A wide-area electromagnetic data processing method, electronic equipment and storage medium

Abstract

The present application belongs to the technical field of data processing, and specifically provides a wide-area electromagnetic data processing method, an electronic device and a storage medium, which comprises the following steps: acquiring original time-domain data collected by a wide-area electromagnetic method; segmenting and denoising the original time-domain data to obtain segmented and denoised time-domain data; performing segmented discrete Fourier transform analysis on the segmented and denoised time-domain data to obtain each segment spectrum; calculating the signal-to-noise ratio of a target frequency point in each segment spectrum and screening out a high-quality spectrum segment; performing spectrum estimation on the amplitude of each target frequency point in the screened high-quality spectrum segment to obtain spectrum estimation results of multiple robust statistical methods; performing bias and consistency analysis on the spectrum estimation results of the multiple robust statistical methods, respectively, and outputting a final spectrum. The present application extracts main periodic signal components through segmented analysis, thereby more effectively suppressing the interference of non-periodic noise on spectrum estimation and improving the extraction stability of weak signals in a strong noise background.

Description

Technical Field

[0001] This invention belongs to the field of data processing technology and relates to a wide-area electromagnetic data processing method, electronic device and storage medium. Background Technology

[0002] Wide-area electromagnetic sounding (CSEM) is a method for detecting underground electrical structures using artificially controlled electromagnetic field sources. In this method, a transmitter typically excites electromagnetic signals at multiple discrete frequency points over a wide frequency band, while a receiver records the received time-domain data at the surface. Due to the controlled excitation by an artificial source, CSEM offers a higher signal-to-noise ratio for shallow and deep subsurface exploration compared to natural field magnetotelluric (MT) methods. However, various strong interference noises still exist in the field environment, posing challenges to data processing. Current techniques typically involve directly performing a Fast Fourier Transform (FFT) on the received entire time-domain data to extract the amplitude and phase corresponding to each transmission frequency point. While the full-segment FFT spectrum analysis method is simple to implement, it has significant shortcomings in high-noise environments: unstable spectrum estimation and poor repeatability are common problems. Although multi-period superposition measurements are often used in practical observations to enhance effective data and average random noise, if some periods are contaminated by sudden interference, simple linear superposition will still introduce these anomalous noises into the spectrum estimation, making it difficult to obtain stable and reliable results.

[0003] Noise sources in wide-area electromagnetic reception signals can generally be divided into two main categories: periodic interference and aperiodic interference. For periodic noise at a fixed frequency (such as 50 / 60 Hz power grid interference or instrument periodic drift noise), it can usually be effectively suppressed in the frequency domain using methods such as narrowband notch filtering. However, aperiodic noise often exhibits characteristics of being sudden, short-duration, and wide-bandwidth, with its energy distribution being unconcentrated and unpredictable. In the entire FFT processing, the presence of this type of aperiodic noise leads to extremely unstable amplitude extraction at the target frequency. For example, if a high-amplitude transient pulse interference is superimposed on the received record at a certain moment, it will appear in the FFT spectrum in a wide-bandwidth form, thus obscuring the useful signal peaks at adjacent frequencies, resulting in severe distortion of the spectral amplitude estimation. Because aperiodic noise is highly non-stationary and its spectrum may overlap with the signal band, traditional frequency domain filtering struggles to remove it without damaging the useful signal. Therefore, under strong noise interference conditions, the anti-interference capability and reliability of existing whole-spectrum extraction methods are significantly limited.

[0004] On the other hand, based on the signal characteristics of the active source electromagnetic method, the received time-domain data can be modeled as the sum of a periodic component (PC) and a non-periodic component (NPC). The periodic component includes a stable signal component generated by periodic excitation of the artificial source and periodic background noise superimposed on it, while the non-periodic component consists of randomly occurring transient interference and background Gaussian white noise. Summary of the Invention

[0005] Wide-area electromagnetic methods typically employ a controllable source to emit a wide-bandwidth continuous signal, acquiring a response containing multiple frequency components at the receiving end. The effective signal persists throughout the observation period, while noise appears randomly or intermittently. This invention aims to provide an efficient noise reduction method for wide-area electromagnetic sounding data. By utilizing the characteristics of wide-area electromagnetic methods, principal component analysis is used to extract robust principal components from the data, thereby separating and removing random noise, improving signal quality, and overcoming the shortcomings of existing technologies in handling complex noise.

[0006] This invention provides a wide-area electromagnetic data processing method, comprising the following steps:

[0007] Step 1: Acquire raw time-domain data obtained by wide-area electromagnetic method ;

[0008] Step 2: Process the raw time-domain data The data is segmented, and PCA is used to analyze the original time-domain data after segmentation. Denoising is performed to obtain segmented denoised time-domain data;

[0009] Step 3: Perform piecewise discrete Fourier transform analysis on the piecewise denoised time-domain data to obtain the spectrum of each segment;

[0010] Step 4: Calculate the signal-to-noise ratio (SNR) of the target frequency in each spectrum segment, sort the spectrum segments according to their SNR, and select high-quality spectrum segments.

[0011] Step 5: Use multiple robust statistical methods to estimate the amplitude of each target frequency point in the selected high-quality spectrum bands, and obtain the spectrum estimation results of multiple robust statistical methods;

[0012] Step 6: Perform bias and consistency analysis on the spectrum estimation results of various robust statistical methods, and output the final spectrum.

[0013] Furthermore, the specific process of step two is as follows:

[0014] Original time domain data Divided by length Segment time-domain data, assuming each segment of time-domain data contains There are consecutive sampling points, and there is a certain distance between two adjacent time-domain data segments. The proportions of the sample sizes overlap;

[0015] All time-domain data are used to construct corresponding column vectors and row vectors, and all column vectors and row vectors are stacked into a data matrix. ;

[0016] For data matrix Perform principal component analysis and based on the data matrix The strategy for selecting the number of principal components to retain in the time-domain data based on the principal component analysis results determines the amount of time-domain data to be retained. One principal component;

[0017] Construct the truncated singular value matrix and utilize Reconstruct the data matrix from the principal components. ;

[0018] Data matrix By restoring the data to its time series form, the reconstructed time-domain data can be obtained. For reconstructing time-domain data For the segmented reconstructed time-domain data with overlapping non-Chinese and African components, each segment of the reconstructed time-domain data is directly concatenated sequentially to obtain the segmented denoised time-domain data; for the reconstructed time-domain data... The time-domain data is reconstructed by segmenting overlapping segments. The average of the signals in the overlapping segments is taken and smoothly connected to obtain the segmented denoised time-domain data.

[0019] Furthermore, an index matrix is ​​used to analyze the original time-domain data. Perform segmentation processing to obtain Time-domain data;

[0020] The index matrix is ;

[0021] ;

[0022] in, The starting position, and , segment number and ; Step size; To round down; The total number of segments, and ; This is the total length of the window; This represents the total number of sampling points for the original time-domain data.

[0023] Furthermore, regarding the data matrix The specific method for performing principal component analysis is as follows:

[0024] For data matrix Centralized processing is performed to obtain a centralized data matrix. ;

[0025] For centralized data matrix Perform principal component analysis to decompose the left singular vector matrix. Singular value matrix and right singular vector matrix .

[0026] Furthermore, based on the data matrix The specific method for selecting the number of principal components to retain in the time-domain data based on the principal component analysis results is as follows:

[0027] By analyzing the singular value spectrum, the boundary between signal and noise is determined, thus identifying whether to retain the time-domain data using either a directly specified method or an automatically determined method based on the cumulative variance contribution rate. One principal component.

[0028] Furthermore, the specific process of step three is as follows:

[0029] by The length is used to further divide the segmented, denoised time-domain data into segments. Segmented data in the time domain;

[0030] The discrete Fourier transform of each segment of time-domain data is calculated by multiplying it by the Hanning window function to obtain the segmented spectrum.

[0031] Furthermore, obtain a single-segment spectrum. The specific process is as follows:

[0032] Set the Hanning window function for single-segment time-domain segmented data The windowed time-domain data is obtained. ;

[0033] Based on windowed time-domain data Perform Discrete Fourier Transform calculations to obtain the Discrete Fourier Transform results. ;

[0034] Based on the results of Discrete Fourier Transform Calculate the positive frequency component of the one-sided amplitude spectrum. ; and the recovery coefficient based on the Hanning window For the positive frequency portion of the single-sided amplitude spectrum To perform recovery and obtain the extent of recovery. That is, the true physical amplitude of the Y-axis of a single-segment time-domain spectrum;

[0035] Calculate the frequency axis That is, to calculate the true physical frequency of the X-axis of a single-segment time-domain spectrum; and to set the discrete index points This is mapped to a specific transmission frequency, resulting in a frequency axis. The actual spectrum;

[0036] Based on the magnitude of recovery and frequency axis Calculate the first Segment data on the frequency axis Complete spectrum data at [location] This yields a single-segment spectrum.

[0037] Furthermore, the specific process of step four is as follows:

[0038] Calculate the spectrum of each segment The ratio of signal strength to background noise at the main transmission frequency point is recorded as the signal-to-noise ratio (SNR) of that frequency band. ;

[0039] Signal-to-noise ratio of each spectrum Sort the spectrum from high to low, and set the percentage threshold for signal-to-noise ratio (SNR) filtering to P%. Retain the SNR values ​​of each spectrum in the top P% as high-quality spectrum segments.

[0040] Furthermore, the robust statistical methods in step five include at least the following three:

[0041] The first approach: Define a set of high-quality spectrum bands. Includes Each effective spectrum segment is located at the target frequency. The amplitude at that point is denoted as ; Calculate the high-quality spectrum band at the target frequency amplitude at simple average And based on simple average Spectrum estimation is performed on each effective spectrum segment within the high-quality spectrum band.

[0042] The second approach: Define a set of high-quality spectrum bands. Includes Each effective spectrum segment is located at the target frequency. The amplitude at that point is denoted as ; Calculate the median of each effective spectrum segment. And based on the median Spectrum estimation is performed on each effective spectrum segment within the high-quality spectrum band.

[0043] The third method involves discarding a certain proportion of the largest and smallest samples from the high-quality spectrum segments, and then averaging the remaining high-quality spectrum segments to obtain the truncated mean. ; and based on the truncated mean Spectrum estimation is performed on each effective spectrum segment within the high-quality spectrum band.

[0044] The fourth method: Obtain the IQR mean of the remaining effective spectrum segments using the mean estimation method based on the interquartile range. And based on the IQR mean Spectrum estimation is performed on each effective spectrum segment within the high-quality spectrum band.

[0045] The fifth method applies the Huber robust estimation concept, assigning different weights to each effective spectrum segment and then averaging the results to obtain the Huber weighted mean. And based on Huber's weighted average Spectrum estimation is performed on each effective spectrum segment within the high-quality spectrum band.

[0046] Furthermore, the specific process of step six is ​​as follows:

[0047] Median of the result estimated by median As a reference benchmark, the deviations and relative errors of the first, third to fifth estimation methods relative to the reference benchmark are compared;

[0048] set up Represented as any estimation method, the absolute deviation of that estimation method at each target frequency point is calculated. and relative error Distribution, determine whether this estimation method is similar to the median. To assess the consistency and robustness of the estimation method, the bias analysis results were obtained.

[0049] Based on the deviation analysis results of all estimation methods, the estimation method that is most stable for estimating the spectral amplitude is selected, and its estimation result is output as the final spectrum.

[0050] As a further aspect of the present invention, the present invention also provides an electronic device, including a memory, one or more processes, and one or more programs stored in the memory, said one or more programs including instructions for performing the wide-area electromagnetic data processing method as described above.

[0051] As a further aspect of the present invention, the present invention also provides a storage medium including one or more programs executable by one or more processors of an electronic device, said one or more programs including instructions for performing the wide-area electromagnetic data processing method as described above.

[0052] Compared with the prior art, the present invention has the following beneficial effects:

[0053] (1) Based on the signal feature points of the active source electromagnetic method, the present invention models the received time-domain data as the sum of periodic components (PC) and non-periodic components (NPC), which serves as the theoretical basis for the segmented noise reduction proposed in this invention: Since the periodic component appears repeatedly in each time period and its energy is relatively concentrated, while the non-periodic component exhibits statistical inconsistency in different time periods, this decomposable characteristic of the signal can be utilized to extract the main periodic signal components through segmented analysis, thereby more effectively suppressing the interference of non-periodic noise on spectrum estimation and improving the extraction stability of weak signals in the background of strong noise.

[0054] (2) This invention addresses the processing requirements of wide-area electromagnetic received signals under high background noise environments by combining principal component analysis (PCA) segmented noise reduction, spectral signal-to-noise ratio screening, and various robust statistical estimation methods to achieve highly robust extraction of the signal spectrum. The purpose of this invention is to improve the anti-interference performance and accuracy of spectral amplitude estimation of wide-area electromagnetic received data under strong interference noise backgrounds, so that weak geological response signals can be stably identified and quantified even under complex noise conditions.

[0055] (3) The outstanding feature of this invention is that it combines PCA segmented noise reduction with a variety of robust statistical estimation methods, thereby achieving highly robust extraction of the spectrum of wide-area electromagnetic received data and significantly enhancing the detectability and estimation accuracy of weak signals in the context of strong noise.

[0056] (4) The robust spectrum extraction scheme proposed in this invention is expected to play an important role in applications such as mineral resource exploration and geothermal energy exploration that require extraction of artificial source electromagnetic signals under strong interference background, thereby improving the detection accuracy and reliability of related electromagnetic sounding technologies.

[0057] In addition to the objectives, features, and advantages described above, the present invention has other objectives, features, and advantages. The invention will now be described in further detail with reference to the figures. Attached Figure Description

[0058] The accompanying drawings, which form part of this application, are used to provide a further understanding of the invention. The illustrative embodiments of the invention and their descriptions are used to explain the invention and do not constitute an undue limitation of the invention. In the drawings:

[0059] Figure 1 This is a flowchart illustrating a wide-area electromagnetic data processing method according to Embodiment 1 of the present invention;

[0060] Figure 2This is a schematic diagram of the PCA segmented noise reduction and signal reconstruction process in Embodiment 1 of the present invention;

[0061] Figure 3(a) is a schematic diagram of the spectrum of the original data in Embodiment 1 of the present invention;

[0062] Figure 3(b) is a schematic diagram of the spectrum of PCA noise reduction data in Embodiment 1 of the present invention;

[0063] Figure 4 This is a bar chart illustrating the errors of different statistical methods in Embodiment 1 of the present invention. Detailed Implementation

[0064] To make the above-mentioned objects, features, and advantages of the present invention clearer and easier to understand, the specific embodiments of the present invention will be described in detail below with reference to the accompanying drawings. It should be noted that the accompanying drawings of the present invention are all in a simplified form and use non-precise proportions, and are only used to facilitate and clearly illustrate the implementation of the present invention.

[0065] Example 1:

[0066] See Figure 1 and Figure 2 As shown, the present invention provides a wide-area electromagnetic data processing method, which relates to the field of geophysical exploration data processing technology, specifically to the field of active source wide-field electromagnetic sounding (WFEM) data processing technology. More specifically, it relates to a method combining principal component analysis (PCA) noise reduction and robust spectrum estimation to improve the anti-interference capability and spectrum analysis accuracy of wide-area electromagnetic received signals in high-noise environments. Furthermore, the present invention is also applicable to anti-interference processing of artificial source frequency domain electromagnetic detection data such as controlled-source electromagnetic sounding (CSEM).

[0067] Specifically, the wide-area electromagnetic data processing method includes the following steps:

[0068] S1. Time-domain data structure modeling and data reading;

[0069] The specific process is as follows:

[0070] S1.1 Obtaining raw time-domain data acquired by the wide-area electromagnetic method The raw time-domain data To receive the electric / magnetic field data recorded by the receiving coil;

[0071] S1.2 Preprocess the time-domain data to obtain the one-dimensional time-domain data to be processed;

[0072] S1.3, Based on the decomposition model, the original time-domain data is... Decomposed into periodic components Aperiodic components the sum; that is, ;

[0073] in, Represented as raw time-domain data The data contains useful periodic data and periodic noise, which correspond to signal components and periodic noise that repeat stably in time (interference of known frequency, such as power frequency and its harmonics). Digital filters can be preferentially applied to remove these specific frequency periodic noises. Represented as raw time-domain data The data contains a portion of non-periodic interference noise, which corresponds to interference components that appear instantaneously and whose amplitude is difficult to predict.

[0074] By employing a decomposition model, the original time-domain data is... Decomposed into periodic components Aperiodic components This is to facilitate subsequent processing, targeting the periodic components. Aperiodic components Measures are taken for these two different types of components to achieve the goals of noise reduction and robust estimation.

[0075] S2, PCA (Principal Component Analysis) segmented noise reduction;

[0076] The specific process is as follows:

[0077] S2.1, For the original time domain data Segmentation is performed to construct a multi-sample matrix to extract the main feature structure of the time-domain data; specifically:

[0078] Original time domain data Classified by fixed length Segment time-domain data, assuming each segment of time-domain data contains There are 10 consecutive sampling points, and there is a certain percentage of sample overlap between two adjacent time-domain data segments (the overlap percentage is denoted as 1 / 2). ), to increase the sample size;

[0079] The first The time-domain data forms a column vector and row vector (in the formula) Represents the transpose symbol, which will transpose the first... Column vector of time-domain data Transpose to row vector Ensure the first Each row of the time-domain data matrix corresponds to a time-domain data segment sample, and the row vector... belong 3D real vector space ( (representing real numbers);

[0080] No. Column vector of time-domain data The expression is as follows:

[0081] ;

[0082] The column vectors and row vectors formed by all time-domain data are stacked vertically to form... 3D data matrix ;

[0083] Data Matrix The expression is as follows:

[0084] ;

[0085] Preferably, an index matrix is ​​used to analyze the original time-domain data. Perform segmentation processing to obtain Time-domain data.

[0086] Further preferably, let the index matrix of any segment of time-domain data be... ;

[0087] in, The starting position, and ; Let be the step size, and ; To round down; The total number of segments, and ; This represents the total number of sampling points for the original time-domain data. Specifically, the step size... The setting rule is: use the total window length. Subtract the number of overlapping samples (i.e., the overlap ratio) Multiply by the total window length (and rounded down), thus obtaining the actual number of points the window needs to move forward each time. Outer layer This is to set a lower limit to ensure that even with extremely high overlap rates, the step size is at least one sampling point, thus avoiding program dead loops.

[0088] S2.2, Data Matrix Perform principal component analysis; the specific method is as follows:

[0089] First, the data matrix Centralized processing is performed to obtain a centralized data matrix. To eliminate the effects of data shifting and ensure that principal component analysis is based on the true variation of the data rather than simple shifting;

[0090] Then, the centralized data matrix Perform principal component analysis (PCA) to center the data matrix. Equivalent to singular value decomposition into a left singular vector matrix Singular value matrix and right singular vector matrix .

[0091] Preferred, centralized data matrix The expression is as follows:

[0092] ;

[0093] ;

[0094] in, The sample mean vector For the first A column vector composed of time-domain data.

[0095] Preferably, for centralized data matrices The expression for performing PCA decomposition is equivalent to singular value decomposition:

[0096] ;

[0097] in, The left singular vector matrix (principal component score); It is a singular value matrix. diagonal elements Singular values ​​( The total number of diagonal elements. ,and ); It is a right singular vector matrix (principal component loading).

[0098] S2.3 Principal component selection strategy;

[0099] By analyzing the singular value spectrum to determine the boundary between signal and noise, two strategies are employed to select and retain time-domain data. Principal components:

[0100] Strategy 1 (Direct Specification): Directly specify the time-domain data to be retained. One principal component.

[0101] Strategy 2 (Variance Threshold): Automatically determine the time-domain data to be retained based on the cumulative variance contribution rate. One principal component;

[0102] ;

[0103] in, For the preset contribution rate threshold (e.g.) ), To ensure the selection of the minimum variance contribution rate that meets the calculation requirements The value is used to reconstruct the signal while preserving its main energy. Principal component number.

[0104] S2.4, Truncation Reconstruction and Feature Extraction;

[0105] Construct the truncated singular value matrix and utilize Reconstruct the data matrix from the principal components. ;

[0106] Data Matrix The expression is as follows:

[0107] ;

[0108] ;

[0109] Equivalent, data matrix The way to express PCA output items is as follows:

[0110] ;

[0111] in, Principal component scores, Principal component loadings. Specifically, the data matrix. The construction method is as follows: after performing singular value decomposition (SVD) on the matrix, only the first few elements are extracted. The most important principal component (i.e., the left-singular matrix) The former Columns, singular value matrix The former Square matrix of order and right singular matrix The former Multiply the columns together, and finally add back the mean from the centering process. This allows us to reconstruct the signal after noise removal. is a truncated singular value matrix.

[0112] Data Matrix The method for constructing PCA output items is as follows: that is, directly using the previous... The score matrix of each principal component (Score, Multiply by the corresponding load matrix (Coefficient, The transpose of ).

[0113] Reconstructed data matrix Equivalent to only the main part of the original signal The approximation formed by the principal components contains the common main signal features in each segment, while filtering out the noise-dominant components, thereby achieving the purpose of extracting the principal components (PC) and suppressing the non-principal components (NPC).

[0114] S2.5, Data Matrix By restoring the data to its time series form, the reconstructed time-domain data can be obtained. For reconstructing time-domain data For the segmented reconstructed time-domain data with overlapping non-Chinese and African components, each segment of the reconstructed time-domain data is directly concatenated sequentially to obtain the segmented denoised time-domain data; for the reconstructed time-domain data... The time-domain data is reconstructed by segmenting overlapping segments. The average of the signals in the overlapping segments is taken and smoothly connected to obtain the segmented denoised time-domain data.

[0115] Preferably, reconstruct time-domain data The expression is as follows:

[0116] ;

[0117] ;

[0118] ;

[0119] ;

[0120] in, For the first The starting position of the segment, and ; This is the Hanning window function.

[0121] To more fully extract the stable portion of the periodic signal, Each consecutive sampling point is set as the segment length, selected as the signal length of the fundamental period of the transmitted signal or an integer multiple thereof, so that each segment contains one or several complete periods of useful signal, thereby improving the effectiveness of PCA in extracting periodic components. In this segmented PCA noise reduction process, since the repetitive structure of the signal is consistently captured by multiple segments, while random noise exhibits differences in different segments, PCA can separate the former as the principal component, achieving the technical effect of effectively extracting signal features and suppressing noise interference. Its mathematical essence can be summarized as follows:

[0122] ;

[0123] in, For segmented operations, for Principal component projection, This is an overlap-average reconstruction operation.

[0124] S3. Perform piecewise discrete Fourier transform analysis on the piecewise denoised time-domain data to obtain the spectrum of each segment;

[0125] The specific method is as follows: Perform segmented spectrum calculation on the time-domain data after PCA noise reduction; then divide the noise-filtered time-domain data again into segments according to a set length. The time-domain segmented data (the segment length can be the same as in step S2.1; typically, each segment overlaps by 50% to increase the number of independent spectrum estimates) is multiplied by a window function (e.g., a Hanning window to reduce spectral leakage) and its discrete Fourier transform (FFT) is calculated to obtain the spectrum of each segment.

[0126] Specifically, the calculation process for a single-segment spectrum includes:

[0127] (1) For a length of For time-domain segmented data, set its Hanning window function. ;

[0128] Hanning window function The expression is as follows:

[0129] ;

[0130] in, For the length of the window function, This is the index of the discrete sampling points within the window function, and .

[0131] (2) Based on Hanning window function Windowing is applied to each segment of the time-domain data to obtain the windowed time-domain data. ;

[0132] Windowed time-domain data The expression is as follows:

[0133] ;

[0134] (3) Calculate the Discrete Fourier Transform and obtain the Discrete Fourier Transform result. ;

[0135] Discrete Fourier Transform Results The expression is as follows:

[0136] ;

[0137] in, Let be the number of discrete points, and .

[0138] (4) Calculate the recovery of the amplitude spectrum To compensate for the energy loss caused by the window function, the recovery of the amplitude spectrum is calculated. It needs to be multiplied by the recovery factor;

[0139] Amplitude spectrum recovery The expression is as follows:

[0140] ;

[0141] ;

[0142] in, Let be the recovery coefficient of the Hanning window, and ; This represents the positive frequency portion of the single-sided amplitude spectrum.

[0143] Amplitude spectrum recovery calculates the Y-axis of the spectrum (the true physical amplitude, millivolts mV); it converts the complex modulus calculated by FFT into the actual voltage / electromagnetic field strength.

[0144] In this invention, the core function of amplitude spectrum recovery is to restore the true physical amplitude of the signal. Because a Hanning window was applied to the time-domain data in the preceding steps, the window function suppresses the signal to near zero at both ends, inevitably causing a decrease in the total signal energy (energy loss). Without multiplying by the recovery coefficient... The amplitude calculated by FFT will be much smaller than the amplitude of the actual physical signal.

[0145] (5) Calculate the frequency axis and discrete index points This is mapped to a specific transmission frequency, resulting in a frequency axis. The actual spectrum;

[0146] Frequency axis The expression is as follows:

[0147] ;

[0148] in, Sampling rate, This represents the discrete frequency index of the positive frequency portion of the Discrete Fourier Transform, and its value range is... .

[0149] The frequency axis calculates the X-axis of the spectrum (the actual physical frequency, such as Hertz Hz); it sets the discrete index points of the FFT. It is mapped to a specific transmission frequency.

[0150] This invention obtains the true physical amplitude (i.e., the vertical axis of the spectrum) after eliminating the energy attenuation of the window function through the above-described amplitude spectrum recovery; and by mapping the frequency axis, the discrete points are... This is converted into actual physical frequencies (i.e., the horizontal axis of the spectrum). These two steps combined complete the conversion from a purely mathematical sequence to a spectrum with real physical meaning, allowing the response to be directly extracted by referring to the transmission frequency.

[0151] (6) Calculation of a single-segment spectrum:

[0152] Calculate the frequency points using the real and imaginary parts of the FFT result. Phase at ;

[0153] phase The expression is as follows:

[0154] ;

[0155] The magnitude of recovery in step (4) With frequency axis Calculate the first Segment data at physical frequency Complete spectrum data at [location] This yields a single-segment spectrum.

[0156] Spectrum data The expression is as follows:

[0157] .

[0158] In this way, each segment will yield a series of frequency points (including the transmission frequency and its harmonics) with amplitude and phase values. Due to PCA preprocessing, most of the random noise in the signal has been attenuated, thus significantly improving the signal-to-noise ratio near the target frequency in each segment's FFT spectrum compared to the original signal. However, the residual noise level may still vary across different time periods, leading to abnormally high or low amplitude values ​​at individual frequency points in certain segments. Therefore, further quality assessment and screening of the spectral results for each segment are necessary. It should be noted that traditional spectrum estimation often uses the Welch method to perform arithmetic averaging of multiple spectrum segments to reduce noise fluctuations. This invention, however, chooses to introduce robust statistical estimation in subsequent steps instead of simple averaging, thereby more effectively addressing the interference of a few abnormal segments on the results.

[0159] S4. Calculate the signal-to-noise ratio (SNR) of the target frequency in each spectrum segment, sort the spectrum segments according to the SNR, and select high-quality spectrum segments.

[0160] The specific process is as follows:

[0161] S4.1 For each segment of the spectrum Calculate the ratio of signal strength to background noise at the main transmission frequency points (e.g., fundamental frequency and important harmonic frequencies), and record this ratio as the signal-to-noise ratio of that frequency band. ;

[0162] Signal-to-noise ratio The expression is as follows:

[0163] .

[0164] Signal strength It can be defined as the spectral amplitude of the target frequency point, with background noise. It can be estimated by the root mean square amplitude or local noise floor value of the adjacent no-signal frequency band.

[0165] S4.2, Calculate the signal-to-noise ratio of each spectrum. The spectrum segments are sorted from highest to lowest signal-to-noise ratio (SNR) and a certain percentage of segments with the highest SNR are selected as high-quality spectrum segments. Assuming a SNR threshold of P%, several segments with SNR values ​​in the top P% are retained (e.g., half of the high-quality segments are retained when P=50%). Low SNR segments are discarded because their spectral values ​​are not representative due to severe noise contamination and are not included in subsequent statistical estimation. This SNR-based screening step eliminates outlier segments with significant spectral distortion, retaining only data from periods where the signal is dominant and the spectral values ​​are reliable, thus improving the robustness of spectrum estimation.

[0166] S5. For the selected high-quality spectrum bands, a variety of robust statistical methods are used to estimate the amplitude of each target frequency point to reduce the impact of the remaining outliers on the final result.

[0167] Preferred, robust statistical methods include:

[0168] The first approach: Define a set of high-quality spectrum bands. Includes Each effective spectrum segment is located at the target frequency. The amplitude at that point is denoted as ;

[0169] Each effective spectrum segment at the target frequency amplitude at The expression is as follows:

[0170] ;

[0171] in, This represents the modulus of a complex number.

[0172] Calculate the high-quality spectrum band at the target frequency amplitude at simple average And based on simple average Spectrum estimation is performed on each effective spectrum segment within the high-quality spectrum band.

[0173] Preferred, simple average The expression is as follows:

[0174] ;

[0175] The second approach: Define a set of high-quality spectrum bands. Includes Each effective spectrum segment is located at the target frequency. The amplitude at that point is denoted as ;

[0176] Calculate the median of each effective spectrum segment. And based on the median Spectrum estimation is performed on each effective spectrum segment within the high-quality spectrum band.

[0177] Preferred, median The expression is as follows:

[0178] ;

[0179] Indicated as filtered The set of amplitudes of each high-quality, effective spectral segment at the target frequency.

[0180] Preferred, set Median of each effective spectrum segment Specifically refers to Amplitude of each effective spectral segment The values ​​that are in the middle position after being sorted by size, where if If the number is even, the average of the two middle values ​​is taken.

[0181] Further preferred, due to the use of the median For sets When estimating the effective spectrum segments, the median is insensitive to outliers; however, when more than half of the data is not severely affected by noise, the median estimation method can provide reliable frequency amplitude estimates. The median has the highest robustness (theoretically resistant to nearly 50% of outliers), and is therefore often regarded as a reliable estimator and benchmark under strong noise conditions.

[0182] The third method involves discarding a certain proportion of the largest and smallest samples from the high-quality spectrum segments, and then averaging the remaining high-quality spectrum segments to obtain the truncated mean. ; and based on the truncated mean Spectrum estimation is performed on each effective spectrum segment within the high-quality spectrum band.

[0183] Preferably, the truncated mean is obtained. The specific process is as follows:

[0184] Set truncation ratio (like );

[0185] Cut-off mean The expression is as follows:

[0186] ;

[0187] in, Indicates the sorted order of the first... Each amplitude value.

[0188] By discarding extreme values, the truncated mean It takes into account the simple average value to a certain extent. and median Advantages: It utilizes most of the data to improve statistical efficiency while significantly reducing the impact of outliers on the results.

[0189] The fourth method: Obtain the mean IQR of the remaining effective spectrum segments using the mean estimation method based on the interquartile range (IQR). And based on the IQR mean Spectrum estimation is performed on each effective spectrum segment within the high-quality spectrum band.

[0190] Preferably, the mean IQR The expression is as follows:

[0191] ;

[0192] ;

[0193] ;

[0194] .

[0195] This invention uses statistical rules to automatically remove highly abnormal outliers, so that the calculated average value is not affected by extreme noise.

[0196] Fifthly, to reduce the influence of outliers in the averaging calculation, this invention applies the Huber robust estimation concept, assigning different weights to each effective spectrum segment before averaging to obtain the Huber weighted mean. Based on Huber weighted average By removing highly anomalous outliers from each effective spectral segment, the high efficiency of the mean method can be maintained while significantly reducing the bias of the results caused by anomalous noise.

[0197] Preferred, Huber weighted average The expression is as follows:

[0198] ;

[0199] Among them, weight Determined by the Huber weight function and calculated based on standardized residuals.

[0200] S6. Perform bias and consistency analysis on the spectrum estimation results of various robust statistical methods, and output the final spectrum;

[0201] To evaluate the reliability of the spectral amplitudes obtained by different estimation methods, a bias and consistency analysis was performed on the results of each method. The median of the estimated results was used. As a reference benchmark, the deviations and relative errors of other methods relative to this benchmark are compared.

[0202] make Represented by various robust statistical methods excluding the median Any estimation method other than (Representing the mean, truncated mean, IQR mean, or Huber mean, etc.), the frequency amplitude estimate of this estimation method is calculated, and the absolute deviation, relative error, and average relative error are defined as follows:

[0203] absolute deviation ;

[0204] relative error ;

[0205] Mean relative error ;in, The target frequency number.

[0206] The absolute deviation of each method at each target frequency point was calculated. and relative error The distribution can be used to determine its consistency and robustness with the median benchmark. For example, if a method's relative error is significant across all frequencies... If the relative error of a method remains within a very small range (e.g., within a few percentage points), it indicates that the estimate is very close to the median reference value and is minimally affected by noise; conversely, if the relative error of a method remains within a very small range (e.g., within a few percentage points), it indicates that the estimate is very close to the median reference value and is minimally affected by noise. If the method reaches a high level at certain frequencies, it indicates that the method may be significantly affected by residual anomalies at those frequencies. Based on the deviation analysis results, the method that is most stable in estimating spectral amplitude can be selected, providing a basis for the selection of spectral parameters in practical applications.

[0207] To further verify the practical application effect of the method described in Embodiment 1 of the present invention, the method was applied to wide-area electromagnetic field measurement data in a certain area. This received data not only contained continuous power frequency interference from the power grid, but also irregular, wideband, strong noise caused by the operation of power generation facilities.

[0208] Figures 3(a) and 3(b) illustrate a comparison of the measured original time-domain data and the data spectrum after PCA segmented denoising using this invention. In the upper part of the original spectrum, the overall background noise level is extremely high due to strong aperiodic electromagnetic interference, and the weak geological response signal is almost completely submerged by environmental noise, making accurate extraction difficult. However, the lower part of the PCA-denoised spectrum shows that after PCA segmented feature extraction and reconstruction using this invention, the aperiodic broadband background noise is significantly suppressed, and the signals at each target emission point (corresponding to the red dotted line in the figure) become prominent. The method of this invention has better denoising and signal recovery effects in the low-frequency band, the peak values ​​of weak low-frequency signals become clearer and more stable, and the overall signal-to-noise ratio is significantly improved.

[0209] Furthermore, combined Figure 4 As shown, a bar chart comparing the errors of amplitude estimation at seven different target frequencies using different statistical methods for high-quality spectrum bands is presented. The horizontal axis represents the frequency index, and the vertical axis represents the amplitude estimation error (mV). It is clear from the figure that the traditional simple mean estimation method (corresponding to the blue bars in the figure) is highly susceptible to interference from residual extreme outliers, exhibiting significant error fluctuations of up to 2-3.5 mV at multiple frequency points (such as frequencies 2, 5, and 7). In contrast, the robust statistical methods introduced in this embodiment, such as median estimation (corresponding to the orange bars in the figure) and weighted mean estimation (corresponding to the green bars in the figure), effectively suppress and control the estimation error to below 1 mV at all tested frequency points.

[0210] This fully demonstrates that the processing scheme of "PCA segmented noise reduction combined with multiple robust statistical estimations" proposed in this invention can effectively resist the bias of residual outliers on the true amplitude under strong noise background, significantly improve the accuracy, high consistency and high reliability of spectrum estimation of wide-area electromagnetic bathymetry data, and provide high-quality data support for subsequent construction and inversion of deep three-dimensional geological models.

[0211] Comparing the proposed method with traditional full-segment FFT spectrum processing methods, it can be seen that the proposed method has significant advantages in noise reduction robustness and estimation reliability. The proposed method, through PCA segmented reconstruction, effectively extracts recurring stable components in the received signal and filters out irrelevant random noise in each segment, significantly improving the signal-to-noise ratio. In contrast, traditional full-segment FFT processing directly performs spectrum transformation on the full-duration signal. If strong interference noise exists in the time-domain record, its energy will be evenly distributed across all frequencies of the spectrum, causing the amplitude at certain frequency points to be dominated by noise, resulting in unstable results. The proposed method, by segmenting and filtering the time-domain data, performs spectral feature fusion only on data segments with high signal-to-noise ratios, greatly reducing the impact of abnormal noise on the results. Simultaneously, the introduction of multiple robust statistical estimates ensures that even if a small number of deviations exist in the retained data segments, they will not cause serious deviations in the final spectral amplitude. For example, when the noise contains occasional spike interference, the spectral amplitude obtained by the traditional method may fluctuate significantly between different measurement times; however, with the proposed method, the segments corresponding to the spike interference will be identified and removed, and the final spectral amplitude will remain consistent across all observations. Furthermore, in cases where there is fixed-frequency narrowband interference, this invention can also reduce its impact in PCA noise reduction and subsequent statistical processes, because these narrowband noises are neither the dominant source of signal variation nor are they automatically weakened in robust estimation.

[0212] In summary, after the above processing steps, the spectrum extraction results of wide-area electromagnetic data are more stable and reliable, and the amplitude of weak signals is closer to the true level. Compared with traditional methods, this invention effectively improves the signal extraction accuracy and consistency of wide-area electromagnetic data under strong noise backgrounds, providing a high-quality input data foundation for subsequent underground electrical structure inversion. Overall, the method of this invention significantly improves the noise resistance and reliability of wide-area electromagnetic data processing without increasing additional hardware investment, providing more accurate and robust spectrum input for subsequent underground dielectric electrical parameter inversion. In terms of computational complexity, this invention mainly adds PCA matrix decomposition and multi-segment FFT operations, but these can be efficiently implemented through parallel computing and optimization algorithms, and the processing of conventional data scales can be completed within a reasonable time. Therefore, this method has good engineering practicality and is expected to be widely applied in wide-area electromagnetic data processing software systems.

[0213] Example 2:

[0214] As a further embodiment of the present invention, the present invention also provides an electronic device, comprising:

[0215] One or more processors;

[0216] Storage device for storing one or more programs;

[0217] When the one or more programs are executed by the one or more processors, the one or more processors implement the aforementioned method.

[0218] In practical use, users can interact with servers, which are also electronic devices, via a network to receive or send messages. Terminal devices are generally various electronic devices equipped with a display and used through a human-computer interface, including but not limited to smartphones, tablets, laptops, and desktop computers. Various specific application software can be installed on these terminal devices as needed, including but not limited to web browsers, instant messaging software, social media platforms, and shopping apps.

[0219] Furthermore, a server is a network service that provides various services, such as a backend server that provides corresponding calculation services for wide-area electromagnetic data transmitted from terminal devices, so as to realize the processing of wide-area electromagnetic data processing methods, calculate the final spectrum, and return it to the terminal devices.

[0220] Example 3:

[0221] As a further embodiment of the present invention, the present invention also provides a storage medium including one or more programs executable by one or more processors of an electronic device, the one or more programs including instructions for performing the wide-area electromagnetic data processing method as described above.

[0222] The above description is merely a preferred embodiment of the present invention and is not intended to limit the invention. Various modifications and variations can be made to the present invention by those skilled in the art. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the present invention should be included within the scope of protection of the present invention.

Claims

1. A wide-area electromagnetic data processing method, characterized in that, Includes the following steps: Step 1: Acquire the raw time-domain data obtained by the wide-area electromagnetic method. ; Step 2: Process the raw time-domain data The data is segmented, and PCA is used to analyze the original time-domain data after segmentation. Denoising is performed to obtain segmented denoised time-domain data; Step 3: Perform piecewise discrete Fourier transform analysis on the piecewise denoised time-domain data to obtain the spectrum of each segment; Step 4: Calculate the signal-to-noise ratio (SNR) of the target frequency in each spectrum segment, sort the spectrum segments according to their SNR, and select high-quality spectrum segments. Step 5: Use multiple robust statistical methods to estimate the amplitude of each target frequency point in the selected high-quality spectrum bands, and obtain the spectrum estimation results of multiple robust statistical methods; Step 6: Perform bias and consistency analysis on the spectrum estimation results of various robust statistical methods, and output the final spectrum; The specific process of step two is as follows: Original time domain data Divided by length Segment time-domain data, assuming each segment of time-domain data contains There are consecutive sampling points, and there is a certain distance between two adjacent time-domain data segments. The proportions of the sample sizes overlap; All time-domain data are used to construct corresponding column vectors and row vectors, and all column vectors and row vectors are stacked into a data matrix. ; For data matrix Perform principal component analysis and based on the data matrix The strategy for selecting the number of principal components to retain in the time-domain data based on the principal component analysis results determines the amount of time-domain data to be retained. One principal component; Construct the truncated singular value matrix and utilize Reconstructing the data matrix from principal components ; Data matrix By restoring the data to its time series form, the reconstructed time-domain data can be obtained. For reconstructing time-domain data For the segmented reconstructed time-domain data with overlapping non-Chinese and African components, each segment of the reconstructed time-domain data is directly concatenated sequentially to obtain the segmented denoised time-domain data; for the reconstructed time-domain data... The time-domain data is reconstructed by segmenting overlapping segments. The average of the signals in the overlapping segments is taken and smoothly connected to obtain the segmented denoised time-domain data.

2. The wide-area electromagnetic data processing method according to claim 1, characterized in that, Using an index matrix to analyze the original time-domain data Perform segmentation processing to obtain Time-domain data; The index matrix is ; ； in, The starting position, and , segment number and ; Step size; To round down; The total number of segments, and ; This is the total length of the window; This represents the total number of sampling points for the original time-domain data.

3. The wide-area electromagnetic data processing method according to claim 1, characterized in that, For data matrix The specific method for performing principal component analysis is as follows: For data matrix Centralized processing is performed to obtain a centralized data matrix. ; For centralized data matrix Perform principal component analysis to decompose the left singular vector matrix. Singular value matrix and right singular vector matrix .

4. The wide-area electromagnetic data processing method according to claim 1, characterized in that, Based on the data matrix The specific method for selecting the number of principal components to retain in the time-domain data based on the principal component analysis results is as follows: By analyzing the singular value spectrum, the boundary between signal and noise is determined, thus identifying whether to retain the time-domain data using either a directly specified method or an automatically determined method based on the cumulative variance contribution rate. One principal component.

5. The wide-area electromagnetic data processing method according to any one of claims 2-4, characterized in that, The specific process of step three is as follows: by The length is used to further divide the segmented, denoised time-domain data into segments. Segmented data in the time domain; The discrete Fourier transform of each segment of time-domain data is calculated by multiplying it by the Hanning window function to obtain the segmented spectrum.

6. The wide-area electromagnetic data processing method according to claim 5, characterized in that, Obtain a single-segment spectrum The specific process is as follows: Set the Hanning window function for single-segment time-domain segmented data The windowed time-domain data is obtained. ; Based on windowed time-domain data Perform Discrete Fourier Transform calculations to obtain the Discrete Fourier Transform results. ; Based on the results of Discrete Fourier Transform Calculate the positive frequency component of the one-sided amplitude spectrum. ; and the recovery coefficient based on the Hanning window For the positive frequency portion of the single-sided amplitude spectrum To perform recovery and obtain the extent of recovery. That is, the true physical amplitude of the Y-axis of a single-segment time-domain spectrum; Calculate the frequency axis That is, to calculate the true physical frequency of the X-axis of a single-segment time-domain spectrum; and to set the discrete index points This is mapped to a specific transmission frequency, resulting in a frequency axis. The actual spectrum; Based on the magnitude of recovery and frequency axis Calculate the first Segment data on the frequency axis Complete spectrum data at [location] This yields a single-segment spectrum.

7. The wide-area electromagnetic data processing method according to claim 6, characterized in that, The specific process of step four is as follows: Calculate the spectrum of a single segment The ratio of signal strength to background noise at the main transmission frequency point is recorded as the signal-to-noise ratio (SNR) of that frequency band. ; Signal-to-noise ratio of each spectrum Sort the spectrum from high to low, and set the percentage threshold for signal-to-noise ratio (SNR) filtering to P%. Retain the SNR values ​​of each spectrum in the top P% as high-quality spectrum segments.

8. The wide-area electromagnetic data processing method according to claim 7, characterized in that, The robust statistical methods in step five include at least three of the following five methods, and each of the at least three methods is used to estimate the spectrum of the selected high-quality spectrum bands, resulting in at least three different spectrum estimation results: The first approach: Define a set of high-quality spectrum bands. Includes Each effective spectrum segment is located at the target frequency. The amplitude at that point is denoted as ; Calculate the high-quality spectrum band at the target frequency amplitude at simple average And based on simple average Spectrum estimation is performed on each effective spectrum segment within the high-quality spectrum band. The second approach: Define a set of high-quality spectrum bands. Includes Each effective spectrum segment is located at the target frequency. The amplitude at that point is denoted as ; Calculate the median of each effective spectrum segment. And based on the median Spectrum estimation is performed on each effective spectrum segment within the high-quality spectrum band. The third method involves discarding a certain proportion of the largest and smallest samples from the high-quality spectrum segments, and then averaging the remaining high-quality spectrum segments to obtain the truncated mean. ; and based on the truncated mean Spectrum estimation is performed on each effective spectrum segment within the high-quality spectrum band. The fourth method: Obtain the IQR mean of the remaining effective spectrum segments using the mean estimation method based on the interquartile range. And based on the IQR mean Spectrum estimation is performed on each effective spectrum segment within the high-quality spectrum band. The fifth method applies the Huber robust estimation concept, assigning different weights to each effective spectrum segment and then averaging the results to obtain the Huber weighted mean. And based on Huber's weighted average Spectrum estimation is performed on each effective spectrum segment within the high-quality spectrum band.

9. The wide-area electromagnetic data processing method according to claim 8, characterized in that, The specific process of step six is ​​as follows: Median of the result estimated by median As a reference benchmark, the deviations and relative errors of the first, third to fifth estimation methods relative to the reference benchmark are compared; set up Represented as any estimation method, the absolute deviation of that estimation method at each target frequency point is calculated. and relative error Distribution, determine whether this estimation method is similar to the median. To assess the consistency and robustness of the estimation method, the bias analysis results were obtained. Based on the deviation analysis results of all estimation methods, the estimation method that is most stable for estimating the spectral amplitude is selected, and its estimation result is output as the final spectrum.

10. An electronic device, characterized in that, It includes a memory, one or more processes, and one or more programs stored in the memory, said one or more programs including instructions for performing the wide-area electromagnetic data processing method as described in any one of claims 1-9.

11. A storage medium, characterized in that, It includes one or more programs executable by one or more processors of an electronic device, the one or more programs including instructions for performing the wide-area electromagnetic data processing method as described in any one of claims 1-9.

Images