魏伟， 张竞文， 常思江， 王浩圣， 王建中

(1.北京理工大学 爆炸科学与技术国家重点实验室， 北京 100081；2.中国兵器工业第208研究所 瞬态冲击技术重点实验室， 北京 102202；.南京理工大学 能源与动力工程学院， 江苏 南京 210094)

0 引言

随着弹药技术的发展，为从根本上提高枪炮弹药的设计水平，枪弹飞行状态参数的获取变得越来越重要。枪弹飞行状态信息可通过在弹上安装各种嵌入式传感器或在外部(如枪弹靶道内)安放专用装置测量得到。外弹道理论及大量工程实践表明，枪弹发射及飞行过程中存在各种随机扰动，导致运动状态具有不确定性。定量研究姿态测量装置的测量数据对枪弹运动状态不确定性的灵敏程度，有助于测量信息的数据处理及深入开展弹药技术研究。

要定量分析灵敏程度，需要开展灵敏度分析，即采用一定方法计算出各输入参数的灵敏度值，根据灵敏度值大小筛选对于待研究参数影响较小甚至无影响的参数，设为标称值，而主要关注那些对传感器测量输出影响较大的参数。

近年来，灵敏度分析方法已逐步应用于枪炮、弹药等领域，目前已有相关研究。如文献[1]建立了一个准线性变参数刚体弹道模型，在平射条件下研究了随机扰动因素对加速度计、陀螺仪、地磁仪等装置输出的灵敏度。文献[4]利用局部灵敏度分析技术开展了炮弹气动和结构参数辨识研究。文献[5]对某运载火箭的输入不确定性及模型不确定性开展了灵敏度分析。文献[6]针对自行火炮某些参数存在相关性等情况，应用全局灵敏度方法实现了面向弹丸炮口状态的某自行火炮结构参数全局灵敏度分析。Sobol法是一种基于方差分解的敏感性分析方法，作为典型的全局灵敏度分析方法，目前在火炮弹药领域有一定的应用。如文献[7]针对混合不确定性的灵敏度分析问题，提出了一种基于辅助变量法的改进Sobol灵敏度分析方法，并据此对某电磁轨道炮的射击精度开展研究。文献[8]应用Sobol法研究了导弹分离的稳定性问题。文献[9]应用Sobol全局灵敏度分析方法定量分析了起始扰动对自行火炮射击密集度的影响。此外，Sobol法也已广泛应用于航空飞行器设计、能源系统设计等领域。扩展傅里叶振幅灵敏度检验(EFAST)法是在FAST法基础上结合Sobol法改进而来的全局敏感性分析方法。近年来，作为一种新兴的全局灵敏度分析方法，EFAST法也逐步应用于弹药领域，如文献[12]定量分析了弹、炮、气象等各种因素对密集度的影响程度，应用EFAST法对落点计算模型进行了全局灵敏度分析，该方法广泛应用于生态水文、地质探测、铁道交通等领域，具有良好的应用效果。根据文献[14]的研究，Sobol和EFAST方法是最为可靠、稳定的全局灵敏度分析方法。与Sobol法的不同之处在于EFAST方法可以得到各阶灵敏度指数，用于评价任意若干个指标之间的耦合作用，而Sobol法只能评价各个指标与全部其他指标的耦合作用。两种方法已经广泛应用于经济、能源、社会等领域大型模型中，具有一定的普遍性和很强的适用性。现有文献[7-10,12]表明，两种方法对火炮弹药等军事领域也均具有良好的适用性。由上述可知，灵敏度分析技术的应用，对于提高弹药精度、提升弹药设计水平等具有重要的作用和意义。

本文以某12.7 mm枪弹为对象，因为目前工程上没有相匹配的弹载测量装置，所以在研究中假设利用外部姿态测量手段(如靶道内连续布置高速摄影)来测得枪弹飞行过程中的转动角速度。本文选取非线性简化刚体弹道模型(称为5自由度刚体弹道模型)为分析模型，综合考虑了枪弹在发射及飞行过程中存在的初速跳动、起始扰动、气动系数偏差等随机因素，利用基于蒙特卡洛打靶的改进Sobol全局灵敏度分析方法和EFAST法，研究姿态测量装置的输出对诸扰动因素的灵敏程度，为开展枪弹姿态测量信息与其飞行状态间关系的深入分析提供基础，同时本文将两种方法进行定量比较，以期为更加深入的应用研究和方法选择提供参考。

1 外弹道模型

1.1 需用坐标系

需用坐标系主要包括基准坐标系、弹轴坐标系、弹体坐标系。

1) 基准坐标系:坐标系原点位于弹箭质心，轴沿水平线指向射击方向，轴沿铅锤向上，轴按右手法则确定为垂直于射击面并指向右方，坐标轴随质心一起平动。

2) 弹轴坐标系：坐标系原点位于弹箭质心，轴为弹轴，轴垂直于轴且指向上方，轴按右手法则垂直于平面且指向右方。

3) 弹体坐标系：坐标系原点位于弹箭质心，轴为弹轴，和固连在弹体上并与弹体一同绕纵轴旋转。

将弹轴坐标系向基准坐标系投影，即可得到弹轴坐标系和基准坐标系间的转换关系如下：

(1)

式中：为弹轴方向角；为弹轴高低角。

弹体坐标系与弹轴坐标系的差别在于坐标平面相对于坐标平面转过一个自转角，其转换关系如下：

(2)

式中：表示弹体滚转角。

本文选取姿态传感器测量枪弹转动角速度在弹体坐标系上的分量，因此只要知道枪弹转动角速度在弹轴坐标系上的分量，即可求出其在弹体坐标系三轴上的分量。

1.2 外弹道及传感器输出模型

为了从传感器输出中分析各输入因素不确定性的灵敏度，需使用枪弹运动模型，建立实测弹道与传感器之间的联系。在6自由度弹道方程数值求解过程中，为确保攻角不发散，一般采用较小积分步长(不超过0005 s)，而5自由度弹道方程略去了快圆运动可使得积分步长大为增加(如可取为004 s)，故计算时间大幅降低。本文在计算灵敏度时主要计算的弹丸姿态参数为枪弹转动角速度(滚转角速率、偏航角速率、俯仰角速率)。根据外弹道学理论，这些参数与快圆运动的关联性很小。通过试算表明，采用6自由度弹道方程和5自由度方程计算所得的枪弹转动角速度结果几乎完全相同，考虑到本文需要计算数万次弹道以得到各因素的灵敏度值，为在确保精度的前提下提高计算效率，故选取5自由度刚体弹道模型如下：

(3)

式中：表示弹箭质心的速度；和分别表示弹箭速度高低角和方向角；、、表示弹箭绕质心转动的角速度；、、表示弹箭质心在地面坐标系、、轴上的变化；为弹箭质量；、、为炮弹所受外力在弹道坐标系三轴上的分量；为轴向转动惯量；、、为外力矩在弹轴坐标系三轴上的分量，轴力矩的具体表达式可参见文献[21]。

根据5自由度外弹道模型可求出枪弹转动角速度在弹轴坐标系上的分量、和，再运用11节中坐标系转换关系可得到枪弹转动角速度在弹体坐标系上的分量，则姿态测量装置的输出可表示为

(4)

式中：、、为姿态传感器的输出。

2 两种灵敏度分析方法

灵敏度分析分为局部灵敏度分析和全局灵敏度分析，其中局部灵敏度分析运算简单，但其用于非线性模型时存在较大偏差；而全局灵敏度分析同时考虑了各因素概率密度函数的分布影响及所有因素不同的变动范围，具有较好的应用效果。为便于分析、比较，本文拟采用两种全局灵敏度分析方法，即Sobol方法和EFAST方法。

2.1 基于蒙特卡洛打靶的Sobol方法

Sobol法的核心思想是方差分解。它通过把模型分解为单个参数及参数之间相互组合的函数，计算单个输入参数或输入参数集的方差对输出方差的贡献，进而来分析参数的重要性以及参数之间的交互效应。基于蒙特卡洛打靶的改进Sobol法是利用仿真实验数据来近似单个输入参数或输入参数集的方差，从而最终得到各阶灵敏度系数。

假设可积函数()的定义域为维单元体，=[,,,…,]，Sobol法的中心思想是将函数()分解为子项之和。

(5)

()和1,2,…,(,, …,)的方差为和1,2,…,，则可得

(6)

(7)

(8)

1阶灵敏系数可以表示为

(9)

式中：为参数的1阶灵敏度，它描述了随机因素“独自”对输出总方差的贡献，越大，表明的不确定性对系统的输出影响越大。

全局灵敏度T是指参数各阶灵敏度系数之和，表示为

(10)

当输入变量的T和相差较大时，可认为变量与其他输入变量之间的交互作用对系统输出影响显著，反之可认为变量与其他输入变量之间的交互作用对系统输出影响细微。

2.2 EFAST方法

EFAST法是一种基于方差分解的全局灵敏度分析方法，结合了FAST和Sobol法的优点改进而来。其算法简单介绍如下：

设有如下模型

(11)

(12)

式中：为傅里叶变换参数；、为傅里叶振幅。

(13)

(14)

为参数的振荡频率，=1,2,…,；为每个参数的随机初相位，取[0,2π]；为标量变量，取[-π,π]。

参数变化引起模型输出方差为

(15)

式中：

(16)

由此可得模型总方差为

(17)

模型总方差可分解为

(18)

式中：为参数自身变化引起的方差；为参数通过参数变化引起的方差；12…为参数通过其余-1个参数相互作用引起的方差。

与Sobol法相同，通过归一化处理后，参数的1阶灵敏度可表示其对模型输出总方差的直接贡献，即

(19)

全局灵敏度为

(20)

式中：～为不包括参数的其他所有参数方差之和。

3 仿真分析

3.1 仿真参数

311 初始发射状态参数的不确定量

本文以某127 mm枪弹为研究对象，取初速930 m/s、平射条件下(射角为0841°)进行仿真计算。在枪口初速的变量上假设有1的不确定性范围；对于初始偏航角速率和初始俯仰角速率，在标称值0=0=0 rad/s基础上假设2 rad/s的不确定区间。枪弹的初始姿态对飞行轨迹和落点有直接影响，在初始高低角和初始方向角的标称值上均假设0000 49 rad的不确定度。初始攻角在标称值0 rad上假设0034 9 rad的不确定度。表1为所有状态变量的标称值及其不确定度。

表1 初始状态变量的标称值及不确定性区间Table 1 Nominal values and uncertainty intervals of the initial state variables

312 模型参数的不确定量

升力系数导数′和阻力系数对枪弹运动规律有很大影响，因此将其选为不确定量，′和标称值见表2。根据外弹道工程经验，在其初值上均考虑了±3的不确定区间。模型参数的初值和其不确定区间汇总如表3所示。

表2 升力系数导数和阻力系数值Table 2 Lift coefficient derivatives and drag coefficient values

根据表1和表3中所给出的初值和不确定性范围，使用Sobol法和EFAST法对影响枪弹转动角速度的因素进行灵敏度分析。该分析是以一种时变的方法进行的，每004 s(即弹道方程组的积分步长)计算一次灵敏度，根据各变量灵敏度变化趋势，判断何种输入变量对炮弹转动角速度有较大影响，同时通过计算时间和结果准确度比较两种灵敏度计算方法的优劣。

表3 枪弹气动参数的标称值及不确定性区间Table 3 Nominal values and uncertainty intervals of aerodynamic parameters of bullet

3.2 仿真计算与结果分析

根据表1和表3中的标称条件及不确定性，使用改进的Sobol法和EFAST法进行仿真。

下面对姿态测量装置三轴输出的1阶灵敏度值进行仿真。1阶灵敏度值是指单个输入参数自身对3种传感器输出的影响程度，1阶灵敏度值越大，表明该参数的不确定性对传感器的输出影响越大；而全局灵敏度值与1阶灵敏度值之差越接近于零，表明对应变量与其他输入变量之间的交互作用越小。

321 Sobol法仿真结果分析

由于Sobol法中取样较少时1阶灵敏度系数易受数值计算影响出现负值，本文经过多次取样计算测试稳定性后，假设样本均匀分布，选取取样个数为20 000。仿真采用22节所述步骤，并且由于使用了简化刚体弹道方程组，每004 s计算一次1阶灵敏度，仿真结果如图1～图4所示。

图1 Sobol法计算陀螺仪输出的1阶灵敏度值Fig.1 First-order sensitivity of the gyroscope output calculated by the EFSAT method

图1(a)中，初速对的1阶灵敏度值稳定在1附近，初速对的1阶灵敏度值远大于其余7种因素对的1阶灵敏度值，表明初速是影响输出值的主要因素。由图1(b)可以看出初速和初始攻角对姿态测量装置输出的1阶灵敏度值在0～1间波动，且初速和初始攻角这两种因素对的1阶灵敏度值在0～3 s间均大于其余因素对的1阶灵敏度值，由此推断初速和初始攻角是影响输出的主要因素，二者交替对输出产生影响。图1(c)中，初速和初始攻角对的1阶灵敏度值在0～1间波动，同时初速和初始攻角这两种因素对的1阶灵敏度值在0～3 s间均大于剩余因素对的1阶灵敏度值，表明初速和初始攻角是影响输出值的主要因素。

为研究各输入变量间的交互作用，本文还计算了8种输入变量对姿态测量装置沿弹体径向输出的全局灵敏度与1阶灵敏度的差值，以此观察各变量交互作用程度，计算结果如图2所示。

图2 Sobol法中各输入变量的全局与1阶灵敏度差值Fig.2 Global and first-order sensitivity difference of input variables in Sobol’s method

由图2可知，8种输入变量对姿态输出的全局灵敏度和1阶灵敏度差值在0～01间小幅波动，表明输入变量的独立作用是影响枪弹飞行状态的主要部分，其中初速与其他输入变量交互耦合作用所占比例最大值为7。

322 EFAST法仿真结果分析

依据EFAST法中参数的取样算法，设置取样个数分别为2 312、3 575、4 610、5 768、8 008，经过计算测试后发现EFAST在取样较少(3 575)时计算结果即趋于稳定。因此本文对不确定性参数空间进行采样，假设样本均匀分布，样本大小为=3 575组。仿真采用23节所述步骤，仿真结果如图3和图4所示。

图3 EFAST法计算陀螺仪输出的1阶灵敏度值Fig.3 First-order sensitivity of the gyroscope output calculated by the EFAST method

图4 EFAST法中各输入变量的全局灵敏度与1阶灵敏度差值Fig.2 Global and first-order sensitivity difference of input variables in the EFAST method

由图3(a)可以看出，初速对的1阶灵敏度值稳定在1附近，且初速对的1阶灵敏度值在0～3 s间远大于其余因素对的1阶灵敏度值，由此可知初速是影响输出值的主要因素。由图3(b)、图3(c)可以看出，初速和初始攻角对姿态测量装置输出、的1阶灵敏度值在0～1间波动，同时初速和初始攻角这两种因素对、的 1阶灵敏度值在0～3 s间均大于剩余6种因素对的 1阶灵敏度值，由此可知初速和初始攻角是影响、输出的主要因素。此结论与运用Sobol法的计算结果基本一致。

由图4可以看出，输入变量的独立作用是影响炮弹飞行状态的主要部分，其中初速与其他输入变量交互耦合作用所占比例最大值为8。

通过以上仿真计算可以看出：初速和初始攻角对姿态测量装置沿弹体径向的输出影响较大，初速对姿态测量装置轴向输出影响较大；输入变量的独立作用是影响枪弹飞行状态的主要部分。后续研究可以参考以上结论，开展相关的试验设计及数据处理工作，譬如通过控制枪弹初速和初始攻角得到姿态角的准确测量结果。

3.3 Sobol法和EFAST法结果比较

为直观对比不同取样数对Sobol方法输出结果的影响，本节采用22节所述步骤，应用5自由度弹道方程模型，选取取样个数为20 000和3 575，应用Sobol法，就影响输出的因素1阶灵敏度值进行仿真计算，结果如图5所示。

图5 Sobol法不同取样次数输出1阶灵敏度值Fig.5 First-order sensitivity at different sampling times in Sobol’s method

由图5(a)、图5(b)可以看出，应用Sobol方法，取样3 575次时出现了各因素输出1阶灵敏度值震荡且输出产生负值的情况。由此可知Sobol法取样较少时，1阶灵敏度系数易受数值计算影响出现负值，由此可见在相同较少的取样数时EFAST法较Sobol法输出结果更稳定准确。

从灵敏度计算结果(见图1～图4)上来看，Sobol法和EFAST法这2种方法的分析结果基本一致，即初速和初始攻角对姿态测量装置沿弹体径向的输出影响较大，初速对姿态测量装置轴向输出的影响较大。但是EFAST法较Sobol法的取样数要少得多，相应的计算过程中前者较后者要节省大量的计算时耗。在本文研究中，以姿态测量装置输出为研究对象时，EFAST法较Sobol法计算节省时间约5 min。

综上所述，EFAST法需要的样本数较少且计算高效、稳健，对于多参数模型，可较大程度地降低计算时耗，因此建议当模型参数较多时，可首选EFAST法作为研究算法。

4 结束语

本文以简化刚体弹道方程组为基础，应用Sobol法和EFAST法开展全局灵敏度分析，综合考虑初速、初始偏航角速率、初始俯仰角速率、初始高低角、初始方向角、升力系数导数、阻力系数、初始攻角等随机因素，定量计算得到了各因素对姿态测量装置输出的灵敏度值。通过对仿真结果的分析，得到如下主要结论：

1) 枪弹初速、初始攻角是影响枪弹转动角速度的主要因素。初始攻角和初速在枪弹飞行期间交替对其飞行状态产生影响。枪弹飞行状态受到各输入变量的独立及交叉耦合影响，前者是影响枪弹飞行的主要部分，交叉耦合影响相对较小。

2) 与Sobol法相比，EFAST法需要的样本数较少且计算快速准确，当模型参数较多时，可首选EFAST法开展灵敏度计算。

在后续的试验设计中，可以通过控制初始攻角、初速的散布状态，得到更为精确的枪弹飞行状态参数。