Patent CNB - 基于双星空间三维插值原理的卫星重力反演方法 - Google PatentsCN BGrantCN Dec 26, 2012Jun 3, 2011Jun 3, 2011.X, CN
B, CN B, CN , CN-B-, CN B, CNB, CN, CN.X, , , , , ,
(4) , 基于双星空间三维插值原理的卫星重力反演方法
本发明涉及一种精密测量地球重力场的方法，特别是一种基于双星空间三维插值原理的卫星重力反演方法；基于将位于空间位置相对不规则的重力双星轨道上的地球扰动位(由GRACE卫星K波段测量仪的星间速度、GPS接收机的卫星轨道位置和卫星轨道速度、加速度计的卫星非保守力、以及恒星敏感器的卫星三维姿态数据计算获得)三维插值到空间位置相对规则且均匀网格划分的参考球面上，进而精确和快速反演了地球重力场；该方法卫星重力反演精度高，解算过程物理含义明确，计算速度快，计算机性能要求低，敏感于重力场中高频信号；因此，双星空间三维插值法是解算高精度和高空间分辨率地球重力场的有效方法。
1. 一种基于双星空间三维插值原理的卫星重力反演方法，其特征如下：
步骤一：利用GRACE卫星K波段测量仪获取星间速度、利用GPS接收机获取卫星轨道位置和速度、利用加速度计获取卫星非保守力、以及利用恒星敏感器获取卫星三维姿态，通过上述数据测量计算卫星轨道处的地球扰动位T (X，y, ζ)；
步骤二 ：以地心为球心选择4个等间距且规则的参考球面A, r2，r3和r4，距地心最近参考球面的半径为A = 6830km，球面间隔为Ar = 20km，轨道高度H = 500km，卫星轨道位于!■2和1*4参考球面之间，在每个球面上按照分辨率=纬度Λ θ X经度Λ λ进行均匀网格划分，其中 Λ Θ 取为 180° /2048，Δ λ 取为 360° /4096 ；
步骤三：利用步骤一中的卫星轨道上测量的每个地球扰动位T(x，y, ζ)三维插值得到4个参考球面上对应的单元格子的地球扰动位T Cr, θ , λ );单元格子的分辨率设计为Δ θ X Δ λ X Ar，其中 Λ Θ 取为 180。/2048，Δ λ 取为 360。/4096，Ar 取为 20km ；步骤四：将纬度Θ划分为2048等份，并固定每个规则参考球面上有限的纬度网格划分值；其次，计算每个小区间中点处的Legendre函数值并存储,其它点上的Legendre函数值l(COS^)用已知点上的函数值作泰勒展开逼近获得，并将计算结果存储以备后面调用；步骤五：利用空间三维插值到规则的4个参考球面上的地球扰动位T (r，θ，λ)解算大型超定方程组（I)，最终获得地球引力位系数；
其中，GM表示地球质量M和万有引力常数G之积；R6表示地球的平均半径；r =如2+y2+ Z2表示卫星的地心半径，X, y, Z分别表示卫星位置矢量r的三个分量;Θ和λ分别表示卫星的地心余纬度和地心经度；@OT(C0S^)表示规格化的Legendre函数，I表示阶数，m表示次数山表示地球引力位按球谐函数展开的最大阶数和表示待求的规格化地球引力位系数；
所述步骤三中的插值公式为：
其中，（X，y，Z) ∈ [X2，X3] X [y2，Υ3] X [z2 , z3]，
T(x, y, z)表示卫星轨道上的地球扰动位,（气，表示参考球面上的直角坐标，Wab表示将卫星轨道上的地球扰动位三维插值到规则参考球面上的权重系数，其中a = 1,2,3,4,
，b = χ, y, ζ表示卫星轨道上的三维直角坐标，b2 = x2, y2, Z2表示第二个参考 球面上的三维直角坐标，Aqn表示卫星轨道和第二个参考球面间的距离；将双星空间三维插值公式的权重系数以矩阵表示，则双星空间三维插值公式矩阵形式表示如下:
因此可得经过双星空间三维插值后4个参考球面上的扰动位
所述步骤五中解算大型超定方程组为：
将（I)式以矩阵形式表示为
DtXl — AtXn · XnX I (5)
其中，Dtxi表示基于双星空间三维插值法将卫星轨道处测量的地球扰动位T(X，y, ζ)空间三维插值到参考球面后的列矩阵，t表示观测点的数量；\&&1表示待求地球引力位系数向量(Gm,I)，其中引力位系数按照阶数I排列形成，次数m固定，W = Z^ax +21_ -3表示待求引力位系数的个数，Lfflax表示引力位按球函数展开的最大阶数;AtXn表示t行η列的设计矩阵，所述设计矩阵表示在地球扰动位卫星观测方程（I)右边除了地球引力位系数(G,,,I)外的其它部分；在上式两边同乘C获得待求的地球引力位系数
基于双星空间三维插值原理的卫星重力反演方法
一、技术领域
[0001] 本发明涉及卫星大地测量学、地球物理学、空间科学等交叉技术领域，特别是涉及一种基于将位于空间位置相对不规则的重力双星轨道上的地球扰动位（由GRACE卫星K波段测量仪的星间速度、GPS接收机的卫星轨道位置和卫星轨道速度、加速度计的卫星非保守力、以及恒星敏感器的卫星三维姿态数据计算获得）三维插值到空间位置相对规则且均匀网格划分的参考球面上，进而精确和快速反演地球重力场的方法。
二、背景技术
[0002] 如图I所示，21世纪是人类利用卫星跟踪卫星高低/低低技术（SST-HL/LL)提升对“数字地球”认知能力的新纪元。地球重力场及其随时间的变化反映地球表层及内部物质的空间分布、运动和变化，同时决定着大地水准面的起伏和变化。因此，确定地球重力场的精细结构及其时变不仅是大地测量、地球物理、海洋勘探、地震预报、空间科学、惯性导航、 航天技术、工业应用等的需求，同时也将为全人类寻求资源、保护环境和预测灾害提供重要的信息资源。第一，重力变化（重力异常）决定于地质构造和矿产资源（如煤炭、石油、天然气等）分布的不均匀性。由于我国远景资源储备中有78%的石油和93%的天然气亟待开发，因此精密重力测量有利于我国将来的资源勘探。第二，由于地震（如2005年3月印尼苏门答腊8. 5级地震、2008年5月中国汶川8. O级地震、2010年2月智利8. 8级地震、2011年3月日本福岛9. O级地震等）、火山、海嘯、龙卷风、泥石流等自然灾害发生前后，地球重力场将发生明显的变化和异常，因此精密重力测量能帮助我们预测自然灾害的发生，进而降低人员伤亡和经济损失。
[0003] 卫星重力反演是指通过分析轨道位置r和轨道速度^、星间距离P 12和星间速度Pu、非保守力f、三维姿态q、重力梯度Vij等卫星观测数据和地球引力位系数（地球引力位按球谐函数展开的常数Clm，Slm)的关系，建立并求解卫星运动观测方程，进而恢复地球引力位系数，最终目的是反演高精度和高空间分辨率的地球重力场。对卫星观测方程（大型线性超定方程组）的精确和快速求解是反演高精度地球重力场的关键技术。至今为止，国内外众多学者在地球重力场反演方法方面已开展了广泛研究。
[0004] I、直接最小二乘法（DLSP):通过正规方阵的直接求逆进而解算地球引力位系数。优点是求解过程严格，解算精度较高；缺点是大型矩阵直接求逆较困难，耗时较多，且需高性能计算机支持。
[0005] 在地心惯性系中，卫星观测方程建立如下
[0006] Gtxi = λ txl+etxl = AtxnXnxi (I)
[0007] 其中，Gtxi表示卫星观测值向量（主要包括轨道位置、轨道速度、星间速度、非保守力、三维姿态等），t表示卫星轨道观测点的数量；λ txl表示观测真值向量（观测真值向量中不含测量误差，而观测值向量中含有测量误差）；χηΧ1表示待求的地球引力位系数向量，其中引力位系数按照阶数I排列形成（次数m固定= +2Zmax _3表示待求引力位系数的个数，Lfflax表示地球引力位按球谐函数展开的最大阶数;AtXn表示t行η列的设计矩阵；etxl表示观测值误差向量，如果etxl是各元素相互独立的随机误差向量，卫星观测方程
(I)有如下统计特征
[0008] !DHxi (2)
[0009] 其中，E表示期望值算子，元xl表示Xnxi的真值。由于观测方程（I)是大型线性超定方程组，因此没有常规解，只有最小二乘解。在（I)式两边同乘4TX,,得
[0013] Atxn是一个庞大的长方转换矩阵，存储需占用大量的内存空间，因此直接存储较难实现。正规方阵Nnxn虽小于Atxn但如果直接存储也会占用大量的内存空间（解算120阶地球引力位系数至少耗费3Gb内存）。如果采用30天的卫星观测值，采样间隔为10s，则观测点个数为t = 259200个，若反演120阶地球重力场，则待求引力位系数的个数为n = L2maK+ 2LmsK — 3 = 14637。由⑷式可知，分别求解Nnxn和xnX1约需要tn2和η3量级的双精度运算，基于DLSP需要的总运算量约为tn2+n3，此运算量即使在目前超大型并行机上计算也非常耗时，因此利用DLSP求解大型线性超定方程组较困难。
[0014] 2、预处理共轭梯度迭代法（PCCG):通过优化选取预处理阵避免了大型矩阵的直接求逆，每步迭代的方向选择以误差最小为原则进而解算出地球引力位系数。优点是适当提高了计算速度，降低了对计算机性能的要求，适合于解算中长波重力场；缺点是选取预处理阵做了近似处理，但可以通过适当增加循环迭代的次数降低精度损失。
[0015] PCCG是目前求解大规模线性超定方程组的有效方法之一，在保证计算精度的前提下可显著提高运算速度且只需要几百Mb的内存空间。在PCCG中，Nnxn的性质至关重要，只有Nnxn满足正定满秩的条件，解算中迭代才可能收敛。主体思想如下：第一，每步迭代均对待求参量进行修正，直至达到预期精度为止；第二，每步迭代的方向选择以误差最小为原则；第三，回避最小二乘法的直接矩阵求逆，通过循环迭代求解真值。Nnxn是块对角占优结构的稠密阵，此性质为快速迭代求解提供了有利条件。由于不同的卫星设计具有不同的Nnxn,因此Nnxn的特殊性造成了预处理阵选择的差异性。由于预处理阵选择的优劣是整个大规模线性超定方程组解算的关键技术，因此对观测方程正规矩阵的优化处理是求解大规模线性超定方程组的首要问题。
[0016] PCCG最关键的部分在于预处理阵Mnxn的选取，优化选取有两个标准■o第一，M丄易于计算，可提高地球重力场反演的速度；第二，与尽可能接近，可保证地球重力场反演的精度。由于Nnxn是块对角占优方阵，因此通常选取Nnxn的块对角部分作为预处理阵，形成的Mnxn为主对角线上按次数m排列，其余部分为O的块对角方阵。如此选取不仅保留了 Nnxn主要特征，而且较易于计算。总之，优化选取预处理阵可较大程度地减少PCCG求解地球引力位系数中循环迭代的次数（计算时间较直接最小二乘法至少降低1000倍）。
[0017] 在方程⑷两边同乘得
[0019] 利用PCCG求解观测方程（5)，不需要存储NnXn，因此只需要几百Mb的内存空间，算法思想如下：
[0020]第一步,初始化：Xq = O, r。= y-Nx0, Z0 = M_1r0, d0 = z0,p0 = r0TZ0 ;
[0021] 第二步，循环迭代直到&
[0029] [8]返回到[3]循环迭代。
[0030] 综上所述，DLSP虽然理论上求解精度较高，但据目前国际超大并行计算机的解算速度而言暂时较难推广（仅美国宇航局喷气推进实验室（NASA-JPL)、德国波兹坦地学研究中心（GFZ)等国际研究机构采用），然而随着国际计算机性能的日益提高，将来有望逐步成为广泛应用的标准方法；PCCG是目前国际众多科研机构广泛采用的解算GRACE (GravityRecovery and Climate Experiment)地球重力场的有效方法之一，但其仅适合于求解中低阶地球重力场。
三、发明内容
[0031] 本发明的目的是：基于双星空间三维插值法有效加快解算速度，而且进一步提高地球重力场反演的精度。
[0032] 为达到上述目的，本发明采用了如下技术方案：
[0033] 基于双星空间三维插值原理的精确和快速卫星重力反演方法包含下列步骤：
[0034] I)利用GRACE卫星K波段测量仪获取星间速度、利用GPS接收机获取卫星轨道位置和卫星轨道速度、利用加速度计获取卫星非保守力、以及利用恒星敏感器获取卫星三维姿态观测数据，快速计算卫星轨道处的地球扰动位。在地固系中，地球扰动位T(r，θ，λ)按球谐函数展开的表达式为
[0036] 其中，GM表示地球质量M和万有引力常数G之积；RS表示地球的平均半径；r = yJX2+y2+ z2表示卫星的地心半径，X, y, Z分别表示卫星位置矢量r的三个分量；Θ和λ分别表示卫星的地心余纬度和地心经度；l(COsO表示规格化的Legendre函数，I表示阶数，m表示次数山表示地球引力位按球谐函数展开的最大阶数表示待求的规格化地球弓I力位系数。
[0037] 2)由于卫星轨道观测点的空间位置不规则，因此直接在卫星轨道处反演地球重力场耗费时间较大（基于CPU为3. OGHz,内存为4Gb的PC机，利用时间为30天和采样间隔为IOs的卫星观测数据，解算120阶地球重力场耗时至少20h)。为了克服上述耗时的缺点，本发明基于新型参考球面法快速解算了地球引力位系数。参考球面个数最优选取的原则如下：在插值误差对地球重力场解算精度不产生实质影响的标准下，尽量采用较少的参考球面个数，进而减少整体计算量和降低对计算机性能的要求。本发明基于双星空间三维插值法，利用时间为30天和采样间隔为IOs的卫星观测数据，并分别选取2个、3个、4个、5个和6个参考球面解算了 120阶GRACE地球重力场，结果表明：由于双星空间三维插值法引入的误差较GRACE卫星各关键载荷的测量误差约小10倍，不会对地球重力场反演精度产生实质性影响，因此选取4个参考球面较优。本发明以地心为球心选择4个等间距且规则的参考球面Ir1, r2, r3和r4,距地心最近参考球面的半径为!T1 = 6830km,球面间隔为Ar = 20km,卫星轨道H = 500km位于巧和!^参考球面之间，并在每个球面上按照分辨率=纬度Λ θ X经度Λ λ进行均匀网格划分，其中Λ Θ取为180° /2048，Δ λ取为360° /4096。
[0038] 3)利用卫星轨道上的I个地球扰动位观测值三维插值得到4个参考球面上I个单元格子的扰动位观测值。单元格子的分辨率设计为Λ θΧΛλΧΛι·，其中Λ Θ取为180 ° /2048，Δ λ取为360 ° /4096，Ar取为20km。计算三维插值
[0042]其中，
I, z)表示卫星轨道上的地球扰动位（本发明基于能量守恒法利用GRACE卫星K波段测量仪的星间速度、GPS接收机的卫星轨道位置和卫星轨道速度、加速度计的卫星非保守力、以及恒星敏感器的卫星三维姿态数据计算获得表示参考球面上的直角坐标，Wab=(a = 1,2,3,4 ；b = x, y, z)表示将卫星轨道上的地球扰动位三维插值到规则参考球
面上的权重系数，
表示卫星轨道上的三维直角
坐标，b2 = x2, y2, Z2表示第二个参考球面上的三维直角坐标，Aqn表示卫星轨道和第二个参考球面间的距离。
[0043] 4)由于Legendre函数(COs6&)的计算是耗时较多部分，因此固定每个规则参考球面上有限的纬度网格划分值Θ = 2048份，通过Legendre递推公式按纬圈快速计算^W(C0S6&)，并将计算结果存储以备后面调用，这是运算速度加快的主要原因。Legendre函数的快速求解是利用双星空间三维插值法求解大型线性超定方程组的关键问题之一。为了达到快速求解Legendre函数的目的，本发明采用如下方法：在保证地球重力场反演精度的前提下，由于Legendre函数只与余纬度Θ有关，因此不需要逐个准确计算出每个观测点的Legendre函数,只需计算有限个余纟韦度点的Legendre函数,而其它点上的Legendre函数值用已知点上的函数值作泰勒展开逼近即可。本发明首先将余纬度划分成2048等份；其次，计算每个小区间中点处的Legendre函数值并存储。具体过程如下：
[0044] [I]计算i时刻采样点的余纬度Θ i，并搜索所在的小区间为[Θ j，Θ J+1]；
[0045] [2 假定 Θ j+1/2 是区间[Θ」，Θ J+1]的中点，则 Plm (cos Θ j)的计算可用 Plm (cos Θ J+1/2)
的泰勒展开式逼近，
[0047] 其中，0( δ Θ L+1)表示Plm(cos Θ D按泰勒公式展开的高阶小量，δ Θ L+1 =
[0048] 正规化缔合Legendre函数及一阶导数的计算是影响大型线性超定方程组计算速度的重要因素。本发明将给出一种改进的正规化缔合Legendre函数的递推算法，具体公式如下：
[0049]假设 u = sin Θ , w = cos θ，则Q=I^1= λ/3μ，正规化缔合 Legendre 函数和
对Θ的一阶导数可表示为如下形式，
[0055] 5)利用双星空间三维插值法快速解算大型超定方程组（方程的个数远大于未知数的个数），进而拟合出地球引力位系数。在地心惯性系中，将地球扰动位卫星观测方程
(6)以矩阵形式表示如下
[0056] Dtxi = Atxn · xnXt (11)
[0057] 其中，Dtxi表示将卫星轨道处的地球扰动位观测值T(X，y, ζ)(本发明基于能量守恒法利用GRACE卫星K波段测量仪的星间速度、GPS接收机的卫星轨道位置和卫星轨道速度、加速度计的卫星非保守力、以及恒星敏感器的卫星三维姿态数据计算获得）空间三维插值到参考球面后的列矩阵，t表示观测点的数量；xnX1表示待求地球引力位系数向量(Clm，& J，其中引力位系数按照阶数I排列形成（次数m固定），《 = L2_ + 2Lmax - 3表示待求引力位系数的个数，Lfflax表示引力位按球谐函数展开的最大阶数;AtXn表示t行η列的设计矩阵（在地球扰动位观测方程（6)右边，除地球引力位系数外的其它部分）。在（11)式两边同乘τχ,,可得待求的地球引力位系数
[0061] 在求解大型线性方程组（13)时的快速解算是核心环节。由于对的直接
求解和存储都较困难，因此为了有效解算GRACE-type卫星跟踪卫星模式观测方程进而快速反演高精度和高空间分辨率地球重力场，本发明提出了前述利用双星空间三维插值法精确和快速解算的新思想。
[0062] 本发明是基于双星空间三维插值原理有利于精确和快速反演地球重力场的特点而设计的，优点是：1)地球重力场反演精度较高；2)卫星重力反演过程物理含义明确；3)在保证重力场反演精度的前提下，较大程度提高了计算速度；4)对计算机性能的要求相对较低；5)对地球重力场中短波信号较敏感。
四、附图说明
[0063] 图I表示卫星跟踪卫星高低/低低观测模式的测量原理，不仅以高轨GPS卫星对低轨重力双星实时精密跟踪定轨（SST-HL)，同时以差分原理测定低轨双星之间的相互运动(SST-LL)。
[0064] 图2表示美国宇航局（NASA)和德国航天局（DLR)合作发射的GRACE重力卫星。[0065] 图3表示由卫星轨道上的I个地球扰动位观测值三维插值到4个参考球面上的I个单兀格子。
[0066] 图4表示基于新型双星空间三维插值原理的卫星重力反演方法原理图。
[0067] 图5表示基于新型双星空间三维插值法反演地球重力场的精度，其中横坐标表示地球引力位按球函数展开的阶数，纵坐标表示大地水准面的累积误差（单位：m)。
五、具体实施方式
[0068] 以下结合附图，对本发明的具体实施方式作进一步的说明。
[0069] 基于新型双星空间三维插值原理的精确和快速卫星重力反演方法的应用：
[0070] I)利用GRACE卫星（如图2所示）K波段测量仪获取星间速度、利用GPS接收机获 取卫星轨道位置和卫星轨道速度、利用加速度计获取卫星非保守力、以及利用恒星敏感器获取卫星三维姿态观测数据，快速计算卫星轨道处的地球扰动位。在地固系中，地球扰动位T(r, θ，λ)按球谐函数展开的表达式为 [0071]
[0072] 其中，GM表示地球质量M和万有引力常数G之积；Re表示地球的平均半径；r =」χ2 + y2 + ζ2表示卫星的地心半径，χ，I, ζ分别表示卫星位置矢量r的三个分量；Θ和λ分别表示卫星的地心余纬度和地心经度；l(COS^)表示规格化的Legendre函数，I表示阶数，m表示次数山表示地球引力位按球谐函数展开的最大阶数；C/w^P*S/m表示待求的规格化地球弓I力位系数。
[0073] 2)由于卫星轨道观测点的空间位置不规则，因此直接在卫星轨道处反演地球重力场耗费时间较大（基于CPU为3. OGHz,内存为4Gb的PC机，利用时间为30天和采样间隔为IOs的卫星观测数据，解算120阶地球重力场耗时至少20h)。为了克服上述耗时的缺点，本发明基于新型参考球面法快速解算了地球引力位系数。参考球面个数最优选取的原则如下：在插值误差对地球重力场解算精度不产生实质影响的标准下，尽量采用较少的参考球面个数，进而减少整体计算量和降低对计算机性能的要求。本发明基于双星空间三维插值法，利用时间为30天和采样间隔为IOs的卫星观测数据，并分别选取2个、3个、4个、5个和6个参考球面解算了 120阶GRACE地球重力场，结果表明：由于双星空间三维插值法引入的误差较GRACE卫星各关键载荷的测量误差约小10倍，不会对地球重力场反演精度产生实质性影响，因此选取4个参考球面较优。本发明以地心为球心选择4个等间距且规则的参考球面Ir1, r2, r3和r4,距地心最近参考球面的半径为!T1 = 6830km,球面间隔为Ar = 20km,卫星轨道H = 500km位于巧和!^参考球面之间，并在每个球面上按照分辨率=纬度Λ θ X经度Λ λ进行均匀网格划分，其中Λ Θ取为180° /2048，Δ λ取为360° /4096。
[0074] 3)利用卫星轨道上的I个地球扰动位观测值三维插值得到4个参考球面上I个单元格子的扰动位观测值。如图3所示，单元格子的分辨率设计为Λ θ X Λ λ X Ar，其中Δ Θ取为180° /2048，Δ λ取为360° /4096，Δ λ取为20km。双星空间三维插值公式表示如下：
I, z)表示卫星轨道上的地球扰动位（本发明基于能量守恒法利用GRACE卫星K波段测量仪的星间速度、GPS接收机的卫星轨道位置和卫星轨道速度、加速度计的卫星非保守力、以及恒星敏感器的卫星三维姿态数据计算获得），（气,·^)表示参考球面上的直角坐标，wab(a = 1,2,3,4 ；b = x, y, ζ)表示将卫星轨道上的地球扰动位三维插值到规则参考球面
上的权重系数，Wli /?(l-p)2，w2b = ^~~P2 + 2^3 ,Wu =^?1 + 4Ρ_3/?2) ’
w4i =-^p2 ip-I) ~ ' A 分 2，P = (b_b2) /Δ qn, b = x, y, z 表示卫星轨道上的三维直角
坐标，b2 = x2, y2, Z2表示第二个参考球面上的三维直角坐标，Δ qn表示卫星轨道和第二个参考球面间的距离。将双星空间三维插值公式的权重系数以行矩阵^^乂^^^^^^^^表示，则双星空间三维插值公式矩阵形式表示如下：
[0081] 因此可得经过双星空间三维插值后4个参考球面上的扰动位
[0083] 4)由于Legendre函数(COs6&)的计算是耗时较多部分，因此固定每个规则参考球面上有限的纬度网格划分值Θ = 2048份，通过Legendre递推公式按纬圈快速计算
(cos(9)，并将计算结果存储以备后面调用，这是运算速度加快的主要原因。Legendre函数的快速求解是利用双星空间三维插值法求解大型线性超定方程组的关键问题之一。为了达到快速求解Legendre函数的目的，本发明采用如下方法：在保证地球重力场反演精度的前提下，由于Legendre函数只与余纬度Θ有关，因此不需要逐个准确计算出每个观测点的Legendre函数,只需计算有限个余纟韦度点的Legendre函数,而其它点上的Legendre函数值用已知点上的函数值作泰勒展开逼近即可。本发明首先将余纬度划分成2048等份；其次，计算每个小区间中点处的Legendre函数值并存储。具体过程如下：
[0084] [I]计算i时刻采样点的余纬度Θ i，并搜索所在的小区间为[Θ j，Θ i+J ；
[0085] [2]假定Θ i+1/2是区间[Θ P Θ i+1]的中点，则Plni(cos Θ j)的计算可用Plm (COS Θ i+1/2)的泰勒展开式逼近，[0086]
[0087] 其中，0( δ Θ L+1)表示Plm(cos Θ D按泰勒公式展开的高阶小量，δ Θ L+1 =
(e「0J+1/2)L+1。
[0088] 正规化缔合Legendre函数及一阶导数的计算是影响大型线性超定方程组计算速度的重要因素。本发明将给出一种改进的正规化缔合Legendre函数的递推算法，具体公式如下：
[0089]假设 u = sin θ , w = cos θ，则A。= I丨：=，正规化缔合 Legendre 函数和对Θ的一阶导数可表示为如下形式，
[0095] 5)利用双星空间三维插值法快速解算大型超定方程组（方程的个数远大于未知数的个数），进而拟合出地球引力位系数。在地心惯性系中，将地球扰动位卫星观测方程以矩阵形式表示如下
[0096] Dtxi = Atxn · Xnxi
[0097] 其中，Dtxi表示基于双星空间三维插值法将卫星轨道处的地球扰动位观测值T (X，1，ζ)(本发明基于能量守恒法利用GRACE卫星K波段测量仪的星间速度、GPS接收机的卫星轨道位置和卫星轨道速度、加速度计的卫星非保守力、以及恒星敏感器的卫星三维姿态数据计算获得）空间三维插值到参考球面后的列矩阵，t表示观测点的数量；Xnxl表示待求地球引力位系数向量(Ch,i/m)，其中引力位系数按照阶数I排列形成（次数m固定），^ = L2max + 2Lmax —3表示待:求弓丨力位;系数;的个数;，Lmax表示弓丨力位;按球i皆函数;展开的最大阶数；atxn表示t行η列的设计矩阵（在地球扰动位卫星观测方程右边，除了地球引力位系mclm,slm)^rn^m^)。在±式;两?同乘乂L可得待求的》球弓丨力位系数
[0098] ^xl = <<人)一1 . 4LAx!
[0099]令L = {Λ?Λ,ηΤ1 . ^ln，因此上式可变形为
[0101] 在求解大型线性方程组时，的快速解算是核心环节。如图4所示，由于对
的直接求解和存储都较困难，因此为了有效解算GRACE-type卫星跟踪卫星模式观测方程 进而快速反演高精度和高空间分辨率地球重力场，本发明提出了前述利用双星空间三维插
值法精确和快速解算的新思想。
[0102] 图5中星号线表示德国GFZ公布的120阶EIGEN-GRACE02S全球重力场模型的精度，在120阶处累计大地水准面精度为18. 938虚线表示利用双星空间三维插值法反演GRACE地球重力场的精度，在120阶处累计大地水准面精度为15.421cm。通过2条曲线在各阶处的符合性可知，双星空间三维插值法是解算高精度和高空间分辨率地球重力场的有效方法之一。本发明基于双星空间三维插值法反演GRACE地球重力场的精度在长波部分（LS 50阶）略优于EIGEN-GRACE02S全球重力场模型的精度，而在中长波部分（50 & LS 120阶）与其符合较好，原因分析如下：德国GFZ基于动力学法利用110天和采样间隔5s的GRACE卫星观测数据解算了 120阶地球重力场。据图5中星号线可知，EIGEN-GRACE02S模型长波重力场的精度并不十分理想，其主要原因是由于GRACE-A/B采用双星相互跟踪差分模式（星间距离220±50km)，在抵消掉双星间共同误差的同时，也差分掉了部分信号，因此导致对部分阶次地球引力位系数敏感性降低。另外，由于观测数据较复杂，在有效处理后仍存在较多不可知和不可控误差，因此会导致地球重力场长波部分信噪比较低。本发明基于GRACE实际参数获取了卫星观测数据，主要考虑了关键载荷误差（K波段测量仪、GPS接收机、加速度计、恒星敏感器等）、大气阻力误差、固体潮汐误差、太阳和月球引力误差等，但不可能将观测数据中存在的目前不可知和不可控误差都考虑在内，因此导致了重力场长波信号的精度（LS 50阶）略优于EIGEN-GRACE02S模型的精度。通过本发明的数值模拟检验，可预测将来随着观测数据处理方法和地球重力场反演方法的进一步优化和改进，GRACE地球重力场长波信号的精度仍有可提高的潜力。另外，我们基于半解析法利用美国JPL公布的GRACE卫星观测误差数据估计了 120阶GRACE地球重力场的精度，其结果在30阶内与EIGEN-GRACE02S模型仍存在一定程度的差异，因此地球重力场反演方法的不同也是引起解算结果差异的原因之一。
*中国科学院测量与地球物理研究所重力卫星编队轨道稳定性优化设计和精密反演地球重力场方法 *中国科学院测量与地球物理研究所重力卫星编队轨道稳定性优化设计和精密反演地球重力场方法International ClassificationC06PublicationC10Request of examination as to substanceC14GrantedC17Cessation of patent rightRotate