A method and system for downward extension of gravity data based on compact constraints

Abstract

The present invention relates to a method and system for downward extension of gravity data based on compact constraints, comprising: obtaining gravity data from each measuring point in a study area and preprocessing it; determining the geometric parameters of an equivalent source in the study area; obtaining the density distribution of the equivalent source based on the preprocessed gravity data in combination with a compact constraint term and a smoothing strategy; determining the compact constraint term based on the equivalent source density and the compact factor; obtaining an equivalent source model based on the geometric parameters of the gravity equivalent source and the density of the equivalent source; and performing forward calculations based on the equivalent source model to obtain the downward extension results of the gravity data. By introducing a compact constraint term, the equivalent source physical properties are concentrated, and the smoothing strategy is used to suppress or eliminate divergent signals in the equivalent source physical properties, thereby significantly improving the distortion of the data morphology, improving the accuracy of data amplitude recovery, and ensuring the stability of downward extension at great depths.

Description

Technical Field

[0001] The present invention relates to the fields of gravity and gravity exploration, and in particular to a method and system for downward extension of gravity data based on compact constraints. Background Art

[0002] Downward continuation is a crucial aspect of gravity data processing. Converting observational data from the observation surface to the passive space below it can improve source resolution, enhance weak signals, and meet the demand for gravity field information at different altitudes. Examples include downcontinuing airborne gravity data to the surface, downconverting ship-derived gravity data to the seafloor thousands of meters deep, and downconverting satellite gravity data from hundreds of kilometers altitude to near-surface altitudes. The high-precision and high-resolution gravity data obtained through downward continuation plays a vital role in fields such as geoscience, oceanography, environmental science, resource exploration, and national defense security.

[0003] Stable, accurate, and deep downward continuation has always been a difficult and hot issue in gravity data processing. Stability is the primary factor limiting the maximum distance and accuracy of downward continuation. Downward continuation methods can be categorized into two types: frequency-domain and spatial-domain methods. These include methods based on Taylor series expansion, Fourier transform, analytical continuation, integral-iteration methods, and solutions based on potential field median theory. These methods, either by solving high-order derivatives or by the amplification effect of continuation factors, amplify high-frequency signals, particularly noise (or spurious interference), which distorts the data morphology and can even drown out valid signals, leading to instability in the downward continuation process. Common approaches to improve stability include regularization, filtering, increasing the number of iterations, and multi-source data constraints. While regularization has been shown to suppress noise, it also compromises signal resolution and amplitude. Filtering has similar effects and applications to regularization. Increasing the number of iterations can mitigate over-amplification of high-frequency signals, but it still fails to fundamentally address the issue of noise. Mutual constraints between multi-source data are currently a popular approach, but in many cases, they are often hindered by the inability to obtain high-quality multi-source data. Therefore, while previous research has actively promoted the development of downward continuation methods, the instability of downward continuation still severely restricts the continuation distance and accuracy, necessitating continued in-depth research.

[0004] The equivalent source method is a commonly used technique in gravity data processing. This method utilizes a set of artificially constructed virtual sources (equivalent sources) to replace real sources and generate gravity data in a sourceless space. Specifically, an equivalent source model is established by fitting the observed data, transforming gravity data processing into forward calculations using the equivalent source model. This method has the advantages of directly processing the original point observation data and taking into account the homology between data, demonstrating excellent performance in gravity data processing. However, when using this method for deep downward continuation, data morphological distortion and amplitude loss can also occur. The root cause is the presence of "divergent" signals in the physical properties of the equivalent sources. Therefore, obtaining high-precision and high-resolution downward continuation results places higher demands on the equivalent source model, necessitating the exploration of new technical methods to ensure the stability of deep downward continuation. Summary of the Invention

[0005] The purpose of the present invention is to provide a method and system for downward extension of gravity data based on compact constraints. By introducing compact constraints, the equivalent source properties are concentratedly distributed. At the same time, a smoothing strategy is considered to suppress or eliminate the "divergent" signals in the equivalent source properties, which significantly improves the problem of data morphological distortion, improves the recovery accuracy of data amplitude, and ensures the stability of downward extension at a large depth.

[0006] To achieve the above object, the present invention provides the following solutions:

[0007] A gravity data downward continuation method based on compact constraints, comprising:

[0008] Obtaining gravity data of each measuring point in the study area and preprocessing the gravity data;

[0009] Determine the geometric parameters of the equivalent source in the study area;

[0010] The density distribution of the equivalent source is obtained based on the preprocessed gravity data in combination with a compact constraint term and a smoothing strategy; the compact constraint term is determined according to the equivalent source unit density and the compact factor;

[0011] Obtaining an equivalent source model of the study area according to the geometric parameters of the equivalent source and the density distribution of the equivalent source;

[0012] A forward calculation is performed based on the equivalent source model to obtain a downward extension result of the gravity data.

[0013] Optionally, the density distribution of the equivalent source is obtained based on the preprocessed gravity data in combination with a compact constraint term and a smoothing strategy, specifically including:

[0014] establishing an inversion objective function based on the preprocessed gravity data, the equivalent source forward modeling kernel matrix, and the compact constraint term;

[0015] Obtaining an inversion calculation equation according to the inversion objective function;

[0016] The equivalent source density is obtained by performing inversion calculation according to the inversion calculation equation in combination with the preconditioned conjugate gradient method, and the obtained equivalent source density is smoothed to obtain a smoothed equivalent source density.

[0017] Optionally, the inversion objective function is expressed as:

[0018]

[0019] Where d represents the preprocessed gravity data; G represents the equivalent source forward kernel matrix; m represents the density of the equivalent source to be determined; μ represents the regularization parameter; m T W T Wm represents the compact constraint term; w j Represents the diagonal elements of the diagonal matrix W in the constraint terms; m′ j is the density of the jth equivalent source unit, α is the compact factor, and ε is 10 -7 .

[0020] Optionally, the inversion calculation equation is expressed as:

[0021] P(G T G+μW T W)m=PG T d

[0022] Where P is the preconditioning matrix, P(G T G+μW T W)≈I, the preconditioning matrix takes the form of a diagonal matrix, and the diagonal elements of the preconditioning matrix G ij is the element of the equivalent source forward modeling kernel matrix G; i represents the measurement point code; M represents the number of measurement points.

[0023] Optionally, performing inversion calculation according to the inversion calculation equation of the equivalent source density in combination with a preconditioned conjugate gradient method to obtain an equivalent source density, and smoothing the obtained equivalent source density to obtain a smoothed equivalent source density, specifically includes:

[0024] Let (G T G+μW T W)=A,G T d = B;

[0025] Set the initial index of the equivalent source density to m0 = 0, the preset error value to e, r0 = A T B; n = 0, 1, ...;

[0026] Let y n =Prn ; When n = 0, S n =y n ; When n≠0, S n =y n +β n-1 S n-1 ;

[0027] Determine the iteration step size And according to the formula m n =m n-1 +t n S n Update equivalent source density;

[0028] according to Determine whether the iteration is finished;

[0029] When not satisfied When n=n+1, calculate r n+1 =r n -t n A T AS n , and return to step "Let y n =Pr n ; When n = 0, S n =y n ; When n≠0, S n =y n +β n-1 S n-1 ”;

[0030] When satisfied When , the iteration ends, and the equivalent source density m after the nth iteration is obtained. n ;

[0031] For the equivalent source density m n Smoothing is performed to obtain the smoothed equivalent source density.

[0032] Optionally, the equivalent source density m n Performing a smoothing process to obtain the smoothed equivalent source density specifically includes:

[0033] According to the equivalent source density m n Performing divergent signal detection and smoothing the detected divergent signal;

[0034] The smoothed equivalent source density is used as the initial value m0 of the equivalent source density, and n = 0 is initialized, and the process returns to step "Let y n =Pr n ; When n = 0, S n =y n ; When n≠0, S n=y n +β n-1 S n-1 "; until the divergent signal and (d-Gm) do not exist in the current smoothed equivalent source density T The value of (d-Gm) is less than the preset value, and the final result of the equivalent source density m is obtained. inv .

[0035] The present invention also provides a gravity data downward continuation system based on compact constraints, comprising:

[0036] A data acquisition and processing module is used to acquire gravity data of each measuring point in the study area and pre-process the gravity data;

[0037] Equivalent source geometric parameter determination module, used to determine the geometric parameters of the equivalent source in the study area;

[0038] An equivalent source density calculation module is used to obtain the density distribution of the equivalent source based on the preprocessed gravity data in combination with a compact constraint term and a smoothing strategy; the compact constraint term is determined according to the equivalent source density and the compact factor;

[0039] An equivalent source model construction module, configured to obtain an equivalent source model of the study area according to geometric parameters of the equivalent source and density distribution of the equivalent source;

[0040] A forward calculation module is used to perform forward calculation based on the equivalent source model to obtain a downward continuation result of the gravity data.

[0041] Optionally, the equivalent source density calculation module specifically includes:

[0042] An objective function establishing unit is used to establish an inversion objective function according to the preprocessed gravity data, the equivalent source forward modeling kernel matrix and the compact constraint term;

[0043] An inversion calculation equation construction unit, configured to obtain an inversion calculation equation according to the inversion objective function;

[0044] The inversion and smoothing processing unit is used to perform inversion calculation according to the inversion calculation equation combined with the preconditioned conjugate gradient method to obtain the equivalent source density, and to smooth the equivalent source density to obtain the smoothed equivalent source density.

[0045] Optionally, the objective function of the gravity equivalent source density is expressed as:

[0046]

[0047] Where d represents the preprocessed gravity data; G represents the equivalent source forward kernel matrix; m represents the density of the equivalent source to be determined; μ represents the regularization parameter; m TW T Wm represents the compact constraint term; w j Represents the diagonal elements of the diagonal matrix W in the constraint terms; m′ j is the density of the jth equivalent source unit, α is the compact factor, and ε is 10 -7 .

[0048] Optionally, the inversion calculation equation of the gravity equivalent source density is expressed as:

[0049] P(G T G+μW T W)m=PG T d

[0050] Where P is the preconditioning matrix, P(G T G+μW T W)≈I, the preconditioning matrix takes the form of a diagonal matrix, and the diagonal elements of the preconditioning matrix G ij is the element of the equivalent source forward modeling kernel matrix G; i represents the measurement point code, and M represents the number of measurement points.

[0051] According to the specific embodiments provided by the present invention, the present invention discloses the following technical effects:

[0052] The present invention relates to a method and system for downward extension of gravity data based on compact constraints, comprising: obtaining gravity data from each measuring point in a study area and preprocessing it; determining the geometric parameters of equivalent sources in the study area; obtaining the density distribution of the equivalent sources based on the preprocessed gravity data in combination with compact constraints and a smoothing strategy; determining the compact constraints based on the equivalent source density and the compact factor; obtaining an equivalent source model based on the geometric parameters and density of the equivalent source; and performing forward calculations based on the equivalent source model to obtain the downward extension results of the gravity data. By introducing compact constraints, the physical properties of the equivalent sources are concentrated, and the smoothing strategy is used to suppress or eliminate divergent signals in the equivalent source density, the distortion of the data morphology is significantly improved, the accuracy of data amplitude recovery is improved, and the stability of downward extension at great depths is ensured. BRIEF DESCRIPTION OF THE DRAWINGS

[0053] In order to more clearly illustrate the embodiments of the present invention or the technical solutions in the prior art, the following briefly introduces the drawings required for use in the embodiments. Obviously, the drawings described below are only some embodiments of the present invention. For ordinary technicians in this field, other drawings can be obtained based on these drawings without paying any creative work.

[0054] Figure 1 A flow chart of a method for downward extension of gravity data based on compact constraints provided in Example 1 of the present invention;

[0055] Figure 2 A technical principle diagram of a gravity data downward continuation method based on compact constraints provided in Example 1 of the present invention;

[0056] Figure 3 This is a schematic diagram of an equivalent source provided in Example 1 of the present invention;

[0057] Figure 4 Schematic diagram showing how the number of iterations varies with μ, provided in Example 1 of the present invention;

[0058] Figure 5 The theoretical model and forward gravity anomaly data provided in Example 1 of the present invention;

[0059] Figure 6 A comparison chart of the downward extension results of the 800m altitude data provided in Example 1 of the present invention;

[0060] Figure 7 The noise-containing and original gravity anomaly data at an altitude of 800 m provided in Example 1 of the present invention;

[0061] Figure 8 This is a comparison chart of the downward extension results of the noisy data at an altitude of 800m provided in Example 1 of the present invention. DETAILED DESCRIPTION

[0062] The following will clearly and completely describe the technical solutions in the embodiments of the present invention in conjunction with the accompanying drawings. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without making creative efforts are within the scope of protection of the present invention.

[0063] The purpose of the present invention is to provide a method and system for downward extension of gravity data based on compact constraints. By introducing compact constraints, the equivalent source properties are concentratedly distributed. At the same time, a smoothing strategy is considered to suppress or eliminate the "divergent" signals in the equivalent source properties, which significantly improves the problem of data morphological distortion, improves the recovery accuracy of data amplitude, and ensures the stability of downward extension at a large depth.

[0064] In order to make the above-mentioned objects, features and advantages of the present invention more obvious and easy to understand, the present invention is further described in detail below with reference to the accompanying drawings and specific embodiments.

[0065] Example 1

[0066] like Figure 1 and 2 As shown, this embodiment provides a method for downward extension of gravity data based on compact constraints, including:

[0067] S1: Obtain gravity data of each measuring point in the study area and preprocess the gravity data.

[0068] This embodiment works on Bouguer gravity anomaly data or Bouguer gravity anomaly gradient data; the data information includes: easting, northing and vertical coordinates of the Cartesian coordinate system, data value, and total data accuracy.

[0069] The process of screening (preprocessing) the data is as follows:

[0070] (1) The Bouguer gravity anomaly data or Bouguer gravity anomaly gradient data at discrete measuring points in the survey line or survey area are mapped (contour map or image map), and the sudden jump points or sharp points in the data are eliminated based on the changing trend and basic laws of the data in the map.

[0071] (2) According to the requirements for data resolution (determined in advance by the task requirements), the data after removing the sudden jump points or sharp points are resampled to remove redundant data, that is, the study area is gridded according to the resolution, and one measuring point is randomly retained in each grid through which the survey line passes; the resampled data is used as the final data for calculation (that is, the preprocessed gravity data used for calculation in step S3); the purpose of data resampling is to reduce the amount of calculation, avoid overfitting, and provide a reference for the horizontal size of the equivalent source unit.

[0072] S2: Determine the geometric parameters of the equivalent source in the study area.

[0073] (1) The equivalent source is equivalent to the real field source. The real field source is all underground materials, such as strata, rocks, minerals, etc., which have mass and density. The equivalent source is the earth or part of the earth. The equivalent source is composed of a set of discrete rectangular blocks with uniform density. The rectangular blocks are arranged in close proximity, and the geometric dimensions of the rectangular blocks are determined by the spacing (or average spacing) of the reference data points.

[0074] (2) Depending on the severity of the terrain, the equivalent source can be set on a curved surface or a flat surface within a certain depth range below the surface.

[0075] like Figure 3As shown, the equivalent source unit adopts a rectangular parallelepiped with uniform density, and the horizontal position of the unit body is determined according to the position of the measuring point and its distribution range; in the absence of geological or geophysical prior information, the top surface of the unit body is buried within the range of 2.5 to 6 times the average measuring point spacing below the surface. In actual operation, it is recommended to place the equivalent source as deep as possible (greater than 4 times the point spacing); when prior information is available, the top surface burial depth of the equivalent source can be determined based on the average depth of the field source reflected by the prior information; the horizontal size of the equivalent source unit is determined according to the gravity data resolution requirement in step S1 (2), and the unit body thickness is 3 to 6 times its minimum horizontal size; when the terrain is undulating, the equivalent source is placed on a curved surface parallel to the terrain, otherwise the equivalent source is placed on a plane; the equivalent source layer is expanded outward from the data coverage area by 3 to 5 equivalent source units to suppress the boundary effect.

[0076] S3: Based on the preprocessed gravity data, a compact constraint term and a smoothing strategy are combined to obtain the density distribution of the equivalent source; the compact constraint term is determined according to the equivalent source density and the compact factor.

[0077] Wherein, step S3 specifically includes:

[0078] S31: establishing an inversion objective function according to the preprocessed gravity data, the equivalent source forward modeling kernel matrix and the compact constraint term.

[0079] Specifically,

[0080] According to the geometric parameters of the equivalent source and the spatial position coordinates of the pre-processed data d, the equivalent source forward kernel matrix G = {G ij}, its forward calculation formula adopts the publicly released analytical formula for cuboid forward calculation to establish the inversion objective function:

[0081] φ=φ d +μφ m (1)

[0082] Where, φ d That is, the variance between the observed data and the model calculated value (d-Gm) T (d-Gm), m is the equivalent source density vector to be determined; the second term φ m is a compact constraint term. The introduction of constraint information through this term can affect the distribution of equivalent source density and is expressed in matrix form as:

[0083] φ m =m T W T Wm (2)

[0084] Where W is a diagonal matrix with the following diagonal elements:

[0085]

[0086] Where m′ j is the jth equivalent source unit density, which belongs to the element in the equivalent source density vector m to be determined, α is the compact factor, and ε is a very small number with a value of 10 -7 The equivalent source units w determined by formula (3) j Construct a diagonal matrix W; bring them into the objective function (1), and we can get

[0087] The expression of the inversion objective function is:

[0088]

[0089] Where d represents the preprocessed gravity data; G represents the equivalent source forward kernel matrix; m represents the equivalent source density to be calculated; μ represents the regularization parameter; m T W T Wm represents the compact constraint term φ m ; w j represents the diagonal elements of the diagonal matrix W; m′ j is the density of the jth equivalent source unit, α is the compact factor, and ε is 10 -7 .

[0090] S32: Obtaining an inversion calculation equation according to the inversion objective function.

[0091] Find m for φ T Derivative, and set it to zero, that is, dφ / dm T =0, we can get

[0092] (G T G+μW T W)m=G T d (5)

[0093] Multiply both ends of equation (5) by the preconditioning matrix P to obtain the final equation for inversion calculation:

[0094] That is, the inversion calculation equation for obtaining the equivalent source density is:

[0095] P(G T G+μW T W)m=PG T d (6)

[0096] Where P is the preconditioning matrix, which is used to improve the condition number of the coefficient matrix and should satisfy P(G T G+μW T W)≈I, the preconditioning matrix takes the form of a diagonal matrix, and the diagonal elements of the preconditioning matrix G ijis the element of the equivalent source forward modeling kernel matrix G; i represents the measurement point code; M represents the number of measurement points.

[0097] S33: performing inversion calculation according to the inversion calculation equation in combination with the preconditioned conjugate gradient method to obtain an equivalent source density, and performing smoothing processing on the equivalent source density to obtain a smoothed equivalent source density.

[0098] Among them, S33 specifically includes:

[0099] S331: Let (G T G+μW T W)=A,G T d=B.

[0100] S332: Set the initial value of the equivalent source density to m0 = 0, the preset error value to e, and r0 = A T B. The initial value of n is 0.

[0101] S333: Let y n =Pr n ; When n = 0, S n =y n ; When n≠0, S n =y n +β n-1 S n-1 ; (This step corresponds to the updated search direction in the preconditioned conjugate gradient method)

[0102] S334: Determine the iteration step size And according to the formula m n =m n-1 +t n S n Update the equivalent source density.

[0103] S335: According to Determine whether the iteration is finished.

[0104] When not satisfied When , the next inversion iteration is performed, let n = n + 1, and calculate r n+1 =r n -t n A T AS n , and return to step S333 "let y n =Pr n ; When n = 0, S n =y n ; When n≠0, S n =y n +β n-1 S n-1 ”.

[0105] When satisfied When the inversion iteration ends, the equivalent source density m after the nth iteration is obtained. n .

[0106] S336: For the equivalent source density m n Smoothing is performed to obtain the smoothed equivalent source density.

[0107] Specifically, step S336 includes:

[0108] (1) According to the equivalent source density m n Diverging signals are detected and smoothed. After compaction, valid signals and "diverging" signals are visibly distinguishable. Diverging signals are characterized by a "ripple" shape around valid signals, with a smaller amplitude than the valid signal, and a significant difference between the two. The two can be distinguished by combining signal shape and amplitude (a signal with an amplitude less than a certain value is considered a diverging signal).

[0109] Check the distribution characteristics of the equivalent source density obtained by inversion and determine the distribution range of the residual "divergent" signal.

[0110] According to the determined distribution range of the "divergent" signal, the signal within the distribution range is set to zero or equal to a certain minimum number (ie, smoothed).

[0111] (2) The equivalent source density after smoothing the “divergent” signal is used as the initial value m0 of the equivalent source density, and n=0 is initialized, and the process returns to step S333 “Let y n =Pr n ; When n = 0, S n =y n ; When n≠0, S n =y n +β n-1 S n-1 "; until the current amplitude of the "divergence" signal is less than the preset amplitude and (d-Gm) T The value of (d-Gm) is less than the preset value, and the final result of the equivalent source density m is obtained. inv , that is, the equivalent source density after smoothing.

[0112] In this step, the smoothed equivalent source density is assigned as the initial value of the equivalent source density, and the inversion iteration process of step S33 is returned. The number of iterations starts from 0, and an equivalent source density is inverted again. Then, the smoothing process is performed again. The smoothed equivalent source density is assigned as the initial value of the equivalent source density again, and the inversion iteration process of step S333 is performed again until the amplitude of the divergent signal in the inverted equivalent source density is less than the preset amplitude and (d-Gm) TIf the value of (d-Gm) is less than the preset value, the final equivalent source density is output.

[0113] The "divergent" signal in the equivalent source properties can be significantly weakened by compact constraints, but some residues still remain. After compact constraints, the "divergent" signal is clearly distinguished from the effective signal. Therefore, after each inversion calculation is completed, smoothing measures are used to further eliminate the "divergent" signal, while enhancing the amplitude of the effective signal, that is, the amplitude of the "divergent" signal is returned to zero or less than a certain value. The smoothed equivalent source properties are used as the initial value of the equivalent source density for further inversion calculation, and this is iterated several times until the "divergent" signal is effectively suppressed.

[0114] In order to balance the accuracy and efficiency of the inversion calculation, it is necessary to reasonably determine the values ​​of the regularization parameter μ and the compactness factor α.

[0115] The role of the regularization parameter μ is to balance the influence of the first term (data fitting function) and the second term (constraint term) in the inversion objective function (1) on the equivalent source properties. Under the premise of determining the data fitting error, when μ is small, the emphasis is placed on fitting the data, and the role of the constraint term is not significant. When μ is large, although the role of the constraint information is well played, the number of iterations will increase, thereby reducing the computational efficiency. Therefore, in this embodiment, μ is set near the sudden jump point of the number of iterations, taking into account the constraint effect and computational efficiency.

[0116] The general rule for selecting the compactness factor α is that when its value is small, the number of inversion iterations is large and the distribution of equivalent source properties is smoother. When its value is large, the number of iterations is small, but the equivalent source property signal will experience large fluctuations. Therefore, when prior information is available, a larger α value is selected. When prior information is not available, the iteration starts with a smaller α value and then gradually increases the α value according to a certain rule until it reaches a certain value and remains unchanged. This embodiment adopts a method of giving α a small initial value (0-1) and then gradually increasing the α value with each iteration, for example, increasing it by 0.2 each time until it reaches 4 and remains unchanged.

[0117] Therefore, when determining the optimal value of μ, multiple random values ​​of μ can be determined in advance. These random values ​​are determined based on experience. These random values ​​are respectively substituted into formula (6) and participate in the inversion calculation respectively. That is, each μ is substituted into formula (6), and step S33 is executed to obtain a schematic diagram of the number of iterations changing with μ, as shown in FIG. Figure 4 As shown in FIG, the optimal value of μ is determined according to the sudden jump point of the iteration number in the schematic diagram of the variation of the iteration number with μ.

[0118] S4: Obtaining an equivalent source model of the study area according to the geometric parameters of the equivalent source and the density distribution of the equivalent source.

[0119] The equivalent source method uses a set of simple artificial sources (usually composed of discrete cuboids or point masses) to replace real sources and generate gravity field information in a source-free space. The general approach is to first manually specify the spatial location and size of the equivalent sources, namely the geometric parameters (step S2), at which point the equivalent source properties are unknown. Then, a system of equations is established based on the observed data. By fitting the observed data, the distribution of the equivalent source properties is inverted and calculated (step S3), thus establishing an equivalent source model. This allows the conversion and prediction of gravity data to be converted into forward calculations of the equivalent source model.

[0120] S5: Perform forward calculation based on the equivalent source model to obtain a downward extension result of the gravity data.

[0121] The forward modeling formula used in this step can be the uniform density cuboid forward modeling formula G′m inv = d′, at this time (the spatial coordinates of the data points in G′ and d′ should be replaced by the coordinates of the downward extension position). Where G′ is the forward kernel function used for downward extension, and d′ is the gravity data for downward extension (i.e., the downward extension result to be determined).

[0122] This embodiment, based on the basic principle of potential field equivalent source, optimizes the distribution of equivalent source properties by introducing compact constraints and, in combination with a smoothing strategy, suppresses or eliminates the "divergent" signals in the equivalent source properties while ensuring the amplitude of the effective signal, thereby achieving the goal of high-precision, deep downward extension of gravity data at a certain height in passive space. This invention can obtain a reliable three-dimensional passive space gravity database and establish a gravity field model for the region, thereby providing technical support and basic data for applications such as gravity research, gravity navigation reference maps, gravity exploration, and gravity instrument testing. Therefore, this embodiment has the following advantages:

[0123] (1) Starting from the perspective of constraining the distribution of physical properties of equivalent sources, it directly targets the fundamental problem that causes the instability of downward continuation. Compared with the traditional equivalent source method, its effect on improving the stability of downward continuation is more significant and effective.

[0124] (2) Combining compact constraints and smoothing strategies to extend gravity data downward significantly improves the suppression effect of distorted (or amplified) signals in the data and effectively ensures high-precision recovery of the gravity data amplitude.

[0125] (3) Adding compact constraints in the application can reduce the number of iterations of the inversion calculation.

[0126] To verify the effectiveness of the technology, we designed a two-dimensional theoretical model, in which the ground is set to be horizontal, the altitude is 0m, the observation range is -2000m to 2000m, and the point spacing is 20m; the field source is a square prism with a cross-sectional area of ​​200m×200m and a residual density of 3g / cm 3, the center coordinates (0, -80) in units of m; Figure 5 Shown are the forward gravity anomalies of the model at the 0m and 800m altitude planes.

[0127] An equivalent source model was set up below the ground. The equivalent source model consisted of a set of adjacent 20m x 20m cross-sections with a top depth of 100m. The gravity anomaly at an altitude of 800m was extended downward to the ground using both the traditional equivalent source method and the technology of the present invention. The improvement effect of the present invention on the downward continuation based on the equivalent source was demonstrated by comparing the physical properties of the equivalent source and the difference between the extended results and the theoretical values.

[0128] like Figure 6 As shown, Figure 6 (a) is the downward extension result of each method; Figure 6 (b) is the difference between the theoretical value and the extension result; Figure 6 (c) Obtain equivalent source property distributions for each method;

[0129] Figure 6 (a) shows the calculation results of the traditional equivalent source method and the technology of the present invention extending the 800m height data, as well as the difference with the theoretical gravity anomaly data on the ground. Figure 6 (b) is the difference between the theoretical value and each downward extension result. The figure shows that the traditional method's calculation results suffer from signal shape distortion and severe amplitude loss. After introducing the compact constraint, the signal distortion is significantly improved, and the signal amplitude increases accordingly. Building on the compact constraint and further incorporating a smoothing strategy, not only does the amplitude of the calculated result return to the theoretical data level, but the signal shape distortion is also effectively eliminated.

[0130] Figure 6 (c) shows the equivalent source density corresponding to each calculation result. The equivalent source density distribution obtained by the traditional equivalent source method has a main effective signal in the area corresponding to the source position, and there are also "divergent" signals with large amplitudes on both sides of it, which seriously affect the shape and amplitude of the downward continuation result. The present invention introduces compact constraints to suppress the "divergent" signal, but there is still a residual. Further smoothing measures are taken to reduce the amplitude to less than 4g / cm 3 The signal is equal to 0. After the third inversion calculation, the "divergent" signal is well eliminated, thereby significantly improving the accuracy of downward continuation, which illustrates the feasibility and effectiveness of the technology of the present invention.

[0131] Add 3% noise to the 800m height data ( Figure 7 ), the noisy data is extended downward, and the noisy data at 800m altitude can be extended downward to obtain the comparison results, such as Figure 8 As shown, Figure 8 (a) is the downward extension result of each method; Figure 8 (b) is the difference between the theoretical value and the extension result; Figure 8 (c) Comparison of the equivalent source property distribution obtained by each method. Figure 8 and Figure 6 This can better reflect the improvement effect of the technology of the present invention on the traditional equivalent source method.

[0132] Example 2

[0133] This embodiment provides a gravity data downward continuation system based on compact constraints, including:

[0134] The data acquisition and processing module M1 is used to acquire the gravity data of each measuring point in the study area and pre-process the gravity data.

[0135] The equivalent source geometric parameter determination module M2 is used to determine the geometric parameters of the equivalent source in the study area.

[0136] The equivalent source density calculation module M3 is used to obtain the density distribution of the equivalent source based on the preprocessed gravity data in combination with the compact constraint term and the smoothing strategy; the compact constraint term is determined according to the equivalent source density and the compact factor.

[0137] The equivalent source density calculation module M3 specifically includes:

[0138] The objective function establishing unit M31 is used to establish an inversion objective function according to the preprocessed gravity data, the equivalent source forward modeling kernel matrix and the compact constraint term.

[0139] The expression of the inversion objective function is:

[0140]

[0141] Where d represents the preprocessed gravity data; G represents the equivalent source forward kernel matrix; m represents the density of the equivalent source to be determined; μ represents the regularization parameter; m T W T Wm represents the compact constraint term; w j Represents the diagonal elements of the diagonal matrix W in the constraint terms; m′ j is the density of the jth equivalent source unit, α is the compact factor, and ε is 10 -7 .

[0142] The inversion calculation equation construction unit M32 is used to obtain the inversion calculation equation according to the inversion objective function.

[0143] The inversion calculation equation is expressed as:

[0144] P(G T G+μW T W)m=PGT d

[0145] Where P is the preconditioning matrix, P(G T G+μW T W)≈I, the preconditioning matrix takes the form of a diagonal matrix, and the diagonal elements of the preconditioning matrix G ij is the element of the equivalent source forward modeling kernel matrix G; i represents the measurement point code; M represents the number of measurement points.

[0146] The inversion and smoothing processing unit M33 is used to perform inversion calculation according to the inversion calculation equation in combination with the preconditioned conjugate gradient method to obtain the equivalent source density, and smooth the equivalent source density obtained by inversion to obtain the smoothed equivalent source density.

[0147] The equivalent source model construction module M4 is used to obtain the equivalent source model of the study area according to the geometric parameters of the equivalent source and the density distribution of the equivalent source.

[0148] The forward calculation module M5 is used to perform forward calculation based on the equivalent source model to obtain the downward continuation result of the gravity data.

[0149] Each embodiment in this specification focuses on the differences from other embodiments, and the same or similar parts between the embodiments can be referred to each other. For the system disclosed in the embodiment, since it corresponds to the method disclosed in the embodiment, the description is relatively simple, and the relevant parts can be referred to the method part.

[0150] This document uses specific examples to illustrate the principles and implementation methods of the present invention. The above examples are only intended to help understand the method and core concept of the present invention. At the same time, those skilled in the art will find that the specific implementation methods and application scopes may vary based on the concept of the present invention. In summary, the contents of this specification should not be construed as limiting the present invention.

Claims

1. A method for downward extension of gravity data based on compact constraints, characterized in that: include: Obtaining gravity data of each measuring point in the study area and preprocessing the gravity data; Determine the geometric parameters of the equivalent source in the study area; The density distribution of the equivalent source is obtained based on the preprocessed gravity data in combination with a compact constraint term and a smoothing strategy; the compact constraint term is determined according to the equivalent source density and the compact factor; Obtaining an equivalent source model of the study area according to the geometric parameters of the equivalent source and the density distribution of the equivalent source; A forward calculation is performed based on the equivalent source model to obtain a downward extension result of the gravity data.

2. The method according to claim 1, characterized in that The density distribution of the equivalent source is obtained based on the preprocessed gravity data in combination with the compact constraint term and the smoothing strategy, specifically including: establishing an inversion objective function based on the preprocessed gravity data, the equivalent source forward modeling kernel matrix, and the compact constraint term; Obtaining an inversion calculation equation according to the inversion objective function; The equivalent source density is obtained by performing inversion calculation according to the inversion calculation equation in combination with the preconditioned conjugate gradient method, and the obtained equivalent source density is smoothed to obtain the smoothed equivalent source density.

3. The method according to claim 2, characterized in that The expression of the inversion objective function is: Where d represents the preprocessed gravity data; G represents the equivalent source forward kernel matrix; m represents the density of the equivalent source to be determined; μ represents the regularization parameter; m T W T Wm represents the compact constraint term; w j Represents the diagonal elements of the diagonal matrix W in the constraint terms; m′ j is the density of the jth equivalent source unit, α is the compact factor, and ε is 10 -7 .

4. The method according to claim 3, characterized in that The inversion calculation equation is expressed as: P(G T G+μW T W)m=PG T d Where P is the preconditioning matrix, P(G T G+μW T W)≈I, the preconditioning matrix takes the form of a diagonal matrix, and the diagonal elements of the preconditioning matrix G ij is the element of the equivalent source forward modeling kernel matrix G; i represents the measurement point code; M represents the number of measurement points.

5. The method according to claim 4, characterized in that The inversion calculation is performed according to the inversion calculation equation in combination with the preconditioned conjugate gradient method to obtain the equivalent source density, and the obtained equivalent source density is smoothed to obtain the smoothed equivalent source density, specifically including: Let (G T G + μW T W) = A, G T d = B; Set the initial index of the equivalent source density to m0 = 0, the preset error value to e, r0 = A T B; n = 0, 1, ...; Let y n = Pr n ; when n = 0, S n = y n ; when n ≠ 0, S n = y n + β n-1 S n-1 ; Determine the iteration step size And according to the formula m n =m n-1 +t n S n Update equivalent source density; according to Determine whether the iteration is finished; When not satisfied When n=n+1, calculate r n+1 =r n -t n A T AS n , and return to step "Let y n =Pr n ; When n = 0, S n =y n ; When n≠0, S n =y n +β n-1 S n-1 ”; When satisfied When , the iteration ends, and the equivalent source density m after the nth iteration is obtained. n ; For the equivalent source density m n Smoothing is performed to obtain the smoothed equivalent source density.

6. The method according to claim 5, characterized in that The equivalent source density m n Performing a smoothing process to obtain the smoothed equivalent source density specifically includes: According to the equivalent source density m n Performing divergent signal detection and smoothing the detected divergent signal; The smoothed equivalent source density is used as the initial value m0 of the equivalent source density, and n is initialized to 0, and then return to step "Let y n =Pr n ; When n = 0, S n =y n ; When n≠0, S n =y n +β n-1 S n-1 "; until the amplitude of the divergent signal in the smoothed equivalent source density is less than the preset amplitude and (d-Gm) T The value of (d-Gm) is less than the preset value, and the final result of the equivalent source density m is obtained. inv .

7. A system based on the method according to any one of claims 1 to 6, characterized in that: include: A data acquisition and processing module is used to acquire gravity data of each measuring point in the study area and pre-process the gravity data; Equivalent source geometric parameter determination module, used to determine the geometric parameters of the equivalent source in the study area; An equivalent source density calculation module is used to obtain the density distribution of the equivalent source based on the preprocessed gravity data in combination with a compact constraint term and a smoothing strategy; the compact constraint term is determined according to the equivalent source density and the compact factor; An equivalent source model construction module, configured to obtain an equivalent source model of the study area according to geometric parameters of the equivalent source and density distribution of the equivalent source; A forward calculation module is used to perform forward calculation based on the equivalent source model to obtain a downward continuation result of the gravity data.

8. The system according to claim 7, characterized in that The equivalent source density calculation module specifically includes: An objective function establishing unit is used to establish an inversion objective function according to the preprocessed gravity data, the equivalent source forward modeling kernel matrix and the compact constraint term; An inversion calculation equation construction unit, configured to obtain an inversion calculation equation according to the inversion objective function; The inversion and smoothing processing unit is used to perform inversion calculation according to the inversion calculation equation in combination with the preconditioned conjugate gradient method to obtain the equivalent source density, and smooth the obtained equivalent source density to obtain the smoothed equivalent source density.

9. The system according to claim 8, characterized in that The expression of the inversion objective function is: Where d represents the preprocessed gravity data; G represents the equivalent source forward kernel matrix; m represents the density of the equivalent source to be determined; μ represents the regularization parameter; m T W T Wm represents the compact constraint term; w j Represents the diagonal elements of the diagonal matrix W in the constraint terms; m′ j is the density of the jth equivalent source unit, α is the compact factor, and ε is 10 -7 .

10. The system according to claim 9, characterized in that The inversion calculation equation is expressed as: P(G T G+μW T W)m=PG T d Where P is the preconditioning matrix, P(G T G+μW T W)≈I, the preconditioning matrix takes the form of a diagonal matrix, and the diagonal elements of the preconditioning matrix G ij is the element of the equivalent source forward modeling kernel matrix G; i represents the measurement point code; M represents the number of measurement points.

Images