基于卡尔曼滤波重构GRACE-FO姿态数据

梁磊，闫易浩，王长青，朱紫彤，高铭，钟敏，于锦海，徐焕

1 中国科学院精密测量科学与技术创新研究院大地测量与地球动力学国家重点实验室,武汉 430077 2 地理信息工程国家重点实验室，西安 710054 3 中山大学物理与天文学院，广东珠海 519082 4 中山大学测绘科学与技术学院，广东珠海 519082 5 中国科学院大学地球与行星科学学院，北京 100049

GRACE-FO（GRACE Follow-On）是GRACE（Gravity Recovery and Climate Experiment）的后续卫星，于2018年5月22日成功发射（Kornfeld et al.,2019）.二者是同类型卫星，通过精确测量两颗卫星之间的距离反演地球时变重力场模型，研究地球质量变化（Tapley et al.,2004）.GRACE-FO相较于GRACE而言除了搭载K波段微波测距系统外，还首次实现了激光干涉测距，实现了以nm级的精度测量两颗卫星之间的距离（Goswami et al.,2021）.GRACE与GRACE-FO的数据产品分为Level-0、Level-1A、Level-1B和Level-2，其中Level-0至Level-1B称为载荷数据处理，Level-1B至Level-2称为时变重力场反演（Wu et al.,2006;Case et al.,2010;Wen et al.,2019）.

在载荷数据处理过程中，GNSS（Global Navigation Satellite System）接收机、加速度计和星间测距等数据处理与姿态数据有关，主要为：第一个是计算KBR（K-Band Ranging）天线相位中心改正和GNSS天线相位中心改正，将测量的信息改正到卫星质心上；
第二个是将加速度计（ACC，Accelerometer）测量的非保守力转换到惯性坐标系；
第三个是标定KBR天线相位中心位置和卫星质心位置；
第四个是评估激光测距系统的耦合噪声.姿态数据误差会通过这四个方面传递到重力场模型中，例如Inácio等（2015）分析GRACE星间距变率残差与KBR天线相位中心改正时，指出姿态误差占总误差的18%左右；
Horwath等（2011）通过KBR天线相位中心改正，指出姿态误差导致了星间指向存在偏差，通过对姿态误差进行建模提高了重力场模型解算精度；
Wegener等（2020）通过姿态数据评估了激光反射棱镜交点与卫星质心不重合对激光测距的影响.因此确定高精度的重力卫星姿态数据是原始载荷数据处理的重要内容之一.

GRACE与GRACE-FO均搭载了星敏感器和惯性测量单元测量卫星的姿态，但二者又有较大区别.对于GRACE而言，每颗卫星上均搭载了两颗星敏感器和一个三轴的惯性测量单元.星敏感器固定安装于卫星左右两个面上，三轴的惯性测量单元固定安装于卫星的内部.然而，在发射之初GRACE-A星的惯性测量单元一个轴出现故障，未有冗余观测量，GRACE-B星的惯性测量单元仅在备份模式下开启（Bandikova,2015），所以两颗卫星均未有IMU的测量数据，因此GRACE卫星姿态主要由星敏感器提供.受日月光照的影响GRACE只有一颗星敏感器或者无星敏感器测量卫星姿态，从而导致卫星姿态精度降低或出现大量间断.借鉴于GRACE上的经验，GRACE-FO在设计之初采用了三颗星敏感器和一个四轴光纤陀螺仪设计方案.相较于GRACE而言，三颗星敏感器设计方案减少了姿态数据间断，并且较多的时间存在两颗星敏感器测量卫星的姿态，部分时间有三颗星敏感器测量卫星姿态；
四轴光纤陀螺仪具有更好的测量性能，固定安装于卫星的内部，同时整个科学任务阶段均可观测（Kornfeld et al.,2019）.

融合多个载荷数据，提高最终数据产品处理精度，是目前数据处理常用的方法.对于GRACE而言，除了星敏感器测量卫星姿态数据外，加速度计可以同步测量卫星的角加速度数据，融合星敏感器数据和角加速度数据，可以提高姿态数据精度，弥补星敏感器数据的缺失，减少姿态数据间断.例如Klinger（2018）在反演重力场模型时，首先基于最小二乘配置法融合了这两类数据，构建了高精度的姿态数据产品；
Goswami等（2018）组合了星敏感器与角加速度计数据，提高了姿态数据产品的精度；
官方机构JPL（Jet Propulsion Laboratory）基于Kalman滤波融合了这两类数据，发布了RL03版本的姿态数据（Harvey and Sakumura,2019）.GOCE卫星是另一颗用于地球重力场研究的卫星，其上搭载了引力梯度仪，测量地球引力梯度数据，可同步获得测量频段内的角加速度.由于在处理梯度数据时，需要改正的离心力项与卫星的角速度是非线性关系（Stummer et al.,2011;Pail,2005；
郭泽华等，2021；
赵宇鹏等，2021），这就需要高精度的姿态信息，由此需要融合星敏感器和梯度仪测量的角加速度数据.例如Cesare等（2008）在时域上通过Kalman滤波融合了这两类数据；
Stummer等（2011）和郭泽华等（2021）在频率域上，基于Wiener滤波融合了这两类数据.对于GRACE-FO而言，星敏感器和IMU是测量卫星姿态的主要载荷，其中星敏感器对姿态的低频部分敏感，IMU对姿态的高频部分敏感，融合这两类数据可以获得高精度的姿态数据.在时域上进行数据融合具有一定的便利性，可以同步评估IMU的漂移参数，弥补微小的间断，因此本文选择在时域下进行星敏感器数据与IMU数据融合.

关于姿态Kalman滤波，Lefferts等（1982）基于状态变量（四元数和陀螺仪的漂移参数）的一阶导形式构建了姿态Kalman滤波算法，此时在计算状态转移矩阵和状态转移方程时，其形式了变成了Riccati方程；
Trawny和Roumeliotis（2005）与Lefferts方法类似，只是推导过程将Riccati方程变成级数形式.目前关于这种方法尚未有研究应用于实际的GRACE-FO数据处理中.Harvey和Sakumura（2019）将角速度近似处理转换成四元数形式微小变化量作为Kalman滤波的状态转移矩阵，并基于此方法解算了GRACE和GRACE-FO的姿态数据，这种近似处理过程理论上具有一定的瑕疵.本文将以四元数和陀螺仪的漂移参数作为状态变量，直接构建另一种姿态Kalman滤波融合算法，处理GRACE-FO的姿态数据.鉴于此，本文开展了GRACE-FO高精度的姿态数据处理的研究.在本文的第1节是介绍星敏感器处理算法、IMU处理算法和姿态Kalman滤波算法；
第2节是分析姿态数据处理结果；
第3节是评估姿态数据产品对时变重力场的影响；
第4节是结论和总结.

GRACE-FO搭载的三颗星敏感器和IMU记录的时标是OBC（Onboard Computer）时标，首先需要利用TIM1B和CLK1B数据将OBC时标改正到GPS（Global Position System）时标上.其次对2 Hz采样的多星敏感器数据进行组合并降采样至1 Hz，然后通过姿态Kalman滤波与8 Hz采样的IMU数据进行融合，输出1 Hz姿态数据产品.下面分别介绍星敏感器数据处理、IMU数据处理和姿态Kalman滤波算法，其中星敏感器数据处理和IMU数据处理中的详细细节可参考Yang等（2022）.

1.1 星敏感器数据处理

星敏感器测量的噪声具有各向异性的特点，即绕垂直于星敏感器视轴方向旋转的精度要比绕视轴方向旋转的精度大约高8～10倍（Bandikova and Flury,2014;Stanton,2000）.由于在重力场反演过程中，需要的是SRF （Science Reference Frame）坐标系下的卫星姿态数据.为此，需要将星敏感器坐标系下测量的四元数转换到SRF坐标系中，其中星敏感器坐标系与SRF的之间的转换矩阵（即星敏感器安装矩阵）记录于QSA1B或者SOE文件中.转换过程可看成是一个线性组合的形式，绕星敏感器视轴旋转的精度低会传递到绕科学坐标系Y轴和Z轴旋转中，导致绕这两个轴旋转的精度降低.融合多个星敏感器数据可以降低这一影响.目前，多星敏感器融合主要有两种方法.第一种是根据星敏感器测量的噪声构建权阵，通过加权的形式融合多星敏感器数据（Romans，2003）.第二种是利用星敏感器的视轴构建一个共同参考框架融合多星敏感器数据（Mandea et al.,2010）.这两种方法均已用于多种卫星，例如CHAMP（Challenging Minisatellite Payload）、GRACE、GRACE-FO和GOCE（Gravity field and steady-state Ocean Circulation Explorer）（Bandikova and Flury,2014;Siemes,2011;Stummer et al.,2011）.对于这两种方法而言，融合后的姿态数据精度基本一致（Bandikova,2015）.本文选择使用第一种方法，融合多个星敏感器的测量值，相关算法具体参考Romans（2003）.

四元数运算有别于实数运算，需要采用合适的插值算法保证插值后的姿态具有很好的光滑性.本文采用单位四元素插值样条曲线对四元数进行插值，保证插值后的卫星姿态具有二阶导连续特点（邢燕等，2017），具体公式为

1.2 惯性测量单元数据处理

GRACE-FO的惯性测量单元是由4个四面体形式的光纤陀螺组成，每个光纤陀螺在其自身的陀螺仪框架下记录滤波角（filtered angle），采样率为8 Hz（Wen et al.,2019）.由于IMU记录的每个轴采样时刻并不一致，需要对每个轴的时间进行统一.IMU记录的滤波角的范围是从-5758°到5758°，当记录到5758°时自动调整到-5758°开始，或与之相反.然而这导致对滤波角进行数值差分计算角速度时，存在异常情况，需要对异常情况进行处理.对于数值差分计算，为了保证数值差分计算的精度，本文选择使用Diebel（2006）的方法，具体形式如下：

（1）

由于IMU测量存在冗余观测，需要对冗余观测量进行处理，根据Jafari（2015）指出同时融合冗余观测量，可以提高角速度的精度，处理方法如下：

ω=（HTH）-1HTm，

（2）

其中，H4×3是一个转换矩阵，表示的是4个光纤陀螺仪测量轴在IMU坐标系下的矢量，记录于GRACE-FO Level-1A数据使用手册中（Wen et al.,2019），m是陀螺仪的角速度测量值，ω是IMU坐标系下角速度.

1.3 姿态Kalman滤波方法

四元数运动方程可写成如下形式（秦永元，2014）：

（3）

其中q是四元数矢量，其运算法则按照附录A的形式计算，ω=[ωx,ωy,ωz]，对四元数微分方程组（3）按照角增量的形式进行求解，省略中间求解过程，可得：

q（t+1）=Φ（ωΔt）q（t），

（4）

其中

（5）

（6）

陀螺仪的噪声模型可以写成如下的形式（Lefferts et al.,1982）

（7）

E[εω]=0，

（8）

E[εω（t1）,εω（t2）]=Q1（t）δ（t1-t2），

（9）

（10）

（11）

δ是克莱罗符号，Q1与Q2是相应的协方差.将公式（7）代入到公式（4）中，可以写成

（12）

（13）

将公式（12）和（13）代入到公式（4）中，有

X（t+1）=f（t,X（t））+G（t）W（t），

（15）

那么f（t,X（t））可以写成如下形式：

（16）

根据附录公式（A3），G（t）W（t）可以写成

（17）

对公式（17）进行泰勒展开，保留线性形式有

（18）

（19）

四元数的观测方程来源于星敏感器，有

Z（t）=HX（t）+vk，

（20）

其中H=[I，0]，状态转移矩阵可以写成

（21）

基于扩展Kalman滤波，姿态Kalman滤波算法过程按照以下五步进行：

第一步：状态预测方程

（22）

第二步：协方差矩阵预测

P（t+1|t）=Φ（t+1|t）P（t|t）ΦT（t+1|t）+Q（t+1）；

（23）

第三步：计算姿态Kalman滤波增益

K（t+1）=P（t+1|t）HT（t+1）[H（t+1）P（t+1|t）×HT（t+1）+R（t+1）]-1

（24）

第四步：状态更新

（25）

第五步：协方差更新

P（t+1）=[I-K（t+1）H（t+1）]P（t+1|t）.

（26）

为了对处理后的姿态数据的精度进行评估，本文将与JPL发布的GRACE-FO Level-1B姿态产品进行对比分析，即SCA1B数据.对比分析的方法主要为星间指向（inter-satellite pointing）和卫星角速度矢量.星间指向为KBR天线相位中心矢量与LOS（Light of Sight）矢量的偏差并用角来表示，其中LOS为两颗卫星质心之间的连线方向上的矢量通过轨道数据计算（Bandikova et al.,2012;Bandikova，2015）.星间指向主要是用于卫星在轨运行期间，监测两颗卫星的对准情况，是判断是否需要进行在轨控制的依据.由于直接计算惯性坐标系下卫星的姿态角，会受到较大值的影响难以展示相关细节，因此为了便于直观展示解算的最终结果，本文采用星间指向进行对比分析.卫星角速度矢量用于评估融合后的姿态数据在中高频部分是否完全融合了IMU信息.

星间指向计算过程如下：

（27）

（28）

求解的φ、θ和ψ即为星间指向上的偏差角.

由于官方机构并未公布Level-0数据，仅公布了Level-1A和Level-1B数据，因此本文将从Level-1A数据出发构建GRACE-FO的姿态数据.本文选择以2019年1月1日的数据为例，对星敏感器数据处理以及与IMU数据进行姿态Kalman滤波融合，并对滤波后的结果进行对比分析.

2.1 星间指向

首先选择2019年1月1日中三颗星敏感器同时存在观测数据的时间段进行分析.单个星敏感器解算的卫星姿态结果如图1所示，其中SCA1、SCA2和SCA3表示星敏感器的编号（Kornfeld et al.,2019）.从星间指向的Yaw方向可见，即图1中第一行，三个星敏感器的时间序列变化并不一致但基本维持在±1×10-3rad之间，从功率谱上看三个星敏感器在高频部分的精度基本一致.从星间指向的Pitch方向可见，即图1中第二行，三个星敏感器的时间序列变化基本维持在±1×10-3rad之间，并且SCA2和SCA3基本一致且有较多毛刺，而SCA1变化更加平滑.从功率谱上看SCA2和SCA3在高频部分的精度基本一致，而SCA1具有更低的高频噪声.这主要是由于SCA1安装在卫星的天顶方向，星敏感器视轴方向刚好与科学坐标系的Z轴平行，而SCA2和SCA3安装在卫星左右两个面方向上，从星敏感器坐标系转换到SRF坐标系过程中，受到绕视轴旋转不精确的影响传递到Pitch方向上.从星间指向的Roll方向可见，即图中第三行，三个星敏感器的时间序列变化趋势基本一致，从功率谱上看三个星敏感器在高频部分精度处于同一水平.

融合多个星敏感器计算的星间指向的结果如图2所示，其中SCA12、SCA13、SCA23和SCA123表示融合不同的星敏感器，例如SCA12表示融合SCA1和SCA2两个星敏感器，SCA123表示融合SCA1、SCA2和SCA3三个星敏感器.从星间指向的Yaw、Pitch和Roll方向的时间序列可见，三个时间序列的变化趋势基本一致.从功率谱上看，与图1中的功率谱对比可见，融合多个星敏感器可以提高高频部分的精度，并且融合三颗星敏感器解算的姿态的高频部分的精度要略优于融合两颗星敏感器的结果.此外，从星间指向的Pitch方向可发现一个有趣的现象，融合不同星敏感器计算的Pitch角存在一个偏差，最大相差约3.0×10-4rad.我们检查了相应的计算流程，发现QSA1B中记录的星敏感器坐标系与SRF坐标系的转换四元数仅提供了地面标定值，未提供在轨标定值.对比GRACE的SOE文件可见，在GRACE入轨后进行了数次在轨标定，提供了相应的星敏感器坐标系与SRF坐标系转换四元数.目前，由于缺少相关的遥测数据，无法验证这样的偏差是否来源于QSA1B，因此在本文的后续处理过程中，暂不考虑这样的微小偏差影响.

融合星敏感器数据与IMU数据解算卫星姿态数据，其中首先对星敏感器数据进行融合并降采样至1 Hz，然后再与8 Hz的IMU数据进行姿态Kalman滤波融合，在姿态Kalman滤波过程中仅保留整数秒的姿态数据.融合后的卫星姿态数据结果如图2所示，其中APM表示本文所处理的姿态数据结果，JPL表示官方机构解算的结果.从星间指向的时间序列上看，Yaw、Pitch和Roll的时间序列趋势变化基本一致，仅在Pitch方向上存在一个小的偏差，约为0.6 mrad.从功率谱看上基本一致，仅在截尾部分（即0.2～0.5 Hz处）JPL略低，这里将在后面进一步分析.

图1 单个星敏感器解算的卫星姿态（a） 左列是星间指向的时间序列，（b） 右列是对应的功率谱，其中第一至三行分别表示星间指向的Yaw、Pitch和Roll三个方向.Fig.1 Satellite attitude calculated by a single star camera（a） The left plots are the time series of inter-satellite pointing,and （b） the right plots are the corresponding amplitude spectral density,where the first to third rows represent the Yaw,Pitch,and Roll directions of the inter-satellite pointing,respectively.

图2 融合多个星敏感器解算的卫星姿态（a） 左列是星间指向的时间序列，（b） 右列是对应的功率谱，其中第一至三行分别表示星间指向的Yaw、Pitch和Roll三个方向.Fig.2 Satellite attitude calculated by fusing multiple star cameras（a） The left plots are the time series of inter-satellite pointing,and （b） the right plots are the corresponding amplitude spectral density,where the first to third rows represent the Yaw,Pitch,and Roll directions of the inter-satellite pointing,respectively.

图3 星间指向的分析（a） 左列表示相应星间指向的时间序列，（b） 右列表示相应的功率谱，其中第一至三行分别表示星间指向的Yaw、Pitch和Roll三个方向.Fig.3 Analysis of inter-satellite pointing（a） The left plots are the time series of inter-satellite pointing,and （b） the right plots are the corresponding amplitude spectral density,where the first to third rows represent the Yaw,Pitch,and Roll directions of the inter-satellite pointing,respectively.

图4 卫星角速度矢量的分析（a） 左列表示时间序列，（b） 右列表示相应的功率谱，其中第一至三行分别表示ωx、ωy和ωz三个方向.Fig.4 Analysis of satellite angular velocity vector（a） The left plots are the time series of inter-satellite pointing,and （b） the right plots are the corresponding amplitude spectral density,where the first to third rows represent the ωx,ωy,and ωz.

图5 惯性坐标系下非保守力的比较（a） 左列表示APM与JPL姿态产品计算的非保守力差异的时间序列，（b） 右列是相应的功率，其中第一行至第三行分别表示惯性坐标系下X，Y和Z方向.Fig.5 Comparison of non-conservative forces in an inertial coordinate system（a） The left plots represent the time series of non-conservative force differences calculated by the APM and JPL attitude data,（b） the right plots represent the corresponding amplitude spectral density,where the first to third rows represent the X,Y,and Z directions of the inertial coordinate system.

图6 KBR天线相位中心改正变率（a） 时间序列;（b） 相对应的功率谱.Fig.6 KBR antenna phase correction rate（a） The time-series;（b） The corresponding amplitude spectral density.

图7 时变重力场模型的阶方差（a） 2018年6月;（b） 2019年3月.Fig.7 The degree variance of the time-variable gravity field model（a） The June 2018;（b） The March 2019.

图8 质量变化等效水高（a） 2018年6月；
（b） 2019年3月.Fig.8 The equivalent water height for mass change（a） June 2018；
（b） March 2019.

2.2 卫星角速度

利用公式（3）计算卫星的角速度验证进行比较，结果如图4所示.图4中APM、JPL和IMU1B分别表示本文解算的结果、官方机构结果和IMU测量的角速度.APM-IMU1B、JPL-IMU1B分别表示APM和 JPL计算的角速度与IMU1B测量值的差异.为了更好的展示细节，在图4的左列中仅展示了2019年1月1日前10000秒的时间序列，其中相对应的功率谱则使用完整的一天数据进行计算.从图4的时间序列可见，{ωx,ωy,ωz}变化趋势基本一致，APM与IMU1B的差异在三个轴的标准差分别为1.89×10-7rad·s-1，1.56×10-7rad·s-1，1.63×10-7rad·s-1，JPL与IMU1B的差异在三个轴的标准差分别为6.09×10-7rad·s-1，3.62×10-7rad·s-1，4.01×10-7rad·s-1，可见APM相较于JPL而言精度至少提高了3倍，具有更高的精度.从功率谱上看，APM、JPL和IMU1B变化基本一致，但是从与IMU1B产品的差异上看，APM在0.004～0.2 Hz上具有更低的噪声，特别在0.01～0.1 Hz相较于JPL而言噪声水平低一个量级左右，即图中青色虚线低于品红色虚线.换句话说本文所解算的姿态数据充分融合IMU1B数据，具有更低的噪声水平.在截尾部分（即0.2～0.5 Hz处）可见，APM与JPL都有明显的下降趋势，这与滤波所选取的截止频率有关.

正如引言所述，姿态数据主要通过非保守力转换和KBR天线相位中心改正影响重力场模型解算精度.因此，本节将对比分析APM和JPL姿态数据对非保守力转换以及KBR天线相位中心改正的影响，并进一步反演时变重力场模型分析姿态数据对时变重力场模型精度的影响.时变重力场反演的过程及等效水高计算过程本文参考Liang等（2021）、梁磊等（2019）和Zhou等（2018）.

计算惯性坐标系下非保守力的结果如图5所示.由于只需要评估APM与JPL姿态数据计算惯性坐标系下非保守力的差异，因此在图5左列中的时间序列只展示了转换后的非保守力的差异，在图5的右列中展示了JPL与APM转换后的非保守力及其差异的功率谱.对图5中的时间序列进行统计，结果如表1所示.从图5左列的时间序列和表1可见，通过APM和JPL姿态数据转换后的非保守力差异的标准差在10-10m·s-2左右，非保守力差异的均值在X和Y轴的方向为10-11m·s-2左右，在Z轴的方向均值差异为4.61×10-10m·s-2，也就是说在6 h积分弧长上，非保守力的差异累积不超过2 cm，低于定轨精度.从图5右列的功率谱上看高频部分的差异主要受姿态数据的影响，低频部分主要是集中于1CPR和2CPR.对于1CPR和2CPR的低频信号，在重力场反演过程中可看成是低频误差（Liang et al.,2021;Kim,2000;Zhao et al.,2011;Zhou et al.,2018,2019），通过低频误差处理或者加速度标定进行吸收，并不影响重力场模型解算精度.

表1 转后非保守力差异的标准差和均值Table 1 Standard deviation and mean of the non-conservative force difference

由于在反演重力场过程中，主要使用的是星间距变率作为观测值，因此这里仅计算KBR天线相位中心改正变率进行分析，计算方法参考Bandikova（2015），结果如图6所示.从图中可见，KBR天线相位中心改正的时间变率量级在10-9m·s-1左右，远低于KBR测距精度1 μm·s-1的要求，因此KBR天线相位中心改正的时间变率对重力场反演影响很小.此外，可见GRACE-FO对两颗卫星的对准情况要求要更为严格.

选取2018年6月和2019年3月数据反演时变重力场模型进行比较，其中反演重力场过程中只考虑姿态数据的差异，其他数据则使用JPL Level-1B产品.反演时变重力场模型的阶方差和等效水高结果如图7—8所示.从图7可见，反演时变重力场模型的阶方差基本一致.图8表示的是两个模型计算的等效水高差异，从图中可见在全球范围内利用APM与JPL姿态数据计算的等效水高基本一致，仅在GRACE-FO任务开始阶段（2018年6月），在质量变化较大的亚马逊流域上有略微差异，约占最大信号的5%左右，进入到任务稳定时期（2019年3月），计算的全球等效水高基本一致.

本文研究了GRACE-FO姿态数据Level-1A至Level-1B的处理算法及其对时变重力场的影响.首先根据GRACE-FO载荷设计特点，从四元数运动微分方程出发，以四元数和陀螺仪漂移参数为状态变量，构建了一种新的姿态Kalman滤波融合算法.然后以GRACE-FO数据为基础，分析了单个星敏感及多星敏感器组合对卫星姿态解算的影响，可见融合多星敏感器数据可以抑制高频部分的噪声.另外，星敏感器数据之间的偏差可能来源于官方机构未提供在轨标定的星敏感器安装矩阵，需要进一步利用遥测数据进行分析验证.最后基于本文推导的姿态Kalman滤波算法，处理实际的GRACE-FO Level-1A数据，从星间指向和角速度矢量可见，本文解算的姿态数据充分发挥了IMU载荷的性能，具有更低的噪声水平，优于JPL解算的结果.

对时变重力场模型精度的影响，本文首先分析了姿态数据对转换后的非保守力和KBR天线相位中心改正的影响.从转换后的非保守力看，本文解算的姿态数据与JPL的姿态数据引起的差异在10-10m·s-2左右，从KBR天线相位中心改正上看，本文解算的姿态数据与JPL的姿态数据计算的KBR天线相位中心改正变率量级在10-9m·s-1，远低于KBR星间变率1 μm·s-1的精度要求.进一步解算了2018年6月和2019年3月的时变重力场进行比较，从阶方差和等效水高差异上看，本文解算的姿态数据与JPL的姿态数据对时变重力场模型的精度影响很小.总的来说，本文构建了一种的新的姿态Kalman滤波算法，解算的姿态数据充分发挥了IMU的性能精度优于JPL，但对时变重力场模型解算精度的影响很小，目前限制时变重力场模型的精度主要来源于其他误差.

致谢感谢GFZ网站提供的GRACE-FO Level-1A和Level-1B数据（ftp:∥isdcftp.gfz-potsdam.de/grace-fo/）.感谢审稿专家对本文提出的宝贵修改意见.

附录A 四元数的定义和运算

称q=q0+q1i+q2j+q3k为一个四元数，其中q0,q1,q2,q3∈R，i,j,k是不同的虚数单位，满足i2=j2=k2=-1，ij=-ji=k，jk=-kj=i，ki=-ik=j.四元数也可以表示成（q0,q1,q2,q3），其中q0是标量.

令任意三个四元数为q=（q0,q1,q2,q3），p=（p0,p1,p2,p3）和s=（s0,s1,s2,s3），那么四元数的运算法则定义如下.

加法：

q±p=（q0±p0,q1±p2,q1±p2,q1±p2），

（A1）

点乘：

q·p=q0p0+q1p1+q2p2+q3p3，

（A2）

乘法：

s=qp,

（A3）

共轭：

（A4）

模：

（A5）

若‖q‖=1，则称q是单位四元数.

逆：

（A6）

单位四元数三角函数形式及指数形式：

q=cosθ+nsinθ，

（A7）

其中n∈S2为单位矢量.

猜你喜欢 重力场星间时变 重力场强度在高中物理中的应用科技视界(2020年14期)2020-07-30基于星间链路的导航卫星时间自主恢复策略空间科学学报(2020年3期)2020-07-24基于空间分布的重力场持续适配能力评估方法中国惯性技术学报(2019年6期)2019-03-04星间激光通信终端精跟踪性能地面测试和分析制导与引信(2017年3期)2017-11-02基于时变Copula的股票市场相关性分析智富时代(2017年4期)2017-04-27基于时变Copula的股票市场相关性分析智富时代(2017年4期)2017-04-27非完整性约束的平面多智能体位置时变一致性控制自动化学报(2017年7期)2017-04-18例谈带电粒子在复合场中的运动分类新高考·高一物理(2015年6期)2015-09-28星间测距对小卫星编队相对轨道状态的修正空间控制技术与应用(2015年3期)2015-06-05基于星间测距的编队卫星一致性导航算法空间控制技术与应用(2015年3期)2015-06-05