CN101144861A - 高光谱亚像元目标探测方法及装置 - Google Patents

Abstract

本发明公开了一种高光谱亚像元目标探测方法，该方法包括：建立目标光谱和待测图像像元光谱二维矩阵的逆模型；获取该逆模型的回归系数向量；根据得到的回归系数向量获取各像元的马氏距离；判定马氏距离大于阀值的回归系数所对应的像元为亚像元目标点。本发明还提供了一种高光谱亚像元目标探测装置，采用本发明的亚像元目标探测方法及装置，不需要背景端元光谱的先验信息，对背景光谱的复杂度不敏感，具有较高的目标探测准确度和较快的运算速度。

Description

高光谱亚像元目标探测方法及装置

技术领域

本发明涉及高光谱遥感探测技术领域，尤其涉及一种高光谱亚像元目标探测方法及装置。

背景技术

高光谱遥感是二十世纪末地球观测系统中最重要的技术突破之一，它克服了传统单波段、多光谱遥感在波段数、波段范围、精细信息表达等方面的局限性，以较窄的波段区间、较多的波段数量提供遥感信息，能够从光谱空间中对地物予以细分和鉴别，在资源遥感、环境遥感、生态遥感等领域得到了广泛应用。高光谱遥感数据最主要的特点是将传统的图像维与光谱维信息融合为一体，在获取地表空间图像的同时，得到每个地物的连续光谱信息，该光谱信息能够反映出复杂背景下弱小目标和背景的细微差别，因此可借助丰富的光谱信息发现空间特征难以或无法探测的地面目标。由于目前高光谱遥感空间分辨率有限，单一像元很难只包括一种地物成分，其光谱往往是多种地物光谱的混叠，因此对于此类体积尺寸小于像元空间分辨率的目标探测称为亚像元目标探测。前述的空间分辨率即是指图像上所能辨别的地面物体最小尺寸。高光谱目标探测技术被广泛的应用于军事目标探测、矿物勘探、植被分布评估、环境污染检测等领域。

目前常见的亚像元目标探测方法主要包括：复合光谱最小二乘分离的探测方法、基于正交子空间投影(OSP，Orthogonal Subpixel Projection)的探测方法和基于广义似然比的探测方法等。

其中，复合光谱最小二乘分离的探测方法为：由于待测图像中的任一像元光谱均可以看成多种基本物质的光谱的加权和，权值即对应像元中各组成物质所占的混合比例，如果图像地物中各像元组分的纯光谱可从已知的光谱数据库中得到，则可以利用有约束最小二乘方法计算出目标对应的组分光谱在混合像元光谱中所占的比例，从而探测小于地面像元的目标。该方法可以对像元光谱进行定量分析，但前提是需要各像元组分光谱的先验信息，因此很难应用到背景未知的目标探测领域。

基于OSP的探测方法，主要是将像元光谱向量投影到背景特征的正交子空间，从而消除背景信息，突出目标光谱信息。现有技术中存在一种无监督的正交子空间投影方法对亚像元目标进行检测，该方法首先利用一种迭代的方法找到一组背景地物的端元光谱，也即地物中各组成物质的纯光谱，而后利用这些背景地物端元光谱构造检测算子，动态的对亚像元目标进行检测。该方法可以在背景地物端元未知的情况下构造出OSP检测算子，对背景噪声有一定的抑制作用，但由于计算所得到的背景地物端元光谱精确度不高，导致探测准确度不高，而且计算耗时较长。

现有技术中还存在一种基于广义似然比的探测方法，主要是先假设背景光谱信息满足某种多维分布，然后构造一定的检测算子通过假设检验的方法判定像元中是否含有潜在的目标光谱。这种方法能够给出理论上的虚警率，从而自适应地调节阈值，但该方法只能检测出背景中的奇异点，且需要一定的图像信息，训练样本的好坏对探测结果的影响较大。由于该探测方法需要有背景地物的先验信息，且对背景复杂度较为敏感，不能在背景变化复杂的小目标探测领域进行。

综上所述，现有技术中的亚像元目标探测方法大都需要背景地物端元光谱的先验信息，对背景光谱的复杂度较为敏感，而且目标探测的准确度较低、速度较慢。

发明内容

有鉴于此，本发明的主要目的在于提供一种高光谱亚像元目标探测方法及装置，以解决现有技术中的亚像元目标探测方法需要背景端元光谱先验信息，目标探测准确度较低和速度较慢的缺陷。

为达到上述目的，本发明的技术方案是这样实现的：

本发明一种高光谱亚像元目标探测方法，包括以下步骤：

建立目标光谱和待测图像像元光谱二维矩阵的逆模型；

获取所述逆模型的回归系数向量；

根据所述回归系数向量获取各像元的马氏距离；

判定马氏距离大于阀值的回归系数所对应的像元为亚像元目标点。

其中，所述建立目标光谱和待测图像像元光谱二维矩阵的逆模型，具体包括：

将获取的待测图像像元的三维高光谱数据表示为高光谱反射率的二维矩阵：R_m×n＝[p₁，p₂...p_y×i+j...p_x×y]，0<i≤x,0<j≤y，

或R_m×n＝[p₁，p₂...p_i+x×j...p_x×y]，0<i≤x,0<j≤y,

其中，R_n×m表示像元光谱的二维矩阵，[P₁,p₂...p_y×i+j...p_x×y]和[P₁,p₂...P_i+x×j...p_x×y]表示像元的光谱矢量，m表示波段数，n表示图像像元总数，x表示图像像元行数，y表示图像像元列数，n＝x×y；

建立所述目标光谱和所述待测图像像元光谱二维矩阵的线性关系：

S_I＝R_m×n×c_PLS+E_m×I，

其中，S_t表示目标光谱向量，R_n×m表示像元光谱的二维矩阵，c_PLS表示回归系数，E_m×1表示噪声向量。

其中，所述将三维高光谱数据表示为二维矩阵和建立线性关系之间，还包括：对所述目标光谱和所述待测图像像元光谱二维矩阵进行预处理。

其中，所述预处理包括标准正交变换处理或附加散射校正处理。

其中，所述标准正交变换处理具体包括：

$P_{\mathrm{hk},\mathrm{SNV}}=\frac{P_{\mathrm{hk}}-\stackrel{\&OverBar;}{P_h}}{\sqrt{\underset{k=1}{\overset{m}{\mathrm{\&Sigma;}}}{(P_{\mathrm{hk}}-\stackrel{\&OverBar;}{P_h})}^2}}{(m-1)}^2$

其中，p_hkSNV表示经过正交变换处理后待测图像中第h个像元在第k个波段的反射率值，

表示待测图像中第h个像元在各个波段处反射率的平均值，m表示波段数，m-1表示自由度；

所述附加散射校正处理具体包括：

计算平均光谱矢量： $\stackrel{-}{p}=\frac{1}{n}\underset{h=1}{\overset{n}{\mathrm{\&Sigma;}}}{p}_h;$

对每一个像元光谱进行线性回归： $p_h=m_h\stackrel{-}{p}+b_n;$

进行附加散射校正： $p_{h\left(\mathrm{MSC}\right)}=\frac{\left(p_h{-b}_h\right)}{m_h};$

其中， $\stackrel{\&OverBar;}{P}$ 表示平均光谱矢量， $\underset{h-1}{\overset{n}{\mathrm{\&Sigma;}}}{P}_h$ 表示对所有像元光谱矢量的求和，ph表示第h个像元的光谱矢量，m_h、b_h分别表示第h个像元光谱矢量p_h与所有像元平均光谱的线性回归的斜率与截距，p_h(MSC)表示经过附加散射校正后的像元光谱矢量。

其中，通过偏最小二乘迭代方法获取所述逆模型的回归系数向量，具体包括：

a、根据目标光谱向量获取初始权重向量： $W_n=R{S}_t^T,$ 其中，S_t表示目标光谱向量，w_n表示初始权重向量，R表示像元光谱的二维矩阵；

b、根据所述初始权重向量计算得分向量：t_n＝R^Tw_n，其中，t_n表示得分向量，R表示像元光谱的二维矩阵，w_n表示初始权重向量；

c、根据所述得分向量计算所述目标光谱的载荷向量：q_n＝S_tt_n，其中，q_n表示所述目标光谱的载荷向量，S_t表示目标光谱向量，t_n表示得分向量；

d、根据所述得分向量计算所述像元光谱二维矩阵的载荷向量：p_n＝Rt_n，其中，p_n表示所述像元光谱二维矩阵的载荷向量，R表示像元光谱的二维矩阵，t_n表示得分向量；

e、获取所述回归系数向量：c_PLS＝W(P^TW)^-1Q^T，其中，Q＝{q₁,q₂...q_n}，P＝{p₁,p₂...p_n}，  W＝{w₁,w₂....w_n}；

f、计算残差平方和： $\mathrm{SS}\left(n\right)=\underset{i=1}{\overset{m}{\mathrm{\&Sigma;}}}{(S_t-R_o{c}_{\mathrm{PLS}})}^2,$ 其中，SS(n)表示残差平方和，m表示波段数，R₀表示原始的像元光谱二维矩阵，S_t表示目标光谱向量，设定阀值为G，若SS(n-1)SS(n)≤G，则取该SS(n)对应的c_PLS为回归系数向量；否则，令S_t＝S_t-t_n q_n，R＝R-t_n p_n，返回步骤a重复上述操作，直到SS(n-1)-SS(n)≤G，然后获取回归系数向量c_PLS。

本发明还提供了一种高光谱亚像元目标探测装置，包括：

模型建立单元，用于建立目标光谱和待测图像像元光谱二维矩阵的逆模型；

回归系数向量获取单元，用于获取所述逆模型的回归系数向量；

马氏距离获取单元，用于根据所述回归系数向量获取各像元的马氏距离；

判定单元，用于判定马氏距离大于阀值的回归系数所对应的像元为亚像元目标点。

其中，所述模型建立单元包括：

矩阵生成子单元，用于将待测图像像元的三维高光谱数据表示为高光谱反射率的二维矩阵；

线性关系建立子单元，用于建立所述目标光谱和所述待测图像像元光谱二维矩阵的线性关系。

其中，所述模型建立单元还包括：预处理子单元，用于对所述目标光谱和所述矩阵生成子单元生成的待测图像像元光谱二维矩阵进行预处理并提供给线性关系建立子单元。

其中，所述装置还包括：数据获取单元，用于获取待测图像像元的三维高光谱数据并提供给模型建立单元。

与现有技术相比，本发明具有以下优点：

本发明所用到的逆模型不需要背景端元光谱的先验信息，仅需要先验的目标光谱信息，因此对背景光谱的复杂度不敏感；探测准确度仅仅依赖于目标光谱的精确程度，从而目标探测准确度高；又由于偏最小二乘中的迭代计算方法仅需要少量矩阵求逆运算，且在检测部分采用了速度较快的马氏奇异值检测方法，使得本发明的探测方法具有较快的运算速度。

附图说明

图1为本发明一种高光谱亚像元目标探测方法的流程图；

图2为本发明实施例的原始高光谱图像的一示意图；

图3为本发明实施例的原始高光谱图像的另一示意图；

图4为本发明实施例的目标光谱示意图；

图5为本发明实施例的混入目标光谱后的像元点示意图；

图6为本发明实施例的回归系数三维示意图；

图7为本发明实施例的马氏距离示意图；

图8为本发明实施例的亚像元目标探测结果示意图；

图9为本发明一种高光谱亚像元目标探测装置的组成结构示意图。

具体实施方式

下面结合附图和具体实施例对本发明的技术方案做进一步详细阐述：

如图1所示，图1为本发明一种高光谱亚像元目标探测方法的流程图，主要包括以下步骤：

步骤101，建立目标光谱和待测图像像元光谱二维矩阵的逆模型。

首先，将获取到待测图像像元的三维高光谱数据表示为高光谱反射率的二维矩阵。该待测图像像元的三维高光谱数据是由成像光谱仪获得，成像光谱仪在每个波段可获得一幅感光图像，而由一些列连续的波段的图像叠加而成的立方体即为高光谱图像，该高光谱图像包括图像空间的两维信息和光谱波段方向的一维信息，因此为三维的高光谱数据。获取待测图像像元的三维高光谱数据，并将该三维的高光谱数据表示为如下的二维矩阵：

R_m×n＝[p₁，p₂...p_y×i+j...p_x×y]，0<i≤x,0<j y    (1)

或R_m×n＝[p₁，p₂...p_i+x×j...p_x×y]，0<i≤x,0<j≤y    (2)

其中，(1)式为待测图像像元按行展开的表示形式，(2)式为待测图像像元按列展开的表示形式，R_n×m表示像元光谱的二维矩阵，x表示高光谱图像像元行数，y表示高光谱图像像元列数，m表示波段数，n表示图像像元总数，n＝x×y。p₁，p₂...p_y×i+j...p_x×y]和[p₁，p₂...p_i+x×j.．.p_x×y]表示高光谱图像像元的光谱矢量，在按行展开的R_n×m中，p_y×i+j表示待测图像中第i行第j列对应像元的光谱矢量；同理，在按列展开的R_n×m中，p_y×i+j则表示待测图像中第i行第j列对应像元的光谱矢量。每个像元的光谱矢量包括该像元在各个波段处的反射率值，例如：假设第h个像元的光谱矢量为p_h，则p_h＝[p_h1,p_h2...p_hk...p_hm]^T，  [p_h1,p_h2...p_hk...p_hm]^T代表[p_h1,p_h2...p_hk...p_hm]的转置矩阵，其中p_hk表示该第h个像元在第k个波段处的反射率值。由于成像光谱仪获取到图像像元点在每个波段处的值即是该像元点在该波段处的反射率值，因此该反射率值为已知量。

然后，对目标光谱和待测图像像元光谱二维矩阵进行预处理，以校正因大气散射引起的光谱误差。该预处理可为标准正交变换处理，也可为附加散射校正处理。当然，本发明的预处理方法并不仅仅局限于上述两种处理方法，其他任何能校正因大气散射引起的光谱误差的处理方法也应属于本发明的保护范围。

标准正交变换处理的公式如下：

$P_{\mathrm{hk},\mathrm{SNV}}=\frac{P_{\mathrm{hk}}-{\stackrel{-}{P}}_h}{{\sqrt{\underset{k=1}{\overset{m}{\mathrm{\&Sigma;}}}(P_{\mathrm{hk}}-\stackrel{-}{P_h})}}^2}{(m-1)}^{1/2}$

其中，p_hk,SNF表示经过正交变换处理后待测图像中第h个像元在第k个波段的反射率值，

表示待测图像中第h个像元在各个波段处反射率的平均值，m表示波段数，m-1表示自由度。

附加散射校正处理的过程包括：

首先，计算平均光谱矢量：

$\stackrel{\&OverBar;}{P}=\frac{1}{n}\underset{h=1}{\overset{n}{\mathrm{\&Sigma;}}}{P}_h$

然后，对每一个像元光谱进行线性回归：

$P_h=m_h\stackrel{\&OverBar;}{P}+b_n$

再进行附加散射校正：

$P_{h\left(\mathrm{MSC}\right)}=\frac{(P_h-b_n)}{m_h}$

上述(4)到(6)式中，

表示平均光谱矢量，

表示对所有像元光谱矢量的求和，p_h表示第h个像元的光谱矢量，m_h、b_h分别表示第h个像元光谱矢量p_h与所有像元平均光谱的线性回归的斜率与截距，p_h(MSC)表示经过附加散射校正后的像元光谱矢量。

分别对目标光谱和待测图像中各像元光谱进行预处理之后，则可得到消除大气散射误差之后较精确的光谱信息。

最后，建立目标光谱和待测图像像元光谱二维矩阵的线性关系，该线性关系表示如下：

S_t＝R_m×n×c_PLS+E_m×1    (7)

其中，S₁表示目标光谱向量，R_n×m表示像元光谱的二维矩阵，c_PLS表示回归系数向量，c_PLS值的大小可以反映目标光谱对各像元光谱的贡献程度，c_PLS值越大，则表明目标光谱对像元光谱的贡献越大，E_m×1表示噪声向量，为一维噪声矩阵。

步骤102，获取该逆模型的回归系数向量。

由于R_n×m的秩小于波段数，R_n×m中各像元光谱存在高度相关性，采用通常的最小二乘方法对回归系数进行c_PLS估计时需要对R_m×n ^TR_m×n进行求逆，R^T _m×n代表R_m×n的转置矩阵，而当R_n×m中的变量高度相关时，行列式|R^T _m×nR_m×n|几乎接近于零，对R_m×n ^TR_m×n求逆会产生严重的舍入误差，因此采用通常的最小二乘方法求得的c_PLS会受到舍入误差的较大影响。本发明采用单因变量偏最小二乘方法求取该逆模型的回归系数向量，具体包括：

a、根据目标光谱向量计算初始权重向量：

$w_s=R{S}_t^T$

其中，S_t表示目标光谱向量，S_t为已知量，S₁ ^T代表S₁的转置矩阵，w_n表示初始权重向量，R表示像元光谱的二维矩阵。

b、根据初始权重向量计算得分向量：

t_n＝R^Tw_n    (9)

其中，t_n表示得分向量，R表示像元光谱的二维矩阵，R^T代表R的转置矩阵，w_n表示初始权重向量。

c、根据得分向量计算目标光谱的载荷向量：

q_n＝S_tt_n    (10)

其中，S_t表示目标光谱向量，q_n表示载荷向量，t_n表示得分向量。

d、根据得分向量计算像元光谱二维矩阵的载荷向量：

p_n＝Rt_n    (11)

其中，p_n表示像元光谱二维矩阵的载荷向量，R表示像元光谱的二维矩阵，t_n表示得分向量。

e、计算回归系数向量：

c_PLS＝W(P^TW)^-1 Q^T    (12)

其中，Q＝{q₁,q₂...q_n}，P＝{p₁,p₂...p_n}，W＝{w₁,w₂....w_n}，(P^TW)^-1代表对(P^TW)进行求逆运算。

f、计算残差平方和：

$\mathrm{SS}\left(n\right)=\underset{i=1}{\overset{m}{\mathrm{\&Sigma;}}}{(S_t-R_o{c}_{\mathrm{PLS}})}^2$

其中，SS(n)表示残差平方和，m表示波段数，R₀表示原始的像元光谱二维矩阵，若SS(n-1)-SS(n)≤G，则表明收敛，取该SS(n)对应的c_PLS为回归系数向量；否则，令St＝St-tnqn，  R＝R-tnpn，返回步骤a重复上述操作，直到SS(n-1)-SS(n)≤G，然后获取回归系数向量c_PLS。在实际应用中，G值可根据实际需要进行选取。

步骤103，根据得到的回归系数向量获取各像元的马氏距离。

根据得到的回归系数向量并利用下式计算各像元的马氏距离，公式如下：

$d_h\left(c_h\right)=\sqrt{\frac{{(c_h-\stackrel{\&OverBar;}{c})}^2}{\mathrm{std}{\left(c\right)}^2}},h\&le;n$

其中，c_h表示回归系数向量中对应第h个像元的回归系数值，n表示像元总数，std(c)表示回归系数的标准差。d_h(c_h)表示第h个像元对应的回归系数的马氏距离。前述std(c)可通过下式求得：

$\mathrm{std}\left(c\right)=\sqrt{\frac{1}{n-1}\underset{h=1}{\overset{n}{\mathrm{\&Sigma;}}}{(c_h-\stackrel{\&OverBar;}{c})}^2}---\left(15\right)$

其中，std(c)表示回归系数的标准差，c_h表示第h个像元对应的回归系数值，

表示所有像元回归系数的平均值， $\stackrel{\&OverBar;}{c}=\underset{h=1}{\overset{n}{\mathrm{\&Sigma;}}}{c}_h.$

步骤104，判定马氏距离大于阀值的回归系数所对应的像元为亚像元目标点。

计算出各像元回归系数所对应的马氏距离后，检测阀值可根据如下方式来快速获得：即可以认为回归系数小于0的像元必为背景像元，则该些已知的背景像元必存在一个马氏距离最大的像元点，可取该像元所对应的马氏距离为背景与目标的分割阈值，即马氏距离检测阈值。

下面结合具体实施例对上述本发明的高光谱亚像元目标探测方法进一步详细阐述。本实例所用的高光谱遥感数据来源于机载成像光谱仪，机载成像光谱仪为采用推扫成像方式的成像光谱仪，在0.4微米-2.45微米的波长范围获取224个波长处的空间图像信息，波长间隔为10纳米，当飞机在20千米高空飞行时，图像空间分辨率可达20×20米。本实施例使用的高光谱图像如图2所示，大小为614×512像元，每个像元光谱包括224个波长，波长范围从369.85纳米到2506.81纳米。具体探测过程如下：

A、取图2中所示白色方框内64×64像元的图像为本实施例的待测图像，图3则为该待测图像的放大示意图，待测图像像元总数n＝64×64＝4096。剔除224个波段中的1-6、33、107-114、153-168、222-224波段等35个信噪比较低的坏波段，其余190波段为本实施例的使用波段，也即波段数m＝190。取图2右下角黑色圆圈内箭头所指的反射能量值较高的屋顶光谱为本实施例的目标光谱，该目标光谱的光谱图如图4所示。提取64×64图像中坐标为(10，32)、(10,42)、(10，52)、(32，32)、(32，42)、(32，52)、(42，32)、(42，42)、(42，52)的九个像元点，分别混入5％的目标光谱，混入目标光谱后的像元点如图5所示。如果将图像中的像元按行进行展开，则(10，32)的像元点对应的光谱矢量即为p_64×10+32，其他像元点也同理，在此不再一一描述。

需要指出的是，本发明实施例的待测图像像元点是任意选取的，选取像元点的数量也是任意的，此处同时选取九个像元点是为了表明本发明的探测方法可同时对多个像元点进行探测。对前述的九个像元点混入目标光谱后，该九个像元点中即存在了目标光谱的信息，则通过本发明的探测方法即可将存在目标光谱的该九个点都探测出来，这也是本发明实施例的最终目的。

B、将所获取的九个像元点的三维高光谱数据表示为二维的高光谱反射率矩阵，表示方法与前述相同，在此不再多述。

C、对目标光谱和像元光谱二维矩阵分别进行预处理，处理方法与前述相同，在此也不再多述。然后，建立目标光谱和像元光谱之间的逆模型。

D、利用单因变量偏最小二乘方法对上述逆模型的回归系数向量c_PLS进行求解，在该偏最小二乘方法中设定阀值为10^-3。求取回归系数的三维图如图6所示，图中的X、Y轴分别代表待测图像的横、纵坐标，Z轴代表回归系数的值。从图6中可看出，混有5％目标光谱的像元点与未混入目标光谱的背景像元点的回归系数值差异较大，可见回归系数的大小可以反映目标光谱对像元光谱的贡献程度。

E、计算各像元点回归系数的马氏距离。计算结果如图7所示，图中的横轴代表像元数，纵轴代表马氏距离值。由上述阈值选取方法可得检验阀值为5.1402，从图7中可看出，所选九个像元点回归系数的马氏距离都大于设定阈值。因此，检测结果如图8所示，所选九个像元点都被判定为亚像元目标点，检测结果与实际情况相符。

本发明还提供了一种高光谱亚像元目标探测装置，如图9所示，该装置包括：数据获取单元100、模型建立单元200、回归系数向量获取单元300、马氏距离获取单元400和判定单元500。其中，数据获取单元100，用于获取待测图像像元的三维高光谱数据。模型建立单元200，连接数据获取单元100，用于根据数据获取单元100获取的像元高光谱数据建立目标光谱和待测图像像元光谱的逆模型。回归系数向量获取单元300，连接模型建立单元200，用于获取逆模型的回归系数向量。马氏距离获取单元400，连接回归系数向量获取单元300，用于根据得到的回归系数向量获取各像元的马氏距离。判定单元500，连接马氏距离获取单元400，用于判定马氏距离大于阀值的回归系数所对应的像元为亚像元目标点。

其中，模型建立单元200还包括：矩阵生成子单元210、预处理子单元220和线性关系建立子单元230。矩阵生成子单元210，用于将待测图像像元的三维高光谱数据表示为高光谱反射率的二维矩阵。预处理子单元220，连接矩阵生成子单元210，用于对目标光谱和矩阵生成子单元210生成的待测图像像元光谱二维矩阵进行预处理并提供给线性关系建立子单元230。线性关系建立子单元230，连接预处理子单元220，用于建立目标光谱和待测像元光谱二维矩阵的线性关系。

综上所述，本发明一种高光谱亚像元目标探测的方法及装置，所用到的逆模型不需要背景端元光谱的先验信息，仅需要先验的目标光谱信息，因此对背景光谱的复杂度不敏感；探测准确度仅仅依赖于目标光谱的精确程度，从而目标探测准确度高；又由于偏最小二乘中的迭代计算方法仅需要少量矩阵求逆运算，且在检测部分采用了速度较快的马氏奇异值检测方法，使得本发明的探测方法具有较快的运算速度。

以上所述，仅为本发明的较佳实施例，并非用于限定本发明的保护范围。

Claims (10)

1.一种高光谱亚像元目标探测方法，其特征在于，包括以下步骤：

建立目标光谱和待测图像像元光谱二维矩阵的逆模型；

获取所述逆模型的回归系数向量；

根据所述回归系数向量获取各像元的马氏距离；

判定马氏距离大于阀值的回归系数所对应的像元为亚像元目标点。

2.如权利要求1所述高光谱亚像元目标探测方法，其特征在于，所述建立目标光谱和待测图像像元光谱二维矩阵的逆模型，具体包括：

将获取的待测图像像元的三维高光谱数据表示为高光谱反射率的二维矩阵：    R_m×n＝[p₁，p₂…p_y×i+j…p_x×y]，0＜i≤x，0＜j≤y，

    或R_m×n＝[p₁，p₂…p_i+x×j…p_x×y]，0＜i≤x，0＜j≤y，

其中，R_n×m表示像元光谱的二维矩阵，[p₁，p₂…p_y×i+j…p_x×y]和[p₁，p₂…p_i+x×j…p_x×y]表示像元的光谱矢量，m表示波段数，n表示图像像元总数，x表示图像像元行数，y表示图像像元列数，n＝x×y；

建立所述目标光谱和所述待测图像像元光谱二维矩阵的线性关系：

S_t＝R_m×n×c_PLS+E_m×1，

其中，S_t表示目标光谱向量，R_n×m表示像元光谱的二维矩阵，c_PLS表示回归系数，E_m×1表示噪声向量。

3.如权利要求2所述高光谱亚像元目标探测方法，其特征在于，所述将三维高光谱数据表示为二维矩阵和建立线性关系之间，还包括：对所述目标光谱和所述待测图像像元光谱二维矩阵进行预处理。

4.如权利要求3所述高光谱亚像元目标探测方法，其特征在于，所述预处理包括标准正交变换处理或附加散射校正处理。

5.如权利要求4所述高光谱亚像元目标探测方法，其特征在于，所述标准正交变换处理具体包括：

$p_{\mathrm{hk},\mathrm{SNV}}=\frac{p_{\mathrm{hk}}-{\stackrel{\&OverBar;}{p}}_h}{\sqrt{\underset{k=1}{\overset{m}{\mathrm{\&Sigma;}}}{(p_{\mathrm{hk}}-{\stackrel{\&OverBar;}{p}}_h)}^2}}{(m-1)}^{1/2}$

其中，p_hk，SNV表示经过正交变换处理后待测图像中第h个像元在第k个波段的反射率值，

表示待测图像中第h个像元在各个波段处反射率的平均值，m表示波段数，m-1表示自由度；

所述附加散射校正处理具体包括：

计算平均光谱矢量： $\stackrel{\&OverBar;}{p}=\frac{1}{n}\underset{h=1}{\overset{n}{\mathrm{\&Sigma;}}}{p}_h;$

对每一个像元光谱进行线性回归： $p_h=m_h\stackrel{\&OverBar;}{p}+b_h;$

进行附加散射校正： $p_{h\left(\mathrm{MSC}\right)}=\frac{(p_h-b_h)}{m_h};$

其中，

表示平均光谱矢量，

表示对所有像元光谱矢量的求和，p_h表示第h个像元的光谱矢量，m_h、b_h分别表示第h个像元光谱矢量p_h与所有像元平均光谱的线性回归的斜率与截距，p_h(MSC)表示经过附加散射校正后的像元光谱矢量。

6.如权利要求1所述高光谱亚像元目标探测方法，其特征在于，通过偏最小二乘迭代方法获取所述逆模型的回归系数向量，具体包括：

a、根据目标光谱向量获取初始权重向量： $w_n={\mathrm{RS}}_t^T,$ 其中，S_t表示目标光谱向量，w_n表示初始权重向量，R表示像元光谱的二维矩阵；

b、根据所述初始权重向量计算得分向量：t_n＝R^Tw_n，其中，t_n表示得分向量，R表示像元光谱的二维矩阵，w_n表示初始权重向量；

c、根据所述得分向量计算所述目标光谱的载荷向量：q_n＝S_tt_n，其中，q_n表示所述目标光谱的载荷向量，S_t表示目标光谱向量，t_n表示得分向量；

d、根据所述得分向量计算所述像元光谱二维矩阵的载荷向量：p_n＝Rt_n，其中，p_n表示所述像元光谱二维矩阵的载荷向量，R表示像元光谱的二维矩阵，t_n表示得分向量；

e、获取所述回归系数向量：c_PLS＝W(P^TW)^-1Q^T，其中，Q＝{q₁，q₂…q_n}，P＝{p₁，p₂…p_n}，W＝{w₁，w₂…w_n}；

f、计算残差平方和： $\mathrm{SS}\left(n\right)=\underset{i=1}{\overset{m}{\mathrm{\&Sigma;}}}{(S_t-R_0{c}_{\mathrm{PLS}})}^2,$ 其中，SS(n)表示残差平方和，m表示波段数，R₀表示原始的像元光谱二维矩阵，S_t表示目标光谱向量，设定阀值为G，若SS(n-1)-SS(n)≤G，则取该SS(n)对应的c_PLS为回归系数向量；否则，令S_t＝S_t-t_nq_n，R＝R-t_np_n，返回步骤a重复上述操作，直到SS(n-1)-SS(n)≤G，然后获取回归系数向量c_PLS。

7.一种高光谱亚像元目标探测装置，其特征在于，包括：

模型建立单元，用于建立目标光谱和待测图像像元光谱二维矩阵的逆模型；

回归系数向量获取单元，用于获取所述逆模型的回归系数向量；

马氏距离获取单元，用于根据所述回归系数向量获取各像元的马氏距离；

判定单元，用于判定马氏距离大于阀值的回归系数所对应的像元为亚像元目标点。

8.如权利要求7所述高光谱亚像元目标探测装置，其特征在于，所述模型建立单元包括：

矩阵生成子单元，用于将待测图像像元的三维高光谱数据表示为高光谱反射率的二维矩阵；

线性关系建立子单元，用于建立所述目标光谱和所述待测图像像元光谱二维矩阵的线性关系。

9.如权利要求8所述高光谱亚像元目标探测装置，其特征在于，所述模型建立单元还包括：预处理子单元，用于对所述目标光谱和所述矩阵生成子单元生成的待测图像像元光谱二维矩阵进行预处理并提供给线性关系建立子单元。

10.如权利要求7至9中任一项所述高光谱亚像元目标探测装置，其特征在于，所述装置还包括：数据获取单元，用于获取待测图像像元的三维高光谱数据并提供给模型建立单元。

Images