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

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

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

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

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

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

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

質(zhì)譜峰形常用高斯分布函數(shù)來描述，故將原始重疊峰信號歸一化為概率密度函數(shù)之后可用GMM模型來逼近，即構(gòu)建重疊峰信號高斯混合的參數(shù)化模型，然后通過期望最大(expectation maximization, EM)算法估計(jì)模型中各高斯分布的參數(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)位置提取重疊峰信號，將信號強(qiáng)度值歸一化作為概率密度函數(shù)，構(gòu)建高斯混合模型，然后由離散直接抽樣按照該概率密度函數(shù)產(chǎn)生相應(yīng)的隨機(jī)數(shù)。

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

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

1.1 GMM模型

GMM是一種基于統(tǒng)計(jì)學(xué)的聚類模型，其基本假設(shè)為數(shù)據(jù)是由幾組不同高斯分布的隨機(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(即信號強(qiáng)度值歸一化的結(jié)果)，可計(jì)算出該組變量的累積分布函數(shù)為

(3)

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

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

(4)

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

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

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

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

(1)變量初始化，需要初始化的參數(shù)有簇?cái)?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ù)的初始值或者上一次迭代的估計(jì)值來更新隱變量(其中第i個變量屬于第j簇的概率)，即

(6)

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

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

(7)

第j簇的方差為

(8)

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

2 實(shí)驗(yàn)研究

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 編寫程序進(jìn)行實(shí)驗(yàn)。

2.3 GMM解析重疊峰

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

2.3.1 簇?cái)?shù)的確定

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

(9)

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

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

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

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

EM算法對初值較為敏感。經(jīng)驗(yàn)值或隨機(jī)值可能會導(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點(diǎn)和B點(diǎn)是通過尋峰所得峰頂以及對應(yīng)的橫坐標(biāo)；從B點(diǎn)出發(fā)，沿y=0.05(最大幅值1%)水平軸向右，與曲線的第一個交點(diǎn)定為C點(diǎn)；D點(diǎn)和E點(diǎn)分別為線段AB和AC的中點(diǎn)。根據(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 實(shí)驗(yàn)結(jié)果與討論

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

3.1 模擬重疊峰解析

模擬重疊峰是在敞開式質(zhì)譜儀實(shí)測信號的基礎(chǔ)上設(shè)計(jì)的，已知簇?cái)?shù)為2，按照小峰與大峰的標(biāo)準(zhǔn)偏差分別為2.5和4，峰位置相差10，峰幅值比為1∶3進(jìn)行模擬重疊峰信號，抽樣次數(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算法估計(jì)模型參數(shù)完成解析的方法是可行的。

3.1.1 重疊峰幅值比例的影響

實(shí)驗(yàn)中，在不改變峰寬的前提下，改變模擬大小峰的幅值比來研究其對解析結(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。因此，解析大小峰不同幅值比例的重疊情況時(shí)，誤差較小，結(jié)果可靠。

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

3.1.2 重疊峰分離度的影響

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

(10)

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

疊信號分離度大于1.047。

3.1.3 噪聲的影響

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

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

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

3.2 實(shí)測重疊峰解析

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

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

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

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

4 結(jié)論

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