基于重力梯度误差模型原理的卫星重力反演方法

技术领域

本发明涉及卫星重力梯度学、大地测量学、地球物理学、航空航天等交叉 技术领域，特别是涉及一种通过分析卫星重力梯度垂向张量误差、水平张量误 差和相关张量误差对累计大地水准面精度的联合影响建立新型重力梯度误差模 型，进而精确和快速反演地球重力场的方法。

背景技术

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

不同于GRACE（Gravity Recovery and Climate Experiment）双星高精度感测 地球中长波重力场，欧空局（ESA）提出了专用于地球中短波重力场精密探测的 GOCE（Gravity Field and Steady-State Ocean Circulation Explorer）卫星重力梯度 （SGG）计划。如图1所示，GOCE卫星已于2009年3月17日成功发射升空， 采用近圆（轨道离心率0.001）、极地（轨道倾角96.5°）和太阳同步轨道，经过 3年的飞行计划，轨道高度由250km降为240km。GOCE卫星采用卫星跟踪卫 星高低模式和卫星重力梯度模式的结合（SST-HL/SGG），除基于高轨道的GPS 和GLONASS卫星对低轨道的GOCE进行精密跟踪定位（定轨精度1cm），同 时利用定位于卫星质心处的星载重力梯度仪（测量精度3×10^-12/s²）高精度测量 卫星轨道高度处引力位的二阶导数。GOCE采用了非保守力补偿技术 （Drag-free），首先利用重力梯度仪测量由非保守力（大气阻力、太阳光压、地 球辐射压、轨道高度和姿态控制力等）引起的卫星质心的线性加速度与卫星平 台的角加速度；最后，结合卫星平台姿态测量数据，通过无阻尼离子微推进器 补偿卫星受到的非保守力。由于卫星重力梯度测量数据中的非保守力效应得到 了有效扣除，因此进一步提高了地球重力场反演的精度和空间分辨率。自20世 纪初匈牙利物理学家设计出第一台重力梯度仪（扭称）以来， 重力梯度仪经历了从单轴旋转到三轴定向，从室温到低温（低于4.2k），从扭力、 静电悬浮、超导到冷原子干涉的发展过程，测量精度日益提高。由于地球重力 场信号随卫星轨道高度的增加而急剧衰减(R_e/r)^l，基于分析卫星轨道运动仅适 合于确定地球中长波重力场，而卫星重力梯度是直接测定地球引力位的二次微 分，其结果将球谐系数放大了l²倍，因此可有效抑制地球引力位随高度的衰减 效应，进而高精度感测地球中高频重力场信号。欧空局独立研制的GOCE卫星 原计划于2004年6月发射，由于星载三维静电悬浮重力梯度仪未能达到预期精 度指标3×10^-12/s²（单个加速度计的分辨率超过10^-13m/s²，较GRACE卫星加速 度计分辨率高约3个数量级）以及重力梯度卫星整体系统研制的困难性，因此 距成功发射为止已推迟至少6次之多。

基于GOCE卫星重力梯度测量计划预计于2014年前结束，以及为了进一步 提高地球重力场中短波信号的探测精度，目前国际众多科研机构正积极推动 GOCE Follow-On卫星重力梯度测量计划的成功实施。我国相关研究机构紧跟国 际卫星重力梯度测量的热点和动态，正积极投身于卫星重力梯度测量计划的需 求论证和载荷预研之中。目前，国内外众多学者在基于卫星重力梯度技术反演 地球重力场的理论和方法等方面已开展了广泛研究。不同于前人的研究，本发 明通过分析卫星重力梯度垂向张量误差、水平张量误差和相关张量误差对累计 大地水准面精度的联合影响建立新型重力梯度误差模型，进而精确和快速反演 地球重力场。本发明不仅可为我国下一代地球重力场模型的精确建立提供理论 基础和技术保证，同时对卫星重力梯度反演技术的发展方向具有一定借鉴意义。

在现有技术中已有通过分析一维垂向分量V_zz和三维全张量V_ij重力梯度对 累计大地水准面精度的影响，建立卫星重力梯度解析误差模型，进而精确和快 速估计GOCE地球重力场精度的方法，该方法假设全球均匀分布的重力梯度观 测值有N₀个，并假设N₀个观测值的误差满足正态分布随机特性，大量数据的平 均可有效降低噪声，因此地球引力位系数的方差正比于1/N₀。同时该方法未考 虑垂向重力梯度和水平重力梯度之间所存在的较强相关性，未考虑基于相关性 重力梯度误差反演累计大地水准面精度，因此其反演精度仍未达到预期的要求。

发明内容

本发明的目的是：通过分析卫星重力梯度垂向张量误差、水平张量误差和 相关张量误差对累计大地水准面精度的联合影响，建立新型重力梯度误差模型， 进而精确和快速反演地球重力场。

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

一种基于重力梯度误差模型的卫星重力反演方法，包括如下步骤：

步骤1，通过重力梯度卫星的星载重力梯度仪采集卫星重力梯度误差数据δT_xyz；

步骤2，建立重力梯度误差模型，具体包括：

步骤2.1，在地固系中，将地球扰动位T(r,θ，λ)按球谐函数展开，分别对重 力梯度卫星轨道位置矢量r的三个分量x,y,z进行二阶求导，其中θ和λ分别表示 重力梯度卫星的地心余纬度和地心经度，将地球扰动位T(r,θ，λ)分别对r,θ，λ进 行一阶求导和二阶求导；

步骤2.2，通过地球扰动位T(r,θ，λ)分别对x,y,z的二阶导数和分别对r,θ，λ 一阶导数、二阶导数以及Legendre函数及其一阶导数和二阶导数得到待求规格 化地球引力位系数和与一维垂向重力梯度T_zz的关系，将此关系表示为垂 向重力梯度公式；

步骤2.3，利用所得到的垂向重力梯度公式以及地球引力位系数精度 与累计大地水准面精度的关系获得基于垂向重力梯度误差数据δT_zz反 演累计大地水准面精度的误差模型；

步骤2.4，基于球谐函数的正交性，得到待求规格化地球引力位系数和与一维水平重力梯度T_xx和T_yy的关系，将此关系表示为水平重力梯度公式，利用 所得到的水平重力梯度公式以及地球引力位系数精度与累计大地水准 面精度的关系分别获得基于水平重力梯度误差数据δT_xx和δT_yy反演累计大地水 准面精度的误差模型；

步骤2.5，根据垂向重力梯度T_zz与水平重力梯度T_xx和T_yy三者之间的非相互 独立性，联合垂向重力梯度误差数据δT_zz反演累计大地水准面精度的误差模型与 水平重力梯度误差数据δT_xx和δT_yy反演累计大地水准面精度的误差模型，得到基 于相关性重力梯度误差数据δT_z-x-y反演累计大地水准面精度的误差模型，忽略 非对角张量对地球重力场精度的影响，得到基于卫星重力梯度全张量误差数据 δT_xyz反演累计大地水准面精度的误差模型，以此作为重力梯度误差模型；

步骤3，基于所述重力梯度误差模型反演地球重力场精度，具体包括：

步骤3.1，在地球表面绘制网格；其次，按照重力梯度卫星的轨道在地球表 面的轨迹点位置依次加入卫星重力梯度误差数据δT_xyz，将分布于地球表面的卫 星重力梯度误差数据δT_xyz平均归算于划分的网格点δT_xyz(φ，λ)处；

步骤3.2，将δT_xyz(φ，λ)按球谐函数展开，得到δT_xyz(φ，λ)按球谐函数展开的 系数与卫星重力梯度误差数据在各阶处方差之间的关系，将 卫星重力梯度误差数据在各阶处方差代入所述重力梯度误差模型进而 反演地球重力场。

本发明是基于新型重力梯度误差模型法有利于快速反演高精度和高空间分 辨率地球重力场的特点而设计的，优点为：

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

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

3）易于重力梯度卫星系统误差分析；

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

5）计算机性能要求低。

附图说明

图1表示GOCE卫星重力梯度计划。

图2表示基于新型重力梯度误差模型法反演累计大地水准面精度对比。

具体实施方式

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

基于重力梯度误差模型原理的卫星重力反演方法包含下列步骤：

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

通过重力梯度卫星的星载重力梯度仪采集卫星重力梯度误差数据δT_xyz。

步骤二：重力梯度误差模型建立

在地固系中，地球扰动位T(r,θ，λ)按球谐函数展开的表达式为

$T(r,\theta ,\lambda )=\frac{\mathrm{GM}}{R_e}\underset{l=2}{\overset{L}{\Sigma }}{\left(\frac{R_e}{r}\right)}^{l+1}\underset{m=0}{\overset{l}{\Sigma }}({\bar{C}}_{\mathrm{lm}}\cos \mathrm{m\lambda }+{\bar{S}}_{\mathrm{lm}}\sin \mathrm{m\lambda }){\bar{P}}_{\mathrm{lm}}\left(\cos \theta \right)---\left(1\right)$

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

T(r,θ，λ)分别对x,y,z的二阶导数表示为

$\frac{{\partial }^{2}T}{\partial x\partial y}=\left(\begin{array}{ccc}{T}_{\mathrm{xx}}& T_{\mathrm{xy}}& T_{\mathrm{xz}}\\ T_{\mathrm{yx}}& T_{\mathrm{yy}}& T_{\mathrm{yz}}\\ T_{\mathrm{zx}}& T_{\mathrm{zy}}& T_{\mathrm{zz}}\end{array}\right)---\left(2\right)$

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

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

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

$\left(\begin{array}{c}{T}_r(r,\theta ,\lambda )=-\frac{\mathrm{GM}}{R_e^2}\underset{l=2}{\overset{L}{\Sigma }}(l+1){\left(\frac{R_e}{r}\right)}^{l+2}\underset{m=0}{\overset{l}{\Sigma }}({\bar{C}}_{\mathrm{lm}}\cos \mathrm{m\lambda }+{\bar{S}}_{\mathrm{lm}}\sin \mathrm{m\lambda }){\bar{P}}_{\mathrm{lm}}\left(\cos \theta \right)\\ T_{\theta }(r,\theta ,\lambda )=-\frac{\mathrm{GM}}{R_e}\underset{l=2}{\overset{L}{\Sigma }}{\left(\frac{R_e}{r}\right)}^{l+1}\underset{m=0}{\overset{l}{\Sigma }}({\bar{C}}_{\mathrm{lm}}\cos \mathrm{m\lambda }+{\bar{S}}_{\mathrm{lm}}\sin \mathrm{m\lambda }){\bar{P}}_{\mathrm{lm}}^{\prime }\left(\cos \theta \right)\sin \theta \\ T_{\lambda }(r,\theta ,\lambda )=\frac{\mathrm{GM}}{R_e}\underset{l=2}{\overset{L}{\Sigma }}{\left(\frac{R_e}{r}\right)}^{l+1}\underset{m=0}{\overset{l}{\Sigma }}m(-{\bar{C}}_{\mathrm{lm}}\sin \mathrm{m\lambda }+{\bar{S}}_{\mathrm{lm}}\cos \mathrm{m\lambda }){\bar{P}}_{\mathrm{lm}}\left(\cos \theta \right)\end{array}\right)---\left(4\right)$

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

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

$\left(\begin{array}{c}{\bar{P}}_{\mathrm{lm}}\left(\cos \theta \right)={\gamma }_m{2}^{-l}{\sin }^m\theta \underset{k=0}{\overset{[(l-m)/2]}{\Sigma }}{(-1)}^k\frac{(2l-2k)!}{k!(l-k)!(l-m-2k)!}{\left(\cos \theta \right)}^{l-m-2k}(m\le l)\\ {\bar{P}}_{\mathrm{lm}}^{\prime }\left(\cos \theta \right)=\frac{1}{\sin \theta }[(l+1)\cos \theta {\bar{P}}_{\mathrm{lm}}\left(\cos \theta \right)-(l-m-1){\bar{P}}_{l+1,m}\left(\cos \theta \right)]\\ {\bar{P}}_{\mathrm{lm}}^{\prime \prime }\left(\cos \theta \right)=-l{\bar{P}}_{\mathrm{lm}}\left(\cos \theta \right)+l\cos \theta {\bar{P}}_{l-1,m}^{\prime }\left(\cos \theta \right)+\frac{l}{4}{\cos }^2\theta [{\bar{P}}_{l-1,m+1}^{\prime }\left(\cos \theta \right)-4{\bar{P}}_{l-1,m-1}^{\prime }\left(\cos \theta \right)\end{array}\right)---\left(6\right)$

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

基于球谐函数的正交性，联合公式（3）和公式（5）可得一维垂向重力梯 度公式

$({\bar{C}}_{\mathrm{lm}},{\bar{S}}_{\mathrm{lm}})=\frac{R_e^3}{4\mathrm{\pi GM}}{\left(\frac{r}{R_e}\right)}^{l+3}{(l+1)}^{-1}{(l+2)}^{-1}\underset{\sigma }{\int \int }{T}_{\mathrm{zz}}{\bar{Y}}_{\mathrm{lm}}(\theta ,\lambda )\mathrm{d\sigma }---\left(7\right)$

其中，T_zz表示一维垂向重力梯度，实际计算时需要离散化数值积分。基于等间 隔的Δθ和Δλ在地球表面进行全球经纬网划分，同时将每个网格中的垂向重力 梯度值取平均值其中i,j表示网格的经纬度标号。因此，公式（7）可改写 为

$({\bar{C}}_{\mathrm{lm}},{\bar{S}}_{\mathrm{lm}})=\frac{R_e^3}{4\mathrm{\pi GM}}{\left(\frac{r}{R_e}\right)}^{l+3}{(l+1)}^{-1}{(l+2)}^{-1}\underset{i,j}{\Sigma }{\bar{T}}_{\mathrm{zz}}{|}_{\mathrm{ij}}\underset{{\sigma }_{\mathrm{ij}}}{\int \int }{\bar{Y}}_{\mathrm{lm}}(\theta ,\lambda )d{\sigma }_{\mathrm{ij}}---\left(8\right)$

累计大地水准面精度公式表示为

${\sigma }_N^L=R_e\sqrt{\underset{l=2}{\overset{L}{\Sigma }}\underset{m=0}{\overset{l}{\Sigma }}{\left(\delta {\bar{C}}_{\mathrm{lm}}\right)}^2+{\left(\delta {\bar{S}}_{\mathrm{lm}}\right)}^2}---\left(9\right)$

其中，表示地球引力位系数精度。

联合公式（8）和公式（9），可得基于一维垂向重力梯度误差数据δT_zz反演 累计大地水准面精度的误差模型

${\sigma }_N\left(T_{\mathrm{zz}}\right)\approx \frac{R_e}{\mathrm{GM}/R_e^3}\sqrt{\underset{l=2}{\overset{L}{\Sigma }}\frac{2l+1}{{(l+1)}^2{(l+2)}^2}{\left(\frac{r}{R_e}\right)}^{2(l+3)}{\sigma }_l^2\left({\mathrm{\delta T}}_{\mathrm{zz}}\right)}---\left(10\right)$

基于球谐函数的正交性，联合公式（3）~（5）可得水平方向重力梯度T_xx(yy)公式

$({\bar{C}}_{\mathrm{lm}},{\bar{S}}_{\mathrm{lm}})=\frac{R_e^3}{4\mathrm{\pi GM}}{\left(\frac{r}{R_e}\right)}^{l+3}{(l+1)}^{-1}{(m-l-1)}^{-1}\underset{\sigma }{\int \int }{T}_{\mathrm{xx}\left(\mathrm{yy}\right)}{\bar{Y}}_{\mathrm{lm}}(\theta ,\lambda )\mathrm{d\sigma }---\left(11\right)$

联合公式（9）和公式（11），可得分别基于水平重力梯度误差数据δT_xx和δT_yy反演累计大地水准面精度的误差模型

${\sigma }_N\left(T_{\mathrm{xx}}\right)={\sigma }_N\left(T_{\mathrm{yy}}\right)\approx \frac{R_e}{\mathrm{GM}/R_e^3}\sqrt{\underset{l=2}{\overset{L}{\Sigma }}\frac{2l+1}{\frac{4{(l+1)}^3(l+2)(l+3)}{9(2l+1)}}{\left(\frac{r}{R_e}\right)}^{2(l+3)}{\sigma }_l^2\left({\mathrm{\delta T}}_{\mathrm{xx}}\right)}---\left(12\right)$

由于垂向重力梯度T_zz以及水平重力梯度T_xx和T_yy非相互独立，而具有较强的 相关性，因此，联合公式（10）和公式（12）可得基于相关性重力梯度误差数 据δT_z-x-y反演累计大地水准面精度的误差模型

${\sigma }_N(T_{\mathrm{zz}}-T_{\mathrm{xx}}-T_{\mathrm{yy}})\approx \frac{R_e}{\mathrm{GM}/R_e^3}\sqrt{\underset{l=2}{\overset{L}{\Sigma }}\frac{2l+1}{{[(l+1)(l+2)-2\sqrt{\frac{4{(l+1)}^3(l+2)(2l+3)}{9(2l+1)}}]}^2}{\left(\frac{r}{R_e}\right)}^{2(l+3)}{\sigma }_l^2\left({\mathrm{\delta T}}_{z-x-y}\right)}---\left(13\right)$

如公式（2）所示，在卫星重力梯度的9个张量中，对角张量（垂向分量T_zz和水平分量T_xx,T_yy）是主要分量，非对角张量对地球重力场精度的影响相对于对 角张量基本可忽略。因此，联合公式（10）、（12）和（13），可得基于卫星重力 梯度全张量误差数据δT_xyz反演累计大地水准面精度的误差模型

${\sigma }_N\left(T_{\mathrm{xyz}}\right)\approx \frac{R_e}{\mathrm{GM}/R_e^3}\sqrt{\underset{l=2}{\overset{L}{\Sigma }}\frac{2l+1}{[2{(l+1)}^2{(l+2)}^2+\frac{8{(l+1)}^3(l+2)(2l+3)}{3(2l+1)}\frac{8}{3}\sqrt{\frac{{(l+1)}^5{(l+2)}^3(2l+3)}{(2l+1)}}]}{\left(\frac{r}{R_e}\right)}^{2(l+3)}{\sigma }_l^2\left({\mathrm{\delta T}}_{\mathrm{xyz}}\right)}---\left(14\right)$

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

基于新型卫星重力梯度误差模型法，利用2012年的GOCE-Level-1B卫星重 力梯度误差数据δT_xyz反演累计大地水准面精度的过程如下

第一，首先以0.3°×0.3°为网格分辨率，在地球表面的经度（0°~360°）和纬 度（-90°~90°）范围内绘制网格；其次，按照GOCE卫星轨道在地球表面的轨 迹点位置依次加入δT_xyz；最后，将分布于地球表面的δT_xyz平均归算于划分的网 格点δT_xyz(φ，λ)处。

第二，将δT_xyz(φ，λ)按球谐函数展开为

${\mathrm{\delta T}}_{\mathrm{xyz}}(\phi ,\lambda )=\underset{l=0}{\overset{L}{\Sigma }}\underset{m=0}{\overset{l}{\Sigma }}\left[(C_{{\mathrm{\delta T}}_{\mathrm{lm}}}\cos \mathrm{m\lambda }+S_{{\mathrm{\delta T}}_{\mathrm{lm}}}\sin \mathrm{m\lambda }){\bar{P}}_{\mathrm{lm}}\left(\sin \phi \right)\right]---\left(15\right)$

其中，表示δT_xyz(φ，λ)按球函数展开的系数

$(C_{{\mathrm{\delta T}}_{\mathrm{lm}}},S_{{\mathrm{\delta T}}_{\mathrm{lm}}})=\frac{1}{4\pi }\int \int [{\mathrm{\delta T}}_{\mathrm{xyz}}(\phi ,\lambda )={\bar{Y}}_{\mathrm{lm}}(\phi ,\lambda )\cos \mathrm{\phi d\phi d\lambda }]---\left(16\right)$

δT_xyz在各阶处的方差表示为

${\sigma }_l^2\left({\mathrm{\delta T}}_{\mathrm{xyz}}\right)=\underset{m=0}{\overset{l}{\Sigma }}(C_{{\mathrm{\delta T}}_{\mathrm{lm}}}^2+S_{{\mathrm{\delta T}}_{\mathrm{lm}}}^2)---\left(\text{17}\right)$

将公式（17）代入公式（14），可有效和快速反演地球重力场精度。

图2表示基于新型重力梯度误差模型法反演累计大地水准面精度对比（公 式(14)），虚线表示未加入相关性重力梯度误差δT_z-x-y反演累计大地水准面精度， 实线表示加入相关性重力梯度误差δT_z-x-y（公式(13)）反演累计大地水准面精度。 通过对图2所示结果的对比研究表明：第一，加入相关性重力梯度误差δT_z-x-y反 演累计大地水准面精度（实线）较未加入相关性重力梯度误差δT_z-x-y反演精度 平均提高2～3倍；第二，由于实际上垂向重力梯度T_zz以及水平重力梯度T_xx和T_yy并非相互独立，而是具有较强的相关性，因此，在本发明建立的新型重力梯度 误差模型中加入相关性重力梯度误差δT_z-x-y是进一步提高地球重力场精度的重 要因素；第三，新型重力梯度误差模型法是建立下一代高精度、高空间分辨率 和高阶次地球重力场模型的有效方法。