CN110390658B - 基于光谱形态和Gram-Schmidt变换约束的高光谱影像变分融合方法 - Google Patents
基于光谱形态和Gram-Schmidt变换约束的高光谱影像变分融合方法 Download PDFInfo
- Publication number
- CN110390658B CN110390658B CN201910528784.8A CN201910528784A CN110390658B CN 110390658 B CN110390658 B CN 110390658B CN 201910528784 A CN201910528784 A CN 201910528784A CN 110390658 B CN110390658 B CN 110390658B
- Authority
- CN
- China
- Prior art keywords
- image
- pixel
- resampled
- fused
- hyperspectral
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Expired - Fee Related
Links
- 238000001228 spectrum Methods 0.000 title claims abstract description 55
- 238000007500 overflow downdraw method Methods 0.000 title claims abstract description 25
- 230000009466 transformation Effects 0.000 title claims abstract description 25
- 230000003595 spectral effect Effects 0.000 claims abstract description 53
- 230000004927 fusion Effects 0.000 claims abstract description 36
- 238000000034 method Methods 0.000 claims abstract description 27
- 239000013598 vector Substances 0.000 claims description 19
- 238000012952 Resampling Methods 0.000 claims description 9
- 238000012545 processing Methods 0.000 claims description 7
- 125000004432 carbon atom Chemical group C* 0.000 claims description 6
- 238000011478 gradient descent method Methods 0.000 claims description 6
- 238000004364 calculation method Methods 0.000 claims description 5
- 238000012935 Averaging Methods 0.000 claims description 3
- 230000008569 process Effects 0.000 claims description 3
- 239000004576 sand Substances 0.000 claims description 3
- 230000009286 beneficial effect Effects 0.000 abstract description 3
- 238000005457 optimization Methods 0.000 abstract description 2
- 238000011426 transformation method Methods 0.000 abstract description 2
- 230000000877 morphologic effect Effects 0.000 description 4
- 238000010586 diagram Methods 0.000 description 3
- 238000011158 quantitative evaluation Methods 0.000 description 3
- 238000004800 variational method Methods 0.000 description 2
- 230000003044 adaptive effect Effects 0.000 description 1
- 238000004458 analytical method Methods 0.000 description 1
- 238000013398 bayesian method Methods 0.000 description 1
- 230000005540 biological transmission Effects 0.000 description 1
- 230000015572 biosynthetic process Effects 0.000 description 1
- 238000000354 decomposition reaction Methods 0.000 description 1
- 238000001514 detection method Methods 0.000 description 1
- 230000000694 effects Effects 0.000 description 1
- 238000005516 engineering process Methods 0.000 description 1
- 238000002474 experimental method Methods 0.000 description 1
- 238000009499 grossing Methods 0.000 description 1
- 238000013507 mapping Methods 0.000 description 1
- 239000011159 matrix material Substances 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 238000010606 normalization Methods 0.000 description 1
- 230000003287 optical effect Effects 0.000 description 1
- 238000007781 pre-processing Methods 0.000 description 1
- 238000000513 principal component analysis Methods 0.000 description 1
- 238000013441 quality evaluation Methods 0.000 description 1
- 230000005855 radiation Effects 0.000 description 1
- 238000004088 simulation Methods 0.000 description 1
- 238000006467 substitution reaction Methods 0.000 description 1
- 238000003786 synthesis reaction Methods 0.000 description 1
- 238000012546 transfer Methods 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T5/00—Image enhancement or restoration
- G06T5/50—Image enhancement or restoration using two or more images, e.g. averaging or subtraction
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T2207/00—Indexing scheme for image analysis or image enhancement
- G06T2207/10—Image acquisition modality
- G06T2207/10032—Satellite or aerial image; Remote sensing
- G06T2207/10036—Multispectral image; Hyperspectral image
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T2207/00—Indexing scheme for image analysis or image enhancement
- G06T2207/20—Special algorithmic details
- G06T2207/20212—Image combination
- G06T2207/20221—Image fusion; Image merging
Landscapes
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Engineering & Computer Science (AREA)
- Theoretical Computer Science (AREA)
- Image Processing (AREA)
Abstract
本发明提供了一种基于光谱形态和Gram‑Schmidt变换约束的高光谱影像变分融合方法，包括：基于传统的影像变分融合模型，设计了新的光谱形态约束项和相关性约束项；光谱形态约束项采用邻域像素的光谱形态特征和新的权重分配方法，修正由空间分辨率的变化引起的光谱失真；相关性约束项基于Gram‑Schmidt变换方法的生成影像建立约束，提高融合影像与标准参考影像的相关性。本发明的有益效果是：本发明所提出的技术方案将影像融合问题转化为能量方程的最优化问题用以重建融合影像；与传统的变分融合方法相比，本方法能够有效提高高光谱影像的空间分辨率，并较好的保持影像中原有的光谱信息。
Description
技术领域
本发明涉及遥感图像处理数据融合技术领域，尤其涉及一种基于光谱形态和Gram-Schmidt变换约束的高光谱影像变分融合方法。
背景技术
高光谱影像能够提供多波段且精细的光谱信息，被广泛地应用于森林填图、城市环境模拟和农业质量检测。然而，光学遥感系统受到入射能量、卫星载荷和传输带宽的影响，使得影像的空间分辨率与光谱分辨率难以同时提高。因此，高光谱遥感影像的空间分辨率通常低于其他影像。遥感影像融合技术可以结合同区域的高光谱影像和高分辨率影像生成理想的高分辨率高光谱影像。
在高光谱影像融合中，如何保留原始高光谱影像的光谱信息，并在更精细空间尺寸下生成准确的光谱特征是重要的问题。目前，现有的融合方法被分为成分替换法、多尺度分析法、矩阵分解法和贝叶斯方法。多数融合方法不能较好的实现光谱信息的保真。变分融合方法通过设计影像约束项，建立能量方程并求最优解的方法来重建融合影像。该方法可以生成较好的空间细节，并且平衡影像中的空间信息与光谱信息。但是，变分融合方法需要进一步提高光谱信息保真度，提高融合影像与真实影像的相关性。
发明内容
为了解决上述问题，本发明提供了一种基于光谱形态和Gram-Schmidt变换约束的高光谱影像变分融合方法，一种基于光谱形态和Gram-Schmidt变换约束的高光谱影像变分融合方法，主要包括以下步骤：
S101：获取待融合区域的高分辨率影像及对应的高光谱影像；并将所述高光谱影像进行重采样，使其空间分辨率和所述高分辨率影像相同，进而得到重采样后的高光谱影像；并将所述重采样后的高光谱影像作为第一次迭代的融合影像；
S102：分别计算融合影像和所述高分辨率影像各波段的梯度；
S103：根据计算得到的所述融合影像各波段的梯度和所述高分辨率影像各波段的梯度，建立空间信息保真项Eg，以增强融合影像的空间细节；
S104：分别计算所述重采样后的高光谱影像和融合影像的光谱形态特征向量，并分别计算融合影像中各精细像素在所述重采样后的高光谱影像中对应的粗像素邻域内的权重；
S105：根据所述重采样后的高光谱影像和融合影像的光谱形态特征向量和融合影像中各精细像素在所述重采样后的高光谱影像中对应的粗像素邻域内的权重，计算得到光谱形态约束项Es；
S106：采用Gram-Schmidt变换对所述高分辨率影像和所述重采样后的高光谱影像进行合并处理，得到处理结果Z；并将融合影像与处理结果Z相减，得到相关性约束项Ec；
S107：根据空间信息保真项Eg、光谱形态约束项Es和相关性约束项Ec，建立能量方程E；
S108：通过梯度下降法计算所述能量方程E的最优解，以重建融合影像，得到重建后的融合影像；
S109：计算所述重建后的融合影像和所述重采样后的高光谱影像之间的光谱角R；
S110：判断条件R<r或者w>num是否成立；若是，则到步骤S111；否则，将w更新为w+1，并将所述重建后的融合影像作为下一次迭代的融合影像，并返回步骤S102；其中，r为预设的光谱角阈值，w为迭代次数，w的初始值为1；num为预设的最大迭代次数；
S111：将当前融合影像作为修正后的最终融合影像，并输出。
进一步地，步骤S101中，将未经过重采样的高光谱影像中的像素作为粗像素，重采样过程将每个所述粗像素细分成了多个待修正的精细像素；重采样后的高光谱影像中的精细像素的总数与所述高分辨率影像的像素总数相同；且所述重采样后的高光谱影像和所述高分辨率影像的影像覆盖面积大小和空间分辨率大小均相同，所述高分辨率影像只有一个波段，而所述重采样后的高光谱影像有多个波段。
进一步地，步骤S102中，计算所述高分辨率影像和所述融合影像各波段的梯度时，计算公式如公式(1)所示：
上式中，
和
分别为融合影像第i个波段的梯度值和高分辨率影像的梯度值，i＝1,2,…,n，n为融合影像的总波段数；
和
分别表示融合影像第i个波段在x方向和在y方向上的偏导数；
和
分别表示所述高分辨率影像在x方向和在y方向上的偏导数；ε2为预设的残差值。
进一步地，步骤S103中，根据计算得到的所述融合影像各波段的梯度和所述高分辨率影像的梯度，采用公式(2)建立空间信息保真项Eg：
上式中，
表示高分辨率影像的标准向量场；i＝1,2,…,n，n为融合影像的总波段数；Ω表示全部融合影像区域；
为融合影像第i个波段的梯度值。
进一步地，步骤S104中，采用公式(3)分别计算所述重采样后的高光谱影像和融合影像的光谱形态特征向量：
上式中，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)所涉及的重采样后的高光谱影像中的像素均指重采样后的精细像素；
其中，将一个精细像素周围T×T大小的粗像素范围作为该精细像素的邻域，即一个精细像素的邻域中有
个相邻粗像素，再加上该精细像素所在的粗像素本身，共有
个粗像素组成了T×T的粗像素范围的邻域，t×t为一个粗像素单位大小，T为预设值，且大于0，为t的整数倍；
计算融合影像中某个精细像素x在所述重采样后的高光谱影像中对应的粗像素邻域内的某个粗像素Y的权重，方法为：
首先，将该精细像素x对应的粗像素领域内的该粗像素Y投影到所述高分辨率影像中，并计算该粗像素Y的中心的精细像素到精细像素x的欧式距离；
然后，分别计算粗像素Y内的各精细像素与精细像素x的差的绝对值，并用所有绝对值之和表示精细像素x与粗像素Y中覆盖地物的相似性程度，即粗像素Y的权重；
按照上述方法依次计算精细像素x在所述重采样后的高光谱影像中对应的粗像素邻域内剩余粗像素的权重。
进一步地，步骤S105中，根据所述重采样后的高光谱影像和融合影像的光谱形态特征向量和融合影像中各精细像素在所述重采样后的高光谱影像中对应的粗像素邻域内的权重，采用公式(5)计算得到光谱形态约束项Es：
Es＝∫Ω∫I(Du(x)-Φ·DH(Y))2w(x,Y)dxdY (5)
上式中，Ω表示全部融合影像区域；I表示一个精细像素的邻域范围大小；
为比例系数；这一约束用于修正融合影像中像素的光谱信息，并生成更高空间分辨率尺度下的光谱特征。
进一步地，步骤S106中，相关性约束项Ec的计算如公式(6)所示：
Ec＝∫Ω(u-Z)2 (6)
上式中，Z表示Gram-Schmidt变换结果；融合影像u和Z相减时，u和Z中的对应的像素值相减。
进一步地，步骤S107中，能量方程E，如公式(7)所示：
上式中，γ、α、β和η均为预设的比例系数；
为融合影像第i个波段的梯度值，ui为融合影像第i个波段值；Φ为比例系数；n为融合影像的总波段数；w(x,Y)为精细像素x的邻域粗像素Y的权重；DH(Y)和Du(x)分别表示原始高光谱影像中的粗像素的光谱形态值和融合影像中精细像素的光谱形态特征值；θ表示高分辨率影像的标准向量场。
进一步地，步骤S108中，通过梯度下降法计算所述能量方程E的最优解时，迭代方程如公式(8)所示：
上式中，k为迭代次数，k的初始值为1，取值范围为[1,100]。
进一步地，步骤S109中，光谱角R的计算方法为：首先，采用公式(9)计算各像素点之间的光谱角：
上式中，a为所述重建后的融合影像中的某个精细像素点，b为所述重采样后的高光谱影像中与a对应的精细像素点；采用公式(9)遍历计算所有像素点之间的光谱角；
然后，把计算得到的所有光谱角求平均值，得到所述重建后的融合影像和所述重采样后的高光谱影像之间的光谱角R。
本发明提供的技术方案带来的有益效果是：本发明所提出的技术方案将影像融合问题转化为能量方程的最优化问题用以求解融合影像；与传统的变分融合方法相比，本方法能够有效提高高光谱影像的空间分辨率，并较好的保持影像中原有的光谱信息。
附图说明
下面将结合附图及实施例对本发明作进一步说明，附图中：
图1是本发明实施例中一种基于光谱形态和Gram-Schmidt变换约束的高光谱影像变分融合方法的流程图；
图2是本发明实施例中一种基于光谱形态和Gram-Schmidt变换约束的高光谱影像变分融合方法的执行框架示意图；
图3是本发明实施例中权重计算方法的示意图。
具体实施方式
为了对本发明的技术特征、目的和效果有更加清楚的理解，现对照附图详细说明本发明的具体实施方式。
本发明的实施例提供了一种基于光谱形态和Gram-Schmidt变换约束的高光谱影像变分融合方法。
请参考图1，图1是本发明实施例中一种基于光谱形态和Gram-Schmidt变换约束的高光谱影像变分融合方法的流程图，具体包括如下步骤：
S101：获取待融合区域的高分辨率影像及对应的高光谱影像；并将所述高光谱影像进行重采样，使其空间分辨率和所述高分辨率影像相同，进而得到重采样后的高光谱影像；并将所述重采样后的高光谱影像作为第一次迭代的融合影像；
S102：分别计算融合影像和所述高分辨率影像各波段的梯度；
S103：根据计算得到的所述融合影像各波段的梯度和所述高分辨率影像各波段的梯度，建立空间信息保真项Eg，以增强融合影像的空间细节；
S104：分别计算所述重采样后的高光谱影像和融合影像的光谱形态特征向量，并分别计算融合影像中各精细像素在所述重采样后的高光谱影像中对应的粗像素邻域内的权重；
S105：根据所述重采样后的高光谱影像和融合影像的光谱形态特征向量和融合影像中各精细像素在所述重采样后的高光谱影像中对应的粗像素邻域内的权重，计算得到光谱形态约束项Es；
S106：采用Gram-Schmidt变换对所述高分辨率影像和所述重采样后的高光谱影像进行合并处理，得到处理结果Z；并将融合影像与处理结果Z相减，得到相关性约束项Ec；
S107：根据空间信息保真项Eg、光谱形态约束项Es和相关性约束项ec，建立能量方程E；
S108：通过梯度下降法计算所述能量方程E的最优解，以重建融合影像，得到重建后的融合影像；
S109：计算所述重建后的融合影像和所述重采样后的高光谱影像之间的光谱角R；
S110：判断条件R<r或者w>num是否成立；若是，则到步骤S111；否则，将w更新为w+1，并将所述重建后的融合影像作为下一次迭代的融合影像，并返回步骤S102；其中，r为预设的光谱角阈值，本发明实施例中取值为0.5，w为迭代次数，w的初始值为1；num为预设的最大迭代次数，本发明实施例中取值为100；
S111：将当前融合影像作为修正后的最终融合影像，并输出。
请参阅图2，图2是本发明实施例中一种基于光谱形态和Gram-Schmidt变换约束的高光谱影像变分融合方法的执行框架示意图；
步骤S101中，将未经过重采样的高光谱影像中的像素作为粗像素，重采样过程将每个所述粗像素细分成了多个待修正的精细像素；重采样后的高光谱影像中的精细像素的总数与所述高分辨率影像的像素总数相同；且所述重采样后的高光谱影像和所述高分辨率影像的影像覆盖面积大小和空间分辨率大小均相同，所述高分辨率影像只有一个波段，而所述重采样后的高光谱影像有多个波段。
步骤S102中，计算所述高分辨率影像和所述融合影像各波段的梯度时，计算公式如公式(1)所示：
上式中，
和
分别为融合影像第i个波段的梯度值和高分辨率影像的梯度值，i＝1,2,…,n，n为融合影像的总波段数；
和
分别表示融合影像第i个波段在x方向和在y方向上的偏导数；
和
分别表示所述高分辨率影像在x方向和在y方向上的偏导数；ε2为预设的残差值。
步骤S103中，根据计算得到的所述融合影像各波段的梯度和所述高分辨率影像的梯度，采用公式(2)建立空间信息保真项Eg：
上式中，
表示高分辨率影像的标准向量场；i＝1,2,…,n，n为融合影像的总波段数；Ω表示全部融合影像区域；
为融合影像第i个波段的梯度值。
步骤S104中，采用公式(3)分别计算所述重采样后的高光谱影像和融合影像的光谱形态特征向量：
上式中，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)所涉及的重采样后的高光谱影像中的像素均指重采样后的精细像素；
其中，将一个精细像素周围T×T大小的粗像素范围作为该精细像素的邻域，即一个精细像素的邻域中有
个相邻粗像素，再加上该精细像素所在的粗像素本身，共有
个粗像素组成了T×T的粗像素范围的邻域，t×t为一个粗像素单位大小，T为预设值，且大于0，为t的整数倍；
计算融合影像中某个精细像素x在所述重采样后的高光谱影像中对应的粗像素邻域内的某个粗像素Y的权重，方法为：
首先，将该精细像素x对应的粗像素领域内的该粗像素Y投影到所述高分辨率影像中，并计算该粗像素Y的中心的精细像素到精细像素x的欧式距离；
然后，分别计算粗像素Y内的各精细像素与精细像素x的差的绝对值，并用所有绝对值之和表示精细像素x与粗像素Y中覆盖地物的相似性程度，即粗像素Y的权重；
按照上述方法依次计算精细像素x在所述重采样后的高光谱影像中对应的粗像素邻域内剩余粗像素的权重。
举例说明：
若x是待修正的精细像素；且将Y的值设为3，t的值设为1，如图3所示，则Y1…Y9是x在原始高光谱影像中对应的邻域内的9个粗像素；Y5为x所在的粗像素；
以粗像素Y2的权重为例：首先，将粗像素Y2投影到所述高分辨率影像中，并计算其中心精细像素p5到精细像素x的欧式距离
在图3中由虚线表示；其中，Ax和Bx分别表示精细像素x和中心精细像素p5的横坐标，AY和BY分别表示精细像素x和中心精细像素p5的纵坐标；
然后，分别计算粗像素Y2内的各精细像素p1,…,p9与精细像素x的差的绝对值，用所有绝对值之和表示精细像素x与粗像素Y2中覆盖地物的相似性程度，即粗像素Y2的权重；具体如公式(4)所示：
上式中，x表示重采样后的高光谱影像中的待修正精细像素，Y表示x的邻域内的粗像素(图3所示,Y1、Y2、…、Y9),表示重采样后粗像素Y内的精细像素(图3所示，p1,…,p9)，C(x)＝∫I(∫Y dρ(x,Y)|P(x)-P(p)|dp)-1dY，表示标准化参数；
按照上述方法依次计算精细像素x在所述重采样后的高光谱影像中对应的粗像素邻域内剩余粗像素的权重。
步骤S105中，根据所述重采样后的高光谱影像和融合影像的光谱形态特征向量和融合影像中各精细像素在所述重采样后的高光谱影像中对应的粗像素邻域内的权重，采用公式(5)计算得到光谱形态约束项Es：
Es＝∫Ω∫I(Du(x)-Φ·DH(Y))2w(x,Y)dxdY (5)
上式中，Ω表示全部融合影像区域；I表示一个精细像素的邻域范围大小(图3所示)；
为比例系数；这一约束用于修正融合影像中像素的光谱信息，并生成更高空间分辨率尺度下的光谱特征。
步骤S106中，相关性约束项Ec的计算如公式(6)所示：
Ec＝∫Ω(u-Z)2 (6)
上式中，Z表示Gram-Schmidt变换结果；融合影像u和Z相减时，u和Z中的对应的像素值相减。
步骤S107中，能量方程E，如公式(7)所示：
上式中，γ、α、β、μ和η均为预设的比例系数，本发明实施例中均取值为1；
为融合影像第i个波段的梯度值，ui为融合影像第i个波段值；Φ为比例系数；n为融合影像的总波段数；w(x,Y)为精细像素x的邻域粗像素Y的权重；DH(Y)和Du(x)分别表示原始高光谱影像中的粗像素的光谱形态值和融合影像中精细像素的光谱形态特征值；θ表示高分辨率影像的标准向量场。
步骤S108中，通过梯度下降法计算所述能量方程E的最优解时，迭代方程如公式(8)所示：
上式中，k为迭代次数，k的初始值为1，取值范围为[1,100]。
步骤S109中，光谱角R的计算方法为：首先，采用公式(9)计算各像素点之间的光谱角：
上式中，a为所述重建后的融合影像中的某个精细像素点，b为所述重采样后的高光谱影像中与a对应的精细像素点；采用公式(9)遍历计算所有像素点之间的光谱角；
然后，把计算得到的所有光谱角求平均值，得到所述重建后的融合影像和所述重采样后的高光谱影像之间的光谱角R。
为突出本发明的创造性，进行如下实验对比说明：
实验选择环境1A卫星的高分辨率影像和高光谱影像进行融合实验：影像尺寸均为400×400，高分辨率影像有1个波段，高光谱影像有92个波段。
选择的对比方法包含：
指导滤波主成分分析(guided filter principal component analysis,GFPCA)；
Gram-Schmidt变换方法(Gram-Schmidt adaptive,GSA)；
调制传递函数(modulation transfer function,MTF_GLP)；
基于平滑滤波的强度调制(smoothing filter-based intensity modulation,SFIM)；
经典变分融合方法(classic variational method)；
波段解耦变分融合方法(band-decoupled variational method,NLVD)。
定量评价指数包含：
光谱角(Spectral angle mapper,SAM)；
均方根误差(Root-mean-square error,RMSE)；
全局综合误差(Relative dimensionless global error in synthesis,ERGAS)；
相关系数(Correlation coefficient,CC)；
全局质量指数(Universal image quality index,UIQI)
其中，光谱角计算两幅影像中光谱信息的向量角以评价融合影像的光谱质量，理想值是0；定义是：
其中，a和b分别是两幅影像中的像素值；
均方根误差描述融合影像中辐射畸变的大小，理想值是0；定义是：
其中，‖A‖F是A的Frobenius范数，n是A中的像素个数；
全局综合误差结合每个波段的均方根误差进行影像的全局质量评价，理想值是0；定义是：
其中，d是高分变率影像A与高光谱影像A的空间分辨率之比，m是波段数量，μ是波段均值；
相关系数评价融合影像与标准影像的相关性，理想值是1；定义是：
全局质量指数综合评价了融合影像的相关性损失、亮度畸变和对比度畸变，理想值是0；定义是：
其中，σAB是A和B的协方差，
是A的均值。
实验结果：
用本发明的方法和GSA、GFPCA、MTF_GLP、SFIM、经典变分方法和NLVD的结果进行对比，定量评价指标如表1所示：
表1.定量评价结果
与其他方法相比，本发明方法增强了高光谱影像的空间信息，并在更高空间分辨率下生成了更加准确的光谱信息。本发明方法生成结果的SAM、ERGAS、RMSE、CC、UIQI分别是2.6843、1.5167、270.2535、0.9302和0.8816，都优于参与对比的其他方法。
本发明的有益效果是：通过设计光谱形态约束项和相关性保真项，并与空间信息保真项组合建立能量方程，迭代计算能量方程的最优解以重建融合影像；空间信息保真项计算高分影像和高光谱影像各个波段的影像梯度，并相减以增强融合影像的空间信息；光谱形态约束项计算原始高光谱影像和融合影像的光谱形态特征，并得到融合影像中每个像素在高光谱影像中每个邻域像素的权重，建立光谱形态约束项修正融合影像的光谱信息；相关性约束项使用Gram-Schmidt变换对高分影像和高光谱影像进行预处理，将融合影像与预处理结果相减以提高相关性；本方法有效的增强了高光谱影像中包含的空间信息，同时在更高空间分辨率尺度下生成了准确的光谱特征，提高了融合影像与真实影像的相关性；因此，本方法有相关性高、光谱特征准确、空间信息清晰的优点。
以上所述仅为本发明的较佳实施例，并不用以限制本发明，凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。
Claims (7)
1.一种基于光谱形态和Gram-Schmidt变换约束的高光谱影像变分融合方法，其特征在于，包括以下步骤：
S101：获取待融合区域的高分辨率影像及对应的高光谱影像；并将所述高光谱影像进行重采样，使其空间分辨率和所述高分辨率影像相同，进而得到重采样后的高光谱影像；并将所述重采样后的高光谱影像作为第一次迭代的融合影像；
S102：分别计算融合影像和所述高分辨率影像各波段的梯度；
S103：根据计算得到的所述融合影像各波段的梯度和所述高分辨率影像各波段的梯度，建立空间信息保真项Eg，以增强融合影像的空间细节；所述空间信息保真项Eg的建立方法如下：
根据计算得到的所述融合影像各波段的梯度和所述高分辨率影像的梯度，采用公式(2)建立空间信息保真项Eg：
上式中，
表示高分辨率影像的标准向量场；i＝1，2，...，n，n为融合影像的总波段数；Ω表示全部融合影像区域；
为融合影像第i个波段的梯度值；
S104：分别计算所述重采样后的高光谱影像和融合影像的光谱形态特征向量，并分别计算融合影像中各精细像素在所述重采样后的高光谱影像中对应的粗像素邻域内的权重；
S105：根据所述重采样后的高光谱影像和融合影像的光谱形态特征向量和融合影像中各精细像素在所述重采样后的高光谱影像中对应的粗像素邻域内的权重，计算得到光谱形态约束项Es；所述光谱形态约束项Es的计算方法如下：
根据所述重采样后的高光谱影像和融合影像的光谱形态特征向量和融合影像中各精细像素在所述重采样后的高光谱影像中对应的粗像素邻域内的权重，采用公式(5)计算得到光谱形态约束项Es：
Es＝∫Ω∫I(Du(x)-Φ·DH(Y))2w(x，Y)dxdY (5)
上式中，Ω表示全部融合影像区域；I表示一个精细像素的粗像素邻域范围大小；
为比例系数；这一约束用于修正融合影像中像素的光谱信息，并生成更高空间分辨率尺度下的光谱特征；
S106：采用Gram-Schmidt变换对所述高分辨率影像和所述重采样后的高光谱影像进行合并处理，得到处理结果Z；并将融合影像与处理结果Z相减，得到相关性约束项Ec；
S107：根据空间信息保真项Eg、光谱形态约束项Es和相关性约束项Ec，建立能量方程E；所述能量方程E如公式(7)所示：
上式中，γ、α、β、μ和η均为预设的比例系数；
为融合影像第i个波段的梯度值，ui为融合影像第i个波段值；Φ为比例系数；n为融合影像的总波段数；w(x，Y)为精细像素x的邻域粗像素Y的权重；DH(Y)和Du(x)分别表示原始高光谱影像中的粗像素的光谱形态值和融合影像中精细像素的光谱形态特征值；θ表示高分辨率影像的标准向量场；
S108：通过梯度下降法计算所述能量方程E的最优解，以重建融合影像，得到重建后的融合影像；
S109：计算所述重建后的融合影像和所述重采样后的高光谱影像之间的光谱角R；
S110：判断条件R＜r或者w＞num是否成立；若是，则到步骤S111；否则，将w更新为w+1，并将所述重建后的融合影像作为下一次迭代的融合影像，并返回步骤S102；其中，r为预设的光谱角阈值，w为迭代次数，w的初始值为1；num为预设的最大迭代次数；
S111：将当前融合影像作为修正后的最终融合影像，并输出。
2.如权利要求1所述的一种基于光谱形态和Gram-Schmidt变换约束的高光谱影像变分融合方法，其特征在于：步骤S101中，将未经过重采样的高光谱影像中的像素作为粗像素，重采样过程将每个所述粗像素细分成了多个待修正的精细像素；重采样后的高光谱影像中的精细像素的总数与所述高分辨率影像的像素总数相同；且所述重采样后的高光谱影像和所述高分辨率影像的影像覆盖面积大小和空间分辨率大小均相同，所述高分辨率影像只有一个波段，而所述重采样后的高光谱影像有多个波段。
3.如权利要求2所述的一种基于光谱形态和Gram-Schmidt变换约束的高光谱影像变分融合方法，其特征在于：步骤S102中，计算所述高分辨率影像和所述融合影像各波段的梯度时，计算公式如公式(1)所示：
上式中，
和
分别为融合影像第i个波段的梯度值和高分辨率影像的梯度值，i＝1，2，...，n，n为融合影像的总波段数；
和
分别表示融合影像第i个波段在x方向和在y方向上的偏导数；
和
分别表示所述高分辨率影像在x方向和在y方向上的偏导数；ε2为预设的残差值。
4.如权利要求2所述的一种基于光谱形态和Gram-Schmidt变换约束的高光谱影像变分融合方法，其特征在于：步骤S104中，采用公式(3)分别计算所述重采样后的高光谱影像和融合影像的光谱形态特征向量：
上式中，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)所涉及的重采样后的高光谱影像中的像素均指重采样后的精细像素；
其中，将一个精细像素周围T×T大小的粗像素范围作为该精细像素的邻域，即一个精细像素的邻域中有
个相邻粗像素，再加上该精细像素所在的粗像素本身，共有
个粗像素组成了T×T的粗像素范围的邻域，t×t为一个粗像素单位大小，T为预设值，且大于0，为t的整数倍；
计算融合影像中某个精细像素x在所述重采样后的高光谱影像中对应的粗像素邻域内的某个粗像素Y的权重，方法为：
首先，将该精细像素x对应的粗像素领域内的该粗像素Y投影到所述高分辨率影像中，并计算该粗像素Y的中心的精细像素到精细像素x的欧式距离；
然后，分别计算粗像素Y内的各精细像素与精细像素x的差的绝对值，并用所有绝对值之和表示精细像素x与粗像素Y中覆盖地物的相似性程度，即粗像素Y的权重；
按照上述方法依次计算精细像素x在所述重采样后的高光谱影像中对应的粗像素邻域内剩余粗像素的权重。
5.如权利要求2所述的一种基于光谱形态和Gram-Schmidt变换约束的高光谱影像变分融合方法，其特征在于：步骤S106中，相关性约束项Ec的计算如公式(6)所示：
Ec＝∫Ω(u-Z)2 (6)
上式中，Z表示Gram-Schmidt变换结果；融合影像u和Z相减时，u和Z中的对应的像素值相减。
6.如权利要求1所述的一种基于光谱形态和Gram-Schmidt变换约束的高光谱影像变分融合方法，其特征在于：步骤S108中，通过梯度下降法计算所述能量方程E的最优解时，迭代方程如公式(8)所示：
上式中，k为迭代次数，k的初始值为1，取值范围为[1，100]。
7.如权利要求1所述的一种基于光谱形态和Gram-Schmidt变换约束的高光谱影像变分融合方法，其特征在于：步骤S109中，光谱角R的计算方法为：首先，采用公式(9)计算各像素点之间的光谱角：
上式中，a为所述重建后的融合影像中的某个精细像素点，b为所述重采样后的高光谱影像中与a对应的精细像素点；采用公式(9)遍历计算所有像素点之间的光谱角；
然后，把计算得到的所有光谱角求平均值，得到所述重建后的融合影像和所述重采样后的高光谱影像之间的光谱角R。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201910528784.8A CN110390658B (zh) | 2019-06-18 | 2019-06-18 | 基于光谱形态和Gram-Schmidt变换约束的高光谱影像变分融合方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201910528784.8A CN110390658B (zh) | 2019-06-18 | 2019-06-18 | 基于光谱形态和Gram-Schmidt变换约束的高光谱影像变分融合方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN110390658A CN110390658A (zh) | 2019-10-29 |
| CN110390658B true CN110390658B (zh) | 2021-04-27 |
Family
ID=68285562
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201910528784.8A Expired - Fee Related CN110390658B (zh) | 2019-06-18 | 2019-06-18 | 基于光谱形态和Gram-Schmidt变换约束的高光谱影像变分融合方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN110390658B (zh) |
Families Citing this family (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN111444835B (zh) * | 2020-03-26 | 2023-08-04 | 贵阳欧比特宇航科技有限公司 | 一种基于多源遥感数据提取地物空间分布位置的方法 |
| CN115479906A (zh) * | 2022-09-27 | 2022-12-16 | 同济大学 | 基于rgb和高光谱图像融合的碎塑料和微塑料检测方法 |
Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US6208752B1 (en) * | 1998-03-12 | 2001-03-27 | The United States Of America As Represented By The Secretary Of The Navy | System for eliminating or reducing exemplar effects in multispectral or hyperspectral sensors |
| CN102013093A (zh) * | 2010-12-02 | 2011-04-13 | 南京大学 | 一种基于Gram-Schmidt融合和LEGION的高分辨率遥感影像分割方法 |
| CN106384332A (zh) * | 2016-09-09 | 2017-02-08 | 中山大学 | 基于Gram‑Schmidt的无人机影像与多光谱影像融合方法 |
| CN108846329A (zh) * | 2018-05-23 | 2018-11-20 | 江南大学 | 一种基于波段选择和特征融合的高光谱人脸识别方法 |
-
2019
- 2019-06-18 CN CN201910528784.8A patent/CN110390658B/zh not_active Expired - Fee Related
Patent Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US6208752B1 (en) * | 1998-03-12 | 2001-03-27 | The United States Of America As Represented By The Secretary Of The Navy | System for eliminating or reducing exemplar effects in multispectral or hyperspectral sensors |
| CN102013093A (zh) * | 2010-12-02 | 2011-04-13 | 南京大学 | 一种基于Gram-Schmidt融合和LEGION的高分辨率遥感影像分割方法 |
| CN106384332A (zh) * | 2016-09-09 | 2017-02-08 | 中山大学 | 基于Gram‑Schmidt的无人机影像与多光谱影像融合方法 |
| CN108846329A (zh) * | 2018-05-23 | 2018-11-20 | 江南大学 | 一种基于波段选择和特征融合的高光谱人脸识别方法 |
Non-Patent Citations (2)
| Title |
|---|
| Global and local Gram-Schmidt methods for hyperspectral pansharpening;Mauro Dalla Mura等;《2015 IEEE International Geoscience and Remote Sensing Symposium (IGARSS)》;20150731;全文 * |
| 基于稀疏表达及空间信息的高光谱图像分类方法研究;杨京辉;《中国博士学位论文全文数据库 信息科技辑》;20171115(第11期);全文 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN110390658A (zh) | 2019-10-29 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN109102477A (zh) | 一种基于非凸低秩稀疏约束的高光谱遥感图像恢复方法 | |
| Song et al. | Spatiotemporal satellite image fusion through one-pair image learning | |
| Loncan et al. | Hyperspectral pansharpening: A review | |
| CN110428387B (zh) | 基于深度学习与矩阵分解的高光谱与全色图像融合方法 | |
| CN110390658B (zh) | 基于光谱形态和Gram-Schmidt变换约束的高光谱影像变分融合方法 | |
| US8699790B2 (en) | Method for pan-sharpening panchromatic and multispectral images using wavelet dictionaries | |
| CN114119444B (zh) | 一种基于深度神经网络的多源遥感图像融合方法 | |
| CN107958450B (zh) | 基于自适应高斯滤波的全色多光谱影像融合方法及系统 | |
| CN103886559B (zh) | 一种光谱图像处理方法 | |
| JP2019537151A (ja) | 画像処理装置、画像処理方法及び画像処理プログラム | |
| CN111583330B (zh) | 一种多尺度时空马尔可夫遥感影像亚像元定位方法及系统 | |
| CN116310883B (zh) | 基于遥感图像时空融合的农业灾害预测方法及相关设备 | |
| CN111383203B (zh) | 基于分区域拟合的全色与多光谱遥感图像融合方法 | |
| CN115170406A (zh) | 一种光学影像与sar强度影像高精度融合方法 | |
| McCarthy et al. | Component-separated, CIB-cleaned thermal Sunyaev-Zel’dovich maps from Planck PR4 data with a flexible public needlet ILC pipeline | |
| Todarello et al. | Constraints on the origin of the radio synchrotron background via angular correlations | |
| CN116883799A (zh) | 成分替换模型引导的高光谱图像深度空间光谱融合方法 | |
| CN110926611A (zh) | 一种应用于压缩感知光谱成像系统的噪声抑制方法 | |
| CN111598115A (zh) | 一种基于交叉皮质神经网络模型的sar影像融合方法 | |
| CN101779112B (zh) | 估计光学系统的波前的至少一个变形或者通过所述光学系统观察到的物体的方法以及相关设备 | |
| CN115100075A (zh) | 基于光谱约束和残差注意力网络的高光谱全色锐化方法 | |
| CN111046844B (zh) | 一种基于邻域选取约束的高光谱图像分类方法 | |
| CN109410165B (zh) | 一种基于分类学习的多光谱遥感图像融合方法 | |
| CN110097501B (zh) | 一种基于非局部均值梯度稀疏正则化的ndvi图像融合方法 | |
| Zeng | Low-light image enhancement algorithm based on lime with pre-processing and post-processing |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant | ||
| CF01 | Termination of patent right due to non-payment of annual fee |
Granted publication date: 20210427 |