扩展形态学与正交子空间投影结合的端元自动提取方法

摘要

一种扩展形态学与正交子空间投影结合的端元自动提取方法，步骤如下：(1)高光谱数据读入；(2)确定结构元素尺寸、初始与最大迭代次数；(3)由扩展膨胀、腐蚀操作计算结构元素邻域内混合度最低和最高的像元；(4)由步骤(3)得到的结果计算得到形态离心率指数值；(5)步骤(3)、(4)随迭代次数的增加不断重复，并利用步骤(3)扩展膨胀操作结果更新图像数据，直到达到最大迭代次数；(6)二值化MEI图像，得到端元数据集；(7)采用光谱角匹配方法计算得到第一个端元，并向得到端元的正交子空间投影更新端元数据集；(8)重复步骤(7)，直到满足误差要求。该方法是一种稳定性强、可靠性高、精确度高的高光谱端元自动提取方法。

说明书

技术领域

本发明涉及扩展形态学与正交子空间投影结合的端元自动提取方法，属于高光谱数据处理方法与应用技术领域，适用于高光谱数据端元提取的理论方法以及在农作物精细分类、矿物蚀变信息提取等领域的应用技术研究。

背景技术

高光谱数据处理与传统的遥感数据处理的一个本质的区别在于后者在对地物分类时通常是将单一像元根据某种规则分到某一不同的种类中，从而更加准确地反映遥感对象的状况。混合像元分解不仅是提高地物识别精度的重要方法，更是定量探测地物成分的有效途径，因此，混合像元分解是高光谱遥感数据应用的主要方向之一，该类技术的发展能够大大推动高光谱数据的应用范围与水平，并不断扩展高光谱数据的应用深度和广度。

目前端元提取的方法主要包括纯像元指数(Pure Pixel Index，PPI)、最小体积变换、迭代误差分析、顶点分量分析等，上述方法均只利用了高光谱数据提供的光谱/特征空间信息，针对这一问题，国外提出了基于数学形态学的方法，该方法同时利用了高光谱数据提供的光谱和空间信息，提高了方法的可靠性。但是该类方法存在着一定程度的问题：采用光谱角制图的方法度量光谱之间的相似度，受噪声影响较大；不同类型的端元采用了大津法，对于光谱相似的不同类型端元很难区分；并且算法的性能直接受结构元素形状与尺寸的影响。

发明内容

本发明的目的是针对现有端元提取方法从数据光谱或空间或特征信息的单一方面出发进行混合像元分解、不同类型端元难以区分等不足，提出一种扩展形态学与正交子空间投影结合的端元自动提取方法。

本发明的技术解决方案是：利用扩展形态学操作与正交子空间投影理论实现高光谱数据端元自动提取的方法，该方法主要是基于扩展数学形态学的理论，利用扩展膨胀和腐蚀操作，并通过计算形态离心率指数进行高光谱数据的端元数据集计算，然后利用光谱角匹配方法提取不同类型的端元，并通过向端元正交子空间投影消除已经提取端元的影响，实现高光谱数据所有端元的自动提取。形态偏心率指数是利用扩展膨胀和腐蚀的结果计算数据的正交投影散度，从而将图像中端元提取出来，引入正交子空间投影方法实现不同类型端元的有效区分，并且在本发明中为了避免由于形态学中结构元素的大小对方法性能的限制和降低计算过程中引入的人为误差的影响，采用迭代的方法，利用扩展膨胀操作不断更新参与计算的图像数据来实现。

本发明扩展形态学与正交子空间投影结合的端元自动提取方法，其步骤如下：

(1)高光谱遥感数据读入；

(2)根据高光谱数据特点，确定结构元素的尺寸、初始迭代次数以及最大迭代次数；

(3)由扩展数学形态学膨胀和腐蚀操作计算每一个结构元素邻域内混合度最低和最高的像元；

(4)由步骤(3)得到的混合度最低和最高的像元结果计算得到形态离心率指数值；

(5)步骤(3)、(4)随迭代次数的增加不断重复，并利用步骤(3)得到的混合度最低的结果不断更新图像数据，直到达到最大的迭代次数；

(6)二值化形态离心率指数MEI图像，得到端元数据集；

(7)利用由步骤(6)得到的端元数据集中任一元素作为参考光谱，采用光谱角匹配方法计算得到参考光谱与端元数据集中所有元素光谱角最大的元素，作为第一个端元，并向第一个端元的正交子空间投影，更新端元数据集；

(8)不断重复步骤(7)，直到满足误差要求，从而得到高光谱遥感图像中所有的端元。

其中，步骤(2)中所述的结构元素尺寸为5×5，初始迭代次数t＝1，最大迭代次数t_max＝5；

其中，步骤(3)所述的“由扩展数学形态学膨胀和腐蚀操作计算每一个结构元素邻域内混合度最低和最高的像元”，其扩展数学形态学膨胀和腐蚀操作定义如下：

$d(x,y)=(f\oplus b)(x,y)=\underset{(s,t)\in D_s}{\arg \_\mathrm{Max}}\left\{D(f(x-s,y-t),b)\right\}$

$e(x,y)=(f\otimes b)(x,y)=\underset{(s,t)\in D_s}{\arg \_\mathrm{Min}}\left\{D(f(x+s,y+t),b)\right\}$

其中，d(x，y)、e(x，y)分别表示像元(x，y)在结构元素b(x，y)内的扩展膨胀、腐蚀操作结果；f(x，y)是图像像元(x，y)的值，b(x，y)是结构元素；arg_Max，arg_Mn分别表示使得D达到最大和最小的像元向量；D是为了确定像元混合度的多维向量的排序关系，引入的一个多维向量的度量算子。该度量算子由结构元素内各个像元累加距离计算得到，定义如下：

$D(f(x,y),b)=\underset{s}{\Sigma }\underset{t}{\Sigma }\mathrm{dist}(f(x,y),f(s,t)),$(s，t)∈D_b

其中，(s，t)∈D_b表示(s，t)是属于结构元素b(x，y)内的点，dist是测量N维向量的逐点线性距离，采用正交投影散度(Orthogonal Projection Divergence，OPD)计算该距离，考虑两个N维光谱信号s_i＝[s_i1，s_i2，...，s_iN]^T，s_j＝[s_j1，s_j2，...，s_jN]^T，则N维光谱信号s_i和s_j之间的正交投影散度OPD表示为：

$\mathrm{OPD}(s_i,s_j)={(s_i^T{P}_{s_j}^{\perp }{s}_i+s_j^T{P}_{s_i}^{\perp }{s}_j)}^{1/2}$

$P_{s_k}^{\perp }=I_{N\times N}-s_k{\left(s_k^T{s}_k\right)}^{-1}{s}_k^T,$k＝i，j，并且I_N×N是N＊N维的单位矩阵。

因此，累加距离D能够根据像元的混合程度大小排序结构元素中的向量。通过以上的分析表明，扩展到高光谱图像的膨胀结果得到的是在结构元素内纯度最高的像元，腐蚀结果得到的是在结构元素内混合度最高的像元。

其中，步骤(4)中所述的“由步骤(3)得到的混合度最低和最高的像元结果计算得到形态离心率指数值”，其含义说明如下：通过正交投影散度融合结构元素邻域内膨胀和腐蚀结果，计算得到形态离心率指数MEI，计算公式如下：

MEI(x，y)＝OPD[d(x，y)，e(x，y)]

其中，步骤(5)中所述的“步骤(3)、(4)随迭代次数的增加不断重复，并利用步骤(3)得到的混合度最低的结果不断更新图像数据，直到达到最大的迭代次数”，其含义说明如下：设f_t(x，y)为第t次迭代计算得到的像元点(x，y)的图像数据，f_t-1(x，y)为第i-1次迭代计算得到的(x，y)的图像数据，更新图像数据的目的是避免人为引入误差或结构元素尺寸对提取端元纯度的影响，因此采用扩展膨胀操作实现图像数据更新，即$f_t(x,y)=(f_{t-1}\oplus b)(x,y).$

其中，步骤(6)中所述的“计算迭代后的形态离心率指数MEI图像，并进行二值化，得到端元数据集”，其含义说明如下：利用[m-3σ，m+3σ]进行MEI图像的二值化，得到所有端元的空间位置，其中，m为MEI图像均值，σ为MEI图像的标准差。

其中，步骤(7)中所述的“利用由步骤(6)得到的端元数据集中任一元素作为参考光谱，采用光谱角匹配方法计算得到参考光谱与端元数据集中所有元素光谱角最大的元素，作为第一个端元，并向第一个端元的正交子空间投影，更新端元数据集”，其含义说明如下：基于线性混合模型，利用正交子空间投影(Orthogonal Subspace Projection，OSP)概念进行混合像元分解时，仍可以继续将端元矩阵分解，则混合模型可以表示为：

r＝dα_d+Uα_U+n

其中r为端元数据集，d为得到的端元光谱信号，α_d为端元光谱信号的混合百分比，U为剩余端元光谱信号，α_U为剩余端元光谱信号的混合百分比，n为噪声信号；通过MEI图像中提取的端元数据，利用光谱角匹配方法得到其中一个端元光谱信号d＝e₁，然后利用向其正交子空间投影得到来消除端元数据中的信号d＝e₁，其中$P_d^{\perp }=I-d{\left(d^{T}d\right)}^{-1}{d}^T.$将代入基于OSP的混合模型中得到$P_d^{\perp }r=P_d^{\perp }{\mathrm{U\alpha }}_U+P_d^{\perp }n,$可以看出，不仅消除了d＝e₁的影响而且有效的压制了噪声然后利用更新端元数据集。

其中，步骤(8)中所述的“不断重复步骤(7)，直到满足误差要求，从而得到图像中所有端元”，其含义说明如下：经过m-1次投影之后得到d＝{e₁，e₂，...，e_m}，投影循环结束的条件由分解误差决定，从而得到了m个端元。

本发明与现有技术相比的优点在于：本方法同时利用了扩展形态学与正交子空间投影的理论，从数学机理上有效的克服了现有端元提取方法无法将不同类型端元有效区分的问题，实现了无先验信息情况下高光谱数据端元自动提取。它具有以下的优点：(1)通过扩展形态学膨胀和腐蚀操作实现图像数据中纯度最高和混合度最高的像元提取，并利用正交投影散度融合扩展膨胀和腐蚀结果得到形态离心率指数，实现了高光谱数据端元提取，有效地抑制了背景等因素对端元提取结果的影响，提高了方法的稳定性；(2)采用迭代的方法，利用扩展膨胀操作不断更新图像数据来克服由于形态学中结构元素的尺寸等对方法性能的限制；(3)引入正交子空间投影和光谱角匹配，实现不同类型端元的有效区分，提高了方法的精确性。

附图说明

图1为本发明扩展形态学与正交子空间投影结合的端元自动提取方法的流程图。

具体实施方式

为了更好的说明本发明涉及的扩展形态学与正交子空间投影结合的端元自动提取方法，利用星载Hyperion高光谱成像仪数据，采用本发明方法进行西藏驱龙地区蚀变矿物端元提取。本发明一种扩展形态学与正交子空间投影结合的端元自动提取方法，具体实现步骤如下：

(1)高光谱遥感数据读入：读入Hyperion高光谱数据，数据维数为255×550×196；

(2)根据高光谱数据特点以及应用目的，确定结构元素的尺寸和初始迭代次数：结构元素尺寸为5×5，初始迭代次数t＝1，最大迭代次数t_max＝5；

(3)由扩展数学形态学膨胀和腐蚀操作计算每一个结构元素邻域内混合度最低和最高的像元；

为了实现更加可靠、稳定、精确的高光谱数据端元自动提取，本发明中的方法同时利用了光谱和空间信息，由此引入了扩展数学形态学膨胀和腐蚀操作，定义如下：

$d(x,y)=(f\oplus b)(x,y)=\underset{(s,t)\in D_s}{\arg \_\mathrm{Max}}\left\{D(f(x-s,y-t),b)\right\}$

$e(x,y)=(f\otimes b)(x,y)=\underset{(s,t)\in D_s}{\arg \_\mathrm{Min}}\left\{D(f(x+s,y+t),b)\right\}$

其中，d(x，y)、e(x，y)分别表示像元(x，y)在结构元素b(x，y)内的扩展膨胀、腐蚀操作结果；f(x，y)是图像像元(x，y)的反射率值，b(x，y)是结构元素；arg_Max，arg_Min分别表示使得D达到最大和最小的像元向量；D是为了确定像元混合度的多维向量的排序关系，引入的一个高维向量的度量算子。该度量算子由结构元素内各个像元累加距离计算得到，定义如下：

$D(f(x,y),b)=\underset{s}{\Sigma }\underset{t}{\Sigma }\mathrm{dist}(f(x,y),f(s,t)),$(s，t)∈D_b

其中，(s，t)∈D_b表示(s，t)是属于结构元素b(x，y)内的点，dist是测量N维向量的逐点线性距离。为了有效地利用高光谱数据提供的光谱和空间信息，本发明采用正交投影散度OPD计算该距离，考虑两个196维光谱信号s_i＝[s_i1，s_i2，...，s_i196]^T，s_j＝[s_j1，s_j2，...，s_j196]^T，则196维光谱信号s_i和s_j之间的正交投影散度OPD表示为：

$\mathrm{OPD}(s_i,s_j)={(s_i^T{P}_{s_j}^{\perp }{s}_i+s_j^T{P}_{s_i}^{\perp }{s}_j)}^{1/2}$

$P_{s_k}^{\perp }=I_{196\times 196}-s_k{\left(s_k^T{s}_k\right)}^{-1}{s}_k^T,$k＝i，j，并且I_196×196是196＊196维的单位矩阵。

因此，累加距离D能够根据像元的混合程度大小排序结构元素中的向量。通过以上的分析表明，扩展到高光谱图像的膨胀结果得到的是在结构元素内纯度最高的像元，腐蚀结果得到的是在结构元素内混合度最高的像元。

(4)由步骤(3)得到的混合度最低和最高的像元结果计算得到形态离心率指数值；

为了克服利用单一膨胀或腐蚀操作人为引入的噪声和背景影响，以及为了定量化图像数据像元的纯度与混合度，采用了形态离心率指数。通过正交投影散度融合结构元素邻域内膨胀和腐蚀结果，计算得到形态离心率指数MEI，实现高光谱数据的端元提取。

MEI(x，y)＝OPD[d(x，y)，e(x，y)]

(5)步骤(3)、(4)随迭代次数的增加不断重复，并利用步骤(3)得到的混合度最低的结果不断更新图像数据，直到达到最大的迭代次数；

为了提取图像数据中纯度最高的像元，利用能够得到结构元素内纯度最高像元的扩展膨胀操作不断更新图像数据，从而有效地提高端元提取精确性。

设f_t(x，y)为第t次迭代计算得到的像元点(x，y)的参与计算图像数据，f_t-1(x，y)为第i-1次迭代计算得到的(x，y)的参与计算图像数据，更新图像数据的准则是避免人为引入误差或结构元素尺寸对提取端元纯度的影响，因此采用扩展膨胀操作实现图像数据更新，即$f_t(x,y)=(f_{t-1}\oplus b)(x,y).$

(6)计算迭代后的形态离心率指数MEI图像，并进行二值化，得到端元数据集；

为了自动将形态离心率指数MEI图像二值化，实现所有端元数据的提取，利用[m-3σ，m+3σ]进行MEI图像的二值化，得到所有端元的空间位置，其中，m为MEI图像均值，σ为MEI图像的标准差。

(7)利用由步骤(6)得到的端元数据集中任一元素作为参考光谱，采用光谱角匹配方法计算得到参考光谱与端元数据集中所有元素光谱角最大的元素，作为第一个端元，并向第一个端元的正交子空间投影，更新端元数据集；

为了克服现有基于扩展数学形态学方法无法准确区分不同类型端元的局限，本发明方法提出了基于正交子空间投影的方法，实现了不同类型端元的有效提取。

基于线性混合模型，利用正交子空间投影OSP概念进行混合像元分解时，仍可以继续将端元矩阵分解，则混合模型可以表示为：

r＝dα_d+Uα_U+n

其中，r为端元数据集，d为得到的端元光谱信号，α_d为端元光谱信号的混合百分比，U为剩余端元光谱信号，α_U为剩余端元光谱信号的混合百分比，n为噪声信号；通过MEI图像中提取的端元数据，利用光谱角匹配方法得到其中一个端元光谱信号d＝e₁，然后利用向其正交子空间投影得到来消除端元数据中的信号d＝e₁，其中$P_d^{\perp }=I-d{\left(d^{T}d\right)}^{-1}{d}^T.$将代入基于OSP的混合模型中得到$P_d^{\perp }r=P_d^{\perp }{\mathrm{U\alpha }}_U+P_d^{\perp }n,$可以看出，不仅消除了d＝e₁的影响而且有效的压制了噪声然后利用更新端元数据集。

(8)不断重复步骤(7)，直到满足误差要求，从而得到图像中所有端元。

经过m-1次投影之后得到d＝{e₁，e₂，...，e_m}，投影循环结束的条件由分解误差决定，从而得到了m个端元。分解误差为当前提取的端元与前一次提取的端元光谱变化值，本发明中的循环终止条件为：‖e_i-e_i-1‖₂≤0.01。

以上整个实现过程不需要人为干预或任何先验信息，因此，大大提高了端元提取方法的自动性。通过本发明方法，得到了西藏驱龙地区的五种主要的蚀变矿物端元光谱：贫铝绢云母，富铝绢云母，高岭石化，绿泥石化，褐铁矿化。