基于子集模拟法的卫星碰撞风险评估研究.pdf

JOURNAL OFNORTHUNIVERSITY OFCHINA(NATURAL SCIENCEEDITION)(Sum No.211)(总第2 11期）2023No.5Vol.442023年第44卷第5期中北大学学报（自然科学版）文章编号：16 7 3-3193（2 0 2 3)0 5-0 544-0 6基于子集模拟法的卫星碰撞风险评估研究邢尧，余红英，樊永生（1.中北大学电气与控制工程学院，山西太原0 30 0 51；2.中北大学大数据学院，山西太原0 30 0 51）摘要：卫星碰撞风险评估是卫星进行空间碎片预警和规避机动的基础，传统的蒙特卡洛法计算碰撞概率效率较低。本文将Gibbs采样与子集模拟法相结合，并应用于卫星碰撞概率计算中，通过选取合适的中间失效事件，将碰撞概率转换为一系列较大的中间失效事件条件概率的乘积，同时使用Gibbs抽样方法生成失效域内的样本，降低计算所需样本数。文中给出了子集模拟法求解卫星碰撞概率的具体步骤，并通过算例验证了子集模拟法的可行性和高效性。结果表明，子集模拟法在卫星碰撞概率计算中精度较高，且碰撞概率越低，相对于传统蒙特卡洛法的效率越高。关键词：蒙特卡洛法；子集模拟法；吉布斯采样；卫星碰撞中图分类号：V448文献标识码：Adoi:10.3969/j.issn.1673-3193.2023.05.009Research on Satellite Collision Risk Assessment Based onSubset Simulation MethodXING Yao,YU Hongying,FAN Yongsheng(1.School of Electrical and Control Engineering,North University of China,Taiyuan 030051,China;2.School of Data Science and Technology,North University of China,Taiyuan 030051,China)Abstract:Satellite collision risk assessment is the basis for space debris warning and avoidance maneu-vers by satellites,and the traditional Monte Carlo method to calculate collision probability is inefficient.In this paper,Gibbs sampling was combined with the subset simulation method and applied to the satel-lite collision probability calculation.By selecting suitable intermediate failure events,the collision prob-ability was converted into the product of a series of larger conditional probabilities of intermediate failureevents,and the Gibbs sampling method was used to generate samples in the failure domain to reduce thenumber of samples required for the calculation.The specific steps of the subset simulation method weregiven to solve the satellite collision probability,and the feasibility and efficiency of the subset simulationmethod were verified by arithmetic examples.The results show that the subset simulation method ismore accurate in the satellite collision probability calculation,and the smaller the collision probabilityis,the more efficient it is compared with the traditional Monte Carlo method.Key words:Monte Carlo method;subset simulation method;Gibbs sampling;satellite collision0引言随着空间飞行器和空间垃圾数量地迅速增加，卫星与卫星、卫星与空间碎片的碰撞风险越来越高。卫星碰撞风险评估对空间目标预警和规避机动都非常重要。碰撞概率计算是风险评估的重要判据。目前，相关文献已经提出多种碰撞概率计算方法。Chan1-21提出了等效圆域法，证明了在相同收稿日期：2 0 2 3-0 2-0 2作者简介：邢尧（1996 一），男，硕士生，主要从事仿真推演方面的研究。(1)545尧等）基于子集模拟法的卫星碰撞风险评估研究（邢（总第2 11期）坐标系下，将两个相互独立的空间目标的位置误差协方差矩阵相加，就可得到二者的联合误差协方差矩阵。王华等3 针对空间目标相遇期间的碰撞风险问题提出3种化简方法。1)通过等效的方法利用均匀密度分布替换概率密度函数；2 将二重积分降维处理，转化为曲线积分的形式；3)利用数值求解的方法。这3种方法通过简化碰撞过程来计算碰撞概率，算法简单，易于理解，但计算需要复杂的坐标变换，偏差较大，适用范围有限。Vries和Sabol等4-5 通过蒙特卡洛方法计算碰撞概率。为了减小蒙特卡洛方法的计算量，Dolado6引人了重要性采样，但重要性采样并不适用于多维情况。蒙特卡洛方法可以减小计算偏差，但仍存在计算效率低和所需样本量大等问题。本文针对卫星碰撞风险评估问题，采用蒙特卡洛方法避免计算过程中复杂的坐标变换，在此基础上，将Gibbs采样与子集模拟法相结合，以减少计算碰撞概率所需的样本数，从而提高计算效率。根据卫星运动模型，对子集模拟法获得的碰撞概率和计算时间进行评估，证明该方法在碰撞概率计算方面具有良好的效率和准确性。1碰撞概率计算的基本方法1.1碰撞概率计算中的假设卫星碰撞概率计算不但要考虑接近时刻与接近距离，而且还要考虑轨道误差协方差。两个空间目标发生碰撞时受力相对复杂，空间目标运动状态的改变会导致误差数据发生变化。为了简化碰撞概率计算过程，本文假设1)3)7 1）由于两个空间目标接近过程非常短暂，相对速度非常大，因此，在接近期间将二者的相对运动视为线性运动，速度没有不确定性。2）两个空间目标位置误差不相关，且服从三维高斯分布，误差椭球在相遇期间保持不变。3）两个空间目标均等效为半径已知的球体。当卫星与空间目标相对距离小于二者半径之和，即认为发生碰撞。1.2碰撞概率计算的提出对于轨道存在交集的两个空间目标，可利用轨道外推模型估计二者的接近时刻和接近距离。在计算过程中，综合考虑两个空间目标的状态分布，给出接近期间相对位置的状态概率分布，即联合概率分布。由于位置误差相互独立，且服从高斯分布，两个空间目标各自的协方差矩阵相加，就得到二者的联合误差协方差矩阵。又因为二者状态相互独立，为方便计算，将卫星与空间目标的半径相加，形成的联合包络球可被看作一颗卫星，将两者的联合误差看作空间目标的误差椭球。此时，卫星位置不存在误差，空间目标可看作空间中的一点。碰撞概率计算可以转化为计算空间中服从三维高斯分布的一点落入半径已知球体内的概率。以卫星质心为原点建立坐标系，各轴平行于地心惯性坐标系的相应轴，以避免计算过程中的坐标变换。根据上述分析，接近期间的包络球与误差椭球如图1所示。包络球的半径为卫星与空间目标半径之和，椭球体3个半轴可通过对联合协方差矩阵进行对角化得到。此时，卫星与空间目标的碰撞概率，可通过在半径为R的球体内对联合概率密度积分获得，如式（1)所示。1Pexpdrdydz,2（2 元）23式中：积分区域V为?+y?+0,发生碰撞。在接近时刻，根据三维高斯概率密度函数生成N个样本，如果相对距离小于阈值D，则认为发生碰撞。碰撞发生的样本数N。除以样本总数NT，即得到碰撞概率8 N。P(d*D)(4)NT2子集模拟法计算碰撞概率2.1子集模拟法基本原理子集模拟法是一种计算工程系统中遇到的高维小概率事件的随机模拟方法，其基本原理是通过引人合适的中间失效事件，即本研究中的碰撞事件，将较小概率的卫星碰撞问题转化为一系列较大条件概率的乘积9。由功能函数g（x）可定义失效区域F=(x：g（x)0)，通过引人一系列中间值blbzb3bm=0，将碰撞事件分为一系列具有包含关系的中间失效事件F=(x：g（x）b）（k=1，2,.,m)，此时有FiDF2D.DF=F，依据概率论，碰撞概率计算表达式为m-lP(F)=P(F)=P(FI)IIP(F+I F,),1=1(5)式中：P（F+1IF,)表示失效事件F，发生情况下，事件Ft+1发生的概率。该方法使用蒙特卡洛模拟法进行初始化，生成1级区域的样本。当确定了失效区域F1，并计算出失效概率P（F1），就通过马尔可夫链蒙特卡洛（MarkovChainMonteCarlo，MCMC)法生成区域2 的条件样本。然后，定位中间失效区域F2，并利用MCMC法生成相应样本。重复该过程，直到定位出卫星与空间目标发生碰撞的失效区域。子集模拟法的关键在于中间失效区域的选择10 。本文通过选择一个恒定的概率值po，并在P(Ft+1lF,）=p。的每个失效区域内寻找相应的临界值。b，为1级失效区域的临界值，类似于式（3），会出现3种情况，即bl，x 处于失效区域Fl。由于每次迭代时的条件概率都等于力。，因此，式（5)中的碰撞概率可以变为P(F)=P(Fm)=P(Fm I Fm-1)po-l(7)O2.2Gibbs 采样MCMC法采样常用的方法为M-H采样。M-H采样面向高维问题时存在两个缺点：1）由于需要计算接受率，导致算法收敛时间变长，计算效率较低；2）多数高维问题的条件概率分布容易得到，但联合分布不容易求得。因此，本文采用Gibbs采样。假设一个三维的概率分布元（1，2，3），对该分布进行Gibbs采样，具体过程为1)4)。1）确定目标分布元（1,2,3)或者所有特征对应的条件概率分布，设定状态转移次数n1和需要的样本数n2。2）随机初始化初始状态值（i，r 2，r 9)。3）从条件分布中采样：a）从条件分布P（l 2,9)中采样得到样本i；b）从条件分布P（2 lri,r)中采样得到样本r2；c）从条件分布P(3lzi,r2)中采样得到样本r3。4）步骤3)每重复n1次，输出一个样本，直到样本数达到n2，停止算法。2.3子集模拟法具体步骤1）设置1=1，通过蒙特卡洛法模拟生成N个服从三维高斯分布且相互独立的样本向量x，k=1，N。上标1表示样本处于1级区域12）计算N个样本x的响应值g（x），将响应值从小到大排列，取第（1一p。）N个响应值作为中间事件Fi=（x：g（x）b 1)的临界值，若b10，则跳到第6 步直接计算碰撞概率；否则，P(F1)=P(g()b,)=po。3）从落在F,-1（l=2，3，m）域内的poN个样本出发，使用马尔科夫链生成(1一p。)N个服从密度函数q(x|Fi+1)的条件样本，构成总共N个条件样本(x：k=1,2,N)。4）利用功能函数g（x），可以计算出N个条件样本对应的响应值，将得到的响应值从小到大排列，选择第（1一p。）N个响应值作为中间事件F,=(x：g（x)b l)的临界值bl。根据定义，相关条件概率为P(F,/F,-1)=P(g(x)b/g(x)bl-1)=po,则P(F,)=po。547尧等）基于子集模拟法的卫星碰撞风险评估研究（邢(总第2 11期）5）重复步骤3）和4），直到某一层的第（1一p。）N个响应值大于0，则令bm=0，结束自动分层。利用式（6）判断m层的N个条件样本是否落人失效区域F。6）分层结束后，m层中落人失效区域Fm（即卫星与空间目标发生碰撞)的样本数记为NF。生成的样本总数为N=N+(m-1)(1-p。)N,(8)式中：m是达到碰撞区域所需的层数。根据式（7)可知，碰撞概率为NFP(dD)=o-P(F.|FuI)=1N(9)式中：N是每个层级的样本总数；N是m层中处于失效区域F的样本个数。子集模拟算法的效率取决于条件失效概率力。和每个步骤的样本数量N。Zu e v 等12 对条件失效概率的选择进行了详细分析，表明力。的最佳选择为0.2。要使用的样本数量N取决于问题维度和预期失效概率13，如果失效区域是断开的，则样本必须足够密集。3算例分析本节将对算法的计算效率和准确性进行验证分析。算法采用MATLABR2021b软件编程实现，计算机运行环境为6 4位Windows10操作系统，Intel（R）C o r e（T M）i 7-116 5G 7 C PU 2.80GHz处理器，16 GB内存。3.1算例1测试算例1选择文献14 中的案例5，该案例中两个空间目标在低地球轨道运行，接近期间为线性相对运动，相关数据的处理也和本文相似。其中，包络球半径为10 m，相对位置误差协方差矩阵通过二者协方差矩阵相加，取左上角3X3部分得到。二者接近时刻状态信息和误差协方差矩阵均由该文献提供，计算结果如表1所示。表1中样本数为10 和10 的蒙特卡洛计算结果由文献14 所提供，子集模拟法的样本数表示每层计算使用的样本数和分层层数。经比较可知，基于Gibbs采样的子集模拟法使用较少的样本就可得到较高的精度，且随着样本数的大幅降低，计算效率也得到提升，因此，该方法在卫星碰撞概率计算方面具有较高的应用价值。表1计算结果Tab.1Calculation results方法样本数NT碰撞概率P误差/%1060.0443360000.2952蒙特卡罗法1080.0444989130.071210100.0444672572000.04460.2985子集模拟法720000.044.440.06133.2算例2本小节将对编号为432 0 的卫星进行风险评估，评估过程如图2 所示。首先，从Space-Track获取所有在轨目标数据。然后经过近地点/远地点筛选，发现该卫星可能与3533个空间目标发生碰撞，通过轨道外推模型计算3天内3533个空间目标与该卫星最接近距离，其中有3个空间目标与卫星最接近距离小于5km。以2 0 2 2 年11月10 日5时49 分37 秒为初始时刻，该时刻卫星与空间目标轨道根数如表2 所示。开始卫星数据近地点/远地点筛选轨道外推最接近距离 5km是计算碰撞概率结束图2 卫星风险评估流程图Fig.2Satelliteriskassessmentflowchart近地点/远地点筛选规则为ra+dra+d,(10)式中：ra、r p 分别为卫星的远地点和近地点；ra、r分别为空间目标的远地点和近地点；d为距离预留量，本文设为5km。当空间目标满足式（10）时，表示该空间目标对卫星总是安全的，可以永久剔除。548中北大学学学报（自然科学版）2023年第5期表2 轨道根数Tab.2Orbital elements半长轴/km偏心率倾角/（）升交点赤经/（）近地点幅角/（）平近点角（）卫星7833.9470.0028668101.543153.820574.744592.3534空间目标17963.7190.120039762.029.3144.0433103.228283.7846空间目标27.716.4470.020764582.4588326.0831246.3987323.7862空间目标37862.4510.001750874.0238271.6085290.5714243.1949设置仿真时间为1d，仿真步长为1s，经过轨道外推模型计算卫星与空间目标的接近时刻与接近距离。对接近事件分别采用不同样本数的蒙特卡洛法和子集模拟法计算碰撞概率。计算过程假设3次接近事件在接近时刻沿坐标轴方向的相对位置误差相等，分别为0.5km，1k m 和2 km，卫星与空间目标的联合包络球半径均为R=0.1km。接近信息及碰撞概率计算结果如表3所示。表3接近信息与碰撞概率计算结果Tab.3Proximity information and collision probability calculation results空间日标接近时刻/s接近距离/m方法总样本N计算时间t/s碰撞概率P误差/%1084.9471233.1060E-50.5221蒙特卡洛10956.8224913.1246E-50.0737139581.456.353.71010774.515.0833.1223E-54.9E+50.5553293.1364E-50.4516子集模拟4.9E+64.8987533.1192E-50.09931084.9357626.4660E-51.1660蒙特卡洛10952.3362786.5497E-50.11312118061 674.94411010712.3374626.5423E-54.2E+50.4741546.5627E-50.3118子集模拟4.2E+63.9829266.5379E-50.06731084.8937502.0400E-64.8290蒙特卡洛10965.0952842.2250E-63.802350.5084.686.805.31010735.7968442.143.5E-66.3E+50.5401772.0547E-62.9670子集模拟6.3E+66.894.0772.149.0E-60.2566算例结果显示，随着样本数的增加，蒙特卡洛法和子集模拟法的计算时间也随之增加。子集模拟法使用较少的样本就可得到较高的计算精度，但计算时间远远短于同精度下的蒙特卡洛法。当碰撞概率为10-5时，蒙特卡罗法需要10 个样本可使误差小于1%，而相同精度情况下子集模拟法只需数量级为10 的样本数，约为蒙特卡洛法的千分之一，计算耗时为蒙特卡洛法的十分之一。当碰撞概率为10-6 时，蒙特卡洛法使用10 个样本已不能得出高精度的结果，而样本数为10 1时，计算时间需十几分钟。当计算误差小于1%时，子集模拟法的样本数仍为10 量级，计算时间增长不明显。由此可看出，在计算卫星碰撞概率方面，子集模拟法的计算效率远高于蒙特卡洛法。4结论卫星碰撞概率计算是空间安全的一个重要研究内容，本文采用基于Gibbs采样的子集模拟法计算卫星碰撞概率，避免了碰撞概率计算过程中复杂的坐标变换。利用轨道外推模型获取卫星与空间目标的最接近时刻及接近时刻的各自状态，在子集模拟法计算卫星碰撞概率时，通过Gibbs采样降低计算所需样本数，提高计算效率。通过与蒙特卡洛法进行对比，证明在卫星碰撞风险评估中，基于Gibbs采样的子集模拟法在卫星碰撞概率计算方面具有良好的准确性，且比蒙特卡洛方法所需样本数少，计算效率较高。想要获取更高精度只需增加子集模拟法的样本数量，且在一定范围内计算效率不会出现明显变化。参考文献：1 J CHAN K.Collision probability analyses for earth-or-biting satellites J.Advances in the AstronauticalSciences，19 9 7,9 6:10 33-10 48.2(CHAN K.Analytical expressions for computingspacecraft collision probabilitiesC/The llth An-549尧等）基于子集模拟法的卫星碰撞风险评估研究（邢（总第2 11期）nual AAS/AIAA Space Flight Mechanics Meeting,2001:305-320.3王华，李海阳，唐国金。飞行器碰撞概率计算的一般方法J.国防科技大学学报，2 0 0 6，2 8（4）：2 7-31.WANG Hua,LI Haiyang,TANG Guojin.General meth-od for calculating spacecraft collision probability J.Journal of National University of Defense Technolo-gy,2006,28(4):27-31.(in Chinese)4 VRIES W D,PHILLION D.Monte carlo method forcollision probability calculations using 3D satellitemodelsC/Advanced Maui Optical and Space Sur-veillance Technologies Conference,2010:1119964.5 SA BO L C,BI NZC,SEG ERM A N A,e t a l.Pr o b a-bility of collision with special perturbations dynamicsusing the Monte Carlo methodJ.Advances in theAstronautical Sciences,2012,142:1081-1093.6 D O LA D O J C，LEG END RE P,G A RM I ER R,et al.Satellite collision probability computation forlong term encountersCJ/AAS-AIAA Astrodynam-ics Specialist Conference,2011:142.7 霍俞蓉，李智，韩蕾。基于拉普拉斯变换的空间目标碰撞概率计算方法.北京航空航天大学学报，2018，44(4):8 10-8 19.HUO Yurong,LI Zhi,HAN Lei.Computationmethod of collision probability of space object basedon Laplace transformationJ.Journal of Beijing U-niversity of Aeronautics and Astronautics,2018,44(4):810-819.(inChinese)8 K RU PENEV D,BO YA RK IN D,IA K U BO VSK IID.Improvement in the computational efficiency of atechnique for assessing the reliability of electric pow-er systems based on the Monte Carlo methodJ.Reliability Engineering and System Safety,2020,204:107171.9J MAESTRINI M,LIZIA P D.Guidance for autono-mous inspection of unknown uncooperative residentspace objectsC/2020 AAS/AIAA AstrodynamicsSpecialist Conference.2020:1-13.10宋述芳，吕震宙基于马尔可夫蒙特卡罗子集模拟的可靠性灵敏度分析方法J.机械工程学报，2 0 0 9，45(4):33-38.SONG Shufang,LU Zhenzhou.Structural reliabilitysensitivity analysis method based on markov chainmonte carlo subset simulationJ.Journal of Me-chanical Engineering，2 0 0 9,45(4):33-38.（i n Ch i-nese)11刘岗，封周权，华旭刚，等基于子集模拟法的风机叶片可靠度分析J.计算力学学报，2 0 19，36（5）：636-641.LIU Gang,FENG Zhouquan,HUA Xugang,et al.Reliability computation of wind turbine blades basedon subset simulation algorithmJ.Chinese Journalof Computational Mechanics,2019,36(5):636-641.(in Chinese)12 ZUEV K M,BECK J L,AU S K,et al.Bayesianpost-processor and other enhancements of subsetsimulation for estimating failure probabilities in highdimensionsJJ.Computers&Structures,2012,92:283-296.13王家，金艳琪，陈雅含基于子集模拟的预制构件生产目标均衡优化算法研究J.公路工程，2 0 2 2，47(5):189-196.WANG Jia,JIN Yanqi,CHEN Yahan.Research onoptimal algorithm of prefabricated components pro-duction target balance based on subset simulationJ.Highway Engineering，2 0 2 2，47（5)：18 9-196.（i nChinese)14 ALFANO S.Satellite conjunction Monte Carlo analy-sis J.Advances in the Astronautical Sciences,2009,134:2007-2024.