An equivalent source denoising method based on shape regularization

Abstract

The present application is suitable for the field of geophysical exploration data processing technology, and provides an equivalent source denoising method based on shaping regularization, which comprises the following steps: obtaining multi-component data; constructing an equivalent source forward model; balancing the contribution of each component through robust scale normalization; introducing zero trace relationship as a physical soft constraint into the inversion system; adopting shaping regularization to stably solve the process and adaptively determining the regularization parameter through HKB method; and finally obtaining the denoising result through forward calculation. The present application realizes integrated denoising of multi-component data, significantly improves the component consistency, physical rationality and processing efficiency of the denoising result, and is suitable for fine processing of airborne and ground gravity gradient measurement data.

Description

Technical Field

[0001] This invention belongs to the field of geophysical exploration data processing technology, and particularly relates to an equivalent source denoising method based on shaping regularization. Background Technology

[0002] Full Tensor Gradiometry (FTG) measurements can obtain information about the second derivative of the gravitational potential function, i.e., the components of the gravity gradient tensor. Compared with traditional gravity anomaly data, gravity gradient data is more sensitive to spatial variations, has higher lateral resolution, and stronger response capability to shallow anomalies. Therefore, it is widely used in fields such as airborne gravity gradient measurement, mineral resource exploration, structural identification, and geological mapping.

[0003] However, during the acquisition and processing of FTG data, it is susceptible to various factors such as platform attitude disturbances, navigation and positioning errors, instrument background noise, and environmental disturbances. This results in random and structured noise being superimposed on each component of the data, reducing the physical consistency and interpretability of the data. Therefore, effective denoising processing must be performed before data interpretation and imaging.

[0004] Currently, most common FTG denoising schemes in engineering applications adopt a "component-based independent processing" approach. This approach treats each gradient component as an independent data channel, denoising each component separately before combining them for output. Its main process includes data regularization, single-component filtering (such as moving average, Gaussian filtering, wavelet thresholding, etc.), component-level parameter adjustment, and result output. This type of method is simple to implement and can reduce the noise level of individual components to a certain extent.

[0005] However, since FTG multi-component data originates from the same gravitational potential field and has a strict physical coupling relationship, single-component independent denoising has the following inherent defects:

[0006] (1) Component consistency is difficult to guarantee: Each component of the FTG corresponds to the second derivative combination of the same potential field. Existing technologies apply independent filters and parameters to each component, which can easily change the phase relationship, boundary response and amplitude ratio between components. This leads to inconsistencies in the spatial location, anomalous boundary and other characteristics of the same anomalous body in different components, which reduces the reliability of the tensor data as a whole for interpretation.

[0007] (2) Difficulty in maintaining harmonic properties in the passive region: The gravitational field satisfies the Laplace equation in the passive observation region, and the diagonal component of its gradient tensor should satisfy the zero-trace relation. The single-component denoising process lacks explicit control over this physical constraint. Inconsistency in the filtering intensity will cause a systematic deviation from the zero-trace relation, forming spatially correlated "trace residuals" or non-harmonic components. These components may manifest as pseudo-anomalies or background drift in subsequent interpretations.

[0008] (3) The selection of regularization parameters is blind and costly: Existing technologies mainly rely on manual trial and error or L-curve method to determine regularization parameters. Manual trial and error is highly subjective and inefficient; L-curve method requires complete inversion calculation of multiple parameter values ​​to determine the inflection point, which is extremely costly and difficult to meet the needs of rapid imaging. Summary of the Invention

[0009] The purpose of this invention is to provide an equivalent source denoising method based on shaping regularization, which aims to solve the problems mentioned in the background art.

[0010] The present invention is implemented as follows: an equivalent source denoising method based on integer regularization includes the following steps:

[0011] Step 1: Obtain the full tensor gravity gradient data and express it as a vector;

[0012] Step 2: Construct an equivalent source forward model, discretize the equivalent source layer into multiple prism elements, and establish a linear forward modeling equation from element density to observation data;

[0013] Step 3: Perform robust scaling normalization on each component of the full tensor gravity gradient data, estimate the noise standard deviation based on the median absolute deviation of each component, and construct a diagonal weighted matrix to balance the contribution of each component.

[0014] Step 4: Utilize the zero-trace property of the gravity gradient tensor to introduce the zero-trace condition as a soft constraint into the linear system and construct the augmented system;

[0015] Step 5: Applying shaping regularization to constrain the density model, where the shaping operator is a triangular filter operator to achieve interpretable smoothing of the model space;

[0016] Step 6: Adaptively determine the regularization parameters using the Hoerl-Kennard-Baldwin (HKB) estimator, substitute them into the integer regularization equation to solve, and output the final denoised density model.

[0017] In a further technical solution, step 1 includes the following specific steps:

[0018] Gravity gradiometer for measuring gravitational field Given the three components of the gravitational field, there are a total of nine spatial gradients, which together form a gradient at every point in three-dimensional space. The tensor; let For gravitational potential, gravitational gradient tensor Derived from gravitational potential:

[0019] (1)

[0020] Among them, superscript The gradient operator is represented as a transpose; using a right-handed Cartesian coordinate system, the gravity gradient tensor is expressed as:

[0021] (2)

[0022] in, Pointing to geographical north, Pointing to the geographical east, Pointing downwards, Indicates gravitational potential along The second rate of change of direction, express directional gravitational component along Rate of change of direction express directional gravitational component along Rate of change of direction express directional gravitational component along Rate of change of direction Indicates gravitational potential along The second rate of change of direction, express directional gravitational component along Rate of change of direction express directional gravitational component along Rate of change of direction Indicates gravitational potential along The second rate of change of direction;

[0023] The gradient tensor is symmetric and has a zero trace, i.e.:

[0024] (3)

[0025] in, For coordinate indicators, the value is... .

[0026] In a further technical solution, step 2 includes the following specific steps:

[0027] A layer with lateral density variation is selected as the equivalent source layer at the observation horizon. The equivalent source layer is defined by the algebraic vector of relative density, thickness, burial depth, and horizontal dimension. The equivalent source layer is set at a certain distance below the observation surface, and the density layer is discretized into... Consider a set of rectangular prisms with approximately constant density within a given cell. Assuming each prism has a constant but unknown density value, the principle of linear superposition of gravitational fields is applied at any observation point. place The component gravity gradient can be written as a linear combination of the contributions of each element:

[0028] (4)

[0029] in, Indicates the observation point place Component gravity gradient, For the first The remaining density of each unit, The kernel function represents the density of the prism produced at the observation point. Component response; The six-component data from each observation point are concatenated in a predetermined order to obtain the matrix form of the full tensor forward modeling equation:

[0030] (5)

[0031] in, This is the sensitivity matrix; For density model vectors, denoted as ; Let the observed data vector be represented as ,in For the number of observation points, Represents a data vector The The component, i.e., the first component One observation value, Indicates the first At each observation point value, Indicates the first At each observation point value, Indicates the first At each observation point value, Indicates the first At each observation point value, Indicates the first At each observation point value, Indicates the first At each observation point value, Indicates the first At each observation point value.

[0032] In a further technical solution, step 3 includes the following specific steps:

[0033] Introducing a diagonal weighted matrix To balance the contributions of each component, a robust scaling estimate based on the absolute deviation of the median is used to obtain the noise level:

[0034] (6)

[0035] in, Indicates the first The standard deviation of each component noise This indicates taking the median. Indicates the first The data vector of each gradient component across all observation points; the corresponding weights are taken as follows: ,in Represents the data weighting matrix; after whitening, the noise is approximately assumed to satisfy... ,in This represents the data noise vector, i.e., the error term contained in the observations. Indicates a normal distribution. Indicates the noise variance. Represent a The identity matrix, The length of the observed data vector is represented by the weighted data fitting term. ,in This represents the observation data vector.

[0036] A further technical solution involves incorporating the zero-trace condition as a soft constraint into the linear system in step 4; (Note: The original text contains some inconsistencies and unclear formatting. A more accurate translation would require the full context.) For the sensitivity submatrix of the corresponding component, weighting coefficients are introduced. Controlling constraint strength to construct an augmented system:

[0037] (7).

[0038] In a further technical solution, step 5 includes the following specific steps:

[0039] Traditional Tikhonov regularization stabilizes the inversion by introducing a roughness penalty:

[0040] (8)

[0041] in, Let be the objective function. This represents the weighted and augmented orthogonal matrix. This represents the weighted and augmented vector of observed data. For roughness operators, Here is the regularization parameter; the corresponding normal equation is:

[0042] (9)

[0043] Introducing Integer Operators To make it approximate :

[0044] (10)

[0045] in, Representing the identity matrix, substituting it into the normal equation yields:

[0046] (11).

[0047] Further technical solutions, shaping operators It is implemented using a triangular filter.

[0048] In a further technical solution, step 6 includes the following steps:

[0049] First, obtain the least squares estimate without regularization. And calculate its residual variance estimate:

[0050] (12)

[0051] in To increase the data dimensionality of the system, Let be the number of cuboids and prisms whose density is approximately constant within the cell. This represents the estimated variance of the observation noise; subsequently, the regularization parameter is given according to the HKB formula:

[0052] (13)

[0053] in, This represents the regularization parameter calculated using the HKB formula; ultimately, it will... Substitute the values ​​into the shaping regularization equation to solve the problem, thus completing the inversion process.

[0054] The present invention provides an equivalent source denoising method based on shaping regularization, the beneficial effects of which are as follows:

[0055] (1) Improve the consistency and stability among components: Through robust scale normalization, the differences in the dimensions and noise levels of each component are eliminated, so that multiple components can participate in the inversion under a unified metric, avoiding the dominance of individual components in the results and ensuring the coordination and consistency of the denoising results among the components.

[0056] (2) Enhance the physical rationality of the results: By introducing zero-trace soft constraints, the physical laws of the gravitational field are integrated into the denoising process, which effectively suppresses the incoordination between non-physical components and components induced by noise, reduces the probability of structured residuals and pseudo-anomalies, and improves the physical interpretability of the denoising results.

[0057] (3) Preserve the anomaly shape and amplitude information: Using shaping regularization, interpretable smoothing operators (such as triangular filtering) are used to constrain the shape of the solution space. While suppressing noise, the spatial shape continuity and amplitude information of the anomaly are better preserved, and the loss of detail caused by over-smoothing is reduced.

[0058] (4) Improve processing efficiency and automation level: The regularization parameters are adaptively determined by the HKB method, avoiding tedious manual trial and error and the computationally expensive L-curve method, making the denoising process more efficient and objective, and having good adaptability and repeatability under different noise conditions. Attached Figure Description

[0059] Figure 1 A flowchart of an equivalent source denoising method based on shaping regularization provided in an embodiment of the present invention;

[0060] Figure 2 The original 3D model diagram (where a is the Hxx component, b is the Hxy component, c is the Hxz component, d is the Hyy component, e is the Hyz component, and f is the Hzz component).

[0061] Figure 3 A schematic diagram of the 3D model after adding random noise (where a is the Hxx component, b is the Hxy component, c is the Hxz component, d is the Hyy component, e is the Hyz component, and f is the Hzz component).

[0062] Figure 4 The denoised full tensor gravity gradient data of the 3D model (where a is the Hxx component, b is the Hxy component, c is the Hxz component, d is the Hyy component, e is the Hyz component, and f is the Hzz component).

[0063] Figure 5 The noise profile removed from the 3D model (where a is the Hxx component, b is the Hxy component, c is the Hxz component, d is the Hyy component, e is the Hyz component, and f is the Hzz component).

[0064] Figure 6 A schematic diagram after adding additional noise to the original survey area data (where a is the Hxx component, b is the Hxy component, c is the Hxz component, d is the Hyy component, e is the Hyz component, and f is the Hzz component).

[0065] Figure 7 The full tensor gravity gradient data after denoising the data of the survey area (where a is the Hxx component, b is the Hxy component, c is the Hxz component, d is the Hyy component, e is the Hyz component, and f is the Hzz component).

[0066] Figure 8The noise profile of the survey area data (where a is the Hxx component, b is the Hxy component, c is the Hxz component, d is the Hyy component, e is the Hyz component, and f is the Hzz component). Detailed Implementation

[0067] To make the objectives, technical solutions, and advantages of this invention clearer, the invention will be further described in detail below with reference to the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are merely illustrative and not intended to limit the invention.

[0068] The specific implementation of the present invention will be described in detail below with reference to specific embodiments.

[0069] like Figure 1 As shown, an equivalent source denoising method based on integer regularization is provided in one embodiment of the present invention. The specific implementation steps are as follows:

[0070] Step 1: Acquisition and representation of full tensor gravity gradient data;

[0071] Gravity gradiometer for measuring gravitational field Given the three components of the gravitational field, there are a total of nine spatial gradients, which together form a gradient at every point in three-dimensional space. The tensor of . Let . It is the gravitational potential, the gravitational gradient tensor. Derived from gravitational potential:

[0072] (1)

[0073] Among them, superscript Let represent the transpose of the gradient operator. Equation (1) shows that each component of the gravity gradient tensor is given by the second derivative of . A right-handed Cartesian coordinate system is used, where the x-axis points to geographic north, the y-axis points to east, and the z-axis points downwards. The gravity gradient tensor is expressed as:

[0074] (2)

[0075] in, Pointing to geographical north, Pointing to the geographical east, Pointing downwards, Indicates gravitational potential along The second rate of change of direction, express directional gravitational component along Rate of change of direction express directional gravitational component along Rate of change of direction express directional gravitational component along Rate of change of direction Indicates gravitational potential along The second rate of change of direction, express directional gravitational component along Rate of change of direction express directional gravitational component along Rate of change of direction Indicates gravitational potential along The second rate of change of direction.

[0076] Since the gravitational potential is a smooth function and satisfies the Laplace equation in the passive observation region, the gradient tensor is symmetric and has a zero trace, i.e.:

[0077] (3)

[0078] in, For coordinate indicators, the values ​​can be... .

[0079] Therefore, only five gradient components are linearly independent. Although five independent components can be obtained at each observation point, they all originate from the second derivative of the same potential function, essentially representing observations and projections of the same potential field in different directions. Thus, these components are coupled with each other.

[0080] Step 2: Construct an equivalent source forward model;

[0081] A layer with lateral density variation is selected as the equivalent source layer at the observation horizon. The equivalent source layer is defined by the algebraic vector of relative density, thickness, burial depth, and horizontal dimensions. The thickness, burial depth, and horizontal dimensions can be adaptively determined based on the existing data quality and the characteristics of the problem to be solved. The equivalent source layer is placed at a certain distance below the observation surface, and the density layer is discretized into... Consider a set of rectangular prisms with approximately constant density within a single cell, assuming each prism has a constant but unknown density value. Using the principle of linear superposition of gravitational fields, at any observation point... place The component gravity gradient can be written as a linear combination of the contributions of each element:

[0082] (4)

[0083] in, Indicates the observation point place Component gravity gradient, For the first The remaining density of each unit, The kernel function represents the density of the prism produced at the observation point. Component response. The kernel function is calculated using a prism analytical expression, which avoids the discrete errors and stability problems caused by numerical integration. The six-component data from each observation point are concatenated in a predetermined order to obtain the matrix form of the full tensor forward modeling equation:

[0084] (5)

[0085] in, This is the sensitivity matrix; For density model vectors, denoted as ; Let the observed data vector be represented as ,in The number of observation points. Represents a data vector The The component, i.e., the first component One observation value, Indicates the first At each observation point value, Indicates the first At each observation point value, Indicates the first At each observation point value, Indicates the first At each observation point value, Indicates the first At each observation point value, Indicates the first At each observation point value, Indicates the first At each observation point value.

[0086] Step 3: Robust scalar normalization;

[0087] The amplitude ranges of the components of the full tensor are usually inconsistent with the noise variance. If the residuals are directly and equally weighted in the same objective function, the components with larger amplitudes will have higher weights in the least squares sense; at the same time, the components with larger noise variance may be overfitted due to their larger residuals, thus causing the inversion results to deviate from the true model. Therefore, a diagonal weighting matrix is ​​introduced. To balance the contributions of each component, a robust scaling estimate based on median absolute deviation (MAD) is used to obtain the noise level.

[0088] (6)

[0089] in, Indicates the first The standard deviation of each component noise This indicates taking the median. Indicates the first The data vector of each gradient component over all observation points. The corresponding weights are taken as follows: ,in This represents a data weighting matrix, thereby achieving equal-weighted fusion of different components in a statistical sense and reducing the impact of outliers on the inversion. After whitening, the noise can be approximately considered to satisfy... ,in This represents the data noise vector, i.e., the error term contained in the observations. Indicates a normal distribution. Indicates the noise variance. Represent a The identity matrix, This represents the length of the observed data vector. The weighted data fitting term is written as... ,in Represents a diagonal weighted matrix. This represents the observation data vector.

[0090] Step 4: Introduce zero-trace physical constraints;

[0091] To enhance physical consistency during inversion, this invention incorporates the zero-trace condition as a soft constraint into the linear system. For the sensitivity submatrix of the corresponding component, weighting coefficients are introduced. Controlling constraint strength to construct an augmented system:

[0092] (7)

[0093] This processing is equivalent to penalizing "trace close to zero" at the data level, which helps to suppress the inhomogeneous components caused by noise or inconsistent components and improves the stability of the full tensor joint inversion.

[0094] Step 5: Solving for the shaping regularization constraint and adaptive parameter selection;

[0095] Illness makes direct minimization possible This can easily lead to high-frequency oscillations and non-physical solutions. Traditional Tikhonov regularization stabilizes the inversion by introducing a roughness penalty:

[0096] (8)

[0097] in, The objective function is defined as follows: given a density model vector... Use this formula to measure quality This represents the weighted and augmented orthogonal matrix. This represents the weighted and augmented vector of observed data. For roughness operators (usually the first-order gradient or the second-order Laplace). The regularization parameter controls the trade-off between data fitting and smoothing constraints. The corresponding normal equation is:

[0098] (9)

[0099] In large-scale problems, Explicit construction and solution are costly, and at the same time Choice and The adjustment of roughness often fails to directly correspond to the desired geological morphological features. The idea behind shaping regularization is to replace explicit roughness penalties with "interpretable smooth mappings." This involves introducing a shaping operator. To make it approximate :

[0100] (10)

[0101] in Representing the identity matrix, substituting it into the normal equation yields:

[0102] (11)

[0103] This invention selects triangular filtering as... The implementation of triangular filtering can be viewed as a quadratic window smoothed convolution, characterized by zero phase and controllable smoothness. Its function is equivalent to applying an interpretable correlation length to the model space, thus directly translating the "smoothing level" into a geologically understandable assumption of structural continuity. By adjusting the filtering radius, a more geologically sound compromise can be achieved between suppressing noise and preserving anomaly boundaries.

[0104] Step 6: Output the results;

[0105] Regularization parameters The trade-off between data fitting and model smoothing is a key factor affecting the inversion results. To reduce the subjectivity of manual trial and error, this invention uses the Hoerl-Kennard-Baldwin (HKB) estimator for adaptive determination. This method originates from the parameter estimation idea of ​​ridge regression: it uses the statistical properties of unregularized least squares to derive a reasonable regularization strength.

[0106] First, obtain the least squares estimate without regularization. And calculate its residual variance estimate:

[0107] (12)

[0108] in To increase the data dimensionality of the system, Let be the number of cuboids and prisms whose density is approximately constant within the cell. This represents the estimated variance of the observation noise. The regularization parameter is then given according to the HKB formula:

[0109] (13)

[0110] in, This represents the regularization parameter calculated using the HKB formula. This estimate statistically links the noise energy to the model energy, making... It can adaptively adjust to the level of data noise, thereby suppressing unstable components while preserving as much effective anomaly information as possible. Ultimately, it will... Substitute the values ​​into the shaping regularization equation to solve the problem, thus completing the inversion process.

[0111] The effectiveness of this method was tested using a three-dimensional model consisting of a sphere and a cube. The sphere has a radius of 20 m, its center is located at (1700, 650, 50), and its density is 1500 kg / m³. 3 The cube has a side length of 50m and a density of 1800 kg / m³. 3 Rotating 30° around the z-axis, the distance between the center of the sphere and the center of the cube is 100 m. This invention uses the original data ( Figure 2 Add random noise interference (such as) to the basis of ) Figure 3 This results in a low signal-to-noise ratio in the data, and the denoising result is as follows: Figure 4 As shown in the results, the signal-to-noise ratio of all components is significantly improved, and the anomalous characteristic relationships of the six components are well preserved. From the removed noise profile ( Figure 5 As can be seen in the data, no effective signal appears in the noise profile, which verifies that the method can effectively separate the signal from the noise in low signal-to-noise ratio multi-component gravity gradient data.

[0112] The effectiveness of this method was tested using open-source gravity gradient data. The survey area was Sullivan North (Block 1), Missouri, USA, with its center at approximately 91°13'W, 38°09'N. The area contains three iron oxide deposits and covers an area of ​​approximately 35 × 37 km². 2 The data was collected in March 2014. The north-south survey line spacing was 400 m, the east-west connecting line was 4000 m, and the flight altitude was 80 m. The noise reduction results are as follows... Figure 7As shown in the results, the signal-to-noise ratio of all components is significantly improved, and the denoising profile ( Figure 8 The absence of a valid signal indicates that this invention is an effective signal-to-noise separation method.

[0113] The above description is only a preferred embodiment of the present invention and is not intended to limit the present invention. Any modifications, equivalent substitutions, and improvements made within the spirit and principles of the present invention should be included within the protection scope of the present invention.

Claims

1. An equivalent source denoising method based on integer regularization, characterized in that, Includes the following steps: Step 1: Obtain the full tensor gravity gradient data and express it as a vector; Step 2: Construct an equivalent source forward model, discretize the equivalent source layer into multiple prism elements, and establish a linear forward modeling equation from element density to observation data; Step 3: Perform robust scaling normalization on each component of the full tensor gravity gradient data, estimate the noise standard deviation based on the median absolute deviation of each component, and construct a diagonal weighted matrix to balance the contribution of each component. Step 4: Utilize the zero-trace property of the gravity gradient tensor to introduce the zero-trace condition as a soft constraint into the linear system and construct the augmented system; Step 5: Applying shaping regularization to constrain the density model, where the shaping operator is a triangular filter operator to achieve interpretable smoothing of the model space; Step 6: Use the HKB estimator to adaptively determine the regularization parameters, substitute them into the integer regularization equation to solve, and output the final denoised density model.

2. The equivalent source denoising method based on shaping regularization according to claim 1, characterized in that, Step 1 includes the following specific steps: Gravity gradiometer for measuring gravitational field Given the three components of the gravitational field, there are a total of nine spatial gradients, which together form a gradient at every point in three-dimensional space. The tensor; let For gravitational potential, gravitational gradient tensor Derived from gravitational potential: （1） Among them, superscript The gradient operator is represented as a transpose; using a right-handed Cartesian coordinate system, the gravity gradient tensor is expressed as: （2） in, Pointing to geographical north, Pointing to the geographical east, Pointing downwards, Indicates gravitational potential along The second rate of change of direction, express directional gravitational component along Rate of change of direction express directional gravitational component along Rate of change of direction express directional gravitational component along Rate of change of direction Indicates gravitational potential along The second rate of change of direction, express directional gravitational component along Rate of change of direction express directional gravitational component along Rate of change of direction Indicates gravitational potential along The second rate of change of direction; The gradient tensor is symmetric and has a zero trace, i.e.: （3） in, For coordinate indicators, the value is... .

3. The equivalent source denoising method based on shaping regularization according to claim 2, characterized in that, Step 2 includes the following specific steps: A layer with lateral density variation is selected at the observation horizon as the equivalent source layer. This equivalent source layer is then positioned at a specified distance below the observation surface, and the density layer is discretized. Consider a set of rectangular prisms with approximately constant density within a single cell, assuming each prism has a constant but unknown density value. Using the principle of linear superposition of gravitational fields, at any observation point... place The component gravity gradient can be written as a linear combination of the contributions of each element: （4） in, Indicates the observation point place Component gravity gradient, For the first The remaining density of each unit, The kernel function represents the density of the prism produced at the observation point. Component response; The six-component data from each observation point are concatenated in a predetermined order to obtain the matrix form of the full tensor forward modeling equation: （5） in, This is the sensitivity matrix; For density model vectors, denoted as ; Let the observed data vector be represented as ,in For the number of observation points, Represents a data vector The The component, i.e., the first component One observation value, Indicates the first At each observation point value, Indicates the first At each observation point value, Indicates the first At each observation point value, Indicates the first At each observation point value, Indicates the first At each observation point value, Indicates the first At each observation point value, Indicates the first At each observation point value.

4. The equivalent source denoising method based on shaping regularization according to claim 3, characterized in that, Step 3 includes the following specific steps: Introducing a diagonal weighted matrix To balance the contributions of each component, a robust scaling estimate based on the absolute deviation of the median is used to obtain the noise level: （6） in, Indicates the first The standard deviation of each component noise This indicates taking the median. Indicates the first The data vector of each gradient component across all observation points; the corresponding weights are taken as follows: ,in Represents the data weighting matrix; after whitening, the noise is approximately assumed to satisfy... ,in This represents the data noise vector, i.e., the error term contained in the observations. Indicates a normal distribution. Indicates the noise variance. Represent a The identity matrix, This represents the length of the observed data vector; the weighted data fitting term is written as... ,in This represents the observation data vector.

5. The equivalent source denoising method based on shaping regularization according to claim 4, characterized in that, In step 4, the zero-trace condition is incorporated as a soft constraint into the linear system; let... For the sensitivity submatrix of the corresponding component, weighting coefficients are introduced. Controlling constraint strength to construct an augmented system: （7）。 6. The equivalent source denoising method based on shaping regularization according to claim 5, characterized in that, Step 5 includes the following specific steps: Traditional Tikhonov regularization stabilizes the inversion by introducing a roughness penalty: （8） in, Let be the objective function. This represents the weighted and augmented orthogonal matrix. This represents the weighted and augmented vector of observed data. For roughness operators, Here is the regularization parameter; the corresponding normal equation is: （9） Introducing Integer Operators To make it approximate : （10） in, Representing the identity matrix, substituting it into the normal equation yields: （11）。 7. The equivalent source denoising method based on shaping regularization according to claim 6, characterized in that, Plasticizer It is implemented using a triangular filter.

8. The equivalent source denoising method based on shaping regularization according to claim 6, characterized in that, Step 6 includes the following steps: First, obtain the least squares estimate without regularization. And calculate its residual variance estimate: （12） in To increase the data dimensionality of the system, Let be the number of cuboid prisms whose density is approximately constant within the cell. This represents the estimated variance of the observation noise; subsequently, the regularization parameter is given according to the HKB formula: （13） in, This represents the regularization parameter calculated using the HKB formula; ultimately, it will... Substitute the values ​​into the regularization equation to solve the problem, thus completing the inversion process.

Images