基于不确定性的飞行器分离可靠性建模与分析方法

张海瑞，王 浩，王 尧，洪东跑

(中国运载火箭技术研究院，北京 100076)

0 引 言

飞行器分离是飞行器发射任务的关键环节之一，其功能是在一定的高度将完成预定功能的部分结构分离并抛去，旨在改善总体性能，增大飞行速度，提高运载能力[1-3]。飞行器实施分离时的恶劣环境将导致两体可靠分离困难，严重制约飞行器可靠性水平的提升。在分离过程中存在着诸多不确定性因素[4]，如反推发动机作用时间不同步，下面级残余推力偏斜以及残余推力作用点位置偏差等。在这些不确定性因素的综合影响下，两体分离存在碰撞风险，可能导致分离任务失败。

传统的分离分析方法通常是基于参数偏差进行的，通过偏差组合考虑最严酷的情况[5]，以此判断分离设计方案是否满足要求。然而随着飞行器可靠性要求不断提高以及新技术方法的发展应用，传统方法为了包络最严酷的极限条件，其设计方案往往偏于保守，增加结构重量，在某些情况下甚至严重影响飞行器总体性能，难以满足精细化设计需求。

针对这一问题，目前的研究趋势是综合考虑不确定性因素影响，将不确定性注入分离动力学仿真模型，建立飞行器分离可靠性模型，实现分离方案精细化分析[6-8]。然而，不确定性因素和模型非线性的影响导致分离可靠性求解困难。国内外学者系统研究了利用蒙特卡洛方法分析不确定性因素对分离运动的影响[9-10]，Roshanian等[11]采用蒙特卡洛方法研究了不确定性因素对分离运动的影响，给出了分离运动参数的包络曲线。沙建科等[12]针对某型导弹级间冷分离方案，采用蒙特卡洛方法对随机不确定性进行了分离仿真分析，给出了两体不发生碰撞的分离可靠性。蒙特卡洛方法具有无偏性以及非侵入性，在工程中有着广泛的应用。当前，对如何进行分离可靠性高效分析与精确评估的研究较少，特别是随着飞行器分离可靠性指标的不断提高，蒙特卡洛方法调用分离动力学仿真模型的次数也会随之急剧增加，从而消耗巨大的计算资源，甚至在有限的研制时间内无法给出精确的分离可靠性分析结果。

为进一步提高可靠性分析的效率，国内外学者提出了诸多高效的可靠性分析方法，其中基于代理模型的可靠性分析方法具有简单直接，适用广泛的特点，是可靠性分析研究的发展趋势之一[7]。文献[13]基于Kriging模型的特点提出了高效全局优化(Efficient global optimization, EGO)方法用于求解复杂函数的优化问题。文献[14]基于EGO方法提出了高效全局可靠性分析(Efficient global reliability analysis, EGRA)方法求解复杂问题的可靠度，其中，采用期望可行函数(Expected feasibility function,EFF)序贯增加训练点。文献[15]结合Kriging模型的特点通过U函数序贯增加训练点，提出了主动学习Kriging方法求解隐式功能函数的可靠度。文献[16]提出了风险期望函数(Expected risk function，ERF)序贯加点，求解非线性功能函数的混合可靠度。文献[17]针对小失效概率问题，将主动学习Kriging和子集模拟方法相结合，提出了AK-SS方法求解隐式功能函数的小失效概率问题。目前，可靠性分析方法已经取得了长足的进步，但是面向高超声速飞行器分离任务需求，开展高效分离可靠性分析的研究还较少。

本文以某轴对称式飞行器级间冷分离方案为研究对象，建立了分离动力学仿真模型，为满足精细化设计需求，通过综合考虑分离过程不确定性因素的影响，构建了飞行器分离可靠性模型。结合代理模型和蒙特卡洛方法的特点，提出了改进主动学习Kriging(Improved active kriging, IAK)的分离可靠性分析方法，给出了一种基于不确定性的飞行器分离可靠性建模与分析方法。

1 飞行器分离动力学建模

发射过程中飞行器的正常分离通常包括助推器分离、级间分离、整流罩分离、星箭分离等[1]，在飞行器设计过程中，级间分离技术仍然是需要攻克的一项关键技术，包括级间热分离和级间冷分离两种方法[18]。本文选取典型的轴对称式飞行器级间冷分离方案为研究对象。针对其环境特点与性能要求，建立分离动力学模型。分离过程中，上面级受到重力、分插拔脱力和气动力的影响，下面级受到重力、分插拔脱力、主发动机残余推力、反推发动机推力以及气动力的影响，飞行器级间冷分离示意图如图1所示。由于上面级发动机位于下面级的壳体内，故在分离过程中需要着重关注两者之间的相对距离。

图1 飞行器级间分离示意图Fig.1 Illustration of vehicle stage separation

定义飞行器上面级和下面级的弹体坐标系原点与各自质心重合，x轴与各自弹体纵轴重合；定义分离坐标系与分离初始时刻飞行器组合体的弹体坐标系重合，由于忽略地球自转影响，分离坐标系即为分离惯性参考系。以此为基础，在分离坐标系下分别针对上面级和下面级建立质心运动方程如下

(1)

其中，m是分析对象的质量；vx,vy,vz分别是分析对象的速度矢量在分离坐标系下的速度投影；Fx,Fy,Fz分别是分析对象受到的合力在分离坐标系下的投影。进一步给出分析对象弹体坐标系到分离坐标系的转换矩阵为

GB=

(2)

其中，φ,ψ,γ分别为分析对象的弹体坐标系相对分离坐标系的欧拉角，通过上述转换矩阵将分析对象在弹体坐标系下受到的合力转换为分离坐标系下受到的合力。

进一步，分别在上面级和下面级的弹体坐标系下建立绕质心的运动方程为

(3)

其中，ωx1,ωy1,ωz1为分析对象转动角速度在弹体坐标系中的分量；Ix1,Iy1,Iz1为分析对象相对其弹体坐标系的转动惯量。分别以上面级和下面级为分析对象，补充建立速度与位置、欧拉角与角速度之间的关系，如式(4)所示:

(4)

针对分离过程的具体特点，结合上述方程分别建立组合体、上面级和下面级的六自由度刚体运动及动力学仿真模型。考虑到下面级主发动机关机后仍会产生残余推力，且残余推力随时间变化而逐渐减小，进一步考虑残余推力的偏移偏斜影响，如图2所示。

P是残余推力的等效作用点，X10,Y10,Z10与下面级弹体坐标系方向平行。其中，λ角为残余推力的推力线与下面级弹体中轴线的夹角；θ角为残余推力的推力线在截面Y10PZ10上的投影与Y10轴的夹角，给出残余推力在弹体坐标系下的分量为

(5)

图2 残余推力偏移偏斜示意图Fig.2 Illustration of the offset and deviation of residual thrust

由于偏移偏斜的影响，残余推力的作用点与下面级质心的相对位置为rpc=[xpc,ypc,zpc]T，发动机推力偏移偏斜产生的力矩为

(6)

在飞行器分离过程中，存在多种分离故障失效模式，其中较为典型的故障判据为两体在分离时发生干涉碰撞，这里采用文献[10]中的碰撞检测方案计算上面级喷管底端面中心相对下面级中轴线的最大偏移距离d，根据分离结构的具体特点，给定最大偏移距离的阈值d0，若从开始分离到喷管完全脱离底端面的过程中，最大偏移距离小于给定阈值d0，则判定分离过程中未发生碰撞。中心点偏移示意图如图3所示。

图3 中心点偏移示意图Fig.3 Illustration of the center offset

(7)

2 飞行器分离可靠性建模

结合分离动力学模型，针对轴对称式飞行器级间冷分离方案的特点，综合考虑质量、气动、残余推力以及反推发动机等不确定性因素影响，给出分离过程的不确定性因素如表1所示。

表1 不确定性因素名称及符号Table 1 Names and symbols of uncertainty parameters

通过中心差分方法对上述20个分离不确定性参数进行显著性分析[19-20]，经归一化处理后，给出上面级喷管底端面中心相对下面级中轴线最大偏移距离d影响最大的10个因素如图4所示。

图4 不确定性因素对最大偏移距离影响百分比排列图Fig.4 Effects of uncertainty parameters on maximum offset distance in the separation process

基于上述分析，选取主要不确定性因素包括：反推发动机不同步时间t0(x16)，气动特性偏差Wqd(x20)，反推推力作用点位置xFt(x14)，残余推力的推力线与下面级弹体中轴线的夹角λ(x18)以及残余推力作用点位置xFe(x17)。这5个主要不确定性因素的影响百分比总和超过95%，识别主要不确定性因素，滤除影响较小的次要不确定性因素。

为满足精细化设计需求，在建立分离不确定性模型的基础上，将主要不确定性因素注入分离动力学模型，建立飞行器分离可靠性模型。结合分离碰撞判据，给出可靠性模型的极限状态方程为分离过程上面级喷管底端面的最大偏移距离d与给定阈值d0之差等于零，即：

g(t0,xFt,xFe,λ,Wqd)=d(t0,xFt,xFe,λ,Wqd)-

d0=0

(8)

其中，t0,xFt,xFe,λ,Wqd为不确定性参数，失效域定义为

D={(t0,xFt,xFe,λ,Wqd)g(t0,xFt,xFe,λ,Wqd)≥0}

(9)

因而在分离过程中两体不发生碰撞的分离可靠性为

R=Pr(d(t0,xFt,xFe,λ,Wqd)-d0<0)

(10)

进一步定义可靠性模型极限状态方程的联合概率密度函数为p(t0,xFt,xFe,λ,Wqd)，分离可靠性还可以表示为

(11)

3 飞行器分离可靠性分析

由于不确定性参数为随机变量，通过试验数据和工程经验确定分离不确定参数的分布，进而由式(11)可得分离可靠度。鉴于不确定参数较多，难以获得失效域D与联合概率密度函数p(t0,xFt,xFe,λ,Wqd)的显示表达，直接积分求解难度较大且复杂耗时，因而通常利用蒙特卡洛法进行可靠性分析，此方法精度较高，但随着分离可靠性指标的提高，调用分离仿真模型的次数也会随之急剧增加，导致求解效率较低，难以满足快速设计迭代的工程需求。针对这一问题，文献[15,21]提出了主动学习Kriging的可靠性分析方法，具有很高的效率，在蒙特卡洛样本点中直接选取失效概率最大的样本点作为新增训练点，然而其新增训练点的选择完全取决于蒙特卡洛样本点集中失效概率最大的样本点，无法选取样本点集之外失效概率更大的未知点，从而限制了效率和精度的进一步提升。

结合代理模型和蒙特卡洛方法的特点，提出了一种改进主动学习Kriging的分离可靠性分析方法，通过新的采样策略在失效概率最大样本点附近完成采样，在满足采样点与训练点集相关条件的基础上，给出失效概率更大的新增训练点，使得新增训练点的最大偏移距离更加趋于给定阈值。

记分离不确定性参数向量x=(t0,xFt,xFe,λ,Wqd)，利用Kriging代理模型，将真实最大偏移距离d与给定阈值d0之差表达为如下形式

(12)

其中，f(x)=[f1(x),f2(x),…,fn(x)]T为回归基函数，针对分离动力学模型的具体特点，采用二阶多项式回归模型；β=[β1,β2,…,βn]T为回归系数，z(x)为高斯随机修正过程。在给出未知点x预测最大偏移距离与给定阈值之差的基础上，还可以进一步给出预测值的标准差，其预测值服从如下高斯分布

(13)

其中，μG(x)为Kriging的预测值，σG(x)为预测值的标准差。

(14)

式中:Rt为目标分离可靠性。在样本点集XMC中采用如下学习函数确定失效概率最大的样本点A。

(15)

其中，μ(x)为预测的分离最大偏移距离与给定阈值之差，σ(x)为预测值的标准差。当学习函数的最小值大于2时，所有样本点的最大失效概率仅为Φ(-2)≈0.0228，满足分离可靠性的计算需求，即满足收敛条件。

式(15)的U函数更倾向于选取最大偏移距离在阈值附近或者尚未探索区域的样本点。因而，进一步将A点分为两类：若A点预测最大偏移距离在给定阈值附近，记为A1；若A点预测最大偏移距离远离给定阈值，则属于探索未知区域的情况，记为A2，如图5所示。

图5 采样策略示意图Fig.5 Illustration of the proposed sampling strategy

针对预测最大偏移距离在给定阈值附近的情况，采用新的采样策略确定失效概率更大的新增训练点，采样策略如下：首先通过A1点的符号判断其属于失效点还是可行点，若A1点属于可行点，则选取失效点集为样本点集XB，反之，选取可行点集为样本点集XB。之后，在样本点集XB中给出与A1点距离最近的样本点B，如图5所示。最后，连接A1,B两点，取N+1个采样点将线段AB等分为N段，记采样点集为XP。在采样点集XP中选取失效概率最大的采样点作为新增训练点。

针对预测最大偏移距离远离给定阈值的情况，其目的是探索未知区域，若在给定阈值附近选取失效概率更大的采样点，此时依然无法提高代理模型未知区域的精度，因而直接选取失效概率密度最大的样本点A2为新增训练点。

由于采样点集XP与当前训练点集XT不相关，因而，存在迭代后期新增训练点xnew与训练点集XT相关程度过大而影响Kriging模型精度的情况，需要进一步补充相关条件。记当前训练点集XT={x1,x2,…,xm}，则新增训练点xnew的Kriging修正函数表述如下

z(xnew)=rT(xnew)Q-1(y-F·β)

(16)

其中，y=[g(x1),g(x2),…,g(xm)]T为训练点集的最大偏移距离与给定阈值之差；Q为相关矩阵，rT(x)为新增训练点xnew与训练点集XT之间的相关向量，定义为

r(x)=[R(xnew,x1),R(xnew,x2),…,R(xnew,xm)]

(17)

其中,R为相关函数。补充两者之间的相关条件如下

max{r(xnew,XT)}<ε

(18)

其中,ε为相关程度的阈值，通常ε∈[0.9, 1]。通过上述方法，将新增训练点与训练点集之间的拥挤距离转化为两者之间的相关函数，在不影响Kriging模型精度的基础上，选取了失效概率更大的新增训练点，建立了分离可靠性极限状态方程的局部高精度代理模型，结合蒙特卡洛方法，减少分离模型的调用次数，提高分离可靠性分析效率。

综上所述，改进主动学习Kriging方法的流程图如图6所示。

图6 改进主动学习Kriging流程图Fig.6 Flowchart of the improved active learning Kriging approach

4 实例结果

以某低空高速轴对称式飞行器级间冷分离方案为应用对象，飞行高度40 km，分离冲量由两个对称安装的反推火箭提供。考虑加工及安装误差的影响，主要不确定性因素分布类型及分布参数如表2所示。

根据识别的分离不确定性模型随机选取50000个蒙特卡洛样本点，采用相同的样本点集分别采用主动学习Kriging和改进主动学习Kriging进行级间分离可靠性分析。为进一步验证改进主动学习Kriging的稳健性，给出10次分离可靠性分析结果，如表3所示。

改进主动学习Kriging的分离仿真模型调用次数更少，分析效率更高，其可靠性分析效率平均提高了16.34%。随机选取10次不同的初始训练点集和蒙特卡洛样本点集，主动学习Kriging和改进主动学习Kriging结果的标准差近似相同，校验了改进主动学习Kriging方法的稳健性。给出其中第7次仿真训练点集的分离最大偏移距离与给定阈值之差的历程对比以及U函数迭代历程对比，如图7、图8所示。

表2 主要不确定性因素分布类型及其参数Table 2 Distribution models and parameters of the significant uncertainty parameters

主动学习Kriging(AK)和改进主动学习Kriging (IAK)都趋于选择分离最大偏移距离在给定阈值附近的样本点。因而，两者在极限状态方程附近都具有较高的精度。同时，U函数迭代历程表明，改进主动学习Kriging方法提高了采样策略的加点质量，尤其是早期加点质量，使得U函数更早地进入到上升阶段，提高了可靠性分析效率。进一步，采用蒙特卡洛方法计算分离可靠性，验证方法的准确性和有效性，给出分离可靠性分析对比如表4所示。其中仿真平台的处理器为Intel(R) Core(TM) i7-6700U CPU @3.40 GHz,内存4.00 GB。

表3 飞行器级间分离可靠性分析结果Table 3 Vehicle stage separation reliability analysis results

表4 飞行器级间分离可靠性分析对比Table 4 Comparison of reliability analysis of vehicle stage separation

图7 新增训练点最大偏移距离与给定阈值之差历程对比Fig.7 Iteration history of the difference between maximum offset distance and given threshold for new training points

图8 U函数迭代历程对比Fig.8 Comparison of U function iterations

为进一步提升分离任务的可靠性，考虑不同高度对分离可靠性的影响，采用上述飞行器可靠性建模与分析方法给出33 km～44 km之间飞行器两体不发生碰撞的分离可靠性，如图9所示。

针对低空高速轴对称式飞行器级间冷分离方案，飞行高度对分离可靠性具有显著影响，由于空气密度较大，飞行马赫数较高，飞行动压较大，因而气动力及气动力矩的影响显著。选取适当的分离高度可以提高分离任务的可靠性，保证两体可靠分离，在量化分离不确定性的同时，给出分离可靠性定量分析结果，为提升分离方案的可靠性提供了重要的方法支撑。

图9 不同高度下的分离可靠性Fig.9 Separation reliability at different altitudes

5 结 论

1)本文以轴对称式飞行器级间冷分离方案为研究对象，通过考虑分离过程中不确定性因素的影响，将主要不确定性注入分离动力学仿真模型,建立了飞行器分离可靠性模型，实现了分离任务可靠性的定量分析，提升了分离可靠性的精度。

2)针对飞行器分离可靠性模型的特点，提出了改进主动学习Kriging的分离可靠性分析方法，提升了分离可靠性定量分析的效率。

3)结合实例结果验证了改进主动学习Kriging方法的正确性、高效性和稳健性，并有效减少分离可靠性分析的计算时间。下一步的研究重点是与多点并行加点策略相结合，进一步降低计算时间，提高拟合效率。