基于主成分回歸的低濃度α-HMX太赫茲光譜定量分析

楊潔萍,鄧 琥,劉泉澄

(西南科技大學信息工程學院,綿陽 621000)

HMX一般指奧克托今(Octogen),化學名為環四亞甲基四硝胺,分子式為C4H8N8O8。HMX是目前爆炸性能最好的一種單質猛炸藥,有α、β、γ、δ 4種晶型,實際應用中使用常溫穩定的β-HMX。大量研究表明，在不同的溫度、壓力等作用下,HMX的不同晶型之間會發生相互轉化,化學物理特性也會發生改變,影響HMX正常使用[1-3]。因此,為確保HMX的安全使用,對HMX中低濃度不穩定的 α-HMX 進行定量檢測顯得十分重要。

太赫茲波是頻率范圍為0.1～10.0 THz的電磁波,介于微波與紅外之間,兼具二者的優點,安全、高效、檢測速度快等優勢使其得到廣泛應用[4]，為HMX的相關研究提供了新思路新手段。目前,針對HMX展開的太赫茲光譜研究多為分析HMX在太赫茲頻段的吸收特性[5-6],HMX定量分析以近紅外光譜法為主,利用太赫茲光譜進行HMX含量檢測的報道較少。溫曉燕等[7]利用近紅外光譜技術,采用偏最小二乘法建立了α-HMX雜質晶型的近紅外光譜模型,預測誤差小于0.13%。但是該法并沒有同時預測混合HMX中不同晶型的含量。劉泉澄等[8]通過線性擬合對α/β型HMX混合物的太赫茲吸收譜建立預測方程,實現α/β型HMX的定量檢測,預測誤差最小為0.9%。但并未實現低濃度α-HMX的定量檢測。

現提出將主成分回歸[9](principle component regression,PCR)與太赫茲時域光譜結合起來實現HMX中β-HMX和低濃度α-HMX含量的同時測定。該方法以太赫茲光譜數據與組分濃度作為輸入和輸出,建立吸收系數與濃度之間的預測模型,實現對α/β型HMX混合物中組分含量的快速定量檢測。

1 實驗部分

1.1 實驗裝置

實驗采用Zomega公司生產的太赫茲時域光譜系統(Z-3)以及美國光譜物理公司生產的飛秒激光器(MaiTai HP)。實驗設置進入Z-3系統的平均飛秒脈沖功率為150 mW,溫度為23 ℃,濕度控制在3%之內。

1.2 樣品制備

實驗采用紅外壓片法制作樣品。聚乙烯粉末在太赫茲波段內無吸收,加入進HMX樣品中可以增強樣品黏合性,同時可以減少樣品過量帶來的強烈吸收。將聚乙烯粉末與HMX樣品按照一定比例混合,使用瑪瑙缽充分碾磨混合均勻,并在12 t的壓力下壓成直徑13 mm的薄片。

1.3 數據處理

(1)太赫茲時域光譜中存在的反射峰可能會隱藏一些重要的吸收特征[10],導致后續分析結果不準確,因此采用時域截斷的方法來去除反射峰。

(2)將去除反射峰后的時域光譜數據作快速傅里葉變換,得到包含振幅以及相位信息的樣品太赫茲頻域譜。

(3)由于低頻時噪聲的影響,相位會發生2π的跳變[11],需要校正相位使零頻相位接近于0。相位校正方法為:將0.5～1.2 THz范圍內的相位與頻率作線性關系得到線性方程,此時得到零頻相位φ,將φ除以2π取整得到N,再用所有頻點的相位減2πN。

(4)根據Dorney等[12]提出的太赫茲時域光譜技術提取光學參數的模型計算得到化合物的吸收系數譜。

2 主成分回歸算法

主成分回歸是以主成分為自變量而進行的回歸分析。它的基本思想是先通過主成分分析對數據進行降維得到主成分變量,然后對主成分變量進行回歸分析得到回歸方程,最后將原變量代入回歸方程。

設X為m個樣品在n個頻率點處的吸收系數矩陣,Y為m個樣品中p種組分的濃度矩陣。X可表示為

(1)

對X進行線性變換,得到原始變量的線性組合,即新變量為

(2)

式(2)中：Z1為新變量中方差最大的,稱為第1主成分；Z2為除Z1之外方差最大的,稱為第2主成分。以此類推,Zn為除前n個新變量之外方差最大的,為第n主成分。

一般選取前h個主成分來表征原始數據矩陣。此時建立Y與前h個主成分Z1,Z2,…,Zh的回歸方程為

yi=β0+β1Zi1+β2Zi2+…+βhZih+εi，

i=1,2,…,m

(3)

式(3)中：βi為各項系數；εi為殘差。

令

則回歸模型的矩陣形式為

(4)

采用最小二乘法尋找β=(β0,β1,…,βh)T的估計β′,使其殘差平方和SR最小,即

(Y-Cβ)T(Y-Cβ)

(5)

可得

β′=(CTC)-1CTY

(6)

由此得到Y與Z1,Z2,…,Zh的線性回歸方程。式(2)說明了Z1,Z2,…,Zh與X的對應關系,將其代入回歸方程中,可得到Y與X的線性回歸方程。

3 結果與分析

3.1 單質HMX的吸收光譜

配制了α和β型HMX的單質樣品,以此來考察單組分的吸收系數譜特性。圖1(a)所示為 α-HMX 和β-HMX的原始吸收系數譜,可以看出明顯的周期振蕩現象,因此樣品的特征吸收峰值很難分辨出。去除反射峰后,同一樣品的吸收系數譜如圖1(b)所示,此時可以明顯看出α-HMX的特征峰值大概在0.8 THz及1.5 THz左右,β-HMX的特征峰值在 1.8 THz 左右,因此可以通過選取特征峰值左右的數據點來進行主成分回歸建模。

圖1 單質HMX吸收系數譜Fig.1 HMX absorption coefficient spectrum

3.2 主成分數的確定

主成分貢獻率如圖2所示,同時考察了主成分數與主成分回歸模型總相對預報誤差(total relative prediction error,RPET)之間的關系，如圖3所示,兩者共同來確定主成分數?？傁鄬︻A報誤差公式為

(7)

式(7)中：ρred,i,j、ρreal,i,j分別為第i個混合物中的第j個組分的實際濃度和預測濃度。

從圖2可看出,第1主成分的貢獻率和前兩個主成分的累計貢獻率都很接近100%,可以選擇第1主成分或前兩個主成分來表征原始數據矩陣。從圖3中可以看出,主成分數為1時,RPET最小,主成分數為2時,RPET相對于主成分數為1時大為增加,而后隨著主成分數的增加,RPET逐漸減小,但仍高于主成分數為1時的誤差。因此,綜合考慮,選取第一主成分來表征原始數據矩陣,建立預測矩陣。

圖2 主成分貢獻率Fig.2 Principal component contribution rate

圖3 主成分與總相對預報誤差之間的關系Fig.3 The relationship between principal component number and total relative prediction error

3.3 預測模型的建立及檢驗

共制備了26組α-β HMX混合樣品,其中使用20組樣品作為訓練集,如表1所示,用以建立預測模型；剩余的6組樣品作為預測集,用以驗證所建校正模型的可靠性。

表1 混合樣品訓練集組分含量Table 1 Content of training set of mixed samples

在MATLAB R2016a的環境下建立PCR預測模型,導入訓練集數據。通過計算得出,預測模型的相關系數為0.956,建立的模型具有較好的線性相關性。將測試集數據導入預測模型中,得出樣品中各組分含量的預測值,計算相對誤差和絕對誤差，如表2所示。

表2 混合樣品預測集組分含量預測值及誤差分析Table 2 Prediction value and error analysis of component content in the prediction set of mixed samples

從表2中可以看出，樣品組分含量預測值與實際值之間存在著一定的誤差。其中,最大絕對誤差為1.098%,最小為0.146%,保持在1%范圍內；相對誤差最大為0.082,最小為0.002,保持在0.1范圍內。樣品組分含量預測值和實際值之間的誤差在可接受范圍內,因此該方法可以較好地預測HMX的組分含量。

樣品各組分含量實際值與預測值關系如圖4所示,可以看出兩種組分的預測結果和實際結果的數據點的線性相關性較好。

表2、圖4說明了預測模型能較為準確地識別出不同組分的不同含量,可以對α/β型HMX混合物的組分含量做出可靠的判斷,說明了基于主成分回歸算法的預測模型的可行性以及可靠性。

圖4 混合物各組分實際值與預測值關系Fig.4 The relationship between the actual and predicted values of each component of the mixture

4 結論

配置了不同濃度配比的α/β型HMX混合物的樣本,獲取其太赫茲時域光譜數據,用以定量分析模型的建立以及檢驗；并在MATLAB環境中編寫主成分回歸算法,建立吸收系數與濃度之間的線性預測模型,并使用檢驗集樣本驗證了該模型的精確性。結果表明,實驗建立的線性預測模型的相關系數為0.956,組分含量預測誤差保持在1%范圍內,相對誤差則保持在0.1的范圍內,對α/β型HMX混合物的定量分析具有較好的預測效果。該模型可同時檢測α/β型HMX混合物中兩種晶型含量,為低濃度α-HMX含量快速檢測提供了新方法。