0 引言
火箭发射任务中残骸坠落地面有可能会对人员及设备设施造成损毁和危害,特别是部分残骸还携带有剧毒燃料及火工品,火箭残骸落点预报与回收是一项非常重要工作。随着落区地方经济快速发展,高价值基础设备设施越发增多,落区面临的人员、财产安全及舆情管控形势也愈加严峻。因此,提高火箭残骸落点预报精度,实现残骸精准快速回收是落区工作的紧迫需求。从实际残骸搜索回收结果看,火箭残骸往往与理论落点和预报落点存在较大偏差,对残骸的及时回收和落区现场快速处理造成较大困难。安金霞等结合历史落点数据开展了火箭残骸落点预报精确度分析、火箭残骸预置落点区域准确度分析和火箭残骸落点散布区域特征分析等方法,给出了火箭发射残骸落点“态势一张图”的可视化设计要素[1]。王佩等根据多基站接收到弹箭残骸舱段主动辐射信号的时间差信息,采用Chan算法完成对弹箭残骸舱段位置信息的解算,实现对弹箭残骸轨迹的连续监测[2]。王颍通过蒙特卡洛法对比分析了3自由度和6自由度动力学模型下的整流罩落点散布特性[3]。龚秋武等在质心运动的基础上,对残骸的姿态运动进行了分析与建模,提出火箭子级残骸6自由度落点预报方法,有效缩小了理论落区范围[4]。李瑞鸿等设计分析了栅格舵技术在火箭一子级回收中的运用,对垂直起降重复使用运载火箭的工程推广具有重要支撑作用[5]。张骁等研究了基于分离点信息和蒙特卡诺方法对火箭残骸落点预报的优化方法,取得了较好的改进效果[6-7]。刘海波等基于单台雷达测量数据在算法层对空气阻力系数进行插值同时考虑高空风的影响,在策略层利用预测散布平面的椭圆中心作为预测落点,增强了算法的稳健性[8]。黄普等研究了残骸陨落飞行规律和落点散布规律,建立了一种考虑轨道摄动的高精度陨落弹道计算分析方法,具有一定工程应用价值[9]。王景国等建立了高空风修正的整流罩残骸落点预报模型,使落区搜索面积减少了约90%[10]。徐明兴与贾洪印等分别研究了栅格舵在控制火箭一级残骸落点分布中的运用,研究表明栅格舵机械结构简单,易于工程实现,适用于各类运载火箭的改造,可大幅度减小落区的范围[11-12]。张强等采用基于神经网络特征提取的异源图像匹配定位技术,实现落区目标精确自主定位[13]。陈志会等分析了国内外火箭残骸搜索回收技术,为发展可重复使用运载火箭提供研究发展思路和参考[14]。以上研究针对火箭残骸不同坠落场景,并对常用算法进行了优化改进,同时探索采用新的技术手段实现可控回收,或者尽可能减小残骸落点预报范围,以期提高火箭回收效率,但基于残骸实时跟踪测量数据的高精度落点预报研究目前还未有公开、可靠应用成果。
为此,将引进一台相控阵雷达和机动小光学,对火箭分离残骸的飞行过程进行实时跟踪测量,实现残骸飞行轨迹的精确定位。与之相应,需要建立火箭残骸飞行轨迹的预测模型,为相控阵雷达和机动小光学实时提供实时跟踪引导,便于其搜索捕获和跟踪丢失后的重补,同时基于残骸飞行实测数据的落点计算可大大提升火箭残骸落点的预报精度,实现残骸的精准定位和快速回收。
弹道空间几何模型是残骸无动力飞行弹道计算最为快捷的方法,但计算的坠落飞行弹道与真实情况差异较大,真实环境下残骸坠落的被动段弹道要复杂很多,受空气动力、目标特性等因素的影响。空气动力包括大气风场力和空气阻力,由于火箭在无控坠落过程中速度非常大,因而受空气阻力影响很大,而受大气风场力影响较小[15]。
针对火箭残骸飞行轨迹的特点,文中提出了一种将弹道空间几何模型和空气动力模型组合设计的飞行轨迹预测数学模型,即基于对火箭残骸飞行轨迹的特性分析,将残骸无动力坠落的被动段弹道分解为空间几何段和空气动力段,在空间几何段采用地球扁率修正的椭圆几何模型,在空气动力段采用加入空气阻力的动力学模型。并基于残骸伞降可控回收试验中残骸自定位装置所获取飞行轨迹实时定位实测数据,与预测模型计算的理论弹道进行对比分析,通过优化模型进一步提高模型预测精度。
1 火箭残骸飞行轨迹分析
大气层沿高度分为5层,即对流层、同温层、中间层、电离层和外大气层。对流层在高度11 km以下,大约集中了整个大气质量的75%;同温层高度约为11~25 km,大约集中了整个大气质量的15%;中间层高度约为25~90 km,在中间层的上层80 km以上,可以忽略大气的存在视为真空来处理[11]。
火箭残骸分离高度已处在中间层的上层,空气极度稀薄,空气阻力接近为零。如图1所示,为简化模型,以火箭残骸分离点K为空间几何段的起点,其对称点E为空间几何段的终点,K点与E点之间的残骸飞行弹道,采用几何椭圆模型;空间几何段的终点即是空气动力段的起点,落点C是空气动力段的终点,之间的残骸飞行弹道考虑空气阻力的影响,采用动力学模型。
动力学模型主要考虑空气阻力、地球引力与速度等运动参数关联以及随时间的变化。火箭残骸的飞行速度高度曲线如图2所示,H,V,T分别为飞行高度、速度。当火箭残骸分离后,由于空气非常稀薄,空气阻力接近为零,残骸主要受地球引力作用,飞行轨迹接近几何椭圆,残骸高度增加,速度减小。当残骸飞行过顶点,地球引力的切向分量与速度方向相同,速度增大、高度减小。到弹道点后,随着高度的降低,空气密度逐渐增大,空气阻力也在增大,当空气阻力与地球引力的切向分量相等时,此时的速度达最大值。随着飞行高度的快速下降、空气密度急剧增加,空气阻力也急剧增大、速度急剧减小。之后因速度减小,阻力也减小,同时也因速度倾角的绝对值增大,引力的切向分量增大,加速度又有所回升,但仍为负值,速度仍然不断减少但程度却在持续减缓。
2 飞行轨迹预测模型
2.1 空间几何段轨迹参数计算
其中,rK,VK,θK,σK分别为分离点地心径向距离、绝对速度、速度倾角和方位角。基于残骸飞行轨迹分析需求,残骸分离点的高度HK、星下点经度LK和纬度BK可表示为:
式(2)需要进行迭代计算,根据实际计算,迭代3次以上就可以达到足够高的精度。根据以上残骸分离点的弹道参数,计算空间几何段的椭圆几何模型的轨道参数:
其中,h,VK,P,ε分别为椭圆模型的动量矩、能量比、半通径、偏心率。
$\begin{array}{c}t_{\beta}=\frac{\left(r_{K} v_{K} \cos \theta_{K}\right)^{4}}{\left(1-\varepsilon^{2}\right)^{3 / 2} f_{\mathrm{M}}^{2} h}\left(2 \arctan \left(\sqrt{\frac{1-\varepsilon}{1+\varepsilon}} \tan \frac{f_{K}+\beta_{A}}{2}\right)-\right. \\2 \arctan \left(\sqrt{\frac{1-\varepsilon}{1+\varepsilon}} \tan \frac{f_{K}+\beta_{A}}{2}\right)- \\\left.\varepsilon\left(\frac{\sqrt{1-\varepsilon^{2}} \sin \left(f_{K}+\beta_{A}\right)}{1+\varepsilon \cos \left(f_{K}+\beta_{A}\right)}-\frac{\sqrt{1-\varepsilon^{2}} \sin f_{K}}{1+\varepsilon \cos f_{K}}\right)\right)\end{array}$
其中,ξ,δη是很小的扁率修正参数,可通过复杂公式计算,也可通过实测数据进行反推拟合获得经验值,便于简化模型。顶点至K点的经度之差ΔLKβ,其正弦值可表示为:
sin ΔLKβ=sin βAsin σK/cos φβ
通过以上求出的参数,利用式(2)即可求出轨迹顶点高度Hβ、星下点经度Lβ和纬度Bβ,同时根据椭圆轨道速度计算公式得到Vβ。
以任一点相对于K的射程角Δβ起步进行计算,基于以上计算过程,可以得到任一点的地心径向距离、飞行时间、高度、速度、星下点经纬度等轨迹参数,用于分析椭圆几何模型的精度。
2.2 空气动力段轨迹参数计算
根据残骸飞行轨迹的理论弹道特性分析,考虑空气阻力的影响,火箭飞行速度为V,对当地水平线的倾角为θ,空气阻力为F,空气动力段运动特性如图3所示。
运动微分方程组为[17]:
当火箭以高速进入稠密大气层时,在巨大的空气阻力作用下,残骸将有一个较大的负加速度,速度会急剧下降。火箭所受的空气阻力为:
F=CxSmρV2/2
其中:Cx为空气阻力的阻力系数,可根据试验数据取经验值,一般在0.1~0.2之间;Sm为阻力受力面积,在空间动力段大部分弧段,火箭速度大,Sm趋于最小,在末端弧段,火箭速度小,Sm存在一定不规则变化,但对阻力及坠落飞行轨迹影响很小,因此模型取火箭最小截面积;ρ为大气密度,ρ=ρ0e-αh,α=(29.27T0)-1,T0为海平面的绝对温度,ρ0取理想气体状态下的数值1.225 05。
对于运动微分方程可采用高精度单步积分的四阶“龙格-库塔”法,设时间步长为τ,其基本公式为[7]:
以数据采样点频率为步长进行数值积分计算,取时间步长τ=0.05 s,以为火箭残骸轨迹E点rE,VE,θE为初值xi= ,代入“龙格-库塔”基本公式,计算下一时刻的ti+1=ti+τ的位置与速度。射程角可表示为:
由此即可获得火箭残骸下一时刻的ti+1轨迹参数ri+1,Vi+1,θi+1,βi+1,同时通过式(2)、式(7)获得高度Hi+1、星下点经度Li+1、纬度Bi+1等轨迹参数。
单步计算的结果作为下一步计算的初值,按此策略顺序计算,即可获得火箭残骸在空气动力段飞行轨迹参数数据。
对于残骸落点参数计算,只需计算落区的地心径向距离rC,以此作为“龙格-库塔”计算公式的终止判断:当ri+1<rC时,则计算终止,达到落点位置,积累单步计算时间t2即为空气动力段飞行时间,对应飞行轨迹参数中的经纬度即为落点大地经纬度。
3 模型验证分析
模型验证分析所需的实测数据,是基于典型残骸伞降可控回收试验,通过装载的残骸自定位装置获取的开伞前的残骸飞行轨迹全程自定位数据。
3.1 对比验证
图4 实测与模型计算的残骸飞行轨迹对比曲线Fig.4 Flight trajectory curves of debris measured and calculated by models |
图5 实测与模型计算的残骸飞行高度对比曲线Fig.5 Flight altitude curve of debris measured and calculated by model |
由图4可知,实测与模型计算的飞行轨迹存在一定偏差,但轨迹趋势基本相同,实测和模型计算轨迹均是在下坠到一定高度后,呈垂直坠落趋势。
由图5可知,模型计算的轨迹最高点略高于实测高度,从分离点至残骸下坠到约30 km高度,模型计算的飞行高度轨迹与实测高度基本吻合,在小于这个高度后,模型计算高度与实测高度曲线分离越来越大,即误差不断增大,说明残骸实际坠落的速度远低于模型计算的速度。开伞点是实测数据与模型计算结果进行有效对比分析的最后位置点,高度误差约5.8 km。
由图6可知,自分离点开始,模型计算的星下点轨迹与实测飞行星下点轨迹发生一定偏离并逐渐增大,至开伞点偏差最大。基于球面三角形公式,利用实测与模型计算的星下点经纬度及其偏差数据,计算出开伞点星下点的平面距离偏差约为8.5 km。
3.2 结果分析
综合残骸实测与模型计算的空间飞行轨迹、高度变化和星下点轨迹的对比分析可知,模型计算结果在高度30 km以上与实测结果吻合较好,随后误差逐渐增大,但轨迹趋势基本相同。说明文中提出的空间几何动力组合的火箭残骸飞行轨迹预测模型,设计科学合理,适用性较好。空间几何段椭圆几何模型与实际飞行轨迹吻合度较好,精度较高。空气动力段的动力学模型,高度30 km以下误差逐渐增大,分析主要原因是模型的空气阻力F偏小,以致模型计算的坠落速度要高于实际坠落速度,可根据本次试验数据对空气阻力的阻力系数Cx和大气密度ρ进行近似拟合求出更贴合实际的经验值。星下点轨迹方向与实测轨迹有一定偏离,但偏离不大,可进一步对模型进行修正。
4 模型优化及验证
基于模型的验证分析,可以对模型进行进一步修正优化,并采用另一次伞降可控回收试验的开伞前的实测数据,对修正优化后的模型进行精度验证。
4.1 模型优化修正
对于动力学模型中空气阻力X偏小的情况,基于大气层的分层结构,可通过试验实测数据对阻力系数Cx和大气密度ρ的乘积进行综合分层拟合。基于实测数据通过正交拟合法求得Cxρ在空气动力段大气中间层40 km以上为0.000 267,中间层25~40 km为0.004 35,同温层11~25 km为0.057 2。由于11 km以下的实测数据是残骸开伞状态的数据,与实际自由落体运动轨迹相差很大,可根据经验数据取值为0.215 9。Sm取火箭最小截面积。
再通过式(4)、式(7)求得修正后的C点与分离点的精度之差ΔLKC和地心纬度φC,修正后经度、纬度为:
4.2 精度验证
4.3 结果分析
修正优化后的模型在坠落末端高度较低时,高度误差虽然也是逐渐增大,但误差减少很多,且变化趋势更为接近。星下点轨迹在空气动力段才开始发生一定偏离,但误差也明显减少,在开伞点高度误差最大,高度偏差约为3.1 km,水平距离偏差约为4.0 km左右,空间距离偏差为5.1 km,经计算平均误差约为1.7 km。说明引入的摄动影响修正和采用实测数据拟合求取空气阻力参数,适用有效,计算结果与实测结果吻合度更高,轨迹更为接近,能够有效提高模型预测精度。因为开伞点的高度已处于残骸近乎垂直坠落运动段的中后段,开伞点空间距离偏差可以近似为模型计算的最大误差。
5 结论
文中建立的空间几何模型和空气动力模型组合设计的火箭残骸飞行轨迹预测模型,科学合理,适用有效,在恰当的空间阻力参数和摄动影响的修正优化后,预测精度更加理想。所预测的飞行轨迹能够很好地为火箭残骸跟踪测量设备提供准确实时数据引导,并可提升火箭残骸落点的预报精度,为火箭残骸飞行过程的精准跟踪以及残骸落点的精准定位和快速回收提供理想的预测数据,可大大提升落区应急处置效率,减轻落区工作负担。