基于離散極大值法的太赫茲光譜特征吸收峰提取方法

曹雨齊，康旭升1, *，陳飄云，謝 晨，喻 潔*，黃平捷，侯迪波，張光新

1.浙大城市學(xué)院計(jì)算機(jī)與計(jì)算科學(xué)學(xué)院，浙江 杭州 310015 2.浙江大學(xué)控制科學(xué)與工程學(xué)院，工業(yè)控制技術(shù)國(guó)家重點(diǎn)實(shí)驗(yàn)室，浙江 杭州 310027

引 言

太赫茲波是頻率位于0.1～10 THz(1 THz=1012Hz)的電磁波，其波長(zhǎng)范圍為0.03～3 mm。由于太赫茲波具有低能性、指紋特性、強(qiáng)穿透性等特點(diǎn)，因此太赫茲時(shí)域光譜技術(shù)被廣泛應(yīng)用于食品安全檢測(cè)[1-5]、藥品質(zhì)量檢測(cè)[6-10]、醫(yī)學(xué)診斷[11-16]等多個(gè)應(yīng)用領(lǐng)域。

太赫茲光譜中的特征吸收峰在含峰目標(biāo)物的鑒定研究中起著重要作用[17-20]。浙江大學(xué)張宏建團(tuán)隊(duì)利用四種殺蟲劑獨(dú)特的特征吸收峰強(qiáng)度對(duì)四種殺蟲劑進(jìn)行定量檢測(cè)[17]。三聚氰胺混合物在0.1～3.0 THz頻率范圍內(nèi)被觀察到三個(gè)特征吸收峰，分別位于2.00，2.26和2.60 THz處。該研究根據(jù)混合物樣品的太赫茲吸收系數(shù)譜圖中的特征峰信息對(duì)三聚氰胺進(jìn)行檢測(cè)。結(jié)果表明，太赫茲時(shí)域光譜技術(shù)在食品質(zhì)量檢測(cè)中具有潛在的應(yīng)用價(jià)值[18]。

上述研究中，樣品均具有明顯特征吸收峰，可被肉眼識(shí)別，但部分弱峰可能被忽略，從而影響含峰目標(biāo)物含量的精準(zhǔn)檢測(cè)。因此研究一種由計(jì)算機(jī)自主有效地提取太赫茲光譜中的吸收峰位特征信息的方法，是太赫茲波廣泛應(yīng)用于各領(lǐng)域物質(zhì)識(shí)別的必要性研究。

為實(shí)現(xiàn)藥品的太赫茲弱吸收峰的精準(zhǔn)識(shí)別，本文提出了伴隨拐點(diǎn)法，并將其應(yīng)用于食品添加劑硝基呋喃類化合物的太赫茲吸收峰識(shí)別中。研究結(jié)果表明伴隨拐點(diǎn)法對(duì)弱吸收峰具有較好的識(shí)別作用。

1 實(shí)驗(yàn)部分

1.1 實(shí)驗(yàn)系統(tǒng)與樣品制備

實(shí)驗(yàn)采用Z3太赫茲時(shí)域光譜系統(tǒng)(Zomega Terahertz Corporation, USA)，系統(tǒng)的詳細(xì)介紹可參考文獻(xiàn)[21]。實(shí)驗(yàn)測(cè)試的各硝基呋喃類藥物的純物質(zhì)均為結(jié)晶性粉末。為制備樣品壓片，首先將樣品粉末放置于真空干燥柜中干燥，然后在瑪瑙研缽中充分研磨。最后，將樣品粉末壓制成厚度為1 mm、直徑為13 mm的樣品壓片。由于純硝基呋喃粉末不易壓制成型，需與聚乙烯粉末按1∶1的比例進(jìn)行混合后壓片。

1.2 離散極大值法

連續(xù)極大值法是指根據(jù)連續(xù)函數(shù)曲線f(x)中x點(diǎn)的一階和二階導(dǎo)數(shù)f′(x)和f″(x)數(shù)值符號(hào)來(lái)判斷該點(diǎn)是否為特征吸收峰位。如果f′(x)=0且f″(x)<0, 則將x點(diǎn)判定為f(x)的吸收峰位。

然而由于太赫茲時(shí)域光譜系統(tǒng)采樣不連續(xù)，因此頻域譜圖為離散函數(shù)，則可能存在含峰物質(zhì)的峰位未被采樣，導(dǎo)致沒(méi)有x滿足f′(x)=0。因此，針對(duì)離散數(shù)據(jù)，判別吸收峰位的條件為對(duì)點(diǎn)x，f′(x)的取值需經(jīng)過(guò)從正值到負(fù)值的轉(zhuǎn)換，同時(shí)f″(x)<0。由于實(shí)驗(yàn)系統(tǒng)中Δx為固定的非零常數(shù)，f′(x)能夠利用中心差分公式近似計(jì)算，如式(1)。同理計(jì)算得到f″(x)，如式(2)

(1)

(2)

根據(jù)離散極大值法識(shí)別特征吸收峰的方法如下。設(shè)xi為頻段范圍內(nèi)任意頻率點(diǎn)，分別利用式(1)和式(2)近似計(jì)算吸收系數(shù)譜圖的一階和二階導(dǎo)數(shù)，若f′(xi-1)f′(xi)<0且f″(xi)<0，則判定xi為特征譜的極大值點(diǎn)或者吸收峰位。

1.3 伴隨拐點(diǎn)法

雖然離散極大值法對(duì)強(qiáng)吸收峰具有較強(qiáng)識(shí)別能力，但它對(duì)弱吸收峰的敏感性不強(qiáng)，這是因?yàn)殡x散極大值法僅根據(jù)該點(diǎn)或在鄰域內(nèi)確定峰值位置。若吸收峰位于采樣間隔Δx內(nèi)而沒(méi)有被采樣，則一階導(dǎo)數(shù)的離散曲線可能并不形成過(guò)零點(diǎn)，導(dǎo)致弱吸收峰被忽略。因此，本文基于離散極大值法提出一種能夠同時(shí)測(cè)定強(qiáng)吸收峰和弱吸收峰的吸收峰識(shí)別方法，即伴隨拐點(diǎn)法。

伴隨拐點(diǎn)法利用特征吸收峰的伴隨拐點(diǎn)信息來(lái)識(shí)別對(duì)應(yīng)的特征吸收峰，其示意圖如圖1所示，其中a，b，c，d為吸收峰，a1與a2、b1與b2、c1與c2、d1與d2分別為a，b，c，d的伴隨拐點(diǎn)對(duì)。伴隨拐點(diǎn)法的具體步驟如下：(1)獲取樣品的原始太赫茲吸收系數(shù)譜(Origin，如圖1藍(lán)色曲線所示，其中部分曲線被紅色曲線重疊覆蓋未顯示)，確定拐點(diǎn)；(2)確立拐點(diǎn)中的各伴隨拐點(diǎn)對(duì)(如圖1中的a1與a2，b1與b2)；(3)連接伴隨拐點(diǎn)對(duì)，與原始譜共同確定基線譜(Baseline，如圖1中紅色曲線所示)；(4)將原始譜與基線譜作差，確定原始譜與基線譜的差譜(Difference，如圖1中褐色曲線所示)；(5)利用離散極大值法計(jì)算差譜的吸收峰位，即為目標(biāo)檢測(cè)物的特征吸收峰。

圖1 伴隨拐點(diǎn)法示意圖

伴隨拐點(diǎn)法中伴隨拐點(diǎn)對(duì)的確定方法如下。首先，設(shè)xi為頻段范圍內(nèi)任意頻率點(diǎn)，若滿足f″(xi-1)f″(xi)<0，則判定xi為特征譜的拐點(diǎn)。然后，根據(jù)f″(xi-1)和f″(xi)的符號(hào)將拐點(diǎn)分類，即起勢(shì)拐點(diǎn)和落勢(shì)拐點(diǎn)。定義若f″(xi-1)>0且f″(xi)<0，則xi為起勢(shì)拐點(diǎn)；若f″(xi-1)<0且f″(xi)>0，則xi為落勢(shì)拐點(diǎn)。最后，將連續(xù)出現(xiàn)的一對(duì)起勢(shì)拐點(diǎn)和落勢(shì)拐點(diǎn)作為一對(duì)伴隨拐點(diǎn)。

為避免伴隨拐點(diǎn)法過(guò)敏感，本研究提出差譜中對(duì)應(yīng)峰值強(qiáng)度設(shè)定吸收峰的識(shí)別規(guī)則。具體為：如果利用“離散極大值法”未在差譜中的某伴隨區(qū)間內(nèi)識(shí)別出極大值點(diǎn)，則直接將對(duì)應(yīng)候選吸收峰排除；否則如果在該區(qū)間內(nèi)識(shí)別出極大值點(diǎn)，但在差譜中的峰值強(qiáng)度低于最大吸收峰峰值強(qiáng)度的5%或者在原始吸收譜中對(duì)應(yīng)的吸收系數(shù)在20%分位數(shù)以下，則認(rèn)為其吸收強(qiáng)度過(guò)低，也對(duì)應(yīng)候選吸收峰排除；其余識(shí)別出的極大值點(diǎn)均可認(rèn)定為原始吸收譜的吸收峰位。

2 結(jié)果與討論

2.1 極大值法與伴隨拐點(diǎn)法的比較分析

利用伴隨拐點(diǎn)法對(duì)一種易被濫用的硝基呋喃類藥物呋喃妥因進(jìn)行了特征吸收峰的識(shí)別，結(jié)果如圖2所示。為驗(yàn)證伴隨拐點(diǎn)法的有效性，同樣采用離散極大值法對(duì)該樣品進(jìn)行吸收峰位識(shí)別，并比較分析。黑色曲線為呋喃妥因在THz波段的原始吸收系數(shù)譜圖，藍(lán)色曲線為原始吸收系數(shù)譜圖的一階導(dǎo)數(shù)，紅色曲線為對(duì)應(yīng)二階導(dǎo)數(shù)。

圖2 呋喃妥因的吸收系數(shù)曲線及一、二階導(dǎo)數(shù)曲線

從圖中觀察到，該物質(zhì)至少存在三個(gè)吸收峰，峰位位于0.90，1.25和1.60 THz附近，其中1.25 THz附近為強(qiáng)吸收峰(圖2中B點(diǎn))，而0.90和1.60 THz處的吸收峰較弱(圖2中A和C點(diǎn))。對(duì)于強(qiáng)吸收峰B，在1.25 THz附近的每個(gè)一階導(dǎo)數(shù)曲線均隨頻率的增大由正值轉(zhuǎn)為負(fù)值，即1.25 THz附近恰為所有的一階導(dǎo)數(shù)曲線的“過(guò)零點(diǎn)”。同時(shí)1.25 THz附近所有的二階導(dǎo)數(shù)的取值均小于零，符合離散極大值法判別吸收峰的條件。同樣，對(duì)于弱吸收峰C，其導(dǎo)數(shù)特征也符合離散極大值法的評(píng)判指標(biāo)，能夠被離散極大值法識(shí)別到。但是，對(duì)于弱吸收峰A，其峰位0.90 THz附近，雖然二階導(dǎo)數(shù)均小于零，但在0.90 THz附近的一階導(dǎo)數(shù)均大于零，因此，離散極大值法無(wú)法識(shí)別到該吸收峰。

由此可見，離散極大值法用于判定強(qiáng)吸收峰十分有效，但用于弱吸收峰的判定時(shí)，有效性低，所以在使用此法尋找吸收峰位時(shí)有可能會(huì)錯(cuò)過(guò)一些特征吸收峰，并不利于物質(zhì)的識(shí)別和鑒定。而伴隨拐點(diǎn)法恰好能滿足這些需求。

呋喃妥因的吸收系數(shù)譜圖及伴隨拐點(diǎn)法的峰位識(shí)別結(jié)果如圖3所示，其中藍(lán)色圓點(diǎn)和紅色圓點(diǎn)分別為用于構(gòu)造伴隨區(qū)間的起勢(shì)拐點(diǎn)和落勢(shì)拐點(diǎn)。圖3中伴隨區(qū)間的數(shù)量達(dá)到7個(gè)，對(duì)應(yīng)的起勢(shì)拐點(diǎn)和落勢(shì)拐點(diǎn)的數(shù)量達(dá)到14個(gè)。但是根據(jù)差分光譜的譜峰識(shí)別，最終識(shí)別到4個(gè)吸收峰，分別為A，B，C和D。值得注意的是，伴隨拐點(diǎn)法不僅能識(shí)別弱吸收峰A，而且能識(shí)別到另一個(gè)肉眼無(wú)法識(shí)別的弱吸收峰D。伴隨拐點(diǎn)法計(jì)算得到的差譜中有兩個(gè)明顯的峰和兩個(gè)相對(duì)較弱的峰。因此，根據(jù)差分光譜，兩個(gè)強(qiáng)吸收峰分別對(duì)應(yīng)B和C，兩個(gè)相對(duì)較弱的吸收峰分別對(duì)應(yīng)A和D。差分光譜中的峰位代表了利用伴隨拐點(diǎn)法識(shí)別到的呋喃妥因的吸收峰位。另外，差分光譜中B峰的強(qiáng)度最強(qiáng)，其次是吸收峰C，這與離散極大值法得到的結(jié)果一致。因此，伴隨拐點(diǎn)法不管是對(duì)強(qiáng)吸收峰還是對(duì)弱吸收峰的識(shí)別，均十分有效。

圖3 呋喃妥因的吸收系數(shù)譜圖及伴隨拐點(diǎn)法的峰位識(shí)別結(jié)果

而對(duì)于候選弱吸收峰E，F(xiàn)和G而言，其強(qiáng)度在差譜中不清晰。根據(jù)峰位判別規(guī)則，利用離散極大值法未在差譜中識(shí)別出E和F點(diǎn)，則直接將其排除。同時(shí)，由于G峰值在差譜中強(qiáng)度低于吸收峰B的峰值強(qiáng)度的5%，也排除吸收峰G，其余識(shí)別出的極大值點(diǎn)(A，B，C，D)均可認(rèn)定為原始吸收譜的吸收峰位。

2.2 基于密度泛函理論的呋喃妥因吸收峰的理論計(jì)算

為驗(yàn)證伴隨拐點(diǎn)法對(duì)弱吸收峰的識(shí)別能力，將吸收峰識(shí)別結(jié)果與仿真結(jié)果進(jìn)行對(duì)比分析。本研究利用密度泛函理論針對(duì)呋喃妥因的吸收峰振動(dòng)模式進(jìn)行了理論解析。輸入呋喃妥因分子的三維結(jié)構(gòu)文件，在Chem3D軟件中采用半經(jīng)驗(yàn)分子軌道PM3方法，獲得呋喃妥因分子的初始最小能量結(jié)構(gòu)，然后在Gaussian量子化學(xué)軟件中采用密度泛函理論中的B3LYP方法，選取6-311G基組水平對(duì)此初始分子結(jié)構(gòu)進(jìn)行了幾何優(yōu)化和頻率計(jì)算。計(jì)算結(jié)果沒(méi)有虛頻出現(xiàn)，說(shuō)明幾何優(yōu)化計(jì)算得到了分子的最小能量結(jié)果。

圖4為呋喃妥因在0.2～1.8 THz范圍內(nèi)的實(shí)驗(yàn)光譜與利用密度泛函理論計(jì)算的理論光譜對(duì)照?qǐng)D，藍(lán)色曲線和紅色曲線分別表示呋喃妥因的實(shí)驗(yàn)光譜和理論光譜。理論計(jì)算結(jié)果顯示，呋喃妥因在0.2～1.8 THz范圍內(nèi)共有三個(gè)吸收峰，分別位于0.89，1.30和1.67 THz處。由于理論光譜中的吸收系數(shù)與實(shí)驗(yàn)光譜不在一個(gè)數(shù)量級(jí)上，為便于觀察，理論光譜中的各縱坐標(biāo)采用了同時(shí)放大為1 000倍后的數(shù)據(jù)，其中放大結(jié)果不影響吸收峰位的指認(rèn)。另外，理論光譜中1.30 THz處的吸收峰的峰值與其余兩個(gè)吸收峰的峰值相比，相對(duì)較小，無(wú)法在圖4中觀察到，但我們?nèi)栽趫D4中做了標(biāo)記，方便對(duì)比。

圖4 呋喃妥因的理論吸收譜與實(shí)驗(yàn)吸收譜對(duì)照?qǐng)D

呋喃妥因的太赫茲吸收系數(shù)光譜中利用伴隨拐點(diǎn)法識(shí)別出的在1.25 THz處的強(qiáng)吸收峰，和在0.90和1.60 THz處的兩個(gè)弱吸收峰均有對(duì)應(yīng)的理論計(jì)算結(jié)果，分別對(duì)應(yīng)于理論光譜中1.30，0.89和1.67 THz處的吸收峰。而其中實(shí)驗(yàn)光譜0.90 THz處的弱吸收峰是無(wú)法利用離散極大值法識(shí)別出來(lái)的，理論光譜中0.89 THz的吸收峰進(jìn)一步驗(yàn)證了利用伴隨拐點(diǎn)法識(shí)別實(shí)驗(yàn)光譜中吸收峰的合理性。不過(guò)實(shí)驗(yàn)光譜中利用伴隨拐點(diǎn)法識(shí)別出的在1.48和1.72 THz處的弱吸收峰未在理論光譜中找到對(duì)應(yīng)的吸收峰，這與理論計(jì)算模擬軟件所采用分子模擬條件有關(guān)。

2.3 伴隨拐點(diǎn)法對(duì)不同樣品吸收峰的識(shí)別結(jié)果分析

為進(jìn)一步驗(yàn)證伴隨拐點(diǎn)法的有效性，利用該方法對(duì)呋喃妥因、呋喃西林、呋喃他酮和呋喃唑酮的吸收峰位也進(jìn)行識(shí)別，結(jié)果如圖5所示。每個(gè)樣本吸收系數(shù)譜圖均有多對(duì)伴隨拐點(diǎn)，其中藍(lán)點(diǎn)是起勢(shì)拐點(diǎn)，紅點(diǎn)表示落勢(shì)拐點(diǎn)。以呋喃妥因?yàn)槔?，共識(shí)別出10對(duì)拐點(diǎn)對(duì)，但是最終識(shí)別出強(qiáng)度相對(duì)較大的5個(gè)吸收峰，分別為1.25 THz處的強(qiáng)吸收峰，0.90，1.48，1.60和1.72 THz處的四個(gè)弱吸收峰。其中0.90，1.48和1.72 THz處的弱吸收峰采用離散極大值法均無(wú)法識(shí)別。從結(jié)果中觀察到，結(jié)合伴隨拐點(diǎn)法的吸收峰位識(shí)別規(guī)則，能夠有效避免過(guò)敏感現(xiàn)象，實(shí)現(xiàn)含峰目標(biāo)物吸收峰位的精準(zhǔn)識(shí)別。

圖5 硝基呋喃類藥物吸收系數(shù)的伴隨拐點(diǎn)法結(jié)果

3 結(jié) 論

提出了一種基于離散極大值法的含峰目標(biāo)物的峰位檢測(cè)方法伴隨拐點(diǎn)法，并將其應(yīng)用于識(shí)別和提取硝基呋喃類化合物的太赫茲吸收峰位。同時(shí)比較了離散極大值法和伴隨拐點(diǎn)法提取的樣品峰位信息，觀察到伴隨拐點(diǎn)法能夠有效識(shí)別離散極大值法無(wú)法識(shí)別的弱吸收峰，但同時(shí)具有較好的避免過(guò)敏感的能力。另外，將密度泛函理論計(jì)算得到的理論光譜與太赫茲譜圖的識(shí)別結(jié)果進(jìn)行比較，進(jìn)一步驗(yàn)證了伴隨拐點(diǎn)法對(duì)位于0.89 THz的弱吸收峰位識(shí)別的有效性。該結(jié)果表明，伴隨拐點(diǎn)法在光譜譜圖解析和物質(zhì)檢測(cè)等領(lǐng)域具有潛在的應(yīng)用前景。