一种水蒸气太赫兹吸收谱的测量方法与流程

本发明涉及一种用于太赫兹吸收谱的测量方法，尤其是一种水蒸气太赫兹吸收谱的测量方法。
背景技术：
：太赫兹(terahertz,thz)波是指频率从0.1thz～10thz的电磁辐射。从频率上看，它介于无线电波毫米波和光波红外线之间的频谱过渡区；能量上介于电子和光子之间。由于其独特的物理特性，使得thz波在成像、无标记的基因检测、生物、物理和化学领域，以及宽带通信等场合都有广泛的应用。而近年来，飞秒激光技术的发展和成熟为thz波的应用和研究提供了有效手段。由于飞秒激光技术受热噪声影响小，而且可以在飞秒时间窗内发射/接收探测信号，因此，thz时域光谱技术(thztime-domainspectroscopy,thz-tds)得到了长足发展。thz-tds本质上属于相位相干电磁辐射探测技术，它能直接检测thz波与被测样品分子的相互作用，同时获得thz辐射的振幅和相位，因此可精确测量物质的传输和吸收特性。这种光谱技术，在气体、液体、药品以及生物分子的谱特性研究中显示出鲜明的优势。由于大气中广泛存在水蒸气，研究水蒸气对thz辐射的吸收特性是诸多样品测试和研究的基础。近年来，国内外学者针对水蒸气的thz吸收特性的研究做了大量工作，且大多研究均基于thz-tds技术，如现有研究比较细致地提取了水蒸气在thz波段的吸收谱“线”，但没有全面地揭示thz-tds技术测量水蒸气吸收谱的时间-频率分布特性。现有研究将小波变换引入到thz-tds信号的吸收谱处理中，取得了较好的时间-频率分布特性呈现效果。由于小波变换具有良好的时频分辨率，理论上采用逐渐细化的频域步长，可以分析任意频率细节，但实际应用中，有两方面原因将影响小波变换的时频分析效果：其一，若要分析任意频率细节，则需要非常精细的频域步长，这将带来不可预估的计算量；其二，由于水蒸气的thz-tds测量信号具有十分典型的非平稳特性，即频率的跳变具有突然性，即便是善于边沿检测的小波也很难捕捉这些突变，从而导致时频特性的模糊度增加。因此有必要设计出一种水蒸气太赫兹吸收谱的测量方法，能够在计算量较小的情况下，实现水蒸气的thz吸收谱时频特性分析，进而估计出各吸收峰对应的频率。自适应短时fourier变换是短时fourier变换的优化版本，它通过对样本的处理和各种预设条件，自适应地调整“短时”时间窗，从而更好地捕捉时频特性的细节信息。但是每种自适应短时fourier变换算法或不同场合所采用的自适应短时fourier变换都各不相同，即“自适应”的判断标准和方法因信号特性而各不相同。技术实现要素：本发明的目的在于：提供一种水蒸气太赫兹吸收谱的测量方法，能够在计算量较小的情况下，实现水蒸气的thz吸收谱时频特性分析，进而估计出各吸收峰对应的频率。为了实现上述发明目的，本发明公开了一种水蒸气太赫兹吸收谱的测量方法，包括如下步骤：步骤1，对采样频率为fs的水蒸气thz-tds信号进行自适应短时fourier变换，将每次由自适应调整后的fft结果保存至二维数组φ(k,j)中；步骤2，根据二维数组φ(k,j)利用频率修正法计算吸收峰对应的频率值。进一步地，步骤1中，对采样频率为fs的水蒸气thz-tds信号进行自适应短时fourier变换的具体步骤为：步骤1.1，在采样频率为fs的水蒸气thz-tds信号上从起始位置开始利用gauss窗依次选取各个信号样本片段；步骤1.2，对gauss窗内选取的信号样本片段进行快速fourier变换，检测fft包络中水蒸气的吸收峰个数；步骤1.3，根据吸收峰个数自适应调整gauss窗的宽度，使fft包络中有且仅有一个吸收峰；步骤1.4，滑动gauss窗至下一信号样本片段，gauss窗的宽度为上一自适应调整后的宽度，再重复步骤1.2和1.3直至gauss窗滑动至水蒸气thz-tds信号末端，且最后一个gauss窗内的fft包络中有且仅有一个吸收峰。进一步地，步骤1.1中，gauss窗的宽度小于等于水蒸气thz-tds信号长度的十分之一。进一步地，步骤1.2中，在检测fft包络中水蒸气的吸收峰个数时，检测的频率范围为0.5thz～2.4thz。进一步地，步骤1.2中，对gauss窗内选取的信号样本片段进行快速fourier变换时，fft点数固定设置为n为整个水蒸气thz-tds信号采样个数，且当n小于3072时，fft点数固定设置为1024。进一步地，步骤1.3中，根据吸收峰个数自适应调整gauss窗的宽度的调整策略为：若吸收峰的个数大于1，则逐步减小gauss窗的宽度并重复步骤1.2，直到fft包络中有且仅有一个吸收峰；若吸收峰的个数小于1，则逐步增大gauss窗的宽度并重复步骤1.2，直到fft包络中有且仅有一个吸收峰；若吸收峰的个数等于1，则无需调整gauss的窗宽度。进一步地，在逐步增大gauss窗的宽度并重复步骤1.2时，若gauss窗的宽度大于64，则先逐步调整gauss窗的方差并重复步骤1.2，即调整gauss窗的方差每次减小0.1，调整次数限定为5次，若还未检测到有效的吸收峰，则再增加gauss窗的宽度，再重复上述调整过程逐步减小gauss窗的方差，直到fft包络中有且仅有一个吸收峰。进一步地，在逐步减小gauss窗的宽度并重复步骤1.2时，若gauss窗的宽度小于16，则先逐步调整gauss窗的方差并重复步骤1.2，即调整gauss窗的方差每次增加0.1，调整次数限定为5次，若吸收峰的个数依然大于1，则再减小gauss窗的宽度，再重复上述调整过程逐步增加gauss窗的方差，直到fft包络中有且仅有一个吸收峰。进一步地，步骤1中，二维数组φ(k,j)为k行j列的二维数组，其中，k为每次fft的点数，j为gauss窗的滑动次数。进一步地，步骤2，利用频率修正法计算吸收峰对应的频率值的具体步骤为：步骤2.1，依次取出二维数组φ(k,j)的每一列φi(k)，其中i＝1,2,...,j；步骤2.2，搜索φi(k)中最小值所对应的位置索引号，记为ηi；步骤2.3，取φi(η-2)、φi(η-1)、φi(η)、φi(η+1)以及φi(η+2)五个元素，计算吸收峰对应位置为：步骤2.4，根据下式计算水蒸气thz-tds信号吸收峰对应的频率值为：式中，fs为水蒸气thz-tds信号的采样频率，fs/k表示为fft的频率分辨率，k为每次fft的点数，j为gauss窗的滑动次数。本发明的有益效果在于：本发明对水蒸气thz-tds信号做自适应短时fourier变换，依据频域吸收峰的个数进行自适应调整gauss窗宽度，一方面保证了所有包含吸收峰的频谱细节都能被呈现，另一方面若gauss窗内信号样本没有吸收峰，将自适应地增大gauss窗的宽度，提高算法处理剩余信号样本的效率，降低运算量；采用频率修正法估计吸收峰频率，克服了离散fourier变换带来的频域“栅栏效应”，提高了频率估计精度。附图说明图1为本发明的方法流程图；图2为本发明的实例样本的整体频谱图；图3为本发明的自适应短时fourier变换获得的室温条件下水蒸气的thz-tds信号时间-频率分布图；图4为现有技术中的小波变换获得的室温条件下水蒸气的thz-tds信号时间-频率分布图；图5为现有技术中的小波变换在0～50皮秒区间呈现出水蒸气thz-tds信号的时间-频率分布特征；图6为本发明的自适应短时fourier变换在0～50皮秒区间呈现出水蒸气thz-tds信号的时间-频率分布特征。具体实施方式下面结合附图对本发明的技术方案进行详细说明，所述的实施方式是示例性的，仅用于解释本发明，而不能理解为对本发明的限制。如图1所示，本发明公开了一种水蒸气太赫兹吸收谱的测量方法，包括如下步骤：步骤1，对采样频率为fs的水蒸气thz-tds信号进行自适应短时fourier变换，将每次由自适应调整后的fft结果保存至二维数组φ(k,j)中，其中，水蒸气thz-tds信号的样本为x(τ),τ＝0...n-1；步骤2，根据二维数组φ(k,j)利用频率修正法计算吸收峰(即由于频率成份的缺失而形成的“谷”)对应的频率值。进一步地，步骤1中，对采样频率为fs的水蒸气thz-tds信号进行自适应短时fourier变换的具体步骤为：步骤1.1，在采样频率为fs的水蒸气thz-tds信号上从起始位置开始利用gauss窗依次选取各个信号样本片段，信号样本片段为x(τ′)＝x(τ).*gauss(λ),τ＝0...λ-1，λ为gauss窗的宽度；步骤1.2，对gauss窗内选取的信号样本片段进行快速fourier变换，检测fft包络中水蒸气的吸收峰个数；步骤1.3，根据吸收峰个数自适应调整gauss窗的宽度，使fft包络中有且仅有一个吸收峰；步骤1.4，滑动gauss窗至下一信号样本片段，即相邻两个gauss窗内的样本只有一个不重合，gauss窗每次只滑动一个样本的步长，gauss窗的宽度为上一自适应调整后的宽度，再重复步骤1.2和1.3直至gauss窗滑动至水蒸气thz-tds信号末端，且最后一个gauss窗内的fft包络中有且仅有一个吸收峰，由于thz-tds信号的最后一个样本已经被囊括其中，无法再根据fft包络的吸收峰个数增加gauss窗的宽度，因此无需进行自适应调整。进一步地，步骤1.1中，gauss窗的宽度小于等于水蒸气thz-tds信号长度的十分之一。gauss窗的宽度会影响频率分辨率，因此不能超过水蒸气thz-tds信号长度的十分之一。进一步地，步骤1.2中，在检测fft包络中水蒸气的吸收峰个数时，检测的频率范围为0.5thz～2.4thz，确定fft各谱线对应频率的方法是：若fft运算点数为k，则fft的频率分辨率为fs/k，各谱线对应频率为(n-1)·fs/k，其中n是fft谱线的索引号。水蒸气的吸收峰主要集中在0.5thz～2.4thz频率范围内，因此，在检测fft包络中的吸收峰时，应针对该频率范围，而不是全频域检测，能够提高运行效率。进一步地，步骤1.2中，对gauss窗内选取的信号样本片段进行快速fourier变换时，，fft点数固定设置为n为整个水蒸气thz-tds信号采样个数，且当n小于3072时，fft点数固定设置为1024。本实施例优选为2048个，能够有效提高频率分辨率，且每个gauss窗内的信号样本片段均做相同点数的fft运算。进一步地，步骤1.3中，根据吸收峰个数自适应调整gauss窗的宽度的调整策略为：若吸收峰的个数大于1，则说明窗宽度太大，窗内信号样本片段的频率分量过于丰富，此时需要逐步减小gauss窗的宽度并重复步骤1.2，直到fft包络中有且仅有一个吸收峰；若吸收峰的个数小于1，则说明窗内信号样本片段的频率分量过于单一，需要逐步增大gauss窗的宽度并重复步骤1.2，直到fft包络中有且仅有一个吸收峰；若吸收峰的个数等于1，则无需调整gauss窗的宽度。在进行gauss窗的宽度调整时，进一步地，在逐步增大gauss窗的宽度并重复步骤1.2时，若gauss窗的宽度大于64，则先逐步调整gauss窗的方差并重复步骤1.2，即调整gauss窗的方差每次减小0.1，调整次数限定为5次，若还未检测到有效的吸收峰，则再增加gauss窗的宽度，再重复上述调整过程逐步减小gauss窗的方差，直到fft包络中有且仅有一个吸收峰。即先逐步减小gauss窗的方差，再调整gauss窗的宽度，能够防止gauss窗的宽度过大而影响频率分辨率。在进行gauss窗的宽度调整时，进一步地，在逐步减小gauss窗的宽度并重复步骤1.2时，若gauss窗的宽度小于16，则先逐步调整gauss窗的方差并重复步骤1.2，即调整gauss窗的方差每次增加0.1，调整次数限定为5次，若吸收峰的个数依然大于1，则再减小gauss窗的宽度，再重复上述调整过程逐步增加gauss窗的方差，直到fft包络中有且仅有一个吸收峰。即先逐步增加gauss窗的方差，再调整gauss窗的宽度，能够防止gauss窗的宽度过小而影响时间分辨率。进一步地，步骤1中，二维数组φ(k,j)为k行j列的二维数组，其中，k为每次fft的点数，j为gauss窗的滑动次数，若以变量n·fs/k为x轴(n＝0,1,...,k-1)，m/fs为y轴(m＝0,1,...,j-1)，二维数组φ(k,j)中元素的幅度为z轴，画出二维数组φ(k,j)对应的三维图，并投影到x-y平面，便可得到水蒸气thz-tds信号吸收谱的时频分布图。进一步地，步骤2，利用频率修正法计算吸收峰对应的频率值的具体步骤为：步骤2.1，依次取出二维数组φ(k,j)的每一列φi(k)，其中i＝1,2,...,j；步骤2.2，搜索φi(k)中最小值所对应的位置索引号，记为ηi；步骤2.3，取φi(η-2)、φi(η-1)、φi(η)、φi(η+1)以及φi(η+2)五个元素，计算吸收峰对应位置为：步骤2.4，根据下式计算水蒸气thz-tds信号吸收峰对应的频率值为：式中，fs为水蒸气thz-tds信号的采样频率，fs/k表示为fft的频率分辨率，k为每次fft的点数，j为gauss窗的滑动次数。在本发明的实例中，水蒸气thz-tds信号样本为采用menlosystemsterak15型thz-tds仪在室温环境20.5℃，湿度为64％时测得，采样频率fs为30thz，信号采样个数n为6256，为降低加性高斯白噪声的影响，信号采集系统已经设置为采集1000次后取平均的结果。信号样本的整体频谱如图2所示。从图2可已清晰地看到，在0.5thz～2.4thz频率范围内存在多个吸收峰，为准确测得这些吸收峰，首先对实例中的信号样本进行自适应短时fourier变换，设置初始gauss窗的宽度为27，fft变换点数统一为2048点；通过自适应短时fourier变换获得的水蒸气thz-tds信号时频分布图如图3所示。为了进行更直观的对比，本发明还采用了现有技术中的小波变换方法对该实例中的信号样本进行了处理，获得水蒸气thz-tds信号时频分布图如图4所示。现有技术中的小波变换方法出处为：“邓玉强,邢岐荣,郎利影.thz波的小波变换频谱分析,物理学报,2005,54(11),pp:5224-5227”。对比图3和图4，可以看出自适应短时fourier变换和小波变换都能够准确地呈现出水蒸气thz-tds信号的时间-频率分布特征，但在细节方面，特别是在时间范围为0～50皮秒区间内，本发明的自适应短时fourier变换比小波变换呈现了更多、更细致的时频特性。如图5和图6所示，图5是0～50皮秒区间小波变换呈现出水蒸气thz-tds信号的时间-频率分布特征，图6是0～50皮秒区间本发明自适应短时fourier变换呈现出水蒸气thz-tds信号的时间-频率分布特征，结合图5和图6可以看出时间范围0～50皮秒、频率0.5thz～2.4thz区间的浅色部分更有区分度了。本发明在自适应短时fourier变换的基础上，采用基于五个元素求重心的频率修正方法，估计各吸收峰对应的频率，并与现有技术中测得的水蒸气的thz吸收峰频率进行了对比，对比结果如表1所示。现有技术中测得的水蒸气的thz吸收峰频率的内容出处为：“x.xin,h.altan,a.saint,et.al.terahertzabsorptionspectrumofparaandorthowatervaporsatdifferenthumiditiesatroomtemperature.journalofappliedphysics.2006,100,https://doi.org/10.1063/1.2357412”。从表1可以看出，本发明方法测得的水蒸气典型thz吸收峰频率与现有技术基本一致，考虑到实验中信号采集环境的温度、湿度的些许差别，表1中所测的频率存在微小差别是意料之中的；需要说明的是，在表1中，吸收峰序号为6的位置，本发明方法检测到两个有效的thz吸收峰频率分别为1.2305thz和1.2891thz，而不是现有技术所揭示的只有一个有效吸收峰。本发明方法检测到的两个有效的thz吸收峰这事实上与图2是吻合的，从图2也能大致确定，在大约1.2thz～1.3thz之间，的确存在两个有效的吸收峰，这也从侧面再一次验证了本发明所展示的方法的有效性。表1现有技术方法与本发明方法对比图吸收峰序号现有技术测得的频率(thz)本发明方法测得的频率(thz)10.5570.559420.7530.752730.9890.990941.0981.094751.1641.171961.22641.2305和1.289171.41161.410681.6031.606191.671071.6792101.7181.7178111.79841.8064121.8691.8701131.92071.9236142.07592.0708152.26442.2852162.34592.3438172.39322.3952当前第1页12