張 龍, 張繼軍, 趙建偉, 張寶國, 張東亮
(西北核技術(shù)研究院,陜西 西安 710024)
在密閉空間爆炸實驗中,準(zhǔn)確測量爆后氣體的溫度變化歷程對于分析爆炸實施效果和氣體擴(kuò)散規(guī)律具有十分重要的意義[1,2]。由于爆炸過程具有極強(qiáng)的瞬時性和動態(tài)性,因此對爆后氣體溫度參數(shù)的測量屬于動態(tài)測量范疇。與靜態(tài)測量相比,動態(tài)測量具有波動性、隨機(jī)性等特征,且動態(tài)測量過程容易受到多種未知隨機(jī)誤差的干擾[3]。為獲得更加準(zhǔn)確的爆后氣體溫度參數(shù),希望在現(xiàn)有傳感器的基礎(chǔ)上通過數(shù)學(xué)方法減小被測量動態(tài)測量結(jié)果的隨機(jī)性誤差,并對其估計真值和動態(tài)測量不確定度做出準(zhǔn)確評定。
在動態(tài)測量數(shù)據(jù)處理方法中,非統(tǒng)計理論表現(xiàn)出較強(qiáng)的優(yōu)勢[4],灰色系統(tǒng)理論[5]、Bootstrap方法[6]等算法在動態(tài)測量數(shù)據(jù)的不確定度評定中得到了廣泛應(yīng)用。在測量數(shù)據(jù)序列概率分布未知的條件下,灰色模型GM(1,1)可以準(zhǔn)確預(yù)測被測量瞬時值的大小,但難以估計其置信區(qū)間,因而無法評估在給定置信概率下的動態(tài)測量不確定度[7,8];Bootstrap方法可以通過自助再抽樣操作模擬動態(tài)測量數(shù)據(jù)序列的概率分布,并估計其置信區(qū)間,但計算機(jī)仿真結(jié)果表明,自助再抽樣方法引入了附加的不確定度分量[9],并且由Bootstrap方法獲得的估計區(qū)間小于原始數(shù)據(jù)的實際分布區(qū)間,從而增大了動態(tài)不確定度的估計誤差。
綜上所述,Bootstrap方法和灰色模型GM(1,1)均無法對密閉空間爆炸溫度的動態(tài)不確定度做出準(zhǔn)確評價。本文融合灰色模型GM(1,1)、Bootstrap方法以及不確定度評定理論,將灰自助模型GBM(1,1)運用到密閉空間爆炸溫度動態(tài)測量不確定度的評定中,并選取估計真值、估計區(qū)間、區(qū)間估計可靠度、動態(tài)測量不確定度和平均不確定度等參數(shù)表征其估計結(jié)果。實驗結(jié)果表明,GBM(1,1)模型融合了灰色模型GM(1,1)和Bootstrap方法的優(yōu)勢,能夠準(zhǔn)確模擬動態(tài)測量數(shù)據(jù)序列的概率分布,并實時跟蹤被測量瞬時值的變化趨勢。相比于灰色模型GM(1,1)和Bootstrap方法,灰自助模型GBM(1,1)具有更高的真值估計精度和區(qū)間估計可靠度,其區(qū)間估計可靠度高于90%,估計區(qū)間能夠更加完整地包絡(luò)被測量的動態(tài)波動范圍。在密閉空間爆炸實驗中,該方法能夠?qū)Ρ覝囟鹊膭討B(tài)測量不確定度做出準(zhǔn)確評估,具有較高的實用價值。
設(shè)密閉爆室內(nèi)氣體的溫度參數(shù)為隨機(jī)變量,通過數(shù)據(jù)采集系統(tǒng)可獲得被測量的動態(tài)測量數(shù)據(jù)序列為:
X={x(t)}, (t=1,2,…)
(1)
式中:x(t)表示t時刻的測量數(shù)據(jù)。從X中抽取t時刻之前的m個數(shù)據(jù)(包括t時刻的數(shù)據(jù)),構(gòu)成時刻t的動態(tài)評估子序列Xm:
Xm={xm(n)},
(n=t-m+1,t-m+2,…,t)
(2)
根據(jù)Bootstrap方法,從Xm中有放回地等概率抽取1個數(shù)據(jù),共抽取m次,得到第一個Bootstrap樣本。重復(fù)上述方法B次,得到B個Bootstrap再抽樣樣本,用向量表示為:
YBoot=(Y1,Y2,…,Yb,),
(b=1,2,…,B)
(3)
式中Yb為第b個Bootstrap再抽樣樣本,且有
Yb={yb(n)},
(n=t-m+1,t-m+2,…,t)
(4)
式中:yb(n)為Yb中第n個Bootstrap再抽樣數(shù)據(jù)。
在灰色系統(tǒng)理論中,累加生成是重要的數(shù)據(jù)處理方法。任意非負(fù)波動數(shù)列經(jīng)累加生成操作后均可轉(zhuǎn)化為遞增數(shù)列,從而削弱了原始數(shù)據(jù)序列的隨機(jī)性,突出其變化趨勢,有助于發(fā)現(xiàn)數(shù)據(jù)內(nèi)在的變化規(guī)律[10]。根據(jù)灰色模型GM(1,1),對Yb做一階累加生成運算,其生成序列表示為:
(5)
設(shè)均值生成序列為:
Mb={mb(n)}={0.5xb(n)+0.5xb(n-1)},
(n=t-m+1,t-m+2,…,t)
(6)
在初始條件xb(1)=yb(1)下,累加生成序列的預(yù)測值為:
(7)
式中c1和c2可由下述公式求得:
(c1,c2)T=(DTD)-1DT(Yb)T,
(n=t-m+2,t-m+3,…,t)
(8)
D=(-Mb,I)T
(9)
I=(1,1,…,1)
(10)
由累減生成方法可求得t+1時刻的預(yù)測值為:
(11)
在t+1時刻,有B個數(shù)據(jù),可構(gòu)成序列:
(12)
t時刻xm的頻率函數(shù)可表示為:
ft+1=ft+1(xm)
(13)
式中ft+1稱為GBM(1,1)模型的概率密度函數(shù)。
(1) 估計真值X0
被測量在t時刻的估計真值可由瞬時加權(quán)真值表示:
(14)
式中:X0為用最大概率值表示的t時刻的最終解;Q表示將Xt+1分為Q組;q為第q組;xmq為第q組數(shù)據(jù)的中值;f(t+1)q為對應(yīng)于xmq的灰自助概率。
(2) 估計區(qū)間[XL,XU]
給定置信概率為P,則t時刻被測量真值的估計區(qū)間可表示為:
[XL,XU]=[Xα/2,X1-α/2]
(15)
式中:a為顯著性水平,a∈[0,1];置信概率P=1-a;XL、XU分別表示估計區(qū)間的下限值和上限值;Xa/2為概率是a/2時變量xm對應(yīng)的值;X1-a/2為概率是1-a/2時變量xm對應(yīng)的值。
(3) 動態(tài)不確定度U
t+1時刻被測量的動態(tài)不確定度可表示為:
U=XU-XL
(16)
(4) 區(qū)間估計可靠度PB
設(shè)測量過程中被測量的采樣次數(shù)為T,若有h個數(shù)據(jù)位于估計區(qū)間[XL,XU]之外,則區(qū)間估計可靠度可表示為:
(17)
(5) 平均不確定度Umean
平均不確定度Umean是一個統(tǒng)計量,可以作為被測量隨機(jī)性的評價指標(biāo)。對平均不確定度Umean最理想的評價方法是在PB=100%的條件下,Umean取最小值。Umean越小,則被測量的波動范圍越小。
為準(zhǔn)確評價動態(tài)不確定度的估計效果,定義平均不確定度Umean為:
(18)
由式(18)可知,平均不確定度Umean是m,B,P的函數(shù),為使Umean獲得最理想的評估結(jié)果,必須選擇合適的參數(shù)m,B,P。由計算機(jī)仿真可得,在置信概率為99.7%的前提下,參數(shù)m,B與估計可靠度PB的關(guān)系如表1所示。
表1 參數(shù)m,B與估計可靠度PB的關(guān)系Tab.1 Relationship between the parameters m, B and the estimation reliability PB
由表1可知:若m值較大,則B的值也應(yīng)當(dāng)較大,算法的運算量也隨之增大,由此導(dǎo)致跟蹤被測量變化趨勢的實時性減弱[11];若m值較小,則B也應(yīng)選取較小值,因選取的數(shù)據(jù)量較小,估計結(jié)果的穩(wěn)定性將會變差[12]。在實際問題中,需結(jié)合測量準(zhǔn)確度要求和數(shù)據(jù)序列的樣本量選擇合適的m,B。在本文中,因數(shù)據(jù)序列的樣本量較大且精度要求較高,應(yīng)選取較大的m和B值。本文中選取m=7,B=1 000。
密閉空間爆炸實驗數(shù)據(jù)采集系統(tǒng)的基本結(jié)構(gòu)如圖1所示,主要包括K型熱電偶溫度傳感器、信號調(diào)理模塊、數(shù)據(jù)采集模塊和采集控制主機(jī)等。所述數(shù)據(jù)采集模塊選用Advantech PCI-1710數(shù)據(jù)采集卡,設(shè)定采樣頻率為1 kHz,選取爆炸零時后300 s的溫度數(shù)據(jù)作為實驗分析數(shù)據(jù)。
圖1 數(shù)據(jù)采集系統(tǒng)基本結(jié)構(gòu)Fig.1 Basic structure of data acquisition system
根據(jù)前期實驗數(shù)據(jù)和傳感器性能指標(biāo),選定鎧裝K型熱電偶作為溫度傳感器,其基本性能指標(biāo)為:溫度測量范圍0~1 250 ℃,準(zhǔn)確度±2.2 ℃,響應(yīng)時間0.73 s;實驗所需傳感器均經(jīng)過計量檢定,滿足測量精度要求。
實驗結(jié)束后,分別采用灰色模型GM(1,1)、Bootstrap方法和GBM(1,1)模型對測量所得數(shù)據(jù)進(jìn)行真值估計,對比分析3種方法的真值估計結(jié)果。在給定置信概率為99%,95%,90%的前提下,分別采用Bootstrap方法和GBM(1,1)模型對測量所得數(shù)據(jù)進(jìn)行區(qū)間估計,比較同一置信概率下2種方法的區(qū)間估計可靠度PB和平均不確定度Umean。
3.3.1 真值估計
選取m=7,B=1 000,Q=10,圖2(a)~圖2(c)分別表示灰色模型GM(1,1)、Bootstrap方法和GBM(1,1)模型對爆炸溫度的真值估計結(jié)果。圖中相對時間是相對于實驗起爆時間。
圖2 3種模型的真值估計結(jié)果Fig.2 True value estimation results of the three methods
由圖2可知,在密閉爆室溫度參數(shù)的真值估計中,GBM(1,1)模型的真值估計結(jié)果優(yōu)于Bootstrap方法和灰色模型GM(1,1)。分析認(rèn)為,GBM(1,1)模型融合了灰色模型GM(1,1)和Bootstrap方法的優(yōu)勢,不僅具有對數(shù)據(jù)樣本的擴(kuò)展能力,而且具備對數(shù)據(jù)的預(yù)測機(jī)制,因此可以更加準(zhǔn)確地模擬測量數(shù)據(jù)列的變化趨勢。
為對上述3種模型的真值估計結(jié)果做出更直觀的對比,對3種模型估計誤差的分布區(qū)間、平均值、標(biāo)準(zhǔn)差和最大相對誤差進(jìn)行比較。首先,剔除估計值中的粗大誤差,先用萊以特準(zhǔn)則,以估計誤差序列的3倍標(biāo)準(zhǔn)差為取舍標(biāo)準(zhǔn),判定并剔除估計值中的粗大誤差[13];其次,研究估計誤差的分布規(guī)律,經(jīng)驗證,3種模型的估計誤差均服從正態(tài)分布,其中GBM(1,1)模型估計誤差的正態(tài)分布P-P圖如圖3所示。
圖3 GBM(1,1)模型估計誤差的正態(tài)分布P-P圖Fig.3 The P-P chart of normal distribution for estimation errors of GBM (1,1) model
3種模型估計誤差的分布區(qū)間、平均值、標(biāo)準(zhǔn)差和最大相對誤差的計算結(jié)果如表2所示:由于3種模型在各個采樣點處的估計誤差均近似服從正態(tài)分布,故其均值約為零;GBM(1,1)模型的估計誤差分布區(qū)間、估計誤差的標(biāo)準(zhǔn)差和最大相對誤差均小于灰色模型GM(1,1)和Bootstrap方法,說明GBM(1,1)模型的估計誤差較小,且具有較好的重復(fù)性和穩(wěn)定性。因此認(rèn)為GBM(1,1)模型的真值估計結(jié)果優(yōu)于灰色模型GM(1,1)和Bootstrap方法。
表2 估計誤差的分布區(qū)間和平均值Tab.2 Distribution interval and mean value of estimation error
3.3.2 區(qū)間估計
圖4給出了不同置信概率下Bootstrap方法對爆室溫度動態(tài)測量數(shù)據(jù)的區(qū)間估計結(jié)果(m=7,B=1 000,Q=10)。圖5為對應(yīng)置信概率下GBM(1,1)模型的區(qū)間估計結(jié)果,其中,圖5(a)~圖5(c)分別對應(yīng)置信概率為99%,95%,90%。
圖4 不同置信概率下Bootstrap方法的估計區(qū)間Fig.4 Estimation interval of the Bootstrap method with different confidence levels
圖5 不同置信概率下GBM(1,1)模型的估計區(qū)間Fig.5 Estimation interval of the GBM (1,1) model with different confidence levels
對比圖4和圖5可知,在相同置信概率下,GBM(1,1)模型的估計區(qū)間寬于Bootstrap方法,能夠更加完整地包絡(luò)被測量的動態(tài)波動范圍。在置信概率不同的情況下,P越小,估計區(qū)間越窄,對被測量的動態(tài)波動范圍包絡(luò)越緊密,但區(qū)間估計的可靠度會隨之降低。
為了更直觀地對比Bootstrap方法和GBM(1,1)模型區(qū)間估計結(jié)果,表3列舉了不同置信概率下2種方法的區(qū)間估計可靠度。由表可知:GBM(1,1)模型的區(qū)間估計可靠度高達(dá)90%以上,優(yōu)于Bootstrap方法。
表3 區(qū)間估計可靠度比較Tab.3 Comparison rusult of the interval estimation reliability
3.3.3 動態(tài)不確定度
圖6給出了不同置信概率下Bootstrap方法和GBM(1,1)模型對爆室溫度動態(tài)測量結(jié)果的估計不確定度(m=7,B=1 000,Q=10)。
圖6 不同置信概率下Bootstrap方法和GBM(1,1)模型的動態(tài)不確定度Fig.6 Dynamic uncertainties of the Bootstrap method and the GBM (1,1) model with different confidence levels
不同置信概率下Bootstrap方法和GBM(1,1)模型2種方法的平均不確定度如表4所示。由表可知,GBM(1,1)模型的平均不確定度略大于Bootstrap方法,且置信水平越高,平均不確定度越大。
表4 平均不確定度比較Tab.4 The mean uncertainty comparison between Bootstrap method and GBM (1,1) model ℃
本文運用灰自助模型GBM(1,1)對密閉空間爆炸實驗中爆室氣體溫度動態(tài)測量數(shù)據(jù)序列的不確定性進(jìn)行了研究,并選取一個基于統(tǒng)計特性的評估參數(shù)(平均不確定度)和3個動態(tài)評估參數(shù)(估計真值、估計區(qū)間、動態(tài)不確定度),對GBM(1,1)模型的評估結(jié)果進(jìn)行了分析。
相比于灰色模型GM(1,1)和Bootstrap方法,GBM(1,1)模型包含了冪函數(shù)和線性回歸的特性,改進(jìn)了傳統(tǒng)的差分模型和離散模型的缺陷,具有較強(qiáng)的趨勢性,且累加生成有效減弱了原始測量數(shù)據(jù)的隨機(jī)性誤差。因此,GBM(1,1)模型表現(xiàn)出更加完善的預(yù)測機(jī)制,能夠?qū)崿F(xiàn)對被測量變化趨勢的實時跟蹤。
估計區(qū)間[XL,XU]描述了被測量的動態(tài)變化范圍,隨著置信概率P的增大,區(qū)間寬度逐漸增大,對被測量波動范圍的包絡(luò)越完整。平均不確定度Umean是被測量動態(tài)不確定度的均值,可以作為被測量隨機(jī)性的評價指標(biāo)。平均不確定度Umean最理想的評價方法為在區(qū)間估計可靠度PB=100%的條件下,Umean取最小值;但在實際問題中,難以實現(xiàn)PB=100%,因此需要選擇合適的m和B,使PB盡可能趨近于100%。由表3可知,在m=7,B=1 000,Q=10的前提下,Bootstrap方法和GBM(1,1)模型的區(qū)間估計可靠度PB均無法達(dá)到100%,且Bootstrap方法的區(qū)間估計可靠度遠(yuǎn)低于GBM(1,1)模型,盡管其平均不確定度Umean小于GBM(1,1)模型,但綜合考慮置信概率P和估計可靠度PB等參數(shù),認(rèn)為GBM(1,1)模型的動態(tài)不確定度評估效果優(yōu)于Bootstrap方法。
爆炸過程的測量屬于動態(tài)測量[14],本文將灰自助方法應(yīng)用于密閉空間內(nèi)爆炸氣體溫度的動態(tài)不確定度評估,構(gòu)建了灰自助模型GBM(1,1),并通過估計真值、估計區(qū)間、估計可靠度、平均不確定度等參數(shù)驗證了該模型對動態(tài)不確定度評定的有效性。GBM(1,1)模型融合了灰色模型GM(1,1)和Bootstrap方法的優(yōu)勢,可以在無任何先驗信息的前提下準(zhǔn)確模擬被測量的概率分布,并實時預(yù)測被測量瞬時值的變化趨勢。相比于Bootstrap方法和灰色模型GM(1,1),灰自助模型GBM(1,1)具有更高的真值估計精度和區(qū)間估計可靠度,估計區(qū)間能夠較為完整地包絡(luò)被測量的動態(tài)變化范圍。在密閉空間爆炸實驗中,該方法能夠?qū)Ρ覝囟鹊膭討B(tài)測量不確定度做出準(zhǔn)確評定,具有較高的實用價值。