利用方差-协方差对角张量原理反演地球重力场的方法

摘要

本发明涉及一种精密测量地球重力场的方法，特别是一种利用卫星重力梯度方差-协方差对角张量原理反演地球重力场的方法；基于卫星重力梯度方差-协方差对角张量原理建立累积大地水准面误差模型，通过星载重力梯度仪的卫星重力梯度测量数据精确和快速反演地球重力场，并利用卫星轨道高度和卫星重力梯度仪精度指标开展GOCE-II卫星重力梯度系统需求论证；该方法地球重力场反演精度高，卫星重力梯度反演速度快，卫星观测方程物理含义明确，易于开展卫星重力梯度系统需求分析，对计算机性能要求低。因此，方差-协方差对角张量卫星重力梯度反演法是解算高精度和高空间分辨率地球重力场的有效方法。

说明书

技术领域

本发明涉及卫星重力梯度学、空间大地测量学、地球物理学、宇航学等交 叉技术领域，特别是涉及一种基于卫星重力梯度方差-协方差对角张量原理建立 累积大地水准面误差模型，并利用卫星轨道高度和卫星重力梯度仪精度指标开 展GOCE-II卫星重力梯度系统需求论证研究的方法。

背景技术

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

如图1所示，欧洲空间局（ESA）独立研制的GOCE（Gravity Field and  Steady-State Ocean Circulation Explorer）重力梯度卫星已于2009年3月17日成 功发射升空；采用近圆（轨道离心率0.001）、极地（轨道倾角96.5°）和太阳同 步轨道；经过20个月的飞行计划，轨道高度由250km降为240km；主要用于 精密探测地球重力场的中短波信号。GOCE采用卫星跟踪卫星高低和卫星重力 梯度（SST-HL/SGG）模式的结合，除基于高轨道的GPS/GLONASS卫星对低轨 道的GOCE进行精密跟踪定位（定轨精度1cm），同时利用定位于卫星质心处 的重力梯度仪（测量精度3×10^-12/s²）高精度测量卫星轨道高度处引力位的二阶 导数。GOCE采用了非保守力补偿技术，首先利用重力梯度仪测量由非保守力 （大气阻力、太阳光压、地球辐射压、轨道高度和姿态控制力等）引起的卫星 质心的线性加速度与卫星平台的角加速度；其次，结合卫星平台姿态测量数据， 通过无阻尼离子微推进器补偿卫星受到的非保守力。由于卫星重力梯度观测数 据中的非保守力效应得到了有效扣除，因此进一步提高了重力场反演的精度和 空间分辨率。由于地球重力场信号随卫星轨道高度的增加而呈指数急剧衰减 [R_e/(R_e+H)]^l+1，其中R_e表示地球的平均半径，H表示卫星轨道高度，l表示地 球引力位按球函数展开的阶数。基于分析卫星轨道运动仅适合于确定中长波地 球重力场，而SGG是直接测定地球引力位的二次微分，其结果将球谐系数放大 了l²倍，因此可有效抑制地球引力位随高度的衰减效应，进而高精度感测中高 频地球重力场信号。基于GOCE重力梯度卫星在高精度感测地球中高频重力场 中的优秀表现，以及GOCE卫星重力梯度计划预计于2015年前结束的原因，国 际大地测量、空间科学等交叉研究领域的众多科研机构正积极开展GOCE-II卫 星重力梯度计划的研究论证，旨在进一步提高全球重力场的测量精度以及获得 地球重力场的时变信号。

在众多卫星重力反演方法中，按照卫星观测方程的建立和求解的不同可分 为数值法和解析法。数值法是指将卫星观测数据按时间序列处理，卫星星历值 直接表示成地球引力位系数的函数，由最小二乘法、预处理共轭梯度法等解算 超定方程组，进而获得地球引力位系数。优点是地球重力场求解精度较高；缺 点是不易于误差分析、求解速度较慢、对计算机性能要求较高，难以解算高阶 次地球重力场模型。解析法是指通过分析地球重力场和卫星观测数据的关系建 立卫星观测方程模型，进而估计地球重力场的精度。优点是卫星观测方程物理 含义明确，易于误差分析且可快速求解高阶地球重力场；缺点是在建立卫星观 测方程模型时作了不同程度的近似。不同于已有卫星重力反演法，本发明首次 建立了方差-协方差对角张量卫星重力梯度反演方法。在卫星重力计划可行性研 究的地球重力场需求分析阶段，可通过方差-协方差对角张量卫星重力梯度反演 法有效和快速论证卫星观测模式、卫星轨道参数（轨道高度、轨道倾角、轨道 离心率等）、关键载荷匹配精度指标（重力梯度仪、GPS接收机、非保守力补偿 系统等）等的合理性和最优设计，分析卫星系统各项误差源对地球重力场反演 精度的影响。本发明不仅为卫星重力梯度系统的轨道参数和关键载荷精度指标 的优化设计提供了理论依据和技术支持，同时对国际月球和太阳系火星等行星 卫星重力测量的发展方向具有一定的借鉴意义。

发明内容

为了解决上述技术问题，本发明提供了一种利用卫星重力梯度方差-协方差 对角张量原理反演地球重力场的方法，包含以下步骤：

步骤1：通过重力梯度卫星的星载重力梯度仪采集卫星重力梯度测量数据 δV_xyz；

步骤2：建立方差-协方差对角张量误差模型，使用所述方差-协方差对角张 量误差模型利用卫星重力梯度测量数据δV_xyz反演累积大地水准面误差；其中所 述步骤2包括：

步骤2.1：基于卫星重力梯度测量数据建立卫星重力梯度方差-协方差对 角张量误差模型，包括：

在地固系中按球谐函数展开表达地球引力位V(r,θ，λ)，分别对卫 星位置矢量r的三个分量x,y,z进行二阶求导，将求导结果表达为矩阵形式，得 到转换矩阵A；

基于最小二乘平差法，根据方差-协方差矩阵公式，得到转换矩阵 A的正规方阵逆(A^TA)^-1的主对角元素为地球引力位系数的阶方差；

基于卫星重力梯度对角张量V_xyz建立地球引力位系数的方差估计 模型，分别获得基于卫星重力梯度垂向张量V_zz和水平张量V_xx和V_yy的累积大地 水准面误差模型，将分别基于卫星重力梯度垂向张量V_zz和水平张量V_xx和V_yy的累 积大地水准面误差模型进行联合，得到基于卫星重力梯度对角张量V_xyz的累积大 地水准面联合误差模型，以此作为卫星重力梯度方差-协方差对角张量误差模型；

步骤2.2：基于所述卫星重力梯度方差-协方差对角张量误差模型反演累 积大地水准面误差。

本发明还提供了一种使用利用卫星重力梯度方差-协方差对角张量原理反演 地球重力场的方法来确定GOCE-II重力梯度卫星参数需求的方法，包括如下步 骤：

基于卫星重力梯度方差-协方差对角张量误差模型，利用GOCE卫星在不同 的卫星轨道高度处所获得的测量数据确定卫星轨道高度对地球重力场精度的影 响关系；

基于卫星重力梯度方差-协方差对角张量误差模型，利用GOCE卫星的重力 梯度仪所获得的测量数据确定不同的重力梯度仪精度指标对地球重力场精度的 影响关系；

根据地球科学各相关学科需求确定所需的地球重力场反演精度，通过所述 卫星轨道高度对地球重力场精度的影响关系和所述重力梯度仪精度指标对地球 重力场精度的影响关系，确定GOCE-II重力梯度卫星的关键载荷精度指标和轨 道参数。

本发明提出的卫星重力梯度方差-协方差对角张量法有利于反演高精度和高 空间分辨率地球重力场，其优点是：

1）地球重力场反演精度高；

2）卫星重力梯度反演速度快；

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

4）易于开展卫星重力梯度系统需求分析；

5）对计算机性能要求低。

附图说明

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

图2表示基于不同卫星轨道高度反演GOCE-II地球重力场精度。

图3表示基于不同卫星重力梯度仪精度反演GOCE-II地球重力场精度。

具体实施方式

为了便于本领域普通技术人员理解和实施本发明，下面结合附图及具体实 施方式对本发明作进一步的详细描述。

利用卫星重力梯度方差-协方差对角张量原理反演地球重力场的方法包含下 列步骤：

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

通过重力梯度卫星的星载重力梯度仪采集卫星重力梯度测量数据δV_xyz。

步骤2：建立方差-协方差对角张量误差模型，反演累积大地水准面误差

2.1）基于卫星重力梯度测量数据建立卫星重力梯度方差-协方差对角张量模型

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

$V(r,\theta ,\lambda )=\frac{\mathrm{GM}}{R_e}\underset{l=0}{\overset{L}{\Sigma }}\underset{m=0}{\overset{l}{\Sigma }}{\left(\frac{R_e}{r}\right)}^{l+1}{\bar{Y}}_{\mathrm{lm}}(\theta ,\lambda ){\bar{X}}_{\mathrm{lm}}---\left(1\right)$

其中，表示球函数，${\bar{Y}}_{\mathrm{lm}}(\theta ,\lambda )={\bar{P}}_{l\left|m\right|}\left(\cos \theta \right)Q_m\left(\lambda \right),$$Q_m\left(\lambda \right)=\left(\begin{array}{cc}\cos \mathrm{m\lambda }& m\ge 0\\ \sin \left|m\right|\lambda & m<0\end{array}\right),$GM表示地球质量M和万有引力常数G之积，R_e表示地球平均半径； 表示卫星的地心半径，x,y,z分别表示卫星轨道位置矢量r的三 个分量，θ和λ分别表示卫星的地心余纬度和地心经度，L表示地球引力位按球 谐函数展开的最大阶数；表示规格化Legendre函数，l表示阶数，m表 示次数；表示待求的规格化地球引力位系数。

地球引力位V(r,θ，λ)分别对x,y,z的二阶导数表示为

$\Gamma =\left(\begin{array}{ccc}{V}_{\mathrm{xx}}& V_{\mathrm{xy}}& V_{\mathrm{xz}}\\ V_{\mathrm{yx}}& V_{\mathrm{yy}}& V_{\mathrm{yz}}\\ V_{\mathrm{zx}}& V_{\mathrm{zy}}& V_{\mathrm{zz}}\end{array}\right)---\left(2\right)$

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

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

其中，地球引力位V(r,θ，λ)分别对r,θ，λ的一阶导数表示为

$\left(\begin{array}{c}{V}_r(r,\theta ,\lambda )=-\frac{\mathrm{GM}}{R_e^2}\underset{l=0}{\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)\\ V_{\theta }(r,\theta ,\lambda )=-\frac{\mathrm{GM}}{R_e}\underset{l=0}{\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 \\ V_{\lambda }(r,\theta ,\lambda )=\frac{\mathrm{GM}}{R_e}\underset{l=0}{\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)$

地球引力位V(r,θ，λ)分别对r,θ，λ的二阶导数表示为

$\left(\begin{array}{c}{V}_{\mathrm{rr}}(r,\theta ,\lambda )=\frac{\mathrm{GM}}{R_e^3}\underset{l=0}{\overset{L}{\Sigma }}(l+1)(l+2){\left(\frac{R_e}{r}\right)}^{l+3}\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)\\ V_{\mathrm{\theta \theta }}(r,\theta ,\lambda )=\frac{\mathrm{GM}}{R_e}\underset{l=0}{\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){\sin }^2\theta -{\bar{P}}_{\mathrm{lm}}^{\prime }\left(\cos \theta \right)\cos \theta ]\\ V_{\mathrm{\lambda \lambda }}(r,\theta ,\lambda )=-\frac{\mathrm{GM}}{R_e}\underset{l=0}{\overset{L}{\Sigma }}{\left(\frac{R_e}{r}\right)}^{l+1}\underset{m=0}{\overset{l}{\Sigma }}{m}^2({\bar{C}}_{\mathrm{lm}}\cos \mathrm{m\lambda }+{\bar{S}}_{\mathrm{lm}}\sin \mathrm{m\lambda }){\bar{P}}_{\mathrm{lm}}\left(\cos \theta \right)\\ V_{\mathrm{r\theta }}(r,\theta ,\lambda )=V_{\mathrm{\theta r}}(r,\theta ,\lambda )=\frac{\mathrm{GM}}{R_e^2}\underset{l=0}{\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}}^{\prime }\left(\cos \theta \right)\sin \theta \\ V_{\mathrm{r\lambda }}(r,\theta ,\lambda )=V_{\mathrm{\lambda r}}(r,\theta ,\lambda )=\frac{\mathrm{GM}}{R_e^2}\underset{l=0}{\overset{L}{\Sigma }}(l+1){\left(\frac{R_e}{r}\right)}^{l+2}\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)\\ V_{\mathrm{\theta \lambda }}(r,\theta ,\lambda )=V_{\mathrm{\lambda \theta }}(r,\theta ,\lambda )=\frac{\mathrm{GM}}{R_e}\underset{l=0}{\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}}^{\prime }\left(\cos \theta \right)\sin \theta \end{array}\right)---\left(5\right)$

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}}^{″}\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{1}{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）的矩阵形式表示为

y=A·x      （7）

其中，y表示全张量重力梯度的9个分量，A表示转换矩阵，x表示待求的规格 化地球引力位系数和在公式（7）两边同乘以A^T得

A^Ty＝A^TA·x      （8）

基于最小二乘平差法，根据方差-协方差矩阵公式，正规方阵逆(A^TA)^-1的主 对角元素为地球引力位系数的阶方差

${\sigma }_l^2(\delta {\bar{C}}_{\mathrm{lm}},\delta {\bar{S}}_{\mathrm{lm}}){\sigma }^2\left({\mathrm{\delta V}}_{\mathrm{xyz}}\right)·\mathrm{Diag}{\left(A^{T}A\right)}^{-1}---\left(9\right)$

其中，σ(δV_xyz)表示卫星重力梯度仪的测量精度，表示规格化地球引力 位系数精度。

联合公式（3）和公式（8），基于球函数的正交归一性，垂向重力梯度V_zz正 规方阵的对角元素表示为

$\mathrm{Diag}{\left(A^{T}A\right)}_{V_{\mathrm{zz}}}={\left(\frac{\mathrm{GM}}{R_e^3}\right)}^2{\left(\frac{R_e}{R_e+H}\right)}^{2l+6}{(l+1)}^2{(l+2)}^2---\left(10\right)$

其中，H表示卫星的平均轨道高度。

水平重力梯度V_xx和V_yy正规方阵和的对角元素表示为

$\mathrm{Diag}{\left(A^{T}A\right)}_{V_{\mathrm{xx}}}\approx \mathrm{Diag}{\left(A^{T}A\right)}_{V_{\mathrm{yy}}}={\left(\frac{\mathrm{GM}}{R_e^3}\right)}^2{\left(\frac{R_e}{R_e+H}\right)}^{2l+6}\frac{{(l+1)}^3(l+2)(2l+3)}{9(2l+1)}---\left(11\right)$

据误差原理可知，如果在卫星轨道采样点处重力梯度观测数据有K₀个，那 么地球引力位系数方差正比于1/K₀

K₀=T/Δt      （12）

其中，T表示卫星观测时间，Δt表示卫星观测数据的采样间隔。

在卫星重力梯度的9个张量中，对角张量垂向梯度V_zz和水平梯度V_xx,V_yy是 主要分量，其它非对角张量对地球重力场精度的影响较对角张量可忽略。因此， 联合公式（9）~（11），基于卫星重力梯度对角张量V_xyz估计地球引力位系数的 方差表示为

${\sigma }_l^2{(\delta {\bar{C}}_{\mathrm{lm}},\delta {\bar{S}}_{\mathrm{lm}})}_{V_{\mathrm{xyz}}}=\frac{{\sigma }^2\left({\mathrm{\delta V}}_{\mathrm{xyz}}\right)}{K_0[\mathrm{Diag}{\left(A^{T}A\right)}_{V_{\mathrm{xx}}}+\mathrm{Diag}{\left(A^{T}A\right)}_{V_{\mathrm{yy}}}+\mathrm{Diag}{\left(A^{T}A\right)}_{V_{\mathrm{zz}}}]}---\left(13\right)$

大地水准面方差表示为

${\sigma }_l^2\left(N\right)=R_e^2\underset{l=0}{\overset{L}{\Sigma }}\underset{m=-l}{\overset{l}{\Sigma }}[{\left(\delta {\bar{C}}_{\mathrm{lm}}\right)}^2+{\left(\delta {\bar{S}}_{\mathrm{lm}}\right)}^2]---\left(14\right)$

联合公式（10）~（14），基于卫星重力梯度垂向张量V_zz，累积大地水准面 误差模型表示为

${\sigma }_l{\left(N\right)}_{V_{\mathrm{zz}}}=\frac{R_e^4}{\mathrm{GM}\sqrt{T/\mathrm{\Delta t}}}\sqrt{\underset{l=0}{\overset{L}{\Sigma }}\frac{2l+1}{{\left(\frac{R_e}{R_e+H}\right)}^{2l+6}{(l+1)}^2{(l+2)}^2}{\sigma }^2\left({\mathrm{\delta V}}_{\mathrm{zz}}\right)}---\left(15\right)$

基于卫星重力梯度水平张量V_xx和V_yy，累积大地水准面误差模型表示为

${\sigma }_l{\left(N\right)}_{V_{\mathrm{xx}}}\approx {\sigma }_l{\left(N\right)}_{V_{\mathrm{yy}}}=\frac{R_e^4}{\mathrm{GM}\sqrt{T/\mathrm{\Delta t}}}\sqrt{\underset{l=0}{\overset{L}{\Sigma }}\frac{2l+1}{{\left(\frac{R_e}{R_e+H}\right)}^{2l+6}\frac{{(l+1)}^3(l+2)(2l+3)}{9(2l+1)}}{\sigma }^2\left({\mathrm{\delta V}}_{\mathrm{xx}}\right)}---\left(16\right)$

基于卫星重力梯度对角张量V_xyz，累积大地水准面联合误差模型表示为

${\sigma }_l{\left(N\right)}_{V_{\mathrm{xyz}}}=\frac{R_e^4}{\mathrm{GM}\sqrt{T/\mathrm{\Delta t}}}\sqrt{\underset{l=0}{\overset{L}{\Sigma }}\frac{2l+1}{{\left(\frac{R_e}{R_e+H}\right)}^{2l+6}[{(l+1)}^2{(l+2)}^2+\frac{2{(l+1)}^3(l+2)(2l+3)}{9(2l+1)}]}{\sigma }^2\left(\delta V_{\mathrm{xyz}}\right)}---\left(17\right)$

基于公式（15）和公式（16），可以分别衡量卫星重力梯度垂向张量V_zz和水平张 量V_xx和V_yy三者分别对累积大地水准面误差的影响；公式（17）为综合公式（15） 和公式（16）所得到的累积大地水准面联合误差模型，综合反映了卫星重力梯 度垂向张量V_zz和水平张量V_xx和V_yy三者整体对累积大地水准面误差的影响。

2.2）基于卫星重力梯度方差-协方差对角张量误差模型反演累积大地水准面误差

基于卫星重力梯度方差-协方差对角张量法，利用卫星重力梯度测量数据 δV_xyz反演累积大地水准面误差的过程如下（其中卫星重力梯度测量数据可以采 用GOCE卫星的重力梯度测量数据或者是未来GOCE-II卫星的重力梯度测量数 据，下面以GOCE卫星重力梯度测量数据为例进行说明）：

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

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

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

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

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

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

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

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

步骤3：GOCE-II卫星重力梯度系统的需求论证

3.1）卫星轨道高度影响

图2表示基于卫星重力梯度方差-协方差对角张量法，利用不同的卫星轨道 高度（200km~500km）估计GOCE-II地球重力场精度。研究结果表明：重力 梯度卫星GOCE-II的轨道高度设计为300~400km较优，原因分析如下：

第一，据欧空间局公布的GOCE-Level-1B中GPS导航实测数据可知，GOCE 卫星的轨道高度主要分布在距地面200~300km的空间范围。经过2年多的全 球重力场测量，GOCE卫星已高精度和高空间分辨率地感测了地球中短波静态 重力场。由于不同卫星轨道高度敏感于不同频段的地球重力场信号，因此GOCE 卫星仅能在特定轨道高度区间发挥其优越性，而在轨道覆盖空间范围之外基本 无能为力。如果重力梯度卫星的轨道高度也同样设计在200~300km的空间范 围，除非反演地球重力场的精度高于GOCE卫星，否则其效果仅相当于GOCE 卫星的简单重复测量，对于地球重力场精度的进一步提高没有实质性贡献。因 此，重力梯度卫星的轨道高度应尽可能选择在GOCE的测量盲区，进而与GOCE 形成互补的态势。

第二，利用重力卫星作为传感器进行地球重力场测量的最大弱点是卫星轨 道高度处的地球重力场成指数衰减。随着重力卫星轨道逐步升高，地球重力场 的长波信号衰减幅度较小，中波信号衰减幅度次之，短波信号衰减幅度最大。 因此，较高轨道的重力卫星对地球重力场中波和短波信号的敏感性较弱，不利 于高阶地球重力场反演。为了克服上述缺点进而反演高精度、高空间分辨率和 全频段的地球重力场，目前最有效的办法是适当降低卫星轨道高度。GOCE卫 星为了尽可能抑制地球重力场信号随卫星轨道高度的衰减效应，因此采用了极 低轨设计（平均轨道高度250km），虽然理论上可以提高地球重力场反演的精度 和空间分辨率，但其实际负面效应不容忽视：（1）卫星轨道高度每降低100km， 大气阻力将提高约10倍，为调整卫星轨道高度和姿态需频繁进行轨道机动，不 稳定的卫星平台工作环境将影响关键载荷的测量精度；（2）由于卫星频繁喷气 引起喷气燃料消耗，将导致星体质心和加速度计检验质量质心存在实时偏差； （3）卫星使用寿命极大地缩减，将影响静态和时变地球重力场反演的精度和空 间分辨率。因此，合理选择卫星轨道高度是反演高精度和高空间分辨率地球重 力场的重要保证。为了有效弥补GOCE卫星设计的不足，GOCE-II重力梯度卫 星可将轨道高度设计在300~400km区间，不仅可有效提高卫星的使用寿命进 而获得地球重力场时变信号，而且通过适当升高轨道高度可为卫星关键载荷（冷 原子干涉重力梯度仪、GNSS复合接收机、非保守力补偿系统、以及恒星敏感器 等）的高精度测量提供安静和稳定的卫星平台工作环境。

3.2）卫星重力梯度仪精度影响

图3表示基于卫星重力梯度方差-协方差对角张量法，利用不同的重力梯度 仪精度指标（10^-11/s²~10^-15/s²）反演GOCE-II地球重力场精度。重力梯度卫星 GOCE-II的重力梯度仪精度指标可设计为10^-13/s²~10^-15/s²，具体原因分析如下： 卫星重力梯度仪是一种能直接探测空间重力加速度矢量梯度的传感器。由于重 力梯度可以较好地反映等位面的曲率和力线的弯曲程度，因此敏感于中短波地 球重力场的信号，更能反应地球重力场的精细结构。在地球卫星内的微重力环 境中，由于不同位置点加速度的差异较小，因此不同属性的重力梯度仪通常由1 ~3对属性相同的加速度计按不同的排列方式组合而成，精确测定每对加速度计 检验质量之间的相对位置变化，通过观测重力加速度差进而得到重力梯度张量， 此为卫星重力梯度仪能在微重力环境下直接测量地球重力场参数的主要原因。 目前卫星重力梯度仪主要包括旋转式重力梯度仪、静电悬浮重力梯度仪（GOCE 卫星）、超导重力梯度仪、冷原子干涉重力梯度仪等。国际卫星重力梯度测量工 程的发展方向以采用高精度和新型的冷原子干涉重力梯度仪为主流。由于冷原 子干涉重力梯度仪具有灵敏度高、结构简单、成本低、抗外界干扰能力强、易 于自动化数据采集等优点，同时我国已具有一定的研究基础，因此，GOCE-II 卫星重力梯度计划搭载冷原子干涉重力梯度仪（10^-13/s²～10^-15/s²）较优。

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