基于Sobol和EFAST法的槍彈轉動角速度全局靈敏度仿真分析

魏偉， 張競文， 常思江， 王浩圣， 王建中

(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法開展靈敏度計算。

在后續的試驗設計中，可以通過控制初始攻角、初速的散布狀態，得到更為精確的槍彈飛行狀態參數。