CN108267792B - 全球重力场模型反演方法 - Google Patents

Abstract

本发明公开了一种全球重力场模型反演方法，包括步骤1：采集低轨重力卫星数据并进行预处理；步骤2：通过数值微分计算纯地球引力加速度观测值；步骤3：观测值定权策略构造新的观测方程；步骤4：参数估计得到地球重力场模型位系数。本发明基于参考模型的加速度域残差分析，直接采用自相关序列作为基本元素来构造对称Toeplitz协方差阵，可以有效避免由轨道先验稠密阵经误差传播定律计算加速度值协方差阵时高维大型矩阵的复杂运算过程，减少定权耗时，从而提高重力场模型解算效率。相较于经典方法，本发明提出的算法操作简洁高效，而且整体上能够保证重力场模型反演计算的效率稳定可靠。

Description

全球重力场模型反演方法

技术领域

本发明涉及卫星重力探测技术领域，更具体的涉及一种全球重力场模型反演方法。

背景技术

地球重力场是地球内部物质分布及质量迁移的历史再现，确定地球重力场的精细结构及其时空变化具有重要的科学和现实意义，因此地球重力场的研究一直是大地测量学领域核心和热点问题。在众多地球重力场模型解算方法中，基于高低卫星跟踪卫星的加速度法由于原理简单计算耗时少而被广泛采用。加速度法是由卫星星历通过数值微分计算得到卫星运动载体加速度值，基于牛顿第二运动定律直接建立卫星高度处加速度观测值与重力场模型位系数之间的线性关系式来估算位系数，为获得最优解需要合理定权，目前最为流行的是基于轨道先验稠密阵误差传播定律的经典方法。其中，轨道先验稠密阵是指定轨确定的先验协方差阵，受GNSS星座影响导致矩阵具有稠密性 (非零元素占优)；轨道位置协方差阵经误差传播定律即可获得加速度观测值协方差信息。

采用经典方法处理时，首先需要已知轨道协方差信息，然后基于误差传播定律获得加速度协方差阵。但是当轨道精度较差时，误差传播后的加速度协方差阵信息较差，导致估算的位系数偏差大，反演地球重力场模型分辨率差，同时需解算高维大型矩阵，运算过程较为复杂，导致重力场反演计算效率低下。

发明内容

本发明要解决的技术问题在于针对现有技术中加速度法在构建观测值权阵时，由于矩阵维数较高，处理运算过程复杂，导致重力场模型反演计算效率低下的缺陷，提供一种能够有效减少定权耗时，提高模型解算效率，且稳定可靠的全球重力场模型反演方法。

本发明解决其技术问题所采用的技术方案是：

提高一种全球重力场模型反演方法，包括以下步骤：

步骤1：采集低轨重力卫星数据并进行预处理；

步骤2：数值微分计算

采用移动中心窗口多项式微分法进行微分计算，并在此基础上扩展微分滤波器，抑制高频噪声的放大影响；根据步骤1的预处理结果和摄动力建模计算各项摄动加速度，同时基于星载加速度计确定非保守力；将载体运动加速度减去除地球引力外的保守力摄动加速度和非保守力摄动加速度，得到纯地球引力加速度观测值；

步骤3：观测值定权策略

引入参考重力场模型仿真模拟重力加速度，对其进行坐标变换旋转到惯性系下的重力加速度模拟值；将纯地球引力加速度观测值与重力加速度模拟值作粗差探测进行数理统计确定阈值，进而剔除超出阈值的信号；利用残差分析确定纯地球引力加速度观测值自相关序列，并将其视作基本元素，以此构造观测值自协方差阵；针对观测值有色噪声特性，对其观测值协方差阵Q 进行Cholesky分解Q＝R^TR,令F＝(R^T)^-1则重新构造新的观测方程y^*＝B^*x+Δ^*，使新的观测值y^*＝Fy满足白噪声(Q_y*＝E)，继而通过经典最小二乘准则解算位系数来提高计算效率；其中R和F分别是观测值协方差阵Q的乔列斯基分解矩阵；y*、B*和Δ^*分别表示新的观测量、设计阵和噪声；Q_y*表示新观测值的协方差阵，E为单位阵；

步骤4：参数估计

针对海量卫星观测数据，实际处理时将整个观测数据分割成若干短弧段分配线程，通过各个线程协同分别计算各弧段对应设计矩阵B_i和观测量y_i，组成若干观测方程y_i＝B_ix+Δ_i，经步骤3形成新的观测方程据此建立各弧段法方程组基于并行算法采用多线程同时计算各弧段法方程阵和继而对各弧段的法方程阵进行叠加以减少矩阵的存储空间，缩减计算耗时；顾及位系数以次为主的顺序排列方式具有块对角占优的特性，采用所述参考重力场模型进行移去-恢复技术，利用最小二乘平差最终解算得到地球重力场模型位系数。

接上述技术方案，采集的低轨重力卫星数据包括几何法轨道观测值、简化动力学轨道观测值和星载加速度计观测值。

接上述技术方案，所述预处理包括粗差探测和时变重力场改正两部分：采用数理统计方法对几何法轨道进行粗差探测，超出阈值部分直接用对应的简化动力学轨道填充；采用IERS提供的背景模型分离出低轨重力卫星数据中的时变重力场信号，获取高精度的静态重力场模型。

接上述技术方案，对步骤3中的观测值定权的计算过程进行优化：

重力加速度模拟值IRF:

局部轨道坐标系下加速度模拟值：

加速度观测值的协方差阵：

其中Δ_i分别表示i时刻卫星在惯性系IRF的微分加速度观测矢量、地球引力矢量和观测噪声，和分别表示扣除地球引力外的保守力和非保守力；(r,θ,λ)分别表示地固球坐标的地心向径、地心余纬和地心经度；GM和R分别指地心引力常数和地球平均半径；n和m为球谐展开的阶和次；N为重力场模型的最大阶数；和表示完全规格化的n阶m次球谐位系数；和分别表示完全规格化的n阶m次缔合勒让德函数及其一阶导数；V_x、V_y和V_z分别表示局部轨道坐标系LNRF下沿轨、法向和径向方向处的加速度模拟值；{σ_i}表示参考重力场模型正演的加速度模拟值与轨道经数值微分扣除各项摄动力后的加速度观测值作差所得到的残差序列；

具体包括：

首先基于参考重力场模型进行粗差探测，将先验重力场模型正演的重力加速度模拟值与轨道经数值微分扣除各项摄动力后的加速度观测值作差得到残差序列{σ_i}，数理统计后设置扰动重力阈值；对于超出阈值的历元处，采用参考重力场模型正演的加速度模拟值直接填充对应观测值，此时对应的残差值置零；然后以将残差序列为自相关基本元素，经正规化处理构造观测值的对称Toeplitz协方差阵。

接上述技术方案，移动中心窗口多项式微分法是以计算点为中心两侧各取等点构成窗口函数结合多项式对窗口内的观测点进行最小二乘拟合通过线性组合滤波核函数与观测值来计算中心值；随后移动窗口，接着计算该点的下一点，依次类推，组成平滑后的波谱。

本发明产生的有益效果是：本发明基于参考模型的加速度域残差分析，直接采用自相关序列作为基本元素来构造对称Toeplitz协方差阵，可以有效避免由轨道先验稠密阵经误差传播定律计算加速度值协方差阵时高维大型矩阵的复杂运算过程，减少定权耗时，从而提高重力场模型解算效率。相较于经典方法，改进算法不仅操作简洁高效，而且可以整体上保证全球重力场模型反演计算效率稳定可靠。

附图说明

下面将结合附图及实施例对本发明作进一步说明，附图中：

图1为本发明实施例提供的一种全球重力场模型反演方法流程图；

图2为本发明实施例提供的定权策略。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅用以解释本发明，并不用于限定本发明。

步骤1：观测数据的采集与预处理精化。

数据源包括CHAMP、GRACE、GOCE、Swarm等专用低轨重力卫星数据，采集的数据内容包括几何法轨道观测值、简化动力学轨道观测值和星载加速度计观测值。重力卫星观测数据的预处理包括粗差探测和时变重力场改正两部分：采用数理统计方法对几何法轨道进行粗差探测，超出阈值部分直接用对应的简化动力学轨道填充；采用IERS提供的背景模型，如三体引力、固体潮、海潮、极潮等分离出时变重力场信号，获取高精度的静态重力场模型。

需要说明的是，发布的轨道数据是在地固系下给定，而观测方程在惯性系下建立，所以需要进行坐标变换，在实际实施时直接采用基于春分点的经典坐标转换方法，即依次考虑岁差、章动、地球自转、极移影响。

步骤2：数值微分计算。

载体运动加速度矢量基于几何法轨道观测值时确定需要通过微分计算，而数值微分会导致高频误差放大效应，本发明采用移动中心窗口多项式微分法进行微分计算，并在此基础上基于“降采样”思想设计扩展微分滤波器，用以抑制高频噪声的放大影响；然后通过步骤1的时变重力场改正经过摄动力建模计算各项摄动加速度，同时基于星载加速度计确定非保守力；最后将载体运动加速度减去除地球引力外的摄动加速度(包括保守力和非保守力)，得到纯地球引力加速度观测值。

需要说明的是，移动中心窗口多项式微分法是以计算点为中心两侧各取等点构成窗口函数结合多项式对窗口内的观测点进行最小二乘拟合通过线性组合滤波核函数与观测值来计算中心值；随后移动窗口，接着计算该点的下一点，依次类推，组成平滑后的波谱。为了抑制高频误差影响，在移动窗口多项式微分算法基础上，基于“降采样”思想设计扩展微分滤波器，这样一方面有利于去噪，用以提高信噪比；另一方面有效地保留信号波谱中的有用信息，用以提高载体运动加速度观测值精度。

步骤3：定权策略。

基于参考模型的残差序列实现定权，其原因在于，当轨道精度较差时，误差传播后的加速度协方差阵信息较差，导致估算的位系数偏差大，反演地球重力场模型分辨率差，同时需进行高维大型矩阵运算，大大增加了计算耗时。改进算法直接利用参考重力场模型(比如EGM96、EGM2008模型等国际上已发布的地球重力场模型)通过残差分析构造观测值协方差阵，可以减少繁琐运算，从而提高计算效率。则由步骤2计算的载体运动加速度观测值扣除各项摄动力，同时对参考模型进行正演得到加速度模拟值，其推导过程如下：

地固系下引力位的球谐级数展开式：

重力加速度模拟值IRF:

局部轨道坐标系下引力位梯度简式：

局部轨道坐标系下加速度模拟值：

而局部轨道坐标系与地固系坐标变换：

地固系经过岁差、章动、地球自转及极移即可转换至惯性系(同步骤1)。由此经一系列坐标转换即可计算惯性系下参考模型正演的加速度模拟值。

其中Δ_i分别表示i时刻卫星在惯性系IRF的微分加速度观测矢量，地球引力矢量和观测噪声，和分别表示扣除地球引力外的保守力(主要包括日月引力、固体潮、海潮、极潮和相对论效应等)和非保守力； (r,θ,λ)分别表示地固球坐标的地心向径、地心余纬和地心经度。GM和R分别指地心引力常数和地球平均半径；n和m为球谐展开的阶和次；N为重力场模型的最大阶数；和表示完全规格化的n阶m次球谐位系数；和分别表示完全规格化的n阶m次缔合勒让德函数及其一阶导数；V_x、V_y和V_z分别表示局部轨道坐标系(LNRF)下沿轨、法向和径向方向处的加速度模拟值。Q表示加速度观测值的协方差阵，{σ_i}表示参考重力场模型正演的加速度模拟值与轨道经数值微分扣除各项摄动力后的加速度观测值作差所得到的残差序列。

本步骤中，改进的定权策略采用残差分析方法确定观测值协方差阵(参见图2)：首先基于参考重力场模型进行粗差探测，将先验重力场模型正演的加速度模拟值与轨道经数值微分扣除各项摄动力后的加速度观测值作差即残差序列{σ_i}，数理统计后设置扰动重力阈值；对于超出阈值的历元处，采用参考模型正演的加速度模拟值直接填充对应观测值，此时对应的残差值置零。然后将上述过程生成的残差序列为自相关基本元素经正规化处理构造观测值的对称Toeplitz协方差阵，即最后针对观测值有色噪声特性，对其协方差阵Q进行Cholesky分解Q＝R^TR,令F＝(R^T)^-1则重新构造新的观测方程y^*＝B^*x+Δ^*，使新的观测值y^*＝Fy满足白噪声(Q_y*＝E)，继而通过经典最小二乘准则解算位系数来提高计算效率。

步骤4：参数估计。

针对海量卫星观测数据，实际处理时将整个观测数据分割成若干短弧段分配线程，通过各个线程协同分别计算各弧段对应设计矩阵B_i和观测量y_i，组成若干观测方程y_i＝B_ix+Δ_i，经步骤3形成新的观测方程据此建立各弧段法方程组基于并行算法采用多线程同时计算各弧段法方程阵和继而对各弧段的法方程阵进行叠加 以减少矩阵的存储空间，缩减计算耗时。顾及位系数以次为主的顺序排列方式具有块对角占优的特性，这样有利于矩阵求逆和采用快速算法，同时为了减少计算误差，采用上述参考模型进行移去-恢复技术，利用最小二乘平差最终解算得到地球重力场模型位系数。

需要说明的是，移去-恢复技术是先将观测值中减去已知重力场信息，从而得到观测值的改正量，以此计算位系数的改正量，再加上参考重力场模型位系数即可得求解的重力场模型参数。

最后说明的是，以上优选实施例仅用以说明本发明的技术方案而非限制，尽管通过上述优选实施例已经对本发明进行了详细的描述，但本领域技术人员应当理解，可以在形式上和细节上对其做出各种各样的改变，而不偏离本发明权利要求书所限定的范围。

Claims (5)

1.一种全球重力场模型反演方法，其特征在于，包括以下步骤：

步骤1：采集低轨重力卫星数据并进行预处理；

步骤2：数值微分计算

采用移动中心窗口多项式微分法进行微分计算，并在此基础上扩展微分滤波器，抑制高频噪声的放大影响；根据步骤1的预处理结果和摄动力建模计算各项摄动加速度，同时基于星载加速度计确定非保守力；将载体运动加速度减去除地球引力外的保守力摄动加速度和非保守力摄动加速度，得到纯地球引力加速度观测值；

步骤3：观测值定权策略

引入参考重力场模型仿真模拟重力加速度，对其进行坐标变换旋转到惯性系下的重力加速度模拟值；将纯地球引力加速度观测值与重力加速度模拟值作粗差探测进行数理统计确定阈值，进而剔除超出阈值的信号；利用残差分析确定纯地球引力加速度观测值自相关序列，并将其视作基本元素，以此构造观测值自协方差阵；针对观测值有色噪声特性，对其观测值协方差阵Q进行Cholesky分解Q＝U^TU,令F＝(U^T)^-1则重新构造新的观测方程y^*＝B^*x+Δ^*，使新的观测值y^*＝Fy满足白噪声特性Q_y*＝E，继而通过经典最小二乘准则解算位系数来提高计算效率；其中U和F分别是观测值协方差阵Q的乔列斯基分解矩阵；y*、B*和Δ^*分别表示新的观测值、设计阵和噪声；Q_y*表示新观测值的协方差阵，E为单位阵；

步骤4：参数估计

针对海量卫星观测数据，实际处理时将整个观测数据分割成若干短弧段分配线程，通过各个线程协同分别计算各弧段对应设计矩阵B_i和观测量y_i，组成若干观测方程y_i＝B_ix+Δ_i，经步骤3形成新的观测方程据此建立各弧段法方程组基于并行算法采用多线程同时计算各弧段法方程阵和继而对各弧段的法方程阵进行叠加以减少矩阵的存储空间，缩减计算耗时；顾及位系数以次为主的顺序排列方式具有块对角占优的特性，采用所述参考重力场模型进行移去-恢复技术，利用最小二乘平差最终解算得到地球重力场模型位系数。

2.根据权利要求1所述的全球重力场模型反演方法，其特征在于，采集的低轨重力卫星数据包括几何法轨道观测值、简化动力学轨道观测值和星载加速度计观测值。

3.根据权利要求1所述的全球重力场模型反演方法，其特征在于，所述预处理包括粗差探测和时变重力场改正两部分：采用数理统计方法对几何法轨道进行粗差探测，超出阈值部分直接用对应的简化动力学轨道填充；采用IERS提供的背景模型分离出低轨重力卫星数据中的时变重力场信号，获取高精度的静态重力场模型。

4.根据权利要求1所述的全球重力场模型反演方法，其特征在于，对步骤3中的观测值定权的计算过程进行优化：

重力加速度模拟值IRF:

局部轨道坐标系下加速度模拟值：

加速度观测值的协方差阵：

其中▽V_i(r,θ,λ)、Δ_i分别表示i时刻卫星在惯性系IRF的微分加速度观测矢量、地球引力矢量和观测噪声，和分别表示扣除地球引力外的保守力和非保守力；r,θ,λ分别表示地固球坐标的地心向径、地心余纬和地心经度；GM和R分别指地心引力常数和地球平均半径；n和m为球谐展开的阶和次；N为重力场模型的最大阶数；和表示完全规格化的n阶m次球谐位系数；和分别表示完全规格化的n阶m次缔合勒让德函数及其一阶导数；V_x、V_y和V_z分别表示局部轨道坐标系LNRF下沿轨、法向和径向方向处的加速度模拟值；{σ_i}表示参考重力场模型正演的加速度模拟值与轨道经数值微分扣除各项摄动力后的加速度观测值作差所得到的残差序列；

具体包括：

首先基于参考重力场模型进行粗差探测，将先验重力场模型正演的重力加速度模拟值与轨道经数值微分扣除各项摄动力后的加速度观测值作差得到残差序列{σ_i}，数理统计后设置扰动重力阈值；对于超出阈值的历元处，采用参考重力场模型正演的加速度模拟值直接填充对应观测值，此时对应的残差值置零；然后以将残差序列为自相关基本元素，经正规化处理构造观测值的对称Toeplitz协方差阵。

5.根据权利要求1所述的全球重力场模型反演方法，其特征在于，移动中心窗口多项式微分法是以计算点为中心两侧各取等点构成窗口函数结合多项式对窗口内的观测点进行最小二乘拟合通过线性组合滤波核函数与观测值来计算中心值；随后移动窗口，接着计算该点的下一点，依次类推，组成平滑后的波谱。