面向彈丸炮口狀態的自行火炮結構參數全局靈敏度

羅中峰， 管小榮， 徐誠

(南京理工大學 機械工程學院， 江蘇 南京 210094)

0 引言

研究自行火炮結構參數在給定分布條件下對彈丸炮口狀態的敏感程度，有助于進一步精確降低彈丸炮口狀態波動，提高自行火炮射擊精度。進行自行火炮結構參數對彈丸炮口狀態的敏感程度分析，首先需要建立自行火炮發射動力學模型。陳光宋[1]基于剛柔耦合和彈炮耦合約束方程建立了某自行火炮彈丸和身管耦合動力學理論模型。羅中峰等[2]基于有限段法建立了某自行火炮彈丸/身管相互作用模型。王曉鋒等[3]通過凱恩方程和Huston方法建立了某自行火炮多體動力學模型。劉飛飛等[4]基于多體系統傳遞矩陣法和總傳遞方程自動推導定理，建立了某自行火炮行進間射擊動力學方程。景鵬淵等[5]基于有限元法建立了某火炮發射動力學模型。上述建模方式存在以下問題：只進行彈丸/身管相互作用過程的建模，不能考慮自行火炮車體振動對彈丸膛內運動的影響；完全根據多剛體動力學理論建模，不能考慮彈丸的膛內運動；完全利用有限元理論建立自行火炮發射動力學模型，計算效率太低。

為了在面向彈丸炮口狀態的自行火炮結構參數敏感程度分析過程中，充分考慮自行火炮參數的實際分布情況，需要使用全局靈敏度分析方法進行相關分析。Saltelli[6]經過研究發現:當相關變量按不同順序產生時，該相關變量對目標函數方差的貢獻是不同的。Saltelli等[[7]利用雙層蒙特卡洛技術對單個相關變量的1階和全敏感系數進行了估計。Kucherenko等[8]通過擴展敏感系數定義，推導了單個相關變量的1階和全敏感系數的估計方法。Xu等[9]建議將相關變量的1階敏感系數分解為相關和不相關兩部分；同時還提供了一種當目標函數與輸入變量為線性關系時，相關變量1階敏感系數相關和不相關部分的計算方法。Li等[10]提供了一種當目標函數與輸入變量為非線性關系時，相關變量1階敏感系數相關和不相關部分的計算方法。Mara等[11]利用條件樣本對相關變量組的敏感系數進行了估計。Tarantola等[12]利用傅里葉振幅敏度測試，對單個相關變量的敏感系數進行了估計。Decarlo等[13]根據貝葉斯校準數據，對單個相關變量的敏感系數進行了估計。Zhou等[14]基于樣本間距離，對單個相關變量的敏感系數進行了估計。綜上所述，當前全局靈敏度分析方法存在估計效率低和處理復雜等問題。

為了克服上述問題，本文首先建立能夠同時考慮彈丸膛內運動和車體振動的某自行火炮整體發射動力學模型；然后給出一種基于Rosenblatt轉換理論的全局靈敏度估計方法；最后基于上述模型和方法，進行面向彈丸炮口狀態的某自行火炮結構參數全局靈敏度分析。

1 自行火炮整體發射動力學模型

1.1 彈丸炮口狀態定義

為了獲得完整的彈丸炮口位置和姿態，參照6自由度彈丸外彈道方程[15]的輸入要求，構建彈丸炮口位置和姿態表征方法，如圖1所示。

圖1 彈丸炮口位置和姿態Fig.1 Position and attitude of projectile at muzzle

由于身管和彈丸間相互作用的結束，所以彈丸在炮口時(彈丸尾部到達炮口)可以認為是剛體。為了便于區別，作如下定義：圖1所涉及的彈丸質心稱為炮口質心，其為剛體彈丸質心；對應的彈軸稱為炮口彈軸，其為剛體彈丸彈軸；彈丸在炮口時，炮口質心速度和水平面的夾角與火炮高低射角之差為彈丸炮口高低偏角， 炮口質心速度與射擊面的夾角為彈丸炮口方向偏角，炮口彈軸和水平面的夾角與火炮高低射角之差為彈丸炮口高低擺角, 炮口彈軸與射擊面的夾角為彈丸炮口方向擺角，繞炮口彈軸的彈丸自轉角為彈丸炮口自轉角。

1.2 自行火炮整體發射動力學模型

建模前，作如下假設：不考慮供輸彈動作對車體的激勵作用；自行火炮在水平和靜止狀態下進行射擊；在射擊過程中，所有車輪處于制動狀態，懸架處于開鎖狀態；地面為松軟土壤地面。

為了獲得較準確的彈丸炮口狀態，建立了某自行火炮整體發射動力學模型。該模型由彈丸/身管相互作用模型和自行火炮發射動力學模型組成。彈丸/身管相互作用模型用于考慮彈丸膛內運動；自行火炮發射動力學模型用于考慮發射過程中的車體振動。彈丸/身管相互作用模型首先計算一個積分步，并在結果收斂條件下將獲得的身管與高低機和搖架與耳軸連接處的受力、位移、速度和加速度輸入到自行火炮發射動力學模型；其次自行火炮發射動力學模型以上述輸入為基礎計算一個積分步，并在結果收斂條件下將獲得的身管與高低機和搖架與耳軸連接處的受力和位移、速度，加速度傳遞給彈丸/身管相互作用模型；最后彈丸/身管相互作用模型再次計算。如此循環，直到彈丸膛內運動結束。彈丸膛內運動結束，只運行自行火炮發射動力學模型。彈丸炮口狀態由彈丸/身管相互作用模型計算得到。

1.2.1 彈丸/身管相互作用模型

1.2.1.1 彈丸膛內狀態定義

為了獲得彈丸炮口位置和姿態，參照圖1構建彈丸膛內位置和姿態的表征方法，如圖2所示。

圖2 彈丸膛內位置和姿態Fig.2 Position and attitude of projectile in bore

由于身管和彈丸間的相互作用，彈丸在身管中運動時存在彈性變形，其質心和彈軸是不確定的。為了計算需要，作如下定義：圖2所涉及的彈丸質心稱為膛內質心，它是彈丸為剛體時其質心所在點(該點隨著彈丸的變形會發生位置變化)；相應的彈軸稱為膛內彈軸，其經過膛內質心，垂直于該質心所在彈丸橫截面；彈丸在身管中運動時，膛內質心速度和水平面的夾角與火炮高低射角之差為彈丸膛內高低偏角，膛內質心速度與射擊面的夾角為彈丸膛內方向偏角，膛內彈軸和水平面的夾角與火炮高低射角之差為彈丸膛內高低擺角，膛內彈軸與射擊面的夾角為彈丸膛內方向擺角，繞膛內彈軸的彈丸自轉角為彈丸膛內自轉角。

當彈丸到達炮口時，圖2所示的彈丸膛內位置和姿態就會和圖1所示的彈丸炮口位置和姿態對應相等。原因是：當彈丸到達炮口時，由于身管和彈丸間相互作用的結束，彈丸由柔性體變成剛體，彈丸膛內質心和膛內彈軸就分別會與彈丸炮口質心和炮口彈軸重合?；谏鲜鲈?，在計算過程中通過持續記錄膛內質心和膛內彈軸位置變化情況，就可以計算出彈丸炮口的位置和姿態。

1.2.1.2 彈丸/身管相互作用模型

為了考慮彈丸的膛內運動，根據有限段思想建立了如圖3[2]和圖4[2]所示彈丸/身管相互作用模型。彈丸/身管相互作用模型的建模詳情見文獻[2]。

圖3 身管半剖有限段模型Fig.3 Half sectional mesh model of barrel

圖4 彈丸半剖有限段模型Fig.4 Half sectional mesh model of projectile

1.2.2 自行火炮發射動力學模型

為了考慮發射過程中自行火炮車體振動，基于Adams軟件建立了如圖5所示某自行火炮發射動力學模型。

圖5 某自行火炮發射動力學模型Fig.5 Launch dynamics model of a self-propelled gun

1.2.2.1 模型拓撲結構

圖5所示自行火炮發射動力學模型的拓撲關系如圖6所示。在圖6中，hl(l=1,2,…,15)表示部件與部件之間的連接方式和載荷作用關系。其中：輪胎與地面間通過摩擦副進行連接，二者之間的載荷包括水平摩擦力和垂直壓力；車輪與車身間通過彈簧和阻尼進行連接，二者之間的載荷包括彈簧力和阻尼力；炮塔與車身間通過襯套進行連接，二者之間的載荷包括旋轉扭矩和接觸壓力；搖架與炮塔間通過旋轉副進行連接，二者之間的載荷包括旋轉扭矩和后坐力；反后坐裝置中后坐部分與搖架間通過滑動副進行連接，反后坐裝置中制退機和復進機部分與搖架固聯；身管與反后坐裝置中制退機和復進機固聯在一起；身管與炮口制退器固聯在一起；身管和高低機通過套箍連接在一起，身管和套箍固聯在一起，套箍和高低機通過旋轉副連接；高低機與炮塔間直接通過旋轉副連接。

圖6 某自行火炮發射動力學模型拓撲圖Fig.6 Launch dynamics model topography of a self-propelled gun

1.2.2.2 碰撞力計算

為了考慮身管與高低機連接處因間隙產生的碰撞，上述自行火炮發射動力學模型根據(1)式所示Lankarani-Nikravesh模型[16]進行相關碰撞力的計算。

(1)

1.2.2.3 地面與輪胎間的接觸力計算

為了準確計算地面與輪胎間的接觸力，上述自行火炮發射動力學模型根據(2)式所示Wong-Reece模型[17]進行相關接觸力的計算。

(2)

1.3 模型驗證

在公開文獻中尚沒有某自行火炮彈丸炮口狀態的測量數據。根據現有數據，通過特定條件下某自行火炮底盤和炮口垂直振動位移的測量值和計算值的對比，實現對上述自行火炮整體發射動力學模型的驗證。在表1所示條件下，某自行火炮底盤和炮口垂直振動位移的測量值和計算值分別如圖7和圖8所示(圖7所示測量值源于文獻[18])。由圖7和圖8可知，在表1所示條件下，某自行火炮底盤和炮口垂直振動位移的測量值與其對應的計算值相對誤差較小。所以上述某自行火炮整體發射動力學模型是可信的。

表1 參數列表

圖7 底盤垂直振動位移對比圖Fig.7 Measured and calculated vertical vibration displacements of chassis

圖8 炮口垂直振動位移對比圖Fig.8 Measured and calculated vertical vibration displacements of muzzle

2 基于Rosenblatt轉換理論的全局靈敏度分析方法

為了在自行火炮結構參數對彈丸炮口狀態的敏感程度分析過程中，充分考慮自行火炮參數的實際分布情況，并克服當前全局靈敏度分析方法估計效率低和處理復雜等問題，本節以Sobol敏感系數定義為基礎發展一種全局靈敏度分析方法。

2.1 敏感系數定義

本節所考慮的目標函數如(3)式所示。

Z=F(x1,…,xi-1,xi,xi+1,…,xn),

(3)

式中：x=[x1,…,xi-1,xi,xi+1,…,xn]是一個連續實隨機向量，i=1,2,…,n，其聯合概率密度函數為h(x).w是包含前i個輸入變量的變量組，w=[x1,x2,…,xi]；w-是包含除前i個輸出變量以外所有輸入變量的變量組，w-=[xi+1,xi+2,…,xn]。

假設(3)式所示目標函數擁有有限方差。當(3)式所示目標函數擁有有限方差時，其方差V(Z)可以按(4)式[11]進行分解。

V(Z)=V[E[F(x)|w]]+

E[V[F(x)|(w|w-)]]=

V[E[F(x)|w-]]+

E[V[F(x)|(w-|w)]]，

(4)

式中：E(·)為函數期望；F(x)|w為在給定變量組w條件下F(x)的取值；w|w-為在給定變量組w-的條件下變量組w的取值。

為了衡量變量組w對目標函數方差的影響，有(5)式[11]和(6)式[11]所示定義。(5)式所示比率為變量組w的1階敏感系數。

(5)

式中：hw-|w(w-)是變量組w-在給定w條件下的概率密度函數;hw(w)是變量組w的概率密度函數；R為實數集，R(i)表示本次積分對i個實變量進行積分，R(n-i)表示本次積分對n-i個實變量進行積分。該敏感系數用于記錄由變量組w引發，變量組w和變量xi+1,xi+2,…,xn單獨，對目標函數方差V(Z)的貢獻率。

(6)式所示比率為變量組w的全敏感系數。

(6)

該敏感系數用于記錄由w引發，變量組w和變量xi+1,xi+2,…,xn的單獨項和它們的交叉項，對目標函數方差V(Z)的貢獻率。

當x1,x2,…,xn都是獨立變量時，(4)式可以轉換成(7)式[6]，(5)式可以轉換成(8)式[6]，(6)式可以轉換成(9)式[6]。

V(Z)=V[E[F(x)|w]]+E[V[F(x)|w]]=
V[E[F(x)|w-]]+E[V[F(x)|w-]]，

(7)

(8)

式中：hxk(xk)為變量xk的邊緣概率密度函數，k=1，2，…，n.

(9)

在后續段落中，視x1,x2,…,xn為相關變量。

2.2 相關變量的產生

2.2.1 相關變量的產生

根據Rosenblatt轉換[19]，有(10)式成立。

(10)

式中：ui是均勻分布于[0,1]的獨立實隨機變量；Hxi|u(xi)是在給定變量集u?[x1,…,xi-1,xi+1,…,xn]的條件下，隨機變量xi的分布函數；Hxi|u是單調增函數。

對(10)式所含各式進行求反，有(11)式成立。

(11)

同樣根據Rosenblatt轉換規則[19]，有(12)式成立。

(12)

式中：Gyi(yi)是連續獨立實隨機變量yi的分布函數。

(12)式代入到(11)式，有(13)式成立。

(13)

由(13)式可知：連續實向量x=[x1,…,xi-1,xi,xi+1,…,xn]可以用任意連續獨立實向量y=[y1,…,yi-1,yi,yi+1,…,yn]表示。當向量x中變量的順序發生變化時，變量xi的條件分布函數會發生變化，(13)式所示對應關系也會發生變化。

調整向量x中各個隨機變量的順序，得到x=[xi,x1,…,xi-1,xi+1,…,xn]。參照(13)式，有(14)式成立。

(14)

2.2.2 相關變量產生性質分析

當參照(14)式產生相關變量時，根據概率密度函數的性質[11]，相關變量(組)的分布情況如(15)式所示。

在條件概率密度函數的條件[11]下，向量x的聯合概率密度函數有(16)式所示形式。

對比(15)式中各個變量(組)的概率密度函數和(16)式中向量x的聯合概率密度函數形式，可知： (14)式是參照(16)式所示聯合概率密度函數的一種形式(hw(w)=hxi(xi)hx1|xi(x1)hx2|x1,xi(x2)…hxi-1|x1,x2,…,xi-2,xi(xi-1);hw-|w(w-)=hxi+1|w(xi+1)hxi+2|w,xi+1(xi+2)…hxn|w,xi+1,xi+2,…,xn-1(xn))產生的相關變量。在2.2.1節所示方法的幫助下，通過調整相關變量的產生順序，(14)式也可以按照(16)式所示聯合概率密度函數的其他形式產生相關變量。比如，參照hw(w)=hxi(xi)hx1|xi(x1)hx2|x1,xi(x2)…xi-1|x1,x2,…,xi-2,xi(xi-1)，hw-|w(w-)=hxi+2|w(xi+2)hxi+1|w,xi+2(xi+1)hxi+3|w,xi+1,xi+2(xi+3)…hxn|w,xi+1,xi+2,…,xn-1(xn)，有(17)式成立。

參照(16)式所示聯合概率密度函數的任意一種形式產生w和w-時，均有x～h(x)、w～hw(w)和w-～hw-|w(w-)成立。參照(16)式所示聯合概率密度函數的任意一種形式產生w和w-，根據(5)式和(6)式計算w的1階敏感系數和全敏感系數時，h(x)、hw(w)、hw-|w(w-)、目標函數和積分域都是一致的。所以有如下結論：不管參照(16)式所示聯合概率密度函數的何種形式產生w和w-，并以此為基礎根據(5)式和(6)式計算w的1階敏感系數和全敏感系數，獲得w的1階敏感系數和全敏感系數是唯一的，都是w的1階敏感系數和全敏感系數。

綜上所述，可以參照(16)式所示x聯合概率密度函數的一種形式獲得w和w-，并以此為基礎根據(5)式和(6)式計算w的1階敏感系數和全敏感系數。下面將以(14)式基礎，進行w的1階敏感系數和全敏感系數的估計研究。

(15)

式中：hxi|u(xi)為在給定u?[x1,…,xi-1,xi+1,…,xn]的條件下，相關變量xi的概率密度函數。

h(x)=hw(w)hw-|w(w-),

(16)

其中

hw(w)=hxi(xi)hx1|xi(x1)hx2|x1,xi(x2)…·
hxi-1|x1,x2,…,xi-2,xi(xi-1)=
hx1(x1)hx2|x1(x2)…hxi|x1,x2,…,xi-1(xi-1)…=
hx2(x2)hx1|x2(x1)hx3|x1,x2(x3)…hxi|x1,x2,…,xi-1(xi-1)，
hw-|w(w-)=hxi+1|w(xi+1)hxi+2|w,xi+1(xi+2)…
hxn|w,xi+1,xi+2,…,xn-1,(xn)=hxi+2|w(xi+2)·
hxi+1|w,xi+2(xi+1)hxi+3|w,xi+1,xi+2(xi+3)…
hxn|w,xi+1,xi+2,…,xn-1(xn)…=
hxn|w(xn)hxi+1|w,xn(xi+1)hxi+2|w,xi+1,xn(xi+2)…
hxn-1|w,xi+1,xi+2,…,xn-2,xn(xn-1)，

式中：相關變量組w包含i個變量，hw(w)有i!種表達形式；相關變量組w-包含n-i個變量，hw-|w(w-)有(n-i)!種表達形式；h(x)有i!(n-i)!種表達形式。

(17)

2.3 敏感系數轉化

(14)式代入(3)式，有(18)式成立。

(18)

式中：y為包含獨立隨機變量y1，y2，…，yn的變量集，y=[y1,y2,…,yn];v為包含獨立隨機變量y1，y2，…，yi的變量集，v=[y1,y2,…,yi];v-為包含獨立隨機變量y′i+1，yi+2，…，yn的變量集，v-=[yi+1,yi+2,…,yn]。

由(14)式取反，可得(19)式。

(19)

對(19)式所含各式兩邊同時微分，根據概率密度函數的性質[20]，有(20)式成立。

(20)

式中：gyi(yi)為連續獨立實隨機變量yi的概率密度函數。

(20)式中前i個等式相乘，有(21)式成立。

(21)

根據條件概率密度函數的性質[11]，有(22)式成立。

hx1,x2,…,xi(x1,x2,…,xi)=hxi(xi)hx1|xi(x1)·
hx2|x1,xi(x2)…hxi-1|x1,x2,…,xi-2,xi(xi-1)，

(22)

式中：hx1,x2,…,xi(x1,x2,…,xi)為相關變量x1,x2,…,xi的聯合概率概率函數。

w=[x1,x2,…,xi]代入(22)式，有(23)式成立。

hw(w)=hxi(xi)hx1|xi(x1)hx2|x1,xi(x2)…·
hxi-1|x1,x2,…,xi-2,xi(xi-1).

(23)

因為w=[x1,x2,…,xi]，有(24)式成立。

(24)

(23)式和(24)式代入(21)式，有(25)式成立。

(25)

(20)式中第i+1個到第n個等式相乘，有(26)式成立。

(26)

根據條件概率密度函數的性質[11]，有(27)式成立。

hxi+1,xi+2,…,xn|x1,x2,…,xi(xi+1,xi+2,…,xn)=
hxi+1|x1,x2,…,xi-1,xi(xi+1)…hxn|x1,x2,…,xi,…,xn-1(xn)，

(27)

式中：hxi+1,xi+2,…,xn|x1,x2,…,xi(xi+1,xi+2,…,xn)是給定變量x1,x2,…,xi條件下，變量xi+1,xi+2,…,xn的聯合概率密度函數。

w=[x1,x2,…,xi]和w-=[xi+1,xi+2,…,xn]代入(27)式，有(28)式成立。

hw-w(w-)=hxi+1|x1,x2,…,xi-1,xi(xi+1)…·
hxn|x1,x2,…,xi,…,xn-1(xn).

(28)

因為w-=[xi+1,xi+2,…,xn]，有(29)式成立。

(29)

(28)式、(29)式代入(26)式，可得(30)式。

(30)

2.3.1 相關變量組1階敏感系數轉化

(18)式、(25)式和(30)式代入(5)式，相關變量組w的1階敏感系數定義式可由(5)式轉換成(31)式。

(31)

參考(8)式，目標函數由F變成Fv時，獨立變量組v的1階敏感系數定義式如(32)式所示。

Sv|F=Fv=

(32)

式中：Sv|F=Fv為目標函數由F變成Fv時，獨立變量組v的1階敏感系數。

比較(31)式和(32)式，有(33)式成立。

Sw=Sv|F=Fv.

(33)

(33)式表明：相關變量組w的1階敏感系數可以轉化成，目標函數由F變成Fv條件下，獨立變量組v的1階敏感系數。

2.3.2 相關變量組全敏感系數轉化

(18)式、(25)式和(30)式代入(6)式，相關變量組w全敏感系數的定義式可由(6)式變成(34)式。

(34)

參考(9)式，目標函數由F變成Fv時，獨立變量組v全敏感系數定義式如(35)式所示。

(35)

比較(34)式和(35)式，有(36)式成立。

(36)

(36)式表明：相關變量組w的全敏感系數可以轉化成，目標函數由F變成Fv條件下，獨立變量組v的全敏感系數。

2.3.3 敏感系數性質分析

1) 當w包含任意相關變量時，可以先參照

(10)式～(13)式進行變量轉換(要求w包含的變量要先于其他變量被產生)，并將目標函數依照上述關系進行轉換；轉換后，w所包含相關變量對應獨立變量組的敏感系數就是w的敏感系數(在(10)式～(13)式中，xi和yi是一一對應的)。

2) 當w只包含一個相關變量時，w-包含其余相關變量時，就可以參照上述方法進行單個相關變量的1階敏感系數和全敏感系數的轉化。因此：單個相關變量的1階敏感系數是相關變量組1階敏感系數的一種特殊情況；單個相關變量的全敏感系數是相關變量組全敏感系數的一種特殊情況。

2.4 估計樣本的產生

為了實現對(5)式和(6)式所示敏感參數的估計，只需先產生兩組獨立變量yi的樣本。具體步驟如下：1)根據獨立變量yi的分布，利用蒙特卡洛法[21]，產生獨立變量yi的樣本；2)重復上述步驟，再次產生變量yi的樣本。詳細的取樣結果如2.5節中矩陣A和B所示。

2.5 敏感系數估計

(37)

式中：yri為變量yi第1次取樣的第r個樣本;m為樣本總數;y′ri為變量yi第2次取樣的第r個樣本。

(38)

3 面向彈丸炮口狀態的自行火炮結構參數全局靈敏度分析

為了降低面向彈丸炮口狀態的自行火炮結構參數全局靈敏度分析的計算量，首先進行面向射擊精度的彈丸炮口狀態全局靈敏度分析，找出對射擊精度影響較大的彈丸炮口狀態特征量；然后以上述影響較大的特征量表征彈丸炮口狀態，進行面向彈丸炮口狀態的自行火炮結構參數全局靈敏度分析。

3.1 面向射擊精度的彈丸炮口狀態全局靈敏度分析

面向射擊精度的彈丸炮口狀態全局靈敏度分析過程如下：首先基于表2所示對彈丸炮口狀態影響相對較大參數的分布情況，在第1節所示自行火炮整體發射動力學模型的幫助下，利用蒙特卡洛法統計出最大射程條件下彈丸炮口狀態的分布情況(結果如表3所示)；然后基于6自由度彈丸外彈道方程[15]和表3所示彈丸炮口狀態分布，在第2節所示全局靈敏度分析方法的幫助下，分別進行最大射程條件下面向彈丸落點縱向位移和橫向位移的彈丸炮口狀態全局靈敏度分析(結果見圖9和圖10)。

表2 參數分布

表3 彈丸炮口狀態分布

圖9 彈丸炮口狀態的敏感系數Fig.9 Sensitivity coefficients of projectile muzzle state

圖10 彈丸炮口狀態的敏感系數Fig.10 Sensitivity coefficients of position and attitude of projectile at nmuzzle

由圖9和圖10可知：彈丸炮口狀態所包含各特征量的1階敏感系數都比較小，同時它們的全敏感系數相對較大。造成這種現象的原因是，彈丸炮口狀態各個特征量主要通過它們之間的交叉項作用于最大射程條件下彈丸落點縱向和橫向位移。鑒于上述現象，可以通過圖9和圖10所示彈丸炮口狀態各特征量的全敏感系數，比較它們對最大射程條件下彈丸落點縱向(橫向)位移的敏感程度。根據上述分析，有如下結論：在表2所示條件下，彈丸質心炮口速率、彈丸炮口高低偏角、彈丸炮口方向偏角、彈丸炮口高低擺角速度和彈丸炮口方向擺角速度是影響某自行火炮最大射程條件下射擊精度的主要彈丸炮口狀態特征量。

3.2 面向彈丸炮口狀態的自行火炮結構參數全局靈敏度分析

根據3.1節的分析結果，為了減少計算量，本節將以對自行火炮射擊精度影響較大的彈丸炮口狀態特征量(彈丸質心炮口速率、彈丸炮口高低偏角、彈丸炮口方向偏角、彈丸炮口高低擺角速度和彈丸炮口方向擺角速度)作為彈丸炮口狀態的表征量，并以此為基礎進行面向彈丸炮口狀態的自行火炮結構參數全局靈敏度分析。

鑒于彈丸炮口狀態服從多元正態分布的情況(如表3所示)，根據(39)式計算彈丸炮口狀態的波動情況。

(39)

式中:Fj是對射擊精度影響較大的彈丸炮口狀態特征量；αj是Fj的權重系數；E(Fj)是Fj的均值；V(Fj)是Fj的方差。F越大，彈丸炮口狀態波動越大；F越小，彈丸炮口狀態波動越小。在本例中，計算所需的均值和方差如表3所示。因為最大射程條件下，彈丸落點縱向位移波動相對其橫向位移波動較大，相關彈丸炮口狀態特征量的權重系數直接取其最大射程條件下面向彈丸落點縱向位移的全敏感系數。

以彈丸炮口狀態波動為分析目標，表2所示參數為分析變量，某自行火炮最大射程條件下面向彈丸炮口狀態的全局靈敏度分析結果如圖11和圖12所示。在圖11和圖12中：自行火炮結構參數的1階敏感系數也都比較小，同時它們的全敏感系數也相對較大；相對內彈道參數，火炮整體結構參數和底盤參數的全敏感系數要大一些。造成火炮整體結構結構參數和底盤參數的全敏感系數較大的原因是：內彈道參數加工精度高，公差范圍??；火炮結構和底盤參數不易控制，波動范圍大。

圖11 結構參數敏感系數Fig.11 Sensitivity coefficients of structural parameters

圖12 結構參數敏感系數Fig.12 Sensitivity coefficients of structural parameters

為了進一步揭示自行火炮結構參數變化對最大射程條件下彈丸炮口狀態的敏感規律，按表2所示分類分別計算內彈道參數、火炮整體結構參數和底盤結構參數3個變量組相對彈丸炮口狀態的全敏感系數(結果如圖13所示)。

圖13 變量組敏感系數Fig.13 Sensitivity coefficients of variable sets

由圖13可知：內彈道參數、火炮整體結構參數和底盤結構參數3個變量組的全敏感系數均小于其包含所有變量的全敏感系數之和；內彈道參數、火炮整體結構參數和底盤結構參數3個變量組的全敏感系數之和大于1. 產生這種現象的原因是：自行火炮是一個復雜綜合性系統，不同變量組的全敏感系數間存在重疊的部分。這也進一步說明內彈道參數、火炮整體結構參數和底盤結構參數并不是單獨作用于彈丸炮口狀態，而是通過相互間的協同作用于彈丸炮口狀態。

鑒于上述現象，如果要了解自行火炮若干個變量同時對彈丸炮口狀態的敏感程度，最好將這些變量組合成一個變量組，然后計算這個變量組對彈丸炮口狀態的敏感系數。這樣能夠更加準確地評估該變量組對彈丸炮口狀態的敏感程度。借助變量組敏感系數的計算，可以更加準確地進行自行火炮不同部件間的敏感程度比較，為設計人員提供更加準確的靈敏度數據。

4 結果驗證

文獻[11]提供的全局靈敏度分析方法計算量較大，但是結果是正確的。為了實現對上述方法和相關分析結果的驗證，利用文獻[11]提供的全局靈敏度分析方法，進行同等條件下彈丸炮口狀態的全局靈敏度分析(詳細結果如表4和表5所示)。對比表4和表5所示結果和圖9～圖13所示結果后，發現：根據兩種不同的計算方法，彈丸炮口狀態各特征量的敏感系數是一致的；自行火炮各結構參數的敏感系數是一致的。這說明本文提供的全局靈敏度方法和相關靈敏度分析結果是可信的。同時估計表4和表5所示結果，文獻[15]提供的方法總共用了98 000組數據，本文提供的方法總共用了8 700組數據；估計表6所示結果，文獻[15]提供的方法總共用了70 000組數據，本文提供的方法總共用了6 500組數據。所以本文提供的方法提高了全局靈敏度分析方法的計算效率。

表4 面向彈丸落點縱向位移的敏感系數

表5 面向彈丸落點橫向位移的敏感系數

表6 面向彈丸炮口狀態的敏感系數

5 結論

為了研究給定分布條件下自行火炮結構參數對彈丸炮口狀態的敏感程度，本文建立了某自行火炮整體發射動力學模型，提出一種全局靈敏度分析方法，并應用上述模型和方法進行了實例分析。具體貢獻和結論如下：

1) 建立了能夠同時考慮彈丸膛內運動和車體振動的自行火炮整體發射動力學模型。

2) 提出了一種全局靈敏度分析方法。該方法證明了連續相關實向量可以用任意連續獨立實向量表示；輸入變量為任意連續獨立實變量條件下，通過對獨立變量組1階敏感系數和全敏感系數的估計，實現了相關變量組1階敏感系數和全敏感系數的估計。

3) 根據上述方法，進行了最大射程條件下，面向彈丸炮口狀態的某自行火炮結構數全局靈敏度分析。結果顯示：底盤縱向轉動慣量是影響最大射程條件下某自行火炮彈丸炮口狀態的最大因素。

參考文獻(References)

[1] 陳光宋. 彈炮耦合系統動力學及關鍵參數識別研究[D]. 南京：南京理工大學, 2016.

CHEN G S. The study on the dynamics of the projectile-barrel coupled system and the corresponding key parameters[D]. Nanjing: Nanjing University of Science and Technology, 2016. (in Chinese)

[2] 羅中峰,管小榮,徐亞棟,等.某火炮彈丸在炮口狀態的動態靈敏度分析[J].兵工學報,2017,38(12):2328-2336.

LUO Z F, GUAN X R, XU Y D, et al. Analyses of parameter sensitivity for position and attitude of projectile at muzzle[J]. Acta Armamentrii, 2017, 38(12):2328-2336. (in Chinese)

[3] 王曉鋒,姜興渭. 有限段法在高炮發射多體系統動力學分析中的應用[J]. 南京理工大學學報,2001,25(1):25-27.

WANG X F, JIANG X W. Application of finite segment method to dynamic analysis of antiaircraft gun launch[J]. Journal of Nanjing University of Science and Technology, 2001,25(1):25-27. (in Chinese)

[4] 劉飛飛,芮筱亭,于海龍, 等. 自行火炮行進間發射動力學研究[J]. 振動工程學報, 2016,29(3):380-385.

LIU F F, RUI X T, YU H L, et al. Study on launch dynamics of the self-propelled artillery marching fire [J]. Journal of Vibration Engineering, 2016,29(3):380-385. (in Chinese)

[5] 景鵬淵,顧克秋. 后坐體結構參量對彈丸起始擾動的影響[J]. 兵器裝備工程學報, 2017,38(2):53-56.

JING P Y，GU K Q．Impact of structure parameter of recoiling parts on the initial projectile disturbance[J]．Journal of Ordnance Equipment Engineering，2017，38(2):53-56. (in Chinese)

[6] SALTELLI A. Making best use of model evaluations to compute sensitivity indices[J]. Computer Physics Communications, 2002, 145(2): 280-297.

[7] SALTELLI A, TARANTOLA S. On the relative importance of input factor in mathematical models[J]. Statistical Association, 2002, 97(459): 8-30.

[8] KUCHERENKO S, TARANTOLA S, ANNONI P. Estimation of global sensitivity indices for models with dependent variables[J]. Computer Physics Communications, 2012, 183(4): 937-946.

[9] XU C G, GERTNER G Z. Uncertainty and sensitivity analysis for models with correlated parameters[J]. Reliability Engineering & System Safety, 2008, 93(10): 1563-1573.

[10] LI G Y, RABITZ H, YELVINGTON P E,et al. Global sensitivity analysis for systems with independent and/or correlated inputs[J]. Physical Chemistry A, 2010, 114(19): 6022-6032.

[11] MARA T A, TARANTOLA S, ANNONI P. Non-parametric methods for global sensitivity analysis of model output with dependent inputs [J]. Environmental Modelling & Software, 2015,72:173-183.

[12] TARANTOLA S, MARA T A. Variance-based sensitivity indices of computer models with dependent inputs: the Fourier amplitude sensitivity test[J]. International Journal for Uncertainty Quantification, 2017,7(6):511-523.

[13] DECARLO E C, MAHADEVAN S, SMARSLOK B P. Efficient global sensitivity analysis with correlated variables[J]. Structural and Multidisciplinary Optimization, 2018, 58(6): 2325-2340.

[14] ZHOU Y C, LU Z Z, XIAO S A, et al. Distance correlation-based method for global sensitivity analysis of models with dependent inputs[J]. Structural and Multidisciplinary Optimization, 2019,60(3):1189-1207.

[15] 韓子鵬. 彈箭外彈道學[M], 北京:北京理工大學出版社,2014.

HAN Z P. Exterior ballistics of projectiles and rockets[M]. Beijing: Beijing Institute of Technology Press,2014. (in Chinese)

[16] LANKARANI H M, NIKRAVESH P E. A contact force model with hysteresis damping for impact analysis of multibody system[J]. Journal of Mechanical Design, 1990, 112(3): 369 -376.

[17] WONG J Y, REECE A R. Prediction of rigid wheel performance based on the analysis of soil-wheel stresses part I. performance of driven rigid wheels[J]. Journal of Terramechanics, 1967,4(1):81-98.

[18] 羅中峰. 某特種機械的穩健設計方法研究[D]. 南京：南京理工大學,2020.

LUO Z F.Research on robustness design methods of a special machinery [D]. Nanjing: Nanjing University of Science and Technology, 2020. (in Chinese)

[19] ROSENBLATT M. Remarks on a multivariate transformation [J]. The Annals of Mathematical Statistics, 1952, 43: 470-472.

[20] 盛驟, 謝式千, 潘承毅. 概率論與數理統計[M]. 北京: 高等教育出版社, 2001.

SHENG Z, XIE S Q, PAN C Y.Probability theory and mathematical statistics[M]. Beijing: Higher Education Press, 2001. (in Chinese)

[21] METROPOLIS N, ULAM S. The Monte Carlo method[J]. Journal of the American Statistical Association, 1949, 44(247):335-341.