沈 飛，王 輝，張 皋

（西安近代化學(xué)研究所，陜西 西安 710065）

1 引言

炸藥爆轟釋能是一個(gè)復(fù)雜的物理化學(xué)過(guò)程，其量化研究大多基于試驗(yàn)觀測(cè)數(shù)據(jù)進(jìn)行唯像建模。由于不同炸藥的爆轟過(guò)程差異較大，試驗(yàn)觀測(cè)數(shù)據(jù)的規(guī)律性表征常采用特定的多重累加型函數(shù)（由多個(gè)形式相同的多項(xiàng)式疊加組成）進(jìn)行回歸擬合，如炸藥水下爆炸近場(chǎng)的沖擊波傳播距離隨時(shí)間變化的曲線［1，2］，穩(wěn)態(tài)爆轟的波陣面形狀曲線［3，4］，圓筒試驗(yàn)中的銅管徑向膨脹位移隨時(shí)間變化的曲線等［5，6］。為了進(jìn)一步獲得評(píng)價(jià)炸藥爆轟性能的物理量，常需要對(duì)這些擬合曲線進(jìn)行微分處理，因而建立擬合過(guò)程不確定度與其它物理量不確定度之間的關(guān)系，這對(duì)于整個(gè)爆轟模型測(cè)量或計(jì)算精度的分析十分重要。為了系統(tǒng)地建立該類模型的不確定度分析方法，本研究針對(duì)圓筒試驗(yàn)?zāi)Ｐ瓦M(jìn)行了相應(yīng)的不確定度分析。

圓筒試驗(yàn)是研究炸藥驅(qū)動(dòng)能力的經(jīng)典試驗(yàn)，通過(guò)狹縫高速掃描攝影法獲取定常滑移爆轟驅(qū)動(dòng)下銅管壁的徑向膨脹位移隨時(shí)間的變化曲線，并采用特定的多重累加函數(shù)進(jìn)行擬合，然后進(jìn)一步計(jì)算出銅管膨脹速度、比動(dòng)能等物理量［6，7］。該試驗(yàn)及分析模型主要包括膨脹位移-時(shí)間曲線的數(shù)據(jù)獲取、擬合、微分3個(gè)步驟，一些文獻(xiàn)［5，7］分析了位移或時(shí)間數(shù)據(jù)點(diǎn)測(cè)量不確定度以及曲線擬合參數(shù)不確定度的基本計(jì)算方法。但現(xiàn)有文獻(xiàn)中并未建立各物理量不確定度之間的關(guān)系，且所采用的擬合參數(shù)不確定度計(jì)算方法不能適用于該多重累加擬合函數(shù)，因此未能建立銅管膨脹速度、比動(dòng)能等物理量不確定度的分析方法。

本研究通過(guò)在膨脹位移表達(dá)式中引入估計(jì)值為1的無(wú)量綱參量，獲得了膨脹位移和膨脹速度兩者測(cè)量不確定度的關(guān)系式，并修正了多重累加函數(shù)的擬合參數(shù)不確定度分析模型，建立了較為系統(tǒng)的圓筒試驗(yàn)不確定度評(píng)定及計(jì)算方法。

2 基本原理及原始數(shù)據(jù)的測(cè)量不確定度

2.1 圓筒試驗(yàn)基本原理

炸藥圓筒試驗(yàn)的原理如圖1所示，將炸藥裝填在標(biāo)準(zhǔn)尺寸的銅管內(nèi)，并從銅管的一端起爆炸藥，采用狹縫掃描相機(jī)觀測(cè)銅管的某一截面。當(dāng)相機(jī)的轉(zhuǎn)鏡高速旋轉(zhuǎn)時(shí)，在膠片的不同位置可記錄銅管在不同時(shí)刻的膨脹狀態(tài)。

圖1 圓筒試驗(yàn)原理示意圖Fig.1 Schematic diagram of cylinder test

圓筒試驗(yàn)的原始數(shù)據(jù)為膠片中的圖像邊界，如圖2所示。其中，橫向?yàn)檗D(zhuǎn)鏡掃描光路的線位移，其除以轉(zhuǎn)鏡的線速度即為時(shí)間值；縱向可反映銅管直徑的變化。

圖2 圓筒膨脹圖像的邊界參量示意圖Fig.2 Schematic diagram of boundary parameters of cylinder expansion image

2.2 銅管徑向位移的測(cè)量不確定度

采用高精度掃描儀將膠片圖像轉(zhuǎn)換為電子圖像，然后由專用判讀軟件進(jìn)行邊界識(shí)別（銅管徑向的像位移（yj－y0），對(duì)應(yīng)的橫向距離為（xj－x0），其識(shí)別的判據(jù)為圖像灰度梯度是否達(dá)到極值。

根據(jù)相應(yīng)的算法分析［8，9］及重復(fù)性驗(yàn)證獲知，對(duì)于2 400 dpi分辨率的圖像，采用一階算法時(shí)，邊界線的識(shí)別偏差為0～2個(gè)像素（pixel），即0～0.02 mm，則邊界識(shí)別的不確定度即為0.02 mm。由于確定兩點(diǎn)距離時(shí)需判定兩個(gè)邊界，則u（yj－y0）＝u（xj－x0）≈0.03 mm。

綜上分析可以發(fā)現(xiàn)，銅管外表面的實(shí)際徑向位移Δrej＝（yj－y0）／β的測(cè)量不確定度來(lái)源于圖像邊界識(shí)別和放大比測(cè)量的隨機(jī)誤差。然而將圖像放大比作為一個(gè)固定值輸入到每一個(gè)數(shù)據(jù)點(diǎn)時(shí)，將產(chǎn)生系統(tǒng)誤差，并影響后續(xù)的擬合、求導(dǎo)等過(guò)程。為了便于對(duì)該問(wèn)題進(jìn)行分析，可將銅管實(shí)際徑向位移的表達(dá)式修改為：

即引入?yún)⒘喀?，其估?jì)值為1，u（λ）≈u（β）／β，對(duì)于?25 mm的標(biāo)準(zhǔn)銅管，則u（λ）≈3.5×10－3。

2.3 時(shí)間的測(cè)量不確定度

由于銅管膨脹時(shí)間為tj＝（xj－x0）／v，其中，v為相機(jī)的掃描線速度，則

從式（2）可以看出，對(duì)于同一發(fā)試驗(yàn)，僅有（xj－x0）為變量，u（tj）隨（xj－x0）單調(diào)遞增。文獻(xiàn)［7］研究認(rèn)為，u（v）＝0.001v，則對(duì)于?25 mm的標(biāo)準(zhǔn)圓筒試驗(yàn)，當(dāng)v＝3 mm／μs時(shí)，式（2）可簡(jiǎn)化為：

3 位移曲線擬合參數(shù)的不確定度

3.1 曲線擬合的基本步驟

目前GJB 8381—2015［5］及國(guó)外公開(kāi)報(bào)道的資料［6］中，均首先將銅管外表面膨脹軌跡的判讀點(diǎn)（tj，Δrej）按照公式（4）轉(zhuǎn)換為銅管質(zhì)量中心面處的數(shù)據(jù)點(diǎn)（tj，Δrmj），

式中：re0和ri0分別為銅管內(nèi)外表面的初始半徑。

由公式（4）和不確定度傳遞規(guī)律可知，u（Δrej）≈u（Δrmj），然后按照公式（5）的形式，利用最小二乘法對(duì)數(shù)據(jù)點(diǎn)（tj，Δrmj）進(jìn)行擬合，從而獲得a1、b1、a2、b2、t0等5個(gè)擬合參數(shù)。此外，公式（5）為典型多重累加函數(shù)，為了便于表示，采用f1和f2分別表示兩個(gè)多項(xiàng)式。

擬合曲線的實(shí)驗(yàn)標(biāo)準(zhǔn)偏差可作為擬合曲線的標(biāo)準(zhǔn)不確定度［10］，即

式中：Δrmj和j分別為時(shí)間tj所對(duì)應(yīng)的徑向位移實(shí)驗(yàn)值和擬合值；（n－p）為自由度，n為實(shí)驗(yàn)點(diǎn)的總數(shù)，p為擬合參數(shù)的數(shù)量。采用決定系數(shù)R2值可作為擬合符合程度的度量。

以TNT和JO-159炸藥的某一發(fā)?25 mm圓筒試驗(yàn)數(shù)據(jù)為例進(jìn)行分析，其結(jié)果見(jiàn)表1，由決定系數(shù)R2值可看出擬合精度較高。

表1 銅管膨脹位移曲線擬合參數(shù)Tab.1 Curve-fitting parameters of the expansion displacement of the cylinder wall

3.2 多重累加函數(shù)擬合參數(shù)不確定度的模型修正

由于公式（5）為典型弱非線性函數(shù)，尤其是當(dāng)t不斷增大后，Δrm近似于線性增長(zhǎng)，因此，計(jì)算擬合參數(shù)的不確定度時(shí)，可采用泰勒一階展開(kāi)式對(duì)其進(jìn)行近似表示［10］：

式中：y和γ分別為數(shù)據(jù)點(diǎn)Δrmj及殘差組成的向量；β為p個(gè)擬合參數(shù)所組成的向量；Aij＝

從而可進(jìn)一步計(jì)算出擬合參數(shù)的方差—協(xié)方差矩陣C＝（ATA）－1，矩陣ATA的元素形式可表示為

擬合參數(shù)βi的標(biāo)準(zhǔn)不確定度為

不同擬合參數(shù)之間協(xié)方差為

然而，若將公式（5）中的5個(gè)擬合參數(shù)作為估計(jì)量，即β＝［a1b1a2b2t0］T，則采用上述方法計(jì)算不確定度時(shí)，銅管膨脹軌跡和公式（5）的固有特征會(huì)導(dǎo)致兩類現(xiàn)象出現(xiàn)，從而造成計(jì)算錯(cuò)誤。第一類現(xiàn)象，若b1或b2的值較大（為便于描述，這里假定b1的值較大），f1≈a1（t＋t0），膨脹曲線的非線性特征主要由f2描述，則由公式（5）可推導(dǎo)出

將其代入公式（8）可知，矩陣ATA的第2列及第2行的元素均近似為0，使得該矩陣為奇異矩陣，難以求解其逆矩陣。第二類現(xiàn)象，若b1或b2較為接近，近似于公式（5）中j＝1的情況，使得?f／?a1≈?f／?a2，從而造成ATA矩陣的兩行或兩列數(shù)據(jù)相近，即近似為奇異矩陣，仍難以求解其逆矩陣。

對(duì)于復(fù)雜非線性函數(shù)造成的不確定度計(jì)算困難，雖然從蒙特卡羅方法中可能會(huì)尋求出相應(yīng)的解決途徑［11～13］，但在后續(xù)的求導(dǎo)等計(jì)算過(guò)程中將會(huì)出現(xiàn)新的問(wèn)題，因此，本研究中仍基于傳統(tǒng)GUM不確定度評(píng)定方法進(jìn)行分析。結(jié)合公式（5）所示的多重累加回歸模型的特點(diǎn)可以發(fā)現(xiàn)，雖然f1和f2一般是同時(shí)進(jìn)行擬合計(jì)算的，但本質(zhì)上這兩個(gè)多項(xiàng)式的擬合過(guò)程可以分步進(jìn)行。第一步，采用Δrm＝f1（t）對(duì)實(shí)驗(yàn)數(shù)據(jù)點(diǎn)（tj，Δrmj）進(jìn)行擬合，可獲得a1、b1、t0這3個(gè)參數(shù)值，并計(jì)算出每個(gè)數(shù)據(jù)點(diǎn)的殘差（tj，Δrmjf1（tj））；第二步，將（tj，Δrmj－f1（tj））數(shù)據(jù)點(diǎn)代入Δrm－f1（t）＝f2（t）進(jìn)行擬合，獲得a2、b2兩個(gè)參數(shù)值?？梢钥闯?，最終殘差值由第二步擬合的精度決定，因此，可認(rèn)為a1、b1、t0的值為真值，而a2、b2的值為估計(jì)值，并分析這兩個(gè)擬合參數(shù)的不確定度。設(shè)β＝［a2b2］T，則ATA矩陣的具體形式為

式中：

需要說(shuō)明的是，選擇估計(jì)值參量時(shí)，a1、b1和a2、b2可以互換。

根據(jù)式（9）、式（10）、式（12），并結(jié)合表1所列數(shù)據(jù)，可計(jì)算出兩種炸藥的擬合參數(shù)的不確定度，具體數(shù)值見(jiàn)表2。

表2 擬合參數(shù)的標(biāo)準(zhǔn)不確定度Tab.2 Uncertainty evaluation of curve-fitting parameters

4 銅管速度及比動(dòng)能的不確定度

4.1 銅管徑向膨脹速度的測(cè)量不確定度

將公式（5）對(duì)時(shí)間t求導(dǎo)，并結(jié)合公式（1），可得到銅管的徑向膨脹速度vm：

則vm的標(biāo)準(zhǔn)不確定度為

式中：g1、g2、g3分別表示圖像放大比、擬合參數(shù)、膨脹時(shí)間的相關(guān)多項(xiàng)式。

為了弄清這3類因素對(duì)u（vm）的影響程度，可基于表1中的數(shù)據(jù)對(duì)ω（λ）＝g1／u2（vm）、ω（a2，b2）＝g2／u2（vm）、ω（t）＝g3／u2（vm）的變化規(guī)律進(jìn)行分析，結(jié)果如圖3所示。

圖3 不同因素對(duì)銅管徑向速度不確定度的影響規(guī)律Fig.3 The influence of different factors on the uncertainty of the radial velocity of copper pipe

由圖3可以看出：對(duì)于TNT和JO-159兩種不同作功能力的炸藥，u（λ）（即圖像放大比的相對(duì)標(biāo)準(zhǔn)不確定度u（β）／β）對(duì)u（vm）的影響均最大，且隨著銅管膨脹位移的增大，其影響程度逐漸增大，而擬合參數(shù)和膨脹時(shí)間的不確定度對(duì)其影響程度逐步減弱，尤其是在膨脹后期，膨脹時(shí)間不確定度的影響基本可以忽略；對(duì)于JO-159炸藥，u（λ）對(duì)u（vm）的影響程度高于TNT，而擬合參數(shù)不確定度對(duì)u（vm）的影響程度低于TNT。因此，在圓筒試驗(yàn)測(cè)量過(guò)程中，對(duì)圖像放大比的測(cè)量精度應(yīng)更為重視。

4.2 圓筒比動(dòng)能的測(cè)量不確定度

獲得圓筒比動(dòng)能時(shí)，首先需要計(jì)算出銅管質(zhì)量中心面的質(zhì)點(diǎn)速度vs，即

式中：D為炸藥在銅管中的爆速，可通過(guò)銅管首、尾端設(shè)置的電探針同步測(cè)量，TNT和JO-159的爆速分別為6.830 mm／μs和8.760 mm／μs。

根據(jù)國(guó)軍標(biāo)GJB 772A-1997［14］可計(jì)算出TNT和JO-159爆速的測(cè)量不確定度u（D）分別為0.010 mm／μs和0.016 mm／μs。由公式（15）可獲得vs的標(biāo)準(zhǔn)不確定度計(jì)算式為

式中：

取包含因子k＝2，置信水平約0.95，則圓筒比動(dòng)能E的相對(duì)擴(kuò)展不確定度為：

基于上述分析模型，可對(duì)TNT和JO-159炸藥的?25 mm圓筒徑向膨脹速度vm和比動(dòng)能E的測(cè)量不確定進(jìn)行分析，結(jié)果如圖4所示。

由圖4（a）可以看出：對(duì)于這2種炸藥，u（E）與E近似呈線性關(guān)系，在銅管膨脹后期，u（E）隨E的增長(zhǎng)速率逐漸變緩；在銅管膨脹過(guò)程中（忽略膨脹初期），JO-159炸藥的E及u（E）均高于TNT。圖4（b）顯示了比動(dòng)能E的相對(duì)擴(kuò)展不確定度U（E）／E（k＝2）的變化曲線，由圖可以看出：隨著E的增大，U（E）／E均呈現(xiàn)先增大（主要在膨脹早期）后減小的趨勢(shì)；JO-159的U（E）／E值變化范圍為1.60%至1.75%，而TNT的U（E）／E值變化范圍為1.67%至1.93%。因此，對(duì)于作功能力較強(qiáng)的炸藥，其U（E）／E值低于作功能力低的炸藥。此外，由于這兩種炸藥具有較強(qiáng)的代表性，可以預(yù)測(cè)，對(duì)于大多數(shù)高能炸藥的?25 mm圓筒試驗(yàn)，其圓筒比動(dòng)能的相對(duì)擴(kuò)展不確定度不超過(guò)2%。

圖4 圓筒比動(dòng)能的不確定度變化曲線Fig.4 Uncertainty curve of specific kinetic energy of cylinder

5 結(jié)論

（1）對(duì)于含多重累加擬合函數(shù)的爆轟模型，分析其擬合參數(shù)的不確定度時(shí)，應(yīng)選擇影響最終擬合殘差值的部分參數(shù)作為估計(jì)值，從而可避免因出現(xiàn)奇異矩陣而導(dǎo)致的計(jì)算錯(cuò)誤。

（2）銅管膨脹速度及比動(dòng)能測(cè)量不確定度的最大影響因素是圖像放大比的相對(duì)不確定度，其次為膨脹曲線擬合參數(shù)的不確定度，而膨脹時(shí)間的不確定度對(duì)其影響最小。

（3）隨著圓筒比動(dòng)能的增大，其標(biāo)準(zhǔn)不確定度近似呈線性增長(zhǎng)趨勢(shì)，而相對(duì)不確定度則呈現(xiàn)先增大（主要在膨脹早期）后減小的趨勢(shì)；且對(duì)于作功能量較強(qiáng)的JO-159炸藥，其圓筒比動(dòng)能的相對(duì)不確定度明顯低于TNT。此外，根據(jù)這2種典型炸藥的分析結(jié)果可以預(yù)測(cè)，對(duì)于大多數(shù)高能炸藥的?25 mm圓筒試驗(yàn)，其圓筒比動(dòng)能的相對(duì)擴(kuò)展不確定度（k＝2）不超過(guò)2%。