技术领域
[0001] 本发明涉及图像融合的技术领域，尤其是指一种基于全变分的耦合非负矩阵分解的图像融合方法。
相关背景技术
[0002] 高光谱图像因其丰富的光谱信息在精准农牧业、生物医药、地质灾害预测等领域都扮演着十分重要的角色。但是，由于光学传感器的局限性等众多因素，多光谱图像虽然在空间分辨率上表现出色，但其光谱分辨率相对较低；而相比之下，高光谱图像具有更高的光谱分辨率，能够捕捉更细致的光谱信息，但其空间分辨率通常较低。这极大限制了高光谱图像在许多领域的应用。因此，研究者们通常将同一场景下的高空间分辨率多光谱图像(High Spatial Resolution Multispectral Image，简称：HR‑MSI)和低空间分辨率高光谱图像(Low Spatial Resolution Hyperspectral Image，简称：LR‑HSI)进行融合，获得高空间分辨率高光谱图像(High Spatial Resolution Hyperspectral Image，简称：HRHSI)，以达到提高高光谱图像空间分辨率的目的。
[0003] 近年来，随着越来越多的研究人员对高光谱图像融合问题产生兴趣，学术界相继提出了多种提高高光谱数据空间分辨率的方法。大致可分为如下几类：利用线性光谱解混(LSU)技术的方法、基于非负矩阵分解的方法、基于张量分解的方法和基于深度学习的方法。
[0004] 虽然现在的融合算法表现出了良好的效果，但在某些方面仍存在一定的不足之处。如，基于张量的方法中，虽然CP分解和Tucker分解展现出了卓越的图像恢复能力，并且以二者为基础的算法框架也具有良好的性能。但是这些方法并没有充分挖掘潜在因子的物理含义，从而难以利用高光谱图像中丰富的先验知识来增强融合算法的性能；基于深度学习的方法仍有两个缺点，首先这类方法需要海量数据来训练神经网络，以获得融合网络的最佳参数，其次深度卷积神经网络泛化能力较弱，在一类数据下训练得到的卷积神经网络可能并不适用于另一类型的数据。
具体实施方式
[0061] 下面结合实施例及附图对本发明作进一步详细的描述，但本发明的实施方式不限于此。
[0062] 如图1所示，本实施例公开了一种基于全变分的耦合非负矩阵分解的图像融合方法，包括以下步骤：
[0063] 步骤1：将低分辨率高光谱图像X和高分辨率多光谱图像Y看作是由目标图像Z(即要恢复的高分辨率高光谱图像)在空间和光谱两个维度下采样得到的：
[0064] Z＝WH+N
[0065] X＝ZS+Es
[0066] Y＝RZ+Er
[0067] 其中， 表示实数集，p表示高光谱传感器光谱通道数，q表示高光谱图像的像素数，m表示多光谱传感器光谱通道数，n表示多光谱图像的像素数；S是空间下采样矩阵，每一列向量表示点扩展函数(PSF)从高光谱图像到多光谱图像逐个像素值的变换，每个PSF都是归一化的，R是光谱下采样矩阵，每一行向量表示从高光谱传感器到多光谱传感器的光谱响应函数的变换，S和R都是稀疏矩阵。N，Es，Er分别代表残差矩阵，W是端元矩阵，H是丰度矩阵，则有：
[0068] X≈WHh
[0069] Y≈WmH
[0070] 式中，Hh≈HS，Wm≈RW，Hh表示空间下采样的丰度矩阵，Wm表示光谱下采样的端元矩阵。
[0071] 步骤2：经典的耦合非负矩阵分解模型，将高光谱图像和多光谱图像进行分解，建立以Frobenius平方范数为代价函数的模型，然后通过空间光谱退化模型将两个分解模型相结合，得到耦合的非负矩阵分解模型。但它是一个典型的病态问题，即想要得到的理想图像(高分辨率高光谱图像)可能与它的空间和光谱退化形式(可观测到的高光谱和多光谱图像)不是唯一反演的，这就使得耦合非负矩阵分解是一个不适定的逆问题。因此，为了得到更合理的解空间，需要引入其它的先验信息来约束优化模型。所以像贝叶斯融合准则需要先验信息一样，CNMF也需要正则化约束。本发明将端元最短距离正则、丰度稀疏正则和全变分正则三项正则化添加到模型中，得到联合正则化图像融合模型：
[0072] 上述耦合非负矩阵分解模型为：
[0073]
[0074]
[0075] 式中，表示“定义为”，||·||F表示Frobenius范数，W≥0和H≥0都是非负约束，是和为1的约束，Hi表示第i个端元在单个像素中的丰度分数，K是端元个数。
[0076] 将端元最短距离正则、丰度稀疏正则和全变分正则这三项正则化添加到上述耦合非负矩阵分解模型中，得到联合正则化图像融合模型，表达形式如下：
[0077]
[0078] s.t.W≥0,H≥0
[0079] 式中， 表示端元最短距离正则项，λ1表示正则化参数，P是一个投影矩阵；λ2||H||1表示丰度稀疏正则项，λ2表示正则化参数，||·||1表示1范数；λ3||DxH||1+λ3||DyH||1表示全变分正则项，λ3表示正则化参数，Dx和Dy分别表示沿水平和垂直方向的二维有限差分算子，由下式给出：
[0080]
[0081]
[0082] 步骤3：使用交替方向乘子法ADMM对联合正则化图像融合模型进行优化求解。
[0083] 首先将模型分解为两个关于W和H的子问题，具体如下：
[0084] 关于W的子问题：
[0085]
[0086] 式中， 表示W≥0。
[0087] 关于H的子问题：
[0088]
[0089] 式中， 表示H≥0。
[0090] 然后对两个子问题分别进行求解。
[0091] 求解关于W的子问题：
[0092] 首先引入新的辅助变量U1，令U1＝W，则可将公式转化为：
[0093]
[0094] 推导出上述问题的增广拉格朗日函数为：
[0095]
[0096] 式中，λ1＞0为正则化参数，ρ1＞0为惩罚参数，A1表示拉格朗日乘子，进一步对其优化分解，其给出为：
[0097]
[0098] 式中，Wk+1表示W的第k+1次求解结果，U1k和U1k+1分别表示U1的第k次和第k+1次求解k k+1结果，A1和A1 分别表示A1的第k次和第k+1次求解结果
[0099] 其中包含W的子问题为：
[0100]
[0101] 该问题是一个二次优化问题，它具有一个唯一的解，其解相当于计算以下的Sylvester矩阵方程:
[0102]
[0103] 可以使用共轭梯度法(CG)来解决上述问题。
[0104] 对于H的子问题：引入新的变量V1、V2、V3、V4，令V1＝H,V2＝DxH,V3＝DyH,V4＝H，则可将子问题转化为：
[0105]
[0106] s.t.V1＝H,V2＝DxH,V3＝DyH,V4＝H
[0107] 推导出上述问题的增广拉格朗日函数为：
[0108]
[0109] 式中，λ2＞0,λ3＞0为正则化参数，α1＞0,α2＞0,α3＞0,α4＞0为惩罚参数，B1,B2,B3,B4为拉格朗日乘子，使用ADMM算法来解决上述问题，其解空间为：
[0110]
[0111] 式中，i＝1,2,3,4；Hk+1表示H的第k+1次求解结果，Vik和Vik+1分别表示Vi的第k次和k k+1第k+1次求解结果，Bi和Bi 分别表示Bi的第k次和第k+1次求解结果。
[0112] 其中包含H的子问题为：
[0113]
[0114] 该问题是一个二次优化问题，它具有一个唯一的解，相当于计算如下的矩阵方程：
[0115]
[0116] 式中，·表示矩阵乘法，[]‑1表示矩阵的逆。
[0117] 步骤4：将步骤3中优化求得的W和H使用公式：Z＝WH+N重构得到最终的目标图像，即高分辨率高光谱图像。
[0118] 本发明在验证所提方法的有效性过程中，使用了三个数据集：Pavia University数据集，Washington DC数据集和青铜峡稻田数据集，如图2中(a)和(b)分别表示Pavia University数据集的高光谱图像和多光谱图像，如图2中(c)和(d)分别表示Washington DC数据集的高光谱图像和多光谱图像，如图2中(e)和(f)分别表示青铜峡稻田数据集的高光谱图像和多光谱图像。
[0119] Pavia University数据集是由意大利反射光学系统成像光谱仪(ROSIS)光学传感器在帕维亚大学市区获得的。图像尺寸为610×340×115，空间分辨率为1.3m，我们在去除水蒸气吸收波段后，波段数减少到102个。鉴于下采样的缘故，我们仅将左上角256×256×102的子图像作为参考图像应用于实验中。在Pavia University数据集上的融合结果展示在图3，下表展示了本发明方法的数值评价指标结果：
[0120] Method PSNR SAM CC RMSE ERGAS Time(s)Best ∞ 0 1 0 0 0
HySure 36.2735 3.0369 0.9831 0.0183 2.5525 126.54
FUSE 32.5423 5.6592 0.9632 0.0204 1.9473 2.36
CNMF 36.6728 4.1959 0.9856 0.0163 1.5793 8.58
LRPL 40.0635 3.0091 0.9912 0.0105 1.3094 96.04
LRSR 39.0126 3.7336 0.9897 0.0127 1.6137 58.43
Proposed 42.2519 1.8606 0.9953 0.0097 1.212 8.76
[0121] Washington DC数据集是由HYDICE传感器采集的华盛顿购物广场的图像，截取1280×307大小的图像进行标注得到的。图像空间分辨率为2.5m，共有210个波段，我们截取大小为256×256×191的子图像作为参考图像进行实验。在Washington DC数据集上的融合结果展示在图4，下表展示了本发明方法的数值评价指标结果：
[0122]Method PSNR SAM CC RMSE ERGAS Time(s)
Best ∞ 0 1 0 0 0
HySure 35.5826 3.0416 0.9786 0.0137 2.1348 35.26
FUSE 34.2647 2.8652 0.9676 0.0118 2.6019 11.94
CNMF 36.2071 2.6853 0.9592 0.0103 1.7206 14.08
LRPL 37.8684 2.1863 0.9913 0.0098 1.5492 337.65
LRSR 38.2639 1.9983 0.9923 0.0094 1.4937 310.07
Proposed 38.9984 2.0992 0.9929 0.0089 1.2445 19.11
[0123] 青铜峡稻田数据集是在2021年7月15日，在宁夏回族自治区吴忠市青铜峡市瞿靖镇瞿靖村的实验稻田由大疆无人机搭载高光谱传感器Pika‑L和多光谱传感器Micasence的飞行活动中获取的场景。采集数据时，无人机飞行速度3m/s，飞行高度120m。高光谱传感器获得的图像尺寸大小为5376×1492×150，图像有150个波段，在去除HSI中一些噪声大的波段，将剩下的60个波段作为实验研究的对象，得到一个包含房屋、树木和水稻田等地物信息，尺寸大小为512×512×60的HSI子图像。利用多光谱传感器获得图像尺寸大小为16103×9677×3，其空间分辨率为5.2cm。我们在相同的位置取数据大小为512×512×3的子图像作为MSI。截取出来的图像面积大概有整幅图像的十分之一。在青铜峡稻田数据集上的融合结果展示在图5，下表展示了本发明方法的数值评价指标结果：
[0124]
[0125]
[0126] 从图3、图4、图5和上述数值评价表综合来看，本发明方法在保留高光谱图像空间和光谱信息的同时，显著增强了图像的细节，实现了高分辨率高光谱图像的重构，在提升融合结果质量方面表现良好。
[0127] 上述实施例为本发明较佳的实施方式，但本发明的实施方式并不受上述实施例的限制，其他的任何未背离本发明的精神实质与原理下所作的改变、修饰、替代、组合、简化，均应为等效的置换方式，都包含在本发明的保护范围之内。
