基于经验模态分解提高太赫兹频率分辨率的方法

卢 敏，李小霞，*，尚丽平，邓 琥，

1.西南科技大学信息工程学院，四川 绵阳 621010 2.特殊环境机器人技术四川省重点实验室，四川 绵阳 621010

基于经验模态分解提高太赫兹频率分辨率的方法

卢 敏1，李小霞1，2*，尚丽平2，邓 琥1，2

1.西南科技大学信息工程学院，四川 绵阳 621010 2.特殊环境机器人技术四川省重点实验室，四川 绵阳 621010

太赫兹频率分辨率是影响物质鉴别的重要因素，但是由于太赫兹时域光谱系统中存在器件反射，使得参考信号和测量信号出现多个反射峰，时域信号长度截断导致频率分辨率很低。为了去除反射峰的影响，提出基于经验模态分解去除时域反射峰，从而提高太赫兹频率分辨率的方法。通过与太赫兹主峰的互相关定位时域反射峰，计算反射峰的上下包络和求平均值，获取本征模函数并代替反射峰，增加时域信号有效长度，提高太赫兹频域分辨率。空气中水蒸汽的太赫兹透射谱实验结果表明，该方法具有自适应去除多个反射峰的能力，对太赫兹时域信号修复效果良好，频率分辨率提高了12倍，而且未丢失有用的吸收谱信息，吸收峰的位置和个数与真实谱一致，很好地保留了太赫兹谱的鉴别能力。

太赫兹时域光谱；经验模态分解；频率分辨率；反射峰去除

引 言

太赫兹(Terahertz，THz)通常是指频率范围在0.1～10 THz的电磁波，处于宏观电子学和微观光子学之间的过渡区域[1]。在食品、药品、毒品、含能材料、生物分子探测等方面具有非常重要的科学价值[2-5]。通过傅里叶变换将太赫兹时域信号变换到频域，计算出介电常数、折射率和吸收系数等物理参数，可以表征和鉴别物质种类和质量，并可得到其他手段无法获取的物理和化学信息[6-7]。

频率分辨率是指能够区分两个最近的频率分量的能力。一般来说，数字信号的物理(真实)分辨率与信号的实际长度成反比，信号采样长度越长，分辨率越高[8]。

或采用现代功率谱估计的方法可提高分辨率，但是该方法会造成频谱的失真[9]。由于太赫兹时域光谱系统中的探测晶体和天线基底材料的反射，使得系统中的太赫兹光会发生多次反射，时域信号出现多个反射峰，特别是样品的镜面反射和复杂的系统使得反射峰出现更早或更多，使得实际可利用的时域信号长度有限。若将THz时域光谱信号截断至第一反射峰前，根据傅里叶变换理论，其频率分辨率会降低，频率分辨率过低将不能反映频谱的精细程度，影响后续对吸收峰的辨识和物质的鉴别。因此太赫兹时域反射峰的去除是很多太赫兹应用的必要技术。

经验模态分解(empirical mode decomposition，EMD)是1998年由Huang提出的自适应时频分析方法[12]。相比于传统的时域相关分析、频域谐波分解、时频域小波变换等信号分解方法，EMD在处理非平稳、非线性数据上具有非常明显的优势，特别适合处理含有多级反射峰的太赫兹时域信号[11]。

利用EMD方法去除反射峰的影响，提高太赫兹频率分辨率。通过修复前后的时域谱和吸收谱的对比分析，本方法修复太赫兹时域谱和提高频率分辨率的效果。

1 EMD基本原理

EMD的核心思想是能将任意复杂信号分解为若干个本征模函数(intrinsic mode function，IMF)之和，各IMF分量包含了原信号不同时间尺度的局部特征信号，与傅里叶变换的谐波函数和小波变换的小波基函数不同，IMF没有固定的表达式，依据信号自身的情况自适应地确定。Huang认为，IMF分量必须满足下面两个条件：(1)在全部时间范围内，局部极值点数和过零点数必须相等，或最多相差一个；(2)在任意时刻，上包络线(局部最大值的包络)和下包络线(局部最小值的包络)的平均值必须为零[12]。下面利用EMD将复杂的原始信号分解成瞬时频率不同的若干个IMF的特性去除THz反射峰。

2 基于EMD去除太赫兹反射峰

以空气中水蒸汽的太赫兹透射谱为例，参考样本为空气中湿度为8.5%时的水蒸汽，测量样本为空气中湿度为49%时的水蒸汽。

2.1 反射峰定位

如图1所示，首先利用峰值查找方法得到主峰谷点对应时间tv0和峰值点对应时间tp0，tv0往前推20个点即主峰起始点对应时间ts0，tp0往后推20个点即主峰终止点对应的时间te0。然后将主峰信号与后面的信号进行互相关，可得到第一反射峰起始点对应时间ts1、终止点对应时间te1，第二反射峰起始点对应时间ts2、终止点对应时间te2，以此类推。

图1 带有反射峰的太赫兹时域光谱

最后利用时间除以采样间隔分别找到对应的数据序号，第一反射峰起始点ns1、终止点ne1，第二反射峰起始点ns2、终止点ne2。ns1～ne1之间的信号为第一反射峰信号，ns2～ne2之间的信号为第二反射峰信号，以此类推。参考信号和测量信号采用相同的流程进行定位。

2.2 反射峰去除

利用EMD方法对各个反射峰进行IMF分解。以参考信号第一反射峰为例(其他反射峰所用方法与之相同)，具体步骤如下：

(1)找出反射峰对应的时域信号x(t)的所有局部最小值点和局部最大值点；

(2)分别对局部最大值点和最小值点通过最小二乘拟合得到上下包络线xmax(t)和xmin(t)

xmin(t)≤x(t)≤xmax(t)，t∈[ts1，te1]

(1)

(3)计算上下包络线的均值

(2)

(4)得到反射信号去掉包络均值的新数据序列h11(t)

h11(t)=x(t)-m11(t)

(3)

(5)检查h11(t)是否满足IMF的条件，通常对h11(t)重复上述处理过程，直到满足IMF的定义要求，则

h1k(t)=h1(k-1)(t)-m1k(t)

(4)

第一个IMF分量为

(5)

令

r1(t)=x(t)-c1(t)

(6)

r1(t)为原始信号，重复上述步骤，则可以得到其他IMF：c2(t)，c3(t)，…，cn(t)，即

cn(t)=rn-1(t)-rn(t)

(7)

式(7)中rn(t)表示原始信号的趋势信息，称为余项。至此，第一反射峰的EMD分解结束，原信号可以表示为

(8)

根据上述步骤对参考信号的第一反射峰进行EMD分解，得到的3个IMF如图2所示。

图2 EMD分解得到的本征模函数

分解出的IMF按频率由高到低依次排列，通常情况下，频率高的IMF是噪声。因此，选取分解出的最后一个IMF作为反射峰修复后的信号，并用它替代原来的反射峰时域信号，获得修复后的太赫兹时域光谱信号。该方法对参考信号和测量信号的多个反射峰都具有自适应修复的能力。

3 太赫兹参数计算

(9)

其中c是真空中光速，ω是角频率，d是样品厚度。吸收系数计算公式为[15]

(10)

根据式(9)和式(10)计算EMD处理前后的吸收系数。

4 结果与讨论

在常温下分别采用商用仪器(Zomega，Z-3)和自建系统采集湿度为8.5%(参考信号)和49%(测量信号)水蒸汽的太赫兹时域信号，采样间隔66.7 fs，采样长度为200 ps。利用上述的方法去除THz时域反射峰，分别计算EMD处理前后的吸收系数。

因为商用仪器采用了在光电导天线上加半球形硅透镜等其他技术，使得反射峰出现的时间较晚，为了验证本方法的正确性，截取商用仪器采集的第一反射峰前的时域数据，计算后当作真实的吸收谱。

图3和图4分别为参考信号和测量信号EMD处理前后的时域谱对比图。为了能更清楚地比较修复的效果，将修复后的时域谱向下移动0.000 5 V。通过比较可以明显看出，参考信号和测量信号的前三个反射峰都被很好地去掉了。

图3 参考信号修复前后时域对比图

图4 测量信号修复前后时域对比图

图5为修复前后太赫兹吸收谱和真实吸收谱的对比图。

图5 修复前后及真实的吸收谱对比图

比较EMD处理前后的吸收谱可以看出，在1.2，1.6，1.8及1.9 THz附近吸收峰的个数存在很大的差异。在1.2，1.6和1.8 THz附近应该有一个或两个吸收峰，但修复前都没有吸收峰。在1.9 THz附近EMD处理后有两个吸收峰，而修复前只有一个吸收峰，表明通过经验模态分解提高了吸收谱的分辨率。

比较真实谱和EMD处理后的吸收谱，其吸收峰位置和个数都基本一致，说明EMD处理未丢失有用的吸收峰信息，很好地保留了频谱的鉴别能力。在频率为0.556，1.165，1.230，1.604及1.868 THz处，修复后的吸收系数与真实谱相比偏小，这是因为修复后的吸收峰分裂为两个峰，而真实谱只有一个峰。

数字信号的频率分辨率计算公式

(11)

5 结 论

提出利用EMD自适应去除太赫兹时域反射峰、提高太赫兹频率分辨率的方法。为了验证方法的正确性，以空气中水蒸汽的太赫兹吸收谱为例，比较了EMD处理前后的时域谱，反射峰去除效果良好；通过EMD处理前后的吸收谱比较，发现增加了吸收峰的个数，频率分辨率大大提高；通过与真实谱比较，其吸收峰的位置和个数都基本一致，说明EMD处理未丢失有用的吸收峰信息，很好地保留了频谱的鉴别能力。

[1] Siegel P H.IEEE Transactions on Microwave Theory and Techniques, 2002, 50(3): 910.

[2] Massaouti M, Daskalaki C, Gorodetsky A, et al.Applied Spectroscopy, 2013, 67(11): 1264.

[3] HE Ting, SHEN Jing-ling(和 梃，沈京玲).Spectroscopy and Spectral Analysis(光谱学与光谱分析)，2013，33(9)：2348.

[4] Huhn, A K, Saenz E, De Maagt P, et al.IEEE Transactions on Terahertz Science and Technology, 2013, 3(5): 649.

[5] Conti Nibali, Valeria, Havenith Martina.Journal of the American Chemical Society, 2014, 136(37): 12800.

[6] Brahm A, Tymoshchuk M, Wichmann F, et al.Journal of Infrared Millimeter and Terahertz Waves, 2014, 35(11): 974.

[7] LI Li-long，XIANG Yang, WU Lei, et al(李利龙，向 洋，吴 磊，等).High Power Laser and Particle Beams(强激光与粒子束), 2013, 25(6): 1566.

[8] CHENG Pei-qing(程佩青).Digital Signal Processing Tutorial(数字信号处理教程).Beijing: Tsinghua University Press(北京：清华大学出版社)，2007.100.

[9] ZHANG Xian-da(张贤达).Modern Signal Processing(现代信号处理).Beijing: Tsinghua University Press(北京：清华大学出版社)，2002.121.

[10] Wu Z, Huang N E.Advances in Adaptive Data Analysis, 2009, 1(1): 1.

[11] ZHANG Yu-yan, ZHOU Hang, YAN Mei-su(张玉燕，周 航，闫美素).Acta Physica Sinica(物理学报)，2015，64(5)：054203.

[12] LIU Tong, XU Zheng-wu, WU Yuan-jie, et al(刘 通，徐政五，吴元杰，等).Journal of Signal Processing, 2013, 29(12): 1650.

[13] Jiang L, Li C, Huang L.Advance Journal of Food Science and Technology, 2014, 6(5): 578.

(Received May 11，2015; accepted Sep.19, 2015)

*Corresponding author

Research on the Method of Improving Terahertz Frequency Resolution Based on Empirical Mode Decomposition

LU Min1, LI Xiao-xia1，2*, SHANG Li-ping2, DENG Hu1，2

1.School of Information Engineering, Southwest University of Science and Technology, Mianyang 621010, China 2.Robot Technology Used for Special Environment Key Laboratory of Sichuan Province, Mianyang 621010, China

Terahertz frequency resolution is an important factor affecting substance identification, but the presence of the device reflection within the terahertz time domain spectroscopy systems, causes the presence of a plurality of reflection peaks in the reference and measurement signals with low frequency resolution, because of the length truncation of time domain signal.In order to remove the influence of the reflection peak, this paper proposes a method based on empirical mode decomposition to remove the time domain reflection peak, to improve the terahertz frequency resolution.The time domain reflection peaks are positioned by correlation with the real terahertz peak, calculating the reflection peak upper and lower envelope and an average, obtaining intrinsic mode functions, and alternating reflection peaks with intrinsic mode functions, the effective length of the time domain signal is increased to improve the terahertz frequency resolution.Water vapor in the air terahertz transmission spectroscopy results show that this method can self adaptively remove a plurality of reflection peaks and has good repairing effect on the terahertz time-domain signal; the frequency resolution is increased by 12 times; and the useful information of absorption spectrum is not lost; the absorption peak position and the number is consistent with real spectrum; the terahertz spectrum ability to identify is well preserved.

Terahertz time-domain spectroscopy (THz-TDS)；Empirical mode decomposition (EMD)；Frequency resolution；Reflection peak removing

2015-05-11，

2015-09-19

国家自然科学基金项目(11176032)，国防技术基础项目(Z202D13T001)，四川省教育厅资助科研项目(11ZB106)，极端条件物质特性实验室课题(13zxjk02)资助

卢 敏，女，1990年生，西南科技大学信息工程学院硕士 e-mail：121188438@qq.com *通讯联系人 e-mail：664368504@qq.com

O433.4

A

10.3964/j.issn.1000-0593(2016)09-2732-04