CN105868503B - 地基激光移除空间碎片过程的三维建模与仿真方法 - Google Patents

Abstract

本发明公开了一种地基激光移除空间碎片过程的三维建模与仿真方法，基于地心赤道坐标系中碎片的初始位置和速度矢量，以及激光站的位置矢量，激光的波长、频率等参数，计算模拟地基激光移除空间碎片的物理过程。该方法极大的简化了计算方法，便于参数设置和显示输出结果，提高了其准确度和使用灵活性。

Description

地基激光移除空间碎片过程的三维建模与仿真方法

技术领域

本发明涉及航天技术领域，尤其是涉及一种地基激光移除空间碎片过程的三维建模与仿真方法。

背景技术

空间碎片是指人类空间活动所产生的太空垃圾。目前地球周围的空间碎片总量已经接近7000t，其中有43％主要分布在400-2000km的低地轨道上。它们严重威胁着航天器的安全，减缓空间碎片已势在必行，如何清除空间碎片成为航天领域关注的热点和前沿问题。

激光移除空间微小碎片是国际上一致认可的理想碎片移除方法。该方法采用激光推进的方式，当高强度脉冲激光照射在碎片表面时，使表面材料等离子化产生高温气体，产生一系列类似于火箭反推的热物质射流，给碎片提供速度增量来降低近地点高度，使其在较短时间内落入大气层烧毁。

在激光移除空间碎片过程建模方面，已有学者对其进行了研究。Phipps等建立了二维变轨模型，通过给定的激光和碎片参数能够预测碎片的降轨高度(参见ClaudeR.Phipps,Kevin L.Baker,Stephen B.Libby,et al.Removing orbital debris withlasers.Advances in Space Research,2012,49:1283–1300)。金星等(参见:金星,洪延姬,常浩.地基激光清除椭圆轨道空间碎片特性的计算分析.航空学报,2013,34:2064-2073)建立了以速度增量为初始量的单脉冲和多脉冲碎片变轨二维模型，该模型仅能分析二维平面内空间碎片的降轨过程。韩威华等(参见:韩威华,甘庆波,何洋,杨新.天基激光清理低轨空间碎片的最佳角度分析与过程设计.航空学报,2015,36:749-756)针对高能脉冲激光清除低轨空间碎片的问题，建立了高能脉冲激光清除空间碎片作用过程和轨道演化的二维模型，推导了脉冲激光作用碎片的最佳角度与碎片轨道参数的解析关系。常浩等(参见:常浩,金星,洪延姬.基于激光烧蚀冲量耦合的空间碎片主动变轨建模与仿真.航空学报,2013,34:2325-2332)针对空间碎片主动离轨防止碰撞在轨卫星的情况，基于轨道动力学建立了空间碎片主动变轨物理过程三维模型。该模型的计算方法十分复杂，以碎片轨道根数作为初始量，碎片进入激光作用范围的判断、碎片轨道参数、速度增量等需要在地心惯性坐标系、轨道坐标系和地心赤道坐标系之间进行转换，这无疑增加了系统的复杂程度和计算量，并且忽略了碎片逃逸等实际问题。

激光移除空间碎片是一个包括多脉冲激光冲量耦合作用、碎片变轨、陨落和逃逸等的复杂过程。虽然目前已经建立一些激光清除碎片过程的轨道动力学模型，但是它们大多数为二维模型，并且以速度增量为初始参数，而忽略了激光与物质的冲量耦合作用、碎片轨道参数变化和轨道倾角变化对速度增量的影响。在已建立的三维模型仅针对某种特定的变轨情况，在激光移除空间碎片领域具有一定的局限性，忽略了真实空间环境下的一些实际问题，并且计算方法十分复杂，结果信息输出不够详尽与直观。因此，研究一种具有普适性的、计算模型简单、能够描述激光移除空间碎片的真实物理过程的三维模型对定量分析激光清除空间碎片的过程具有重要意义。

发明内容

鉴于现有技术的上述缺陷，本发明旨在提供一种具有普适性、计算模型简单、能够描述地基激光移除空间碎片真实物理过程的三维仿真方法，建立包括激光冲量耦合作用、碎片变轨、陨落和逃逸等关键物理过程的仿真模型，真实、准确的输出激光移除空间碎片的物理过程，为激光参数选取和碎片移除的可行性判断提供基础数学模型。

本发明需要解决的技术问题包括：(1)如何采用简单的计算方法将地基激光移除空间碎片中的一系列关键问题集成在一个三维模型中；(2)如何避免在不同坐标系之间转换来判断碎片是否进入激光作用范围；(3)如何方便的设置参数和显示输出结果。

本发明采用了如下的解决方案:

一种地基激光移除空间碎片过程的三维建模与仿真方法，包括以下步骤：

(1)建立地心赤道坐标系，以空间碎片和激光站的位置矢量(x,y,z)和速度矢量(v_x,v_y,v_z)作为输入参数，建立M函数rab_check.m，判断碎片是否在激光的作用范围内；其中，基于地心赤道坐标系(单坐标系)，采用激光站地表切面，对是否具备清除条件进行判断，假设激光站坐标为(x₀,y₀,z₀),碎片的位置矢量为(x₁,y₁,z₁)，地球球面方程为：x²+y²+z²＝6378²，根据几何方法确定激光站点在地球表面的切面方程：x₀x+y₀y+z₀z＝6378²；只有在碎片移动到激光站上方时(x₁x₀+y₁y₀+z₁z₀-6378²>0)，才具有清除条件，此时向量夹角即为天顶角；

(2)建立计算轨道根数的M函数coe_from_sv，将采用地心赤道坐标系中的位置矢量和速度矢量作为输入参数，根据二体运动方程和牛顿定律转换为轨道根数；

(3)建立轨道与时间关系的M函数rv_from_r0v0，对碎片运行轨迹进行计算和信息记录，利用此函数实时跟踪碎片的移动过程，更新碎片的轨道信息；

(4)确定激光参数，基于激光在大气中的传输模型和激光与物质相互作用机理，建立碎片速度增量矢量M函数dv，其中包含常见空间碎片材料的相关参数，对常见空间碎片材料进行计算，确定激光脉冲作用后碎片的速度V_to_dV，对激光作用次数进行计数，计算激光作用后靶的质量m＝m₀-μE，μ为激光烧蚀率，E为靶上的激光能量；

(5)激光作用判断，判断激光脉冲作用后碎片轨道近地点是否降低，如果降低则激光作用计数+1，更新碎片速度矢量和轨道根数；反则，碎片按照时间步长继续运动；

(6)每个时间步长判断碎片是否移动出激光作用范围，如果移动出激光作用范围，则计算时间步长增加十倍，以提高计算效率，直至碎片再次进入天顶角90°警戒范围时，时间步长恢复，开始下一个计算周期，并且碎片过顶次数+1；

(7)当碎片近地点高度降至≤200km时，认为碎片坠入大气层烧毁，即成功移除；

(8)在计算结果显示输出方面，模型计算和储存每个时间步长的所有参数信息，用户选择输出所需参数和格式，基于Matlab GUI平台建立了图像显示输出模块，通过该模块可实现碎片变轨过程图形可视化输出，从不同视角采用动画的方式呈现碎片的移除过程，并且对关键位置进行突出化显示。

其中，天顶角的警戒范围为45-90°。

其中，只有当天顶角在-45°～45°范围内才能使从地面发射的激光辐照在空间碎片上。

其中，在三维空间中，存在轨道倾角，激光作用产生的速度增量分量，使碎片轨道倾角增大，逃逸出激光站可管辖的区域，进一步增加了对碎片逃逸情况的判断。

与现有技术相比，本发明首次提出了一种地基激光移除空间碎片的三维仿真方法。该方法打破了以往仿真方法的局限性和传统思路，极大的简化了计算方法，便于参数设置和显示输出结果，提高了其准确度和使用灵活性。基于Matlab建立了物理过程三维模型，该模型能够根据输入激光和碎片参数，计算碎片的变轨过程，根据需要实时输出碎片的轨道信息和清除结果，并且实现图形化输出。考虑到速度增量分量对轨道倾角的影响，该模型增加了碎片逃逸情况的判断，具有更广泛的适用性。

附图说明

图1是本发明的地基激光移除空间碎片过程的三维建模与仿真方法中地心赤道坐标系和轨道根数示意图；

图2是本发明的地基激光移除空间碎片过程的三维建模与仿真方法中三维仿真计算方法流程；

图3是本发明的地基激光移除空间碎片过程的三维建模与仿真方法中计算结果显示输出示意图。

具体实施方式

下面结合附图对本发明的地基激光移除空间碎片三维仿真方法进行进一步说明，该说明仅仅是示例性的，并不旨在限制本发明的保护范围。

参见图1，图1是本发明的地基激光移除空间碎片过程的三维建模与仿真方法中地心赤道坐标系和轨道根数示意图，其中，本发明的方法，基于地心赤道坐标系中碎片的初始位置和速度矢量，以及激光站的位置矢量，激光的波长、频率等参数，计算模拟地基激光移除空间碎片的物理过程。

首先，需要判断碎片是否具有清除条件，即：碎片的位置是否在激光的作用范围内。如果不在作用范围内，碎片按照脉冲时间间隔t开始移动，直到碎片进入激光作用范围内，开始计算脉冲激光与碎片的相互作用。依据相互作用机理，得到单脉冲变轨后碎片的状态矢量，计算变轨后的轨道根数；计算时间t后，碎片具有的位置和速度矢量，开始下一轮计算；重复计算，直到近地点高度小于设定值(即清除成功)，或者碎片超出激光作用范围。对于后者，根据碎片的移动，重新判断碎片是否具有清除条件，开始二次或多次过顶计算直到碎片被清除，或者逃逸出具备清除条件的轨道，无法继续清除。

2.1参数输入

本发明以空间碎片和激光站的位置矢量(x,y,z)和速度矢量(v_x,v_y,v_z)作为输入参数，调用M函数rab_check_new对输入参数的合理性进行判断。首先基于轨道力学对碎片位置矢量和速度矢量的是否匹配进行判断；再对碎片轨道进行分析，判断是否需要进行移除操作；然后对激光站位置矢量进行分析，判断激光站是否在地球表面；最后，分析碎片的初始位置是否在激光作用的区域内，并判断其是否具备移除条件。

2.2轨道根数的确定

在地心赤道坐标系中(参见图2，图2是本发明的地基激光移除空间碎片过程的三维建模与仿真方法中三维仿真计算方法流程)，X轴指向春分点的方向，XY平面为地球的赤道平面，Z轴与地球的旋转轴一致，且指向北。单位矢量I,J和K满足右手定则。在给定时刻，空间碎片的初始状态向量速度v₀和位置r₀，在地心赤道坐标系中状态向量可表示为：

r₀＝x₀I+y₀J+z₀K (1)

v₀＝v_x0I+v_y0J+v_z0K (2)

在该点受到激光辐照后，获得速度增量

Δv＝Δv_xI+Δv_yJ+Δv_zK (3)

变轨后速度矢量和位置矢量分别为：

v＝v₀+Δv＝(v_x0+Δv_x)I+(v_y0+Δv_y)J+(v_z0+Δv_z)K (4)

r＝r₀ (5)

根据二体运动方程和牛顿定律可以计算出变轨后的轨道根数：比角动量的模h、轨道倾角i、升交点赤经Ω、偏心率e、近地点幅角ω和真近点角θ。

2.3轨道与时间的函数

如果已知t₀时刻的位置r₀和速度v₀，可由拉格朗日系数f和g及其一阶导数，根据下述表达式求出任意时刻的位置r和速度v：

r＝fr₀+gv₀ (6)

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

其中长半轴的导数α为：

对于椭圆轨道，α＞0.

如果已知Δt、r₀、v₀和α，可从全局开普勒方程中解出全局近点角χ，具体步骤如下：

计算χ₀合理的初始估计值：

以χ₀为初始数据进行迭代计算：

其中vr₀为v₀在r₀方向上的投影,vr₀＝(r₀·v₀)/r₀。

算出比值ratioi＝f(χ_i)/f′(χ_i)。如果|ratio_i|超出精度范围10^-8，则按χ_i+1＝χ_i－ratio_i重新选取χ估计值，计算式(14)(15)直到|ratio_i|＜10^-8，则接受χ_i作为解。这样就可以根据式(6)(8)(9)得到Δt时刻的位置矢量r，根据式(7)(10)(11)得到Δt时刻的速度矢量v。

2.4判定激光作用条件

在地基清除系统中，由于受到大气的影响，理论上只有当天顶角在-45°～45°范围内才能使从地面发射的激光辐照在空间碎片上。但是在计算中，向量夹角的范围在0～180°，需要根据碎片的坐标和象限来计算天顶角，特别是对于任意激光站位置和空间碎片轨道的情况，这无疑增加了计算的复杂程度。本文采用激光站地表切面的方法，对是否具备清除条件进行判断。

假设激光站坐标为(x₀,y₀,z₀),碎片的位置矢量为(x₁,y₁,z₁)，地球球面方程为：x²+y²+z²＝6378²，根据几何方法可确定激光站点在地球表面的切面方程：x0x+y₀y+z₀z＝6378²。

只有在碎片移动到激光站上方时(x₁x₀+y₁y₀+z₁z₀-6378²>0)，才具有清除条件，此时向量夹角即为天顶角。这样就可对碎片的天顶角进行精确判断，而无需考虑碎片和激光站的象限问题。实现方式为：

elseif(R(1)*RL(1)+R(2)*RL(2)+R(3)*RL(3)-6378^2>0)

2.5激光与碎片相互作用

在激光传输方面，本发明采用简化模型进行分析。为了在较大距离上将激光传输到靶上，需要数KJ的激光脉冲能量和大型发射镜来克服光的衍射：

d_s＝aM²λz/D_eff (20)

其中d_s为靶上光斑直径；M²为光束质量因子(≥1，1代表最好)；Deff为计算衍射时输出孔径D中的有效照射光直径；a为一个乘数因子，对于高斯光束a＝4/π。

在确定速度增量时，激光与材料的动量耦合系数是主要参数：

C_m＝pτ/Φ＝mΔv/E(N·s/J) (21)

其中，激光作用后靶的质量m＝m₀-μE；μ为激光烧蚀率；p为一个强度为I、脉宽τ的激光脉冲在靶上产生的烧蚀压；激光通量为Ф；E为靶上的激光能量。目前已经通过实验获得了多种材料的C_m值。在计算模型中，我们对激光辐照次数进行记录，每次辐照后碎片的质量变化都能够计算，从而获得精确的速度增量。具体实现方式如下：

Fl＝(W*(Deff^2)*4*Teff)/(pi*(M^4)*((a*lamda*z)^2))；％作用在碎片上的激光通量％

E＝Fl*(pi*((dt/2)^2))；％作用在碎片上的能量％

m0＝density*4/3*pi*((dt/2)^3)；％初始质量％

m＝m0-(count+1)*u*E；％作用count次后的质量％

dv＝Cm*E/m/1000；％速度增量％。

2.6碎片逃逸

在三维轨道计算中，速度增量可能与碎片轨道不在同一个平面上，这样就会产生速度增量改变碎片的轨道倾角，当轨道倾角增加到一定程度时，碎片就逃逸出激光的作用范围，无法继续作用。在地基系统中，如果碎片沿当前轨道移动一周后仍然不具备激光辐照的条件，则认为碎片逃逸。

2.7计算结果显示输出

模型记录每个时间步长下的碎片参数，根据用户需要将计算结果以TXT文本格式输出，方便后期处理。调用plot函数对计算结果进行图形化输出，并对关键点进行标注，具体实施方式如下：

标注碎片初始位置：plot3(rn0(1),rn0(2),rn0(3),'r.',

'MarkerEdgeColor','r','markerFaceColor','r','markersize',20)；

绘制碎片轨道：plot3(rx,ry,rz,'b')；hold on

选择视角：view([-27 17])；

实现动画输出：comet3(rx,ry,rz)；

表1为计算实例的输入参数，表2为输出的选择性输出的关键点信息，可见本发明所建立的仿真方法实现了预期的功能，表3为五种常见空间碎片的计算结果，实现了碎片逃逸功能的预测。图3为是本发明的地基激光移除空间碎片过程的三维建模与仿真方法中计算结果显示输出示意图。

表1地基激光系统和空间碎片相关参数

表2直径1cm碎片的清除过程(地基)

表3

尽管上文对本发明的具体实施方式给予了详细描述和说明，但是应该指明的是，我们可以依据本发明的构想对上述实施方式进行各种等效改变和修改，其所产生的功能作用仍未超出说明书及附图所涵盖的精神时，均应在本发明的保护范围之内。

Claims (4)

1.地基激光移除空间碎片过程的三维建模与仿真方法，包括以下步骤：

（1）建立地心赤道坐标系，以空间碎片的位置矢量(x₁,y₁,z₁)和速度矢量(v_x,v_y,v_z)，以及激光站的位置矢量(x₀,y₀,z₀)作为输入参数，建立M函数rab_check.m，判断碎片是否在激光的作用范围内；其中，基于地心赤道坐标系，采用激光站地表切面，对是否具备清除条件进行判断，假设激光站坐标为 (x₀,y₀,z₀), 碎片的位置矢量为(x₁,y₁,z₁)，地球球面方程为：x²+y²+z²=6378²，根据几何方法确定激光站点在地球表面的切面方程：x₀x+y₀y+z₀z=6378²；只有在碎片移动到激光站上方即x₁x₀+y₁y₀+z₁z₀-6378²>0时，才具有清除条件，此时向量夹角即为天顶角；

（2）建立计算轨道根数的M函数coe_from_sv，将采用地心赤道坐标系中的位置矢量和速度矢量作为输入参数，根据二体运动方程和牛顿定律转换为轨道根数；

（3）建立轨道与时间关系的M函数rv_from_r0v0，对碎片运行轨迹进行计算和信息记录，利用此函数实时跟踪碎片的移动过程，更新碎片的轨道信息；

（4）确定激光参数，基于激光在大气中的传输模型和激光与物质相互作用机理，建立碎片速度增量矢量M函数dv，其中包含已知空间碎片材料的相关参数，对已知空间碎片材料进行计算，确定激光脉冲作用后碎片的速度V_to_dV，对激光作用次数进行计数，计算激光作用后靶的质量m=m₀-μE，μ为激光烧蚀率，E为靶上的激光能量；

（5）激光作用判断，判断激光脉冲作用后碎片轨道近地点是否降低，如果降低则激光作用计数+1，更新碎片速度矢量和轨道根数；反则，碎片按照时间步长继续运动；

（6）每个时间步长判断碎片是否移动出激光作用范围，如果移动出激光作用范围，则计算时间步长增加十倍，以提高计算效率，直至碎片再次进入天顶角90°警戒范围时，时间步长恢复，开始下一个计算周期，并且碎片过顶次数+1；

（7）当碎片近地点高度降至≤200km时，认为碎片坠入大气层烧毁，即成功移除；

（8）在计算结果显示输出方面，模型计算和储存每个时间步长的所有参数信息，用户选择输出所需参数和格式，基于Matlab GUI平台建立了图像显示输出模块，通过该模块可实现碎片变轨过程图形可视化输出，从不同视角采用动画的方式呈现碎片的移除过程，并且对关键位置进行突出化显示。

2.如权利要求1所述的方法，其中，天顶角的警戒范围为45-90°。

3.如权利要求1所述的方法，其中，只有当天顶角在-45°~45°范围内才能使从地面发射的激光辐照在空间碎片上。

4.如权利要求1所述的方法，其中，在三维空间中，存在轨道倾角，激光作用产生的速度增量分量，使碎片轨道倾角增大，逃逸出激光站可管辖的区域，进一步增加了对碎片逃逸情况的判断。