1、引　言

空间分辨力的先天不足，使得高光谱图像中单一像元往往包含多种实际地物类型，形成所谓“混合像元”(Xu 等，2001)。混合像元的“解混”是提高地物分类与目标识别精度，促进高光谱遥感定量化应用的关键问题之一。

近年来，非负矩阵分解NMF (Non-negative Matrix Factorization)技术逐渐成为高光谱解混领域的研究热点。NMF由Lee和Seung(1999)在Nature杂志上正式提出，根据绝大多数物质的“非负”属性，以及由“非负”导致的“全加性”约束和“局部构成整体”思想，将高阶信号的解混转化为两个相乘的低阶非负矩阵的分解，显著降低了信号解混的难度与计算量(黄春海，2014)。NMF适合处理高维海量数据，同时，作为一种代表性的盲信号处理技术，能够降低解混过程对先验知识的依赖程度，符合数据量大且缺少先验知识的高光谱解混的基本需求。但通用性和局部极小等问题导致NMF的解混性能并不令人满意(Geng 等，2016)，需要引入符合具体应用特征的约束条件，弥补上述问题带来的不利影响。例如，Pauca等人(2006)提出了加入端元光谱平滑性和丰度稀疏性约束的NMF像元解混算法，Jia和Qian(2009)也提出了分段平滑和稀疏约束的NMF PSNMFSC(Piecewise Smoothness NMF with Sparseness Constraints)，这两种算法都需要预先指定稀疏度，实际应用中实现难度较大；Zymnis等人(2007)提出了交互投影子梯度的非负矩阵分解算法，但算法复杂度较高；Miao和Qi(2007)提出了最小体积约束的非负矩阵分解像元解混算法MVCNMF (Minimum Volume Constrained Nonnegative Matrix Factorization)，问题在于约束条件值域为无穷大，容易导致算法发散。因此，有必要从其他角度继续探索NMF解混算法中约束条件的构造原理和实现方式，使之更加合理和有效。

相关关系普遍存在于世间万物中，高光谱图像也不例外。近年来，相关性分析正越来越多地出现在高光谱图像的不同应用领域，分析对象涉及高光谱图像的光谱特征、空间特征、统计参数及各种物理参数等(Liu 等，2014)，基于空间相关或光谱相关特征的高光谱解混研究也取得了大量学术成果。例如，詹锡兰和吴波(2011)利用高斯马尔可夫随机场模型来描述组分空间相关特性，提出了一种基于高斯马尔可夫随机场模型的混合像元分解方法；杜世强等人(2012)提出了一种基于图正则化的半监督NMF算法(GSNMF)，克服了NMF忽略样本数据局部几何结构的缺陷；刘建军等人(2012)利用训练数据构造结构化字典，添加空间相关性约束项和训练数据的空间信息，提出了一种新的高光谱解混方法；Chen等人(2013)利用先验概率密度函数表达两个相邻区域的空间相关程度，提出了一种区域相关的NMF解混算法；王楠等人(2014)通过联合最小化光谱相关度函数和NMF误差函数，提出了一种最小化光谱相关度约束的NMF方法；姜小燕等人(2016)提出了一种基于稀疏约束和图正则化的半监督NMF方法，在进行低维非负分解时，能够保持数据几何结构；许宁等人(2016)基于图像中存在的邻域以及非局部相似等图像空间特征和联合稀疏解混思想，提出一种基于高光谱图像空间相关性度量的多任务联合稀疏解混方法。

上述算法主要存在以下两方面的问题：一是没有对图像在空间和谱间两方面的典型相关性特征进行充分利用，对应文献包括詹锡兰和吴波(2011)，王楠等人(2014)，许宁等人(2016)等；二是不属于非监督算法，所需的大量先验知识通常难以获取，导致算法实用性不够理想，对应文献包括杜世强等人(2012)，刘建军等人(2012)，Chen等(2013)，姜小燕等人(2016)等。

因此，有必要深入挖掘并充分利用高光谱遥感图像中蕴含的各类相关性特征，探索其在降低先验知识依赖程度，改善解混算法性能方面的应用潜力。

本文在线性光谱混合模型基础上，针对NMF在高光谱解混中存在的通用性和局部极小问题，综合利用马尔可夫随机场MRF (Markov Random Field)模型(Chen 等，2014)和信号预测度技术(尤传雨 等，2016)，提出一种空间与谱间相关性联合约束的NMF解混算法NMFSSC (NMF hyperspectral unmixing algorithm based on Spatial and Spectral Correlation)。

2、NMF解混原理及其在像元解混中存在的主要问题     (2.1) 线性光谱混合模型

进行NMF像元解混之前，首先需要确定光谱混合模型SMM (Spectral Mixture Model)。SMM可以简单分为两大类，即线性光谱混合模型LSMM (Linear Spectral Mixture Model)和非线性光谱混合模型NLSMM (Nonlinear Spectral Mixture Model)。相对于NLSMM，LSMM形式简单，物理意义明确，应用也更广泛(Zhao 等，2015)，本文算法主要基于LSMM进行设计。

LSMM数学形式如式(1)

${{X}} = \sum\limits_{i = 1}^p {{{{M}}_i}{{{S}}_i}} + {{N}} = {{MS}} + {{N}}$

(1)

式中， ${{X}}$ 是 $l$ ( $l$ 为波段数量)维的已知高光谱图像； ${{M}}$ 是 $l \times p$ ( $p$ 为端元数目)的端元矩阵，每一列代表一个端元的光谱向量；向量 ${{S}}$ 是像元中各端元的丰度； ${{N}}$ 是 $l$ 维的高斯随机噪声。如果把高光谱图像的全部 $c$ (假设 $c$ 为像元个数)个像元均考虑进式(1)，则式(1)可扩展为

${X} = {MS} + {N}$

(2)

式中， ${{X}} \in {{{R}}^{l \times c}}$ ， ${{R}}$ 代表非负矩阵，列向量为每个像元的光谱值； ${{S}} \in {{{R}}^{p \times c}}$ 表示丰度矩阵或称混合矩阵。 ${{S}}$ 存在两个约束条件，一是非负性约束( ${s_{i, j}} \geqslant 0, s = 1, 2, \cdots, p, j = 1, 2, \cdots, c$ )；二是全加性约束( $\sum\limits_{i = 1}^p {{s_{ij}}} = 1, j = 1, 2, \cdots, c$ )。

    (2.2) NMF算法原理

NMF通过把一个高阶非负矩阵分解为两个低阶非负矩阵的乘积，达到信号分解、维数约减和提取数据隐含特征等目的，具体过程为：已知非负矩阵 ${{{X}}_1} \in {{{R}}^{m \times n}}$ ( ${{R}}$ 表示非负有理数)和正整数 $p(p < \min (m, n))$ ，同时寻找元素均为非负的矩阵 ${{{X}}_2} \in {{{R}}^{m \times p}}$ 和 ${{{X}}_3} \in {{{R}}^{p \times n}}$ ，使其满足

${{{X}}_1} \approx {{{X}}_2}{{{X}}_3}$

(3)

或等价表示为

${{{X}}_1}_{_j} \approx {{{X}}_2}{{{X}}_3}_{_j}$

(4)

式中， ${{{X}}_1}_{_j}$ 是 ${{{X}}_1}$ 的第 $j$ 列列向量， ${{{X}}_3}_{_j} \in {{{R}}^p}$ ，其中 $p$ 为矩阵 ${{{X}}_2}$ 的期望秩，一般可由先验信息或 ${{{X}}_1}$ 本身确定。

解NMF问题的常用方法是通过最小化 ${{{X}}_1}$ 和 ${{{X}}_2}{{{X}}_3}$ 之间的欧氏距离，寻求最优分解。

$\min f({{{X}}_2}, {{{X}}_3}) = \frac{1}{2}\left\| {{{{X}}_1} - {{{X}}_2}{{{X}}_3}} \right\|_{{F}}^2$

(5)

式中，符号 $\left\| {\,\,} \right\|_{{F}}^2$ 代表Frobenius模。

$\left\| {{{{X}}_1} - {{{X}}_2}{{{X}}_3}} \right\|_{{F}}^2 = \sum\limits_{i = 1}^m {\sum\limits_{j = 1}^n {{{({{{X}}_1}_{_{i, j}} - {{({{{X}}_2}{{{X}}_3})}_{i, j}})}^2}} } $

(6)

对于目标函数式(6)，一般采用乘法准则进行迭代。

    (2.3) NMF在像元解混中存在的主要问题

对比式(2)和式(3)可知，NMF与LSMM的数学形式相似，且所有变量均设定为非负，但NMF存在着诸如分解维度p确定、初始化、并行计算等难题。具体到高光谱解混应用，NMF的问题主要体现在两方面：一是目标函数因非凸导致的局部极小，二是通用性导致对高光谱图像重要特性的忽略。上述两点将显著损害NMF解混算法的精度和稳定性，必须予以处理。

首先分析局部极小问题。NMF的求解过程可理解为通过迭代求目标函数最小值的过程，理想情况下目标函数(6)应为凸函数。凸函数的一个重要性质是“任何局部极小同时也是全局最小”，如果目标函数为凸且取得局部极小收敛，就说明得到了函数的全局最小值，亦即全局最优解。但Hoyer(2004)的研究表明，式(6)同时对于 ${{{X}}_2}$ 和 ${{{X}}_3}$ 是非凸函数。也就是说，NMF目标函数存在大量局部极小，可能得到关于 ${{{X}}_2}$ 和 ${{{X}}_3}$ 的很多对不唯一解，这就是所谓NMF的“局部极小”问题。该问题会降低NMF解混的收敛速度，增大结果的不确定性，并显著降低解混精度。

接下来分析通用性问题。NMF是一种抽象且宽泛的信号分析技术，将自然界中的任意非负混合信号都视为两个非负矩阵的乘积。由于NMF的通用性，在处理具体问题时容易忽略问题本身的个性特征，导致实际效果不佳。具体到高光谱解混，由于高光谱图像是通过光谱特征和空间特征(或称图像特征)共同描述的，相邻波段的光谱相关特征，以及相邻像元的空间相关特征，均与混合像元联系紧密，对于解混具有重要意义。然而，在NMF解混过程中，只将高光谱图像视为光谱测量记录的无规律数学矩阵，忽略了空间图像属性，导致结果容易出现较大偏差。

因此，需要深入挖掘高光谱图像中蕴含的相关性特征，建立对应的必要约束条件，引导NMF解混过程向符合地物真实情况的方向靠近，消除或弱化上述问题带来的不利影响。

3、结合空间与谱间相关分析的NMF解混

高光谱遥感图像中蕴含丰富的相关特征，其中与解混相关的包括空间相关特征和谱间相关特征等。前者包括端元组成的种类和数量随混合像元空间位置的变化而变化，以及相邻像元取值相近等表现形式；后者包括端元光谱曲线分段平滑，以及相邻波段间存在大量光谱信息冗余等表现形式。

论文提出的NMFSSC算法，针对NMF的通用性问题，在NMF目标函数外部加入基于MRF的相邻像元空间相关约束，通过体现高光谱图像的固有空间特征，修正标准NMF解混过程中地物空间能量分布的畸变，提高解混精度；针对NMF的局部极小问题，在NMF目标函数内部加入基于信号预测度技术的相邻波段谱间相关(光谱分段平滑)约束项，通过体现高光谱图像的固有谱间特征，改善标准NMF解混过程中的目标函数收敛性能，进一步提高解混精度。

下面分别讨论两种相关性约束的相关概念和算法设计。

    (3.1) 基于MRF的空间相关特征约束

在高光谱遥感图像中，地物的空间分布往往呈现两种特征: 一是都有各自的主导区域;二是在地表空间上分布连续(汤毅 等，2014)。上述特点导致相邻像元观测数值(DN值、辐亮度或地表反射率等)之间存在不同程度的近似性和相互依赖性，既是空间相关特征的主要表现形式之一，也是地物蕴含于混合像元及其形成过程的重要物理属性。论文选择以前文提到的MRF模型为载体，在解混算法中显式且定量地描述相邻像元的空间相关特征。

设 ${q} = \{ {q_1}, {q_2}, \cdots, {q_n}\} $ 是随机变量的集簇 ${{Q}} = \{ {{{Q}}_1}, {{{Q}}_2}, \cdots, {{{Q}}_n}\} $ 的一个实现，所有可能的实现组成一个空间 ${{\varPhi }}$ ， ${{N}} = \{ {{{N}}_\theta }, \theta \in {{\varPhi }}\} $ 为定义在 ${{\varPhi }}$ 上的邻域系统，若随机场 ${{Q}}$ 满足条件：

(1) ${{P}}({{Q}} = q) > 0, q \in {{\varPhi }}$

(2) ${{P}}({{{Q}}_\theta } = {q_\theta }|{{{Q}}_\tau } = {q_\tau }, \theta \ne \tau) = {{P}}({{{Q}}_\theta } = {q_\theta }|{{{Q}}_\tau } = {q_\tau },$ $ \tau \in {{{N}}_\theta })$

称 ${{Q}}$ 为以 ${{N}}$ 为邻域系统的MRF，条件(2)被称为MRF的局部特性。

MRF包含“Markov性质”和“随机场”两个要素，可简单解释为具有Markov性质的随机场。其中，Markov性质也称“无后效性”，是指过程(或系统)在时刻t₀所处的状态为已知的条件下，过程在时刻t>t₀所处状态的条件分布，与过程在时刻t₀之前所处的状态无关的特性；随机场是指来自同一个样本空间的随机变量(随机变量之间可能有依赖关系)的集合。2维数字图像可看作随机场，空间相关性质也可类比于Markov性质(仅相邻时刻间的状态相关)，则空间相关特征显著的高光谱遥感图像可用MRF近似表示。

图像的空间能量大小与地物的空间变化频率和幅度成正比。对于空间相关特征显著的高光谱图像，解混结果中空间能量越小，越可能接近地面空间能量分布的真实情况。因此，论文以Deng和Clausi(2004)于MRF的图像分割模型为基础，建立反映图像空间相关特征的能量函数模型如下

$E = {E_{{F}}} + {E_{{R}}}$

(7)

模型将能量函数分为 ${E_{{F}}}$ 和 ${E_{{R}}}$ 两部分，能够反映像元本身以及与相邻像元之间两种特征。

${E_{{F}}}$ 描述像元本身特征，只和观测值分布有关。用F表示从图像( ${{X}}$ =x)中提取出的特征向量，假设观测图像含有高斯噪声，可得能量函数 ${E_{{F}}}$

${E_{{F}}} = \sum\limits_{i = 1}^B {\left({\frac{1}{2}\ln {{({\sigma _i}^k)}^2} + \frac{{{{({x_i} - \mu _i^k)}^2}}}{{2{{({\sigma _i}^k)}^2}}}} \right)} $

(8)

式中，B为波段数， $\mu _i^k$ 和 ${\sigma _i}^k$ 分别表示 ${{F}}$ 中第 $k$ 个特征的均值和标准差。

${E_{{R}}}$ 描述像元与邻域间的相关性。为简化计算，假设模型具有均匀性和各向同性，使得势函数与集簇的位置和方向无关，可得能量函数 ${E_{{R}}}$ 。

${E_{{R}}} = \sum\limits_{i = 1}^c {\sum\limits_{r = 1}^d {[{\psi _r}\delta ({x_i}, {x_{i: + r}})]} } $

(9)

式中，c为图像中像素个数，当 $m = n$ 时， $\delta (m, n) = - 1$ ；否则 $\delta (m, n) = 0$ ； $d = 2$ 表示一阶邻域模型， $d = 4$ 表示二阶邻域模型； $\{ {\psi _1}, {\psi _2}, \cdots, {\psi _c}\} $ 是权重系数。

方便起见，下文中假设所有高光谱图像行列数均相等，行列数不相等的图像可作为多个行列数相等子图像的组合。基于上述假设，下文中所有丰度矩阵 ${{S}}$ 均为方阵，即存在逆矩阵。

设 ${{U}}$ 为分离矩阵，即丰度矩阵 ${{S}}$ 的逆， ${{Y}}$ 为端元矩阵 ${{M}}$ 的估计，由式(3)可得

${{U}} = {{{S}}^{ - 1}}$

(10)

${{Y}} = {{XU}}$

(11)

式中, ${{X}}$ 为高光谱图像矩阵。论文以 ${{U}}$ 为求解对象，通过NMF和MRF约束交替运行的方式计算 ${{U}}$ 。

为简化推导过程，这里仅考虑一阶邻域系统相关性。把所有权重系数 $\psi $ 设置为1，用 ${{W}}$ 表示特征向量的均值矩阵，同时假设特征向量中各个特征的标准差 ${\sigma ^k}$ 相等，结合式(8)和式(9)，可以推出 ${{E}}({{U}})$ 如式(12)。

$\begin{aligned}{{E}}({{U}}) = & {{{E}}_{{R}}}({{U}}) + {{{E}}_{{F}}}({{U}}) = \sum\limits_{i = 1}^P {\sum\limits_{r = 1}^c {[\delta ({x_i}, {x_{i + r}})]} } + \\ & \frac{N}{2}\ln {\sigma ^2} + \frac{1}{{2{\sigma ^2}}}({{{U}}^{{T}}}{{{R}}^{{T}}}{{RU}} - 2{{WRU}} + {{W}}{{{W}}^{{T}}})\end{aligned}$

(12)

能量函数 ${{E}}({{U}})$ 取最小值时，MRF模型描述的“空间能量”最小，空间分布情况最理想，此时式(12)中 ${{E}}({{U}})$ 关于 ${{U}}$ 的偏导数应为零，据此可最终推导出 ${{U}}$ 的解如式(13)。

${{U}} = {({{{R}}^{{T}}} \cdot {{R}})^{ - 1}} \cdot {({{W}} \cdot {{R}})^{{T}}}$

(13)

为减小与已有约束(NMF目标函数的欧氏距离约束或KL散度约束，非负约束、全加约束等)的相互干扰，新的相邻像元空间相关约束采用了与标准NMF目标函数交替运行的独立函数形式，即迭代解混过程由NMF步骤和MRF步骤交替组成。

    (3.2) 基于预测度技术的谱间相关约束

位于NMF标准目标函数外部的空间相关约束，无法影响函数本身性能，要改善其非凸性和收敛性能，必须从函数内部进行改进，即在目标函数内部加入新约束项。论文选择利用相邻波段的谱间相关特征在光谱曲线形状方面的表现形式—分段平滑，作为构建新约束项的依据。

谱间相关性源于传感器的频谱交叠，一般而言，光谱分辨率越高，谱间相关越强。分段平滑特征即光谱曲线在几个特定的光谱范围(光谱范围主要依据可见光、近红外、短波红外、中红外等进行划分，各光谱范围对应的波段随成像光谱仪的不同而变化)内分别具有非常明显的平滑性，在各光谱范围的边界处则平滑特征减弱，存在不同程度的突变。该相关特征反映了谱线的结构性和冗余性，即根据某波长位置的数值可大致预测出附近其他波长位置的数值。真实高光谱遥感数据中，源于谱间相关的分段平滑特征非常显著和普遍。

信号预测度也可称为复杂度映射 (complexity pursuit)，利用信号的可预测性(predictability)度量信号复杂度，通过在混合信号的变换中寻找具有最小复杂度的信号来达到恢复原始信号的目的，其中预测度反映了利用信号本身某些部分预测其他部分的难易程度。

这里用 $F({{M}})$ 表示高光谱图像端元光谱信号 ${{M}}$ 的预测度，其计算公式为

$F({{M}}) = \sum\limits_{p = 1}^P {\ln \frac{{\sum\limits_{l = 1}^L {{{({{{\bar{ M}}}_p} - {{{M}}_{lp}})}^2}} }}{{\sum\limits_{l = 1}^L {{{({{{\tilde{ M}}}_p} - {{{M}}_{lp}})}^2}} }}} = \sum\limits_{p = 1}^P {\ln \frac{{{{{V}}_p}}}{{{{{U}}_p}}}} $

(14)

${{\tilde{ M}}_{lp}} = {\lambda _s}{{\tilde{ M}}_{(l - 1)p}} + (1 - {\lambda _s}){{{M}}_{(l - 1)p}}, 0 \leqslant {\lambda _s} \leqslant 1$

(15)

式中，分子 ${{{V}}_p}$ 表示信号 ${{\bar{ M}}_p}$ 的方差( ${{\bar{ M}}_p}$ 是 ${{{M}}_{lp}}$ 的均值)，它保证信号序列不会变为常数序列；分母 ${{{U}}_p}$ 表示信号的“粗糙度”，信号越粗糙， ${{{U}}_p}$ 越大；相反地，信号越平滑， ${{{U}}_p}$ 越小。实际上， ${{{U}}_p}$ 反映了先前信号序列的短期移动平均值 ${{\tilde{ M}}_{lp}}$ 对 ${{{M}}_{lp}}$ 的预测程度。 ${{\tilde{ M}}_{lp}}$ 是之前信号序列的指数级加权之和， ${\lambda _{{s}}}$ 则表示预测率。

最大化 ${{{V}}_p}/{{{U}}_p}$ 需要同时满足两个条件：

(1) ${{{M}}_{lp}}$ 存在一定的变化范围；

(2) ${{{M}}_{lp}}$ 的值在时间上变化“缓慢”。

由于 ${{{M}}_{lp}}$ 代表高光谱图像的端元光谱信号，绝大多数的端元光谱都具有分段平滑特征，即同时具备变化性和平滑性，很好地满足了上述两个条件。因此，可利用预测度函数定量描述端元光谱的分段平滑特征，作为NMF目标函数的新约束项。加入光谱预测度约束后的NMF目标函数为

$\begin{array}{l}D({{M}}, {{S}}) = {{Euc}}({{M}}, {{S}}) + \frac{\zeta }{2}{\left\| {{\tilde{ M}} - {{M}}} \right\|^2}\\ = {\left\| {{{R}} - {{MS}}} \right\|^2} + \frac{\zeta }{2}{\left\| {{\tilde{ M}} - {{M}}} \right\|^2}\end{array}$

(16)

式中， $\zeta $ 为规则化参数， ${{M}}$ 为 $L \times P$ 维的矩阵，在位置 $\left({l, p} \right)$ 处的取值为 ${{M}} \leftarrow {{M}}({{R}}{{{S}}^{{T}}} + c{\tilde{ M}})./({{MS}}{{{S}}^{{T}}} +$ $ c{{M}})$ 。

按照梯度下降算法构造更新规则，得到更新规则如下

${{M}} \leftarrow {{M}} \cdot ({{R}}{{{S}}^{{T}}} + c{\tilde{ M}})/({{MS}}{{{S}}^{{T}}} + c{{M}})$

(17)

${{S}} \leftarrow {{S}} \cdot {{{M}}^{{T}}}{{R}}/({{{M}}^{{T}}}{{MS}})$

(18)

此外，为满足ASC，S需要被归一化，归一化计算公式如下

${{{S}}_{pb}} \leftarrow \frac{{{{{S}}_{pb}}}}{{\sum\limits_{p = 1}^P {{{{S}}_{pb}}} }}, 1 \leqslant p \leqslant P, 1 \leqslant b \leqslant B$

(19)

    (3.3) NMFSSC算法流程

综上所述，基于MRF的空间相关特征约束，以及基于预测度技术的谱间相关约束，针对NMF解混存在的问题，分别引入了反映高光谱图像特征的不同约束条件，但二者作用目标和侧重点所有不同。前者由于与NMF目标函数独立、交替运行，主要侧重于修正NMF流程在端元空间能量分布方面的偏差，通过改善“通用性问题”，提高解混精度；而后者存在于NMF目标函数内部，在改善“通用性问题”的同时，主要侧重于修正NMF目标函数的非凸性，通过改善“局部极小问题”，确保并加快NMF目标函数的收敛。

用3.2节中包含谱间相关约束的NMF替换3.1节中的NMF，可得NMFSSC算法步骤如下：

(1) 利用基于最小误差的高光谱信号识别法HySIME (Hyperspectral Signal Identification by Minimum Error)估算端元数量P；

(2) 初始化端元矩阵 ${{M}}$ 和丰度矩阵 ${{S}}$ ；

(3) 根据更新规则式(17)和式(18)分别计算端元矩阵 ${{M}}$ 和丰度矩阵 ${{S}}$ 的迭代结果，并利用 ${{M}}$ 和 ${{S}}$ 计算分离矩阵 ${{U}}$ 的迭代结果；

(4) 归一化分离矩阵 ${{U}}$ 的每一列，同时估计像元对应特征向量的均值矩阵 $W$ ，然后利用式(13)计算 ${{U}}$ ；

(5) 重复步骤(3)、(4)，直到各自的停止准则同时满足，得到一个估计的成分，将其转换为矩阵，即获得一个端元分布。继续迭代，直至满足阈值条件，得到端元分布的估计结果。

其中，端元数量估计算法HySime是一种估计高光谱信号子空间的方法，计算过程复杂，但不需任何参数，具有自适应性，估计准确度较高。文献(Bioucas和Nascimento，2008)给出了该算法的原理与实现过程。

算法步骤(3)中NMF迭代的目标函数收敛性证明可参见Lee和Seung(2001)的文献；步骤(3)收敛保证了分离矩阵U的迭代结果收敛，则步骤(4)中对U的每一列进行归一化得到的特征向量均值矩阵W也是收敛的；容易证明，式(13)关于W收敛，因此步骤(4)中的目标函数收敛。由于步骤(3)和步骤(4)共同组成了一个带有初始边界值的并行交替迭代法，该算法可抽象为一种D-N(Dirichlet-Neuman)交替迭代法，吴金彪(2002)对于D-N交替迭代法及其收敛性的分析结果表明，可以证明该算法的收敛性。

需要说明的是，本文算法主要针对空间相关程度较高，或者空间相关程度一般但谱间相关程度较高的高光谱遥感场景；反之也成立。也就是说，只有空间、谱间两种相关特征同时不显著时，才会严重影响算法性能，多数自然场景下的真实高光谱数据一般不会出现上述情况。下文实验结果能够证明上述结论。

4、实验与分析     (4.1) 实验数据1

实验数据包括两部分，首先选择美国Purdue大学遥感实验田的AVIRIS高光谱数据集，其伪彩色图像(R:40、G:25、B:10)如图1，主要的地表覆盖类型包括：农作物(玉米、大豆、小麦等)、自然植被(草地、树林等)和人工建筑(高速公路、房屋、高塔等)。图1行列数均为145，空间分辨率为17 m，原始波段数量为220，波段范围0.4—2.5 μm，光谱分辨率为10 nm。由于水汽吸收或低信噪比等原因，本文将波段1—4，103—111，149—165和214—220共计37个波段的数据去除，利用剩余共183个波段进行解混实验。

该数据可在Purdue大学网站下载(http://cobweb.ecn.purdue.edu/~biehl/MultiSpec/[2016-07-07])，同时可下载的还有该地区的实地调查报告及配套实测数据。实地调查报告把该区域划分为16类，根据式(20)计算16种端元光谱A与B两两间的光谱角距离(SAD)，其均值为0.1461，说明各波段间的谱间相关性较高，冗余信息多。

${{SAD}}({{A}}, {{B}}) = {\cos ^{ - 1}}\left({\frac{{{{{A}}^{{T}}} \times {{B}}}}{{\left\| {{A}} \right\| \times \left\| {{B}} \right\|}}} \right)$

(20)

根据式(21)计算整幅图像的全局Moran’I空间自相关指数，结果为0.4174，说明各像元间的空间相关特征也较为显著。

$I = \frac{{n\sum\limits_{i = 1}^n {\sum\limits_{j = 1}^n {{w_{ij}}} } \left( {{x_i} - \bar x} \right)\left( {{x_j} - \bar x} \right)}}{{\sum\limits_{i = 1}^n {\sum\limits_{j = 1}^n {{w_{ij}}} } \sum\limits_{i = 1}^n {{{\left( {{x_i} - \bar x} \right)}^2}} }}$

(21)

式中，n为样本量，即空间位置的个数。x_i、x_j是空间位置i和j的观察值，w_ij表示空间位置i和j的邻近关系，当i和j为邻近的空间位置时，w_ij=1；反之，w_ij=0。全局Moran指数I的取值范围为[–1, 1]。

上述计算结果表明，实验数据同时具有较显著的空间相关特征和谱间相关特征，这也是大部分真实高光谱遥感图像的共同特点。

对于配套实测数据，利用式(22)分别计算端元光谱之间的相关系数。对光谱相关度高(相关系数在0.8以上)的端元类型进行合并，同时舍弃大多数小目标(即把实际分布面积比例很小的几类地物类型视为噪声，假设它们不存在，所有像元中都只存在下文中6种主要地物)，可得如下6类主要地物：玉米(A)、大豆(B)、干草堆(C)、小麦(D)、树木I和草地(F)。第i波段和第j波段之间的相关系数为

${R_{ij}} = \frac{{\sum\limits_{k = 1}^N {(x_k^{(i)} - {{\bar x}^{(i)}})(x_k^{(j)} - {{\bar x}^{(j)}})} }}{{\sqrt {\sum\limits_{k = 1}^N {{{(x_k^{(i)} - {{\bar x}^{(i)}})}^2} \cdot \sum\limits_{k = 1}^N {{{(x_k^{(j)} - {{\bar x}^{(j)}})}^2}} } } }}$

(22)

式中，N为全部像元个数， $x_k^{(i)}$ 为第i波段数据的第k个像元值， ${\bar x^{(i)}}$ 为第i波段的平均值。

解混实验主要针对上述6类地物进行。

实验中用到的解混算法，除NMFSSC外，还包括前文中提到的分段平滑和稀疏约束的非负矩阵分解(PSNMFSC)、交互投影子梯度的非负矩阵分解(APSNMF)和最小体积约束的非负矩阵分解(MVCNMF)3种算法。其中，MVCNMF的目标函数包括两部分：一是估计观测数据与端元和丰度重建数据之间的近似误差，二是最小体积限制。PSNMFSC、APSNMF和MVCNMF均为目前比较有代表性的、基于NMF的高光谱解混算法。

图2列出了根据实地调查报告和地面实测数据得出的6种主要地物的真实丰度分布图，其中的白色区域表示该类地物的真实空间覆盖范围。

图3、图4、图5和图6分别列出了该数据的MVCNMF、PSNMFSC、APSNMF和NMFSSC盲分解丰度估计结果。论文依据实地调查报告及配套实测数据确定真值，利用式(20)计算SAD代表端元光谱估计精度，利用式(23)计算均方根误差(RMSE)计算丰度估计精度。

${{RMSE}}_{i, j}^k = \sqrt {\frac{1}{L}\sum\limits_{j = 1}^L {{{(\varepsilon _{i, j}^k)}^2}} } $

(23)

式中， $\varepsilon _{i, j}^k = S_{i, j}^k - \hat S_{i, j}^k$ 是第k个端元在图像中位置为(i, j)的像元的丰度余差， ${{RMSE}}_{i, j}^k$ 是第k个端元在图像中位置为(i, j)的像元的均方根误差，求全部像元的均值可得第k个端元的丰度均方根误差，即丰度估计精度。

端元光谱和丰度解混的精度分析结果分别如表1和表2，其中粗体代表每一列的精度最高项。可以看出：4种解混算法均能以较为理想的精度有效分解出6种主要地物类型，其中本文提出的NMFSSC算法综合性能最好，以端元光谱分解结果为准，NMFSSC比其他3种算法的分解精度分别提高了11.4%，12.6%和8.91%；以丰度估计结果为准，NMFSSC比其他3种算法的分解精度分别提高了10.8%，11.1%和7.12%。其余3种算法精度接近，其中APSNMF(交互投影子梯度的非负矩阵分解算法)略好一些，而MVCNMF和PSNMFSC算法精度基本相当。

上述结果表明，NMFSSC算法中空间和谱间相关性特征约束项的联合，对于克服NMF的通用性问题和局部极小问题，进一步提高解混精度具有积极意义和较理想能力。

    (4.2) 实验数据2

上述实验图像的自然场景以耕地里成片分布的自然场景为主，可以预见其空间相关性较高。为进一步验证算法对空间相关程度一般的场景的解混效果，实验数据2选择获取于1995年7月美国Nevada州Cuprite采矿区的AVIRIS高光谱数据(图7)。图7中的A、B、C、K、M等字样分别代表该地区5种广泛分布的矿物Alunite、Buddingtointe、Calcite、Kaolinite和Muscovite的大致分布位置。

该图像大小为400列，350行，空间分辨率为20 m，波长范围为1.99—2.48 μm，光谱分辨率为10 nm，共50波段，依次为AVIRIS原始谱段中的第172—221波段。该地区位于美国Nevada州南部，地表多为裸露矿物，基本无植被覆盖。相对于上一实验数据，该数据的空间相关特征不明显，几种主要矿物空间分布的随机性较大。

由于AVIRIS为机载光谱成像仪，空间分辨率(20 m)相对较高，可合理假设图像中每种端元都含有一定数量的纯像元。因此，本文通过在原图像中人工选择参考点的方式，收集每种端元(地物)的光谱作为参考值，利用全约束最小二乘法计算端元的丰度参考值。

下面重复数据1的实验环节，以此对比和验证算法针对不同场景时的性能表现。需要说明的是，实验只针对上述5种主要矿物类型进行解混实验，对于其余小目标代表的地物类型的假设和处理与数据1完全相同，具体情况不再赘述。

图8为本文算法NMFSSC对应的Cuprite采矿区AVIRIS数据丰度估计结果。由于篇幅限制，MVCNMF等其他3种算法的结果图不再一一列出。

表3和表4分别给出了4种解混算法的端元光谱和丰度解混的精度分析结果。可以看出：4种解混算法均能有效分解出5种主要矿物类型，其中本文提出的NMFSSC算法性能仍然最好，但与其余3种算法均比较接近。以端元光谱分解结果为准，NMFSSC比其他3种算法的分解精度分别提高了5.53%，7.56%和6.03%；以丰度估计结果为准，NMFSSC比其他4种算法的分解精度分别提高了6.18%，6.76%和7.94%。其余3种算法精度基本相当。

上述结果表明，本文算法在实验数据空间相关程度显著减小的情况下，相对于其他几种NMF代表性算法，解混精度有所下降，但仍然具有一定优势。分析其原因，实验数据空间相关程度的降低，会不可避免地导致算法基于空间相关特征的约束部分出现一定程度的性能下降；但由于算法还包括了谱段间的相关性约束，能够在较大程度上弥补空间相关约束部分解混性能降低带来的负面影响，在整体上仍保持较为理想的解混性能。

最后，为了更全面地分析算法解混性能，给出了上述4种算法在matlab7.0环境下分别解混上述两组空间相关程度不同的实验数据时运行时间的统计与比较结果(表5)，单位为秒。其中，为便于和遥感试验田数据进行比较，截取了Cuprite采矿区数据中相同行列数(145行，145列)的子区作为实验数据；波段数也都统一选定为50(从数据1中随机选取50个波段)，两组试验数据的感兴趣端元数量均设定为5。

从表5中可以看出，论文算法NMFSSC的运行时间要少于其他3种算法，运算效率较高；同时该算法在处理空间相关程度较低的数据时耗时有所增大。分析其原因，主要是由于空间相关程度的降低使得算法迭代次数增大。

5、结　论

论文研究结果表明，在解混过程中利用MRF的改进模型与复杂度映射技术，分别引入高光谱图像的空间相关特征(相邻像元灰度近似)与谱间相关特征(端元光谱分段平滑)，能够有效遏制NMF的局部极小问题和通用性问题，引导NMF解混过程及结果不断趋近于真实高光谱图像特征，相对其余几种代表性NMF解混算法，能够进一步改善解混精度。

相对于前人在利用相关性分析改进高光谱解混方面的研究成果，论文主要在如下两个方面进行了改进和创新：一是针对NMF影响解混性能的“局部极小”和“通用性”两大问题，在解混模型中同时对图像的空间与谱间相关性特征进行了描述，较全面地弱化了NMF本身特性对解混的不利影响；二是NMFSSC属于非监督算法，能够降低对于难以获取的先验知识的依赖，相对于部分对先验知识依赖度较高的现有算法，在一定程度上改善了算法实用性。

最后应当指出，尽管相关特征在高光谱图像中普遍存在，但并非全部真实数据都同时具有显著的(相邻像元)空间相关和(光谱平滑)谱间相关特征。对于空间相关与谱间相关特征均不显著的极端情况，论文算法可能出现精度显著降低或性能不稳定等问题。如何进一步拓展算法适用范围，稳定算法性能，将是下一步研究工作的主要内容。