孫銘澤，王艷平，羅太明，饒?jiān)迄i，張騰月，曾 丹

(1.中國(guó)兵器工業(yè)火炸藥工程與安全技術(shù)研究院， 北京 100053；2.中國(guó)兵器科學(xué)研究院， 北京 100089)

引 言

靜電放電一直是火炸藥安全生產(chǎn)領(lǐng)域重點(diǎn)關(guān)注的安全風(fēng)險(xiǎn)因素之一，其放電火花可引起火炸藥燃燒或爆炸，繼而對(duì)火炸藥生產(chǎn)、儲(chǔ)運(yùn)等過(guò)程造成潛在風(fēng)險(xiǎn)和重大安全隱患[1-5]。因此，火炸藥的靜電感度，或者靜電放電引燃的最小能量(Emin)，是表征火炸藥安全性能的重要參數(shù)之一。

由于靜電放電最小引燃能量(Emin)難以直接且準(zhǔn)確地測(cè)量，加之環(huán)境因素影響具有隨機(jī)性和偶然性，故當(dāng)前火炸藥的靜電安全特性主要通過(guò)靜電感度來(lái)表征，即火炸藥在靜電放電火花的作用下，發(fā)生燃燒或爆炸的難易程度[2,4]。

當(dāng)前，靜電感度試驗(yàn)主要基于數(shù)理統(tǒng)計(jì)方法設(shè)定試驗(yàn)程序，開(kāi)展系列靜電發(fā)火試驗(yàn)，計(jì)算得到靜電放電能量與火炸藥響應(yīng)概率的關(guān)系[7-12]?，F(xiàn)行火炸藥?kù)o電感度試驗(yàn)標(biāo)準(zhǔn)[13-17]都采用升降法設(shè)定試驗(yàn)程序，準(zhǔn)確度較差，只能得到被試樣品50%響應(yīng)概率所對(duì)應(yīng)的靜電放電能量(E50)，不能準(zhǔn)確測(cè)算極小概率(如百萬(wàn)分之一)響應(yīng)點(diǎn)，一定程度上影響了火炸藥?kù)o電安全風(fēng)險(xiǎn)評(píng)估的準(zhǔn)確性和有效性。

除升降法外，典型的感度試驗(yàn)數(shù)理統(tǒng)計(jì)方法還包括蘭利法和最優(yōu)化D法。蘭利法是一種變步長(zhǎng)的升降法[9]，被應(yīng)用于部分炸藥短脈沖起爆感度試驗(yàn)方法和標(biāo)準(zhǔn)中，可快速得到50%概率感度數(shù)據(jù)；最優(yōu)化D法屬于最優(yōu)化方法[7-8]，操作和計(jì)算流程較為復(fù)雜，但試驗(yàn)準(zhǔn)確性較高，已被應(yīng)用于美國(guó)感度試驗(yàn)軍用標(biāo)準(zhǔn)。但因缺少相關(guān)研究驗(yàn)證，蘭利法和最優(yōu)化D法均未被應(yīng)用于國(guó)內(nèi)的火炸藥?kù)o電感度試驗(yàn)方法和標(biāo)準(zhǔn)之中。

因此，亟需開(kāi)展基于不同數(shù)理統(tǒng)計(jì)方法的靜電感度試驗(yàn)準(zhǔn)確性研究，優(yōu)選或改進(jìn)現(xiàn)有數(shù)理統(tǒng)計(jì)方法，在最少試驗(yàn)次數(shù)內(nèi)，獲得更準(zhǔn)確地火炸藥?kù)o電感度分布，以提高對(duì)百萬(wàn)分之一概率發(fā)火能量(E1ppm)等極小概率響應(yīng)點(diǎn)的測(cè)算準(zhǔn)確性，也能更準(zhǔn)確地估計(jì)Emin等火炸藥安全性能參數(shù)。

本研究利用蒙特卡羅方法進(jìn)行計(jì)算機(jī)模擬試驗(yàn)對(duì)比，分析基于典型數(shù)理統(tǒng)計(jì)方法進(jìn)行靜電感度試驗(yàn)的特點(diǎn)，并對(duì)現(xiàn)有的蘭利法進(jìn)行優(yōu)化和改進(jìn)，開(kāi)展升降法、蘭利法、最優(yōu)化D法以及優(yōu)化和改進(jìn)后的蘭利法4種方法的準(zhǔn)確性研究，以期為優(yōu)化火炸藥?kù)o電感度試驗(yàn)方法和標(biāo)準(zhǔn)提供理論依據(jù)，也為獲得更準(zhǔn)確的火炸藥?kù)o電最低引燃能量提供方法和途徑。

1 典型數(shù)理統(tǒng)計(jì)方法及對(duì)比分析

國(guó)內(nèi)現(xiàn)行的感度試驗(yàn)數(shù)理統(tǒng)計(jì)方法標(biāo)準(zhǔn)[18]中推薦了升降法和蘭利法等方法，美國(guó)軍用標(biāo)準(zhǔn) MIL-DTL-23659D于2013年將最優(yōu)化D法作為火炸藥感度試驗(yàn)數(shù)理統(tǒng)計(jì)方法。因此，升降法、蘭利法和最優(yōu)化D法都是目前國(guó)內(nèi)外最具代表性的感度試驗(yàn)方法，本研究將對(duì)這3種典型數(shù)理統(tǒng)計(jì)方法展開(kāi)研究分析。

1.1 升降法

升降法最早由 W. J. Dixon 和 A. M. Mood 在1948年提出[5,10,12,19]，其流程簡(jiǎn)單，是國(guó)內(nèi)最常用的火炸藥?kù)o電感度試驗(yàn)數(shù)理統(tǒng)計(jì)方法。GJB/Z 377A-94對(duì)升降法的數(shù)據(jù)處理給出的均值計(jì)算公式為：

(1)

(2)

(3)

(4)

1.2 蘭利法

蘭利法由H. J. Langlie于1963年提出，并由Guillaume 在2010年補(bǔ)充了操作細(xì)節(jié)[20]，是一種變步長(zhǎng)升降法，可快速得到50%概率感度數(shù)據(jù)，具有可操作性高的特點(diǎn)，試驗(yàn)前只需提供全不響應(yīng)刺激量xL作為試驗(yàn)下限、全響應(yīng)刺激量xU作為試驗(yàn)上限，每次試驗(yàn)點(diǎn)的計(jì)算只涉及簡(jiǎn)單的四則運(yùn)算。首次試驗(yàn)時(shí)選擇上下限平均值作為初始靜電刺激量x0，隨著試驗(yàn)次數(shù)增加，測(cè)點(diǎn)集中在感度分布均值處，可快速收斂獲得靜電感度分布均值，但不能準(zhǔn)確測(cè)算分布方差，繼而影響了如E1ppm等依賴(lài)于方差的極小概率靜電引燃能量的測(cè)算準(zhǔn)確性。

1.3 最優(yōu)化D法

最優(yōu)化D法是由Barry T. Neyer在1994年提出的一種基于Fisher信息矩陣的靜電感度試驗(yàn)數(shù)理統(tǒng)計(jì)方法[7,8,21,22]，每次試驗(yàn)前需先計(jì)算靜電感度分布，然后優(yōu)化新的刺激量，使試驗(yàn)后靜電感度分布的整體測(cè)算誤差最小，但每次試驗(yàn)需要求解復(fù)雜的非線(xiàn)性?xún)?yōu)化問(wèn)題，程序設(shè)定和計(jì)算十分復(fù)雜，操作難度大。此外，最優(yōu)化D法需在預(yù)試驗(yàn)后進(jìn)行，Barry T. Neyer 最初提出的方法被稱(chēng)為Neyer-D法，基于折半查找法進(jìn)行預(yù)試驗(yàn)[22]。華東師范大學(xué)的吳蔚[23]對(duì)其進(jìn)行改進(jìn)提出蘭利-最優(yōu)化D法，可同時(shí)獲得較為準(zhǔn)確的感度分布均值和方差，并對(duì)E1ppm等極小概率靜電引燃能量進(jìn)行準(zhǔn)確度較高的測(cè)算。因此，本研究選擇蘭利-最優(yōu)化D法作為最優(yōu)化D法的代表。

基于最優(yōu)化D法進(jìn)行靜電感度試驗(yàn)時(shí)，若試驗(yàn)次數(shù)較少時(shí)，則會(huì)因樣本不足以覆蓋極端情況而使測(cè)算的方差整體偏小，可通過(guò)引入大于1的偏量修正因子提高方差測(cè)算的準(zhǔn)確性：

σC=β×σ

(5)

式中：β為修正因子，是試驗(yàn)次數(shù)N的函數(shù)；σ為修正前方差；σc為修正后方差。本研究基于最優(yōu)化D法開(kāi)展不同次數(shù)的模擬試驗(yàn)，計(jì)算得到靜電感度分布方差與設(shè)定真值之比作為偏量修正因子βi：

(6)

式中：Ni為模擬試驗(yàn)次數(shù)；βi為Ni次模擬試驗(yàn)對(duì)應(yīng)的偏量修正因子；σi為Ni次模擬試驗(yàn)后靜電感度分布方差測(cè)算值；σ0為靜電感度分布方差設(shè)定真值。將試驗(yàn)次數(shù)Ni與對(duì)應(yīng)的偏量修正因子與組合為數(shù)據(jù)對(duì)(Ni,βi)，并基于最小二乘法擬合得到最優(yōu)化D法的偏量修正因子β與試驗(yàn)次數(shù)N之間的數(shù)學(xué)表達(dá)式。

本研究分別基于修正前后的最優(yōu)化D法進(jìn)行不同次數(shù)的模擬試驗(yàn)，不僅能驗(yàn)證修正因子的有效性，也能分析基于該方法開(kāi)展靜電感度試驗(yàn)的特性。

2 蘭利法的優(yōu)化和改進(jìn)

針對(duì)蘭利法方差測(cè)算準(zhǔn)確性不高的問(wèn)題，本研究突破了原有方法測(cè)點(diǎn)集中在50%概率響應(yīng)點(diǎn)附近的模式，探索性地提出將測(cè)點(diǎn)分散在30%概率和70%概率響應(yīng)點(diǎn)周?chē)?，以尋找被測(cè)樣品靜電感度分布的離散程度，從而更準(zhǔn)確地捕捉分布方差。

基于改進(jìn)的蘭利法開(kāi)展靜電感度試驗(yàn)的具體做法如下：

(1)試驗(yàn)前設(shè)定全不響應(yīng)刺激量xL作為試驗(yàn)下限、全響應(yīng)刺激量xU作為試驗(yàn)上限，首次試驗(yàn)選擇試驗(yàn)上下限的算數(shù)平均值作為測(cè)點(diǎn)刺激量；

(2)在計(jì)算出的刺激量下進(jìn)行靜電感度試驗(yàn)，記錄第i次試驗(yàn)的刺激量xi和響應(yīng)vi后，從本次試驗(yàn)開(kāi)始往前查找歷史試驗(yàn)記錄。

(4)若查找到j(luò)=1仍未滿(mǎn)足要求，則當(dāng)vi=1時(shí)取定中間變量為xL，反之vi=0則取定中間變量為xU。

(6)重復(fù)步驟(2)～(5)，直至試驗(yàn)次數(shù)滿(mǎn)足設(shè)定要求后，基于最大似然法計(jì)算靜電感度分布均值和方差以及不同響應(yīng)概率臨界刺激量等結(jié)果。

當(dāng)試驗(yàn)次數(shù)足夠多時(shí)，符合條件的xj會(huì)接近樣品50%概率響應(yīng)點(diǎn)，而中間變量取為其和試驗(yàn)上限或下限的均值，則可有效地將測(cè)點(diǎn)分散在較小或較大概率響應(yīng)點(diǎn)周?chē)?，可比傳統(tǒng)蘭利法更準(zhǔn)確地測(cè)算靜電感度方差。同時(shí)，改進(jìn)后的蘭利法雖然操作上略復(fù)雜于原方法，但其每次測(cè)點(diǎn)建議值依然可通過(guò)簡(jiǎn)單的四則運(yùn)算得到，計(jì)算復(fù)雜度遠(yuǎn)低于最優(yōu)化D法。

基于改進(jìn)的蘭利法進(jìn)行靜電感度試驗(yàn)時(shí)，也可根據(jù)(5)式和(6)式引入偏量修正因子減小方差的測(cè)算誤差。本研究通過(guò)模擬試驗(yàn)，擬合得到改進(jìn)蘭利法的偏量修正因子，并分別基于修正前后的改進(jìn)蘭利法進(jìn)行不同次數(shù)的模擬試驗(yàn)，不僅能驗(yàn)證修正因子的有效性，也能分析基于該方法開(kāi)展靜電感度試驗(yàn)的特性。

3 蒙特卡羅試驗(yàn)研究方法與流程

蒙特卡羅方法是一種對(duì)隨機(jī)過(guò)程進(jìn)行模擬后統(tǒng)計(jì)結(jié)果的計(jì)算機(jī)模擬試驗(yàn)方法。本研究將通過(guò)蒙特卡羅法開(kāi)展計(jì)算機(jī)數(shù)值模擬試驗(yàn)，分析基于升降法、蘭利法、最優(yōu)化D法、改進(jìn)的蘭利法4種不同數(shù)理統(tǒng)計(jì)方法獲得結(jié)果的準(zhǔn)確性，具體研究方法如下：

(1)設(shè)定初始刺激量x0=μ，步長(zhǎng)d=σ，分別基于勘誤前后的升降法進(jìn)行模擬試驗(yàn)，對(duì)比勘誤前后結(jié)果的準(zhǔn)確性；分別設(shè)定d=σ和d=σ/2進(jìn)行模擬試驗(yàn)，分析不同步長(zhǎng)對(duì)升降法準(zhǔn)確性的影響。

(2)設(shè)定試驗(yàn)上下限為[xU,xL]=[μ+3σ,μ-3σ]，基于蘭利法進(jìn)行不同次數(shù)的模擬試驗(yàn)。

(3)基于最優(yōu)化D法進(jìn)行模擬試驗(yàn)，擬合偏量修正因子；設(shè)定試驗(yàn)上下限為[μ+3σ,μ-3σ]，對(duì)比有無(wú)偏量修正因子時(shí)，基于最優(yōu)化D法進(jìn)行模擬試驗(yàn)的準(zhǔn)確性；分別設(shè)定試驗(yàn)上下限為[μ-3σ,μ+3σ]和[μ-2σ,μ+σ]兩種情況，基于最優(yōu)化D法進(jìn)行模擬試驗(yàn)，對(duì)比分析試驗(yàn)上下限對(duì)其準(zhǔn)確性的影響。

(4)基于優(yōu)化和改進(jìn)后的蘭利法進(jìn)行模擬試驗(yàn)，擬合偏量修正因子；設(shè)定試驗(yàn)上下限為[μ+3σ,μ-3σ]，對(duì)比有無(wú)偏量修正因子時(shí)，基于改進(jìn)的蘭利法進(jìn)行模擬試驗(yàn)的準(zhǔn)確性；分別設(shè)定試驗(yàn)上下限為[μ-3σ,μ+3σ] 和[μ-2σ,μ+σ]兩種情況，基于改進(jìn)的蘭利法進(jìn)行模擬試驗(yàn)，對(duì)比分析試驗(yàn)上下限對(duì)其準(zhǔn)確性的影響。

(5)對(duì)比分析基于上述4種方法開(kāi)展模擬試驗(yàn)的測(cè)點(diǎn)選擇的有效性，靜電感度分布均值和方差的測(cè)算準(zhǔn)確性，以及小概率響應(yīng)點(diǎn)的測(cè)算準(zhǔn)確性。

上述研究?jī)?nèi)容中開(kāi)展模擬試驗(yàn)的具體流程如下：

(1)設(shè)定靜電感度分布的真值，將其設(shè)為均值μ0=10，方差σ0=1的正態(tài)分布；

(2)按方法計(jì)算每次試驗(yàn)的測(cè)點(diǎn)刺激量xi，并基于靜電感度分布真值計(jì)算響應(yīng)概率，并生成試驗(yàn)結(jié)果vi；

(3)設(shè)定模擬試驗(yàn)的次數(shù)N，完成一組N次試驗(yàn)后，根據(jù)測(cè)點(diǎn)刺激量x1,x2,…,xN和響應(yīng)情況v1,v2,…,vN，基于所選方法計(jì)算靜電感度分布均值和方差的測(cè)算值μα、σα，其中α=1,2,...標(biāo)記不同組模擬試驗(yàn)；

(4)基于測(cè)算出的均值和方差計(jì)算1%、0.01%、10-4%等小概率響應(yīng)點(diǎn)臨界刺激量：

xα(p)=μα+upσα

(7)

式中：μα、σα為第α組模擬試驗(yàn)得到的靜電感度分布均值和方差；xα(p)為概率p響應(yīng)點(diǎn)靜電臨界刺激量；up為概率p對(duì)應(yīng)的標(biāo)準(zhǔn)配位數(shù)，例如u1%≈-2.33，u1ppm≈-4.75；

(5)重復(fù)步驟(2)～(4)進(jìn)行M組模擬試驗(yàn)后，統(tǒng)計(jì)分析靜電感度分布測(cè)算值與真值之間誤差的均值與標(biāo)準(zhǔn)差：

ΔXα=Xα-X0

(8)

(9)

(10)

式中：Xα為第α組模擬試驗(yàn)測(cè)算結(jié)果(可為感度均值、方差或小概率響應(yīng)點(diǎn))；X0為對(duì)應(yīng)設(shè)定真值；ΔXα為第α組模擬試驗(yàn)測(cè)算誤差；ME(ΔXα)為第α組模擬試驗(yàn)測(cè)算誤差均值；SD(ΔXα)為第α組模擬試驗(yàn)測(cè)算誤差標(biāo)準(zhǔn)差；M為模擬試驗(yàn)組數(shù)；

(6)計(jì)算模擬試驗(yàn)測(cè)算結(jié)果誤差在約68.26%置信度(單σ)下的置信區(qū)間，并繪制誤差線(xiàn)圖：

ΔXU=ME(ΔX)+Z(1-p)/2SD(ΔX)

(11)

ΔXD=ME(ΔX)+Z(1-p)/2SD(ΔX)

(12)

式中：p為置信度；ΔXU為置信區(qū)間上限；XD為置信區(qū)間下限；ME(ΔX)為結(jié)果X誤差的均值；SD(ΔX)為結(jié)果X誤差的標(biāo)準(zhǔn)差;；Z為正態(tài)分布標(biāo)準(zhǔn)分?jǐn)?shù)，為計(jì)算方便本研究將其取1對(duì)應(yīng)約68.26%的置信度。

當(dāng)靜電感度分布均值和方差測(cè)算誤差的均值和標(biāo)準(zhǔn)差相對(duì)較小時(shí)，小概率響應(yīng)點(diǎn)臨界刺激量測(cè)算誤差的均值和標(biāo)準(zhǔn)差可由其計(jì)算得到：

ME[Δx(p)]=ME(Δμ)+upME(Δσ)

(13)

SD[Δx(p)]=SD(Δμ)+upSD(Δσ)

(14)

式中：μ、σ、x(p)和up的定義同式(7)；ME為其均值；SD為其標(biāo)準(zhǔn)差。因此，本研究在通過(guò)數(shù)值模擬試驗(yàn)分析4種方法特點(diǎn)時(shí)只分析了靜電感度均值和方差。

試驗(yàn)中需令M?N(將M取為100×N)，使單組模擬試驗(yàn)不確定性帶來(lái)的誤差降低到可忽略，保證統(tǒng)計(jì)的有效性。此外，取不同的試驗(yàn)次數(shù)N進(jìn)行模擬試驗(yàn)，研究試驗(yàn)次數(shù)對(duì)結(jié)果準(zhǔn)確性的影響，并將最大試驗(yàn)次數(shù)取為1000。蒙特卡羅計(jì)算機(jī)模擬試驗(yàn)的具體技術(shù)途徑如圖1所示。

圖1 蒙特卡羅計(jì)算機(jī)模擬試驗(yàn)流程圖Fig.1 Flow chart of Monte Carlo computer simulation test

4 結(jié)果對(duì)比與準(zhǔn)確性分析

針對(duì)升降法、蘭利法、最優(yōu)化D法和改進(jìn)的蘭利法進(jìn)行蒙特卡羅模擬試驗(yàn)的結(jié)果進(jìn)行對(duì)比分析，具體見(jiàn)圖2。

4.1 升降法

模擬試驗(yàn)結(jié)果顯示，按照前文1.1中介紹的勘誤模式1的升降法進(jìn)行模擬試驗(yàn)，當(dāng)試驗(yàn)次數(shù)達(dá)1000次時(shí)，絕對(duì)誤差的統(tǒng)計(jì)平均值仍有約+0.32；而經(jīng)前文1.1中勘誤模式2的升降法對(duì)靜電感度分布均值和方差的測(cè)算都較準(zhǔn)確，誤差統(tǒng)計(jì)平均值幾乎為0，且測(cè)算誤差隨試驗(yàn)次數(shù)增加而降低，如圖2(a)所示，圖2和圖3中的所有子圖均以誤差線(xiàn)的形式展現(xiàn)測(cè)算誤差，因此，本研究采用模式2對(duì)GJB/Z 377A-94中介紹的升降法勘誤；此外，當(dāng)步長(zhǎng)設(shè)定為靜電感度分布方差真值時(shí)，升降法的準(zhǔn)確性尚可，但當(dāng)步長(zhǎng)設(shè)定偏離方差真值時(shí)，其在試驗(yàn)次數(shù)較少時(shí)，對(duì)方差的測(cè)算誤差顯著增大，如圖2(b)所示。

因此，基于升降法進(jìn)行試驗(yàn)可對(duì)靜電感度分布的均值有效測(cè)算，得到較準(zhǔn)確的50%響應(yīng)點(diǎn)，但對(duì)方差測(cè)算誤差較大，所以無(wú)法獲得準(zhǔn)確的極小概率響應(yīng)點(diǎn)；而且升降法測(cè)試結(jié)果的準(zhǔn)確性依賴(lài)于初始步長(zhǎng)設(shè)定，穩(wěn)定性較低。

4.2 蘭利法

模擬試驗(yàn)結(jié)果顯示，基于蘭利法進(jìn)行靜電感度試驗(yàn)時(shí)，在約40次試驗(yàn)內(nèi)可以快速獲得較準(zhǔn)確的靜電感度分布均值；但隨著試驗(yàn)次數(shù)增加，試驗(yàn)中的靜電刺激量會(huì)聚集在測(cè)算出的50%概率靜電響應(yīng)點(diǎn)處，無(wú)法獲取其他響應(yīng)點(diǎn)信息，導(dǎo)致測(cè)算誤差隨試驗(yàn)次數(shù)增加反而增大，如圖2(c)和(d)所示。此外，蘭利法和升降法類(lèi)似，只對(duì)靜電感度分布均值測(cè)算較準(zhǔn)，而對(duì)方差的測(cè)算誤差較大。

4.3 最優(yōu)化D法

根據(jù)模擬試驗(yàn)結(jié)果，可擬合最優(yōu)化D法偏量修正因子表達(dá)式：

(15)

式中：β為修正因子；N為試驗(yàn)次數(shù)；a為擬合參數(shù)。模擬試驗(yàn)結(jié)果顯示，基于最優(yōu)化D法進(jìn)行靜電感度試驗(yàn)可同時(shí)準(zhǔn)確獲得其分布均值和方差，且偏量修正因子會(huì)顯著降低100次試驗(yàn)以?xún)?nèi)的方差的測(cè)算誤差，如圖2(e)所示，只有100次試驗(yàn)內(nèi)修正前后的感度分布方差測(cè)算誤差有所區(qū)別，其他上下限條件下的結(jié)果幾乎在圖像上重合，如圖2(f)所示，說(shuō)明上下限取值對(duì)最優(yōu)化D法獲得的靜電感度分布測(cè)算誤差無(wú)明顯影響。

圖2 四種數(shù)理統(tǒng)計(jì)方法在不同情況下試驗(yàn)結(jié)果測(cè)算誤差對(duì)比Fig.2 Comparison of the test result error under different conditions by four different mathematical statistical methods

因此，引入偏量修正因子后，最優(yōu)化D法對(duì)靜電感度分布的均值和方差可同時(shí)準(zhǔn)確測(cè)算，且結(jié)果準(zhǔn)確性不依賴(lài)于試驗(yàn)上下限的選取，穩(wěn)定性較高。

4.4 改進(jìn)的蘭利法

根據(jù)模擬試驗(yàn)結(jié)果，可擬合改進(jìn)的蘭利法偏量修正因子表達(dá)式：

(16)

式中：β為修正因子；N為試驗(yàn)次數(shù)；a,b為擬合參數(shù)。模擬試驗(yàn)結(jié)果顯示，基于改進(jìn)的蘭利法進(jìn)行靜電感度試驗(yàn)，也可同時(shí)準(zhǔn)確地測(cè)算分布均值和方差，修正因子可有效改善100次試驗(yàn)內(nèi)方差測(cè)算的整體偏差，如圖2(g)所示，只有100次試驗(yàn)內(nèi)修正前后的感度分布方差測(cè)算誤差有所區(qū)別，其他情況下修正前后的數(shù)據(jù)線(xiàn)在圖中重合。此外，上下限選擇會(huì)影響改進(jìn)蘭利法的測(cè)算結(jié)果，特別是方差測(cè)算結(jié)果準(zhǔn)確性，但其穩(wěn)定性依然優(yōu)于原始的蘭利法，如圖2(h)所示。

因此，引入偏量修正因子后，改進(jìn)的蘭利法可憑借遠(yuǎn)低于最優(yōu)化D法的操作難度和計(jì)算復(fù)雜度，對(duì)靜電感度分布均值和方差同時(shí)進(jìn)行準(zhǔn)確測(cè)算，只是穩(wěn)定性略差于最優(yōu)化D法。

4.5 方法準(zhǔn)確性對(duì)比與分析

首先，對(duì)比4種不同方法試驗(yàn)的測(cè)點(diǎn)分布：升降法測(cè)點(diǎn)集中在靜電感度分布均值及其兩側(cè)；蘭利法測(cè)點(diǎn)聚集在均值附近；最優(yōu)化D法的測(cè)點(diǎn)分散在約15%和約85%概率響應(yīng)點(diǎn)兩處；而改進(jìn)的蘭利法測(cè)點(diǎn)均勻地分散在約30%概率和約70%概率響應(yīng)點(diǎn)周?chē)?，如圖3(a)所示。

圖3 不同數(shù)理統(tǒng)計(jì)方法測(cè)點(diǎn)對(duì)比分析Fig.3 Comparison and analysis of measuring points with different mathematical statistics methods

利用Fisher信息矩陣的行列式，即D函數(shù)，可評(píng)估試驗(yàn)中測(cè)點(diǎn)選擇的有效性。在測(cè)點(diǎn)位置進(jìn)行一次新的試驗(yàn)后，靜電感度分布測(cè)算誤差矩陣的行列式等于D函數(shù)的倒數(shù)，可反映測(cè)點(diǎn)對(duì)靜電感度分布準(zhǔn)確測(cè)算的有效貢獻(xiàn)。結(jié)果顯示，從蘭利法、升降法、改進(jìn)的蘭利法到最優(yōu)化D法，D函數(shù)遞增，測(cè)點(diǎn)選擇有效性依次增加，如圖3(b)所示。

其次，對(duì)比不同方法對(duì)靜電感度分布均值與方差的測(cè)算誤差。結(jié)果顯示，與均值為10和方差為1的設(shè)定真值比較，升降法、最優(yōu)化D法和改進(jìn)的蘭利法對(duì)靜電感度分布的測(cè)算誤差均隨試驗(yàn)次數(shù)增加而降低；升降法對(duì)均值的測(cè)算誤差尚可，但對(duì)方差的測(cè)算誤差最大；蘭利法能快速測(cè)算出較準(zhǔn)確的均值，但試驗(yàn)次數(shù)增多后測(cè)算誤差反而增大；優(yōu)化和改進(jìn)的蘭利法對(duì)均值與方差的測(cè)算誤差均尚可；最優(yōu)化D法表現(xiàn)最佳，對(duì)均值和方差的測(cè)算誤差均最小，如圖4(a)和(b)所示。

圖4 不同數(shù)理統(tǒng)計(jì)方法均值和方差誤差對(duì)比Fig.4 Comparison of measurement errors of mean and variance by different mathematical statistics methods

僅從均值和方差的測(cè)算上考慮，優(yōu)化和改進(jìn)后的蘭利法和最優(yōu)化D法的準(zhǔn)確性難分伯仲。因此，進(jìn)一步分析其對(duì)小概率響應(yīng)點(diǎn)的測(cè)算準(zhǔn)確性，與百萬(wàn)分之一、萬(wàn)分之一和百分之一3個(gè)小概率響應(yīng)點(diǎn)的設(shè)定真值(約5.245、6.28、7.67)比較，各方法的測(cè)算誤差如圖5所示，典型試驗(yàn)次數(shù)下百萬(wàn)分之一概率響應(yīng)點(diǎn)的測(cè)算結(jié)果的測(cè)算誤差對(duì)比如表1所示。

圖5 不同數(shù)理統(tǒng)計(jì)方法小響應(yīng)概率點(diǎn)結(jié)果比較Fig.5 Comparison of results of small response probability points of different mathematical statistics methods

表1 典型試驗(yàn)次數(shù)下10-4%概率響應(yīng)點(diǎn)測(cè)算誤差68.26%置信區(qū)間對(duì)比Table 1 Comparison of 68.26% confidence interval of 10-4% probability response points error under typical test times

結(jié)果顯示，最優(yōu)化D法對(duì)小概率響應(yīng)點(diǎn)的測(cè)算準(zhǔn)確性?xún)?yōu)于改進(jìn)的蘭利法，其在較少的250次試驗(yàn)內(nèi)，可將E1ppm的測(cè)算誤差的68.26%置信區(qū)間控制在0.456以?xún)?nèi)，使其相對(duì)誤差小于10%。綜上所述，最優(yōu)化D法是四種數(shù)理統(tǒng)計(jì)方法中最適合靜電感度試驗(yàn)的方法；優(yōu)化和改進(jìn)的蘭利法具有操作便捷的優(yōu)點(diǎn)，試驗(yàn)中的計(jì)算復(fù)雜度遠(yuǎn)低于最優(yōu)化D法，其將測(cè)點(diǎn)集中在較小和較大概率響應(yīng)點(diǎn)的改進(jìn)思路，可為進(jìn)一步優(yōu)化靜電感度試驗(yàn)數(shù)理統(tǒng)計(jì)方法提供參考。

5 結(jié) 論

(1)通過(guò)蒙特卡羅法計(jì)算機(jī)模擬試驗(yàn)，分析了基于現(xiàn)有典型數(shù)理統(tǒng)計(jì)方法進(jìn)行靜電感度試驗(yàn)的特點(diǎn)：升降法操作最簡(jiǎn)單，可對(duì)靜電感度分布均值和方差進(jìn)行較準(zhǔn)確測(cè)算，但對(duì)初值和步長(zhǎng)選擇敏感，穩(wěn)定性差；蘭利法也具有操作簡(jiǎn)單的特點(diǎn)，可快速獲得分布均值，但對(duì)方差的測(cè)算誤差較大；最優(yōu)化D法操作難度和計(jì)算復(fù)雜度較高，但對(duì)感度分布均值和方差測(cè)算準(zhǔn)確性和穩(wěn)定性都最高；增加試驗(yàn)次數(shù)可提高升降法和最優(yōu)化D法的測(cè)算準(zhǔn)確性。

(2)對(duì)比研究發(fā)現(xiàn)：最優(yōu)化D法可在較少的試驗(yàn)次數(shù)內(nèi)最準(zhǔn)確且穩(wěn)定地測(cè)算靜電感度分布，在250次試驗(yàn)內(nèi)可將E1ppm的測(cè)算相對(duì)誤差的68.26%置信區(qū)間控制在10%以?xún)?nèi)，能準(zhǔn)確表征火炸藥?kù)o電最低引燃能量，為火炸藥?kù)o電安全風(fēng)險(xiǎn)的定量化分析提供基礎(chǔ)數(shù)據(jù)。但其使用過(guò)程需要進(jìn)行大量復(fù)雜數(shù)學(xué)計(jì)算，試驗(yàn)時(shí)對(duì)人員技術(shù)能量要求較高，不適合作為常規(guī)靜電感度試驗(yàn)的方法，未來(lái)可研發(fā)專(zhuān)用試驗(yàn)程序軟件，自動(dòng)完成復(fù)雜的測(cè)點(diǎn)刺激量計(jì)算和試驗(yàn)結(jié)果計(jì)算等過(guò)程，輔助操作人員基于最優(yōu)化D法開(kāi)展靜電感度試驗(yàn)。

(3)優(yōu)化和改進(jìn)后的蘭利法，突破了傳統(tǒng)靜電感度試驗(yàn)方法測(cè)點(diǎn)集中在50%概率響應(yīng)點(diǎn)的思路模式，將測(cè)點(diǎn)分散在約30%和約70%概率響應(yīng)點(diǎn)周?chē)?，在程序設(shè)定和計(jì)算復(fù)雜度遠(yuǎn)低于最優(yōu)化D法的同時(shí)，達(dá)到與其相似的準(zhǔn)確性。

(4)提出的分散測(cè)點(diǎn)的思路可有效測(cè)算方差的準(zhǔn)確性，可作為未來(lái)進(jìn)一步研究靜電感度試驗(yàn)方法和程序的基礎(chǔ)，以期建立一種新的感度試驗(yàn)數(shù)理方法，在提高極小概率響應(yīng)點(diǎn)測(cè)算準(zhǔn)確性的同時(shí)，具有較為簡(jiǎn)單、可操作性強(qiáng)的特點(diǎn)，方便用于常規(guī)靜電感度試驗(yàn)。