利用不同轨道倾角卫星联合反演地球重力场的方法

摘要

本发明涉及一种地球重力场精密测量方法，特别是一种利用不同轨道倾角卫星联合反演地球重力场的方法；基于扰动星间加速度法，利用GRACE卫星K波段测距仪的星间距离、星间速度和星间加速度测量数据、GPS接收机的卫星轨道位置测量数据和卫星轨道速度测量数据、以及加速度计的非保守力测量数据，通过不同轨道倾角卫星最优组合，进而精确和快速反演地球重力场；该方法卫星重力反演精度高，同时敏感于引力位带谐项和田谐项系数精度，地球重力场解算速度快，卫星观测方程物理含义明确，计算机性能要求低；因此，不同轨道倾角卫星联合法是反演高精度和高空间分辨率地球重力场的优选方法。

说明书

技术领域

本发明涉及卫星大地测量学、卫星重力学、空间科学等交叉技术领域，特 别是涉及一种基于扰动星间加速度法，利用GRACE卫星K波段测距仪的星间 距离、星间速度和星间加速度测量数据、GPS接收机的卫星轨道位置测量数据 和卫星轨道速度测量数据、以及加速度计的卫星非保守力测量数据，通过不同 轨道倾角卫星最优组合，进而精确和快速反演地球重力场的方法。

背景技术

自伽利略于16世纪末第一次进行地球重力场测量以来，国内外的许多科研 机构在全球范围内的陆地、海洋和空间采用多种技术和方法进行了大量的地球 重力场测量。目前卫星重力测量技术的实现是继美国全球定位系统（GPS）星座 成功构建之后在大地测量领域的又一项创新和突破。它不仅促进了大地测量学 本身的发展，同时带动了地球物理学、大气研究、海洋探测、空间技术等领域 的发展。卫星重力测量技术之所以被国际大地测量学界公认为是当前地球重力 场探测研究中最高效、最经济和最有发展潜力的方法之一，是因为它既不同于 传统的车载、船载和机载测量，也不同于轨道摄动分析和卫星测高技术，而是 通过卫－卫跟踪技术（SST）和卫星重力梯度技术（SGG）反演高精度和高空间 分辨率地球重力场。

如图1所示，美国宇航局（NASA）和德国波兹坦地学研究中心（GFZ）合 作研制并于2002年3月17日成功发射的GRACE双星是继2000年7月15日升 空的CHAMP单星之后的又一项专用于地球重力场中长波精密测量的卫星计划。 GRACE编队飞行计划是卫星重力测量中卫－卫跟踪技术（SST）的再次成功实 现，在10年的飞行任务中采用变轨技术（500km~300km）感测不同波段的地 球重力场。为了满足GRACE整体卫星系统测量地球重力场的要求，双星轨道设 计为近极轨模式（轨道倾角89°）。对于反演120阶GRACE地球重力场而言， 由于89°轨道倾角在地球南北极形成的极沟区（未覆盖区）2×|90°-I₁|=2°小于 对应的空间分辨率360°/L_max＝3°，因此，该模式的优点是不仅可达到卫星近似全 球覆盖的目的，同时可忽略极沟区对地球重力场反演精度的影响。

由于适当增大卫星的轨道倾角有利于提高地球引力位带谐项系数反演的精 度，适当降低卫星的轨道倾角有利于提高地球引力位田谐项系数反演的精度， 因此，采用多颗不同轨道倾角卫星联合测量可互相取长补短，进而共同反演高 精度和高空间分辨率的地球重力场。至目前为止，国内外众多学者在卫星不同 轨道倾角如何影响地球引力位系数反演精度的方面开展了卓有成效的研究工 作。GRACE卫星采用高轨道倾角89°的设计可有效提高地球引力位带谐项系数 的精度，但对地球引力位田谐项系数的敏感度较低，因此可采用第二组较低轨 道倾角的卫星高精度测量地球引力位田谐项系数。数值模拟结果表明：第二组 GRACE卫星采用轨道倾角82°~84°是较优选择，该设计可有效弥补单组89°轨 道倾角卫星对地球引力位田谐项系数敏感度较低的不足。综上所述，不同轨道 倾角卫星的联合测量是反演高精度和高空间分辨率地球重力场的有效途径。

本发明基于扰动星间加速度法利用不同轨道倾角卫星的组合反演了120阶 GRACE地球重力场，阐述了卫星不同轨道倾角和地球引力位带谐项、扇谐项和 田谐项系数反演精度的关系，得到了两组GRACE双星分别采用89°和83°轨道 倾角联合反演累计大地水准面的精度较单组89°轨道倾角的精度平均提高约2倍 的结论，分析了两组GRACE双星分别采用89°和82°~84°轨道倾角反演地球重 力场是较优组合的原因。

发明内容

本发明的目的是：基于扰动星间加速度法，通过不同轨道倾角卫星联合进 一步提高地球重力场反演精度。

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

一种利用不同轨道倾角卫星联合反演地球重力场的方法，包括如下步骤：

步骤1：采集GRACE双星的关键载荷数据从而获得双星测量数据，通过 GRACE双星的星载K波段测距仪采集星间距离ρ₁₂、星间速度和星间加速度 通过星载GPS接收机采集双星轨道位置（r₁,r₂）和双星轨道速度通过星载加速度计获取作用于双星的非保守力（f₁，f₂）；

步骤2：通过残余地心引力差、残余地球扰动引力差、残余保守力差和残余 非保守力差对扰动轨道加速度差的影响，建立扰动星间加速度的观测方程，将 步骤1中所获得的双星测量数据代入扰动星间加速度观测方程，使用最小二乘 法反演获得地球引力位系数和

步骤3：基于扰动星间加速度观测方程，利用不同轨道倾角卫星的测量数据 进行地球重力场反演，将反演得到的地球引力位系数精度进行对比，进而分析 不同轨道倾角的GRACE双星观测数据分别对地球引力位带谐项、扇谐项、田谐 项系数精度的影响；

步骤4：基于不同轨道倾角的双星测量数据联合反演大地水准面精度，其中 不同轨道倾角的双星测量数据至少包括一组89°高轨道倾角的GRACE双星测量 数据和另一组82°~84°低轨道倾角的GRACE双星测量数据。

本发明是利用不同轨道倾角卫星联合有利于精确反演地球重力场的特点而 设计的，优点是：

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

2）同时敏感于引力位带谐项和田谐项系数精度；

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

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

5）计算机性能要求低。

附图说明

图1为GRACE双星测量原理图。

图2a-2d表示基于不同轨道倾角卫星反演地球引力位系数精度对比；

图2a为基于不同轨道倾角卫星反演地球引力位带谐项系数精度对比；

图2b为基于不同轨道倾角卫星反演地球引力位扇谐项系数精度对比；

图2c为基于不同轨道倾角卫星反演地球引力位田谐项系数精度对比

图2d为基于不同轨道倾角卫星反演地球综合引力位系数精度对比。

图3表示基于89°轨道倾角和不同轨道倾角组合反演引力位系数阶误差之比 的平均值。

图4表示基于两组不同轨道倾角卫星联合反演累计大地水准面精度。

具体实施方式

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

利用不同轨道倾角卫星联合反演地球重力场方法的应用：

步骤1：卫星关键载荷数据采集

（1）通过星载K波段测距仪采集一组GRACE-A/B双星的星间距离ρ₁₂、 星间速度和星间加速度通过星载GPS接收机采集双星轨道位置（r₁,r₂） 和双星轨道速度通过星载加速度计获取作用于双星的非保守力（f₁,f₂）。

（2）利用9阶Runge-Kutta线性单步法和12阶Adams-Cowell线性多步法 数值模拟公式获取双星参考轨道位置和双星参考轨道速度

（3）参考星间距离参考星间速度和参考星间加速度通过参考轨 道位置和参考轨道速度计算获得。

（4）参考非保守力通过DTM2000阻力温度模型计算获得。

（5）通过国际公布模型DE-405、IERS96和CSR4.0联合计算获取作用于 双星的保守力（F₁,F₂）和参考保守力

其中步骤（4）、（5）的计算方法已在【郑伟,许厚泽,钟敏,员美娟,周旭 华,彭碧波.卫星跟踪卫星测量模式中星载加速度计高低灵敏轴分辨率指标优 化设计论证.地球物理学报,2009,52(11):2712-2720.】和【Tapley B,Ries J, Bettadpur S,Chambers D,Cheng M,Condi F,Gunter B,Kang Z,Nagel P，Pastor R, Pekker T,Poole S,Wang F.GGM02-An improved Earth gravity field model from  GRACE.Journal of Geodesy,2005,79(8):467–478.】中公开。

步骤2：扰动星间加速度观测方程建立

在地心惯性系中，GRACE-A/B的星间距离ρ₁₂表示如下

ρ₁₂=r₁₂·e₁₂         （1） 其中，r₁₂=r₂-r₁表示GRACE-A/B的相对轨道位置矢量，r₁和r₂分别表示双星的 绝对轨道位置矢量；e₁₂=r₁₂/|r₁₂|表示由GRACE-A指向GRACE-B的单位矢量。

GRACE-A/B的参考星间距离表示如下

其中，表示GRACE-A/B的相对参考轨道位置矢量，和分别表示 双星的绝对参考轨道位置矢量。

通过公式（1）-公式（2），GRACE-A/B的扰动星间距离δρ₁₂表示如下

δρ₁₂=δr₁₂·e₁₂         （3）

其中，

在（3）式两边同时对时间t求导数，可得GRACE-A/B的扰动星间速度

$\delta {\stackrel{·}{\rho }}_{12}=\delta {\stackrel{·}{r}}_{12}·e_{12}+\delta r_{12}·{\stackrel{·}{e}}_{12}---\left(4\right)$

其中，表示GRACE-A/B的相对扰动轨道速度矢量，和表示相 对轨道速度和相对参考轨道速度矢量，表示垂直于GRACE-A/B连线的单位 矢量

${\stackrel{·}{e}}_{12}=\frac{\delta {\stackrel{·}{r}}_{12}-\delta {\stackrel{·}{\rho }}_{12}{e}_{12}}{\delta {\rho }_{12}}---\left(5\right)$

因为所以（4）式可简化为

$\delta {\stackrel{·}{\rho }}_{12}=\delta {\stackrel{·}{r}}_{12}·e_{12}---\left(6\right)$

在（6）式两边同时对时间t求导数，可得GRACE-A/B的扰动星间加速度

$\delta {\stackrel{··}{\rho }}_{12}=\delta {\stackrel{··}{r}}_{12}·e_{12}+\delta {\stackrel{·}{r}}_{12}·{\stackrel{·}{e}}_{12}---\left(7\right)$

其中，表示GRACE-A/B的相对扰动轨道加速度矢量

$\delta {\stackrel{··}{r}}_{12}=\delta g_{12}+\delta T_{12}+\delta F_{12}+\delta f_{12}---\left(8\right)$

其中，δT₁₂表示作用于双星的残余地球扰动引力差；表示除地球引力之外的残余保守力差，F₁和F₂表示作用于双星的保守力，和 表示参考保守力；表示残余非保守力差，f₁和f₂表示作用于双星的非保守力，和表示参考非保守力； 表示残余地心引力差，g₁和g₂表示双星的地心引力， 和表示参考地心引力

其中，GM表示地球质量M和万有引力常数G的乘积，分别表示双星的地心半径，x₁₍₂₎,y₁₍₂₎,z₁₍₂₎表示轨道位置矢量r₁₍₂₎的3个分量。

基于公式（5）～公式（8），扰动星间加速度观测方程表示如下

其中，表示残余地球扰动位的一阶梯度，V₁和V₂表 示地球扰动位，和表示参考地球扰动位

$V(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)$

其中，r,θ，λ分别表示地心半径、地心余纬度和地心经度，R_e表示地球平均半径， L表示地球引力位按球谐函数展开的最大阶数；表示正规化的缔合 Legendre函数，l表示阶数，m表示次数；和表示待估的地球引力位系数。

最后，通过将“步骤1：卫星关键载荷数据采集”中的双星测量数据代入扰动 星间加速度观测方程（10），利用最小二乘法反演获得地球引力位系数和

步骤3：分析不同轨道倾角卫星对地球引力位系数精度的影响

图2a~2d表示基于扰动星间加速度法，采用GRACE卫星公布的关键载荷精 度指标（星间速度1×10^-6m/s、轨道位置3×10^-2m、轨道速度3×10^-5m/s、非保守 力3×10^-10m/s²），利用不同轨道倾角（85°、87°和89°）卫星的测量数据反演地 球引力位系数精度对比。图2a表示反演引力位带谐项系数（l≠0,m=0）的精 度；图2b表示反演引力位扇谐项系数（l＝m≠0）的精度；图2c表示反演引力 位田谐项系数（l≠m≠0）的精度；图2d表示反演引力位系数（带谐、扇谐和 田谐项综合贡献）的精度。星号线表示德国地学研究中心（GFZ）公布的120 阶EIGEN-GRACE02S地球重力场模型的引力位系数的实测精度；实线表示基于 89°轨道倾角反演引力位系数的精度；圆圈线表示基于87°轨道倾角反演引力位 系数的精度；虚线表示基于85°轨道倾角反演引力位系数的精度。

经图2a~2d对比，模拟结果表明：

第一，据图2a可知，随着轨道倾角逐渐增加（85°、87°和89°），反演地球 引力位带谐项系数的精度依次提高。原因分析如下：地球引力位带谐项系数反 演的精度决定于反演地球重力场空间分辨率和地球两极的极沟尺寸的比值 比值越大反演引力位带谐项系数的精度越高。Mackenzie基于Kaula 提出的轨道倾角函数阐述了随着单颗卫星轨道倾角逐渐增加（0°~ 89°），地球引力位带谐项系数对应的轨道倾角函数的幅值依次增大，进而地球引 力位带谐项系数对应的地球引力位信号依次增强，因此，反演地球引力位带谐 项系数精度依次提高。综上所述，适当增大卫星的轨道倾角有利于提高地球引 力位带谐项系数反演的精度。

第二，据图2b可知，随着轨道倾角逐渐增加（85°、87°和89°），反演地球 引力位扇谐项系数的精度无显著变化。原因分析如下：地球引力位扇谐项系数 反演的精度决定于卫星重力测量数据的空间分辨率，D＝20000/L_max。对于反演 120阶GRACE地球重力场而言，由于三种不同轨道倾角设计方案中卫星重力测 量数据的空间分辨率均相同，因此反演地球引力位扇谐项系数的精度基本相同。 为保证反演地球引力位扇谐项系数的精度，卫星绕地球飞行的总圈数N_r应至少 大于反演地球重力场最高阶数L_max的2倍（N_r>2L_max）。综上所述，卫星轨道倾 角的变化对反演地球引力位扇谐项系数精度的影响较小。

第三，据图2c可知，随着轨道倾角逐渐增加（85°、87°和89°），反演地球 引力位田谐项系数的精度依次降低。原因分析如下：地球引力位田谐项系数的 精度决定于卫星轨道在地球表面覆盖面积内测量数据的密度。当卫星测量数据 长度和采样间隔均相同时，如果卫星轨道倾角较小将导致轨道覆盖面积内测量 数据的密度较大，因此反演地球引力位田谐项系数的精度也较高。综上所述， 适当降低卫星的轨道倾角有利于提高地球引力位田谐项系数反演的精度。

第四，据图2d可知，随着轨道倾角逐渐增加（85°、87°和89°），反演引力 位系数（带谐、扇谐和田谐项综合贡献）的精度整体呈升高趋势。但是，由于 不同卫星轨道倾角敏感于不同阶l和次m的引力位系数，因此采用多颗不同轨 道倾角卫星联合测量可互相取长补短，进而共同反演高精度和高空间分辨率地 球重力场。

步骤4：基于不同轨道倾角卫星联合反演大地水准面精度

两组GRACE双星分别采用89°和0°~89°轨道倾角运行，根据两组GRACE 双星测量数据联合反演120阶GRACE地球重力场的引力位系数表示如下

$P_{\mathrm{lm}}^1{x}_{\mathrm{lm}}^{89}+P_{\mathrm{lm}}^2{x}_{\mathrm{lm}}^i=x_{\mathrm{lm}}^{89+i}---\left(11\right)$

其中，和分别表示基于89°轨道倾角和0°~89°轨道倾角反演的地球引力 位系数向量；和分别表示地球引力位系数和的 权系数向量，和分别表示基于89°轨道倾角和0°~89°轨道倾角反演地球 引力位系数的方差向量，可由各自的最小二乘协方差阵的对角线元素得到

$D\left(x\right)={\sigma }^2{\left(A_{t\times n}^T{A}_{t\times n}\right)}^{-1}---\left(12\right)$

其中，表示卫星测量数据单位权方差的无偏估值，表 示卫星观测方程的残差向量，G_T×1表示卫星测量数据向量，A_t×n表示t行n列的 测量数据和待求引力位系数之间的转换矩阵，t表示卫星轨道测量数据的数量， 表示待求引力位系数的个数，x_n×1表示待求的引力位系数向 量。正规矩阵是一个块对角占优结构的稠密阵，对于反演120阶 GRACE地球重力场而言，正规矩阵N_n×n为14637行×14637列的超大型方阵， 因此直接求逆较为困难。基于N_n×n的块对角占优性，本发明选取N_n×n的块对角 部分作为预处理阵M_n×n，形成的M_n×n阵为主对角线上按次数m排列，其余部分 为0的块对角方阵，如此选取不仅保留了N_n×n阵的主要特征，而且较易于计算。因此，（12）式可表示为

$D\left(x\right)={\sigma }^2{M}_{n\times m}^{-1}---\left(13\right)$

图3表示基于89°轨道倾角和不同轨道倾角组合89°和0°~89°反演引力位 系数阶误差之比的平均值，其中，表示基于89° 轨道倾角反演地球引力位系数的阶误差向量，表示基于89°和0°~89° 轨道倾角组合反演引力位系数的阶误差向量。图4表示在各阶处基于不同轨道 倾角反演地球累计大地水准面精度对比。星号线表示德国地学研究中心（GFZ） 公布的EIGEN-GRACE02S地球重力场模型的累计大地水准面的实测精度。采用 GRACE卫星公布的关键载荷精度指标，虚线表示基于89°轨道倾角反演累计大 地水准面的精度；实线表示两组GRACE双星分别采用89°和83°轨道倾角联合 反演累计大地水准面的精度。

据图3和图4中数值模拟结果可知：两组GRACE双星分别采用89°和82°~ 84°轨道倾角反演120阶GRACE地球重力场是较优组合。在120阶内，两组 GRACE双星分别采用89°和83°轨道倾角联合反演累计大地水准面的精度较单 组89°轨道倾角的精度平均提高约2倍。原因分析如下：由于不同卫星轨道倾角 敏感于不同阶l和次m的引力位系数，因此采用多颗不同轨道倾角卫星联合测 量可互相取长补短共同反演高精度和高空间分辨率的地球重力场。

本发明所提出的基于不同轨道倾角卫星联合反演大地水准面精度的方法， 可以采用多组GRACE双星分别在不同轨道倾角的轨道上运行以获得不同轨道 倾角的卫星测量数据，优选采用两组GRACE双星分别在高低不同的轨道倾角的 轨道上运行以获得不同轨道倾角的卫星测量数据，也可以采用一组GRACE双星 分别在不同轨道倾角的轨道上运行获得不同轨道倾角卫星测量数据。目前现有 的GRACE系统采用双星工作方式运行于89°轨道倾角，即处于高轨道倾角位置， 可以发射另一组GRACE双星运行于低轨道倾角，从而实现基于不同轨道倾角卫 星的联合反演。

本发明第一组GRACE双星将轨道倾角设计为89°，对于反演120阶地球重 力场而言，89°轨道倾角在地球南北极形成的极沟区2×|90°-I₁|=2°小于对应的 空间分辨率360°/L_max＝3°，因此，该模式不仅可达到卫星近似全球覆盖的目的， 同时可保证120阶地球重力场反演的精度。第一组GRACE双星采用高轨道倾角 89°的设计可有效提高地球引力位带谐项系数的精度，但对地球引力位田谐项系 数的敏感度较低，因此需要第二组较低轨道倾角的卫星高精度测量地球引力位 田谐项系数，以弥补单组89°轨道倾角卫星反演120阶GRACE地球重力场的不 足。如图3所示，本发明通过权系数和的合理引入，基于89°轨道倾角和 0°~89°轨道倾角反演的每组引力位系数均按误差最小原则进行优化组合，因此， 有效提高了地球重力场反演的精度。如果第二组GRACE卫星的轨道倾角设计较 高85°~89°，其结果只相当于高轨道倾角卫星测量数据的简单重复，因此不能 有效提高地球引力位田谐项系数的精度；如果第二组GRACE卫星的轨道倾角设 计较低0°~81°，虽然适当提高了地球引力位田谐项系数的精度，但同时也急剧 降低了地球引力位带谐项系数的精度，其结果使地球引力位带谐项系数精度降 低的幅度超过了地球引力位田谐项系数精度提高的幅度，因此不能有效提高地 球重力场反演的精度；第二组GRACE卫星采用轨道倾角82°~84°是较优选择， 82°~84°轨道倾角的设计可有效弥补单组89°轨道倾角卫星对地球引力位田谐项 系数敏感度较低的不足。综上所述，两组GRACE双星分别采用89°和82°~84° 轨道倾角反演120阶GRACE地球重力场是较优组合。

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