BP-GA 算法確定未反應(yīng)炸藥的JWL 狀態(tài)方程參數(shù)

崔 浩，郭 銳，宋 浦，顧曉輝，周 昊，楊永亮，江 琳，俞旸暉

（1. 南京理工大學(xué)機(jī)械工程學(xué)院，江蘇 南京 210094；2. 西安近代化學(xué)研究所 燃燒與爆炸技術(shù)重點(diǎn)實(shí)驗(yàn)室，陜西 西安 710065；3. 南京理工大學(xué) 國(guó)家特種超細(xì)粉體工程技術(shù)研究中心，江蘇 南京 210094）

1 引言

炸藥穩(wěn)定爆轟時(shí)，前導(dǎo)沖擊波壓縮炸藥從初始狀態(tài)沿著C-J 爆轟Rayleigh 線突躍到未反應(yīng)Hugoniot 線的馮諾依曼峰值（Von Neuman Spike，VNS），隨后化學(xué)反應(yīng)使炸藥的狀態(tài)沿C-J 爆轟Rayleigh 線回到爆轟Hugoniot 曲線上的C-J 狀態(tài)，C-J 爆轟Rayleigh 線和爆轟Hugoniot 曲線相切是爆轟波穩(wěn)定傳播的條件［1］。當(dāng)沖擊波強(qiáng)度較低時(shí)，從初始點(diǎn)出發(fā)的Rayleigh 線的斜率低于C-J 爆轟Rayleigh 線的斜率，導(dǎo)致炸藥無(wú)法穩(wěn)定爆轟，內(nèi)部存在未反應(yīng)的部分。未反應(yīng)炸藥的狀態(tài)方程即描述前導(dǎo)沖擊波陣面過(guò)后未反應(yīng)炸藥的壓力、比容和內(nèi)能等物理量之間關(guān)系的函數(shù)，在數(shù)值仿真軟件中未反應(yīng)炸藥的狀態(tài)方程通常采用Jones-Wilkings-Lee（JWL）狀態(tài)方程［2-3］形式。未反應(yīng)炸藥JWL 狀態(tài)方程參數(shù)的確定是進(jìn)行反應(yīng)度流場(chǎng)計(jì)算、點(diǎn)火增長(zhǎng)模型擬合和炸藥裝藥沖擊起爆仿真等工作的基礎(chǔ)。

基于JWL 或Grüneisen 狀態(tài)方程以及波陣面上的守恒關(guān)系可推導(dǎo)得到未反應(yīng)炸藥JWL 形式的Hugoniot 曲線表達(dá)式，據(jù)此可擬合出JWL 狀態(tài)方程中的未知參數(shù)，因此炸藥沖擊Hugoniot 曲線的測(cè)量是確定未反應(yīng)炸藥狀態(tài)方程的前提［4-5］。關(guān)于炸藥沖擊Hugoniot 關(guān)系的測(cè)量國(guó)內(nèi)外學(xué)者進(jìn)行了大量的研究，主要通過(guò)平面準(zhǔn)一維沖擊試驗(yàn)測(cè)量炸藥內(nèi)部的沖擊波速度和粒子速度關(guān)系而獲得，包括采用激光干涉測(cè)速儀或電磁速度計(jì)直接測(cè)量粒子速度［6-9］和采用錳銅壓力計(jì)測(cè)量炸藥內(nèi)部壓力并結(jié)合壓力對(duì)比法間接計(jì)算得到粒子速度［10-11］等方式。通過(guò)試驗(yàn)數(shù)據(jù)獲得粒子速度后再根據(jù)動(dòng)量守恒關(guān)系或沖擊波速度擬合法推導(dǎo)得到?jīng)_擊波速度，基于多組粒子速度-沖擊波速度試驗(yàn)數(shù)據(jù)即可擬合出未反應(yīng)炸藥的Hugoniot 關(guān)系。此外，眾多學(xué)者也開(kāi)展了基于Grüneisen 狀態(tài)方程和Hugoniot 關(guān)系標(biāo)定未反應(yīng)炸藥JWL 狀態(tài)方程的研究。其中，王延飛等［12］根據(jù)試驗(yàn)獲得的沖擊Hugoniot關(guān)系對(duì)三項(xiàng)式狀態(tài)方程進(jìn)行了標(biāo)定，確定了JOB-9003 炸藥的JWL 狀態(tài)方程參數(shù)；黃奎邦等［13］基于Grüneisen 狀態(tài)方程推導(dǎo)得到了JWL 形式的Hugoniot 曲線表達(dá)式，并提出了用最小二乘法擬合JWL 參數(shù)的方法。上述研究給出了未反應(yīng)炸藥JWL 參數(shù)的推導(dǎo)過(guò)程和參數(shù)擬合流程，但目前尚未出現(xiàn)利用智能算法計(jì)算未反應(yīng)炸藥JWL 參數(shù)的研究。遺傳算法（Genetic Algorithm，GA）是目前最常用的尋優(yōu)算法之一，單一的GA 算法通常用于求解非線性系統(tǒng)的最優(yōu)解，但有些非線性系統(tǒng)無(wú)法建立精確的數(shù)學(xué)模型，而B(niǎo)P 神經(jīng)網(wǎng)絡(luò)-遺傳算法（BP-GA 算法）不僅可以尋求到這些非線性系統(tǒng)的最優(yōu)解［14］，還可直接預(yù)測(cè)染色體的適應(yīng)度值，提高求解效率［15］。

為此，本研究提出了一種利用BP-GA 算法計(jì)算未反應(yīng)炸藥JWL 狀態(tài)方程參數(shù)的方法。此方法基于JWL狀態(tài)方程和波陣面上的守恒方程推導(dǎo)得到的JWL 形式的Hugoniot 曲線表達(dá)式，利用3 個(gè)約束方程并選取未反應(yīng)炸藥沖擊Hugoniot 曲線上的5 個(gè)特定點(diǎn)對(duì)試驗(yàn)數(shù)據(jù)確定的Hugoniot 曲線進(jìn)行擬合，最終可確定未反應(yīng)炸藥的JWL 參數(shù)。文中采用這種方法確定了8 種炸藥的JWL 參數(shù)，結(jié)果表明通過(guò)該方法獲得的曲線和試驗(yàn)數(shù)據(jù)相吻合，能夠滿足工程應(yīng)用的需要。

2 沖擊Hugoniot 關(guān)系和JWL 狀態(tài)方程

2.1 沖擊Hugoniot 關(guān)系

在沖擊試驗(yàn)中，炸藥初始態(tài)的粒子速度和壓力通常為0，因此波陣面上的質(zhì)量、動(dòng)量和能量守恒方程分別為［16］：

式中，ρ0和ρH分別為炸藥初始和波后密度，g·cm-3；US為波陣面速度，km·s-1；UP為波后粒子速度，km·s-1；pH為波后炸藥內(nèi)部壓力，GPa；e0和eH分別為炸藥初始和波后比內(nèi)能，kJ·g-1；VH0和VH為別為炸藥初始和波后比容，cm3·g-1。

通常情況下，材料內(nèi)的沖擊波速度US和粒子速度UP呈線性關(guān)系，即：

式中，v=VH/VH0，為相對(duì)比容。由于Hugoniot 系數(shù)C0和S根據(jù)試驗(yàn)數(shù)據(jù)確定，因此式（9）即為試驗(yàn)數(shù)據(jù)確定的未反應(yīng)炸藥p-v形式的沖擊Hugoniot 曲線。

2.2 JWL 狀態(tài)方程

JWL 狀態(tài)方程定義的壓力為［2-3］：

將初始相對(duì)比容v=1 和初始溫度T0=298.15 K 代入JWL 形式的內(nèi)能表達(dá)式可得未反應(yīng)炸藥的初始內(nèi)能E0，將E0代入式（12）可得完整的JWL 形式的Hugoniot 曲線表達(dá)式為［13，17］：

式中，CV為比定容熱容，GPa·K-1。

為了減少未反應(yīng)炸藥JWL 狀態(tài)方程中未知參數(shù)的個(gè)數(shù)以及獲得一組封閉的JWL 狀態(tài)方程參數(shù)［18］。根據(jù)未反應(yīng)炸藥的沖擊Hugoniot 關(guān)系，JWL 參數(shù)之間需滿足如下3 個(gè)約束方程：

（1）初始態(tài)時(shí)，v=1 時(shí)，pHJWL=0，可得［13］：

（2）當(dāng)v=vVNS時(shí)，pHJWL=pVNS，即：

（3）此外，為了模擬炸藥延遲爆轟時(shí)內(nèi)部真實(shí)張力水平，通常設(shè)置R1和R2有如下簡(jiǎn)單的線性關(guān)系［19］：

上述3 個(gè)約束方程使得JWL 狀態(tài)方程中的6 個(gè)參數(shù)只有3 個(gè)是獨(dú)立的，只需確定其中任意3 個(gè)參數(shù)便可確定剩余參數(shù)。對(duì)于一組JWL 參數(shù)，當(dāng)式（9）與式（13）的誤差值在［vVNS，1］區(qū)間上均等于或接近0 時(shí)，可認(rèn)為此組JWL參數(shù)符合炸藥真實(shí)的JWL參數(shù)，為了提高擬合精度，選取多個(gè)特定點(diǎn)進(jìn)行計(jì)算?？紤]到不同炸藥的［vVNS，1］區(qū)間相差較多，不方便選取具體相對(duì)比容值，因而提出以vVNS為參考點(diǎn)，均勻選取vVNS+0.05、vVNS+0.10、vVNS+0.15、vVNS+0.20 和vVNS+0.25 這5 個(gè)點(diǎn)為特定點(diǎn)的方法。為了定量評(píng)價(jià)一組JWL 參數(shù)的適用性，基于選取的5 個(gè)特定點(diǎn)定義兩條曲線的相似度函數(shù)為：

觀察式（17）可發(fā)現(xiàn)，函數(shù)值F取值區(qū)間為［0，1］，當(dāng)函數(shù)值越接近1 時(shí)，可認(rèn)為兩條曲線越相似，吻合程度越高。

3 BP-GA 算法和計(jì)算結(jié)果

BP（Back Propagation）神經(jīng)網(wǎng)絡(luò)是一種按誤差反向傳播訓(xùn)練的多層前饋網(wǎng)絡(luò)，其基本思想是梯度下降法，利用梯度搜索技術(shù)，以期使網(wǎng)絡(luò)的實(shí)際輸出值和期望輸出值的誤差均方誤差為最小。BP 算法主要包括信息的正向傳播和誤差的反向傳播兩個(gè)過(guò)程。網(wǎng)絡(luò)正向傳播時(shí)，輸入信號(hào)通過(guò)隱含層作用于輸出節(jié)點(diǎn)，若實(shí)際輸出與期望輸出相差較大，則轉(zhuǎn)入反向傳播，并將輸出誤差通過(guò)隱含層向輸入層逐層反傳，通過(guò)調(diào)整各層的神經(jīng)元閾值，使誤差沿梯度方向下降，經(jīng)過(guò)反復(fù)學(xué)習(xí)訓(xùn)練可確定與最小誤差相對(duì)應(yīng)的網(wǎng)絡(luò)參數(shù)。

為了尋求非線性系統(tǒng)（不同的JWL 參數(shù)組合代表的F值組成的非線性系統(tǒng)）中的最優(yōu)解，通常需要借助GA 優(yōu)秀的搜尋能力，然而由于JWL 參數(shù)的取值范圍較大，這要求GA 的種群規(guī)模要足夠大，且需要求解3 個(gè)約束方程才可計(jì)算染色體的適應(yīng)度值（本研究染色體適應(yīng)度值為式（17）的F）。種群規(guī)模的擴(kuò)大、種群進(jìn)化次數(shù)的延長(zhǎng)使得染色體適應(yīng)度值的計(jì)算次數(shù)大大增多，直接降低了求解過(guò)程的計(jì)算效率。為了提高計(jì)算效率，本研究對(duì)BP 神經(jīng)網(wǎng)絡(luò)進(jìn)行訓(xùn)練，使其可以擬合非線性系統(tǒng)，從而直接預(yù)測(cè)每條染色體的適應(yīng)度值。

3.1 BP 神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)及訓(xùn)練過(guò)程

由前述可知，只需確定JWL 參數(shù)中的任意3 個(gè)便可通過(guò)求解3 個(gè)約束方程得到剩余參數(shù)，從而可確定F值?？紤]到JWL 參數(shù)中的R1、R2、ω和CV取值范圍相對(duì)較小，且R1和R2有固定的線性關(guān)系，選取R1、ω和CV為自變量。國(guó)內(nèi)外眾多研究表明大多數(shù)炸藥的R1、ω和CV的取值范圍分別集中在6.2～14.1、0.8～1.251 和1.68×10-3～2.78×10-3GPa·K-1［20-24］。為了尋求JWL參數(shù)組合的多樣性和最優(yōu)解，本研究將R1、ω和CV的取值范圍分別放大為5～16、0.6～2.6和1×10-3～3×10-3GPa·K-1。

文獻(xiàn)［25］研究發(fā)現(xiàn)雙隱含層BP 神經(jīng)網(wǎng)絡(luò)的訓(xùn)練誤差低于單隱含層，但隨著隱含層數(shù)的增加，網(wǎng)絡(luò)會(huì)更加復(fù)雜，網(wǎng)絡(luò)權(quán)值和閾值的訓(xùn)練時(shí)間也會(huì)增加。為了達(dá)到較高的訓(xùn)練精度以及縮短訓(xùn)練時(shí)間，本研究采用雙隱含層的網(wǎng)絡(luò)結(jié)構(gòu)，如圖1 所示，包括輸入層、雙隱含層和輸出層，輸入層節(jié)點(diǎn)數(shù)為3（自變量R1、ω、CV）、輸出層節(jié)點(diǎn)數(shù)為1（染色體適應(yīng)度F）。由于自變量的取值范圍較大，在自變量取值范圍內(nèi)隨機(jī)生成10000個(gè)樣本進(jìn)行訓(xùn)練，通過(guò)最佳隱含層計(jì)算公式［26］確定第一隱含層節(jié)點(diǎn)數(shù)為40，并采用試湊法［27］確定第二隱含層的節(jié)點(diǎn)數(shù)為30。

圖1 BP 神經(jīng)網(wǎng)絡(luò)示意圖Fig.1 Schematic diagram of BP neural network

為了獲得高精度的神經(jīng)網(wǎng)絡(luò)，設(shè)置當(dāng)網(wǎng)絡(luò)訓(xùn)練集的均方誤差（Mean Square Error，MSE）低于1×10-10時(shí)停止訓(xùn)練，訓(xùn)練步驟具體如下［15］：

（1）網(wǎng)絡(luò)初始化，網(wǎng)絡(luò)的權(quán)值和閾值隨機(jī)賦初值。

（2）信息正向傳播，計(jì)算各層的輸入和輸出。設(shè)神經(jīng)網(wǎng)絡(luò)有P個(gè)訓(xùn)練樣本，對(duì)任意訓(xùn)練樣本Xk［xk1，xk2，xk3］，期望輸出為T［tk］。計(jì)算各層神經(jīng)元節(jié)點(diǎn)輸入、輸出分別為：

式中，η為學(xué)習(xí)步長(zhǎng)。

（5）記憶訓(xùn)練

重復(fù)步驟（2）～（4），正向傳播和反向傳播交替進(jìn)行，不斷逐層修正權(quán)值和閾值直至訓(xùn)練集的均方誤差（Mean Square Error，MSE）低于1×10-10。

本研究訓(xùn)練集、驗(yàn)證集和測(cè)試集的比例（依次為0.7∶0.15∶0.15）采用Matlab 軟件默認(rèn)比例，其各自的均方誤差隨訓(xùn)練次數(shù)的變化曲線如圖2 所示，圖中縱坐標(biāo)為了觀察方便采用對(duì)數(shù)形式。觀察圖2 可發(fā)現(xiàn)，訓(xùn)練集的MSE 在335 代時(shí)低于預(yù)期設(shè)置的1×10-10，網(wǎng)絡(luò)訓(xùn)練結(jié)束，而驗(yàn)證集和測(cè)試集的MSE 也接近1×10-10，表明訓(xùn)練后的BP 神經(jīng)網(wǎng)絡(luò)誤差較小、預(yù)測(cè)精度較高。此外，訓(xùn)練結(jié)束后的訓(xùn)練集、驗(yàn)證集和測(cè)試集的決定系數(shù)R2均為1。

圖2 均方誤差隨訓(xùn)練次數(shù)的變化曲線Fig.2 The change of MSE versus training times

3.2 GA 流程及參數(shù)設(shè)置

GA 的染色體編碼方式、選擇機(jī)制、交叉方式和變異算子等均同文獻(xiàn)［28］，具體步驟如下：

（1）選擇，從種群中選擇一定數(shù)量的染色體；

（2）交叉，選擇出的染色體兩兩交叉產(chǎn)生下一代；

（3）變異，交叉生成的染色體以一定概率發(fā)生變異；

（4）精英保存，用上一代種群中適應(yīng)性最好的染色體替換掉本代適應(yīng)性最差的個(gè)體。

（5）種群替換，經(jīng)過(guò)步驟（1）～（4）后種群進(jìn)化了一代，重復(fù)此步驟直到種群進(jìn)化成熟。

由于自變量取值范圍較大，設(shè)置種群中染色體的數(shù)量為200 個(gè)，同時(shí)為了尋求到全局最優(yōu)解，設(shè)置種群連續(xù)50 代無(wú)進(jìn)化為種群進(jìn)化成熟標(biāo)準(zhǔn)。在種群進(jìn)化過(guò)程中，選擇機(jī)制、交叉方式、變異算子和精英保存等基因操作均需要計(jì)算染色體的適應(yīng)度值，這些適應(yīng)度值采用訓(xùn)練好的BP 神經(jīng)網(wǎng)絡(luò)進(jìn)行直接預(yù)測(cè)。

3.3 BP-GA 算法流程

圖3 為BP-GA 算法確定未反應(yīng)炸藥JWL 狀態(tài)方程參數(shù)的計(jì)算流程，BP-GA 算法主要包括BP 神經(jīng)的訓(xùn)練過(guò)程以及隨后的GA 尋優(yōu)過(guò)程，具體流程為：（1）BP神經(jīng)網(wǎng)絡(luò)的訓(xùn)練過(guò)程，輸入炸藥的初始密度、爆速、Hugoniot 系數(shù)C0和S后隨機(jī)生成10000 個(gè)樣本進(jìn)行神經(jīng)網(wǎng)絡(luò)的訓(xùn)練，當(dāng)訓(xùn)練樣本的MSE 小于設(shè)定的1×10-10時(shí)停止訓(xùn)練，此時(shí)訓(xùn)練好的BP 神經(jīng)網(wǎng)絡(luò)即可預(yù)測(cè)任意一組JWL 參數(shù)對(duì)應(yīng)的F值；（2）GA 的尋優(yōu)過(guò)程，初始化種群，設(shè)置種群數(shù)量200 個(gè)，經(jīng)過(guò)選擇、交叉、變異和精英保存等基因操作后種群進(jìn)化一代，循環(huán)多次直到種群進(jìn)化成熟，進(jìn)化成熟后種群中適應(yīng)度值最大的染色體攜帶的遺傳信息即為JWL 參數(shù)。

圖3 BP-GA 算法計(jì)算未反應(yīng)炸藥JWL 參數(shù)流程圖Fig.3 Flow chart of calculating JWL parameters of unreacted explosive by BP-GA algorithm

3.4 計(jì)算結(jié)果和試驗(yàn)對(duì)比

表1 炸藥材料參數(shù)Table 1 Parameters of explosives

表2 炸藥未反應(yīng)JWL 狀態(tài)方程參數(shù)Table 2 Parameters of JWL equation of state for unreacted explosives

參數(shù)所確定p-v曲線和試驗(yàn)數(shù)據(jù)確定的p-v曲線的對(duì)比如圖4 所示，并采用決定系數(shù)R2定量分析兩者的相似性。觀察圖4 可發(fā)現(xiàn)，除TATB 炸藥的R2為0.9995外，其余7 種炸藥的R2均為1.0，表明吻合程度較好，證明了本研究提出的采用BP-GA 算法確定未反應(yīng)炸藥JWL 狀態(tài)方程參數(shù)的精度高。

圖4 試驗(yàn)和BP-GA 算法確定未反應(yīng)炸藥沖擊的p-v 形式的Hugoniot 曲線對(duì)比Fig.4 Comparison of Hugoniot curves in p-v form of unreacted explosive determined by experiment and BP-GA Algorithm

4 結(jié)論

提出了一種基于BP-GA 算法和沖擊Hugoniot 關(guān)系確定未反應(yīng)炸藥的JWL 狀態(tài)方程參數(shù)的方法，采用這種方法確定了8 種炸藥的JWL 狀態(tài)方程參數(shù)，主要結(jié)論如下：

（1）基于JWL 狀態(tài)方程和波陣面上的守恒方程可推導(dǎo)得到JWL 形式的Hugoniot 曲線表達(dá)式，利用3 個(gè)約束方程和Hugoniot 曲線上的5 個(gè)特定點(diǎn)并采用BP-GA 算法對(duì)試驗(yàn)數(shù)據(jù)確定的Hugoniot 曲線進(jìn)行擬合，可確定未反應(yīng)炸藥的JWL 參數(shù)。

（2）已知某種炸藥的初始密度、爆速、Hugoniot系數(shù)C0和S，便可利用BP-GA 算法確定其JWL 參數(shù)，結(jié)果顯示BP-GA 算法確定的8 種未反應(yīng)炸藥的p-v曲線和試驗(yàn)數(shù)據(jù)確定的p-v曲線相一致，且8 條p-v曲線的R2均不低于0.9995。