一种重力卫星质心在轨标定算法

辛宁 邱乐德 张立华 丁延卫 王大雷

（1航天东方红卫星有限公司，北京100094）（2中国空间技术研究院，北京100094）

1 引言

重力场探测卫星的核心有效载荷为测量非保守力的静电加速度计，用于确定星间距离的K波段微波星间双向测距仪（K Band Ranging System，KBR）。当卫星质心发生漂移时，静电加速度计内部悬浮质量块质心与卫星质心无法重合，严重影响了加速度计的测量精度［1］，同时卫星质心的偏差对KBR相位中心的几何标定也有一定的影响［2］。然而，由于卫星在发射前质心校准存在误差，发射中经受振动环境，卫星入轨后姿态调整飞行中燃料损失等原因均会引起卫星质心发生漂移。为了使卫星在轨飞行测量中，达到悬浮质量块质心和卫星质心的最大重合，提高静电加速度计和KBR的测量精度，进而反演出高精度的地球重力场信息，必须进行重力卫星质心的在轨标定研究［3］。

文献［4］提出了一种利用卫星自旋晃动，并通过静电加速度计实现卫星质心在轨标定的方案。但是该在轨自旋方案控制难度较大，工程上难以实现。文献［5］从工程实现的角度出发，提出了一种利用周期性磁力矩使卫星产生姿态机动，利用经典的Kalman滤波算法对卫星质心进行检测的方案。此方案将卫星质心在轨标定问题转换为对卫星角加速度的求解，并进行了尝试性的仿真研究。但该算法需假设卫星加速度模型近似线性，且过程噪声为高斯白噪声，如果实际的角速度响应中存在着非线性、非平稳过程，则会降低质心估计精度。文献［6］提出了一种利用最小模型误差准则（MME）估计卫星角加速度，利用EKF估计卫星质心的算法。由于该算法的测量元件中增加了陀螺仪，因此不需要再进行卫星角速度的估计，节省了一定的计算量。MME算法虽然无需利用假设的系统模型即可获得角速度估计，但是MME估计器是一种离线估计算法，无法对质心进行实时估计。为了解决实时性问题，在经典的卡尔曼滤波和MME估计器的基础上，文献［7］提出了一种应用于无陀螺角速度确定预测滤波算法，这种算法不再局限于系统模型误差只是高斯白噪声的情况，能够对空间飞行姿态进行实时估计，有效解决了存在显著动态模型误差情况下的非线性滤波问题。

在上述文献的基础上，提出了一种使用静电加速度计和陀螺数据估计卫星质心的算法，其将预测滤波与扩展卡尔曼滤波相结合，以卫星角加速度和质心偏移量为状态估计量，将卫星角加速度作为系统模型误差可以实时进行估计。该算法中的滤波器不需要卫星的控制输入以及动力学特性信息，针对卫星质心及角加速度能够实时估计，具有快速的收敛性能和良好的稳定精度。

2 卫星质心在轨标定算法

卫星质心的在轨标定基本原理如图1所示，根据地球磁场模型及卫星的轨道位置，利用磁力矩器给卫星施加一定频率的力矩，使卫星产生相应的转动运动。采用加速度计、陀螺仪作为检测传感器，对陀螺仪观测矢量进行预测滤波处理，估计出角加速度；对静电加速度观测矢量进行扩展卡尔曼滤波处理，分离出具有线性特征的非重力加速度分量，由此来确定卫星质心的位置。

图1 重力卫星质心在轨标定原理Fig.1 Principle of on-orbit calibration of center of mass for gravity satellite

2.1 系统模型

静电加速度计在卫星本体坐标系下输出的观测矢量表达式可表示为［5］

式中Aobs为静电加速度计在卫星本体坐标系下的输出；RAB为安装误差矩阵；s1为加速度计的标度因数；Ab为零偏误差；An为加速度计的噪声误差；aout为静电加速度计的实际输出：

式中ω为卫星角速度观测数据；为卫星角加速度；l为卫星质心偏差；gg为由于质心偏差引起的重力梯度加速度；ang为非保守力加速度。

陀螺仪的观测数据可表示为：

式中k1为标度因数；ωs为加速度计输出；ωb为零偏误差；ωn为陀螺的噪声误差。

由文献［5］可知，由重力梯度和非保守力因素引起的加速度分量在较短的标定时间t内基本保持线性关系，因此质心在轨标定问题的测量方程可表示为

式中α和β为aout、gg以及Ab线性化近似后的斜率估计量和截距估计量；为等效矩阵：

显然，在构建上述测量方程时，需要利用卫星的角加速度信息。刚体卫星角加速度为

式中Jp为转动惯量；M为作用于卫星上的合力矩。由于Jp和M在质心标定过程中存在不确定性，因此利用公式（4）无法实时准确的得到卫星角加速度。为了弥补上述公式的缺点，本文首先将陀螺数据代入预测滤波器中，实时获取到角加速度信息，然后利用扩展卡尔曼滤波对质心位置进行标定。

将式（4）右边的状态信息全部定义为模型误差d（t），由此获得的系统状态方程为

式中G为d（t）的系数矩阵。

观测方程为

式中为陀螺仪输出；为静电加速度计输出。观测值

2.2 角加速度预测滤波估计

预测滤波器的物理概念：设已经获得了tk时刻的状态估值，采样间隔时间为Δt，在接到tk＋Δt时刻的测量信息后，利用模型最小误差逼近可预测 ［tk，tk＋1］间隔内的模型误差d（t），然后将d（t）的预测值代入系统状态方程和观测输出方程中，将状态估值传播到tk＋1时刻，从而得到较高的估计精度。

将观测方程（7）泰勒级数展开得到［7］

将式（9）代入式（10）后求导，得到最小模型误差：

其中

式中ci为系统观测模型函数c的第i行；Λ（Δt）为对角矩阵，其对角元素为λii＝Δtpi／pi！，pi是对公式（7）微分时中首次出现d（t）中元素的最低阶数；是ci关于f的k阶李导数；si为S的第i行元素；j为d（t）的行索引；gj是矩阵G的第j列元素；R为量测噪声；W为半正定模型误差加权矩阵。

利用上述迭代公式便可实时精确地获取卫星角加速度˙ω（t）的估计值。

2.3 质心的扩展卡尔曼滤波估计

在卫星角速度及角加速度估计值已知的条件下，根据系统状态方程式（6）及观测方程式（8）对X1＝［lαβ］进行估计。

待估状态变量一步预测均值估计为

式中Φ（tk）＝Ⅰ9×9为状态转移矩阵。

一步估计均方误差阵为P（tk＋1／k）＝Φ（tk）P（tk）ΦT（tk）

反馈增益矩阵为K（tk＋1）＝P（tk＋1／k）H（tk＋1）T（H（tk＋1）P（tk＋1／k）H（tk＋1）＋Q（tk＋1））-1，其 中，为量测转移矩阵；Q（tk＋1）为量测噪声矩阵。

估计均方差阵为P（tk＋1）＝（Ⅰ-K（tk＋1）H（tk＋1））P（tk＋1／k）

状态估计值：

上述质心标定过程的实质是先利用预测滤波预测下一时刻的模型误差校正量d（t），然后再对状态估值进行传播估计，从而有效地提高了角加速度的估计精度。然后将角速度及高精度的角加速度估计值代入中更新量测转移矩阵H，根据公式（12）实现了卫星质心的高精度估计。

3 数值仿真

卫星初始时刻：2015年1月1日正午（UTC）。

卫星动力学参数［8］：10阶EMG96重力场模型；DTM大气模型，阻力系数Cd＝2.2；太阳光压反照系数η＝0.5；地球磁场模型，13阶IGRF2005模型。

卫星仿真参数［9］：质量m＝500kg；初始质心位置三轴偏差磁力矩器力偶mm＝ ［0 1 -sgn（ByBz）］T·30sin（2π×0.1t）（A·m2），By、Bz为地球磁场在卫星本体坐标系下的分量；静电加速度计标度因数s1＝ ［1＋0.01% 1＋0.01% 1＋0.01%］T；零偏误差Ab＝ ［2×10-62×10-6-2×10-6］Tm／s2，测量误差An＝ ［1×10-101×10-101×10-10］Tm／s2；陀螺仪 标 度 因 数k1＝ 1＋0.1% ，三轴零偏误差ωb＝［0.3 0.3 0.3］／［（°）／h］，测量误差ωn＝［0.1 0.1 0.1］／［（°）／h］。

加速度计与陀螺仪的输出周期为0.1s，预测估计滤波算法的周期为0.1s，为了检验所设计的磁力矩器力偶信号的持续时间与质心标定精度的关系，预测估计滤波算法的运行时间分别设定为1min、3min、5min、10min。为了分析得到数据采集的轨道相位与质心标定精度的关系，分别选在赤道附近及南极附近进行卫星质心标定。

在赤道附近进行姿态机动，图2为地球磁场强度仿真结果，可以看出三轴磁场强度在短时间内可认为是近似线性变化。图3为施加在卫星上的周期性磁力矩仿真结果，三轴磁力矩的最大幅值分别为0.5mNm、0.75mNm和0.75mNm，形式为正弦，频率为0.1Hz。线性加速度的仿真时序如图4所示，其中X方向的线性加速度近似效果相对较差，主要原因为X方向上大气阻力的线性度比Y方向和Z方向要差很多，因此线性度较差。图5为陀螺仪输出的角速度数据，由图5中的角速度数据计算的角加速度估计值如图6所示，预测滤波算法具有很好的估计精度和可靠性，三轴估计精度分别为2.939×10-6rad／s2，0.975×10-6rad／s2，0.858×10-6rad／s2；该算法仅利用了陀螺提供的角速度便能实时地估计出高精度的卫星角加速度，南极上空的现象与赤道上空一致，这里不再赘述。通过图7中的静电加速度计的输出仿真值，利用扩展卡尔曼滤波得到的质心标定结果如图8所示，由于X方向的线性加速度线性度较其他两个方向相对较差，因此质心X方向的估计时间较长，Y方向和Z方向的估计时间较短。卫星质心在轨标定的轨道相位、精度和标定时间的关系如表1所示。

图2 赤道上空三轴磁场强度分布Fig.2 Magnetic flux density over equator

图3 赤道上空三轴磁力矩Fig.3 Magnetic torque simulation over equator

图4 赤道上空线性加速度Fig.4 Linear acceleration over equator

图5 赤道上空陀螺仪观测值Fig.5 Gyro observation over equator

图6 赤道上空角加速度估计值Fig.6 Estimation of angle acceleration over equator

图7 赤道上空静电加速度计观测值Fig.7 Accelerometer observation over equator

图8 赤道上空质心卡尔曼滤波估计值Fig.8 Kalman estimation of center of mass over equator

表1 卫星质心在轨标定误差精度Tab.1 Center of mass calibration error of satellite mm

由表1可以看出，X轴的最佳标定精度为0.062 8mm，Y轴的最佳标定精度为0.032 4 mm，Z轴的最佳精度为0.041 4mm。本文提出的质心标定算法满足标定精度需求。在相同条件下，在赤道附近进行姿态机动，质心标定结果要略优于在南极附近机动的结果。其主要原因是赤道附近激励起的三轴角加速度较为均匀，这样利于辨识精度的提高；而南极激励的三轴角加速度幅值差异较大。通过不同的激励时长看，3min或5min的总机动时间是可以接受的，再长的激励时间，则质心测量精度反而降低。其主要原因是，辨识中假设非重力加速度项在辨识期间保持线性关系，而随着辨识时间的增加，这种假设的误差就越大，反而不利于质心辨识精度的提高。

4 结束语

综上所述，目前关于重力卫星质心的标定方法均为离线估计算法，无法实现质心的高精度实时估计，针对这一问题，本文通过设计估计角加速度的预测滤波器与质心扩展卡尔曼滤波器实现了卫星质心的实时估计，仿真结果表明该算法具有较高的精度，鲁棒性能强，可为质心在轨修正及KBR天线相位中心在轨标定提供设计依据。同时考虑该算法对陀螺仪的精度要求较高，下一步将研究在陀螺仪精度不高的情况下，如何确保算法的精度。

［1］周泽兵，白彦峥，祝竺，等.卫星重力测量中加速度计在轨参数校准方法研究 ［J］.中国空间科学技术，2009，29（6）：74-80.ZHOU ZEBING，BAI YANZHENG，ZHU ZHU，et al.In-orbit calibration methods of accelerometer parameters on satellite-borne gravimetry ［J］.Chinese Space Science and Technology，2009，29（6）：74-80.

［2］佘世刚，王锴，周毅，等.高精度星间微波测距技术 ［J］.宇航学报，2006，27（3）：403-406.SHE SHIGANG，WANG KAI，ZHOU YI，et al.The technology of high accuracy inter-satellite microwave ranging ［J］.Journal of Astronautics，2006，27（3）：403-406.

［3］陈光锋，唐富荣，薛大同.重力卫星在轨质心修正原理 ［J］.宇航学报，2005，26（5）：567-570.CHEN GUANGFENG，TANG FURONG，XUE DATONG.The trim principle of center of mass of gravity satellite during orbit flight［J］.Journal of Astronautics，2005，26（5）：567-570.

［4］TANGYIN S，WILLIAMS T.Mass property estimation using coasting maneuvers［J］.Journal of Guidance，Control and Dynamics，1997，20（4）：625-632.

［5］WANG F，TAPLEY B D.Precise determination of accelerometer proof mass ［C］.AAS／AIAA Space Flight Mechanics Meeting，Texas，2002：739-752.

［6］王本利，廖鹤，韩毅.基于MME／EKF算法的卫星质心在轨标定 ［J］.宇航学报，2010，31（9）：2150-2156.WANG BENLI，LIAO HE，HAN YI.On-orbit calibration of satellite center of mass based on MME／EKF algorithm ［J］.Journal of Astronautics，2010，31（9）：2150-2156.

［7］CRASSIDIS J L.Predictive filtering for attitude estimation without rate sensors ［J］.Journal of Guidance，Control and Dynamics，1997，20（3）：522-527.

［8］KIM J R.Simulation study of a low-low satellite-to-satellite tracking mission ［D］.Austin：University of Texas at Austin，2000.

［9］WANG FURUN.Study on center of mass calibration and K-band ranging system calibration of the GRACE mission［D］.Austin：University of Texas at Austin，2003.