发布时间: 2017-04-16
摘要点击次数:
全文下载次数:
DOI: 10.11834/jig.20170414
2017 | Volume 22 | Number 4

遥感图像处理

结合空间光谱预处理和约束非负矩阵分解的高光谱图像混合像元分解

彭倩¹, 张兵¹, 孙旭¹, 高连如¹, 于文博^1,2

1. 中国科学院遥感与数字地球研究所, 北京 100094;

2. 郑州大学信息工程学院, 郑州 450001

收稿日期: 2016-10-24; 修回日期: 2017-01-03

基金项目: 国家自然科学基金项目（41325004，41571349）

第一作者简介: 彭倩 (1993-), 女, 中国科学院遥感与数字地球研究所信号与信息处理专业硕士研究生, 主要研究方向为高光谱混合像元分解。E-mail:pengqian@radi.ac.cn

中图法分类号: TP751

文献标识码: A

文章编号: 1006-8961(2017)04-0542-09

摘要

目的混合像元问题在高光谱遥感图像处理分析中普遍存在，非负矩阵分解的方法被引入到高光谱图像解混中。本文提出结合空间光谱预处理和约束非负矩阵分解的混合像元分解流程。方法结合空间光谱预处理的约束非负矩阵分解，如最小体积约束、流行约束等，通过加入邻域的空间和光谱信息进行预处理获得更优的预选端元，从而对非负矩阵分解的解混结果进行优化。结果在5组不同信噪比的模拟数据实验中，空间预处理（SPP）和空间光谱预处理（SSPP）均能够有效提高约束非负矩阵分解（最小体积约束的非负矩阵分解和图正则非负矩阵分解）的解混结果，其中SPP在不同信噪比的情况下都能优化约束非负矩阵分解的结果，而SSPP在低信噪比的情况下，预处理效果更佳。利用美国内华达州Cuprite矿区数据进行真实数据实验，SPP提高了约束非负矩阵分解的解混精度，而SSPP在复杂场景下，解混精度更佳。模拟数据和真实数据的实验均表明，空间光谱预处理能够有效地提高约束非负矩阵分解的解混精度，特别是对于信噪比较低的情况下，融合空间和光谱信息对噪声有很好的鲁棒性。结论本文对约束非负矩阵分解的解混算法添加空间光谱预处理，利用高光谱遥感数据的空间和光谱信息，优化预选端元，加入空间光谱预处理的非负矩阵解混实验流程，在复杂场景情况下，对噪声具有较好的鲁棒性。

关键词

高光谱图像; 混合像元分解; 非负矩阵分解; 空间光谱预处理; 光谱混合分析

Hyperspectral unmixing based on spatial and spectral preprocessing prior and constrained non-negative matrix factorization

Peng Qian¹, Zhang Bing¹, Sun Xu¹, Gao Lianru¹, Yu Wenbo^1,2

1. Institute of Remote Sensing and Digital Earth, Chinese Academy of Sciences, Beijing 100094, China;

2. School of Information Engineering, Zhengzhou University, Zhengzhou 450001, China

Supported by: Supported by: National Natural Science Foundation of China (41325004, 41571349)

Abstract

Objective The problem of mixed pixels is common in hyperspectral remote sensing image processing and analysis, and the non-negative matrix factorization (NMF) method has been introduced into hyperspectral unmixing. This study proposes a new spatial and spectral preprocessing technique for constrained NMF, thereby making the unmixing process robust to noise. Method The spatial and spectral preprocessing technique is based on constrained NMF, such as minimum volume constrained NMF and graph-regularized NMF. The unmixing result of constrained NMF is improved by obtaining better endmember candidates from the preprocessing within the spatial and spectral information of the neighborhood. Result Spatial preprocessing and spatial-spectral preprocessing can both improve the unmixing accuracy of constrained NMF, such as minimum volume constrained NMF and graph-regularized NMF, in five groups of simulation data with different signal-to-noise ratios (SNRs). Spatial preprocessing can enhance the results of constrained NMF in all SNRs, whereas spatial-spectral preprocessing effectively optimizes accuracy, particularly for conditions with low SNR. Real data experiment based on the well-known hyperspectral data captured by the Airborne Visible/Infrared Imaging Spectrometer over Cuprite, Nevada is conducted. Spatial preprocessing improves the unmixing results of constrained NMF in the real data set, and spatial-spectral preprocessing is better than spatial preprocessing. The experiments on both simulated data and real data show that spatial and spectral preprocessing can efficiently improve the unmixing accuracy of constrained NMF, particularly in low SNR condition. This finding indicates robustness to noise with spatial and spectral information. Conclusion This study introduces spatial and spectral preprocessing into constrained NMF. The unmixing process becomes more robust to noise with constrained NMF by optimizing the endmember candidate using the spatial and spectral information of hyperspectral remote sensing data in complex remote sensing scenes.

Key words

hyperspectral image; pixel unmixing; non-negative matrix factorization; spatial and spectral preprocessing; spectral mixture analysis

0 引言

高光谱遥感技术同时探测地物的2维空间信息和1维光谱信息，高光谱波段通常连续且多达十个甚至数百个，因此广泛应用到精细农业、地质调查、军事目标识别、环境监测、灾害评估等行业中^[1-2]，受到广泛重视。由于遥感器空间分辨率的限制和地物的复杂多样性，高光谱遥感数据的空间分辨率较低，混合像元成为高光谱遥感数据处理和分析中的重要问题。光谱信息反映地物的理化信息，可以分解混合像元光谱获得其组成端元和比例^[3-5]。

光谱解混算法主要依赖于场景混合模型的建立，现有的模型主要有线性和非线性模型两种^[5]。线性混合模型是迄今为止高光谱图像处理中应用最多的模型。线性模型假设太阳入射辐射只与一种地物表面发生作用，物体间没有相互作用^[2]。

基于线性解混模型的算法主要包括基于纯像元、基于最小体积模型和基于统计模型的方法^[6-7]。在基于统计模型的方法中，将无监督线性高光谱混合像元分解问题视为盲源分离问题 (BSS)。非负矩阵分解 (NMF) 作为一种盲源分离的方法，其分解前后的矩阵都具有非负性，能够满足高光谱数据的端元矩阵和丰度矩阵的非负性要求，因此被运用到混合像元分解问题中^[8-15]。

但是在实际高光谱混合像元分解问题中，仅有非负矩阵分解中的非负约束是远远不够的，为了得到理想的分解值，很多约束条件被加入到分解中，其中用端元来约束非负矩阵分解的有最小体积约束的非负矩阵分解算法 (MVCNMF)，用上下文约束的有流行正则非负矩阵分解算法 (GNMF)^[9-11]等。这些算法都是将原始数据加入到约束非负矩阵分解的系列算法中，通过对算法本身的约束调整获得最优解，对异常端元敏感，在信噪比较低的情况下，解混结果受噪声影响较大。

本文对原始高光谱数据进行融合空间和光谱的预处理，考虑相邻像元间的空间关系和光谱相似性^[16-17]，然后将预处理之后的光谱数据运用到约束非负矩阵分解中，预处理对预选端元进行优化，为非负矩阵分解带来更好的初始值，从而在迭代之后获得更好的解混结果，并且对噪声有更好的鲁棒性。对模拟数据和真实数据的实验均表明，加入空间光谱预处理的约束非负矩阵分解解混实验流程，能够对有效的提高解混精度。

1 预处理方法

1.1 空间预处理

空间预处理 (SPP) 算法，在传统的基于光谱信息的端元提取算法中加入空间信息，根据像元与相邻像元的光谱相似程度提出一个比例因子，作为该像元在空间上下文中光谱信息的重要程度的权重^[16]。

假设高光谱数据中包含纯像元的均质区域，在其中寻找最纯的光谱作为端元。对每个像元，计算它与邻域像元的空间相似性程度，并作为该像元光谱信息的空间权重。以像元${\boldsymbol{R}\left( {i, j} \right)}$为中心建一个$ w \times w$($w $为奇数) 的正方形邻域窗口，定义比例权重

$ \alpha \left( {i, j} \right) = \sum\limits_{r = i-d}^{r = i + d} {\;\sum\limits_{s = j-d}^{s = j + d} {\beta \left( {r-i, s - j} \right)} } \cdot \gamma \left( {r - i, s - j} \right) $

(1)

$ \left\{ \begin{array}{l} d = \left( {w-1} \right)/2\\ \gamma \left( {r-i, s-j} \right) = \gamma \left( {\boldsymbol{R}\left( {r, s} \right), \boldsymbol{R}\left( {i, j} \right)} \right) \end{array} \right. $

(2)

式中，$ \gamma $是中心像元${\boldsymbol{R}\left( {i, j} \right)} $和邻域像元$ {\boldsymbol{R}\left( {r, s} \right)}$的相似性度量，即计算两个像元间的光谱角距离^[2]，$\beta $是相似性度量值的权重，根据窗口内的位置来调节权重，离中心像元越近的位置越重要，权重也越大。

用$ \rho \left( {i, j} \right) = {\left( {1 + \sqrt {\alpha \left( {i, j} \right)} } \right)^2}$定义空间权重因子，按比例缩小原始图像的像元光谱值，即

$ \boldsymbol{R}'\left( {i,j} \right) = \frac{1}{{\rho \left( {i,j} \right)}}\left( {\boldsymbol{R}\left( {i,j} \right) - \boldsymbol{\overline R} } \right) + \boldsymbol{\overline R} $

(3)

式中，$\boldsymbol{R}'\left( {i, j} \right) $是加权后的光谱特征值，${\boldsymbol{\overline R} } $是原始图像所有像元的光谱均值。以上过程获得修正后的单形体，完成空间预处理。

空间预处理根据邻域像元的光谱相似程度所提出的比例因子对原始单形体进行修正，使得单形体顶点更有可能是位于空间均质区域的纯像元，从而在后续的端元提取中更好的获得端元位置，从而为非负矩阵分解提供更好的初始值。

1.2 空间光谱预处理

空间光谱预处理 (SSPP) 算法同时运用了空间和光谱两种信息进行预处理，从而能够获得空间同质且光谱纯度高的像元用于后续的端元提取过程。SSPP算法能够有效提高端元提取的精度，不会增加提取的复杂度，而且能够有效去除噪声点干扰^[17]。

空间光谱预处理主要包括四步：多尺度高斯滤波、计算空间同质性、光谱聚类以及空间与光谱信息融合。设原始高光谱图像为$\boldsymbol{R} $，分别对$ L$个波段分别进行高斯滤波处理并组成新的图像，相当于对图像进行平滑处理。计算多尺度高斯滤波处理前后两个图像，将每个像元多个波段像元值的均方根误差定义为原始图像中每个像元的空间同质性指数，均方根误差越小则说明原始图像中该像元与其邻域中像元的相似度越高。

对原始高光谱图像进行基于光谱信息的聚类，本文采用K-均值聚类。然后在各类别中根据每个像元的空间同质指数排序，按照比例${\rho _1} \in \left( {0, 1} \right) $选择同质指数最小的部分像元。同时对原始高光谱图像进行纯像元指数 (PPI)^[2]计算，即生成单位随机变量$ {\boldsymbol{\varphi }_k} $，并计算图像在${\boldsymbol{\varphi }_k}$方向上的投影$ \boldsymbol{R}{'_k} = \boldsymbol{\varphi }_k^{\rm{T}}\boldsymbol{\boldsymbol{R}}$，找到$\boldsymbol{R}{'_k} $中最大最小的位置，多次操作并更新所有像元的纯像元指数，按照比例$ {\rho _2} \in \left( {0, 1} \right)$选择光谱纯度指数最高的部分像元。然后由空间同质性和光谱纯度指数筛选获得的像元组成候选端元集，此时完成空间光谱预处理过程。

2 解混与非负矩阵分解

2.1 线性光谱混合模型

线性光谱混合模型 (LSMM)，假设太阳入射辐射只与一种地物表面发生作用，物体间没有相互作用。LSMM是应用最广的解混模型^[2]，其矩阵形式为

$ \boldsymbol{\boldsymbol{R = }}\boldsymbol{EA + }\varepsilon $

(4)

式中，$\boldsymbol{R} $表示图像矩阵 (每列一个像元)，$ \boldsymbol{E}$表示端元矩阵 (每列一个端元)，$\boldsymbol{A} $表示丰度矩阵，$\varepsilon $表示误差矩阵^[17]。丰度矩阵中的每个元素$ A\left( {i, j} \right) \ge 0$且每列之和为1。因此高光谱混合像元分解问题就是在“非负”约束与“和为1”约束下，将一个高维矩阵分解为两个低维矩阵的过程。

2.2 非负矩阵分解

NMF的数学模型为：找到两个非负矩阵$ \boldsymbol{Q} \in {{\bf{R}}^{L \times P}}, \boldsymbol{S} \in {{\bf{R}}^{P \times N}}$近似表达已知的非负矩阵$\boldsymbol{X} \in {{\bf{R}}^{L \times N}} $，可表示为^[18]

$ \boldsymbol{X} \approx \boldsymbol{QS} $

(5)

式中，$P < {\rm{min}}\left( {L, N} \right) $。

一种常见的目标函数是基于欧氏距离

$ f\left( {\boldsymbol{Q}, \boldsymbol{S}} \right) = \left\| {\boldsymbol{X}-\boldsymbol{QS}} \right\|_{\rm{F}}^2 = \sum\limits_{i, j = 1}^{L, N} {} {\left( {{\boldsymbol{X}_{ij}}-{{\left( {\boldsymbol{QS}} \right)}_{ij}}} \right)^2} $

(6)

式中，${\left\| {\; \cdot \;} \right\|_{\rm{F}}} $是Frobenius范数。由于目标函数的非凸性，标准NMF算法只能解决目标函数局部最小解问题，因此越来越多的学者在标准NMF的基础上加入一定的约束条件以实现NMF的改进。

2.3 最小体积约束的非负矩阵分解

最小体积约束的非负矩阵分解算法 (MVCNMF) 进行高光谱混合像元分解^[9]，将几何模型中的体积限制加入到NMF算法中，强调光谱数据非负性的同时，要求端元构成的单形体在离散数据空间中体积最小。

MVCNMF算法的目标函数为

$ f\left( {\boldsymbol{E}, \boldsymbol{A}} \right) = \frac{1}{2}\left\| {\boldsymbol{R}-\boldsymbol{EA}} \right\|_{\rm{F}}^2 + \lambda J\left( \boldsymbol{E} \right) $

(7)

式中，$ J\left( \boldsymbol{E} \right)$为单形体体积惩罚函数，$ \lambda $是正则化参数。目标函数由重建数据和原始数据间的均方根误差以及最小体积限制构成。均方根误差使得端元向点云外部移动，最小体积限制又使得端元靠近，所以最小体积限制不必定义为单形体的体积，能起到约束作用即可。

由文献[9]知MVCNMF算法的目标函数式 (7) 可表示

$ \begin{array}{l} \mathop {{\rm{min}}}\limits_{E, A} f\left( {\boldsymbol{E}, \boldsymbol{A}} \right) = \frac{1}{2}\left\| {\boldsymbol{R}- \boldsymbol{EA}} \right\|_{\rm{F}}^2 + \frac{\tau }{2}{\rm{de}}{{\rm{t}}^2}\left( {\boldsymbol{P} + \boldsymbol{Q}{\boldsymbol{U}^{\rm{T}}}\left( {\boldsymbol{E}- \boldsymbol{\overline R} \cdot {\bf{1}}_m^{\rm{T}}} \right)} \right)\\ {\rm{s}}{\rm{.t.}}\;\;\;\boldsymbol{E} \in {\bf{R}}_ + ^{L \times m}, \boldsymbol{A} \in {\bf{R}}_ + ^{m \times n}, {\bf{1}}_m^{\rm{T}}\boldsymbol{A} = {{\bf{1}}_n}, \\ \boldsymbol{P} = \left[\begin{array}{l} {\bf{1}}_m^{\rm{T}}\\ {{\bf{0}}_{\left( {m-1} \right) \times m}} \end{array} \right], \boldsymbol{Q} = \left[\begin{array}{l} {\bf{0}}_{m-1}^{\rm{T}}\\ {\boldsymbol{I}_{m-1}} \end{array} \right] \end{array} $

(8)

式中，$\boldsymbol{U} \in {\boldsymbol{R}^{L \times \left( {m-1} \right)}} $为原图像矩阵$ \boldsymbol{R}$主成分分析 (PCA) 变换后所得矩阵，$ {\boldsymbol{\overline R} } $为像元均值，$\tau = \lambda /\left( {m-1} \right)! $，${{\bf{1}}_m} $是一个$ m$维的列向量且元素都为1，$ {{\bf{0}}_{\left( {m-1} \right) \times m}}$是一个$_{\left( {m-1} \right) \times m} $的0矩阵，${{\bf{0}}_{\left( {m - 1} \right)}}$是一个$\left( {m-1} \right) $维列向量且元素都为0，$ {\boldsymbol{I}_{m-1}}$是一个$_{\left( {m-1} \right) \times \left( {m-1} \right)} $的单位矩阵。

2.4 图正则非负矩阵分解

根据高光谱数据的结构特性，若邻近像元属于同一种类别，则其光谱向量以及对应丰度向量的欧氏距离都很接近。根据流形学习理论，对高光谱数据$\boldsymbol{R} = \left[{{\boldsymbol{R}_1}, {\boldsymbol{R}_2}, \cdot \cdot \cdot, {\boldsymbol{R}_N}} \right] \in {{\bf{R}}^{L \times N}} $的所有像元$\left\{ {{\boldsymbol{R}_N}} \right\}_{n = 1}^N $构造一个结构图$\boldsymbol{W} \in {{\bf{R}}^{N \times N}} $。本文构造$\boldsymbol{W} $使用热核函数

$ {\boldsymbol{W}_{ij}} = {{\rm{e}}^{-\frac{{{{\left\| {{R_i}-{R_j}} \right\|}^2}}}{\sigma }}} $

(9)

当光谱向量${\boldsymbol{R}_i} $和${\boldsymbol{R}_j} $越相似，权值$ {\boldsymbol{W}_{ij}}$越大，丰度分量$ {\boldsymbol{A}_i}$和$ {\boldsymbol{A}_j}$也越相似。

光谱图正则项为

$ \frac{1}{2}\sum\limits_{i, j = 1}^N {{{\left\| {{\boldsymbol{A}_i}-{\boldsymbol{A}_j}} \right\|}^2}{\boldsymbol{W}_{ij}} = } \sum\limits_{i = 1}^N {\boldsymbol{A}_i^{\rm{T}}} {\boldsymbol{A}_i}{\boldsymbol{D}_{ii}}-\sum\limits_{i, j = 1}^N {\boldsymbol{A}_i^{\rm{T}}} {\boldsymbol{A}_j}{\boldsymbol{W}_{ij}} = tr\left( {\boldsymbol{AD}{\boldsymbol{A}^T}} \right)-tr\left( {\boldsymbol{AW}{\boldsymbol{A}^{\rm{T}}}} \right) = tr\left( {\boldsymbol{AL}{\boldsymbol{A}^{\rm{T}}}} \right) $

(10)

式中，${\boldsymbol{D}_{ii}} = \sum\limits_{ j } {{\boldsymbol{W}_{ij}}} $是一个对角矩阵，$\boldsymbol{L} = \boldsymbol{D}-\boldsymbol{W} $。

将光谱图正则加入非负矩阵分解中，得到图正则非负矩阵分解 (GNMF) 的目标函数^[10]为

$ f\left( {\boldsymbol{E}, \boldsymbol{A}} \right) = \frac{1}{2}\left\| {\boldsymbol{R}-\boldsymbol{EA}} \right\|_{\rm{F}}^2 + \lambda {\rm{tr}}\left( {\boldsymbol{AL}{\boldsymbol{A}^{\rm{T}}}} \right) $

(11)

式中，$\lambda $是正则化参数。

3 结合空间光谱预处理和约束非负矩阵分解的高光谱图像解混流程

非负矩阵分解前后矩阵元素都满足非负限制，而其非凸性会导致局部最优问题，各种约束条件的加入使得非负矩阵分解算法能够获得满足要求的最优解。本文通过对高光谱遥感数据的预处理，融入空间和光谱信息，使得在约束非负矩阵分解中有更优的初始化条件，实验表明加入预处理对于约束非负矩阵解混算法的结果有优化作用。预处理与非负矩阵分解算法结合的解混实验方案流程框架如图 1所示。

图 1 结合空间光谱预处理和约束非负矩阵分解的高光谱图像混合像元分解实验方案流程图

Fig. 1 Processing flow chart of hyperspectral unmixing experiment based on spatial and spectral preprocessing prior and constrained non-negative matrix factorization

在获得高光谱图像后，使用虚拟维数 (VD)^[19]或者最小误差信号子空间识别 (HySIME)^[20]算法确认该图像的端元数量。然后分别按照空间预处理 (SPP)、空间光谱预处理 (SSPP) 两种预处理方法对高光谱图像进行预处理获得待选端元集。

在对约束NMF算法进行初始化的时候，本文使用顶点成分分析 (VCA) 算法^[21]作为初始对比算法，所得的初始端元作为初始值，分别使用MVCNMF，GNMF进行处理，同时获得最终的端元光谱和丰度分量，最后进行精度评定。

根据光谱角距离 (SAD) 和原始高光谱图像与重构图像的均方根误差 (RMSE)^[2]来进行精度评价。

光谱角距离 (SAD) 用来表示光谱之间的差异程度，SAD越小，对应光谱差异越小，其表达公式为

$ {\rm{SAD}}\left( {\boldsymbol{a}, \boldsymbol{b}} \right) = arccos\left( {\frac{{{\boldsymbol{a}^{\rm{T}}}\boldsymbol{b}}}{{{{\left\| \boldsymbol{a} \right\|}_2} \cdot {{\left\| \boldsymbol{b} \right\|}_2}}}} \right) $

(12)

式中，$\boldsymbol{a}, \boldsymbol{b} \in {{\bf{R}}^L} $分别为算法提取的端元光谱和参考端元光谱。参考端元一般可以从实测光谱库中获得。

均方根误差 (RMSE) 用来评价两个图像的差异程度，表达公式为

$ {\rm{RMSE}} = \frac{1}{n}\sum\limits_{i = 1}^n {\sqrt {\frac{1}{L}\left\| {{\boldsymbol{R}_i}-\boldsymbol{R}{\prime _i}} \right\|_2^2} } $

(13)

式中，${\boldsymbol{R}_i}, \boldsymbol{R}{\prime _i} $分别是原始图像和重构图像的像元。RMSE越小，说明端元提取和丰度反演的结果越好。

4 实验与分析

4.1 模拟数据实验

4.1.1 实验设计

模拟实验是从USGS的veg_1dry光谱库中选择5种地物光谱曲线，分别记为$ {r_1}, {r_2}, {r_3}, {r_4}, {r_5}, $光谱范围为0.4 μm至2.5 μm，波段数量为206个，像元数量200×200，${r_1}, {r_2}, {r_3}, {r_4} $分布在图像的四角，$ {r_5}$在图像的中心位置，每个端元的丰度分别从端元位置向周围减小，丰度值分布在[0, 1]之间^[22]。

端元的丰度图和光谱曲线如图 2所示, 将端元光谱按照丰度图进行混合，得到无噪声的模拟图像。同时为了检验算法对不同程度的噪声的鲁棒性，按照20 :1、40:1、60:1、80:1、100:1的信噪比在无噪声的模拟图像中加入高斯白噪声获得五组模拟数据。

图 2 模拟端元丰度图和光谱曲线

Fig. 2 Simulation abundance distributions and corresponding reflectance spectra ((a) abundance distributions; (b) reflectance spectra)

4.1.2 实验结果分析

利用所得模拟实验数据，分别进行了空间预处理、空间光谱预处理，并将输出结果后续应用于MVCNMF和GNMF中，最后将实验结果与直接应用VCA、MVCNMF、GNMF对原始数据进行解混的结果进行对比，比较结果如图 3所示，图 3(a)是所提取出来的端元与参考光谱曲线的平均光谱角距离，图 3(b)是重构出来的图像与原图的均方根误差。如图 3(a)所示，在信噪比很低的时候 (SNR=20:1)，预处理之后的结果有显著提高，特别是SSPP同时结合了空间光谱信息，而在信噪比较高的时候 (SNR≥60:1)，噪声较小，不论是否加入预处理均能得到较好解混效果。如表 1、表 2中，将SPP或SSPP处理之后的约束非负矩阵分解算法与未经预处理之前的算法结果对比。对于信噪比很低的情况下，很明显两种预处理都能够有效的提高解混精度，对于信噪比较高的情况下，SPP依然能够优化结果，SSPP处理之后的精度虽没有得到提高，但也能达到直接处理所获得的解混精度。针对高信噪比下的这种情况，在实验中通过调整SSPP中对空间同质指数和光谱纯度指数设置的阈值${\rho _1}, {\rho _2} $发现，在信噪比低的时候，选择较大的阈值筛选出更小比例更高纯度的像元，能够达到更好的实验效果；相反，在信噪比高时，选择较小的阈值，获得更大比例的待选像元，能达到更佳效果。对阈值调整的实验分析得出，可能是在高信噪比的情况下，原本的端元提取算法就能达到很好的效果，通过阈值选择反而在初选时去掉了部分可能的端元。图 3(b)也支持了以上论点，在噪声干扰较大的情况下进行预处理，能够有效地提高解混结果。比较SPP和SSPP的实验效果发现，在噪声较高的时候，加入光谱信息的SSPP能够达到更好效果，而在噪声较低的时候，更适合加入邻域空间信息的SPP。表 1、表 2是对应图 3的解混结果。

图 3 不同信噪比下预处理前后所得结果的平均光谱角距离和均方根误差

Fig. 3 The mean SAD and RMSE compare of preprocessing NMF and original NMF under different SNR ((a) mean SAD; (b) RMSE)

表 1 不同信噪比下预处理后的NMF解混结果与无预处理结果对比平均光谱角距离
Table 1 The mean SAD compare of preprocessing NMF and original NMF under different SNR

下载CSV

/rad
SNR	VCA	MVCNMF	GNMF	SPP-MVCNMF	SPP-GNMF	SSPP-MVCNMF	SSPP-GNMF
20	0.254	0.257	0.288	0.229	0.217	0.162	0.164
40	0.031	0.029	0.030	0.025	0.038	0.029	0.030
60	0.023	0.023	0.022	0.024	0.023	0.022	0.022
80	0.025	0.025	0.025	0.016	0.016	0.025	0.029
100	0.025	0.025	0.025	0.015	0.016	0.025	0.034
注：加粗字体为预处理之后的结果更优或者能够达到未经预处理所获得的解混精度。

表 2 不同信噪比下预处理后的NMF解混结果与无预处理结果对比均方根误差
Table 2 The RMSE compare of preprocessing NMF and original NMF under different SNR

下载CSV

SNR	VCA	MVCNMF	GNMF	SPP-MVCNMF	SPP-GNMF	SSPP-MVCNMF	SSPP-GNMF
20	5.95×10^-2	5.93×10^-2	5.95×10^-2	2.94×10^-2	2.91×10^-2	5.94×10^-2	5.94×10^-2
40	1.08×10^-2	7.67×10^-3	1.07×10^-2	7.38×10^-3	7.38×10^-3	1.07×10^-2	1.07×10^-2
60	1.13×10^-3	1.12×10^-3	3.46×10^-3	9.76×10^-4	1.04×10^-3	1.28×10^-3	1.25×10^-3
80	1.13×10^-4	1.12×10^-4	3.02×10^-4	1.51×10^-4	2.13×10^-4	1.28×10^-4	8.28×10^-4
100	1.13×10^-5	1.12×10^-5	2.30×10^-4	9.40×10^-5	1.60×10^-4	2.64×10^-3	1.12×10^-3
注：加粗字体为预处理之后的结果更优或者能够达到未经预处理所获得的解混精度。

4.2 真实数据实验

本文使用的真实数据是美国内华达州Cuprite地区1997年成像的0.4 μm到2.5 μm波段的AVIRIS (airborne visible infrared imaging spectrometer)^[23]数据。截取350×350个像元的图像 (如图 4) 进行实验，去除低信噪比和水汽吸收波段 (1~3、105~115、150~170)，剩余189个波段，并将原始数据乘10^－4并归一化到0, 1。这个区域中有多种矿物质以及较为典型的纯像元存在，在高光谱图像混合像元分解的研究中广泛应用，且其中矿物质都在USGS (united states geological survey) 光谱库中。综合VD和HySime进行端元估计的结果，将端元数量设置为$p = 19 $^[17]，将实验所得端元与USGS光谱库中的光谱进行对比，找到对应最小光谱角距离的光谱，去除重复提取的端元获得8个端元，在这个区域中主要包括明矾石 (Alunite)、高岭石 (Kaolinite)、玉髓 (Chalcedony)、白云母 (Muscovite)、蒙脱石 (Montmorillonite)、黄钾铁矾 (Jarosite)、方解石 (Calcite)，还有少量的水铵长石 (Buddingtonite) 等其他矿物^[2]。表 3中展示的是提取出的端元与USGS光谱库中对比所得最小的光谱角以及对应的矿物。表 3中标粗数值是SPP、SSPP预处理后与所对应的原始处理解混结果相比得到了优化的情况，可以看出在真实数据中，场景比较复杂，加入SPP、SSPP的非负矩阵分解MVCNMF和GNMF解混结果都有提高，其平均光谱角距离和均方根误差的精度更高。在复杂场景下，加入SSPP预处理的解混结果明显优势较强，SPP预处理的解混结果虽在数值上略有优势，但是有部分端元没有提取出来，而SSPP综合了空间光谱特征，相比SPP只考虑空间信息，解混精度最佳。

图 4 Cuprite地区AVIRIS高光谱数据 (假彩色)

Fig. 4 AVIRIS hyperspectral image at Cuprite (false color composition)

表 3 Cuprite数据预处理后的NMF解混结果与无预处理结果对比平均光谱角距离和均方根误差
Table 3 The mean SAD and RMSE compare of preprocessing NMF and original NMF on Cuprite data

下载CSV

Cuprite	光谱角距离/rad									均方根误差
Cuprite	明矾石	水铵长石	白云母	方解石	蒙脱石	高岭石	黄钾铁矾	玉髓	平均	均方根误差
VCA	0.10	0.07	0.14	0.11	0.17	0.19	0.15	0.11	0.13	0.012 428
MVCNMF	-	0.09	0.15	0.12	0.19	0.21	0.15	0.12	0.15	0.007 776
GNMF	0.10	0.16	0.16	0.10	0.16	0.14	0.20	0.11	0.14	0.006 618
SPP-MVCNMF	-	0.09	0.10	-	0.15	0.12	0.11	0.06	0.11	0.004 088
SPP-GNMF	-	0.09	0.10	-	0.14	-	0.14	0.07	0.11	0.004 357
SSPP-MVCNMF	0.09	0.11	0.11	0.09	0.11	0.12	0.17	0.09	0.11	0.006 133
SSPP-GNMF	0.07	0.11	0.11	0.10	0.12	0.12	0.16	0.07	0.11	0.006 632
注：加粗字体为预处理之后的结果更优或者能够达到未经预处理所获得的解混精度, 无法提取端元的结果用“-”表示。

5 结论

通过结合空间光谱预处理和约束非负矩阵分解来进行高光谱图像混合像元分解，对非负矩阵分解中的初始端元进行优化，能够有效地提高解混精度。特别是在处理低信噪比的数据时，预处理之后的算法更具有优势，而且考虑邻域光谱信息比只考虑空间信息更具优势。在高信噪比数据中，预处理之后同样能够达到高精度的解混结果，且只考虑空间信息能够在信噪比变高的情况下，一直保持较好的实验结果。通过对模拟数据和真实数据的实验，利用光谱角距离和均方根误差评价指标进行比较，均证明了本文所提出的结合空间光谱预处理和约束非负矩阵分解的解混实验流程，对于解混结果有一定的改进。其中无论信噪比高低，SPP-MVCNMF和SPP-GNMF均能获得较好的实验结果，而SSPP-MVCNMF和SSPP-GNMF更适合于信噪比较低的情况。

但是在本文所做的空间预处理只是简单获得像元邻域信息，适合有较大均质区域的图像，没有考虑存在较小地物，所以在较大噪声情况下，空间预处理 (SPP) 效果不佳，这也是下一步需要研究的内容。

参考文献

[1] Tong Q X, Zhang B, Zheng L F. Hyperspectral Remote Sensing:Principle, Technology and Application[M]. Beijing, China: Higher Education Press, 2006: 1-3. [ 童庆禧, 张兵, 郑兰芬. 高光谱遥感:原理、技术与应用[M]. 北京: 高等教育出版社, 2006: 1-3.]

[2] Zhang B, Sun X. Hyperspectral Image Spectral Unmixing[M]. Beijing: Science Press, 2015: 2-20. [ 张兵, 孙旭. 高光谱图像混合像元分解[M]. 北京: 科学出版社, 2015: 2-20.]

[3] Keshava N, Mustard J F. Spectral unmixing[J]. IEEE Signal Processing Magazine, 2002, 19(1): 44–57. [DOI:10.1109/79.974727]

[4] Plaza A, Martín G, Plaza J, et al. Recent developments in endmember extraction and spectral unmixing[M]//Prasad S, Bruce L M, Chanussot J. Optical Remote Sensing:Advances in Signal Processing and Exploitation Techniques. Berlin Heidelberg:Springer, 2011:235-267.[DOI:10.1007/978-3-642-14212-3_12]

[5] Plaza A, Benediktsson J A, Boardman J W, et al. Recent advances in techniques for hyperspectral image processing[J]. Remote Sensing of Environment, 2009, 113(S1): S110–S122. [DOI:10.1016/j.rse.2007.07.028]

[6] Bioucas-Dias J M, Plaza A, Dobigeon N, et al. Hyperspectral unmixing overview:geometrical, statistical, and sparse regression-based approaches[J]. IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing, 2012, 5(2): 354–379. [DOI:10.1109/JSTARS.2012.2194696]

[7] Parente M, Plaza A. Survey of geometric and statistical unmixing algorithms for hyperspectral images[C]//Proceedings of the 2nd Workshop on Hyperspectral Image and Signal Processing:Evolution in Remote Sensing. Reykjavik, Iceland:IEEE, 2010:1-4.[DOI:10.1109/WHISPERS.2010.5594929]

[8] Jia S, Qian Y T. Constrained nonnegative matrix factorization for hyperspectral unmixing[J]. IEEE Transactions on Geoscience and Remote Sensing, 2009, 47(1): 161–173. [DOI:10.1109/TGRS.2008.2002882]

[9] Miao L D, Qi H R. Endmember extraction from highly mixed data using minimum volume constrained nonnegative matrix factorization[J]. IEEE Transactions on Geoscience and Remote Sensing, 2007, 45(3): 765–777. [DOI:10.1109/TGRS.2006.888466]

[10] Cai D, He X F, Han J W, et al. Graph regularized nonnegative matrix factorization for data representation[J]. IEEE Transactions on Pattern Analysis and Machine Intelligence, 2011, 33(8): 1548–1560. [DOI:10.1109/TPAMI.2010.231]

[11] Lu X Q, Wu H, Yuan Y, et al. Manifold regularized sparse NMF for hyperspectral unmixing[J]. IEEE Transactions on Geoscience and Remote Sensing, 2013, 51(5): 2815–2826. [DOI:10.1109/TGRS.2012.2213825]

[12] Wang N, Du B, Zhang L P. An endmember dissimilarity constrained non-negative matrix factorization method for hyperspectral unmixing[J]. IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing, 2013, 6(2): 554–569. [DOI:10.1109/JSTARS.2013.2242255]

[13] Qian Y T, Jia S, Zhou J, et al. L_1/2 Sparsity constrained nonnegative matrix factorization for hyperspectral unmixing[C]//Proceedings of 2010 International Conference on Digital Image Computing:Techniques and Applications. Sydney, Australia:IEEE, 2010:447-453.[DOI:10.1109/DICTA.2010.82]

[14] Huang C H. Researches on hyperspectral unmixing based on nonnegative matrix factorization[D]. Xi'an:Xidian University, 2014. [黄春海. 基于非负矩阵分解的高光谱图像解混技术研究[D]. 西安: 西安电子科技大学, 2014.]

[15] Yu X L. The theory of non negative matrix factorization and its application in the highspectrual unmixing[D]. Chengdu:Chengdu University of Technology, 2015. [余肖玲. 非负矩阵分解理论及其在高光谱解混中的应用[D]. 成都: 成都理工大学, 2015.]

[16] Zortea M, Plaza A. Spatial preprocessing for endmember extraction[J]. IEEE Transactions on Geoscience and Remote Sensing, 2009, 47(8): 2679–2693. [DOI:10.1109/TGRS.2009.2014945]

[17] Martin G, Plaza A. Spatial-spectral preprocessing prior to endmember identification and unmixing of remotely sensed hyperspectral data[J]. IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing, 2012, 5(2): 380–395. [DOI:10.1109/JSTARS.2012.2192472]

[18] Pauca V P, Piper J, Plemmons R J. Nonnegative matrix factorization for spectral data analysis[J]. Linear Algebra and Its Applications, 2006, 416(1): 29–47. [DOI:10.1016/j.laa.2005.06.025]

[19] Chang C I, Du Q. Estimation of number of spectrally distinct signal sources in hyperspectral imagery[J]. IEEE Transactions on Geoscience and Remote Sensing, 2004, 42(3): 608–619. [DOI:10.1109/TGRS.2003.819189]

[20] Bioucas-Dias J M, Nascimento J M P. Hyperspectral subspace identification[J]. IEEE Transactions on Geoscience and Remote Sensing, 2008, 46(8): 2435–2445.

[21] Nascimento J M P, Dias J M B. Vertex component analysis:a fast algorithm to unmix hyperspectral data[J]. IEEE Transactions on Geoscience and Remote Sensing, 2005, 43(4): 898–910. [DOI:10.1109/TGRS.2005.844293]

[22] Plaza A, Martínez P, Pérez R, et al. A quantitative and comparative analysis of endmember extraction algorithms from hyperspectral data[J]. IEEE Transactions on Geoscience and Remote Sensing, 2004, 42(3): 650–663. [DOI:10.1109/TGRS.2003.820314]

[23] Green R O, Eastwood M L, Sarture C M, et al. Imaging spectroscopy and the airborne visible/infrared imaging spectrometer (AVIRIS)[J]. Remote Sensing of Environment, 1998, 65(3): 227–248. [DOI:10.1016/S0034-4257(98)00064-9]