茶籽含油量高光谱检测技术研究

茶籽含油量高光谱检测技术研究摘 要：为了研究油茶籽含油量高光谱检测技术，本试验采用了神经网络、主成分分析法、偏最小二乘法等算法，分析处理了油茶籽的高光谱数据，提取相应特征值（面积参数） 。对高光谱数据进行建模，通过模型预测油茶籽含油量的面积参数。通过预测面积与实际面积的偏差值，评价所建模型的预测能力。研究包括使用 FieldSpec HH 2 光谱仪测出 27 组油茶籽光谱，然后使用气相色谱法对样品进行标定，从而得到高光谱模型和实际含油量的对应关系，并进行外部验证。其具备对预测集油茶籽含油量的检测。关键词：近红外；ANN 神经网络；主成分分析；偏最小二乘法2Research on Detection of Camellia Oilfera Abel Seed Oleaginousness by Hyperspectral TechnologyAbstract: Artificial neural network, principal component analysis, partial least squares etc. were used to analyze the hyper spectral data of oil tea camellia for the research on detection of Camellia Oilfera Abel seeds oleaginousness by hyper spectral technology, extracting the appropriate characteristics（area parameters）. As a result, a model was built to predict the area parameters of oil tea camellia seed oleaginousness successfully and to study the deviation of the forecast area and actual area, to evaluate the predictive ability of the model. The result proved the model was capable to achieve the detection of camellia seed oil content This approach was composed of four major procedures: Measuring spectrum corresponding to the 27 group of seed by the Field Spec HH2 hyper spectral instrument, Quantitative analysis of the sample by gas chromatography technology, Building a model of the correspondence between oleaginousness and spectral data and recalibrating externally. Especially, a novel methodology, model based on spectral analysis method, was proposed to accomplish measuring the obscured oleaginousness of Camellia Oilfera Abel seed by utilizing the property of model.Key words: NIR (Near Infrared Spectroscopy); ANN (Artificial Neural Network); PCA (Principal Component Analysis); PLS (Partial Least Squares)3目 录摘要1关键词11 前言21.1 课题背景21.1.1 油茶籽的简介21.1.2 含油量检测技术的目的和意义21.2 研究进展和国内外现状21.3 研究内容与技术路线31.3.1 研究内容31.3.2 技术路线42 总体方案设计42.1 近红外光谱定量分析理论依据 42.1.1 近红外光谱产生机理42.1.2 油茶籽分子光谱特性42.2 试验原理52.3 试验方案63 试验方法与过程63.1 试验材料63.2 样品的光谱采集73.3 采集的部分高光谱数据93.4 高光谱数据处理的必要性114 油茶籽光谱特性124.1 分析软件和方法124.2 光谱数据预处理1344.2.1 平滑处理提高信噪比134.2.2 峰值及峰位134.3 光谱数据降维174.3.1 PCA 降维174.3.2 NEWC 降维224.3.3 反射峰积分264.4 BP 神经网络建模274.5 外部验证315 结束语32参考文献33致谢34附录3511 前言 1.1 课题背景1.1.1 油茶籽的简介茶油全称是油茶籽油，是油茶树(Camellia oleiferaAbel )所产富含脂肪的种子经压榨或浸提得到的植物油脂。茶油的高营养价值与其脂肪酸组成有关，茶油主要由脂肪酸组成，即油酸(C18:1)74%87%、亚油酸(C18:2)7.0%14%、饱和酸 7%11% 1。故采用油茶籽油的特征指标即脂肪酸组成，作为油茶籽油真实属性的主要判定依据。茶油中油酸的含量为天然植物油中之冠，含有的亚油酸、亚麻酸等多不饱和脂肪酸为人体生长和生理活动不可缺少的必需脂肪酸，茶油食用后人体的消化吸收率达 97%。优质茶油富含生理活性物质甾醇、生育酚、角鲨烯等，对提高人体抗病能力，延缓衰老等有重要作用 2，其市场价格大约是普通食用油的 36 倍 3。1.1.2 含油量检测技术的目的和意义随着科技的进步和人们健康意识的增强，相比食品的色香味，人们更看重食品的营养价值和安全性。消费者普遍关心的是市售品牌油、散装油、以及餐饮用油的品质。近年来出口和国内市场的食用油脂中掺伪现象也日趋严重。茶油因其种类不同、营养价值不同而价格差异很大。一些生产经营者为了获取暴利，在高价茶油籽中掺入廉价的茶油籽。甚至还有的厂家将过期变质油茶籽掺入合格油茶籽中以次充好，但茶农和消费者很难从外观上识别掺杂的油茶籽。掺伪茶油籽不仅影响品质和营养成分，而且危害到茶农的利益和消费者的健康。为了保护合法生产经营者和消费者的利益，能否找到快速准确检测油茶籽品质和进行油茶籽种类鉴别的方法成为当务之急。传统的检测油茶籽含油量主要有理化检验方法、气相色谱法、高效液相色谱法及同位素比值法等方法。传统化学方法往往需要多种化学仪器和试剂，样品需要预处理，操作繁琐，耗时费力。因此，寻求出一种能快速、准确、在线、无损检测油茶籽品质的新方法对提高油茶籽市场实时质量监控能力，规范油茶籽市场，保障茶农权利等方面具有十分重要的意义。1.2 研究进展和国内外现状国内外许多学者利用电学特性、声学特性、光学特性、声波振动特性、太赫兹、核磁共振、X 射线、电子鼻、机器视觉等技术对种子的品质进行无损检测 4。种子品质无损检测是在不破坏被检测种子的情况下，应用一定的检测技术和分析方法对其外部特征和内部品质加以测定，并按一定的标准对其做出评价的过程。油茶籽含油量用2传统的检测方法进行测量费时、费力，且内部组织成分的测定需依靠破坏性检测方法，受人工影响较大，这严重影响了对油茶籽进行快速、无损、准确的评价。因此研究快速、无损、准确的油茶籽含油量的检测技术，对推动农业的健康、持续发展，提高油茶籽的交易价格，增加果农收入，具有十分重要的现实意义。众所周知，光谱分析是自然科学中一种重要的研究手段，光谱技术能检测到被测物体的物理结构、化学成分等指标。高光谱遥感技术在植被生化成分定量提取与作物农学参数估算研究上的逐步深入，高光谱数据与叶绿素含量( Chlorophyll Content, CH.C) 、叶绿素密度(Chlorophyll Density, CH.D)等农学参数之间的相关关系的研究已日渐增多 5。吴长山等通过微分光谱与水稻、玉米 CH.D 的相关分析，选出 CH.D 相关性高且受植被种类影响很小的特征波段，建立了统一的线性模型，模型精度达 80.6%6。国内有关含油率 NIR 模型分析的研究有油菜、玉米、黄豆、棉籽等植物油种类，如李钧等研究的油菜籽含油量NIR 模型的内部交叉验证的决定系数和均方差分别为 0.98 为 0.7557。方彦等建立的玉米 NIR 模型其定标集、 检验集的相关系数分别为 0.958 和 0.9578。这些模型都有很好的预测结果，为近红外光谱分析技术在植物油含量的快速检测奠定良好的理论基础。原姣姣等人 1通以索氏提取法分析了 30 份油茶籽的含油量,通过多种预处理方法和回归方法建立了较精准的分析模型。其校正相关系数(RC)0.92567、外部验证系数(r)0.8978 其最优波段为 4200.204088.35 cm -1 和 4666.894639.89 cm -1。这模型为高光谱检测技术的发展起到了积极的意义。近红外光谱定量和定性分析在农副产品及食品的研究方面有较广泛的研究，而关于高光谱在油茶籽品质检测与掺假鉴别分析尚无较深入的研究。1.3 研究内容与技术路线1.3.1 研究内容该课题研究目标是：采用不同的高光谱预处理，降维方法建模，得出最佳模型。具体的研究内容包括：（1）收集油茶籽样品，其中包括了含油量差别较大的品种。使茶油籽含油量的变幅尽可能包含待检测茶油籽的含油量，便于选出有代表性的样品参加标定计算。（2）采集油茶籽的高光谱。（3）对高光谱进行预处理，提高信噪比。（4）降维处理，便于软件分析。（5）特征值提取，对高光谱进行主成分分析，NEWC 和积分面积处理，提取其特3征值。 （6）使用 BP 网络对高光谱进行建模。（7）外部验证，验证模型的可行性。1.3.2 技术路线（1）选择有代表性的校正集样本，并测量其近红外光谱。（2）采用气相色谱法进行化学成分标定。（3）根据测量的光谱和化学数据通过合理的化学计量学方法建立校正模型。（4）根据校正模型对预测集样品的组成及性质进行测定。2 总体方案设计 2.1 近红外光谱定量分析理论依据 2.1.1 近红外光谱产生机理在分子中，如果成键的两个原子的正负电荷中心不重合，这两个原子就形成一对电偶极子 9。电偶极子以一定的频率振动时，能吸收光子，因而会对吸收谱产生影响。若某物质中的原子完全非极性且完全对称，该物质在不同波长处的吸收就相等，吸收谱就会是一条直线 10。实际中的分子并非如此，其吸收谱中的每一波峰或每一波谷都反映了分子的组成结构特征。这就是光谱分析的物理基础。近红外波段的吸收一般是由质量最小的氢原子 H 引起。氢原子不仅能引起分子的基频吸收，还会产生非简谐振动并引起较强的倍频吸收 11。H 基频吸收在28003700cm -1，绝大多数倍频在 450025000cm -1，(约 4002200nm)内，40005500cm -1(约 1800 2500nm)为合频区域。H 的基频和倍频与分子中的其它振动会产生了遍布整个近红外区域的合频吸收 12。每个分子都会有许多个吸收带，吸收带的强弱还受分子浓度的影响，这是利用近红外光谱分析样品的浓度的化学基础 13。2.1.2 油茶籽分子光谱特性查找相关文献和化学分析方法，确定油茶籽中的主要成份为油酸 1。根据红外光谱的波数范围为 1280010 cm -1（0.751000um） 。近红外区处于可见区和中红外之间，吸收带主要是由低能电子跃进、含氢原子团（如 O-H、N-H、C-H ）伸缩振动的倍频及组合频吸收。基于 O-H 伸缩振动的第一泛音吸收带出现在 7100 cm-1，可以测定各种试剂中的水分，以及甘油、肼、发烟硝酸等 14，可以定量测定酚、醇、有机酸等。基于羟基伸缩的第一泛音吸收带出现在 33303600 cm -1，可以测定酯、酮。它的测量准确度与紫4外、可见吸收光谱相当 15。另外，漫反射可测定未处理的固体和液体试样，例如油茶籽 16。在通常情况下，分子都处于基态振动，一般极性分子吸收红外光主要属于从基态到第一激发态之间的跃迁，其振动能量变化为（1）/2khE对应的谱带称为基频反射带或基本振动谱带，若用波数表示，可以表示成(cm-1) （2）/1370其中由公式（2）可见键力常数 k 应采用 N.cm-1 为单位，而 应采用原子质量单位(u) 为单位。由公式（2）可见理论上油茶籽的油酸 C18H34O2 的基频为：/k1370=1370*sqrt(5.1/(18*12+34+32)*(18+34+2)=1353.873167 cm-1根据这个公式可以测量各种类型的化学键力常数 k，一般来说单键常数的平均值约为 5N.cm-1，而双键和三键的键力常数分别大约是此值的两倍和三倍 17。相反可以利用实验得到的键力常数的平均值估算吸收频率应为 2993cm-1 而实际光谱实际测值为2885cm-1。由公式（ 2）可见化学键的力常数 k 越大，原子折合的质量 就越小，则原子间的振动频率就越高，吸收峰将出现在高波数区；反之，将出现在低波数区 18。例如，叁键双键单键，吸收峰出现的位置不同；叁键（2222cm -1）双键（1667cm -1）单键（1429cm -1）此外多原子分子的振动，不仅包括双原子分子沿核-核的伸缩运动，还有键角渗入的各种可能的变形振动：伸缩运动和变形运动 19。综上所述，油茶籽分子这种多原子分子，由于组成原子数目增多，加之分子原子排布情况的不同，及组成分子的键或基团和空间结构的不同，其振动光谱远比双原子复杂得多 20， 。故其光谱数据具有很多的变量。因此在分析油茶籽高光谱数据时，需对其进行优化，分离出相应分子的高光谱值，才能建立相应的模型。2.2 试验原理本试验主要是检测油茶籽油脂肪酸含量。采用 NIR 数据结合气相色谱法测定的含量，建立高光谱与含油量的模型，从而达到快速用 NIR 测量油茶籽的含油量，甄别油茶籽的优劣属性 20。高光谱的数据与油茶籽分子振动有密切的关系。分子的振动自由度可以通过红外5光谱的吸收峰来体现。从原则上讲，每一个振动自由度相当于红外区的一个吸收峰，但实际的红外吸收峰的数目常少于振动自由度的数目。因为不伴随偶极变化的振动没有红外吸收峰；振动频率相同的不同振动形式会重叠。在建立油茶籽模型之前，高光谱数据需要进行预处理，然后再利用优化之后光谱，分析出其相应的特性，建立高光谱数据与含油量的模型。检测未知油茶籽含油量，测量其光谱数据，并优化处理，进行模式识别即可得到含油量。2.3 试验方案选取 27 种相同质量，不同品种的油茶籽并编号；在相同的光照、湿度和温度下，利用高光谱分析仪测量其对应的高光谱数据并记录；使用 ViewSpec 软件中的修正功能对数据进行修正，并生成.mat 文件；用 SPSS、SAS、MATLAB 等软件进行PCA、 NEWC、积分面积等分析，得到油茶籽的特性；使用气相色谱法对其进行含油量定标；结合实际油茶籽的含油量和高光谱数据特征值，用 BP 神经网络等进行建模分析，建立起高光谱与含油量的相关数据处理模型，从而实现利用高光谱对样品进行含油量的检测。图 1 总体方案框图Fig. 1 Diagram of the overall program3 试验方法与过程3.1 试验材料（1）样品筛选。27 组油茶籽样品(产自江西南昌、湖北武汉、江西平南、徐州、淮安等，部分样品产于长沙、常德市、益阳市、等地网购)。由于不同地区的生长条件诸如土壤、光照、温湿度、肥力、季节等因数的不同，所搜集的样品包括了含油量差别较大的品种。这样才能使茶油籽含油量的变幅尽可能包含待检测茶油籽的含油量，便于选出有代表性的样品参加标定计算。表 1 油茶籽的种类分布表Table 1 Camellia species distribution table编号 油茶籽产地 重量（500g） 编号 油茶籽产地 重量（500g）用已建的模型计算预测集样本的含油量降维后的光谱和含油量建模气相色谱法测量含油量光谱数据物理特性27 组油茶籽样品光谱使用SPSS、 MATLAB 等软件进行降维处理模型的评价61 江西南昌 1 15 湖北观音湖 22 湖北武汉 1 16 湖北仙居定顶 23 广西平南县 1 17 湖北木兰湖 24 徐州 1 18 湖北省阳新县 25 淮安 1 19 湖北大悟山 2编号 油茶籽产地 重量（500g） 编号 油茶籽产地 重量（500g）6 山东日照 1 20 湖北小悟山 27 江苏宿迁 1 21 上海 18 广东沐阳 1 22 长沙县 149 广东梅州 1 23 长沙葵花山庄 410 江苏 1 24 长沙黄花机场 1411 湖北随州 1 25 湖南绥宁 412 浏阳小河中学 1 26 广西柳州 113 江苏花卉种植园 1 27 福建 114 湖北悟峰 2（2）样品制备。由于环境温度，特别是湿度对光谱测量有较大的影响，所以必须把环境温湿度控制在一定的范围内，以保证采集光谱过程中免受温湿度变化的影响。油茶籽中的水分对其他成分的标定的影响不可忽视，因为水的基频吸收带为 3650cm-1、3750cm -1（ -OH 的对称与不对称伸缩振动）与 1595cm-1，而-OH 的合频吸收带为5160cm-1，故应使水分严格保持一致，才能获得稳定的光谱。其具体的实施办法是将待测的各组油茶籽在测量之前同时进行干燥。3.2 样品的光谱采集试验采用 FieldSpec HandHeld 2 光谱仪（图 2）进行实验。图 2 FieldSpec HandHeld 2 光谱仪外表结构Fig. 2 Feature of the FieldSpec HandHeld表 2 光谱仪的外部个接口Table 2 Import and outport of the hyper spectrum instrument序号 名称 序号 名称1 25前视场角入光口 8 远处控制接口2 激光瞄准器 9 电源接口3 三角架连接口 10 操作屏（可转动）续表 174 扳机连接口 11 电池槽，需 4 节 AA 电池5 USB 迷你-B 接口，用于电脑控制光谱仪操作 12 选配瞄准镜连接处6 USB 接口 13 电源接口7 USB 接口 14 操作屏（可转动）把光谱仪放到被测样品上方（注意：光谱仪的方位应当与采集白板参比光谱时相同） ，此时界面上显示的就是相对反射光谱；按拇指键存储当前的光谱曲线。保存完成后能够听到提示音；注意：保存数据时界面是否提示出现饱和。如果出现饱和则存储数据前必须重新优化。将光谱仪对准白板，并使得白板已经被最佳照明并且在镜头视场内只有白板。点击 OPT 图标优化光谱仪的积分时间。注意每间隔 15-20 分钟或者照明条件以及环境条件（比如云层覆盖、湿度变化、太阳移动等）改变时，重新优化对数据质量是有益处的。事实上，每次更换油茶籽采集光谱数据之前都重新优化；点击 DC/WR 图标，仪器会自动重新采集暗电流，几秒钟之后界面上显示一条反射率数值为 1.00 的平直线，同时界面左上角显示 R 字符如图 3； 图 3 高光谱仪器界面Fig. 3 Workplace of the fieldSpec handheld 28图 4 油茶籽含油量检测的现场图Fig4. Scene of measuring the Camellia seeds oil content3.3 采集的部分高光谱数据表 3 五组高光谱 900930nm 数据Table 3 Five groups of hyper spectral data through 900 to 930nm波长 nm第 1 组反射率（%） 第 2 组反射率（%） 第 3 组反射率（%） 第 4 组反射率（%） 第 5 组反射率（%）900 0.467 0.539 0.532 0.579 0.527901 0.47 0.534 0.532 0.581 0.527902 0.474 0.53 0.532 0.588 0.525903 0.475 0.534 0.534 0.591 0.527904 0.473 0.537 0.529 0.589 0.529905 0.471 0.538 0.522 0.586 0.53906 0.476 0.537 0.53 0.589 0.527907 0.477 0.537 0.536 0.591 0.525908 0.472 0.539 0.535 0.591 0.525909 0.475 0.539 0.536 0.588 0.529910 0.476 0.538 0.536 0.587 0.531911 0.471 0.539 0.535 0.591 0.527912 0.475 0.54 0.532 0.589 0.526913 0.482 0.54 0.529 0.583 0.526914 0.48 0.539 0.535 0.582 0.528915 0.479 0.539 0.536 0.584 0.527916 0.48 0.542 0.531 0.587 0.525917 0.478 0.542 0.531 0.589 0.523918 0.478 0.541 0.532 0.588 0.525919 0.48 0.538 0.531 0.583 0.531920 0.48 0.541 0.535 0.584 0.531921 0.48 0.543 0.54 0.587 0.528922 0.487 0.535 0.54 0.587 0.5319923 0.487 0.534 0.538 0.586 0.533924 0.482 0.541 0.536 0.586 0.533925 0.482 0.539 0.53 0.587 0.526926 0.486 0.537 0.529 0.588 0.521波长 nm第 1 组反射率（%） 第 2 组反射率（%） 第 3 组反射率（%） 第 4 组反射率（%） 第 5 组反射率（%）927 0.489 0.538 0.534 0.588 0.526928 0.489 0.541 0.534 0.59 0.529929 0.489 0.54 0.531 0.591 0.527930 0.474 0.53 0.532 0.588 0.525注：其中光谱的反射率，为相对反射率，即测量的光谱值相对于白板校验光谱值的百分比。输入 MATLAB 后对其中的 4 组和 7 组直接描点。具体的程序如下： load(BP.mat) plot(P(:,4),:,LineWidth,2) hold on plot(P(:,7)%P 为优化后的 75127 矩阵的光谱数据程序运行的结果如图 5 的左图所示。 plot(P(:,5),:,LineWidth,2) hold on plot(P(:,6),LineWidth,3) hold on plot(P(:,8)程序运行如下图 5 的右图所示（其中横轴代表波长，纵轴代表相对反射率） 。续表 310图 5 两组不同的高光谱数据Fig. 5 Two group of spectral data由图 5 的左图可见，第 4 组为实线表示的图形，而第 7 组数据为虚线表示的图形。可以看到高光谱子在 8001000nm 范围内，不同的波长段有不同的反射峰，具有很大的敏感度，能够利用相应的光谱特性来建立相应的模型。再由图 5 的右图可见，第 5 组为 2 号实线表示的图形，而第 6 组数据为 3 号实线表示的图形。第 8 组为最细的实线显示的图形，也可以看到高光谱子在 800900nm范围内，有一个很明显的波峰，能够利用相应的此光谱特性来建立相应的模型。3.4 高光谱数据处理的必要性1）在收集样品近红外光谱时，一些主要来自高频随机噪音、基线漂移、信号本底、样品颗粒大小和光散射等的噪音信息会随之带入光谱中，这将会严重干扰近红外反射强度与样品中各项指标的关系，并直接影响到近红外回归分析方程的可靠性和准确性。因此，在光谱分析时，应综合光谱预处理和数学处理来滤除这些噪音。本课题的油茶籽高光谱分析软件中，包括光谱预处理和数学处理，其中光谱预处理包括趋势变换法、标准正态变量转换（Transformation of Standard Normal Variable，简称 SNV）法等。趋势变换法是用一个人工模拟的模型（一般为多项式）来描述基线，然后将整个波长区每一点的光密度减去该点的基线值，用以消除基线的漂移。标准正态变量转换（SNV）即每个波长点处的光密度减去整个光谱区光密度的平均值，再除以整个波长区光密度的标准差，该转换方法可以减少非线性弯曲作用产生的偏差。数学处理主要有导数处理和平滑处理两种方法。其中导数处理可以提高光谱的分辨率以及减小基线的漂移，而平滑处理可以去掉高频噪音对光谱信号的干扰。对第 4 组的油茶籽数据的导数处理如图 6。图 6 一组高光谱数据的微分11Fig. 6 Derivative of the reflection由图 6 可见，光谱在短波和长波段，具有很够的离散性，具有很多光谱的特性。并且其导数也可以直接作为光谱的处理的原始数据，其抑制了对基线的漂移。2）光谱数据较大，例如每组油茶籽，采集 30 次光谱数据，每次采集 751 个数据，就是每组油茶籽对应 30751=22430 个值，每个值为四位浮点数。一般软件无法完成油茶籽高光谱数据计算，需要对其进行降维处理。本研究采用的算法有主成分分析方法、自竞争神经网络、积分面积等。采用降维处理后的数据，能够表征原谱图的特性，能够直接快速地使用这种特性进行分析研究。4 油茶籽光谱特性4.1 分析软件和方法近红外技术和其他分析技术一样有其长处，也有一定的局限性。例如因为近红外的反射率较低，不适宜做含量过低的样品、微量样品的显微分析。因为近红外光谱谱峰复杂重叠不宜进行官能团定性分析。因此，应注意近红外光谱分析技术的综合应用。例如为了实现对复杂物的综合品质性状的分析，需要运用多种仪器分析手段进行分析。对此可以运用“ 数据融合”（date fusion）技术，对多种仪器分析得到的信息在一定准则下加以综合、分析、过滤、相关与合成。本课题研究了 MLR（multiple linear regression, 多元线性回归） 、逐步回归(stepwise regression, SMR） 、主成分分析（PCA） 、主成分回归（principal component regression, PCR）与偏最小二乘法（PLS ） 、人工神经网络（ANN）和拓扑（Topological）等。MLR 和 SMR 发在分析样品时只采用了一些特征波长点的光谱信息，其他点的信息被丢失，易产生模型的过适应性（over fitting） 。PCR 和 PLS 的显著特点就是利用了全部的光谱信息，可以压缩所需要样品数量，将高度相关的波长点归于一个独立变量中，根据位数不多的独立变量建立回归方程，通过内部验证（cross validation）来防止过模型现象，比 MLR 和 SMR 分析精度高。油茶籽数据结果采用 SPSS 统计软件，其分析包括描述性统计、均值比较、一般线性模型、相关分析、回归分析、对数线性模型、聚类分析、数据简化、生存分析、时间序列分析、多重响应等几大类，每类中又分好几个统计过程，比如回归分析中又分线性回归分析、曲线估计、Logistic 回归、Probit 回归、加权估计、两阶段最小二乘法、非线性回归等多个统计过程。SPSS 也有专门的绘图系统，可以根据数据绘制各种图形。很容易对光谱数据进行初步的处理，研究相关算法。因为油茶籽光谱的数据量很大，故使用 MATLAB 分析软件，使用 Simulink 12Toolbox、Signal Processing Toolbox、神经网络工具箱等工具箱对光谱的数据进行分析处理。本研究是采用全波段进行回归分析构建回归方程，因此，需用一些标准来衡量回归方程的优劣。这些参数主要包括定标标准差、预测标准差、定标相关系数、预测相关系数、定标决定系数、预测决定系数。构建回归方程时，一般要求相关系数与决定系数很高同时具有较小的各类标准差。4.2 光谱数据预处理4.2.1 平滑处理提高信噪比平滑技术是降低谱图噪声的一种方法，通过对数据做一定的数据平均计算，使误差重新分配，从而达到降低噪声的目的，尤其是降低高频噪声。 可用于平滑技术的数学方法很多，本课题采用的是（Savitzky-Golay）法。它是根据最小二乘法而采用的多项式近似法，可选择不同的平滑点数。平滑点数又称平滑级数，必须是奇数值，如 5,7,9,11以 5 点平滑为例，它取第 15 个数据点的 5 个 Y值的平均值作为第 3 点的 Y 轴值；取第 26 点的 Y 值的平均值为第 4 点的 Y 值；取第 37 点的 Y 值的平均值为第 5 点的 Y 值；依次类推。这种平滑方法是将各数据点完全等同地进行处理，并不区分在峰尖与峰谷处的数据点的差别。这种数据处理方法可能会使谱峰变形，使峰值发生位移。平滑处理时所取级数愈高，信噪比得到改善的程度愈大，但是谱峰变形也愈严重，降低了光谱分辨率。、图 7 一组油茶籽光谱数据不同的平滑处理Fig. 7 Different scales of smoothing图 7 的左图和右图是第四组油茶籽样品，分别采用的平滑级数为 5 和平滑级数为15 的时候对光谱的处理结果。从图中可以看出，左图光谱较右图显示更为精细的结构，平滑级数愈高，反射峰的失真愈多。但是平滑级数为 15 时同时反映了基本相同的特13征信息，而且光谱较平滑，数据量较少。所以在相同有用信息的情况下，采用级数为15，即右图的光谱。4.2.2 峰值及峰位（1）matlab 的拟合系数数目的选择。理论上一个振动的自由度，在红外光谱上相应产生一个基频反射带。实际上，绝大多数的化合物在红外光谱图上出现的波峰数，远小于理论上的波峰数，原因如下：没有偶极距变化的振动，不产生红外吸收，及非红外活性；相同频率的振动会重叠；仪器不能区别那些频率十分相近的振动，或因反射带很弱，仪器检测不出。峰值及反射峰峰尖的横坐标值，峰尖与峰谷的位置可通过多种数学方法计算出来。峰尖与峰谷的斜率变化均为 0，它们的位置可由一阶导数值确定，计算机由左到右对谱图逐步求导就可得到反射峰的峰尖和峰谷的位置。在红外光谱中常采用的坐标单位有横坐标以波数为单位的吸收峰、横坐标以波数为单位的投射谱、横坐标以波长为单位的吸收谱和横坐标以波长为单位的投射光谱四种情况。横坐标以波数为单位的吸收谱的横坐标值从左到右的波数值逐渐减小（例如中红外谱为 4000400cm -1） ，当计算机沿横坐标由左至右依次计算 dA/d时，d 总是负值，在尖峰附近，dA 由正值变化为负值，一阶导数值相应的由负值变化为正值，其间导数值为 0 时相应横坐标值为吸收峰的峰值。当一阶导数值由正值变化为负值时，期间为 0 值的相应横坐标值为峰谷的位置。横坐标以波数为单位的投射谱从左到右的d总是负值，当 dT/ d由正值变为负值，经过 0 值时所对应的横坐标值为峰位，一阶导数由负变正经 0 值时所对应的横坐标为峰谷位置。横坐标以波长为单位的吸收谱图从左到右其波长值逐渐增加，所以 d总是正值，当一阶导数 dA/ d由正值变为负值其间经过 0 值时所对应的横坐标为峰位，一阶导数由负值变正值经过 0 值时相应的横坐标值为峰谷位置。横坐标以波长为单位的投射谱的 d总是正值。当一阶导数由负值变化为正值，经过导数值为 0 处相应的横坐标值为峰位；一阶导数由正值变化为负值，经过导数值为 0 处相应的横坐标为峰谷峰位。这些用一阶导数求峰值的方法对于连续函数是精确的，但是对于离散数据点组成的谱图就会有误差。数据点取得愈少，误差就愈大，即对于分辨率低的谱图，求出的峰值离真实值偏差越大。在油茶籽光谱处理中，计算机通过二次差值法找出峰值。对于一条油茶籽的光谱（采取平均采样间隔） ，用峰值附近的三个数据点的相应的强度 、 与 就可1iai1i14以计算出峰强与峰位：)2(811iiiMaa（3）iii4)(1（4）由式（3） （4）可见 i 是相应于 ai 处的波长值，计算机通过每次取三个数据点、逐步右移一个数据点的方法比较每次求出的 aM 值，最大的 aM 值为峰强，其相应的波数值 M 值为峰值（峰值位） 。由于计算机不能区分噪声与反射峰，因此需输入一个阀值，凡大于该阀值的峰才被标出，这样可以避免噪声的干扰。例如 CH3(CH2)7CH=CH(CH2)7COOH 油酸为油茶籽的主要成份 21。油酸 C18H34O2 即 3(18+34+2)-6=318 个波峰，所以针对油茶籽的高光谱拟合的系数数目不能太少。（2）原始高光谱数据波峰的计算。编写程序计算单个原始数据波峰数目： for k=0u=1:749%Y 中为油茶籽光谱的数据b1=Y(u,1); b2=Y(u+1,1);b3=Y(u+2,1);%计算峰值点if b2-b10&b3-b2k=k+1;endendkk =68（3）优化之后的波峰计算。Matlab 针对一个优化后的油茶籽高光谱数据进行波峰数目分析，得到的函数图像15为图 8。相应的 matlab 程序为： clear load data x=325:1075; y=mean(P);%P 为优化之后的光谱矩阵 size(y) clear bm=1;n=1;for k=1:747if y(1,k)y(1,k+2)&k-n5b(m)=y(1,k+1);n=k;m=m+1;plot(x(k+1),y(1,k+1),o);hold on;endend plot(x,y,r)图 8 两组油茶籽波长和波峰图16Fig. 8 Peak of two groups spectral data图 8 为用近红外光谱扫描出的两组油茶籽样品的优化后的光谱图波峰分布图（波峰用图中的“O”表示，横坐标为光谱范围，纵坐标为反射强度） 。从图中可见，在全光谱范围内，油菜籽在不同的波长段有不同的反射峰，这表明反射强度与所测定成分的含量成正比。因此，油茶籽的近红外光谱图可作为定量分析的依据。m=40 即有 m-1=39 个波峰，对应至少 39 种，优化后的 27 组光谱数据，其波峰数目的处理结果（表 4）表 4 优化之后的每组油茶籽的光谱波峰数据Table 4 Number of peaks after optimizing组编号 波峰数目 组编号 波峰数目 组编号 波峰数目1 53 10 48 19 432 53 11 45 20 513 47 12 42 21 404 56 13 40 22 525 55 14 35 23 526 48 15 46 24 577 48 16 47 25 438 51 17 54 26 549 56 18 54 27 48由表 4 可知，27 组高光谱的光谱图的波峰数目大致相同，其主要原因是因为油茶籽的品质不同，其成分种类相差很大，从而导致不同的波峰数目。故可以从样品的波峰数据，进行相同的优化处理，分析其波峰数，就可以将光谱进行降维。综上所述，波峰数目反映了平滑处理后的数据特性。其表示了每组光谱的曲线拟合之后的，一阶导数的取值为的数目（见附录图 1） ，即待降维的光谱的基本特性。油茶籽的近红外光谱图可作为定量分析的依据。4.3 光谱数据降维经典分析化学中利用标准曲线法进行浓度测定的过程仅以单点数据为测量点，如光谱以最大反射波长，而色谱则以面积等来表征仪器数据，这种方法失去很多有用数据。而且这种以单点测定的方法对油茶籽的矩阵型数据往往无能为力，对于较复杂的对象也往往是一筹莫展。另一方面，由于光谱数据变得越来越复杂，但是同时要求快速准确的定性，定量及结构分析的结果。多元回归分析（multivariate calibration）提供了现代分析仪器所提供的大量数据进行解析的方法。因此数学分析方法诸如主成份分析、因子分析、聚17类分析等引入，用于解决所面临的新问题。多元回归方法能解决如下问题：使分析对象数结构简化，用简单的方法形象地表示所研究的复杂对象，把所涉及的分析对象按其测量性质进行分类或对各种变量进行分组并进行变量之间相互关系的分析。另外提高信噪比、提高测量精密度、改善分析选择性、拓展应用范围也是多元回归包括的内容。4.3.1 PCA 降维主成分分析是设法将原来众多具有一定相关性（比如 P 个指标）重新组合成一组新的互相无关的综合指标来代替原来的指标。通常数学上的处理就是将原来 P 个指标作线性组合，作为新的综合指标。最经典的做法就是用 F1（选取的第一个线性组合，即第一个综合指标）的方差来表达，即 Var(F1)越大，表示 F1 包含的信息越多。因此在所有的线性组合中选取的 F1 应该是方差最大的，故称 F1 为第一主成分。如果第一主成分不足以代表原来 P 个指标的信息，再考虑选取 F2 即选第二个线性组合，为了有效地反映原来信息，F1 已有的信息就不需要再出现在 F2 中，用数学语言表达就是要求 cov(F1,F2)=0，则称 F2 为第二主成分，依此类推可以构造出第三、第四，第 P 个主成分。现在油茶籽的光谱数据矩阵 Xn p，n 表示测量的光谱数据即 3251075 nm,p 为测量的光谱组数目 27。对 Xn p 奇异值分解（singular value decomposition, SVD）这种矩阵算法用于数据阵 X 的分解可得X=USVT （5）由公式（5）可见：U 为 n p 阶行正交矩阵；V 为 p r 阶列正交矩阵，而 S 为r r 阶对角矩阵； r 为维数，它的最大值可为 n 或 p 中的较小值。对角阵 S 中的对角元素为正值，而非对角元素均为零。SVD 分解还可用矩阵 U 和 S 的列向量来表示：（6）r1iTkTrT2T1X由于矩阵 U 和 V 中列向量的正交性，即：UTU=VTV=Ir （7）其中 Ir 为 r r 阶单位矩阵，该式表示矩阵 U 和 V 各列元素平方和为 1，而列向量间的内积为零，向量的这种性质称为正交性。由以上所述的主成分分析方法得到的特征值个数（或称因子数） ，一般与数据矩阵 x 中行数和列数中较小的值相等，或小于这个数。若列数 p 小于行数 n，则有特征值个数， 。由于任何实际的化学测量都存在实验误差或任意一个数据矩阵总npr有重要因子和次要因子，因此不能认为任何从数学上求出的特征值个数就是重要因子18数或化学测量体系中存在的主要组分数或更一般地说主成分数。当然如果不存在实验误差或不存在次要因子，那么矩阵的特征值个数就是试样中的组成分或测量数据中的重要因子或矩阵的主成分数。实际中油茶籽的光谱矩阵 X 的响应信号来源于两方面：一是油茶籽中的主要组分（待测物质和共存干扰物质） ；二是测量误差、仪器噪声以及含量甚微的次要组分。在测量矩阵分解得到的特征值中，只有前面 r 个特征值是由于主要组分的存在引起的。按照前面对油茶籽的高光谱的奇异值分解（SVD）可将矩阵 X 进行分解求得 p 个特征值 ，在迭代算法中个特征值是按照由大到小的次序排列为：（8）p1rkr1kp1k由公式（8）可见中右边第一项是主要因子（主要成分）对特征值的贡献，后面一项是误差及其他微量成分等因素对特征值的贡献。前面已指出， 的大小与原始数据矩阵所进行预处理有关，因而一般不能直接依据特征值绝对值的大小来判别主因子数或组分数，而应该根据其相对大小及其变化规律加以判别。常用于判别组分数的判别指标有比率法（RATIO）和（RSD ）函数法。比率法从组分数等于 1 开始计算 r 取不同值时的 RATIO 值，当该组分百分数值大于某确定值（根据情况可取 85%95%）时的 r 值为主成分数。（9）%10)/(RATIOr1kpk上式（9）中仅表示列数 p 小于行数 n 的情况，否则可将累加号中的 p 以