CN103163562B - 基于滤波原理的卫星重力梯度反演方法 - Google Patents

Abstract

本发明涉及一种地球重力场精密测量方法，特别是一种基于滤波原理的卫星重力梯度反演方法；通过滤波原理建立新型卫星重力梯度观测方程，进而精确和快速反演地球重力场；该方法卫星重力梯度反演精度高，地球重力场计算速度快，易于卫星重力梯度系统敏感度分析，卫星观测方程物理含义明确，计算机性能要求低；滤波卫星重力梯度反演法是解算高精度和高空间分辨率地球重力场的有效方法。

Description

基于滤波原理的卫星重力梯度反演方法

技术领域

本发明涉及卫星重力梯度学、空间大地测量学、地球物理学、航空航天等交叉技术领域，特别是涉及一种通过滤波原理建立新型卫星重力梯度观测方程，进而精确和快速反演地球重力场的方法。

背景技术

地球重力场及其时变反映地球表层及内部物质的空间分布、运动和变化，同时决定着大地水准面的起伏和变化，因此确定地球重力场的精细结构及其时变不仅是大地测量学、地震学、海洋学、空间科学、国防建设等的需求，同时也将为全人类寻求资源、保护环境和预测灾害提供重要的信息资源。

GOCE（Gravity Field and Steady-State Ocean Circulation Explorer）重力梯度卫星由欧洲空间局（ESA）独立研制，已于2009年3月17日发射升空，主要用于地球中短波重力场的精密测量。GOCE采用近圆（轨道离心率0.001）、近极地（轨道倾角96.5°）和太阳同步轨道，经过20个月的飞行计划，轨道高度由250km降为240km。GOCE采用卫星跟踪卫星高低（SST-HL）和卫星重力梯度（SGG）模式的结合，除基于高轨道的GPS/GLONASS卫星对低轨道的GOCE卫星进行精密跟踪定位（1cm）外，利用定位于质心处的重力梯度仪（3×10^-12/s²）高精度测量卫星轨道高度处引力位的二阶导数，同时基于无阻尼离子微推进器补偿卫星受到的非保守力。早在20世纪80年代，国外便开始制定国际卫星重力梯度计划。由于地球重力场信号随卫星轨道高度的增加而急剧衰减(R_e/r)^l，基于分析卫星轨道运动仅适合于精密确定中长波地球重力场，而卫星重力梯度技术可直接测定地球引力位的二次微分，进而在一定程度上抑制了地球重力场信号的衰减效应，因此卫星重力梯度测量有利于高精度感测地球中短波重力场信号。

在利用卫星重力观测数据反演地球重力场的众多方法中，时域法的优点是直接对卫星观测数据进行处理，不需作任何近似，求解精度较高且能有效处理色噪声；缺点是随着卫星观测数据的增多，观测方程数量剧增，极大地增加了计算量。为了满足下一代卫星重力梯度测量计划中精密和快速解算高阶地球重力场的要求，本发明紧跟国际卫星重力梯度测量的热点和动态，利用新型滤波卫星重力梯度法，精确和快速反演了250阶GOCE地球重力场。

在现有技术中已有利用卫星重力梯度对角张量（V_xx,V_yy,V_zz），通过4个参考球面，在未考虑卫星重力梯度数据滤波情况下，恢复地球重力场的数值模拟方法。由于该方法未考虑卫星重力梯度交叉张量的贡献，而且未对卫星重力梯度观测数据进行有效滤波处理，因此地球重力场反演精度仍未达到预期要求。不同于已有技术，本发明基于卫星重力梯度全张量（V_xx,V_yy,V_zz,V_xy,V_xz,V_yz），通过5个参考球面，采用卫星重力梯度数据滤波处理技术，在不过多增加反演计算量的前提下，进一步提高了GOCE地球重力场的反演精度和计算速度。

发明内容

本发明的目的是：通过滤波原理，建立新型卫星重力梯度观测方程，进而精确和快速反演地球重力场。

为达到上述目的，本发明采用了如下技术方案：

一种基于滤波原理的卫星重力梯度反演方法，包括如下步骤：

步骤1：采集重力梯度卫星观测数据，其中通过重力梯度卫星的星载重力梯度仪采集卫星重力梯度全张量观测数据V_xx,V_yy,V_zz,V_xy,V_xz,V_yz，通过重力梯度卫星的星载GPS/GLONASS复合接收机采集卫星轨道位置观测数据r；

步骤2：建立卫星重力梯度观测模型，具体包括：

在地固系中，按球谐函数展开地球扰动位T(r,θ,λ)，θ和λ分别表示地心余纬度和地心经度，并分别对卫星位置矢量r的三个分量x,y,z进行二阶求导，在地心惯性系中，以矩阵方式表达地球扰动位T(r,θ,λ)的球谐函数展开式对三个分量x,y,z的二阶导数，以此建立卫星重力梯度观测方程作为卫星重力梯度观测模型，其中，y_g×1表示卫星轨道处的重力梯度观测数据，g表示重力梯度观测数据的个数；Γ_g×n表示g行n列的设计矩阵，n=L²+2L-3，L表示球谐函数展开的最大阶数；表示n×1列的待求地球引力位系数矩阵；

步骤3：使用滤波技术将卫星重力梯度观测模型处理为卫星重力梯度滤波观测模型，利用预处理共轭梯度迭代法求解卫星重力梯度滤波观测模型，进而反演地球重力场，具体包括：

以地心为球心选择若干个等间距且规则的参考球面，并在每个参考球面上进行均匀网格划分，在所划分的网格上利用卫星轨道上的重力梯度观测数据插值得到各个参考球面上对应单元格子的重力梯度观测数据；

对卫星重力梯度观测模型进行滤波变换，方程等式两边同时左乘滤波因子C_d的逆矩阵和设计矩阵的转置Γ^T以及预处理阵P_n×n的逆矩阵，得到卫星重力梯度滤波观测模型 $P_{n\×n}^1{G}_{n\×1}=P_{n\×n}^{-1}{N}_{n\×n}\·{\stackrel{\‾}{x}}_{n\×1},$ 其中 $G={\mathrm{\Γ}}^T{C}_d^{-1}y,$ $n={\mathrm{\Γ}}^T{C}_d^{-1}\mathrm{\Γ};$

求解获得滤波因子C_d；

利用预处理共轭梯度迭代法，快速解算卫星重力梯度滤波观测模型，获

得地球引力位系数和进而完成地球重力场反演。

本发明是基于滤波卫星重力梯度反演法有利于快速反演高精度和高空间分辨率地球重力场的特点而设计的，优点是：

1）卫星重力梯度反演精度高；

2）地球重力场计算速度快；

3）易于卫星重力梯度系统敏感度分析；

4）卫星观测方程物理含义明确；

5）计算机性能要求低。

附图说明

图1是GOCE重力梯度卫星的示意图。

图2表示基于滤波卫星重力梯度反演法解算地球重力场的原理图。

图3表示正规方阵N_n×n的块对角占优特性，色标代表矩阵元素数值的大小，色标条采用以10为底的对数表示。

图4表示基于滤波卫星重力梯度法反演地球引力位系数精度。

图5表示基于滤波卫星重力梯度法反演累计大地水准面精度。

具体实施方式

以下结合附图，对本发明的具体实施方式作进一步的说明。

基于滤波原理的卫星重力梯度反演方法包含下列步骤：

步骤一：重力梯度卫星观测数据采集

（1）通过重力梯度卫星的星载重力梯度仪采集卫星重力梯度全张量观测数据（V_xx,V_yy,V_zz,V_xy,V_xz,V_yz）。

（2）通过重力梯度卫星的星载GPS/GLONASS复合接收机采集卫星轨道位置观测数据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(1\right)$

其中，GM表示地球质量M和万有引力常数G之积，R_e表示地球的平均半径，L表示球函数展开的最大阶数；表示卫星的地心半径，x,y,z分别表示卫星位置矢量r的三个分量，θ和λ表示地心余纬度和地心经度；表示规格化的Legendre函数，l表示阶数，m表示次数；和表示待求的规格化地球引力位系数。

地球扰动位T(r,θ,λ)分别对x,y,z的二阶导数表示如下

其中，地球扰动位二阶导数是对称张量，同时在真空情况下满足Laplace方程表现为无迹性，T_xx+T_yy+T_zz=0，因此在9个重力梯度分量中有5个是独立的。全张量重力梯度的9个分量表示如下

$\left\{\begin{array}{c}{T}_{\mathrm{xx}}(r,\mathrm{\θ},\mathrm{\λ})=\frac{1}{r}{T}_r(r,\mathrm{\θ},\mathrm{\λ})+\frac{1}{r^2}{T}_{\mathrm{\θ\θ}}(r,\mathrm{\θ},\mathrm{\λ})\\ T_{\mathrm{yy}}(r,\mathrm{\θ},\mathrm{\λ})=\frac{1}{r}{T}_r(r,\mathrm{\θ},\mathrm{\λ})+\frac{1}{r^2}\cot {\mathrm{\θT}}_{\mathrm{\θ}}(r,\mathrm{\θ},\mathrm{\λ})+\frac{1}{r^2{\sin }^2\mathrm{\θ}}{T}_{\mathrm{\λ\λ}}(r,\mathrm{\θ},\mathrm{\λ})\\ T_{\mathrm{zz}}(r,\mathrm{\θ},\mathrm{\λ})=T_{\mathrm{rr}}(r,\mathrm{\θ},\mathrm{\λ})\\ T_{\mathrm{xy}}(r,\mathrm{\θ},\mathrm{\λ})=T_{\mathrm{yx}}(r,\mathrm{\θ},\mathrm{\λ})=\frac{1}{r^2\sin \mathrm{\θ}}[-\cot {\mathrm{\θT}}_{\mathrm{\λ}}(r,\mathrm{\θ},\mathrm{\λ})+T_{\mathrm{\θ\λ}}(r,\mathrm{\θ},\mathrm{\λ})]\\ T_{\mathrm{xz}}(r,\mathrm{\θ},\mathrm{\λ})=T_{\mathrm{zx}}(r,\mathrm{\θ},\mathrm{\λ})=\frac{1}{r^2}{T}_{\mathrm{\θ}}(r,\mathrm{\θ},\mathrm{\λ})-\frac{1}{r}{T}_{\mathrm{r\θ}}(r,\mathrm{\θ},\mathrm{\λ})\\ T_{\mathrm{yz}}(r,\mathrm{\θ},\mathrm{\λ})=T_{\mathrm{zy}}(r,\mathrm{\θ},\mathrm{\λ})=\frac{1}{r\sin \mathrm{\θ}}[\frac{1}{r}{T}_{\mathrm{\λ}}(r,\mathrm{\θ},\mathrm{\λ})-T_{\mathrm{r\λ}}(r,\mathrm{\θ},\mathrm{\λ})]\end{array}\right.---\left(3\right)$

其中，地球扰动位T(r,θ,λ)分别对r,θ,λ的一阶导数表示为

$\left\{\begin{array}{c}{T}_r(r,\mathrm{\θ},\mathrm{\λ})=-\frac{\mathrm{GM}}{R_e^2}\underset{l=2}{\overset{L}{\mathrm{\Σ}}}(l+1){\left(\frac{R_e}{r}\right)}^{l+2}\underset{m=0}{\overset{1}{\mathrm{\Σ}}}({\stackrel{\‾}{C}}_{\mathrm{lm}}\cos \mathrm{m\λ}+{\stackrel{\‾}{S}}_{\mathrm{lm}}\sin \mathrm{m\λ}){\stackrel{\‾}{P}}_{\mathrm{lm}}\left(\cos \mathrm{\θ}\right)\\ T_{\mathrm{\θ}}(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)\sin \mathrm{\θ}\\ T_{\mathrm{\λ}}(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{\Σ}}}m(-{\stackrel{\‾}{C}}_{\mathrm{lm}}\sin \mathrm{m\λ}+{\stackrel{\‾}{S}}_{\mathrm{lm}}\cos \mathrm{m\λ}){\stackrel{\‾}{P}}_{\mathrm{lm}}\left(\cos \mathrm{\θ}\right)\end{array}\right.---\left(4\right)$

地球扰动位T(r,θ,λ)分别对r,θ,λ的二阶导数表示为

$\left\{\begin{array}{c}{T}_{\mathrm{rr}}(r,\mathrm{\θ},\mathrm{\λ})=\frac{\mathrm{GM}}{R_e^3}\underset{l=2}{\overset{L}{\mathrm{\Σ}}}(l+1)(l+2){\left(\frac{R_e}{r}\right)}^{l+3}\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)\\ T_{\mathrm{\θ\θ}}(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){\sin }^2\mathrm{\θ}-{\stackrel{\‾}{P}}_{\mathrm{lm}}^{\′}\left(\cos \mathrm{\θ}\right)\cos \mathrm{\θ}]\\ T_{\mathrm{\λ\λ}}(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{\Σ}}}{m}^2({\stackrel{\‾}{C}}_{\mathrm{lm}}\cos \mathrm{m\λ}+{\stackrel{\‾}{S}}_{\mathrm{lm}}\sin \mathrm{m\λ}){\stackrel{\‾}{P}}_{\mathrm{lm}}\left(\cos \mathrm{\θ}\right)\\ R_{\mathrm{r\θ}}(r,\mathrm{\θ},\mathrm{\λ})=T_{\mathrm{\θr}}(r,\mathrm{\θ},\mathrm{\λ})=\frac{\mathrm{GM}}{R_e^2}\underset{l=2}{\overset{L}{\mathrm{\Σ}}}(l+1){\left(\frac{R_e}{r}\right)}^{l+2}\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)\sin \mathrm{\θ}\\ T_{\mathrm{r\λ}}(r,\mathrm{\θ},\mathrm{\λ})=T_{\mathrm{\λr}}(r,\mathrm{\θ},\mathrm{\λ})=\frac{\mathrm{GM}}{R_e^2}\underset{l=2}{\overset{L}{\mathrm{\Σ}}}(l+1){\left(\frac{R_e}{r}\right)}^{l+2}\underset{m=0}{\overset{l}{\mathrm{\Σ}}}m({\stackrel{\‾}{C}}_{\mathrm{lm}}\sin \mathrm{m\λ}-{\stackrel{\‾}{S}}_{\mathrm{lm}}\cos \mathrm{m\λ}){\stackrel{\‾}{P}}_{\mathrm{lm}}\left(\cos \mathrm{\θ}\right)\\ T_{\mathrm{\θ\λ}}(r,\mathrm{\θ},\mathrm{\λ})=T_{\mathrm{\λ\θ}}(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{\Σ}}}m({\stackrel{\‾}{C}}_{\mathrm{lm}}\sin \mathrm{m\λ}-{\stackrel{\‾}{S}}_{\mathrm{lm}}\cos \mathrm{m\λ}){\stackrel{\‾}{P}}_{\mathrm{lm}}^{\′}\left(\cos \mathrm{\θ}\right)\sin \mathrm{\θ}\end{array}\right.---\left(5\right)$

Legendre函数及一阶导数和二阶导数表示为

$\left\{\begin{array}{c}{\stackrel{\‾}{P}}_{\mathrm{lm}}\left(\cos \mathrm{\θ}\right)={\mathrm{\γ}}_m{2}^{-l}{\sin }^m\mathrm{\θ}\underset{k=0}{\overset{[(l-m)/2]}{\mathrm{\Σ}}}{(-1)}^k\frac{(2l-2k)!}{k!(l-k)!(l-m-2k)!}{\left(\cos \mathrm{\θ}\right)}^{l-m-2k}(m\≤l)\\ {\stackrel{\‾}{P}}_{\mathrm{lm}}^{\′}\left(\cos \mathrm{\θ}\right)=\frac{1}{\sin \mathrm{\θ}}[(l+1)\cos \mathrm{\θ}{\stackrel{\‾}{P}}_{\mathrm{lm}}\left(\cos \mathrm{\θ}\right)-(l-m-1){\stackrel{\‾}{P}}_{l+1,m}\left(\cos \mathrm{\θ}\right)]\\ {\stackrel{\‾}{P}}_{\mathrm{lm}}^{\′\′}\left(\cos \mathrm{\θ}\right)=-l{\stackrel{\‾}{P}}_{\mathrm{lm}}\left(\cos \mathrm{\θ}\right)+l\cos \mathrm{\θ}{\stackrel{\‾}{P}}_{l-1,m}^{\′}\left(\cos \mathrm{\θ}\right)+\frac{l}{4}{\cos }^2\mathrm{\θ}[{\stackrel{\‾}{P}}_{l-1,m+1}^{\′}\left(\cos \mathrm{\θ}\right)-4{\stackrel{\‾}{P}}_{l-1,m-1}^{\′}\left(\cos \mathrm{\θ}\right)]\end{array}\right.---\left(6\right)$

其中， ${\mathrm{\γ}}_m=\left\{\begin{array}{cc}\sqrt{2(2l+1)\frac{(l-|m\left|\right)!}{(l+|m\left|\right)!}}& (m\≠0)\\ \sqrt{2l+1}& (m=0)\end{array}\right..$

在地心惯性系中，将公式（3）的矩阵形式表示为卫星重力梯度观测方程

$y_{g\×1}={\mathrm{\Γ}}_{g\×n}{\stackrel{\‾}{x}}_{n\×1}---\left(7\right)$

其中，y_g×1表示卫星轨道处的重力梯度观测数据，g表示重力梯度观测数据的个数；Γ_g×n表示g行n列的设计矩阵，n=L²+2L-3；表示n×1列的待求地球引力位系数矩阵。

步骤三：卫星重力梯度反演

如图2所示，本发明基于滤波卫星重力梯度技术反演250阶GOCE地球重力场的主要思想如下：

第一，本发明以地心为球心选择5个等间距且规则的参考球面r₁，r₂，r₃，r₄和r₅，距地心最近参考球面的半径为r₁=6588km，球面间隔为Δr=20km，使得卫星轨道H=250km位于r₂和r₄参考球面之间，并在每个参考球面上按照分辨率＝纬度Δθ（180°/4096）×经度Δλ（360°/8192）进行均匀网格划分。参考球面个数最优选取的原则如下：在重力梯度插值误差对地球重力场解算精度不产生实质影响的标准下，尽可能采用较少的参考球面个数，进而减少整体计算量和降低对计算机性能的要求。本发明利用时间长度为180天和采样间隔为5s的卫星观测数据，并分别选取4个、5个、6个和7个参考球面解算了250阶GOCE地球重力场，结果表明：由于插值计算引入的误差较GOCE卫星关键载荷的测量精度约小10倍，不会对地球重力场反演精度产生实质影响，因此选取5个参考球面较优。

第二，利用卫星轨道上的1个重力梯度观测数据插值得到5个参考球面上1个单元格子的重力梯度观测数据，单元格子的分辨率设计为Δθ（0.05°）×Δλ（0.05°）×Δr（20km）。

第三，由于卫星重力梯度观测数据量非常大，而且C_d不是对角阵以及C_d的阶数很大，因此的计算非常困难。本发明将C_d分为只与观测到的某个重力梯度张量的分量对应的子矩阵，且每个子矩阵都是N_p×N_p的对称Toeplitz阵，N_p表示观测点的总数，因此，滤波处理过程得以较大程度简化。

第四，利用预处理共轭梯度迭代法，快速解算大型超定方程组（7）式，获得地球引力位系数和进而完成地球重力场反演。

卫星观测方程（7）是大型线性超定方程组，因此没有精确解，只有最小二乘解。在方程两边同时左乘得到

${\mathrm{\Γ}}^T{C}_d^{-1}y={\mathrm{\Γ}}^T{C}_d^{-1}\mathrm{\Γ}\·\stackrel{\‾}{x}---\left(8\right)$

其中，C_d为滤波因子。

令 $G={\mathrm{\Γ}}^T{C}_d^{-1}y,$ $N={\mathrm{\Γ}}^T{C}_d^{-1}\mathrm{\Γ},$ 则方程（8）可变形为

$G_{n\×1}=N_{n\×n}\·{\stackrel{\‾}{x}}_{n\×1}---\left(9\right)$

在方程（9）两边同时左乘得

$P_{n\×n}^{-1}{G}_{n\×1}=P_{n\×n}^{-1}{N}_{n\×n}\·{\stackrel{\‾}{x}}_{n\×1}---\left(10\right)$

其中，P_n×n表示预处理阵。

Γ_g×n是一个庞大的长方矩阵，存储需占用大量的内存空间，因此直接存储较难实现，而正规方阵N_n×n虽较Γ_g×n缩小了许多，但如果直接存储也会占用大量的内存空间（约占12Gb）。如图3所示，N_n×n是一个块对角占优的方阵，此性质为本发明迭代求解的加速提供了有利条件。

所谓滤波处理即是计算通常情况下，卫星重力梯度观测数据包含色噪声，因此C_d不是对角阵。由于卫星重力梯度观测数据量非常大，因此C_d的阶数很大，进而的计算非常困难。但如果对观测模型适当处理，则的计算便可看成是一个相对容易处理的滤波过程，此即的计算被称为滤波处理的原因。如果所有卫星重力梯度观测数据按采样时间的先后顺序计数，同时满足如下条件：（a）张量的不同分量互不相关；（b）噪声相对稳定，即观测点i和观测点j的相关性只由|i-j|来决定；（c）所有相邻观测点的采样间隔相同。在此前提下，C_d可被分为只与观测到的某个重力梯度张量的分量对应的子矩阵，且每个子矩阵都是N_p×N_p的对称Toeplitz阵，N_p表示观测点的总数，因此，滤波处理过程得以较大程度简化。另外，如果对应于|i-j|的相关性与对应于|i-j±N_p|的相关性等价，则上述每个子矩阵将变成更为简单的循环Toeplitz阵。循环Toeplitz阵有两个重要性质：一是它可仅由矩阵的第一行来生成，二是它的逆阵仍为循环Toeplitz阵，且原阵和逆阵的元素在Fourier域内满足如下简单关系

$F\left(C_d^{(..)-1}\right)=\frac{1}{N_p^{2}F\left(C_d^{(..)}\right)}---\left(11\right)$

其中，F(·)表示循环Toeplitz阵第一行元素的Fourier变换。由于基于卫星重力梯度观测数据可获得噪声功率谱密度的均方根信息因此，利用公式（11）在Fourier域内可直接获得滤波因子C_d。

由于卫星轨道位置和重力梯度观测数据中含有色噪声，而在使用最小二乘求解大型线性超定方程组的现有技术中并未考虑色噪声的影响（即相当于将公式（8）中的滤波因子矩阵简化为单位阵来进行计算），由于随着阶数l的增加，Γ^TΓ的病态性急剧增强，进而严重影响了地球重力场反演的精度。为了解决此技术问题，现有技术中提出了采用正则化处理来降低Γ^TΓ病态性的方法，例如采用Kaula正则化的方式来进行处理，这些正则化处理方法虽然在一定程度上缓解了Γ^TΓ的病态性随着阶数l急剧增加的问题，但是仍然降低了卫星重力梯度反演精度，而且增加了反演计算量。不同于已有技术中应对卫星轨道位置和重力梯度观测数据所包含的色噪声影响的方法，本发明使用基于滤波处理技术的求解方法，通过使用公式（8）中的滤波因子对这些色噪声进行有针对性的去除，在不过多增加反演计算量的前提下，进一步改善了地球重力场的反演精度和计算速度。

预处理共轭梯度迭代法（PCCG）是目前求解大型线性方程组的有效方法之一，主体思想如下：第一，每一步迭代均对待求参量进行修正，直到达到预期精度为止；第二，每一步迭代的方向选择以误差最小为原则；第三，回避最小二乘法的直接矩阵求逆，通过循环迭代求解真值。利用PCCG求解观测方程（10），不需要直接存储N_n×n，总运算量只需要900Mb的内存空间。PCCG最关键的部分是P_n×n的选取，标准如下：第一，与越接近越好，保证了大型线性方程组求解的精度；第二，易于计算，提高了大型线性方程组求解的速度。本发明选取N_n×n的块对角部分作为预处理阵，形成的P_n×n为主对角线上按次数m排列且其余部分为0的块对角方阵，如此选取不仅保留了N_n×n的主要特征，而且较易于计算。总而言之，适当选取预处理阵可极大地减少PCCG求解引力位系数中循环迭代的次数，较直接最小二乘法计算速度至少提高1000倍。本发明基于PCCG反演了250阶GOCE地球重力场，经过55步迭代，总计耗时为32小时。

图4表示经过55步迭代，基于滤波卫星重力梯度反演法，利用重力梯度全张量（V_xx,V_yy,V_zz,V_xy,V_xz,V_yz）反演250阶GOCE地球引力位系数的精度（由上而下），其中轨道误差为1cm，卫星重力梯度误差为3×10^-12/s²。在250阶处，反演地球引力位系数的精度为4.668×10^-11（第55步迭代）。

如图5所示，实线表示基于滤波卫星重力梯度反演法，利用重力梯度全张量（V_xx,V_yy,V_zz,V_xy,V_xz,V_yz），在第55步迭代处反演250阶GOCE累计大地水准面精度；在250阶处，累计大地水准面精度为5.163cm。

以上具体实施方式仅为本发明的一种实施示例，其描述较为具体和详细，但不能因此而理解为对本发明专利范围的限制。其具体实施步骤顺序和模型参数可根据实际需要进行相应的调整。应当指出的是，对于本领域的普通技术人员来说，在不脱离本发明构思的前提下，还可以做出若干变形和改进，这些都属于本发明的保护范围。

Claims (6)

1.一种基于滤波原理的卫星重力梯度反演方法，其特征在于包括如下步骤：

步骤1：采集重力梯度卫星观测数据，其中通过重力梯度卫星的星载重力梯度仪采集卫星重力梯度全张量观测数据V_xx,V_yy,V_zz,V_xy,V_xz,V_yz，通过重力梯度卫星的星载GPS/GLONASS复合接收机采集卫星轨道位置观测数据r；

步骤2：建立卫星重力梯度观测模型，具体包括：

在地固系中，按球谐函数展开地球扰动位T(r,θ,λ)，θ和λ分别表示地心余纬度和地心经度，并分别对卫星位置矢量的三个分量x,y,z进行二阶求导，在地心惯性系中，以矩阵方式表达地球扰动位T(r,θ,λ)的球谐函数展开式对三个分量x,y,z的二阶导数，以此建立卫星重力梯度观测方程作为卫星重力梯度观测模型，其中，y_g×1表示卫星轨道处的重力梯度观测数据，g表示重力梯度观测数据的个数；Γ_g×n表示g行n列的设计矩阵，n＝L²+2L-3，L表示球谐函数展开的最大阶数；表示n×1列的待求地球引力位系数矩阵；

步骤3：使用滤波技术将卫星重力梯度观测模型处理为卫星重力梯度滤波观测模型，利用预处理共轭梯度迭代法求解卫星重力梯度滤波观测模型，进而反演地球重力场，具体包括：

以地心为球心选择若干个等间距且规则的参考球面，并在每个参考球面上进行均匀网格划分，在所划分的网格上利用卫星轨道上的重力梯度观测数据插值得到各个参考球面上对应单元格子的重力梯度观测数据；

对卫星重力梯度观测模型进行滤波变换，方程等式两边同时左乘滤波因子C_d的逆矩阵和设计矩阵的转置Γ^T以及预处理阵P_n×n的逆矩阵，得到卫星重力梯度滤波观测模型 $P_{n\×n}^{-1}{G}_{n\×1}=P_{n\×n}^{-1}{N}_{n\×n}\·{\stackrel{\‾}{x}}_{n\×1},$ 其中 $G={\mathrm{\Γ}}^T{C}_d^{-1}y,$ $N={\mathrm{\Γ}}^T{C}_d^{-1}\mathrm{\Γ};$

求解获得滤波因子C_d；

利用预处理共轭梯度迭代法，快速解算卫星重力梯度滤波观测模型，获得地球引力位系数和进而完成地球重力场反演。

2.如权利要求1所述的基于滤波原理的卫星重力梯度反演方法，其特征在于：所述步骤3中求解获得滤波因子C_d的方法为：

将卫星重力梯度全张量观测数据V_xx,V_yy,V_zz,V_xy,V_xz,V_yz按采样时间的先后顺序计数，同时满足如下条件：(a)张量的不同分量互不相关；(b)噪声相对稳定，即观测点i和观测点j的相关性只由|i-j|来决定；(c)所有相邻观测点的采样间隔相同；将C_d分为只与观测到的某个重力梯度张量的分量对应的子矩阵，且每个子矩阵都是N_p×N_p的对称Toeplitz阵，N_p表示观测点的总数。

3.如权利要求2所述的基于滤波原理的卫星重力梯度反演方法，其特征在于：所述步骤3中求解获得滤波因子C_d的方法为：

使得观测点i和观测点j的对应于|i-j|的相关性与对应于|i-j±N_p|的相关性等价，则C_d分为的每个子矩阵变为循环Toeplitz阵，利用循环Toeplitz阵在Fourier域的性质，在Fourier域内直接求解获得滤波因子C_d。

4.如权利要求1-3任意一项所述的基于滤波原理的卫星重力梯度反演方法，其特征在于：所述步骤3中选择若干个参考球面的方法为：

以地心为球心选择5个等间距且规则的参考球面r₁，r₂，r₃，r₄和r₅，距地心最近参考球面的半径为r₁＝6588km，球面间隔为Δr＝20km，使得卫星轨道H＝250km位于r₂和r₄参考球面之间。

5.如权利要求4所述的基于滤波原理的卫星重力梯度反演方法，其特征在于：所述步骤3中在每个参考球面上进行均匀网格划分，利用卫星轨道上的重力梯度全张量观测数据插值得到各个参考球面上每个单元格子的重力梯度观测数据的方法为：

在每个参考球面上进行均匀网格划分；将卫星轨道上的重力梯度全张量观测数据插值得到每个参考球面上对应单元格子的重力梯度观测数据，单元格子的分辨率为Δθ(0.05°)×Δλ(0.05°)×Δr(20km)。

6.如权利要求1-5任意一项所述的基于滤波原理的卫星重力梯度反演方法，其特征在于：所述重力梯度卫星为GOCE卫星。