多目標優化算法在彈頭侵徹明膠運動模擬中的應用*

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

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

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

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

1 彈道模型

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

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

(1)

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

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

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

(2)

令彈軸在全局坐標系中的方向余弦為(λ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)

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

(6)

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

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

(7)

根據推廣的歐拉動力學方程，彈頭的彈道方程可表示為:

(8)

(9)

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

(10)

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

2 多目標優化模型與算法

2.1 目標函數組

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

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

(11)

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

(12)

式中:N為實驗數據的總量。

2.2 優化變量

2.3 優化算法

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

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

2.4 決策策略

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

(1) 將多目標函數的解構成決策矩陣及規范化決策矩陣；

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

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

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

(5) 選出最優的多目標函數解[18]。

3 算 例

3.1 實驗

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

3.2 約束條件

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

表1 待優化力學系數的取值范圍Table 1 Bounds of undetermined mechanical coefficients

3.3 計算結果

設定Gamultiobj算法的具體參數為：種群的大小，200；初始種群的生成，隨機法；變異操作，自適應可行策略；交叉操作，中間策略；優化的終止條件，迭代次數達600次或平均適應值函數的相對誤差小于0.5×10-4。由于多目標函數組(12)中的各個目標之間相互獨立，令TOPSIS算法中各個目標的權重相同均為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節中Gamultiobj算法通過167次迭代獲得了優化變量的最優解，由于種群的大小為200，相當于進行了33 400次彈道模型的計算。若采用搜索算法對目標函數組(12)的12個優化變量進行優化，即使每個變量取10個數值，模型的運算次數也高達1012次，遠大于Gamultiobj的計算次數，可見Gamultiobj算法能夠大幅度縮減目標函數組優化所需的時間。

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

表3 最優解對應的目標函數值Table 3 Values of objective functions corresponding to optimal solutions

5 結 論

在彈頭侵徹明膠彈道模型的基礎上，以彈頭空間位置和姿態理論均方根誤差作為目標函數，對7.62 mm步槍彈侵徹明膠的運動進行了模擬分析，得到以下結論：

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

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

(3) 優化變量的最優解與初始種群的構成有關，但所得的最優解變異系數較小，結果具有較強的可信性。