1、引　言

异常探测是高光谱遥感的重要应用之一(Chang，2003)。利用高光谱图像数据进行异常探测需求解压制后的背景统计信息，即样本协方差矩阵K(Reed和Yu，1990)或自相关矩阵R(Chang和Chiang，2002)的逆矩阵K^–1或R^–1。由于高光谱数据波段多，数据处理耗时较长，传统异常探测算法无法满足实时处理的时效性要求，如火灾火点与移动目标的实时探测。实时高光谱异常探测能够在成像光谱仪获取数据的过程中，对当前时刻所获取的待测样本立即进行快速处理，无须等待研究区所有数据获取完毕，并且实时处理的时间复杂度恒定且较短，不随数据量的增长而增长，有效提高高光谱目标探测的在线处理与应用能力。

针对摆扫型成像光谱仪(如AVIRIS)逐像元采集数据的方式，目前高光谱异常探测的实时处理方法主要为：(1) 基于背景矩阵QR分解的方法，可避免直接求解K^–1或R^–1，在保持数值计算精度的同时降低异常探测的计算复杂度(Chang 等，2001；Chang和Hsueh，2006；Du和Zhang，2011；Molero 等，2013)。(2) 基于卡尔曼滤波器的思想，利用Woodbury矩阵引理(Du和Nekovei，2009；Acito 等，2011；Chen 等，2014；Rossi 等，2014；Chang 等，2015)，结合上一时刻K^–1或R^–1与当前时刻所获取的像元，递归求解当前时刻所需K^–1或R^–1，该方法避免采用所有像元进行重复计算，降低背景压制过程的时间消耗。需要注意的是，基于多核处理器平台(FPAG，多DSP，GPU)的并行计算方法(Paz和Plaza，2010；Guo 等，2014；Zhao 等，2014)并不是实时处理，而是快速处理。

传统的实时高光谱异常探测算法存在时间复杂度高以及数值稳定性差的问题。一方面，由于计算机舍入误差的存在，基于Woodbury矩阵引理的递归算子对K^–1或R^–1进行低阶修正时易造成严重的数值精度损失(Stewart，1998；http://infoscience. epft.ch/record/161458[2016-11-21])，在机载或星上成像光谱仪不断获取数据并进行实时处理的过程中，该损失会一直存在，并且在后续探测过程中不断累积，造成严重的数值计算错误，同时，该递归算子涉及多次矩阵与向量之间的乘积操作，空间、时间等消耗依然较大。另一方面，基于矩阵QR分解的实时处理算法虽然避免了矩阵求逆运算，但该方法在实时处理过程中，每一次对当前时刻K或R进行分解所需时间复杂度为O(n³) (n为高光谱数据波段数量)，由于高光谱数据波段较多，该方法计算代价仍较高。

基于以上问题，提出一种新型的逐像元实时高光谱异常探测算法，论证高光谱遥感图像背景统计矩阵的正定对称性结构，利用Cholesky分解及其分解因子的一阶修正方法，将背景压制过程转换为下三角线性系统求解过程，快速更新当前时刻所需背景统计量，保证数值计算的稳定性，降低实时处理的时间与空间复杂度，且探测精度保持不变。

2、传统异常探测

经典高光谱异常探测(AD)算法RX(Reed和Yu，1990)以全局样本像元的协方差矩阵K作为背景统计量，记作K-AD。以全局样本像元的自相关矩阵R(Chang和Chiang，2002)为背景统计信息的异常探测算法考虑到样本一阶统计量的影响，记作R-AD。设 ${{r}}_i^{\rm{T}} = {({r_{i1}},{r_{i2}}, \cdots ,{r_{iB}})^{\rm{T}}}$ 为B×1维样本光谱矢量，B代表高光谱数据波段数量。

    (2.1) K-AD和R-AD

K-AD采用背景样本协方差矩阵，可表示为

${\delta ^{{\rm{K}} - {\rm{AD}}}}\left( {{{{r}}_i}} \right) = {\left( {{{{r}}_i} - {{\mu }}} \right)^{\rm{T}}}{{{K}}^{ - 1}}\left( {{{{r}}_i} - {{\mu }}} \right) \lessgtr \eta$

(1)

式中， ${\delta ^{{\rm{K - AD}}}}({{{r}}_i})$ 表示待测像元r_i与背景之间的马氏距离，若 ${\delta ^{{\rm{K - AD}}}}({{{r}}_i})$ 大于预设阈值η，则可判定r_i为异常像元； ${{\mu }} = (1/N)\sum\limits_{i{\rm{ = 1}}}^N {{{{r}}_i}} $ 为B×1维全局背景样本均值向量；N为像元总数； ${{K}} = (1/N)\sum\limits_{i = 1}^N {({{{r}}_i} - {{\mu }}){{({{{r}}_i} - {{\mu }})}^{\rm{T}}}} $ ，K^–1可视作背景压制(Chang，2016)。

R-AD (Chang和Chiang，2002)采用样本自相关矩阵

${\delta ^{{\rm{R}} - {\rm{AD}}}}\left( {{{{r}}_i}} \right) = {{r}}_i^{\rm{T}}{{{R}}^{ - 1}}{{{r}}_i} \lessgtr \eta$

(2)

式中， ${{R}} = (1/N)\sum\limits_{i = 1}^N {{{{r}}_i}{{r}}_i^{\rm{T}}} $ 。

    (2.2) 因果式R-AD与K-AD

由于R-AD与K-AD需要探测区域所有样本数据参与背景估计，而逐像元实时异常探测要求背景统计信息随当前时刻所采集的样本进行实时更新。需要注意的是，术语“因果”强调了在实时处理过程中，仅能采用当前时刻已获取的数据进行异常探测，未被获取的像元不能参与计算。

假设r_n为成像光谱仪当前时刻所获取的第n个样本，R(n)为当前时刻已获取的n个样本的自相关矩阵，即 ${{R}}(n) = (1/n)\sum\limits_{i{\rm{ = 1}}}^n {{{{r}}_i}{{r}}_i^{\rm{T}}} $ 。则因果式R-AD (CR-AD)可表示为

${\delta ^{{\rm{CR}} - {\rm{AD}}}}\left( {{{{r}}_n}} \right) = {{r}}_n^{\rm{T}}{{{R}}^{ - 1}}\left( n \right){{{r}}_n}$

(3)

式中， ${\delta ^{{\rm{CR - AD}}}}({{{r}}_n})$ 为当前像元r_n的异常响应值。

基于样本协方差矩阵的异常探测算法K-AD需要计算样本一阶统计量，即样本均值向量μ。假设μ(n)为当前时刻已获取的n个样本的均值向量，即 ${{\mu }}(n) = (1/n)\sum\limits_{i{\rm{ = 1}}}^n {{{{r}}_i}} $ ；K(n)为该样本集所对应的协方差矩阵，即 ${{K}}(n) = (1/n)\sum\limits_{i{\rm{ = 1}}}^n {[{{{r}}_i} - {{\mu }}(n)]{{[{{{r}}_i} - {{\mu }}(n)]}^{\rm{T}}}} $ 。则在获取r_n的情况下，CK-AD可表示为

${\delta ^{{\rm{CK}} - {\rm{AD}}}}\left( {{{{r}}_n}} \right) = {\left[ {{{{r}}_n} - {{\mu }}\left( n \right)} \right]^{\rm{T}}}{{{K}}^{ - 1}}\left( n \right)\left[ {{{{r}}_n} - {{\mu }}\left( n \right)} \right]$

(4)

随着所采集数据量的不断增大，对背景统计矩阵的估计与压制时间将不断增长，限制了当前时刻对新获取像元进行实时处理的速度(Zhao 等，2015)。

3、实时异常探测

通过对上一时刻背景统计矩阵的Cholesky分解因子进行一阶修正，在恒定时间复杂度内递归计算当前时刻背景统计信息，可有效提高实时高光谱异常探测效率。

    (3.1) CR-AD实时处理          3.1.1. 线性系统输出异常

为提高实时处理速度及数值稳定性，基于背景自相关矩阵的数据结构特性，对异常探测过程进行优化。首先，式(2)中背景自相关矩阵为对称矩阵；其次，假设x为任意B×1维非零矢量，可得

$\begin{array}{l}{{{x}}^{\rm{T}}}{{R}}(n){{x}} = {{{x}}^{\rm{T}}}\left[ {\left( {1/n} \right)\sum\limits_{i{\rm{ = 1}}}^n {{{{r}}_i}{{r}}_i^{\rm{T}}} } \right]{{x}}\\ = \left( {1/n} \right)\sum\limits_{i{\rm{ = 1}}}^n {{{\left( {{{r}}_i^{\rm{T}}{{x}}} \right)}^{\rm{T}}}\left( {{{r}}_i^{\rm{T}}{{x}}} \right)} \\ = \left( {1/n} \right)\sum\limits_{i{\rm{ = 1}}}^n {{{\left\| {{{r}}_i^{\rm{T}}{{x}}} \right\|}^{\rm{2}}}} \geqslant 0\end{array}$

(5)

式(5)表明自相关矩阵R(n)为半正定对称矩阵，其特征值 $({e_{\rm{1}}},{e_{\rm{2}}}, \cdots ,{e_B})$ 均为非负。进一步由式(2)可知，R(n)具有可逆的必要条件。一般地，在进行实时异常探测之前，需对背景自相关矩阵进行初始化，须有足够数量的初始样本保证R(n)为非奇异满秩矩阵。因此R(n)的所有特征值必须为非零值，即

${e_i} > 0\;\;\left( {i = 1,2, \cdots ,B} \right)$

(6)

因此，R(n)为对称正定矩阵。充分利用该特性将有效提高CR-AD的实时处理效率与数值计算精度。

进一步地，基于R(n)的正定对称性，利用Cholesky分解，可将R(n)快速分解为一个下三角矩阵与其转置矩阵的乘积

${{R}}\left( n \right) = {{{C}}_n}{{C}}_n^{\rm{T}}$

(7)

式中，C_n为B×1维下三角矩阵。该分解过程能够最大限度保持数值稳定性(Stewart，1998；Seeger，2004)，减少对后续实时异常探测精度的影响。

为避免当前时刻待测像元r_n对背景估计的影响，防止r_n被背景R(n)所压制，每一次获取新像元之后，采用上一时刻背景统计信息R(n–1)对r_n进行探测处理(Chang，2016)。对CR-AD进行改写，可得实时CR-AD，记作RCR-AD

$\begin{array}{l}{\delta ^{{\rm{RCR - AD}}}}\left( {{{{r}}_n}} \right) = {{r}}_n^{\rm{T}}{{{R}}^{ - 1}}\left( {n - 1} \right){{{r}}_n}\\[7pt] = {{r}}_n^{\rm{T}}{\left[ {{{{C}}_{n - 1}}{{C}}_{n - 1}^{\rm{T}}} \right]^{ - 1}}{{{r}}_n}\\[7pt] = {\left( {{{C}}_{n - 1}^{ - 1}{{{r}}_n}} \right)^{\rm{T}}}\left( {{{C}}_{n - 1}^{ - 1}{{{r}}_n}} \right)\\[7pt] = {\left\| {{{{z}}_n}} \right\|^2}\end{array}$

(8)

式中， ${{{z}}_n} = {{C}}_{n - 1}^{ - 1}{{{r}}_n}$ 为B×1维向量。

基于式(8)，求解出z_n即可输出当前时刻r_n的异常响应值。值得注意的是，为保证数值稳定性以及加速求解过程，并不需要对下三角矩阵C_n进行求逆运算，改写 ${{{z}}_n} = {{{r}}_n}/{{C}}_{n - 1}$

${{{C}}_{n - 1}}{{{z}}_n} = {{{r}}_n}$

(9)

式(9)为下三角线性方程组系统，利用前向迭代可快速求解出z_n，具体求解过程如下

${{{z}}_n}\left( k \right) = \frac{{{{{r}}_n}\left( k \right) - \sum\limits_{l = 1}^{k - 1} {{{{C}}_{n - 1}}\left( {k,l} \right){{{z}}_n}\left( l \right)} }}{{{{{C}}_{n - 1}}\left( {k,k} \right)}},\;k = \left( {1,2, \cdots ,B} \right)$

(10)

式中，z_n(k)与r_n(k)分别表示向量z_n与r_n中第k个元素， ${{{C}}_{n - 1}}(k{\rm{,}}l)$ 表示 ${{{C}}_{n - 1}}$ 中第k行第l列元素。

         3.1.2. 背景优化与实时更新

为实时处理下一时刻所获取样本，需在恒定且较短时间内更新当前时刻背景统计量R(n)

$\begin{aligned}{{R}}(n) = & \left( {1/n} \right)\sum\limits_{i = 1}^n {{{{r}}_i}{{r}}_i^{\rm{T}}} = & \left( {1 - 1/n} \right){{R}}(n - 1) + \left( {1/n} \right){{{r}}_n}{{r}}_n^{\rm{T}}\end{aligned}$

(11)

式中，R(n–1)为当前时刻已知的前(n–1)个像元的样本自相关矩阵，r_n为当前时刻新获取的样本。

由式(8)中RCR-AD可知，每一次对r_n进行实时处理时，仅需提供上一状态下R(n–1)的Cholesky分解因子 ${{{C}}_{n - 1}}$ ，而R(n–1)无需显式地参与运算。因此，为加速下一个像元的实时处理速度，将实时异常探测所需背景统计信息由R(n)转化为其对应的Cholesky分解因子C_n。在逐像元实时处理过程中，直接更新C_n，无需计算R(n)后对其重新进行Cholesky分解。

假设 ${\hat{{A}}} = {{A}} + {{x}}{{{x}}^{\rm{T}}}$ ，其中 ${\hat{{A}}}$ 与A均为正定对称矩阵，其Cholesky分解分别为： ${\hat{{A}}} = {\hat{{P}}}{{\hat{{P}}}^{\rm{T}}}$ 和 ${{A}} = {{P}}{{{P}}^{\rm{T}}}$ 。若已知Cholesky分解因子P与一维向量x，存在对Cholesky分解因子的一阶修正方法快速更新得到 ${\hat{{P}}}$ (Gill 等，1974；Stewart，1998)，无需由A与x更新 ${\hat{{A}}}$ 后重复对 ${\hat{{A}}}$ 进行复杂的Cholesky分解。

如式(11)，R(n)由上一时刻R(n–1)以及当前时刻r_n更新所得。由于矩阵 ${{{r}}_n}{{r}}_n^{\rm{T}}$ 为对称一阶矩阵(Gill 等，1974)，对R(n–1)的Cholesky分解因子 ${{{C}}_{n - 1}}$ 进行一阶修正，可快速计算当前时刻R(n)的Cholesky分解因子C_n。改写式(11)可得

$\begin{aligned}{{R}}(n) = & \left( {1 - 1/n} \right)\left[ {{{R}}(n - 1) + \left( {1/(n - 1)} \right){{{r}}_n}{{r}}_n^{\rm{T}}} \right] =\\ & \left( {1 - 1/n} \right)[{{R}}(n - 1) + {{{\tilde{{r}}}}_n}{\tilde{{r}}}_n^{\rm{T}}] =\\ & \left( {1 - 1/n} \right)[{{{C}}_{n - 1}}{{C}}_{n - 1}^{\rm{T}} + {{{\tilde{{r}}}}_n}{\tilde{{r}}}_n^{\rm{T}}] =\\ & \left( {1 - 1/n} \right){{{\tilde{{C}}}}_n}{\tilde{{C}}}_n^{\rm{T}} = {{{C}}_n}{{C}}_n^{\rm{T}}\end{aligned}$

(12)

式中， ${{\tilde{{r}}}_n} = \sqrt {1/(n - 1)} {{{r}}_n}$ ， ${{{C}}_n} = \sqrt {1 - 1/n} {{\tilde{{C}}}_n}$ 。在实时处理r_n时， ${{{C}}_{n - 1}}$ 为已知的上一时刻背景自相关矩阵R(n–1)所对应的Cholesky分解因子。 ${{\tilde{{r}}}_n}$ 可理解为对 ${{{C}}_{n - 1}}$ 进行更新的一阶修正因子。本研究采用文章(Stewart，1998)中所述的矩阵“混合旋转”算法， ${{{C}}_{n - 1}}$ 与 ${{\tilde{{r}}}_n}$ 作为输入，可快速更新求解 ${{\tilde{{C}}}_n}$ 。进一步地， ${{{C}}_n} = \sqrt {1 - 1/n} {{\tilde{{C}}}_n}$ ，得到下一时刻所需背景统计信息。

    (3.2) CK-AD实时处理

根据3.1.1中证明过程，结合样本协方差矩阵的构造方式，容易得出K(n)同样为对称正定矩阵。设K(n)的Cholesky分解为 ${{K}}(n) = {{{L}}_n}{{L}}_n^{\rm{T}}$ ，实时CK-AD (记作RCK-AD)，可表示为

$\begin{aligned}{\delta ^{{\rm{RCK - AD}}}}\left( {{{r}}_n} \right) = & {\left[ {{{{r}}_n} - {{\mu}} \left( n \right)} \right]^{\rm{T}}}{{{K}}^{ - 1}}\left( {n - 1} \right)\left[ {{{{r}}_n} - {{\mu}} \left( n \right)} \right] = \\ & \tilde {{r}}_n^{\rm{T}}{\left[ {{{{L}}_{n - 1}}{{L}}_{n - 1}^{\rm{T}}} \right]^{ - 1}}{{\tilde {{r}}}_n} = \\ & {\left( {{{L}}_{n - 1}^1{{\tilde {{r}}}_n}} \right)^{\rm{T}}}\left( {{{L}}_{n - 1}^{ - 1}{{\tilde {{r}}}_n}} \right) = {\left\| {{{{z}}_n}} \right\|^2}\end{aligned}$

(13)

式中， ${{\tilde{{r}}}_n} = [{{{r}}_n} - {{\mu }}(n)]$ 。采用如式(10)所示前向迭代过程，可快速求解z_n，输出 ${\delta ^{{\rm{RCK}} - {\rm{AD}}}}\left( {{{{r}}_n}} \right)$ 。

如式(13)，RCK-AD的实时处理方式比RCR-AD更复杂。当前时刻获取像元数据r_n后，首先需要对样本均值μ(n)进行更新，μ(n)由上一时刻均值μ(n–1)与当前时刻r_n计算所得

$\begin{aligned}{{\mu }}(n) = & \left( {1/n} \right)\sum\limits_{i = 1}^n {{{{r}}_i}} {\rm{ = }}\left( {1/n} \right)\sum\limits_{i = 1}^{n - 1} {{{{r}}_i}} + \left( {1/n} \right){{{r}}_n} = \\ & \left[ {1 - \left( {1/n} \right)} \right]\left[ {1/\left( {n - 1} \right)\sum\limits_{i = 1}^{n - 1} {{{{r}}_i}} } \right] + \left( {1/n} \right){{{r}}_n} = \\ & \left[ {1 - \left( {1/n} \right)} \right]{{\mu }}(n - 1) + \left( {1/n} \right){{{r}}_n}\end{aligned}$

(14)

进一步，对背景协方差矩阵K(n)进行改写

$\begin{aligned}{{K}}\left( n \right) = & \frac{1}{n}\sum\limits_{i = 1}^n {\left[ {{{{r}}_i} - {{\mu }}\left( n \right)} \right]{{\left[ {{{{r}}_i} - {{\mu }}\left( n \right)} \right]}^{\rm{T}}}} = \\ & \frac{1}{n}\sum\limits_{i = 1}^n {{{{r}}_i}{{r}}_i^{\rm{T}} - {{\mu }}\left( n \right){{{\mu }}^{\rm{T}}}\left( n \right)} = \\ &{{R}}\left( n \right) - {{\mu }}\left( n \right){{{\mu }}^{\rm{T}}}\left( n \right)\end{aligned}$

(15)

考虑到 ${{K}}(n) = {{{L}}_n}{{L}}_n^{\rm{T}}$ ，因此

$\begin{aligned}{{K}}\left( n \right) = & {{{L}}_n}{{L}}_n^{\rm{T}}{\rm{ = }}{{R}}\left( n \right) - {{\mu }}\left( n \right){{{\mu }}^{\rm{T}}}\left( n \right) =\\ & {{{C}}_n}{{C}}_n^{\rm{T}} - {{\mu }}\left( n \right){{{\mu }}^{\rm{T}}}\left( n \right)\end{aligned}$

(16)

式中， ${{\mu }}(n){{{\mu }}^{\rm{T}}}(n)$ 为对称一阶矩阵，L_n与C_n分别为K(n)与R(n)所对应的Cholesky分解因子，通过3.1.2中所述Cholesky分解因子的一阶修正方法，由C_n与μ(n)可快速计算得出L_n。

图1为RCR-AD以及RCK-AD的实时处理流程。

    (3.3) 时间复杂度分析

实时异常探测对算法计算性能要求苛刻，对RCR-AD及RCK-AD进行时间复杂度分析有助于评价探测器的计算性能。时间复杂度以处理过程中的计算机浮点数运算次数(flops)作为衡量指标。

不考虑背景样本自相关矩阵的初始Cholesky分解过程，RCR-AD算法可分为两个部分：

(1) 获取待测像元r_n后，结合已知的R(n–1)的Cholesky分解因子 ${{{C}}_{n - 1}}$ ，在下三角线性方程组系统中求解 ${{{z}}_{n - 1}}$ ，根据式(10)前向迭代过程，所需时间复杂度为B² flops，而向量 ${{{z}}_{n - 1}}$ 的求模运算仅涉及2B flops，可忽略不计。

(2) 为下一时刻探测所需，对 ${{{C}}_{n - 1}}$ 进行一阶修正，更新得到R(n)所对应的Cholesky分解因子C_n，所需时间复杂度为2.5B² flops(Bojanczyk 等，1987)。

因此，RCR-AD实时处理每一个当前时刻像元所需时间复杂度总计为3.5B² flops。

RCK-AD在RCR-AD处理过程的基础上，需首先去除一阶统计量μ(n)，如式(15)，该过程需要2.5B² flops。因此，RCK-AD算法实时处理当前时刻新获取的像元数据所需时间复杂度总计为6B² flops。

值得注意的是，若不采用Cholesky分解因子一阶修正方法，对每一次更新所得背景统计矩阵R(n)或K(n)重新进行Cholesky分解，仅该分解过程所需时间复杂度即为O(B³) flops(Press 等，2007)。因此，本算法大大降低实时异常探测的时间消耗。同时，因背景统计矩阵的重复求逆过程易引入严重的舍入误差(Stewart，1998；Seeger，2004)，本算法无需进行求逆运算，保持了数据处理的数值稳定性。

4、实验结果与分析

为验证所提出算法的时效性，采用实验室采集的高光谱图像数据进行计算机仿真测试。测试环境为64位Windows 10系统，Intel Core^TM i7-6770HQ 2.59 GHz CPU，仿真实验平台为Visual Studio C++ 2015。

    (4.1) 实验室FISS数据

采用的FISS数据由地面成像光谱辐射测量系统FISS在实验室采集，该系统由中国科学院遥感与数字地球研究所于2008年自主研发完成。其光谱范围为379—870 nm，共344个连续波段，光谱分辨率优于5 nm，空间分辨率为1 mrad。经过数据预处理，去除信噪比较低的波段(1—34，325—334)，剩余290个波段进行仿真实验。

图2为FISS测试数据，其中图2(a)为该数据真彩色合成图，共18个不同尺寸(1×1 cm，0.5×0.5 cm，0.2×0.2 cm)的人工目标均匀排列在一块黑色纺织布料上。目标中第1—4行为纸质目标，分别为红色(PRed)、黑色(PBla)、蓝色(PBlu)和黄色(PYel)，第5—6行为纺织布料目标，分别为淡蓝色(FLGr)和深蓝色(FDGr)。图2(b)为标记目标位置的地面验证数据。

    (4.2) 实时异常探测结果

利用实时异常探测算法(RCR-AD、RCK-AD)对上述FISS数据进行仿真处理，图3和图4为实时异常探测过程示意图。为便于目视判读，探测结果经过非线性拉伸，灰度值转换到dB空间，即灰度值 $x = 20{\lg}(x)$ ，并采用伪彩色显示。

需要注意的是，为保证背景自相关矩阵R(n)或协方差矩阵K(n)非奇异，初始状态下用于估计R(n)或K(n)的像元数量须大于R(n)或K(n)的秩，即该高光谱数据的波段数量。

为进一步评价实时异常探测的精度，利用全局异常探测器(R-AD，K-AD)进行测试，图5为4种探测器(RCR-AD，RCK-AD，R-AD和K-AD)最终的探测结果灰度图。

    (4.3) 探测精度评价

从图3和图4知，在RCR-AD和RCK-AD进行实时异常探过程中，异常目标出现之前，探测结果伪彩色图各像元强度值对比度较小；当第一个异常像元出现时，图像背景立即变暗，而异常像元则高亮显示；在所有6行目标中，第2行的异常响应最弱，与背景对比度最小。通过计算光谱相关系数，6种目标与背景之间的光谱相关性分别为0.898 (PRed)、0.962 (PBla)、0.873 (PBlu)、0.543 (PYel)、0.904 (FLGr)、0.919 (FDGr)，可见第2行黑色纸质目标(PBla)与背景相关性最大，相对于背景而言，异常程度最低。因此，第2行目标在该场景中可理解为弱信号异常目标，而其他5行目标的信号强度则较大。

进一步观察发现，图5所示4种探测器的最终探测结果灰度图中，全局性异常探测器(R-AD，K-AD)均无法提取出第2行弱信号目标，而实时异常探测器(RCR-AD，RCK-AD)在渐进式处理过程中能够有效识别出该异常信号，防止其被后续探测到的强信号压制而无法识别。

利用图2(b)地面验证数据，采用3维接收器曲线(3D-ROC)(Chang，2013)对4种探测器的探测精度进行定量分析。首先将图5探测结果图的灰度归一化至0—1，在0—1范围内变换决策阈值τ，可将该探测结果图二值化为目标和背景，结合图2(b)地面验证数据，可得到特定阈值τ所对应的探测率P_D，虚警率P_F。因此，通过变化决策阈值τ，将得到一组(τ，P_F，P_D)用于生成3D-ROC。进一步地，从3D-ROC中可提取出3条2维ROC (2D-ROC)，即2D-ROC (P_D，P_F)、2D-ROC (τ，P_D)、2D-ROC (τ，P_F)。

3D-ROC能够提供探测率P_D和虚警率P_F随决策阈值τ的变化趋势。此外，由于实际应用中异常探测为非监督过程，异常目标的识别主要依赖于对探测结果图的目视判断，因此，背景的压制效果对异常目标的决策影响较大(Chen 等，2014)。在这种情况下，2D-ROC (τ，P_F)可以提供背景压制随τ的变化而变化的信息，而2D-ROC (τ，P_D)则体现了探测率与决策阈值的变化关系。

图6为4种探测器(RCR-AD，RCK-AD，R-AD和K-AD)的ROC曲线图。为进一步定量化评价FISS数据的探测精度，表1为不同探测器的ROC曲线下面积AUC，探测性能最优的值为阴影显示，其中，A(τ，P_D)和A(P_F，P_D)越大，或A(τ，P_F)越小，探测精度越高，该异常探测器的性能越好。

从图6及表1可得出以下结论：(1) 由A(τ，P_F)可知，实时异常探测器相较于全局性异常探测器而言，其背景压制更强，但是前者探测率A(τ，P_D)比后者低，过度的背景压制易导致较低的探测率。(2) 虽然全局性异常探测器(R-AD，K-AD)具有较低的A(P_F，P_D)，但其A(τ，P_D)高于实时异常探测器(RCR-AD，RCK-AD)。因此，A(P_F，P_D)并非与A(τ，P_D)正相关，仅利用传统的2D-ROC并不能有效地评价一个异常探测算法的探测性能。而3D-ROC能够从探测率、虚警率、背景压制等多方面综合评价一个探测器的性能。

    (4.4) 算法时效性分析

算法时效性是衡量实时异常探测性能最重要的指标之一。为验证本研究所提出的实时异常探测器(RCR-AD，RCK-AD)的计算性能，对比分析了目前广泛使用的实时处理算法的计算性能并进行实验测试，主要包括基于背景统计矩阵QR分解的方法(QRR-AD、QRK-AD)(Chang 等，2001；Chang和Hsueh，2006；Hsueh，2007；Du和Zhang，2011；Molero 等，2013)；基于Woodbury引理递归更新背景统计量的方法(WR-AD、WK-AD)(Du和Nekovei，2009；Acito 等，2011；Rossi 等，2012；Chen 等，2014；Rossi 等，2014；Chang 等，2015)。

基于Woodbury引理直接计算背景自相关矩阵的逆矩阵，其递归更新公式为

${{{R}}^{ - 1}}\left( n \right) = {{{\tilde{{R}}}}^{ - 1}}\left( {n - 1} \right) - \frac{{\left[ {{{{\tilde{{R}}}}^{ - 1}}\left( {n - 1} \right){{{\tilde{{r}}}}_n}} \right]\left[ {{\tilde{{r}}}_n^{\rm{T}}{{{\tilde{{R}}}}^{ - 1}}\left( {n - 1} \right)} \right]}}{{1 + {\tilde{{r}}}_n^{\rm{T}}{{{\tilde{{R}}}}^{ - 1}}\left( {n - 1} \right){{{\tilde{{r}}}}_n}}}$

(17)

式中， ${{\tilde{{R}}}^{ - 1}}(n - 1) \! = \! {[(1 \! - \! 1/n){{R}}(n \! - \! 1)]^{ \! - \! 1}}$ ， ${{\tilde{{r}}}_n} \! = \! (1/\! \sqrt n ){{{r}}_n}$ 。

根据式(17)，计算 ${{\tilde{{R}}}^{ - 1}}(n - 1)$ 所需时间复杂度为B² flops，第2项分子部分 ${{\tilde{{R}}}^{ - 1}}(n - 1)$ 与 ${{\tilde{{r}}}_n}$ 的乘积需2B² flops，而 $[{{\tilde{{R}}}^{ - 1}}(n - 1){{\tilde{{r}}}_n}]$ 与其转置的乘积则涉及B²flops，分母部分的计算所需时间复杂度较小可忽略不计。分子与分母相除以及第1项减去第2项各涉及时间复杂度为B² flops。最后由式(3)求解当前待测像元r_n的异常响应值涉及2B² flops。因此，基于Woodbury引理递归求解的实时处理算法共所需时间复杂度为8B² flops。根据更新所得背景自相关矩阵的逆矩阵 ${{{R}}^{ - 1}}(n)$ 以及样本均值向量μ(n)可根据Woodbury引理进一步更新得到背景协方差矩阵的逆矩阵 ${{{K}}^{ - 1}}(n)$ ，该过程与式(17)类似，所需时间复杂度为5B² flops。因此，基于背景协方差矩阵进行计算，共所需浮点数运算次数为5B² flops。

基于背景自相关矩阵QR分解的方法，其步骤为：(1)对背景自相关矩阵R(n–1)进行QR分解 ${{R}}(n - 1) = {\tilde{{Q}}}{\tilde{{R}}}$ ，所需时间复杂度为2B³ flops(Golub和Van Loan，2012)；(2)基于式(3)，利用线性方程组系统求解中间变量 ${{y}}(n) = {{{R}}^{ - 1}}(n - 1){{{r}}_n}$ ，即 ${{\tilde{{Q}}}^{\rm{T}}}{{{r}}_n} = $ $ {\tilde{{R}}{{y}}}(n)$ ，该过程时间复杂度为3B² flops。(3)输出 ${{r}}_n^{\rm{T}}{{y}}(n)$ 所需时间复杂度为2B flops，可忽略不计。(4)根据式(11)，利用待测像元r_n更新背景自相关矩阵得到R(n)涉及2B² flops。因此，基于背景自相关矩阵QR分解的实时异常探测方法，处理每一个待测像元的时间复杂度总和为 $(2{B^3} + 5{B^2})$ flops。根据式(15)可由背景自相关矩阵R(n)及均值向量μ(n)计算得到协方差矩阵K(n)，该过程所需时间复杂度为2B² flops。因此，若基于背景协方差矩阵进行实时异常探测，共涉及运算量为 $(2{B^3} + 7{B^2})$ flops。

表2为上述各实时异常探测算法的理论时间复杂度，与本文所提出算法进行时效性对比分析。

对上述各算法分别进行仿真实验测试。表3为各算法实时处理单个像元、整行像元的平均时间，及处理完整幅图像所花费的时间。可以看出，基于QR分解的实时处理算法涉及时间复杂度为O(B³)，虽然避免了背景统计矩阵的求逆运算，但计算代价巨大；而基于Woodbury矩阵引理的方法需求解高维背景统计矩阵的逆，易造成数值不稳定(Stewart，1998；Seeger，2004)，且由于其多次矩阵操作，涉及时间复杂度依然较高。本文所提出的基于Cholesky分解因子一阶修正的实时处理方法所需时间复杂度最低，能够在恒定且最短时间内输出当前时刻所获取像元的异常响应值，有效提高实时高光谱异常探测器的计算性能。

5、结　论

本文提出一种基于Cholesky分解一阶修正的逐像元实时高光谱异常探测算法，定量讨论了算法的计算性能与探测精度，利用地面光谱成像辐射测量系统FISS所采集的室内高光谱数据进行仿真实时处理，得出以下结论：

(1) 通过对背景统计矩阵的Cholesky分解因子进行一阶修正，可快速更新下一时刻背景统计信息。结合下三角线性系统求解过程，实时处理时间缩短为基于QR分解算法的0.4%—0.65%，或基于Woodbury矩阵引理算法的27%—33%，且数值稳定性得以保持。

(2) 实时处理有助于观察背景压制效果的实时变化，防止弱异常信号被后续强信号所压制；另外，3D-ROC曲线有助于从背景压制、探测率、虚警率等方面综合评价探测器的性能。

本文所提出的实时高光谱异常探测算法在保持数值稳定性以及探测精度的前提下，有效提高实时处理的速度。但本文算法仅进行了计算机仿真实验，未利用机载高光谱成像仪进行真实场景下的实时处理测试，下一步将着重于无人机机载实时高光谱异常探测系统的研究。