基于多线性判别分析的玉米种子高光谱图像多特征转换方法

技术领域

本发明涉及一种玉米种子高光谱图像多特征转换方法，尤其是一种基于多线性判别分析的玉米种子高光谱图像多特征转换方法。

背景技术

玉米是世界总产量最高的粮食作物和主要的畜牧业饲料来源，同时也是工业生产中的重要原料，因而对于减少种子的混杂，保证农业生产的顺利进行具有重要的价值。目前，对玉米种子的外部品质检测主要采用可见光图像检测和红外图像检测等较成熟的技术，但是这些检测技术无法获取玉米种子内部品质有效信息。近红外光谱技术可在线实时无损检测种子的内部品质，而近红外光谱技术只提供对检验个体一个小区域的检测。由于种子的品质在空间上存在差异，所以该方法还存在一定的局限。高光谱图像技术是一种广泛应用的农产品品质无损检测方法，可以获得种子的几何形态学特征和光谱特征等优点，在种子品种识别中得到了广泛的应用。

种子的品种识别从本质上来说是一个模式识别问题，分类特征的充分挖掘获取是识别模型精度和鲁棒性的保证。尽管高光谱图像技术可获得种子的外部形态学特征和内部的化学成分特征，但是现有的高光谱图像识别中，多是利用单一特征(如光谱特征等)，存在着分类特征信息丢失的可能性。较充分的特征信息有利于分类模型的学习，从而提高模型的识别能力。然而，过多的特征学习也会存在着过拟合的风险，高光谱图像波段数目众多，同时波段之间含有大量的冗余信息，不仅影响识别的精度，而且不利于高光谱图像技术在种子品质检测在线应用的实时性，因而有必要对数据进行降维处理。

特征提取和波段选择是高光谱数据降维的两种主要途径。特征提取是将原始特征转换为一组具有明显物理意义、统计意义或核的特征降维方法，如主成分分析(Principal Component Analysis，PCA)、线性判别分析(Linear discriminant analysis，LDA)等。波段选择是寻找最能代表原始光谱信息波段的过程，以此减少和优化特征空间且不改变原始的光谱信息，如无信息变量消除法(Uninformative Variable Eliminate，UVE)、连续投影算法(Successive Projections Algorithm，SPA)等。特征提取和波段选择的降维方法已经被广泛应用到农产品品质的无损检测中。一方面，这些方法利用的是对于分类有利的特征，在包含更多对于分类的无用特征时，这些特征提取方法效果难以保证。另一方面，由于这些多特征波段选择方法多是先利用均值光谱特征进行波段选择，然后将选择的波段应用到其他特征，导致对其他特征信息考虑不足，不能确保选择的最优波段对于其他特征的适用性。在如何充分提取高光谱图像的有用信息，并同时满足高光谱图像无损检测技术在线应用的实时性和高分类精度上仍存在改进的空间。

发明内容

本发明的目的是克服现有技术中存在的不足，提供一种基于多线性判别分析的玉米种子高光谱图像多特征转换方法，其能够实现玉米种子高光谱图像多特征条件下特征选择，操作简单，快速有效，并具有较高的鲁棒性等优点。

按照本发明提供的技术方案，所述基于多线性判别分析的玉米种子高光谱图像多特征转换方法，其特征是，包括：

a、将待识别的Num粒玉米种子样本放置在高光谱图像采集系统中，采集并获取该M粒玉米种子样本在L个波段下的L个高光谱图像，L为自然数；

b、利用阈值分割获取玉米种子的轮廓曲线，提取轮廓下玉米种子的均值光谱特征和纹理特征，将L个波段下获得的共L个均值光谱特征和L个纹理特征进行特征融合，将融合特征作为待识别玉米种子的特征参数X；

c、利用步骤b所得到的玉米种子特征参数X，基于多线性判别分析(MLDA)算法对玉米种子特征参数X进行特征提取，并得到特征提取后的融合特征矩阵计算每个波段下的原始特征对第一主成分的贡献率W，设定阈值，选择较大W值所对应的波段作为最优波段，输出最优波段集合Θ，将最优波段集合Θ应用到特征提取后的融合特征矩阵得到特征转换后的特征S；

d、建立预测模型，对转换后的特征进行评价。

进一步的，在所述步骤b中，获得玉米种子的特征参数X的具体步骤包括：

玉米种子高光谱图像的均值光谱特征和纹理特征分别为

(其中，X¹表示均值光谱特征；X²～X¹⁴表示纹理特征，且X²～X⁶分别表示为纹理特征中一阶数据统计的：能量(energy)、熵(entropy)、偏度(skew)、标准差(standard>7～X¹⁴分别为纹理特征中二阶数据统计在0°、45°、90°和135°四个方向上的对比度均值(contrast>i^j为第i粒玉米种子的第j个特征，Num个样本在L个波段下构成一个L×Num的矩阵；首先对特征进行归一化处理，使融合特征达到数量级上的一致性；再对每一个样本在第p个波段下分别进行特征组合(为第j个特征)，得到14个特征在L个波段下Num个样本的融合特征矩阵X∈{K^{14×14×L×Num}}，其中，p＝1,2,...,L。

进一步的，在所述步骤c中，获得玉米种子的特征参数转换特征的具体步骤包括：

根据MLDA计算最佳投影向量将14×14×L×Num维度的原始融合特征空间投影到J₁×J₂×J₃×J₄的特征空间，得到投影后的特征矩阵其中，Y_m的计算公式为：

$Y_m=X\times {{}_{1}U}^{{\left(1\right)}^T}\times {{}_{2}U}^{{\left(2\right)}^T}\times {{}_{3}U}^{{\left(3\right)}^T}\times {{}_{4}U}^{{\left(4\right)}^T}---\left(1\right);$

不同n-mode的投影方向U⁽ⁿ⁾计算公式为：

$\begin{array}{c}{U}^{\left(n\right)}{|}_{n=1}^4=\underset{U^{\left(1\right)},U^{\left(2\right)},U^{\left(3\right)},U^{\left(4\right)}}{\arg \max }\frac{{\Sigma }_{c=1}^C{n}_c\left|\right|({\bar{X}}_c-\bar{X})\times {{}_{1}U}^{{\left(1\right)}^T}\times {{}_{2}U}^{{\left(2\right)}^T}\times {{}_{3}U}^{{\left(3\right)}^T}\times {{}_{4}U}^{{\left(4\right)}^T}|{|}^2}{{\Sigma }_{i=1}^{Num}{\Sigma }_{c=1}^C\left|\right|(X_i-{\bar{X}}_c)\times {{}_{1}U}^{{\left(1\right)}^T}\times {{}_{2}U}^{{\left(2\right)}^T}\times {{}_{3}U}^{{\left(3\right)}^T}\times {{}_{4}U}^{{\left(4\right)}^T}|{|}^2}\\ =\underset{U^{{\left(n\right)}^T}{U}^{\left(n\right)}=I}{\arg \max }\frac{trace(U^{{\left(n\right)}^T}·S_{B\left(n\right)}·U^{{\left(n\right)}^T})}{trace(U^{{\left(n\right)}^T}·S_{W\left(n\right)}·U^{{\left(n\right)}^T})}\end{array}---\left(2\right);$

公式(2)中，共有C类样本，n_c为第c类的样本数，是第c类样本融合特征的平均值，为全部样本融合特征的平均值，X_i表示第i个样本；与分别表示为按照l-mode张量展开的类内离散方差和类间离散方差，其中，U⁽³⁾设置为L×L的正交单位矩阵，U⁽⁴⁾为Num×Num的单位矩阵；

根据线性判别分析(Linear Discriminant Analysis,LDA)中主成分与原始特征之间满足线性关系t＝WXⁱ的思想，推广到MLDA中，所求取的第一主成分t₁与原始融合特征X∈{K^{14×14×L×Num}}也满足线性关系，U⁽³⁾的第一行可反映为每个波段下的特征信息对第一主成分的贡献率W，W对应的值越大，则该波段下的特征信息对主成分的贡献率越大，以此选出最优波段集合Θ；将最优波段集合应用到特征选择后的融合特征矩阵，得到特征转换后的矩阵

进一步的，所述步骤c中MLDA算法首先对融合特征进行特征提取，得到特征提取后的矩阵，根据每个波段下的特征信息在第一主成分上的贡献率作为波段选择的依据，将较大贡献率所对应的波段作为最优波段，并应用到特征提取后的矩阵得到转换特征，具体包括：

S1、计算整体样本均值第c类样本均值初始化投影矩阵根据公式(2)和公式(3)对融合特征进行特征选择；

S2、通过最优化迭代法得到第t次S_W(n)^(t)和S_B(n)^(t)的l-mode张量展开分别为：

$\begin{array}{c}{S_{B\left(n\right)}}^{\left(t\right)}={\Sigma }_{c=1}^C{n}_c{({\bar{X}}_c-\bar{X})}_{\left(n\right)}·{\Phi }_{\bar{n}}^{{(t-1)}^T}·{\Phi }_{\bar{n}}^{(t-1)}·{({\bar{X}}_c-\bar{X})}_{\left(n\right)}^T\\ S_{W\left(n\right)}={\Sigma }_{i=1}^{Num}{\Sigma }_{c=1}^C(X_i-{\bar{X}}_c)·{\Phi }_{\bar{n}}^T·{\Phi }_{\bar{n}}·{({\bar{X}}_c-\bar{X})}_{\left(n\right)}^T\end{array}---\left(3\right);$

S3、根据广义特征值分解计算(S_W(n)^(t))^-1S_B(n)^(t)＝UΛU^T中最大广义特征值所对应的单位广义特征向量

S4、当达到设定的阈值或最大迭代次数K时则停止迭代，输出最佳投影方向向量U⁽ⁿ⁾和特征选择后的融合特征矩阵否则返回步骤S2；

S5、计算每个波段下的原始特征对第一主成分的贡献率W，设定阈值，选择较大W值所对应的波段作为最优波段，输出最优波段集合Θ；

S6、将最优波段集合Θ应用到特征选择后的特征矩阵得到转换特征S。

本发明所述基于多线性判别分析的玉米种子高光谱图像多特征转换方法，通过对玉米种子高光谱图像进行特征转换，其能够实现多特征条件下特征转换，操作简单，快速有效，并具有较高鲁棒性等优点。

附图说明

图1为玉米种子高光谱图像的图像分割和特征提取过程示意图。

图2为本发明所述基于多线性判别分析的玉米种子高光谱图像多特征转换方法的流程示意图。

具体实施方式

下面结合具体附图对本发明作进一步说明。

本发明所述基于多线性判别分析的玉米种子高光谱图像多特征转换方法，包括：

a、如图1所示：首先选择待识别玉米种子轮廓最清晰的图像对应的波段(在700.1nm处)，利用自适应阈值分割法，获得该波段下的待识别玉米种子的轮廓曲线。将该轮廓曲线投射到L个波段上，提取L个波段下玉米种子轮廓图像的均值光谱特征和纹理特征(一阶统计量和二阶统计量)，L为自然数；

b、如图2所示：对均值光谱特征和纹理特征进行归一化处理，消除特征之间在数量级上的差异；

c、玉米种子高光谱图像的均值光谱特征和纹理特征分别为

(其中，X¹表示均值光谱特征；X²～X¹⁴表示纹理特征，且X²～X⁶分别表示为纹理特征中一阶数据统计的：能量(energy)、熵(entropy)、偏度(skew)、标准差(standard>7～X¹⁴分别为纹理特征中二阶数据统计在0°、45°、90°和135°四个方向上的对比度均值(contrast>14×14×L×Num}，其中，p＝1,2,...,L；

d、根据MLDA(多重线性判别分析，Multiple Linear Discrminant Analysis)寻找具有最大类间散度和最小类内散度(广义瑞利商)的投影方向的思想计算最佳投影向量最佳投影向量将14×14×L×Num维度的原始融合特征空间投影到J₁×J₂×J₃×J₄的特征空间，得到投影后的特征矩阵其中，Y_m的计算公式为：

$Y_m=X\times {{}_{1}U}^{{\left(1\right)}^T}\times {{}_{2}U}^{{\left(2\right)}^T}\times {{}_{3}U}^{{\left(3\right)}^T}\times {{}_{4}U}^{{\left(4\right)}^T}---\left(1\right);$

假设共有C类样本，n_c为第c类的样本数，是第c类样本融合特征的平均值，为全部样本融合特征的平均值，X_i表示第i个样本。不同n-mode的投影方向U⁽ⁿ⁾计算公式为：

$\begin{array}{c}{U}^{\left(n\right)}{|}_{n=1}^4=\underset{U^{\left(1\right)},U^{\left(2\right)},U^{\left(3\right)},U^{\left(4\right)}}{\arg \max }\frac{{\Sigma }_{c=1}^C{n}_c\left|\right|({\bar{X}}_c-\bar{X})\times {{}_{1}U}^{{\left(1\right)}^T}\times {{}_{2}U}^{{\left(2\right)}^T}\times {{}_{3}U}^{{\left(3\right)}^T}\times {{}_{4}U}^{{\left(4\right)}^T}|{|}^2}{{\Sigma }_{i=1}^{Num}{\Sigma }_{c=1}^C\left|\right|(X_i-{\bar{X}}_c)\times {{}_{1}U}^{{\left(1\right)}^T}\times {{}_{2}U}^{{\left(2\right)}^T}\times {{}_{3}U}^{{\left(3\right)}^T}\times {{}_{4}U}^{{\left(4\right)}^T}|{|}^2}\\ =\underset{U^{{\left(n\right)}^T}{U}^{\left(n\right)}=I}{\arg \max }\frac{trace(U^{{\left(n\right)}^T}·S_{B\left(n\right)}·U^{{\left(n\right)}^T})}{trace(U^{{\left(n\right)}^T}·S_{W\left(n\right)}·U^{{\left(n\right)}^T})}\end{array}---\left(2\right);$

公式(2)中，与分别表示为按照l-mode张量展开的类内离散方差和类间离散方差，

最佳投影向量U^(l)可以转化为计算矩阵S_W(l)^-1S_B(l)求取最大特征值所对应的特征向量的问题，通过k阶最优化迭代法来计算U^(l)。其中，U⁽³⁾设置为L×L的正交单位矩阵，U⁽⁴⁾为Num×Num的单位矩阵，从而实现对特征的降维而不改变波段的数目。

e、根据线性判别分析(Linear Discriminant Analysis,LDA)中主成分与原始特征之间满足线性关系t＝WXⁱ的思想，推广到MLDA中，所求取的第一主成分t₁与原始融合特征X∈{K^{14×14×L×Num}}也满足线性关系，U⁽³⁾的第一行可反映为每个波段下的特征信息对第一主成分的贡献率W，W对应的值越大，则该波段下的特征信息对主成分的贡献率越大，以此选出最优波段集合Θ；

f、将最优波段集合应用到特征选择后的融合特征矩阵，得到特征转换后的矩阵

g、对特征转换后的特征组合建立分类识别模型，根据分类精度进行评价。