为了得到准确的天宫二号动力学质心位置，需要在精密定轨的同时完成对动力学质心的估计。另外，质心估计还有其他重要作用，如结果可用于评估大型航天器内部任意位置点的微重力水平[5]。经测算，距离动力学质心1 m处的微重力水平为1.35×10-6 m/s2，距离动力学质心10 m处微重力水平将下降一个数量级，即1.35×10-5 m/s2，某些实验对微重力水平和质心位置精度要求较高。当天宫二号和天舟一号组合体运行时，天舟一号上的加速度计无法安装在质心处，由质心偏移而引起的离心力加速度和线性加速度则干扰了加速度计对非保守力的测量。因此当完成组合体动力学质心估计后，可以有效剥离加速度计数据中由于质心偏移安装所引入的离心力加速度和线性加速度干扰项，精确评估加速度计对非保守力加速度的测量水平和精度。航天器质心的轨道状态还可用于大气密度的反演[6]。获得航天器精确的动力学质心后，可以利用GNSS测量数据和非质心安装加速度计数据对航天器的姿态角进行估计[7]。
目前，国内外已经针对航天器的质心标定问题开展了研究[8]，多数文献采用对敏感器数据进行连续处理，数据量大、计算复杂，且实验对象多为重力卫星这类小型航天器。文献[9]提出了一种利用卫星自旋晃动，并通过静电加速度计实现卫星质心在轨标定的方案，通过该方案计算得到的质心精度达到毫米级，但是该自旋方案控制难度较大，工程上难以实现。文献[10-11]针对重力卫星质心在轨标定问题，提出了一种利用周期性磁力矩使卫星产生姿态机动，利用经典的Kalman滤波算法对卫星质心进行检测的方案，该方案将卫星质心在轨标定问题转换为对卫星角加速度的求解，并进行了尝试性的仿真研究，使质心位置的三轴最佳标定精度达到毫米级，但该算法需要假设加速度模型近似线型，且过程噪声为高斯白噪声，如果实际的角速度响应中存在非线型、非平稳过程，则会降低质心估计精度。文献[12-13]利用陀螺仪和加速度计的测量数据，提出了基于指数加权递归的最小二乘算法来估计航天器的质心位移，计算得到质心位置的估计精度达到约6 mm水平。以上方法大都不适合天宫二号，未搭载加速度计等敏感器。
本文分析了处于偏航状态下的天宫二号质心位移估计的可行性，提出了基于GNSS测量数据和简化动力学最小二乘批处理(Batch Least-Squares, Batch LSQ)的质心位移估计方法。大多数航天器都会配置GNSS接收机，将GNSS测量数据用于航天器质心估计具有普遍意义。同时，针对天宫二号结构和姿态运动特点，构建了较为精细的结构模型Cylinder-wing[14]，实时计算大气阻力和辐射光压面质比，提高非保守力建模精度。由于天宫二号在轨期间，其真实质心位置无法确定，缺乏直接验证质心位移估计精度的手段，只能通过定轨结果的内外符合来间接验证。因此，本文从后验残差、重叠轨道检验和地面激光数据检核等方面对质心位移校正前后的精密轨道进行了对比，从而评估该质心估计方法的有效性。
1 天宫二号结构及姿态动力学模型 1.1 天宫二号结构 天宫二号由实验舱和资源舱2个舱段构成，其总质量约为8 600 kg[3]，全长10.4 m，最大直径3.35 m，太阳能翼宽约18.4 m，结构如图 1所示。资源舱后端面与运载火箭连接。从资源舱后端沿空间实验室纵轴向实验舱看，按顺时针方向分为4个象限，由象限线Ⅰ、Ⅱ、Ⅲ、Ⅳ分割。各象限线间隔90°，正常在轨飞行时Ⅰ象限对地，如图 1所示。
 |
| 图 1 天宫二号结构和象限线示意图 Fig. 1 Illustration of Tiangong-2 structure and quadrant line definition |
| 图选项 |
天宫二号的整体坐标系和本体坐标系[15]如图 1所示。整体坐标系原点Od取在资源舱与运载火箭连接端面的几何中心；OdXd轴沿飞行器的纵轴，指向对接机构；OdYd轴与OdXd轴垂直，指向Ⅲ象限线；OdZd轴与OdXd和OdYd轴垂直，指向Ⅳ象限线。本体坐标系原点Ob为空间实验室质心，ObYb轴指向空间实验室的Ⅳ象限线，ObZb轴指向空间实验室的Ⅰ象限线，ObXb轴、ObYb轴和ObZb轴满足右手螺旋法则。
轨道坐标系OoXoYoZo的坐标原点位于空间实验室质心处。+Zo轴指向地心，+Xo轴在轨道平面内与Zo轴垂直并指向与空间实验室运动速度呈锐角的方向，Yo轴根据Xo轴和Zo轴的指向，通过右手螺旋法则确定，垂直于轨道面，如图 2所示。轨道坐标系的各坐标轴单位矢量的数学表达式为
 | (1) |
 |
| 图 2 轨道坐标系示意图 Fig. 2 Illustration of orbital coordinate system |
| 图选项 |
式中：r为从地心指向卫星质心Oo的位置矢量；v为卫星的运动速度矢量。
天宫二号在运行期间，开展了精密定轨任务，测量载荷包括双频GNSS接收机和激光反射器(Laser Reroreflector Array, LRA)，其在整体坐标系下的坐标如表 1所示。天线安装角为34°，天线法线与整体坐标系+Zd轴的夹角为14°，如图 3所示。激光发射器安装在天宫二号象限线Ⅰ方向上，指向地心。
表 1 测量载荷与发射前参考质心的位置 Table 1 Positions of measured loads and reference COM before launch
| 位置 | Xd/cm | Yd/cm | Zd/cm |
| 天线相位中心 | 428.66 | 158.5 | -101.31 |
| 激光角反射镜 | 414.88 | -159.57 | -0.09 |
| 发射前参考质心 | 416.03 | 0.63 | 0.53 |
表选项
 |
| 图 3 天线安装位置及天线法线示意图 Fig. 3 Illustration of position and normal direction of antenna |
| 图选项 |
1.2 天宫二号姿态运动 天宫二号航天器有2种姿态模式，分别为“三轴对地稳定”和“连续偏航”。当太阳高度角β>25°时，天宫二号将执行连续偏航模式。在这种模式下，卫星本体坐标系与轨道坐标系在本体坐标系的XbYb平面内存在夹角，称作姿态偏航角ψ(定义参见附录A)，如图 4所示。同时姿态偏航角ψ连续发生变化，保持太阳能电池板垂直于太阳方向，即太阳能帆板法向n与太阳方向一致，可以获得最大的太阳能。姿态偏航角ψ和姿态偏航角速度
 | (2) |
 | (3) |
 |
| 图 4 天宫二号偏航姿态运动示意图 Fig. 4 Illustration of Tiangong-2 yaw-steeringattitude mode |
| 图选项 |
如果λ>0，则
 | (4) |
否则
 | (5) |
式中：ω为轨道角速度绝对值；Sx、Sy和Sz为太阳矢量表示在轨道坐标系的单位向量。
当太阳高度角β < 25°时，天宫二号将转向三轴对地稳定模式。在三轴对地稳定模式下，卫星本体坐标系与轨道坐标系重合，太阳能电池板将尽量保持最大的横截面对太阳，所以太阳能帆板法向n与太阳方向在轨道平面上的投影一致。太阳能帆板法向n随着飞行方向(整体坐标系OdXd轴)的变化而不断变化[15]。
以一个完整的轨道周期为例，姿态偏航角ψ随轨道升交角距的变化如图 5所示。可知，随着升交角距的增大，姿态偏航角在40°~140°之间连续变化。
 |
| 图 5 姿态偏航角变化示意图 Fig. 5 Illustration of yaw angle change |
| 图选项 |
1.3 动力学精确建模 由天宫二号的姿态运动可知，天宫二号处于不同姿态模式时，大气阻力截面积不一样，需要进行精细的动力学建模。在动力学建模时，将天宫二号的卫星结构模型简化为Cylinder-wing模型，该结构模型由1个圆柱形的舱体和2个长方形的太阳翼构成。圆柱高H=10.4 m，底面积直径R=3.35 m，l为太阳能帆板的宽度，h为太阳能帆板的长度，约18.4 m，其侧视图和俯视图如图 6所示。
 |
| 图 6 Cylinder-wing模型示意图 Fig. 6 Illustration of model of Cylinder-wing |
| 图选项 |
根据天宫二号的姿态运动推导出2种模式下的大气阻力截面积的计算公式。
1) 计算大气速度的单位方向向量
在地心惯性坐标系[16]下，卫星的运动速度为V，卫星质心的位置向量为Rs, 则大气速度的单位方向向量为
 | (6) |
2) 计算圆柱体的大气阻力截面积
三轴对地稳定姿态下，圆柱体的大气阻力截面积为
 | (7) |
偏航姿态下，假设姿态偏航角为ψ，圆柱体的大气阻力截面积为
 | (8) |
3) 计算太阳能帆板的大气阻力截面积
在地心惯性坐标系下，太阳的坐标向量为Rsun，则卫星指向太阳的向量Rlos为
 | (9) |
将Rlos投影到地心惯性坐标系的XY面上，同时单位化得到Rvlos：
 | (10) |
三轴对地稳定姿态下，太阳能帆板的大气阻力截面积为
 | (11) |
偏航姿态下，太阳能帆板的大气阻力截面积为
 | (12) |
可知，总大气阻力截面积为
 | (13) |
据此即可计算出飞行器的大气阻力面质比，图 7为不同姿态模式下，大气阻力面质比在一个轨道周期内随升交角距的变化情况。大气阻力面质比随升交角距连续变化，变化范围约为0.001~0.006 m2/kg。
 |
| 图 7 偏航模式和三轴对地稳定模式下，大气阻力面质比在不同升交角距下的变化规律 Fig. 7 Law of the change of atmospheric drag mass of atmospheric drag in attitude mode of yaw steering and in attitude mode of three-axis Earth-pointing stabilization under different argument of latitude |
| 图选项 |
2 质心位移估计方法 2.1 质心位移建模与估计 简化动力学最小二乘批处理方法是对动力学模型参数和观测模型参数进行迭代平差，使得给定时间段上的观测量和观测模型值的差的平方和最小。根据待估参数初值条件，卫星在给定时间段上进行数值积分，得到观测模型值，与观测量做差，得到残差值。建立观测模型关于待估参数的偏导数矩阵，组成超定线性法方程组，求解法方程得到初值条件的校正量，多次迭代直至待估参数的解收敛[17]。该方法每次迭代完成给定时间段上所有观测数据的处理，观测量必须事先可用，满足事后轨道确定的条件。
最小二乘批处理方法的具体推导过程可参考文献[15]，在此基础上增加质心偏移估计参数。将质心偏移加入动力学待估参数向量Y中，即
 | (14) |
式中：roffsets为质心位移量; y0为卫星在t0时刻的位置与速度；CD和CR分别为大气阻力系数和太阳辐射光压系数；a0, a1, …, ana－1表示t0~tna－1时刻的逐段经验加速度，在定轨过程中同其他参数一同估计，用于吸收动力学模型误差，其量级大小可表示非保守力建模误差[18]。
GNSS伪距消电离测量模型[17]为
 | (15) |
式中：r为飞行器；s为某一颗GNSS星座卫星；t为信号接收时间；c为光速; ρrs(t)为接收机天线相位中心和GNSS卫星天线相位中心的几何距离；δtr(t)为接收机在t时刻的钟差；δts(t-τrs(t))为GNSS卫星在信号发射时刻t-τrs(t)的钟差，τrs(t)为信号传输时间；εRIF(t)为测量噪声，假设为均值为0的随机变量，精度在分米级; MRIF(t)包括天线相位中心校正、重力场延迟等。
载波相位消电离测量模型为
 | (16) |
式中：λIF为波长；NIF为载波相位的整周数；εLIF(t)为测量噪声；MLIF(t)为其他系统偏差。
对ρrs(t)进行线性化，得到
 | (17) |
式中：ρr0s(t)为GNSS卫星质心与天宫二号参考质心的距离；yr0(t)为天宫二号参考质心在地心地固下的位置；ys(t)为GNSS卫星质心在地心地固系下的位置；rpco, m为天宫二号动力学质心到接收机天线相位中心的向量；rpco为天宫二号参考质心到接收机天线相位中心的向量；roffsets=rpco, m－rpco为质心位移偏差，即作为待估参数，在批处理过程中假设是常量；ers(t)为天宫二号指向GNSS卫星的单位向量，如图 8所示。
 |
| Fig. 8 Positional relationship among satellite reference COM, dynamic COM and antenna phase center |
| 图选项 |
采用微分修正的方式对质心位移量进行估计，需计算测量值对质心位移偏差的偏导数。
 | (18) |
2.2 质心估计可行性分析 通过定性分析和定量分析2个方面来说明连续偏航模式下进行质心估计的可行性。
2.2.1 定性分析 通过以下假设来定性分析是否可以基于GNSS测量实现对卫星质心位移的估计。假设卫星内部有2个位置点：一个是真正的动力学质心，一个是参考质心。因参考质心位置已知，便可以和GNSS天线相位中心建立对应关系，或者说，GNSS“间接”测量的是航天器参考质心的位置。而动力学质心的运动反映的是航天器真实的运动状态，不受GNSS测量的影响。若GNSS测量值(测量参考质心的运动)能与动力学质心的运动完全匹配，即参考质心和动力学质心的轨迹一致，则GNSS测量是无法区分这2个位置点的，也就是不具备质心估计可行性。若GNSS测量值与动力学质心运动不匹配，则可以通过GNSS测量进行质心位移估计。
1) 假设这2个质点仅在轨道坐标系中的Zo方向(即轨道径向)存在位置差，则两质点由于轨道高度不相同造成轨道周期存在差异，随着时间增加，GNSS测量的参考质心与动力学质心的位置差越来越大。因此，轨道径向上的质心位移是可以估计的。
2) 假设这2个质点仅在轨道坐标系中的Yo方向(即轨道法向)上存在位置差，则两质点运动的轨道面存在一定夹角，位置差发生周期性变化。因此，轨道法向上的质心位移是可以估计的。
3) 假设这2个质点仅在轨道坐标系中的Xo方向(即轨道切向)上存在位置差，在两质点与地心的方向矢量形成的夹角足够小的情况下(由于航天器尺寸限制，该条件一般都能满足)，2个质点的轨迹是相同的，其相对位置关系固定不变，因此，轨道切向上的质心位移是无法估计的。而且由于2个质点的位置差很小，非保守力也不会有明显差异。
以上分析表明，轨道径向和法向的质心位移可观测性比较强，而轨道切向上的质心位移可观测性较差。燃料消耗引起的天宫二号质心位移本体坐标系Xb轴方向变化显著。根据上述定性分析，在三轴对地稳定姿态下，本体坐标系Xb轴与轨道切向重合，难以通过GNSS测量正确本体坐标系Xb轴的质心位移。而在连续偏航模式下，本体坐标系Xb轴在轨道法向上有较大分量。因此，本文利用连续偏航模式下的GNSS观测数据，完成对天宫二号卫星的质心确定。
2.2.2 定量分析 利用GNSS实测数据计算本体坐标系下不同方向的质心偏移误差对精密轨道的影响。三轴对地稳定情况下的定量仿真分析已有文献给出[19]。因此，本文选取天宫二号卫星处于偏航姿态下的时期(姿态偏航角变化范围为40°~140°)，即利用天宫二号卫星2018年1月8日—1月10日的观测量。在不考虑质心位移的前提下，采用简化动力学精密定轨方法计算得到轨道产品Orbit-noOFFSET。分别在天宫二号卫星本体坐标系的Xb、Yb、Zb 3个方向上添加10 cm的质心位置位移，即roffsetsX=(10, 0, 0)cm、roffsetsY=(0, 10, 0)cm和roffsetsZ=(10, 0, 0)cm。仍采用相同的定轨方法对上述日期获得的观测数据进行解算，获得轨道产品Orbit-REF，并将其与Orbit-noOFFSET进行比对，分析不同方向的质心位移变化对定轨结果的影响，结果如表 2~表 4所示。表中：Rx、Ry、Rz分别表示Orbit-REF在本体坐标系3个方向的位置残差(与Orbit-noOFFSET比较)。
表 2 本体坐标系Xb方向10 cm的质心位移变化对定轨结果的影响分析 Table 2 Analysis of influence of 10 cm COM offset in Xb-axis direction on orbit determination results
| 日期 | Rx/cm | Ry/cm | Rz/cm |
| 2018-01-08 | 6.39 | 6.81 | 2.58 |
| 2018-01-09 | 5.69 | 7.28 | 1.67 |
| 2018-01-10 | 5.42 | 7.63 | 1.85 |
| 平均值 | 5.83 | 7.24 | 2.03 |
表选项
表 3 本体坐标系Yb方向10 cm的质心位移变化对定轨结果的影响分析 Table 3 Analysis of influence of 10 cm COM offset in Yb-axis direction on orbit determination results
| 日期 | Rx/cm | Ry/cm | Rz/cm |
| 2018-01-08 | 7.8 | 6.3 | 1.00 |
| 2018-01-09 | 7.49 | 6.52 | 0.85 |
| 2018-01-10 | 7.83 | 6.58 | 1.04 |
| 平均值 | 7.71 | 6.47 | 0.96 |
表选项
表 4 本体坐标系Zb方向10 cm的质心偏移变化对定轨结果的影响分析 Table 4 Analysis of influence of 10 cm COM offset in Zb-axis direction on orbit determination results
| 日期 | Rx/cm | Ry/cm | Rz/cm |
| 2018-01-08 | 3.18 | 0.61 | 8.14 |
| 2018-01-09 | 2.75 | 0.53 | 8.86 |
| 2018-01-10 | 3.08 | 0.44 | 8.65 |
| 平均值 | 3.00 | 0.53 | 8.55 |
表选项
本体坐标系Xb方向加入10 cm的质心位移roffsetsX后，得到的定轨结果与Orbit-noOFFSET进行对比，发现本体坐标系Xb方向定轨残差均值为5.83 cm，Yb方向为7.24 cm, Zb方向为2.03 cm，与质心位移roffsetsX对比有明显区别，说明质心位移影响到了GNSS测量模型值，导致测量模型值与实际测量值不再匹配，从而影响了定轨结果。因此，偏航状态下，GNSS测量数据可以用于本体坐标系Xb方向上的质心位移估计。
由于天宫二号在轨机动的燃料消耗，本体坐标系Xb方向上的质心位移最大，远高于另外2个正交方向的质心位移。在三轴对地定向姿态模式中，本体坐标系与轨道坐标系重合，本体坐标系Xb方向质心位移与轨道切向重合。谷德峰等[19]已证明轨道切向的质心位移不可估计。但在连续偏航模式下，本体坐标系Xb方向质心位移在轨道切向和轨道法向上都有分量。现假定Xb方向存在10 cm的质心位移，由于姿态的偏航运动，该质心位移在轨道切向和法向上的分量如图 9所示，轨道法向的质心位移分量在6~10 cm之间变化，大部分情况下，质心位移在法向上的分量比在切向上的分量要大。此时本体坐标系Xb方向质心位移是可以估计的[20]。
 |
| 图 9 偏航模式下，本体坐标系Xb轴10 cm质心位移在轨道切向和法向上的分量变化 Fig. 9 Component changes of Xb-axis 10 cm-offsets of COM in along-track and cross-track directions in attitude mode of yaw steering |
| 图选项 |
本体坐标系其他2个方向情况大致类似，质心位移对定轨结果都有影响。通过质心估计定量分析结果看，姿态偏航模式下本体坐标系3个方向均具备质心位移估计的可行性。
2.3 质心估计校验方法 本文利用GNSS实测数据在偏航模式下估计天宫二号本体坐标系3个方向的质心位移，并将定轨结果与未进行质心估计的定轨结果做比较，从4个方面评估质心位移估计精度，前3个为内符合精度校验，最后1个为外符合精度校验。
1) 载波相位后验残差精度校验。载波相位后验残差精度表示定轨结果与观测资料之间的拟合程度。计算若干历元定轨结果对应的模型值与观测量之间的载波相位后验残差，得到均值和均方根(Root Mean Square，RMS)。若均值接近于0，RMS值符合观测量类型的噪声水平，说明内符合精度较高，定轨模型与GNSS观测量是一致的，轨道确定的迭代过程是收敛的。
2) 弧段重叠精度校验。相同的定轨软件通过2次独立解算得到的2段轨道中如果包含有重叠弧段，那么尽管重叠段的观测数据相同，但是这2段重叠定轨结果是相对独立的，可以认为不相关，轨道重叠部分的符合程度反映出定轨精度(采用重叠部分之差的RMS值)。
3) 经验加速度估计水平校验。同时段的观测资料独立解算得到的经验加速度估计幅值的大小反映出定轨过程中采用的动力学模型的精度。经验加速度幅值越大，动力学建模精度越差，间接反映质心位移估计精度不高，反之则说明质心位移估计精度较高。
4) 激光测距(Satellites Laser Ranging，SLR)外部校验。利用激光测距观测值对轨道数据进行外部检验是一种独立的轨道检验方法，其原理是：比较地面激光测距站直接测量的站星距(测站和卫星轨道间的几何距离)与基于GNSS测量值求解的轨道与测站的站星距的差，并进行统计分析。
3 质心估计结果和分析 发射前，天宫二号在地面进行了参考质心的标定。天宫二号燃料消耗几乎只影响其本体坐标系Xb轴方向的质心变化。表 5列出了天宫二号发射时刻和3.1节里的实验方案中所涉及时间段的质心位置(称作参考质心，表示在整体坐标系中)和天线相位中心位置(表示在本体坐标系中)，这2个位置是根据燃料消耗和质心变化推算的。因为天宫二号内部结构较为复杂，数学推算的质心位置精度比较差，所以有必要利用在轨测量数据进行精确估计。
表 5 根据燃料消耗推算的天宫二号质心位置和天线相位中心位置 Table 5 Position of COM and antenna phase center of Tiangong-2 computed by fuel consumption
| 阶段 | 发生时间(UTC) | 质量/kg | 质心位置(整体坐标系)/cm | 天线相位中心位置(本体坐标系)/cm | |
| 天宫二号整器 | 剩余燃料 | ||||
| 轨控 | 2016-09-16 08:58:55 | 8 406.4 | 692.4 | (416.03, 0.63, 0.53) | (12.63, -101.84, -157.87) |
| 阶段A~阶段C | 2018-01-01 00:00:00 | 8 488.97 | 889.88 | (393.29, 1.55, 0.53) | (35.37, -101.84, -156.95) |
表选项
3.1 实验方案 本文基于天宫二号实测GNSS双频伪距和载波相位测量数据，在定轨过程采用了2种方案，并对定轨结果进行对比。方案1为不考虑质心位移估计的精密轨道确定；方案2为考虑质心位移估计的精密轨道确定[21]。通过这样的对比实验来验证质心估计对天宫二号轨道精度的影响。
考虑到GNSS观测数据充足，激光数据受观测时段和测站的影响而有限，因此，在GNSS观测数据连续且充足的时段进行质心估计，在有激光数据的时段进行激光校验[20]。具体时段设置如表 6所示。
表 6 所选实验时段的具体说明 Table 6 Specific instructions for selected experimental period
| 时段名称 | 时段 | 时段内天宫二号姿态模式 | 时段内是否有激光观测数据 | 具体实验方法 |
| 阶段A | 2018-01-01—2018-01-10 | 连续偏航 | 无 | 对本体坐标系Xb、Yb、Zb 3个方向同时进行质心位移估计，给出质心位移，并对比质心估计前后的定轨结果 |
| 阶段B | 2018-08-20—2018-08-24 | 三轴对地稳定 | 有 | 利用阶段A估计的质心位移进行质心校正后，完成精密定轨，并对轨道进行激光外符合校验 |
| 阶段C | 2018-08-27—2018-08-30 | 连续偏航 | 有 | 利用阶段A估计的质心位移进行质心校正后，完成精密定轨，并对轨道进行激光外符合校验 |
表选项
当天宫二号卫星处于阶段A、阶段B、阶段C这3个时段时，天宫二号卫星均没有燃料消耗，可视作质心没有发生变化。
3.2 结果与分析 选取来自实际任务的天宫二号卫星多天的在轨观测数据，如表 6所示。GNSS卫星星历以及钟差来自于IGS(International GNSS Service)，为已知量且具有较高精度。GNSS星历位置利用8阶拉格朗日插值以获得GNSS卫星在任意时刻的位置，误差在5 cm量级以下[1]，对GNSS卫星钟差利用线性插值可获得1 cm以下的精度。
定轨相关参数设置为：批处理一天(24 h)的观测数据，观测历元间隔30 s，待估钟差参数为nT=2 880，载波相位模糊度参数nB≈450~500。CD和CR在整个区间上是待估常量，设置经验加速度相关时间te=600 s，动力学参数nY≈443。待估参数总数约为3 800个。计算法方程的权重矩阵，设置消电离伪距噪声标准差为0.4 m，消电离载波相位噪声标准差为0.01 m。动力学轨道积分采用变步长RKF78阶方法，可以满足高精度的数值积分要求。
利用2018年1月1日—2018年1月10日(阶段A)的观测数据进行简化动力学最小二乘批处理得到的质心位移估计量如表 7所示。在本体坐标系下，Xb方向质心位移估计量的平均值为-15.84 cm，标准差为0.46 cm，Zb方向的质心位移估计量均值为1.04 cm，标准差为0.16 cm，这2个方向的质心位移估计精度较高。
表 7 偏航姿态模式下，本体坐标系3个方向的质心位移估计值及不确定度 Table 7 Estimates and uncertainties of COM offset in three directions in satellite-body-fixed coordinate system in attitude mode of yaw steering
| 日期 | Xb/cm | Yb/cm | Zb/cm |
| 2018-01-01 | -16.54±0.35 | -7.49±0.74 | 1.20±0.13 |
| 2018-01-02 | -16.09±0.36 | -5.68±0.86 | 1.21±0.13 |
| 2018-01-03 | -16.26±0.41 | -4.69±1.07 | 1.04±0.11 |
| 2018-01-04 | -16.02±0.43 | -6.76±1.38 | 1.11±0.14 |
| 2018-01-05 | -15.74±0.41 | -6.53±1.34 | 1.29±0.16 |
| 2018-01-06 | -16.13±0.44 | -9.57±2.19 | 1.02±0.13 |
| 2018-01-07 | -15.30±0.48 | -7.19±3.03 | 0.94±0.14 |
| 2018-01-08 | -15.54±0.47 | -3.84±3.36 | 0.80±0.14 |
| 2018-01-09 | -15.79±0.53 | -11.25±5.09 | 0.87±0.16 |
| 2018-01-10 | -15.03±0.43 | -10.87±3.27 | 0.88±0.19 |
| 平均值 | -15.84 | -7.39 | 1.04 |
| 标准差 | 0.46 | 2.49 | 0.16 |
表选项
尽管本文应用的经验加速度参数吸收了部分动力学建模误差，其与质心位移也存在一定耦合，对质心位移估计产生一定影响，但由于定轨过程中对经验加速度参数进行了约束，使得估计后的经验加速度幅值的大小被约束在合理的范围内。计算估计后的经验加速度的时间序列如图 10所示。经计算，本体坐标系Xb方向的经验加速度对该方向的质心位移的影响为2~3 cm，其大小约为本体坐标系Xb方向质心估计结果(15 cm)的1/5(见表 7)，不会产生根本性影响。
 |
| 图 10 偏航模式下，质心位移估计后本体坐标系下3个方向的经验加速度 Fig. 10 Empirical acceleration after COM offset estimation in three directions in satellite-body-fixed coordinate system in attitude mode of yaw steering |
| 图选项 |
Yb方向的质心位移估计量均值为-7.39 cm，标准差为2.49 cm。因为天宫二号卫星偏航模式下的姿态偏航角主要在40°~140°之间变化，所以多数情况本体坐标系Yb方向质心位移在轨道切向上的分量比在轨道法向上的分量要大，该方向的质心估计结果相对较差，表 7中Yb方向的质心位移估计标准差比其他2个方向差一些。
3.2.1 内符合分析 为了检验定轨结果与观测量的一致性(内符合)，给出了阶段A时段内的载波相位后验残差的标准差结果，如图 11所示。残差平均值为0，定轨结果没有系统性偏差，未进行质心位移估计时相位后验残差标准差的均值为0.90 cm，进行质心位移估计后相位后验残差标准差的均值为0.86 cm，在进行质心位移估计后，相位后验残差标准差均值明显减小，说明内符合结果更好，轨道精度更高。
 |
| 图 11 2018年1月1日—2018年1月10日的轨道后验残差标准差变化 Fig. 11 Changes of standard deviation of post-fit residuals from January 1st, 2018 to January 10th, 2018 |
| 图选项 |
3.2.2 经验加速度估计水平分析 仍以2018年1月1日—2018年1月10日(阶段A)的观测数据为例，分析经验加速度的补偿水平情况。图 12比较了未进行和进行质心位移估计定轨过程的经验加速度补偿水平。可以看出，与未进行质心位移估计相比，进行质心位移估计后，轨道径向、切向和法向经验加速度分别降低62%、50%和65%，表明质心得到校正后，轨道动力学模型精度提高，相应的经验加速度的补偿水平才会降低。
 |
| 图 12 经验加速度的补偿水平在轨道径向、切向和法向上的变化 Fig. 12 Changes of radial, along-track and cross-track directions in compensation level of empirical acceleration |
| 图选项 |
3.2.3 重叠弧段分析 图 13为2018年1月1日—2018年10日(阶段A)的重叠弧段三维位置均方根(3D Root Mean Square, 3D RMS)在质心位移估计前后的变化。定轨弧长为27 h，采用前一弧段最后3 h与当前弧段的前3 h重叠轨道差异以评价轨道一致性。与未进行质心位移估计的结果进行对比，质心位移估计后，天宫二号卫星轨道重叠弧段3D RMS由4.33 cm降低到2.64 cm，定轨精度较好。
 |
| 图 13 重叠弧段精度变化 Fig. 13 Change of accuracy of overlapping arc |
| 图选项 |
3.2.4 激光校验结果分析 阶段B为三轴对地稳定模式，不适合本体坐标系Xb、Yb方向质心位移估计，但该阶段存在地面激光数据。将阶段A计算的质心位移量应用到这个阶段的精密定轨过程中。分为2种情形对比：一个是不应用质心位移校正结果(利用参考质心定轨，不做质心位移计算)；另一个是应用质心位移校正结果(即利用阶段A计算的质心位移更新参考质心，但不做质心位移计算)。选取太阳高度角在45°以上的激光标准点数据，分别与上述2种情形的基于GNSS计算的精密轨道进行对比。在天宫二号卫星处于三轴对地稳定姿态下(见表 8，阶段B)时，是否应用质心位移校正对定轨结果影响不大，这与之前的分析是一致的。
表 8 三轴对地稳定姿态下，卫星激光测距比对结果 Table 8 SLR comparison results in attitude mode of three-axis Earth-pointing stabilization
| 日期 | 未应用质心位移量的定轨结果与激光测距比较残差/cm | 应用质心位移量的定轨结果与激光测距比较残差/cm |
| 2018-08-20 | 3.14 | 3.01 |
| 2018-08-21 | 3.29 | 3.17 |
| 2018-08-22 | 2.96 | 2.76 |
| 2018-08-23 | 3.03 | 2.95 |
| 2018-08-24 | 3.27 | 3.08 |
| 平均值 | 3.14 | 2.99 |
表选项
阶段C为连续偏航模式，适合做质心位移估计且该阶段存在地面激光数据。将阶段A计算的质心位移量应用到这个阶段的精密定轨过程中。分为2种情形对比：一个是不应用质心位移校正结果(利用参考质心定轨，不做质心位移计算)；另一个是应用质心位移校正结果(即利用阶段A计算的质心位移更新参考质心，但不做质心位移计算)。选取太阳高度角在45°以上的激光标准点数据，分别与上述2种情形的基于GNSS计算的精密轨道进行对比。在天宫二号处于连续偏航姿态下(见表 9，阶段C)时，和未应用质心位移校正的定轨结果与激光测距比较残差相比，应用质心位移校正的定轨结果与激光测距比较残差降低了约0.86 cm，表明应用质心位移校正的定轨结果的轨道精度更高。
表 9 连续偏航姿态下，卫星激光测距比对结果 Table 9 SLR comparison results in attitude mode of yaw steering
| 日期 | 未应用质心位移量的定轨结果与激光测距比较残差/cm | 应用质心位移量的定轨结果与激光测距比较残差/cm |
| 2018-08-27 | 3.45 | 2.73 |
| 2018-08-28 | 3.12 | 2.12 |
| 2018-08-29 | 2.93 | 2.15 |
| 2018-08-30 | 2.82 | 1.87 |
| 平均值 | 3.08 | 2.22 |
表选项
4 结论 本文基于天宫二号卫星的GNSS观测数据，给出了利用简化动力学最小二乘批处理方法来确定天宫二号卫星动力学质心位移的方法。
1) 定性和定量分析的结果表明，基于GNSS测量的定轨结果对处于轨道切向上的质心位移不敏感，而对轨道径向和法向上的质心位移较为敏感。
2) 利用实测数据进行轨道确定的结果表明，当天宫二号卫星处于连续偏航姿态下时，本体坐标系下的Xb和Zb方向的质心偏移估计精度较高，Yb方向质心估计精度相对较差。①由内符合分析可知，在偏航姿态下，与未进行质心位移校正得到的定轨结果相比，进行质心位移校正后得到的相位后验残差标准差降低了0.04 cm，经验加速度在轨道径向、切向和法向上分别降低62%、50%和65%，重叠弧段轨道精度提高1.69 cm。②由外符合分析可知，激光校核的结果表明天宫二号卫星本体坐标系Xb方向质心位移校正对三轴对地稳定姿态下的定轨精度改善不明显，但该方向质心位移校正对偏航姿态的定轨精度有明显改善，与激光数据比对的轨道精度提高0.86 cm。
附录A: 卫星姿态和姿态角定义 卫星的姿态一般通过卫星本体坐标系相对于轨道坐标系(或者惯性坐标系)之间的转动关系来描述，转动的角度即为卫星的姿态角。姿态偏航角、姿态滚动角和姿态俯仰角定义为本体坐标系相对轨道坐标系的3个欧拉角。初始状态下，卫星的本体坐标系ObXbYbZb与轨道坐标系OoXoYoZo重合(红)，经过偏航转动至OX′oY′oZo(蓝)；经过滚动转动至OX′oYbZ′o(紫)，俯仰转动至OXbYbZb，转动顺序为Z-X-Y，即3-1-2转序，如图 A1所示。
 |
| 图 A1 姿态角示意图 Fig. A1 Illustration of attitude angle |
| 图选项 |
在图 A1中，姿态偏航角为ψ，姿态滚动角为φ，姿态俯仰角为θ。卫星本体坐标系坐标轴相对轨道坐标系坐标轴旋转满足右手螺旋法则，即逆时针旋转时角度为正，顺时针方向旋转为负，姿态角单位为(°)。
参考文献
| [1] | MONTENBRUCK O, VAN H T, KROES R, et al. Reduced dynamic orbit determination using GPS code and carrier measurements[J]. Aerospace Science and Technology, 2005, 9(3): 261-271. DOI:10.1016/j.ast.2005.01.003 |
| [2] | 秦建, 郭金运, 孔巧丽, 等. Jason-2卫星星载GPS数据cm级精密定轨[J]. 武汉大学学报(信息科学版), 2014, 39(2): 137-141. QIN J, GUO J Y, KONG Q L, et al. Precise orbit determination of Jason-2 with precision of centimeters based on satellite-borne GPS technique[J]. Geomatics and Information Science of Wuhan University, 2014, 39(2): 137-141. (in Chinese) |
| [3] | LI K, ZHOU X, WANG W, et al. Centimeter-level orbit determination for TG02 spacelab using onboard GNSS data[J]. Sensors, 2018, 18(8): 2671. DOI:10.3390/s18082671 |
| [4] | 张强. 采用GPS与北斗的低轨卫星及其编队精密定轨关键技术研究[D]. 武汉: 武汉大学, 2018: 22-35. ZHANG Q. Research on the key technologies of precise orbit determination for low earth orbit satellites and their formation using GPS and BDS[D]. Wuhan: Wuhan University, 2018: 22-35(in Chinese). |
| [5] | 秦显平. 星载GPS低轨卫星定轨理论及方法研究[J]. 测绘科学与工程, 2010, 30(1): 77-78. QIN X P. Research on precision orbit determination theory and method of low earth orbiter based on GPS technique[J]. Science of Surveying and Mapping, 2010, 30(1): 77-78. (in Chinese) |
| [6] | BRUINSMA S, TAMAGNAN D, BIANCALE R. Atmospheric densities derived from CHAMP/STAR accelerometer observations[J]. Planetary and Space Science, 2004, 52(4): 297-312. DOI:10.1016/j.pss.2003.11.004 |
| [7] | 戴小蕾. 基于平方根信息滤波的GNSS导航卫星实时精密定轨理论与方法[D]. 武汉: 武汉大学, 2016: 18-20. DAI X L. Real-time precise GNSS satellite orbit determination using the SRIF method: Theory and implementation[D]. Wuhan: Wuhan University, 2016: 18-20(in Chinese). |
| [8] | 刘伟, 俞洁, 杨立峰, 等. GEO卫星在轨横向质心快速估算方法[J]. 航天器工程, 2016, 25(5): 39-44. LIU W, YU J, YANG L F, et al. Method of in-orbit lateral centroid fast estimation of GEO satellite[J]. Spacecraft Engineering, 2016, 25(5): 39-44. DOI:10.3969/j.issn.1673-8748.2016.05.006 (in Chinese) |
| [9] | TANYGIN S, WILLIAMS T. Mass property estimation using coasting maneuvers[J]. Journal of Guidance, Control, and Dynamics, 1997, 20(4): 625-632. DOI:10.2514/2.4099 |
| [10] | 郭正勇, 张增安, 汪礼成, 等. 一种基于推力器控制的卫星质心在轨估算方法研究[J]. 上海航天, 2017, 34(5): 76-82. GUO Z Y, ZHANG Z A, WANG L C, et al. Study of on-orbit estimation method of satellite'centroid based on thrust control[J]. Aerospace Shanghai, 2017, 34(5): 76-82. (in Chinese) |
| [11] | 谭沧海, 梁翠娜, 薛宏伟. 伪距定位算法中天线相位中心偏差的修正及误差分析[J]. 现代导航, 2017, 8(5): 328-333. TAN C H, LIANG C N, XUE H W. Analysis and correction of antenna phase center offsets on pseudo-range positioning algorithm[J]. Modern Navigation, 2017, 8(5): 328-333. (in Chinese) |
| [12] | 辛宁, 邱乐德, 张立华, 等. 一种重力卫星质心在轨标定算法[J]. 中国空间科学技术, 2013, 33(4): 9-15. XIN N, QIU L D, ZHANG L H, et al. Study on on-orbit calibration of center of mass for gravity satellite[J]. Chinese Space Science and Technology, 2013, 33(4): 9-15. (in Chinese) |
| [13] | 李洪银, 屈少波, 白彦峥, 等. 静电悬浮加速度计在轨质心位置的最小二乘估计[J]. 地球物理学报, 2017, 60(3): 897-902. LI H Y, QU S B, BAI Y Z, et al. Least squares estimation of in-orbit mass center position of the electrostatic accelerometer[J]. Chinese Journal of Geophysics, 2017, 60(3): 897-902. (in Chinese) |
| [14] | 陈光锋, 唐富荣, 薛大同. 重力卫星在轨质心修正原理[J]. 宇航学报, 2005, 26(5): 567-570. CHEN G F, TANG F R, XUE D T. The trim principle of center of mass of gravity satellite during orbit flight[J]. Journal of Astronautics, 2005, 26(5): 567-570. (in Chinese) |
| [15] | WANG W, LIU J, SHUA L, et al. Precise orbit determination of large-scale spacecraft in low earth orbit: Preliminary results[C]//Proceedings of the 68th International Astronautical Congress (IAC), 2017: 3-12. |
| [16] | SHAO K, GU D, JU B, et al. Analysis of Tiangong-2 orbit determination and prediction using onboard dual-frequency GNSS data[J]. GPS Solutions, 2020, 24(1): 11. DOI:10.1007/s10291-019-0927-y |
| [17] | KROES R. Precise relative positioning of formation flying spacecraft using GPS[D]. The Netherland: Delft University of Technology, 2006: 56-63. |
| [18] | WANG W, GAO Y. Effective empirical acceleration modeling and its application to enhanced-accuracy orbit prediction[J]. Transactions of the Japan Society for Aeronautical and Space Sciences, Aerospace Technology Japan, 2016, 14(30): 39-45. |
| [19] | JIA T U, GU D F, YI W U, et al. Phase residual estimations for PCVs of spaceborne GPS receiver antenna and their impacts on precise orbit determination of GRACE satellites[J]. Chinese Journal of Aeronautics, 2012, 25(4): 631-639. |
| [20] | RIM H J, YOON S, SCHUTZ B, et al. Effect of center of mass position error on icesat precision orbit determination[C]//AIAA/AAS Astrodynamics Specialist Conference. Reston: AIAA, 2010: 77-79. |
| [21] | CHOI K R. Jason-1 precision orbit determination using GPS combined with SLR and DORIS tracking data[D]. Austin: The University of Texas at Austin, 2003: 43-47. |
