高光谱亚像元分解预测花生中的黄曲霉毒素B1

韩仲志 刘 杰

（青岛农业大学理学与信息科学学院 山东青岛266109）

黄曲霉毒素（Aflatoxin）是一种剧毒、强致癌物，其毒性为砒霜的68 倍，是目前发现的最强的I类化学致癌物质。其广泛存在于花生、玉米籽粒表面及其制品中。中美两国都对其进行了强制限量：食品级和饲料级的限量标准分别为20 μg/kg 和100 μg/kg[1-2]。目前对黄曲霉毒素的检测主要是生化方法，包括薄层层析法、高效液相色谱法、微柱法、酶联免疫吸附法等[3]。虽然检测精度高，但检测手段和检测仪器复杂，检测速度慢，价格昂贵，不能在线检测。

近年来，光谱成像技术作为一项新的化学计量学手段，广泛应用于农产品及食品检测中[4-5]。黄曲霉毒素具有紫外荧光特性和表面浅表面分布特性。高光谱技术为一种图谱合一的新兴手段，在黄曲霉毒素检测中受到普遍关注。M.Atas 等[6]研究了市场上购买的辣椒粉黄曲霉毒素污染的高光谱成像分析技术，提出基于神经网络权值优化的特征选择方法，并找到对检测起关键作用的特征波长，在420 nm 处留一法识别率达到85%；H.Yao[7]团队主要集中在玉米的黄曲霉毒素污染的高光谱研究，他们通过人工种植的方法发现荧光强度随着黄曲霉素含量的增加而减少，并且存在峰移现象，还基于此两点发现申请了专利[8]，然而这种设备的检测还是手工进行，并不能做到自动检测。王伟[9]研究组使用美国农业部（USDA）的高光谱数据研究黄曲霉毒素的检测问题，指出根据普通CCD 成像，黄曲霉毒素自然污染籽粒的检出率为87.5%，他们采用Fisher 优选了700～1 100 nm的5 个高光谱波长，对人工污染黄曲霉毒素浓度的预测正确率达到88.3%[10]。

上述研究忽略了一个关键因素，黄曲霉毒素在籽粒表面分布并不均匀，这是由于黄曲霉毒素是黄曲霉菌的代谢物，呈颗粒状不均匀分布。通过光学探测器得到的像元往往是特定分辨率下各种物质平均像元，得到的光谱是各种物质（包括黄曲霉毒素）的混合光谱。单个像素内黄曲霉毒素丰度是定量计算整个籽粒含量的重要基础，有必要对单像素内的黄曲霉毒素丰度进行解析。本研究拟借用遥感领域像元分解的方法，探讨花生籽粒表面黄曲霉毒素的丰度，进而对黄曲霉毒素进行定量反演。

1 试验材料

1.1 样品制备

试验所用花生样品为市场上购买4 粒红小花生，经青岛市海润农大检测中心检测未检出黄曲霉毒素，挑选1 g 左右表面光滑250 籽粒备用。在青岛农业大学食品学院实验室对购买的黄曲霉毒素B1标准样品，通过乙氰配比20，50，100，500 和1 000 μg/kg 共5 种黄曲霉毒素溶液，然后使用移液器将1 μL 的黄曲霉毒素溶液滴到花生表面，自然风干。如表1，共制备4 粒红花生5 个浓度，每个浓度50 粒，共250 粒花生，对这250 个籽粒以此编号为1～250，正面朝上，放置在无荧光纸板上，冷藏备用。

表1 试验材料浓度与编号Table 1 Concentration and number of experimental materials

1.2 图像采集

图像采集在中国科学院青岛光电所光谱实验室的暗室环境下进行。图1（a）列出了试验中使用的采集设备，试验仪器为便携式成像仪，该成像仪主要由液晶可调式滤波器（Liquid Crystal Tunable Filter，LCTF）和一台普通CCD 相机组成，谱段范围为400～720 nm，光谱带宽为10 nm，使用一台UV 365 nm 大功率LED 紫外灯 （美国陆阳LUYOR-3404 台式紫外灯，样品处照度7 000 流明）作为光源。高光谱图像的分析采用ENVI4.7（ITT Visual Information Solutions，Boulder，Colo.），Matlab2012b（Math Works，USA）完成。计算机配置为：CPU 为Intel-E4600，主频2.39 GHz，内存为DDR3-3.25 GB。高光谱成像仪配套采集软件可根据图像整体亮度情况自动调整每个波长下的曝光时间从而避免过饱和和欠饱和现象。

紫外灯下，将花生籽粒污染面朝上每9 个一组，每次采集9 颗花生籽粒，依此采完250 个籽粒的图像，每次获得的高光谱图像数据块大小为1 392×1 040×33。图1（a）是基于RGB 图像合成原理，通过33，5，16 三个波段（720，440，550 nm）合成的假彩色图像。

图1 图像采集设备图Fig.1 Image acquisition equipment

图2 假彩色合成图像Fig.2 Composited color image

2 研究方法

2.1 算法

高光谱图像像素的混合过程通常通过线性混合模型 （linear spectral random mixture model，LSRMM）[11]来描述，高光谱某个混合像元的辐射值：

其中，ei，j为第j个纯像元的光谱曲线，j=1，2，…m，m 为纯像元个数；其中i=1，2，…L，L 为波段数；aj为第j 个像元的丰度，aj满足非负约束和为一约束。

根据线性混合模型，需要确定混合像元的端元波谱，N-FINDR 端元提取算法[12]是一种从已获得的高光谱图像自动化的获取端元的方法。高光谱的所有像素在高维空间形成一个凸体，每个端元对应于凸体的顶点，非端元则分布在凸体的内部、棱或面上，该算法的迭代思想是，对经PCA 或MNF 降维后的数据，随机选择一组像素作为初始端元，然后将图像中每个像素替换现有的端元向量中的端元计算替换后的体积，若体积增大则用新的端元向量替换原端元向量。重复遍历所有像元，最终的端元向量即是该图像的端元。

著名的科学杂志《Nature》刊登了两位科学家D.D.Lee 和H.S.Seung[13]对数学中非负矩阵研究的突出成果。该文提出了一种新的矩阵分解思想——非负矩阵分解（Non-negative Matrix Factorization，NMF）算法，即NMF 是在矩阵中所有元素均为非负数约束条件之下的矩阵分解方法。

丰度图像的灰度分布直方图可能提供一些有价值的信息。可用于黄曲霉素含量的预测中，对一副灰度图像的灰度量化特征表示为：

式中：n=1，2，…，N 表示灰度值量化为灰度阶的数量。I（x，y）——坐标（x，y）处灰度值。将这些直方图量化特征作为支持向量机回归的输入特征。

支持向量机（Support Vector Machine，SVM）是Corinna Cortes 和Vapnik 等[14]首先提出的，它在解决小样本、非线性及高维模式识别中表现出许多特有的优势，并能够推广应用到函数拟合等其它机器学习问题中。支持向量机算法可描述为：

式中：αk——拉格朗日乘子；K（x，xk）——核函数，本研究选自RBF 核函数；b--偏置系数。

2.2 研究框架

黄曲霉毒素在花生表面分布呈微小颗粒状不均匀分布，通过光学探测器得到的像元是一定区域（一个像素）内各种物质平均像元。将单个像元进行亚像元分解，分解为花生种皮光谱、黄曲霉素光谱和背景光谱各个组分的丰度，通过黄曲霉素组分丰度图像可更为精确地反演黄曲霉含量。图3是本文提出的方法流程图。

图3 方法流程图Fig.3 Method flow chart

对通过高光谱仪获得的高光谱图像，首先进行图像预处理和ROI 提取，将单个花生籽粒高光谱图像块提取出来，然后通过N-FINDR 算法，提取花生种皮、黄曲霉毒素及背景的端元光谱，反复比较提取较为纯净的3 个端元光谱，并观察3 个端元光谱波形，确定黄曲霉素的光谱，然后使用非负矩阵分解（NMF）方法对籽粒高光谱图像进行非负矩阵分解，得到黄曲霉毒素丰度图像，对丰度图像求取直方图量化特征，使用该特征进行回归预测输入参数，将该丰度与实际数据进行支持向量机回归，为比较预测效果，与偏最小二乘回归预测效果进行比较，进行精度评价和检测方法的有效性。

3 结果与分析

3.1 图像预处理

背景区域对图像分析造成干扰，需要去除，采用RGB 图像分量运算（R-1.2G）去除背景，并将背景置零；根据籽粒区域平均亮度一致原则将籽粒区域进行光照补偿，然后标记籽粒区域，将籽粒区域定义为感兴趣区域（ROI），根据籽粒左上角坐标与籽粒区域长宽，在33 个波段上将籽粒提取出来，顺序存储，将单籽粒高光谱数据存为.mat 文件，对应的伪彩色图像存为.jpg 文件，共得到250个单个籽粒（ROI）高光谱图像块，预处理过程通过Matlab 编程自动实现。图4（a）是分割得到其中一个单籽粒高光谱图像块ROI 对应的假彩色图像，可见图像上黄曲霉素区域的荧光现象比较明显。

3.2 N-FINDR 纯端元提取

端元提取就是在高光谱图像上找到纯像元（端元），图4（b，c）是基于N-FINDR 对图4（a）提取的端元所在位置图和端元曲线。图4（b）绿色圈点表示寻找到的端元的位置。图4（c）从端元曲线来看端元1 在400～420 波段内有一明显的波峰，这于黄曲霉素的荧光曲线相符合，可见端元1 即为黄曲霉素端元，端元3 为纯黑色很显然是背景的端元，因为预处理是已将背景置为0，那么端元2 则为花生种皮光谱端元。可以看出，基本将3 种端元寻找出来。根据经验，多个籽粒反复比较得到黄曲霉毒素、花生种皮和背景的纯像元光谱。

图4 端元及其光谱曲线Fig.4 Endmember and spectrum

由于高光谱数据量巨大，在端元提取之前首先要进行最大噪声分离（MNF）[15]以减少数据量，这样可以用前几个波段来代替整个高光谱影像，可以有效加快算法的执行效率。图5是最大噪声分离后的信噪比分布，可以看出，前面MNF 波段信噪比高，5 个以后信噪比已经低于1，前2 个MNF 分量的分布相当集中（图5b）。

图5 MNF 前两个端元分布与等高线图Fig.5 First two end-members distribution and contour map

3.3 NMF 丰度提取

使用非负矩阵分解（NMF）进行亚像元分解，就是得到各个端元的丰度，图6是各个端元的丰度，可以明显看出，图6（a）是黄曲霉毒素的丰度，图6（b）是花生种皮的丰度，图6（c）是背景的丰度，通过黄曲霉素丰度图（图6a）像素值求和即可计算出每个籽粒黄曲霉素的丰度和。

图6 各端元的丰度图Fig.6 Endmember abundance image

对黄曲霉毒素丰度图像进行求和。图7（a）是基于MNF 进行丰度分解后得到的250 个籽粒黄曲霉毒素的丰度和分布，单从图像上看并无规律。可见仅通过丰度求和较难准确预测黄曲霉毒素含量。丰度概率分布可能蕴含着大量的信息，对丰度图像进行直方图量化是体现丰度概率分布的有效手段。

3.4 直方图量化

直方图量化就是将丰度图像（0～1）转化为灰度图像（0～255），然后将灰度值直方图分成不同的区间（如0～50，50～100，100～150 等），统计灰度值落入每个区间的像素个数。图7（b）是对图6（a）黄曲霉毒素丰度图的直方图量化为12 个区间得到的直方图量化结果，通过量化可以发现，灰度主要集中在低灰度区域，这与整副图像偏黑有关。各个区域直方图分布代表了不同丰度黄曲霉毒素的分布，这些分布信息与黄曲霉毒素含量直接相关，可将该量化特征（12 个）作为回归模型的输入特征。

图7 黄曲霉毒素丰度统计及直方图量化Fig.7 Aflatoxin abundance and quantification histograms

3.5 支持向量机回归

将直方图量化特征作为黄曲霉毒素回归模型的输入参数，进行含量预测，与黄曲霉毒素的实际含量比较得到预测误差。图8是比较了支持向量机回归（SVR）与传统的偏最小二乘回归（PLS）得到的误差图，前200 个样本用来训练，后5 个样本用来测试，可以看出SVR 效果较PLS 效果好，且训练集误差小于测试集误差。

表2是采用五折交叉验证法得到的两种模型的预测结果，可以看出SVR 训练集误差只有0.89%。平均相对误差为12.16%，由于黄曲霉素限量很低，只有20～100 μg/kg，1 μg/kg 相当于1 t 里的1 g[15]，能够达到这样一个精度，可以在很大程度上将高黄曲霉毒素含量的籽粒鉴别出来，这是实际生产应用可以接受的。

图8 SVR 和PLS 两种方法预测误差Fig.8 Predicted error of SVR and PLS methods

表2 两种模型的预测结果Table 2 Predicted result of the two models

4 结论

环境中的灰尘、背景板在一定程度上都具有荧光现象，在高光谱伪彩色合成图像上表现为高亮黄绿荧光区域，这些因素的存在对黄曲霉素的准确预测造成了干扰，实用的技术中需要对这些因素进行预先的去除，比如加装必要的除尘设备和配置没有荧光现象的背景板，特氟龙材质的背景板是一种很好的选择。

基于亚像元分解的方法被广泛应用遥感影像的定量解析中，本研究基于MNF 的亚像元分解，提出了一种能在食品安全检测领域应用的黄曲霉丰度检测方法，并构建了基于丰度图像的直方图量化特征进行非线性回归，进而预测黄曲霉素精度，该方法最优总体预测误差以降低到12.16%。由于黄曲霉毒素含量在mg/kg 数量级，也就是百万分之一数量级，这一精度，已能够较为准确地识别花生是否被黄曲霉毒素污染，进而将污染籽粒剔除。相比较生化的检测方法，虽然预测精度较低，但该方法预测速度快，可以在瞬间完成，对黄曲霉的在线快速检测及其相关便携式装备的研发具有积极意义。