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

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 kind of airborne gravity downward continuation method based on integral iteration algorithm

Technical field

The present invention is mainly concerned with the technical field of airborne gravimetry, refers in particular to a kind of airborne gravity downward continuation method based on integral iteration algorithm.

Background technology

The basic physical field of one of earth gravity field to be earth gravity field the be earth, has Major Strategic and basic meaning in the fields such as basic science, military affairs and national security.Airborne gravimetry take aircraft as the method that carrier determines region and Local Gravity Field, and its measuring speed is fast, scope is wide, cost is low, is the most effective means obtaining high precision, middle high-resolution gravity field information in regional extent.

What airborne gravimetry obtained is the gravity anomaly data on enroute altitude, and these measurement data are even going up the aerial of km apart from earth surface hundreds of rice usually.Thus obtained data cannot directly be used on earth surface or geoid surface, so to need airborne gravity measurement data by certain method downward continuation to earth surface or geoid surface, and then more apply.This process of earth surface or geoid surface that airborne gravity measurement data calculated is exactly the downward continuation of airborne gravity data.

Downward Continuation of Airborne Gravity Data is the non-stationary process that signal amplifies, and very little observation noise can cause the error that parameter to be asked is larger, belongs to solving of astable problem.At present, the downward continuation of airborne gravity measurement data is still solve according to sphere Poisson integration, and more common method mainly contains quasi-point mass method, gradient method, least squares collocation, directly method of representatives and regularization method.

Virtual point quality model method is simple, but larger by the impact of boundary effect.Gradient method can not combine existing ground data, and is difficult in actual applications accurately try to achieve due to the VG (vertical gradient) of gravity anomaly.Least squares collocation can be taken the error of observed quantity into account and calculate error co-variance matrix to be estimated, but this method needs the inverse matrix solving ultra-large type covariance matrix, implements more difficult.Direct method of representatives needs surveys landform high data meticulousr in district as with reference to value, cannot realize quick continuation to the area lacking terrain data.The regularization algorithm of downward continuation can effectively suppress observation noise and gravity anomaly high fdrequency component on the impact of continuation precision, but the method key is to solve suitable regularization parameter, and calculates more complicated.

Summary of the invention

The technical problem to be solved in the present invention is just: the technical matters existed for prior art, the invention provides simple, easy and simple to handle, the easy popularization of a kind of principle, can improve the airborne gravity downward continuation method based on integral iteration algorithm of data precision and degree of confidence.

For solving the problems of the technologies described above, the present invention by the following technical solutions:

Based on an airborne gravity downward continuation method for integral iteration algorithm, when carrying out downward continuation to survey district airborne gravity measurement data, airborne gravity measurement data is 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.

As a further improvement on the present invention: the described downward continuation degree of depth is no more than 20 times of observation data point distance.

As a further improvement on the present invention: its concrete steps are:

(1) by δ g _r(x, y) as the initial value in face to be asked, namely

(2) according to following formula, will upward continuation, to airborne gravity measurement data height, obtains

${\mathrm{\δg}}_r(x,y)=F_{-1}^2\{\tilde{K}(u,v)\·F_2\{{\mathrm{\δg}}_R(x,y)\left\}\right\}$

Wherein, ²f with ²f ^-1represent two-dimension fourier transform and inverse fourier transform respectively; $\tilde{K}(u,v)=\frac{H}{2\mathrm{\π}}\underset{-\∞}{\overset{\∞}{\∫}}\underset{-\∞}{\overset{\∞}{\∫}}\frac{e^{-i(\mathrm{ux}+\mathrm{vy})}}{{(x^2+y^2+H^2)}^{3/2}}\mathrm{dxdy}=e^{-H\sqrt{u^2+v^2}},$ Integral kernel $K(x,y)=\frac{1}{2\mathrm{\π}}\frac{H}{{(x^2+y^2+H^2)}^{3/2}},$ H is the height on distance ground;

(3) to the gravity anomaly data that upward continuation obtains with airborne gravity measurement data δ g _r(x, y) compares, and judges whether the termination condition meeting integral iteration process, that is:

When ${\left|\right|{\mathrm{\&delta;g}}_r^k(x,y)-{\mathrm{\&delta;g}}_r(x,y)\left|\right|}_2<\mathrm{\&epsiv;}$ Time, finishing iteration process, skips to step (6);

When ${\left|\right|{\mathrm{\&delta;g}}_r^k(x,y)-{\mathrm{\&delta;g}}_r(x,y)\left|\right|}_2\&GreaterEqual;\mathrm{\&epsiv;}$ Time, continue next step;

(4) select suitable iteration step length s, and treat according to following integral iteration formula and ask the initial value in face to upgrade, obtain new initial value

${\mathrm{\&delta;g}}_R^{k+1}(x,y)={\mathrm{\&delta;g}}_R^k(x,y)+s[{\mathrm{\&delta;g}}_r^k(x,y)-{\mathrm{\&delta;g}}_r(x,y)]$

Wherein k is iterations, and s is iteration step length; represent the downward continuation result obtained through k iteration, be upward continuation is to the result of corresponding height;

(5) make k=k+1, repeat step (2) ~ (4).

(6) continuation result be: ${\mathrm{\&delta;g}}_R(x,y)={\mathrm{\&delta;g}}_R^k(x,y)+s[{\mathrm{\&delta;g}}_r^k(x,y)-{\mathrm{\&delta;g}}_r(x,y)].$

Compared with prior art, the invention has the advantages that: the airborne gravity downward continuation method based on integral iteration algorithm of the present invention, principle is simple, easy and simple to handle, easily promote, have in the numerical precision calculated, degree of confidence and efficiency etc. and significantly improve.The present invention can overcome traditional downward continuation method antijamming capability weak, be easily subject to the shortcoming such as observation noise and the impact of gravimetry high fdrequency component, and then improve the data precision that obtains of continuation, reduce error, improve its degree of confidence.

Accompanying drawing explanation

Fig. 1 is the principle schematic of the inventive method in application example.

Fig. 2 is the schematic flow sheet of the inventive method in application example.

Embodiment

Below with reference to Figure of description and specific embodiment, the present invention is described in further details.

A kind of airborne gravity downward continuation method based on integral iteration algorithm of the present invention, for: when carrying out downward continuation to survey district airborne gravity measurement data, airborne gravity measurement data is 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 ultimate principle of the inventive method is as follows:

According to the Dirichlet problem of First Boundary, if known certain function V meets harmonic equation in Spherical Surface S outside, so function V of the outer any point of ball _ecan be expressed as:

$V_e=\frac{R(r^2-R^2)}{4\mathrm{\&pi;}}\underset{S}{\&Integral;}\frac{1}{l^3}\mathrm{VdS}---\left(1\right)$

Wherein, R is the mean radius of this spheroid, and r is the distance of the outer unknown point of ball to the centre of sphere, and l is certain any distance on unknown point and sphere.Suppose that r is the distance of a certain data point of airborne gravimetry and earth center, R is the mean radius of the earth.

Make δ g _rwith δ g _rto represent on earth surface certain certain any gravimetric measurements (being also gravity anomaly) a bit and on airborne survey circuit respectively, then r δ g _rmeet harmonic equation, can obtain according to formula (1):

${\mathrm{r\&delta;g}}_r=\frac{R(r^2-R^2)}{4\mathrm{\&pi;}}\underset{S}{\&Integral;}\frac{1}{l^3}{\mathrm{R\&delta;g}}_R\mathrm{dS}$

Then

${\mathrm{\&delta;g}}_r=\frac{R(r^2-R^2)}{4\mathrm{\&pi;r}}\underset{S}{\&Integral;}\frac{1}{l^3}{\mathrm{R\&delta;g}}_R\mathrm{dS}---\left(2\right)$

Formula (2) is launched under spherical coordinates can obtain

${\mathrm{\&delta;g}}_r({\mathrm{\&theta;}}_i,{\mathrm{\&lambda;}}_j)=\frac{R(r^2-R^2)}{4\mathrm{\&pi;r}}\underset{{\mathrm{\&theta;}}_i}{\overset{\mathrm{\&pi;}}{\&Integral;}}\underset{{\mathrm{\&lambda;}}_j}{\overset{2\mathrm{\&pi;}}{\&Integral;}}\frac{1}{l^3}{\mathrm{\&delta;g}}_R({\mathrm{\&theta;}}_i,{\mathrm{\&lambda;}}_j)\sin {\mathrm{\&theta;}}_i\mathrm{d\&theta;d\&lambda;}---\left(3\right)$

Make integral kernel then formula (3) can regard following convolution process as

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

The Fourier transform of known above-mentioned integral kernel K (x, y) is:

$\tilde{K}(u,v)=\frac{H}{2\mathrm{\&pi;}}\underset{-\&infin;}{\overset{\&infin;}{\&Integral;}}\underset{-\&infin;}{\overset{\&infin;}{\&Integral;}}\frac{e^{-i(\mathrm{ux}+\mathrm{vy})}}{{(x^2+y^2+H^2)}^{3/2}}\mathrm{dxdy}=e^{-H\sqrt{u^2+v^2}}$

Then the downward continuation of airborne gravity measurement data can be expressed as:

${\mathrm{\&delta;g}}_R(x,y)=F_{-1}^2\left\{\frac{F_2\left\{{\mathrm{\&delta;g}}_r(x,y)\right\}}{\tilde{K}(u,v)}\right\}---\left(5\right)$

Wherein ²f with ²f ^-1represent two-dimension fourier transform and inverse fourier transform respectively.

If known ground abnormal data δ g _r(x, y), obtain distance floor level is the gravity anomaly of H, i.e. ground data upward continuation, can be expressed as:

${\mathrm{\&delta;g}}_r(x,y)=F_{-1}^2\{\tilde{K}(u,v)\&CenterDot;F_2\{{\mathrm{\&delta;g}}_R(x,y)\left\}\right\}---\left(6\right)$

From formula (5), when carrying out downward continuation, signal amplifies by integral kernel K (x, y).This just means, if containing noise in observation data, so observation noise will be exaggerated, thus cause the decline of continuation precision.From formula (6), the upward continuation of ground gravity exception is the smoothing process to signal, can suppress the high fdrequency component in observational error and gravimetric data.

As shown in Figure 1, be an embody rule example of the airborne gravity downward continuation method based on integral iteration algorithm of the present invention, in figure, " 1 " represents the course line of aircraft flight when carrying out airborne gravimetry.Before carrying out downward continuation, need the gravity survey data on many course lines to carry out graticule mesh.Dotted line frame " 2 " represents integral iteration process.Wherein, " 2.1 " are assignment and the renewal of iteration initial value, during first time iteration, and can using measurement data as initial value." 2.2 " represent by waiting to ask towards enroute altitude continuation process, are the steps necessarys of integral iteration method.

Based on above principle, as shown in Figure 2, be the schematic flow sheet of the present invention in embody rule example, wherein δ g _r(x, y) is the GRAVITY ANOMALIES obtained by airborne gravimetry, is known quantity; δ g _r(x, y) represents the GRAVITY ANOMALIES of downward continuation, is unknown quantity.Idiographic flow of the present invention is:

(1) the gravity anomaly form of expression that airborne gravimetry obtains is generally survey line or area grid.The present invention, when applying integral iteration method, needs the data of these two kinds of forms to carry out graticule mesh according to the height of downward continuation.

Under normal circumstances, the continuation degree of depth is no more than 20 times of observation data point distance, obtains airborne gravimetry graticule mesh data δ g thus _r(x, y).And continuation depth value is larger, continuation precision is poorer, and degree of confidence is lower.

(2) by δ g _r(x, y) as the initial value in face to be asked, namely

(3) according to following formula, will upward continuation, to airborne gravity measurement data height, obtains

${\mathrm{\&delta;g}}_r(x,y)=F_{-1}^2\{\tilde{K}(u,v)\&CenterDot;F_2\{{\mathrm{\&delta;g}}_R(x,y)\left\}\right\}$

(4) to the gravity anomaly data that upward continuation obtains with airborne gravity measurement data δ g _r(x, y) compares, and judges whether the termination condition meeting integral iteration process, namely

When ${\left|\right|{\mathrm{\&delta;g}}_r^k(x,y)-{\mathrm{\&delta;g}}_r(x,y)\left|\right|}_2<\mathrm{\&epsiv;}$ Time, finishing iteration process, skips to step (7);

When ${\left|\right|{\mathrm{\&delta;g}}_r^k(x,y)-{\mathrm{\&delta;g}}_r(x,y)\left|\right|}_2\&GreaterEqual;\mathrm{\&epsiv;}$ Time, continue next step.

(5) select suitable iteration step length s, and treat according to following integral iteration formula and ask the initial value in face to upgrade, obtain new initial value

${\mathrm{\&delta;g}}_R^{k+1}(x,y)={\mathrm{\&delta;g}}_R^k(x,y)+s[{\mathrm{\&delta;g}}_r^k(x,y)-{\mathrm{\&delta;g}}_r(x,y)]---\left(7\right)$

Wherein k is iterations, and s is iteration step length. represent the downward continuation result obtained through k iteration, be upward continuation is to the result of corresponding height.

(6) make k=k+1, repeat step (3) ~ (5).

(7) continuation result be:

${\mathrm{\&delta;g}}_R(x,y)={\mathrm{\&delta;g}}_R^k(x,y)+s[{\mathrm{\&delta;g}}_r^k(x,y)-{\mathrm{\&delta;g}}_r(x,y)]$

Below be only the preferred embodiment of the present invention, protection scope of the present invention be not only confined to above-described embodiment, all technical schemes belonged under thinking of the present invention all belong to 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 principles of the present invention, should be considered as protection scope of the present invention.

Claims (2)

1. based on an airborne gravity downward continuation method for integral iteration algorithm, it is characterized in that, 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; Concrete steps are:

(1) by δ g _r(x, y) as the initial value in face to be asked, namely wherein, δ g _r(x, y) is the GRAVITY ANOMALIES obtained by airborne gravimetry, is known quantity; δ g _r(x, y) represents the GRAVITY ANOMALIES of downward continuation, is unknown quantity;

(2) according to following formula, will upward continuation, to airborne gravity measurement data height, obtains

${\mathrm{\&delta;g}}_r(x,y){=}^2{F}^{-1}\{\tilde{K}(u,v)\&CenterDot;F_2\{{\mathrm{\&delta;g}}_R(x,y)\left\}\right\}$

Wherein, ²f with ²f ^-1represent two-dimension fourier transform and inverse fourier transform respectively; $\tilde{K}(u,v)=\frac{H}{2\mathrm{\&pi;}}\underset{-\mathrm{\&infin;}}{\overset{\mathrm{\&infin;}}{\&Integral;}}\underset{-\mathrm{\&infin;}}{\overset{\mathrm{\&infin;}}{\&Integral;}}\frac{e^{-i(ux+vy)}}{{(x^2+y^2+H^2)}^{3/2}}dxdy=e^{-H\sqrt{u^2+v^2}},$ Integral kernel $K(x,y)=\frac{1}{2\mathrm{\&pi;}}\frac{H}{{(x^2+y^2+H^2)}^{3/2}},$ H is the height on distance ground;

(3) to the gravity anomaly data that upward continuation obtains with airborne gravity measurement data δ g _r(x, y) compares, and judges whether the termination condition meeting integral iteration process, and ε is the threshold value of termination condition setting, that is:

When time, finishing iteration process, skips to step (6);

When $\left|\right|{\mathrm{\&delta;g}}_r^k(x,y)-{\mathrm{\&delta;g}}_r(x,y)|{|}_2\&GreaterEqual;\mathrm{\&epsiv;}$ Time, continue next step;

(4) select suitable iteration step length s, and treat according to following integral iteration formula and ask the initial value in face to upgrade, obtain new initial value

${\mathrm{\&delta;g}}_R^{k+1}(x,y)={\mathrm{\&delta;g}}_R^k(x,y)+s\&lsqb;{\mathrm{\&delta;g}}_r^k(x,y)-{\mathrm{\&delta;g}}_r(x,y)\&rsqb;$

Wherein k is iterations, and s is iteration step length; represent the downward continuation result obtained through k iteration, be upward continuation is to the result of corresponding height;

(5) make k=k+1, repeat step (2) ~ (4);

(6) continuation result be: ${\mathrm{\&delta;g}}_R(x,y)={\mathrm{\&delta;g}}_R^k(x,y)+s\&lsqb;{\mathrm{\&delta;g}}_r^k(x,y)-{\mathrm{\&delta;g}}_r(x,y)\&rsqb;.$

2. the airborne gravity downward continuation method based on integral iteration algorithm according to claim 1, is characterized in that, the described downward continuation degree of depth is no more than 20 times of observation data point distance.