顶点成分分析算法(VCA)的原理与求解

顶点成分分析算法(Vertex Component Analysis, VCA)是一种快速无监督提取高光谱图像端元的算法。该算法的主要思想是利用端元在高维特征空间中形成数据云的顶点的性质,通过迭代寻找位于数据云边界的端元,从而实现端元提取。下面将详细介绍该算法的数学原理和求解过程。

问题描述

设有一幅 \(L\) 波段的高光谱图像,每个像素的光谱向量为 \(\mathbf{r} \in \mathbb{R}^{L}\),整个图像可表示为 \(\mathbf{R} = [\mathbf{r}_1, \mathbf{r}_2, \dots, \mathbf{r}_N] \in \mathbb{R}^{L \times N}\),其中 \(N\) 为像素个数。根据线性混合模型,每个光谱向量 \(\mathbf{r}_i\) 可表示为 \(p\) 个端元 \(\mathbf{m}_1, \mathbf{m}_2, \dots, \mathbf{m}_p \in \mathbb{R}^{L}\) 的线性组合:

\[\mathbf{r}_i = \sum_{j=1}^{p} a_{ij} \mathbf{m}_j + \mathbf{n}_i\]

其中 \(a_{ij} \geq 0\) 为满足条件 \(\sum_{j=1}^{p} a_{ij} = 1\) 的端元丰度系数, \(\mathbf{n}_i\) 为高斯白噪声。将上式写成矩阵形式:

\[\mathbf{R} = \mathbf{M} \mathbf{A} + \mathbf{N}\]

其中 \(\mathbf{M} = [\mathbf{m}_1, \mathbf{m}_2, \dots, \mathbf{m}_p] \in \mathbb{R}^{L \times p}\) 为端元矩阵, \(\mathbf{A} \in \mathbb{R}^{p \times N}\) 为丰度矩阵, \(\mathbf{N} \in \mathbb{R}^{L \times N}\) 为噪声矩阵。

VCA算法的目标是从高光谱数据 \(\mathbf{R}\) 中无监督地提取出端元矩阵 \(\mathbf{M}\)。

VCA算法原理

数据云顶点性质

假设存在 \(p\) 个端元, 那么所有光谱向量都落在由这 \(p\) 个端元张成的 \(p-1\) 维单纯形内。在没有噪声的情况下,位于单纯形边界的光谱向量必定是某个端元。因此,VCA算法的核心思想是寻找位于数据云边界的光谱向量,将其作为提取的端元。

投影方法

为了寻找边界光谱向量, VCA算法利用投影技术将高维数据投影到低维空间,从而使边界光谱向量的几何特征更加明显。

首先利用奇异值分解(SVD)计算 \(p\) 阶投影矩阵 \(\mathbf{U}_d \in \mathbb{R}^{L \times d}\) (\(d=p\) 或 \(d=p-1\),取决于信噪比SNR),并将数据 \(\mathbf{R}\) 投影到 \(d\) 维空间:

\[\mathbf{X} = \mathbf{U}_d^T \mathbf{R}\]

其中 \(\mathbf{X} \in \mathbb{R}^{d \times N}\) 为投影后的数据。

接下来需要从 \(\mathbf{X}\) 中找到边界光谱向量。观察到在 \(p\) 维单纯形中,点到边界的距离大于点到内部的距离。因此可以通过最大化投影矩阵中的行向量 \(\mathbf{x}_i^T \in \mathbb{R}^{1 \times N}\) 与所有点的夹角余弦值,从而找到距离最远的边界点。

\[\max_{\|\mathbf{y}\|=1} \frac{\|\mathbf{y}^T \mathbf{x}_i\|}{\sqrt{\sum_{j=1}^{N} (\mathbf{y}^T \mathbf{x}_j)^2}}\]

根据Cauchy-Schwarz不等式,上式可化简为求 \(\mathbf{x}_i\) 的最大绝对值分量:

\[\max_{1 \leq j \leq N} |\mathbf{x}_{ij}|\]

因此,VCA算法采取以下迭代方式求边界光谱向量:

1. 对任意非零向量 \(\mathbf{w} \in \mathbb{R}^{p}\),令 \(\mathbf{f} = \mathbf{w} - \mathbf{A} (\mathbf{A}^T \mathbf{A})^{-1} \mathbf{A}^T \mathbf{w}\) 为 \(\mathbf{w}\) 在当前求解向量的补空间的投影,并归一化 \(\mathbf{f}\)。
2. 对归一化后的 \(\mathbf{f}\),计算 \(\mathbf{v} = \mathbf{f}^T \mathbf{Y}\),其中 \(\mathbf{Y} = [\mathbf{X}; \mathbf{c}]\) (\(\mathbf{c}\) 为尺度常数)。
3. 找到 \(\mathbf{v}\) 的最大绝对值分量对应的下标 \(k\),将 \(\mathbf{A}\) 的第 \(i\) 列赋值为 \(\mathbf{Y}\) 的第 \(k\) 列。

重复以上步骤,最终可获得 \(p\) 个顶点向量,即提取的端元矩阵 \(\mathbf{U}\)。

算法求解过程

根据上述原理,VCA算法的具体求解步骤如下:

1. 计算数据均值向量 \(\overline{\mathbf{r}} = \frac{1}{N} \sum_{i=1}^{N} \mathbf{r}_i\),并对数据进行中心化: \(\widetilde{\mathbf{R}} = \mathbf{R} - \overline{\mathbf{r}} \mathbf{1}_N^T\)。

2. 计算噪声协方差矩阵的 \(p\) 阶投影矩阵 \(\mathbf{U}_d\):

\[\mathbf{U}_d = \mathrm{svd}(\widetilde{\mathbf{R}} \widetilde{\mathbf{R}}^T / N, p)\]

3. 将数据投影到 \(d\) 维空间: \(\mathbf{X} = \mathbf{U}_d^T \widetilde{\mathbf{R}}\)。

4. 根据信噪比 SNR 判断是否需要尺度约束:

   • 若 \(\mathrm{SNR} > \mathrm{SNR}_\mathrm{th}\),则 \(\mathbf{Y} = \mathbf{X}\)。
   • 若 \(\mathrm{SNR} \leq \mathrm{SNR}_\mathrm{th}\),则 \(\mathbf{Y} = [\mathbf{X}; c]\),其中 \(c = \max_i \|\mathbf{x}_i\|_2\)。

5. 初始化 \(p \times p\) 矩阵 \(\mathbf{A}\):

   • 第一列赋值为 \([1, 0, \dots, 0]^T\)。
   • 其余列置零。

6. 对每个端元进行迭代:

   • 生成随机向量 \(\mathbf{w}\)。
   • 计算 \(\mathbf{f} = (\mathbf{I} - \mathbf{A} (\mathbf{A}^T \mathbf{A})^{-1} \mathbf{A}^T) \mathbf{w}\),并归一化 \(\mathbf{f}\)。
   • 计算 \(\mathbf{v} = \mathbf{f}^T \mathbf{Y}\),找到 \(\mathbf{v}\) 的最大绝对值分量对应的下标 \(k\)。
   • 将 \(\mathbf{A}\) 的当前列赋值为 \(\mathbf{Y}\) 的第 \(k\) 列。

7. 估计端元矩阵 \(\mathbf{U}\):

   • 若 \(\mathrm{SNR} > \mathrm{SNR}_\mathrm{th}\),则 \(\mathbf{U} = \mathbf{U}_d \mathbf{X}_\mathrm{indicies}\)。
   • 若 \(\mathrm{SNR} \leq \mathrm{SNR}_\mathrm{th}\),则 \(\mathbf{U} = \mathbf{U}_d \mathbf{X}_\mathrm{indicies} + \overline{\mathbf{r}} \mathbf{1}_p^T\)。

其中 \(\mathbf{X}_\mathrm{indicies}\) 为由提取端元对应的列向量组成的矩阵。

以上就是VCA算法的完整求解过程。该算法巧妙地利用了数据投影和顶点特性,通过迭代方式无监督地提取出高光谱图像的端元,具有很好的效率和稳健性。

代码实现

function [ U, indicies ] = vca( R, p ) %-------------------------------------------------------------------------- 
% 顶点成分分析(Vertex Component Analysis)算法
%
% 用法
% [ U, indicies ] = vca( R, p )
% 输入
% R : 高光谱图像数据(bands,pixels)  
% p : 端元个数
% 输出
% U : 端元矩阵(bands,p)
% indicies : 端元在R中的索引
%
% 参考文献
% J. M. P. Nascimento and J. M. B. Dias, "Vertex component analysis: A
% fast algorithm to unmix hyperspectral data," IEEE Transactions on
% Geoscience and Remote Sensing, vol. 43, no. 4, pp. 898 - 910, Apr. 2005.
%-------------------------------------------------------------------------- 
[L, N]=size(R); % 获取图像波段数L和像素数N
r_m = mean(R,2); % 计算每个波段的均值
R_o = R - repmat(r_m,[1 N]); % 对数据进行去均值
[Ud,~,~] = svds(R_o*R_o'/N,p); % 计算p主成分投影矩阵Ud
x_p = Ud' * R_o; % 将去均值数据投影到p维空间
P_y = sum(R(:).^2)/N; % 计算原始数据的能量
P_x = sum(x_p(:).^2)/N + r_m'*r_m; % 计算投影数据的能量
SNR = abs(10*log10( (P_x - p/L*P_y)/(P_y- P_x) )); % 估计信噪比SNR

%fprintf('SNR estimate [dB]: %g\n', SNR); % 输出估计的SNR(注释掉)

% 确定使用哪种投影
SNRth = 15 + 10*log(p) + 8; % 设置SNR阈值 
%SNRth = 15 + 10*log(p); % 原论文中提出的阈值(注释掉)

if (SNR > SNRth) % 若SNR高于阈值
    d = p; 
    [Ud,~,~] = svds(R*R'/N,d); % 计算d主成分投影矩阵Ud
    X = Ud'*R; % 将原始数据投影到d维空间
    u = mean(X, 2); % 计算每个投影分量的均值u
    Y = X ./ repmat( sum( X .* repmat(u,[1 N]) ) ,[d 1]); % 对X进行归一化得到Y
else % 若SNR低于阈值
    d = p-1;
    r_m = mean(R,2); % 重新计算均值
    R_m = repmat(r_m,[1 N]); % 生成均值矩阵 
    R_o = R - R_m; % 对数据去均值
    [Ud,~,~] = svds(R_o*R_o'/N,d); % 计算d主成分投影矩阵
    X = Ud'*R_o; % 将去均值数据投影到d维空间
    X = X(1:d,:); % 只保留前d行
    c = max(sum(X.^2,1))^0.5; % 计算X每列的L2范数的最大值c  
    Y = [X ; c*ones(1,N)] ; % 将X和c组合形成Y
end

e_u = zeros(p, 1); e_u(p) = 1; % 生成单位矩阵的最后一列
A = zeros(p, p); A(:, 1) = e_u; % 初始化矩阵A
indicies = zeros(1,p); % 初始化indicies

% 使用投影搜索找到p个顶点
for i=1:p
    w = rand(p, 1); % 生成随机向量w
    f = w - A*pinv(A)*w; % 计算w在A列空间的正交分量f
    f = f / sqrt(sum(f.^2)); % 将f归一化
    v = f'*Y; % 将Y投影到f方向
    [~, indicies(i)] = max(abs(v)); % 找到投影的最大绝对值对应的索引
    A(:,i) = Y(:,indicies(i)); % 更新A的第i列
end

% 根据顶点索引计算端元矩阵U
if (SNR > SNRth)
    U = Ud*X(:,indicies); % 若SNR高,直接利用投影数据
else
    U = Ud*X(:,indicies) + repmat(r_m,[1 p]); % 若SNR低,加上均值
end
end