高光谱解混

第9章 稀疏约束的广义双线性高光谱遥感影像解混
解混是高光谱图像处中一项非常重要的任务，指从高光谱图像的混合像素中分解得到端元及丰度的过程，一般认为，端元代表图像中存在的纯物质；丰度代表某个像素中每个端元所占的百分比。线性混合模型（linear mixing model, LMM)假设观察到的光谱是一些端元的线性组合，它是一个简化的光谱模型，只考虑一阶散射光子，而忽略多个光子间的相互作用。尽管基于LMM的混合端元分解方法能得到有物理意义的结果，但是光谱混合模型中的非线性分量已在很多工作中被指出。基于线形混合模型的光谱解混理论方法已得到了广泛的研究与应用，但由于受实际地物间复杂关系以及大气散射的影响，光谱混合可能是非线性的，这就使得应用传统的基于线性光谱混合模型的解混结果难以满足精度要求。
9.1 高光谱中的解混问题
通过搭载在航空航天设备上的成像光谱仪可以得到拥有地物电磁波谱幅值信息的幅值数据立方体，但是这个数据立方体并不能直接应用于解混研究，成像光谱仪所获得的高光谱遥感图像在某些波段受到大气中的水蒸气和成像时仪器的抖动等因素的影响，使得数据需要经过一些技术上的处理，如大气校正与降维等，才能应用到具体的研究实际当中去。
9.1.3 解混
高光谱图像解混的主要目的就是要获得混合像元中包含的端元，即端元提取及其在像元中所占的比例即丰度估计。端元提取包括端元识别和端元生成两种手段，端元识别是从光谱数据中选择有代表意义或是满足一定条件的像元集作为整幅图像的端元集，并没有改变谱数据的原有信息；而端元生成的手段是生成原光谱数据中可能不存在的谱特征作为端元，可能导致得到的光谱特征并没有实际的物理意义。在给定高光谱数据和端元集谱特征情况下，丰度估计是一个组合优化的逆问题，可通过有效的方法予以解决。
目前高光谱图像解混主要应用到的方法可以分为三个类别，分别是基于几何、统计和稀疏的方法。基于几何方法的光谱解混主要思想是利用了光谱数据的所有像元分布在高维单形体空间或是在一个正的凸锥区域；基于统计方法的光谱解混旨在利用参数估计的技术来决定端元集和丰度矩阵；基于稀疏回归方法的光谱解混主要是利用光谱库作为字典这一先验知识，作为一种半监督的方式将分离转化为一个稀疏回归的问题，类似于压缩感知的稀疏问题；基于稀疏编码的光谱解混主要是直接利用原光谱数据作为字典，然后对字典不断进行更新迭代除去了对光谱库的依赖，但随之而来的问题是关于字典和光谱数据之间的校正问题。

根据文献袁静等人的分析，我们将介绍NMF、原型分析、贝叶斯和稀疏解混方法的优缺点。基于 NMF的解混方法具有较明确的物理意义， 较好的数学操作性， 克服了几何方法的缺陷， 能够处理高混合像元的情况， 能够为某些特定的谱变异情况建模。但其存在的问题是目标函数非凸且易产生虚假端元， 谱变异模型泛化能力较弱。为了克服 NMF 的非凸问题， 近两年出现了采用原型策略的光谱解混算法， 称之为基于原型分析的光谱解混算法。然而， 在高混合像元存在的场景中， 如何调整松弛因子和尺度参数使得其能够找到需要的端元依然存在难度，且原型分析方法中没有考虑谱变异情况。而基于贝叶斯理论的解混方法能有效地将光谱、端元和丰度的变异和不确定性融入解混模型， 并通过增加合理先验信息改善解混效果， 其代价是需要大量的计算。进一步分析， 由于地面的复杂性， 通过概率分布对复杂场景下的端元和丰度进行建模亦存在难。上述三种方法中容易出现的问题是不可靠的端元提取大大地降低丰度估计的精度。为了避免不可靠的端元提取带来的负面影响， 稀疏解混越来越受到业界的关注。其存在的问题是所建立的端元谱库的可迁移性较差， 导致采用端元谱库逼近真实光谱图像时无法避免真实环境的误差影响。

9.2 高光谱解混模型
9.2.1 线性光谱混合模型
混合像元解混的前提是需要明确和了解构成混合像元的机理，通过分析场景中的物质对光的反射和散射作用，得到描述混合像元的基本模型，再利用各种解混算法进行端元提取和丰度估计。
线性混合模型（LMM）是从宏观上把握物质对光的反射作用，认为反射面是地物成规则混合的，并且是平整的，不考虑物质之间对光的相互作用，混合像元的形成主要是由于成像光谱仪的低分辨率造成的。图9.2形象地阐述了线性光谱混合过程，在地面的三种物质分别为m1、m2和m3，所占区域的比例分别为a1、a2和a3，太阳光照射在地物上时它们各自的反射光并没有射光的干扰，而是平行进入成像光谱仪的传感器中形成一个混合像素点，混合现象发生在传感器中，得到的混合像元的反射率为y=m1a1+m2a2+m3*a3，是各个光谱特征的线性叠加。

Hyperspectral unmixing

For hyperspectral unmixing, the $L_{2,1}$-norm regularizar $f(X)={\sum }_{i=1}^M\Vert X^i{\Vert }_2$, where, the abundance matrix $X\in R^{M\times N}$, $X^i$ represents the ith row of the abundance matrix. This $L_{2,1}$-norm regularizar promotes the row sparsity of $X$, this is to say, each material exists in the finite zone for a specific scene.
Spatial Correlation: Considering the fact that only a few endmembers play a dominant role in a homogeneous region, the corresponding abundance for every homogeneous region should be sparse.

基于张量分解的光谱解混

高光谱数据压缩、降维、特征提取和去噪任务主要关注数据最小重构误差。但是解混关注的是分解得到的因子是否和物理混合模型一致。尽管CPD和Tucker分解已经被用于高光谱解混中，他们和线性混合模型之间关系却没有像矩阵分解那么明确。CPD分解的秩是定义为最小秩。1张量个数，其数目很难确定，是一个NP问题。在解混中，端元个数在大多数情况下可以由统计、几何或领域知识得到，但是通过这些方法得到的秩不能被看作是张量的秩。这是因为给定张量$\mathcal{Y}\in {\mathbb{R}}^{I\times J\times K}$，在实际应用中，它的CPD秩通常设置$\max (I,J,K)$或更大数值。这个数值比端元个数大的多，因此CPD分解不能直接应用于高光谱解混。事实上，CPD分解把张量所有模同等对待，与高光谱数据特性不符合。尽管高光谱数据可以看做是一个三阶张量，它的各个模可以根据空谱结构区分，即两个模表示空间信息，一个模表示光谱信患。Tucker分解可以根据数据的物理意义区分张量不同模。虽然光谱维度的秩可以直接等价为端元个数，但是Tucker分解的正交性和端元的线性物理混合模型不一致。另外，Tucker分解的核张量和各模的因子矩阵强关联性质，使得很难把核张量转化为丰度矩阵和端元形式，即Tucker分解和线性光谱混合模型之间的关系同样不明确。因此，CPD分解和Tucker分解都不适合高光谱解混，我们有必要提出另外一个基于张量分解并且和线性光谱混合模型一致的解混算法。

综述

在高光谱遥感影像中，由于大部分像元是由小部分端元混合而成的，因此对于一个像元，它的丰度对应于丰度矩阵的列，而且是列稀疏的。对于稀疏解混算法，由于光谱库足够大，而场景中的端元个数远远小于给定光谱库中光谱个数，所以丰度系数必然从全局上表现出行稀疏的特性。陈允杰等人得出行稀疏性比 $l_1$ 稀疏更能有效刻画丰度系数的稀疏模式，具有更强的稀疏性，所以他们为了刻画丰度矩阵的行稀疏性而在CLSUnSAL-TV中使用了联合稀疏正则项$\Vert \mathbf{X}{\Vert }_{2,1}$，其中$\Vert \mathbf{X}{\Vert }_{2,1}={\sum }_{i=1}^m\Vert {\mathbf{x}}_i{\Vert }_2$。

高光谱遥感影像局部相邻像元通常包含相似的物质,其光谱特征具有相似性，利用高光谱影像像元间的这种相关性，可以有效提高解混的准确性。进一步分析像元间的空间相关性，可以发现这种相关性通常表现为像元间的平滑性，由此SUnSAL-TV算法和CLSUnSAL-TV加入了TV正则项。

在SUnSAL算法中，Bioucas-Dias和Figueiredo引入了正则项$\Vert \mathbf{X}{\Vert }_{1,1}$。SUnSAL算法的优化目标函数为
min ⁡ X 1 2 ∥ Y − A X ∥ F 2 + λ ∥ X ∥ 1 , 1 + ι { 1 } ( 1 T X ) + ι R + ( X ) \min_\mathbf{X}\frac{1}{2}\|\mathbf{Y-AX}\|_F^2+\lambda\|\mathbf{X}\|_{1,1}+\iota_{\{1\}}(\mathbf{1}^T\mathbf{X})+\iota_{R+}(\mathbf{X}) Xmin​21​∥Y−AX∥F2​+λ∥X∥1,1​+ι{1}​(1TX)+ιR+​(X)

对于$\mathbf{X}\in {\mathbb{R}}^{m,n}$，$\Vert \mathbf{X}{\Vert }_{1,1}={\sum }_{i=1}^n\Vert x_i{\Vert }_1$，其中${\mathbf{x}}_i$为$\mathbf{X}$的第$i$列。
定义${\mathbf{H}}_h:{\mathbb{R}}^{m\times n}\to {\mathbb{R}}^{m\times n}$是计算$\mathbf{X}$水平方向差分的一个线性算子，即${\mathbf{H}}_h\mathbf{X}=[{\mathbf{d}}_1,{\mathbf{d}}_2,\cdots ,{\mathbf{d}}_n]$，其中${\mathbf{d}}_i={\mathbf{x}}_i-{\mathbf{x}}_{i_h}$，$i_h$是的$i$邻接像元。定义${\mathbf{H}}_v:{\mathbb{R}}^{m\times n}\to {\mathbb{R}}^{m\times n}$是计算$\mathbf{X}$垂直方向差分的一个线性算子，即${\mathbf{H}}_v\mathbf{X}=[{\mathbf{v}}_1,{\mathbf{v}}_2,\cdots ,{\mathbf{v}}_n]$，其中${\mathbf{v}}_i={\mathbf{x}}_i-{\mathbf{x}}_{i_v}$，$i_v$是的$i$邻接像元。于是，$\mathbf{H}\mathbf{X}=\left[\begin{array}{c}{\mathbf{H}}_h\mathbf{X}\\ {\mathbf{H}}_v\mathbf{X}\end{array}\right]$。在SUnSAL-TV算法中，Iordache等人引入了正则项$\Vert \mathbf{H}\mathbf{X}{\Vert }_{1,1}$。SUnSAL-TV算法的优化目标函数为
$$\underset{\mathbf{X}}{\min }\frac{1}{2}\Vert \mathbf{Y}-\mathbf{A}\mathbf{X}{\Vert }_F^2+\lambda \Vert \mathbf{X}{\Vert }_{1,1}+{\lambda }_{\text{TV}}\Vert \mathbf{H}\mathbf{X}{\Vert }_{1,1}\Vert +{\iota }_{R+}(\mathbf{X})$$
其中，${\iota }_{R+}(\mathbf{X})={\sum }_{i=1}^n{\iota }_{R+}({\mathbf{x}}_i)$是指示函数，${\mathbf{x}}_i$是$\mathbf{X}$的第$i$列。如果${\mathbf{x}}_i$属于非负象限和$+\infty$，则${\mathbf{x}}_i$非零。

在CLSUnSAL算法中，Iordache等人使用了联合稀疏正则项$\Vert \mathbf{X}{\Vert }_{2,1}$，其中$\Vert \mathbf{X}{\Vert }_{2,1}={\sum }_{i=1}^m\Vert {\mathbf{x}}_i{\Vert }_2$。CLSUnSAL算法的优化目标函数为
$$\underset{\mathbf{X}}{\min }\Vert \mathbf{Y}-\mathbf{A}\mathbf{X}{\Vert }_F^2+\lambda \Vert \mathbf{X}{\Vert }_{2,1}+{\iota }_{R+}(\mathbf{X})$$

陈允杰等人在CLSUnSAL算法优化目标函数之上加入了SUnSAL-TV算法所使用的，得到了一种基于协同稀疏和全变差的高光谱线性解混方法，即CLSUnSAL-TV。CLSUnSAL-TV的优化目标函数为
KaTeX parse error: No such environment: matrix* at position 8: \begin{̲m̲a̲t̲r̲i̲x̲*̲}̲[l] &\min_\math…

MV-NTF算法没有利用空间结构信息，高光谱影像$\mathbf{Y}\in {\mathbb{R}}^{R\times I\times J}$的优化目标函数为：
$$\underset{\mathbf{A},\mathbf{B},\mathbf{C}}{\min }\Vert \mathcal{Y}-\sum _{r=1}^R({\mathbf{A}}_r{\mathbf{B}}_r^T)\cdot {\mathbf{c}}_r{\Vert }_F^2+\lambda \Vert 1_{I\times J}-\mathbf{A}{\mathbf{B}}^T{\Vert }_F^2$$
其中，$R$为端元个数。

对于TV-RSNMF算法的目标函数，首先添加了加权稀疏正则项$\Vert \mathbf{W}\odot \mathbf{X}{\Vert }_1$，其中$\odot$为元素乘算符。对于灰度图像$\mathbf{x}\in {\mathbb{R}}^{m\times n}$，Beck和Teboulle给出的各向异性TV范数为
$$\begin{gathered} \Vert x{\Vert }_{\text{TV}}=\sum _{i=1}^{m-1}\sum _{j=1}^{n-1}|x_{i,j}-x_{i+1,j}|+|x_{i,j}-x_{i,j+1}| \\ +\sum _{i=1}^{m-1}|x_{i,n}-x_{i+1,n}|+\sum _{j=1}^{n-1}|x_{m,j}-x_{m,j+1}| \end{gathered}$$
当处理丰度时，将丰度矩阵$\mathbf{X}\in {\mathbb{R}}^{L\times N}$转换为丰度张量$\mathcal{X}\in {\mathbb{R}}^{L\times m\times n}$，则$\Vert \mathbf{X}{\Vert }_{\text{HTV}}={\sum }_{j=1}^L\Vert {\mathcal{X}}_{j,:,:}{\Vert }_{\text{TV}}$。于是，TV-RSNMF算法的优化目标函数为
$$\underset{\mathbf{A},\mathbf{X}}{\min }\Vert \mathbf{Y}-\mathbf{A}\mathbf{X}{\Vert }_F^2+\lambda \Vert \mathbf{W}\odot \mathbf{X}{\Vert }_1+\tau \Vert \mathbf{X}{\Vert }_{\text{HTV}}$$

参考文献

陈允杰, 葛魏东 and 孙乐, 2018. 一种基于协同稀疏和全变差的高光谱线性解混方法. 自动化学报, 44(1), pp.116-128.

Iordache, M. D., Bioucas-Dias, J. M., & Plaza, A. (2012). Total variation spatial regularization for sparse hyperspectral unmixing. IEEE Transactions on Geoscience and Remote Sensing, 50(11), 4484-4502.

Bioucas-Dias, J. M., & Figueiredo, M. A. (2010, June). Alternating direction algorithms for constrained sparse regression: Application to hyperspectral unmixing. In 2010 2nd Workshop on Hyperspectral Image and Signal Processing: Evolution in Remote Sensing (pp. 1-4). IEEE.

He, W., Zhang, H., & Zhang, L. (2017). Total variation regularized reweighted sparse nonnegative matrix factorization for hyperspectral unmixing. IEEE Transactions on Geoscience and Remote Sensing, 55(7), 3909-3921

Beck, A., & Teboulle, M. (2009). Fast gradient-based algorithms for constrained total variation image denoising and deblurring problems. IEEE transactions on image processing, 18(11), 2419-2434.

Qian, Y., Xiong, F., Zeng, S., Zhou, J., & Tang, Y. Y. (2016). Matrix-vector nonnegative tensor factorization for blind unmixing of hyperspectral imagery. IEEE Transactions on Geoscience and Remote Sensing, 55(3), 1776-1792.

Iordache, M. D., Bioucas-Dias, J. M., & Plaza, A. (2013). Collaborative sparse regression for hyperspectral unmixing. IEEE Transactions on geoscience and remote sensing, 52(1), 341-354.