遥感图像异常检测是一类无需先验信息监督的广义目标检测技术。它不仅包含了感兴趣目标，同时可以发掘可能感兴趣的区域，与环境监测、国防安全、矿藏勘测、公共营救等领域密切相关，具有重要的研究意义^[1-2]。

随着高分辨率传感设备的投入使用，地表地物构成与纹理结构特征的精细描述成为现实。同时，背景地物的高精细描述使得丰富边缘纹理信息等高频特性与异常目标特征形成较强的相关性，仅仅依据线性的分割平面无法准确区分背景与异常。寻找能将复杂背景特征与异常目标成功区分的非线性分割平面是目前异常检测算法研究的热点。

目前，面向复杂背景的异常检测比较主流的算法有基于复杂背景模型^[3]、机器学习及非线性空间变换^[4]等。然而，由于涉及背景特征的分类、学习以及建模，通常具有较高的计算复杂度且消耗大量的计算时间，较难满足遥感图像解译对实时性的要求。

针对以上问题，无需背景建模的新型算法不断涌现。2005年Gaucel等人^[5]提出基于空域相关性滤波器的异常检测算法(WSCF)。2012年Cohen等人^[6]提出的背景纹理提前抑制及后处理方法。2014年Romano等人^[7]提出利用形态学分析，通过增强异常目标特征或者抑制背景特征有效降低异常检测的虚警率。2015年Yuan等人^[8]提出一种高阶2维交叉滤波器(2DCAD)，利用异常目标显著边缘特征作为检测依据，实现快速超光谱异常检测的算法。该算法无需背景建模，在实时性能上体现出较强的优势。

然而对于存在复杂背景特征的遥感图像，相比于高亮度的背景地物，小尺寸异常目标在空域上的显著性明显下降。仅从单波段空域特征入手，小尺寸异常目标会由于缺乏显著空域特征而被漏检。另外，背景地物的高频空域特征与目标空域特征高度相关，致使线性滤波无法区分高频信息的来源，从而导致检测结果虚警率偏高。

生物视觉系统通过神经细胞的逐级滤波处理，可以从海量视觉数据中提取有用信息，实现复杂环境下微小目标快速、准确的检测与识别^[9]，其高效的自适应特性以及超敏锐信息感知能力均优于现有的计算机视觉系统。其中，生物视觉通过非线性滤波综合时、空域信息，提取高维特征信息实现目标检测与识别，不仅有效抑制背景与目标在时空域上的冗余，同时增强感兴趣目标的特征信息，具有较强的抗噪能力。

基于此，受到生物视觉非线性滤波机理的启发，提出构建联合光谱域与空域信息的相关型仿生非线性滤波器。该滤波器利用波段之间光谱变化特征的差异增强小尺寸异常目标特征，同时抑制复杂背景特征。该算法弥补了线性滤波器仅从空域分析无法区分背景特征和目标特征的缺点，实现快速有效的异常检测。仿生非线性滤波器计算复杂度较低，无需背景建模，具有较好的实时性与普适性。仿真实验将从检测算法的虚警率及计算时间等方面验证所提算法的优越性。

1956年Richard设计了一种相关型的视觉信息提取模型，并命名为初级运动检测器模型(EMD)。研究证实，该模型是目前描述生物视觉系统初级信息提取机理最成功的数学模型^[10]。EMD模型模拟生物视觉系统联合时、空域信息提取高维特征信息的生物学机理，在冗余信息压缩、有用信息增强方面具有计算原理简单，抗噪性能强等特点。

EMD模型通过空域特征在时域上的变化情况，判断是否存在随着时间变化而变化的空域特征。不随时间变化的空域特征(即背景特征)将被抑制；随时间变化而变化的空域特征(目标的边缘特征)则被增强。

以水平方向为例，将不同时刻不同空间位置上的亮度信息相乘作为EMD水平方向的输入，然后通过差分计算获得初级运动检测器的输出，即

$\begin{gathered} {E_h}=\left({i, j, t} \right)={I_t}\left({i, j} \right) \times {I_{t + \tau }}\left({i + \Delta s, j} \right)-\hfill \\ {I_t}\left({i + \Delta s, j} \right) \times {I_{t + \tau }}\left({i, j} \right) \hfill \\ \end{gathered} $

(1)

同理, 垂直方向输出E_v(i, j, t)为

$\begin{gathered} {E_v}=\left({i, j, t} \right)={I_t}\left({i, j} \right) \times {I_{t + \tau }}\left({i, j + \Delta v} \right)-\hfill \\ {I_t}\left({i, j + \Delta v} \right) \times {I_{t + \tau }}\left({i, j} \right) \hfill \\ \end{gathered} $

(2)

式中, I_t(i, j)和I_t+τ(i+Δs, j)分别为第t时刻图像上第(i, j)像素的亮度值和第t+τ时刻第(i+Δs, j)像素位置上的亮度值。Δs、Δv分别为水平、垂直方向上像素之间的空间距离，τ为时间间隔。

EMD模型通过特定空间区域Δs内，特定时间间隔内空域特征的相关性检测该区域中是否存在因目标运动而产生的随时间变化而变化的高维特征。若不同时刻上的空域特征不相关，说明在该时间间隔内有运动目标通过，其边缘以及表面纹理特征在局部空域内与背景不具备相关性。

EMD模型获得的高维特征综合了时域、空域双重信息，大大消除了背景及目标在空域、时域上的相关性冗余。另外，高维特征不依赖于空域信息或时域信息，即便目标的空域特征不显著，也会因其在时域上的变化而被EMD感知。这类不同属性信息交互作用的检测模式较好地克服了依赖单一特征(时域特征或空域特征)检测运动目标的片面性，有利于获得更为可靠和稳定的检测结果。

在遥感图像中，地物属性的差异反映在谱变化特征的差异性上，因此，判断地物属性通常需要综合多个甚至所有波段的光谱数据。异常目标是一类具有特殊光谱特性的地物，通常因其本身地物属性的稀缺性及特殊性，在与周围背景地物交接处容易形成较明显的光谱亮度跳变，在空域上则表现为目标的边缘特征。因此，从空域上的光谱亮度差异分析，采用线性滤波等方法可以快速区分异常目标和背景。

然而，现实情况下比较容易出现的两类问题是：背景特征复杂，异常目标边缘特征的高频属性将降低，利用空域滤波提取高频信息的方法不利于异常目标检测；另外，当在背景具有较高亮度光谱值的背景地物时，异常目标的光谱异常性则不显著。此时，单独分析单个波段光谱变化特征不能准确检测异常。虽然，单波段空域特征差异不显著，但比较多个波段光谱变化特性仍存在较强的差异性。

基于以上分析，如何联合空域与光谱域信息挖掘异常目标在光谱波形上的差异是实现复杂背景下异常目标检测的关键。

基于以上分析，基于异常目标光谱特性在波谱变化特征的异常性，利用仿生非线性滤波有效抑制复杂背景特征对异常检测准确性的干扰，同时联合空域信息和光谱信息挖掘异常目标，算法流程如图 1所示。该算法主要分为仿生非线性滤波、二值化处理及基于双窗模型的异常检测3部分。

1) 仿生非线性滤波。基于生物视觉非线性滤波机理通过空域与光谱域的交互分析，将局部区域内反复出现的背景纹理抑制，同时增强异常目标的边缘特征。

2) 二值化处理。带通滤波获得图像的中频成分。通过中频成分将高频与低频数据分开，高于中频的成分用“1”表示，相反则用“0”表示。

3) 基于双窗口模式的异常检测。采用双窗口模式，对待检测像元进行局部范围内的异常显著性统计作为最后的检测结果。保护带的存在有效阻断了目标与背景的混淆，从而避免了较大尺寸目标中心区域的异常性丢失问题。

考虑到遥感图像多波段提供的地物波谱特性将有助于提取场景中的边缘纹理特征。按照空域信息与光谱信息交叉相关的思想，选择空间相邻两个像元的光谱强度值作为特征检测器模型的输入，采用对称交叉相乘的方式计算空域特征相关性如图 2所示。

图 2中白色圆点和灰色圆点分别标记了对应相乘的像元，I_k1(i, j)和I_k2(i, j)分别为k₁波段与k₂波段在相同空间位置(i, j)上的光谱强度值，I_k1(i+Δø, j)为k₁波段在水平方向上与第(i, j)个像元间隔Δø距离的像元光谱强度值，I_k2(i+Δø, j)同理.对于第(i, j)个像元，输出R_h(i, j)中的高亮度像素代表局部区域内具有较强的特殊光谱特征。

$\begin{gathered} {R_h}=\left({i, j, } \right)={I_{k1}}\left({i, j} \right) \times {I_{k2}}\left({i + \Delta \phi, j} \right)-\hfill \\ {I_{k1}}\left({i + \Delta \phi, j} \right){I_{k2}}\left({i, j} \right) \hfill \\ \end{gathered} $

(3)

$\begin{gathered} {R_v}=\left({i, j, } \right)={I_{k1}}\left({i, j} \right) \times {I_{k2}}\left({i, j + \Delta \phi } \right)-\hfill \\ {I_{k1}}\left({i, j + \Delta \phi } \right){I_{k2}}\left({i, j} \right) \hfill \\ \end{gathered} $

(4)

式中，(R_h(i, j), R_v(i, j))为第(i, j)个像元在局部区域内的光谱波形特征向量。参数Δø为空域尺度参数，较小的Δø表示较精细的观测尺度；相反较大的Δø表示较粗糙的观测尺度。通过欧氏距离，计算第(i, j)相关的光谱特征强度值

${M_{{k_1}, {k_{2, \Delta \phi }}}}\left({i, j} \right)=\sqrt {{R_h}{{\left({i, j} \right)}^2} + {R_v}{{\left({i, j} \right)}^2}} $

(5)

M_{k1, k2, Δø}(i, j)值大小反映了第(i, j)像元周围邻近区域内光谱强度变化的剧烈程度，即该空间位置上是否存在显著的光谱特征差异。

考虑到不同波段组合下的光谱特征不同，为了统计整个波谱范围内的光谱变化特性，叠加不同波段组合获得的特征信息M_{k1, k2, Δø}(i, j)获得第(i, j)像元的光谱特征值

${M_{\Delta \phi }}\left({i, j} \right)=\sum\limits_{{k_1}, {k_2}} {{M_{{k_1}, {k_2}, \Delta \phi }}\left({i, j} \right)} $

(6)

式中，k₁、k₂为波段数，其取值范围为[1, K]，K为波段总数。图 3选取Landset-7卫星获得的多光谱遥感数据，采用仿生非线性滤波提取高维特征信息，提取结果如图 3(b)所示。其中波段组合分别为(1, 2)、(3, 4)、(5, 6)，叠加不同波段组合下获得的提取结果可得到如图 3(b)所示的边缘纹理特征信息。其中空域尺度为Δø=1，波段组合选择为相邻波段组合。

异常目标突变的边缘信息是较为可靠的检测依据。2007年Urhan^[11]提出采用图像的中频成分作为分割阈值，将图像中的高频信息与低频信息区分开，分别采用“1”和“0”的二值化符号标记异常与背景。该中频信息为像元周围邻域光频谱强度的均值，可以较好地克服自定义阈值容易导致分类误差的问题。

对于高维特征图像数据M，中频成分$\hat{\boldsymbol{M}}$采用带通滤波器获得，即

$\hat{\boldsymbol{M}}=\boldsymbol{M} \otimes \boldsymbol{K}$

(7)

式中，$ \otimes $表示卷积运算。K为卷积核函数，其大小为17×17，且定义为

$K\left({k, l} \right)=\left\{ \begin{gathered} 1/25\;\; \;\; k, l \in \left[{1, 4, 8, 12, 16} \right] \hfill \\ 0\;\; \;\; \;\; \;\; \;其他\hfill \\ \end{gathered} \right.$

(8)

式中，k、l分别为2维卷积核矩阵的下标。

根据中频信息$\hat M$，将原始遥感数据进行二值化处理。图像中的高频与低频部分为

$\boldsymbol{M}'=\boldsymbol{M}-\hat{\boldsymbol{M}}$

(9)

该二值化处理过程可以表示为

$G\left({i, j} \right)=\left\{ \begin{gathered} 1\;\; \;\; M\left({i, j} \right)-\hat M\left({i, j} \right){\text{ > }}0 \hfill \\ 0\;\; \;\; \;\; \;\; \;其他\hfill \\ \end{gathered} \right.$

(10)

式中，G(i, j)为第(i, j)像元二值化后的结果。二值化处理将灰度图像转换为仅包含0与1的二值化符号集合，基于二值化图像的异常检测算法计算量小实时性能可以得到较大的提升。

考虑到异常目标与其周围背景环境存在显著的特征差异，且具有封闭的边缘特征，而背景纹理及边缘特征则无封闭的边缘。Yuan等人^[8]对二值化高维特征图G通过局部区域内像元级的加权异或分析提取具有显著高维光谱特征的像元。最后采用局部内外窗口的均值差异判定内窗口包围的像元是否为异常目标。

然而，当异常目标仅仅占据单个像元时，通过待测像元与其局部区域内背景像元的异或运算可以得到较满意的检测结果。但当异常目标的尺寸较大时，由于目标本身光谱特征的相似性，目标中心像元与其边缘像元的异或容易产生空圈现象。导致内窗中异常目标光谱特征差异均值降低。另外，双窗口尺寸选择过小时，可能导致外窗口中包含异常像元，从而降低了内外窗口的均值差异，不利于异常目标检测。

为了避免以上情况的发生，将双窗口模式中的内窗口设定为保护区域，用于阻断异常像元与周围背景区域的混淆问题，从而更加准确地计算像元高维特征的异常显著性。双窗口模式如图 4所示。将N(i, j)划分为内窗口与外窗口两部分。外部灰色区域为外窗口，表示背景区域N_out(i, j)，紧紧围绕在待检测像元的白色区域为内窗口N_in(i, j)起到隔离目标与背景的作用。窗口尺寸的选择要根据异常目标尺寸的大小而定。

在交叉分析的异或计算中，将内窗口作为待测像元的保护带，通过比较待测像元与其外窗口N_out(i, j)之间的差异性判断该像元的属性，避免了异常目标内部像元之间的异常运算。计算公式为

$D\left({i, j} \right)=\sum\limits_{\left({m, n} \right) \in {N_{out}}\left({i, j} \right)} {\left({G\left({i, j} \right) \oplus G\left({m, n} \right)} \right){w_{m, n}}} $

(11)

式中，D(i, j)为第(i, j)个待测像元交叉分析后的结果。N_out(i, j)为像元(i, j)周围相邻的外窗口区域，待测像元G(i, j)与其背景像元G(m, n)通过⊕异或运算，寻找高维特征的相关性。若待测像元G(i, j)为异常目标，经过异或运算，(i, j)坐标对应“1”；相反，若G(i, j)为背景像元，则经交叉分析后，其值为“0”。w_{m, n}为加权系数，定义为

${w_{m, n}}=\frac{1}{{{z_{i, j}}}}\exp \left({-\frac{1}{{2{\sigma ^2}}}{{\left\| {{p_{i, j}}-{p_{m, n}}} \right\|}^2}} \right)$

(12)

式中，z_{i, j}为正则化系数，用于确保加权系数w_{m, n}在领域N_out(i, j)之中的和为1。σ=log(N)，用于控制差分计算的范围。N为局部区域N_out(i, j)内像素总数。p_{i, j}与p_{m, n}分别为第(i, j)个像元与第(m, n)个像元的坐标。加权系数w_{m, n}的值反应了相邻像素之间的交互作用是随着空间距离增加而减小，$\left\| \cdot \right\|$为L₂范数距离。

通过交叉分析，可以有效地消除背景纹理特征，在最后的异常检测阶段，对于待检测第(i, j)像元，采用单个像元的光谱亮度值与外窗口背景灰度值进行比较。由于是以单个像元为检测对象，双窗口的尺寸选择相对原方法受到的限制更小。考虑到尺寸在3×3个像元以内的异常目标，内窗口的尺寸大于等于5×5即可。

$\varphi \left({i, j} \right)=\frac{{{{\left\| {D\left({i, j} \right)-mean\left({{N_{out}}\left({i, j} \right)} \right)} \right\|}_2}}}{{mean\left({{N_{out}}\left({i, j} \right)} \right)}}$

(13)

式中，N_out(i, j)为外窗口。φ(i, j)的值大小反应了异常显著性的强弱。mean (·)为求均值函数。

采用高光谱遥感数据(图 5)与高光谱仿真遥感数据(图 6)对本文算法进行实验验证。

图 5(a)为HYDICE高光谱数据，该数据的光谱分辨率为10 nm，共160个波段，空间分辨率为1 m。截取大小为80×100像素的区域作为目标区域，选择该区域中若干尺寸较小的目标作为异常。目标的准确位置参考图如图 5(b)所示。

采用较为经典的异常检测算法：WSCF算法, 2DCAD算法及RSAD(random selection anomaly detector)算法^[12]与本文所提仿生检测算法进行比较。实验结果如图 5所示。图 5(c)为本文算法的检测结果。异常目标均被检测到。图 5(d)为2DCAD的检测结果，由于目标尺寸较小，基本均未被检测到。图 5(f)为RSAD检测结果，其中RSAD算法中涉及到的参数：查表可得自由度为df的卡方分布值为χ_0.25²=117.679和自由度df=160。RSAD算法虽然检测到部分目标，但是运行时间较长。图 5(e)为WSCF算法的检测结果，同样采用与2DCAD算法相同的5×5与7×7的内外窗口。可见该算法对小目标的敏感性较低漏检率较高。

图 6为环境一号卫星HJ-1A超光谱成像仪HSI(hyper spectral imager)数据，截取2011年7月30日拍摄的上海地区海陆交接部分数据作为实验对象。该数据的光谱分辨率达到4 nm，空间分辨率为100 m，共有115个波段。模拟存在异常目标的仿真数据如图 6(a)所示，该区域背景环境具有丰富的纹理细节特征。设置的异常目标大小分别为3×3与2×2两种尺寸。异常目标的光谱特征同样选择图 3的做法，选取该区域中具有特殊光谱的目标光谱变化特性，插入到各个波段数据中。标准的参考结果如图 6(b)所示。图 6(c)—(h)分别采用了2DCAD算法、WSCF算法、MSAD(multiple-window anomaly detector)算法^[13]、RSAD算法、SVDD算法^[14](support vector-based anomaly detector)进行检测。2DCAD算法与WSCF算法的内外窗口尺寸为5×5与7×7。MSAD算法、RSAD算法的参数为：自由度为df=116，开方表查阅阈值χ_0.25 ²=88.084。其检测结果可以准确的探测到异常目标，但是虚警率较高；SVDD算法的运算时间长，不能满足实时性要求。图 6(g)为本文算法的提取结果。

采用公开数据库blind test dataset提供的超光谱数据。该数据来源于HyMap传感器空间分辨率为3 m，光谱波段为126。拍摄于2008年7月英国Cookey市周围小镇，截取其中大小为100×200像素大小的图像块如图 7(a)所示。图 7中2DCAD算法与WSCF算法的内外窗口尺寸为5×5与7×7。考虑到该数据中，异常目标的尺寸变换范围较大，而2DCAD算法与WSCF算法均需要设定固定的内外窗口进行异常分析，因此得到的检测结果不够理想。

表 1中列举了不同数据通过不同算法实现异常检测的运行时间。可见对于本文算法，其时效性要好于RSAD以及SVDD算法。由于需要背景建模以及背景模型的不断优化，RSAD以及SVDD算法需要较长的运算时间。2DCAD与WSCF为两种实时性较好的算法。本文算法的时效性接近2DCAD算法，具有较理想的实时性。

采用ROC曲线比较不同算法检测结果的虚警率。具体数据见图 8与图 9所示。图 8中，绿色曲线代表的本文算法检测结果均优于其他算法，效果较理想。图 9中绿色曲线代表的本文算法检测结果明显优于SVDD算法、2DCAD算法、WSCF算法。RSAD算法虽然取得较大的AUC(area under curve)面积，但是结合检测结果可知，RSAD算法错误地将水体背景检测为异常。综合考察，算法具有较为理想的检测效果。

针对异常检测算法的实时性要求高，以异常目标高维光谱特征存在异常性的特点，受生物视觉非线性滤波机理的启发，提出一种快速有效的仿生异常检测算法。该算法中仿生非线性滤波器通过空域与光谱域上光谱变化特性的联合分析，有效改善现有基于滤波器的异常检测算法难以从高频信息中分离复杂背景与异常目标的不足。克服了线性滤波器无法应对复杂背景特征干扰的问题，在存在丰富纹理特征的背景区域，仍然可以获得较好的检测效果。由于不涉及到背景建模，本文算法具有较好的实时性能。然而对于非线性滤波器的优化，及波段选择问题还有待进一步的研究。