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

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

引 言

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

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

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

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

1 實驗部分

1.1 實驗系統(tǒng)與樣品制備

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

1.2 離散極大值法

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

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

(1)

(2)

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

1.3 伴隨拐點法

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

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

圖1 伴隨拐點法示意圖

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

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

2 結(jié)果與討論

2.1 極大值法與伴隨拐點法的比較分析

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

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

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

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

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

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

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

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

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

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

圖4 呋喃妥因的理論吸收譜與實驗吸收譜對照圖

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

2.3 伴隨拐點法對不同樣品吸收峰的識別結(jié)果分析

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

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

3 結(jié) 論

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