多目標(biāo)優(yōu)化算法在彈頭侵徹明膠運動模擬中的應(yīng)用*

莫根林，金永喜，李忠新，吳志林

(1.江蘇大學(xué)先進制造研究院，江蘇 鎮(zhèn)江 212013；2.中國兵器工業(yè)第208研究所瞬態(tài)沖擊技術(shù)重點實驗室，北京 102202；3.南京理工大學(xué)機械工程學(xué)院，江蘇 南京 210094)

為揭示彈頭侵徹人體組織的致傷機理，通常采用彈道明膠作為人體組織的替代物[1-2]。程可[3]、S.S.Liu等[4]和Y.Wen等[5]、溫垚珂等[6]建立了步槍彈侵徹明膠的數(shù)值模擬模型；D.Datoc[7]、Y.Wang等[8]、G.H.Yoon等[9]建立了手槍彈侵徹明膠的數(shù)值模擬模型。劉坤等[10]、Liu等[11]基于彈頭所受合力和合力矩的經(jīng)驗假設(shè)建立了彈頭侵徹明膠的平面運動模型；莫根林等[12]基于彈頭表面微元的受力假設(shè)建立了彈頭的六自由彈道模型。研究表明這些模型均能較好地模擬對應(yīng)投射物的運動規(guī)律，但模型中的力學(xué)系數(shù)和彈頭的一些初始運動參數(shù)一般不是由實驗直接得到的[1]，而是由試算的方法進行推測的。由于待定參數(shù)較多，為獲得較好的彈頭運動模擬結(jié)果，試算過程會消耗大量的時間。

待定參數(shù)的試算過程可看作是一個多目標(biāo)問題的求解過程。傳統(tǒng)求解此類問題的方法一般是利用線性加權(quán)法、約束法等將其轉(zhuǎn)化為單目標(biāo)問題，但其存在探索未知空間的能力不強、容易陷入局部極值點等問題，不能很好地處理復(fù)雜、非線性問題[13]。而多目標(biāo)進化算法如多目標(biāo)遺傳算法(MOGA)、Pareto存檔進化策略(PAES)、非劣分類遺傳算法(NSGA)以及基于遺傳算法的多目標(biāo)優(yōu)化算法(Gamultiobj)等，具有全局性和多向性的優(yōu)點，可以很好地解決此類問題[14]。為此，本文中以待定參數(shù)作為優(yōu)化變量建立彈頭侵徹明膠運動模擬的多目標(biāo)優(yōu)化模型，并運用Gamultiobj算法和決策策略對優(yōu)化變量進行求解分析。

1 彈道模型

建立全局坐標(biāo)系O′-x′y′z′，其中x′軸平行于明膠塊底面，y′軸豎直向上，z′軸指向明膠塊內(nèi)部，如圖1所示。在彈頭上建立連體坐標(biāo)系O-xyz，其中z軸和彈軸重合，O點位于彈頭的質(zhì)心位置。采用歐拉角描述彈頭的空間姿態(tài)，即全局坐標(biāo)系依次繞z軸旋轉(zhuǎn)ψ(進動角)得到輔助坐標(biāo)系O′-x″y″z″；再繞x″軸旋轉(zhuǎn)θ(章動角)得到輔助坐標(biāo)系O′-ξηζ，最后繞ζ軸旋轉(zhuǎn)φ(自轉(zhuǎn)角)即可得到坐標(biāo)系O-xyz,如圖2所示。

根據(jù)文獻[10]，彈頭侵徹明膠過程中所受到的阻力FD、升力FL以及翻滾力矩M的大小可表示為：

(1)

式中：ρ為明膠密度、v為彈頭速度、A0為彈頭的橫截面積；δ為彈頭速度方向和彈軸方向的夾角(攻角)；CD0和C1為阻力系數(shù)，CL0為升力系數(shù)，CM0為彈頭的轉(zhuǎn)矩系數(shù)。阻力FD方向與彈頭速度方向相反，升力FL方向垂直于速度方向且與FD和彈軸共面，翻滾力矩M垂直于FD和FL。

令彈頭速度方向余弦為(λ1，λ2，λ3)，則彈頭的水平速度分量vz、豎直速度分量vy和側(cè)向速度分量vx可以表示為：

vz=vλ3,vy=vλ2,vx=vλ1

(2)

令彈軸在全局坐標(biāo)系中的方向余弦為(λ4，λ5，λ6)，可由歐拉角表示為[15]：

λ4=sinθsinψ,λ5=sinθcosψ,λ6=cosθ

(3)

令升力FL的方向余弦為(λ7，λ8，λ9)，翻滾力矩M的方向余弦為(λ10，λ11，λ12)，由定義可知，它們可由彈頭速度的方向余弦和彈軸的方向余弦表示：

(4)

λ10=λ2λ6-λ3λ5,λ11=λ3λ4-λ1λ6,λ12=λ1λ5-λ2λ4

(5)

在連體坐標(biāo)系中，令M的方向余弦為(λ13，λ14，λ15)，它與(λ10，λ11，λ12)的關(guān)系可表示為[15]：

(6)

因此，M在連體坐標(biāo)系中的分量Mx、My和Mz可以表示為：

Mx=Mλ13,My=Mλ14,Mz=Mλ15

(7)

根據(jù)推廣的歐拉動力學(xué)方程，彈頭的彈道方程可表示為:

(8)

(9)

式中：Jx、Jy和Jz分別為彈頭繞x、y、z軸的轉(zhuǎn)動慣量；ωx、ωy和ωz為彈頭繞x、y、z軸的角速度分量，它們與歐拉角角速率的關(guān)系為：

(10)

將式(1)～(7)和式(9)～(10)帶入式(8)，并通過龍格-庫塔法求解，可得彈頭運動位移和空間歐拉角隨時間的變化規(guī)律。

2 多目標(biāo)優(yōu)化模型與算法

2.1 目標(biāo)函數(shù)組

令彈軸在x′z′平面上的投影與z′軸的夾角(偏航角)為α′、彈軸在y′z′平面上投影與z′軸的夾角(俯仰角)為β′，它們可由彈頭的歐拉角計算得到：

tanα′=cosψtanθ, tanβ′=sinψtanθ

(11)

令彈頭質(zhì)心坐標(biāo)的實驗值為(xO，yO，zO)，彈頭偏航角的實驗值為α，俯仰角的實驗值為β。以彈頭質(zhì)心坐標(biāo)、偏航角和俯仰角的理論值均方根誤差表示彈道模型模擬彈頭運動的好壞程度，則目標(biāo)函數(shù)組可表示為：

(12)

式中:N為實驗數(shù)據(jù)的總量。

2.2 優(yōu)化變量

2.3 優(yōu)化算法

Gamultiobj算法的基本流程如圖3所示。初始種群一般由隨機產(chǎn)生的多個個體組成，每個個體由待求解的優(yōu)化變量組成；該算法首先判斷初始種群是否滿足用戶設(shè)定條件，若滿足則得到Pareto解集(每個解至少存在一個目標(biāo)值優(yōu)于集合中所有其他的解)；否則將使種群進化一代，并再一次對其判斷，循環(huán)往復(fù)直至滿足終止條件。

種群的進化流程如圖4所示。首先，基于個體的序值和擁擠距離，采用錦標(biāo)賽算法從上一種群(父種群)中選擇一定數(shù)量的個體，其中序值小的個體先被選中；序值相同時，擁擠距離大的個體被選中。其次對選擇得到的個體進行交叉、變異操作，以產(chǎn)生一個由新個體的集合即子種群。其次，將子種群和父種群合并成一個新的種群(此時新種群的大小是父種群的兩倍)。最后，通過修剪算法從合并種群中選擇與父種群大小相同的新種群，其中修剪算法是基于序值和擁擠距離進行的，序值高、擁擠距離小的個體被剔除[16]。

2.4 決策策略

為從Pareto解集中挑選出最優(yōu)解，引入TOPSIS綜合評價方法對解集中的解進行綜合評價[17]。其決策過程可概括為：

(1) 將多目標(biāo)函數(shù)的解構(gòu)成決策矩陣及規(guī)范化決策矩陣；

(2) 確定近似理想解，一般由解集中同一目標(biāo)的最小值構(gòu)成；

(3) 求各個解與理想解的近似度。一般首先求出各個解與理想解的距離，其次對多個目標(biāo)賦予權(quán)重，最后將各個權(quán)重與各個距離的乘積之和作為近似度；

(4) 根據(jù)近似度的大小，對解集中的解進行排序；

(5) 選出最優(yōu)的多目標(biāo)函數(shù)解[18]。

3 算 例

3.1 實驗

實驗所用明膠塊的制作方法為：將明膠顆粒和水按照1∶9的質(zhì)量比例混合，并放入水浴爐中加熱至60 ℃。保溫2 h后將澄清的凝膠液體注入300 mm×300 mm×300 mm的模具中。待模具中的明膠液體冷卻至室溫時，將其放入4C的恒溫箱中保溫24 h，經(jīng)測量其密度為1 060 kg/m3。實驗時，以81式7.62 mm自動步槍瞄準(zhǔn)明膠塊中心位置射擊56式7.62 mm普通彈，發(fā)射時槍口距離明膠塊約2.2 m。實驗示意圖如圖5所示。鏡面放置在明膠塊上方，明膠塊側(cè)面的高速攝像機同時拍攝彈頭在明膠塊中的運動和彈頭在鏡面中像的運動。高速攝像機的采樣頻率為30 000 s-1、曝光時間為2s、圖像分辨率為640×480。實驗所得彈頭侵徹明膠的運動過程如圖6所示。通過Phantom Camera Control軟件的測量功能從拍攝的圖片上測量彈頭運動的實驗值xO、yO、zO、α和β。

3.2 約束條件

對于彈道模型中的力學(xué)系數(shù)，在文獻[10]的基礎(chǔ)上通過試算確定它們的上、下邊界,如表1所示。對于彈頭的初始運動參數(shù)，通過對水平位移zO實驗值的差分運算可得vz0=680 m/s，考慮差分誤差的影響，限定vz0取值范圍為650～730 m/s；由于實驗中通過瞄準(zhǔn)使彈頭的入射方向基本垂直于明膠塊的入射面，因此假定彈頭速度方向和水平方向的夾角小于3.5。再由vx、vy和vz的關(guān)系可知，vx和vy的取值范圍均為-30～30 m/s；由圖6可知彈頭的初始章動角θ0較小，考慮誤差的影響，設(shè)定其取值范圍為-5～5；考慮到彈頭軸線方向可能為空間的任意取值，設(shè)定彈頭進動角ψ0的范圍為0～180；由文獻[19]可知，彈頭出槍口時的自轉(zhuǎn)角速率φt0=1 357 rad/s、進動角速率ψt0=17 997 rad/s，以它們2～3倍的數(shù)值估計它們的取值范圍,即0<φt0<3 000 rad/s、0<ψt0<40 000 rad/s。由于一般認(rèn)為彈頭的章動角速率θt0較小，設(shè)定其取值范圍為0～2 000 rad/s。

表1 待優(yōu)化力學(xué)系數(shù)的取值范圍Table 1 Bounds of undetermined mechanical coefficients

3.3 計算結(jié)果

設(shè)定Gamultiobj算法的具體參數(shù)為：種群的大小，200；初始種群的生成，隨機法；變異操作，自適應(yīng)可行策略；交叉操作，中間策略；優(yōu)化的終止條件，迭代次數(shù)達(dá)600次或平均適應(yīng)值函數(shù)的相對誤差小于0.5×10-4。由于多目標(biāo)函數(shù)組(12)中的各個目標(biāo)之間相互獨立，令TOPSIS算法中各個目標(biāo)的權(quán)重相同均為0.2。

序號C1CD0CL0CM0vx0/(m·s-1)vy0/(m·s-1)vz0/(m·s-1)θ0/(°)ψ0/(°)θt0/(rad·s-1)ψt0/(rad·s-1)φt0/(rad·s-1)17.490.058 60.6360.064 7-22.4-14.207131.91068.81 0201 22015 40028.890.066 20.6690.082 915.3-17.407030.11568.81 4701 91017 60037.720.051 30.7440.091 123.54.107200.28181.41 2202 20019 100410.900.064 20.7230.071 918.37.737170.72881.91 4101 84012 80054.830.047 80.8130.062 815.6-16.907150.68270.51 7502 26016 900

4 討 論

3.3節(jié)中Gamultiobj算法通過167次迭代獲得了優(yōu)化變量的最優(yōu)解，由于種群的大小為200，相當(dāng)于進行了33 400次彈道模型的計算。若采用搜索算法對目標(biāo)函數(shù)組(12)的12個優(yōu)化變量進行優(yōu)化，即使每個變量取10個數(shù)值，模型的運算次數(shù)也高達(dá)1012次，遠(yuǎn)大于Gamultiobj的計算次數(shù)，可見Gamultiobj算法能夠大幅度縮減目標(biāo)函數(shù)組優(yōu)化所需的時間。

由Gamultiobj算法的基本原理可知，Pareto解集中的解與初始種群的構(gòu)成有關(guān)[16]。表2中第2～5組最優(yōu)解為采用隨機生成法生成初始種群時其他4次優(yōu)化得到的；它們對應(yīng)的目標(biāo)函數(shù)值如表3所示?？梢婋m然每次優(yōu)化得到的最優(yōu)解各不相同，但各最優(yōu)解對應(yīng)的均方根誤差均較小，表明它們均能較好地模擬彈頭的運動規(guī)律。由計算可知，最優(yōu)解中力學(xué)系數(shù)C1、CD0、CL0、CM0的變異系數(shù)分別為0.278、0.138、0.095 8和0.162，彈頭初始參數(shù)vx0、vy0、vz0、θ0、ψ0、θt0、ψt0和φt0對應(yīng)的變異系數(shù)分別為1.83、1.67、0.009、0.945、0.091 1、0.200、0.224和0.146。除了由于vx0、vy0、θ0的均值接近零，導(dǎo)致其變異系數(shù)較大外，其他優(yōu)化變量的變異系數(shù)均較小，表明所得最優(yōu)解具有較強可信性。

表3 最優(yōu)解對應(yīng)的目標(biāo)函數(shù)值Table 3 Values of objective functions corresponding to optimal solutions

5 結(jié) 論

在彈頭侵徹明膠彈道模型的基礎(chǔ)上，以彈頭空間位置和姿態(tài)理論均方根誤差作為目標(biāo)函數(shù)，對7.62 mm步槍彈侵徹明膠的運動進行了模擬分析，得到以下結(jié)論：

(1) Gamultiobj算法和TOPSIS分析方法獲得的彈頭初始運動參數(shù)和彈道模型力學(xué)系數(shù)，能較好地模擬步槍彈侵徹明膠的空間運動過程，有助于揭示彈頭對人體組織的致傷機理。

(2) Gamultiobj算法和TOPSIS分析方法能夠快速獲得優(yōu)化變量的多組最優(yōu)解。

(3) 優(yōu)化變量的最優(yōu)解與初始種群的構(gòu)成有關(guān)，但所得的最優(yōu)解變異系數(shù)較小，結(jié)果具有較強的可信性。