高光谱图像具有高分辨率、图谱合一等特点，可提供区别地物光谱微小差别的能力，在目标检测领域具有广泛应用^[1]。受到光照、大气等影响，很难获取准确的先验信息，因此不需要先验知识、具有较强实用性的异常检测得到了更多的关注。目前，对于异常目标还没有明确的定义，但普遍认为一个目标在一定程度上偏离了背景分布，即光谱信息明显不同于背景信息，则这个目标就是异常^[2]。

较为经典的异常检测算法是由Reed等提出的RX(reed-xiaoli)检测算法^[3]。在渐进意义下，其检测器原理为求待检测像元到周围背景像元均值之间的马氏距离。该算法是在高斯背景统计特性和空间白化的前提假设下进行的，而实际地物分布是复杂多变的，因此RX算法的检测性能受到限制。Kwon等提出KRX(kernel RX)算法，通过引入高斯径向基核函数，将高光谱数据映射到高维的特征空间进行处理，挖掘高光谱数据所隐含的非线性信息，显著提升检测精度^[4]。然而，高斯径向基核对于同类地物光谱强度的变化较为敏感，容易引起误判，且该算法仍然以RX基本假设作为前提，具有一定的局限性，且计算量较大。Li等提出基于协同表示的异常检测(collaborative representation based detector, CRD)算法，得到较好的检测性能，但算法易受到背景中存在的异常目标和噪声的干扰，从而引起误判^[5]。光谱角匹配算法可避免对背景模型的假设，经不断发展，现已应用于无监督的异常检测中^[6-7]。现有的基于光谱角累加的异常检测中，由于需要进行波段选择去噪的预处理过程，丢失了部分目标信息。同时，由于光谱角度量未充分利用像元间的非线性特性、对阈值的设置非常敏感，这些问题都会影响算法性能，因此该算法存在局限性。目前，异常检测算法更注重高光谱图像中丰富的光谱信息，而忽略了像元间的空间相关性。然而，为了更加准确、稳定的进行高光谱图像异常检测，综合考虑空间与光谱信息是很重要的。

随着GPU通用计算架构的不断发展，设计并行处理算法来提升效率成为研究热点^[8-11]。Antonio Plaza等利用GPU优化自动目标检测和分类算法(automatic target detection and classification algorithm, ATDCA)^[8]。Sanchez等提出了基于GPU并行计算架构的混合像元分解优化算法，并取得了较好的加速效果^[11]。通过利用GPU的并行计算优势，设计更优化的并行处理系统，并充分挖掘算法潜在的并行特性，可显著提升算法的执行效率。

为了去除对背景模型的假设，避免上述利用光谱角进行异常检测算法的弊端，本文提出一种基于核光谱角累加的异常检测算法(spatial-spectral kernel spectral angle mapping, SSKSAM)。该算法充分利用高光谱图像的空间信息，与光谱信息有效结合，通过计算待检测像元与其邻域内像元的核光谱角累加和，得到每个像元的异常程度。利用扩展形态学，对得到的异常程度结果进行修正，有效抑制噪声的同时，降低虚警率，从而得到最终的异常检测结果。考虑到实际应用中，处理的高光谱图像数据量巨大，因此利用GPU的并行计算能力，设计符合算法特性的并行系统架构。

1 基于核光谱角匹配的异常检测 1.1 核光谱角匹配

光谱角匹配(spectral angle mapping，SAM)是用来衡量两个光谱向量间差异程度的方法^[12]。如图 1所示，当两个光谱向量均为类内向量时，其对应的光谱角β相对较小。当两个光谱向量为类间向量时，其对应的光谱角α相对较大。设两个光谱向量为x=[x₁ x₂…x_L]与y=[y₁ y₂…y_L]，L为波段数。则光谱角距离公式可表示为

$ \begin{array}{l} {S_A}\left( {\boldsymbol{x},\boldsymbol{y}} \right) = \arccos \left( {\frac{{\sum\limits_{i = 1}^L {{x_i}{y_i}} }}{{\sqrt {\sum\limits_{i = 1}^L {x_i^2} } \sqrt {\sum\limits_{i = 1}^L {y_i^2} } }}} \right) = \\ \;\;\;\;\;\;\;\arccos \left( {\frac{{\left\langle {\boldsymbol{x},\boldsymbol{y}} \right\rangle }}{{\left\| \boldsymbol{x} \right\| \cdot \left\| \boldsymbol{y} \right\|}}} \right) \end{array} $

(1)

由式(1) 可知光谱角仅与向量的方向有关，与向量本身的长度无关，故光谱角能够有效地度量两个光谱形状上的相似性，降低因太阳入射角、地势、观测角度等因素的影响。此外，光谱角匹配具有良好的抵抗乘性干扰的能力，因此不受光照强度变化的影响。由于光谱角匹配未充分考虑波段间的非线性信息及潜在的空间信息，因此直接用于异常检测得到的检测效果并不理想。

为了避免光谱角匹配度量算法的不足，本文采用光谱角匹配核函数形式进行异常检测^[13]。通过非线性映射函数Φ可将原始空间数据r映射到高维空间中：

$ \boldsymbol{\mathit{\Phi}} \left( \boldsymbol{r} \right) = [\boldsymbol{\mathit{\Phi}} ({r_1})\;\;\boldsymbol{\mathit{\Phi}} ({r_2}) \ldots \;\;\boldsymbol{\mathit{\Phi}} ({r_L})] $

(2)

相应地得到特征空间中光谱角距离表达式为

$ {S_A}\left( {\boldsymbol{\mathit{\Phi}} \left( \boldsymbol{x} \right),\boldsymbol{\mathit{\Phi}} \left( \boldsymbol{y} \right)} \right) = \arccos \frac{{ < \boldsymbol{\mathit{\Phi}} \left( \boldsymbol{x} \right),\boldsymbol{\mathit{\Phi}} \left( \boldsymbol{y} \right) > }}{{\left\| {\boldsymbol{\mathit{\Phi}} \left( \boldsymbol{x} \right)} \right\| \cdot \left\| {\boldsymbol{\mathit{\Phi}} \left( \boldsymbol{y} \right)} \right\|}} $

(3)

由于特征空间的维数非常高，且非线性映射函数Φ是未知的，若直接在特征空间中计算将非常困难甚至不能实现。而核函数可将特征空间的点积运算转化为原始空间数据的核函数实现，因此利用核函数隐式映射可有效避免上述问题，核函数的一般形式为

$ k\left( {\boldsymbol{x},\boldsymbol{y}} \right) = \left( {\boldsymbol{\mathit{\Phi}} \left( \boldsymbol{x} \right),\boldsymbol{\mathit{\Phi}} \left( \boldsymbol{y} \right)} \right) $

(4)

高斯径向基核(radial basis function, RBF)是较为常用的正定核函数，表示为

$ k\left( {\boldsymbol{x},\boldsymbol{y}} \right) = \exp( - {\left\| {\boldsymbol{x} - \boldsymbol{y}} \right\|^2}/c) $

(5)

式中：c为常数，表示高斯径向基核的宽度。其原理是计算向量间的欧式距离，由式(5) 可知RBF核只与向量间差值大小有关，与其绝对位置无关，因而具有平移不变性。与光谱角距离结合可有效弥补RBF核对光照强度变化等敏感的缺陷，且当x=y时，k(x, y)=1，则式(3) 中的分母为1，可有效降低计算复杂度。将式(4) 和式(5) 代入式(3) 可得核光谱角距离(kernel spectral angle mapping, KSAM):

$ \begin{array}{c} {\rm{KSAM}}\left( {\boldsymbol{x},\boldsymbol{y}} \right) = {S_A}\left( {\boldsymbol{\mathit{\Phi}} \left( \boldsymbol{x} \right),\boldsymbol{\mathit{\Phi}} \left( \boldsymbol{y} \right)} \right) = \\ \arccos \left( {\frac{{ < \boldsymbol{\mathit{\Phi}} \left( \boldsymbol{x} \right),\boldsymbol{\mathit{\Phi}} \left( \boldsymbol{y} \right) > }}{{\left\| {\boldsymbol{\mathit{\Phi}} \left( \boldsymbol{x} \right)} \right\| \cdot \left\| {\boldsymbol{\mathit{\Phi}} \left( \boldsymbol{y} \right)} \right\|}}} \right) = \\ \arccos \left( {\frac{{k\left( {\boldsymbol{x},\boldsymbol{y}} \right)}}{{\left\| {k\left( {\boldsymbol{x},\boldsymbol{x}} \right)} \right\| \cdot \left\| {k\left( {\boldsymbol{y},\boldsymbol{y}} \right)} \right\|}}} \right) = \\ \arccos \cdot \\ \left( {\frac{{\exp( - {{\left\| {\boldsymbol{x} - \boldsymbol{y}} \right\|}^2}/c)}}{{\left\| {\exp( - {{\left\| {\boldsymbol{x} - \boldsymbol{x}} \right\|}^2}/c)} \right\| \cdot \left\| {\exp( - {{\left\| {\boldsymbol{y} - \boldsymbol{y}} \right\|}^2}/c)} \right\|}}} \right) = \\ \arccos \left( {\frac{{\exp( - {{\left\| {\boldsymbol{x} - \boldsymbol{y}} \right\|}^2}/c)}}{{\left\| {\exp\left( 0 \right)} \right\| \cdot \left\| {\exp\left( 0 \right)} \right\|}}} \right) = \\ \arccos (\exp( - {\left\| {\boldsymbol{x} - \boldsymbol{y}} \right\|^2}/c)) \end{array} $

(6)

式中：x为待检测像元向量，y为待检测像元x的邻域像元向量。在异常检测中，通常通过对比待检测像元与周围像元之间的差异程度，来判定待检测像元的异常程度。KSAM通过计算两个光谱向量之间的核光谱角距离评价光谱向量之间的相似程度，因此本文通过计算待检测像元与周围邻域内像元的核光谱角匹配累加结果衡量待检测像元的异常程度，核光谱角累加和表达式为

$ \begin{array}{c} R = \sum\limits_{i = 1}^n {{\rm{KSAM}}(\boldsymbol{x},{\boldsymbol{y}_i})} = \\ \sum\limits_{i = 1}^n {\arccos \left( {\exp\left( { - {{\left\| {\boldsymbol{x} - {\boldsymbol{y}_i}} \right\|}^2}/c} \right)} \right)} \end{array} $

(7)

式中：n为待检测像元邻域内像元个数。其核光谱角累加计算的邻域像元区域示意图如图 2所示。因待检测像元与自身的光谱角距离为0，所以在计算过程中无需剔除待检测像元。

1.2 扩展形态学

数学形态学是一种分析像元间非线性空间信息的处理技术^[14]。在形态学处理灰度图像时，将像元的数值大小作为排序标准进行最小或最大的灰度值计算，若将形态学用于高光谱图像中，因其像元为多维向量，所以不能简单地进行大小比较。Antonio Plaza等将形态学引入到了高光谱图像处理中，通过改进的扩展形态学操作有效地解决了上述不能简单直观比较像元大小的问题^[15]。在扩展形态学中，为了确定高光谱图像中多维光谱向量间的排序关系，引入一个度量算子，该度量算子是将高光谱图像中每个光谱向量与其结构元素内其他光谱向量的距离之和作为排序准则，表达式为

$ C\left( {\boldsymbol{F}\left( {x,y} \right)} \right) = \sum\limits_{\left( {s,t} \right) \in s} {{\rm{dist}}\left( {\boldsymbol{F}\left( {x,y} \right),\boldsymbol{F}\left( {s,t} \right)} \right)} $

(8)

式中：F(x, y)为N维光谱向量，N为波段数。F(s, t)为结构元素内其他像元光谱向量。dist()为N维向量的线性距离测度，用来表征光谱向量间的相似性。S为结构元素，其大小可根据异常目标的大小来设置。根据上述构造的度量因子，定义扩展形态学的膨胀和腐蚀操作：

$ \begin{array}{l} d\left( {x,y} \right) = \left( {\boldsymbol{F\varTheta S}} \right)\left( {x,y} \right) = \\ \arg {\max _{\left( {s,t} \right) \in s}}\left\{ {C\left( {\boldsymbol{F}\left( {x - s,y - t} \right),S} \right)} \right\} \end{array} $

(9)

$ \begin{array}{l} e\left( {x,y} \right) = \left( {\boldsymbol{F} \oplus S} \right)\left( {x,y} \right) = \\ \arg {\min _{\left( {s,t} \right) \in S}}\left\{ {C\left( {\boldsymbol{F}\left( {x + s,y + t} \right),S} \right)} \right\} \end{array} $

(10)

式中：argmax和argmin表示为度量因子达到最大或最小的像元向量。由地物分布的空间连续性可知，不同地物对应的端元点不会是相邻的像元，即端元点邻近像元通常为混合像元点。高光谱数据在特征空间中为单形体结构，混合像元点位于单形体内部，端元点位于顶点处或顶点附近，所以度量因子最小的必然为混合像元点，度量因子最大的更可能为端元点。通过上述分析可知，扩展形态学的膨胀结果将得到在结构元素内与背景信息差异较大的像元，腐蚀结果将得到在结构元素内与背景信息相似的像元，其示意图如图 3所示。

由于扩展形态学的腐蚀操作将得到在结构元素内部与背景相似的像元信息，因此可利用腐蚀操作对异常结果进行修正，以此来去除背景噪声的干扰，降低虚警率。因上节内容中，已讨论将核光谱角累加和作为衡量像元异常程度的标准，因此在进行扩展形态学腐蚀操作时，可利用核光谱角累加和作为扩展形态学的度量因子，在此基础之上进行扩展形态学的腐蚀操作，最后将由核光谱角累加和得到的像元异常程度与腐蚀结果进行差值计算，对异常结果进行修正，得到最终的异常检测结果。其表达式为：

$ \boldsymbol{r}\left( {x,y} \right) = \left\| {\boldsymbol{R} - \boldsymbol{e}\left( {x,y} \right)} \right\| $

(11)

式中：r(x, y)为异常检测最终结果像元，R为通过核光谱角累加计算的像元，e(x, y)为扩展形态学腐蚀结果像元。

图 4所示为SS-KSAM算法设计流程图。

2 SS-KSAM算法并行优化模型

实际应用中，需处理的高光谱数据量巨大的，而异常检测算法的执行效率会随着高光谱图像的增大而降低。因此，为了更符合实际需求，本节将所提出的SS-KSAM算法进行并行优化，设计适合算法的并行系统架构。

近些年，随着GPU的不断发展，其强大的通用计算能力为提升异常检测算法效率提供了良好的平台。随着NVIDIA公司推出统一设备架构CUDA模型，GPU在高光谱图像处理领域得到了更为广泛的应用。GPU由多个流多处理器(stream multiprocessors，SM)阵列组成，每个SM由多个标量处理器核心(streaming processor，SP)及少量其他计算单元构成。CUDA编程模型为CPU+GPU异构模式。CPU端为主机端，主要负责整体流程的控制及串行部分的计算。GPU端为设备端，其负责的并行部分计算被封装成为Kernel，根据并行执行的任务不同将Kernel配置成不同大小的网格，网格中包含多个线程块，而每个线程块又包含多个线程，通过索引将每个计算任务分配给不同的线程去完成，实现计算任务的并行化处理过程。

通过封装核函数的设计，可完成某一并行任务。在执行核函数任务时，通过线程资源的合理分配，并进行二维线程索引，将每个线程明确定位，从而实现线程的统一调度。图 5所示为Kernel_Dis核函数并行实现待检测像元与邻域像元二范数平方‖x-y‖²运算的设计理念示意图。通过分配不同的线程块，把多个经窗口滑动形成的像元组实现并行分块处理，再通过块内多个线程的系统调度设计，将每个具体任务划分给各独立线程去并行执行，从而实现用一组运算执行时间完成整个运算过程的高效目的。其中，在累加过程中可利用归约思想进行计算，图 6所示为归约计算流程示意图，假定共有N个元素参与归约计算，共需进行lb N轮运算。有效降低计算复杂度，进一步提升运算效率。本文设计多个不同的封装函数用于对SS-KSAM算法进行并行优化操作，充分挖掘算法潜在的并行特性，将算法并行度达到更优化。因GPU端与CPU端通过PCI-E总线进行通信，其传输带宽相对较慢，因此在并行设计时应该尽量避免不必要的数据通信，以此减少传输带宽对并行效果的影响，同时也可起到增大并行与串行部分运算比例的作用。

并行优化处理具体流程如下：

1) 构造扩展图像，因需要邻域窗口滑动实现像元异常程度的计算，所以需要将图像扩展来完成边缘像元的计算过程。扩展宽度由窗口大小来决定。

2) 将扩展后的高光谱图像通过CudaMemcpy()函数传送给设备端的全局存储器，并在设备端分配相应内存空间以供并行处理。

3) 调用设备端函数Kernel_Dis，并合理分配网格中的线程块和线程，设计全面系统的线程定位方式，将每个像元计算任务分配给某一线程块中的确定线程，从而完成‖x-y‖²的计算过程。

4) 计算待检测像元与其邻域像元的核光谱角距离，通过Kernel_Ang封装核函数，将核光谱角距离计算分配给多个线程进行独立并行实现，并将结果累加得到核光谱角累加和，表达式如式(7)。

5) 为便于图像边缘像元的腐蚀处理，将得到的异常程度结果通过设计的Kernel_Pad封装核函数进行多线程扩展，扩展宽度由结构元素大小确定。

6) 根据扩展形态学的定义，将上述像元异常程度作为度量因子进行腐蚀操作。利用GPU端核函数Kernel_Min，寻找待检测像元邻域内的极值，得到腐蚀结果，腐蚀操作见式(10)。

7) 通过上述步骤，得到了每个像元的异常程度结果及腐蚀结果，根据式(11)，调用Kernel函数Kernel_Com，通过线程定位设计将每个像元的差值运算过程分配给多个线程并行执行。

8) 最后，通过CudaMemcpy()函数将最终异常检测结果传送给主机端，完成异常检测过程。

3 实验结果与分析

为了验证本文所提的SS-KSAM算法的有效性及高效性，分别利用真实数据与模拟数据进行实验仿真操作。实验仿真环境：CPU处理器为Intel Xeon E7-4820，主频为2.0 GHZ，内存为DDR3 96 GB。GPU处理器为NVIDIA Tesla K40 m。实验平台为Microsoft Visual Studio 2010集成开发坏境。

3.1 算法有效性

采用一幅真实高光谱数据来验证本文所提算法的有效性，该高光谱图像由PaciaC(意大利北部的帕维亚中心)截取的像素大小为115×115的数据，背景主要由桥、水波、影子组成，去除水吸收波段及信噪比较低的波段后，剩余102个波段用于本次实验。图 7所示为此高光谱数据第20波段的灰度图像及真实地物分布的二值图像。

利用经典RX算法、KRX算法、CRD算法、SAM累加算法(简称SAM)和KSAM算法(该算法为本文算法的第一步，即无空间-光谱约束情况下直接利用KSAM累加)作为对比算法。通过交叉验证，RX、KRX及CRD算法的内窗设为3×3像素，外窗设为11×11像素，其他3种算法的累加邻域设为11×11像素。KRX及两种核光谱角距离算法的高斯核参数为10，SS-KSAM算法的形态学结构元素大小为3×3像素。由于CRD对参数λ不敏感，因此直接将其设为1^[5]。图 8为各算法异常检测结果灰度图，图 9为各算法异常检测结果灰度峰度图。

ROC曲线描述检测概率和虚警概率间的变化关系，可用来对异常检测算法的性能进行定量分析。图 10为上述6种算法的ROC曲线图。

由异常检测结果图及ROC曲线可以看出，因传统RX算法及SAM算法均未充分利用高光谱图像波段间的非线性信息，导致其检测性能低于其他三种算法，且由于RX算法需要背景模型假设条件的局限，导致其在虚警概率为0.03时，还未准确检测到全部异常点。KRX算法通过引入核函数，充分利用了波段间的非线性信息，因而具有较强的非线性检测性能。根据ROC曲线可知，由于SS-KSAM算法利用数据空谱信息且具有较强的光谱形状分辨能力，使其在虚警概率高于0.003时，检测概率超过KRX算法，且在虚警概率为0.006时即可检测出全部的异常点。但由于核光谱角匹配对像元光谱形状具有较高的敏感性，使得在虚警概率低于0.003时，SS-KSAM算法的检测概率稍微低于KRX算法，但KRX算法仍需在背景假设模型的局限下进行，且算法的计算复杂度要明显高于SS-KSAM算法，同时，通过图 8与图 9中(b)和(f)的对比可看出，KRX算法对局部噪声的抑制效果要低于SS-KSAM算法。因此，SS-KSAM算法检测性能优于KRX算法。因CRD算法易受噪声影响，所以检测精度稍低于SS-KSAM算法。与KSAM算法相比，SS-KSAM算法检测概率有了一定提高，且在虚警概率高于0.002 8时，性能优势更为明显，并且由灰度图对比可看出，SS-KSAM算法具有更好的抑制噪声的能力。综上所述，SS-KSAM算法的检测性能优于其他五种算法。

3.2 算法有效性

采用一组模拟数据进行算法的高效性验证实验。此模拟数据由AVIRIS高光谱遥感中的多种目标进行目标插入得到，共有4×4=16个目标，合成数据大小为128×128像素，共128个波段。每行目标由同种地物构成。其中第一列目标由纯像元构成，第二列目标包含75%目标像元光谱特性和25%背景光谱特性，第三列目标包含50%目标像元光谱特性和50%背景光谱特性，第四列目标包含25%目标像元光谱特性和75%背景光谱特性。图 11为模拟数据第20波段灰度图像及其目标分布二值图像。

为了充分展现随着高光谱数据量的增大，算法的并行实现效果，在保证目标的密度不变的情况下，将上述模拟数据扩展成256×256像素和512×512像素，分别利用这三组空间大小不同的模拟数据进行实验。表 1为三组数据下，SS-KSAM算法串行及并行实现检测时间加速比情况。

为了去除计算机本身的脉冲误差，表 1中的算法检测时间均为100次实验结果的均值。由表 1可看出，在数据大小相同的条件下，由GPU并行优化的算法检测时间明显快于串行算法。模拟数据为128×128大小时，与CPU串行算法相比，加速比为4.294倍，当模拟数据为512×512大小时，加速比可达到14.776倍，可见随着数据量的增大，GPU的加速效果也更加凸显。因此利用GPU/CUDA架构进行算法并行优化处理更加符合实际应用海量数据高效处理需求。上述实验证明，通过并行优化处理的异常检测算法的执行效率得到显著提高，并随着高光谱图像的增大，其加速效果会愈加明显。

4 结论

1) 随着高光谱图像空间分辨率的不断增大，亚像元、混合像元及背景噪声等因素对检测性能的影响不可忽视，因此为适应当前高光谱图像的特性，充分利用空间信息是解决问题的重要方向之一。

2) 充分利用GPU并行特性，设计并行系统架构，提升了算法执行效率，且随数据量的增大，加速效果愈加凸显。

3) 通过理论分析与实验结果证明，该算法具有较强的实用性，在保证较高的检测精度的同时，实现算法的高效检测过程。

下一步将对算法参数的自适应选取进行深入研究。