基于高斯多峰法的密閉空間爆炸特性曲線擬合

謝 威， 蔣新生， 徐建楠， 周 毅， 王 冬

(陸軍勤務學院 油料系，重慶 401311)

在石油、化工安全防護中，研究易燃、易爆介質產生的爆炸超壓等特性有著極其重要的意義。通過對其發展規律的研究探索，能有效地預測爆炸產生的危害，以便對爆炸防護裝置的強度等參數進行設計、計算、分析。對于爆炸超壓的預測和計算，國內外許多學者都已有較多研究，主要包括空氣中爆炸超壓[1-4]到密閉容器、有約束和無約束的泄壓容器爆炸超壓的預測計算研究[5]。由于受限空間爆炸的傳播受距離、擾動、障礙物等的影響會加速，甚至可能形成爆轟，造成更大的破壞，引起了研究者們的重視[6]。而且，不同尺度的受限空間，其爆炸規律也不盡相同。所使用的研究手段主要有三種：①基于試驗數據的信號處理[7-8]；②基于爆炸理論的理論推導；③基于數值計算與爆炸理論相結合的數值模擬。

針對密閉空間或容器爆炸的研究，王淑蘭等[9]對密閉容器內中心弱點火條件下的烴類氣體的燃爆特性進行了計算研究；趙衡陽總結了密閉容器爆炸的三種簡化模型：等溫爆炸模型、絕熱爆炸模型和一般模型，并將密閉容器氣體或粉塵爆炸的理論和試驗結果進行比較分析，驗證了模型的可靠性；嚴清華等[10]利用數值模擬的手段，考慮流體動力學和化學反應動力學，對大型球型密閉容器內的可燃氣體爆炸過程進行了數值模擬，彌補了解析解模型如等溫爆炸模型不能對濃度、密度、溫度等流場進行定量描述的缺點。但是，由于密閉容器爆炸的理論研究、數值模擬研究往往與實際環境、實際需求相差較大，得到的結論再投入工程應用的部分比較少，因此，利用試驗數據結合理論分析來反推和修正理論數值進而投入應用，是有效解決問題的途徑之一。

本文提出用高斯多峰擬合法對爆炸特性曲線進行擬合，利用多個高斯單峰函數的線性疊加趨近真實曲線，提高擬合精度，通過比較不同數量分峰時的殘差分布，選擇一個能較好反映原曲線總體特征的簡化表達式，對進一步量化分析試驗特性參量，如初始濃度、溫度、壓力等與爆炸的聯系具有重要意義。

1 高斯多峰法應用于爆炸特性曲線擬合

高斯多峰法應用到信號波的處理上的研究比較廣泛，比如Giorgos等[11]將其應用到波形激光雷達的信號處理上,Samjin等[12]將其應用于心臟頻譜的分割，Gabriel等[13]將其應用于復雜色譜分析等等，而國內目前比較常見的應用領域是統計學與光譜分析[14-16]，在其他領域的應用比較少見，但其本身卻是很好的求曲線擬合的一種方法。而且高斯多峰的疊加性質可以很好契合爆炸過程中各種爆炸波(泛指爆炸過程中所有具有可疊加性質、震蕩特性的爆炸表征)，比如，壓力波、速度波、火焰波、聲波等等的疊加，尤其是在密閉空間中，由于容器等約束的限制，引起爆炸壓力波的徑向傳播受阻，而軸向傳播得到加速，后產生的爆炸壓力波追上之前的爆炸壓力波，不斷疊加，在前方產生更強的疊加壓縮波，當傳播距離足夠長時，最終將形成爆轟波。

在爆炸特征參量-時間(f-t)曲線上，可以將高斯單峰函數看作某一波動的波動能量隨時間的變化過程，多個單峰曲線的疊加，就等于同一時刻不同波動能量的疊加，意味著某一測點某一時刻受多個波動的影響。而在爆炸特征參量-距離(f-x)曲線上，可以將高斯單峰函數看作某時刻某一波動的波動能量隨位置的變化過程，多個單峰曲線的疊加，就等于同一位置不同波動能量的疊加，也意味著某一時刻某一位置受多個波的影響。此外，爆炸過程總伴隨著振動頻率先不斷增大后逐漸降低的壓力波振蕩現象，如需更精確的描述爆炸特性曲線，可以對曲線進行修正

f(t)=M(t)+Di(t)+Fi(t)(i=1,2,3,4,...,n)

(1)

式中：f(t)為修正值；M(t)為擬合值；Di(t)為修正爆炸波高頻區間的振動函數；Fi(t)為修正低頻區間的脈動函數，且Fi(t)的脈動變化是先增強后慢慢變弱，利用頻譜圖分析高頻和低頻段的頻率和振幅，找到對應的函數。

2 高斯多峰擬合法

高斯單峰的函數表達式為

(2)

式中：y0為函數的基準值；W為控制峰寬的峰寬函數；A為高斯函數的面積函數，可正可負；Xc為高斯函數的峰軸。

標準的高斯函數曲線就是正態分布曲線，其基值為0，中心軸為x軸，積分面積恒為1，如圖1所示。

圖1 標準的高斯函數曲線Fig.1 Standard Gaussian function curve

高斯多峰擬合法就是利用高斯單峰函數在不同位置的疊加來擬合出需要的近似曲線，并得到一個函數表達式。高斯多峰的函數表達式為

(3)

當y0=0時，記

(4)

則有

(5)

(6)

也就是說原函數的積分等于分峰后的所有單峰面積求和。

3 密閉容器爆炸試驗

試驗主體為φ700×3 000 mm，耐壓10 MPa的長密閉圓管，點火位置都在圖2左側法蘭中心，壓力測點位于末端中心，可燃氣體為汽油蒸汽，采用真空泵循環，如圖2所示。輔助設備包括：WGDH-5型高能點火器；TST6300壓力數據采集系統；GXH-1050碳氫分析儀——用于測定試驗系統初始狀態時內部的可燃氣體體積分數。為保證預混均勻，在停止充氣后，需繼續循環一段時間然后靜置10 min。若各測量點碳氫分析儀的讀數變動在0.3%范圍內，即說明試驗容器中的可燃氣體混合均勻。試驗共進行5次，以保證試驗的可重復性，并從所得的數據中選取具有代表性而且振動特性比較明顯的一次試驗數據作為試驗結果。

圖2 試驗布置示意圖Fig.2 Experiment layout diagram

測得的密閉空間內末端的油氣爆炸典型壓力-時間曲線，如圖3所示。由于密閉容器爆炸產生的壓力振蕩前期得到不斷加強，當振動頻率足夠大時，會引發容器振動和傳感器干擾，使得試驗采集到的數據存在一定的干擾誤差。所以，在進行數據擬合前，可以先利用傅里葉變換進行適當的濾波預處理，使曲線便于擬合計算，如圖4所示。

圖3 可燃氣體爆炸試驗實測p-t曲線Fig.3 Measured p-t curve of combustible gas explosion

圖4 實測數據濾波后的p-t曲線Fig.4 Data curve of p-t after filtering

由于擬合所采用的數據是基于瞬態壓力數據采集系統采集到的壓力數據，而且傳感器為壓阻式，因而可視其為一維壓力波動數據，方向為垂直于傳感器方向，所以采用高斯多峰法對該壓力-時間曲線的擬合，僅是對該點的壓力波在一維尺度上的近似分析。

4 密閉空間爆炸壓力特性曲線擬合及結果分析

將高斯多峰函數法具體應用到爆炸壓力特性曲線p-t上，對比油氣體積分數為2.0%(其爆炸極限為1.4%～7.6%)時管道爆炸試驗測得的曲線圖與不同峰數擬合得到的曲線圖。具體的操作步驟，如圖5所示。先對試驗曲線進行濾波處理，濾掉高頻的一些干擾和振動，得到一條相對光滑的曲線(見圖4)。然后，利用Origin軟件對曲線進行峰分析和以高斯函數為基函數的多峰曲線擬合[17]。分析濾波后的曲線，取一較小峰數(n=6)對曲線擬合，得到的圖形如圖6所示。但是，由于Origin軟件目前的最大迭代步數限制為400步，若能繼續迭代，將得到更適合的表達形式。

圖5 詳細操作步驟Fig.5 Operation details

通過擬合可以看出，分峰數為6時，擬合的曲線基本可以體現出爆炸壓力發展的整體過程，其擬合的確定數COD(R2)，也就是校正后的相關系數平方為0.995 17，如圖6所示。由圖6可知，利用高斯多峰擬合法對爆炸壓力曲線進行分析時，將所得Fit Peak按照峰中心出現時間及幅值進行排序，可得到兩類峰：一類是構成超壓曲線主體，稱為主體峰，其表征了整個超壓曲線的總體變化趨勢；另一類是構成超壓曲線出現最大值時的振蕩峰，稱為特征峰，其表征了超壓振蕩的壓力波動細節，幅值、峰寬一般都較小。對于分峰數為6時的擬合曲線來說，峰1、峰3、峰4、峰5、峰6為主體峰，峰2為特征峰。可見超壓最大時的壓力波動細節體現明顯不足，其殘差波動也比較大，因此只能滿足對爆炸超壓總體趨勢的估計，其相關參數，如表1所示。

圖6 分峰數為6的擬合情況Fig.6 The case of six fit peaks

表1 擬合曲線的相關參數Tab.1 Parameters of the fitting curve

實際上，6個高斯函數擬合后的曲線是對原壓力曲線的比較粗糙地光滑，為了更逼近原始曲線，嘗試采用采用12個高斯函數疊加來對其進行擬合，得到的結果，如圖7所示。

通過擬合可以看出，分峰數為12時，擬合的曲線更加接近爆炸壓力發展過程，其擬合的確定數COD(R)為0.998 46。其中峰1、峰2、峰5、峰6、峰8、峰9、峰10、峰11、峰12為主體峰，峰3、峰4、峰7為特征峰。隨著主體峰數量的增加，壓力的整體變化趨勢愈加精確。而特征峰增加不明顯，雖然超壓最大時的壓力波動細節比6峰時要好，殘差波動也比6峰時小，但是壓力波動細節的體現仍然略顯不足。其相關參數如表2所示。

為了更加體現壓力波動細節，達到增加特征峰的目的，在超壓曲線達到最大值時多選取幾個分峰點，并采用20個高斯函數對其進行擬合，得到的結果如圖8所示。

表2 擬合曲線的相關參數Tab.2 Parameters of the fitting curve

圖7 分峰數為12的擬合情況Fig.7 The case of twelve fit peaks

圖8 分峰數為20的擬合情況Fig.8 The case of twenty fit peaks

表3 擬合曲線的相關參數Tab.3 Parameters of the fitting curve

分峰數為20時，擬合的曲線幾乎與原爆炸壓力曲線沒有多少區別，其擬合的確定數COD(R2)為0.999 4，得到的擬合曲線表達式為

式中：Fit Peak3在該擬合段的值都為0，可以忽略，峰1、峰2、峰3、峰4、峰12、峰14、峰15、峰17、峰18、峰19、峰20為主體峰，峰5～峰11、峰13、峰16為特征峰，因此其壓力波動細節明顯比12峰時更好，殘差波動也比12峰時小。但遺憾的是，由于Origin軟件迭代步數的限制，使得結果仍不能完全達到預期效果，出現了負向波動的情況。若能繼續迭代的話，可得到更加精確的結果，將原曲線精確得劃分為主體峰與特征峰。擬合曲線的相關參數，如表3所示。

4.1 分峰數量與擬合效果的關系

通過擬合結果可以看出，最終擬合曲線與尋峰的峰軸位置及數目相關；在滿足迭代步數足夠多的前提下，隨著分峰數目的增加，擬合曲線的殘差值的范圍會縮小，擬合精度會越高，分出來的主體峰與特征峰會更加精確。當峰的數量較多時，高斯單峰之間的間距越少，其主體峰的擬合方式越接近連續變換，其峰值、峰軸、峰寬也可進行二次擬合得到面積函數A(t)、峰寬函數w(t)、峰軸函數Xc(t)。也就是說，只要高斯單峰的基數量n足夠多，則擬合出來的曲線可以無限接近試驗測得的爆炸特性曲線，其表達式可以寫成

(7)

式中：pi(t)為特征峰函數；Pz(t)為主體峰復擬合函數，其積分的具體表達式為

(8)

但是，通過試驗發現，隨著分峰數的增加，相鄰峰軸越來越接近，其產生負向峰的數量也會增加，擬合的成功率會降低，所需要的計算時間及迭代步數會大大增加，擬合的優化效果變化越來越不明顯，越難得到一個理想的擬合結果。因此，我們只需要在滿足精確度的前提下，找到一個合理的峰數即可。

4.2 分峰參數與爆炸壓力變化過程的對應關系

高斯單峰函數本身的幾個參數都具有顯著的幾何含義，其中，峰的高度Height計算為

(9)

(10)

半峰寬的計算為

(11)

峰的高度在壓力-時間曲線上表示的含義是單峰的最大超壓值；半峰寬的大小則與單峰的升壓速率相關，同時也表征了擬合峰的作用時間。

如圖9所示，以6個高斯函數擬合為例，采用高斯函數對爆炸壓力-時間曲線進行擬合，就是對爆炸模型的簡化，將爆炸的整個過程看作是多個不同時刻的不同幅度的單個壓力脈沖的疊加。也就是說，采用脈沖疊加的方式可以模擬出接近實際爆炸的超壓等效過程。

圖9 6高斯函數擬合的三維圖Fig.9 The three-dimensional curves of the case of six Gaussian peaks

4.3 單峰積分與爆炸沖量的關系

對每一個高斯函數進行積分，可以得到高斯峰面積。由于基線從0時刻開始，小于0時刻的積分幾乎可以忽略不計，因此高斯峰總面積可以近似表示為

(12)

也就是說，爆炸產生的超壓總沖量近似等于單個峰的超壓沖量之和，而高斯函數的積分是比較容易求得的

(13)

因此，利用高斯多峰擬合法可以非常快速的求出單峰的超壓沖量，進而估計出整個爆炸的總沖量。

5 結 論

為了更加細致、量化地研究可燃氣體爆炸超壓振動的主要特性，本文從試驗的角度出發，采用高斯多峰擬合法對試驗所得的爆炸超壓曲線進行了分析，得到如下結論:

(1)擬合結果的精確度與尋峰的峰軸位置及數目相關。當用12個高斯函數去擬合時，擬合曲線的主體變化趨勢比6個高斯函數擬合時更精確；當增加超壓曲線達到最大值時的分峰數時，擬合曲線的波動細節體現更明顯，其精確度高達99.94%。

(2)高斯多峰擬合法能得到比較精確的密閉空間爆炸超壓函數表達式，其對密閉空間爆炸超壓的預測及規律研究具有重要意義。

(3)高斯多峰擬合能得到爆炸超壓波動的等效模型，采用多個不同時刻不同幅值的脈沖疊加，可以等效實際爆炸的超壓過程。

(4)高斯多峰擬合法得到的單峰的面積積分之和就是整個爆炸超壓的沖量大小，采用此方法能快速地估計出整個爆炸超壓的總沖量。

在實際應用時，應根據實際需要適當調整計算時間及迭代步數等參數，提高擬合效率。