不同于红外图像和多光谱图像，高光谱图像 (hyperspectral image，HSI) 具有数百个窄带，包含一段连续光谱范围内地物的辐射值，光谱分辨率可达nm级. 基于这一特性，高光谱图像在目标检测方面具有强大优势，被广泛应用于智能农业、矿产勘探和大气科学等多个领域^[1]. 作为高光谱图像处理的重要组成部分，异常检测旨在以无监督的方式检测出高维数据中的异常信息. 其中，异常是指在空间上稀疏分布的少量像素或者混合像元，具有和背景像素不同的光谱特征.

目前，Molero等^[2]提出各种用于高光谱异常检测的方法. 最经典的是RX算法（Reed-Xiaoli，RX），它假设背景像素遵循正态分布统计模型，通过计算测试像素与相邻背景像素之间的马氏距离来检测异常. 在实际应用中，由于场景地物类型复杂，以整幅图像的协方差矩阵估计背景统计特征不准确.

根据高光谱数据的异常稀疏性、背景低秩性的特点，一些基于非高斯分布的检测模型^[3-5]被提出. 鲁棒主成分分析（robust principal component analysis，RPCA）^[6]将图像分解为低秩的背景矩阵和稀疏的异常矩阵，通过矩阵分解重构残差. 低秩表示模型（low-rank representation，LRR）^[7]改进RPCA方法中单一子空间表示模型，将复杂背景用多个子空间线性混合表示. Xu等^[8]对低秩系数添加局部约束的稀疏项，提出低秩稀疏表示模型（low rank and sparse representation，LRASR）. 低秩稀疏表示的检测模型忽略了噪声和异常对构建背景字典过程中的污染问题. 只关注高光谱图像的光谱维信息，忽略图像的空间维信息. 这两点在一定程度上限制了最终的检测精度.

针对上述问题，本文提出基于FrFT（fractional Fourier transform，FrFT）变换与全变分正则化的异常检测模型. 通过FrFT变换提取光谱维的时频域特征，在每个子空间内构造FrFT-RX算子，剔除潜在的异常原子，得到纯净的背景字典. 在低秩稀疏模型中引入全变分正则化项，表示经过变换后中间域内背景部分的空间平滑性，综合利用图像的空谱信息. 采用交替方向乘子法（alternating direction method of multipliers，ADMM）^[9]将模型求解转换成多个子问题的求解问题，得到异常部分. 计算异常部分的 $ {{l}_{2,1}} $范数，得到异常检测结果.

基于背景的低秩特性和异常的稀疏特性，低秩稀疏表示模型把图像分为背景部分和异常部分. 对于任意高光谱图像 $\boldsymbol{Y}=[{\boldsymbol{y}}_1,{\boldsymbol{y}}_2,\cdots , {\boldsymbol{y}}_N]\in {\bf{R}}^{B\times N}$，表示模型如下：

(1) $ \boldsymbol{Y}=\boldsymbol{A}\boldsymbol{X}+\boldsymbol{S} . $

式中： $\boldsymbol{Y}\in {\bf{R}}^{B\times N}$为高光谱的二维矩阵形式，其中 $ B $为波段数, $ N $为像元数； $ \boldsymbol{A} $为背景字典矩阵， $\boldsymbol{A}= [{\boldsymbol{a}}_{1},{\boldsymbol{a}}_{2},\cdots ,{\boldsymbol{a}}_{{M}}]$，其中 $ M\mathrm{为}\mathrm{子}\mathrm{空}\mathrm{间}\mathrm{数}；\boldsymbol{X} $为背景字典的低秩系数矩阵， $\boldsymbol{X}=[{{\boldsymbol{x}}}_{1},{{\boldsymbol{x}}}_{2},\cdots ,{{\boldsymbol{x}}}_{N}]$； $ \boldsymbol{S} $为稀疏矩阵， $\boldsymbol{S}=[{{\boldsymbol{s}}}_{1},{{\boldsymbol{s}}}_{2},\cdots ,{{\boldsymbol{s}}}_{N}]$. 背景像元能被字典原子线性表示，对应系数存在低秩性；异常像元占比小，表现为稀疏性. 上述问题等价于求解如下问题：

(2) $ \left. \begin{array}{*{20}{c}} {{{\min }_{{\boldsymbol{X,S}}}}\;{{\left\| {\boldsymbol{X}} \right\|}_*} + \beta {{\left\| {\boldsymbol{X}} \right\|}_1} + \lambda {{\left\| {\boldsymbol{S}} \right\|}_{2,1}};}\\ {{\rm{s}}.{\rm{t}}.\;\;{\boldsymbol{Y}} = {\boldsymbol{AX}} +{\boldsymbol{ S}}.} \end{array} \right\}$

式中：‖·‖_*为矩阵的核范数，表示矩阵奇异值的和，用于约束矩阵的低秩性；‖·‖₁和‖·‖_2,1分别为矩阵的 ${l _1}$范数和 ${l _{2,1}}$范数； $ \;\beta $和 $ \lambda $分别为低秩部分和稀疏部分的平衡参数. 采用自适应惩罚的线性化交替方向法（linearized alternating direction method with adaptive penalty，LADMAP）^[10]，求解式（2）. 计算异常部分的 ${l _{2,1}}$范数,得到检测图.

针对低秩稀疏模型中背景字典易被污染、空间信息利用不足的问题，提出基于FrFT变换的全变分正则化模型，整体流程的示意图如图1所示. 通过k-means聚类映射至 $ M $个子空间， $Y^i$表示第 $ i $个子空间. 通过FrFT变换及马氏距离检测器,挑选 $ q $个原子构成背景字典. 在FrFT变换的中间域内，建立全变分正则化约束的低秩稀疏模型,计算得到模型最优解和异常检测结果.

FrFT变换可以看作信号在时频面内绕坐标轴原点逆时针旋转任意角度后构成的分数阶傅里叶域（即中间域）上的表示方法^[11]. 它通过阶次的变化，获得信号在时域、频域及时频域的特征；相似信号在最优阶次域中具有很强的能量聚集性，对抑制噪声、目标检测^[12]均具有很好的效果. 高光谱图像中存在大量光谱曲线相似的背景信号，因此，可以通过最优阶次的选取，在中间域上提取高光谱图像的内在特征，增强异常目标与背景的光谱差异性，抑制噪声干扰. 对于高光谱图像上的任意一点 ${{\boldsymbol{y}}}_{i}\in {\bf{R}}^{B\times N}$，FrFT变换表示如下：

(3) $ \tilde{{{\boldsymbol{y}}}_{i}}\left(v\right)=\frac{1}{B}\sum\nolimits _{u=1}^{B}{{\boldsymbol{y}}}_{i}\left(u\right){K}_{p}\left(u,v\right), $

(4) $ {K}_{p}\left(u,v\right)={A}_{\alpha }\mathrm{e}\mathrm{x}\mathrm{p}\;\left\{{\rm{j}}\text{π} \right({u}^{2}{\rm{cot}}\;\alpha -2uv{\rm{csc}}\;\alpha + {v}^{2}{\rm{cot}}\;\alpha \left)\right\}    , $

(5) $ {A}_{\alpha }=\frac{\mathrm{e}\mathrm{x}\mathrm{p}\;[-{\rm{j}}\text{π} {\rm{sgn}}\;({\rm{sin}}\;\alpha )/4+{\rm{j}}{\text{α}} /2]}{{\left|{\rm{sin}}\;\alpha \right|}^{1/2}} . $

$ \mathrm{式}\mathrm{中}：\alpha $为旋转角， $\alpha =p \text{π} /2；{A}_{\alpha }$表示变换系数； $ {K}_{p}\left(u,v\right) $表示核函数，其中0 $ < p < 1，u $表示时域轴， $ v $表示频域轴； ${{{\boldsymbol{y}}}_{i}}\left(v\right)$为像元点经FrFT变换后的幅值.

FrFT变换对高光谱图像中异常像元和背景像元的光谱曲线的变换效果如图2所示. 图中，R为反射值，A_m为光谱幅值.

FrFT变换的本质是对光谱曲线进行特征提取. 实验表明，当选择合适的 $ p $时，经过FrFT变换，背景和异常像元的差异性增大. 构建背景字典来表示变换后异常像元和背景像元的光谱差异性，全变分正则化项表示变换后两者的空间差异性.

背景字典的构建在异常检测中非常重要. 纯净的背景字典将进一步提高异常检测的精度，降低虚警率. 传统的基于表示的检测方法^[7]常用整个高光谱图像或者部分图像构建背景字典，这样会引入异常像元的污染，影响后续检测的精度. 尽管LRR模型引入多个子空间概念，每个子空间代表相似地物类型的像元，但可能包含异常混合像元. 考虑上述因素，在采用多个子空间的基础上，结合FrFT变换剔除异常像元的干扰. 采用k-means聚类算法，将光谱信息相似的像元映射到一个子空间中. FrFT变换的表达式如下：

(6) $ \left.\begin{array}{c}\boldsymbol{Y}={\cap }_{i=1,\cdots ,M}{\boldsymbol{Y}}^{i}.\\ {\boldsymbol{Y}}^{i}\cap {\boldsymbol{Y}}^{j}=\mathrm{\varnothing }\text{，}i\ne j.\end{array} \right\}$

式中： $ {\boldsymbol{Y}}^{i} $表示第 $ i $个子空间； $ M $为子空间个数，子空间之间无交集. $ \boldsymbol{A} $初始化为空集. 在每个子空间内，开展FrFT变换，将变换后像元的幅值作为新的反射值，在子空间中间域上计算每一个原子与局部均值的马氏距离:

(7) $ {\rm{MD}}\left(\tilde{{{\boldsymbol{y}}}_{i}}\right)={\left(\tilde{{\boldsymbol{y}}_{i}}-\boldsymbol{\mu }\right)}^{{\rm{T}}}{\boldsymbol{\varSigma }}^{-1}\left(\tilde{{\boldsymbol{y}}_{i}}-\boldsymbol{\mu }\right);\; j=\mathrm{1,2},\cdots ,{N}_{i} . $

式中： $ {N}_{i} $为第 $ i $个子空间的原子个数， $ \boldsymbol{\mu } $为子空间均值. 选取距离最小的 $ q $个原子，作为该子空间的背景原子 ${\boldsymbol{A}}^{i}={[{{\boldsymbol{y}}}_{1},{{\boldsymbol{y}}}_{2},\cdots ,{{\boldsymbol{y}}}_{q}]}^{{\rm{T}}}$. 当该子空间内包含的像元数量少于 $ q $个，则直接跳过. 选取 $ qM $个背景原子，构成背景字典 $ \boldsymbol{A}=\boldsymbol{A}\bigcup {\boldsymbol{A}}^{i} $.

在字典学习中，当 $ p $较小时，在每个子空间内经FrFT变换后保留较多的时域信息，挑选原子时易受到噪声异常信息的干扰. 当 $ p $较大时，保留较多的频域信息，部分背景原子特征被抑制，可能无法选择最合适的背景原子描述图像的背景特征. 选择合适的 $ p $值，剔除异常原子的干扰，得到完整的背景字典. 为了保证字典学习和模型优化中FrFT映射的中间域的一致性，两者的 $ p $相同.

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

由于空间相邻像素 ${{\boldsymbol{y}}}_{i}$和 ${{\boldsymbol{y}}}_{j}$可以表示为背景原子的线性组合，空间正则化约束可以通过线性组合系数，即低秩系数 ${{\boldsymbol{x}}}_{i}$和 ${{\boldsymbol{x}}}_{j}$来表示^[14]. 对于高光谱图像上的任意点 ${{\boldsymbol{y}}}_{i}\in {{\bf{R}}}^{B\times N}$， $ \mathrm{T}\mathrm{V} $ 正则化器表示如下：

(8) $ \mathrm{T}\mathrm{V}\left(\boldsymbol{X}\right)=\sum\nolimits _{\{i,j\}\in \epsilon }{\left\| {{{\boldsymbol{x}}}_{i}-{{\boldsymbol{x}}}_{j}} \right\|}_{1}. $

式中： $ \epsilon $为每个像素点空间邻域点的集合. 引入水平方向和垂直方向的线性算子 $ {\boldsymbol{H}}_{{\rm{h}}} $和 ${\boldsymbol{H}}_{{\rm{v}}}$^[15]. ${\boldsymbol{H}}_{{\rm{h}}}$表示水平方向相邻像素对应向量之间的差异， ${\boldsymbol{H}}_{{\rm{h}}}\boldsymbol{X}= [{\boldsymbol{d}}_{1},{\boldsymbol{d}}_{2},\cdots ,{\boldsymbol{d}}_{N}]$， ${\boldsymbol{d}}_{i}={\boldsymbol{x}}_{i,j}-{\boldsymbol{x}}_{i+1,j}$. ${\boldsymbol{H}}_{{\rm{v}}}$表示垂直方向相邻像素对应向量之间的差异， ${\boldsymbol{H}}_{{\rm{v}}}\boldsymbol{X}=[{\boldsymbol{v}}_{1},{\boldsymbol{v}}_{2},\cdots ,{\boldsymbol{v}}_{N}]$， ${\boldsymbol{v}}_{i}={\boldsymbol{x}}_{i,j}-{\boldsymbol{x}}_{i,j+1}$. 式（8）转化为

(9) $ \mathrm{T}\mathrm{V}\left(\boldsymbol{X}\right)={\left\| {\left[\begin{array}{c}{\boldsymbol{H}}_{{\rm{h}}}\boldsymbol{X}\\ {\boldsymbol{H}}_{{\rm{v}}}\boldsymbol{X}\end{array}\right]} \right\|}_{\mathrm{1,1}}={\left\| {\boldsymbol{H}\boldsymbol{X}} \right\|}_{\mathrm{1,1}}. $

式中： ${{l }_{1,1}}$范数定义为矩阵每列的 ${{l }_{1}}$范数之和， ${\left\| {\boldsymbol{X}} \right\|}_{\mathrm{1,1}}=\displaystyle\sum\nolimits _{i=1}^{n}{\left\| {{\boldsymbol{x}}_{i}} \right\|}_{1}$.

经过FrFT变换后，在中间域上引入全变分正则化项约束，模型表示为

(10) $ \left. \begin{array}{*{20}{c}} {{\rm{min}}_{{\boldsymbol{X,S}}}}\;{\left\| {\boldsymbol{X}} \right\|}_{*}+\lambda {\left\| {\boldsymbol{S}} \right\|}_{\mathrm{2,1}}+{\beta \mathrm{T}\mathrm{V}\left(\boldsymbol{X}\right)+\gamma \left\| {\boldsymbol{X}} \right\|}_{1},\\ {\rm{s.t.}} \;\tilde{\boldsymbol{Y}}=\boldsymbol{A}\boldsymbol{X}+\boldsymbol{S}.\end{array}\right\} $

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

(11) $ {{\rm{min}}_{{\boldsymbol{X,S}}}}\;{\left\|\boldsymbol{X}\right\|}_{*}+\lambda {\left\|\boldsymbol{S}\right\|}_{\mathrm{2,1}}+{\beta {\left\|\boldsymbol{H}\boldsymbol{X}\right\|}_{\mathrm{1,1}}+\gamma \left\|\boldsymbol{X}\right\|}_{1}. $

式中：模型第3项是全变分正则化项,第4项是低秩系数的稀疏性约束. 该模型结合背景表示系数的低秩全局性和空间平滑性以及稀疏部分的局部信息，可以更好地对异常进行检测. 式（11）对每个变量都是凸函数，目标函数是可分离的，因此可以通过ADMM解决，迭代解是收敛的^{[10, 16]}. 引入辅助变量，将式（11）转化成以下等价问题：

(12) $ \left. \begin{array}{*{20}{c}} \underset{{\boldsymbol{V}}_{1},{\boldsymbol{V}}_{3},{\boldsymbol{V}}_{4},\boldsymbol{S}}{\mathrm{min}}\;{\left\| {{\boldsymbol{V}}_{1}} \right\|}_{*}+\lambda {\left\| {\boldsymbol{S}} \right\|}_{\mathrm{2,1}}+{\beta {\left\| {{\boldsymbol{V}}_{3}} \right\|}_{\mathrm{1,1}}+\gamma \left\| {{\boldsymbol{V}}_{4}} \right\|}_{1};\\ {\rm{s.t.}}\;\tilde{\boldsymbol{Y}}=\boldsymbol{A}\boldsymbol{X}+\boldsymbol{S},{\boldsymbol{V}}_{1}=\boldsymbol{X},{\boldsymbol{V}}_{2}=\boldsymbol{X},{\boldsymbol{V}}_{3}=\boldsymbol{H}{\boldsymbol{V}}_{2},{\boldsymbol{V}}_{4}=\boldsymbol{X}.\end{array} \right\}$

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

(13) $ \begin{split} &\mathcal{L}\left({\boldsymbol{V}}_{1},{\boldsymbol{V}}_{2},{\boldsymbol{V}}_{3},{\boldsymbol{V}}_{4},\boldsymbol{X},\boldsymbol{S},{\boldsymbol{D}}_{1},{\boldsymbol{D}}_{2},{\boldsymbol{D}}_{3},{\boldsymbol{D}}_{4},{\boldsymbol{D}}_{5},\tau \right) =\\ &{\left\| {{\boldsymbol{V}}_{1}} \right\|}_{*}+\lambda {\left\| {\boldsymbol{S}} \right\|}_{\mathrm{2,1}}+{\beta {\left\| {{\boldsymbol{V}}_{3}} \right\|}_{\mathrm{1,1}}+\gamma \left\| {{\boldsymbol{V}}_{4}} \right\|}_{1} +\\ &\frac{\tau }{2}\left({\left\| {\tilde{\boldsymbol{Y}}-\boldsymbol{A}\boldsymbol{X}-\boldsymbol{S}-{\boldsymbol{D}}_{1}} \right\|}_{{\rm{F}}}^{2}+{\left\| {{\boldsymbol{V}}_{1}-\boldsymbol{X}-{\boldsymbol{D}}_{2}} \right\|}_{{\rm{F}}}^{2} \right. +\\ &\left.{\left\| {{\boldsymbol{V}}_{2}-\boldsymbol{X}-{\boldsymbol{D}}_{3}} \right\|}_{{\rm{F}}}^{2}+{\left\| {{\boldsymbol{V}}_{3}-\boldsymbol{H}{\boldsymbol{V}}_{2}-{\boldsymbol{D}}_{4}} \right\|}_{{\rm{F}}}^{2} + {\left\| {{\boldsymbol{V}}_{4}-\boldsymbol{X}-{\boldsymbol{D}}_{5}} \right\|}_{{\rm{F}}}^{2}\right). \end{split} $

式中： ${\boldsymbol{D}}_{1}、{\boldsymbol{D}}_{2}、{\boldsymbol{D}}_{3}、{\boldsymbol{D}}_{4}、{\boldsymbol{D}}_{5}$为拉格朗日乘子， $ \tau $为惩罚项系数. ADMM算法通过固定其他变量，迭代优化每一个变量的目标函数，由于增广拉格朗日函数对每个变量是凸函数，可以保证迭代后解的全局最优性. 该问题可以分为以下子问题.

1）固定变量（ ${{\boldsymbol{V}}_{1}、\boldsymbol{V}}_{2}、{\boldsymbol{V}}_{3}、{\boldsymbol{V}}_{4}、\boldsymbol{S}$），优化变量 $ \boldsymbol{X} $.

(14) $ \begin{split} {\boldsymbol{X}}^{(k+1)}=\;&{\rm{arg}}\underset{\boldsymbol{X}}{\mathrm{m}\mathrm{i}\mathrm{n}}\frac{{\tau }^{\left(k\right)}}{2}\left({\left\| {\tilde{\boldsymbol{Y}}-\boldsymbol{A}\boldsymbol{X}-{\boldsymbol{S}}^{\left(k\right)}-{\boldsymbol{D}}_{1}^{\left(k\right)}} \right\|}_{{\rm{F}}}^{2} +\right.\\ &{\left\| {{\boldsymbol{V}}_{1}^{\left(k\right)}-\boldsymbol{X}-{\boldsymbol{D}}_{2}^{\left(k\right)}} \right\|}_{{\rm{F}}}^{2}+{\left\| {{\boldsymbol{V}}_{2}^{\left(k\right)}-\boldsymbol{X}-{\boldsymbol{D}}_{3}^{\left(k\right)}} \right\|}_{{\rm{F}}}^{2} +\\ &\left.{\left\| {{\boldsymbol{V}}_{4}^{\left(k\right)}-\boldsymbol{X}-{\boldsymbol{D}}_{5}^{\left(k\right)}} \right\|}_{{\rm{F}}}^{2}\right) . \end{split} $

式中： $ k $为迭代次数. 式（14）中的目标函数为多个 ${{l }_{2}}$范数的和，令目标函数微分右式为零：

(15) $ \begin{split} &{\boldsymbol{X}}^{(k+1)}={({\boldsymbol{A}}^{{\rm{T}}}\boldsymbol{A}+3\boldsymbol{I})}^{-1}\left[{\boldsymbol{A}}^{{\rm{T}}}\left(\boldsymbol{Y}-{\boldsymbol{S}}^{\left(k\right)}-{\boldsymbol{D}}_{1}^{\left(k\right)}\right) +\right.\\ &\left.\left({\boldsymbol{V}}_{1}^{\left(k\right)}-{\boldsymbol{D}}_{2}^{\left(k\right)}\right)+\left({\boldsymbol{V}}_{2}^{\left(k\right)}-{\boldsymbol{D}}_{3}^{\left(k\right)}\right)+\left({\boldsymbol{V}}_{4}^{\left(k\right)}-{\boldsymbol{D}}_{5}^{\left(k\right)}\right)\right]. \end{split} $

2）固定变量（ $\boldsymbol{X}、{\boldsymbol{V}}_{2}、{\boldsymbol{V}}_{3}、{\boldsymbol{V}}_{4}、\boldsymbol{S}$），优化变量 $ {\boldsymbol{V}}_{1} $.

(16) $ {\boldsymbol{V}}_{1}^{(k+1)}={\rm{arg}}\underset{{\boldsymbol{V}}_{1}}{{\rm{min}}}\;{\left\| {{\boldsymbol{V}}_{1}} \right\|}_{*} +\frac{{\tau }^{\left(k\right)}}{2}{\left\| {{\boldsymbol{V}}_{1}-{\boldsymbol{X}}^{(k+1)}-{\boldsymbol{D}}_{2}^{\left(k\right)}} \right\|}_{{\rm{F}}}^{2} . $

上述问题为核范数和 $ {{l}}_{2} $范数的优化问题，引入奇异值算子（singular value thresholding，SVT）^[17]解决，其中 $\varTheta$为奇异值阈值算子（SVT）：

(17) $ {\boldsymbol{V}}_{1}^{(k+1)}={{\varTheta }}_{1/{\tau }^{\left(k\right)}}({\boldsymbol{X}}^{\left(k+1\right)}+{\boldsymbol{D}}_{2}^{\left(k\right)}). $

3）固定变量（ $ {\boldsymbol{X}、\boldsymbol{V}}_{1}、{\boldsymbol{V}}_{3}、{\boldsymbol{V}}_{4}、\boldsymbol{S} $），优化变量 $ {\boldsymbol{V}}_{2} $.

(18) $ \begin{split} {\boldsymbol{V}}_{2}^{(k+1)}=\; &{\rm{arg}}\underset{{\boldsymbol{V}}_{2}}{\mathrm{min}}\frac{{\tau }^{\left(k\right)}}{2}\left(\left\| {{\boldsymbol{V}}_{2}-{\boldsymbol{X}}^{(k+1)}-{\boldsymbol{D}}_{3}^{\left(k\right)}} \right\|_{{\rm{F}}}^{2} +\right.\\ &\left.{\left\| {{\boldsymbol{V}}_{3}^{\left(k\right)}-\boldsymbol{H}{\boldsymbol{V}}_{2}-{\boldsymbol{D}}_{4}^{\left(k\right)}} \right\|}_{{\rm{F}}}^{2}\right). \end{split} $

对式(18)中的目标函数求导，等式右边取零，得到

(19) $ {\boldsymbol{V}}_{2}^{(k+1)}={({\boldsymbol{H}}^{{\rm{T}}}\boldsymbol{H}+\boldsymbol{I})}^{-1} \left[{\boldsymbol{X}}^{\left(k+1\right)}+{\boldsymbol{D}}_{3}^{\left(k\right)}+\boldsymbol{H}^{{\rm{T}}}\left({\boldsymbol{V}}_{3}^{\left(k\right)}+{\boldsymbol{D}}_{4}^{\left(k\right)}\right)\right]. $

式（19）可以通过离散傅里叶变换对角化^[15] ，进行有效计算.

4）固定变量（ ${\boldsymbol{X}、\boldsymbol{V}}_{1}、{\boldsymbol{V}}_{2}、{\boldsymbol{V}}_{4}、\boldsymbol{S}$），优化变量V₃.

(20) $ {\boldsymbol{V}}_{3}^{(k+1)}={\rm{arg}}\underset{{\boldsymbol{V}}_{3}}{\mathrm{m}\mathrm{i}\mathrm{n}}\;\;\beta {\left\| {{\boldsymbol{V}}_{3}} \right\|}_{\mathrm{1,1}} +\frac{{\tau }^{\left(k\right)}}{2}{\left\| {{\boldsymbol{V}}_{3}-\boldsymbol{H}{\boldsymbol{V}}_{2}^{(k+1)}-{\boldsymbol{D}}_{4}^{\left(k\right)}} \right\|}_{{\rm{F}}}^{2}. $

引入软阈值收缩运算符 $ \mathcal{S} $（soft-thresholding shrinkage operator）^[18]求解如下：

(21) $ {\boldsymbol{V}}_{3}^{(k+1)}={\mathcal{S}}_{\beta /{\tau }^{\left(k\right)}}(\boldsymbol{H}{\boldsymbol{V}}_{2}^{\left(k+1\right)}+{\boldsymbol{D}}_{4}^{\left(k\right)}) . $

5）固定变量（ ${\boldsymbol{X}、\boldsymbol{V}}_{1}、{\boldsymbol{V}}_{2}、{\boldsymbol{V}}_{3}、\boldsymbol{S}$），优化变量 $ {\boldsymbol{V}}_{4} $.

(22) $ {\boldsymbol{V}}_{4}^{(k+1)}={\rm{arg}}\underset{{\boldsymbol{V}}_{4}}{\mathrm{m}\mathrm{i}\mathrm{n}}\;\;\gamma {\left\| {{\boldsymbol{V}}_{4}} \right\|}_{1} +\frac{{\tau }^{\left(k\right)}}{2}{\left\| {{\boldsymbol{V}}_{4}-{\boldsymbol{X}}^{(k+1)}-{\boldsymbol{D}}_{5}^{\left(k\right)}} \right\|}_{{\rm{F}}}^{2}, $

(23) $ {\boldsymbol{V}}_{4}^{(k+1)}={\mathcal{S}}_{\gamma /{\tau }^{\left(k\right)}}({\boldsymbol{X}}^{(k+1)}+{\boldsymbol{D}}_{5}^{\left(k\right)}) . $

6）固定变量（ ${\boldsymbol{X}、\boldsymbol{V}}_{1}、{\boldsymbol{V}}_{2}、{\boldsymbol{V}}_{3}、{\boldsymbol{V}}_{4}$），优化变量 $ \boldsymbol{S} $：

(24) $ {\boldsymbol{S}}^{\left(k+1\right)}={\rm{arg}}\underset{\boldsymbol{S}}{\mathrm{m}\mathrm{i}\mathrm{n}}\;\lambda {\left\| {\boldsymbol{S}} \right\|}_{\mathrm{2,1}} +\frac{{\tau }^{\left(k\right)}}{2}{\left\| {‖\tilde{\boldsymbol{Y}}-\boldsymbol{A}{\boldsymbol{X}}^{\left(k+1\right)}-\boldsymbol{S}-{\boldsymbol{D}}_{1}^{\left(k\right)}} \right\|}_{{\rm{F}}}^{2} .$

引入 ${{l }_{2,1}}$范数最小化算子^[7]，求解如下：

(25) $ {\boldsymbol{S}}^{\left(k+1\right)}={\mathrm{\varOmega }}_{\lambda /{\tau }^{\left(k\right)}}(\tilde{\boldsymbol{Y}}-\boldsymbol{A}{\boldsymbol{X}}^{\left(k+1\right)}-{\boldsymbol{D}}_{1}^{\left(k\right)}). $

式中： $\mathrm{\varOmega }$为 ${{l }_{2,1}}$范数最小化算子.

7）迭代拉格朗日乘子 $ {\boldsymbol{D}}_{1}{\mathbf{、}\boldsymbol{D}}_{2}{\mathbf{、}\boldsymbol{D}}_{3}{\mathbf{、}\boldsymbol{D}}_{4}{\mathbf{、}\boldsymbol{D}}_{5} $和惩罚项系数 $ \tau $：

(26) $ {\boldsymbol{D}}_{1}^{\left(k+1\right)}={\boldsymbol{D}}_{1}^{\left(k\right)}-(\tilde{\boldsymbol{Y}}-\boldsymbol{A}{\boldsymbol{X}}^{\left(k+1\right)}-{\boldsymbol{S}}^{\left(k+1\right)}). $

(27) $ {\boldsymbol{D}}_{2}^{\left(k+1\right)}={\boldsymbol{D}}_{2}^{\left(k\right)}-({\boldsymbol{V}}_{1}^{(k+1)}-{\boldsymbol{X}}^{(k+1)}) . $

(28) $ {\boldsymbol{D}}_{3}^{\left(k+1\right)}={\boldsymbol{D}}_{3}^{\left(k\right)}-({\boldsymbol{V}}_{2}^{(k+1)}-{\boldsymbol{X}}^{(k+1)}). $

(29) $ {\boldsymbol{D}}_{4}^{\left(k+1\right)}={\boldsymbol{D}}_{4}^{\left(k\right)}-({\boldsymbol{V}}_{3}^{(k+1)}-\boldsymbol{H}{\boldsymbol{V}}_{2}^{(k+1)}). $

(30) $ {\boldsymbol{D}}_{5}^{\left(k+1\right)}={\boldsymbol{D}}_{5}^{\left(k\right)}-({\boldsymbol{V}}_{4}^{(k+1)}-{\boldsymbol{X}}^{(k+1)}). $

(31) $ {\tau }^{(k+1)}=\mathrm{m}\mathrm{i}\mathrm{n}\;\{\rho {\tau }^{\left(k\right)},{\tau }_{{\rm{max}}}\}. $

当迭代条件满足

(32) $ \begin{split} &{\left\| {\tilde{\boldsymbol{Y}}-\boldsymbol{A}\boldsymbol{X}-\boldsymbol{S}} \right\|}_{{\rm{F}}}^{2}+{\left\| {{\boldsymbol{V}}_{1}-\boldsymbol{X}} \right\|}_{{\rm{F}}}^{2}+{\left\| {{\boldsymbol{V}}_{2}-\boldsymbol{X}} \right\|}_{{\rm{F}}}^{2} +\\ &\qquad{\left\| {{\boldsymbol{V}}_{3}-\boldsymbol{H}{\boldsymbol{V}}_{2}} \right\|}_{{\rm{F}}}^{2}+{\left\| {{\boldsymbol{V}}_{4}-\boldsymbol{X}} \right\|}_{{\rm{F}}}^{2}\leqslant \delta \text{，} \end{split} $

或最大迭代次数达到要求时，迭代结束，得到低秩矩阵和稀疏矩阵. 由于FrFT变换具有噪声抑制的效果，图像经过FrFT变换后，在中间域上引入全变分正则化约束进行迭代，得到的稀疏矩阵剔除噪声干扰，可以表示为异常部分.

当目标函数达到最优解时，异常部分可以求解如下：

(34) $ T\left(\tilde{{{\boldsymbol{y}}}_{i}}\right)={\left\| {{\left[{\boldsymbol{S}}^{\boldsymbol{*}}\right]}_{:,i}} \right\|}_{2}. $

式中： $\tilde{{\boldsymbol{y}}}_i$为经过分数阶傅里叶变换后的图像矩阵上的点； $T (\tilde{{\boldsymbol{y}}}_i)$为每一点对应的异常值，若异常值大于设定阈值，则该像素属于异常点.

式（33）表示 $ {\boldsymbol{S}}^{*} $异常矩阵最优解的第 $ i $列的 ${{l }_{2}}$范数，本质上是异常矩阵的 $ {{l }_{2,1}} $范数. 若 $T\left(\tilde{{{\boldsymbol{y}}}_{i}}\right)$大于设定的阈值，则该像素点属于异常点.

通过论述和公式推导，得到基于FrFT变换和全变分正则化的异常检测算法流程如下.

算法1：基于FrFT变换和全变分正则化异常检测算法

输入：高光谱图像 $\boldsymbol{Y}\in {{\bf{R}}}^{B\times N}$，参数 $ \lambda > 0 $， $ \;\beta > 0 $， $ \gamma > 0 $，变换参数 $ p $，聚类类别数 $ M $, 每个聚类的原子数 $ q $.

1）利用2.2节的内容，构建背景字典 $ \boldsymbol{A} $.

2） 利用FrFT变换，将高光谱数据映射至中间域，得到 $ \tilde{\boldsymbol{Y}}\in {\mathbf{R}}^{B\times N} $.

3）初始化变量： $\delta =10^{-6}$， $ k=0 $， $ {k}_{\mathrm{m}\mathrm{a}\mathrm{x}}=500 $， ${\tau }_{0}=10^{-4}$， ${\tau }_{\mathrm{m}\mathrm{a}\mathrm{x}}=10^{10}$， $ \rho =1.5 $；初始化 $ \boldsymbol{X} $、 $ {\boldsymbol{V}}_{1} $、 $ {\boldsymbol{V}}_{2} $、 $ {\boldsymbol{V}}_{3} $、 $ {\boldsymbol{V}}_{4} $、 $ \boldsymbol{S} $、 $ {\boldsymbol{D}}_{1} $、 $ {\boldsymbol{D}}_{2} $、 $ {\boldsymbol{D}}_{3} $、 $ {\boldsymbol{D}}_{4} $、 $ {\boldsymbol{D}}_{5} $为零矩阵.

4）模型求解. 当模型不收敛或未达到迭代次数时，执行以下步骤.

a）根据式（15）更新 $ \boldsymbol{X} $.

b) 根据式（17）、（19）、（21）、（23），分别更新 $ {\boldsymbol{V}}_{1} $、 $ {\boldsymbol{V}}_{2} $、 $ {\boldsymbol{V}}_{3} $、 $ {\boldsymbol{V}}_{4} $.

c）根据式（25）更新 $ \boldsymbol{S} $.

d) 根据式（26）~（31），分别更新 $ {\boldsymbol{D}}_{1} $、 $ {\boldsymbol{D}}_{2} $、 $ {\boldsymbol{D}}_{3} $、 $ {\boldsymbol{D}}_{4} $、 $ {\boldsymbol{D}}_{5} $、 $ \tau $.

e) 计算是否满足收敛要求或达到最大迭代次数.

输出：稀疏矩阵 $ {\boldsymbol{S}}^{\boldsymbol{*}} $、异常检测图.

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

第1组数据由机载成像光谱仪（hyperspectral digital image collection experiment，HYDICE）采集. 它的光谱分辨率是10 nm，空间分辨率是1.5 m，共有224个光谱波段，波长为400~2 500 nm. 去除水汽吸收和低信噪比波段（1~4、76、87、101~111、136~153和198~210），剩下162个波段作为实验数据. 实验中使用尺寸为80×120像素的场景，包含21个异常目标，假彩色图像和目标真值图如图3所示.

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

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

本文方法的参数主要有FrFT变换参数 $ p $，聚类数 $ M $,每类原子数 $ q $，模型参数 $ \lambda $、 $\; \beta $和 $ \gamma $. 在字典构建和模型优化中，两者的FrFT变换参数相同. 在聚类参数中 $ ，M $为子空间数， $ q $为每个子空间选择的背景原子数. 文献[12]中的分析分别取 $ M=15 $， $ q=20 $. 以AUC（area under curve，AUC）指标评价算法的检测率. 改变 $ p $，其余参数固定不变，分别记录3个数据集上不同 $ p $对检测率的影响，结果如图6所示. 图6描述了在不同数据集上变换参数取值对AUC的影响.

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

模型参数包括 $ \lambda $、 $\; \beta $和 $ \gamma $，如图7所示，当 $ p=0.6 $时，固定2个变量，分析参数对检测结果的影响. 从图7可知，当 $ \lambda $和 $ \beta $分别取0.5和0.05时，3个数据集的AUC达到最大. 对于 $ \gamma $参数，Urban数据集的AUC在 $ \gamma $ = 0.1时达到最大，其余2个数据集的AUC $ \mathrm{在}\gamma $ = 0.05时达到最大. Urban数据集和Pavia数据集对参数在某段区间上的变化不敏感， Hyperion数据集检测精度随着参数的变化有较大的变化. 结合上述分析，实验取 $ \lambda =0.5 $， $\; \beta =0.05 $， $ \gamma =0.05 $.

为了验证本文方法的异常检测效果，与5种经典算法GRX、RPCA、FRFE、LSMAD、LRASR进行对比. 其中，基于分数傅里叶熵的高光谱异常检测算法（fractional Fourier transform，FRFE）^[19]在RX算法中引入FrFT变换，GRX是经典的RX系列算法，RPCA、LSMAD和LRASR属于稀疏表示类算法，具有较好的检测效果. 不同算法的定量评价指标采用ROC曲线（receiver operating characterstic，ROC）和ROC曲线下的面积AUC. 其中，ROC曲线越靠近左上方，表示检测效果越好. AUC越高，表示检测效果越好.

不同方法在3个数据集上的检测效果分别如图8~10所示. 以真值图作为比较，检测图与真值图越接近，定性认为该方法的检测效果越好. 从图8~10可以看出，本文方法不仅在一定程度上能够抑制背景的影响，而且能够准确检测出目标的异常区域及异常点位置. LSMAD和LRASR算法都能够较好地显示异常点位置，但对于一些混合像元，在利用本文方法得到的异常检测图中，混合像元与背景像元的可分性更高. 此外，FRFE算法能够较好地抑制背景的干扰，但FRFE的检测率不如本文方法.

为了进一步分析本文方法的检测精度，如图11所示为不同算法的ROC曲线. 图中，P_d为检测概率，P_f为虚警概率. 不同算法在3个数据集上的AUC如表1所示.

在3个数据集中，本文方法的AUC最高，分别是0.999 2、0.999 0和0.997 7. 对于Urban数据集和Hyperion数据集，ROC曲线处于左上角的位置，这表明本文方法在抑制背景与检测精度上均优于其他算法. 对于Pavia数据集，当虚警率约为 $ {10}^{-4} $时，LRASR的背景抑制效果较好. 随着虚警率的增大，本文方法的检测效果优于LRASR算法.

为了进一步表示算法的背景抑制效果和目标检测能力，引入 $ {(P}_{\mathrm{f}},\eta ) $和 $ {(P}_{\mathrm{d}},\eta ) $^[20]，其中 $ \eta $为阈值. 引入上述2个参数值实际上是对ROC曲线的细化，它们表示不同的 $ \eta $对应的虚警概率和检测概率. 以 $ \eta $为横坐标， $ {P}_{\mathrm{f}} $、 $ {P}_{\mathrm{d}} $为纵坐标，计算2条曲线下的面积，得到 $ {(P}_{\mathrm{f}},\eta ) $和 $ {(P}_{\mathrm{d}},\eta ) $. $ {(P}_{\mathrm{f}},\eta ) $表示算法的背景抑制效果，数值越小，背景抑制能力越强； $ {(P}_{\mathrm{d}},\eta ) $表示算法的目标检测效果，数值越大，检测能力越好. 如表2、3所示为6种算法的背景抑制效果和检测能力. 在3个数据集中，本文方法的 $ {(P}_{\mathrm{d}},\eta ) $最大，这表明本文方法具有较强的检测能力. 对于 $ {(P}_{\mathrm{f}},\eta ) $，在Urban数据集中所提方法最小，在其余2个数据集中数值略大于LSMAD算法. 虽然本文方法的 $ {(P}_{\mathrm{f}},\eta ) $都低于LRASR算法，但在一定程度上改善了LRASR算法高检测率带来的高虚警率问题. 这表示算法的背景抑制能力可以进一步提高.

综合而言，本文方法融合了FrFT变换和全变分正则化两者的优势，具有较高的检测率和较低的虚警率. 如表4所示为不同算法的运行时间t_r. 由于利用图像的空间信息，算法的运行时间较长.

针对传统异常检测算法背景字典污染和空间信息利用不足而导致检测精度低的问题，本文提出基于FrFT变换和全变分正则化的异常检测算法. 利用FrFT变换能够有效地提取光谱维的时频域特征，在构造背景字典时有利于剔除潜在的异常原子，得到较纯净的背景字典. 将高光谱数据通过FrFT变换映射至中间域，在低秩稀疏模型中引入全变分正则化约束，以表示中间域内背景部分的空间平滑性，综合利用了图像的空谱信息. 通过交替方向乘子法，将模型求解转换成求解多个子问题的最优解，得到稀疏部分. 计算稀疏部分的范数，得到最终的检测结果. 通过定性和定量分析，本文算法在提高检测率和降低虚警率方面均取得了一定的效果. 由于联合空谱信息，算法的运行时间较长，在未来工作中，将进一步优化算法，在更多数据集上进行测试，提高算法的鲁棒性.