陳 華，周海兵，劉國昭，孫占峰，張樹道

（1．北京應(yīng)用物理與計(jì)算數(shù)學(xué)研究所，北京100094；2．中國工程物理研究院流體物理研究所沖擊波物理與爆轟物理實(shí)驗(yàn)室，四川綿陽621999）

圓筒試驗(yàn)JWL狀態(tài)方程參數(shù)的貝葉斯標(biāo)定＊

陳 華1，周海兵1，劉國昭1，孫占峰2，張樹道1

（1．北京應(yīng)用物理與計(jì)算數(shù)學(xué)研究所，北京100094；2．中國工程物理研究院流體物理研究所沖擊波物理與爆轟物理實(shí)驗(yàn)室，四川綿陽621999）

在研究數(shù)值模擬的輸入?yún)?shù)引入的不確定性時，通常需要人為給定每個輸入?yún)?shù)的概率分布，且輸入?yún)?shù)概率分布的選擇可能會對分析結(jié)果產(chǎn)生直接影響。利用貝葉斯方法標(biāo)定了圓筒試驗(yàn)JWL狀態(tài)方程參數(shù)，得到了標(biāo)定參數(shù)的估計(jì)值和后驗(yàn)分布，并研究了不同統(tǒng)計(jì)模型假設(shè)對標(biāo)定參數(shù)的估計(jì)值和后驗(yàn)分布的影響。貝葉斯后驗(yàn)分布融合了基準(zhǔn)試驗(yàn)的試驗(yàn)數(shù)據(jù)的信息，因此將其作為不確定度量化分析時輸入?yún)?shù)的初始概率分布，可以盡量減少分布選擇引入的認(rèn)知不確定性。

不確定度量化；參數(shù)標(biāo)定；后驗(yàn)概率分布；貝葉斯方法；JWL狀態(tài)方程；圓筒試驗(yàn)；高斯過程

隨著計(jì)算機(jī)技術(shù)的發(fā)展，數(shù)值模擬已經(jīng)成為復(fù)雜系統(tǒng)研究及設(shè)計(jì)的重要手段。以往，研究人員常常將材料參數(shù)、幾何尺寸等作為確定性的輸入?yún)?shù)，代入程序進(jìn)行數(shù)值模擬，從而得到確定性的結(jié)果。然而實(shí)際上，數(shù)值模擬過程中存在各種不確定性來源。近十幾年，不確定度量化（uncertainty quantification，UQ）方法發(fā)展迅速，逐漸成為數(shù)值模擬研究的熱點(diǎn)［12］。不確定度量化主要關(guān)注如何表征、量化、減小數(shù)值模擬過程中的各種不確定性來源，是實(shí)現(xiàn)高置信度數(shù)值模擬的有效手段。在進(jìn)行不確定度量化分析時，為了研究輸入?yún)?shù)引入的不確定性，每個輸入?yún)?shù)常常需要人為給定概率分布。輸入?yún)?shù)概率分布的選擇很可能會對分析結(jié)果有直接影響，為了減少分布選擇的主觀性，輸入?yún)?shù)的概率分布最好能夠由試驗(yàn)數(shù)據(jù)提供。

目前，圓筒試驗(yàn)是評估炸藥的做功能力、確定爆轟產(chǎn)物狀態(tài)方程參數(shù)的基準(zhǔn)試驗(yàn)，已經(jīng)得到了廣泛應(yīng)用［38］。標(biāo)準(zhǔn)的圓筒試驗(yàn)方法主要采用電離探針測量炸藥爆速，用高速掃描相機(jī)記錄定?；票Z驅(qū)動下圓筒壁的徑向膨脹距離隨時間的變化關(guān)系，并利用經(jīng)驗(yàn)公式計(jì)算圓筒壁的膨脹速度和比動能等表征炸藥做功能力的特征量，最后利用這些試驗(yàn)結(jié)果標(biāo)定JWL（Jones－Wilkins－Lee）狀態(tài)方程的參數(shù)［910］。本文中，利用貝葉斯方法對圓筒試驗(yàn)JWL狀態(tài)方程參數(shù)進(jìn)行標(biāo)定，給出標(biāo)定參數(shù)的估計(jì)值，得到標(biāo)定參數(shù)的后驗(yàn)分布，然后將該后驗(yàn)分布作為不確定度量化分析時輸入?yún)?shù)的初始分布，以減小輸入?yún)?shù)分布選擇引入的不確定性。

1 圓筒試驗(yàn)

圓筒試驗(yàn)的試驗(yàn)裝置示意圖如圖1所示，高速相機(jī)VISAR從狹縫掃描記錄的是圓筒壁外表面的徑向膨脹過程，通過對照片底片的處理獲得圓筒壁的徑向膨脹位移和時間的關(guān)系。得到的觀測數(shù)據(jù)如圖2所示。炸藥爆轟產(chǎn)物JWL狀態(tài)方程形式和等熵條件如下：式中：p為壓力，下標(biāo)“S”表示熵，V為相對比容，e為初始體積能量，A、B、C、R1、R2和ω為待標(biāo)定參數(shù)。利用JWL狀態(tài)方程在C－J點(diǎn)的特性，由圓筒試驗(yàn)和相關(guān)試驗(yàn)的觀測數(shù)據(jù)得到爆壓、爆速、密度、內(nèi)能和參數(shù)ω，再構(gòu)造代數(shù)方程組，給定R1和R2，聯(lián)立求解A、B和C［11］。

圖1 圓筒試驗(yàn)的試驗(yàn)裝置示意圖Fig．1 Device of cylinder test

圖2 圓筒試驗(yàn)徑向半徑和時間的關(guān)系Fig．2 Radius vs．time in cylinder test

2 貝葉斯標(biāo)定方法

M．Kennedy和A．OHagan［12］提出用貝葉斯方法進(jìn)行參數(shù)標(biāo)定和預(yù)測。由此發(fā)展出來一系列方法，統(tǒng)稱為KOH方法［1315］。KOH方法用高斯過程（Gaussian process，GP）對計(jì)算結(jié)果和觀測結(jié)果進(jìn)行建模，同時對模型偏差進(jìn)行建模。KOH方法的主要思想是用高斯過程建立代理模型，然后用代理模型計(jì)算后驗(yàn)分布和標(biāo)定參數(shù)。下面簡單介紹KOH方法的框架，更多技術(shù)細(xì)節(jié)參見文獻(xiàn)［12］。

2．1 定義和統(tǒng)計(jì)模型

令X＝（X（1），…，X（qX））為qX維控制變量，θ＝（θ（1），…，θ（qθ））為qθ維標(biāo)定參數(shù)，Y為數(shù)值模擬計(jì)算的輸出結(jié)果，Z為試驗(yàn)觀測結(jié)果。KOH模型如下：

式中：fC（·）為計(jì)算模型，b（·）為模型偏差，εj為觀測誤差，ρ為回歸系數(shù)，M 為數(shù)值模擬樣本量，N為試驗(yàn)樣本量，x＊i為控制變量在第i次數(shù)值模擬中的取值，xj為控制變量在第j次實(shí)驗(yàn)中的取值（按照數(shù)理統(tǒng)計(jì)中約定，采用大寫字母表示隨機(jī)變量，相應(yīng)的小寫字母表示該隨機(jī)變量的取值，本文中，“X”、“Y”、“Z”、“T”、“R1”、“R2”均采用該約定）。為了得到數(shù)值模擬的結(jié)果，需要給標(biāo)定參數(shù)賦值，令s表示數(shù)值模擬時輸入的標(biāo)定參數(shù)，θ表示標(biāo)定參數(shù)的真實(shí)值，以便于區(qū)分。

假設(shè)εj服從均值為零、方差為λ的正態(tài)分布；fC（·）用高斯過程建模，均值函數(shù)為m1（x，s），協(xié)方差函數(shù)為c1（·，·）；模型偏差b（x）也用高斯過程建模，均值函數(shù)為m2（x），協(xié)方差函數(shù)為c2（·，·）。假設(shè)均值函數(shù)和協(xié)方差函數(shù)的形式如下：

式中：hi（·）為統(tǒng)計(jì)建模時任意給定的連接函數(shù)，σ2i、ωXi、ωθi為協(xié)方差函數(shù)建模中需要估計(jì)的未知參數(shù)。令ψ＝（ψ1，ψ2），ψi是協(xié)方差函數(shù)ci（·，·）中的超參數(shù)（σ2i，ωXi，ωθi），i＝1，2。KOH模型的超參數(shù)包括β＝（β1，β2）和φ＝（ρ，λ，ψ），則全部要估計(jì)的參數(shù)為（θ，β，φ）。

2．2 數(shù)據(jù)結(jié)構(gòu)

數(shù)據(jù)包括觀測數(shù)據(jù)z和數(shù)值模擬結(jié)果y，即dT＝（zT，yT）。數(shù)值模擬過程中輸入變量的取值記為D1＝｛（x＊1，s1），…，（x＊M，sM）｝，試驗(yàn)數(shù)據(jù)控制變量的取值記為D2＝｛x1，…，xN｝，與D1相對應(yīng)，可以將D2擴(kuò)充為D2（θ）＝｛（x1，θ），…，（xN，θ）｝。

2．3 后驗(yàn)分布

令先驗(yàn)分布P（β）∝1。根據(jù)貝葉斯理論，在給定先驗(yàn)分布的條件下，后驗(yàn)分布為：

式中：f（d；md（θ），Vd（θ））是多元正態(tài)分布N（md（θ），Vd（θ））的密度函數(shù)，且

式中：H1、H2為由連接函數(shù)hi構(gòu)成的向量，V1（D1）為由集合D1中數(shù)據(jù)構(gòu)成的方差矩陣，V2（D2）為由集合D2中數(shù)據(jù)構(gòu)成的方差矩陣，C1為D1和D2中數(shù)據(jù)構(gòu)成的協(xié)方差矩陣，IN為N維單位矩陣。

2．4 計(jì)算步驟

直接計(jì)算式（6）中的后驗(yàn)分布較困難，因此將計(jì)算過程分為4個步驟：步驟1和步驟2是建立GP模型，步驟3是用建立的GP模型計(jì)算后驗(yàn)分布，步驟4是利用步驟3得到的后驗(yàn)分布標(biāo)定參數(shù)。

步驟1：建立GP代理模型，估計(jì)β1和c1（·，·）的超參數(shù)ψ1。用數(shù)值模擬所得數(shù)據(jù)y估計(jì)β1和ψ1，得到估計(jì)值步驟2：建立觀測數(shù)據(jù)的GP代理模型，估計(jì)ρ、λ、β2和c2（·，·）的超參數(shù)ψ2。在給定的條件下，用觀測數(shù)據(jù)z估計(jì)ρ、λ、β2和c2（·）的超參數(shù)ψ2，從而得到估計(jì)值步驟3：計(jì)算標(biāo)定參數(shù)θ的后驗(yàn)分布。由于，因此通過計(jì)算即可得到θ的后驗(yàn)分布。步驟4：利用θ的后驗(yàn)分布標(biāo)定參數(shù)。標(biāo)定參數(shù)的估計(jì)可以定義為后驗(yàn)分布的均值＾θ（mean）、中位數(shù)＾θ（mid）或最大概率值＾θ（max），其中＾θ（max）也稱為最大后驗(yàn)估計(jì)（maximum posterior estimates，MPE）。

3 JWL狀態(tài)方程參數(shù)的標(biāo)定

采用KOH貝葉斯方法標(biāo)定圓筒試驗(yàn)JWL狀態(tài)方程參數(shù)。將標(biāo)定參數(shù)記為（R1，R2）；將時間視為設(shè)計(jì)變量，記為T；同時令M＝30，N＝10。數(shù)值計(jì)算過程中，R1和R2的抽樣區(qū)間分別為［4．3，5．5］和［0．8，1．5］，用拉丁超立方抽樣（latin hypercube sample，LHS）方法得到M組參數(shù)并代入模擬程序進(jìn)行計(jì)算，得到M條徑向半徑隨時間變化的曲線。對設(shè)計(jì)變量T進(jìn)行LHS抽樣得到M個樣本。將T和（R1，R2）的樣本隨機(jī)配對得到模擬計(jì)算的輸入數(shù)據(jù)D1＝｛（t＊1，R11，R21），…，（t＊M，R1M，R2M）｝。D1對應(yīng)的計(jì)算輸出結(jié)果為y＝｛y1，…，yM｝。需要注意的是，每次數(shù)值計(jì)算的結(jié)果都是徑向半徑隨時間變化的曲線，這里只是從曲線上隨機(jī)取一個抽樣點(diǎn)作為輸出結(jié)果。對時間變量T隨機(jī)抽樣得到N 個樣本，即D2＝｛t1，…，tN｝。將D2擴(kuò)充為D2（r1，r2）＝｛（t1，r1，r2），…，（tN，r1，r2）｝，其中（r1，r2）是標(biāo)定參數(shù)的真實(shí)值。D2（r1，r2）對應(yīng)的徑向半徑觀測數(shù)據(jù)為z＝｛z1，…，zN｝。假設(shè)（R1，R2）的先驗(yàn)分布為均勻分布，即對（R1，R2）沒有任何知識；超參數(shù)ψ的先驗(yàn)分布為標(biāo)準(zhǔn)正態(tài)分布。數(shù)值模擬使用的計(jì)算程序?yàn)橛杀本?yīng)用物理與計(jì)算數(shù)學(xué)研究所自主研制的相容拉氏流體動力學(xué)工具包CHAP（compatible hydrodynamics analysis program）［16］。

下面分別采用兩個模型進(jìn)行建模，其主要區(qū)別在于模型中是否包含偏差。模型1是KOH模型的通用形式，如式（2）～（3）所示；模型2是KOH模型的特殊形式，相當(dāng)于在模型1的基礎(chǔ)上令ρ＝1，b（xj）≡0，此時式（3）可以簡化為：

按照第2．4節(jié)中的4個步驟依次進(jìn)行計(jì)算。第1步是建立數(shù)值模擬的GP代理模型。GP代理模型中參數(shù)的估計(jì)值分別為為檢驗(yàn)GP代理模型預(yù)測的能力，重新隨機(jī)抽樣4組（R1，R2）樣本代入數(shù)值模擬程序進(jìn)行計(jì)算，將計(jì)算得到的結(jié)果和采用GP模型預(yù)測的結(jié)果進(jìn)行比較，如圖3所示。從圖3可以看出：代理模型預(yù)測曲線和數(shù)值計(jì)算的曲線基本重合，說明代理模型的預(yù)測足夠精確。第2步是建立觀測數(shù)據(jù)的GP代理模型。模型1中模型參數(shù)的估計(jì)值分別為由于模型2忽略了模型偏差并且ρ＝1，所以該步驟中只需要估計(jì)λ，估計(jì)值為第3步是計(jì)算標(biāo)定參數(shù)的后驗(yàn)分布。圖4給出了標(biāo)定參數(shù)（R1，R2）的后驗(yàn)概率分布等高線，其中后驗(yàn)概率共分為9個概率區(qū)間，分別由9種顏色進(jìn)行區(qū)分，從小到大依次用Pk（k＝1，…，9）表示，P9表示最大概率區(qū)間的值域。圖4（a）為采用模型1得到的后驗(yàn)概率分布，可以看出，概率分布從右下角向左上角逐級遞增，概率最大區(qū)域P9位于左上角；圖4（b）為采用模型2得到的后驗(yàn)概率分布，可以看出，概率分布從右下角和左邊向中間區(qū)域逐級遞增，P9近似位于中間位置。第4步是標(biāo)定參數(shù)。由圖4可知：模型1的MPE為兩個模型的MPE均近似位于P9區(qū)域的中心位置。

圖3 GP代理模型預(yù)測的結(jié)果和數(shù)值計(jì)算結(jié)果的比較Fig．3 Comparison of predictions of the GP surrogate model with numerical simulations

圖4 兩模型的后驗(yàn)概率等高線和最大后驗(yàn)估計(jì)Fig．4 Posterior distribution and MPE with different models

圖5為模型1和模型2的預(yù)測效果，其中標(biāo)定結(jié)果（紅線）是將MPE代入數(shù)值模擬程序計(jì)算得到，預(yù)測結(jié)果（綠線）采用觀測數(shù)據(jù)的GP代理模型計(jì)算得到，試驗(yàn)結(jié)果用黑線表示。圖5顯示：兩個模型的預(yù)測結(jié)果均與試驗(yàn)結(jié)果符合很好，說明建立的代理模型都能夠比較精確地描述試驗(yàn)數(shù)據(jù)。模型2中試驗(yàn)結(jié)果和標(biāo)定結(jié)果符合得比較好，表明標(biāo)定的參數(shù)值是最佳擬合的估計(jì)；但是，當(dāng)t＞39μs時，模型1的標(biāo)定結(jié)果和試驗(yàn)結(jié)果差異逐漸增大。由于數(shù)值模擬和基于觀測數(shù)據(jù)代理模型的預(yù)測（即和）都足夠精確，所以兩者的差異即為模型偏差的預(yù)測＾b（x），說明模型1標(biāo)定的參數(shù)值是“調(diào)和”了模型偏差和標(biāo)定參數(shù)的最佳估計(jì)。

圖5 徑向半徑的預(yù)測、試驗(yàn)和模擬結(jié)果的比較Fig．5 Comparison of radius histories obtained by prediction，experiment and simulation

4 結(jié) 論

（1）用KOH貝葉斯方法標(biāo)定圓筒試驗(yàn)JWL狀態(tài)方程的參數(shù)R1和R2，得到了標(biāo)定參數(shù)的估計(jì)值和后驗(yàn)分布，并將后驗(yàn)分布作為子系統(tǒng)級（或系統(tǒng)級）復(fù)雜系統(tǒng)數(shù)值模擬不確定度量化分析的輸入?yún)?shù)的初始分布，以盡可能地減少分布選擇引入的認(rèn)知不確定性。（2）分析了兩種KOH統(tǒng)計(jì)模型假設(shè)對標(biāo)定參數(shù)估計(jì)值和后驗(yàn)分布的影響。從參數(shù)標(biāo)定的角度分析，模型2得到的估計(jì)值更好；從不確定度量化的角度分析，由于模型1還分析了模型偏差引入的不確定性，因此模型1更合適。

［1］ Oberkampf W L，Roy C J．Verification and validation in scientific computing［M］．Cambridge：Cambridge University Press，2010．

［2］ Smith R C．Uncertainty quantification：Theory，implementation，and applications［M］．Philadelphia：Society for Industrial and Applied Mathematics，2013．

［3］ Kury J W，Honig H C，Lee E L，et al．Metal acceleration by chemical explosives［C］∥4th Symposium on Detonation．Maryland：White Oak，1966：3－13．

［4］ Lan I F，Hung S C，Chen C Y，et al．An improved simple method of deducing JWL parameters from cylinder expansion test［J］．Propellants，Explosives，Pyrotechnics，1993，18（1）：18－24．

［5］ 于川，劉文翰，李良忠，等．鈍感炸藥圓筒試驗(yàn)與爆轟產(chǎn)物JWL狀態(tài)方程研究［J］．高壓物理學(xué)報，1997，11（3）：227－233．Yu Chuan，Liu Wenhan，Li Liangzhong，et al．Studies on cylinder test and JWL equation of state of detonation product for insensitive high explosive［J］．Chinese Journal of High Pressure Physics，1997，11（3）：227－233．．

［6］ 陳朗，黃毅民，馮長根．含鋁炸藥圓筒試驗(yàn)及爆轟產(chǎn)物JWL狀態(tài)方程研究［J］．火炸藥學(xué)報，2001，24（3）：13－15．Chen Lang，Huang Yimin，F(xiàn)eng Changgen．The cylinder test and JWL equation of state detontion product of aluminized explosives［J］．Chinese Journal of Explosives and Propellants，2001，24（3）：13－15．

［7］ 陳朗，李澤仁．激光速度干涉儀測量法在炸藥圓筒試驗(yàn)中的應(yīng)用［J］．爆炸與沖擊，2001，21（3）：229－232．Chen Lang，Li Zeren．The cylinder tests measured by VISAR interferometer［J］．Explosion and Shock Waves，2001，21（3）：229－232．

［8］ 徐輝，孫占峰，李慶忠．標(biāo)準(zhǔn)圓筒試驗(yàn)數(shù)據(jù)處理和不確定度評定方法［J］．北京理工大學(xué)學(xué)報，2010，30（5）：626－630．Xu Hui，Sun Zhanfeng，Li Qingzhong．Study on data processing and uncertainty evaluation of standard cylinder test［J］．Transactions of Beijing Institute of Technology，2010，30（5）：626－630．

［9］ Lee E L，Tarver C M．Phenomenological model of shock initiation in heterogeneous explosives［J］．Physics of Fluids，1980，23（12）：2362－2372．

［10］ Wescott B L，Stewart D S，Davis W C．Equation of state and reaction rate for condensed－phase explosives［J］．Journal of Applied Physics，2005，98（5）：053514．

［11］ 孫承緯，衛(wèi)玉章，周之奎．應(yīng)用爆轟物理［M］．北京：國防工業(yè)出版社，2000．

［12］ Kennedy M，OHagan A．Bayesian calibration of computer models（with discussion）［J］．Journal of the Royal Statistical Society Series B：Statistical Methodology，2001，68：425－464．

［13］ Williams B，Higdon D，Gattiker J，et al．Combining experimental data and computer simulations，with an application to flyer plate experiments［J］．Bayesian Analysis，2006，1（1）：765－792．

［14］ Higdon D，Gattiker J，Williams B，et al．Computer model calibration using high－dimensional output［J］．Journal of the American Statistical Association，2008，103（482）：570－583．

［15］ Bayarri M J，Berger J O，Rui P，et al．A framework for validation of computer models［J］．Technometrics，2007，49（2）：138－154．

［16］ 張樹道，周海兵，熊俊．相容拉氏流體動力學(xué)CHAP程序研制及其應(yīng)用研究：GF－A0091846G［R］．2005．

Bayesian calibration for parameters of JWL equation of state in cylinder test

Chen Hua1，Zhou Haibing1，Liu Guozhao1，Sun Zhanfeng2，Zhang Shudao1

（1．Institute of Applied Physics and Computational Mathematics，Beijing100094，China；2．Laboratory for Shock Wave and Detonation Physics Research，Institute of Fluid Physics，Chinese Academy of Engineering Physics，Mianyang621999，Sichuan，China）

Since probability distributions of input parameters are usually assigned subjectively for uncertainty quantification（UQ）in numerical simulations，the selection of distributions may have significant effect on analysis results of UQ．In this paper，to calibrate more precisely the parameters of JWL equation of state in the cylinder test，we proposed to adopt Bayesian methods that provide estimators and posterior distributions of calibration parameters．Further，we investigated the effects of model assumptions on estimators and posterior distributions of calibration parameters．Our study shows that，owing to the information of cylinder tests they contain，the posterior distributions can be used as initial probability distributions of the input parameters in UQ to reduce the epistemic uncertainty．

uncertainty quantification；calibration；posterior distribution；Bayesian methods；JWL equation of state；cylinder test；Gaussian process

O381國標(biāo)學(xué)科代碼：13035

A

10．11883／1001－1455（2017）04－0585－06

（責(zé)任編輯 王玉鋒）

2015－12－30；

2016－04－25

國家自然科學(xué)基金項(xiàng)目（11472060，11671413）；中國工程物理研究院科學(xué)技術(shù)發(fā)展基金項(xiàng)目（2013A0101004，2015B0201037）；北京應(yīng)用物理與計(jì)算數(shù)學(xué)研究所所長青年基金項(xiàng)目（ZYSZ1518－13）

陳 華（1979－ ），女，博士，副研究員，chen＿h(yuǎn)ua＠iapcm．a(chǎn)c．cn。