CN106919952A - 基于结构稀疏表示和内部聚类滤波的高光谱异常目标检测方法 - Google Patents

Abstract

本发明公开了一种基于结构稀疏表示和内部聚类滤波的高光谱异常目标检测方法，用于解决现有高光谱异常目标检测方法目标检测效率低的技术问题。技术方案是在选择初始背景像元后，利用主成分分析字典学习法学习得到鲁棒性的背景字典。稀疏向量求解和图像重建过程中，引入重加权拉普拉斯先验，提高稀疏向量求解精度。计算原始图像与重建图像之间的误差得到稀疏表示误差。利用内部聚类滤波表示高光谱数据的空谱特性，通过计算待测像元与它同类其他像元线性表示结果的误差，得到内部聚类误差，最后联合稀疏表示误差与内部聚类误差的线性加权和实现异常目标的精确提取。经测试，本发明相对于背景技术在恒虚警率的前提下检测率提高了10％～15％。

Description

基于结构稀疏表示和内部聚类滤波的高光谱异常目标检测 方法

技术领域

本发明涉及一种高光谱异常目标检测方法，特别是涉及一种基于结构稀疏表示和内部聚类滤波的高光谱异常目标检测方法。

背景技术

高光谱异常目标检测技术是一种不需要提供待测目标先验光谱信息的前提下，仅利用图像像元之间的光谱差异性进行目标检测的高光谱目标检测技术，在实际应用中具有较强的实用性。

相对于整个高光谱图像，异常目标不仅出现的概率较低，而且所占的比例也较少，传统的异常目标检测算法一般假设图像背景服从高斯分布。在这种假设下可以利用全局或局部的统计特性来检测异常目标。然而，在实际应用中由于受到空间分辨率限制，图像背景具有较高的复杂性，仅通过假设背景服从高斯分布来进行异常目标检测，不能取得很好的检测效果。

近年来，基于稀疏表示的高光谱异常目标检测算法受到越来越多的关注。该类算法假设图像中的像元可以通过背景和目标的字典线性稀疏表示，在求得稀疏向量后，利用重建误差的大小提取异常目标。在实际应用中具有很好的灵活性和可扩展性。文献“Li J,Zhang H,Zhang L,et al.Hyperspectral anomaly detection by the use ofbackground joint sparse representation[J].Selected Topics in Applied EarthObservations and Remote Sensing,IEEE Journal of,2015,8(6):2523-2533.”公开了一种基于背景联合稀疏表示的高光谱异常目标检测方法。该方法首先以全部图像像元作为字典，通过稀疏表示模型求解稀疏向量，根据所求的非零稀疏向量位置确定局部背景像元，然后根据局部背景像元建立局部正交子空间提取异常目标。文献所述方法在选取局部背景像元时每次都使用了大规模的全局像元作为字典，给求解稀疏向量带来较大的计算负担；在求解稀疏向量时没有考虑稀疏向量内部的结构稀疏信息，造成求解稀疏向量不准确；未有效考虑高光谱图像中像元之间的空间约束，目标检测效率不高。

发明内容

为了克服现有高光谱异常目标检测方法目标检测效率低的不足，本发明提供一种基于结构稀疏表示和内部聚类滤波的高光谱异常目标检测方法。该方法在稀疏表示理论框架下，首先基于局部RX算法选择初始背景像元后，利用主成分分析字典学习法学习得到鲁棒性的背景字典。稀疏向量求解和图像重建过程中，引入重加权拉普拉斯先验，提高稀疏向量求解精度。计算原始图像与重建图像之间的误差得到稀疏表示误差。此外，利用内部聚类滤波表示高光谱数据的空谱特性，通过计算待测像元与它同类其他像元线性表示结果的误差，得到内部聚类误差，最后联合稀疏表示误差与内部聚类误差的线性加权和实现异常目标的精确提取。在真实的高光谱卫星图像AVIRIS、HYDICE和仿真的高光谱数据集上的试验结果表明，本发明获得的检测结果相对于已有的异常目标检测方法在恒虚警率的前提下检测率提高了10％～15％。

本发明解决其技术问题所采用的技术方案是：一种基于结构稀疏表示和内部聚类滤波的高光谱异常目标检测方法，其特点是包括以下步骤：

步骤一、输入一幅高光谱图像，选用局部双窗RX算法获取背景像元集合。局部RX算法计算公式如下：

D_LocalRX(x_i)＝(x_i-μ_local(i))^T(∑_local(i))^-1(x_i-μ_local(i)) (1)

式中，x_i表示第i个像元，μ_local(i)表示x_i像元所对应的局部背景的均值，(∑_local(i))^-1表示x_i像元所对应的局部背景的协方差矩阵的逆矩阵，利用公式(1)遍历所有高光谱像元后选取阈值进行背景分割，小于阈值的对应像元为背景像元，得到背景像元矩阵Z_local，大小为n_b×N_Z，n_b表示高光谱数据波段个数，N_Z表示局部背景像元的个数。所述阈值取0.007～0.01；

步骤二、根据主成分分析学习算法，对步骤(一)得到的背景像元矩阵Z_local转置后求其协方差矩阵，定义为CovZ_local，大小为n_b×n_b。然后求协方差矩阵CovZ_local的特征值V和特征向量P，将特征值按照从大到小排列后，相应的特征向量也依次排列，得到的新的特征向量矩阵即为最终学习到的背景字典D。

步骤三、在稀疏表示框架下，高光谱图像像元能被学习到的字典D稀疏表示，即

X＝DY+N (2)

式中，X为高光谱图像，D为步骤(二)学习到的背景字典，Y＝[y₁,…,y_n]为稀疏表示矩阵，N为稀疏表示误差和图像噪声。将重加权拉普拉斯稀疏先验引入Y，得到基于重加权拉普拉斯稀疏先验的结构稀疏表示优化模型。

步骤四、采用坐标下降法求解优化模型得到最优稀疏信号Y_rec，然后通过公式得到重建后的高光谱图像X_rec。根据原图像与重建图像得到重建误差，重建误差计算公式如下：

r_sprase(X)＝||X-X_rec||₂ (3)

其中，X、X_rec分别表示原图像和重建后图像，r_sprase(X)表示基于稀疏表示的重建误差，||·||₂表示向量的2范数。

步骤五、采用k-means++算法对输入的高光谱数据进行聚类。假设输入图像X被聚为M个类，表示第m个聚类包含的像元集合，n_b代表图像波段数量，n_m代表第m个聚类有n个像元。

步骤六、计算类内像元之间的相似权值，计算公式如下：

式中，表示第m个聚类里第i、j个像元，h是尺度因子，为归一化因子。

步骤七、针对每一个图像像元，用与之同类型的其他像元在相似权值的基础之上进行线性表示，表示公式如下：

式中，表示第m个聚类中第i个元素被同一聚类中其他像元线性加权表示后的结果。

步骤八、计算原来像元与线性表示后的像元之间的表示误差。如果该像元是背景像元，则能够被该类其他背景像元线性表示，表示误差较小；否则，若该像元为异常像元，则表示误差较大。据此能够提取异常目标。表示误差的计算公式如下：

式中，X、分别是原图像和聚类表示后的图像，r_cluster(X)表示基于内部聚类的重建误差，||·||₂表示向量的2范数。

步骤九、加权求和稀疏重建误差和聚类误差，最终得到异常目标检测结果。计算公式如下：

r_final＝0.5×r_sparse+0.5×r_cluster (7)

在两个重建误差加权求和所得的灰度图基础之上，选取分割阈值δ，图像中大于δ的像素标记为1表示目标，小于等于δ的像素标记0表示背景，δ的取值范围为0.1～0.3。最终得到只有0和1的二值结果图，完成异常目标检测。

本发明的有益效果是：该方法在稀疏表示理论框架下，首先基于局部RX算法选择初始背景像元后，利用主成分分析字典学习法学习得到鲁棒性的背景字典。稀疏向量求解和图像重建过程中，引入重加权拉普拉斯先验，提高稀疏向量求解精度。计算原始图像与重建图像之间的误差得到稀疏表示误差。此外，利用内部聚类滤波表示高光谱数据的空谱特性，通过计算待测像元与它同类其他像元线性表示结果的误差，得到内部聚类误差，最后联合稀疏表示误差与内部聚类误差的线性加权和实现异常目标的精确提取。在真实的高光谱卫星图像AVIRIS、HYDICE和仿真的高光谱数据集上的试验结果表明，本发明获得的检测结果相对于已有的异常目标检测方法在恒虚警率的前提下检测率提高了10％～15％。

下面结合具体实施方式对本发明作详细说明。

具体实施方式

本发明基于结构稀疏表示和内部聚类滤波的高光谱异常目标检测方法具体步骤如下：

假设输入的高光谱图像是一个三维数据立方体，包含n_b个波段，每个波段是一幅n_row行和n_col列大小的图像。为了计算方便，将每个波段拉伸成为一个行向量，所有行向量组成一个二维矩阵X，其中，X的每一列表示每一个像素对应的光谱，该方向为光谱维；X的每一行对应一个波段的所有像素值(即n_p＝n_row×n_col)，该方向为空间维。具体如下：

1、基于主成分分析的鲁棒性背景字典学习。

(1)采用双窗局部RX算法获取背景像素集合。

根据输入高光谱图像分辨率选择窗口尺寸都为奇数的矩形外窗Outer和内窗Inner(Outer＝n_outer×n_outer,Inner＝n_inner×n_inner,n_outer>n_innner,n_outer、n_inner分别表示外窗和内窗的尺寸，都为奇数)。以每一个输入像素x_i∈X,i＝1,…,n_p为中心，根据外窗和内窗尺寸获取局部背景区域n_local＝Outer-Inner，n_local为n_b×N_local大小的局部背景矩阵，然后计算n_local的均值和协方差矩阵分别得到μ_local＝[μ₁,…,μ_N]^T和∑_local。局部RX计算公式如下：

D_LocalRX(x_i)＝(x_i-μ_local(i))^T(∑_local(i))^-1(x_i-μ_local(i)) (1)

式中，μ_local(i)表示第i个输入像元所对应的局部背景的均值，(∑_local(i))^-1表示第i输入像元所对应的局部背景的协方差矩阵的逆矩阵，利用公式(1)遍历所有高光谱像元后选取阈值进行背景分割，小于阈值的对应像元为背景像元，得到背景像元矩阵Z_local，大小为n_b×N_Z，N_Z表示局部背景像元的个数。为了得到纯净的背景像元，选择较小的阈值较为合适，建议取值为0.007～0.01。

(2)基于主成分分析的鲁棒性背景字典学习。

将背景像元矩阵Z_local转置后求其协方差矩阵，定义为CovZ_local，大小为n_b×n_b。然后求协方差矩阵CovZ_local的特征值V和特征向量P(每一个特征值对应一个特征向量)，将特征值按照从大到小排列后，相应的特征向量也依次排列，得到的新的特征向量矩阵即为最终学习到的背景字典D。

2、建立基于重加权拉普拉斯稀疏先验的结构稀疏表示模型。

在稀疏表示框架下，高光谱图像像元可以被学习到的字典D稀疏表示，即

X＝DY+N (2)

D为已学习到的背景字典，Y＝[y₁,…,y_n]为稀疏表示矩阵，N为稀疏表示误差和图像噪声。

假设N服从的矩阵正态分布，则对应的似然函数为

其中，Σ_n＝diag(λ)是以λ元素为对角线元素的对角矩阵，用于表示误差和噪声的强度。表示Q矩阵的加权迹范数。

为了表示稀疏向量内部的结构稀疏性，本发明将重加权拉普拉斯稀疏先验引入Y。首先假设Y服从如下分布

其中，Σ_y＝diag(γ)表示以γ的元素为对角线元素的对角矩阵，控制Y中每一行的稀疏度，γ_i＝0表示Y的第i行为0。假设则其中任一列y_i服从高斯分布。假设超参数γ服从以下的伽马分布，

以上两级先验，等价于重加权拉普拉斯分布，因为对于y_i有

其中，为保证求解便利，本发明中不直接使用重加权拉普拉斯先验，而是采用级联先验，而且λ，γ和κ均为待估计参数。

由于λ，γ和κ未知，无法采用最大后验估计(Maximum a posterior estimation，MAP)直接对稀疏信号Y进行估计。因此，本发明根据经验贝叶斯框架，先基于输入的高光谱数据X利用MAP估计未知参数λ，γ和κ，如下

其中，λ_opt,γ_opt,κ_opt分别表示最优的λ，γ，κ。通过积分，并引入-2log运算，容易得知式子(7)等价于最小化如下的式子

其中，tr(·)表示迹范数，Σ_by＝Σ_n+DΣ_yD^T，为代价函数。

通过变形式子(8)，得到稀疏信号Y的非分离稀疏约束模型。首先，对式子(8)的第一份部分进行变形

然后，将式子(9)带入到式子(8)中，得

接着，引入新的代价方程如下

显然，而且可以证明，最小化式子(8)再对稀疏信号Y进行MAP估计，与直接最小化式子(11)得到的λ，γ和κ相同，关于Y的解仅相差一个常量因此，式子(11)可以看作是关于稀疏信号Y的正则化回归模型，其中为稀疏信号的非分离稀疏约束。该约束不能拆分成对于Y中每一行的独立约束，因此该约束能同时约束稀疏信号中非零元素，潜在地考量非零元素之间的相关性。此外，Σ_by中包含了表征噪声强度的λ，因此，得到的稀疏约束可以随着估计的噪声强度自适应的变化，具有噪声鲁棒性。

最终得到了如下的基于重加权拉普拉斯稀疏先验的结构稀疏表示优化模型：

3、模型求解并重建高光谱图像。

已知待检测高光谱数据X，本发明采用坐标下降法求解式子(12)，每次迭代中仅优化一个变量而固定剩余的其他变量。λ⁰，γ⁰，κ⁰分别表示初始的数值，t记录迭代次数，η表示更新阈值，λ^t，γ^t，κ^t，Y^t分别表示迭代到第t次时的数值，λ^t+1，γ^t+1，κ^t+1，Y^t+1分别表示迭代到第t+1次时的数值。

具体步骤如下：

①初始化，λ⁰，γ⁰，κ⁰均初始化为对应长度的全1向量，计数变量t＝0；

②更新中间变量Σ_n＝diag(λ^t)，Σ_y＝diag(γ^t)，Σ_by＝Σ_n+DΣ_yD^T；

③固定λ^t，γ^t和κ^t，根据式子(11)得到关于Y的优化形式，如下

求解得到Y的更新规则如下，

④固定Y^t+1，λ^t和κ^t，得到关于γ的优化形式，如下

求解得到如下的更新形式：

其中，为γ^t+1的第i个元素，代表V^T+Y^t+1(Y^t+1)^T的对角线元素组成的向量，z_i为z的第i个元素；

⑤固定Y^t+1，γ^t+1和κ^t，得到关于λ的优化形式，如下

求解得到如下的更新形式：

其中，根号运算表示向量每一个元素开方后组成的向量，./运算代表两个向量对应元素相除后组成的向量，代表对角线元素组成的向量。

⑥固定Y^t+1，γ^t+1和λ^t+1，得到关于κ的优化形式，如下:

求解得到如下的更新形式：

上式子中的加法和除法运算均作用在向量的每一个元素上，得到一个新向量，d＝10^-6的引入是为了确保γ中出现0时，式子(20)依然有意义。

⑦计算稀疏信号Y更新前后的差异，如下

其中，表示对Y^t+1内的每一个元素乘以||·||_F表示弗罗贝尼乌斯范数(Frobenius norm),如果计数器t>5或者更新差异η<10^-3，则退出循环；否则t+1，循环执行②至⑦。

⑧假设上述循环结束得到的最优稀疏信号为Y_rec，则待重建的高光谱图像X_rec可以通过如下方式得到：

4、计算重建误差。

显然X_rec是基于背景字典重建的结果，由此可知原高光谱图像中背景像素可以被精确重建，而异常像素不能被精确重建。重建误差计算公式如下：

r_sprase(X)＝||X-X_rec||₂ (23)

其中，X、X_rec分别表示原图像和重建后图像，都为n_b×n_p二维矩阵，||·||₂表示向量的2范数。

5、基于内部聚类滤波算法的高光谱图像聚类。

(1)利用k-means++算法对图像聚类。

针对输入高光谱图像X，设置聚类个数K的值(根据不同图像数值不同，通常范围为30～50)，然后进行以下具体步骤：

①随机选择一个像元作为第一个聚类中心；

②对于图像中的每一个像元x，计算它与最近聚类中心(指已选择的聚类中心)的距离D(x)；

③选择一个新的像元作为新的聚类中心，选择原则是D(x)较大的点，被选取作为下一个聚类中心；

④重复步骤②和③直到K个聚类中心被选出来；

⑤遍历每个像元计算它与k个聚类中心的距离，该像元与哪个聚类中心最近就将它划分到该类，这样将所有图像像元聚集为K类；

⑥分别计算由⑤得到的K个聚类的质心，得到K个新的聚类中心；

⑦重复步骤⑤和⑥，直到质心不再变化，得到将原图像X聚为K个类的结果。假设输入图像X被聚为M个类表示第m个聚类包含的像元集合，n_b代表图像波段数量，n_m代表第m个聚类有n个像元。；

(2)计算类内像元之间的相似权值，计算公式如下：

式中表示第m个聚类里第i、j个像元，h是尺度因子，为归一化因子。

(3)针对每一个图像像元，用与之同类型的其他像元在相似权值的基。

础之上进行线性表示，表示公式如下：

式中，表示第m个聚类中第i个元素被同一聚类中其他像元线性加权表示后的结果。

(4)计算原来像元与线性表示后的像元之间的表示误差。

如果该像元是背景像元，则可以被该类其他背景像元线性表示，表示误差较小；否则，若该像元为异常像元，则表示误差较大。据此可以提取异常目标。表示误差的计算公式如下：

式中，X、分别是原图像和聚类表示后的图像。

6、加权求和稀疏重建误差和聚类误差，最终得到异常目标检测结果。通过上述计算，分别得到了结构稀疏表示和内部聚类误差，通过加权

计算两者误差之和，最终分割出异常目标，得到检测结果。计算公式如下：

r_final＝0.5×r_sparse+0.5×r_cluster (27)

根据灰度图r_final选取分割阈值δ，图像中大于δ的像素标记为1表示目标，小于等于δ的像素标记0表示背景，δ的建议取值范围为0.1～0.3。最终得到只有0和1的二值结果图，完成异常目标检测，结束。

Claims (1)

1.一种基于结构稀疏表示和内部聚类滤波的高光谱异常目标检测方法，其特征在于包括以下步骤：

步骤一、输入一幅高光谱图像，选用局部双窗RX算法获取背景像元集合；局部RX算法计算公式如下：

D_LocalRX(x_i)＝(x_i-μ_local(i))^T(∑_local(i))^-1(x_i-μ_local(i)) (1)

式中，x_i表示第i个像元，μ_local(i)表示x_i像元所对应的局部背景的均值，(∑_local(i))^-1表示x_i像元所对应的局部背景的协方差矩阵的逆矩阵，利用公式(1)遍历所有高光谱像元后选取阈值进行背景分割，小于阈值的对应像元为背景像元，得到背景像元矩阵Z_local，大小为n_b×N_Z，n_b表示高光谱数据波段个数，N_Z表示局部背景像元的个数；所述阈值取0.007～0.01；

步骤二、根据主成分分析学习算法，对步骤(一)得到的背景像元矩阵Z_local转置后求其协方差矩阵，定义为CovZ_local，大小为n_b×n_b；然后求协方差矩阵CovZ_local的特征值V和特征向量P，将特征值按照从大到小排列后，相应的特征向量也依次排列，得到的新的特征向量矩阵即为最终学习到的背景字典D；

步骤三、在稀疏表示框架下，高光谱图像像元能被学习到的字典D稀疏表示，即X＝DY+N(2)

式中，X为高光谱图像，D为步骤(二)学习到的背景字典，Y＝[y₁,…,y_n]为稀疏表示矩阵，N为稀疏表示误差和图像噪声；将重加权拉普拉斯稀疏先验引入Y，得到基于重加权拉普拉斯稀疏先验的结构稀疏表示优化模型；

步骤四、采用坐标下降法求解优化模型得到最优稀疏信号Y_rec，然后通过公式得到重建后的高光谱图像X_rec；根据原图像与重建图像得到重建误差，重建误差计算公式如下：

r_sprase(X)＝||X-X_rec||₂ (3)

其中，X、X_rec分别表示原图像和重建后图像，r_sprase(X)表示基于稀疏表示的重建误差，||·₂表示向量的2范数；

步骤五、采用k-means++算法对输入的高光谱数据进行聚类；假设输入图像X被聚为M个类，表示第m个聚类包含的像元集合，n_b代表图像波段数量，n_m代表第m个聚类有n个像元；

步骤六、计算类内像元之间的相似权值，计算公式如下：

${\mathrm{\&omega;}}_{ij}=\frac{1}{W_i}\exp (-||x_m^i-x_m^j|{|}_2^2/h)---\left(4\right)$

式中，表示第m个聚类里第i、j个像元，h是尺度因子，为归一化因子；

步骤七、针对每一个图像像元，用与之同类型的其他像元在相似权值的基础之上进行线性表示，表示公式如下：

${\hat{x}}_m^i=\underset{j\&NotEqual;i}{\overset{n_m}{\&Sigma;}}{\mathrm{\&omega;}}_{ij}{x}_m^j---\left(5\right)$

式中，表示第m个聚类中第i个元素被同一聚类中其他像元线性加权表示后的结果；

步骤八、计算原来像元与线性表示后的像元之间的表示误差；如果该像元是背景像元，则能够被该类其他背景像元线性表示，表示误差较小；否则，若该像元为异常像元，则表示误差较大；据此能够提取异常目标；表示误差的计算公式如下：

$r_{cluster}\left(X\right)=\left|\right|X-\hat{X}|{|}_2---\left(6\right)$

式中，X、分别是原图像和聚类表示后的图像，r_cluster(X)表示基于内部聚类的重建误差，||·||₂表示向量的2范数；

步骤九、加权求和稀疏重建误差和聚类误差，最终得到异常目标检测结果；计算公式如下：

r_final＝0.5×r_sparse+0.5×r_cluster (7)

在两个重建误差加权求和所得的灰度图基础之上，选取分割阈值δ，图像中大于δ的像素标记为1表示目标，小于等于δ的像素标记0表示背景，δ的取值范围为0.1～0.3；最终得到只有0和1的二值结果图，完成异常目标检测。