CN102262248B - Satellite gravity inversion method based on double-satellite spatial three-dimensional interpolation principle - Google Patents

Abstract

The invention relates to a method for precisely measuring an earth gravitational field, in particular to a satellite gravity inversion method based on a double-satellite spatial three-dimensional interpolation principle. An earth disturbing potential (obtained by computing an inter-satellite speed of a GRACE satellite K wave band measuring instrument, a satellite orbit position and a satellite orbit speed of a global positioning system (GPS) receiver, a satellite nonconservative force of an accelerometer and a satellite three-dimensional attitude data of a fixed star sensitive device) which is positioned on a gravity double-satellite orbit with a relatively irregular spatial position is interpolated three-dimensionally onto a reference spherical surface with a relatively regular spatial position and uniform grid division, so that the earth gravitational field is precisely and quickly inversed. By the method, satellite gravity inversion precision is high, physical meaning in the resolving process is clear, the computing speed is high, requirements on performance of a computer are low, and the method is sensitive to high-frequency signals in the gravitational field, so that the double-satellite spatial three-dimensional interpolation method is an effective method for resolving the earth gravitational field with high precision and high spatial resolution.

Description

Satellite gravity inversion method based on double star space three-dimensional interpolation principle

One, technical field

The present invention relates to interleaving techniques fields such as satellite geodesy, geophysics, space science; Particularly relate to a kind of based on will being positioned at earth disturbing potential on the irregular relatively gravity double star track of locus (satellite orbital position and satellite orbit speed, the satellite nonconservative force of accelerometer and the satellite 3 d pose data computation of Star Sensor by speed, GPS receiver between the star of GRACE satellite K wave band measuring instrument obtain) three-dimensional interpolation to the reference sphere that rule and uniform grid are divided relatively of locus, and then accurately and the method for fast inversion earth gravity field.

Two, background technology

As shown in Figure 1,21 century is that human use's Satellite Tracking satellite height/low low technical (SST-HL/LL) promoted " new era of digital earth cognitive ability.Earth gravity field and reflect space distribution, motion and the variation of epigeosphere and inner material is over time determining the fluctuating and the variation of geoid surface simultaneously.Therefore; Confirm the fine structure of earth gravity field and become the demand of being not only geodetic surveying, geophysics, seafari, earthquake prediction, space science, inertial navigation, spationautics, commercial Application etc. at that time, also will important information resources be provided simultaneously for the whole mankind seeks resource, protection environment and prediction disaster.The first, Gravity changer (gravity anomaly) is decided by the unevenness that tectonic structure and mineral resources (like coal, oil, rock gas etc.) distribute.Because have in China's distant view resources reserve 78% oil and 93% rock gas to demand exploitation urgently, therefore accurate gravimetry helps China's resource exploration in the future.Second; Because before and after disasteies such as earthquake (like the 8.5 grades of earthquakes in March, 2005 Indonesia Sumatera, in May, 2008 China's 8.0 grades of earthquakes in Wenchuan, 8.8 grades of earthquakes of in February, 2010 Chile, 9.0 grades of earthquakes of in March, 2011 Fukushima, Japan etc.), volcano, tsunami, wind spout, rubble flow take place; Obvious variation will take place with unusual in earth gravity field; Therefore accurate gravimetry can help us to predict the generation of disaster, and then reduces casualties and economic loss.

The Satellite gravity inverting is meant through analyzing orbital position r and orbital velocity

Interstellar distance ρ ₁₂And speed between star

Nonconservative force f, 3 d pose q, gravity gradient V _Ij(gravitation potential of earth is by the constant C of spherical-harmonic expansion Deng moonscope data and terrestrial gravitation potential coefficient _Lm, S _Lm) relation, set up and find the solution the satellite motion observation equation, and then recover the terrestrial gravitation potential coefficient, final purpose is the earth gravity field of inverting high precision and high spatial resolution.Accurate and rapid solving to moonscope equation (large-scale linear overdetermined equation group) is the inverting gordian technique of gravity field accurately.Up to now, the lot of domestic and foreign scholar is carrying out broad research aspect the earth gravity field inversion method.

1, direct least square method (DLSP): directly inverting and then resolve the terrestrial gravitation potential coefficient through regular square formation.Advantage is that solution procedure is strict, and calculation accuracy is higher; Shortcoming is the large-scale matrix difficulty of directly inverting, and is consuming time more, and needs the high-performance computer support.

In Earth central inertial system, the moonscope equation is set up as follows

G _t×1＝λ _t×1+e _t×1＝A _t×nx _n×1 (1)

Wherein, G _{T * 1}Expression moonscope value vector (mainly comprising speed, nonconservative force, 3 d pose etc. between orbital position, orbital velocity, star), t representes the quantity of satellite orbit observation station; λ _{T * 1}Expression observation true value vector (do not contain measuring error in the observation true value vector, and contain measuring error in the observation vector); x _{N * 1}Represent gravitation potential of earth coefficient vector to be asked, wherein the gravitation potential coefficient is arranged according to exponent number l and is formed (number of times m fixes),

Expression waits to ask the number of gravitation potential coefficient, L _MaxThe expression gravitation potential of earth is by the maximum order of spherical-harmonic expansion; A _{T * n}The design matrix of the capable n row of expression t; e _{T * 1}If expression observed reading error vector is e _{T * 1}Be the separate stochastic error vector of each element, moonscope equation (1) has following statistical nature

$\left\{\begin{array}{c}E\left({\mathrm{\λ}}_{t\×1}\right)=A_{t\×n}{\stackrel{\‾}{x}}_{n\×1}\\ E\left(e_{t\×1}\right)=0\end{array}\right.---\left(2\right)$

Wherein, E representes expected value operator,

Expression x _{N * 1}True value.Because observation equation (1) is large-scale linear overdetermined equation group, therefore there is not routine to separate, have only least square solution.Get with taking advantage of

on (1) formula both sides

$A_{t\×n}^T{G}_{t\×1}{A}_{t\×n}^T{A}_{t\×n}{x}_{n\×1}---\left(3\right)$

Order $y_{n\×1}=A_{t\×n}^T{G}_{t\×1},$ $N_{n\×n}=A_{t\×n}^T{A}_{t\×n},$ Then

y _n×1＝N _n×nx _n×1 (4)

A _{T * n}Be a huge rectangular transition matrix, storage need take a large amount of memory headrooms, and therefore directly storage is difficult realizes.Regular square formation N _{N * n}Though less than A _{T * n}If but directly storage also can take a large amount of memory headroom (resolve 120 terrace ball gravitation potential coefficients and expend the 3Gb internal memory at least).If adopt 30 days moonscope value, the SI is 10s, and then the observation station number is t=259200, if inverting 120 terrace gravity fields then wait to ask the number of gravitation potential coefficient to do

Can know by (4) formula, find the solution N respectively _{N * n}And x _{N * 1}Need tn approximately ²And n ³The double-precision arithmetic of magnitude, the total operand that needs based on DLSP is about tn ²+ n ³Even this operand calculates also very consuming time on present ultra-large type parallel machine, therefore utilize DLSP to find the solution difficulty of large-scale linear overdetermined equation group.

2, pre-service conjugate gradient process of iteration (PCCG): choose the pre-service battle array through optimization and avoided directly inverting of large-scale matrix, the direction of per step iteration is selected to be principle and then to calculate the terrestrial gravitation potential coefficient with the error minimum.Advantage is suitably to have improved computing velocity, has reduced the requirement to computing power, long wave gravity field in being suitable for resolving; Shortcoming is to choose the pre-service battle array to have done approximate processing, but can reduce loss of significance through the number of times that suitably increases loop iteration.

PCCG is one of effective ways of finding the solution at present extensive linear overdetermined equation group, at the memory headroom that guarantees can to significantly improve arithmetic speed under the prerequisite of computational accuracy and only to need hundreds of Mb.In PCCG, N _{N * n}Character most important, have only N _{N * n}Satisfy the condition of positive definite full rank, iteration just possibly restrain in resolving.Juche idea is following: the first, and per step iteration is all treated and is asked parameter correction, till reaching anticipate accuracy; The second, it is principle that the direction of per step iteration is selected with the error minimum; The 3rd, avoid the direct matrix of least square method and invert, find the solution true value through loop iteration.N _{N * n}Be the dense battle array of piece diagonal dominance structure, this character is that iteratively faster is found the solution advantage is provided.Because different design of satellites has different N _{N * n}, so N _{N * n}The singularity otherness that caused the pre-service battle array to select.Because the quality that the pre-service battle array is selected is the gordian technique that whole extensive linear overdetermined equation group is resolved, therefore the optimization process to the observation equation normal matrix is a matter of utmost importance of finding the solution extensive linear overdetermined equation group.

The part of PCCG most critical is pre-service battle array M _{N * n}Choose, optimize two standards of having chosen: the first, Be easy to calculate, can improve the speed of earth gravity field inverting; The second,

With

Approaching as far as possible, can guarantee the precision of earth gravity field inverting.Because N _{N * n}Be piece diagonal dominance square formation, therefore choose N usually _{N * n}Piece diagonal angle part as the pre-service battle array, the M of formation _{N * n}Arrange for pressing number of times m on the principal diagonal, remainder is 0 piece diagonal angle square formation.So choose and not only kept N _{N * n}Principal character, and

Be easy to calculate.In a word, optimize and to choose the pre-service battle array and can reduce the number of times (computing time, directly least square method reduced by 1000 times at least) that PCCG finds the solution loop iteration in the terrestrial gravitation potential coefficient largely.

Get with taking advantage of

on equation (4) both sides

$M_{n\×n}^{-1}{y}_{n\×1}=M_{n\×n}^{-1}{N}_{n\×n}{x}_{n\×1}---\left(5\right)$

Utilize PCCG to find the solution observation equation (5), need not store N _{N * n}, therefore only needing the memory headroom of hundreds of Mb, algorithm thought is following:

The first step, initialization: x ₀=0, r ₀=y-Nx ₀, z ₀=M ^-1r ₀, d ₀=z ₀,

In second step, loop iteration is up to

[1] $q_i={\mathrm{\α}}_i^T\left({\mathrm{\α}}_i{d}_i\right)$ (i＝0，1，2，…)，

[2] ${\mathrm{\α}}_i={\mathrm{\ρ}}_i/d_i^T{q}_i,$

[3]x _i×1＝x _i+α _id _i，

[4]r _i+1＝r _i-α _iq _i，

[5]z _i+1＝M ^-1r _i+1，

[6] ${\mathrm{\ρ}}_{i+1}=r_{i+1}^T{z}_{i+1},$

[7]d _i+1＝z _i+1+ρ _i+1/ρ _id _i，

[8] turn back to [3] loop iteration.

In sum; Though DLSP solving precision in theory is higher; But according to temporary transient difficult popularization of the computing speed of present international super large parallel computer (only NASA jet propulsion laboratory (NASA-JPL), German Potsdam earth science research center international research mechanisms such as (GFZ) adopt); Yet, be expected to progressively become in the future the standard method of widespread use along with the raising day by day of international computer performance; PCCG is one of effective ways that resolve GRACE (Gravity Recovery and Climate Experiment) earth gravity field of the extensive employing of international numerous scientific research institutions at present, but it only is suitable for finding the solution middle low order earth gravity field.

Three, summary of the invention

The objective of the invention is: effectively accelerate computing speed based on double star space three-dimensional method of interpolation, and further improve the precision of earth gravity field inverting.

For achieving the above object, the present invention has adopted following technical scheme:

Accurate and quick Satellite gravity inversion method based on double star space three-dimensional interpolation principle comprises the following step:

1) utilizes GRACE satellite K wave band measuring instrument to obtain speed between star, utilize the GPS receiver to obtain satellite orbital position and satellite orbit speed, utilize accelerometer to obtain the satellite nonconservative force and utilize Star Sensor to obtain satellite 3 d pose observation data, calculate the earth disturbing potential at satellite orbit place fast.In ground was admittedly, (λ) expression formula by spherical-harmonic expansion did earth disturbing potential T for r, θ

$T(r,\mathrm{\θ},\mathrm{\λ})=\frac{\mathrm{GM}}{R_e}\underset{l=2}{\overset{L}{\mathrm{\Σ}}}{\left(\frac{R_e}{r}\right)}^{l+1}\underset{m=0}{\overset{l}{\mathrm{\Σ}}}({\stackrel{\‾}{C}}_{\mathrm{lm}}\cos \mathrm{m\λ}+{\stackrel{\‾}{S}}_{\mathrm{lm}}\sin \mathrm{m\λ}){\stackrel{\‾}{P}}_{\mathrm{lm}}\left(\cos \mathrm{\θ}\right)---\left(6\right)$

Wherein, GM representes that earth quality M and gravitational constant G are long-pending; R _eThe mean radius of the expression earth;

The earth's core radius of expression satellite, x, y, z represent three components of satellite position vector r respectively; θ and λ represent the geocentric colatitude degree and the geocentric longitude of satellite respectively;

Represent normalized Legendre function, l representes exponent number, and m representes number of times; L representes that gravitation potential of earth presses the maximum order of spherical-harmonic expansion;

With Represent normalization terrestrial gravitation potential coefficient to be asked.

2) because the locus of satellite orbit observation station is irregular; Therefore directly expending time at satellite orbit place inverting earth gravity field, (based on CPU is 3.0GHz greatly; In save as the PC of 4Gb; Utilizing the time is to be the moonscope data of 10s with the SI in 30 days, resolves 120 terrace gravity fields 20h at least consuming time).In order to overcome above-mentioned shortcoming consuming time, the present invention is based on novel separating fast and calculated the gravitation potential of earth coefficient with reference to the sphere method.The principle of choosing with reference to sphere number optimum is following: the earth gravity field calculation accuracy is not produced under the standard of materially affect in interpolation error, adopt less reference sphere number as far as possible, and then reduce overall calculation amount and the requirement that reduces computing power.The present invention is based on double star space three-dimensional method of interpolation; Utilizing the time is to be the moonscope data of 10s with the SI in 30 days; And choose 2,3,4,5 and 6 respectively and resolved 120 rank GRACE earth gravity fields with reference to sphere; The result shows: because the error that double star space three-dimensional method of interpolation is introduced is approximately little 10 times than the measuring error of each crucial load of GRACE satellite, can not produce substantial effect to the earth gravity field inversion accuracy, so choose 4 more excellent with reference to sphere.The present invention is that the centre of sphere is selected 4 equidistant and regular reference sphere r with the earth's core ₁, r ₂, r ₃And r ₄, be r with reference to the radius of sphere recently apart from the earth's core ₁=6830km, sphere are spaced apart Δ r=20km, and satellite orbit H=500km is positioned at r ₂And r ₃Between sphere, and on each sphere, carry out uniform grid according to resolution=latitude Δ θ * longitude Δ λ and divide, wherein Δ θ is taken as 180 °/2048, and Δ λ is taken as 360 °/4096.

3) utilize 1 earth disturbing potential observed reading three-dimensional interpolation on the satellite orbit to obtain 4 disturbing potential observed readings with reference to 1 unit cell on the sphere.The resolution of unit cell is designed to Δ θ * Δ λ * Δ r, and wherein Δ θ is taken as 180 °/2048, and Δ λ is taken as 360 °/4096, and Δ r is taken as 20km.The Calculation of Three Dimensional interpolation $T(x,y,z)=\underset{i_x=1}{\overset{4}{\mathrm{\Σ}}}\underset{i_y=1}{\overset{4}{\mathrm{\Σ}}}\underset{i_z=1}{\overset{4}{\mathrm{\Σ}}}{w}_{i_{x}x}{w}_{i_{y}y}{w}_{i_{z}z}T(x_{i_x},y_{i_y},z_{i_z})$ (x, y, z) ∈ [x ₂, x ₃] * [y ₂, y ₃] * [z ₂, z ₃]

$=\underset{i_x=1}{\overset{4}{\mathrm{\Σ}}}{w}_{i_{x}x}\underset{i_y=1}{\overset{4}{\mathrm{\Σ}}}{w}_{i_{y}y}\underset{i_z=1}{\overset{4}{\mathrm{\Σ}}}[w_{1z}T(x_{i_x},y_{i_y},z_1)+w_{2z}T(x_{i_x},y_{i_y},z_2)+w_{3z}T(x_{i_x},y_{i_y},z_3)+w_{4z}T(x_{i_x},y_{i_y},z_4)]---\left(7\right)$

$=\underset{i_x=1}{\overset{4}{\mathrm{\Σ}}}{w}_{i_{x}x}[w_{1y}B(x_{i_x},y_1)+w_{2y}B(x_{i_x},y_2)+w_{3y}B(x_{i_x},y_3)+w_{4y}B(x_{i_x},y_4)]$

$=w_{1x}C\left(x_1\right)+w_{2x}C\left(x_2\right)+w_{3x}C\left(x_3\right)+w_{4x}C\left(x_4\right)$

Wherein, $C\left(x_{i_x}\right)=w_{1y}B(x_{i_x},y_1)+w_{2y}B(x_{i_x},y_2)+w_{3y}B(x_{i_x},y_3)+w_{4y}B(x_{i_x},y_4),$ $B(x_{i_x},y_{i_y})=w_{1z}T(x_{i_x},y_{i_y},z_1)+w_{2z}T(x_{i_x},y_{i_y},z_2)+w_{3z}T(x_{i_x},y_{i_y},z_3)+w_{4z}T(x_{i_x},y_{i_y},z_4),$ T (x; Y; Z) earth disturbing potential of expression on the satellite orbit (the present invention is based on conservation of energy utilizes satellite orbital position and satellite orbit speed, the satellite nonconservative force of accelerometer and the satellite 3 d pose data computation of Star Sensor of speed, GPS receiver between the star of GRACE satellite K wave band measuring instrument to obtain)

Expression is with reference to the rectangular coordinate on the sphere, w _Ab=(a=1,2,3,4; B=x, y, z) expression with the earth disturbing potential three-dimensional interpolation on the satellite orbit to rule with reference to the weight coefficient on the sphere, $w_{1b}=-\frac{1}{2}p{(1-p)}^2,$ $w_{2b}=1-\frac{5}{2}{p}^2+\frac{3}{2}{p}^3,$ $w_{3b}=\frac{1}{2}p(1+4p-3p^2),$ $w_{4b}=\frac{1}{2}{p}^2(p-1),$ $p=\frac{b-b_2}{\mathrm{\Δ}{q}_n},$ P=(b-b ₂)/Δ q _n, b=x, y, z represent the three-dimensional rectangular coordinate on the satellite orbit, b ₂=x ₂, y ₂, z ₂Represent second with reference to the three-dimensional rectangular coordinate on the sphere, Δ q _nSatellite orbit and second are with reference to the distance between sphere in expression.

4) because the calculating of Legendre function

is consuming time than many parts; Therefore fixing each rule is with reference to latitude grid dividing value θ limited on the sphere=2048 parts; By latitude circle calculates

fast and result of calculation storage is called in order to the back, this is the main cause that arithmetic speed is accelerated through the Legendre recursion formula.The rapid solving of Legendre function is to utilize double star space three-dimensional method of interpolation to find the solution one of key issue of large-scale linear overdetermined equation group.In order to reach the purpose of rapid solving Legendre function; The present invention adopts following method: under the prerequisite that guarantees the earth gravity field inversion accuracy; Because the Legendre function is only relevant with colatitude θ; Therefore do not need accurately to calculate one by one the Legendre function of each observation station, only need to calculate the Legendre function of limited colatitude point, and the Legendre functional value on other aspect approaches as Taylor expansion with the functional value on the known point and gets final product.The present invention at first is divided into 2048 equal portions with colatitude; Secondly, calculate the Legendre functional value and the storage of each minizone midpoint.Detailed process is following:

[1] calculates the i colatitude θ of sampled point constantly _i, and the minizone at search place is [θ _j, θ _J+1];

[2 supposition θ _J+1/2Be interval [θ _j, θ _J+1] mid point, P then _Lm(cos θ _i) the available P of calculating _Lm(cos θ _J+1/2) Taylor expansion approach,

$P_{\mathrm{lm}}\left(\cos {\mathrm{\θ}}_i\right)=P_{\mathrm{lm}}\left(\cos {\mathrm{\θ}}_{j+1/2}\right)+\underset{k=1}{\overset{L}{\mathrm{\Σ}}}\left[\frac{1}{k!}{P}_{\mathrm{lm}}^{\left(k\right)}\left(\cos {\mathrm{\θ}}_{j+1/2}\right){({\mathrm{\θ}}_i-{\mathrm{\θ}}_{j+1/2})}^k\right]+O\left(\mathrm{\δ}{\mathrm{\θ}}^{L+1}\right)---\left(8\right)$

Wherein, O (δ θ ^L+1) expression P _Lm(cos θ _i) high-order that launches by Taylor's formula in a small amount, δ θ ^L+1=(θ _i-θ _J+1/2) ^L+1

The calculating of normalization association Legendre function and first order derivative is the key factor of the large-scale linear overdetermined equation set of calculated speed of influence.The present invention will provide a kind of recursive algorithm of improved regular association Legendre function, and concrete formula is following:

Suppose u=sin θ; W=cos θ; regular association Legendre function and can be expressed as following form then to the first derivative of θ

$\left\{\begin{array}{cc}\frac{{\stackrel{\‾}{P}}_{\mathrm{lm}}\left(\cos \mathrm{\θ}\right)}{u^m}=\frac{1}{\sqrt{k}}[g_{\mathrm{lm}}w\frac{{\stackrel{\‾}{P}}_{l,m+1}\left(\cos \mathrm{\θ}\right)}{u^{m+1}}-h_{\mathrm{lm}}{u}^2\frac{{\stackrel{\‾}{P}}_{l,m+2}\left(\cos \mathrm{\θ}\right)}{u^{m+2}}]& (l>m)\\ \frac{{\stackrel{\‾}{P}}_{\mathrm{mm}}\left(\cos \mathrm{\θ}\right)}{u^m}=\sqrt{3}\underset{i=2}{\overset{m}{\mathrm{\Π}}}\sqrt{\frac{2i+1}{2i}}& (m\≥1)\\ \frac{{\stackrel{\‾}{P}}_{\mathrm{lm}}^{\′}\left(\cos \mathrm{\θ}\right)}{u^{m-1}}=\mathrm{mw}\frac{{\stackrel{\‾}{P}}_{\mathrm{lm}}\left(\cos \mathrm{\θ}\right)}{u^m}-e_{\mathrm{lm}}{u}^2\frac{{\stackrel{\‾}{P}}_{l,m+1}\left(\cos \mathrm{\θ}\right)}{u^{m+1}}& (l\≥m)\end{array}\right.---\left(9\right)$

Wherein, $k=\left\{\begin{array}{cc}2& (m=0)\\ 1& (m>0)\end{array}\right.,$ $g_{\mathrm{lm}}=\frac{2(m+1)}{\sqrt{(l-m)(l+m+1)}},$ $h_{\mathrm{lm}}=\sqrt{\frac{(l+m+2)(l-m-1)}{(l-m)(l+m+1)}},$

$e_{\mathrm{lm}}=\sqrt{\frac{(l+m+1)(l-m)}{k}}.$

Order ${\tilde{P}}_{\mathrm{Lm}}\left(\mathrm{Cos}\mathrm{\θ}\right)=\frac{{\stackrel{\‾}{P}}_{\mathrm{Lm}}\left(\mathrm{Cos}\mathrm{\θ}\right)}{u^m}$ With ${\tilde{P}}_{\mathrm{Lm}}^{\′}\left(\mathrm{Cos}\mathrm{\θ}\right)=\frac{{\stackrel{\‾}{P}}_{\mathrm{Lm}}^{\′}\left(\mathrm{Cos}\mathrm{\θ}\right)}{u^{m-1}},$ Therefore (9) formula is rewritten as,

$\left\{\begin{array}{cc}{\tilde{P}}_{\mathrm{lm}}\left(\cos \mathrm{\θ}\right)=\frac{1}{\sqrt{k}}[g_{\mathrm{lm}}w{\tilde{P}}_{l,m+1}\left(\cos \mathrm{\θ}\right)-h_{\mathrm{lm}}{u}^2{\tilde{P}}_{l,m+2}\left(\cos \mathrm{\θ}\right)]& (l>m)\\ {\tilde{P}}_{\mathrm{mm}}\left(\cos \mathrm{\θ}\right)=\sqrt{3}\underset{i=2}{\overset{m}{\mathrm{\Π}}}\sqrt{\frac{2i+1}{2i}}& (m\≥1)\\ {\tilde{P}}_{\mathrm{lm}}^{\′}\left(\cos \mathrm{\θ}\right)=\mathrm{mw}{\tilde{P}}_{\mathrm{lm}}\left(\cos \mathrm{\θ}\right)-e_{\mathrm{lm}}{u}^2{\tilde{P}}_{l,m+1}\left(\cos \mathrm{\θ}\right)& (l\≥m)\end{array}\right.---\left(10\right)$

5) utilize double star space three-dimensional method of interpolation to resolve large-scale overdetermined equation group (number of equation is much larger than the number of unknown number) fast, and then simulate the terrestrial gravitation potential coefficient.In Earth central inertial system, earth disturbing potential moonscope equation (6) is represented as follows with matrix form

D _t×1＝A _t×n·x _n×t (11)

Wherein, D _{T * 1}Expression is with the earth disturbing potential observed reading T (x at satellite orbit place; Y; Z) (the present invention is based on conservation of energy utilizes the satellite orbital position of speed, GPS receiver between the star of GRACE satellite K wave band measuring instrument and satellite orbit speed, the satellite nonconservative force of accelerometer and the satellite 3 d pose data computation of Star Sensor to obtain) space three-dimensional is interpolated into reference to the column matrix behind the sphere, and t representes the quantity of observation station; x _{N * 1}Expression waits to ask the gravitation potential of earth coefficient vector Wherein the gravitation potential coefficient is arranged according to exponent number l and is formed (number of times m fixes),

Expression waits to ask the number of gravitation potential coefficient, L _MaxExpression gravitation position is by the maximum order of spherical-harmonic expansion; A _{T * n}The design matrix of the capable n row of expression t (on earth disturbing potential observation equation (6) the right, removes the terrestrial gravitation potential coefficient

Other outer part).On (11) formula both sides with take advantage of

can be to be asked the terrestrial gravitation potential coefficient

$x_{n\×1}={\left(A_{t\×n}^T{A}_{t\×n}\right)}^{-1}\·A_{t\×n}^T{D}_{t\×1}---\left(12\right)$

Order ${\tilde{A}}_{t\×n}={\left(A_{t\×n}^T{A}_{t\×n}\right)}^{-1}\·A_{t\×n}^T,$ Therefore (12) formula deformable does

$x_{n\×1}={\tilde{A}}_{t\×n}{D}_{t\×1}---\left(13\right)$

When finding the solution large-scale system of linear equations (13), resolving fast of

is core link.Because directly finding the solution and storing

is all difficult; Therefore in order effectively to resolve GRACE-type Satellite Tracking satellite mode observation equation and then fast inversion high precision and high spatial resolution earth gravity field, the present invention proposes the aforementioned new thought of utilizing double star space three-dimensional method of interpolation accurately and fast to resolve .

The present invention is based on that double star space three-dimensional interpolation principle helps accurately and the characteristics of fast inversion earth gravity field design, and advantage is: 1) the earth gravity field inversion accuracy is higher; 2) Satellite gravity refutation process physical meaning is clear and definite; 3) under the prerequisite that guarantees the gravity field inversion accuracy, improved computing velocity largely; 4) to computing power require relatively low; 5) responsive to earth gravity field medium short wave signal.

Four, description of drawings

Fig. 1 representes the measuring principle of Satellite Tracking satellite height/low observation mode, not only with high rail gps satellite low rail gravity double star real-time accurate is followed the tracks of orbit determination (SST-HL), measures the mutual motion (SST-LL) between the low rail double star with differential principle simultaneously.

Fig. 2 representes the GRACE Gravity Satellite of NASA (NASA) and German space agency (DLR) cooperation emission.

Fig. 3 representes by 1 earth disturbing potential observed reading three-dimensional interpolation to 4 on the satellite orbit with reference to 1 unit cell on the sphere.

Fig. 4 representes the Satellite gravity inversion method schematic diagram based on three-dimensional interpolation principle between new type double starry sky.

Fig. 5 representes the precision based on three-dimensional interpolation method inverting earth gravity field between new type double starry sky, and wherein horizontal ordinate representes that gravitation potential of earth presses the exponent number that spheric function is launched, and ordinate is represented the cumulative errors (unit: m) of geoid surface.

Five, embodiment

Below in conjunction with accompanying drawing, specific embodiments of the invention is further described.

Application based on the accurate and quick Satellite gravity inversion method of three-dimensional interpolation principle between new type double starry sky:

1) utilizes GRACE satellite (as shown in Figure 2) K wave band measuring instrument to obtain speed between star, utilize the GPS receiver to obtain satellite orbital position and satellite orbit speed, utilize accelerometer to obtain the satellite nonconservative force and utilize Star Sensor to obtain satellite 3 d pose observation data, calculate the earth disturbing potential at satellite orbit place fast.In ground was admittedly, (λ) expression formula by spherical-harmonic expansion did earth disturbing potential T for r, θ

$T(r,\mathrm{\θ},\mathrm{\λ})=\frac{\mathrm{GM}}{R_e}\underset{l=2}{\overset{L}{\mathrm{\Σ}}}{\left(\frac{R_e}{r}\right)}^{l+1}\underset{m=0}{\overset{l}{\mathrm{\Σ}}}({\stackrel{\‾}{C}}_{\mathrm{lm}}\cos \mathrm{m\λ}+{\stackrel{\‾}{S}}_{\mathrm{lm}}\sin \mathrm{m\λ}){\stackrel{\‾}{P}}_{\mathrm{lm}}\left(\cos \mathrm{\θ}\right)$

Wherein, GM representes that earth quality M and gravitational constant G are long-pending; R _eThe mean radius of the expression earth;

The earth's core radius of expression satellite, x, y, z represent three components of satellite position vector r respectively; θ and λ represent the geocentric colatitude degree and the geocentric longitude of satellite respectively; Represent normalized Legendre function, l representes exponent number, and m representes number of times; L representes that gravitation potential of earth presses the maximum order of spherical-harmonic expansion;

With Represent normalization terrestrial gravitation potential coefficient to be asked.

2) because the locus of satellite orbit observation station is irregular; Therefore directly expending time at satellite orbit place inverting earth gravity field, (based on CPU is 3.0GHz greatly; In save as the PC of 4Gb; Utilizing the time is to be the moonscope data of 10s with the SI in 30 days, resolves 120 terrace gravity fields 20h at least consuming time).In order to overcome above-mentioned shortcoming consuming time, the present invention is based on novel separating fast and calculated the gravitation potential of earth coefficient with reference to the sphere method.The principle of choosing with reference to sphere number optimum is following: the earth gravity field calculation accuracy is not produced under the standard of materially affect in interpolation error, adopt less reference sphere number as far as possible, and then reduce overall calculation amount and the requirement that reduces computing power.The present invention is based on double star space three-dimensional method of interpolation; Utilizing the time is to be the moonscope data of 10s with the SI in 30 days; And choose 2,3,4,5 and 6 respectively and resolved 120 rank GRACE earth gravity fields with reference to sphere; The result shows: because the error that double star space three-dimensional method of interpolation is introduced is approximately little 10 times than the measuring error of each crucial load of GRACE satellite, can not produce substantial effect to the earth gravity field inversion accuracy, so choose 4 more excellent with reference to sphere.The present invention is that the centre of sphere is selected 4 equidistant and regular reference sphere r with the earth's core ₁, r ₂, r ₃And r ₄, be r with reference to the radius of sphere recently apart from the earth's core ₁=6830km, sphere are spaced apart Δ r=20km, and satellite orbit H=500km is positioned at r ₂And r ₃Between sphere, and on each sphere, carry out uniform grid according to resolution=latitude Δ θ * longitude Δ λ and divide, wherein Δ θ is taken as 180 °/2048, and Δ λ is taken as 360 °/4096.

3) utilize 1 earth disturbing potential observed reading three-dimensional interpolation on the satellite orbit to obtain 4 disturbing potential observed readings with reference to 1 unit cell on the sphere.As shown in Figure 3, the resolution of unit cell is designed to Δ θ * Δ λ * Δ r, and wherein Δ θ is taken as 180 °/2048, and Δ λ is taken as 360 °/4096, and Δ λ is taken as 20km.Double star space three-dimensional interpolation formula is represented as follows:

$T(x,y,z)=\underset{i_x=1}{\overset{4}{\mathrm{\Σ}}}\underset{i_y=1}{\overset{4}{\mathrm{\Σ}}}\underset{i_z=1}{\overset{4}{\mathrm{\Σ}}}{w}_{i_{x}x}{w}_{i_{y}y}{w}_{i_{z}z}T(x_{i_x},y_{i_y},z_{i_z})$ (x，y，z)∈[x ₂，x ₃]×[y ₂，y ₃]×[z ₂，z ₃]

$=\underset{i_x=1}{\overset{4}{\mathrm{\Σ}}}{w}_{i_{x}x}\underset{i_y=1}{\overset{4}{\mathrm{\Σ}}}{w}_{i_{y}y}\underset{i_z=1}{\overset{4}{\mathrm{\Σ}}}[w_{1z}T(x_{i_x},y_{i_y},z_1)+w_{2z}T(x_{i_x},y_{i_y},z_2)+w_{3z}T(x_{i_x},y_{i_y},z_3)+w_{4z}T(x_{i_x},y_{i_y},z_4)]$

$=\underset{i_x=1}{\overset{4}{\mathrm{\Σ}}}{w}_{i_{x}x}[w_{1y}B(x_{i_x},y_1)+w_{2y}B(x_{i_x},y_2)+w_{3y}B(x_{i_x},y_3)+w_{4y}B(x_{i_x},y_4)]$

$=w_{1x}C\left(x_1\right)+w_{2x}C\left(x_2\right)+w_{3x}C\left(x_3\right)+w_{4x}C\left(x_4\right)$

Wherein, $C\left(x_{i_x}\right)=w_{1y}B(x_{i_x},y_1)+w_{2y}B(x_{i_x},y_2)+w_{3y}B(x_{i_x},y_3)+w_{4y}B(x_{i_x},y_4),$ $B(x_{i_x},y_{i_y})=w_{1z}T(x_{i_x},y_{i_y},z_1)+w_{2z}T(x_{i_x},y_{i_y},z_2)+w_{3z}T(x_{i_x},y_{i_y},z_3)+w_{4z}T(x_{i_x},y_{i_y},z_4),$ T (x; Y; Z) earth disturbing potential of expression on the satellite orbit (the present invention is based on conservation of energy utilizes satellite orbital position and satellite orbit speed, the satellite nonconservative force of accelerometer and the satellite 3 d pose data computation of Star Sensor of speed, GPS receiver between the star of GRACE satellite K wave band measuring instrument to obtain)

Expression is with reference to the rectangular coordinate on the sphere, w _Ab(a=1,2,3,4; B=x, y, z) expression with the earth disturbing potential three-dimensional interpolation on the satellite orbit to rule with reference to the weight coefficient on the sphere, $w_{1b}=-\frac{1}{2}p{(1-p)}^2,$ $w_{2b}=1-\frac{5}{2}{p}^2+\frac{3}{2}{p}^3,$ $w_{3b}=\frac{1}{2}p(1+4p-3p^2),$ $w_{4b}=\frac{1}{2}{p}^2(p-1),$ $p=\frac{b-b_2}{\mathrm{\Δ}{q}_n},$ P=(b-b ₂)/Δ q _n, b=x, y, z represent the three-dimensional rectangular coordinate on the satellite orbit, b ₂=x ₂, y ₂, z ₂Represent second with reference to the three-dimensional rectangular coordinate on the sphere, Δ q _nSatellite orbit and second are with reference to the distance between sphere in expression.The weight coefficient of double star space three-dimensional interpolation formula is represented that with row matrix

then double star space three-dimensional interpolation formula matrix form is represented as follows:

$T(x,y,z)=W(w_{i_{x}x},w_{i_{y}y},w_{i_{z}z})\·T(x_{i_x},y_{i_y},z_{i_z})$

Therefore can pass through after the double star space three-dimensional interpolation 4 with reference to the disturbing potential on the sphere

$T(w_{i_x},w_{i_y},w_{i_z})={\left(W^{T}W\right)}^{-1}{W}^{T}T(x,y,z)$

4) because the calculating of Legendre function

is consuming time than many parts; Therefore fixing each rule is with reference to latitude grid dividing value θ limited on the sphere=2048 parts; By latitude circle calculates

fast and result of calculation storage is called in order to the back, this is the main cause that arithmetic speed is accelerated through the Legendre recursion formula.The rapid solving of Legendre function is to utilize double star space three-dimensional method of interpolation to find the solution one of key issue of large-scale linear overdetermined equation group.In order to reach the purpose of rapid solving Legendre function; The present invention adopts following method: under the prerequisite that guarantees the earth gravity field inversion accuracy; Because the Legendre function is only relevant with colatitude θ; Therefore do not need accurately to calculate one by one the Legendre function of each observation station, only need to calculate the Legendre function of limited colatitude point, and the Legendre functional value on other aspect approaches as Taylor expansion with the functional value on the known point and gets final product.The present invention at first is divided into 2048 equal portions with colatitude; Secondly, calculate the Legendre functional value and the storage of each minizone midpoint.Detailed process is following:

[1] calculates the i colatitude θ of sampled point constantly _i, and the minizone at search place is [θ _j, θ _I+1];

[2] supposition θ _I+1/2Be interval [θ _i, θ _I+1] mid point, P then _Lm(cos θ _i) the available P of calculating _Lm(cos θ _I+1/2) Taylor expansion approach,

$P_{\mathrm{lm}}\left(\cos {\mathrm{\θ}}_i\right)=P_{\mathrm{lm}}\left(\cos {\mathrm{\θ}}_{j+1/2}\right)+\underset{k=1}{\overset{L}{\mathrm{\Σ}}}\left[\frac{1}{k!}{P}_{\mathrm{lm}}^{\left(k\right)}\left(\cos {\mathrm{\θ}}_{j+1/2}\right){({\mathrm{\θ}}_i-{\mathrm{\θ}}_{j+1/2})}^k\right]+O\left(\mathrm{\δ}{\mathrm{\θ}}^{L+1}\right)$

Wherein, O (δ θ ^L+1) expression P _Lm(cos θ _i) high-order that launches by Taylor's formula in a small amount, δ θ ^L+1=(θ _i-θ _J+1/2) ^L+1

The calculating of normalization association Legendre function and first order derivative is the key factor of the large-scale linear overdetermined equation set of calculated speed of influence.The present invention will provide a kind of recursive algorithm of improved regular association Legendre function, and concrete formula is following:

Suppose u=sin θ; W=cos θ;

regular association Legendre function and can be expressed as following form then to the first derivative of θ

$\left\{\begin{array}{cc}\frac{{\stackrel{\‾}{P}}_{\mathrm{lm}}\left(\cos \mathrm{\θ}\right)}{u^m}=\frac{1}{\sqrt{k}}[g_{\mathrm{lm}}w\frac{{\stackrel{\‾}{P}}_{l,m+1}\left(\cos \mathrm{\θ}\right)}{u^{m+1}}-h_{\mathrm{lm}}{u}^2\frac{{\stackrel{\‾}{P}}_{l,m+2}\left(\cos \mathrm{\θ}\right)}{u^{m+2}}]& (l>m)\\ \frac{{\stackrel{\‾}{P}}_{\mathrm{mm}}\left(\cos \mathrm{\θ}\right)}{u^m}=\sqrt{3}\underset{i=2}{\overset{m}{\mathrm{\Π}}}\sqrt{\frac{2i+1}{2i}}& (m\≥1)\\ \frac{{\stackrel{\‾}{P}}_{\mathrm{lm}}^{\′}\left(\cos \mathrm{\θ}\right)}{u^{m-1}}=\mathrm{mw}\frac{{\stackrel{\‾}{P}}_{\mathrm{lm}}\left(\cos \mathrm{\θ}\right)}{u^m}-e_{\mathrm{lm}}{u}^2\frac{{\stackrel{\‾}{P}}_{l,m+1}\left(\cos \mathrm{\θ}\right)}{u^{m+1}}& (l\≥m)\end{array}\right.$

Wherein, $k=\left\{\begin{array}{cc}2& (m=0)\\ 1& (m>0)\end{array}\right.,$ $g_{\mathrm{lm}}=\frac{2(m+1)}{\sqrt{(l-m)(l+m+1)}},$ $h_{\mathrm{lm}}=\sqrt{\frac{(l+m+2)(l-m-1)}{(l-m)(l+m+1)}},$

$e_{\mathrm{lm}}=\sqrt{\frac{(l+m+1)(l-m)}{k}}.$

Order ${\tilde{P}}_{\mathrm{Lm}}\left(\mathrm{Cos}\mathrm{\θ}\right)=\frac{{\stackrel{\‾}{P}}_{\mathrm{Lm}}\left(\mathrm{Cos}\mathrm{\θ}\right)}{u^m}$ With ${\tilde{P}}_{\mathrm{Lm}}^{\′}\left(\mathrm{Cos}\mathrm{\θ}\right)=\frac{{\stackrel{\‾}{P}}_{\mathrm{Lm}}^{\′}\left(\mathrm{Cos}\mathrm{\θ}\right)}{u^{m-1}},$ Therefore following formula can be rewritten as,

$\left\{\begin{array}{cc}{\tilde{P}}_{\mathrm{lm}}\left(\cos \mathrm{\θ}\right)=\frac{1}{\sqrt{k}}[g_{\mathrm{lm}}w{\tilde{P}}_{l,m+1}\left(\cos \mathrm{\θ}\right)-h_{\mathrm{lm}}{u}^2{\tilde{P}}_{l,m+2}\left(\cos \mathrm{\θ}\right)]& (l>m)\\ {\tilde{P}}_{\mathrm{mm}}\left(\cos \mathrm{\θ}\right)=\sqrt{3}\underset{i=2}{\overset{m}{\mathrm{\Π}}}\sqrt{\frac{2i+1}{2i}}& (m\≥1)\\ {\tilde{P}}_{\mathrm{lm}}^{\′}\left(\cos \mathrm{\θ}\right)=\mathrm{mw}{\tilde{P}}_{\mathrm{lm}}\left(\cos \mathrm{\θ}\right)-e_{\mathrm{lm}}{u}^2{\tilde{P}}_{l,m+1}\left(\cos \mathrm{\θ}\right)& (l\≥m)\end{array}\right..$

5) utilize double star space three-dimensional method of interpolation to resolve large-scale overdetermined equation group (number of equation is much larger than the number of unknown number) fast, and then simulate the terrestrial gravitation potential coefficient.In Earth central inertial system, earth disturbing potential moonscope equation is represented as follows with matrix form

D _t×1＝A _t×n·x _n×1

Wherein, D _{T * 1}Expression is based on the earth disturbing potential observed reading T (x of double star space three-dimensional method of interpolation with the satellite orbit place; Y; Z) (the present invention is based on conservation of energy utilizes the satellite orbital position of speed, GPS receiver between the star of GRACE satellite K wave band measuring instrument and satellite orbit speed, the satellite nonconservative force of accelerometer and the satellite 3 d pose data computation of Star Sensor to obtain) space three-dimensional is interpolated into reference to the column matrix behind the sphere, and t representes the quantity of observation station; x _{N * 1}Expression waits to ask the gravitation potential of earth coefficient vector

Wherein the gravitation potential coefficient is arranged according to exponent number l and is formed (number of times m fixes),

Expression waits to ask the number of gravitation potential coefficient, L _MaxExpression gravitation position is by the maximum order of spherical-harmonic expansion; a _{T * n}The design matrix of the capable n row of expression t is (on earth disturbing potential moonscope equation the right, except the terrestrial gravitation potential coefficient Other outer part).On the following formula both sides with take advantage of

can be to be asked the terrestrial gravitation potential coefficient

$x_{n\×1}={\left(A_{t\×n}^T{A}_{t\×n}\right)}^{-1}\·A_{t\×n}^T{D}_{t\×1}$

Order ${\tilde{A}}_{t\×n}={\left(A_{t\×n}^T{A}_{t\×n}\right)}^{-1}\·A_{t\×n}^T,$ Therefore the following formula deformable does

$x_{n\×1}={\tilde{A}}_{t\×n}{D}_{t\×1}.$

When finding the solution large-scale system of linear equations, resolving fast of

is core link.As shown in Figure 4; Because directly finding the solution and storing

is all difficult; Therefore in order effectively to resolve GRACE-type Satellite Tracking satellite mode observation equation and then fast inversion high precision and high spatial resolution earth gravity field, the present invention proposes the aforementioned new thought of utilizing double star space three-dimensional method of interpolation accurately and fast to resolve .

The asterisk line is represented the precision of the 120 rank EIGEN-GRACE02S whole world gravity field model that German GFZ announces among Fig. 5, and accumulative total geoid surface precision is 18.938cm at place, 120 rank; Dotted line representes to utilize the precision of double star space three-dimensional method of interpolation inverting GRACE earth gravity field, and accumulative total geoid surface precision is 15.421cm at place, 120 rank.Can know that in the accordance at each place, rank double star space three-dimensional method of interpolation is one of effective ways of resolving high precision and high spatial resolution earth gravity field through 2 curves.The precision that the present invention is based on double star space three-dimensional method of interpolation inverting GRACE earth gravity field slightly is superior to the precision of EIGEN-GRACE02S whole world gravity field model in long wave part (L≤50 rank); And meet with it better in middle long wave part (50＜L≤120 rank), the analysis of causes is following: German GFZ utilized 110 days based on KINETIC METHOD and the GRACE moonscope data solver of SI 5s 120 terrace gravity fields.Can know according to asterisk line among Fig. 5; The precision of EIGEN-GRACE02S model long wave gravity field is not very good; Its main cause is because GRACE-A/B adopts double star to follow the tracks of difference modes (interstellar distance 220 ± 50km) each other; Balancing out between double star in the common error, also difference has been fallen part signal, therefore causes part order gravitation potential of earth coefficient susceptibility is reduced.In addition,, after effectively handling, still exist and morely can not know and uncontrollable error, therefore can cause earth gravity field long wave part signal to noise ratio (S/N ratio) lower because observation data is complicated.The present invention is based on the GRACE actual parameter and obtained the moonscope data; Crucial load error (K wave band measuring instrument, GPS receiver, accelerometer, Star Sensor etc.), atmospherical drag error, solid tide error, the sun and lunar gravitation error etc. have mainly been considered; But can not can not knowing at present with uncontrollable error of existing in the observation data all be taken into account, therefore cause the precision (L≤50 rank) of gravity field long-wave signal slightly to be superior to the precision of EIGEN-GRACE02S model.Through numerical simulation check of the present invention, measurable future, the precision of GRACE earth gravity field long-wave signal still had the potentiality that can improve along with the further optimization and the improvement of observation data disposal route and earth gravity field inversion method.In addition; The GRACE moonscope error information that we utilize U.S. JPL to announce based on semi analytical method has been estimated the precision of 120 rank GRACE earth gravity fields; Still there is difference to a certain degree in its result with the EIGEN-GRACE02S model in 30 rank, so the difference of earth gravity field inversion method also is to cause one of reason of resolving result difference.

Claims (1)

1. Satellite gravity inversion method based on double star space three-dimensional interpolation principle, its characteristic is following:

Step 1: utilize GRACE satellite K wave band measuring instrument to obtain speed between star, utilize the GPS receiver to obtain satellite orbital position and speed, utilize accelerometer to obtain the satellite nonconservative force and utilize Star Sensor to obtain the satellite 3 d pose; Calculate the earth disturbing potential T (x at satellite orbit place through above-mentioned DATA REASONING; Y, z);

Step 2: with the earth's core is that the centre of sphere is selected 4 equidistant and regular reference sphere r ₁, r ₂, r ₃And r ₄, be r with reference to the radius of sphere recently apart from the earth's core ₁=6830km, sphere are spaced apart Δ r=20km, orbit altitude H=500km, and satellite orbit is positioned at r ₂And r ₄Between sphere, on each sphere, carry out uniform grid and divide according to resolution=latitude Δ θ * longitude Δ λ, wherein Δ θ is taken as 180 °/2048, and Δ λ is taken as 360 °/4096;

Step 3: utilize each earth disturbing potential T of measuring on the satellite orbit in the step 1 (x, y, z) three-dimensional interpolation obtain 4 earth disturbing potential T with reference to corresponding cells grid on the sphere (r, θ, λ); The resolution of unit cell is designed to Δ θ * Δ λ * Δ r, and wherein Δ θ is taken as 180 °/2048, and Δ λ is taken as 360 °/4096, and Δ r is taken as 20km;

Step 4: latitude θ is divided into 2048 equal portions, and fixing each rule is with reference to latitude grid dividing value limited on the sphere; Secondly; Calculate the Legendre functional value and the storage of each minizone midpoint; Legendre functional value on other aspect approaches acquisition with the functional value on the known point as Taylor expansion, and the result of calculation storage is called in order to the back;

Step 5: utilize space three-dimensional be interpolated into the rule 4 with reference to the earth disturbing potential T on the sphere (r, θ λ) resolve large-scale overdetermined equation group (1), finally obtain the terrestrial gravitation potential coefficient;

$T(r,\mathrm{\θ},\mathrm{\λ})=\frac{\mathrm{GM}}{R_e}\underset{l=2}{\overset{L}{\mathrm{\Σ}}}{\left(\frac{R_e}{r}\right)}^{l+1}\underset{m=0}{\overset{l}{\mathrm{\Σ}}}({\stackrel{\‾}{C}}_{\mathrm{lm}}\cos \mathrm{m\λ}+{\stackrel{\‾}{S}}_{\mathrm{lm}}\sin \mathrm{m\λ}){\stackrel{\‾}{P}}_{\mathrm{lm}}\left(\cos \mathrm{\θ}\right)---\left(1\right)$

Wherein, GM representes that earth quality M and gravitational constant G are long-pending; R _eThe mean radius of the expression earth;

The earth's core radius of expression satellite, x, y, z represent three components of satellite position vector r respectively; θ and λ represent the geocentric colatitude degree and the geocentric longitude of satellite respectively;

Represent normalized Legendre function, l representes exponent number, and m representes number of times; L representes that gravitation potential of earth presses the maximum order of spherical-harmonic expansion;

With

Represent normalization terrestrial gravitation potential coefficient to be asked;

Interpolation formula in the said step 3 is:

$T(x,y,z)=\underset{i_x=1}{\overset{4}{\mathrm{\Σ}}}\underset{i_y=1}{\overset{4}{\mathrm{\Σ}}}\underset{i_z=1}{\overset{4}{\mathrm{\Σ}}}{w}_{i_{x}x}{w}_{i_{y}y}{w}_{i_{z}z}T(x_{i_x},y_{i_y},z_{i_z})$

$=\underset{i_x=1}{\overset{4}{\mathrm{\Σ}}}{w}_{i_{x}x}\underset{i_y=1}{\overset{4}{\mathrm{\Σ}}}{w}_{i_{y}y}\underset{i_z=1}{\overset{4}{\mathrm{\Σ}}}[w_{1z}T(x_{i_x},y_{i_y},z_1)+w_{2z}T(x_{i_x},y_{i_y},z_2)+w_{3z}T(x_{i_x},y_{i_y},z_3)+w_{4z}T(x_{i_x},y_{i_y},z_4)]---\left(2\right)$

$=\underset{i_x=1}{\overset{4}{\mathrm{\Σ}}}{w}_{i_{x}x}[w_{1y}B(x_{i_x},y_1)+w_{2y}B(x_{i_x},y_2)+w_{3y}B(x_{i_x},y_3)+w_{4y}B(x_{i_x},y_4)]$

$=w_{1x}C\left(x_1\right)+w_{2x}C\left(x_2\right)+w_{3x}C\left(x_3\right)+w_{4x}C\left(x_4\right)$

Wherein, (x, y, z) ∈ [x ₂, x ₃] * [y ₂, y ₃] * [z ₂, z ₃],

$C\left(x_{i_x}\right)=w_{1y}B(x_{i_x},y_1)+w_{2y}B(x_{i_x},y_2)+w_{3y}B(x_{i_x},y_3)+w_{4y}B(x_{i_x},y_4),$

$B(x_{i_x},y_{i_y})=w_{1z}T(x_{i_x},y_{i_y},z_1)+w_{2z}T(x_{i_x},y_{i_y},z_2)+w_{3z}T(x_{i_x},y_{i_y},z_3)+w_{4z}T(x_{i_x},y_{i_y},z_4),$

T (x, y, the z) earth disturbing potential on the expression satellite orbit,

Expression is with reference to the rectangular coordinate on the sphere, w _AbExpression with the earth disturbing potential three-dimensional interpolation on the satellite orbit to rule with reference to the weight coefficient on the sphere, a=1 wherein, 2,3,4, $w_{1b}=-\frac{1}{2}p{(1-p)}^2,$ $w_{2b}=1-\frac{5}{2}{p}^2+\frac{3}{2}{p}^3,$ $w_{3b}=\frac{1}{2}p(1+4p-3p^2),$ $w_{4b}=\frac{1}{2}{p}^2(p-1),$ $p=\frac{b-b_2}{\mathrm{\Δ}{q}_n},$ B=x, y, z represent the three-dimensional rectangular coordinate on the satellite orbit, b ₂=x ₂, y ₂, z ₂Represent second with reference to the three-dimensional rectangular coordinate on the sphere, Δ q _nSatellite orbit and second are with reference to the distance between sphere in expression; With the weight coefficient of double star space three-dimensional interpolation formula with matrix

Expression, then double star space three-dimensional interpolation formula matrix form is represented as follows:

$T(x,y,z)=W(w_{i_{x}x},w_{i_{y}y},w_{i_{z}z})\·T(x_{i_x},y_{i_y},z_{i_z})---\left(3\right)$

Therefore can pass through after the double star space three-dimensional interpolation 4 with reference to the disturbing potential on the sphere

$T(x_{i_x},y_{i_y},z_{i_z})={\left(W^{T}W\right)}^{-1}{W}^{T}T(x,y,z)---\left(4\right)$

Resolving large-scale overdetermined equation group in the said step 5 is:

(1) formula is expressed as with matrix form

D _t×1＝A _t×n·x _n×1 (5)

Wherein, D _{T * 1}(z) space three-dimensional is interpolated into reference to the column matrix behind the sphere earth disturbing potential T that expression is measured the satellite orbit place based on double star space three-dimensional method of interpolation for x, y, and t representes the quantity of observation station; x _{N * 1}Expression waits to ask the gravitation potential of earth coefficient vector

Wherein the gravitation potential coefficient is arranged according to exponent number l and is formed, and number of times m fixes,

Expression waits to ask the number of gravitation potential coefficient, L _MaxThe maximum order that launch by spheric function expression gravitation position; A _{T * n}The design matrix of the capable n row of expression t, said design matrix are illustrated in earth disturbing potential moonscope equation (1) the right except the terrestrial gravitation potential coefficient

Other outer part; Take advantage of together on the following formula both sides

The terrestrial gravitation potential coefficient that acquisition is to be asked

$x_{n\×1}={\left(A_{t\×n}^T{A}_{t\×n}\right)}^{-1}\·A_{t\×n}^T{D}_{t\×1}---\left(6\right).$

Images