基于光谱形态和Gram-Schmidt变换约束的高光谱影像变分融
合方法
技术领域
背景技术
[0002] 高光谱影像能够提供多波段且精细的光谱信息，被广泛地应用于森林填图、城市环境模拟和农业
质量检测。然而，光学遥感系统受到入射
能量、卫星
载荷和传输带宽的影响，使得影像的空间
分辨率与光谱分辨率难以同时提高。因此，高光谱遥感影像的空间分辨率通常低于其他影像。遥感影像融合技术可以结合同区域的高光谱影像和高分辨率影像生成理想的高分辨率高光谱影像。
[0003] 在高光谱影像融合中，如何保留原始高光谱影像的光谱信息，并在更精细空间尺寸下生成准确的光谱特征是重要的问题。目前，现有的融合方法被分为成分替换法、多尺度分析法、矩阵分解法和贝叶斯方法。多数融合方法不能较好的实现光谱信息的保真。变分融合方法通过设计影像约束项，建立能量方程并求最优解的方法来重建融合影像。该方法可以生成较好的空间细节，并且平衡影像中的空间信息与光谱信息。但是，变分融合方法需要进一步提高光谱信息保真度，提高融合影像与真实影像的相关性。
发明内容
[0004] 为了解决上述问题，本发明提供了一种基于光谱形态和Gram-Schmidt变换约束的高光谱影像变分融合方法，一种基于光谱形态和Gram-Schmidt变换约束的高光谱影像变分融合方法，主要包括以下步骤：
[0005] S101：获取待融合区域的高分辨率影像及对应的高光谱影像；并将所述高光谱影像进行重
采样，使其空间分辨率和所述高分辨率影像相同，进而得到重采样后的高光谱影像；并将所述重采样后的高光谱影像作为第一次
迭代的融合影像；
[0006] S102：分别计算融合影像和所述高分辨率影像各波段的梯度；
[0007] S103：根据计算得到的所述融合影像各波段的梯度和所述高分辨率影像各波段的梯度，建立空间信息保真项Eg，以增强融合影像的空间细节；
[0008] S104：分别计算所述重采样后的高光谱影像和融合影像的光谱形态
特征向量，并分别计算融合影像中各精细
像素在所述重采样后的高光谱影像中对应的粗像素邻域内的权重；
[0009] S105：根据所述重采样后的高光谱影像和融合影像的光谱形态特征向量和融合影像中各精细像素在所述重采样后的高光谱影像中对应的粗像素邻域内的权重，计算得到光谱形态约束项Es；
[0010] S106：采用Gram-Schmidt变换对所述高分辨率影像和所述重采样后的高光谱影像进行合并处理，得到处理结果Z；并将融合影像与处理结果Z相减，得到相关性约束项Ec；
[0011] S107：根据空间信息保真项Eg、光谱形态约束项Es和相关性约束项Ec，建立能量方程E；
[0012] S108：通过
梯度下降法计算所述能量方程E的最优解，以重建融合影像，得到重建后的融合影像；
[0013] S109：计算所述重建后的融合影像和所述重采样后的高光谱影像之间的光谱
角R；
[0014] S110：判断条件R＜r或者w＞num是否成立？若是，则到步骤S111；否则，将w更新为w+1，并将所述重建后的融合影像作为下一次迭代的融合影像，并返回步骤S102；其中，r为预设的光谱角
阈值，w为迭代次数，w的初始值为1；num为预设的最大迭代次数；
[0015] S111：将当前融合影像作为修正后的最终融合影像，并输出。
[0016] 进一步地，步骤S101中，将未经过重采样的高光谱影像中的像素作为粗像素，重采样过程将每个所述粗像素细分成了多个待修正的精细像素；重采样后的高光谱影像中的精细像素的总数与所述高分辨率影像的像素总数相同；且所述重采样后的高光谱影像和所述高分辨率影像的影像
覆盖面积大小和空间分辨率大小均相同，所述高分辨率影像只有一个波段，而所述重采样后的高光谱影像有多个波段。
[0017] 进一步地，步骤S102中，计算所述高分辨率影像和所述融合影像各波段的梯度时，计算公式如公式(1)所示：
[0018]
[0019] 上式中， 和 分别为融合影像第i个波段的梯度值和高分辨率影像的梯度值，i＝1，2，...，n，n为融合影像的总波段数； 和 分别表示融合影像第i个波段在x方向和2
在y方向上的偏导数； 和 分别表示所述高分辨率影像在x方向和在y方向上的偏导数；ε为预设的残差值。
[0020] 进一步地，步骤S103中，根据计算得到的所述融合影像各波段的梯度和所述高分辨率影像的梯度，采用公式(2)建立空间信息保真项Eg：
[0021]
[0022] 上式中， 表示高分辨率影像的标准向量场；i＝1，2，...，n，n为融合影像的总波段数；Ω表示全部融合影像区域； 为融合影像第i个波段的梯度值。
[0023] 进一步地，步骤S104中，采用公式(3)分别计算所述重采样后的高光谱影像和融合影像的光谱形态特征向量：
[0024] DHi(xj)＝Hi(xj)-mean(H(xj))
[0025] Dui(xj)＝ui(xj)-mean(u(xj)) (3)
[0026] 上式中，DHi(xj)和Dui(xj)分别表示重采样后的高光谱影像中第j个像素在第i个波段的光谱形态特征值和融合影像中第j个像素在第i个波段的光谱形态特征值；Hi(xj)和ui(xj)分别表示重采样后的高光谱影像中第j个像素在第i个波段的值和融合影像中第j个像素在第i个波段的值； 表示所述采样后的高光谱影像中第j个像素在所有波段的均值， 表示所述融合影像中第j个像素在
所有波段的均值，i＝1，2，...，n，j＝1，2，...，m；n和m分别为融合影像的总波段数和像素总个数；公式(3)所涉及的重采样后的高光谱影像中的像素均指重采样后的精细像素；
[0027] 其中，将一个精细像素周围T×T大小的粗像素范围作为该精细像素的邻域，即一个精细像素的邻域中有 个相邻粗像素，再加上该精细像素所在的粗像素本身，共有个粗像素组成了T×T的粗像素范围的邻域，t×t为一个粗像素单位大小，T为预设值，且大于0，为t的整数倍；
[0028] 计算融合影像中某个精细像素x在所述重采样后的高光谱影像中对应的粗像素邻域内的某个粗像素Y的权重，方法为：
[0029] 首先，将该精细像素x对应的粗像素领域内的该粗像素Y投影到所述高分辨率影像中，并计算该粗像素Y的中心的精细像素到精细像素x的欧式距离；
[0030] 然后，分别计算粗像素Y内的各精细像素与精细像素x的差的绝对值，并用所有绝对值之和表示精细像素x与粗像素Y中覆盖地物的相似性程度，即粗像素Y的权重；
[0031] 按照上述方法依次计算精细像素x在所述重采样后的高光谱影像中对应的粗像素邻域内剩余粗像素的权重。
[0032] 进一步地，步骤S105中，根据所述重采样后的高光谱影像和融合影像的光谱形态特征向量和融合影像中各精细像素在所述重采样后的高光谱影像中对应的粗像素邻域内的权重，采用公式(5)计算得到光谱形态约束项Es：
[0033] Es＝∫Ω∫I(Du(x)-Φ·DH(Y))2w(x，Y)dxdY (5)
[0034] 上式中，Ω表示全部融合影像区域；I表示一个精细像素的邻域范围大小；为比例系数；这一约束用于修正融合影像中像素的光谱信息，并生成更高空间分辨率尺度下的光谱特征。
[0035] 进一步地，步骤S106中，相关性约束项Ec的计算如公式(6)所示：
[0036] Ec＝∫Ω(u-Z)2 (6)
[0037] 上式中，Z表示Gram-Schmidt变换结果；融合影像u和Z相减时，u和Z中的对应的像素值相减。
[0038] 进一步地，步骤S107中，能量方程E，如公式(7)所示：
[0039]
[0040] 上式中，γ、α、β和η均为预设的比例系数； 为融合影像第i个波段的梯度值，ui为融合影像第i个波段值；Φ为比例系数；n为融合影像的总波段数；w(x，Y)为精细像素x的邻域粗像素Y的权重；DH(Y)和Du(x)分别表示原始高光谱影像中的粗像素的光谱形态值和融合影像中精细像素的光谱形态特征值；θ表示高分辨率影像的标准向量场。
[0041] 进一步地，步骤S108中，通过梯度下降法计算所述能量方程E的最优解时，迭代方程如公式(8)所示：
[0042]
[0043] 上式中，k为迭代次数，k的初始值为1，取值范围为[1，100]。
[0044] 进一步地，步骤S109中，光谱角R的计算方法为：首先，采用公式(9)计算各像素点之间的光谱角：
[0045]
[0046] 上式中，a为所述重建后的融合影像中的某个精细像素点，b为所述重采样后的高光谱影像中与a对应的精细像素点；采用公式(9)遍历计算所有像素点之间的光谱角；
[0047] 然后，把计算得到的所有光谱角求平均值，得到所述重建后的融合影像和所述重采样后的高光谱影像之间的光谱角R。
[0048] 本发明提供的技术方案带来的有益效果是：本发明所提出的技术方案将影像融合问题转化为能量方程的最优化问题用以求解融合影像；与传统的变分融合方法相比，本方法能够有效提高高光谱影像的空间分辨率，并较好的保持影像中原有的光谱信息。
附图说明
[0049] 下面将结合附图及
实施例对本发明作进一步说明，附图中：
[0050] 图1是本发明实施例中一种基于光谱形态和Gram-Schmidt变换约束的高光谱影像变分融合方法的
流程图；
[0051] 图2是本发明实施例中一种基于光谱形态和Gram-Schmidt变换约束的高光谱影像变分融合方法的执行
框架示意图；
[0052] 图3是本发明实施例中权重计算方法的示意图。
具体实施方式
[0053] 为了对本发明的技术特征、目的和效果有更加清楚的理解，现对照附图详细说明本发明的具体实施方式。
[0054] 本发明的实施例提供了一种基于光谱形态和Gram-Schmidt变换约束的高光谱影像变分融合方法。
[0055] 请参考图1，图1是本发明实施例中一种基于光谱形态和Gram-Schmidt变换约束的高光谱影像变分融合方法的流程图，具体包括如下步骤：
[0056] S101：获取待融合区域的高分辨率影像及对应的高光谱影像；并将所述高光谱影像进行重采样，使其空间分辨率和所述高分辨率影像相同，进而得到重采样后的高光谱影像；并将所述重采样后的高光谱影像作为第一次迭代的融合影像；
[0057] S102：分别计算融合影像和所述高分辨率影像各波段的梯度；
[0058] S103：根据计算得到的所述融合影像各波段的梯度和所述高分辨率影像各波段的梯度，建立空间信息保真项Eg，以增强融合影像的空间细节；
[0059] S104：分别计算所述重采样后的高光谱影像和融合影像的光谱形态特征向量，并分别计算融合影像中各精细像素在所述重采样后的高光谱影像中对应的粗像素邻域内的权重；
[0060] S105：根据所述重采样后的高光谱影像和融合影像的光谱形态特征向量和融合影像中各精细像素在所述重采样后的高光谱影像中对应的粗像素邻域内的权重，计算得到光谱形态约束项Es；
[0061] S106：采用Gram-Schmidt变换对所述高分辨率影像和所述重采样后的高光谱影像进行合并处理，得到处理结果Z；并将融合影像与处理结果Z相减，得到相关性约束项Ec；
[0062] S107：根据空间信息保真项Eg、光谱形态约束项Es和相关性约束项Ec，建立能量方程E；
[0063] S108：通过梯度下降法计算所述能量方程E的最优解，以重建融合影像，得到重建后的融合影像；
[0064] S109：计算所述重建后的融合影像和所述重采样后的高光谱影像之间的光谱角R；
[0065] S110：判断条件R＜r或者w＞num是否成立？若是，则到步骤S111；否则，将w更新为w+1，并将所述重建后的融合影像作为下一次迭代的融合影像，并返回步骤S102；其中，r为预设的光谱角阈值，本发明实施例中取值为0.5，w为迭代次数，w的初始值为1；num为预设的最大迭代次数，本发明实施例中取值为100；
[0066] S111：将当前融合影像作为修正后的最终融合影像，并输出。
[0067] 请参阅图2，图2是本发明实施例中一种基于光谱形态和Gram-Schmidt变换约束的高光谱影像变分融合方法的执行框架示意图；
[0068] 步骤S101中，将未经过重采样的高光谱影像中的像素作为粗像素，重采样过程将每个所述粗像素细分成了多个待修正的精细像素；重采样后的高光谱影像中的精细像素的总数与所述高分辨率影像的像素总数相同；且所述重采样后的高光谱影像和所述高分辨率影像的影像覆盖面积大小和空间分辨率大小均相同，所述高分辨率影像只有一个波段，而所述重采样后的高光谱影像有多个波段。
[0069] 步骤S102中，计算所述高分辨率影像和所述融合影像各波段的梯度时，计算公式如公式(1)所示：
[0070]
[0071] 上式中， 和 分别为融合影像第i个波段的梯度值和高分辨率影像的梯度值，i＝1，2，...，n，n为融合影像的总波段数； 和 分别表示融合影像第i个波段在x方向和在y方向上的偏导数； 和 分别表示所述高分辨率影像在x方向和在y方向上的偏导数；ε2为预设的残差值。
[0072] 步骤S103中，根据计算得到的所述融合影像各波段的梯度和所述高分辨率影像的梯度，采用公式(2)建立空间信息保真项Eg：
[0073]
[0074] 上式中， 表示高分辨率影像的标准向量场；i＝1，2，...，n，n为融合影像的总波段数；Ω表示全部融合影像区域； 为融合影像第i个波段的梯度值。
[0075] 步骤S104中，采用公式(3)分别计算所述重采样后的高光谱影像和融合影像的光谱形态特征向量：
[0076] DHi(xj)＝Hi(xj)-mean(H(xj))
[0077] Dui(xj)＝ui(xj)-mean(u(xj)) (3)
[0078] 上式中，DHi(xj)和Dui(xj)分别表示重采样后的高光谱影像中第j个像素在第i个波段的光谱形态特征值和融合影像中第j个像素在第i个波段的光谱形态特征值；Hi(xj)和ui(xj)分别表示重采样后的高光谱影像中第j个像素在第i个波段的值和融合影像中第j个像素在第i个波段的值； 表示所述采样后的高光谱影像中第j个像素在所有波段的均值， 表示所述融合影像中第j个像素在
所有波段的均值，i＝1，2，...，n，j＝1，2，...，m；n和m分别为融合影像的总波段数和像素总个数；公式(3)所涉及的重采样后的高光谱影像中的像素均指重采样后的精细像素；
[0079] 其中，将一个精细像素周围T×T大小的粗像素范围作为该精细像素的邻域，即一个精细像素的邻域中有 个相邻粗像素，再加上该精细像素所在的粗像素本身，共有个粗像素组成了T×T的粗像素范围的邻域，t×t为一个粗像素单位大小，T为预设值，且大于0，为t的整数倍；
[0080] 计算融合影像中某个精细像素x在所述重采样后的高光谱影像中对应的粗像素邻域内的某个粗像素Y的权重，方法为：
[0081] 首先，将该精细像素x对应的粗像素领域内的该粗像素Y投影到所述高分辨率影像中，并计算该粗像素Y的中心的精细像素到精细像素x的欧式距离；
[0082] 然后，分别计算粗像素Y内的各精细像素与精细像素x的差的绝对值，并用所有绝对值之和表示精细像素x与粗像素Y中覆盖地物的相似性程度，即粗像素Y的权重；
[0083] 按照上述方法依次计算精细像素x在所述重采样后的高光谱影像中对应的粗像素邻域内剩余粗像素的权重。
[0084] 举例说明：
[0085] 若x是待修正的精细像素；且将Y的值设为3，t的值设为1，如图3所示，则Y1...Y9是x在原始高光谱影像中对应的邻域内的9个粗像素；Y5为x所在的粗像素；
[0086] 以粗像素Y2的权重为例：首先，将粗像素Y2投影到所述高分辨率影像中，并计算其中心精细像素p5到精细像素x的欧式距离 在图3中由虚线表示；其中，Ax和Bx分别表示精细像素x和中心精细像素p5的横坐标，AY和BY分别表示精细像素x和中心精细像素p5的纵坐标；
[0087] 然后，分别计算粗像素Y2内的各精细像素p1，...，p9与精细像素x的差的绝对值，用所有绝对值之和表示精细像素x与粗像素Y2中覆盖地物的相似性程度，即粗像素Y2的权重；具体如公式(4)所示：
[0088]
[0089] 上式中，x表示重采样后的高光谱影像中的待修正精细像素，Y表示x的邻域内的粗像素(图3所示，Y1、Y2、...、Y9)，p表示重采样后粗像素Y内的精细像素(图3所示，p1，...，p9)，C(x)＝∫I(∫Ydρ(x，Y)|P(x)-P(p)|dp)-1dY，表示标准化参数；
[0090] 按照上述方法依次计算精细像素x在所述重采样后的高光谱影像中对应的粗像素邻域内剩余粗像素的权重。
[0091] 步骤S105中，根据所述重采样后的高光谱影像和融合影像的光谱形态特征向量和融合影像中各精细像素在所述重采样后的高光谱影像中对应的粗像素邻域内的权重，采用公式(5)计算得到光谱形态约束项Es：
[0092] Es＝∫Ω∫I(Du(x)-Φ·DH(Y))2w(x，Y)dxdY (5)
[0093] 上式中，Ω表示全部融合影像区域；I表示一个精细像素的邻域范围大小(图3所示)； 为比例系数；这一约束用于修正融合影像中像素的光谱信息，并生成更高空间分辨率尺度下的光谱特征。
[0094] 步骤S106中，相关性约束项Ec的计算如公式(6)所示：
[0095] Ec＝∫Ω(u-Z)2 (6)
[0096] 上式中，Z表示Gram-Schmidt变换结果；融合影像u和Z相减时，u和Z中的对应的像素值相减。
[0097] 步骤S107中，能量方程E，如公式(7)所示：
[0098]
[0099] 上式中，γ、α、β、μ和η均为预设的比例系数，本发明实施例中均取值为1； 为融合影像第i个波段的梯度值，ui为融合影像第i个波段值；Φ为比例系数；n为融合影像的总波段数；w(x，Y)为精细像素x的邻域粗像素Y的权重；DH(Y)和Du(x)分别表示原始高光谱影像中的粗像素的光谱形态值和融合影像中精细像素的光谱形态特征值；θ表示高分辨率影像的标准向量场。
[0100] 步骤S108中，通过梯度下降法计算所述能量方程E的最优解时，迭代方程如公式(8)所示：
[0101]
[0102] 上式中，k为迭代次数，k的初始值为1，取值范围为[1，100]。
[0103] 步骤S109中，光谱角R的计算方法为：首先，采用公式(9)计算各像素点之间的光谱角：
[0104]
[0105] 上式中，a为所述重建后的融合影像中的某个精细像素点，b为所述重采样后的高光谱影像中与a对应的精细像素点；采用公式(9)遍历计算所有像素点之间的光谱角；
[0106] 然后，把计算得到的所有光谱角求平均值，得到所述重建后的融合影像和所述重采样后的高光谱影像之间的光谱角R。
[0107] 为突出本发明的创造性，进行如下实验对比说明：
[0108] 实验选择环境1A卫星的高分辨率影像和高光谱影像进行融合实验：影像尺寸均为400×400，高分辨率影像有1个波段，高光谱影像有92个波段。
[0109] 选择的对比方法包含：
[0110] 指导滤波主成分分析(guided filter principal component analysis，GFPCA)；
[0111] Gram-Schmidt变换方法(Gram-Schmidt adaptive，GSA)；
[0112] 调制传递函数(modulation transfer function，MTF_GLP)；
[0113] 基于平滑滤波的强度调制(smoothing filter-based intensity modulation，SFIM)；
[0114] 经典变分融合方法(classic variational method)；
[0115] 波段解耦变分融合方法(band-decoupled variational method，NLVD)。
[0116] 定量评价指数包含：
[0117] 光谱角(Spectral angle mapper，SAM)；
[0118] 均方根误差(Root-mean-square error，RMSE)；
[0119] 全局综合误差(Relative dimensionless global error in synthesis，ERGAS)；
[0120] 相关系数(Correlation coefficient，CC)；
[0121] 全局质量指数(Universal image quality index，UIQI)
[0122] 其中，光谱角计算两幅影像中光谱信息的向量角以评价融合影像的光谱质量，理想值是0；定义是：
[0123]
[0124] 其中，a和b分别是两幅影像中的像素值；
[0125] 均方根误差描述融合影像中
辐射畸变的大小，理想值是0；定义是：
[0126]
[0127] 其中，||A||F是A的Frobenius范数，n是A中的像素个数；
[0128] 全局综合误差结合每个波段的均方根误差进行影像的全局质量评价，理想值是0；定义是：
[0129]
[0130] 其中，d是高分变率影像A与高光谱影像A的空间分辨率之比，m是波段数量，μ是波段均值；
[0131] 相关系数评价融合影像与标准影像的相关性，理想值是1；定义是：
[0132]
[0133] 全局质量指数综合评价了融合影像的相关性损失、
亮度畸变和
对比度畸变，理想值是0；定义是：
[0134]
[0135] 其中，σAB是A和B的协方差，是A的均值。
[0136] 实验结果：
[0137] 用本发明的方法和GSA、GFPCA、MTF_GLP、SFIM、经典变分方法和NLVD的结果进行对比，定量评价指标如表1所示：
[0138] 表1.定量评价结果
[0139]
[0140] 与其他方法相比，本发明方法增强了高光谱影像的空间信息，并在更高空间分辨率下生成了更加准确的光谱信息。本发明方法生成结果的SAM、ERGAS、RMSE、CC、UIQI分别是2.6843、1.5167、270.2535、0.9302和0.8816，都优于参与对比的其他方法。
[0141] 本发明的有益效果是：通过设计光谱形态约束项和相关性保真项，并与空间信息保真项组合建立能量方程，迭代计算能量方程的最优解以重建融合影像；空间信息保真项计算高分影像和高光谱影像各个波段的影像梯度，并相减以增强融合影像的空间信息；光谱形态约束项计算原始高光谱影像和融合影像的光谱形态特征，并得到融合影像中每个像素在高光谱影像中每个邻域像素的权重，建立光谱形态约束项修正融合影像的光谱信息；相关性约束项使用Gram-Schmidt变换对高分影像和高光谱影像进行预处理，将融合影像与预处理结果相减以提高相关性；本方法有效的增强了高光谱影像中包含的空间信息，同时在更高空间分辨率尺度下生成了准确的光谱特征，提高了融合影像与真实影像的相关性；
因此，本方法有相关性高、光谱特征准确、空间信息清晰的优点。
[0142] 以上所述仅为本发明的较佳实施例，并不用以限制本发明，凡在本发明的精神和原则之内，所作的任何
修改、等同替换、改进等，均应包含在本发明的保护范围之内。
