CN103399350B

CN103399350B - A kind of airborne gravity downward continuation method based on integral iteration algorithm

Info

Publication number: CN103399350B
Application number: CN201310322111.XA
Authority: CN
Inventors: 吴美平; 赵磊; 张开东
Original assignee: National University of Defense Technology
Current assignee: National University of Defense Technology
Priority date: 2013-07-29
Filing date: 2013-07-29
Publication date: 2016-02-24
Anticipated expiration: 2033-07-29
Also published as: CN103399350A

Abstract

The invention discloses a kind of airborne gravity downward continuation method based on integral iteration algorithm, when carrying out downward continuation to survey district airborne gravity measurement data, airborne gravity measurement data being carried out upward continuation as the initial value of integral iteration; Upward continuation result upward continuation obtained and airborne gravity measurement data are asked after difference as the feedback quantity of downward continuation result; Feedback quantity is superposed with initial value and obtains new iteration initial value and carry out next iteration.The present invention has simple, easy and simple to handle, the easy popularization of principle, can improve the advantage such as data precision and degree of confidence.

Description

A downward continuation method of aviation gravity based on integral iterative algorithm

技术领域technical field

本发明主要涉及到航空重力测量的技术领域，特指一种基于积分迭代算法的航空重力向下延拓方法。The invention mainly relates to the technical field of airborne gravity measurement, in particular to an airborne gravity downward continuation method based on an integral iterative algorithm.

背景技术Background technique

地球重力场是地球重力场是地球的一种基本物理场，在基础科学、军事和国家安全等领域中的具有重大战略和基础性意义。航空重力测量是以飞机为载体确定区域和局部重力场的方法，其测量速度快、范围广、成本低，是区域范围内获取高精度、中高分辨率重力场信息的最有效手段。The earth's gravitational field is a basic physical field of the earth, which has great strategic and fundamental significance in the fields of basic science, military and national security. Aeronautical gravity measurement is a method of determining the regional and local gravity field with the aircraft as the carrier. Its measurement speed is fast, the range is wide, and the cost is low.

航空重力测量得到的是在航线高度上的重力异常数据，这些测量数据通常在距地球表面几百米甚至上千米的空中。由此获得的数据无法直接运用在地球表面或者大地水准面上，所以需要将航空重力测量数据通过一定的方法向下延拓至地球表面或者大地水准面，进而进行更多的应用。将航空重力测量数据推算到地球表面或者大地水准面的这一过程就是航空重力数据的向下延拓。Aeronautical gravity measurement obtains gravity anomaly data at flight altitude, and these measurement data are usually in the air hundreds of meters or even thousands of meters above the earth's surface. The data thus obtained cannot be directly applied to the earth's surface or the geoid, so it is necessary to extend the airborne gravity measurement data down to the earth's surface or the geoid by a certain method, and then carry out more applications. The process of extrapolating airborne gravity measurement data to the earth's surface or geoid is the downward extension of airborne gravity data.

航空重力测量数据向下延拓是信号放大的非平稳过程，很小的观测噪声会引起待求参数较大的误差，属于非稳定问题的求解。目前，航空重力测量数据的向下延拓依然是根据球面Poisson积分进行求解，比较常见的方法主要有虚拟点质量法、梯度法、最小二乘配置法、直接代表法和正则化方法。The downward continuation of airborne gravity measurement data is a non-stationary process of signal amplification, and small observation noise will cause large errors in the parameters to be obtained, which belongs to the solution of non-stable problems. At present, the downward continuation of airborne gravity measurement data is still solved according to the spherical Poisson integral. The more common methods mainly include virtual point mass method, gradient method, least squares configuration method, direct representation method and regularization method.

虚拟点质量模型方法简单，但受边界效应的影响较大。梯度法不能联合已有的地面数据，而且由于重力异常的垂直梯度在实际应用中难以精确求得。最小二乘配置法能顾及观测量的误差和计算待估量的误差协方差矩阵，但这一方法需要求解超大型协方差矩阵的逆矩阵，实现起来比较困难。直接代表法需要有测区内较精细的地形高数据作为参考值，对缺少地形数据的地区无法实现快速延拓。向下延拓的正则化算法可以有效抑制观测噪声以及重力异常高频分量对延拓精度的影响，但是该方法关键在于求解合适的正则化参数，而且计算比较复杂。The method of the virtual point mass model is simple, but it is greatly affected by the boundary effect. The gradient method cannot be combined with the existing ground data, and it is difficult to accurately obtain the vertical gradient due to the gravity anomaly in practical applications. The least squares configuration method can take into account the error of the observed quantity and calculate the error covariance matrix to be estimated, but this method needs to solve the inverse matrix of a very large covariance matrix, which is difficult to implement. The direct representation method needs finer topographic height data in the survey area as a reference value, and it cannot achieve rapid extension for areas lacking topographic data. The regularization algorithm of downward continuation can effectively suppress the influence of observation noise and gravity anomaly high-frequency components on the continuation accuracy, but the key to this method is to find the appropriate regularization parameters, and the calculation is relatively complicated.

发明内容Contents of the invention

本发明要解决的技术问题就在于：针对现有技术存在的技术问题，本发明提供一种原理简单、操作简便、易推广、可提高数据精度和置信度的基于积分迭代算法的航空重力向下延拓方法。The technical problem to be solved by the present invention is that: aiming at the technical problems existing in the prior art, the present invention provides an air gravity downlink algorithm based on an integral iterative algorithm that is simple in principle, easy to operate, easy to popularize, and can improve data accuracy and confidence. continuation method.

为解决上述技术问题，本发明采用以下技术方案：In order to solve the problems of the technologies described above, the present invention adopts the following technical solutions:

一种基于积分迭代算法的航空重力向下延拓方法，在对测区航空重力测量数据进行向下延拓时，将航空重力测量数据作为积分迭代的初始值进行向上延拓；将向上延拓得到的向上延拓结果与航空重力测量数据求差后作为向下延拓结果的反馈量；将反馈量与初始值叠加得到新的迭代初始值进行下一次迭代。An airborne gravity downward continuation method based on the integral iterative algorithm. When the airborne gravity measurement data of the survey area is extended downward, the airborne gravity measurement data is used as the initial value of the integral iteration to carry out upward continuation; the upward continuation The difference between the obtained upward continuation result and the airborne gravity measurement data is used as the feedback amount of the downward continuation result; the feedback amount and the initial value are superimposed to obtain a new iteration initial value for the next iteration.

作为本发明的进一步改进：所述向下延拓深度不超过观测数据点距的20倍。As a further improvement of the present invention: the downward extension depth does not exceed 20 times of the observation data point pitch.

作为本发明的进一步改进：其具体步骤为：As a further improvement of the present invention: its concrete steps are:

（1）将δg_r(x,y)作为待求面的初始值，即 (1) Take δg _r (x, y) as the initial value of the surface to be found, namely

（2）根据下式，将向上延拓至航空重力测量数据高度，得到 (2) According to the following formula, the Extending upward to the height of the airborne gravity measurement data, we get

${δg δg}_{r r} ((x x,, y the y)) = = {F f}_{- - 11}^{22} {{\overset{~ ~}{K K} ((u u,, v v)) \cdot &Center Dot; {F f}_{22} {{{δg δg}_{R R} ((x x,, y the y))}}}}$

其中，²F与²F^-1分别表示二维傅立叶变换与傅立叶反变换； $\tilde{K} (u, v) = \frac{H}{2 π} {&Integral;}_{- \infty}^{\infty} {&Integral;}_{- \infty}^{\infty} \frac{e^{- i (ux + vy)}}{{(x^{2} + y^{2} + H^{2})}^{3 / 2}} dxdy = e^{- H \sqrt{u^{2} + v^{2}}},$ 积分核 $K (x, y) = \frac{1}{2 π} \frac{H}{{(x^{2} + y^{2} + H^{2})}^{3 / 2}},$ H为距离地面的高度；Wherein, ² F and ² F ^-1 represent two-dimensional Fourier transform and inverse Fourier transform respectively; $\tilde{K} (u, v) = \frac{h}{2 π} {&Integral;}_{- \infty}^{\infty} {&Integral;}_{- \infty}^{\infty} \frac{e^{- i (ux + vy)}}{{(x^{2} + {the y}^{2} + h^{2})}^{3 / 2}} dxdy = e^{- h \sqrt{u^{2} + v^{2}}},$ Integral kernel $K (x, the y) = \frac{1}{2 π} \frac{h}{{(x^{2} + {the y}^{2} + h^{2})}^{3 / 2}},$ H is the height from the ground;

（3）对向上延拓得到的重力异常数据与航空重力测量数据δg_r(x,y)进行比较，并判断是否满足积分迭代过程的结束条件，即：(3) Gravity anomaly data obtained by upward continuation Compare with the airborne gravity measurement data δg _r (x, y), and judge whether the end condition of the integral iteration process is met, namely:

当 ${| | {δg}_{r}^{k} (x, y) - {δg}_{r} (x, y) | |}_{2} < ϵ$ 时，结束迭代过程，跳至步骤（6）；when ${| | {δ g}_{r}^{k} (x, the y) - {δg}_{r} (x, the y) | |}_{2} < ϵ$ , end the iterative process and skip to step (6);

当 ${| | {δg}_{r}^{k} (x, y) - {δg}_{r} (x, y) | |}_{2} &GreaterEqual; ϵ$ 时，继续下一步；when ${| | {δ g}_{r}^{k} (x, the y) - {δ g}_{r} (x, the y) | |}_{2} &Greater Equal; ϵ$ , proceed to the next step;

（4）选择合适的迭代步长s，并按照下列积分迭代公式对待求面的初始值进行更新，得到新的初始值 (4) Select an appropriate iteration step size s, and update the initial value of the surface to be sought according to the following integral iteration formula to obtain a new initial value

${δg δg}_{R R}^{k k + + 11} ((x x,, y the y)) = = {δg δg}_{R R}^{k k} ((x x,, y the y)) + + s the s [[{δg δg}_{r r}^{k k} ((x x,, y the y)) - - {δg δg}_{r r} ((x x,, y the y))]]$

其中k是迭代次数，s是迭代步长；表示经过k次迭代得到的向下延拓结果，是向上延拓至相应高度的结果；Where k is the number of iterations and s is the iteration step size; Indicates the downward continuation result obtained after k iterations, yes The result of extending upwards to the corresponding height;

（5）令k=k+1，重复步骤（2）～（4）。(5) Let k=k+1, repeat steps (2)~(4).

（6）延拓结果为： ${δg}_{R} (x, y) = {δg}_{R}^{k} (x, y) + s [{δg}_{r}^{k} (x, y) - {δg}_{r} (x, y)] .$ (6) The extension result is: ${δg}_{R} (x, the y) = {δg}_{R}^{k} (x, the y) + the s [{δ g}_{r}^{k} (x, the y) - {δ g}_{r} (x, the y)] .$

与现有技术相比，本发明的优点在于：本发明的基于积分迭代算法的航空重力向下延拓方法，原理简单、操作简便、易推广，在计算的数值精度、置信度和效率等方面具有大幅改进。本发明能够克服传统向下延拓方法抗干扰能力弱、容易受到观测噪声和重力测量高频分量影响等缺点，进而提高延拓得到的数据精度，降低误差，提高其置信度。Compared with the prior art, the present invention has the advantages of: the air gravity downward continuation method based on the integral iterative algorithm of the present invention has simple principle, simple and convenient operation, and is easy to popularize. with substantial improvements. The invention can overcome the shortcomings of the traditional downward continuation method, such as weak anti-interference ability and being easily affected by observation noise and high-frequency components of gravity measurement, and further improve the accuracy of data obtained by continuation, reduce errors, and increase its confidence.

附图说明Description of drawings

图1是本发明方法在应用实例中的原理示意图。Fig. 1 is a schematic diagram of the principle of the method of the present invention in an application example.

图2是本发明方法在应用实例中的流程示意图。Fig. 2 is a schematic flow chart of the method of the present invention in an application example.

具体实施方式detailed description

以下将结合说明书附图和具体实施例对本发明做进一步详细说明。The present invention will be further described in detail below in conjunction with the accompanying drawings and specific embodiments.

本发明的一种基于积分迭代算法的航空重力向下延拓方法，为：在对测区航空重力测量数据进行向下延拓时，将航空重力测量数据作为积分迭代的初始值进行向上延拓；将向上延拓得到的向上延拓结果与航空重力测量数据求差后作为向下延拓结果的反馈量；将反馈量与初始值叠加得到新的迭代初始值进行下一次迭代。A method for downward continuation of aerial gravity based on the integral iterative algorithm of the present invention is as follows: when the aerial gravity measurement data of the survey area is extended downward, the aerial gravity measurement data is used as the initial value of the integral iteration to carry out upward continuation ; The difference between the upward continuation result obtained by the upward continuation and the airborne gravity measurement data is used as the feedback amount of the downward continuation result; the new iteration initial value is obtained by superimposing the feedback amount and the initial value for the next iteration.

本发明方法的基本原理如下：The basic principle of the inventive method is as follows:

根据第一边值的Dirichlet问题，如果已知某函数V在球面S外部满足调和方程，那么球外任意一点的函数V_e可以表示为：According to the Dirichlet problem of the first boundary value, if it is known that a certain function V satisfies the harmonic equation outside the sphere S, then the function V _e at any point outside the sphere can be expressed as:

${V V}_{e e} = = \frac{R R (({r r}^{22} - - {R R}^{22}))}{44 π π} \underset{S S}{&Integral; &Integral;} \frac{11}{{l l}^{33}} VdS VdS - - - - - - ((11))$

其中，R是该球体的平均半径，r是球外待求点到球心的距离，l是待求点与球面上某一点的距离。假设r为航空重力测量某一数据点与地球中心的距离，R是地球的平均半径。Among them, R is the average radius of the sphere, r is the distance from the point outside the sphere to the center of the sphere, and l is the distance from the point to be found to a certain point on the sphere. Assume that r is the distance from a certain data point of airborne gravimetry to the center of the earth, and R is the average radius of the earth.

令δg_r与δg_R分别表示地球表面上某一点和航空测量线路上某一点的重力测量值（也即重力异常），则rδg_r满足调和方程，根据式（1）可得：Let δg _r and δg _R denote the gravity measurement value (that is, the gravity anomaly) of a certain point on the earth’s surface and a certain point on the aerial survey line respectively, then rδg _r satisfies the harmonic equation, according to formula (1):

${rδg rδg}_{r r} = = \frac{R R (({r r}^{22} - - {R R}^{22}))}{44 π π} \underset{S S}{&Integral; &Integral;} \frac{11}{{l l}^{33}} {Rδg Rδg}_{R R} dS wxya$

则but

${δg δ g}_{r r} = = \frac{R R (({r r}^{22} - - {R R}^{22}))}{44 πr πr} \underset{S S}{&Integral; &Integral;} \frac{11}{{l l}^{33}} {Rδg Rδg}_{R R} dS wxya - - - - - - ((22))$

将式（2）在球坐标下展开可得Expand equation (2) in spherical coordinates to get

${δg δ g}_{r r} (({θ θ}_{i i},, {λ λ}_{j j})) = = \frac{R R (({r r}^{22} - - {R R}^{22}))}{44 πr πr} {&Integral; &Integral;}_{{θ θ}_{i i}}^{π π} {&Integral; &Integral;}_{{λ λ}_{j j}}^{22 π π} \frac{11}{{l l}^{33}} {δg δ g}_{R R} (({θ θ}_{i i},, {λ λ}_{j j})) sin sin {θ θ}_{i i} dθdλ dθdλ - - - - - - ((33))$

令积分核则式（3）可以看作如下的卷积过程Integral kernel Equation (3) can be regarded as the following convolution process

δg_r(x,y)=K(x,y)*δg_R(x,y)（4）δg _r (x,y)=K(x,y)*δg _R (x,y) (4)

已知上述积分核K(x,y)的傅立叶变换为：It is known that the Fourier transform of the above integral kernel K(x,y) is:

$\overset{~ ~}{K K} ((u u,, v v)) = = \frac{H h}{22 π π} {&Integral; &Integral;}_{- - \infty \infty}^{\infty \infty} {&Integral; &Integral;}_{- - \infty \infty}^{\infty \infty} \frac{{e e}^{- - i i ((ux ux + + vy vy))}}{{(({x x}^{22} + + {y the y}^{22} + + {H h}^{22}))}^{33 / / 22}} dxdy dxdy = = {e e}^{- - H h \sqrt{{u u}^{22} + + {v v}^{22}}}$

则航空重力测量数据的向下延拓可以表示为：Then the downward continuation of airborne gravity measurement data can be expressed as:

${δg δg}_{R R} ((x x,, y the y)) = = {F f}_{- - 11}^{22} {{\frac{{F f}_{22} {{{δg δg}_{r r} ((x x,, y the y))}}}{\overset{~ ~}{K K} ((u u,, v v))}}} - - - - - - ((55))$

其中²F与²F^-1分别表示二维傅立叶变换与傅立叶反变换。Among them, ² F and ² F ^-1 represent two-dimensional Fourier transform and inverse Fourier transform respectively.

如果已知地面异常数据δg_R(x,y)，要得到距离地面高度为H的重力异常，即地面数据向上延拓，可以表示为：If the ground anomaly data δg _R (x, y) is known, the gravity anomaly at a height H from the ground is to be obtained, that is, the ground data is extended upwards, which can be expressed as:

${δg δg}_{r r} ((x x,, y the y)) = = {F f}_{- - 11}^{22} {{\overset{~ ~}{K K} ((u u,, v v)) \cdot &Center Dot; {F f}_{22} {{{δg δg}_{R R} ((x x,, y the y))}}}} - - - - - - ((66))$

由式（5）可知，当进行向下延拓时，积分核K(x,y)将信号进行了放大。这就意味着，如果观测数据中含有噪声，那么观测噪声将被放大，从而导致延拓精度的下降。由式（6）可知，地面重力异常的向上延拓是对信号的平滑过程，可以抑制观测误差和重力数据中的高频分量。It can be seen from formula (5) that when extending downward, the integral kernel K(x, y) amplifies the signal. This means that if the observation data contains noise, the observation noise will be amplified, resulting in a decrease in continuation accuracy. It can be seen from Equation (6) that the upward continuation of the ground gravity anomaly is a smoothing process for the signal, which can suppress observation errors and high-frequency components in the gravity data.

如图1所示，为本发明的基于积分迭代算法的航空重力向下延拓方法的一个具体应用实例，图中“1”表示进行航空重力测量时飞机飞行的航线。进行向下延拓之前，需要将多条航线上的重力测量数据进行格网化。虚线框“2”表示积分迭代过程。其中，“2.1”是迭代初始值的赋值与更新，第一次迭代时，可以将测量数据作为初始值。“2.2”表示由待求面向航线高度延拓过程，是积分迭代方法的必要步骤。As shown in Figure 1, it is a specific application example of the air gravity downward continuation method based on the integral iterative algorithm of the present invention, and "1" in the figure indicates the flight route of the aircraft when the air gravity measurement is carried out. Before the downward continuation, the gravity measurement data on multiple routes need to be gridded. The dotted box "2" represents the integral iterative process. Among them, "2.1" is the assignment and update of the initial value of the iteration. In the first iteration, the measured data can be used as the initial value. "2.2" indicates the continuation process from the altitude of the flight line to be sought, which is a necessary step of the integral iteration method.

基于以上原理，如图2所示，为本发明在具体应用实例中的流程示意图，其中δg_r(x,y)是由航空重力测量得到的重力异常值，是已知量；δg_R(x,y)表示向下延拓的重力异常值，是未知量。本发明的具体流程为：Based on the above principles, as shown in Figure 2, it is a schematic flow chart of the present invention in a specific application example, wherein δg _r (x, y) is a gravity anomaly obtained by airborne gravity measurement, which is a known quantity; δg _R (x ,y) represents the gravity anomaly extending downward, which is an unknown quantity. Concrete flow process of the present invention is:

（1）航空重力测量得到的重力异常表现形式一般是测线或者区域网格。本发明在应用积分迭代方法时，需要将这两种形式的数据按照向下延拓的高度进行格网化。(1) The gravity anomalies obtained by airborne gravity survey are generally in the form of survey lines or regional grids. When the present invention applies the integral iterative method, the two forms of data need to be gridded according to the height of the downward extension.

通常情况下，延拓深度不超过观测数据点距的20倍，由此得到航空重力测量格网化数据δg_r(x,y)。而且，延拓深度值越大，延拓精度越差，置信度越低。Usually, the continuation depth does not exceed 20 times the point distance of the observation data, and thus the airborne gravimetry gridded data δg _r (x, y) is obtained. Moreover, the larger the value of the continuation depth, the worse the continuation accuracy and the lower the confidence.

（2）将δg_r(x,y)作为待求面的初始值，即 (2) Take δg _r (x, y) as the initial value of the surface to be found, namely

（3）根据下式，将向上延拓至航空重力测量数据高度，得到 (3) According to the following formula, the Extending upward to the height of the airborne gravity measurement data, we get

${δg δ g}_{r r} ((x x,, y the y)) = = {F f}_{- - 11}^{22} {{\overset{~ ~}{K K} ((u u,, v v)) \cdot &Center Dot; {F f}_{22} {{{δg δg}_{R R} ((x x,, y the y))}}}}$

（4）对向上延拓得到的重力异常数据与航空重力测量数据δg_r(x,y)进行比较，并判断是否满足积分迭代过程的结束条件，即(4) Gravity anomaly data obtained by upward continuation Compare with the airborne gravity measurement data δg _r (x, y), and judge whether the end condition of the integral iteration process is satisfied, that is,

当 ${| | {δg}_{r}^{k} (x, y) - {δg}_{r} (x, y) | |}_{2} < ϵ$ 时，结束迭代过程，跳至步骤（7）；when ${| | {δ g}_{r}^{k} (x, the y) - {δ g}_{r} (x, the y) | |}_{2} < ϵ$ , end the iterative process and skip to step (7);

当 ${| | {δg}_{r}^{k} (x, y) - {δg}_{r} (x, y) | |}_{2} &GreaterEqual; ϵ$ 时，继续下一步。when ${| | {δ g}_{r}^{k} (x, the y) - {δg}_{r} (x, the y) | |}_{2} &Greater Equal; ϵ$ , continue to the next step.

（5）选择合适的迭代步长s，并按照下列积分迭代公式对待求面的初始值进行更新，得到新的初始值 (5) Select an appropriate iteration step size s, and update the initial value of the surface to be sought according to the following integral iteration formula to obtain a new initial value

${δg δg}_{R R}^{k k + + 11} ((x x,, y the y)) = = {δg δg}_{R R}^{k k} ((x x,, y the y)) + + s the s [[{δg δg}_{r r}^{k k} ((x x,, y the y)) - - {δg δg}_{r r} ((x x,, y the y))]] - - - - - - ((77))$

其中k是迭代次数，s是迭代步长。表示经过k次迭代得到的向下延拓结果，是向上延拓至相应高度的结果。where k is the number of iterations and s is the iteration step size. Indicates the downward continuation result obtained after k iterations, yes The result of extending upwards to the corresponding height.

（6）令k=k+1，重复步骤（3）～（5）。(6) Let k=k+1, repeat steps (3)~(5).

（7）延拓结果为：(7) The extension result is:

${δg δg}_{R R} ((x x,, y the y)) = = {δg δg}_{R R}^{k k} ((x x,, y the y)) + + s the s [[{δg δg}_{r r}^{k k} ((x x,, y the y)) - - {δg δg}_{r r} ((x x,, y the y))]]$

以上仅是本发明的优选实施方式，本发明的保护范围并不仅局限于上述实施例，凡属于本发明思路下的技术方案均属于本发明的保护范围。应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理前提下的若干改进和润饰，应视为本发明的保护范围。The above are only preferred implementations of the present invention, and the protection scope of the present invention is not limited to the above-mentioned embodiments, and all technical solutions under the idea of the present invention belong to the protection scope of the present invention. It should be pointed out that for those skilled in the art, some improvements and modifications without departing from the principle of the present invention should be regarded as the protection scope of the present invention.

Claims

1. The aviation gravity downward continuation method based on the integral iteration algorithm is characterized in that when aviation gravity measurement data of a measurement area is subjected to downward continuation, the aviation gravity measurement data is used as an initial value of integral iteration to be subjected to upward continuation; calculating the difference between the upward continuation result obtained by upward continuation and the aviation gravity measurement data to be used as the feedback quantity of the downward continuation result; overlapping the feedback quantity with the initial value to obtain a new iteration initial value for next iteration; the method comprises the following specific steps:

(1) g is prepared from_r(x, y) as initial values of the plane to be solved, i.e.Wherein, g_r(x, y) is a gravity anomaly value from an airborne gravity measurement, which is a known quantity; g_R(x, y) represents a downwardly extending gravity anomaly value, which is an unknown quantity;

(2) according to the following formulaExtending upwards to the height of the aerial gravity measurement data to obtain

{δg}_{r} (x, y) =^{2} F^{- 1} {\tilde{K} (u, v) \cdot F_{2} {{δg}_{R} (x, y)}}

Wherein,²f and²F^-1respectively representing two-dimensional Fourier transform and inverse Fourier transform;

\tilde{K} (u, v) = \frac{H}{2 π} {&Integral;}_{- \infty}^{\infty} {&Integral;}_{- \infty}^{\infty} \frac{e^{- i (u x + v y)}}{{(x^{2} + y^{2} + H^{2})}^{3 / 2}} d x d y = e^{- H \sqrt{u^{2} + v^{2}}},

integral kernel

K (x, y) = \frac{1}{2 π} \frac{H}{{(x^{2} + y^{2} + H^{2})}^{3 / 2}},

H is the height from the ground;

(3) gravity anomaly data obtained by extending upwardAnd aerial gravity measurement data g_r(x, y) comparing, and judging whether an end condition of the integration iteration process is met, wherein the threshold set for the end condition is as follows:

when in useWhen the iteration process is finished, jumping to the step (6);

when in use

| | {δg}_{r}^{k} (x, y) - {δg}_{r} (x, y) | |_{2} &GreaterEqual; ϵ

If so, continuing the next step;

(4) selecting a proper iteration step length s, and updating the initial value of the to-be-solved surface according to the following integral iteration formula to obtainNew initial value

{δg}_{R}^{k + 1} (x, y) = {δg}_{R}^{k} (x, y) + s [{δg}_{r}^{k} (x, y) - {δg}_{r} (x, y)]

Where k is the number of iterations and s is the iteration step;representing the downward continuation results obtained through k iterations,is thatExtending upwards to a corresponding height;

(5) repeating steps (2) to (4) with k being k + 1;

(6) the continuation result is:

{δg}_{R} (x, y) = {δg}_{R}^{k} (x, y) + s [{δg}_{r}^{k} (x, y) - {δg}_{r} (x, y)] .

2. the aviation gravity downward continuation method based on the integral iterative algorithm as claimed in claim 1, wherein the downward continuation depth is not more than 20 times of the observation data point distance.