覆盖分析模块是航天动力学软件^[1]中的重要内容，可为遥感、侦察、预警等卫星的任务规划提供基础数据，而星下点的精确和实时求解能为其提供有力支撑^[2]。传统的星下点计算方法较为成熟^[3-4]，但地球在空间中的运动很复杂，受岁差、章动和极移等因素的影响^[5]，应用该方法时需要将轨道六要素转换到真赤道春分点坐标系中，如果直接在卫星工作的J2000.0平赤道坐标系中使用会产生误差。本文结合仿真实例给出该误差的大小，结果表明，进行坐标转换是必要的。

1 传统的星下点计算方法

将地球视为圆球体，利用球面三角形知识，通过以下公式得到t时刻卫星的星下点坐标：

$ \left\{\begin{aligned} \varphi &=\arcsin (\sin i \sin u) \\ \lambda &=\mathit{\Omega }+\arctan (\cos i \tan u)-\\ & \theta_{g 0}-\omega_{c}\left(t-t_{0}\right) \\ u &=\omega+f \end{aligned}\right. $

(1)

式中，φ为星下点纬度，λ为星下点经度，i为轨道倾角，u为t时刻卫星与升交点的角距，Ω为升交点赤经，θ_g0为t₀时刻对应的格林尼治恒星时，ω_c为地球自转角速度，ω为近地点幅角，f为t时刻的真近点角。

考虑J₂摄动时的星下点计算公式为：

$ \left\{\begin{aligned} \varphi &=\arcsin (\sin i \sin u) \\ \lambda &=\mathit{\Omega }_{0}+\arctan (\cos i \tan u)-\\ & \theta_{g0}-\left(\omega_{c}-\dot{\mathit{\Omega }}\right)\left(t-t_{0}\right) \end{aligned}\right. $

(2)

式中，$\mathit{\dot \Omega } $为受到J₂摄动作用引起的升交点赤经变化率。可以看出，传统的星下点计算需要在真赤道春分点坐标系中进行，否则由于地球运动中的岁差、章动及极移等现象，会使星下点求解产生较大的偏差。

2 基于GCRS-ITRS变换关系的星下点计算方法

通常使用轨道六要素来表示卫星在空间中的运动状态^[6]，将轨道六要素投影到近焦点坐标系后，利用3次旋转变换，得到卫星的位置矢量r_J2000.0。

2.1 GCRS和ITRS坐标变换

J2000.0平赤道坐标系与地心天球坐标系(GCRS)存在参考架偏差，通过公式完成相互转换：

$ \boldsymbol{r}_{\mathrm{GCRS}}=\boldsymbol{B}^{-1} \cdot \boldsymbol{r}_{\mathrm{J} 2000.0} $

(3)

$ \boldsymbol{B}=\boldsymbol{R}_{1}\left(-\eta_{0}\right) \boldsymbol{R}_{2}\left(\xi_{0}\right) \boldsymbol{R}_{3}\left(\mathrm{d} \alpha_{0}\right) $

(4)

$ \left\{\begin{array}{l}{\xi_{0}=-16.617 \pm 0.010 \mathrm{ mas}} \\ {\eta_{0}=-6.819 \pm 0.010 \mathrm{ mas}} \\ {\mathrm{d} \alpha_{0}=-14.6 \pm 0.5 \mathrm{ mas}}\end{array}\right. $

(5)

式中，r_J2000.0和r_GCRS分别为J2000.0平赤道坐标系与地心天球坐标系下的位置矢量，R₁、R₂、R₃为初等旋转矩阵。

刘林等^[7]认为参考架偏差很小，通常无需考虑，因此本文不区分J2000.0平赤道坐标系与GCRS坐标系。GCRS-ITRS变换即地心天球坐标系(接近J2000.0平赤道坐标系)与国际地球参考系(地固坐标系)之间的相互转换，其精度决定本文星下点计算结果的可靠性。

根据IAU决议，坐标变换有基于春分点和格林尼治恒星时(GST)的经典变换及基于天球中间零点(CIO)和地球自转角(ERA)的变换2种方法，其在微角秒的精度上是一致的^[5]。本文采用基于春分点和格林尼治恒星时(GST)的经典变换，使用IERS2010规范^[8]。

2.1.1 基于春分点和格林尼治恒星时坐标转换

假设[GCRS]和[ITRS]中卫星位置矢量分别为x₁、x_2、x₃和y₁、y₂、y₃，则：

$ \left[\begin{array}{l}{\boldsymbol{y}_{1}} \\ {\boldsymbol{y}_{2}} \\ {\boldsymbol{y}_{3}}\end{array}\right]=[\mathrm{ITRS}]^{-1}[\mathrm{GCRS}]\left[\begin{array}{l}{\boldsymbol{x}_{1}} \\ {\boldsymbol{x}_{2}} \\ {\boldsymbol{x}_{3}}\end{array}\right] $

(6)

$ [\mathrm{ITRS}]^{-1}[\mathrm{GCRS}]=\boldsymbol{W}(t) \cdot \boldsymbol{R}(t) \cdot \boldsymbol{M}(t) $

(7)

式中，W(t)为极移矩阵，R(t)为地球自转矩阵，M(t)为岁差-章动矩阵(其形式较为复杂，详细推导及结果见文献[5])。为方便表述，[ITRS]^-1[GCRS]为GCRS到ITRS的旋转矩阵，下文中类似的表述代表同样的意义。

2.1.2 利用SOFA程序库完成GCRS-ITRS变换

IAU基本天文学标准(SOFA)程序库提供的IAU 2006/2000A变换矩阵算法程序库，可用来完成GCRS和ITRS之间的坐标变换。在使用前，需要根据IERS(国际地球自转和参考系服务机构)提供的公告，将地球自转参数UT1-UTC、极移分量x_p和y_p导入程序，这3个值分别在10^-1 s、10^-1″和10^-1″量级，甚至更小，但若看成定值将引起数10 m的转换误差。图 1给出SOFA标准库的调用流程，其具体过程如下：

1) 给定TT时间的儒略日JD_TT，调用函数iauNut06a和iauPn06，得到[CIRS]^-1[GCRS]旋转矩阵；

2) 结合IERS公告数据，通过函数iauUtcut1得到对应的JD_UT1，调用函数iauGST06得到格林尼治视恒星时GST，给出[TIRS]^-1[CIRS]旋转矩阵;

3) 调用函数iauSp00和iauPom00得到极移矩阵[ITRS]^-1[TIRS]，最后调用函数iauC2teqx，得到[ITRS]^-1[GCRS]旋转矩阵，结合式(6)完成坐标转换。

STK是航天领域专业仿真软件，内含高精度数学模型，可提供高精度的转换结果^[9]。通过STK中的仿真数据检验本文的转换结果，设定检验仿真时间为2013-11-30 03：39：1.818，给定天球坐标系下位置矢量(-5 420 499.167, 7 425 896.670, 367 855.733)m。本文及通过STK得到的地固系下位置矢量结果见表 1，结果显示，两者之间的偏差在3 m以内，满足精度要求，证明了模型的准确性。

2.2 空间直角坐标与大地坐标的转换

通过GCRS-ITRS变换，得到卫星在地固坐标系下的空间直角坐标(X, Y, Z)，该点亦可表示在大地坐标系中，记为(B, L, H)，其中B为大地纬度，L为大地经度，H为大地高，(B, L)即为星下点坐标。

2.2.1 大地纬度的计算

采用经典的迭代解法计算大地纬度B^[10]：

$ \begin{array}{c}{\tan B_{i+1}=} \\ {\frac{1}{\sqrt{X^{2}+Y^{2}}}\left(Z+\frac{a \cdot e^{2} \tan B_{i}}{\sqrt{1+\left(1-e^{2}\right) \tan ^{2} B_{i}}}\right)}\end{array} $

(8)

式中，a为地球长半径，e为地球第一偏心率。取迭代初值$\tan B_{1}=\frac{Z}{\sqrt{X^{2}+Y^{2}}}\left(1+e^{2}\right) $，通常2~3次即可达到10^-4″的精度^[10]。

2.2.2 大地经度的计算

在大地坐标系中，大地经度L可表示为：

$ \tan L=\frac{Y}{X} $

(9)

$ L=\arctan \left(\frac{Y}{X}\right)+[1-\operatorname{sign}(X)] \cdot \frac{\pi}{2} $

(10)

由式(10)可以确定L的大小及象限，其中sign为符号函数。

3 算法仿真结果

结合仿真实例对2种方法得到的星下点进行分析。设置仿真初始时刻为2013-11-30 03 :39 :1.818，对应的轨道六要素(J2000.0平赤道坐标系下)见表 2，通过仿真得到1 200 min内的星下点坐标序列。其中，采用基于GCRS-ITRS变换方法得到的部分星下点坐标序列见表 3，直接在J2000.0平赤道坐标系下采用传统方法得到的部分星下点坐标序列见表 4。

图 2为1 200 min内使用两种方法得到的星下点轨迹，图 3和4分别为1 200 min内使用两种方法得到的纬度之差与经度之差。由图可知，两种方法得到的星下点之差近似呈周期性变化，其中纬度之差最大为0.143 9°，经度之差最大为0.214 4°，而0.2°的误差将引起地面星下点22.263 4 km的偏差。因此，在J2000.0平赤道坐标系下，直接采用传统方法会产生较大误差，必须进行相应的坐标转换。

4 结语

本文针对J2000.0平赤道坐标系下星下点的求解进行探讨，综合考虑地球存在的岁差、章动和极移等因素，基于GCRS-ITRS变换关系，使用成熟的SOFA标准库，在大地椭球模型的基础上得到了精度较高的星下点坐标。结合仿真实例，给出在J2000.0平赤道坐标系下直接使用传统方法进行求解的误差，其中纬度方向最大误差为0.143 9°，经度方向最大误差为0.214 4°，而0.2°的误差将引起地面星下点22.263 4 km的偏差，所以必须进行相应的坐标转换。本文实现了J2000.0平赤道坐标系下星下点的实时高精度求解，可为航天动力学软件的编写、覆盖分析等研究提供一定参考。