天基激光移除空间碎片仿真平台研究与开发

杨武霖 ，陈川，余谦，李明，龚自正

1．北京卫星环境工程研究所，北京 100094

2．北京卫星环境工程研究所 可靠性与环境工程重点实验室，北京 100094

3．中国空间技术研究院，北京 100094

尺寸在1～10cm的厘米级空间碎片，因其难以跟踪编目和被动防护，对航天器的在轨安全构成了严重威胁。空间碎片之间的相互撞击则会产生更多的碎片，造成碎片数量的急剧增加。空间碎片主动移除技术，是应对厘米级碎片对航天器和空间环境严重影响的最有效途径。在众多的主动移除技术中，天基激光移除空间碎片技术因其激光光斑可调、作用距离远、反应迅速、可多次使用的优点，获得广泛关注。其主要技术原理是激光辐照空间碎片表面熔融、气化、电离产生等离子体，等离子体等气化产物反喷使碎片获得速度增量，改变碎片的轨道，经过多次对碎片的探测跟踪、烧蚀驱动和变轨评估过程，使碎片进入稠密大气层烧毁，实现碎片移除。

天基激光移除空间碎片仿真主要用于移除方案的效果分析和不同因素下对移除效果影响分析。在移除方案的效果分析方面，始终致力于天基激光移除空间碎片技术研究的Phipps在2016年提出了采用紫外激光移除低轨碎片的天基移除方案［1］。方案的移除效果仿真分析采用单脉冲共面变轨过程中的能量守恒原理，变轨过程中的激光所需能量即碎片移除前后轨道上的能量变化，不考虑碎片和移除平台具体的变轨过程。国内仿真分析了不同影响因素下的降轨效果。王成林等采用单脉冲共面变轨模式，通过分析半长轴和偏心率的变化趋势，对碎片轨道不同位置上的速度分量的变化对其新轨道的近地点高度的影响进行了研究［2］。康博琨等将天基激光驱动空间碎片的驱动力当做碎片轨道摄动力，仿真分析了驱动力作用下碎片轨道六要素的变化趋势，分析了特定碎片群的移除效果［3］。方英武等也将激光的驱动力当做碎片轨道摄动力，仿真分析了激光脉冲数对碎片半长轴、偏心率和倾角的变化趋势［4］。张品亮等基于轨道力学和激光与物质相互作用理论建立了激光移除空间碎片的三维变轨模型，在考虑速度增量分量对轨道倾角影响的情况下，仿真分析了不同材料的降轨效果［5］。温泉在单脉冲共面变轨的理论基础上，通过近地点高度的变化情况，仿真分析了碎片旋转情况下的降轨过程［6］。

目前国内外开展的数值仿真研究，主要针对特定轨道或区域的空间碎片，采用简单的变轨模型和kJ级高能激光器。但是部署在太空的移除平台在运行过程中并非固定不动，需要根据碎片和移除平台的实时相对位置判定是否能够烧蚀驱动特定碎片，规划移除平台的工作窗口，且并非所有特定的碎片都能够被激光驱动或移除。其次，当前的变轨模型大都以霍曼转移、大速度下的单脉冲共面变轨或从能量守恒的角度出发，只关注初始状态和结束状态时的轨道信息，忽略了激光烧蚀驱动过程中的轨道变化。此外，当前的方案设计中大都采用KJ及高能激光器，这要求单次脉冲下产生的速度增量在80～100m/s之间，但当前的天基激光器技术水平无法满足要求，需要对低能量激光器的可行性进行研究。开发一种面向所有选定碎片、考虑碎片和移除平台的实时状态和过程变化、基于实际硬件水平的三维数值仿真平台，对定量研究移除过程的碎片和移除平台的状态变化、从整体角度评估移除效率和效果具有重要意义。

针对以上问题，本文首先基于轨道动力学和激光烧蚀驱动空间碎片力学效应建立了天基激光移除空间碎片三维变轨模型。其次从总体设计出发，分析了三维数值仿真平台的需求和总体框架，设计各模块功能和数据接口，采用C++/Qt开发了三维数值仿真平台。最后基于当前的激光器技术水平，分析计算了其烧蚀驱动的动力学效果。基于当前碎片分布特性，输入多个碎片和移除平台数据，展示了多个移除平台协同工作的碎片移除过程和效果，对三维数值仿真平台的有效性进行了验证。

该仿真系统的研究与开发，解决了常用的单脉冲共面变轨和轨道摄动过于理想、无法分析碎片与移除平台所有相对位置关系的局限性。其次，该系统计算和分析碎片和移除平台的位置和速度，克服了从能量角度无法精细刻画移除过程的劣势。最后，仿真系统的平台和碎片输入参数可调，可对低能量激光器和近距离的激光有效作用距离进行有效的仿真分析，为发展移除卫星星座奠定基础，同时可在当前硬件技术条件下加速工程化推进。

1 三维变轨模型

三维变轨模型是根据碎片和移除平台的轨道六要素初始信息，移除平台的激光有效作用距离、频率和单脉冲能量等参数，模拟计算碎片和移除平台在激光烧蚀驱动过程中的实时状态变化，并在移除任务结束后判定各碎片的最终降轨状态和平台的效率及效果。

首先需要根据碎片和移除平台的轨道六要素信息，判定二者最短相对距离Dmin与激光有效作用距离R的大小。激光有效作用距离即碎片在该距离范围内能够被激光辐照并烧蚀驱动产生速度增量，在该距离范围之外则激光无法烧蚀驱动碎片。对于最短相对距离小于激光有效作用距离的碎片，在二者的运动过程的某一时刻，碎片将会进入激光烧蚀驱动的区域。当二者的实时相对距离D小于激光有效作用距离时，激光烧蚀并驱动碎片产生速度增量，改变碎片的速度。碎片在激光烧蚀区域内时，启动激光持续烧蚀碎片以改变碎片的速度和位置，直至二者的实时相对距离再次大于激光有效作用距离。当碎片运行一个周期后再次进入激光烧蚀驱动区域，重复上述烧蚀过程，直至烧蚀后的碎片与移除平台的最短相对距离大于激光有效作用距离，或者烧蚀后碎片近地点满足移除条件。因此，对于所有输入的碎片，其结果是被移除或此次任务无法移除两种结果。图1为模型的计算流程。

图1 变轨模型流程Fig.1 Procedures of orbital transfer model

2 三维变轨模型

2．1 碎片和移除平台的初始状态

当前关于空间碎片的轨道信息，包括国内外在碎片移除系统方案设计和仿真中采用的碎片轨道数据，大都来自于美国NORAD卫星星历，又称双行轨道数据(Two-Line Orbital Element，TLE)。该数据主要用于描述空间目标的发射信息和轨道信息。格式如下:

1)30323U 07003A 07067．68277059 ．00069181 13771-5 44016-2 0 587

2)30323025．0330 358．9828 7594216 197．8808 102．7839 01．92847527 650

其中，轨道要素信息的数据如表1所示。

表1 TLE数据中的轨道要素信息Table 1 Orbital parameters of the TLE data

为方便轨道六要素和位置速度之间的相互转换，简化后续的位置判定和速度增量计算，常采用的轨道六要素为动量矩、偏心率、轨道倾角、升交点赤经、近地点幅角和真近点角。

(1)TLE数据中每天环绕地球的圈数n转换为轨道六要素中的动量矩h

首先由圈数求解空间目标的运行周期T:

式中:T0为地球自转周期，T0=86164．09053s。随后求解空间目标的半长轴a:

最后由

获得空间目标动量矩h。

式中:μ为地球引力常数，μ=3．986×105km3/s2;e为轨道偏心率，由TLE数据可直接获取。

(2)平近点角M转换为真近点角f

首先将平近点角M转换为偏近点角E，平近点角到偏近点角常用迭代法进行解算［7］。具体做法是把开普勒方程写成E=M+esinE形式，迭代的起步是将方程右边的E用M代替，求出偏近点角的第一次近似值，然后再进行若干次迭代求出满足要求的近似值。

式中:e为轨道偏心率。由式(4)将偏近点角E转换为真近点角f。

经过上述转换关系，可将TLE数据中的轨道信息转换为轨道六要素，用于记录三维数值仿真系统中空间碎片和移除平台的轨道信息。

2．2 最短相对距离Dmin

在已知空间碎片和移除平台初始轨道六要素的条件下，比较碎片的运行周期TD和移除平台的运行周期TL，以长周期的运行时间作为最短相对距离判定的采样周期T=Max{TD，TL}。以固定的时间步长Δts计算在采样周期内碎片和移除平台的相对距离，选取相对距离最小值作为二者运行过程中的最短相对距离Dmin。

轨道六要素转换为地心赤道惯性系下位置r={rx，ry，rz} 和速度 v={vx，vy，vz} 信息的关系如下。

首先将轨道六要素转换为近焦点坐标系下的速度v珔X和位置r珔X信息:

近焦点坐标系中的位置r珔X和速度v珔X信息转换为地心赤道惯性坐标系下的位置r和速度v信息为:

式中:Q珋xX为转移矩阵［8］。

2．3 碎片和移除平台轨道随时间的变化关系

碎片和移除平台的位置和速度随时间的变化关系可由泰勒级数展开式近似计算［9］。如果已知t0时刻的位置r0和速度v0，可由拉格朗日系数f和g及其一阶导数，根据下述表达式求出任意时刻的位置r和速度v:

利用全局变量χ和斯达姆夫函数C(Z)与S(Z)所表示的拉格朗日系数如下:

其中半长轴的倒数α为:

对于椭圆轨道，α ＞ 0。如果已知Δt、r0、v0和α，可从全局开普勒方程中解出全局近点角χ，进而求解得出Δt时刻后空间目标的位置和速度信息。

2．4 实时相对距离D

已知t时刻碎片在地心赤道惯性坐标系中的位置矢量rD为rD={xD，yD，zD}，移除平台此时的位置矢量rL为rL={xL，yL，zL}，则此时二者的相对距离D为:

2．5 激光烧蚀驱动碎片

激光辐照空间碎片表面时，熔融气化产生的反喷物质运动使碎片获得速度增量，冲量耦合系数反映了烧蚀过程中激光能量转化为碎片冲量的能力。在确定激光烧蚀驱动碎片获得的速度增量时，激光与碎片的冲量耦合系数Cm是主要参数，

式中:m为碎片质量，Δv为速度增量。

在已知碎片质量、激光能量和烧蚀过程的冲量耦合系数条件下，烧蚀过程中产生的速度增量Δv为:

速度改变后的轨道六要素变化，可根据二体运动和轨道六要素的定义计算［10］。

2．6 碎片移除判据

当激光多次烧蚀驱动碎片，使碎片再也无法进入激光烧蚀驱动区域时，根据此时的轨道要素计算其近地点高度H。

式中:RE为地球半径。当近地点高度低于200km时，碎片将快速坠入大气层烧毁，实现碎片移除［11］。

3 仿真平台设计

3．1 仿真平台需求分析

仿真平台的定位是采用真实碎片数据模拟仿真碎片移除过程中碎片与移除平台的状态变化，为满足不同需求的移除任务的规划和优化提供支持。仿真平台的基本功能是:在输入初始碎片和移除平台信息条件下，采用激光烧蚀驱动变轨模型，针对特定碎片和移除平台采用三维效果模拟驱动变轨过程;针对分布在不同区域的大量碎片和多个移除平台，模拟多个碎片区域的移除效果和多移除平台的移除能力。

因碎片和移除平台的数据数量相对较大，仿真平台应具有从文件或数据库中直接导入相关初始信息的功能。平台针对显示移除过程细节和移除后的轨道和移除平台能力评估的需求，需要仿真平台的视角随时间变化逐渐切换。对于移除过程中进入驱动范围、正在驱动和驱动后已经降轨移除的碎片，采用不同的标记方法在三维场景中进行展示。

3．2 系统框架设计

根据需求分析，采用4层构架模型对天基激光移除空间碎片三维数值仿真平台进行系统框架设计。如图2所示，平台分为通用层、服务层、业务层和表示层。通用层为服务层、业务层及表示层组件提供开发测试和运行环境，属于公共服务功能域，包括操作系统、开发工具(QT、C/C++、OSG)、测试工具等商业软件。服务层为业务层和表示层提供所需基本数据采集、协议转换、协同分布处理等服务，包括仿真推演组件。业务层为表示层提供专业计算分析服务，包括轨道计算分析组件。表示层为满足不同场景表达需求而提供多种场景视觉呈现和用户交互服务，包括三维可视化组件。

图2 三维数值仿真平台体系结构Fig.2 Structure of the 3D simulation platform

系统主要包括两大模块:空间碎片态势推演模块和轨道计算分析模块。模块功能与工作流程如图3所示。

图3 模块功能与工作流程Fig.3 Functions and workflows of the models

(1)空间碎片态势推演模块

该模块主要根据导入的天基激光器和空间碎片轨道信息，显示基本空间态势，推演并显示碎片被激光照射后的轨道演化过程。主要功能如下:

1)批量空间碎片导入:以文件方式导入以上空间碎片的轨道参数(碎片初始化参数为轨道六要素、碎片的冲量耦合系数Cm和质量m)。

2)批量天基激光器导入:以文件方式导入激光器参数(天基激光器的参数为轨道六要素、有效作用距离R和激光能量E)。

3)基本空间态势显示:显示导入空间碎片的空间运动轨迹，显示空间运动轨迹及其作用范围，显示地球、星空等背景。

4)接近判定和清除过程显示:实时判定并高亮显示位于天基激光作用范围内的空间碎片;进入天基激光作用范围内的碎片，激光每秒钟作用空间碎片1次(次数可在控制面板上手动设置)。判定作用后的碎片是否低于200km(该判据可在控制面板上手动设置)，当碎片高度低于200km时自动消失。

5)移除结果输出显示:在面板上用表格显示目标的状态信息(位置和速度信息)，激光平台、尚在运行的碎片、已移除碎片采用不同的颜色和字体在表格中显示和标注，同时可以输出清除结果的文本文件。

(2)轨道计算分析模块

该模块主要根据导入的天基激光器和空间碎片轨道信息，完成空间目标星历预报，解算分析空间目标相互位置关系，解算分析在外力作用下(由用户提供作用力大小和方向)空间碎片的轨道变化。主要功能如下:

1)空间目标星历预报:计算J2000．0惯性系位置速度、地固系位置速度、J2000惯性系开普勒要素预报。

2)空间目标相互位置:解算空间目标的相互位置关系。

3)空间碎片轨道演化:能够完成在外力作用下空间碎片轨道变化的解算。

4)计算结果支持文件输出。

5)具备参数配置完成情况下自动运行的能力，同时具备人工干预能力。

4 仿真平台开发及仿真运行效果

采用 C++/Qt作为开发语言，基于 Visual Studio 2016进行了软件开发。首先根据仿真平台的总体设计，将3．2节中描述的功能转化为C++代码描述。整个仿真系统的界面由多个窗体组成，利用可视化工具可分别创建不同场景，并对场景中的元素进行定制化设计。

基于该仿真平台，选取两个场景进行仿真模拟:1)分别从碎片和移除平台的角度，刻画天基激光移除空间碎片的实时过程。上述仿真可用于碎片移除方案的优化。2)多个移除平台构成的星座移除给定区域内空间碎片的过程及其效果。

上述仿真可用于区域内碎片移除策略及优化分析。

(1)天基激光移除特定空间碎片的实时过程

在仿真模拟天基激光移除特定空间碎片的实时过程的示例中，移除平台与碎片位于不同的倾角和轨道高度上，且移除平台轨道略高于碎片轨道。假设激光驱动过程中碎片获得的速度增量沿移除平台指向碎片方向。移除平台和碎片参数如表2～表5所示。

表2 碎片冲量耦合特性Table 2 Physical properties of the debris

表3 碎片轨道参数Table 3 Orbital elements of the debris

表4 移除平台激光特性Table 4 Laser properties of the removal platform

表5 移除平台轨道参数Table 5 Orbital elements of the removal platform

仿真表明:从初始位置开始，第1179 s时碎片进入移除平台的激光有效作用范围内，激光烧蚀驱动空间碎片长达35s后，碎片离开激光有效作用范围。184s后，碎片进入大气层烧毁，碎片被移除。整个过程如图4所示，蓝色圈表示激光有效作用范围。

(2)针对特定区域的多个碎片采用多个移除平台的移除过程

针对碎片分布密集的轨道高度在800 km的区域，任意选择了698个碎片，并在其附近部署了21个天基激光移除平台，激光特性如表6所示。仿真结果显示在1个月的任务周期内，被移除的碎片数目仅为1个。分析表明，碎片移除的成功与否与移除平台的部署位置和激光特性密切相关，采用霍曼转移等简单的变轨分析并不适用于特定区域或者整个空域的空间碎片等多个碎片移除场景。仿真过程如图5所示。

图4 天基激光移除特定空间碎片的实时过程Fig.4 Real-time steps for specific debris removal

表6 激光特性参数Table 6 Input laser parameters

5 结束语

图5 针对特定区域的多个碎片采用多个移除平台的移除过程Fig.5 Real-time steps for multiple debris removal with multiple removal platforms for specific area

针对当前天基激光移除空间碎片仿真中变轨模型过于理想、变轨过程描述不足和当前硬件技术无法满足系统设计参数的问题，研究和设计了一种基于真实的可观测空间碎片数据的一套天基激光移除空间碎片三维数值仿真平台，仿真平台运行达到预期目标。研究的不足之处在于碎片移除过程中仅考虑激光烧蚀驱动过程，忽略了移除平台的探测、跟踪过程;同时激光驱动碎片变轨模型的速度增量方向固定，未考虑碎片形状和旋转对速度增量的影响。上述问题将在后续的研究中展开，不断丰富和优化仿真系统的性能。