基于精确星历最小二乘拟合的轨道快速递推方法.pdf

2023 年第 5 期计算机与数字工程收稿日期：2022年11月2日，修回日期：2022年12月12日作者简介：车征，男，硕士，工程师，研究方向：航天器轨道动力学。1引言高精度的轨道快速递推模型可为卫星敏感器、控制器等部件的在轨自主工作提供位置输入信息。由于星上计算资源和运算效率限制，无法使用高精度数值积分方式进行轨道预报。传统的摄动分析方法采用一阶或二阶摄动理论12，虽然能够进行较快速的平均根数与瞬时根数的转换和轨道递推，但需要计算大量的公式推导工作，且在摄动模型项较多时运算仍然较为复杂。GPS星历外推使用了一种参数拟合的方式，将摄动规律转为有限个参数描述的近似项35。参数个数从最初的16参数、发展出18参数、22参数等69。这种方法具有很高的拟合精度，但缺点是仅仅适用于短时间内的位置解算，当拟合时段变长则误差急剧增大，且没有包含速度信息，不适合用于长时间轨道递推1013。还有一种方法是多项式插值方式1415，在拟合区间内能够实现很高精度，但多项式阶数较高时计算量较大，且用于预报存在发散趋势。本文从瞬根数变化规律中提取轨道根数各分量的频率信息，使用含少量待定系数的三角函数和线性函数组合关系式描述轨道根数的变化趋势，使用精确动力学预报结果通过最小二乘拟合的方式确定系数值。避免了复杂的解析公式运算，同时保证了与精确动力学模型预报结果的吻合度。仿真算例表明拟合结果可以保证24h平均预报误差小于400m，且计算耗时小。基于精确星历最小二乘拟合的轨道快速递推方法车征（陕西星邑空间技术有限公司西安710043）摘要提出了一种基于第一类无奇点根数的低轨卫星瞬时根数拟合方法，根据精确星历预报结果规律，对每个根数分量设置不同的周期项和时间变化项，使用最小二乘法解算对应的系数。同时根据周期函数的特点设计了快速获得系数初值的解析计算式，减小了迭代计算量。方法可用于已知精确星历时的快速轨道递推计算；还可以快速分离出短周期和长周期以及线性时变项，实现瞬根数到平根数的转换。关键词星历预报；最小二乘；快速轨道递推中图分类号V19DOI：10.3969/j.issn.1672-9722.2023.05.016Fast Orbit Predict Method Based on the Least Square Fitting ofEphemeris DataCHE Zheng（Xingyi Space Technology Co.，Ltd.，Xian710043）AbstractAn ephemeris fitting method of LEO satellite based on the non-singular orbital elements is proposed.According tothe results of high precision orbit propagation，different periodic terms and linear terms are set to approximate the six elements.Then the least square method is applied to calculate the optimal coefficients to make the fitting function as close as the real values.The fast calculating method of initial values for the least square is also proposed to reduce the time cost.The method in this papercan be well used to quickly predict orbit elements.It can also separate the long terms out from the short terms and thus can be usedto quickly convert osculating elements into mean elements.Key Wordsephemeris propagation，least square method，quick orbit predictionClass NumberV19总第 403期2023 年第 5期计算机与数字工程Computer&Digital EngineeringVol.51No.51065第 51 卷2拟合模型基本思想为根据瞬根周期规律进行傅里叶分解，找出每个根数变化的短周期频率，分别指定拟合函数，对每个轨道根数单独拟合。以第一类无奇点根数为变量，定义各个分量的变化规律如下：1）a的周期为T/2，还包括随时间线性变化项。2）ex和ey的周期为T/3和T。3）i的周期为T/2和1/2天。4）的周期为T/2，还包括随时间线性变化项。5）的周期为T和1/2天，还包括随时间线性变化项。T为平均轨道周期。因此设置拟合函数如下式：a=a0+a1sin(tT4+a2)+a3tex=ex0+ex1sin(tT6+ex2)+ex3sin(tT2+ex4)ey=ey0+ey1sin(tT6+ey2)+ey3sin(tT2+ey4)i=i0+i1sin(tT4+i2)+i3sin(t864002+i4)=0+1sin(tT4+2)+t3=0+1sin(tT2+2)+t3+4sin(t432002+5)+6sin(tT4+7)（1）一共有30个参数。对六根数分别拟合，设待拟合的参数符号为 x（代表 a，ex，ey，i，中的一个），待拟合参数的初值为x0到xn，则：Dx=xx0Dx0+xxnDxn=xx0 xxnDx0.Dxn（2）利用不同时间共 m组 x的精确动力学模型预报数据，有L=xx0(t1)xxn(t1)xx0(tm)xxn(tm)Dx0Dxn=BDx0Dxn（3）其中L=xpre-x为按照式（1）和时间序列预报的xpre与真实动力学积分得到的x之间的偏差。则根据L可以反算出x0到xn应修正量和修正后的待拟合变量为Dx0Dxn=(BTB)-1BTLx0 xn=x0 xn+Dx0Dxn（4）其中（a，ex，ey，i，）对各自待拟合变量的偏导数为AA0=1AA1=sin(tT4+A2)AA2=A1cos(tT4+A2)exex0=1exex1=sin(tT6+ex2)exex2=ex1cos(tT6+ex2)exex3=sin(tT2+ex4)exex4=ex3cos(tT2+ex4)eyey0=1eyey1=sin(tT6+ey2)eyey2=ex1cos(tT6+ey2)eyey3=sin(tT2+ey4)eyey4=ey3cos(tT2+ey4)ii0=1ii1=sin(tT4+i2)ii2=i1cos(tT4+i2)ii3=sin(t864002+i4)ii4=i3cos(t864002+i4)0=11=sin(tT4+2)2=1cos(tT4+2)3=t0=11=sin(tT4+2)2=1cos(tT4+2)3=t4=sin(t432002+5)5=4cos(t432002+5)6=sin(tT4+7)7=6cos(tT4+7)（5）重复式（3）（6）直到计算的x0到xn应修正量接近0，即表示待拟合变量收敛。由于最小二乘方法需要一个较为准确的初值才能稳定收敛，下节给出初值的计算方法。3拟合初值计算方法3.1a的系数初值计算对于给定的数值积分获得的一系列时刻的 a（ti），i=1 到 N，首先根据大气衰减较小的假设，取a3=0。然后根据式（1），可知0tfa(ti)dt=a0tf-Ta12(cos(tfT2+a2)-cosa2)（6）当tf远大于T时，可以认为a0=0tfa(ti)dttf=i=1Na(ti)N（7）然后设轨道角速度n=a30+Dn，可得预测的卫星在最后时刻纬度幅角：(tN)=(t1)+(N-1)Dt(a30+Dn)（8）其中Dt为数据时间间隔，且(N-1)DtDn是一个小量，则车征：基于精确星历最小二乘拟合的轨道快速递推方法10662023 年第 5 期计算机与数字工程Dn=Mod(tN)-(t1)2)(N-1)DtT=2 a30+Dn（9）T作为常值用于随后各个轨道根数系数的拟合过程中。根据式（1）的性质，可知：a(0)=a0+a1sina2a(T8)=a0+a1sin(2+a2)=a0+a1cosa2（10）则：a1=(a(0)-a0)2+(a(T8)-a0)2a2=arctan2(a(0)-a0a(T8)-a0)（11）其中a(T8)对应a（ti）中的i=int(T8Dt)，其中Dt为数据时间间隔，后面类似符号（如T4T2等）同。3.2exey的系数初值计算类似式（9）可知ex0=i=1Nex(ti)N，根据式（1），有ex(0)=ex0+ex1sin(ex2)+ex3sin(ex4)ex(T3)=ex0+ex1sin(ex2)+ex3sin(23)cos(ex4)+cos(23)sin(ex4)ex(2T3)=ex0+ex1sin(ex2)+ex3sin(43)cos(ex4)+cos(43)sin(ex4)（12）则ex(T3)-ex(0)=ex3sin(23)cos(ex4)+(cos(23)-1)sin(ex4)ex(2T3)-ex(0)=ex3sin(43)cos(ex4)+(cos(43)-1)sin(ex4)（13）求解方程组得：ex3cos(ex4)=(ex(T3)-ex(0)-ex(2T3)+ex(0)3ex3sin(ex4)=(ex(T3)-ex(0)+ex(2T3)-ex(0)-3ex3=(ex3cos(ex4)2+(ex3sin(ex4)2ex4=arctan2(ex3sin(ex4)ex3ex3cos(ex4)ex3)（14）同理利用ex(T12)=ex0+ex1cosex2+ex3sin(6+ex4)（15）可计算：ex1cosex2=ex(T12)-ex0-ex3sin(6+ex4)ex1sinex2=ex(0)-ex0-ex3sin(ex4)ex1=(ex1cos(ex2)2+(ex1sin(ex2)2ex2=arctan2(ex1sin(ex2)ex1ex1cos(ex2)ex1)（16）ey的初值计算与式（15）到式（18）完全相同。3.3i的系数初值计算同 样 地，i0初 值 直 接 取 i 的 平 均 值，i0=i=1Ni(ti)N。根据i(T2)和i(T)的特性，有：i3sini4=-(i(T2)-i(0)sin(2T43200)+(i(T)-i(0)sin(T43200)(-cos(T43200)-1)sin(2T43200)+(cos(2T43200)-1)sin(T43200)i3cosi4=i(T2)-i(0)-i3sini4(cos(T43200)-1)sin(T43200)i3=(i3cosi4)2+(i3sini4)2i4=arctan2(i3sini4i3i3cosi4i3)（17）然后根据i(T8)的特性，有：i1sini2=i(0)-i0-i3sini4i1cosi2=i(T8)-i(0)-i3sin(i4+14T43200)i1=(i1cosi2)2+(i1sini2)2i2=arctan2(i1sini2i1i1cosi2i1)（18）3.4的系数初值计算根据：(0)=0+1sin2(T2)=0+1sin2+T23(T4)=0-1sin2+T43(T8)=0+1cos2+T83（19）可得：3=(T2)-(0)T 20=(T2)+(T4)-3T4321sin2=(0)-01cos2=(T8)-0-T831=(1cos2)2+(1sin2)22=arctan2(1sin211cos21)（20）3.5的系数初值首先3初值取为式（10）中的轨道角速度 n。1067第 51 卷根据下式：(0)=0+1sin(2)+4sin(5)+6sin(7)(T)=0+1sin(2)+T3+4sin(T432002+5)+6sin(7)(2T)=0+1sin(2)+2T3+4sin(2T432002+5)+6sin(7)（21）得：4sin5=(T)-(0)sin(2T432002)-(2T)-(0)sin(T432002)(cos(T432002)-1)sin(2T432002)-(cos(2T432002)-1)sin(T432002)4cos5=(T)-(0)-4sin5(cos(T432002)-1)sin(T432002)4=(4sin5)2+(4cos5)25=arctan2(4sin544cos54)（22）再根据下式：(T2)=0-1sin(2)+T23+4sin(T43200+5)+6sin(7)(T4)=0+1cos(2)+T43+4sin(T432002+5)-6sin(7)(3T4)=0-1cos(2)+3T43+4sin(3T4320024+5)-6sin(7)（23）得：1sin(2)=(0)-(T2)-4sin(5+T43200)-T231cos(2)=(3T4)-(T4)+T23-4sin(3T4320024+5)+4sin(T432002+5)20=(0)+(T4)-T43-1sin(2)-1cos(2)-4sin(5)-4sin(5+T432002)26sin(7)=(0)-0-1sin(2)-4sin(5)6cos(7)=(T8)-0-1sin(2+4)-4sin(5+T432004)6=(6sin7)2+(6cos7)27=arctan2(6sin766cos76)1=(1sin2)2+(1cos2)22=arctan2(1sin211cos21)（24）4测试结果初始瞬根为 7130km，0，98，0，0，360。则拟合一天的参数为a=7120.9688+9.0790sin(2t/T+1.5751)-3.3803e-08tex=-0.00046978+0.0007439sin(6t/T+1.5752)+0.0002941sin(2t/T-1.5749)ey=2.2291e-5+0.0007435sin(6t/T+0.004545)+0.0009327sin(2t/T-3.1340)i=1.7105+8.9871e-05sin(4t/T-1.5839)+2.7833e-05sin(2t/43200-0.6091)=2.0126e-05+9.0666e-05sin(4t/T-3.1393)+1.9014e-07t=0.0001174+0.0009410sin(2t/T-3.0928)+0.001049t+9.7748e-05sin(2t/43200-2.1751)+9.4583e-05sin(4t/T+0.008591)对应的残差换算到三维位置偏差见图1，平均总误差小于400m。设去除3倍周期以内的短周期项后的剩余长期项（常数项、长周期项和时间线性项）为平根数，与瞬根的对比如图2图7。图1位置残差图2半长轴长期项提取图3ex长期项提取车征：基于精确星历最小二乘拟合的轨道快速递推方法10682023 年第 5 期计算机与数字工程图4ey长期项提取图5i长期项提取图图6长期项提取图7长期项提取可见拟合获得的平根能够很好地反映轨道根数长期变化。5结语本文使用三角函数和时间线性函数来逼近低轨近圆轨道的轨道根数变化趋势，给出了最小二乘拟合函数系数的算法和流程，并设计了快速获得拟合系数初值的计算方法。仿真实例表明，拟合方法能够快速获得较高精度的轨道根数递推函数，24h数据拟合精度达到400m，且计算量很小。算法可用于星上快速轨道递推和瞬平根数转换。参 考 文 献1刘林，王彦荣.卫星轨道预报的一种分析方法 J.天文学报，2005，46（03）：307-313.LIU Lin，WANG Yanrong.An Analytical Method for Satellite Orbit ForecastingJ.ACTA Astronomica Sinica，2005，46（03）：307-313.2杨志涛，刘静，刘林.轨道分析解的改进方法及其应用J.系统工程与电子技术，2020，42（02）：427-433.YANG Zhitao，LIU Jing，LIU Lin.Improved method of orbit analytical solution and its application J.Systems Engineering and Electronics，2020，42（02）：427-433.3马开锋，彭碧波，洪樱.基于卫星轨道特征的低轨卫星星历参数拟合法 J.大地测量与地球动力学，2007（01）：85-90.MA Kaifeng，PENG Bibo，HONG Ying.LEO Satelliteephemeris parameters fitting method based on orbit character of satellite J.Journal of Geodesy and Geodynamics，2007，27（01）：85-90.4张如伟，刘根友.低轨卫星轨道拟合及预报方法研究J.大地测量与地球动力学，2008（04）：115-120.ZHANG Ruwei，LIU Genyou.Discussion on orbit fittingand orbit forecasting of low earth satellites J.Journal ofGeodesy and Geodynamics，2008，28（04）：115-120.5柳青松，刘峰，巨涛.小偏心率低轨卫星的星历参数拟合法 J.航天器工程，2010，19（04）：17-22.LIU Qingsong，BAI Muyi，JU Tao.Ephemeris ParametersFitting Method for Low Eccentricity LEO SatelliteJ.Spacecraft Engineering，2010，19（04）：17-22.6温渊，白沐炎，易灵，等.低轨卫星高精度轨道递推算法研究 J.上海航天，2019，36（2）：14-17.WEN Yuan，BAI Muyan，YI Ling，et al.LEO High Accuracy Orbit Extrapolation Algorithm J.Aerospace Shanghai，2019，36（2）：14-17.7方善传，杜兰，高云鹏，等.LEO卫星轨道根数型星历参数与接口设计 J.测绘学报，2019，48（02）：198-206.FANG Shanchuan，DU Lan，GAO Yunpeng，et al.Orbital elements ephemerides and interfaces design of LEO satellitesJ.Journal of Geodesy and Geoinformation Science，2019，48（02）：198-206.8关梅倩，焦文海，贾小林，等.基于轨道根数的LEO广播星历参数设计 J.大地测量与地球动力学，2018，38（09）：948-952.GUAN Meiqian，JIAO Wenhai，JIA Xiaolin，et al.Broadcast Ephemeris Design of LEO Satellite Based on the Orbital Elements J.Journal of Geodesy and Geodynamics，2018，38（09）：948-952.（下转第1095页）10692023 年第 5 期计算机与数字工程J.Review of Educational Research，2010（80）：272-295.9Grave B.The effect of student time allocation on academicachievementJ.Education Economics，2011（19）：291-310.10Wald A，Muennig P A，Oconnell K A，et al.Associations between healthy lifestyle behaviors and academicperformance in US undergraduates：a secondary analysisof the American College Health Associations NationalCollege Health Assessment IIJ.American Journal ofHealth Promotion，2013（28）：298-305.11Taylor D J，Vatthauer K E，Bramoweth A D，et al.Therole of sleep in predicting college academic performance：is it a unique predictor？J.Behavioral SleepMedicine，2013（11）：159-172.12Brinton C G，Chiang M.MOOC performance predictionvia clickstream data and social learning networks A.InProceedings of the 2015 IEEE Conference on ComputerCommunicationsC.New York（NY）：IEEE Press，2015：2299-2307.13Wang R，Harari G，Hao P，et al.SmartGPA：how smartphones can assess and predict academic performance ofcollege studentsA.In Proceedings of the 2015 ACMInternational Joint Conference on Pervasive and Ubiquitous ComputingC.New York（NY）：ACM Press，2016：295-306.14Cao Y，Gao J，Lian D，et al.Orderliness predicts academic performance：behavioural analysis on campus lifestyle J.Journal of The Royal Society Interface，2018：1-8.15Xu P，Yin L，Yue Z，et al.On predictability of time series EB/OL.https：/arxiv.org/abs/1806.03876.16Spearman C.The proof and measurement of associationbetween two things J.International Journal of Epidemiology，2010（39）：1137-1150.17王鼎，任虹宇，王倩.综合素质评价结果有效性研究基于学业成绩与综合素质评价指标的相关性分析 J.教育发展研究，2019，39（20）：18-28.WANG Ding，REN Hongyu，WANG Qian.Research onEffectiveness of Comprehensive Quality Evaluation Results：Correlation Analysis Based on Academic Performance and Comprehensive Quality Evaluation Index J.Research in Educational Development，2019，39（20）：18-28.9黄华，何峰，刘林.广播星历参数物理意义分析与相关性研究 J.宇航学报，2014，35（02）：171-176.HUANG Hua，HE Feng，LIU Lin.Research on PhysicalMeaning and Correlation of Broadcast Ephemeris Parameters J.Journal of Astronautics，2014，35（02）：171-176.10李兆铭，杨文革，丁丹，等.带摄动力拟合的低轨卫星实时定轨 STCKF 算法 J.电波科学学报，2016，31（05）：843-850.LI Zhaoming，YANG Wenge，DING Dan，et al.An STCKF algorithm for LEO satellite orbit determination withdisturbing fitting function J.Chinese Journal of RadioScience，2016，31（05）：843-850.11毛悦，宋小勇，胡小工，等.自主导航星上快速星历拟合算法研究 J.测绘科学，2017，42（08）：1-7.MAO Yue，SONG Xiaoyong，HU Xiaogong，et al.Research on on-board broadcast ephemeris fast fitting algorithm used on autonomous navigation satelliteJ.Science of Surveying and Mapping，2017，42（08）：1-7.12冯来平，崔先强，贾小林，等.GPS 16参数广播星历改进参数的变化规律分析 J.空间科学学报，2017，37（02）：214-221.FENG Laiping，CUI Xianqiang，JIA Xiaolin，et al.Analysis of the Changing Regularity of Improved Parameters of 16-parameter GPS Broadcast Ephemeris J.Chinese Journal of Space Science，2017，37（02）：214-221.13刘金健.基于拉格朗日插值法的GPS卫星轨道位置拟合 J.科技创新与生产力，2018（07）：11-13.LIU Jinjian.Fitting of GPS Satellite Orbit Location Basedon Lagrange Interpolation MethodJ.Taiyuan Scienceand Technology，2018（07）：11-13.14李振昌，李仲勤.滑动式切比雪夫多项式拟合法在BDS精密星历内插中的应用 J.测绘工程，2019，28（04）：49-53.LI Zhenchang，LI Zhongqin.Application of sliding Chebyshev polynomial fitting method to Beidou satellite precise ephemeris interpolation J.Engineering of Surveying and Mapping，2019，28（04）：49-53.15袁晓波，张超，杜兰.基于第一类无奇点根数的LEO历书参数设计 J.测绘科学技术学报，2016，33（06）：572-576.YUAN Xiaobo，ZHANG Chao，DU Lan.Parameter Design of LEO Almanac Ephemeris Based on the FirstClass of Non-Singularity Variables J.Journal of Geomatics Science and Technology，2016，33（06）：572-576.（上接第1069页）1095