高光谱解混（三）——基于几何和统计的线性光谱解混方法

GEOMETRICAL

 基于几何的方法基本可以分为两大类：基于纯像元（PP）的方法和基于最小体积（MV）的方法。此外，还有一些其他方法也将被讨论。 在Section Ⅱ已经简单提到：通过将最小体积（MV）（最小面积）单纯形拟合到数据来完成端元估计（给定点集寻找构成凸包的点），可以作为几种基于几何的解混算法的理论基础 。

基于纯像元的几何解混方法

 基于纯像元的算法仍属于最小体积MV方法，但假设数据中每个端元至少有一个相应的纯像元。这意味着数据单纯形的每个顶点上至少有一个光谱向量（可以理解为：单纯形内部的点是混合像元，单纯形的顶点是纯净像元—即，端元）。尽管从计算的角度来看，这种假设能够设计出非常高效的算法，但这是一个强大的必要条件，在许多数据集中可能并不适用。在任何情况下，这些算法都能找到数据中最纯的像元集。
 基于纯像元假设的线性几何解混算法的计算复杂度低和理论意义明确，因此在高光谱解混算法中得到广泛应用。

像元纯度指数(PPI)

 PPI算法使用最大噪声分量 MNF作为预处理步骤来降维，提高信噪比。

 理论依据：如果将高光谱影像所有像元组成一个特征空间，那么以影像的端元作为顶点所组成的单形体可以将高光谱影像中所有的像元都包含在内。将高光谱数据的光谱向量投影到每个测试向量上的时候，端元只会投影到向量的两端，而混合像元则会投影到向量的中部。如下图所示，随机产生了几个穿过单纯形的测试向量$X_1,X_2,X_3$，端元$a_i$在测试向量上的投影始终位于最外端。混合像元的投影位于测试向量上两个端元投影点的中间。

 为了减少模型和噪声对算法性能的影响，在开始的时候先给每个像元定义一个纯净像元指数（初始值为零)，用以记录每个像元投影到测试向量端点的次数。每个像元的纯像元指数越大，意味着该像元成为端元的可能性越大。

 实施步骤：考虑原始高光谱图像中噪声的存在，为了降低噪声的影响，先将高光谱图像使用MNF或者PCA进行降维（一般降至q-1维，q为端元数量）和去噪，然后再使用PPI 算法。

• 利用NMF或者PCA变换将原始高光谱图像进行降维和去噪，获取低维数据。
• 在进行变换后得到的低维数据空间中随机生成大量的测试向量。
• 将该数据空间中所有的数据投影到测试向量中，统计每个点的纯净像元指数。
• 从所有点中选择纯净像元指数最大的q个点作为端元。

 缺点：

• 需要预先设置随机单位向量的个数与选择端元时所需要的阈值q(即需要预先知道端元个数)；
• 向量的随机性使得每次求解的结果不一致；
• 每次使用的测试向量的数量比较大，因此PPI算法的速度不是很快

 2006年chang等人基于奈曼一皮尔逊引理检测理论的阈值方法(HFC)来确定高光谱数据中光谱端元的数量——虚拟维数(VD)，然后将其作为随机变量的个数，避免了阈值的设定。

N-FINDR

 理论依据：在光谱维度中，由最纯像元形成的单纯形定义的体积大于由任何其他像元组合定义的任何其他体积。

 基本步骤：如下图所示，从高光谱图像中随机选择一定数量的像元（假设为p )作为光谱端元，计算由这些像元所构成的单形体的体积。固定这些像元中的p-1个像元，使用其余像元依次替代第p个像像元，如果单形体体积增加，则替代原来的像素点，否则无效则淘汰。依次循环，直到单形体的体积没有变化或者变化很小为止。同理，考虑原始高光谱图像中噪声的存在，为了降低噪声的影响，需要先对高光谱图像进行降维处理。

 缺点：该算法便于理解，低噪声情况下识别效果好。但也存在如下缺点：

• 受搜索策略的限制难以收敛到全局最优解，且像元个数较多时时间复杂度高;
• 需要预先知道端元个数;
• 最终获取的端元提取结果与初始选定的作为端元的像元有关。
  
  

VCA顶点成分分析

 理论依据：1）端元为单形体的顶点；2）单形体经仿射变换后仍为单形体。

 基本步骤：顶点成分分析(vertex component analysis, VCA)算法迭代地将高光谱像元向量投影到《已经确定的端元张成的子空间》正交的方向上，然后选定投影最大的像元作为新的端元加入已知端元集。该算法迭代直到成功的选择所有端元。

 缺点：每次寻找新端元时都要遍历所有的像元点，因此VCA算法具有较高的时间复杂度。初始端元的选取对最终端元提取的结果也有一定的影响。

迭代误差分析（IEA）算法

迭代误差分析（IEA）算法使用一系列线性约束解混，每次都选择使得纯净图像中剩余误差最小化的像元作为端元 。

SMACC

 序列最大角凸锥(SMACC)算法基于凸锥来表示谱向量。该算法从一个端元开始，并在维数上递增。根据端元与现有圆锥的夹角来识别新的端元。选择与现有圆锥最大夹角的数据向量作为下一个端元，扩大端元集。在一定的公差范围内，当所有的像元向量都在凸锥内时，算法终止。

基于N-FINDR

• SGA单形体增长：SGA算法是由Chang等人在总结N-FINDR算法局限性的基础上提出的。理论依据：高光谱数据位于端元所构成的单形体中，当数据中存在每种物质的纯像元时，这些纯像元相较于其他混合像元所构成的单形体体积最大。基本步骤：SGA算法相对于N-FINDR算法具有以下几点进步。首先，SGA算法采用虚拟维度、信号子空间估计等方法获取高光谱图像所包含的端元数目，从而不需要设置算法的阈值终止条件；其次，该算法采用逐步增加单形体维度的方式逐个识别目标端元，极大的降低了算法执行的时间复杂度；最后，该算法专门设计了起始端元的初始化方法，极大保障了算法在重复运行时结果的一致性。
• AVMAX:受N-FINDR启发，交替体积最大化(AVMAX)以循环的方式最大化由端元定义的单形的体积，每次针对一个端元。AVMAX与N-FINDR的SC-N-FINDR变体非常相似。

基于VCA

• 连续体积最大化（SVMAX）：连续体积最大化(SVMAX)类似于VCA。主要的区别在于像元向量向已经确定的端元张成的子空间正交的方向上进行投影的方式。VCA考虑的是这些子空间中的一个随机方向，而SVMAX考虑的是完整的子空间。

Other

• 协作凸框架(collaborative convex framework)将像元向量矩阵Y分解为一个非负混合矩阵M和一个稀疏且非负丰度矩阵S。混合矩阵M的列被约束为数据Y的列。
• 晶格联想记忆（LAM） 将光谱集建模为部分有序实值向量的晶格元素。 格操作用于非线性构造 LAMS。 通过从光谱像素构造所谓的最小和最大 LAM 来找到端元。 这些 LAM 包含光谱像素的最大和最小坐标（经过适当的附加缩放）并且是候选端元。 端元是使用仿射独立性和相似性度量的概念从 LAMS 中选择的，例如光谱角、相关性、互信息或切比雪夫距离。

基于最小体积的几何解混算法

基本概念

 MV方法寻求一个混合矩阵M，该混合矩阵能够最小化由其列向量构成的单形体$conv(\mathrm{M})$的体积。同时该单形体遵循如下约束：$conv(\mathrm{M})$中包含观测到的光谱向量。这个约束可以是软约束，也可以是硬约束。此处不再强制遵循纯像元假设，因此解混问题也变成了一个非凸优化问题。下图对包含所有观测数据的最小单形体进行抽象说明。估计的混合矩阵$\hat{M}=[\widehat{m_1},\widehat{m_2},\widehat{m_3}]$与真正的混合矩阵略有不同，因为每个面（边）上没有足够的数据点($p-1$个面/边)来定义真正的单纯形。

 假设数据集投影到子空间$S$，$S$的维度为$p$，$S$的列向量$m_i\in {\mathbb{R}}^p，i=1,...,p$是仿射独立的（或者说$m_i-m_1,i=2,...,p$是线性独立的）。因此，单纯形$conv(\mathrm{M})$的维数是p，在${\mathbb{R}}^p$中$conv(\mathrm{M})$的体积是0。要得到非零体积，通常考虑包含原点的扩展单纯形${\mathrm{M}}_0\equiv [0,\mathrm{M}]$。${\mathrm{M}}_0$的体积/凸包被定义为（det代表求行列式，行列式的物理意义就是${\mathbb{R}}^n$维空间下n维空间几何的体积）

 对上式，如果将数据集平移到原点，并在维度等于$p-1$的子空间内工作时，单纯形的体积可以表示为

MVT：最小体积变换

 理论依据：Craig于1994年提出了MVT算法。不同于单形体体积增长类方法，MVT更适合于像元高度混合的高光谱图像解混方法。因为在像元高度混合的高光谱图像中，纯像元假设通常不再成立，基于单形体体积增长类方法所检索到的端元仍然可能是混合像元。相较之下，MVT致力于寻找包围所有高光谱数据点的最小单形体，该单形体的顶点更接近于真实端元。

个人理解：
 N-FINDR类的单形体体积增长方法将认为所有的观测谱向量中必有纯净像元且纯净像元是单形体的顶点。所以该类方法以像元向量为出发点，使选中的几个像元向量构成的单形体包含所有的观测谱向量，固定p-1个像元，替换第p个像元以寻求可以组成的最大体积单形体。
 实际上，高光谱图像可能并不存在纯净像元，也就没有端元是单形体的顶点这一概念，所以MVT类算法则是以单形体为出发点，直接使拟合单形体去包围所有的观测谱向量，固定拟合的一个面/一条边，迭代求解其他的面/边，从而寻找到最小体积的单形体。

 具体步骤：该算法每次迭代地改变单纯形的一个面/一条边，保持其他面/边不变，从而保证体积最小化且所有的谱向量都包含在形成的单形体中（ANC约束和ASC约束：$\alpha ={\hat{\mathbf{M}}}^{-1}{\mathbf{y}}_i⪰0\text{ and }{1}_p^T{\hat{\mathbf{M}}}^{-1}{\mathbf{y}}_i=1$）。MVT算法中单形体的体积可以表示为

算法核心部分可以用公式描述为

MVSA/SISAL、MVES

 在MVT思想的基础上，Li等人提出的最小体积单纯形分析(MVSA)和J. M. Bioucas-Dias提出的基于变量分裂和增强拉格朗日(SISAL)的单纯形识别算法实现了基于MV的鲁棒几何解混。MVSA 以 VCA算法的结果进行初始化,可以摒除初始单形体内的所有搜索区域，然后采用序列二次规划方法搜索最终端元。鲁棒性是通过允许违反非负性约束引入的【将非负性约束ANC变为软约束】。为了理解这种修改的相关性，下图描绘了有噪声的光谱观测数据。虚线单纯形表示包含所有数据所需的最小体积单纯形；通过允许违反非负性约束，MVSA和SISAL算法产生了一个非常接近真实单纯形的单纯形。由于噪声或任何其他扰动源的存在，谱向量可能位于真正的数据单纯形之外。使用MV算法会产生虚线估计，虚线估计形成的单纯形与原始真正单纯形有一定的区别。

 为了更准确地估计端元，MVSA/SISAL允许对非负性约束的违反。这种弹性违反约束使用惩罚函数$hinge(x)=0,ifx\ge 0and-xifx<0$来表示。MVSA/SISAL将数据投影到信号子空间，矩阵M是方形的，理论上是可逆的（病态可能会使数值计算逆矩阵变得困难）。此外，逆问题求解丰度向量$x$可以表示为：

 MVSA/SISAL旨在解决以下优化问题：

其中，$-\lambda 1_p^T\mathrm{hinge}(\mathbf{Q}\mathbf{Y})1_n$权衡了非负性约束ANC的违反程度。当正则化参数$\lambda$趋近于$\infty$时，软约束趋向于硬约束。

 最小体积封闭单纯形（MVES）旨在解决优化问题【上图中（10）】在$\lambda$=$\infty$的情况，即硬约束情况。MVES使用线性规划（LPs）实现循环最小化。虽然优化问题（10）是非凸的，但在（10）中证明了纯像元的存在是MVES识别真正端元的充分条件。

MVC-NMF

 MVC-NMF在非负矩阵分解算法的框架中采用最小化单形体体积作为正则约束项，解决了直接应用于原始数据集的无降维优化问题：

其中，$S\equiv [{\alpha }_1,...,{\alpha }_n]\in {\mathbb{R}}^{p\times n}$是包含丰度向量的矩阵，$\Vert \mathbf{A}{\Vert }_F^2\equiv \mathrm{tr}\left({\mathbf{A}}^T\mathbf{A}\right)$是矩阵$\mathbf{A}$的F范数，$\lambda$是正则化参数。(11)最小化优化一个两项目标函数，前项衡量近似误差，后项衡量$\mathbf{M}$的列向量构成的单纯形体积的平方。(个人理解：前项代表最小化噪声误差，后项代表最小化单形体体积，正则化参数平衡最小化单形体体积和最小化重建误差)。MVC-NMF针对$\mathbf{S}$（二次规划问题）和$\mathbf{M}$（非凸规划问题）实现了一系列交替最小化。 MVC-NMF和MVSA/SISAL/MVES算法之间的主要区别在于，后者允许违反ANC，从而使SU逆问题具有鲁棒性，而前者不允许。

ICE

 迭代约束端元(ICE)算法旨在解决一个与MVC-NMF类似的优化问题，其中单纯形的体积被一个更易于处理的近似值所取代——所有单纯形顶点之间距离的平方和。ICE通过优化所有单形体顶点之间的距离平方和来获取端元。ICE中的正则化项是二次的，在任何环境维度和退化单形中都能被很好的定义，这也是相较于$\Vert det(\mathbf{M})\Vert$正则化项的一个优势（$\Vert det(\mathbf{M})\Vert$正则化项是非凸的，当HU问题条件恶劣或端元的数量不准确时，容易出现并发问题 ）。由于使用了二次正则化项，ICE相对于MVC-NMF的一个优点是，在前者中，最小化是一个二次规划问题，它是可以解析求解的最小二乘问题。而在MVC-NMF中，关于M的优化是一个非凸问题。

 ICE算法的变体是SPICE。SPICE结合了稀疏促进的先验知识，旨在寻找正确的端元数量。SPICE在二次目标函数添加了一个与一个端元相关的所有比例的线性项。线性项对应于指数先验。初始化过程中使用了大量端元，迭代过程中优先倾向于将与特定端元相关的所有比例值变为零。如果与一个端元对应的所有比例都为零，则可以丢弃该端元。由于最小化仍然涉及一个二次规划，因此增加稀疏性促进先验并不会增加模型的复杂性。
 ICE和SPICE算法中使用的二次体积正则化项在允许数据点位于单纯形$conv\mathbf{M}$之外的条件上也提供了一定的鲁棒性。ICE最小化目标函数定义为：

其中，$\mu \in [0,1]$是一个正则化参数，用于控制误差和更小的单形体体积之间的权衡。如果$\mu =1$，最优解是将所有端元收缩到一个点，这样所有数据都将位于单纯形之外。 如果$\mu =0$，最优解是不考虑单纯形的大小的无误解。${\sum }_{\mathbf{M}}$表示端元的样本协方差矩阵，$trace({\sum }_{\mathbf{M}})$表示所有单纯形顶点之间距离的平方和。

$L_{1/2}$-NMF

有一部分因为赶时间没有仔细看和翻译，先占个坑~

STATISTICAL

 当光谱混合程度较高时，由于单形体的面上没有足够的谱向量且不一定存在纯像元，基于最小体积的几何的方法效果较差（无法找到准确的拟合单形体）。在这些情况下，统计方法是一种强大的替代方法，但通常需要付出代价:与基于几何的方法相比，它的计算复杂度更高。统计方法还提供了一个自然的框架来表示端元的可变性。在统计框架下，谱分解是一个统计推理问题。
 在大多数情况下，端元的数量和它们的反射率是未知的，因此高光谱解混问题也归属于一种盲源分解问题。独立分量分析(Independent Component Analysis, ICA)是一种有名的盲源分离技术，目前已提出了基于独立分量分析的高光谱数据盲解混的几种算法[155]- [157]。然而，ICA基于“源相互独立的假设”，但高光谱数据的丰度向量满足丰度和为一约束，即源之间存在统计依赖性。这种依赖性削弱了ICA对高光谱数据的适用性。实际上，ICA通过将光谱向量与解混矩阵相乘寻找端元特征，该解混矩阵最小化通道间的互信息。如果源是独立的，ICA能够正确的解混。因为通道间互信息的最小值是与独立的源相对应的。对于高光谱数据所依赖的丰度向量而言，这已不再成立。然而，有些端元可能是近似纯净的。
 贝叶斯方法能够为统计可变性建模，并施加先验，将所求解限制在物理意义范围内，并将所求解正则化。后一性质通常被认为是解决不适定问题的必要条件。推理机采用贝叶斯框架，是待估计随机量的后验密度。假设未知混合矩阵M和丰度系数矩阵S是先验独立的，Bayes范式允许M和S的联合后验使用如下方式计算：

 $p_A$和$p_{A|B}$分别代表A和给定B的情况下A的概率密度函数pdf。上式中，$p_{Y|M,S}$是依赖于观测模型和先验分布$p_M(\mathbf{M})\&p_s(\mathbf{S})$的似然函数。先验分布总结了关于这些未知参数的先验知识。
 一个最常用的的贝叶斯估计是联合最大后验估计：

 在线性混合模型下，假设随机噪声向量$\mathbf{w}$服从高斯分布且协方差为${\sigma }^2\mathbf{I}$。由此则有：

 由此，ICE/SPICE和MVC-NMF算法，可以被归类为几何算法，也可以被归类为统计算法，从而产生(15)中的联合MAP估计。在所有这些算法中，端元估计都是通过最小化一个两项目标函数来获得的。其中$-lo{g}_{p_{Y|M,S}}(\mathbf{Y}|\mathbf{M},\mathbf{S})$用于数据拟合，${\log }_{p_M}(\mathbf{M})-{\log }_{p_S}(\mathbf{S})$组成惩罚项。相反，从贝叶斯的角度来看，分别指定端元混合矩阵M和丰度矩阵A的先验分布是确保观测模型固有物理约束的一种便捷的方法。
 参考文献[160]引入了一种贝叶斯方法，其中假设线性解混模型是零均值高斯白噪声，协方差为${\sigma }^2\mathbf{I}$。且丰度系数均匀分布在单纯形上，混合矩阵M的先验是一个自回归模型。[160]以迭代的方式实现负对数后验分布的最大化。丰度系数的最大化被表述为具有线性约束的n个加权最小二乘问题分别求解。与混合矩阵M相关的优化使用梯度下降法进行求解。

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.