Electromagnetic interference detection method for earthquake monitoring

Abstract

The invention belongs to the technical field of geophysical monitoring, and relates to an electromagnetic interference detection method for earthquake monitoring, which comprises the following steps: deploying a multi-dimensional electromagnetic field sensing unit and an earthquake signal acquisition unit and realizing high-precision synchronization; synchronously acquiring multi-dimensional electromagnetic field and seismic physical signal data; performing multi-feature and noise feature extraction on the two; inputting the features into a pre-trained machine learning / deep learning model for identification, and outputting information such as interference type, intensity and azimuth angle; and carrying out influence evaluation and quality marking on the seismic data according to the output, and selectively inhibiting or compensating. Through the above scheme, multi-dimensional perception and characterization, high-precision synchronization and correlation analysis, intelligent identification and spatial positioning are realized, and the reliability and accuracy of earthquake monitoring data are significantly improved.

Description

Technical Field

[0001] This invention belongs to the field of geophysical monitoring technology, and more specifically, relates to an electromagnetic interference detection method for earthquake monitoring. Background Technology

[0002] With the continuous improvement of global earthquake monitoring networks, advancements in earthquake monitoring technology have become a core issue in Earth science research and disaster prevention and mitigation. Against this backdrop, the accuracy and reliability of data acquisition from various sensors and supporting equipment used for collecting and analyzing seismic wave signals are crucial. However, with the acceleration of industrialization and the widespread adoption of wireless communication technologies, the electromagnetic environment is becoming increasingly complex, and the resulting electromagnetic interference (EMI) is having a significant impact on highly sensitive earthquake monitoring equipment. EMI can intrude into earthquake monitoring systems in various forms, leading to signal distortion, increased noise, and even data loss, thereby severely weakening the effectiveness of monitoring data and the accuracy of subsequent analysis, directly affecting the quality of earthquake early warning and scientific research. Therefore, developing an efficient and accurate electromagnetic interference detection method to ensure the purity and reliability of earthquake monitoring data has become a critical problem urgently needing to be solved in the field of earthquake monitoring technology.

[0003] To address the challenges of electromagnetic interference in various application scenarios, multiple solutions have emerged in existing technologies. For example, Chinese Patent Publication No. CN107204822B discloses an unmanned surface vessel (USV) with electromagnetic interference detection capabilities and a method for dealing with electromagnetic interference. The core of this solution lies in the USV's onboard electromagnetic interference detection antenna, which senses the intensity of electromagnetic interference in the surrounding environment in real time. Once interference is detected, the USV intelligently adjusts its position and the direction of its directional receiving antenna according to a preset strategy. This aims to shield or avoid electromagnetic interference sources through spatial position changes and antenna directivity optimization, thereby re-establishing and maintaining a stable communication link. This technical solution effectively suppresses electromagnetic interference in specific fields, especially for maritime communication or remote sensing platforms, through proactive physical avoidance, demonstrating its unique advantages in dynamic and maneuverable scenarios.

[0004] For example, Chinese Patent Publication No. CN112964956B discloses a method and apparatus for detecting electromagnetic interference in vehicles. This solution focuses on the complex electronic systems inside vehicles, using the number of erroneous frames in the vehicle under test as a characterization of electromagnetic interference. Furthermore, it innovatively employs a diagnostic mode that controls multiple onboard components to sequentially enter a non-operating state. Through this step-by-step elimination method, the solution can effectively locate and determine the specific component generating electromagnetic interference. This method requires no additional on-site testing equipment installation, enabling rapid diagnosis of electromagnetic interference inside vehicles, greatly improving the efficiency and convenience of troubleshooting onboard electronic systems, and making a significant contribution to ensuring the stable operation of vehicle electronic systems.

[0005] However, as earthquake monitoring technology evolves towards higher precision, wider coverage, and longer-term observation series, and with the increasing demands for data reliability, some inherent characteristics of the aforementioned existing technologies at the principle level have gradually revealed their inherent limitations and deep-seated contradictions when applied to fixed earthquake monitoring scenarios. Specifically, the unmanned surface vessel (USV) method described in CN107204822B relies heavily on "mobility" for its electromagnetic interference avoidance strategy. Shielding interference sources by adjusting the USV's spatial position and antenna directionality is its fundamental approach to solving the problem. However, earthquake monitoring equipment, especially seismographs in reference stations or regional networks, is typically deployed in a fixed manner. Once their location is selected, it cannot be changed arbitrarily, and the direction of their antennas (such as satellite communication antennas) must also be precisely fixed to ensure data transmission. This fundamental contradiction between "fixedness" and the "mobility" of the USV solution makes this method unsuitable for direct application to fixed seismic stations, and its core interference response mechanism completely fails in non-mobile scenarios. Furthermore, in complex electromagnetic environments, this scheme is inadequate in its ability to identify multi-source interference from different directions, frequencies, and characteristics, making it difficult to meet the needs of earthquake monitoring for refined, multi-dimensional electromagnetic field analysis.

[0006] Correspondingly, the vehicle electromagnetic interference detection method described in CN112964956B fundamentally relies on the "number of error frames" as an indicator to determine and trace the source of interference. Error frames are a direct manifestation of data transmission damage in digital communication systems and are suitable for diagnosing the integrity of communication links. However, the core task of earthquake monitoring equipment is to acquire weak vibration signals from the Earth's interior with high precision. These signals are typically extremely low-frequency, broadband, and weak analog signals or high-sampling-rate digital signals. The impact of electromagnetic interference on earthquake signals is not simply to cause "error frames" or communication interruptions, but rather manifests as the superposition of indistinguishable noise into the target signal, or to cause nonlinear responses in the sensor itself, leading to deviations or distortions in the measurement results. Therefore, this judgment logic based on digital communication "error frames" is fundamentally incompatible with the unique signal characteristics and interference manifestations of earthquake monitoring equipment. Furthermore, this scheme focuses on identifying components inside the vehicle that generate interference, rather than locating interference sources and their spatial distribution in the external electromagnetic environment. In earthquake monitoring applications, interference sources are often external and variable, such as urban industrial activities, power lines, communication base stations, and even atmospheric ionospheric disturbances, and the intensity and direction of the interference are dynamically changing. This method lacks the ability to accurately locate and characterize the direction, intensity, and spectral characteristics of electromagnetic interference sources. This makes it difficult to effectively distinguish electromagnetic noise from real seismic signals in complex electromagnetic environments, and also fails to provide a reliable basis for taking targeted external interference suppression measures. This lack of diagnostic logic and location capability results in the inability to effectively identify, separate, and quantitatively analyze complex electromagnetic noise in the external environment, thus severely restricting the reliability and accuracy of earthquake monitoring data.

[0007] The aforementioned issues collectively demonstrate that existing electromagnetic interference detection technologies, in terms of principles, application scope, and problem-solving focus, exhibit significant fundamental differences and technological gaps compared to the needs of seismic monitoring equipment operating in fixed, high-precision environments with multi-source external interference. This is particularly evident in fixed seismic monitoring scenarios lacking mobility and unable to rely on communication error frames as the primary basis for judgment, highlighting the challenge of accurately identifying, effectively characterizing, and distinguishing complex, multi-source electromagnetic interference from seismic signals. Summary of the Invention

[0008] To achieve the aforementioned objectives, namely, to provide an electromagnetic interference detection method for earthquake monitoring, this invention aims to improve detection accuracy, enhance multi-source interference identification capabilities, and adapt to the application scenarios of earthquake monitoring equipment, thereby improving the reliability and accuracy of earthquake monitoring data. This invention provides an electromagnetic interference detection method capable of high-precision, multi-dimensional, real-time synchronous perception, characterization, location, and classification of electromagnetic interference in the environment where fixed-deployment earthquake monitoring equipment is located, and performing correlation analysis with earthquake monitoring data, thereby achieving effective identification of electromagnetic interference and accurate differentiation of earthquake signals.

[0009] The electromagnetic interference detection method for earthquake monitoring provided by this invention includes the following steps:

[0010] First, at least one multi-dimensional electromagnetic field sensing unit and at least one seismic signal acquisition unit are deployed at the location of the earthquake monitoring equipment or in a predetermined area around it, and the time synchronization accuracy between the multi-dimensional electromagnetic field sensing unit and the seismic signal acquisition unit is better than 1 microsecond.

[0011] Second, through the multi-dimensional electromagnetic field sensing unit, multi-dimensional electromagnetic field data within the region are collected synchronously. The multi-dimensional electromagnetic field data includes, but is not limited to, electric field strength data, magnetic field strength data, electromagnetic wave polarization data, and electromagnetic wave angle of arrival data. The multi-dimensional electromagnetic field sensing unit has wideband response characteristics and can cover the electromagnetic spectrum range from extremely low frequency (ELF) to ultra-high frequency (UHF), for example, from 1Hz to 3GHz.

[0012] Third, the seismic signal acquisition unit synchronously acquires seismic physical signal data within the area, including but not limited to ground motion velocity, ground motion acceleration, or ground strain; the seismic signal acquisition unit has high sensitivity and wide dynamic range characteristics, and can respond to seismic signal frequency bands from 0.01Hz to 100Hz.

[0013] Fourth, multi-feature extraction and analysis are performed on the multi-dimensional electromagnetic field data to generate multi-dimensional electromagnetic interference feature vectors; the multi-feature extraction includes, but is not limited to, time-domain feature extraction, frequency-domain feature extraction, time-frequency-domain feature extraction, polarization-domain feature extraction, and spatial-domain feature extraction.

[0014] Fifth, noise features are extracted from the earthquake physical signal data to generate earthquake signal noise feature vectors; the noise feature extraction includes, but is not limited to, analyzing the energy distribution, instantaneous amplitude changes, correlation fluctuations, and abnormal pulses caused by non-earthquake sources in the earthquake physical signal data within a specific frequency band;

[0015] Sixth, the multi-dimensional electromagnetic interference feature vector and the seismic signal noise feature vector are input into a pre-trained electromagnetic interference identification model; the electromagnetic interference identification model is constructed based on machine learning or deep learning algorithms and trained using historical, labeled electromagnetic interference data samples and clean seismic signal data samples to learn the mapping relationship and discrimination boundary between electromagnetic interference features and seismic signal noise features.

[0016] Seventh, through the electromagnetic interference identification model, output the probability of the presence of electromagnetic interference, type classification, intensity estimation, dominant frequency, temporal correlation and possible spatial azimuth information;

[0017] Eighth, based on the output of the electromagnetic interference identification model, the seismic physical signal data is assessed for electromagnetic interference impact and quality is marked, and electromagnetic interference suppression or compensation processing can be selectively performed.

[0018] Further, the multi-dimensional electromagnetic field sensing unit in the first step of this invention specifically includes: a triaxial orthogonal electric field sensor array for acquiring the electric field intensity data; a triaxial orthogonal magnetic field sensor array for acquiring the magnetic field intensity data; a broadband radio frequency receiving front end for receiving and digitizing the electromagnetic field data; and a high-precision time synchronization module for achieving time alignment between the electromagnetic field data and the seismic physical signal data. The triaxial orthogonal electric field sensor array consists of three mutually perpendicular electric dipole antennas, each with a symmetrical arm of 100 mm to 500 mm in length, a gain flatness maintained within ±3 dB in the range of 1 Hz to 3 GHz, an input impedance of 50 ohms, and a noise figure of less than 3 dB. The triaxial orthogonal magnetic field sensor array consists of three mutually perpendicular single-turn or multi-turn loop coil antennas, each with a diameter of 50 mm to 200 mm. Its induced voltage is proportional to the rate of change of magnetic flux, with a response frequency range of 1 Hz to 100 MHz, a sensitivity of no less than 10 mV / nT at 10 kHz, and an inherent noise of less than 1 pT / √Hz. The broadband RF receiving front-end includes a low-noise amplifier (LNA), a mixer, an intermediate frequency filter, and a high-speed analog-to-digital converter (ADC). The LNA has a noise figure of less than 2 dB and a gain of 20 dB to 40 dB. The high-speed ADC employs 14-bit or 16-bit resolution and a sampling rate of 200 MSPS to 5 GSPS to ensure sufficient sampling of broadband electromagnetic signals. The high-precision time synchronization module uses a local high-precision crystal oscillator or atomic clock (e.g., a rubidium atomic clock) calibrated by a Global Navigation Satellite System (GNSS) receiver to achieve nanosecond or sub-nanosecond time synchronization with other units through Network Time Protocol (NTP) or Precision Time Protocol (PTP), ensuring the consistency of timestamps across all data acquisition channels.

[0019] Further, the seismic signal acquisition unit in the third step of this invention specifically includes: a triaxial broadband velocity-type seismic sensor for acquiring ground motion velocity data, with a frequency response range of 0.01Hz to 100Hz, a sensitivity of not less than 2000V / (m / s), and a self-noise level lower than the Low Noise Model for Seismology (NLNM); a high-resolution analog-to-digital converter for digitizing the seismic physical signal data; and a data preprocessing module for filtering and gain adjustment of the digitized seismic physical signal data. The high-resolution analog-to-digital converter uses a 24-bit or 26-bit Σ-Δ ADC with an effective bit depth of not less than 22 bits, a sampling rate configurable to 100SPS, 200SPS, or 500SPS, and a total harmonic distortion plus noise (THD+N) of less than -110dB. The data preprocessing module integrates a hardware digital filter, which can realize low-pass, high-pass, band-pass, or band-stop filtering functions. For example, to remove power line interference, it can be configured as a band-stop filter with a center frequency of 50Hz or 60Hz and a stopband attenuation depth greater than 80dB.

[0020] Furthermore, the multi-feature extraction and analysis in the fourth step of this invention specifically includes:

[0021] First, time-domain feature extraction: Instantaneous power, peak amplitude, root mean square (RMS) amplitude, zero-crossing rate, pulse width, rise time, and fall time are calculated from the multi-dimensional electromagnetic field data. The instantaneous power is obtained by calculating the square of the instantaneous amplitude of the electric or magnetic field. The peak amplitude represents the maximum absolute value of the signal within a specified time window. The RMS amplitude represents the effective value of the signal power within a specified time window. The zero-crossing rate represents the number of times the signal crosses the zero level within a specified time window, used to characterize the frequency variation characteristics of the signal. The pulse width, rise time, and fall time are used to characterize the duration and waveform steepness of instantaneous pulse-type interference.

[0022] Second, frequency domain feature extraction: Fast Fourier Transform (FFT) or Welch spectrum estimation is performed on the multi-dimensional electromagnetic field data to obtain features such as power spectral density (PSD), dominant frequency, harmonic components, bandwidth, spectral centroid, and peak-to-sidelobe ratio; the calculation window length of the power spectral density curve can be set to 2 to 10 seconds, with an overlap rate of 50% to 75%; the dominant frequency is the frequency point corresponding to the peak value of the power spectral density; the harmonic component identification is achieved by detecting energy peaks that appear at integer multiples of the fundamental frequency;

[0023] The bandwidth is the frequency range at which the signal power drops to half (-3dB) or one-tenth (-10dB) of the peak power.

[0024] Third, time-frequency domain feature extraction: performing a short-time Fourier transform on the multi-dimensional electromagnetic field data.

[0025] (STFT) or wavelet transform is used to obtain time-varying spectral characteristics, instantaneous frequency and energy concentration, etc., to identify transient, time-varying or non-stationary electromagnetic interference; the STFT adopts Hanning window or Gaussian window, the window length can be configured from 20 milliseconds to 500 milliseconds, and the window overlap rate is 50%.

[0026] The wavelet transform employs the Morlet wavelet or Daubechies wavelet basis to analyze the local features of the signal at different scales.

[0027] Fourth, polarization domain feature extraction: Based on the output of the triaxial orthogonal electric field sensor array or magnetic field sensor array, the Stokes parameters of the electromagnetic wave are calculated to obtain features such as linear polarization degree, circular polarization degree, and polarization angle, which are used to identify electromagnetic interference sources of specific modulation types; the Stokes parameters are obtained by calculating the cross-correlation function and the autocorrelation function, for example, I = E_x^2 + E_y^2, Q = E_x^2 - E_y^2, U = 2Re(E_x E_y^), V = 2Im(E_x E_y^).

[0028] Fifth, spatial domain feature extraction: If multiple multi-dimensional electromagnetic field sensing units are deployed, multi-station time difference of arrival (TDOA), angle of arrival (AoA), or electromagnetic beamforming techniques are used to locate and track electromagnetic interference sources, thereby obtaining features such as the spatial azimuth, elevation angle, distance, and movement trajectory of the interference sources; the AoA technique can employ MUSIC (Multiple Signal Classification) or ESPRIT (Estimation of Signal Parameters via Rotational Invariance).

[0029] The Techniques algorithm estimates the direction of arrival of a signal by analyzing the phase difference between signals from multiple receiving antennas; the TDOA technology calculates the location of the interference source by measuring the time difference between the arrival of the same interference signal from different sensors, requiring at least three synchronous sensors for two-dimensional positioning or at least four synchronous sensors for three-dimensional positioning.

[0030] Furthermore, the noise feature extraction in the fifth step of this invention specifically includes:

[0031] First, the seismic physical signal data undergoes preprocessing operations such as detrending, DC removal, and instrument response removal. The detrending is achieved by fitting the signal using the least squares method and removing linear or polynomial trend terms. The DC removal is achieved using a high-pass filter with a cutoff frequency of, for example, 0.001 Hz. The instrument response removal is achieved by deconvolution or frequency domain division to remove the influence of the seismograph's transfer function.

[0032] Second, within a specific time window, power spectral density analysis is performed on the preprocessed seismic physical signal data to identify anomalous energy concentration frequency bands that exceed the background seismic noise spectrum; the power spectral density analysis employs FFT or Welch spectral estimation methods similar to the electromagnetic interference feature extraction, and the window length matches the electromagnetic field data analysis window.

[0033] Third, the energy ratio or signal-to-noise ratio decrease of the seismic physical signal data in preset electromagnetic interference sensitive frequency bands (e.g., 50Hz, 60Hz, FM broadcast band, cellular communication band) is calculated to quantify the potential electromagnetic interference impact; the sensitive frequency bands are identified by analyzing historical electromagnetic environment data to identify common electromagnetic interference frequencies that affect seismic data.

[0034] Fourth, analyze the instantaneous pulses or spikes in the earthquake physical signal data, characterized by steep rising edges, narrow pulse widths, and high instantaneous amplitudes, and distinguish them from typical earthquake signal waveforms;

[0035] The instantaneous pulse can be detected using methods based on amplitude thresholds or specific waveform template matching. Fifth, through cross-correlation analysis or coherence analysis, the correlation between noise components in the seismic physical signal data and specific frequency bands or types of electromagnetic interference characteristics in the multi-dimensional electromagnetic field data is evaluated.

[0036] The cross-correlation analysis calculates the time delay and correlation coefficient between seismic noise and electromagnetic interference signals; the coherence analysis calculates the linear correlation between the two in the frequency domain to determine whether electromagnetic interference is the main source of seismic signal noise.

[0037] Furthermore, the electromagnetic interference identification model in the sixth step of this invention specifically includes:

[0038] First, the input layer receives the multi-dimensional electromagnetic interference feature vector (e.g., feature values ​​in the time domain, frequency domain, time-frequency domain, polarization domain, and spatial domain, with a total dimension of up to hundreds of dimensions) and the seismic signal noise feature vector (e.g., power spectral density, energy ratio, pulse characteristics, and cross-correlation coefficients, with a total dimension of up to tens of dimensions), and merges these feature vectors into a high-dimensional input vector.

[0039] Second, the feature fusion layer: a convolutional neural network (CNN) or a multilayer perceptron (MLP) is used to perform preliminary feature learning and dimensionality reduction on the merged high-dimensional input vector to extract higher-level abstract features; the CNN layer may contain multiple convolutional kernels (e.g., 16 to 64, with a size of 3x3 or 5x5), the activation function is ReLU, and pooling layers may be connected for dimensionality reduction; the MLP layer may contain 2 to 5 fully connected layers, with 64 to 512 neurons in each layer.

[0040] Third, the sequence processing layer: Considering the temporal characteristics of electromagnetic interference or its time-varying impact on seismic data, recurrent neural network (RNN) structures such as Long Short-Term Memory (LSTM) or Gated Recurrent Unit (GRU) can be introduced to capture the dependencies of feature vectors over time. The LSTM layer can contain 1 to 3 layers, each containing 64 to 256 units, for processing continuous feature sequences. Fourth, the classification decision layer: The Softmax function is used as the activation function of the output layer to output the probability distribution of different electromagnetic interference types (e.g., power line harmonic interference, broadcast signal interference, communication base station interference, industrial equipment transient interference, lightning discharge interference, background electromagnetic noise), and to identify the state of clean seismic signals without electromagnetic interference. The classification decision layer is a fully connected layer, and the number of neurons is equal to the number of predefined electromagnetic interference types plus one (used to represent the interference-free state).

[0041] Fifth, training optimization: The model is trained using gradient descent optimization algorithms (such as Adam optimizer or RMSprop optimizer) and cross-entropy loss function. During training, techniques such as batch normalization and dropout are used to prevent overfitting. The training dataset contains electromagnetic interference samples and clean seismic signal samples collected from different earthquake monitoring stations under different electromagnetic environments, which are manually or semi-automatically labeled, with a sample size of at least tens of thousands.

[0042] Furthermore, the electromagnetic interference impact assessment and quality marking in the eighth step of this invention specifically includes:

[0043] First, based on the output probability of the electromagnetic interference identification model, when the probability of a specific type of electromagnetic interference exceeds a preset threshold (e.g., 0.8), it is determined that the seismic data for that time period is affected by that type of electromagnetic interference.

[0044] Second, by combining the intensity estimate, dominant frequency, and time correlation of the electromagnetic interference, the percentage decrease in data purity or the decibel reduction in signal-to-noise ratio (SNR) caused by the electromagnetic interference to the seismic physical signal data is calculated; the decibel reduction in SNR is calculated by comparing the power spectrum or root mean square value of the seismic signal during the interference period and the interference-free period.

[0045] Third, for seismic physical signal data segments affected by electromagnetic interference, a standardized quality flag or label is added. For example, the quality control word defined in the SEED (Standard for the Exchange of Earthquake Data) data format is used to identify the data segment as "affected by electromagnetic interference", "high noise", "pending review", etc. The quality flag may include the interference type code, intensity level, start and end points of the affected period, and suggested follow-up processing measures (e.g., suggesting specific frequency filtering or skipping the data segment).

[0046] Fourth, the electromagnetic interference suppression or compensation processing can be selectively implemented using digital signal processing methods such as adaptive noise cancellation algorithms, wavelet domain denoising, or notch filtering. The adaptive noise cancellation algorithm uses LMS (Least Mean Square) or RLS (Recursive Least Squares) algorithms, taking the multi-dimensional electromagnetic field data as a reference input, and subtracting the relevant electromagnetic interference components from the seismic physical signal data. The center frequency and bandwidth of the notch filter are dynamically adjusted according to the identified dominant frequency and interference bandwidth.

[0047] The present invention achieves significant technical effects through the above-described technical solution, specifically in the following aspects:

[0048] First, this invention achieves "multi-dimensional" perception and characterization of electromagnetic interference, overcoming the limitations of existing technologies that only focus on a single indicator (such as communication error frames). This invention not only collects the electric and magnetic field components of the electromagnetic field but also analyzes its time-domain, frequency-domain, time-frequency-domain, polarization-domain, and spatial-domain characteristics, constructing a comprehensive portrait of the electromagnetic environment. This multi-dimensional feature set can more precisely distinguish different types of electromagnetic interference sources, such as power frequency harmonic interference, radio frequency broadcast signals, impulse noise, or partial discharge noise, significantly improving the accuracy and classification capability of electromagnetic interference identification.

[0049] Secondly, it achieves "high-precision synchronous acquisition and correlation analysis" of electromagnetic interference data and seismic signal data, solving the problem that electromagnetic interference detection and seismic signal acquisition are independent and difficult to correlate effectively in existing technologies. Through nanosecond or sub-nanosecond precise time synchronization technology, it ensures the time alignment of electromagnetic interference events with corresponding noise segments in the seismic signal. This close temporal correlation makes subsequent cross-analysis possible, enabling a clear determination of the scope and extent of the impact of specific electromagnetic interference events on seismic monitoring data, thus providing a solid foundation for data cleaning and quality assessment.

[0050] Third, an "electromagnetic interference identification model" based on machine learning or deep learning was constructed, solving the problem that traditional threshold discrimination or simple filtering methods are unable to effectively distinguish between electromagnetic noise and weak seismic signals in complex electromagnetic environments. This model can learn and identify complex patterns of electromagnetic interference in a multi-dimensional feature space, achieving intelligent classification and accurate identification of electromagnetic interference even in situations with low signal-to-noise ratios or the simultaneous presence of multiple interference sources. Through training on a large amount of labeled data, the model possesses the ability to generalize to novel or unseen interference patterns, significantly improving the robustness of detection.

[0051] Fourth, it provides the capability of "spatial location and azimuth estimation" of electromagnetic interference sources, solving the dilemma of being unable to avoid interference sources in fixed seismic monitoring scenarios. By deploying multiple electromagnetic field sensing units and combining them with multi-station positioning algorithms, this invention can estimate the direction and distance of interference sources in real time. This capability not only helps to gain a deeper understanding of the physical source of interference (e.g., whether it is a nearby power line, communication tower, or distant industrial equipment), but also provides data support for the site selection optimization, shielding design, or directional interference suppression of seismic stations. Although the location of the seismograph is fixed, understanding the direction of the interference source can guide the modification of the surrounding environment of the station, or allow for focused attention on interference from specific directions during data processing.

[0052] Fifth, this invention enables "electromagnetic interference impact assessment and quality labeling" of earthquake monitoring data, significantly improving the reliability and accuracy of the data. The invention can output quantitative indicators of the type, intensity, dominant frequency, and specific impact of electromagnetic interference events on the signal-to-noise ratio of earthquake data, and attach standard quality flags to affected data segments. This allows subsequent data users to clearly identify data that may be contaminated by electromagnetic interference, avoiding misleading earthquake event interpretation and focal mechanism research. It also provides targeted guidance for data preprocessing, such as selectively denoising specific frequency bands or directly removing contaminated data, thereby ensuring the purity of earthquake monitoring data and the effectiveness of scientific analysis.

[0053] Sixth, the method proposed in this invention possesses "universality and adaptability," and can be applied to various fixed-deployment earthquake monitoring scenarios, including regional seismic networks, deep-well seismic stations, and ocean bottom seismometers. Its core lies in the refined perception and intelligent analysis of the electromagnetic environment, unrestricted by the physical mobility of monitoring equipment, and independent of non-seismic signal characteristics such as communication error rates. Therefore, it effectively overcomes the inherent contradictions and limitations of existing technologies in fixed, high-precision earthquake monitoring scenarios. Attached Figure Description

[0054] To more clearly illustrate the technical solutions in the embodiments of the present invention, the accompanying drawings used in the description of the embodiments will be briefly introduced below. Obviously, the accompanying drawings described below are only some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0055] Figure 1 This is a flowchart illustrating the electromagnetic interference detection method for earthquake monitoring according to the present invention. Detailed Implementation

[0056] This invention provides an electromagnetic interference detection method for earthquake monitoring. It aims to achieve effective identification of electromagnetic interference and accurate differentiation of weak seismic signals by performing multi-dimensional, high-precision, real-time synchronous sensing, characterization, localization, and classification of the electromagnetic field in the environment where the earthquake monitoring equipment is located, and then conducting deep correlation analysis with synchronously acquired seismic physical signal data. This method not only improves the accuracy and multi-source identification capability of electromagnetic interference detection, but more importantly, it provides solid technical support for the quality assessment, contamination identification, and subsequent processing of earthquake monitoring data, significantly enhancing the reliability and accuracy of earthquake monitoring data. The system architecture and execution flow of this method can be found in the appendix. Figure 1 .

[0057] In one specific embodiment, the electromagnetic interference detection method for earthquake monitoring provided by the present invention includes the following key steps in its operation process.

[0058] First, at least one multi-dimensional electromagnetic field sensing unit 10 and at least one seismic signal acquisition unit 20 are precisely deployed within the geographical location of the fixedly deployed seismic monitoring equipment or within a predetermined radius therearound. The selection of deployment location is crucial, requiring comprehensive consideration of the potential distribution of electromagnetic interference sources, the propagation characteristics of seismic signals, and the power supply and communication conditions of the equipment. Preferably, the multi-dimensional electromagnetic field sensing unit 10 and the seismic signal acquisition unit 20 should be deployed as close as possible to ensure a high spatial correlation between the electromagnetic environment they sense and the seismic physical signals, thereby avoiding a decrease in data correlation due to excessive spatial separation. To ensure the accuracy and effectiveness of subsequent data analysis, the time synchronization accuracy between the multi-dimensional electromagnetic field sensing unit 10 and the seismic signal acquisition unit 20 is strictly controlled, requiring a precision better than 1 microsecond. Specifically, this high-precision time synchronization is achieved through a high-precision time synchronization module 14 built into the multi-dimensional electromagnetic field sensing unit 10. This module typically uses a timing signal provided by a Global Navigation Satellite System (GNSS) receiver (e.g., a receiver supporting multi-constellation positioning functions such as GPS, BeiDou, and GLONASS) to calibrate a local high-precision crystal oscillator or atomic clock. The atomic clock can be a rubidium atomic clock with high robustness and extremely low drift rate, with a long-term stability better than 5 x 10^-11 / day and a short-term stability better than 1 x 10^-12 over a 1-second average time. In actual deployment, the GNSS receiver receives satellite signals through an external antenna and sends timing pulses (e.g., 1 PPS, i.e., pulses per second) to the local clock module for synchronization calibration. The calibrated local clock broadcasts time information to all connected sensing and acquisition units via a high-precision time protocol (e.g., IEEE 1588 Precision Time Protocol, PTPv2) or Network Time Protocol (NTP), ensuring that all acquired data packets are appended with nanosecond or sub-nanosecond timestamps. For example, using the PTP protocol, the clock offset between units can be kept within 50 nanoseconds, which is far better than the 1 microsecond accuracy requirement, laying a solid time foundation for subsequent accurate correlation analysis of electromagnetic interference and seismic signals.

[0059] Second, the multi-dimensional electromagnetic field sensing unit 10 synchronously and continuously collects multi-dimensional electromagnetic field data within the region. This multi-dimensional electromagnetic field data constitutes a comprehensive electromagnetic environment profile, including but not limited to electric field strength data, magnetic field strength data, electromagnetic wave polarization data, and electromagnetic wave angle of arrival data. This multi-dimensional sensing capability is key to the invention's ability to precisely distinguish complex electromagnetic interference sources. The multi-dimensional electromagnetic field sensing unit 10 features a wide-bandwidth response characteristic in its hardware design, capable of covering a broad electromagnetic spectrum range from extremely low frequency (ELF) to ultra-high frequency (UHF), such as from 1Hz power frequency noise and ground current disturbances to high-frequency interference such as 3GHz wireless communication and microwave radiation. Specifically, the multi-dimensional electromagnetic field sensing unit 10 is subdivided into: a triaxial orthogonal electric field sensor array 11 for accurately collecting electric field strength data; a triaxial orthogonal magnetic field sensor array 12 for capturing magnetic field strength data; a broadband radio frequency receiving front-end 13 responsible for receiving and digitizing electromagnetic field data; and the aforementioned high-precision time synchronization module 14. The triaxial orthogonal electric field sensor array 11 consists of three physically perpendicular electric dipole antennas. The symmetrical arm length of each electric dipole antenna is preferably designed to be 300 mm to balance low-frequency response and high-frequency characteristics. These antennas maintain a gain flatness within ±3 dB over a wide frequency range of 1 Hz to 3 GHz, have an input impedance of 50 ohms, and a noise figure strictly controlled to less than 3 dB to ensure low signal distortion and high fidelity. The triaxial orthogonal magnetic field sensor array 12 consists of three physically perpendicular single-turn or multi-turn loop coil antennas. In a preferred embodiment, each coil antenna is wound with 50 turns of enameled wire, has a diameter of 100 mm, and its induced voltage is strictly proportional to the rate of change of magnetic flux. The response frequency range covers from 1 Hz to 100 MHz, achieving a sensitivity of no less than 10 mV / nT at 10 kHz, with an inherent noise level of less than 1 pT / √Hz. This is crucial for capturing weak magnetic field changes (such as ground currents or leakage magnetic fields from industrial equipment). The broadband RF receiver front-end 13 is the core of the signal link, comprising a low-noise amplifier (LNA), a mixer, an intermediate frequency (IF) filter, and a high-speed analog-to-digital converter (ADC). The LNA has a noise figure better than 2dB, and its gain can be flexibly configured between 20dB and 40dB to accommodate electromagnetic signals of varying intensities, preventing weak signals from being drowned out by noise. The mixer down-converts the RF signal to the IF frequency, and then the IF filter removes out-of-band noise. The high-speed ADC employs 14-bit resolution and a sampling rate up to 2GSPS, enabling sufficient Nyquist sampling of electromagnetic signals within a 3GHz bandwidth to ensure complete signal digitization.

[0060] Third, the seismic signal acquisition unit 20 synchronously and continuously acquires seismic physical signal data within the region. This seismic physical signal data encompasses various physical manifestations of Earth's internal activity and surface motion, including but not limited to ground motion velocity, ground motion acceleration, and ground strain. The seismic signal acquisition unit 20 is designed to possess high sensitivity and a wide dynamic range, effectively responding to near-field high-frequency seismic signals ranging from ultra-long-period seismic waves of 0.01Hz to near-field high-frequency seismic signals of 100Hz. Specifically, the seismic signal acquisition unit 20 is subdivided into: a three-axis broadband velocity-type seismic sensor 21 for acquiring ground motion velocity data; a high-resolution analog-to-digital converter 22 for digitizing the seismic physical signal data; and a data preprocessing module 23 for performing preliminary filtering and gain adjustment on the digitized seismic physical signal data. The three-axis broadband velocity seismic sensor 21, such as the Guralp CMG-3ESPC broadband seismograph, has a frequency response range accurately covering 0.01Hz to 100Hz, a sensitivity of no less than 2000V / (m / s), and a self-noise level far below the New Low Noise Model (NLNM), ensuring its ability to capture extremely weak ground vibrations. The high-resolution analog-to-digital converter 22 employs a 26-bit Σ-Δ ADC chip (e.g., Analog Devices AD7768-1), with an effective bit depth of 22 or more, configurable sampling rates including 100SPS, 200SPS, or 500SPS, and a total harmonic distortion plus noise (THD+N) of less than -110dB, ensuring ultra-high fidelity digitization of the seismic signal. The data preprocessing module 23 integrates configurable digital filters at the hardware level, enabling low-pass, high-pass, band-pass, or band-stop filtering functions. For example, to effectively suppress common power line AC interference (such as 50Hz or 60Hz), it can be configured as a narrowband bandstop filter with a center frequency of 50Hz, and its stopband attenuation depth can reach more than 80dB, thereby initially purifying seismic data at the hardware level and reducing the burden of subsequent software processing.

[0061] Fourth, multi-feature extraction and analysis are performed on the multi-dimensional electromagnetic field data to generate multi-dimensional electromagnetic interference feature vectors. This step is the cornerstone for identifying different types of electromagnetic interference by uncovering the inherent laws of electromagnetic signals from different dimensions. Multi-feature extraction methods cover time-domain feature extraction, frequency-domain feature extraction, time-frequency-domain feature extraction, polarization-domain feature extraction, and spatial-domain feature extraction.

[0062] Specifically, time-domain feature extraction includes calculating features such as instantaneous power, peak amplitude, root mean square (RMS) amplitude, zero-crossing rate, pulse width, rise time, and fall time from electromagnetic field data. Instantaneous power is obtained by calculating the square of the instantaneous amplitude of the electric or magnetic field, for example, P(t) = E(t)^2 or P(t) = H(t)^2. Peak amplitude represents the maximum absolute value of the signal within a specified time window (e.g., 10 milliseconds), reflecting the intensity of the instantaneous interference. RMS amplitude represents the effective value of the signal power within a specified time window, for example, RMS = sqrt(mean(signal^2)), which provides information about the average energy of the signal. Zero-crossing rate represents the number of times the signal crosses the zero level within a specified time window, used to characterize the frequency variation characteristics of the signal; a high zero-crossing rate usually indicates high-frequency components. Pulse width, rise time, and fall time are used to accurately characterize the duration of instantaneous pulse-type interference and the steepness of its waveform; these features are crucial for distinguishing transient discharges, switching noise, etc.

[0063] Frequency domain feature extraction involves performing Fast Fourier Transform (FFT) or Welch spectrum estimation on multi-dimensional electromagnetic field data to obtain features such as power spectral density (PSD), dominant frequency, harmonic components, bandwidth, spectral centroid, and peak-to-sidelobe ratio. The calculation window length for the power spectral density curve can be set to 5 seconds with an overlap rate of 75% to balance frequency resolution and temporal smoothness. The dominant frequency is the frequency point corresponding to the energy peak in the power spectral density curve. Harmonic component identification is achieved by detecting energy peaks that occur at integer multiples of the fundamental frequency, such as identifying power frequency harmonics like 50Hz, 100Hz, and 150Hz. The bandwidth is defined as the frequency range where the signal power drops to half of the peak power (-3dB), or, depending on the specific application requirements, the frequency range where it drops to one-tenth (-10dB). The spectral centroid indicates the average position of the signal energy in the frequency domain.

[0064] Time-frequency domain feature extraction involves performing Short-Time Fourier Transform (STFT) or Wavelet Transform on multi-dimensional electromagnetic field data to obtain time-varying spectral characteristics, instantaneous frequency, and energy concentration. This method is particularly suitable for identifying transient, time-varying, or non-stationary electromagnetic interference. STFT employs a Hanning window with a configurable window length of 100 milliseconds and a window overlap rate of 50%, achieving a balance between time and frequency resolution. Wavelet Transform uses Morlet wavelet basis functions; its excellent time-frequency localization properties enable effective analysis of signal local features at different scales, making it particularly suitable for capturing the fine structure of transient interference (such as lightning discharges and switching noise). Through time-frequency analysis, the evolution of interference events in time and frequency can be clearly displayed, for example, for identifying frequency sweeping or frequency hopping interference.

[0065] Polarization domain feature extraction is based on the output of the triaxial orthogonal electric field sensor array 11 or magnetic field sensor array 12. Stokes parameters of the electromagnetic wave are calculated to obtain features such as linear polarization degree, circular polarization degree, and polarization angle, which are used to identify electromagnetic interference sources of specific modulation types. The Stokes parameters include four components: I, Q, U, and V. These are obtained by calculating the cross-correlation and autocorrelation functions of the three orthogonal electric field components (Ex, Ey, Ez). For example, for a plane wave, I = Ex^2 + Ey^2, Q = Ex^2 - Ey^2, U = 2Re(Ex Ey*), V = 2Im(Ex Ey*), where E_x and E_y represent the complex amplitudes of the two orthogonal electric field components. These parameters can be used to quantify the polarization state of electromagnetic waves, such as linear polarization (sqrt(Q^2+U^2) / I) and circular polarization (|V| / I). This information is of great significance for distinguishing interference from different radiation sources (such as vertically polarized antennas, horizontally polarized antennas, or circularly polarized antennas).

[0066] Spatial domain feature extraction, if multiple multi-dimensional electromagnetic field sensing units are deployed, utilizes Time Difference of Arrival (TDOA), Angle of Arrival (AoA), or beamforming techniques to locate and track electromagnetic interference sources, obtaining features such as the spatial azimuth, elevation, distance, and trajectory of the interference source. AoA technology employs the MUSIC (Multiple Signal Classification) algorithm, estimating the direction of arrival of a signal by analyzing the phase difference between signals from multiple receiving antennas, achieving an accuracy of 0.5 degrees. This method requires an array of at least two or more (e.g., four) electric or magnetic field sensors. TDOA technology calculates the location of the interference source by measuring the time difference between arrivals of the same interference signal from different sensors. This method requires extremely high time synchronization accuracy, necessitating at least three synchronized sensors for two-dimensional positioning, or at least four synchronized sensors for three-dimensional positioning, using the hyperbolic positioning principle to solve for the interference source coordinates. These spatial features provide important clues for analyzing the physical origin of the interference source and its relative position to seismic stations.

[0067] Fifth, noise features are extracted from the seismic physical signal data to generate a seismic signal noise feature vector. This step aims to identify noise components introduced by non-seismic sources (especially electromagnetic interference) from the seismic data. The specific process of noise feature extraction includes the following aspects.

[0068] First, the seismic physical signal data undergoes preprocessing operations such as detrending, DC removal, and instrument response removal. Detrending involves fitting and removing linear or polynomial trend terms from the signal using the least squares method to eliminate long-term drift. DC removal is achieved using a high-pass filter, with a cutoff frequency typically set to 0.001 Hz, to remove DC bias from the signal. Instrument response removal uses deconvolution or frequency domain division to remove the influence of the seismograph's transfer function, converting the sensor output into true ground motion and ensuring that subsequent analyses are based on the actual changes in physical quantities.

[0069] Secondly, within a specific time window, such as a 5-second window matching the electromagnetic field data analysis window, power spectral density analysis is performed on the preprocessed seismic physical signal data to identify anomalous energy concentration frequency bands exceeding the background seismic noise spectrum. The power spectral density analysis employs FFT or Welch spectral estimation methods similar to those used for electromagnetic interference feature extraction, with a window overlap rate of 75%. By comparing the real-time calculated power spectrum with a pre-defined background seismic noise model (e.g., NLNM or NHNM, New High Noise Model), anomalously high-energy frequency bands can be identified; these anomalies are typically manifestations of electromagnetic interference or other non-seismic noise.

[0070] Furthermore, the energy proportion or signal-to-noise ratio (SNR) decrease of seismic physical signal data within pre-defined electromagnetic interference sensitive frequency bands is calculated to quantify the potential impact of electromagnetic interference. Sensitive frequency bands are identified through long-term analysis of historical electromagnetic environment data, revealing common electromagnetic interference frequencies that affect seismic data, such as power frequencies (50Hz / 60Hz), FM broadcasting bands (88-108MHz), and cellular communication bands (800 / 900 / 1800 / 2100MHz). By calculating the proportion of signal energy to total energy within these frequency bands, or by comparing the SNR changes in these bands before and after interference, the degree of electromagnetic interference contamination of seismic data can be accurately quantified.

[0071] Furthermore, the analysis of transient pulses or spike-like noise in seismic physical signal data reveals their characteristics: steep rise edges, narrow pulse widths, and high instantaneous amplitudes, distinguishing them from typical seismic signal waveforms. Transient pulse detection can employ amplitude threshold-based methods (e.g., instantaneous peak values ​​exceeding six times the root mean square amplitude) or methods based on specific waveform template matching (e.g., cross-correlation matching with predefined electromagnetic pulse waveforms). Typical seismic signals (such as P-waves and S-waves) typically have relatively gentle rise edges and wider durations, significantly different from the characteristics of electromagnetic pulse noise.

[0072] Finally, cross-correlation analysis or coherence analysis is used to assess the correlation between noise components in seismic physical signal data and specific frequency bands or types of electromagnetic interference characteristics in multidimensional electromagnetic field data. Cross-correlation analysis calculates the time delay and correlation coefficient between seismic noise and electromagnetic interference signals. A high correlation coefficient (e.g., absolute value greater than 0.7) and a time delay close to zero strongly indicate that both originate from the same event. Coherence analysis calculates the linear correlation (between 0 and 1) between the two in the frequency domain to determine whether electromagnetic interference is the primary source of noise in seismic signals. If electromagnetic interference signals in a specific frequency band exhibit high coherence with the noise components of seismic physical signals in that frequency band (e.g., coherence greater than 0.8), it can be inferred that the electromagnetic interference is the primary contributor to noise in the seismic data of that frequency band.

[0073] Sixth, the multi-dimensional electromagnetic interference feature vector and the seismic signal noise feature vector are input into the pre-trained electromagnetic interference identification model 30. The electromagnetic interference identification model 30 is constructed based on machine learning or deep learning algorithms and is rigorously trained using historical, labeled electromagnetic interference data samples and clean seismic signal data samples to learn the complex mapping relationship and discrimination boundary between electromagnetic interference features and seismic signal noise features. The core of this model lies in its multi-layered architecture.

[0074] Specifically, the electromagnetic interference identification model 30 includes the following components.

[0075] First, the input layer 31 receives multi-dimensional electromagnetic interference feature vectors (e.g., feature values ​​in the time domain, frequency domain, time-frequency domain, polarization domain, and spatial domain, with a total dimension of up to 300) and seismic signal noise feature vectors (e.g., power spectral density, energy ratio, impulse characteristics, and cross-correlation coefficients, with a total dimension of up to 50). These feature vectors are merged into a high-dimensional input vector at the input layer, forming a unified data representation for model processing.

[0076] Second, the feature fusion layer 32: A Convolutional Neural Network (CNN) or a Multilayer Perceptron (MLP) is used to perform preliminary feature learning and dimensionality reduction on the merged high-dimensional input vector to extract higher-level abstract features. The CNN layer can contain three convolutional kernel layers. For example, the first convolutional layer contains 32 kernels (5x1 size), the second contains 64 kernels (3x1 size), and the third contains 128 kernels (3x1 size). All convolutional layers use the ReLU activation function, and each convolutional layer is followed by a max-pooling layer with a 2x1 kernel size for dimensionality reduction. This structure helps to automatically extract discriminative local patterns from the original features. The MLP layer can contain three fully connected layers, with 256, 128, and 64 neurons per layer, respectively. It also uses the ReLU activation function to achieve non-linear feature transformation.

[0077] Third, sequence processing layer 33: Considering the temporal characteristics of electromagnetic interference or the time-varying nature of its impact on seismic data, recurrent neural network (RNN) structures such as Long Short-Term Memory (LSTM) networks or gated recurrent units (GRUs) can be introduced to capture the dependencies of feature vectors over time. An LSTM layer can contain two layers, each with 128 units, for processing continuous feature sequences. For example, when analyzing a continuous 5-second feature sequence, LSTM can understand the correlation between the interference features at the current time point and the interference features at the previous second or even longer time points, thereby improving the ability to identify intermittent or periodic interference.

[0078] Fourth, the classification decision layer 34: The Softmax function is used as the activation function for the output layer to output the probability distribution of different electromagnetic interference types (e.g., power line harmonic interference, broadcast signal interference, communication base station interference, industrial equipment transient interference, lightning discharge interference, and background electromagnetic noise), and to identify the state of clean seismic signals without electromagnetic interference. The classification decision layer is a fully connected layer with the number of neurons equal to the predefined number of electromagnetic interference types (e.g., 6 types) plus one (used to represent the interference-free state), i.e., 7 neurons.

[0079] Fifth, training optimization: The model is trained using gradient descent optimization algorithms (e.g., the Adam optimizer) and cross-entropy loss function. The Adam optimizer, with its adaptive learning rate adjustment, can accelerate model convergence and improve training stability. Batch normalization is used during training to accelerate training and improve the model's generalization ability, while techniques such as Dropout (e.g., a Dropout rate of 0.3) are combined to effectively prevent overfitting. The training dataset is crucial for building the model, containing electromagnetic interference samples and clean seismic signal samples collected from different earthquake monitoring stations in China and Germany under different electromagnetic environments, labeled manually or semi-automatically. The sample size is at least tens of thousands, for example, containing 50,000 labeled 5-second data segments, including electromagnetic interference events of different intensities and types, as well as clean seismic background noise, ensuring the model's robustness and generalization ability.

[0080] Seventh, the output of the electromagnetic interference identification model 30 provides the probability of electromagnetic interference presence, type classification, intensity estimation, dominant frequency, temporal correlation, and possible spatial azimuth information. The probability of electromagnetic interference presence is a value between 0 and 1, representing the likelihood that the current data segment is affected by electromagnetic interference. Type classification clearly indicates the nature of the interference, such as "power frequency interference" or "GSM communication interference." Intensity estimation provides the quantified intensity of the interference signal, such as dBm or nV / m. The dominant frequency indicates the frequency point where the interference energy is most concentrated. Temporal correlation illustrates the precise temporal correspondence between the interference event and the noise segment in the seismic data. If multiple sensors are deployed, the model can also output the inferred spatial azimuth of the interference source, such as 30 degrees east of north, and the potential distance range. These outputs collectively constitute a comprehensive diagnostic report of electromagnetic interference events.

[0081] Eighth, based on the output of the electromagnetic interference identification model, the electromagnetic interference impact of the seismic physical signal data is assessed and quality is marked, and electromagnetic interference suppression or compensation can be selectively performed. This step is the final step in this method to improve the quality of seismic data.

[0082] First, based on the output probability of the electromagnetic interference identification model, when the probability of a specific type of electromagnetic interference exceeds a preset threshold (e.g., 0.8), the seismic data for that time period is determined to be significantly affected by that type of electromagnetic interference. This threshold can be adjusted according to the actual application scenario and the requirements for data purity.

[0083] Secondly, by combining the intensity estimation, dominant frequency, and temporal correlation of electromagnetic interference (EMI), the percentage decrease in data purity or the reduction in signal-to-noise ratio (SNR) in dB caused by EMI to seismic physical signal data is calculated. The SNR reduction in dB is calculated by comparing the power spectrum or root mean square value of the seismic signal during the interfered period and the interference-free period (or based on a background noise model estimation). For example, SNR_loss_dB = 10*log10(P_signal_clean / P_signal_noisy), where P_signal_clean represents the clean signal power and P_signal_noisy represents the interfered signal power. This provides a quantitative indicator of the impact of EMI.

[0084] Furthermore, standardized quality flags are added to seismic physical signal data segments affected by electromagnetic interference. For example, quality control words defined in the SEED (Standard for the Exchange of Earthquake Data) format are used to identify affected data segments as "affected by electromagnetic interference," "high noise," or "pending review." Quality flags can further include interference type codes (e.g., "EM_POWERLINE," "EM_RF"), intensity levels (e.g., "low," "medium," "high"), the start and end points of the affected period (accurate to milliseconds), and suggested follow-up processing measures (e.g., recommending specific frequency filtering, adaptive noise reduction, or skipping the use of this data segment in critical analyses). These flags ensure that data users can clearly identify and avoid contaminated data, thereby guaranteeing the reliability of subsequent seismological research.

[0085] In a preferred embodiment of the present invention, electromagnetic interference suppression or compensation can be selectively implemented using digital signal processing methods such as adaptive noise cancellation algorithms, wavelet domain denoising, or notch filtering. The adaptive noise cancellation algorithm employs LMS (Least Mean Square) or RLS (Recursive Least Squares) algorithms, using multi-dimensional electromagnetic field data as a reference input signal to subtract relevant electromagnetic interference components from the seismic physical signal data. For example, when 50Hz power frequency interference is identified, the 50Hz component in the electric or magnetic field data can be used as a reference signal input to the LMS adaptive filter. This filter minimizes the output of 50Hz noise in the seismic data by adjusting its weights, thereby achieving dynamic suppression of related interference. The center frequency and bandwidth of the notch filter are dynamically adjusted according to the identified dominant frequency and interference bandwidth. For example, if broadcast signal interference with a dominant frequency of 140MHz is identified, a notch filter with a center frequency of 140MHz and a bandwidth of 1MHz is configured to accurately remove interference in this frequency band while minimizing the impact on the target seismic signal. Wavelet domain denoising utilizes wavelet transform to decompose the signal into multiple scales, identifies and thresholds the noise-related parts of the wavelet coefficients, and then reconstructs the denoised signal through inverse wavelet transform.

[0086] In a specific experimental embodiment, the method of this invention was deployed at a seismic station in Sichuan Province, China, surrounded by high-voltage power lines and mobile communication base stations. The monitoring period was 30 days, during which a total of 120TB of electromagnetic field data and seismic physical signal data were collected. The seismic signal acquisition unit used a Guralp CMG-3ESPC broadband seismograph with a sampling rate of 200 SPS. The multi-dimensional electromagnetic field sensing unit consisted of customized triaxial electric and magnetic field sensors, with the electric field sensor having a response frequency up to 2GHz and the magnetic field sensor up to 100MHz, and a sampling rate of 1GSPS. All sensors achieved nanosecond-level time synchronization via the PTPv2 protocol.

[0087] Example: Identification and Suppression of Power Frequency Interference in High-Voltage Transmission Lines

[0088] Between 2:30 AM and 3:15 AM on a certain day, seismic physical signal data showed abnormally high energy at 50 Hz and its harmonics, with a PSD value 15 dB higher than the background noise. Traditional methods, relying solely on fixed thresholds or a single low-pass filter, cannot effectively distinguish weak seismic signals from strong power frequency interference in this frequency band, or may lead to excessive attenuation of the target seismic signal.

[0089] The method of this invention first extracts frequency domain features from multi-dimensional electromagnetic field data. By performing FFT on the magnetic field data, significant narrowband energy peaks at 50Hz, 100Hz, and 150Hz are identified, with peak power spectral density 30dB higher than surrounding frequency bands. Simultaneously, electric field data also shows a clear energy concentration in the same frequency band. Time-domain analysis reveals that the signal is periodic, consistent with the periodic characteristics of alternating current. Polarization-domain analysis shows that the electromagnetic field in this frequency band mainly exhibits linear polarization characteristics, and the polarization direction is roughly parallel to the direction of the high-voltage transmission line. Spatial-domain analysis (using AoA technology with multiple electromagnetic field sensing units) estimates the azimuth angle of the interference source to be approximately 45 degrees northeast of the station, which matches the actual location of the high-voltage transmission line.

[0090] These multi-dimensional electromagnetic interference feature vectors (including PSD peak, dominant frequency, harmonic components in the frequency domain, periodicity and RMS value in the time domain, linear polarization degree in the polarization domain, and azimuth angle in the spatial domain) and the synchronous seismic signal noise feature vectors (including PSD anomalies, signal-to-noise ratio decrease, and cross-correlation coefficients with electromagnetic field data in the 50Hz band, where the cross-correlation coefficient is as high as 0.85 at 50Hz) are input together into the pre-trained electromagnetic interference identification model.

[0091] The model output shows that the probability of "power line harmonic interference" is as high as 0.98, with an estimated intensity of moderate to high. The dominant frequency is 50Hz, exhibiting extremely strong temporal correlation and a high degree of correlation with the 50Hz noise component in the seismic data. An impact assessment of the seismic data for this period shows a signal-to-noise ratio decrease of approximately 8.5dB. Based on the model's recommendations, this data segment is labeled "affected by electromagnetic interference_high noise," and adaptive noise cancellation (LMS algorithm) is recommended for processing.

[0092] In LMS adaptive noise cancellation processing, the 50Hz frequency and its harmonic components in the multi-dimensional electromagnetic field data are used as reference input to adaptively subtract relevant interference from the contaminated seismic signal. After processing, the PSD value of the seismic signal in the 50Hz band is reduced by approximately 12dB. At the same time, the accuracy of P-wave initial motion interpretation during this period is significantly improved, and the analysis of high-frequency attenuation is also more accurate.

[0093] The method proposed in this invention is universal and adaptable, applicable to various fixed-deployment earthquake monitoring scenarios, including regional seismic networks, deep-well seismic stations, and ocean bottom seismometers. Its core lies in the refined perception and intelligent analysis of the electromagnetic environment, unrestricted by the physical mobility of monitoring equipment and independent of non-seismic signal characteristics such as communication error rates. Therefore, it effectively overcomes the inherent contradictions and limitations of existing technologies in fixed, high-precision earthquake monitoring scenarios. Through the detailed description of the above embodiments, the technical solution of this invention and its significant technical effects are fully demonstrated. It can significantly improve the reliability and accuracy of earthquake monitoring data, providing cleaner and more reliable observational data for earthquake scientific research.

Claims

1. An electromagnetic interference detection method for earthquake monitoring, characterized in that, Includes the following steps: First, at least one multi-dimensional electromagnetic field sensing unit (10) and at least one seismic signal acquisition unit (20) are deployed at the location of the earthquake monitoring equipment or in a predetermined area around it, and the time synchronization accuracy between the multi-dimensional electromagnetic field sensing unit (10) and the seismic signal acquisition unit (20) is better than 1 microsecond. Second, through the multi-dimensional electromagnetic field sensing unit (10), multi-dimensional electromagnetic field data in the region are collected synchronously. The multi-dimensional electromagnetic field data includes electric field strength data, magnetic field strength data, electromagnetic wave polarization data and electromagnetic wave angle of arrival data. The multi-dimensional electromagnetic field sensing unit (10) has wideband response characteristics and can cover the electromagnetic spectrum range from 1Hz to 3GHz. Third, the seismic signal acquisition unit (20) synchronously acquires seismic physical signal data in the area. The seismic physical signal data includes ground motion velocity, ground motion acceleration or ground strain. The seismic signal acquisition unit (20) has high sensitivity and wide dynamic range characteristics and can respond to the seismic signal frequency band from 0.01Hz to 100Hz. Fourth, multi-feature extraction and analysis are performed on the multi-dimensional electromagnetic field data to generate multi-dimensional electromagnetic interference feature vectors. The multi-feature extraction includes time-domain feature extraction, frequency-domain feature extraction, time-frequency-domain feature extraction, polarization-domain feature extraction, and spatial-domain feature extraction. Fifth, noise features are extracted from the earthquake physical signal data to generate an earthquake signal noise feature vector; Sixth, the multi-dimensional electromagnetic interference feature vector and the seismic signal noise feature vector are input into the pre-trained electromagnetic interference identification model (30). The electromagnetic interference identification model (30) is constructed based on machine learning or deep learning algorithms and trained using historical, labeled electromagnetic interference data samples and clean seismic signal data samples to learn the mapping relationship and discrimination boundary between electromagnetic interference features and seismic signal noise features. Seventh, through the electromagnetic interference identification model (30), output the probability of the presence of electromagnetic interference, type classification, intensity estimation, dominant frequency, time correlation and possible spatial azimuth information; Eighth, based on the output of the electromagnetic interference identification model (30), the seismic physical signal data is subjected to electromagnetic interference impact assessment and quality labeling, and electromagnetic interference suppression or compensation processing is selectively performed.

2. The electromagnetic interference detection method for earthquake monitoring according to claim 1, characterized in that, The multi-dimensional electromagnetic field sensing unit (10) mentioned in the first step specifically includes: A triaxial orthogonal electric field sensor array (11) is used to collect the electric field strength data. The triaxial orthogonal electric field sensor array (11) consists of three mutually perpendicular electric dipole antennas. Each electric dipole antenna has a symmetrical arm with a length of 100 mm to 500 mm. Its gain flatness is maintained within ±3dB in the range of 1 Hz to 3 GHz. Its input impedance is 50 ohms and its noise figure is less than 3dB. A triaxial orthogonal magnetic field sensor array (12) is used to collect the magnetic field strength data. The triaxial orthogonal magnetic field sensor array (12) consists of three mutually perpendicular single-turn or multi-turn loop coil antennas. The diameter of each coil antenna is 50 mm to 200 mm. Its induced voltage is proportional to the rate of change of magnetic flux. The response frequency range is 1 Hz to 100 MHz. The sensitivity at 10 kHz is not less than 10 mV / nT. The inherent noise is less than 1 pT / √Hz. A broadband radio frequency receiving front-end (13) is used to receive and digitize the electromagnetic field data. The broadband radio frequency receiving front-end (13) includes a low-noise amplifier, a mixer, an intermediate frequency filter, and a high-speed analog-to-digital converter. The noise figure of the low-noise amplifier is less than 2dB, and the gain is 20dB to 40dB. The high-speed analog-to-digital converter adopts 14-bit or 16-bit resolution and the sampling rate can reach 200MSPS to 5GSPS. And a high-precision time synchronization module (14) is used to realize the time alignment between the electromagnetic field data and the seismic physical signal data. The high-precision time synchronization module (14) adopts a local high-precision crystal oscillator or atomic clock calibrated by a global navigation satellite system receiver, and realizes nanosecond or sub-nanosecond time synchronization with other units through network time protocol or precise time protocol.

3. The electromagnetic interference detection method for earthquake monitoring according to claim 1, characterized in that, The earthquake signal acquisition unit (20) in the third step specifically includes: a three-axis broadband velocity type seismic sensor (21) for acquiring ground motion velocity data, with a frequency response range of 0.01Hz to 100Hz, a sensitivity of not less than 2000V / (m / s), and a self-noise level lower than that of the seismic low-noise model; A high-resolution analog-to-digital converter (22) is used to digitize the seismic physical signal data. The high-resolution analog-to-digital converter (22) is a 24-bit or 26-bit Σ-Δ analog-to-digital converter with an effective number of bits not less than 22 bits, a sampling rate configured as 100SPS, 200SPS or 500SPS, and a total harmonic distortion plus noise of less than -110dB. And a data preprocessing module (23) for filtering and gain adjustment of the digitized seismic physical signal data. The data preprocessing module (23) integrates a hardware digital filter that can realize low-pass, high-pass, band-pass or band-stop filtering functions. The band-stop filter is configured with a center frequency of 50Hz or 60Hz and a stopband attenuation depth greater than 80dB.

4. The electromagnetic interference detection method for earthquake monitoring according to claim 1, characterized in that, The fourth step, multi-feature extraction and analysis, specifically includes: The multidimensional electromagnetic field data is subjected to time-domain feature extraction, which includes calculating the instantaneous power, peak amplitude, root mean square amplitude, zero-crossing rate, pulse width, rise time, and fall time characteristics of the multidimensional electromagnetic field data. The instantaneous power is obtained by calculating the square of the instantaneous amplitude of the electric or magnetic field; the peak amplitude represents the maximum absolute value of the signal within a specified time window; the root mean square amplitude represents the effective value of the signal power within a specified time window; the zero-crossing rate represents the number of times the signal crosses the zero level within a specified time window; and the pulse width, rise time, and fall time are used to characterize the duration and waveform steepness of the instantaneous pulse-type interference.

5. The electromagnetic interference detection method for earthquake monitoring according to claim 4, characterized in that, The fourth step, multi-feature extraction and analysis, specifically includes: Frequency domain feature extraction is performed on the multi-dimensional electromagnetic field data. The frequency domain feature extraction includes performing fast Fourier transform or Welch spectrum estimation on the multi-dimensional electromagnetic field data to obtain power spectral density, dominant frequency, harmonic components, bandwidth, spectral centroid and spectral peak-sidelobe ratio features. The calculation window length of the power spectral density curve is set to 2 to 10 seconds, and the overlap rate is 50% to 75%; the dominant frequency is the frequency point corresponding to the peak of the power spectral density; the harmonic component identification is achieved by detecting the energy peak that appears at an integer multiple of the fundamental frequency; the bandwidth is the frequency range when the signal power drops to half or one-tenth of the peak power.

6. The electromagnetic interference detection method for earthquake monitoring according to claim 5, characterized in that, The fourth step, multi-feature extraction and analysis, specifically includes: Time-frequency domain feature extraction is performed on the multi-dimensional electromagnetic field data. The time-frequency domain feature extraction includes performing short-time Fourier transform or wavelet transform on the multi-dimensional electromagnetic field data to obtain time-varying spectrum features, instantaneous frequency and energy concentration features. The short-time Fourier transform uses a Hanning window or a Gaussian window with a window length configured from 20 milliseconds to 500 milliseconds and a window overlap rate of 50%; the wavelet transform uses a Morlet wavelet or a Daubechies wavelet basis. And to extract polarization domain features from the multi-dimensional electromagnetic field data, the polarization domain feature extraction is based on the output of the triaxial orthogonal electric field sensor array (11) or the triaxial orthogonal magnetic field sensor array (12), to calculate the Stokes parameters of the electromagnetic wave in order to obtain the linear polarization degree, circular polarization degree and polarization angle features; The Stokes parameters are obtained by calculating the cross-correlation function and the autocorrelation function.

7. The electromagnetic interference detection method for earthquake monitoring according to claim 6, characterized in that, The fourth step, multi-feature extraction and analysis, specifically includes: When multiple multi-dimensional electromagnetic field sensing units (10) are deployed, multi-station time difference positioning technology, angle of arrival technology or electromagnetic beamforming technology are used to locate and track electromagnetic interference sources to obtain the spatial azimuth, elevation, distance and movement trajectory characteristics of the interference sources; wherein, the angle of arrival technology adopts the MUSIC algorithm or ESPRIT algorithm to estimate the arrival direction of the signal by analyzing the phase difference of the signals between multiple receiving antennas; the multi-station time difference positioning technology calculates the location of the interference source by measuring the time difference of different sensors arriving at the same interference signal; the multi-station time difference positioning technology requires at least three synchronous sensors for two-dimensional positioning or at least four synchronous sensors for three-dimensional positioning.

8. The electromagnetic interference detection method for earthquake monitoring according to claim 1, characterized in that, The noise feature extraction in the fifth step specifically includes: The seismic physical signal data is preprocessed by detrending, DC removal, and removal of instrument response. Within a specific time window, power spectral density analysis is performed on the preprocessed seismic physical signal data to identify anomalous energy concentration frequency bands that exceed the background seismic noise spectrum; the energy ratio or signal-to-noise ratio decrease of the seismic physical signal data within a preset electromagnetic interference sensitive frequency band is calculated. Analyze the instantaneous pulse or spike noise in the earthquake physical signal data and distinguish it from typical earthquake signal waveforms; The correlation between noise components in the seismic physical signal data and electromagnetic interference characteristics of specific frequency bands or types in the multidimensional electromagnetic field data is assessed through cross-correlation analysis or coherence analysis. Specifically, the detrending is achieved by fitting and removing linear or polynomial trend terms from the signal using the least squares method; the DC removal is achieved using a high-pass filter with a cutoff frequency of 0.001 Hz; and the removal of instrument response is achieved by deconvolution or frequency domain division to remove the influence of the seismograph's transfer function. The power spectral density analysis employs either the Fast Fourier Transform or Welch spectral estimation method, with the window length matching the electromagnetic field data analysis window. The sensitive frequency bands are identified by analyzing historical electromagnetic environment data, which reveals common electromagnetic interference frequencies that affect seismic data. The instantaneous pulse is detected using a method based on an amplitude threshold or based on a specific waveform template matching. The cross-correlation analysis calculates the time delay and correlation coefficient between seismic noise and electromagnetic interference signals; the coherence analysis calculates the linear correlation between the two in the frequency domain.

9. The electromagnetic interference detection method for earthquake monitoring according to claim 1, characterized in that, The electromagnetic interference identification model (30) in the sixth step specifically includes: an input layer (31), which is used to receive the multi-dimensional electromagnetic interference feature vector and the seismic signal noise feature vector, and merge the two into a high-dimensional input vector; The feature fusion layer (32) is used to perform preliminary feature learning and dimensionality reduction on the merged high-dimensional input vector using a convolutional neural network or a multilayer perceptron, so as to extract higher-level abstract features; The sequence processing layer (33) is used to capture the dependencies of the feature vectors on the time series by employing a long short-term memory network or a gated recurrent unit recurrent neural network structure. The classification decision layer (34) is used to use the Softmax function as the activation function of the output layer to output the probability distribution of different electromagnetic interference types and identify the state of pure seismic signals without electromagnetic interference. The training optimization step is used to train the model using the gradient descent optimization algorithm and the cross-entropy loss function. During the training process, batch normalization and Dropout techniques are used to prevent overfitting. The convolutional neural network layer contains multiple convolutional kernels, uses ReLU as the activation function, and can be connected to pooling layers for dimensionality reduction; the multilayer perceptron layer contains 2 to 5 fully connected layers, with 64 to 512 neurons in each layer. The Long Short-Term Memory (LSTM) network layer comprises 1 to 3 layers, each containing 64 to 256 units; The classification decision layer is a fully connected layer, and the number of its neurons is equal to the number of predefined electromagnetic interference types plus one. The training dataset contains electromagnetic interference data samples and clean seismic signal data samples collected from different earthquake monitoring stations under different electromagnetic environments, which have been manually or semi-automatically labeled. The number of samples is at least tens of thousands.

10. The electromagnetic interference detection method for earthquake monitoring according to claim 1, characterized in that, The electromagnetic interference impact assessment and quality labeling in step eight specifically includes: According to the output probability of the electromagnetic interference identification model (30), when the probability of a specific type of electromagnetic interference exceeds a preset threshold, it is determined that the seismic data in that time period is affected by that type of electromagnetic interference. Based on the intensity estimation, dominant frequency, and time correlation of the electromagnetic interference, calculate the percentage decrease in data purity or the decibel decrease in signal-to-noise ratio caused by the electromagnetic interference to the seismic physical signal data. For the seismic physical signal data segments affected by electromagnetic interference, add standardized quality markers or flags; The electromagnetic interference suppression or compensation process is selectively implemented through adaptive noise cancellation algorithms, wavelet domain denoising, or notch filtering digital signal processing methods. The reduction in signal-to-noise ratio in decibels is calculated by comparing the power spectrum or root mean square value of the seismic signal during the period of interference and the period of no interference. The quality label includes the interference type code, intensity level, start and end points of the affected period, and suggested follow-up measures; The adaptive noise cancellation algorithm uses the LMS algorithm or the RLS algorithm, takes the multi-dimensional electromagnetic field data as a reference input, and subtracts the relevant electromagnetic interference components from the seismic physical signal data. The center frequency and bandwidth of the notch filter are dynamically adjusted according to the identified dominant frequency and interference bandwidth.

Images