A Robust Minimum Volume Based Algorithm with Automatically Estimating Regularization Parameters for Hyperspectral Unmixing
-
摘要: 对高光谱图像解混的目的在于从低空间分辨率的高光谱图像中找到端元与对应的丰度.本文根据解混算法中的最小体积准则，提出了一种自适应鲁棒最小体积高光谱解混算法（Robust minimum volume based algorithm with automatically estimating regularization parameters for hyperspectral unmixing，RMVHU）.本算法通过引入负数惩罚正则项，替换了同类算法中的丰度非负性约束（Non-negativity constraint，ANC），使算法对图像中的噪声与异常值具有更强的鲁棒性；采用循环最小化方法，将非凸优化问题分解为凸优化子问题，然后应用交替方向乘子法解决随着像素点个数增大带来的求解困难问题；对于正则项系数，本算法提出了一种自适应调整策略，提高了算法的收敛性，并且通过定性分析，说明了该调整方法的合理性.将算法应用于合成数据与实际数据，实验结果表明，与同类算法相比，本文提出的算法能够取得更为优秀的效果.Abstract: Hyperspectral unmixing aims at finding hidden endmembers and their corresponding abundances from hyperspectral images with low spatial resolution. Based on the well-known minimum volume (MV) rule in geometrical based approaches, a robust minimum volume based algorithm with automatically estimating regularization parameters for hyperspectral unmixing (RMVHU) is proposed. In this algorithm, the ANC constraint is replaced with a negative number punished regularizer which may lead to a more robust result to outliers and noise. A cyclic minimization algorithm is used to split the nonconvex RMVHU problem into convex subproblems, and ADMM is referred to sovle the large scale optimization problem with the increasing number of pixels in the image. To improve the convergence of the algorithm, a strategy to estimate the regularization parameters of the regularizer automatically is proposed. Compared with some existing geometrical based methods, experimental results show the superiority of the RMVHU algorithm on both synthetic datasets and real datasets.
-
随着遥感传感器的飞速发展, 高光谱图像由于包含丰富的光谱间信息与空间信息, 在分类、探测等方面得到了更加广泛的应用.与普通图像不同, 在高光谱图像中, 每一个像素点(像元)包含了成百上千个波段的反射率信息[1-2].高光谱的高谱间分辨率, 为像素级乃至亚像素级的图像分类与探测提供了条件, 例如文献[3]就利用高光谱的空谱特性进行了有效的异常探测.在光谱图像中, 若某一像元中只存在一种物质, 我们将这样的像元称为纯像元[4].但是高光谱图像的空间分辨率较低, 图像中存在着大量混合像元[4].混合像元广泛存在于高光谱图像中, 是影响遥感分类精度和目标探测效果的重要因素之一.为了进一步利用高光谱数据, 我们需要通过分解混合像元得到图像中一系列基本物质(这样的基本物质被称作端元)的光谱信息, 同时还要求取端元在混合像元中的占比(丰度), 而上述通过原始高光谱图像得到端元与丰度的过程就是高光谱图像的解混过程[5].
近年来, 学者们提出了一系列基于线性混合模型的解混算法.在线性混合模型中, 像元可以由一系列不相关的端元线性表示, 且丰度满足"非负性" (ANC)与"和为1" (ASC)两个约束条件[6].在线性模型的基础上, 基于凸面几何学的解混算法被广泛研究.该类算法的主要思想是:根据线性模型的全约束条件, 所有数据点均被包含在某一类单形体中, 在这些单形体中, 由端元作为顶点构成的单形体体积是最小的.同时也有较多学者研究基于非负矩阵的解混算法[7-10]与稀疏解混算法[11-15].相较于其他算法, 凸面几何学类解混算法运算速度快, 解混精度较高, 因此本文将继续深入研究凸面几何学的算法.
在基于凸面几何学的算法中, 有一类基于纯像元假设的算法, 如VCA[16]、PPI[17]、NFINDR[18]等.虽然这些算法在图像满足纯像元假设时[5]有不俗的表现, 但是若图像中存在大量高度混合的数据, 解混精度将大大下降.为解决上述问题, 学者们提出了一类基于最小体积变换的算法, 该类算法能够在不满足纯像元假设的情况下取得较好的端元提取效果, 如MVES[19]、MVSA[20]、SISAL[21]等.但是由于MVES、MVSA这类算法须严格满足丰度非负约束, 观测数据的异常值将会对解混精度造成很大的影响; 而SISAL虽然通过惩罚项放宽了约束条件, 但是在求解优化问题时, SISAL存在两个问题: 1)将目标函数做了二次近似, 将导致近似的目标函数由于高阶项的缺失, 在某些情况下会大尺度地偏离原目标函数; 2)二次项的Hesse矩阵采用对角矩阵近似而非目标函数的Hesse矩阵.上述问题的存在将给解混精度的提高带来一定的限制.
为了改善上述算法出现的问题, 本文提出了一种自适应鲁棒最小体积解混算法, 该算法有如下优点:
1) 将约束条件放宽, 通过引入负数惩罚正则项替代丰度非负性约束, 获得了更强的抗干扰性能.
2) 求解的目标函数根据行列式的形式展开, 变换之后的目标函数与原函数完全等价, 不再是二次近似, 解决了二次逼近带来的误差.
3) 为了能够解决惩罚正则项带来的大规模凸优化求解困难问题, 应用交替方向乘子法[22] (ADMM), 获得了理想的计算精度与速度.
4) 创造性地提出了一种自适应的正则系数调整方法, 提高了算法的收敛性与稳定性, 并通过定性的分析说明了此自适应参数调整方法的合理性.
5) 对于ADMM中的惩罚系数, 应用文献[22]中自适应调节的方法, 进一步加快了算法的收敛速度.
文章结构主体如下所述.在第1节中详细介绍线性混合模型, 凸面几何学的数学基础以及最小体积解混模型; 在第2节中讲述自适应鲁棒最小体积高光谱解混算法(RMVHU)模型, 参数的自适应调整方法, RMVHU算法步骤以及算法复杂度与收敛性; 在第3节中讲述实验用的数据来源, 并将数据应用于算法的实验与分析; 在第4节中, 总结实验的分析结果, 得出结论.
1. 几何解混相关工作
1.1 高光谱线性混合模型
通常情况下, 在线性混合模型中, 高光谱图像中的每个像元都可被近似认为是图像中各个端元的线性组合:
$ \boldsymbol{ y}=\sum\limits_{i=1}^{p}\boldsymbol{a}_{i}s_{i}+\boldsymbol{e}=A \boldsymbol{ s}+\boldsymbol{ e} $
(1) $ s_{i}\geq{0}, \forall{i}\leq{p} $
(2) $ \sum\limits_{i=1}^{p}s_{i}=1 $
(3) 式中, $p$为端元个数, $\boldsymbol{y}\in {\bf R}^{{L}}$为任意像元的维光谱向量(L为图像波段数), $A=\big[{matrix} \boldsymbol{a}_{1}&\boldsymbol{a}_{2}& \cdots & \boldsymbol{a}_{p} {matrix}\big]$为大小是$L\times p$的端元矩阵, $\boldsymbol{a}_{i}$为端元的向量表示. $\boldsymbol{s}=\big[{matrix} {s}_{1}&{s}_{2}& \cdots & {s}_{p} {matrix}\big]^{\rm T}$为丰度的向量表示, ${s}_{i}$表示像元$\boldsymbol{y}$中端元$\boldsymbol{a}_{i}$所占的比例, $\boldsymbol{e}$为误差项.
线性混合模型一般可分为3种情形:式(1)为无约束的线性混合模型; 加上约束条件(2)则为非负约束混合模型; 再加上约束条件(3)则为全约束混合模型.线性解混就是提取端元, 同时求出各个端元在像元中所占的比例, 得到端元丰度的过程.
在研究几何学高光谱解混算法之前, 本文将依次介绍相关数学符号的定义、求解凸优化问题的著名算法框架——ADMM, 以及仿射集与凸包的数学概念.
本文应用到的符合及其意义如表 1所示:
表 1 数学符号及其意义Table 1 Mathematical notations and their meaning符号 符号意义 ${\bf R}$ 实数 ${\bf R}^{{L}}$ ${L}$维向量 ${\bf R}^{{L}\times {\rm N}}$ ${L}\times {\rm N}$维矩阵 ${\bf R}_{+}$ 非负实数 ${\bf R}_+^{L}$ 非负${L}$维向量 ${\bf R}_{+}^{{L}\times {N}}$ 非负${L}\times {N}$维矩阵 $1_{N}$ ${N}\times 1$维全1列向量 $I_{N}$ ${N}\times {N}$维单位矩阵 $X^†$ $X$的广义逆矩阵 $\boldsymbol{ e}_i$ 第$i$个元素为1, 其他元素为0的单位列向量 $\boldsymbol{x}\succeq {a}, {\boldsymbol{x}}\in {\bf R}^{{L}}$ $L$维列向量$\boldsymbol{x}$中所有元素均不小于$a$ $l_{\bf R_{+}}(\boldsymbol{x})$ $l_{\bf R_{+}}(\boldsymbol{x})={cases} 0, &\boldsymbol{x}\succeq 0 \\ + \infty, & \mbox{其他} {cases}$ ${\rm sgn} (y)$ ${\rm sgn}(y)={cases} 1, &y > 0 \\ 0, & y=0\\-1, &y < 0 {cases}$ 1.1.1 交替方向乘子法(ADMM)
考虑如下线性约束优化问题：
$ \left\{ {\begin{array}{*{20}{l}} {\mathop {\min }\limits_{\boldsymbol{x}, \boldsymbol{z}} f({\rm}\boldsymbol{x}) + g({\rm}\boldsymbol{z})}\\ {{\rm{s}}.{\rm{t}}.A{\rm}\boldsymbol{x} + B{\rm}\boldsymbol{z} = {\rm}\boldsymbol{c}} \end{array}} \right. $
(4) 式中, $\boldsymbol{x}\in {\bf R}^{{n}}$, $\boldsymbol{z}\in {\bf R}^{{m}}$, 矩阵${A}\in {\bf R}^{{l} \times {n}}$, ${B}\in {\bf R}^{{l} \times {m}}$, $\boldsymbol{c}\in {\bf R}^{{p}}$.根据文献[22], 若函数$f(\boldsymbol{ x})$与$g(\boldsymbol{ z})$是凸函数且问题(4)有最优解, 那么在式(5)$\sim$(7)的迭代步骤下, 序列$\big\{\boldsymbol{x}^k \big\}$则将收敛到最优解上, 否则序列$\big\{ \boldsymbol{z}^k\big\}$或者$\big\{ \boldsymbol{d}^k \big\}$将发散.
$ \begin{array}{l} {{\boldsymbol{x}}^{k + 1}} = \arg \;\mathop {\min }\limits_\boldsymbol{ x} f({\rm{ }}\boldsymbol{x}) + \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\\ \frac{\mu }{2}\left\| {A{\rm{ }}\boldsymbol{x} + B{\rm{ }}{\boldsymbol{z}^k} - {\rm{ }}\boldsymbol{c} - {\rm{ }}{\boldsymbol{d}^k}} \right\|_2^2 \end{array} $
(5) $ \begin{align} \,\boldsymbol{ z}^{k + 1} = \arg\;\mathop {\min }\limits_\boldsymbol{ z} f({\rm{ }}\boldsymbol{x}) +\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\nonumber\\ \frac{\mu}{2}{\left\| {A \boldsymbol{ x}^{k+1} + B{\boldsymbol{ z}} - \boldsymbol{ c} - {\boldsymbol{ d}^k}} \right\|}_2^2 \end{align} $
(6) $ \begin{align} \boldsymbol{ d}^{k + 1} = \boldsymbol{ d}^{k}- \left( { A\boldsymbol{ x}^{k+1} + B\boldsymbol{ z}^{k+1} - \boldsymbol{ c} } \right) \end{align} $
(7) 式中, $\mu$为与收敛速度相关的常值, $\boldsymbol{d}$为尺度对偶变量(Scaled dual variable)[22].
采用ADMM的框架, 能够在一步步的迭代过程中逼近凸问题的最优解.
1.1.2 仿射集与凸包
有向量集$\big\{ \boldsymbol{ a}_1, \boldsymbol{ a}_2, \cdots, \boldsymbol{ a}_{{p}} \big\} \subset {{\bf R}^{\ L}}$, 则向量集$\big\{ \boldsymbol{ a}_1, \boldsymbol{ a}_2, \cdots, \boldsymbol{ a}_{p} \big\}$的仿射包定义为:
$ \begin{align}%tag(8) {\rm aff}\{ {\boldsymbol{ a}}_1, {\boldsymbol{ a}}_2, \cdots, {\boldsymbol{ a}}_{{p}} \big\} = \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\nonumber\\ \Big\{ \boldsymbol{ x} = \sum\limits_{i = 1}^{p} {{{\theta}_i}{\boldsymbol{ a}_i}} \Big| {1_{p}^{\rm T}\boldsymbol{ \theta} = 1, \boldsymbol{ \theta} \in {\bf R}^{p} } \Big. \Big\} \end{align} $
(8) 仿射包是一个仿射集, 所以也能被表示为:
$ \begin{align}%tag(9) &{\rm aff}\left\{ \boldsymbol{ a}_1, \cdots, \boldsymbol{ a}_{{p}} \right\} =\left\{ \boldsymbol{ x}= C\boldsymbol{ \alpha} + \boldsymbol{ d} \left| \boldsymbol{ \alpha}\in{\bf R}^{{k}} \right. \right\}:=\nonumber\\ &\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; \Lambda(C, \boldsymbol{ d}) \end{align} $
(9) 式中, $\Lambda(C, \boldsymbol{ d})$表示由参数$(C, \boldsymbol{ d})$构成的仿射集, $C \in {{\bf R}^{{{L}} \times {{k}}}}$, $\boldsymbol{ d} \in {{\bf R}^{{L}} }$, 参数$C$需满足秩条件——rank$(C)={ k}$. ${\ k}$为仿射包的仿射维数, 其必须满足${\ k}\leq{ p}-1$的约束条件.若向量集$\big\{ \boldsymbol{ a}_1, \cdots, \boldsymbol{ a}_{p} \big\}$仿射无关(即${\ p}-1$个向量$\big\{ \boldsymbol{ a}_1-\boldsymbol{ a}_{\ p}, \cdots, \boldsymbol{ a}_{{\ p}-1}-\boldsymbol{ a}_{p} \big\}$是线性无关的), 则${\ k}={\ p}-1$.
在文献[23]中, 通过求解式(10)描述的问题, 根据向量集$\big\{ \boldsymbol{ a}_1, \cdots, \boldsymbol{ a}_{p} \big\}$与$\ k$, 得到一组参数$(C, \boldsymbol{ d})$来表示仿射集.
$ \begin{align}%tag(10) (C, \boldsymbol{ d}) = \arg\;\underset{{C, \boldsymbol{ d}}\atop{C^{\rm T}C={I_{ k} } } } {\min}\sum\limits_{i=1}^{p}{e_{\Lambda(C, \, \boldsymbol{ d})}(\ \boldsymbol{a}_i)} \end{align} $
(10) 式中, 约束$C^{\rm T}C={I_{\ k} } $是为了满足rank$(C)={ k}$的秩约束条件, $e_{\Lambda(C, \, \boldsymbol{ d})}(\ \boldsymbol{ a}_i)$是$\boldsymbol{ a}_i$投影到仿射集$\Lambda(C, \, \boldsymbol{ d})$上的误差, $e_{\Lambda(C, \, \boldsymbol{ d})}(\ \boldsymbol{ a}_i)$定义为:
$ \begin{align}%tag(11) e_{\Lambda(C, \, \boldsymbol{ d})}(\ a_i) =\underset{\boldsymbol{ a}\in{\Lambda(C, \, \boldsymbol{ d})}}{\min}{\left\| \boldsymbol{ a}_i - \boldsymbol{ a} \right\|^2} \end{align} $
(11) 问题(10)有如下闭合解的形式[23]:
$ \begin{align} \begin{cases} %\tag(12) \begin{array}{lll} C& = &[{\begin{array}{ccc} {\boldsymbol{ q}_1(U^{\rm T}U)}&{\cdots }&{\boldsymbol{ q}_{\ k}(U^{\rm T}U)} \end{array}}]\\ \boldsymbol{ d}&=&\frac{1}{\ p} \sum\limits_{i = 1}^{p}\boldsymbol{ a}_i\\ U&= &[{\begin{array}{ccc} {\boldsymbol{ a}_1-\boldsymbol{ d}}&{\cdots }&{\boldsymbol{ a}_{p}-\boldsymbol{ d}} \end{array}}] \end{array} \end{cases} \end{align} $
(12) 式中, ${\boldsymbol{ q}_i(U^{\rm T}U)}$代表$U^{\rm T}U$的第$i$个主成分.
下面在仿射集的基础上介绍凸集的概念.给定向量集$\big\{ \boldsymbol{ a}_1, \boldsymbol{ a}_2, \cdots, \boldsymbol{ a}_{{p}} \big\} \subset {{\bf R}^{ L}}$, 则由其表示的凸包的定义为:
$ \begin{align} \,{\rm conv}\big\{ \boldsymbol{ a}_1, \boldsymbol{ a}_2, \cdots, \boldsymbol{ a}_{p} \big\} = \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\nonumber\\ \Big\{ \boldsymbol{ x} = \sum\limits_{i = 1}^{p} {{{\theta}_i}{\boldsymbol{ a}_i}} \Big| {1_{p}^{\rm T}\boldsymbol{ \theta} = 1, \boldsymbol{ \theta} \in {\bf R}^{p}_{+} } \Big. \Big\} \end{align} $
(13) 凸包中的点$\boldsymbol{x}$若不能被表示为严格凸组合的形式, 即:若${1_{p}^{\rm T}\boldsymbol{ \theta} = 1, \boldsymbol{ \theta}\in {\bf R}^{p}_{+}}$且$\boldsymbol{ \theta}\neq {\pmb e}_i, \forall i=1, \cdots, {\ p}$, $\boldsymbol{ x} \neq {\sum_{i=1}^{p}}{\theta}_i\boldsymbol{ a}_i$, 称$\boldsymbol{x}$为凸包的顶点.若${\ L}={\ p}-1$且$\big\{ \boldsymbol{ a}_1, \boldsymbol{ a}_2, \cdots, \boldsymbol{ a}_{{p}} \big\}$是仿射无关的, 则称凸包$\rm{conv}\big\{ \boldsymbol{ a}_1, \boldsymbol{ a}_2, \cdots, \boldsymbol{ a}_{{p}} \big\}$为单形体, 向量集$\big\{ \boldsymbol{ a}_1, \boldsymbol{ a}_2, \cdots, \boldsymbol{ a}_{{p}} \big\}$为单形体顶点的集合.
仿射集、凸包的概念如下图所示:
从图 1中可以看到, 由三个顶点描述的仿射集是一个平面, 而凸包则是在仿射集平面上, 由顶点$\boldsymbol{ a}_1, \boldsymbol{ a}_2, \boldsymbol{ a}_3$组成的三角形内的区域, 这个三角形也被称为二维单形体.
图 1 顶点个数为3时, 凸包与仿射集的概念说明Fig. 1 Illustration of convex hull and affine hull when the number of vertices is three1.2 最小体积解混模型
一般来说, 式(1)中的端元矩阵是列满秩的, 即各个端元%向量线性无关.又由于式(3) "和为1"的假设, 文献[19]指出, 由观测数据集$\big\{ \boldsymbol{ y}[1], \boldsymbol{ y}[2], \cdots, \boldsymbol{ y}[{{N}}] \big\}$组成的仿射包与由端元$\big\{ \boldsymbol{ a}_1, \boldsymbol{ a}_2, \cdots, \boldsymbol{ a}_{{p}} \big\}$组成的仿射包是同一个.因此${\rm aff} \big\{ \boldsymbol{ y}[1], \boldsymbol{ y}[2], \cdots, \boldsymbol{ y}[{N}] \big\}$与${\rm aff} \big\{ \boldsymbol{ a}_1, \boldsymbol{ a}_2, \cdots, \boldsymbol{ a}_{p} \big\}$都能由同一组参数$(C, \boldsymbol{ d})$来表示, 其中:
$ \begin{align} \begin{cases} \begin{array}{l} C= [{\begin{array}{ccc} {\boldsymbol{ q}_1(U^{\rm T}U)}&{\cdots }&{\boldsymbol{ q}_{{\ p}-1}(U^{\rm T}U)} \end{array}}]\\ \boldsymbol{ d} \; = \frac{1}{\ N} \sum\limits_{n = 1}^{\ N}\boldsymbol{ y}[n]\\ U = [{\begin{array}{ccc} {\boldsymbol{ y}[1]-\boldsymbol{ d}}&{\cdots }&{\boldsymbol{ y}[{\ N}]-\boldsymbol{ d}} %C& = &[{\begin{array}{ccc} {\boldsymbol{ q}_1(U^{\rm T}U)}&{\cdots }&{\boldsymbol{ %q}_{{\ p}-1}(U^{\rm T}U)} %\end{array}}]\\ %\boldsymbol{ d}&=&\frac{1}{\ N} \sum\limits_{n = 1}^{\ N}\boldsymbol{ y}[n]\\ %U&= &[{\begin{array}{ccc} {\boldsymbol{ y}[1]-\boldsymbol{ d}}&{\cdots }&{\boldsymbol{ %y}[{\ N}]-\boldsymbol{ d}} \end{array}}] \end{array} \end{cases} \end{align} $
(14) 由于观测点$\boldsymbol{ y}[n]$是在仿射集${\Lambda(C, \boldsymbol{ d})}$中, 所以我们能够写成如下形式:
$ \begin{align} \boldsymbol{ y}[n]=C\tilde{\boldsymbol{ y}}[n]+\boldsymbol{ d} \end{align} $
(15) $\tilde{y}[n]$能够被表示为:
$ \begin{align} \tilde{\boldsymbol{ y}}[n]=C^{†}(\boldsymbol{ y}[n]-\boldsymbol{ d})\in{\bf R}^{{ p}-1} \end{align} $
(16) 式中, $\tilde{\boldsymbol{ y}}[n]$也是高光谱原始数据$\boldsymbol{ y}[n]$的低维表示.将式(3)代入式(16), 则有:
$ \begin{align} %\tag(17) \tilde{\boldsymbol{ y}}[n]=\sum\limits_{i=1}^{p}{s_i}[n]\boldsymbol{ \alpha}_i \end{align} $
(17) 式中, $\boldsymbol{ \alpha}_i = {C^† }({\boldsymbol{ a}_i} -d) \in { {\bf R}^{{\ p} -1}} $.
根据凸包的定义以及式(2)、(3)能够发现, ${\tilde{\boldsymbol{ y}}[n]} \in {\rm conv}\big\{ \boldsymbol{ a}_1, \boldsymbol{ a}_2, \cdots, \boldsymbol{ a}_{p} \big\}$, 且此凸包是单形体; 该单形体与由端元%向量作为顶点构成的凸包具有一一对应的关系.所以如果得到了凸包${\tilde{\boldsymbol{ y}}[n]} \in {\rm conv}\big\{ \boldsymbol{ a}_1, \boldsymbol{ a}_2, \cdots, \boldsymbol{ a}_{p} \big\}$的顶点, 那么通过$\boldsymbol{ \alpha}_i$与$\boldsymbol{ a}_i$的关系, 最终就能够得到端元矩阵.
考虑如下优化问题:
$ \begin{align} \begin{cases} %\tag(18) \underset{ \boldsymbol{ \beta}_1, \cdots, \boldsymbol{ \beta}_{p} }{\min} V(\boldsymbol{ \beta}_1, \cdots, \boldsymbol{ \beta}_{p})\\ {\rm s.t. }{\tilde{\boldsymbol{ y}}[n]} \in {\rm conv}\big\{ \boldsymbol{ \beta}_1, \cdots, \boldsymbol{ \beta}_{p} \big\} \end{cases} \end{align} $
(18) 式中, $V(\boldsymbol{ \beta}_1, \cdots, \boldsymbol{ \beta}_{p})$代表由$\boldsymbol{ \beta}_1, \cdots, \boldsymbol{ \beta}_{p}$作为顶点构成的单形体的体积, $\tilde{\boldsymbol{ y}}[n]$代表光谱数据的低维表示.文献[19]指出, 在纯像元假设下, 问题(18)的最优解如下:
$ \begin{align} %\tag(19) {[\boldsymbol{ \beta}_1^*, \cdots, \boldsymbol{ \beta}_{p}^*]}^{\rm T}= {[\boldsymbol{ \alpha}_1, \cdots, \boldsymbol{ \alpha}_{p}]}^{\rm T} \end{align} $
(19) 同时文献[19]指出, 在纯像元假设不严格成立的情况下, 采用式(19)表示的可行解也接近问题(18)的最优解.根据文献[24], 单形体的体积能够被表示为:
$ \begin{align} \begin{cases} %\tag(20) V(\boldsymbol{ \beta}_1, \cdots, \boldsymbol{ \beta}_{p}) =\frac{1}{(p-1)!} \left| {\det(B)} \right|\\ B=[\boldsymbol{ \beta}_1-\boldsymbol{ \beta}_{p}, \cdots, \boldsymbol{ \beta}_{{ p}-1}-\boldsymbol{ \beta}_{p}] \end{cases} \end{align} $
(20) 根据凸包定义, ${\tilde{\boldsymbol{ y}}[n]} \in {\rm conv}\big\{ \boldsymbol{ \beta}_1, \cdots, \boldsymbol{ \beta}_{p} \big\}$能够被表示为:
$ \begin{align} \begin{cases} %\tag(21) \tilde{\boldsymbol{ y}}[n]=\sum\limits_{j=1}^{p}{s_j}[n]\boldsymbol{ \beta}_j=\boldsymbol{ \beta}_{p}+B\boldsymbol{ s}'[n]\\ {\rm s.t.}s_{p}[n]=1-1_{ {\ p} - 1}^{\rm T}\boldsymbol{ s}'[n]\geq0\\ \boldsymbol{ s}'[n]=[s_{1}[n], \cdots, s_{ {\ p} -1}[n]]^{\rm T}\succeq0 \end{cases} \end{align} $
(21) 将式(20)与式(21)代入式(18), 式(18)描述的优化问题能够被等价表示为:
$ \begin{align}%\tag(22) \begin{cases} \underset{ B, \boldsymbol{ \beta}_{p} }{\min}\left| {\det(B)} \right|\\ {\rm s.t.}\boldsymbol{ s}'[n]\succeq0, 1_{ {\ p} - 1}^{\rm T}\boldsymbol{ s}'[n]\leq1, \\ \tilde{\boldsymbol{ y}}[n]=\boldsymbol{ \beta}_{p}+B\boldsymbol{ s}'[n], \forall n=1, \cdots, {\ N} \end{cases} \end{align} $
(22) 再令$H=B^{-1}, \boldsymbol{ g}=B^{-1}\boldsymbol{ \beta}_{p}$, 则可以使问题(22)的非线性约束转换为线性约束的等价问题:
$ \begin{align}%\tag(23) \begin{cases} \underset{ H, \boldsymbol{ g} }{\max}\left| {\det(H)} \right|\\ {\rm s.t.}1_{ {\ p} - 1}^{\rm T}(H\tilde{\boldsymbol{ y}}[n]-\boldsymbol{ g})\leq1, \\ H\tilde{\boldsymbol{ y}}[n]-\boldsymbol{ g}\succeq0, \forall n=1, \cdots, {\ N} \end{cases} \end{align} $
(23) 式(23)即为最小体积解混算法模型的数学描述.
2. 自适应鲁棒最小体积高光谱解混算法
2.1 RMVHU算法模型
在凸面几何学的最小体积类算法中, 我们通常要求单形体能够包围所有的数据点集, 即必须满足丰度系数非负性条件.然而, 在实际情况下, 在图像中很有可能存在异常点, 或者由于噪声的影响导致像元分布在单行体之外.若在这些情况下仍强行满足丰度非负性条件, 则可能出现为了将这些单形体之外的点包入估计的单形体中, 所估计的单形体顶点与真实的顶点相差甚远的情况.因此, 为了容忍异常点给端元估计带来的影响, 我们需要构造一个负数惩罚正则项替换非负性约束, 从而在尽可能保证单形体体积最小的情况下, 增强算法对异常值、噪声值的鲁棒性.同时通过引入此正则项, 减少在求解非凸优化问题时, 陷入局部最小解的概率.
考虑引入负数惩罚正则项:
$ \begin{align}%\tag(24) \left\| {S} \right\|_{h}=\sum\limits_{i, j}h(s_{ij}), h(s_{ij})=\max({0, -s_{ij}}) \end{align} $
(24) 式中, $s_{ij}$为矩阵$S$的元素.该正则项能够对负系数进行惩罚, 而对于非负系数则没有任何影响.
式(24)中的$h(s_{ij})$, 可以写成如下形式:
$ \begin{align}%\tag(25) h(s_{ij})=\frac{1}{2}\left| {{s_{ij}}} \right| - \frac{1}{2}{s_{ij}} \end{align} $
(25) 将此正则项代替优化问题(23)中的非负约束项, 建立如下优化模型:
$ \begin{align}%tag(26) &\mathop {\min }\limits_ {H, \boldsymbol{ g}}-\left| {\det(H)} \right|+\nonumber\\ &\qquad2\lambda \Big( { \sum\limits_{n=1}^{\ N} }\left\| {H\tilde{\boldsymbol{ y}}[n]-\boldsymbol{ g}} \right\|_{h} + \nonumber\\ &\qquad{ \sum\limits_{n=1}^{\ N} }\left\| {1- 1_{{\ p} - 1}^{\rm T}({H\tilde{\boldsymbol{ y}}[n]-\boldsymbol{ g}}) }\right\|_{h} \Big) \end{align} $
(26) 将式(25)代入式(26), 上述问题等价为如下无约束非凸最优化问题:
$ \begin{align}%tag(27) &\mathop {\min }\limits_ {H, \boldsymbol{ g}}-\left| {\det(H)} \right|+\nonumber\\ &\qquad\lambda \Big( { \sum\limits_{n=1}^{\ N} }\left\| {H\tilde{\boldsymbol{ y}}[n]-\boldsymbol{ g}} \right\|_{1} + \nonumber\\ &\qquad{ \sum\limits_{n=1}^{\ N} }\left| {1- 1_{{\ p} - 1}^{\rm T}({H\tilde{\boldsymbol{ y}}[n]-\boldsymbol{ g}}) }\right| \Big) \end{align} $
(27) 式(27)描述的问题即为RMVHU算法的问题模型.但是式(27)中的目标函数非凸, 很难求解得到全局最优解, 因此, 可以采用循环最小化思想, 逐行更新矩阵变量$H$与$\boldsymbol{g}$的值, 以将此问题转换为带有绝对值的凸优化问题, 随后通过代数余子式展开行列式, 且脱去$\left| {\det(H)} \right|$项的绝对值, 将凸优化问题转换为两个等价的凸优化子问题.
具体的, 考虑更新矩阵变量$H$与$\boldsymbol{g}$的第$i$行值, 得到RMVHU算法的凸优化子问题$p^{*}$与$q^{*}$:
$ \begin{align} \begin{cases} %\tag(28) p^{*}=\mathop {\min }\limits_\boldsymbol{ x}\lambda \left\| {A\boldsymbol{ x}+\boldsymbol{ b}} \right\|_{1} + \boldsymbol{ c}^{\rm T}\boldsymbol{ x}\\ {\rm s.t.} \boldsymbol{ c}^{\rm T} \boldsymbol{ x} < 0 \end{cases} \end{align} $
(28) $ \begin{align} \begin{cases} %\tag(29) q^{*}=\mathop {\min }\limits_\boldsymbol{ x}\lambda \left\| {A\boldsymbol{ x}+\boldsymbol{ b}} \right\|_{1} - \boldsymbol{ c}^{\rm T}\boldsymbol{ x}\\ {\rm s.t.} \boldsymbol{ c}^{\rm T} \boldsymbol{ x} > 0 \end{cases} \end{align} $
(29) 式中, $\boldsymbol{ x}=\big[{{array} {cc} { \boldsymbol{ h}_i }& {g_i} {array}} \big]^{\rm T}$, ${ \boldsymbol{ h}_i }$与$ {g_i}$分别为$H$与$\ g$的第$i$行元素, 常数矩阵$A\in{\bf R}^{2{\ N}\times{p}}$与向量$\boldsymbol{ b}\in{\bf R}^{2{\ N}}, \boldsymbol{ c}\in{\bf R}^{p}$被表示为:
$ \begin{align} %\tag(30) A=\left[{\begin{array}{cccccc} {\tilde{\boldsymbol{ y}}[1]}& \cdots&{\tilde{\boldsymbol{ y}}[\ N]}&{-\tilde{\boldsymbol{ y}}[1]}& \cdots&{-\tilde{\boldsymbol{ y}}[\ N]}\\ -1 &\cdots&-1 &1 & \cdots& 1 \end{array}} \right]^{\rm T} \end{align} $
(30) $ \begin{align} %\tag(31) c=\left[{\begin{array}{cccc} {{H_{i, 1}}}& \cdots&{H_{i, { p}-1}} &0 \end{array}} \right]^{\rm T} \end{align} $
(31) $ \begin{align} \begin{cases} %\tag(32) b=\left[{\begin{array}{cccccc} {0} & \cdots & 0 & b_{1} &\cdots& b_{\ N} \end{array}} \right]^{\rm T}\\ b_{n}=1-\sum\limits_{j = 1, j \ne i}^{{{p}} - 1} \left( { \boldsymbol{ h}_j^{\rm T} {\tilde y}[n] -g_j } \right) \end{cases} \end{align} $
(32) 式(31)中, $H_{i, j}$为$H$中元素$h_{i, j}$的代数余子式.
2.2 RMVHU算法子问题的求解
问题(28)与(29)为带有$l_1$范数的最优化问题.下面介绍一种求解此类的方法.
定理1. 问题(33)的最优解与问题(34)的最优解相同:
$ \begin{align} \begin{cases} %\tag(33) f=\mathop {\min }\limits_\boldsymbol{ x}\lambda \left\| {A\boldsymbol{ x}+\boldsymbol{ b}} \right\|_{1} + \boldsymbol{ c}^{\rm T}\boldsymbol{ x}\\ {\rm s.t.} B \boldsymbol{ x} < \boldsymbol{ q} \end{cases} \end{align} $
(33) $ \begin{align} \begin{cases} %\tag(34) f=\mathop {\min }\limits_\boldsymbol{ x, \ y, \ z}1_{\ N}^{\rm T}{\left( { \boldsymbol{ y} + \boldsymbol{ z} }\right) } + \boldsymbol{ c}^{\rm T}\boldsymbol{ x}\\ {\rm s.t.} B \boldsymbol{ x} < \boldsymbol{ q}\\ \boldsymbol{ y} \succeq 0, \boldsymbol{ z} \succeq 0, \boldsymbol{ y}-\boldsymbol{ z}=A\boldsymbol{ x}+\boldsymbol{ b} \end{cases} \end{align} $
(34) 证明. 假设当前的可行解是$\big(\boldsymbol{ x}_0, \boldsymbol{ y}_0, \boldsymbol{ z}_0 \big)$, $\boldsymbol{ y}_0$与$\boldsymbol{ z}_0$为所有元素大于0的列向量, $f_0=1_{\ N}^{\rm T}{\big({ \boldsymbol{ y}_0 + \boldsymbol{ z}_0 }\big) }+\boldsymbol{ c}^{\rm T}\boldsymbol{ x}_0$.此时, 存在一个$\alpha \succ 0$, 使得$\boldsymbol{ y}_1=\boldsymbol{ y}_0-\boldsymbol{ \alpha}\succ 0, \boldsymbol{ z}_1=\boldsymbol{ z}_0-\boldsymbol{ \alpha} \succ 0 $.
则目标函数值有:
$ \begin{align} %\tag(35) f_1=\, &1_{\ N}^{\rm T}{\big( { \boldsymbol{ y}_1 + \boldsymbol{ z}_1 }\big) }+\boldsymbol{ c}^{\rm T}\boldsymbol{ x}_0=\nonumber\\ &1_{\ N}^{\rm T}{\big( { \boldsymbol{ y}_0 + \boldsymbol{ z}_0 }\big) }+\boldsymbol{ c}^{\rm T}\boldsymbol{ x}_0-2\cdot1_{\ N}^{\rm T}\boldsymbol{ \alpha}<\nonumber\\ &f_0 \end{align} $
(35) 只有当$\boldsymbol{ y}_0$的第$i$个元素$y_0(i)=0$或者$\boldsymbol{ z}_0$的第$i$个元素$z_0(i)=0$时, 能够满足在当前可行解$\boldsymbol{ x}_0$下, 目标函数取得最小值.根据此结论, 结合约束条件$\boldsymbol{ y}-\boldsymbol{ z}=A\boldsymbol{ x}+\boldsymbol{ b}$, 有:
$ \begin{align}%\tag(36) \begin{array}{l} y(i) + z(i) = \left\{ {\begin{array}{*{20}{l}} {\boldsymbol{ a}_i\boldsymbol{ x} + {b_i}, \boldsymbol{ a}_i\boldsymbol{ x} + {b_i} \ge 0}\\ { - \big( \boldsymbol{ a}_i\boldsymbol{ x} + {b_i} \big), \boldsymbol{ a}_i\boldsymbol{ x} + {b_i} < 0} \end{array}} \right.=\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; \big| { \boldsymbol{ a}_i \boldsymbol{ x} + {b_i} } \big| \end{array} \end{align} $
(36) 式中$\boldsymbol{ a}_i$为矩阵$A$的行向量, ${b}_i$为向量$\boldsymbol{ b}$的第$i$个元素.所以, 问题(34)可以等价为问题(33)的形式.
利用定理1, 可以将凸优化子问题(28)与(29)转换为线性规划问题求解.
但是此种解法仅仅适用于图像像素点较少的情况, 原因在于$A\in{\bf R}^{2{\ N}\times{\ p}}$, 转换后新增的变量$\boldsymbol{ y}\in{\bf R}^{2{\ N}}, \boldsymbol{ z}\in{\bf R}^{2{\ N}}$, 变量的个数为${p}+4{\ N}$, 约束条件的个数为$1+6{\ N}$.当像素点个数达到${N}=10\, 000$的时候, 所需解决的就是一个中规模的线性规划问题.而这仅仅是算法中更新矩阵变量$H$与$\boldsymbol{g}$第$i$行元素的子步骤, 若要完成算法的一个完整的迭代过程, 将带来很大的时间开销.因此需要能够快速求解子问题(28)与(29)的新方法.
2010年后, ADMM在求解大型凸优化问题的领域展现了其强大的能力, 也给求解此类凸优化问题带来了契机.
以子问题(28)为例, 给出采用ADMM算法解决此问题的步骤.
问题(28)有如下等价形式:
$ \begin{align} \begin{cases} %\tag(37) \underset{\boldsymbol{ x}, z_1, \boldsymbol{ z_2}}{\min}z_1+\lambda \left\| \boldsymbol{ z_2} \right\|_{1} \\ {\rm s.t.} z_1=\boldsymbol{ c}^{\rm T} \boldsymbol{ x}, z_1 < 0, \boldsymbol{ z_2}={A\boldsymbol{ x}+\boldsymbol{ b}} \end{cases} \end{align} $
(37) 此问题的ADMM形式为:
$ \begin{align}%\tag(38) L(\boldsymbol{ x}, {z_1}, \boldsymbol{ z}_2, d_1, \boldsymbol{ d}_2)=z_1+\lambda\left\| \boldsymbol{ z_2} \right\|_{1} + l_{ {\bf R}_{+} }(-z_1)+\nonumber\\ \frac{\mu}{2}\left\| { \boldsymbol{ c}^{\rm T} \boldsymbol{ x} -z_1 - d_1 } \right\|_{2}^{2}+ \frac{\mu}{2}\left\| { A \boldsymbol{ x} +\boldsymbol{ b} -\boldsymbol{ z}_2 - \boldsymbol{ d}_2 } \right\|_{2}^{2} \end{align} $
(38) 根据ADMM算法的框架, 采用表 2所描述的算法流程, 依次更新$\boldsymbol{ x}, z_1, \boldsymbol{ z}_2, d_1, \boldsymbol{ d}_2$, 求解子问题(28)中的$p^{*}$.
表 2 RMVHU中求解子问题$p^{*}$的算法步骤Table 2 Steps for solving subproblem $p^{*}$ in RMVHU算法名称: 求解问题$p^{*}$的ADMM算法步骤 步骤1: 设初值$k=0, \mu > 0, \boldsymbol{ x}_0, {z}_1^0, {d}_1^0, \boldsymbol{ z}_2^0, \boldsymbol{ d}_2^0 $ 步骤2: 更新$\boldsymbol{x}$: $\boldsymbol{ x}^{k+1} = \arg\; \mathop {\min }\limits_\boldsymbol{ x}L(\boldsymbol{ x}, {z_1^k}, \boldsymbol{ z}_2^k, d_1^k, \boldsymbol{ d}_2^k) $ 步骤3: 更新$z_1$: ${z_1^{k+1} }= \arg\; \underset{z_1}{\min}L(\boldsymbol{ x}^{k+1}, {z_1}, \boldsymbol{ z}_2^k, d_1^k, \boldsymbol{ d}_2^k) $ 步骤4: 更新$\boldsymbol{ z}_2$: $\boldsymbol{ z}_2^{k+1} = \arg\; \underset{\boldsymbol{ z}_2}{\min}L(\boldsymbol{ x}^{k+1}, {z_1^{k+1}}, \boldsymbol{ z}_2, d_1^k, \boldsymbol{ d}_2^k) $ 步骤5: 更新$d_1$: ${d_1^{k+1} }= d_1^k -(\boldsymbol{ c}^{\rm T} \boldsymbol{ x}^{k+1} -z_1^{k+1})$ 步骤6: 更新$\boldsymbol{ d}_2$: 步骤7: 判断是否满足迭代停止条件, 若满足, 则转步骤8, 否则$k=k+1$, 转步骤2. 步骤8: 输出$\boldsymbol{ x}^{k+1}$: 下面将介绍如何求解步骤2至步骤4的优化子问题.
步骤2的子问题为:
$ \begin{align}%\tag(39) \boldsymbol{ x}^{k+1}=\, &\arg\; \mathop {\min }\limits_\boldsymbol{ x} \frac{1}{2}\left\| { \boldsymbol{ c}^{\rm T} \boldsymbol{ x} -z_1^{k} - d_1^{k} }\right\|_{2}^{2}+\nonumber\\ &\frac{1}{2}\left\| { A \boldsymbol{ x} +\boldsymbol{ b} -\boldsymbol{ z}_2^k - \boldsymbol{ d}_2^k } \right\|_{2}^{2} \end{align} $
(39) 式(39)的问题为简单的无约束二次规划问题, 目标函数梯度为0的极值点即为最优解.令问题(39)目标函数的梯度为0, 得到步骤2的解析最优解:
$ \begin{align}%\tag(40) \begin{cases} %\tag(40) \boldsymbol{ x}^{k+1}={\Psi}^{-1}\boldsymbol{ c} ( {z_1^k+d_1^k}) -\\ \qquad\;\;\;\;\;{\Psi}^{-1}A^{\rm T} (\boldsymbol{ b}-\boldsymbol{ z}_2^k-\boldsymbol{ d}_2^k) \\ {\Psi} = \boldsymbol{ c}\boldsymbol{ c}^{\rm T}+{A}^{\rm T}A \end{cases} \end{align} $
(40) 步骤3的子问题为:
$ \begin{align}%\tag(41) {z}_1^{k+1}=\, &\arg\; \underset{{z}_1}{\min} \frac{\mu}{2}\left\| { z_1 - ( \boldsymbol{ c}^{\rm T} \boldsymbol{ x}^{k+1} -d_1^{k}) }\right\|_{2}^{2}+\nonumber\\ &{z}_1+l_{ {\bf R}_{+} }(-z_1) \end{align} $
(41) 该问题是对变量$z_1$解耦的不等式约束最优化问题, 因此我们能够直接给出问题(41)的最优解形式.更方便的是, $z_1$其实就是一个实数.式(41)最优解的解析形式表示为:
$ \begin{align}%\tag(42) {z}_1^{k+1}={\min}\big( {0, \boldsymbol{ c}^{\rm T} \boldsymbol{ x}^{k+1} -d_1^{k}-\frac{1}{\mu}} \big) \end{align} $
(42) 步骤4描述的问题为:
$ \begin{align}%\tag(43) \boldsymbol{ z}_2^{k+1}=\, &\arg\; \underset{\boldsymbol{ z}_2}{\min}\lambda\left\| \boldsymbol{ z_2} \right\|_{1}+\nonumber\\ &\frac{\mu}{2}\left\| { \boldsymbol{ z}_2 - (A \boldsymbol{ x} +\boldsymbol{ b} - \boldsymbol{ d}_2^k ) } \right\|_{2}^{2} \end{align} $
(43) 因此, 问题的解可以采用著名的软阈值[25]给出:
$ \begin{align}%\tag(44) \boldsymbol{ z}_2^{k+1}=\rm{soft}\big( {{A \boldsymbol{ x}^{k+1} +\boldsymbol{ b} - \boldsymbol{ d}_2^k}, \frac{\lambda}{\mu}} \big) \end{align} $
(44) 式中, $\rm{soft}\big({\Omega}, \nu\big)$即为著名的软阈值函数. $\rm{soft}\big({\Omega}, \nu\big)$表示对矩阵$\Omega$中每个元素${\Omega}_{i, j}$作如下映射: ${\Omega}_{i, j}\mapsto{\rm sgn} ({\Omega}_{i, j})\max\{|{\Omega}_{i, j}|-\nu, 0\}$.
在矩阵的运算中, 最为耗时的是求解矩阵的逆.我们可以发现, 在问题(28)的求解中, 除了采用式(40)更新$\boldsymbol{ x}^k$时需要进行矩阵的逆运算, 其余变量的更新只要进行简单的矩阵加减乘除与元素大小的比较操作, 因此着重关注式(40)中矩阵$\Psi$的大小.根据式(30)与(31), $A\in{\bf R}^{2{\ N}\times{\ p}}, \boldsymbol{ c}\in{\bf R}^{p}$, 因此$\Psi\in{\bf R}^{ {\ p} \times {\ p} }$, $\ p$为图像中的端元个数.在高光谱图像中, 端元个数$p$几乎都在几十以内, 因此, 求解$\Psi^{-1}$的运算量很小, 并且在子问题的求解中, $\Psi^{-1}$为常值, 意味着在迭代过程中只需要计算一次$\Psi^{-1}$, 这意味着从时间与空间开销上考虑, 求逆运算对子问题的求解并不会产生很大的影响.求解问题(29)的方法与表 2中的步骤类似, 下面直接给出应用表 2中步骤求解时, 问题(29)对应的步骤2至步骤4的最优解表达形式:
$ \begin{align}%\tag(45) \begin{cases} %\tag(40) \boldsymbol{ x}^{k+1}={\Psi}^{-1}\boldsymbol{ c} ( {z_1^k+d_1^k}) -\\ \qquad\;\;\;\;\;{\Psi}^{-1}A^{\rm T} (\boldsymbol{ b}-\boldsymbol{ z}_2^k-\boldsymbol{ d}_2^k) \\ {\Psi} = \boldsymbol{ c}\boldsymbol{ c}^{\rm T}+{A}^{\rm T}A \end{cases} \end{align} $
(45) $ \begin{align}%\tag(46) {z}_1^{k+1}={\max}\big( {0, \boldsymbol{ c}^{\rm T} \boldsymbol{ x}^{k+1} -d_1^{k}+\frac{1}{\mu}} \big) \end{align} $
(46) $ \begin{align}%\tag(47) \boldsymbol{ z}_2^{k+1}=\rm{soft}\big( {{A \boldsymbol{ x}^{k+1} +\boldsymbol{ b} - \boldsymbol{ d}_2^k}, \frac{\lambda}{\mu}} \big) \end{align} $
(47) 关于表 2步骤7中的迭代停止条件, 当满足如下两种情形中的一种时, 停止迭代.
情形1. 迭代步数$k$到达最大迭代步数.
情形2. 原始残差与对偶残差(见下文第2.3.2节中定义)小于预设阈值.
特别要说明的是关于两个优化问题的优化值取舍的问题.如果$p^{*} < q^{*}$, 则采用子问题(28)的解, 否则采用子问题(29)的解.
最后, 假设RMVHU算法最终求解问题(27)得到的最终解为$H^*$与$\boldsymbol{ g}^*$, 则可由式(48)得到端元矩阵$A$:
$ \begin{align}%\tag(48) \begin{cases} A=CE+d1_{p}^{\rm T}\\ E=\big[{ B, \boldsymbol{ \beta}_{p} } \big]\\ B=H^{*-1}, \boldsymbol{ \beta}_{p}=B\boldsymbol{ g}^{*}\\ \end{cases} \end{align} $
(48) 式中, $C$与$\boldsymbol{d}$由式(14)得到.
根据式(21), 可得到第$n$个像元的丰度$\hat{\boldsymbol{s}}[n]$:
$ \begin{align}%\tag(49) \hat{\boldsymbol{ s}}[n]=\left[{ \begin{array}{cc} ( {H^*} \tilde{\boldsymbol{ y}}[n] -\boldsymbol{ g}^*)^{\rm T}&1-1_{ {\ p} -1}^{\rm T}( {H^*} \tilde{\boldsymbol{ y}}[n] -\boldsymbol{ g}^*) \end{array} } \right]^{\rm T} \end{align} $
(49) 式中, $\ p$为端元个数.
2.3 RMVHU算法参数的自适应调节策略
2.3.1 正则项系数的自适应调节策略
正则项系数的大小控制着正则项在整个优化函数中的权值, 也影响着整个算法的解混效果, 当正则项所占的比例较小时, 其起的作用相应较小, 反之亦然.在RMVHU算法中, 若正则项系数$\lambda$不变, 则在每一步的迭代中, 由于$\boldsymbol{ c}^{\rm T} \boldsymbol{ x}$项每次都会变化, 其所起的作用将会随每一次的更新而改变, 不利于整体优化问题的收敛.因此有必要固定其在每一步子优化问题中的占比.受文献[26]启发, 采用如下策略自适应调整参数$\lambda$:
$ \begin{align}%\tag(50) \lambda=\omega\frac{\left|\boldsymbol{ c}^{\rm T} \boldsymbol{ x}_0 \right|}{ \left\|{A \boldsymbol{ x}_{0} +\boldsymbol{ b} }\right\|_1 } \end{align} $
(50) 式中, $\boldsymbol{ x}_0$为$\boldsymbol{ x}$在迭代之前的初值. $\omega$为常值比例系数, 实验中发现其值在$30\, \sim\, 50$之间算法效果较好.
下面通过定性的分析, 表明采用上述的自适应参数调整法是行之有效的.
首先单独考虑体积约束项$\big|\boldsymbol{ c}^{\rm T} \boldsymbol{ x}_0 \big|$.若$\big|\boldsymbol{ c}^{\rm T} \boldsymbol{ x}_0 \big|$很大, 则表示端元构成的单形体体积已经很小了, 这时, 我们需要加重负数正则惩罚项的作用才能防止过多的像素点不满足非负条件, 进而给解混带来的较大误差.所以在体积约束项较大时, 负数惩罚正则项的系数应该随之增大.反之, 若$\big|\boldsymbol{ c}^{\rm T} \boldsymbol{ x}_0 \big|$很小, 则表明端元构成的单形体体积仍然很大, 负数正则惩罚项的作用还不是那么重要, 若我们相应地减小系数, 则可以加快算法的收敛速度.
其次单独考虑负数惩罚正则项${ \big\|{A \boldsymbol{ x}_{0} +\boldsymbol{ b} }\big\|_1 } $.若${ \big\|{A \boldsymbol{ x}_{0} +\boldsymbol{ b} }\big\|_1 } $较大, 则表示在数据集中, 很多像素点已经不满足非负性条件了, 实际显然不可能发生这样的情况, 所以此时需要减弱非负惩罚项的作用, 使单形体能够包含更多的像素点.而当${ \big\|{A \boldsymbol{ x}_{0} +\boldsymbol{ b} }\big\|_1 } $较小时, 说明单形体中基本包含了全部的像素点, 这时为了能够剔除出一些异常点, 需要增加正则惩罚项的权重, 也就需要相应增大正则项系数.
以上两点原因, 解释了为何在正则项的系数与体积约束项大小成正比, 与负数惩罚正则项大小成反比时, RMVHU算法能够取得很好的效果.
2.3.2 惩罚系数${\boldsymbol{\mu}}$的自适应调整
根据文献[22], 子问题ADMM形式(38)中的惩罚项系数$\mu$在如下简单自适应调节的策略下, 使算法的收敛性得到提高:
$ \begin{align}%\tag(51) \mu^{k+1} = \begin{cases} \tau\mu^k,& \left\| { \boldsymbol{ r}^k} \right\|_2 > \gamma\left\| { \boldsymbol{ s}^k} \right\|_2 \\ \frac{\mu^k}{\tau},& \left\| { \boldsymbol{ s}^k} \right\|_2 < \gamma\left\| { \boldsymbol{ r}^k} \right\|_2 \\ \mu^k,& \mbox{其他} \end{cases} \end{align} $
(51) 式中, $\boldsymbol{ r}^k$为原始残差, $\boldsymbol{ s}^k$为对偶残差.根据原始残差与对偶残差在文献[22]中的定义, 在RMVHU算法中, 原始残差与对偶残差能够由以下公式进行计算:
$ \begin{align}%\tag(52) { \boldsymbol{ r}^k} = \left[{ \begin{array}{c} {\boldsymbol{ c}^{\rm T}} \\ A \end{array} }\right] \boldsymbol{ x}^k - \left[{ \begin{array}{c} {z}_1^{k} \\ \boldsymbol{ z}_2^k \end{array} }\right] - \left[{ \begin{array}{c} 0 \\-\boldsymbol{ b} \end{array} } \right] \end{align} $
(52) $ \begin{align}%\tag(53) { \boldsymbol{ s}^k} = -\mu\left[{ \begin{array}{c} \boldsymbol{ c}^{\rm T} \\ A \end{array} }\right]^{\rm T} \left( \left[{ \begin{array}{c} {z}_1^{k} \\ \boldsymbol{ z}_2^k \end{array} }\right] - \left[{ \begin{array}{c} {z}_1^{k-1} \\ \boldsymbol{ z}_2^{k-1} \end{array} }\right] \right) \end{align} $
(53) 当惩罚系数$\mu$改变后, ${d}_1^{k}$与$\boldsymbol{ d}_2^k$也要相应地进行改变.当$\mu$翻倍时, ${d}_1^{k}$与$\boldsymbol{ d}_2^k$要相应减半, 反之亦然.关于参数$\mu$与$\tau$的选择对算法收敛性的影响, 将在第3.2节的实验中分析.
算法的主要流程如下表所示:
在表 3描述的算法流程中, 在更新矩阵变量$H$与$\boldsymbol{g}$时, 一次完整的迭代步骤为:采用步骤3至步骤5将矩阵变量的每一行元素都更新一遍.也就是说, 在算法的整个求解过程中, 分别最多求解${\ M}\times ({\ p} -1)$个子问题(28)与(29).算法在采用表 2描述的ADMM算法求解子问题后是相对高效的, 因此, RMVHU算法必然能够在有限时间内给出可行解.
表 3 RMVHU算法步骤Table 3 Procedure for solving RMVHU算法名称: RMVHU 输入: 端元个数$\ p$, 高光谱图像数据集$Y\in{\bf R}^{ {\ L} \times {\ N}}$, 收敛阈值$\varepsilon$, 最大迭代步数值$\ M$. 输出: 端元矩阵$A$, 丰度矩阵$S$. 步骤1: 根据式(14)与(16)对数据集$Y$进行降维, 得到降维后的光谱数据集$\tilde{Y}\in{\bf R}^{ ({\ p} -1) \times {\ N}}$. 步骤2: 设置当前迭代步数$i=0$, 并寻找满足如下条件的矩阵$H$与$\ g$作为初始迭代值: } 步骤3: 对$k=1, 2, \cdots, {\ p}-1$, 执行步骤4与步骤5. 步骤4: 使用表 2中的ADMM算法与式(51)自适应更新, 求解子问题(28)与(29), 得到单步最优解$(\boldsymbol{ h}_{k}^{-}, g_{k}^{-})$与$(\boldsymbol{ h}_{k}^{+}, g_{k}^{+})$及其对应的目标函数值$p^*$与$q^*$. 步骤5: 如果$p^* < q^*$, 则, 否则$(\boldsymbol{ h}_{k}, g_{k}) = (\boldsymbol{ h}_{k}^{+}, g_{k}^{+})$, 更新矩阵$H$与$\boldsymbol{ g}$相应的元素得到$H^*$与$\boldsymbol{ g}^*$. 步骤6: 如果$\left| {\left| {\det \left({{H^*}} \right)} \right| -\rho }\right|/\rho < \varepsilon$, 或者$i > {\ M}$, 则转步骤7, 否则$i=i+1$, , 应用式(50)更新正则系数$\lambda $, 转步骤3. 步骤7: 根据式(48)、(49)得到端元矩阵$A$与丰度矩阵$S$. 2.4 RMVHU算法时间复杂度与收敛性分析
为进一步研究RMVHU算法的效率, 本节将详细分析算法的复杂度.由于算法的本质为求解${\ M}\times ({\ p} -1)$个子问题(28)与(29), 因此, 将以求解子问题(28)的时间复杂度为分析的着眼点.
当采用表 2所描述的算法求解子问题(28)时, 更新$\boldsymbol{ x}$需做$2{Np}+4{\ N}+{\ p}$次加法或者乘法运算, 更新$z_1$需做${p}+2$次加法或者乘法运算, 更新$\boldsymbol{ z}_2$需做$2{\ Np}+12{N}$次加法或者乘法运算, 更新$d_1$需做${\ p}+2$次加法或者乘法运算, 更新$\boldsymbol{ d}_2$需做$2{\ Np}+6{\ N}$次加法或者乘法运算, 在计算$\big\| { \boldsymbol{ s}^k} \big\|_2$时, 需做$2{\ Np}+6{N}+{\ p}+2$次加法或者乘法运算, 在计算$\big\| { \boldsymbol{ r}^k} \big\|_2$时, 需做$2{\ Np}+4{\ N}+{\ p}+8$次加法或者乘法运算.因此, 在表 2描述的算法中, 循环一次需要做$10{\ Np}+32{\ N}+5{p}+8$次加法或者乘法运算, 算法的时间复杂度为相应的为O$({Np}+{\ N}+{\ p})$, 可以发现, 时间复杂度与图像像素点个数${N}$成正比, 与端元个数${\ p}$成正比, 当${\ p}\ll{N}$时, 算法的时间复杂度能够被近似表示为O$({\ N})$.
当采用定理1的方法, 将问题(28)转化为线性规划问题求解时, 变量个数为${4N}+{\ p}$.采用Karmarkar算法求解时, 时间复杂度[27]为O$({N}^{3.5})$.
显然, 采用ADMM算法求解在时间效率上具有显然的优势, 尤其当$ {N}$很大的时候.
为了分析RMVHU算法的收敛性, 本文将从如下两个方面进行阐述.
1) 在式(27)描述的优化目标函数中, 当固定正则项系数$\lambda$不变, 采用循环最小化的方法, 逐行更新$H$与$\boldsymbol{g}$时, 所求解的两个子问题都是凸问题, 并且ADMM算法在求解凸问题时具有收敛性, 因此, 式(27)的目标函数值必将逐渐减少, 说明在采用循环最小化的方法, 执行步骤4至步骤5时, 算法是收敛的.
2) 当执行完步骤4与步骤5, 算法仍未收敛时, 将根据式(50)的正则系数自适应调整策略, 重新计算$\lambda$, 这时, 根据第2.3.1节中的定性分析, 算法将在此种策略下加速收敛.
3. 仿真校验
本部分由两块内容组成, 首先介绍实验数据的来源以及算法评价准则, 随后给出实验结果并进行分析.
3.1 算法实验数据来源与评价准则
算法的实验数据分为合成光谱数据与真实光谱数据.
3.1.1 合成高光谱数据
在USGS光谱库[28]中选择$\ p$个光谱向量作为端元. 图 2显示了库中的3条光谱曲线.
根据文献[16]的方法, 产生像元的丰度信息.像元的丰度信息满足狄利克雷分布, 狄利克雷分布表示为:
$ \begin{align}%\tag(54) p({\alpha _1}, {\alpha _2}, \cdots, {\alpha _{p} })=\, &\frac{ \Gamma ({\mu _1} + {\mu _2} + \cdots + {\mu _{\ {p}} }) } {\Gamma ({\mu _1}) \Gamma ({\mu _2}) \cdots \Gamma ({\mu _{p}}) } \times \nonumber\\ &{\alpha _1^{\mu_1-1}}{\alpha _2^{\mu_2-1}}\cdots {\alpha _{p}^{ {\mu _{\ {p}} } -1 } } \end{align} $
(54) 式中, 第$i$个端元的丰度$\alpha _i$满足$0 \leq \alpha _i \leq 1$, $\sum_{i=1}^{p} \alpha_i = 1$, 且丰度期望${\rm E}[\alpha_i]={\mu _i}/ \sum_{k=1}^{p}{\mu _k}$. $\Gamma(\cdot)$表示Gamma函数.狄利克雷的密度分布不仅能够保证丰度符合非负与和为1约束, 还能够根据不同的参数${\mu_i}$, 产生不同分布特征的丰度系数.为了产生高度混合的数据, 对于丰度${\alpha_i}>\rho$的像元, 将其丰度重新匹配, 改为$1/{\ p}$. $\rho$为表示像元混合度的参数且满足$\rho\in{(0, 1)}$, $\rho$越小, 图像混合度越高.
为了使数据更具有真实性, 加入独立同分布的零均值高斯加性噪声.信噪比SNR表示为:
$ \begin{align}%\tag(55) {\rm SNR} =\frac{\sum\limits_{n=1}^{\ N} {\left\| {\boldsymbol{ y}[n]} \right\|}_2} {{\sigma ^2}{{LN}}} \end{align} $
(55) 式中, $\boldsymbol{ y}[n]$为真实的光谱数据, $\sigma ^2$代表高斯噪声的方差, $\ L$为波段维数, $\ N$为像素点个数.
在仿真数据中, 还应考虑异常点对算法的影响, 因此需要生成异常点.假设异常点的个数为$k$, 则重新生成第1至$k$个像素点的丰度.第$i\ (0 \leq i \leq k)$个点的丰度由以下公式重新产生:
$ \begin{align}%\tag(56) \boldsymbol{ s}(i) = \begin{cases} {{s_j}(i) = 1 + 0.2\delta },& j = k\\ {s_m}(i) = \frac{ {s_m}(i) } { \sum\limits_{m \ne j}{s_m}(i) }- \frac{1+0.2\delta}{{\ p}-1},&m \ne k \end{cases} \end{align} $
(56) 式中, $\delta$为常数, $k={\rm randperm}({{p}})$, 代表在1与$\ p$之间满足均匀分布的随机整数, $\boldsymbol{ s}(i)$为第$i$个像元丰度的向量表示, ${s_m}(i)$代表$\boldsymbol{ s}(i)$的第$m$个元素.
3.1.2 实际数据的来源
本文使用的实际数据是美国内华达州Cuprite地区1997年成像的224波段0.4 ${\rm \mu m}$至2.5 ${\rm \mu m}$, 大小为$250 \times 191$的子图[29].这组数据在高光谱图像端元提取的研究中广泛应用, 具有很好的代表性.由于水汽吸收的干扰和低信噪比的原因, 第1$\, \sim\, $2, 104$\, \sim\, $113, 148$\, \sim\, $167以及221$\, \sim\, $224波段的数据被剔除.伪彩色子图如图 3所示.
图 3 Cuprite地区高光谱图像伪彩色子图, R: 2.109 ${\rm \mu m}$, G: 2.209 ${\rm \mu m}$, B: 2.308 ${\rm \mu m}$Fig. 3 The Pseudo color subimage of the AVIRIS Cuprite dataset, R: 2.109 ${\rm \mu m}$, G: 2.209 ${\rm \mu m}$, B: 2.308 ${\rm \mu m}$该地区的矿物分布图由Tricorder 3.3软件[30]提供, 如图 4所示.需要注意的是此物质分布在1995年就已问世, 而Cuprite地区的数据是在1997年采集.因此该物质分布图仅作为解混效果好坏的参考.
3.1.3 算法评价准则
为了描述解混算法的效果好坏, 通常需要研究解混后端元与真实端元的相似度以及解混后估计丰度与真实丰度的相似程度.可以定义光谱角距离(SAD)和误差平方根(RMSE)来评估实验所得的端元和丰度.
对于第$i$个端元, 假设实际的端元的向量表示形式为$\boldsymbol{ a}_i$, 采用解混算法得到的端元估计值的向量表示为$\hat{\boldsymbol{ a}}_i $, 第$i$个端元的光谱角距离表示为$\phi_i$:
$ \begin{align}%\tag(57) {\phi_i} = {\cos ^{ - 1}}\left( { \frac{ \hat{\boldsymbol{ a}}_i^{\rm T}{\boldsymbol{ a}} } {\left\| {{{\hat {\boldsymbol{ a}}}_i}} \right\|\left\| {{\boldsymbol{ a}_i}} \right\|} } \right) \end{align} $
(57) 对于第$i$个端元, 假设端元对应的真实丰度表示为$\boldsymbol{ s}_i$, 采用解混算法得到与端元对应的丰度估计值表示为$\hat{\boldsymbol{ s}}_i$, $N$为图像像素点个数, 第$i$个端元的均方根误差表示为$\theta_i$:
$ \begin{align}%\tag(58) {\theta_i} = \left( \frac{1}{\ N} \left\| \boldsymbol{ s}_i-\hat{\boldsymbol{ s}}_i \right\|^{2}_{2} \right)^{\frac{1}{2}} \end{align} $
(58) 为全局地评价解混算法的好坏, 将采用如下两个平均指标:平均光谱角距离$\varepsilon_\phi$与平均误差平方根$\varepsilon_\theta$进行算法效果的衡量:
$ \begin{align}%\tag(59) {\varepsilon _\phi } = \sum\limits_{i = 1}^{p} \frac{\phi_i}{ p} \end{align} $
(59) $ \begin{align}%\tag(60) {\varepsilon_\theta } = \sum\limits_{i = 1}^{p} \frac{\theta_i}{\ p} \end{align} $
(60) 式(59)与(60)中, $\ p$为图像中端元个数.
3.2 实验及结果分析
实验分为合成数据解混与真实数据解混两组实验.在采用合成数据的实验时, 首先研究$\tau$与$\gamma$的取值对ADMM算法收敛的影响, 研究第2.3.1节中正则系数自适应调整的策略的有效性, 随后在端元数${p}=3$且存在异常点的数据集中, 将RMVHU算法与MVES算法进行比较, 阐述RMVHU算法的优点; 其次将RMVHU算法与被广泛使用的VCA、MVES、MVSA、SISAL四种凸面几何学算法进行比较, 分别研究真实端元个数, 噪声大小, 数据混合度$\rho$, 异常点个数, 像元个数以及错估端元个数时对不同算法解混效果的影响.最后, 对于真实数据, 给出RMVHU算法的解混结果.
3.2.1 合成数据的仿真
实验1. RMVHU惩罚系数的调整实验
本实验将进一步研究第2.3.2节中惩罚系数$\gamma$与$\tau$的取值对算法收敛性的影响.实验中, $\gamma$取值为1, 10, 100, 200, $\tau$取值为0.5, 1.5, 2, 4.
图 5与图 6为求解子问题(28)时, 目标函数值在不同$\gamma$与$\tau$下的变化趋势.从图 5可以发现, 当$\gamma$较小时, 函数收敛较快, 但波动较大, 当$\gamma$增大时, 目标函数虽然波动较小, 但收敛速度变慢.产生上述结果的原因在于, $\gamma$较小时, 只要原始残差与对偶残差存在差别, 那么惩罚系数$\mu$便会频繁的变化, 虽然这样的变化能使收敛速度提高, 但会带来收敛稳定性的不足; 而当$\gamma$较大时, 原始残差与对偶残差的差异性将很难影响到惩罚系数$\mu$的改变, 而合适的$\mu$将对收敛速度产生较大的影响, 因此, 虽然ADMM算法具有收敛性, 但是收敛速度将放缓.在图 5中可以发现, $\gamma=10$时, 函数收敛较快, 波动较小.
从图 6可以发现, 当$\tau < 1$时, 函数不收敛, 随着$\tau$的增大时, 函数收敛速度加快, 将很快收敛于极值, 但是波动将会增加, 其原因在于, 当对偶残差与原始残差较大时, 若$\tau < 1$, 根据调整策略, 惩罚系数$\mu$将增大, 这将导致对偶残差更大, 势必造成算法的不收敛; 当$\tau>1$, 则根据调整策略, 惩罚系数$\mu$将减小, 这将缩小对偶残差与原始残差的大小, 使得算法收敛性得到提升, 但是过大的$\tau$将使对偶残差波动较大, 带来收敛稳定性的不足.在图 6中可以发现$\tau=2$时, 函数收敛较快, 波动较小.
因此在接下来的实验中, 取$\gamma=10, \tau=2$.
实验2. 正则项系数的自适应调节对算法的影响
为了进一步研究第2.3.1节的正则项系数的自适应调节策略对算法的影响, 将进行相关实验进行对比验证.实验中将比较随着端元个数的变化时, 两种方法的解混精度.固定系数的算法取参数$\lambda=0.002$, 采用第2.3.1节中描述的策略的实验, 取参数$\omega = 40$.
实验中, 模拟数据的像素点个数${\ N}=10\, 000$, 波段维数${L}=220$, 端元数${\ p} $从3至8变化, SNR = 30 dB, 数据混合度$\rho = 0.8$, 异常点个数为25, 估计端元个数无错估.实验结果如图 7与图 8所示.
图 7 正则项系数固定或自适应变化时, 各算法平均光谱角距离Fig. 7 The average spectral angle distance of different algorithms when the regularization parameter is fixed or changes automatically
图 8 正则项系数固定或自适应变化时, 各算法平均均方根误差Fig. 8 The average root mean square error of different algorithms when the regularization parameter is fixed or changes automatically图 7与图 8表明, 当正则项系数固定时, 实验效果将会随着端元个数的改变发生很大的变化.而采用第2.3.1节策略的算法则能够稳定地估计真实端元与丰度.由此可以说明, 采用第2.3.1节策略将使算法具有更强的鲁棒性.
实验3. RMVHU的鲁棒性验证实验
本实验主要通过与一种流行的凸面几何学解混算法MVES进行比较, 直观地表明RMVHU算法在存在异常点情况下的鲁棒性.
模拟数据的像素点个数${\ N}=10\, 000$, 波段维数${\ L}=220$, 端元个数${\ p}=3$ (选取的端元光谱曲线见图 2), SNR = 30 dB, 数据混合度$\rho=0.8$, 异常点个数为25 (像素点总数的0.25%), 估计端元个数为3 (端元个数无错估).
RMVHU与MVES端元估计结果如图 9所示.在图 9中, 实线代表真实的端元光谱曲线, 虚线代表解混算法提取出的端元曲线.可以发现, 在少量异常点的干扰下, RMVHU算法的端元结果远远优于MVES算法.其主要原因在于, RMVHU算法引入了正则项, 该正则项能够容忍丰度值为负的情形, 使得噪声与异常点对最小体积逼近的影响不那么明显, 而MVES算法需要严格满足非负的条件, 因此MVES计算的单形体顶点需严格包含所有观测数据, 使异常点对最小体积逼近的影响尤为突出, 导致了端元提取效果的剧烈下降.
图 9 RMVHU与MVES端元估计结果, (a)、(b)、(c)为RMVHU的端元估计结果, (d)、(e)、(f)为MVES的端元估计结果Fig. 9 Estimated results of endmembers via two algorithms: RMVHU and MVES ((a)、(b)、(c) are the results of RMVHU and (d)、(e)、(f) are the results of MVES)图 10显示了数据在二维空间的分布以及RMVHU、MVES算法提取的端元在二维空间的分布, 此图亦印证了上述缘由.
图 10 数据集与端元在$p-1$ (二)维子空间的投影Fig. 10 Two dimensional subspace projection of datasets, estimated endmembers and real endmembers表 4与表 5分别记录了RMVHU与MVES算法的光谱角距离, 平均光谱角距离与均方根误差, 平均均方根误差值.结果最优值在表中加粗.从表 4与表 5中可以发现, 若图像数据中存在异常值, 则在进行图像的解混时, RMVHU算法无论是在端元提取精度上, 还是丰度反演精度上, 都远远优于MVES算法.随后进行的合成数据仿真实验中, 将分别研究真实端元个数, 噪声大小, 数据混合度, 有无异常点, 像素点个数以及错估端元个数对不同解混算法的影响.
表 4 RMVHU与MVES算法的光谱角距离与平均光谱角距离Table 4 Spectral angle distance and average spectral angle distance via RMVHU and MVES端元名称 RMVHU MVES Spessartine 0.30 1.86 Nontronite 0.26 18.29 Arsenopyrite 2.03 120.60 $\varepsilon_\theta$ (°) 0.87 46.92 表 5 RMVHU与MVES算法的均方根误差与平均均方根误差Table 5 Root mean square error and average root mean square error via RMVHU and MVES端元名称 RMVHU MVES Spessartine 0.012 0.239 Nontronite 0.008 0.151 Arsenopyrite 0.013 0.165 $\varepsilon_\theta$ (°) 0.011 0.185 实验4. 真实端元个数变化的实验
在每一幅高光谱图像中, 端元个数会根据不同的拍摄场景而变化, 因此, 研究端元个数对解混结果的影响很有必要.
本实验中, 生成端元个数${\ p}$从3$\sim$10变化的8组数据, 其他参数设置为:模拟数据的像素点个数${N}=10\, 000$, 波段维数${\ L}=220$, SNR = 30 dB, 数据混合度$\rho=0.8$, 异常点个数为25 (像素点总数的0.25%), 端元个数无错估.
平均光谱角距离结果如图 11所示.平均均方根误差结果如图 12所示.
图 11 端元个数变化时, 各算法平均光谱角距离Fig. 11 The average spectral angle distance of different algorithms when the number of endmembers changes
图 12 端元个数变化时, 各算法平均均方根误差Fig. 12 The average root mean square error of different algorithms when the number of endmembers changes图 11与图 12表明, 相较于其他算法, 端元个数的变化对RMVHU算法的影响并不很明显. RMVHU算法的解混效果相对于基于纯像元假设的算法——VCA得到的结果优秀很多.相比于需要严格满足非负性条件的最小体积变化算法: MVES与MVSA, RMVHU提取的{端元}的精度要优秀很多.对于同样可以不满足非负条件, 但是采用二次逼近体积项的SISAL算法来说, 由于RMVHU算法子问题中的体积项采用等价的行列式展开公式, 因此其求解的精度也在大多数情况下比SISAL算法高.
实验5. 抗噪能力的测试实验
设计本实验的目的在于测试在不同大小的噪声下, RMVHU算法的稳定性.
本实验中, 信噪比SNR分别设置为15 dB、20 dB、25 dB、30 dB、35 dB、40 dB, 获得6组数据, 其他参数设置为:模拟数据的像素点个数${\ N}=10\, 000$, 波段维数${L}=220$, 端元个数${\ p}=6$, 数据混合度$\rho=0.8$, 异常点个数为25 (像素点总数的0.25%), 端元个数无错估.
在不同噪声大小下, 平均光谱角距离与平均均方根误差结果如图 13与图 14所示.
图 13 在不同噪声大小下, 各算法平均光谱角距离Fig. 13 The average spectral angle distance of different algorithms when SNR changes
图 14 在不同噪声大小下, 各算法平均均方根误差Fig. 14 The average root mean square error of different algorithms when SNR changes从图 13与图 14中可以发现, 信噪比越高, RMVHU算法的解混效果越好.由于VCA中采用了根据不同噪声大小进行不同的数据降维方法, 因此VCA在低信噪比下有较突出的表现.而本文提出的RMVHU算法效果在低信噪比下的端元提取精度仅次于VCA, 侧面印证了RMVHU算法正则项及自适应的正则算子调节方法有较强的抗噪能力.相较于其他几何学解混算法, 在信噪比大于20 dB的情况下, RMVHU算法的表现出色, 能够得到最低的平均光谱角距离与平均均方根误差值. RMVHU算法较强的抗噪声能力得到了实验结果的支持.
实验6. 数据混合度变化的实验
设计本实验的目的在于:验证在不同混合度的场景下, RMVHU算法也有较为突出的表现, 即能够胜任不同混合程度的高光谱图像的解混.
在本实验中, 数据混合度$\rho$分别设置为0.65、0.7、0.75、0.8、0.85、0.9、0.95、1, 获得8组数据.其他参数设置如下:模拟数据的像素点个数${\ N}=10\, 000$, 波段维数${\ L}=220$, 真实端元个数${\ p}$为6, SNR为30 dB, 异常点个数为25 (像素点总数的0.25%), 端元个数无错估.
在不同混合度下, 平均光谱角距离与平均均方根误差结果如图 15与图 16所示.
图 15 在不同混合度下, 各算法平均光谱角距离Fig. 15 The average spectral angle distance of different algorithms when the purity of pixels changes
图 16 在不同混合度下, 各算法平均均方根误差Fig. 16 The average root mean square error of different algorithms when the purity of pixels changes由图 15与图 16发现, RMVHU在不同的混合度下, 都能取得最为优异的结果.同时, 由于异常值的存在, 使得本应该在混合度较低, 即$\rho$接近1的情况下能够获得更好效果的纯像元提取类算法VCA变得精度很差.但是本文提出的解混算法RMVHU在异常点存在的情况下, 对数据的混合度大小并不敏感.同时由于实验4中所述的原因, RMVHU相比于SISAL也更加优秀.本实验的成功意味着该算法在对不同混合度的高光谱图像解混时, 能够取得很好的效果.
实验7. 高于异常点个数变化的实验
设计本实验的目的在于:验证在异常点个数发生变化时, RMVHU算法有较强的稳定性, 并且与其他算法一起, 比较异常点的有无及数量对算法解混精度的影响.
实验中, 异常点个数分别设置为像素点总数的0%(无异常数据), 0.25%、0.5%、0.75%、1%、1.25%、1.5%、1.75%、2%, 获得9组数据.其他参数设置如下:模拟数据的像素点个数${\ N}=10\, 000$, 波段维数${\ L}=220$, 真实端元个数${\ p}$为6, SNR为30 dB, 数据混合度$\rho$设置为0.8, 端元个数无错估.
在异常点个数改变的情形下, 平均光谱角距离与平均均方根误差结果如图 17与图 18所示.
图 17 在不同异常点个数下, 各算法平均光谱角距离Fig. 17 The average spectral angle distance of different algorithms when the number of outliers changes
图 18 在不同异常点个数下, 各算法平均均方根误差Fig. 18 The average root mean square error of different algorithms when the number of outliers changes由图 17与图 18可以发现, 图像中无异常数据时, 基于最小体积变化的算法RMVHU、MVES、MVSA及SISAL算法都能够取得较好端元提取与解混的效果; 但是当出现异常数据时, 除了RMVHU算法能够取得最为优异的解混结果, SISAL算法勉强能够较为准确地获得图像端元与丰度信息外, 其他算法获得的结果都与真实值有较大偏差.同时, 相比于SISAL算法的结果随着异常点个数增加而变差, RMVHU算法的解混误差基本不变, 维持在最低的水平.本实验进一步说明了RMVHU算法的稳定性.
实验8. 关于像素点个数的变化的实验
设计本实验的目的在于研究图像像素点个数对RMVHU算法解混精度的影响以及算法的时间复杂度(即能否在较短的时间内得到较为精确的解混结果).
在本实验中, 设置像素点个数为2 000、4 000、6 000、8 000、10 000、12 000、14 000、16 000、18 000、20 000, 共得到10组数据.其他参数设置如下:真实端元个数${p}$为6, 波段维数${\ L}=220$, SNR为30 dB, 数据混合度$\rho$设置为0.8, 异常点个数为像素点总数的0.25%, 端元个数无错估.计算机为联想S20工作站.
在不同大小的高光谱图像中, 平均光谱角距离与平均均方根误差结果如图 19与图 20所示, 运行时间如图 21所示.
图 19 在不同像素点个数下, 各算法平均光谱角距离Fig. 19 The average spectral angle distance of different algorithms when the number of pixels changes
图 20 在不同像素点个数下, 各算法平均均方根误差Fig. 20 The average root mean square error of different algorithms when the number of pixels changes
图 21 在不同像素点个数下, 各算法运行时间Fig. 21 Computational time of different algorithms when the number of pixels changes由图 19与图 20可以发现, 相较于其他算法, RMVHU算法解混偏差并不会像SISAL、MVES、MVSA、VCA等算法随着像素点的个数变化出现巨大的波动, 它的表现一直很稳定.当像素点个数超过4 000时, RMVHU算法的解混结果是所有算法中最优秀的.此实验验证了本文提出的RMVHU算法能够在各种大小的高光谱图像解混中, 取得稳定且优秀的解混效果的说法.
在图 21中, 可以发现随着像素点个数的增多, RMVHU算法运行的时间线性增长, 这也与第2.4节中算法的复杂度分析相吻合.在像素点个数较大时, 由于不用求解大规模的线性规划问题, 其时间消耗反而要比MVES算法少.所以, 在有限时间内, 应用RMVHU算法, 定能得到理想的解混效果.相较于SISAL与MVSA等算法, 虽然其时间消耗较多, 但若更注重解混的精度与结果的稳定性, 那么在时间上的一些牺牲是值得的.
实验6. 关于错误估计端元个数的实验
在解混算法的实际应用中, 大多数解混算法需要预先估计端元的个数, RMVHU算法也不例外.但是可能由于混合模型的非线性偏差, 端元个数的估计有可能发生错误, 因此, 设计实验研究算法在错误估计端元个数时的表现也很有必要.
在本实验中, 真实端元的个数${ p}$在3$\sim$8间变化.估计的端元个数为${\ p} -1$, 模拟端元个数错误估计的情况.其他实验参数设置如下:模拟数据的像素点个数${\ N}=10\, 000$, 波段维数${\ L}=220$, SNR为30 dB, 数据混合度$\rho$设置为0.8, 异常点个数设置为0.
在对被错估端元个数的高光谱图像解混后, 平均光谱角距离与平均均方根误差结果如图 22与图 23所示.
图 22 在错估端元个数时, 各算法平均光谱角距离Fig. 22 The average spectral angle distance of different algorithms when the number of endmembers is incorrectly estimated
图 23 在错估端元个数时, 各算法平均均方根误差Fig. 23 The average root mean square error of different algorithms when the number of endmembers is incorrectly estimated由图 22发现在绝大多数端元个数错估的情况下, RMVHU算法相比其他的最小体积变换算法: MVSA、MVES及SISAL, 能够取得更为优秀的端元估计结果.
在对丰度值的估计方面, 由图 23可以发现, RMVHU算法在不同的端元个数下都获得了不错的效果, 且随着真实端元个数的增加, 精度也有提高的趋势.对RMVHU算法而言, 端元提取误差与丰度反演的误差随着真实端元个数的增加, 也有降低的趋势, 这也说明了RMVHU算法在错误估计端元个数的情况下仍然具有较强的鲁棒性.
3.2.2 真实数据的实验
采用Cuprite子图, 剔除干扰波段后, 文献[8]的研究结果表明, 这片区域中有14种矿物.由于有些同种类不同化学成分的矿物的光谱只存在很微小的差别, 因此在解混时, 将端元的个数减少至11.解混后各物质丰度图如图 24所示.
在图 24的Cuprite地区估计丰度图中, (a)为Desert Varnish的丰度图, (b)为Sphene的丰度图, (c)为Nontronite的丰度图, (d)为Kaolinite的丰度图, (e)为Dumortierite的丰度图, (f)为Chalcedony的丰度图, (g)为Pyrope的丰度图, (h)为Andradite的丰度图, (i)为Montmorillonite的丰度图, (j)为Muscovite的丰度图, (k)为Alunite的丰度图.
同时为了比较端元提取效果, 将RMVHU、VCA、MVES、MVSA、SISAL算法的端元提取结果与光谱库中的标准数据对比, 得到各物质的光谱角距离, 如表 6所示.针对所有物质, 各算法得到的光谱角距离的最小值已在表中加粗.
表 6 RMVHU、VCA、MVES、MVSA、SISAL算法提取端元与真实端元的光谱角距离Table 6 The spectral angle distance of real datasets via unmixing algorithms: RMVHU, VCA, MVES, MVSA and SISAL物质 RMVHU VCA MVES MVSA SIS (a) Desert Varnish 15.08 21.48 16.86 16.02 22.35 (b) Sphene 86.03 28.07 62.25 59.24 18.33 (c) Nontronite 33.57 18.22 35.51 33.13 33.80 (d) Kaolinite 20.32 31.95 29.44 37.59 14.90 (e) Dumortierite 14.55 24.94 23.75 21.60 64.74 (f) Chalcedony 15.15 16.26 38.04 39.76 28.07 (g) Pyrope 23.74 29.89 44.67 54.42 63.03 (h) Andradite 17.18 32.38 45.62 60.30 27.50 (i) Montmorillonite 21.24 46.39 35.85 41.08 21.20 (j) Muscovite 20.66 37.47 65.06 78.46 24.64 (k) Alunite 15.72 37.55 38.57 29.36 12.03 ${\varepsilon_\phi}$ (°) 25.75 29.51 39.60 42.81 30.05 将图 24与图 4相比, 可以发现, RMVHU算法确实能够将高光谱图像中的物质信息有效地提取出来.
从表 6中可以发现, RMVHU算法在对如下6种物质: Desert Varnish、Dumortierite、Chalcedony、Pyrope、Andradite以及Muscovite的端元提取中, 表现最为优秀.对于其他的物质, 除了Sphene, 另外四种物质的端元提取效果也较为优秀, RMVHU算法的端元提取结果也非常接近由其他算法所获得的最优秀的结果.最后可以发现, RMVHU算法的平均光谱角误差是所有算法中最小的, 证明了本算法在实际解混应用中性能的优越性.
4. 结论
本文通过引入负数惩罚正则项, 提出了自适应鲁棒最小体积高光谱解混算法(RMVHU), 在求解时采用循环最小化方法拆分成可以求解的凸优化子问题, 并成功应用ADMM框架对子问题进行求解.此外, 文章对正则项的系数$\lambda$提出了一种自适应的调节策略, 加快了收敛速度, 提高了算法的稳定性.
仿真实验表明相较于其他解混算法, RMVHU算法在任意端元个数下均拥有优秀且稳定的解混表现; 此外, 算法对数据噪声大小, 数据异常点个数以及图像像素点个数亦很不敏感, 尤其在端元数估计错误的情况下也有比较稳定的表现; 在不同混合度的场景下RMVHU也有很优秀的表现.
对Cuprite真实数据的解混结果表明, RMVHU算法提取出的端元及估算出的丰度图与实际地物组成基本一致.
需要注意的是, RMVHU算法的计算时间相对于诸如VCA等纯像元提取算法较长, 因此在下一步的研究中将尝试采用并行化技术, 以提高运算的效率.
致谢: 作者感谢审稿专家在论文修改过程中提出的宝贵建议和意见.作者对哈尔滨工业大学航天学院检测实验室的研究支持表示感谢.
-
图 1 顶点个数为3时, 凸包与仿射集的概念说明
Fig. 1 Illustration of convex hull and affine hull when the number of vertices is three
图 3 Cuprite地区高光谱图像伪彩色子图, R: 2.109 ${\rm \mu m}$, G: 2.209 ${\rm \mu m}$, B: 2.308 ${\rm \mu m}$
Fig. 3 The Pseudo color subimage of the AVIRIS Cuprite dataset, R: 2.109 ${\rm \mu m}$, G: 2.209 ${\rm \mu m}$, B: 2.308 ${\rm \mu m}$
图 7 正则项系数固定或自适应变化时, 各算法平均光谱角距离
Fig. 7 The average spectral angle distance of different algorithms when the regularization parameter is fixed or changes automatically
图 8 正则项系数固定或自适应变化时, 各算法平均均方根误差
Fig. 8 The average root mean square error of different algorithms when the regularization parameter is fixed or changes automatically
图 9 RMVHU与MVES端元估计结果, (a)、(b)、(c)为RMVHU的端元估计结果, (d)、(e)、(f)为MVES的端元估计结果
Fig. 9 Estimated results of endmembers via two algorithms: RMVHU and MVES ((a)、(b)、(c) are the results of RMVHU and (d)、(e)、(f) are the results of MVES)
图 10 数据集与端元在$p-1$ (二)维子空间的投影
Fig. 10 Two dimensional subspace projection of datasets, estimated endmembers and real endmembers
图 11 端元个数变化时, 各算法平均光谱角距离
Fig. 11 The average spectral angle distance of different algorithms when the number of endmembers changes
图 12 端元个数变化时, 各算法平均均方根误差
Fig. 12 The average root mean square error of different algorithms when the number of endmembers changes
图 13 在不同噪声大小下, 各算法平均光谱角距离
Fig. 13 The average spectral angle distance of different algorithms when SNR changes
图 14 在不同噪声大小下, 各算法平均均方根误差
Fig. 14 The average root mean square error of different algorithms when SNR changes
图 15 在不同混合度下, 各算法平均光谱角距离
Fig. 15 The average spectral angle distance of different algorithms when the purity of pixels changes
图 16 在不同混合度下, 各算法平均均方根误差
Fig. 16 The average root mean square error of different algorithms when the purity of pixels changes
图 17 在不同异常点个数下, 各算法平均光谱角距离
Fig. 17 The average spectral angle distance of different algorithms when the number of outliers changes
图 18 在不同异常点个数下, 各算法平均均方根误差
Fig. 18 The average root mean square error of different algorithms when the number of outliers changes
图 19 在不同像素点个数下, 各算法平均光谱角距离
Fig. 19 The average spectral angle distance of different algorithms when the number of pixels changes
图 20 在不同像素点个数下, 各算法平均均方根误差
Fig. 20 The average root mean square error of different algorithms when the number of pixels changes
图 21 在不同像素点个数下, 各算法运行时间
Fig. 21 Computational time of different algorithms when the number of pixels changes
图 22 在错估端元个数时, 各算法平均光谱角距离
Fig. 22 The average spectral angle distance of different algorithms when the number of endmembers is incorrectly estimated
图 23 在错估端元个数时, 各算法平均均方根误差
Fig. 23 The average root mean square error of different algorithms when the number of endmembers is incorrectly estimated
表 1 数学符号及其意义
Table 1 Mathematical notations and their meaning
符号 符号意义 ${\bf R}$ 实数 ${\bf R}^{{L}}$ ${L}$维向量 ${\bf R}^{{L}\times {\rm N}}$ ${L}\times {\rm N}$维矩阵 ${\bf R}_{+}$ 非负实数 ${\bf R}_+^{L}$ 非负${L}$维向量 ${\bf R}_{+}^{{L}\times {N}}$ 非负${L}\times {N}$维矩阵 $1_{N}$ ${N}\times 1$维全1列向量 $I_{N}$ ${N}\times {N}$维单位矩阵 $X^†$ $X$的广义逆矩阵 $\boldsymbol{ e}_i$ 第$i$个元素为1, 其他元素为0的单位列向量 $\boldsymbol{x}\succeq {a}, {\boldsymbol{x}}\in {\bf R}^{{L}}$ $L$维列向量$\boldsymbol{x}$中所有元素均不小于$a$ $l_{\bf R_{+}}(\boldsymbol{x})$ $l_{\bf R_{+}}(\boldsymbol{x})=\begin{cases} 0, &\boldsymbol{x}\succeq 0 \\ + \infty, & \mbox{其他} \end{cases}$ ${\rm sgn} (y)$ ${\rm sgn}(y)=\begin{cases} 1, &y > 0 \\ 0, & y=0\\-1, &y < 0 \end{cases}$ 
下载: 导出CSV
表 2 RMVHU中求解子问题$p^{*}$的算法步骤
Table 2 Steps for solving subproblem $p^{*}$ in RMVHU
算法名称: 求解问题$p^{*}$的ADMM算法步骤 步骤1: 设初值$k=0, \mu > 0, \boldsymbol{ x}_0, {z}_1^0, {d}_1^0, \boldsymbol{ z}_2^0, \boldsymbol{ d}_2^0 $ 步骤2: 更新$\boldsymbol{x}$: $\boldsymbol{ x}^{k+1} = \arg\; \mathop {\min }\limits_\boldsymbol{ x}L(\boldsymbol{ x}, {z_1^k}, \boldsymbol{ z}_2^k, d_1^k, \boldsymbol{ d}_2^k) $ 步骤3: 更新$z_1$: ${z_1^{k+1} }= \arg\; \underset{z_1}{\min}L(\boldsymbol{ x}^{k+1}, {z_1}, \boldsymbol{ z}_2^k, d_1^k, \boldsymbol{ d}_2^k) $ 步骤4: 更新$\boldsymbol{ z}_2$: $\boldsymbol{ z}_2^{k+1} = \arg\; \underset{\boldsymbol{ z}_2}{\min}L(\boldsymbol{ x}^{k+1}, {z_1^{k+1}}, \boldsymbol{ z}_2, d_1^k, \boldsymbol{ d}_2^k) $ 步骤5: 更新$d_1$: ${d_1^{k+1} }= d_1^k -(\boldsymbol{ c}^{\rm T} \boldsymbol{ x}^{k+1} -z_1^{k+1})$ 步骤6: 更新$\boldsymbol{ d}_2$: $\boldsymbol{ d}_2^{k+1} = \boldsymbol{ d}_2^{k} -(A \boldsymbol{ x}^{k+1} +\boldsymbol{ b} -\boldsymbol{ z}_2^{k+1})$ 步骤7: 判断是否满足迭代停止条件, 若满足, 则转步骤8, 否则$k=k+1$, 转步骤2. 步骤8: 输出$\boldsymbol{ x}^{k+1}$:
下载: 导出CSV
表 3 RMVHU算法步骤
Table 3 Procedure for solving RMVHU
算法名称: RMVHU 输入: 端元个数$\ p$, 高光谱图像数据集$Y\in{\bf R}^{ {\ L} \times {\ N}}$, 收敛阈值$\varepsilon$, 最大迭代步数值$\ M$. 输出: 端元矩阵$A$, 丰度矩阵$S$. 步骤1: 根据式(14)与(16)对数据集$Y$进行降维, 得到降维后的光谱数据集$\tilde{Y}\in{\bf R}^{ ({\ p} -1) \times {\ N}}$. 步骤2: 设置当前迭代步数$i=0$, 并寻找满足如下条件的矩阵$H$与$\ g$作为初始迭代值: $ {H} \tilde{\boldsymbol{ y}}[n] -\boldsymbol{ g} \succeq 0, 1_{ {\ p} -1}^{\rm T}({H} \tilde{\boldsymbol{ y}}[n] -\boldsymbol{ g})\leq 1, \forall n=1, \cdots, {\ N} $} 步骤3: 对$k=1, 2, \cdots, {\ p}-1$, 执行步骤4与步骤5. 步骤4: 使用表 2中的ADMM算法与式(51)自适应更新, 求解子问题(28)与(29), 得到单步最优解$(\boldsymbol{ h}_{k}^{-}, g_{k}^{-})$与$(\boldsymbol{ h}_{k}^{+}, g_{k}^{+})$及其对应的目标函数值$p^*$与$q^*$. 步骤5: 如果$p^* < q^*$, 则$(\boldsymbol{ h}_{k}, g_{k}) = (\boldsymbol{ h}_{k}^{-}, g_{k}^{-})$, 否则$(\boldsymbol{ h}_{k}, g_{k}) = (\boldsymbol{ h}_{k}^{+}, g_{k}^{+})$, 更新矩阵$H$与$\boldsymbol{ g}$相应的元素得到$H^*$与$\boldsymbol{ g}^*$. 步骤6: 如果$\left| {\left| {\det \left({{H^*}} \right)} \right| -\rho }\right|/\rho < \varepsilon$, 或者$i > {\ M}$, 则转步骤7, 否则$i=i+1$, $\rho = {\det \left({{H^*}} \right)}$, 应用式(50)更新正则系数$\lambda $, 转步骤3. 步骤7: 根据式(48)、(49)得到端元矩阵$A$与丰度矩阵$S$.
下载: 导出CSV
表 4 RMVHU与MVES算法的光谱角距离与平均光谱角距离
Table 4 Spectral angle distance and average spectral angle distance via RMVHU and MVES
端元名称 RMVHU MVES Spessartine 0.30 1.86 Nontronite 0.26 18.29 Arsenopyrite 2.03 120.60 $\varepsilon_\theta$ (°) 0.87 46.92
下载: 导出CSV
表 5 RMVHU与MVES算法的均方根误差与平均均方根误差
Table 5 Root mean square error and average root mean square error via RMVHU and MVES
端元名称 RMVHU MVES Spessartine 0.012 0.239 Nontronite 0.008 0.151 Arsenopyrite 0.013 0.165 $\varepsilon_\theta$ (°) 0.011 0.185
下载: 导出CSV
表 6 RMVHU、VCA、MVES、MVSA、SISAL算法提取端元与真实端元的光谱角距离
Table 6 The spectral angle distance of real datasets via unmixing algorithms: RMVHU, VCA, MVES, MVSA and SISAL
物质 RMVHU VCA MVES MVSA SIS (a) Desert Varnish 15.08 21.48 16.86 16.02 22.35 (b) Sphene 86.03 28.07 62.25 59.24 18.33 (c) Nontronite 33.57 18.22 35.51 33.13 33.80 (d) Kaolinite 20.32 31.95 29.44 37.59 14.90 (e) Dumortierite 14.55 24.94 23.75 21.60 64.74 (f) Chalcedony 15.15 16.26 38.04 39.76 28.07 (g) Pyrope 23.74 29.89 44.67 54.42 63.03 (h) Andradite 17.18 32.38 45.62 60.30 27.50 (i) Montmorillonite 21.24 46.39 35.85 41.08 21.20 (j) Muscovite 20.66 37.47 65.06 78.46 24.64 (k) Alunite 15.72 37.55 38.57 29.36 12.03 ${\varepsilon_\phi}$ (°) 25.75 29.51 39.60 42.81 30.05
下载: 导出CSV
-
[1] 潘宗序, 禹晶, 肖创柏, 孙卫东.基于光谱相似性的高光谱图像超分辨率算法.自动化学报, 2014, 40(12):2797-2807 http://www.aas.net.cn/CN/abstract/abstract18558.shtmlPan Zong-Xu, Yu Jing, Xiao Chuang-Bai, Sun Wei-Dong. Spectral similarity-based super resolution for hyperspectral images. Acta Automatica Sinica, 2014, 40 (12):2797-2807 http://www.aas.net.cn/CN/abstract/abstract18558.shtml [2] 倪鼎, 马洪兵.基于近邻协同的高光谱图像谱-空联合分类.自动化学报, 2015, 41 (2):273-284 http://www.aas.net.cn/CN/abstract/abstract18606.shtmlNi Ding, Ma Hong-Bing. Spectral-spatial classification of hyperspectral images based on neighborhood collaboration. Acta Automatica Sinica, 2015, 41 (2):273-284 http://www.aas.net.cn/CN/abstract/abstract18606.shtml [3] Du B, Zhang L P. A discriminative metric learning based anomaly detection method. IEEE Transactions on Geoscience and Remote Sensing, 2014, 52 (11):6844-6857 doi: 10.1109/TGRS.2014.2303895 [4] Lin C H, Chi C Y, Wang Y H, Chan T H. A fast hyperplane-based minimum-volume enclosing simplex algorithm for blind hyperspectralunmixing. IEEE Transactions on Signal Processing, 2016, 64 (8):1946-1961 doi: 10.1109/TSP.2015.2508778 [5] Bioucas-Dias J M, Plaza A, Dobigeon N, Parente M, Du Q, Gader P, Chanussot J. Hyperspectralunmixing overview:geometrical, statistical, and sparse regression-based approaches. IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing, 2012, 5 (2):354-379 doi: 10.1109/JSTARS.2012.2194696 [6] Lillesand T, Kiefer R W, Chipman J. Remote Sensing and Image Interpretation (Seventh Edition). New York:John Wiley and Sons, 2014. 23-61 [7] Yuan Y, Fu M, Lu X Q. Substance dependence constrained sparse NMF for hyperspectralunmixing. IEEE Transactions on Geoscience and Remote Sensing, 2015, 53 (6):2975-2986 doi: 10.1109/TGRS.2014.2365953 [8] Wang W H, Qian Y T, Tang Y Y. Hypergraph-regularized sparse NMF for hyperspectralunmixing. IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing, 2016, 9 (2):681-694 doi: 10.1109/JSTARS.2015.2508448 [9] Qu Q, Nasrabadi N M, Tran T D. Subspace vertex pursuit:a fast and robust near-separable nonnegative matrix factorization method for hyperspectralunmixing. IEEE Journal of Selected Topics in Signal Processing, 2015, 9 (6):1142-1155 doi: 10.1109/JSTSP.2015.2419184 [10] Li J, Bioucas-Dias J M, Plaza A, Liu L. Robust collaborative nonnegative matrix factorization for hyperspectralunmixing. IEEE Transactions on Geoscience and Remote Sensing, 2016, 54 (10):6076-6090 doi: 10.1109/TGRS.2016.2580702 [11] Feng R Y, Zhong Y F, Zhang L P.An improved nonlocal sparse unmixing algorithm for hyperspectral imagery. IEEE Geoscience and Remote Sensing Letters, 2015, 12 (4):915-919 doi: 10.1109/LGRS.2014.2367028 [12] Iordache M D, Bioucas-Dias J M, Plaza A. Collaborative sparse regression for hyperspectralunmixing. IEEE Transactions on Geoscience and Remote Sensing, 2014, 52 (1):341-354 doi: 10.1109/TGRS.2013.2240001 [13] Meyer T R, Drumetz L, Chanussot J, Bertozzi A L, Jutten C. Hyperspectralunmixing with material variability using social sparsity. In:Proceedings of 2016 IEEE International Conference on Image Processing. Phoenix, USA:IEEE, 2016. 2187-2191 [14] Giampouras P V, Themelis K E, Rontogiannis A A, Koutroumbas K D. Simultaneously sparse and low-rank abundance matrix estimation for hyperspectral image unmixing. IEEE Transactions on Geoscience and Remote Sensing, 2016, 54 (8):4775-4789 doi: 10.1109/TGRS.2016.2551327 [15] Zheng C Y, Li H, Wang Q, Chen C L P. Reweighted sparse regression for hyperspectralunmixing. IEEE Transactions on Geoscience and Remote Sensing, 2016, 54 (1):479-488 doi: 10.1109/TGRS.2015.2459763 [16] Nascimento J M P, Dias J M B. Vertex component analysis:a fast algorithm to unmixhyperspectral data. IEEE Transactions on Geoscience and Remote Sensing, 2005, 43 (4):898-910 doi: 10.1109/TGRS.2005.844293 [17] Boardman J W. Automating spectral unmixing of AVIRIS data using convex geometry concepts. In:Proceedings of Summaries of the 4th Annual JPL Airborne Geoscience Workshop. Arlington, Virginia:JPL, 1993. 11-14 [18] Winter M E. N-FINDR:an algorithm for fast autonomous spectral end-member determination in hyperspectral data. In:Proceedings of 1999 SPIE's International Symposium on Optical Science, Engineering, and Instrumentation. Denver, USA:SPIE, 1999. 266-275 [19] Chan T H, Chi C Y, Huang Y M, Ma W K. A convex analysis-based minimum-volume enclosing simplex algorithm for hyperspectralunmixing. IEEE Transactions on Signal Processing, 2009, 57 (11):4418-4432 doi: 10.1109/TSP.2009.2025802 [20] Li J, Agathos A, Zaharie D, Bioucas-Dias J M, Plaza A, Li X. Minimum volume simplex analysis:a fast algorithm for linear hyperspectralunmixing. IEEE Transactions on Geoscience and Remote Sensing, 2015, 53 (9):5067-5082 doi: 10.1109/TGRS.2015.2417162 [21] Bioucas-Dias J M. A variable splitting augmented Lagrangian approach to linear spectral unmixing. In:Proceedings of the 1st Workshop on Hyperspectral Image and Signal Processing:Evolution in Remote Sensing. Grenoble, France:IEEE, 2009. 1-4 [22] Boyd S, Parikh N, Chu E, Peleato B, Eckstein J. Distributed optimization and statistical learning via the alternating direction method of multipliers. Foundations and Trends in Machine Learning, 2010, 3 (1):1-122 doi: 10.1561/2200000016 [23] Chan T H, Ma W K, Chi C Y, Wang Y. A convex analysis framework for blind separation of non-negative sources. IEEE Transactions on Signal Processing, 2008, 56 (10):5120-5134 doi: 10.1109/TSP.2008.928937 [24] StrangG. Linear Algebra and Its Applications. San Diego, CA:Thomson, 2005. http://math.ecnu.edu.cn/~zhan/papers/HZ2015a.pdf [25] Chen S S, Donoho D L, Saunders M A. Atomic decomposition by basis pursuit. SIAM Review, 2001, 43 (1):129-159 doi: 10.1137/S003614450037906X [26] Shi Z W, An Z Y, Tan X Y, Zhu Z X, Jiang Z G. Hyperspectralunmixing using non-negative matrix factorization with automatically estimating regularization parameters. In:Proceedings of the 7th International Conference on Natural Computation. Shanghai, China:IEEE, 2011. 1836-1840 [27] 陈宝林.最优化理论与算法.第2版.北京:清华大学出版社, 2005. 180-193Chen Bao-Lin. Optimization Theories and Algorithms (Second Edition). Beijing:Tsinghua University Press, 2005. 180-193 [28] Clark R N, Swayze G A, Wise R, Livo E, Hoefen T, Kokaly R, Sutley S J. USGS digital spectral library[Online], available:http://speclab.cr.usgs.gov/spectral-lib.html , September 13, 2016. [29] AVIRIS. AVIRIS data-ordering free AVIRIS standard data products[Online], available:http://aviris.jpl.nasa.gov/html/aviris.freedata.html , September 13, 2016. [30] Clark R N, Swayze G A, Livo K E, Kokaly R F, Sutley S J, Dalton J B, McDougal R R, Gent C A. Imaging spectroscopy:earth and planetary remote sensing with the usgstetracorder and expert systems[Online], available:http://speclab.cr.usgs.gov/PAPERS/tetracorder , September 13, 2016. 期刊类型引用(6)
1. 王丽，王威，刘勃妮. 顶点成分分析法应用于高光谱图像端元提取. 信息技术. 2022(01): 37-43 . 
百度学术2. 黄鸿，唐玉枭，段宇乐. 面向高光谱分类的局部几何稀疏保持嵌入. 自动化学报. 2022(10): 2496-2507 . 本站查看3. 闫钧华，黄伟，张寅，许祯瑜，苏恺. 基于空间约束的高光谱图像解混算法. 仪器仪表学报. 2019(03): 188-195 . 百度学术4. 刘颖，梁楠楠，李大湘，杨凡超. 基于光谱距离聚类的高光谱图像解混算法. 计算机应用. 2019(09): 2541-2546 . 百度学术5. 许宁，胡玉新，耿修瑞. 凸体几何光谱解混研究进展及若干问题浅析. 遥感技术与应用. 2019(05): 1028-1039 . 百度学术6. 邹丽，蔡希彪，孙静，孙福明. 基于双图正则的半监督NMF混合像元解混算法. 计算机科学. 2018(12): 251-254+278 . 百度学术其他类型引用(7)
-
下载:
计量
- 文章访问数: 3898
- HTML全文浏览量: 278
- PDF下载量: 416
- 被引次数: 13
