【《农产品产地溯源的近红外光谱数据的去噪案例分析》7900字（论文）】

农产品产地溯源的近红外光谱数据的去噪案例分析目录TOC\o"1-3"\h\u23224农产品产地溯源的近红外光谱数据的去噪案例分析 186651.1引言 1323951.2光谱数据的预处理方法 241301.3基于EMD算法的光谱去噪方法 11272731.3.1EMD算法介绍 11233591.3.2光谱数据的去噪处理 14166541.3.3EMD算法的去噪性能评价 211871.4去噪结果与分析 24引言水稻中含有各种营养素，营养素中又含有许多与氢基有关的基团，这些基团会吸收包含在近红外光波中的同频率光波能量，产生与大米各种成分含量有关的近红外光谱图，大米的各成分含量与大米的产地相互关联，判断出不同结构的光谱图即可判断各个大米的产地类型，这就实现了大米产地溯源的定量分析。基于近红外光谱的大米产地溯源检测的原理即为：近红外光被氢基团吸收后，近红外光衰减后被探测器接收，而不同区域的水稻，对应着不同种类和数量的氢基团，它们在近红外光波段的吸收能力也各不相同，从而根据光谱数据建立起数学模型。在建立的数学模型的基础上，利用傅里叶变换近红外光谱仪采集待测产地未知水稻样品的近红外光谱数据，并将其代入上述已经建立好的数学模型，最后得到具体的结果。然而，水稻近红外波段比较复杂，混叠严重。此外，由于水稻自身、傅里叶变换光谱仪和实验测量环境的影响，采集到的近红外光谱信号存在不同程度的噪声，吸收的光谱信号强度较弱。因此，在建立水稻产地溯源模型之前，需要对光谱数据进行去噪，借此提高信号的信噪比同时，还可以获得能够真实反映水稻成分含量的有用信息。本章节的主要研究内容是大米近红外光谱数据的去噪处理。首先以大米光谱数据为背景，对近红外大米光谱数据所涉及到的预处理算法原理、使用方法进行阐述；其次由于大米近红外光谱数据存在噪声，进而导致大米产地溯源模型检测精度低的问题，提出基于EMD算法的大米近红外光谱数据去噪方法，该方法通过EMD算法对不同大米的近红外光谱数据进行IMF的自适应分类，并通过分析采用恰当的区分方式对包含有有用信号的相关模态进行合并重构；最后通过对比试验等分析EMD算法的有效性。光谱数据的预处理方法傅里叶变换红外光谱的数据主要是由被测大米的红外光谱特征主导的，同时采集又受到外界环境的影响较大，众多不确定因素造就了采集得到的近红外光谱信号的复杂性。在近红外光谱数据中，除了样品的自身有用信息外，不可避免会含有其它无关噪声，这些噪声的存在会影响数据的分析和模型的建立。为了消除各种噪声对光谱信号分析的影响，需要通过适当的预处理手段对光谱数据进行降噪，以便减少各种后续操作的错误率，提高光谱的信噪比，进而提升所建立近红外光谱模型的准确性和预测精度，改善光谱分析信号的性能。但光谱数据去噪只是光谱数据预处理中最基本的方法之一。近红外光谱信号的常用的预处理方法主要还有以下几种：1、主线校正主线校正顾名思义即通过设置不同的主线使被测得到的光谱曲线消除偏移的方法，可以利用全光谱的数据减去其中的最低值，即所谓借助常偏移量进行消除。还可以利用全光谱数据按照最小二乘法原理拟合一条直线，然后用全部数据减去该直线，进而完成主线校正。该类型方法通常已经封装在光谱仪内部，得到的结果也已进行过相应处理，故该类型方法通常不体现在外部数据处理的MATLAB代码中。图3-1原始数据光谱图2、坐标图转换由于光谱采集软件默认使用透过率作为纵坐标，在计算时为了与文献保持一致，需要先将透过率转化为吸光度，再进行其他光谱前处理。MATLAB语句是：x=2-log10(x)。图3-2透光率转吸光度后光谱图限定范围使用的是MATLAB内自带的函数mapminmax，这是一种会规定函数上限与下限的函数，在训练模型的时候，能够避免因数据值过大或过小造成的误差，所以该方法使用非常广泛，同时MATALB自带该算法语句，侧面也说明了使用的频繁程度。图3-3限定范围后光谱图3、数据增强型处理采用偏最小二乘法和支持向量机建立近红外光谱分类模型时，光谱数据的变化程度会与待测大米产地或组成成分的波动进行互联。因此，该类方法的主要目的是扩大不同产地大米近红外光谱数据之间的差异，同时还可以削减无关信号的影响作用，进而提高后续模型的预测精度。常用的数据增强型处理方法有归一化、标准化、平均化以及正态变量交化等。归一化的原理是用向量自身去除以向量模长，即使原向量同比例缩小到1个单位长度；标准化与归一化类同，不同的是标准化目的是使原向量的均值变为0；而正态变量交化的目的是为了让数据符合正态分布的形式。他们的本质意义都是要让纵坐标（透过率）幅度上的微小差异不会被单次测量之间的误差所掩盖。平均化同样也是一种使用广泛的预处理方法，即将多个样品数据的值求平均值。这种预处理手段能够有效避免因操作失误或者是环境波动所引起的个别测量误差，但这样做的缺点是会导致训练数据数量的下降。平均化以及标准化这两种方法在数据预处理领域使用频繁，值得注意的是，在对近红外光谱数据进行预处理的同时，这两种方法也会对被测数据的性质或组成进行同样的变换。正态变量交化是在1998年首次由SvanteWold等人提出，该方法原理是利用正交化方法去除原始数据中与测量无关的数据，是目前使用最广的光谱数据预处理方法之一。图3-4平均化（上左）、标准化（上右）、归一化（下左）、正态变量交化（下右）后的光谱图4、平滑处理由于近红外光谱采集仪器设备本身的问题，或是周围环境的波动等原因，会导致光谱信号中含有不同程度的噪声。噪声的存在会影响后续模型的预测精度，所以需要进行去噪处理，信号平滑就是常见的去噪手段，最终达到提高信号信噪比的目标。滑动平均平滑法中主要的是如何选择一个合适的窗口宽度，常用的有五点平滑法，五点平滑法的理论依据是利用五个连续的数据点作为一个滑动窗口，将窗口包含的五个数据进行特殊处理来获得一个新的数据点，五点平滑的算法是可以根据自我需求进行自行修改的。需要注意的是，由于原信号的前2个点和后2个点处，不足以凑够5个点，所以需要使用特殊且规定好的平滑规则；而从第三个点到倒数第三个点，使用的则是另外一种算法。例如，采用如下处理方式，其中b(k)为各数据点，其前两个数据点和中间各点以及后两个数据点的规定算法如下图3-5所列：图3-5五点平滑算法内部程序语句下面分别用10次、100次来平滑数据，得出结果如图所示，蓝色为原始数据、红色为平滑10次、绿色为平滑100次，观察可以发现，数据随处理次数增加变得更加平滑。图3-6五点平滑后光谱图但是，滑动平均平滑法正是由于其本质只是简单的加权平均算法这样的处理原理，会导致以下一些问题：如果滑动窗口包含的数据太少，就会使得平滑效果较差，去噪不够明显；如果滑动窗口包含的数据太多，就会导致处理过度，缺失有用信息，导致失真。Savitsky-Golay卷积平滑法的本质也是一种利用多数据点平滑的方法，这一点与滑动平均平滑法一致，不同的是，前者需要通过对滑动窗口内的数据，进行多项式拟合平滑，但后者只是简单的求均值，所以SG卷积平滑法更看重中心点的作用。Savitsky-Golay卷积平滑算法主要涉及到两个参数：窗口宽度和多项式次数。不同的移动窗口宽度及多项式次数会对平滑产生不同效果，如下图所示，利用SG平滑法对采集得到的数据进行处理，其中蓝色为原始的近红外光谱图，黑色为窗口宽度20、多项式次数2的处理后光谱数据，绿色为窗口宽度20、多项式次数5的处理后光谱数据，黄色为窗口宽度50、多项式次数2的处理后光谱数据，红色为窗口宽度50、多项式次数5的处理后光谱数据。图3-7Savitsky-Golay卷积平滑后光谱图由上图分析可得：多项式次数一定时，窗口宽度越大，所包含的数据越多，平滑起来的效果也越好，但窗口的选择也不能过大，过大会使得包含的数据过多，造成原本没必要消减的信息也被平滑掉；窗口宽度一定时，多项式次数越多，处理后的数据越贴合原始数据，对噪声的抵抗能力越弱，但多项式次数也不宜过大，否则容易发生过拟合。5、小波变换小波变换（Wavelettransform,WT）是被较少使用到的一种预处理手段，它同时也是一种高效的去噪处理方法，其去噪效果远比前述的平滑去噪算法要好的多，它可以方便明确地区分开高频噪声。但是，在信号已经足够平滑的情况下，小波变换去噪不会有明显效果。小波变换的原理就在于将原信号根据不同的频率小波分解为不同的小块信号，再分析各小块信号的特性，最后按照一定的阈值对其区分并重构，从而实现平滑去噪。小波变换信号去噪的一般步骤分为如下3个步骤：（1）确定要变换的小波类型，同时选择需要分解的层次，最后分解计算。（2）确定用去区分系数的阈值定义。（3）进行具体的小波重构。从上面的3个步骤中可以看出，小波变换的重点就在于如何选取合适的区分阈值，使得信号被很好地区分开来，因此信号去噪后表现优劣直接取决于阈值的选取。通常，利用MATLAB解决小波变换问题中，阈值的选取方法有三种：①系统自带阈值；②人为给定阈值；③强制给定。上述的三种阈值选取方式中需要使用到的MATLAB系统自带函数有：ddencmp——分解信号的默认阈值wdencmp——去噪处理。wthresh——阈值量化处理。针对本文研究所使用的大米近红外光谱数据，分别采用以上三种方式进行小波变换，变换方式是：db5小波+3层分解，处理后的数据光谱图如下。图3-8小波变换后光谱图根据图像可以得出如下结论：针对本次测量使用的大米近红外光谱数据而言，方法①<方法②<方法③。另外，除了阈值选取会影响结果外，分解所用到的小波函数的选择也会使得小波变换去噪处理后的结果变得不一样。因此，使用小波变换时要根据具体问题选择不同的小波函数，现阶段大部分使用的方法是利用大量实验验证或是依靠经验来选择最佳的小波函数。由于小波变换涉及到的优化参数较多，处理过程相对复杂难懂，于是，针对本文所要研究的大米光谱数据而言，采用下面将要介绍的EMD算法对光谱数据进行去噪处理，具有更好的性能和效果。基于EMD算法的光谱去噪方法EMD算法介绍经验模态分解（EmpiricalModeDecomposition,EMD）经常被用于解决非线性处理问题，他的发明者Huang教授认为任何信号都是由若干个IMF组成的，IMF即本征模态函数（IntrinsicModeFunction,IMF）。而EMD算法的主要任务就是找到IMF分量，并对它们进行相应的分析区分。EMD算法从根本上解决了小波变换中小波基和分解尺度选择的难题。在本文的研究中，考虑到大米光谱数据的信号特征，使用本文规定好的阈公式对信号进行分析分解。EMD算法具体步骤如下，其中是由两个信号tone（有用信号）和chirp（噪声信号）叠加而成，如下图所示，对于信号而言。图3-9信号组成光谱示意图利用三次样条插值法顺序连接局部最大值数据点，可以获得上包络线；同理，连接局部最小值数据点，可以获得下包络线。将上下包络线的均值定义为。提取瞬时局部振荡值。重复步骤(1)-(3)直到接近于零，则即为一个IMF，定义为。计算余项。将余项赋值于，重复步骤(1)-(5)直到下一个余项为单调函数。因此，原始信号可以重构为：(3-1)式中为第个IMF，为余项。通过EMD算法对实际信号进行分析时，迭代次数可能就无限进行，这就需要对IMF的迭代次数进行限定，公式如下：(3-2)其中，为信号长度，当小于某值时（通常为0.2~0.3间），筛选动作停止进行。图3-10信号原始数据及各IMF分量数据图图3-11前n个IMF重构信号图图3-12后n个IMF重构信号图从图中分解后的IMF分量可以看出，IMF2分量和有用信号在形状上最为相近，而从代表倒序（从后向前重组）的后n个IMF分量重构后的光谱图中可以看到，c2f6重构信号与原本数据中的有用信号基本一致，即IMF2与IMF3、IMF4、IMF5、IMF6以及余项RES组合重构后的c2f6即为去噪结果。光谱数据的去噪处理本论文中大米数据如下图所示，绿色为五常、洋红色为响水、蓝色为越光、红色为银水，从图中可以看出大米近红外光谱数据普遍存在噪声，且噪声集中分布，个别波段区间几乎不含噪声，光谱数据混叠严重，受噪声影响较大，分析产生该现象的原因如下：由于本课题研究中使用的近红外光谱仪器是偏中红外波段的仪器，所以对应的光谱数据图左端也即波数较小处（靠近8000cm-1处）的噪声较大，主要是由仪器本身造成，后期特征波段选取时的意义不大（不过在下一章也证实了该区间波段用于产地分类的用处不大）；而光谱图中间区段（即5200cm-1左右处），通过前述光谱分析可知，该区间有一较明显的吸收峰，可用于表征水分含量，所以该处噪声产生的主要原因是水分含量不稳定造成的。本文中使用的是MATLABR2016a版本，该版本软件中本身不包含EMD算法，首先需要下载EMD算法工具箱以及时频工具箱，将下载好的package_emd和tftb-0.2文件添加到MATLAB的toolbox文件夹中，并将其添加到工作路径中，完成编译。图3-13大米数据原始光谱图根据上一节文章内容阐述可知，被处理信号为，其中是对于后续数据处理有用的部分，而为噪声干扰的部分。EMD算法将被测信号分解成n个本征模态分量和一个余项式，分解后得到的结果是将所有的IMF分量按照频率高低进行顺序排列的结果。通常认为噪声是一种高频信号，因此构成噪声信号的IMF分量主要集中在高频部分；相反地，对后续处理数据有用的信号主要集中在低频部分。因此，只要找到一个临界阈值序号，将IMF分量分为低频和高频两部分，有用信号和噪声信号分别由这两部分重构组成，如下式所示：(3-3)所以想要获得去噪后的有用信号，主要的目的就是寻找分解过后的各相关模态的分界阈值序号，以此来可以区分噪声信号模态和纯净信号模态，其中EMD算法中的取值可以通过MATLAB语句获得，语句如下：imf=emd(y);[m,n]=size(imf);k=ceil(m/2)+1;通常曲线波形的突变主要集中在分解序数的一半以右，分界序号可以用来区分噪声模态和信号模态，以单独的一个五常大米近红外光谱数据来看，通过EMD算法进行分解并重构。其中，f2cn表示将分解后的前n个IMF分量进行组合重构，c2fn表示将将分解后的后n个IMF分量进行组合重构，EMD分解后的各分量以及按照前向和反向重构的可视化程序语句为：emd_visu(y,1:length(y),imf)。具体的处理结果如下图所示：图3-14单个大米数据的EMD处理可视化由于有用数据集中在低频区域，所以可以通过c2fn的组合看出EMD算法分解重构后何时才能能够达到最佳的状态，为了方便观察和处理重构后的数据，将上述单个大米数据EMD分解后的后n个IMF分量重构的c2fn单独处理出来，得如下图像矩阵，其中红色代表重构后的光谱数据，绿色表示原始的大米光谱数据。图3-15单个大米的c2f8（左）、c2f7（右）重构图图3-16单个大米的c2f6（左）、c2f5（右）重构图图3-17单个大米的c2f4（左）、c2f3（右）重构图图3-18单个大米的c2f2（左）、c2f1（右）重构图通过上图可以观察得出结论：该大米样本数据的最佳去噪组合为c2f5，即IMF4分量及以后分量与余项的组合重构数据；的选择越小，处理后的数据越贴合原始数据，但去噪效果越差；的选择越大，处理后的数据越平滑，但在原本噪声含量不大的波段处偏移程度越大；虽然c2f5为该大米样本数据全波段最佳的重构去噪组合，但仔细观察可以发现，在靠近横坐标左端即波数数值较大处（此处噪声较大），c2f5在该区域的去噪效果不如c2f4好，但c2f4在中间区域（噪声较小，比左侧区域大，比左侧区域小）以及右侧区域（几乎不含噪声）存在去噪严重的现象，其处理后的数据与原始数据偏移较大、平滑过渡，数据无法使用。针对上述出现的问题，对大米光谱数据进行分段处理，如下图所示。将光谱数据分为三个区间，大噪声区间：7800~5970cm-1（对应数据点：1~950）、小噪声区间：5600~5110cm-1（对应数据点：1142~1396）、无噪声区间（除上述两区间以外区间）。将前两个噪声区间分开进行EMD算法分解重构，获得去噪后的数据，将去噪后的数据与无噪声区间数据结合，得到最终的全波段大米近红外光谱去噪数据。图3-19大米数据的分区段光谱图基于EMD算法的大米近红外光谱数据去噪具体步骤如下：(1)导入全波段光谱数据信号。(2)截取选定的两个噪声较大的特殊波，即波数为7800~5970cm-1（对应数据点1~950）、5600~5110cm-1（对应数据点1142~1396）这两部分。(3)利用EMD算法分别对两波段进行分解，获得一系列IMF。(4)每个特殊波段的相关模态序号值由下式确定：(3-4)对原始信号进行重构，得如下：(3-5)对每个大米光谱数据进行上述处理，并将处理出来的结果保存到.mat文件中。EMD算法的去噪性能评价本节通过对比前述的不同光谱去噪方法，对基于EMD算法分解重构后的光谱信号进行去噪程度方面的评价。对比的方法包括有：小波变换法、五点平滑、Savitsky-Golay卷积平滑法、EMD去噪这四种方法。仿真实验基于R2016a版本的MATLAB计算平台，计算机参数为主频：2.5GHz，IntelCorei5-7200U处理器，内存为4GB。其中，小波变换：db5小波、3层分解、给定软阈值去噪；五点平滑次数：100；SG平滑：多项式次数5、窗口宽度20，步骤如下：（1）由前文分析可知，选取噪声较大的波段7800~5970cm-1为实验波段；（2）对于实际采集的光谱信号来说，由于不含噪声的信号未知，因此不能通过信噪比，但根据信噪比的定义，我们可以仿照信噪比的原理，提出与之类似的概念，即“信差比”，公式与信噪比公式一致即下式（3-6）所示：(3-6)其中为原始信号，为重构信号。信噪比的定义中，首先基于的条件是：不含噪声的纯净信号是已知的。公式中，分子中的原始信号指的是纯净信号（即信号中不包含任何噪声的信号），分母为原始信号（不含噪声）和重构信号之差的平方值，所以当重构信号包含的噪声越小时，分母中原始信号（不含噪声）与重构信号差值的平方值就越小，所以，在分子原始信号（不含噪声）不变的情况下，包含噪声越小，分母越小，信噪比就越大。然而，不同于信噪比数值越大噪声越小的道理，由于不含噪声的纯净信号未知，所以本设计中的“信差比”定义是：分子为原始信号的平方值，但这里的原始信号指的是大米近红外光谱的原始数据，即本身就包含噪声的光谱数据，当去噪效果好时，分母中的原始信号（含噪声）与重构信号差值的平方越大，在分子原始信号（含噪声）不变的情况下，去噪的效果越好，分母就越大，“信差比”就越小。所以采用这种方法的结果中，“信差比”越小的方法，去噪效果越好。（3）将不同分解方法的重构结果与原始包含噪声的数据进行对比，采用输出“信差比”的评价准则从数值角度验证不同分解方法信号重构结果的优劣，得出具体数值。（4）作出四种去噪方法所获得的信号图像，利用数形结合的方法分析各方法性能。根据以上步骤进行实验验证，得到的数据和图像如下表3-1及下图3-20：表3-1各去噪方法性能评价表去噪方法“信差比”/dB小波变换28.4711五点平滑29.5761SG卷积平滑法29.5425EMD算法去噪27.9900图3-20各去噪方法重构信号效果图表3-1和图3-20为针对不同的去噪算法重构信号结果的输出信差比，从图表中结果可以得出如下结论：针对大米近红外光谱中高噪声波段