Gravity anomaly upward continuation calculation method based on central grid data block separation method

Abstract

The invention discloses a gravity anomaly upward continuation calculation method based on a central grid data block separation method, belongs to the technical field of geophysical exploration data processing, and solves the problems that when an existing method is used for processing an ultra-low altitude upward continuation problem, singularity occurs when an integral kernel function coincides with a calculation point and a data point; the method comprises the following steps: stripping a central data grid where a calculation point is located from an integral domain based on a block separation method, and carrying out practical modification processing on a global integral model based on the block separation method to obtain a modified model without integral singularity; according to the method, a practical modification scheme for converting a gravity anomaly upward continuation global integral model into a local numerical integral model is provided, numerical calculation verification is carried out on a kernel function singularity solution and the practical modification scheme of a calculation model by utilizing an ultrahigh-order bit model EGM2008, and model internal coincidence precision superior to 1mGal can be obtained; therefore, the high-precision earth external gravitational field assignment requirement is met.

Description

Technical Field

[0001] The present invention belongs to the technical field of geophysical exploration data processing, and in particular relates to a gravity anomaly upward extension calculation method based on a central grid data block separation method. Background Art

[0002] In the field of geophysical exploration and data processing, the upward extension of gravity anomalies is one of the core technical means to solve problems such as surface and shallow geological structure inversion and mineral resource exploration. Its core goal is to extend the gravity anomaly data collected on the sphere to the external space of the sphere through the potential field theoretical model, providing key data support for subsequent gravity field separation, terrain correction and deep structural analysis. At present, the classical solution based on the Poisson integral formula has long been widely used in the numerical calculation of the upward extension of gravity anomalies because it strictly satisfies the spherical solution characteristics of the first boundary value problem of the potential field (Dirichlet problem) in theory and has the advantages of computational stability and efficiency.

[0003] According to the theory of gravity potential field, the global integral model for calculating the external gravity anomaly of the sphere by extending the gravity anomaly distributed on the sphere upward can be derived from the integral formula of the spherical solution of the first boundary value problem of the potential field, namely the Dirichlet problem:

[0004]

[0005] Where Δg p is the calculation point outside the sphere with a height of h The gravity anomaly at are the geocentric radius, geocentric latitude and longitude of the spherical external space coordinates respectively; R is the average radius of the earth ellipsoid; r = R + h, Δg q Integral flow point for the sphere The known gravity anomaly at the location; σ represents the unit sphere; dσ is the integral area element of the unit sphere; K(r,ψ) is the integral kernel function; l represents the calculation point To points flow point The spatial distance between pq is the spherical angular distance between the calculation point and the flow point; P n (cosψ) is the nth-order Legendre function.

[0006] The Poisson integral formula has strict physical meaning and mathematical rigor in theory. However, when dealing with the problem of upward extension at ultra-low altitude, the integral kernel function in the formula will have serious singularity problems. Specifically, it can be seen from equations (1)-(3) that when r→R and ψ≠0, K(r,ψ)→0, and the corresponding integral term is zero; and when the calculation point approaches the data point, that is, when r→R and ψ→0, the denominator term l→0 will appear. At this time, the integral kernel function K(r,ψ) becomes singular, resulting in uncontrollable uncertainty problems in the numerical calculation of the integral formula (1), which manifests as violent oscillation of the calculation results, sharp error amplification, and even algorithm crash. Therefore, when implementing the numerical calculation of the integral formula (1), special treatment must be made for the singularity problem of the integral kernel function.

[0007] It can be seen from this that in the existing technical solutions, although the Poisson integral formula can achieve high-precision numerical calculations at conventional altitudes, its applicability to ultra-low-altitude scenarios (such as airborne gravity measurements, satellite gravity gradient inversion, etc.) is severely limited. Although some studies have attempted to alleviate the influence of singularities by truncating high-order terms, introducing regularization operators, or adjusting grid partitioning strategies, none of them have fundamentally solved the essential singularity problem of the integral kernel function when r→R. Specifically, when the calculation points and data points coincide with each other due to discretization errors or actual overlap in numerical calculations, the traditional method cannot avoid the zero value problem of the denominator term, resulting in a loss of numerical stability, which in turn affects the reliability of the upward extension results. Therefore, in response to the singularity problem of the integral kernel function in ultra-low-altitude scenarios in the upward extension of gravity anomalies, we proposed a gravity anomaly upward extension calculation method based on the central grid data block separation method. Summary of the Invention

[0008] The purpose of the present invention is to address the shortcomings of the existing technology and provide a gravity anomaly upward extension calculation method based on the central grid data block separation method, which solves the problem that the existing method will cause numerical instability and affect the calculation accuracy when dealing with ultra-low altitude upward extension problems due to the singularity of the integral kernel function when the calculation point coincides with the data point.

[0009] The present invention is implemented as follows: a gravity anomaly upward continuation calculation method based on a central grid data block separation method, the gravity anomaly upward continuation calculation method based on a central grid data block separation method comprises:

[0010] S10, based on the block separation method, the central data grid where the calculation point is located is separated from the integral domain, and the integral kernel function in the small block is processed using a plane approximate analytical method, and a global integral model is outputted based on the central grid data block separation method to remove the singularity of the integral kernel function;

[0011] S20, loading the constructed global integration model based on the block separation method, dividing the global integration domain into a near zone and a far zone, performing practical modification processing on the global integration model based on the block separation method, and obtaining a modified model with removed integration singularities.

[0012] The method for separating the central data grid where the calculation point is located from the integration domain based on the block separation method includes:

[0013] S101, performing a planar approximation process on the integral kernel function within the central grid, and using polar coordinate expansion to approximate the integral kernel function;

[0014] S102, deriving an integral expression for the approximate integral kernel function and outputting a local analytical integral;

[0015] S103 approximates that the integral kernel function is rewritten as the sum of the separated main integral term and the contribution value of the central block. In view of the drastic change in gravity anomaly of the central data grid, an additional correction term is introduced to output a global integral model based on the central grid data block separation method to remove the singularity of the integral kernel function.

[0016] When the integral kernel function is approximated as a plane within the central grid:

[0017] The integral kernel function is approximated in the center grid to satisfy r=R+h;R 2 dσ≈sdsdαAt this time, the defined integral kernel function is approximated by polar coordinate expansion:

[0018]

[0019] When deriving the integral expression of the approximate integral kernel function:

[0020] Assuming that the radius of the central data grid is s0, the contribution of the gravity anomaly of the spherical central data grid to the calculated value of the external gravity anomaly at height h is expressed as:

[0021]

[0022] Within the central data grid, Δg q As a constant, Δg q =Δg Rp , Δg Rp To calculate the known spherical gravity anomaly of the data grid where the point is located, complete the integration of formula (6) and obtain:

[0023]

[0024] When h = 0, equation (7) is simplified to:

[0025] Δg p0 =(Δg Rp ) (8).

[0026] The approximate integral kernel function is rewritten as the sum of the separated main integral term and the contribution value of the central block. Taking into account the separation effect of the contribution value of the central block, the integral kernel function formula is rewritten as:

[0027]

[0028] When the gravity anomaly of the grid where the calculation point is located changes dramatically and cannot be regarded as a constant value, the additional impact brought about by this is shown in the following formula:

[0029]

[0030] Taking into account the compensation effect of formula (10), formula (9) can be rewritten as:

[0031]

[0032] Formula (11) is applicable to all height segments where h≥0. In actual calculation, when calculating the integral radius s0, it is calculated according to formula (12) and formula (13);

[0033]

[0034] Where, is the geodetic latitude of the calculation point; and Δλ are the data longitude and latitude grid spacings, respectively.

[0035] The method for improving the practicality of the global integration model based on the block separation method includes:

[0036] S201, loading the constructed global integration model based on the block separation method, dividing the global integration domain into a near zone and a far zone, calculating the far zone integral using a high-order potential model, introducing a potential model reference field, and removing the recovery reference field;

[0037] S202, using the truncated form of the Wong-Gore kernel function to remove the corresponding kernel function spherical harmonic expansion of the same order as the potential model reference field to suppress the propagation of observation errors;

[0038] S203: Using a high-order spherical harmonic potential model to compensate for the far-zone integral term, a modified model with no integral singularity is obtained.

[0039] When the far-zone integral is calculated using a high-order potential model, the high-order potential model is expressed as:

[0040]

[0041] Where GM is the product of the universal gravitational constant and the mass of the Earth;

[0042] When removing the restoration reference field:

[0043]

[0044] N is the highest order of the reference field introduced into the potential model; L is the highest order of the ultra-high-order gravity potential model; P n (cosψ) is the Legendre function; To fully normalize the associated Legendre function, and is the fully normalized perturbation coefficient.

[0045] The Wong-Gore kernel function truncation form is used to remove the corresponding kernel function spherical harmonic expansion of the same order as the potential model reference field, and its calculation formula is:

[0046]

[0047] Where K WG (r, ψ) is called the truncated kernel function of K(r, ψ). After partitioning the global integration domain and introducing the potential model reference field and the above truncated kernel function, Equation (11) can be rewritten as:

[0048]

[0049] Where Δg qref is the gravity anomaly of the potential model reference field on the spherical surface; Δg pref is the potential model reference field gravity anomaly at the calculation point outside the sphere; its calculation formulas are:

[0050]

[0051] At this time, Δg in formula (17) p0 Formula (7) is rewritten as:

[0052]

[0053] When the high-order spherical harmonic potential model is used to compensate for the far-zone integral term, the far-zone integral term is compensated using the EGM2008 high-order spherical harmonic potential model.

[0054] Compared with the prior art, the embodiments of the present application have the following beneficial effects:

[0055] In this embodiment of the present invention, planar approximation is applied to both the kernel function and the integral domain, simplifying the computational process. A practical modification scheme is proposed for converting the global integral model for upward continuation of gravity anomalies into a local numerical integral model. The corresponding expressions for the far-zone effect, removal of the restored reference field, and kernel function truncation are presented. Numerical verification of the kernel function singularity solution and the practical modification scheme for the computational model is performed using the ultra-high-order potential model EGM2008. This demonstrates that the proposed scheme for upward continuation of gravity anomalies can achieve an in-model accuracy better than 1 mGal, thus meeting the requirements for high-precision Earth external gravity field assignment.

[0056] In the embodiment of the present invention, the block separation method is used to remove the singular area, and the computational difficulty of ultra-low altitude scenes (near-zone singularity) is converted into a local problem that can be processed analytically, making the method stably applicable to ultra-low altitude scenes below 1 km. By isolating the singular area + local analytical processing, the core strategy fundamentally solves the singularity problem of the integral kernel function in the upward extension of ultra-low altitude gravity anomalies, significantly improving numerical stability, computational accuracy and scenario applicability, while optimizing engineering implementation efficiency. This method provides key technical support for ultra-low altitude scenes such as aerial gravity measurement and low-orbit satellite data processing, and promotes the in-depth application of gravity field detection technology in resource exploration, geological disaster monitoring and other fields.

[0057] In an embodiment of the present invention, when the global integral model based on the block separation method is modified for practicality, a collaborative mechanism of long-wave separation, error suppression, and global compensation is adopted, and the position model reference field is introduced and the recovery reference field is removed. This can effectively separate gravity field signals of different sources and scales, highlight long-wave signals, reduce the interference of the background field on local anomalies, and make the calculation more focused on the characteristics of local gravity anomalies. This solves the singularity, error propagation, and data dependence problems of the traditional Poisson integral model in ultra-low-altitude scenarios, and significantly improves the accuracy, stability, and computational efficiency of the upward extension of gravity anomalies. DETAILED DESCRIPTION

[0058] Unless otherwise defined, all technical and scientific terms used herein have the same meaning as commonly understood by those skilled in the art to which this application belongs; the terms used in the specification of the application herein are only for the purpose of describing specific embodiments and are not intended to limit this application.

[0059] The upward continuation of gravity anomalies is a key method in geophysical exploration data processing. The Poisson integral formula is a classic solution, and its computational stability and efficiency have always been widely recognized. However, when dealing with ultra-low altitude upward continuation problems, the existing methods will cause numerical instability due to the singularity of the integral kernel function when the calculation point and the data point coincide, which affects the calculation accuracy. To address the above problems, we propose a gravity anomaly upward continuation calculation method based on the central grid data block separation method. In short, when the method is implemented, the central data grid where the calculation point is located is first separated from the integral domain based on the block separation method, and the integral kernel function in the small block is processed using a plane approximation analysis. The global integral model based on the central grid data block separation method is output to remove the singularity of the integral kernel function. Then, the global integral model constructed based on the block separation method is loaded, and the global integral domain is divided into a near zone and a far zone. The global integral model based on the block separation method is practically modified to obtain a modified model with the integral singularity removed. In the embodiment of the present invention, by performing plane approximation processing on both the kernel function and the integral domain, the calculation process is relatively simple. At the same time, a practical modification scheme for converting the global integral model of upward extension of gravity anomaly into a local numerical integral model is proposed. The calculation expressions for far-zone effect, removal of recovery reference field and kernel function truncation modification corresponding to the local integral model are given respectively. The kernel function singularity solution and the practical modification scheme of the calculation model are numerically verified using the ultra-high-order potential model EGM2008. It is confirmed that the use of this scheme to realize the upward extension calculation of gravity anomaly can achieve an in-model accuracy better than 1mGal, thus meeting the demand for high-precision assignment of the Earth's external gravity field.

[0060] An embodiment of the present invention provides a method for calculating the upward extension of gravity anomalies based on a central grid data block separation method. The method for calculating the upward extension of gravity anomalies based on a central grid data block separation method specifically includes:

[0061] S10, based on the block separation method, the central data grid where the calculation point is located is separated from the integral domain, and the integral kernel function in the small block is processed using a plane approximate analytical method, and a global integral model is outputted based on the central grid data block separation method to remove the singularity of the integral kernel function;

[0062] In this embodiment of the present invention, a "block separation method" is proposed and implemented for removing singularities in the Poisson integral kernel function during the upward continuation calculation of gravity anomalies. This method removes the central data grid where the calculation point is located from the integral domain and applies a planar approximate analytical expression to its local kernel function. This method avoids the numerical divergence and instability issues that arise at integral singular points in traditional methods. This method is particularly suitable for ultra-low-altitude gravity continuation scenarios below 1 km.

[0063] S20, loading the constructed global integration model based on the block separation method, dividing the global integration domain into a near zone and a far zone, performing practical modification processing on the global integration model based on the block separation method, and obtaining a modified model with removed integration singularities.

[0064] In this embodiment of the present invention, planar approximation is applied to both the kernel function and the integral domain, simplifying the computational process. A practical modification scheme is proposed for converting the global integral model for upward continuation of gravity anomalies into a local numerical integral model. The corresponding expressions for the far-zone effect, removal of the restored reference field, and kernel function truncation are presented. Numerical verification of the kernel function singularity solution and the practical modification scheme for the computational model is performed using the ultra-high-order potential model EGM2008. This demonstrates that the proposed scheme for upward continuation of gravity anomalies can achieve an in-model accuracy better than 1 mGal, thus meeting the requirements for high-precision Earth external gravity field assignment.

[0065] In an embodiment of the present invention, the method for separating the central data grid where the calculation point is located from the integration domain based on the block separation method includes:

[0066] S101, performing a planar approximation process on the integral kernel function within the central grid, and using polar coordinate expansion to approximate the integral kernel function;

[0067] In this embodiment, when the integral kernel function is subjected to planar approximation processing within the central grid:

[0068] The integral kernel function is approximated in the center grid to satisfy r=R+h;R 2 dσ≈sdsdαAt this time, the defined integral kernel function is approximated by polar coordinate expansion:

[0069]

[0070] S102, deriving an integral expression for the approximate integral kernel function and outputting a local analytical integral;

[0071] In this embodiment, when the integral expression is derived for the approximate integral kernel function:

[0072] Assuming that the radius of the central data grid is s0, the contribution of the gravity anomaly of the spherical central data grid to the calculated value of the external gravity anomaly at height h is expressed as:

[0073]

[0074] Within the central data grid, Δg q As a constant, we can take Δg q =ΔgRp , Δg Rp To calculate the known spherical gravity anomaly of the data grid where the point is located, complete the integration of formula (6) and obtain:

[0075]

[0076] When h = 0 (i.e., r = R), equation (7) is simplified to:

[0077] Δg p0 =(Δg Rp ) (8).

[0078] S103 approximates that the integral kernel function is rewritten as the sum of the separated main integral term and the contribution value of the central block. In view of the drastic change in gravity anomaly of the central data grid, an additional correction term is introduced to output a global integral model based on the central grid data block separation method to remove the singularity of the integral kernel function.

[0079] In this embodiment, when the approximate integral kernel function is rewritten as the sum of the separated main integral term and the contribution value of the central block, the integral kernel function formula is rewritten as follows after taking into account the separation effect of the contribution value of the central block:

[0080]

[0081] When the gravity anomaly of the grid where the calculation point is located changes dramatically and cannot be regarded as a constant value, the additional impact brought about by this should also be taken into account. The impact is shown in the following formula:

[0082]

[0083] Taking into account the compensation effect of formula (10), formula (9) should be rewritten as:

[0084]

[0085] Formula (11) is applicable to all height segments where h≥0. In actual calculation, when calculating the integral radius s0, it is calculated according to formula (12) and formula (13);

[0086]

[0087] Where, is the geodetic latitude of the calculation point; and Δλ are the data latitude and longitude grid spacings, respectively. Here, the process of processing the kernel function singularity according to Equation (11) is called the solution for removing the integral kernel function singularity based on the center grid data block separation method. It should be noted that when using this solution, the error caused by the plane approximation of the integral kernel function (see Equation (5)) will also increase as the extension calculation height increases. Therefore, the application scope of this solution must be limited to the ultra-low altitude segment, such as below the 1km extension height. At altitudes above 1km, since the kernel function no longer has singularity problems, the original integral model can restore the normal process to complete accurate and effective numerical calculations. Therefore, the above processing scheme is sufficient to ensure that we obtain stable and reliable gravity anomaly upward extension solution results at all altitudes.

[0088] In the embodiment of the present invention, the block separation method is used to remove the singular area, and the computational difficulty of ultra-low altitude scenes (near-zone singularity) is converted into a local problem that can be processed analytically, making the method stably applicable to ultra-low altitude scenes below 1 km. By isolating the singular area + local analytical processing, the core strategy fundamentally solves the singularity problem of the integral kernel function in the upward extension of ultra-low altitude gravity anomalies, significantly improving numerical stability, computational accuracy and scenario applicability, while optimizing engineering implementation efficiency. This method provides key technical support for ultra-low altitude scenes such as aerial gravity measurement and low-orbit satellite data processing, and promotes the in-depth application of gravity field detection technology in resource exploration, geological disaster monitoring and other fields.

[0089] In a further preferred embodiment of the present invention, a method for improving the practicality of a global integral model based on a block separation method is provided. The method for improving the practicality of a global integral model based on a block separation method specifically includes:

[0090] S201, loading the constructed global integration model based on the block separation method, dividing the global integration domain into a near zone and a far zone, calculating the far zone integral using a high-order potential model, introducing a potential model reference field, and removing the recovery reference field;

[0091] It should be noted that the gravity field is composed of the superposition of long waves and short waves. The traditional method directly uses the measured data integration, which makes the long waves easily contaminated by the measured short wave noise. However, this embodiment introduces the potential model reference field and removes the recovery reference field, which can effectively separate the gravity field signals of different sources and scales, highlight the long wave signal, reduce the interference of the background field on the local anomaly, and make the calculation more focused on the local gravity anomaly characteristics. At the same time, with the help of the high-order potential model to compensate for the far-zone integral term, the result accuracy is further improved.

[0092] S202, using the truncated form of the Wong-Gore kernel function to remove the corresponding kernel function spherical harmonic expansion of the same order as the potential model reference field to suppress the propagation of observation errors;

[0093] In the embodiment of the present invention, when the far-zone integral is calculated using a high-order bit model, the high-order bit model is expressed as:

[0094]

[0095] Where GM is the product of the universal gravitational constant and the mass of the Earth;

[0096] When removing the restoration reference field:

[0097]

[0098] N is the highest order of the reference field introduced into the potential model; L is the highest order of the ultra-high-order gravity potential model; P n (cosψ) is the Legendre function; To fully normalize the associated Legendre function, and is the fully normalized perturbation coefficient.

[0099] S203: Using a high-order spherical harmonic potential model to compensate for the far-zone integral term, a modified model with no integral singularity is obtained.

[0100] When the high-order spherical harmonic potential model is used to compensate for the far-zone integral term, the far-zone integral term is compensated using the EGM2008 high-order spherical harmonic potential model.

[0101] Among them, the Wong-Gore kernel function truncation form is used to remove the corresponding kernel function spherical harmonic expansion of the same order as the potential model reference field, and its calculation formula is:

[0102]

[0103] Where K WG (r, ψ) is called the truncated kernel function of K(r, ψ). After partitioning the global integration domain and introducing the potential model reference field and the above truncated kernel function, Equation (11) can be rewritten as:

[0104]

[0105] Where Δg qref is the gravity anomaly of the potential model reference field on the spherical surface; Δg pref is the potential model reference field gravity anomaly at the calculation point outside the sphere; its calculation formulas are:

[0106]

[0107] At this time, Δg in formula (17) p0 Formula (7) is rewritten as:

[0108]

[0109] In an embodiment of the present invention, when the global integral model based on the block separation method is modified for practicality, a collaborative mechanism of long-wave separation, error suppression, and global compensation is adopted, and the position model reference field is introduced and the recovery reference field is removed. This can effectively separate gravity field signals of different sources and scales, highlight long-wave signals, reduce the interference of the background field on local anomalies, and make the calculation more focused on the characteristics of local gravity anomalies. This solves the singularity, error propagation, and data dependence problems of the traditional Poisson integral model in ultra-low-altitude scenarios, and significantly improves the accuracy, stability, and computational efficiency of the upward extension of gravity anomalies.

[0110] In order to verify the effectiveness of the modified model, the experimental area of ​​Mount Everest, where the gravity anomaly field changes drastically, was selected as the experimental area during the numerical calculation verification and analysis of the model. The specific coverage area is: 6°×6°( 25°N-31°N; λ:84°E-90°E).

[0111] First, the potential model EGM2008 truncated to 360th order is selected as the reference field, that is, N = 360. Then, the truncated potential model EGM2008 of 361~2160th order is selected as the standard field for calculation verification, that is, L = 2160. The model similar to Equation (18) is used to calculate the “true value” of gravity anomaly on the sphere (Δg t ); then select r i =R+h i , R = 6371 km, using the same truncated position model as above, the theoretical “true value” (Δg pti The nine height planes were calculated at altitudes of 0 km, 0.01 km, 0.05 km, 0.1 km, 0.3 km, 1 km, 3 km, 5 km, and 10 km. The coordinates of the calculation points coincide with the data grid points, and each height plane corresponds to 180 × 180 = 32,400 calculation points. Table 1 lists the statistical results of the theoretical "true values" of gravity anomalies on five of these height planes, calculated using the 361-2160 order EGM2008 model. As shown in Table 1, the gravity anomaly in the test area varies dramatically.

[0112] Table 1 Statistics of gravity anomalies on five altitude surfaces calculated by the EGM2008 model (unit: mGal)

[0113] Height / km Maximum Minimum average value RMS value Standard deviation 0 189.40 -162.67 -1.34 53.88 53.87 0.1 186.53 -160.49 -1.32 53.21 53.19 1 164.12 -142.44 -1.18 47.66 47.64 3 127.30 -110.70 -0.91 37.98 37.96 10 65.47 -51.53 -0.38 19.40 19.40

[0114] In order to compare and analyze the calculation effect of the modified model, in this embodiment, the zero-height surface, that is, the 2′×2′ grid gravity anomaly “true value” Δg on the spherical surface is used. tAs the observed quantity, two modified models are used to simultaneously extend the 2′×2′ grid gravity anomaly on the 9 height planes of the corresponding block. Among them, the first model refers to the use of the practical modified gravity anomaly upward extension model for numerical calculation, and the influence of the central data grid that coincides with the calculation point is eliminated for all extended height segments, that is, the central data grid does not participate in the integral summation calculation; the second model refers to the modified model corresponding to Equation (17). This model calculates the contribution Δg of the central data grid to the external gravity anomaly of the sphere according to Equations (20) and (10) at all heights. p0 and Δg p01 ; Compare the calculated values ​​of the two modified models with the theoretical “true value” of gravity anomaly at the corresponding height Δg pti By comparing the different modified models, we can obtain information on the calculation accuracy assessment. The specific comparison results are listed in Table 2. The integration radius is uniformly set to ψ0 = 2°. Considering the influence of the integral edge effect, Table 2 only lists the data comparison results within the 2° × 2° range of the test center area.

[0115] Table 2 Comparison of the upward extension values ​​of gravity anomalies calculated by different modified models and the “true values” (unit: mGal)

[0116]

[0117] From the comparative calculation results given in Table 2, it can be seen that the first model successfully avoids the integral singularity problem by eliminating the influence of the central data grid that coincides with the calculation point at all altitude segments, but this also introduces a non-negligible calculation error. In the ultra-low altitude segment, the root mean square value of this error exceeds 5 mGal. The second model is a simplified version of the modified model that eliminates the integral singularity in this application. Compared with the first model, the calculation accuracy of the second model has been greatly improved. The improvement in calculation accuracy is particularly significant in the ultra-low altitude segment of h ≤ 1 km. As mentioned above, the calculation error of this model increases slightly with the increase of the continuation altitude. This is due to the increase in the planar approximation error of the integral kernel function. In this case, if a segmented calculation model change strategy is adopted, that is, for example, in the altitude segment of h greater than 1 km, the second model is returned to the original model and the normal calculation process (so that the integral singularity modification process is no longer required), then, by combining the calculation results of the second model at altitudes of h ≤ 1 km with the calculation results of the original model at altitudes greater than 1 km, we can obtain stable and reliable continuation calculation results for all altitude segments.

[0118] The “true value” of the gravity anomaly of the potential model Δg as the observed quantity tError interference, ±3mGal random noise and 0.5mGal system deviation were added to the 9 height planes. Then, according to the same calculation scheme and process as before, the modified model was used to complete the external gravity anomaly extension calculation of the 2′×2′ grid sphere on 9 height planes. Finally, the calculation results were compared and evaluated with the “true value” of the position model at the corresponding height. The specific comparison results are shown in Table 3.

[0119] Table 3 Comparison of the calculation results of different modified models under data error interference and the "true value" (unit: mGal)

[0120]

[0121] The statistical results in Table 3 also show that after adding the interference of systematic deviations in data observations, the deviation of the gravity anomaly continuation calculation results from the model “true value” does not change significantly, indicating that under the coupling of random noise interference in data observations and systematic deviations, the combined effect of the error has a certain degree of uncertainty.

[0122] In summary, the present invention provides a method for calculating the upward extension of gravity anomalies based on the central grid data block separation method. In the embodiments of the present invention, planar approximation is applied to both the kernel function and the integration domain, simplifying the calculation process. A practical modification scheme for converting the global integral model for the upward extension of gravity anomalies into a local numerical integral model is also proposed. The corresponding calculation expressions for the far-zone effect, removal of the recovery reference field, and kernel function truncation are given. Numerical verification of the kernel function singularity solution and the practical modification scheme for the calculation model are performed using the ultra-high-order bit model EGM2008. This demonstrates that the proposed scheme for calculating the upward extension of gravity anomalies can achieve an in-model accuracy better than 1 mGal, thus meeting the requirements for high-precision Earth external gravity field assignment.

[0123] It should be noted that for the aforementioned embodiments, for simplicity of description, they are all expressed as a series of action combinations. However, those skilled in the art should be aware that the present invention is not limited by the order of the actions described, because according to the present invention, certain steps may be performed in other orders or simultaneously. Secondly, those skilled in the art should also be aware that the embodiments described in this specification are all preferred embodiments, and the actions and modules involved are not necessarily required by the present invention.

[0124] The above embodiments are only used to illustrate the technical solutions of the present invention, rather than to limit the scope of protection of the invention. Obviously, the embodiments described are only some embodiments of the present invention, rather than all embodiments. Based on these embodiments, all other embodiments obtained by ordinary technicians in this field without making creative work are within the scope of protection of the present invention. Although the present invention has been described in detail with reference to the above embodiments, ordinary technicians in this field can still combine, add, delete or make other adjustments to the features in the various embodiments of the present invention according to the circumstances without conflict, without making creative work, so as to obtain different other technical solutions that do not deviate from the concept of the present invention in essence, and these technical solutions also fall within the scope of protection of the present invention.

Claims

1. A gravity anomaly upward extension calculation method based on the central grid data block separation method is characterized by: The method comprises: S10, based on the block separation method, the central data grid where the calculation point is located is separated from the integral domain, and the integral kernel function in the small block is processed using a plane approximate analytical method, and a global integral model is outputted based on the central grid data block separation method to remove the singularity of the integral kernel function; S20, loading the constructed global integration model based on the block separation method, dividing the global integration domain into a near zone and a far zone, performing practical modification processing on the global integration model based on the block separation method, and obtaining a modified model with removed integration singularities.

2. The method for calculating gravity anomaly upward extension based on the central grid data block separation method according to claim 1, characterized in that: The method for separating the central data grid where the calculation point is located from the integration domain based on the block separation method includes: S101, performing a planar approximation process on the integral kernel function within the central grid, and using polar coordinate expansion to approximate the integral kernel function; S102, deriving an integral expression for the approximate integral kernel function and outputting a local analytical integral; S103 approximates that the integral kernel function is rewritten as the sum of the separated main integral term and the contribution value of the central block. In view of the drastic change in gravity anomaly of the central data grid, an additional correction term is introduced to output a global integral model based on the central grid data block separation method to remove the singularity of the integral kernel function.

3. The method for calculating gravity anomaly upward extension based on the central grid data block separation method according to claim 2, characterized in that: When the integral kernel function is approximated as a plane within the central grid: The integral kernel function is approximated in the center grid to satisfy r=R+h;R 2 dσ≈sdsdαAt this time, the defined integral kernel function is approximated by polar coordinate expansion:

4. The method for calculating gravity anomaly upward extension based on the central grid data block separation method according to claim 3, characterized in that: When deriving the integral expression of the approximate integral kernel function: Assuming that the radius of the central data grid is s0, the contribution of the gravity anomaly of the spherical central data grid to the calculated value of the external gravity anomaly at height h is expressed as: Within the central data grid, Δg q As a constant, Δg q =Δg Rp , Δg Rp To calculate the known spherical gravity anomaly of the data grid where the point is located, complete the integration of formula (6) and obtain: When h = 0, equation (7) is simplified to: Δg p0 =(Δg Rp ) (8)。 5. The method for calculating gravity anomaly upward extension based on the central grid data block separation method according to claim 4, characterized in that: The approximate integral kernel function is rewritten as the sum of the separated main integral term and the contribution value of the central block. Taking into account the separation effect of the contribution value of the central block, the integral kernel function formula is rewritten as: When the gravity anomaly of the grid where the calculation point is located changes dramatically and cannot be regarded as a constant value, the additional impact brought about by this is taken into account. The impact is expressed as follows: Taking into account the compensation effect of formula (10), formula (9) can be rewritten as: Formula (11) is applicable to all height segments where h≥0. In actual calculation, when calculating the integral radius s0, it is calculated according to formula (12) and formula (13); Where, is the geodetic latitude of the calculation point; and Δλ are the data longitude and latitude grid spacings, respectively.

6. The method for calculating gravity anomaly upward extension based on the central grid data block separation method according to claim 5, characterized in that: The method for improving the practicality of the global integration model based on the block separation method includes: S201, loading the constructed global integration model based on the block separation method, dividing the global integration domain into a near zone and a far zone, calculating the far zone integral using a high-order potential model, introducing a potential model reference field, and removing the recovery reference field; S202, using the truncated form of the Wong-Gore kernel function to remove the corresponding kernel function spherical harmonic expansion of the same order as the potential model reference field to suppress the propagation of observation errors; S203: Using a high-order spherical harmonic potential model to compensate for the far-zone integral term, a modified model with no integral singularity is obtained.

7. The method for calculating gravity anomaly upward extension based on the central grid data block separation method according to claim 6, characterized in that: When the far-zone integral is calculated using a high-order potential model, the high-order potential model is expressed as: Where GM is the product of the universal gravitational constant and the mass of the Earth; When removing the restoration reference field: N is the highest order of the reference field introduced into the potential model; L is the highest order of the ultra-high-order gravity potential model; P n (cosψ) is the Legendre function; To fully normalize the associated Legendre function, and is the fully normalized perturbation coefficient.

8. The method for calculating gravity anomaly upward extension based on the central grid data block separation method according to claim 7, characterized in that: The Wong-Gore kernel function truncation form is used to remove the corresponding kernel function spherical harmonic expansion of the same order as the potential model reference field, and its calculation formula is: Where K WG (r, ψ) is called the truncated kernel function of K(r, ψ). After partitioning the global integration domain and introducing the potential model reference field and the above truncated kernel function, Equation (11) can be rewritten as: Where Δg qref is the gravity anomaly of the potential model reference field on the spherical surface; Δg pref Calculate the potential model reference field gravity anomaly at points outside the sphere; The calculation formulas are: At this time, Δg in formula (17) p0 Formula (7) is rewritten as:

9. The method for calculating gravity anomaly upward extension based on the central grid data block separation method according to claim 8, characterized in that: When the high-order spherical harmonic potential model is used to compensate for the far-zone integral term, the far-zone integral term is compensated using the EGM2008 high-order spherical harmonic potential model.

Images