CN110532585A - 快速解算GOCE卫星重力场模型的Torus方法和系统 - Google Patents

Abstract

本申请实施例提供一种快速解算GOCE卫星重力场模型的Torus方法和系统，其中方法包括：获取GOCE卫星引力梯度观测值，将所述观测值归算至名义轨道，并进行内插，生成Torus圆环面上的格网观测值；根据先验重力场信息获取伪观测值，利用快速傅里叶变换算法对所述伪观测值进行处理，生成Torus圆环面上的格网模拟值；根据所述格网观测值和所述格网模拟值的差异确定球谐系数改正量；根据所述球谐系数改正量确定最终球谐系数，生成最终球谐模型。本申请实施例的快速解算GOCE卫星重力场模型的Torus方法和系统，在解算卫星重力场模型时，能兼顾精度和计算难度，提高计算效率。

Description

快速解算GOCE卫星重力场模型的Torus方法和系统

技术领域

本申请涉及卫星重力场技术领域，尤其涉及一种快速解算GOCE卫星重力场模型的Torus方法和系统。

背景技术

地球重力场作为地球的一个基本物理场，它反映了地球物质空间分布、运动和变化，地球重力场及其精细结构的确定是现代大地测量的主要科学目标之一，对大地测量学、固体地球物理、海洋学、冰川学等方面的研究具有重要的科学意义。确定地球重力场，主要是建立全球重力位模型，确定区域或全球高分辨率、高精度的大地水准面模型。

目前，官方发布的GOCE卫星重力场模型有GOSG01S、IGGT_R1、 IfE_GOCE05s、GO_CONS_GCF_2_SPW_R5、NULP-02s、ITU_GGC16、GGM05G、 GO_CONS_GCF_2_TIM_R5、GO_CONS_GCF_2_DIR_R5、JYY_GOCE04S、 GOGRA04S、JYY_GOCE02S、ITG-Goce02等，其解算方法主要有直接法(direct， DIR)、空域法(space-wise，SPW)和时域法(time-wise，TIM)。

由于直接法没有任何近似和假设，所以精度较高，但对计算机要求也较高。空域法通过归算和内插大大减少了观测量，且无需考虑数据中断等问题，但同时也引入了误差。时域法不需要任何其他先验重力场信息，所以可实现“纯”GOCE 卫星重力场模型的确定，但计算效率难以提高。

发明内容

有鉴于此，本申请的目的在于提出一种快速解算GOCE卫星重力场模型的 Torus方法和系统，来解决现有卫星重力场模型的解算方法不能兼顾精度和计算难度的技术问题。

基于上述目的，在本申请的第一个方面，提出了一种快速解算GOCE卫星重力场模型的Torus方法，包括：

获取GOCE卫星引力梯度观测值，将所述观测值归算至名义轨道，并进行内插，生成Torus圆环面上的格网观测值；

根据先验重力场信息获取伪观测值，利用快速傅里叶变换算法对所述伪观测值进行处理，生成Torus圆环面上的格网模拟值；

根据所述格网观测值和所述格网模拟值的差异确定球谐系数改正量；

根据所述球谐系数改正量确定最终球谐系数，生成最终球谐模型。

在一些实施例中，所述确定球谐系数改正量，包括：确定公式

中为集总系数，又称伪观测值，利用二维快速傅里叶变换得到，是转换系数和球谐系数α_nm、β_nm的线性组合：

其中，加和下标n＝n_min[2]，n_min＝max(|k|,m)+δ，δ取值取决于k和max(|k|,m)的奇偶性关系，奇偶性相同时，δ为0，不同时δ为1。

在一些实施例中，公式中的u和Λ为定义在两个不同方向的圆周，其变化范围均为[0,2π)，这两个圆周可形成一个封闭的圆环。

在一些实施例中，所述格网观测值分布在所述封闭的圆环上。

在一些实施例中，所述格网观测值利用归算后的卫星引力梯度观测值内插得到，内插方法为Shepard方法。

在一些实施例中，所述归算后的卫星引力梯度观测值利用泰勒级数展开式得到。

在一些实施例中，所述名义轨道轨道半径与长半轴相等，没有近地点，且真近点角和平近点角相等，地球的扁率会引起以下三个轨道根数的长周期摄动；

其中，根据开普勒第三定律，卫星平均运动角速度n与轨道半径r 满足n²r³＝μ；其中，

其中，r₀、I₀均为常数，分别是名义轨道的半径和倾角。

在一些实施例中，球谐系数α_nm、β_nm满足以下约束条件：

其中，为完全规格化球谐系数。

基于上述目的，在本申请的第二个方面，还提出了一种快速解算GOCE卫星重力场模型的Torus系统，包括：

观测值获取模块，用于获取GOCE卫星引力梯度观测值，将所述观测值归算至名义轨道，并进行内插，生成Torus圆环面上的格网观测值；

模拟值获取模块，用于根据先验重力场信息获取伪观测值，利用快速傅里叶变换算法对所述伪观测值进行处理，生成Torus圆环面上的格网模拟值；

球谐系数改正量确定模块，用于根据所述格网观测值和所述格网模拟值的差异确定球谐系数改正量；

模型生成模块，用于根据所述球谐系数改正量确定最终球谐系数，生成最终球谐模型。

在一些实施例中，所述球谐系数改正量确定模块具体用于：

确定公式

中的的值，

其中为集总系数，是转换系数和球谐系数α_nm、β_nm的线性组合：

其中，加和下标n＝n_min[2]，n_min＝max(|k|,m)+δ，δ取值取决于k和max(|k|,m)的奇偶性关系，奇偶性相同时，δ为0，不同时δ为1。

本申请实施例提供一种快速解算GOCE卫星重力场模型的Torus方法和系统，其中方法包括：获取GOCE卫星引力梯度观测值，将所述观测值归算至名义轨道，并进行内插，生成Torus圆环面上的格网观测值；根据先验重力场信息获取伪观测值，利用快速傅里叶变换算法对所述伪观测值进行处理，生成Torus圆环面上的格网模拟值；根据所述格网观测值和所述格网模拟值的差异确定球谐系数改正量；根据所述球谐系数改正量确定最终球谐系数，生成最终球谐模型。本申请实施例的快速解算GOCE卫星重力场模型的Torus方法和系统，在解算卫星重力场模型时，能兼顾精度和计算难度，提高计算效率。

附图说明

通过阅读参照以下附图所作的对非限制性实施例所作的详细描述，本申请的其它特征、目的和优点将会变得更明显：

图1是本申请实施例一的快速解算GOCE卫星重力场模型的Torus方法的流程图；

图2是本申请实施例二的快速解算GOCE卫星重力场模型的Torus系统的结构示意图；

图3是本申请实施例一的快速解算GOCE卫星重力场模型的Torus方法的具体流程图；

图4是不同模型相对于模型GO_TIM_R5的阶误差RMS。

具体实施方式

下面结合附图和实施例对本申请作进一步的详细说明。可以理解的是，此处所描述的具体实施例仅仅用于解释相关发明，而非对该发明的限定。另外还需要说明的是，为了便于描述，附图中仅示出了与有关发明相关的部分。

需要说明的是，在不冲突的情况下，本申请中的实施例及实施例中的特征可以相互组合。下面将参考附图并结合实施例来详细说明本申请。

如图1所示，是本申请实施例一的快速解算GOCE卫星重力场模型的Torus 方法的流程图。从图1中可以看出，本实施例的快速解算GOCE卫星重力场模型的Torus方法，可以包括以下步骤：

S101：获取GOCE卫星引力梯度观测值，将所述观测值归算至名义轨道，并进行内插，生成Torus圆环面上的格网观测值。

所述将所述观测值利用泰勒级数展开归算至名义轨道，并实现内插，使得生成的格网观测值步长相同。

S102：根据先验重力场信息获取伪观测值，利用快速傅里叶变换算法对所述伪观测值进行处理，生成Torus圆环面上的格网模拟值。

具体地，确定公式

其中为集总系数，是伪观测值,也是转换系数和球谐系数α_nm、β_nm的线性组合：

其中，加和下标n＝n_min[2]，n_min＝max(|k|,m)+δ，δ取值取决于k和max(|k|,m)的奇偶性关系，奇偶性相同时，δ为0，不同时δ为1。

所述名义轨道轨道半径与长半轴相等，没有近地点，且真近点角和平近点角相等，地球的扁率会引起以下三个轨道根数的长周期摄动。

其中，根据开普勒第三定律，卫星平均运动角速度n与轨道半径r 满足n²r³＝μ；其中，

其中，r₀、I₀均为常数，分别是名义轨道的半径和倾角。

公式中的u和Λ为定义在两个不同方向的圆周，其变化范围均为[0,2π)，这两个圆周可形成一个封闭的圆环。所述格网观测值分布在所述封闭的圆环上。

S103：根据所述格网观测值和所述格网模拟值的差异确定球谐系数改正量。

具体地，确定公式

S104：根据所述球谐系数改正量确定最终球谐系数，生成最终球谐模型。

根据卫星轨道摄动理论，局部轨道坐标系(local orbit reference famce，LORF)下引力位与卫星6个轨道根数之间的关系可表示为：

其中，a为卫星轨道长半轴；I为轨道倾角；u为升交角距，u＝ω+M，ω和M分别为近地点角距和平近点角；Λ为升交点经度，Λ＝Ω－θG，其中Ω为升交点赤经，θ_G为格林尼治恒星时角；μ为地心引力常数，即万有引力常数与地球质量的乘积；R为地球平均半径；n、m为球谐函数的阶和次，第3 个加和下标中的[2]表示k的变化步长为2。F_nmk(I)和Gnkq(e)分别为正规化的倾角函数和偏心率函数，e为偏心率，ψ_mkq为u、Λ和M的函数，ψ_mkq＝ku+mΛ+qM。球谐系数α_nm和β_nm随n－m的奇偶性取不同的值，具体关系如下：

其中，为完全规格化球谐系数。

对于e≈0的近圆轨道，上式可以简化为：

式中，r为地心向径；ψ_mk＝ku+mΛ，此时，u＝ω+f，f为真近点角。上式中级数求和的n应取至无穷大，但实际计算中常截止到某最高阶N。

经过符号代换，可写为：

式中，A_mk、B_mk为集总系数，是转换系数H^V _nmk和球谐系数(α_nm和β_nm) 的线性组合。

式中，加和下标n＝n_min[2]，n_min＝max(|k|,m)+δ，若k与max(|k|,m)奇偶性相同，则δ＝0，否则δ＝1。

公式是SA方法的基本表达式，是傅里叶级数的形式，集总系数A_mk和B_mk为傅里叶系数。SA方法要求使用名义轨道上的观测值，名义轨道具有3个特点：轨道为圆形、轨道倾角为常量、轨道受J2项影响而长周期进动，同时，为便于 FFT快速计算，要求等间隔采样的观测值。公式中变量u和Λ分别为定义在两个不同方向的圆周，其变化范围均为[0，2π)，这两个圆周可形成一个封闭的圆环，若观测值均匀分布在该圆环面上，则可以利用2D-FFT方法解算集总系数，因此2D-FFT对应的SA方法称为Tours方法。采用Tours方法由GOCE沿轨卫星引力梯度观测值确定地球重力场模型的具体流程如图3所示，其中a_i表示第i 次迭代的集总系数，K_i+1为第i次迭代后的球谐系数。

本实施例的快速解算GOCE卫星重力场模型的Torus方法，在解算卫星重力场模型时，能兼顾精度和计算难度，提高计算效率。

如图2所示，是本申请实施例二的快速解算GOCE卫星重力场模型的Torus 系统的结构示意图。本实施例的快速解算GOCE卫星重力场模型的Torus系统，包括：

观测值获取模块201，用于获取GOCE卫星引力梯度观测值，将所述观测值归算至名义轨道，并进行内插，生成Torus圆环面上的格网观测值。

模拟值获取模块202，用于根据先验重力场信息获取伪观测值，利用快速傅里叶变换算法对所述伪观测值进行处理，生成Torus圆环面上的格网模拟值。

球谐系数改正量确定模块203，用于根据所述格网观测值和所述格网模拟值的差异确定球谐系数改正量。

模型生成模块204，用于根据所述球谐系数改正量确定最终球谐系数，生成最终球谐模型。

其中，所述球谐系数改正量确定模块203具体用于：

确定公式

中的的值，

其中为集总系数，是转换系数和球谐系数α_nm、β_nm的线性组合：

其中，加和下标n＝n_min[2]，n_min＝max(|k|,m)+δ，δ取值取决于k和max(|k|,m)的奇偶性关系，奇偶性相同时，δ为0，不同时δ为1。

公式中的u和Λ为定义在两个不同方向的圆周，其变化范围均为[0,2π)，这两个圆周可形成一个封闭的圆环。所述格网观测值分布在所述封闭的圆环上。所述将所述观测值归算至名义轨道，并进行内插，包括：所述将所述观测值利用泰勒级数展开归算至名义轨道，并利用Shepard方法实现内插。具体地：先将所述观测值利用泰勒级数展开式进行归算，然后利用Shepard方法实现内插，使得生成的格网观测值步长相同。

本实施例的快速解算GOCE卫星重力场模型的Torus系统，在解算卫星重力场模型时，能兼顾精度和计算难度，提高计算效率。

利用2009年11月～2011年8月GRF下的GOCE卫星引力梯度观测值，选择 EGM2008前200阶次的系数为参考模型，采用Torus方法解算200阶次的GOCE 卫星重力场模型CASM_Torus_GOCE。

不同模型相对于模型GO_TIM_R5的阶误差RMS见图4。所有模型的阶误差 RMS均小于模型GO_TIM_R5信号的大小。其中，模型GOSG01C同样采用了2009 年11月～2011年8月的卫星引力梯度观测值(仅Vxx、Vyy、Vzz三个分量)，还联合采用了2009.11.01～2010.07.05的GOCE卫星跟踪卫星数据，解算方法为直接法，最高阶次为220，没有使用任何先验模型，且只对带谐项和近带谐项采用了正则化(Xu et al.，2017)。与GOCE_Torus0模型(71天实测数据解算的模型) 相比，130阶以后，CASM_Torus_GOCE(即图中的GOCE_Torus1)模型精度提高明显。受格网化误差和参考模型的影响，与同期模型相比，100～160阶内， CASM_Torus_GOCE模型精度略低，但明显高于EGM2008的精度。由于高阶采用了Kaula约束，所以170阶以后，CASM_Torus_GOCE模型的精度优于GOSG01C。

利用美国和中国的GPS/水准点对各重力场模型进行了外部检核，具体结果见表1和表2。在美国地区，CASM_Torus_GOCE模型与GO_TIM_R5、EGM2008 和其他同期模型的精度相当，差别仅在4mm以内。在中国地区， CASM_Torus_GOCE模型的精度与同期模型相当。

表1：地球重力场模型(截断至200阶次)确定的大地水准面高与美国GPS/ 水准的差值(单位：m)

表2：地球重力场模型(截断至200阶次)确定的大地水准面高与中国GPS/ 水准的差值(单位：m)

利用EGM2008模型200阶次以上的系数对各模型进行了补充，利用GPS/水准数据进行精度检核，结果见表3。各模型在美国地区的精度均相当。由于缺乏地面重力数据，EGM2008模型在中国的误差仍然最大，达0.240m。受EGM2008 模型的影响，CASM_Torus_GOCE模型的精度略低于同期模型的精度，但相比模型EGM2008，CASM_Torus_GOCE模型误差明显减小，为0.182m。

表3：地球重力场模型确定的大地水准面高与GPS/水准的差值(单位：m)， 200阶次以上的系数用EGM2008模型的系数进行了补充。

以上描述仅为本申请的较佳实施例以及对所运用技术原理的说明。本领域技术人员应当理解，本申请中所涉及的发明范围，并不限于上述技术特征的特定组合而成的技术方案，同时也应涵盖在不脱离上述发明构思的情况下，由上述技术特征或其等同特征进行任意组合而形成的其它技术方案。例如上述特征与本申请中公开的(但不限于)具有类似功能的技术特征进行互相替换而形成的技术方案。

Claims (9)

1.快速解算GOCE卫星重力场模型的Torus方法，其特征在于，包括：

获取GOCE卫星引力梯度观测值，将所述观测值归算至名义轨道，并进行内插，生成Torus圆环面上的格网观测值；

根据先验重力场信息获取伪观测值，利用快速傅里叶变换算法对所述伪观测值进行处理，生成Torus圆环面上的格网模拟值；

根据所述格网观测值和所述格网模拟值的差异确定球谐系数改正量；

根据所述球谐系数改正量确定最终球谐系数，生成最终球谐模型。

2.根据权利要求1所述的方法，其特征在于，所述确定球谐系数改正量，包括：确定公式

中为集总系数，又称伪观测值，利用二维快速傅里叶变换得到，是转换系数和球谐系数α_nm、β_nm的线性组合：

其中，加和下标n＝n_min[2]，n_min＝max(|k|，m)+δ，δ取值取决于k和max(|k|，m)的奇偶性关系，奇偶性相同时，δ为0，不同时δ为1。

球谐系数α_nm、β_nm满足以下约束条件：

其中，为完全规格化球谐系数。

3.根据权利要求2所述的方法，其特征在于，观测值应位于名义轨道，其轨道半径与长半轴相等，没有近地点，且真近点角和平近点角相等，地球的扁率会引起以下三个轨道根数的长周期摄动。

其中，根据开普勒第三定律，卫星平均运动角速度n与轨道半径r满足n²r³＝μ；其中，

其中，r₀、I₀均为常数，分别是名义轨道的半径和倾角。

4.根据权利要求3所述的方法，其特征在于，u和Λ为定义在两个不同方向的圆周，其变化范围均为[0，2π)，这两个圆周可形成一个封闭的圆环。

5.根据权利要求4所述的方法，其特征在于，格网观测值分布在所述封闭的圆环上。

6.根据权利要求5所述的方法，其特征在于，格网观测值利用归算后的卫星引力梯度观测值内插得到，内插方法为Shepard方法。

7.根据权利要求6所述的方法，其特征在于，所述归算后的卫星引力梯度观测值是利用泰勒级数展开式得到的。

8.快速解算GOCE卫星重力场模型的Torus系统，其特征在于，包括：

观测值获取模块，用于获取GOCE卫星引力梯度观测值，将所述观测值归算至名义轨道，并进行内插，生成Torus圆环面上的格网观测值；

模拟值获取模块，用于根据先验重力场信息获取伪观测值，利用快速傅里叶变换算法对所述伪观测值进行处理，生成Torus圆环面上的格网模拟值；

球谐系数改正量确定模块，用于根据所述格网观测值和所述格网模拟值的差异确定球谐系数改正量；

模型生成模块，用于根据所述球谐系数改正量确定最终球谐系数，生成最终球谐模型。

9.根据权利要求8所述的系统，其特征在于，所述球谐系数改正量确定模块具体用于：

确定公式

中的的值，

其中为集总系数，是转换系数和球谐系数α_nm、β_nm的线性组合：

其中，加和下标n＝n_min[2]，n_min＝max(|k|，m)+δ，δ取值取决于k和max(|k|，m)的奇偶性关系，奇偶性相同时，δ为0，不同时δ为1。