基于光学测角数据的风云四号同步轨道卫星精密定轨

刘思语，黄勇，毛银盾，贾耀红，鲁文强,4，黄乘利,3，郑景辉，杨鹏

1. 中国科学院上海天文台，上海 200030 2. 中国科学院大学 天文与空间科学学院，北京 100049 3. 上海科技大学 物质学院，上海 201210 4. 国家卫星气象中心，北京 100081

1 引言

风云四号气象卫星于2016年12月11日由西昌卫星发射中心成功发射系列第一颗卫星(以下称FY-4A)并于2017年9月25日正式交付用户使用，是中国第二代静止轨道气象卫星，用于在地球同步轨道带上进行高效、稳定、连续的气象观测[1]。

FY-4A目前采用无线电测距技术实现对卫星的精密定轨，岁差章动模型参考IAU1976，N体摄动模型参考DE405历表，采用数值求解二阶微分方程组的KSG定步长线性多步积分器求解动力学方程，重叠1 h轨道的内符合精度优于20 m。FY-4A卫星定点于东经99°，地面五个跟踪站(北京，佳木斯，广州，乌鲁木齐，腾冲)同时进行观测，数据采样率为1 s，观测模式为每小时或每半小时连续观测10 min，每日轨控一次，获取得到S波段与L波段的测距信号，对测距数据经过系统差标校后，利用双频改正方法进行电离层延迟改正，每个测站配备自动气象站，为大气折射修正提供气压、温度、湿度等参数。

对于大部分地球同步轨道卫星(GEO)定轨，国内现阶段除了光学测角手段之外，还包括无线电测距、激光测距(SLR)、雷达等技术[2]。宋叶志等基于多台站双频双程测距数据计算风云四号气象卫星轨道，通过轨道重叠分析，在非变轨时期获得20 m的定位精度[3]；刘凯等采用一发多收的转发式测定轨方法测量北斗导航卫星的轨道，精度达到10 m[4]；卿芸等利用星载GNSS和SLR数据对北斗C08、C10地球同步轨道卫星进行联合定轨，定轨精度优于0.1 m[5]；Zongbo HUYAN等基于星间链路测量对北斗导航卫星进行定轨，48 h内轨道预报精度优于2.4 m[6]；张志斌等使用口径1 m的“微型天线”在上海、都匀和乌鲁木齐三个测站组网，监测“亚太六号”地球同步轨道卫星，轨道精度可达百米量级[7]；Yong Huang等利用地基VLBI数据对地球同步轨道卫星定轨，VLBI时延精度可达3.6 ns，定轨精度优于20 m[8]。

地基光学观测手段相比上述测量技术具有运行成本低、结果更直观等优势，并且作为一种无源观测手段，在探测非合作目标时拥有更广泛的应用范围[9]。目前地面的光学观测设备主要成像终端为CCD或CMOS相机，通过照相天体测量的方式获取目标星的位置。照相天体测量技术需要预知参考星的理想坐标、参考星和目标星在CCD或CMOS像面上的量度坐标，利用每一组参考星的理想坐标和量度坐标构建底片模型，再通过底片模型和目标星的量度坐标反解出目标星的理想坐标和天球坐标。GEO相对地面测站几乎静止不动，受到地球自转影响，地球同步轨道卫星与背景恒星存在运动差异(大约15＂/s)，因此相机采用凝视模式观测并获得GEO圆形星象的同时，背景恒星的星象图像会拖长，无法精确确定参考星的位置中心。于涌等利用CCD漂移扫描技术解决星象拖长的问题[10]。望远镜对目标精密跟踪的情况下，星像固定于CCD视场。星像所产生的电荷在同一像素区域内逐渐积累。当积累的能量明显超过周围的背景噪声时，信噪比高的星象会显现出来。CCD漂移扫描技术实现电荷在像素间的行转移，当电荷转移速度与目标天体视运动速度在CCD像面上的投影同步且目标运动方向与电荷转移方向一致时，能够在固定望远镜指向下实现对运动目标的电荷跟踪，得到圆形星象。严丹等利用相似的方式观测某颗GEO和木卫六，分别在赤经赤纬上取得0.32＂、0.28＂和0.100＂、0.105＂的观测精度[11]。2011年西班牙皇家天文台基于光学测角数据对H1D GEO卫星精密定轨，取得100 m的外符合精度[12]。

本文利用光学测角数据对FY-4A数据进行定轨分析，分析了仅利用测角数据定轨以及利用测距测角数据联合定轨的精度，并评估了仅利用测距数据的实时轨道精度。

2 观测实验

2019年12月27～31日，上海天文台GEO全球联合观测网50 cm望远镜(以下称50 cm望远镜)对FY-4A进行了光学观测，观测样例见图1和图2。利用CCD相机观测得到FY-4A的测角数据，直接测量量为目标星的量度坐标，再转换为定轨所需要的赤经赤纬或高度角与方位角。

图1 利用CCD漂移扫描技术观测恒星背景实例Fig.1 An example of CCD drift scanning technology for sidereal stars

50 cm望远镜安装了支持时延积分功能的Alta-F9000相机和GPS/北斗双模时统，具备CCD漂移扫描功能具体参数见表1。2019年12月27～31日该望远镜对FY-4A进行观测，CCD相机每一轮观测拍摄7帧图像，首末2帧以漂移扫描模式拍摄背景恒星，中间5帧以凝视模式跟踪FY-4A。相机快门开启时时间系统记录开始时刻，快门关闭时记录结束时刻。每2帧图像之间的时间间隔约为6.7 s，其中5 s为曝光时间，1.7 s为CCD电荷读取时间。5 s的曝光时间能保证望远镜可以观测15星等的GEO目标。光学测角要求测站处于阴影且被测目标处于可见光照射，保证良好天光背景，因此观测时间多集中于天气状况良好的晴朗夜间[13]，具体观测弧段见表2。

图2 利用凝视模式观测GEO卫星实例Fig.2 An example of stare-mode for GEO satellites

表1 望远镜参数及测站位置

表2 光学测角弧段

初始图像经过中值滤波去除高频信号，保留低频信号，抑制背景的不均匀。再通过连通域星象法检测图像中星象，若某星象的信噪比小于设置门限，则去除该星象。同时基于二维修正矩定心法消除背景噪声的影响，计算星象中心。最后利用相邻帧图像差分和航迹关联法识别GEO卫星[14]。

3 测量模型

3.1 照片天体测量与底片模型

图3 赤经赤纬与理想坐标几何关系Fig.3 The geometric relationship between ideal coordinate and ascension,declination

照相天体测量技术基于底片上目标星的量度坐标推算其理想坐标，并计算其赤经赤纬，赤经赤纬与理想坐标的几何关系见图3。根据观测结果所计算的赤经和赤纬均在站心赤道坐标系下表示。

在处于焦平面上的照相底片取一点C，以C为原点，以焦平面为坐标平面构建笛卡尔坐标系C-ξη，η为赤纬圈在底片上的投影，以此坐标系为理想坐标系。CO轴为望远镜的光轴，C'(α0,δ0)为CO轴延长线与天球的交点,其中α0和δ0为该点的赤经和赤纬，|CO|为焦距F。现有一颗目标星在焦平面上的投影点为S,其对应的在天球上的点位为S'，P为北极所在位置。若设焦距F为单位变量，在三角形OSC中则有：

(1)

式中：φ为OC和OS夹角；θ为CS与η轴夹角。在球面三角形C'S'P中，存在以下转换关系：

(2)

联立(1)(2)有：

(3)

量度坐标系设立于CCD图像，其坐标平面为CCD图像平面，即CCD芯片所在平面。第一轴与第二轴为矩形像元的两条边，第三轴垂直于坐标平面。理想坐标系的坐标平面即为望远镜焦平面，X轴和Y轴分别平行于赤经圈与赤纬圈，Z轴垂直于坐标平面。设备安装中无法保证CCD芯片与焦平面平行，导致整套光学系统不共轴。理想坐标系的原点为望远镜的光心，量度坐标系的原点为相关量度设备在照相底片上的位置中心，二者无法完全一致。不同的望远镜焦距和底片尺寸会导致理想坐标系和量度坐标系的坐标轴比例尺无法严格相等。CCD像元由于制作工艺的因素可能导致像元不是严格的矩形，量度坐标系的两轴无法严格正交，可能与理想坐标系的两轴不平行，除此之外量度坐标系受到光学系统像差等影响，因此通过多项式建立量度坐标(x,y)和理想坐标(ξ,η)的关系[15]：

(4)

底片模型参数aij和bij通过背景恒星的理想坐标和量度坐标建立观测方程，基于最小二乘法求解获得。由(4)可建立赤经赤纬与量度坐标系的转换关系(5)[16]：

(5)

本文定轨所使用观测资料为高度角与方位角，赤经赤纬与高度角方位角存在以下转换关系，A和E为目标星的方位角和高度角，S为测站所在恒星时。

3.2 高度角与方位角

方位角和高度角是基于测站坐标系下的概念，测站坐标系的原点位于观测站站心，即测量设备跟踪天线的旋转中心。站心所在位置的地平面为基本平面，主方向由站心指向正北方向。笛卡尔系下的测站坐标系，YS轴在基本平面内指向东方，XS轴指向主方向，ZS轴垂直于基本平面指向上方。

现设测站坐标系下有一目标位置矢量为rtp(xtp,ytp,ztp),方位角A表现为目标位置矢量在基本平面上的投影与Xs轴的夹角[17]：

方位角A对卫星状态的偏导数：

(6)

式中：r为卫星在J2000.0天球坐标系中的位矢；rb为卫星在地固坐标系中的位矢。式(6)具体展开为：

(7)

式中：HG为J2000.0天球坐标系到地固坐标系的转换矩阵，M为地固坐标系到测站坐标系的转换矩阵。HG包括岁差、章动、自转和极移的修正，λ和φ分别为站心的经度和大地纬度,式(7)中M为：

高度角E表现为目标位置矢量与Zs轴的夹角：

高度角E对卫星状态参数的偏导数：

3.3 测距

测距型数据获取的是测站与卫星相对距离关系，可表示为：

式中：ρ为测站相位中心和星载天线相位中心的距离；ρu为上行信号传播距离；ρd为下行信号传播距离。ρu和ρd具体可表示为：

式中：t为测站接收下行信号的时刻；Δt1为上行光行时；Δt2为下行光行时。r(t-Δt2)为卫星发送下行信号时，卫星的位置矢量。R(t-Δt1-Δt2)为测站发送上行信号时，测站的位置矢量。R(t)为测站接受下行信号时的位置矢量。Δρ是误差修正，包括相对论修正、电离层修正、大气折射修正等。εu和εd分别为上行信号传播距离和下行信号传播距离的测量误差。上行光行时Δt1和下行光行时Δt2无法直接获得，可先迭代计算下行光行时，推导得到卫星发送下行信号的时刻t-Δt2，再经过一次计算上行光行时，推导得到发送上行信号的时刻t-Δt1-Δt2。具体的迭代过程如下：

3.4 定轨原理

在对FY-4A进行光学观测的过程中。假设在ti时刻获得一组观测量Yi:

(8)

设x(t)=X(t)-X*(t)，由(8)式得：

4 数据分析及小结

4.1 测角数据定轨及精度分析

基于50 cm望远镜的测角数据，对2019年12月28～31日的FY-4A轨道进行定轨，FY-4A每日17～18时进行变轨，本文仅分析变轨前后的轨道。统计多天的定轨结果，高度角残差rms为0.25＂，方位角残差rms为0.45＂，具体定轨残差见表3。

表3 高度角与方位角拟合残差

将有测角数据的定轨弧段和参考轨道重叠比较位置差异，参考轨道基于FY-4A测距数据精密定轨获得，测距数据由北京、佳木斯、乌鲁木齐、广州以及腾冲的测站提供，采样率为一秒一个点。在有测角数据的弧段，轨道精度总体优于50 m。具体定轨精度见表4。

表4 存在测角数据弧段的轨道精度

在定轨弧段内，包含测角数据弧段的定轨精度优于其它弧段，结果见图4。

T1、T2-Period with angel-data图4 12月29～30日测角数据定轨精度Fig.4 Precision of orbit determination with optical angle measurement data from 29th December to 30th December

以本次动量轮卸载后到下次动量轮卸载前的弧段为一个定轨周期，在一个定轨周期内，随着测角数据的增加，定轨结果与参考轨道的位置差异会呈现下降的趋势，具体结果见图5。

图5 3个单天的测角数据定轨精度Fig.5 Precision of orbit determination with optical angle measurement data of three single day

以FY-4A运控轨道为参考轨道，选择12月29日和12月30日有测角数据的弧段，计算测角残差O-C，结果见图6、图7。高度角与方位角残差精度在绝大多数的弧段优于0.5＂，高度角平均O-C为0.23＂，方位角平均O-C为0.34＂，具体结果见图6和图7。

图6 12月29日测角残差O-CFig.6 Residual O-C of orbit determination for FY-4A on 29th December

图7 12月30日测角残差O-CFig.7 Residual O-C of orbit determination for FY-4A on 30th December

基于测距数据的轨道精度约为30 m，本次实验取得的测角残差O-C约为0.2＂，故无法利用测角残差评估本次定轨结果。

4.2 测角和测距联合数据定轨及精度分析

基于测距数据和50 cm望远镜的测角数据对2019年12月29日和12月30日的FY-4A进行联合定轨。利用12月28日18时-12月29日6时的测距和测角数据(12 h定轨弧段)，以及12月29日5时-12月29日17时的测距和测角数据分别定轨(12 h定轨弧段)，比较12月29日5时-6时的重叠轨道，轨道精度约为15 m。利用12月29日18时-12月30日6时的测距和测角数据(12 h定轨弧段)，以及12月30日5时-12月30日17时的测距和测角数据分别定轨(12 h定轨弧段)，比较12月30日5时-6时的重叠轨道，轨道精度约为13 m，具体定轨结果见图8和图9。

图8 12月29日1 h轨道重叠精度Fig.8 Precision of orbit determination of 1 h overlap on 29th December

图9 12月30日1 h轨道重叠精度Fig.9 Precision of orbit determination of 1 h overlap on 30th December

以测距+测角数据定轨结果为参考轨道，对12月29日和12月30日运控测距轨道进行精度评估，比较位置差异。评估结果表明，在每日卫星机动(机动时间约为UTC17:30)之后，控后半小时定轨预报半小时轨道精度约为30～45 m，此后随着测距数据的增加，定轨精度稳步提升，控后6 h至下次机动，轨道精度稳定在20 m以内，具体结果见图10和图11。

图10 12月29日测距和测角联合定轨精度Fig.10 Precision of orbit determination with optical angel measurement data and ranging data for FY-4A on 29th December

图11 12月30日测距和测角联合定轨精度Fig.11 Precision of orbit determination with optical angel measurement data and ranging data for FY-4A on 30th December

5 结论

基于上海天文台GEO全球联合观测网50 cm望远镜的光学测角数据对FY-4A进行定轨分析，剔除动量轮卸载变轨的弧段，单天测角数据定轨的高度角残差rms优于0.25＂，方位角残差rms优于0.45＂。与基于测距数据提供的运控轨道相比，在有测角数据的弧度，轨道精度优于50 m。在单天定轨弧段内，随着测角数据量的增加，轨道的精度呈现上升趋势。利用测角和测距数据对FY-4A进行联合定轨，以多组12 h的数据为定轨弧段，两组数据存在1 h的相同弧段，重叠弧段的精度优于15 m。以一个单天的数据为定轨弧段，高度角平均O-C为0.23＂，方位角平均O-C为0.34＂。在FY-4A动量轮卸载后，控制后半小时轨道精度在30-45 m，随着测量数据的增加，定轨精度稳定提升，控后6 h至下次变轨前，轨道精度优于20 m。