基于差分進(jìn)化算法的X熒光重疊峰的分解

廖先莉, 黃進(jìn)初， 賴萬(wàn)昌, 辜潤(rùn)秋, 王廣西, 唐 琳, 翟 娟

1. 成都理工大學(xué)核技術(shù)與自動(dòng)化工程學(xué)院， 四川 成都 610059 2. 成都大學(xué)信息科學(xué)與工程學(xué)院， 四川 成都 610106

引 言

對(duì)于譜峰重疊問題， 一般用數(shù)學(xué)解析法進(jìn)行重疊譜分解， 譜峰重疊數(shù)學(xué)分解方法的研究， 對(duì)熒光譜進(jìn)一步的定量、 定性分析都有十分重要的意義， 現(xiàn)階段已有不少相關(guān)的研究報(bào)告[1-3]。 其中楊熙等提出了GMM-SDR模型和粒子群算法相結(jié)合的重疊譜的解譜方法[1]； 胡耀垓等運(yùn)用曲線擬合完成了光譜重疊峰解析方法[2]; 徐喜榮等提出了一種基于小波變換和連續(xù)Hopfield神經(jīng)網(wǎng)絡(luò)的譜圖重疊峰解析策略[3]； 目前還沒有一種算法被公認(rèn)為沒有局限性的數(shù)學(xué)解譜方法， 比如曲線擬合度不夠高， 易陷入局部收斂， 使用限制條件不易滿足。

本文在高斯混合統(tǒng)計(jì)模型的基礎(chǔ)上， 提出兩種情況下的參數(shù)模型， 利用差分進(jìn)化算法全局搜索優(yōu)勢(shì)， 得到了重疊譜的最優(yōu)分解模型。 兩種模型下的解譜結(jié)果誤差范圍內(nèi)都是有效的， 但是解譜精度卻不相同， 為類似數(shù)學(xué)解譜方法提供參考， 同時(shí)該方法的搜索速度快， 種群規(guī)模對(duì)尋優(yōu)結(jié)果的影響比較小， 結(jié)合了光譜的隨機(jī)物理特性， 保證的原譜數(shù)據(jù)的“零損失”。

1 算法與原理

1.1 統(tǒng)計(jì)模型GMM

一種融合了參數(shù)估計(jì)法和非參數(shù)估計(jì)法優(yōu)點(diǎn)的修正模型， 模型為

(1)

其中ai為各分支的權(quán)重， 且滿足

(2)

式(1)中，M為分支數(shù)，ui和σi表示第i分支的均值和標(biāo)準(zhǔn)差， 由各個(gè)分支的權(quán)重、 均值、 標(biāo)準(zhǔn)差構(gòu)成了差分進(jìn)化算法尋優(yōu)體參數(shù)， 所得的最優(yōu)解， 即為重疊峰分解后各個(gè)小峰的參數(shù)。

如果不能提前確定GMM模型中參數(shù)間的關(guān)系， 認(rèn)為其是相互獨(dú)立的， 建立模型為GMM參數(shù)獨(dú)立模型， 模型參數(shù)為θ=[a1，a2， …，aM；u1，u2， …，uM；σ1，σ2， …，σM]。

如果能夠知道重疊譜各個(gè)小峰參數(shù)間的關(guān)系， 建立模型為GMM參數(shù)關(guān)聯(lián)模型， 可以通過這種關(guān)系減少差分進(jìn)化算法尋優(yōu)個(gè)體的參數(shù)個(gè)數(shù)， 比如均值和標(biāo)準(zhǔn)差間存在線性關(guān)系σi=σ1ui/u1， 模型參數(shù)為θ=[a1，a2， …，aM；u1，u2， …，uM；σ1]， 下文中均以均值和標(biāo)準(zhǔn)差間存在的線性關(guān)系建立GMM參數(shù)關(guān)聯(lián)模型。

1.2 GMM模型參數(shù)的差分進(jìn)化算法估計(jì)

差分進(jìn)化算法是模擬自然界物種進(jìn)化原理的一種尋優(yōu)算法， 算法通過對(duì)父代個(gè)體進(jìn)行變異、 交叉操作， 生成新一代個(gè)體， 選擇子代和父代個(gè)體中滿足目標(biāo)條件且適應(yīng)度值最優(yōu)的個(gè)體作為本次尋優(yōu)的結(jié)果， 經(jīng)過多級(jí)迭代， 使搜索結(jié)果逐漸趨近全局最優(yōu)解[4]。 具體方法如下：

(1) 初始化種群

GMM參數(shù)獨(dú)立模型個(gè)體的維數(shù)為3M, 個(gè)體中元素ai，ui和σi在有限區(qū)間上隨機(jī)生成， 種群中粒子規(guī)模為NP。

GMM參數(shù)關(guān)聯(lián)模型個(gè)體維數(shù)2M+1， 個(gè)體中元素ai，ui和σi在有限區(qū)間上隨機(jī)生成， 種群中粒子規(guī)模為NP。 整個(gè)尋優(yōu)過程需滿足條件σi=σ1ui/u1， 其中i=2, 3， …，M。

為了滿足GMM模型中式(2)的條件， 首先要對(duì)模型參數(shù)ai做歸一化處理， 使得參數(shù)ai的取值范圍在0到1之間， 方便后續(xù)的運(yùn)算。 在對(duì)尋優(yōu)得到的參數(shù)進(jìn)行還原， 方便得到的模型曲線和實(shí)際譜線對(duì)比。

(2) 適應(yīng)度值的評(píng)估

種群中的每一個(gè)個(gè)體對(duì)應(yīng)一個(gè)適應(yīng)度值， 適應(yīng)度值由式(3)算得， 每一代種群中適應(yīng)度值最小的個(gè)體為這一代種群中的最優(yōu)個(gè)體， 多代尋優(yōu)后適應(yīng)度值不變時(shí)， 搜索結(jié)束， 得到的參數(shù)θ為全局最優(yōu)個(gè)體， 即為我們所求的分解模型參數(shù)。

(3)

由式(3)可知，g(i)為隨機(jī)序列(道址)對(duì)應(yīng)的計(jì)數(shù)值， 由隨機(jī)序列x(i)帶入GMM模型得到的值取對(duì)數(shù)后值為負(fù)數(shù)，P(x(j)/θ)越大f(θ)越大，Y(θ)的值越小， 所得個(gè)體越優(yōu)， 該方法充分使用了測(cè)量所得的所有隨機(jī)數(shù)據(jù)， 保證了原譜數(shù)據(jù)的“零損失”。

(3) 變異策略

采用式(4)的變異策略來更新種群個(gè)體， 其中，V(i)表示變異后的第i個(gè)個(gè)體，X(i)為變異前的第i個(gè)個(gè)體，X(r(1))為變異前的第r(1)個(gè)個(gè)體，X(r(2))為變異前的第r(2)個(gè)個(gè)體，r(1)和r(2)是[0, NP]范圍內(nèi)產(chǎn)生了不等于i的整數(shù)，F(xiàn)為變異縮放因子。

V(i)=X(i)+F(bestX-X(i))+F(X(r(1))-X(r(2)))

(4)

(4) 交叉策略

如式(5)所示， 個(gè)體中的X(i,j)經(jīng)過變異策略后變?yōu)閭€(gè)體V(i,j), 若[0, 1]之間生成的隨機(jī)數(shù)正好等于交叉概率CR， 或者1到Dim之間隨機(jī)生成的整數(shù)jrand等于j時(shí)， 新生成個(gè)體中的U(i,j)等于V(i,j), 否則U(i,j)為變異前的X(i,j), 該交叉策略確保了最新產(chǎn)生的個(gè)體中至少有一個(gè)U(i,j)=V(i,j)。

(5)

(5) 選擇策略

如式(6)所示， 當(dāng)經(jīng)過交叉、 變異后新產(chǎn)生的種群個(gè)體所對(duì)應(yīng)的適應(yīng)度函數(shù)值小于原來個(gè)體所對(duì)應(yīng)的適應(yīng)度函數(shù)值時(shí)， 下一代種群這個(gè)位置上的個(gè)體變?yōu)樽钚庐a(chǎn)生的種群個(gè)體， 否則， 下一代種群中這個(gè)位置上的個(gè)體保持不變。

(6)

(6) 終止條件

迭代次數(shù)滿， 或者最優(yōu)適應(yīng)度值Y(θ)連續(xù)多次不變時(shí) ， 算法終止， 否則返回(3)繼續(xù)搜索。

一直到滿足終止條件得到最優(yōu)個(gè)體， 得到各個(gè)分支高斯函數(shù)的權(quán)值、 均值和均方差， 完成重疊峰分解 。

2 重疊峰的分解

采用離散直接抽樣方法產(chǎn)生隨機(jī)數(shù)值x(1) ，x(2) ， …，x(N)， 由這些隨機(jī)數(shù)的統(tǒng)計(jì)分布構(gòu)成X熒光重疊峰， 重疊峰由幾個(gè)峰位十分接近的高斯峰重疊而成[5-7]， 本設(shè)計(jì)以重疊譜解譜為重點(diǎn)， 本底計(jì)數(shù)已扣除， 下面將依次介紹兩類重疊峰的分解過程。

2.1 三峰分解

圖1為三峰重疊譜， 峰位為210， 200和190， 面積分別為25 000， 15 000和10 000， 即歸一化后權(quán)重分別為0.5， 0.3和0.2， 橫坐標(biāo)為道址， 縱坐標(biāo)為計(jì)數(shù)。

采用GMM參數(shù)獨(dú)立的模型， 設(shè)種群的個(gè)體為X=[a1,a2,a3;u1,u2,u3;σ1,σ2,σ3], 各參數(shù)的取值下限[0.01, 0.01, 0.01; 160, 160, 160; 4, 4, 4], 取值上限[1, 1, 1; 250, 250, 250; 8, 8, 8], 選擇取值范圍內(nèi)的隨機(jī)數(shù)初始化種群， 種群規(guī)模50， 最大迭代次數(shù)1000， 縮放因子CR為0.1和交叉概率F為0.1。 采用GMM參數(shù)關(guān)聯(lián)的模型， 則提前知曉均值和標(biāo)準(zhǔn)差間的線性關(guān)系σi=σ1ui/u1， 種群個(gè)體為X=[a1,a2,a3;u1,u2,u3;σ1]， 取值范圍從[0.01, 0.01, 0.01; 160, 160, 160; 4]到[1, 1, 1; 250, 250, 250; 8]， 選擇取值范圍內(nèi)的隨機(jī)數(shù)初始化種群， 種群規(guī)模50， 最大迭代次數(shù)1 000， 縮放因子CR為0.1和交叉概率F為0.4。

圖1 三峰重疊譜

如表1可知， 參數(shù)獨(dú)立模型和參數(shù)關(guān)聯(lián)模型分別得到的權(quán)重最大誤差為8.15%和2%， 峰位最大誤差為0.30%和0.06%， 標(biāo)準(zhǔn)差的最大誤差為7.5%和1.35%。

表1 權(quán)重、 峰位和標(biāo)準(zhǔn)差的比較

圖2 原始譜、 GMM參數(shù)獨(dú)立曲線、 分解峰

圖2和圖3分別為兩種模型下原始重疊譜和分解譜的擬合情況。 已經(jīng)能夠明顯觀察到運(yùn)用GMM參數(shù)關(guān)聯(lián)模型對(duì)重疊譜進(jìn)行分解的精度比運(yùn)用GMM參數(shù)獨(dú)立模型對(duì)重疊譜進(jìn)行分解的精度高。

圖3 原始譜、 參數(shù)關(guān)聯(lián)GMM譜、 分解峰

2.2 四峰分解

如圖4所示， 為四峰重疊譜， 各子峰的峰位分別為200， 210， 225和240， 峰面積分別為20 000， 30 000， 30 000和20 000， 即歸一化后權(quán)重分別為0.2， 0.3， 0.3， 0.2， 重疊嚴(yán)重。

圖4 四峰重疊譜

設(shè)置種群個(gè)體50， 如果選擇GMM獨(dú)立參數(shù)的模型， 則每個(gè)種群個(gè)體有12個(gè)參數(shù)， 它們的取值范圍從[0.01， 0.01， 0.01， 0.01； 175， 175， 175， 175； 4， 4， 4， 4]到[1， 1， 1， 1； 265， 265， 265， 265； 8， 8， 8， 8]， 如果能夠提前知道重疊譜參數(shù)間的關(guān)系， 得到GMM的相關(guān)聯(lián)參數(shù)模型， 比如均值和標(biāo)準(zhǔn)差間的線性關(guān)系δi=δiui/u1， 則每個(gè)個(gè)體參數(shù)10個(gè)， 它們的取值范圍從[0.01， 0.01， 0.01， 0.01； 175， 175， 175， 175； 4]到[1， 1， 1， 1； 265， 265， 265， 265； 8]， 兩種模型尋優(yōu)的迭代次數(shù)為1 000， 交叉概率F=0.8, 縮放因子CR=0.8。

如表2可知， 參數(shù)獨(dú)立模型和參數(shù)關(guān)聯(lián)模型分別得到的權(quán)重最大誤差為8.3%和4.3%， 峰位最大誤差為0.12%和0.13%， 標(biāo)準(zhǔn)差的最大誤差為5.04%和0.45%。

表2 權(quán)重、 峰位、 均方差的比較

如圖5和圖6所示， 分別為通過獨(dú)立GMM模型和相關(guān)聯(lián)GMM模型， 對(duì)原始重疊峰進(jìn)行分解后的各峰擬合情況。 已經(jīng)能夠明顯觀察到運(yùn)用GMM參數(shù)關(guān)聯(lián)模型對(duì)重疊譜進(jìn)行分解的精度比運(yùn)用GMM參數(shù)獨(dú)立模型對(duì)重疊譜進(jìn)行分解的精度高。

圖5 原始譜、 GMM參數(shù)獨(dú)立曲線、 分解峰

圖6 原始譜、 參數(shù)關(guān)聯(lián)GMM譜、 分解譜

每次尋優(yōu)分解的效率和搜索結(jié)果的精度不盡相同， 略有差異， 相對(duì)而言， GMM參數(shù)關(guān)聯(lián)模型運(yùn)用差分進(jìn)化算法搜索結(jié)果更加穩(wěn)定， 精度也相對(duì)較高。

2.3 重疊譜峰的分解實(shí)例

取深井中的巖屑作為基樣， 加入鏑粉， 磨勻配置鏑元素含量分別為20， 5和2 μg·g-1的三種樣品， 由于巖屑中本身含有Mn和Fe元素， 當(dāng)加入Dy元素后， 三種元素的部分特征X射線會(huì)發(fā)生譜重疊現(xiàn)象， 如圖7所示， 為三種樣品用X熒光儀測(cè)得的全譜圖。

圖7 樣品能譜圖

Mn元素的Kβ系的特征X射線能量為6.49 keV， Fe元素Kα系的特征X射線能量為6.403 keV， Dy元素Lα系的特征X射線能量為6.495 keV， 由于Dy元素含量比較低， 為了能夠更好的分析譜重疊情況， 對(duì)圖進(jìn)行局部放大如圖8所示， 為三種元素的重疊峰。

圖8 實(shí)測(cè)三峰重疊譜

將X熒光儀實(shí)測(cè)譜線重疊峰相關(guān)數(shù)據(jù)取出扣除本底后導(dǎo)入matlab中， 根據(jù)建立的獨(dú)立GMM模型， 對(duì)重疊譜進(jìn)行分解， Dy元素含量為20 μg·g-1的樣品， 分解后各峰的情況如圖9所示。

獨(dú)立GMM模型對(duì)三種樣品的重疊譜進(jìn)行分解， 得到Dy元素的含量分別為19.10， 5.27和2.27 μg·g-1, 相對(duì)誤差分別為-4.5%， 5.4%和13.4%。 可見用該種方法進(jìn)行重疊峰的分解時(shí)， 隨著元素含量降低， 分解得到的結(jié)果相對(duì)誤差增大。

在已知所求元素種類的情況下， 可以提前已知峰位間的關(guān)系， 建立相關(guān)聯(lián)GMM模型， Dy元素和Mn元素的特征X射線能量差不多， Dy元素Lα系的特征X射線能量為6.495 keV， 可得式(1)中的u1=u2和u3=0.986 6u1， 模型中參數(shù)個(gè)數(shù)由9個(gè)變成了7個(gè)， 根據(jù)建立的相關(guān)聯(lián)GMM模型， 對(duì)重疊譜進(jìn)行分解， Dy元素含量為20 μg·g-1的樣品， 分解后各峰的情況如圖10所示。

圖9 原始譜、 GMM參數(shù)獨(dú)立曲線、 分解峰

相關(guān)聯(lián)GMM模型對(duì)三種樣品的重疊譜進(jìn)行分解， 計(jì)算得到Dy元素的含量分別為20.18, 4.85和2.15 μg·g-1, 相對(duì)誤差分別為0.9%， -3%， 7.5%， 可見該方法解譜得到的結(jié)果相對(duì)誤差比獨(dú)立GMM模型解譜得到結(jié)果的相對(duì)誤差小一些， 元素含量越低， 脈沖計(jì)數(shù)越少， 分解得到的結(jié)果相對(duì)誤差越大。

由實(shí)驗(yàn)可知， 兩種方法進(jìn)行重疊譜分解時(shí)， 能夠知道各峰間的關(guān)系， 建立相關(guān)聯(lián)GMM模型， 減少參數(shù)個(gè)數(shù)， 分解后的精度更高， 但是隨著元素含量的降低， 分解測(cè)量的精度降低了。

圖10 原始譜、 參數(shù)關(guān)聯(lián)GMM譜、 分解譜

3 結(jié) 論

運(yùn)用文中方法對(duì)重疊峰的分解結(jié)果表明， 兩種模型下， 均能實(shí)現(xiàn)重疊譜的分解， 從模擬仿真可知， 針對(duì)相對(duì)復(fù)雜的重疊峰兩種模型得到的分解結(jié)果精度都較高， 但是GMM參數(shù)關(guān)聯(lián)模型的分解精度比GMM參數(shù)獨(dú)立模型的分解精度明顯高一些。 從三峰重疊的分解實(shí)例來看， 由于實(shí)際測(cè)量過程中對(duì)測(cè)量結(jié)果的影響因素相對(duì)復(fù)雜一些， 分解計(jì)算得到的結(jié)果較仿真計(jì)算結(jié)果相對(duì)誤差大一些， 兩種方式的結(jié)果表明， 如果能夠提前得到各個(gè)相互重疊小峰之間的關(guān)系， 建立相關(guān)聯(lián)GMM模型， 減少尋優(yōu)參數(shù)個(gè)數(shù)， 對(duì)提高復(fù)雜峰的分解精度是非常重要的。