基于樣本優(yōu)化和主成分分析的多通道拉曼光譜重建及其快速成像

范賢光，劉 龍，支瑜亮，康哲銘，夏 宏，張佳杰，王 昕*

1. 廈門大學(xué)航空航天學(xué)院儀器與電氣系，福建 廈門 361005 2. 傳感技術(shù)福建省高等學(xué)校重點(diǎn)實(shí)驗(yàn)室，福建 廈門 361005

引 言

拉曼光譜作為一種無損快速檢測技術(shù)，在材料、 化工、 生物醫(yī)學(xué)、 環(huán)保等領(lǐng)域已有著廣泛的應(yīng)用[1-3]。然而，大多數(shù)生物、 藥材等樣本中的自發(fā)拉曼信號(hào)都較弱且易受熒光背景和噪聲的干擾，因此采譜時(shí)間往往較長，導(dǎo)致常規(guī)拉曼成像的動(dòng)態(tài)性能較差，嚴(yán)重阻礙了拉曼光譜成像技術(shù)在快速動(dòng)態(tài)體系中的應(yīng)用。

目前，常用的拉曼成像技術(shù)有兩大類[4-5]：自發(fā)拉曼成像和非線性拉曼成像。自發(fā)拉曼成像主要包括掃描模式和寬場模式。掃描模式可獲得完整的三個(gè)維度信息，但是掃描時(shí)間過長，難以觀察快速變化的動(dòng)態(tài)過程；寬場模式可獲取更好的動(dòng)態(tài)性能，但是只能收集兩個(gè)空間維度的單波長光譜。非線性拉曼成像主要包括表面增強(qiáng)拉曼散射成像和針尖增強(qiáng)拉曼光譜成像。這兩種成像技術(shù)具備高空間分辨率和快速成像的優(yōu)勢，但是其對樣品有一定的損傷且儀器結(jié)構(gòu)復(fù)雜、 造價(jià)昂貴，一般由實(shí)驗(yàn)室自行搭建，用于前沿科學(xué)研究，商品化儀器相對較少，推廣難度大。

因此，開發(fā)低成本且快速成像的技術(shù)至關(guān)重要。多通道拉曼成像技術(shù)[6]是解決這一問題的有效途徑之一，其基本原理如圖1所示。硬件上，通過寬場模式激發(fā)樣品，將散射光同時(shí)經(jīng)過多個(gè)不同參數(shù)的濾光通道，耦合于面陣CCD的不同感光區(qū)域，從而得到多通道窄帶數(shù)據(jù)；軟件算法上，基于事先標(biāo)定好的轉(zhuǎn)移矩陣(不同樣本的轉(zhuǎn)移矩陣不同)將多通道數(shù)據(jù)重建為完整的拉曼光譜，從而實(shí)現(xiàn)拉曼成像。

由圖1可知，多通道拉曼光譜成像的核心是拉曼光譜標(biāo)定-重建算法。目前，常用光譜重建算法有偽逆法[7-8]、 Wiener估計(jì)算法[9]等?，F(xiàn)有的這些算法來源于光譜反射率重建領(lǐng)域，并非專門的多通道拉曼重建算法，要么沒有考慮非線性因素的影響，要么沒有對訓(xùn)練樣本進(jìn)行優(yōu)化處理且計(jì)算量過大。所以，本文提出一種專用的多通道拉曼光譜重建算法，首先對訓(xùn)練樣本進(jìn)行優(yōu)化處理，排除壞樣本的干擾；然后對得到的多通道窄帶測量值進(jìn)行多項(xiàng)式回歸拓展，降低非線性因素的影響；最后引入歸一化處理，并依據(jù)主成分分析[10-11](principal component analysis，PCA)簡化計(jì)算，完成拉曼光譜快速重建。同時(shí)，在仿真與實(shí)驗(yàn)中對偽逆法、 Wiener估計(jì)算法等進(jìn)行了驗(yàn)證，并與本文算法進(jìn)行了比較。

圖1 多通道拉曼成像原理圖Fig.1 The theory of multi-channel Raman spectroscopy imaging

1 算 法

1.1 光譜重建問題闡述

假設(shè)待測樣本的拉曼光譜為向量r1(N×l維)，N取決于拉曼光路的分光光柵參數(shù)和線陣CCD的像素點(diǎn)個(gè)數(shù)。通過濾光片組產(chǎn)生相應(yīng)窄帶測量值由u1表示，二者之間的關(guān)系如式(1)

u1=Tr1+n

(1)

其中，n表示(M×1維)窄帶測量中的隨機(jī)噪聲，T(M×N維，其中M為濾光片數(shù)量)表示濾光片組的光譜響應(yīng)函數(shù)矩陣，u1是M×1向量。

快速重建拉曼光譜過程就是建立r1和u1之間的逆向數(shù)學(xué)模型，即轉(zhuǎn)移矩陣W(N×M維)。通過多通道拉曼成像系統(tǒng)，獲得實(shí)測窄帶測量值u，進(jìn)而重建完整拉曼光譜r，如式(2)

r=Wu

(2)

1.2 基于PCA的拉曼光譜重建

針對圖1中測試樣品A，在開始建立重建算法前，借助常規(guī)拉曼光譜儀完成對樣品A的多次測量，獲取訓(xùn)練樣本集R={r1,r2, …,ri, …,rk}(N×k維，其中k為樣本數(shù)量)。根據(jù)式(1)，不計(jì)隨機(jī)噪聲，完成窄帶通濾光片從整個(gè)拉曼光譜中的特征峰波段內(nèi)選擇性地檢測拉曼信號(hào)，得樣本模擬窄帶測量值集U={u1,u2, …,ui, …,uk}(M×k維)。至此，求解轉(zhuǎn)移矩陣W的準(zhǔn)備工作已完成，具體求解如下。

1.2.1 預(yù)處理階段

依據(jù)Wiener估計(jì)思想，待測樣品A的訓(xùn)練樣本R數(shù)目越多，重建完整拉曼光譜的精度會(huì)相應(yīng)提高。但是，訓(xùn)練樣本增多，不可避免的會(huì)有差異過大的樣本存在，故將樣本R進(jìn)行優(yōu)化處理，具體操作如下：

(3)

ui1ui4,ui2ui2, …,ui2ui4, …,ui3ui4,ui4ui4}

(4)

1.2.2 校準(zhǔn)階段

(5)

其中

(6)

式(6)中，上標(biāo)“T”表示矩陣的轉(zhuǎn)置，上標(biāo)“-1”表示方陣的逆矩陣。

聯(lián)立式(3)，式(5)和式(6)，得

(7)

1.2.3 測試階段

根據(jù)校準(zhǔn)階段求得轉(zhuǎn)移矩陣W，聯(lián)立式(2)，完整由窄帶測量值快速重建完整拉曼光譜，整體流程如圖2所示。

圖2 基于PCA的拉曼光譜重建流程圖Fig.2 The process of Raman spectrum reconstruction based on PCA

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

本文以有機(jī)玻璃(polymethyl methacrylate, PMMA)作為實(shí)驗(yàn)樣品，借助常規(guī)拉曼光譜儀(532 nm激光光源)，完成對樣品的多次測量，獲的取訓(xùn)練樣本集R={r1,r2, …,ri, …,rk}，其中k代表測量次數(shù)。原始拉曼光譜帶有熒光背景，為訓(xùn)練樣本優(yōu)化做準(zhǔn)備，對其作熒光去除。本文采用B樣條去除熒光背景，處理結(jié)果如圖3所示。

圖3 PMMA原始和去熒光背景拉曼光譜Fig.3 Original spectra and eliminating background fluorescence of PMMA

圖4 四款濾光片光譜響應(yīng)函數(shù)Fig.4 The spectral response functions of four filters

圖5 三種方法對PMMA拉曼光譜的仿真重建效果對比Fig.5 Comparison among spectra reconstruction of PMMA based on three algorithms in theory

(8)

表1 三種算法重建光譜與訓(xùn)練樣本的RMSETable 1 RMSE comparison between reconstructed spectra and training samples by three methods

現(xiàn)為驗(yàn)證本文提出的拉曼光譜重建算法，基于實(shí)驗(yàn)室搭建的多通道拉曼窄帶圖像采集系統(tǒng)，獲得PMMA的四組濾光片對應(yīng)的拉曼窄帶圖像，如圖6所示。

圖6 PMMA多通道窄帶圖像(a)：560/10 nm濾光片測得窄帶圖像；(b)：570/10 nm濾光片測得窄帶圖像；(c)：580/10 nm濾光片測得窄帶圖像；(d)：590/10 nm濾光片測得窄帶圖像Fig.6 The multi-channel narrow-band image of PMMA(a)：Narrow-band image taken by the filter 560/10 nm; (b)：Narrow-band image taken by the filter 570/10 nm; (c)：Narrow-band image taken by the filter 580/10 nm; (d)：Narrow-band image taken by the filter 590/10 nm

由圖6分析可知，鑒于PMMA的最高特征峰主要集中在波長555～565 nm之間(如圖3可知)，所以560/10 nm濾光片測得窄帶圖像最亮，其他順次變暗。由于測試樣本是PMMA固體粉末顆粒，顆粒密度分布不均，所以得到的窄帶圖像呈現(xiàn)出凹凸感。

根據(jù)圖6，求得任一像素點(diǎn)對應(yīng)的四通道窄帶測量值u。聯(lián)立式(2)和式(7)，基于訓(xùn)練樣本優(yōu)化的PCA算法，即可實(shí)現(xiàn)拉曼光譜的快速重建。同時(shí)，依據(jù)偽逆法，Wiener估計(jì)算法，也完成拉曼光譜重建，如圖7所示。

圖7 三種方法對PMMA拉曼光譜的實(shí)驗(yàn)重建效果對比Fig.7 Comparison among spectra reconstruction of PMMA based on three algorithms in experiment

首先，定性分析三條重建光譜：三種算法均完成了PMMA拉曼光譜的重建，但是偽逆法和Wiener估計(jì)算法在完成光譜重建時(shí)，重建光譜在波長568 nm的峰值出現(xiàn)偏離，在波長587 nm的峰值拔高。而本文提出的樣本優(yōu)化的PCA重建算法，不會(huì)出現(xiàn)上述問題。

然后，定量分析三條重建光譜：基于式(8)，分別求解出三條重建光譜的RMSE，如表2所示。由表2可知，訓(xùn)練樣本優(yōu)化的PCA算法的RMSE最小。對比表1可知，實(shí)驗(yàn)重建拉曼光譜的RMSE，比理論RMSE偏大，其原因應(yīng)為實(shí)測多通道窄帶數(shù)據(jù)中含有雜散光和噪聲的干擾。

表2 實(shí)驗(yàn)重建光譜與訓(xùn)練樣本的RMSETable 2 RMSE between reconstructed spectra and training samples in experiment

最后，根據(jù)本文重建算法，完成對圖6多通道拉曼窄帶圖像某一指定區(qū)域重建拉曼圖像，具體如圖8(a)所示。該區(qū)域包含131×131(1.43×1.43 mm2)個(gè)像素點(diǎn)。重建該區(qū)域拉曼圖像的步驟，具體如下：

(1) 提取該區(qū)域內(nèi)131×131個(gè)像素點(diǎn)的對應(yīng)的多通道窄帶響應(yīng)值；

(2) 根據(jù)本文提出的重建算法，聯(lián)立式(2)和式(7)，得到131×131條完整的重建拉曼光譜；

(3) 對上述重建光譜在波長560/10 nm內(nèi)，執(zhí)行積分，得到一個(gè)131×131的矩陣，即為重建后的拉曼圖像，如圖8(d)所示。

對比驗(yàn)證本文算法重建拉曼圖像的質(zhì)量，依據(jù)偽逆法和Wiener估計(jì)算法重建拉曼圖像，如圖8(b)和8(c)所示。實(shí)驗(yàn)中，圖8(a)的中心位置PMMA的濃度僅比周圍略高，圖8(b)和(c)在該位置明顯過亮，其原因在于激發(fā)激光在該處的干擾。由圖8分析可知，基于樣本優(yōu)化的PCA重建拉曼圖像圖8(d)中間亮度相對較弱，在一定程度上克服了激光干擾，更好地展現(xiàn)了PMMA的拉曼圖像。同時(shí)，本文提出的重建算法，實(shí)現(xiàn)131×131個(gè)像素點(diǎn)重建拉曼圖像所需的時(shí)間遠(yuǎn)低于常規(guī)的點(diǎn)掃描和線掃描成像模式(如表3所示)。

表3 不同拉曼成像形式所需成像時(shí)間Table 3 The time of Raman imaging among different imaging formats

3 結(jié) 論

基于樣本優(yōu)化和主成分分析，提出了一種快速拉曼光譜重建方法。該方法首先通過優(yōu)化訓(xùn)練樣本，剔除差樣本干擾；然后，引入多項(xiàng)式回歸，降低非線性因素影響；最后利用PCA，經(jīng)歸一化處理，完成拉曼光譜重建。仿真和實(shí)驗(yàn)結(jié)果表明，本文算法比傳統(tǒng)算法更好地完成了光譜重建，不會(huì)出現(xiàn)重建光譜波峰的偏移和突變。同時(shí)，借助自行開發(fā)的多通道拉曼成像儀器，利用本文算法完成了PMMA的拉曼光譜重建和成像。因此，本文的研究成果，有力地完善了多通道拉曼成像技術(shù)，進(jìn)而推進(jìn)了多通道拉曼成像技術(shù)在快速動(dòng)態(tài)體系中的應(yīng)用。