CN116184519B - A method and system for downward extension of gravity data based on compact constraints - Google Patents
A method and system for downward extension of gravity data based on compact constraintsInfo
- Publication number
- CN116184519B CN116184519B CN202211414981.5A CN202211414981A CN116184519B CN 116184519 B CN116184519 B CN 116184519B CN 202211414981 A CN202211414981 A CN 202211414981A CN 116184519 B CN116184519 B CN 116184519B
- Authority
- CN
- China
- Prior art keywords
- equivalent source
- density
- compact
- inversion
- equivalent
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V7/00—Measuring gravitational fields or waves; Gravimetric prospecting or detecting
- G01V7/02—Details
- G01V7/06—Analysis or interpretation of gravimetric records
Landscapes
- Physics & Mathematics (AREA)
- Life Sciences & Earth Sciences (AREA)
- General Life Sciences & Earth Sciences (AREA)
- General Physics & Mathematics (AREA)
- Geophysics (AREA)
- Complex Calculations (AREA)
- Management, Administration, Business Operations System, And Electronic Commerce (AREA)
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
The invention relates to the field of gravity mechanics and gravity exploration, in particular to a gravity data downward continuation method and system based on compact constraint.
Background
Downward continuation is one of important contents in gravity data processing, and the observation data is converted into a passive space below the observation surface from the observation surface, so that the resolution of the observation surface to a field source can be improved, weak signals can be enhanced, and the requirement on gravity field information of different heights can be met. Such as extending down aerial gravity data to the ground, extending down shipboard gravity data to the seafloor thousands of meters deep, or extending down satellite gravity data at heights of hundreds of kilometers to near-surface heights. The high-precision and high-resolution gravity data obtained by downward continuation plays an important role in fields such as geology, ocean science, environmental science, resource exploration, national defense safety and the like.
The downward continuation of stability, accuracy and depth has been a difficult and hot spot problem in gravity data processing, where stability is a major factor limiting the maximum distance and accuracy of the downward continuation. The downward continuation is divided into two types of methods, namely a frequency domain and a spatial domain, such as a method based on Taylor series expansion, a method based on Fourier transformation, a analytic continuation method, an integral-iteration method, a solution method based on a bit field median theory and the like. The above method causes distortion of data morphology, even "drowning" the effective signal, and instability of the down-delay process, because of high frequency signal amplified, especially noise (or false interference) signal, either by solving the high order derivative or by amplifying the extension factor. To improve stability, the usual measures taken are regularization, filters, increasing the number of iterations, and multi-source data constraints. The regularization thought has positive effect on suppressing noise signals, but also sacrifices the resolution and amplitude of the signals, and the filter has similar effect to the regularization. Increasing the number of iterations may improve over-amplification of the high frequency signal, but still fail to fundamentally address the problem of noise. The mutual constraint among multi-source data is a popular practice at present, but in many cases, the problem that high-quality multi-source data cannot be obtained is faced. Therefore, although the former research work actively promotes the development of the downward extension method, the instability problem of the downward extension still severely limits the extension distance and the precision, and the intensive research work needs to be continuously carried out.
The equivalent source method is a relatively common technology in gravity data processing, and the method utilizes a group of artificially constructed virtual field sources (equivalent sources) to replace real field sources to generate gravity data in a passive space, namely an equivalent source model is built by fitting observation data, and the gravity data processing is converted into equivalent source model forward calculation. The method has the advantages of directly processing the original point observation data and considering the homology among the data, and has good performance in gravity data processing, but the method can also generate the problems of data form distortion and amplitude loss when used for carrying out large-depth downward continuation, and the root cause is a divergent signal existing in equivalent source physical properties. Therefore, obtaining the high-precision and high-resolution extension results puts higher demands on the equivalent source model, and a new technical method needs to be explored to ensure the stability of the large-depth downward extension.
Disclosure of Invention
The invention aims to provide a gravity data downward continuation method and system based on compact constraint, which lead equivalent source physical properties to be intensively distributed by introducing compact constraint items, and simultaneously consider a smoothing strategy to suppress or eliminate divergent signals in the equivalent source physical properties, thereby remarkably improving the problem of data morphological distortion, improving the recovery precision of data amplitude and ensuring the stability of large-depth downward continuation.
In order to achieve the above object, the present invention provides the following solutions:
A gravity data downward continuation method based on compact constraints, comprising:
acquiring gravity data of each measuring point in a research area, and preprocessing the gravity data;
determining geometric parameters of equivalent sources of the research area;
obtaining density distribution of an equivalent source based on the preprocessed gravity data in combination with a compact constraint term and a smoothing strategy, wherein the compact constraint term is determined according to the density of an equivalent source unit and a compact factor;
obtaining an equivalent source model of the research area according to the geometric parameters of the equivalent source and the density distribution of the equivalent source;
and forward calculation is carried out based on the equivalent source model to obtain a downward continuation result of the gravity data.
Optionally, the obtaining the density distribution of the equivalent source based on the preprocessed gravity data and the compact constraint term and the smoothing strategy specifically includes:
Establishing an inversion objective function according to the preprocessed gravity data, the equivalent source forward kernel matrix and the compact constraint term;
Obtaining an inversion calculation equation according to the inversion objective function;
And carrying out inversion calculation according to the inversion calculation equation and combining a precondition conjugate gradient method to obtain equivalent source density, and carrying out smoothing treatment on the obtained equivalent source density to obtain smoothed equivalent source density.
Optionally, the expression of the inversion objective function is:
Wherein d represents preprocessed gravity data, G represents an equivalent source forward kernel matrix, m represents the density of an equivalent source to be solved, mu represents a regularization parameter, and m TWT Wm represents a compact constraint term; w j represents the diagonal element of the diagonal matrix W in the constraint term, m' j is the j-th equivalent source unit density, α is a compact factor, and ε has a value of 10 -7.
Optionally, the expression of the inversion calculation equation is:
P(GTG+μWTW)m＝PGTd
wherein P is a precondition matrix, P (G TG+μWT W) is approximately equal to I, the precondition matrix adopts a diagonal matrix form, and diagonal elements of the precondition matrix G ij is an element of an equivalent source forward kernel matrix G, i represents the code of the measuring points, and M represents the number of the measuring points.
Optionally, the inversion calculation according to the inversion calculation equation of the equivalent source density combines with a precondition conjugate gradient method to perform inversion calculation to obtain the equivalent source density, and performs smoothing treatment on the obtained equivalent source density to obtain the smoothed equivalent source density, which specifically includes:
Let (G TG+μWTW)＝A,GT d=b;
Setting the initial index of equivalent source density to m 0 =0, the preset error value to e, r 0＝AT B, n=0, 1.;
Let y n＝Prn, when n= 0,S n＝yn, when n+. 0,S n＝yn+βn-1Sn-1;
Determining iteration step And updating the equivalent source density according to the formula m n＝mn-1+tnSn;
According to Judging whether iteration is finished;
when not meeting Let n=n+1 and calculate r n+1＝rn-tnATASn and return to step "let y n＝Prn; when n= 0,S n＝yn; when n+. 0,S n＝yn+βn-1Sn-1";
When meeting the requirements When the iteration is finished, obtaining the equivalent source density m n of inversion after the nth iteration;
and carrying out smoothing treatment on the equivalent source density m n to obtain the smoothed equivalent source density.
Optionally, the smoothing processing is performed on the equivalent source density m n to obtain the smoothed equivalent source density, which specifically includes:
Detecting a divergent signal according to the equivalent source density m n, and smoothing the detected divergent signal;
and taking the smoothed equivalent source density as an initial value m 0 of the equivalent source density, initializing n=0, returning to the step of enabling y n＝Prn, when n= 0,S n＝yn, when n is equal to 0,S n＝yn+βn-1Sn-1', obtaining an equivalent source density final result m inv until no divergence signal exists in the current smoothed equivalent source density and the value of (d-Gm) T (d-Gm) is smaller than a preset value.
The invention also provides a gravity data downward continuation system based on compact constraint, comprising:
The data acquisition and processing module is used for acquiring the gravity data of each measuring point in the research area and preprocessing the gravity data;
The equivalent source geometric parameter determining module is used for determining the geometric parameters of the equivalent source of the research area;
the equivalent source density calculation module is used for obtaining the density distribution of an equivalent source based on the preprocessed gravity data and combining a compact constraint item and a smoothing strategy, wherein the compact constraint item is determined according to the equivalent source density and a compact factor;
The equivalent source model construction module is used for obtaining an equivalent source model of the research area according to the geometric parameters of the equivalent source and the density distribution of the equivalent source;
And the forward calculation module is used for carrying out forward calculation based on the equivalent source model to obtain a downward continuation result of the gravity data.
Optionally, the equivalent source density calculation module specifically includes:
the objective function building unit is used for building an inversion objective function according to the preprocessed gravity data, the equivalent source forward kernel matrix and the compact constraint item;
the inversion calculation equation construction unit is used for obtaining an inversion calculation equation according to the inversion objective function;
and the inversion and smoothing processing unit is used for carrying out inversion calculation according to the inversion calculation equation and combining a precondition conjugate gradient method to obtain equivalent source density, and carrying out smoothing processing on the equivalent source density to obtain the smoothed equivalent source density.
Optionally, the expression of the objective function of the gravity equivalent source density is:
Wherein d represents preprocessed gravity data, G represents an equivalent source forward kernel matrix, m represents the density of an equivalent source to be solved, mu represents a regularization parameter, and m TWT Wm represents a compact constraint term; w j represents the diagonal element of the diagonal matrix W in the constraint term, m' j is the j-th equivalent source unit density, α is a compact factor, and ε has a value of 10 -7.
Optionally, the expression of the inversion calculation equation of the gravity equivalent source density is:
P(GTG+μWTW)m＝PGTd
wherein P is a precondition matrix, P (G TG+μWT W) is approximately equal to I, the precondition matrix adopts a diagonal matrix form, and diagonal elements of the precondition matrix G ij is an element of an equivalent source forward kernel matrix G, i represents the measurement point code, and M represents the number of measurement points.
According to the specific embodiment provided by the invention, the invention discloses the following technical effects:
The invention relates to a gravity data downward continuation method and system based on compact constraint, and the method comprises the steps of obtaining gravity data of each measuring point of a research area, preprocessing, determining geometric parameters of an equivalent source of the research area, obtaining density distribution of the equivalent source based on the preprocessed gravity data in combination with compact constraint terms and a smoothing strategy, determining the compact constraint terms according to the density of the equivalent source and compact factors, obtaining an equivalent source model according to the geometric parameters of the equivalent source and the density of the equivalent source, and performing forward calculation based on the equivalent source model to obtain a downward continuation result of the gravity data. The compact constraint item is introduced to enable the equivalent source physical properties to be intensively distributed, meanwhile, a smoothing strategy is considered to suppress or eliminate divergent signals in the equivalent source density, so that the distortion of the data form is remarkably improved, the recovery precision of the data amplitude is improved, and the stability of large-depth downward continuation is ensured.
Drawings
In order to more clearly illustrate the embodiments of the present invention or the technical solutions in the prior art, the drawings that are needed in the embodiments will be briefly described below, and it is obvious that the drawings in the following description are only some embodiments of the present invention, and other drawings may be obtained according to these drawings without inventive effort for a person skilled in the art.
FIG. 1 is a flow chart of a gravity data downward continuation method based on compact constraints provided in embodiment 1 of the present invention;
FIG. 2 is a technical schematic diagram of a gravity data downward continuation method based on compact constraint according to embodiment 1 of the present invention;
FIG. 3 is a schematic diagram of an equivalent source provided in embodiment 1 of the present invention;
FIG. 4 is a graph showing the variation of the number of iterations with μ according to embodiment 1 of the present invention;
FIG. 5 is a theoretical model and forward gravity anomaly data provided in embodiment 1 of the present invention;
FIG. 6 is a graph showing the comparison of the downward continuation result of the 800m height data provided in example 1 of the present invention;
FIG. 7 is a graph of 800m high noise and raw gravity anomaly data provided in example 1 of the present invention;
FIG. 8 is a graph showing the comparison of the downward continuation results of 800m high noise data provided in example 1 of the present invention.
Detailed Description
The following description of the embodiments of the present invention will be made clearly and completely with reference to the accompanying drawings, in which it is apparent that the embodiments described are only some embodiments of the present invention, but not all embodiments. All other embodiments, which can be made by those skilled in the art based on the embodiments of the invention without making any inventive effort, are intended to be within the scope of the invention.
The invention aims to provide a gravity data downward continuation method and system based on compact constraint, which lead equivalent source physical properties to be intensively distributed by introducing compact constraint items, and simultaneously consider a smoothing strategy to suppress or eliminate divergent signals in the equivalent source physical properties, thereby remarkably improving the problem of data morphological distortion, improving the recovery precision of data amplitude and ensuring the stability of large-depth downward continuation.
In order that the above-recited objects, features and advantages of the present invention will become more readily apparent, a more particular description of the invention will be rendered by reference to the appended drawings and appended detailed description.
Example 1
As shown in fig. 1 and 2, the present embodiment provides a gravity data downward continuation method based on compact constraint, including:
S1, acquiring gravity data of each measuring point in a research area, and preprocessing the gravity data.
The embodiment works against the Bragg gravity anomaly data or the Bragg gravity anomaly gradient data, wherein the data information comprises east, north and vertical coordinates, data values and total data accuracy of a Cartesian coordinate system.
The data were screened (pre-processed) as follows:
(1) And (3) mapping the Bragg gravity anomaly data or the Bragg gravity anomaly gradient data on the discrete measuring points in the measuring line or the measuring region (contour map or image map), and eliminating the jump points or sharp points in the data according to the change trend and the basic rule of the data in the map.
(2) According to the requirement of the data resolution (pre-determined by task requirements), resampling the data after eliminating the jump points or the sharp points, removing redundant data, namely, dividing a research area into grids according to the resolution, randomly reserving 1 measuring point in each grid through which a measuring line passes, taking the resampled data as the data finally participating in calculation (namely, the preprocessed gravity data participating in calculation in the step S3), and resampling the data to reduce the calculated amount, avoid overfitting and provide reference for the horizontal size of an equivalent source unit.
S2, determining the geometric parameters of the equivalent source of the research area.
(1) Equivalent sources are equivalent to real field sources. The real field source is all underground substances, strata, rocks, minerals and the like, has mass and density, and is equivalent to the earth or a part of the earth. The equivalent source consists of a set of discrete, uniformly dense cuboids arranged in close proximity, with the spacing (or average spacing) of the reference data points determining the geometry of the cuboid.
(2) Depending on the intensity of the terrain relief, the equivalent source may optionally be positioned below the surface over a curved or planar surface within a depth range.
The method comprises the steps of (1) obtaining a horizontal dimension of an equivalent source unit according to the gravity data resolution requirement of the step (2), obtaining a value of the minimum dimension of the unit when the prior information is available, obtaining the equivalent source unit according to the field source average depth reflected by the prior information when the prior information is available, obtaining the equivalent source unit horizontal dimension of the unit according to the gravity data resolution requirement of the step (2), obtaining the equivalent source unit thickness of 3-6 times of the minimum dimension of the unit when the relief of the topography is severe, obtaining the equivalent source unit on a curved surface parallel to the topography when the relief of the topography is severe, and obtaining the equivalent source unit outside the data coverage area by expanding the edge by 3-5 equivalent source units when the equivalent source layer is located on the plane to suppress boundary effects.
And S3, obtaining the density distribution of the equivalent source based on the preprocessed gravity data and combining a compact constraint term and a smoothing strategy, wherein the compact constraint term is determined according to the density of the equivalent source and a compact factor.
The step S3 specifically includes:
and S31, establishing an inversion objective function according to the preprocessed gravity data, the equivalent source forward kernel matrix and the compact constraint term.
In particular, the method comprises the steps of,
Establishing an equivalent source forward kernel matrix G= { G ij } according to the equivalent source geometric parameter and the space position coordinates of the preprocessed data d, wherein a forward calculation formula adopts a published cuboid forward calculation analytic formula, and an inversion objective function is established:
φ=φd+μφm (1)
Wherein phi d is the variance (d-Gm) T (d-Gm) between the observed data and the model calculated value, m is the equivalent source density vector to be solved, the second term phi m is a compact constraint term, the distribution of the equivalent source density can be influenced by introducing constraint information through the term, and the distribution is expressed as a matrix:
φm＝mTWTWm (2)
Where W is a diagonal matrix, the diagonal elements are:
Wherein m' j is the j-th equivalent source unit density, belongs to an element in the vector m of the equivalent source density to be solved, alpha is a compact factor, epsilon is an extremely small number, and the value is 10 -7. The equivalent source units W j determined by the formula (3) form a diagonal matrix W, and are all brought into the target function formula (1) to obtain
The expression of the inversion objective function is:
Wherein d represents preprocessed gravity data, G represents an equivalent source forward kernel matrix, m represents equivalent source density to be solved, mu represents regularization parameters, and m TWT Wm represents a compact constraint term phi m; W j denotes the diagonal element of the diagonal matrix W, m' j is the j-th equivalent source unit density, α is a compact factor, and ε has a value of 10 -7.
S32, obtaining an inversion calculation equation according to the inversion objective function.
Taking the derivative of m T for phi and letting it be zero, i.e. dphi/dm T =0, we can obtain
(GTG+μWTW)m＝GTd (5)
Multiplying both ends of the formula (5) by a precondition matrix P to obtain an equation finally used for inversion calculation:
That is, the expression of the inversion calculation equation for obtaining the equivalent source density is:
P(GTG+μWTW)m＝PGTd (6)
Wherein P is a precondition matrix for improving the condition number of the coefficient matrix, and P (G TG+μWT W) is about I, the precondition matrix adopts a diagonal matrix form, and diagonal elements of the precondition matrix G ij is an element of an equivalent source forward kernel matrix G, i represents the code of the measuring points, and M represents the number of the measuring points.
And S33, carrying out inversion calculation according to the inversion calculation equation and combining a precondition conjugate gradient method to obtain equivalent source density, and carrying out smoothing treatment on the equivalent source density to obtain the smoothed equivalent source density.
Wherein S33 specifically includes:
Let (G TG+μWTW)＝A,GT d=b) in formula (6).
S332, setting the initial value of the equivalent source density to m 0 =0, and setting the preset error value to e, r 0＝AT B. The initial value of n is 0.
S333 let y n＝Prn, when n= 0,S n＝yn, when n+. 0,S n＝yn+βn-1Sn-1, the step corresponds to the updated search direction in the precondition conjugate gradient method
S334, determining an iteration stepAnd updates the equivalent source density according to equation m n＝mn-1+tnSn.
S335 according toAnd judging whether the iteration is ended or not.
When not meetingIf so, then the next inversion iteration is performed, let n=n+1, and r n+1＝rn-tnATASn is calculated, and step S333 "let y n＝Prn" is returned, when n= 0,S n＝yn, when n+. 0,S n＝yn+βn-1Sn-1 ".
When meeting the requirementsAnd when the inversion iteration is finished, obtaining the inversion equivalent source density m n after the nth iteration.
And S336, performing smoothing treatment on the equivalent source density m n to obtain the smoothed equivalent source density.
Specifically, step S336 includes:
(1) And detecting a divergence signal according to the equivalent source density m n, and smoothing the detected divergence signal. After the compact constraint, the effective signal and the divergent signal have macroscopic distinction, the divergent signal is the signal with the wave shape characteristic around the effective signal, the amplitude of the divergent signal is smaller than that of the effective signal, and the two signals have larger difference, and the two signals can be distinguished by combining the shape and the amplitude of the signals (the signals with the amplitude smaller than a certain value can be regarded as divergent signals).
And checking the equivalent source density distribution characteristics obtained by inversion, and determining the distribution range of the residual divergent signal.
Depending on the determined "divergent" signal distribution range, the signal within that distribution range is zeroed out or equal to some minimum (i.e., smoothed).
(2) And taking the smoothed equivalent source density of the divergent signal as an initial value m 0 of the equivalent source density, initializing n=0, returning to step S333 to enable y n＝Prn, when n= 0,S n＝yn, when n is not equal to 0,S n＝yn+βn-1Sn-1, until the current amplitude of the divergent signal is smaller than a preset amplitude and the value of (d-Gm) T (d-Gm) is smaller than a preset value, and obtaining an equivalent source density final result m inv, namely the smoothed equivalent source density.
The step of assigning the smoothed equivalent source density as an initial value of the equivalent source density, returning to the inversion iteration process of the step S33, inverting again to an equivalent source density from 0 for the iteration times, performing smoothing again, assigning the smoothed equivalent source density as the initial value of the equivalent source density again, performing the inversion iteration process of the step S333 again until the amplitude of the divergence signal in the equivalent source density obtained by inversion is smaller than a preset amplitude and the value of (d-Gm) T (d-Gm) is smaller than a preset value, and outputting the final equivalent source density.
The "divergent" signal in the equivalent source property can be significantly attenuated by the compact constraint, but there is still a residual; the method comprises the steps of carrying out compact constraint, carrying out inversion calculation by taking the physical property of an equivalent source after smoothing as an initial value of the density of the equivalent source, and iterating for a plurality of times until the divergent signal is effectively pressed.
In order to compromise the accuracy and efficiency of the inversion calculation, the values of the regularization parameter μ and the compaction factor α need to be reasonably determined.
The regularization parameter mu has the functions of balancing the influence degree of the first term (data fitting function) and the second term (constraint term) on equivalent source physical properties in the inversion target function formula (1), on the premise that the data fitting difference is determined, when mu is smaller, the effect of the constraint term is biased on the fitting of data, and when mu is larger, the effect of constraint information is well exerted, but the iteration times are increased, so that the calculation efficiency is reduced, and therefore, the embodiment takes mu near the jump point of the iteration times in consideration of the constraint effect and the calculation efficiency.
The general rule of the compact factor alpha selection is that when the value of the compact factor alpha is smaller, the inversion iteration times are more, the equivalent source physical property distribution is smoother, when the value of the compact factor alpha is larger, the iteration times are less, but the equivalent source physical property signal has larger fluctuation, so that a larger alpha value is selected under the condition of prior information, iteration is started from smaller alpha under the condition of no prior information, and then the alpha value is gradually increased according to a certain rule until the alpha value is increased to a certain value and then is kept unchanged. In this embodiment, the method of giving a smaller initial value (0-1) to α and then gradually increasing the α value at each iteration, for example, increasing by 0.2 each time, until increasing by 4, is kept unchanged.
Therefore, when determining the optimal value of μ, a plurality of random values of μmay be predetermined, and these random values are empirically determined, and are respectively brought into equation (6) and respectively participate in inversion calculation, that is, each μ is brought into equation (6), step S33 is performed, and a schematic diagram of the number of iterations with μ is obtained, and as shown in fig. 4, the optimal value of μ is determined according to the jump point of the number of iterations in the schematic diagram of the number of iterations with μ.
And S4, obtaining an equivalent source model of the research area according to the geometric parameters of the equivalent source and the density distribution of the equivalent source.
Equivalent source methods provide a set of simple artificial field sources (typically consisting of discrete cuboids or particles) instead of real field sources to generate gravitational field information in a passive space. The general method comprises the steps of firstly manually setting the space position and the size of an equivalent source, namely geometric parameters (namely the process of step S2), wherein the physical properties of the equivalent source are undetermined, then establishing an equation set with observed data, and inverting and calculating the physical property distribution of the equivalent source (namely the process of step S3) through fitting the observed data, namely establishing an equivalent source model, so that the conversion and prediction processing of gravity data can be converted into forward calculation of the equivalent source model.
And S5, forward calculation is carried out based on the equivalent source model to obtain a downward continuation result of the gravity data.
The forward formula applied in this step may be a uniform density cuboid forward formula G'm inv = d', where the spatial position coordinates of the data points in G 'and d' should be exchanged for coordinates of the down-extending positions. Wherein G 'is forward kernel function for downward continuation, and d' is gravity data of downward continuation (i.e. downward continuation result to be solved)
According to the basic principle of the potential field equivalent source, the embodiment optimizes the physical property distribution of the equivalent source by introducing a compact constraint term and combines a smoothing strategy to suppress or eliminate the divergent signal in the physical property of the equivalent source and simultaneously ensure the amplitude of the effective signal, thereby realizing the purpose of carrying out high-precision and large-depth downward continuation on the gravity data at a certain height in the passive space. The invention can obtain a reliable three-dimensional passive space gravity database, and establishes a gravity field model of the area, thereby providing technical support and basic data for applications such as gravity research, gravity navigation reference diagram, gravity exploration, gravity instrument test and the like. Therefore, the present embodiment has the following advantages:
(1) From the aspect of restraining the equivalent source physical property distribution, the method directly aims at the fundamental problem of causing instability of downward continuation, and compared with the traditional equivalent source method, the method has more remarkable and effective effect of improving the stability of downward continuation.
(2) The gravity data is extended downwards by combining the compact constraint and smoothing strategy, so that the suppression effect on distortion (or amplification) signals in the data is obviously improved, and the high-precision recovery of the gravity data amplitude is effectively ensured.
(3) Adding compact constraint terms to the application may reduce the number of iterations of the inversion calculation.
In order to verify the effectiveness of the technology, a two-dimensional theoretical model is designed, wherein the ground is set to be horizontal, the altitude is 0m, the observation range is-2000 m, the point distance is 20m, the field source is a square cylinder with the cross section of 200m multiplied by 200m, the residual density is 3g/cm 3, the center coordinate (0, -80) unit m, and the model is forward gravity anomaly on the height plane of 0m and 800m as shown in fig. 5.
An equivalent source model is arranged under the ground, the equivalent source model consists of a group of adjacent cross sections of 20m multiplied by 20m, and the depth of the top surface is 100m. The gravity anomaly at the height of 800m is extended downwards to the ground by using the traditional equivalent source method and the technology of the invention, and the improvement effect of the invention on the downward extension based on the equivalent source is illustrated by comparing the physical property of the equivalent source and the difference between the extension result and the theoretical value.
As shown in fig. 6, fig. 6 (a) shows the downward continuation result of each method, fig. 6 (b) shows the difference between the theoretical value and the continuation result, and fig. 6 (c) shows the equivalent source physical property distribution obtained by each method;
Fig. 6 (a) shows the calculation result of the conventional equivalent source method and the height data of 800m extended by the technology of the present invention, and the difference value from the theoretical gravity anomaly data on the ground, and fig. 6 (b) shows the difference value between the theoretical value and each downward extended result. The figure shows that the signal shape distortion exists in the calculation result of the traditional method, the amplitude loss is serious, the signal distortion condition is well improved after the compact constraint is introduced, the signal amplitude is correspondingly increased, the smoothing strategy is further combined on the basis of the compact constraint, the amplitude of the calculation result is restored to the level of theoretical data, and the signal shape distortion is well eliminated.
FIG. 6 (c) shows the equivalent source density corresponding to each calculation result, wherein the equivalent source density distribution obtained by the traditional equivalent source method has a main effective signal in the area corresponding to the field source position, and the 'divergent' signals with larger amplitude are also arranged at the two sides of the effective signal, so that the form and the amplitude of the downward continuation result are seriously influenced.
Adding 3% of noise into 800m height data (figure 7), carrying out downward extension on the noisy data to obtain the downward extension result comparison of the noisy data of 800m height, as shown in figure 8, wherein figure 8 (a) is the downward extension result of each method, figure 8 (b) is the difference value between the theoretical value and the extension result, figure 8 (c) is the equivalent source physical property distribution obtained by each method, and the improvement effect of the technology of the invention on the traditional equivalent source method can be more embodied compared with figures 8 and 6.
Example 2
The present embodiment provides a gravity data downward continuation system based on compact constraints, comprising:
The data acquisition and processing module M1 is used for acquiring the gravity data of each measuring point in the research area and preprocessing the gravity data.
And the equivalent source geometric parameter determination module M2 is used for determining the geometric parameters of the equivalent source of the research area.
The equivalent source density calculation module M3 is used for obtaining the density distribution of the equivalent source based on the preprocessed gravity data and combining a compact constraint term and a smoothing strategy, wherein the compact constraint term is determined according to the equivalent source density and a compact factor.
The equivalent source density calculating module M3 specifically includes:
And the objective function establishing unit M31 is used for establishing an inversion objective function according to the preprocessed gravity data, the equivalent source forward kernel matrix and the compact constraint term.
The expression of the inversion objective function is:
Wherein d represents preprocessed gravity data, G represents an equivalent source forward kernel matrix, m represents the density of an equivalent source to be solved, mu represents a regularization parameter, and m TWT Wm represents a compact constraint term; w j represents the diagonal element of the diagonal matrix W in the constraint term, m' j is the j-th equivalent source unit density, α is a compact factor, and ε has a value of 10 -7.
And the inversion calculation equation construction unit M32 is used for obtaining an inversion calculation equation according to the inversion objective function.
The expression of the inversion calculation equation is:
P(GTG+μWTW)m＝PGTd
wherein P is a precondition matrix, P (G TG+μWT W) is approximately equal to I, the precondition matrix adopts a diagonal matrix form, and diagonal elements of the precondition matrix G ij is an element of an equivalent source forward kernel matrix G, i represents the code of the measuring points, and M represents the number of the measuring points.
And the inversion and smoothing processing unit M33 is used for carrying out inversion calculation according to the inversion calculation equation and combining a precondition conjugate gradient method to obtain equivalent source density, and carrying out smoothing processing on the equivalent source density obtained by inversion to obtain the smoothed equivalent source density.
And the equivalent source model construction module M4 is used for obtaining an equivalent source model of the research area according to the geometric parameters of the equivalent source and the density distribution of the equivalent source.
And the forward computing module M5 is used for performing forward computation based on the equivalent source model to obtain a downward continuation result of the gravity data.
In this specification, each embodiment is mainly described in the specification as a difference from other embodiments, and the same similar parts between the embodiments are 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 points refer to the description of the method section.
The principles and embodiments of the present invention have been described herein with reference to specific examples, which are intended to facilitate an understanding of the principles and concepts of the invention and are to be varied in scope and detail by persons of ordinary skill in the art based on the teachings herein. In view of the foregoing, this description should not be construed as limiting the invention.
Claims (10)
1. A gravity data downward continuation method based on compact constraints, comprising:
acquiring gravity data of each measuring point in a research area, and preprocessing the gravity data;
determining geometric parameters of equivalent sources of the research area;
Obtaining density distribution of an equivalent source based on the preprocessed gravity data in combination with a compact constraint term and a smoothing strategy, wherein the compact constraint term is determined according to the equivalent source density and a compact factor;
obtaining an equivalent source model of the research area according to the geometric parameters of the equivalent source and the density distribution of the equivalent source;
and forward calculation is carried out based on the equivalent source model to obtain a downward continuation result of the gravity data.
2. The method according to claim 1, wherein the deriving the density distribution of the equivalent sources based on the preprocessed gravity data in combination with a compact constraint term and a smoothing strategy comprises:
Establishing an inversion objective function according to the preprocessed gravity data, the equivalent source forward kernel matrix and the compact constraint term;
Obtaining an inversion calculation equation according to the inversion objective function;
and carrying out inversion calculation according to the inversion calculation equation and combining a precondition conjugate gradient method to obtain equivalent source density, and carrying out smoothing treatment on the obtained equivalent source density to obtain smoothed equivalent source density.
3. The method of claim 2, wherein the inversion objective function is expressed as:
Wherein d represents preprocessed gravity data, G represents an equivalent source forward kernel matrix, m represents the density of an equivalent source to be solved, mu represents a regularization parameter, and m TWT Wm represents a compact constraint term; w j represents the diagonal element of the diagonal matrix W in the constraint term, m' j is the j-th equivalent source unit density, α is a compact factor, and ε has a value of 10 -7.
4. A method according to claim 3, wherein the expression of the inversion calculation equation is:
P(GTG+μWTW)m＝PGTd
wherein P is a precondition matrix, P (G TG+μWT W) is approximately equal to I, the precondition matrix adopts a diagonal matrix form, and diagonal elements of the precondition matrix G ij is an element of an equivalent source forward kernel matrix G, i represents the code of the measuring points, and M represents the number of the measuring points.
5. The method according to claim 4, wherein the inversion calculation is performed according to the inversion calculation equation in combination with a preconditioned conjugate gradient method to obtain an equivalent source density, and the obtained equivalent source density is smoothed to obtain a smoothed equivalent source density, and the method specifically comprises:
Let (G TG+μWTW)＝A,GT d=b;
Setting the initial index of equivalent source density to m 0 =0, the preset error value to e, r 0＝AT B, n=0, 1.;
Let y n＝Prn, when n= 0,S n＝yn, when n+. 0,S n＝yn+βn-1Sn-1;
Determining iteration step And updating the equivalent source density according to the formula m n＝mn-1+tnSn;
According to Judging whether iteration is finished;
when not meeting Let n=n+1 and calculate r n+1＝rn-tnATASn and return to step "let y n＝Prn; when n= 0,S n＝yn; when n+. 0,S n＝yn+βn-1Sn-1";
When meeting the requirements When the iteration is finished, obtaining the equivalent source density m n of inversion after the nth iteration;
and carrying out smoothing treatment on the equivalent source density m n to obtain the smoothed equivalent source density.
6. The method of claim 5, wherein the smoothing the equivalent source density m n to obtain the smoothed equivalent source density specifically includes:
Detecting a divergent signal according to the equivalent source density m n, and smoothing the detected divergent signal;
And taking the smoothed equivalent source density as an initial value m 0 of the equivalent source density, initializing n=0, returning to the step of ' y n＝Prn ', when n= 0,S n＝yn, when n is equal to 0,S n＝yn+βn-1Sn-1 ', and obtaining an equivalent source density final result m inv when the divergence signal amplitude in the smoothed equivalent source density is smaller than a preset amplitude and the value of (d-Gm) T (d-Gm) is smaller than a preset value.
7. A system based on the method of any one of claims 1 to 6, comprising:
The data acquisition and processing module is used for acquiring the gravity data of each measuring point in the research area and preprocessing the gravity data;
The equivalent source geometric parameter determining module is used for determining the geometric parameters of the equivalent source of the research area;
The equivalent source density calculation module is used for obtaining the density distribution of an equivalent source based on the preprocessed gravity data and combining a compact constraint item and a smoothing strategy, wherein the compact constraint item is determined according to the equivalent source density and a compact factor;
The equivalent source model construction module is used for obtaining an equivalent source model of the research area according to the geometric parameters of the equivalent source and the density distribution of the equivalent source;
And the forward calculation module is used for carrying out forward calculation based on the equivalent source model to obtain a downward continuation result of the gravity data.
8. The system of claim 7, wherein the equivalent source density calculation module specifically comprises:
the objective function building unit is used for building an inversion objective function according to the preprocessed gravity data, the equivalent source forward kernel matrix and the compact constraint item;
the inversion calculation equation construction unit is used for obtaining an inversion calculation equation according to the inversion objective function;
And the inversion and smoothing processing unit is used for carrying out inversion calculation according to the inversion calculation equation and combining a precondition conjugate gradient method to obtain equivalent source density, and carrying out smoothing processing on the obtained equivalent source density to obtain the smoothed equivalent source density.
9. The system of claim 8, wherein the expression of the inversion objective function is:
Wherein d represents preprocessed gravity data, G represents an equivalent source forward kernel matrix, m represents the density of an equivalent source to be solved, mu represents a regularization parameter, and m TWT Wm represents a compact constraint term; w j represents the diagonal element of the diagonal matrix W in the constraint term, m' j is the j-th equivalent source unit density, α is a compact factor, and ε has a value of 10 -7.
10. The system of claim 9, wherein the expression of the inversion calculation equation is:
P(GTG+μWTW)m＝PGTd
wherein P is a precondition matrix, P (G TG+μWT W) is approximately equal to I, the precondition matrix adopts a diagonal matrix form, and diagonal elements of the precondition matrix G ij is an element of an equivalent source forward kernel matrix G, i represents the code of the measuring points, and M represents the number of the measuring points.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202211414981.5A CN116184519B (en) | 2022-11-11 | 2022-11-11 | A method and system for downward extension of gravity data based on compact constraints |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202211414981.5A CN116184519B (en) | 2022-11-11 | 2022-11-11 | A method and system for downward extension of gravity data based on compact constraints |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN116184519A CN116184519A (en) | 2023-05-30 |
| CN116184519B true CN116184519B (en) | 2025-09-12 |
Family
ID=86431414
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202211414981.5A Active CN116184519B (en) | 2022-11-11 | 2022-11-11 | A method and system for downward extension of gravity data based on compact constraints |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN116184519B (en) |
Families Citing this family (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN119126249B (en) * | 2024-09-10 | 2025-03-25 | 中国地震局地球物理研究所 | A method for optimizing mineral resource exploration based on gravity forward and inversion |
Citations (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN112363236A (en) * | 2020-10-15 | 2021-02-12 | 中国地质大学（武汉） | Gravity field data equivalent source continuation and data type conversion method based on PDE |
| CN114154111A (en) * | 2021-11-26 | 2022-03-08 | 兰州大学 | Bit field data downward continuation method for frequency domain continuous-fraction expansion |
Family Cites Families (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US11487036B2 (en) * | 2017-01-12 | 2022-11-01 | Cgg Services Sas | Reflection full waveform inversion methods with density and velocity models updated separately |
-
2022
- 2022-11-11 CN CN202211414981.5A patent/CN116184519B/en active Active
Patent Citations (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN112363236A (en) * | 2020-10-15 | 2021-02-12 | 中国地质大学（武汉） | Gravity field data equivalent source continuation and data type conversion method based on PDE |
| CN114154111A (en) * | 2021-11-26 | 2022-03-08 | 兰州大学 | Bit field data downward continuation method for frequency domain continuous-fraction expansion |
Also Published As
| Publication number | Publication date |
|---|---|
| CN116184519A (en) | 2023-05-30 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN111323776B (en) | Method for monitoring deformation of mining area | |
| CN111665556B (en) | Construction Method of Formation Acoustic Velocity Model | |
| WO2018129844A1 (en) | Seismic diffracted wave separation method and device | |
| CN116466402B (en) | Electromagnetic inversion method based on geological information and electromagnetic data combined driving | |
| CN116774292B (en) | Seismic wave travel time determining method, system, electronic equipment and storage medium | |
| CN107783185B (en) | A kind of processing method and device for chromatographic static correction | |
| CN102866421A (en) | Scattered wave pre-stack imaging method for identifying small-fault throw breakpoints | |
| CN111856598B (en) | Magnetic measurement data multilayer equivalent source upper extension and lower extension method | |
| CN109884710A (en) | For the micro logging chromatography imaging method of excitation well depth design | |
| CN116184519B (en) | A method and system for downward extension of gravity data based on compact constraints | |
| CN120446906A (en) | A method and system for measuring and determining forest tree height based on lidar point cloud data | |
| CN113109793A (en) | Adaptive resolution water depth curved surface filtering method and device | |
| CN109655918A (en) | Ground shallow well micro-seismic monitoring measuring platform station location determines method and system | |
| CN114167511A (en) | Continuous-fraction expansion downward continuation-based bit field data rapid inversion method | |
| CN115032694B (en) | A VSP first arrival travel time tomography method and system | |
| CN114722590B (en) | Design optimization method of random acquisition observation system based on geophysical model | |
| CN114966878B (en) | A 3D Gravity Field Inversion Method Based on Mixed Norm and Cross-correlation Constraints | |
| CN115826038A (en) | VTI medium adjoint state method travel time multi-parameter tomography method and system | |
| CN115598700B (en) | A method, apparatus, storage medium, and electronic device for seismic profile imaging. | |
| CN121186774A (en) | A method and system for slope monitoring based on unmanned aerial vehicles (UAVs) | |
| Li et al. | Modified Bott-Parker method for gravimetric Moho modeling | |
| CN109444973A (en) | Gravity forward modeling accelerated method under a kind of spherical coordinate system | |
| CN116794735B (en) | Aeromagnetic vector gradient data equivalent source multi-component joint denoising method and device | |
| CN114740540B (en) | Method and system for constructing magnetic anomaly map of middle ridge area in ocean based on direction constraint | |
| CN113568033A (en) | Design method and device of three-dimensional irregular sampling seismic acquisition observation system |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |