CN107895361A - 一种基于局部密度纯化背景的高光谱异常检测方法 - Google Patents

Abstract

本发明公开了一种基于局部密度纯化背景的高光谱异常检测方法，首先，利用同心双层窗模型获取当前被检测像元对应的初始背景。然后，计算初始背景各像元的局部密度。接着，设定最大异常比例，按该比例选取最小局部密度对应的像元。再利用最大类间方差法对其进一步分割。最后，结合LRXD异常检测方法对高光谱图像进行检测。设定阈值，将检测值大于阈值的像元标记为异常点。本发明通过纯化背景去除背景中的异常数据，有助于分析目标与背景中的差异，从而有效的降低虚警率。

Description

一种基于局部密度纯化背景的高光谱异常检测方法

技术领域

本发明涉及高光谱异常检测方法，尤其涉及一种基于局部密度纯化背景的高光谱异常检测方法。

背景技术

高光谱图像是一种图谱合一的三维新型遥感数据，其光谱分辨率高，具有许多连续的光谱波段。这为目标检测提供了丰富的判别信息。根据目标检测是否利用先验信息(一般为目标光谱)，将其分为两类：有监督目标检测与无监督目标检测。无监督目标检测也称为异常目标检测，异常检测在未知先验信息的情况下，直接利用背景与目标的光谱差异，检测图像异常点。异常检测方法为背景和目标分布未知的检测识别问题提供了解决方式。

经典的异常检测方法是Reed等人利用广义似然比检验推导出的RX检测方法。RX检测方法根据选取背景方式不同，分为局部RX检测方法(LRXD)和全局RX检测方法(GRXD)。GRXD选取高光谱图像的所有像元作为背景，而LRXD基于同心双层窗模型选择检测像元的局部近邻对背景统计信息进行估计。RX检测方法假设背景服从多维高斯分布模型，且异常目标在高光谱图像中是低概率存在的。然而，在实际情况中，高光谱图像一般存在多种地物，同时受阴影、光照、大气干扰的影响，获取的高光谱图像中背景复杂多变。因此多维高斯分布模型不能全面体现真实高光谱图像地物的分布特性。只要存在0.5％的数据污染就会产生协方差失真的问题。对于高光谱图像，背景的复杂性造成了统计信息的失真，这使得RX检测方法的虚警率较高。针对这一问题，通过纯化背景将潜在异常目标去除，使背景更为符合多维高斯分布模型。再利用背景纯化方法得到更为准确的背景统计信息(背景协方差与均值向量)，降低RX检测方法的虚警率。

发明内容

本发明目的是要提供一种基于局部密度纯化背景的高光谱异常检测方法，它能有效地降低RX检测方法的虚警率。

技术方案：本发明方法先获取纯化后背景，然后增强背景与目标差异程度，从而降低检测方法虚警率，包括如下步骤：

步骤1，输入高光谱图像；

步骤2，遍历图像中的每个像元，当前像元为被检测像元，通过同心双层窗获取对应的初始背景；

步骤3，计算初始背景中的每个像元的局部密度；

步骤4，进行双重分割，得到纯化后背景；

步骤5，根据纯化后背景，利用异常检测方法对高光谱图像的被检测像元进行检测；

步骤6，高光谱图像检测完毕，根据检测值得到异常点，输出结果。

步骤2中，同心双层窗模型为具有相同中心的两个矩形窗的叠加，分内窗与外窗。双层窗尺寸根据感兴趣目标大小和检测全图大小确定，一般内窗尺寸略大于目标尺寸，外窗尺寸要保证背景集合具有一定数量的像元。采用同心双层窗模型获取对应的初始背景，被检测像元为同心双层窗的中心像元，同心双层窗窗口尺寸表示为(w_in,w_out)，w_in为内窗尺寸，w_out为外窗尺寸。可以计算出初始背景像元数：N＝w_out×w_out-w_in×w_in。这里w_in设置为5，w_out设置为11。初始背景矩阵表示为X_b＝[x₁,…,x_j,…,x_N]，像元x_j＝[x_j(1),x_j(2),…,x_j(L)]^T，其中，第j个像元表示为x_j＝[x_j(1),x_j(2),…,x_j(L)]^T，x_j(L)为第L波段的光谱辐射度或反射率，L为高光谱图像波段数，由成像仪的特性决定，被检测像元表示为y＝[y(1),y(2),…,y(L)]^T，y(L)表示第L波段的被检测像元。

步骤3中，对于像元p，其局部密度DEN(p)计算公式如下所示：

DEN(p)＝|{q|distance(q,p)<d,d>0}|

式中，p和q为初始背景X_b中的像元，distance(·)表示像元距离，这里使用欧氏距离进行度量，在多维空间中，以一点为中心，一定距离内的空间为超球体，这里p为球心，q为超球体中的像元，d表示像元超球体的区域半径，|·|表示计算超球体中像元的数目。这里设定d的取值范围为[0.01,0.05]。根据上式，计算初始背景X_b各个像元x_j的局部密度den_j，得到初始背景局部密度向量den＝[den₁,…,den_N]。局部密度越小，像元异常程度越高。

步骤4包括：

步骤4-1，进行第一重分割，设定最大异常比例，按该比例选取最小局部密度值及对应的像元：将最大异常比例设置为20％，分割过程如下所示：

对应的：

其中，den_l为den中80％最大元素，X_l为X_b中对应的80％局部密度最大的像元；den_s为den中80％最大元素，X_s为X_b中对应的80％局部密度最大的像元；

步骤4-2，进行第二重分割，利用最大类间方差法，对den_s进行处理：灰度化为den_sg。den_sg＝round(255×den_s/max(den_s))，

式中，max(·)表示其中元素的最大值，round(·)表示根据四舍五入取整。

然后遍历灰度值，取值范围为[0,255]，以灰度值为界将den_sg分成异常和背景两部分，计算两部分间的方差，当两部分间的方差最大时记录该灰度值th：

式中，w₀表示den_sg中不大于th的元素所占比例，w₁表示den_sg中大于th的元素所占比例；μ₀、μ₁、μ分别为den_sg中不大于th的元素的均值、den_sg中大于th的元素的均值、den_sg的均值。当存在多个灰度值满足条件时，取它们的平均值作为th。接着，将th去灰度化，得到分割阈值th_f。

th_f＝max(den_s)×th/255

最后，纯化后背景X_rb表示为：

X_rb＝{x_j|den_j>th_f,den_j＝DEN(x_j)},j＝1,2,…,N

步骤5中，采用的异常检测方法为局部RX算法(local Reed-Xiaoli Detector,LRXD)，被检测像元y＝[y(1),y(2),...,y(L)]^T的检测值计算公式如下所示：

LRXD(y)＝(y-μ)^TK^-1(y-μ)，

式中，μ为纯化后背景X_rb的均值向量，K为X_rb的协方差矩阵。

其中，σ_ij＝(1/N)∑_k(x_ik-m_i)(x_jk-m_j),k＝1,...,N，x_ik表示第k个像元第i个波段取值，m_i，m_j分别为第i个波段的均值和第j个波段的均值。

步骤6中，通过设定置信度的方法设定判别阈值，将检测值大于判别阈值的像元标记为异常点，输出结果。

步骤6中，还可以设定一定量数目，数目由实验者根据经验选择，将检测值较大的像元判为异常点，输出结果。

有益效果：

1、本发明的纯化背景方法作为一种预处理手段，可以在异常检测方法执行之前将初始背景中的潜在异常像元去除。这样高光谱异常检测方法可免受潜在异常干扰从而增强异常目标与背景间的差异程度，降低异常检测方法的虚警率。

2、本发明的纯化背景方法基于局部密度概念与双重分割实现。像元局部密度表征像元近邻域内像元数目。像元数目越多，局部密度越大。也就是说局部密度越大的像元，其异常程度越低，属于背景的概率越大。第一重分割，设定最大异常比例，按该比例选取最小局部密度对应的像元。加入第一重分割是为了防止最大类间方差法将背景过分割。而第二重分割利用最大类间方差法进一步分割，得到分割阈值。将局部密度大于分割阈值的像元提取为纯化背景。

附图说明

下面结合附图和具体实施方式对本发明做更进一步的具体说明，本发明的上述或其他方面的优点将会变得更加清楚。

图1为一幅AVIRIS高光谱图像；

图2a为AVIRIS高光谱图像第60波段图像；

图2b用于仿真实验截取子图像

图2c地面目标分布图；

图3为同心双层窗模型；

图4为本文方法(Density Background Refinement Based Anomaly Detector,DBRAD)在d取不同值时的ROC曲线比较

图5为本文方法同其它方法的接收机工作特性(receiver operatorcharacteristic,ROC)曲线比较，对比方法包括LRXD、基于随机选取的方法(Random-Selection-Based Anomaly Detector,RSAD)、基于概率统计的异常检测方法(Probabilistic Anomaly Detector,PAD)。

图6a为LRXD检测方法检测结果图。

图6b为RSAD检测方法检测结果图。

图6c为PAD检测方法检测结果图。

图6d为DBRAD检测方法检测结果图。

具体实施方式

下面结合附图及实施例对本发明做进一步说明。

本发明包括如下步骤：

步骤(1)输入高光谱图像，如图2a和图2b所示，图2b为图1的子图像，空间大小100×100。图2c给出了目标分布位置，用于判断检测出的像元是否为真实目标所在位置。

步骤(2)对当前高光谱图像的被检测像元，利用同心双层窗模型(如图3所示)获取对应的初始背景。同心双层窗模型为具有相同中心的两个矩形窗的叠加，分内窗与外窗。双层窗尺寸根据感兴趣目标大小和检测全图大小确定，一般内窗尺寸略大于目标尺寸，外窗尺寸要保证背景集合具有一定数量的像元。窗口尺寸表示为(w_in,w_out)，w_in为内窗尺寸，w_out为外窗尺寸。可以计算出初始背景像元数：N＝w_out×w_out-w_in×w_in。这里w_in设置为5，w_out设置为11。初始背景矩阵表示为X_b＝[x₁,…,x_j,…,x_N]，像元x_j＝[x_j(1),x_j(2),…,x_j(L)]^T，L为高光谱图像波段数。

步骤(3)计算初始背景中的每个像元的局部密度。对于像元p，其局部密度计算公式如下所示：

DEN(p)＝|{q|distance(q,p)<d,d>0}|

式中，p和q为初始背景X_b中的像元，distance(·)表示像元距离，这里使用欧氏距离进行度量，d表示像元超球体的区域半径，|·|表示计算超球体中像元的数目。这里设定d的取值范围为[0.01,0.05]。根据上式，计算初始背景X_b各个像元x_j的局部密度den_j，得到初始背景局部密度向量den＝[den₁,…,den_N]。局部密度越小，像元异常程度越高。d取值不同，最终检测结果也会受到影响，如图4所示。图4给出了d取不同值(0.05、0.03或0.01)时DBRAD的ROC曲线比较情况。可以看出当d取值为0.01时，检测效果最好，在相同虚警率下具有较高的检测率。ROC曲线、虚警率和检测率为方法性能评价的通用指标。

步骤(4)双重分割：

(4-1)第一重分割，设定最大异常比例，按该比例选取最小局部密度值及对应的像元；

(4-2)第二重分割，利用最大类间方差法对(4-1)中选取部分进行处理，得到一个分割阈值。将局部密度大于该阈值的像元提取为纯化后背景。

将最大异常比例设置为20％，分割过程如下所示：

对应的：

利用最大类间方差法，对den_s进行处理。首先，灰度化为den_sg。

den_sg＝round(255×den_s/max(den_s))

式中，max(·)表示其中元素的最大值，round(·)表示根据四舍五入取整。

然后遍历灰度值，取值范围为[0,255]，以灰度值为界将den_sg分成异常和背景两部分，计算两部分间的方差，当两部分间的方差最大时记录该灰度值th：

式中，w₀表示den_sg中不大于th的元素所占比例，w₁表示den_sg中大于th的元素所占比例；μ₀、μ₁、μ分别为den_sg中不大于th的元素的均值、den_sg中大于th的元素的均值、den_sg的均值。当存在多个灰度值满足条件时，取它们的平均值作为th。接着，将th去灰度化，得到分割阈值th_f。

th_f＝max(den_s)×th/255

最后，纯化后背景X_rb表示为：

X_rb＝{x_j|den_j>th_f,den_j＝DEN(x_j)},j＝1,2,…,N

步骤(5)根据纯化后背景，利用LRXD异常检测方法对高光谱图像的被检测像元进行检测。采用的异常检测方法为LRXD，被检测像元y＝[y(1),y(2),…,y(L)]^T的检测值计算公式如下所示：

LRXD(y)＝(y-μ)^TK^-1(y-μ)

式中，μ为纯化后背景X_rb的均值向量，K为X_rb的协方差矩阵。

步骤(6)高光谱图像检测完毕，设定判别阈值，将检测值大于判别阈值的像元标记为异常点，输出结果。

将本发明提供的检测方法与当前主流异常检测方法效果进行比较，如图5、图6a～图6d所示。图5给出了DBRAD(d＝0.01)同LRXD，RSAD，PAD算法的ROC曲线比较。可以看出DBRAD的ROC曲线变化更快，检测率高于其它算法。图6a～图6d给出了检测结果二值图。可以看出，相对PAD、RSAD、LRXD有些目标检测结果不够清晰完整，DBRAD可以有效的检测出目标，虚警点更少，细节体现更为全面。

本发明提供了一种基于局部密度纯化背景的高光谱异常检测方法，具体实现该技术方案的方法和途径很多，以上所述仅是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以做出若干改进和润饰，这些改进和润饰也应视为本发明的保护范围。本实施例中未明确的各组成部分均可用现有技术加以实现。

Claims (7)

1.一种基于局部密度纯化背景的高光谱异常检测方法，其特征在于，包括如下步骤：

步骤1，输入高光谱图像；

步骤2，遍历图像中的每个像元，当前像元为被检测像元，通过同心双层窗获取对应的初始背景；

步骤3，计算初始背景中的每个像元的局部密度；

步骤4，进行双重分割，得到纯化后背景；

步骤5，根据纯化后背景，利用异常检测方法对高光谱图像的被检测像元进行检测；

步骤6，高光谱图像检测完毕，根据检测值得到异常点，输出结果。

2.根据权利要求1所述的方法，其特征在于，步骤2中，采用同心双层窗模型获取对应的初始背景，被检测像元为同心双层窗的中心像元，同心双层窗窗口尺寸表示为(w_in,w_out)，w_in为内窗尺寸，w_out为外窗尺寸，计算初始背景像元数N：

N＝w_out×w_out-w_in×w_in，

初始背景X_b表示为：

X_b＝[x₁,…,x_j,…,x_N]，

其中，第j个像元表示为x_j＝[x_j(1),x_j(2),…,x_j(L)]^T，x_j(L)为第L波段的光谱辐射度或反射率，L为高光谱图像波段数，由成像仪的特性决定，被检测像元表示为y＝[y(1),y(2),…,y(L)]^T，y(L)表示第L波段的被检测像元。

3.根据权利要求2所述的方法，其特征在于，步骤3中，对于像元p，其局部密度DEN(p)计算公式如下所示：

DEN(p)＝|{q|distance(q,p)<d,d>0}|，

其中，p和q为初始背景X_b中的像元，distance(·)表示像元距离，这里使用欧氏距离进行度量，d表示像元超球体的区域半径，|·|表示计算超球体中像元的数目，根据上式，计算初始背景X_b各个像元x_j的局部密度den_j，得到初始背景局部密度向量den＝[den₁,…,den_N]。

4.根据权利要求3所述的方法，其特征在于，步骤4包括：

步骤4-1，进行第一重分割，设定最大异常比例，按该比例选取最小局部密度值及对应的像元：

将最大异常比例设置为20％，分割过程如下所示：

对应的：

其中，den_l为den中80％最大元素，X_l为X_b中对应的80％局部密度最大的像元；den_s为den中80％最大元素，X_s为X_b中对应的80％局部密度最大的像元；

步骤4-2，进行第二重分割，利用最大类间方差法，对den_s进行处理：

将den_s灰度化为den_sg：

den_sg＝round(255×den_s/max(den_s))，

式中，max(·)表示其中元素的最大值，round(·)表示根据四舍五入取整；

遍历灰度值，取值范围为[0,255]，以灰度值为界将den_sg分成异常和背景两部分，计算两部分间的方差，当两部分间的方差最大时记录该灰度值th：

<mrow> <mi>t</mi> <mi>h</mi> <mo>=</mo> <mi>arg</mi> <mi> </mi> <munder> <mrow> <msub> <mi>maxw</mi> <mn>0</mn> </msub> </mrow> <mrow> <mi>t</mi> <mi>h</mi> </mrow> </munder> <msup> <mrow> <mo>(</mo> <msub> <mi>&amp;mu;</mi> <mn>0</mn> </msub> <mo>-</mo> <mi>&amp;mu;</mi> <mo>)</mo> </mrow> <mn>2</mn> </msup> <mo>+</mo> <msub> <mi>w</mi> <mn>1</mn> </msub> <msup> <mrow> <mo>(</mo> <msub> <mi>&amp;mu;</mi> <mn>1</mn> </msub> <mo>-</mo> <mi>&amp;mu;</mi> <mo>)</mo> </mrow> <mn>2</mn> </msup> <mo>,</mo> </mrow>

其中，w₀表示den_sg中不大于th的元素所占比例，w₁表示den_sg中大于th的元素所占比例；μ₀、μ₁、μ分别为den_sg中不大于th的元素的均值、den_sg中大于th的元素的均值、den_sg的均值，当存在两个以上灰度值满足条件时，取它们的平均值作为th；

将th去灰度化，得到分割阈值th_f：

th_f＝max(den_s)×th/255

纯化后背景X_rb表示为：

X_rb＝{x_j|den_j>th_f,den_j＝DEN(x_j)},j＝1,2,…,N。

5.根据权利要求4所述的方法，其特征在于，步骤5中，采用的异常检测方法为局部RX算法，被检测像元y＝[y(1),y(2),…,y(L)]^T的检测值LRXD(y)计算公式如下所示：

LRXD(y)＝(y-μ)^TK^-1(y-μ)，

其中，μ为纯化后背景X_rb的均值向量，K为X_rb的协方差矩阵，

<mrow> <mi>&amp;mu;</mi> <mo>=</mo> <mfrac> <mn>1</mn> <mi>N</mi> </mfrac> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>j</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>N</mi> </munderover> <msub> <mi>x</mi> <mi>j</mi> </msub> <mo>,</mo> <msub> <mi>x</mi> <mi>j</mi> </msub> <mo>&amp;Element;</mo> <msub> <mi>X</mi> <mrow> <mi>r</mi> <mi>b</mi> </mrow> </msub> </mrow>

<mrow> <mi>K</mi> <mo>=</mo> <mfenced open = "[" close = "]"> <mtable> <mtr> <mtd> <msub> <mi>&amp;sigma;</mi> <mn>11</mn> </msub> </mtd> <mtd> <msub> <mi>&amp;sigma;</mi> <mn>12</mn> </msub> </mtd> <mtd> <mo>...</mo> </mtd> <mtd> <msub> <mi>&amp;sigma;</mi> <mrow> <mn>1</mn> <mi>L</mi> </mrow> </msub> </mtd> </mtr> <mtr> <mtd> <mo>...</mo> </mtd> <mtd> <mrow></mrow> </mtd> <mtd> <mrow></mrow> </mtd> <mtd> <mrow></mrow> </mtd> </mtr> <mtr> <mtd> <msub> <mi>&amp;sigma;</mi> <mrow> <mi>L</mi> <mn>1</mn> </mrow> </msub> </mtd> <mtd> <msub> <mi>&amp;sigma;</mi> <mrow> <mi>L</mi> <mn>2</mn> </mrow> </msub> </mtd> <mtd> <mo>...</mo> </mtd> <mtd> <msub> <mi>&amp;sigma;</mi> <mrow> <mi>L</mi> <mi>L</mi> </mrow> </msub> </mtd> </mtr> </mtable> </mfenced> </mrow>

其中，σ_ij＝(1/N)∑_k(x_ik-m_i)(x_jk-m_j),k＝1,...,N，x_ik表示第k个像元第i个波段取值，m_i，m_j分别为第i个波段的均值和第j个波段的均值。

6.根据权利要求5所述的方法，其特征在于，步骤6中，通过设定置信度的方法设定判别阈值，将检测值大于判别阈值的像元标记为异常点，输出结果。

7.根据权利要求5所述的方法，其特征在于，步骤6中，设定一定量数目，将检测值较大的像元判为异常点，输出结果。