专利名称：基于双星空间三维插值原理的卫星重力反演方法
技术领域：
本发明涉及卫星大地测量学、地球物理学、空间科学等交叉技术领域，特别是涉及一种基于将位于空间位置相对不规则的重力双星轨道上的地球扰动位(由GRACE卫星K波段测量仪的星间速度、GPS接收机的卫星轨道位置和卫星轨道速度、加速度计的卫星非保守力、以及恒星敏感器的卫星三维姿态数据计算获得)三维插值到空间位置相对规则且均勻网格划分的参考球面上，进而精确和快速反演地球重力场的方法。
背景技术：
如
图1所示，21世纪是人类利用卫星跟踪卫星高低/低低技术(SST-HL/LL)提升对“数字地球”认知能力的新纪元。地球重力场及其随时间的变化反映地球表层及内部物质的空间分布、运动和变化，同时决定着大地水准面的起伏和变化。因此，确定地球重力场的精细结构及其时变不仅是大地测量、地球物理、海洋勘探、地震预报、空间科学、惯性导航、 航天技术、工业应用等的需求，同时也将为全人类寻求资源、保护环境和预测灾害提供重要的信息资源。第一，重力变化(重力异常)决定于地质构造和矿产资源(如煤炭、石油、天然气等)分布的不均勻性。由于我国远景资源储备中有78%的石油和93%的天然气亟待开发，因此精密重力测量有利于我国将来的资源勘探。第二，由于地震(如2005年3月印尼苏门答腊8. 5级地震、2008年5月中国汶川8. 0级地震、2010年2月智利8. 8级地震、2011 年3月日本福岛9.0级地震等)、火山、海啸、龙卷风、泥石流等自然灾害发生前后，地球重力场将发生明显的变化和异常，因此精密重力测量能帮助我们预测自然灾害的发生，进而降低人员伤亡和经济损失。卫星重力反演是指通过分析轨道位置r和轨道速度广、星间距离P 12和星间速度 Pu、非保守力f、三维姿态q、重力梯度Vij等卫星观测数据和地球引力位系数(地球引力位按球谐函数展开的常数Clm，Slm)的关系，建立并求解卫星运动观测方程，进而恢复地球引力位系数，最终目的是反演高精度和高空间分辨率的地球重力场。对卫星观测方程(大型线性超定方程组)的精确和快速求解是反演高精度地球重力场的关键技术。至今为止，国内外众多学者在地球重力场反演方法方面已开展了广泛研究。1、直接最小二乘法(DLSP)通过正规方阵的直接求逆进而解算地球引力位系数。 优点是求解过程严格，解算精度较高；缺点是大型矩阵直接求逆较困难，耗时较多，且需高性能计算机支持。在地心惯性系中，卫星观测方程建立如下Gtxi = λ txl+etxl = AtxnXnxi(1)其中，Gtxi表示卫星观测值向量(主要包括轨道位置、轨道速度、星间速度、非保守力、三维姿态等)，t表示卫星轨道观测点的数量；λ txl表示观测真值向量(观测真值向量中不含测量误差，而观测值向量中含有测量误差)；Xnxi表示待求的地球引力位系数向量， 其中引力位系数按照阶数1排列形成(次数m固定Km = Z^ix+2Zmax-3表示待求引力位系数的个数，Lfflax表示地球引力位按球谐函数展开的最大阶数；Atxn表示t行η列的设计矩阵；etxl表示观测值误差向量，如果etxl是各元素相互独立的随机误差向量，卫星观测方程 (1)有如下统计特征 = (2) 1取〗)=0 其中，E表示期望值算子，&χ1表示Ixi的真值。由于观测方程(1)是大型线性超定方程组，因此没有常规解，只有最小二乘解。在α)式两边同乘得 _ο] Α βιΛ = AlnAtxnXnxl(3)
_ 1 ]令凡X1 = ^lnGm ,Nnxn = A^xnAtxn，则ynX1 = NnxnXnxi(4)Atxn是一个庞大的长方转换矩阵，存储需占用大量的内存空间，因此直接存储较难实现。正规方阵Nnxn虽小于Atxn但如果直接存储也会占用大量的内存空间(解算120 阶地球引力位系数至少耗费3( 内存)。如果采用30天的卫星观测值，采样间隔为10s， 则观测点个数为t = 259200个，若反演120阶地球重力场，则待求引力位系数的个数为 n = L2maK+ 2LmsK - 3 = 14637。由⑷式可知，分别求解Nnxn和^cnxi约需要tn2和η3量级的双精度运算，基于DLSP需要的总运算量约为tn2+n3，此运算量即使在目前超大型并行机上计算也非常耗时，因此利用DLSP求解大型线性超定方程组较困难。2、预处理共轭梯度迭代法(PCCG)通过优化选取预处理阵避免了大型矩阵的直接求逆，每步迭代的方向选择以误差最小为原则进而解算出地球引力位系数。优点是适当提高了计算速度，降低了对计算机性能的要求，适合于解算中长波重力场；缺点是选取预处理阵做了近似处理，但可以通过适当增加循环迭代的次数降低精度损失。PCCG是目前求解大规模线性超定方程组的有效方法之一，在保证计算精度的前提下可显著提高运算速度且只需要几百Mb的内存空间。在PCCG中，Nnxn的性质至关重要， 只有Nnxn满足正定满秩的条件，解算中迭代才可能收敛。主体思想如下第一，每步迭代均对待求参量进行修正，直至达到预期精度为止；第二，每步迭代的方向选择以误差最小为原则；第三，回避最小二乘法的直接矩阵求逆，通过循环迭代求解真值。Nnxn是块对角占优结构的稠密阵，此性质为快速迭代求解提供了有利条件。由于不同的卫星设计具有不同的 Nnxn，因此Nnxn的特殊性造成了预处理阵选择的差异性。由于预处理阵选择的优劣是整个大规模线性超定方程组解算的关键技术，因此对观测方程正规矩阵的优化处理是求解大规模线性超定方程组的首要问题。PCCG最关键的部分在于预处理阵Mnxn的选取，优化选取有两个标准■ 第一，M丄易于计算，可提高地球重力场反演的速度；第二，Μ:1"与Λ ^尽可能接近，可保证地球重力场反演的精度。由于Nnxn是块对角占优方阵，因此通常选取Nnxn的块对角部分作为预处理阵，形成的Mnxn为主对角线上按次数m排列，其余部分为0的块对角方阵。如此选取不仅保留了 Nnxn主要特征，而且Mn^较A^1n易于计算。总之，优化选取预处理阵可较大程度地减少 PCCG求解地球引力位系数中循环迭代的次数(计算时间较直接最小二乘法至少降低1000 倍)。在方程(4)两边同乘Ma^得
5_ 8] M'nlnynxX = M=NrixnXrixl(5 )利用PCCG求解观测方程(5)，不需要存储Nnxn，因此只需要几百Mb的内存空间，算法思想如下第一步，初始化 = 0，r。= y-Nx0, Z0 = M^r0, d0 = z0,p0 = Tq1Z0 ；第二步，循环迭代直到^^ < ε，[1] q( =^iaiCii) (i = 0,1,2,…)，[2] a. = PiId^qi ,[3]xiX1 = Xi+α jdi,[4]ri+1 = r-a iqi,[5]zi+1 = M_1ri+1,[6] pi+1 = r^zl+l，[7]di+1 = zi+1+p i+1/p idi,[8]返回到[3]循环迭代。综上所述，DLSP虽然理论上求解精度较高，但据目前国际超大并行计算机的解算速度而言暂时较难推广(仅美国宇航局喷气推进实验室(NASA-JPL)、德国波兹坦地学研究中心(GFZ)等国际研究机构采用)，然而随着国际计算机性能的日益提高，将来有望逐步成为广泛应用的标准方法；PCCG是目前国际众多科研机构广泛采用的解算GRACE (Gravity Recovery and Climate Experiment)地球重力场的有效方法之一，但其仅适合于求解中低阶地球重力场。
发明内容
本发明的目的是基于双星空间三维插值法有效加快解算速度，而且进一步提高地球重力场反演的精度。为达到上述目的，本发明采用了如下技术方案基于双星空间三维插值原理的精确和快速卫星重力反演方法包含下列步骤1)利用GRACE卫星K波段测量仪获取星间速度、利用GPS接收机获取卫星轨道位置和卫星轨道速度、利用加速度计获取卫星非保守力、以及利用恒星敏感器获取卫星三维姿态观测数据，快速计算卫星轨道处的地球扰动位。在地固系中，地球扰动位T(r，θ，λ) 按球谐函数展开的表达式为
GM 1 (R V+1 1 ———
Τ{ν,θ,λ) = -- \ 丄 X{Clm cosπιλ + 4 sinmX)Plm(cosΘ)(6)
Ae 1=2 \ r J m=0其中，GM表示地球质量M和万有引力常数G之积武表示地球的平均半径； r = yJx2+y2 + z2表示卫星的地心半径，X, y，Z分别表示卫星位置矢量r的三个分量；θ 和λ分别表示卫星的地心余纬度和地心经度ilCcos^)表示规格化的Legendre函数，1 表示阶数，m表示次数；L表示地球引力位按球谐函数展开的最大阶数表示待求的规格化地球引力位系数。2)由于卫星轨道观测点的空间位置不规则，因此直接在卫星轨道处反演地球重力场耗费时间较大(基于CPU为3. OGHz,内存为4( 的PC机，利用时间为30天和采样间隔为IOs的卫星观测数据，解算120阶地球重力场耗时至少20h)。为了克服上述耗时的缺点， 本发明基于新型参考球面法快速解算了地球引力位系数。参考球面个数最优选取的原则如下在插值误差对地球重力场解算精度不产生实质影响的标准下，尽量采用较少的参考球面个数，进而减少整体计算量和降低对计算机性能的要求。本发明基于双星空间三维插值法，利用时间为30天和采样间隔为IOs的卫星观测数据，并分别选取2个、3个、4个、5个和 6个参考球面解算了 120阶GRACE地球重力场，结果表明由于双星空间三维插值法引入的误差较GRACE卫星各关键载荷的测量误差约小10倍，不会对地球重力场反演精度产生实质性影响，因此选取4个参考球面较优。本发明以地心为球心选择4个等间距且规则的参考球面r1 r2, r3和r4，距地心最近参考球面的半径为Γι = 6830km，球面间隔为Ar = 20km, 卫星轨道H = 500km位于巧和！^参考球面之间，并在每个球面上按照分辨率=纬度Δ θ X 经度Δ λ进行均勻网格划分，其中Δ θ取为180° /2048，Δ λ取为360° /4096。3)利用卫星轨道上的1个地球扰动位观测值三维插值得到4个参考球面上1个单元格子的扰动位观测值。单元格子的分辨率设计为Δ θΧΔλΧΔι·，其中Δ θ取为180 ° /2048，Δ λ取为360 ° /4096，Ar取为20km。计算三维插值
4 4 4
Τ(χ^^ζ) = ΣΣΣ Viyywi7(\‘^‘z^) (X，y，ζ) e [X2，X3] X [y2，y3] X [z2 z3] ',=1^=1 ‘=1
444
二ΣmVΣmVΣ，yiy,a),yjy，。)+,a) + w4zT(xix，^, )](7)
'^=I ',=1v ' y
4= ΣiwlyB(\,乃)+ W2yBixi ,y2) + W3yBixi，JV3) +,y4)]
',=1 * ‘ =WlxC(X1) + W2xC(X2) + W3xC(X3) + W4xC(X4)其中，c(Xix)= wlyB(xix,乃)+ w2yB(xir,y2) + w3yB(xix,y3) + w4yB(x、,y4),
,yly) = ，A) + ，^^2) + ，;ν4) + ,·>νA)，τ (x，
y，z)表示卫星轨道上的地球扰动位(本发明基于能量守恒法利用GRACE卫星K波段测量仪的星间速度、GPS接收机的卫星轨道位置和卫星轨道速度、加速度计的卫星非保守力、以及恒星敏感器的卫星三维姿态数据计算获得)，(\,少,2表示参考球面上的直角坐标，Wab =(a = 1,2,3,4 ；b = χ, y，ζ)表示将卫星轨道上的地球扰动位三维插值到规则参考球
面上的权重系数，W16 =-去/7(1 -ρ)2 ,w2b = 1 -昏ρ2 + 昏ρ3， =+ 4p-3p2),
w --P2(P-I) == (b-b2)/Aqn,b = χ, y，ζ 表示卫星轨道上的三维直角
坐标，b2 = x2, y2, &表示第二个参考球面上的三维直角坐标，Aqn表示卫星轨道和第二个参考球面间的距离。4)由于Legendre函数^m(COs6>)的计算是耗时较多部分，因此固定每个规则参考球面上有限的纬度网格划分值θ = 2048份，通过Legendre递推公式按纬圈快速计算 PjcosΘ)，并将计算结果存储以备后面调用，这是运算速度加快的主要原因。Legendre函数的快速求解是利用双星空间三维插值法求解大型线性超定方程组的关键问题之一。为了达到快速求解Legendre函数的目的，本发明采用如下方法在保证地球重力场反演精度的前提下，由于Legendre函数只与余纬度θ有关，因此不需要逐个准确计算出每个观测点的 Legendre函数，只需计算有限个余纬度点的Legendre函数，而其它点上的Legendre函数值用已知点上的函数值作泰勒展开逼近即可。本发明首先将余纬度划分成2048等份；其次， 计算每个小区间中点处的Legendre函数值并存储。具体过程如下[1]计算i时刻采样点的余纬度θ i，并搜索所在的小区间为[θ ρ θ J+1]；[2 假定 θ j+1/2 是区间[θ」，θ J+1]的中点，则 Plm(cos θ J 的计算可用 Plm(cos θ J+1/2)
的泰勒展开式逼近，
P1X) = Plm(cos0J+y2) + P^{cose^2Xei -6J+lll)k] +(8)其中，0( δ θ L+1)表示Plm(cos θ J按泰勒公式展开的高阶小量，δ θ L+1 =
(θ i_ θ j+1/2) L+1。正规化缔合Legendre函数及一阶导数的计算是影响大型线性超定方程组计算速度的重要因素。本发明将给出一种改进的正规化缔合Legendre函数的递推算法，具体公式如下假设u = sin θ，w = C0s θ，则巧。=1 j = ^u，正规化缔合 Legendre 函数和
对θ的一阶导数可表示为如下形式，
权利要求
1.一种基于双星空间三维插值原理的卫星重力反演方法，其特征如下步骤一利用GRACE卫星K波段测量仪获取星间速度、利用GPS接收机获取卫星轨道位置和速度、利用加速度计获取卫星非保守力、以及利用恒星敏感器获取卫星三维姿态，通过上述数据计算卫星轨道处的地球扰动位T (x，y，z)；步骤二 以地心为球心选择4个等间距且规则的参考球面r1; r2，r3和r4，距地心最近参考球面的半径为T1 = 6830km，球面间隔为Ar = 20km，轨道高度H = 500km，卫星轨道位于 ! 2和1~4参考球面之间，在每个球面上按照分辨率=纬度Δ θ X经度Δ λ进行均勻网格划分，其中 Δ θ 取为 180° /2048，Δ λ 取为 360° /4096 ；步骤三利用卫星轨道上测量的每个地球扰动位T(x，y, ζ)三维插值得到4个参考球面上对应的单元格子的地球扰动位T(r，θ , λ)；单元格子的分辨率设计为 Δ θ X Δ λ X Δι·，其中 Δ θ 取为 180° /2048，Δ λ 取为 360° /4096，Ar 取为 20km ；步骤四将纬度θ划分为2048等份，并固定每个规则参考球面上有限的纬度网格划分值；其次，计算每个小区间中点处的Legendre函数值并存储，其它点上的Legendre函数值 gm(cos60用已知点上的函数值作泰勒展开逼近获得，并将计算结果存储以备后面调用；步骤五利用空间三维插值到规则的4个参考球面上的地球扰动位T (r，θ , λ)解算大型超定方程组(1)，最终获得地球引力位系数；
2.如权利要求1所述的卫星重力反演方法，其特征在于步骤三中的插值公式为
3.如权利要求1所述的卫星重力反演方法，其特征在于在步骤五中解算大型超定方程组为将(1)式以矩阵形式表示为DtXl — AtXn · ΧηΧ 1(^)其中，Dtxi表示基于双星空间三维插值法将卫星轨道处的地球扰动位观测值T(X，y， ζ)(本发明基于能量守恒法利用GRACE卫星K波段测量仪的星间速度、GPS接收机的卫星轨道位置和卫星轨道速度、加速度计的卫星非保守力、以及恒星敏感器的卫星三维姿态数据计算获得)空间三维插值到参考球面后的列矩阵，t表示观测点的数量；Xnxl表示待求地球引力位系数向量(C/w,$w)，其中引力位系数按照阶数1排列形成(次数m固定)， n^ L2max+ Hmax - 3表示待求引力位系数的个数，Lfflax表示引力位按球函数展开的最大阶数；Atxn表示t行η列的设计矩阵(在地球扰动位卫星观测方程(1)右边，除了地球引力位系数外的其它部分)；在上式两边同乘可得待求的地球引力位系数
全文摘要
本发明涉及一种精密测量地球重力场的方法，特别是一种基于双星空间三维插值原理的卫星重力反演方法；基于将位于空间位置相对不规则的重力双星轨道上的地球扰动位(由GRACE卫星K波段测量仪的星间速度、GPS接收机的卫星轨道位置和卫星轨道速度、加速度计的卫星非保守力、以及恒星敏感器的卫星三维姿态数据计算获得)三维插值到空间位置相对规则且均匀网格划分的参考球面上，进而精确和快速反演了地球重力场；该方法卫星重力反演精度高，解算过程物理含义明确，计算速度快，计算机性能要求低，敏感于重力场中高频信号；因此，双星空间三维插值法是解算高精度和高空间分辨率地球重力场的有效方法。
文档编号G01V7/00GK102262248SQ20111015003
公开日2011年11月30日 申请日期2011年6月3日 优先权日2011年6月3日
发明者刘成恕, 熊熊, 许厚泽, 郑伟, 钟敏 申请人:中国科学院测量与地球物理研究所