CN104952050A - 基于区域分割的高光谱图像自适应解混方法 - Google Patents

Abstract

本发明公开了一种基于区域分割的高光谱图像自适应解混方法，考虑线性混合和双线性混合共存的情况，其实现有：输入高光谱图像；用基于最小错误高光谱信号识别法估计端元数目；用顶点成分分析算法提取端元矩阵；用K均值聚类方法对高光谱数据聚类，将图像分割为匀质区域和细节区域；匀质区域采用线性模型，用稀疏约束的非负矩阵分解方法解混，细节区域采用广义双线性模型，用稀疏约束的半非负矩阵分解方法解混。本发明结合高光谱数据光谱和丰度的特点，获得高光谱图像更精确的表示，提高解混的准确率。在丰度上添加稀疏约束条件，克服了半非负矩阵分解算法易陷入局部最小值的缺点，得到更加准确的丰度，用于高光谱图像的地物识别。

Description

基于区域分割的高光谱图像自适应解混方法

技术领域

本发明属于图像处理技术领域，主要涉及无监督目标识别方法，具体是一种基于区域分割的高光谱图像自适应解混方法，可应用于高光谱图像的地物识别。

背景技术

高光谱遥感是利用测谱学原理，即在电磁波谱的紫外、可见光、近红外和中红外区域，利用成像光谱仪获取的许多非常窄且光谱连续的图像数据。成像光谱仪以像元为单位来获取的地面反射或发射光谱信号。图像中每个像元对应的地物空间区域，往往包含有着不同的光谱特征的不同物质。若该像元中仅包含一种物质或者该物质所占的比例很高，则称为纯像元，也叫作“端元”；若该像元包含不止一种物质，该像素点被称作混合像元。混合像元所对应的光谱曲线由几种不同物质的光谱曲线混合而成。

混合像元分解方法按照所采用的分解模型，大致可以分为基于线性光谱混合模型的分解方法和基于非线性光谱混合模型的分解方法。

光谱解混提高了目标识别的精度，是高光谱图像处理中一项非常重要的任务，线性模型(LMM)假设观察到的光谱是一些端元的线性组合，它是一个简化的光谱模型，只考虑一阶散射光子，而忽略多个光子间的相互作用。尽管基于LMM的混合端元分解方法能得到有物理意义的结果，光谱混合模型中的非线性分量已在很多工作中被指出。基于线性混合模型的光谱解混理论和方法已得到了广泛的研究与应用，但由于受实际地物间复杂关系以及大气散射的影响，光谱混合都是非线性的，这就使得应用传统的基于线性光谱混合模型的解混结果难以满足高光谱图像解混的精度要求。近年来，非线性解混是基于非线性模型的，在遥感图像处理中非常受重视。多层场景中，在不同的层中有多重散射，通常发生在植被和土壤之间。根据地物散射的特点，双线性模型(BMM)考虑了光子在不同材料之间的二阶散射。双线性比较适合于高光谱图像中灌木丛、植被等会产生折射现象的地类目标。双线性模型是非线性解混中的常用模型。

非负矩阵分解(NMF)将一个非负矩阵分解为两个非负矩阵相乘，是目前处理线性混合模型的非常有效的方法。半非负矩阵分解(semi-NMF)将一个无约束矩阵分解为一个无约束矩阵和一个非负矩阵相乘,可用来优化基于广义双线性模型(GBM)的解混。然而高光谱图像往往不是单纯只存在线性光谱混合或者只存在双线性光谱混合，大多数高光谱数据需要考虑线性混合和双线性混合模型共存的情况，目前还没有有效的方法来处理这样的图像。而且目前现有的基于广义双线性模型GBM的解混方法没有考虑丰度矩阵的稀疏性，容易限于局部最优解。

发明内容

为了克服现有的高光谱图像解混方法不适用于处理线性混合和双线性混合模型共存的高光谱数据，本发明提供一种考虑线性混合和双线性混合共存的基于区域分割的高光谱图像自适应解混方法。

本发明的基于区域分割的高光谱图像自适应解混方法，其实现包括有如下步骤：

(1)输入高光谱图像，高光谱图像用数据Y∈R^L×N表示，其中L表示高光谱数据的波段数，N表示高光谱数据样本总个数，R表示实数域；

(2)使用基于最小错误高光谱信号识别法估计高光谱数据Y∈R^L×N的信号子空间，得到信号子空间维数K，即高光谱数据的端元数目；

(3)用顶点成分分析算法提取高光谱图像的端元矩阵A∈R^L×K，其每行向量表示高光谱图像第i个端元的光谱曲线；

(4)用K均值聚类方法(kmeans)对高光谱数据Y∈R^L×N聚类，聚为K+2类，聚类类标分别为1,2,3...,K+2，根据类标，用2*2的窗口扫描全图，若窗口中至少有一个类标不同，则把该2*2窗口对应的像素点归为高光谱图像的细节区域数据其余像素点归为高光谱图像的匀质区域数据其中N＝N₁+N₂，N₁表示高光谱图像的匀质区域数据样本个数，N₂表示高光谱图像的细节区域数据样本个数；本发明在该步骤中把高光谱图像分割为细节区域和匀质区域分别表示，在之后的处理步骤中分别进行处理；

(5)将高光谱图像匀质区域数据Y₁和高光谱图像端元矩阵A作为线性光谱混合模型的输入信号，用L_1/2约束的非负矩阵分解(L_1/2-NMF)方法得到匀质区域数据的系数矩阵即一阶丰度矩阵，其中每一列向量第n个像素的丰度向量；本发明在高光谱图像匀质区域采用线性模型，用稀疏约束的非负矩阵分解(L_1/2-NMF)方法解混。

(6)将高光谱图像细节区域数据Y₂和高光谱图像端元矩阵A作为双线性光谱混合模型的输入信号，用稀疏约束的半非负矩阵分解方法得到细节区域数据Y₂对应的一阶丰度矩阵其中每一列向量表示第m个像素的丰度向量，和细节区域数据Y₂对应的二阶丰度矩阵其中每一列向量表示第m个像素的双线性丰度向量；本发明在高光谱图像细节区域采用广义双线性模型(GBM)，用稀疏约束的半非负矩阵分解(L_1/2-Semi-NMF)方法解混。

(7)将高光谱图像匀质区域Y₁的一阶丰度矩阵和高光谱图像细节区域Y₂的一阶丰度矩阵合并为得到为整个高光谱数据的丰度矩阵，完成高光谱图像解混。

本发明中高光谱图像考虑了线性混合和双线性混合模型共存的情况，通过K均值聚类方法把高光谱图像分割为细节区域和匀质区域分别表示。其中匀质区域采用线性模型，用稀疏约束的非负矩阵分解(L_1/2-NMF)方法解混；细节区域采用广义双线性模型(GBM)，用稀疏约束的半非负矩阵分解(L_1/2-Semi-NMF)方法解混，在GBM模型中考虑了丰度的稀疏信息，使得解混更加稳定，不易陷于局部最优解。

本发明的实现还在于(5)所述的用L_1/2约束的非负矩阵分解方法得到匀质区域数据Y₁的一阶丰度矩阵X₁，包括有如下步骤：

(5a)根据高光谱成像理论，在高光谱图像匀质区域数据的丰度矩阵X₁中加入L_1/2范数，得到稀疏约束表达式作为丰度矩阵X₁的稀疏约束项，其中x_1n(k)是高光谱图像匀质区域数据Y₁中在第n个像元的第k个端元的丰度；

(5b)将步骤(5a)得到的稀疏约束项添加到以欧氏距离为基础的非负矩阵分解算法的目标函数中，构成新的目标函数：

$min\frac{1}{2}\left|\right|Y_1-{\mathrm{AX}}_1|{|}_2^2+\mathrm{\λ}|\left|X_1\right|{|}_{1/2}$

条件：X₁≥0，1^TX₁＝1^T

其中λ是稀疏约束正则化参数，X₁≥0，1^TX₁＝1^T是高光谱图像匀质区域丰度矩阵的“非负”和“和为一”约束；

(5c)对步骤(5b)得到的目标函数用迭代乘法进行优化，得到匀质区域数据的端元矩阵A的更新公式和匀质区域数据的一阶丰度矩阵X₁的更新公式其中(·)^T表示矩阵的转置，·*和·/分别表示按元素的乘法和除法；

(5d)重复执行步骤(5c)得到的高光谱图像匀质区域数据的端元矩阵A的更新公式和匀质区域数据的一阶丰度矩阵X₁的更新公式，直到迭代次数达到设定次数c。

目前正则化方法通常用来约束丰度的稀疏性，因为大多数像素只存在少数的端元的混合。稀疏约束可以通过正则化方法实现，如L₀正则和L₁正则。前者可以指定的非零元素的数目，但难以求解，而后者只能控制稀疏正则化参数，但不是精确控制非零数。本发明采用L_1/2约束的非负矩阵分解方法得到匀质区域数据Y₁的一阶丰度矩阵X₁，提供了比使用L₁更稀疏，更准确的结果。

本发明的实现还在于步骤(6)用L_1/2-Semi-NMF方法得到细节区域数据Y₂对应的一阶丰度矩阵X₂和二阶丰度矩阵E，按照如下步骤进行：

(6a)在高光谱图像细节区域Y₂采用双线性模型表示如下

Y₂＝AX₂+BE+M

其中，表示双线性区域数据对应的一阶丰度矩阵，其中每一列向量表示第n个像素的丰度向量，是双线性端元矩阵,是双线性区域数据对应的二阶丰度矩阵，其中每一列向量表示第n个像素的双线性丰度向量，表示噪声矩阵；

(6b)在高光谱图像匀质区域数据的丰度矩阵X₂中加入L_1/2范数，得到稀疏约束表达式作为丰度矩阵X₂的稀疏约束项，其中x_2n(k)是高光谱图像细节区域数据Y₂中在第n个像元的第k个端元的丰度；

(6c)将步骤(6b)得到的稀疏约束项添加到步骤(6a)所述双线性模型上，得到双线性模型的最小化目标函数

$\min \left|\right|Y_2-{\mathrm{AX}}_2-BE|{|}_F^2+\mathrm{\λ}|\left|X_2\right|{|}_{1/2}$

条件：

其中(X^*)_(i,j)n＝(X₂)_in(X₂)_jn(n∈{1,2,...,N₂}，||·||_F表示F范数；

(6d)引入Y₂₁＝Y₂-AX₂，Y₂₂＝Y₂-BE，可将步骤(6a)中的双线性模型写成如下形式Y₂₁＝AX₂+M和Y₂₂＝ΒΕ+M，然后分别关于A和B迭代的优化原优化问题可以通过迭代优化算法转化为优化两个局部优化问题；

(6e)因为物理约束，A，X₂，B和E都是非负的。用迭代乘法优化步骤(6c)得到的局部优化问题，得到细节区域数据Y₂对应的一阶丰度矩阵X₂的更新公式和细节区域数据Y₂对应的二阶丰度矩阵E的更新公式 $E^T.*\sqrt{{\left(Y_{22}^{T}B\right)}^{+}+E^T{\left(B^{T}B\right)}^{-}/{\left(Y_{22}^{T}B\right)}^{-}+E^T{\left(B^{T}B\right)}^{+}};$

(6f)重复执行步骤(6e)得到的细节区域数据Y₂对应的一阶丰度矩阵X₂的更新公式和细节区域数据Y₂对应的二阶丰度矩阵E的更新公式，直到迭代次数达到设定次数c。

大多数现有的非线性解混算法没有考虑数据的稀疏的信息，而高光谱数据的高度相关性导致了丰度的稀疏性。直接用半非负矩阵分解方法来解混GBM的线性部分，这导致收敛速度缓慢，而且容易陷入局部最优解。为了克服这些缺点，本发明通过添加稀疏约束来拓展GBM，以提高解混的稳定性和结果的正确性。

本发明与现有的技术相比具有以下优点：

1、本发明由于采用了基于区域分割的高光谱图像自适应解混方法，考虑了线性混合和双线性混合模型共存的情况，通过K均值聚类方法把高光谱图像分割为细节区域和匀质区域分别表示，其中匀质区域采用线性模型，用稀疏约束的非负矩阵分解方法解混；细节区域采用广义双线性模型，用稀疏约束的半非负矩阵分解方法解混，很好的保持了双线性丰度的边缘信息，从而获得高光谱图像更精确的表示，提高解混的准确率。

2、本发明由于在细节区域数据的双线性混合模型中考虑了丰度矩阵的稀疏性，收敛速度加快且不易限于局部最优解，提高了解混的稳定性和结果的正确性。

对比实验表明，本发明有效提高了高光谱遥感图像的解混准确率。

附图说明

图1是本发明的流程示意图；

图2是本发明仿真采用的真实图像第30波段的rgb图；

图3是本发明根据GBM模型的莫菲特场用本发明方法估计的丰度图；

图4是本发明根据GBM模型的莫菲特场用本发明方法估计的二阶丰度图。

具体实施方式

无监督光谱解混是高光谱遥感信号处理中一个重要的技术领域，要解决的技术问题是识别所拍摄的场景中存在的物质以及各个物质所占的组分。高光谱遥感技术采集图像过程中利用成像光谱仪纳米级的光谱分辨率克服了传统的遥感成像波段范围窄、光谱分辨率低等局限性，以较窄的波段区间、较多的波段数量提供遥感信息，使得解混能从光谱空间中对地物予以细分和鉴别，在海洋、植被、地质、大气、环境、军事、医学等方面都有广泛的应用和更广阔的应用前景。

实施例1

本发明是一种基于区域分割的高光谱图像自适应解混方法，参照图1，本发明的具体实施步骤包括：

(1)输入高光谱图像，高光谱图像用数据Y∈R^L×N表示，其中L表示高光谱数据的波段数，N表示高光谱数据样本总个数，R表示实数域；本例中输入高光谱图像第30波段的rgb图如图2所示，图中包含了植被、水和土壤三种物质。高光谱解混就是要得到这三种不同的物质分别占的百分比，即丰度图。

(2)使用基于最小错误高光谱信号识别法估计高光谱数据Y∈R^L×N的信号子空间，得到信号子空间维数K，即高光谱数据的端元数目，本例中端元数目是3。

(3)用顶点成分分析算法提取高光谱图像的端元矩阵A∈R^L×K，其每行向量表示高光谱图像第i个端元的光谱曲线，分别为植被、水和土壤这3个端元的光谱曲线。

以上是对高光谱图像的预处理。

(4)用K均值聚类方法(kmeans)对高光谱数据Y∈R^L×N聚类，聚为K+2类，类标分别为1,2,3...,K+2，根据类标，用2*2的窗口扫描全图，若窗口中类标全部相同，则把该2*2窗口对应的像素点归为高光谱图像的匀质区域数据其余像素点归为高光谱图像细节区域数据其中N＝N₁+N₂。本例中举例说明，若2*2窗口对应的像素点的类标分别为1，1，1，2，或1，1，2，3，则将这2*2窗口对应的4个像素点归为高光谱图像的细节区域；若2*2窗口对应的像素点的类标分别为1，1，1，1，或2，2，2，2，则将这2*2窗口对应的4个像素点归为高光谱图像的匀质区域。在本领域中，匀质区域也称作平滑区域。

本发明在该步骤中把高光谱图像分割为细节区域和匀质区域分别表示，并且在之后的处理步骤中对不同的技术区域分别进行不同的处理，以获得更加精确的解混效果。

(5)将高光谱图像匀质区域数据Y₁和高光谱图像端元矩阵A作为线性光谱混合模型的输入信号，用L_1/2约束的非负矩阵分解(L_1/2-NMF)方法得到匀质区域数据的系数矩阵即一阶丰度矩阵，其中每一列向量表示第n个像素的丰度向量。

(6)将高光谱图像细节区域数据Y₂和高光谱图像端元矩阵A作为双线性光谱混合模型的输入信号，用稀疏约束的半非负矩阵分解(L_1/2-Semi-NMF)方法得到细节区域数据Y₂对应的一阶丰度矩阵其中每一列向量表示第n个像素的丰度向量，和细节区域数据Y₂对应的二阶丰度矩阵其中每一列向量表示第n个像素的双线性丰度向量；

(7)将高光谱图像匀质区域Y₁的一阶丰度矩阵和高光谱图像细节区域Y₂的一阶丰度矩阵合并为得到为整个高光谱数据的丰度矩阵，完成高光谱图像解混，得到植被、水和土壤的一阶丰度图。植被、水和土壤之间的相互二阶丰度图即为细节区域的二阶丰度矩阵

本发明考虑了线性混合和双线性混合模型共存的情况，通过K均值聚类方法把高光谱图像分割为细节区域和匀质区域分别表示，其中匀质区域采用线性模型，用稀疏约束的非负矩阵分解方法解混；细节区域采用广义双线性模型，用稀疏约束的半非负矩阵分解方法解混，很好的保持了双线性丰度的边缘信息，从而获得高光谱图像更精确的表示，提高解混的准确率。

实施例2

基于区域分割的高光谱图像自适应解混方法同实施例1，其中步骤(5)所述的用L_1/2约束的非负矩阵分解方法得到匀质区域数据Y₁的一阶丰度矩阵X₁，包括有如下步骤：

(5a)根据高光谱成像理论，在高光谱图像匀质区域数据的丰度矩阵X₁中加入L_1/2范数，得到稀疏约束表达式作为丰度矩阵X₁的稀疏约束项，其中x_1n(k)是高光谱图像匀质区域数据Y₁中在第n个像元的第k个端元的丰度。

(5b)将步骤(5a)得到的稀疏约束项添加到以欧氏距离为基础的非负矩阵分解算法的目标函数中，构成新的目标函数：

$min\frac{1}{2}\left|\right|Y_1-{\mathrm{AX}}_1|{|}_2^2+\mathrm{\λ}|\left|X_1\right|{|}_{1/2}$

条件：X₁≥0，1^TX₁＝1^T

其中λ是稀疏约束正则化参数，X₁≥0，1^TX₁＝1^T是高光谱图像匀质区域丰度矩阵的“非负”和“和为一”约束；本例中正则化参数λ的取值范围一般为(0.001,0.5)。

(5c)对步骤(5b)得到的目标函数用迭代乘法进行优化，得到匀质区域数据的端元矩阵A的更新公式和匀质区域数据的一阶丰度矩阵X₁的更新公式其中(.)^T表示矩阵的转置，·*和·/分别表示按元素的乘法和除法；高光谱图像端元丰度的和为一约束可以采用在高光谱图像匀质区域数据矩阵Y₁和高光谱图像端元矩阵A上加一行常量的方法，定义如下：

$\begin{array}{cc}{Y}_{1f}=\left[\begin{array}{c}{Y}_1\\ \mathrm{\δ}{1}_N^T\end{array}\right]& A_f=\left[\begin{array}{c}A\\ \mathrm{\δ}{1}_K^T\end{array}\right]\end{array}$

其中δ控制和为一约束对最终端元丰度的影响。δ越大，X₁的列的总和越接近一。所以用Y_1f和A_f分别替代Y₁和A，作为X₁的考虑和为一约束的更新公式 $X_1.*A_f^T{Y}_{1f}/(A_f^T{A}_f{X}_1+\frac{\mathrm{\λ}}{2}{X_1}^{-\frac{1}{2}});$ 本例中参数δ的取值为25。

(5d)重复执行步骤(5c)得到的高光谱图像匀质区域数据的端元矩阵A的更新公式和匀质区域数据的一阶丰度矩阵X₁的更新公式，直到迭代次数达到设定次数c，本例中迭代次数范围是(100,1000)，因为如果迭代次数太大，计算时间长，迭代次数太小，精度不高，所以一般取400。

本发明采用L_1/2约束的非负矩阵分解方法得到匀质区域数据Y₁的一阶丰度矩阵X₁，提供了比使用L₁更稀疏，更准确的结果。

实施例3

基于区域分割的高光谱图像自适应解混方法同实施例1-2,其中步骤(6)用L_1/2-Semi_NMF方法得到细节区域数据Y₂对应的一阶丰度矩阵X₂和二阶丰度矩阵E，按照如下步骤进行：

(6a)在高光谱图像细节区域Y₂采用双线性模型表示如下

Y₂＝AX₂+BE+M

其中，表示双线性区域数据对应的一阶丰度矩阵，其中每一列向量表示第n个像素的丰度向量，是双线性端元矩阵,是双线性区域数据对应的二阶丰度矩阵，其中每一列向量表示第n个像素的双线性丰度向量，表示噪声矩阵；

(6b)在高光谱图像匀质区域数据的丰度矩阵X₂中加入L_1/2范数，得到稀疏约束表达式作为丰度矩阵X₂的稀疏约束项，其中x_2n(k)是高光谱图像细节区域数据Y₂中在第n个像元的第k个端元的丰度；

(6c)将步骤(6b)得到的稀疏约束项添加到步骤(6a)所述双线性模型上，得到双线性模型的最小化目标函数

$\min \left|\right|Y_2-{\mathrm{AX}}_2-BE|{|}_F^2+\mathrm{\λ}|\left|X_2\right|{|}_{1/2}$

条件： $X_2\≥0,\underset{i=1}{\overset{K}{\Σ}}{\left(X_2\right)}_{il}=1,0\≤E\≤X^{*}$

其中(X^*)_(i,j)n＝(X₂)_in(X₂)_jn(n∈{1,2,...,N₂}，||·||_F表示F范数；

(6d)引入Y₂₁＝Y₂-AX₂，Y₂₂＝Y₂-BE，可将步骤(6a)中的双线性模型写成如下形式Y₂₁＝AX₂+M和Y₂₂＝ΒΕ+M，然后分别关于A和B迭代的优化原优化问题可以通过迭代优化算法转化为优化两个局部优化问题；

(6e)用迭代乘法优化步骤(6c)得到的局部优化问题，得到细节区域数据Y₂对应的一阶丰度矩阵X₂的更新公式)和细节区域数据Y₂对应的二阶丰度矩阵E的更新公式 $E^T.*\sqrt{{\left(Y_{22}^{T}B\right)}^{+}+E^T{\left(B^{T}B\right)}^{-}/{\left(Y_{22}^{T}B\right)}^{-}+E^T{\left(B^{T}B\right)}^{+}};$

(6f)重复执行步骤(6e)得到的细节区域数据Y₂对应的一阶丰度矩阵X₂的更新公式和细节区域数据Y₂对应的二阶丰度矩阵E的更新公式，直到迭代次数达到设定次数c，本例中迭代次数范围是(100,1000)，因为如果迭代次数太大，计算时间长，迭代次数太小，精度不高，所以一般取400。

本发明通过添加稀疏约束来拓展GBM，收敛速度加快且不易限于局部最优解，提高了解混的稳定性和结果的正确性。

实施例4

基于区域分割的高光谱图像自适应解混方法，同实施例1-3，包括有如下步骤：

(1)输入高光谱图像，高光谱图像用数据Y∈R^L×N表示，其中L表示高光谱数据的波段数，N表示高光谱数据样本总个数，R表示实数域；

(2)使用基于最小错误高光谱信号识别法估计高光谱数据Y∈R^L×N的信号子空间，得到信号子空间维数K，即高光谱数据的端元数目；

(3)用顶点成分分析算法提取高光谱图像的端元矩阵A∈R^L×K，其每行向量表示高光谱图像第i个端元的光谱曲线；

(4)用K均值聚类方法(kmeans)对高光谱数据Y∈R^L×N聚类，聚为K+2类，类标分别为1,2,3...,K+2，根据类标，用2*2的窗口扫描全图，若窗口中至少有一个类标不同，则把该2*2窗口对应的像素点归为高光谱图像的细节区域数据其余像素点归为高光谱图像匀质区域数据其中N＝N₁+N₂，N₁表示高光谱图像的匀质区域数据样本个数，N₂表示高光谱图像的细节区域数据样本个数；

(5)将高光谱图像匀质区域数据Y₁和高光谱图像端元矩阵A作为线性光谱混合模型的输入信号，用L_1/2约束的非负矩阵分解(L_1/2-NMF)方法得到匀质区域数据的系数矩阵即一阶丰度矩阵，其中每一列向量表示第n个像素的丰度向量；包括有如下步骤：

(5a)在高光谱图像匀质区域数据的丰度矩阵X₁中加入L_1/2范数，得到稀疏约束表达式作为丰度矩阵X₁的稀疏约束项，其中x_1n(k)是高光谱图像匀质区域数据Y₁中在第n个像元的第k个端元的丰度；

(5b)将步骤(5a)得到的稀疏约束项添加到以欧氏距离为基础的非负矩阵分解算法的目标函数中，构成新的目标函数：

条件：X₁≥0，1^TX₁＝1^T

其中λ是稀疏约束正则化参数，X₁≥0，1^TX₁＝1^T是高光谱图像匀质区域丰度矩阵的“非负”和“和为一”约束；

(5c)对步骤(5b)得到的目标函数用迭代乘法进行优化，得到匀质区域数据的端元矩阵A的更新公式和匀质区域数据的一阶丰度矩阵X₁的更新公式其中(.)^T表示矩阵的转置，·*和·/分别表示按元素的乘法和除法；

(5d)重复执行步骤(5c)得到的高光谱图像匀质区域数据的端元矩阵A的更新公式和匀质区域数据的一阶丰度矩阵X₁的更新公式，直到迭代次数达到设定次数c。

(6)将高光谱图像细节区域数据Y₂和高光谱图像端元矩阵A作为双线性光谱混合模型的输入信号，用稀疏约束的半非负矩阵分解(L_1/2-Semi-NMF)方法得到细节区域数据Y₂对应的一阶丰度矩阵其中每一列向量表示第n个像素的丰度向量，和细节区域数据Y₂对应的二阶丰度矩阵其中每一列向量表示第n个像素的双线性丰度向量；

(6a)在高光谱图像细节区域Y₂采用双线性模型表示如下

Y₂＝AX₂+BE+M

其中，表示双线性区域数据对应的一阶丰度矩阵，其中每一列向量表示第n个像素的丰度向量，是双线性端元矩阵,是双线性区域数据对应的二阶丰度矩阵，其中每一列向量表示第n个像素的双线性丰度向量，表示噪声矩阵；

(6b)在高光谱图像匀质区域数据的丰度矩阵X₂中加入L_1/2范数，得到稀疏约束表达式作为丰度矩阵X₂的稀疏约束项，其中x_2n(k)是高光谱图像细节区域数据Y₂中在第n个像元的第k个端元的丰度；

(6c)将步骤(6b)得到的稀疏约束项添加到步骤(6a)所述双线性模型上，得到双线性模型的最小化目标函数

$\min \left|\right|Y_2-{\mathrm{AX}}_2-BE|{|}_F^2+\mathrm{\λ}|\left|X_2\right|{|}_{1/2}$

条件： $X_2\≥0,\underset{i=1}{\overset{K}{\Σ}}{\left(X_2\right)}_{il}=1,0\≤E\≤X^{*}$

其中(X^*)_(i,j)n＝(X₂)_in(X₂)_jn(n∈{1,2,...,N₂}，||·||_F表示F范数；

(6d)引入Y₂₁＝Y₂-AX₂，Y₂₂＝Y₂-BE，可将步骤(6a)中的双线性模型写成如下形式Y₂₁＝AX₂+M和Y₂₂＝ΒΕ+M，然后分别关于A和B迭代的优化

(6e)用迭代乘法优化步骤(6c)得到的局部优化问题，得到细节区域数据Y₂对应的一阶丰度矩阵X₂的更新公式)和细节区域数据Y₂对应的二阶丰度矩阵E的更新公式 $E^T.*\sqrt{{\left(Y_{22}^{T}B\right)}^{+}+E^T{\left(B^{T}B\right)}^{-}/{\left(Y_{22}^{T}B\right)}^{-}+E^T{\left(B^{T}B\right)}^{+}};$

(6f)重复执行步骤(6e)得到的细节区域数据Y₂对应的一阶丰度矩阵X₂的更新公式和细节区域数据Y₂对应的二阶丰度矩阵E的更新公式，直到迭代次数达到设定次数c。

(7)将高光谱图像匀质区域Y₁的一阶丰度矩阵和高光谱图像细节区域Y₂的一阶丰度矩阵合并为得到为整个高光谱数据的丰度矩阵，完成高光谱图像解混。

本发明通过K均值聚类方法把高光谱图像分割为细节区域和匀质区域分别表示。其中匀质区域用稀疏约束的非负矩阵分解方法解混；细节区域用稀疏约束的半非负矩阵分解方法解混，在GBM模型中考虑了丰度的稀疏信息，使得解混更加稳定，不易陷于局部最优解。

实施例5

基于区域分割的高光谱图像自适应解混方法，同实施例1-4，本发明的效果通过以下仿真实验进一步说明：

线性模拟数据由美国地质调查局(USGS)光谱库(共498种光谱信号)随机选择的3种光谱。为了产生合成数据，生成类似地面实况的丰度矩阵。步骤如下：先把z²×z²大的图像分成z×z个区域，每个区域用相同类型的端元初始化，即随机选择其中的一个端元填充。然后对每个像素用(z+1)×(z+1)的低通滤波器来生成混合数据，并且使产生的像素变化匀质。然后加入双线性丰度生成基于广义双线性模型GBM图像。

仿真实验在CPU为Intel Core(TM)2Duo、主频2.33GHz，内存为2G的WINDOWS 7系统上用MATLAB R2012b软件进行。

仿真内容及分析：

使用本发明与现有技术中三种方法对合成数据进行解混，现有三种方法分别是：全约束最小二乘(FCLS)，L_1/2约束的非负矩阵分解(L_1/2-NMF)，半非负矩阵分解(Semi-NMF)。本发明设置参数z为5，步骤(5b)和(6c)中参数λ设置为0.01，步骤(5d)和(6f)中迭代次数c设置为400。表中RE为重建误差，RMSE为均方根误差。

表1

	FCLS	L1/2-NMF	Semi-NMF	本发明
					RE	0.0031	0.0594	0.0805	0.0625
RMSE	0.021037	0.04795	0.020344	0.0038636

表1显示了本发明与现有技术中三种方法对三端元的合成数据进行解混性能的比较：解混的RE值和丰度的RMSE值。解混的RE值和丰度的RMSE值越小，解混效果越好。由重建误差表明FCLS解混效果较好，本发明的方法优于Semi-NMF方法的解混效果。由均方根误差表明本发明的RMSE值最小，且比其他三种方法小一个数量级。所以综合对比解混的RE值和丰度的RMSE值，明显看出本发明的解混效果较好。

实施例6

基于区域分割的高光谱图像自适应解混方法，同实施例1-4，本发明的效果通过以下仿真实验进一步说明：

真实的高光谱数据摄于1997年，是莫菲特场，位于加利福尼亚州的旧金山湾南端。这个图像包含了AVIRIS传感器采集的从400至2500纳米波长光谱，有224波段。因为原始图像太大，会导致非常昂贵的计算成本，所以从原图中截取160×250个像素的区域。除去水蒸汽和大气影响的波段后，剩有189波段。具体地讲，有三个端元在该高光谱图像图像内，即“土壤”，“植被”和“水”。本发明步骤(5b)和(6c)中参数λ设置为0.01，步骤(5d)和(6f)中迭代次数c设置为400。

通过本发明的解混方法对上述图解混，其中图3显示了根据GBM模型的莫菲特场用本方法估计的一阶丰度图。图3(a)显示了植被的一阶丰度图，图3(b)显示了水的一阶丰度图，图3(c)显示了土壤的一阶丰度图。

图4显示了根据GBM模型的莫菲特场用本方法估计的二阶丰度图。在图4中，图4(a)表明水体和土壤之间的相互作用发生在沿海地区，图4(b)表明植被和土壤之间的相互作用位于土壤部分,图4(c)表明植被和水之间的相互作用大多发生于植被区域。

通过图3、图4，表明本发明能够精确的显示匀质区域和细节区域的丰度图。

综上，本发明的基于区域分割的高光谱图像自适应解混方法。考虑了线性混合和双线性混合模型共存的情况，其实现主要有：输入高光谱图像；用基于最小错误高光谱信号识别法估计高光谱图像的端元数目；用顶点成分分析算法提取高光谱图像的端元矩阵；用K均值聚类方法对高光谱数据聚类，将图像分割为匀质区域和细节区域；匀质区域采用线性模型，用稀疏约束的非负矩阵分解方法解混，细节区域采用广义双线性模型，用稀疏约束的半非负矩阵分解方法解混。本发明结合了高光谱数据光谱和丰度的特点，获得高光谱图像更精确的表示，提高解混的准确率。并在丰度上添加稀疏性约束条件，克服了半非负矩阵分解算法容易陷入局部最小点的缺点，得到更加准确的丰度，用于高光谱图像的地物识别。

Claims (3)

1.一种基于区域分割的高光谱图像自适应解混方法，其特征在于，包括有如下步骤：

(1)输入高光谱图像，高光谱图像用数据用Y∈R^L×N表示，其中L表示高光谱数据的波段数，N表示高光谱数据样本总个数，R表示实数域；

(2)使用基于最小错误高光谱信号识别法估计高光谱数据Y∈R^L×N的信号子空间，得到信号子空间维数K，即高光谱数据的端元数目；

(3)用顶点成分分析算法提取高光谱图像的端元矩阵A∈R^L×K，其每行向量表示高光谱图像第i个端元的光谱曲线；

(4)用K均值聚类方法对高光谱数据Y∈R^L×N聚类，聚为K+2类，类标分别为1,2,3...,K+2，根据类标，用2*2的窗口扫描全图，若窗口中至少有一个类标不同，则把该2*2窗口对应的像素点归为高光谱图像的细节区域数据其余像素点归为高光谱图像匀质区域数据其中N＝N₁+N₂，N₁表示高光谱图像的匀质区域数据样本个数，N₂表示高光谱图像的细节区域数据样本个数；

(5)将高光谱图像匀质区域数据Y₁和高光谱图像端元矩阵A作为线性光谱混合模型的输入信号，用L_1/2约束的非负矩阵分解方法得到匀质区域数据的系数矩阵即一阶丰度矩阵，其中每一列向量表示第n个像素的丰度向量；

(6)将高光谱图像细节区域数据Y₂和高光谱图像端元矩阵A作为双线性光谱混合模型的输入信号，用稀疏约束的半非负矩阵分解方法得到细节区域数据Y₂对应的一阶丰度矩阵其中每一列向量表示第m个像素的丰度向量，和细节区域数据Y₂对应的二阶丰度矩阵其中每一列向量表示第m个像素的双线性丰度向量；

(7)将高光谱图像匀质区域Y₁的一阶丰度矩阵和高光谱图像细节区域Y₂的一阶丰度矩阵合并为得到为整个高光谱数据的丰度矩阵，细节区域数据Y₂对应的二阶丰度矩阵E即为整个高光谱数据的二阶丰度矩阵，完成高光谱图像解混。

2.根据权利要求1所述的基于区域自适应分割的高光谱图像解混方法，其特征在于，步骤5所述的用L_1/2约束的非负矩阵分解方法得到匀质区域数据Y₁的一阶丰度矩阵X₁，包括有如下步骤：

(5a)求解在高光谱图像匀质区域数据的丰度矩阵X₁的中加入L_1/2范数，得到稀疏约束表达式作为丰度矩阵X₁的稀疏约束项，其中x_1n(k)是高光谱图像匀质区域数据Y₁中在第n个像元对应的第k个端元的丰度；

(5b)将步骤(5a)得到的稀疏约束项添加到非负矩阵分解算法的目标函数中，构成新的目标函数：

$min\frac{1}{2}\left|\right|Y_1-{\mathrm{AX}}_1|{|}_2^2+\mathrm{\λ}|\left|X_1\right|{|}_{1/2}$

约束条件：X₁≥0，1^TX₁＝1^T

其中，λ是稀疏约束正则化参数，X₁≥0，1^TX₁＝1^T是对高光谱图像匀质区域丰度矩阵的“非负”和“和为一”约束；

(5c)对步骤(5b)得到的目标函数用迭代乘法进行优化，匀质区域数据的一阶丰度矩阵X₁的更新公式为：其中(·)^T表示矩阵的转置，.*和./分别表示按元素的乘法和除法；

(5d)重复执行步骤(5c)中匀质区域数据的一阶丰度矩阵X₁的更新公式，直到迭代次数达到设定次数c。

3.根据权利要求1所述的基于区域自适应分割的高光谱图像解混方法，其中步骤(6)用稀疏约束的半非负矩阵分解方法得到细节区域数据Y₂对应的一阶丰度矩阵X₂和二阶丰度矩阵E，按照如下步骤进行：

(6a)在高光谱图像细节区域Y₂采用双线性模型表示如下

Y₂＝AX₂+BE+M

其中，表示双线性区域数据对应的一阶丰度矩阵，其中每一列向量表示第m个像素的丰度向量，是双线性端元矩阵，其中每列向量b_(i,j)＝a_i⊙a_j(i,j∈{1,2,...,K}表示第i个和第j个端元之间的二阶丰度向量，计算如下：

其中⊙表示Hadamard乘积；是双线性区域数据对应的二阶丰度矩阵，其中每一列向量表示第m个像素的双线性丰度向量，表示噪声矩阵；

(6b)求解高光谱图像细节区域数据的丰度矩阵X₂的L_1/2范数，用稀疏约束表达式作为X₂的稀疏约束项，其中x_2m(k)是高光谱图像细节区域数据Y₂中在第m个像元的第k个端元的丰度；

(6c)将步骤(6b)中得到的稀疏约束项添加到步骤(6a)所述双线性模型中，得到双线性模型的最小化目标函数

$\min \left|\right|Y_2-{\mathrm{AX}}_2-BE|{|}_F^2+\mathrm{\λ}|\left|X_2\right|{|}_{1/2}$

约束条件：X₂≥0,1^TX₂＝1^T,0≤E≤X^*

其中的每个元素用(X^*)_(i,j)n＝(X₂)_in(X₂)_jn(n∈{1,2,...,N₂}计算，||·||_F表示F范数，X₂≥0，1^TX₂＝1^T分别是对高光谱图像细节区域丰度矩阵的“非负”和“和为一”约束；

(6d)引入两个变量Y₂₁＝Y₂-AX₂，Y₂₂＝Y₂-BE，可将步骤(6a)中的双线性模型写成如下形式Y₂₁＝AX₂+M和Y₂₂＝ΒΕ+M，然后分别关于A和B迭代的优化和

(6e)用迭代乘法优化步骤(6c)中的和得到细节区域数据Y₂对应的一阶丰度矩阵X₂的更新公式和细节区域数据Y₂对应的二阶丰度矩阵E的更新公式

(6f)重复执行步骤(6e)得到的细节区域数据Y₂对应的一阶丰度矩阵X₂的更新公式和细节区域数据Y₂对应的二阶丰度矩阵E的更新公式，直到迭代次数达到设定次数c。