CN102314685B - 一种基于随机投影的高光谱图像稀疏解混方法 - Google Patents

Abstract

一种基于随机投影的高光谱图像稀疏解混方法，它有四大步骤：一、计算机在MATLAB R2008b环境下读取数据；二、计算机对高光谱图像数据和高光谱库数据进行随机投影；三、构建稀疏解混的目标函数，使用分裂Bregman算法优化目标函数求极值，直至达到收敛停止条件。四、设定合适的阈值处理丰度矩阵，获得最终的丰度图和端元。本发明利用了高光谱数据库来选择端元，克服了以往算法所求出的端元与标准高数据库中的纯物质光谱无法严密对应的缺点；并使用随机投影技术对原始数据进行降维，从而达到了节省内存和减少计算量的效果。本发明实现了对高光谱图像快速定量分析，它在高光谱遥感图像分析领域里具有实用价值和广阔的应用前景。

Description

一种基于随机投影的高光谱图像稀疏解混方法

(一)技术领域：

本发明涉及一种基于随机投影的高光谱图像稀疏解混方法，属于高光谱遥感图像分析技术领域。

(二)背景技术：

在过去三十年间，随着成像光谱技术(imaging spectroscopy)的不断发展，在飞机或卫星平台上搭载的成像光谱仪(imaging spectrometer)采集得到的遥感图像(remote sensingimages)包含了越来越丰富的空间、辐射和光谱信息，从而为地表物质的信息提取和目标检测提供了一个强有力的手段。光谱分辨率作为遥感技术发展的一个重要指标，目前已由多光谱发展到高光谱，正向超高光谱发展。

一般来说，在同一环境下，由同一光谱传感器检测得到的物质光谱具有不同物质不同光谱特征的特点。根据这一特点，可以对高光谱图像的每一个像素点进行物质分析。然而由于遥感探测仪器较低的空间分辨率限制以及自然界地物的复杂多样性，使得单个像元点的光谱特征反映的不一定只是一种物质的特性，而可能是几种不同物质光谱的混合，这样的特征像素点被称为混合像元(mixed pixel)。为了提高遥感应用的精度，就必须解决混合像元的分解问题，使遥感应用由像元级达到亚像元级。进入像元内部，将混合像元分解为不同的“基本组成单元”，或称“端元”(endmember)，并求得这些基本组分所占的比例，或称“丰度”(abundancefractions)，这就是所谓的“光谱解混”过程(spectral unmixing)(童庆禧，张兵，郑兰芬.高光谱遥感——原理、技术与应用[M].高等教育出版社2006)。光谱解混分析实现了高光谱数据的定量分析。

传统的高光谱图像解混方法是：首先利用端元选择(endmember selection)技术获取存在于高光谱图像中的成份光谱；然后通过混合像元分解(mixed pixel decomposition)技术确定各个端元在混合像元中所占的比例。由于这类方法是以端元光谱为条件建立的，因此必须先以监督或非监督的方法找到所谓的端元光谱。近年来，随着盲信号分离(blind source separation，BSS)和非负矩阵分解(Nonnegative Matrix Factorization，NMF)技术的兴起和发展，非监督的光谱解混技术逐渐受到遥感学者们的重视。它是指在端元信息完全未知的情况下，直接从遥感图像本身入手，根据混合像元的光谱模型以及约束条件等信息，利用非监督的信号处理方法得到端元光谱及其组分信息。目前，非监督的光谱解混算法主要集中在两个方向：基于凸集理论的单纯形几何学方法和基于盲信号分离(BSS)技术的独立成分分析(Independent ComponentAnalysis，ICA)及非负矩阵分解(nonnegative matrix factorization，NMF)方法。这两种算法在混合像元分解上都可以获得较为精确的分解结果，但同时也面临一些需要解决的问题。例如，单纯形方法往往需要假设高光谱图像数据中的每一个端元均存在纯像元，否则分解精度会受到显著影响。非负矩阵分解方法存在很多局部极小，这会造成分解的不唯一，极大的影响了算法的性能。

随着高光谱成像仪的广泛应用，高光谱数据库得到不断完善，越来越多纯物质的光谱特征得到测量，一种基于高光谱数据库的半监督解混算法成为研究的热点，即高光谱图像稀疏解混方法。它将高光谱数据库所有的物质光谱作为假想的端元，真正的端元在光谱库中所占的比例是很小的。高光谱图像稀疏解混的目的就是要找到每个像元光谱在高光谱数据库中的稀疏表达。可以利用假想端元的丰度稀疏性作为正则化项，求解所有假想端元所对应的丰度，并在光谱库中挑选出真正的端元。稀疏解混方法能够准确直接的找到端元，避免了前两种方法求解端元不精确的问题。

三种高光谱图像解混技术面对的共同问题是海量的高光谱数据造成解混过程中计算机内存不足以及计算时间较长。针对上述问题，本发明采用随机投影算法将海量的高光谱数据投影在较低维的空间，以投影后的较少量数据代替原始数据进行稀疏解混处理，并采用分裂Bregman算法对高光谱图像所有的像素点进行整体处理，来代替传统的逐个像素进行解混的方式，进一步减少时间消耗，取得了较好的效果。

(三)发明内容：

1、目的：本发明的目的是提供一种基于随机投影的高光谱图像稀疏解混方法，此方法能够对原始海量高光谱数据进行降维处理，减少内存消耗和计算时间；并利用高光谱数据库，为丰度加以稀疏性约束，进而求得高光谱图像每个像元的组成物质以及其含量百分比，实现高光谱图像的定量分析。

2、技术方案：本发明是通过以下技术方案实现的：

本发明一种基于随机投影的高光谱图像稀疏解混方法，包含以下几个步骤：

步骤一：用计算机读取数据：计算机在MATLAB R2008b环境下读取高光谱图像数据，数据来源于成像光谱仪采集到的遥感图像，得到数据立方体，高光谱图像数据应去除被水汽吸收的波段和信噪比较低的波段。将高光谱图像数据逐像素点排列，得到原始的高光谱图像矩阵V。假设高光谱图像有m个波段，共有n个像素点，V＝[V₁，v₂，...，v_n]，v_i是高光谱图像第i个像素点的光谱列向量，是一个m维向量。读取现有的高光谱数据库数据，选择光谱库中大量的纯物质光谱数据构建光谱库矩阵S。假设纯物质光谱的数量为q，S＝[s₁，s₂，...，s_q]，s_i是光谱库中第i个纯物质的光谱列向量，同样是一个m维向量。

步骤二：高光谱数据随机投影：用计算机读取数据后，对高光谱图像数据和高光谱库数据进行随机投影，就是用一个包含的所有元素符合正态分布的矩阵R_d与矩阵V，S分别相乘。假设矩阵R_d的大小为d×m，d＜m，得到降维后的数据R_dV和R_dS，其矩阵大小分别为d×n，d×q，代替原始的高光谱图像数据V和高光谱数据库S参与下述步骤三的计算，其中d是随机投影后光谱波段的个数。

步骤三：高光谱图像稀疏解混：数据进行随机投影后，接下来就要对降维后的高光谱图像数据进行光谱解混。假设高光谱数据库中含有足够丰富的纯物质光谱，那么高光谱图像中含有的端元只占极少一部分，也就是说高光谱图像中每个像素点的光谱曲线是由高光谱数据中少数的纯物质的光谱曲线以某种方式组合构成，这体现了端元的稀疏性表达。端元在高光谱数据库中是稀疏的，他们所对应的丰度也具有稀疏性，即丰度矩阵是稀疏的，以丰度的稀疏性作为正则化项，符合实际的物理意义。在这里，采用线性混合模型，即将一条检测得到的光谱曲线，分解为光谱库中纯物质光谱线性组合的形式，其系数即为相应丰度。假设丰度矩阵为H，其大小为q×n，其所有元素满足非负性。线性混合模型如下表示：

V＝SH

对原始高光谱图像稀疏解混的目标函数为：

$\begin{array}{cc}\underset{H}{\min }{\left|\right|H\left|\right|}_1& s.t.V=\mathrm{SH},H\≥0\end{array}$

||H||₁表示H包含的所有元素绝对值的和。因为丰度具有非负性，||H||₁越小，丰度矩阵H含有的非零元素值将有可能越少，其表现出的稀疏性就较强。将上式转变为无约束优化问题：

$\underset{H}{\min }\frac{1}{2}{\left|\right|V-\mathrm{SH}\left|\right|}_2^2+\mathrm{\λ}{\left|\right|H\left|\right|}_1$

在上式目标函数中，l₂项是拟合项，体现的是解混的结果与原高光谱图像的拟合程度；l₁即稀疏性正则化项；λ是正则化参数，调节着拟合项与正则化项的权重，其大小需要根据实际数据进行设置。因为V，S数据量很大，在迭代求解的过程中占用内存较大，并且计算时间较长。本发明采用随机投影算法，能够有效保持上式中的l₂拟合项，实现数据量压缩。随机投影后的目标函数为：

$\underset{H}{\min }\frac{1}{2}{\left|\right|R_{d}V-R_d\mathrm{SH}\left|\right|}_2^2+\mathrm{\λ}{\left|\right|H\left|\right|}_1$

针对上述目标函数，采用分裂Bregman算法求解上述函数极值，从而得到该优化模型的最优解——最优丰度矩阵H。

引入辅助变量D，并且D＝H，上式目标函数转化为：

$\begin{array}{cc}\underset{H,d}{\min }\frac{1}{2}{\left|\right|R_{d}V-R_d\mathrm{SH}\left|\right|}_2^2+\mathrm{\λ}{\left|\right|D\left|\right|}_1& s.t.D=H\end{array}$

转化为无约束形式为：

$\underset{H,d}{\min }{\mathrm{\λ}{\left|\right|D\left|\right|}_1+\frac{\mathrm{\μ}}{2}\left|\right|R_{d}V-R_d\mathrm{SH}\left|\right|}_2^2+\frac{1}{2}{\left|\right|D-H\left|\right|}_2^2$

其中μ是正则化参数，其大小根据实际的高光谱图像数据设置。引入误差反馈变量B，上述优化问题转变为迭代求解过程，如下：

$\left\{\begin{array}{c}{H}^{k+1}=\arg \min \{\frac{\mathrm{\μ}}{2}{\left|\right|R_{d}V-R_d\mathrm{SH}\left|\right|}_2^2+\frac{1}{2}{\left|\right|D^k-H-B^k\left|\right|}_2^2\}\\ D^{k+1}=\arg \min \{\mathrm{\λ}{\left|\right|D\left|\right|}_1+\frac{1}{2}{\left|\right|D-H^{k+1}-B^k\left|\right|}_2^2\}\end{array}\right.$

具体求解步骤如下：

(1)初始化丰度矩阵H以及辅助变量矩阵B，D，V_d，

H⁰＝B⁰＝D⁰＝0；V_d ⁰＝R_dV

(2)迭代变量H，D，B，V_d，交替迭代过程如下：

$\left\{\begin{array}{c}{H}^{k+1}={[\mathrm{\μ}{\left(R_{d}S\right)}^T\left(R_{d}S\right)+I]}^{-1}[\mathrm{\μ}{\left(R_{d}S\right)}^T{V_d}^k-B^k+D^k]\\ D^{k+1}=\max (H^{k+1}+B^k-\mathrm{\λ},0)\\ B^{k+1}=B^k+H^{k+1}-D^{k+1}\\ {V_d}^{k+1}={V_d}^k+R_{d}V-R_d{\mathrm{SH}}^{k+1}\end{array}\right.$

这里k为迭代次数，V_d ^k+1的迭代过程可以加速收敛。

(3)停止迭代条件。当丰度矩阵H相邻两次迭代的相对误差不大时，迭代停止。迭代停止条件为：

$\frac{{\left|\right|H^{k+1}-H^k\left|\right|}_2^2}{{\left|\right|H^k\left|\right|}_2^2}<\mathrm{tol}$

本发明中，取tol＝10^-6，若迭代步数少于10次就能收敛，使用迭代次数控制迭代过程，设迭代次数为1000步。

步骤四：获取丰度图及真正的端元：在获得丰度矩阵后，设置合适的阈值，丰度矩阵中小于阈值的元素设为零，不小于阈值的元素不做处理。找出丰度矩阵中含有非零元素的行标号，丰度矩阵的行标号与光谱库的列标号对应。取出光谱矩阵中列标号所对应的列，即为真正的端元光谱；取出丰度矩阵中行标号所对应的行，即为端元所对应的丰度图。在丰度图中颜色较亮的区域表示端元所占的比例较大，较暗的区域对应的比例小，颜色为黑色的区域表示不含有该端元。

3、优点及功效：本发明的优点是：它利用现有的光谱数据库，通过对丰度的稀疏性正则化约束，能够同时实现端元的选择和丰度图的获取。在计算过程中，对原始的海量高光谱图像数据和拥有大量纯物质光谱的光谱库数据进行随机投影，大大减少了后续计算的数据维数，从而达到了节省内存以及减少计算量的效果；在对目标函数求极值的过程中，采用分裂Bregman算法对矩阵进行操作，而非传统算法中对向量进行重复性操作，进一步减少了计算耗时。

(四)附图说明：

图1本发明所述方法流程图

图2(a)-(h)仿真哈勃望远镜的高光谱图像包含的8种端元

图3(a)-(h)对仿真纯像元哈勃望远镜的高光谱图像进行解混所求得的端元

图4(a)-(h)对仿真混合像元哈勃望远镜的高光谱图像进行解混所求得的端元

仿真的高光谱图像来源于美国航空航天局约翰逊航天中心(NASAJohnson Space Center)的哈勃望远镜3D模型，其包含的8种端元光谱(如图2(a)-(h)所示)来自NASA高光谱数据库，物质名称依次是螺栓、剥离铜、太阳能电池、哈勃铝、哈勃绿胶水、哈勃蜂窝顶、哈勃蜂窝边以及黑色橡胶边。仿真高光谱图像大小为177*193*100，波段范围为0.4μm-2.5μm，共100波段。构建高光谱数据库的数据除了8种端元外，其余的物质光谱来自USGS数据库，包括矿石、人造物、植物等591物质。纯像元仿真的高光谱图像每个像素至多含有一种物质，其可视化图和八种组成物质所对应的空间分布如图3所示。将纯像元高光谱图像的空间分辨率降低为原来的1/4，即可得到混合像元的仿真高光谱图像。

本发明分别对纯像元和混合像元仿真高光谱图像进行稀疏解混，对于解混的结果——端元和丰度分别采用谱信息散度(spectral information divergence，SID)和均方根误差(root meansquare error， RMSE)来度量与真实端元和丰度的偏差。

对于端元光谱提取的准确性，本发明采用SID来度量第p个端元光谱W_p(W即由仿真图像中所有组成物质的光谱构成，W_p是端元矩阵W的第p列)和它的估计

之间的相似性，SID对两个光谱之间整体形状的相似性度量很有效。SID越小说明端元提取的结果越准确。

${\mathrm{SID}}_p=D\left(W_p\right|\left|{\hat{W}}_p\right)+D\left({\hat{W}}_p\right|\left|W_p\right)$

其中

是

对于W_o的相关熵，定义如下：

$D\left(W_p\right|\left|{\hat{W}}_p\right)=\underset{l=1}{\overset{L}{\mathrm{\&Sigma;}}}{p}_l\log \left(\frac{p_l}{q_l}\right)$

其中 $p_l=W_{\mathrm{lp}}/\underset{j=1}{\overset{L}{\mathrm{\&Sigma;}}}{W}_{\mathrm{jp}},$ $q_l={\hat{W}}_{\mathrm{lp}}/\underset{j=1}{\overset{L}{\mathrm{\&Sigma;}}}{\hat{W}}_{\mathrm{jp}}.$

度量丰度矩阵准确性，本发明采取均方根误差(root mean square error，RMSE)作为指标，RMSE越小，说明估计的丰度矩阵越接近理论值。其定义如下：

${\mathrm{RMSE}}_p=\sqrt{\frac{1}{K}\underset{k=1}{\overset{K}{\mathrm{\&Sigma;}}}{(H_{\mathrm{pk}}-{\hat{H}}_{\mathrm{pk}})}^2}$

其中，K为仿真图像中像素点个数。

利用所有端元的SID和RMSE的中值来度量算法性能，分别记作SID和RMSE。

(五)具体实施方式：

为了更好地理解本发明的技术方案，以下结合附图对本发明的实施方式作进一步描述：

本发明在MATLAB R2008b语言环境下实现。计算机读取高光谱遥感图像数据后，获取的是数据立方体，其计算机配置采用：Intel(R)Core(TM)2Duo CPU E73002.66GHz。

本发明一种基于随机投影的高光谱图像稀疏解混方法，其流程图见图1所示，该解混方法包括如下步骤：

步骤一：用计算机读取数据。计算机在MATLAB R2008b环境下读取高光谱图像数据，数据来源于成像光谱仪采集到的遥感图像，得到数据立方体，高光谱图像数据应去除被水汽吸收的波段和信噪比较低的波段。将高光谱图像数据逐像素点排列，得到原始的高光谱图像矩阵V。假设高光谱图像有m个波段，共有n个像素点，V＝[v₁，v₂，...，v_n]，v_i是高光谱图像第i个像素点的光谱列向量，是一个m维向量。读取现有的高光谱数据库数据，选择光谱库中大量的纯物质光谱数据构建光谱库矩阵S。假设纯物质光谱的数量为q，S＝[s₁，s₂，...，s_q]，s_i是光谱库中第i个纯物质的光谱列向量，同样是一个m维向量。

步骤二：高光谱数据随机投影。用计算机读取数据后，对高光谱图像数据和高光谱库数据进行随机投影，就是用一个所包含的元素符合正态分布的矩阵R_d与矩阵V，S分别相乘。假设矩阵R_d的大小为d×m，d＜m，得到降维后的数据R_dV和R_dS，其矩阵大小分别为d×n，d×q，代替原始的高光谱图像数据V和高光谱数据库S参与下述步骤三的计算，其中d是随机投影后光谱波段的个数。

步骤三：高光谱图像稀疏解混。数据进行随机投影后，接下来就要对降维后的高光谱图像数据进行光谱解混。假设高光谱数据库中含有足够丰富的纯物质光谱特征，那么高光谱图像中含有的端元只占极少一部分，也就是说高光谱图像中每个像素点的光谱曲线是由高光谱数据中少数的纯物质的光谱曲线线性组合构成，这体现了端元的稀疏性表达。端元在高光谱数据库中是稀疏的，他们所对应的丰度也具有稀疏性，即丰度矩阵是稀疏的，以丰度的稀疏性作为正则化项，符合实际的物理意义。在这里，采用线性混合模型，即将一条检测得到的光谱曲线，分解为光谱库中纯物质光谱线性组合的形式，其系数即为相应丰度。假设丰度矩阵为H，其大小为q×n，其所有元素满足非负性。线性混合模型如下表示：

V＝SH

对原始高光谱图像稀疏解混的目标函数为：

$\underset{H}{\min }{\left|\right|H\left|\right|}_1$     s.t.V＝SH，H≥0

||H||₁表示H包含的所有元素绝对值的和。因为丰度具有非负性，||H||₁越小，丰度矩阵H含有的非零元素值将有可能越少，其表现出的稀疏性就较强。将上式转变为无约束优化问题：

$\underset{H}{\min }\frac{1}{2}{\left|\right|V-\mathrm{SH}\left|\right|}_2^2+\mathrm{\&lambda;}{\left|\right|H\left|\right|}_1$

在上式目标函数中，l₂项是拟合项，体现的是解混的结果与原高光谱图像的拟合程度；l₁即稀疏性正则化项；λ是正则化参数，调节着拟合项与正则化项的权重，其大小根据实际的高光谱图像数据进行设置。因为V，S数据量很大，在迭代求解的过程中占用内存较大，并且计算时间较长。本发明采用随机投影算法，能够有效保持上式中的l₂距离项，实现数据量压缩。随机投影后的目标函数为：

$\underset{H}{\min }\frac{1}{2}{\left|\right|R_{d}V-R_d\mathrm{SH}\left|\right|}_2^2+\mathrm{\&lambda;}{\left|\right|H\left|\right|}_1$

针对上述目标函数，采用分裂Bregman算法求解上述函数极值，从而得到该优化模型下的最优解——最优丰度矩阵H。

引入辅助变量D，并且D＝H，上式目标函数转化为：

$\begin{array}{cc}\underset{H,d}{\min }\frac{1}{2}{\left|\right|R_{d}V-R_d\mathrm{SH}\left|\right|}_2^2+\mathrm{\&lambda;}{\left|\right|D\left|\right|}_1& s.t.D=H\end{array}$

无约束形式为：

$\underset{H,d}{\min }{\mathrm{\&lambda;}{\left|\right|D\left|\right|}_1+\frac{\mathrm{\&mu;}}{2}\left|\right|R_{d}V-R_d\mathrm{SH}\left|\right|}_2^2+\frac{1}{2}{\left|\right|D-H\left|\right|}_2^2$

其中，μ是正则化参数。引入误差反馈变量B，上述优化问题转变为迭代求解过程，如下：

$\left\{\begin{array}{c}{H}^{k+1}=\arg \min \{\frac{\mathrm{\&mu;}}{2}{\left|\right|R_{d}V-R_d\mathrm{SH}\left|\right|}_2^2+\frac{1}{2}{\left|\right|D^k-H-B^k\left|\right|}_2^2\}\\ D^{k+1}=\arg \min \{\mathrm{\&lambda;}{\left|\right|D\left|\right|}_1+\frac{1}{2}{\left|\right|D-H^{k+1}-B^k\left|\right|}_2^2\}\end{array}\right.$

具体求解步骤如下：

(1)初始化丰度矩阵H以及辅助变量矩阵B，D，V_d，H⁰＝B⁰＝D⁰＝0；V_d ⁰＝R_dV

(2)迭代变量H，D，B，V_d，交替迭代过程如下：

$\left\{\begin{array}{c}{H}^{k+1}={[\mathrm{\&mu;}{\left(R_{d}S\right)}^T\left(R_{d}S\right)+I]}^{-1}[\mathrm{\&mu;}{\left(R_{d}S\right)}^T{V_d}^k-B^k+D^k]\\ D^{k+1}=\max (H^{k+1}+B^k-\mathrm{\&lambda;},0)\\ B^{k+1}=B^k+H^{k+1}-D^{k+1}\\ {V_d}^{k+1}={V_d}^k+R_{d}V-R_d{\mathrm{SH}}^{k+1}\end{array}\right.$

这里k为迭代次数，V_d ^k+1的迭代过程可以加速收敛。

(3)停止迭代条件。当丰度矩阵H相邻两次迭代的相对误差不大时，迭代停止。迭代停止条件为：

$\frac{{\left|\right|H^{k+1}-H^k\left|\right|}_2^2}{{\left|\right|H^k\left|\right|}_2^2}<\mathrm{tol}$

本发明中，取tol＝10^-6，若迭代步数少于10次就能收敛，使用迭代次数控制迭代过程，设迭代次数为1000步。

步骤四：获取丰度图及真正的端元。在获得丰度矩阵后，设置合适的阈值，丰度矩阵中小于阈值的元素设为零，不小于阈值的元素不做处理。找出丰度矩阵中含有非零元素的行标号，丰度矩阵的行标号与光谱库的列标号对应。取出光谱矩阵中列标号所对应的列，即为真正的端元光谱；取出丰度矩阵中行标号所对应的行，即为端元所对应的丰度图。在丰度图中颜色较亮的区域表示端元所占的比例较大，较暗的区域对应的比例小，颜色为黑色的区域表示不含有该端元。

有益效果：

实验结果：为了本发明方法的有效性，我们采用仿真的纯像元和混合像元的哈勃望远镜高光谱图像进行试验，通过图3(a)-(h)、图4(a)-(h)与图2(a)-(h)的对比，可以看出，本发明能够准确选择出纯像元和混合像元高光谱图像中的端元，在表1中的SID＝0，也能说明这一点。列表1中，RMSE能够说明本发明求丰度解的精确性，所求解的丰度值与理论丰度的偏差在10^-3数量级，其精度比较高。

表1本发明使用仿真高光谱图像进行性能检测

从实验结果可以看出，本发明对高光谱图像进行光谱解混能够取得较好的效果，能够准确提取端元，并得到精度较高的丰度图，实现了对高光谱图像定量分析，它在高光谱遥感图像分析技术领域里具有实用价值和广阔的应用前景。

Claims (1)

1.一种基于随机投影的高光谱图像稀疏解混方法，其特征在于：该方法具体步骤如下：

步骤一：用计算机读取数据：计算机在MATLAB R2008b环境下读取高光谱图像数据，数据来源于成像光谱仪采集到的遥感图像，得到数据立方体，高光谱图像数据去除被水汽吸收的波段和信噪比较低的波段，将高光谱图像数据逐像素点排列，得到原始的高光谱图像矩阵V；设高光谱图像有m个波段，共有n个像素点，V=[v₁,v₂,...,v_n]，v_i是高光谱图像第i个像素点的光谱列向量，是一个m维向量；读取现有的高光谱数据库数据，选择光谱库中纯物质光谱数据构建光谱库矩阵S；设纯物质光谱的数量为q，S=[s₁,s₂,...,s_q]，s_i是光谱库中第i个纯物质的光谱列向量，同样是一个m维向量；

步骤二：高光谱数据随机投影：用计算机读取数据后，对高光谱图像数据和高光谱库数据进行随机投影，就是用一个所包含的元素符合正态分布的矩阵R_d与矩阵V,S分别相乘；设矩阵R_d的大小为d×m,d<m，得到降维后的数据R_dV和R_dS，其矩阵大小分别为d×n,d×q，代替原始的高光谱图像数据V和高光谱数据库S参与下述步骤三的计算，其中d是随机投影后光谱波段的个数；

步骤三：高光谱图像稀疏解混：数据进行随机投影后，就对降维后的高光谱图像数据进行光谱解混；设高光谱数据库中含有足够丰富的纯物质光谱特征，那么高光谱图像中含有的端元只占极少一部分，也就是说高光谱图像中每个像素点的光谱曲线是由高光谱数据中少数的纯物质的光谱曲线线性组合构成，这体现了端元的稀疏性表达；端元在高光谱数据库中是稀疏的，他们所对应的丰度也具有稀疏性，即丰度矩阵是稀疏的，以丰度的稀疏性作为正则化项，符合实际的物理意义；在这里，采用线性混合模型，即将一条检测得到的光谱曲线，分解为光谱库中纯物质光谱线性组合的形式，其系数即为相应丰度；设丰度矩阵为H，其大小为q×n，其所有元素满足非负性，线性混合模型如下表示：

V=SH

对原始高光谱图像稀疏解混的目标函数为：

$\begin{array}{cc}\underset{H}{\min {\left|\right|H\left|\right|}_1}& s.t.V=\mathrm{SH},H\&GreaterEqual;0\end{array}$

||H||₁表示H包含的所有元素绝对值的和，因为丰度具有非负性，||H||₁越小，丰度矩阵H含有的非零元素值将越少，其表现出的稀疏性就较强；将上式转变为无约束优化问题：

$\underset{H}{\min }\frac{1}{2}{\left|\right|V-\mathrm{SH}\left|\right|}_2^2+\mathrm{\&lambda;}{\left|\right|H\left|\right|}_1$

在上式目标函数中，||H||₁项是拟合项，体现的是解混的结果与原高光谱图像的拟合程度；||H||₁即稀疏性正则化项；λ是正则化参数，调节着拟合项与正则化项的权重，其大小根据实际高光谱数据设置；因为V,S数据量很大，在迭代求解的过程中占用内存较大，并且计算时间较长；采用随机投影算法，能够有效保持上式中的

距离项，实现数据量压缩，随机投影后的目标函数为:

$\underset{H}{\min }\frac{1}{2}{\left|\right|R_{d}V-R_d\mathrm{SH}\left|\right|}_2^2+\mathrm{\&lambda;}{\left|\right|H\left|\right|}_1$

针对上述目标函数，采用分裂Bregman算法求解上述函数极值，从而得到该优化模型下的最优解——最优丰度矩阵H；

引入辅助变量D，并且D=H，上式目标函数转化为：

$\begin{array}{cc}\underset{H,d}{\min }\frac{1}{2}{\left|\right|R_{d}V-R_d\mathrm{SH}\left|\right|}_2^2+\mathrm{\&lambda;}{\left|\right|D\left|\right|}_1& s.t.D=H\end{array}$

转化无约束形式为：

$\underset{H,d}{\min }\mathrm{\&lambda;}{\left|\right|D\left|\right|}_1+\frac{\mathrm{\&mu;}}{2}{\left|\right|R_{d}V-R_d\mathrm{SH}\left|\right|}_2^2+\frac{1}{2}{\left|\right|D-H\left|\right|}_2^2$

其中μ是正则化参数，其大小根据实际高光谱数据设置，引入误差反馈变量B，上述优化问题转变为迭代求解过程，如下：

$\left\{\begin{array}{c}{H}^{k+1}=\arg \min \{\frac{\mathrm{\&mu;}}{2}{\left|\right|R_{d}V-R_d\mathrm{SH}\left|\right|}_2^2+\frac{1}{2}{\left|\right|D^k-H{B}^k\left|\right|}_2^2\}\\ D^{k+1}=\arg \min \{\mathrm{\&lambda;}{\left|\right|D\left|\right|}_1+\frac{1}{2}{\left|\right|{D-H}^{k+1}-B^k\left|\right|}_2^2\}\end{array}\right.$

具体求解步骤如下:

(1)初始化丰度矩阵H以及辅助变量矩阵B,D,V_d，H⁰=B⁰=D⁰=0;V_d ⁰=R_dV

(2)迭代变量H,D,B,V_d,交替迭代过程如下：

$\left\{\begin{array}{c}{H}^{k+1}={[\mathrm{\&mu;}{\left(R_{d}S\right)}^T\left(R_{d}S\right)+I]}^{-1}[\mathrm{\&mu;}{\left(R_{d}S\right)}^T{V_d}^k-B^k+D^k]\\ D^{k+1}=\max (H^{k+1}+B^k-\mathrm{\&lambda;},0)\\ B^{k+1}=B^k+H^{k+1}-D^{k+1}\\ {V_d}^{k+1}={V_d}^k+R_{d}V-R_d{\mathrm{SH}}^{k+1}\end{array}\right.$

这里k为迭代次数，

的迭代过程可以加速收敛；I表示单位矩阵；

(3)停止迭代条件：当丰度矩阵H相邻两次迭代的相对误差不大时，迭代停止；迭代停止条件为:

$\frac{{{\left|\right|H}^{k+1}{-H}^k\left|\right|}_2^2}{{{\left|\right|H}^k\left|\right|}_2^2}<\mathrm{tol}$

这里，取tol=10^-6，若迭代步数少于10次就能收敛，使用迭代次数控制迭代过程，设迭代次数为1000步；

步骤四：获取丰度图及真正的端元：在获得丰度矩阵后，设置合适的阈值，丰度矩阵中小于阈值的元素设为零，不小于阈值的元素不做处理；找出丰度矩阵中含有非零元素的行标号，丰度矩阵的行标号与光谱库的列标号对应；取出光谱矩阵中列标号所对应的列，即为真正的端元光谱；取出丰度矩阵中行标号所对应的行，即为端元所对应的丰度图；在丰度图中颜色较亮的区域表示端元所占的比例较大，较暗的区域对应的比例小，颜色为黑色的区域表示不含有该端元。

Images