基于密度泛函理論與自舉軟縮減法的酒石酸太赫茲光譜特征譜區(qū)分析指認(rèn)

唐 鑫， 周勝靈， 祝詩平， 馬羚凱， 鄭 權(quán)， 普 京

西南大學(xué)工程技術(shù)學(xué)院， 重慶 402160

引 言

太赫茲波(Terahertz, THz)頻率介于0.1～10 THz之間， 作為一門銜接了微觀量子理論和經(jīng)典電磁波理論的新興科學(xué)， 展現(xiàn)了許多獨(dú)特的性質(zhì)， 如安全性、 透視性和較強(qiáng)的分辨能力。 THz光譜對固態(tài)分子的排列、 振動(dòng)和轉(zhuǎn)動(dòng)非常敏感[1]， 可以通過振幅和相位信息， 對材料的成分及其結(jié)構(gòu)的細(xì)微變化進(jìn)行有效的分析。

目前， 國內(nèi)外許多學(xué)者針對THz光譜信息篩選特征譜區(qū)， 簡化模型展開了研究。 Hu等[2]利用競爭性自適應(yīng)加權(quán)采樣(CARS-PLS)對苯甲酸在1.0～3.0 THz頻段建立光譜模型， 得到相關(guān)系數(shù)為0.9953， 較傳統(tǒng)方法， 檢測精度大幅度提高。 Li等[3]利用蒙特卡洛無信息變量消除法(MC-UVE-PLS)提取近紅外光譜信息， 用于測定棉籽中棉酚的含量， 結(jié)果表明， MC-UVE-PLS算法可以通過提取有效信息和剔除無關(guān)變量達(dá)到特征譜區(qū)篩選的目的。 Jiang等[4]通過間隔區(qū)間偏最小二乘法(iPLS)實(shí)現(xiàn)對0.2～1.6 THz的山梨酸和山梨酸鉀樣品THz光譜的特征譜區(qū)篩選， 得到了很好的預(yù)測精度。 自舉軟縮減法(bootstrapping soft shrinkage， BOSS)是由Deng于2016年提出[5]， 該算法基于隨機(jī)抽樣統(tǒng)計(jì)技術(shù)(bootstrap sampling， BSS)和加權(quán)引導(dǎo)采樣技術(shù)(weight bootstrap sampling, WBS)， 用于生成變量的隨機(jī)組合并構(gòu)建子模型， 對不同的變量賦予不同的權(quán)重。 基于生成的權(quán)重模型， 采用模型集群分析(MPA)結(jié)合PLS提取信息并建模分析。 對上述步驟進(jìn)行反復(fù)的變量迭代提取， 簡化變量空間， 權(quán)重較大的變量有更大概率得以保留。 通過迭代計(jì)算中最小的交叉驗(yàn)證均方根誤差(RMSECV)， 確定最優(yōu)變量集。 對變量中存在共線性信息導(dǎo)致的模型不準(zhǔn)確問題有良好的處理效果。

上述方法， 多是從提升預(yù)測結(jié)果出發(fā)， 通過迭代優(yōu)化參數(shù)進(jìn)行特征譜區(qū)篩選， 但對篩選的特征譜區(qū)是否合理缺乏理論支撐， 因此本文將通過量子化學(xué)方法對特征峰在分子結(jié)構(gòu)振動(dòng)的歸屬予以說明， 對篩選的特征譜區(qū)正確性加以驗(yàn)證。 隨著計(jì)算機(jī)性能的顯著提高， 利用量子化學(xué)方法對分子振動(dòng)模式進(jìn)行仿真， 分析光譜指紋特征的產(chǎn)生機(jī)理， 得到了廣泛應(yīng)用。 密度泛函理論(density functional theory， DFT)通過粒子密度函數(shù)來描述原子和分子等基態(tài)的特征， 常用于計(jì)算和預(yù)測分子的物理結(jié)構(gòu)和化學(xué)性質(zhì)， 基于其顯著的“魯棒性”， 以簡便有效的運(yùn)算為分子特征提供合理準(zhǔn)確的預(yù)測[6-7]。 Zhou等[8]利用DFT對卡馬西平與煙酰胺、 糖精和富馬酸的單個(gè)組分的晶體振動(dòng)進(jìn)行了理論研究并得到THz光譜。 實(shí)驗(yàn)結(jié)果顯示， 仿真的THz光譜成功再現(xiàn)了實(shí)驗(yàn)光譜中所有的晶體特征， 證實(shí)了DFT理論在光譜指紋特性仿真計(jì)算中具有很高的計(jì)算效率和精度。 Zhang等[9]通過THz-TDS對尿囊素晶體的THz光譜吸收和色散特性(弱相互作用分析)進(jìn)行解碼， 結(jié)果表明了DFT理論可以仿真晶體的THz光譜指紋特征。

綜上所述， 以L-酒石酸(L-TA)為研究對象， 建立并分析了基于DFT理論的分子模型及THz特征頻譜， 從理論計(jì)算角度驗(yàn)證了實(shí)驗(yàn)吸收峰來源判斷的正確性。 在此基礎(chǔ)上， 分別采用CARS-PLS， MC-UVE-PLS， iPLS和BOSS共4種方法對L-TA在0.2～1.6 THz頻段進(jìn)行特征譜區(qū)的篩選， 建立了定量分析模型。

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

1.1 裝置

采用EKSPLA公司的T-SPEC太赫茲時(shí)域光譜設(shè)備。 裝置采用低溫砷化鎵(LT-GaAs)作為光導(dǎo)天線， FF50飛秒激光器作為超短脈沖激光光源， 中心波長為1 064 nm， 輸出脈寬為150 fs， 發(fā)射器和探測器之間的光程約為62.5 cm。 光路通過光學(xué)透鏡控制和校準(zhǔn)， 系統(tǒng)原理如圖1所示。 飛秒脈沖通過半波片， 由分光鏡將其分為強(qiáng)脈沖和弱脈沖。 強(qiáng)脈沖經(jīng)過斬波器后， 由反射鏡引導(dǎo)經(jīng)可變延遲線入射到LT-GaAs光導(dǎo)天線上， 產(chǎn)生THz電磁輻射脈沖， 然后將該THz脈沖聚焦在待測試的樣品上。 弱脈沖作為探測光， 待測樣品發(fā)射的太赫茲脈沖與微弱脈沖進(jìn)行合并， 合并后的信號(hào)送入鎖相放大器放大處理， 得到含有待測樣品信息的THz時(shí)域頻譜。 本實(shí)驗(yàn)在室溫下進(jìn)行， 整個(gè)光譜處于充入氮?dú)獾姆忾]箱， 相對濕度控制在5%以下， 以減少空氣中水蒸氣對THz波的吸收影響。

圖1 T-SPEC太赫茲時(shí)域光譜設(shè)備原理圖

1.2 樣品制備與數(shù)據(jù)獲取

實(shí)驗(yàn)樣品L-TA和聚乙烯固體粉末均購于Signa-aldrich公司， 樣品純度大于98%， 使用時(shí)未進(jìn)一步提純處理。 將L-TA和聚乙烯混合， 得到L-TA濃度分別為10%， 20%， 40%， 50%， 60%和80%的混合粉末， 將其放入研缽中研磨以減少顆粒引起的散射效應(yīng)， 保證兩者混合均勻。 分別取各濃度適量粉末， 置于15 MPa壓力下持續(xù)壓制2 min， 壓制成約直徑為13 mm， 厚度1.0 mm表面光滑無裂痕、 上下面平行的圓形片狀樣本。 將樣片置于T-SPEC太赫茲時(shí)域光譜設(shè)備中， 通過改變固定樣片的二位平移臺(tái)和THz波匯聚點(diǎn)的相對位置， 每隔0.2 mm采集一個(gè)樣本點(diǎn)的光譜數(shù)據(jù)。

1.3 光譜數(shù)據(jù)處理

1.3.1 光學(xué)參數(shù)

THz波的原始信號(hào)是時(shí)域信號(hào)， 經(jīng)快速傅里葉變換(FFT)后， 將參考信號(hào)r(t)、 樣品信號(hào)s(t)、 樣品厚度d和光速c作為輸入， 得到透射函數(shù)H(ω)和吸收系數(shù)αs(ω)[10]表達(dá)式如式(1)和式(2)

(1)

(2)

式(1)和式(2)中，Ks(ω)為消光系數(shù)，ω為角頻率，n0和ns分別為真空和被測樣品復(fù)折射率。

1.3.2 原始光譜及光譜預(yù)處理

剔除明顯異常數(shù)據(jù)， 每個(gè)濃度獲得57組， 共計(jì)342組有效數(shù)據(jù)。 由FFT得到吸收譜， 如圖2(a)所示。 實(shí)驗(yàn)光譜探測區(qū)間為0～2.3 THz， 根據(jù)信噪比選擇0.2～1.6 THz區(qū)間進(jìn)行分析。 經(jīng)多種預(yù)處理方法的比較， 采用Savtizky-Golay平滑(SG)及均值中心化(mean centering)對數(shù)據(jù)預(yù)處理。 處理后吸收譜如圖2(b)所示， 1.1 THz附近出現(xiàn)明顯的特征波峰， 這與論文資料一致[11]。 各濃度光譜曲線有輕微的偏移， 且相互疊加重合， 掩蓋了光譜的指紋特性。 實(shí)驗(yàn)所得THz光譜特征波峰出現(xiàn)的偏移可能是因?yàn)閷?shí)驗(yàn)過程中環(huán)境溫度的波動(dòng)所致。

1.3.3 樣本選擇方法

為了有效地覆蓋多維向量空間， 增加樣本間的差異性和代表性， 提高模型穩(wěn)定性， 采用SPXY方法[12]將樣品分為訓(xùn)練集和測試集， 如表1所示。

表1 L-TA樣本劃分?jǐn)?shù)據(jù)統(tǒng)計(jì)表

2 理論計(jì)算分析

2.1 DFT仿真計(jì)算

使用GAUSSIAN 16W在Windows 10系統(tǒng)環(huán)境中進(jìn)行仿真計(jì)算。 由劍橋晶體數(shù)據(jù)中心獲得單分子構(gòu)型[13]如圖3所示， 采用DFT理論B3LYP方法， 對氧原子和氫原子分別添加d和p軌道函數(shù)， 以6-31G*(d， p)為基組對單分子結(jié)構(gòu)進(jìn)行優(yōu)化迭代計(jì)算[14-16]。 相關(guān)研究證明， 理論計(jì)算的頻率是諧振頻率， 由于忽視了非諧振效應(yīng)， 基于此計(jì)算基頻等可能會(huì)有一定誤差。 因此添加矯正因子(0.960)進(jìn)行頻率矯正。 經(jīng)優(yōu)化結(jié)構(gòu)的迭代計(jì)算與頻率的計(jì)算， 最終得到的L-TA計(jì)算結(jié)果沒有虛頻， 說明得到了分子能量最穩(wěn)定的結(jié)構(gòu)。

圖3 實(shí)驗(yàn)及仿真計(jì)算在1.11 THz分子振動(dòng)示意圖

在0.2～1.6 THz對L-TA單分子進(jìn)行仿真計(jì)算。 隨機(jī)選取一組實(shí)驗(yàn)數(shù)據(jù)進(jìn)行擬合。 由圖3可以看出， 圖中黑色曲線表示高斯仿真計(jì)算光譜， 分子中C1相連羥基的轉(zhuǎn)動(dòng)帶動(dòng)碳鏈C1—C2—C3—C4進(jìn)行集體振動(dòng)， 同時(shí)C2—H4， C3—H5和O13—H14形成的羥基有輕微擺動(dòng)， 造成了L-TA在1.1 THz處特征峰的出現(xiàn)。 圖3中紅色曲線表示L-TA實(shí)驗(yàn)光譜擬合曲線， 與仿真計(jì)算所得光譜的比較發(fā)現(xiàn)， 實(shí)測光譜吸收峰較仿真計(jì)算光譜吸收峰向前偏移0.02 THz， 其輕微的譜峰偏移可能是由于實(shí)驗(yàn)測量是在室溫下進(jìn)行而理論仿真是基于絕對零度的計(jì)算下得出的。

2.2 基于BOSS算法的L-TA光譜的特征譜區(qū)篩選

2.2.1 變量的選取

在變量空間中利用BBS生成K個(gè)子集， 對每個(gè)子集提取BBS選擇的變量， 并剔除重復(fù)變量以保持唯一性， 對剔除后的變量賦予相同的權(quán)重(w)。

對提取的K個(gè)子集建立PLS模型， 以交叉驗(yàn)證均方根誤差(RMSECV)為模型判斷標(biāo)準(zhǔn)， 選取RMSECV最小值的模型為最佳模型， 采用決定系數(shù)(R2)對模型進(jìn)行評價(jià)。

(3)

(4)

將回歸向量中的元素變?yōu)榻^對值并歸一化處理， 再通過回歸向量累加值結(jié)果賦予新的權(quán)重。

(5)

式(5)中，K是子模型的數(shù)量，bi, k是第k個(gè)子模型中變量i的歸一化回歸系數(shù)絕對值。

根據(jù)變量新的權(quán)重選取新的子集， 保證最佳子模型中回歸系數(shù)較大的變量有更大的概率進(jìn)入最優(yōu)模型。 重復(fù)以上步驟， 直至新子集變量為1。

3 結(jié)果與討論

3.1 BOSS算法對L-TA分析

采用BOSS算法對L-TA的THz光譜數(shù)據(jù)進(jìn)行分析， 由圖4(a)可以看出， RMSECV值在1～7次迭代中穩(wěn)定降低， 由0.030 0減小為0.026 7， 在第7～12次迭代有持續(xù)上升趨勢， 由0.0267升為0.028 7， 隨后增長趨勢明顯， 在第7次迭代時(shí)有唯一最小的RMSECV值0.026 7。 因此， 選取迭代數(shù)為7時(shí)的變量集進(jìn)行分析。

圖4 RMSECV值隨迭代次數(shù)變化圖(a)和迭代次數(shù)為7時(shí)的權(quán)重圖(b)

提取第7次迭代中變量集賦予權(quán)重的值， 如圖4(b)所示， 在0.4～0.6 THz有一個(gè)權(quán)重波峰， 但賦予權(quán)重較小。 在0.98～1.36 THz頻段， 有非常明顯的權(quán)重波峰， 遠(yuǎn)大于0.4～0.6 THz的權(quán)重值。 因此， 根據(jù)權(quán)重值， 將0.98～1.36 THz選為特征譜區(qū)。

3.2 不同建模方法比較

為了進(jìn)一步驗(yàn)證BOSS算法對L-TA光譜數(shù)據(jù)特征譜區(qū)的篩選效果， 將該算法與iPLS， CARS-PLS和MC-UVE-PLS三種經(jīng)典的特征譜區(qū)選擇方法及全譜PLS模型進(jìn)行比較。 BOSS， CARS-PLS和MC-UVE-PLS模型參數(shù)均設(shè)置為5折交叉驗(yàn)證， 采樣次數(shù)為500次， 最大潛在變量的數(shù)量為10個(gè)， 迭代次數(shù)為50次， 建模前均中心化預(yù)處理。 對iPLS模型參數(shù)設(shè)置為5折交叉驗(yàn)證， 建模前中心化預(yù)處理， 主成分?jǐn)?shù)量(principal components, PCs)經(jīng)PCA回歸確定為5個(gè)， 子區(qū)間數(shù)范圍為2～60。

通過不同選擇方法對L-TA的THz光譜進(jìn)行特征譜段篩選。 iPLS特征譜區(qū)篩選效果受子區(qū)間數(shù)量設(shè)置的影響較大， 子區(qū)間數(shù)為4時(shí)有最小的RMSECV和最大的R2。 由圖5可以看出， iPLS選擇譜區(qū)較為集中， 但誤差選擇范圍也大于其他算法， 結(jié)果在表2得到驗(yàn)證， 其RMSECV和R2僅高于全譜PLS模型， 低于其他三種算法。 BOSS， CARS-PLS和MC-UVE都選擇了0.45， 0.93和1.1～1.24 THz區(qū)域， 這對應(yīng)了DFT模型計(jì)算中L-TA分子O—H鍵轉(zhuǎn)動(dòng)的影響。 CARS-PLS和MC-UVE-PLS傾向于選擇更多的變量， CARS-PLS在0.33和1.5 THz附近選擇了譜區(qū)， MC-UVE-PLS在0.2～0.4和0.92 THz也有譜區(qū)的選擇， 但這些譜區(qū)未體現(xiàn)指紋特征， 為無效信息譜區(qū)。 相較于CARS-PLS和MC-UVE-PLS， BOSS的特征譜區(qū)選擇更為集中， 除0.5 THz附近有篩選的譜區(qū)， 其他篩選譜區(qū)均在0.93～1.34 THz頻段， 與DFT仿真計(jì)算光譜特征譜區(qū)有良好的重合度， 較完整地提取了特征波峰所在譜區(qū)， 對特征譜區(qū)體現(xiàn)更加充分。

圖5 多方法對L-TA THz光譜的譜區(qū)篩選

結(jié)果表明， BOSS算法相較于其他三種經(jīng)典譜區(qū)算法， 與DFT仿真光譜特征譜區(qū)有良好的一致性， 譜區(qū)篩選效果更好， 從而使模型具有更好的預(yù)測性能和穩(wěn)定性能， 這一研究也為其他復(fù)雜成分實(shí)驗(yàn)樣品的定量分析提供了參考。

4 結(jié) 論