两阶段决策引导双字典构建的高光谱图像异常检测方法

1.本发明属于图像处理技术领域，进一步涉及高光谱图像异常检测，具体为一种两阶段决策引导双字典构建的高光谱图像异常检测方法，可用于智能农业、战场侦察和空间探索。


背景技术：

2.高光谱图像异常检测指在没有先验知识的情况下，找出与邻近像素光谱差异较大的像素。同高光谱图像目标检测技术相比，由于其不需要光谱先验信息，因而广泛应用于智能农业、战场侦察、空间探索等领域。
3.在文献[huyan n,zhang x,zhou h,et al.hyperspectral anomaly detection via background and potential anomaly dictionaries construction[j].ieee transactions on geoscience and remote sensing,2018,57(4):2263-2276]中，huyan等人提出了一种基于背景和潜在异常字典的高光谱图像异常检测方法，该方法通过利用联合稀疏表示模型在局部区域内的残差定义像素异常等级，并以此构建潜在异常字典，然而在构建潜在异常字典过程中，潜在异常字典容易受到背景像素的干扰，这在一定程度上影响了异常的刻画效果。


技术实现要素：

[0004]
本发明的目的是针对上述现有技术的不足，提出一种两阶段决策引导双字典构建的高光谱图像异常检测方法，用于解决现有方法中构建的潜在异常字典纯度低的问题。首先，对高光谱图像进行超像素分割，以获取的超像素为内窗，通过设计外窗，实现显著性检测，并对显著性检测结果进行后处理，得到粗糙二值图；其次，构建背景估计网络，以获取的粗糙二值图指导挑选背景像素，并将其用于背景估计网络的训练，待训练完成后，将高光谱图像送入训练好的背景估计网络中，得到重构图像，对输入的高光谱图像与重构图像的残差图像进行背景建模，计算每个像素的马氏距离，以此作为异常判别标准，即可得到异常响应图，并对异常响应图进行后处理，得到精细二值图；然后，利用获取的粗糙二值图、精细二值图以及超像素共同构建背景字典和潜在异常字典；最后，利用获取的背景字典和潜在异常字典构建低秩和稀疏表示模型，以获取的异常部分计算检测结果。本发明能够获取纯净的潜在异常字典，有利于更好描述异常特性，进而提升高光谱图像异常检测效果。
[0005]
本发明实现上述目的具体步骤，包括如下：
[0006]
(1)生成粗糙二值图：
[0007]
(1.1)对高光谱图像进行超像素分割，得到超像素集合s＝{si|i∈[1,k
*
]}，其中h、w和b分别表示高光谱图像x的高、宽和波段数，表示第i个超像素，ni表示第i个超像素中的像素数，k
*
表示超像素总数；
[0008]
(1.2)将si作为内窗，以si的最小封闭矩形框四条边为基准，分别向外扩展距离r来获取外窗，对si中第t个像素其中t∈[1,ni]，分别计算与其内外窗间所有像素的显著值，
以其平均值作为该像素的显著值遍历所有超像素得到显著性检测结果d1；
[0009]
(1.3)对显著性检测结果d1进行后处理，得到粗糙二值图
[0010]
(2)生成精细二值图：
[0011]
(2.1)搜索粗糙二值图中数值为0的元素，并记录其空间位置sp；
[0012]
(2.2)在高光谱图像x中挑选与sp相同的像素，构成背景像素集合并将其作为训练样本，其中表示cb中第i'个背景像素，q表示cb中包含的背景像素总数；
[0013]
(2.3)构建背景估计网络，并利用cb进行训练，待训练完成后，将x送入训练好的背景估计网络中，得到重构图像x
rec
；
[0014]
(2.4)对高光谱图像x与重构图像x
rec
的残差图像x
res
进行背景建模，并计算像素的马氏距离即残差图像x
res
中第i”个像素的异常判别结果：
[0015][0016]
其中，n表示高光谱图像x的像素总数；表示残差图像x
res
中第i”个像素，μ和σ分别表示x
res
的均值向量和协方差矩阵；遍历所有像素，得到检测结果
[0017]
(2.5)对检测结果d2进行后处理，并将其转化为2维矩阵，得到精细二值图
[0018]
(3)构建双字典：
[0019]
(3.1)构建背景字典：将粗糙二值图和精细二值图的和作为背景索引图b
idx
，并将b
idx
按照超像素集合s进行划分，得到划分结果其中表示b
idx
中第i个超像素，计算中所有像素之和，判断其是否为0；若为0，则记录中质心的空间位置spb，并根据spb在高光谱图像x中选择对应像素作为原子，构建背景字典；若不为0，则对应超像素不参与背景字典构建；当遍历所有超像素后，即可获得背景字典db；
[0020]
(3.2)构建潜在异常字典：将粗糙二值图和精细二值图的点乘作为潜在异常索引图a
idx
，搜索a
idx
中数值为1的元素，记录其空间位置spa，并根据spa在高光谱图像x中选择对应像素，以此作为原子，构建潜在异常字典da；
[0021]
(4)根据背景字典db和潜在异常字典da，构建基于双字典的低秩和稀疏表示模型：
[0022]
x＝dbw+das+e，
[0023]
其中，w和s分别表示背景项和异常项的表示系数，e表示噪声项；
[0024]
(5)目标函数优化求解：
[0025]
根据背景的低秩特性以及异常的稀疏特性，建立目标函数：
[0026][0027]
采用交替方向乘子法优化求解目标函数，得到w、s、e，通过对das的每一列向量求解l2范数，即可得到检测结果；其中β和λ表示平衡各项的系数，s.t.表示约束条件，||
·
||
*
表示矩阵的核范数，||
·
||1表示矩阵的l1范数，||
·
||
2,1
表示矩阵的l
2,1
范数。
[0028]
本发明与现有技术相比，具有以下优点：
[0029]
第一、由于本发明采用了基于超像素的显著性检测方法，从而能够有效缓解内外窗间异常像素对检测结果的影响，使得检测结果更为精确；
[0030]
第二、本发明采用两阶段决策联合引导潜在异常字典的构建，能够获取更加纯净的潜在异常字典，从而有效描述异常特性；同时，在两阶段决策联合引导背景字典构建过程中，通过引入超像素，确保背景字典覆盖所有的地物类别，并且地物越大，用于构建背景字典的像素越多，从而实现对背景特性的有效描述。
附图说明
[0031]
图1为本发明方法的实现流程图；
[0032]
图2为本发明方法在具体实例中的应用示意图；
[0033]
图3为本发明中使用的数据集；其中(a)表示hydice数据集及其对应标签，(b)表示pavia数据集及其对应标签，(c)是los angeles数据集及其对应标签，(d)表示san diego-i数据集及其对应标签，(e)表示san diego-ii数据集及其对应标签，(f)表示texas coast数据集及其对应标签；
[0034]
图4为本发明方法与现有方法主观效果对比图；其中(i)-(vi)分别表示hydice、pavia、los angeles、san diego-i、san diego-ii和texas coast数据集，(a)-(g)分别表示现有对比方法rx、crd、lrasr、lsmad、pab-dc、lsdm-mog和kifd，(h)表示本发明方法，(i)表示数据集对应的标签；
[0035]
图5为本发明方法与现有方法在hydice数据集上的roc曲线与背景-异常分离图；其中(a)表示roc曲线；(b)表示背景-异常分离图；
[0036]
图6为本发明方法与现有方法在pavia数据集上的roc曲线与背景-异常分离图；其中(a)表示roc曲线；(b)表示背景-异常分离图；
[0037]
图7为本发明方法与现有方法在los angeles数据集上的roc曲线与背景-异常分离图；其中(a)表示roc曲线；(b)表示背景-异常分离图；
[0038]
图8为本发明方法与现有方法在san diego-i数据集上的roc曲线与背景-异常分离图；其中(a)表示roc曲线；(b)表示背景-异常分离图；
[0039]
图9为本发明方法与现有方法在san diego-ii数据集上的roc曲线与背景-异常分离图；其中(a)表示roc曲线；(b)表示背景-异常分离图；
[0040]
图10为本发明方法与现有方法在texas coast数据集上的roc曲线与背景-异常分离图；其中(a)表示roc曲线；(b)表示背景-异常分离图。
具体实施方式
[0041]
下面结合附图和具体实施例，对本发明作进一步的详细描述。
[0042]
参照图1和图2，本发明提出的一种两阶段决策引导双字典构建的高光谱图像异常检测方法，具体按照如下步骤实现：
[0043]
步骤1：生成粗糙二值图：
[0044]
(1.1)对高光谱图像进行超像素分割，得到超像素集合s＝{si|i∈[1,k
*
]}，其中h、w和b分别表示高光谱图像x的高、宽和波段数，表示第i个超像素，ni表
示第i个超像素中的像素数，k
*
表示超像素总数。对高光谱图像x进行超像素分割的步骤包括：
[0045]
(1.1.1)设定预期超像素个数k，以作为初始超像素边长，得到初始超像素集合并挑选所有超像素的中心，以其作为初始聚类中心，得到初始聚类中心集合其中表示对
·
向下取整，表示第i0个初始超像素，k0表示初始超像素个数，且其中h和w分别表示高光谱图像的高和宽；
[0046]
(1.1.2)以第i0个超像素的初始聚类中心为中心，分别计算与其2s
×
2s范围内像素的距离其中包含空间距离d
spa
和光谱距离d
spe
，具体如下：
[0047][0048][0049][0050]
其中，表示以为中心2s
×
2s范围内的第j个像素，表示第i0个超像素中心的空间坐标，(xj,yj)表示以为中心2s
×
2s范围内第j个像素的空间坐标；＜
·
，
·
＞表示两个向量的内积，||
·
||表示向量的模，ω表示空间-光谱距离加权系数；
[0051]
(1.1.3)将所有像素分配至距离其最近的聚类，然后计算每个聚类中所有像素的均值，并以该聚类中距离均值最近的像素作为新的聚类中心；
[0052]
(1.1.4)重复步骤(1.1.2)-(1.1.3)直至聚类中心不再变化，得到超像素分割结果；
[0053]
(1.1.5)对步骤(1.1.4)中得到的超像素分割结果进行处理，即将过小的超像素合并到邻近超像素中，[参照文献：r.achanta,a.shaji,k.smith,a.lucchi,p.fua and s.s
ü
sstrunk,"slic superpixels compared to state-of-the-art superpixel methods,"in ieee transactions on pattern analysis and machine intelligence,vol.34,no.11,pp.2274-2282]，形成最终的超像素集合s＝{si|i∈[1,k
*
]}。
[0054]
(1.2)将si作为内窗，以si的最小封闭矩形框四条边为基准，分别向外扩展距离r来获取外窗，对si中第t个像素其中t∈[1,ni]，分别计算与其内外窗间所有像素的显著值，以其平均值作为该像素的显著值遍历所有超像素得到显著性检测结果d1；显著性检测的步骤包括：
[0055]
(1.2.1)计算待测像素与其内外窗间像素xm的显著值，具体如下：
[0056][0057]
[0058][0059]
其中，表示待测像素与其内外窗间第m个像素xm的显著值，和分别表示待测像素与其内外窗间像素xm的光谱距离和空间距离，c是一个常数，由于其对结果影响不大，为方便起见，此处将c设为1；＜
·
，
·
＞表示两个向量的内积，||
·
||表示向量的模，和(xm,ym)分别为待测像素与其内外窗间像素xm的空间坐标；
[0060]
(1.2.2)计算待测像素与其内外窗间所有像素显著值的均值，具体如下：
[0061][0062]
其中，m表示待测像素的内外窗间像素个数。
[0063]
(1.3)对显著性检测结果d1进行后处理，得到粗糙二值图具体是：首先利用otsu算法[otsu n.a threshold selection method from gray-level histograms[j].ieee transactions on systems,man,and cybernetics,1979,9(1):62-66]对显著性检测结果d1进行处理，得到初步二值图再将该二值图中较大的连通域去除，得到粗糙二值图
[0064][0065]
其中，area_filtering(
·
；τ)表示连通域过滤操作，当连通域面积超过指定阈值τ时，本实施例中设置该阈值为n/100，将对应区域的数值置0。
[0066]
步骤2：生成精细二值图：
[0067]
(2.1)搜索粗糙二值图中数值为0的元素，并记录其空间位置sp；
[0068]
(2.2)在高光谱图像x中挑选与sp相同的像素，构成背景像素集合并将其作为训练样本，其中表示cb中第i'个背景像素，q表示cb中包含的背景像素总数；
[0069]
(2.3)构建背景估计网络，并利用cb进行训练，待训练完成后，将x送入训练好的背景估计网络中，得到重构图像x
rec
；这里构建背景估计网络及训练细节具体如下：
[0070]
构建背景估计网络：背景估计网络包含3个子网络，分别是编码器、解码器和判别器，其中编码器和解码器共同构成自编码器，编码器和判别器共同构成生成式对抗网络。所有子网络均由4个全连接层构成，3个子网络的节点数分别为{b,512,512,nz},{nz,512,512,b}和{nz,512,512,1}，其中b表示高光谱图像的波段数，在每个全连接层后，分别包含斜率为0.2的leakyrelu激活函数。在每层激活函数后，使用采样率为0.5的dropout策略，以防止过拟合。此外，在判别器之后，加上sigmoid函数，确保输出在[0,1]。
[0071]
模型训练细节：在训练背景估计网络时，编码器和解码器的学习率均为0.001，判别器的学习率为0.0001，3个子网络的优化器均为adam，迭代次数为1000。背景估计网络的损失函数l包含两部分，分别是自编码器的损失l
ae
和生成式对抗网络的损失l
gan
，具体如下：
[0072]
l＝l
ae
+α
·
l
gan
[0073][0074][0075]
其中，x
bac
和分别表示挑选的背景像素和重构的背景像素，||
·
||2表示向量的l2范数，z
′
表示从n(0,i)中采样得到的真实样本，en(
·
)表示编码操作，d(
·
)表示判别操作，p
data
表示挑选的背景像素的分布，α表示权衡自编码器损失和生成式对抗网络损失的系数，此处α设置为1。
[0076]
(2.4)对高光谱图像x与重构图像x
rec
的残差图像x
res
进行背景建模，并计算像素的马氏距离即残差图像x
res
中第i”个像素的异常判别结果：
[0077][0078]
其中，n表示高光谱图像x的像素总数；表示残差图像x
res
中第i”个像素，μ和σ分别表示x
res
的均值向量和协方差矩阵；遍历所有像素，得到检测结果
[0079]
(2.5)对检测结果d2进行后处理，并将其转化为2维矩阵，得到精细二值图该步骤中对检测结果d2进行后处理的原理同步骤(1.3)，具体是：首先利用otsu算法对检测结果d2进行处理，得到初步二值图再将该二值图中较大的连通域去除，得到粗糙二值图
[0080][0081]
步骤3：构建双字典：
[0082]
(3.1)构建背景字典：将粗糙二值图和精细二值图的和作为背景索引图b
idx
，并将b
idx
按照超像素集合s进行划分，得到划分结果其中表示b
idx
中第i个超像素，计算中所有像素之和，判断其是否为0；若为0，则记录中质心的空间位置spb，并根据spb在高光谱图像x中选择对应像素作为原子，构建背景字典；若不为0，则对应超像素不参与背景字典构建；当遍历所有超像素后，即可获得背景字典db；
[0083]
(3.2)构建潜在异常字典：将粗糙二值图和精细二值图的点乘作为潜在异常索引图a
idx
，搜索a
idx
中数值为1的元素，记录其空间位置spa，并根据spa在高光谱图像x中选择对应像素，以此作为原子，构建潜在异常字典da；
[0084]
步骤4：根据背景字典db和潜在异常字典da，构建基于双字典的低秩和稀疏表示模型：
[0085]
x＝dbw+das+e，
[0086]
其中，w和s分别表示背景项和异常项的表示系数，e表示噪声项；
[0087]
步骤5：目标函数优化求解：
[0088]
根据背景的低秩特性以及异常的稀疏特性，建立目标函数：
[0089]
[0090]
采用交替方向乘子法优化求解目标函数，得到w、s、e，通过对das的每一列向量求解l2范数，即可得到检测结果；其中β和λ表示平衡各项的系数，s.t.表示约束条件，||
·
||
*
表示矩阵的核范数，||
·
||1表示矩阵的l1范数，||
·
||
2,1
表示矩阵的l
2,1
范数。
[0091]
所述目标函数优化求解，具体如下：
[0092]
引入辅助变量j和l，目标函数转化为：
[0093][0094]
s.t.x＝dbw+das+e,w＝j,s＝l
[0095]
对上述目标函数构建增广拉格朗日函数：
[0096][0097]
其中，y1、y2和y3表示拉格朗日乘子，μ表示惩罚参数，表示矩阵frobenius范数的平方；采用交替方向乘子法求解，即每次更新一个变量，固定其余变量：
[0098]
1)固定w、s、e和l，更新j：
[0099][0100]
其中，u,σ＝(σ1，σ2，
···
，σr)
t
,v分别表示对矩阵w+y2/μ进行奇异值分解所得的左奇异矩阵，对角矩阵的主对角元素和右奇异矩阵；
[0101]
2)固定w、s、e和j，更新l：
[0102][0103]
3)固定w、s、l和j，更新e：
[0104][0105]
其中，
[0106]
4)固定s、e、l和j，更新w：
dictionaries construction[j].ieee transactions on geoscience and remote sensing,2018,57(4):2263-2276]、lsdm-mog：[参照文献：l.lu,l.wei,d.qian,and t.ran,“low-rank and sparse decomposition with mixture of gaussian for hyperspectral anomaly detection,”ieee trans.cybern.,vol.51,no.9,pp.4363-4372,feb.2020]和kifd：li s,zhang k,duan p,et al.hyperspectral anomaly detection with kernel isolation forest[j].ieee transactions on geoscience and remote sensing,2019,58(1):319-329)从定性和定量两个角度进行分析，其中定性分析主要包括检测图、roc曲线、背景-异常分离图、定量分析采用auc值。
[0125]
图4为本发明方法与七种对比方法在六种数据集上的检测图对比结果，其中(i)-(vi)分别为hydice、pavia、los angeles、san diego-i、san diego-ii和texas coast数据集，(a)-(h)分别为rx、crd、lrasr、lsmad、pab-dc、lsdm-mog、kifd和本发明方法，通过与标签(i)对比，可以看出本发明方法效果最好。
[0126]
图5(a)、图6(a)、图7(a)、图8(a)、图9(a)和图10(a)分别为本发明方法与七种对比方法在hydice、pavia、los angeles、san diego-i、san diego-ii和texas coast数据集上的roc曲线，曲线下的面积越大，表示检测效果越好，可见本发明方法效果最好。
[0127]
图5(b)、图6(b)、图7(b)、图8(b)、图9(b)和图10(b)分别为本发明方法与七种对比方法在hydice、pavia、los angeles、san diego-i、san diego-ii和texas coast数据集上的背景-异常分离图，背景和异常间的距离越大，表明分离效果越好，即检测效果越好，可以看出本发明方法的背景与异常间距最大，即分离效果最好。
[0128]
下表(表1)列出了本发明方法与七种对比方法在六种高光谱数据集上的auc值对比结果，auc值在0～1之间，值越大表示效果越好。
[0129]
表1本发明方法与七种对比方法在不同数据集上的auc值对比(其中，最佳结果以加粗方式表示，次佳结果以阴影方式表示)
[0130][0131]
由上表可以看出，在不同数据集上本发明方法与现有方法相比均体现出了明显优势，采用本发明方法进行检测的auc值最大，即检测效果最佳。
[0132]
上述实验结果证明了本发明所提方法的正确性与有效性。
[0133]
本发明未详细说明部分属于本领域技术人员公知常识。
[0134]
以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，显然对于本领域的专业人员来说，在了解了本发明内容和原理后，都可能在不背离本发明原理、结构的情况下，进行形式和细节上的各种修正和改变，但是这些基于本发明思想的修正和改变仍在本发明的权利要求保护范围之内。