A method for directly obtaining the deep fine structure of underground space by using gravity

Abstract

The present invention relates to a method for directly obtaining the deep fine structure of underground space by using gravity, comprising: S1, obtaining Bouguer gravity anomaly data; S2, obtaining the regularization downward continuation stability coefficient by using the sliding window numerical deviation statistical accumulation estimation method; S3, obtaining the Fourier cosine series of the Bouguer gravity anomaly; S4, obtaining the downward continuation section by using the adaptive stable regularization downward continuation algorithm; S5, obtaining the high-order vertical differential difference according to the downward continuation section; S6, obtaining the low-frequency difference map by combining low-pass filtering; S7, performing depth correction on the low-frequency difference map to obtain the deep fine structure of underground space. The method for directly obtaining the deep fine structure of underground space provided by the present invention solves the stability problem of the downward continuation technology at any depth and the problem of separating the fine structure of underground space including layered structure, deep structure, and small structures, and realizes the direct acquisition of the fine structure and tectonic information of underground space based on gravity.

Description

Technical Field

[0001] The present invention relates to the technical field of geophysical exploration, and particularly to a method for directly obtaining the deep fine structure of the underground space by using gravity. Background Art

[0002] The geophysical exploration work covering the whole country mainly focuses on gravity and magnetic methods. Taking gravity as an example, after nearly 70 years of accumulation, most parts of China have basically completed the 1:1,000,000 gravity coverage below the snow line. The main ore-forming belts have basically achieved the 1:200,000 regional gravity coverage, and the vast areas in the east and the main ore concentration areas in the northwest have basically achieved the 1:50,000 gravity coverage. The coverage accuracy of airborne and ground magnetic surveys is even higher. Applying these data to major national needs such as energy resource exploration, geological disaster prevention, and national defense construction has always been the key work in the gravity field. However, due to the relatively low coverage of deep exploration and the insufficient ability of gravity to directly obtain geological information of the underground space, these gravity data have not played a major role in the exploration of the underground space structure.

[0003] Applying gravity to directly obtain the fine structure of the underground space is a worldwide difficult and hot issue in the field of gravity exploration. The traditional methods for obtaining the underground space structure information by gravity can be roughly divided into two categories: The first category is the inversion of gravity data, including main methods such as automatic inversion and man-machine interactive inversion. It often requires prior information to construct the initial model of the underground space structure. The automatic inversion can only adapt to the situation with relatively simple structures and has a shallow inversion depth. The man-machine interactive inversion can obtain the density structure characteristics of the relatively complex underground space, but it needs to refer to the interpretation results of electrical methods, seismic methods, etc. as a reference, edit a more detailed initial model, and lacks the ability of direct imaging during the input inversion process. At the same time, whether it is the automatic inversion or the man-machine interactive inversion method, both have a large amount of calculation and are difficult to adapt to the acquisition of the underground space structure in a large range and large area. The second category is the downward continuation method. Due to the poor stability of the traditional downward continuation method, it is difficult to obtain the information of geological bodies at greater depths, and its ability to depict the deep fine structure including small structures and hidden structures is very weak, thus seriously affecting the application of the gravity method in the geological field and its own development.

[0004] Therefore, there is an urgent need for a method that can directly obtain the deep fine structure of the underground space by using gravity. Summary of the Invention

[0005] (1) Technical Problems to be Solved

[0006] In view of the above-mentioned shortcomings and deficiencies of the prior art, the present invention provides a method for directly obtaining the deep fine structure of the underground space by using gravity.

[0007] (2) Technical Solutions

[0008] To achieve the above object, the main technical solutions adopted by the present invention include:

[0009] In a first aspect, an embodiment of the present invention provides a method for directly obtaining the deep fine structure of an underground space by using gravity. The core is based on the gravity data (mainly Bouguer gravity anomaly data) obtained for the corresponding area of the underground space to be measured, and an algorithm for obtaining the underground space structure by downward continuation. In the downward continuation calculation, a sliding window numerical deviation statistical accumulation estimation method is used to obtain an adaptive regularization downward continuation stability coefficient, and further, a high-order vertical differential method is used to separate the fine structure of the downward continuation section. Combining low-pass filtering, the fine structure and tectonic information of the underground space are directly obtained based on gravity.

[0010] The method includes:

[0011] S1. Obtain the Bouguer gravity anomaly data of the corresponding area of the underground space to be measured;

[0012] S2. According to the Bouguer gravity anomaly data, use the sliding window numerical deviation statistical accumulation estimation method to obtain the regularization downward continuation stability coefficient;

[0013] S3. According to the Bouguer gravity anomaly data, obtain the Fourier cosine series of the Bouguer gravity anomaly;

[0014] S4. According to the Bouguer gravity anomaly data, the Bouguer gravity anomaly cosine series, and the regularization downward continuation stability coefficient, use the adaptive stable regularization downward continuation algorithm to control the depth input value z to start from the shallowest downward continuation depth Z0, and increase by a preset depth step Z each time 步长 until z is greater than the maximum downward continuation control depth Z max to obtain the downward continuation profile at each depth, and store all the downward continuation profiles by depth as a data set as the downward continuation section;

[0015] S5. According to the downward continuation section, calculate the K-th order vertical differential difference as the high-order vertical differential difference, where K is an integer and the value range of K is 1 ≤ K ≤ 3;

[0016] S6. According to the high-order vertical differential difference, obtain a low-frequency difference map through combined low-pass filtering;

[0017] S7. Perform depth correction on the low-frequency difference map to obtain the deep fine structure of the underground space;

[0018] Among them, the Bouguer gravity anomaly data includes a Bouguer gravity anomaly data profile Δg with N data points, a length of L, and a point spacing of Δx.

[0019] Optionally, S2 includes:

[0020] The regularization downward continuation stability coefficient is calculated by using the statistical accumulation estimation method of the sliding window numerical deviation, and the formula is as follows:

[0021]

[0022] Where α γ is the regularization downward continuation stability coefficient, N is the number of data points in the Bouguer gravity anomaly data, Δg is the Bouguer anomaly profile in the Bouguer gravity anomaly data, γ is an adjustment parameter and its value range is 0.9 ≤ γ ≤ 1.1, and m is an integer and m = 1, 2, 3, …, 10.

[0023] Optionally, S3 includes: calculating the Fourier cosine series A0, A1, A2, …, A of the Bouguer gravity anomaly according to the following formula n ,

[0024]

[0025] Where n is an integer and the value range of n is n = 1, 2, 3, …, L, N is the number of data points in the Bouguer gravity anomaly data, Δg is the Bouguer anomaly profile in the Bouguer gravity anomaly data, L is the length of the Bouguer gravity anomaly data, and Δx is the point spacing in the Bouguer gravity anomaly data.

[0026] Optionally, the calculation formula of the adaptive stable regularization downward continuation algorithm is

[0027]

[0028]

[0029] Where U(x, z) is the downward continuation profile with spatial coordinates x and z, α γ is the regularization downward continuation stability coefficient, A0, A n is the Fourier cosine series of the Bouguer gravity anomaly, N is the number of data points in the Bouguer gravity anomaly data, Δx is the point spacing in the Bouguer gravity anomaly data, L is the length of the Bouguer gravity anomaly data, and w is an intermediate variable.

[0030] Optionally, S5 includes:

[0031] S51. When K = 1, calculate the first-order vertical differential difference according to the downward continuation section, and the formula is:

[0032]

[0033] Where u(i, j) is the first-order vertical differential difference, is the downward continuation section dh is an intermediate variable;

[0034] S52. When K > 1, Iterate the formula five continuously for K times to obtain the K-th order vertical differential difference as the high-order vertical differential difference.

[0035] Optionally, the S6 includes:

[0036] S61. Perform a moving average filter on the high-order vertical differential difference to obtain first filtered data, where the moving average filter uses the following formula:

[0037]

[0038] Where, is the first filtered data, m1 is the width of the half window of the moving average filter, is the high-order vertical differential difference;

[0039] S62. Perform M times of Gaussian low-pass filtering on the first filtered data to obtain second filtered data as the low-frequency difference map, where M is a positive integer and the value range is 10 ≤ M ≤ 15;

[0040] The Gaussian low-pass filter uses the normalized optimization algorithm of two-dimensional Gaussian filtering, and the formula is as follows:

[0041]

[0042]

[0043] Where, u G is the second filtered data, u G x is the intermediate result of Gaussian low-pass filtering in the x direction, u is the first filtered data, x1, z1 are the filter center positions, x, z are the coordinates of the data points within the filter window, sum is the accumulation function, σ is the adjustment coefficient and the value range is 0 < σ < 1.0, G(x), G(z) are the filter functions of Gaussian low-pass filtering in the x and z directions.

[0044] Optionally, the S7 includes:

[0045] S71. Use the low-frequency difference map to obtain the downhole viewing depth H d ;

[0046] S72. Use the following formula to perform depth correction on the low-frequency difference map to obtain the deep fine structure of the underground space. The calculation formula for the depth correction is:

[0047] H r = βH dFormula Nine;

[0048] where H r is the depth after depth correction, H d is the downward extended apparent depth, β is the conversion coefficient and 1.02 < β < 1.08.

[0049] Optionally, the value of K is K = 3.

[0050] According to the gravity data of the corresponding area of the underground space to be measured obtained, based on the downward continuation technology, the above method of the present invention uses the sliding window numerical deviation statistical accumulation estimation method to obtain the adaptive regularization downward continuation stability coefficient, solving the stability problem of the downward continuation technology at any depth; further using the high-order vertical differential method to separate the fine structure of the downward continuation section, solving the problem of separating the fine structure of the underground space including layered structure, deep structure, and small structures; and using the low-pass filtering combination technology to eliminate the high-frequency interference caused by the high-order vertical differential difference calculation, solving the problem of extracting the complex underground space structure, and thus forming a brand-new large-depth and refined gravity underground space fine imaging method, realizing the direct acquisition of the fine structure and tectonic information of the underground space based on gravity. This method can be applied to fields such as regional gravity underground space fine structure extraction and identification, deep exploration, regional ore-forming prediction, deep and marginal ore prospecting in mining areas, and gravity three-dimensional modeling in large basins, basin-mountain junctions, orogenic belts, etc., and supports the construction of national and regional three-dimensional geophysical-geological models.

[0051] In a second aspect, the present invention provides a computer, characterized in that it includes a memory, a processor, and a computer program stored in the memory and executable on the processor, and when the processor executes the computer program, it implements the method of directly obtaining the deep fine structure of the underground space by gravity as described in any one of the above first aspects.

[0052] In a third aspect, the present invention provides a computer-readable storage medium, characterized in that a computer program is stored on the computer-readable storage medium, and when the computer program is executed by a processor, it implements the method of directly obtaining the deep fine structure of the underground space by gravity as described in any one of the above first aspects.

[0053] (III) Beneficial Effects

[0054] Compared with the prior art, the above method of the present invention, based on the acquired gravity data of the corresponding area of the underground space to be measured, uses the downward continuation technology and the sliding window numerical deviation statistical accumulation estimation method to obtain the adaptive regularization downward continuation stability coefficient, solving the stability problem of the downward continuation technology at any depth; further using the high-order vertical differential method to separate the fine structure of the downward continuation section, solving the problem of separating the fine structure of the underground space including layered structures, deep structures, and small structures; and using the low-pass filtering combination technology to eliminate the high-frequency interference caused by the high-order vertical differential difference calculation, solving the problem of extracting the complex underground space structure, thus forming a brand-new large-depth and refined gravity underground space fine imaging method, realizing the direct acquisition of the fine structure and tectonic information of the underground space based on gravity. This method can be applied to fields such as regional gravity underground space fine structure extraction and identification, deep exploration, regional ore-forming prediction, deep-edge ore prospecting in mining areas, and gravity three-dimensional modeling in large basins, basin-mountain junctions, orogenic belts, etc., and supports the construction of national and regional three-dimensional geophysical-geological models. BRIEF DESCRIPTION OF THE DRAWINGS

[0055] Figure 1 FIG. is a schematic diagram of the main process of the method for directly obtaining the deep fine structure of the underground space using gravity provided by an embodiment of the present invention;

[0056] Figure 2 FIG. is a schematic diagram of the data flow of the method for directly obtaining the deep fine structure of the underground space using gravity provided by an embodiment of the present invention;

[0057] Figure 3 FIG. is the comprehensive test effect diagram of the single-fault model of layered media;

[0058] Figure 4 FIG. is the comprehensive test effect diagram of the combined model of multi-layer and multi-fault media;

[0059] Figure 5 FIG. is the comparison effect diagram of the directly extracted fine structure of the underground space by gravity and the seismic interpretation result of the north-south seismic profile in the Wugong area of the Guanzhong Basin;

[0060] Figure 6 FIG. is the detailed comparison diagram of the north-south seismic profile in the Wugong area of the Guanzhong Basin and the directly extracted fine structure of the underground space by gravity;

[0061] Figure 7 FIG. is the comprehensive comparison diagram of the effect of the gravity-extracted underground space structure and the long-period magnetotelluric inversion interpretation of the same profile in the deep exploration of the Gonghe Basin. DETAILED DESCRIPTION OF THE EMBODIMENTS

[0062] To better understand the above technical solution, the exemplary embodiments of the present invention will be described in more detail below with reference to the accompanying drawings. Although the exemplary embodiments of the present invention are shown in the drawings, it should be understood that the present invention can be implemented in various forms and should not be limited by the embodiments set forth herein. On the contrary, these embodiments are provided to enable a clearer and more thorough understanding of the present invention and to fully convey the scope of the present invention to those skilled in the art.

[0063] To better understand the present invention, some terms used in the present invention are explained below.

[0064] Downward continuation technique (or downward continuation technique of gravity and magnetic fields): It is widely used in the processing and interpretation of gravity and magnetic data, such as downward continuation of aeronautical data potential fields, flattening of curved surfaces, inversion of potential field observation planes, and potential field source seeking. A series of classical algorithms have been developed, such as spatial domain interpolation method, conventional downward continuation FFT method, regularization downward continuation method, generalized inverse method, spline function method, integral iteration method, Milin method, and hybrid iteration method, etc.

[0065] Regularization downward continuation method, regularization continuation coefficient, and regularization downward continuation stability coefficient: The regularization downward continuation method is a method in the downward continuation technique. It has the fastest calculation speed and the strongest adaptability. This method requires an appropriate regularization continuation coefficient. Currently, the commonly used methods for selecting the regularization continuation coefficient include the L-curve method, C-norm method, generalized cross-validation (GCV) method, radial spectrum method, etc. The common feature of these methods is that first, a range of regularization parameters is determined, a common ratio less than or greater than 1 is set, and the regularization parameters are decreased or increased within this range. Then, the downward continuation results of all regularization parameters within this range are calculated, and finally, the optimal parameter is determined based on a certain criterion, such as the inflection point or maximum curvature of the L-curve method, the minimum value of the C-norm curve and GCV curve, the inflection point of the radial spectrum method, etc. Therefore, these methods for selecting regularization parameters require a long calculation time to obtain accurate optimal regularization parameters, and their ability to solve practical problems (stability and broad-spectrum) is insufficient. The regularization downward continuation stability coefficient is a simple, practical, and fast-calculating method for calculating the regularization downward continuation stability coefficient proposed by the present invention starting from the regularization downward continuation stability condition and combining the characteristics of the original data.

[0066] Embodiment 1

[0067] As Figure 1 and Figure 2 shown, this embodiment provides a method for directly obtaining the deep fine structure of the underground space using gravity, including the following steps:

[0068] S1. Obtain the Bouguer gravity anomaly data of the corresponding area of the underground space to be measured;

[0069] Among them, the Bouguer gravity anomaly data includes a Bouguer gravity anomaly data profile Δg with the number of data points being N, the length being L, and the point spacing being Δx.

[0070] S2. According to the Bouguer gravity anomaly data, a regularization downward continuation stability coefficient is obtained by using a sliding window numerical deviation statistical accumulation estimation method;

[0071] Specifically, the regularization downward continuation stability coefficient is calculated by using a sliding window numerical deviation statistical accumulation estimation method, and the formula is as follows:

[0072]

[0073] where α γ is the regularization downward continuation stability coefficient, N is the number of data points in the Bouguer gravity anomaly data, Δg j and Δg n (obtained from Δg) is the Bouguer anomaly profile in the Bouguer gravity anomaly data, γ is an adjustment parameter and its value range is 0.9 ≤ γ ≤ 1.1, m is the half-window width value of the sliding window and m = 1, 2, 3,..., 10, and n and j are intermediate variables for accumulation control.

[0074] It should be noted that the above method for obtaining the regularization downward continuation stability coefficient of the present invention starts from the regularization downward continuation stability condition and combines the characteristics of the original data. Its theoretical basis is that the value condition and range equation of the regularization coefficient satisfy the Schwarz inequality:

[0075]

[0076] Among them, δ is the error between the result of upward continuation after downward continuation at depth z and the original value, C0 is a constant, L is the profile length, z is the depth of downward continuation, is the data deviation degree, k is the order, and n is the number of points;

[0077] Taking k = 1 and solving formula 1-1 simultaneously, we get:

[0078]

[0079] Since always holds, so always holds, and always satisfies formula 1-2, so we can take:

[0080]

[0081] In the present invention, in order to better standardize the relationship between the regularization coefficient and the original data and at the same time meet the requirements of different actual data processing, formula 1-3 is improved, and the improved regularization coefficient α γThe calculation formula is as follows:

[0082]

[0083] Among them, for the Bouguer anomaly profile Δg that is distributed along the x - direction, with the number of data points being N, the length being L, and the point - spacing being Δx, a numerical deviation - square cumulative - sum calculation is carried out. The half - width of the moving - average window is m, and the calculation formula is as follows:

[0084]

[0085] According to Formula 1 - 4 and Formula 1 - 5, the first formula of the regularization downward - continuation stability coefficient can be obtained.

[0086] It should be noted that in specific implementation, by adjusting the size of the moving - average window (half - width is m), the downward - continuation problem at large depths can be solved. The larger the moving - average window, the greater the stable downward - continuation depth, but the ability to distinguish fine structures in the later stage is also weakened; the smaller the moving - average window, the smaller the stable downward - continuation depth, and the stronger the ability to distinguish fine structures in the later stage.

[0087] S3. Obtain the Fourier cosine series of the Bouguer gravity anomaly according to the Bouguer gravity anomaly data;

[0088] Specifically, the Fourier cosine series \(A_0\), \(A_1\), \(A_2\), …, \(A_n\) of the Bouguer gravity anomaly can be calculated according to the following formula n ,

[0089]

[0090] Among them, n is an integer and the value range of n is \(n = 1,2,3,\cdots,N\), N is the number of data points in the Bouguer gravity anomaly data, \(\Delta g_0\), \(\Delta g_1\), \(\Delta g_2\) (obtained from \(\Delta g\)) are the Bouguer anomaly profiles in the Bouguer gravity anomaly data, L is the length of the Bouguer gravity anomaly data, \(\Delta x\) is the point - spacing in the Bouguer gravity anomaly data, and i is an intermediate variable. L , \(\Delta g_1\) i (obtained from \(\Delta g\)) are the Bouguer anomaly profiles in the Bouguer gravity anomaly data, L is the length of the Bouguer gravity anomaly data, \(\Delta x\) is the point - spacing in the Bouguer gravity anomaly data, and i is an intermediate variable.

[0091] S4. According to the Bouguer gravity anomaly data, the Bouguer gravity anomaly cosine series, and the regularization downward - continuation stability coefficient, use the adaptive - stability regularization downward - continuation algorithm to control the depth input value z starting from the shallowest downward - continuation depth \(Z_0\), and increasing by a preset depth step \(Z\) 步长 each time until z is greater than the maximum downward - continuation control depth \(Z\) max , obtain the downward - continuation profile at each depth, and store all the downward - continuation profiles by depth as a data set as the downward - continuation section;

[0092] The calculation formula of the adaptive stable regularization downward continuation algorithm is as follows:

[0093]

[0094]

[0095] where U(x, z) is the downward continuation profile with spatial coordinates x and z. is the regularization downward continuation operator, and α γ is the regularization downward continuation stability coefficient, A0, A n is the Fourier cosine series of the Bouguer gravity anomaly. n is an integer and the value range of n is n = 1, 2, 3, …, N, where N is the number of data points in the Bouguer gravity anomaly data, Δx is the point spacing in the Bouguer gravity anomaly data, L is the length in the Bouguer gravity anomaly data, w is an intermediate variable, and e is the natural constant.

[0096] It should be noted that formula three of the above adaptive stable regularization downward continuation algorithm is constructed based on the downward continuation formula of even extension in the frequency domain. The downward continuation calculation formula of even extension in the frequency domain is:

[0097]

[0098] where U(x, z) is the downward continuation profile with spatial coordinates x and z, φ(w, z) is the regularization downward continuation operator, A0, A n is the Fourier cosine series of the Bouguer gravity anomaly. n is an integer and the value range of n is n = 1, 2, 3, …, N, where N is the number of data points in the Bouguer gravity anomaly data, Δx is the point spacing in the Bouguer gravity anomaly data, and w is an intermediate variable that satisfies formula four.

[0099] The calculation formula of the regularization downward continuation operator φ(w, z) is:

[0100]

[0101] where α γ is the regularization downward continuation stability coefficient, w is an intermediate variable that satisfies formula four, e is the natural constant, and z is the downward continuation depth (i.e., the spatial coordinate in the z direction).

[0102] Substituting formula 3-2 into formula 3-1 and combining with the regularization downward continuation stability coefficient obtained according to formula one, formula three of the adaptive stable regularization downward continuation algorithm can be obtained.

[0103] As Figure 1 and Figure 2 shown, the specific steps to obtain the downward continuation profile dataset are as follows:

[0104] Set the shallowest depth of the ground or downward extension as z0, and the maximum controlled depth of downward extension as Z max , and preset the depth step as Z 步长 (or expressed as dh),

[0105] A11. Set the initial depth z = z0;

[0106] A12. Obtain the downward continuation profile U(x, z) of the depth z through Formula 3 and store it in the data set where x and z are the spatial coordinates of the profile;

[0107] A13. Control the depth z = z + dh, and repeat the operation of A12 and A13 until z > Z max ;

[0108] The data set of the downward continuation profile obtained through the above steps A11 - A13 is the downward continuation section that fully covers the range of the downward continuation depth (z0, z max ).

[0109] S5. Calculate the Kth-order vertical differential difference as the high-order vertical differential difference according to the downward continuation section, where K is an integer and the value range of K is 1 ≤ K ≤ 3;

[0110] Specifically, in the implementation process, the vertical differential difference is obtained according to the downward continuation section formed in step S4, that is, for Since the higher the order, the higher the fineness of the extracted structure, but the smaller the difference value, and the smaller the numerical difference between each other. At the same time, the higher the order, the greater the high-order dispersion and high-frequency interference signals. Generally speaking, in practice, the differential order can be selected as 1 - 3 orders according to needs.

[0111] To better illustrate step S5, the following will illustrate step S5 through sub-steps S51 - S52:

[0112] S51. When K = 1, calculate the first-order vertical differential difference according to the downward continuation section, and the formula is:

[0113]

[0114] where u(i, j) is the first-order vertical differential difference (i is the row number, j is the column number), (obtained from , i, i + 1 are the row numbers, j is the column number) is the downward continuation section (can also be abbreviated as ), dh, i, j are intermediate variables;

[0115] S52. When K > 1, The formula five can be iterated continuously for K times to obtain the K-th order vertical differential difference as the high-order vertical differential difference.

[0116] Generally, the high-order vertical differential difference can be expressed as u diff (i, j) or u diff .

[0117] Preferably, for example, in this embodiment, the order is selected such that K = 3. Therefore, in this embodiment, the obtained high-order vertical component difference is the difference of the third-order vertical component.

[0118] S6. According to the high-order vertical differential difference, a low-frequency difference map is obtained through combined low-pass filtering;

[0119] It should be explained that the high-order vertical differential difference obtained in step S5 contains a large amount of dispersion and high-frequency interference signals attached during the high-order differential calculation process, and various structural features therein are still masked. It is necessary to perform combined low-pass filtering to remove the interference in order to better obtain the spatial distribution characteristics and geological representations of the high-order vertical differential difference.

[0120] It should be noted that in this embodiment, the underground space structure information is implemented according to the adaptive regularization downward continuation algorithm. Therefore, the obtained underground space and the subsequent high-order vertical differential difference extraction structure information also have the characteristics of downward continuation, that is, mainly vertical strip signals with large intensity, while the layered and deep signals have weak intensity and are often suppressed and masked by the former. At the same time, there are certain dispersion and high-order interference waves. This requires protecting the horizontal signals as much as possible, suppressing the high-frequency interference signals, and not over-suppressing the vertical signals during the filtering process. According to this data characteristic, the present invention uses moving average filtering and Gaussian low-pass filtering for combined low-pass filtering. Among them, the moving average filtering is only implemented in the x direction, mainly to enhance the continuity of the horizontal signals, play a role in strengthening the protection of the layered formation signals, and at the same time eliminate local horizontal interference; the Gaussian filtering is to eliminate dispersion and high-frequency interference signals in both the x and z directions and highlight the mid-low frequency signals.

[0121] To better illustrate step S6, S6 will be decomposed into sub-steps S61 - S62 for description below.

[0122] S61. Perform moving average filtering on the high-order vertical differential difference to obtain first filtered data, where the moving average filtering uses the following formula:

[0123]

[0124] where, is the first filtered data (i is the row number, j is the column number), m1 is the width of the half window of the moving average filtering, (Obtained from u diff (where k is the number of rows and j is the number of columns) is the high-order vertical differential difference, and i, j, and k are intermediate variables;

[0125] S62. Perform M times of Gaussian low-pass filtering on the first filtered data to obtain second filtered data as the low-frequency difference map, where M is a positive integer and the value range is 10 ≤ M ≤ 15;

[0126] The Gaussian low-pass filtering adopts a normalized optimization algorithm of two-dimensional Gaussian filtering. Its working process is to first perform Gaussian low-pass filtering on all data in the x direction, and then based on this, perform Gaussian low-pass filtering in the z direction. Its advantage is that compared with the traditional two-dimensional Gaussian low-pass filtering, this method is faster, and has stronger stability and adaptability.

[0127] In this embodiment, the formula of the Gaussian low-pass filtering is as follows:

[0128]

[0129]

[0130] where u G is the second filtered data, u G x is the intermediate result of Gaussian low-pass filtering in the x direction, u (i.e., in Formula 6 ) is the first filtered data, x1, z1 are the filtering center positions, x, z are the coordinates of the data points within the filtering window, sum is the accumulation function, σ is the adjustment coefficient and the value range is 0 < σ < 1.0, G(x), G(z) are the filter functions of Gaussian low-pass filtering in the x and z directions, and e is the natural constant.

[0131] It can be understood that in this embodiment, according to the low-pass filtering of the above steps, first perform moving average filtering with 1 filtering time, and then perform Gaussian low-pass filtering with 10 - 15 filtering times, which can achieve a good effect of removing high-frequency interference signals.

[0132] Preferably, as Figure 2 shown, in actual implementation, according to the signal interference elimination situation after low-pass filtering of the above steps, if the interference has not been completely removed, the steps S61 - S62 can be further used for iterative filtering (usually can continue to iterate 2 - 4 times) to achieve a better filtering effect.

[0133] In specific implementation, for example, the Gaussian filtering window width can be selected from 5 groups of filtering forms such as 3 points, 5 points, 7 points, 9 points, and 11 points, and the moving average filtering can be selected from 3 groups of filtering forms such as 5 points, 7 points, and 9 points.

[0134] S7. Perform depth correction on the low-frequency differential map to obtain the fine structure of the deep underground space;

[0135] S71. Use the low-frequency differential map to obtain the downward extended apparent depth H d ;

[0136] S72. Perform depth correction on the low-frequency differential map using the following formula to obtain the fine structure of the deep underground space. The calculation formula for the depth correction is:

[0137] H r = βH d Formula Nine;

[0138] where H r is the depth after depth correction, H d is the downward extended apparent depth, and β is the conversion coefficient with 1.02 < β < 1.08.

[0139] It should be noted that the value range of the conversion coefficient β, 1.02 < β < 1.08, roughly reflects that the downward continuation technology adopted in this embodiment can reach a depth between 0.93 - 0.98 for approximating the true field source, which is higher than the existing level of 0.9 and is mainly related to the specific geological conditions of the gravity data acquisition area.

[0140] After the above correction process, a vertical differential difference map based on downward continuation can be obtained. Through theoretical model testing and comparison, and repeated comparison with the results of geophysical methods such as seismic and electromagnetic methods, it can be determined that this difference map has a good one-to-one correspondence with the actual geological boundary of the underground space, and it is a new type of gravity imaging method with fast calculation speed and high reliability.

[0141] The above method of the present invention, based on the gravity data of the corresponding area of the underground space to be measured obtained, adopts the sliding window numerical deviation statistical accumulation estimation method to obtain the adaptive regularization downward continuation stability coefficient based on the downward continuation technology, so as to solve the stability problem of the downward continuation technology at any depth; further use the high-order vertical differential method to separate the fine structure of the downward continuation section and solve the problem of separating the fine structure of the underground space including layered structure, deep structure, and small structures; and adopt the low-pass filtering combination technology to eliminate the high-frequency interference brought by the high-order vertical differential difference calculation and solve the problem of extracting the complex underground space structure, thereby forming a brand-new large-depth and refined gravity underground space fine imaging method, realizing the direct acquisition of the fine structure and tectonic information of the underground space based on gravity. This method can be applied to fields such as large basins, basin-mountain junctions, orogenic belts, etc. for regional gravity underground space fine structure extraction and identification, deep exploration, regional ore-forming prediction, deep and marginal ore prospecting in mining areas, and gravity three-dimensional modeling, and supports the construction of national and regional three-dimensional geophysical-geological models.

[0142] Example 2

[0143] According to the method for directly obtaining the deep fine structure of underground space using gravity in Example 1 of the present invention, comprehensive tests on layered single-fault medium models, comprehensive tests on multi-layer multi-fault medium combined models, extraction of the deep structure of the gravity profile in the Wugong area of the Guanzhong Basin and fine structure comparison between it and the seismic interpretation results of the same profile, exploration of the deep gravity structure of the Gonghe Basin and structure comparison with the long-period magnetotelluric sounding profile of the same profile, etc. have been carried out, and good results have been achieved.

[0144] The comprehensive test effect of the layered medium single-fault model is as Figure 3 shown. The stratification boundary of the model and the exact position of the spatial fault ( Figure 3 the second part in Figure 3 ), by applying forward modeling of the model to obtain the Bouguer gravity anomaly ( Figure 3 the first part in

[0145] ), it can be seen that the method of the present invention can perform relatively accurate inversion.

[0145] At the same time, in order to better test the adaptability of the present invention to more complex models, the comprehensive test of the multi-layer multi-fault medium model ( Figure 4 the second part in Figure 4 ), the test effect is as Figure 4 shown. First, forward simulation is carried out to obtain the Bouguer anomaly curve ( Figure 4 the first part in

[0146] ), after directly extracting the underground space structure, through systematic comparison, the extraction results can effectively extract all the fine information such as the spatial position of the interface of the layered medium, the occurrence and spatial position of the fault, etc., indicating the effectiveness of the method and its ability to adapt to complex media. Generally speaking, through model testing, this method can effectively adapt to layered media, and at the same time can more precisely extract the real structure of underground space, thus breaking through the key technology of automatically extracting the structure of layered media by gravity, and providing support for regional gravity underground space fine imaging.

[0146] As Figure 5 shows the application effect of the north-south profile in the Wugong area of the Guanzhong Basin (obtained by using the method of Example 1), and as a comparison, as Figure 6 shown, the seismic exploration interpretation results of the same profile. It can be seen that there is a good corresponding relationship between the two in the main horizons, the spatial position and dip of the fault structure, etc. At the same time, due to the influence of the Qinling uplift in the south, there is a bedding collapse on the side of the Guanzhong Basin in the later stage, thus forming obvious bedding structure characteristics, and the north side has obvious thrust nappe structure characteristics. The above two characteristics are better than the seismic effect. By comparison, it can be known that the results obtained by the method of the present invention are comparable to those of seismic, indicating that the method has a high degree of fineness in extracting the underground space structure and can be applied to high-precision gravity imaging scenarios.

[0147] As Figure 7 shown,Figure 7 It shows the results of gravity underground space extraction from the deep exploration profile of the Gonghe Basin. The structural characteristics of the gravity underground space extracted by gravity are highly similar to the tectonic morphology inverted by long-period magnetotellurics, indicating the effectiveness of the gravity method. The results show that the shallow part of the basin presents low-value gravity anomalies, which coincide with the basin range and are consistent with the basin range reflected by the characteristics of low resistivity anomalies. There is a phenomenon of "uplift" of high-density anomalies below 25 km in the basin towards the shallow part, and the uplift amplitude exceeds 30 km. The electrical characteristics also show that the deep high resistivity uplifts upwards to the area above the Moho surface, presenting an overall "M"-type anomaly characteristic with low values on both the north and south sides and high values in the middle. It is speculated that there is an upwelling of mantle materials or a solidified upwelling mantle in the deep part of the basin. The application effect proves that the method of the present invention can apply gravity to the field of deep exploration, and further solve the problem of extracting the fine density structure above the Moho surface, greatly expanding the application scope of gravity in the field of deep exploration.

[0148] Embodiment III

[0149] This embodiment provides a computer, which is characterized by including a memory, a processor, and a computer program stored in the memory and executable on the processor. When the processor executes the computer program, it implements the method for directly obtaining the deep fine structure of the underground space as described in any one of the above Embodiment I;

[0150] In addition, this embodiment also provides a computer-readable storage medium, which is characterized in that a computer program is stored on the computer-readable storage medium, and when the computer program is executed by the processor, it implements the method for directly obtaining the deep fine structure of the underground space as described in any one of the above Embodiment I.

[0151] Embodiment IV

[0152] This embodiment provides an apparatus for obtaining the deep fine structure of the underground space, including:

[0153] A gravity information acquisition unit for acquiring Bouguer gravity anomaly data of the corresponding area of the underground space to be measured;

[0154] A stability coefficient processing unit for obtaining a regularization downward continuation stability coefficient by using the sliding window numerical deviation statistical accumulation estimation method according to the Bouguer gravity anomaly data;

[0155] A Fourier cosine series processing unit for obtaining the Fourier cosine series of the Bouguer gravity anomaly according to the Bouguer gravity anomaly data;

[0156] A cross-section processing unit, configured to control the depth input value z to start from the shallowest downward continuation depth Z0 and increase by a preset depth step Z each time according to the Bouguer gravity anomaly data, the Bouguer gravity anomaly cosine series, and the regularization downward continuation stability coefficient, by using an adaptive stability regularization downward continuation algorithm, until z is greater than the maximum downward continuation control depth Z, obtain the downward continuation profiles at each depth, and store all the downward continuation profiles as a data set according to the depth as the downward continuation cross-section. 步长 , until z is greater than the maximum downward continuation control depth Z max , obtain the downward continuation profiles at each depth, and store all the downward continuation profiles as a data set according to the depth as the downward continuation cross-section;

[0157] A high-order vertical differential processing unit, configured to calculate the Kth-order vertical differential difference according to the downward continuation cross-section as the high-order vertical differential difference, where K is an integer and the value range of K is 1 ≤ K ≤ 3;

[0158] A low-pass filtering unit, configured to obtain a low-frequency difference map by combining low-pass filtering according to the high-order vertical differential difference;

[0159] A correction unit, configured to perform depth correction on the low-frequency difference map to obtain the fine structure of the deep underground space.

[0160] It should be noted that in the description of this specification, the descriptions of terms such as "one embodiment", "some embodiments", "embodiment", "example", "specific example", or "some examples" mean that the specific features, structures, materials, or characteristics described in connection with the embodiment or example are included in at least one embodiment or example of the present invention. In this specification, the schematic representations of the above terms do not necessarily refer to the same embodiment or example. Moreover, the specific features, structures, materials, or characteristics described can be combined in a suitable manner in any one or more embodiments or examples. In addition, without contradiction, those skilled in the art can combine and combine the different embodiments or examples described in this specification and the features of different embodiments or examples.

[0161] Although the preferred embodiments of the present invention have been described, those skilled in the art can make additional changes and modifications after learning the basic creative concepts. Obviously, those skilled in the art can make various modifications and variations to the present invention without departing from the spirit and scope of the present invention.

Claims

1. A method for directly obtaining the deep fine structure of underground space by using gravity, characterized in that The method includes: S1. Obtain Bouguer gravity anomaly data of the corresponding area of the underground space to be measured; S2. According to the Bouguer gravity anomaly data, use the sliding window numerical deviation statistical accumulation estimation method to obtain the regularization downward continuation stability coefficient; S3. According to the Bouguer gravity anomaly data, obtain the Fourier cosine series of the Bouguer gravity anomaly; S4. According to the Bouguer gravity anomaly data, the cosine series of the Bouguer gravity anomaly, and the regularization downward continuation stability coefficient, using the adaptive stable regularization downward continuation algorithm, control the depth input value z to start from the shallowest downward continuation depth Z0, and increase it by a preset depth step Z each time 步长 until z is greater than the maximum downward continuation control depth Z max , obtain the downward continuation profile at each depth, and store all the downward continuation profiles by depth as a data set as the downward continuation section; S5. According to the downward continuation section, calculate the K-th order vertical differential difference as the high-order vertical differential difference, where K is an integer and the value range of K is 1≤K≤3; S6. According to the high-order vertical differential difference, obtain a low-frequency difference map through combined low-pass filtering; S7. Perform depth correction on the low-frequency difference map to obtain the deep fine structure of the underground space; Among them, the Bouguer gravity anomaly data includes: a Bouguer gravity anomaly data profile Δg with the number of data points N, length L, and point spacing Δx; The S2 includes: Calculate the regularization downward continuation stability coefficient using the sliding window numerical deviation statistical accumulation estimation method. The formula is as follows: where α γ is the regularization downward continuation stability coefficient, N is the number of data points in the Bouguer gravity anomaly data, Δg is the Bouguer anomaly profile in the Bouguer gravity anomaly data, γ is an adjustment parameter and its value range is 0.9 ≤ γ ≤ 1.1, m is an integer and m = 1, 2, 3, …, 10.

2. The method according to claim 1, wherein The S3 includes: calculating the Fourier cosine series A0, A1, A2, …, A of the Bouguer gravity anomaly according to the following formula n , Where n is an integer and the value range of n is n = 1, 2, 3,..., N, N is the number of data points in the Bouguer gravity anomaly data, Δg is the Bouguer anomaly profile in the Bouguer gravity anomaly data, L is the length of the Bouguer gravity anomaly data, and Δx is the point spacing in the Bouguer gravity anomaly data.

3. The method according to claim 2, wherein The calculation formula of the adaptive stability regularization downward continuation algorithm is Among them, U(x,z) is the downward continuation profile with spatial coordinates x and z, and α γ is the regularization downward continuation stability coefficient, A0, A n is the Fourier cosine series of the Bouguer gravity anomaly, N is the number of data points in the Bouguer gravity anomaly data, Δx is the point spacing in the Bouguer gravity anomaly data, L is the length in the Bouguer gravity anomaly data, and w is an intermediate variable.

4. The method according to claim 1, wherein The S5 includes: S51. When K = 1, according to the downward continuation section, calculate the first-order vertical differential difference. The formula is: where u(i, j) is the first-order vertical differential difference quantity, for the downward continuation section dh is an intermediate variable; S52. When K > 1, Iterate formula five continuously for K times to obtain the K-th order vertical differential difference as the high-order vertical differential difference.

5. The method according to claim 4, wherein The S6 includes: S61. Perform moving average filtering on the high-order vertical differential difference to obtain the first filtered data. The moving average filtering uses the following formula: Among them, is the first filtered data, and m1 is the width of the half window of the moving average filter, is the high-order vertical differential difference; S62. Perform M times of Gaussian low-pass filtering on the first filtered data to obtain the second filtered data as the low-frequency difference map, where M is a positive integer and the value range is 10≤M≤15; The Gaussian low-pass filtering uses the normalized optimization algorithm of two-dimensional Gaussian filtering. The formula is as follows: where, u G is the second filtered data, u G x is the intermediate result of Gaussian low-pass filtering in the x direction, u is the first filtered data, x1, z1 are the filtering center positions, x, z are the coordinates of the data points within the filtering window, sum is the accumulation function, σ is the adjustment coefficient and its value range is 0 < σ < 1.0, G(x), G(z) are the filter functions of Gaussian low-pass filtering in the x and z directions.

6. The method according to claim 1, wherein The S7 includes: S71. Obtain the downward-looking apparent depth H by using the low-frequency differential map d ; S72. Use the following formula to perform depth correction on the low-frequency difference map to obtain the deep fine structure of the underground space. The calculation formula of the depth correction is: H r = βH d Formula IX; where H r is the depth after depth correction, H d is the downward-looking depth, β is the conversion coefficient and 1.02 < β < 1.

08.

7. The method according to claim 1, characterized in that, The value of K is K = 3.

8. A computer, characterized in that, It includes a memory, a processor, and a computer program stored in the memory and executable on the processor. When the processor executes the computer program, it implements the method for directly obtaining the deep fine structure of the underground space using gravity as described in any one of claims 1 to 7.

9. A computer-readable storage medium, characterized in that, A computer program is stored on the computer-readable storage medium. When the computer program is executed by the processor, it implements the method for directly obtaining the deep fine structure of the underground space using gravity as described in any one of claims 1 to 7.

Images