CN104198397A - 氮磷钾交互作用下温室作物营养含量检测的方法 - Google Patents

Abstract

本发明氮磷钾交互作用下番茄营养含量检测的方法，涉及温室作物营养含量亏缺诊断技术领域。建立包含氮磷钾交互作用的光谱表达式方程用于定量分析番茄营养胁迫的情况，在数学上表现为足够的组合变量和特征来同时反演作物几种养分的水平、以及养分之间的交互作用等胁迫信息，重点开发了在考虑氮、磷、钾三者交互作用下的番茄营养含量检测模型的建立，获得精度更高精度的检测模型。

Description

氮磷钾交互作用下温室作物营养含量检测的方法

技术领域

本发明涉及温室作物营养含量亏缺诊断技术领域；特指一种氮磷钾交互作用下诊断番茄营养含量的方法。

背景技术

番茄(LycopersiconesculentumMill.)是我国温室栽培的主要蔬菜作物之一。氮、磷、钾是番茄生长必需的也是最重要的营养元素，氮、磷、钾营养亏缺会使其生理发生变化，直接影响到产量的高低和口感，进而影响经济效益。植物生理学的研究已经证实，叶片是对营养状况反映最敏感的部位之一，营养亏缺能引起叶片颜色、纹理、粗糙度和气孔等特征发生显著变化，因此，以作物叶片为研究对象诊断营养状况成为近年来研究的热点。

国内外科研学者在影响偏振反射光谱监测植物对象特性的单独因素方面开展研究的同时，还考虑到各个因素之间的主次顺序和参数优选。试验因素都不是孤立的，对目标结果的影响具有主次顺序和轻重关系。有研究表明，基于高光谱技术和偏振光谱技术的多信息融合并结合恰当的建模方法可以较大程度的提高模型的精度。但是氮、磷、钾之间是否存在互作关系，即发生氮元素营养胁迫的情况下引起了缺磷或者缺钾的状态，以此类推，在钾元素或者磷元素发生营养胁迫时，是否间接引起了另外两种元素的亏缺。近几年来国内外一些学者主要是将高光谱技术、偏振光谱技术应用于农产品质量以及作物病害的检测中，但未见利用高光谱图像技术、偏振光谱技术来诊断作物氮磷钾交互作用的检测。

发明内容

本发明建立包含氮磷钾交互作用的光谱表达式方程用于定量分析番茄营养胁迫的情况，在数学上表现为足够的组合变量和特征来同时反演作物几种养分的水平、以及养分之间的交互作用等胁迫信息，重点开发了在考虑氮、磷、钾三者交互作用下的番茄营养含量检测模型的建立，以期获得精度更高精度的检测模型。

与以往单独建立某一种元素亏缺的预测模型相比较，综合考虑氮、钾、磷交互作用情况下的营养水平诊断。通过分析氮、磷与钾之间的交互关系，消除互作之间的影响。同时反演作物几种养分的水平以及养分之间的交互作用等胁迫信息时，光谱变量间有谐波关系，在数学上表现为需要足够的组合变量，即从机理上讲方程之间是不相互独立的，从而通过若干特征波长的光谱反射率组合求解和反演得出的氮、磷、钾营养信息对于不同品种、不同生育期及不同营养条件下的检测模型。

本发明氮磷钾交互作用下温室作物营养含量检测的方法，按照下述步骤进行：

(一)、建立氮(N)、磷(P)、钾(K)营养胁迫试验样本，每个营养元素分为五个水平进行处理，

(二)、分别建立基于图像特征、光谱特征、偏振度特征的单一特征模型，

(三)、定义氮磷钾交互作用的方程式表达，

(四)、权重系数矩阵的计算，

(五)、确定交互影响系数矩阵，

(六)、交互模型的建立，用该模型检测作物是否发生氮磷钾亏缺。

其中步骤(一)中每个营养元素分为不同水平进行处理，是指按照在正常配方中氮、磷、钾的正常含量的25％～150％(质量比)。

其中步骤(二)中所述的分别建立基于图像特征、光谱特征、偏振度特征的单一特征模型，按照下述步骤进行：(1)偏振光谱采集，指利用偏振光谱采集系统采集温室番茄叶片的偏振光谱。

(2)偏振度特征提取，根据斯托克斯公式计算步骤(1)中番茄叶片的的偏振光谱的偏振度。

(3)高光谱图像采集，利用高光谱图像采集系统采集番茄叶片的高光谱图像。

(4)图像预处理，对步骤(3)中的高光谱图像进行波段筛选、滤波和利用掩膜对图像进行分割。

(5)图像纹理特征的提取，对步骤(4)中经预处理后的高光谱图像，首先通过主成分分析得到氮、磷、钾的敏感波长，然后在敏感波长下进行基于二阶概率统计滤波的纹理特征提取。

(6)光谱特征的提取，对步骤(5)中提取到的纹理特征，利用区间偏最小二乘法-遗传算法优选特征变量，对光谱特征变量进行提取。

(7)模型建立，采用支持向量机建立番茄开花期氮、磷、钾营养含量预测模型，具体按照下述步骤进行：(a)对样本进行归一化预处理，确定输入特征的数量；

(b)分别基于网格搜索法(GS)和粒子群优化算法(PSO)对进行支持向量机回归(SVR)参数进行寻优；

(c)将偏振度特征、纹理特征和光谱特征采用SVR的特征层融合的方法建立番茄开花期氮、磷、钾营养含量的预测模型。

其中步骤(三)中所述的定义氮磷钾交互作用的方程式表达是指利用提取的光谱特征、图像特征、偏振度特征建立交互作用下的氮磷钾检测模型。

其中步骤(四)中所述的权重系数矩阵的计算是指对建立的交互作用下的氮磷钾检测模型中的权重系数矩阵进行求解。

其中步骤(五)中所述的确定交互影响系数矩阵是指对建立的交互作用下的氮磷钾检测模型中的交互影响系数矩阵进行求解。

其中步骤(六)中所述的交互模型的建立是指将权重系数矩阵和交互影响系数矩阵代入交互作用下的氮磷钾检测模型。

本发明所述的温室作物优选番茄、生菜、黄瓜等。

本发明的有益效果：

本发明通过权重系数矩阵和交互影响系数矩阵对模型进一步修正，研究建立了在考虑氮、磷、钾三者交互作用下的番茄营养含量检测模型，预测结果表明该方法能够获得精度更高的检测模型，N、P、K元素的交互模型比单一特征类型模型精度分别平均提高了4.3％，9.0％，3.3％，模型的预测值更接近作物营养含量的真值，为交互作用下作物营养快速无损检测提供新的方法和思路。

附图说明

图1.高光谱成像系统，

其中：1.光箱；2.光源；3.控制器；4.计算机；5.近红外相机；6.成像光谱仪；7.步进电机；8.玻璃光纤灯。

图2.偏振光谱采集系统，

其中：1.光源；2.光源光纤灯；3.探测光纤；4.旋转测量机构；5.样品台；6.光谱仪；7.微电流计；8.计算机。

图3.番茄开花期氮素交互模型实测值和预测值相关图。

图4.番茄开花期磷素交互模型实测值和预测值相关图。

图5.番茄开花期钾素交互模型实测值和预测值相关图。

图6.番茄结果中期氮素交互模型实测值和预测值相关图。

图7.番茄结果中期磷素交互模型实测值和预测值相关图。

图8.番茄结果中期钾素交互模型实测值和预测值相关图。

具体实施方式

下面以番茄为例，结合附图对本发明进行进一步详细描述。

本发明具体实施方式中所采用的高光谱图像采集系统参阅图1。利用图1所示的高光谱图像采集系统采集温室番茄叶片高光谱图像，其包括近红外相机5(XEVA-FPA-1.7-320，XenICs，Leuven，Belgium)，光谱范围900-1700nm，成像光谱仪6(ImspectorN17E，Spectral ImagingLtd.，Finland)，分辨率为5nm，150W卤钨灯的直流可调光源2(2900-ER+9596-E，Illumination Technologies，Inc.，EastSyracuse，NY，USA)，位移单元由步进电机7(MTS120，北京光学仪器厂，北京，中国)和控制器3(SC100，北京光学仪器厂，北京，中国)组成，成像光谱仪可以采集图像传输到计算机4(DELLInspiron 530s，USA)中，玻璃光纤灯8为图像的采集提供必要的照明。近红外相机5、成像光谱仪6、步进电机7和玻璃光纤灯8位于光箱1。本发明2012年3月至2012年9月在江苏大学现代农业装备与技术教育部重点实验室的Venlo型温室中进行实验。培育品种为辽宁省农业科学院蔬菜研究所培育的L-402。为保证前期的基础性研究能够对番茄分的有效特征进行准确提取，本发明采用无土栽培技术进行样本培育。在保证其他营养元素均衡的情况下，对氮、磷、钾进行精确控制，以获取纯正的营养胁迫样本。营养液pH值为6-6.5、EC值为1.2ms/cm。首次栽植后浇灌正常营养液，为避免珍珠岩中营养成分残留，在番茄生长至壮苗期转入缺素培育阶段时，进行二次移栽。

营养胁迫试验样本分为氮(N)、磷(P)、钾(K)三组，每组中各个营养元素均分为五个水平进行处理，按照在正常配方中氮、磷、钾的正常含量的25％～150％，分别形成五种不同营养水平的样本，依次是重度胁迫25％、中度胁迫50％、轻度胁迫75％、适量100％、过量150％。(均为质量比)

化学值含量的测定工作与光谱试验同步进行，培育的样本在编号的自封袋中排序，并放进专业植物保鲜箱，立即带回光谱实验室，开始偏振反射光谱测量实验和高光谱图像采集。采集结束后将叶片放入烘箱，80℃烘干至恒重置于干燥器内。采用凯氏定氮法(GB/T 5009.5-1985)测定样本的全氮含量，仪器为英国SEAL公司生产的Auto Analyzer 3型连续流动分析仪(Seal Analytical Instruments Co.，Ltd，England)。采用钼锑抗分光光度法(GB 11893-1989)测定样本中磷的含量，仪器为美国瓦里安公司紫外可见分光光度计(Varian Inc.，Palo Alto，USA；Model Cary 100)。采用火焰光度法(GB/T 18633-2002)测定样本中钾的含量，仪器为BWB-XP多元素火焰光度计(BWB Co.，British)。化学值含量测定后用于下面的模型建立及校正。

偏振光谱采集：

偏振反射光谱测量分析系统为课题组自主研发，该仪器测量波长范围350-1000nm，如图2所示，光源1接出的光源光纤灯2安装在左侧的旋转测量机构4，光谱仪6接出的探测光纤3安装在右侧的旋转测量机构4，样品台5用于放置待测样品，微电流计7分别连接光谱仪6及操作电脑8。

偏振反射光谱和高光谱图像采集保证在尽量短的时间内依次快速进行，以保证样品的统一性。数据采集前，两套仪器均需进行预热以及黑场和白场标定，以消除环境因素引起的系统误差，每片叶片测量3次，取平均值作为最终测量结果。光谱实验样本为氮、磷、钾每组各120个，采集时间均为早晨8点，并统一选择每株倒七叶，叶片的叶宽应大于2cm，一经采下迅速装入自封袋封好并进行现场编号，放进专业植物保鲜箱，立即带回光谱实验室，开始偏振反射光谱测量实验和高光谱图像采集。实验中为防止外界环境光的干扰，偏振反射光谱测量在暗室中进行；高光谱图像则在光箱中采集。测量结束后将叶片放入烘箱，80℃烘干至恒重置于干燥器内以备化学值测定之用。

高光谱图像采集：

高光谱图像数据的采集是基于SpectralCube(Spectral Imaging Ltd.，Finland)软件平台；实际采集的光谱范围为871.6～1766.3nm，空间分辨率为62.5um，采样间隔为3.5nm，一次采集可获取采样光谱范围内以3.5nm为间隔的256幅独立的高光谱图像。

确定近红外相机的曝光时间以保证图像的清晰，同时确定位移台的速度以避免图像尺寸和空间分辨率的失真。经过分析比较确定曝光时间为为20ms，位移台的移动速度为1.25mm/s。数据采集时，首先进行黑场和白场标定，设定反射率范围，进而利用二阶巴特沃茨滤波器进行数字滤波，去除噪声干扰。

图像预处理：

为了避免传感器暗电流以及光源的照度在各波段下分布不均而使得图像含有较大的噪声，或是造成不同波长下的较大亮度值差异，先对所原始样本图像进行标定。全白的标定图像W通过扫描硫酸钡标准白板得到；全黑的标定图像B则是通过盖上摄像机的镜头后采集得到。按照公式(1)，将采集得到的绝对图像I变成相对图像R。

$R=\frac{I-B}{W-B}---\left(1\right)$

采集的番茄叶片高光谱的光谱曲线的范围是390～1050nm。光谱曲线在450nm以下和950nm以上区域存在着明显的噪声，因此在后期的数据处理过程中，选取450～950nm范围内，共388个波段的高光谱数据进行后续分析研究。

本发明采用5×5窗口的中值滤波法对图像进行滤波，5×5窗口的中值滤波能够既降噪又不失真，对后续特征提取更有利。将在700nm处将图像分割得到的二值图像建立掩膜，掩膜是ENVI中一种特殊的图像，是一个由0和1组成的二进制图像。当有掩膜参与到一个高光谱图像处理时，1值区域被处理，0值区域被屏蔽，即背景为黑色，掩膜的部分将不参与后续的计算处理，白色区域为分割出来的番茄叶片区域，参与到后续处理中，极大的减少了背景对特征提取的影响。

图像纹理特征的提取：

主成分分析法(PCA)是高光谱图像降维中最常用的方法。它的目标是寻求一种变换，把原始数据映射到一个新空间。分别求出经过PCA转换后氮、磷、钾的前五个主成分的权重系数曲线图，绘制出的氮、磷、钾高光谱前五个主成分的权重系数曲线图，根据提取权重系数曲线主要波峰和波谷处所对应的波长作为敏感波长。通过比对后其中N素的特征波长分别为：464.91nm，566.29nm，696.28nm，724.66nm；P素的特征波长分别为：474.85nm，567.54nm，693.71nm，738.89nm；K素的特征波长分别为：565.03nm，691.14nm，733.71nm，766.14nm。

在特征波长下采用基于二阶概率统计滤波的纹理特征提取中值、协方差、同质性、熵、相异性、二阶矩、对比度和相关性。

(1)中值(Mean)：

这里的中值运算类似于卷积，但计算的不是加权求和，而是把邻域中的像素按灰度级进行排序，然后选择该组的中间值作为输出，其主要功能是让与周围像素灰度值的差比较大的像素改取与周围的像素值接近的值。

(2)协方差(Variance)： $\mathrm{VAR}=\underset{i=0}{\overset{L-1}{\mathrm{\Σ}}}\underset{j=0}{\overset{L-1}{\mathrm{\Σ}}}(i-u_x)(j-u_y)p(i,j)---\left(3\right)$

式中：

一般情况下协方差直观上表示的是两个变量总体误差的期望。这里协方差用于衡量两个变量的总体误差。

(3)同质性(Homogeneity)：

又称逆差距，它能够度量图像纹理局部变化的多少。其值大则说明图像纹理在不同区域间变化较少，局部非常均匀，值小说明局部分布不均匀。

(4)熵：

图像所具有的信息量的度量。它反映了图像的混乱程度和无序程度，代表了图像中纹理的非均匀性和复杂性，熵值越大纹理越复杂，熵值越小则纹理越均匀。

(5)相异性(Dissimilarity)： $\mathrm{DIS}=\underset{i=0}{\overset{L-1}{\mathrm{\Σ}}}\underset{j=0}{\overset{L-1}{\mathrm{\Σ}}}\frac{p(i,j)}{1+{|i-j|}^2}---\left(6\right)$

相异性大小用于考量元素在行的方向或列的方向上所有元素之间的相异程度，改指标是相对于其他行或者其他列而言的，如果图像的某方向上差异性较大，则该方向的DIS值将大于其他方向的DIS值。

(6)二阶矩(AngularSecondMoment)：

又称为能量，它是通过灰度共生矩阵元素值平方和的计算，反映图像灰度分布的纹理粒子间粗细的程度。二阶矩的值越高说明纹理较粗，而二阶矩的值越低说明纹理较细。

(7)对比度(Contrast)：

改指标对比了纹理沟纹的深浅程度及图像清晰程度。假如沟纹越深，CON值越大，视觉效果越清晰且图像越明显；CON值越小，则表明沟纹浅，效果模糊。

(8)相关(Correlation)：

式中：

该指标用于衡量灰度共生矩阵中的元素在行方向上，或者列方向上的相似度。它能够反映图像局部灰度的相关性。如果图像的行方向上或者列方向上的COR值越高说明元素的相似程度越高，反之则越低。

然后将统计出的这八个特征参数分别与番茄叶片N、K、P含量的实测值进行相关性分析，结果见表1，表2和表3。

表1 敏感波长下纹理特征与番茄N含量的相关系数

表2 敏感波长下纹理特征与番茄P含量的相关系数

表3 敏感波长下纹理特征与番茄K含量的相关系数

四个敏感波长下纹理特征的中值特征与N、P、K三者含量的相关性都较弱，说明其是无效特征，应舍弃。根据表1，表2和表3，优选相关性高的特征作为用于模型建立的特征变量，为了方便后续的研究，从相关系数中找出六个与氮、磷、钾的相关性都较高的特征，作为氮磷钾的共同高光谱特征，分别为：VAR_693.71、CON_566.29、DIS_693.71、ENT_733.71、ASM_566.29、COR_733.71。再从剩余的特征中分别选择两个相关性最高的特征作为氮磷钾的特有的高光谱特征，依次分别为：氮ASM_464.91、COR_464.91；磷HOM_693.71、ENT_474.85；钾HOM_762.24、ENT_762.24。

建立基于图像特征的模型：

本研究分别采用MLR、PCR和PLS建模方法建立番茄N素胁迫叶片高光谱的模型，模型对应的校正和预测的结果如表4所示。表4中，Rc表示校正模型的相关系数；Rp表示预测模型的相关系数；RMSECV和RMSEP分别交互验证均方根误差和预测均方根误差；MRE表示预测模型的平均相对误差(下同)。分别采用MLR、PCR和PLS建模方法建立番茄P素胁迫叶片高光谱的模型，模型对应的校正和预测的结果如表5所示。分别采用MLR、PCR和PLS建模方法建立番茄K素胁迫叶片高光谱的模型，模型对应的校正和预测的结果如表6所示。

表4 基于高光谱图像特征的N元素模型校正和预测结果

表5 基于高光谱图像特征的P元素模型校正和预测结果

表6 基于高光谱图像特征的K元素模型校正和预测结果

光谱特征的提取：

遗传算法(GA)模拟了对生物界自然遗传机制和自然选择，希望通过模拟来解决变量间的优化问题。GA比较适合用于解决的复杂的、非线性优化问题，而传统搜索方法在解决此类问题时效果欠佳，遗传算法目前已经被广泛用于优选近红外光谱特征波数点。将区域偏最小二乘优选特征谱区和遗传算法优选特征波数点的思想相结合，采用间隔偏最小二乘法-遗传算法(iPLS-GA)优选光谱的特征谱区。首先用iPLS截选出建模精度最佳的波段，再采用GA的方法优选出能够代表此波段的若干个变量后再进行建模，希望通过尽量少的特征变量代替全谱数据同时又能获得较好的模型精度。

iPLS-GA模型在iPLS模型的基础上优选出4个变量，即N元素的敏感波长，分别是741.48nm，755.74nm，767.44nm，784.37nm。iPLS-GA模型在iPLS模型的基础上优选出8个变量，即P元素的敏感波长，分别是770.04nm，779.16nm，813.12nm，824.92nm。优选出4个变量，即K元素的敏感波长，分别是618.23nm，630.97nm，645.00nm，705.30nm。基于光谱特征的氮磷钾iPLS-GA模型的结果如表7所示。

表7 基于光谱特征的氮磷钾模型结果

偏振度计算：

偏振度P为全偏振分量的强度与该光波总强度的比值，由此证明光束的偏振度P可以通过公式(10)求得：

在自然大气背景及目标物体对光源入射的偏振效应中，圆偏振的分量极少，圆偏振分量在仪器可以检测的范围内更小，相对于误差来说可以忽略。对于番茄叶片表面的反射光而言，其圆偏振分量对结果的影响予以忽略，故其斯托克斯矢量中的V分量假设约为0，仅用到了I，Q，U三个参量，此时的偏振度可以写成：

因而，只需测出光线在0°、90°、+45°、-45°方向上的线偏振分量的光强，就可以通过公式(11)完全确定一束光线的偏振状态，根据上一节的研究内容，在最优的组合条件下，分别提取探测器处偏振片在0°、90°、+45°、-45°方向上每个番茄样本叶片区域的平均偏振光谱。对不同氮磷钾营养水平番茄叶片样本的偏振度数据逐点进行相关性分析，获得350～1000nm波段范围番茄的N、P、K含量的相关性曲线，根据相关分析的结果选择与N、P、K含量相关性高的偏振度特征。

由于相邻波段的偏振度特征存在极强的相关性，假设从相邻波长随意选择必定存在重复，使得选出的特征波长或者特征波长的组合不能具备代表性，因此需从相关性极其显著的敏感波段内选取相关性较弱的波长作为敏感波长，也就是说敏感波长必须满足以下条件：(1)在相关分析得出的波段范围内；(2)确保选出来的波段之间具有弱相关性；(3)为了满足后续研究的需要，选出N、P、K共有的敏感波长，尽量能够选出对N、P、K都比较敏感的波长。具体实施步骤为：首先将相关分析得出的敏感波段分成若干个小波段，每个小波段含有10个波长；然后根据相关系数的大小将小波段进行排序，组成波段子集Ui；将相关系数最大的波段放到已选择波段集合U_s；按照前后顺序依次从小波段子集Ui中选择不在U_s中的波段，如果它满足与U_s中所有波段的相关系数R都小于0.8，则放入U_s中，否则予以剔除。根据筛选条件，选出最能代表偏振度与番茄叶片N、P、K含量有显著关系的特征波长，其中655.41nm，744.48nm，850.58nm为N、P、K共有敏感波长；而N、P、K特有的敏感波长分别为380.49nm，914.56nm，488.42nm。

建立基于偏振度特征的模型：

表8、表9、表10分别为氮、磷、钾元素的采用MLR、PCR、PLS建立的基于偏振度模型的预测结果。

表8 番茄叶片N含量偏振度模型校正和预测结果

表9 番茄叶片P含量偏振度模型校正和预测结果

表10 番茄叶片K含量偏振度模型校正和预测结果

(1)定义氮磷钾交互作用的光谱方程式的表达：

设{X}为多信息特征向量，即{X}＝{T_g，T_t，T_p}^T，其中T_g，T_t，T_p分别为光特征、图像特征、偏振度特征，角标g、t、p分别代表光谱、图像、偏振度(下同)；[A]为权重系数矩阵；{Y}为氮、磷、钾检测值向量，即：{Y}＝{T_j，T_j，T_j}^T，则：式中的系数矩阵如果看作是变量的权重，即[A]矩阵表示对{X}向量的权重系数矩阵。那么一般线性模型可以写成(12)式：

{Y}＝[A]{X}   (12)

其中： $\left[\begin{array}{c}A\end{array}\right]=\left[\begin{array}{ccc}{a}_{11}& a_{12}& a_{13}\\ a_{21}& a_{22}& a_{23}\\ a_{31}& a_{32}& a_{33}\end{array}\right]$

式中：a₁₁、a₁₂、a₁₃表示对氮素而言强度光谱特征、图像特征、偏振度特征分别在模型建立过程中所占的比重；

a₂₁、a₂₂、a₂₃表示对磷素而言强度光谱特征、图像特征、偏振度特征分别在模型建立过程中所占的比重；

a₃₁、a₃₂、a₃₃表示对钾素而言强度光谱特征、图像特征、偏振度特征分别在模型建立过程中所占的比重。

将氮磷钾的检测值代入向量{Y}＝{T_j，T_j，T_j}^T，使得{Y}＝{N_j，P_j，K_j}^T，则(12)式可以写成：

$\left[\begin{array}{c}{N}_j\\ P_j\\ K_j\end{array}\right]=\left[\begin{array}{ccc}{a}_{11}& a_{12}& a_{13}\\ a_{21}& a_{22}& a_{23}\\ a_{31}& a_{32}& a_{33}\end{array}\right]\left[\begin{array}{c}{T}_g\\ T_t\\ T_p\end{array}\right]---\left(13\right)$

由(13)式可知，通过等号右边的光谱、图像、偏振度的相应特征值作为输入，计算可以得到等式左边的N、P、K的检测值。而(13)式也能够看作是涵盖所有的线性融合模型的通用方程。在实际的植物生长过程中，往往不只是仅发生缺某一种元素的情况，大多数情况下会同时受到两种或三种主要元素的营养胁迫。也就是说作物缺氮同时可能引起缺磷或者缺钾等。当我们考虑氮磷钾的两两交互作用时，我们在(13)式的基础上引入系数矩阵[B]，将矩阵[B]定义为：

$\left[\begin{array}{c}B\end{array}\right]=\left[\begin{array}{ccc}{b}_{11}& b_{12}& b_{13}\\ b_{21}& b_{22}& b_{23}\\ b_{31}& b_{32}& b_{33}\end{array}\right]$

式中：b₁₁表示氮对氮的影响，b₁₂表示磷对氮的影响，b₁₃表示钾对氮的影响；

b₂₁表示氮对磷的影响，b₂₂表示磷对磷的影响，b₂₃表示钾对磷的影响；

b₃₁表示氮对钾的影响，b₃₂表示磷对钾的影响，b₃₃表示钾对钾的影响。

有了[B]矩阵对N、P、K之间的交互作用影响作进一步修正，使之成为营养元素含量的真实值和检测值之间的一个修正矩阵。因此，营养元素的真实值和检测值这两者与[B]矩阵之间的关系可以用(14)式来进行表示(其中N_z，P_z，K_z分别表示氮、磷、钾的真实值)：

$\left[\begin{array}{c}{N}_j\\ P_j\\ K_j\end{array}\right]=\left[\begin{array}{ccc}{b}_{11}& b_{12}& b_{13}\\ b_{21}& b_{22}& b_{23}\\ b_{31}& b_{32}& b_{33}\end{array}\right]\left[\begin{array}{c}{N}_z\\ P_z\\ K_z\end{array}\right]---\left(14\right)$

因此通过(13)式和(14)式的等价交换可以推导出：

$\left[\begin{array}{ccc}{b}_{11}& b_{12}& b_{13}\\ b_{21}& b_{22}& b_{23}\\ b_{31}& b_{32}& b_{33}\end{array}\right]\left[\begin{array}{c}{N}_z\\ P_z\\ K_z\end{array}\right]=\left[\begin{array}{ccc}{a}_{11}& a_{12}& a_{13}\\ a_{21}& a_{22}& a_{23}\\ a_{31}& a_{32}& a_{33}\end{array}\right]\left[\begin{array}{c}{T}_g\\ T_t\\ T_p\end{array}\right]---\left(15\right)$

将(15)式中[B]矩阵移到等号右边，转化得到(16)式：

$\left[\begin{array}{c}{N}_z\\ P_z\\ K_z\end{array}\right]={\left[\begin{array}{ccc}{b}_{11}& b_{12}& b_{13}\\ b_{21}& b_{22}& b_{23}\\ b_{31}& b_{32}& b_{33}\end{array}\right]}^{-1}\left[\begin{array}{ccc}{a}_{11}& a_{12}& a_{13}\\ a_{21}& a_{22}& a_{23}\\ a_{31}& a_{32}& a_{33}\end{array}\right]\left[\begin{array}{c}{T}_g\\ T_t\\ T_p\end{array}\right]---\left(16\right)$

在(16)式等号右边的[A]、[B]矩阵均已求得，而分别代表光谱、图像、偏振度特征的特征向量{T_g，T_t，T_p}^T的对应取值也由实验获取的前提下，由此可以求出等式左边的{N_z，P_z，K_z}^T的真实值。

(16)式即为考虑氮磷钾交互作用下番茄营养含量的预测模型，它较好的描述了在考虑交互作用下的真实值和仪器检测值之间满足的关系，涵盖了光谱、图像和偏振度的特征向量{X}＝{T_g，T_t，T_p}^T各自的比重，在考虑三类不同来源的特征的权重的基础上，引入了交互影响系数矩阵[B]作为单独缺素模型的一个有效的补充，方程不仅考虑了仅某一种元素营养胁迫的情况，还考虑了N、P、K元素之间交互作用影响的情况，有望进一步提高常规方法所建立的融合模型的预测精度。

(2)权重系数矩阵的计算：

阐明了氮磷钾交互作用下的番茄营养含量的预测模型方程式表达的定义，接下来就是要对权重系数矩阵[A]和交互影响系数矩阵[B]进行求解。本节先计算权重系数矩阵[A]，在求之前先要知道权重的概念，所谓权重，是指某指标在整体评价中的相对重要程度。权重越大则该指标的重要性越高，对整体的影响就越高。每个指标的权重在0～1之间，并且满足所有指标之和为1。

本发明中需要求解的权重是针对某一种元素而言，其强度光谱特征、图像特征、偏振度特征这三种不同来源的特征被同时用于建立融合模型的过程中分别所占的比重；而这就必须考虑到这三种不同来源的特征在单独建立模型时的精度和误差。简而言之，假设基于图像特征的番茄叶片的N元素预测模型相对于基于强度光谱特征、偏振度特征建立的预测模型的精度更高，且误差更小，那么建立融合模型时，图像特征理应获得更高的权重。

分别基于强度光谱特征、图像特征和偏振度特征建立了N、P、K元素的预测模型，获得了相应的最佳模型。这里尝试根据以最佳模型的预测相关系数Rp和预测均方根误差RMSEP为考量来计算三种不同来源的特征的权重系数，理由是Rp和RMSEP是评价预测模型精度的两个最主要的指标，也是最重要的指标；并且这两个指标中，Rp是反映预测值和实测值之间相关关系的密切程度，RMSEP是对一组测量中的特小或者特大的误差反映十分敏感，它能够很好地反映出测量的精密度。因此，将这两种指标作为求解权重的考量标准是最具有说服力的。

下面进行权重系数矩阵[A]的求解，具体权重系数统计和求解过程如下：

第一步将N、P、K基于强度光谱特征T_g、图像特征T_t和偏振度特征T_p建立的最佳模型的Rp和RMSEP的值进行汇总，汇总结果见表11。

表11 最佳单一特征来源模型预测结果汇总表

第二步进行各特征的权重的计算。Rp和RMSEP在被用于评价模型的预测精度和稳定性时，衡量标准刚好相反，Rp的值是越接近1越好，说明预测值与实测值之间的相关程度越好；而RMSEP的值却是越小越好，越小的值越说明模型对外部样本的预测能力越高。基于此，在计算权重之前，为了便于计算，先要将这两者的进行数学上的统一，即用(1-RMSEP)后的值与Rp的值以各占一半的比例进行计算。例如N元素的光谱特征的权重的具体计算公式为：

$W_{\mathrm{Ng}}=\frac{{\mathrm{Rp}}_g\×0.5+(1+{\mathrm{RMSEP}}_g)\×0.5}{[{\mathrm{Rp}}_g\×0.5+(1-{\mathrm{RMSEP}}_g)\×0.5]+[{\mathrm{Rp}}_t\×0.5+(1-{\mathrm{RMSEP}}_t)\×0.5]+[{\mathrm{Rp}}_p\×0.5+(1-{\mathrm{RMSEP}}_p)\×0.5]}---\left(17\right)$

式中：W_Ng表示N元素的光谱特征的权重；

Rp_g、Rp_t、Rp_p分别表示基于光谱、图像、偏振度特征建立的最佳模型的预测相关系数；

RMSEP_g、RMSEP_t、RMSEP_p分别表示基于光谱、图像、偏振度特征建立的最佳模型的预测均方根误差。

以此类推，各元素的各类特征的权重均可以参考(17)式进行计算，例如计算N元素的强度光谱特征T_g的权重：

特征T_g的权重＝[0.9048*0.5+(1-0.6320)*0.5]/{[0.9048*0.5+(1-0.6320)*0.5]

+[0.9205*0.5+(1-0.4486)*0.5]+[0.9145*0.5+(1-0.7299)*0.5]}

＝0.6364/1.9647＝0.3239；

由此方法逐一计算后得出[A]矩阵的值为：

$\left[\begin{array}{c}A\end{array}\right]=\left[\begin{array}{ccc}{a}_{11}& a_{12}& a_{13}\\ a_{21}& a_{22}& a_{23}\\ a_{31}& a_{32}& a_{33}\end{array}\right]=\left[\begin{array}{ccc}0.3239& 0.3746& 0.3015\\ 0.4059& 0.3961& 0.1980\\ 0.3042& 0.3830& 0.3128\end{array}\right]$

从求得的[A]矩阵的结果可知，对N元素而言，图像特征的权重相对略高，P元素的光谱特征权重最高，K元素的图像特征权重最高。

(3)确定交互影响系数矩阵：

当已经求得[A]矩阵时，上一节中的(13)式的线性方程形式可以写成：

N_j＝a₁₁f_(Ng)+a₁₂f_(Nt)+a₁₃f_(Np)

P_j＝a₂₁f_(Pg)+a₂₂f_(Pt)+a₂₃f_(Pp)   (18)

K_j＝a₃₁f_(Kg)+a₃₂f_(Kt)+a₃₃f_(Kp)

这里(13)式中特征向量T_g，T_t，T_p分别写成了(18)中的光谱、图像和偏振度特征的线性模型的方程式的形式。(18)式等号右边的矩阵[A]已经计算得出，而通过实验，样本的光谱、图像、偏振度特征的相应的特征变量取值也已经统计出，代入线性融合模型的方程式f，可以计算出具体的某一种来源特征的检测值，如此等号右边均为已知，那么等号左边的氮磷钾三种不同特征变量信息来源融合后的检测值{N_j，P_j，K_j}^T就可以计算求得。

上一节中的(14)式的方程形式可以写成

N_j＝b₁₁Nz+b₁₂Pz+b₁₃Kz

P_j＝b₂₁Nz+b₂₂Pz+b₂₃Kz   (19)

K_j＝b₃₁Nz+b₃₂Pz+b₃₃Kz

将(18)式的计算结果代入(19)式，(19)式右边的化学方法测定的真实值也由化学实验得出，因此可以计算得出矩阵[B]。而这里的[B]矩阵的值根据生长周期的不同，并不是固定不变的。

下面就以番茄开花期和结果中期的数据样本为例进行[B]矩阵的求解。其中代入的真实值为开花期样本的平均值。

将开花期的数据为例进行计算，先联列方程组求得[B]矩阵第一行：

3.51b₁₁+4.34b₁₂+4.70b₁₃＝3.11

5.72b₁₁+2.83b₁₂+4.66b₁₃＝5.32

5.85b₁₁+4.28b₁₂+2.26b₁₃＝5.75

解出b₁₁＝1.0096；b₁₂＝0.0216；b₁₃＝-0.1120。接着联列个方程组求得[B]矩阵第二行和第三行。

由此得出开花期的交互系数矩阵[B]的值为：

$\left[\begin{array}{c}B\end{array}\right]=\left[\begin{array}{ccc}{b}_{11}& b_{12}& b_{13}\\ b_{21}& b_{22}& b_{23}\\ b_{31}& b_{32}& b_{33}\end{array}\right]=\left[\begin{array}{ccc}1.0096& 0.0216& -0.1120\\ -0.0858& 1.1468& -0.0304\\ -0.0418& 0.0792& 1.0049\end{array}\right]$

同理，将结果中期的数据为例联列方程组进行计算，得出结果中期的交互系数矩阵[B]的值为：

$\left[\begin{array}{c}B\end{array}\right]=\left[\begin{array}{ccc}{b}_{11}& b_{12}& b_{13}\\ b_{21}& b_{22}& b_{23}\\ b_{31}& b_{32}& b_{33}\end{array}\right]=\left[\begin{array}{ccc}1.0108& 0.0407& -0.0926\\ -0.0765& 1.1205& -0.0364\\ -0.0391& 0.0662& 1.0087\end{array}\right]$

从开花期和结果中期求得的两组[B]矩阵的值可以发现：①矩阵中的数值有正数和负数，正数表明对植株的营养元素有增强的作用，负数则表明对营养元素具有抑制作用。②在矩阵对角线上表示氮对氮影响的b₁₁，表示磷对磷影响的b₂₂和表示钾对钾影响的b₃₃数值均在1左右，且均大于1的正数，说明元素对自身的影响均为增强作用。③两组[B]矩阵中出现负号的数值位置一致，即发生抑制作用的元素是一致的，氮和磷之间互为抑制作用，氮和钾对磷均为抑制作用；而磷对氮和磷对钾均为增强作用。④从数值高低也可以看出氮对磷影响的数值均为两组数值中最高，说明氮对磷的影响比较明显。

(4)交互模型的建立，用该模型检测作物是否发生氮磷钾亏缺。

下面将番茄开花期和结果中期对应数据求得的[B]矩阵分别代入(16)式进行模型的验证，以预测相关系数Rp、预测均方根误差RMSEP和平均相对误差MRE为预测模型的评价指标。

(1)首先开花期模型的验证，将对应数据求得的[B]矩阵代入(16)式得到：

$\left[\begin{array}{c}{N}_y\\ P_y\\ K_y\end{array}\right]={\left[\begin{array}{ccc}1.0096& 0.0216& -0.1120\\ -0.0858& 1.1468& -0.0304\\ -0.0418& 0.0792& 1.0049\end{array}\right]}^{-1}\left[\begin{array}{ccc}0.3239& 0.3746& 0.3015\\ 0.4059& 0.3961& 0.1980\\ 0.3042& 0.3830& 0.3128\end{array}\right]\left[\begin{array}{c}{T}_g\\ T_t\\ T_p\end{array}\right]---\left(20\right)$

采用24个番茄叶片样本进行预测，将样本的光谱、图像、偏振的特征值代入(9)式，计算后得出N、P、K的预测值与真实值的预测相关系数Rp、预测均方根误差RMSEP和平均相对误差MRE如下表12所示，N、P、K预测结果的相关图分别见图3、图4、图5。

表12 开花期交互模型预测结果汇总表

由表12可知，采用开花期的交互模型进行预测，其中N元素的Rp＝0.9585，RMSEP＝0.2609％，MRE＝5.596％；P元素的Rp＝0.9201，RMSEP＝0.1739％，MRE＝5.346％；K元素的Rp＝0.9194，RMSEP＝0.2263％，MRE＝9.244％。N元素的预测相关系数最高，而P元素的预测均方根误差和平均相对误差最小，但N、P、K的预测相关系数均达到0.9以上，证明了通过权重系数矩阵[A]和交互影响系数矩阵[B]建立考虑三元素交互作用下的预测模型是可行的。

(2)下面进行结果中期模型的验证，将对应数据求得的[B]矩阵代入(16)式得到：

$\left[\begin{array}{c}{N}_y\\ P_y\\ K_y\end{array}\right]={\left[\begin{array}{ccc}1.0108& 0.0407& -0.0926\\ -0.0765& 1.1205& -0.0364\\ -0.0391& 0.0662& 1.0087\end{array}\right]}^{-1}\left[\begin{array}{ccc}0.3239& 0.3746& 0.3015\\ 0.4059& 0.3961& 0.1980\\ 0.3042& 0.3830& 0.3128\end{array}\right]\left[\begin{array}{c}{T}_g\\ T_t\\ T_p\end{array}\right]---\left(21\right)$

采用24个番茄叶片样本进行预测，将样本的光谱、图像、偏振的特征值代入(21)式，计算后得出N、P、K的预测值与真实值的预测相关系数Rp、预测均方根误差RMSEP和平均相对误差MRE如下表13所示，N、P、K预测结果的相关图分别见图6、7、8。

由表13可知，采用结果中期的交互模型进行预测，其中N元素的Rp＝0.9461，RMSEP＝0.2452％，MRE＝6.910％；P元素的Rp＝0.9183，RMSEP＝0.2616％，MRE＝8.387％；K元素的Rp＝0.9144，RMSEP＝0.2436％，MRE＝10.287％。结果中期的N、P、K的预测相关系数也均达到了0.9以上。

表13 结果中期交互模型预测结果汇总表

由表12和表13的交互模型的预测结果证明了通过权重系数矩阵和交互影响系数矩阵建立考虑三元素交互作用下的预测模型的可行性。将单一特征类型模型、信息融合模型和交互模型的预测结果进行汇总，结果见表14。

表14 不同建模方法模型预测结果汇总表

由表14可以看出将三元素的交互作用模型与单一信息来源的模型分别进行比较，N元素的基于图像特征的检测模型相对与基于光谱特征、偏振度特征的检测模型而言精度较高，预测集的Rp＝0.9205，RMSEP＝0.4486％；P元素的基于光谱特征的检测模型相对其他两种特征的模型精度较高，预测集的Rp＝0.8740，RMSEP＝0.5120％；K元素的基于偏振度特征的检测模型精度较高，预测集的Rp＝0.9009，RMSEP＝0.7982％，而从开花期和结果中期数据建立的交互模型的验证结果来看，无论是N、P还是K元素，预测精度和模型的稳定性均有提高。而与采用BP-ANN以及SVR方法建立非线性融合模型相比较，N元素的BP-ANN模型预测集Rp＝0.9400，RMSEP＝0.1995％；P元素的SVR-PSO模型预测集Rp＝0.8998，RMSEP＝0.1912％；K元素的SVR-GS模型预测集Rp＝0.9101，RMSEP＝0.1417％，三种元素的考虑交互作用的预测模型的预测精度和模型的稳定性也有一定程度的提高。通过计算N、P、K元素的信息融合模型比单一特征类型模型的精度分别平均提高了2.1％，3.5％，0.9％，而N、P、K元素的交互模型比单一特征类型模型精度分别平均提高了4.3％，9.0％，3.3％，其中P元素提升幅度相对较大。

在前人研究某种营养元素胁迫的过程中，通常情况是假设这种元素亏缺，且其他元素不缺，在这样的前提下建立检测模型，而这也就无形中设定了该元素与其他元素不发生任何互作的情形。而实际的作物生长过程中这种情况是比较少的，通常元素之间存在着交互作用，[B]矩阵的求解证明了三种元素之间确实存在着交互作用和影响，交互模型的预测结果也表明通过[B]矩阵的求解对N、P、K交互作用的综合解耦，对模型作进一步修正可以提高多信息融合模型的预测精度。与以往单独建立某一种元素亏缺的预测模型相比较，该方法不仅融入了光谱、图像、偏振度三种不同信息来的特征，实现优势互补，还综合考虑了N、P、K交互作用情况来展开营养的定量诊断。通过交互影响系数矩阵的求解对番茄N、P、K之间的互作情况进行了量化，将交互影响矩阵作为N、P、K诊断模型的修正因子，消除互作之间的影响，同时反演作物几种营养成分的含量以及营养之间的交互作用的信息，这样模型的预测值更接近作物营养含量的真值。该方法能有效提高数学模型的预测精度，为交互作用下作物营养快速无损检测建模提供新的思路。

Claims (7)

1.氮磷钾交互作用下温室作物营养含量检测的方法，其特征在于按照下述步骤进行：

（一）、建立氮（N）、磷（P）、钾（K）营养胁迫试验样本，每个营养元素分为五个水平进行处理，

（二）、分别建立基于图像特征、光谱特征、偏振度特征的单一特征模型，

（三）、定义氮磷钾交互作用的方程式表达，

（四）、权重系数矩阵的计算，

（五）、确定交互影响系数矩阵，

（六）、交互模型的建立，用该模型检测作物是否发生氮磷钾亏缺。

2.根据权利要求1所述的氮磷钾交互作用下番茄营养含量检测的方法，其特征在于其中步骤（一）中每个营养元素分为不同水平进行处理，是指按照在正常配方中氮、磷、钾的正常含量的25%~150％（质量比），

根据权利要求1所述的氮磷钾交互作用下番茄营养含量检测的方法，其特征在于其中步骤（二）中所述的建立基于图像特征、光谱特征、偏振度特征的单一特征模型，按照下述步骤进行：

（1）偏振光谱采集，指利用偏振光谱采集系统采集温室番茄叶片的偏振光谱；

（2）偏振度特征提取，根据斯托克斯公式计算步骤（1）中番茄叶片的的偏振光谱的偏振度；

（3）高光谱图像采集，利用高光谱图像采集系统采集番茄叶片的高光谱图像；

（4）图像预处理，对步骤（3）中的高光谱图像进行波段筛选、滤波和利用掩膜对图像进行分割；

（5）图像纹理特征的提取，对步骤（4）中经预处理后的高光谱图像，首先通过主成分分析得到氮、磷、钾的敏感波长，然后在敏感波长下进行基于二阶概率统计滤波的纹理特征提取；

（6）光谱特征的提取，对步骤（5）中提取到的纹理特征，利用区间偏最小二乘法-遗传算法优选特征变量，对光谱特征变量进行提取；

（7）模型建立，采用支持向量机建立番茄开花期氮、磷、钾营养含量预测模型，具体按照下述步骤进行：（a）对样本进行归一化预处理，确定输入特征的数量；

（b）分别基于网格搜索法（GS）和粒子群优化算法（PSO）对进行支持向量机回归（SVR）参数进行寻优；

（c）将偏振度特征、纹理特征和光谱特征采用SVR的特征层融合的方法建立番茄开花期氮、磷、钾营养含量的预测模型。

3.根据权利要求1所述的氮磷钾交互作用下番茄营养含量检测的方法，其特征在于其中步骤（三）中所述的定义氮磷钾交互作用的光谱方程式的表达是指利用提取的光谱特征、图像特征、偏振度特征建立交互作用下的氮磷钾检测模型。

4.根据权利要求1所述的氮磷钾交互作用下番茄营养含量检测的方法，其特征在于其中步骤（四）中所述的权重系数矩阵的计算是指对建立的交互作用下的氮磷钾检测模型中的权重系数矩阵进行求解。

5.根据权利要求1所述的氮磷钾交互作用下番茄营养含量检测的方法，其特征在于其中步骤（五）中所述的确定交互影响系数矩阵是指对建立的交互作用下的氮磷钾检测模型中的交互影响系数矩阵进行求解。

6.根据权利要求1所述的氮磷钾交互作用下番茄营养含量检测的方法，其特征在于其中步骤（六）中所述的交互模型的建立是指将权重系数矩阵和交互影响系数矩阵代入交互作用下的氮磷钾检测模型。

7.根据权利要求1所述的氮磷钾交互作用下番茄营养含量检测的方法，其特征在于所述的温室作物为番茄。