炮弹飞行运动状态的全局灵敏度分析

张竞文,邵渤涵,常思江,谢 飞

(1.南京理工大学 能源与动力工程学院,江苏 南京 210094;2.中国船舶集团有限公司系统工程研究院,北京 100094;3.辽沈工业集团有限公司,辽宁 沈阳 110045)

随着信息化弹药技术的不断深入和发展,炮弹飞行状态参数的获取变得越来越重要和迫切。炮弹飞行状态信息主要通过在弹上安装各种嵌入式传感器测量得到,而炮弹发射及飞行过程中存在各种随机扰动,导致运动状态具有不确定性。定量研究嵌入式传感器数据对炮弹运动状态不确定性的灵敏程度,有助于测量信息的数据处理及深入开展信息化弹药技术的研究。所谓灵敏度分析,主要是采用一定的方法,定量计算出各输入参数的灵敏度值,根据灵敏度值的大小筛选出对于待研究参数没有影响或者影响较小的参数,将其设为标称值,而将主要研究精力或关注点集中于讨论或设计对传感器输出影响较大的参数上。

近年来,灵敏度分析方法已逐步应用于弹药、火炮等领域,目前已有一些相关的应用研究。文献[1]建立了一个准线性变参数刚体弹道模型,据此在平射条件下研究了一些随机扰动因素对弹载测量装置输出的灵敏度;文献[6]引入局部灵敏度分析技术,开展了炮弹的气动和结构参数辨识研究;文献[7]提出一种解析目标级联方法,对某运载火箭的输入不确定性和模型不确定性进行了全局灵敏度分析;文献[8]则针对混合不确定性的灵敏度分析问题,提出了一种基于辅助变量法的改进Sobol’灵敏度分析方法,并据此对某电磁轨道炮的射击精度开展研究。文献[9]针对自行火炮某些参数存在相关性等情况,提出了一种基于Rosenblatt转换理论的全局灵敏度分析方法,实现了面向弹丸炮口状态的某自行火炮结构参数全局灵敏度分析。文献[10]为评估、设计、提高自行火炮的射击密集度,应用Sobol’全局灵敏度分析方法定量分析了起始扰动对自行火炮射击密集度的影响。

本文以某大口径炮弹为研究对象,假设弹上安装了陀螺仪、加速度计及磁力计3种传感器,分别用于测量炮弹转动角速度、加速度和姿态角。选取非线性简化刚体弹道模型(称为五自由度刚体弹道模型)为分析模型,综合考虑了炮弹在发射及飞行过程中存在的初速跳动、起始扰动、弹丸质量偏差、气动系数偏差等随机因素,利用基于蒙特卡洛打靶的改进Sobol’全局灵敏度分析方法,研究陀螺仪、加速度计及磁力计的三轴输出对诸扰动因素的灵敏程度,为开展弹载传感器测量信息与炮弹飞行状态间关系的深入分析提供基础,与文献[1]相比,本文较为全面地考虑了扰动因素(如起始攻角和滚转阻尼力系数等),并扩展了分析条件(如射角和飞行时间),以期为深入的应用研究提供参考。

1 动力学模型

炮弹姿态角的定义和建立炮弹弹道方程组所使用的坐标系及其转换关系等密切相关,姿态角的解算则依赖于不同坐标系的转换关系。

1.1 坐标系及其相互转换

本文用到的坐标系包括:地面坐标系、基准坐标系、弹轴坐标系、弹体坐标系。

①地面坐标系:坐标系原点通常选取在射出点,轴通常是弹道面与水平面交线,指向目标为正;轴沿垂线向上,轴与其他两轴垂直并构成右手坐标系。

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

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

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

将弹轴坐标系向基准坐标系投影,即可得到弹轴坐标系和基准坐标系间的转换关系:

(1)

式中:为弹轴方向角;为弹轴高低角。

弹体坐标系与弹轴坐标系的差别只在于坐标平面相对于坐标平面转过一个自转角,转换关系为

(2)

式中:为弹体滚转角。

本文选取三轴陀螺仪测量炮弹转动角速度在弹体坐标系上的分量,三轴磁力计测量炮弹在弹体坐标系上的姿态角,三轴加速度计测量炮弹在弹体坐标系上的加速度。只要知道炮弹转动角速度、地磁要素、弹箭加速度在弹轴坐标系上的分量,即可转换求出其在弹体坐标系三轴上的分量。

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

为了从3种传感器输出中分析各输入因素不确定性的灵敏度,需使用炮弹运动模型,建立实测弹道与陀螺仪、加速度计、磁力计输出之间的联系。由于普通六自由度弹道方程的求解步长较小(一般为0.002 s),而五自由度弹道方程由于略去了快圆运动可使得积分步长大为增加(如可取为0.04 s),所以计算时间几乎只有普通六自由度弹道方程的1/10。考虑到本文以分析炮弹飞行状态为主,且需计算数十万次弹道以得到各因素的灵敏度值,为在确保计算精度的前提下提高计算效率,选用如下五自由度刚体弹道模型:

(3)

式中:为弹箭质心的速度;和分别为弹箭速度高低角和方向角;和分别为弹轴高低角和方位角;为弹体滚转角;,,为弹箭绕质心转动的角速度;,,为弹箭质心在地面坐标系的位置分量;为弹箭质量;,,为炮弹所受外力在弹道坐标系三轴上的分量;为轴向转动惯量;,,为外力矩在弹轴坐标系三轴上的分量;轴力矩的具体表达式可参见文献[11]。

根据五自由度外弹道模型可求出炮弹转动角速度在弹轴坐标系上的分量、和,再运用第11节中坐标系转换关系可得到炮弹转动角速度在弹体坐标系上的分量,则三轴陀螺仪的输出可表示为

(4)

式中:,,为陀螺仪三轴上的输出。

由五自由度外弹道模型可求出弹箭质心的速度,运用坐标系转换关系可得到炮弹运动加速度在弹体坐标系上的分量,即加速度计输出可表示为

(5)

式中:,,为三轴加速度计的输出;为加速度计的采样间隔。

(6)

式中:为地磁强度;、、为地磁强度在基准坐标系上的分量;为射向角;为磁倾角;为磁偏角。以南京地区为例,2010年南京地区的地磁要素=49 6671 nT,磁偏角=5′12″,磁倾角=48′24″。具体应用时可根据时期和地理区域加以选取。

运用第11节中坐标系转换关系,可得到地磁要素在弹体坐标系上的分量,即三轴磁力计的输出可表示为

(7)

式中:,,为地磁强度在弹体坐标系上的分量,即磁力计的三轴输出。

2 灵敏度分析方法

灵敏度分析可分为局部灵敏度分析和全局灵敏度分析2种方式。一般来说,局部灵敏度分析方法运算简单,计算效率高,但用其处理非线性模型时偏差较大,并且所有因素均需满足相同的误差范围。与之相反,全局灵敏度分析既考虑了各因素概率密度函数的分布影响,又满足了计算分析时所有因素不同的变动范围。

2.1 Sobol’法基本原理

Sobol’法的核心思想是方差分解。它通过把模型分解为单个参数及参数之间相互组合的函数,计算单个输入参数或输入参数集的方差对输出方差的贡献,进而来分析参数的重要性以及参数之间的交互效应。

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

(8)

对于求解函数1,2,…,可由多重积分及输出在不同条件下的期望表示,本文利用如下方程求解:

(9)

(10)

式中:为第个因子,～为矩阵中除去的所有因子,～(|)表示在因子不变而其他因子发生改变时输出的期望,类似地可求出高阶项。称()和1,2,…,(,,…,)的方差为和1,2,…,,则可得:

(11)

(12)

即

(13)

定义灵敏系数1,2,…,,表示为1,2,…,=1,2,…,。

为衡量输入参数作用的灵敏度,由式(13)可知:

式中:为参数的一阶灵敏度,它描述了随机因素“独自”对输出的总方差的贡献,越大,表明的不确定性对系统的输出影响越大;1,2,…,为参数的阶灵敏度,它描述了随机因素交互作用对输出总方差的贡献。

2.2 基于蒙特卡洛打靶的改进Sobol’法

基于蒙特卡洛打靶的改进Sobol’法的原理是利用仿真数据来近似式(11)中各项因式的方差,从而最终得到各阶灵敏度系数,一般应用步骤可描述如下。

①在输入因子设定范围内生成2个相互独立的矩阵×,×,其中为样本数。

③对式(11)用蒙特卡洛打靶法变换为

④一阶灵敏度可表示为

⑤同理,全局灵敏度可表示为

利用以上步骤,可以针对不同输入变量对系统输出的影响程度进行定量计算和分析研究。与传统Sobol’法相比,基于蒙特卡洛的改进Sobol’法在取值时,考虑了系统变量的具体散布特性(如选取弹道散布源)进行取值,更为接近所研究对象的实际状况。

3 仿真分析

3.1 仿真参数

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

本文以某大口径榴弹为研究对象,取初速900 m/s、射角45°进行仿真计算。在炮口初速的变量上假设有1%的不确定性范围;对于初始偏航角速率和初始俯仰角速率,在标称值0=0=0基础上假设了3 rad/s的不确定区间。炮弹的初始姿态对飞行轨迹和落点有直接的影响,在初始高低角和初始方向角的标称值上均假设了0.000 49 rad的不确定度。初始攻角在标称值0上假设了0.034 9 rad的不确定度。表1为所有状态变量的标称值及其不确定度。

表1 初始状态变量的标称值及不确定性区间

3.1.2 模型参数的不确定量

升力系数导数′、阻力系数和滚转阻尼力系数′对炮弹运动规律有很大的影响,因此将其选为不确定量。′、和′标称值见表2。根据外弹道工程经验,在其初值上均考虑了±5的不确定区间。其次,炮弹质量也具有不确定性,本文以45 kg炮弹为研究对象,根据经验在其初值上考虑了-0.25%～0.25%的不确定区间。模型参数的初值和其不确定区间汇总如表3所示。

表2 升力系数导数和阻力系数值

表3 模型参数的标称值及不确定性区间

根据表1和表3中所给出的初值和不确定性范围,使用Sobol’法对影响炮弹转动角速度的因素进行灵敏度分析。该分析是以一种时变的方法进行的,每0.04 s(即弹道方程组的积分步长)计算一次灵敏度,根据各变量灵敏度变化趋势,判断何种输入变量对炮弹飞行状态有较大影响。

3.2 仿真计算与结果分析

根据表1和表3中的标称条件及不确定性,使用改进的Sobol’法进行模拟。以下仿真在全参数和简化参数2种条件下进行。全参数仿真是指在考虑初速、初始偏航角速率、初始俯仰角速率、初始高低角、初始方向角、初始攻角、质量、升力系数导数、阻力系数、滚转阻尼力系数这10种参数扰动的情况下进行仿真;简化参数是指除去以上10种参数中对传感器输出影响较大的参数,对剩余参数的灵敏度值进行研究。

下面对陀螺仪、磁力计、加速度计三轴输出的一阶灵敏度值进行仿真。一阶灵敏度值是指单个输入参数自身对3种传感器输出的影响程度,一阶灵敏度值越大,表明该参数的不确定性对传感器的输出影响越大;而全局灵敏度值与一阶灵敏度值之差越接近于0,表明对应变量与其他输入变量之间的交互作用越小。

3.2.1 基于改进Sobol’法的全参数仿真的结果分析

对不确定性参数空间进行采样,假设样本均匀分布,样本大小为=20 000。仿真采用第2.2节所述步骤,并且由于使用了简化刚体弹道方程组,故每0.04 s计算一次一阶灵敏度,仿真结果如图1～图4所示。

图1 陀螺仪输出的一阶灵敏度值

图2 磁力计输出的一阶灵敏度值

图3 加速度计输出的一阶灵敏度值

图4 各传感器初始偏航角速率的全局与一阶灵敏度差值

由图1(a)和图1(b),炮弹初始偏航角速率、初始俯仰角速率对陀螺仪输出的一阶灵敏度值在0～1区间波动,由此推断初始偏航角速率和初始俯仰角速率是影响输出的主要因素,而初始攻角、初速和滚转阻尼力系数的灵敏度值在0～0.2间,表明初始攻角、初速和滚转阻尼力系数是影响输出的次要因素。图1(c)中,初速对的一阶灵敏度值稳定在1附近,这表明初速是影响输出值的主要因素。由图1(d)可知,炮弹初始偏航角速率和初始俯仰角速率是影响输出的主要因素;而初速在炮弹运动后期(出炮口约2 s后)对输出值产生小幅的影响。

根据图2(a)和图2(b),炮弹初始偏航角速率、初始俯仰角速率对磁力计输出的一阶灵敏度值在0～1区间波动,而初始攻角、初速和滚转阻尼力系数的灵敏度值在0～0.2区间波动,这表明,炮弹初始偏航角速率和初始俯仰角速率是影响输出的主要因素,而初始攻角、初速和滚转阻尼力系数对输出产生影响较小。由图2(c)、图2(d)可知,初速是影响、输出的主要因素;炮弹初始偏航角速率和初始俯仰角速率对、输出的影响随时间逐渐减小。

由图3(a)和图3(b)可知,阻力系数是影响输出的主要因素,而炮弹初始偏航角速率和初始俯仰角速率在炮弹运动初期对输出产生的影响较小。由图3(c)、图3(d)可知,炮弹初始偏航角速率和初始俯仰角速率是影响、输出的主要因素;初速在炮弹运动约2.8 s后对、输出产生较大影响。

为研究各输入变量间的交互作用,本文还计算了10种输入变量的全局灵敏度与一阶灵敏度的差值。计算结果表明,与其他9种输入变量相比,初始偏航角速率对3种传感器输出的一阶灵敏度值影响较大,故仅以初始偏航角速率为例,给出输入变量之间的交互作用计算结果,如图4所示。由图4可知,初始偏航角速率对3种传感器输出的全局灵敏度和一阶灵敏度差值在0～0.1间小幅波动,这表明初始偏航角速率的独立作用是影响炮弹飞行状态的主要部分,图4中输出的全局灵敏度系数与一阶灵敏度系数差值在9种传感器输出中最大,最大值为0.083 42,由此可知,在传感器输出的9个量中,元素交互耦合作用所占比例最大值约为8%。

通过以上仿真计算可以看出,初始偏航角速率和初始俯仰角速率对陀螺仪沿弹体径向的输出、加速度计沿弹体径向的输出和磁力计轴向输出的影响较大,初速对陀螺仪轴向输出和磁力计沿弹体径向的输出的影响较大,阻力系数对加速度计轴向输出的影响较大。后续研究可以参考以上结论,开展相关的试验设计及数据处理工作,譬如通过控制炮弹初速得到姿态角的准确测量结果。

本节结合灵敏度计算结果,辨识出了炮弹初速、初始偏航角速率、初始俯仰角速率、初始攻角和滚转阻尼力系数对炮弹飞行状态的影响程度。以上所述5项输入参数计算得到的一阶灵敏度数值较大,相比之下气动力系数等参数对炮弹飞行状态的影响可忽略不计。为将其余5项输入参数的影响程度进行分析排序,需要开展进一步的仿真计算。

3.2.2 结果验证

文献[16]提供了一个傅里叶幅值灵敏度测试法(Fourier amplitude sensitivity test,FAST),其计算复杂度较高,但经过验证是正确的。为对基于蒙特卡洛打靶的改进Sobol’法所得结果进行验证,利用文献[16]方法,在第3.1.1节所示条件下开展全局灵敏度分析,并将结果与第3.2.1节的结果共同列于表4。为便于比较,表4所列为0～2 s时段影响陀螺仪输出变量的一阶灵敏度变化趋势,表中的表达方式参考文献[1]。其中,符号“d”表示“影响程度极大”,“di”表示“影响程度随时间增加”,“ds”表示“影响程度极大且呈波动状变化”,“dl”表示“影响程度极小”,“ni”表示“无影响”。根据表4所示结果可知,采用两种不同计算方法得到的影响陀螺仪各轴输出的变量灵敏度分析结果完全一致,表明本文所用方法及相关分析结果是正确、可信的。

表4 影响陀螺仪输出变量的灵敏度(0～2 s)

3.2.3 简化参数仿真的结果分析

本节将考察3.2.1节中无法辨识出影响程度的5种输入,即初始高低角、初始方位角、炮弹质量、阻力系数和升力系数导数对炮弹飞行状态的影响。由图5(a)可知,阻力系数对的一阶灵敏度值稳定在0.8～0.9,质量的一阶灵敏度值稳定在0.1～0.2,则阻力系数和质量是影响输出的元素,其中阻力系数为主要影响因素。由图5(b)、图5(c)知,升力系数导数和阻力系数是影响和输出的主要因素。

图5 陀螺仪输出的一阶灵敏度值

由图6(a)可知,初始高低角是影响输出的最主要因素,同时升力系数导数和阻力系数也对输出值产生影响。由图6(b)知,初始方向角在炮弹飞行初期是影响输出的主要因素,后期影响输出的因素主要为升力系数导数和阻力系数。由图6(c)可知,在炮弹飞行期间,初始高低角、升力系数导数和阻力系数先后为影响输出的主要因素。

图6 磁力计输出的一阶灵敏度值

由图7(a)可以看出,阻力系数对的一阶灵敏度值稳定在0.97,其他4种元素灵敏度值几乎为0,则阻力系数是影响输出的主要因素。由图7(b)和图7(c)可知,升力系数导数和阻力系数是影响和输出的主要因素;升力系数导数和阻力系数交替对炮弹加速度产生影响。

图7 加速度计输出的一阶灵敏度值

4 结束语

本文以简化刚体弹道方程组为基础,应用基于蒙特卡洛打靶的改进Sobol’法,开展全局灵敏度分析,综合考虑初速、初始偏航角速率、初始俯仰角速率、初始高低角、初始方向角、质量、升力系数导数、阻力系数、初始攻角、滚转阻尼力系数等随机因素对3种传感器三轴输出的影响,定量计算得到了各因素对各传感器输出的灵敏度值。通过对仿真结果的分析,可得到如下结论:

①在全参数条件下,炮弹初始偏航角速率、初始俯仰角速率和初速是影响炮弹飞行状态的主要因素。初始偏航角速率和初始俯仰角速率在炮弹飞行期间交替对其飞行状态产生影响。

②从简化参数条件的分析结果看。阻力系数和升力系数导数是影响3种传感器输出的主要因素,质量对炮弹飞行状态的影响可忽略不计。

③炮弹的飞行状态受到各输入变量的独立影响和交叉耦合影响,其中前者是影响炮弹飞行的主要部分,交叉耦合影响在实际应用中可不考虑。

在后续的试验设计中,可以通过控制初始偏航角速率、初始俯仰角速率和初速的散布状态,得到更为精确的炮弹飞行状态参数。