Optimization calculation method for downward continuation of airborne gravity data

Abstract

The invention relates to an optimization calculation method for downward continuation of airborne gravity data, which comprises the following steps of: calculating a Poisson kernel truncation coefficient according to a rigorous method, solving and correcting unknown parameters by utilizing the truncation coefficient, and correcting a Poisson kernel function; calculating the influence of a far region, and discretizing a near region to form an observation equation; and performing regularization processing on the equation, and calculating a least square solution based on singular value decomposition. By adopting the technical scheme of the invention, noise interference can be effectively suppressed through kernel function correction and regularization processing, and the stability of a continuation result is improved; and the least square solution is calculated by adopting singular value decomposition, so that the calculation efficiency is improved.

Description

Technical Field

[0001] The invention relates to the technical field of aviation gravity measurement, in particular to an optimization calculation method for downward extension of aviation gravity data. Background Art

[0002] In geodesy, the Earth's gravity field is used to determine the Earth's shape and elevation benchmark. Gravity measurement is a discrete sampling of the continuous gravity field. Aerial gravity measurement has obvious advantages in difficult areas (such as mountainous areas and polar regions) that are difficult to reach by ground gravity measurement, as well as in gravity blank areas at the junction of land and sea. It can quickly obtain medium and short wave gravity information. Aerial gravity data is the gravity value on the flight survey line. For most geodetic applications, it is usually necessary to extend the gravity data in the air downward to the ground or other reference surfaces for multi-source gravity data fusion or direct calculation of the geoid.

[0003] The downward continuation of aviation gravity is an ill-posed inverse problem. The continuation process amplifies the observation error, and it is difficult to obtain a stable and reliable continuation solution. At present, the main method to solve this problem is to solve the inverse Poisson integral equation. Based on the Poisson integral equation, an iterative continuation model can be obtained after transformation, but the iteration causes the observation noise to accumulate continuously, and it is replaced by the least squares continuation model of direct inversion.

[0004] In order to improve the stability of the extension, the "removal-restore" of the potential coefficient model is used to change the spectral characteristics of the gravity data. By correcting the Poisson kernel function, the distortion of long waves can be avoided, the influence of the far-field effect can be reduced, and the accuracy of the extension calculation can be improved. Commonly used modified Poisson kernel function methods include the spherical Poisson kernel based on the WG idea and the modified spherical Poisson kernel based on the Molodensky method. The truncation coefficient of the spherical Poisson kernel is greater than that of the modified spherical Poisson kernel, so the modified spherical Poisson kernel has more advantages.

[0005] Airborne gravity measurement data is affected by various errors. When the data interval is small or the extended height is large, the normal equation is ill-conditioned and regularization is required during the least squares solution. Selecting an appropriate regularization parameter is the key to the entire regularization algorithm. If the regularization parameter is too small, the data observation noise cannot be effectively suppressed. If the regularization parameter is too large, the solution will be too smooth. In order to find the optimal regularization parameter, it is usually necessary to loop within a certain range. However, this process involves large matrix multiplication and inversion operations, which greatly limits the computational efficiency.

[0006] In summary, due to the ill-conditioning of the downward continuation observation equation of aviation gravity, when the data interval is small or the continuation height is large, the degree of ill-conditioning is aggravated, and it is difficult to obtain a stable and reliable continuation solution. In addition, the calculation of Poisson kernel truncation coefficients usually adopts numerical integration method. When the integration step size is reduced, the numerical calculation accuracy is improved, but the calculation speed is slowed down. Summary of the invention

[0007] The purpose of the present invention is to provide an optimization calculation method for downward extension of aerial gravity data, to establish a rigorous calculation method for Poisson kernel cutoff coefficients, and an optimization method for downward extension of aerial gravity based on kernel correction, regularization, and singular value decomposition, so as to improve the accuracy of downward extension data of aerial gravity while taking into account the calculation efficiency.

[0008] To achieve the above object, the present invention provides the following technical solutions:

[0009] The optimization calculation method for downward extension of aerial gravity data includes the following steps:

[0010] Step 1: Calculate the Poisson kernel cutoff coefficient according to the rigorous formula of Poisson kernel cutoff coefficient;

[0011] Step 2: Use the Poisson kernel cutoff coefficient to solve and correct the unknown parameters, modify the Poisson kernel function, calculate the influence of the far area, and discretize the near area to form the observation equation;

[0012] Step 3: Regularize the observation equation and calculate the extension result based on singular value decomposition.

[0013] In the above method, the implementation method of step one is:

[0014] First use formula (1) to calculate the intermediate quantity R n,k (t) and R n,n (t):

[0015]

[0016] Where ψ0 represents the integral radius of the spherical cap, let t = cosψ0 represent the cosine of the integral radius of the spherical cap, P n (t) represents the nth-order Legendre polynomial, P k (t) represents the k-order Legendre polynomial, n represents the order of the Poisson kernel cutoff coefficient to be calculated, and k represents the loop index, ranging from 0 to n;

[0017] Then use the recursive formula (4) to calculate the intermediate quantity U n (t,h) and V n (t,h):

[0018] The initial values ​​are formula (2) and formula (3):

[0019]

[0020] The recursive formula (4) is:

[0021]

[0022] In the formula, h = R / r, R is the radius of the earth, and r is the diameter from the center of the earth at the calculation point;

[0023] Finally, the Poisson kernel cutoff coefficient Q is directly calculated according to formula (5): n (r,ψ0,R):

[0024]

[0025] The implementation method of step 2 is:

[0026] First, calculate the modified Poisson kernel function: According to formula (6), L+1 linear equations can be formed. Solve the equations to calculate the modified unknown parameter t m (r,ψ0,R):

[0027]

[0028] Where L represents the correction order, n represents the Poisson kernel truncation coefficient order to be calculated, ranging from 0 to L, and m represents the correction unknown parameter index, ranging from 0 to L;

[0029] Modified Poisson kernel function K m (r,ψ,R) is expressed as formula (7)

[0030]

[0031] Where ψ represents the spherical angular distance between the calculation point and the flow point;

[0032] Then, calculate the far zone effects:

[0033] Modified Poisson kernel cutoff coefficient It can be expressed as formula (8):

[0034]

[0035] Gravity anomaly nth-order Laplace harmonic function Δg n (θ,λ) is calculated by the bit coefficient model formula (9):

[0036]

[0037] Where GM is the Earth's gravitational constant, P nm(cosθ) is the fully normalized associated Legendre function, and is the fully normalized perturbation coefficient;

[0038] The series form of the calculation of far-zone impact is formula (10):

[0039]

[0040] Where N max represents the maximum order of the model, θ represents the pole distance of the calculation point, and λ represents the longitude of the calculation point;

[0041] Finally, the observation equation is formed:

[0042] The off-diagonal element A in the i-th row and j-th column of the coefficient matrix ij and the diagonal element A in row i and column i ii They are formula (11) and formula (12) respectively:

[0043]

[0044] in,

[0045]

[0046] Where, d m (r i ,ψ0,R) represents the intermediate temporary variable, N c Indicates the number of points to be found within the integral radius;

[0047] The observation equation is formula (14):

[0048]

[0049] In the formula, Δg i (r,θ,λ) represents the airborne gravity observation value, Δg j (R,θ′,λ′) represents the extended value to be calculated,

[0050] represents the far zone influence, and N represents the number of extended values ​​to be calculated;

[0051] Assuming there are M airborne gravity observations, then according to formula (14), the matrix form is formed, that is, formula (15):

[0052] b=Ax (15)

[0053] Where b represents the column vector formed by the M airborne gravity observations on the left side of formula (14) minus the influence of the distant area, A represents the coefficient matrix of M rows and N columns calculated according to formulas (11), (12) and (13), and x represents the column vector composed of N extended values ​​to be calculated;

[0054] The implementation method of step three is:

[0055] First, use formula (16) to perform singular value decomposition (SVD) on the coefficient matrix A:

[0056] A=UΛV T (16)

[0057] Where U and V are composed of the left and right eigenvectors respectively, U = [u1 u2…u M ],V=[v1 v2…v N ],Λ=diag(λ1,λ2,…,λ N ) is a diagonal matrix, λ i is a singular value;

[0058] Then, the generalized cross validation GCV method is used to determine the regularization parameter. The range and step size of the regularization parameter α are set, and the traversal calculation is performed. The regularization solution of the downward extension of the aerial gravity data is converted into formula (17):

[0059]

[0060] In the formula, represents the uncalculated extension value when the regularization parameter is α, u i represents the i-th left eigenvector;

[0061] Use formula (18) to calculate the intermediate quantity F α

[0062]

[0063] Judgment Criteria G α Formula (19):

[0064]

[0065] Take the regularization parameter corresponding to the minimum value of formula (19) as the final regularization parameter, and substitute it into formula (17) to calculate the least squares extension solution.

[0066] The beneficial effects of the present invention are:

[0067] Aiming at the ill-posed problem of downward extension of aeronautical gravity, the present invention establishes a rigorous calculation method for Poisson kernel truncation coefficient, which is used to calculate and correct unknown parameters and far-zone influences. The downward extension of aeronautical gravity is optimized based on three methods: kernel correction, regularization, and singular value decomposition, which significantly improves the stability and accuracy of the extension results and enhances the calculation efficiency. BRIEF DESCRIPTION OF THE DRAWINGS

[0068] Figure 1It is a schematic diagram of the working process of the present invention. DETAILED DESCRIPTION

[0069] In order to make the purpose, technical solution and advantages of the present invention more clearly understood, the present invention is further described in detail below in conjunction with the accompanying drawings. It should be understood that the specific embodiments described herein are only used to explain the present invention and are not used to limit the present invention.

[0070] refer to Figure 1 The purpose of the present invention is to provide an optimization calculation method for downward extension of aviation gravity data, comprising the following steps:

[0071] Step 1: Calculate the Poisson kernel cutoff coefficient according to the direct calculation formula of the Poisson kernel cutoff coefficient;

[0072] Step 2: Use the Poisson kernel cutoff coefficient to solve and correct the unknown parameters, modify the Poisson kernel function, calculate the influence of the far area, and discretize the near area to form the observation equation;

[0073] Step 3: Regularize the observation equation and calculate the extension result based on singular value decomposition.

[0074] The implementation method of step one is:

[0075] The conventional method for calculating the Poisson kernel truncation coefficient is numerical integration, but it affects the calculation efficiency. By establishing a direct calculation method for the Poisson kernel truncation coefficient, the calculation results can be obtained quickly without losing accuracy.

[0076] First use formula (1) to calculate the intermediate quantity R n,k (t) and R n,n (t):

[0077]

[0078] Where ψ0 represents the integral radius of the spherical cap, let t = cosψ0 represent the cosine of the integral radius of the spherical cap, P n (t) represents the nth-order Legendre polynomial, P k (t) represents the k-order Legendre polynomial, n represents the order of the Poisson kernel cutoff coefficient to be calculated, and k represents the loop index, ranging from 0 to n;

[0079] Then use the recursive formula (4) to calculate the intermediate quantity U n (t,h) and V n (t,h):

[0080] The initial values ​​are formula (2) and formula (3):

[0081]

[0082] The recursive formula (4) is:

[0083]

[0084] In the formula, h = R / r, R is the radius of the earth, and r is the diameter from the center of the earth at the calculation point;

[0085] Finally, the Poisson kernel cutoff coefficient Q is directly calculated according to formula (5): n (r,ψ0,R):

[0086]

[0087] The implementation method of step 2 is:

[0088] In actual calculations, the observed data are discrete values ​​within a certain area, and the integration area is generally divided into two parts: the near area and the far area. In order to improve the stability of the continuation, when performing the downward continuation calculation, the low-order model gravity anomaly can be removed first, the remaining gravity anomaly can be extended, and then the model gravity anomaly can be restored. In order to avoid the distortion of the long wave by the Poisson kernel function, the Poisson kernel function needs to be corrected. The corrected Poisson kernel function can further reduce the influence of the far-zone effect and improve the accuracy of the continuation calculation.

[0089] First, calculate the modified Poisson kernel function: According to formula (6), L+1 linear equations can be formed. Solve the equations to calculate the modified unknown parameter t m (r,ψ0,R):

[0090]

[0091] Where L represents the correction order, n represents the Poisson kernel truncation coefficient order to be calculated, ranging from 0 to L, and m represents the correction unknown parameter index, ranging from 0 to L;

[0092] Modified Poisson kernel function K m (r,ψ,R) is expressed as formula (7)

[0093]

[0094] Where ψ represents the spherical angular distance between the calculation point and the flow point;

[0095] Then, calculate the far zone effects:

[0096] Modified Poisson kernel cutoff coefficient It can be expressed as formula (8):

[0097]

[0098] Gravity anomaly nth-order Laplace harmonic function Δg n (θ,λ) is calculated by the bit coefficient model formula (9):

[0099]

[0100] Where GM is the Earth's gravitational constant, P nm (cosθ) is the fully normalized associated Legendre function, and is the fully normalized perturbation coefficient;

[0101] The series form of the calculation of far-zone impact is formula (10):

[0102]

[0103] Where N max represents the maximum order of the model, θ represents the pole distance of the calculation point, and λ represents the longitude of the calculation point;

[0104] Finally, the observation equation is formed:

[0105] The off-diagonal element A in the i-th row and j-th column of the coefficient matrix ij and the diagonal element A in row i and column i ii They are formula (11) and formula (12) respectively:

[0106]

[0107] in,

[0108]

[0109] Where, d m (r i ,ψ0,R) represents the intermediate temporary variable, N c Indicates the number of points to be found within the integral radius;

[0110] The observation equation is formula (14):

[0111]

[0112] In the formula, Δg i (r,θ,λ) represents the airborne gravity observation value, Δg j (R,θ′,λ′) represents the extended value to be calculated,

[0113] represents the far zone influence, and N represents the number of extended values ​​to be calculated;

[0114] Assuming there are M airborne gravity observations, then according to formula (14), the matrix form is formed, that is, formula (15):

[0115] b=Ax (15)

[0116] Where b represents the column vector formed by the M airborne gravity observations on the left side of formula (14) minus the influence of the distant area, A represents the coefficient matrix with M rows and N columns calculated according to formulas (11), (12) and (13), and x represents the column vector composed of N extended values ​​to be calculated.

[0117] The implementation method of step three is:

[0118] Normally, the observation equation of downward extension is ill-conditioned, and the extension process will amplify the observation error and unmodeled signal. Regularizing the equation can suppress the influence of various errors on the extended estimate. However, the regularized solution process faces multiple large matrix multiplication and inversion operations, which greatly affects the calculation efficiency. The spectral decomposition form can avoid large matrix multiplication and inversion operations, saving calculation time.

[0119] First, use formula (16) to perform singular value decomposition SVD on the coefficient matrix A:

[0120] A=UΛV T (16)

[0121] Where U and V are composed of the left and right eigenvectors respectively, U = [u1 u2…u M ],V=[v1 v2…v N ],Λ=diag(λ1,λ2,…,λ N ) is a diagonal matrix, λ i is a singular value;

[0122] Then, the generalized cross validation GCV method is used to determine the regularization parameter: the regularization parameter α is set to select the range and step size, and the traversal calculation is performed. The regularization solution of the downward extension of the aerial gravity data is converted into formula (17):

[0123]

[0124] In the formula, represents the uncalculated extension value when the regularization parameter is α, u i represents the i-th left eigenvector;

[0125] Use formula (18) to calculate the intermediate quantity F α

[0126]

[0127] Judgment Criteria G α Formula (19):

[0128]

[0129] Take the regularization parameter corresponding to the minimum value of formula (19) as the final regularization parameter, and substitute it into formula (17) to calculate the least squares extension solution.

[0130] The gravity data simulated by the gravity field potential coefficient model are used for test analysis. The 3°×3° range (28°N—31°N, 86°E—89°E) in the Qinghai-Tibet Plateau is taken as the study area, with an average extension height of 5km. The 2190th order of the EIGEN-6C4 potential model is used to simulate and calculate the 2′×2′ resolution grid gravity anomaly data at the ground and 10km flight altitude, and the standard deviation of 2mGal random noise is used to simulate the measurement noise. Experiments are carried out according to the two types of gravity anomaly data without error and adding 2mGal error. The least squares adjustment solution, regularized solution, and regularized solution of the Poisson kernel function modified by removing the 360th order model value are used in turn for downward extension. The regularized solution is calculated according to the singular value decomposition, and the integral radius is selected as 1°. The ground gravity anomaly value is used as a reference for comparison. In order to reduce the edge effect, the range of the statistical results is reduced by 0.5° from the boundary of the study area to the inside.

[0131] The far-zone influence of the study area is calculated, and the statistical results are shown in Table 1. It can be seen that the far-zone influence of the uncorrected Poisson kernel function is significant and exceeds the accuracy of airborne gravity measurement. This influence should be considered in the continuation calculation; while the far-zone influence of the corrected Poisson kernel function is significantly reduced and can be ignored in the calculation.

[0132] Table 1 Far zone impact statistics (unit: mGal)

[0133] Correction Type average value Standard Deviation Minimum Maximum Uncorrected 1.495 0.493 -0.290 2.370 Correction 0.000 0.001 -0.003 0.003

[0134] Table 2 shows the difference statistics of the ground gravity anomaly minus the downward continuation calculation results. It can be seen that without random noise, the continuation accuracy of the three methods does not change much, and the standard deviation is close to 1mGal. The modified kernel regularization method improves the regularization solution by 35% compared with the uncorrected one. When 2mGal random noise is added, the least squares solution becomes extremely unstable, and the difference change even exceeds the gravity anomaly itself. After regularization, the standard deviation is reduced from 97.888mGal to 7.006mGal, and the result is relatively stable. On this basis, the accuracy is improved by 18% after the Poisson kernel function is corrected. In general, the accuracy is reduced after adding random noise.

[0135] Table 2 Statistics of the difference between ground gravity anomaly and downward extension value (unit: mGal)

[0136]

[0137] In order to analyze the computational time of the general regularization solution method and the singular value decomposition method using the matrix inversion method, the computational time of the regularization solution of the uncorrected integral kernel using the inversion method and the SVD decomposition method is recorded, and the results are shown in Table 3. It can be seen that the two methods use the same set of equations for regularization calculations, and the two calculations of the set of equations are close in time. The main difference lies in the loop traversal of the regularization parameter. The calculation speed of the SVD method is 16 times that of the inversion method.

[0138] Table 3 Regularization calculation time statistics (unit: minutes)

[0139] Calculation method Group equation Regularization solution Total time Inverse 0.035 65.454 65.494 SVD 0.036 3.921 3.959

[0140] The above experiments selected difficult areas with complex terrain, which pose great challenges to the downward extension of aerial gravity data. The experimental results further prove the effectiveness of the technical solution provided by the present invention.

[0141] What is described above is only a preferred specific implementation manner of the present invention, but the protection scope of the present invention is not limited thereto. Any technician familiar with the technical field can make equivalent replacements or changes according to the technical scheme and inventive concept of the present invention within the technical scope disclosed by the present invention, which should be covered by the protection scope of the present invention.

Claims

1. An optimization calculation method for downward extension of aerial gravity data, characterized in that: The following steps are involved: Step 1: Calculate the Poisson kernel cutoff coefficient according to the rigorous formula of Poisson kernel cutoff coefficient; Step 2: Use the Poisson kernel cutoff coefficient to solve and correct the unknown parameters, modify the Poisson kernel function, calculate the influence of the far area, and discretize the near area to form the observation equation; Step 3: Regularize the observation equation and calculate the extension result based on singular value decomposition.

2. The optimization calculation method for downward extension of aerial gravity data according to claim 1 is characterized in that: The implementation method of step one is: First use formula (1) to calculate the intermediate quantity R n,k (t) and R n,n (t): Where ψ0 represents the integral radius of the spherical cap, let t = cosψ0 represent the cosine of the integral radius of the spherical cap, P n (t) represents the nth-order Legendre polynomial, P k (t) represents the k-order Legendre polynomial, n represents the order of the Poisson kernel cutoff coefficient to be calculated, and k represents the loop index, ranging from 0 to n; Then use the recursive formula (4) to calculate the intermediate quantity U n (t,h) and V n (t,h): The initial values ​​are formula (2) and formula (3): The recursive formula (4) is: In the formula, h = R / r, R is the radius of the earth, and r is the diameter from the center of the earth at the calculation point; Finally, the Poisson kernel cutoff coefficient Q is directly calculated according to formula (5): n (r,ψ0,R): The implementation method of step 2 is: First calculate the modified Poisson kernel function: Formula (6) can form L+1 linear equations, and the unknown parameters t can be corrected by solving the equations. m (r,ψ0,R): Where L represents the correction order, n represents the Poisson kernel truncation coefficient order to be calculated, ranging from 0 to L, and m represents the correction unknown parameter index, ranging from 0 to L; Modified Poisson kernel function K m (r,ψ,R) is expressed as formula (7): Where ψ represents the spherical angular distance between the calculation point and the flow point; Then calculate the far zone impact: Modified Poisson kernel cutoff coefficient It can be expressed as formula (8): Gravity anomaly nth-order Laplace harmonic function Δg n (θ,λ) is calculated by the bit coefficient model formula (9): Where GM is the Earth's gravitational constant, P nm (cosθ) is the fully normalized associated Legendre function, and is the fully normalized perturbation coefficient; The series form of the calculation of far-zone impact is formula (10): Where N max represents the maximum order of the model, θ represents the pole distance of the calculation point, and λ represents the longitude of the calculation point; Finally, the observation equation is formed: The off-diagonal element A in the i-th row and j-th column of the coefficient matrix ij and the diagonal element A in row i and column i ii They are formula (11) and formula (12) respectively: in, Where, d m (r i ,ψ0,R) represents the intermediate temporary variable, N c Indicates the number of points to be found within the integral radius; The observation equation is formula (14): In the formula, Δg i (r,θ,λ) represents the airborne gravity observation value, Δg j (R,θ′,λ′) represents the extended value to be calculated, represents the far zone influence, and N represents the number of extended values ​​to be calculated; Assuming there are M airborne gravity observations, then according to formula (14), the matrix form is formed, that is, formula (15): b=Ax (15) Where b represents the column vector formed by the M airborne gravity observations on the left side of formula (14) minus the influence of the distant area, A represents the coefficient matrix of M rows and N columns calculated according to formulas (11), (12) and (13), and x represents the column vector composed of N extended values ​​to be calculated; The implementation method of step three is: First, use formula (16) to perform singular value decomposition SVD on the coefficient matrix A: A=UΛV T (16) Where U and V are composed of the left and right eigenvectors respectively, U = [u1 u2…u M ], V=[v1 v2…v N ],Λ=diag(λ1,λ2,…,λ N ) is a diagonal matrix, λ i is a singular value; Then the generalized cross validation GCV method is used to determine the regularization parameter: that is, the selection range and step size of the regularization parameter α are set, and traversal calculation is performed; The regularized solution of the downward extension of the aerial gravity data is expressed as formula (17): In the formula, represents the uncalculated extension value when the regularization parameter is α, u i represents the i-th left eigenvector; Use formula (18) to calculate the intermediate quantity F α Judgment Criteria G α Formula (19): Take the regularization parameter corresponding to the minimum value of formula (19) as the final regularization parameter, and substitute it into formula (17) to calculate the least squares extension solution.

Images