0 引言
雷达的优点是白天黑夜均能探测远距离的目标，且不受雾、云和雨的阻挡，具有全天候、全天时的特点，并有一定的穿透能力。雷达已广泛应用于社会经济发展(如气象预报、资源探测、环境监测等[1-2])和科学研究(天体研究、大气物理、电离层结构研究等[3])中，但其应用最广的还是军事领域，如防空警戒、反导预警、空间监视、指挥引导等[4-5]。
空间目标监视雷达是空间监视领域的骨干装备，与常规防空雷达相比，主要增加了空间目标轨道确定功能，轨道改进实现该功能的核心模块，轨道改进模块的研制涉及非常复杂的航天动力学知识，这对从事雷达系统设计的非轨道力学专业的研究人员来说带来极大困难。虽有一些开源定轨软件可供借鉴[6-10]，但这些软件主要用于处理光学观测数据，并不适合处理雷达观测数据。为此，作者着手开发RadarOrbDet函数库[11]来辅助非轨道力学专业的人员进行空间目标监视雷达的系统设计，目前已完成坐标转换模块[12]、初轨确定模块[13]和摄动分析模块[14-16]的开发，本文研究轨道改进模块的开发。
轨道改进是空间目标轨道确定的重点和难点，它涉及的内容较多，包括摄动分析、偏微分方程的变分法分析、最小二乘优化等理论[17]，本文聚焦空间目标雷达定轨中的批处理轨道改进问题，着重介绍了最小二乘轨道改进原理，讲解了模型中的偏导数求解方法，描述了程序实现流程，最后用ODTK软件验证了本文软件设计的有效性。
1 最小二乘轨道改进原理
最小二乘轨道确定的基本原理是找到一条轨道，使得理论观测值和实际观测值之间的残差平方和(即损耗函数)最小，即求解x0，使得式(1)的值最小。
J(x0)=[z-h(x0)]T[z-h(x0)]
(1)
式中:x0=[r(t0);v(t0)]是卫星在t0时刻的位置和速度构成的状态参量，由x0通过轨道预报可以得到一条参考轨道；z=[z1;z2;…;zN]是N次雷达实际观测矢量，每次观测矢量zn包含距离、方位和俯仰三个分量(ρn,αn,βn)；h(x0)是由参考轨道计算得到的理论观测值。
由于h是一个非线性函数，式(1)描述的是一个非线性最小二乘问题，求解异常困难，通常采用泰勒展开将其转化为一个线性最小二乘问题，然后通过迭代方式求解它的线性最小二乘近似解，表达式为
(2)
式中：
是一个由观测误差构成的对角方阵，起加权作用。迭代初值为
它是x0的一个初始估计，通常可由初轨确定得到[17]。迭代过程直到连续两次损耗函数的相对改变量小于某个给定的门限时终止。雅克比矩阵H为
(3)
2 模型中的偏导数求解
式(2)的求解关键在于偏导数雅克比矩阵H的计算，不失一般性，以某一时刻的雷达观测量为例。为书写方便，除t0时刻外，省略其他时刻索引。
利用求导链式法则，H转换为
(4)
式中，矩阵A是当前时刻观测量对当前状态参量的偏导数矩阵，Φ是当前时刻状态参量对初始时刻状态参量的偏导数矩阵，也称为状态误差传递矩阵。
2.1 观测偏导数矩阵求解
在忽略光行差等因素后，雷达观测量只与卫星当前的位置有关，而与其速度无关。因此，矩阵A可拆分成
(5)
式中，rECI是卫星在惯性坐标系(ECI坐标系)下的位置矢量。
雷达观测是基于站心地平坐标进行的(如南-东-天坐标系，SEZ)，直接计算观测量对ECI坐标系下的位置矢量的偏导数比较困难，因此，A矩阵的计算需进一步运用链式求导法则进行分解：
(6)
式中，ρSEZ是卫星在SEZ坐标系下的直角坐标，用(ρS,ρE,ρZ)表示；ρECEF是雷达对卫星的观测矢量在地固坐标系(ECEF坐标系)下的表示；rECEF是卫星在ECEF坐标系的位置矢量。
注意到
ρECEF= rECEF- rsiteECEF
(7)
式中rsiteECEF是雷达在ECEF坐标系下的位置矢量，是个常矢量。
因此
(8)
根据坐标变换有
(9)
式中
为SEZ坐标系到ECEF坐标系的变换矩阵。
根据式(9)有
(10)
同理有
(11)
将式(8)、(10)和(11)代入式(6)有
(12)
式中
为SEZ坐标系到ECI坐标系的变换矩阵。
在SEZ坐标系下，由直角坐标和极坐标的关系，易得
(13)
将式(13)代入式(12)，有
(14)
计算式(14)中的两个矩阵时，都需要在当前观测时刻进行。最后，将式(14)代入式(5)后，可得矩阵A。
2.2 状态误差传递矩阵求解
设卫星的运动方程为
(15)
则
(16)
其中，矩阵F为
(17)
由于采用迭代方式进行最小二乘问题求解，因此对Φ矩阵的求解精度要求不高，从而对F矩阵的求解精度也要求不高。因此，F矩阵中的加速度可以只考虑地球中心引力和低阶的非球形引力摄动，以加快计算机求解速度[18]。
在只考虑地球引力(包含非球形摄动)情况下，加速度只与位置有关，而与速度无关，并且它的计算通常是在ECEF坐标系下进行，具体表达式可以参考文献[18]，篇幅所限，这里不再赘述。
最后还需将ECEF坐标系下的偏导数矩阵转换到ECI坐标系下。在忽略科里奥利力和离心力情况下，有如下关系式：
(18)
计算得到当前时刻的F矩阵后，采用数值微分方程方法求解Φ。
3 程序设计
图1是最小二乘轨道改进算法的程序设计流程图。
图1 最小二乘轨道改进算法的程序设计流程图
图中每点雷达观测数据
包括距离、方位和俯仰(ρ,α,β)三个分量；利用数值法求解参考轨道时，摄动因素要结合问题求解的精度尽可能考虑全面；而求解状态误差传递矩阵时可以只考虑地球中心引力和低阶非球形摄动。
程序分为三层循环，第一层循环是数值微分方程求解的内部循环，输出结果是第i时刻的参考轨道和状态误差传递矩阵；第二层循环是遍历每个观测值，根据最小二乘原理求解单次轨道改进量；第三层循环是执行多次最小二乘轨道改进，使最终结果趋近于实际轨道状态。
4 仿真实验
4.1 接口调用示例
采用C语言按照第3部分的程序设计原理实现了最小二乘轨道改进的接口函数rodLeastSquare，并集成到RadarOrbDet库中。rodAccelHarmonic接口函数的软件调用方法如图2所示。
图2 最小二乘轨道改进接口函数说明
4.2 实验比较
在STK11.2中仿真一颗卫星，其历元轨道参数设置如表1所示，运动模式设置为HPOP(高精度轨道预报)，考虑21阶非球形摄动。雷达站址的经纬高为(120°,30°,0)。选择2021-09-07 19:15:01至2021-09-07 19:15:10间隔1秒的10点数据，对距离、方位和俯仰分别加上100 m、0.1°和0.1°的观测误差。
表1 卫星轨道根数
历元时刻半长轴偏心率轨道倾角近地点幅角升交点赤经真近点角2021-09-07 04:00:006678.14km028.5‱0‱0‱0‱
采用本文软件和AGI公司开发的一款先进的轨道确定与分析软件[19]——Orbit Determination Tool Kit(ODTK)进行定轨比较，ODTK与本文软件计算结果如表2所示，采用式(19)所描述的位置和速度相对误差进行定量比较。
(19)
从表2可以看出，定轨的位置误差为0.059 6%，速度误差为0.35%，表明本文软件算法是有效的，其误差主要来源为计算程序的坐标与时间转换、计算截断误差等。
表2 定轨结果对比表(ODTK与本文软件)
类别定轨结果(历元:19:15:00;单位:km, km/s)xyzvxvyvzODTK5189.103690.492003.17软件5191.773692.332005.47相对误差0.0596%-4.84075.28982.8745-4.86145.27872.86100.35%
5 结束语
本文针对空间目标雷达定轨中的批处理轨道改进问题，介绍了最小二乘轨道改进的原理，描述了模型中的偏导数求解，给出了程序设计思路，仿真表明了软件设计的有效性。此外，本文开发的RadarOrbDet函数库可以很好地辅助空间目标监视雷达系统设计人员进行定轨指标的快速设计和验证。
[1] SOKOL Z. The Role of Weather Radar in Rainfall Estimation and Its Application in Meteorological and Hydrological Modelling—A Review[J]. Remote Sensing, 2021, 13(3):351.
[2] GUTIERREZ R E, KUMJIAN M R. Environmental and Radar Characteristics of Gargantuan Hail-Producing Storms[J]. Monthly Weather Review, 2021(8):2523-2538.
[3] DING Chunyu, XIAO Zhiyong, WU Bo, et al. Rock Fragments in Shallow Lunar Regolith: Constraints by the Lunar Penetrating Radar Onboard the Chang’E-4 Mission[J].Journal of Geophysical Research:Planets, 2021,126(9):1-8.
[4] LUO Ru, CHEN Lifu, XING Jin, et al. A Fast Aircraft Detection Method for SAR Images Based on Efficient Bidirectional Path Aggregated Attention Network[J]. Remote Sensing, 2021, 13(15):2940.
[5] CATALDO D, GENTILE L, GHIO S, et al. Multibistatic Radar for Space Surveillance and Tracking[J]. IEEE Aerospace and Electronic Systems Magazine, 2020, 35(8):14-30.
[6] University of Pisa. The OrbFit Software Package[EB/OL].(2017-04-20)[2020-10-20]. http:∥adams.dm.unipi.it/～orbmaint/orbfit/.
[7] Pasquale Tricarico. Orbit Reconstruction, Simulation and Analysis ORSA[EB/OL]. (2018-09-14) [2020-10-20]. http:∥orsa.sourceforge.net/index.html.
[8] European Space Agency. Asteroid and Comet Trajectory Propagator[EB/OL]. (2015-11-23) [2020-10-20]. http:∥neo.ssa.esa.int/neo-propagator.
[9] The United States Naval Observatory. NOVAS[EB/OL]. (2016-11-02)[2020-10-20]. https:∥pypi.org/project/novas/.
[10] National Aeronautics and Space Administration. Orbit Determination Toolbox[EB/OL]. (2019-05-15) [2020-10-20]. https:∥sourceforge.net/projects/odtbx/.
[11] 黄晓斌. 雷达定轨函数库(RadarOrbDet) [EB/OL]. (2021-09-01)[2021-09-01]. https:∥pan.baidu.com/s/19kYFKZy77awg6pNnL-Y73Q.
[12] 黄晓斌, 张燕, 肖锐. 空间目标的雷达定轨坐标转换问题[J]. 导航定位学报,2020,8(4):39-43.
[13] 张燕, 黄晓斌, 肖锐, 等. 空间目标的雷达初轨确定问题探讨[C]∥2020年全国航空电子信息技术高端论坛论文集(上),北京: 中国电子科技集团有限公司, 2020:598-600.
[14] 黄晓斌, 张燕, 肖锐, 等. 空间目标的地球引力加速度分析及工程实现[J]. 测绘与空间地理信息,2021,44(9):200-204.
[15] HUANG Xiaobin, ZHANG Yan, XIAO Rui, et al. Analysis of Sun-Moon Gravitational Acceleration and Engineering Realization of Space Target[C]∥ 2021 International Conference on Computer, Communication, Control, Automation and Robotics,Shanghai:IOP Publishing Ltd, 2021:012036.
[16] HUANG Xiaobin, ZHANG Yan, XIAO Rui. Analysis of Atmospheric Drag Acceleration and Engineering Realization of Space Target[C]∥ 2021 7th International Conference on Mechanical Engineering and Automation Science (ICMEAS 2021), Seoul： IEEE, 2021：199-203.
[17] VALLADO D A. Fundamentals of Astrodynamics and Applications[M]. 4th ed. New York：McGraw-Hill, 2013.
[18] MONTENBRUCK O, GILL E. 卫星轨道——模型、方法和应用[M]. 王家松, 祝开建, 胡小工,译. 北京:国防工业出版社, 2012.
[19] VALLADO D A, HUJSAK R S, JOHNSON T M, et al. Orbit Determination Using ODTK Version 6[C]∥Madrid, Spain:ESA/ESAC Astronomy Centre,2010.https:∥www.researchgate.net/publication/228532272.
