最大距离(MaxD)算法提取高光谱图像端元

最大距离(MaxD)算法是一种无监督的端元提取算法,用于从高光谱图像数据中识别出纯净的端元光谱。本文将详细介绍MaxD算法的数学原理和MATLAB实现。

算法原理

MaxD算法基于以下几个关键思想:

1. 最大距离准则: 端元光谱对应于数据样本在特征空间中最远的点。这是因为端元光谱代表了纯净的地物光谱,而混合光谱位于端元光谱之间。

2. 逐步提取: 算法按照端元光谱与其他光谱的距离从大到小的顺序,逐步提取出所有端元光谱。

3. 向量投影: 在提取新的端元光谱时,利用已知端元光谱的线性组合构造一个向量空间,新的端元光谱在该空间的complementary空间中具有最大投影长度。

设有 \(N\) 个波段的高光谱数据 \(\mathbf{X} = \{\mathbf{x}_1, \mathbf{x}_2, \ldots, \mathbf{x}_n\}\),其中 \(\mathbf{x}_i \in \mathbb{R}^N\) 为第 \(i\) 个样本光谱,我们需要从中找出 \(p\) 个端元光谱 \(\mathbf{E} = \{\mathbf{e}_1, \mathbf{e}_2, \ldots, \mathbf{e}_p\}\)。算法步骤如下:

1. 初始化: 计算所有样本光谱的范数 \(\|\mathbf{x}_i\|_2\),选取范数最大的光谱作为第一个端元光谱 \(\mathbf{e}_1\):

   \[\mathbf{e}_1 = \mathbf{x}_{i_1}, \quad i_1 = \arg\max_i \|\mathbf{x}_i\|_2\]

2. 提取第二个端元光谱: 对于剩余的每个光谱 \(\mathbf{x}_j (j \neq i_1)\),计算其与第一个端元光谱 \(\mathbf{e}_1\) 的距离 \(d_j = \|\mathbf{x}_j - \mathbf{e}_1\|_2\),选取距离最大的光谱作为第二个端元光谱 \(\mathbf{e}_2\):

   \[\mathbf{e}_2 = \mathbf{x}_{i_2}, \quad i_2 = \arg\max_j d_j\]

3. 提取第三个端元光谱: 对于剩余的每个光谱 \(\mathbf{x}_k (k \neq i_1, i_2)\),计算其在由 \(\mathbf{e}_1\) 和 \(\mathbf{e}_2\) 张成的平面的补空间中的投影长度 \(d_k^{(3)}\),选取投影长度最大的光谱作为第三个端元光谱 \(\mathbf{e}_3\):

   \[\begin{aligned} \mathbf{v}_{12} &= \mathbf{e}_2 - \mathbf{e}_1 \\ \mathbf{v}_{k} &= \mathbf{x}_k - \mathbf{e}_2 \\ \theta_k &= \arccos\left(\frac{\mathbf{v}_k^T\mathbf{v}_{12}}{\|\mathbf{v}_k\|_2\|\mathbf{v}_{12}\|_2}\right) \\ d_k^{(3)} &= \|\mathbf{v}_k\|_2 \sin\theta_k \\ \mathbf{e}_3 &= \mathbf{x}_{i_3}, \quad i_3 = \arg\max_k d_k^{(3)} \end{aligned}\]

4. 提取第四个端元光谱: 对于剩余的每个光谱 \(\mathbf{x}_l (l \neq i_1, i_2, i_3)\),计算其在由 \(\mathbf{e}_1\)、\(\mathbf{e}_2\) 和 \(\mathbf{e}_3\) 张成的空间的complementary空间中的投影长度 \(d_l^{(4)}\),选取投影长度最大的光谱作为第四个端元光谱 \(\mathbf{e}_4\):

   \[\begin{aligned} \mathbf{E}_3 &= [\mathbf{e}_1, \mathbf{e}_2, \mathbf{e}_3] \\ r &= \mathrm{rank}(\mathbf{E}_3) \\ \mathbf{G} &= \mathrm{null}(\mathbf{E}_3, r) \\ \mathbf{n} &= \mathbf{G}(:,1) \\ d_l^{(4)} &= \frac{|\mathbf{v}_l^T\mathbf{n}|}{\|\mathbf{n}\|_2} \\ \mathbf{e}_4 &= \mathbf{x}_{i_4}, \quad i_4 = \arg\max_l d_l^{(4)} \end{aligned}\]

5. 提取剩余端元光谱: 对于第 \(m\) 个端元光谱 \((m > 4)\),重复上一步的过程,但使用前 \(m-1\) 个端元光谱计算complementary空间和投影长度。算法重复进行,直到没有新的端元光谱可提取为止。

通过上述步骤,MaxD算法能够从高光谱数据中逐步提取出所有的端元光谱。需要注意的是,该算法对初始光谱的顺序不敏感,但需要预先设定端元个数 \(p\)。

MATLAB实现

以下是MaxD算法在MATLAB中的实现代码:

function [Endmember]=MaxD(X,num_End)
   % 将输入数据X赋值给mixed
   mixed = X;
   % 计算每个样本光谱的范数
   D_O = sum(mixed.^2,2).^(1/2);
   % 找出范数最大的样本光谱的索引，作为第一个端元
   [~,E_lab1]=max(D_O);
   % 将第一个端元存储在E的第一行
   E(1,:)=E_lab1(1,1);
   %% The second endmember
   % 计算其余样本光谱与第一个端元的距离
   D_1p = sum((mixed(E_lab1(1,1),:)-mixed).^2,2).^(1/2);
   % 找出距离最大的样本光谱的索引，作为第二个端元
   [~,E_lab2] = max(D_1p);
   % 将第二个端元存储在E的第二行
   E(2,:)=E_lab2;
   %% The third endmember
   % 计算第二个端元与第一个端元的向量差
   v_12=(mixed(E_lab2,:)-mixed(E_lab1(1,1),:));
   % 计算其余样本光谱与第二个端元的向量差
   v_2p = mixed-mixed(E_lab2,:);
   % 计算其余样本光谱在由v_12张成的平面的补空间中的投影长度的余弦值
   cos = v_2p*v_12'./((sum(v_2p.^2,2).^(1/2))*(sum(v_12.^2,2).^(1/2)));
   % 计算v_2p的长度
   norm_v_p2 = sum(v_2p.^2,2).^(1/2);
   % 计算其余样本光谱在由v_12张成的平面的补空间中的投影角度
   angle = acos((cos));
   % 计算其余样本光谱在由v_12张成的平面的补空间中的投影长度
   D_p_12 = norm_v_p2.*sin(angle(:));
   % 找出投影长度最大的样本光谱的索引，作为第三个端元
   [~,E_lab3]=max(D_p_12);
   % 将第三个端元存储在E的第三行
   E(3,:)=E_lab3;
   %% The forth endmember
   % 提取前三个端元构成的矩阵
   E_3 = mixed(E,:);
   % 计算E_3的秩
   r =rank(E_3);
   % 计算E_3的null空间的基
   Gen_sol = null(E_3,r);
   % 取null空间的第一个基向量
   n = Gen_sol(:,1);
   % 计算其余样本光谱在由n张成的空间中的投影长度
   [~,E_lab4] = max(abs(v_2p*n)/norm(n));
   % 将投影长度最大的样本光谱索引作为第四个端元
   E(4,:) = E_lab4;
   %% The next endmember
   % 对于剩余的端元
   for i = 5:num_End
       % 提取前i-1个端元构成的矩阵
       A5 = mixed(E,:);
       % 计算A5的秩
       r5 = rank(A5);
       % 计算A5的null空间的基
       Gen_sol5 = null(A5,r5);
       % 取null空间的第一个基向量
       n5 = Gen_sol5(:,1);
       % 取第一个端元对应的样本光谱
       x=mixed(E(1,1),:);
       % 计算x在由n5张成的空间中的投影长度的常数项
       b5 = abs((x*n5)/norm(n5));
       % 计算其余样本光谱在由n5张成的空间中的投影长度
       DL5=mixed*n5+b5;
       % 计算投影长度与n5长度的比值
       dis5 = (sum(abs(DL5).^2,2).^(1/2))/norm(n5);
       % 将dis5从大到小排序，获取排序后的值和对应索引
       [rol5,cow5]=sort(dis5,'descend');
       % 如果最大投影长度比值大于阈值0.00001
       if rol5(1,1)>0.00001
           % 则将对应的样本光谱索引作为第i个端元
           E(i,:)=cow5(1,1);
           a = rol5(1,1);
       else
           % 否则退出循环
           break
       end
   end
   % 根据端元索引从mixed中提取端元光谱
   Endmember = mixed(:,E);

end

该函数的输入参数包括:

• X: 高光谱数据矩阵,每一行对应一个光谱样本
• num_End: 需要提取的端元个数

函数的输出是一个矩阵 Endmember,每一列对应一个提取出的端元光谱。

代码主要按照前面介绍的算法步骤实现了端元提取过程。值得注意的是,在提取第四个及之后的端元时,代码使用了矩阵的null空间计算,以得到complementary空间的基向量。

通过这个MATLAB函数,我们可以方便地将MaxD算法应用于实际的高光谱数据分析中,获得样本的端元成分。

总的来说,MaxD算法是一种重要的无监督端元提取方法,具有数学理论简单、易于实现的特点。它在高光谱遥感、矿物分析等领域有着广泛的应用。