一种高光谱图像端元自动提取方法

摘要

本发明涉及一种高光谱图像端元自动提取方法，利用高光谱数据具有近似于低维流形数据的特性，采用正交投影的原理，通过将提取的端元扩充至正交投影算子中，分析投影后各个像元向量的p范数值，进行端元数目自动估计，从而提取端元。本发明自动提取端元，在一定程度上提高混合像元分解的自动化程度，减少人工介入，该方法有较强的鲁棒性，在信噪比不高的情况下仍能正确地估计端元数目，提取端元。

说明书

技术领域

本发明涉及高光谱遥感应用领域，尤其涉及一种高光谱图像端元自动提取方法。

背景技术

高光谱遥感是遥感技术发展的重要趋势，其光谱通道数多，通常有数十甚至数百个，这 种特性被成功地应用在地质勘探、农林业调查和环境监测等领域的研究，取得了引人瞩目的 成果。高光谱具有图谱合一的特点，能够得到各个像元的光谱曲线和图像信息。随着高光谱 遥感技术的迅速发展，与之相辅的高光谱图像端元自动提取技术也得到了人们越来越多的重 视。

端元数目的正确估计是端元提取的先决条件，端元提取的结果在一定程度上直接决定了 混合像元分解的精度，目前，学者们提出了很多有效的端元提取算法，但大多数算法主要是 根据判读经验来估计端元数目，这往往会因为判读人员的个人经验或者光谱数据中存在较高 噪声而降低端元数目估计的精度。因此，自动且具有较高鲁棒性的端元数目估计是混合像元 分解中端元提取的关键环节。

发明内容

有鉴于此，本发明的主要目的在于提出一种高光谱图像端元自动提取方法，自动提取端 元，在一定程度上提高混合像元分解的自动化程度，减少人工介入，该方法在低信噪比情况 下仍能正确地估计端元数目，并且可以正确地提取端元。

为达到上述目的，本发明的技术方案是这样实现的：

一种高光谱图像端元自动提取方法，包括如下步骤，如图1所示：

步骤A、获取高光谱图像数据；

获取高光谱图像数据X^M×N，其中M是光谱的波段数目，N是所有像元的数目。

步骤B、高光谱图像数据预处理；

为了简化计算以及降低噪声对结果的影响，要将获取的高光谱数据X^M×N进行降维。

步骤B1、求取高光谱数据X^M×N的平均光谱d：

$d=\frac{1}{N}\underset{n=1}{\overset{N}{\Sigma }}x\left[n\right]---\left(1\right)$

其中x[n]表示高光谱数据X^M×N中的一条光谱,n取1，2，…N；

步骤B2、将高光谱数据X^M×N进行均值零化处理，处理后的光谱数据为U：

U＝[x[1]-d,x[2]-d,…,x[N]-d]∈R^M×N(2)

步骤B3、求取仿射变换算子C：

C＝[q₁(UU^T),q₂(UU^T),…,q_pmax-1(UU^T)](3)

其中q_i(UU^T)表示求取UU^T第i个特征值所对应的特征向量，pmax最大端元数；

步骤B4、将均值零化处理后的数据仿射变换，求取降维后的像元向量

$\tilde{x}\left[n\right]=C^T\left(x\right[n]-d)\in R^{pmax-1}---\left(4\right)$

步骤B5、将降维后的像元向量添加一维，获取增广降维光谱数据

$\bar{x}\left[n\right]=\left(\begin{array}{c}\tilde{x}\left[n\right]\\ 1\end{array}\right)---\left(5\right)$

步骤C、求取第一个端元；

根据公式(6)，求取增广降维后的一个端元${\bar{m}}_1:$

$\begin{array}{cc}{\bar{m}}_1=\tilde{x}\left[l_1\right]& l_1\in \arg \underset{n=1,\mathrm{...},N}{\max }\left\{\left|\right|\bar{x}\left[n\right]\left|\right|\right\}\end{array}---\left(6\right)$

步骤D、将求取的降维后端元组成矩阵投影算子Q：

$Q=\left[{\bar{m}}_1\right]\in R^{pmax\times 1}---\left(7\right)$

pmax为最大端元数；

步骤E、根据矩阵Q不断迭代提取端元,并将提取的端元添加到Q中，最后获得端元光谱； 步骤E1、求取Q的正交投影算子

其中I_pmax是单位矩阵，维数为pmax，表示矩阵Q的伪逆。

步骤E2、求取不同于的端元

$\begin{array}{cc}{\bar{m}}_j=\bar{x}\left[l_j\right]& l_j\in \arg \underset{n=1,...,N}{max}\left\{\right|\left|P_Q^{\perp }\bar{x}\right[n\left]\right|{|}_p\}\end{array}---\left(9\right)$

步骤E3、将新求取的端元添加到矩阵Q中，组成新的矩阵Q：

$Q=\left[\begin{array}{cc}Q& {\bar{m}}_j\end{array}\right]=[{\bar{m}}_1,...,{\bar{m}}_k]\in R^{pmax\times k}$

步骤E4、求取的中位数计算 $D_j=\underset{n=1,\mathrm{...},N}{\max }\left\{\right|\left|P_Q^{\perp }\bar{x}\right[n\left]\right|{|}_p\}-\underset{n=1,\mathrm{...},N}{median}\{\left|\right|P_Q^{\perp }\bar{x}\left[n\right]\left|{|}_p\right\},$p可以取任意正整数，通常p取2；

步骤E5、k未达到最大端元数pmax时，重复第E1步，达到最大端元数pmax时，计算 判定值$T_i=\frac{D_i}{\underset{i=1,\mathrm{...},P\max }{\min }\left\{D_i\right\}};$

步骤E6、设定阈值T_threshold，估计端元数目，比较判定值T与阈值的关系，当T_i≤T_threshold时，i便是要进行估计的端元数目；

步骤E7、输出端元光谱集合Q。

本发明与现有技术相比的优点在于：

本发明利用高光谱数据具有近似于低维流形数据的特性，采用正交投影的原理，将混合 像元投影至正交投影算子中，通过寻找投影向量p范数最大的点即单形体的顶点来筛选端元， 将筛选出的端元再次添加到投影算子中，不断迭代。本方法无需设定初始端元向量，避免了 常规方法中随机设定端元初始向量带来的端元提取结果的不确定性，提高了端元提取的精度 和自动化程度。

在端元筛选的过程中，本发明提出了自动估计端元数目的判别函数，通过分析投影后各 个像元向量的p范数值，将每一次投影中的p范数最大值与p范数中位数作差，并将其与所 有次投影中p范数最大值与p范数中位数之差的最小项作商，通过多次投影的商值与判定值 阈值进行比较，从而实现端元数目自动估计。该方法解决了目前大多数端元提取算法首先根 据经验人为确定端元数目导致的端元提取精度降低的问题。该方法简洁，容易实现，精度高， 并且具有较强的鲁棒性。

附图说明

图1为本发明方法的流程图；

图2为实验采用的USGS库中的5种端元光谱图；

图3为实验采用高光谱图像仿真图；

图4为提取的端元顺序图；

图5为判定值T与投影次数的关系。

具体实施方式

为使本发明的目的、技术方案和优点更加清楚明白，以下举实施例并参照附图，对本发 明进一步详细说明。

为了明确地探究端元数据估计方法的优劣性，将使用仿真图像作为实验数据。选取USGS (美国地质勘探局)光谱库中的5种材料光谱(记作P₁、P₂、P₃、P₄、P₅)按照一定比例构 成仿真图像，波长范围从350nm-1000nm。仿真图像包含25个方形区域，所在位置如图3所 示，图像中其它区域像元称为背景像元，大小为5种物质光谱的均值(简写为B)，相应的 背景光谱特征在图2中画出。这25个方形区域静仿真按如下的思路：1-5行第一列为个4×4 的纯像元，1-5行第二列为2×2的纯像元，1-5行第三列为2×2的混合像元，1-5行第四列为 1×1的子像元，1-5行第五列为1×1的子像元。引入第三列的混合像元和第四、第五列的子 像元的目的是为了研究并分析在五种端元不同混合程度时算法的性能。具体的混合信息如表 1所示。在表1中，P代表区域中所有像元，P_ij代表区域中第i行、第j列。在仿真图的第一 列和第二列中总共有100个纯像元，第三列中有20个混合像元，第四列和第五列中有10个 子像元。此外，在仿真图像中分别添加30dB的高斯白噪声。

表1仿真合成图像的像元组成

参照图1，该实施例高光谱图像端元自动提取实现过程如下：

步骤A、获取高光谱图像数据；

获取高光谱图像数据X^651×130，其中651是光谱的波段数目，130是所有像元的数目。

步骤B、高光谱图像数据预处理；

为了简化计算以及降低噪声对结果的影响，要将获取的高光谱数据X^651×130进行降维。

步骤B1、求取高光谱数据的平均光谱d：

$d=\frac{1}{130}\underset{n=1}{\overset{130}{\Sigma }}x\left[n\right]---\left(1\right)$

步骤B2、将高光谱数据X^651×130进行均值零化处理，处理后的光谱数据为U：

U＝[x[1]-d,x[2]-d,…,x[130]-d]∈R^651×130(2)

步骤B3、设定最大端元数目pmax＝20，求取仿射变换算子C：

C＝[q₁(UU^T),q₂(UU^T),…,q₁₉(UU^T)](3)

步骤B4、将均值零化处理后的数据仿射变换，求取降维后的像元向量

$\tilde{x}\left[n\right]=C^T\left(x\right[n]-d)\in R^{19}---\left(4\right)$

步骤B5、获取增广降维光谱数据

$\bar{x}\left[n\right]=\left(\begin{array}{c}\tilde{x}\left[n\right]\\ 1\end{array}\right)---\left(5\right)$

步骤C、求取第一个端元；

求取增广降维后的一个端元m₁：

$\begin{array}{cc}{\bar{m}}_1=\tilde{x}\left[l_1\right]& l_1\in \arg \underset{n=1,\mathrm{...},N}{\max }\left\{\left|\right|\bar{x}\left[n\right]\left|\right|\right\}\end{array}---\left(6\right)$

步骤D、将求取的增广降维端元组成矩阵Q：

$Q=\left[{\bar{m}}_1\right]\in R^{19\times k}---\left(7\right)$

步骤E、不断提取端元，将提取的端元添加到Q；

步骤E1、求取Q的正交投影算子

其中表示矩阵Q的伪逆。

步骤E2、求取不包含在Q中端元

$\begin{array}{cc}{\bar{m}}_j=\bar{x}\left[l_j\right]& l_j\in \arg \underset{n=1,...,130}{max}\left\{\right|\left|P_Q^{\perp }\bar{x}\right[n\left]\right|{|}_2\}\end{array}---\left(9\right)$

步骤E3、将新求取的端元添加到矩阵Q中，组成新的矩阵Q：

$Q=\left[\begin{array}{cc}Q& {\bar{m}}_j\end{array}\right]=[{\bar{m}}_1,...,{\bar{m}}_k]\in R^{pmax\times k}$

步骤E4、求取的中位数计算 $D_j=\underset{n=1,\mathrm{...},N}{\max }\left\{\right|\left|P_Q^{\perp }\bar{x}\right[n\left]\right|{|}_p\}-\underset{n=1,\mathrm{...}N}{median}\{\left|\right|P_Q^{\perp }\bar{x}\left[n\right]\left|{|}_p\right\},$p取值为2；

步骤E5、未达到最大迭代次数20时，重复第E1步，达到最大迭代次数20时，计算判 定值$T_i=\frac{D_i}{\underset{i=1,\mathrm{...},20}{\min }\left\{D_i\right\}};$

步骤E6、设定阈值T_threshold＝2估计端元数目；

T值与投影次数的关系如图5所示，T₅≤T_threshold时，估计端元数目为5。

步骤E7、输出端元光谱集合，其在仿真图像中的位置如图4所示。

提供以上实施例仅仅是为了描述本发明的目的，而并非要限制本发明的范围。本发明的 范围由所附权利要求限定。不脱离本发明的精神和原理而做出的各种等同替换和修改，均应 涵盖在本发明的范围之内。