出版日期: 2018-05-25
点击次数:
下载次数:
DOI: 10.11834/jrs.20187284
2018 | Volumn22 | Number 3
上一篇  |  下一篇


技术方法 
列式鲁棒主成分分析的高光谱遥感异常探测
expand article info 孙伟伟1,2 , 李飞1 , 杨刚1 , 张殿发1
1. 宁波大学 地理与空间信息技术系,宁波 315211
2. 武汉大学 测绘遥感信息工程国家重点实验室,武汉 430079

摘要

传统的基于鲁棒主成分分析的高光谱异常探测模型中,稀疏异常矩阵假设为非低秩且其非零元素满足随机分布条件。这导致稀疏矩阵的非零元素影响低秩背景矩阵的估计,进而制约背景信息和异常信息的有效分离。提出列式鲁棒主成分分析的异常探测方法,改进异常矩阵为列稀疏条件来解决上述问题。该方法分解高光谱影像2维矩阵为低秩背景矩阵,列稀疏异常矩阵和噪声矩阵,松弛目标方程为凸优化问题,并采用非精确增强拉格朗日乘子算法来求解得到列稀疏异常矩阵的最优估计。最后,对稀疏异常矩阵中所有列的L2范数值进行阈值分割来探测得到异常像元。利用两个高光谱影像数据集,对比5种主流的异常探测方法来验证提出方法的有效性。实验结果表明,列式鲁棒主成分分析方法优于包括传统鲁棒主成分分析模型在内的5种异常探测方法,且计算效率适中。

关键词

高光谱遥感, 异常探测, 鲁棒主成分分析, 列稀疏, 非精确增强拉格朗日乘子

Hyperspectral anomaly detection using column-wise robust principal component analysis
expand article info SUN Weiwei1,2 , LI Fei1 , YANG Gang1 , ZHANG Dianfa1
1.Department of Geography and Spatial Information Techniques, Ningbo University, Ningbo 315211, China
2.State Key Laboratory of Information Engineering in Surveying, Mapping and Remote Sensing, Wuhan 430079, China

Abstract

Hyperspectral imagery (HSI) collects the detailed spectral response of ground objects on the Earth’s surface by using hundreds of narrow bands and presents a great potential for use in detecting anomalies (i.e., small and low-probability ground objects) from the main background. Sparsity theory has recently attracted increasing interest because of intelligent processing in the hyperspectral field. Many sparsity-based anomaly detection methods have been proposed in literature, and robust principal component analysis (RPCA)-based detectors are typical examples. However, regular RPCA-based anomaly detectors show that the sparse anomaly matrix should not be of low rank, and its nonzero entries should be randomly scattered in the matrix. This condition negatively affects the estimation of a low-rank background matrix and seriously impacts the effective separation between the background and anomalies in the image scene. Moreover, RPCA-based approaches involve many iteration procedures of primal variables, which leads to high computational complexity. Therefore, this study proposed column-wise RPCA (CWRPCA) to resolve these problems and improve the detection results. CWRPCA assumes that background information exists in low-dimensional randomized column subspace and possesses low-rank properties and that the anomalies are sparse and do not lie in the column subspace of the background. CWRPCA aims to decompose the HSI data matrix into the sum of a low-rank background matrix, a column-sparse anomaly matrix with small portions of nonzero columns, and a noise matrix. In this study, when the column subspace of the background was determined, anomalies were estimated from the nonzero columns of the anomaly matrix. The problem of determining the background and anomaly matrices was formulated into a convex optimization program in (2). The inexact augmented Lagrange multiplier algorithm was implemented to optimize the objective function. This algorithm introduces two auxiliary variables into the objective function and iteratively updates primal variables by fixing other variables. The low-rank background and sparse anomaly matrices were obtained when the iteration procedure terminated. The detection result was achieved by segmenting the L2 norms of column vectors in the anomaly matrix. Two HSI datasets were used to verify the detection performance of CWRPCA. The Receiver Operating Characteristic (ROC) curve and Area Under Curve (AUC) were utilized to evaluate detection performance. ROC describes the probabilities of detection and false alarm. AUC quantifies the area under the ROC curve and shows how far the ROC curve is from the baseline. The detection results of the proposed method were compared with those of five state-of-the-art anomaly detection methods, namely, low-rank and sparse matrix decomposition-based anomaly detection method, global Reed-Xiaoli method, dual window-based eigen separation transform, collaborative representation-based detector, and low-rank and sparse representation. Experimental results showed that the proposed CWRPCA outperformed the five state-of-the-art anomaly detection methods in temrs of ROC curve and AUC with a moderate computation cost. CWRPCA is better in detecting anomalies than other methods and can be a good alternative for hyperspectral anomaly detection.

Key words

hyperspectral imagery, anomaly detection, robust principal component analysis, column sparse, inexact augmented Laplacian multiplier

1 引 言

高光谱遥感依托成像光谱仪设备,利用数百个波段来采集地表地物的细微光谱响应信息(童庆禧 等,2016张良培和李家艺,2016)。获取的高光谱遥感影像利用目标探测技术来区分低概率且细小地物与主要背景信息,并用于军事侦察、公共安全和生态监测等实践中(王跃明 等,2016杜培军 等,2016)。作为非监督的目标探测技术,异常探测不需要目标地物的光谱特征信息和复杂的影像预处理工作(如大气校正和辐射定标等),因此在高光谱遥感领域中应用更为广泛(Zhang 等,2010Zhao 等,2014Wu 等,2016)。异常探测假定影像场景中的主要地物为背景信息,细小且低概率的地物为异常信息,目的是有效区分异常信息与背景信息所占的像元。

常规的异常探测方法大多利用背景像元和异常像元的统计特性差异来探测异常信息。经典的RX(Reed-Xiaoli)方法假设背景信息满足多变量正态分布的统计特性,利用背景信息的概率分布密度来确定某一像元是否属于背景或异常信息(Reed和Yu,1990)。然而高光谱遥感影像并不一定满足正态分布假设,尤其在较为复杂的地物分布场景中(Nasrabadi,2014)。后来,学者提出许多针对RX方法的改进模型包括核RX、加权RX和子空间RX等(Nasrabadi,2014),以及一些更为先进的基于统计的异常探测算法,如基于双窗口的本征分离转换DWEST(Dual Window-based Eigen Separation Transform)(Kwon 等,2003)和多重窗口异常探测器MWAD(Multiple-Window Anomaly Detection)(Liu和Chang,2013)。以上这些方法通过压制背景信息来更好实现背景和异常信息的有效区分。然而,这些方法在估计背景统计特性时并未排除异常信息的影响,导致探测结果可能会混杂一些背景像元进而影响结果的准确性和可靠性(Sun 等,2014)。

近年来,随着压缩感知和稀疏理论的兴起,鲁棒主成分分析方法逐渐被应用于解决高光谱的异常探测问题。鲁棒主成分分析RPCA(Robust Principal Component Analysis)于2009年提出,目的为更好解决传统主成分分析中背景信息容易受到噪声和粗差影响的问题(Candes 等,2009),其目前广泛应用于视频监控、人脸识别和潜在语义索引等方面。在高光谱遥感领域,鲁棒主成分分析将高光谱遥感影像数据看作2维矩阵,采用矩阵分解的思想将其分解为一个低秩的背景矩阵和一个稀疏异常组分矩阵(以下简称稀疏矩阵),其中稀疏矩阵的非零元素包含影像的异常信息(Sun 等,2014)。目前已有学者利用鲁棒主成分分析来研究高光谱遥感的异常探测问题。如低秩和稀疏分解方法LSD(Low-rank and Sparse Decomposition)将高光谱影像划分为一系列的影像子块来分别构建鲁棒主成分分析模型,通过排列所有像元估计得到的稀疏系数的范数值来探测异常像元(Cui 等,2014)。低秩稀疏矩阵分解方法LRaSMD(Low-rank and sparse decomposition based anomaly detector)考虑高光谱遥感中的噪声问题来改进传统的鲁棒主成分分析模型,利用矩阵分解将原始高光谱遥感影像数据分解为低秩背景组分,噪声组分和稀疏异常组分,进而排列所有像元的稀疏异常组分的二阶范数值来探测异常信息(Sun 等,2014)。随后马氏距离被考虑来更好融合背景信息到LRaSMD模型中以进一步提升其探测结果的可靠性(Zhang 等,2015)。后来,低秩稀疏表达LRASR(Low-Rank And Sparse Representation)(Xu 等,2016)方法引入低秩表达模型来表达背景信息的低秩特性,利用重构背景组分的稀疏残差信息来探测异常像元。相比常规的基于统计的异常探测方法,基于RPCA的模型无须背景信息的任何统计假设,能够更好应用于高光谱遥感影像来实现异常信息的有效估计。然而,上述RPCA模型假定稀疏矩阵为非低秩且非零元素在稀疏矩阵中呈现随机分布特性,非零元素没有明显的结构特征(Candes 等,2009)。这导致稀疏矩阵的非零元素同样影响低秩背景矩阵的列元素估计(Rahmani和Atia,2015Xu 等,2010),进而降低背景信息估计与异常信息的探测结果。

因此,本文提出列式鲁棒主成分分析方法CWRPCA(Column-Wise RPCA)来尝试解决传统RPCA模型应用于高光谱遥感异常探测时存在的上述问题,采用非精确增强拉格朗日乘子算法IALM(Inexact Augmented Lagrange Multipliers)来优化求解CWRPCA模型,进而有效探测得到异常信息。采用PaviaC(Pavia City)和Cri两个高光谱数据集,对比5种主流的异常探测方法来验证方法的有效性。

2 高光谱遥感的CWRPCA模型

传统RPCA模型认为高光谱影像背景信息构成的矩阵具有低秩特性,异常信息在稀疏矩阵中具有随机分布特性。CWRPCA改进传统RPCA模型,认为异常信息在稀疏矩阵中为少数非零列即稀疏矩阵呈现列稀疏的特性,即稀疏矩阵中每一非零列对应于异常像元。另一方面,CWRPCA假设所有背景像元存在于一个低维的背景列子空间中,背景列子空间可以通过随机采样背景像元来构建得到;同时假设稀疏矩阵中的异常信息对应的非零列不存在于该子空间中。稀疏矩阵中的非零列与部分背景像元对应列不相关,这使得低秩背景矩阵的部分列元素的估计并不受到稀疏矩阵中非零列的影响,进而改善背景信息与异常信息的分离结果(Rahmani和Atia,2015)。假设高光谱影像数据形成2维矩阵 ${ Y} = \left\{ {{{{y}}_i}} \right\}_{i = 1}^N \in {{\bf{R}} ^{D \times N}}$ ,式中,D为波段数量,N为像素个数, ${{{y}}_i}$ 为任一像元的光谱向量。考虑到高光谱遥感成像过程中的高斯随机噪声,CWRPCA将原始高光谱2维矩阵 ${ Y}$ 分解为一个低秩背景矩阵 ${ B} \in {{\bf{R}} ^{D \times N}}$ ,一个噪声矩阵 ${ E} \in {{\bf{R}} ^{D \times N}}$ 和一个列稀疏异常矩阵 ${ S} \in {{\bf{R}}^{D \times N}}$

${ Y} = { B} + { S} + { E},{\rm{s}}.{\rm{t}}.\left\{ \begin{array}{l}{\rm{rank}}({ B}) = r\\[3pt]{\left\| { S} \right\|_{{\rm{column}},0}} \leqslant K\end{array} \right.$ (1)

式中, ${\rm{rank}}({ B}) = r > 0$ 是为限定背景矩阵的低秩特性, ${\left\| { S} \right\|_{{\rm{column}},0}} \leqslant K$ 是为保证列稀疏异常矩阵 $ S$ 中非零列的数量不超过K。可以看出,稀疏列矩阵中 $ S$ 的的非零列对应于高光谱影像中的异常信息,同时这些异常信息在背景矩阵 $ B$ 的对应位置应该为0。式(1)可转换为以下目标函数的优化求解问题

$\begin{array}{l}\mathop {\arg \min }\limits_{{{B}},{{S}}} \;{\rm{rank}}( B) + {\textit{λ}} {\left\| { S} \right\|_{0,2}}\\\;\;{\rm{s}}.{\rm{t}}.{\left\| {{ Y} - { B} -{ S}} \right\|_{\rm F}} \leqslant \varepsilon \end{array}$ (2)

式中, ${\left\| \cdot \right\|_{\rm F}}$ 代表弗罗贝尼乌斯范数(Frobenius范数), ${\left\| { S} \right\|_{0,2}}$ 为稀疏矩阵中非零列的约束项, ${\textit{λ}} $ 为正则化参数以平衡低秩约束和列稀疏约束的贡献大小, $\varepsilon $ 为设定的逼近误差值或噪声范数约束值。

3 基于CWRPCA的高光谱异常探测

3.1 CWRPCA模型的优化求解

式(2)是一个NP-hard的非凸优化问题,因此较难求解得到低秩背景矩阵 $ B$ 和列稀疏矩阵 $ S$ 的全局最优解。理论证明,核范数能够很好地逼近低秩约束项 ${\rm{rank}}({ B})$ (Sun 等,2017),同时 ${{\rm{L}}_{1,2}}$ 混合范数可以作为矩阵 ${{\rm{L}}_{0,2}}$ 范数的最佳近似(Liu 等,2013)。因此式(2)被松弛为以下凸优化问题

$\begin{array}{l}\mathop {\arg \min }\limits_{{{B}},{\rm{S}}} \;{\left\| { B} \right\|_ * } + {\textit{λ}} {\left\| { S} \right\|_{1,2}}\\\;\;{\rm{s}}.{\rm{t}}.{\left\| {{ Y} -{ B} - { S}} \right\|_{\rm F}} \leqslant \varepsilon \end{array}$ (3)

式中, ${\left\| \cdot \right\|_ * }$ 为矩阵核范数运算符, ${\left\| { S} \right\|_{1,2}} \!=\!\! \sum\limits_{i = 1}^N {({{\left\| {{ S}(:,i)} \right\|}_2})} $ ${{\rm{L}}_{1,2}}$ 混合范数,表示S中每一列向量 ${{\rm{L}}_2}$ 的范数之和。前项 ${\left\| { B} \right\|_ * }$ 限定背景矩阵的低秩特性,后项 ${\left\| { S} \right\|_{1,2}}$ 保证异常矩阵的列稀疏特性。式(3)的求解可以采用离群值追踪算法(Xu 等,2010)、调制器追踪算法(Chen 等,2011)和加速邻近梯度算法(Lin 等,2009)等。本文中,考虑到计算简单和方便,采用非精确增强拉格朗日乘子(IALM)算法(Lin 等,2010)来求解得到式(3)的全局优化解。IALM通过引入辅助变量来增强式(3)的约束条件,转换式(3)为一系列优化主变量和拉格朗日乘子的目标子方程来依次进行迭代求解。首先,引入辅助变量 ${ J} = { B}$ ,式(3)转换为

$\begin{array}{l}\mathop {\arg \min }\limits_{{{B}},{{S}}} \;{\left\| { J} \right\|_ * } + {\textit{λ}} {\left\| { S} \right\|_{1,2}}\\\;\;{\rm{s}}.{\rm{t}}.{\left\| {{ Y} - { B} - { S}} \right\|_{\rm{F}}} \leqslant \varepsilon ,{ J} = { B}\end{array}$ (4)

进一步,引入两个拉格朗日乘子 ${{ Z}_1}$ ${{ Z}_2}$ ,式(4)可以转换为如下问题

$\begin{array}{l}\mathop {\arg \min }\limits_{{{B}},{{S}},{{J}},{{{Z}}_1},{{{Z}}_2},\beta } \;{\left\| { J} \right\|_ * } + {\textit{λ}}{\left\| { S} \right\|_{1,2}} + {\rm{tr}}({ Z}_1^{\rm T}({ Y} - { B} - { S})) + \\{\rm{tr}}({ Z}_2^{\rm T}({ J} - { B})) + \displaystyle\frac{\beta }{2}(\left\| {{ Y} - { B} - { S}} \right\|_{\rm F}^2 + \left\| {{ J} - { B}} \right\|_{\rm F}^2)\end{array}$ (5)

式中, ${\rm{tr}}( \cdot )$ 为矩阵的迹运算符, $\displaystyle\frac{\beta }{2}(\left\| {{ Y} - { B} - { S}} \right\|_{\rm F}^2 + $ $ \left\| {{ J} - { B}} \right\|_{\rm F}^2)$ 为惩罚项以增强目标函数的凸特性且并不影响其全局最优解估计, $\, \beta $ 为设定的惩罚正则化参数,其随着迭代过程而发生改变。假设初始条件为 ${{ J}^{(0)}} = {{ B}^{(0)}} = 0$ ${{ S}^{(0)}} = 0$ ${ Z}_1^{(0)} = 0$ ${ Z}_2^{(0)} = 0$ $\, {\beta ^{(0)}} = {10^{ - 6}}$ 。在t+1次迭代中,当固定其他变量时,变量 ${{ J}^{(t + 1)}}$ 可以通过求解以下目标子式(6)来实现

${{ J}^{(t + 1)}} = \arg \min \frac{1}{{{\beta ^{(t)}}}}{\left\| {{J}} \right\|_*} + \frac{1}{2}\left\| {{ J} - ({{ B}^{(t)}} + { Z}_2^{(t)}/{\beta ^{(t)}})} \right\|_{\rm F}^2$ (6)

进一步,当固定其他变量时,变量 ${B^{(t + 1)}}$ 可以通过求解以下目标子式(7)来实现

$\begin{aligned}{{ B}^{(t + 1)}} = & \arg \min \displaystyle\frac{1}{2}\left\| {{ B} - \left( {{ Y} - {{ S}^{(t)}} + { Z}_1^{(t)}/{\beta ^{(t)}}} \right)} \right\|_{\rm F}^2 +\\ &\displaystyle\frac{1}{2}\left\| {{ B} - \left( {{{ J}^{(t + 1)}} - { Z}_2^{(t)}/{\beta ^{(t)}}} \right)} \right\|_{\rm F}^2\end{aligned}$ (7)

接下来,固定其他变量,变量 ${{ S}^{(t + 1)}}$ 可以通过求解以下目标子式(8)来实现

$\begin{aligned}{{ S}^{(t + 1)}} = & \arg \min \displaystyle\frac{{\textit{λ}} }{\beta }{\left\| { S} \right\|_{1,2}} +\\ &\displaystyle\frac{1}{2}\left\| {{ S} - ({ Y} - {{ B}^{(t + 1)}} + { Z}_1^{(t)}/{\beta ^{(t)}})} \right\|_{\rm F}^2\end{aligned}$ (8)

最后,两个拉格朗日乘子 ${ Z}_1^{(t{\rm{ + }}1)}$ ${ Z}_2^{(t{\rm{ + }}1)}$ 及惩罚正则化参数 ${\beta ^{(t{\rm{ + }}1)}}$ 分别通过式(9)—(11)来依次进行更新

${ Z}_1^{(t + 1)} = { Z}_1^{(t)} + {\beta ^{(t)}}({ Y} - {{ B}^{(t + 1)}} - {{ S}^{(t + 1)}})$ (9)
${ Z}_2^{(t + 1)} = { Z}_2^{(t)} + {\beta ^{(t)}}({{ B}^{(t + 1)}} - {{ J}^{(t + 1)}})$ (10)
${\beta ^{(t + 1)}} = \min (1.1{\beta ^{(t)}},{10^{10}})$ (11)

上述迭代过程持续更新直至满足以下条件: ${\left\| {{{ B}^{(t + 1)}} - {{ J}^{(t + 1)}}} \right\|_\infty } < \tau $ ${\left\| {{ Y} - {{ B}^{(t + 1)}} - {{ S}^{(t + 1)}}} \right\|_\infty } < \tau $ ;或者最大迭代次数 ${t_{\max }} > {10^6}$ ${\left\| \cdot \right\|_\infty }$ 表示无穷范数, $\tau $ 为设定的逼近误差。迭代结束,最终得到低秩背景矩阵和列稀疏异常矩阵的最优解 ${ B} = {{ B}^{(t + 1)}}$ ${ S} = $ ${{ S}^{(t + 1)}}$

3.2 基于列稀疏矩阵S的高光谱异常信息估计

由式(1)可知,求解得到的列稀疏异常矩阵S非零列对应于原始高光谱影像中的异常信息,其每一列对应一个异常像元。然而,实际计算过程中,由于噪声等因素的影响,非零列的数量会往往大于实际异常像元的个数。因此,采用列稀疏矩阵S的每一列向量的L2范数来求解得到各个像元的异常值,通过设定阈值η来分割得到异常像元集合 ${ Y}(:,C)$

${ Y}(:,C) = \left\{ {C\left| {{{\left\| {{{ S}_{(:,C)}}} \right\|}_2} > \eta } \right.} \right\},C \in [1,N],\eta > 0$ (12)

在此基础上得到CWRPCA的2维异常探测结果。

3.3 基于CWRPCA的异常探测流程

基于CWRPCA的高光谱异常探测流程如图1所示,包括以下步骤:

图 1 CWRPCA用于异常探测的流程
Fig. 1 The process of hyperspectral anomaly detection using CWRPCA

(1) 将3维高光谱影像数据立方体转换为2维矩阵 ${ Y} \in {{\bf{R}} ^{D \times N}}$ ,式中D为波段数量,N为像素个数;

(2) 利用式(2)来构造CWRPCA模型,并通过范数松弛转换为凸优化问题(3);

(3) 采用IALM算法,利用式(5)—(11)迭代更新各个变量,求解得到背景矩阵和列稀疏异常矩阵的最优解 ${ B} = {{ B}^{(t + 1)}}$ ${ S} = {{ S}^{(t + 1)}}$

(4) 采用式(12)来求解稀疏矩阵S中各个列向量的L2范数并通过阈值分割得到异常探测图。

4 实验和分析

4.1 实验数据

PaviaC数据由ROSIS光学系统成像光谱仪采集得到,影像覆盖意大利北部的帕维亚大学区域,具有准确的真实地物空间分布信息(Sun 等,2014Zhang 等,2015)。PaviaC数据的原始波段数为115,光谱范围为430—860 nm,空间分辨率为1.3 m。实验采用原始影像中较小的一块区域,大小为108×120像素;对原始115波段进行预处理,移除13个低信噪比的波段,剩余102波段。图2(a)中3种主要地物构成影像场景的背景信息:桥、水体和阴影。代表桥上车辆的47个像元被认为是场景中的异常信息,因为车辆的光谱特征曲线与背景地物差异较大。图2(b)展示异常像元的实际分布情况,图2(c)展示异常像元和背景像元的光谱响应曲线信息。

图 2 PaviaC数据信息
Fig. 2 The information of PaviaC dataset

Cri数据由Nuance Cri成像光谱仪采集得到,具有46个光谱波段,覆盖光谱区间为650—1100 nm(Zhang 等,2015)。实验采用原始影像中较小的一块区域,大小为400×400像素(图3(a))。10块石头散落在草地中,占据图像2216个像元。这些石块与主要背景地物的光谱特征差异很大,被认为是场景中的异常信息。图3(b)展示异常像元的实际分布情况,图3(c)展示异常像元和背景像元的光谱响应曲线信息。

图 3 Cri数据信息
Fig. 3 The information of Cri dataset

4.2 实验分析

本节利用以上两个高光谱影像数据集来设计两组实验以验证提出的CWRPCA异常探测方法的有效性。第一组实验探讨不同的正则化参数λ对CWRPCA的异常探测性能的影响,试图为后续应用中选取合适的正则化参数λ提供依据。第二组实验对比CWRPCA与5种主流的异常探测方法的探测结果。对比方法包括全局RX算法GRX(Global RX)、两种局部方法双窗口的本征分离转换算法(DWEST)(Kwon 等,2003)和协同表达探测算法CRD(Collaborative Representation based Detector)(Li和Du,2015),基于传统RPCA模型的低秩稀疏矩阵分解的异常探测算法(LRaSMD)(Sun 等,2014),以及低秩稀疏表达方法(LRASR)(Xu 等,2016)。采用受试者工作特征ROC(Receiver Operating Characteristic)曲线和ROC曲线的线下面积AUC(Area Under Curve)来定量评价异常探测结果。ROC曲线能够同时描述探测概率和虚警概率,而AUC能够定量评价ROC曲线偏离其理想ROC曲线(AUC=1)的程度。

(1)正则化参数λ对CWRPCA的影响分析

本实验通过改变不同的正则化参数λ来分析其对CWRPCA的探测性能的影响。PaviaC和Cri数据中,正则化参数λ的选择区间人为设定为[0.001, 0.005, 0.01, 0.02, 0.05]。

图4展示不同正则化参数λ情况下的CWRPCA在PaviaC和Cri数据中得到的ROC曲线。从图4(a)可以看出,随着λ从0.001增大至0.05,CWRPCA的ROC曲线所表现的探测性能逐渐下降。当λ=0.05时ROC曲线表现最差,达到100%的探测率时具有最大的虚警率(f)。类似地,图4(b)中,ROC曲线所表现出来的CWRPCA探测性能总体呈现下降趋势。当λ从0.001升至0.005时,ROC曲线性能明显提升;在λ=0.005时表现最佳,达到100%的探测率时具有最小的虚警率(f)。随着λ持续增大至0.05时,CWRPCA的探测性能逐渐下降。当λ=0.05时,ROC曲线表现最差。另一方面,表1列出不同λ情况下CWRPCA的ROC曲线的AUC值。可以看出,AUC随着λ的增加整体呈下降趋势,与ROC曲线的结果变化情况一致。因此,可以得出如下结论:正则化参数λ的选取对CWRPCA的探测性能影响较大;较小的λ能够保证CWRPCA能够得到较好的探测结果。

图 4 不同正则化参数λ情况下的ROC曲线
Fig. 4 The ROC curves of CWRPCA with different choices of regularized parameter λ

表 1 不同正则化参数λ情况下的AUC值
Table 1 The AUC values of CWRPCA with different choices of λ

下载CSV 
高光谱影像 线下面积AUC
λ=0.001 λ=0.005 λ=0.01 λ=0.02 λ=0.05
PaviaC数据 0.9997 0.9996 0.9994 0.9992 0.9878
Cri数据 0.8890 0.9836 0.9626 0.8724 0.4999

(2)CWRPCA的探测性能分析

本实验对比CWRPCA与其他5种主流异常探测方法的探测性能。利用交叉验证方法,PaviaC和Cri数据中,DWEST方法的最大窗口大小设置分别为9×9和21×21。利用交叉验证方法,CRD方法中,PaviaC数据的外窗和内窗大小分别设置为5×5和3×3;Cri数据的外窗和内窗大小分别设置为17×17和3×3。利用交叉验证方法,LRaSMD方法中,PaviaC数据的秩数和正则化参数分别设置为2和0.4;Cri数据的秩数和正则化参数分别设置为3和0.6。PaviaC和Cri数据中,LRASR的正则化参数都设置为0.1。利用交叉验证,PaviaC和Cri数据中,CWRPCA的正则化参数分别设置为0.001和0.005。

图5为CWRPCA和其他5种方法在PaviaC和Cri数据中得到的ROC曲线。图5(a)中,CWRPCA的ROC曲线明显优于其他5种方法,在同样的虚警率(f)条件下始终具有最高的探测率。LRaSMD方法的ROC优于DWEST、CRD、GRX和LRASR。GRX的ROC曲线优于两种局部方法DWEST和CRD。DWEST和CRD表现不佳的原因解释如下。47个像元的异常信息位于细长的桥梁中部或边界,邻域窗口内背景地物较为复杂,导致其并不满足高斯空间分布假设。6种异常探测方法中,LRASR方法表现最差,探测性能明显低于其他5种方法。图5(b)中,CWRPCA同样优于其他5种方法,在达到100%的探测率时具有最低的虚警率(f)。DWEST方法的ROC曲线优于GRX、LRASR和LRaSMD。表2列出6种异常探测方法在PaviaC和Cri数据中的AUC值。表2中可以看出,CWRPCA方法获得最大的AUC值;LRaSMD方法的AUC大于GRX和LRASR。尤其CWRPCA的ROC曲线和AUC优于LRaSMD,这说明CWRPCA模型能够改善背景信息与异常信息的估计,进而提升异常探测的性能。图6图7列出6种方法在PaviaC和Cri数据中得到的阈值分割前的2维探测图。

图 5 6种异常探测方法的ROC曲线对比
Fig. 5 The contrast in ROC curves from six anomaly detection methods

表 2 6种异常探测方法的线下面积AUC对比
Table 2 The contrast in AUC from six anomaly detection methods

下载CSV 
数据集 异常探测方法
GRX DWEST CRD LRaSMD LRASR CWRPCA
PaviaC数据 0.9946 0.9871 0.9822 0.9979 0.9643 0.9997
Cri数据 0.9135 0.9206 0.9612 0.9158 0.7057 0.9836
图 6 PaviaC数据的阈值分割前的2维探测图
Fig. 6 The detection maps on PaviaC data without thresholding segmentation from GRX, DWEST, CRD, LRaSMD, LRASR, CWRPCA
图 7 Cri数据的阈值分割前的2维探测图
Fig. 7 The detection maps on Cri data without thresholding segmentation from GRX, DWEST, CRD, LRaSMD, LRASR, CWRPCA

此外,表3列出6种异常探测方法在PaviaC和Cri数据上的运行时间。6种异常探测方法都通过MALTAB 2014a编程实现,运算环境为Inter(R) CORE(TM) i5-3230M双核处理器,内存为8 GB和Windows 10操作系统。可以看出,GRX的计算效率最高,所需的计算时间最短。LRaSMD计算效率高于CWRPCA、DWEST、LRASR和CRD。LRASR由于包含过多迭代过程,导致计算效率最低。DWEST和CRD的计算效率较低,这是由于遍历所有像元的局部背景信息统计所需的计算时间较长。因此,可以得出如下结论:5种异常探测方法中,CWRPCA的探测性能最优,优于基于传统RCPA的LRaSMD方法;CWRPCA的计算效率适中,LRASR的计算效率最低。

表 3 6种异常探测方法的计算时间对比
Table 3 The contrast in computation times from six anomaly detection methods

下载CSV 
数据集 异常探测方法
GRX DWEST CRD LRaSMD LRASR CWRPCA
PaviaC数据 0.144 87.418 29.802 14.116 578.52 60.369
Cri数据 1.902 1271.722 2371.417 85.243 3242.503 364.591

5 结 论

本文提出列式鲁棒主成分分析CWRPCA模型来解决传统RPCA用于异常探测时的不足,降低稀疏异常矩阵中非零元素对背景信息估计的影响,增强背景信息和异常信息的分离程度以提升异常探测效果。CWRPCA将非凸优化问题松弛为拉格朗日函数的优化求解问题,采用非精确增强拉格朗日乘子算法IALM来求解得到列稀疏的异常矩阵,进而通过其列向量的L2范数值的阈值分割来实现异常探测。基于PaviaC和Cri数据来设计两组实验以验证CWRPCA的探测性能。实验结果表明,CWRPCA的探测性能明显优于基于LRaSMD和其他4种方法GRX,DWEST,LRASR与CRD,且所需的计算时间适中。此外,一个较小的正则化参数λ能够使CWRPCA具有较好的异常探测结果。

参考文献(References)

  • Candes E J, Li X D, Ma Y and Wright J. 2009. Robust principal com-ponent analysis? arXiv:0912.3599
  • Chen Y D, Xu H, Caramanis C and Sanghavi S. 2011. Robust matrix completion and corrupted columns//Proceedings of the 28th International Conference on Machine Learning. Bellevue, WA, USA: ICML: 873–880.
  • Cui X Q, Tian Y, Weng L B and Yang Y P. 2014. Anomaly detection in hyperspectral imagery based on low-rank and sparse decomposition//Proceedings of SPIE Volume 9069, Fifth International Conference on Graphic and Image Processing. Hong Kong, China: SPIE: 90690R [DOI: 10.1117/12.2050229]
  • Du P J, Xia J S, Xue C H, Tan K, Su H J and Bao R. 2016. Review of hyperspectral remote sensing image classification. Journal of Remote Sensing, 20 (2): 236–256. [DOI: 10.11834/jrs.20165022] ( 杜培军, 夏俊士, 薛朝辉, 谭琨, 苏红军, 鲍蕊. 2016. 高光谱遥感影像分类研究进展. 遥感学报, 20 (2): 236–256. [DOI: 10.11834/jrs.20165022] )
  • Kwon H, Der S Z and Nasrabadi N M. 2003. Adaptive anomaly detection using subspace separation for hyperspectral imagery. Optical Engineering, 42 (11): 3342–3351. [DOI: 10.1117/1.1614265]
  • Li W and Du Q. 2015. Collaborative representation for hyperspectral anomaly detection. IEEE Transactions on Geoscience and Remote Sensing, 53 (3): 1463–1474. [DOI: 10.1109/TGRS.2014.2343955]
  • Lin Z, Ganesh A, Wright J, Wu L, Chen M and Ma Y. 2009. Fast convex optimization algorithms for exact recovery of a corrupted low-rank matrix. Proceedings in International workshops on Computational Advances in Multi-Sensor Adaptive Processing (CAMSAP), 1–18 [DOI: 10.1.1.231.4169]
  • Lin Z C, Chen M M and Ma Y. 2010. The augmented Lagrange multiplier method for exact recovery of corrupted low-rank matrices. arXiv:1009.5055
  • Liu G C, Lin Z C, Yan S C, Sun J, Yu Y and Ma Y. 2013. Robust recovery of subspace structures by low-rank representation. IEEE Transactions on Pattern Analysis and Machine Intelligence, 35 (1): 171–184. [DOI: 10.1109/TPAMI.2012.88]
  • Liu W M and Chang C I. 2013. Multiple-window anomaly detection for hyperspectral imagery. IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing, 6 (2): 644–658. [DOI: 10.1109/JSTARS.2013.2239959]
  • Nasrabadi N M. 2014. Hyperspectral target detection: an overview of current and future challenges. IEEE Signal Processing Magazine, 31 (1): 34–44. [DOI: 10.1109/MSP.2013.2278992]
  • Rahmani M and Atia G. 2015. Randomized robust subspace recovery for high dimensional data matrices. arXiv:1505.05901
  • Reed I S and Yu X. 1990. Adaptive multiple-band CFAR detection of an optical pattern with unknown spectral distribution. IEEE Transactions on Acoustics, Speech, and Signal Processing, 38 (10): 1760–1770. [DOI: 10.1109/29.60107]
  • Sun W W, Liu C, Li J L, Lai Y M and Li W Y. 2014. Low-rank and sparse matrix decomposition-based anomaly detection for hyperspectral imagery. Journal of Applied Remote Sensing, 8 (1): 083641 [DOI: 10.1117/1.JRS.8.083641]
  • Sun W W, Yang G, Du B, Zhang L F and Zhang L P. 2017. A sparse and low-rank near-isometric linear embedding method for feature extraction in hyperspectral imagery classification. IEEE Transactions on Geoscience and Remote Sensing, 55 (7): 4032–4046. [DOI: 10.1109/TGRS.2017.2686842]
  • Tong Q X, Zhang B and Zhang L F. 2016. Current progress of hyperspectral remote sensing in China. Journal of Remote Sensing, 20 (5): 689–707. [DOI: 10.11834/jrs.20166264] ( 童庆禧, 张兵, 张立福. 2016. 中国高光谱遥感的前沿进展. 遥感学报, 20 (5): 689–707. [DOI: 10.11834/jrs.20166264] )
  • Wang Y M, Jia J X, He Z P and Wang J Y. 2016. Key technologies of advanced hyperspectral imaging system. Journal of Remote Sensing, 20 (5): 850–857. [DOI: 10.11834/jrs.20166206] ( 王跃明, 贾建鑫, 何志平, 王建宇. 2016. 若干高光谱成像新技术及其应用研究. 遥感学报, 20 (5): 850–857. [DOI: 10.11834/jrs.20166206] )
  • Wu Z B, Li Y L, Plaza A, Li J, Xiao F and Wei Z H. 2016. Parallel and distributed dimensionality reduction of hyperspectral data on cloud computing architectures. IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing, 9 (6): 2270–2278. [DOI: 10.1109/JSTARS.2016.2542193]
  • Xu H, Caramanis C and Sanghavi S. 2010. Robust PCA via outlier pursuit//Proceedings of the 23rd International Conference on Neural Information Processing Systems. Vancouver, British Columbia, Canada: Curran Associates Inc.: 2496–2504
  • Xu Y, Wu Z B, Li J, Plaza A and Wei Z H. 2016. Anomaly detection in hyperspectral images based on low-rank and sparse representation. IEEE Transactions on Geoscience and Remote Sensing, 54 (4): 1990–2000. [DOI: 10.1109/TGRS.2015.2493201]
  • Zhang L P, Du B and Zhong Y F. 2010. Hybrid detectors based on selective endmembers. IEEE Transactions on Geoscience and Remote Sensing, 48 (6): 2633–2646. [DOI: 10.1109/TGRS.2010.2040284]
  • Zhang L P and Li J Y. 2016. Development and prospect of sparse representation-based hyperspectral image processing and analysis. Journal of Remote Sensing, 20 (5): 1091–1101. [DOI: 10.11834/jrs.20166050] ( 张良培, 李家艺. 2016. 高光谱图像稀疏信息处理综述与展望. 遥感学报, 20 (5): 1091–1101. [DOI: 10.11834/jrs.20166050] )
  • Zhang Y X, Du B and Zhang L P. 2015. A sparse representation-based binary hypothesis model for target detection in hyperspectral images. IEEE Transactions on Geoscience and Remote Sensing, 53 (3): 1346–1354. [DOI: 10.1109/TGRS.2014.2337883]
  • Zhao R, Du B and Zhang L P. 2014. A robust nonlinear hyperspectral anomaly detection approach. IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing, 7 (4): 1227–1234. [DOI: 10.1109/JSTARS.2014.2311995]