High-precision density inversion method and system for gravity anomalies based on regular, stable, and fast solution

Abstract

The present invention discloses a method and system for high-precision density inversion of gravity anomalies based on regular stable fast solution. [Solution] Based on the above source regularization solution, we propose using data containing spurious anomalies at the source imaging depth as boundary value conditions. Then, through the upper half-space solution of the first-class boundary value problem of Laplace's equation, we calculate the gravity anomaly on the ground, and perform a least-squares optimization iteration with the actual ground anomaly to obtain high-precision anomaly imaging data at the source depth without spurious anomalies caused by the Gibbs effect. This depth is used as the top interface of the source, and the top interface is divided into grid prisms to construct a source configuration. Density inversion is performed on each prism top surface in the frequency domain to obtain the density distribution of the source configuration. Because this algorithm is implemented in the frequency domain and the inversion is performed at a known source depth, the present invention provides a gravity density inversion method with high speed, high precision, and high resolution, and provides a new concept for solving gravity source density parameters.

Description

[Technical Field]

[0001] The present invention relates to the technical field of mineral exploration, and in particular to a method and system for high-precision density inversion of gravity anomalies based on regular, stable, and fast solution. [Background technology]

[0002] Gravity surveying is a common geophysical exploration technique with wide applications in a variety of fields, including energy and mineral resource exploration, geological mapping, and engineering construction. Gravity anomaly data processing and transformation are an important part of gravity data interpretation theory. Solving the Laplace equation in the upper half-space (region away from the field source) with surface gravity anomalies as the first-class boundary value condition allows for a regular, stable solution. Solving the Laplace equation in the lower half-space (region approaching the field source) provides higher-resolution gravity anomaly information, but the lower half-space is subject to the boundary value constraints of the Laplace equation, making it impossible to obtain gravity anomaly information at the field source. In the process of solving the lower half-space, the potential field solved at the field source is not singular. Therefore, to obtain the complete field distribution in the lower half-space, it is necessary to optimize and improve the method for solving the gravity field in the lower half-space.

[0003] The distribution characteristics of the gravitational and magnetic potential fields in the lower half-space are highly intuitive and provide an important basis for studying the field source. Using ground-based gravity anomalies as the first-class boundary condition, a regularized stable solution of the Laplace lower half-space can provide high-resolution anomaly data at the source depth. However, due to the Gibbs cutoff effect of the frequency domain solution, weak anomaly zones with alternating positive and negative "eight" shapes exist on both sides of the main anomaly zone, centered on the distribution area of ​​the main source. These anomalies can be isolated or overlap with the main anomalies of other field sources to form compound anomalies, making the interpretation of the field source more difficult. Summary of the Invention [Problem to be solved by the invention]

[0004] To solve the above technical problems, the present invention provides a method and system for high-precision density inversion of gravity anomalies based on regularized, stable, and fast solution. Based on the source regularization solution, data containing spurious anomalies at the source imaging depth are used as boundary value conditions. The gravity anomaly on the ground is calculated through the upper half-space solution of the first-class boundary value problem of Laplace's equation. Then, using the actual ground anomaly and least-squares optimization iterative solution, high-precision anomaly imaging data at the source depth, free of spurious anomalies due to the Gibbs effect, is obtained. Using this depth as the top interface of the source, the top interface is divided into grid prisms to construct a source configuration. Density inversion is performed on each prism top surface in the frequency domain to obtain the density distribution of the source configuration. Since this algorithm is implemented entirely in the frequency domain and the inversion is performed at a known source depth, the present invention provides a high-speed, high-precision, and high-resolution gravity density inversion method, providing a new concept for comprehensively interpreting field sources from both physical properties and geometric parameters. [Means for solving the problem]

[0005] The present invention discloses a high-precision density inversion method for gravity anomaly based on regular stable fast solution, the method comprising: S1: obtaining a stable Fourier series solution under general conditions through the Laplace equation of gravity field and corresponding boundary conditions, i.e., gravity data from ground observations; S2, in the lower half space, select a stable Fourier series with regular factors to describe the gravity solution in the lower half space, and make the center of gravity anomaly correspond to the center of the field source; S3: Using the data containing the false anomaly at the center of the field source depth in the lower half space as the initial boundary value condition of the upper half space, calculating the ground gravity anomaly data based on the boundary value condition through the solution of the upper half space, and performing the least squares optimization iteration solution with the actual ground gravity data to obtain the lower half space gravity anomaly information that does not include the Gibbs effect at the center of the field source depth; S4: The field source central depth is set as the top interface of the field source, the top interface is divided into grid prisms to construct a field source shape, and density inversion is performed on the top surface of each prism in the frequency domain to obtain the density distribution at the gravitational field source.

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

number

number

number

[0007] Preferably, in S2, the step of selecting a stable Fourier series with a regular factor in the lower half space to describe the gravity solution in the lower half space and making the center of the gravity anomaly correspond to the center of the field source includes the following equation:

number

number

[0008] Preferably, in S3, the step of using the data including the false anomaly at the center of the field source depth in the lower half space as the initial boundary value condition of the upper half space, calculating the ground gravity anomaly data based on the boundary value condition by solving the upper half space, and performing least squares optimization iteration with the actual ground gravity data to obtain the lower half space gravity anomaly information without the Gibbs effect at the center of the field source depth is as follows: S31, Suppose there are two planes ΨO and ΨP, which represent the planes at a certain depth above and below ground, respectively.

number

number

number

number

number

number

number

number

number

number

number

number

number

number

number

number

number

number

number

number

number

number

Number

Number

Number

Number

Number

[0009] Preferably, in the above S4, taking the center depth of the field source as the top surface of the field source, dividing the top surface into grid prisms to construct the field source form, and performing density inversion on the top surface of each prism in the frequency domain to obtain the density distribution in the gravity field source, the steps are: Gravity anomaly

Number

number

number

number

number

number

number

number

number

number

number

number

number

number

number

[0010] The present invention further provides a high-precision density inversion system for gravity anomaly based on regular stable fast solution, which is used to realize any one of the above methods, and includes 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 gravity field and corresponding boundary conditions, i.e., gravity data of ground observation; The second calculation module is used to select a stable Fourier series with regular factors in the lower half space to describe the gravity solution in the lower half space, and make the center of the gravity anomaly correspond to the field source center; The iterative module uses the data containing the false anomaly at the center of the field source depth in the lower half space as the initial boundary value condition of the upper half space, calculates the ground gravity anomaly data based on the boundary value condition through the solution of the upper half space, and performs least squares optimization iterative solution with the actual ground gravity data to obtain the gravity anomaly information of the lower half space that does not contain the Gibbs effect at the center of the field source depth; The inversion module uses the central depth of the field source as the top interface of the field source, divides the top interface into grid prisms to construct a field source shape, and performs density inversion on each prism top surface in the frequency domain to obtain the density distribution at the gravitational field source.

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

number

number

number

[0012] Preferably, in the second calculation module, the step of selecting a stable Fourier series with regular factors in the lower half space to describe the gravity solution in the lower half space and making the center of gravity anomaly correspond to the field source center includes the following equation:

number

number

[0013] Compared with the prior art, the beneficial effects of the present invention are as follows:

[0014] This invention provides a high-precision density inversion method for gravity anomalies based on regularized stable fast solution. Compared with traditional inversion methods, this invention first uses ground gravity anomalies as the first-class boundary value condition to perform regularized stable solution of the Laplace lower half-space, and then uses a function with a regularization factor to correct the series term frequency of the gravity field value solved in the lower half-space, effectively avoiding divergence of the results and obtaining high-resolution anomaly data around the source depth. Based on the above field regularized solution, it proposes using data with spurious anomalies at the source imaging depth as the boundary value condition, and then calculates the gravity anomaly on the ground through solving the first-class boundary value problem of the Laplace equation in the upper half-space. Then, it performs least-squares optimization iteration with the actual ground anomaly to obtain high-precision anomaly imaging data at the source depth without spurious anomalies caused by the Gibbs effect.

[0015] To separate anomalies generated on the surface by sources at different depths underground, it is often necessary to solve the lower half-space up to the top of different depth layers. The stability and depth of the solution directly affect the quality of the final results. Research has shown that to address the instability of potential field solutions in the lower half-space and the depth limitations of traditional FFT methods for solving the lower half-space, a least-squares optimization iterative solution method in the spatial domain has been adopted. This method is based on a stable FFT upper half-space solution, is simple in principle, and is stable downward. Assuming there are two planes, ΨO and ΨP, representing the surface and underground planes at a certain depth, and there is no field source distribution between them, i.e., a source-free space, the least-squares optimization iterative solution flowchart can be used to iteratively solve the lower half-space based on the known depth of the upper half-space solution, thereby obtaining high-precision anomaly imaging data at the source depth without false anomalies caused by the Gibbs effect. Furthermore, this depth is used as the top interface of the field source, and the top interface is divided into grid prisms to construct a field source configuration, and density inversion is performed on each prism top surface in the frequency domain, thereby obtaining the density distribution of the field source configuration.In this invention, the gravity data of the field source is Fourier transformed to obtain its frequency domain signal, that is, this algorithm is all realized in the frequency domain, and the inversion is performed at a known field source depth, so the present invention is a gravity density inversion method with high speed, high accuracy, and high resolution characteristics.

[0016] The present invention employs a high-precision density inversion method for gravity anomalies based on regularization, stable, and fast solution, to obtain high-precision anomaly imaging data of the field source depth without false anomalies caused by the Gibbs effect. Furthermore, since the algorithm is fully implemented in the frequency domain and the inversion is performed at a known field source depth, the calculation speed of this method is very fast and the stability and accuracy of the calculation results are effectively improved. Theoretically, if the depth determined by the regularization method is accurate, it is possible to obtain high-precision true physical property parameters of the gravity field, providing a new concept for solving the gravity field source density parameters. [Brief explanation of the drawings]

[0017] In order to more clearly explain the technical solutions of the present invention, the following briefly describes the drawings that need to be used in the embodiments. Obviously, the drawings described below are only some embodiments of the present invention, and those skilled in the art can obtain other drawings based on these drawings without any creative efforts. [Figure 1] 1 is a technology roadmap of the present invention. [Figure 2] FIG. 10 is a diagram showing specific parameters of a single sphere model according to an embodiment of the present invention. [Figure 3] 1 is a schematic diagram of a model space position of an embodiment of the present invention; [Figure 4] FIG. 2 is a spatial schematic diagram of a least squares optimization iterative solution according to an embodiment of the present invention; [Figure 5] 1 is a flowchart of a least squares optimization iterative solution according to an embodiment of the present invention. [Figure 6] FIG. 1 is a diagram showing gravity anomaly imaging of a spherical model according to an embodiment of the present invention. [Figure 7] 1 is a cross-sectional view of a principal cross section curve and its Laplace lower half-space regularized stable solution according to an embodiment of the present invention; [Figure 8] This is a diagram of gravity anomaly imaging at the top interface of the field source center depth (left) and a diagram of apparent density inversion results (right) of an embodiment of the present invention. DETAILED DESCRIPTION OF THE INVENTION

[0018] The technical solutions in the embodiments of the present invention will be clearly and completely described below in conjunction with the drawings in the embodiments of the present invention. Obviously, the described embodiments are only some embodiments of the present invention, and are not all embodiments. Based on the embodiments of the present invention, all other embodiments that can be obtained by those skilled in the art without creative labor shall fall within the protection scope of the present invention.

[0019] Unless otherwise defined, technical or scientific terms used in the embodiments of the present disclosure should have the ordinary meaning understood by a person of ordinary skill in the field to which the present disclosure belongs. The terms "first," "second," and similar terms used in the embodiments of the present disclosure do not indicate any order, quantity, or importance, but are merely used to distinguish between different components. Similar terms such as "comprise" or "include" mean that the element or object appearing before the term covers the elements or objects listed thereafter and their equivalents, and do not exclude other elements or objects. Similar terms such as "connection" or "connected" are not limited to physical or mechanical connections, but can also include electrical connections, whether direct or indirect. Terms such as "upper," "lower," "left," and "right" only indicate relative positional relationships, and if the absolute position of the described object changes, these relative positional relationships may also change accordingly.

[0020] In order to make the above objects, features and advantages of the present invention clearer and easier to understand, the present invention will be described in more detail below with reference to the drawings and specific embodiments.

[0021] Example 1

[0022] As shown in FIG. 1, an embodiment of the present invention provides a high-precision density inversion method for gravity anomaly based on regular stable fast solution, the method comprising: S1: First, perform a Fourier series transformation of the observed gravity anomaly data, and then calculate the Laplace equation of the gravity field and its boundary conditions (gravity data from ground observations).

number

number

[0023] In this embodiment, in the step S1, the step of 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 of the ground observation, is as follows: Laplace's equation of gravity and its boundary conditions (gravity data from ground observations) include:

number

number

number

number

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

number

number

number

number

number

[0025] There is a relationship between frequency and depth of the lower half-space solution, so shallower information is higher in the frequency domain, i.e.

number

number

[0026] In this embodiment, in S2, the step of selecting a stable Fourier series with regular factors in the lower half space to describe the gravity solution in the lower half space and making the center of the gravity anomaly correspond to the field source center is as follows: regularization factor term

number

number

number

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

number

[0028] This allows us to obtain the gravity anomaly at the depth of the field source. In practice, by adding a regularization factor, the potential field passing through the field source in the lower half-space solution process is not singular, and not only can we obtain the complete field value distribution in the lower half-space, but the Gibbs effect is also not particularly pronounced.

[0029] In this embodiment, in S3, the data including the false anomaly at the center of the field source depth in the lower half space is used as the initial boundary value condition of the upper half space, and the ground gravity anomaly data based on the boundary value condition is calculated by solving the upper half space, and the actual ground gravity data is used for least squares optimization iteration to obtain the lower half space gravity anomaly information without the Gibbs effect at the center of the field source depth. Taking the gravity anomaly at the source depth as the boundary condition, the upper half-space calculation method is used to calculate the ground gravity anomaly data, and the calculated ground gravity data and the ground measured gravity data are subjected to least squares optimization iteration to solve the data, and the gravity data at the source is corrected through continuous data fitting, so that high-precision anomaly image data at the source depth can be obtained that does not contain false anomalies caused by the Gibbs effect. Specifically, as shown in Figures 4 and 5 of the instruction manual, S31, Suppose there are two planes ΨO and ΨP, which represent the planes at a certain depth above and below ground, respectively.

number

number

number

number

number

number

number

number

number

number

number

number

number

number

number

number

number

Number

Number

Number

Number

Number

Number

Number

Number

Number

Number

Number

number

number

number

[0030] In this embodiment, in S4, the field source center depth is the top interface of the field source, the top interface is divided into grid prisms to construct the field source form, and the density inversion of each prism top surface in the frequency domain is performed to obtain the density distribution of the gravitational field source.

[0031] gravity anomaly

number

number

number

number

number

number

number

number

number

number

number

number

number

number

number

[0032] In this embodiment, in step S5, the physical property parameter constraint inversion of the single model gravity anomaly is carried out.

[0033] The present invention designs a 101x101 grid, and the point interval and line interval are all 1 (unit: meter), that is, the size of the measurement range is 100x100m 2 The specific parameters of the single sphere model are shown in Figure 2, and the spatial position of the model is shown in Figure 3.

[0034] Based on the forward calculation of the above model, the corresponding gravity anomaly is obtained, as shown in Figure 6. Figure 7 shows the main cross-sectional curve of the gravity anomaly and cross-sectional views of each depth layer obtained by the regularized stable solution of the Laplace lower half-space. It can be seen that the center of the anomaly coincides with the specified model depth. However, due to the Gibbs cutoff effect of the frequency domain solution, there is an "eight"-shaped, alternating positive and negative weak anomaly band centered on the distribution region of the main source body. This reduces the lateral resolution of the field source imaging and increases the difficulty of utilizing the spatial distribution characteristics of the gravitational potential field.

[0035] Therefore, based on the above-mentioned source regularization solution, the present invention proposes to use the data with spurious anomalies at the source imaging depth as the boundary value condition, and again calculate the gravity anomaly on the ground through solving the upper half-space of the first-class boundary value problem of Laplace's equation, and then perform least-squares optimization iteration with the actual ground anomaly to obtain high-precision anomaly imaging data at the source depth that does not contain spurious anomalies caused by the Gibbs effect. The spatial conceptual diagram of the processing is shown in Figure 4, and the specific flow is shown in Figure 5.

[0036] The obtained field source central depth is used as the top interface of the field source, and the top interface is divided into grid prisms to construct the field source shape. Density inversion is then performed on each prism top surface in the frequency domain, thereby obtaining the density distribution of the field source shape. From the apparent density inversion results, it is easy to see that they are nearly identical to the assumed single spherical model density. Figure 8 shows the gravity anomaly (left) and apparent density inversion results (right) for the top interface of the field source central depth. Since this algorithm is entirely implemented in the frequency domain and inversion is performed at a known field source depth, the present invention provides a gravity density inversion method with high speed, high accuracy, and high resolution.

[0037] In summary, the high-precision density inversion method for gravity anomalies based on regular stable fast solution provided by the present invention has the following advantages over conventional conversion methods:

[0038] 1) The accuracy of gravity anomaly parameter inversion is related to the accuracy of the constraint depth, and complex anomalies may reduce the accuracy of the constraint depth. However, by using the ground gravity anomaly as the first-class boundary value condition to perform a regularized stable solution of the Laplace lower half-space and correcting the series term frequency of the gravity field value in the lower half-space solution using a function with a regularization factor, divergence of the results can be effectively avoided, and the high-resolution information obtained near the apex interface still ensures very high lateral resolution of the physical property parameters.

[0039] 2) Based on the above-mentioned field source regularization solution, we propose to use data with spurious anomalies at the field source imaging depth as the boundary value condition, and again calculate the gravity anomaly on the ground through solving the upper half-space of the first-class boundary value problem of Laplace's equation, and then perform a least-squares optimization iterative solution with the actual ground anomaly, thereby obtaining high-precision anomaly imaging data at the field source depth that does not contain spurious anomalies caused by the Gibbs effect.

[0040] 3) This algorithm is implemented entirely in the frequency domain, and the inversion is performed at a known field source depth. This method has a very fast calculation speed and effectively improves the stability and accuracy of the calculation results. In theory, if the field source center depth determined by the regularized stable solution of the Laplace lower half-space is accurate, it is possible to obtain highly accurate true physical parameters of the gravitational field.

[0041] 4) The high-precision density inversion method for gravity anomalies based on regular stable high-speed solution is beneficial to improving the accuracy of depth-constrained physical property parameter inversion for complex geological models, which can improve the effectiveness and accuracy of deep mineral resource exploration, increase the exploration success rate, and reduce exploration costs, providing important technical guarantees for China's resource development.

[0042] It should be noted that the method of the present disclosure may be performed by a single device, such as a computer or server, or may be implemented in a distributed scenario where multiple devices cooperate with each other to execute the method. In such a distributed scenario, only one of the multiple devices may execute one or more steps of the method of the present disclosure, and the multiple devices interact with each other to complete the method.

[0043] The above describes several exemplary embodiments of the present disclosure. Other exemplary embodiments are also within the scope of the claims. In some cases, it should be understood that the magnitude of the sequence numbers of each step in the above exemplary embodiments does not imply an execution order, and the execution order of each process is determined by its function and inherent logic, and does not impose any limitations on the implementation process of the exemplary embodiments of the present invention. The operations or steps described in the claims may be performed in an order different from that of the above exemplary embodiments and still achieve the expected results. Furthermore, the processes depicted in the figures do not necessarily require the specific order or sequential order shown, and may still achieve the expected results. In some embodiments, multitasking and parallel processing may also be possible or advantageous.

[0044] Example 2

[0045] According to the same inventive concept, and corresponding to the method of any of the above embodiments, the present invention further provides a high-precision density inversion system for gravity anomaly based on regular stable fast solution, the system is used to realize any one of the methods, and includes: 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 gravity field and corresponding boundary conditions, i.e., gravity data of ground observation; The second calculation module selects a stable Fourier series with regular factors in the lower half space, which is used to describe the gravity solution in the lower half space and make the center of the gravity anomaly correspond to the field source center; The iterative module uses the data containing the false anomaly at the center of the field source depth in the lower half space as the initial boundary value condition of the upper half space, calculates the ground gravity anomaly data based on the boundary value condition through the solution of the upper half space, and performs least squares optimization iterative solution with the actual ground gravity data to obtain the lower half space gravity anomaly information that does not contain the Gibbs effect at the center of the field source depth; The inversion module uses the central depth of the field source as the top interface of the field source, divides the top interface into grid prisms to construct the field source shape, and performs density inversion on each prism top surface in the frequency domain to obtain the density distribution at the gravitational field source.

[0046] In this embodiment, in the first calculation module, the step of 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 of the ground observation, includes the following equation:

number

number

number

[0047] In this embodiment, in the second calculation module, the step of selecting a stable Fourier series with regular factors in the lower half space to describe the gravity solution in the lower half space and making the center of gravity anomaly correspond to the field source center includes the following formula:

number

[0048] The system of the above embodiment can be used to realize a high-precision density inversion method for gravity anomaly based on the corresponding regular stable fast solution in any of the above embodiments, and has the beneficial effects of the corresponding method embodiment, so the description will be omitted here.

[0049] The gravity anomaly high-precision density inversion system based on the regular stable high-speed solution is realized in the form of a functional unit. The "module" here can be realized in software and / or hardware form, but is not limited thereto.

[0050] For example, a "module" may be a software program, a hardware circuit, or a combination thereof that implements the described functionality. The hardware circuit may include an application-specific integrated circuit (ASIC), electronic circuitry, a processor (e.g., a shared processor, dedicated processor, or group processor) and memory that executes one or more software or firmware programs, integrated logic circuitry, and / or other suitable components that support the described functionality.

[0051] The embodiments of the present disclosure are intended to cover all such substitutions, amendments and variations that fall within the broad scope of the claims, and therefore, any omissions, amendments, equivalent substitutions, improvements, etc. within the spirit and scope of the embodiments of the present disclosure shall fall within the scope of protection of the present disclosure.

Claims

1. A high-precision density inversion method for gravity anomaly based on regular stable fast solution, the method comprising: S1: obtaining a stable Fourier series solution under general conditions through the Laplace equation of the gravitational field and corresponding boundary conditions, i.e., the gravity data of ground observations; S2, in the lower half space, select a stable Fourier series with regular factors to describe the gravity solution in the lower half space, and make the center of gravity anomaly correspond to the center of the field source; S3: Using the data containing the false anomaly at the center of the field source depth in the lower half space as the initial boundary value condition of the upper half space, calculating the ground gravity anomaly data based on the boundary value condition through the solution of the upper half space, and performing least squares optimization iteration solution with the actual ground gravity data to obtain the lower half space gravity anomaly information that does not include the Gibbs effect at the center of the field source depth; S4. The field source central depth is the top interface of the field source, and the top interface is divided into grid prisms to construct a field source shape, and density inversion is performed on the top surface of each prism in the frequency domain to obtain the density distribution at the gravitational field source; In the step S1, the step of 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 of the ground observation, includes the following equation: [Equation 1] where: [Equation 2] is the gravity field value at the measurement point (x, y, w), w is the altitude of the observation point, and n is a discrete series order from 0 to N. [Equation 3] is the frequency, μ ｎ represents a coefficient, and e represents a natural constant. In the step S2, the step of selecting a stable Fourier series with regular factors in the lower half space to describe the gravity solution in the lower half space and making the center of the gravity anomaly correspond to the center of the field source includes the following equation: [Equation 4] where η ｎ is the regularization factor, r is a multiple of the measurement point interval, and increases as the series order increases. [Equation 5] is the gravitational field value at the measurement point (r, w), μ ｎ A high-precision density inversion method for gravity anomalies based on regular, stable, and fast solution, characterized by the fact that is a coefficient.

2. In the step S3, the data including the false anomaly at the center of the field source depth in the lower half space is used as the initial boundary value condition of the upper half space, and the ground gravity anomaly data based on the boundary value condition is calculated by solving the upper half space, and the actual ground gravity data is used for least squares optimization iteration to obtain the lower half space gravity anomaly information that does not include the Gibbs effect at the center of the field source depth. S31, assuming there are two planes ΨO and ΨP representing the planes at a certain depth above ground and underground, respectively, an anomaly on the surface ΨO [Equation 6] The kth abnormal value in [Equation 7] is placed at the vertical projection point k on the plane ΨP, and as an outlier on the plane ΨP, [Equation 8] Let us say, [Equation 9] , the spectrum of the outliers on the plane ΨP is obtained by Fourier transformation. [Equation 10] and S32, when there is no field source distribution between the two planes ΨO and ΨP, the gravitational field value satisfies Laplace's equation, and according to the solution formula in the frequency domain of the upper half space of the first class boundary value problem of Laplace's equation, [0011] Outliers on the plane ΨP [0012] Spectrum of [0013] of [0014] Substitute into Anomalies of k-points on the ground ΨO [Equation 15] When w=d, the following equation is obtained: [0016] where s and t are spatial frequencies, and represent x and y wave numbers, respectively. [Equation 17] is the gravity anomaly value on the ground ΨO (w = 0) [Equation 18] represents the Fourier transform of －１ is the inverse Fourier transform, d is the step representing the height of the upper half-space solution, and S33, [Equation 19] Using the difference between [Equation 20] Correct the new outlier [0000] Obtained, [Equation 22] Here, q represents the step length, which is generally set to 1, and S34, repeat S32 and S33 to obtain the least squares optimization iterative solution formula; [Equation 23] [0000] in the case of, [Equation 25] from [Equation 26] where δ represents a predetermined number close to zero, [0000] satisfies the design accuracy, that is, if u<U, stop the iteration, set the number of iterations u to 30 to 50, and set the following equation: [0000] S35, gravity anomaly values ​​at all points on the plane ΨP [0000] Calculate the equation, and use the solution formula in the frequency domain of the upper half space of the first class boundary value problem of Laplace's equation. [Equation 30] Perform an inverse Fourier transform, [Equation 31] Gravity anomaly at the depth d of the source center in the lower half space [Equation 32] and obtaining:

3. In the step S4, the field source center depth is set as the top surface of the field source, the top surface is divided into grid prisms to construct the field source shape, and density inversion is performed on the top surface of each prism in the frequency domain to obtain the density distribution of the gravitational field source; gravity anomaly [Equation 33] Divide the irregular geological body causing the above into A B equal-sized vertical small prisms of finite extension, where A B = A × B; The coordinates of the center point of the small prism are (x ０ｉ , y ０ｉ ), the residual density is θi, where i = 1, 2, 3, ..., AB, the side lengths are a and b, the prism top surface burial depth is d, and the height is cd. In the frequency domain, each prism is located at a known field source center depth at the top interface (x ｄ , y ｄ Spectrum of gravity anomalies generated at [Equation 34] We obtain the formula [Equation 35] Based on the property of superposition, the forward inversion formula of the gravity anomaly of the field source of known physical parameters, that is, the entire depth layer is known, the field source center depth top interface (x ｄ , y ｄ Spectrum of gravity anomaly G(x, y, 0) generated at [Equation 36] Obtained, [Equation 37] where s and t are spatial frequencies, and x and y are wave numbers, respectively, and G is the gravitational constant. [Number 38] When the gravity anomaly G(x, y, 0) generated by the corresponding top interface of the deep layer is known, d=0; if d≠0, i.e., the buried depth is not zero, calculate the gravity anomaly data of the top interface of the deep layer by solving the least-squares optimization iteration in the lower half space; and [0.39] According to this, the following formula is set: [Equation 40] [Equation 41] and [0.001] From this, we obtain the following equation: [Equation 43] [Equation 44] According to this, the following formula is derived: [Equation 45] [Equation 46] , and perform an inverse Fourier transform on the depth layer top interface at each point (x ｄ , y ｄ , 0) to obtain the residual density values, i.e., to estimate the inversion equations of the relevant physical parameters from the known field source gravity anomalies; [Equation 47] 3. The method of claim 2, further comprising estimating an inversion equation for relevant physical parameters from known field source gravity anomalies and calculating apparent density values ​​corresponding to the gravity anomalies.

4. A high-precision density inversion system for gravity anomaly based on regular stable fast solution, the system being used for implementing the method according to any one of claims 1 to 3, comprising: 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 gravity field and corresponding boundary conditions, i.e., gravity data of ground observation; The second calculation module is used to select a stable Fourier series with regular factors in the lower half space to describe the gravity solution in the lower half space, and make the center of the gravity anomaly correspond to the field source center; The iterative module uses the data containing the false anomaly at the center of the field source depth in the lower half space as the initial boundary value condition of the upper half space, calculates the ground gravity anomaly data based on the boundary value condition through the solution of the upper half space, and performs least squares optimization iterative solution with the actual ground gravity data to obtain the gravity anomaly information of the lower half space that does not contain the Gibbs effect at the center of the field source depth; The inversion module is used to obtain the density distribution at the gravitational field source by dividing the top surface of the field source into grid prisms at the center depth of the field source and performing density inversion at each prism top surface in the frequency domain.

5. In the first calculation module, the step of 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 of the ground observation, includes the following equation: [Number 48] where: [Number 49] is the gravity field value at the measurement point (x, y, w), w is the altitude of the observation point, and n is a discrete series order from 0 to N. [Number 50] is the frequency, μ ｎ 5. The system of claim 4, wherein: is a coefficient, and e is a natural constant.

6. In the second calculation module, the step of selecting a stable Fourier series with regular factors in the lower half space to describe the gravity solution in the lower half space and making the center of gravity anomaly correspond to the field source center includes the following formula: [0.51] where η ｎ is the regularization factor, r is a multiple of the measurement point interval, and increases as the series order increases. [Number 52] is the gravitational field value at the measurement point (r, w), μ ｎ 6. The system of claim 5, wherein: is a coefficient.

Images