CN114820741A - 一种高光谱影像全波段超分重建方法 - Google Patents

Abstract

本发明公开了一种高光谱影像全波段超分重建方法，涉及遥感图像处理领域；通过分析获取的HS遥感图像和MS遥感图像的光谱范围，划分波段；定义光谱范围相同的波段为匹配波段，不同的波段为非匹配波段；对于波段匹配部分，采用小波变换方法进行融合处理，融合后的重建图像的纹理细节等方面更清晰；对于非匹配部分，利用匹配波段的重建结果与原始HS遥感图像之间像元值的前后变化为训练知识，构建BP神经网络。以非匹配波段HS遥感图像的像元作为网络的输入，结合HS遥感图像全波段归一化的信息熵以及像元的邻域空间关系对模拟结果进行优化，实现HS遥感图像全波段超分重建，重建后的HS遥感图像具有MS遥感图像的高空间分辨率，还保留了本身高光谱分辨率。

Description

一种高光谱影像全波段超分重建方法

技术领域

本发明涉及遥感图像处理领域，尤其涉及一种高光谱影像全波段超分重建方法。

背景技术

高光谱(Hyper-Spectral,HS)图像收集了数百个相邻波段，提供了不同地物的精细光谱细节。 HS遥感图像以其高光谱分辨率在地物识别与地面目标分类、环境监测、土地利用与土地覆盖 变化检测等多个领域具有良好的应用前景。而HS影像空间分辨率低，严重限制了其空间精 细解译能力。通常，提高其空间分辨率可以从遥感设备和图像处理两个方面来实现。与后者 相比，前者在技术上总是困难得多，成本也昂贵得多。

国内外研究者在HS遥感图像超分重建领域做出了许多工作，有力推动了遥感图像重建方 法的发展。HS遥感图像重建起源于全色(Panchromatic,PAN)图像锐化。PAN图像锐化指的是 多光谱(Multi-Spectral,MS)图像与相应的高分辨率PAN图像融合来增强MS遥感图像的空间 分辨率。基于图像融合的HS遥感图像超分重建不同于PAN图像泛锐化。HS遥感图像比PAN 图像具有更丰富的光谱信息，在HS遥感图像泛锐化的过程中，会出现光谱失真等问题。因 此，一些学者对典型的融合算法进行了改进，并将其应用于HS遥感图像和MS遥感图像融 合。目前，典型的HS遥感图像和MS遥感图像空谱融合方法主要分为四类：基于矩阵变换 的方法、基于混合光谱解混理论的方法、基于概率统计理论的方法和基于深度学习的方法。

(1)基于矩阵变换的方法

该方法大体上包括成分替换以及频域滤波。成分替换法中最经典的方法是 GS(Gram-Schmidt)变换，HS遥感图像和MS遥感图像分别进行GS变换得到各自的GS分量， 用MS遥感图像的第一GS分量替换HS遥感图像的第一GS分量，新组合的HS遥感图像GS 分量进行GS逆变换得到高空间分辨率HS遥感图像。频域滤波法采用滤波提取HS遥感图像 的低频信息以及MS遥感图像的高频信息，对提取出的信息通过反向滤波等操作获得融合图 像。Nardecchia等人对HS遥感图像的每个波段进行平稳小波变换。平稳小波变换利用固定的滤波器将HS遥感图像在三个方向上分解为低频小波系数和高频小波系数。用MS遥感图像代替低频小波系数，参与小波重构，得到融合图像。基于频域滤波的超分辨率重建过于理想化，无法应对噪声等因素造成的图像分辨率降低。

(2)基于混合光谱解混理论的方法

该方法利用高空间分辨率数据的空间信息，辅助对HS遥感图像进行光谱解混，得到高空 间尺度下的像元光谱。Zhang等人对HS遥感图像和MS遥感图像分别进行混合像元解混得到 对应的丰度矩阵和端元矩阵，对端元矩阵进行最小化约束以及对丰度矩阵进行稀疏性约束， 将约束后的HS遥感图像丰度矩阵与MS遥感图像的端元矩阵混合得到高空间分辨率HS遥感 图像。然而，该类方法存在处理时间长，且实时性较差的特点。

(3)基于概率统计理论的方法

这类方法主要包含两个方面：能量函数的构建和优化求解。将HS遥感图像超分重建视作 一个线性回归问题，根据理想的高空间分辨率HS遥感图像和表示退化数据的HS遥感图像之 间的观测模型设计一个能量函数，然后对能量函数进行优化求解。优化求解通常采用迭代的 方式，如梯度下降算法、共轭梯度算法等。Song等人根据光谱的低维特性，将MS遥感图像 和HS遥感图像划分为几个部分，分别对每个部分的理想图像进行近似估计。对于大多数优 化解，都是通过迭代求解得到的融合图像。这类方法虽然效果较好，但其计算效率过低，严 重影响了这类方法的应用。

(4)基于深度学习的方法

复杂多样的神经网络模型能够较好地提取HS遥感图像的光谱特征和MS遥感图像的空间 特征，实现HS遥感图像的空谱融合。Palsson等人提出使用卷积神经网络(Convolutional Neural Networks,CNN)实现超分重建。该算法通过PCA提取HS遥感图像的几个主成分和MS遥感 图像将作为CNN的输入，利用MS遥感图像与HS遥感图像的几个主成分之间的分辨率因子 提取MS遥感图像的空间特征，并采用双三次插值法，得到融合图像。融合后的图像通过逆 变换生成高空间分辨率HS遥感图像。Han等人提出了一种新的基于多分支BP(Back Propagation)神经网络的聚类图像融合算法。该算法的核心是找到HS遥感图像和MS遥感图 像之间的光谱映射。通过将对应的MS遥感图像和HS遥感图像输入到BP神经网络中，得到 地物的非线性光谱映射。根据该光谱映射，将MS遥感图像作为该BP神经网络的输入得到 高空间分辨率HS遥感图像。尽管基于深度学习的融合方法因其优异的空间信息获取能力而 得到广泛的认可，但是网络训练需要丰富的数据集，而HS遥感图像样本获取困难，公开的 训练数据集少。并且，想要获取与HS遥感图像匹配的MS遥感图像难度高，限制了深度学 习在HS融合领域的发展。

近些年，一些学者考虑到HS遥感图像的时间分辨率和覆盖范围不足，采用宽覆盖和高时 间分辨率数据融合弥补了HS遥感图像的缺陷，时-空-谱四维数据一体化融合及光谱信息的时 空范围拓展逐渐得到重视。Sun等人认为HS遥感图像中像元与MS遥感图像的对应像元存在 线性关系，该线性关系可通过最小二乘估计表达为一组变换矩阵。然后，根据该变换矩阵， 将MS遥感图像逐像元地转化为模拟图像。模拟图像的光谱分辨率与HS遥感图像相同，并 且空间分辨率也得到了提高。Zhao等人将HS遥感图像的光谱划分为多个区域，并采用小波 变换在每个区域内融合HS遥感图像和MS遥感图像，采用比值重采样的方式对HS遥感图像 中非多光谱波段覆盖的缺失数据进行插值。这些方法不仅使得HS遥感图像空间分辨率得到 了提升，时间分辨率和覆盖范围也与MS遥感图像一致。

综上所述，大部分基于融合的超分重建都是针对多光谱波段范围之内的HS遥感图像，且 算法种类繁多、成熟。而对于多光谱波段范围之外的HS遥感图像的重建或者模拟的研究工 作相对缺乏。故本发明采用融合和模拟策略对HS遥感图像进行超分重建研究与实现，该方 法能有效地对HS遥感图像的空间分辨率进行改进，且HS遥感图像的光谱信息保存性也相对 较好。

发明内容

针对现有技术的不足，本发明提出一种高光谱影像全波段超分重建方法；采用小波变换 和BP神经网络，实现了高光谱影像的全波段超分重建。

所采用的技术方案为：

一种高光谱影像全波段超分重建方法，包括以下步骤：

步骤1：获取HS遥感图像和MS遥感图像；

步骤2：HS遥感图像和MS遥感图像匹配波段的小波融合；

对HS遥感图像进行上采样，使HS遥感图像与MS遥感图像具有相同的空间域和像元尺 寸；对MS遥感图像进行小波分解,得到匹配波段的HS遥感图像的小波分解；对于所有的分 解层，采用低频融合规则融合第p波段的HS遥感图像H_p ^(RS_m)和第q波段的MS遥感图像M_q的低频系数，高频融合规则融合第p波段的HS遥感图像H_p ^(RS_m)和第q波段的MS遥感图像 M_q的三个方向的高频系数；从尺度最大的层次逐层重构到尺度最小的层次进行小波重构，直 到最后一层，最后得到融合图像H^(RC_m)；

步骤3：BP神经网络重建不匹配波段；

根据遥感图像中重建后的高空间分辨率HS遥感图像H^(RC)和重采样后的HS遥感图像H^(RS)的像元，匹配波段的融合图像和HS遥感图像对应像元集合{H^(RS_m),H^(RC_m)}作为训练样本， 其中，H^(RS_m)为H^(RS)的匹配波段图像，H^(RC_m)为H^(RC)的匹配波段图像；训练BP神经网络重 建非匹配波段图像H^(RC_u)；选取BP神经网络参数；网络结构设置为三层；通过像元邻域空间 关系结合波段信息熵的形式优化H^(RC_u)'的所有波段；信息熵归一化后表示为n_IE^(RS)；结合 H^(RC_m)和H^(RC_u)得到超分辨率重建图像H^(RC)。

所述步骤2的具体过程，包括以下步骤：

步骤2.1：HS遥感图像，如式(1)所示：

H＝{H_i(x'_i,y'_i),i＝1,2,…,I} (1)

其中，i为像元索引；I为像元总数；

为像元i的网格坐标，其中， D_i ^(H)为像元i的分辨率单位面积，D^(H)＝D₁ ^(H)∪…∪D_i ^(H)…∪D_I ^(H)是H的空间域；H_i＝(H_ip；p ∈Ω)为像元i的光谱测度矢量，p为HS遥感图像的波段索引，Ω＝{1,2,…,P}，P为波段总 数，H_ip为像元i在p波段的光谱值；

MS遥感图像为：

M＝{M_j(x'_j,y'_j),j＝1,2,…,J} (2)

其中，j为像元索引，J为像元总数，

是像元j的网格坐标，D_j ^(M)为 像元j的分辨率单位面积，D^(M)＝D₁ ^(M)∪…∪D_j ^(M)…∪D_J ^(M)是M的空间域；M_j＝(M_jq；q∈Λ)为像元j的光谱测度矢量，q为MS遥感图像的波段索引，Λ＝{1,2,…,Q}，Q是波段总数， M_jq为像元j在q波段的光谱值；

fr_p ^(H)为HS遥感图像第p波段的中心频率；B_p ^(H)为HS遥感图像第p波段的光谱范围；fr_q ^(M)为MS遥感图像第q波段的中心频率；B_q ^(M)为MS遥感图像第q波段的光谱范围；HS遥感图像中光谱范围B_p ^(H)位于MS遥感图像光谱范围B_q ^(M)之内的波段为匹配波段，如式(3)所示：

HS遥感图像的光谱范围与MS遥感图像的光谱范围不重叠的波段为非匹配波段，如式(4) 所示：

Ω_u＝Ω\Ω_m (4)

重建后的高空间分辨率HS遥感图像为：

H^(RC)＝{H_j ^(RC)(x'_i,y'_i)，j＝1,2,…,J} (5)

其中，H_j ^(RC)＝(H_jp ^(RC)；p∈Ω)是像元j的光谱测度矢量，H_jk ^(RC)为重建后的光谱值；

所述H^(RC)具有与M相同的空间结构与H相同的光谱结构；H的超分重建通过M估计H在匹配和非匹配波段的H^(RC)的光谱值；对于匹配波段

H的超分重建图像表示为：

H^(RC_m)＝{H_j ^(RC_m)(x'_i,y'_i),j＝1,2,…,J} (6)

H_j ^(RC_m)＝(H_jp ^(RC)；p∈Ω_m) (7)

对于非匹配波段Ω_u＝Ω\Ω_m，H的超分重建图像表示为：

H^(RC_u)＝{H_j ^(RC_u)(x'_i,y'_i),j＝1,2,…,J} (8)

H_j ^(RC_u)＝(H_jp ^(RC)；p∈Ω_u) (9)

步骤2.2：对HS遥感图像H＝{H_i(x'_i,y'_i)；i＝1,2,…,I}进行上采样，如图2所示，使HS 遥感图像与MS遥感图像具有相同的空间域和像元尺寸；

重采样后的HS遥感图像为：

H^(RS)＝{H_j ^(RS)(x'_j,y'_j)；j＝1,2,…,J} (10)

H_j ^(RS)＝(H_jp ^(RS)；p∈Ω) (11)

根据Ω_u和Ω_m，H^(RS)的匹配波段图像和非匹配波段图像表示为：

H^(RS_m)＝{H_j ^(RS_m)(x'_j,y'_j)；j＝1,2,…,J} (12)

H^(RS_u)＝{H_j ^(RS_u)(x'_j,y'_j)；j＝1,2,…,J} (13)

H_j ^(RS_m)＝(H_jp ^(RS)；p∈Ω_m)andH_j ^(RS_u)＝(H_jp ^(RS)；p∈Ω_u) (14)

步骤2.3：MS遥感图像是索引为(x'_j,y'_j)的二维信号，其中，M_j(x'_j,y'_j)表示MS遥感图像 第x'_j行y'_j列的像元值；M_j(x'_j,:)为x'_j行的像元集合，M_j(:,y'_j)为y'_j列的像元集合；如图3所 示，MS遥感图像的小波分解为：

其中，ξ_gg为小波分解变换，表示采用低通滤波器g先垂直方向后水平方向的下采样；ξ_hg为小波分解变换，表示先采用低通滤波器g垂直方向下采样，后采用高通滤波器h水平方向 下采样；ξ_gh为小波分解变换，表示先采用高通滤波器h垂直方向下采样，后采用低通滤波器 g水平方向下采样；ξ_hh为小波分解变换，表示采用高通滤波器h先垂直方向后水平方向的下 采样；M_q ^(m)表示第q波段的MS遥感图像M_q经过尺度m次小波分解后的低频系数；Mh_q ^(m+1)表示M_q ^(m)小波分解后水平方向的高频系数；Mv_q ^(m+1)表示M_q ^(m)小波分解后垂直方向的高频系数；Md_q ^(m+1)表示M_q ^(m)小波分解后对角方向的高频系数；

匹配波段的HS遥感图像的小波分解为，

其中，H_p ^(RS_m)(m)表示第p波段的HS遥感图像H_p ^(RS_m)经过尺度m次小波分解后的低频系 数；Hh_p ^(RS_m)(m+1)表示H_p ^(RS_m)(m)小波分解后水平方向的高频系数；Hv_p ^(RS_m)(m+1)表示H_p ^{(RS _m)(m)}小波分解后垂直方向的高频系数；Hd_p ^(RS_m)(m+1)表示H_p ^(RS_m)(m)小波分解后对角方向的高频系数；

步骤2.4：对于所有的分解层，采用低频融合规则融合H_p ^(RS_m)和M_q的低频系数，采用高 频融合规则融合H_p ^(RS_m)和M_q的三个方向的高频系数；尺度m下的单波段HS-MS遥感图像的低频信息融合表示为：

其中，a'是权重系数，0<a'<1，WF_p ^(RC)(m)表示低频小波系数WF_p ^(RC)尺度m下的小波低频系数；

尺度m下的单波段HS-MS遥感图像的高频信息绝对值最大的融合表示为：

其中，WFh_p ^(RC)(m)表示WF_p ^(RC)尺度m下水平方向的小波高频系数；WFv_p ^(RC)(m)表示WF_p ^(RC)尺度m下垂直方向的小波高频系数；WFd_p ^(RC)(m)表示WF_p ^(RC)尺度m下对角方向的小波高频系数；

步骤2.5：小波重构是小波分解的逆过程，从尺度大的层次逐层重构到尺度小的层次； MS遥感图像的小波重构可表示为：

其中，ζ_gg为小波重构变换，表示采用低通滤波器g先水平方向后垂直方向的上采样；ζ_hg为小波重构变换，表示先采用高通滤波器h水平方向上采样，后采用低通滤波器g垂直方向 上采样；ζ_gh为小波重构变换，表示先采用低通滤波器g水平方向上采样，后采用高通滤波器 h垂直方向上采样；ζ_hh为小波重构变换，表示采用高通滤波器h先水平方向后垂直方向的上 采样；M_q ^(m-1)是m层小波系数进行小波逆变换得到的低频小波系；M_q ^(m-1)作为m层小波逆变换 的低频小波系数，参与(m-1)层小波逆变换，当m＝1时，得到重构小波图像M_q ⁽⁰⁾；

HS遥感图像的小波重构可表示为：

从尺度m的分解层开始进行小波重构，m层的一个低频系数和三个方向的高频系数参与 重构，重构得到的图像作为m-1层低频小波系数参与m-1层的重构，依次递减，直到最后一 层，最后得到融合图像H^(RC_^m)。

所述步骤3的具体过程，包括以下步骤：

步骤3.1：根据遥感图像中H^(RC)和H^(RS)的像元，像元j及其邻域像元在p波段上的光谱 测量值{H_jp ^(RC)}和{H_jp ^(RS)}有相同的空间关系，匹配波段的融合图像和HS遥感图像对应像元 集合{H^(RS_m),H^(RC_m)}作为训练样本训练BP神经网络重建H^(RC_u)；

网络正向传播时的网络输入为：

Int_ip＝(Int_ipι,ι＝1,...,Θ)＝(H_jp ^(RS_m)(x'_j,y'_j) (21)

(x'_j,y'_j)∈D_i ^(H)；p∈Ω_m) (22)

网络正向传播时的网络输出为：

Out_ip＝(Out_ipυ,υ＝1,...,Θ)＝(H_jp ^(RC_m)(x'_j,y'_j) (23)

(x'_j,y'_j)∈D_i ^(H)；p∈Ω_m) (24)

其中，Θ为图像空间域D_i ^(H)个数和波段p的乘积，BP神经网络的样本集{H^(RS_m),H^{(RC _m)}} 表示为{Int_ip,Out_ip}；

步骤3.2：选取BP神经网络参数，包括以下步骤：

步骤3.2.1：选取网络的层数，以Int_ip作为输入，Out_ip作为输出，选取三层BP神经网络， 包括输入层、单层隐含层和输出层；

步骤3.2.2：选取隐含层神经元个数，BP神经网络的输入层和输出层的神经元个数为Θ， 隐含层神经元个数为：

其中，Θ'为隐含层神经元个数；c为取值1～10的常数；隐含层神经元个数与输入参数和 输出参数的设置有关；

隐含层输出，如式(26)所示：

步骤3.2.3：设定训练次数，训练次数设置为100次；

步骤3.2.4：设定学习效率，学习效率是在数据反向传递的过程中起到调整权值的作用， 取值为0.01～0.1之间；

步骤3.2.5：设定期望最小误差，设定网络期望最小误差为0.00001，当训练结果满足该 精度，停止训练；

通过设定BP神经网络的初始值后，以Int_ip作为网络输入，以Out_ip作为网络期望输出， 随着信号的正向传播，误差反向传播，网络权值的不断调整，直至网络的最终输出与Out_ip的相对误差小于设定的阈值N₁；

步骤3.3：网络结构设置为三层，重建后的H^(RC_u)'的像元邻域空间关系表示为：

H^(RC_u)'＝{H_j ^(RC_u)'(x'_i,y'_i)；j＝1,2,…,J} (27)

H_j ^(RC_u)'＝(H_jp ^(RC)'；p∈Ω_u) (28)

利用H^(RS)的信息熵和H_jp ^(RC_m)的像元邻域关系对重建图像H^(RC_^u)'的像元拟合现象进行优 化；

计算每个波段的重建图像H^(RC_m)在空间域D_i ^(H)内的邻域像元集合H_jp ^(RC_m)与其均值的距 离为：

其中，mean_ip表示空间域D_i ^(H)内邻域像元的均值；dist_jp是空间域D_i ^(H)内像元j与均值mean_ip的差，表示像元j在空间域D_i ^(H)内的空间关系；计算dist_j'p在匹配波段p∈Ω_m中波段均值 mean_dist_j'，如式(31)所示：

其中，P'表示匹配波段总数；

信息熵表征图像单波段的整体信息，通过像元邻域空间关系结合波段信息熵的形式优化 H^(RC_u)'的所有波段为：

其中，H_jp ^(RC),p∈Ω_u表示优化后的光谱值，IE^(RS)为H^(RS)的每个波段的信息熵，IE^(RS)＝ (IE_p ^(RS)；p∈Ω)，信息熵归一化后表示为：

n_IE^(RS)＝(n_IE_p ^(RS)；p∈Ω) (33)

结合H^(RC_m)和H^(RC_u)得到超分辨率重建图像H^(RC)。

本发明有益技术效果：

1、本发明提出了一种结合小波变换和BP神经网络的高光谱影像全波段超分重建方法。 通过分析获取的HS遥感图像和MS遥感图像的光谱范围，对HS遥感图像和MS遥感图像的 波段进行划分。定义光谱范围相同的波段为匹配波段，定义光谱范围不同的波段为非匹配波 段。对于HS遥感图像和MS遥感图像波段匹配部分，采用小波变换方法进行融合处理，融 合后的重建图像的纹理细节等方面更清晰。

2、对于HS遥感图像和MS遥感图像波段非匹配部分，利用匹配波段的重建结果与原始 HS遥感图像之间像元值的前后变化为训练知识，构建BP神经网络。以非匹配波段HS遥感图像的像元作为该网络的输入，结合HS遥感图像全波段归一化的信息熵以及像元的邻域空间关系对模拟结果进行优化，实现HS遥感图像全波段超分重建，重建后的HS遥感图像不仅具有MS遥感图像的高空间分辨率，还保留了本身高光谱分辨率的特点，并结合多种指标对非匹配波段重建图像进行分析，发现整体上具有较好的重建效果。

附图说明

图1为本发明实施例提供的一种高光谱影像全波段超分重建方法流程图；

图2为本发明实施例提供的HS-MS遥感图像波段匹配示意图；

图3为本发明实施例提供的小波融合的步骤流程图；

图4为本发明实施例提供的BP神经网络重建过程示意图；

图5为本发明实施例提供的Sentinel-2B各波段对应的图像示意图；

其中，图5(a)为Sentinel-2B 1波段对应的图像示意图，图5(b)为Sentinel-2B 2波段对应 的图像示意图，图5(c)为Sentinel-2B 3波段对应的图像示意图，图5(d)为Sentinel-2B 4波段 对应的图像示意图，图5(e)为GaoFen-530波段对应的图像示意图，图5(f)为GaoFen-540波 段对应的图像示意图，图5(g)为GaoFen-565波段对应的图像示意图，图5(h)为GaoFen-5100 波段对应的图像示意图；

图6为本发明实施例提供的Sentinel-2B各波段对应的图像示意图；

其中，图6(a)为Sentinel-2B 1波段和GaoFen-530波段的小波融合图像示意图，图6(b) 为Sentinel-2B 2波段和GaoFen-540波段的小波融合图像示意图，图6(c)为Sentinel-2B 3波 段和GaoFen-565波段的小波融合图像示意图，图6(d)为Sentinel-2B 4波段和GaoFen-5100 波段的小波融合图像示意图；

图7为本发明实施例提供的GaoFen-5各波段对应的图像示意图；

其中，图7(a)为GaoFen-515波段对应的图像示意图，图7(b)为GaoFen-535波段对应的 图像示意图，图7(c)为GaoFen-555波段对应的图像示意图，图7(d)为GaoFen-580波段对应 的图像示意图，图7(e)为GaoFen-5160波段对应的图像示意图，图7(f)为基于小波变换的BP 神经网络重建图像H^(RC)的15波段对应的图像示意图，图7(g)为基于小波变换的BP神经网络 重建图像H^(RC)的35波段对应的图像示意图，图7(h)为基于小波变换的BP神经网络重建图像 H^(RC)的55波段对应的图像示意图，图7(i)为基于小波变换的BP神经网络重建图像H^(RC)的80 波段对应的图像示意图，图7(j)为基于小波变换的BP神经网络重建图像H^(RC)的160波段对 应的图像示意图；

图8为本发明实施例提供的像元大小为300×300的GaoFen-5的重建三维立体图；

图9为本发明实施例提供的H和H^(RC)间相关系数和结构相似性示意图；

其中，图9(a)为H和H^(RC)之间的CC示意图，图9(b)为H和H^(RC)之间的SSIM示意图；

图10为本发明实施例提供的H和H^(RC)的信息熵示意图；

图11为本发明实施例提供的附图11为H和H^(RC)的光谱梯度变化示意图；

其中，图11(a)为H和H^(RC)水体的光谱一阶微分示意图，图11(b)为H和H^(RC)建筑的光谱一阶微分示意图，图11(c)为H和H^(RC)植被的光谱一阶微分示意图，图11(d)为H和H^(RC)土壤的光谱一阶微分示意图。

具体实施方式

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

本实施例提出了一种高光谱影像全波段超分重建方法，采用小波变换和BP神经网络， 实现了高光谱影像的全波段超分重建；如图1所示，包括以下步骤：

一种高光谱影像全波段超分重建方法，包括以下步骤：

步骤1：获取HS遥感图像和MS遥感图像；

步骤2：HS遥感图像和MS遥感图像匹配波段的小波融合；

对HS遥感图像进行上采样，使HS遥感图像与MS遥感图像具有相同的空间域和像元尺 寸；对MS遥感图像进行小波分解,得到匹配波段的HS遥感图像的小波分解；对于所有的分 解层，采用低频融合规则融合第p波段的HS遥感图像H_p ^(RS_m)和第q波段的MS遥感图像M_q的低频系数，高频融合规则融合第p波段的HS遥感图像H_p ^(RS_m)和第q波段的MS遥感图像 M_q的三个方向的高频系数；从尺度最大的层次逐层重构到尺度最小的层次进行小波重构，直 到最后一层，最后得到融合图像H^(RC_m)，包括以下步骤：

步骤2.1：HS遥感图像，如式(1)所示：

H＝{H_i(x'_i,y'_i),i＝1,2,…,I} (1)

其中，i为像元索引；I为像元总数；

为像元i的网格坐标，其中， D_i ^(H)为像元i的分辨率单位面积，D^(H)＝D₁ ^(H)∪…∪D_i ^(H)…∪D_I ^(H)是H的空间域；H_i＝(H_ip；p ∈Ω)为像元i的光谱测度矢量，p为HS遥感图像的波段索引，Ω＝{1,2,…,P}，P为波段总 数，H_ip为像元i在p波段的光谱值；

MS遥感图像为：

M＝{M_j(x'_j,y'_j),j＝1,2,…,J} (2)

其中，j为像元索引，J为像元总数，

是像元j的网格坐标，D_j ^(M)为 像元j的分辨率单位面积，D^(M)＝D₁ ^(M)∪…∪D_j ^(M)…∪D_J ^(M)是M的空间域；M_j＝(M_jq；q∈Λ)为像元j的光谱测度矢量，q为MS遥感图像的波段索引，Λ＝{1,2,…,Q}，Q是波段总数， M_jq为像元j在q波段的光谱值；

fr_p ^(H)为HS遥感图像第p波段的中心频率；B_p ^(H)为HS遥感图像第p波段的光谱范围；fr_q ^(M)为MS遥感图像第q波段的中心频率；B_q ^(M)为MS遥感图像第q波段的光谱范围；HS遥感图像中光谱范围B_p ^(H)位于MS遥感图像光谱范围B_q ^(M)之内的波段为匹配波段，如式(3)所示：

HS光谱范围与MS遥感图像不重叠的波段为非匹配波段，如式(4)所示：

Ω_u＝Ω\Ω_m (4)

重建后的高空间分辨率HS遥感图像为：

H^(RC)＝{H_j ^(RC)(x'_i,y'_i)，j＝1,2,…,J} (5)

其中，H_j ^(RC)＝(H_jp ^(RC)；p∈Ω)是像元j的光谱测度矢量，H_jk ^(RC)为重建后的光谱值；

所述H^(RC)具有与M相同的空间结构与H相同的光谱结构；H的超分重建通过M估计H在匹配和非匹配波段的H^(RC)的光谱值；对于匹配波段

，H的超分重建图像表示为：

H^(RC_m)＝{H_j ^(RC_m)(x'_i,y'_i),j＝1,2,…,J} (6)

H_j ^(RC_m)＝(H_jp ^(RC)；p∈Ω_m) (7)

对于非匹配波段Ω_u＝Ω\Ω_m，H的超分重建图像表示为：

H^(RC_u)＝{H_j ^(RC_u)(x'_i,y'_i),j＝1,2,…,J} (8)

H_j ^(RC_u)＝(H_jp ^(RC)；p∈Ω_u) (9)

步骤2.2：对HS遥感图像H＝{H_i(x'_i,y'_i)；i＝1,2,…,I}进行上采样，使HS遥感图像与 MS遥感图像具有相同的空间域和像元尺寸；

重采样后的HS遥感图像为：

H^(RS)＝{H_j ^(RS)(x'_j,y'_j)；j＝1,2,…,J} (10)

H_j ^(RS)＝(H_jp ^(RS)；p∈Ω) (11)

根据Ω_u和Ω_m，H^(RS)的匹配波段图像和非匹配波段图像表示为：

H^(RS_m)＝{H_j ^(RS_m)(x'_j,y'_j)；j＝1,2,…,J} (12)

H^(RS_u)＝{H_j ^(RS_u)(x'_j,y'_j)；j＝1,2,…,J} (13)

H_j ^(RS_m)＝(H_jp ^(RS)；p∈Ω_m)andH_j ^(RS_u)＝(H_jp ^(RS)；p∈Ω_u) (14)

步骤2.3：MS遥感图像是索引为(x'_j,y'_j)的二维信号，其中，M_j(x'_j,y'_j)表示MS遥感图像 第x'_j行y'_j列的像元值；M_j(x'_j,:)为x'_j行的像元集合，M_j(:,y'_j)为y'_j列的像元集合；MS遥感图 像的小波分解为，

其中，ξ_gg为小波分解变换，表示采用低通滤波器g先垂直方向后水平方向的下采样；ξ_hg为小波分解变换，表示先采用低通滤波器g垂直方向下采样，后采用高通滤波器h水平方向 下采样；ξ_gh为小波分解变换，表示先采用高通滤波器h垂直方向下采样，后采用低通滤波器 g水平方向下采样；ξ_hh为小波分解变换，表示采用高通滤波器h先垂直方向后水平方向的下 采样；M_q ^(m)表示第q波段的MS遥感图像M_q经过尺度m次小波分解后的低频系数；Mh_q ^(m+1)表示M_q ^(m)小波分解后水平方向的高频系数；Mv_q ^(m+1)表示M_q ^(m)小波分解后垂直方向的高频系数；Md_q ^(m+1)表示M_q ^(m)小波分解后对角方向的高频系数；

匹配波段的HS遥感图像的小波分解为：

其中，H_p ^(RS_m)(m)表示第p波段的HS遥感图像H_p ^(RS_m)经过尺度m次小波分解后的低频系 数；Hh_p ^(RS_m)(m+1)表示H_p ^(RS_m)(m)小波分解后水平方向的高频系数；Hv_p ^(RS_m)(m+1)表示H_p ^{(RS _m)(m)}小波分解后垂直方向的高频系数；Hd_p ^(RS_m)(m+1)表示H_p ^(RS_m)(m)小波分解后对角方向的高频系数；

步骤2.4：对于所有的分解层，采用低频融合规则融合H_p ^(RS_m)和M_q的低频系数，采用高 频融合规则融合H_p ^(RS_m)和M_q的三个方向的高频系数；尺度m下的单波段HS-MS遥感图像的低频信息融合表示为：

其中，a'是权重系数，0<a'<1，WF_p ^(RC)(m)表示低频小波系数WF_p ^(RC)尺度m下的小波低频系数；

尺度m下的单波段HS-MS遥感图像的高频信息绝对值最大的融合表示为：

其中，WFh_p ^(RC)(m)表示WF_p ^(RC)尺度m下水平方向的小波高频系数；WFv_p ^(RC)(m)表示WF_p ^(RC)尺度m下垂直方向的小波高频系数；WFd_p ^(RC)(m)表示WF_p ^(RC)尺度m下对角方向的小波高频系数；

步骤2.5：小波重构是小波分解的逆过程，从尺度大的层次逐层重构到尺度小的层次； MS遥感图像的小波重构可表示为：

其中，ζ_gg为小波重构变换，表示采用低通滤波器g先水平方向后垂直方向的上采样；ζ_hg为小波重构变换，表示先采用高通滤波器h水平方向上采样，后采用低通滤波器g垂直方向 上采样；ζ_gh为小波重构变换，表示先采用低通滤波器g水平方向上采样，后采用高通滤波器h垂直方向上采样；ζ_hh为小波重构变换，表示采用高通滤波器h先水平方向后垂直方向的上 采样；M_q ^(m-1)是m层小波系数进行小波逆变换得到的低频小波系；M_q ^(m-1)作为m层小波逆变换 的低频小波系数，参与(m-1)层小波逆变换，当m＝1时，得到重构小波图像M_q ⁽⁰⁾；

HS遥感图像的小波重构可表示为：

从尺度m的分解层开始进行小波重构，m层的一个低频系数和三个方向的高频系数参与 重构，重构得到的图像作为m-1层低频小波系数参与m-1层的重构，依次递减，直到最后一 层，最后得到融合图像H^(RC_m)；

步骤3：BP神经网络重建不匹配波段；

根据遥感图像中重建后的高空间分辨率HS遥感图像H^(RC)和重采样后的HS遥感图像H^(RS)的像元，匹配波段的融合图像和HS遥感图像对应像元集合{H^(RS_m),H^(RC_m)}作为训练样本， 其中，H^(RS_m)为H^(RS)的匹配波段图像，H^(RC_m)为H^(RC)的匹配波段图像；训练BP神经网络重 建非匹配波段图像H^(RC_u)；选取BP神经网络参数；网络结构设置为三层；通过像元邻域空间 关系结合波段信息熵的形式优化H^(RC_u)'的所有波段；信息熵归一化后表示为n_IE^(RS)；结合 H^(RC_m)和H^(RC_u)得到超分辨率重建图像H^(RC)，包括以下步骤：

步骤3.1：根据遥感图像中H^(RC)和H^(RS)的像元，像元j及其邻域像元在p波段上的光谱 测量值{H_jp ^(RC)}和{H_jp ^(RS)}有相同的空间关系，匹配波段的融合图像和HS遥感图像对应像元 集合{H^(RS_m),H^(RC_m)}作为训练样本训练BP神经网络重建H^(RC_u)；

网络正向传播时的网络输入为：

Int_ip＝(Int_ipι,ι＝1,...,Θ)＝(H_jp ^(RS_m)(x'_j,y'_j) (21)

(x'_j,y'_j)∈D_i ^(H)；p∈Ω_m) (22)

网络正向传播时的网络输出为：

Out_ip＝(Out_ipυ,υ＝1,...,Θ)＝(H_jp ^(RC_m)(x'_j,y'_j) (23)

(x'_j,y'_j)∈D_i ^(H)；p∈Ω_m) (24)

其中，Θ为图像空间域D_i ^(H)个数和波段p的乘积，BP神经网络的样本集{H^(RS_m),H^{(RC _m)}} 表示为{Int_ip,Out_ip}；

步骤3.2：选取BP神经网络参数，如图4所示，包括以下步骤：

步骤3.2.1：选取网络的层数，单层隐含层可以拟合任何“包含从一个有限空间到另一个 有限空间的连续映射”的函数，而本发明以Int_ip作为输入，Out_ip作为输出，三层BP神经网络 可以高效运行；因此，网络包括输入层、单层隐含层和输出层；

步骤3.2.2：选取隐含层神经元个数，BP神经网络的输入层和输出层的神经元个数为Θ， 隐含层神经元个数为：

其中，Θ'为隐含层神经元个数；c取值为1～10的常数；隐含层神经元个数与输入参数和 输出参数的设置有关；隐含层输出，如式(26)所示：

步骤3.2.3：设定训练次数，训练次数设置为100次；

步骤3.2.4：设定学习效率，学习效率是在数据反向传递的过程中起到调整权值的作用， 取值为0.01～0.1之间，本实施例中，设置为0.03；

步骤3.2.5：设定期望最小误差，设定网络期望最小误差为0.00001，当训练结果满足该 精度，停止训练；

通过设定BP神经网络的初始值后，以Int_ip作为网络输入，以Out_ip作为网络期望输出， 随着信号的正向传播，误差反向传播，网络权值的不断调整，直至网络的最终输出与Out_ip的相对误差小于设定的阈值N₁；

步骤3.3：网络结构设置为三层，重建后的H^(RC_u)'的像元邻域空间关系

H^(RC_u)'＝{H_j ^(RC_u)'(x'_i,y'_i)；j＝1,2,…,J} (27)

H_j ^(RC_u)'＝(H_jp ^(RC)'；p∈Ω_u (28)

利用H^(RS)的信息熵和H_jp ^(RC_m)的像元邻域关系对重建图像H^(RC_u)'的像元拟合现象进行优 化；

计算每个波段的重建图像H^(RC_m)在空间域D_i ^(H)内的邻域像元集合H_jp ^(RC_m)与其均值的距 离为：

其中，mean_ip表示空间域D_i ^(H)内邻域像元的均值；dist_jp是空间域D_i ^(H)内像元j与均值mean_ip的差，表示像元j在空间域D_i ^(H)内的空间关系；计算dist_j'p在匹配波段p∈Ω_m中波段均值 mean_dist_j'，如式(31)所示：

其中，P'表示匹配波段总数；

信息熵表征图像单波段的整体信息，通过像元邻域空间关系结合波段信息熵的形式优化 H^(RC_u)'的所有波段为：

其中，H_jp ^(RC),p∈Ω_u表示优化后的光谱值，IE^(RS)为H^(RS)的每个波段的信息熵，IE^(RS)＝ (IE_p ^(RS)；p∈Ω)，信息熵归一化后表示为：

n_IE^(RS)＝(n_IE_p ^(RS)；p∈Ω) (33)

结合H^(RC_m)和H^(RC_u)得到超分辨率重建图像H^(RC)。

为对本发明作进一步详细说明，进行了如下实验：

其中，实验数据与参数设置：

本实验使用的HS数据为GaoFen-5 AHSI获取的空间分辨率为30m的HS遥感图像，该图像具有330个波段，包括可见光(Visible Light,VL)、近红外(Near-Infrared,NIF)和短波红外 (Short Wave Infrared,SWIF)波段，其中193-200和246-262波段的灰度值为0，视为无效波段。 去除这些无效波段后，HS遥感图像保留305个波段。使用的MS数据为哨兵二号(Sentinel-2B) 获取的空间分辨率为10m的MS遥感图像。Sentinel-2B MS遥感图像具有4个波段，包括 VL和NIF。表1和表2分别给出了本发明采用的HS遥感图像和MS遥感图像的波段信息。

表1 GaoFen-5 AHSI信息表

Table1 AHSI information table of GaoFen-5.

表2 Sentinel-2B信息表

Table2 Information table of Sentinel-2B.

如表3所示，Sentinel-2B MS遥感图像和GaoFen-5 HS遥感图像的波段匹配情况。

表3 GaoFen-5和Sentinel-2B波段匹配表

Table3 Band matching table of GaoFen-5 and Sentinel-2B.

实验结果与分析：

本发明选取了尺度为300×300像元的Sentinel-2B图像和GaoFen-5图像进行融合；其中， 图5(a)为Sentinel-2B 1波段对应的图像示意图，图5(b)为Sentinel-2B 2波段对应的图像示意 图，图5(c)为Sentinel-2B 3波段对应的图像示意图，图5(d)为Sentinel-2B 4波段对应的图像 示意图，图5(e)为GaoFen-530波段对应的图像示意图，图5(f)为GaoFen-540波段对应的图 像示意图，图5(g)为GaoFen-565波段对应的图像示意图，图5(h)为GaoFen-5100波段对应 的图像示意图；

对图5(a)和图5(e)使用小波变换进行融合得到附图6(a)的小波融合图像示意图，对附图 5(b)的和图5(f)的GaoFen-5图像使用小波变换进行融合得到附图6(b)的小波融合图像示意图， 对附图5(c)的和图5(g)使用小波变换进行融合得到附图6(c)的小波融合图像示意图，对附图 5(d)和图5(h)使用小波变换进行融合得到附图6(d)的小波融合图像示意图。

图7为Sentinel-2B图像和GaoFen-5图像的BP神经网络重建图像H^(RC_u)；考虑到Sentinel-2B和GaoFen-5的波段匹配信息，图7(a)为H^(RS_u)的15波段对应的图像示意图，附图7(b)为H^(RS_u)的35波段对应的图像示意图，附图7(c)为H^(RS_u)的55波段对应的图像示意图， 附图7(d)为H^(RS_u)的80波段对应的图像示意图，附图7(e)为H^(RS_u)的160波段对应的图像示 意图，附图7(f)为BP神经网络重建图像H^(RC)的15波段对应的图像示意图，附图7(g)为BP 神经网络重建图像H^(RC)的35波段对应的图像示意图，附图7(h)为BP神经网络重建图像H^(RC)的55波段对应的图像示意图，附图7(i)为BP神经网络重建图像H^(RC)的80波段对应的图像示意图，附图7(j)为BP神经网络重建图像H^(RC)的160波段对应的图像示意图，从附图7(f)-(j) 能够明显看出，本发明的算法生成的非匹配波段的重建图像具有更丰富的纹理细节，但是仅 从主观评价的角度，无法细致地区分本发明算法的优劣。因此，对BP神经网络重建图像进 行客观分析。

对像元大小为300×300的GaoFen-5匹配波段的重建图像H^(RC_m)和非匹配波段的重建图像H^(RC_u)进行组合，图8为了本发明算法重构的三维立体图，其波段总数为305，包括VL波段、 NIF波段和SWIF波段，重建后的空间分辨率为10m。

本发明的实验结果综合考虑了主观和客观因素，采用以下四个指标对性能进行评价：

(1)信息熵(Information Entropy,IE)，用来测量图像的复杂性。图像越复杂，它包含的纹 理信息就越丰富，反之越少，信息熵的定义如下：

IE＝-∑ip(gray)logip(gray) (34)

其中，ip(gray)为像元光谱值gray出现在HS遥感图像中的概率。

(2)光谱一阶微分(First-Order Differential of the Spectrum,FODS)，反映HS遥感图像中 光谱曲线在某一波长处的细微变化，有助于分析HS遥感图像地物光谱变换。光谱曲线的一 阶微分定义如下：

其中，fods为光谱一阶微分，p为波段索引，λ_p为波段p的波长，ε为步长，R'_p+ε和R'_p-ε 分别为波长λ_p+ε和λ_p-ε处光谱曲线的反射率。当ε取不同值时，公式(35)能反映任意波长尺度下 光谱曲线的变化。

(3)相关系数(Correlation Coefficient,CC)，反映HS遥感图像与MS遥感图像之间空间信 息的相关性程度，也可以用来描述超分重建图像与HS遥感图像之间空间信息的相关性程度， 相关系数值越大，说明超分重建图像效果越好。相关系数的定义如下：

其中，CC为相关系数，i为HS遥感图像的像元索引，I为HS遥感图像的像元总数，(x'_i,y'_i) 是HS遥感图像的像元i的网格坐标，

是HS遥感图像的像元均值。j为MS遥感图像的像 元索引，J为MS遥感图像的像元总数，(x'_j,y'_j)是MS遥感图像的像元j的网格坐标

是MS 遥感图像的像元均值。

(4)结构相似性(Structural Similarity Index Measure,SSIM)，是一种符合人类直觉的图像质 量评价标准。结构相似性从图像组成的角度将结构信息定义为不依赖于亮度和对比度，反映 场景中目标结构的属性。并且，对于纹理信息变化较为剧烈的区域更敏感，而对于纹理信息 变化缓慢区域的敏感程度较差；结构相似性定义如下：

其中，SSIM为结构相似性，μ为亮度估计的均值，σ为对比度估计的方差，σ_HM用来衡量HS 遥感图像和MS遥感图像结构相似性的协方差，c₁和c₂分别为常数，以避免由于分母为零而 引起的系统误差。

统计HS遥感图像和H^(RC)的CC和SSIM，结果如附图9所示。由附图9(a)-(b)可以看出，结合两项客观评价指标，波段180的重建效果较好。有MS数据参与重构的H^(RC)比HS 遥感图像具有更丰富的空间信息和更高的相似性，在180～240波段H^(RC)与之间的SSIM显著 降低。波段240后的重构效果随着HS遥感图像本身像素灰度的减小而减小，而且重构时模 拟的像元灰度值也会有较大的误差导致衰减。重构结果较差主要体现在附图9(b)中SSIM的波谷。在无效波段附近，传感器的能量采集性能较弱。去除无效波段之后，无效波段附近的波段关联较弱，因此，H^(RC)的重建效果与HS遥感图像在结构上不同。

如图10所示，H和H^(RC)的信息熵示意图，H^(RC)在1-135、180-305波段的信息熵与H基本相同，说明H^(RC)在这些波段的空间信息保留效果较好。H^(RC)在135-180波段具有较弱的信息丰富度，其IE低于H。

结合SSIM，H^(RC)的重建效果大致分为三个阶段：有MS数据参与的高效重建阶段、没有 MS数据参与的中等效率重建阶段和能量损失的低效重建阶段。

如图11所示，HS遥感图像和H^(RC)的光谱梯度变化。对比这些FODS曲线图像，建筑物、 植被和土壤在135波段附近的光谱差异较大，而水体的光谱差异主要体现在NIF波段和160-190波段。

从FODS图可以看出，NIF波段地物H^(RC)的光谱保留度低于其他波段。其中水体的FODS 曲线变化最明显，说明H^(RC)对水体的光谱保留性能较差，其余几种地物的光谱保留性能较好。

HS遥感图像因其具有丰富的光谱信息而被广泛应用于各个领域中，但是由于当前高光谱 传感器本身的局限性，HS遥感图像在具备高光谱分辨率的同时，通常具有较低的空间分辨率， 因此采用超分辨率重建技术改善HS遥感图像的空间分辨率具有重要意义。为了实现HS遥感 图像的全波段的超分重建，本发明提出了一种HS遥感图像和MS遥感图像匹配波段融合、 非匹配波段模拟的重建方法；结论如下：

(1)通过分析Sentinel-2B以及GaoFen-5数据的光谱范围，对HS遥感图像和MS遥感图 像的波段进行划分。定义光谱范围相同的波段为匹配波段，包括VL波段和NIF波段。定义光谱范围不同的波段为非匹配波段，包括VL波段，NIF波段和SWIF波段。重建结果具有GaoFen-5数据相同的光谱范围，且其空间分辨率也得到了提高。

(2)对小波变换理论进行研究，采用小波变换对HS遥感图像匹配波段进行重建。采用 Sentinel-2B MS遥感图像与GaoFen-5 HS遥感图像进行融合重建，通过对比原始HS遥感图像， 能够明显看出纹理细节等方面更清晰。

Claims (7)

1.一种高光谱影像全波段超分重建方法，其特征在于：包括以下步骤：

步骤1：获取HS遥感图像和MS遥感图像；

步骤2：HS遥感图像和MS遥感图像匹配波段的小波融合；

对HS遥感图像进行上采样，使HS遥感图像与MS遥感图像具有相同的空间域和像元尺寸；对MS遥感图像进行小波分解,得到匹配波段的HS遥感图像的小波分解；对于所有的分解层，采用低频融合规则融合第p波段的HS遥感图像H_p ^(RS_m)和第q波段的MS遥感图像M_q的低频系数，高频融合规则融合第p波段的HS遥感图像H_p ^(RS_m)和第q波段的MS遥感图像M_q的三个方向的高频系数；从尺度最大的层次逐层重构到尺度最小的层次进行小波重构，直到最后一层，最后得到融合图像H^(RC_m)；

步骤3：BP神经网络重建不匹配波段；

根据遥感图像中重建后的高空间分辨率HS遥感图像H^(RC)和重采样后的HS遥感图像H^(RS)的像元，匹配波段的融合图像和HS遥感图像对应像元集合{H^(RS_m),H^(RC_m)}作为训练样本，其中，H^(RS_m)为H^(RS)的匹配波段图像，H^(RC_m)为H^(RC)的匹配波段图像；训练BP神经网络重建非匹配波段图像H^(RC_u)；选取BP神经网络参数；网络结构设置为三层；通过像元邻域空间关系结合波段信息熵的形式优化H^(RC_u)'的所有波段；信息熵归一化后表示为n_IE^(RS)；结合H^(RC_m)和H^(RC_u)得到超分辨率重建图像H^(RC)。

2.如权利要求1所述的高光谱影像全波段超分重建方法，其特征在于：所述步骤2的具体过程，包括以下步骤：

步骤2.1：HS遥感图像，如式(1)所示：

H＝{H_i(x'_i,y'_i),i＝1,2,…,I} (1)

其中，i为像元索引；I为像元总数；

为像元i的网格坐标，其中，D_i ^(H)为像元i的分辨率单位面积，D^(H)＝D₁ ^(H)∪…∪D_i ^(H)…∪D_I ^(H)是H的空间域；H_i＝(H_ip；p∈Ω)为像元i的光谱测度矢量，p为HS遥感图像的波段索引，Ω＝{1,2,…,P}，P为波段总数，H_ip为像元i在p波段的光谱值；

MS遥感图像为：

M＝{M_j(x'_j,y'_j),j＝1,2,…,J} (2)

其中，j为像元索引，J为像元总数，

是像元j的网格坐标，D_j ^(M)为像元j的分辨率单位面积，D^(M)＝D₁ ^(M)∪…∪D_j ^(M)…∪D_J ^(M)是M的空间域；M_j＝(M_jq；q∈Λ)为像元j的光谱测度矢量，q为MS遥感图像的波段索引，Λ＝{1,2,…,Q}，Q是波段总数，M_jq为像元j在q波段的光谱值；

fr_p ^(H)为HS遥感图像第p波段的中心频率；B_p ^(H)为HS遥感图像第p波段的光谱范围；fr_q ^(M)为MS遥感图像第q波段的中心频率；B_q ^(M)为MS遥感图像第q波段的光谱范围；HS遥感图像中光谱范围B_p ^(H)位于MS遥感图像光谱范围B_q ^(M)之内的波段为匹配波段，如式(3)所示：

HS遥感图像的光谱范围与MS遥感图像的光谱范围不重叠的波段为非匹配波段，如式(4)所示：

Ω_u＝Ω\Ω_m (4)

重建后的高空间分辨率HS遥感图像为：

H^(RC)＝{H_j ^(RC)(x'_i,y'_i)，j＝1,2,…,J} (5)

其中，H_j ^(RC)＝(H_jp ^(RC)；p∈Ω)是像元j的光谱测度矢量，H_jk ^(RC)为重建后的光谱值；

所述H^(RC)具有与M相同的空间结构与H相同的光谱结构；H的超分重建通过M估计H在匹配和非匹配波段的H^(RC)的光谱值；对于匹配波段

H的超分重建图像表示为：

H^(RC_m)＝{H_j ^(RC_m)(x'_i,y'_i),j＝1,2,…,J} (6)

H_j ^(RC_m)＝(H_jp ^(RC)；p∈Ω_m) (7)

对于非匹配波段Ω_u＝Ω\Ω_m，H的超分重建图像表示为：

H^(RC_u)＝{H_j ^(RC_u)(x'_i,y'_i),j＝1,2,…,J} (8)

H_j ^(RC_u)＝(H_jp ^(RC)；p∈Ω_u) (9)

步骤2.2：对HS遥感图像H＝{H_i(x'_i,y'_i)；i＝1,2,…,I}进行上采样，使HS遥感图像与MS遥感图像具有相同的空间域和像元尺寸；

重采样后的HS遥感图像为：

H^(RS)＝{H_j ^(RS)(x'_j,y'_j)；j＝1,2,…,J} (10)

H_j ^(RS)＝(H_jp ^(RS)；p∈Ω) (11)

根据Ω_u和Ω_m，H^(RS)的匹配波段图像和非匹配波段图像表示为：

H^(RS_m)＝{H_j ^(RS_m)(x'_j,y'_j)；j＝1,2,…,J} (12)

H^(RS_u)＝{H_j ^(RS_u)(x'_j,y'_j)；j＝1,2,…,J} (13)

H_j ^(RS_m)＝(H_jp ^(RS)；p∈Ω_m)andH_j ^(RS_u)＝(H_jp ^(RS)；p∈Ω_u) (14)

步骤2.3：MS遥感图像是索引为(x'_j,y'_j)的二维信号，其中，M_j(x'_j,y'_j)表示MS遥感图像第x'_j行y'_j列的像元值；M_j(x'_j,:)为x'_j行的像元集合，M_j(:,y'_j)为y'_j列的像元集合；MS遥感图像的小波分解为：

其中，ξ_gg为小波分解变换，表示采用低通滤波器g先垂直方向后水平方向的下采样；ξ_hg为小波分解变换，表示先采用低通滤波器g垂直方向下采样，后采用高通滤波器h水平方向下采样；ξ_gh为小波分解变换，表示先采用高通滤波器h垂直方向下采样，后采用低通滤波器g水平方向下采样；ξ_hh为小波分解变换，表示采用高通滤波器h先垂直方向后水平方向的下采样；M_q ^(m)表示第q波段的MS遥感图像M_q经过尺度m次小波分解后的低频系数；Mh_q ^(m+1)表示M_q ^(m)小波分解后水平方向的高频系数；Mv_q ^(m+1)表示M_q ^(m)小波分解后垂直方向的高频系数；Md_q ^{(m +1)}表示M_q ^(m)小波分解后对角方向的高频系数；

匹配波段的HS遥感图像的小波分解为，

其中，H_p ^(RS_m)(m)表示第p波段的HS遥感图像H_p ^(RS_m)经过尺度m次小波分解后的低频系数；Hh_p ^(RS_m)(m+1)表示H_p ^(RS_m)(m)小波分解后水平方向的高频系数；Hv_p ^(RS_m)(m+1)表示H_p ^(RS_m)(m)小波分解后垂直方向的高频系数；Hd_p ^(RS_m)(m+1)表示H_p ^(RS_m)(m)小波分解后对角方向的高频系数；

步骤2.4：对于所有的分解层，采用低频融合规则融合H_p ^(RS_m)和M_q的低频系数，采用高频融合规则融合H_p ^(RS_m)和M_q的三个方向的高频系数；尺度m下的单波段HS-MS遥感图像的低频信息融合表示为：

其中，a'是权重系数，0<a'<1，WF_p ^(RC)(m)表示低频小波系数WF_p ^(RC)尺度m下的小波低频系数；

尺度m下的单波段HS-MS遥感图像的高频信息绝对值最大的融合表示为：

其中，WFh_p ^(RC)(m)表示WF_p ^(RC)尺度m下水平方向的小波高频系数；WFv_p ^(RC)(m)表示WF_p ^(RC)尺度m下垂直方向的小波高频系数；WFd_p ^(RC)(m)表示WF_p ^(RC)尺度m下对角方向的小波高频系数；

步骤2.5：小波重构是小波分解的逆过程，从尺度大的层次逐层重构到尺度小的层次；MS遥感图像的小波重构可表示为：

其中，ζ_gg为小波重构变换，表示采用低通滤波器g先水平方向后垂直方向的上采样；ζ_hg为小波重构变换，表示先采用高通滤波器h水平方向上采样，后采用低通滤波器g垂直方向上采样；ζ_gh为小波重构变换，表示先采用低通滤波器g水平方向上采样，后采用高通滤波器h垂直方向上采样；ζ_hh为小波重构变换，表示采用高通滤波器h先水平方向后垂直方向的上采样；M_q ^(m-1)是m层小波系数进行小波逆变换得到的低频小波系；M_q ^(m-1)作为m层小波逆变换的低频小波系数，参与(m-1)层小波逆变换，当m＝1时，得到重构小波图像M_q ⁽⁰⁾；

HS遥感图像的小波重构可表示为：

从尺度m的分解层开始进行小波重构，m层的一个低频系数和三个方向的高频系数参与重构，重构得到的图像作为m-1层低频小波系数参与m-1层的重构，依次递减，直到最后一层，最后得到融合图像H^(RC_m)。

3.如权利要求1所述的高光谱影像全波段超分重建方法，其特征在于：所述步骤3的具体过程，包括以下步骤：

步骤3.1：根据遥感图像中H^(RC)和H^(RS)的像元，像元j及其邻域像元在p波段上的光谱测量值{H_jp ^(RC)}和{H_jp ^(RS)}有相同的空间关系，匹配波段的融合图像和HS遥感图像对应像元集合{H^(RS_m),H^(RC_m)}作为训练样本训练BP神经网络重建H^(RC_u)；

网络正向传播时的网络输入为：

Int_ip＝(Int_ipι,ι＝1,...,Θ)＝(H_jp ^(RS_m)(x'_j,y'_j) (21)

(x'_j,y'_j)∈D_i ^(H)；p∈Ω_m) (22)

网络正向传播时的网络输出为：

Out_ip＝(Out_ipυ,υ＝1,...,Θ)＝(H_jp ^(RC_m)(x'_j,y'_j) (23)

(x'_j,y'_j)∈D_i ^(H)；p∈Ω_m) (24)

其中，Θ为图像空间域D_i ^(H)个数和波段p的乘积，BP神经网络的样本集{H^(RS_m),H^(RC_m)}表示为{Int_ip,Out_ip}；

步骤3.2：选取BP神经网络参数；

步骤3.3：网络结构设置为三层，重建后的H^(RC_u)'的像元邻域空间关系表示为：

H^(RC_u)'＝{H_j ^(RC_u)'(x'_i,y'_i)；j＝1,2,…,J} (27)

H_j ^(RC_u)'＝(H_jp ^(RC)'；p∈Ω_u) (28)

利用H^(RS)的信息熵和H_jp ^(RC_m)的像元邻域关系对重建图像H^(RC_u)'的像元拟合现象进行优化；

计算每个波段的重建图像H^(RC_m)在空间域D_i ^(H)内的邻域像元集合H_jp ^(RC_m)与其均值的距离为：

其中，mean_ip表示空间域D_i ^(H)内邻域像元的均值；dist_jp是空间域D_i ^(H)内像元j与均值mean_ip的差，表示像元j在空间域D_i ^(H)内的空间关系；计算dist_j'p在匹配波段p∈Ω_m中波段均值mean^_dist_j'，如式(31)所示：

其中，P'表示匹配波段总数；

信息熵表征图像单波段的整体信息，通过像元邻域空间关系结合波段信息熵的形式优化H^(RC_u)'的所有波段为：

其中，H_jp ^(RC),p∈Ω_u表示优化后的光谱值，IE^(RS)为H^(RS)的每个波段的信息熵，IE^(RS)＝(IE_p ^(RS)；p∈Ω)，信息熵归一化后表示为：

n_IE^(RS)＝(n_IE_p ^(RS)；p∈Ω) (33)

结合H^(RC_m)和H^(RC_u)得到超分辨率重建图像H^(RC)。

4.如权利要求3所述的高光谱影像全波段超分重建方法，其特征在于：所述选取BP神经网络参数具体包括以下步骤：

步骤3.2.1：选取网络的层数，以Int_ip作为输入，Out_ip作为输出，选取三层BP神经网络，包括输入层、单层隐含层和输出层；

步骤3.2.2：选取隐含层神经元个数，BP神经网络的输入层和输出层的神经元个数为Θ，隐含层神经元个数为：

其中，Θ'为隐含层神经元个数；c为取值1～10的常数；隐含层神经元个数与输入参数和输出参数的设置有关；

隐含层输出，如式(26)所示：

Hid_ip＝(Hid_ipθ,θ＝1,...,Θ') (26)

步骤3.2.3：设定训练次数；

步骤3.2.4：设定学习效率，学习效率是在数据反向传递的过程中起到调整权值的作用，取值为0.01～0.1；

步骤3.2.5：设定期望最小误差，设定网络期望最小误差为0.00001，当训练结果满足该精度，停止训练；

通过设定BP神经网络的初始值后，以Int_ip作为网络输入，以Out_ip作为网络期望输出，随着信号的正向传播，误差反向传播，网络权值的不断调整，直至网络的最终输出与Out_ip的相对误差小于设定的阈值N₁。

5.如权利要求4所述的高光谱影像全波段超分重建方法，其特征在于：训练次数设置为100次。

6.如权利要求4所述的高光谱影像全波段超分重建方法，其特征在于：设定学习效率为0.01。

7.如权利要求4所述的高光谱影像全波段超分重建方法，其特征在于：设定网络期望最小误差为0.00001。

Images