炮彈飛行運動狀態(tài)的全局靈敏度分析

張競文,邵渤涵,常思江,謝 飛

(1.南京理工大學 能源與動力工程學院,江蘇 南京 210094;2.中國船舶集團有限公司系統(tǒng)工程研究院,北京 100094;3.遼沈工業(yè)集團有限公司,遼寧 沈陽 110045)

隨著信息化彈藥技術(shù)的不斷深入和發(fā)展,炮彈飛行狀態(tài)參數(shù)的獲取變得越來越重要和迫切。炮彈飛行狀態(tài)信息主要通過在彈上安裝各種嵌入式傳感器測量得到,而炮彈發(fā)射及飛行過程中存在各種隨機擾動,導致運動狀態(tài)具有不確定性。定量研究嵌入式傳感器數(shù)據(jù)對炮彈運動狀態(tài)不確定性的靈敏程度,有助于測量信息的數(shù)據(jù)處理及深入開展信息化彈藥技術(shù)的研究。所謂靈敏度分析,主要是采用一定的方法,定量計算出各輸入?yún)?shù)的靈敏度值,根據(jù)靈敏度值的大小篩選出對于待研究參數(shù)沒有影響或者影響較小的參數(shù),將其設(shè)為標稱值,而將主要研究精力或關(guān)注點集中于討論或設(shè)計對傳感器輸出影響較大的參數(shù)上。

近年來,靈敏度分析方法已逐步應(yīng)用于彈藥、火炮等領(lǐng)域,目前已有一些相關(guān)的應(yīng)用研究。文獻[1]建立了一個準線性變參數(shù)剛體彈道模型,據(jù)此在平射條件下研究了一些隨機擾動因素對彈載測量裝置輸出的靈敏度;文獻[6]引入局部靈敏度分析技術(shù),開展了炮彈的氣動和結(jié)構(gòu)參數(shù)辨識研究;文獻[7]提出一種解析目標級聯(lián)方法,對某運載火箭的輸入不確定性和模型不確定性進行了全局靈敏度分析;文獻[8]則針對混合不確定性的靈敏度分析問題,提出了一種基于輔助變量法的改進Sobol’靈敏度分析方法,并據(jù)此對某電磁軌道炮的射擊精度開展研究。文獻[9]針對自行火炮某些參數(shù)存在相關(guān)性等情況,提出了一種基于Rosenblatt轉(zhuǎn)換理論的全局靈敏度分析方法,實現(xiàn)了面向彈丸炮口狀態(tài)的某自行火炮結(jié)構(gòu)參數(shù)全局靈敏度分析。文獻[10]為評估、設(shè)計、提高自行火炮的射擊密集度,應(yīng)用Sobol’全局靈敏度分析方法定量分析了起始擾動對自行火炮射擊密集度的影響。

本文以某大口徑炮彈為研究對象,假設(shè)彈上安裝了陀螺儀、加速度計及磁力計3種傳感器,分別用于測量炮彈轉(zhuǎn)動角速度、加速度和姿態(tài)角。選取非線性簡化剛體彈道模型(稱為五自由度剛體彈道模型)為分析模型,綜合考慮了炮彈在發(fā)射及飛行過程中存在的初速跳動、起始擾動、彈丸質(zhì)量偏差、氣動系數(shù)偏差等隨機因素,利用基于蒙特卡洛打靶的改進Sobol’全局靈敏度分析方法,研究陀螺儀、加速度計及磁力計的三軸輸出對諸擾動因素的靈敏程度,為開展彈載傳感器測量信息與炮彈飛行狀態(tài)間關(guān)系的深入分析提供基礎(chǔ),與文獻[1]相比,本文較為全面地考慮了擾動因素(如起始攻角和滾轉(zhuǎn)阻尼力系數(shù)等),并擴展了分析條件(如射角和飛行時間),以期為深入的應(yīng)用研究提供參考。

1 動力學模型

炮彈姿態(tài)角的定義和建立炮彈彈道方程組所使用的坐標系及其轉(zhuǎn)換關(guān)系等密切相關(guān),姿態(tài)角的解算則依賴于不同坐標系的轉(zhuǎn)換關(guān)系。

1.1 坐標系及其相互轉(zhuǎn)換

本文用到的坐標系包括:地面坐標系、基準坐標系、彈軸坐標系、彈體坐標系。

①地面坐標系:坐標系原點通常選取在射出點,軸通常是彈道面與水平面交線,指向目標為正;軸沿垂線向上,軸與其他兩軸垂直并構(gòu)成右手坐標系。

②基準坐標系:坐標系原點位于彈箭質(zhì)心,軸沿水平線指向射擊方向,軸沿鉛錘向上,軸按右手法則確定為垂直于射擊面并指向右方,坐標軸隨質(zhì)心一起平動。

③彈軸坐標系:坐標系原點位于彈箭質(zhì)心,軸為彈軸,軸垂直于軸且指向上方,軸按右手法則垂直于平面且指向右方。

④彈體坐標系:坐標系原點位于彈箭質(zhì)心,軸為彈軸,軸和軸固連在彈體上并與彈體一同繞縱軸旋轉(zhuǎn)。

將彈軸坐標系向基準坐標系投影,即可得到彈軸坐標系和基準坐標系間的轉(zhuǎn)換關(guān)系:

(1)

式中:為彈軸方向角;為彈軸高低角。

彈體坐標系與彈軸坐標系的差別只在于坐標平面相對于坐標平面轉(zhuǎn)過一個自轉(zhuǎn)角,轉(zhuǎn)換關(guān)系為

(2)

式中:為彈體滾轉(zhuǎn)角。

本文選取三軸陀螺儀測量炮彈轉(zhuǎn)動角速度在彈體坐標系上的分量,三軸磁力計測量炮彈在彈體坐標系上的姿態(tài)角,三軸加速度計測量炮彈在彈體坐標系上的加速度。只要知道炮彈轉(zhuǎn)動角速度、地磁要素、彈箭加速度在彈軸坐標系上的分量,即可轉(zhuǎn)換求出其在彈體坐標系三軸上的分量。

1.2 外彈道及傳感器輸出模型

為了從3種傳感器輸出中分析各輸入因素不確定性的靈敏度,需使用炮彈運動模型,建立實測彈道與陀螺儀、加速度計、磁力計輸出之間的聯(lián)系。由于普通六自由度彈道方程的求解步長較小(一般為0.002 s),而五自由度彈道方程由于略去了快圓運動可使得積分步長大為增加(如可取為0.04 s),所以計算時間幾乎只有普通六自由度彈道方程的1/10?？紤]到本文以分析炮彈飛行狀態(tài)為主,且需計算數(shù)十萬次彈道以得到各因素的靈敏度值,為在確保計算精度的前提下提高計算效率,選用如下五自由度剛體彈道模型:

(3)

式中:為彈箭質(zhì)心的速度;和分別為彈箭速度高低角和方向角;和分別為彈軸高低角和方位角;為彈體滾轉(zhuǎn)角;,,為彈箭繞質(zhì)心轉(zhuǎn)動的角速度;,,為彈箭質(zhì)心在地面坐標系的位置分量;為彈箭質(zhì)量;,,為炮彈所受外力在彈道坐標系三軸上的分量;為軸向轉(zhuǎn)動慣量;,,為外力矩在彈軸坐標系三軸上的分量;軸力矩的具體表達式可參見文獻[11]。

根據(jù)五自由度外彈道模型可求出炮彈轉(zhuǎn)動角速度在彈軸坐標系上的分量、和,再運用第11節(jié)中坐標系轉(zhuǎn)換關(guān)系可得到炮彈轉(zhuǎn)動角速度在彈體坐標系上的分量,則三軸陀螺儀的輸出可表示為

(4)

式中:,,為陀螺儀三軸上的輸出。

由五自由度外彈道模型可求出彈箭質(zhì)心的速度,運用坐標系轉(zhuǎn)換關(guān)系可得到炮彈運動加速度在彈體坐標系上的分量,即加速度計輸出可表示為

(5)

式中:,,為三軸加速度計的輸出;為加速度計的采樣間隔。

(6)

式中:為地磁強度;、、為地磁強度在基準坐標系上的分量;為射向角;為磁傾角;為磁偏角。以南京地區(qū)為例,2010年南京地區(qū)的地磁要素=49 6671 nT,磁偏角=5′12″,磁傾角=48′24″。具體應(yīng)用時可根據(jù)時期和地理區(qū)域加以選取。

運用第11節(jié)中坐標系轉(zhuǎn)換關(guān)系,可得到地磁要素在彈體坐標系上的分量,即三軸磁力計的輸出可表示為

(7)

式中:,,為地磁強度在彈體坐標系上的分量,即磁力計的三軸輸出。

2 靈敏度分析方法

靈敏度分析可分為局部靈敏度分析和全局靈敏度分析2種方式。一般來說,局部靈敏度分析方法運算簡單,計算效率高,但用其處理非線性模型時偏差較大,并且所有因素均需滿足相同的誤差范圍。與之相反,全局靈敏度分析既考慮了各因素概率密度函數(shù)的分布影響,又滿足了計算分析時所有因素不同的變動范圍。

2.1 Sobol’法基本原理

Sobol’法的核心思想是方差分解。它通過把模型分解為單個參數(shù)及參數(shù)之間相互組合的函數(shù),計算單個輸入?yún)?shù)或輸入?yún)?shù)集的方差對輸出方差的貢獻,進而來分析參數(shù)的重要性以及參數(shù)之間的交互效應(yīng)。

假設(shè)可積函數(shù)()的定義域為維的單元體,=(…),Sobol’法的中心思想是將函數(shù)()分解為子項之和。

(8)

對于求解函數(shù)1,2,…,可由多重積分及輸出在不同條件下的期望表示,本文利用如下方程求解:

(9)

(10)

式中:為第個因子,～為矩陣中除去的所有因子,～(|)表示在因子不變而其他因子發(fā)生改變時輸出的期望,類似地可求出高階項。稱()和1,2,…,(,,…,)的方差為和1,2,…,,則可得:

(11)

(12)

即

(13)

定義靈敏系數(shù)1,2,…,,表示為1,2,…,=1,2,…,。

為衡量輸入?yún)?shù)作用的靈敏度,由式(13)可知:

式中:為參數(shù)的一階靈敏度,它描述了隨機因素“獨自”對輸出的總方差的貢獻,越大,表明的不確定性對系統(tǒng)的輸出影響越大;1,2,…,為參數(shù)的階靈敏度,它描述了隨機因素交互作用對輸出總方差的貢獻。

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

基于蒙特卡洛打靶的改進Sobol’法的原理是利用仿真數(shù)據(jù)來近似式(11)中各項因式的方差,從而最終得到各階靈敏度系數(shù),一般應(yīng)用步驟可描述如下。

①在輸入因子設(shè)定范圍內(nèi)生成2個相互獨立的矩陣×,×,其中為樣本數(shù)。

③對式(11)用蒙特卡洛打靶法變換為

④一階靈敏度可表示為

⑤同理,全局靈敏度可表示為

利用以上步驟,可以針對不同輸入變量對系統(tǒng)輸出的影響程度進行定量計算和分析研究。與傳統(tǒng)Sobol’法相比,基于蒙特卡洛的改進Sobol’法在取值時,考慮了系統(tǒng)變量的具體散布特性(如選取彈道散布源)進行取值,更為接近所研究對象的實際狀況。

3 仿真分析

3.1 仿真參數(shù)

3.1.1 初始發(fā)射狀態(tài)參數(shù)的不確定量

本文以某大口徑榴彈為研究對象,取初速900 m/s、射角45°進行仿真計算。在炮口初速的變量上假設(shè)有1%的不確定性范圍;對于初始偏航角速率和初始俯仰角速率,在標稱值0=0=0基礎(chǔ)上假設(shè)了3 rad/s的不確定區(qū)間。炮彈的初始姿態(tài)對飛行軌跡和落點有直接的影響,在初始高低角和初始方向角的標稱值上均假設(shè)了0.000 49 rad的不確定度。初始攻角在標稱值0上假設(shè)了0.034 9 rad的不確定度。表1為所有狀態(tài)變量的標稱值及其不確定度。

表1 初始狀態(tài)變量的標稱值及不確定性區(qū)間

3.1.2 模型參數(shù)的不確定量

升力系數(shù)導數(shù)′、阻力系數(shù)和滾轉(zhuǎn)阻尼力系數(shù)′對炮彈運動規(guī)律有很大的影響,因此將其選為不確定量?！?span id="j5i0abt0b" class="subscript">、和′標稱值見表2。根據(jù)外彈道工程經(jīng)驗,在其初值上均考慮了±5的不確定區(qū)間。其次,炮彈質(zhì)量也具有不確定性,本文以45 kg炮彈為研究對象,根據(jù)經(jīng)驗在其初值上考慮了-0.25%～0.25%的不確定區(qū)間。模型參數(shù)的初值和其不確定區(qū)間匯總?cè)绫?所示。

表2 升力系數(shù)導數(shù)和阻力系數(shù)值

表3 模型參數(shù)的標稱值及不確定性區(qū)間

根據(jù)表1和表3中所給出的初值和不確定性范圍,使用Sobol’法對影響炮彈轉(zhuǎn)動角速度的因素進行靈敏度分析。該分析是以一種時變的方法進行的,每0.04 s(即彈道方程組的積分步長)計算一次靈敏度,根據(jù)各變量靈敏度變化趨勢,判斷何種輸入變量對炮彈飛行狀態(tài)有較大影響。

3.2 仿真計算與結(jié)果分析

根據(jù)表1和表3中的標稱條件及不確定性,使用改進的Sobol’法進行模擬。以下仿真在全參數(shù)和簡化參數(shù)2種條件下進行。全參數(shù)仿真是指在考慮初速、初始偏航角速率、初始俯仰角速率、初始高低角、初始方向角、初始攻角、質(zhì)量、升力系數(shù)導數(shù)、阻力系數(shù)、滾轉(zhuǎn)阻尼力系數(shù)這10種參數(shù)擾動的情況下進行仿真;簡化參數(shù)是指除去以上10種參數(shù)中對傳感器輸出影響較大的參數(shù),對剩余參數(shù)的靈敏度值進行研究。

下面對陀螺儀、磁力計、加速度計三軸輸出的一階靈敏度值進行仿真。一階靈敏度值是指單個輸入?yún)?shù)自身對3種傳感器輸出的影響程度,一階靈敏度值越大,表明該參數(shù)的不確定性對傳感器的輸出影響越大;而全局靈敏度值與一階靈敏度值之差越接近于0,表明對應(yīng)變量與其他輸入變量之間的交互作用越小。

3.2.1 基于改進Sobol’法的全參數(shù)仿真的結(jié)果分析

對不確定性參數(shù)空間進行采樣,假設(shè)樣本均勻分布,樣本大小為=20 000。仿真采用第2.2節(jié)所述步驟,并且由于使用了簡化剛體彈道方程組,故每0.04 s計算一次一階靈敏度,仿真結(jié)果如圖1～圖4所示。

圖1 陀螺儀輸出的一階靈敏度值

圖2 磁力計輸出的一階靈敏度值

圖3 加速度計輸出的一階靈敏度值

圖4 各傳感器初始偏航角速率的全局與一階靈敏度差值

由圖1(a)和圖1(b),炮彈初始偏航角速率、初始俯仰角速率對陀螺儀輸出的一階靈敏度值在0～1區(qū)間波動,由此推斷初始偏航角速率和初始俯仰角速率是影響輸出的主要因素,而初始攻角、初速和滾轉(zhuǎn)阻尼力系數(shù)的靈敏度值在0～0.2間,表明初始攻角、初速和滾轉(zhuǎn)阻尼力系數(shù)是影響輸出的次要因素。圖1(c)中,初速對的一階靈敏度值穩(wěn)定在1附近,這表明初速是影響輸出值的主要因素。由圖1(d)可知,炮彈初始偏航角速率和初始俯仰角速率是影響輸出的主要因素;而初速在炮彈運動后期(出炮口約2 s后)對輸出值產(chǎn)生小幅的影響。

根據(jù)圖2(a)和圖2(b),炮彈初始偏航角速率、初始俯仰角速率對磁力計輸出的一階靈敏度值在0～1區(qū)間波動,而初始攻角、初速和滾轉(zhuǎn)阻尼力系數(shù)的靈敏度值在0～0.2區(qū)間波動,這表明,炮彈初始偏航角速率和初始俯仰角速率是影響輸出的主要因素,而初始攻角、初速和滾轉(zhuǎn)阻尼力系數(shù)對輸出產(chǎn)生影響較小。由圖2(c)、圖2(d)可知,初速是影響、輸出的主要因素;炮彈初始偏航角速率和初始俯仰角速率對、輸出的影響隨時間逐漸減小。

由圖3(a)和圖3(b)可知,阻力系數(shù)是影響輸出的主要因素,而炮彈初始偏航角速率和初始俯仰角速率在炮彈運動初期對輸出產(chǎn)生的影響較小。由圖3(c)、圖3(d)可知,炮彈初始偏航角速率和初始俯仰角速率是影響、輸出的主要因素;初速在炮彈運動約2.8 s后對、輸出產(chǎn)生較大影響。

為研究各輸入變量間的交互作用,本文還計算了10種輸入變量的全局靈敏度與一階靈敏度的差值。計算結(jié)果表明,與其他9種輸入變量相比,初始偏航角速率對3種傳感器輸出的一階靈敏度值影響較大,故僅以初始偏航角速率為例,給出輸入變量之間的交互作用計算結(jié)果,如圖4所示。由圖4可知,初始偏航角速率對3種傳感器輸出的全局靈敏度和一階靈敏度差值在0～0.1間小幅波動,這表明初始偏航角速率的獨立作用是影響炮彈飛行狀態(tài)的主要部分,圖4中輸出的全局靈敏度系數(shù)與一階靈敏度系數(shù)差值在9種傳感器輸出中最大,最大值為0.083 42,由此可知,在傳感器輸出的9個量中,元素交互耦合作用所占比例最大值約為8%。

通過以上仿真計算可以看出,初始偏航角速率和初始俯仰角速率對陀螺儀沿彈體徑向的輸出、加速度計沿彈體徑向的輸出和磁力計軸向輸出的影響較大,初速對陀螺儀軸向輸出和磁力計沿彈體徑向的輸出的影響較大,阻力系數(shù)對加速度計軸向輸出的影響較大。后續(xù)研究可以參考以上結(jié)論,開展相關(guān)的試驗設(shè)計及數(shù)據(jù)處理工作,譬如通過控制炮彈初速得到姿態(tài)角的準確測量結(jié)果。

本節(jié)結(jié)合靈敏度計算結(jié)果,辨識出了炮彈初速、初始偏航角速率、初始俯仰角速率、初始攻角和滾轉(zhuǎn)阻尼力系數(shù)對炮彈飛行狀態(tài)的影響程度。以上所述5項輸入?yún)?shù)計算得到的一階靈敏度數(shù)值較大,相比之下氣動力系數(shù)等參數(shù)對炮彈飛行狀態(tài)的影響可忽略不計。為將其余5項輸入?yún)?shù)的影響程度進行分析排序,需要開展進一步的仿真計算。

3.2.2 結(jié)果驗證

文獻[16]提供了一個傅里葉幅值靈敏度測試法(Fourier amplitude sensitivity test,FAST),其計算復雜度較高,但經(jīng)過驗證是正確的。為對基于蒙特卡洛打靶的改進Sobol’法所得結(jié)果進行驗證,利用文獻[16]方法,在第3.1.1節(jié)所示條件下開展全局靈敏度分析,并將結(jié)果與第3.2.1節(jié)的結(jié)果共同列于表4。為便于比較,表4所列為0～2 s時段影響陀螺儀輸出變量的一階靈敏度變化趨勢,表中的表達方式參考文獻[1]。其中,符號“d”表示“影響程度極大”,“di”表示“影響程度隨時間增加”,“ds”表示“影響程度極大且呈波動狀變化”,“dl”表示“影響程度極小”,“ni”表示“無影響”。根據(jù)表4所示結(jié)果可知,采用兩種不同計算方法得到的影響陀螺儀各軸輸出的變量靈敏度分析結(jié)果完全一致,表明本文所用方法及相關(guān)分析結(jié)果是正確、可信的。

表4 影響陀螺儀輸出變量的靈敏度(0～2 s)

3.2.3 簡化參數(shù)仿真的結(jié)果分析

本節(jié)將考察3.2.1節(jié)中無法辨識出影響程度的5種輸入,即初始高低角、初始方位角、炮彈質(zhì)量、阻力系數(shù)和升力系數(shù)導數(shù)對炮彈飛行狀態(tài)的影響。由圖5(a)可知,阻力系數(shù)對的一階靈敏度值穩(wěn)定在0.8～0.9,質(zhì)量的一階靈敏度值穩(wěn)定在0.1～0.2,則阻力系數(shù)和質(zhì)量是影響輸出的元素,其中阻力系數(shù)為主要影響因素。由圖5(b)、圖5(c)知,升力系數(shù)導數(shù)和阻力系數(shù)是影響和輸出的主要因素。

圖5 陀螺儀輸出的一階靈敏度值

由圖6(a)可知,初始高低角是影響輸出的最主要因素,同時升力系數(shù)導數(shù)和阻力系數(shù)也對輸出值產(chǎn)生影響。由圖6(b)知,初始方向角在炮彈飛行初期是影響輸出的主要因素,后期影響輸出的因素主要為升力系數(shù)導數(shù)和阻力系數(shù)。由圖6(c)可知,在炮彈飛行期間,初始高低角、升力系數(shù)導數(shù)和阻力系數(shù)先后為影響輸出的主要因素。

圖6 磁力計輸出的一階靈敏度值

由圖7(a)可以看出,阻力系數(shù)對的一階靈敏度值穩(wěn)定在0.97,其他4種元素靈敏度值幾乎為0,則阻力系數(shù)是影響輸出的主要因素。由圖7(b)和圖7(c)可知,升力系數(shù)導數(shù)和阻力系數(shù)是影響和輸出的主要因素;升力系數(shù)導數(shù)和阻力系數(shù)交替對炮彈加速度產(chǎn)生影響。

圖7 加速度計輸出的一階靈敏度值

4 結(jié)束語

本文以簡化剛體彈道方程組為基礎(chǔ),應(yīng)用基于蒙特卡洛打靶的改進Sobol’法,開展全局靈敏度分析,綜合考慮初速、初始偏航角速率、初始俯仰角速率、初始高低角、初始方向角、質(zhì)量、升力系數(shù)導數(shù)、阻力系數(shù)、初始攻角、滾轉(zhuǎn)阻尼力系數(shù)等隨機因素對3種傳感器三軸輸出的影響,定量計算得到了各因素對各傳感器輸出的靈敏度值。通過對仿真結(jié)果的分析,可得到如下結(jié)論:

①在全參數(shù)條件下,炮彈初始偏航角速率、初始俯仰角速率和初速是影響炮彈飛行狀態(tài)的主要因素。初始偏航角速率和初始俯仰角速率在炮彈飛行期間交替對其飛行狀態(tài)產(chǎn)生影響。

②從簡化參數(shù)條件的分析結(jié)果看。阻力系數(shù)和升力系數(shù)導數(shù)是影響3種傳感器輸出的主要因素,質(zhì)量對炮彈飛行狀態(tài)的影響可忽略不計。

③炮彈的飛行狀態(tài)受到各輸入變量的獨立影響和交叉耦合影響,其中前者是影響炮彈飛行的主要部分,交叉耦合影響在實際應(yīng)用中可不考慮。

在后續(xù)的試驗設(shè)計中,可以通過控制初始偏航角速率、初始俯仰角速率和初速的散布狀態(tài),得到更為精確的炮彈飛行狀態(tài)參數(shù)。