CN106646644A - Two-step integral inverse method for determining geoidal surface based on band-limited aerial vector gravity - Google Patents

Abstract

The present invention relates to a two-step integral inverse method for determining a geoidal surface based on band-limited aerial vector gravity. The method includes the following steps that: based on a generalized horizontal boundary value theory, band-limited aerial vector gravity at the altitude of a flight line is adopted to calculate band-limited aeronautical disturbing potential; and an inverse Poisson integral model is adopted to extend downwards the band-limited aeronautical disturbing potential to a sea surface, so that band-limited sea surface disturbing potential can be obtained, and the band-limited sea surface disturbing potential is converted into the geoidal surface through using the Bruns formula. The method of the invention is reasonable in design. With the method adopted, a new solution is provided for calculating the geoidal surface based on the band-limited aerial vector gravity, a calculation function for determining the geoidal surface by using the band-limited aerial vector gravity is realized, and the accuracy of the obtained geoidal surface can satisfy engineering application requirements.

Description

Aviation vector gravitational is limited based on band and determines that two steps of geoid integrate anti-solution

Technical field

The invention belongs to airborne vector gravimetry technical field, especially a kind of to be determined based on band limit aviation vector gravitational Two steps of geoid integrate anti-solution.

Background technology

Airborne vector gravimetry technology is with aircraft as carrier, with inertial navigation system and global navigation satellite positioning system The inertial acceleration information than force information and carrier is united as sensor and obtains respectively, by the joint such as Coordinate Conversion and filtering Calculate the vector gravitational on aviation height.Geoid is an occluding surface for representing the figure of the earth, be defined as with entirely Ball without the optimal closely sealed terrestrial gravitation equipotential surface of the static mean sea level of tide, while be also can reflect earth's internal structure with it is close The physical surface of degree distribution characteristics.Determine that geoid is one of main purpose of airborne vector gravimetry.

At present, determine that the research of geoid is adopted similar to astronmical leveling principle based on aviation vector gravitational, by The horizontal component of vector gravitational is converted into survey line section integration the disturbing potential on course line, by disturbing potential downward continuation to big ground water Behind quasi- face, then undulation of the geoid is calculated using Bruns formula, but, what said method was obtained is level surface relative to the earth, Need to add and be only geoid after benchmark, it is difficult to meet the needs of engineering application.

The content of the invention

It is an object of the invention to overcome the deficiencies in the prior art, there is provided one kind is determined big based on band limit aviation vector gravitational Two steps of ground-level integrate anti-solution, realize that aviation vector gravitational determines the computing function of geoid and the earth for obtaining Level surface precision meets engineering application demand.

The present invention solves existing technical problem and takes technical scheme below to realize：

It is a kind of that the anti-solution of two steps integration that aviation vector gravitational determines geoid is limited based on band, comprise the following steps：

It is step 1, theoretical based on the horizontal boundary values of broad sense, calculate band limit boat using band limit aviation vector gravitational on enroute altitude Empty disturbing potential；

Step 2, using inverse Poisson integral models, by the band limit aviation disturbing potential downward continuation obtained in step 1 to extra large On face, band limit sea disturbance position is obtained, and geoid will be converted into limit sea disturbance position by Bruns formula.

It is as follows that the step 1 calculates the formula with limit aviation disturbing potential：

In formula：T^b(r, θ, λ) is band limit aviation disturbing potential, and r is the earth's core at aviation calculating point to footpath, and θ and λ is aviation meter The colatitude and longitude at point is calculated, b is band limit order, and GM is Gravitational coefficient of the Earth, and R is earth radius, and L is far field truncation funcation Maximum order, l is the reference to gravitational field model order removed, C_n(H,ψ₀) it is band limit aviation Vector operation band limit aviation disturbing potential Far field truncation funcation, H be aviation height, ψ₀It is to integrate radius, T_n(θ, λ) is the Laplace harmonic functions of disturbing gravity position, π It is pi, N is measuring point number in integration radius,WithRespectively band limit aviation vector gravitational measuring point j North-south component and thing component,It is the north-south component with limit aviation vector gravitational measuring point j and thing point Amount calculates the integration kernel function with limit aviation disturbing potential, ψ_jIt is the spherical angle that calculates between point of aviation measuring point and aviation away from Δ σ_jIt is Integration unit area.

Further, the north-south component with limit aviation vector gravitational measuring point j and thing component calculate band limit aviation disturbance The integration kernel function of positionComputing formula be：

θ in formula_jAnd λ_jIt is the colatitude and longitude at aviation measuring point j, P_n(cosψ_j) it is Legendre function.

Further, the far field truncation funcation C that aviation disturbing potential is limited with limit aviation Vector operation band_n(H,ψ₀) calculating Formula is：

R in formula_nm(ψ₀) be Legendre function integral function, be expressed as：

Further, the inverse Poisson integral models of the step 2 are：

In formula：T^b(R, θ ', λ ') is the band limit sea disturbance position of point, (θ ', λ ') be respectively the colatitude of point and Longitude, K^b(R, ψ, r) be Poisson integral models kernel function, its computing formula is：

Inverse Poisson integral models are that computing formula is such as to carrying out process of inverting after Poisson integral model discretizations Under：

T^b(R)=(A^TA)^-1A^TT^b(r)

In formula：T^b(R) it is the expression matrix with limit sea disturbance position, T^b(r) be with limit aviation disturbing potential expression matrix, A It is the expression matrix of Poisson integral models.

Further, the step 2 is based on the formula that Bruns formula calculate geoid：

In formula：N^b(R, θ, λ) is geoid, and γ is normal gravity.

Advantages of the present invention and good effect are：

The present invention is first depending on the horizontal boundary values theory of broad sense, calculates band limit aviation based on band limit aviation vector gravitational and disturbs Position；Then according to inverse Poisson integral models, band limit sea disturbance position is obtained based on band limit disturbing potential downward continuation, and is passed through Bruns formula calculate geoid, and to calculate geoid based on aviation vector gravitational new solution is provided, and lead to Experimental verification is crossed, the precision of the geoid obtained using the method can meet the demand of engineer applied.

Description of the drawings

Fig. 1 a are 2km band limit aviation disturbing gravities north and south horizontal component schematic diagram；

Fig. 1 b are 2km band limit aviation disturbing gravity thing horizontal component schematic diagrames；

Fig. 2 is geoid schematic diagram；

Fig. 3 is the difference value schematic diagram of the geoid that the present invention is obtained and standard value.

Specific embodiment

The embodiment of the present invention is further described below in conjunction with accompanying drawing.

The present invention is a kind of new method that geoid is calculated based on aviation vector gravitational data, main to include in following Hold：It is theoretical according to the horizontal boundary values of broad sense, band limit aviation disturbing potential is calculated based on band limit aviation vector gravitational；According to inverse Poisson Integral model, based on band limit disturbing potential downward continuation band limit sea disturbance position is obtained, and calculates the earth level by Bruns formula Face.

To make the objects, technical solutions and advantages of the present invention become more apparent, develop simultaneously embodiment pair below in conjunction with accompanying drawing Present invention elaborates.

Using the air strips limit aviation vector gravitational north and south of the global order of high-order gravity field model EGM2008 simulation calculations 2160 Horizontal componentWith thing horizontal componentBased on data, the flying height of the aviation vector gravitational is 2000 meters, area The latitude in domain is 36 ° to 41 °, and longitude is 247 ° to 252 °, and resolution ratio is 2.0 '.In order to simulate airborne vector gravimetry data The error of result, band limit aviation vector gravitational horizontal component is introduced respectively observation error white noise error (σ=± 3mGal).As shown in table 1：

Table 1 is with limit aviation horizontal component statistical form/mGal

Fig. 1 a and Fig. 1 b sets forth 2km bands limit aviation disturbing gravity north and south horizontal component and 2km band limit aviation disturbances Gravity thing horizontal component.

In order to check effectiveness of the invention, corresponding geoid N is calculated based on EGM2008 gravity field models^b (R, θ, λ), zoning is 38 ° to 39 °, and longitude is 249 ° to 250 °, and resolution ratio is 2.0 '.

The geoid of table 2 statistics/centimetre

Type	Minimum of a value	Maximum	Mean value	Standard deviation	Middle error	Number
							Geoid	-46.516	50.115	2.878	19.533	19.744	900

Fig. 2 gives geoid schematic diagram.

This aviation vector gravitational data are calculated using the present invention calculate concretely comprising the following steps for geoid：

It is step 1, theoretical based on the horizontal boundary values of broad sense, calculate band limit boat using band limit aviation vector gravitational on enroute altitude Empty disturbing potential, specific formula for calculation is：

In formula：T^b(r, θ, λ) is band limit aviation disturbing potential, and r is the earth's core at aviation calculating point to footpath, and θ and λ is aviation meter The colatitude and longitude at point is calculated, b is band limit order, and GM is Gravitational coefficient of the Earth, and R is earth radius, and L is far field truncation funcation Maximum order, l is the reference to gravitational field model order removed, C_n(H,ψ₀) it is band limit aviation Vector operation band limit aviation disturbing potential Far field truncation funcation, H be aviation height, ψ₀It is to integrate radius, T_n(θ, λ) is the Laplace harmonic functions of disturbing gravity position, π It is pi, N is measuring point number in integration radius,WithRespectively band limit aviation vector gravitational measuring point j North-south component and thing component,It is the north-south component with limit aviation vector gravitational measuring point j and thing point Amount calculates the integration kernel function with limit aviation disturbing potential, ψ_jIt is the spherical angle that calculates between point of aviation measuring point and aviation away from Δ σ_jIt is Integration unit area.

Wherein, the north-south component and thing component with limit aviation vector gravitational measuring point j calculates the product with limit aviation disturbing potential Pyrene functionComputing formula be：

θ in formula_jAnd λ_jIt is the colatitude and longitude at aviation measuring point j, P_n(cosψ_j) it is Legendre function.

Wherein, the far field truncation funcation C of aviation disturbing potential is limited with limit aviation Vector operation band_n(H,ψ₀) computing formula be：

R in formula_nm(ψ₀) be Legendre function integral function, be expressed as：

Step 2, based on inverse Poisson integral models, by the band limit aviation disturbing potential downward continuation obtained in step 1 to extra large On face, band limit sea disturbance position is obtained, and geoid will be converted into limit sea disturbance position by Bruns formula, specifically Computing formula is：

Poisson integral models are：

In formula：T^b(R, θ ', λ ') is the band limit sea disturbance position of point, (θ ', λ ') be respectively the colatitude of point and Longitude, K^b(R, ψ, r) be Poisson integral models kernel function, its computing formula is：

Inverse Poisson integral models are that computing formula is such as to carrying out process of inverting after Poisson integral model discretizations Under：

T^b(R)=(A^TA)^-1A^TT^b(r)

In formula：T^b(R) it is the expression matrix with limit sea disturbance position, T^b(r) be with limit aviation disturbing potential expression matrix, A It is the expression matrix of Poisson integral models.

It is based on the formula of Bruns formula calculating geoid：

In formula：N^b(R, θ, λ) is geoid, and γ is normal gravity.

Step 2 result of calculation is compared with the geoid of standard, statistics such as following table：

Table 3 based on the geoid that two steps integrate anti-solution compare with standard value statistical form/centimetre

	Minimum of a value	Maximum	Mean value	Standard deviation	Middle error	Number
							Fiducial value	-13.856	10.875	-1.085	5.274	5.384	900

Fig. 3 gives the difference value that the geoid with standard value of anti-solution are integrated based on two steps, by table 3 and Fig. 3 As can be seen that on 2000 meters of aviation height, integrating anti-solution using two steps and being based on (σ=± 3mGal) containing white noise error Band limit aviation vector gravitational to calculate the precision of geoid be 5.384 centimetres, fully meet the requirement of engineering application.

It is emphasized that embodiment of the present invention is illustrative, rather than it is determinate, therefore present invention bag The embodiment for being not limited to described in specific embodiment is included, it is every by those skilled in the art's technology according to the present invention scheme The other embodiment for drawing, also belongs to the scope of protection of the invention.

Claims (6)

1. it is a kind of aviation vector gravitational is limited based on band to determine that two steps of geoid integrate anti-solution, it is characterised in that include with Lower step：

It is step 1, theoretical based on the horizontal boundary values of broad sense, calculate band limit aviation using band limit aviation vector gravitational on enroute altitude and disturb Dynamic position；

Step 2, using inverse Poisson integral models, the band obtained in step 1 is limited into aviation disturbing potential downward continuation to sea On, band limit sea disturbance position is obtained, and geoid will be converted into limit sea disturbance position by Bruns formula.

It is 2. according to claim 1 that the anti-solution of two steps integration that aviation vector gravitational determines geoid is limited based on band, It is characterized in that：It is as follows that the step 1 calculates the formula with limit aviation disturbing potential：

$\begin{array}{c}{T}^b(r,\mathrm{\θ},\mathrm{\λ})-\frac{GM}{2R}\underset{n=l}{\overset{L}{\Σ}}n(n+1)C_n(H,{\mathrm{\ψ}}_0)T_n(\mathrm{\θ},\mathrm{\λ})\\ =\frac{GM}{4\mathrm{\π}r}\underset{j=1}{\overset{N}{\Σ}}\[{\mathrm{\δg}}_N^b(r,{\mathrm{\θ}}_j,{\mathrm{\λ}}_j),{\mathrm{\δg}}_E^b(r,{\mathrm{\θ}}_j,{\mathrm{\λ}}_j)\]{\[G_{NT}^b\left({\mathrm{\ψ}}_j\right),G_{ET}^b\left({\mathrm{\ψ}}_j\right)\]}^T{\mathrm{\Δ\σ}}_j,L\≤l+b\end{array}$

In formula：T^b(r, θ, λ) is band limit aviation disturbing potential, and r is the earth's core at aviation calculating point to footpath, and θ and λ is that aviation calculates point The colatitude and longitude at place, b is band limit order, and GM is Gravitational coefficient of the Earth, and R is earth radius, and L is the maximum of far field truncation funcation Order, l is the reference to gravitational field model order removed, C_n(H,ψ₀) it is with the remote of limit aviation Vector operation band limit aviation disturbing potential Area's truncation funcation, H be aviation height, ψ₀It is to integrate radius, T_n(θ, λ) is the Laplace harmonic functions of disturbing gravity position, and π is round Frequency, N is measuring point number in integration radius,WithRespectively with the south of limit aviation vector gravitational measuring point j Northern component and thing component,It is the north-south component with limit aviation vector gravitational measuring point j and thing component meter Calculate the integration kernel function with limit aviation disturbing potential, ψ_jIt is the spherical angle that calculates between point of aviation measuring point and aviation away from Δ σ_jIt is integration Unit area.

It is 3. according to claim 2 that the anti-solution of two steps integration that aviation vector gravitational determines geoid is limited based on band, It is characterized in that：The north-south component with limit aviation vector gravitational measuring point j and thing component are calculated with limit aviation disturbing potential Integration kernel functionComputing formula be：

$\left\{\begin{array}{c}{G}_{NT}^b\left({\mathrm{\ψ}}_j\right)=-\underset{n=l}{\overset{l+b}{\Σ}}\frac{2n+1}{n(n+1)}\frac{\∂P_n\left({\mathrm{cos\ψ}}_j\right)}{\∂{\mathrm{\θ}}_j}\\ G_{ET}^b\left({\mathrm{\ψ}}_j\right)=\frac{1}{{\mathrm{sin\θ}}_j}\underset{n=l}{\overset{l+b}{\Σ}}\frac{2n+1}{n(n+1)}\frac{\∂P_n\left({\mathrm{cos\ψ}}_j\right)}{\∂{\mathrm{\λ}}_j}\end{array}\right.$

θ in formula_jAnd λ_jIt is the colatitude and longitude at aviation measuring point j, P_n(cosψ_j) it is Legendre function.

It is 4. according to claim 2 that the anti-solution of two steps integration that aviation vector gravitational determines geoid is limited based on band, It is characterized in that：The far field truncation funcation C that aviation disturbing potential is limited with limit aviation Vector operation band_n(H,ψ₀) computing formula For：

$C_n(r,{\mathrm{\ψ}}_0)=\underset{m=l}{\overset{l+b}{\Σ}}\frac{2m+1}{m(m+1)}{R}_{nm}\left({\mathrm{\ψ}}_0\right),l\≤n\≤L$

R in formula_nm(ψ₀) be Legendre function integral function, be expressed as：

$R_{nm}\left({\mathrm{\ψ}}_0\right)={\∫}_{x={\mathrm{\ψ}}_0}^{\mathrm{\π}}{P}_n\left(\cos x\right)P_m\left(\cos x\right)\sin xdx.$

It is 5. according to claim 1 that the anti-solution of two steps integration that aviation vector gravitational determines geoid is limited based on band, It is characterized in that：The inverse Poisson integral models of the step 2 are：

$T^b(r,\mathrm{\θ},\mathrm{\λ})=\frac{1}{4\mathrm{\π}}\underset{\mathrm{\Θ}}{\∫\∫}{T}^b(R,{\mathrm{\θ}}^{\′},{\mathrm{\λ}}^{\′})K^b(R,\mathrm{\ψ},r)d({\mathrm{\θ}}^{\′},{\mathrm{\λ}}^{\′}),H>0$

In formula：T^b(R, θ ', λ ') is the band limit sea disturbance position of point, and (θ ', λ ') is respectively the colatitude and longitude of point, K^b(R, ψ, r) be Poisson integral models kernel function, its computing formula is：

$K^b(R,\mathrm{\ψ},r)=\underset{n=l}{\overset{l+b}{\Σ}}(2n+1){\left(\frac{R}{r}\right)}^{N+1}{P}_n\left(cos\mathrm{\ψ}\right)$

Inverse Poisson integral models are that computing formula is as follows to carrying out process of inverting after Poisson integral model discretizations：

T^b(R)=(A^TA)^-1A^TT^b(r)

In formula：T^b(R) it is the expression matrix with limit sea disturbance position, T^bR () is the expression matrix with limit aviation disturbing potential, A is The expression matrix of Poisson integral models.

It is 6. according to claim 1 that the anti-solution of two steps integration that aviation vector gravitational determines geoid is limited based on band, It is characterized in that：The step 2 is based on the formula of Bruns formula calculating geoid：

$N^b(R,\mathrm{\θ},\mathrm{\λ})=\frac{T^b(R,\mathrm{\θ},\mathrm{\λ})}{\mathrm{\γ}}$

In formula：N^b(R, θ, λ) is geoid, and γ is normal gravity.