基于原型分析的高光谱图像处理

原型分析简介

在高光谱图像处理中,原型分析(Archetypal Analysis, AA)是一种无监督的端元提取和丰度估计方法。与传统的端元提取方法不同,原型分析并不假设端元位于单纯形的顶点,而是寻找原始数据点的最小凸包。

给定高光谱数据矩阵 \(\mathbf{X}=\left\{\mathbf{x}_{1}, \mathbf{x}_{2}, \ldots, \mathbf{x}_{n}, \ldots \mathbf{x}_{N}\right\} \in \mathbf{R}^{M \times N}\),其中 \(M\) 为谱带数, \(N\) 为像素数。原型分析的目标是从 \(\mathbf{X}\) 中估计出 \(D(D\ll N)\) 个原型 \(\mathbf{Z}=\left\{\mathbf{z}_{1}, \mathbf{z}_{2}, \ldots, \mathbf{z}_{d}, \ldots \mathbf{z}_{D}\right\} \in \mathbf{R}^{M \times D}\),使得每个原始数据 \(\mathbf{x}_{n}\) 可以表示为原型的线性组合:

\[\mathbf{x}_{n}=\mathbf{Z} s_{n}\]

相应地,每个原型 \(\mathbf{z}_{d}\) 也可以表示为原始数据的线性组合:

\[\mathbf{z}_{d}=\mathbf{X} \mathbf{c}_{d}\]

其中 \(\mathbf{C} \in \mathbf{R}^{N \times D}\) 为系数矩阵, \(\mathbf{S} \in \mathbf{R}^{D \times N}\) 为丰度矩阵。

原型分析的目标函数可以表示为:

\[ \begin{array}{ll} \min _{\mathbf{C}, \mathbf{S}} & f(\mathbf{X} \mid \mathbf{X C S})=\|\mathbf{X}-\mathbf{X C S}\|_{F}^{2} \\ \text { s.t. } & \mathbf{C} \geq 0, \mathbf{S} \geq 0 \\ & \left\|\mathbf{c}_{d}\right\|_{1}=1,\left\|\mathbf{s}_{n}\right\|_{1}=1, \end{array} \]

其中 \(\|\cdot\|_F\) 表示矩阵的范数。约束条件保证了系数矩阵 \(\mathbf{C}\) 和丰度矩阵 \(\mathbf{S}\) 的非负性,以及各行(列)元素和为 1。

为了提高原型分析的鲁棒性,可以引入松弛因子 \(\boldsymbol{\alpha}\),使得原型不必严格位于原始数据点的凸包内。此时目标函数变为:

\[ \begin{array}{rl} \min _{\boldsymbol{\alpha}, \mathbf{C}, \mathbf{S}} & f(\mathbf{X} \mid \mathbf{X} \mathbf{C} \operatorname{diag}(\boldsymbol{\alpha}) \mathbf{S}) \\ \text { s.t. } & 1-\delta \leq \alpha_{d} \leq 1+\delta \\ & \mathbf{C} \geq 0, \mathbf{S} \geq 0 \\ & \left|\mathbf{c}_{d}\right|_{1}=1,\left|\mathbf{s}_{n}\right|_{1}=1 \end{array} \]

其中 \(\delta\) 是控制原型偏离程度的参数。

基于原型分析的高光谱图像融合

高光谱图像融合是将低空间分辨率的高光谱图像与高空间分辨率的多光谱图像融合,获得同时具有高光谱和高空间分辨率的图像。假设理想的高光谱图像为 \(\mathbf{Z}\),低空间分辨率的高光谱图像为 \(\mathbf{Y}\),高空间分辨率的多光谱图像为 \(\mathbf{X}\),它们可以表示为:

\[ \begin{aligned} \mathbf{Z} &= \mathbf{Z}\mathbf{C}\mathbf{S} \\ \mathbf{Y} &= \mathbf{Z}\mathbf{P} = \mathbf{Z}\mathbf{C}\mathbf{S}\mathbf{P} \\ \mathbf{X} &= \mathbf{R}\mathbf{Z} = \mathbf{R}\mathbf{Z}\mathbf{C}\mathbf{S} \end{aligned} \]

其中 \(\mathbf{P}\) 为空间下采样矩阵, \(\mathbf{R}\) 为光谱响应矩阵。根据最小二乘准则,融合目标函数可以表示为:

\[ \begin{aligned} \min_{\mathbf{Z},\mathbf{C},\mathbf{S}}&\frac{1}{2} \left\|\mathbf{Y} - \mathbf{Z}\mathbf{C}\mathbf{S}\mathbf{P}\right\|_F^2 + \frac{1}{2} \left\|\mathbf{X} -\mathbf{R}\mathbf{Z}\mathbf{C}\mathbf{S} \right\|_F^2 \end{aligned} \]

引入松弛因子 \(\boldsymbol{\alpha}\) 后,目标函数变为:

\[ \begin{aligned} \min_{\mathbf{Z},\mathbf{C},\alpha, \mathbf{S}}&\frac{1}{2} \left\|\mathbf{Y} - \mathbf{Z}\mathbf{C}\operatorname{diag}(\alpha)\mathbf{S}\mathbf{P}\right\|_F^2 + \frac{1}{2} \left\|\mathbf{X} -\mathbf{R}\mathbf{Z}\mathbf{C}\operatorname{diag}(\alpha)\mathbf{S} \right\|_F^2\\ \text { s.t. } & 1-\delta \leq \alpha_{d} \leq 1+\delta \\ \end{aligned} \]

由于 \(\mathbf{Z}\) 是未知量,因此该模型无法直接求解。一种常见的做法是先使用端元提取算法(如 PPI、FPPI 或 VCA)从 \(\mathbf{Y}\) 估计出过完备的端元矩阵 \(\mathbf{Q}\),然后利用 \(\mathbf{Q}\) 代替 \(\mathbf{Z}\),目标函数简化为:

\[ \begin{aligned} \min_{\mathbf{C},\alpha, \mathbf{S}}&\frac{1}{2} \left\|\mathbf{Y} - \mathbf{Q}\mathbf{C}\operatorname{diag}(\alpha)\mathbf{S}\mathbf{P}\right\|_F^2 + \frac{1}{2} \left\|\mathbf{X} -\mathbf{R}\mathbf{Q}\mathbf{C}\operatorname{diag}(\alpha)\mathbf{S} \right\|_F^2\\ \text { s.t. } & 1-\delta \leq \alpha_{d} \leq 1+\delta \\ \end{aligned} \]

为了获得更稀疏的解,可以引入 \(\ell_1\) 范数正则化项:

\[ \begin{aligned} \min_{\mathbf{C},\alpha, \mathbf{S}}&~\frac{1}{2} \left\|\mathbf{Y} - \mathbf{Q}\mathbf{C}\operatorname{diag}(\alpha)\mathbf{S}\mathbf{P}\right\|_F^2 + \frac{1}{2} \left\|\mathbf{X} -\mathbf{R}\mathbf{Q}\mathbf{C}\operatorname{diag}(\alpha)\mathbf{S} \right\|_F^2+\left\|\mathbf{S} \right\|_1\\ \text { s.t. } &~1-\delta \leq \alpha_{d} \leq 1+\delta \\ &\mathbf{C}\succeq \mathbf{0}\\ & \mathbf{S} \succeq \mathbf{0}\\ \end{aligned} \]

该目标函数可以通过交替优化的方式求解,具体更新步骤如下:

\[ \begin{aligned} \mathbf{C}^{k + 1}\in &\mathop{argmin}\limits_{\mathbf{C}\succeq \mathbf{0}}\frac{1}{2} \left\|\mathbf{Y} - \mathbf{Q}\mathbf{C}\operatorname{diag}(\alpha)\mathbf{S}\mathbf{P}\right\|_F^2 + \frac{1}{2} \left\|\mathbf{X} -\mathbf{R}\mathbf{Q}\mathbf{C}\operatorname{diag}(\alpha)\mathbf{S} \right\|_F^2\\ \alpha^{k + 1}\in &\mathop{argmin}\limits_{1-\delta \leq \alpha_{d} \leq 1+\delta}\frac{1}{2} \left\|\mathbf{Y} - \mathbf{Q}\mathbf{C}^{k + 1}\operatorname{diag}(\alpha)\mathbf{S}\mathbf{P}\right\|_F^2 + \frac{1}{2} \left\|\mathbf{X} -\mathbf{R}\mathbf{Q}\mathbf{C}^{k + 1}\operatorname{diag}(\alpha)\mathbf{S} \right\|_F^2\\ \mathbf{S}^{k + 1}\in~ &\mathop{argmin}\limits_{\mathbf{S}\succeq \mathbf{0}}\frac{1}{2} \left\|\mathbf{Y} - \mathbf{Q}\mathbf{C}^{k + 1}\operatorname{diag}(\alpha^{k + 1})\mathbf{S}\mathbf{P}\right\|_F^2 + \frac{1}{2} \left\|\mathbf{X} -\mathbf{R}\mathbf{Q}\mathbf{C}^{k + 1}\operatorname{diag}(\alpha^{k + 1})\mathbf{S} \right\|_F^2+\left\|\mathbf{S} \right\|_1\\ \end{aligned} \]

其中 \(\mathbf{C}\) 和 \(\boldsymbol{\alpha}\) 的更新采用投影梯度法, \(\mathbf{S}\) 的更新利用 ADMM 算法求解。

基于原型分析的高光谱图像解混

在高光谱图像解混中,原型分析也可以用于估计端元矩阵 \(\mathbf{A}\) 和丰度矩阵 \(\mathbf{S}\),使得:

\[ \begin{aligned} &\min \frac{1}{2}\|\mathbf{R}-\mathbf{R} \mathbf{C} \mathbf{S}\|_{F}^{2}\\ &\text { s.t. } \mathbf{C} \geq \mathbf{0}, \mathbf{1}_{n}^{T} \mathbf{C}=\mathbf{1}_{p}^{T}, \mathbf{S} \geq \mathbf{0}, \mathbf{1}_{p}^{T} \mathbf{S}=\mathbf{1}_{n}^{T} \end{aligned} \]

其中 \(\mathbf{R}\) 为高光谱数据矩阵, \(n\) 为像素数, \(p\) 为端元数。

考虑到光谱易变性,可以先利用 PPI 算法从 \(\mathbf{R}\) 中提取出过完备的数据点集 \(\mathbf{Y}\),然后利用约束原型分析估计端元和丰度:

\[ \begin{aligned} \arg \min &\frac{1}{2}|| \mathbf{Y} \mathbf{C} \operatorname{diag}(\boldsymbol{\alpha}) \mathbf{S}-\mathbf{R}\left\|_{F}^{2}+\lambda|| \mathbf{S}\right\|_{1}\\ &\text { s.t. } \mathbf{C} \geq \mathbf{0}, 1-\delta \leq \alpha_{i} \leq 1+\delta\\ &\mathbf{S} \geq \mathbf{0},\left\|c_{i}\right\|_{1}=1,\left\|s_{j}\right\|_{1} \leq 1 \end{aligned} \]

其中 \(\lambda\) 为正则化参数,用于控制丰度矩阵的稀疏性。该优化问题可以通过交替优化的方式求解,具体更新步骤如下:

\[ \begin{aligned} \boldsymbol{C}^{k+1} &=\underset{\boldsymbol{C} \geq 0,\left\|c_{i}\right\|_{1}=1}{\operatorname{argmin}} F\left(\boldsymbol{C}^{k}, \boldsymbol{\alpha}, \boldsymbol{S}\right) \\ \boldsymbol{\alpha}^{k+1} &=\underset{1-\delta \leq \alpha_{i} \leq 1+\delta}{\operatorname{argmin}} F\left(\boldsymbol{C}, \boldsymbol{\alpha}^{k}, \boldsymbol{S}\right) \\ \boldsymbol{S}^{k+1} &=\underset{\boldsymbol{S} \geq \mathbf{0},\left\|s_{j}\right\|_{1} \leq 1}{\operatorname{argmin}} F\left(\boldsymbol{C}, \boldsymbol{S}^{k}\right) \end{aligned} \]

其中 \(F\) 为目标函数:

\[ F=\frac{1}{2}\left\|\mathbf{Y} \widehat{\mathbf{C}} \mathbf{S}-\mathbf{R}\right\|_{F}^{2}+\lambda\|\mathbf{S}\|_{1} \]

\(\widehat{\mathbf{C}}\) 表示 \(\tilde{\mathbf{C}} \operatorname{diag}(\boldsymbol{\alpha})\)。

\(\mathbf{C}\) 和 \(\boldsymbol{\alpha}\) 的更新采用投影梯度法:

\[ \begin{aligned} \mathbf{C}^{k+1}&=\max \left(0, \mathbf{C}^{k}-\mu_{C}^{k} \nabla_{C} F\left(\mathbf{C}^{k}, \mathbf{S}^{k}\right)\right)\\ \boldsymbol{\alpha}^{k+1}&=P_{\alpha}\left[\boldsymbol{\alpha}^{k}-\mu_{\alpha}^{k} \nabla_{\alpha} F\left(\boldsymbol{\alpha}^{k}\right)\right] \end{aligned} \]

其中 \(\mu_C^k\) 和 \(\mu_\alpha^k\) 为步长,通过线搜索获得。\(P_\alpha\) 为投影算子,将 \(\boldsymbol{\alpha}\) 约束在 \([1-\delta, 1+\delta]\) 范围内。

\(\mathbf{S}\) 的更新采用 ADMM 算法求解:

\[ \begin{aligned} \mathbf{z}^{k+1}&=\operatorname{argmin} \lambda\left\|\mathbf{z}^{k}\right\|_{1}+\frac{\mu}{2}\left\|\mathbf{S}^{k}-\mathbf{z}^{k}-\mathbf{d}^{k}\right\|_{F}^{2}\\ \mathbf{S}^{k+1}&=\underset{\mathbf{S} \geq \mathbf{0},\left\|s_{j}\right\|_{1} \leq 1}{\operatorname{argmin}} \frac{1}{2}\left\|\mathbf{A} \mathbf{S}-\mathbf{R}\right\|_{F}^{2}+\frac{\mu}{2}\left\|\mathbf{S}-\mathbf{z}^{k+1} \mathbf{d}\right\|_{F}^{2}\\ \mathbf{d}^{k+1}&=\mathbf{d}^{k}-\left(\mathbf{S}^{k+1}-\mathbf{z}^{k+1}\right) \end{aligned} \]

其中 \(\mathbf{A}=\mathbf{Y} \widehat{\mathbf{C}}\), \(\mathbf{z}\) 和 \(\mathbf{d}\) 为辅助变量。

\(\mathbf{z}\) 的更新可以通过软阈值算子高效完成:

\[ \mathbf{z}^{k+1}=\operatorname{soft}\left(\mathbf{z}^{k+1}-\mathbf{d}, \frac{\lambda}{\mu}\right) \]

\(\mathbf{S}\) 的更新借鉴了快速梯度映射(FGM)算法,具体步骤如下:

Input: A matrix: $\boldsymbol{A} \in R^{L \times p}$, a matrix $\boldsymbol{R} \in R^{L \times n}$ and the initial $\boldsymbol{S}=\boldsymbol{S}_{0}$.

Output: An approximate solution $\boldsymbol{S} \approx \underset{\boldsymbol{S} \in \Delta^{p}}{\operatorname{argmin}} f(\boldsymbol{A} \boldsymbol{S}-\boldsymbol{R})$.

1. $\alpha_{0} \in(0,1) ; \boldsymbol{Q}=\boldsymbol{S} ; L^{*}=\left\|L_{f}\right\|_{F}^{2}$  

2. for $k=1:$ maxtier do

3.    $\boldsymbol{S}^{*}=\boldsymbol{S}$.

4.    $\boldsymbol{S}=P_{\Delta}\left(\boldsymbol{Q}-\frac{1}{L^{*}} \nabla f(\boldsymbol{A} \boldsymbol{Q}-\boldsymbol{R})\right)$ 

5.    $\mathbf{Q}=\boldsymbol{S}+\beta_{k}\left(\boldsymbol{S}-\boldsymbol{S}^{*}\right)$, 其中 $\beta_{k}=\frac{\alpha_{k}\left(1-\alpha_{k}\right)}{\alpha_{k}^{2}+\alpha_{k+1}}$, $\alpha_{k+1} \geq 0$ 使得 $\alpha_{k+1}^{2}=\left(1-\alpha_{k+1}\right) \alpha_{k}^{2}$

6. end for

其中 \(\Delta^p\) 为单位简单体, \(P_\Delta\) 为投影算子,将 \(\mathbf{S}\) 约束在 \(\Delta^p\) 内。 \(L_f\) 为目标函数 \(f\) 的利普希茨常数。

总的算法流程如下:

Input: The hyperspectral data $\boldsymbol{R} \in R^{L \times n}$, the number of endmembers $p$, and the parameter $\delta$ that regulates the size of the archetypal in (5).

Output: Endmember matrix $\boldsymbol{A}$, abundance matrix $\boldsymbol{S}$.  

Pre-processing: Extracting master data points set $\boldsymbol{Y}$ by PPI. Then, compute the sparse regular term parameter $\lambda$.

Initialize $\boldsymbol{C}$ by SNPA [19] and $S$ by applying the updating rule $B$.  

repeat
    1. Update $\boldsymbol{C}$ and $\boldsymbol{\alpha}$ by the updating rule $A$.
    2. Update $\boldsymbol{S}$ by applying the updating rule B. 
until Converge.

3. Calculate endmembers matrix $\boldsymbol{A}: \boldsymbol{A}=\boldsymbol{Y} \boldsymbol{C}$.

以上就是基于原型分析的高光谱图像处理算法的主要数学原理和求解过程。该算法具有很强的鲁棒性,可以有效处理噪声和异常值,在高光谱图像融合和解混等任务中表现出色。