假设高光谱图像的波段数为B和像元数为N, 则可以使用X_o={x_i }i=1N∈ R^{B× N}表示待检测高光谱图像上的所有像元即观测像元。 高光谱图像上的像元在大概率上为背景像元, 而异常像元出现的概率为小概率事件。

假设背景像元可以用其他的背景像元来线性组合表示, 且在实际中用于表示的组合系数往往会是低秩的。 因此引入Latent LRR(latent low-rank representation)^[9]模型解决该问题, 可以得到如下的表达式(1)

minZrank(Z) s.t Xo=[Xo, Xh]Z(1)

这里[X_o, X_h]表示字典, Z∈ R^{N× N}是表示矩阵, X_h表示不存在于待检测高光谱图像上的像元, 而是通过光谱组成原理可合成的潜在像元即未观测像元, rank(· )表示矩阵的秩。 通过Liu的理论3.1^[9]可知观测像元和未观测像元同样来自于低秩的空间。 并推导出通过利用两者的信息获得更准确的背景表示, 获得表达式(2)

minZ, Trank(Z)+rank(T)s.t.X=XZ+TX(2)

这里简化方程X=X_o, T∈ R^{B× B}。 其中XZ为背景信息的主要特征, TX为背景信息中的显著特征。 通过Liu的论文^[9]可以证明rank(Z)和rank(T)这两部分是自动平衡的, 无需进行参数调整。 对于高光谱图像, 这意味着同时利用光谱相似的像元和光谱维度上典型波段进行背景像元表示。

由于异常像元属于无法用线性组合表示的部分, 通过l_{2, 1}范数引入异常项‖ E‖ _{2, 1}= ∑j=1n∑i=1mEij2, 其中E∈ R^{B× N}。 该项假设异常像元是被损坏的正常像元, 这将会激励E的列朝0收缩导致异常的稀疏性, 意味着异常像元是发生的小概率事件。 求解rank(· )是一个“ NP难(NP-hardness)” 问题, 可以通过求解核范数‖ · ‖ _* 获得近似的解。 进一步获得Latent LRR的表达式(3)

minZ, T, E‖Z‖* +‖T‖* +τ‖E‖2, 1s.t. X=XZ+TX+E(3)

这里τ 作为平衡异常像元部分的系数, 通过调整它的大小来控制异常的概率。

尽管Latent LRR模型通过引入观测数据和未观测数据[X_o, X_h]作为字典, 以全局的信息对数据样本进行表示。 但是在实际中, 空间中的数据样本往往可以由全局中某些典型的数据样本线性组合而成, 而大部分的数据样本在表示中仅仅占据少量成分。 对应于高光谱图像, 像元光谱特征往往是由少数典型的像元近似地表示^[8]。 所以引入l₀范数约束背景表示矩阵Z, 这会使系数矩阵在所有列上都表现出稀疏的特性, 将l₀范数放松到l₁范数近似求解。

minZ, T, E‖Z‖* +‖T‖* +β‖Z‖1+τ‖E‖2, 1s.t. X=XZ+TX+E(4)

这里β 作为平衡背景表示矩阵Z的稀疏性。

当假设邻近像元邻近关系为线性的时候, 高光谱像元上的每个像元都来自潜在的子流形M中。 这些像元的局部几何结构可以用他周围邻近像元重构而成。 本工作将会构造一个带有n个顶点的k最邻近图G, 高光谱每个像元表示顶点。 通过构造对称的权重矩阵W_ij, 他的值表示第i像元同第j像元的相连的权重系数。 构成一个无向图G=(V, E; W), 其中V表示顶点的集合, E表示边的集合, 权重矩阵W的取值如式(5)。

Wij=1, xi∈Nknn(xj)orxj∈Nknn(xi)0, otherwise(5)

这里N_knn(x_i)表示像元x_i的最近的knn个邻像元。

根据高光谱像元在数据分布的内部几何体越靠近, 则它们在重新投影到新的空间中也将会相互靠近的。 所以通过保护顶点对之间的仿射变换的低维向量来描述图G的顶点, 他们之间的相似性可以通过边的权重来衡量, 可以表示为

min∑ij‖zi-zj‖2Wij(6)

这里z_i和z_j表示像元x_j和x_j通过某些变换的投影。 拉普拉斯矩阵可以通过L=D-W得到, 其中度矩阵D是一个对角矩阵, 它的对角线上的值为d_ij= ∑jW_ij。 可以推导出

min∑ij‖zi-zj‖2Wij=mintr(ZTLZ)(7)

将拉普拉斯正则项加入模型中

minZ, T, E‖Z‖* +‖T‖* +β‖Z‖1+τ‖E‖2, 1+αtr(ZTLZ)s.t. X=XZ+TX+E(8)

其中α 为控制拉普拉斯约束的权重。 称该模型为拉普拉斯约束低秩表示。 运用该模型进行高光谱图像异常检测。

从数据中恢复低秩信息的优化问题有很多解决方法, 交替方向方法ADM(alternating direction method)^[11]是其中最有效的工具。 为了高效解决该模型, 运用自适应惩罚的线性交替方向算法LADMAP(linearized alternating direction method with adaptive penalty)^[12]来解目标方程。

在这个目标方程中, 引入辅助变量来帮助分离目标。 优化问题可以重新表述为

minZ, J, T, E‖Z‖* +‖T‖* +β‖J‖1+τ‖E‖2, 1+αtr(ZTLZ)s.t. X=XZ+TX+E, Z=J(9)

该方程是约束优化问题, 引入两个拉格朗日乘子Y₁和Y₂将这两个线性等式约束移除。

f(Z, T, J, E, Y1, Y2)=‖Z‖* +‖T‖* +β‖J‖1+τ‖E‖2, 1+αtr(ZTLZ)+< Y1, X-XZ-TX-E> +< Y2, Z-J> +μ2(‖X-XZ-TX-E‖F2+‖Z-J‖F2)(10)

其中μ 作为惩罚参数, 根据策略进行调整。 主要通过固定其他变量求解一个变量, 交替迭代进行。

(1)固定J, E, 求解Z, T

minZ, T‖Z‖* +‖T‖* +q(Z, T, J, E, Y1, Y2)(11)

其中定义

q(Z, T, J, E, Y1, Y2)=αtr(ZTLZ)+μ2‖X-XZ-TX-E+1μY1‖F2+‖Z-J+1μY2‖F2(12)

算法1 高光谱异常检测算法

输入: X, 参数β > 0, τ > 0, α > 0, knn> 0

输出: 异常程度分布图D(· )

初始化: Z₀=J₀=Y_{1, 0}=0∈ R^{N× N}, E₀=Y_{2, 0}=0∈ R^{B× N}, T₀=0∈ R^{B× B}, μ ₀=0.05, μ _max=10¹⁰, ρ ₀=1.2, ε ₁=2× 10^-4, ε ₂=1× 10^-5, k₌0, η _Z=η _T=‖ X ‖22

1. 通过计算每个x_i的knn最近邻来构造W, 得到拉普拉斯L=D-W。

2. 循环开始, 当‖ X-XZ_k-T_kX-E_k‖ _F/‖ X‖ _F≥ ε ₁或者μ _kmax(‖ Z_k-Z_k+1‖ _F, ‖ T_k-T_k+1‖ _F, ‖ J_k-J_k+1‖ _F, ‖ E_k-E_k+1‖ _F)/‖ X‖ _F≥ ε ₂。

3. 根据式(17)更新Z_k+1

4. 根据式(18)更新T_k+1

5. 根据式(19)更新J_k+1

6. 根据式(21)更新E_k+1

7. 更新拉格朗日乘子Y_{1, k+1}和Y_{2, k+1}

Y1, k+1=Y1, k+μk(X-XZk+1-Tk+1X-Ek+1)Y2, k+1=Y2, k+μk(Zk+1-Jk+1)

8. 更新μ _k+1

9. μ _k+1=min(μ _max, ρ μ _k)这里

ρ=ρ0如果μkmax(‖Zk-Zk-1‖F, ‖Tk-Tk-1‖F, ‖Jk-Jk-1‖F, ‖Ek-Ek-1‖F)/‖X‖F≤ε21, otherwise

10. 更新迭代次数k=k+1

11. 结束循环

12. 输出

D(xi)=∑j([E* ]ij)2

该问题可以分成两个问题求解

minZ‖Z‖* +ηZμk2‖Z-(Zk+1ηZqZ)‖F2(13)minT‖T‖* +ηTμk2‖T-(Tk+1ηTqT)‖F2(14)

其中q_Z, q_T为跟Z和T相关的子梯度。

qZ=-αμk(L+LT)Zk+XTX-XZk-TkX-Ek+Y1, kμk-Zk-Jk+Y2, kμk(15)qT=X-XZk+1-TkX-Ek+Y1, kμkXT(16)

使用奇异值阈值算子(SVT)Θ (· )可以获得闭合解。

Zk+1=Θ(ηZμk)-1(Zk+qZ/ηZ)(17)Tk+1=Θ(ηZμk)-1(Tk+qT/ηT)(18)

(2)固定Z, T, E, 求解J。

同样固定其余变量, 来优化求解 minJβ ‖ J‖ ₁+ μk2‖ J-(Z_k+1+ 1μkY2, k) ‖F2。

Jk+1=Sβ/μk(Zk+1+Y2, k/μk)(19)

求这个子问题可以有如下近似解。 其中S(· )是收缩算子(shrinkage)。

(3)固定Z, T, J, 求解E。

同样固定其余变量, 来优化求解

minEτ‖E‖2, 1+μk2‖E-(X-XZk+1-Tk+1X+1μkY1, k)‖F2(20)

求这个子问题可以有如式(21)近似解

Ek+1=Ωτ/μkX-XZk+1-Tk+1X+1μkY1, k(21)

其中Ω (· )是l_{2, 1}最小化X算子。

算法1总结了本方法的步骤, 其中初始化参数时, Z₀, J₀, Y_{1, 0}, E₀, Y_{2, 0}以及T₀矩阵作为优化问题的初始点, 可以设为零矩阵或者随机矩阵。 μ ₀和μ _max是控制拉格朗日算子更新的速率, ρ ₀则是控制μ 的变化, 这三个参数根据迭代速度设置。 ε ₁和ε ₂决定迭代的误差的大小, 根据需要检出异常的量级来设置大小。 k表示算法迭代的次数。 η _Z和η _T是归一化数。 算法1将会输出图像中异常程度分布图。

使用仿真数据和真实数据来验证方法。

仿真数据来自于罗切斯特理工大学(RIT)在美国纽约大都会区罗切斯特市附近飞行采集的实验数据集, 波长为400~2 500 nm, 共360个波段, 全部波段用于实验。 空间分辨率为1 m。 用于合成仿真图像的数据分为背景数据和异常数据:背景数据取长为100像元, 宽为100像元, 波段为360的图像块, 如图1(a)所示; 异常数据为五种材料, 分别是黄色的毛毡(Yellow Felt), 黄色的棉布(Yellow Cotton), 粉色的毛毡(Pink Felt), 金色的毛毡(Gold Felt)和蓝色的毛毡(Blue Felt)的光谱, 如图1(c)所示。 合成仿真数据过程中, 将利用异常的光谱以不同的成分大小同背景数据进行线性融合, 如式(22)

S(i, j)=f×Sa+(1-f)×Sbf∈(0.1, 0.3, 0.5, 0.8, 1)(22)

如图1(b)所示, 以每行对应不同的材料光谱S_a和图中从左到右的(i, j)位置光谱的S_b, 按每列为异常占的比重f的大小, 从小到大进行合成S(i, j)后放置于同一位置(i, j)。

图1

Fig.1

	Figure Option ViewDownloadNew Window
	图1 (a)仿真数据的伪彩图; (b)真实分布图; (c)异常的光谱曲线Fig.1 (a) False color image of the simulation data; (b) Ground-truth map; (c) The spectral profile of anomaly

第一份真实数据同样来自罗切斯特理工大学(RIT)数据。 取其中长宽为100像元, 波段360的图像块, 如图2(a)。 图中的异常目标为事先铺设好的布, 异常点占所有像元的1.35%, 他的真实分布图如图2(b); 第二份真实数据是来自AVRIS上采集的美国加利福尼亚州的圣地亚哥市飞机场。 图像块由波长为370~2 510 nm采集的224个波段组成的。 空间分辨率为3.5 m。 在去除掉对应的水汽吸收区域、 信噪比低和损坏的波段后(1— 9, 101— 115, 151— 179和217— 224), 剩余183个波段用于实验。 空间上取长为100像元, 宽为100像元的图像块进行测试, 如图3(a)所示。 图中飞机作为异常数据进行检测, 异常点占所有像元的0.62%, 它的真实分布图如图3(b); 第三份真实数据是来自机载的HYDICE光谱数据。 该数据的空间分辨率为每像元1米。 波长分辨率为10 nm的210个波段数据, 去除水汽吸收区域、 低信噪比以及损坏的波段(1— 4, 76, 87, 101— 111, 136— 153, 198— 210)后剩下的长宽分别为100像元、 80像元的162个波段, 如图4(a)所示。 图中的异常为汽车, 占全部像元的0.21%, 他的真实分布图如图4(b)。

图2

Fig.2

	Figure Option ViewDownloadNew Window
	图2 (a)罗切斯特理工大学数据的伪彩图和(b)真实分布图Fig.2 (a) False color image of the RIT data; (b) Ground-truth map

图3

Fig.3

	Figure Option ViewDownloadNew Window
	图3 (a)圣地亚哥数据的伪彩图和(b)真实分布图Fig.3 (a) False color image of the Sandiego data; (b) Ground-truth map

图4

Fig.4

	Figure Option ViewDownloadNew Window
	图4 (a)都市数据的伪彩图和(b)真实分布图Fig.4 (a) False color image of the Urban data; (b) Ground-truth map

实验将分析和对比不同的六种算法在四类数据中的检测性能。 由检测的目视结果图定性的评价算法的效果。 图5中色彩条表明在目视结果图中, 色彩越接近蓝色端判定为背景概率越高, 而色彩越接近红色端则意味着判定为异常概率越高。 除此之外, 还分别绘制出Global-RX(Global Reed-Xiaoli), DWEST, Seg-RX, RPCA-RX, LRASR和本算法在各个数据集上的ROC(receiver operating characteristic)曲线, 并统计出各自AUC(area under the curve)值进行定量的对比。 结果表明提出的算法明显突出异常像元的分布, 且在同样的误检率(false alarm rate)下, 将获得更高的检测率(probability of detection)。

图5

Fig.5

	Figure Option ViewDownloadNew Window
	图5 异常色彩条Fig.5 Colorbar for anomaly

对于仿真数据, 算法对比的目视结果图如6(a)— (f)所示。 可以看出随着异常成分所占比率的升高, 本算法在准确的位置检出更高的异常值, 这表明异常更显著。 并且本算法在异常值低于0.3以下时依然可以检出异常, 较其他的五种算法相比明显效果更佳。 不同算法的ROC曲线如图6(g)所示, 可以看出本算法的ROC曲线包裹住其他方法的曲线即说明新算法在性能上更加优异。 特别是在0.01~0.3的误检率之间, 本方法获得了更好的检测率, 远远高于另外五种算法。 新算法在误检率为0.2左右时, 能检出所有的异常值。

图6

Fig.6

	Figure Option ViewDownloadNew Window
	图6 不同方法的对仿真数据的检测结果和ROC曲线Fig.6 The detection results and ROC curves by different methods for simulation data

对于罗切斯特理工大学的数据, 图7(a)— (f)所示为六种算法对于他的目视结果图。 相较于Global-RX, DWEST, Seg-RX传统算法, 基于低秩原理的另外其他三种方法获得更好的检测效果; 本算法通过未观测数据表示和低维流形的约束, 对于罗切斯特理工大学的数据中房子部分的抑制获得更好的效果。 减少房子的检出将会在同样误检率下提高检测率。 ROC曲线如图7(g)所示也表明了本算法的效果。

图7

Fig.7

	Figure Option ViewDownloadNew Window
	图7 不同方法的对RIT数据的检测结果和ROC曲线Fig.7 The detection results and ROC curves by different methods for RIT data

圣地亚哥数据集中, 六种算法的目视结果图如图8(a)— (f)所示。 由于地物背景种类很多导致RPCA-RX效果不明显, 是因为RPCA-RX本身是基于单一子空间的算法, 对于位于多子空间上的数据无法获得较好的检测率。 而LRASR和本算法是基于低秩表示的, 对于多子空间中的数据能更好的检测出异常。 LRASR和本算法在目视结果中没有太多的差别。 通过他们的ROC曲线分析, 如图8(g)所示。 在0.02~0.14误检率之间新算法检测率明显较LRASR更高。 主要因为LRASR依赖于字典的构建, 他无法避免字典构建过程中, 异常的意外引入导致的字典污染, 致使不正确的表示。 而新算法通过全部高光谱图像上观测数据和潜在的未观测数据表示, 无需进行额外的字典训练过程, 并引入拉普拉斯约束限制异常像元参与待检测像元的表示, 所以效果更好。

图8

Fig.8

	Figure Option ViewDownloadNew Window
	图8 不同方法的对圣地亚哥数据的检测结果和ROC曲线Fig.8 The detection results and ROC curves by different methods for Sandiego data

图9(a)— (f)所示的是六种算法对于城市数据的效果图和图9(g)所示是其ROC曲线图。 图中异常的比例较低, 本算法依然有效的提高了异常的检出率。 图中可以看出公路与草地等背景均具有低秩的特点, 在新算法中将会抑制具有低秩属性背景来突出图中异常像元。

图9

Fig.9

	Figure Option ViewDownloadNew Window
	图9 不同方法的对城市数据的检测结果和ROC曲线Fig.9 The detection results and ROC curves by different methods for Urban data

表1是对六种算法对不同数据在ROC曲线下的面积AUC对比。 其中AUC的值越接近百分百, 说明检测器的性能更加优异。 本算法将检测率稳定在99%以上, 几乎达到了100%的程度, 相对于其他的算法提高了大约1%2%的AUC。 本检测算法比其他五个传统算法优异。

表1

Table 1

表1(Table 1)

表1 对于不同数据用不同算法获得的AUC(%) Table 1 AUC obtained by different methods for different datasets(%) 算法 仿真 RIT 圣地亚哥 城市 GlobalRX 88.75 98.75 87.98 98.92 DWEST 91.36 88.38 95.23 91.70 SegRX 79.84 76.30 93.25 98.65 RPCARX 95.71 99.02 96.91 98.70 LRASR 93.06 98.20 97.15 99.56 LRLLRR 99.18 99.75 99.04 99.88

表1 对于不同数据用不同算法获得的AUC(%) Table 1 AUC obtained by different methods for different datasets(%)