CN103617424B

CN103617424B - 一种高光谱图像端元个数自动估计的端元提取方法

Info

Publication number: CN103617424B
Application number: CN201310601834.3A
Authority: CN
Inventors: 许毅平; 田岩; 谭志敏
Original assignee: Huazhong University of Science and Technology
Current assignee: Huazhong University of Science and Technology
Priority date: 2013-11-25
Filing date: 2013-11-25
Publication date: 2016-05-25
Anticipated expiration: 2033-11-25
Also published as: CN103617424A

Abstract

本发明公开了一种高光谱图像端元个数自动估计的端元提取方法，首先使用PPI算法计算高光谱图像中所有像元的PPI值，根据像元的PPI值确定初始候选端元集，并对初始候选端元按PPI值由大到小进行排序；然后依次剔除初始候选端元集中独立性最弱的候选端元，同时获得随端元个数递减的候选端元集的重构精度曲线和端元独立性曲线，通过寻找这两条曲线变化的显著性位置确定端元个数的取值范围；最后利用SGA算法获取不同端元个数下的二次选择候选端元集，根据这些候选端元集所构成的重构精度曲线，确定最终的端元个数，同时获得最终的端元提取结果。本发明不需要事先给定端元的个数，有利于端元提取的自动化处理，同时提取的端元准确性高。

Description

一种高光谱图像端元个数自动估计的端元提取方法

技术领域

本发明属于高光谱遥感应用领域，具体涉及一种高光谱图像端元个数自动估计的端元提取方法。

背景技术

在高光谱图像处理中，端元提取是高光谱图像地物分类和识别的关键技术之一。依靠所提取到端元进行光谱匹配，可以准确有效地判断地物类别。目前常用的端元提取方法在进行端元提取时，需要事先给出端元的个数，而端元个数的估计往往需要人工根据经验对场景进行分析来确定，影响了高光谱图像自动化处理的实现。

发明内容

针对现有技术的不足，本发明的目的在于提出一种高光谱图像端元个数自动估计的端元提取方法，减少人工介入，自动提取端元。

本发明解决其技术问题所采用的技术方案是：

一种高光谱图像端元个数自动估计的端元提取方法，包括以下步骤：

S1、使用PPI算法计算每个像元的PPI值，将PPI值大于0的像元作为初始候选端元，并用矩阵E_L×N表示，其中L为数据维数，N为候选端元的个数；

S2、对初始候选端元按PPI值从大到小进行排序，根据候选端元间的相关性逐个剔除端元，并获得候选端元集的光谱重构精度曲线R和独立性曲线I；

S3、根据光谱重构精度曲线和独立性曲线随端元个数减少的变化率，确定端元个数的取值范围；

S4、在步骤S3确定的端元个数取值范围内，利用SGA算法提取不同端元个数所对应的候选端元集，并计算这些候选端元集的光谱重构精度曲线，选择曲线变化最显著的位置值作为最终的端元个数值，并将该个数值对应应的候选端元集作为最终的端元集。

本发明的有益效果为：

（1）提取端元时，不需要事先给定端元的个数，减少了人工的介入，有利于端元提取的自动化处理。

（2）提取端元时，综合考虑了端元具有像元纯净值高和相关性弱的特点，通过候选端元的迭代剔除，使提取的端元准确性更高，抗干扰能力更强。

附图说明

参照下面的说明，结合附图，可以对本发明有最佳的理解。在附图中，相同的部分可由相同的标号表示。

图1为本发明较佳实施例的端元提取流程图；

图2为真实光谱数据；

图3为五种目标端元与背景1的光谱曲线图；

图4为五种目标端元与背景2的光谱曲线图；

图5为构造的四组高光谱图像测试数据，其中图5(a)采用背景BKG，以TI方式插入小子块；图5(b)采用背景BKG2，以TI方式插入小子块；图5(c)采用背景BKG，以TE方式插入小子块；图5(d)采用背景BKG2，以TE方式插入小子块；

图6为将图5的4个测试图像分别作为输入数据得到的端元提取结果。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及示例性实施例，对本发明进行进一步详细说明。应当理解，此处所描述的示例性实施例仅用以解释本发明，并不用于限定本发明的适用范围。

本发明的基本思想为：结合纯像元指数算法（PixelPureIndex，PPI）和单形体增长算法（SimplexGrowingAlgorithm，SGA），自适应地估计端元个数并提取端元光谱。由于在地物光谱的遥感测量过程中受各种因素的影响，所获取的地物光谱信息包含随机噪声，这些噪声的存在会影响光谱曲线的特性，使光谱曲线之间存在较大的相关性，阻碍了端元个数的估计。因此，在端元个数的估计过程中，采取逐渐逼近的方法，通过逐个排除非端元像元，并分析候选端元集的光谱重构能力和端元间的独立性，根据重构精度曲线和相关性曲线的变化规律，选择曲线拐点位置所对应的端元数作为最终的端元个数，并提取相应的端元。

如图1所示，本发明较佳实施例的一种高光谱图像端元个数自动估计的端元提取方法，包括以下步骤：

S1、使用PPI算法计算每个像元的PPI值，将PPI值大于0的像元作为初始候选端元，并用矩阵E_L×N表示，其中L为数据维数，N为候选端元的个数。

S2、对初始候选端元按PPI值从大到小进行排序，根据候选端元间的相关性逐个剔除端元，并获得候选端元集的光谱重构精度曲线R和独立性曲线I。

其中，步骤S2具体为：

（a）令p＝N，S_p＝E_L×p；

（b）将候选像元根据其PPI值从大到小排序，结果存入S_p中；

（c）计算p个候选端元两两间的光谱角距离，用距离矩阵C_p×p表示，该矩阵为下三角矩阵，对角线上的值设为0；

（d）根据式(1)得到I(p)，根据式(2)～(4)得到R(p)；

I(p)＝c_i,j(1)

α_{p} = {(S_{p}^{T} \cdot S_{p})}^{- 1} \cdot S_{p}^{T} \cdot γ - - - (2)

{\hat{γ}}_{p} = S_{p} \cdot α_{p} - - - (3)

R (p) = \cos < γ, {\hat{γ}}_{p} > = \frac{γ \cdot \hat{γ}}{| γ | | \hat{γ} |} - - - (4)

其中，c_i,j为距离矩阵C_p×p中的最大值，γ为图像数据的均值向量；

（e）删除矩阵C_p×p中的第i行和第j列；

（f）令p＝p-1；如果p＞0，则跳至（d）继续，否则退出S2。

S3、根据光谱重构精度曲线和独立性曲线随端元个数减少的变化率，确定端元个数的取值范围。

具体而言，按以下步骤估计端元个数的取值范围：

（a）根据计算式(5)和(6)分别计算光谱重构精度曲线R和独立性曲线I的二阶向后差分的绝对值，分别记为Diff_R(p)和Diff_I(p)；

Diff_R(p)＝|[R(p)-R(p-1)]-[R(p-1)-R(p-2)]|(5)

Diff_I(p)＝|[I(p)-I(p-1)]-[I(p-1)-I(p-2)]|(6)

（b）寻找Diff_R和Diff_I的前两个极小值，分别用n_R1、n_R2、n_I1和n_I2表示，则端元个数的取值范围为[n₁,n₂]，其中n₁＝max(n_R1,n_I1)、n₂＝max(n_R2,n_I2)。

S4、在步骤S3确定的端元个数取值范围内，利用SGA算法提取不同端元个数所对应的候选端元集，并计算这些候选端元集的光谱重构精度曲线，选择曲线变化最显著的位置值作为最终的端元个数值，并将该值对应应的候选端元集作为最终的端元集。

其中，步骤S4具体为：

（a）取步骤S2中p为2×n₂时的候选端元集作为SGA算法的输入数据；

（b）设定端元个数从1变化到2×n₂。每设定一次端元个数，用SGA算法提取相应的端元，形成端元矩阵，记为并用式(7)和(8)计算单形体体积，其中，s_i为端元向量，p表示端元个数。另外根据计算式(2)～(4)得到第二条重构精度曲线R^SGA；

E = {|\begin{matrix} s_{1} \cdot s_{1} & s_{1} \cdot s_{2} & \cdot \cdot \cdot & s_{1} \cdot s_{p} \\ \cdot \cdot \cdot & \cdot \cdot \cdot & \cdot \cdot \cdot & \cdot \cdot \cdot \\ s_{p} \cdot s_{1} & s_{p} \cdot s_{2} & \cdot \cdot \cdot & s_{p} \cdot s_{p} \end{matrix}|}^{\frac{1}{2}} - - - (7)

V (E) = \frac{1}{p!} E - - - (8)

（c）根据计算式(9)，得到曲线Diff_R^SGA。在p的取值范围[n₁,n₂]内，寻找Diff_R^SGA下降幅度最大的点，该点的p值为最终的端元个数估计值，记为Npure。

Diff_R^SGA(p)＝10lg(|[R^SGA(p)-R^SGA(p-1)]-

[R^SGA(p-1)-R^SGA(p-2)]|)(9)

（d）选择端元个数为Npure的端元集作为最终的端元提取结果。

以下利用仿真测试图像数据来验证本方法的有效性。

图2所示为真实的高光谱图像，从图中选取5种不同目标地物的纯像元作为目标端元（图中分别标记为+1、+2、+3、+4和+5），另外选取图中BKG标识区域的像元的平均值作为待构建数据的背景端元，标记为BKG。图3为6个端元的光谱曲线图，从图3中可以看到背景像元和4号目标端元的光谱曲线很相似，两者的光谱角距离为0.9979。因此，可以认为，在构造的图像数据中只存在5种端元。

利用上述6种光谱数据，构建大小为191×64×64的仿真测试图像数据，并在图像中插入5×5个子块，子块的构造方法如下：

第一列子块分别为5个纯端元子块，大小为3×3；

第二列子块分别为5个纯端元子块，大小为2×2；

第三列子块由表1确定，每个像元中都是对应的纯端元占50%，其他任一纯端元占50%的混合像元子块，大小为2×2；

第四和第五列小子块由表2确定。第四列中的像元是对应的纯端元占50%，背景占50%；第五列的像元是对应的纯端元占25%，背景占75%。

表1第三列小子块的像元

表2第四列和第五列小模块的像元

	第四列	第五列
			第一行	0.5A+0.5BKG	0.25A+0.75BKG
第二行	0.5B+0.5BKG	0.25B+0.75BKG
			第三行	0.5C+0.5BKG	0.25C+0.75BKG
第四行	0.5D+0.5BKG	0.25D+0.75BKG
			第五行	0.5E+0.5BKG	0.25E+0.75BKG

设定前两列小子块的目的是为了验证本发明在进行端元提取时，能不能排除处于不同空间位置的相同纯端元对端元提取的干扰。设置后三列中子块的目的是验证本发明会不会漏掉纯端元而将混合像元作为最终的端元输出，即验证本发明能不能提取到数据凸集的顶点。

由于真实的光谱数据中都是存在噪声的，所以在构建的数据中插入信噪比为20:1的高斯白噪声，使得数据更逼近真实光谱数据。构造数据的过程是先构造好背景，然后加入高斯噪声，最后插入含纯端元的小子块。

若称待插入的小子块为目标，在构建数据时，目标的插入方法有两种：

（1）目标植入（TI：TargetImplantation）：该方法是首先在背景中加入高斯噪声（信噪比为20:1），然后将待插入的小子块替换背景中相应的像元。按该方法构造的子块在后续实验中以TI为前缀命名。

（2）目标嵌入（TE：TargetEmbeddedness）：该方法是把待插入的小子块与相应背景中的像元进行叠加而得，按该方法构造的子块在后续实验中以TE为前缀命名。

在构造数据时，除了选择图1所示的BKG背景端元外，另外还选择了与5种目标端元差异较大的光谱作为背景像元（记为BKG2）。BKG2光谱曲线以及5个端元的光谱曲线如图4所示。以BKG2作为背景时，由于BKG2的光谱与其他5个端元差异较大，因此，可以认为由这些端元构建的场景中有6种不同的端元。由此构建四组不同的测试图像，分别命名为TI_1,TI_2,TE_1,TE_2（TI_1表示采用背景BKG，以TI方式插入小子块；TI_2表示采用背景BKG2，以TI方式插入小子块；TE_1表示采用背景BKG，以TE方式插入小子块；TE_2表示采用背景BKG2，以TE方式插入小子块），其中TI_1和TE_1图像中的端元个数为5个，TI_2和TE_2图像中的端元个数为6个。图5为构造的四个高光谱测试图像第32号波段的单波段图像。

将图5所示的4个测试图像分别作为本发明的输入数据，可分别得到不同图像的端元提取结果，如图6所示。

从图6的实验结果看，虽然端元个数自动估计的端元提取算法对端元个数的估计值比真实值偏大，但是四组数据的结果显示，应该被提取到的纯端元都包含在了结果中，并且这些端元都是优先提取到的。即，如果数据中真实端元个数为5个，而自适应端元个数的端元提取算法的估计个数为7，那么，提取端元时，提取的前5个端元为数据中的真实端元，后两个为多提取的。根据该构造数据的实验结果，多提取的端元有可能会与前5个端元有重复。这是因为自适应端元个数的端元提取算法运用到了构成最大体积的顶点即为端元的假设。在这种假设下，算法认为当出现这种重复时得到的体积才能是最大的，所以结果中可能会出现重复现象。另一方面，由于只有在这种构造的数据中，其形成的数据凸集才是完全符合线性假设模型的，而真实的光谱数据中则不会，所以这种重复提取的现象在输入真实高光谱数据时并不会出现。

以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。

Claims

1.一种高光谱图像端元个数自动估计的端元提取方法，包括以下步骤：

S4、在步骤S3确定的端元个数取值范围内，利用SGA算法提取不同端元个数所对应的候选端元集，并计算这些候选端元集的光谱重构精度曲线，选择曲线变化最显著的位置值作为最终的端元个数值，并将该个数值对应的候选端元集作为最终的端元集。

2.根据权利要求1所述的端元提取方法，其中，步骤S2具体为：

(a)令p＝N，S_p＝E_L×p；

(b)将初始候选像元根据其PPI值从大到小排序，结果存入S_p中；

(c)计算p个候选端元两两间的光谱角距离，用距离矩阵C_p×p表示，该矩阵为下三角矩阵，对角线上的值设为0；

(d)根据式(1)得到I(p)，根据式(2)～(4)得到R(p)；

I(p)＝c_i,j(1)

α_{p} = {(S_{p}^{T} \cdot S_{p})}^{- 1} \cdot S_{p}^{T} \cdot γ - - - (2)

{\hat{γ}}_{p} = S_{p} \cdot α_{p} - - - (3)

R (p) = \cos < γ, {\hat{γ}}_{p} > = \frac{γ \cdot \hat{γ}}{| γ | | \hat{γ} |} - - - (4)

其中，p为中间变量，c_i,j为距离矩阵C_p×p中的最大值，γ为图像数据的均值向量；

(e)删除矩阵C_p×p中的第i行和第j列；

(f)令p＝p-1；如果p＞0，则跳至(d)继续，否则退出S2。

3.根据权利要求2所述的端元提取方法，其中，步骤S3具体为：

(a)根据计算式(5)和(6)分别计算光谱重构精度曲线R和独立性曲线I的二阶向后差分的绝对值，分别记为Diff_R(p)和Diff_I(p)；

Diff_R(p)＝|[R(p)-R(p-1)]-[R(p-1)-R(p-2)]|(5)

Diff_I(p)＝|[I(p)-I(p-1)]-[I(p-1)-I(p-2)]|(6)

(b)寻找Diff_R和Diff_I的前两个极小值，分别用n_R1、n_R2、n_I1和n_I2表示，则端元个数的取值范围为[n₁,n₂]，其中n₁＝max(n_R1,n_I1)、n₂＝max(n_R2,n_I2)。

4.根据权利要求3所述的端元提取方法，其中，步骤S4具体为：

(a)取步骤S2中p为2×n₂时的候选端元集作为SGA算法的输入数据；

(b)设定端元个数从1变化到2×n₂；每设定一次端元个数，用SGA算法提取相应的端元，形成端元矩阵，记为并用式(7)和(8)计算单形体体积，其中，s_i为端元向量，p表示端元个数；另外根据计算式(2)～(4)得到第二条重构精度曲线R^SGA；

E = {|\begin{matrix} s_{1} \cdot s_{1} & s_{1} \cdot s_{2} & ... & s_{1} \cdot s_{p} \\ ... & ... & ... & ... \\ s_{p} \cdot s_{1} & s_{p} \cdot s_{2} & ... & s_{p} \cdot s_{p} \end{matrix}|}^{\frac{1}{2}} - - - (7)

V (E) = \frac{1}{p!} E - - - (8)

(c)根据计算式(9)，得到曲线Diff_R^SGA，在p的取值范围[n₁,n₂]内，寻找Diff_R^SGA下降幅度最大的点，该点的p值为最终的端元个数估计值，记为Npure；

Diff_R^SGA(p)＝10lg(|[R^SGA(p)-R^SGA(p-1)]-[R^SGA(p-1)-R^SGA(p-2)]|)(9)

(d)选择端元个数为Npure的端元集作为最终的端元提取结果。