地基激光移除空间碎片过程的三维建模与仿真方法
【专利摘要】本发明公开了一种地基激光移除空间碎片过程的三维建模与仿真方法，基于地心赤道坐标系中碎片的初始位置和速度矢量，以及激光站的位置矢量，激光的波长、频率等参数，计算模拟地基激光移除空间碎片的物理过程。该方法极大的简化了计算方法，便于参数设置和显示输出结果，提高了其准确度和使用灵活性。
【专利说明】
地基激光移除空间碎片过程的三维建模与仿真方法
技术领域
[0001]本发明涉及航天技术领域，尤其是涉及一种地基激光移除空间碎片过程的三维建 模与仿真方法。
【背景技术】
[0002] 空间碎片是指人类空间活动所产生的太空垃圾。目前地球周围的空间碎片总量已 经接近7000t，其中有43 %主要分布在400-2000km的低地轨道上。它们严重威胁着航天器的 安全，减缓空间碎片已势在必行，如何清除空间碎片成为航天领域关注的热点和前沿问题。
[0003] 激光移除空间微小碎片是国际上一致认可的理想碎片移除方法。该方法采用激光 推进的方式，当高强度脉冲激光照射在碎片表面时，使表面材料等离子化产生高温气体，产 生一系列类似于火箭反推的热物质射流，给碎片提供速度增量来降低近地点高度，使其在 较短时间内落入大气层烧毁。
[0004]在激光移除空间碎片过程建模方面，已有学者对其进行了研究。Phipps等建立了 二维变轨模型，通过给定的激光和碎片参数能够预测碎片的降轨高度（参见Claude R.Phipps ,Kevin L. Baker, Stephen B. Libby,et al.Removing orbital debris with lasers.Advances in Space Research,2012,49:1283-1300)。金星等（参见：金星，洪延姬， 常浩.地基激光清除椭圆轨道空间碎片特性的计算分析.航空学报，2013，34:2064-2073)建 立了以速度增量为初始量的单脉冲和多脉冲碎片变轨二维模型，该模型仅能分析二维平面 内空间碎片的降轨过程。韩威华等(参见:韩威华，甘庆波，何洋，杨新.天基激光清理低轨空 间碎片的最佳角度分析与过程设计.航空学报，2015，36:749-756)针对高能脉冲激光清除 低轨空间碎片的问题，建立了高能脉冲激光清除空间碎片作用过程和轨道演化的二维模 型，推导了脉冲激光作用碎片的最佳角度与碎片轨道参数的解析关系。常浩等(参见:常浩， 金星，洪延姬.基于激光烧蚀冲量耦合的空间碎片主动变轨建模与仿真.航空学报，2013， 34:2325-2332)针对空间碎片主动离轨防止碰撞在轨卫星的情况，基于轨道动力学建立了 空间碎片主动变轨物理过程三维模型。该模型的计算方法十分复杂，以碎片轨道根数作为 初始量，碎片进入激光作用范围的判断、碎片轨道参数、速度增量等需要在地心惯性坐标 系、轨道坐标系和地心赤道坐标系之间进行转换，这无疑增加了系统的复杂程度和计算量， 并且忽略了碎片逃逸等实际问题。
[0005]激光移除空间碎片是一个包括多脉冲激光冲量耦合作用、碎片变轨、陨落和逃逸 等的复杂过程。虽然目前已经建立一些激光清除碎片过程的轨道动力学模型，但是它们大 多数为二维模型，并且以速度增量为初始参数，而忽略了激光与物质的冲量耦合作用、碎片 轨道参数变化和轨道倾角变化对速度增量的影响。在已建立的三维模型仅针对某种特定的 变轨情况，在激光移除空间碎片领域具有一定的局限性，忽略了真实空间环境下的一些实 际问题，并且计算方法十分复杂，结果信息输出不够详尽与直观。因此，研究一种具有普适 性的、计算模型简单、能够描述激光移除空间碎片的真实物理过程的三维模型对定量分析 激光清除空间碎片的过程具有重要意义。
【发明内容】
[0006] 鉴于现有技术的上述缺陷，本发明旨在提供一种具有普适性、计算模型简单、能够 描述地基激光移除空间碎片真实物理过程的三维仿真方法，建立包括激光冲量耦合作用、 碎片变轨、陨落和逃逸等关键物理过程的仿真模型，真实、准确的输出激光移除空间碎片的 物理过程，为激光参数选取和碎片移除的可行性判断提供基础数学模型。
[0007] 本发明需要解决的技术问题包括：（1)如何采用简单的计算方法将地基激光移除 空间碎片中的一系列关键问题集成在一个三维模型中；（2)如何避免在不同坐标系之间转 换来判断碎片是否进入激光作用范围；（3)如何方便的设置参数和显示输出结果。
[0008] 本发明采用了如下的解决方案：
[0009] 一种地基激光移除空间碎片过程的三维建模与仿真方法，包括以下步骤：
[0010] (1)建立地心赤道坐标系，以空间碎片和激光站的位置矢量（x，y，z)和速度矢量 (vx，Vy，vz)作为输入参数，建立Μ函数rab_check.m，判断碎片是否在激光的作用范围内；其 中，基于地心赤道坐标系（单坐标系），采用激光站地表切面，对是否具备清除条件进行判 断，假设激光站坐标为(xo，yo，zo)，碎片的位置矢量为( X1，yi，Z1)，地球球面方程为:x2+y2+z 2 = 63782，根据几何方法确定激光站点在地球表面的切面方程:XQX+yoy+Z()Z = 63782;只有在 碎片移动到激光站上方时（xixQ+yiyQ+zizo-esTS2^))，才具有清除条件，此时向量夹角即为 天顶角；
[0011] (2)建立计算轨道根数的Μ函数C〇e_fr〇m_SV，将采用地心赤道坐标系中的位置矢 量和速度矢量作为输入参数，根据二体运动方程和牛顿定律转换为轨道根数；
[00?2] (3)建立轨道与时间关系的Μ函数rv_f rom_r0v0，对碎片运行轨迹进行计算和信息 记录，利用此函数实时跟踪碎片的移动过程，更新碎片的轨道信息；
[0013] (4)确定激光参数，基于激光在大气中的传输模型和激光与物质相互作用机理，建 立碎片速度增量矢量Μ函数dv，其中包含常见空间碎片材料的相关参数，对常见空间碎片材 料进行计算，确定激光脉冲作用后碎片的速度V_to_dV，对激光作用次数进行计数，计算激 光作用后靶的质量πι=π?()-μΕ，μ为激光烧蚀率，E为靶上的激光能量；
[0014] (5)激光作用判断，判断激光脉冲作用后碎片轨道近地点是否降低，如果降低则激 光作用计数+1，更新碎片速度矢量和轨道根数;反则，碎片按照时间步长继续运动；
[0015] (6)每个时间步长判断碎片是否移动出激光作用范围，如果移动出激光作用范围， 则计算时间步长增加十倍，以提高计算效率，直至碎片再次进入天顶角90°警戒范围时，时 间步长恢复，开始下一个计算周期，并且碎片过顶次数+1;
[0016] (7)当碎片近地点高度降至< 200km时，认为碎片坠入大气层烧毁，即成功移除；
[0017] (8)在计算结果显示输出方面，模型计算和储存每个时间步长的所有参数信息，用 户选择输出所需参数和格式，基于Matlab GUI平台建立了图像显示输出模块，通过该模块 可实现碎片变轨过程图形可视化输出，从不同视角采用动画的方式呈现碎片的移除过程， 并且对关键位置进行突出化显示。
[0018] 其中，天顶角的警戒范围为45-90°。
[0019]其中，只有当天顶角在-45°~45°范围内才能使从地面发射的激光辐照在空间碎 片上。
[0020] 其中，在三维空间中，存在轨道倾角，激光作用产生的速度增量分量，使碎片轨道 倾角增大，逃逸出激光站可管辖的区域，进一步增加了对碎片逃逸情况的判断。
[0021] 与现有技术相比，本发明首次提出了一种地基激光移除空间碎片的三维仿真方 法。该方法打破了以往仿真方法的局限性和传统思路，极大的简化了计算方法，便于参数设 置和显示输出结果，提高了其准确度和使用灵活性。基于Matlab建立了物理过程三维模型， 该模型能够根据输入激光和碎片参数，计算碎片的变轨过程，根据需要实时输出碎片的轨 道信息和清除结果，并且实现图形化输出。考虑到速度增量分量对轨道倾角的影响，该模型 增加了碎片逃逸情况的判断，具有更广泛的适用性。
【附图说明】
[0022]图1是本发明的地基激光移除空间碎片过程的三维建模与仿真方法中地心赤道坐 标系和轨道根数示意图；
[0023]图2是本发明的地基激光移除空间碎片过程的三维建模与仿真方法中三维仿真计 算方法流程；
[0024]图3是本发明的地基激光移除空间碎片过程的三维建模与仿真方法中计算结果显 示输出示意图。
【具体实施方式】
[0025]下面结合附图对本发明的地基激光移除空间碎片三维仿真方法进行进一步说明， 该说明仅仅是示例性的，并不旨在限制本发明的保护范围。
[0026]参见图1，图1是本发明的地基激光移除空间碎片过程的三维建模与仿真方法中地 心赤道坐标系和轨道根数示意图，其中，本发明的方法，基于地心赤道坐标系中碎片的初始 位置和速度矢量，以及激光站的位置矢量，激光的波长、频率等参数，计算模拟地基激光移 除空间碎片的物理过程。
[0027] 首先，需要判断碎片是否具有清除条件，即：碎片的位置是否在激光的作用范围 内。如果不在作用范围内，碎片按照脉冲时间间隔t开始移动，直到碎片进入激光作用范围 内，开始计算脉冲激光与碎片的相互作用。依据相互作用机理，得到单脉冲变轨后碎片的状 态矢量，计算变轨后的轨道根数;计算时间t后，碎片具有的位置和速度矢量，开始下一轮计 算;重复计算，直到近地点高度小于设定值（即清除成功），或者碎片超出激光作用范围。对 于后者，根据碎片的移动，重新判断碎片是否具有清除条件，开始二次或多次过顶计算直到 碎片被清除，或者逃逸出具备清除条件的轨道，无法继续清除。
[0028] 2.1参数输入
[0029] 本发明以空间碎片和激光站的位置矢量(X，y，z)和速度矢量(Vx，Vy，Vz)作为输入 参数，调用Μ函数rab_check_new对输入参数的合理性进行判断。首先基于轨道力学对碎片 位置矢量和速度矢量的是否匹配进行判断;再对碎片轨道进行分析，判断是否需要进行移 除操作;然后对激光站位置矢量进行分析，判断激光站是否在地球表面;最后，分析碎片的 初始位置是否在激光作用的区域内，并判断其是否具备移除条件。
[0030] 2.2轨道根数的确定
[0031]在地心赤道坐标系中（参见图2,图2是本发明的地基激光移除空间碎片过程的三 维建模与仿真方法中三维仿真计算方法流程），χ轴指向春分点的方向，XY平面为地球的赤 道平面，Ζ轴与地球的旋转轴一致，且指向北。单位矢量I，J和Κ满足右手定则。在给定时刻， 空间碎片的初始状态向量速度V0和位置Π )，在地心赤道坐标系中状态向量可表示为：
[0032] r〇 = xoI+yoJ+zoK (1)
[0033] vo = Vx〇I+Vy〇J+vz〇K (2)
[0034] 在该点受到激光辐照后，获得速度增量
[0035] Δ ν= Δ νχΙ+ Δ vyJ+ Δ vzK (3)
[0036] 变轨后速度矢量和位置矢量分别为：
[0037] ν = νο+ Δ ν= (νχ〇+ Δ vx)I+(vyo+ Δ vy)J+(vz〇+ Δ νζ)Κ (4)
[0038] r = r〇 (5)
[0039] 根据二体运动方程和牛顿定律可以计算出变轨后的轨道根数:比角动量的模h、轨 道倾角i、升交点赤经Ω、偏心率e、近地点幅角ω和真近点角Θ。
[0040] 2.3轨道与时间的函数
[0041]如果已知to时刻的位置Π )和速度νο,可由拉格朗日系数f和g及其一阶导数，根据下 述表达式求出任意时刻的位置r和速度v:
[0042] r = fr〇+gvo (6)
[0043] v = fr0 + gv0 (7)
[0044] 利用全局变量x和斯达姆夫函数C(z)与S(z)所表示的拉格朗日系数如下：
[0051 ] 对于椭圆轨道，α>〇·
[0052] 如果已知At、r〇、V()和α，可从全局开普勒方程中解出全局近点角X，具体步骤如下：
[0053]计算XQ合理的初始估计值：
[0055]以x〇为初始数据进行迭代计算：
[0062]算出比值= (Xi)。如果 I ratioi I 超出精度范围10-8,则按Xi+i = Xi - ratioi重新选取X估计值，计算式（14) (15)直到| ratioi | <10-8,则接受Xi作为解。这样就可 以根据式(6)(8)(9)得到At时刻的位置矢量r，根据式(7)(10)(11)得到At时刻的速度矢 量v〇
[0063] 2.4判定激光作用条件
[0064]在地基清除系统中，由于受到大气的影响，理论上只有当天顶角在-45°~45°范围 内才能使从地面发射的激光辐照在空间碎片上。但是在计算中，向量夹角的范围在〇~ 180°，需要根据碎片的坐标和象限来计算天顶角，特别是对于任意激光站位置和空间碎片 轨道的情况，这无疑增加了计算的复杂程度。本文采用激光站地表切面的方法，对是否具备 清除条件进行判断。
[0065]假设激光站坐标为(χο,γο,ζο)，碎片的位置矢量为(xi，yi，zi)，地球球面方程为:x2 +y2+z2 = 63782，根据几何方法可确定激光站点在地球表面的切面方程：xOx+yoy+zoz = 63782。
[0066]只有在碎片移动到激光站上方时（XiXo+yiyo+ZiZo-eSTS2〉。），才具有清除条件，此 时向量夹角即为天顶角。这样就可对碎片的天顶角进行精确判断，而无需考虑碎片和激光 站的象限问题。实现方式为：
[0067] elseif(R(l)*RL(l)+R(2)*RL(2)+R(3)*RL(3)-6378~2>0)
[0068] 2.5激光与碎片相互作用
[0069]在激光传输方面，本发明采用简化模型进行分析。为了在较大距离上将激光传输 到靶上，需要数KJ的激光脉冲能量和大型发射镜来克服光的衍射：
[0070] ^3ΕΙΜτ (卜)）
[0071] ds = aM2Az/Deff (20)
[0072] 其中ds为靶上光斑直径;Μ2为光束质量因子U 1，1代表最好）；Deff为计算衍射时 输出孔径D中的有效照射光直径;a为一个乘数因子，对于高斯光束a = 4A。
[0073] 在确定速度增量时，激光与材料的动量耦合系数是主要参数：
[0074] Cm=px/Φ =mA v/E(N · s/J) (21)
[0075] 其中，激光作用后靶的质量πι = ΠΚ)-μΕ;μ为激光烧蚀率;p为一个强度为I、脉宽τ的 激光脉冲在靶上产生的烧蚀压;激光通量为Φ ;Ε为靶上的激光能量。目前已经通过实验获 得了多种材料的U值。在计算模型中，我们对激光辐照次数进行记录，每次辐照后碎片的质 量变化都能够计算，从而获得精确的速度增量。具体实现方式如下：
[0076] Fl = (W*(Deff~2)*4*Teff)/(pi*(]T4M(a*lamda*zr2));% 作用在碎片上的激 光通量％
[0077] E=Fl*(pi*((dt/2r2));%作用在碎片上的能量 ％
[0078] m0 = density*4/3*pi*( (dt/2) ~3); %初始质量％
[0079] m=m〇-(count+l)*u*E; %作用count次后的质量％
[0080] dv = Cm*E/m/1000;% 速度增量％。
[0081 ] 2.6碎片逃逸
[0082]在三维轨道计算中，速度增量可能与碎片轨道不在同一个平面上，这样就会产生 速度增量改变碎片的轨道倾角，当轨道倾角增加到一定程度时，碎片就逃逸出激光的作用 范围，无法继续作用。在地基系统中，如果碎片沿当前轨道移动一周后仍然不具备激光辐照 的条件，则认为碎片逃逸。
[0083] 2.7计算结果显示输出
[0084] 模型记录每个时间步长下的碎片参数，根据用户需要将计算结果以TXT文本格式 输出，方便后期处理。调用plot函数对计算结果进行图形化输出，并对关键点进行标注，具 体实施方式如下：
[0085] 标注碎片初始位置:plot3(rn0(l)，rn0(2)，rn0(3)，'r. '，
[0086] 'MarkerEdgeColor'，'r'，'markerFaceColor'，'r'，'markersize'，20);
[0087] 绘制碎片轨道:plot3(rx，ry，rz,，b，）；hold on
[0088] 选择视角：view([_27 17]);
[0089] 实现动画输出：comet3(rx，ry，rz);
[0090] 表1为计算实例的输入参数，表2为输出的选择性输出的关键点信息，可见本发明 所建立的仿真方法实现了预期的功能，表3为五种常见空间碎片的计算结果，实现了碎片逃 逸功能的预测。图3为是本发明的地基激光移除空间碎片过程的三维建模与仿真方法中计 算结果显示输出示意图。
[0091] 表1地基激光系统和空间碎片相关参数
[0092]
[0093] 表2直径lcm碎片的清除过程(地基）
[0094]
[0095] 表 3
[0096]
[0097] 尽管上文对本发明的【具体实施方式】给予了详细描述和说明，但是应该指明的是， 我们可以依据本发明的构想对上述实施方式进行各种等效改变和修改，其所产生的功能作 用仍未超出说明书及附图所涵盖的精神时，均应在本发明的保护范围之内。
【主权项】
1. 一种地基激光移除空间碎片过程的三维建模与仿真方法，包括以下步骤： (1) 建立地心赤道坐标系，以空间碎片和激光站的位置矢量(X，y，Z)和速度矢量(Vx，Vy， vz)作为输入参数，建立Μ函数rab_check.m，判断碎片是否在激光的作用范围内；其中，基于 地心赤道坐标系(单坐标系），采用激光站地表切面，对是否具备清除条件进行判断，假设激 光站坐标为(XQ，yo，ZQ)，碎片的位置矢量为(XI，yi，Z1)，地球球面方程为:x 2+y2+z2 = 63782， 根据几何方法确定激光站点在地球表面的切面方程:x〇x+y()y+z()z = 63782;只有在碎片移动 到激光站上方时(xix〇+yiy〇+ziz()-63782>0)，才具有清除条件，此时向量夹角即为天顶角； (2) 建立计算轨道根数的Μ函数C〇e_fr〇m_SV，将采用地心赤道坐标系中的位置矢量和 速度矢量作为输入参数，根据二体运动方程和牛顿定律转换为轨道根数； (3) 建立轨道与时间关系的Μ函数rv_from_;rOvO，对碎片运行轨迹进行计算和信息记 录，利用此函数实时跟踪碎片的移动过程，更新碎片的轨道信息； (4) 确定激光参数，基于激光在大气中的传输模型和激光与物质相互作用机理，建立碎 片速度增量矢量Μ函数dv，其中包含常见空间碎片材料的相关参数，对常见空间碎片材料进 行计算，确定激光脉冲作用后碎片的速度V_to_dV，对激光作用次数进行计数，计算激光作 用后靶的质量πι=π?()-μΕ，μ为激光烧蚀率，E为靶上的激光能量； (5) 激光作用判断，判断激光脉冲作用后碎片轨道近地点是否降低，如果降低则激光作 用计数+1，更新碎片速度矢量和轨道根数;反则，碎片按照时间步长继续运动； (6) 每个时间步长判断碎片是否移动出激光作用范围，如果移动出激光作用范围，则计 算时间步长增加十倍，以提高计算效率，直至碎片再次进入天顶角90°警戒范围时，时间步 长恢复，开始下一个计算周期，并且碎片过顶次数+1; (7) 当碎片近地点高度降至< 200km时，认为碎片坠入大气层烧毁，即成功移除； (8) 在计算结果显示输出方面，模型计算和储存每个时间步长的所有参数信息，用户选 择输出所需参数和格式，基于Matlab GUI平台建立了图像显示输出模块，通过该模块可实 现碎片变轨过程图形可视化输出，从不同视角采用动画的方式呈现碎片的移除过程，并且 对关键位置进行突出化显示。2. 如权利要求1所述的方法，其中，天顶角的警戒范围为45-90°。3. 如权利要求1所述的方法，其中，只有当天顶角在-45°~45°范围内才能使从地面发 射的激光辐照在空间碎片上。4. 如权利要求1所述的方法，其中，在三维空间中，存在轨道倾角，激光作用产生的速度 增量分量，使碎片轨道倾角增大，逃逸出激光站可管辖的区域，进一步增加了对碎片逃逸情 况的判断。
【文档编号】G06T17/00GK105868503SQ201610262000
【公开日】2016年8月17日
【申请日】2016年4月25日
【发明人】张品亮, 龚自正, 牟永强, 曹燕, 徐坤博, 杨武霖, 刘国青, 李明
【申请人】北京卫星环境工程研究所