2. 西安科技大学 通信与信息工程学院, 陕西 西安 710054
高光谱遥感是在电磁波谱的可见光、近红外、中红外和远红外波段范围内，获取许多非常窄的连续光谱影像数据的技术[1]。高光谱目标检测是高光谱图像处理中的重要研究方向之一，它利用丰富的光谱信息，通过对目标和背景信息的分析和处理，可以有效地将弱小的地物目标从背景图像中探测出来[2]。通常目标检测算法建立在一定先验信息的基础上，而实际应用中，先验信息的获取十分困难。因此，利用无监督方式来对异常目标进行检测，具有广泛的实用价值，其研究和发展也受到了大量的关注。异常检测实际上是一种二分类问题，将高光谱图像分为目标和背景，根据异常目标光谱和周围背景光谱特征的差异(这种差异常常是微弱的)进行信息提取。基于这种差异的高光谱图像异常检测，在森林火灾检测、矿产勘探、医学肿瘤检查、军事目标监测等领域具有重要应用价值[3-4]。
异常检测主要依靠全局或局部统计变化来检测异常。异常检测算法主要包括RX(Reed-Xiaoli)算法[5-8]、核RX(KRX)算法[9-10]、基于支持向量描述(SVDD)的算法[11-13]、基于稀疏表示的算法[14-15]以及基于协同表示的算法[16]。近年来深度学习技术在图像分类中取得了前所未有的卓越性能，Li等[17]在迁移学习基础上，建立了基于CNN(convolution neural network，卷积神经网络)框架的高光谱异常检测算法，通过比对测试点和近邻域像素来获取异常目标，但其只考虑了光谱差异。
随着成像技术的发展，成像光谱仪的光谱分辨率和空间分辨率不断提高。高光谱数据不仅包含了地物的光谱辐射信息还包含了地物的空间分布、形状以及纹理等空间信息。异常检测的目的是寻找与背景有明显差异的目标，这种差异不仅体现在光谱信息上还体现在空间结构上(如纹理特性)。然而传统异常检测算法大多只考虑异常点与背景像素的光谱差异、忽略了二者在空间结构上的差异，这在一定程度上造成了高光谱数据的浪费，也使得检测性能未能得到充分地发挥。
鉴于此，本文提出一种基于局部空间结构差异的空谱联合高光谱图像异常检测方法。该方法，不仅计算待测像素与背景窗内像素在光谱维度上的差异程度，还度量了内窗与背景窗在空间结构上的差异。由于异常检测中，目标形状、纹理、大小等信息都是未知的，存在尺寸很小甚至为亚像元的目标，传统的检测器对小目标的检测效果较差。因此本文又构建了自适应多层异常检测框架，较传统的双窗模型(全局RX、局部RX和KRX)更好地利用了局部空间结构和光谱维度信息，降低了虚警率，提升了对像素少的异常目标的检测效果。
1 局部空间结构差异的异常检测算法 1.1 算法模型异常检测算法多基于Reed和Xiaoli Yu提出的RX双窗模型[5]，通过检测当前像素与背景像素的光谱差异来发现异常点，但却往往忽略了目标和背景在空间结构上的关联性。因此，本文提出一种基于局部空间结构差异的空谱联合异常检测算法。该算法不仅考虑光谱信息，还通过内窗和背景窗的空间结构相似性计算异常值，同时调整空间异常值和光谱异常值的权值，得到最终的异常指数。该算法模型如图 1所示。
|
| 图 1 局部空间结构差异的空谱联合算法模型 |
在检测过程中, 采用相同中心点的双窗模型。图 1中, 红色代表待测像素; 绿色内窗大小为rin×rin, 可以看作目标区域, 起到保护待测像素的作用; 灰色部分的外窗大小为rout×rout, 可以看作背景区域。通常rin和rout取奇数, rout=K·rin(K∈R+)。为了更好地描述待测像素与背景区域间的异常性, 设置异常指数(anomaly index, AI)进行评估。
1.2 光谱异常值将高光谱遥感图像各个像元每条光谱曲线看作是光谱空间的一个矢量, 通过计算待测像元光谱矢量与背景光谱矢量之间的角度大小来表示两者之间的相似度, 进一步度量光谱异常值(spectral anomaly value, SAV)。
待测像素r(ik, jk)={r1, r2, …rL}T, 背景窗像素rb(ik, jk)={rb1, rb2, …rbL}T, 均表示为L×1维列向量(L为总波段数), 两者之间的相似函数定义为
|
(1) |
式中，||·||为2范数, 用于求解向量的长度。由高光谱数据某一像素点在所有波段的成像性质可知, sim(r, rb)∈[0, 1], 待测像素与背景像素的光谱波形越相似, 夹角越小, 相似函数值越大。由于异常值需要计算待测像元和背景像元间的“相异程度”, 对(1)式变形得到
|
(2) |
f(r, rb)的值越大则待测像素和背景像素的光谱差异越大。计算待测像素与背景窗口内所有像素点之间的相异程度, 取均值即可得到该待测像素与外窗内背景像元之间的光谱异常值
|
(3) |
本文定义了局部空间结构的异常值(locally spatial structure anomaly value, LSSAV), 将待测像素同内窗看作一个整体, 以体现出待测点的空间结构信息。
如图 1所示, 在当前第k波段中, 将以待测像素位置(ik, jk)为中心点的局部空间结构窗口表示为Sik, jk, 像素归一化后的灰度值表示为g(Sik, jk), 该窗口可以理解为传统RX算法中的内窗结构(图 1中绿色和红色方块组成部分)。搜索窗中心点位置为
使用l2距离对图像中2个空间结构相似性进行度量[18]。定义不同窗口的空间结构之间的异常值为
|
(4) |
计算完当前窗口之后, 搜索窗按照图 1所示的搜索路径, 以l为步长进行滑动。每次滑动后的计算结果作为当前2个空间结构的异常值, 最终取当前像素(ik, jk)所有空间结构异常值中的最小值作为待测像素内窗空间结构异常值
|
(5) |
将待测像素内窗与最相似结构的欧氏距离作为异常指数的评价准则, 充分考虑了异常目标“小”的特性, 减弱了利用均值时大面积的背景对评价准则的影响。δLSSAV越大, 当前像素(ik, jk)的局部空间结构异常度越高; 反之, 当前像素的空间结构异常度越低。
1.4 异常指数第k波段待测像素(ik, jk)的异常指数δ(ik, jk)由光谱异常值δSAV和局部空间结构的异常指数δLSSAV组成
|
(6) |
式中，w为[0, 1]之间的常数, 作为权重调整空间异常值和光谱异常值对异常检测结果的影响, 待测像素遍历整幅图像即可得到当前第k波段的异常检测结果。对于波段数为L的高光谱数据, 将各个波段的异常指数进行累加, 最终得出该算法当前像素(i, j)的异常检测结果为
|
(7) |
在空谱联合检测的基础上本文又构建了基于背景抑制的自适应多层结构的异常检测算法框架(anomaly detection algorithm based on adaptive multi-layers structure, ADAML), 进一步提高异常检测的性能。ADAML框架主要由异常检测层、背景抑制层和判决停止层构成, 如图 2所示。
|
| 图 2 基于自适应多层结构的异常检测算法框架 |
大小为M×N×L的高光谱数据, 每个波段内像元总数为P(P=M×N), 总波段数为L。该数据可以看作一个P×L的矩阵X=[x1, x2, …xL], 每行代表该波段的所有数据。这里使用上文提出的空谱联合异常检测器作为异常检测层, 采用图 1所示的双窗模型, 第n层检测器的输出为
|
(8) |
计算出的yn是[0, 1]上的实值, 大小为M×N, 表示各个像素点的异常指数, 即为当前的检测结果。
2.2 背景抑制层第n层异常检测层的输出yn, 经过非负函数激活后与该层的输入xn相乘, 在背景抑制层对背景进行抑制, 构造出xn+1作为第n+1层检测器的输入值
|
(9) |
非线性激活函数q(·)表示如下
|
(10) |
式中, λ为正常数, 通过λ数值可以调整q(·)的形状, 同时改变对于xn的约束力度。对于背景光谱, yn∈[0, 1)并接近于0, 每个光谱向量xn与(yn)λ相乘后会快速下降; 对于目标光谱, yn接近于1, xn与(yn)λ相乘后基本保持不变。因此, 通过对背景光谱的非线性约束, 背景光谱和目标光谱之间的差异逐渐扩大, 经过背景抑制层处理后, 结果如下
|
(11) |
在多层约束之后, 背景被抑制后的结果仅包含值为0的向量, 异常目标谱仍然为原始的光谱矢量值。虽然该算法包含多层异常检测器, 但是随着层数的增加, 数据的稀疏性也随之增加, 除了第一层以外, 计算速度逐层加快, 具有较高的计算效率。由图 2可以看出背景抑制层的作用。
2.3 判决停止层如果人为设置异常检测的层数n, 不仅需要不断调整n值使之达到较优结果, 而且易陷入局部最优解, 计算效率较低。因此, 需要设置一个判决函数, 使该算法能够自适应确定层数n, 在达到某种情况时, 停止计算, 同时输出当前层的异常检测结果。本文中的判决函数使用第n层和第n+1层的平均能量之差
|
(12) |
在理想情况中, η=0, 但在实际实验过程中, η为趋近于0的一个正值。当平均能量之差大于η时, 继续下一层迭代; 否则, 停止检测。本文通过设计了判决停止层, 自适应确定层数n, 从而达到良好的异常检测效果。
算法流程如图 3所示:
|
| 图 3 ADAML算法流程图 |
实验数据1和实验数据2来自AVIRIS传感器获取的San Diego海军机场图像, 该高光谱数据集在光谱范围370~2 510 nm之间共有224个波段, 将低信噪比和水汽吸收波段剔除后还有189个波段。原始图像大小为400×400, 为使检测效果更直观, 截取左上角100×100的子图像作为实验数据1, 灰度图如图 4a)所示, 目标分布如图 4b)所示, 共57个目标像素点。选取原始图像中间大小为200×200的子图作为实验数据2, 灰度图像如图 4c)所示, 目标分布如图 4d)所示, 共有134个目标像素点。实验数据3来自ROSIS传感器拍摄的Pavia Center, 在光谱范围430~860 nm之间共有115个波段, 去除低信噪比和水汽吸收波段后, 还有102个波段。原始图像大小为1 096×715, 截取的子图像横坐标范围在223~322之间, 纵坐标范围在1~96之间, 该子图像大小为100×96, 其灰度图像如图 4e)表示, 目标分布如图 4f)所示, 共有33个目标像素点。
|
| 图 4 实验数据图示 |
针对本文提出的空谱联合算法设计了多组对比实验, 调整空域异常值权重(w), 得到ROC曲线如图 5所示。由图 5可以看出, 当空域异常值权重在[0.3, 0.5]之间效果较好, 越远离这个区间, 效果越差。当空域异常值所占权重为0.9时, 检测的结果最差。由此说明: 高光谱数据空间和谱域的信息都很重要, 不能仅考虑单一方面的信息, 要有效结合空域和谱域信息, 从而得到更优的检测结果。
|
| 图 5 不同w取值下ROC曲线对比(实验数据1) |
内窗作为待测像素的空间结构窗, 内窗和外窗之间的区域作为待测像素的背景。外窗既要能反映像素的局部信息, 又要能体现背景的结构信息, 因此不能太大也不能太小。综合考虑, 令rout=3×rin。
调整内窗边长rin, 通过计算ROC的面积AUC, 验证窗口大小对检测效果的影响, 其中w=0.4, 步长l设置为1, 实验结果如表 1所示。
| rin | 实验数据1 | 实验数据2 | 实验数据3 | |||
| AUC | 时间/s | AUC | 时间/s | AUC | 时间/s | |
| 3×3 | 0.990 1 | 15.92 | 0.862 9 | 20.07 | 0.989 9 | 12.01 |
| 5×5 | 0.998 3 | 17.33 | 0.886 9 | 23.43 | 0.989 1 | 13.38 |
| 7×7 | 0.995 2 | 19.09 | 0.924 5 | 26.96 | 0.989 0 | 14.47 |
| 9×9 | 0.994 3 | 23.77 | 0.904 4 | 29.32 | 0.989 0 | 15.42 |
| 11×11 | 0.993 9 | 26.98 | 0.893 3 | 34.44 | 0.988 7 | 16.91 |
| 13×13 | 0.993 3 | 30.38 | 0.892 9 | 39.49 | 0.988 0 | 19.02 |
可见, 不同大小的窗口影响异常检测的结果, 如果窗口大小与异常目标尺寸相近或略大于目标尺寸时, 异常检测的性能较好, 窗口过大或过小, 都会影响异常检测的效果。综合AUC和运算时间, 3组数据的内窗大小分别设置为5×5, 7×7, 3×3。
3.3 空谱融合异常检测算法的仿真实验结果图 6至8分别显示了采用G-RX、L-RX、KRX和本文提出的空谱联合4种异常检测算法对实验数据1、实验数据2和实验数据3的异常检测效果二值图。各异常检测算法的运行时间如表 2所示。
|
| 图 6 不同算法在实验数据1上的检测结果对比 |
|
| 图 7 不同算法在实验数据2上的检测结果对比 |
|
| 图 8 不同算法在实验数据3上的检测结果对比 |
| 算法 | 实验数据1 | 实验数据2 | 实验数据3 |
| GRX | 2.567 7 | 3.564 9 | 1.783 2 |
| LRX | 29.970 7 | 41.114 3 | 11.734 5 |
| KRX | 395.897 2 | 576.825 0 | 265.940 1 |
| 空谱联合 (本文算法) |
17.33 | 26.96 | 12.01 |
上述实验结果表明: G-RX算法的检测结果虚警率较高, 目标更易受到背景干扰; L-RX算法能在一定程度上减弱背景的影响, 使得异常目标更为显著, 但是计算效率较低; KRX方法将非线性的像素投影到更高维度的空间中, 将其转变为线性可分的情况, 有效利用了高光谱数据的高阶统计信息, 能够较好地检测出目标的轮廓, 在实验数据1和实验数据3中都有较好的表现, 但是在实验数据2中, 却并没有将背景和目标进行很好分离, 而该方法的计算量比L-RX算法更大, 需要更多的运算时间; 本文提出的空谱联合算法可以更好地利用空间结构信息和光谱域的信息, 更加准确地检测异常目标, 并且虚警率较低, 无需计算逆矩阵, 运行时间较短。
3.4 ADAML算法的仿真实验由图 6至8可以看出, 空谱联合算法虽然可以检测出异常目标, 但是对尺寸较小的异常目标检测时仍存在较大的虚警率。利用自适应多层结构框架ADAML分别在实验数据1、实验数据2和实验数据3上进行异常检测仿真实验, 仿真结果的ROC曲线如图 9所示。实验中, η=10-5, λ=3。
|
| 图 9 空谱联合算法和ADAML算法的ROC曲线对比 |
由图 9的ROC曲线可以看出, 在使用ADAML结构后, ADAML的ROC曲线处于仅使用单层空谱联合算法结构的上方, 说明本文提出的自适应多层结构对图像背景有抑制作用, 可以更好地分离目标和背景, 降低虚警率。而在相同虚警率的条件下, 提升了检测率。
4 结论本文针对尺寸较小的异常目标, 提出了一种自适应多层结构的异常检测新方法。该方法主要由异常检测层、背景抑制层和判决停止层构成。其中, 异常检测层采用了空谱融合的异常检测器, 将局部空间结构异常值与光谱域异常值相结合, 构造出新的空谱联合异常指数; 背景抑制层通过结合上一层检测器的输出结果, 非线性地约束背景向量, 增大目标与背景的差异; 判决停止层设置算法停止条件, 自适应完成背景抑制, 避免人为设置参数导致性能不稳定和计算效率低等问题。实验结果表明, ADAML算法通过抑制背景信号, 提高了异常检测效果。
| [1] |
童庆禧, 张兵, 郑兰芬. 高光谱遥感: 原理、技术与应用[M]. 北京: 高等教育出版社, 2006.
TONG Qingxi, ZHANG Bin, ZHENG Lanfen. Hyperspectral remote sensing: the principle, technology and application[M]. Beijing: Higher Education Press, 2006. (in Chinese) |
| [2] | MANOLAKIS D, SHAW G. Detection algorithm for hyperspectral imaging applications[J]. IEEE Signal Processing Magazine, 2002, 19(1): 29-43. DOI:10.1109/79.974724 |
| [3] | STEIN D, BEAVEN S, HOFF L, et al. Anomaly detection from hyperspectral imagery[J]. IEEE Signal Processing Magazine, 2002, 19(1): 58-69. DOI:10.1109/79.974730 |
| [4] | MATTEOLI S, DIANI M, CORSINI G. A tutorial overview of anomaly detection in hyperspectral images[J]. IEEE Aerospace and Electronic Systems Magazine, 2010, 25(7): 5-28. DOI:10.1109/MAES.2010.5546306 |
| [5] | REED I, YU X. Adaptive multiple-band CFAR detection of an optical pattern with unknown spectral distribution[J]. IEEE Trans on Acoustics Speech and Signal Processing, 1990, 38(10): 1760-1770. DOI:10.1109/29.60107 |
| [6] | LI W, DU Q. Decision fusion for dual-window-based hyperspectral anomaly detector[J]. Journal of Applied Remote Sensing, 2015, 9(1): 097297-097297. DOI:10.1117/1.JRS.9.097297 |
| [7] | CHEN S, WANG W, WU C, et al. Real-time causal processing of anomaly detection for hyperspectral imagery[J]. IEEE Trans on Aerospace and Electronic Systems, 2014, 50: 1511-1534. DOI:10.1109/TAES.2014.130065 |
| [8] | ZHAO C, WANG Y, QI B, et al. Global and local real-time anomaly detectors for hyperspectral remote sensing imagery[J]. Remote Sensing, 2015, 7(4): 3966-3985. DOI:10.3390/rs70403966 |
| [9] | KWON H, NASRABADI N. Kernel RX-algorithm: a nonlinear anomaly detector for hyperspectral imagery[J]. IEEE Trans on Geoscience and Remote Sensing, 2005, 43(2): 388-397. DOI:10.1109/TGRS.2004.841487 |
| [10] | HE M, MEI H, WU Y, YAN H. Weighted kernel-based signature subspace projection for hyperspectral target detection[C]//9th Workshop on Hyperspectral Image and Signal Processing: Evolution in Remote Sensing, 2018: 1-5 |
| [11] | BANERJEE A, BURLINA P, DIEHL C. A support vector method for anomaly detection in hyperspectral imagery[J]. IEEE Trans on Geoscience and Remote Sensing, 2006, 44(8): 2282-2291. DOI:10.1109/TGRS.2006.873019 |
| [12] | KHAZAI S, HOMAYOUNI S, SAFARI A, et al. Anomaly detection in hyperspectral images based on an adaptive support vector method[J]. IEEE Geoscience and Remote Sensing Letters, 2011, 8(4): 646-650. DOI:10.1109/LGRS.2010.2098842 |
| [13] | GURRAM P, KWON H, HAN T. Sparse kernel-based hyperspectral anomaly detection[J]. IEEE Geoscience and Remote Sensing Letters, 2012, 9(5): 943-947. DOI:10.1109/LGRS.2012.2187040 |
| [14] | CUI X, TIAN Y, WENG L, et al. Anomaly detection in hyperspectral imagery based on low-rank and sparse decomposition[C]//International Society for Optics and Photonics International Conference on Graphic and Image Processing, 2014: 9069 |
| [15] | ZHANG Y, DU B, ZHANG L, et al. A low-rank and sparse matrix decomposition-based mahalanobis distance method for hyperspectral anomaly detection[J]. IEEE Trans on Geoscience and Remote Sensing, 2016, 54(3): 1376-1389. DOI:10.1109/TGRS.2015.2479299 |
| [16] | LI W, DU Q. Collaborative representation for hyperspectral anomaly detection[J]. IEEE Trans on Geoscience and Remote Sensing, 2015, 53(3): 1463-1474. |
| [17] | LI W, WU G, DU Q. Transferred deep learning for anomaly detection in hyperspectral imagery[J]. IEEE Geoscience and Remote Sensing Letters, 2017, 14(5): 597-601. |
| [18] | EFROS A, LEUNG T, et al. Texture synthesis by non-parametric sampling[C]//Proceedings of the Seventh IEEE International Conference on Computer Vision, 1999: 1033-1038 |
2. School of Communication and Information Engineering, Xi'an University of Science and Technology, Xi'an 710054, China
