基于 N-FINDER 的端元提取算法

N-FINDER 算法是一种用于从高光谱数据立方体中提取端元特征的算法。该算法的核心思想是在数据集中寻找表示最大体积的单纯形,其顶点即为所需的端元特征。算法的数学原理和具体步骤如下:

算法数学原理

设有 \(M \times N \times C\) 维的高光谱数据立方体 \(\mathbf{X}\),我们需要从中找出 \(p\) 个端元特征 \(\mathbf{e}_1, \mathbf{e}_2, \ldots, \mathbf{e}_p\),这些特征向量组成的矩阵记为 \(\mathbf{E} = [\mathbf{e}_1, \mathbf{e}_2, \ldots, \mathbf{e}_p]^T \in \mathbb{R}^{C \times p}\)。

根据线性混合模型,数据立方体 \(\mathbf{X}\) 中的每个像素 \(\mathbf{x}_i\) 可表示为端元特征的线性组合:

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

其中 \(a_{ij}\) 为相应的丰度分数,满足约束条件 \(\sum_{j=1}^p a_{ij} = 1\) 和 \(a_{ij} \geq 0\),\(\mathbf{n}_i\) 为加性噪声项。

N-FINDER 算法的目标是找到一组端元特征 \(\mathbf{E}\),使得由这些特征向量构成的 \(p\) 维单纯形具有最大体积。该单纯形的体积可由端元特征矩阵 \(\mathbf{E}\) 的行列式的绝对值表示:

\[\mathrm{vol}(\mathbf{E}) = |\det(\mathbf{E})|\]

算法通过迭代不断更新端元特征矩阵 \(\mathbf{E}\),直到找到最大化上述行列式的解。

算法步骤

1. 初始化: 从数据立方体中随机选择 \(p\) 个像素作为初始端元特征向量,构成初始端元特征矩阵 \(\mathbf{E}^{(0)}\),计算其行列式的绝对值 \(\mathrm{vol}(\mathbf{E}^{(0)})\)。
2. 迭代:
   • 对每个端元特征向量 \(\mathbf{e}_i^{(t)}\),用数据立方体中的其他像素 \(\mathbf{x}_j\) 代替,构成新的端元特征矩阵 \(\mathbf{E}_{ij}^{(t+1)}\)。
   • 计算 \(\mathrm{vol}(\mathbf{E}_{ij}^{(t+1)})\),若其大于 \(\mathrm{vol}(\mathbf{E}^{(t)})\),则更新 \(\mathbf{E}^{(t+1)} = \mathbf{E}_{ij}^{(t+1)}\)。
   • 重复上述步骤,直到无法再提高体积为止。
3. 输出: 得到的最终端元特征矩阵 \(\mathbf{E}\) 即为所求解。

代码实现

function endmemSig = nfindr(cube, num, option)
% NFINDR 通过在数据立方体上找到单纯形来提取端元特征
%   ENDMEMSIG = NFINDR(CUBE, NUM) 通过在降维后的数据上找到单纯形的体积来估计端元特征
%   NUM 描述了 CUBE 中存在的端元特征数量
%   CUBE 是一个包含高光谱数据或高光谱立方体对象的 M-by-N-by-C 数值数组
%   M 和 N 分别是高光谱立方体的行数和列数
%   C 描述了高光谱立方体中的波段数量
%   NUM 是一个大于 1 的正整数值
%   输出 ENDMEMSIG 是一个 C-by-NUM 的数值矩阵,它提供了端元特征

%   ENDMEMSIG = NFINDR(___, Name, Value) 使用选定的名称-值对在降维后的数据上找到最大体积单纯形来估计端元特征

%   'NumIterations'  算法运行的次数,如果没有收敛,默认值为: 3*NUM

%   'ReductionMethod' - 用于执行维数减少的方法。选择如下:
%
%       'PCA'  基于最大数据变化,主成分分析 (PCA) 将大量数据转换为较小的数据
%              同时保持大部分信息
%
%       'MNF'  基于噪声分数,最大噪声分数 (MNF) 将大量数据转换为较小的数据
%              同时保持大部分信息。MNF 是一种根据图像质量对分量进行排序的方法
%              默认方法是 MNF

%   注释
%   -----
%   1) NFINDR 算法使用默认的随机种子来生成随机单位向量
%
%   2) 对于较大数据集,NFINDR 算法计算量很大

%   类支持
%   -------------
%   输入数组 CUBE 必须是以下类之一: uint8、uint16、uint32、int8、int16、int32、single、double、uint64、int64 或 hypercube 对象
%   它必须是实数且非稀疏矩阵
%   NUM 和 NUMITER 是正整数标量
%   输出 ENDMEMSIG 是与输入立方体相同类的 C-by-NUM 数值矩阵

%   参考
%   ---------
%   M.E. Winter, "N-FINDR: an algorithm for fast autonomous spectral
%   end-member determination in hyperspectral data", Proc. SPIE 3753,
%   Imaging Spectrometry V, 1999.

%   示例 - 1 : 找到端元特征
%   -----------
%   % 加载 indian pines 数据集
%   cube = load('indian_pines.mat');
%
%   % 估计端元数量
%   num = countEndmembersHFC(cube.indian_pines);
% 
%   % 估计并可视化端元特征
%   endmemSig = nfindr(cube.indian_pines, num);
%   plot(endmemSig);
%   xlabel('Band index')
%   ylabel('Signatures')
%   title('Endmembers')

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

%   Copyright 2020-2021 The MathWorks, Inc.

arguments
    cube {mustBeA(cube,'hypercube')} % 强制 cube 是 hypercube 对象
     num (1,1) {mustBeNumeric, mustBeNonNan, mustBeFinite, mustBeNonsparse, ...
        mustBeNonempty, mustBeReal, mustBeInteger, mustBePositive, ...
        mustBeLessSize(num, cube)} % 强制 num 是正整数且小于波段数
    option.NumIterations (1,1) {mustBeNumeric, mustBeNonNan, mustBeFinite, ...
        mustBeNonsparse, mustBeNonempty, mustBeReal, mustBeInteger,...
        mustBePositive} = 3*num % 默认迭代次数为 3*num
    option.ReductionMethod (1,1) string ...
        {validatestring(option.ReductionMethod,{'MNF', 'PCA'})} = "MNF"  % 默认降维方法为 MNF
end

% 用于编译器测试
if isdeployed
    rootDir = ctfroot;
else
    rootDir = matlabroot;
end
% 注册资源/hyperspectral
matlab.internal.msgcat.setAdditionalResourceLocation(rootDir);

% 输入解析
if isobject(cube)
    cube = cube.DataCube; % 如果 cube 是对象,则提取数据立方体
end

maxItr = option.NumIterations; % 最大迭代次数
method = validatestring(option.ReductionMethod, {'MNF', 'PCA'}); % 降维方法

% 立方体维度
[rows, cols, channels] = size(cube);
numOfSamples = rows*cols; % 样本数量

% 重塑原始数据立方体
originalCube = reshape(cube, numOfSamples, channels);

if isinteger(cube)
    cube = single(cube); % 将整数类型转换为 single
    num = single(num);
end

if strcmpi(method, 'PCA')
    % 对高光谱立方体进行主成分分析 (PCA)
    cube = hyperpca(cube, num-1);
    
else
    % 使用最大噪声分数 (MNF) 变换来减少计算要求
    cube = hypermnf(cube, num-1);
end

% 重塑降维后的立方体
cube = reshape(cube, [numOfSamples num-1]);

% 初始化矩阵,并将第一行设置为 1,每一列向量包含端元的光谱
signatures = zeros(num);
signatures(1,:) = 1;

% 在 1 到 numOfSamples 范围内生成随机样本
prevS = rng('default');
randSamplesInd = randi([1 numOfSamples], 1, num);

% 清理随机状态
cleanup = onCleanup(@() rng(prevS));
% 选择随机像素并将其设置为端元
signatures(2:num, :) = cube(randSamplesInd, :)';

% 对于该算法,除以 (factorial(p-1)) 是一个不变量,因此跳过了
actualVolume = abs(det(signatures));

% 初始化迭代次数
itr = 1;

% 初始化体积 v1 和 v2
initVol = -itr;
maxVol = actualVolume;

% 计算所有像素的体积,如果替换像素会导致体积增加,则将替换像素视为端元
% 该过程重复进行,直到没有更多端元替换
while maxVol > initVol
    for sigCount=1:num
         
         % 找到体积的单纯形
         [randSamplesInd, actualVolume, signatures] = findSimplex(signatures, cube, ...
             randSamplesInd, actualVolume, numOfSamples, num, sigCount); 
    end
    
    % 如果达到最大迭代次数,则中断
    if itr == maxItr
        break;
    end
    
    % 更新迭代次数、v1 和 v2
    itr = itr+1;
    initVol = maxVol;
    maxVol = actualVolume;
end

% 根据样本索引从体积中提取端元特征
endmemSig = originalCube(randSamplesInd, :)';
end

function [randSamplesInd, actualVolume, signatures] = findSimplex(signatures, reducedCube, ...
    randSamplesInd, actualVolume, numOfSamples, num, sigCount)
for sampleCount=1:numOfSamples
    
    % 检查降维后的每个样本
    signatures(2:num, sigCount) = reducedCube(sampleCount, :);

    % 使用端元估计形成的单纯形的体积与特征矩阵的行列式成正比
    volume = abs(det(signatures));
    if volume > actualVolume
        actualVolume = volume;
        randSamplesInd(sigCount) = sampleCount;
    end
end

% 用最新样本更新端元特征
signatures(2:num, sigCount) = reducedCube(randSamplesInd(sigCount), :);

end

function mustBeA(hypercube,classIn)

if isobject(hypercube)
    validateattributes(hypercube,{classIn},{'nonempty'},'nfindr','hypercube');
    cube = hypercube.DataCube;
else
    cube = hypercube;
end

validateattributes(cube, ...
    {'numeric'}, {'nonsparse', 'nonempty', 'real','nonnan', 'finite', ...
    'ndims',3}, mfilename, 'hypercube', 1);

end

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

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

该算法的优点是思路清晰,能够有效提取出最大体积簇对应的端元特征。但缺点是计算量很大,对于大数据集会耗费大量时间。因此在实际应用中,通常会先利用 PCA 或 MNF 等方法对数据进行降维,再在降维后的数据上运行 N-FINDER 算法,以缩减计算量。

总的来说,N-FINDER 算法通过最大化端元特征向量组成的单纯形体积,来有效提取出数据中的端元特征,是一种较为经典的无监督端元提取算法。