面向彈丸炮口狀態(tài)的自行火炮結(jié)構(gòu)參數(shù)全局靈敏度

羅中峰， 管小榮， 徐誠

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

0 引言

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

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

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

1 自行火炮整體發(fā)射動力學(xué)模型

1.1 彈丸炮口狀態(tài)定義

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

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

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

1.2 自行火炮整體發(fā)射動力學(xué)模型

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

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

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

1.2.1.1 彈丸膛內(nèi)狀態(tài)定義

為了獲得彈丸炮口位置和姿態(tài)，參照圖1構(gòu)建彈丸膛內(nèi)位置和姿態(tài)的表征方法，如圖2所示。

圖2 彈丸膛內(nèi)位置和姿態(tài)Fig.2 Position and attitude of projectile in bore

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

當(dāng)彈丸到達炮口時，圖2所示的彈丸膛內(nèi)位置和姿態(tài)就會和圖1所示的彈丸炮口位置和姿態(tài)對應(yīng)相等。原因是：當(dāng)彈丸到達炮口時，由于身管和彈丸間相互作用的結(jié)束，彈丸由柔性體變成剛體，彈丸膛內(nèi)質(zhì)心和膛內(nèi)彈軸就分別會與彈丸炮口質(zhì)心和炮口彈軸重合?；谏鲜鲈颍谟嬎氵^程中通過持續(xù)記錄膛內(nèi)質(zhì)心和膛內(nèi)彈軸位置變化情況，就可以計算出彈丸炮口的位置和姿態(tài)。

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

為了考慮彈丸的膛內(nèi)運動，根據(jù)有限段思想建立了如圖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 自行火炮發(fā)射動力學(xué)模型

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

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

1.2.2.1 模型拓撲結(jié)構(gòu)

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

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

1.2.2.2 碰撞力計算

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

(1)

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

為了準(zhǔn)確計算地面與輪胎間的接觸力，上述自行火炮發(fā)射動力學(xué)模型根據(jù)(2)式所示W(wǎng)ong-Reece模型[17]進行相關(guān)接觸力的計算。

(2)

1.3 模型驗證

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

表1 參數(shù)列表

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

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

2 基于Rosenblatt轉(zhuǎn)換理論的全局靈敏度分析方法

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

2.1 敏感系數(shù)定義

本節(jié)所考慮的目標(biāo)函數(shù)如(3)式所示。

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

(3)

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

假設(shè)(3)式所示目標(biāo)函數(shù)擁有有限方差。當(dāng)(3)式所示目標(biāo)函數(shù)擁有有限方差時，其方差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(·)為函數(shù)期望；F(x)|w為在給定變量組w條件下F(x)的取值；w|w-為在給定變量組w-的條件下變量組w的取值。

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

(5)

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

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

(6)

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

當(dāng)x1,x2,…,xn都是獨立變量時，(4)式可以轉(zhuǎn)換成(7)式[6]，(5)式可以轉(zhuǎn)換成(8)式[6]，(6)式可以轉(zhuǎn)換成(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的邊緣概率密度函數(shù)，k=1，2，…，n.

(9)

在后續(xù)段落中，視x1,x2,…,xn為相關(guān)變量。

2.2 相關(guān)變量的產(chǎn)生

2.2.1 相關(guān)變量的產(chǎn)生

根據(jù)Rosenblatt轉(zhuǎn)換[19]，有(10)式成立。

(10)

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

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

(11)

同樣根據(jù)Rosenblatt轉(zhuǎn)換規(guī)則[19]，有(12)式成立。

(12)

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

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

(13)

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

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

(14)

2.2.2 相關(guān)變量產(chǎn)生性質(zhì)分析

當(dāng)參照(14)式產(chǎn)生相關(guān)變量時，根據(jù)概率密度函數(shù)的性質(zhì)[11]，相關(guān)變量(組)的分布情況如(15)式所示。

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

對比(15)式中各個變量(組)的概率密度函數(shù)和(16)式中向量x的聯(lián)合概率密度函數(shù)形式，可知： (14)式是參照(16)式所示聯(lián)合概率密度函數(shù)的一種形式(hw(w)=hxi(xi)hx1|xi(x1)hx2|x1,xi(x2)…h(huán)xi-1|x1,x2,…,xi-2,xi(xi-1);hw-|w(w-)=hxi+1|w(xi+1)hxi+2|w,xi+1(xi+2)…h(huán)xn|w,xi+1,xi+2,…,xn-1(xn))產(chǎn)生的相關(guān)變量。在2.2.1節(jié)所示方法的幫助下，通過調(diào)整相關(guān)變量的產(chǎn)生順序，(14)式也可以按照(16)式所示聯(lián)合概率密度函數(shù)的其他形式產(chǎn)生相關(guān)變量。比如，參照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)…h(huán)xn|w,xi+1,xi+2,…,xn-1(xn)，有(17)式成立。

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

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

(15)

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

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)…h(huán)xi|x1,x2,…,xi-1(xi-1)…=
hx2(x2)hx1|x2(x1)hx3|x1,x2(x3)…h(huán)xi|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)，

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

(17)

2.3 敏感系數(shù)轉(zhuǎn)化

(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)式所含各式兩邊同時微分，根據(jù)概率密度函數(shù)的性質(zhì)[20]，有(20)式成立。

(20)

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

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

(21)

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

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

(22)

式中：hx1,x2,…,xi(x1,x2,…,xi)為相關(guān)變量x1,x2,…,xi的聯(lián)合概率概率函數(shù)。

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)

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

hxi+1,xi+2,…,xn|x1,x2,…,xi(xi+1,xi+2,…,xn)=
hxi+1|x1,x2,…,xi-1,xi(xi+1)…h(huán)xn|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的聯(lián)合概率密度函數(shù)。

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 相關(guān)變量組1階敏感系數(shù)轉(zhuǎn)化

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

(31)

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

Sv|F=Fv=

(32)

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

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

Sw=Sv|F=Fv.

(33)

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

2.3.2 相關(guān)變量組全敏感系數(shù)轉(zhuǎn)化

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

(34)

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

(35)

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

(36)

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

2.3.3 敏感系數(shù)性質(zhì)分析

1) 當(dāng)w包含任意相關(guān)變量時，可以先參照

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

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

2.4 估計樣本的產(chǎn)生

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

2.5 敏感系數(shù)估計

(37)

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

(38)

3 面向彈丸炮口狀態(tài)的自行火炮結(jié)構(gòu)參數(shù)全局靈敏度分析

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

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

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

表2 參數(shù)分布

表3 彈丸炮口狀態(tài)分布

圖9 彈丸炮口狀態(tài)的敏感系數(shù)Fig.9 Sensitivity coefficients of projectile muzzle state

圖10 彈丸炮口狀態(tài)的敏感系數(shù)Fig.10 Sensitivity coefficients of position and attitude of projectile at nmuzzle

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

3.2 面向彈丸炮口狀態(tài)的自行火炮結(jié)構(gòu)參數(shù)全局靈敏度分析

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

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

(39)

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

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

圖11 結(jié)構(gòu)參數(shù)敏感系數(shù)Fig.11 Sensitivity coefficients of structural parameters

圖12 結(jié)構(gòu)參數(shù)敏感系數(shù)Fig.12 Sensitivity coefficients of structural parameters

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

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

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

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

4 結(jié)果驗證

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

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

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

表6 面向彈丸炮口狀態(tài)的敏感系數(shù)

5 結(jié)論

為了研究給定分布條件下自行火炮結(jié)構(gòu)參數(shù)對彈丸炮口狀態(tài)的敏感程度，本文建立了某自行火炮整體發(fā)射動力學(xué)模型，提出一種全局靈敏度分析方法，并應(yīng)用上述模型和方法進行了實例分析。具體貢獻和結(jié)論如下：

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

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

3) 根據(jù)上述方法，進行了最大射程條件下，面向彈丸炮口狀態(tài)的某自行火炮結(jié)構(gòu)數(shù)全局靈敏度分析。結(jié)果顯示：底盤縱向轉(zhuǎn)動慣量是影響最大射程條件下某自行火炮彈丸炮口狀態(tài)的最大因素。

參考文獻(References)

[1] 陳光宋. 彈炮耦合系統(tǒng)動力學(xué)及關(guān)鍵參數(shù)識別研究[D]. 南京：南京理工大學(xué), 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] 羅中峰,管小榮,徐亞棟,等.某火炮彈丸在炮口狀態(tài)的動態(tài)靈敏度分析[J].兵工學(xué)報,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] 王曉鋒,姜興渭. 有限段法在高炮發(fā)射多體系統(tǒng)動力學(xué)分析中的應(yīng)用[J]. 南京理工大學(xué)學(xué)報,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] 劉飛飛,芮筱亭,于海龍, 等. 自行火炮行進間發(fā)射動力學(xué)研究[J]. 振動工程學(xué)報, 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] 景鵬淵,顧克秋. 后坐體結(jié)構(gòu)參量對彈丸起始擾動的影響[J]. 兵器裝備工程學(xué)報, 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] 韓子鵬. 彈箭外彈道學(xué)[M], 北京:北京理工大學(xué)出版社,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] 羅中峰. 某特種機械的穩(wěn)健設(shè)計方法研究[D]. 南京：南京理工大學(xué),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] 盛驟, 謝式千, 潘承毅. 概率論與數(shù)理統(tǒng)計[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.