«上一篇
文章快速检索     高级检索
下一篇»
  应用科技  2017, Vol. 44 Issue (6): 20-26  DOI: 10.11991/yykj.201610010
0

引用本文  

王霖郁, 刘一博. 基于光谱角匹配加权的高光谱图像异常检测[J]. 应用科技, 2017, 44(6), 20-26. DOI: 10.11991/yykj.201610010.
WANG Linyu, LIU Yibo. Anomaly detection for hyperspectral image based on weighted spectral angle match[J]. Applied Science and Technology, 2017, 44(6), 20-26. DOI: 10.11991/yykj.201610010.

通信作者

刘一博,E-mail:413974717@qq.com

作者简介

王霖郁(1977−),女,副教授;
刘一博(1992−),男,硕士研究生

文章历史

收稿日期:2016-10-25
网络出版日期:2017-03-31
基于光谱角匹配加权的高光谱图像异常检测
王霖郁, 刘一博    
哈尔滨工程大学 信息与通信工程学院,黑龙江 哈尔滨 150001
摘要:针对高光谱背景中存在异常和噪声的问题,提出了一种基于光谱角匹配(SAM)加权的核RX异常检测算法。首先对图像背景像元进行K-均值聚类,得到不同类背景对应的聚类中心,然后计算背景像元与聚类中心的光谱角余弦,选出较纯净的背景作为新背景,最后新背景中的每个像元将自己的光谱角信息作为权值,构造加权核RX异常检测算子,通过加权削弱了残留其中的异常和噪声的干扰。为验证算法的有效性,利用真实的AVIRIS和ROSIS-03遥感器采集高光谱数据进行了仿真实验,结果表明与对比算法相比,所提算法对潜在的异常具有较强的抑制能力,提高了检测精度。
关键词高光谱图像    K-均值聚类    加权    核RX    光谱角匹配    异常检测    光谱角余弦    背景净化    
Anomaly detection for hyperspectral image based on weighted spectral angle match
WANG Linyu, LIU Yibo    
College of Information and Communication Engineering, Harbin Engineering University, Harbin 150001, China
Abstract: In order to overcome the problem that hyperspectral image background samples contain anomalous pixels and noise, a kernel RX anomaly detection algorithm based on weighted spectral angle match (SAM) was proposed. Firstly, k-means clustering was performed on the background pixels of the image to obtain the cluster centers, then the spectral angle cosine of the background pixels and the cluster centers was calculated, the pure background was selected as the new background. Each pixel in the new background will own its spectral angle information as the weight, which is given to every background pixel to construct weighted kernel RX anomaly detector to weaken the interference of the residual outliers and noise. To validate the effectiveness of the proposed algorithm, experiments were conducted on real hyperspectral data from AVIRIS and ROSIS-03 remote sensor. The results show that by comparison with the compared algorithms, the proposed algorithm has strong suppression ability against potential outliers and can improve the detection accuracy.
Key words: hyperspectral image    K-means clustering    weighted    kernel RX    spectral angle match    anomaly detection    spectral angle cosine    background optimization    

高光谱图像是一种具有“图谱合一”特性的遥感图像,具有很高的光谱分辨率,在目标检测与识别方面与多光谱相比有很大优势[1]。异常检测能在没有目标先验信息的情况下,将不符合背景统计特性的异常检测出来,在军事、农业、地质勘测等领域有很大的利用价值[2-3]

常用的高光谱异常检测算法是由Reed和Yu[4]提出的RX算法,是一种基于广义似然检验的恒虚警比率检测算法,RX算法仅仅利用了高光谱数据的低阶统计特性且假定背景信息符合高斯分布,认为检测窗内所有像元均为背景,忽略了潜在的异常像元对背景均值以及协方差矩阵估计的影响。在实际应用中高光谱图像常包含多种地物[5],背景复杂多变,所以对图像的检测效果不是很理想。为充分利用高光谱数据的非线性信息,Kown等[6]利用核函数将原始空间数据映射到高维特征空间,构造了核RX(kernel RX,KRX)算法,但KRX算法同样未对包含在背景中的异常抑制,对背景协方差估计不准确,检测结果中仍含有较多虚警点。文献[7]提出了一种核加权RX(weighted kernel RX,WKRX)算法,利用背景像元到质心的欧氏距离,自适应地赋予背景像元权值,从而削减异常信息在协方差矩阵中的影响,取得了良好的检测效果。但对于欧氏距离相同的两组数据,不一定具有很高的光谱相似度[8],以此作为加权准则难以保证异常像元与质心距离的最大化。文献[9]提出了阻塞式自适应高效计算异常探测器(blocked adaptive computationally efficient outlier nominators,BACON)算法,通过RXD迭代方式净化背景,并引入校正因子增强背景净化的稳健性,取得了不错的检测效果。文献[10]在BACON的基础上,提出了随机子集估计的异常检测(random-selection-based anomaly detector,RSAD)算法,随机选取一部分像元作为初始背景,降低异常在初始背景中的比例,同样采用迭代方式更新背景,最终获得稳健的背景,有效检测出异常目标。

对背景中潜在的异常进行抑制,有利于准确估计背景统计信息,提高检测精度。本文提出了一种基于光谱角匹配加权的异常检测算法,采用光谱角余弦(spectral angle cosine,SAC)作为光谱角匹配准则,使用双窗口模型,窗口大小由感兴趣的目标大小决定[10]。首先通过K-均值聚类得到不同类背景相对应的聚类中心,求取背景像元与聚类中心的SAC,选出较为纯净的背景作为新背景,然后新背景中的每个像元将自己的光谱角余弦信息作为权值结合KRX算法进行加权KRX异常检测,最后利用真实高光谱数据进行仿真实验,通过与RX、KRX、加权RX(weighted RX,WRX)和WKRX以及BACON进行对比,说明了所提算法的有效性。

1 核RX算法

经典的RX算法是由Reed等人于1990年提出的,其本质是求样本到总体均值向量的马氏距离。为了充分利用高光谱数据非线性信息,Kwon等人将原始空间数据通过非线性映射函数ϕ映射到高维特征空间中,利用核函数构造了基于核方法的KRX算法,将N个光谱向量 ${{\mathit{\boldsymbol{X}}}_b} = [{{\mathit{\boldsymbol{x}}}_1},{{\mathit{\boldsymbol{x}}}_2} \cdots ,{{\mathit{\boldsymbol{x}}}_n}]$ 映射到高维特征空间得 $\phi \left( {{{\mathit{\boldsymbol{X}}}_b}} \right) = \left[ {\phi \left( {{{\mathit{\boldsymbol{x}}}_1}} \right),\phi \left( {{{\mathit{\boldsymbol{x}}}_2}} \right), \cdots ,\phi \left( {{{\mathit{\boldsymbol{x}}}_n}} \right)} \right]$ 。KRX异常检测算子可表示为

${\rm KRX}\left( {\phi \left( {\mathit{\boldsymbol{r}}} \right)} \right) = {\left( {\phi \left( {\mathit{\boldsymbol{r}}} \right) - {{\mathit{\boldsymbol{\mu }}}_\phi }} \right)^{\rm{T}}}{{\hat{C}}_{b\phi }}^{ - 1}\left( {\phi \left( {\mathit{\boldsymbol{r}}} \right) - {{\mathit{\boldsymbol{\mu }}}_\phi }} \right)$

式中 ${{\mathit{\boldsymbol{\mu }}}_\phi }$ ${{\mathit{\boldsymbol{C}}}_{b\phi }}$ 分别为特征空间背景均值向量和协方差矩阵的估计,表达式为

$\begin{array}{c}{{\mathit{\boldsymbol{\mu }}}_\phi } = \displaystyle\frac{1}{N}\sum\limits_{i = 1}^N {\phi \left( {{{\mathit{\boldsymbol{x}}}_i}} \right)} \\{{{\hat{C}}}_{b\phi }} = \displaystyle\frac{1}{N}\sum\limits_{i = 1}^N {(\phi \left( {{{\mathit{\boldsymbol{x}}}_i}} \right) - {{\mathit{\boldsymbol{\mu }}}_\phi }){{(\phi \left( {{{\mathit{\boldsymbol{x}}}_i}} \right) - {{\mathit{\boldsymbol{\mu }}}_\phi })}^{\rm{T}}}} \end{array}$

为了便于计算,通常采用原始数据空间核函数来表示高维特征空间数据的内积:

$k({{\mathit{\boldsymbol{x}}}_i},{{\mathit{\boldsymbol{x}}}_j}) = < \phi ({{\mathit{\boldsymbol{x}}}_i}),\phi ({{\mathit{\boldsymbol{x}}}_j}) > {\rm{ = }}\phi {({{\mathit{\boldsymbol{x}}}_i})^{\rm{T}}} \cdot \phi ({{\mathit{\boldsymbol{x}}}_j})$

KRX算法能充分利用高光谱数据的非线性信息,本文加权选择在KRX基础上进行,在实验部分与WRX算法进行对比。根据构造背景方式不同,异常检测算法分为局部异常检测算法和全局异常检测算法。全局异常检测算法采用图像的所有像元作为背景,而局部异常检测算法基于滑动同心双层窗模型进行检测,外窗内的像元集合作为背景,内窗作为保护窗将待检测像元与背景分离。由于局部异常检测能够更有效地检测出图像中小异常目标[11],因此本文将使用双层窗模型进行加权KRX异常检测。

2 基于光谱角匹配加权的异常检测算法 2.1 K-均值聚类

K-均值聚类算法是模式识别中的经典算法,因其计算简便、能够动态聚类、并且聚类效果明显,被广泛应用于各个领域。高光谱图像中异常目标占有的比例非常低,通常不能形成单独的类,可以通过K-均值聚类将不同类的背景像元归属到相应的类中。K-均值聚类算法的基本步骤如下:

1) 在数据中选取K个向量作为初始聚类中心,记为 ${{\mathit{\boldsymbol{z}}}_1},{{\mathit{\boldsymbol{z}}}_2}, \cdots ,{{\mathit{\boldsymbol{z}}}_k}$

2) 计算数据中的每个向量 ${{\mathit{\boldsymbol{x}}}_i}$ 到每个聚类中心的最小距离,按最小距离准则将每个向量划分到相应的类中,即如果 ${d_{ij}} = \min {\left\| {{{\mathit{\boldsymbol{x}}}_i} - {{\mathit{\boldsymbol{z}}}_j}} \right\|_2}$ ,其中 $i = {\rm{1,2,}} \cdots ,N$ N为数据个数, $j = {\rm{1,2,}} \cdots ,k$ ,则 ${{\mathit{\boldsymbol{x}}}_i}$ 属于以 ${{\mathit{\boldsymbol{z}}}_j}$ 为中心的 ${{\mathit{\boldsymbol{w}}}_j}$ 类中。

3) 重新计算各类的中心:

${{\mathit{\boldsymbol{z}}}_j} = \frac{1}{n}\sum\limits_{{{\mathit{\boldsymbol{x}}}_i} \in {{\mathit{\boldsymbol{w}}}_j}} {{{\mathit{\boldsymbol{x}}}_i}} ,j = 1,2, \cdots ,k$

4) 如果类中心不再改变,则算法结束;否则转至步骤2)继续进行聚类。

通过聚类得到的聚类中心用于描述不同类的背景像元,从而进行后续的光谱角匹配,相比于仅仅使用检测窗内的背景均值作光谱角匹配,在处理检测窗内含有多类背景的情况时,将聚类中心作为参考对象进行光谱角匹配,得到的结果更精确。

2.2 光谱角匹配加权

光谱角匹配(spectral angle match, SAM)对于光谱向量的相似性描述效果比较理想,在高光谱分类、解混等领域运用广泛,相比于欧氏距离描述的光谱相似性,利用SAM描述更具有准确性[8],它不仅反应了2个光谱向量在数值上的差异,也能体现出光谱向量在光谱曲线、光谱形状上的差异。实际中通常通过计算SAC代替光谱角的计算[12],SAC的计算公式为

${S_{{\rm{SAC}}}}({\mathit{\boldsymbol{x}}},{\mathit{\boldsymbol{y}}}) = \frac{{\langle {\mathit{\boldsymbol{x}}},{\mathit{\boldsymbol{y}}}\rangle }}{{\sqrt {\langle {\mathit{\boldsymbol{x}}},{\mathit{\boldsymbol{x}}}\rangle } \sqrt {\langle {\mathit{\boldsymbol{y}}},{\mathit{\boldsymbol{y}}}\rangle } }}$

式中:xy为像元的光谱向量, $\langle \cdot , \cdot \rangle $ 为内积运算。光谱角是2个光谱向量的夹角,光谱角越小,计算出的光谱角余弦越接近1,两光谱向量就越相似。

通常情况下,我们认为检测窗口内像元就是背景像元,直接用其来计算协方差矩阵和均值,但当像元中存在异常目标时,计算出来的协方差矩阵和均值就不能准确地表征背景数据的分布,因此需要通过加权来削弱异常目标在检测窗口内所占的比重。

通过K-均值聚类方法得到相应的类中心 ${{\mathit{\boldsymbol{z}}}_j}$ ,然后计算每个背景向量 ${{\mathit{\boldsymbol{x}}}_i}$ 与各类中心的SAC记为 ${T_{{\rm{SAC}}}}_{_{ij}}$ ,取 $\max [{T_{{\rm{SAC}}}}_{_{ij}}]$ 作为 ${{\mathit{\boldsymbol{x}}}_i}$ 的光谱角余弦值,记为 ${S_{{\rm{SAC}}}}({{\mathit{\boldsymbol{x}}}_i},{{\mathit{\boldsymbol{z}}}_j})$ 并按升序排列。将阈值η设置为第M个相似度量值,因为异常目标面积比较小,所占像元数目较少,可以先假定M为像元总数目的5%[13],剔除 ${S_{{\rm{SAC}}}}({{\mathit{\boldsymbol{x}}}_i},{{\mathit{\boldsymbol{z}}}_j})$ 小于η的背景像元,大于等于η的所有像元构成较纯净的新背景,将新背景像元的 ${S_{{\rm{SAC}}}}({{\mathit{\boldsymbol{x}}}_i},{{\mathit{\boldsymbol{z}}}_j})$ 记为 ${s_i}\left( {i = {\rm{1,2,}} \cdots ,N} \right)$ ,将si的最大值和最小值记为smaxsmin,最大误差 $d = {s_{\max }} - {s_{\min }}$ ,将sismin相减差值与d的比值作为权值赋予每个背景像元,为了避免权值等于0导致加权后协方差矩阵不可逆,对每个背景像元的SAC加上一个极小值 $\varDelta $ ,这样背景像元权值可表示为

${w_i} = ({s_i} - {s_{\min }}{\rm{ + }}\varDelta )/\left( {d{\rm{ + }}\varDelta } \right)$ (1)

通过光谱角匹配加权,与背景越接近的像元,其对应的权值越大。将权值赋予背景中的每个像元,抑制背景中残留的异常目标对背景协方差矩阵和均值计算造成的干扰,从而更有效地表征出背景数据特征分布。

2.3 加权KRX异常检测算子

进行加权之后的协方差矩阵与均值向量分别表示为

$\begin{array}{c}{{{\hat{C}}}_{w\phi }} = \displaystyle\sum\limits_{i = 1}^N {{w_i}(\phi \left( {{{\mathit{\boldsymbol{x}}}_i}} \right) - {{\mathit{\boldsymbol{\mu }}}_{w\phi }}){{(\phi \left( {{{\mathit{\boldsymbol{x}}}_i}} \right) - {{\mathit{\boldsymbol{\mu }}}_{w\phi }})}^{\rm{T}}}} /\displaystyle\sum\limits_{i = 1}^N {{w_i}} \\{{\mathit{\boldsymbol{\mu }}}_{w\phi }} = \displaystyle\sum\limits_{i = 1}^N {{w_i}\phi \left( {{{\mathit{\boldsymbol{x}}}_i}} \right)} /\displaystyle\sum\limits_{i = 1}^N {{w_i}} \end{array}$

则加权KRX异常检测算子表示为

${W_{{\rm{KRX}}}}\left( {\phi \left( {\mathit{\boldsymbol{r}}} \right)} \right) = {\left( {\phi \left( {\mathit{\boldsymbol{r}}} \right) - {{\mathit{\boldsymbol{\mu }}}_{w\phi }}} \right)^{\rm{T}}}{{\hat{C}}_{w\phi }}^{{\rm{ - }}1}\left( {\phi \left( {\mathit{\boldsymbol{r}}} \right) - {{\mathit{\boldsymbol{\mu }}}_{w\phi }}} \right)$

根据文献[7]推导可得:

$\begin{split}& {W_{{\rm{KRX}}}}\left( {\phi \left( {\mathit{\boldsymbol{r}}} \right)} \right) = {\left( {\phi \left( {\mathit{\boldsymbol{r}}} \right) - {{\mathit{\boldsymbol{\mu }}}_{w\phi }}} \right)^{\rm{T}}}{{{\hat{C}}}_{w\phi }}^{{\rm{ - }}1}\left( {\phi \left( {\mathit{\boldsymbol{r}}} \right) - {{\mathit{\boldsymbol{\mu }}}_{w\phi }}} \right) = \\& W{\left( {\phi \left( {\mathit{\boldsymbol{r}}} \right) - {{\mathit{\boldsymbol{\mu }}}_{w\phi }}} \right)^{\rm{T}}}{{\mathit{\boldsymbol{X}}}_{w\phi }}{\mathit{\boldsymbol{K}}}_w^{ - 1}{\mathit{\boldsymbol{X}}}_{w\phi }^{\rm{T}}\left( {\phi \left( {\mathit{\boldsymbol{r}}} \right) - {{\mathit{\boldsymbol{\mu }}}_{w\phi }}} \right)\end{split}$ (2)

式中: $W{\rm{ = }}\displaystyle\sum\limits_{i = 1}^N {{w_i}} $ ${{\mathit{\boldsymbol{K}}}_w}$ 为特征空间加权中心化数据核矩阵, ${\phi _w}({{\mathit{\boldsymbol{x}}}_i}) \!=\!\! \sqrt {{w_i}} [\phi ({{\mathit{\boldsymbol{x}}}_i}) \!-\! {{\mathit{\boldsymbol{\mu }}}_{w\phi }}]$ ${{\mathit{\boldsymbol{X}}}_{w\phi }} \!=\! [{\phi _w}({{\mathit{\boldsymbol{x}}}_1}),{\phi _w}({{\mathit{\boldsymbol{x}}}_2}), \cdots, {\phi _w}({{\mathit{\boldsymbol{x}}}_N})]$ 。由于映射函数ϕ未知,无法求得 ${{\mathit{\boldsymbol{K}}}_w}$ ,因此需用输入空间未中心化的数据核矩阵 ${{\mathit{\boldsymbol{K}}}_0}[{({{\mathit{\boldsymbol{K}}}_0})_{ij}} = \langle \phi ({{\mathit{\boldsymbol{x}}}_i}),\phi ({{\mathit{\boldsymbol{x}}}_j})\rangle ]$ 来计算 ${{\mathit{\boldsymbol{K}}}_w}$ 。根据 ${\phi _w}({{\mathit{\boldsymbol{x}}}_i}) = \sqrt {{w_i}} [\phi ({{\mathit{\boldsymbol{x}}}_i}) - {{\mathit{\boldsymbol{\mu }}}_{w\phi }}]$ 及文献[14]中的核矩阵运算法则可得

$\begin{array}{l}{{\mathit{\boldsymbol{K}}}_w} \!\!=\! {\mathit{\boldsymbol{B}}}({{\mathit{\boldsymbol{K}}}_0} \!-\! {{\mathit{\boldsymbol{K}}}_0}{\mathit{\boldsymbol{w}}}{{\mathit{\boldsymbol{I}}}^{\rm{T}}}/W \!-\! {\mathit{\boldsymbol{I}}}{{\mathit{\boldsymbol{w}}}^{\rm{T}}}{{\mathit{\boldsymbol{K}}}_0}/W \!+\! {\rm{ }}{{\mathit{\boldsymbol{w}}}^{\rm{T}}}{{\mathit{\boldsymbol{K}}}_0}{\mathit{\boldsymbol{w}}}/{W^2}){\mathit{\boldsymbol{B}}} \!\!=\!\! {\mathit{\boldsymbol{B}}}{{\mathit{\boldsymbol{K}}}_x}{\mathit{\boldsymbol{B}}}\end{array}$

式中:I为元素都为1的列向量, ${\mathit{\boldsymbol{w}}} = {[{w_1},{w_2}, \cdots ,{w_N}]^{\rm{T}}}$ 是一个列向量, ${\mathit{\boldsymbol{B}}} = {\rm{diag}}\left( {\sqrt {{w_1}} ,\sqrt {{w_2}} , \cdots ,\sqrt {{w_N}} } \right)$ 为一个对角阵。

由式(2)可知,仍需要将在特征空间运算的 $\phi {({\mathit{\boldsymbol{r}}})^{\rm{T}}}{{\mathit{\boldsymbol{X}}}_{w\phi }}$ $\,{\mu _{w\phi }}^{\rm{T}}{{\mathit{\boldsymbol{X}}}_{w\phi }}$ 转化为低维空间的核运算:

$\begin{aligned}\phi {({\mathit{\boldsymbol{r}}})^{\rm{T}}}{{\mathit{\boldsymbol{X}}}_{w\phi }} = & \phi {({\mathit{\boldsymbol{r}}})^{\rm{T}}}\{ [\phi ({{\mathit{\boldsymbol{x}}}_1}),\phi ({{\mathit{\boldsymbol{x}}}_2}), \cdots ,\phi ({{\mathit{\boldsymbol{x}}}_N})] - {\mu _{w\phi }}\} {\mathit{\boldsymbol{B}}} = \\& [K({{\mathit{\boldsymbol{X}}}_b},{\mathit{\boldsymbol{r}}}) - K({{\mathit{\boldsymbol{X}}}_b},{\mathit{\boldsymbol{r}}}){\mathit{\boldsymbol{w}}}/W]{\mathit{\boldsymbol{B}}} = {\mathit{\boldsymbol{K}}}_r^{\rm{T}}\\\mu _{w\phi }^{\rm{T}}{{\mathit{\boldsymbol{X}}}_{w\phi }} \!\!=\!\! & \sum\limits_{i = 1}^N {{w_i}\phi \left( {{{\mathit{\boldsymbol{x}}}_i}} \right)} /W\{ [\phi ({{\mathit{\boldsymbol{x}}}_1}),\phi ({{\mathit{\boldsymbol{x}}}_2}), \cdots \phi ({{\mathit{\boldsymbol{x}}}_N})] - {\mu _{w\phi }}\} {\mathit{\boldsymbol{B}}}\!\!=\\& ({{\mathit{\boldsymbol{w}}}^{\rm{T}}}{{\mathit{\boldsymbol{K}}}_0}/W - {{\mathit{\boldsymbol{w}}}^{\rm{T}}}{{\mathit{\boldsymbol{K}}}_0}{\mathit{\boldsymbol{w}}}/{W^2}){\mathit{\boldsymbol{B}}} = {\mathit{\boldsymbol{K}}}_\mu ^{\rm{T}}\end{aligned}$

由此可得加权核RX异常检测算子可表示为

${W_{{\rm{KRX}}}}({\mathit{\boldsymbol{r}}}) = W({\mathit{\boldsymbol{K}}}_r^{\rm{T}} - {\mathit{\boldsymbol{K}}}_\mu ^{\rm{T}}){\mathit{\boldsymbol{K}}}_w^{ - 1}{({\mathit{\boldsymbol{K}}}_r^{\rm{T}} - {\mathit{\boldsymbol{K}}}_\mu ^{\rm{T}})^{\rm{T}}}$ (3)

本文提出的基于光谱角匹配的异常检测算法(weighted kernel RX anomaly detection algorithm based on spectral angle match,SWKRX)具体步骤如下:

1) 首先利用图像中全部像元构造均值和协方差矩阵,计算每个像元的马氏距离,然后选出马氏距离最小的前80%作为背景数据库,利用K-均值聚类将背景数据库中不同类的背景像元归属到相应的类中,得到每个类的聚类中心 ${{\mathit{\boldsymbol{z}}}_j}$

2) 通过计算 ${S_{{\rm{SAC}}}}({{\mathit{\boldsymbol{x}}}_i},{{\mathit{\boldsymbol{z}}}_j})$ 得到较纯净的背景,将新背景每个像元的si带入式(1)计算权值,削弱残留的异常目标的干扰。

3) 将权值wi带入式(3),得到像元判别值,构成一幅灰度图像,设定阈值,将灰度图像转换为二值图像。

3 实验结果与分析

为了验证本文算法的有效性,利用两幅真实高光谱图像进行仿真实验。实验在Intel Core i5-480M,主频2.67 GB CPU,4 GB内存,MATLAB 2014a环境下进行。

第一幅图像为50像元×50像元的高光谱图像,是美国圣地亚哥机场的一部分,图像覆盖了从可见光到近红外的连续光谱范围,共有224波段,波长范围为0.2~2.4 μm,去除受水汽影响和低信噪比的波段后,剩下的126个波段用于仿真实验,第60波段图像及地面目标分布图如图1所示,包含3架飞机。

图 1 第60波段图像及地面目标分布

窗口大小采用外窗13像元×13像元,内窗7像元×7像元,空间核函数选择常用的高斯径向基核函数,核参数σ选取10。为了验证本算法的检测效果,选择RX、KRX、采用本文加权方法的WRX以及WKRX和BACON作为对比算法,检测结果如图2所示。

图 2 不同算法在相同阈值下的检测结果

可以看出,RX算法直接假设背景符合高斯分布,并未对背景中潜在的异常目标做任何处理,检测结果中飞机形状不完整,虚警点多;WRX对于右上角的飞机检测出的目标像元多于RX;BACON经过背景净化能检测出比较完整的飞机形状;SWKRX、WKRX和KRX算法都充分挖掘了光谱间高阶统计特性,检测结果要好于RX和WRX,其中SWKRX检测出的目标像元多于KRX和WKRX,检测出来的飞机形状与地面目标分布图最接近,这充分说明了利用光谱角信息进行加权比欧氏距离更有效,更有利于削弱潜在异常目标对背景信息统计的影响。

接收机操作特性(receiver operating characteristic, ROC)曲线是常见的衡量检测特性的重要指标,用于描述检测概率Pd与虚警概率Pf之间的对应关系,能够为算法的检测性能提供定量分析。将检测到的真实目标像素个数与实际目标个数比值定义为检测率Pd,将检测到的虚警像素个数同整幅高光谱图像像素总数的比值定义为虚警率Pf

图3给出了6种算法的ROC曲线。可以看出,在相同虚警概率下,由于利用核函数的算法充分发掘了高光谱数据的非线性信息,所以检测概率明显高于RX算法,WKRX算法在虚警概率大于0.01后基本与KRX算法的检测概率相同。因为根据光谱角匹配加权相比于欧氏距离加权能更好地抑制目标信息对背景矩阵的干扰,所以SWKRX算法的ROC曲线变化最快,检测概率一直高于其他5种算法的检测概率。

图 3 6种算法的ROC曲线

表1给出了6种算法的ROC曲线下面积(area under the ROC curve, AUC)。通过比较AUC大小,可以明显看出SWKRX算法要优于其他5种算法,具有很好的检测性能。

表 1 6种算法的曲线下面积

第二幅图像选取了由ROSIS-03遥感器采集的意大利Pavia大学图像,尺寸为610像元×340像元,原始波段数为115个,波长范围为0.43~0.96 μm,去除受水汽影响及低信噪比波段后,剩下的102个波段用于仿真实验。截取了原图像的一部分用于实验,尺寸为115像元×115像元,第60波段原始图像及地面目标分布图如图4所示,含有7个待检测目标。

图 4 第60波段图像及地面目标分布

根据感兴趣的目标大小,将外窗大小设为11像元×11像元,内窗大小设为3像元×3像元,核函数选择为高斯径向基核函数,核函数的σ值经过多次实验对比,取值为40,6种算法检测结果如图5所示。因为RX算法仅仅利用了高光谱数据的低阶特性,忽略了波段间的相关性,所以RX算法只检测出5个目标并且有较多的虚警点;经过加权处理,WRX检测出了6个目标;SWKRX和WKRX抑制了在背景矩阵中的目标信息,检测结果都好于KRX,BACON检测出了全部目标且虚警个数少于WKRX。将图5(a)(b)(c)进行比较可以看出,SWKRX检测结果中虚警数目最少,SWKRX的效果要好于WKRX和BACON。

图 5 不同算法在相同阈值下的检测结果

图6给出了6种算法的ROC曲线。表2给出了6种算法的AUC。通过比较ROC曲线及AUC可以明显看出,在相同的虚警概率下,RX算法的检测概率最低,SWKRX算法检测概率最高,当虚警概率为0.008时检测概率已达到1。通过表2可以清楚地看出,SWKRX对应的AUC值也大于其他5种算法,充分体现出了SWKRX的有效性。

图 6 6种算法的ROC曲线
表 2 6种算法的曲线下面积

第三幅图像仍为美国圣地亚哥机场的一部分,大小为60像元×60像元,包含19个点目标,第60波段图像及地面目标分布图如图7所示。

图 7 第60波段图像及地面目标分布

窗口大小外窗为11像元×11像元,内窗3像元×3像元,高斯径向基核参数σ值为10。6种算法的检测结果如图8所示。由于这幅图像中异常目标较多且比较密集,RX算法仅仅利用低阶统计特性表征背景数据分布,在左下角以及右下角由于异常目标与背景像元的光谱较相似,RX无法有效检测出目标,具有较高的虚警率;WRX的检测出的目标像元个数多于RX;BACON和KRX检测出了18个目标,检测出的目标像元个数少于SWKRX和WKRX;WKRX以及本文提出的SWKRX都检测出了全部异常目标,相比较而言,SWKRX检测结果虚警点更少,体现了算法对于潜在异常目标的抑制作用。

图 8 不同算法在相同阈值下的检测结果

图9表3分别给出了6种算法的ROC曲线和AUC,通过比较可以看出,在相同虚警概率下,SWKRX算法的检测概率最高,对应的AUC也最大。

图 9 6种算法的ROC曲线

K-均值聚类方法中聚类中心个数K的确定与图像本身具有的类数相关,K的取值要大于等于真实的类数,常用的判别准则为最小误差平方和准则,但经实验验证,其只适用于低维且呈正团状分布数据自适应聚类数目的确定,而高光谱属于高维数据且不一定呈正团状分布,并不一定适合用最小误差平方和准则进行类数的确定。因此,本文中K值通过实验验证得到。如图10所示,给出了3次实验中K取值为2~10时对应的AUC,可见,当实验1中K取5,实验2中K取3,实验3中K取3时,SWKRX对应的AUC最大。

表 3 6种算法的曲线下面积
图 10 聚类个数K与AUC关系
4 结论

本文提出的基于光谱角匹配加权的高光谱图像异常检测算法。

1) 通过K-均值聚类得到高光谱图像中不同类背景对应的聚类中心,利用SAM得到较为纯净的背景;

2) 自适应地赋予每个新背景像元相应的权值,从而减小残留异常目标信息对背景特性估计的影响。

实验结果表明,利用光谱角匹配信息作为权值对背景像元进行加权处理,能够有效削弱异常目标对背景协方差矩阵和均值计算的影响。在相同的虚警概率下,与对比算法作比较,本文所提算法的检测概率最高,对应的接收机操作特性曲线下面积最大,具有良好的检测性能。但对高斯径向基核参数的选取是通过大量实验获得的,如何优化核参数的选择将有待进一步的研究。

参考文献
[1] 杨桄, 张俭峰, 赵波, 等. 基于PCA和KRX算法的高光谱异常检测[J]. 应用科技, 2014, 41(5): 11-13. (0)
[2] CHUDNOVSKY A, KOSTINSKI A, HERRMANN L, et al. Hyperspectral spaceborne imaging of dust-laden flows: anatomy of Saharan dust storm from the bodélé depression[J]. Remote sensing of environment, 2011, 115(4): 1013-1024. (0)
[3] TIWARI K C, ARORA M K, SINGH D. An assessment of independent component analysis for detection of military targets from hyperspectral images[J]. International journal of applied earth observations & sgeoinformation, 2011, 13(5): 730-740. (0)
[4] REED I S, YU X. Adaptive multiple-band CFAR detection of an optical pattern with unknown spectral distribution[J]. IEEE transactions on acoustics speech & signal processing, 1990, 38(10): 1760-1770. (0)
[5] SAFA K, ABDOLREZA S. Anomaly detection in hyperspectral images based on an adaptive support vector method[J]. IEEE geoscience and remote sensing letters, 2011, 8(4): 646-650. (0)
[6] KWON H, NASRABIDI N M. Kernel RX-algorithm: a nonlinear anomaly detector for hyperspectral imagery[J]. IEEE Transactions on geoscience and remote sensing, 2005, 43(2): 388-397. (0)
[7] 赵春晖, 李杰, 梅锋. 核加权RX高光谱图像异常检测算法[J]. 红外与毫米波学报, 2011, 29(5): 378-382. (0)
[8] 王玉磊, 赵春晖, 齐滨. 基于光谱相似度量的高光谱图像异常检测[J]. 吉林大学学报: 工学版, 2013, 43(增): 149-153. (0)
[9] BILLOR N, HADI A S, VELLEMAN P F. BACON: blocked adaptive computationally efficient outlier nominators[J]. Computational statistics & data analysis, 2000, 34(99): 279-298. (0)
[10] DU B, ZHANG L. Random-selection-based anomaly detector for hyperspectral imagery[J]. IEEE transactions on geoscience & remote sensing, 2011, 49(5): 1578-1589. (0)
[11] MOLERO J M, GARZON E M, GARCIA I, et al. Analysis and optimizations of global and local versions of the RX algorithm for anomaly detection in hyperspectral data[J]. IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing, 2013, 6(2): 801-814. (0)
[12] JAMES NORMAN SWEET. The spectral similarity scale and its application to the classification of hyperspectral remote sensing data[D]. New York: State university of New York, College of Environmental Science and Forestry, 2002: 40-54. (0)
[13] LOU C, ZHAO H. Local density-based anomaly detection in hyperspectral image[J]. Journal of applied remote sensing, 2015, 9(1): 095070. (0)
[14] JOHN SHAWE-TAYLOR, NELLO CRISTIANINI. Kernel methods for pattern analysis[M]. Cambridge: Cambridge University, 2004: 152-198. (0)