CN103679210A - 基于高光谱图像解混的地物识别方法 - Google Patents

Abstract

本发明公开了一种基于高光谱图像解混的地物识别方法，主要解决现有方法判断混合像素点所属地物类别不准确的问题。其实现步骤为：输入一幅高光谱图像，将该高光谱图像中的混合像素点按列排成一个矩阵，组成数据矩阵；用数据矩阵的流形约束，丰度矩阵的稀疏约束和端元矩阵的平滑约束构成的约束项，加入到NMF算法的目标函数里，构成新的目标函数；对新的目标函数进行优化解混，得到该高光谱图像解混之后的端元矩阵和丰度矩阵；根据解混之后的端元矩阵和丰度矩阵判断该高光谱图像中所有混合像素点的地物类别。本发明能提高解混得到的端元值和丰度值的精度，从而提高了高光谱图像地物识别的精度，可用于目标跟踪。

Description

基于高光谱图像解混的地物识别方法

技术领域

本发明属于遥感图像处理技术领域，是一种基于高光谱图像解混的地物识别方法，该方法可用于高光谱图像的分析，将一个混合像素点分解为端元和对应的丰度值。 

背景技术

高光谱图像是利用成像光谱仪对同一地表区域的几十乃至几百个波段同时成像而获得的三维图像，由二维空间信息和一维光谱信息组成。利用这些丰富的光谱信息对地物进行细分和鉴别，在多领域得到了广泛应用。因为高光谱图像的光谱传感器的光谱分辨率高，所以形成的谱段多，但是每个谱段接受的能量小，所以只能提高接受光谱的地面面积，即降低空间分辨率。所以由于高光谱图像的空间分辨率不高还有自然界地物的复杂性从而形成混合像素点。混合像素点的普遍存在不仅影响地物的识别和分类精度，而且是遥感技术向定量化发展的重要障碍。因此如何有效地进行混合像素点分解是高光谱遥感图像应用的关键问题之一。 

高光谱图像中混合像素点的模型一般采用线型混合模型，它的优点是算法简单，物理意义明确。该模型的数学过程简要描述如下：一个具有L个谱段的像素点表示为X_ij∈R^L×1，具有P个端元的端元矩阵表示为M∈R^L×P，M对应的丰度矩阵表示为S_ij∈R^P×1。则有：X_ij＝MS_ij+n，其中n为噪声。在实际环境中该模型受两个条件的限制：①M_up≥0,(1≤u≤L,1≤p≤P)② 

上面的两个式子分别表示光谱的能量不存在负值和混合的能量的大小是一定的，不可能无限大。上述的模型和限制条件都满足非负矩阵分解（Nonnegative Matrix Factorization，NMF）的数学模型，所以可以用NMF算法来进行解混。 

目前提出的基于非负矩阵分解（Nonnegative Matrix Factorization，NMF）的高光谱图像的解混方法都是在NMF算法的目标函数上加上正则项。因为该方法 没有充分考虑高光谱图像的特性，使得解混效果差，从而导致地物识别的精度不高。 

发明内容

本发明的目的在于针对现有方法的不足，提出一种基于高光谱图像解混的地物识别方法，本方法将高光谱图像的端元矩阵的平滑性，丰度矩阵的稀疏性和数据矩阵的流形假设结构结合在一起，使得高光谱图像解混效果更好，从而使得地物识别的精度更高。 

实现本方法的技术方案是：输入一幅高光谱图像，将该高光谱图像中的混合像素点按列排成一个矩阵，组成数据矩阵，用数据矩阵的流形约束，丰度矩阵的稀疏约束和端元矩阵的平滑约束构成的约束项，加入到NMF算法的目标函数里，构成新的目标函数，然后对这个新的目标函数进行优化解混，得到该高光谱图像的端元矩阵和丰度矩阵，然后根据该解混出的端元矩阵和丰度矩阵判断该高光谱图像中所有混合像素点的地物类别。具体步骤包括如下： 

（1）输入一幅高光谱图像X∈R^M×N×L，并将该高光谱图像中混合像素点X_ij∈R^1×L按列排列，构成数据矩阵Z∈R^L×B，其中M和N为二维图像的行和列，i和为j为二维图像的横坐标和纵坐标，L为谱段数，B为高光谱图像中混合像素点总数，B＝M×N，R表示实数集合； 

（2）根据流形假设理论，构造数据矩阵Z的流形约束项： 

$\frac{1}{2}\underset{i,j=1}{\overset{B}{\mathrm{\Σ}}}{\left|\right|s_i-s_j\left|\right|}^2{W}_{\mathrm{ij}}=\underset{i=1}{\overset{B}{\mathrm{\Σ}}}{s}_i^T{s}_i{D}_{\mathrm{ii}}-\underset{i,j=1}{\overset{B}{\mathrm{\Σ}}}{s}_i^T{s}_j{W}_{\mathrm{ij}}$

$=\mathrm{Tr}\left({\mathrm{SYS}}^T\right)$

其中z_i是Z的第i列，z_j是Z的第j列，且z_i是z_j的k个近邻中的一个，S是丰度矩阵，s_i是S的第i列，W是Z的权值矩阵，W_ij是W的一个元素， 

为z_i和z_j的权值，Tr(·)表示矩阵的迹，T表示矩阵的转置，D是Z的对角线权值矩阵，D_ii是D对角线上的一个元素，D_ii＝Σ_jW_ij，Y是流 形系数矩阵，Y＝D-W； 

（3）根据高光谱图像成像理论，在丰度矩阵S中加入L_1/2范数，得到稀疏约束表达式‖S‖_1/2，以作为丰度矩阵S的稀疏约束项； 

（4）根据高光谱图像成像理论，在端元矩阵M中加入Frobenius范数，得到平滑约束表达式 

以作为端元矩阵M的平滑约束项； 

（5）将步骤(2)-(4)得到的三个约束项添加到NMF算法的目标函数 

里，以构成新的目标函数： 

$f\text{'}(M,S)=\frac{1}{2}{\left|\right|Z-\mathrm{MS}\left|\right|}_F^2+\mathrm{\α}{\left|\right|S\left|\right|}_{1/2}+\mathrm{\β}{\left|\right|M\left|\right|}_F^2+\mathrm{\γTr}\left({\mathrm{SYS}}^T\right)$

其中，α为丰度矩阵S的稀疏约束正则参数，β为端元矩阵M的平滑约束正则参数，γ为数据矩阵Z的流形约束正则参数； 

（6）对步骤（5）得到的目标函数f’（M，S）用迭代乘法进行优化求解，得到高光谱图像X∈R^M×N×L的端元矩阵M和丰度矩阵S； 

（7）将上述高光谱图像X∈R^M×N×L中混合像素点X_ij用步骤（6）求解得到的端元矩阵M和丰度向量s_i表示，即混合像素点X_ij＝Ms_i； 

（8）根据高光谱图像统计分布理论，由步骤（7）中的丰度向量s_i对混合像素点X_ij进行地物类别判断，即当max(s_i)=s_ai时，则判混合像素点X_ij属于第a类，得到该混合像素点的类别标签为v_ij=a，其中max(·)表示取向量中的最大值，a＝1,2,...,P表示该高光谱图像中相应的地物类别编号，P表示该高光谱图像中地物类别总数，s_ai是s_i的第a个元素； 

（9）对上述高光谱图像X∈R^M×N×L中所有混合像素点用步骤（8）的操作进行地物类别判断，得到该高光谱图像X∈R^M×N×L的地物类别矩阵V∈R^M×N。 

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

1.本发明相比现有的方法更充分的考虑了高光谱图像的结构特性，将数据矩阵的流形约束加入到解混，与未加入流形约束的NMF算法，CNMF算法和分段光滑NMF算法相比，得到的丰度矩阵的精度更高，使得判断混合像素点所属地物类别更准确，从而达到较高的地物识别精度。 

2.本发明同时将端元矩阵的平滑约束和丰度矩阵的稀疏约束加入到解混，与现有的GLNMF算法相比，更充分考虑了高光谱图像特性，大大提高了端元矩阵的精度，使得计算出的地物类别更接近真实地物类别，更加符合实际应用的要求。 

附图说明

图1是本发明的实现流程图； 

图2是本发明使用的高光谱图像示意图； 

图3是本发明与现有几种算法解混出来的端元的光谱曲线和实际光谱曲线的对比示意图； 

图4是本发明与现有几种算法解混出的端元的光谱角距离SAD和解混出的丰度的均方根误差RMSE的平均30次误差棒对比图。 

具体实施方式

参照图1，本发明的具体实施步骤如下： 

步骤1：输入高光谱图像，构建数据矩阵，获得该高光谱图像的真实的地物类别矩阵，真实的端元矩阵和真实的丰度矩阵。 

1.1)输入如图2所示的高光谱图像，该图像大小为145×145，共有16类地物，该图像中的每一个混合像素点可看做是由200个谱段的光谱信息构成的光谱向量； 

1.2)将该高光谱图像X∈R^M×N×L中混合像素点X_ij∈R^1×L按列排列，构成数据矩阵Z∈R^L×B，其中，M和N为二维图像的行和列，i和j分别为二维图像的横坐标和纵坐标，L为谱段数，P为地物类别数，B为高光谱图像中混合像素点总数，B＝M×N，R表示实数集合，在本实例中L取值为200，B取值为21025，P取值为16； 

1.3)获得该高光谱图像的真实类别矩阵为V∈R^M×N，真实端元为M∈R^L×P， 真实丰度为S∈R^P×B。 

步骤2：构造约束项： 

1.1）根据流形假设理论，构造数据矩阵Z的流形约束项： 

1.1a）用热核函数计算数据矩阵Z的第i列z_i和数据矩阵Z的第j列z_j的权值W_ij： 

$W_{\mathrm{ij}}=e^{-\frac{{\left|\right|z_i-z_j\left|\right|}^2}{\mathrm{\σ}}},$

其中，σ为热核参数，取值为1； 

1.1b）用L₂范数计算丰度矩阵S的第i列s_i和丰度矩阵S的第j列s_j的距离： 

‖s_i-s_j‖²； 

1.1c）根据流形假设理论，将步骤1.1a）的式子和步骤1.1b）的式子结合起来构成流形约束项： 

$\frac{1}{2}\underset{i,j=1}{\overset{B}{\mathrm{\Σ}}}{\left|\right|s_i-s_j\left|\right|}^2{W}_{\mathrm{ij}}=\underset{i=1}{\overset{B}{\mathrm{\Σ}}}{s}_i^T{s}_i{D}_{\mathrm{ii}}-\underset{i,j=1}{\overset{B}{\mathrm{\Σ}}}{s}_i^T{s}_j{W}_{\mathrm{ij}}$

$=\mathrm{Tr}\left({\mathrm{SYS}}^T\right)$

其中z_i是z_j的k个近邻中的一个，W是Z的权值矩阵，Tr(·)表示矩阵的迹，T表示矩阵的转置，D是Z的对角线权值矩阵，D_ii是D对角线上的一个元素，D_ii＝Σ_jW_ij，Y是流形系数矩阵，Y＝D-W； 

1.2）根据高光谱图像成像理论，在丰度矩阵S中加入L_1/2范数，得到稀疏约束表达式‖S‖_1/2，以作为丰度矩阵S的稀疏约束项； 

1.3）根据高光谱图像成像理论，在端元矩阵M中加入Frobenius范数，得到平滑约束表达式 

以作为端元矩阵M的平滑约束项。 

步骤3：构造新的目标函数 

1.1）将步骤2得到的三个约束项添加入到NMF算法的目标函数 

里，以构成新的目标函数： 

$f\text{'}(M,S)=\frac{1}{2}{\left|\right|Z-\mathrm{MS}\left|\right|}_F^2+\mathrm{\α}{\left|\right|S\left|\right|}_{1/2}+\mathrm{\β}{\left|\right|M\left|\right|}_F^2+\mathrm{\γTr}\left({\mathrm{SYS}}^T\right),$

其中，α为丰度矩阵S的稀疏约束正则参数，取值为2，β为端元矩阵M的平滑约束正则参数，取值为1，γ为数据矩阵Z的流形约束正则参数，取值为0.6。 

步骤4：对步骤3得到的新的目标函数f’（M，S）用迭代乘法进行优化求解，得到端元矩阵M和丰度矩阵S。 

1.1）用[0,1]之间随机数初始化端元矩阵M和丰度矩阵S，将S的每列进行归一化处理； 

1.2）输入端元矩阵M和丰度矩阵S； 

1.3）通过对新的目标函数f’（M，S）中的端元矩阵M和丰度矩阵S求导，分别得到端元矩阵M和丰度矩阵S的迭代乘法的计算公式： 

M’=M.*(ZS^T-βM)./MSS^T， 

$S.*(M^{T}Z+2\mathrm{\γSW})./(M^T\mathrm{MS}+\frac{1}{2}\mathrm{\α}{S}^{-\frac{1}{2}}+2\mathrm{\γSD})$

其中，（.*）表示矩阵的元素乘法，（./）表示矩阵的元素除法； 

1.4）用步骤1.3）的式子计算出的新的丰度矩阵S’和新的端元矩阵M’，将新的丰度矩阵S’和新的端元矩阵M’作为步骤1.2）的输入； 

1.5）重复执行步骤1.2）～1.4）共n次，输出最终需要求出的端元矩阵M和丰度矩阵S，计算结束，其中，n为执行次数，取值为3000。 

步骤5：判断步骤1输入的高光谱图像X∈R^M×N×L中所有混合像素点的地物类别。 

1.1）对步骤1输入的高光谱图像X∈R^M×N×L中的混合像素点X_ij用步骤3解混得到的端元矩阵M和丰度向量s_i表示，即X_ij=Ms_i； 

1.2）根据高光谱图像统计分布理论，由步骤1.1）中得到的丰度向量s_i对 混合像素点X_ij进行地物类别判断，即当max(s_i)=s_ai时，则判混合像素点X_ij属于第a类，得到该混合像素点的类别标签为v_ij=a，其中max(·)表示取向量中的最大值，a＝1,2,...,P表示该高光谱图像中相应的地物类别编号，P表示该高光谱图像中地物类别总数，s_ai是s_i的第a个元素； 

1.3）对该高光谱图像X∈R^M×N×L中所有混合像素点用步骤1.2）的操作进行地物类别判断，得到该高光谱图像X∈R^M×N×L的解混之后的地物类别矩阵V∈R^M×N。 

本发明的效果通过以下仿真实验进一步说明： 

（1）实验仿真条件： 

本实验使用的高光谱图像是典型的AVIRIS高光谱图像：取自1992年6月拍摄的美国印第安纳州西北部印第安遥感试验区，地物类别共计16类，图像的大小为145×145。原始数据共有220个谱段，除去被噪声污染和水域污染的20个谱段，仅保留剩下的200个谱段。本实验在CPU为Intel(R)Core(TM)i5-2450、主频2.5GHz，内存为4G的WINDOWS7***上采用软件MATLAB2009a进行仿真。 

（2）评价指标： 

1.1）对于端元，用光谱角距离SAD来比较真实端元M_t和估计端元M_t的相似性。光谱角距离SAD的值越小则两个端元的光谱曲线越接近。对于丰度，用均方根误差RMSE来比较真实丰度S_t和估计丰度S_t之间的差异。均方根误差RMSE的值越小则两个丰度值越接近。上述两个评价标准的公式分别为： 

${\mathrm{SAD}}_t={\cos }^{-1}\left(\frac{M_t^T{M}_t}{\left|\right|M_t^T\left|\right|\left|\right|M_t\left|\right|}\right),$

${\mathrm{RMSE}}_t={\left(\frac{1}{B}{|S_t-S_t|}^2\right)}^{\frac{1}{2}},$

其中，M_t为真实端元矩阵M的第t列，M_t为估计端元矩阵M的第t列，S_t为真实丰度矩阵S的第t列，S_t为估计丰度矩阵S的第t列，t表示地物的类别， t=1,2,...,P，P为地物类别总数； 

（3）实验仿真内容： 

实验一 

利用本发明对步骤（1）中所述高光谱图像进行优化解混，得到端元矩阵M和丰度矩阵S，再分别用NMF算法，CNMF算法，GLNMF算法和分段光滑NMF算法对上述高光谱图像进行优化解混，与本发明解混得到的结果进行对比，实验结果如图3所示，其中: 

图3（a）是用NMF算法解混出的一种地物氯黄晶光谱曲线与实际地物氯黄晶光谱曲线的对比示意图； 

图3（b）是用CNMF算法解混出的一种地物氯黄晶光谱曲线与实际地物氯黄晶光谱曲线的对比示意图； 

图3（c）是用GLNMF算法解混出的一种地物氯黄晶光谱曲线与实际地物氯黄晶光谱曲线的对比示意图； 

图3（d）是用分段光滑NMF算法解混出的一种地物氯黄晶光谱曲线与实际地物氯黄晶光谱曲线的对比示意图； 

图3（e）是用本发明的算法解混出的一种地物氯黄晶光谱曲线与实际地物氯黄晶光谱曲线的对比示意图； 

由图3可见，本发明相比于其他现有方法，得到的地物的光谱曲线更接近真实地物地物的光谱曲线。 

实验二 

利用步骤（2）对实验一得到的端元矩阵M和丰度矩阵S分别计算端元矩阵M的光谱角距离SAD值和丰度矩阵S的均方根误差RMSE值，实验结果如图4所示，其中: 

图4（a）是本发明与上述4种现有算法解混出的端元的光谱角距离SAD值的平均30次结果误差棒对比图； 

图4（b）是本发明与上述4种算法解混出的丰度的均方根误差RMSE值的平均30次结果误差棒对比图； 

由图4可见，本发明相比于其他现有方法，解混得到的端元值和丰度值的精度更高。 

实验三 

利用实验一得到的端元矩阵M和丰度矩阵S对该高光谱图像中所有混合像素点进行类别判断，得到解混之后的类别矩阵V，再利用支撑矢量机SVM对解混之后得到的类别矩阵V和真实类别矩阵V进行地物识别精度的测量，将本发明测量所得的地物识别精度与上述4种现有方法测量所得的地物识别精度进行对比，如表1所示。 

表1 地物识别精度数值指标对比 

从表1可见，利用本发明进行地物识别得到的地物识别精度已经明显的优于上述4种现有算法的地物识别精度。 

综上所述，本发明可以大幅度提高高光谱图像中端元值和丰度值的精度，从而更好地提高了地物识别的精度，所以本发明作为一种基于高光谱图像解混的地物识别方法在高光谱图像识别领域具有广阔的应用前景。 

Claims (3)

1.一种基于高光谱图像解混的地物识别方法，包括如下步骤： 

（1）输入一幅高光谱图像X∈R^M×N×L，并将该高光谱图像中混合像素点X_ij∈R^1×L按列排列，构成数据矩阵Z∈R^L×B，其中M和N为二维图像的行和列，i和为j为二维图像的横坐标和纵坐标，L为谱段数，B为高光谱图像中混合像素点总数，B＝M×N，R表示实数集合； 

（2）根据流形假设理论，构造数据矩阵Z的流形约束项： 

其中z_i是Z的第i列，z_j是Z的第j列，且z_i是z_j的k个近邻中的一个，S是丰度矩阵，s_i是S的第i列，s_j是S的第j列，，W是Z的权值矩阵，W_ij是W的一个元素，   

为z_i和z_j的权值，Tr(·)表示矩阵的迹，T表示矩阵的转置，D是Z的对角线权值矩阵，D_ii是D对角线上的一个元素，D_ii＝Σ_jW_ij，Y是流形系数矩阵，Y＝D-W； 

（3）根据高光谱图像成像理论，在丰度矩阵S中加入L_1/2范数，得到稀疏约束表达式‖S‖_1/2，以作为丰度矩阵S的稀疏约束项； 

（4）根据高光谱图像成像理论，在端元矩阵M中加入Frobenius范数，得到平滑约束表达式   

以作为端元矩阵M的平滑约束项； 

（5）将步骤(2)-(4)得到的三个约束项添加到NMF算法的目标函数    

里，以构成新的目标函数： 

其中，α为丰度矩阵S的稀疏约束正则参数，β为端元矩阵M的平滑约束正则参数，γ为数据矩阵Z的流形约束正则参数； 

（6）对步骤（5）得到的目标函数f’（M，S）用迭代乘法进行优化求解，得到高光谱图像X∈R^M×N×L的端元矩阵M和丰度矩阵S； 

（7）将上述高光谱图像X∈R^M×N×L中混合像素点X_ij用步骤（6）求解得到的端元矩阵M和丰度向量s_i表示，即混合像素点X_ij＝Ms_i； 

（8）根据高光谱图像统计分布理论，由步骤（7）中的丰度向量s_i对混合像素点X_ij进行地物类别判断，即当max(s_i)=s_ai时，则判混合像素点X_ij属于第a类，得到该混合像素点的类别标签为v_ij=a，其中max(·)表示取向量中的最大值，a＝1,2,...,P表示该高光谱图像中相应的地物类别编号，P表示该高光谱图像中地物类别总数，s_ai是s_i的第a个元素； 

（9）对上述高光谱图像X∈R^M×N×L中所有混合像素点用步骤（8）的操作进行地物类别判断，得到该高光谱图像X∈R^M×N×L的地物类别矩阵V∈R^M×N。 

2.根据权利要求书1所述的基于高光谱图像解混的地物识别方法，其中步骤（2）所述的根据流形假设理论，构造数据矩阵Z的流形约束项，按如下步骤进行： 

2a）用热核函数计算数据矩阵Z的第i列z_i和数据矩阵Z的第j列z_j的权值W_ij： 

其中，σ为热核参数，取值为1； 

2b）用L₂范数计算丰度矩阵S的第i列s_i和丰度矩阵S的第j列s_j的距离： 

‖s_i-s_j‖²； 

2c）根据流形假设理论，将步骤2a）的式子和步骤2b）的式子结合起来构成流形约束项：   

其中，B为高光谱图像中混合像素点的总数。 

3.根据权利要求书1所述的基于高光谱图像解混的地物识别方法，其中所述步骤（6）中用迭代乘法优化求解目标函数f’（M，S），按如下步骤进行： 

3a)用（0,1）之间随机数初始化端元矩阵M和丰度矩阵S，将丰度矩阵S的每列进行归一化处理；输入端元矩阵M和丰度矩阵S； 

3b)对目标函数f’（M，S）=0中的丰度矩阵S进行求导，得到丰度矩阵S的迭代乘法计算公式： 

其中（.*）表示矩阵的元素乘法，（./）表示矩阵的元素除法，W是Z的权值矩阵，D是Z的对角线权值矩阵，Z为数据矩阵，T表示矩阵的转置，α为丰度矩阵S的稀疏约束正则参数，γ为数据矩阵Z的流形约束正则参数； 

3c)对目标函数f’（M，S）=0中的端元矩阵M进行求导，得到端元矩阵M的迭代乘法计算公式： 

M.*(ZS^T-βM)./MSS^T

其中，β为端元矩阵M的平滑约束正则参数； 

3d)用步骤3b)和步骤3c)的式子计算出的新的丰度矩阵S’和新的端元矩阵M’，将丰度矩阵S’和新的端元矩阵M’作为步骤3a)的输入； 

3e)重复执行步骤3a)～3d)n次，输出最终需要求出的端元矩阵M和丰度矩阵S，则计算结束，其中，n为执行次数，取值为3000。 

Images