A method for processing gravity and magnetic data based on wavelet decomposition and multifractal

Abstract

A gravity and magnetic data processing method based on wavelet decomposition and multifractal analysis is presented, belonging to the field of geophysical exploration technology. This method aims to address the problems of traditional filtering techniques, such as over-smoothing weak anomalies and difficulty in distinguishing signal and noise singularities. The steps include: acquiring raw gravity and magnetic data and performing normalization; selecting a wavelet basis and decomposition level, performing L-level wavelet decomposition on the normalized data to obtain approximation coefficients and detail coefficients; calculating the adjusted Holder exponent based on the detail coefficients to quantify signal singularities; designing an adaptive threshold by combining the Holder exponent and the number of sampling points, and filtering the detail coefficients through triple mask fusion; reconstructing the filtered coefficients using inverse wavelet transform, followed by adaptive moving average smoothing and denormalization to obtain the final filtering result. This invention integrates the multi-scale decomposition capability of wavelet transform with the singularity identification advantages of multifractals, achieving precise noise suppression and complete preservation of useful signals, improving the reliability of gravity and magnetic data, and providing technical support for deep geological exploration, mineral exploration, and other tasks.

Description

Technical Field

[0001] This invention belongs to the field of geophysical exploration technology, specifically relating to a gravity and magnetic data processing method based on a combination of wavelet decomposition and multifractal methods. Background Technology

[0002] In airborne geophysical exploration, core tasks such as deep geological exploration, strategic mineral exploration, and oil and gas reservoir evaluation place extremely high demands on the accuracy and reliability of exploration data. Superconducting gravity gradiometers, with their 1E-level sensitivity, can capture subtle gravity gradient changes caused by density differences in geological bodies. They are a core support for in-depth data processing in this field (anomaly separation, noise suppression, etc.), laying the foundation for high-precision geological interpretation. However, raw airborne gravity gradient data is susceptible to interference from multiple sources of noise: high-frequency random noise introduced by flight platform attitude jitter and airflow disturbances, as well as instrument electronic noise and external electromagnetic interference. This noise, mixed with weak and locally singular useful gravity gradient anomaly signals, reduces the data signal-to-noise ratio, interferes with subsequent tensor inversion and geological body boundary positioning, and further restricts the identification of weak gradient anomalies caused by small deep ore bodies and concealed fault zones, becoming a bottleneck for the field's refined and in-depth development. Traditional filtering methods have limitations: although Gaussian filtering can suppress noise, it is prone to over-smoothing useful weak anomalies and losing geological details; traditional wavelet denoising relies on coefficient amplitude threshold screening, which makes it difficult to accurately distinguish noise from signal singularities, and its robustness is insufficient at low signal-to-noise ratios. Summary of the Invention

[0003] To address the aforementioned problems, the purpose of this invention is to provide a gravity and magnetic data processing method based on a combination of wavelet decomposition and multifractal methods, thereby solving the problems.

[0004] To achieve the above objectives, the technical solution adopted by the present invention includes: A gravity and magnetic data processing method based on wavelet decomposition and multifractal combination includes the following steps: S1, obtain the length as N Raw gravity and magnetic data X and the original gravity and magnetic data X Normalization is performed to obtain normalized data. Xnorm ; N Raw gravity and magnetic data X The number of sampling points.

[0005] S2. Selecting wavelet basis and decomposition level L Normalized data obtained from S1 Xnorm conduct L Wavelet decomposition of level 1 yields the 1st wavelet. L Approximation coefficients of the layer A L and from the 1st floor to the L Layer detail factorD j ,in j = 1, 2, ..., L ; S3. Using the detail coefficients Dj obtained from S2, calculate the global Holder exponent for each layer. and the local window energy of each layer's detail coefficient at position k And through the global Holder index at each layer and the local window energy of each layer's detail coefficient at position k The adjusted Holder index for each layer was calculated. The adjusted Holder index for each layer Including the adjusted Holder exponent at position k in layer j ; S4, the adjusted Holder index obtained from S3 and number of sampling points N For detail coefficients D j Perform filtering to obtain the filtered detail coefficients. D after j and approximation coefficients A L As the approximation coefficients after filtering A after L ; S5, the filtered detail coefficients obtained from S4 D after j and the approximate coefficients after filtering A after L Perform inverse wavelet transform to reconstruct the time-domain signal. Xnorm after Then, for the time-domain signal... Xnorm after Adaptive moving average smoothing is performed, and the normalized data obtained from S1 is used. Xnorm Perform denormalization to obtain the final filtering result. X result .

[0006] Preferably, wavelet bases and decomposition levels are selected in S2. L Specifically, based on the signal length N and the filter length of the selected wavelet basis dec_len Determine the maximum decomposition level ,in, This represents the floor function for rounding down; and takes...L= min (L input ,L max ) This represents the actual decomposition level, where L input This is the preset level.

[0007] Preferred, preset level L input Take 3.

[0008] Preferably, the adjusted Holder index at position k in the j-th layer of S3. As shown in the following formula: ;in, Let be the average energy of the detail coefficients at the j-th layer.

[0009] Preferably, S4 specifically includes: S4.1, For all detail coefficients D j The noise standard deviation σ is estimated using the median absolute deviation method. S4.2. Based on the noise standard deviation and the number of sampling points N obtained in S4.1, calculate the basic threshold. T base ; S4.3, Based on the fully adjusted Holder index The base threshold obtained from S4.2 T base Dynamic adjustments are made to obtain an adaptive threshold. T The dynamic adjustment is as follows: ; S4.4 Mask Generation and Fusion: (a) For detail coefficients D j Each coefficient in D j（k） Determine its absolute value | D j（k） | Is it greater than or equal to the adaptive threshold obtained in S4.3? T If so, the generated first mask holder_mask at position k indicates that the coefficient satisfies the continuity condition; where k represents the coefficient at any detail coefficient. D j (a) Position index in the middle; (b) For detail coefficients D j Each coefficient in D j（k） Determine its absolute value |D j（k）| Is it greater than or equal to the reinforcement threshold γ? σ; if so, the generated second mask amplitude_mask indicates at position k that the coefficient satisfies the amplitude condition; where γ is a preset coefficient greater than 1; (c) for detail coefficients D j Each coefficient in D j（k） Check if there are other coefficients within its preset neighborhood that are indicated by the first mask holder_mask or the second mask amplitude_mask as satisfying the condition; if so, the generated third mask continuity_mask indicates at position k that the coefficient satisfies the continuity condition; (d) Generate a combined mask combined_mask from the first mask holder_mask, the second mask amplitude_mask and the third mask continuity_mask; for the coefficient at position k, the combined mask combined_mask indicates that the coefficient is retained only if the first mask holder_mask indicates that it satisfies the Holder exponent condition, and the second mask amplitude_mask indicates that it satisfies the amplitude condition or the third mask continuity_mask indicates that it satisfies the continuity condition; S4.5. Preserve detail coefficients according to the instructions of the combined_mask. D j The coefficients that were instructed to be retained are filtered out, and the coefficients that were not instructed to be retained are removed to obtain the filtered detail coefficients. D after j .

[0010] Preferably, in S5, the adaptive moving average smoothing specifically involves: calculating the time-domain signal. Xnorm after The local variance is normalized to the [0,1] interval to obtain the normalized variance value. According to the rule of dynamically selecting the moving average window size, the moving average window size is dynamically selected for smoothing.

[0011] Preferably, the rule for dynamically selecting the size of the moving average window is as follows: when the normalized variance is greater than 0.6, a 3-point window is used; when the normalized variance is between 0.3 and 0.6, a 5-point window is used; and when the normalized variance is less than 0.3, a 7-point window is used.

[0012] Preferred raw gravity and magnetic data X This includes airborne gravity gradient data, ground gravity anomaly data, or magnetic anomaly data.

[0013] A computer-readable storage medium storing a computer program, which, when executed by a processor, implements the gravity and magnetic data processing method based on wavelet decomposition and multifractal combination disclosed in this application.

[0014] A computer program product includes a computer program / instructions that, when executed by a processor, implement the gravity and magnetic data processing method based on wavelet decomposition and multifractal combination disclosed in this application.

[0015] Compared with the prior art, the advantages of the present invention are: (1) Accurately distinguish singular components: By using multifractal theory to quantitatively describe the distribution of signal singularity, the noise with "high singularity and low amplitude" can be effectively separated from the useful geological anomaly signal with "low singularity and high amplitude", thus overcoming the limitation of traditional methods that are difficult to distinguish between the two. (2) Both noise suppression and detail preservation: The wavelet transform’s “time-frequency” localization capability enables multi-scale signal decomposition. Combined with adaptive threshold and triple mask fusion filtering, it can effectively suppress multi-source noise such as flight attitude jitter and electromagnetic interference, while fully preserving the subtle anomaly details corresponding to weak geological bodies such as deep small ore bodies and hidden fault zones. (3) Improve data reliability and interpretation accuracy: By dynamically adjusting the Holder exponent threshold and adaptive moving average smoothing, the signal-to-noise ratio of gravity and magnetic data is significantly improved, providing high-precision data support for subsequent tensor inversion, geological body boundary positioning and other work, and helping airborne geophysical exploration to develop towards refinement and depth; (4) Strong robustness and applicability: It is compatible with various gravity and magnetic data such as airborne gravity gradient data, ground gravity anomaly data, and magnetic anomaly data. Furthermore, the wavelet basis selection is flexible and the decomposition level adapts to the signal length. It maintains a stable filtering effect even in low signal-to-noise ratio scenarios, which is superior to the over-smoothing of Gaussian filtering and the insufficient robustness of traditional wavelet denoising. Attached Figure Description

[0016] The accompanying drawings are provided to further illustrate the invention and form part of the specification. They are used together with the following detailed description to explain the invention, but do not constitute a limitation thereof. In the drawings: Figure 1 It is a comparison between the original signal and the noisy signal; Figure 2 This is a graph of wavelet decomposition coefficients; Figure 3 This is a comparison chart of the reconstructed signal and the smoothed signal; Figure 4 This is a comparison chart of the final filtered result and the original signal. Detailed Implementation

[0017] The invention is not limited to the specific embodiments described below. All equivalent modifications made based on the technical solutions of this application fall within the protection scope of this invention. Unless otherwise specified, all components and devices in this invention utilize components and devices known in the prior art.

[0018] This application discloses a gravity and magnetic data processing method based on a combination of wavelet decomposition and multifractal methods, comprising the following steps: S1, obtain the length as N Raw gravity and magnetic data X and the original gravity and magnetic data X Normalization is performed to obtain normalized data. Xnorm ; N Raw gravity and magnetic data X The number of sampling points.

[0019] Original gravity and magnetic data in this embodiment X This includes airborne gravity gradient data, ground gravity anomaly data, or magnetic anomaly data.

[0020] S2. Selecting wavelet basis and decomposition level L Normalized data obtained from S1 Xnorm conduct L Wavelet decomposition of level 1 yields the 1st wavelet. L Approximation coefficients of the layer A L and from the 1st floor to the L Layer detail factor D j ,in j = 1, 2, ..., L ; Selection of wavelet basis and decomposition level L Specifically: Based on signal length N and the filter length of the selected wavelet basis dec_len Determine the maximum decomposition level ,in, The floor function represents rounding down; and take L= min (L input ,L max ) This represents the actual decomposition level, where L input As a preset level, this embodiment uses level 3.

[0021] Experiments have verified that the wavelet basis has a relatively small impact on this embodiment. The selection of the wavelet basis does not require analysis and selection; commonly used ones are sufficient. In this embodiment, the filter length dec_len of the wavelet basis is determined by the principle of discrete wavelet decomposition itself, and the signal length becomes 1 / 2 of the original after one decomposition.

[0022] S3. Using the detail coefficients Dj obtained from S2, calculate the global Holder exponent for each layer. and the local window energy of each layer's detail coefficient at position k And through the global Holder index at each layer and the local window energy of each layer's detail coefficient at position k The adjusted Holder index for each layer was calculated. ; Adjusted Holder index for each layer Including the adjusted Holder exponent at position k in layer j As shown in the following formula:

[0023] in, Let be the average energy of the detail coefficients at the j-th layer.

[0024] S4, the adjusted Holder index obtained from S3 and number of sampling points N For detail coefficients D j Perform filtering to obtain the filtered detail coefficients. D after j and approximation coefficients A L As the approximation coefficients after filtering A after L ; S4 specifically includes: S4.1, For all detail coefficients D j The noise standard deviation σ is estimated using the median absolute deviation method. S4.2. Based on the noise standard deviation and the number of sampling points N obtained in S4.1, calculate the basic threshold. T base ; S4.3, Based on the fully adjusted Holder index The base threshold obtained from S4.2 T base Dynamic adjustments are made to obtain an adaptive threshold. T The dynamic adjustment is as follows: ; S4.4 Mask Generation and Fusion: (a) Regarding detail factor D j Each coefficient in D j（k） Determine its absolute value |D j（k） | Is it greater than or equal to the adaptive threshold obtained in S4.3? T If so, the generated first mask holder_mask indicates at position k that the coefficient satisfies the continuity condition; Where k represents the coefficient in any detail coefficient D j Position index in; (b) Regarding detail factor D j Each coefficient in D j（k） Determine its absolute value |D j（k） | Is it greater than or equal to the reinforcement threshold γ? σ; if so, the generated second mask amplitude_mask indicates at position k that the coefficient satisfies the amplitude condition; Where γ is a preset coefficient greater than 1, and in this embodiment, 3 is preferred; (c) Regarding detail factor D j Each coefficient in D j（k） Check if there are other coefficients in its preset neighborhood that are indicated by the first mask holder_mask or the second mask amplitude_mask to meet the condition; if so, the generated third mask continuity_mask indicates at position k that the coefficient meets the continuity condition. (d) Generate a combined mask, combined_mask, from the first mask holder_mask, the second mask amplitude_mask, and the third mask continuity_mask; For the coefficient at position k, the combined mask indicates that the coefficient is retained only if the first mask holder_mask indicates that it satisfies the Holder exponent condition, the second mask amplitude_mask indicates that it satisfies the amplitude condition, or the third mask continuity_mask indicates that it satisfies the continuity condition. S4.5. Preserve detail coefficients according to the instructions of the combined_mask. D j The coefficients that were instructed to be retained are filtered out, and the coefficients that were not instructed to be retained are removed to obtain the filtered detail coefficients. D after j .

[0025] S5, the filtered detail coefficients obtained from S4 D after j and the approximate coefficients after filtering A after L Perform inverse wavelet transform to reconstruct the time-domain signal. Xnorm after Then, for the time-domain signal... Xnorm after Adaptive moving average smoothing is performed, and the normalized data obtained from S1 is used. Xnorm Perform denormalization to obtain the final filtering result. X result .

[0026] The adaptive moving average smoothing in this embodiment is specifically as follows: Calculate time-domain signals Xnorm after The local variance is normalized to the [0,1] interval to obtain the normalized variance value. According to the rule of dynamically selecting the moving average window size, the moving average window size is dynamically selected for smoothing.

[0027] The rule for dynamically selecting the size of the moving average window in this embodiment is as follows: When the normalized variance is greater than 0.6, use a 3-point window; Use a 5-point window when the normalized variance is between 0.3 and 0.6; When the normalized variance is less than 0.3, a 7-point window is used.

[0028] This application also discloses a computer-readable storage medium storing a computer program, which, when executed by a processor, implements the gravity and magnetic data processing method based on wavelet decomposition and multifractal combination disclosed in this application.

[0029] This application also discloses a computer program product, including a computer program / instruction, which, when executed by a processor, implements the gravity and magnetic data processing method based on wavelet decomposition and multifractal combination disclosed in this application.

[0030] Example This embodiment uses a one-dimensional analog signal as the processing object, and the specific steps are as follows: Raw data preparation and preprocessing: The input consists of raw data composed of multiple superimposed sine and cosine signals with different amplitudes and phases (e.g., ...). Figure 1 (As shown), and Gaussian noise with a maximum value of 100% is added to simulate the disturbed gravity and magnetic data scenario; the original data is then normalized to obtain normalized data. Xnorm Meanwhile, the normalization parameters (Xmax and Xmin) are saved for subsequent denormalization to recover the signal.

[0031] Wavelet decomposition parameter selection and decomposition execution: The db4 wavelet is selected as the wavelet basis, with a preset decomposition level Linput=3; based on the signal length N and the filter length dec_len of the db4 wavelet basis, the maximum decomposition level Lmax is calculated using the formula, and the actual decomposition level L=min(3,Lmax)=3 is finally determined; a 3-level wavelet decomposition is performed on the normalized data Xnorm to obtain the 3rd level approximation coefficients AL and the 1st-3rd level detail coefficients Dj (j=1,2,3), as shown in the figure. Figure 2 As shown.

[0032] Holder index calculation: Based on the detail coefficients Dj of each layer, the global Holder index is first calculated through cross-scale analysis. Then calculate the local window energy at each position k. and the average energy of the detail factor of that layer Using formulas The adjusted Holder index at position k in the j-th layer is obtained. The Holder index after adjustment at each position is obtained. .

[0033] Adaptive filtering: The median absolute deviation method is used to estimate the noise standard deviation σ for all detail coefficients Dj; the base threshold is calculated by combining the number of sampling points N. T base Then through the formula Based on the fully adjusted Holder index The adaptive threshold T is obtained through adjustment; three types of masks (holder_mask, amplitude_mask, continuity_mask) are generated and fused into a combined mask (combined_mask). Only coefficients satisfying the "Holder exponential condition + (amplitude condition or continuity condition)" are retained to obtain the filtered detail coefficients. D after j Meanwhile, the approximation coefficient AL is retained as A after L .

[0034] Signal reconstruction and smoothing optimization: A after L and D after j Perform inverse wavelet transform to reconstruct the time-domain signal. Xnorm after (like Figure 3 The reconstructed signal is shown in the figure); calculation Xnorm after The local variance was normalized to [0,1], and then adaptively smoothed using a moving average according to the rule of "variance > 0.6 using a 3-point window, 0.3~0.6 using a 5-point window, and < 0.3 using a 7-point window" (e.g., ...). Figure 3 (The smoothed signal is shown in the image). The normalization parameters saved in step 1 are used for denormalization to obtain the final filtered result. X result Its comparison with the original signal is as follows: Figure 4 As shown.

[0035] The preferred embodiments of this disclosure have been described in detail above with reference to the accompanying drawings. However, this disclosure is not limited to the specific details of the above embodiments. Within the scope of the technical concept of this disclosure, various simple modifications can be made to the technical solutions of this disclosure, and these simple modifications all fall within the protection scope of this disclosure.

[0036] It should also be noted that the various specific technical features described in the above specific embodiments can be combined in any suitable manner without contradiction. In order to avoid unnecessary repetition, this disclosure will not describe the various possible combinations separately.

[0037] Furthermore, the various implementation methods disclosed in this solution can be combined in any way, as long as they do not violate the spirit of this disclosure, they should also be regarded as the content invented by this disclosure.

Claims

1. A gravity and magnetic data processing method based on a combination of wavelet decomposition and multifractal methods, characterized in that, Includes the following steps: S1, obtain the length as N Raw gravity and magnetic data X and the original gravity and magnetic data X Normalization is performed to obtain normalized data. Xnorm ; The N Raw gravity and magnetic data X The number of sampling points. S2. Selecting wavelet basis and decomposition level L Normalized data obtained from S1 Xnorm conduct L Wavelet decomposition of level 1 yields the 1st wavelet. L Approximation coefficients of the layer A L and from the 1st floor to the L Layer detail factor D j ,in j = 1, 2, ..., L ; S3. Using the detail coefficients Dj obtained from S2, calculate the global Holder exponent for each layer. and the local window energy of each layer's detail coefficient at position k And through the global Holder index of each layer and the local window energy of each layer's detail coefficient at position k The adjusted Holder index for each layer was calculated. ; The adjusted Holder index corresponding to each layer Including the adjusted Holder exponent at position k in layer j ; S4, the adjusted Holder index obtained from S3 and the number of sampling points N For the detail coefficients D j Perform filtering to obtain the filtered detail coefficients. D after j and the approximation coefficients A L As the approximation coefficients after filtering A after L ; S5, the filtered detail coefficients obtained from S4 D after j and the approximate coefficients after filtering A after L Perform inverse wavelet transform to reconstruct the time-domain signal. Xnorm after Then, the time-domain signal... Xnorm after Adaptive moving average smoothing is performed, and the normalized data obtained from S1 is used. Xnorm Perform denormalization to obtain the final filtering result. X result .

2. The gravity and magnetic data processing method based on wavelet decomposition and multifractal combination as described in claim 1, characterized in that, The selection of wavelet bases and decomposition levels described in S2 L Specifically: Based on signal length N and the filter length of the selected wavelet basis dec_len Determine the maximum decomposition level ,in, The floor function represents rounding down; and take L= min (L input ,L max ) This represents the actual decomposition level, where L input This is the preset level.

3. The gravity and magnetic data processing method based on wavelet decomposition and multifractal combination as described in claim 2, characterized in that, The preset level L input Take 3.

4. The gravity and magnetic data processing method based on wavelet decomposition and multifractal combination as described in claim 1, characterized in that, The adjusted Holder index at position k in layer j as described in S3 As shown in the following formula: in, Let be the average energy of the detail coefficients at the j-th layer.

5. The gravity and magnetic data processing method based on wavelet decomposition and multifractal combination as described in claim 1, characterized in that, S4 specifically includes: S4.1, For all detail coefficients D j The noise standard deviation σ is estimated using the median absolute deviation method. S4.

2. Based on the noise standard deviation obtained in S4.1 and the number of sampling points N, calculate the basic threshold. T base ; S4.3, Based on the fully adjusted Holder index The base threshold obtained from S4.2 T base Dynamic adjustments are made to obtain an adaptive threshold. T The dynamic adjustment is as follows: ； S4.4 Mask Generation and Fusion: (a) Regarding the aforementioned detail factor D j Each coefficient in D j（k） Determine its absolute value |D j（k） | Is it greater than or equal to the adaptive threshold obtained in S4.3? T If so, the generated first mask holder_mask indicates at position k that the coefficient satisfies the continuity condition; Where k represents the coefficient in any detail coefficient D j Position index in; (b) Regarding the aforementioned detail factor D j Each coefficient in D j（k） Determine its absolute value |D j（k） | Is it greater than or equal to the reinforcement threshold γ? σ; if so, the generated second mask amplitude_mask indicates at position k that the coefficient satisfies the amplitude condition; Where γ is a preset coefficient greater than 1; (c) Regarding the aforementioned detail factor D j Each coefficient in D j（k） Check if there are other coefficients in its preset neighborhood that are indicated by the first mask holder_mask or the second mask amplitude_mask to meet the condition; if so, the generated third mask continuity_mask indicates at position k that the coefficient meets the continuity condition; (d) Generate a combined mask, combined_mask, from the first mask holder_mask, the second mask amplitude_mask, and the third mask continuity_mask; For the coefficient at position k, the combined mask indicates that the coefficient is retained only if the first mask holder_mask indicates that it satisfies the Holder exponent condition, the second mask amplitude_mask indicates that it satisfies the amplitude condition, or the third mask continuity_mask indicates that it satisfies the continuity condition. S4.

5. Based on the indication of the combined_mask, preserve the detail coefficients. D j The coefficients that were instructed to be retained are filtered out, and the coefficients that were not instructed to be retained are removed to obtain the filtered detail coefficients. D after j .

6. The gravity and magnetic data processing method based on wavelet decomposition and multifractal combination as described in claim 1, characterized in that, In S5, the adaptive moving average smoothing specifically refers to: Calculate the time-domain signal Xnorm after The local variance is normalized to the [0,1] interval to obtain the normalized variance value. According to the rule of dynamically selecting the moving average window size, the moving average window size is dynamically selected for smoothing.

7. The gravity and magnetic data processing method based on wavelet decomposition and multifractal combination as described in claim 6, characterized in that, The rule for dynamically selecting the size of the moving average window is as follows: When the normalized variance is greater than 0.6, use a 3-point window; Use a 5-point window when the normalized variance is between 0.3 and 0.6; When the normalized variance is less than 0.3, a 7-point window is used.

8. The gravity and magnetic data processing method based on wavelet decomposition and multifractal combination as described in any one of claims 1-7, characterized in that, The original gravity and magnetic data X This includes airborne gravity gradient data, ground gravity anomaly data, or magnetic anomaly data.

9. A computer-readable storage medium, characterized in that, The computer-readable storage medium stores a computer program that, when executed by a processor, implements the gravity and magnetic data processing method based on wavelet decomposition and multifractal combination as described in any one of claims 1-8.

10. A computer program product, characterized in that, Includes a computer program / instruction, which, when executed by a processor, implements the gravity and magnetic data processing method based on wavelet decomposition and multifractal combination as described in any one of claims 1-8.

Images