一种基于太赫兹反射光谱提取固体薄片复折射率的方法

摘要

本申请公开了一种基于太赫兹反射光谱提取固体薄片复折射率的方法，制备同一种材料的薄厚不同的两个固体薄片样品，其中薄样品为被测样品，厚样品作为时域信号基线校正时的辅助样品，测量二者的太赫兹时域反射信号；用被测样品和辅助样品的太赫兹时域反射信号之差作为基线，将被测样品的第一、二个峰分别作为基线校正后的参考信号和样品信号；对基线校正后的样品信号和参考信号分别进行复合数字滤波，以去除系统回波和随机噪声；根据滤波后的样品信号和参考信号，用自参考方法提取被测样品的复折射率。

说明书

技术领域

本申请涉及光子学技术与电子学技术领域，尤其涉及一种基于太赫兹反射光谱提取固体薄片复折射率的方法。

背景技术

太赫兹波通常指的是频率在0.1THz～10THz(1THz＝1012Hz)，波长在3mm～30μm范围内的电磁辐射。太赫兹波在电磁波谱上位于毫米波和红外线之间，属于远红外波段，太赫兹频段在电磁波谱中的位置如图1所示，是光子学技术与电子学技术、宏观与微观的过渡区域。太赫兹辐射具有非常重要的科学价值，目前已成为世界性的研究热点，是2l世纪科学研究最前沿的领域之一。太赫兹频段所具有的独特性质决定了它具有广阔的应用前景。

利用一些光谱技术可以对药物进行检测，目前已采用拉曼散射光谱、红外光谱对药物进行了一些研究。相比于这些技术，太赫兹时域光谱技术具有独特的优势，例如，它对化合物晶型有很高的灵敏度，能够反映整个分子的低频振动和分子间的弱相互作用，能量低不会造成药物发生化学变化等。因此，利用太赫兹光谱在化学药品的生产、储存和流通环节进行质量的检查等，成为太赫兹辐射最具前景的应用之一。

化学药品通常以固体薄片、颗粒和液态形式存在。通过对化学药品薄片的太赫兹时域信号进行快速傅立叶变换，可以同时获得它的太赫兹信号的振幅和相位信息，从而提取到物质在太赫兹频段的复折射率，即获得药品的“指纹谱”，从而达到分析和鉴别药品的目的。

现有技术中，通常采用全反射镜参考方法提取固体薄片的复折射率。图2(a)所示为太赫兹脉冲垂直入射到固体薄片样品上时光线的传播情况，图中表示被测固体薄片放置在干燥空气中。若忽略空气对太赫兹辐射的吸收，则空气的折射率为实数，且近似为1；设薄片的复折射率为其中n_s为样品的折射率，k_s为消光系数。

现有的全反射镜参考方法的处理流程如图3所示，包括如下步骤：

步骤301：基于反射式太赫兹时域光谱系统，测量全反射镜的时域信号并将其作为参考信号对其做傅里叶变换得到参考信号的频谱

步骤302：将全反射镜取下换成样品，并使得二者的前表面位置完全重合，测量样品的时域信号同样再对其做傅里叶变换，得到样品的频谱

步骤303：根据公式(6)计算得到样品的折射率n_s和消光系数k_s，从而获得复折射率

假设全反射镜无损耗反射，则全反射镜的反射信号等于入射信号以图2(a)中作为样品信号，将作为参考信号。则参考信号和样品信号的频谱可以分别表示为

${\tilde{E}}_r\left(\omega \right)={\tilde{E}}_0\left(\omega \right),---\left(1\right)$

${\tilde{E}}_s\left(\omega \right)={\tilde{E}}_1\left(\omega \right)={\tilde{E}}_0\left(\omega \right){\tilde{r}}_{as},---\left(2\right)$

其中，可由菲涅尔(Fresnel)公式得到，通常采用P分量探测，则

${\tilde{r}}_{as}=\frac{{\tilde{n}}_s\left(\omega \right)-1}{{\tilde{n}}_s\left(\omega \right)+1}.---\left(3\right)$

用样品信号频谱除以参考信号频谱可得到传递函数的表达式

$\tilde{H}\left(\omega \right)=\frac{{\tilde{E}}_s\left(\omega \right)}{{\tilde{E}}_r\left(\omega \right)}=\frac{{\tilde{E}}_0\left(\omega \right){\tilde{r}}_{as}}{{\tilde{E}}_0\left(\omega \right)}={\tilde{r}}_{as}.---\left(4\right)$

联立公式(3)和(4)，得到样品复折射率和传递函数的关系

${\tilde{n}}_s\left(\omega \right)=\frac{1+\tilde{H}\left(\omega \right)}{1-\tilde{H}\left(\omega \right)}.---\left(5\right)$

是复数，令则由公式(5)可以推导出样品的折射率和消光系数的表达式

$\left(\begin{array}{c}{n}_s\left(\omega \right)=\frac{1+\rho \left(\omega \right)}{1+\rho \left(\omega \right)-2\sqrt{\rho \left(\omega \right)}\cos \phi \left(\omega \right)}\\ k_s\left(\omega \right)=\frac{2\sqrt{\rho \left(\omega \right)}\sin \phi \left(\omega \right)}{1+\rho \left(\omega \right)-2\sqrt{\rho \left(\omega \right)}\cos \phi \left(\omega \right)}\end{array}\right).---\left(6\right)$

吸收系数α_s＝2ωk_s/c，其中c是真空中的光速。传递函数的振幅ρ(ω)和相位φ(ω)都是可测量的量，因此，由公式(6)便可以获得折射率n_s和消光系数k_s，从而获得复折射率${\tilde{n}}_s\left(\omega \right)=n_s+{\mathrm{jk}}_s.$

这种全反射镜参考方法存在如下问题：首先，传递函数的相位φ(ω)不可避免地会存在一个机械相移误差。如图2(b)所示，样品表面和全反射镜的表面很难实现完全复位，二者的位置会存在ΔL的光程差。因此，参考信号与样品信号会存在一个固有的相移使得提取的固体薄片的复折射率不够精确；其次，该方法要求测量系统的机构精密度足够高，因此，增加了系统硬件制作的难度；最后，由于需要分别测量参考与样品的时域信号，因此，测量效率低。

发明内容

本申请提供了一种基于太赫兹反射光谱提取固体薄片复折射率的方法，可以提高测量得到的复折射率的精确度，并提高测量效率。

本申请实施例提供的一种基于太赫兹反射光谱提取固体薄片复折射率的方法，包括：

A、制备同一种材料的薄厚不同的两个固体薄片样品，其中薄样品为被测样品，厚样品作为时域信号基线校正时的辅助样品，测量二者的太赫兹时域反射信号；

B、用被测样品和辅助样品的太赫兹时域反射信号之差作为基线，将被测样品的第一、二个峰分别作为基线校正后的参考信号和样品信号；

C、对基线校正后的样品信号和参考信号分别进行复合数字滤波，以去除系统回波和随机噪声；

D、根据滤波后的样品信号和参考信号，用自参考方法提取被测样品的复折射率，其中包括：采用遗传算法对公式

进行优化求解，以传递函数为目标函数，进行至少N次种群迭代，最终得到优化后的复折射率的精确解；其中，为传递函数的相位；为传递函数的振幅；菲涅尔系数r_as、t_as、r_sa、t_sa可由菲涅耳公式得到：

$r_{as}=\frac{{\tilde{n}}_s\left(\omega \right)-1}{{\tilde{n}}_s\left(\omega \right)+1},t_{as}=\frac{2}{1+{\tilde{n}}_s\left(\omega \right)},r_{sa}=\frac{1-{\tilde{n}}_s\left(\omega \right)}{1+{\tilde{n}}_s\left(\omega \right)},t_{sa}=\frac{2{\tilde{n}}_s\left(\omega \right)}{1+{\tilde{n}}_s\left(\omega \right)};$表示复折射率，c是真空中的光速，d表示被测样品厚度，下标as表示由空气到样品的界面，下标sa表示由样品到空气的界面。

较佳地，步骤B包括：

B1、将被测样品和辅助样品的太赫兹时域反射信号的第一个峰在时间轴上对齐；

B2、用被测样品的太赫兹时域反射信号减去辅助样品的太赫兹时域反射信号，得到样品信号；

B3、将被测样品的太赫兹时域反射信号中的第二个峰置为前后数据点的平均值，将辅助样品的太赫兹的时域反射信号中的第一个峰置为前后数据点的平均值，再用被测样品的太赫兹时域反射信号减去辅助样品的太赫兹的时域反射信号，得到参考信号。

较佳地，所述复合数字滤波包括解卷积和小波变换。

较佳地，所述解卷积包括：

通过实验得到第一个回波的振幅A₁和第一个回波与主峰的时间延迟(t₁-t₀)；

依据公式滤波，获得主峰的频谱；其中，A(ω)和A₀(ω)分别为原始时域信号和主峰的频谱。

较佳地，所述小波变换包括：将基线校正后的样品信号和参考信号多尺度分解，使噪声按频率区域化分；采用Matlab中的wdencmp函数，对高频系数进行阈值量化去噪处理，最后根据小波分解的各层系数和经过量化处理后的高频系数，进行一维信号的重构，得到去除噪声后的样品信号和参考信号。

较佳地，所述N大于或等于100。

较佳地，由公式计算得到的n₀(ω)、k₀(ω)作为遗传算法的初代种群。

较佳地，所述被测样品和辅助样品的厚度差的绝对值等于或大于0.8mm。

从以上技术方案可以看出，与目前的全反射镜参考方法相比，本发明的方法具有以下优点：

(1)仅需一次性制备同种材料的厚样品，以后再测相同材料时可重复使用作为基线校正的辅助样品，从而提高了测量效率；

(2)采用数字滤波技术，可以通过编程灵活地调节参数，不受硬件的限制，并且可以采用最先进的数字滤波算法；

(3)无需用全反射镜作为参考，避免了样品表面与参考物不共面所带来的机械相移误差，提高了复折射率的测量精度；

(4)将被测固体薄片的前表面反射的光作为参考信号，将后表面反射的光作为样品信号，光经过样品内部后将携带更多的样品信息，便于提取样品的特征吸收峰。

附图说明

图1为太赫兹频段在电磁波谱中的位置示意图；

图2(a)为太赫兹脉冲垂直入射到固体薄片上的光线传播示意图，图2(b)为相应的反射相移示意图；

图3为现有的全反射镜参考方法的处理流程示意图；

图4为本申请提供的基于太赫兹反射光谱提取固体薄片复折射率的方法的流程示意图；

图5所示为全反射镜与厚度为2.588mm、直径为13mm的葡萄糖薄片的太赫兹时域信号；

图6所示为采用现有的全反射镜参考方法提取的葡萄糖薄片的折射率(a)和消光系数(b)；

图7所示为厚度分别为2.588mm和3.750mm的葡萄糖薄片的太赫兹反射时域信号；

图8所示为厚度为2.588mm的葡萄糖薄片去除基线前后的时域信号：(a)样品信号；(b)参考信号；

图9所示为厚度为2.588mm的葡萄糖薄片第一个峰的频谱复合滤波前后的比较：(a)测量的原始频谱；(b)小波去噪后的频谱；(c)解卷积后的频谱；(d)复合滤波后的频谱；

图10所示为厚度为2.588mm的葡萄糖薄片的折射率和吸收系数随频率变化曲线。

具体实施方式

为使本申请技术方案的技术原理、特点以及技术效果更加清楚，以下结合具体实施例对本申请技术方案进行详细阐述。

图4为本申请提供的基于太赫兹反射光谱提取固体薄片复折射率的方法的流程示意图，该流程包括：

步骤401：制备同一种材料的薄厚不同的两个固体薄片样品，测量二者的太赫兹时域反射信号，其中薄样品为被测样品，厚样品作为时域信号基线校正时的辅助样品，薄厚样品的厚度差的绝对值应等于或大于0.8mm；

步骤402：用被测样品和辅助样品的太赫兹时域反射信号之差作为基线，将被测样品的第一、二个峰分别作为基线校正后的参考信号和样品信号；

步骤403：对基线校正后的样品信号和参考信号分别进行复合数字滤波，以去除系统回波和随机噪声；

步骤404：根据滤波后的样品信号和参考信号，用自参考方法提取被测样品的复折射率，其中包括：采用遗传算法进行优化求解，以传递函数为目标函数，进行至少N次种群迭代，最终得到优化后的复折射率的精确解。其中，N大于等于100。

其中，步骤402的作用是消除基线噪声。具体可以包括如下子步骤：

子步骤402-1：将被测样品和辅助样品的太赫兹时域反射信号的第一个峰在时间轴上对齐；

子步骤402-2：用被测样品的太赫兹时域反射信号减去辅助样品的太赫兹时域反射信号，得到样品信号；

子步骤402-3：将被测样品的太赫兹时域反射信号中的第二个峰置为前后数据点的平均值，将辅助样品的太赫兹的时域反射信号中的第一个峰置为前后数据点的平均值，再用被测样品的太赫兹时域反射信号减去辅助样品的太赫兹的时域反射信号，得到参考信号。

步骤403中的复合数字滤波包括解卷积和小波变换。

太赫兹波在砷化镓(GaAs)基底和碲化锌(ZnTe)探测晶体中的多次反射会产生系统回波，给太赫兹光谱带来系统误差。本申请方案中利用解卷积算法去除系统回波。

假定考虑到测量是在干燥空气条件下进行，认为THz脉冲在真空中传播的过程中无耗散，因此系统回波脉冲等于主脉冲的相对强度乘以时间延迟：

$A\left(t\right)=A_0\left(t\right)\otimes \underset{i=0}{\overset{\infty }{\Sigma }}{A}_i\delta \left(t_i\right),---\left(8\right)$

其中，A(t)和A₀(t)分别为原始时域信号和主脉冲；A_i和t_i分别表示第i个回波相对于主脉冲的幅值和延迟时间，当i＝0表示主脉冲；为卷积运算符。理论上说，系统回波的数目是无限的，其振幅依次衰减，计算时考虑的回波数目越多计算结果越精确，但是计算过程越复杂。本申请实施例中，只考虑较强的第一个系统回波，即取i＝1。将时域信号表达式(8)进行傅里叶变换，并取i＝1，可以得到

$A\left(\omega \right)=\frac{A_0\left(\omega \right)}{1+A_1{e}^{-i\omega (t_1-t_0)}},---\left(9\right)$

其中，A(ω)和A₀(ω)分别为原始时域信号和主峰的频域表达式；(t₁-t₀)为系统的第一个回波与主峰的时间延迟。由公式(9)可以看出，A(ω)相当于在A₀(ω)基础上加正弦振荡。因此，只要通过实验得到第一个回波的振幅A₁和其与主峰的时间延迟(t₁-t₀)，就可以去除系统回波噪声。

本申请方案利用小波变换去除随机噪声，小波变换是把函数在特定的空间内按照小波基展开和逼近，可以同时进行时域和频域的分析，并具有很强的表征信号局部特征的能力。因为小波变换去噪对于具有尖峰和突变特性的太赫兹时域信号保真度高。根据小波变换理论将噪声信号进行多尺度分解，使噪声按频率区域化分，再分别用滤波器对信号进行滤波，用得到的近似系数和精确系数进行小波重构，最后得到去除随机噪声的真实时域信号。对于信号f(t)，其小波变换的形式为

${\omega }_f(a,b)=\int f\left(t\right){\psi }_{(a,b)}\left(t\right)dt=\frac{1}{\sqrt{\left|a\right|}}\int f\left(t\right)\psi \left(\frac{t-b}{a}\right)dt,---\left(10\right)$

其中，ψ_(a,b)(x)为母小波ψ(x)经过伸缩和平移得到的小波函数，a表示伸缩因子，b表示平移因子。小波变换ω_f(a,b)体现的是原信号在x＝b处、在|x-b|≤a范围内对信号进行分析，具有时频局域化和多分辨的特性。

一个含噪声的一维时域信号可以认为是真实信号叠加了噪声信号，以如下方式表示：

S(t)＝f(t)+N(t),t＝0,1,...,n-1,(11)

其中，S(t)为测量信号，f(t)为真实信号，N(t)为噪声信号。在实际应用中，噪声多表现为高频信号。因此，在用小波分析时，将一维时域信号S(t)(基线校正后的时域样品信号和时域参考信号)进行分解后，采用Matlab中的wdencmp函数，对高频系数进行阈值量化去噪处理，最后根据小波分解的各层系数和经过量化处理后的高频系数，进行一维信号的重构，得到去除噪声后的时域样品信号和时域参考信号。

以下给出步骤404所述的一种提取固体薄片太赫兹频段复折射率的自参考方法的具体实现过程，该方法以图2中的作为时域参考信号，而以作为时域样品信号，将二者分别进行傅里叶变换，得到各自的频谱复折射率的求解公式推导过程如下

${\tilde{E}}_r\left(\omega \right)={\tilde{E}}_1\left(\omega \right)={\tilde{E}}_0\left(\omega \right)r_{as},---\left(8\right)$

${\tilde{E}}_s\left(\omega \right)={\tilde{E}}_2\left(\omega \right)={\tilde{E}}_0\left(\omega \right)t_{as}{p}_s[\omega ,L_s\left(\omega \right)]r_{sa}{p}_s[\omega ,L_s\left(\omega \right)]t_{sa}.---\left(9\right)$

其中，P分量菲涅尔系数r_as、t_as、r_sa、t_sa可由菲涅耳公式得到

$r_{as}=\frac{{\tilde{n}}_s\left(\omega \right)-1}{{\tilde{n}}_s\left(\omega \right)+1},t_{as}=\frac{2}{1+{\tilde{n}}_s\left(\omega \right)},r_{sa}=\frac{1-{\tilde{n}}_s\left(\omega \right)}{1+{\tilde{n}}_s\left(\omega \right)},t_{sa}=\frac{2{\tilde{n}}_s\left(\omega \right)}{1+{\tilde{n}}_s\left(\omega \right)}.---\left(10\right)$

下标as表示由空气到样品的界面，下标sa表示由样品到空气的界面。

传播因子p_s[ω,L_s(ω)]的表达式为，

$\left(\begin{array}{c}{p}_s[\omega ,L_s\left(\omega \right)]=\exp \left(\frac{-{\mathrm{jL}}_s\left(\omega \right)\omega }{c}\right)\\ L_s\left(\omega \right)={\tilde{n}}_s\left(\omega \right)·d\end{array}\right).---\left(11\right)$

d表示被测样品厚度。

传递函数为

$\stackrel{···}{H}\left(\omega \right)=\frac{{\tilde{E}}_s\left(\omega \right)}{{\tilde{E}}_r\left(\omega \right)}=\frac{t_{as}{r}_{sa}{t}_{sa}}{r_{as}}\exp \left(\frac{-2j{\tilde{n}}_s\left(\omega \right)d\omega }{c}\right).---\left(12\right)$

将代入上式，并化简得到

$\tilde{H}\left(\omega \right)=\frac{t_{as}{r}_{sa}{t}_{sa}}{r_{as}}·\exp (-2k_s\left(\omega \right)·\omega d/c)·\exp (-2{\mathrm{jn}}_s\left(\omega \right)·\omega d/c).---\left(13\right)$

若令则有

$\rho \left(\omega \right)=abs\left(\frac{t_{as}{r}_{sa}t}{r_{as}}\right)·\exp (-\frac{2k_s\left(\omega \right)·\omega d}{c}),---\left(15\right)$

由上述等式可推导出

abs()表示对括号内求绝对值，ln[]表示对括号内求自然对数。

对于弱吸收介质，即n_s>>k_s，可以忽略样品复折射率的虚部k_s，则$\rho \left(\omega \right)=\frac{t_{as}{r}_{sa}{t}_{sa}}{r_{as}^s}·\exp (-\frac{2k_s\left(\omega \right)·\omega d}{c}),r_{as}=\frac{n_s\left(\omega \right)-1}{n_s\left(\omega \right)+1},t_{as}=\frac{2}{1+n_s\left(\omega \right)},r_{sa}=\frac{1-n_s\left(\omega \right)}{1+n_s\left(\omega \right)},t_{sa}=\frac{2n_s\left(\omega \right)}{1+n_s\left(\omega \right)},$从而得到样品复折射率的近似解

利用公式(16)，采用遗传算法，通过MATLAB编程计算，便可以提取出样品的复折射率。

较佳地，将由公式(17)计算得到的n₀(ω)、k₀(ω)作为遗传算法的初代种群，这样便于提高计算效率，但是，最终的优化结果与赋给的初值无关。

该自参考方法不仅避免了反射太赫兹时域光谱系统的固有相移问题，而且由于光程增加而提取了更强的物质属性信息。另外，作为基线校正的厚样品一次性制备完成后，或者当测量基线完好，并不需要校正时，仅需测量一次样品的时域信号，便可以获得复折射率的值。因此，实验效率提高了。这对于太赫兹光谱技术的实际应用具有重要意义。本发明的方法不局限于固体薄片，也可以推广到液态样品复折射率的提取。若将颗粒样品压制成薄片，同样可以应用本发明方法提取颗粒样品的复折射率。

以下以厚度为2.588mm、直径为13mm的无水D-葡萄糖粉末的固体薄片的实验为例，说明本申请方法与传统方法的区别。

I.现有的全反射镜参考方法

现有的全反射镜参考方法要求参考全反射镜表面与样品前表面的位置完全重合。图5为全反射镜与厚度为2.588mm、直径为13mm的葡萄糖薄片的太赫兹时域信号，小图为第一个峰放大后的曲线。如图5中的小图可见，从波峰位置看到，金镜和各样品的表面并没有完全重合。利用公式(1)-(6)提取的葡萄糖薄片的折射率和消光系数，如图6所示。

II.本发明方法

1)光谱相减法基线校正

在制备了被测样品(厚度为2.588mm、直径为13mm的葡萄糖粉末的固体薄片)后，还制备了厚度为3.750mm、直径为13mm的同一种葡萄糖粉末的固体薄片，用于时域信号基线校正时的辅助样品。

图7展示的分别是纯葡萄糖不同厚度薄片的两条时域信号，粗线、细线分别是厚度为3.750mm和2.588mm的薄片的时域信号，图中的虚线为零线。可以明显看到，随着测量时间的增加，谱线向下漂移，表现出基线存在误差。本发明对原始数据进行如下基线校正：

(1)将厚度不同的纯葡萄糖薄片的时域信号的第一个峰位通过调节时间手动重合；然后，根据数据点统一将其时间截到60ps；最后，用薄样品的时域信号减去厚样品的时域信号，结果如图8(a)中细线所示，粗线为原始时域信号；

(2)将薄样品时域信号中的第二个峰置为前后数据点的平均值，将厚样品时域信号中的第一个峰置为前后数据点的平均值；然后，用薄样品的时域信号减去厚样品的时域信号，即得到参考信号。如图8(b)中的细线所示，粗线为测量的原始时域信号。

综上所述，采取不同厚度的样品时域谱来消除基线噪声方法的具体使用步骤可以归纳为：首先，将不同厚度样品的两条时域谱的第一个峰手动对齐；然后，用薄样品的时域信号减去厚样品的时域信号，得到样品信号；最后，将薄样品时域信号中的第二个峰置为前后数据点的平均值，将厚样品时域信号中的第一个峰置为前后数据点的平均值，用薄样品的时域信号减去厚样品的时域信号，得到参考信号。

由图4可见，虽然经过基线校正后，系统回波信号幅度有所减弱，但是依然存在。因此，还要考虑消除系统回波。

2)复合滤波法去噪

图9是厚度为2.588mm厚的葡萄糖薄片的复合滤波前后的比较。按照上述解卷积和小波变换的原理进行滤波去噪，得到如图9所示的第一个峰的傅里叶变换结果。

3)自参考法提取复折射率

将参考信号和样品信号进行傅里叶变换获得各自的频谱图，两式相除获得传递函数，按照自参考的原理计算获得厚度为2.588mm厚的葡萄糖薄片的折射率和吸收系数，结果如图10所示。可以看出，葡萄糖的折射率在1.3THz后呈现反常色散，且吸收系数随着频率的增加而增大，并存在特征吸收峰。

III.本发明方法与传统全反射镜参考方法结果比较

图6和图10分别为用传统的全反射镜参考方法和本发明方法得到的葡萄糖薄片的折射率和消光系数。由图6(a)可见，折射率变化比较平缓，图6(b)中的消光系数也没有表现出明显的特征吸收，二者没有满足的介电模型。然而，从图10(a)中可见，折射率有多处突变，且突变位置与图10(b)中的消光系数的峰位相关，峰位在图中标出，二者满足洛伦兹介电模型。这些吸收峰是无水D-葡萄糖的特征吸收峰，从特征吸收峰我们可以鉴别是否为无水D-葡萄糖这种物质。

以上所述仅为本申请的较佳实施例而已，并不用以限制本申请的保护范围，凡在本申请技术方案的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本申请保护的范围之内。