基于纯像元指数(PPI)的端元提取算法

PPI算法用于从高光谱数据立方体中提取端元光谱。该算法通过在降维后的数据上投影随机生成的单位向量,来识别出表现出最大投影值的端元像元。

算法原理

1. 对高光谱数据立方体进行降维处理,可选择PCA(主成分分析)、MNF(最大噪声分数)或不降维三种方式。
2. 生成一定数量的随机单位向量(默认10000个)。
3. 将降维后的数据投影到每个随机单位向量上,记录每个像元点获得的最大投影值。
4. 根据最大投影值对像元点排序,选取投影值最大的N个像元点作为端元光谱。

算法实现

这里介绍如何使用MATLAB实现PPI算法。

  function endmemSig = ppi(cube, num, option)
%PPI 使用像素纯度指数提取端元签名。
%   ENDMEMSIG = PPI(CUBE, NUM) 通过在随机生成的单位向量上投影降维数据来估计端元签名。
%   NUM 描述了CUBE中存在的端元签名数量。CUBE是一个包含高光谱数据或hypercube对象的
%   M-by-N-by-C的数值数组。M和N分别是高光谱立方体的行和列，C是高光谱立方体中的波段数量。
%   ENDMEMSIG是一个数值矩阵，提供了维度为C-by-NUM的端元签名。
%
%   ENDMEMSIG = PPI(___, Name, Value) 通过在选定NV对的降维数据上投影随机生成的单位向量来
%   估计端元签名。
%
%   'NumVectors'    随机生成的单位向量数量。随着数量增加，算法在找到端元的同时也增加了计算复杂度。
%                   默认值为10^4。
%
%   'ReductionMethod': 执行降维的方法。选项包括：
%
%       'PCA'  基于最大数据变异量的主成分分析（PCA），通过保留大部分信息将大数据集转换为较小的数据集。
%
%       'MNF'  基于噪声分数的最大噪声分数（MNF），同样是通过保留大部分信息来缩减数据集大小。MNF是一种根据图像质量排序组件的方法。默认方法为MNF。
%
%       'NONE' 在随机生成的单位向量上投影原始数据。
%
%   注：
%   1) PPI算法使用默认的随机种子生成随机单位向量。
%
%   类支持：
%   输入数组CUBE必须是以下类之一：uint8, uint16, uint32, int8, int16, int32, single, double, uint64, int64或hypercube对象。
%   它必须是实数且非稀疏的。NUM和NumVectors是正整数值的数值标量。输出ENDMEMSIG是与输入cube同类的数值矩阵，大小为C-by-NUM。
%
%   参考：
%   J.W  Boardman, F.A. Kruse 和 R.O. Green, "通过部分解混AVIRIS数据映射目标签名", 技术报告, 加利福尼亚, 美国, 1995。
%
%   示例 - 1 : 寻找端元签名。
%   % 加载indian pines数据集
%   cube = load('indian_pines.mat');
%
%   % 估计端元数量。
%   num = countEndmembersHFC(cube.indian_pines);
%
%   % 估计并可视化端元签名
%   endmemSig = ppi(cube.indian_pines, num);
%   plot(endmemSig);
%   xlabel('波段索引')
%   ylabel('签名')
%   title('端元')
%
%   示例 - 2 : 使用'PCA'降维方法寻找端元签名。
%   % 加载高光谱立方体。
%   obj = hypercube('paviaU.hdr');
%
%   % 估计端元数量。
%   num = countEndmembersHFC(obj);
%
%   % 使用PPI估计并可视化端元签名。
%   endmemSig = ppi(obj, num, 'ReductionMethod', 'pca', 'NumVectors', 10^3);
%   plot(endmemSig);
%   xlabel('波段索引')
%   ylabel('签名')
%   title('端元')



arguments
    cube {mustBeA(cube,'hypercube')}
    num (1,1) {mustBeNumeric, mustBeNonNan, mustBeFinite, mustBeNonsparse, ...
        mustBeNonempty, mustBeReal, mustBeInteger, mustBePositive, ...
        mustBeLessSize(num, cube)}
    option.NumVectors (1,1) {mustBeNumeric, mustBeNonNan, mustBeFinite, ...
        mustBeNonsparse, mustBeNonempty, mustBeReal, mustBeInteger,...
        mustBePositive} = 10^4
    option.ReductionMethod (1,1) string ...
        {validatestring(option.ReductionMethod,{'MNF', 'PCA', 'NONE'})} = "MNF"

end

% For compiler tests
if isdeployed
    rootDir = ctfroot;
else
    rootDir = matlabroot;
end
% Register resources/hyperspectral.
matlab.internal.msgcat.setAdditionalResourceLocation(rootDir);

% Input parsing.
if isobject(cube)
    cube = cube.DataCube;
end

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

% Volume dimensions.
[rows, cols, channels] = size(cube);
numOfSamples = rows*cols;

originalCube = reshape(cube, numOfSamples, channels);

% Cast the input data to single data type
if isinteger(cube)
    cube = single(cube);
    num = single(num);
end

if strcmpi(method, 'PCA')
    % Principal Component Analysis (PCA) on hyperspectral cube.
    cube = reshape(hyperpca(cube, num),numOfSamples, num);

elseif strcmpi(method, 'MNF')
    % Maximum noise fraction (MNF) transformation is used to reduce the
    % computational requirement.
    cube = reshape(hypermnf(cube, num), numOfSamples, num);

else
    cube = reshape(cube, numOfSamples, channels);
end

cube = cube';

% Determine the maximum projection.
maxProjections = cast(zeros(size(cube, 2), 1), class(cube));

% Generate random unit vectors.
prevS = rng('default');
unitVec = randn(numOfUnitVec, size(cube,1));

% Clean up the random state.
cleanup = onCleanup(@() rng(prevS));
for numVect = 1:numOfUnitVec
    [~, idx] =  max(abs(unitVec(numVect,:) *cube));
    maxProjections(idx) = maxProjections(idx) + 1;
end

% Find the maximum projected votes to decide the endmember signatures.
[~, idx] = sort(maxProjections, 'descend');

% Extract endmembers form the volume.
endmemSig = originalCube(idx(1:num),:)';
end

function mustBeA(hypercube,classIn)

if isobject(hypercube)
    validateattributes(hypercube,{classIn},{'nonempty'}, 'ppi','hypercube');
    cube = hypercube.DataCube;
else
    cube = hypercube;
end
validateattributes(cube, ...
    {'numeric'}, {'nonsparse', 'nonempty', 'real','nonnan', ...
    'finite', 'ndims', 3}, mfilename, 'hypercube', 1);

end

% Validate the numComponents are less than number of bands.
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

总的来说，PPI算法是一种基于几何的无监督端元提取方法，通过将降维数据投影到随机单位向量上并找到具有最大投影值的像素来识别端元像素。它通常用于高光谱图像的端元提取和非线性解混分析。