CN115937302A - 结合边缘保持的高光谱图像亚像元定位方法 - Google Patents

Abstract

本发明涉及结合边缘保持的高光谱图像亚像元定位方法，包括：S1、利用高斯模糊核和原始高光谱图像进行非盲去卷积来达到削减原始高光谱图像中点扩展函数效应；S2、使用域变换递归滤波器对步骤S1的结果进行滤波来起到边缘保持和减少噪声的作用；S3、通过线性光谱解混技术获取丰度图像；S4、对丰度图像使用插值方法进行上采样以获得各类别亚像元软类值；S5、利用类别分配策略为各亚像元分配类别标签，得到最终的亚像元定位结果。本发明能够处理有效地处理点扩展函数效应以及图像中的纹理细节等噪声，从而使亚像元定位精度得以提高。

Description

结合边缘保持的高光谱图像亚像元定位方法

技术领域

本发明涉及一种结合边缘保持的高光谱图像亚像元定位方法，属于高光谱图像亚像元定位技术领域。

背景技术

由于受到高光谱传感器瞬时视场的制约而使得高光谱遥感图像空间分辨率低，以及地面真实情况的复杂性，这导致混合像元大量存在，传统的分类方法以像元为单位来为其赋予类别标签，该种分类往往不能反映地面的真实情况。并且光谱解混技术是用来获取混合像元内端元种类以及各端元对应的丰度，但是无法明确各种地物在空间上的分布。而亚像元定位(subpixel mapping,SPM)的目的就是通过一定的算法或模型来确定每个像元内不同地物的空间位置，传统的两步SPM直接对原始高光谱图像进行处理，首先将原始高光谱图像进行光谱解混得到各类别丰度矩阵，然后按照某一尺度分别对各类别丰度矩阵进行上采样依次得到亚像尺度下的各类别软类值，最后使用类别分配策略获取最终的亚像元定位结果。然而，两步SPM没有考虑到高光谱遥感成像系统的点扩展函数效应带给原始高光谱图像的模糊，以及远距离测量时不可避免的会受到云雾等所带来的无关纹理细节对亚像元定位精度的影响,为了解决该问题，本文提出了一种结合边缘保持(edge-preserving，EP)的高光谱亚像元定位方法。

发明内容

为了克服现有研究的不足，本发明提供了一种结合边缘保持的高光谱图像亚像元定位方法，用来缓解原始高光谱易受到来自于成像系统的点扩展函数效应的影响的现象，同时用来减少原始高光谱图像上无关纹理细节。

一种结合边缘保持的高光谱图像亚像元定位方法的具体步骤如下：

S1：利用高斯模糊核和原始高光谱图像进行非盲去卷积来削减原始高光谱图像中点扩展函数效应；

S2：使用域变换递归滤波器滤波来起到边缘保持和减少噪声的作用；

S3：通过光谱解混技术得到丰度图像；

S4：对丰度图像进行上采样以获得各类别亚像元软类值；

S5：利用类别分配策略为各亚像元分配标签，得到亚像元定位结果。

作为优选方案，在步骤S1中，受到点扩展函数效应影响的高光谱成像模型表示为：

B＝I*K+N

其中，B∈R^m×n×b为高光谱遥感传感器获取得到的模糊图像，K∈R^k×k为点扩展函数所对应的模糊核，I∈R^m×n×b为所求清晰图像，N∈R^m×n×b表示加性噪声，*表示卷积运算，其中m×n表示空间分辨率，b表示波段数，k为点扩展函数对应的卷积核阶数。削弱点扩展函数导致对高光谱图像的影响的过程就是利用模糊图像和模糊核来得到清晰图像。

设B＝(B₁,B₂,...,B_b)^T，I＝(I₁,I₂,...,I_b)^T，i＝1,2,...,b，(·)^T表示转置。假设B_i满足泊松噪声分布，则有：

B_i＝Poisson((I_i*K)(x))

其中，该式表示成像系统获取到的高光谱图像与待求清晰图像的各波段灰度图用卷积核K卷积后的结果满足泊松分布。则可将待求清晰图像I的似然概率函数表示为：

其中，！表示阶乘运算，K的表示如下：

其中，σ表示高斯函数标准差，V(i,j)表示以(i,j)为中心的局部窗口的空间范围，对I的似然函数取对数得到其能量函数，结果如下；

E(I_i)＝∑{(I_i*K)-B_iln(I_i*K)}

即求其最大似然解

为：

令点扩展函数对应归一化的模糊核K满足∫K(x)dx＝1，则对上式求导得到迭代式：

其中，K^*为K的伴随矩阵，t为迭代次数，一般选择I_i＝B_i作为初始条件进行迭代，通过设置合适的迭代次数进而得到所求图像近似值I。

作为优选方案，在步骤S2中，非盲去卷积得到的清晰图像I使用域变换递归滤波器进行滤波，具体实现为：

将经过步骤S1获得的清晰高光谱图像按空间维度展开，即I＝(I₁,I₂,...,I_b)^T，i＝1,2,...,b，将任意波段I_i看成二维灰度图，将其分别沿着图像的水平和垂直方向转换为多组一维信号。对于任意给定一维信号S，其域变换可被定义为：

其中，S₀表示一维信号初始时刻的值，δ_s和δ_r为控制滤波器平滑度的两个常数，S_j-1和S_j表示第j-1时刻和第j时刻一维信号的值，·表示绝对值，U_i表示域变换后的信号强度。则递归保边滤波器可以在变换域中定义为：

J[n]＝(1-α^d)S[n]+α^dJ[n-1]

其中，

为反馈系数，d＝U_n-U_n-1表示变换域中相邻两信号间的距离，S_n表示第n时刻一维信号的强度值，J[n-1]和J[n]分别表示第n-1和第n时刻一维信号域变换递归滤波的值。对于二维图像按行展开成一维信号执行域变换递归滤波后重新组合成二维信号，再将上一步得到的二维信号按列展开成一维信号执行域变换递归滤波后重新组合成二维信号，按照该步骤递归三次得到滤波后的结果。

作为优选方案，在步骤S3中，光谱解混技术具体为：利用基于线性解混模型的光谱解混方法来处理域变换递归滤波器滤波后的图像，以获得各类别丰度图像。线性解混模型为：线性解混模型是指混合像元内代表地物类别的端元与其所占的比例成线性关系，则设

为待解混的高光谱数据，d为波段数，n为像元数，

为端元矩阵，p为端元个数，

为丰度系数矩阵，

为噪声矩阵，则有，

其中，

表示1行p列的全1向量，

表示1行n列的全1向量，两个约束条件分别代表丰度系数矩阵的非负约束和每个像元丰度系数和为1约束。

作为优选方案，在步骤S4中，上采样方法使用的径向基函数插值，通过径向基函数插值对步骤S3获得的各地物类别丰度图像进行处理，得到每个亚像元的各类别的软类值。

作为优选方案，在步骤S5中,类别分配策略选用的是以类为单元的类别分配策略(UOC)，具体为：假设场景中含有C类地物，P_j表示第j个像元，F_k(P_j)表示第j个像元中第k类地物的丰度系数，则有，

s.t.:F_k(P_j)≥0

其中，该式表示步骤S3得到的每个像元的丰度系数和为一，约束条件表示该像元各地物对应丰度系数非负。则该像元通过尺度因子s进行上采样使得该像元划分为s²个亚像元，则该像元内第k类地物所占亚像元个数N_k(P_j)可通过下式获取：

N_k(P_j)＝round[F_k(P_j)·s²]

其中，round[·]表示四舍五入运算。即在类别分配时，必须满足每个亚像元仅仅被分配为一类地物，并且每个像元内各类地物对应亚像元的数量与丰度系数成比例。

最后借助步骤S4获取到的各亚像元软类值，并依据衡量空间相关性的全局莫兰指数从高到低来确定地物类别分配次序，然后结合各像元内每类亚像元数量逐像元地按照类别分配次序和软类值中从大到小来为各亚像元赋予类别，分配结束后即可得到最终的亚像元定位结果。

与现有技术相比，本发明的有益效果在于：

传统的两步亚像元定位方法没有考虑点扩展函数效应对于亚像元定位精度的制约，直接使用未处理的原始高光谱图像进行亚像元定位，而本发明使用非盲去卷积来处理原始高光谱图像，可以有效的缓解点扩展函数效应对亚像元定位精度的影响。针对高光谱图像中的无关纹理细节(如伪影等)，本发明使用域变换递归滤波器进行空间维滤波来起到对无关纹理细节的减少，且可以起到到边缘保持作用，从而提高亚像元定位精度。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。

图1是本发明实施例的一种结合边缘保持的高光谱图像亚像元定位方法的流程图。

图2是本发明实施例的实验部分数据集的高分辨率灰度图像与其降采样得到的低分辨率灰度图像，(a)高分辨率灰度图像，(b)低分辨率灰度图像。

图3是本发明实施例的实验部分参考分类图与使用本发明方法获得的亚像元定位结果图，(a)参考分类图，(b)定位结果图。

具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

实施例：

本发明提出的一种结合边缘保持的高光谱图像亚像元定位方法，其实施例流程图如图1所示。针对原始低分辨率高光谱图像B∈R^m×n×b，

表示实数域，m×n表示空间分辨率，b表示波段数，设尺度因子为s，则最后亚像元定位获取到的结果图的空间分辨率为ms×ns。

本实施例方法具体步骤如下：

S1：利用高斯模糊核和原始高光谱图像进行非盲去卷积来削弱原始高光谱图像中点扩展函数效应；

S2：使用域变换递归滤波器滤波来起到边缘保持和减少噪声的作用；

S3：通过光谱解混技术得到丰度图像；

S4：对丰度图像进行上采样以获得各类别亚像元软类值；

S5：利用类别分配策略为各亚像元分配标签，得到亚像元定位结果。

对各步骤更具体的描述如下：

在步骤S1中，受到点扩展函数效应影响的高光谱成像模型表示为：

B＝I*K+N

其中，B∈R^m×n×b为高光谱遥感传感器获取得到的模糊图像，K∈R^k×k为点扩展函数所对应的模糊核，I∈R^m×n×b为所求清晰图像，N∈R^m×n×b表示加性噪声，*表示卷积运算，其中m×n表示空间分辨率，b表示波段数，k为点扩展函数对应的卷积核阶数。削弱点扩展函数导致对高光谱图像的影响的过程就是利用模糊图像和模糊核来得到清晰图像。

设B＝(B₁,B₂,...,B_b)^T，I＝(I₁,I₂,...,I_b)^T，i＝1,2,...,b，(·)^T表示转置。假设B_i满足泊松噪声分布，则有：

B_i＝Poisson((I_i*K)(x))

其中，该式表示成像系统获取到的高光谱图像与待求清晰图像的各波段灰度图用卷积核K卷积后的结果满足泊松分布。则可将待求清晰图像I的似然概率函数表示为：

其中，！表示阶乘运算，K的表示如下：

其中，σ表示高斯函数标准差，V(i,j)表示以(i,j)为中心的局部窗口的空间范围，对I的似然函数取对数得到其能量函数，结果如下；

E(I_i)＝∑{(I_i*K)-B_iln(I_i*K)}

即求其最大似然解

为：

其中，

表示目标函数在取得最小值时I_i的取值。令点扩展函数对应归一化的模糊核K满足∫K(x)dx＝1，则对上式求导得到迭代式：

其中，K^*为K的伴随矩阵，t为迭代次数，一般选择I_i＝B_i作为初始条件进行迭代，通过设置合适的迭代次数进而得到所求图像近似值I。因为随着迭代次数的增加，振铃效应也就越严重，故此方法使用一次迭代获取到所求图像近似值I。

在步骤S2中，非盲去卷积得到的清晰图像I使用域变换递归滤波器进行滤波，具体实现为：

将S1中获得的清晰高光谱图像按空间维度展开，即I＝(I₁,I₂,...,I_b)^T，i＝1,2,...,b，将任意波段I_i看成二维灰度图，将其分别沿着图像的水平和垂直方向转换为多组一维信号。对于任意给定一维信号S，其域变换可被定义为：

其中，S₀表示一维信号初始时刻的值，δ_s和δ_r为控制滤波器平滑度的两个常数，S_j-1和S_j表示第j-1时刻和第j时刻一维信号的值，·表示绝对值，U_i表示域变换后的信号强度。则递归保边滤波器可以在变换域中定义为：

J[n]＝(1-α^d)S[n]+α^dJ[n-1]

其中，

为反馈系数，d＝U_n-U_n-1表示变换域中相邻两信号间的距离，S_n表示第n时刻一维信号的强度值，J[n-1]和J[n]分别表示第n-1和第n时刻一维信号域变换递归滤波的值。对于二维图像按行展开成一维信号执行域变换递归滤波后重新组合成二维信号，再将上一步得到的二维信号按列展开成一维信号执行域变换递归滤波后重新组合成二维信号，按照该步骤递归三次得到滤波后的结果。

在步骤S3中，光谱解混技术具体为：利用基于线性解混模型的光谱解混方法来处理域变换递归滤波器滤波后的图像，以获得各类别丰度图像。线性解混模型为：线性解混模型是指混合像元内代表地物类别的端元与其所占的比例成线性关系，设步骤S2处理后的结果

d为波段数，n为像元数，

为端元矩阵，p为端元个数，

为丰度系数矩阵，

为噪声矩阵，则有，

其中，

表示1行p列的全1向量，

表示1行n列的全1向量，两个约束条件分别代表丰度系数矩阵的非负约束和每个像元丰度系数和为1约束。

在步骤S4中，上采样方法使用的径向基函数插值，通过径向基函数插值对步骤S3获得的各地物类别丰度图像进行处理，得到每个亚像元的各类别的软类值。

在步骤S5中,类别分配策略选用的是以类为单元的类别分配策略(UOC)，具体为：假设场景中含有C类地物，P_j表示第j个像元，F_k(P_j)表示第j个像元中第k类地物的丰度系数，则有，

s.t.:F_k(P_j)≥0

其中，该式表示S3得到的每个像元的丰度系数和为一，约束条件表示该像元各地物对应丰度系数非负。则该像元通过尺度因子s进行上采样使得该像元划分为s²个亚像元，则该像元内第k类地物所占亚像元个数N_k(P_j)可通过下式获取：

N_k(P_j)＝round[F_k(P_j)·s²]

其中，round[·]表示四舍五入运算。即在类别分配时，必须满足每个亚像元仅仅被分配为一类地物，并且每个像元内各类地物对应亚像元的数量与丰度系数成比例。

最后借助步骤S4获取到的各亚像元软类值，并依据衡量空间相关性的全局莫兰指数从高到低来确定地物类别分配次序，然后结合各像元内每类亚像元数量逐像元地按照类别分配次序和软类值中从大到小来为各亚像元赋予类别，分配结束后即可得到最终的亚像元定位结果。

以下通过具体实验数据来解释本发明方法的优越性：

选择JasperRidge高光谱数据集作为验证数据集进行实验。JasperRidge高光谱数据集的拥有512×614个像元，包含224个波段，光谱分辨率可达9.46nm，去除噪声波段后还留有198波段，由于其地面真实情况过于复杂，对于尺度因子为2和4时截取的100×100像元，尺度因子为3时截取99×99个像元。首先将数据集使用下采样矩阵降采样得到的低分辨率图像并以此作为亚像元定位的原始低分辨率高光谱图像，尺度因子为2时，低分辨率图像为50×50；尺度因子为3时，低分辨率图像为33×33；尺度因子为4时，低分辨率图像为25×25。其中，以尺度因子为2为例，高分辨率图像如图2(b)，其对应的低分辨率图像如图2(b)所示；对高分辨率图像采用支持向量机分类算法得到参考分类图，如图3(a)所示，类别为4，分别包括树、水、土壤和道路；本实施例所提方法进行亚像元定位得到的结果图如图3(b)所示。为反映最终的定位效果，本发明使用生产者精度(PCC)、总体精度(OA)、平均精度(AA)和Kappa系数(KC)四个指标进行评估。本实施例在尺度因子为2、3和4时生产者精度、总体精度、平均精度以及Kappa系数对比见表1，表中性能更优者采用加粗表示。其中，RBF-EP对应本实施例所提方法。

表1实施例在各尺度下性能对比

从表中可以看出，本实施例中所提的方法在OA、AA以及KC上均优于未加所提的边缘保持亚像元定位方法，并且就PCC而言，除了土壤外，其他三类也有提升，这证明了本实施例所提方法的有效性。

以上结合附图对本发明的实施方式作了详细说明，但本发明不限于所描述的实施方式。对于本领域的技术人员而言，在不脱离本发明原理和精神的情况下，对这些实施方式进行多种变化、修改、替换和变型，仍落入本发明的保护范围内。

Claims (6)

1.结合边缘保持的高光谱图像亚像元定位方法，其特征在于：包括以下步骤：

S1：利用高斯模糊核和原始高光谱图像进行非盲去卷积来削弱原始高光谱图像中点扩展函数效应；

S2：使用域变换递归滤波器滤波来起到边缘保持和减少噪声的作用；

S3：通过光谱解混技术得到丰度图像；

S4：对丰度图像进行上采样以获得各类别亚像元软类值；

S5：利用类别分配策略为各亚像元分配标签，得到亚像元定位结果。

2.根据权利要求1所述的结合边缘保持的高光谱图像亚像元定位方法，其特征在于：在步骤S1中，受到点扩展函数效应影响的高光谱成像模型表示为：

B＝I*K+N

其中，B∈R^m×n×b为高光谱遥感传感器获取得到的模糊图像，K∈R^k×k为点扩展函数所对应的模糊核，I∈R^m×n×b为所求清晰图像，N∈R^m×n×b表示加性噪声，*表示卷积运算，其中m×n表示空间分辨率，b表示波段数，k为点扩展函数对应的卷积核阶数；

设B＝(B₁,B₂,...,B_b)^T，I＝(I₁,I₂,...,I_b)^T，i＝1,2,...,b，假设B_i满足泊松噪声分布，将待求清晰图像I的似然概率函数表示为：

其中，！表示阶乘运算，K为高斯模糊核，最大似然解出现在似然函数偏导数为零的情况下，则得到迭代式：

其中，K^H为K的共轭转置矩阵，n为迭代次数，选择I_i＝B_i作为初始条件进行迭代，通过设置合适的迭代次数进而得到所求图像近似值I。

3.根据权利要求1所述的结合边缘保持的高光谱图像亚像元定位方法，其特征在于：在步骤S2中，非盲去卷积得到的清晰图像I使用域变换递归滤波器进行滤波，具体实现为：令I＝(I₁,I₂,...,I_b)^T，i＝1,2,...,b，将任意波段I_i看成二维灰度图，将其分别沿着图像的水平和垂直方向转换为多组一维信号，对于任意给定一维信号S，其域变换可被定义为：

其中，S₀表示一维信号初始时刻的值，δ_s和δ_r为控制滤波器平滑度的两个常数，S_j-1和S_j表示第j-1时刻和第j时刻一维信号的值，|·|表示绝对值，U_i表示域变换后的信号强度，递归保边滤波器可以在变换域中定义为：

J[n]＝(1-α^d)S[n]+α^dJ[n-1]

其中，α为反馈系数，d＝U_n-U_n-1表示变换域中相邻两信号间的距离，S_n表示第n时刻一维信号的强度值，J[n-1]和J[n]分别表示第n-1和第n时刻一维信号域变换递归滤波的值。

4.根据权利要求1所述的结合边缘保持的高光谱图像亚像元定位方法，其特征在于：在步骤S3中，光谱解混技术具体为：利用基于线性解混模型的光谱解混方法来处理域变换递归滤波器滤波后的图像，以获得各类别丰度图像，线性解混模型为：线性解混模型是指混合像元内代表地物类别的端元与其所占的比例成线性关系，则设Y为待解混的高光谱数据，A为端元矩阵，X为丰度矩阵，N₁为噪声矩阵，则有

Y＝AX+N₁

其中，丰度矩阵X满足非负约束以及和为1约束。

5.根据权利要求1所述的结合边缘保持的高光谱图像亚像元定位方法，其特征在于：在步骤S4中，通过上采样方法对步骤S3获得的各地物类别丰度图像进行处理，得到各亚像元的每种类别的软类值。

6.根据权利要求1所述的结合边缘保持的高光谱图像亚像元定位方法，其特征在于：在步骤S5中，类别分配策略选用的是以类为单元的类别分配策略，具体为：

假设场景中含有C类地物，P_j表示第j个像元，F_k(P_j)表示第j个像元中第k类地物的丰度系数，则该像元通过尺度因子s进行上采样使得该像元划分为s²个亚像元，则该像元内第k类地物所占亚像元个数N_k(P_j)可通过下式获取：

N_k(P_l)＝round[F_k(P_j)·s²]

其中，round[·]表示四舍五入运算。

Images