1 引言
随着我国航天事业的飞速发展，越来越多的卫星进入太空，在导航、定位、通信、遥感等领域发挥着重要的作用。近年来各大国对空间资源的竞争愈发激烈，使得轨道资源愈发紧张，尤其是拥有重要战略地位的地球同步轨道资源。因此，为更加合理地规划轨道资源、避免卫星碰撞、保障太空安全，各国都对空间目标的精密定轨提出了更高的要求，除了要对已知的合作目标进行精密跟踪，还需要对非合作的空间目标进行监测。
传统的地基监测系统在时间、空间的监测上均存在监视盲区，无法满足我国当前对空间态势感知能力的高要求。自20世纪末，美国就致力于天基监测系统的建设，试图将地基系统与天基观测相结合，实现对天区的全面覆盖。相比美国，中国对天基系统的建设起步较晚，且在地基台站的建设上同时受到政治、地理等因素的限制，在空间覆盖上还需要其他的手段来进行补充。海基观测可以在陆地以外的海域进行，包括我国的领海以及公海，所以海基观测在一定程度上可以弥补地基观测的不足，我国天基观测系统建成后，海基观测依然可以作为一个重要组成部分，为航天测控提供有力的支持。
目前地基测控研究已经有了诸多成果，于涌等基于上海天文台佘山站1.56 m望远镜，利用单站CCD漂移扫描技术对北斗同步卫星进行定轨，单圈精度为数千米量级，多圈精度可达50 m[1]。宋叶志等利用多台站，采用双频双程测距模式，对风云四号进行定轨，非变轨期间精度优于20 m[2]。YongHuang等利用VLBI测轨数据以及C波段转发式测距数据对GEO进行定轨分析，VLBI时延测量精度为3.6 ns，定轨精度约为20 m[3]。美国从1996年发射第一个携带天基光学传感器的低轨卫星，论证了单个平台对高轨带监视的可行性[4]开始，在此后的20多年的时间里，不断发射卫星，构建SSBS天基太空监视系统[5-7]。国内的天基监视系统尚未建设完成，周海银等设计SBSS系统覆盖区域的简易算法，对天基监视系统进行仿真，得到低轨目标精度10 m，高轨目标精度500 m左右的结论[8]。周庆勇等基于天基光学监视平台的空间目标可视性以及轨道改进算法，对多平台联合定轨进行仿真[9]。海基测控在国内起步较晚，康德勇等在2005年提出，海上移动站精密测轨的方法，并对低轨卫星进行了定轨[10]。在2010年，分析了船位误差对船载外测数据的影响，随着距离的增加对定轨结果的影响逐渐减弱[11]。朱伟康等建立测速数据的船姿船速修正模型等一系列优化途径，使得初轨计算半长轴外符合精度提高3倍[12]。陈红英等使用最小二乘法解算双站目标坐标的数据处理方法，提高了弹道状态估计精度[13]。倪兴等提出在航天测控数据处理中使用简化的大气折射公式，得到的轨道精度满足要求，同时满足实时数据处理数据的计算速度[14]。
目前航天测量船的测控手段多为雷达测控，且观测对象多为低轨目标，应用于高轨光学监测的相对较少。本研究使用海基单站测量设备采集到的光学测角数据，利用数值方法，在观测弧段较短且不连续的情况下，对GEO卫星进行轨道改进。
2 海基天文光学测量方法
海基天文光学测量方法是利用CCD成像中的空间目标和背景恒星的相对位置，通过计算其理想坐标，得到其赤经赤纬。
以原点C建立直角坐标系C-ξη,C是CCD底片上的一点，η为赤纬圈投影，ξ与η垂直，为赤经圈投影，正方向取赤经赤纬增加的方向，该坐标系为理想坐标系。望远镜视场平面与天球的交点为(α0,δ0)。空间目标在该坐标系中的投影坐标(ξ,η)与其赤道坐标(α,δ)的转换公式为：
(1)
理想坐标系和量度坐标系由于设备参数的原因在中心原点、坐标轴比例以及坐标轴平行等方面存在差异，可以通过多项式建立二者之间的关联。
底片模型参数是通过理想坐标和量度坐标之间建立方程，利用最小二乘原理求得。利用底片模型参数，可以得到观测目标在理想坐标系下的坐标，进而得到其观测坐标，即赤经、赤纬
(2)
3 海基测定轨原理
3.1 海基测定轨涉及的坐标系统
用于定轨的观测资料是在惯导地平坐标系Ol-XlYlZl(舰船地理坐标系)下描述的，其坐标原点是舰船摇摆运动中心在舰船甲板上的投影，通过惯导三轴(北向、东向、方位)定义基准面。OlXl指向正东，OlYl指向正北，OlZl指向由OlXl,OlYl,按右手螺旋规则确定。此外还有惯导甲板坐标系Od-XdYdZd(舰船甲板坐标系)，其中心与Od与Ol重合，OdYd指向与艏艉线一致，OdXd指向垂直OdYd，平行于甲板面，以右舷为正，OdZd指向由OdXd和OdYd按照右手螺旋规则确定；原始的观测资料是在测量坐标系Or-XrYrZr(设备直角坐标系)下描述，该坐标系的中心位于设备中心，三轴指向与Od-XdYdZd一致[15]。
将原始观测资料从Or-XrYrZr下转换到Ol-XlYlZl需要借助船载的形变测量设备和惯导测量设备。假设某个时刻，形变测量设备的观测值为(xdelt,ydelt,zdelt)，该值是设备中心Or到Ol(Od)之间几何矢量的形变；惯导设备的观测量为(K,φ,θ)，该值分别是舰船的艏摇角、纵摇角和横摇角。
图1 船载坐标系示意图
Fig.1 Schematic diagram of shipborne coordinate system
由Or-XrYrZr转换到Od-XdYdZd，考虑船体形变，进行如下变换：
(3)
再由Od-XdYdZd转换到Ol-XlYlZl，考虑船体姿态，进行如下变换：
(4)
其中旋转矩阵RM具体表达如下：
3.2 光学测角精密定轨方法
空间目标在围绕地球进行运行的过程中会受到许多力的作用。这些力大致可以分为两类，一类是保守力，诸如中心天体的引力，其他大型天体的引力等。另一类是耗散力，则系统能量损耗，对于卫星而言，主要是大气阻力、太阳光压、地球红外辐射等。
在惯性系中，空间目标的运动方程可以表示为：
(5)
式中：r为空间目标在某时刻t的位置矢量；
为空间目标的加速度矢量，方程右端为作用在卫星单位质量上的力的总和。其中， f0为中心引力， fε为各种摄动加速度。
为参考时刻t0时空间目标的状态矢量。一般参考时刻的状态量无法精确获得，只能对其参考值不断地观测精化，以取得高精度的初始状态。
观测量为赤经赤纬，其观测方程为：
(6)
Δx=x-xs
Δy=y-ys
Δz=z-zs
(7)
式中：(x,y,z),(xs,ys,zs)分别为目标卫星和观测平台在天球参考系下的坐标； ξα,ξδ为大气折射改正；ηα,ηδ为周年光行差改正；να,νδ为周日光行差改正。
设在ti时刻有一组观测量Yi：
(8)
式中：
为空间目标在ti时刻的状态矢量；
为空间目标的速度矢量；p为待估参数，包括大气阻力系数、太阳光压系数等；F(X,t)是和状态量相关的力学参数以及速度；G(Xi,ti)是ti时刻观测量的真值； Vi为测量误差。由于观测方程，即式(6)为非线性方程，可以将其线性化，假设某一时刻参考状态为X*(t)，在参考状态X*(t)展开，并令x(t)=X(t)-X*(t)，略去二次以上高阶项有：
(9)
对式(9)进行变形可得：
(10)
式中：y为观测值与理论值之差，即残差(O-C，Observation-Calculation)；x0为待估状态量X0的改正值；H矩阵中的
是测量矩阵，具体表达如下：
(11)
式中： ρ为目标到台站的距离；
为状态转移矩阵，具体公式如下：
(12)
式中，m为待估参数的个数，n=m+6。
若令
(13)
可得到观测方程，即条件方程：
y=Hx0+V
(14)
卫星精密定轨的工作即在上述基础上进行，由大量观测采样数据Yj求解条件方程，可给出待估状态量X0的改正值
从而得到改进的历元状态量
定轨过程是一个迭代过程，将每次得到
作为估计量
重复前面的过程，直至满足精度为止。
利用协方差分析理论可以研究不同误差源对估值精度的影响，将观测方程表达如下：
y=Hxx+Hcc+V
(15)
式中：x为求解参数；c为考察参数；Hx和Hc分别待估参数和考察参数的偏导矩阵。定义灵敏矩阵为：
(16)
(17)
式中：P为计算协方差；W为权阵。另外定义摄动矩阵Γ=Sπ1/2，其意义是考察参数有1σ的先验误差时，该误差对待估参数的影响。摄动矩阵是评价解算的外符合精度的重要统计量。其中：
(18)
利用灵敏度矩阵和摄动矩阵，考察协方差阵可以表述为：
Pc=P+SπST=P+ΓΓT
(19)
因为ΓΓT为正定矩阵，式(19)意味着考察参数的误差会降低待估参数的精度，并且协方差分析可定量的计算这种影响的大小。
4 数据处理与分析
4.1 仿真分析
对海基单站光学跟踪GEO卫星定轨进行仿真模拟。初始历元为2021年6月5日14时，台站(船载设备)的初始位置以及目标GEO卫星的初始星历如表1所示(船载设备的位置以经纬度、高程表示，GEO卫星的初轨以无奇点轨道根数表示)。表2所示为台站的初始位置。
表1 GEO目标的初始星历
Table 1 Initial ephemeris of GEO
a/kmhkGEO目标42 165.2730.000 4-0.002ξ/10-5η/10-5λ/(°)GEO目标-2.046-9.788235.893
无奇点轨道根数与开普勒轨道根数的关系式如下：
(20)
式中：a为轨道半长轴；i为轨道倾角；e为轨道偏心率；ω为近地点幅角；Ω为升交点赤经；M为平近地点角。
表2 台站的初始位置
Table 2 Initial position of station
Lat/(°)Long/(°)Height/m台站17.970133.21619.560
仿真的观测噪声为3 ″，该噪声为随机误差，服从高斯分布，采样率为3 s，模拟的观测弧段为2021年6月5日14∶00—15∶00以及16∶00—17∶30，观测不连续。
4.1.1 台站位置噪声为50 m
将台站位置在地固系下3个方向上的分量依次加上50 m噪声，使用时间序列化的台站位置和模拟的观测数据进行定轨，定轨得到的残差如图2所示。
图2 位置噪声为50 m的定轨残差示意图
Fig.2 POD residuals with noise amplitude of 50 meters
改进后的轨道和仿真参考轨道在观测弧段内的重叠精度如图3所示。
图3 改进轨道与仿真轨道曲线
Fig.3 Overlap of POD and simulation orbit
4.1.2 台站位置噪声为10 m
将台站位置在地固系下3个方向上的分量依次加上10 m噪声，同样使用时间序列化的台站位置和模拟的观测数据进行定轨，定轨得到的残差如图4所示。
图4 位置噪声为10 m的定轨残差示意图
Fig.4 POD residuals with noise amplitude of 10 meters
改进后的轨道和参考轨道在观测弧段内的重叠精度如图5所示。
对比台站位置噪声振幅50 m和10 m的定轨结果可以看出：不同的噪声振幅对于定轨残差的影响较小，这是由较小的台站位置偏移和台站到目标的长距离所共同决定；不同的噪声振幅对于轨道重叠的精度影响有所差异，但是基本处于相同的量级，相比较于其他观测手段在相同条件下，精度也基本符合。
图5 改进轨道与仿真轨道曲线
Fig.5 Overlap of POD and simulation orbit
4.2 实测数据处理
本文利用某次海上测控任务采集的数据进行分析，观测目标主要针对静止轨道卫星。选取观测弧段为2021年6月5日14∶00—15∶00以及16∶00—17∶30内的四颗GEO卫星进行数据解算，根据GEO卫星的NORAD 编号查询其对应的两行根数，转换成位置速度作为轨道改进的初始值，4颗GEO卫星在历元为2021年6月5日14时的初轨如表3所示(初轨使用无奇点根数表示)。
表3 4颗GEO目标的初始星历
Table 3 Initial ephemeris of four GEOs
IDa/kmh/10-5kξ/10-5η/10-5λ/°4104342 166.60118.402 0-0.002 08.332 7-8.213 1592.493 93780442 165.0947.647 7-0.001 29.395 6-1.970 3593.802 03723442 165.705103.030 4-0.031 7-4.687 2-32.766 4230.366 03276742 165.27340.238 1-0.002 0-2.230 4-10.669 3235.893 0
4.2.1 NORAD编号41043目标的定轨结果
NORAD编号为41043号的卫星名是LAOSAT-1，为老挝的通信卫星。测量设备从2021年6月5日14点09分16秒开始至2021年6月5日14点50分49秒，再由2021年6月5日16点29分12秒至2021年6月5日17点03分01秒对其进行观测，具体定轨的残差如图6所示。
从残差的结果可以看出，对41043目标的定轨残差基本在2 ″以内，这个船载的光学测量设备的观测噪声基本符合。
4.2.2 NORAD编号37804目标的定轨结果
NORAD编号为37804号的卫星名是CHINASAT-1A，为中国的烽火二号卫星。测量设备从2021年6月5日14点09分16秒开始至2021年6月5日14点50分56秒，再由2021年6月5日16点29分12秒至2021年6月5日17点27分53秒对其进行观测，具体定轨的残差如图7所示。
从残差的结果可以看出，对37804目标定轨的赤经残差基本在2 ″以内；赤纬残差大多数在5 ″以内，少数在5～10 ″，出现这种原因多数是由于惯导、形变测量在某些时刻未能完全修正船姿。但是整体来看依旧与船载的光学测量设备的观测噪声基本符合。
4.2.3 NORAD编号37234目标的定轨结果
NORAD编号为37234号的卫星名是CHINASAT-20A，为中国的神通一号卫星。测量设备从2021年6月5日14点09分16秒开始至2021年6月5日14点50分56秒，再由2021年6月5日16点17分42秒至2021年6月5日17点03分35秒对其进行观测，具体定轨的残差如图8所示。
图6 41043定轨残差示意图
Fig.6 POD residuals of 41043
图7 37804定轨残差示意图
Fig.7 POD residuals of 37804
图8 37234定轨残差示意图
Fig.8 POD residuals of 37234
从残差的结果可以看出，对37234目标定轨的赤经赤纬残差基本在5 ″以内，这个结果与船载的光学测量设备的观测噪声基本符合。
4.2.4 NORAD编号32767目标的定轨结果
NORAD编号为32767号的卫星名是VINASAT-1，为越南的通信卫星。测量设备从2021年6月5日14点09分16秒开始至2021年6月5日14点50分56秒，再由2021年6月5日16点29分12秒至2021年6月5日17点28分36秒对其进行观测，具体定轨的残差如图9所示。
图9 32767定轨残差示意图
Fig.9 POD residuals of 32767
从残差的结果可以看出，对32767目标定轨的赤经赤纬残差大多数在2 ″以内，少数在2～4 ″，这个结果与船载的光学测量设备的观测噪声基本符合。
四颗目标的实测定轨残差和仿真定轨残差基本吻合，仿真的观测弧段基本按照实测数据弧段进行筛选。对比仿真的定轨结果，海基单站实测目标的定轨精度在本文的观测条件下约为数百米，在观测情况较差的情况下，精度约在公里量级。
5 结论
利用海基单站对四颗非合作的GEO目标进行跟踪，使用光学测角数据，在观测弧段较短且数据不连续的情况下对其轨道进行改进。同时类比实测场景，进行仿真分析，对比实测结果和仿真结果可得：
1) 通过海基单站光学测量静止轨道卫星，精密定轨残差典型在角秒量级与测量设备标称值基本一致，反映了定轨策略的合理性。
2) 受海上测控资源限制，由于缺少高精度的已知目标标校比对，仿真相同场景下的静止轨道卫星数据处理，通过残差与轨道分析，海基单站测量对静止卫星定轨精度在百米至公里量级。实际任务中，由于受设备性能和测量弧长等因素影响，这一数据可能会有一定误差。
海基观测是航天测控网重要的一部分，通过海基测控对非合作目标的快速确认以及轨道改进，可为后续的其他测控手段甚至战略攻防提供参考。
[1] 于涌，罗浩，毛银盾，等.基于单站CCD漂移扫描光电设备的同步轨道目标测定轨[J].空间科学学报，2019,39(03):365-372.
Yu Y,Luo H,Mao Y D,et al.Orbit determination for geostationary earth orbit targets based on the CCD drift-scan photoelectric telescope at one single station[J].Chinese Journal of Space Science，2019，39(03):365-372.
[2] 宋叶志，黄勇，杨建华，等.风云四号卫星双程测距系统精密轨道确定[J].宇航学报,2020,41(03):270-275.
Song Y Z,Huang Y,Yang J H,et al.Orbit determination for FY-4 satellite using two-way ranging system[J].Journal of Astronautics,2020,41(03):270-275.
[3] Huang Y,Hu X G,Zhang X Z,et al.Improvement of orbit determination for geostationary satellites with VLBI tracking[J].Chinese Science Bulletin,2011,56(26):2765-2722.
[4] HUEBSCHMAN R K.The MSX spacecraft system design[J].Johns Hopkins APL Technical Digest,1996,17(01):41-48.
[5] STOKES G H,BRAUN C V,SRIDHARAN R,et al.The space-based visible program[J].Lincoln Laboratory Journal,1998,11(02):205-238.
[6] SHARMA J,STOKES G H,BRAUN C V,et al.Toward operational space-based space surveillance[J].Lincoln Laboratory Journal,2002,13(02):309-334.
[7] 宫经刚，宁宇，吕楠.美国高轨天基态势感知技术发展与启示[J].空间控制技术与应用，2021,47(01):1-7.
Gong J G,Ning Y,Lv N.Development and enlightenment of space based situational awareness technology for high orbit in the united states[J].Aerospace Control and Application,2021,47(01):1-7.
[8] 周海银，潘晓刚，李董辉.基于天基空间目标监视系统的定轨技术研究[J].系统仿真学报，2008,20(13):3538-3547.
Zhou H Y,Pan X G,Li D H.Research on orbit determination algorithm based on SBSS[J].Journal of System Simulation,2008,20(13):3538-3547.
[9] 周庆勇，杜兰，蓝朝桢.天基光学监视系统空间目标定轨精度仿真分析[J].系统仿真学报，2010,22(11):2660-2663.
Zhou Q Y,Du L,Lan Z Z.Analysis of orbit determination accuracy for space objects based on space-based optical surveillance system[J].Journal of System Simulation,2010,22(11):2660-2663.
[10] 康德勇，徐劲，茅永兴.海上移动站精密测轨方法[J].天文学报，2005(03):314-321.
Kang D Y,Xu J,Mao Y X.Precise orbit determination method for maritime floating observation station[J].Acta Astronomica Sinica,2005(03):314-321.
[11] 康德勇，李晓勇，王旭良，等.船位误差对外弹道测量及定轨精度的影响[J].电讯技术，2010,50(09):106-109.
Kang D Y,Li X Y,Wang X L,et al.Influence of ship’s position error on exterior trajectory measurement and orbit determination[J].Telecommunication Engineering,2010,50(09):106-109.
[12] 朱伟康，李辉芬，陈德明，等.提高航天测量船定轨精度的途径[J].电讯技术，2012,52(05):658-662.
Zhu W K,Li H F,Chen D M,et al.Ways to improve orbit determination precision of space TT&C ships[J].Telecommunication Engineering,2012,52(05):658-662.
[13] 陈红英，郭才发，向颉，等.基于最小二乘的海基双站测量数据综合处理[J].电讯技术，2015,55(12):1413-1416.
Chen H Y,Guo C F,Xiang J,et al.Data processing based on least squares for sea-based Bi-station measurement[J].Telecommunication Engineering,2015,55(12):1413-1416.
[14] 倪兴，李晓勇，杨磊，等.近似公式在航天测控中的应用[J].兵器装备工程学报，2016,37(06):174-176.
Ni X,Li X Y,Yang L,et al.Application of approximate formula on TT&C system[J].Journal of Ordnance Equipment Engineering,2016,37(06):174-176.
[15] 季勤超，赵建军，贺林波，等.舰船姿态测量误差对雷达探测精度的影响分析[J].计算机与数字工程，2019,47(02):330-333.
Ji Q C,Zhao J J,He L B,et al.Analysis of effects on radar detection precision caused by ship attitude measurement errors[J].Computer & Digital Engineering,2019,47(02):330-333.