专利名称：扩展形态学与正交子空间投影结合的端元自动提取方法
技术领域：
本发明涉及扩展形态学与正交子空间投影结合的端元自动提取方法，属 于高光谱数据处理方法与应用技术领域，适用于高光谱数据端元提取的理论 方法以及在农作物精细分类、矿物蚀变信息提取等领域的应用技术研究。
背景技术：
高光谱数据处理与传统的遥感数据处理的一个本质的区别在于后者在 对地物分类时通常是将单一像元根据某种规则分到某一不同的种类中，从而 更加准确地反映遥感对象的状况。混合像元分解不仅是提高地物识别精度的 重要方法，更是定量探测地物成分的有效途径，因此，混合像元分解是高光 镨遥感数据应用的主要方向之一 ，该类技术的发展能够大大推动高光谱数据 的应用范围与水平，并不断扩展高光谱数据的应用深度和广度。
目前端元提取的方法主要包括纯像元指数(Pure Pixel Index, PPI)、
最小体积变换、迭代误差分析、顶点分量分析等，上述方法均只利用了高光 谱数据提供的光谱/特征空间信息，针对这一问题，国外提出了基于数学形 态学的方法，该方法同时利用了高光谱数据提供的光谱和空间信息，提高了 方法的可靠性。但是该类方法存在着一定程度的问题采用光谱角制图的方 法度量光谱之间的相似度，受噪声影响较大；不同类型的端元采用了大津法， 对于光谱相似的不同类型端元很难区分；并且算法的性能直接受结构元素形 状与尺寸的影响。
发明内容
本发明的目的是针对现有端元提取方法从数据光谱或空间或特征信息 的单一方面出发进行混合像元分解、不同类型端元难以区分等不足，提出一
种扩展形态学与正交子空间投影结合的端元自动提取方法。本发明的技术解决方案是利用扩展形态学操作与正交子空间投影理论 实现高光谱数据端元自动提取的方法，该方法主要是基于扩展数学形态学的 理论，利用扩展膨胀和腐蚀操作，并通过计算形态离心率指数进行高光谱数 据的端元数据集计算，然后利用光谱角匹配方法提取不同类型的端元，并通 过向端元正交子空间投影消除已经提取端元的影响，实现高光谱数据所有端 元的自动提取。形态偏心率指数是利用扩展膨胀和腐蚀的结果计算数据的正 交投影散度，从而将图像中端元提取出来，引入正交子空间投影方法实现不 同类型端元的有效区分，并且在本发明中为了避免由于形态学中结构元素的 大小对方法性能的限制和降低计算过程中引入的人为误差的影响，采用迭代 的方法，利用扩展膨胀操作不断更新参与计算的图像数据来实现。
本发明扩展形态学与正交子空间投影结合的端元自动提取方法，其步骤 如下
(1 )高光谱遥感数据读入；
(2)根据高光谱数据特点，确定结构元素的尺寸、初始迭代次数以及 最大迭代次数；
(3 )由扩展数学形态学膨胀和腐蚀操作计算每一个结构元素邻域内混 合度最低和最高的像元；
(4) 由步骤(3)得到的混合度最低和最高的像元结果计算得到形态离 心率指数值；
(5) 步骤(3) 、 (4)随迭代次数的增加不断重复，并利用步骤(3) 得到的混合度最低的结果不断更新图像数据，直到达到最大的迭代次数；
(6) 二值化形态离心率指数MEI图像，得到端元数据集；
(7) 利用由步骤(6)得到的端元数据集中任一元素作为参考光谱，采 用光谱角匹配方法计算得到参考光谱与端元数据集中所有元素光谱角最大 的元素，作为第一个端元，并向第一个端元的正交子空间投影，更新端元数 据集；(8)不断重复步骤(7)，直到满足误差要求，从而得到高光谱遥感图 像中所有的端元。
其中，步骤(2)中所述的结构元素尺寸为5x5，初始迭代次数,=1,最 大迭代次数^=5;
其中，步骤(3)所述的"由扩展数学形态学膨胀和腐蚀操作计算每一 个结构元素邻域内混合度最低和最高的像元"，其扩展数学形态学膨胀和腐 蚀操作定义如下
"0,力=(/ 6)0，力=arg一 Mxc {D(/0 /)， 6)}
力=(/ 6)0，力=arg— Af/" {/)(/0 + 5，^ + 0， 6)}
其中，"(x,力、e(x,力分别表示像元(;c,力在结构元素ZKx，力内的扩展膨胀、腐 蚀操作结果；/(x,力是图像像元(x,力的值，6(x,力是结构元素；arg—Ma;c， arg—M/" 分别表示使得"达到最大和最小的像元向量；D是为了确定像元混合度的多 维向量的排序关系，引入的一个多维向量的度量算子。该度量算子由结构元 素内各个像元累加距离计算得到，定义如下
其中，(W)eA表示0,/)是属于结构元素6(x,力内的点，力W是测量7V维向量的 逐点线性距离，采用正交投影散度(Orthogonal Projection Divergence, OPD)
计算该距离，考虑两个^维光谱信号^=^^,2,...,^]", ~=[^a2,...,v:t,
则iV维光谱信号&和、之间的正交投影散度OPD表示为
,一)=(《々/+《尸 /2
s)=^xw-A(《&)—、〖，"U'，并且/M是^N维的单位矩阵。
因此，累加距离D能够根据像元的混合程度大小排序结构元素中的向 量。通过以上的分析表明，扩展到高光谱图像的膨胀结果得到的是在结构元 素内纯度最高的像元，腐蚀结果得到的是在结构元素内混合度最高的像元。其中，步骤(4)中所述的"由步骤(3)得到的混合度最低和最高的像 元结果计算得到形态离心率指数值"，其含义说明如下通过正交投影散度 融合结构元素邻域内膨胀和腐蚀结果，计算得到形态离心率指数MEI,计算 公式如下
感/(x，力=O尸岸(x，力，e(x， y)]
其中，步骤(5)中所述的"步骤(3) 、 (4)随迭代次数的增加不断 重复，并利用步骤(3)得到的混合度最低的结果不断更新图像数据，直到 达到最大的迭代次数，，，其含义说明如下设/(x,力为第/次迭代计算得到的 像元点的图像数据，(,(x，力为第/-1次迭代计算得到的(x，力的图像数据， 更新图像数据的目的是避免人为引入误差或结构元素尺寸对提取端元纯度 的影响，因此釆用扩展膨胀操作实现图像数据更新，即/(x,力-(乂一e6)(x,力。
其中，步骤(6)中所述的"计算迭代后的形态离心率指数MEI图像， 并进行二值化，得到端元数据集"，其含义说明如下利用[w-3j，w + 3(j]进 行MEI图像的二值化，得到所有端元的空间位置，其中，w为MEI图像均 值，cx为MEI图像的标准差。
其中，步骤(7)中所述的"利用由步骤(6)得到的端元数据集中任一 元素作为参考光谱，采用光语角匹配方法计算得到参考光谱与端元数据集中
所有元素光谱角最大的元素，作为笫一个端元，并向第一个端元的正交子空 间投影，更新端元数据集"，其含义说明如下基于线性混合模型，利用正 交子空间才殳影(Orthogonal Subspace Projection, OSP)扭克念进4亍混合像元分解时， 仍可以继续将端元矩阵分解，则混合模型可以表示为
其中r为端元数据集，d为得到的端元光谱信号，^为端元光谱信号的混合 百分比，C/为剩余端元光谱信号，c^为剩余端元光镨信号的混合百分比，《
为噪声信号；通过MEI图像中提取的端元数据，利用光谱角匹配方法得到其
中一个端元光谱信号d=A ，然后利用向其正交子空间投影得到p/来消除端元数据中的信号^/ = ^其中尸/=/-^fc/)-7。将尸/代入基于OSP的混合模 型中得到P/r =尸/"^ +尸/"，可以看出，不仅消除了 d = q的影响而且有效的 压制了噪声P/n。然后利用/^更新端元数据集。
其中，步骤(8)中所述的"不断重复步骤(7),直到满足误差要求， 从而得到图像中所有端元"，其含义说明如下经过m-l次投影之后得到 d = {ei,e2,..., },投影循环结束的条件由分解误差决定，从而得到了 m个端元。
本发明与现有技术相比的优点在于本方法同时利用了扩展形态学与正 交子空间投影的理论，从数学机理上有效的克服了现有端元提取方法无法将 不同类型端元有效区分的问题，实现了无先验信息情况下高光谱数据端元自 动提取。它具有以下的优点(1 )通过扩展形态学膨胀和腐蚀操作实现图 像数据中纯度最高和混合度最高的像元提取，并利用正交投影散度融合扩展 膨胀和腐蚀结果得到形态离心率指数，实现了高光谱数据端元提取，有效地 抑制了背景等因素对端元提取结果的影响，提高了方法的稳定性；(2)采 用迭代的方法，利用扩展膨胀操作不断更新图像数据来克服由于形态学中结 构元素的尺寸等对方法性能的限制；(3)引入正交子空间投影和光谱角匹 配，实现不同类型端元的有效区分，提高了方法的精确性。
图1为本发明扩展形态学与正交子空间投影结合的端元自动提取方法的流程图。
具体实施例方式
为了更好的说明本发明涉及的扩展形态学与正交子空间投影结合的端 元自动提取方法，利用星载Hyperion高光谱成像仪数据，采用本发明方法 进行西藏驱龙地区蚀变矿物端元提取。本发明 一种扩展形态学与正交子空间 投影结合的端元自动提取方法，具体实现步骤如下
(1)高光谱遥感数据读入读入Hyperion高光i普数据，数据维数为
255x550x196 ;(2) 根据高光谱数据特点以及应用目的，确定结构元素的尺寸和初始 迭代次数结构元素尺寸为5x5，初始迭代次数/ = 1，最大迭代次数^^5;
(3) 由扩展数学形态学膨胀和腐蚀操作计算每一个结构元素邻域内混
合度最低和最高的像元；
为了实现更加可靠、稳定、精确的高光谱数据端元自动提取，本发明中 的方法同时利用了光谱和空间信息，由此引入了扩展数学形态学膨胀和腐蚀
操作，定义如下
d(x,力=(/ e 。(x，力=arg— {D(/(x - ^ _y —,), 6)}
"',')efl、
力=(/ 。(A力=arg_ M/" {￡ (/(x + s， y + f)，的}
(.、■,')",、■
其中，J(x,j;)、 e(x，力分别表示像元O，;;)在结构元素Z)O，力内的扩展膨胀、腐 蚀操作结果；/(x，力是图像像元(x,力的反射率值，6(;c,力是结构元素；arg_Mzx, arg_M>7分别表示使得"达到最大和最小的像元向量；D是为了确定像元混合 度的多维向量的排序关系，引入的一个高维向量的度量算子。该度量算子由 结构元素内各个像元累加距离计算得到，定义如下 D(/(xj),6) = Z2>、,(/"})，/(w)) , Cv)e A
其中，0力eA表示0力是属于结构元素6(;c,力内的点，^,是测量7V维向量的 逐点线性距离。为了有效地利用高光谱数据提供的光谱和空间信息，本发明 采用正交投影散度opd计算该距离，考虑两个196维光谱信号
s,[^，^,…,^J, 、=[^a2，...a.196]、则196维光谱信号s,和^之间的正交投 影散度opd表示为
《=/19函-A(")X， 并且/19&<196是196*196维的单位矩阵。
因此，累加距离"能够根据像元的混合程度大小排序结构元素中的向 量。通过以上的分析表明，扩展到高光镨图像的膨胀结果得到的是在结构元 素内纯度最高的像元，腐蚀结果得到的是在结构元素内混合度最高的像元。
10(4) 由步骤(3)得到的混合度最低和最高的像元结果计算得到形态离 心率指数值；
为了克服利用单一膨胀或腐蚀操作人为引入的噪声和背景影响，以及为 了定量化图像数据像元的纯度与混合度，采用了形态离心率指数。通过正交 投影散度融合结构元素邻域内膨胀和腐蚀结果，计算得到形态离心率指数 MEI，实现高光谱数据的端元提取。
ME/0，力=C^D[^0,力，e(x， y)]
(5) 步骤(3) 、 (4)随迭代次数的增加不断重复，并利用步骤(3) 得到的混合度最低的结果不断更新图像数据，直到达到最大的迭代次数；
为了提取图像数据中纯度最高的像元，利用能够得到结构元素内纯度最 高像元的扩展膨胀操作不断更新图像数据，从而有效地提高端元提取精确 性。
设乂(x，力为第/次迭代计算得到的像元点(;c,力的参与计算图像数据， ;,(x，力为第/-1次迭代计算得到的的参与计算图像数据，更新图像数据 的准则是避免人为引入误差或结构元素尺寸对提取端元纯度的影响，因此采 用扩展膨胀操作实现图像数据更新，即/0c,力^y;一④6)(;c,力。
(6) 计算迭代后的形态离心率指数MEI图像，并进行二值化，得到端 元数据集；
为了自动将形态离心率指数MEI图像二值化，实现所有端元数据的提
取，利用[m-3c7，W + 3f7]进行MEI图像的二值化，得到所有端元的空间位置，
其中，m为MEI图像均值，cr为MEI图像的标准差。
(7) 利用由步骤(6)得到的端元数据集中任一元素作为参考光谱，采 用光谱角匹配方法计算得到参考光谱与端元数据集中所有元素光谱角最大 的元素，作为第一个端元，并向第一个端元的正交子空间投影，更新端元数 据集；
为了克服现有基于扩展数学形态学方法无法准确区分不同类型端元的局
ii限，本发明方法提出了基于正交子空间投影的方法，实现了不同类型端元的 有效提取。
基于线性混合模型，利用正交子空间投影OSP概念进行混合像元分解时， 仍可以继续将端元矩阵分解，则混合模型可以表示为
其中，,为端元数据集，"为得到的端元光谱信号，《(,为端元光谱信号
的混合百分比，f/为剩余端元光语信号，c^为剩余端元光谱信号的混合百分 比，"为噪声信号；通过MEI图像中提取的端元数据，利用光谱角匹配方法 得到其中 一个端元光谱信号"=^,然后利用向其正交子空间投影得到来 消除端元数据中的信号dw"其中尸/一-^/V)-V7。将尸/代入基于0SP的 混合模型中得到尸/,=户/^(,+户/",可以看出，不仅消除了 d^的影响而且有 效的压制了噪声尸/"。然后利用尸>更新端元数据集。
(8)不断重复步骤(7)，直到满足误差要求，从而得到图像中所有端元。
经过m-l次投影之后得到^ = {21^2,..., }，投影循环结束的条件由分解误 差决定，从而得到了 m个端元。分解误差为当前提取的端元与前一次提取的 端元光谱变化值，本发明中的循环终止条件为-^,1^0.01。
以上整个实现过程不需要人为干预或任何先验信息，因此，大大提高了 端元提取方法的自动性。通过本发明方法，得到了西藏驱龙地区的五种主要 的蚀变矿物端元光谱贫铝绢云母，富铝绢云母，高岭石化，绿泥石化，褐 铁矿化。
权利要求
1、一种扩展形态学与正交子空间投影结合的端元自动提取方法，其特征在于它包含以下步骤(1)高光谱遥感数据读入；(2)根据高光谱数据特点，确定结构元素的尺寸、初始迭代次数以及最大迭代次数；(3)由扩展数学形态学膨胀和腐蚀操作计算每一个结构元素邻域内混合度最低和最高的像元；(4)由步骤(3)得到的混合度最低和最高的像元结果计算得到形态离心率指数值；(5)步骤(3)、(4)随迭代次数的增加不断重复，并利用步骤(3)得到的混合度最低的结果不断更新图像数据，直到达到最大的迭代次数；(6)二值化形态离心率指数MEI图像，得到端元数据集；(7)利用由步骤(6)得到的端元数据集中任一元素作为参考光谱，采用光谱角匹配方法计算得到参考光谱与端元数据集中所有元素光谱角最大的元素，作为第一个端元，并向第一个端元的正交子空间投影，更新端元数据集；(8)不断重复步骤(7)，直到满足误差要求，从而得到高光谱遥感图像中所有的端元。
2、 根据权利要求1所述的扩展形态学与正交子空间投影结合的端元自动提取方法，其特征在于步骤(3)中所述的由扩展数学形态学膨胀和腐 蚀操作计算每一个结构元素邻域内混合度最低和最高的像元，计算公式如下力=(/ W(A力=arg— Max {D(/0 — _y — 0, e(x,力=(/ 6)0，力=arg一 Mz'" { "(/O + ， _y + 0,其中，雖，力、e(x,力分别表示像元0，力在结构元素6(;c，力内的扩展膨胀、腐 蚀操作结果；/0,力是图^f象^(象元(;c，力的值，^c,力是结构元素；arg—Max, arg—M'" 分别表示使得。达到最大和最小的像元向量；D是一个高维向量的度量算子，计算公式如下<formula>formula see original document page 3</formula>其中，(W)eA表示GM)是属于结构元素^c,力内的点，采用正交投影散度计 算两个7V维光谱信号& =， &,... Av]7'和~ =， & ，..., > f之间的必r距离，正交 投影散度OPD表示为其中，《=/歸1(")乂， "^',并且/歸是ITN维的单位矩阵。
3、 根据权利要求1所述的扩展形态学与正交子空间投影结合的端元自 动提取方法，其特征在于其步骤(4)中所述的计算形态离心率指数MEI 值计算公式如下細/(x,力=OPD,，力，e(x, y)] 其中，OPD表示正交投影散度，d(x，力表示像元(;c,力在结构元素6(;c，W内的 扩展膨胀操作结果，e(;c,力表示像元(x,少)在结构元素60c,力内的扩展腐蚀操作 结果。
4、 根据权利要求1所述的扩展形态学与正交子空间投影结合的端元自 动提取方法，其特征在于步骤(5)中所述的更新图像数据采用扩展膨胀 操作实现其中，X0c,W为第,次迭代计算得到的像元点0c，力的图像数据，；,(x，力为第 f-1次迭代计算得到的(x，力的图像数据，6(x,力为结构元素。
5、 根据权利要求1所述的扩展形态学与正交子空间投影结合的端元自 动提取方法，其特征在于步骤(6)中所述的二值化形态离心率指数MEI图像方法如下利用[w-3c7，m + 3cr]进行MEI图像的二值化，得到所有端元的空间位置，其中，m为MEI图像均值，cr为MEI图像的标准差。
6、根据权利要求1所述的扩展形态学与正交子空间投影结合的端元自动提取方法，其特征在于步骤(7)中所述的正交子空间投影计算如下利用正交子空间投影进行混合像元分解模型表示为其中，r为端元数据集，"为得到的端元光谱信号，c^为端元光谱信号的混合百分比，f/为剩余端元光谱信号， 为剩余端元光谱信号的混合百分比，"为噪声信号；然后利用向其正交子空间投影得到尸/ = /-^^)-7，将p/代入基于正交子空间投影的混合模型中得到p》=尸/t/ + p/"实现端元数据集的更新。
全文摘要
一种扩展形态学与正交子空间投影结合的端元自动提取方法，步骤如下(1)高光谱数据读入；(2)确定结构元素尺寸、初始与最大迭代次数；(3)由扩展膨胀、腐蚀操作计算结构元素邻域内混合度最低和最高的像元；(4)由步骤(3)得到的结果计算得到形态离心率指数值；(5)步骤(3)、(4)随迭代次数的增加不断重复，并利用步骤(3)扩展膨胀操作结果更新图像数据，直到达到最大迭代次数；(6)二值化MEI图像，得到端元数据集；(7)采用光谱角匹配方法计算得到第一个端元，并向得到端元的正交子空间投影更新端元数据集；(8)重复步骤(7)，直到满足误差要求。该方法是一种稳定性强、可靠性高、精确度高的高光谱端元自动提取方法。
文档编号G01J3/00GK101504315SQ20091007844
公开日2009年8月12日 申请日期2009年2月23日 优先权日2009年2月23日
发明者州 徐, 娜 李, 辉 蔡, 贾国瑞, 赵慧洁 申请人:北京航空航天大学