CN106778680A

CN106778680A - 一种基于关键波段提取的高光谱图像波段选择方法及装置

Info

Publication number: CN106778680A
Application number: CN201710009535.9A
Authority: CN
Inventors: 黄珍; 赵辽英; 张文强; 厉小润
Original assignee: Hangzhou Dianzi University
Current assignee: Hangzhou Dianzi University
Priority date: 2017-01-06
Filing date: 2017-01-06
Publication date: 2017-05-31
Anticipated expiration: 2037-01-06
Also published as: CN106778680B

Abstract

本发明公开了一种基于关键波段提取的高光谱图像波段选择方法及装置。本发具体实现如下：对高光谱图像确定端元个数并提取端元光谱；对每个端元光谱利用基于三点矢量夹角和变化幅度分析的方法提取关键点子集，合并所有关键点子集构造侯选波段子集；根据波段间相似性具有块状分布的特点，构造局部信息散度矩阵的可视图，确定最佳波段数范围；最后，通过信息量度量法和最优子集选择准则法两类方法中的一种方法优中选优，确定最优波段子集。本发明中的端元光谱的关键波段是不同地物间最具区分性的特征波段，并从关键波段子集中选取最优波段，能缩短后续波段选择的时间，因此，本发明提出的波段选择方法能够提高波段选择的效率。

Description

一种基于关键波段提取的高光谱图像波段选择方法及装置

技术领域

本发明涉及图像处理领域，尤其涉及一种基于关键波段提取的高光谱图像波段选择方法及装置。

背景技术

高光谱遥感是近年发展起来的一种新型遥感探测技术，具有广泛应用前景。高光谱遥感图像一般由上百个波段组成，包含了丰富的空间、辐射和光谱信息。但大量的波段增加了高光谱图像异常检测的时间，波段间的相关性降低了检测精度。因此，有效利用高光谱数据的前提，是选择合适的特征对高光谱数据进行降维，现有两种方法实现降维：一是特征提取，二是波段选择。利用特征提取的方法会导致提取的特征失去原有的物理意义，算法复杂，且不利于图像解释；波段选择的方法比特征提取的方法更简单，更关键的是选择后的波段仍然保留了原始波段的物理含义以及地物的光谱特性。

国内外现有的波段选择算法大体上可以分为4类：(1)基于信息度量的波段选择方法；(2)基于最优子集选择准则方法(3)聚类后分组选择的方法；(4)基于波段重构误差最小的方法。这些方法直接从原始图像的所有波段中选取最优波段子集，波段选择的速度和精度都有进一步提升空间。快速的关键波段子集预提取是提升现有波段选择方法效率的有效手段。但现有的关键波段提取方法存在两个主要问题：(1)基于图像聚类的类平均光谱提取关键波段，高光谱图像混合像元的存在使得基于聚类的硬分类结果不能准确反映地物真实分布；(2)现有的关键波段提取方法只考虑所有的波峰波谷点，结果导致漏提一些明显的关键点或提取出较多不明显的关键点，起不到提升波段选择效率的作用。

发明内容

针对现有技术中的不足，本发明的目的是提供一种基于关键波段提取的高光谱图像波段选择方法及装置，提高高光谱图像关键波段提取效果，有效确定最佳波段数，提升现有多种波段选择方法的效率，解决因波段间的相关性导致高光谱图像异常检测精度低的问题。

为实现上述目的，本发明提供了如下技术方案：

一种基于关键波段提取的高光谱图像波段选择方法，具体实现步骤如下：

步骤1、确定高光谱图像典型地物种类，即端元个数p；

步骤2、对高光谱图像提取p个端元光谱；

步骤3、对每个端元光谱利用基于三点矢量夹角分析的方法提取关键点子集，合并所有关键点子集构造侯选波段子集；

步骤4、根据波段间相似性具有块状分布的特点，构造波段间局部信息散度矩阵的可视图，确定最佳波段数范围；

步骤5、根据所述最佳波段数范围设定待选择的波段数，采用信息量度量法中的任意一种方法或最优子集选择准则法中的任意一种方法，确定最优波段子集。

所述信息量度量的方法是指采用不同信息量度量指标的一大类方法，信息量度量指标包括波段图像的信噪比、偏度、峰度和波段间的K-L散度、互信息等。

所述最优子集选择准则法是指采用不同最优子集选择准则的一大类方法，最优子集选择准则包括丰度估计协方差最小、协方差矩阵行列式最大和噪声最小且协方差矩阵行列式最大等。

进一步地，步骤3所述的基于三点矢量夹角分析的方法具体包括：

3-1.获取第i条端元光谱曲线上第q-1、q、q+1三点的光谱值，即曲线纵坐标值z_q-1、z_q、z_q+1，其中q＝2,3,...,L-1；其中L为波段数；

3-2.计算q-1、q、q+1三点形成的夹角θ，其中0<θ≤π；

3-3.根据判断准则判断第q点是否为关键点，若是，添加到当前第i个关键波段集合KeyBand_i中。

进一步地，所述的步骤3-2的夹角θ的计算方法为：

进一步地，所述的步骤3-3的判断准则如下：

当所述q-1、q、q+1三点的纵坐标值z_q-1、z_q、z_q+1满足{(((z_q-1<z_q)∩(z_q+1<z_q))∪((z_q-1>z_q)∩(z_q+1>z_q)))∩(θ>α)}∪(θ>β)，并且abs(z_q+1-z_q)+abs(z_q+1-z_q)>τ，则q点为关键点，其中，夹角θ由所述q-1、q、q+1三点形成，α、β为决定提取出的关键点个数的预设参数，且满足0<α<β≤π，τ为反映变化幅度的参数。

进一步地，步骤4所述的最佳波段数范围的确定方法具体如下：

4-1.根据波段间相似性具有块状分布的特点，计算各波段与左右相邻三分之一波段数范围内的波段间的相似性度量值；

4-2.根据相似性度量值，描绘所述局部信息散度矩阵的可视图，根据可视图对角线为中心的大小方块数确定最佳波段数范围。

进一步地，步骤4-1所述的波段间相似性度量值的度量指标选用K-L散度；N个像元的高光谱图像第i个波段数据表示为B_i＝[b_i1,b_i2,…b_iN]^T，将波段B_i归一化后得到计算第j个波段图像相对于第i个波段图像的K-L散度第j个波段和第i个波段间的相似性度量值为

本发明还提供了一种基于关键波段提取的高光谱图像波段选择装置，所述装置包括端元数确定模块、端元提取模块、关键波段子集构造模块、相似性度量矩阵可视图绘制模块以及最优波段子集确定模块；

所述的端元数确定模块，用于确定高光谱图像典型地物种类，即端元个数p；

所述的端元提取模块，用于对高光谱图像提取p个端元光谱；

所述的关键波段子集构造模块，用于对每个端元光谱利用基于三点矢量夹角分析的方法提取关键点子集，合并所有关键点子集构造侯选波段子集；

所述的相似性度量矩阵可视图绘制模块，用于根据波段间相似性具有块状分布的特点，构造波段间局部信息散度矩阵的可视图，确定最佳波段数范围；

所述的最优波段子集确定模块，根据所述最佳波段数范围设定待选择的波段数，采用信息量度量法中的任意一种方法或最优子集选择准则法中的任意一种方法，确定最优波段子集。

进一步地，所述高光谱图像波段选择装置还包括应用模块，所述应用模块使用波段选择结果来进行高光谱图像目标检测或分类。

进一步地，所述关键波段子集构造模块构造侯选波段子集的具体步骤包括：

(1)获取第i条端元光谱曲线上第q-1、q、q+1三点的光谱值，即曲线纵坐标值z_q-1、z_q、z_q+1，其中q＝2,3,...,L-1；其中L为波段数；

(2)计算q-1、q、q+1三点形成的夹角θ，其中0<θ≤π；计算方法为：

(3)根据判断准则判断第q点是否为关键点，若是，添加到当前第i个关键波段集合KeyBand_i中；所述的判断准则如下：

(4)利用每一条光谱曲线的关键波段集合KeyBand₁,KeyBand₂,...,KeyBand_p构造候选波段集合KeySet，所述候选波段集合满足KeySet＝KeyBand₁∪KeyBand₂∪KeyBand₃∪...∪KeyBand_p。

本发明的有益效果在于：

(1)针对现有技术中直接对原始图像的所有波段中选取最优波段子集，导致波段选择效率不高的问题，本发明采用先对高光谱图像的原始数据进行关键波段提取，在此基础上再应用信息度量方法或最优子集选择准则方法，优中选优，从而提高了波段选择效率；

(2)本发明提出基于三点矢量夹角分析的关键波段提取新方法，能够尽可能地提取所有明显的关键点，同时避免提取出噪声引起的关键点或不明显的关键点，得到的结果更好。

(3)本发明提出对端元光谱曲线提取关键波段子集，减少了高光谱图像混合像元的普遍存在对关键波段提取效果的影响。

(4)本发明提出基于波段相似性度量矩阵可视图的最佳波段数范围确定新方法，并且只需计算各波段与左右相邻三分之一波段数范围内的波段间的相似性度量值，减少了计算量。

附图说明

图1为本发明高光谱图像波段选择方法实施例的基本步骤流程图。

图2为本发明高光谱图像波段选择装置的结构示意图。

图3为某海军机场AVIRIS高光谱图像

图4中，(a)为图3截取的100×100子图的第72波段，(b)为改子图中38个飞机目标的真实位置分布。

图5为从图4(a)中提取的7条端元光谱曲线及关键波段示意图。

图6中，(a)为关键波段间相似性可视图，(b)原始波段间相似性可视图。

图7中，(a)-(d)分别为SFS_MNBS、Key_SFS_MNBS、MI和Key_MI四种方法结合RX异常检测的二值化检测结果。

图8为MI波段选择+RX异常检测所得结果的ROC曲线。

图9为Key_MI波段选择+RX异常检测所得结果的ROC曲线。

图10为SFS_MNBS波段选择+RX异常检测所得结果的ROC曲线。

图11为Key_SFS_MNBS波段选择+RX异常检测所得结果的ROC曲线。

具体实施方式

为了使本发明的目的、技术方案和优点更加清楚，下面结合具体的实施例详细说明本发明。下面描述了具体实施例以简化本发明。但是需要认识到，本发明不局限于所说明的实施例，并且在不脱离基本原理的前提下，本发明的各种修改是可能的，这些等价形式同样落于本申请所附权利要求书所限定的范围。

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

如图1所示，为所发明的高光谱图像波段选择方法实施例的基本步骤流程图。如图2所示，为所发明的高光谱图像波段选择装置的结构示意图。

实施例一

步骤一：使用HFC(Harsanyi-Farrand-Chang)方法确定端元个数。

(1)计算得到图像数据的协方差矩阵K_L×L和自相关矩阵R_L×L。

(2)分别求出协方差矩阵和自相关矩阵的特征值集合，记为{λ₁≥λ₂≥…λ_L}和其中L为光谱波段数。

(3)求得光谱图像的第第l个波段的近似噪声方差值其中M×N表示图像中元素的个数。

(4)计算概率密度函数

(5)给定虚警概率P_F，根据和求得τ_l值

(6)满足的特征值个数，就是所求的波段数。

步骤二：使用FNSGA(Fast NewSimplexGrowingAlgorithm)单形体增长算法实现端元提取，得到端元光谱曲线。

(1)对图像中每个像元r_i(i＝1,2,3,...,M×N)，计算其空间像素纯度指数(spatialpixelpurityindex，SPPI)，且像元r_i的空间像素纯度指数用P_i表示，那么将具有最小值的P_i所对应的像元作为第一个端元，即e₁＝r_id(1)，其中此时计数器n＝1；

(2)对图像中每一个像元r_i(i＝1,2,...,M×N)，计算其得到的单形体体积将得到最大单形体体积V的像元作为第二个端元e₂，即：e₂＝r_id(2)，此时n＝2；

(3)循环步骤：当n≥2时，每增加一个像元r_i，构造矩阵A_n-1＝[e₂-e₁,e₃-e₁,...,e_n-e₁]及A_n ⁽ⁱ⁾＝[A_n-1 r_i-e₁]，并且计算体积其中，det(A^TA)根据下式计算：

寻找具有最大单形体体积V的像元作为第n个端元，即，e_n＝r_id(n)，

(4)若n<p，则n＝n+1，跳至循环(3)，否则终止循环，得到最终目标端元{e₁,e₂,e₃,..,e_p}。

步骤三：对每个端元的光谱曲线使用提出的基于三点的关键点提取算法，提取关键波段，构造关键波段集合。

光谱曲线个数为p，每条曲线具有L个数据点，关键波段的下标集合为

(1)计数器i＝1；

(2)输入第i条光谱曲线(实际为一维行向量)，e_n＝[d₁,d₂,d₃,..,d_L],d_j为光谱曲线上第j个数据点，令计数器j＝1；

(3)取d_j、d_j+1、d_j+2的数据存入z_q-1、z_q、z_q+1，计算三点所成的夹角θ，计算表达式如下：

(4)根据下述判断准则判断第j点是否为关键点，若是则KeyS＝KeyS∪{j}，若否，则不做改变；

{(((z_q-1<z_q)∩(z_q+1<z_q))∪((z_q-1>z_q)∩(z_q+1>z_q)))∩(θ>α)}∪(θ>β)，并且abs(z_q+1-z_q)+abs(z_q+1-z_q)>τ，其中，θ为所述q-1、q、q+1三点形成的夹角，α、β为决定提取出的关键点个数的预设参数，且满足0<α<β≤π，τ为反映变化幅度的参数。

(5)若j<L，则j＝j+1，跳至第(3)步，否则进入第(6)步

(6)若i<p，则i＝i+1，跳至第(2)步，否则终止循环；

(7)合并各条端元曲线提取的关键点，得到最终关键波段号集合KeyS。

步骤四：根据关键波段子集中各波段间的相似性度量可视图，确定最佳波段数范围。

假设关键波下标集合KeyS的大小为k，关键波段集合记作Φ_2d＝{B₁,B₂,…,B_k}∈R^MN×k，其中MN＝M×N。

(1)将波段B_i归一化后得到，

(2)t＝k/3，并取整，初始化信息散度矩阵C为k×k的全0矩阵。

(3)计算

(4)计算

(5)画C矩阵的值对应的二维图像，根据图像对角线为中心的大小方块数确定最佳波段数范围。

步骤五：对关键波段子集KeyS对应的数据，使用基于MNBS准则的SFS波段选择方法，进行波段选择。

假设关键波段子集中有k个波段，k个波段数据记作Φ＝{D₁,D₂,…,D_k}∈R^M×N×k，任意设定要选择的波段数，使其属于步骤四的波段数范围，记为num(num<k)。

(1)计算Φ对应的噪声矩阵N，对于Φ中每个波段D_i∈R^M×N×1，其对应噪声矩阵为noise_i∈R^M*N*1，noise_i中各项通过下式计算：

其中D_i(m,n)和noise_i(m,n)分别代表第i波段图像和对应的噪声图像中第m行第n列的元素，有m＝1,2,3...,M，n＝1,2,3,...,N。

(2)得到Φ和N的二维形式，即将一副M×N的图像矩阵存入MN×1的一维向量中，得到新的数据Φ_2d＝[B₁,B₂,…,B_k]∈R^MN×k和N_2d＝[n₁,n₂,…,n_k]∈R^MN×k，计算原数据Φ和噪声数据N的协方差矩阵Σ_Φ∈R^k×k和Σ_N∈R^k×k，即Σ_Φ＝Φ_2d ^T*Φ_2d，Σ_N＝N_2d ^T*N_2d。设初始已选波段为令计数器j＝1。

(3)对于所有k个波段，计算每个波段的图像质量指数Q，第i个波段band_i的图像质量指数其中将具有最大图像质量指数的波段作为最优波段，对应波段号加入已选波段集合S，即：S＝S∪{id(1)}，令j＝2。

(4)当j≥2时，计算所有k个波段的图像质量指数Q。对于第i个波段band_i的图像质量指数Q_i，根据下式计算：

其中，

将具有最大图像质量指数的波段作为最优波段，对应波段号加入已选波段集合S，即：S＝S∪{id(j)}，

(5)若j<n，j＝j+1，并且跳至第(4)步。否则，终止循环，得到最终选择出的波段集合KeyS(S)。

实施例二

实施例二的步骤一至四与实施例一完全相同。

步骤五：对关键波段子集对应的数据，使用基于最大信息量的波段选择方法，进行波段选择。

假设关键波段集合中有k个波段，k个波段数据记作Φ_2d＝{B₁,B₂,…,B_k}∈R^MN×k，其中MN＝M×N。波段选择所需波段数量设为num(num<k)。

对步骤四得到的C矩阵

(1)计算maxC＝max(C)，定义集合

(2)计算C矩阵每一行元素的最小值，记为minR_i，(i＝1,2,3,...,k)；

(3)计算并调整SET＝SET∪g；

(4)修改C矩阵的元素使C(g,i)＝C(g,i)＝maxC；

(5)如果SET中元素数目少于k-num，跳至第(3)步，否则跳至第(6)步；

(6)从关键波段子集KeySet中移除SET中的元素，剩余的即为最后选择的波段。

下面以真实高光谱图像为例说明具体的实施方式。

采用的图像是加利福尼亚州-圣地亚哥地区某海军机场的AVIRIS高光谱图像，如图3所示，大小是400×400像元,空间分辨率3.5m,除去水汽吸收波段，保留189个波段。用ENVI软件对原图像进行截取，选取100×100的子图，其第72个波段的灰度图如图4(a)所示，选取后的区域中包含了38个待检测目标，图4(b)为目标的真实位置分布。

根据实施例步骤一，得到100×100子图中端元个数为7，图5为经过实施例步骤二和三得到的端元光谱曲线和关键波段示意图，其中曲线上标记‘o’的点横坐标值为关键波段号。

图6为实施例步骤四得到的波段相似性度量矩阵的可视图，其中图6(a)由关键波段子集计算得到，图6(b)为189个原始波段计算所得。从图6可以看出波段间的相似性具有块状分布的特点，统计对角线上的典型方块数目，可以确定最佳的波段数范围为4-14。

为便于描述，下面先给出各种算法的简称。基于MNBS准则的前向搜索算法用于原始波段进行波段选择，对应算法简称为SFS_MNBS，实施例一所描述的完整算法简称为Key_SFS_MNBS。基于最大信息量的波段选择算法用于原始波段进行波段选择，对应算法简称为MI，实施例二所描述的完整算法简称为Key_MI。用SFS_MNBS、Key_SFS_MNBS、MI和Key_MI四种算法对100×100子图进行波段选择后，再进行RX异常检测。

图7给出了四种波段选择方法结合RX异常检测的二值化结果图，其中(a)-(d)分别对应SFS_MNBS、Key_SFS_MNBS、MI和Key_MI方法。图7结果表明Key_SFS_MNBS的结果比SFS_MNBS的结果稍微好一点，同样Key_MI的结果也略优于MI的结果。

为了更加直观的比较各种波段选择对RX异常检测结果的改进，接下来采用接收者操作特性曲线(Receiver Operating Characteristic Curve，简称ROC曲线)来比较它们的检测性能。ROC曲线是由很多检测概率和虚警概率对应点所构成的曲线，它反映了检测概率(Pd)随虚警概率(Pf)变化而变化的关系，能够提供对算法检测性能的定量分析。ROC曲线中两个衡量算法总体检测性能的重要参数：ROC曲线下面积(area under the ROC curve,AUC)和ROC曲线上距离理想点(0,1)最小距离的平方Δγ，它们反映了检测算法的鲁棒性和平均灵敏性。图8-11分别给出了所选波段数分别为2，4，6，8，10，12，14，20和30时，SFS_MNBS、Key_SFS_MNBS、MI和Key_MI四种算法结合RX异常检测所得结果的ROC曲线。表1给出了对应ROC曲线的AUC和指标比较。图8-11和表1的结果表明：(1)关键波段提取能提升波段选择的效果；(2)波段间相似性度量的可视图可以比较准确地确定最佳波段数范围。

表1不同波段选择方法的RX检测性能比较

Claims

1.一种基于关键波段提取的高光谱图像波段选择方法，其特征在于具体实现步骤如下：

步骤1、确定高光谱图像典型地物种类，即端元个数p；

步骤2、对高光谱图像提取p个端元光谱；

2.根据权利要求1所述的一种基于关键波段提取的高光谱图像波段选择方法，其特征在于步骤3所述的基于三点矢量夹角分析的方法具体包括：

3-2.计算q-1、q、q+1三点形成的夹角θ，其中0＜θ≤π；

3.根据权利要求2所述的一种基于关键波段提取的高光谱图像波段选择方法，其特征在于所述的步骤3-2的夹角θ的计算方法为：

θ = \arccos | \frac{1 + (z_{q + 1} - z_{q}) (z_{q} - z_{q - 1})}{\sqrt{1 + {(z_{q} - z_{q - 1})}^{2}} \sqrt{1 + {(z_{q + 1} - z_{q})}^{2}}} | .

4.根据权利要求2所述的一种基于关键波段提取的高光谱图像波段选择方法，其特征在于所述的步骤3-3的判断准则如下：

当所述q-1、q、q+1三点的纵坐标值z_q-1、z_q、z_q+1满足{(((z_q-₁＜z_q)∩(z_q+1＜z_q))∪((z_q-₁＞z_q)∩(z_q+1＞z_q)))∩(θ＞α)}∪(θ＞β)，并且abs(z_q+1-z_q)+abs(z_q+1-z_q)＞τ，则q点为关键点，其中，夹角θ由所述q-1、q、q+1三点形成，α、β为决定提取出的关键点个数的预设参数，且满足0＜α＜β≤π，τ为反映变化幅度的参数。

5.根据权利要求2所述的一种基于关键波段提取的高光谱图像波段选择方法，其特征在于步骤4所述的最佳波段数范围的确定方法具体如下：

6.根据权利要求5所述的一种基于关键波段提取的高光谱图像波段选择方法，其特征在于步骤4-1所述的波段间相似性度量值的度量指标选用K-L散度；N个像元的高光谱图像第i个波段数据表示为B_i＝[b_i1，b_i2，…b_iN]^T，将波段B_i归一化后得到计算第j个波段图像相对于第i个波段图像的K-L散度第j个波段和第i个波段间的相似性度量值为

7.根据权利要求1所述的一种基于关键波段提取的高光谱图像波段选择方法的实现装置，其特征在于包括端元数确定模块、端元提取模块、关键波段子集构造模块、相似性度量矩阵可视图绘制模块以及最优波段子集确定模块；

所述的端元提取模块，用于对高光谱图像提取p个端元光谱；

8.根据权利要求7所述的一种基于关键波段提取的高光谱图像波段选择方法的实现装置，其特征在于所述高光谱图像波段选择装置还包括应用模块，所述应用模块使用波段选择结果来进行高光谱图像目标检测或分类。

9.根据权利要求7所述的一种基于关键波段提取的高光谱图像波段选择方法的实现装置，其特征在于所述关键波段子集构造模块构造侯选波段子集的具体步骤包括：

(1)获取第i条端元光谱曲线上第q-1、q、q+1三点的光谱值，即曲线纵坐标值z_q-1、z_q、z_q+1，其中q＝2，3，...，L-1；其中L为波段数；

(2)计算q-1、q、q+1三点形成的夹角θ，其中0＜θ≤π；计算方法为：

θ = \arccos | \frac{1 + (z_{q + 1} - z_{q}) (z_{q} - z_{q - 1})}{\sqrt{1 + {(z_{q} - z_{q - 1})}^{2}} \sqrt{1 + {(z_{q + 1} - z_{q})}^{2}}} |;

当所述q-1、q、q+1三点的纵坐标值z_q-1、z_q、z_q+1满足{(((z_q-1＜z_q)∩(z_q+1＜z_q))∪((z_q-1＞z_q)∩(z_q+1＞z_q)))∩(θ＞α)}∪(θ＞β)，并且abs(z_q+1-z_q)+abs(z_q+1-z_q)＞τ，则q点为关键点，其中，夹角θ由所述q-1、q、q+1三点形成，α、β为决定提取出的关键点个数的预设参数，且满足0＜α＜β≤π，τ为反映变化幅度的参数；

(4)利用每一条光谱曲线的关键波段集合KeyBand₁，KeyBand₂，...，KeyBand_p构造候选波段集合KeySet，所述候选波段集合满足KeySet＝KeyBand₁∪KeyBand₂∪KeyBand₃∪...∪KeyBand_p。