传统寻找信号峰值方法, 一般都通过求导来实现。 但由于噪声的影响, 求导法很容易获得许多假峰。 文献[6, 7]中介绍了利用连续小波变换(CWT)的办法得到小波系数, 再通过脊线搜索的办法确定峰位置的算法。 但实际使用的效果不佳, 甚至对纯净的信号进行处理时, 峰位置也会存在偏差。

本文中使用阈值循环迭代的办法交替寻找信号中的极小、 极大值。 在寻找极小值点时, 比较当前点值与该点之前到上一次寻找到的极值点之间所有点的最小值之差。 如果该值大于设定的阈值即判断该点的前一个序列点为一个局部极小值点, 可用式(1)来表示

x(k+1)-min{x(p), …, x(k)}> Δ(1)

其中, x(i)为信号序列; x(p)为上一个极值点或者起始点, 若该式成立, 即判定x(k)为一个极小值点。 Δ 为预先设定的迭代阈值。

寻找极大值点时, 比较当前点值与该点之前到上一次寻找到的极值点之间所有点的最大值之差。 可用式(2)来表示

max{x(p), …, x(k)}-x(k+1)> Δ(2)

式(2)成立, 则判定x(k)为一个极大值点。 调整合理的迭代阈值可以准确的找出信号中所有感兴趣的极小和极大值点。

对于拉曼信号而言, 所有的极大值点位置即被确认为拉曼特征峰的位置。 每个峰的区间宽度的估计可以用式(3)来表示

Wi=2×(Pi-Vi)(3)

式(3)中, W_i为第i个峰区间的宽度, P_i为第i个峰的位置, V_i为该峰位置之前相邻的一个极小值的位置。

匹配追踪(matching pursuit)是一种适用于信号稀疏逼近的贪婪迭代算法, 广泛应用在图像、 语音处理, 压缩感知等领域^{[8, 9]}。 所谓稀疏逼近, 指在通过在过完备字典内选择尽可能少的原子的线性组合来表示原信号。 其收敛法则有稀疏度约束和误差约束法则。 对于给定的字典D, 含噪信号y的稀疏重构可由式(4)来描述^{[10, 11]}

γ=argminγ‖γ‖0s.t‖y-Dγ‖22≤ε(4)

其中, γ 为信号y在字典D上的稀疏表示; ε 为收敛误差; Dγ 为重构信号。

设给定的字典D=[d₁, d₂, …, d_m]包含m个原子, 当信号y为长度n的向量时, 原子为长度为n的单位列向量。 设γ ∈ R^m为稀疏系数向量, θ 为迭代残余部分。 每次迭代更新后的稀疏系数向量γ 以及相应的残余部分θ 均加右上标来表示。 则用匹配追踪算法重构含噪信号y的具体步骤如下:

(1) 设定初始稀疏系数向量γ ⁰=[0, 0, …, 0]^T∈ R^m, θ ⁰=y, 信号y可由式(5)表示

y=Dγ0+θ0(5)

(2) 计算第p-1次迭代残余部分θ ^p-1与字典D中各原子的内积的最大值, 并将所得系数加入稀疏系数向量。 见式(6)

αip=maxi∈M(θp-1, di)(6)

其中i为所选原子在字典D中的索引。

(3) 用更新的稀疏系数向量计算第P次迭代的稀疏逼近值y^p=Dγ ^p, 计算迭代残余θ ^p, 见式(7)

θp=θp-1-αip×di(7)

(4) 根据稀疏度约束条件p≤ K或误差约束条件(‖ θ ^p ‖22≤ ε ), 判断是否停止迭代。 其中K为设定的最大稀疏度, ε 为设定的最大残差。 若未达到终止条件, 则转到第2步, 继续搜索匹配的原子。

(5) 迭代终止, 获得含噪信号y基于字典D的稀疏逼近: y˙=Dγ 。

本文提出的算法首先通过循环迭代寻峰的办法在平均谱上找出特征峰的位置和区间。 据此, 可以将噪声谱分割为有信号的区间和无信号的区间。 如式(8)所示

Noise=∑iSi+∑jNj(8)

其中, Noise为噪声信号, S_j为有信号区间, N_j为无信号区间。 对于无信号区间, 其值为0; 对于有信号区间, 先用多项式拟合的办法对其进行基线的扣除, 然后使用匹配追踪的办法来寻求真实拉曼信号的稀疏逼近。 重构后的信号可由式(9)表示。

Noise∧=∑is˙i+∑j0j(9)

重构信号的准确性在一定程度上取决于字典D的设计。 字典的来源一般有两个^[11]: 提前设定和基于学习算法生成字典。 本文中字典是根据特征峰区间以及位置信息来设计的。 单个拉曼峰的分布可以近似认为正态分布^[12], 因此本文中采用高斯密度函数来设计原子, 如式(10)所示

d(x)=12πσe-(x-μ)22σ2(10)

其中μ 为峰位置σ 为峰宽因子。

至此, 可以将利用匹配追踪算法重构拉曼光谱信号的具体步骤描述如下:

(1) 利用阈值循环迭代的方法在平均谱上找出特征峰的位置和区间等信息。

(2) 利用特征峰的位置和区间等信息将待处理的噪声谱划分区间。

(3) 对每个有信号区间光谱分别进行多项式基线校正, 将光谱无信号区间均置0。

(4) 利用匹配追踪算法, 寻求每个有信号区间光谱的稀疏逼近解。

(5) 合成重构后的拉曼光谱信号。

本文实验仪器为nanophoton公司生产的第三代显微拉曼成像系统Raman-11, 该仪器具有较好的拉曼成像性能。 本文实验样品为头孢呋辛酯片粉末, 采用线激光作为激发光, 每次扫描曝光时间为0.5 s。

由于曝光时间短, 激发光能量低以及荧光背景等原因, 单个拉曼光谱的噪声很大。 为验证算法的性能, 分别使用了本文提出的算法以及常规算法对实测的头孢呋辛酯片光谱进行了去噪处理, 实验结果如图5和图6所示。

图5

Fig.5

	Figure Option ViewDownloadNew Window
	图5 (a) 头孢的平均拉曼信号; (b) 头孢的带噪拉曼信号; (c) 本方法获得的去噪信号Fig.5 (a) Average Raman signal of cefuroxime axetil; (b) Noisy Raman signal of cefuroxime axetil; (c) Denoised signal of the Noisy signal by algorithm proposed

图6

Fig.6

	Figure Option ViewDownloadNew Window
	图6 头孢的拉曼信号 (a): 小波软阈值法的去噪结果; (b): 小波硬阈值法的去噪结果; (c): S-G滤波法的去噪结果; (d): EMD算法去噪结果Fig.6 Raman signal of cefuroxime axetil (a): Denoised by wavelet using soft threshold; (b): Denoised by wavelet using hard threshold; (c): Denoised by S-G filter; (d): Denoised by EMD

图5(a)所示为头孢呋辛酯的平均拉曼光谱信号。 图5(b)所示为选取的扫描区域内某一点的光谱信号, 拉曼光谱信号几乎均被噪声淹没。 图5(c)所示为使用本文算法重构后的信号。 参照图5(a)所示平均谱信号可以看出, 重构后的信号不论在谱峰位置上还是峰区间上, 均与平均谱信号较为吻合, 且没有基线的影响无须进一步处理。 表明算法较好的提取了高噪声背景中的弱拉曼信号。

图6(a)所示为采用小波软阈值对图5(b)所示的含噪光谱信号去噪后的结果; 图6(b)所示为采用小波硬阈值去噪后的结果; 图6(c)所示为采用S-G滤波去噪后的结果; 图6(d)所示为采用EMD去噪后的结果。 对比图5(a)所示平均谱信号, 可以发现, 以上四种常规算法在处理低信噪比的拉曼光谱信号时均不理想。 最后, 为了衡量计算速度, 在相同的软硬件环境下(软件: Matlab R2014b, 硬件: Intel Core M-5Y10C, 4GB RAM (DDR3L), 给出了以上各种算法处理头孢呋辛酯片拉曼光谱计算时间, 如表3所示。 由于本文算法同时对光谱的基线和噪声进行了处理, 需要相对较长的处理时间, 但仍保持在0.05s以内, 完全可以满足实际用户需求。

表3

Table 3

表3(Table 3)

表3 不同方法的运行时间 Table 3 Processing times in the experiment of different algorithms Algorithm Processing time/s S-G 0.001 Wavelet using soft threshold 0.012 Wavelet using hard threshold 0.010 EMD 0.017 Algorithm proposed 0.034

表3 不同方法的运行时间 Table 3 Processing times in the experiment of different algorithms