CN113537062B

CN113537062B - 一种基于FrFT变换和全变分正则化的异常检测方法

Info

Publication number: CN113537062B
Application number: CN202110805375.5A
Authority: CN
Inventors: 孙菲; 厉小润
Original assignee: Zhejiang University ZJU
Current assignee: Zhejiang University ZJU
Priority date: 2021-07-16
Filing date: 2021-07-16
Publication date: 2024-01-09
Anticipated expiration: 2041-07-16
Also published as: CN113537062A

Abstract

本发明属于遥感图像处理技术领域，具体为一种基于FrFT变换和全变分正则化的异常检测方法。首先，通过聚类算法将图像高维数据映射至多个子空间，构造FrFT‑RX算子，增大背景和异常的可分性，得到较为纯净的背景字典。然后，将高光谱数据通过FrFT变换映射至中间域，在低秩稀疏模型中引入全变分正则化约束表示中间域内背景部分的空间平滑性，建立异常检测模型；最后，通过交替方向乘子法把模型求解转换成求解多个子问题的最优解，得到稀疏部分，计算稀疏部分的范数得到最终的检测结果。本发明综合利用字典学习概念，降低异常原子的干扰，算法具有更低的虚警率；同时联合图像空谱信息，建立新的异常检测模型，算法具有更高的检测率。

Description

一种基于FrFT变换和全变分正则化的异常检测方法

技术领域

本发明属于遥感图像处理技术领域，具体设计一种基于FrFT变换和全变分正则化的异常检测方法，用于高光谱异常检测。

背景技术

不同于红外图像和多光谱图像，高光谱图像具有数百个窄带，它包含一段连续光谱范围内地物的辐射值，其光谱分辨率可达纳米级。作为高光谱图像处理的重要组成部分，异常检测旨在以无监督的方式识别高维数据中的异常信息。其中，异常是指在空间上稀疏分布的少量像素或者混合像元，具有和背景像素不同的光谱特征。

目前，已经提出了各种用于高光谱异常检测的方法。RX(Reed-Xiaoli，RX)算法假设背景像素遵循多元正态分布统计模型，通过计算测试像素与相邻背景像素之间的马氏距离来检测异常。但在实际应用中，由于场景地物类型复杂，不符合多元正态分布模型的条件。此外，受噪声和异常点的影响，以整幅图像的协方差矩阵估计背景统计特征并不准确。

为避免传统方法的上述缺陷，一些没有基于高斯分布的检测模型也被提出。它们主要围绕着高光谱数据的异常稀疏性、背景低秩性，提出更加符合实际场景的数据模型。基于表示的检测模型把高光谱图像看作由背景和异常两部分组成，基于背景的低秩性和目标的稀疏性通过矩阵分解重构残差。鲁棒主成分分析(Robust Principal ComponentAnalysis，RPCA)将图像分解为低秩的背景矩阵和稀疏的异常矩阵，通过适当松弛求解凸优化问题。低秩稀疏矩阵分解(Low-Rank and Sparse Matrix Decomposition Detector，LRaSMD)考虑噪声部分，把高光谱数据看做由低秩部分、噪声部分和稀疏部分组成，通过Godec算法快速迭代最小化噪声部分分解原始数据。随后zhang等人结合马氏距离检测器融合背景低秩部分进一步提高异常检测的准确性。在此基础上，低秩表示模型(Low-RankRepresentation，LRR)改进RPCA单一子空间表示模型，根据高维数据分布特点提出多子空间模型，将复杂背景用多个子空间线性混合表示。Yang等人对低秩系数添加局部约束的稀疏项，提出低秩稀疏表示模型(Low Rank and Sparse Representation，LRASR)。但是，基于表示的检测模型忽略了噪声和异常对构建背景字典过程中的污染问题；它只关注高光谱图像光谱维信息，忽略图像的空间信息，这两点在一定程度上限制最终的检测精度。

发明内容

针对低秩稀疏表示的高光谱异常检测算法中背景字典易被污染、空间信息利用不足等问题，本发明提出一种基于FrFT变换和全变分正则化的高光谱异常检测方法来解决此问题。首先，通过FrFT变换提取高光谱图像的异常特征，抑制噪声干扰；其次，构造FrFT-RX算子剔除异常原子，得到纯化的背景字典；然后，引入全变分正则化项约束，建立全变分正则化约束项的低秩稀疏模型，以表示经过FrFT变换后中间域内图像的空间信息；通过交替方向乘子算法最优化目标函数。最终，计算异常部分的l_2,1范数得到异常检测结果。与传统的异常检测算法相比，本发明在多个数据集上具有更高的检测率和更低的虚警率，具有一定的鲁棒性。

为了实现上述目的，本发明具体采用了如下技术方案：

一种基于FrFT变换和全变分正则化的高光谱异常检测方法，步骤如下：

1)获得纯净的背景字典：

对于高光谱图像通过k-means聚类算法，将图像高维数据中光谱信息相似的像元映射到一个子空间中：

Y＝∩_i＝1,…,MYⁱ

其中，Y∈R^B×N表示高光谱二维矩阵形式，Yⁱ表示第i个子空间，M表示子空间个数，两两子空间之间无交集；y_j表示像元点，B表示图像波段数，N表示每个波段的像元点数；

构造FrFT-RX算子：

K_p(u,v)＝A_αexp{jπ(u²cotα-2uv cscα+v²cotα)}

其中，y_j表示像元点，表示子空间上每个像元点经FrFT变换后的幅值，A_α表示变换系数，K_p(u,v)表示变换核函数，u表示时域轴，v表示FrFT变换后的频域轴，两者范围均属于(1,B)；α＝pπ/2表示旋转角，p为变换阶次，0<p<1；N_i是第i个子空间的原子个数；/>表示变换后每个像元点对应的马氏距离，μ为子空间均值，上角标T表示转置；

计算每一个子空间的原子的MD距离，选取距离最小的q个原子作为该子空间的背景原子Aⁱ＝[y₁,y₂,…,y_q]^T，其中Aⁱ表示第i个子空间的背景原子集合，q表示选取的背景原子个数，当该子空间内包含的像元数量少于q个，则直接舍弃；最终得到由q*M个背景原子构成的背景字典A＝A¹∪A²∪…∪A^M。

2)将高光谱数据通过FrFT变换映射至中间域，引入全变分正则化约束，建立基于全变分正则化的稀疏低秩模型作为异常检测模型，表示如下：

其中，表示经过FrFT变换后的在中间域上的变换矩阵，每个像元点为变换后的幅值；A表示背景字典，X表示背景字典的低秩系数，S为异常部分，即稀疏矩阵；‖·‖_*表示矩阵的核函数，λ、β、γ分别为各项的平衡系数，‖·‖_2,1表示l_2,1范数，H表示全变分正则化线性算子。

3)采用交替方向乘子法，把异常检测模型求解问题转换成求解多个子问题的最优解，得到稀疏部分；计算稀疏部分的范数得到最终的检测结果：

其中，表示每个像元点对应的异常度量值，[S^*]_:,i表示求解得到稀疏矩阵的行向量，其本质为求解稀疏矩阵的l_2,1范数；若/>大于设定的阈值，则该像素点属于异常点。

本发明的优点在于：从高光谱图像的低秩特性出发，建立全变分正则化约束的低秩稀疏模型，充分利用高光谱图像的光谱维和空间维信息，同时引入学习字典，将高光谱数据映射至多个子空间，并通过FrFT变换将异常信息从背景字典中剔除。本发明提出的算法对初始参数有较好鲁棒性，可以作为异常检测的有效手段，在地面异常小目标检测与识别方面有着重要的应用价值。

实际高光谱数据实验表明，基于FrFT的字典学习和全变分正则化有效提升了异常的检测率，抑制背景干扰。与另外传统的GRX、RPCA、FRFE、LSMAD、LRASR算法相比，本发明具有更好的探测效果，该方法的实际应用对于高光谱图像的异常检测问题的解决有着重要意义。

附图说明

图1算法示意图；

图2高光谱数据集。其中，(a)Urban数据集；(b)Pavia数据集；(c)Hyperion数据集；

图3不同p值对检测率的影响图；

图4不同参数对检测率的影响图；

图5Urban数据集不同算法的异常检测图；

图6Pavia数据集不同算法的异常检测图；

图7Hyperion数据集不同算法的异常检测图；

图8不同数据集的ROC曲线图。其中，(a)Urban(b)Pavia(c)Hyperion。

具体实施方式

下面结合附图，对本发明做进一步地说明。

具体内容介绍如下：

一、分数阶傅里叶变换(FrFT)

FrFT变换可看做信号在时频面内绕坐标轴原点逆时针旋转任意角度后构成的分数阶傅里叶域(即中间域)上的表示方法。它通过阶次的变化获得信号在时域、频域以及时频域的特征，相似信号在最优阶次域中具有很强的能量聚集性，对抑制噪声、目标检测等均具有很好的效果。高光谱图像中存在大量光谱曲线相似的背景信号，因此，本发明通过选取最优阶次，在中间域上提取高光谱图像的内在特征，增强异常目标与背景的光谱差异性，同时抑制噪声干扰。对于高光谱图像上任意一点y_i∈R^B×N，其FrFT变换表示如下：

K_p(u,v)＝A_αexp{jπ(u²cotα-2uv cscα+v²cotα)} (2)

其中，α＝pπ/2表示旋转角，0<p<1。K_p(u,v)表示核函数，u表示时域轴，v表示频域轴。表示像元点经FrFT变换后的幅值。

二、字典学习

背景字典的构建在异常检测中非常重要。字典原子能线性表示背景像元，而不能表示异常像元。纯净的背景字典将进一步提高异常检测的精度，降低虚警率。过去一些检测方法常用整个高光谱图像或者部分图像构建背景字典，这样会引入异常像元的污染，影响后续检测的精度。LRR模型引入多个子空间概念，每个子空间代表相似地物类型的像元，但其中也可能包含异常混合像元。考虑上述因素，本发明采用多个子空间的概念，通过k-means聚类算法，将光谱信息相似的像元映射到一个子空间中。其表达式如下：

其中，Yⁱ表示第i个子空间，M表示子空间个数，两两子空间之间无交集。A表示背景字典。在每个子空间内，FrFT变换抑制噪声和异常的干扰，增大背景原子和异常的可分性。把FrFT变换后像元的幅值作为新的反射值，并在子空间中间域上计算每一个原子与局部均值的马氏距离:

其中，N_i是第i个子空间的原子个数，μ为子空间均值。选取距离最小的q个原子作为该子空间的背景原子Aⁱ＝[y₁,y₂,…,y_q]^T。当该子空间内包含的像元数量少于q个，则直接跳过。最终，选取q*M个背景原子构成背景字典。

三、全变分正则化

空间特征是指当高光谱图像中两个像素在空间维度上相邻，那么它们对应的表示系数也相似。全变分正则化器(total variation，TV)在保留边缘信息和促进分段平滑性方面非常强大。高光谱图像通过FrFT变换后，背景像元与异常像元的空间差异性增强，即图像在背景区域是趋于平滑，在异常点邻域内存在梯度变化。因此，本发明引入全变分正则化描述中间域内高光谱图像的空间信息，以改进背景部分的估计。

由于空间相邻像素y_i和y_j可以表示为背景原子的线性组合，因此空间正则化约束可通过线性组合系数即低秩系数x_i和x_j来表示。对于高光谱图像上的任意点y_i∈R^B×N，TV正则化器表示如下：

TV(X)＝∑_{i,j}∈ε||x_i-x_j||₁ (7)

ε表示每个像素点空间邻域(水平和垂直方向)点的集合。引入水平方向和垂直方向的线性算子H_h和H_v，H_h表示水平方向相邻像素对应向量之间的差异，H_hX＝[d₁,d₂,…,d_N]，d_i＝|x_i,j-x_i+1,j|。H_v表示垂直方向相邻像素对应向量之间的差异，H_vX＝[v₁,v₂,…,v_N]，v_i＝|x_i,j-x_i,j+1|。公式(7)转化成：

其中，l_1,1范数定义为矩阵每列的l₁范数之和，

四、建模与优化

经过FrFT变换后，在中间域上引入全变分正则化项约束，本发明提出的模型表示如下：

其中，表示经过FrFT变换后的在中间域上的变换矩阵，每个像元点为变换后的幅值。A表示背景字典，X表示背景字典的低秩系数，S为异常部分。将全变分正则化项的范数形式代入：

min_X,S‖X‖_*+λ‖S‖_2,1+β‖HX‖_1,1+γ‖X‖₁ (10)

其中，模型第三项是全变分正则化项,第四项是低秩系数的稀疏性约束。该模型结合背景表示系数的低秩全局性和空间平滑性，以及稀疏部分的局部信息，可以更好地对异常进行检测。引入辅助变量：

上述目标函数是可分离的，优化问题可以通过ADMM解决，构造增广拉格朗日函数如下：

其中，D₁,D₂,D₃,D₄,D₅是拉格朗日乘子，τ是惩罚项系数。由于增广拉格朗日函数对每个变量是凸函数，可以保证迭代后解的全局最优性。该问题可分为多个变量的子问题。

固定变量(V₁,V₂,V₃,V₄,S)，优化变量X：

其中，k为迭代次数。

固定变量(X,V₂,V₃,V₄,S)，优化变量V₁：

其中，Θ为奇异值阈值算子(SVT)。

固定变量(x,V₁,V₃,V₄,S)，优化变量V₂：

固定变量(X,V₁,V₂,V₄,S)，优化变量V₃

其中，为软阈值收缩运算符。

固定变量(X,V₁,V₂,V₃,S)，优化变量V₄：

固定变量(X,V₁,V₂,V₃,V₄)，优化变量S：

其中，Ω表示l_2,1范数最小化算子。

迭代拉格朗日乘子D₁,D₂,D₃,D₄,D₅和惩罚项系数μ

τ^(k+1)＝min(ρτ^(k),τ_max) (24)

迭代要求：

五、算法流程框架

图1描述了所提算法的流程图，其主要步骤如下：

输入：高光谱图像Y∈R^B×N，参数λ>0，β>0，γ>0，变换参数p，聚类参数M,q

1.构建背景字典A；

2.FrFT变换把高光谱数据映射至中间域，得到

3.初始化变量：δ＝1e^-6，k＝0，k_max＝500，τ₀＝1e^-4，τ_max＝1e¹⁰，ρ＝1.5，其余矩阵变量初始化为零矩阵；

4.模型求解；

While模型不收敛或未达到迭代次数时do

4.1根据式(13)更新X

4.2根据式(14)，(15)，(16)，(17)分别更新V₁，V₂，V₃，V₄

4.3根据式(18)更新S

4.4根据式(19)—(24)分别更新D₁，D₂，D₃，D₄，D₅，τ

4.5计算是否达到收敛要求或迭代次数

End while

输出：稀疏矩阵S^*，异常检测图

上述算法得到稀疏矩阵S^*，得到最终探测结果：

它表示S^*异常矩阵最优解的第i列的l₂范数，若大于设定的阈值，则该像素点属于异常点。

实施例

为验证所提算法的有效性，本发明在3组真实高光谱数据集上进行实验。

一、实验数据

第一组数据由机载成像光谱仪(Hyperspectral Digital Image CollectionExperiment，HYDICE)采集，它的光谱分辨率是10nm，空间分辨率是1.5m，共有224个光谱波段，波长范围为400nm—2500nm。去除水吸收和低信噪比波段(1-4、76、87、101-111、136-153和198-210)，剩下162个波段作为实验数据。实验中使用尺寸为80*120的场景，包含21个异常目标，假彩色图像和目标真值图如图2(a)所示。

第二组数据是覆盖意大利北部帕维亚市中心区域(Pavia)，几何分辨率为1.3m。原始数据的空间大小为1096*715个像素，102个波段，波长范围为430—860nm，用于实验的场景覆盖面积为108*120个像素，总共有61个异常像元。数据集的假彩色图像和目标真值图如图2(b)所示。

第三组数据是由EO-1Hyperion传感器在奥卡万戈三角洲上空收集。光谱分辨率为10nm，空间分辨率为30m。去除水汽吸收和低信噪比波段(10-55、82-97、102-119、134-164和187-220)，使用其中145个波段。该场景的空间大小为1476*256像素，采用100*100像素大小的区域，共有32个异常像素。假彩色图像和目标真值图如图2(c)所示。

二、参数分析

本发明中，主要有以下参数需要预先设置：FrFT变换参数p，聚类参数M,q，模型参数λ、β和γ。在字典构建和模型优化中两者的FrFT变换参数相同，即映射至同一个中间域内。聚类参数中M表示子空间数，q表示每个子空间选择的背景原子数，分别取M＝15,q＝20。以AUC(Area Under Curve，AUC)指标评价算法的检测率。先对变换参数p进行实验，其余参数固定不变，分别记录在三个数据集上不同p值对检测率的影响。结果如图3所示。

当p值取0.6时，三个数据集的AUC值均达到最大，这说明当p取0.6时既保留数据原有的时域特征，又提取了傅里叶频域中的特征，因此取p＝0.6。而当p值大于0.6时，检测率下降，说明FrFT变换效果优于全局傅里叶变换，部分时域特征有必要被保留。

模型参数包括λ、β和γ，如图4所示，当p＝0.6时，固定两个变量，分析参数对检测结果的影响。由图4可知，当λ值和β值分别取0.5和0.05时，三个数据集的AUC值达到最大。对于γ参数，Urban数据集在0.1时达到最大，其余两个数据集在0.05时达到最大。Urban数据集和Pavia数据集对参数在某段区间上的变化不敏感，而Hyperion数据集检测精度随着参数的变化有较大的变化。结合上述分析，本发明模型参数取λ＝0.5，β＝0.05，γ＝0.05。

三、对比分析

为验证本发明所提方法的异常检测效果，与五种经典算法GRX,RPCA,FRFE,LSMAD,LRASR进行对比。其中，FRFE算法首次在RX算法中引入FrFT变换，GRX和RPCA是经典的RX系列算法，LSMAD和LRASR属于稀疏表示类算法，具有较好的检测效果。不同算法的定量评价指标采用ROC曲线(Receiver Operating Characterstic，ROC)和ROC曲线下的面积值AUC。ROC曲线越靠近左上方，检测效果越好。AUC值越高，表示检测效果越好。

图5、图6、图7分别表示不同方法在三个数据集上的检测效果。以真值图作为比较，当检测图与真值图越接近，定性认为该方法的检测效果越好。从图中可以看出本发明提出的方法(图5(f)、图6(f)和图7(f))不仅在一定程度上能抑制背景的影响，而且能准确检测出目标的异常区域以及异常点位置。LSMAD和LRASR算法都能较好显示异常点位置，但对于一些混合像元，在发明所提方法中混合像元与背景像元的可分性更高。此外，FrFE算法能较好抑制背景的干扰，但FrFE的检测率不如本发明所提方法。

为进一步定量分析所提方法的检测精度，图8表示不同算法的ROC曲线。表1显示不同算法在三个数据集上的AUC值。

在三个数据集中，本发明所提算法的AUC值最高，分别是0.9992,0.9990和0.9977。其中对于Urban数据集和Hyperion数据集，ROC曲线处于左上角的位置，这表明所提算法在抑制背景与检测精度上均优于其他算法。对于Pavia数据集，当虚警率为10^-4左右，LRASR的背景抑制效果较好。随着虚警率的增大，所提方法的检测效果优于LRASR算法。综合而言，所提方法融合了FrFT变换和全变分正则化两者的优势，具有较高的检测率和较低的虚警率。在不同的数据集上表现出较好的检测效果，算法具有一定的鲁棒性。

表1不同算法在三个数据集上的AUC值

本发明实施例所示的附图说明，可使本发明的目的、技术方案及优点介绍得更加清楚明白。应当说明，此处所描述的具体实施例仅仅用于解释本发明，并不用于限定本发明。凡在本发明提供的方法思路和原则之内所作的等同替换、改进等，均应包含在本发明的保护范围之内。

Claims

1.一种基于FrFT变换和全变分正则化的高光谱异常检测方法，其特征在于，步骤如下：

1)获得纯净的背景字典：

构造FrFT-RX算子，计算每一个子空间的原子的MD距离，选取距离最小的若干原子作为该子空间的背景原子，构建背景字典；在每个子空间内，FrFT-RX算子表示为：

K_p(u,v)＝A_αexp{jπ(u²cotα-2uvcscα+v²cotα)} (3)

2)将高光谱数据通过FrFT变换映射至中间域，引入全变分正则化约束，建立异常检测模型；所述中间域为信号在时频面内绕坐标轴原点逆时针旋转任意角度后构成的分数阶傅里叶域；

所述的步骤2)具体为：

对于高光谱图像上的任意点y_i，全变分正则化约束表示如下：

TV(X)＝∑_{i,j}∈ε||x_i-x_j||₁ (6)

其中，ε表示每个像素点空间邻域点的集合；x_i和x_j表示相邻像元点y_i和y_j的低秩系数，||·||₁表示l₁范数，TV(X)表示全变分正则化约束，X表示低秩系数矩阵；

引入水平方向和垂直方向的线性算子H_h和H_v，将公式(6)转化为：

其中，H_h表示水平方向相邻像素对应向量之间的差异，H_hX＝[d₁，d₂，...，d_N]，d_i＝|x_i，j-x_i+1，j|；H_v表示垂直方向相邻像素对应向量之间的差异，H_vX＝[v₁，v₂，...，v_N]，v_i＝|x_i，j-x_i，j+1|；d_i表示水平方向相邻像素低秩系数差异的列向量，x_i，j表示像元点对应的低秩系数，v_i表示垂直方向相邻像素低秩系数差异的列向量，|·|表示绝对值，||.||_1，1定义为矩阵每列的l₁范数之和；

经过FrFT变换后，在中间域上引入全变分正则化项约束，得到基于全变分正则化的稀疏低秩模型，表示如下：

其中，表示经过FrFT变换后的在中间域上的变换矩阵，每个像元点为变换后的幅值；A表示背景字典，X表示背景字典的低秩系数，S为异常部分，即稀疏矩阵；||·||_*表示矩阵的核函数，λ、β、γ分别为各项的平衡系数，||·||_2，1表示l_2，1范数，H表示全变分正则化线性算子；

3)采用交替方向乘子法，把异常检测模型求解问题转换成求解多个子问题的最优解，得到稀疏部分；计算稀疏部分的范数得到最终的检测结果。

2.根据权利要求1所述的基于FrFT变换和全变分正则化的高光谱异常检测方法，其特征在于，将选取的子空间的背景原子表示为Aⁱ＝[y₁，y₂，...，y_q]^T，其中Aⁱ表示第i个子空间的背景原子集合，q表示选取的背景原子个数，当该子空间内包含的像元数量少于q个，则直接舍弃；最终得到由q*M个背景原子构成的背景字典。

3.根据权利要求1所述的基于FrFT变换和全变分正则化的高光谱异常检测方法，其特征在于，所述的步骤3)具体为：

3.1)引入辅助变量，将基于全变分正则化的稀疏低秩模型转成以下等价问题：

其中，V₁，V₂，V₃，V₄是辅助变量；

3.2)构造增广拉格朗日函数如下：

其中，D₁，D₂，D₃，D₄，D₅是拉格朗日乘子，τ是惩罚项系数，表示矩阵的l₂范数；

3.3)将等价问题分为多个变量的子问题：

固定变量(V₁，V₂，V₃，V₄，S)，优化变量X：

其中，k为迭代次数，上角标T表示转置；

固定变量(X，V₂，V₃，V₄，S)，优化变量V₁：

其中，Θ为奇异值阈值算子；

固定变量(X，V₁，V₃，V₄，s)，优化变量V₂：

其中，I为单位矩阵；

固定变量(X，V₁，V₂，V₄，S)，优化变量V₃：

其中，为软阈值收缩运算符，τ^(k)表示第k次迭代得到的惩罚项系数；

固定变量(X，V₁，V₂，V₃，S)，优化变量V₄：

固定变量(X，V₁，V₂，V₃，V₄)，优化变量S：

其中，Ω表示l_2，1范数最小化算子；

迭代拉格朗日乘子D₁，D₂，D₃，D₄，D₅和惩罚项系数τ：

τ^(k+1)＝min(ρτ^(k)，τ_max) (23)

循环公式(12)-公式(23)，迭代计算参数X、V₁、V₂、V₃、V₄、S、D₁、D₂、D₃、D₄、D₅、τ；迭代要求：

其中，ρ表示惩罚项迭代系数，δ表示收敛值；

3.4)根据步骤3.3)获得的稀疏矩阵S^*计算探测结果：

其中，表示每个像元点对应的异常度量值，[S^*]_：，i表示求解得到稀疏矩阵的行向量，其本质为求解稀疏矩阵的l_2，1范数；若/>大于设定的阈值，则该像元点属于异常点。