基于离散极大值法的太赫兹光谱特征吸收峰提取方法

曹雨齐，康旭升1, *，陈飘云，谢 晨，喻 洁*，黄平捷，侯迪波，张光新

1.浙大城市学院计算机与计算科学学院，浙江 杭州 310015 2.浙江大学控制科学与工程学院，工业控制技术国家重点实验室，浙江 杭州 310027

引 言

太赫兹波是频率位于0.1～10 THz(1 THz=1012Hz)的电磁波，其波长范围为0.03～3 mm。由于太赫兹波具有低能性、指纹特性、强穿透性等特点，因此太赫兹时域光谱技术被广泛应用于食品安全检测[1-5]、药品质量检测[6-10]、医学诊断[11-16]等多个应用领域。

太赫兹光谱中的特征吸收峰在含峰目标物的鉴定研究中起着重要作用[17-20]。浙江大学张宏建团队利用四种杀虫剂独特的特征吸收峰强度对四种杀虫剂进行定量检测[17]。三聚氰胺混合物在0.1～3.0 THz频率范围内被观察到三个特征吸收峰，分别位于2.00，2.26和2.60 THz处。该研究根据混合物样品的太赫兹吸收系数谱图中的特征峰信息对三聚氰胺进行检测。结果表明，太赫兹时域光谱技术在食品质量检测中具有潜在的应用价值[18]。

上述研究中，样品均具有明显特征吸收峰，可被肉眼识别，但部分弱峰可能被忽略，从而影响含峰目标物含量的精准检测。因此研究一种由计算机自主有效地提取太赫兹光谱中的吸收峰位特征信息的方法，是太赫兹波广泛应用于各领域物质识别的必要性研究。

为实现药品的太赫兹弱吸收峰的精准识别，本文提出了伴随拐点法，并将其应用于食品添加剂硝基呋喃类化合物的太赫兹吸收峰识别中。研究结果表明伴随拐点法对弱吸收峰具有较好的识别作用。

1 实验部分

1.1 实验系统与样品制备

实验采用Z3太赫兹时域光谱系统(Zomega Terahertz Corporation, USA)，系统的详细介绍可参考文献[21]。实验测试的各硝基呋喃类药物的纯物质均为结晶性粉末。为制备样品压片，首先将样品粉末放置于真空干燥柜中干燥，然后在玛瑙研钵中充分研磨。最后，将样品粉末压制成厚度为1 mm、直径为13 mm的样品压片。由于纯硝基呋喃粉末不易压制成型，需与聚乙烯粉末按1∶1的比例进行混合后压片。

1.2 离散极大值法

连续极大值法是指根据连续函数曲线f(x)中x点的一阶和二阶导数f′(x)和f″(x)数值符号来判断该点是否为特征吸收峰位。如果f′(x)=0且f″(x)<0, 则将x点判定为f(x)的吸收峰位。

然而由于太赫兹时域光谱系统采样不连续，因此频域谱图为离散函数，则可能存在含峰物质的峰位未被采样，导致没有x满足f′(x)=0。因此，针对离散数据，判别吸收峰位的条件为对点x，f′(x)的取值需经过从正值到负值的转换，同时f″(x)<0。由于实验系统中Δx为固定的非零常数，f′(x)能够利用中心差分公式近似计算，如式(1)。同理计算得到f″(x)，如式(2)

(1)

(2)

根据离散极大值法识别特征吸收峰的方法如下。设xi为频段范围内任意频率点，分别利用式(1)和式(2)近似计算吸收系数谱图的一阶和二阶导数，若f′(xi-1)f′(xi)<0且f″(xi)<0，则判定xi为特征谱的极大值点或者吸收峰位。

1.3 伴随拐点法

虽然离散极大值法对强吸收峰具有较强识别能力，但它对弱吸收峰的敏感性不强，这是因为离散极大值法仅根据该点或在邻域内确定峰值位置。若吸收峰位于采样间隔Δx内而没有被采样，则一阶导数的离散曲线可能并不形成过零点，导致弱吸收峰被忽略。因此，本文基于离散极大值法提出一种能够同时测定强吸收峰和弱吸收峰的吸收峰识别方法，即伴随拐点法。

伴随拐点法利用特征吸收峰的伴随拐点信息来识别对应的特征吸收峰，其示意图如图1所示，其中a，b，c，d为吸收峰，a1与a2、b1与b2、c1与c2、d1与d2分别为a，b，c，d的伴随拐点对。伴随拐点法的具体步骤如下：(1)获取样品的原始太赫兹吸收系数谱(Origin，如图1蓝色曲线所示，其中部分曲线被红色曲线重叠覆盖未显示)，确定拐点；(2)确立拐点中的各伴随拐点对(如图1中的a1与a2，b1与b2)；(3)连接伴随拐点对，与原始谱共同确定基线谱(Baseline，如图1中红色曲线所示)；(4)将原始谱与基线谱作差，确定原始谱与基线谱的差谱(Difference，如图1中褐色曲线所示)；(5)利用离散极大值法计算差谱的吸收峰位，即为目标检测物的特征吸收峰。

图1 伴随拐点法示意图

伴随拐点法中伴随拐点对的确定方法如下。首先，设xi为频段范围内任意频率点，若满足f″(xi-1)f″(xi)<0，则判定xi为特征谱的拐点。然后，根据f″(xi-1)和f″(xi)的符号将拐点分类，即起势拐点和落势拐点。定义若f″(xi-1)>0且f″(xi)<0，则xi为起势拐点；若f″(xi-1)<0且f″(xi)>0，则xi为落势拐点。最后，将连续出现的一对起势拐点和落势拐点作为一对伴随拐点。

为避免伴随拐点法过敏感，本研究提出差谱中对应峰值强度设定吸收峰的识别规则。具体为：如果利用“离散极大值法”未在差谱中的某伴随区间内识别出极大值点，则直接将对应候选吸收峰排除；否则如果在该区间内识别出极大值点，但在差谱中的峰值强度低于最大吸收峰峰值强度的5%或者在原始吸收谱中对应的吸收系数在20%分位数以下，则认为其吸收强度过低，也对应候选吸收峰排除；其余识别出的极大值点均可认定为原始吸收谱的吸收峰位。

2 结果与讨论

2.1 极大值法与伴随拐点法的比较分析

利用伴随拐点法对一种易被滥用的硝基呋喃类药物呋喃妥因进行了特征吸收峰的识别，结果如图2所示。为验证伴随拐点法的有效性，同样采用离散极大值法对该样品进行吸收峰位识别，并比较分析。黑色曲线为呋喃妥因在THz波段的原始吸收系数谱图，蓝色曲线为原始吸收系数谱图的一阶导数，红色曲线为对应二阶导数。

图2 呋喃妥因的吸收系数曲线及一、二阶导数曲线

从图中观察到，该物质至少存在三个吸收峰，峰位位于0.90，1.25和1.60 THz附近，其中1.25 THz附近为强吸收峰(图2中B点)，而0.90和1.60 THz处的吸收峰较弱(图2中A和C点)。对于强吸收峰B，在1.25 THz附近的每个一阶导数曲线均随频率的增大由正值转为负值，即1.25 THz附近恰为所有的一阶导数曲线的“过零点”。同时1.25 THz附近所有的二阶导数的取值均小于零，符合离散极大值法判别吸收峰的条件。同样，对于弱吸收峰C，其导数特征也符合离散极大值法的评判指标，能够被离散极大值法识别到。但是，对于弱吸收峰A，其峰位0.90 THz附近，虽然二阶导数均小于零，但在0.90 THz附近的一阶导数均大于零，因此，离散极大值法无法识别到该吸收峰。

由此可见，离散极大值法用于判定强吸收峰十分有效，但用于弱吸收峰的判定时，有效性低，所以在使用此法寻找吸收峰位时有可能会错过一些特征吸收峰，并不利于物质的识别和鉴定。而伴随拐点法恰好能满足这些需求。

呋喃妥因的吸收系数谱图及伴随拐点法的峰位识别结果如图3所示，其中蓝色圆点和红色圆点分别为用于构造伴随区间的起势拐点和落势拐点。图3中伴随区间的数量达到7个，对应的起势拐点和落势拐点的数量达到14个。但是根据差分光谱的谱峰识别，最终识别到4个吸收峰，分别为A，B，C和D。值得注意的是，伴随拐点法不仅能识别弱吸收峰A，而且能识别到另一个肉眼无法识别的弱吸收峰D。伴随拐点法计算得到的差谱中有两个明显的峰和两个相对较弱的峰。因此，根据差分光谱，两个强吸收峰分别对应B和C，两个相对较弱的吸收峰分别对应A和D。差分光谱中的峰位代表了利用伴随拐点法识别到的呋喃妥因的吸收峰位。另外，差分光谱中B峰的强度最强，其次是吸收峰C，这与离散极大值法得到的结果一致。因此，伴随拐点法不管是对强吸收峰还是对弱吸收峰的识别，均十分有效。

图3 呋喃妥因的吸收系数谱图及伴随拐点法的峰位识别结果

而对于候选弱吸收峰E，F和G而言，其强度在差谱中不清晰。根据峰位判别规则，利用离散极大值法未在差谱中识别出E和F点，则直接将其排除。同时，由于G峰值在差谱中强度低于吸收峰B的峰值强度的5%，也排除吸收峰G，其余识别出的极大值点(A，B，C，D)均可认定为原始吸收谱的吸收峰位。

2.2 基于密度泛函理论的呋喃妥因吸收峰的理论计算

为验证伴随拐点法对弱吸收峰的识别能力，将吸收峰识别结果与仿真结果进行对比分析。本研究利用密度泛函理论针对呋喃妥因的吸收峰振动模式进行了理论解析。输入呋喃妥因分子的三维结构文件，在Chem3D软件中采用半经验分子轨道PM3方法，获得呋喃妥因分子的初始最小能量结构，然后在Gaussian量子化学软件中采用密度泛函理论中的B3LYP方法，选取6-311G基组水平对此初始分子结构进行了几何优化和频率计算。计算结果没有虚频出现，说明几何优化计算得到了分子的最小能量结果。

图4为呋喃妥因在0.2～1.8 THz范围内的实验光谱与利用密度泛函理论计算的理论光谱对照图，蓝色曲线和红色曲线分别表示呋喃妥因的实验光谱和理论光谱。理论计算结果显示，呋喃妥因在0.2～1.8 THz范围内共有三个吸收峰，分别位于0.89，1.30和1.67 THz处。由于理论光谱中的吸收系数与实验光谱不在一个数量级上，为便于观察，理论光谱中的各纵坐标采用了同时放大为1 000倍后的数据，其中放大结果不影响吸收峰位的指认。另外，理论光谱中1.30 THz处的吸收峰的峰值与其余两个吸收峰的峰值相比，相对较小，无法在图4中观察到，但我们仍在图4中做了标记，方便对比。

图4 呋喃妥因的理论吸收谱与实验吸收谱对照图

呋喃妥因的太赫兹吸收系数光谱中利用伴随拐点法识别出的在1.25 THz处的强吸收峰，和在0.90和1.60 THz处的两个弱吸收峰均有对应的理论计算结果，分别对应于理论光谱中1.30，0.89和1.67 THz处的吸收峰。而其中实验光谱0.90 THz处的弱吸收峰是无法利用离散极大值法识别出来的，理论光谱中0.89 THz的吸收峰进一步验证了利用伴随拐点法识别实验光谱中吸收峰的合理性。不过实验光谱中利用伴随拐点法识别出的在1.48和1.72 THz处的弱吸收峰未在理论光谱中找到对应的吸收峰，这与理论计算模拟软件所采用分子模拟条件有关。

2.3 伴随拐点法对不同样品吸收峰的识别结果分析

为进一步验证伴随拐点法的有效性，利用该方法对呋喃妥因、呋喃西林、呋喃他酮和呋喃唑酮的吸收峰位也进行识别，结果如图5所示。每个样本吸收系数谱图均有多对伴随拐点，其中蓝点是起势拐点，红点表示落势拐点。以呋喃妥因为例，共识别出10对拐点对，但是最终识别出强度相对较大的5个吸收峰，分别为1.25 THz处的强吸收峰，0.90，1.48，1.60和1.72 THz处的四个弱吸收峰。其中0.90，1.48和1.72 THz处的弱吸收峰采用离散极大值法均无法识别。从结果中观察到，结合伴随拐点法的吸收峰位识别规则，能够有效避免过敏感现象，实现含峰目标物吸收峰位的精准识别。

图5 硝基呋喃类药物吸收系数的伴随拐点法结果

3 结 论

提出了一种基于离散极大值法的含峰目标物的峰位检测方法伴随拐点法，并将其应用于识别和提取硝基呋喃类化合物的太赫兹吸收峰位。同时比较了离散极大值法和伴随拐点法提取的样品峰位信息，观察到伴随拐点法能够有效识别离散极大值法无法识别的弱吸收峰，但同时具有较好的避免过敏感的能力。另外，将密度泛函理论计算得到的理论光谱与太赫兹谱图的识别结果进行比较，进一步验证了伴随拐点法对位于0.89 THz的弱吸收峰位识别的有效性。该结果表明，伴随拐点法在光谱谱图解析和物质检测等领域具有潜在的应用前景。