基于N-FINDR端元提取算法的线性光谱分解

混合像元光谱分解模型用于分离和提取各像元组分，通常有线性光谱混合模型和非线性光谱混合模型两种类型。无论哪种类型的光谱分解模型，都是把像元的反射率表示为端元（endmember）组分的光谱特征和它们丰度的函数。因此，端元提取是光谱分解中重要一步，端元提取的质量直接影响光谱分解结果精度。

常见的端元提取算法有以下几种：

1. 人工选取端元。通常有两条途径，一是实地测量或直接从光谱库中获取，二是从影像分析中获得。
2.  像元纯净指数法：基于凸面几何理论，即认为在理想情况下高光谱影像的所有N维空间数据点由以图像中所有存在的纯净像元为顶点的单形体所包围。该算法将所有N维数据点反复投影在一个随机向量上，端元更有可能出现在向量的两端。因此，反复多次投影，统计各点出现在向量的两端的次数（即像元纯净指数PPI），PPI越大说明该像元为端元的可能性更高。
3. N-FINDR算法：也是基于凸面几何理论。与像元纯净指数法不同的是，N-FINDR算法的基本思想是将N维空间中体积最大的单形体的顶点作为端元。

---

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

第一步：利用PCA、MNF等手段对原始影像进行降维，保留降维过后的前维，并选取个像元作为初始端元，构成的增广矩阵：

其中为端元的光谱向量，为端元数。

第二步：按公式计算这个端元构成的单形体体积：

第三步：用不同的像元以此替换不同端元，如果单形体体积增大，则替换该端元。

第四步：获取最大单形体对应的端元以及该组端元对应的光谱向量，作为输出结果进行后续操作。

从上述步骤中不难看出，算法的关键在于获取最大单形体。然而，由于像元数量众多，穷举每组端元并计算单形体体积效率太低；而其他方法端元提取算法的结果可能不是不完备最大体积单形体。基于此有许多改进的N-FINDR算法，如改善端元的选取方式：串型N-FINDR算法(SQ N-FINDR)及并型N-FINDR算法(SC N-FINDR)等。

---

代码实现

MATLAB实现遍历循环获取非完备最大单形体对应端元及光谱向量：

clear,close all

% MFN file
MFN_fileName='dataMNF';
[MFN_image,MFN_p,MFN_t]=freadenvi(MFN_fileName);
MFN_image=MFN_image';

n=12;%端元个数
% 选取前12个像素作为初始端元集合
select_index=1:n;
select_pixel_bands=MFN_image(1:n-1,select_index);
E=cat(1,ones(1,n),select_pixel_bands);
VE=abs(det(E))/factorial(n-1);

% 固定端元中的n-1个像元，使用其余像元依次替代第i个像元
% 如果单形体体积增加，则替代原来的像素点
times=0;
while 1
    VE_n=VE;
    times=times+1;
    for i=1:n
        for j=1:400*350
            tmp_index=select_index;
            if ismember(j,select_index)==0
                tmp_index(1,i)=j;
                tmp_bands=MFN_image(1:n-1,tmp_index);
                tmp_E=cat(1,ones(1,n),tmp_bands);
                tmp_VE=abs(det(tmp_E))/factorial(n-1);
                if tmp_VE>VE
                    select_index=tmp_index;
                    VE=tmp_VE;
                end
            end
        end
    end

    if abs(VE_n-VE)<1 || times>10
        break;
    end
end

VE_max_pixel_bands=MFN_image(1:n-1,select_index);
B=[];
for i=1:400*350
    pixel_bands=MFN_image(1:n-1,i);
    B=cat(1,B,(VE_max_pixel_bands\pixel_bands)');
end

其中freadenvi函数用于获取HDR图像信息。代码思路是将固定个端元，按所有像元排列顺序不断替换其中一个端元（若该像元已在端元中则跳过）。由于一次遍历获取单形体的可能不是体积最大单形体，迭代计算直到单形体体积不再变化。

获取到端元光谱向量后，逐像元计算端元丰度值。按照线性模型，其中为端元光谱向量，为某像元光谱向量，则该像元端元丰度。

---

实验结果

MNF图像（第一分量）：

端元丰度图（仅展示部分）：

可以看出，不同端元的丰度在图像上可较好地显示出来。然而，由于未增加丰度非负值及丰都总和为1的限制条件，实验结果可能不符合实际物理意义。此外，MATLAB中Image Processing Toolbool Hyperspectral Imaging Library中有一个nfindr函数，可直接用于端元提取。

MFN_fileName='dataMNF';
[MFN_image,MFN_p,MFN_t]=freadenvi(MFN_fileName);
MFN_image=reshape(MFN_image,MFN_p);
MFN_image=MFN_image(:,:,1:12);%保留前12个分量
endmembers=nfindr(MFN_image,12);

效果如下：