力学计算与仿真（测绘学报阮仁桂）
力学计算与仿真（测绘学报阮仁桂）阮仁桂 贾小林 冯来平 等. 北斗三号MEO卫星非保守力建模[J]. 测绘学报，2022，51(9)：1862-1869. DOI: 10.11947/j.AGCS.2022.20210319关键词：北斗三号全球卫星导航系统 太阳光压 地球辐射 天线辐射 非保守力 激光测距残差 1. 西安测绘研究所 陕西 西安 710054;2. 地理信息工程国家重点实验室 陕西 西安 710054基金项目：国家自然科学基金(42074025;41774012);地理信息工程国家重点实验室基金(SKLGIE2019-ZZ-3)
本文内容来源于《测绘学报》2022年第9期（审图号GS京（2022）0824号）
北斗三号MEO卫星非保守力建模
阮仁桂1 2
贾小林1 2 冯来平1 2 王龙1 2 张奋1 2
1. 西安测绘研究所 陕西 西安 710054;
2. 地理信息工程国家重点实验室 陕西 西安 710054
基金项目：国家自然科学基金(42074025;41774012);地理信息工程国家重点实验室基金(SKLGIE2019-ZZ-3)
摘要：非保守力模型精度不高是制约BDS-3卫星定轨精度的主要因素之一。本文针对BDS-3 MEO卫星构建了地球辐射、天线辐射和箱体-两翼(BW)太阳光压模型 对典型的经验光压模型(ECOM1和ECOM2)进行补偿得到多个非保守力模型 收集全球观测网的数据进行定轨试验 通过轨道重叠互差和激光测距残差分析比较不同轨道模型的优劣。试验结果表明 经验光压模型是影响轨道精度的主要因素 在名义偏航模式下 ECOM2具有更好的表现 但ECOM1对卫星的姿态模式更不敏感。地球辐射和天线辐射会引起北斗卫星轨道径向约3 cm的系统性偏差 对二者建模后 几乎可以完全消除卫星C29和C30的激光残差系统偏差 但卫星C20和C21的系统偏差反而增大。此外 增加box-wing模型对于提高轨道精度也是有益的。关键词：北斗三号全球卫星导航系统 太阳光压 地球辐射 天线辐射 非保守力 激光测距残差
阮仁桂 贾小林 冯来平 等. 北斗三号MEO卫星非保守力建模[J]. 测绘学报，2022，51(9)：1862-1869. DOI: 10.11947/j.AGCS.2022.20210319
RUAN Rengui JIA Xiaolin FENG Laiping et al. Modeling non-conservative forces for BDS-3 MEO satellites[J]. Acta Geodaetica et Cartographica Sinica 2022 51(9): 1862-1869. DOI: 10.11947/j.AGCS.2022.20210319
阅读全文：http://xb.chinasmp.com/article/2022/1001-1595/20220903.htm
引 言北斗三号全球导航卫星系统(BDS-3)已于2020年7月31日正式开通运行。得益于独有的全星座星间链路，BDS-3的空间信号精度和定位服务性能已跻身主要GNSS的前列，其中广播轨道的径向精度优于10 cm[1-2]。可以预见，BDS-3将在高精度定位、时间频率传递、参考框架建立与维持等领域发挥越来越大的作用[3]。
高精度的卫星轨道是实现BDS-3高精度应用的前提。精准确定卫星轨道，需要精细地考虑卫星受到的各种摄动力[4]。BDS-3卫星运行于中、高轨道，力模型的不确定性主要来自太阳光压、地球辐射等非保守力摄动。为提高定轨精度，文献[5—7]对BDS-3卫星的轨道模型进行了探讨，主要工作集中于改进BDS-3卫星的太阳光压模型。针对BDS-3卫星平台为长方体的特点，文献[5]基于可估计参数的箱体-两翼(box-wing BW)模型[8-9]和频谱分析方法构建了两个先验补偿模型以增强五参数ECOM(extend CODE orbit model)[10]模型(ECOM1)。文献[6]同样基于可估参数的BW模型，针对不同厂家的BDS-3卫星建立了分析型光压模型。为了支持高精度应用[11]，中国卫星导航系统管理办公室发布了北斗卫星的元数据[12]，包括几何和物理参数。文献[7]根据BDS-3卫星元数据建立的BW模型联合ECOM1可以获得优于九参数的扩展ECOM模型(D4B1)[13]的定轨结果。与太阳光压相比，地球辐射和天线辐射产生的摄动力要小得多，研究表明二者会引起导航卫星轨道径向的系统性偏差[14-16]，但它们对BDS-3定轨的影响尚未被涉及。
文献[6—7 17]研究表明，即使考虑了BW模型补偿，BDS-3卫星轨道的SLR残差仍然有不可忽视的系统性偏差，这说明BDS-3卫星的非保守力模型还有待进一步完善。本文基于BDS-3卫星元数据和天线辐射功率信息构建地球辐射、天线辐射和BW光压模型对典型的经验型太阳光压模型进行补偿，考查各非保守力因素或模型对精密定轨的影响，以期优选出适合BDS-3 MEO卫星的非保守力建模策略。
1 BDS-3卫星的非保守力建模BDS-3卫星受到的非保守力包括[4]：太阳光压、地球辐射、天线辐射和卫星体热辐射。其中，卫星体热辐射与卫星体表面辐射率、绝对温度有关，难以准确建模且部分被光压模型吸收[13]，在这里不单独考虑。
1.1 太阳光压摄动
太阳光压是BDS-3卫星承受的最大的非保守力。针对导航卫星的太阳光压建模，国内外学者已经积累了非常丰富的经验[10 13 18-21]，大体上可以分为分析法和经验法。两种方法各有优劣，为了获得最佳效果，联合采用经验法和分析法对太阳光压进行建模逐渐成为一种趋势。因此，太阳光压引起的摄动加速度aSRP可表示为
式中，aapr表示先验的光压摄动加速度，通常采用分析模型进行计算，如ROCK模型[18]、BW模型[9]等，其中的模型参数需事先测定或计算得到；aemp表示经验模型，如GSPM[20]、ECOM系列[10]、可估计参数的BW模型[8]。本文将采用BW模型和ECOM系列模型进行太阳光压建模。
1.1.1 BW模型
BW模型的核心思想是将卫星视为一个箱体和两翼的太阳帆板所构成，在给定各表面的几何、物理参数及卫星姿态信息的条件下，可以计算出太阳光压在各个表面产生的摄动加速度，将各个面的加速度叠加就得到太阳光压引起的总加速度。由于卫星体和太阳帆板表面材料的热辐射特性不同，分别用式(1)和式(2)计算卫星体各表面元(i)和太阳帆板(sp)的太阳光压加速度[9]
(1)
(2)
式中，a*(*=i sp，下同)表示加速度矢量；A*表示面积；e⊙表示太阳方向矢量；en *表示表面元(*)的法向矢量；θ*为e⊙和en *的夹角；α*、δ*和ρ*表示表面元(*)的吸收、漫反射和镜面反射系数；m表示卫星质量；c为光速；Φ⊙表示太阳在卫星所处位置的辐射通量，计算公式为
式中，Φ0≈1367 W/m2表示太阳在地球附近的长期平均辐射通量；1AU表示太阳到地球的平均距离；r⊙和r分别表示太阳和卫星的位置矢量。
1.1.2 ECOM模型
ECOM模型是欧洲定轨中心(Center of Orbit Determination atEurope CODE)针对GPS卫星提出的经验型太阳光压模型[10]，现已被广泛应用于各种导航卫星的精密定轨。其核心思想是将太阳光压加速度在3个正交方向eD、eY和eB上进行表示
(3)
式中，eD为卫星指向太阳的矢量；eY为“日-星-地”平面法向矢量；eB=eD×eY，构成右手系；D(u)、Y(u)和B(u)为对应3个方向的加速度，采用截断傅里叶级数表示，其具体计算公式有多个不同的版本[10]，最新的扩展(或改进的)ECOM模型为
(4)
式中，[D0 Y0 B0]、[D2i c B2i－1 c]、[D2i s B2i－1 s]上进行(i=1 2 …)是扩展ECOM模型参数，需作为待估参数进行解算；Δu=u－u⊙ u⊙和u分别表示太阳在轨道面的投影和卫星相对于轨道升交点的夹角。
式(4)中eD方向的加速度仅考虑了偶次项，eB方向仅考虑了奇次项。研究表明，对于卫星体为正方体的GPS卫星，仅考虑{D0 Y0 B0 Bc Bs}五参数的简化ECOM模型，即ECOM1模型，可获得最优的轨道和有关大地参数的解算结果[10]。对于非正方体的卫星，文献[13]推荐七参数{D0 D2 c D2 s Y0 B0 Bc Bs}的D2B1模型(下文用ECOM2表示)或九参数{D0 D2 c D2 s D4 c D4 s Y0 B0 Bc Bs}的D4B1模型[13 21]。
1.2 地球辐射摄动
地球辐射摄动是指地球表面受到太阳光照射后直接反射或红外辐射对卫星引起的辐射压力。文献[15]建立了一系列的地球辐射通量模型，包括分析型、数值型和数据型[15]。其中，将地球视为球形朗伯辐射体的分析模型计算简单，与基于CERES(clouds and Earth's radiant energy system)数据的结果相比，最大误差不超过15%，平均误差为1.6%。这也说明准确计算地球辐射通量非常困难。分析模型计算直接反射Φrefl和红外辐射Φinfr通量的公式分别为
(5)
式中，RE为地球半径；AE=πRE2，为地球截面积；αalb为地球表面漫反射系数，通常取0.3；ε为星-地-日夹角；h为卫星的轨道高度。
用Φrefl和Φinfr代替Φ⊙，用地球方向矢量e⊕代替e⊙，就可用式(1)和式(2)计算地球辐射引起的卫星体各表面和太阳帆板的摄动加速度。
1.3 天线辐射摄动
天线辐射摄动是卫星的导航信号天线持续对地发射射频信号产生的与辐射方向相反的摄动力。天线辐射加速度aAT可由式(6)计算
(6)
式中，P为天线信号发射功率。
1.4 卫星信息
中国卫星导航系统管理办公室提供了BDS-3卫星表面几何和光学参数文件[12]。各类卫星的参数见表 1，其中用CAST和SECM区分中国空间技术研究院和上海微小卫星工程中心研制的卫星。由于文献[12]仅给出了各表面的吸收系数，没有提供镜面反射率。本文根据文献[7]，假设各卫星各个表面的镜面反射率为1，仅CAST卫星 X和-Z面的镜面反射率为0；同时假设太阳帆板反面的光学参数与正面相同，且所有表面材料对红外辐射的镜面反射、漫反射和吸收系数分别为0.1、0.1和0.8。GNSS卫星的天线辐射功率可以用大口径天线进行测量[14]，本文进行天线辐射建模时采用IGS(international GNSS service)提供的数据，SECM卫星和CAST卫星的天线辐射功率分别为280、310 W。
表 1 BDS-3卫星的几何和物理参数
Tab. 1 Geometric and physical parameters of BDS-3 satellites
姿态信息是建立分析型光压模型和地球辐射模型的另一关键信息。BDS-3 MEO和IGSO卫星在大部分时间里采用名义偏航(NYS)模式：即 Z轴指向地心， X轴位于“日-星-地”平面并指向太阳侧。当卫星进入深度地影区(太阳相对于轨道面的高度角小于3°)且日、星、地三者接近共线时进入控制偏航(CYS)模式： X轴与日-星-地平面存在一定夹角，具体的姿态模型见文献[7]。根据上文的公式和有关数据可得：CAST和SECM卫星的太阳光压加速度分别约为130、70 nm/s2，地球辐射加速度分别约为1.5、0.5 nm/s2，天线辐射加速度约为1.5 nm/s2。
2 试验和分析为了分析不同非保守力因素对BDS-3卫星精密定轨的影响，收集全球观测网的数据，构建多个非保守力模型进行定轨试验。
2.1 数据及试验设计
收集2021年年积日001—110期间全球约75个监测站的观测数据，测站分布如图 1所示。采用西安测绘研究所的SPODS软件[22-23]进行定轨解算，定轨弧长为30 h，以每天的12∶00为中心。按照IGS或iGMAS分析中心的日常数据处理模式，同时解算轨道、钟差、测站坐标、对流层延迟等参数，具体策略见表 2。
图 1 本文所用的IGS和iGMAS测站分布Fig. 1 Distribution of the selected IGS and iGMAS stations
图选项
表 2 数据处理算法、模型和策略Tab. 2 Algorithms models and strategies for data processing
本文首先关注经验型太阳光压模型ECOM1和ECOM2的优劣，在此基础上考查地球辐射、天线辐射和先验光压模型的影响。初步分析发现地球辐射与天线辐射主要作用于径向方向，对定轨结果的影响相似。为此，设计了14组不同的非保守力模型配置，见表 3。模型的命名采用以下规则：“B”代表BW模型; “Sn”(n=1、2)分别代表ECOM1和ECOM2；“E”代表地球辐射；“A”代表天线辐射；“*”表示“E”和“A”仅用于SECM卫星。在不至于混淆的情况下，下文用ECOM1和ECOM2分别代表以ECOM1和ECOM2为经验光压模型的一类非保守力模型。
表 3非保守力模型配置表
Tab. 3 Configuration of the non-conservative force models
2.2 结果分析
通过考查相邻两天定轨结果6 h重叠段的互差进行轨道精度的内部评估；分析激光数据残差进行外部检核。考虑到IGSO和GEO卫星的区域覆盖特性，几何观测条件较差，本文不对其结果进行分析。
2.2.1 轨道重叠互差分析
在试验数据弧段内，C20、C21所在的轨道和C29、C30所在轨道分别在年积日033—040和083—094期间进入深度地影区，轨道面内的卫星会进入CYS模式。图 2统计了NYS和CYS模式下不同模型得到的BDS-3 MEO卫星轨道重叠互差在R、T和N方向的RMS。由图 2可以看出，在NYS模式下，不同力模型同一轨道方向的精度差异不大，3个方向的RMS分别在3、5和5.5 cm左右，只是BS1和BS1EA*在轨道T和N方向的精度似乎较其他模型稍差。力模型相同时，与NYS模式的结果相比，CYS模式下轨道各方向的误差都明显增大，增幅普遍超过1 cm；其中ECOM2在R方向的增幅超过3 cm，这说明，ECOM2对姿态模式更敏感。还可看出，考虑地球辐射、天线辐射和BW模型不能显著降低轨道重叠互差。
图 2 不同模型得到的轨道重叠互差RMSFig. 2 RMS of overlap orbit differences with different force models
图选项
2.2.2 激光残差分析
BDS-3卫星都安装了激光反射棱镜，国际卫星激光测距组织对4颗BDS-3卫星(C20、C21、C29和C30)的激光测距数据可以用于评估卫星轨道径向精度。在进行激光残差计算时，各颗卫星的激光反射器的质心偏差参数见表 4[12]。值得注意的是，表 4中的数值与文献[5—7 17]所采用的数值不完全一致。
表 4 BDS-3卫星的激光反射器质心偏差参数
Tab. 4 Center of mass corrections of the laser array for BDS-3 satellites
图 3绘出有代表性的6个力模型得到的激光残差随星-地-日夹角ε的变化，红色数值为红色拟合直线的斜率。可以看出，力模型相同时，卫星C20和C21的激光残差具有更明显的随ε变化的趋势。在经验光压模型基础上，加入地球辐射和天线辐射模型不能减弱激光残差与ε的相关性；而加入BW模型有助于降低CAST卫星激光残差与ε的相关性。对比图 3(a)、(d)可见，相比于ECOM1，采用ECOM2可明显减弱激光残差与ε的相关性。
图 3 激光残差随ε的变化Fig. 3 Variation of SLR residuals withε
图选项
同一厂家卫星的激光残差具有相似的统计特性，图 4统计了两类卫星激光残差的STD。由图 4可以看出，在NYS模式下，不同模型的激光残差STD基本都小于4.5 cm；ECOM2对应的STD更小，基本不超过3.5 cm。在CYS模式下，ECOM1对应的STD普遍略有增大，但都不超过5 cm；而ECOM2对应的STD显著增大，增幅为2~3 cm。可见，在经验光压模型基础上，增加BW模型似可降低激光残差STD，而加入地球辐射和天线辐射则未显示出显著的有益影响。
图 4 卫星轨道的SLR残差STDFig. 4 STD of SLR residuals of orbits with different force models
图选项
图 5统计了激光残差的系统偏差。可以看出，在NYS模式下，仅考虑经验光压模型时，不同厂家卫星的系统偏差符号相反；在经验光压模型基础上，考虑地球辐射和天线辐射，可显著降低SECM卫星的系统偏差，但却导致CAST卫星的系统偏差显著增大。例如，在S2EA的结果中，SECM卫星的系统偏差非常接近于0，而CAST卫星的系统偏差则超过7 cm。此外，增加BW模型或采用ECOM2，有助于降低激光残差的系统偏差。与S1的结果相比，对于CAST和SECM卫星，BS1EA*的系统偏差分别降低13%和88%，S2的系统偏差分别降低15%和45%，BS2EA*的系统偏差分别降低26%和100%。
图 5 卫星轨道的SLR残差的系统偏差Fig. 5 SLR residuals biases of orbits with different force models
图选项
与NYS模式的结果相比，在CYS模式下，对于CAST卫星，各模型得到的系统偏差都普遍增大；对于SECM卫星，采用ECOM2模型对应的系统偏差显著增大，ECOM1模型的结果则未见显著恶化。
3 结论
本文基于卫星几何、光学参数及天线辐射功率等信息，对BDS-3卫星非保守力进行建模，采用全球观测网的数据进行定轨试验，分析在不同的经验光压模型(ECOM1和ECOM2)基础上，增加BW光压模型、地球辐射和天线辐射模型对定轨结果的影响。结果表明，经验光压模型是影响轨道精度的主要因素，在名义偏航模式下，ECOM2模型具有更好的表现，但ECOM1模型对卫星的姿态模式更不敏感。在控制偏航模式下，ECOM2得到的轨道径向重叠互差RMS和激光残差STD显著增大。在经验光压模型基础上，增加地球辐射或天线辐射，可以显著降低SECM卫星激光残差系统偏差；对CAST卫星则适得其反，其原因还有待进一步分析。此外，增加BW模型对于提高轨道精度也是有益的。为了进一步满足BDS-3高精度应用需求，构建更加准确的非保守力模型，需要卫星厂家提供更详细和准确的卫星表面几何和光学参数等信息。
作者简介
第一作者简介：阮仁桂(1983—) 男 博士 助理研究员 研究方向为精密定轨、定位和时间同步。E-mail: rrg2002me@163.com
初审：张 琳复审：宋启凡终审：金 君资讯
