A high-precision density inversion method and system for stably and rapidly solving gravity anomalies based on regularities

Abstract

The present invention discloses a high-precision density inversion method and system for stably and rapidly solving gravity anomalies based on regularization. On the basis of the above-mentioned regularization solution of the field source, it is proposed to use the data carrying false anomalies at the imaging depth of the field source as the boundary value condition, and then calculate the gravity anomaly at the ground through the solution of the upper half-space of the first boundary value problem of the Laplace equation, and perform least-squares optimization iteration with the actual ground anomaly to obtain high-precision anomaly imaging data at the depth of the field source without Gibbs effect false anomalies. Taking this depth as the top interface of the field source, the field source morphology is constructed by dividing the grid prisms at the top interface, and the density inversion is carried out on the top surface of each prism in the frequency domain to obtain the density distribution of the field source morphology. Since this algorithm is implemented in the frequency domain and the inversion is carried out at the known depth of the field source, the present invention is a fast, high-precision, and high-resolution gravity density inversion method, providing a new idea for realizing the solution of the density parameters of the gravity field source.

Description

Technical Field

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

[0002] Gravity exploration, as a commonly used geophysical exploration method, has a wide range of applications in aspects such as energy and mineral resource exploration, geological mapping, and engineering construction. The processing and conversion of gravity anomaly data are the key parts of the theory of gravity data interpretation. Among them, taking the ground gravity anomaly as the first-kind boundary condition and solving the Laplace equation in the upper half-space (i.e., far from the field source) can obtain a regular and stable solution. When solving the Laplace equation in the lower half-space (i.e., close to the field source), on the one hand, higher-resolution gravity anomaly information can be obtained, but the lower half-space is restricted by the boundary value of the Laplace equation, and the gravity anomaly information at the field source cannot be obtained. To ensure that the potential field solved at the field source is non-singular during the solution process of the lower half-space and obtain the complete field value distribution of the lower half-space, it is necessary to optimize and improve the gravity field solution method in the lower half-space.

[0003] The distribution characteristics of the lower half-space of the gravity and magnetic potential fields are highly intuitive and are an important basis for studying the field source body. By taking the ground-observed gravity anomaly as the first-kind boundary condition and performing regularized and stable solution of the Laplace lower half-space, high-resolution anomaly data around the field source depth can be obtained. However, due to the influence of the Gibbs truncation effect in the frequency-domain solution, on both sides of the main anomaly band, there are often weak anomaly bands with positive and negative alternations in a "figure-eight" shape centered on the distribution area of the main source body. It may be a single anomaly or may be superimposed with the main anomaly of another source body to form a composite anomaly, thus increasing the difficulty of interpreting the field source body. Summary of the Invention

[0004] To solve the above technical problems, the present invention provides a high-precision density inversion method and system for stably and rapidly solving gravity anomalies based on regularization. On the basis of the regularized solution of the field source, it is proposed to use the data carrying false anomalies at the imaging depth of the field source as the boundary condition, and then through the solution of the first-kind boundary value problem of the Laplace equation in the upper half-space, calculate the gravity anomaly at the ground, and perform least-squares optimization iteration with the actual ground anomaly to obtain high-precision anomaly imaging data at the field source depth without the false anomalies of the Gibbs effect. Taking this depth as the top interface of the field source, the field source morphology is constructed by dividing the grid prisms at the top interface, and the density inversion is performed on the top surface of each prism in the frequency domain, thereby obtaining the density distribution of the field source morphology. Since this algorithm is implemented in the frequency domain and the inversion is performed at the known field source depth, the present invention is a fast, high-precision, and high-resolution gravity density inversion method, providing a new idea for the comprehensive interpretation method of the field source body from both physical property and geometric parameter aspects.

[0005] The present invention discloses a high-precision density inversion method for stably and rapidly solving gravity anomalies based on regularization, and the method includes:

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

[0007] S2. In the lower half-space, select a stable Fourier series with a regularization factor 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;

[0008] S3. Use the data containing false anomalies at the center depth of the field source in the lower half-space as the starting boundary value condition for the upper half-space. Through the solution of the upper half-space, calculate the ground gravity anomaly data based on this boundary value condition, and perform least squares optimization iteration with the actual ground gravity data to obtain the gravity anomaly information in 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, construct the field source morphology from the top interface by dividing the grid prisms, and perform density inversion on the top surface of each prism in the frequency domain to obtain the density distribution at the gravity field source.

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

[0011]

[0012] Wherein, is the gravitational field value at the measurement point (x, y, w), w is the elevation of the observation point; n represents the discrete series order taken from 0 to N; is the frequency; μ n is the coefficient; e represents the natural constant.

[0013] Preferably, in S2, in the lower half-space, selecting a stable Fourier series with a regularization factor 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, includes:

[0014]

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

[0016] Preferably, in S3, the data including false anomalies at the depth center of the lower half-space field source are used as the starting boundary value conditions for the upper half-space. Through the solution of the upper half-space, the ground gravity anomaly data based on these boundary value conditions are calculated, and the least squares optimization iteration solution is carried out 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, including:

[0017] S31. Suppose there are two planes ΨO and ΨP, representing the ground surface and a plane at a certain depth underground respectively. The anomaly value at the k-th point in is placed at the vertical projection point k of the plane ΨP as the anomaly value on the plane ΨP, denoted as That is: The spectrum of the anomaly value on ΨP is obtained through 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 of boundary value problem of the Laplace equation in the upper half-space in the frequency domain:

[0019]

[0020] Substitute the spectrum of the anomaly value on the plane ΨP into to obtain the anomaly value at the k-th point on the ground ΨO. If w = d, that is:

[0021]

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

[0023] S33. Use the difference between and to correct the anomaly value on the plane ΨP to obtain a new anomaly value

[0024]

[0025] where q represents the step size, generally taken as 1;

[0026] S34. Repeat S32 and S33 to obtain the least squares optimization iteration solution formula:

[0027]

[0028] When it is, from obtain 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 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

[0032] Preferably, in the 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 by dividing the grid prisms, and performing density inversion on the top surface of each prism in the frequency domain to obtain the density distribution at the gravity field source including:

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

[0034] Let the centroid point coordinates of the small prism 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 prism be d, and the height be cd, to obtain the spectrum of the gravity anomaly generated by each prism at the known top interface (x d , y d , 0) of the center depth of the field source in the frequency domain Formula:

[0035]

[0036] According to the superposition property, the forward modeling 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 known top interface (x d , y d , 0) of the center depth of the field source in the entire depth layer

[0037]

[0038] where s and 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 at the corresponding top interface of this depth layer is known, then d = 0. If d ≠ 0, that is, the buried depth is not zero, the gravity anomaly data of the top interface of this depth layer is calculated by least squares optimization iteration in the lower half space, and then according to

[0040] Let:

[0041]

[0042] From

[0043] and we get:

[0044]

[0045] According to it is deduced that:

[0046]

[0047] For performing the Fourier inverse transform, the residual density values of each point (x d , y d , 0) on the top interface of this depth layer are obtained, that is, the inversion formula for calculating relevant physical property parameters from the known field source gravity anomaly:

[0048]

[0049] According to the inversion formula for calculating relevant physical property parameters from the known field source gravity anomaly, the apparent density value corresponding to the gravity anomaly is calculated.

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

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

[0052] The second calculation module is used to select a stable Fourier series with a regularization factor 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;

[0053] The iterative module is used to take the data with false anomalies in the depth center of the field source in the lower half space as the starting boundary value condition for the upper half space. Through the solution of the upper half space, the ground gravity anomaly data based on this boundary value condition is calculated, and the least squares optimization iteration solution is carried out 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;

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

[0055] Preferably, in the first calculation module, through the Laplace equation of the gravity field and the corresponding boundary conditions, that is, the gravity data observed on the ground, the stable Fourier series solution under general conditions is obtained, including:

[0056]

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

[0058] Preferably, in the second calculation module, in the lower half space, 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 field source center, including:

[0059]

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

[0061] Compared with the prior art, the beneficial effects of the present invention are:

[0062] A high-precision density inversion method for regular and stable rapid solution of gravity anomalies provided by the present invention. Compared with traditional transformation methods, the present invention first uses the ground gravity anomaly as the first-kind boundary condition to perform regularized stable solution of the lower half space of Laplace, and adopts a function carrying a regularization factor to correct the series terms of the gravity field value solved in the lower half space by frequency, 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 regularized solution of the field source, it is proposed to use the data carrying false anomalies at the imaging depth of the field source as the boundary condition, and then calculate the gravity anomaly at the ground through the solution of the first-kind boundary value problem of the upper half space of the Laplace equation, and perform least-square optimization iteration with the actual ground anomaly to obtain high-precision anomaly imaging data at the field source depth without Gibbs effect false anomalies.

[0063] Since separating the anomalies generated by layer sources at different depths underground on the ground often requires solving to the top of different depth layers in the lower half space, and the stability and depth of the solution will directly affect the quality of the final result. According to research on the instability of the solution of the potential field lower half space and the limited depth of the solution of the classical FFT method in the lower half space, at this time, a spatial domain least-square optimization iteration solution method is adopted, which is realized based on the stable FFT upper half space solution. The principle is simple and the downward continuation is stable. By assuming the existence of two planes ΨO and ΨP, representing the ground surface and a plane at a certain depth underground respectively, and there is no field source distribution between the two planes, that is, a source-free space, then the lower half space iteration solution based on the known depth of the upper half space solution can be completed according to the least-square optimization iteration solution flow chart, and high-precision anomaly imaging data at the field source depth without Gibbs effect false anomalies is obtained. Then, taking this depth as the top interface of the field source, the field source morphology is constructed by dividing the grid prisms at the top interface, and density inversion is performed on the top surface of each prism in the frequency domain, and thus the density distribution of the field source morphology can be obtained. The present invention performs Fourier transform on the field source gravity data to obtain its frequency domain signal, that is, the algorithm is all implemented in the frequency domain, and the inversion is performed at the known field source depth, so the present invention is a fast, high-precision, high-resolution characteristic gravity density inversion method.

[0064] The present invention adopts a high-precision density inversion method for regular and stable rapid solution of gravity anomalies, obtains high-precision anomaly imaging data at the field source depth without Gibbs effect false anomalies, and because the algorithm is all implemented in the frequency domain and the inversion is performed at the known field source depth, the calculation speed of this method is very fast, effectively enhancing the stability and accuracy of the calculation result. Theoretically, as long as the depth determined by the regularization method is accurate, relatively high-precision true physical property parameters of the gravity field can be obtained, providing a new idea for realizing the solution of the density parameters of the gravity field source. Description of the Drawings

[0065] To more clearly illustrate the technical solution of the present invention, the accompanying drawings required in the embodiments are briefly introduced below. Obviously, the accompanying drawings in the following description are only some embodiments of the present invention. For those of ordinary skill in the art, without creative efforts, other accompanying drawings can be obtained based on these drawings.

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

[0067] Figure 2 are the specific parameters of the single sphere model in the embodiments of the present invention;

[0068] Figure 3 is the schematic diagram of the model spatial position in the embodiments of the present invention;

[0069] Figure 4 is the schematic diagram of the least squares optimization iterative solution space in the embodiments of the present invention;

[0070] Figure 5 is the flowchart of the least squares optimization iterative solution in the embodiments of the present invention;

[0071] Figure 6 is the gravity anomaly imaging map of the sphere model in the embodiments of the present invention;

[0072] Figure 7 is the main profile curve and the regularized stable solution section diagram of its Laplace lower half space in the embodiments of the present invention;

[0073] Figure 8 are the gravity anomaly imaging map (left) of the top interface of the field source center depth and the apparent density inversion result map (right) in the embodiments of the present invention. Detailed implementation manners

[0074] Next, the technical solutions in the embodiments of the present invention will be clearly and completely described in conjunction with the accompanying drawings in the embodiments of the present invention. Obviously, the described embodiments are only some embodiments of the present invention, rather than all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those of ordinary skill in the art without creative efforts belong to 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 have the ordinary meanings understood by those with ordinary skills in the field to which the present disclosure belongs. The "first", "second" and similar terms used in the embodiments of the present disclosure do not denote any order, quantity or importance, but are only used to distinguish different components. Words such as "including" or "comprising" 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. Words such as "connected" or "linked" are not limited to physical or mechanical connections, but may include electrical connections, whether direct or indirect. "Upper", "lower", "left", "right", etc. are only used to represent relative positional relationships, and when the absolute position of the object being described changes, the relative positional relationship may also change accordingly.

[0076] To make the above objects, features and advantages of the present invention more obvious and understandable, the present invention will be further described in detail below with reference to the accompanying drawings and specific embodiments.

[0077] Embodiment 1

[0078] As Figure 1 shown, the embodiments of the present invention disclose a high-precision density inversion method for regular stable and fast solution of gravity anomalies, and the method includes:

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

[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. To overcome the oscillation effect caused by the 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 regularization factor contains Gibbs effect false anomalies. To eliminate this false anomaly, the data containing false anomalies at the center of the field source depth in the lower half-space is used as the starting boundary value condition in the upper half-space. Through the solution in the upper half-space, the ground gravity anomaly data based on this boundary value condition is calculated, and the least squares optimization iteration is solved with the actual ground gravity data. After multiple fittings of the above two data, the gravity anomaly information in the lower half-space without the Gibbs effect at the center depth of the field source is obtained;

[0082] S4. Finally, taking the depth of the field source center as the top interface of the field source, constructing the field source morphology by dividing the grid prisms from the top interface, and performing density inversion on the top surface of each prism in the frequency domain, from which the density distribution at the gravity field source can be quickly obtained;

[0083] S5. Finally, perform physical property parameter constrained inversion of the single-model gravity anomaly.

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

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

[0086]

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

[0088]

[0089] Let D be the number of measurement points, and Δr be the spacing between measurement points. There is:

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

[0091]

[0092] where r is the multiple of the measurement point spacing and increases as the series order increases.

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

[0094]

[0095] where is the gravity field value at the measurement point (x, y, w), is the gravity field value observed on the ground; w is the elevation of the observation point, which is the data observed on the ground at this time, i.e., w0 = 0; n represents the discrete series order taken from 0 to N; is the frequency; μ n is the coefficient; e represents the natural constant, is the gravity field value at the measurement 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, i.e., The larger the value is, this will cause the value of the gravitational field to be infinite, resulting in the divergence of the obtained gravity anomaly data. Secondly, when solving the lower half-space, the greater the depth, that is, the larger the w value, will also make the field value infinite, leading to data divergence. That is The exponential part of

[0097] In this embodiment, in S2, in the lower half-space, a stable Fourier series with a regularization factor is selected to describe the gravity solution in the lower half-space. The center of the gravity anomaly corresponding to the center of the field source includes:

[0098] Add a regularization factor term It is a function related to the frequency and the depth (w) and coefficient (μ) in the space from below the observation surface to above the field source, etc. At this time, and w as the denominator part of the regularization factor can effectively suppress the information with the above-mentioned increasing exponent, and can be used for effective correction of the frequency components of the series terms of the potential field in the lower half-space;

[0099] Therefore, the stable Fourier series corrected by the regularization factor is:

[0100]

[0101] where η 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 it ensure that the potential field is not singular when passing through the field source during the solution process of the lower half-space, and obtain the complete field value distribution in the lower half-space, but also the Gibbs effect can not be particularly significant.

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

[0104] Taking the gravity anomaly at the depth of the field source as the boundary condition, using the calculation method of the upper half-space, calculating the ground gravity anomaly data, and carrying out the least squares optimization iteration solution with the calculated ground gravity data and the measured ground gravity data. Through continuous fitting of the data, the gravity data at the field source can be corrected, and high-precision anomaly imaging data without the Gibbs effect false anomaly at the depth of the field source can be obtained, as shown in Figure 4 and Figure 5 of the specification drawings, and its specific steps are as follows:

[0105] S31. Suppose there are two planes ΨO and ΨP, representing the earth's surface and a plane at a certain depth underground respectively. Place the anomaly value of the k-th point in the anomaly values on the ground ΨO at the vertical projection point k of the plane ΨP as the anomaly value on the plane ΨP, denoted as That is: Obtain the spectrum of the anomaly value on ΨP through Fourier transform

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

[0107]

[0108] Among them, represents the spectrum expression solved in the frequency domain.

[0109] Substitute the spectrum of the anomaly value on the plane ΨP into to find the anomaly value of the k-th point 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; represents the Fourier transform of the gravitational anomaly value on the ground ΨO (w = 0); F -1 represents the inverse Fourier transform; d represents the solution height of the upper half space;

[0112] S33. Use the difference between and to correct the anomaly value on the plane ΨP to obtain a new anomaly value

[0113]

[0114] Among them, q represents the step size, generally taken as 1;

[0115] S34. Repeat S32 and S33 to obtain the least squares optimization iteration solution formula:

[0116]

[0117] Among them, 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, then 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 field source center in the lower half space to be solved

[0122] In this embodiment, in S4, taking the depth of the field source center 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 of 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 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. Obtain the spectrum of the gravity anomaly generated by each prismatic body at the known top interface (x d , y d , 0) of the field source center depth in the frequency domain Formula:

[0125]

[0126] According to the superposition property, obtain the forward 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 dThe spectrum of the gravity anomaly G(x, y, 0) generated at (x, y, 0)

[0127]

[0128] where s and 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 at the corresponding top interface of this depth layer is known, then d = 0. If d ≠ 0, that is, the buried depth is not zero, the gravity anomaly data of the top interface of this depth layer is calculated by least squares optimization iteration in the lower half space, and then according to

[0130] Let:

[0131]

[0132] From

[0133] and we get:

[0134]

[0135] According to it is deduced that:

[0136]

[0137] For performing Fourier inverse transform, the residual density values of each point (x d , y d , 0) of the top interface of this depth layer are obtained, that is, the inversion formula for calculating relevant physical property parameters from the known field source gravity anomaly:

[0138]

[0139] According to the inversion formula for calculating relevant physical property parameters from the known field source gravity anomaly, the apparent density values corresponding to the gravity anomaly are calculated. Since these data all contain the corresponding depth, therefore, fast density inversion can be performed on them to obtain high-precision anomaly imaging data of the field 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, with the point distance and line distance both being 1 (unit: meter), that is, the size of the survey area is 100×100m 2 , and the specific parameters of a single sphere model are as Figure 2As shown, the spatial position of the model is as Figure 3 shown.

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

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

[0144] Taking the obtained field source center depth as the top interface of the field source, the field source morphology is constructed by dividing the grid prisms at the top interface. The density inversion is performed on the top surface of each prism in the frequency domain. From this, the density distribution of the field source morphology can be obtained. It is easy to see from the apparent density inversion result diagram that it is basically consistent with the density of the assumed single sphere model. Figure 8 is the gravity anomaly (left) and the apparent density inversion result (right) at the top interface of the field source center depth. Since the algorithm is implemented in the frequency domain and the inversion is performed at the known field source depth, the present invention is a fast, high-precision, and high-resolution gravity density inversion method.

[0145] To sum up: A high-precision density inversion method for gravity anomaly based on regular and stable fast solution provided by the present invention has the following advantages compared with the traditional transformation method:

[0146] 1) The accuracy of the gravity anomaly parameter inversion is related to the accuracy of the constrained depth. Complex anomalies will cause the accuracy of the constrained depth to decrease to some extent. However, since the ground gravity anomaly is used as the first boundary value condition for the regularized stable solution of the Laplace lower half space, and a function carrying a regularization factor is adopted to correct the series terms of the gravity field value solved in the lower half space by frequency, the divergence of the result is effectively avoided. The high-resolution information near the top interface obtained in this way can still ensure that the physical property parameters have a very high lateral resolution;

[0147] 2) Based on the above solution of field source regularization, it is proposed to use the data carrying false anomalies at the imaging depth of the field source as the boundary value condition. Then, by solving the upper half space of the first boundary value problem of the Laplace equation again, the gravity anomaly at the ground is calculated, and the least squares optimization iteration is carried out with the actual ground anomaly to obtain high-precision anomaly imaging data at the field source depth without Gibbs effect false anomalies;

[0148] 3) This algorithm is implemented in the frequency domain, and the inversion is carried out at the known field source depth. The calculation speed of this method is very fast, effectively enhancing the stability and accuracy of the calculation results. Theoretically, as long as the center depth of the field source determined by the regularized stable solution of the lower half space of the Laplace equation is accurate, relatively high-precision true physical property parameters of the gravity field can be obtained;

[0149] 4) The high-precision density inversion method based on regular and stable fast solution of gravity anomalies is conducive to improving the inversion accuracy of physical property parameters with depth constraints for complex geological models, can effectively improve the effectiveness and accuracy of deep mineral resources exploration, increase the prospecting success rate, reduce the exploration cost, and provide important technical support for resource development in our country.

[0150] It should be noted that the method of the embodiments of the present disclosure can be executed by a single device, such as a computer or a server. The method of this embodiment can also be applied to a distributed scenario and completed by multiple devices cooperating with each other. In this case of a distributed scenario, one of the multiple devices can only execute one or more steps of the method of the embodiments of the present disclosure, and these 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 appended claims. In some cases, it should be understood that the size of the sequence numbers of the steps in the above embodiments does not mean the order of execution. The execution order of each process should be determined by its function and internal logic, and should not constitute any limitation to the implementation process of the embodiments of the present invention. The actions or steps recorded in the claims can be executed in a different order from that in the above embodiments and still achieve the desired results. Additionally, the processes depicted in the drawings do not necessarily require the specific order or continuous order shown to achieve the desired results. In certain embodiments, multitasking and parallel processing are also possible or may be advantageous.

[0152] Embodiment 2

[0153] Based on the same inventive concept, corresponding to the method of any of the above embodiments, the present disclosure also provides a high-precision density inversion system based on regular and stable fast solution of gravity anomalies. The system is used to implement the method of any one of the above, and includes: 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 gravitational field and the corresponding boundary conditions, i.e., the gravitational data observed on the ground.

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

[0156] The iteration module is used to take the data containing false anomalies at the center depth of the field source in the lower half-space as the starting boundary value condition in the upper half-space. Through the solution in the upper half-space, the ground gravitational anomaly data based on this boundary value condition is calculated, and the least squares optimization iteration is performed with the actual ground gravitational data to obtain the gravitational anomaly information in the lower half-space without Gibbs effect at the center depth of the field source.

[0157] The inversion module is used to take the center depth of the field source as the top interface of the field source, construct the field source morphology from the top interface meshed prismatic bodies, and perform density inversion on the top surface of each prismatic body in the frequency domain to obtain the density distribution at the gravitational field source.

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

[0159]

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

[0161] In this embodiment, in the second calculation module, selecting a stable Fourier series with a regularization factor in the lower half-space to describe the gravitational solution in the lower half-space, and the center of the gravitational anomaly corresponds to the center of the field source includes:

[0162]

[0163] Wherein, η n is the regularization factor.

[0164] The system of the above embodiment is used to implement the corresponding high-precision density inversion method for stably and quickly solving gravitational anomalies in any one of the foregoing embodiments, and has the beneficial effects of the corresponding method embodiments, which will not be elaborated here.

[0165] It should be noted that the above high-precision density inversion system for stably and rapidly solving gravity anomalies based on regular expressions is embodied in the form of functional units. The term "module" here can be implemented in the form of software and / or hardware, and no specific limitation is made thereto.

[0166] For example, the "module" can be a software program, a hardware circuit, or a combination of both that implements the above functions. The hardware circuit may include an application specific integrated circuit (ASIC), an electronic circuit, a processor for executing one or more software or firmware programs (such as a shared processor, a proprietary processor, or a group of processors, etc.), and a memory, a merging 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 principle of the embodiments of the present disclosure shall be included within the protection scope of the present disclosure.

Claims

1. A high-precision density inversion method for stably and rapidly solving gravity anomalies based on regularization, characterized in that, The method includes: S1. Obtain the 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; S2. In the lower half-space, select the stable Fourier series with a regularization factor 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; S3. Use the data containing false anomalies at the center of the field source depth in the lower half-space as the starting boundary value condition for the upper half-space. Through the solution of the upper half-space, calculate the ground gravity anomaly data based on this boundary value condition, and perform least squares optimization iterative solution with the actual ground gravity data to obtain the gravity anomaly information in the lower half-space without Gibbs effect at the center depth of the field source; S4. Take the center depth of the field source as the top interface of the field source, construct the field source morphology by dividing the grid prisms at the top interface, and perform density inversion on the top surface of each prism in the frequency domain to obtain the density distribution at the gravity field source; In the above S1, obtaining the 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: Among them, is the gravity field value at the measurement point (x, y, w), where w is the elevation of the observation point; n represents the discrete series order taken from 0 to N; ζ n is the frequency; μ n is the coefficient; e represents the natural constant; In the above S2, selecting the stable Fourier series with a regularization factor 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 includes: Among them, η n is the regularization factor, r is the multiple of the measuring point spacing, and it increases as the series order increases. is the gravity field value at the measuring point (r, w), and μ0 is the coefficient.

2. The method according to claim 1, characterized in that, In the above S3, using the data containing false anomalies at the center depth of the field source in the lower half-space as the starting boundary value condition for the upper half-space. Through the solution of the upper half-space, calculate the ground gravity anomaly data based on this boundary value condition, and perform least squares optimization iterative solution with the actual ground gravity data to obtain the gravity anomaly information in the lower half-space without Gibbs effect at the center depth of the field source includes: S31. Assume there are two planes ΨO and ΨP, representing the surface and a plane at a certain depth underground respectively. Place the anomaly value of the k-th point in the anomaly values on the ground ΨO at the vertical projection point k on the plane ΨP as the anomaly value on the plane ΨP, denoted as That is: Obtain the spectrum of the anomaly value on ΨP through Fourier transform S32. When there is no field source distribution between 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: Substitute the outliers on the plane ΨP into the spectrum of Find the anomaly value at point k on the ground ΨO If w = d, that is: where s and t are spatial frequencies, representing the wave numbers of x and y respectively; is expressed as the gravity anomaly value on the ground ΨO(w = 0); is the Fourier transform of; F -1 represents the inverse Fourier transform; d represents the solution height in the upper half space; S33. Use the and difference to correct the outlier on the plane ΨP to obtain a new outlier where q represents the step size, and generally it is taken as 1; S34. Repeat S32 and S33 to obtain the least squares optimization iterative solution formula: When it is, from obtained wherein, δ represents a given number approaching zero; if meets the design accuracy, i.e., when u < U, the iteration is stopped, 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 in the frequency domain for the first boundary value problem of the Laplace equation in the upper half space: Obtain by performing the inverse Fourier transform The gravity anomaly at the depth d of the field source center in the lower half-space to be solved 3. The method according to claim 2, wherein In the above S4, taking the center depth of the field source as the top interface of the field source, constructing the field source morphology by dividing the grid prisms at the top interface, and performing density inversion on the top surface of each prism in the frequency domain to obtain the density distribution at the gravity field source includes: will cause gravity anomalies The irregular geological body is divided into AB vertically small prisms with the same size and different centroid positions, where AB = A × B; Let the centroid coordinates of the small prism be (x 0i , y 0i ), the remaining density be θ i , i = 1, 2, 3, …, AB, with side lengths a and b respectively, the buried depth of the top surface of the prism be d, and the height be cd. Obtain the spectrum of the gravity anomaly generated by each prism at the top interface (x d , y d , 0) at the known depth of the field source center in the frequency domain Formula: According to the superposition principle, the forward formula of the field-source gravity anomaly with known physical property parameters is obtained, that is, the spectrum of the gravity anomaly G(x, y, 0) generated by the entire depth layer at the top interface of the known field-source center depth (x d , y d , 0) where s and 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 this depth layer at the corresponding top interface is known, then d = 0. If d ≠ 0, that is, the buried depth is not zero, calculate the gravity anomaly data of the top interface of this depth layer through the least squares optimization iterative solution in the lower half-space, and then according to Let: Derived from and we get: According to It is deduced that there is: Pair Perform the inverse Fourier transform to obtain the residual density values of each point (x d , y d , 0) on the top interface of this depth layer, which is the inversion formula for calculating relevant physical property parameters from the known field source gravity anomaly: According to the inversion formula for calculating relevant physical property parameters from the known gravity anomaly of the field source, calculate the apparent density value corresponding to the gravity anomaly.

4. A high-precision density inversion system for stably and rapidly solving gravity anomalies based on regularization, which is used to implement the method described in any one of claims 1-3, and is characterized in that, It includes: The first calculation module, the second calculation module, the iteration module and the inversion module; The first calculation module is used to obtain the 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 with a regularization factor in the lower half-space to describe the gravity solution in the lower half-space, where the center of the gravity anomaly corresponds to the center of the field source; The iteration module is used to take the data containing false anomalies at the center depth of the field source in the lower half-space as the starting boundary value condition for the upper half-space. Through the solution of the upper half-space, the ground gravity anomaly data based on this boundary value condition is calculated, and the least squares optimization iteration is performed with the actual ground gravity data to obtain the gravity anomaly information in the lower half-space without Gibbs effect at the center depth of the field source; The inversion module is used to take the center depth of the field source as the top interface of the field source, construct the field source morphology from the top interface dissected grid prisms, and perform density inversion on the top surface of each prism in the frequency domain to obtain the density distribution at the gravity field source.

5. The system according to claim 4, wherein In the first calculation module, through the Laplace equation of the gravity field and the corresponding boundary conditions, i.e., the gravity data observed on the ground, the stable Fourier series solution under general conditions includes: Among them, is the gravity field value at the measurement point (x, y, w), where w is the elevation of the observation point; n represents the discrete series order taken from 0 to N; ζ n is the frequency; μ n is the coefficient; e represents the natural constant.

6. The system according to claim 5, characterized in that In the second calculation module, in the lower half-space, selecting a stable Fourier series with a regularization factor to describe the gravity solution in the lower half-space, where the center of the gravity anomaly corresponds to the center of the field source includes: where η n is the regularization factor, r is the multiple of the measuring point spacing, and increases as the series order increases. is the gravity field value at the measuring point (r, w), and μ0 is the coefficient.

Images