CN103645154A - 一种利用太赫兹光谱信号提取材料光学常数的方法 - Google Patents

Abstract

本发明属于光电子学与材料学应用技术领域，具体涉及一种利用太赫兹光谱信号提取材料光学常数的方法，包括以下步骤：步骤一，太赫兹波垂直入射到放置于空气中厚度为l的平板样品中，根据Fabry-Perot腔模型获得样品信号和参考信号太赫兹频谱的理论比值A(ω)；步骤二，利用模拟退火优化算法获得样品信号和参考信号太赫兹频谱理论比值A(ω)和实验比值A_meas(ω)之间的全局最小值，获得最小值处对应的实折射率n(ω))和消光系数κ(ω)。该方法采用精确公式结合模拟退火优化算法提取材料光学常数，可以取消弱吸收近似和材料边界约束条件，适用于不同厚度的材料，采用该方法提取的材料光学常数精确度较高。

Description

一种利用太赫兹光谱信号提取材料光学常数的方法

技术领域

本发明属于光电子学与材料学应用技术领域，具体涉及一种利用太赫兹光谱信号提取材料光学常数的方法。

背景技术

目前，各种光谱技术已经广泛的应用于物质分析，采用太赫兹光谱信号提取材料光学常数的方法已经成为一种专业领域内的常用方法。现有的太赫兹时域光谱系统的光路如图1所示，该系统用于测量样品透射信号和参考信号的太赫兹光谱。常用提取材料光学常数的步骤如下：

1.通过实验测得样品信号和参考信号的太赫兹频谱，将两者相比表示为：

A_meas(ω)=E_sam(ω)/E_ref(ω)=T_meas(ω)exp[-jΔφ_meas(ω)]     (1)

式中E_sam(ω)和E_ref(ω)为频域中的复电场，分别表示样品信号和参考信号。T_meas(ω)为透过样品的太赫兹电场，Δφ_meas(ω)为相位变化。

2.根据实验测得的T_meas(ω)和Δφ_meas(ω)，可确定实折射率n(ω)、消光系数κ(ω)和吸收系数α(ω)分别为：

$n\left(\mathrm{\ω}\right)=1+\mathrm{\Δ}{\mathrm{\φ}}_{\mathrm{meas}}\left(\mathrm{\ω}\right)\frac{c}{\mathrm{\ωl}}---\left(2\right)$

$\mathrm{\κ}\left(\mathrm{\ω}\right)=\ln \left\{\frac{4n\left(\mathrm{\ω}\right)}{T_{\mathrm{meas}}\left(\mathrm{\ω}\right){[1+n\left(\mathrm{\ω}\right)]}^2}\right\}\frac{c}{\mathrm{\ωl}}---\left(3\right)$

$\mathrm{\α}\left(\mathrm{\ω}\right)=\frac{2\mathrm{\ω\κ}\left(\mathrm{\ω}\right)}{c}=\frac{2}{l}\ln \left\{\frac{4n\left(\mathrm{\ω}\right)}{T\left(\mathrm{\ω}\right){[1+n\left(\mathrm{\ω}\right)]}^2}\right\}---\left(4\right)$

其中，l为样品厚度，c为真空中光速。

综上所述，常用技术流程如图2所示。现有技术的利用太赫兹光谱信号提取材料光学常数的方法存在几个问题：其一，进行了弱吸收近似（n>>κ）；其二，对所测试的材料有一定的边界约束，如较高的折射率和/或较厚的样品；其三，这些算法的精度还受到被测信号的不确定性限制。

发明内容

本发明的目的在于针对要解决的技术问题，提供一种利用太赫兹光谱信号提取材料光学常数的方法，该方法可以取消弱吸收近似和材料边界约束条件，获得的材料光学常数精确度较高。

一种利用太赫兹光谱信号提取材料光学常数的方法，包括以下步骤：步骤一，太赫兹波垂直入射到放置于空气中厚度为l的平板样品中，根据Fabry-Perot腔模型获得样品信号和参考信号太赫兹频谱的理论比值A(ω)；具体为：

$A\left(\mathrm{\ω}\right)=\frac{4{\tilde{n}}^2\left(\mathrm{\ω}\right)}{{(\tilde{n}\left(\mathrm{\ω}\right)+1)}^2}\exp [-j(\tilde{n}\left(\mathrm{\ω}\right)-1)\mathrm{\ωl}/c]\mathrm{FP}\left(\mathrm{\ω}\right)---\left(a\right)$

其中，为复折射率，n(ω)为实折射率、κ(ω)为消光系数，c为真空中光速，FP(ω)为Fabry-Perot系数，FP(ω)的表达式为：

$\mathrm{FP}\left(\mathrm{\ω}\right)=\underset{p=0}{\overset{m}{\mathrm{\Σ}}}{\left[\frac{(\tilde{n}\left(\mathrm{\ω}\right)-1)}{(\tilde{n}\left(\mathrm{\ω}\right)+1)}\exp (-\mathrm{j\ω}\tilde{n}l)\right]}^{2p}---\left(b\right)$

其中，p为整数，表示反射回波级次，p=0表示主反射回波，m表示最大反射回波级次；步骤二，利用模拟退火优化算法获得样品信号和参考信号太赫兹频谱理论比值A(ω)和实验比值A_meas(ω)之间的全局最小值，获得最小值处对应的实折射率n(ω))和消光系数κ(ω)；具体步骤如下：

（21）构建精确频谱和测量频谱之间的误差函数：

根据样品厚度l、样品信号与参考信号的时间延迟Δt、时间窗口大小t_c，获得估计折射率n_g和最大反射回波级次m值：

$n_g=\frac{\mathrm{c\Δt}}{l}+n_0---\left(c\right)$

$m\≤\frac{c{t}_c}{2n_{g}l}-\frac{1}{2}---\left(d\right)$

其中，时间窗口大小是指从参考信号入射到测量结束的时间，n₀为太赫兹波在空气中的折射率，c为太赫兹波空气中的速度；

A(ω)=T(ω)exp[-jΔφ(ω)]          （e）

同理：A_meas(ω)=T_meas(ω)exp[-jΔφ_meas(ω)]       （f）

根据样品信号峰值|E_s,max|和参考信号峰值|E_r,max|可以求出消光系数κ_g的近似表达式：

${\mathrm{\κ}}_g=\frac{c}{\mathrm{\ωl}}\ln \frac{\left|E_{r,\max }\right|}{\left|E_{s,\max }\right|}---\left(g\right)$

结合公式（e）、（f）、（g）构建一个精确频谱和测量频谱之间的误差函数：

δ(n,κ)={ln[T(ω)]-ln[T_meas(ω)]}²+{ln[Δφ(ω)]-ln[Δφ_meas(ω)]}²     （h）

（22）利用模拟退火算法求无约束条件二元函数的全局最小值，获得δ(n,κ)的全局最小值，进一步得到最小值处对应的实折射率n(ω)和消光系数κ(ω)：

（221）输入初始变量：选取初始点x₀，使用估计折射率n_g和估计消光系数κ_g，估计最小值区间的上、下边界u和l，选定最大迭代数k_max>0，退火因子1>q>0，以及相应的函数值浮动误差容许值ε_f；

(222)求解基本变量：x=x₀，x⁰=x，f⁰=f(x⁰)；

(223)进行迭代计算，找出近似全局最小值，即从k=0到k=k_max执行以下循环块：

生成一个在[-1,1]上服从均匀分布的N×1的随机向量y，y与x大小相同，利用μ^-1定理可以求出Δx：

Δx=(u-l)·g_μ ^-1(y)， ${g_{\mathrm{\μ}}}^{-1}\left(y\right)=\frac{{(1+\mathrm{\μ})}^{\left|y\right|}-1}{\mathrm{\μ}}\mathrm{sign}\left(y\right),$ μ=10¹⁰⁰(k/k_max)^q；

在(l,u)区间内确定x₁=x+Δx；

如果Δf=f(x₁)-f(x)<0，或者在[0,1]上生成的均匀分布的随机数z<p=exp[-(k/k_max)^qΔf/|f(x)|/ε_f]，则令x=x₁，f(x)=f(x₁)。

如果f(x)<f⁰，则令x⁰=x，f⁰=f(x⁰)；

令k=k+1，继续循环；

(224)根据步骤（223）求出近似全局最小值x⁰后，假定x⁰为初值，进而求出准确的全局最小值；

其中，x₀，u，l，x，x⁰，x₁，Δx和g_μ(y)均为矩阵向量。

其中，还包括步骤三，利用Nelder-Mead单纯形算法求一元函数的局部最小值，获得实折射率n(ω)和消光系数κ(ω)对样品厚度变化的最小值，确定样品厚度l，并修正实折射率n(ω)和消光系数κ(ω)；具体为：

（31）将样品从厚度范围l_min到l_max分为I等份，i表示的厚度为

构造一元误差函数：

$\mathrm{\&epsiv;}\left(l\right)=\underset{i=1}{\overset{I}{\mathrm{\&Sigma;}}}\left[\right|n(i-1)-n\left(i\right)|+|\mathrm{\&kappa;}(i-1)+\mathrm{\&kappa;}\left(i\right)\left|\right]---\left(i\right)$

其中，n(i)，κ(i)分别表示厚度为l_i的实折射率和消光系数；

（32）准确样品厚度确定法：利用Nelder-Mead单纯形算法求公式（i）的一元函数的最小值，准确确定样品厚度l，并修正各频率点的实折射率n(ω)和消光系数κ(ω)；具体求解过程如下：

(321)根据不同的厚度值l_i，按模拟退火法求出对应的实折射率n(i)和消光系数κ(i)；

(322)用Nelder-Mead单纯形算法求出一元误差函数ε(l)的最小值及其对应的厚度l；

(323)代入l修正各频率点的实折射率n(ω)和消光系数κ(ω)。

其中，所述步骤(224)求出准确的全局最小值采用的方法为最速下降法、牛顿迭代法或Nelder-Mead单纯形算法中的任意一种。

本发明的有益效果在于：本发明的利用太赫兹光谱信号提取材料光学常数的方法，采用精确公式结合模拟退火优化算法提取材料光学常数，可以取消弱吸收近似和材料边界约束条件，适用于不同厚度的材料，采用该方法提取的材料光学常数精确度较高。

附图说明

图1为本发明的太赫兹时域光谱系统的光路图。

图2为现有技术的流程图。

图3为本发明的流程图。

图4为本发明的参考信号（虚线）和经过GaAs晶体后的样品信号（实线）。

图5为傅立叶变换后得到的参考信号振幅（虚线）和样品信号振幅（实线）。

图6为傅立叶变换后得到的参考信号相位（虚线）和样品信号相位（实线）。

图7为GaAs晶体在0.2—1.6THz波段的折射率。

图8为GaAs晶体在0.2—1.6THz波段的消光系数。

图9为GaAs晶体在0.2—1.6THz波段的吸收系数。

具体实施方式

为了使发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。

如图3和图9所示，一种利用太赫兹光谱信号提取材料光学常数的方法，包括以下步骤：

一.根据Fabry-Perot腔模型获得样品信号和参考信号太赫兹频谱的理论比值A(ω)；

二.利用模拟退火优化算法寻找样品信号和参考信号太赫兹频谱理论比值A(ω)和实验比值A_meas(ω)之间的全局最小值，获得最小值处对应的实折射率n(ω))和消光系数κ(ω)。

三.准确样品厚度确定法：利用Nelder-Mead单纯形算法求一元函数的局部最小值，获得实折射率n(ω)和消光系数κ(ω)对样品厚度变化的最小值，准确确定样品厚度l，并修正实折射率n(ω)和消光系数κ(ω)。

步骤一至步骤三的具体实施方式为：

1.假定太赫兹波垂直入射到放置于空气中厚度为l的平板样品，考虑太赫兹波在具有平行表面样品中传播的Fabry-Perot效应，可以获得样品信号和参考信号太赫兹频谱的理论比值A(ω)：

$A\left(\mathrm{\&omega;}\right)=\frac{4{\tilde{n}}^2\left(\mathrm{\&omega;}\right)}{{(\tilde{n}\left(\mathrm{\&omega;}\right)+1)}^2}\exp [-j(\tilde{n}\left(\mathrm{\&omega;}\right)-1)\mathrm{\&omega;l}/c]\mathrm{FP}\left(\mathrm{\&omega;}\right)---\left(5\right)$

其中，

为复折射率，n(ω)为实折射率、κ(ω)为消光系数，c为真空中光速，ω为角频率，即频率f的2π倍，j为虚数单位，FP(ω)为Fabry-Perot系数，其精确表达式为

$\mathrm{FP}\left(\mathrm{\&omega;}\right)=\underset{p=0}{\overset{m}{\mathrm{\&Sigma;}}}{\left[\frac{(\tilde{n}\left(\mathrm{\&omega;}\right)-1)}{(\tilde{n}\left(\mathrm{\&omega;}\right)+1)}\exp (-\mathrm{j\&omega;}\tilde{n}l)\right]}^{2p}---\left(6\right)$

其中，p为整数，表示反射回波级次，p=0表示主反射回波，m表示最大反射回波级次。

2.根据样品厚度l、样品信号与参考信号的时间延迟Δt、时间窗口大小t_c（时间窗口大小即从参考信号入射到测量结束的时间），获得估计折射率n_g和最大反射回波级次m值：

$n_g=\frac{\mathrm{c\&Delta;t}}{l}+n_0---\left(7\right)$

$m\&le;\frac{c{t}_c}{2n_{g}l}-\frac{1}{2}---\left(8\right)$

其中，n₀为太赫兹波在空气中的折射率，c为太赫兹波空气中的速度。由m值可以获得式(5)(6)表示的A(ω)。此外，根据样品信号峰值|E_s,max|和参考信号峰值|E_r,max|可以近似估计消光系数：

${\mathrm{\&kappa;}}_g=\frac{c}{\mathrm{\&omega;l}}\ln \frac{\left|E_{r,\max }\right|}{\left|E_{s,\max }\right|}---\left(9\right)$

3.已经获得A(ω)和A_meas(ω)，表示为：

A(ω)=T(ω)exp[-jΔφ(ω)]            (10)

A_meas(ω)=T_meas(ω)exp[-jΔφ_meas(ω)]        (11)

A_meas(ω)由实验测定；峰值信号|E_s,max|和参考峰值信号|E_r,max|也是由实验测出，如图4中的虚线和实线的最大值。

构建一个精确频谱和测量频谱之间的误差函数：

δ(n,κ)={ln[T(ω)]-ln[T_meas(ω)]}²+{ln[Δφ(ω)]-ln[Δφ_meas(ω)]}²   (12)

4.利用模拟退火算法求无约束条件二元函数的全局最小值，获得δ(n,κ)的全局最小值，进一步得到最小值处对应的实折射率n(ω)和消光系数κ(ω)。

模拟退火法是通过退火过程与最小化过程之间的比较，并根据概率确定局部最小点的取舍，从而跳过局部最小点，获得全局最小值，其求解过程如下：

(1)选取初始变量：选取初始点x₀，使用估计折射率n_g和估计消光系数κ_g，估计最小值区间的上下边界u和l，选定最大迭代数k_max>0，退火因子1>q>0，退火速度较快时，因子较大接近1，以及相应的函数值浮动误差容许值ε_f。x₀为一个复数，实部为折射率n，虚部为消光系数κ；计算时，该复数的实部为折射率n_g，虚部为消光系数κ_g；

(2)求解基本变量：x=x₀，x⁰=x，f⁰=f(x⁰)。x为迭代中间量，f(x)为公式（12）式表示的误差函数，是根据x₀求出误差函数（12）式的值f⁰，上、下边界u和l根据x₀和材料进行选定。

(3)进行迭代计算，找出近似全局最小值，即从k=0到k=k_max执行以下循环块：

生成一个在[-1,1]上服从均匀分布的N×1的随机向量y，y与x大小相同，利用μ^-1定理可以求出Δx：

Δx=(u-l)·g_μ ^-1(y)， ${g_{\mathrm{\&mu;}}}^{-1}\left(y\right)=\frac{{(1+\mathrm{\&mu;})}^{\left|y\right|}-1}{\mathrm{\&mu;}}\mathrm{sign}\left(y\right),$ μ=10¹⁰⁰(k/k_max)^q。

在(l,u)区间内确定x₁=x+Δx。

如果Δf=f(x₁)-f(x)<0，或者在[0,1]上生成的均匀分布的随机数z<p=exp[-(k/k_max)^qΔf/|f(x)|/ε_f]，则令x=x₁，f(x)=f(x₁)。

如果f(x)<f⁰，则令x⁰=x，f⁰=f(x⁰)。

令k=k+1，继续循环。

(4)根据上述步骤求出近似全局最小值x⁰后，假定x⁰为初值，并运用任意一种（最速下降法、牛顿迭代法或Nelder-Mead单纯形算法）求局部最小值的方法来寻找准确的全局最小值。

注意上述步骤中x₀，u，l，x，x⁰，x₁，Δx和g_μ(y)均为矩阵向量，因此，上述算法适用于多元函数求全局最小值。此外，通过调整最大迭代数k_max，退火因子q，以及相应的函数值浮动误差容许值ε_f，可以获得足够精确的收敛的全局最小值。

5.以上求出的实折射率n(ω))和消光系数κ(ω)对于频率来说不是光滑的曲线，为了获得光滑曲线，必须精确选定样品厚度。将样品厚度范围l_min到l_max分为I等份，i表示的厚度为

构造一元误差函数：

$\mathrm{\&epsiv;}\left(l\right)=\underset{i=1}{\overset{I}{\mathrm{\&Sigma;}}}\left[\right|n(i-1)-n\left(i\right)|+|\mathrm{\&kappa;}(i-1)+\mathrm{\&kappa;}\left(i\right)\left|\right]---\left(13\right)$

其中，n(i)，κ(i)分别表示厚度为l_i的实折射率和消光系数。

6.准确样品厚度确定法：利用Nelder-Mead单纯形算法求上述一元函数的最小值，准确确定样品厚度l，并修正各频率点的实折射率n(ω)和消光系数κ(ω)。

由于厚度l在l_min到l_max一个较小的范围内变化，可以利用Nelder-Mead单纯形算法求上述一元函数的最小值，准确确定样品厚度l，并修正各频率点的实折射率n(ω)和消光系数κ(ω)。求解过程如下：

(1)根据不同的厚度值l_i，可以通过不同的I值调整厚度l_i的步长，按模拟退火法求出对应的实折射率n(i)和消光系数κ(i)，此时角频率ω取固定值；

(2)用Nelder-Mead单纯形算法求出一元误差函数ε(l)的最小值及其对应的厚度l；

(3)代入l修正各频率点的实折射率n(ω)和消光系数κ(ω)，此时角频率ω在整个测试频谱范围内取值。

具体实施例一：

1.实验测试GaAs半导体样品信号和参考信号，样品厚度为0.63mm。图4为测量得到的时域上的参考信号（虚线）和样品信号（实线），图5、图6为经傅里叶变换后频域上的参考信号（虚线）和样品信号（实线）的振幅和相位。

2.根据样品厚度l、样品信号与参考信号的时间延迟Δt、时间窗口大小t_c（从参考信号入射到测量结束的时间），获得估计折射率n_g、最大反射回波级次m值和消光系数κ_g。

从图4可以看出放置样品前后测量得到的THz波形的振幅和形状都有变化，这是因为THz脉冲经过样品GaAs晶体时，信号波形携带有样品的色散信息，同时样品信号波形与参考波形相比延迟了Δt=5.47ps，这是由于样品的折射率比空气的折射率大，引起了附加的光程差，根据这个延迟时间可以得到样品的估计折射率n_g，见公式（7）有：

$n_g=\frac{3\&times;{10}^8\&times;5.47\&times;{10}^{-12}}{0.63\&times;{10}^{-3}}+1=3.603$

求得GaAs样品的估计折射率为3.603，这个值可用作数据处理初值，也可以与数据处理结果进行相比较。

由(8)式，最大反射回波级次m值满足

$m\&le;\frac{3\&times;{10}^8\&times;20\&times;{10}^{-12}}{2\&times;3.603\&times;0.63\&times;{10}^{-3}}-\frac{1}{2}=0.82$

所以，m=0，只有0级反射回波（主反射回波、主脉冲）存在。

同理，参见图5，利用(9)式，按1THz估算，消光系数为

${\mathrm{\&kappa;}}_g=\frac{3\&times;{10}^8}{2\mathrm{\&pi;}\&times;{10}^{12}\&times;0.63\&times;{10}^{-3}}\ln \frac{0.048}{0.033}=0.0284$

3.利用模拟退火算法求无约束条件二元函数（式(12)）的全局最小值，求得最小值处对应的实折射率n(ω)和消光系数κ(ω)，注意实折射率n(ω)和消光系数κ(ω)有一个波动范围，该波动是由于样品厚度的不确定性引起。

4.利用准确样品厚度确定法，即利用Nelder-Mead单纯形算法求上述一元函数（式(13)）的最小值，准确确定样品厚度l=0.68mm。

5.利用准确的样品厚度，再次应用模拟退火算法获得各频率点的实折射率n(ω)和消光系数κ(ω)，并进一步导出吸收系数α(ω)（式(4)）。实折射率n(ω)如图7，消光系数κ(ω)如图8，吸收系数α(ω)如图9。

6.图7、图8和图9关于GaAs晶体的折射率、消光系数和吸收系数的结果与以前报道的测量结果及理论计算值基本一致。

本发明的利用太赫兹光谱信号提取材料光学常数的方法，采用精确公式结合模拟退火优化算法和准确样品厚度确定法提取材料光学常数，可以取消弱吸收近似和材料边界约束条件，适用于不同厚度的材料，采用该方法提取的材料光学常数精确度较高。该发明在基础科学研究、食品安全、生物医学、国家安全等领域中具有重要的现实意义和经济价值。

以上内容仅为本发明的较佳实施例，对于本领域的普通技术人员，依据本发明的思想，在具体实施方式及应用范围上均会有改变之处，本说明书内容不应理解为对本发明的限制。

Claims (3)

1.一种利用太赫兹光谱信号提取材料光学常数的方法，其特征在于：包括以下步骤：步骤一，太赫兹波垂直入射到放置于空气中厚度为l的平板样品中，根据Fabry-Perot腔模型获得样品信号和参考信号太赫兹频谱的理论比值A(ω)；具体为：

$A\left(\mathrm{\&omega;}\right)=\frac{4{\tilde{n}}^2\left(\mathrm{\&omega;}\right)}{{(\tilde{n}\left(\mathrm{\&omega;}\right)+1)}^2}\exp [-j(\tilde{n}\left(\mathrm{\&omega;}\right)-1)\mathrm{\&omega;l}/c]\mathrm{FP}\left(\mathrm{\&omega;}\right)---\left(a\right)$

其中，

为复折射率，n(ω)为实折射率、κ(ω)为消光系数，c为真空中光速，FP(ω)为Fabry-Perot系数，FP(ω)的表达式为：

$\mathrm{FP}\left(\mathrm{\&omega;}\right)=\underset{p=0}{\overset{m}{\mathrm{\&Sigma;}}}{\left[\frac{(\tilde{n}\left(\mathrm{\&omega;}\right)-1)}{(\tilde{n}\left(\mathrm{\&omega;}\right)+1)}\exp (-\mathrm{j\&omega;}\tilde{n}l)\right]}^{2p}---\left(b\right)$

其中，p为整数，表示反射回波级次，p=0表示主反射回波，m表示最大反射回波级次；步骤二，利用模拟退火优化算法获得样品信号和参考信号太赫兹频谱理论比值A(ω)和实验比值A_meas(ω)之间的全局最小值，获得最小值处对应的实折射率n(ω))和消光系数κ(ω)；具体步骤如下：

（21）构建精确频谱和测量频谱之间的误差函数：

根据样品厚度l、样品信号与参考信号的时间延迟Δt、时间窗口大小t_c，获得估计折射率n_g和最大反射回波级次m值：

$n_g=\frac{\mathrm{c\&Delta;t}}{l}+n_0---\left(c\right)$

$m\&le;\frac{c{t}_c}{2n_{g}l}-\frac{1}{2}---\left(d\right)$

其中，时间窗口大小是指从参考信号入射到测量结束的时间，n₀为太赫兹波在空气中的折射率，c为太赫兹波空气中的速度；

A(ω)=T(ω)exp[-jΔφ(ω)]              （e）

同理：A_meas(ω)=T_meas(ω)exp[-jΔφ_meas(ω)]     （f）

根据样品信号峰值|E_s,max|和参考信号峰值|E_r,max|可以求出消光系数κ_g的近似表达式：

${\mathrm{\&kappa;}}_g=\frac{c}{\mathrm{\&omega;l}}\ln \frac{\left|E_{r,\max }\right|}{\left|E_{s,\max }\right|}---\left(g\right)$

结合公式（e）、（f）、（g）构建一个精确频谱和测量频谱之间的误差函数：

δ(n,κ)={ln[T(ω)]-ln[T_meas(ω)]}²+{ln[Δφ(ω)]-ln[Δφ_meas(ω)]}²      （h）

（22）利用模拟退火算法求无约束条件二元函数的全局最小值，获得δ(n,κ)的全局最小值，进一步得到最小值处对应的实折射率n(ω)和消光系数κ(ω)：

（221）输入初始变量：选取初始点x₀，使用估计折射率n_g和估计消光系数κ_g，估计最小值区间的上、下边界u和l，选定最大迭代数k_max>0，退火因子1>q>0，以及相应的函数值浮动误差容许值ε_f；

(222)求解基本变量：x=x₀，x⁰=x，f⁰=f(x⁰)；

(223)进行迭代计算，找出近似全局最小值，即从k=0到k=k_max执行以下循环块：

生成一个在[-1,1]上服从均匀分布的N×1的随机向量y，y与x大小相同，利用μ^-1定理可以求出Δx：

Δx=(u-l)·g_μ ^-1(y)， ${g_{\mathrm{\&mu;}}}^{-1}\left(y\right)=\frac{{(1+\mathrm{\&mu;})}^{\left|y\right|}-1}{\mathrm{\&mu;}}\mathrm{sign}\left(y\right),$ μ=10¹⁰⁰(k/k_max)^q；

在(l,u)区间内确定x₁=x+Δx；

如果Δf=f(x₁)-f(x)<0，或者在[0,1]上生成的均匀分布的随机数z<p=exp[-(k/k_max)^qΔf/|f(x)|/ε_f]，则令x=x₁，f(x)=f(x₁)。

如果f(x)<f⁰，则令x⁰=x，f⁰=f(x⁰)；

令k=k+1，继续循环；

(224)根据步骤（223）求出近似全局最小值x⁰后，假定x⁰为初值，进而求出准确的全局最小值；

其中，x₀，u，l，x，x⁰，x₁，Δx和g_μ(y)均为矩阵向量。

2.根据权利要求1所述的一种利用太赫兹光谱信号提取材料光学常数的方法，其特征在于：还包括步骤三，利用Nelder-Mead单纯形算法求一元函数的局部最小值，获得实折射率n(ω)和消光系数κ(ω)对样品厚度变化的最小值，确定样品厚度l，并修正实折射率n(ω)和消光系数κ(ω)；具体为：

（31）将样品从厚度范围l_min到l_max分为I等份，i表示的厚度为

构造一元误差函数：

$\mathrm{\&epsiv;}\left(l\right)=\underset{i=1}{\overset{I}{\mathrm{\&Sigma;}}}\left[\right|n(i-1)-n\left(i\right)|+|\mathrm{\&kappa;}(i-1)+\mathrm{\&kappa;}\left(i\right)\left|\right]---\left(i\right)$

其中，n(i)，κ(i)分别表示厚度为l_i的实折射率和消光系数；

（32）准确样品厚度确定法：利用Nelder-Mead单纯形算法求公式（i）的一元函数的最小值，准确确定样品厚度l，并修正各频率点的实折射率n(ω)和消光系数κ(ω)；具体求解过程如下：

(321)根据不同的厚度值l_i，按模拟退火法求出对应的实折射率n(i)和消光系数κ(i)；

(322)用Nelder-Mead单纯形算法求出一元误差函数ε(l)的最小值及其对应的厚度l；

(323)代入l修正各频率点的实折射率n(ω)和消光系数κ(ω)。

3.根据权利要求2所述的一种利用太赫兹光谱信号提取材料光学常数的方法，其特征在于：所述步骤(224)求出准确的全局最小值采用的方法为最速下降法、牛顿迭代法或Nelder-Mead单纯形算法中的任意一种。

Images