High-precision density inversion method and system for stably and rapidly solving gravity anomaly based on regularization

Abstract

The invention discloses a high-precision density inversion method and system for stably and rapidly solving gravity anomalies based on regularization, on the basis of regularization solving of a field source, data carrying false anomalies at the imaging depth of the field source is used as a boundary value condition, and the first kind of boundary value problem of the Laplacian equation is solved in the upper half space, so that the high-precision density inversion of the gravity anomalies is realized. And calculating the gravity anomaly at the ground, and carrying out least square optimization iteration solution on the gravity anomaly and the actual ground anomaly to obtain high-precision anomaly imaging data at the field source depth without Gibbs effect false anomaly. And taking the depth as a top interface of the field source, dividing the grid prisms by the top interface to construct a field source form, and performing density inversion on the top surfaces of the prisms in a frequency domain to obtain the density distribution of the field source form. The algorithm is realized in a frequency domain, and inversion is carried out at a known field source depth, so that the method is a rapid, high-precision and high-resolution characteristic gravity density inversion method, and a new idea is provided for solving gravity field source density parameters.

Description

Technical Field

[0001] The present invention relates to the technical field of mineral exploration and detection methods, and in particular to a high-precision density inversion method and system for rapidly solving gravity anomalies based on regularization and stability. Background Art

[0002] As a commonly used geophysical exploration method, gravity exploration is widely used in energy and mineral resource exploration, geological mapping, engineering construction and other aspects. Gravity anomaly data processing and conversion is a key part of gravity data interpretation theory. Taking the ground gravity anomaly as the first-class boundary condition, the upper half space (i.e., far away from the field source) of the Laplace equation is solved to obtain a regular and stable solution. On the other hand, solving the lower half space (i.e., close to the field source) of the Laplace equation can obtain higher-resolution gravity anomaly information, but the lower half space will be restricted by the boundary value of the Laplace equation, and the gravity anomaly information at the field source cannot be obtained. In order to ensure that the potential field solved at the field source is not singular during the solution of the lower half space, and to obtain the complete field value distribution of the lower half space, it is necessary to optimize and improve the gravity field solution method of the lower half space.

[0003] The lower half-space distribution characteristics of the gravity and magnetic potential field are highly intuitive and are an important basis for studying the field source body. By taking the ground-observed gravity anomaly as the first-class boundary condition and performing a regularized stable solution of the Laplace lower half-space, high-resolution anomaly data at the periphery of the field source depth can be obtained. However, due to the influence of the Gibbs truncation effect of the frequency domain solution, there are often "eight-shaped" alternating positive and negative weak anomaly belts on both sides of the main anomaly belt, centered on the main source body distribution area. It may be a single anomaly, or it may be superimposed with the main anomaly of another source body to form a composite anomaly, which increases the difficulty of interpreting the field source body. Summary of the invention

[0004] In order to solve the above technical problems, the present invention provides a high-precision density inversion method and system based on regular stable rapid solution of gravity anomaly. On the basis of regularized solution of field source, it is proposed to use the data carrying false anomalies at the imaging depth of the field source as boundary conditions, and again calculate the gravity anomaly at the ground through the solution of the first-class boundary value problem of Laplace equation in the upper half space, and perform least square optimization iteration solution with the actual ground anomaly to obtain high-precision anomaly imaging data at the depth of the field source without Gibbs effect false anomaly, and use the depth as the top interface of the field source, and construct the field source form by dividing the grid prism from the top interface, and perform density inversion on the top surface of each prism in the frequency domain, thereby obtaining the density distribution of the field source form. Since the algorithm is implemented in the frequency domain and the inversion is performed at the known field source depth, the invention is a fast, high-precision, high-resolution characteristic gravity density inversion method, which provides a new idea for comprehensive interpretation of the field source body from both physical properties and geometric parameters.

[0005] The present invention discloses a high-precision density inversion method for quickly solving gravity anomalies based on regularization stability, the method comprising:

[0006] S1. Through the Laplace equation of gravity field and the corresponding boundary conditions, i.e. the gravity data observed on the ground, a stable Fourier series solution under general conditions is obtained;

[0007] S2. In the lower half space, a stable Fourier series with a regularization factor is selected to describe the gravity solution of the lower half space, and the center of the gravity anomaly corresponds to the center of the field source;

[0008] S3, taking the data containing false anomalies at the depth center of the field source in the lower half space as the starting boundary condition of the upper half space, calculating the ground gravity anomaly data based on the boundary condition by solving the upper half space, and performing least square optimization iterative solution with the actual ground gravity data to obtain the gravity anomaly information of the lower half space without Gibbs effect at the depth of the field source center;

[0009] S4. Take the depth of the field source center as the top interface of the field source, divide the top interface into grid prisms to construct the field source shape, perform density inversion on the top surface of each prism in the frequency domain, and obtain the density distribution at the gravity field source.

[0010] Preferably, in S1, the stable Fourier series solution under general conditions is obtained by using the Laplace equation of the gravity field and the corresponding boundary conditions, i.e., the gravity data observed on the ground, including:

[0011]

[0012] in, is the gravity field value at the measuring point (x, y, w), w is the elevation of the measuring point; n represents the discrete series order from 0 to N; is the frequency; μ n is the coefficient; e represents a natural constant.

[0013] Preferably, in S2, in the lower half space, a stable Fourier series carrying a regularization factor is selected to describe the gravity solution of the lower half space, and the center of the gravity anomaly corresponds to the field source center including:

[0014]

[0015] Among them, η n is the regularization factor, r is the multiple of the measurement point spacing, and increases with the increase of the series order. is the gravity field value at the measuring point (r, w), and μ0 is the coefficient.

[0016] Preferably, in S3, the data containing false anomalies at the depth center of the field source in the lower half space is used as the starting boundary condition of the upper half space, and the ground gravity anomaly data based on the boundary condition is calculated by solving the upper half space, and the least squares optimization iterative solution is performed with the actual ground gravity data to obtain the lower half space gravity anomaly information without Gibbs effect at the depth of the field source center, including:

[0017] S31. Assume that there are two planes ΨO and ΨP, representing the surface and a plane at a certain depth underground respectively. The outlier value of the kth point in Put it on the vertical projection point k of plane ΨP as the outlier on plane ΨP, let it be Right now: The spectrum of the outliers on ΨP is obtained by Fourier transform

[0018] S32. When there is no field source distribution between the two planes ΨO and ΨP, the gravity field value satisfies the Laplace equation. According to the solution formula of the first kind boundary value problem of the Laplace equation in the upper half space in the frequency domain:

[0019]

[0020] The outliers on the plane ΨP Spectrum Substitution Find the outlier value of point k on the ground ΨO If w=d, that is:

[0021]

[0022] Among them, s and t are spatial frequencies, representing the wave number respectively; Expressed as the gravity anomaly on the ground ΨO (w = 0) Fourier transform of -1 represents the inverse Fourier transform; d represents the height of the upper half space solution;

[0023] S33. Utilization and The difference between the outliers on the plane ΨP Correction is performed to obtain new outliers

[0024]

[0025] Among them, q represents the step length, which is generally taken as 1;

[0026] S34, repeat S32 and S33 to obtain the least squares optimization iterative solution formula:

[0027]

[0028] When it is obtained that wherein, represents a given number approaching zero; if meets the design accuracy, that is, when u < U, the iteration stops, and the number of iterations u is taken as 30 to 50 times,

[0029] S35. Calculate the gravity anomaly values of all points on the plane ΨP and use the solution formula of the first boundary value problem of the Laplace equation in the upper half space in the frequency domain:

[0030]

[0031] Perform Fourier inverse transform to obtain The gravity anomaly at the center depth d of the field source in the lower half space to be solved

[0032] Preferably, in the S4, taking the center depth of the field source as the top interface of the field source, constructing the field source shape from the top interface meshed prismatic bodies, and performing density inversion on the top surface of each prismatic body in the frequency domain to obtain the density distribution at the gravity field source, including:

[0033] Dividing the irregular geological body that causes the gravity anomaly into AB vertically small prismatic bodies with the same size and different centroid positions, where AB = A × B;

[0034] Let the centroid point coordinates of the small prismatic body be (x 0i , y 0i ), the remaining density be θ i , i = 1, 2, 3,..., AB, the side lengths be a and b respectively, the top surface burial depth of the prismatic body be d, and the height be cd, obtaining in the frequency domain the spectrum of the gravity anomaly generated by each prismatic body at the top interface (x d , y d , 0) of the known center depth of the field source Formula:

[0035]

[0036] According to the superposition property, the forward formula of the gravity anomaly of the field source with known physical property parameters is obtained, that is, the spectrum of the gravity anomaly G(x, y, 0) generated at the top interface (x d , y d , 0) of the known center depth of the field source for the entire depth layer

[0037]

[0038] Among them, s, t are spatial frequencies, representing the wave numbers of x and y respectively; G is the gravitational constant;

[0039] If the gravity anomaly G(x, y, 0) generated by the depth layer at the corresponding top interface is known, then d = 0. If d ≠ 0, that is, the burial depth is not zero, the gravity anomaly data of the top interface of the depth layer is calculated by the least squares optimization iterative solution in the lower half space, and then according to

[0040] make:

[0041]

[0042] Depend on

[0043] and It turns out that:

[0044]

[0045] according to Launched:

[0046]

[0047] right Perform inverse Fourier transform to obtain the top interface points (x d ,y d , 0) is the inversion formula for calculating the relevant physical parameters based on the known field source gravity anomaly:

[0048]

[0049] The apparent density value corresponding to the gravity anomaly is calculated based on the inversion formula for calculating relevant physical parameters based on the known field source gravity anomaly.

[0050] The present invention also provides a high-precision density inversion system for rapidly solving gravity anomalies based on regular stability, the system being used to implement any one of the methods described, comprising: a first calculation module, a second calculation module, an iteration module and an inversion module;

[0051] The first calculation module is used to obtain a stable Fourier series solution under general conditions through the Laplace equation of the gravity field and the corresponding boundary conditions, i.e., the gravity data observed on the ground;

[0052] The second calculation module is used to select a stable Fourier series carrying a regularization factor to describe the gravity solution of the lower half space in the lower half space, and the center of the gravity anomaly corresponds to the center of the field source;

[0053] The iterative module is used to use the data containing false anomalies at the depth center of the field source in the lower half space as the starting boundary condition of the upper half space, calculate the ground gravity anomaly data based on the boundary condition by solving the upper half space, and perform the least square optimization iterative solution with the actual ground gravity data to obtain the lower half space gravity anomaly information without Gibbs effect at the depth of the field source center;

[0054] The inversion module is used to use the center depth of the field source as the top interface of the field source, divide the top interface into grid prisms to construct the field source shape, perform density inversion on the top surface of each prism in the frequency domain, and obtain the density distribution at the gravity field source.

[0055] Preferably, in the first calculation module, obtaining a stable Fourier series solution under general conditions through the Laplace equation of the gravity field and the corresponding boundary conditions, i.e., the gravity data observed on the ground, includes:

[0056]

[0057] in, is the gravity field value at the measuring point (x, y, w), w is the elevation of the measuring point; n represents the discrete series order from 0 to N; is the frequency; μ n is the coefficient; e represents a natural constant.

[0058] Preferably, in the second calculation module, in the lower half space, a stable Fourier series carrying a regularization factor is selected to describe the gravity solution of the lower half space, and the center of the gravity anomaly corresponds to the field source center and includes:

[0059]

[0060] Among them, η n is the regularization factor, r is the multiple of the measurement point spacing, and increases with the increase of the series order. is the gravity field value at the measuring point (r, w), and μ0 is the coefficient.

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

[0062] The present invention provides a high-precision density inversion method based on regularized stable and fast solution of gravity anomaly. Compared with the traditional transformation method, the present invention first uses the ground gravity anomaly as the first-class boundary value condition to perform regularized stable solution of the Laplace lower half space, and adopts a function carrying a regularization factor to correct the frequency distribution of the series terms of the gravity field value solved in the lower half space, effectively avoiding the divergence of the result, and obtaining high-resolution anomaly data at the periphery of the field source depth. On the basis of the above-mentioned field source regularization solution, it is proposed to use the data carrying false anomalies at the field source imaging depth as the boundary value condition, and again solve the upper half space of the first-class boundary value problem of the Laplace equation to calculate the gravity anomaly at the ground, and perform the least squares optimization iterative solution with the actual ground anomaly to obtain high-precision anomaly imaging data at the field source depth without Gibbs effect false anomalies.

[0063] Due to the anomalies generated by the layer sources at different depths underground, the lower half-space solution is often required to reach the top of the layers at different depths. The stability and depth of the solution will directly affect the quality of the final result. According to the research, the instability of the lower half-space solution of the potential field and the limited depth of the lower half-space solution of the classical FFT method are solved. At this time, the spatial domain least squares optimization iterative solution method is adopted. It is based on the stable FFT upper half-space solution. The principle is simple and the extension is stable. By assuming that there are two planes ΨO and ΨP, representing the planes at the surface and a certain depth underground respectively, and there is no source distribution between the two planes, that is, the passive space, the lower half-space iterative solution can be completed based on the known depth of the upper half-space solution according to the least squares optimization iterative solution flow chart, and high-precision anomaly imaging data at the source depth without Gibbs effect false anomalies can be obtained. Then, this depth is used as the top interface of the source, and the grid prism is divided by the top interface to construct the source morphology. The density inversion of the top surface of each prism is performed in the frequency domain, thereby obtaining the density distribution of the source morphology. The present invention performs Fourier transform on the field source gravity data to obtain its frequency domain signal, that is, the algorithm is implemented in the frequency domain, and the inversion is performed at a known field source depth. Therefore, the invention is a fast, high-precision, high-resolution characteristic gravity density inversion method.

[0064] The present invention adopts a high-precision density inversion method based on regularization to quickly solve gravity anomalies in a stable manner, and obtains high-precision anomaly imaging data of the source depth without Gibbs effect false anomalies. Moreover, since the algorithm is implemented in the frequency domain and the inversion is performed at a known source depth, the calculation speed of this method is very fast, which effectively enhances the stability and accuracy of the calculation results. In theory, as long as the depth determined by the regularization method is accurate, the true physical parameters of the gravity field with higher accuracy can be obtained, which provides a new idea for solving the density parameters of the gravity field source. BRIEF DESCRIPTION OF THE DRAWINGS

[0065] In order to more clearly illustrate the technical solution of the present invention, 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 creative labor.

[0066] Figure 1 It is the technical roadmap of the present invention;

[0067] Figure 2 are specific parameters of a single sphere model according to an embodiment of the present invention;

[0068] Figure 3 A schematic diagram of the spatial position of a model according to an embodiment of the present invention;

[0069] Figure 4 A schematic diagram of the least squares optimization iterative solution space of an embodiment of the present invention;

[0070] Figure 5 This is a flowchart of the iterative solution of least squares optimization according to an embodiment of the present invention;

[0071] Figure 6 This is a gravity anomaly imaging diagram of a spherical model according to an embodiment of the present invention;

[0072] Figure 7 The main profile curve of the embodiment of the present invention and the regularized stable solution cross-section diagram of its Laplace lower half space;

[0073] Figure 8 This is a gravity anomaly imaging diagram of the top interface at the depth of the field source center (left) and an apparent density inversion result diagram (right) of an embodiment of the present invention. DETAILED DESCRIPTION

[0074] The following will be combined with the drawings in the embodiments of the present invention to clearly and completely describe the technical solutions in the embodiments of the present invention. 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 creative work are within the scope of protection of the present invention.

[0075] It should be noted that, unless otherwise defined, the technical terms or scientific terms used in the embodiments of the present disclosure should be understood by people with ordinary skills in the field to which the present disclosure belongs. The "first", "second" and similar words used in the embodiments of the present disclosure do not indicate any order, quantity or importance, but are only used to distinguish different components. "Including" or "comprising" and similar words mean that the elements or objects appearing before the word cover the elements or objects listed after the word and their equivalents, without excluding other elements or objects. "Connect" or "connected" and similar words are not limited to physical or mechanical connections, but can include electrical connections, whether direct or indirect. "Up", "down", "left", "right" and the like are only used to indicate relative positional relationships. When the absolute position of the described object changes, the relative positional relationship may also change accordingly.

[0076] 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.

[0077] Embodiment 1

[0078] like Figure 1 As shown, the embodiment of the present invention provides a high-precision density inversion method for quickly solving gravity anomalies based on regular stability. The method includes:

[0079] S1. First, the Fourier series transformation of the observed gravity anomaly data is performed, and the Laplace equation of the gravity field and its boundary conditions (gravity data observed on the ground ), we can obtain the stable Fourier series solution under general conditions

[0080] S2, in the lower half space, as the calculation approaches and reaches the field source, its boundary value conditions become unsatisfied, and the solution in S1 becomes ill-posed. In order to overcome the oscillation effect caused by ill-posedness in the solution process, a stable Fourier series with a regularization factor is selected to describe the gravity solution in the lower half space, and the center of the gravity anomaly corresponds to the center of the field source;

[0081] S3. The stable solution of the lower half space obtained by using the above-mentioned regularization factor contains a false anomaly of the Gibbs effect. In order to eliminate this false anomaly, the data containing the false anomaly at the depth center of the field source in the lower half space is used as the starting boundary condition of the upper half space. By solving the upper half space, the ground gravity anomaly data based on the boundary condition is calculated, and the least squares optimization iterative solution is performed with the actual ground gravity data. After multiple fittings of the above two data, the gravity anomaly information of the lower half space without the Gibbs effect at the depth of the field source center is obtained;

[0082] S4. Finally, the center depth of the field source is used as the top interface of the field source. The top interface is used to divide the grid prisms to construct the field source shape. The density of the top surface of each prism is inverted in the frequency domain, so that the density distribution at the gravity field source can be quickly obtained.

[0083] S5. Finally, the physical parameter constrained inversion of the single model gravity anomaly is implemented.

[0084] In this embodiment, in S1, the stable Fourier series solution under general conditions is obtained by using the Laplace equation of the gravity field and the corresponding boundary conditions, i.e., the gravity data observed on the ground, including:

[0085] From the Laplace equation of gravity field and its boundary conditions (gravity data observed on the ground):

[0086]

[0087] The Taylor expansion equation for the gravity field value is as follows:

[0088]

[0089] Let D be the number of measuring points, Δr be the spacing between measuring points, have:

[0090] r=(n-1)·Δr

[0091]

[0092] Among them, r is a multiple of the distance between the measuring points, which increases with the increase of the order of the series.

[0093] Using the separation of variables method and boundary conditions, we can obtain a stable Fourier series solution under general conditions:

[0094]

[0095] in, is the gravity field value at the measuring point (x, y, w), is the ground-observed gravity field value; w is the elevation of the observation point, which is the ground-observed data at this time, that is, w0=0; n represents the discrete series order from 0 to N; is the frequency; μ n is the coefficient; e represents the natural constant, is the gravity field value at the measuring point (r, w).

[0096] Since the frequency is related to the depth of the solution in the lower half space, the shallower the information, the higher the frequency domain, that is, The larger the value, the infinite value of the gravity field will be, resulting in divergence of the acquired gravity anomaly data. Secondly, when solving the lower half space, the greater the depth, that is, the greater the w value, the infinite value of the field will be, resulting in data divergence. The exponential part of will have a great influence on the solution.

[0097] In this embodiment, in S2, in the lower half space, a stable Fourier series carrying a regularization factor is selected to describe the gravity solution of the lower half space, and the center of the gravity anomaly corresponds to the field source center including:

[0098] Add regularization factor It is related to frequency Calculate the function related to a certain depth (w) and coefficient (μ) in the space from below the observation surface to above the field source. At this time, and w as the denominator of the regularization factor, can effectively suppress the information of the above-mentioned exponential increase, and can be used to effectively correct the frequency division of the series term of the potential field in the lower half space;

[0099] Therefore, the stable Fourier series after correction with regularization factor is:

[0100]

[0101] Among them, η n is the regularization factor and μ0 is the coefficient.

[0102] Thus, the gravity anomaly at the depth of the field source is obtained. Practice has proved that after adding the regularization factor, not only can the potential field be guaranteed to be non-singular when passing through the field source during the solution of the lower half space, and the complete field value distribution of the lower half space can be obtained, but the Gibbs effect can also be not particularly significant.

[0103] In this embodiment, in S3, the data containing false anomalies at the depth center of the field source in the lower half space is used as the starting boundary condition of the upper half space. By solving the upper half space, the ground gravity anomaly data based on the boundary condition is calculated, and the least squares optimization iterative solution is performed with the actual ground gravity data to obtain the lower half space gravity anomaly information without Gibbs effect at the depth of the field source center, including:

[0104] Taking the gravity anomaly at the depth of the field source as the boundary condition, the calculation method of the upper half space is used to calculate the ground gravity anomaly data, and the calculated ground gravity data and the ground measured gravity data are iteratively solved by least square optimization. Through continuous data fitting, the gravity data at the field source can be corrected to obtain high-precision anomaly imaging data at the depth of the field source without Gibbs effect false anomalies, as shown in the attached figure of the specification. Figure 4 and Figure 5 As shown, it specifically includes the following steps:

[0105] S31. Assume that there are two planes ΨO and ΨP, representing the surface and a plane at a certain depth underground respectively. The outlier value of the kth point in Put it on the vertical projection point k of plane ΨP as the outlier on plane ΨP, let it be Right now: The spectrum of the outliers on ΨP is obtained by Fourier transform

[0106] S32. When there is no field source distribution between the two planes ΨO and ΨP, the gravity field value satisfies the Laplace equation. According to the solution formula of the first kind boundary value problem of the Laplace equation in the upper half space in the frequency domain:

[0107]

[0108] in, express Spectral expression solved in the frequency domain.

[0109] The outliers on the plane ΨP Spectrum Substitution Find the outlier value of point k on the ground ΨO If w=d, that is:

[0110]

[0111] Among them, s and t are spatial frequencies, representing the wave numbers of x and y respectively; Expressed as the gravity anomaly on the ground ΨO (w = 0) Fourier transform of -1 represents the inverse Fourier transform; d represents the height of the upper half space solution;

[0112] S33. Utilization and The difference between the outliers on the plane ΨP Correction is performed to obtain new outliers

[0113]

[0114] Among them, q represents the step length, which is generally taken as 1;

[0115] S34, repeat S32 and S33 to obtain the least squares optimization iterative solution formula:

[0116]

[0117] in, Denote the outlier of the k-th point on the plane ΨP at the u-th iteration; Denote the outlier of the k-th point on the ground ΨO at the u-th iteration; U represents the maximum number of iterations.

[0118] When , from obtain where δ represents a given number approaching zero; if meets the design accuracy, that is, when u < U, stop the iteration. The number of iterations u is generally taken as 30 - 50 times.

[0119] S35. Calculate the gravity anomaly values of all points on the plane ΨP And use the solution formula in the frequency domain for the first boundary value problem of the Laplace equation in the upper half space:

[0120]

[0121] Perform the inverse Fourier transform to obtain The gravity anomaly at the depth d of the center of the field source in the lower half space to be solved

[0122] In this embodiment, in S4, taking the depth of the center of the field source as the top interface of the field source, constructing the field source morphology from the top interface meshed prismatic bodies, and performing density inversion on the top surface of each prismatic body in the frequency domain to obtain the density distribution at the gravity field source including:

[0123] Divide the irregular geological body causing the gravity anomaly into AB vertically small prismatic bodies with the same size and different centroid positions, where AB = A × B;

[0124] Let the centroid point coordinates of the small prismatic body be (x 0i , y 0i ), the residual density be θ i , i = 1, 2, 3,..., AB, the side lengths be a and b respectively, the top surface burial depth of the prismatic body be d, and the height be cd. Obtain the spectrum of the gravity anomaly generated by each prismatic body at the known top interface (x d , y d , 0) of the center depth of the field source in the frequency domain Formula:

[0125]

[0126] According to the superposition property, obtain the forward modeling formula of the gravity anomaly of the field source with known physical property parameters, that is, the entire depth layer at the known top interface (x d , y d, 0) is the frequency spectrum of the gravity anomaly G(x, y, 0)

[0127]

[0128] Among them, s, t are spatial frequencies, representing the wave numbers of x and y respectively; G is the gravitational constant;

[0129] If the gravity anomaly G(x, y, 0) generated by the depth layer at the corresponding top interface is known, then d = 0. If d ≠ 0, that is, the burial depth is not zero, the gravity anomaly data of the top interface of the depth layer is calculated by the least squares optimization iterative solution in the lower half space, and then according to

[0130] make:

[0131]

[0132] Depend on

[0133] and It turns out that:

[0134]

[0135] according to Launched:

[0136]

[0137] right Perform inverse Fourier transform to obtain the top interface points (x d ,y d , 0) is the inversion formula for calculating the relevant physical parameters based on the known field source gravity anomaly:

[0138]

[0139] According to the inversion formula of the known source gravity anomaly to infer the relevant physical parameters, the apparent density value corresponding to the gravity anomaly is calculated. Since the data all contain the corresponding depth, it can be quickly inverted to obtain high-precision anomaly imaging data of the source depth without Gibbs effect false anomalies.

[0140] In this embodiment, in S5, physical property parameter constrained inversion of single model gravity anomaly is implemented.

[0141] The present invention designs a 101×101 grid, and the point spacing and line spacing are both 1 (unit: meter), that is, the measurement area size is 100×100m 2 , the specific parameters of a single sphere model are as follows Figure 2As shown, the spatial position of the model is Figure 3 shown.

[0142] According to the forward calculation of the above model, the corresponding gravity anomaly is obtained, such as Figure 6 As shown; Figure 7 It is the main profile curve of the gravity anomaly and the regularized stable solution of the lower half space of Laplace to obtain the information section of each depth layer. It can be seen that the depth of the anomaly center is consistent with the given model. However, due to the influence of the Gibbs truncation effect of the frequency domain solution, there is a weak anomaly belt with alternating positive and negative in the shape of an "eight" with the main source body distribution area as the center, resulting in poor horizontal resolution of field source body imaging, which increases the difficulty of using the spatial distribution characteristics of the gravity potential field;

[0143] In view of this, based on the above-mentioned field source regularization solution, the present invention proposes to use the data carrying false anomalies at the field source imaging depth as the boundary value condition, and again solve the upper half space of the first-kind boundary value problem of the Laplace equation to calculate the gravity anomaly at the ground, and perform the least squares optimization iterative solution with the actual ground anomaly to obtain high-precision abnormal imaging data at the field source depth without Gibbs effect false anomalies. The processed spatial schematic diagram is shown in Figure 4 The specific process is as shown in Figure 5 As shown;

[0144] The obtained source center depth is used as the top interface of the source, and the grid prism is divided by the top interface to construct the source shape. The density inversion of the top surface of each prism is performed in the frequency domain, thereby obtaining the density distribution of the source shape. It is easy to see from the apparent density inversion result diagram that the density is basically consistent with the assumed single sphere model. Figure 8 The gravity anomaly (left) and apparent density inversion result (right) of the top interface at the depth of the field source center. Since the algorithm is implemented in the frequency domain and the inversion is performed at a known field source depth, the invention is a fast, high-precision, high-resolution characteristic gravity density inversion method.

[0145] In summary, the invention provides a high-precision density inversion method for rapidly solving gravity anomalies based on regularization, stability, and has the following advantages compared with the traditional transformation method:

[0146] 1) The accuracy of gravity anomaly parameter inversion is related to the accuracy of the constrained depth. Complex anomalies will reduce the accuracy of the constrained depth. However, since the ground gravity anomaly is used as the first-class boundary condition, the regularized stable solution of the Laplace lower half space is performed, and the frequency distribution of the series terms of the gravity field values ​​solved in the lower half space is corrected by using a function carrying a regularization factor, which effectively avoids the divergence of the results. The high-resolution information near the top interface obtained in this way can still ensure that the physical parameters have a high lateral resolution;

[0147] 2) Based on the above-mentioned source regularization solution, it is proposed to use the data carrying false anomalies at the source imaging depth as boundary conditions, and again calculate the gravity anomaly at the ground through the solution of the first-kind boundary value problem of the Laplace equation in the upper half space, and perform the least squares optimization iterative solution with the actual ground anomaly, so as to obtain high-precision anomaly imaging data at the source depth without Gibbs effect false anomalies;

[0148] 3) The algorithm is implemented in the frequency domain, and the inversion is performed at a known source depth. This method has a fast calculation speed and effectively enhances the stability and accuracy of the calculation results. In theory, as long as the source center depth determined by the regularized stable solution of the Laplace lower half space is accurate, the true physical parameters of the gravity field can be obtained with high accuracy.

[0149] 4) The high-precision density inversion method based on regular stable and rapid solution of gravity anomalies is conducive to improving the inversion accuracy of depth-constrained physical parameters of complex geological models, which can effectively improve the effectiveness and accuracy of deep mineral resource exploration, increase the success rate of prospecting, reduce exploration costs, and provide important technical support for my country's resource development.

[0150] It should be noted that the method of the embodiment of the present disclosure can be performed by a single device, such as a computer or a server. The method of the present embodiment can also be applied in a distributed scenario and completed by multiple devices cooperating with each other. In the case of such a distributed scenario, one of the multiple devices can only perform one or more steps in the method of the embodiment of the present disclosure, and the multiple devices will interact with each other to complete the described method.

[0151] It should be noted that some embodiments of the present disclosure are described above. Other embodiments are within the scope of the attached claims. In some cases, it should be understood that the size of the serial number of each step in the above embodiment does not mean the order of execution, and the execution order of each process should be determined by its function and internal logic, and should not constitute any limitation on the implementation process of the embodiment of the present invention. The actions or steps recorded in the claims can be performed in an order different from that in the above embodiment and still achieve the desired results. In addition, the process depicted in the accompanying drawings does not necessarily require the specific order or continuous order shown to achieve the desired results. In some embodiments, multitasking and parallel processing are also possible or may be advantageous.

[0152] Embodiment 2

[0153] Based on the same inventive concept, corresponding to any of the above-mentioned embodiments and methods, the present disclosure further provides a high-precision density inversion system for rapidly solving gravity anomalies based on regular stability, the system being used to implement any of the above-mentioned methods, comprising: a first calculation module, a second calculation module, an iteration module and an inversion module;

[0154] The first calculation module is used to obtain a stable Fourier series solution under general conditions through the Laplace equation of the gravity field and the corresponding boundary conditions, i.e., the gravity data observed on the ground;

[0155] The second calculation module is used to select a stable Fourier series carrying a regularization factor to describe the gravity solution of the lower half space in the lower half space, and the center of the gravity anomaly corresponds to the center of the field source;

[0156] The iteration module is used to use the data containing false anomalies at the depth center of the lower half-space field source as the starting boundary condition of the upper half-space, calculate the ground gravity anomaly data based on the boundary condition by solving the upper half-space, and perform the least square optimization iterative solution with the actual ground gravity data to obtain the lower half-space gravity anomaly information without Gibbs effect at the depth of the field source center;

[0157] The inversion module is used to take the center depth of the field source as the top interface of the field source, divide the top interface into grid prisms to construct the field source shape, perform density inversion on the top surface of each prism in the frequency domain, and obtain the density distribution at the gravity field source.

[0158] In this embodiment, in the first calculation module, the stable Fourier series solution under general conditions is obtained through the Laplace equation of the gravity field and the corresponding boundary conditions, i.e., the gravity data observed on the ground, including:

[0159]

[0160] in, is the gravity field value at the measuring point (x, y, w), w is the elevation of the measuring point; n represents the discrete series order from 0 to N; is the frequency; μ n is the coefficient; e represents a natural constant.

[0161] In this embodiment, in the second calculation module, in the lower half space, a stable Fourier series carrying a regularization factor is selected to describe the gravity solution of the lower half space, and the center of the gravity anomaly corresponds to the field source center including:

[0162]

[0163] Among them, η n is the regularization factor.

[0164] The system of the above embodiment is used to implement a high-precision density inversion method for quickly solving gravity anomalies based on regularization stability in any of the above embodiments, and has the beneficial effects of the corresponding method embodiments, which will not be repeated here.

[0165] It should be noted that the above-mentioned high-precision density inversion system for rapidly solving gravity anomalies based on regular stable formula is embodied in the form of functional units. The term "module" here can be implemented in the form of software and / or hardware, and is not specifically limited to this.

[0166] For example, a "module" may be a software program, a hardware circuit, or a combination of the two that implements the above functions. The hardware circuit may include an application specific integrated circuit (ASIC), an electronic circuit, a processor (such as a shared processor, a dedicated processor, or a group processor, etc.) and a memory for executing one or more software or firmware programs, a combined logic circuit, and / or other suitable components that support the described functions.

[0167] The embodiments of the present disclosure are intended to cover all such substitutions, modifications and variations that fall within the broad scope of the appended claims. Therefore, any omissions, modifications, equivalent substitutions, improvements, etc. made within the spirit and principles of the embodiments of the present disclosure should be included in the scope of protection of the present disclosure.

Claims

1. A high-precision density inversion method for rapidly solving gravity anomalies based on regularization stability, characterized in that: The method comprises: S1. Through the Laplace equation of gravity field and the corresponding boundary conditions, i.e. the gravity data observed on the ground, a stable Fourier series solution under general conditions is obtained; S2. In the lower half space, a stable Fourier series with a regularization factor is selected to describe the gravity solution of the lower half space, and the center of the gravity anomaly corresponds to the center of the field source; S3, taking the data containing false anomalies at the depth center of the field source in the lower half space as the starting boundary condition of the upper half space, calculating the ground gravity anomaly data based on the boundary condition by solving the upper half space, and performing least square optimization iterative solution with the actual ground gravity data to obtain the gravity anomaly information of the lower half space without Gibbs effect at the depth of the field source center; S4. Take the depth of the field source center as the top interface of the field source, divide the top interface into grid prisms to construct the field source shape, perform density inversion on the top surface of each prism in the frequency domain, and obtain the density distribution at the gravity field source.

2. The method according to claim 1, characterized in that In S1, the stable Fourier series solutions under general conditions are obtained through the Laplace equation of the gravity field and the corresponding boundary conditions, i.e., the gravity data observed on the ground, including: in, is the gravity field value at the measuring point (x, y, w), w is the height of the measuring point; n represents the discrete series order from 0 to N; ζ n is the frequency; μ n is the coefficient; e represents a natural constant.

3. The method according to claim 2, characterized in that In S2, in the lower half space, a stable Fourier series carrying a regularization factor is selected to describe the gravity solution of the lower half space, and the center of the gravity anomaly corresponds to the field source center including: Among them, η n is the regularization factor, r is the multiple of the measurement point spacing, which increases with the increase of the series order. is the gravity field value at the measuring point (r, w), and μ0 is the coefficient.

4. The method according to claim 3, characterized in that In S3, the data containing false anomalies at the depth center of the field source in the lower half space is used as the starting boundary condition of the upper half space. By solving the upper half space, the ground gravity anomaly data based on the boundary condition is calculated, and the least squares optimization iterative solution is performed with the actual ground gravity data to obtain the lower half space gravity anomaly information without Gibbs effect at the depth of the field source center, including: S31. Assume that there are two planes ΨO and ΨP, representing the surface and a plane at a certain depth underground respectively. The outlier value of the kth point in Put it on the vertical projection point k of plane ΨP as the outlier on plane ΨP, let it be Right now: The spectrum of the outliers on ΨP is obtained by Fourier transform S32. When there is no field source distribution between the two planes ΨO and ΨP, the gravity field value satisfies the Laplace equation. According to the solution formula of the first kind boundary value problem of the Laplace equation in the upper half space in the frequency domain: The outliers on the plane ΨP Spectrum Substitution Find the outlier value of point k on the ground ΨO If w=d, that is: Among them, s and t are spatial frequencies, representing the wave numbers of x and y respectively; Expressed as the gravity anomaly on the ground ΨO (w = 0) Fourier transform of -1 represents the inverse Fourier transform; d represents the height of the upper half space solution; S33. Utilization and The difference between the outliers on the plane ΨP Make corrections to get new outliers Among them, q represents the step length, which is generally taken as 1; S34, repeat S32 and S33 to obtain the least squares optimization iterative solution formula: When then, from obtain wherein, δ represents a given number approaching zero; if meets the design accuracy, that is, when u < U, then stop the iteration, and the number of iterations u is taken as 30 to 50 times, S35. Calculate the gravity anomaly values ​​of all points on the plane ΨP And use the solution formula of the first kind of boundary value problem of Laplace equation in the upper half space in the frequency domain: Perform inverse Fourier transform to obtain The gravity anomaly at the depth d of the source center in the lower half space solved 5. The method according to claim 4, characterized in that In S4, the depth of the field source center is used as the top interface of the field source, the top interface is divided into grid prisms to construct the field source shape, and the density of the top surface of each prism is inverted in the frequency domain to obtain the density distribution at the gravity field source, including: Will cause gravity anomaly The irregular geological body is divided into AB small vertical prisms of finite extension with the same size and different centroid positions, where AB = A × B; Assume that the coordinates of the center of mass of the small prism are (x 0i ,y 0i ), the residual density is θ i , i = 1, 2, 3, ..., AB, the side lengths are a and b respectively, the top surface of the prism is buried at a depth of d and the height is cd, and the top interface (x d ,y d , 0) produces a frequency spectrum of gravity anomaly formula: According to the superposition property, the forward formula of the gravity anomaly of the field source with known physical parameters is obtained, that is, the top interface (x d ,y d , 0) is the frequency spectrum of the gravity anomaly G(x, y, 0) Among them, s, t are spatial frequencies, representing the wave numbers of x and y respectively; G is the gravitational constant; If the gravity anomaly G(x, y, 0) generated by the depth layer at the corresponding top interface is known, then d = 0. If d ≠ 0, that is, the burial depth is not zero, the gravity anomaly data of the top interface of the depth layer is calculated by the least squares optimization iterative solution in the lower half space, and then according to make: Depend on and It turns out that: according to Launched: right Perform inverse Fourier transform to obtain the top interface points (x d ,y d , 0) is the inversion formula for calculating the relevant physical parameters based on the known field source gravity anomaly: The apparent density value corresponding to the gravity anomaly is calculated based on the inversion formula for calculating relevant physical parameters based on the known field source gravity anomaly.

6. A high-precision density inversion system based on regular stable rapid solution of gravity anomaly, the system is used to implement the method according to any one of claims 1 to 5, characterized in that: include: A first calculation module, a second calculation module, an iteration module and an inversion module; The first calculation module is used to obtain a stable Fourier series solution under general conditions through the Laplace equation of the gravity field and the corresponding boundary conditions, i.e., the gravity data observed on the ground; The second calculation module is used to select a stable Fourier series carrying a regularization factor to describe the gravity solution of the lower half space in the lower half space, and the center of the gravity anomaly corresponds to the center of the field source; The iterative module is used to use the data containing false anomalies at the depth center of the field source in the lower half space as the starting boundary condition of the upper half space, calculate the ground gravity anomaly data based on the boundary condition by solving the upper half space, and perform the least square optimization iterative solution with the actual ground gravity data to obtain the lower half space gravity anomaly information without Gibbs effect at the depth of the field source center; The inversion module is used to use the center depth of the field source as the top interface of the field source, divide the top interface into grid prisms to construct the field source shape, perform density inversion on the top surface of each prism in the frequency domain, and obtain the density distribution at the gravity field source.

7. The system according to claim 6, characterized in that In the first calculation module, the stable Fourier series solution under general conditions is obtained through the Laplace equation of the gravity field and the corresponding boundary conditions, i.e., the gravity data observed on the ground, including: in, is the gravity field value at the measuring point (x, y, w), w is the height of the measuring point; n represents the discrete series order from 0 to N; ζ n is the frequency; μ n is the coefficient; e represents a natural constant.

8. The system according to claim 7, characterized in that In the second calculation module, in the lower half space, a stable Fourier series carrying a regularization factor is selected to describe the gravity solution of the lower half space, and the center of the gravity anomaly corresponds to the field source center including: Among them, η n is the regularization factor, r is the multiple of the measurement point spacing, which increases with the increase of the series order. is the gravity field value at the measuring point (r, w), and μ0 is the coefficient.

Images