高光谱解混（二）——线性混合模型（section Ⅱ）和信号子空间识别技术（section Ⅲ）

线性混合模型

基础模型

 忽略不同端元之间的多次散射且利用丰度系数对地表进行划分时，每个像元的光谱可以由相关丰度系数加权的端元光谱的线性混合结果逼近表示。在这种情况下，给定像元在$B$个波段下的线性混合模型（LMM）可以表示为
$$y_i=\sum _{j=1}^p{\rho }_{ij}{\alpha }_j+w_i$$
其中，${\rho }_{ij}\ge 0$代表端元$j$（ j ∈ { 1 , . . . , p } j\in\left\{1,...,p\right\} j∈{1,...,p}）在第$i$个波段上的光谱测量值，${\alpha }_j\ge 0$代表端元$j$的丰度系数，$w_i$代表加性干扰（噪声和模型误差），$p$代表端元总数。对于给定像素，丰度系数${\alpha }_j$代表第$j$个端元在地表所占的比例区域。因此，它应满足如下两个约束

其中，丰度向量$\alpha =[{\alpha }_1,{\alpha }_2,...,{\alpha }_p{]}^T$位于标准$(p-1)$-单纯形或者概率$(p-1)$-单纯形中。在高光谱解混术语中，上述两个约束分别被称为丰度非负性约束(ANC)和丰度和为一约束(ASC)。

 假设$y$是一个B维（波段数）的列向量，$m_j=[{\rho }_{1j},{\rho }_{2j},...,{\rho }_{Bj}]$代表第$j$个端元的光谱特征。那么将LMM模型以向量模式可以表示为
$$y=M\alpha +w$$其中，$y$代表高光谱图像中一个像元在$B$个光谱波段上的混合光谱特征向量，$M=[m_1,m_2,...,m_p]$是覆盖区域上包含的端元的混合光谱特征矩阵（每一列代表一个端元在该覆盖区域上的光谱特征向量。个人觉得：一个高光谱图像所有的像元对应一个混合光谱矩阵，但是丰度向量是和每一个混合像元一一对应的），$w$也是一个B维的列向量。
 如果$M$的列是仿射独立(线性无关的基础上再附加组合系数和为1这一条件就构成仿射独立)的，则有单纯形C

即，M的列向量是分布在${\mathbb{R}}^B$空间内的点，所有的点（列向量）构成的凸包可以认为是${\mathbb{R}}^B$空间内的一个(p-1)单纯形。（因为对很多数学概念不清楚，所以此处我暂时将这句话理解为：对于在B个波段上具有p个端元的高光谱图像数据，如果高光谱特征混合矩阵M的列是仿射独立的，那么就可以构成(p-1)-单纯形C，构成这个凸包的几个点就是端元）

 下图演示了一个含有3个端元的虚拟混合矩阵M构成的2-单纯形$C$。绿色的点表示谱向量，红色的点是单纯形的顶点并且对应于端元。注意，混合矩阵M的推理等价于单纯形C顶点的识别任务。（个人初步理解该部分是将解混问题转化成了凸包求解的问题：凸包问题是给定点集，求构成凸包的点。解混问题是给定混合矩阵，求构成单纯性的顶点）这一几何观点促进了许多解混算法的诞生，具体在后文详细讲述。

 假设噪声可以忽略不计，下图是包含三个数据集数据的最小体积单纯形概念的说明。根据$y=M\alpha$生成的谱向量位于单纯形内部（绿色点），单纯形的顶点对应端元（红色点）。三个数据集有不同的特点：

• 最左侧的数据集包含纯像元。即对于p个端元，都至少有一个像元只包含该某一端元对应的材料。
• 中间的数据集不含纯像元。但在单纯形的每个面（边）上至少包含p-1个谱向量。

 从图中可以看出，在左侧和中间两个数据集中，可以通过将最小体积（MV）（最小面积）单纯形拟合到数据来完成端元估计，这可以作为几种基于几何的解混算法的理论基础 。而对于最右侧的数据集，MV单纯形（虚线）比真实单纯形小。 这种情况对应于一个高度混合的数据集，其中单纯形面（边）附近没有光谱向量。对于这类问题，我们通常求助于统计框架，其中混合矩阵和丰度系数的估计通过对所涉及的变量和参数（即丰度系数和混合矩阵）采用合适的概率模型表述为一个统计推断问题。

高光谱解混逆问题的特征

 给定包含n个像元的B维（B个波段）向量的数据集$Y=[y_1,y_2,...y_n]\in {\mathbb{R}}^{B\times n}$，线性高光谱解混问题可以描述为估计像元混合矩阵M和丰度向量系数${\alpha }_i,i=1,..,n$。这是一个困难的逆问题，因为光谱特征往往是强相关的，特别是对于产生条件差的混合矩阵。即，高光谱解混估计问题对噪声高度敏感。
 如下图所示，端元$m_2,m_3$非常接近，这样就得到了一个条件恶劣的矩阵。图中将噪声的影响用不确定区域表示。

 为了表征线性解混逆问题，我们使用信噪比（SNR）

其中${\mathrm{R}}_x,{\mathrm{R}}_w$分别代表信号（$x\equiv M\alpha$）和噪声的相关矩阵，$\mathbb{E}$代表期望值。

 除了信噪比之外，还引入了信噪比谱分布（SNR-SD）

其中$({\lambda }_{i,x},e_{i,x})$是根据${\lambda }_{i,x}$值对${\mathrm{R}}_x$进行递减排序后的第$i$个特征向量。$\mathrm{S}\mathrm{N}\mathrm{R}-\mathrm{S}\mathrm{D}(i)$沿信号方向$e_{i,x}$产生SNR。因此为了获取可接受的解混结果，需要保证$\mathrm{S}\mathrm{N}\mathrm{R}-\mathrm{S}\mathrm{D}(i)\gg 1，i=1,...,p$。否则，信号子空间中会存在明显被噪声破坏的方向。（这一段不是很懂）

 下图绘制了$i=1,..,50$区间内针对以下三个数据集的$\mathrm{S}\mathrm{N}\mathrm{R}-\mathrm{S}\mathrm{D}(i)$。

• SudP5SNR40：模拟数据集。从区间[0,1]内均匀分布的随机变量中采样生成混合矩阵M；$p=5,n=5000$; 丰度系数均分分布在概率4-单纯形上;SNR=40dB。
• SusgsP5SNR40：模拟数据集；从USGS光谱库中采样构成混合矩阵M；$p=5,n=5000$; 丰度系数均分分布在概率4-单纯形上;SNR=40dB。
• Rcuprite：真实数据集；著名的AVIRIS铜质数据立方体的子集，大小为250行、191列、188个波段(去除噪声波段)。

 信号和噪声的相关矩阵通过HySime分发的算法和代码得到。从下图可以看出，对于SudP5SNR40，$i<5$时$\mathrm{S}\mathrm{N}\mathrm{R}-\mathrm{S}\mathrm{D}(i)\gg 1$；$i>5$时，$\mathrm{S}\mathrm{N}\mathrm{R}-\mathrm{S}\mathrm{D}(i)\ll 1$，这表明信号子空间的信噪比较高。对于SusgsP5SNR40，由于USGS光谱特征的高相关性，混合矩阵奇异值衰减更快。然而，“大图”与SudP5SNR40数据集的“大图”是相似的。Rcuprite数据集产生了更困难的逆问题，因为$\mathrm{S}\mathrm{N}\mathrm{R}-\mathrm{S}\mathrm{D}(i)$“接近凸形”的值接近1的速度较为迟缓。这是一个条件恶劣的逆问题的明显迹象。

信号子空间识别（降维）

 该部分主要是对高光谱图像采用特征提取（线性空间变换）的方法进行降维，从而最大限度保留信号和压缩噪声，以便后续处理。关于主成分分析PCA和最大噪声分数变换的具体实现步骤可以参看这篇博文（博文中还对光谱解混和观光谱分类进行了概述，是一篇很好的文章）。对于向量空间和子空间的理解参看博文1和博文2。（个人理解：（1）子空间必须过零向量（2）子空间的维度小于等于原始空间的维度）

概述

 给定场景中存在的端元数量通常远小于波段 B 的数量。因此，假设线性模型是对高光谱图像数据的一个很好的近似值，谱向量位于或非常接近低维线性子空间。 子空间的识别能够实现频谱向量的低维准确的表示，投影到信号子空间需要大量的计算和存储空间，但同时也会带来信噪比增益。 对信号子空间中表示的数据进行操作通常是有利的，有时是必要的。 因此，需要将信号子空间识别算法作为混合像元分解前的预处理步骤。

 无监督子空间的识别方法有很多。波段选择或波段提取，顾名思义，就是利用相邻波段之间存在的高相关性，在信噪比较高的波段中选择少数几个光谱成分。投影技术通过优化目标函数来寻找表示数据的最佳子空间。例如：

• 主成分分析(PCA)也称为K.L变换，它通过优化最小均方误差得到的最佳正交线性变换。首先求出原始图像的协方差矩阵，然后求出协方差矩阵的特征值，根据特征值进行排序，求解出相应特征值对应的特征向量；使用由前几个特征向量组成的正交变换矩阵，对原图像进行线性变换，然后得到新的图像。通过主成分变换，可以在保留原图像大部分信息的前提下对原图像数据进行了压缩，提取图像特征。
• 奇异值分解(SVD)使功率最大化。与PCA方法类似：已知子空间维数为p，如果噪声为加性白噪声，则在经验相关矩阵的前 p 个特征向量上的投影就是信号子空间的最大化似然。奇异值分解特别适合于波段高度相关的高光谱数据，此时PCA有可能得到很差的结果。
• 最大噪声分量 (MNF) 和噪声调节主成分 (NAPC)最小化噪声功率与信号功率的比率从而去除噪声并完成降维。NAPC在数学上等同于MNF并且可以解释为两个主成分变换的序列：第一个适用于噪声，第二个适用于变换后的数据集。在白噪声的情况下，这两个指标是等效的。这两种变换可以确定图像的固有特征分布范围，分离与平衡高光谱影像数据中的噪声，减少后续的计算复杂度。经过MNF变换后，可以实现噪声的去相关与标准化。特征值大小可以直接反映了输出分量方差的大小，当特征值越接近1的时候，表明图像上噪声越多，当特征值为1的时候，表明图像上只有噪声存在。可以对每个分量的特征值进行排序输出，大于1的特征值的数量即为有效端元个数。

 光学实时自适应光谱识别系统(ORASIS)框架由美国海军研究实验室开发，旨在进行实时数据处理，已被用于降维和端元提取。这个框架由几个模块组成。在这些模块中，通过识别场景中的一个能够表达场景中的可变性的范例像元子集来实现降维。从场景中收集的每个新像元与每个范例像元之间利用角度进行度量比较。如果新像元与现有的每个范例像元之间有充分的不一致，则将其添加到范例子集。使用改进的Gram-Schmidt方法，可以利用当前的范例子集周期性地创建一个正交基。

 信号子空间的识别是一个模型阶次推理问题，出现了最小描述长度(MDL)或Akaike信息准则(AIC)等信息论准则。这些标准实际上已被用于高光谱应用，采用了Wax和Kailath在[97]中提出的方法。反过来，Harsanyi、Farrand和Chang提出了一种基于奈曼一皮尔逊引理检测理论的阈值方法(HFC)来确定高光谱数据中光谱端元的数量，在[96]中称为虚拟维数(VD)。HFC方法基于一个建立在样本相关和协方差矩阵的特征值上的检测器。一种被称为噪声白化HFC (NWHFC)的修正版本，包括一个噪声白化步骤。HySime (hyperspectral signal identification by minimum error)采用一种基于最小均方误差的方法来推断信号子空间。该方法是基于特征分解的，无监督的，全自动的(即，它不依赖于任何调谐参数)。它首先估计信号和噪声的相关矩阵，然后选择在最小二乘误差含义上最能代表信号子空间的特征值子集。当频谱混合是非线性的时候，线性情况下的低维子空间经常被低维流形所代替，这是拓扑数学主题中定义的一个概念。

信号子空间投影

 假设信号子空间S已经用上面提到的方法之一识别出来，并且$\mathrm{E}\equiv [e_1,e_2,...,e_p]$的列是S的一个标准正交基，其中$e_i\in {\mathbb{R}}^B,i=1,...,p$。一个光谱向量y在S上的关于基E的正交投影的坐标定义为$y_S={\mathrm{E}}^{T}y\in {\mathbb{R}}^p$。将$y$替换为LMM模型后可以表示为

如前所述，投影到信号子空间需要大量的计算和存储空间，但同时也会带来信噪比增益。前两个指标是因为在大多数情况下$p\ll B$。假设噪声w均值为0，协方差为${\sigma }^2\mathrm{I}$，则投影噪声项${\mathrm{E}}^{T}w$的平均功率为$\mathbb{E}||{\mathrm{E}}^{T}w|{|}^2={\sigma }^{2}p$ ($\mathbb{E}$为均值)。易看出，经过投影所带来的噪声功率的相对衰减为$p/B$。噪声能量得到了衰减，因此信噪比有一定的增益。

 下图说明了将数据集投影到信号子空间上的优点。用HySime估计噪声和信号子空间。左边的图显示了来自模拟数据集SusgsP5SNR30的噪声和相应的投影光谱。子空间维数被正确识别（p=5）。投影数据集的信噪比为46.6dB，比噪声数据集的信噪比高16.6dB($\approx (B/p)$)。右边的图显示了来自Rcuprite数据集的噪声和相应的投影光谱。被识别的子空间维数为18。预测数据集的信噪比为47.5dB，比噪声数据集高出5dB。加性噪声的有色性质解释了$(B/p)dB-5dB\simeq 5dB$的差异。

 注意，尽管将数据集投影到信号子空间通常会去除很大一部分噪声，但这并不能改善HU逆问题的条件，因为这种投影不会改变信号子空间特征分量的$\mathrm{S}\mathrm{N}\mathrm{R}-\mathrm{S}\mathrm{D}(i)$值 。

 进一步降低信号子空间噪声的一种可能的方法是利用光谱和空间上下文信息。我们在空间域中给出了一个简短的说明。下图左上侧显示了Rcuprite数据集上第18号特征图像（即在$e_i^T{y}_s,i=18$的情况下获得的图像）。信号子空间的基利用HySime算法得到。右上角显示了使用BM3D的滤波版本。在这个例子中，去噪算法是非常有效的。正如左下角图像所示，在噪声估计中(噪声和去噪图像之间的差异)没有结构。这种有效性也可以从右下角图中显示的有噪声(蓝点)和去噪(绿点)特征图像17和18的散点图中看出。去噪后的图像所对应的散点图更加密集，反映了更低的方差。

仿射集投影

 从现在开始，我们假设观测到的数据集已经投影到信号子空间上，为了简单起见，我们仍然将投影向量表示为$y=\mathrm{M}\alpha +w$。其中$y,w\in {\mathbb{R}}^p,\mathrm{M}\in {\mathbb{R}}^{p\times p}$（p是降维后的维度）。由于$\mathrm{M}$的列属于向量子空间，原始的混合矩阵就可以表示为$\mathrm{E}\mathrm{M}$（$\mathrm{E}\equiv [e_1,e_2,...,e_p]$的列是子空间S的一个标准正交基，其中$e_i\in {\mathbb{R}}^B,i=1,...,p$）。
 在一些解混算法中，已经研究并考虑了光谱特征的可变性，包括所有将端元视为分布的统计算法。这种变异性是基于振幅的，因此主要以光谱形状不变性为特征。也就是说，虽然端元的光谱形状相当一致，但它们的振幅是可变的。这意味着端元特征受到一个正的比例因子的影响，该因子因像元而异。因此，端元光谱特征矩阵不是一整个场景的一个端元光谱矩阵，而是每个像元的端元光谱矩阵

在这种情况且在没有噪声的情况下，观测到的谱向量不再位于一个由固定端元集定义的单纯形中，而是在这个集合中

具体如下图所示：观测数据在超平面上的投影：（a）在超平面上的正交投影（投影向量发生旋转）；（b） 透视投影（缩放$y/(y^{T}u)$将它们带到由$y^{\prime T}u=1$定义的超平面）

 因此，端元光谱的系数$s_i^t\alpha$不必和为一，尽管它们仍然是非负的。需要对数据进行转换以提高模型与现实的匹配度。 如果可以找到从辐射单位到反射率的真实映射，那么这种转换就足够了。 然而，估计映射可能是困难的问题或不可能的。 可以应用其他方法来确保 sum-to-one 约束是更好的模型，例如：

• a)正交投影:使用主成分分析(PCA)来识别在最小二乘意义上最能代表观测数据的仿射集，然后计算观测向量在该集合上的正交投影。这一投影如上图所示。
• b)透视投影:这是在[81]中提出的所谓的暗点固定变换(dark point fixed transform, dpft)。对于给定的观测向量y，这个投影，如上图所示，相当于根据$y/(y^{T}u)$缩放y，其中u对于数据集中的每一个y都满足$y^{T}u>0$。对于任意的$v\in {\mathbb{R}}^p$，包含投影向量的超平面定义为$v^{T}u=1$。

 值得注意的是，正交投影会改变光谱向量的方向，而透视投影不会。另一方面，透视投影引入了大尺度因子，对于接近与u正交的光谱向量，这些因子可能会变为负值。此外，具有不同角度的向量u产生非平行仿射集，从而产生不同的丰度系数，这意味着u的选择对于精确估计来说是一个关键问题。

  对于Rterrain数据集，效果如下图所示。Rterrain数据集是一个公开的高光谱数据立方体，由美国陆军工程兵团陆军地理空间中心分发，由高光谱图像数据收集实验（HYDICE）收集。它的尺寸是307像素乘500行，210个光谱带。左侧的图形绘制了未投影向量和正交投影向量之间的角度，作为未投影向量范数的函数。较小范数的向量通常对应于阴影区域，其角度较高，约为1–7。右边的图绘制了投影向量的范数，作为未投影向量范数的函数。相应的比例系数大约在1/3和10之间变化。

 减轻投影误差的一种可能方法是丢弃有问题的投影：在透视投影的情况下，这些有问题的投影可以认为是投影向量和未投影向量之间的角度大于给定小阈值的向量；在正交投影的情况下，这些有问题的投影可以认为是具有非常小或负比例因子$y^{T}u$的向量。

Reference

[1] 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.