Gravity data downward continuation method, device, medium and product based on discretization model of backward limit Poisson integral

Abstract

The application discloses a gravity data downward continuation method based on a discrete model of a reverse band-limit Poisson integral, equipment, a medium and a product, relates to the technical field of gravity data processing, and performs gridization on discrete gravity point data, determines a grid limit frequency, a truncation order of grid gravity anomaly data and a type of a discrete model of a reverse band-limit Poisson integral, and when a direct discrete model in a removal recovery mode is removed, determines a coefficient matrix, a parameter vector and a constant vector based on the grid gravity anomaly data, a continuation height and an earth gravity field model; determines a target calculation formula based on the coefficient matrix, the parameter vector and the constant vector; when an indirect discrete model in a near-zone integral mode is removed, determines a target calculation formula of the discrete model based on the grid gravity anomaly data and the continuation height; and obtains a downward continuation result based on the target calculation formula of the discrete model. The application can stably and highly accurately perform downward continuation on gravity data.

Description

TECHNICAL FIELD

[0001] The present application relates to the field of gravity data processing, in particular to a gravity data downward continuation method based on a discretization model of a reverse band-limited Poisson integral, equipment, medium and product. BACKGROUND

[0002] Gravity data is an important potential field data in physical geodesy, and its acquisition methods include absolute gravity measurement, relative gravity measurement, satellite altimetry technology, shipborne gravity measurement, airborne gravity measurement, etc. When the measurement height surface of gravity data is inconsistent with the actual application height surface, the gravity data needs to be converted to the required height surface through a continuation algorithm.

[0003] Poisson integral is an important theoretical tool for potential field data continuation and is widely used in practice. Poisson integral is a continuous integral, but since the gravity data being integrated is grid data, the Poisson integral needs to be discretized when downward continuation. Reverse band-limited Poisson integral has filtering properties and stable numerical operation, and has gained more and more attention in the field of physical geodesy. Research has found that the discretization method has an important influence on the downward continuation accuracy of the reverse band-limited Poisson integral. In order to obtain stable and high-precision downward continuation results, it is urgent to develop a discretization method for reverse band-limited Poisson integral that can efficiently continue according to the characteristics of the frequency spectrum information contained in the gravity data. SUMMARY

[0004] The purpose of the present application is to provide a gravity data downward continuation method based on a discretization model of a reverse band-limited Poisson integral, equipment, medium and product, which can stably and accurately downward continue the gravity data.

[0005] To achieve the above purpose, the present application provides the following solutions:

[0006] In a first aspect, the present application provides a gravity data downward continuation method based on a discretization model of a reverse band-limited Poisson integral, which comprises:

[0007] obtaining discrete gravity point data to be continued, grid resolution and continuation height;

[0008] performing grid processing on the discrete gravity point data according to the grid resolution to obtain grid gravity anomaly data;

[0009] determining the grid frequency limit based on the grid resolution, and determining the truncation order of the grid gravity anomaly data based on the grid frequency limit and the discrete gravity point data;

[0010] determine a type of a discretization model of the inverse Poisson integral based on the grid limited frequency and the grid gravity anomaly data; the type of the discretization model is a direct discretization model in a remove-restore mode or an indirect discretization model in a near zone integral mode;

[0011] when the type of the discretization model is the direct discretization model in the remove-restore mode:

[0012] obtain a reference earth gravity field model in the remove-restore mode;

[0013] determine a coefficient matrix of the discretization model based on the grid gravity anomaly data and the continuation height; the coefficient matrix includes diagonal elements and non-diagonal elements;

[0014] determine a parameter vector and a constant vector of the discretization model based on the grid gravity anomaly data, the continuation height and the earth gravity field model;

[0015] determine a target calculation formula of the discretization model based on the coefficient matrix, the parameter vector and the constant vector of the discretization model;

[0016] when the type of the discretization model is the indirect discretization model in the near zone integral mode:

[0017] determine a coefficient matrix and a parameter vector of the discretization model based on the grid gravity anomaly data and the continuation height;

[0018] determine a target calculation formula of the discretization model based on the coefficient matrix and the parameter vector of the discretization model;

[0019] obtain the downward continuation result based on the target calculation formula of the corresponding discretization model.

[0020] In a second aspect, a computer device is provided, which includes a memory, a processor, and a computer program stored in the memory and executable on the processor, and the processor executes the computer program to implement the gravity data downward continuation method based on the inverse Poisson integral discretization model.

[0021] In a third aspect, a computer readable storage medium is provided, which stores a computer program executable on a processor to implement the gravity data downward continuation method based on the inverse Poisson integral discretization model.

[0022] In a fourth aspect, a computer program product is provided, which includes a computer program executable on a processor to implement the gravity data downward continuation method based on the inverse Poisson integral discretization model.

[0023] According to the specific embodiments provided in the application, the application has the following technical effects:

[0024] The application discloses a gravity data downward continuation method based on a discrete model of a band-limited Poisson integral, a device, a medium and a product. First, discrete gravity point data are gridded according to grid resolution; based on the grid resolution, grid frequency limitation is determined, and based on the grid frequency limitation and the discrete gravity point data, the truncation order of grid gravity anomaly data is determined; based on the grid frequency limitation and the truncation order, the type of the discrete model of the band-limited Poisson integral is determined; the type of the discrete model is a direct discrete model in a removal recovery mode or an indirect discrete model in a near-zone integral mode; second, for different discrete models, a target calculation formula of the corresponding discrete model is determined; finally, based on the target calculation formula of the discrete model, a downward continuation result is obtained. According to different situations, the application adopts different discrete models to process the band-limited Poisson integral of gravity data to be continued, so that the precision of the gravity data obtained after downward continuation is significantly improved. BRIEF DESCRIPTION OF DRAWINGS

[0025] In order to more clearly illustrate the technical solutions in the embodiments of the application or the prior art, the following will briefly introduce the drawings needed in the embodiments. Obviously, the drawings in the following description are only some embodiments of the application, and for those skilled in the art, other drawings can also be obtained from these drawings without creative labor.

[0026] Figure 1 An application environment diagram of the gravity data downward continuation method based on the discrete model of the band-limited Poisson integral in an embodiment of the application;

[0027] Figure 2 A flowchart of the gravity data downward continuation method based on the discrete model of the band-limited Poisson integral provided in an embodiment of the application;

[0028] Figure 3 A structural diagram of a computer device provided in an embodiment of the application.

[0029] Reference signs:

[0030] Terminal-102, server-104. DETAILED DESCRIPTION

[0031] With reference to the drawings, the technical solutions in the embodiments of the present application will be fully described below, obviously, the described embodiments are only a part of the embodiments of the present application, rather than all the embodiments. Based on the embodiments in the present application, all other embodiments obtained by those skilled in the art without creative labor are within the scope of protection of the present application.

[0032] In order to make the above-mentioned purposes, features and advantages of the present application more obvious and easy to understand, the present application will be further described in detail below with reference to the drawings and specific embodiments.

[0033] The gravity data downward continuation method based on the discretization model of the reverse band-limit Poisson integral provided by the embodiments of the present application can be applied in the application environment as shown in Figure 1 The terminal 102 communicates with the server 104 through the network. The data storage system can store the data required to be processed by the server 104. The data storage system can be separately arranged, or integrated on the server 104, or placed on the cloud or other servers. The terminal 102 can send the discrete gravity point data to be continued, the grid resolution and the continuation height to the server 104, and the server 104 receives the discrete gravity point data to be continued, the grid resolution and the continuation height. For the discrete gravity point data to be continued, the grid resolution and the continuation height, the server 104 performs gravity data downward continuation based on the discretization model of the reverse band-limit Poisson integral. The server 104 can feed back the obtained downward continuation result to the terminal 102. In addition, in some embodiments, the gravity data downward continuation method based on the discretization model of the reverse band-limit Poisson integral can also be realized by the server 104 or the terminal 102 alone, for example, the terminal 102 can directly perform gravity data downward continuation based on the discretization model of the reverse band-limit Poisson integral for the discrete gravity point data to be continued, the grid resolution and the continuation height, or the server 104 can obtain the discrete gravity point data to be continued, the grid resolution and the continuation height from the data storage system, and perform gravity data downward continuation based on the discretization model of the reverse band-limit Poisson integral for the discrete gravity point data to be continued, the grid resolution and the continuation height.

[0034] The terminal 102 can be, but is not limited to, various desktop computers, notebook computers, smart phones and tablet computers. The server 104 can be realized by an independent server or a server cluster composed of multiple servers, and can also be a cloud server.

[0035] In an exemplary embodiment, as Figure 2As shown, a method for downward continuation of gravity data based on a discretized model using a rewind-limited Poisson integral is provided. This method is executed by a computer device, specifically by a terminal or server alone, or by both a terminal and a server. In this embodiment, the method is applied to... Figure 1 Taking server 104 as an example, the explanation includes the following steps S1 to S13. Wherein:

[0036] Step S1: Obtain the discrete gravity point data, grid resolution, and extension height of the data to be extended.

[0037] Step S2: Grid the discrete gravity point data according to the grid resolution to obtain grid gravity anomaly data.

[0038] Step S3: Based on the grid resolution, determine the grid frequency limit, and based on the grid frequency limit and discrete gravity point data, determine the truncation order of the grid gravity anomaly data.

[0039] Step S4: Determine the type of discretization model for the rewind-limited Poisson integral based on the truncation order of the grid frequency limit and grid gravity anomaly data; the type of discretization model is either a direct discretization model in the removal recovery mode or an indirect discretization model in the near-field integration mode.

[0040] When the discretization model is a direct discretization model under the removal and recovery mode, execute steps S5 to S8:

[0041] Step S5: Obtain the Earth's gravity field model with the recovery mode reference removed.

[0042] Step S6: Based on the grid gravity anomaly data and the extended height, determine the coefficient matrix of the discretized model; the coefficient matrix includes diagonal elements and off-diagonal elements.

[0043] Step S7: Based on grid gravity anomaly data, extended height, and Earth's gravity field model, determine the parameter vector and constant vector of the discretized model.

[0044] Step S8: Based on the coefficient matrix, parameter vector, and constant vector of the discretized model, determine the target calculation formula for the discretized model.

[0045] When the discretization model is an indirect discretization model under near-zone integration mode, execute steps S9 to S10:

[0046] Step S9: Based on the grid gravity anomaly data and the extended height, determine the coefficient matrix and parameter vector of the discretized model.

[0047] Step S10: Based on the coefficient matrix and parameter vector of the discretized model, determine the target calculation formula of the discretized model.

[0048] Step S11, based on the target calculation formula of the corresponding discretization model, the downward continuation result is obtained.

[0049] As an optional implementation, step S3 specifically comprises:

[0050] Step S31, based on the grid resolution, the grid limit frequency is calculated by the formula M' = π / Δ; wherein, M' represents the grid limit frequency; Δ represents the grid resolution in radians.

[0051] Step S32, based on the grid limit frequency and the discrete gravity point data, the truncation order of the grid gravity data is calculated by the formula M = M'·min{n1 / n2,1}; wherein, M represents the truncation order of the grid gravity data; n1 represents the number of discrete gravity points; n2 represents the number of grid gravity data.

[0052] As an optional implementation, step S4 specifically comprises:

[0053] Step S41, it is judged whether the difference between the grid limit frequency and the truncation order of the grid gravity data is greater than or equal to a preset threshold value, and a judgment result is obtained. Specifically, if the grid resolution is 3', the difference threshold value is 900 orders; if the grid resolution is 2', the difference threshold value is 2200 orders, the higher the resolution, the greater the difference threshold value.

[0054] Step S42, if the judgment result is yes, it is determined that the type of the discretization model of the backward limit Poisson integral is the direct discretization model in the remove restore mode.

[0055] Step S43, if the judgment result is no, it is determined that the type of the discretization model of the backward limit Poisson integral is the indirect discretization model in the near zone integral mode.

[0056] As an optional implementation, when the discretization model is the direct discretization model in the remove restore mode, in step S6, the calculation formula of the non-diagonal element in the coefficient matrix of the discretization model is

[0057]

[0058] wherein, represents the element in the i-th row and the j-th column of the coefficient matrix of the direct discretization model in the remove restore mode; s j represents the area of the j-th integral grid, s j = Δ 2 sinθ j , θ j is the co-latitude of the center point of the j-th integral grid, Δ represents the grid resolution in radians; N 1，1N represents the low-order truncation order of the back-limited Poisson kernel function when the discretization model is a direct discretization model in the remove-restore mode 1，1 L, L represents the truncation order of the reference Earth gravity field model; N 1,2 N represents the high-order truncation order of the back-limited Poisson kernel function when the discretization model is a direct discretization model in the remove-restore mode 2,1 M; R represents the mean radius of the Earth; r represents the geocentric radial distance of the calculation point in the spherical approximation; P n represents the Legendre polynomial of order n; ψ ij represents the angular distance between the center of the i-th calculation grid and the center of the j-th integral grid; ψ1 represents the radius of the near-zone integral domain.

[0059] In step S6, the calculation formula of the diagonal elements in the coefficient matrix of the discretization model is

[0060]

[0061] wherein, represents the i-th row diagonal element in the coefficient matrix of the direct discretization model in the remove-restore mode; I n represents the value of the n-th order function with parameter ψ0; ψ0 represents the radius of the circular approximation region of the calculation grid. Wherein, ψ0 can be determined according to the principle that the area of the calculation grid is equal to that of its circular approximation region, and the specific calculation formula is

[0062]

[0063] wherein, θ represents the co-latitude of the calculation point.

[0064] I n The calculation formula of ψ0 is

[0065]

[0066] In step S7, the calculation formula of the parameter vector of the discretization model is

[0067]

[0068] wherein, represents the i-th element in the parameter vector; Δg i represents the i-th grid gravity anomaly data; represents the reference gravity anomaly of the i-th grid calculated by using the Earth gravity field model.

[0069] In step S7, the calculation formula of the constant vector of the discretization model is

[0070]

[0071] where b i represents the i-th element of the constant vector; N represents the highest order of the high-frequency far-zone effect; Ω i represents the spherical coordinate of the i-th calculation point; represents the n-th order gravity anomaly calculated by using the earth gravity field model; represents the n-th order far-zone truncation coefficient of the deconvolution Poisson kernel function.

[0072] where, The calculation formula of is

[0073]

[0074] where K IBL represents the deconvolution Poisson kernel function, and when the discretization model is a direct discretization model in the remove-restore mode, the calculation formula is

[0075]

[0076] In step S8, the target calculation formula of the discretization model is

[0077] Δg High (R) = A IBL1 Δg High (r) + b (9)

[0078] where Δg H igh(R) represents the output vector of the direct discretization model in the remove-restore mode; A IBL1 represents the coefficient matrix of the direct discretization model in the remove-restore mode; Δg High (r) represents the parameter vector; and b represents the constant vector.

[0079] As an optional implementation, when the type of the discretization model is an indirect discretization model in the near-zone integral mode, the calculation formula of the coefficient matrix of the discretization model is

[0080] In step S9, the calculation formula of the non-diagonal element in the coefficient matrix of the discretization model is

[0081]

[0082] where, represents the i-th row and j-th column element in the coefficient matrix of the indirect discretization model in the near-zone integral mode; N 2，1 represents the low-order truncation order of the deconvolution Poisson kernel function when the type of the discretization model is the indirect discretization model in the near-zone integral mode, N 2，1 =-1; N 2,2N represents the high-order truncation order of the deconvolution limited Poisson kernel function when the type of the discretization model is the indirect discretization model under the near-zone integral mode 2,2 = M, i.e. N 2,2 The value of N IBL2 IBL2 The diagonal element in the coefficient matrix of the discretization model in step S9 is calculated by the formula IBL2 IBL2 IBL1 High G wherein, IBL1 The diagonal element in the coefficient matrix of the indirect discretization model under the near-zone integral mode is represented by i; J represents the number of the near-zone grids. G IBL At this time, the parameter vector Δg(r) of the discretization model is constructed by using the gridded gravity anomaly data. n High In step S10, the target calculation formula of the discretization model is High n+2 Δg n (R) = A n+2 Δg(r) (12) n n+2 wherein, Δg n (R) represents the output vector of the indirect discretization model under the near-zone integral mode, A n+2 represents the coefficient matrix of the indirect discretization model under the near-zone integral mode; and Δg(r) represents the parameter vector constructed by the gridded gravity anomaly data. n n+2 As an optional implementation, when the type of the discretization model is the direct discretization model under the remove-restore mode, step S11 specifically includes: n In step S1111, the deconvolution limited Poisson integral processing is performed based on the target calculation formula of the direct discretization model under the remove-restore mode, so as to obtain the output vector of the discretization model.

[0001] In step S1112, the output vector after the reference model value is restored is obtained by using the formula Δg (R) = Δg

[0002] (R) + Δg (R) based on the output vector of the discretization model and the earth gravity field model; wherein, Δg

[0003] (R) represents the output vector after the reference model value is restored; and Δg (R) is the reference gravity anomaly vector on the boundary surface with the radius R calculated by using the earth gravity field model.

[0004] Step S1113, based on the output vector after the reference model value is recovered, eliminate the data of the edge of the calculation area to obtain the downward continuation result.

[0094] As an optional implementation, when the type of the discretization model is an indirect discretization model in the near zone integral mode, step S11 specifically includes:

[0095] Step S1121, based on the target calculation formula of the direct discretization model in the recovery mode, perform the backward limit Poisson integral processing to obtain the output vector of the discretization model;

[0096] Step S1122, based on the output vector of the discretization model, eliminate the data of the edge of the calculation area to obtain the downward continuation result.

[0097] Backward limit Poisson integral

[0098] The backward limit Poisson integral takes the gravity data of a certain height of the measurement surface as the integral physical quantity, and its expression is

[0099]

[0100] Wherein, Ω0 represents the global integral domain, dΩ' is the spherical differential unit, Δg(R, Ω) is the gravity anomaly at the calculation point on the boundary surface, Δg(r, Ω') is the gravity anomaly at the flow integral point on the measurement surface, K IBL represents the backward limit Poisson kernel function, and its calculation formula is

[0101]

[0102] Wherein, N1 is the low-order truncation order, N2 is the high-order truncation order, ψ is the central angle distance between the calculation point and the integral point, and n is the spherical harmonic order.

[0103] Since the boundary surface gravity data can be obtained without inverse operation when the downward continuation is performed by using the backward limit Poisson integral, the downward continuation process based on the backward limit Poisson integral has numerical stability. When the backward limit Poisson integral is applied, discretization processing needs to be performed. If the central grid center point kernel function value is taken as the approximate average value of the kernel function in the grid during discretization, a large discretization error will be caused. In order to avoid a large central area discretization error, a direct or indirect discretization scheme can be used in actual application. The direct discretization method is a calculation mode for directly calculating the central grid integral value by using the analytical expression of the continuous integral in the central grid, and the indirect discretization method is a calculation mode for avoiding direct calculation of the central area integral value by using the modified integral expression that makes the central grid integral value 0.

[0104] The past Poisson integral is usually applied to the non-removal recovery mode, and in fact, the past Poisson integral can also be applied to the non-removal recovery mode, and the non-removal recovery mode can be further divided into a near-zone integral mode and an integral mode considering the influence of a far zone. A discretization method of the past Poisson integral has an important influence on the continuation accuracy. A discretization algorithm of the traditional past Poisson integral is essentially an indirect discretization algorithm containing high-frequency far-zone influence in the non-removal recovery mode. In order to further improve the downward continuation accuracy of the past Poisson integral, the present application first derives direct and indirect discretization formulas of the past Poisson integral in the near-zone integral mode, the integral mode considering the influence of the far zone and the non-removal recovery mode, and then verifies the effectiveness of the discretization algorithm in different calculation modes through numerical experiments.

[0105] Near-zone integral mode

[0106] Direct discretization method:

[0107] Theoretically, the Poisson integral is a global integral, but due to the large amount of calculation of the global integral, the Poisson integral is usually used in the near-zone integral in practical application.

[0108] In order to derive the analytical formula of the past Poisson integral of the central grid, the central grid is approximated as a circular region C0 with the center point of the grid as the origin. The radius ψ0 of C0 can be determined according to the principle that the area of C0 is equal to the area of the central grid (i.e. ). Through derivation, it can be obtained that the diagonal element of the discretization coefficient matrix of the near-zone past Poisson integral when using the direct discretization method is

[0109]

[0110] In the formula, ψ1 represents the radius of the near-zone integral domain. n The calculation formula of I (ψ0) is

[0111]

[0112] The non-diagonal element of the discretization coefficient matrix of the past Poisson integral is

[0113]

[0114] In the formula, ψ1 represents the radius of the near-zone integral domain.

[0115] Indirect discretization method: ​

[0116] The low-order truncation order N1 in the backward continuation Poisson kernel function is an integer greater than or equal to -1. By using the orthogonality of Legendre polynomials, we have

[0117]

[0118] The spherical area integral of the backward continuation Poisson kernel function is

[0119]

[0120] By multiplying the both ends of equation (19) by r / R·Δg(r,Ω) and then subtracting the backward continuation Poisson integral, the modified form of the backward continuation Poisson integral is as follows:

[0121]

[0122] The near-zone backward continuation Poisson integral in the modified form is obtained by changing the integral domain in equation (20) to the near-zone integral domain. According to equation (20), the diagonal elements of the indirect discretization coefficient matrix of the near-zone backward continuation Poisson integral are

[0123]

[0124] where J represents the number of near-zone grids.

[0125] By comparing equations (20) and (13), it can be seen that the non-diagonal elements of the indirect discretization coefficient matrix of the near-zone backward continuation Poisson integral are the same as those of the direct discretization coefficient matrix. In addition, because the near-zone integral mode ignores the far-zone integral value, the direct and indirect discretization formulas under the near-zone integral mode do not contain constant vectors.

[0126] Integral mode considering far-zone influence

[0127] Direct discretization method:

[0128] Because direct far-zone Poisson integral requires far-zone measured gravity data and has a large amount of calculation, the far-zone Poisson integral value is usually calculated by using the Earth's gravity field model in practical applications. The expression of the backward continuation Poisson integral considering far-zone influence is

[0129]

[0130] where C1 represents the near-zone integral domain, L' is the cutoff order of far-zone influence, is the n-order model gravity anomaly at the calculation point.

[0131] It can be seen that the difference between the direct discretization formula of the backward continuation Poisson integral considering far zone effect and the direct discretization formula of the backward continuation Poisson integral in near zone is only that a constant term composed of far zone effect is added.

[0132] Indirect discretization method:

[0133] Through derivation, the modified formula of the backward continuation Poisson integral considering far zone effect is

[0134]

[0135] It can be seen that the diagonal elements of the indirect discretization coefficient matrix of the backward continuation Poisson integral considering far zone effect are

[0136]

[0137] By comparing formula (22) and (23), it can be seen that the constant vector of the indirect discretization formula of the Poisson integral considering far zone effect is the same as the constant vector of the direct discretization formula of the Poisson integral considering far zone effect. In addition, the calculation formula of the non-diagonal elements of the direct and indirect discretization coefficient matrix of the backward continuation Poisson integral considering far zone effect is the same as the calculation formula of the non-diagonal elements of the backward continuation Poisson integral in near zone.

[0138] Remove and restore mode

[0139] Direct discretization method:

[0140] When the remove and restore technique is applied, the low frequency information of the gravity data to be continued is first removed, so that the gravity data to be continued in the Poisson integral only contains high frequency information. In this way, the gravity data obtained by the Poisson integral only contains high frequency information, and the low frequency information can be restored to obtain the final continuation result. Among them, the low frequency information of the gravity data is calculated by the earth gravity field model. The expression of the backward continuation Poisson integral in the remove and restore mode is

[0141]

[0142] In the formula, Δg High (R,Ω) is the high frequency gravity anomaly of the projection position of the calculation point on the boundary sphere, Δg High (r,Ω') is the high frequency gravity anomaly at the flow integral point on the measurement sphere.

[0143] If the low-order truncation order of the high frequency gravity data in the remove and restore mode is lower than the highest order of the gravity field model, the high frequency far zone effect can be further increased. The expression of the backward continuation Poisson integral considering high frequency far zone effect in the remove and restore mode is

[0144]

[0145] where N denotes the highest order of the far-zone high-frequency effects, and L denotes the truncated degree of the reference Earth gravity field model in the remove-restore mode. According to the formulas (25) and (26), the calculation formulas of the diagonal and off-diagonal elements of the direct-discretized coefficient matrix of the back azimuth limited Poisson integral in the remove-restore mode are the same as those of the near-zone back azimuth limited Poisson integral.

[0146] Indirect-discretized method:

[0147] The modified form of the back azimuth limited Poisson integral in the remove-restore mode is

[0148]

[0149] The modified back azimuth limited Poisson integral in the remove-restore mode containing the far-zone high-frequency effects is

[0150]

[0151] According to the formulas (27) and (28), the calculation formulas of the diagonal elements of the indirect-discretized coefficient matrix of the back azimuth limited Poisson integral in the remove-restore mode not containing the far-zone high-frequency effects are the same as those of the near-zone back azimuth limited Poisson integral; and the calculation formulas of the diagonal elements of the indirect-discretized coefficient matrix of the back azimuth limited Poisson integral in the remove-restore mode containing the far-zone high-frequency effects are the same as those of the back azimuth limited Poisson integral considering the far-zone effects. For the off-diagonal elements of the discretized coefficient matrix of the back azimuth limited Poisson integral, the calculation formulas are the same no matter what kind of discretized algorithm and calculation mode is used.

[0152] Numerical experiment

[0153] The actual measured discrete gravity point values can contain full-band gravity information, but the spectral content that can be expressed by the gridded gravity data is limited by the spatial distribution of the discrete gravity points and the highest order number of the spectral domain that can be accommodated by the grid. The highest order number that can be accommodated by the grid is related to the grid resolution. In order to conduct the downward continuation test of the band-limited Poisson integral on gravity data containing different spectral domain information under different resolution conditions, the present application first generates model gravity anomaly data of different resolutions and truncated to different orders on the 0m and 4000m height surfaces using the XGM2019e model. The test area is a mountainous area ranging from 106°E to 109°E and 33°N to 36°N, and the terrain of the test area with a resolution of 2'x2' has a minimum of 365.24m, a maximum of 3423.56m, an average of 1300.84m, and a root mean square value of 1383.16m.

[0154] Table 1 Statistics of model gravity anomaly data truncated to different orders under different resolutions (unit: mGal)

[0155]

[0156]

[0157] According to Table 1, the signal strength of gravity data mainly depends on the truncation order of gravity data. The higher the truncation order, the stronger the signal of gravity data on the same height surface. In the downward continuation test, the model gravity anomaly data on the 4000m height without error and with 2mGal random error are respectively taken as the gravity data to be continued, and the model gravity anomaly data on the 0m height is taken as the check data. In the test, the Poisson integral radius is 1°, and to avoid the influence of edge effect, the data within 1° range from the edge of the test area is removed after continuation. The error of the continuation using the discretized Poisson integral contains two parts: the error of the algorithm itself and the error caused by the measurement noise. The error of the continuation of the gravity data without noise actually reflects the error of the algorithm itself. In order to compare and illustrate the downward continuation effect of the discretized band-limited Poisson integral, the present application first conducts a downward continuation test using the standard Poisson integral considering the far zone influence. Table 2 lists the standard deviation statistics of the continuation error, where M represents the truncation order of the gravity anomaly data (the same below).

[0158] Table 2 Downward continuation error of standard Poisson integral (unit: mGal)

[0159]

[0160] The downward continuation process amplifies the noise of gravity data. According to Table 2, the amplification of the measurement noise is very close when the standard Poisson integral is used to downward continue the gravity data with the same resolution but different truncation orders. This is because the noise of gravity data is amplified to the limit frequency corresponding to the grid resolution when the standard Poisson integral with no filtering property is used to downward continue. According to the Nyquist sampling theorem, the highest order (limit frequency) that the grid gravity data with 5', 4', 3' and 2' resolution can accommodate is 2160, 2700, 3600 and 5400, respectively. The higher the resolution of gravity data, the more serious the amplification of the measurement noise. The experiment in Table 2 shows that the standard Poisson integral is not suitable for downward continuation of the noisy high resolution gravity data. Table 3 lists the standard deviation of the downward continuation error of the discretized cut-off Poisson integral with the non-removal recovery mode for the gravity data with different resolutions. In Table 3, N1 is taken as -1 and N2 is taken as M. The cut-off order of the far zone effect is taken as 360 in the integral mode considering the far zone effect.

[0161] Table 3 Downward continuation error of the cut-off Poisson integral with the non-removal recovery mode (unit: mGal)

[0162]

[0163] It can be seen from Table 3 and Table 2 that the downward continuation accuracy of the cut-off Poisson integral with N2 taken as M for the non-full order gravity data with the measurement noise is obviously higher than that of the standard Poisson integral. This is because the measurement noise of gravity data is amplified to the grid limit frequency when the standard Poisson integral is used for downward continuation, while the noise of gravity data in the downward continuation value obtained by the cut-off Poisson integral with the filtering property is at most amplified to N2 order. The downward continuation accuracy of the cut-off Poisson integral with the direct discretization method for the non-full order gravity data is obviously higher than that for the full order gravity data with the same resolution, which reflects that the diagonal element formula of the direct discretization coefficient matrix is 0.5∑(2n+1)(r / R) n+2 I n The discretization error of the (ψ0) term is very large. The difference between the downward continuation errors of the near zone and the far zone considering integral with the cut-off Poisson integral is very small when the direct discretization method is used, which shows that the (ψ0) term has little effect on the downward continuation solution. The term has little effect on the downward continuation solution. However, the downward continuation accuracy of the cut-off Poisson integral with the indirect discretization method has obvious difference between the two calculation modes, from which it can be inferred that the diagonal element formula of the indirect discretization coefficient matrix is 0.5∑(2n+1)(r / R) n+2 I nThe term (ψ1) has a significant impact on the accuracy of downward continuation. Since the magnitude of this term is independent of the grid resolution, the indirect discretization algorithm for the rewind-limited Poisson integral, which takes into account far-field effects, is not very sensitive to resolution changes in gravity data truncated to the same order. 0.5∑(2n+1)(r / R) n+2 I n The order error of the (ψ1) term increases with the order, reaching its maximum at the highest order (N2). Therefore, when using the indirect discretization method for rewound-limited Poisson integration that considers the far-field effect, the difference in downward continuation error calculated by the kernel function with the same N2 value under different resolution conditions is relatively small. As shown in Table 3, when using the indirect discretized rewound-limited Poisson integration for downward continuation, the continuation accuracy of the near-field integration mode is always significantly higher than that of the integration mode that considers the far-field effect. This phenomenon indicates that when using the indirect discretization method for rewound-limited Poisson integration, the 0.5∑(2n+1)(r / R) of the integration mode that considers the far-field effect... n+2 I n The error of term (ψ1) far exceeds the far-field truncation error caused by ignoring this term in the near-field integration mode. According to the experimental results in Table 3, for non-full-order gravity data with a cutoff order significantly lower than that of grid-limited gravity, the direct discretization formula of the inverted-band-limited Poisson integral, which takes into account the far-field effect, can obtain a high-precision downward continuation solution in the non-removal recovery mode. However, for full-order or near-full-order gravity data, the indirect discretization formula of the near-field inverted-band-limited Poisson integral yields the most ideal continuation effect. Finally, Table 4 statistically analyzes the standard deviation of the continuation error of the inverted-band-limited Poisson integral for discretizing gravity data containing 2mGal random noise under different resolution conditions in the removal recovery mode. Here, N is 720, L is 360, N1 equals L, and N2 is M.

[0164] Table 4. Continuation error of Poisson integral with rewind limit removed in recovery mode (unit: mGal)

[0165]

[0166] As shown in Table 4, the direct discretization method can obtain high-precision downward continuation solution when the non-full rank gravity data with the cut-off order much smaller than the grid limit frequency is subjected to the band-limited Poisson integral in the removal recovery mode, as in the non-removal recovery mode. Considering that the kernel function with N1 equal to L has low-order spectrum leakage, the band-limited Poisson integral with high-frequency far-zone influence should be performed at this time. When the indirect discretization method is used to perform the band-limited Poisson integral in the removal recovery mode, for the gravity data with the cut-off order M equal to 2160, the downward continuation precision with high-frequency far-zone influence is slightly higher than that without high-frequency far-zone influence; for the gravity data with M exceeding 2700, the downward continuation precision without high-frequency far-zone influence is obviously higher than that with high-frequency far-zone influence. This shows that when N2 is large, the error of the (ψ1) term exceeds the influence of the spectrum leakage on the downward continuation precision. n+2 I n The error of the (ψ1) term exceeds the influence of the spectrum leakage on the downward continuation precision. According to the test results in Tables 3 and 4, when the full rank or nearly full rank gravity data is subjected to downward continuation, the indirect discretization formula with the near-zone band-limited Poisson integral has the best continuation effect; and for the non-full rank gravity data with the cut-off order significantly lower than the grid limit frequency, the direct discretization formula with the band-limited Poisson integral considering high-frequency far-zone influence is recommended to be used for downward continuation in the removal recovery mode. According to the test results in Tables 3 and 4, for the 3′ resolution grid, when the difference between the grid limit frequency and the cut-off order of the grid gravity data is greater than or equal to 900 orders, the removal recovery mode should be used, and the direct discretization method considering high-frequency far-zone influence should be used; for the 2′ resolution grid, the difference threshold is about 2200 orders. The higher the grid resolution is, the larger the difference threshold is.

[0167] The traditional discretization algorithm of the band-limited Poisson integral is the indirect discretization formula with high-frequency far-zone influence in the removal recovery mode. According to the test of the present application, when the non-full rank gravity data with 2mGal random noise and truncated to 2700 orders under the condition of 2′ resolution is subjected to downward continuation, the downward continuation precision of the band-limited Poisson integral with high-frequency far-zone influence using the direct discretization method in the removal recovery mode is 68.27% higher than that of the traditional algorithm. For the non-full rank gravity data with 2mGal random noise and truncated to 3600 orders under the condition of 2′ resolution, the downward continuation precision of the near-zone band-limited Poisson integral using the indirect discretization method is 59.86% higher than that of the traditional algorithm. When the full rank gravity data with 2mGal random noise under the conditions of 5′, 4′, 3′ and 2′ resolutions is subjected to downward continuation, the downward continuation precision of the near-zone band-limited Poisson integral using the indirect discretization method is 47.52%, 45.48%, 52.04% and 51.95% higher than that of the traditional algorithm, respectively.

[0168] The application further provides an application scenario of the gravity data downward continuation method based on the discretization model of the reverse band-limited Poisson integral. Specifically, the gravity data downward continuation method based on the discretization model of the reverse band-limited Poisson integral provided in the embodiment can be applied in aerial gravity data processing. Aerial gravity data is obtained, and the aerial gravity data is input into the discretization model based on the reverse band-limited Poisson integral to obtain ground gravity data.

[0169] In an exemplary embodiment, a computer device is provided, which includes a memory, a processor, and a computer program stored in the memory and executable on the processor, and the processor executes the computer program to implement the gravity data downward continuation method based on the discretization model of the reverse band-limited Poisson integral.

[0170] In an exemplary embodiment, a computer device is provided, which can be a server or a terminal, and an internal structure diagram thereof can be as shown in Figure 3 The computer device includes a processor, a memory, an input / output interface (I / O), and a communication interface. The processor, the memory, and the input / output interface are connected through a system bus, and the communication interface is connected to the system bus through the input / output interface. The processor of the computer device is configured to provide computing and control capabilities. The memory of the computer device includes a non-volatile storage medium and an internal memory. The non-volatile storage medium stores an operating system, a computer program, and a database. The internal memory provides an environment for the operating system and the computer program in the non-volatile storage medium to run. The input / output interface of the computer device is configured to exchange information between the processor and external devices. The communication interface of the computer device is configured to communicate with external terminals through a network connection. The computer program is executed by the processor to implement a gravity data downward continuation method based on a discretization model of a reverse band-limited Poisson integral.

[0171] Those skilled in the art can understand that Figure 3 The structure shown in the above

[0172] It should be noted that the user information (including but not limited to user device information, user personal information, etc.) and data (including but not limited to data for analysis, stored data, displayed data, etc.) involved in the present application are all information and data authorized by the user or authorized by all parties, and the collection, use, and processing of the relevant data need to comply with relevant regulations.

[0173] Those skilled in the art can understand that all or part of the processes in the above-mentioned embodiment methods can be completed by instructing the relevant hardware through a computer program. The computer program can be stored in a non-volatile computer readable storage medium, and when executed, can include the processes of the above-mentioned embodiment methods. Any reference to memory, database or other medium used in the embodiments provided in the present application can include at least one of non-volatile and volatile memory. Non-volatile memory can include read-only memory (ROM), magnetic tape, floppy disk, flash memory, optical storage, high-density embedded non-volatile memory, resistive memory (ReRAM), magnetoresistive random access memory (MRAM), ferroelectric memory (FRAM), phase change memory (PCM), graphene memory, etc. Volatile memory can include random access memory (RAM) or external cache memory, etc. As an illustration but not limitation, RAM can be in various forms, such as static random access memory (SRAM) or dynamic random access memory (DRAM), etc.

[0174] The database involved in the embodiments provided in the present application can include at least one of a relational database and a non-relational database. The non-relational database can include a distributed database based on a blockchain, etc., without being limited thereto. The processor involved in the embodiments provided in the present application can be a general-purpose processor, a central processing unit, a graphics processing unit, a digital signal processor, a programmable logic device, a data processing logic device based on quantum computing, etc., without being limited thereto.

[0175] The technical features of the above embodiments can be combined arbitrarily. In order to make the description simple, all possible combinations of the technical features in the above embodiments are not described, however, as long as the combinations of the technical features do not exist contradictory, they should be considered as the scope of the present application.

[0176] The principles and implementation manners of the present application are described in the specific examples in the present application, and the above example descriptions are only used to help understand the method of the present application and its core idea; meanwhile, for those skilled in the art, according to the idea of the present application, the specific implementation manners and application ranges will all have changes. In conclusion, the content of the present specification should not be understood as a limitation on the present application.

Claims

1. A gravity data downward continuation method based on a discretized model of the back-laplace Poisson integral, characterized in that, The gravity data downward continuation method based on the discretization model of the reverse limited Poisson integral comprises: acquiring discrete gravity point data to be continued, grid resolution and continuation height; performing grid processing on the discrete gravity point data according to the grid resolution to obtain grid gravity anomaly data; determining grid frequency limitation based on the grid resolution, and determining the truncation order of the grid gravity anomaly data based on the grid frequency limitation and the discrete gravity point data; determining the type of the discretization model of the reverse limited Poisson integral based on the grid frequency limitation and the truncation order of the grid gravity anomaly data; specifically comprising: judging whether the difference between the grid frequency limitation and the truncation order of the grid gravity data is greater than or equal to a preset threshold to obtain a judgment result; the type of the discretization model is a direct discretization model in a remove restore mode or an indirect discretization model in a near zone integral mode; when the type of the discretization model is the direct discretization model in the remove restore mode: acquiring a reference earth gravity field model in the remove restore mode; determining a coefficient matrix of the discretization model based on the grid gravity anomaly data and the continuation height; the coefficient matrix comprises diagonal elements and non-diagonal elements; determining a parameter vector and a constant vector of the discretization model based on the grid gravity anomaly data, the continuation height and the earth gravity field model; determining a target calculation formula of the discretization model based on the coefficient matrix, the parameter vector and the constant vector of the discretization model; when the type of the discretization model is the indirect discretization model in the near zone integral mode: determining the coefficient matrix and the parameter vector of the discretization model based on the grid gravity anomaly data and the continuation height; determining the target calculation formula of the discretization model based on the coefficient matrix and the parameter vector of the discretization model; obtaining a downward continuation result based on the target calculation formula of the corresponding discretization model; when the type of the discretization model is the direct discretization model in the remove restore mode, obtaining the downward continuation result based on the target calculation formula of the corresponding discretization model, specifically comprising: performing reverse limited Poisson integral processing based on the target calculation formula of the direct discretization model in the remove restore mode to obtain an output vector of the discretization model; Based on the output vector of the discretized model and the earth gravity field model, the formula Δg IBL1 (R) = Δg High (R) + Δg G (R) is used to obtain the output vector after the reference model value is recovered; wherein Δg IBL1 (R) represents the output vector after the reference model value is recovered; Δg High (R) represents the output vector after the direct discretized model in the recovery mode is removed; Δg G (R) is the reference gravity anomaly vector on the boundary surface with a radius of R calculated by the earth gravity field model. eliminating data at the edge of a calculation region based on the output vector after the reference model value is restored to obtain the downward continuation result.

2. The gravity data downward continuation method based on the discretized model of the deconvolution limited Poisson integral of claim 1, wherein, determining the grid frequency limitation based on the grid resolution, and determining the truncation order of the grid gravity anomaly data based on the grid frequency limitation and the discrete gravity point data, specifically comprising: calculating the grid frequency limitation by using a formula M' = π / Δ based on the grid resolution; wherein M' represents the grid frequency limitation; and Δ represents the grid resolution in radians; calculating the truncation order of the grid gravity data by using a formula M = M'·min{n1 / n2, 1} based on the grid frequency limitation and the discrete gravity point data; wherein M represents the truncation order of the grid gravity data; n1 represents the number of the discrete gravity points; and n2 represents the number of the grid gravity data.

3. The gravity data downward continuation method based on the discretized model of the deconvolution limited Poisson integral of claim 2, wherein, determining the type of the discretization model of the reverse limited Poisson integral based on the grid frequency limitation and the truncation order of the grid gravity anomaly data, further comprises: If the judgment result is yes, the type of the discretization model of the backward continuation Poisson integral is determined as a direct discretization model in the remove restore mode; If the judgment result is no, the type of the discretization model of the backward continuation Poisson integral is determined as an indirect discretization model in the near zone integral mode.

4. The gravity data downward continuation method based on the discretized model of the deconvolution limited Poisson integral of claim 3, wherein, When the discretization model is the direct discretization model in the remove restore mode: The calculation formula of the non-diagonal elements in the coefficient matrix of the discretization model is wherein denotes the element in the i-th row and j-th column of the coefficient matrix of the direct discretization model in the remove-restore mode; s j denotes the area of the j-th integration grid; N 1，1 denotes the low-order truncation order of the back-to-back Poisson kernel function when the discretization model is the direct discretization model in the remove-restore mode, N 1，1 = L, L denotes the truncation order of the reference Earth gravity field model; N 1,2 denotes the high-order truncation order of the back-to-back Poisson kernel function when the discretization model is the direct discretization model in the remove-restore mode, N 1,2 = M; R denotes the mean radius of the Earth; r denotes the geocentric radius of the calculation point in the spherical approximation; P n denotes the Legendre polynomial of order n; ψ ij denotes the angular distance between the center point of the i-th calculation grid and the center point of the j-th integration grid; ψ1 denotes the radius of the near-zone integration domain; The calculation formula of the diagonal elements in the coefficient matrix of the discretization model is wherein denotes the diagonal element of the i-th row of the coefficient matrix of the direct discretization model in the recovery mode; I n (ψ0) denotes the value of the n-th order function with parameter ψ0; ψ0denotes the radius of the circular approximation region of the calculation grid; The calculation formula of the parameter vector of the discretization model is wherein, denotes the i-th element of the parameter vector; Δg i denotes the i-th grid gravity anomaly data; denotes the reference gravity anomaly of the i-th grid calculated using the Earth's gravitational field model; The calculation formula of the constant vector of the discretization model is where b i represents the i-th element of the constant vector; N represents the highest order of the high-frequency far-zone effect; Ω i represents the spherical coordinates of the i-th calculation point; represents the n-order gravity anomaly calculated using the Earth's gravity field model; represents the n-order far-zone truncation coefficient of the back-laplace limited Poisson kernel function; The target calculation formula of the discretization model is Δg High (R) = A IBL1 Δg High (r) + b; where Δg High (R) denotes the output vector of the direct discretized model in the remove and recover mode; A IBL1 denotes the coefficient matrix of the direct discretized model in the remove and recover mode; Δg High (r) denotes the parameter vector; b denotes the constant vector.

5. The gravity data downward continuation method based on the discretized model of the deconvolution limited Poisson integral of claim 3, wherein, When the type of the discretization model is the indirect discretization model in the near zone integral mode: The calculation formula of the non-diagonal elements in the coefficient matrix of the discretization model is wherein, represents the element in the i-th row and j-th column of the coefficient matrix of the indirect discretization model under the near-zone integral mode; R represents the mean radius of the Earth; r represents the geocentric radial distance of the calculation point under the spherical approximation; s j represents the area of the j-th integral grid; N 2，1 represents the low-order truncation order of the deconvolution-limited Poisson kernel function when the type of the discretization model is the indirect discretization model under the near-zone integral mode, N 2，1 = -1; N 2,2 represents the high-order truncation order of the deconvolution-limited Poisson kernel function when the type of the discretization model is the indirect discretization model under the near-zone integral mode, N 2,2 = M; P n represents the n-th Legendre polynomial; ψ ij represents the angular distance between the center point of the i-th calculation grid and the center point of the j-th integral grid; ψ1 represents the radius of the near-zone integral domain; The calculation formula of the diagonal elements in the coefficient matrix of the discretization model is wherein represents the diagonal element of the i-th row of the coefficient matrix of the indirect discretization model in the near zone integration mode; J represents the number of near zone grids; represents the element of the i-th row and j-th column of the discretization coefficient matrix of the rewind limit Poisson integral; The target calculation formula of the discretization model is Δg IBL2 (R) = A IBL2 Δg(r); where Δg IBL2 (R) represents the output vector of the indirect discretization model under the near-zone integral mode, A IBL2 represents the coefficient matrix of the indirect discretization model under the near-zone integral mode; and Δg(r) represents the parameter vector constructed by the grid gravity anomaly data.

6. The gravity data downward continuation method based on the discretized model of the deconvolution limited Poisson integral of claim 1, wherein, When the type of the discretization model is the indirect discretization model in the near zone integral mode, the downward continuation result is obtained based on the target calculation formula of the corresponding discretization model, and specifically includes: Based on the target calculation formula of the direct discretization model in the remove restore mode, the backward continuation Poisson integral processing is performed to obtain the output vector of the discretization model; Based on the output vector of the discretization model, the data at the edge of the calculation region is eliminated to obtain the downward continuation result.

7. A computer device comprising: A memory, a processor and a computer program stored in the memory and executable on the processor, characterized in that the processor executes the computer program to implement the gravity data downward continuation method based on the discretization model of the backward continuation Poisson integral according to any one of claims 1-6.

8. A computer-readable storage medium having stored thereon a computer program, characterized in that, The computer program is executed by the processor to implement the gravity data downward continuation method based on the discretization model of the backward continuation Poisson integral according to any one of claims 1-6.

9. A computer program product comprising a computer program, characterized in that, The computer program is executed by the processor to implement the gravity data downward continuation method based on the discretization model of the backward continuation Poisson integral according to any one of claims 1-6.

Images