CN101807301B - 一种基于高阶统计量的高光谱图像目标检测方法 - Google Patents

Abstract

本发明一种基于高阶统计量的高光谱图像目标检测方法，它有四大步骤：步骤一：计算机在MATLAB R2008b环境下读取高光谱图像数据；步骤二：计算机对数据进行预处理，即对数据进行去均值化和白化处理；步骤三：构造检测滤波器在对目标光谱增益为1的约束下最小化输出数据的高阶统计量，求解检测滤波器最佳权向量；步骤四：设定合适的阈值，获取检测结果图像。本发明克服了现有技术的不足，充分利用数据的高阶统计量，取得了较好的检测效果，特别是能够提高低虚警率情况下的检测概率，它在高光谱遥感图像目标检测技术领域里具有实用价值和广阔的应用前景。

Description

一种基于高阶统计量的高光谱图像目标检测方法

(一)技术领域：

本发明涉及一种基于高阶统计量的高光谱图像目标检测方法，属于高光谱遥感图像目标检测技术领域。

(二)背景技术：

随着高光谱成像技术在过去三十年的迅猛发展，高光谱图像(Hyperspectral Images)的分析和处理方法成为当前国内外遥感图像领域的研究热点之一。高光谱图像和传统遥感图像相比，突出特点是光谱分辨率高，可获得观测对象的几十个甚至上百个光谱波段的图像信息。成像光谱***获得的连续波段宽度一般都在10nm以内，因此这种数据能以足够的光谱分辨率区分出那些具有诊断性光谱特征的地表物质，这也是高光谱技术用于探测甚至识别地物的物质基础(耿修瑞.高光谱遥感图像目标探测与分类技术研究[D].北京：中国科学院遥感应用研究所，2005)。高光谱图像数据是一个数据立方体，其一层图像对应一个光谱波段；每个像素点位置对应一条光谱曲线(贺霖，潘泉等.高光谱图像目标检测研究进展[J].电子学报，2009，37(9)：2016-2024)。因此，高光谱图像具有图谱合一的特性，即同时包含被观测物体的空间信息和光谱信息。

正是由于高光谱图像的高光谱分辨率，使得许多原来利用普通光学图像信息不能解决的问题，可以通过分析高光谱数据得到解决。高光谱图像目标检测是利用已知的目标光谱信息在高光谱图像中对感兴趣的目标进行检测、定位的技术(刘颖.高光谱图像目标识别新技术研究[D].哈尔滨：哈尔滨工业大学，2007)。高光谱图像目标检测技术在军事和民用领域中都有重要的应用价值。在军事领域可用于对飞机、坦克等军事目标进行检测、定位。在民用领域也可用于房屋等目标的检测。

目前存在的高光谱图像目标检测方法，主要是基于统计方法。其中经典的有匹配滤波器(MF，Matched Filter)，能量约束最小化(CEM，Constrained Energy Minimization)，自适应匹配滤波(AMF，Adaptive Matched Filter)，自适应一致性估计(ACE，Adaptive Coherence Estimator)等。MF假设目标像素和背景像素分别服从均值不同但协方差相同的高斯分布，然后利用Fisher线性判别准则构造检测算子。CEM在对目标光谱增益为1的约束下，最小化输出数据的能量，求解出检测算子。AMF假设背景服从零均值高斯分布，包含背景的像素也服从高斯分布，协方差和背景的协方差相同，但均值由目标子空间决定。ACE假设背景服从零均值高斯分布，包含背景的像素也服从高斯分布，协方差和背景的协方差结构相同但方差不同，均值由目标子空间决定。AMF和ACE都是基于广义似然比检验的方法，根据像素服从的概率分布的假设进行广义似然比检验，得出检测结果。

目前存在的这些目标检测方法主要利用数据的二阶统计量进行计算，主要涉及到了数据的协方差阵或相关阵，而没有利用数据的高阶统计量。而数据的高阶统计量往往包含着数据更重要的特征信息，对检测来说是非常有用的。在复杂背景下，目标所占的像素往往很少，高阶统计量能够更好的描述目标的统计特征。本发明针对上述情况，提出一种基于高阶统计量的目标检测方法，充分利用了数据的高阶统计信息，取得了较好的效果。

(三)发明内容：

1、目的：本发明的目的是提供一种基于高阶统计量的高光谱图像目标检测方法，它充分利用数据的高阶统计量，取得了较好的检测效果，特别是能够提高低虚警率情况下的检测概率，这点对于实际应用是十分有利的。

2、技术方案：本发明是通过以下技术方案实现的：

本发明一种基于高阶统计量的高光谱图像目标检测方法，它包括如下步骤：

步骤一：用计算机读取数据。计算机在MATLAB R2008b环境下读取高光谱图像数据，数据来源于成像光谱仪采集到的遥感图像，得到数据立方体。高光谱图像数据应去除被水汽吸收的波段和信噪比较低的波段。目标光谱来自已知的光谱库或对目标所在像素的光谱曲线取平均值得到。假设高光谱图像共有M个波段，像素的光谱曲线用向量形式表示x₀＝[x₀₁，x₀₂，...，x_0M]^T，x_0i为该像素第i个波段的值。已知的目标光谱也用向量形式表示s₀＝[s₀₁，s₀₂，...，s_0M]^T，s_0i为目标光谱第i个波段的值。

步骤二：数据预处理。用计算机读取数据后需要对数据进行预处理，即对数据进行去均值化和白化。

(1)去均值化。对高光谱图像去均值化，每个像素的光谱曲线都减去所有像素光谱曲线的平均值，去均值后使得整个高光谱图像数据的均值为零。另外，已知的目标光谱也要减去均值。

去均值化可通过下式进行：

x′＝x₀-E{x₀}

其中E{}表示取数学期望，可用图像各像素光谱曲线的平均值近似代替，即

其中x_0i表示第i个像素的光谱曲线，N等于图像的像素个数。相应的已知的目标光谱也需要进行相应的处理，减去图像各像素光谱曲线的期望即：

s′＝s₀-E{x₀}

(2)数据白化。对图像进行去均值化以后，还需要对高光谱图像数据进行白化。白化是去相关的过程，白化之后像素的光谱曲线各分量之间是不相关的，即像素光谱曲线的协方差矩阵是单位阵。另外，目标光谱也要进行相应的处理，即用白化矩阵乘以目标光谱。白化可通过下式进行：

x＝Γ^-1/2x′

其中Γ为数据的协方差矩阵，因为已进行了去均值化，协方差矩阵等于自相关阵，即

x_i′表示去均值化后的第i个像素的光谱曲线，T表示转置。相应的目标光谱也需要进行相应的处理即：

s＝Γ^-1/2s′

数据经过去均值化和白化处理后均值为零，协方差矩阵是单位阵即：

E{x}＝0

E{xx^T}＝I

步骤三：求解检测滤波器最佳权向量。数据进行去均值化和白化后，接下来的检测过程可以看作是一个滤波的过程。检测滤波器可以写为：

y＝w^Tx

x是每个像素的光谱曲线，y是滤波器的输出，w＝[w₁，w₂，...w_M]^T是滤波器的权向量。接下来就是要在对目标光谱增益为1的约束下，最小化输出数据的高阶统计量，求出最佳权向量。这样，最佳权向量w的求解问题可以写为：

$\left\{\begin{array}{c}\min E\left\{G\left(y\right)\right\}\\ s.t.w^{T}s=1\end{array}\right.$

其中，s为经过预处理后的目标光谱，y＝w^Tx，G为高阶统计量函数。本发明选取了四种高阶统计量，即：G(y)＝y³，G(y)＝y⁴，G(y)＝tanh(y)，G(y)＝sign(y)。

求解最佳权向量本质上是条件极值的求解问题，可用拉格朗日乘子法转化为无条件极值的问题，即求下列函数的极值：

J(w)＝E{G(w^Tx)}-λ(w^Ts-1)

等价于求解下列方程组：

$\left\{\begin{array}{c}\▿J=E\left\{G^{\′}\left(w^{T}x\right)x\right\}-\mathrm{\λs}=0\\ w^{T}s-1=0\end{array}\right.$

可从方程组解得λ＝E{(w^Tx)G′(w^Tx)}。本发明采用梯度下降法求解函数J(w)的极值，从而求解出最佳权向量w。具体步骤如下：

(1)、初始化w。w的初始值可以随机给定，然后进行归一化，即：

w⁺＝rand(w)

w＝w⁺/‖w⁺‖

(2)、迭代w。利用梯度下降法对w进行迭代，迭代规则为：

λ＝E{(w^Tx)G′(w^Tx)}

w⁺＝w-μ▽_J＝w-μ{E{G′(w^Tx)x}-λs}

w＝w⁺/‖w⁺‖

其中，μ为步长，本发明取μ＝10^-4。

(3)、停止迭代条件。当两次迭代的权向量w在相邻两次迭代的过程中变化不大时停止迭代，本发明中停止条件为：

‖w-w_old‖＜tol

w_old表示上一次迭代中w的值。本发明中取tol＝10^-4。

步骤四：获取检测结果图像。在求解出最佳权向量后，将每个像素的光谱曲线通过检测滤波器，得出输出数据。然后设定合适的阈值，滤波器输出大于阈值的判定该像素存在目标，小于阈值的判定该像素不存在目标。将输出大于阈值所对应的像素灰度值设为255，输出小于阈值所对应的像素灰度值设为0，就得到了检测结果的二值图像。在二值图像中目标对应的区域为白色，非目标对应的区域为黑色，从而完成了对目标的检测和定位。

3、优点及功效：本发明的优点是：针对目前高光谱图像目标检测方法没有利用数据的高阶统计量，只是利用了二阶统计量的情况，提供了一种基于高阶统计量的高光谱图像目标检测的新方法，充分利用了数据的高阶统计信息，可以获得较好的检测效果，特别是能够提高低虚警率条件下的检测概率，这点对于实际应用是十分有益的。因为在高虚警率条件下，目标会淹没在大片的虚警中，使得目标无法有效分离出来。

(四)附图说明：

图1本发明所述方法流程图

图2本发明所述方法检测过程

图3(a)本发明实验选用y³作为高阶统计量所得的接收机操作特性曲线

图3(b)本发明实验选用y⁴作为高阶统计量所得的接收机操作特性曲线

图3(c)本发明实验选用tanh(y)作为高阶统计量所得的接收机操作特性曲线

图3(d)本发明实验选用作sign(y)为高阶统计量所得的接收机操作特性曲线

本发明实验所用的数据来自AVIRIS传感器获取的美国圣地亚哥机场高光谱图像。实验图像从原图像中截取获得，实验图像从原图像的第110行、第125列所对应的像素点开始向右向下分别截取200个像素点，得到的实验图像共200×200个像素点。原图像共有224个波段，去除被水汽吸收的波段和信噪比较低的波段(1-6，33-35，97，107-113，153-166，221-224波段)后剩下189个波段。本发明中的实验主要针对图像中的飞机目标。

(五)具体实施方式：

为了更好地理解本发明的技术方案，以下结合附图对本发明的实施方式作进一步描述：

本发明在MATLAB R2008b语言环境下实现。计算机读取高光谱遥感图像数据后，获取的是数据立方体，首先对数据进行去均值化使得数据的均值为零，然后对数据进行白化，去除数据的相关性。检测过程可以看作是一个滤波的过程。去均值化和白化后的每个像素所对应的光谱曲线x＝[x₁，x₂，...，x_M]^T作为滤波器的输入，滤波器权向量w＝[w₁，w₂，...，w_M]^T和输入x的乘积w^Tx作为输出。设定输出数据的高阶统计量作为目标函数，寻找最佳权向量w，使得对已知目标光谱的增益为1的约束下，最小化输出数据的高阶统计量，而高阶统计量的选择可以有多种。这样w的求解问题实际上是一个条件极值的求解问题，可以用拉格朗日乘子法转化为无条件极值的求解问题，然后利用各种最优化的方法求得最佳权向量，本发明采用梯度下降法求取最佳权向量。求出最佳权向量后，将每个像素的光谱曲线通过滤波器得到输出数据，设定合适的阈值，判定输出值大于阈值所对应的像素存在目标，小于阈值所对应的像素不存在目标。

本发明提供一种基于高阶统计量的高光谱图像目标检测方法，该方法的流程图见图1所示，其计算机配置采用：Intel(R)Core(TM)2 Duo CPU E73002.66GHz。该检测方法包括如下步骤：

步骤一、用计算机读取数据。计算机在MATLAB R2008b环境下读取高光谱遥感图像数据。数据来源于成像光谱仪采集到的遥感图像，获取的是数据立方体。获取到的高光谱图像数据应去除被水汽吸收的波段和信噪比较低的波段。目标光谱可以来自已知的光谱库或对目标所在像素的光谱曲线取平均值得到。假设高光谱图像共有M个波段，像素的光谱曲线用向量形式表示x₀＝[x₀₁，x₀₂，...，x_0M]^T，x_0i为该像素第i个波段的值。已知的目标光谱也用向量形式表示s₀＝[s₀₁，s₀₂，...，s_0M]^T，s_0i为目标光谱第i个波段的值。

步骤二、数据预处理。得到遥感图像数据后，需要对数据进行预处理，即进行去均值化和白化。

(1)、去均值化。需要对高光谱图像数据进行去均值化，使得整个高光谱图像的均值为零。去均值化可通过下式进行：

x′＝x⁰-E{x₀}

其中E{}表示取数学期望，可用图像各像素光谱曲线的平均值近似代替，即

其中x_0i表示第i个像素的光谱曲线，N等于图像的像素个数。相应的已知的目标光谱也需要进行相应的处理，减去图像各像素光谱曲线的期望即：

s′＝s₀-E{x₀}

(2)、数据白化。进行去均值化之后还需要对数据进行白化处理，去除数据的相关性。白化可通过下式进行：

x＝Γ^-1/2x′

其中Γ为数据的协方差矩阵，因为已进行了去均值化，协方差矩阵等于自相关阵，即

x_i′表示去均值化后的第i个像素的光谱曲线，T表示转置。相应的目标光谱也需要进行相应的处理即：

s＝Γ^-1/2s′

数据经过去均值化和白化处理后均值为零，协方差矩阵是单位阵即：

E{x}＝0

E{xx^T}＝I

步骤三、求解检测滤波器最佳权向量。在进行数据预处理后，检测过程可以看作是一个滤波过程，如图2所示，检测滤波器的输入是经过去均值化和白化后的每个像素的光谱曲线x，输出y＝w^Tx。w＝[w₁，w₂，...w_M]^T为滤波器的权向量。若输出y＞η，则判定该像素存在目标，若输出y＜η，则判定该像素不存在目标。而滤波器的设计过程即为最佳权向量w的求解过程。

基于高阶统计量的检测方法的基本思想为：在对目标光谱增益为1的约束下，最小化输出数据的高阶统计量。用公式表示如下：

$\left\{\begin{array}{c}\min E\left\{G\left(y\right)\right\}\\ s.t.w^{T}s=1\end{array}\right.$

其中，s为经过预处理后的目标光谱，y＝w^Tx，G为高阶统计量函数。本发明选取了四种高阶统计量，即：G(y)＝y³，G(y)＝y⁴，G(y)＝tanh(y)，G(y)＝sign(y)。

求解最佳权向量本质上是条件极值的求解问题，可用拉格朗日乘子法转化为无条件极值的问题，即求下列函数的极值：

J(w)＝E{G(w^Tx)}-λ(w^Ts-1)

等价于求解下列方程组：

$\left\{\begin{array}{c}\▿J=E\left\{G^{\′}\left(w^{T}x\right)x\right\}-\mathrm{\λs}=0\\ w^{T}s-1=0\end{array}\right.$

可从方程组解得λ＝E{(w^Tx)G′(w^Tx)}。本发明采用梯度下降法求解函数J(w)的极值，从而求解出最佳权向量w。具体步骤如下：

(1)、初始化w。w的初始值可以随机给定，然后进行归一化，即：

w⁺＝rand(w)

w＝w⁺/‖w⁺‖

(2)、迭代w。利用梯度下降法对w进行迭代，迭代规则为：

λ＝E{(w^Tx)G′(w^Tx)}

w⁺＝w-μ▽_J＝w-μ{E{G′(w^Tx)x}-λs}

w＝w⁺/‖w⁺‖

其中，μ为步长，本发明取μ＝10^-4。

(3)、停止迭代条件。当两次迭代的权向量w在相邻两次迭代的过程中变化不大时停止迭代，本发明中停止条件为：

‖w-w_old‖＜tol

w_old表示上一次迭代中w的值。本发明中取tol＝10^-4。

步骤四、获取检测结果图像。在求得最佳权向量w后，如图2所示将每个像素的光谱曲线x输入检测滤波器，得到输出数据y＝w^Tx。设定合理的阈值η，若y＞η则判定该像素存在目标，若y＜η则判定该像素不存在目标。将判定存在目标的像素点灰度值设为255，判定不存在目标的像素点灰度值设为0，则可得到检测结果的二值图像，在二值图像中目标对应的区域为白色，非目标对应的区域为黑色，从而完成了对目标的检测和定位。

有益效果：

实验结果：为了验证本发明方法的有效性，我们使用该方法进行实验，并取得了较好的检测效果。本发明实验所用的数据来自AVIRIS传感器获取的美国圣地亚哥机场高光谱图像。实验图像从原图像中截取获得，实验图像从原图像的第110行、第125列所对应的像素点开始向右向下分别截取200个像素点，得到的实验图像共200×200个像素点。原图像共有224个波段，去除被水汽吸收的波段和信噪比较低的波段(1-6，33-35，97，107-113，153-166，221-224波段)后剩下189个波段。本发明中的实验主要针对图像中的飞机目标。

图3(a)、图3(b)、图3(c)、图3(d)是分别选用y³、y⁴、tanh(y)、sign(y)作为高阶统计量函数得到的检测结果的接收机操作特性曲线即ROC(Receiver Operation Curve)曲线。ROC曲线的横坐标是虚警率即非目标像素误检测为目标像素的概率，纵坐标是检测概率即目标像素正确识别为目标像素的概率，ROC曲线下的面积越大表明在同样的虚警率下取得的检测概率越大，说明检测器的性能越好。

从实验结果可以看出，我们发明的方法能取得较好的检测结果，特别是在低虚警率情况下检测概率较高，此方法能够应用到军用和民用领域的高光谱图像目标检测中，具有广阔的应用前景和价值。

Claims (1)

1.一种基于高阶统计量的高光谱图像目标检测方法，其特征在于：该方法具体步骤如下：

步骤一：用计算机读取数据；计算机在MATLAB R2008b环境下读取高光谱图像数据，数据来源于成像光谱仪采集到的遥感图像，得到数据立方体；高光谱图像数据应去除被水汽吸收的波段和信噪比较低的波段，目标光谱来自已知的光谱库或对目标所在像素的光谱曲线取平均值得到；假设高光谱图像共有M个波段，像素的光谱曲线用向量形式表示x₀＝[x₀₁，x₀₂，...，x_0M]^T，x_0i为该像素第i个波段的值，已知的目标光谱也用向量形式表示s₀＝[s₀₁，s₀₂，...，s_0M]^T，s_0i为目标光谱第i个波段的值；

步骤二：数据预处理；用计算机读取数据后需要对数据进行预处理，即对数据进行去均值化和白化；

(1)去均值化；对高光谱图像去均值化，每个像素的光谱曲线都减去所有像素光谱曲线的平均值，去均值后使得整个高光谱图像数据的均值为零；另外，已知的目标光谱也要减去均值；去均值化通过下式进行：

x′＝x₀-E{x₀}

其中：E{}表示取数学期望，用图像各像素光谱曲线的平均值近似代替，即 

x_0i表示第i个像素的光谱曲线，N等于图像的像素个数；相应的已知的目标光谱也需要进行相应的处理，减去图像各像素光谱曲线的期望，即：

s′＝s₀-E{x₀}

(2)数据白化；对图像进行去均值化以后，还需要对高光谱图像数据进行白化；白化是去相关的过程，白化之后像素的光谱曲线各分量之间是不相关的，即像素光谱曲线的协方差矩阵是单位阵；另外，目标光谱也要进行相应的处理，即用白化矩阵乘以目标光谱；

白化通过下式进行：

x＝Γ^-1/2x′

其中：Γ为数据的协方差矩阵，因为已进行了去均值化，协方差矩阵等于自相关阵，即 

x_i′表示去均值化后的第i个像素的光谱曲线，T表示转置；相应的目标光谱也需要进行相应的处理，即： 

s＝Γ^-1/2s′

数据经过去均值化和白化处理后均值为零，协方差矩阵是单位阵，即：

E{x}＝0

E{xx^T}＝I

步骤三：求解检测滤波器最佳权向量；数据进行去均值化和白化后，接下来的检测过程看作是一个滤波的过程，检测滤波器写为：

y＝w^Tx

其中，x是每个像素的光谱曲线，y是滤波器的输出，w＝[w₁，w₂，...w_M]^T是滤波器的权向量；

接下来就是要在对目标光谱增益为1的约束下，最小化输出数据的高阶统计量，求出最佳权向量；

这样，最佳权向量w的求解问题写为：

其中，s为经过预处理后的目标光谱，y＝w^Tx，G为高阶统计量函数；本发明选取了四种高阶统计量，即：G(y)＝y³，G(y)＝y⁴，G(y)＝tanh(y)，G(y)＝sign(y)；

求解最佳权向量本质上是条件极值的求解问题，用拉格朗日乘子法转化为无条件极值的问题，即求下列函数的极值：

J(w)＝E{G(w^Tx)}-λ(w^Ts-1)

等价于求解下列方程组：

从方程组解得λ＝E{(w^Tx)G′(w^Tx)}；本发明采用梯度下降法求解函数J(w)的极值，从而求解出最佳权向量w；其具体步骤如下：

(1)、初始化w；w的初始值随机给定，然后进行归一化，即：

w⁺＝rand(w)

w＝w⁺/||w⁺||

(2)、迭代w；利用梯度下降法对w进行迭代，迭代规则为： 

λ＝E{(w^Tx)G′(w^Tx)}

w＝w⁺||w⁺||

其中，μ为步长，本发明取μ＝10^-4；

(3)、停止迭代条件；当两次迭代的权向量w在相邻两次迭代的过程中变化不大时停止迭代，本发明中停止条件为：

||w-w_old||＜tol

w_old表示上一次迭代中w的值，本发明中取tol＝10^-4；

步骤四：获取检测结果图像；在求解出最佳权向量后，将每个像素的光谱曲线通过检测滤波器，得出输出数据；然后设定合适的阈值，滤波器输出大于阈值的判定该像素存在目标，小于阈值的判定该像素不存在目标；将输出大于阈值所对应的像素灰度值设为255，输出小于阈值所对应的像素灰度值设为0，就得到了检测结果的二值图像；在二值图像中目标对应的区域为白色，非目标对应的区域为黑色，从而完成了对目标的检测和定位。 

Images