基于高斯混合模型的敞開式質(zhì)譜重疊峰解析方法

黃安貽， 繆東升， 聞路紅， 郭 榮， 胡舜迪， 洪歡歡， 吳 帥

(1.武漢理工大學(xué)機電工程學(xué)院， 武漢 430070； 2.寧波大學(xué)高等技術(shù)研究院， 寧波 315211)

敞開式質(zhì)譜是一種無需或僅需簡單樣品前處理過程，可在敞開環(huán)境下直接對樣品實現(xiàn)離子化并進行分析的新型質(zhì)譜技術(shù)，該技術(shù)能滿足實時、快速的分析需求，同時兼具傳統(tǒng)質(zhì)譜分析高靈敏度、高通量等特點。敞開式質(zhì)譜在開放環(huán)境工作時分子離子反應(yīng)較復(fù)雜[1]、易受環(huán)境因素或復(fù)雜進樣中其他碎片離子的影響[2]，由于儀器分辨率較低，在質(zhì)譜圖中常出現(xiàn)重疊峰現(xiàn)象，為了快速、準(zhǔn)確識別待測物特征峰，解析敞開式質(zhì)譜重疊峰非常必要。

在質(zhì)譜圖的重疊峰解析中，常采用色譜和質(zhì)譜聯(lián)用技術(shù)或者高分辨率儀器來完成，而通過信號處理技術(shù)來解析敞開式質(zhì)譜重疊峰的方法則有限[3]。解析重疊峰的方法主要有兩大類[4]：一類是利用數(shù)學(xué)變換的方式對原始信號進行處理來分離重疊峰，如傅里葉自去卷積，小波變換等。傅里葉自去卷積技術(shù)中合適的去卷積函數(shù)和截至函數(shù)不僅能有效抑制負(fù)旁瓣效應(yīng)的產(chǎn)生，同時能提高解析后的信噪比和分辨率[5]，此外，還能較好地解析近乎完全重疊信號，但是該方法在重疊峰不對稱時構(gòu)建去卷積函數(shù)困難而解析效果較差[4]；基于小波變換的重疊峰分析方法突破傅里葉分析時頻單一的局限，充分利用基函數(shù)的時頻局部特性在分離信號中噪聲和有用信息的同時有效提高了解析后的分辨率[6-8]，而且還能完成不同重疊程度的信號解析，但是與原始峰形相比解析后的峰強有所增大、峰寬變窄[3]，從而影響提取離子流(extracted ion chromatogram，EIC)及檢出限的設(shè)定。另一類是建立重疊子峰的數(shù)學(xué)模型，通過估計各子峰最優(yōu)參數(shù)來實現(xiàn)分離，如曲線擬合等。曲線擬合技術(shù)基于最小二乘法原理使合成信號與實測重疊信號的誤差平方和達(dá)到最小，實現(xiàn)徹底分離重疊峰的同時可獲取各子峰的峰形特征，但是噪聲干擾容易導(dǎo)致含糊的解析結(jié)果或者無法得到收斂的解[9-10]。

作為傳統(tǒng)的質(zhì)譜檢測方法之一，EIC表示一定寬度內(nèi)所有信號強度之和，常用來計算信噪比。EIC對譜峰的峰強和峰寬等峰形特征十分敏感，為了提高質(zhì)譜檢測的準(zhǔn)確性和信噪比，基于高斯分布函數(shù)是描述質(zhì)譜峰常用函數(shù)之一[11]，以及高斯混合模型(gaussian mixture model, GMM)可逼近源信號的概率密度函數(shù)[12-13]、受譜峰重疊程度影響較小等優(yōu)點，提出基于GMM的敞開式質(zhì)譜重疊峰的解析方法，改變重疊峰的幅值比、分離度和噪聲后進行解析，結(jié)果未改變峰形特征，并提高了信噪比。

1 高斯混合模型解析重疊峰理論

質(zhì)譜峰形常用高斯分布函數(shù)來描述，故將原始重疊峰信號歸一化為概率密度函數(shù)之后可用GMM模型來逼近，即構(gòu)建重疊峰信號高斯混合的參數(shù)化模型，然后通過期望最大(expectation maximization, EM)算法估計模型中各高斯分布的參數(shù)[14-15]，從而完成重疊峰的解析。采用GMM解析重疊峰流程如圖1所示。

圖1 重疊峰解析流程Fig.1 Flow of overlapping peaks analysis

GMM解析重疊峰的具體描述如下。

(1)將測量所得敞開式質(zhì)譜數(shù)據(jù)作基線校正扣除背景噪聲，獲得較為干凈的譜圖。

(2)在目標(biāo)質(zhì)荷比(mass-to-charge ratio,m/z)位置提取重疊峰信號，將信號強度值歸一化作為概率密度函數(shù)，構(gòu)建高斯混合模型，然后由離散直接抽樣按照該概率密度函數(shù)產(chǎn)生相應(yīng)的隨機數(shù)。

(3)采用EM算法對隨機數(shù)進行迭代運算，求出GMM模型的各個參數(shù)，即完成重疊峰的解析。

在該方法中，需要解決的關(guān)鍵問題是EM算法估計模型參數(shù)時初始值的設(shè)置，包括簇數(shù)、每一簇的均值、標(biāo)準(zhǔn)偏差和權(quán)重。有效選取初始值可以縮小EM算法的搜索空間、避免算法的局部收斂[12, 16]，使迭代計算更為高效準(zhǔn)確。

1.1 GMM模型

GMM是一種基于統(tǒng)計學(xué)的聚類模型，其基本假設(shè)為數(shù)據(jù)是由幾組不同高斯分布的隨機變量組合而成，它能準(zhǔn)確地逼近任意形狀的密度分布[17-18]。若有數(shù)據(jù)集X是來自多個高斯分布的混合體，則其概率分布模型為

(1)

式(1)中：αi為各高斯分布的權(quán)重，αi≥ 0且各高斯分布的權(quán)重和為1；θi= (μi,σi2)為均值μi、方差σi2的向量表示；k為模型中符合高斯分布的分支個數(shù)；pi(X,θi)為概率密度函數(shù)，表達(dá)式為

(2)

1.2 離散直接抽樣

直接抽樣方法是對任意給定的分布函數(shù)，產(chǎn)生其樣本的一種抽樣方法。若有離散分布的變量x1,x2, …,xn(即質(zhì)譜圖中的質(zhì)荷比)，已知對應(yīng)概率p1,p2, …,pn(即信號強度值歸一化的結(jié)果)，可計算出該組變量的累積分布函數(shù)為

(3)

式(3)中：p0=0, ∑pi=1。抽樣時直接產(chǎn)生服從[0，1]均勻分布的隨機數(shù)μ，求滿足表達(dá)式的k值，即

F(Xk-1)<μ≤F(Xk)

(4)

離散變量的第k個值xk即為欲抽取的值。

1.3 EM算法估計模型參數(shù)

EM算法是一種求解似然估計的迭代最優(yōu)化算法[19]，通常用來估計GMM中各個高斯函數(shù)的參數(shù)。EM算法是當(dāng)數(shù)據(jù)存在缺失問題時，在模型中引入隱變量之后再計算似然函數(shù)，交替迭代至對數(shù)似然函數(shù)收斂時停止，即可求解出各高斯分布的參數(shù)。

EM算法估計模型參數(shù)的具體步驟如下：

(1)變量初始化，需要初始化的參數(shù)有簇數(shù)k、每一簇的均值μ和方差σ2，以及隱變量W。在k、μ和σ2初值設(shè)定較為合理的情況下，對隱變量初值的要求大幅度降低，一般Wi,j設(shè)為1/k；而第j簇的權(quán)重根據(jù)隱變量求得，即

(5)

則權(quán)重初始值αj= 1/k。

(2)E步驟(期望)：根據(jù)均值、方差和權(quán)重參數(shù)的初始值或者上一次迭代的估計值來更新隱變量(其中第i個變量屬于第j簇的概率)，即

(6)

再根據(jù)式(5)可更新每一簇的權(quán)重αj。

(3)M步驟(最大化)：針對對數(shù)似然函數(shù)的期望值進行極大化估計，根據(jù)E步驟得到的隱變量值來更新均值和方差，第j簇的均值為

(7)

第j簇的方差為

(8)

(4)E步驟和M步驟交替迭代，直至收斂，即完成GMM的參數(shù)估計。

2 實驗研究

2.1 儀器和試劑

Craiv-110質(zhì)譜儀：寧波市華儀寧創(chuàng)智能科技有限公司；試劑：冰毒樣品，濃度100 μg/L。

2.2 質(zhì)譜條件

離子化能量70 eV；掃描速度4 000 amu/s；質(zhì)量范圍50～700 amu(1 amu為碳12原子質(zhì)量的1/12)；線性離子阱質(zhì)量分析器；分辨率：半峰全寬，亦稱半峰寬(full width at half mzxima，F(xiàn)WHM)為0.5；使用Python 3.7 編寫程序進行實驗。

2.3 GMM解析重疊峰

采用GMM解析重疊峰需要解決的關(guān)鍵問題：模型初始值的設(shè)置，包括簇數(shù)、每一簇的均值和標(biāo)準(zhǔn)偏差。對于實測數(shù)據(jù)，選取目標(biāo)質(zhì)譜峰信號附近存在干擾峰的重疊峰數(shù)據(jù)，本文選取冰毒碎片離子(m/z=119)的質(zhì)譜峰進行重疊峰的解析實驗。目標(biāo)峰的峰形一般均可用高斯分布函數(shù)描述，則實驗中采用多個高斯函數(shù)來模擬構(gòu)建不同重疊形式的信號，從而驗證GMM解析重疊峰的可行性、局限性以及抗噪聲能力。

2.3.1 簇數(shù)的確定

在GMM模型中，簇數(shù)k需提前確定。本文采用手肘法從數(shù)據(jù)本身出發(fā)來確定簇數(shù)k。手肘法的核心指標(biāo)是誤差平方和(sum of the squared errors，SSE)，定義為

(9)

式(9)中：Ci是第i簇；p是Ci中的元素點；mi是第i簇的均值。隨著簇數(shù)k的增大，每個簇的聚合程度會逐漸提高，則SSE會逐漸變小。當(dāng)k小于真實簇數(shù)時，由于k的增大會大幅增加每個簇的聚合程度，則SSE的下降幅度會較大；而當(dāng)k到達(dá)真實簇數(shù)后，再增加k所得到的聚合度變化幅度會變小，即SSE的下降幅度會驟減直至趨于平緩，也就是說SSE和k的關(guān)系圖是一個手肘的形狀，稱為手肘圖，而這個肘部對應(yīng)的k就是數(shù)據(jù)的真實聚類數(shù)。

通過上述手肘法確定重疊信號中的真實簇數(shù)k，結(jié)果如圖 2所示。選擇多種不同重疊情況的質(zhì)譜信號，手肘圖中肘部對應(yīng)的位置平均簇數(shù)k=2處，故模擬重疊峰數(shù)定為2。

圖2 實測重疊峰與手肘圖Fig.2 Measured overlapping peaks and elbow diagram

2.3.2 均值和標(biāo)準(zhǔn)偏差的確定

EM算法對初值較為敏感。經(jīng)驗值或隨機值可能會導(dǎo)致EM算法的局部收斂。為了避免該問題，可根據(jù)譜圖在化學(xué)量測中的物理意義來設(shè)置初值，即質(zhì)譜圖中的信號峰位、峰高和峰寬與高斯函數(shù)的均值、幅值和標(biāo)準(zhǔn)偏差有關(guān)。如圖3所示，采用模擬重疊峰數(shù)據(jù)給出求初值示意圖，均值即為峰幅值處對應(yīng)的橫坐標(biāo)，可通過尋峰的方式確定初始均值；50%峰高處的寬度(半高全峰寬)等于2.35倍的標(biāo)準(zhǔn)偏差，則根據(jù)圖3中DE段可確定標(biāo)準(zhǔn)偏差初值。其中，A點和B點是通過尋峰所得峰頂以及對應(yīng)的橫坐標(biāo)；從B點出發(fā)，沿y=0.05(最大幅值1%)水平軸向右，與曲線的第一個交點定為C點；D點和E點分別為線段AB和AC的中點。根據(jù)三角形中位線定理得DE等于BC的一半。同理，可設(shè)置另一峰的初值。

圖3 確定均值和標(biāo)準(zhǔn)偏差初值示意圖Fig.3 Schematic diagram of determining the initial value of the mean and standard deviation

3 實驗結(jié)果與討論

為了評價GMM模型結(jié)合EM算法解析重疊峰的效果，引入相對誤差、相關(guān)系數(shù)R2和信噪比三個指標(biāo)。由于已知模擬重疊峰數(shù)據(jù)的均值和標(biāo)準(zhǔn)偏差等參數(shù)，故本文采用解析前后各參數(shù)的相對誤差和R2來評價模擬數(shù)據(jù)解析的效果；對于實測數(shù)據(jù)，解析前的均值和標(biāo)準(zhǔn)偏差均未知，故采用信號和噪聲的EIC比值作為信噪比以及R2來評價實測數(shù)據(jù)解析結(jié)果。

3.1 模擬重疊峰解析

模擬重疊峰是在敞開式質(zhì)譜儀實測信號的基礎(chǔ)上設(shè)計的，已知簇數(shù)為2，按照小峰與大峰的標(biāo)準(zhǔn)偏差分別為2.5和4，峰位置相差10，峰幅值比為1∶3進行模擬重疊峰信號，抽樣次數(shù)N=5×104，根據(jù)2.2節(jié)的方法設(shè)置初始均值分別為10.345和19.655，初始標(biāo)準(zhǔn)偏差分別為3和5.5，解析前后的結(jié)果如圖4所示。從圖4(b)中可以看出，解析前后的兩個單峰基本重合，相關(guān)系數(shù)R2均大于等于0.99，沒有引起較大峰寬等峰形特征的變化。由此可知，GMM描述重疊質(zhì)譜信號，然后利用EM算法估計模型參數(shù)完成解析的方法是可行的。

3.1.1 重疊峰幅值比例的影響

實驗中，在不改變峰寬的前提下，改變模擬大小峰的幅值比來研究其對解析結(jié)果的影響，解析前后參數(shù)的相對誤差如表1所示。從表1中可以看出，采用GMM解析不同重疊情形的模擬數(shù)據(jù)，解析前后的相對誤差指標(biāo)：均值的絕對值均不大于0.4%，標(biāo)準(zhǔn)偏差的絕對值均不大于2%，即解析后對原始數(shù)據(jù)峰形影響較小，文中采取的初值設(shè)置方法可避免EM算法的局部收斂；另外，解析前后各對應(yīng)曲線的相關(guān)系數(shù)都能達(dá)到0.99。因此，解析大小峰不同幅值比例的重疊情況時，誤差較小，結(jié)果可靠。

表1 重疊峰不同幅值比解析結(jié)果Table 1 Analysis results of different amplitude ratios of overlapping peaks

3.1.2 重疊峰分離度的影響

實驗中，在大小峰幅值比為3∶1、峰寬不變的前提下，改變重疊峰的分離度研究其影響。結(jié)合質(zhì)譜理論將色譜中的分離度定義為

(10)

表2 重疊峰不同分離度解析結(jié)果Table 2 Analysis results of different resolution of overlapping peaks

疊信號分離度大于1.047。

3.1.3 噪聲的影響

實驗中，在大小峰幅值比為3∶1、峰位置為20和10、標(biāo)準(zhǔn)偏差為4和2.5的模擬重疊峰基礎(chǔ)上，增加不同強度的隨機噪聲(均值為0，改變方差)，觀察了噪聲對解析結(jié)果的影響，解析結(jié)果如圖5所示，

圖5 加入噪聲后解析結(jié)果Fig.5 Analysis results after adding noise

第一列為加入噪聲前后的重疊信號；第二列為加入噪聲重疊信號的解析結(jié)果。由結(jié)果可知，本文所述方法解析重疊峰方法具有一定的抗噪能力；但是，隨著噪聲的增強，嚴(yán)重影響原始信號的峰形特征，解析重疊峰能力降低。噪聲的干擾會影響離散直接抽樣產(chǎn)生的隨機數(shù)據(jù)，同時峰形的改變影響初值的設(shè)定，進而導(dǎo)致EM迭代計算出現(xiàn)局部收斂問題，無法得到準(zhǔn)確的解析結(jié)果。

3.2 實測重疊峰解析

對于冰毒碎片離子在m/z=119處的重疊峰信號，獲取常見的不同重疊形式的數(shù)據(jù)進行解析實驗。抽樣次數(shù)N=5×104，初始均值、標(biāo)準(zhǔn)偏差設(shè)置分別通過尋峰、峰形的物理意義分別確定，解析后的結(jié)果如圖6所示。從圖6中可以看出，對于這3種不同的實測重疊峰信號，重構(gòu)后所得GMM曲線與原始實測信號基本重合，相關(guān)系數(shù)R2均大于0.99，即解析后沒有引起峰形的變化，不會影響EIC的計算，解析效果較好。

圖6 不同重疊程度的實測數(shù)據(jù)解析結(jié)果圖Fig.6 Analysis results of measured data with different degrees of overlap

對于實測數(shù)據(jù)無法得知構(gòu)建原始重疊峰的單峰均值和方差等信息，故無法計算得到其相對誤差來量化結(jié)果。質(zhì)譜儀常用判斷檢出的條件是通過信號與噪聲的EIC比值即信噪比來設(shè)定閾值。未解析重疊峰之前，為了提高結(jié)果的準(zhǔn)確性，一般計算EIC時選取的隔離寬度較小(如隔離寬度設(shè)為1，目標(biāo)峰在m/z=119，選取的隔離范圍為119±0.5)。通過GMM解析重疊峰之后可獲取完整的單峰情況，計算EIC時可擴大隔離寬度來增大該值。實測數(shù)據(jù)解析重疊峰后信噪比結(jié)果如表3所示，其中計算EIC的隔離寬度為1。根據(jù)表3中解析前后的信噪比會發(fā)現(xiàn)，解析后的目標(biāo)峰EIC在數(shù)值上增大，信噪比提高的幅度隨著不同重疊形式的信號有所區(qū)別，最高可提高10.20%。

表3 實測重疊峰解析后信噪比Table 3 Signal-to-noise ratio after analysis of measured overlapping peak

4 結(jié)論

提出了基于高斯混合模型的解析敞開式質(zhì)譜重疊峰方法。采用手肘法和質(zhì)譜圖在化學(xué)量測中各參數(shù)的意義結(jié)合三角形中位線定理設(shè)置合理的初始參數(shù)，然后改變模擬重疊信號的幅值比和分離度、在信號中加入不同強度的隨機噪聲，對模擬數(shù)據(jù)和實測不同重疊情況的數(shù)據(jù)進行解析，通過較好的解析效果驗證了該方法可解決敞開式質(zhì)譜重疊峰問題。相對于傳統(tǒng)的信號處理技術(shù)，本文方法在分離出目標(biāo)譜峰信號時不受譜峰對稱性的影響，不改變峰形特征，并且該方法具有一定的抗噪聲干擾能力。