Gravity and magnetic field imaging method and system
By integrating the secondary vertical derivative method and Pad's approximate expansion downward extension technology, combined with median filtering, the problems of low-pass filtering and depth-weighting complexity in existing heavy magnetic field imaging technologies are solved, achieving higher imaging accuracy and calculation simplification.
Patent Information
- Application Number
- CN202210474484.8
- Authority / Receiving Office
- CN · China
- Patent Type
- Patents(China)
- Current Assignee / Owner
- Filing Date
- 2022-04-29
- Publication Date
- 2025-05-02
- Estimated Expiration
- 2042-04-29
AI Technical Summary
The existing heavy magnetic field imaging technology has low-pass filtering properties in data processing, suppressing high-frequency information, resulting in insufficient extraction of effective information of modeled data, and complex settings of depth weighting parameters, affecting the accuracy of imaging results.
The vertical derivatives of each order of the heavy magnetic field data are calculated by integrating the secondary vertical derivative method, and the downward extension result of the heavy magnetic field data is calculated by Pad's approximate expansion, and finally median filtering is performed to improve the accuracy of the imaging results.
It improves the accuracy of heavy magnetic field imaging, reduces the dependence on depth weighting, reduces the computational complexity, avoids the tedious operation of early field source distribution prediction, and improves the accuracy of heavy magnetic field exploration data modeling.
Smart Images
Figure CN114779357B_ABST
Abstract
Description
Technical Field
[0001] The present invention relates to the field of earth exploration and information technology, and in particular to a gravity and magnetic field imaging method and system. Background Art
[0002] Gravity and magnetic field imaging technology is the main means to serve the efficient modeling of the gravity and magnetic field inside the earth. The overall requirement for establishing the density and magnetic field model inside the earth is high accuracy, and the requirement for gravity and magnetic field imaging to achieve modeling is also high accuracy. However, there are certain problems with the existing gravity and magnetic field imaging related technologies that have been developed: First, since the current gravity and magnetic field data imaging method is based on upward continuation, this technology has a low-pass filter property that can ensure the stability and accuracy of data processing, but due to this property, it will also suppress the effective information of the high-frequency part of the data, resulting in insufficient extraction of effective information from the data used for modeling and insufficient accuracy; secondly, in order to avoid the results being concentrated in the shallow part, the imaging method uses depth weighting in the processing process, and the depth weighting parameter setting directly affects the accuracy of the imaging results. In theory, to obtain the type of gravity and magnetic field target in advance, it is necessary to know it in advance through other channels, which complicates the imaging processing and contradicts the advantage of gravity and magnetic field imaging technology in modeling.
[0003] Gravity and magnetic field imaging technology is a simple, efficient and stable means of modeling the gravity and magnetic field inside the earth. Its results are comparable to three-dimensional discrete physical property inversion. Gravity and magnetic field imaging technology is based on the upward extension of gravity and magnetic field data. The equivalent gravity and magnetic field distribution inside the earth is obtained through the depth weighting function, and then the corresponding geological body is imaged. The upward extension calculation of gravity and magnetic field data is stable and the results are accurate (Fedi, 2007). Depth weighting can improve the vertical distribution of gravity and magnetic field imaging. Therefore, gravity and magnetic field data imaging technology has been widely developed and applied in the field of geophysical modeling.
[0004] In 2007, Fedi et al. proposed a method to determine the depth of the extreme point (DEXP) based on the obvious power function law between the upward extension field and the upward extension height. This power function law is determined by the scale function r, which can be expressed as
[0005]
[0006] Among them, r represents the gravity and magnetic field, z represents the height of the upward extension, z0 represents the unknown source depth, ζ is the estimated source depth, and N is the structural index (a unitless quantity related to the morphology of the source body, generally 0, 1, 2). The structural index N here requires a specific method to estimate. Here, the morphological information of the source needs to be estimated in advance, which will cause trouble for imaging calculations. In addition, when calculating the unknown source depth z0, it is necessary to continuously change ζ until r(z,ζ) becomes a straight line with a slope of approximately 0, and the structural index N can be estimated. At this time, the DEXP method can be expressed as
[0007]
[0008] Among them, Ω(r,z i ) indicates that at z i The gravity and magnetic field of the DEXP scale at this altitude, T(r,z i ) indicates that at z i The gravity and magnetic field at this height is extended upward. represents the power exponential function of the height of the DEXP scale, is a quantity related to the field source. The results show that Ω(r,z i ) corresponds to the position of the field source, that is, the depth position of the field source is given by the inverse of the height of the extreme point. So we can replace z in formula (2) with i Replace with h i ,
[0009]
[0010] Among them, h i =|z i |.
[0011] It is not difficult to see that the key to imaging the gravity and magnetic field using the DEXP method based on upward continuation is depth weighting. This key quantity requires a lot of preliminary work, which affects the rapid imaging of gravity and magnetic fields. In addition, this method is based on upward extension. Although the general upward extension method is stable and accurate, its inherent low-pass filtering and smoothing effect will lose the detailed information in the gravity and magnetic field data. Summary of the invention
[0012] The technical problem to be solved by the present invention is to provide a gravity-magnetic field imaging method and system, which can focus on improving accuracy and reducing dependence on depth weighting under the premise of ensuring stability.
[0013] In order to solve the above technical problems, the present invention provides a gravity and magnetic field imaging method, the method comprising: using the integrated quadratic vertical derivative method to calculate the vertical derivatives of gravity and magnetic field data {T n}; Substitute the measured first-order vertical derivative into the integrated second-order vertical derivative method, and the integrated second-order vertical derivative method degenerates into a method that only needs to be calculated in the spatial domain, without the need to calculate the gravity and magnetic field integral in the wave number domain; Use the obtained vertical derivatives of each order {T n}Calculate the coefficient equation of the downward extension of the Padé approximate expansion And solve the coefficient {p l}, {q m}; Using the coefficient {p l}, {q m Substitute into the Padé approximation expansion Calculate the downward extension result T(x,y,z0+h) of the gravity and magnetic field data; perform median filtering on the obtained downward extension result, and directly image the filtered result to obtain the equivalent gravity and magnetic field imaging result.
[0014] In some embodiments, the vertical derivatives of each order of gravity and magnetic field data are calculated using the integrated second-order vertical derivative method {T n}, including: measured first-order vertical derivatives; using the second-order and higher-order finite difference method to calculate the horizontal derivatives of the gravity potential UT(x,y,z0+h); using the Laplace equation to obtain the vertical derivatives of the gravity potential UT(x,y,z0+h).
[0015] In some embodiments, the horizontal derivatives of the gravity potential UT(x, y, z0+h) are calculated using a second-order or higher-order finite difference method, including: calculating the horizontal second-order derivative of the gravity potential UT(x, y, z0+h) according to the following formula:
[0016]
[0017] Wherein, Δx, Δy represent the sampling interval in the horizontal direction.
[0018] In some embodiments, the method is exemplified by second order.
[0019] In some embodiments, the Laplace equation is used to obtain the vertical derivatives of the gravity potential UT(x, y, z0+h), including: calculating the vertical second derivative of the gravity potential according to the following formula:
[0020] T(x,y,z0+h)=UT zz (x,y,z0+h)=-UT xx (x,y,z0+h)-UT yy (x,y,z0+h).
[0021] In some embodiments, the vertical derivatives of each order {T n}Calculate the coefficient equation of the downward extension of the Padé approximate expansion And solve the coefficient {p l}, {q m},include:
[0022]
[0023]
[0024] Among them, a n is the coefficient of Taylor expansion, {p l}, {q m} are the coefficients of the Padé approximate expansion.
[0025] In some embodiments, the obtained downward extension result is subjected to median filtering, and the filtered result is directly imaged to obtain an equivalent gravity and magnetic field imaging result, including: defining a K*K window with an odd length, K=2N+1, N is a positive integer; assuming that at a certain position, the equivalent gravity and magnetic field data elements in the window are T(xK,yk,z0+h), ..., T(x,y,z0+h), ..., T(x+K,y+k,z0+h), where T(x,y,z0+h) is the equivalent gravity and magnetic field data element value located at the center of the window; after arranging these K*K equivalent gravity and magnetic field data values in ascending order, the value T middle Replace the value T(x,y,z0+h) at the position (x,y,z0+h), and the result is defined as the output of the median filter.
[0026] In some embodiments, the method uses a window of K=3.
[0027] In addition, the present invention also provides a gravity-magnetic field imaging system, which includes: one or more processors; a storage device for storing one or more programs; when the one or more programs are executed by the one or more processors, the one or more processors implement the gravity-magnetic field imaging method described above.
[0028] After adopting such a design, the present invention has at least the following advantages:
[0029] The present invention belongs to earth exploration and information technology in the field of geological exploration technology and engineering. The present invention can perform equivalent gravity and magnetic field imaging based on downward extension of Padé expansion on data measured by satellite, aviation, ground and ocean. Due to the introduction of measured first-order vertical derivatives, the integrated second-order vertical derivative method degenerates into a method for calculating only in the space domain, avoiding the reduction of convergence space caused by wave number domain calculation; since the Padé expansion downward extension method is stable, accurate, and deep, and the Padé expansion downward extension can directly extend downward through the field source body, the stability and precision of gravity and magnetic field imaging are guaranteed; and since depth weighting is not used in the gravity and magnetic field imaging process, the complexity of calculation is reduced, avoiding cumbersome operations such as early estimation of field source distribution. This is of great significance to the improvement of the accuracy of gravity and magnetic field exploration data modeling. BRIEF DESCRIPTION OF THE DRAWINGS
[0030] The above is only an overview of the technical solution of the present invention. In order to more clearly understand the technical means of the present invention, the present invention is further described in detail below in conjunction with the accompanying drawings and specific implementation methods.
[0031] Figure 1 is a flow chart of the method;
[0032] Figure 2 It is the gravity anomaly caused by four spherical gravity field sources with different underground distribution and density;
[0033] Figure 3 It is a gravity field imaging that extends downwards by 5[1 / 4] order Padre expansion. DETAILED DESCRIPTION
[0034] The preferred embodiments of the present invention are described below in conjunction with the accompanying drawings. It should be understood that the preferred embodiments described herein are only used to illustrate and explain the present invention, and are not used to limit the present invention.
[0035] Although a large number of studies on gravity and magnetic field imaging have been carried out and serve modeling work, the overall requirement for establishing the Earth's internal density and magnetic field model is high accuracy. As the main means of modeling, the existing gravity and magnetic field imaging technology faces obvious difficulties:
[0036] First of all, since the current gravity and magnetic field data imaging technology is based on upward continuation, this technology has a low-pass filtering property that can ensure the stability and accuracy of data processing. However, due to this property, it will also suppress the effective information in the high-frequency part of the data, resulting in insufficient extraction of effective information from the data used for modeling and low resolution of the imaging results.
[0037] Secondly, in order to avoid the results being concentrated in the shallow part, imaging technology adopts depth weighting in the processing process. The depth weighting parameter setting directly affects the accuracy of the imaging results. In theory, to obtain the type of gravity and magnetic field target in advance, it is necessary to know it in advance through other means. This complicates the imaging processing and is not suitable for work areas with little preliminary work.
[0038] Therefore, gravity and magnetic field data imaging technology should focus on improving accuracy and reducing dependence on depth weighting while ensuring stability. In view of the above situation, the present invention proposes a gravity and magnetic field imaging technology based on the downward extension of the spatial domain Padé expansion of the integrated vertical derivative.
[0039] In order to solve the shortcomings of the existing gravity and magnetic field imaging technology and achieve the above-mentioned efficient and precise gravity and magnetic field imaging purpose, the present invention proposes a gravity and magnetic field imaging technology based on the downward extension of the spatial domain Padé approximation of the integrated vertical derivative. Figure 1 , the specific steps of the technical solution are as follows:
[0040] S11, using the integrated secondary vertical derivative method (ISVD method) to calculate the vertical derivatives of gravity and magnetic field data {T n};
[0041] S12, the measured first-order vertical derivative is substituted into the integrated second-order vertical derivative method, and the integrated second-order vertical derivative method is degenerated into a method that only needs to be calculated in the space domain, without the need to calculate the gravity and magnetic field integral in the wave number domain;
[0042] S13, using the vertical derivatives of each order {T n}Calculate the coefficient equation of the downward extension of the Padé approximate expansion And solve for the coefficient p L ,q M ;
[0043] S14, using the coefficient p obtained in step 2 L ,q M Substitute the Padé approximation into Calculate the downward extension result T(x,y,z0+h) of gravity and magnetic field data;
[0044] S15, performing median filtering on the downward extension result obtained in step 3, the parameters of the median filter can be selected as 3*3 or 5*5, and directly imaging the filtered result can obtain an equivalent gravity and magnetic field imaging result.
[0045] The present invention belongs to earth exploration and information technology in the field of geological exploration technology and engineering. The present invention can perform equivalent gravity and magnetic field imaging based on downward extension of Padé expansion on data measured by satellite, aviation, ground and ocean. Due to the introduction of measured first-order vertical derivatives, the integrated second-order vertical derivative method degenerates into a method for calculating only in the space domain, avoiding the reduction of convergence space caused by wave number domain calculation; since the Padé expansion downward extension method is stable, accurate, and deep, and the Padé expansion downward extension can directly extend downward through the field source body, the stability and precision of gravity and magnetic field imaging are guaranteed; and since depth weighting is not used in the gravity and magnetic field imaging process, the complexity of calculation is reduced, avoiding cumbersome operations such as early estimation of field source distribution. This is of great significance to the improvement of the accuracy of gravity and magnetic field exploration data modeling.
[0046] In the above step S11, {T n} is the vertical derivative of gravity and magnetic field at the height z0 (also the height of the observation surface) (not the actual observed data). Its general expression is:
[0047]
[0048] {T n The calculation method of} is the integrated secondary vertical derivative method (ISVD method), which is a mature method. The specific calculation process is as follows:
[0049] First, the gravity field data conversion method in the wave number domain is used to calculate the gravity potential UT(x,y,z0+h) of the gravity field T(x,y,z0+h):
[0050]
[0051] Where F represents Fourier transform, F -1 represents the inverse Fourier transform, k x ,k y represents the wave number domain coordinates corresponding to the spatial domain coordinates x, y, represents the vertical integration operator in the wavenumber domain.
[0052] Then, the finite difference method is used to calculate the second-order derivative of the gravity potential UT(x,y,z0+h) in the horizontal direction:
[0053]
[0054] Wherein, Δx, Δy represent the sampling interval in the horizontal direction.
[0055] Finally, using the Laplace equation, we can obtain the second-order derivative of the gravity potential UT(x,y,z0+h) in the vertical direction, that is, the first-order derivative T1(x,y,z0+h) of the gravity and magnetic field data T(x,y,z0+h) in the vertical direction:
[0056] T(x,y,z0+h)=UT zz (x,y,z0+h)=-UT xx (x,y,z0+h)-UT yy (x,y,z0+h) (7)
[0057] Similarly, other vertical derivatives higher than the first order can also be obtained by using the vertical derivatives obtained in the previous step through equations (5) and (6) to obtain other vertical derivatives two orders higher.
[0058] In the above step S12, the measured first-order vertical derivative is substituted into the integrated second-order vertical derivative method, eliminating the wave number domain calculation of the integrated second-order vertical derivative method, that is, omitting the calculation of formula (4), and only needing to repeatedly calculate formulas (5) and (6) in the space domain to complete the calculation of each order of vertical derivatives. The integrated second-order vertical derivative method degenerates into a method that only needs to be calculated in the space domain, and there is no need to calculate the gravity and magnetic field integral in the wave number domain.
[0059] The coefficients {p L}, {q M}, from (3), we can get the Taylor expansion The coefficients are as follows:
[0060]
[0061] Among them, h is the depth of downward extension, and the Taylor expansion and Padre expansion have the following corresponding relationship
[0062]
[0063] in, {p L}, {q M} are the coefficients of the Padé expansion, so,
[0064]
[0065] Therefore, we can get the following linear equations
[0066]
[0067] After determining the order L and M, the Taylor expansion coefficients {a n Substituting equation (7) into the linear equation system (10), solving the linear equation system (10) can obtain the Padé approximate expansion downward continuation equation of the gravity and magnetic field: The coefficients {p L}, {q M}.
[0068] In the above step S14, {pL}, {q M The obtained T(x,y,z0+h) is the downward extension of the gravity and magnetic field based on the Padre expansion.
[0069] The median filtering in the above step S15 refers to the filtering process of an equivalent gravity and magnetic field data sequence T(x, y, z0+h) obtained after downward extension. First, a K*K window with an odd length is defined, K=2N+1, and N is a positive integer. Suppose at a certain position, the equivalent gravity and magnetic field data elements in the window are T(xK, yK, z0+h), ..., T(x, y, z0+h), ..., T(x+K, y+K, z0+h), where T(x, y, z0+h) is the equivalent gravity and magnetic field data element value located at the center of the window. After arranging these K*K equivalent gravity and magnetic field data values in ascending order, the value T middle Replace the value T(x,y,z0+h) at the position (x,y,z0+h), and the result is defined as the output of the median filter.
[0070] The present invention uses the fifth-order [5=1 / 4] Padé approximation expansion to image the gravity field of the model, so the expression used is:
[0071]
[0072] The calculation of each order vertical derivative in formula (3) is obtained by integrating the second order vertical derivative method (ISVD method).
[0073] Substituting the calculated vertical derivatives of each order into {a5}, and then solving the linear equation determined by equation (10), the coefficients of the 5[1 / 4]th order Padé expansion are obtained as follows:
[0074]
[0075] The gravity anomalies caused by four spherical gravity field sources with different underground distribution and density (such as Figure 2 ) Perform a 5[1 / 4]-order Padre expansion downward extension of gravity field imaging, and obtain the following results Figure 3 .
[0076] See also Figure 2 The white lines are the projections of the four spheres on the observation surface. The density difference between the four spheres and the surrounding rocks is 0.3×10 3 kg / m 2 , 0.8×10 3 kg / m 2 , 1.0×10 3 kg / m 2 and 1.8×10 3 kg / m 2The center burial depths are 20m, 50m, 40m, and 70m respectively; along the x-axis, they are -120m, -40m, 40m to 120m respectively, along the y-axis, they are all 0m, and the z-axis is positive downwards, which means it represents depth.
[0077] See also Figure 3 , the white lines are the projections of the four spheres on the vertical slices, and the vertical slices are along the y-axis; the areas with abnormal displays are the equivalent gravity fields obtained.
[0078] right Figure 2 and Figure 3 Added description and instructions: Figure 2 :Gravity field anomalies are generally caused by the density difference between underground rocks, ore bodies or structures and the surrounding geological bodies. Now suppose there are four spheres underground, and their densities are different. The density of the surrounding rocks outside the spheres is uniform. The density difference between the four spheres and the surrounding rocks is 0.3×10 3 kg / m 2 , 0.8×10 3 kg / m 2 , 1.0×10 3 kg / m 2 and 1.8×10 3 kg / m 2 The central burial depths are 20m, 50m, 40m, and 70m respectively; along the x-axis, they are -120m, -40m, 40m, and 120m respectively, and along the y-axis, they are all 0m. Due to the density difference between these four spheres and the surroundings, anomalies in the gravity field are caused. Figure 2 The following are the results of the gravity field anomaly plane contour line of the 0m depth measurement surface (unit is mGal). Here, the 0m depth measurement surface is assumed to be the ground or general observation surface. The white circles are the positions and geometric projections of four spheres with different buried depths on the 0m depth measurement surface. The x-axis and y-axis can be regarded as the coordinates of the measurement surface in the horizontal direction.
[0079] Figure 3 : Figuratively speaking, it can be considered Figure 3 is Figure 2 The gravity field anomaly is the result of a cut from the straight line y=0 downward to the ground. Gravity field anomalies exist not only on the ground and above, but also underground. The downward extension of the gravity field anomaly described in the present invention can complete the conversion of the underground gravity field anomaly. However, due to certain uncertainties and theoretical difficulties in the converted underground gravity field anomaly, we call the gravity field anomaly below the ground obtained by downward extension the equivalent gravity field anomaly. Figure 3 The position is the vertical slice obtained by cutting downward (along the positive direction of the z-axis) perpendicular to the xy plane along the line y=0 parallel to the x-axis. Figure 3The equivalent gravity field anomaly result (in mGal) obtained by downward extension conversion in this paper on the vertical plane obtained by vertical slicing is shown in Figure 1. The white circles are the positions and geometric shapes of four spheres at different underground depths on the vertical plane obtained by vertical slicing. The depth axis (z axis in the three-dimensional rectangular coordinate system) is the new vertical axis, and the x axis is still the horizontal axis.
[0080] The above description is only a preferred embodiment of the present invention and does not limit the present invention in any form. Those skilled in the art may make some simple modifications, equivalent changes or modifications using the technical contents disclosed above, which all fall within the protection scope of the present invention.
Claims
1. A gravity magnetic field imaging method, characterized in that: include: The vertical derivatives of gravity and magnetic field data are calculated by integrating the second-order vertical derivative method {T n },include: Measured first-order vertical derivatives; Use the second-order and higher-order finite difference method to calculate the horizontal derivatives of the gravity potential UT(x,y,z0+h); Using the Laplace equation, we can obtain the vertical derivatives of the gravity potential UT(x,y,z0+h); Substituting the measured first-order vertical derivative into the integrated second-order vertical derivative method, the integrated second-order vertical derivative method degenerates into a method that only needs to be calculated in the space domain, without the need to calculate the gravity and magnetic field integrals in the wave number domain; Using the obtained vertical derivatives of various orders {T n }Calculate the coefficient equation of the downward extension of the Padé approximate expansion And solve the coefficient {p l }, {q m },include: Among them, a n is the coefficient of Taylor expansion, {p l }, {q m } are the coefficients of Padé's approximate expansion; Using the coefficient {p l }, {q m Substitute into the Padé approximation expansion Calculate the downward extension result T(x,y,z0+h) of gravity and magnetic field data; Perform median filtering on the downward extension results, and directly image the filtered results to obtain equivalent gravity and magnetic field imaging results, including: Define a K*K window with an odd length, K=2N+1, N is a positive integer; Assume that at a certain position, the equivalent gravity and magnetic field data element in the window is T(xK, yk, z0+h), ..., T(x, y, z0+h), ..., T(x+K, y+k, z0+h), where T(x, y, z0+h) is the equivalent gravity and magnetic field data element value located at the center of the window; After arranging these K*K equivalent gravity and magnetic field data values from small to large, the value T middle Replace the value T(x,y,z0+h) at the position (x,y,z0+h), and the result is defined as the output of the median filter.
2. The gravity magnetic field imaging method according to claim 1, characterized in that: The second-order and higher-order finite difference methods are used to calculate the horizontal derivatives of the gravity potential UT(x,y,z0+h), including: The second-order derivative of the gravity potential UT(x,y,z0+h) in the horizontal direction is calculated according to the following formula: Wherein, Δx, Δy represent the sampling interval in the horizontal direction.
3. The gravity magnetic field imaging method according to claim 2, characterized in that: Computations are performed using second-order and higher-order finite difference methods.
4. The gravity magnetic field imaging method according to claim 2, characterized in that: Using the Laplace equation, we can obtain the vertical derivatives of the gravity potential UT(x,y,z0+h), including: The second derivative of the vertical direction of the gravity potential is calculated according to the following formula: T(x,y,z0+h)=UT zz (x,y,z0+h)=-UT xx (x,y,z0+h)-UT yy (x,y,z0+h).
5. The gravity magnetic field imaging method according to claim 1, characterized in that: This method uses a window of K=3.
6. A gravity and magnetic field imaging system, characterized in that: include: one or more processors; A storage device for storing one or more programs; When the one or more programs are executed by the one or more processors, the one or more processors implement the gravity and magnetic field imaging method according to any one of claims 1 to 5.