基于GMM-ICSA的光譜重疊峰分解研究

朱丹丹，孫世成，楊惠，張燦，朱奇光，李康，史彥新，陳穎

(1.燕山大學(xué)電氣工程學(xué)院河北省測(cè)試計(jì)量技術(shù)及其儀器重點(diǎn)實(shí)驗(yàn)室，河北秦皇島066004； 2.燕山大學(xué)信息科學(xué)與工程學(xué)院河北省特種光纖與光纖傳感器重點(diǎn)實(shí)驗(yàn)室，河北秦皇島066004； 3.自然資源部水文地質(zhì)環(huán)境地質(zhì)調(diào)查中心, 地質(zhì)環(huán)境監(jiān)測(cè)工程技術(shù)創(chuàng)新中心，河北保定071051)

1 引 言

近年來(lái)，X射線熒光光譜法(X-ray fluorescence spectroscopy, XRF)作為一種準(zhǔn)確的快速檢測(cè)物質(zhì)成分的技術(shù)[1]被許多行業(yè)廣泛應(yīng)用，但由于元素能級(jí)本身的寬度和探測(cè)器的分辨率不足等影響，使元素激發(fā)出的特征峰為類高斯分布的譜峰[2]，且由XRF的原理可得，某些元素所對(duì)應(yīng)的XRF特征峰位比較接近，因此會(huì)出現(xiàn)相互重疊的現(xiàn)象，這將對(duì)組分的定性辨識(shí)與定量分析造成嚴(yán)重的影響。

目前，處理光譜重疊峰的問題大多采取2種方法：幾何法[3]與信息處理法[4]。幾何法可以在重疊峰位不密集時(shí)，從原始光譜直接獲得波峰、波谷、切線、垂線等幾何信息來(lái)完成分峰過程，整體比較直觀且快速，但當(dāng)譜峰重疊嚴(yán)重時(shí)，計(jì)算誤差可能較大。信息處理法主要是借助小波變換、神經(jīng)網(wǎng)絡(luò)[5]、因子分析[6]、元啟發(fā)式算法等來(lái)完成分峰過程，計(jì)算誤差雖然較小，但計(jì)算速度較慢且可能遇到算法不收斂或陷入局部最優(yōu)的缺點(diǎn)。在重疊峰較為嚴(yán)重同時(shí)考慮模型精度的前提下，廣大學(xué)者使用信息處理法偏多，并對(duì)此問題進(jìn)行了深入的研究。黃洪全等[7]優(yōu)化了高斯混合模型(Gaussian mixture model, GMM)并與遺傳算法相結(jié)合進(jìn)行了XRF重疊峰的分解；陶維亮等[8]利用連續(xù)小波變換與神經(jīng)網(wǎng)絡(luò)的結(jié)合，分離出重疊峰的峰位和峰寬進(jìn)而分解了重疊峰；曾國(guó)強(qiáng)等[9]將遺傳算法與免疫算法的優(yōu)勢(shì)結(jié)合建立出的種群算法應(yīng)用于重疊峰的分解上獲得了良好的效果。

在此基礎(chǔ)上，在考慮重疊峰的峰形意義的前提下，選用了常用來(lái)處理多個(gè)高斯函數(shù)問題的GMM模型，同時(shí)選用了參數(shù)較少模型較為簡(jiǎn)單的烏鴉搜索算法(crow search algorithm, CSA)。但CSA的魯棒性較弱，無(wú)法與各種各樣的問題有一個(gè)較好的銜接，且面對(duì)高維問題的時(shí)候容易陷入局部最優(yōu)的缺陷[10]，因此對(duì)算法進(jìn)行了一定的改造。

本文對(duì)CSA做出了一些改進(jìn)：在種群中引入“烏鴉反哺”的特性，使得模型之間銜接更平滑。即在初始化的過程中引入常用于處理XRF重疊峰的GMM模型的結(jié)果值，將其設(shè)定為種群中的年老烏鴉，意識(shí)概率設(shè)為0使其更容易被其他種群跟蹤；引入了梯度型的意識(shí)概率，使種群迭代更具有多樣性，增強(qiáng)了算法全局優(yōu)化的能力，并且減少了一個(gè)可調(diào)參數(shù)使得算法的復(fù)雜度降低；改進(jìn)了全局的優(yōu)化策略，使算法更快速與穩(wěn)定。將改進(jìn)型烏鴉算法(improved the crow search algorithm, ICSA)與其他同類型的4種算法相比，其迭代時(shí)間與重構(gòu)峰后的均方誤差均優(yōu)于其他算法。

2 算法與原理

2.1 GMM模型及其參數(shù)估計(jì)

以混合土壤的X射線熒光光譜為例，其由高斯峰擬合效果如圖1所示。其中，橫坐標(biāo)為激發(fā)能量，單位為keV；縱坐標(biāo)為計(jì)數(shù)率，單位為s-1,即每秒收到的被激發(fā)熒光次數(shù)。高斯函數(shù)擬合誤差較低，可以用于替代原始光譜的離散點(diǎn)，所以可將高斯函數(shù)的性質(zhì)與原始光譜結(jié)合起來(lái)，即高斯函數(shù)的位置參數(shù)和離散程度可與光譜的峰位與峰寬對(duì)應(yīng)起來(lái)。同時(shí)由Moseley定律可得，特征X射線能量與原子序數(shù)成正比[1]：

v=Q(Z-σ)2

(1)

式中：v為對(duì)應(yīng)的X射線能量；Q為比例常數(shù)；Z為原子序數(shù)；σ為屏蔽系數(shù)。

由式(1)可得，As元素與Pb元素激發(fā)特征峰的能量峰位分別為：能量為10.542 keV的As元素峰；能量為10.506 keV的As元素峰和能量為10.448 keV的Pb元素峰，如圖1所示。因此光譜重疊峰分解的問題的關(guān)鍵點(diǎn)可以轉(zhuǎn)化為如何正確的估計(jì)在重疊峰中每個(gè)類高斯峰的參數(shù)。

圖1 高斯函數(shù)擬合光譜效果圖Fig.1 Gaussian function fitting spectrum effect diagram

由于重疊峰屬于高度重合狀態(tài)，使用信息處理法較為合理而高斯混合模型是處理多個(gè)混合高斯函數(shù)問題的常用手段，因此可以應(yīng)用于重疊峰的分解問題。GMM模型[11]可表示為

(2)

p(x=xi)=

(3)

(4)

GMM模型的參數(shù)常使用最大期望算法(Expectation Maximum, EM)來(lái)估計(jì)，EM算法初始化參數(shù)模型后，會(huì)運(yùn)用最大似然與Jensen不等式的原理來(lái)進(jìn)行參數(shù)估計(jì)，具體分為E步與M步，E步為求取數(shù)據(jù)j對(duì)于每個(gè)高斯峰i的最佳分布概率[12]：

(5)

M步為更新參數(shù)的步驟，應(yīng)找到使得似然函數(shù)最大化的模型參數(shù)，具體迭代方法為

(6)

(7)

(8)

EM算法用于估計(jì)GMM模型的參數(shù)十分方便，但EM算法也存在一些問題，如EM算法的性能依賴初始的先驗(yàn)知識(shí)且容易收斂到局部極大值[13]，所以后續(xù)考慮用元啟發(fā)式算法進(jìn)行進(jìn)一步的優(yōu)化。

2.2 ICSA算法

CSA算法的靈感即來(lái)自于烏鴉存儲(chǔ)與盜竊食物的行為[14]。本文選用CSA作為優(yōu)化模型的原因在于CSA中要設(shè)置的參數(shù)個(gè)數(shù)較少，特定參數(shù)只有意識(shí)概率Ap與飛行距離fl。較少的參數(shù)使得整體算法的復(fù)雜度較低，方便使用者應(yīng)用與調(diào)試，但CSA的不平衡搜索策略和過早收斂問題[15]使得優(yōu)化算法容易陷入局部最優(yōu)值，基于此后續(xù)對(duì)CSA做出了一些改進(jìn)，并提出了ICSA。

ICSA沿用了原算法對(duì)烏鴉種群的假設(shè),同時(shí)ICSA做出了一些改進(jìn)：

(1) 在種群初始化中引入“烏鴉反哺”的特性，將GMM模型的參數(shù)估計(jì)結(jié)果作為種群中的年老烏鴉，將其Ap值設(shè)為0，使得年老烏鴉更容易被跟蹤，即可在GMM模型的最優(yōu)值附近產(chǎn)生新的值，從而實(shí)現(xiàn)跳出局部最優(yōu)值的效果，且減少了一定的隨機(jī)初始化值使得算法更容易收斂。

(2) 將CSA中的固定值A(chǔ)p設(shè)為梯度型的Ap，使得后續(xù)計(jì)算種群適應(yīng)度時(shí)，可以將較優(yōu)的烏鴉個(gè)體分配給較低的Ap值，使得較好的優(yōu)化結(jié)果所代表的烏鴉更容易被跟蹤，從而提高產(chǎn)生更優(yōu)結(jié)果的可能性。

(3) 優(yōu)化了位置矩陣更新的方式，現(xiàn)在如果烏鴉j發(fā)現(xiàn)烏鴉i在跟蹤自己，烏鴉j則會(huì)選擇當(dāng)前的全局記憶最優(yōu)值的附近而不是產(chǎn)生一個(gè)隨機(jī)的地點(diǎn)，使得算法不容易陷入局部最優(yōu)值。

具體步驟如下：

步驟1 初始化種群參數(shù)。設(shè)立種群大小N、最大迭代次數(shù)nmax、飛行距離fl和階梯型的Ap矩陣。Ap矩陣的設(shè)立方式為:

(9)

步驟3 設(shè)立種群中年老烏鴉的初始化位置與記憶，并將對(duì)應(yīng)的Ap值設(shè)為0。

步驟4 計(jì)算種群的適應(yīng)度。本文關(guān)于重疊峰的適應(yīng)度函數(shù)描述為：將烏鴉對(duì)應(yīng)的高斯參數(shù)帶入高斯函數(shù)中，并由子峰疊加獲得重構(gòu)峰，從而對(duì)比原重疊峰與重構(gòu)峰的均方誤差，在限制參數(shù)的范圍下，該數(shù)值越小即表示烏鴉的適應(yīng)度越好。

步驟5 由搜索策略更新位置矩陣。搜索規(guī)則如下：

(10)

式中：r1為均勻分布的隨機(jī)數(shù)，本文限制在0至0.01之間；mbest為當(dāng)前種群中最優(yōu)烏鴉存儲(chǔ)食物的記憶地點(diǎn)，c1表示如下：

(11)

式中：c1的靈感來(lái)自Aljarah等[16]學(xué)者的研究，該數(shù)值隨著迭代次數(shù)的增加而逐漸減小，主要用于防止算法的過早收斂，防止其掉入局部最優(yōu)值。

步驟6 檢測(cè)新位置的可行性，并由下式更新記憶矩陣。同時(shí)由記憶矩陣的值按小至大排序的序列更新自我反跟蹤意識(shí)矩陣SA的順序，以此來(lái)保證適應(yīng)度較好的烏鴉獲得較低的對(duì)應(yīng)Ap值。

(12)

式中：f(x)為計(jì)算優(yōu)化參數(shù)可行度的函數(shù)；mj,n+1表示當(dāng)烏鴉j現(xiàn)在的位置好于原來(lái)存儲(chǔ)食物的位置時(shí)，將會(huì)更新存儲(chǔ)食物的位置，否則不更新其位置。

步驟7 檢查最佳適應(yīng)度是否小于最佳迭代預(yù)設(shè)值否則返回步驟5并繼續(xù)迭代至nmax。

3 結(jié)果與討論

本文以土壤污染中較為常見的重金屬As和重金屬Pb的X射線熒光光譜重疊峰問題為例，通過對(duì)原始光譜進(jìn)行預(yù)處理后，將光譜用GMM模型與EM算法進(jìn)行參數(shù)估計(jì)，得出的結(jié)果再通過ICSA優(yōu)化，用以尋求全局最優(yōu)值從而分解重疊峰。同時(shí)以均方誤差與迭代時(shí)間為評(píng)判標(biāo)準(zhǔn)，與同類算法：烏鴉搜索CSA算法、粒子群優(yōu)化算法(particle swarm optimization, PSO)、蝙蝠算法(bat algorithm, BA)、和聲搜索算法(harmony search, HS)做出對(duì)比，驗(yàn)證了ICSA的優(yōu)越性。

3.1 實(shí)驗(yàn)流程

As與Pb是土壤中存在的典型重金屬元素，且兩者的特征峰位較為接近，使得在分析2種元素時(shí)經(jīng)常受到重疊峰問題的干擾，如果不去進(jìn)行重疊峰分解方面的研究，會(huì)大大降低重金屬含量預(yù)測(cè)模型的準(zhǔn)確性。通過對(duì)重疊峰的分解來(lái)獲得不同元素的凈峰面積或峰高等特征來(lái)進(jìn)行模型的建立，會(huì)使得后續(xù)分析變得更加準(zhǔn)確。由于儀器獲得的原始光譜會(huì)有失真或噪聲等影響，因此首先需對(duì)原始光譜數(shù)據(jù)進(jìn)行預(yù)處理。其步驟如下：

1) 通過Savitzky-Golay法[17]來(lái)進(jìn)行去噪處理；

2) 通過Isolation-Forest算法[18]進(jìn)行異常值的剔除；

3) 通過線性本底法來(lái)去除光譜基底。并由Moseley等基礎(chǔ)定律[19]可確定各個(gè)元素的峰位，得出重金屬As和重金屬Pb的重疊峰問題大致在道址684至740之間，理論計(jì)算得出的3個(gè)重疊峰為：能量為10.542 KeV的As元素特征峰；能量為10.506 KeV的As元素特征峰和能量為10.448 KeV的Pb元素特征峰。

以As元素含量800 mg/kg和Pb元素含量100 mg/kg混合土壤的激發(fā)重疊峰為例，原始光譜和預(yù)處理后的光譜對(duì)比如圖2所示。

圖2 原始光譜與預(yù)處理光譜對(duì)比Fig.2 Comparison of original spectrum and preprocessed spectrum

GMM模型可以通過EM算法來(lái)進(jìn)行參數(shù)估計(jì)，進(jìn)而達(dá)到分解重疊峰的作用，但GMM模型在面對(duì)數(shù)據(jù)量較少且重疊峰較為嚴(yán)重時(shí)，分解效果并不理想，由圖1的預(yù)處理光譜放入GMM模型進(jìn)行參數(shù)估計(jì)的結(jié)果如圖3所示。其中重構(gòu)峰是由3個(gè)子峰疊加而成，重構(gòu)峰的擬合效果并不理想，均方誤差為0.661 6。

圖3 GMM模型分解重疊峰Fig.3 GMM model decomposition of overlapping peaks

將GMM模型估計(jì)的參數(shù)引入ICSA模型中作為年老烏鴉。其中種群大小N設(shè)為40，同時(shí)將GMM模型迭代出的3個(gè)子峰分別建立種群，即將3個(gè)子峰分成3個(gè)種群在一起迭代，方便更好地設(shè)置不同峰位的限制同時(shí)提高了迭代的速度。飛行距離fl設(shè)為0.002，通過MSE來(lái)確定最大迭代次數(shù)nmax的選擇最終為100，如圖4所示。ICSA的最終迭代結(jié)果如圖5所示。優(yōu)化模型結(jié)果的對(duì)比如表1所示。經(jīng)ICSA優(yōu)化后的不同元素的峰位與整體的峰面積均更加準(zhǔn)確。

圖4 迭代次數(shù)的選擇Fig.4 Selection of the number of iterations

圖5 ICSA優(yōu)化后的分解重疊峰Fig.5 Decomposition overlapping peaks after ICSA optimization

表1 不同模型結(jié)果對(duì)比Tab.1 Comparison of results of different models

3.2 實(shí)驗(yàn)結(jié)果對(duì)比與分析

元啟發(fā)式新算法有很多，有關(guān)學(xué)者均提出或改進(jìn)該類算法[20]，將本文提出的ICSA同經(jīng)典的同類算法做出對(duì)比實(shí)驗(yàn)，對(duì)比算法有CSA、PSO、HS、BA。實(shí)驗(yàn)過程為將不同算法分段均迭代至200次，以MSE和迭代時(shí)間為衡量算法優(yōu)劣的標(biāo)準(zhǔn)。其中不同算法需要調(diào)整的參數(shù)不同，PSO中應(yīng)調(diào)的參數(shù)有慣性權(quán)重、自我學(xué)習(xí)因子、群體學(xué)習(xí)因子等；HS中應(yīng)調(diào)的參數(shù)有和聲庫(kù)的大小、帶寬的大小、音調(diào)調(diào)節(jié)率等；BA中應(yīng)調(diào)的參數(shù)有音量衰減系數(shù)、搜索頻率系數(shù)、脈沖發(fā)射率等；CSA中只需要調(diào)整飛行距離與意識(shí)概率；ICSA中則只需調(diào)整飛行距離即可。

對(duì)比需調(diào)節(jié)參數(shù)個(gè)數(shù)可知，5種算法中ICSA的算法復(fù)雜度最小。為保證對(duì)比實(shí)驗(yàn)的準(zhǔn)確性，所有算法均引入了GMM模型的參數(shù)估計(jì)值，對(duì)比同類算法的結(jié)果如圖6所示。

圖6 實(shí)驗(yàn)對(duì)比圖Fig.6 Experimental comparison chart

由圖6可知，ICSA在速度與精度上均優(yōu)于其他4種同類型的算法，重疊峰的分解精度提高了4.93%。

4 結(jié) 論

針對(duì)XRF中As和Pb的重疊峰問題，提出將GMM模型與ICSA模型相結(jié)合來(lái)提高重疊峰的分解精度，通過縱向?qū)Ρ葘?shí)驗(yàn)得出使用ICSA模型使得重疊峰的分解精度提高了4.93%。同時(shí)為了進(jìn)一步驗(yàn)證ICSA的可行性，將其與同類4種算法做橫向?qū)Ρ葘?shí)驗(yàn)，得出ICSA的迭代速度與精度均表現(xiàn)較好，可以應(yīng)用于重疊峰的分解問題中，為后續(xù)XRF的物質(zhì)分析工作提供一定的參考。