基于太赫兹光谱用于糖类分析的稀疏表示分类方法与流程

1.本发明涉及一种利用太赫兹光谱进行糖类检测的方法，具体涉及基于太赫兹光谱用于糖类分析的稀疏表示分类方法，属于太赫兹光谱技术领域。


背景技术：

2.太赫兹波(thz)是一种电磁波，频率在0.1thz到10thz之间，波长在3mm～0.03mm之间。thz波所处的位置独特，兼具微波和红外两者的一些特性，具有良好的透视性、安全性和光谱分辨率，并且具有分子指纹特性，能量低，穿透性强等优势。thz技术已经应用于许多领域，如航空航天、生物医学、安全检测，材料性能分析等；目前在农产品和食品安全领域的应用也在不断增多，如农产品和食品含水量检测、食品内部品质检测和种子识别等。
3.糖类是粮食中重要的组成成分，粮食中含有的糖类主要有葡萄糖、蔗糖和乳糖等，在粮食储藏的过程中糖类的变化会引起粮食品质发生改变。且糖类在太赫兹光谱中有明显的吸收峰，以蔗糖、果糖和乳糖作为糖类样品，对糖类进行定性分析，可以为储粮品质的检测奠定一定的基础。
4.针对太赫兹时域光谱对糖类的定性分析，文献“基于太赫兹光谱技术的d
‑
无水葡萄糖定性定量分析研究，光谱学与光谱分析，2017，37(07):2165
‑
2170”中采用多元线性回归和偏最小二乘回归的方法进行定性分析。由于thz光谱数据不仅含有线性关系，且有非线性关系，光谱维度高，存在大量的冗余信息，直接作为定性分析模型的输入会影响模型的预测精度；同时数据量过大会增加模型计算复杂度，导致计算量大、耗时长。目前已有的算法在对thz光谱数据进行处理时，往往需要大量的计算，且模型分类时间长，复杂度高，因此光谱降维成为光谱处理中不可缺少的一部分。


技术实现要素：

5.本发明的目的是：解决现有太赫兹时域光谱对糖类定性分析的算法，在处理数据时需要大量计算，且模型分类时间长、复杂度高的问题，提供基于太赫兹光谱用于糖类分析的稀疏表示分类方法，将稀疏表示引入到thz光谱领域对数据进行降维处理，使thz光谱在建模分析时更加快速准确，且数据尽可能简洁的表示出来。
6.为实现上述目的，本发明采用了以下技术方案：基于太赫兹光谱用于糖类分析的稀疏表示分类方法，包括以下步骤：
7.s1、获取糖类样品的thz吸收光谱特性：选取蔗糖、果糖和乳糖三种糖类作为样品，采用太赫兹时域光谱分析系统thz
‑
tds
‑
spectrometer先测量出一个时域光谱作为参考光谱，再对每类样品多次测量计算出平均时域光谱，并利用时域光谱计算出频域光谱，通过对比参考光谱最终计算出每类样品的吸收光谱；
8.参考光谱由thz
‑
tds系统直接测量的时域光谱信号得到，定义公式为f
r
(t)，其中t为频率，r为f(t)的小标；
9.平均时域光谱由重复测量多次每一类样品的时域光谱信号求平均值得到，定义公
式为f
s
(t)，其中t为频率，s为f(t)的小标；
10.频域光谱由平均时域光谱采用傅里叶变换得到，变换的计算公式为：其中f(ω)为频域光谱信号，f(t)为时域光谱信号；
11.吸收光谱由每类样品的频域光谱和参考光谱的频域光谱对比得到，计算公式为：其中α为吸收光谱，d为样品压片的厚度，a
r
为参考光谱的频域光谱的幅度，a
s
为糖类样品的频域光谱的幅度。
12.s2、对糖类thz吸收光谱数据进行预处理：选择快速傅氏变换fft滤波，对s1中计算出的每类样品的thz吸收光谱数据进行预处理，去除thz光谱数据中频率高于截止频率的噪声信号，得到滤波后的thz吸收光谱数据的数据集，得到每类样品的thz吸收光谱数据更加平滑的波形；
13.利用快速傅氏变换fft滤波的具体步骤为：
14.s21：对s1中计算出的每类样品的thz吸收光谱数据，计算前1％数据点和后1％数据点的平均值；
15.s22：通过s21中计算出的两个平均值点构造一条基线，并用这条基线减去每类样品的所有thz吸收光谱数据；
16.s23：对s22中得到的数据差进行fft变换；
17.s24：将s23中变换后获得的thz吸收光谱数据进行低通滤波，去除频率高于截止频率的傅里叶分量，截止频率为：其中n为指定的窗口点数，n设为15，t为两个相邻数据点之间的时间间隔；
18.s25：对滤波后的频谱进行快速傅立叶反变换ifft计算，整理后获得每类样品thz吸收光谱数据的数据集；
19.s26：将s22中构造的基线添加到s25中获取的数据集中。
20.s3、对预处理后的糖类thz吸收光谱数据进行稀疏表示，具体步骤为：
21.s31、利用预处理后的糖类thz吸收光谱数据的数据集构建稀疏表示算法模型：假设一个糖类样品的thz吸收光谱数据预处理后的数据集y表示为原始样本矩阵y∈r
m
×
n
，其中原始糖类样品吸收光谱信号y为数据集y的其中一类糖类样品的thz吸收光谱数据，则稀疏表示模型为：
[0022][0023]
其中d为过完备字典矩阵，d
i
为原子，d＝{d1，d2，...d
k
}，x为光谱稀疏表示系数且只有有限个非零元素，x＝{x1，x2，...x
n
}；
[0024]
求解得出：s.t.||x||0＜t，其中t为光谱稀疏表示系数中非零分量的数目上限；
[0025]
s32、采用正交匹配追踪omp算法求解光谱稀疏系数矩阵：
[0026]
设置稀疏度k得到迭代次数，通过输入字典矩阵d、原始糖类样品吸收光谱信号y和稀疏度k，求解得出稀疏系数x，具体算法步骤为：
[0027]
初始化：初始化字典矩阵d及稀疏度k，残差r0＝y，索引集λ＝o，t＝1；
[0028]
a1：找到残差r与字典矩阵d最匹配的原子，即找到与字典矩阵d中内积最大的某一列d
i
，其对应的下标为λ，λ
t
＝argmax
j＝1
，...n|<r
t
‑1，d
j
>|；
[0029]
a2：更新索引集λ
t
＝λ
t
‑1∪{λ
t
}和字典矩阵d中的原子
[0030]
a3：采用最小二乘法带入步骤s31中的求解公式计算，得到糖类样品吸收光谱信号的稀疏系数x
t
＝argmax||y
‑
d
t
x
t
||；
[0031]
a4：更新糖类样品吸收光谱信号的残差r
t
＝y
‑
d
t
x
t
，t＝t 1；
[0032]
a5：判断t是否大于k，若t＞k，算法停止，得到糖类样品吸收光谱的稀疏系数矩阵；若t≤t 1，则重新跳转到步骤a1执行算法。
[0033]
s33、采用k
‑
svd算法更新字典矩阵：
[0034]
通过输入字典矩阵d，原始糖类样品吸收光谱信号y和步骤s32中计算出的稀疏系数矩阵，得到更新后的字典矩阵d，具体算法步骤为：
[0035]
初始化：随机在原始糖类样品吸收光谱信号y中取k个列向量或者从原始样本矩阵y的左奇异矩阵中取前k个向量作为初始化字典的原子，初始化字典d0∈r
m
×
k
，j＝0，得到初始字典d
j
；
[0036]
b1：利用得到的字典d
j
进行光谱稀疏编码，求得糖类样品吸收光谱的稀疏系数矩阵x
j
∈r
k
×
n
；
[0037]
b2：逐列更新字典矩阵d
j
，字典的列d
k
∈{d1，d2，...，d
k
}，计算误差矩阵
[0038]
b3：取糖类样品吸收光谱的稀疏矩阵x
j
第k个不为零的行向量的索引的集合
[0039]
b4：从e
k
中取出对应索引ω
k
不为零的列，得到e
′
k
；
[0040]
b5：对e
′
k
作奇异值分解e
k
＝uδv
t
，取u的第1列更新字典的第k列，d
k
＝u(
·
，1)，令将更新到原处，j＝j 1；
[0041]
b6：判断是否达到指定的迭代步数和容许误差，若是则停止算法，得到三种糖类样品吸收光谱信号的稀疏系数矩阵x；若为否，则继续进入omp算法进行循环计算直至满足迭代条件为止。
[0042]
本发明的有益效果是：
[0043]
1)本发明方法以糖类为研究对象，针对thz光谱数据量大、维数高、结构复杂及模型计算时间长等问题，采用稀疏表示对预处理后的thz光谱数据进行降维约简；使thz光谱在建模分析时更加快速准确，且数据尽可能简洁的表示出来。
[0044]
2)本发明方法中选用omp算法计算稀疏系数，k
‑
svd更新字典，实现太赫兹光谱数据的稀疏表示，并与主成分分析、线性判别分析、局部线性嵌入和等度量映射传统降维方法对比，建立svm分类模型；相对于其它降维方法，构建的太赫兹稀疏表示模型的运算时间大大缩减，实验结果验证了稀疏表示非线性降维方法提高模型效率的有效性，提升了太赫兹光谱数据质量和模型检测精度，丰富和发展了太赫兹检测理论和方法。
附图说明
[0045]
图1为三种糖类样品的thz吸收光谱图；
[0046]
图2为三种糖类样品的thz吸收光谱进行fft滤波处理后的平滑效果图；
[0047]
图3为稀疏表示的算法图；
[0048]
图4为thz吸收光谱稀疏表示的算法流程图；
[0049]
图5为稀疏表示与常见降维方法的对比结果表；
[0050]
图6为thz光谱稀疏表示模型与常见降维方式模型的运算时间柱状图。
具体实施方式
[0051]
下面结合附图和具体实施例对本发明作进一步的解释说明。
[0052]
实施例：如图1
‑
6所示，本发明所述的基于太赫兹光谱用于糖类分析的稀疏表示分类方法，包括以下步骤：
[0053]
s1、获取糖类样品的thz吸收光谱特性：
[0054]
选取蔗糖、果糖和乳糖三种糖类作为样品，采用太赫兹时域光谱分析系统thz
‑
tds
‑
spectrometer先测量出一个时域光谱作为参考光谱，再对每类样品多次测量计算出平均时域光谱，并利用时域光谱计算出频域光谱，通过对比参考光谱最终计算出每类样品的吸收光谱；如图1给出了不同糖类样品的吸收光谱图。
[0055]
在使用时域光谱仪进行测量时域光谱数据时会对同一个样品连续测量多次再求平均得到平均时域光谱，可以减少实验过程中的系统误差，在实际实验过程一般测量10次求平均。
[0056]
参考光谱由thz
‑
tds系统直接测量的时域光谱信号得到，定义公式为f
r
(t)，其中t为频率，r为f(t)的小标；
[0057]
平均时域光谱由重复测量多次每一类样品的时域光谱信号求平均值得到，定义公式为f
s
(t)，其中t为频率，s为f(t)的小标；
[0058]
频域光谱由平均时域光谱采用傅里叶变换得到，变换的计算公式为：其中f(ω)为频域光谱信号，f(t)为时域光谱信号；
[0059]
吸收光谱由每类样品的频域光谱和参考光谱的频域光谱对比得到，计算公式为：其中α为吸收光谱，d为样品压片的厚度，a
r
为参考光谱的频域光谱的幅度，a
s
为糖类样品的频域光谱的幅度。
[0060]
s2、对糖类thz吸收光谱数据进行预处理：
[0061]
选择快速傅氏变换fft滤波，对s1中计算出的每类样品的thz吸收光谱数据进行预处理，去除thz光谱数据中频率高于截止频率的噪声信号，得到滤波后的thz吸收光谱数据的数据集，使得到每类样品的thz吸收光谱数据的波形更加平滑；处理结果如图2所示。
[0062]
利用快速傅氏变换fft滤波的具体步骤为：
[0063]
s21：对s1中计算出的每类样品的thz吸收光谱数据，计算前1％数据点和后1％数据点的平均值；
[0064]
s22：通过s21中计算出的两个平均值点构造一条基线，并用这条基线减去每类样
品的所有thz吸收光谱数据；
[0065]
s23：对s22中得到的数据差进行fft变换；
[0066]
s24：将s23中变换后获得的thz吸收光谱数据进行低通滤波，去除频率高于截止频率的傅里叶分量，截止频率为：其中n为指定的窗口点数，n设为15，t为两个相邻数据点之间的时间间隔；
[0067]
s25：对滤波后的频谱进行快速傅立叶反变换ifft计算，整理后获得每类样品thz吸收光谱数据的数据集；
[0068]
s26：将s22中构造的基线添加到s25中获取的数据集中。
[0069]
s3、对预处理后的糖类thz吸收光谱数据进行稀疏表示，流程图如图4所示，具体步骤为：
[0070]
s31、利用预处理后的糖类thz吸收光谱数据的数据集构建稀疏表示算法模型：
[0071]
假设一个糖类样品的thz吸收光谱数据预处理后的数据集y表示为原始样本矩阵y∈r
m
×
n
，其中原始糖类样品吸收光谱信号y为数据集y的其中一类糖类样品的thz吸收光谱数据，则稀疏表示模型为：
[0072][0073]
其中d为过完备字典矩阵，d
i
为原子，d＝{d1，d2，...d
k
}，x为光谱稀疏表示系数且只有有限个非零元素，x＝{x1，x2，...x
n
}；
[0074]
求解得出：s.t.||x||0＜t，其中t为光谱稀疏表示系数中非零分量的数目上限；稀疏表示的算法图如图3所示。
[0075]
s32、采用正交匹配追踪omp算法求解光谱稀疏系数矩阵：
[0076]
设置稀疏度k得到迭代次数，通过输入字典矩阵d、原始糖类样品吸收光谱信号y和稀疏度k，求解得出稀疏系数x，具体算法步骤为：
[0077]
初始化：初始化字典矩阵d及稀疏度k，残差r0＝y，索引集λ＝o，t＝1；
[0078]
a1：找到残差r与字典矩阵d最匹配的原子，即找到与字典矩阵d中内积最大的某一列d
i
，其对应的下标为λ，λ
t
＝argmax
j＝1
，...n|<r
t
‑1，d
j
>|；
[0079]
a2：更新索引集λ
t
＝λ
t
‑1∪{λ
t
}和字典矩阵d中的原子
[0080]
a3：采用最小二乘法带入步骤s31中的求解公式计算，得到糖类样品吸收光谱信号的稀疏系数x
t
＝argmax||y
‑
d
t
x
t
||；
[0081]
a4：更新糖类样品吸收光谱信号的残差r
t
＝y
‑
d
t
x
t
，t＝t 1；
[0082]
a5：判断t是否大于k，若t＞k，算法停止，得到糖类样品吸收光谱的稀疏系数矩阵；若t≤t 1，则重新跳转到步骤a1执行算法。
[0083]
s33、采用k
‑
svd算法更新字典矩阵：
[0084]
通过输入字典矩阵d，原始糖类样品吸收光谱信号y和步骤s32中计算出的稀疏系数矩阵，得到更新后的字典矩阵d，具体算法步骤为：
[0085]
初始化：随机在原始糖类样品吸收光谱信号y中取k个列向量或者从原始样本矩阵y的左奇异矩阵中取前k个向量作为初始化字典的原子，初始化字典d0∈r
m
×
k
，j＝0，得到初
始字典d
j
；
[0086]
b1：利用得到的字典d
j
进行光谱稀疏编码，求得糖类样品吸收光谱的稀疏系数矩阵x
j
∈r
k
×
n
；
[0087]
b2：逐列更新字典矩阵d
j
，字典的列d
k
∈{d1，d2，...，d
k
}，计算误差矩阵
[0088]
b3：取糖类样品吸收光谱的稀疏矩阵x
j
第k个不为零的行向量的索引的集合
[0089]
b4：从e
k
中取出对应索引ω
k
不为零的列，得到e
′
k
；
[0090]
b5：对e
′
k
作奇异值分解e
k
＝uδv
t
，取u的第1列更新字典的第k列，d
k
＝u(
·
，1)，令将更新到原处，j＝j 1；
[0091]
b6：判断是否达到指定的迭代步数和容许误差，若是则停止算法，得到三种糖类样品吸收光谱信号的稀疏系数矩阵x；若为否，则继续进入omp算法进行循环计算直至满足迭代条件为止。
[0092]
在k
‑
svd算法中代入omp算法计算出的稀疏系数矩阵，并利用稀疏系数矩阵进行字典更新得到更新字典矩阵，并将更新后的字典矩阵d重新代入正交匹配追踪omp算法计算稀疏系数矩阵，直至达到设定的迭代次数停止计算，最终得到最后一次计算出的稀疏系数矩阵x。
[0093]
以蔗糖、果糖和乳糖三种样品预处理后的吸收光谱数据为原始矩阵样本，初始化字典时将原始矩阵样本中的部分太赫兹光谱数据作为初始字典，然后设置迭代次数，其中最大迭代次数为300次，采用omp算法求解thz吸收光谱的稀疏系数矩阵，采用k
‑
svd算法更新字典矩阵，计算误差矩阵；当达到设置的迭代条件时训练结束，此时可以得到三种糖类样品吸收光谱信号的稀疏系数矩阵。
[0094]
为了验证稀疏表示在thz光谱数据处理上的高效性，将其与传统的降维方式主成分分析(pca)、线性判别分析(lda)、局部线性嵌入(lle)和等度量映射(isomap)等进行对比，选取了蔗糖、果糖和乳糖样品预处理后的吸收光谱为研究对象，以其稀疏表示后计算得到的稀疏系数矩阵作为模型的输入，构建支持向量机(support vector machine，svm)分类模型，对比结果如图5所示。
[0095]
从图6中可以看到thz光谱稀疏表示后作为svm分类模型的输入相比其他三种降维方式模型的运算时间大大缩减；因此选用稀疏表示作为thz光谱数据的降维方法可以有效的降低模型的运算时间。
[0096]
本发明方法针对thz光谱数据量大、维数高、结构复杂及模型计算时间长等问题，采用稀疏表示对预处理后的thz光谱数据进行降维约简；使thz光谱在建模分析时更加快速准确，且数据尽可能简洁的表示出来。
[0097]
选用omp算法计算稀疏系数，k
‑
svd更新字典，实现太赫兹吸收光谱数据的稀疏表示，并与主成分分析、线性判别分析、局部线性嵌入和等度量映射四种传统降维方法对比，建立svm分类模型。相对于其它降维方法，构建的太赫兹稀疏表示模型耗时缩减为0.132s，实验结果验证了稀疏表示非线性降维方法提高模型效率的有效性，提升了太赫兹光谱数据
质量和模型检测精度，丰富和发展了太赫兹检测理论和方法。
[0098]
以上所述仅用以说明本发明的技术方案而非限制，本领域普通技术人员对本发明的技术方案所做的其他修改或者等同替换，只要不脱离本发明技术方案的精神和范围，均应涵盖在本发明的权利要求范围当中。