改进的快速迭代纯像素指数提取端元特征

FIPPI算法简介

改进的快速迭代纯像素指数(FIPPI)算法是一种从高光谱遥感数据中提取端元谱签的有效方法。它是基于纯像素指数(PPI)算法的迭代改进版本,旨在通过合理的终止条件和精心设计的迭代规则,提高计算效率并获得更好的端元估计结果。

算法数学原理

降维

在FIPPI算法中,首先对原始高维数据立方体进行降维,以降低计算复杂度。常用的降维方法有:

1. 主成分分析(PCA):将高维数据投影到主成分空间,保留主要信息。
2. 最大噪声分数变换(MNF):根据图像质量对分量进行排序,保留高质量分量。

对于原始数据矩阵 \(\boldsymbol{X} \in \mathbb{R}^{m \times n}\),经过降维后得到 \(\boldsymbol{X}^{*} \in \mathbb{R}^{m \times p}\),其中 \(p < n\)。

自动目标生成过程(ATGP)

ATGP算法用于根据正交子空间投影估计初始端元集合。给定降维后的数据 \(\boldsymbol{X}^{*}\),包含 \(p\) 个特征向量,ATGP算法如下:

1. 找到第一个端元向量索引 \(i_1\),使 \(\sum_{j=1}^{p} (\boldsymbol{x}_j^*)^2\) 最大,即 \(\boldsymbol{u}_1 = \boldsymbol{x}_{i_1}^*\)。
2. 对于第 \(k(k \geq 2)\) 个端元向量,求解: \[\boldsymbol{u}_k = \arg\max_{\boldsymbol{x}_j^*} \left\|\left(\boldsymbol{I}_p - \boldsymbol{U}_{k-1}\boldsymbol{U}_{k-1}^+\right)\boldsymbol{x}_j^*\right\|^2\] 其中 \(\boldsymbol{U}_{k-1} = \begin{bmatrix}\boldsymbol{u}_1 & \boldsymbol{u}_2 & \cdots & \boldsymbol{u}_{k-1}\end{bmatrix}\), \(\boldsymbol{U}_{k-1}^+\) 为 \(\boldsymbol{U}_{k-1}\) 的伪逆。

算法重复上述过程直至获得 \(p\) 个初始端元向量。

FIPPI迭代规则

使用ATGP算法获得初始端元集合 \(\boldsymbol{U} \in \mathbb{R}^{p \times p}\) 后,FIPPI算法按照如下迭代规则进行优化:

1. 将降维数据投影到当前端元空间: \(\boldsymbol{Y} = \boldsymbol{U}^T\boldsymbol{X}^*\)
2. 对每个样本,找到其在投影空间中的最大绝对投影值及其对应的样本索引 \(i_j\): \[i_j = \arg\max_{1 \leq i \leq m} \left|\boldsymbol{y}_i\right|\]
3. 构建极值集合 \(\boldsymbol{X}^{ext} = \left\{\boldsymbol{x}_{i_j}^*\right\}_{j=1}^q\),其中 \(q\) 为极值数量。
4. 检查终止条件: 如果 \(\boldsymbol{X}^{ext} \subseteq \boldsymbol{U}\),则终止迭代;否则将 \(\boldsymbol{X}^{ext}\) 中不在 \(\boldsymbol{U}\) 的向量添加到 \(\boldsymbol{U}\),并返回步骤1继续迭代。

经过多次迭代,算法最终输出 \(\boldsymbol{U}\) 作为估计的端元集合。

MATLAB 实现

以下是FIPPI算法在MATLAB中的实现,详细说明请参考代码注释:

function endMembers = fippi(cube, num, option)
% fippi 函数用于从高光谱数据立方体中提取端元谱签,使用快速迭代纯像素指数算法。

%   ENDMEMBERS = FIPPI(CUBE, NUM) 通过将降维后的数据投影到初始估计的端元上,来估计端元谱签。
%   NUM 表示存在于 CUBE 中的端元谱签数量。CUBE 是一个包含高光谱数据的 M-by-N-by-C 维数值数组或者
%   一个高光谱数据立方体对象。M 和 N 分别表示高光谱数据立方体的行和列。C 表示高光谱数据立方体的波段数。
%   NUM 是一个非零正整数值。输出 ENDMEMBERS 是一个数值矩阵,它提供了 C-by-K 维的端元谱签,其中 K 是
%   根据终止条件估计得到的最终端元集合的大小。

%   ENDMEMBERS = FIPPI(___,Name, Value) 通过将降维后的数据投影到初始估计的端元上来估计端元谱签。
%   数据降维使用了两种方法,参数如下:
%
%   'ReductionMethod' - 用于执行维数降低的方法,可选如下:
%
%       'PCA' 基于最大数据方差的主成分分析(PCA)方法,通过保留大部分信息,
%             将大规模数据转换为较小规模。
%
%       'MNF' 基于噪声分数的最大噪声分数(MNF)变换方法,通过保留大部分信息,
%             将大规模数据转换为较小规模。MNF是一种根据图像质量对分量进行排序的方法。
%             默认使用 MNF 方法。

%   注意事项
%   -----
%   1) FIPPI 是一种迭代算法,其中迭代规则被设计用于改进每次迭代,直到达到最终的端元集合。
%   2) FIPPI 的输出不一定是 NUM 个端元,它会根据 FIPPI 的终止条件输出 K 个端元。
%   FIPPI 有一个替换规则用于迭代地细化每次迭代,并在最终的探针集合不再包含新的端元时终止。
%   3) NUM 输入参数有助于使用自动目标生成过程(ATGP)算法估计初始端元集合。

%   类支持
%   -------------
%   输入数组 CUBE 必须是以下类型之一: uint8, uint16, uint32, int8, int16,
%   int32, single, double, uint64, int64 或者是一个高光谱数据立方体对象。
%   它必须是实数且非稀疏矩阵。输出 ENDMEMBERS 是一个 C-by-K 维的数值矩阵。

%   参考文献
%   ----------
%   CI Chang, A. Plaza, "A fast iterative algorithm for implementation of
%   pixel purity index", IEEE Geoscience and Remote Sensing Letters,
%   volume-3, pages: 63-67, 2006.

%   示例 1: 查找端元谱签
%   -----------
%   % 加载 indianpines 数据集
%   cube = load('indian_pines.mat');
%
%   % 估计端元数量
%   num = countEndmembersHFC(cube.indian_pines);
%
%   % 使用 FIPPI 算法估计并可视化端元谱签
%   endMembers = fippi(cube.indian_pines, num);
%   plot(endMembers);
%   xlabel('Band index')
%   ylabel('Signatures')
%   title('Endmembers')

%   示例 2: 使用 'PCA' 维数降低方法查找端元谱签
%   -----------
%   % 读取高光谱数据立方体
%   obj = hypercube('indian_pines.hdr');
%
%   % 估计端元数量
%   num = countEndmembersHFC(obj);
%
%   % 使用 FIPPI 算法估计端元谱签
%   endMembers = fippi(obj, num, 'ReductionMethod', 'pca');
%   plot(endMembers);
%   xlabel('Band index')
%   ylabel('Signatures')
%   title('Endmembers')
%
%   另请参阅 ppi, nfindr, hypercube, countEndmembersHFC。

%   版权所有 2020-2021 The MathWorks, Inc.

arguments
    cube {mustBeA(cube,'hypercube')}
    % cube 必须是一个高光谱数据立方体
    num (1,1) {mustBeNumeric, mustBeNonNan, mustBeFinite, mustBeNonsparse, ...
        mustBeNonempty, mustBeReal, mustBeInteger, mustBePositive, ...
        mustBeLessSize(num, cube)}
    % num 必须是一个正整数,且小于数据立方体的波段数
    option.ReductionMethod (1,1) string ...
        {validatestring(option.ReductionMethod,{'MNF', 'PCA'})} = "MNF"
    % option.ReductionMethod 是一个字符串,用于指定维数降低方法,默认为 'MNF'
end

% 编译测试所需
if isdeployed
    rootDir = ctfroot;
else
    rootDir = matlabroot;
end
% 注册资源/高光谱数据
matlab.internal.msgcat.setAdditionalResourceLocation(rootDir);

% 输入解析
if isobject(cube)
    cube = cube.DataCube;
end

method = validatestring(option.ReductionMethod, {'MNF', 'PCA'});

% 获取数据立方体的大小
[rows, cols, channels] = size(cube);
numofPixels = rows*cols;

% 重塑数据立方体用于提取谱签
originalCube = reshape(cube, numofPixels, channels);

if isinteger(cube)
    cube = single(cube);
    num = single(num);
end

if strcmpi(method, 'PCA')
    % 对高光谱数据立方体执行主成分分析 (PCA)
    cube = hyperpca(cube, num);
    
else 
    % 最大噪声分数 (MNF) 变换用于确定图像数据的内在维数,去除数据中的噪声,并降低计算需求
    cube = hypermnf(cube, num);
end

% 使用 ATGP 算法初始化探针
initializeSkewers = atgp(cube,numofPixels, num);

% FIPPI 算法
stoppingCriteria = false;
cube = reshape(cube, numofPixels, num);

% 终止条件用于寻找最终端元
while ~stoppingCriteria
    
    % 将降维后的数据投影到初始估计的端元上
    projPoints = abs(initializeSkewers*cube');
    
    % 找到极值位置以形成极值集合
    [~, maxLoc] = max(projPoints,[], 2);
    maxLoc = unique(maxLoc);
    
    % 检查新端元集合
    stoppingCriteria = all(ismember(cube(maxLoc,:), initializeSkewers), 'all');
    initializeSkewers = unique([cube(maxLoc, :);initializeSkewers],'sorted','rows');
end

endMembers = originalCube(maxLoc, :)';
end

function U = atgp(dimRed, numVariables, num)
% 自动目标生成过程 (ATGP) 算法基于正交子空间投影估计初始签名集合

% 重塑降维后的数据立方体
dimRed = reshape(dimRed, numVariables, num)';

% 找到第一个签名索引
[~, idx] = max(sum(dimRed.^2));

% 初始化并获取第一个端元签名
U = zeros(num, class(dimRed));
U(:,1) = dimRed(:,idx);

for n=1:num-1
    % 最大绝对正交子空间投影
    [~, idx] = max(sum(((eye(num) - (U * pinv(U)))*dimRed).^2)); 
    U(:,n+1) = dimRed(:,idx);
end
end

function mustBeA(hypercube,classIn)

if isobject(hypercube)
    validateattributes(hypercube,{classIn},{'nonempty'},'fippi','hypercube');
else
    validateattributes(hypercube, ...
    {'uint8','int8','uint16','int16','uint32','int32','single','double', 'uint64', 'int64'},...
    {'nonsparse', 'nonempty', 'real','nonnan', 'finite', 'ndims',3}, mfilename, 'hypercube', 1);
end

end

% 验证 numComponents 小于波段数
function mustBeLessSize(a, b)
if isobject(b)
    cube = b.DataCube;
else
    cube = b;
end
sz = size(cube);
validateattributes(a, {'numeric'},...
    {'>', 0, '<=',sz(3),'scalar'},mfilename, 'num',2);

if (a >= sz(1)*sz(2))
    error(message('hyperspectral:hyperpca:notGreaterEqual'));    
end
end

使用示例:

% 载入数据
cube = load('indian_pines.mat');

% 估计端元数量
num = countEndmembersHFC(cube.indian_pines);

% 使用FIPPI算法估计端元谱签
endMembers = fippi(cube.indian_pines, num);

% 可视化端元谱签
plot(endMembers);
xlabel('Band index');
ylabel('Signatures');
title('Endmembers');

总结

FIPPI算法通过合理的数学模型和精心设计的迭代规则,有效地从高光谱遥感数据中提取端元谱签。它克服了传统PPI算法的一些缺陷,提高了计算效率和端元估计精度,是遥感端元提取的有力工具。