反向匹配法是以纯物质的标准拉曼谱图数据库中的谱图为基准, 逐张去对比待测谱图的特征信息, 而不是用待测谱图去检索数据库, 其匹配流程如图1所示^[10]。 首先待测样品和数据库中纯物质的拉曼谱图都需要经过数据预处理, 如滤波、 校正和基线扣除等; 然后通过寻峰算法^[14]对谱图进行谱峰信息提取, 以获得各个谱峰的拉曼位移、 峰强或峰高等信息。 由于谱图中存在大量噪声, 导致获取的谱峰信息中可能含有噪声峰, 而噪声峰的存在会影响反向匹配方法的准确性, 通常以一个固定的峰强或峰高作为滤噪阈值, 将小于阈值的谱峰作为噪声滤除, 从而获得有效的特征峰信息; 紧接着用标谱数据库中第i张谱图的每个特征峰P_s(n)去比对待测谱图是否在对应位置存在特征峰P_(n), 若存在则P_s(n)被标记为命中。 然后将标谱中所有命中的特征峰的峰强或峰高之和与标谱中所有特征峰的峰强或峰高之和求比值, 从而求得该标谱与待测谱图的相似度Ri。 遍历标谱数据库中所有谱图与待测谱图进行对比, 分别计算出相似度; 最后通过设定的阈值(如80%)将相似度大于阈值的物质标谱标记为匹配。

图1

Fig.1

	Figure Option ViewDownloadNew Window
	图1 反向匹配的基本流程图Fig.1 Basic flow diagram of RMM

反向匹配方法以标准谱图数据为基础, 反向搜索待测谱图的谱峰数据, 这在很大程度上可以避免待测谱图中包含的其他物质的谱峰干扰。 但是, 强峰与弱峰的强度可能会有几个数量级的差异, 这将导致弱峰的权重可能被强峰完全掩盖, 容易造成误匹配。 如图2所示, 含乙酸异丁酯样品中误检出苯, 反向匹配计算出的相似度为92%, 强峰(992 cm^-1)与弱峰(607, 1 586和1 606 cm^-1)的强度存在数量级差异导致了弱峰权重被掩盖, 从而造成错误的匹配结果。 如果以匹配的谱峰数量为基数计算相似度, 又会过度增大弱峰的权重, 因为强峰基本不会消失, 而弱峰容易在采集中丢失, 从而容易造成漏匹配的问题。

图2

Fig.2

	Figure Option ViewDownloadNew Window
	图2 含乙酸异丁酯样品中误检出苯Fig.2 Sample of isobutyl acetate mistakenly detected benzene

其次, 由于不同物质的拉曼散射效应差异很大, 可以通过样品处理、 调整谱图采集参数和使用滤波等方法来提高谱图质量, 其中控制采集谱图的积分时间是常用的一种方法^[15]。 但是随着积分时间的增加, 谱图中噪声峰的强度也会随之增大。 反向匹配通常以一个固定的峰强或峰高作为滤噪阈值, 而这并不能精准地将噪声峰滤除, 甚至可能会滤除有效的特征峰。 如图3所示, 图3(a)和(b)分别为硫酸氢钾和二苯基甲烷的拉曼谱图, 受积分时间的影响, 二者的噪声幅度有较大的不同, 若以2 500为固定的滤噪阈值, 从图3(a)可以看出硫酸氢钾仍有大量噪声峰未被过滤, 而图3(b)二苯基甲烷却有部分有效特征峰被过滤, 由此可以看出用传统的反向匹配方法在噪声峰的过滤上存在较大的问题, 会直接影响匹配的准确性。

图3

Fig.3

	Figure Option ViewDownloadNew Window
	图3 不同噪声强度的谱图的滤噪情况 (a): 硫酸氢钾; (b): 二苯基甲烷Fig.3 Noise filtering effect of spectra under different noise intensities (a): Potassium bisulfate; (b): Diphenylmethane

针对传统反向匹配方法中, 存在的强峰和弱峰的权重合理设置问题, 本研究提出了一种基于衰减函数的权重优化方法。 从反向匹配方法存在的问题分析中可以看出, 物质的强拉曼特征峰在采样中不容易消失, 而弱拉曼特征峰可能会因为样品浓度低或仪器灵敏度差等原因导致部分消失, 因此在待测样品的拉曼谱图的比对中如果物质的强特征峰不匹配, 则表示有很大的概率不含有该物质。 相比之下, 如果弱特征峰不存在, 则表示有相对较小的概率不含有该物质。 从这个角度分析, 在计算相似度的时候强特征峰的权重应该比弱特征峰高, 但是如果直接以峰强或峰高作为权重因子则有可能导致强弱峰的权重相差太大, 此时比对结果基本由强特征峰决定, 而弱特征峰的作用基本被忽略, 从而容易导致误匹配, 尤其是在物质的拉曼特征峰个数较少的情况下该表现更为明显。 为此, 需要对强峰进行衰减, 对弱峰进行适当增强, 使得强峰和弱峰的权重在可比的数量级范围内。 因此需要构造一个关于特征峰强度(或峰高)的权重曲线, 使得强峰权重能够得到抑制, 弱峰权重能够适当提高, 并且随着特征峰个数的变化加重或减轻对谱峰权重的影响效果。 为此, 本研究根据对数函数和指数函数的特征, 构造了一种谱峰的权重衰减函数[式(1)], 以实现特征峰的权重随强度和数量的变化而有效增强或衰减。

W=(logSmaxS)λNmatch (1)

式(1)中, W为调整后的特征峰权重, S为特征峰的强度, S_max为特征峰中最大的强度, N_match为待测谱图和标准谱图相匹配的特征峰个数, λ 为在0~1.0之间的可调比例系数。

为验证该函数对特征峰权重的调节效果, 以匹配的特征峰个数N_match分别为2, 5, 10, 20为例, 权重衰减函数曲线如图4所示。 由于特征峰的强度S为大于1的数, 且S_max≥ S, 因此对数函数在以S_max底数和以S为真数值时, 为斜率小于1的增函数, 从而确保了谱峰的权重与谱峰强度(或高度)为衰减关系, 并且衰减后的权重W将被缩放到1内, 最大谱峰的权重始终为1, 起到了归一化的作用。 其次, 从图中单条曲线可以看出, 随着S的增加, 斜率逐渐由大变小, 从而实现了弱峰权重衰减较小, 强峰权重衰减较大的功能。 此外, 从多条曲线的趋势可以得出, 匹配峰个数较少时(N_match较小), 曲线的斜率变化较大, 即弱峰的权重衰减较缓和, 而强峰的权重衰减更严重。 这是因为当物质的匹配峰的个数较少时, 比较容易出现误判, 因此要提高弱峰的权重, 以确保很少的特征峰中, 弱峰的贡献不被忽略。

图4

Fig.4

	Figure Option ViewDownloadNew Window
	图4 不同谱峰命中个数的权重衰减函数曲线Fig.4 The function curves of weight for the different hit number of peaks

经权重衰减函数处理后, 优化后的反向匹配的相似度计算公式如式(2)所示

R%=∑i=1Nmatch(logSmaxSi)λNmatch∑j=1N(logSmaxSj)λN(2)

式(2)中, S_j为标准谱图中各个特征峰的谱峰强度, N为标准谱图的特征峰个数, S_i为标准谱图中与待测谱图相匹配的各个特征峰的谱峰强度, N_match为相匹配的特征峰个数, S_max为特征峰中最大谱峰强度, λ 为可调比例系数。

经过衰减函数优化后的反向匹配对比效果如图5(a)所示, 权重未优化前含3-叔丁基-4-羟基苯甲醚样品误匹配碘乙烷, 具有较高的匹配度(R=91%), 这是因为碘乙烷标谱中被匹配到的两个特征峰(499和1 201 cm^-1)的信号掩盖了未被匹配上的弱峰信号, 从而导致误判。 而图5(b)为衰减函数优化后的结果, 当λ 取值不同时, 反向匹配中相似度的变化情况。 从图可以看出, 经过权重衰减优化后, 错误匹配的碘乙烷的相似度从原来的91%快速降到了40%附近, 而正确匹配的3-叔丁基-4-羟基苯甲醚相似度保持在80%以上。 这主要是因为误匹配的碘乙烷的相似度主要由强峰(499和1 201 cm^-1)决定, 其他占比小的弱峰(950, 1 050和1 432 cm^-1)都未被匹配, 当通过权重衰减函数调整后, 强峰的权重被减弱, 未匹配的弱峰权重增加, 使得弱峰对相似度的影响变大, 所以碘乙烷的相似度变化较大; 而正确匹配的3-叔丁基-4-羟基苯甲醚只有少数几个弱峰(526, 1 074和1 394 cm^-1)没有被匹配, 所以即使弱峰权重被增强, 对相似度的影响也较小。 从图中还可以得出随着λ 值减少, 弱峰权重逐渐增加, 而谱图中未被匹配的都是弱峰, 所以相似度呈现下降趋势。 如图5(c)和(d)所示的含乙酸异丁酯误匹配到苯的情况与图5(a)和(b)相似, 经衰减函数处理后误匹配的苯的相似度快速下降。 实际匹配过程中, k的取值可根据匹配的谱峰个数进行动态调整, 使得衰减函数的效果达到最优。 引入权重衰减函数处理后的误匹配谱图的相似度出现了明显的降低, 而正确匹配的谱图变化较小, 说明该权重衰减函数能够有效实现强弱特征峰的权重平衡, 提高反向匹配方法的识别准确性。

图5

Fig.5

Figure Option ViewDownloadNew Window

图5 衰减函数效果对比 (a): 权重优化前含3-叔丁基-4-羟基苯甲醚样品中误检出碘乙烷; (b): 权重优化前后含3-叔丁基-4-羟基苯甲醚样品误匹配碘乙烷相似度变化; (c): 权重优化前含乙酸异丁酯样品中误检出苯; (d) 权重优化前后含乙酸异丁酯样品误匹配苯相似度变化Fig.5 Comparison of different attenuation function effect (a): Sample of 3-tert-butyl-4-hydroxyanisole mistakenly detected iodoethane before weight optimization; (b): Similiarity change of sample of 3-tert-butyl-4-hydroxyanisole mistakenly detecting iodoethane while weight optimization; (c): Sample of isobutyl acetate mistakenly detected benzene before weight optimization; (d): Similarity changes of sample of isobutyl acetate mistakenly detected enzene while weight optimization

针对反向匹配的特征峰过滤中, 采用固定阈值过滤噪声峰方法容易导致谱图特征峰提取错误的问题, 本文根据噪声基本满足正态概率分布的特征, 提出了动态阈值滤除噪声峰的方法, 以实现噪声峰的自适应过滤^[16]。 由于谱图中的噪声分布是随机的, 所产生的噪声峰的强度会随着积分时间的改变而增强或减弱, 但其统计学分布基本都集中在一个区域内, 且呈现出正态分布的特征。 因此, 可以通过划定噪声峰强度分布的集中区域, 计算出用于过滤的阈值, 实现符合统计学分布的绝大部分噪声峰的过滤。 设计的动态阈值自适应滤除噪声峰方法的原理如图6所示, 首先对含有大量噪声峰的拉曼谱图数据[图6(a)]求一阶导数, 从而确定谱图中所有极大值点的位置[图6(b)]。 极大值点的位置表示可能含有的谱峰, 因此通过极大值点的位置反向搜索并计算出其在谱图中的对应的真实峰高h[图6(c)], 然后利用直方图对所有峰高进行统计, 再用正态分布函数对直方图进行拟合分析[图6(d)]。 由于绝大部分噪声的强度具有正态分布的特点, 并且噪声峰密集, 而特征峰数量相比于噪声峰很少, 因此将直方图正态分布拟合的谱峰密集的一定区域判定为噪声峰集中区, 而区域外的谱峰判定为特征峰。 最后根据划分的区域边界, 确定动态滤噪的阈值, 实现在不同等级噪声强度下的自适应滤噪。

图6

Fig.6

	Figure Option ViewDownloadNew Window
	图6 动态阈值自适应滤噪方法实现原理 (a): 原始拉曼谱图; (b): 一阶导数谱; (c): 基线扣除后的谱图; (d): 对扣除基线后的峰高进行直方图统计和分布拟合Fig.6 Implementation principle of dynamic threshold adaptive noise filtering method (a): Original Raman spectra; (b): Spectrum after first-order derivative; (c): Spectrum after baseline subtraction; (d): Histogram statistics and fitting for peaks height after baseline subtraction

为验证动态阈值滤噪的效果, 对图3的两份谱图数据进行基线扣除^[17]等预处理后, 进行动态自适应滤噪, 将所有谱峰进行直方图统计和正态分布拟合后, 设定拟合区域的95%面积左右为噪声集中区域, 从而计算出图3(a)的动态滤噪阈值为4 985, 图3(b)为494。 由结果可知, 图3(a)的噪声峰很强, 导致固定的峰高阈值未能将较强的噪声峰滤除, 而采用动态峰高阈值的操作能准确将噪声峰滤除从而保留有效的特征峰信息; 相比之下, 图3(b)的噪声强度较弱, 部分特征峰也较弱, 固定阈值时容易将这些弱特征峰错误地滤除, 而采用动态阈值的方法, 滤噪阈值随噪声强度而变化, 自动调整下来, 从而保证噪声过滤的同时能将弱特征峰有效地保留下来。

首先, 将本文提出的改进的反向匹配方法与传统的反向匹配方法以及命中质量指数方法分别应用于实验室制作的混合物表面增强拉曼谱图的测试集(支持信息表S2)与数据库(支持信息表S1)的比对分析中, 验证各方法的识别准确性。 由于待测谱图和标准谱图来源于不同仪器, 因此在使用命中质量指数(HQI)的计算前进行了相对强度校准^[19]。 各方法的相似度计算范围均为0~1.0, 本文将匹配阈值设置为0.8, 在与数据库中所有标准谱图的相似度比对中, 相似度大于等于阈值时为匹配, 小于阈值时为不匹配。 支持信息表S2中记录了命中质量指数(HQI)、 传统的反向匹配方法(RMM)和改进后的反向匹配方法(IRMM)对于常见的酒水饮料中单独添加不同浓度的苯二氮卓类和毒品类药物组成的混合物样品所采集的拉曼谱图在数据库中的匹配情况, 其中部分物质的识别结果如表1所示。 由表1所示部分物质的识别结果可知, HQI在纯物质的识别中具有较好的匹配效果, 而对于混合物的识别, 由于无法排除其他物质的拉曼信号干扰, 因此无法给出可靠的识别结果。 而RMM方法容易出现多匹配, 即匹配结果中包含正确匹配的物质和误匹配的物质。 相比之下, IRMM方法可以很好地识别混合物中的成分, 可以准确地排除传统的RMM方法中误匹配的物质。

表1

Table 1

表1(Table 1)

表1 支持信息表S2中部分物质(SERS)的匹配结果 Table 1 Matching results of some substances (SERS) in supporting information (S2) 被添加的 酒水饮料 添加的 药物 浓度 /ppm 匹配的物质 HQI RMM IRMM 相似度 /% 匹配 结果 相似度 /% 匹配 结果 相似度 /% 匹配 结果 纯水 吗啡 20 吗啡 82 正确 87 正确 87 正确 摇头丸 20 摇头丸 54 错误 90 正确 82 正确 劳拉西泮 20 劳拉西泮 环戊基甲酰芬太尼 63 28 错误 89 83 存在 多匹 89 57 正确 褪黑素 20 褪黑素 5-甲氧基-2-甲基色胺盐酸盐 87 87 存在 多匹 99 87 存在 多匹 99 79 正确 三唑仑 20 三唑仑 氟阿普唑仑 芬太尼类 91 84 31 存在 多匹 99 74 83 存在 多匹 99 74 66 正确 ORGMA 吗啡 20 吗啡 海洛因 54 65 错误 97 87 存在 多匹 97 69 正确 RIO 摇头丸 20 摇头丸 苯甲酰爱康宁 43 29 错误 97 84 存在 多匹 96 78 正确 豹发力 氯硝西泮 50 氯硝西泮 氯胺酮 64 14 错误 96 83 存在 多匹 96 59 正确 红牛 氯胺酮 200 氯胺酮 芬太尼类 73 42 错误 93 83 存在 多匹 83 61 正确 威士忌 阿普唑仑 50 阿普唑仑 3-甲基硫代芬太尼 15 20 错误 86 81 存在 多匹 86 66 正确

表1 支持信息表S2中部分物质(SERS)的匹配结果 Table 1 Matching results of some substances (SERS) in supporting information (S2)

根据表1中的匹配结果存在的准确匹配和多匹配的情况, 将准确识别的结果定义为SCIR, 对存在多匹配的结果定义为MCIR, 对其余结果即无匹配或错误匹配情况定义为NIR^[20], 则对支持信息表S2识别的统计结果如表2所示。 由表2中的汇总结果可知, HQI在混合物的表面增强拉曼谱图识别中准确识别率仅有14.81%, 这是由于混合物的拉曼谱图由多种物质的谱峰信号叠加组成, 在HQI的计算中存在较大的干扰。 而RMM方法以标准谱图的谱峰为基准反向搜索待测谱图, 这在一定程度上避免了其他物质的拉曼信号干扰, 识别准确率达到44.44%, 但多匹配情况较为严重, 达到55.56%。 而IRMM方法的准确率高达96.30%, 多匹配情况基本没有, 识别准确率大幅提升。

表2

Table 2

表2(Table 2)

表2 支持信息表S2识别结果统计 Table 2 The statistics of matching results of S2 Algorithms SCIR MCIR NIR HQI 14.81% 7.41% 77.78% RMM 44.44% 55.56% 0% IRMM 96.30% 0% 3.70%

表2 支持信息表S2识别结果统计 Table 2 The statistics of matching results of S2

为了进一步验证本文提出的改进IRMM算法的识别性能, 进行了大量的常规拉曼和表面增强拉曼样本数据的匹配测试, 将三种识别方法分别应用于收集自不同仪器的待测谱图(支持信息表S3)的比对分析中, 验证各方法对于实际样本的识别效果。 由于支持信息表S3包含的样本数据量很大, 因此实验数据的统计结果中未列出多匹配中错误匹配的物质名称, 只统计待测谱图是否正确匹配物质和是否存在多匹配情况, 实验的具体结果见支持信息表S3, 统计结果如表3所示。

表3

Table 3

表3(Table 3)

表3 支持信息表S3识别结果统计(单位: %) Table 3 The statistics of matching results of S3 (Unit: %) algorithms 常规拉曼n=355 表面增强拉曼n=152 综合准确率n=507 SCIR MCIR NIR SCIR MCIR NIR SCIR MCIR NIR HQI 60.56 8.73 30.70 28.95 9.87 61.18 51.08 9.07 39.84 RMM 4.79 94.37 0.85 44.08 55.92 0.00 16.57 82.84 0.59 IRMM 90.14 7.89 1.97 94.74 5.26 0.00 91.52 7.10 1.38

表3 支持信息表S3识别结果统计(单位: %) Table 3 The statistics of matching results of S3 (Unit: %)

首先, 由表3的统计结果可知, 基于标准谱图数据库的应用中, IRMM方法对常规拉曼和表面增强拉曼谱图的识别准确率(SCIR)分别为90.14%和94.74%, 说明该方法不仅在常规拉曼谱图的识别中表现良好, 而且对表面增强拉曼的识别同样有效, 并且相比于RMM方法和HQI方法的识别准确率均有大幅度的提升。 而HQI方法在常规拉曼(60.56%)和表面增强拉曼(28.95%)谱图的识别中有较大的差异。 传统的RMM方法在常规拉曼谱图的识别中准确识别的概率只有4.79%, 而多匹配的情况非常严重, 占94.37%。 这主要是由于常规拉曼标准谱图数据库中的标准物质数量有上千种, 并且很多物质特征峰比较尖锐, 存在大量的特征峰重叠问题, 使得匹配的准确度急剧下降到4.79%, 而在表面增强拉曼的标准数据库谱图数量较少, 多匹配的情况为55.92%。 相比之下, 本文提出的IRMM方法在常规拉曼和表面增强拉曼中多匹配的问题有较大的改善, 仅为7.89%和5.26%, 并且综合匹配结果的准确性从16.57%提升到了91.52%, 因此匹配的准确度有大幅度的提升。