范贤光，刘 龙，支瑜亮，康哲铭，夏 宏，张佳杰，王 昕*

1. 厦门大学航空航天学院仪器与电气系，福建 厦门 361005 2. 传感技术福建省高等学校重点实验室，福建 厦门 361005

引 言

拉曼光谱作为一种无损快速检测技术，在材料、 化工、 生物医学、 环保等领域已有着广泛的应用[1-3]。然而，大多数生物、 药材等样本中的自发拉曼信号都较弱且易受荧光背景和噪声的干扰，因此采谱时间往往较长，导致常规拉曼成像的动态性能较差，严重阻碍了拉曼光谱成像技术在快速动态体系中的应用。

目前，常用的拉曼成像技术有两大类[4-5]：自发拉曼成像和非线性拉曼成像。自发拉曼成像主要包括扫描模式和宽场模式。扫描模式可获得完整的三个维度信息，但是扫描时间过长，难以观察快速变化的动态过程；宽场模式可获取更好的动态性能，但是只能收集两个空间维度的单波长光谱。非线性拉曼成像主要包括表面增强拉曼散射成像和针尖增强拉曼光谱成像。这两种成像技术具备高空间分辨率和快速成像的优势，但是其对样品有一定的损伤且仪器结构复杂、 造价昂贵，一般由实验室自行搭建，用于前沿科学研究，商品化仪器相对较少，推广难度大。

因此，开发低成本且快速成像的技术至关重要。多通道拉曼成像技术[6]是解决这一问题的有效途径之一，其基本原理如图1所示。硬件上，通过宽场模式激发样品，将散射光同时经过多个不同参数的滤光通道，耦合于面阵CCD的不同感光区域，从而得到多通道窄带数据；软件算法上，基于事先标定好的转移矩阵(不同样本的转移矩阵不同)将多通道数据重建为完整的拉曼光谱，从而实现拉曼成像。

由图1可知，多通道拉曼光谱成像的核心是拉曼光谱标定-重建算法。目前，常用光谱重建算法有伪逆法[7-8]、 Wiener估计算法[9]等。现有的这些算法来源于光谱反射率重建领域，并非专门的多通道拉曼重建算法，要么没有考虑非线性因素的影响，要么没有对训练样本进行优化处理且计算量过大。所以，本文提出一种专用的多通道拉曼光谱重建算法，首先对训练样本进行优化处理，排除坏样本的干扰；然后对得到的多通道窄带测量值进行多项式回归拓展，降低非线性因素的影响；最后引入归一化处理，并依据主成分分析[10-11](principal component analysis，PCA)简化计算，完成拉曼光谱快速重建。同时，在仿真与实验中对伪逆法、 Wiener估计算法等进行了验证，并与本文算法进行了比较。

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

1 算 法

1.1 光谱重建问题阐述

假设待测样本的拉曼光谱为向量r1(N×l维)，N取决于拉曼光路的分光光栅参数和线阵CCD的像素点个数。通过滤光片组产生相应窄带测量值由u1表示，二者之间的关系如式(1)

u1=Tr1+n

(1)

其中，n表示(M×1维)窄带测量中的随机噪声，T(M×N维，其中M为滤光片数量)表示滤光片组的光谱响应函数矩阵，u1是M×1向量。

快速重建拉曼光谱过程就是建立r1和u1之间的逆向数学模型，即转移矩阵W(N×M维)。通过多通道拉曼成像系统，获得实测窄带测量值u，进而重建完整拉曼光谱r，如式(2)

r=Wu

(2)

1.2 基于PCA的拉曼光谱重建

针对图1中测试样品A，在开始建立重建算法前，借助常规拉曼光谱仪完成对样品A的多次测量，获取训练样本集R={r1,r2, …,ri, …,rk}(N×k维，其中k为样本数量)。根据式(1)，不计随机噪声，完成窄带通滤光片从整个拉曼光谱中的特征峰波段内选择性地检测拉曼信号，得样本模拟窄带测量值集U={u1,u2, …,ui, …,uk}(M×k维)。至此，求解转移矩阵W的准备工作已完成，具体求解如下。

1.2.1 预处理阶段

依据Wiener估计思想，待测样品A的训练样本R数目越多，重建完整拉曼光谱的精度会相应提高。但是，训练样本增多，不可避免的会有差异过大的样本存在，故将样本R进行优化处理，具体操作如下：

(3)

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

(4)

1.2.2 校准阶段

(5)

其中

(6)

式(6)中，上标“T”表示矩阵的转置，上标“-1”表示方阵的逆矩阵。

联立式(3)，式(5)和式(6)，得

(7)

1.2.3 测试阶段

根据校准阶段求得转移矩阵W，联立式(2)，完整由窄带测量值快速重建完整拉曼光谱，整体流程如图2所示。

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

2 实验部分

本文以有机玻璃(polymethyl methacrylate, PMMA)作为实验样品，借助常规拉曼光谱仪(532 nm激光光源)，完成对样品的多次测量，获的取训练样本集R={r1,r2, …,ri, …,rk}，其中k代表测量次数。原始拉曼光谱带有荧光背景，为训练样本优化做准备，对其作荧光去除。本文采用B样条去除荧光背景，处理结果如图3所示。

图3 PMMA原始和去荧光背景拉曼光谱Fig.3 Original spectra and eliminating background fluorescence of PMMA

图4 四款滤光片光谱响应函数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 三种算法重建光谱与训练样本的RMSETable 1 RMSE comparison between reconstructed spectra and training samples by three methods

现为验证本文提出的拉曼光谱重建算法，基于实验室搭建的多通道拉曼窄带图像采集系统，获得PMMA的四组滤光片对应的拉曼窄带图像，如图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固体粉末颗粒，颗粒密度分布不均，所以得到的窄带图像呈现出凹凸感。

根据图6，求得任一像素点对应的四通道窄带测量值u。联立式(2)和式(7)，基于训练样本优化的PCA算法，即可实现拉曼光谱的快速重建。同时，依据伪逆法，Wiener估计算法，也完成拉曼光谱重建，如图7所示。

图7 三种方法对PMMA拉曼光谱的实验重建效果对比Fig.7 Comparison among spectra reconstruction of PMMA based on three algorithms in experiment

首先，定性分析三条重建光谱：三种算法均完成了PMMA拉曼光谱的重建，但是伪逆法和Wiener估计算法在完成光谱重建时，重建光谱在波长568 nm的峰值出现偏离，在波长587 nm的峰值拔高。而本文提出的样本优化的PCA重建算法，不会出现上述问题。

然后，定量分析三条重建光谱：基于式(8)，分别求解出三条重建光谱的RMSE，如表2所示。由表2可知，训练样本优化的PCA算法的RMSE最小。对比表1可知，实验重建拉曼光谱的RMSE，比理论RMSE偏大，其原因应为实测多通道窄带数据中含有杂散光和噪声的干扰。

表2 实验重建光谱与训练样本的RMSETable 2 RMSE between reconstructed spectra and training samples in experiment

最后，根据本文重建算法，完成对图6多通道拉曼窄带图像某一指定区域重建拉曼图像，具体如图8(a)所示。该区域包含131×131(1.43×1.43 mm2)个像素点。重建该区域拉曼图像的步骤，具体如下：

(1) 提取该区域内131×131个像素点的对应的多通道窄带响应值；

(2) 根据本文提出的重建算法，联立式(2)和式(7)，得到131×131条完整的重建拉曼光谱；

(3) 对上述重建光谱在波长560/10 nm内，执行积分，得到一个131×131的矩阵，即为重建后的拉曼图像，如图8(d)所示。

对比验证本文算法重建拉曼图像的质量，依据伪逆法和Wiener估计算法重建拉曼图像，如图8(b)和8(c)所示。实验中，图8(a)的中心位置PMMA的浓度仅比周围略高，图8(b)和(c)在该位置明显过亮，其原因在于激发激光在该处的干扰。由图8分析可知，基于样本优化的PCA重建拉曼图像图8(d)中间亮度相对较弱，在一定程度上克服了激光干扰，更好地展现了PMMA的拉曼图像。同时，本文提出的重建算法，实现131×131个像素点重建拉曼图像所需的时间远低于常规的点扫描和线扫描成像模式(如表3所示)。

表3 不同拉曼成像形式所需成像时间Table 3 The time of Raman imaging among different imaging formats

3 结 论

基于样本优化和主成分分析，提出了一种快速拉曼光谱重建方法。该方法首先通过优化训练样本，剔除差样本干扰；然后，引入多项式回归，降低非线性因素影响；最后利用PCA，经归一化处理，完成拉曼光谱重建。仿真和实验结果表明，本文算法比传统算法更好地完成了光谱重建，不会出现重建光谱波峰的偏移和突变。同时，借助自行开发的多通道拉曼成像仪器，利用本文算法完成了PMMA的拉曼光谱重建和成像。因此，本文的研究成果，有力地完善了多通道拉曼成像技术，进而推进了多通道拉曼成像技术在快速动态体系中的应用。