摘要
以低成本探空火箭为研究对象，基于不同学科仿真模块建立探空火箭多学科仿真模型，实现多学科耦合的火箭飞行性能仿真。针对探空火箭飞行性能的不确定性传播问题，提出一种基于动态扩充采样和代理模型的不确定性传播分析方法。基于物理分析建立了火箭飞行性能的不确定性偏差模型，通过有界扩充拉丁超立方设计方法实现偏差参数的动态采样，通过逆累积分布变化法获得满足指定分布的偏差样本。利用改进增广径向基函数模型进行飞行性能特征参数近似建模，通过调用少数样本点建立探空火箭飞行性能特征参数的近似预测模型。通过与蒙特卡罗模拟方法得到的火箭飞行性能特征参数进行对比，验证了文中方法可以在指定分布的偏差模型下通过调用少数仿真样本实现飞行性能参数统计值的快速、准确预测。
关键词
Abstract
Taking low-cost sounding rockets as the research object, a multi-disciplinary simulation model for sounding rockets was established on the basis of various disciplinary simulation modules, enabling the simulation of rocket flight performance with multi-disciplinary coupling. An uncertainty propagation analysis method based on dynamic expansion sampling and surrogate models was proposed to address the uncertainty propagation issues related to the flight performance of the sounding rocket. An uncertainty deviation model for rocket flight performance was established on the basis of physical analysis. The deviation parameters were dynamically sampled using the bounded augmented Latin hypercube design method, and deviation samples that satisfied the specified distribution were obtained through the inverse cumulative distribution transformation method. An improved expanded radial basis function model was employed for the approximate modeling of flight performance characteristic parameters, and an approximate prediction model for the flight performance characteristic parameters of sounding rockets was established by utilizing a minimal number of sample points. The flight performance characteristic parameters of rocket obtained from the proposed method were compared with those from the Monte Carlo simulation method. The results validate that the proposed method can achieve rapid and accurate statistical prediction of flight performance parameters by utilizing a minimal number of simulation samples under the specified distribution deviation model.
低成本探空火箭的研制对准确获取临近空间大气环境原位探测数据、提升中远程导弹武器打击精度具有极其重要的意义。当前临近空间大气环境原位探测手段较为缺乏[1]，自膨胀落球技术取得突破后，落球探测方法得以迅速发展，并催生了低成本的标枪型火箭[2]。得益于原位探测功能，标枪型探空火箭成为临近空间大气探测不可或缺的重要手段[3]。
然而，标枪型探空火箭多为无控自稳定火箭，在执行飞行任务时可能受大气风场变化、设计缺陷及制造误差的影响，发生失稳或飞行高度不满足设计指标的现象。因此，有必要对探空火箭飞行性能参数进行不确定性量化分析。通过计算不同飞行环境和火箭参数偏差对探空火箭实际飞行性能的影响，评估在实际飞行过程中的不确定性影响下，火箭的弹道最大高度、飞行过程中的最大攻角、侧滑角和最大动压等飞行性能参数是否满足设计要求，从而为实际飞行任务提供可靠参考。
探空火箭飞行性能不确定性量化问题可具体分为不确定性建模和不确定性分析。其中，不确定性建模包括系统参数建模和不确定性参数建模，不确定性分析包括不确定性传播和灵敏度分析等[4]。本研究针对探空火箭飞行时大气环境和发动机参数偏差，开展不确定性传播分析方法的研究。传统不确定性传播分析方法包括采样法、局部展开法等。采样法主要指蒙特卡罗模拟（Monte Carlo simulation，MCS）技术及其变体[5]；局部展开法应用混沌多项式（polynomial chaos expansion，PCE）表示随机变量[6]，但该种方法的计算量随维度上升会急剧增加。
不确定性传播需要量化输入不确定性对输出产生的影响。探空火箭系统具有复杂的结构和强耦合性，使用仿真模型计算其飞行性能会导致巨大的计算负担，因此常常使用效率更高的近似模型来替代仿真模型以提高计算效率。根据实施方式的不同，可以将不确定性传播的近似方法分为嵌入式和非嵌入式方法[7]。嵌入式方法是指将近似模型嵌入系统方程中，通过修改原模型并加入不确定性项来进行建模。非嵌入式方法将系统视为黑箱模型，在不修改系统的情况下，通过分析模型的输入和输出来进行计算。非嵌入式方法避免了方法与对象之间的耦合关系，具有更广阔的应用前景[8]。
综上所述，评估探空火箭在不确定性影响下飞行性能可否满足任务设计指标具有现实意义。同时，传统的分析方法会面临大量迭代重复、计算量庞大、“维度诅咒”以及仿真高耗时等问题[9]。对此，本研究提出了基于动态扩充采样和代理模型的探空火箭性能不确定性传播分析方法。该方法基于动态扩充采样减少高耗时仿真模型的调用次数，并通过近似建模技术提升了低成本探空火箭飞行性能不确定性传播分析的效率。
1 探空火箭多学科耦合计算模型
1.1 几何构型模型
探空火箭采用一级动力和一级载荷方案，从前到后依次为载荷级（标枪级）、级间段、固体火箭发动机和尾舱。动力级和标枪级采用气动分离方案，在发动机工作结束后，采用动力级和标枪级截面差产生的气动阻力差完成分离。标枪级依靠惯性继续飞行至弹道顶点后释放膨胀落球。根据总体指标和低成本需求，确定标枪型探空火箭总体布局如图1所示。
图1探空火箭总体布局
Fig.1Overall layout of sounding rocket
火箭利用翼柱型药柱的单室单推固体发动机提供动力，依靠静稳定尾翼实现全程零攻角飞行。在动力级燃药耗尽时，依靠气动阻力差实现标枪级和动力级的分离，减少火箭在惯性段飞行时气动阻力造成的速度损失，同时简化两级的分离结构、降低火箭的造价。探空火箭的主要设计参数如表1所示。
表1探空火箭主要参数
Tab.1Main parameters of a sounding rocket
1.2 气动特性计算模型
气动特性对探空火箭总体性能具有重要影响。根据飞行弹道估算，确定分离前的弹道高度范围为0~10 km，马赫数范围为0~6.0；标枪级火箭的弹道高度在5 km以上，马赫数范围为0~6.0。针对所得飞行工况范围，利用计算流体力学（computational fluid dynamics，CFD）进行火箭流场仿真，仿真计算的方法和条件设置如表2所示。
表2CFD参数设置
Tab.2Parameters setting of CFD
采用笛卡儿网格并局部加密，计算得到火箭不同高度、不同马赫数以及不同姿态角时的气动力系数和气动力矩系数，网格结构如图2所示。流体计算域的前后边界距离箭体头部均为5倍箭长，即整个流场的计算区间设置为箭体头部的前方5倍弹长到弹体尾部后方4倍弹长。为展示网格结构，图2对网格进行了放大处理。为了保障网格的质量和计算精度，在弹体附近区域采用了O形网格设计，并对模型表面特别是头部和尾翼附近区域进行了精密加密处理，其中第一层厚度为0.1 mm，网格总数约为300万。
图2结构网格示意图
Fig.2Schematic diagram of the structural grid
1.3 动力系统计算模型
本研究中，低成本探空火箭采用一级固体火箭发动机作为动力系统，药柱类型为翼柱型，推进剂类型为三组元丁羟复合推进剂，其结构示意如图3所示。
图3固体发动机结构示意图
Fig.3Schematic diagram of solid rocket motor structure
固体发动机工作过程涉及燃烧流动和传热，内部工作过程十分复杂，为了准确描述这些复杂的过程，需要采用合适的发动机内弹道数学模型。已有的数学模型包括零维、一维、二维和三维内弹道方程。为满足快速仿真要求，本研究采用如式（1）所示的零维内弹道方程组对动力系统性能进行描述[10]。
(1)
其中：Vc为燃烧室总自由容积；为燃气平均密度；ρP为推进剂密度；S为装药总燃面面积；为平均燃速；φ2为流量修正系数；Γ为比热比的函数；为燃烧室平均压强；At为喉部面积；R为气体常数；为平均温度；χ为散热系数；k为比热比；Tp为推进剂定压燃烧温度；a为燃速系数；n为燃速压强指数。
根据所构建的零维内弹道计算模型求解发动机喷管质量流量，其计算公式如式（2）所示。
(2)
式中，χRT0是考虑了散热系数的点火药力，为质量流率损失系数，p0为喷管中的滞止压强，c*为特征速度。推力的计算如式（3）所示。
(3)
式中，CF为推力系数，为推力系数效率。推力系数CF由式（4）计算：
(4)
式中，γ为比热比，ζe为膨胀比，πe为喷管出口压强与燃烧室平均压强的比值，πa为环境压强与燃烧室平均压强的比值。其中燃烧室平均压强根据零维内弹道方程（1）确定。
1.4 六自由度弹道仿真模型
本研究采用半速度坐标系中六自由度无控弹道仿真动力学模型进行低成本探空火箭弹道计算，其质心动力学方程表示为式（5）。
(5)
其中，m、v、θ、σ、ν、α和β分别为探空火箭质量、速度、速度倾角、航迹偏航角、倾侧角、攻角和侧滑角，Pe为发动机推力，Cx、Cy和Cz为三个方向的气动力系数，q=0.5ρv2为动压，g为重力加速度。
弹体坐标系中的绕质心动力学方程为式（6）。
(6)
其中，、和分别为三个方向的转动惯量，、和分别为绕三轴的角速度，、和分别为三个方向的力矩。
质心运动学方程如式（7）所示。
(7)
绕质心运动学方程如式（8）所示。
(8)
式中，φ为箭体俯仰角，ψ为箭体偏航角，ξ为箭体滚转角。
补充方程如式（9）所示。
(9)
1.5 多学科模块集成与仿真
低成本探空火箭总体方案涉及几何、气动、动力和弹道等多个学科。其中几何构型为CFD气动仿真提供参数基础，动力学科采用固体发动机零维内弹道计算模型，弹道学科采用六自由度无控弹道仿真模型。对此，本研究整合并调用多学科专业计算工具，构建探空火箭多学科联合仿真流程，实现探空火箭飞行性能的快速、多学科耦合仿真计算，其仿真流程如图4所示。
图4探空火箭多学科仿真流程
Fig.4Multidisciplinary simulation process of sounding rocket
根据1.1节的火箭总体参数，对探空火箭进行多学科联合仿真。其中，基于CFD气动仿真软件和探空火箭气动构型，结合火箭飞行工况计算得到的火箭零攻角阻力系数Cd随大地高度H、马赫数变化如图5、图6所示，标枪级和全箭的升力、力矩系数如表3、表4所示。考虑到40 km以上大气密度不到海平面大气密度的3‰，所以40 km高度以上的升、阻力系数和力矩系数通过三次样条法对数据进行外插计算。
图5全箭阻力系数
Fig.5Drag coefficient of entire rocket
图6标枪级阻力系数
Fig.6Drag coefficient of javelin stage
表3全箭升力、力矩系数
Tab.3Lift and moment coefficients of entire rocket
表4标枪级升力、力矩系数
Tab.4Lift and moment coefficients of javelin stage
探空火箭采用一级动力方案，固体发动机为翼柱型药柱的单室单推力型，其不同温度时推力仿真曲线如图7所示。
图7发动机推力-时间曲线
Fig.7Thrust-time curve of engine
基于探空火箭总体参数，利用多学科联合仿真模型进行火箭六自由度弹道仿真，其大地高度、马赫数六自由度仿真曲线如图8所示，攻角、侧滑角六自由度仿真曲线如图9所示。
图8大地高度、马赫数六自由度仿真曲线
Fig.8Geodetic altitude and Mach number six-degree-of-freedom simulation curve
图9攻角、侧滑角六自由度仿真曲线
Fig.9Angle of attack and sideslip angle six-degree-of-freedom simulation curve
2 基于BRELHD和ARBF模型的不确定性传播分析方法
基于近似模型的不确定性分析方法中，不确定性传播通过量化随机输入系统对输出的影响，描述输出的统计分布。现有传播分析方法难以判断算法收敛性，且大多不适用于随机输入存在多种分布的情形。对此，本研究提出一种基于动态扩充采样和代理模型的不确定性传播分析方法，统一不同输入变量的概率分布形式，有效判断算法收敛程度，显著提高飞行器不确定性传播分析效率。
2.1 有界扩充拉丁超立方实验设计方法
近似建模过程中，广泛采用拉丁超立方设计（Latin hypercube design，LHD）作为实验设计方法。然而，现有LHD方法大多难以兼顾设计效果和计算效率，且多为样本数量需事先给定的一次性采样方法。近似模型的精度与训练数据规模密切相关，样本冗余会导致资源浪费和过度训练，样本不足则难以保证模型预测精度。因此，针对如何动态生成均匀性良好的实验设计问题，有学者[11-12]提出一种递归演化拉丁超立方设计（recursive evolution permutation Latin hypercube design，RELHD）方法，实现样本规模动态线性扩充，解决了最优样本规模难以提前预估的难题，提出了面向不同随机变量均匀采样的有界序列扩充改进方案，从而实现不同随机变量训练数据高效获取。
RELHD方法中采样点被设置为可以到达设计空间的边界，该方法的2维采样扩充示例如图10（a）所示。图中的黑色圆点和红色圆点分别表示初始设计和演化设计中的样本点。给定初始样本数为p，此时填充演化空间所需的种子数量为p-1。
然而，针对飞行器不确定性传播分析问题，输入变量往往不符合均匀分布，难以通过简单设置使样本达到设计边界。因此，对RELHD中的扩充方法进行改进，使其更适用于不确定性传播问题。扩充过程中，初始的p个样本会将整个设计空间划分为p+1份。上述情况下，要完成整个演化空间的填充则需要p+1个种子，具体示意如图10（b）所示，可见由于初始空间未到达边界，整个采样空间划分数相较于图10（a）中空间划分数增加了1。
图10不确定条件下演化空间对比
Fig.10Evolutionary space comparison under uncertain conditions
对于飞行器不确定性传播分析问题的样本点扩充方式，其样本空间一维投影如图11所示，由于每轮扩充后的总长度仅比上一轮多出间隔的1倍，这种扩充方式的总区间长度是有边界的，因此将其叫作有界扩充LHD方法。
图11有界扩充一维投影示意
Fig.11One-dimensional projection schematic of bounded expansion
针对有界扩充LHD方法的有界性展开证明：设初始样本点个数为p，样本初始间隔为δ1，则初始样本包含p-1个间隔，即设计域宽为（p-1）×δ1。扩充k次后，由图11可知，第k次的总宽度相较于第k-1次只增加δk-1，则扩充k次后的总长度与第k-1次扩充后的总长度关系为：
(10)
又因，式（10）转化为：
(11)
通过等比数列求和可得：
(12)
则：
(13)
极限状态下：
(14)
结合上述证明可知，所提出的扩充方式有界，将其命名为RELHD的有界扩充型（bounded RELHD，BRELHD）。
应用数值算例，比较增强随机进化算法（enhanced stochastic evolutionary algorithm，ESE）、平移传播算法（translational propagation algorithm，TPA）和lhsdesign三种常见流行算法与BRELHD算法在样本空间均匀性和计算效率上的差异。以上算法均在相同配置的计算机上独立运行。为降低随机影响，以每种算法独立运行50次的平均值作为最终结果。
表5中考虑了Viana等[13]提出的采样设置问题。BRELHD以小样本ESE的结果作为初始设计并进行序列扩充，更适合大中型样本设计。因此，主要研究维度d范围为2~8的大中型采样问题以实现BRELHD的快速合理评估。
表5LHD的采样设置
Tab.5Sampling settings of LHD
图12给出了大中型采样问题空间均匀性能的比较结果。可以看出：lhsdesign的性能最差且随机性较强；TPA在维度较低且样本数量较大时表现良好，随着维度的增加，其竞争力下降；在大多数情况下，ESE与BRELHD有着相似的空间均匀性，且明显优于其他两种方法；BRELHD在绝大多数问题上有着较优的表现，其优势随着样本数和维度的增加而更加明显。
图12大中型采样问题对比
Fig.12Comparison of medium and large-scale sampling problems
由于BRELHD和ESE的空间均匀质量在大多情况下较为接近，且BRELHD方法的初始样本是通过ESE生成，因此需对BRELHD和ESE的计算效率进行对比，以验证BRELHD性能。
图13显示了不同维度下（2，5，10，20）、样本数为100和200时计算效率的比较结果。从图中可以看出，得益于LHD的直接构造，TPA和lhsdesign相较于BRELHD耗时更少。然而，TPA和lhsdesign的空间均匀性往往并不令人满意。在满足良好均匀性的前提下，BRELHD的计算效率要明显优于ESE，且随着维数的增加，BRELHD的优势愈发明显。
图13不同维度下算法耗时的比较
Fig.13Comparison of algorithm execution time in different dimensions
2.2 改进增广径向基函数近似模型
不同基函数应用于不同问题表现差异明显，模型基函数的选取很大程度上决定了模型的泛化性能。因此，为了利用不同基函数特点提高模型预测精度，构造基于PCE的增广径向基函数[14]（augmented radial basis function，ARBF）如下：
(15)
式中，λj和ωi分别为PCE和RBF的系数，Qj（x）和φi（x）为对应基函数。注意到，φi（x）为Gauss基函数，如
(16)
式中，ci为第i个基函数的形状参数。
通过求解式（17）来计算系数λ和ω：
(17)
式中：Φ、Q和P分别为φi（x）、Qj（x）和Pij组成的系数矩阵，Pij为Qi（x）与Qj（x）的加权内积项；O为元素均为0的矩阵。
由于样本的空间分布不均匀，基于样本密度选择合适的形状参数可以有效提高元模型的精度。采用样本的局部密度来表征样本的稀疏性，其表示如下：
(18)
其中，与形状参数类似，表示样本对局部密度的影响。将训练样本的影响体积Vij归一化，得到每个样本和维度对应的形状参数参考核宽度为：
(19)
最终的核宽度由比例因子λj乘以得到：
(20)
因此，式（7）中的各向异性表达式可构造为：
(21)
最终各向异性方法改进的增广径向基函数模型为：
(22)
采用快速交叉验证方法对模型的泛化误差进行评估，该方法的预测误差Δy如式（23）所示。
(23)
式中：Ypre，r和Yr分别为r个训练样本的预测值和实际输出；ωr为对应权重系数；L*r×r为由Φ、Q、O和P组成并变换得到的系数矩阵，具体推导见文献[14]。
为验证近似模型的有效性，分别用三个标准算例进行测试，并与目前常用方法进行比较。算例信息如表6所示。为揭示问题维度对模型精度的影响，将训练样本数设置为设计变量的10倍和20倍。所有测试近似模型的训练样本均由RELHD方法生成。为全面评估模型准确性，将测试样本数设置为问题维度的104倍。采用测试样本的均方根误差（root mean square error，RMSE）评估不同近似模型性能：
(24)
式中，xi为第i个样本点，N为测试样本大小，y（xi）为样本xi的真实响应，为近似模型获得的预测输出。
表6三个基准测试函数
Tab.6Three benchmark test functions
统计重复10次测试的RMSE，以降低随机误差影响，计算其平均值，测试结果如表7和图14所示，将ARBF与其他近似模型测试的RMSE值进行比较。随着训练样本量的增长，大多数近似模型更接近真实模型，ARBF方法在不同基准函数的精度和稳定性方面比其他常用近似模型表现出更好性能。在高度非线性基准函数中，大多数近似模型无法逼近真实模型，而ARBF可得到较好结果，表明ARBF近似模型在逼近非线性模型时具有十分优越的性能。
表7测试函数结果
Tab.7Test function results
图14不同近似模型的RMSE结果箱型图
Fig.14Box plot of RMSE results for different approximation models
2.3 算法流程
目前的实验设计和近似建模方法通常对样本集的分布类型有一定的要求，例如通常要求样本来自常见的均匀分布。然而，在不确定性问题的研究中，很难满足这种要求。为了提高基于近似模型的不确定性传播方法的适用性，引入了逆累积分布变换方法，该方法可以将来自任意区间均匀分布的样本转换为符合不同概率分布的样本。
逆累积分布变换方法通常被认为是一种非均匀概率密度的采样方法，用于生成符合指定概率密度分布的样本集。下面以标准正态分布为例来说明该方法的基本流程。对于服从概率密度函数的随机变量x，其累积分布函数（cumulative distribution function，CDF）如式（25）所示：
(25)
使用累积分布函数G（x），可以将服从指定概率密度函数的随机变量集映射为均匀分布在区间[0，1]中的随机变量集。相反，逆累积分布函数（inverse cumulative distribution function，ICDF）G-1（x）可以将区间[0，1]中的随机变量集映射到服从指定概率密度函数的任意区间的随机变量集。基于这种转换，可以采用常规实验设计方法生成近似模型所需的样本集。归一化和ICDF变换采样集的过程如图15所示。
图15采样集的变换过程
Fig.15Transformation process of the sampling set
基于ARBF近似建模的不确定性传播分析方法（ARBF-based uncertainty propagation，ARBF-UP）的算法流程如下。
步骤1：使用BRELHD方法生成初始实验设计，采样点数设为p，则根据式（14）和样本归一化值，初始的设计点间隔δ取值为。
步骤2：通过ICDF将初始设计变换为符合实际概率分布的采样点。
步骤3：将采样点代入仿真模型并获得输出，生成训练样本。
步骤4：基于训练样本构建ARBF近似模型，计算样本点处近似模型的预测值，与仿真进行对比，计算仿真值与预测值间的相对均方根误差ε。当ε＜0.01时，认为模型精度满足需求。ε的计算公式为：
(26)
其中，EN为新扩充的样本集，xk和yk分别为新扩充样本集第k个输入向量和仿真值，为未扩充样本前建立的近似模型，为样本集中仿真值的均值。
步骤5：若精度要求不满足，则基于BRELHD对训练样本集进行扩充，扩充过程中，p个样本会将整个设计空间划分为p+1份，其余的扩充机制与递归演化拉丁超立方实验设计方法[11-12]保持一致。扩充完成后，转至步骤2，直至精度满足。
步骤6：若满足精度需求，则基于所构建的ARBF近似模型进行不确定性传播分析。
在这个流程中，首先通过ICDF将初始设计样本转换为符合实际概率分布的采样点，可以更准确地反映真实系统的特性。然后，通过使用ARBF近似模型来快速预测输出，并根据需要对训练样本进行扩充以提高模型精度。最后，基于构建的ARBF近似模型，使用蒙特卡罗仿真方法进行不确定性传播分析，以评估系统输出的不确定性。整个算法的流程如图16所示。
图16基于增广径向基函数的不确定传播分析方法
Fig.16Uncertainty propagation analysis method based on augmented radial basis functions
3 探空火箭飞行性能不确定性传播分析
在火箭系统的分析和设计过程中，面对大气环境的变化和发动机参数偏差等因素的影响，飞行器的性能指标可能会发生变化，从而导致轨迹偏差、故障甚至失效。因此针对上述不确定性参数偏差，提出一种高效的不确定性传播分析方法对火箭飞行性能的不确定性评估具有重要意义。
3.1 探空火箭不确定参数偏差建模
在火箭发射和飞行过程中，火箭飞行性能主要受大气密度、大气风场和发动机参数不确定性偏差的影响，因此本研究针对这些参数进行偏差建模。
3.1.1 大气密度偏差建模
大气密度不是一成不变的，相反，随着时间的变化，大气密度会有一定程度的改变。而大气密度改变直接影响着火箭所受气动力的大小，进而会影响火箭的飞行弹道和飞行姿态。查阅相关文献得知，同一高度下大气密度在不同时间的变化范围一般在10%左右[15]，因此在六自由度弹道程序的大气密度插值模块中引入大气密度变化系数。其取值服从均值为0、标准差为0.033的正态分布，以此近似模拟10%的大气密度随机偏差，进而可以通过大气密度变化系数的改变来模拟六自由度仿真中大气密度场的改变。
3.1.2 基于实测数据的大气风场模型
本研究中大气风场模型为基于实测数据建立的大气风场分段插值模型。共收集了4年内1 310组不同时刻发射场附近大气风速、风向在不同高度的观测数据。选取了5组时间跨度较大的观测数据，如图17所示。同时，统计了1 310组大气观测数据不同高度处风速、风向的均值及标准差，如图18所示。由观测数据可知，虽然观测时刻不同，但大气风速随高度的分布表现出较为一致的规律：在0~12 km，随着高度增大，大气风速呈增大趋势；在12~20 km，随着高度增大，大气风速呈减小趋势；在20~30 km，大气风速呈增大趋势但增速较小。经分析，这与大气对流层、平流层的分布有关。大气风向均值在150°以上，且在0~20 km，风向均值大于200°，表明观测地点多为西南风。
图17大气风场变化样本观测数据
Fig.17Observational data of atmospheric wind field variation samples
图18大气风场观测数据统计值
Fig.18Statistical values of atmospheric wind field observation data
考虑到高空大气稀薄，对火箭飞行产生的气动力影响减小，当高度大于30 km时，大气风速、风向取值与30 km处观测数据相同，为常值风场。利用上述基于大气风场观测数据的风场模型，在进行探空火箭六自由度弹道仿真时，可随机加入基于某观测时刻的大气风场，以模拟现实环境中大气风场的随机变化。
3.1.3 发动机参数偏差建模
在实际飞行过程中，发动机某些参数的不确定变化会对发动机推力变化产生较大的影响，进而影响火箭实际的飞行弹道，因此有必要对发动机的不确定参数进行偏差建模以评估发动机参数变化对火箭飞行弹道的影响。基于固体发动机原理可知，发动机加工总装完成后，由于发射时环境状态的不同，装药的燃速、特征速度、喷管的烧蚀速率会产生较小的随机变化，导致发动机推力改变。因此本研究针对装药燃速ι、特征速度C*x和喷管烧蚀速率v*进行偏差建模。假设上述参数符合正态分布，取当前修正后的发动机装药燃速、特征速度、喷管烧蚀速率为均值，取参数值的5%为三倍标准差，则各参数具体分布情况如表8所示。
表8发动机偏差参数分布
Tab.8Distribution of engine deviation parameters
3.2 探空火箭弹道性能不确定性传播分析
下面基于文中提出的方法对探空火箭飞行性能进行不确定性传播分析。所得结果与MCS方法对比，以验证提出方法的有效性。对于总体方案设定的初始条件，将3.1节的不确定性耦合到多学科联合仿真模型中，应用MCS方法进行5 000次随机采样，对样本进行六自由度弹道仿真。由于引入不确定性，弹道仿真结果会出现一定偏差，得到弹道最大高度、最大攻角、最大侧滑角和最大动压飞行性能分布如图19所示。
表9给出了仿真结果的分布特征，包括均值、标准差和变异系数cv。使用变异系数cv表征数据的离散程度，cv计算公式如下：
(27)
式中，σdata和μdata分别表示数据的标准差和均值。
本研究中的探空火箭为无控自稳定火箭，其设计高度要求不低于115 km。从分布直方图可以看出，火箭弹道最大高度在110~140 km之间，且主要分布在120~130 km之间。经统计，99.7%的样本点仿真弹道最大高度大于115 km，因此认为探空火箭的飞行高度可以满足设计需求。此外，其最大攻角和最大侧滑角均分布在0~5°之间，特别是最大侧滑角不会大于1°，说明火箭在飞行过程中可以满足稳定性需求。注意到最大攻角和最大侧滑角的变异系数分别为68.9%和8.16%，说明其分布受不确定性影响较大。火箭最大动压分布在1.47~1.71 MPa之间，表明火箭在飞行过程中对结构强度和防热有较高的要求。经气动热仿真显示，两级火箭在4 km、Ma=6分离点时面临最严苛的气动热环境。图20展示了两级火箭在分离点时的表面温度分布，箭体温度均在1 800 K以上，其中相对高温区主要集中在标枪级火箭头部（1 930 K）、级间连接扩张群段（1 960 K）、动力级火箭箭体表面（1 950 K）、标枪级火箭尾翼根部（1 930 K）和助推级火箭根部（1 970 K）。考虑到该工况工作时间较短，将采用绝热涂层的方式在不改变总体质量的情况下，提高火箭的热防护能力。
图19探空火箭飞行性能分布
Fig.19Distribution of sounding rocket flight performance
表9输出参数分布特征
Tab.9Distribution characteristics of output parameters
注：表中均值和标准差均只保留了三位小数，变异系数是由原始数据计算而来。
图20分离点箭体表面温度
Fig.20Separation point body surface temperature
应用MCS方法获得5 000个样本点进行探空火箭六自由度高耗时弹道仿真，计算不同采样点数下弹道最大高度、最大攻角、最大侧滑角和最大动压的均值和标准差，如图21和图22所示。可见在采样点数大于1 000时，弹道特征参数的均值和标准差基本收敛，足以衡量探空火箭飞行性能的动态分布。为确保精度，选择MCS采样点数为5 000时的统计结果作为真值与近似建模结果进行对比。
基于BRELHD方法进行动态扩充采样，利用ARBF模型进行参数建模，初始采样点数为40。发现当采样点数为200时，近似模型的相对均方根误差ε满足收敛条件。为验证BRELHD方法和逆累积分布变换法的有效性，以大气密度正态偏差的样本为例，200组BRELHD方法采样点逆累积分布变换后统计直方图如图23所示，可见采样点取值较好地符合正态分布规律，且由于不会在样本边界取值，逆变换后的样本值不会取到正负无穷，满足了生成正态分布偏差的仿真需求。
图21不同采样数下弹道参数均值变化
Fig.21Variation of mean ballistic parameters under different sample sizes
为验证方法有效性，针对探空火箭飞行性能参数统计值，将扩充过程中样本点数为40、80、120、160和200时建立的ARBF近似模型预测结果与5 000次MCS方法仿真结果进行对比。最终得到不同样本点数下弹道最大高度、最大攻角、最大侧滑角和最大动压飞行性能参数预测结果和MCS方法仿真结果，其对比如图24、图25所示，具体的统计值大小（保留三位小数）如表10所示。从图表中可以看出，基于200组采样点建立的ARBF近似模型的预测结果统计值与MCS方法进行5 000次飞行仿真所得到结果较为接近。从图中的误差带和表10中的数据对比可知，飞行性能参数统计均值的预测精度可以达到1%甚至1‰量级，标准差的预测精度也可以控制在10%以内，表明利用200组飞行仿真样本建立的近似模型可实现探空火箭飞行性能参数统计值的高精度预测，其精度可达到基于5 000次蒙特卡罗仿真所获得的结果。随着样本数的增加，ARBF近似模型通过BRELHD方法实现动态更新，便于判断预测结果的可行性。至此验证了本研究算法的准确性与高效性。
图22不同采样数下弹道参数标准差变化
Fig.22Variation of ballistic parameter standard deviation under different sample sizes
图23采样点取值分布
Fig.23Distribution of sampling point values
图24MCS与近似建模方法飞行性能统计均值对比
Fig.24Comparison of flight performance statistical means between MCS and approximation modeling methods
图25MCS与近似建模方法飞行性能统计标准差对比
Fig.25Comparison of flight performance statistical standard deviations between MCS and approximation modeling methods
4 结论
本研究通过分析火箭各学科间的耦合关系和各学科仿真特点，以探空火箭为研究对象建立了低成本探空火箭的多学科联合仿真流程，通过梳理各仿真模块间的工作流程，实现探空火箭飞行性能的高效计算。
而后，针对探空火箭飞行性能不确定性传播分析问题，开展了基于动态扩充采样和代理模型的不确定性传播分析方法的研究。针对如何动态生成均匀性良好的实验设计问题以及在飞行器不确定性传播分析过程中输入变量往往不符合均匀分布的问题，提出了面向不同分布随机变量均匀采样的有界序列扩充拉丁超立方实验设计方案，实现了不同随机变量训练数据高效获取。
进一步，为了提升近似模型在面对不同问题时的泛化性能，兼顾PCE和RBF模型的优势，建立了改进增广径向基混合近似模型，并将其应用到探空火箭的不确定性传播分析问题中。通过确定影响探空火箭飞行性能的偏差参数建立不确定性偏差模型，基于BRELHD方法对偏差参数进行动态扩充采样并建立改进增广径向基混合近似模型。
表10探空火箭飞行性能不确定性传播分析结果比较
Tab.10Comparison of uncertainty propagation analysis results for the sounding rocket
最后，通过近似模型对火箭飞行性能特征参数的均值、标准差进行预测并与传统MCS方法的结果进行算法精度对比。仿真结果表明，利用文中提出的方法，调用200次飞行仿真模型建立ARBF近似模型，其预测精度与传统蒙特卡罗方法调用5 000次飞行仿真模型获得的结果精度相近，证明了本研究提出的探空火箭飞行性能不确定性传播方法的有效性。
