恒溫法測爆熱的快速系統辨識方法*

楊 杰，賀元吉，趙宏偉，徐明利，高洪泉

（中國人民解放軍96901部隊，北京 100094）

爆熱是評價炸藥作功能力的重要指標[1]，常用測量方法包括恒溫法、絕熱法與水下爆炸法[2-3]。絕熱法和恒溫法為直接且常用測量方法，兩種方法均需獲取測量結束時的內桶水溫度（終點溫度），區別在于，絕熱法外桶溫度需跟蹤內桶溫度、修正溫升（爆炸產物向內桶系統傳熱引起內桶系統的溫升）等于內桶溫升，而恒溫法外桶溫度恒定、修正溫升等于補償后的內桶溫升[4]；水下爆炸法為間接測量方法，需將沖擊波能與氣泡能折算至炸藥爆熱，誤差較大[5]。

恒溫法由于可測藥量大（可達100 g）、測量精度高、硬件簡單可靠，而成為測量炸藥爆熱的首選[1]。GJB 772A—97《炸藥試驗方法》中方法701.1中，恒溫法的初期（外桶水升溫結束至炸藥起爆時間段）、末期（爆炸產物基本完成放熱之后時間段）時間均為11 min，研究表明，如初期、主期（炸藥起爆至爆炸產物基本完成放熱時間段）、末期分別延長至21、80和21 min 時，測量誤差將大幅減小[4]。

然而，恒溫法測量時間比較長，如在主期或末期突遇異常溫控、突然斷電、軟件崩潰等系統故障，會導致無法獲取內桶水準確終點溫度，此時常規測量方法宣告失敗，造成很大的人力與物力的損失。因此，爆熱的估算也是爆熱測量的重要研究內容。俞統昌等[6]提出了混合炸藥爆熱的經驗估算公式；韓早[7]提出了半經驗公式法估算含鋁炸藥的爆熱；此外，還有量子力學方法、神經網絡方法、支持向量機方法等爆熱估算方法[8-10]。可是，這些估算均為事前估算，目前尚未開展測量過程中爆熱實時估算的研究。因此，本文中擬開展基于故障前內桶水溫數據辨識炸藥爆熱的研究，提出爆熱辨識算法，通過實驗測試算法收斂性能，為爆熱的實時估計提供理論依據。

1 量熱計結構及恒溫法測爆熱經典算法

量熱計（見圖1）測爆熱原理為：炸藥起爆后，爆炸產物向內桶系統（由爆熱彈金屬殼、內桶水及其外殼組成）傳熱，根據內桶水的溫升計算爆炸產物的放熱量，再除以炸藥質量即得到炸藥的爆熱。

根據GJB 772A—97中方法701.1，恒溫法測量炸藥爆熱的經典計算公式為：

式中：TI(τ0)、TI(τn)、n分別為主期的初溫（℃）、末溫（℃）、溫度讀數個數；線性擬合初期、末期內桶水溫升曲線得到的起點溫度（℃）分別為θ0、θn，溫升速率（℃/min）分別為V0、Vn；Δτ 為采樣間隔(min)；Δθ、ΔtJ、Δt為補償溫度、經典方法修正溫升、所有方法修正溫升(℃)；C為量熱計的系統熱容（即內桶系統熱容）(J/℃)；M1、M2為被測炸藥、傳爆藥的質量，g；q為雷管爆熱（J）；Q1、Q2分別為被測炸藥、傳爆藥的爆熱（J/g）。

2 爆熱測量過程的傳熱分析

2.1 外桶系統向內桶系統的傳熱分析

由圖1，內外桶之間的傳熱主要有熱傳導與熱輻射兩種方式，內外桶溫差較小時（若要求輻射傳熱量線性計算相對誤差不超過2%，外桶溫度為25℃時，內外桶溫差須不超過4℃），這兩種方式的傳熱量均與內外桶溫差為線性關系[11]：

式中：ΦI2(τ)、QI2(τ)分別為外桶系統在時刻τ 向內桶系統傳熱的熱流量、熱量（為便于計算，QI2(τ)在初期、主末期階段的計算起點分別為初期起始時刻τS、起爆時刻τM），TI(τ)、TW(τ)分別為內桶水、外桶水在時刻τ 的溫度（對于恒溫法，外桶水溫度在初始調溫、初期、主期、末期等階段均為常量），kNW為綜合內外桶熱傳導與熱輻射的總導熱系數。

在初始調溫階段，內桶水溫須高于外桶水溫（內外桶水溫差記為TC）才能將攪拌產生的熱量傳導出去，進而維持內外桶溫度恒定，因此根據能量守恒定律與式(4)得到攪拌產熱速率為：

式中：QC(τ)為攪拌熱，且產熱速率恒定（QC(τ)、QI2(τ)在初始調溫階段的計算起點均為該階段起始時刻）。

在恒溫法中，外桶水溫度TW進入初期階段后是一個常量。而在絕熱法中，外桶水溫度TW(τ)需跟蹤內桶水溫度的變化，且與內桶水溫度的關系為：

外桶水溫度表達式為TW時，表示該溫度為常量；表達式為TW(τ)時，表示該溫度為變量。

2.2 爆炸產物向內桶系統的傳熱分析

爆炸產物向內桶系統的傳熱分為兩個階段，第一階段是炸藥爆轟與爆炸產物飛散過程對內桶系統輻射傳熱，第二階段是爆炸產物充滿爆熱彈后對內桶系統傳熱。

宋浦等[12]通過含鋁炸藥光輻射實驗研究指出，炸藥爆轟及無氧燃燒反應階段的熱輻射基本為可見光輻射，而可見光輻射總能量僅占總爆熱的0.02%，因此，第一階段傳熱對總爆熱的影響可忽略不計。

第二階段中：爆炸產物充滿爆熱彈后的極短時間內，它在爆熱彈內劇烈振蕩（壓強振蕩時間不超過20 ms[2]），此階段爆炸產物與爆熱彈殼體的傳熱方式為強對流傳熱；隨著爆熱彈內壓強振蕩趨穩，爆炸產物與爆熱彈殼體的傳熱方式轉為弱對流傳熱；當爆炸產物完全穩定后，爆炸產物與爆熱彈殼體之間的傳熱方式變為熱傳導。因此，弱對流傳熱與熱傳導是爆炸產物向內桶系統傳熱的主要方式。

由于爆熱彈殼體（鋼材質）的導熱系數高、殼體較薄，弱對流與熱傳導狀態的爆炸產物導熱系數相對較小，內桶水對流強度不高，因此，爆熱彈殼體的溫度場沿殼體徑向近似為線性[11]。

由于爆熱彈殼體導熱良好且較薄，爆炸產物初溫非常高[13]，因此，局部升溫階段時間很短。而外桶系統傳熱與攪拌產熱速率很小，因此，該階段外桶系統傳熱量與攪拌產熱量均近似為零，內桶水及其外殼溫度在該階段近似恒定，進而得到該階段爆熱彈殼體溫度場分布：

圖2 爆炸產物的傳熱過程Fig.2 Heat transfer process of explosive product

式中：t(x,τ)為在時刻τ 距離爆熱彈殼體內壁（徑向）x處的溫度，TE(τ)為爆炸產物在時刻τ 的溫度，L為爆熱彈殼體的平均厚度。

由式(7)，局部升溫階段爆熱彈殼體吸收的熱量為：

式中：QN1(τ)為時刻τM至τ 爆炸產物的傳熱量，cG、mG分別為爆熱彈殼體的比熱容、質量。

炸藥爆轟后的氧化、燃燒等二次反應時間極短（無氧環境不超過0.3s[2]），因此，二次反應階段爆炸產物因傳熱導致熱量損失可忽略，可視為瞬態升溫過程，隨后爆炸產物一直處于放熱過程，因此有：

式中：cZ、mZ分別為爆炸產物的比熱容、質量。

在全局升溫階段，根據傅里葉傳熱定律，爆炸產物向內桶系統傳熱的熱流量ФN1(τ)正比于爆熱彈殼體的溫度場梯度：

3 量熱計傳熱模型

3.1 初期階段的傳熱模型

初期階段（絕熱法沒有初期階段）只有外桶傳熱與攪拌熱引起內桶系統的升溫，由于外桶傳熱與攪拌產熱速率低，因此，可認為此階段爆熱彈殼體、內桶水及其外殼同步升溫：

式中：τS為初期階段起始時刻；cT、mT分別為內桶水及其外殼的比熱容、質量，cN、mN分別為內桶系統的比熱容、質量。cN、mN滿足：

將式(5)、(12)代入式(4)，得到初期階段量熱計的傳熱學方程：

解式(14)，再結合式(12)，得到內桶水在初期階段的理論溫升曲線：

3.2 主期與末期階段的傳熱模型

由于局部升溫階段時間很短，因此全局升溫階段可近似為主期與末期階段。根據能量守恒定律，該階段爆炸產物傳熱、攪拌產熱、外桶系統傳熱共同引起爆熱彈殼體、內桶水及其外殼的溫升：

式中：PG為爆熱彈殼體熱容占內桶系統熱容的比例。PG為：

根據式(18)、(20)～(21)易證，λ2＜λ1＜0。

聯合式(4)～(6)、(9)、(11)、(16)，得到絕熱法在全局升溫階段的傳熱方程：

解式(23)，再結合式(9)、(10)、(16)，得到絕熱法在全局升溫階段內桶水理論溫升曲線：

4 快速系統辨識方法的提出

系統辨識方法是通過系統的輸入輸出來反推系統結構或辨識系統參數，參數辨識的目的是，通過參數估計得到的模型是在選定的模型類中最好的?；谠摲椒?，提出恒溫法測爆熱的快速系統辨識方法（以下簡稱快辨識法）的思路如下：首先，分別基于內桶水在初期階段、主末期階段的溫升曲線辨識參數kNW/(cNmN)、u1、λ1、λ2；然后，基于隔離易受干擾參數λ1的思路快速辨識修正溫升與爆熱。

4.1 各測量階段的參數辨識

量熱計如在初期階段出現故障，由于炸藥未起爆，尚有調整與修復的機會，因此，初期階段的全部測量數據可用于參數辨識；如在主末期階段出現故障，由于炸藥已起爆，很難有調整與修復的機會，因此，主末期階段只有故障前的測量數據可用于參數辨識。

考慮到內桶系統比熱容很大，而內外桶之間絕熱措施良好，因此有kNW/(cNmN)≈0，根據泰勒展開，內桶水在初期階段的理論溫升曲線（式(15)）可近似為線性方程：

對內桶水在初期階段的實測溫升曲線進行線性擬合，得到擬合方程：

式中：k1、k2、ε(τ)為擬合方程的一次項系數、零次項、擬合殘差。

對比式(26)～(27)可得，基于初期階段測量數據辨識參數kNW/(cNmN)為：

利用內桶水在主末期階段理論溫升曲線（式（19））擬合實測溫升曲線，以擬合殘差均方根最小化為目標來辨識系統參數，得到基于故障前測量數據辨識中間參數u1、λ1、λ2的目標函數為：

式中：τj為起爆后的第j個采樣時間點，k為起爆至測量時間τ 的采樣點個數。

注意，考慮局部升溫階段時間很短，因此，式(29)中將τ1L替換為τM對辨識精度影響可忽略。

由式(19)、(21)可知，λ1為內桶水溫升曲線慢變特征量，且kNW/cNmN≈0時，有λ1≈0，即內桶水溫升曲線慢變特征不顯著，因此，λ1的辨識結果易受到快變過程與溫度波動的干擾。多次模擬表明，基于目標函數（式(29)）最小化辨識λ1時，如測量時間較短且λ1無約束，則內桶水溫度辨識曲線的全局擬合精度較差，目標函數衰減速度很慢。因此，為了加快算法收斂速度，需對λ1的取值范圍進行約束。考慮到kNW＜＜kZT，再根據式(20)～(21)可得：

再考慮式(30)中kNW＜＜kZT、cZmZ＜＜cNmN，本文中給出一個相對寬松的約束條件：

因此，基于系統辨識的思想，本文中提出通過求解目標函數（式(29)）在約束條件（式(31)）下最小值的策略來辨識中間參數u1、λ1、λ2。

4.2 修正溫升與爆熱的快速辨識

首先介紹絕熱法修正溫升的辨識方法，然后借鑒其思路給出恒溫法修正溫升的辨識方法。

根據式(23)計算結果可以證明，絕熱法中爆炸產物向內桶系統的總傳熱量為：

再根據修正溫升的定義，結合式(25)，可以證明，絕熱法修正溫升在測量時刻τ 的辨識值為：

注意，本文中用(*)τ表示參數*在時刻τ 的辨識值。

因此，絕熱法辨識思路為，利用內桶水理論溫升曲線（式(24)）擬合實測溫升曲線辨識參數α，再根據式(33)、式(3)計算修正溫升、爆熱。由于參數α 反映溫升快變項的幅值，因此，α 收斂很快，爆熱辨識值收斂也很快。

而對恒溫法測量數據多次模擬表明，增加約束條件（式(31)）雖有利于提高溫升曲線的全局擬合精度，但λ1辨識值受溫度快變過程與溫度波動的影響仍會出現頻繁的振蕩現象，進而降低了修正溫升與爆熱辨識值的收斂速度；然而，爆炸產物放熱結束時，如內外桶溫差較大，則內桶水溫度將經歷一個長時間的慢變過程，慢變過程的累積影響也是不可忽略的，不可簡單地按絕熱過程（λ1=0）辨識內桶水的修正溫升。

綜上所述，基于隔離參數λ1振蕩影響，兼顧考慮恒溫法測量中普遍存在的內桶水溫度長時間慢變過程，參考絕熱法修正溫升辨識公式(33)，提出恒溫法修正溫升在測量時刻τ 快速辨識算法：

將式(34)代入式(3)，即可得到炸藥爆熱的辨識值。

4.3 辨識誤差分析

考慮kNW＜＜kZT，將式(22)代入式(34)可得，恒溫法修正溫升的辨識值近似收斂值為：

比較式(33)與式(35)可知，恒溫法與絕熱法修正溫升辨識值的收斂值一致。

根據爆熱的經典計算公式(1)～(2)，結合主末期階段的溫升曲線（式(19)），及起點溫度與溫升速率的理論辨識值，可證明經典方法修正溫升理論值為：

根據快辨識法與經典方法得到的修正溫升的相對誤差，來評價爆熱辨識值計算精度，令：

代入式(35)～(36)，計算可得：

由于TE(τM)＞＞TI(τM)、cNmN＞＞cZmZ，因此有P≈0，即快辨識法與經典方法的爆熱計算精度相當。

從式(34)～(35)、圖3可知，恒溫法修正溫升辨識值收斂特性不受參數λ1的直接影響，且參數u1反映溫度快變幅值，受干擾小。因此，綜合前述分析，爆熱的快速辨識方法具有較高的收斂速度、精度與穩定性。

圖3 實驗模擬的各參數辨識值Fig.3 Identified values of each parameter in an experiment simulation

5 實驗分析與討論

爆熱測量的經典算法與系統辨識算法雖基于相同的傳熱學模型，但誤差機制卻有較大的不同。經典算法本質上計算爆炸產物向內桶系統傳熱量的積分，其計算精度取決于爆炸產物放熱的完全性，因此，需要盡可能長的測量時間和盡可能小的外界累積干擾；而系統辨識算法本質上是基于爆炸產物向內桶系統傳熱特征量估計的傳熱量預測，其計算精度取決于傳熱模型的吻合性與特征量辨識的精確性。一般認為，測量時間足夠長時，經典算法得到的爆熱值能較好地收斂于真實值。本節以經典方法為基準，通過實驗數據檢驗快辨識法的性能，并提出可用于實際測量中收斂時間點判斷的實驗判據。

利用恒溫法測量設備進行爆熱測量。實驗測量真空條件下炸藥爆熱，實驗步驟參照GJB 772A—97中方法701.1，但為了降低測量誤差，初期、主期、末期3個階段分別延長為21、100、20 min?？紤]炸藥爆熱值主要在4～9 kJ/g 范圍，為了檢驗辨識方法的普適性，選取了8個基本均勻分布在該爆熱范圍內的典型炸藥（見表1）進行測量與模擬分析，其中，樣本1～4為高爆熱炸藥，樣本5～8為含TATB的炸藥。

表1 爆熱辨識值的收斂情況Table1 Convergent statusof system identify value of explosion heat

為了便于評價爆熱辨識值收斂效果，定義如下參數。

（1）爆熱辨識值相對誤差為：

式中：Q1為測量結束后經典方法的爆熱計算值，(Q1)τ為測量時刻τ 系統辨識方法的爆熱計算值。

（2）相對誤差上限水平、相對誤差平均水平為時刻τ 后爆熱辨識值相對誤差的最大值、平均值。

（3）極限收斂水平取相對誤差上限水平最小值。

選取表1中樣本1的實驗數據進行具體分析，如圖4所示。初始調溫階段內外桶溫差TC=0.552℃（TC取決于攪拌產熱速率），外桶升溫結束后的溫度TW=27.019℃。炸藥于85 min 起爆，起爆后，內桶水溫度經歷了一個快速上升的過程，理論上講，測量時間足夠長時，內桶水溫度將達到TW+TC。

圖4 內外桶溫度隨采樣時間的變化Fig.4 Variationsof inner and outer barrel temperature with sampling time

樣本1各參數收斂情況如圖3所示，參數λ2、u1收斂穩定，參數λ1振蕩明顯，但修正溫升Δt有效地隔離了參數λ1的振蕩影響，起爆40 min 后變化量僅為0.11℃，從而檢驗了本算法的收斂穩定性。

將各參數辨識值代入式(19)，即得到主末期內桶水溫升的擬合曲線，擬合殘差均方根變化情況如圖5所示。隨著辨識時間的增加，擬合殘差均方根迅速而穩定地減小，在起爆后27 min 即達到0.2℃量級，從而檢驗了本算法對內桶水溫升變化的預測能力。

圖5 內桶水主末期溫度曲線的擬合殘差均方根Fig.5 RMSof fitting residualsof inner barrel water temperature curve at main and end stage

樣品1的模擬表明，以測量結束后經典方法爆熱計算值7.765 kJ/g 為參照，爆熱辨識值在起爆后40、43、51 min 快速穩定地收斂至3.5%（7.498 kJ/g）、3%（7.537 kJ/g）、2%（7.613 kJ/g）的相對誤差上限水平，并最終達到了0.024 5%（7.764 kJ/g）的極限收斂水平，如圖6～7所示，從而檢驗了本算法的收斂速度與精度。

圖6 爆熱系統辨識值和相對誤差隨測量時間的收斂Fig.6 Convergence rule of system identify valueand relativeerror of explosive heat with measurement time

圖7 相對誤差上限水平和平均水平隨測量時間的變化Fig.7 Variationsof upper limit level and average level of relative error with measurement time

對炸藥的恒溫法爆熱測量數據應用快辨識法進行模擬分析，各炸藥爆熱的經典值、辨識值收斂至3.5%相對誤差上限水平所需最小時間、極限收斂水平等參數計算結果，見表1。

由表1可見：快辨識法得到的爆熱值能快速地（一般不超過40 min，即主末期總時間的1/3）收斂至經典值3.5%的相對誤差上限水平內；除樣品5極限收斂水平為3.221 6%，其他樣本還能在起爆后50 min內進入2%的相對誤差上限水平。這驗證了本算法的普適性。根據式(38)，樣品5收斂相對較差的原因可能是炸藥爆熱值較低，對應爆溫也較低，導致收斂精度略差。

為了提高本方法的可操作性，給出測量過程中爆熱辨識值的收斂情況，根據修正溫升辨識值收斂時的慢變特性，并且研究發現，線性擬合測量時刻τ 及前14 min（共15個測量點）內的修正溫升辨識值-時間曲線的直線斜率絕對值（記為J1）、擬合殘差均方根（記為J2）特征量能較好地反映時刻τ 的數據變化速率與波動性。因此，本文中利用統計方法分析了各樣本的這兩個特征量的分布規律，得到了實驗判據：如J1≤0.008℃/min，且J2≤0.006℃，則爆熱辨識值在測量時刻τ 的相對誤差上限水平小于3.5%。

模擬結果表明，本實驗判據得到的收斂時間略大于理論收斂時間且判斷誤差不超過20 min，具有較好的可靠性。但考慮本實驗判據所統計的樣本量較小，還需要開展大樣本實驗進一步完善。

6 結 論

快辨識法在確保收斂精度的前提下，有效地壓縮了測量時間，降低了系統故障導致測量失敗的風險；或者說，當出現系統故障時，經典測量方法失效的情況下，快辨識法仍能得到一個不算差的估算結果。因此，快辨識法可作為經典方法的備用方案。通過理論與實驗研究，主要獲得以下結論。

（1）快辨識法只辨識中間參數kNW/(cNmN)、u1、λ1、λ2，而不辨識目標參數TE(τM)、kZT/(cNmN)、kZT/(cZmZ)，避免了參數λ1振蕩通過目標參數傳導至修正溫升辨識值，算法的收斂速度與穩定性有了顯著提高，但由于收斂值與經典值存在偏差，收斂精度會略有損失。

（2）當kNW/(cNmN)=0時，快辨識法的修正溫升表達式退化為絕熱法修正溫升表達式，因此，快辨識法可拓展至絕熱法的爆熱辨識。

（3）對樣本1 模擬結果表明，中間參數λ2、u1收斂穩定，中間參數λ1雖存在明顯振蕩，但修正溫升Δt有效地隔離了λ1的振蕩影響，并在起爆后40 min 收斂于0.11℃波動水平，驗證了快辨識法的收斂穩定性。

（4）對樣本1模擬結果表明，對內桶水在主末期溫升曲線擬合殘差均方根在起爆后27 min 即達到0.2℃量級，爆熱辨識值在起爆后40 min 即收斂至3.5%的相對誤差上限水平，并最終達到了0.0245%的極限收斂水平，驗證了快辨識法的收斂速度與精度。

（5）對爆熱值分布在4～9 kJ/g 的8 個炸藥樣本模擬結果表明，爆熱辨識值可快速地（一般不超過40 min，即主末期總時間的1/3）收斂于經典值3.5%的相對誤差水平內，驗證了快辨識法的普適性。

（6）提出了在測量過程中實時判斷爆熱辨識值收斂情況的實驗判據，且收斂時刻判斷誤差不超過20 min。