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

Aviation gravity downward continuation method based on integral iteration algorithm

Technical Field

The invention mainly relates to the technical field of aviation gravity measurement, in particular to an aviation gravity downward continuation method based on an integral iteration algorithm.

Background

The earth gravitational field is a basic physical field of the earth, and has great strategic and fundamental significance in the fields of basic science, military, national security and the like. The aviation gravity measurement is a method for determining an area and a local gravity field by taking an airplane as a carrier, has high measurement speed, wide range and low cost, and is the most effective means for acquiring high-precision and medium-high resolution gravity field information in an area range.

Airborne gravimetry results in gravity anomaly data at the altitude of the flight line, which is typically airborne several hundred meters or even kilometers from the earth's surface. The data obtained by the method cannot be directly applied to the earth surface or the ground level, so that the aeronautical gravity measurement data needs to be extended downwards to the earth surface or the ground level by a certain method, and further more applications are carried out. The process of extrapolating the airborne gravimetry data to the earth's surface or level is a downward continuation of the airborne gravimetry data.

The downward continuation of the aerial gravity measurement data is a non-stable process of signal amplification, and very small observation noise can cause a large error of a parameter to be solved, belonging to the solving of a non-stable problem. At present, the downward continuation of the aviation gravity measurement data is still solved according to the spherical Poisson integral, and the common methods mainly include a virtual point mass method, a gradient method, a least square configuration method, a direct representation method and a regularization method.

The virtual point quality model method is simple, but is greatly influenced by the boundary effect. The gradient method cannot be combined with the existing ground data, and the vertical gradient due to gravity anomaly is difficult to accurately obtain in practical application. The least square configuration method can take the errors of observed quantities into consideration and calculate the covariance matrix of the errors to be estimated, but the method needs to solve the inverse matrix of the ultra-large covariance matrix and is difficult to realize. The direct representation method needs finer terrain height data in a measuring area as a reference value, and cannot realize rapid continuation for areas lacking terrain data. The downward continuation regularization algorithm can effectively inhibit observation noise and the influence of gravity abnormal high-frequency components on continuation precision, but the method has the key point of solving proper regularization parameters and is complex in calculation.

Disclosure of Invention

The technical problem to be solved by the invention is as follows: aiming at the technical problems in the prior art, the invention provides the aviation gravity downward continuation method based on the integral iteration algorithm, which has the advantages of simple principle, simplicity and convenience in operation and easiness in popularization and can improve the data precision and the confidence coefficient.

In order to solve the technical problems, the invention adopts the following technical scheme:

an aviation gravity downward continuation method based on an 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; and overlapping the feedback quantity with the initial value to obtain a new iteration initial value for next iteration.

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

As a further improvement of the invention: 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.

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

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

Wherein,²f and²F^-1respectively representing two-dimensional Fourier transform and inverse Fourier transform; $\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 from the ground;

(3) gravity anomaly data obtained by extending upwardAnd aerial gravity measurement data g_r(x, y) and judging whether the end condition of the integration iteration process is met, namely:

when in use ${\left|\right|{\mathrm{\&delta;g}}_r^k(x,y)-{\mathrm{\&delta;g}}_r(x,y)\left|\right|}_2<\mathrm{\&epsiv;}$ When the iteration process is finished, jumping to the step (6);

when in use ${\left|\right|{\mathrm{\&delta;g}}_r^k(x,y)-{\mathrm{\&delta;g}}_r(x,y)\left|\right|}_2\&GreaterEqual;\mathrm{\&epsiv;}$ 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 obtain a 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)]$

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) let k = k +1, and repeat steps (2) to (4).

(6) The continuation result is: ${\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 the prior art, the invention has the advantages that: the aviation gravity downward continuation method based on the integral iterative algorithm has the advantages of simple principle, simplicity and convenience in operation and easiness in popularization, and has great improvement in the aspects of numerical accuracy, confidence coefficient, efficiency and the like of calculation. The method can overcome the defects that the traditional downward continuation method is weak in anti-interference capability and is easily influenced by observation noise and gravity measurement high-frequency components, and the like, so that the data precision obtained by continuation is improved, the error is reduced, and the confidence coefficient is improved.

Drawings

FIG. 1 is a schematic diagram of the principle of the method of the present invention in an application example.

FIG. 2 is a schematic flow chart of the method of the present invention in an application example.

Detailed Description

The invention will be described in further detail below with reference to the drawings and specific examples.

The invention discloses an aviation gravity downward continuation method based on an integral iterative algorithm, which comprises the following steps: when the aviation gravity measurement data of the measurement area is extended downwards, the aviation gravity measurement data is used as an initial value of integral iteration to be extended upwards; 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; and overlapping the feedback quantity with the initial value to obtain a new iteration initial value for next iteration.

The basic principle of the method of the invention is as follows:

according to the Dirichlet problem of the first edge value, if a certain function V is known to satisfy the harmonic equation outside the sphere S, the function V at any point outside the sphere S_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 average radius of the sphere, R is the distance from the point to be solved outside the sphere to the center of the sphere, and l is the distance from the point to be solved to a certain point on the sphere. Let R be the distance of a certain data point from the earth's center for aeronautical gravimetry, and R be the average radius of the earth.

Let g_rAnd g_RRepresenting a gravity measurement (i.e. gravity anomaly) at a point on the earth's surface and at a point on the aerial survey line, respectively, rg_rSatisfy blendThe equation, which can be derived from equation (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)$

The formula (2) can be developed under the spherical coordinate

${\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)$

Let integral kernelEquation (3) can be viewed as a convolution process as follows

g_r(x,y)=K(x,y)*g_R(x,y)（4）

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

$\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}}$

the downward continuation of the airborne gravimetry 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 and²F^-1representing a two-dimensional fourier transform and an inverse fourier transform, respectively.

If the ground abnormal data g is known_R(x, y), to obtain a gravity anomaly with a height H from the ground, i.e., a ground data up 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)$

as can be seen from equation (5), when downward continuation is performed, the integrating kernel K (x, y) amplifies the signal. This means that if noise is contained in the observation data, the observation noise is amplified, resulting in a decrease in continuation accuracy. From the equation (6), the upward continuation of the ground gravity anomaly is a smoothing process for the signal, and can suppress the observation error and the high-frequency component in the gravity data.

Referring to fig. 1, a specific application example of the aviation gravity downward continuation method based on the integral iterative algorithm of the present invention is shown, wherein "1" represents a flight path of an airplane during aviation gravity measurement. Before proceeding downward continuation, the gravity measurement data on a plurality of routes needs to be gridded. The dashed box "2" represents the integration iteration process. Wherein, 2.1 is the assignment and update of the iteration initial value, and the measurement data can be used as the initial value in the first iteration. "2.2" indicates the continuation process of the flight path height from the sought flight path, and is a necessary step of the integral iteration method.

Based on the above principle, as shown in FIG. 2, a schematic flow chart of the present invention in a specific application example is shown, 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. The specific process of the invention is as follows:

(1) the abnormal expression form of gravity obtained by aviation gravity measurement is generally a measuring line or an area grid. When the method is applied to the integral iteration method, the data in the two forms needs to be gridded according to the height of downward continuation.

Normally, the continuation depth does not exceed 20 times of the observation data point distance, thereby obtaining the aeronautical gravity measurement gridding data g_r(x, y). Moreover, the larger the extension depth value, the worse the extension accuracy, and the lower the confidence.

(2) G is prepared from_r(x, y) as initial values of the plane to be solved, i.e.

(3) According to the following formulaExtending upwards to the height of the aerial gravity measurement data to obtain

${\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) Gravity anomaly data obtained by extending upwardAnd aerial gravity measurement data g_r(x, y) are compared and it is determined whether the end condition of the integration iteration process is met, i.e.

When in use ${\left|\right|{\mathrm{\&delta;g}}_r^k(x,y)-{\mathrm{\&delta;g}}_r(x,y)\left|\right|}_2<\mathrm{\&epsiv;}$ When the iteration process is finished, jumping to the step (7);

when in use ${\left|\right|{\mathrm{\&delta;g}}_r^k(x,y)-{\mathrm{\&delta;g}}_r(x,y)\left|\right|}_2\&GreaterEqual;\mathrm{\&epsiv;}$ Then, the next step is continued.

(5) 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 obtain a 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)$

Where k is the number of iterations and s is the iteration step.Representing the downward continuation results obtained through k iterations,is thatAnd extending upwards to the corresponding height.

(6) Let k = k +1, and repeat steps (3) to (5).

(7) The continuation result is:

${\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)]$

the above is only a preferred embodiment of the present invention, and the protection scope of the present invention is not limited to the above-mentioned embodiments, and all technical solutions belonging to the idea of the present invention belong to the protection scope of the present invention. It should be noted that modifications and embellishments within the scope of the invention may be made by those skilled in the art without departing from the principle of the invention.

Claims (2)

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

${\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 and²F^-1respectively representing two-dimensional Fourier transform and inverse Fourier transform; $\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 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 $\left|\right|{\mathrm{\&delta;g}}_r^k(x,y)-{\mathrm{\&delta;g}}_r(x,y)|{|}_2\&GreaterEqual;\mathrm{\&epsiv;}$ 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

${\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;$

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: ${\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 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.