张 骁,刘丙杰,王瑞臣
(海军潜艇学院，山东 青岛 266199)
在火箭或弹道导弹的发射过程中，其各级火箭发动机在燃料耗尽或者到达预定位置后会与本体分离，在自身分离速度、地球引力和空气阻力等多种复杂因素的共同作用下，经过一定时间后以高速冲击地面。随着航天和导弹技术的发展应用，产生的发动机残骸越来越多，对地面人员和设施带来的威胁也越来越大。当前各火箭发射基地在选址过程中均对发动机残骸的落点位置进行了考虑，一般能够避开人员密集区域，但仍然需要找到可靠方法计算临时选址或者机动式发射前提下的火箭发动机残骸落点范围[1]。
在以往的发射任务中，指挥人员可以使用简易模型与经验公式对残骸落点进行预测。但随着实际任务数量的增大、火箭射程的增加、目标点位的变化、发射机动范围的扩大、发射条件的放宽，残骸可能的落点范围越发难以预测，给任务负责人员带来了较大的决策压力。因此为进一步减小对地面人员和设施的威胁，避免发生严事故，减小后续残骸回收的难度，需探索建立一种更加精准、可靠、方便的残骸落点计算方法[2]。
目前已有相关文献阐述残骸落点计算的方法，但算法侧重点和面向的任务特点各不相同。文献[3]提出了基于实时定位的火箭残骸落点计算模型，精度较高，但需要残骸具有自主定位和信号发射的能力，无法应用于当前大多数火箭发射任务；
文献[4]提出了非定常 CFD 数值方法，主要针对具有复杂形态的大型残骸，特点是能够实现了对目标气动发热和分离解体过程的模拟，但在落点计算的结果上未能明显降低误差；
文献[5]提出了容积卡尔曼滤波(CKF)和时变自回归(TVAR)模型融合的助推器落点预测方法，在误差上比单一时间序列的预测方法要小，但精度仍不如数值积分方法。
本研究基于发射任务单位的实际信息获取能力和精度要求，以能够获得残骸与火箭分离点在发射坐标系中的具体坐标、各速度分量、发射点的经、纬度、高程、射向、残骸质量(假定不变)和有效参考面积为前提，构建一种新的残骸落点计算方法，用以解决任务单位的实际问题。
火箭发动机残骸在与火箭主体分离后不再受主动力的影响，将完全处于被动段飞行，除地球引力、柯氏惯性力和牵引惯性力等必须考虑的因素外，空气阻力在其被动段飞行的后半部分具有相当大的影响，是残骸落地速度和位置的重要因素。在实际过程中，各力对残骸的具体作用与残骸的初始状态都是十分复杂的，为兼顾模型精确度和简洁实用性，首先对残骸受力过程进行适当的简化[6-7]。
1.1 残骸弹道特性分析
1.1.1 升降两端等高度残骸速度模型
如图1所示，假设残骸由点q运动至相同高度的对应点q′的情况下，根据理论力学相关知识可列公式：
(1)
式中,F为从q点到q′的所有外力，即引力mg和空气阻力X的和，通过对作用距离s进行积分得到外力F对残骸所做总功，其值应当等于从q到q′的动能变化量。由于q和q′位于等高度，即两点地心矢径的数值相同，因此引力势能相同，即在此过程中引力做功总量为零。而空气阻力X的方向一定与速度方向相反，对残骸做功只能是负值。因此引力mg和空气阻力X在这一过程中对残骸的总做功为负值，即公式(1)右边的积分一定为负值，由此可推出Vg′ 图1 残骸受力分析图 1.1.2 残骸速度极值点位置分析 在不考虑空气阻力的情况下，只有引力对残骸做功，而弹道顶点为引力势能最大的点，因此动能一定最小，即速度最小。但考虑空气阻力后，残骸速度变化就复杂很多，不仅最小值点将很难确定，而且会出现包括速度极大值点在内的多个极值点，但无论如何，极值点的切向加速度等于零为必要条件，即dv/dt=0。 残骸在其运动的切线方向只受两个力的影响，空气阻力X和引力mg在切线方向的分量mgsinθ(θ为速度与水平方向的夹角)。切向加速度为零则意味着X与mgsinθ大小相等方向相反。根据图1可知，在点q时的X和mgsinθ的方向均与速度方向相反，残骸一直处于减速状态，只有在到达点q′时X和mgsinθ才可能出现方向相反的情况，因此速度极值点必然仅存在于下降段[8]。 一般情况下存在两个残骸速度极值点：点一是速度极小值点，位于残骸刚经过弹道顶点的时候，此时由于空气极为稀薄，空气阻力X的值与引力mg相比很小，因此只需要很小的sinθ值即可达到切向加速度为零的条件，即刚过弹道顶点时速度向下且与水平方向成某一小角度时。点二是速度极大值点，位于残骸以极高速度进入稠密大气层的时候，此时的空气阻力X受残骸速度和大气密度的双重影响而极速增大，直至在某一位置与mgsinθ达到平衡，即此时的残骸虽然处于弹道的下降段，但停止加速，转为减速。 上述两个速度极值点一定存在，但在实际中很可能不仅仅有两个极值点。残骸在进入稠密大气层后的减速阶段后期，sinθ趋近于-1，即弹道接近垂直，在残骸速度的减小的同时，大气密度却不再有明显变化，空气阻力X逐渐减小后与引力mgsinθ缓慢趋向于平衡，加速度趋近于0，可视为第三个速度极值点。如果考虑实际更复杂的情况，火箭发动机残骸一般具有较为复杂且变化的气动外形，其等效参考面积较大，更容易受到空气阻力的影响，甚至在下降段可能存在空气抬升力和下降力的复杂影响，因此有更多极值点出现的可能性。 1.1.3 残骸射程影响因素分析 将式(1)中做功为零的mg去掉并适当变形可得： (2) 残骸的坠落过程中，其有效面积S和质量m并不是一成不变的，在超高速飞行时，空气对残骸产生的热效应和剥离效应将直接影响S与m的值。对这一过程进行量化的最佳方法是建立S、m与残骸速度V、残骸所处大气状态C的函数关系，且必须考虑残骸材料和形状的影响。具体是否复杂化，需综合考量模型使用单位对精度和效率的要求。 了解上述残骸弹道特性后有助于我们在建立落点预测模型的过程中的整体把握，减少明显理论错误发生的可能性，且有助于后期定性的验证模型建立的准确性。 为充分说明典型火箭发动机残骸的坠落过程，根据某次试射任务的观测数据，制作图2～4所示。 图2 V-t图 图3 H-t图 图4 sinθ-t图 上述主要运动参数的分析根据分离点的高低、分离速度的大小和方向改变而稍有变化，但仍可在后续模型建立、修改和验证的过程中起到有效指导作用。 残骸再入大气后的受力情况十分复杂，不仅受地球引力作用，还受空气阻力和空气动力矩的影响，在气动加热的情况下，还会使残骸质量、外形发生较大变化，这些都会在一定程度上影响残骸落点。为使问题简化，此处不考虑气动加热和残骸外形的影响。残骸受分离力的影响，还存在一定程度的翻滚，在10 km高度以上时由于空气阻力微弱，其翻滚造成的影响可以忽略不计，在进入稠密大气后残骸的翻滚会很快受到阻力限制而趋向以某有效面积与速度方向保持稳定角度。因此认为全程忽略空气动力矩与翻滚的影响，将残骸视为一个质点仍能够保证模型的有效精确。 残骸无动力段的运动形式与火箭被动段的运动相似，在模型建立过程中可充分借鉴远程火箭弹道学中关于被动段的成熟理论。在将残骸看成一个质点的情况下，残骸主要受力为地心引力与空气阻力。火箭与残骸分离点较高，飞行时间较长，不能忽视因地球转动所产生的牵连惯性力和柯氏惯性力的影响。根据动力学基本理论，在发射坐标系中可得到火箭残骸无动力段的矢量表达式为: (3) 为方便计算，在对弹道进行研究的时候一般使用发射坐标系，以发射点为原点，以射面与水平面的交线为x轴，以过原点并与发射点位置的水平面相垂直的轴为y轴，z轴与x轴和y轴成右手螺旋关系。式(3)在发射坐标系中的投影形式为: (4) 又根据运动学可得方程: (5) (6) 在抬升力、侧向力和空气动力矩忽略不计后，空气给残骸带来的影响只剩空气阻力X。该阻力方向恒定与速度方向相反，在速度坐标系中可表示为： X=CxdtqdtSmdt (7) 式中,Cxdt为残骸阻力系数，qdt=1/2ρV2为头部速度头,ρ为空气密度，Smdt为残骸的有效面积(根据经验可取最大横截面积)。 即使通过上述简化后，想要对X的值进行精确计算仍然是一个非常复杂的问题，难以仅通过理论推导和计算来完成。目前比较科学可靠的方法是使用实验与模型相结合的方法，即通过空气动力学理论尽可能详细的建立接近现实变化的数学模型，再通过空气动力学实验的方法对模型中的关键系数与其他要素的变化规律进行观察记录。即使如此也仅能获得当前实验条件下的数据，环境稍作变化就可能出现不同的情况。开发人员在火箭的研制过程中必然会根据其外形和材料进行空气动力学研究，但由此获得的图表和数据只能在一定程度上代表火箭本身，无益于此处对发动机残骸空气动力系数的研究。 Cxdt在残骸运动的过程中变化最为复杂，在以往案例中通常使用经验数值来表示。其影响因素很多，主要包括残骸的外形和材料、残骸在空中飞行的姿态、残骸速度、飞行高度等。在掌握上述影响因素与Cxdt的关系后才可能使用数学表达式对Cxdt进行近似表示。 残骸的外形和材料对空气动力的影响很大，具体作用过程也最为复杂。简单来说，长细比越大的物体，其阻力系数越小，表面材料越光滑且凸起越少的物体，其阻力系数越小。残骸速度对Cxdt的影响主要通过改变残骸周边气体的可压缩程度来实现，但要找到Cxdt与速度V之间的直接关系是极为复杂的，此处速度需要用马赫数ma来表示，音速u与高度H之间本就存在复杂关系，且ma与气体可压缩程度的关系也并不是简单线性，而是根据研究对象的形状和材料不同表现出不同的曲线，因此在实际操作中往往需通过查阅实验数据表格获得。同时Cxdt与大气的各项基本属性例如密度、压力、温度等均密切相关，而大气中这些属性的最主要决定者就是高度H，因此Cxdt与高度H之间存在着复杂关系。综上可得，在研究对象确定的情况下，对空气动力系数Cxdt影响最大的因素便是飞行马赫数ma和高度H。 同理，ρ受其他因素影响较小，可认为仅与高度H相关。 鉴于Cxdt的详细估算十分复杂，因此在以往计算中常使用经验值代替，也有文献显示落点计算仍可达到一定精度。但如此则将其过分简化，必然无法反映出飞行物或残骸在大气中的真实运行轨迹[3,6,11]，尤其是在残骸对大气层的入射角度较小时，有阻力段在全射程中的占比增加，其误差必然增大。 本研究根据国际标准大气表(1976美国)的参数对空气密度ρ和不同高度的音速u与几何高度H的关系进行了探索，使用matlab2013b进行函数拟合。因空气密度ρ的数值从H=0到H=80 km的范围内跨越6个数量级，导致用一个连续函数进行表达时必然会出现误差平方和极低，但中高空拟合结果与实际数值出现1至2个数量级差距的问题。推测此处问题的本质原因可能在于地球大气的复杂分层导致无法使用一个连续函数进行合理表示。因此建议采用分段拟合的方式，可得到误差平方和在10-5级别的表达式ρ=f(H)和u=h(H)。根据某火箭发动机残骸Cxdt系数与马赫数ma的关系表格： 表1 Cxdt与ma的关系 对Cxdt与ma的函数关系进行拟合，得到Cxdt=g(ma)。联立三式可得Cxdt与H、V的函数关系式: Cxdt=I(H,V) (8) 将式(8)代入式(7)即可得到在速度坐标系中的空气阻力X，需将其转换至发射坐标系才能够代入式(4)进行计算。根据速度坐标系与发射坐标系的转换关系可得: (9) 其中:Xi(i=x,y,z)为空气阻力X在发射坐标系各轴上的分量，θ是指残骸速度矢量V与发射坐标系oxz平面即发射点所在水平面所成的角度，上为正，下为负，σ是指残骸速度矢量V与发射坐标系oxy平面即射面所成的角度，从坐标原点沿发射坐标系x轴正方向看时，右为正，左为负。由定义可见，角θ变化范围较大，角σ在主动段由于控制系统的作用，数值很小，在被动段其对射面的偏离主要来自分离时z方向的初始速度和横向高空风的影响，在以往建模中一般忽略不计。但正如1.1.3中所作分析，由于残骸的S与m比值相对火箭本身较大，更容易受偏离力的影响，所以此处不可忽略不计。由θ和σ的定义可推出: (10) 在将地球看做正常地球椭球体时，可通过下列公式计算得到地球引力在发射坐标系各轴上的分量[12]： (11) (12) (13) (14) B-φs=690.309″/3600/180*π*sin2B (15) (16) (17) (18) 式中,ω、fM、μ均为接近地球实际的总地球椭球体基本数据，可通过查表获得，r为残骸到地心的矢径，gr为引力加速度在地心坐标系矢径r上的分量，gω为引力加速度在地球自转角速度矢ω上的分量。R0i(i=x,y,z)为发射点地心矢径R0在发射坐标系各轴上的投影。B为残骸所处地理纬度，φs为残骸所处的地心纬度，B0和φs0分别指发射点的地理纬度和地心纬度，北纬定义为正，南纬定义为负，C=B0-φs0为发射坐标系中oy轴与矢径R0的夹角，即发射点地理纬度与发射点地心纬度的差值,可通过常用经验式(15)计算得。ωi(i=x,y,z)为地球自转角速度矢ω在发射坐标系各坐标轴上的分量，B0、A0分别为发射点地理纬度和发射方位角[7-8]。 φs的计算方法除使用上述式(14)外，还可采用由发射坐标系向地心直角坐标系转换，再向地心极坐标系转换的方法。在需要得到地理纬度B时可采用式(15)进行计算。 根据定义可知，发射坐标系是固连在发射原点上的，因此在随地球相对惯性坐标系做平移运动的前提下，还绕地球两极做圆周运动。根据理论力学知识，导弹残骸在运动过程中始终受柯氏惯性力的影响，其在发射坐标轴上的加速度分量可分别表示为[7]： (19) 因发射坐标系随地球相对惯性参考系运动，因此必须考虑施加于导弹残骸上的牵引惯性力，其在发射坐标轴上的加速度分量可分别表示为[7]： (20) 其中: (21) 联立上述(4)至(21)式，在给定初始x、y、z和Vx、Vy、Vz的情况下通过对时间积分即可得到残骸弹道数据。该结果数据为残骸在发射坐标系中的坐标和速度，要得到常用地理经、纬度坐标，还需要进行由发射坐标系至地心直角坐标系，再到地心极坐标系的转换，然后使用式(15)得到地理经、纬度。 本研究使用Matlab2013b进行仿真实验，根据上述公式进行建模并以固定步长进行积分，直至满足条件高度H小于等于设定高度H0时停止积分，输出落点经纬度及其它运动参数。 为评价模型精度，首先可采用实际落点(LMDs，Bs)与预测落点(LMDc，Bc)之间距离Δs来表示误差的大小。因两点距离较近，可忽略地球扁率影响得[13]： cos(Bs)cos(Bc)cos(LMDs-LMDc)] (22) 此处Δs的值即为预测值与实际值的绝对误差，虽然能够在一定程度上说明模型精度，但仅使用Δs无法全面比较不同落点计算模型之间的精度优劣。不同模型所使用的数据均来自不同的任务或实验，发射点和目标点各异，残骸飞行距离差别巨大，甚至达到上千公里。当残骸射程相对增大时，落点计算模型所积累的误差也随之增大，Δs的绝对值也必然相应增大。因此为将不同射程的实测数据用来对比，模型精度的评价参数除Δs外还应当包括其与残骸射程的比值Δs/L[14-15]。 使用某火箭位一固定地点向固定方向进行多次实射所采集的数据进行验证，取步长dt=0.01 s，落点高度H0=0结果如表2和图5所示。 根据对比计算的结果显示，在残骸射程由235.928 km增大至1 469.388 km的五次实验中，仿真结果的绝对误差分别为192.373 m、1 081.079 m、3 274.104 m、5 508.066 m、4 292.573 m，与3.1节中的推论相应和。随着残骸射程的增加绝对误差Δs有明显的增大，在射程达到1 469.388 km时，误差最大达到5 508.066 m，误差百分比达到了0.389 9%。综合来看，五次仿真的平均相对精度a(Δs/L)=0.220 6%。 考虑到模型建立的主要目的，千分之二左右的误差能够满足预定的落点预测功能，足以实现安全预警和便于回收等一系列目的，能够为发射指挥人员的判断提供可靠依据。 论证该落点计算模型的精度，不能仅依据现有数据，还应当能够在与同类型预测模型的对比中保持一定的优势地位。目前所阅相关文献中均对所建模型的精度进行了一定的论证，但达到足够精度并能够给出足够数据用以支撑模型之间横向对比的文献却并不多见。决定采用文献[3]和文献[4]中的模型与本落点计算模型进行对比。 表2 精度分析结果 图5 3号仿真轨迹图示 文献[3]中使用火箭残骸实时定位信息为基础，在残骸落地前的40 km至20 km处测量得到残骸具体运动参数，用以计算残骸在高度10 km处的位置。型号一火箭残骸在高度相差10 km时平均偏差为0.63 km，高度相差20 km时平均偏差为0.72 km，高度相差30 km时平均偏差为1.36 km。型号二火箭残骸在高度相差10 km时平均偏差为0.51 km，高度相差20 km时平均偏差为0.74 km，高度相差30 km时平均偏差为1.27 km。综合来看该模型能够将绝对误差控制在1.36 km至0.51 km之间，但相对误差却达到4.533%与5.1%之间。该方法虽能够控制绝对精度，但其只对残骸落地前的一小部分有空气阻力弹道进行预测，且需要残骸的实时定位信息为基础数据，要求在残骸上内置定位装置，对观测能力或者信息传输能力要求较高，适用性受到一定限制。 文献[4]通过建立火箭不同仰角和侧滑角下的气动数据库，对初始仰角、侧滑角与射程、侧位移的关系进行了详细研究，使用非定常CFD方法能够较为精准的预测火箭再入大气层后的形变、解体和落点预测。进行了火箭发动机残骸自30 km高度陨落的飞行试验，绘图结果显示在射程约52 km的情况下误差在1 km以内，但因并未给出具体数据，所以无法进一步与本文方法比较。该方法预测的绝对误差较小，但需要占用的计算资源较多，且不具有普适性，主要针对结构较为复杂或形态巨大的特定目标。 本研究所用方法从分离点的运动参数入手，能够对发动机残骸运动的全过程进行预测，并且不依靠后续观测数据的矫正，占用资源较少，计算速度较快，绝对误差控制在任务单位的搜索能力范围内，相对误差较小，能够适用于结构和材料较为简单的所有发动机残骸落点预测，符合执行任务单位的能力要求。 本研究缺少残骸的真实形状、材料、剩余燃料等重要数据，因此对空气动力的模拟只局限于速度坐标系下沿x轴负方向的阻力，忽略了空气升力和侧向力及其力矩带来的旋转，不能完全反映空气对残骸的影响，同时也无法估算出分离点的实际残骸质量和后续因气动发热造成的质量损失[16-17]。限于信息获取能力有限，代入的初始状态为火箭本体与发动机分离时的运动数据，忽略了分离过程对残骸运动状态的影响。也未能充分考虑高空风和扰动引力在残骸近地飞行时对轨迹带来的扰动[18-19]。后续还可使用神经网络等智能预测方法对同课题进行研究，尝试在残骸基本属性未知的情况下对落点进行预测[20]。 根据某火箭实际使用和训练任务的需求推动，在理论力学的基础上建立动力学方程组，将地球大气密度、空气阻力系数和音速a与高度H的关系进行拟合，代入方程组得到满足精度要求的残骸落点计算模型，以分离点数据为残骸初始状态，对时间积分即可得到残骸全射程运动参数，进而实现对落点的准确预测。 通过与任务实际数据进行对比，证明该方法误差在合理范围内，能够满足指挥员的辅助决策需求。通过与其他方法的横向对比，证明该方法在当前任务需求下具有计算速度快、信息需求少及经济性好的优势。后续仍可对残骸受力的过程进一步细化，逐步加入分离力、气动升力和侧向力、扰动引力及质量损失、形态变化和分离解体等因素的影响。
1.2 残骸动态分析
2.1 残骸所受空气阻力X的计算
2.2 残骸所受地球引力加速度g的计算
2.3 残骸所受柯氏加速度的计算
2.4 残骸所受牵连加速度的计算
3.1 精度分析方法
3.2 精度分析结果
3.3 横向精度对比
3.4 模型问题分析