基于双星空间三维插值原理的卫星重力反演方法技术

本发明专利技术涉及一种精密测量地球重力场的方法，特别是一种基于双星空间三维插值原理的卫星重力反演方法；基于将位于空间位置相对不规则的重力双星轨道上的地球扰动位(由GRACE卫星K波段测量仪的星间速度、GPS接收机的卫星轨道位置和卫星轨道速度、加速度计的卫星非保守力、以及恒星敏感器的卫星三维姿态数据计算获得)三维插值到空间位置相对规则且均匀网格划分的参考球面上，进而精确和快速反演了地球重力场；该方法卫星重力反演精度高，解算过程物理含义明确，计算速度快，计算机性能要求低，敏感于重力场中高频信号；因此，双星空间三维插值法是解算高精度和高空间分辨率地球重力场的有效方法。

【技术实现步骤摘要】

本专利技术涉及卫星大地测量学、地球物理学、空间科学等交叉
，特别是涉及一种基于将位于空间位置相对不规则的重力双星轨道上的地球扰动位(由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是目前本文档来自技高网...

【技术保护点】
１．一种基于双星空间三维插值原理的卫星重力反演方法，其特征如下：步骤一：利用ＧＲＡＣＥ卫星Ｋ波段测量仪获取星间速度、利用ＧＰＳ接收机获取卫星轨道位置和速度、利用加速度计获取卫星非保守力、以及利用恒星敏感器获取卫星三维姿态，通过上述数据计算卫星轨道处的地球扰动位Ｔ（ｘ，ｙ，ｚ）；步骤二：以地心为球心选择４个等间距且规则的参考球面ｒ１，ｒ２，ｒ３和ｒ４，距地心最近参考球面的半径为ｒ１＝６８３０ｋｍ，球面间隔为Δｒ＝２０ｋｍ，轨道高度Ｈ＝５００ｋｍ，卫星轨道位于ｒ２和ｒ４参考球面之间，在每个球面上按照分辨率＝纬度Δθ×经度Δλ进行均匀网格划分，其中Δθ取为１８０°／２０４８，Δλ取为３６０°／４０９６；步骤三：利用卫星轨道上测量的每个地球扰动位Ｔ（ｘ，ｙ，ｚ）三维插值得到４个参考球面上对应的单元格子的地球扰动位Ｔ（ｒ，θ，λ）；单元格子的分辨率设计为Δθ×Δλ×Δｒ，其中Δθ取为１８０°／２０４８，Δλ取为３６０°／４０９６，Δｒ取为２０ｋｍ；步骤四：将纬度θ划分为２０４８等份，并固定每个规则参考球面上有限的纬度网格划分值；其次，计算每个小区间中点处的Ｌｅｇｅｎｄｒｅ函数值并存储，其它点上的Ｌｅｇｅｎｄｒｅ函数值用已知点上的函数值作泰勒展开逼近获得，并将计算结果存储以备后面调用；步骤五：利用空间三维插值到规则的４个参考球面上的地球扰动位Ｔ（ｒ，θ，λ）解算大型超定方程组（１），最终获得地球引力位系数；（ｍａｔｈ）??（ｍｒｏｗ）?（ｍｉ）Ｔ（／ｍｉ）?（ｍｒｏｗ）?（ｍｏ）（（／ｍｏ）?（ｍｉ）ｒ（／ｍｉ）?（ｍｏ），（／ｍｏ）?（ｍｉ）＆ｔｈｅｔａ；（／ｍｉ）?（ｍｏ），（／ｍｏ）?（ｍｉ）＆ｌａｍｂｄａ；（／ｍｉ）?（ｍｏ））（／ｍｏ）?（／ｍｒｏｗ）?（ｍｏ）＝（／ｍｏ）?（ｍｆｒａｃ）?（ｍｉ）ＧＭ（／ｍｉ）?（ｍｓｕｂ）?（ｍｉ）Ｒ（／ｍｉ）?（ｍｉ）ｅ（／ｍｉ）?（／ｍｓｕｂ）?（／ｍｆｒａｃ）?（ｍｕｎｄｅｒｏｖｅｒ）?（ｍｉ）＆Ｓｉｇｍａ；（／ｍｉ）?（ｍｒｏｗ）?（ｍｉ）ｌ（／ｍｉ）?（ｍｏ）＝（／ｍｏ）?（ｍｎ）２（／ｍｎ）?（／ｍｒｏｗ）?（ｍｉ）Ｌ（／ｍｉ）?（／ｍｕｎｄｅｒｏｖｅｒ）?（ｍｓｕｐ）?（ｍｒｏｗ）?（ｍｏ）（（／ｍｏ）?（ｍｆｒａｃ）?（ｍｓｕｂ）?（ｍｉ）Ｒ（／ｍｉ）?（ｍｉ）ｅ（／ｍｉ）?（／ｍｓｕｂ）?（ｍｉ）ｒ（／ｍｉ）?（／ｍｆｒａｃ）?（ｍｏ））（／ｍｏ）?（／ｍｒｏｗ）?（ｍｒｏｗ）?（ｍｉ）ｌ（／ｍｉ）?（ｍｏ）＋（／ｍｏ）?（ｍｎ）１（／ｍｎ）?（／ｍｒｏｗ）?（／ｍｓｕｐ）?（ｍｕｎｄｅｒｏｖｅｒ）?（ｍｉ）＆Ｓｉｇｍａ；（／ｍｉ）?（ｍｒｏｗ）?（ｍｉ）ｍ（／ｍｉ）?（ｍｏ）＝（／ｍｏ）?（ｍｎ）０（／ｍｎ）?（／ｍｒｏｗ）?（ｍｉ）ｌ（／ｍｉ）?（／ｍｕｎｄｅｒｏｖｅｒ）?（ｍｒｏｗ）?（ｍｏ）（（／ｍｏ）?（ｍｓｕｂ）?（ｍｏｖｅｒ）?（ｍｉ）Ｃ（／ｍｉ）?（ｍｏ）＆ＯｖｅｒＢａｒ；（／ｍｏ）?（／ｍｏｖｅｒ）?（ｍｉ）ｌｍ（／ｍｉ）?（／ｍｓｕｂ）?（ｍｉ）ｃｏｓ（／ｍｉ）?（ｍｉ）ｍ＆ｌａｍｂｄａ；（／ｍｉ）?（ｍｏ）＋（／ｍｏ）?（ｍｓｕｂ）?（ｍｏｖｅｒ）?（ｍｉ）Ｓ（／ｍｉ）?（ｍｏ）＆ＯｖｅｒＢａｒ；（／ｍｏ）?（／ｍｏｖｅｒ）?（ｍｉ）ｌｍ（／ｍｉ）?（／ｍｓｕｂ）?（ｍｉ）ｓｉｎ（／ｍｉ）?（ｍｉ）ｍ＆ｌａｍｂｄａ；（／ｍｉ）?（ｍｏ））（／ｍｏ）?（／ｍｒｏｗ）?（ｍｓｕｂ）?（ｍｏｖｅｒ）?（ｍｉ）Ｐ（／ｍｉ）?（ｍｏ）＆ＯｖｅｒＢａｒ；（／ｍｏ）?（／ｍｏｖｅｒ）?（ｍｉ）ｌｍ（／ｍｉ）?（／ｍｓｕｂ）?（ｍｒｏｗ）?（ｍｏ）（（／ｍｏ）?（ｍｉ）ｃｏｓ（／ｍｉ）?（ｍｉ）＆ｔｈｅｔａ；（／ｍｉ）?（ｍｏ））（／ｍｏ）?（／ｍｒｏｗ）?（ｍｏ）－（／ｍｏ）?（ｍｏ）－（／ｍｏ）?（ｍｏ）－（／ｍｏ）?（ｍｒｏｗ）?（ｍｏ）（（／ｍｏ）?（ｍｎ）１（／ｍｎ）?（ｍｏ））（／ｍｏ）?（／ｍｒｏｗ）?（／ｍｒｏｗ）?（／ｍａｔｈ）其中，ＧＭ表示地球质量Ｍ和万有引力常数Ｇ之积；Ｒｅ表示地球的平均半径；表示卫星的地心半径，ｘ，ｙ，ｚ分别表示卫星位置矢量ｒ的三个分量；θ和λ分别表示卫星的地心余纬度和地心经度；表示规格化的Ｌｅｇｅｎｄｒｅ函数，ｌ表示阶数，ｍ表示次数；Ｌ表示地球引力位按球谐函数展开的最大阶数；和表示待求的规格化地球引力位系数。...

【技术特征摘要】
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...

【专利技术属性】
技术研发人员：郑伟，许厚泽，熊熊，钟敏，刘成恕，
申请(专利权)人：中国科学院测量与地球物理研究所，
类型：发明
国别省市：83