一种基于新型多元偏态t分布模型的高光谱图像异常检测方法

1.本发明属于遥感图像处理技术邻域，具体涉及一种基于自编码器和新型多元偏态t分布模型的高光谱图像异常检测方法。


背景技术：

2.高光谱图像在传统遥感图像的基础上加入了光谱维，是包含着目标空间特征和光谱信息的三维数据结构。高光谱图像检测技术一般可以分为目标检测和异常检测。由于不需要提前知晓目标物体光谱信息便可以检测出异常目标的特点，高光谱异常检测被广泛应用于农业、气象、军事等诸多领域。高光谱异常检测的研究变得愈发重要。
3.基于自编码器的高光谱图像异常检测主要包括图像重建和异常检出两个步骤。图像重建通过寻找高光谱图像中背景信号的方式来重新建立图像，异常检出则依据异常评分机制评估各个像元是否属于背景，不属于背景的即为异常像元。异常评分机制的好坏直接影响到最终的异常检测结果。同时，由于目前基于自编码器的异常检测算法大都是基于偏差的方法，将重建误差认定为服从l1、l2范数形式的拉普拉斯分布和正态分布，实际应用过程中效果并不好，缺乏鲁棒性。


技术实现要素：

4.为解决上述现有技术中存在的缺陷，本发明的目的在于一种基于新型多元偏态t分布模型的高光谱图像异常检测方法，基于自编码器和新型多元偏态t分布模型，实现高光谱图像异常检测。
5.为实现上述目的，本发明采取的技术方案如下：
6.一种基于高光谱图像的异常检测方法，包括如下步骤：
7.步骤1，对原始高光谱图像数据x∈r
mn
×b，使用rx算法进行预处理，得到其背景图像部分；其中m、n表示图像的空间维度，mn即像元数，b表示光谱维度；
8.步骤2，将背景图像部分作为基于光谱角度余弦的栈式降噪自编码器网络的输入，对网络进行训练以获得训练好的栈式降噪自编码器；
9.步骤3，将原始高光谱图像作为训练好的栈式降噪自编码器的输入，选择输出层结果输出，得到重建的高光谱图像数据x
′
∈r
mn
×b；
10.步骤4，计算重建的高光谱图像中第i个像元的重建光谱向量xi′
与原始光谱向量xi的差值di，从而得到光谱向量差值图像，对光谱向量差值图像使用局部标准偏差过滤器，从而通过滤波获得异常更为突出的重建误差ri∈rb，i＝1,2,
…
,mn；
11.步骤5，将重建误差输入到mvskt分布模型中，计算得到异常评分ad(ri)，根据异常评分ad(ri)判断待测像元类别，输出异常检测结果。
12.所述步骤1中，背景图像部分的获取过程如下：
13.步骤1.1，设原始高光谱图像中第i个像元的光谱向量为xi∈rb；
14.步骤1.2，设定原始高光谱图像中的背景像元服从高斯分布，通过构造目标与背景窗口，在图像局部区域利用马氏距离计算待测像元与背景像元之间的差异值，并根据差异值大小来判断待测像元类型，根据下式计算得出判断结果：
[0015][0016]
式中，μb表示背景的均值，cb表示背景的协方差矩阵，η表示rx算法判定阈值，如果rx(xi)《η则初步判断待测像元为背景部分，否则，初步判断待测像元为异常部分；
[0017]
步骤1.3，针对包含超过图像总像素数百分之一的大异常目标图像，将rx算法阈值设定为500，采用固定值替换法得到背景图像；针对包含少于图像总像素数百分之一的小异常目标图像，将rx算法阈值设定为350，采用周围像元替换法得到背景图像；最终将原始高光谱图像根据rx算法的检测结果分为背景图像部分与异常图像部分。
[0018]
所述步骤2中，栈式降噪自编码器的训练过程如下：
[0019]
步骤2.1，采用栈式降噪自编码器网络，在训练时使用基于光谱角度余弦cos(θ
sa
)的重建损失函数，其中cos(θ
sa
)的计算式如下：
[0020][0021]
式中a和b表示自编码器网络的输出重建光谱向量与输入光谱向量，o表示光谱的波段数；
[0022]
将光谱角度纳入到重建损失函数e
sa
中有如下公式：
[0023][0024]
式中，f(z
(l)
)代表了自编码器的输出重建光谱向量a，y代表了输入光谱向量b，l表示自编码器的层数，f(z
(l)
)中的z
(l)
表示为：
[0025]z(l)
＝w
(l-1)a(l-1)
b
(l-1)
[0026]
式中，a
(l)
＝f(z
(l)
)表示经过自编码器l层处理过后的光谱向量，a
(1)
＝x为原始光谱向量，l＝l,l-1,l-2,
…
,2，l表示自编码器层数，w
(l-1)
、b
(l-1)
分别表示自编码器第l-1层的权重矩阵与偏置矩阵；
[0027]
将正则项包括在内，得到最终的重建损失函数e
sa
如下式，
[0028][0029]
式中，w、b分别表示权重矩阵与偏置矩阵，y
(m)
代表了第m个像元的输入光谱向量，λ为正则化参数，i、j分别代表第l层和l 1层中的像元数量；
[0030]
步骤2.2，栈式降噪自编码器网络中使用sigmoid函数作为激活函数，在其隐藏层中，设定隐藏层1到5层基于隐藏层3神经元数量呈对称分布，第1个到第3个自编码器中的神经元数量逐层下降；
[0031]
步骤2.3，在输入背景图像后，首先进行预训练，第一步训练从输入层到隐藏层1和隐藏层5到输出层的参数，实现对输入层添加高斯噪声的图像的重建，训练完成后，固定训练所得到的参数；第二步训练从隐藏层1到隐藏层2和隐藏层4到隐藏层5的参数，实现对隐藏层1的输出图像的重建，训练完成后，固定训练所得到的参数；第三步训练最中间隐藏层3的输入和输出参数，实现对隐藏层2的输出图像的重建，训练完成后，固定训练所得到的参数；
[0032]
步骤2.4，经过预训练，网络中的各个参数值均初始化到了合适的数值，再将训练得到的参数的隐藏层整合到一起进行训练，实现对输入层添加高斯噪声的图像的重建，对各层之前训练得到参数进行小范围的调整，最终得到训练好的栈式降噪自编码器。
[0033]
所述步骤2.2，网络学习率设定为10-3
，epoch设定为100，批大小设定为400。
[0034]
所述步骤4，重建误差ri的获得步骤如下：
[0035]
步骤4.1，对于重建高光谱图像x
′
∈r
mn
×b，按照下式计算每个像元的重建光谱向量xi′
与原始光谱向量xi的差值di，从而获得光谱向量差值图像；
[0036]di
＝x
i-xi′
[0037]
步骤4.2，对获得的光谱向量差值图像使用局部标准偏差过滤器预处理，经过过滤器滤波后，异常像元的局部标准偏差高于背景像元的局部标准偏差，从而通过滤波获得异常更为突出的重建误差ri。
[0038]
所述步骤5，异常评分ad(ri)的计算过程如下：
[0039]
步骤5.1，使用多元偏态t分布(multivariate skewed t-distribution，mvskt)模型模拟重建误差ri，在该模型中，ri表示为下式：
[0040][0041]
式中，wi表示高斯随机向量，z表示广义逆高斯分布的随机变量，m为均值向量，b为偏置向量，t为逆协方差矩阵，设定b》1,α《0,γ
→
∞，且m和t满足正态-wishart先验分布，如下式所示：
[0042][0043]
式中，m0、λ、v0和ψ0为超参数；
[0044]
则多元偏态t分布模型的概率分布函数如下式所示：
[0045][0046]
式中ri＝(r
i-m)
t
t(r
i-m)，k
α
(
·
)为第二类α阶修正贝塞尔(bessel)函数,为第二类阶修正贝塞尔(bessel)函数；
[0047]
步骤5.2，当参数确定后，计算概率分布函数的负对数作为异常评分从而实现异常检测，则其异常评分如下式所示，
[0048][0049]
步骤5.3，将步骤4中所得到的重建误差ri输入到异常评分公式中计算，即得到最终异常评分ad(ri)；
[0050]
步骤5.4，利用阈值分割检测方法根据异常评分ad(ri)判断待测像元类别，得到最终异常检测结果。
[0051]
与现有技术相比，本发明对参数的要求更低，算法的适应性更好，并在此基础上实现了良好的异常检测，具有良好的鲁棒性。
附图说明
[0052]
图1为本发明原始高光谱图像的伪彩色图。
[0053]
图2为本发明流程图。
[0054]
图3为rx算法预处理结果背景图像部分图。
[0055]
图4为栈式降噪自编码器网络结构图。
[0056]
图5为栈式降噪自编码器训练过程auc结果图。
[0057]
图6为本发明实例中异常检测结果对比图，其中(a)为mvc方法检测结果，(b)为mvl方法检测结果，(c)为mvskt方法检测结果。
[0058]
图7为本发明实例中异常检测结果对比三维曲线图，其中(a)为mvc方法检测结果，(b)为mvl方法检测结果，(c)为mvskt方法检测结果。
[0059]
图8为本发明实例中异常检测结果对比roc图。
[0060]
表1为本发明实例中异常检测结果对比auc图。
具体实施方式
[0061]
为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。
[0062]
如前所述，现有的高光谱图像异常检测方法中，由于评分机制不完善和重建误差不精准等方面原因，使得检测的效果欠佳，尤其是缺乏鲁棒性。因此，建立一个新型的评分机制并保证高准确性和高鲁棒性是高光谱异常检测领域的一个重点。本发明即是对上述问题作出的一系列改进。
[0063]
为论证本发明方法的有效性，使用真实高光谱图像来进行异常检测。实验使用abu-airport图像，如图1所示，数据来源于机载可见光/红外成像光谱仪(aviris)采集的airport-beach-urban(abu)高光谱数据集，图像的空间分辨率为一个像元17.2m，光谱波段从365到2500nm。实验所用的图像大小为100
×
100，波段数为205个，并将图像中的两架飞机作为异常目标。参考图2，针对abu-airport图像数据的异常检测基本流程如下：
[0064]
步骤1，对原始高光谱图像x∈r
mn
×b使用rx算法进行预处理，得到属于背景部分的高光谱图像，即其背景图像部分。其中m、n表示图像的空间维度，mn即像元数，b表示光谱维
度。具体步骤可表述如下：
[0065]
步骤1.1，原始高光谱图像是一个100
×
100
×
205的三维矩阵，共有100
×
100个空间像素点，包含地面目标在205个波段内的光谱数据，与每一个空间像素相对应1
×1×
205大小的一维向量表示对应的光谱信息。设原始高光谱图像中第i个像元的光谱向量为xi∈rb。
[0066]
步骤1.2，假定原始高光谱图像中的背景像元服从高斯分布，通过构造目标与背景窗口，在图像局部区域利用马氏距离计算待测像元与背景像元之间的差异值，并根据差异值大小来判断待测像元类型，根据下式计算得出判断结果。
[0067][0068]
式中，μb表示背景的均值，cb表示背景的协方差矩阵，η表示rx算法判定阈值，此公式的含义为：如果rx(xi)《η则初步判断待测像元为背景部分，否则，初步判断待测像元为异常部分。
[0069]
步骤1.3，对飞机图像采用rx算法进行异常检测。
[0070]
rx算法的异常检测情况对于后续自编码器网络的训练至关重要，而rx算法是根据设定阈值判定异常，所以需要对阈值的选择进行实验研究。为了便于对比观察，将检测到的异常像元值替换为固定值0，在图像中表示为纯黑色像元。本发明中，将包含超过图像总像素数百分之一的异常目标图像称为大异常目标图像，将包含少于图像总像素数百分之一的异常目标图像称为小异常目标图像。针对大异常目标图像，将rx算法阈值设定为500，采用固定值替换法得到背景图像；针对小异常目标图像，将rx算法阈值设定为350，采用周围像元替换法得到背景图像。示例地，在本实施例中，飞机目标较大，显然属于大目标，因此将rx算法阈值设定为500，并采用固定值替换法得到背景图像。最终将原始高光谱图像根据rx算法的检测结果分为背景图像部分与异常图像部分，共两部分，处理结果背景图像部分如图3所示。
[0071]
步骤2，将背景图像部分作为基于光谱角度余弦的栈式降噪自编码器网络的输入，对网络进行训练以获得训练好的栈式降噪自编码器。具体的训练过程如下：
[0072]
步骤2.1，本发明采用栈式降噪自编码器网络，其将包含五个隐藏层的栈式自编码器与添加随机高斯噪声的降噪自编码器相结合，网络的具体结构如图4所示。
[0073]
在训练时，使用基于光谱角度余弦cos(θ
sa
)的重建损失函数，其中cos(θ
sa
)的计算式如下：
[0074][0075]
式中a和b表示自编码器网络的输出重建光谱向量与输入光谱向量，o表示光谱的波段数；
[0076]
将光谱角度纳入到重建损失函数e
sa
中有如下公式：
[0077][0078]
式中，f(z
(l)
)代表了自编码器的输出重建光谱向量a，y代表了输入光谱向量b，l表示自编码器的层数，f(z
(l)
)中的z
(l)
表示为：
[0079]z(l)
＝w
(l-1)a(l-1)
b
(l-1)
[0080]
式中，a
(l)
＝f(z
(l)
)表示经过自编码器l层处理过后的光谱向量，a
(1)
＝x为原始光谱向量，l＝l,l-1,l-2,
…
,2，l表示自编码器层数，w
(l-1)
、b
(l-1)
分别表示自编码器第l-1层的权重矩阵与偏置矩阵；
[0081]
将正则项包括在内，得到最终的重建损失函数e
sa
如下式，
[0082][0083]
式中，w、b分别表示权重矩阵与偏置矩阵，y
(m)
代表了第m个像元的输入光谱向量，λ为正则化参数，i、j分别代表第l层和l 1层中的像元数量。
[0084]
步骤2.2，栈式降噪自编码器网络中使用sigmoid函数作为激活函数，网络学习率设定为10-3
，epoch设定为100，批大小设定为400。在其隐藏层中，设定隐藏层1到5层基于隐藏层3神经元数量呈对称分布，考虑到网络数据特征提取的目的，第1个到第3个自编码器中的神经元数量逐层下降，而为了避免变量过多，本实施例中设定第1个和第2个编码器中的编码单元数量分别为100和50，即隐藏层1和隐藏层5中神经元数量为100，隐藏层2和隐藏层4中神经元数量为50，而最中间隐藏层3的神经元数量，由于其直接关系到最终输出图像的光谱维数，分别设定5、10、15、20、25、30、35共七种数值作为隐藏层3的神经元数量进行实验，分别用三种自编码器网络mse-sdae、csa-sdae、sid-sdae对pavia桥梁图像和urban道路图像进行训练，并结合协同表示算法进行异常检测，根据auc结果进行评判，选择30作为隐藏层3的神经元数量，auc结果如图5所示。
[0085]
步骤2.3，在输入背景图像后，首先进行预训练，第一步训练从输入层到隐藏层1和隐藏层5到输出层的参数，实现对输入层添加高斯噪声的图像的重建，训练完成后，固定训练所得到的参数；第二步训练从隐藏层1到隐藏层2和隐藏层4到隐藏层5的参数，实现对隐藏层1的输出图像的重建，训练完成后，固定训练所得到的参数；第三步训练最中间隐藏层3的输入和输出参数，实现对隐藏层2的输出图像的重建，训练完成后，固定训练所得到的参数。
[0086]
步骤2.4，经过前面三步的预训练，网络中的各个参数值均初始化到了合适的数值。再将前三步训练得到的参数的隐藏层整合到一起进行训练，实现对输入层添加高斯噪声的图像的重建，对各层之前训练得到参数进行小范围的调整，最终得到训练好的栈式降噪自编码器。
[0087]
步骤3，对训练好的栈式降噪自编码器，输入原始高光谱图像进行图像重建。具体步骤包括：将原始高光谱图像数据x∈r
mn
×b作为输入，传入到基于光谱角度余弦的栈式降噪自编码器中，选择输出层结果输出，得到重建的高光谱图像数据x
′
∈r
mn
×b。
[0088]
步骤4，计算重建的高光谱图像中第i个像元的重建光谱向量xi′
与原始光谱向量xi的差值di，从而得到光谱向量差值图像，再对光谱向量差值图像使用局部标准偏差过滤器，从而通过滤波获得异常更为突出的重建误差ri∈rb(i＝1,2,
…
,mn)。具体步骤可描述如下：
[0089]
步骤4.1，对于自编码器输出的重建高光谱图像x
′
∈r
mn
×b，按照下式计算每个像元的重建光谱向量xi′
与原始光谱向量xi的差值di，从而获得光谱向量差值图像。
[0090]di
＝x
i-xi′
[0091]
步骤4.2，对获得的光谱向量差值图像使用局部标准偏差过滤器预处理，经过过滤器滤波后，异常像元的局部标准偏差要高于背景像元的局部标准偏差，从而通过滤波获得异常更为突出的重建误差ri∈rb(i＝1,2,
…
,mn)。
[0092]
步骤5，将重建误差ri输入到mvskt分布模型中，在不同的阈值下计算得到异常评分ad(ri)，根据异常评分ad(ri)判断待测像元类别，输出异常检测结果。并根据检测结果中的虚警概率和检测概率绘制相应roc曲线、auc曲面图，通过与其他检测模型比较分析mvskt分布模型的性能。具体步骤可描述如下：
[0093]
步骤5.1，本发明使用多元偏态t分布(multivariate skewed t-distribution，mvskt)模型来模拟重建误差ri，在该模型中ri可表示为下式，
[0094][0095]
式中，wi表示高斯随机向量，z表示广义逆高斯分布的随机变量，m为均值向量，b为偏置向量，t为逆协方差矩阵，设定b》1,α《0,γ
→
∞，且m和t满足正态-wishart先验分布，如下式所示，
[0096][0097]
式中，m0、λ、ν0和ψ0为超参数。
[0098]
则多元偏态t分布模型的概率分布函数如下式所示，
[0099][0100]
式中ri＝(r
i-m)
t
t(r
i-m)，k
α
(
·
)为第二类α阶修正贝塞尔(bessel)函数,为第二类阶修正贝塞尔(bessel)函数。
[0101]
步骤5.2，计算模型中的各个参数，采用vb方法，首先需要考虑最大边际似然估计问题，如下式所示。
[0102][0103]
其中，θ表示mvskt分布模型的参数集合。对潜在变量φ＝{m,t,z}进行积分，使用近似后验分布q(φ)，式中的边际对数似然可表示为
[0104]
logf(r
1:n
∣θ)＝f(q,θ) d
kl
(f(φ∣r
1:n
,θ)‖q)
[0105]
式中f(q,θ)和d
kl
(f||q)的求解分别如下式所示。
[0106]
f(q,θ)＝∫log((f(r
1:n
,φ∣θ))/q(φ))q(φ)dφ
[0107]dkl
9f‖q)＝-∫log((f(φr
1:n
,θ))/q(φ))q(φ)dφ
[0108]
其中，d
kl
(f||q)表示f(φ|r
1:n
,θ)与q(φ)之间的kl散度。
[0109]
vb方法使用精确后验值的因式分解来近似表示后验分布q(φ))，如下式，
[0110]
f(φ∣r
1:n
,θ)≈q(φ)＝q(m)q(t)q(z)
[0111]
这种近似方法方便计算期望积分，因此在vb方法中，需要首先找到q(m)、q(t)和q(z)，然后以迭代的方式更新参数集。
[0112]
对于q(m)，其分布应满足下式：
[0113][0114]
其中const.表示常量。
[0115]
求解得到的q(m)应服从多元正态分布，即其中
[0116][0117]cm
＝((e[z-1
]n λ)e[t])-1
[0118][0119]
对于q(t)，其分布应满足下式，
[0120][0121]
求解得到的q(t)应服从wishart分布，即其中
[0122][0123]
vn＝(τ n 1)
[0124]
对于q(z)，其分布应满足下式
[0125][0126]
求解得到的q(z)应服从广义逆高斯分布，即其中
[0127]
c＝ntr(e(t)bb
t
)
[0128][0129]
其中的期望值计算式如下，
[0130]
e[m]＝μm[0131][0132][0133][0134]
初始参数设置完毕后，接下来就是迭代更新参数，通过将边际对数似然表示最大化，即最大化f(q,θ)的下限来寻找参数，其中m0、ψ0、β、λ的下限分别由下式给出。
[0135][0136][0137][0138][0139]
计算最大化的估计量，其中的变量及参数为
[0140]
m0＝e[m]
[0141]
ψ0＝(τ/τ-l-1)(e[t])-1
[0142]
β＝-(2α/(e[z-1
]))
[0143][0144]
其中一些参数如下所示，
[0145]
α＝l/2
[0146]
τ＝(((n-l-1)p l 1)/1-p)
[0147]
p＝((τ-l-1)/(n τ-l-1))
[0148]
0≤p≤1
[0149]
步骤5.3，当参数确定后，计算概率分布函数的负对数作为异常评分从而实现异常检测，则其异常评分ad如下式所示，
[0150][0151]
步骤5.4，将步骤4中所得到的重建误差ri∈rb(i＝1,2,
…
,mn)输入到异常评分公式中计算，即得到最终异常评分ad(ri)。利用阈值分割检测方法根据异常评分ad(ri)判断待测像元类别，得到最终异常检测结果。
[0152]
步骤5.5，绘制相应最终异常评分ad(ri)结果伪彩图，如图6所示，绘制相应最终异常评分ad(ri)结果三维曲线图，如图7所示，相较于其他模板对异常有更准确的凸显。
[0153]
步骤5.6，依据模板进行虚警率和准确率估计绘制成roc曲线图和auc面积图。roc曲线图如图8所示，auc统计结果在考虑运行时间的基础上与其他模型进行比较，结果如下表所示。
[0154]
算法名称auc时间(s)grx0.84030.0693lrx0.905241.1269cr0.90696.6130mvc0.91156.9031mvj0.91054.7637mvl0.91014.7676mvst0.91054.8084mvn0.84580.1245本章mvskt0.93536.7372
[0155]
根据比对结果可知本发明提出的mvskt算法在abu-airport图像上的异常检测auc值明显高于其它算法，grx与mvn算法的auc值最低，除此之外，其它几种基于自编码器的异常检测算法效果都要好于传统算法，可见对于复杂场景的异常检测，基于自编码器的正态混合均值-方差分布方法的检测效果要比传统方法更加适用些，而相比其它几种多元正态混合均值-方差分布，mvskt算法则是进一步提升了自编码器方法的检测性能。在算法的运行时间上，本发明所提出算法的耗时位于中间水平，与几种多元分布模型算法的耗时较为接近。与所使用的传统算法相比，基于自编码器的多元分布算法对于参数的要求更低，算法的适应性相对更好。