CN111223049B - 一种基于结构-纹理分解的遥感图像变分融合方法 - Google Patents

Abstract

本发明公开了一种基于结构‑纹理分解的多光谱与全色图像变分融合方法，基于遥感领域中PAN图像和HRMS图像之间结构分量和纹理分量的相似性，分别对全局信息和局部空间细节(高阶信息)进行保留。通过梯度稀疏度来确保PAN图像和HRMS图像的结构相似度，根据低轶属性来描述其纹理相似度，提出我们的变分模型。然后通过交替方向乘子法求解该模型，从而得到具有高分辨率的多光谱图像。实验表明：本发明提出的方法考虑到PAN图像和HRMS图像之间结构分量和纹理分量的相似性，使用不同的模拟数据集和真实数据集进行融合实验，与对比方法相比，我们的融合图像的各项指标均最接近于理想值，视觉上最接近于参考图像。

Description

一种基于结构-纹理分解的遥感图像变分融合方法

技术领域

本发明属于图像融合领域，涉及一种基于结构-纹理(cartoon-texture)分解的多光谱与全色图像变分融合方法，适用于由低分辨率多光谱图像与高分辨率全色图像融合得到高分辨率多光谱图像的应用场景。

背景技术

在遥感领域，通过卫星获取的地球图像的数据具有重要的应用价值，例如应用于土地覆盖分类，灾害检测和植被绘图等。由于卫星的存储与传输负担等限制，通常使用不同的传感器分别获得同一画面的光谱信息丰富但空间分辨率低的多光谱图像(LRMS)和光谱信息单一但空间分辨率高的全色图像(PAN)，并通过对其进行融合来获得同时具有丰富的光谱信息和高分辨率的多光谱图像(HRMS)，以便有利于对信息的分析和理解，该图像融合技术也被称为全锐化技术。

基于不同的融合原理，全锐化方法通常被分为4类，包括成分替换方法(componentsubstitution，CS)，多分辨率分析法(multi-resolution analysis，MRA)，基于深度卷积神经网络(deep convolutional neuron network-based，DCNN)的全锐化方法，以及变分法。CS方法将MS图像映射到一个新的空间，用PAN图像替换其中的空间分量后逆变换到原空间。MRA方法是通过多分辨率分解提取出PAN图像的细节然后将其注入MS波段中，这两类全锐化方法均可以得到分辨率较高的多光谱图像，但通常遭受严重的光谱失真或局部伪像。随着深度学习理论的不断发展，研究人员们开始探索基于深度学习的全锐化方法(deeplearning-based pansharpening method，DLP)，并展现出优异的性能。但这类方法需要大量的计算资源和训练数据，因此在资源受限的情况下缺乏灵活性。近年来，变分方法依靠其优异的融合效果成为了一个广受欢迎的融合方法，这类方法的核心为基于某些假设和先验信息建立目标能量函数。例如，C.Chen等提出引入L_2,1约束来鼓励动态梯度稀疏性和群体稀疏性的SIRF方法，较好的提升了融合质量。X Fu等提出一种基于局部梯度约束(LocalGradient Constraint，LGC)的全锐化变分模型，该方法考虑PAN和HRMS图像在不同局部块和光谱波段中的梯度差异来进行局部梯度约束，实现了高质量的锐化。

在已提出的众多变分方法中，考虑到PAN图像和MS图像取自同一场景，它们具有共同的几何特征，因此通常利用梯度相似性来对取自同一场景的PAN图像和MS图像的几何结构一致性进行约束。然而梯度算子只能获得图像的一阶边缘信息，因此融合过程中全色图像的高阶空间信息并未得到有效的利用。

发明内容

本发明针对现有变分方法的不足，提供一种新的将卡通-纹理分解与二阶导变分模型结合的变分全锐化方法。

本发明所采用的技术方案是：一种结构-纹理分解与二阶导变分模型结合的变分全锐化方法。图1展示了一张全色图像的梯度图像以及结构-纹理分量，其中结构分量包含了全局结构信息，而纹理分量中包含着一些梯度图像中不存在的图案细节(即高阶信息)，这对于全锐化过程中的细节信息提取更有利。基于PAN图像和HRMS图像之间结构分量和纹理分量的相似性，我们分别使用梯度稀疏性和低轶性进行先验约束来保留空间信息，通过数据保真项保留MS图像的光谱信息，通过交替方向乘子法(alternative directionmultiplier method，ADMM)制定并求解变分模型，从而获得具有高分辨率的多光谱图像。

步骤1，选取同一时间同一画面的经过配准处理的多光谱图像和全色图像，低分辨率多光谱图像的记为M，全色图像记为P，融合得到的高分辨率多光谱图像记为F，全色图像与多光谱图像之间的分辨率之比为c；

步骤2，构建数据保真项J(F)来保留低分辨率多光谱图像M的光谱信息，构建先验约束项K(F)以确保融合图像F与全色图像P之间的空间结构一致性，从而获得目标能量函数：

argmin_FJ(F)+K(F) (1)

步骤3：基于步骤2的目标能量函数，通过交替方向乘子法来求解优化问题，从而得到融合图像。

进一步的，步骤2的具体实现包括以下子步骤，

步骤2.1，构建数据保真项J(F)，来保留低分辨率多光谱图像M的光谱信息，低分辨率多光谱图像M被视为高分辨率多光谱图像F的下采样模糊版本，即M＝ΨF+n，其中Ψ表示下采样运算符，n表示均值为0的高斯噪声；因此，为了保留LR MS图像的光谱信息，构造保真项为：

其中

表示Frobenius范数；

步骤2.2，基于全色图像和相应的高分辨率多光谱图像之间结构分量和纹理分量的相似性，在结构-纹理空间构建先验约束项K(F)；具体原理如下：

融合图像F表示为纹理分量R和结构分量S的组合：

F＝R+S (3)

其中结构分量S代表的是图像中的全局结构信息，即大面积平滑区域和锐利的边缘，纹理分量R代表的是图像中的小尺寸细节内容；

选择使用核范数来对纹理分量进行约束，因此，先验约束项K(F)：

其中，

和

表示将U和V分别复制到b个波段，U和V为全色图像P的结构分量和纹理分量，由图像分解算法获得；R和S分别为F的结构和纹理分量，

为梯度算子，||·||_2,1表示群体稀疏；||·||_*表示核范数，ρ和τ为平衡各项的参数；

为块映射函数，用于将矩阵

无重叠的分割为r×r的块并将每个块重新排列成长度为r²的列向量；

步骤2.3：基于子步骤2.1和2.2，提出的目标能量函数模型为：

进一步的，步骤3的具体实现包括以下子步骤，

步骤3.1：对于步骤2得到的目标能量函数模型，通过交替方向乘子法来求解该优化问题，通过引入辅助变量Y，(5)的增广拉格朗日函数表示为：

其中λ为拉格朗日参数，(6)式的优化问题被分解成4个子问题，分别进行求解；

步骤3.2：求解F子问题：

其中下标t表示迭代次数，上述问题相当于：

其中I是单位矩阵，等式(8)通过高斯消去法求解；

步骤3.3：求解R子问题：

等式(9)通过快速迭代收缩阈值算法来解决；

步骤3.4：求解S子问题：

等式(10)是低轶优化问题，通过奇异值阈值算法解决；

步骤3.5：求解Y子问题，该子问题由下式直接更新：

Y_t+1＝Y_t+λ(F_t-R_t-S_t) (11)。

本发明主要针对由低分辨率多光谱图像和高分辨率全色图像融合得到高分辨率多光谱图像的应用需求。考虑到在图像全锐化过程中，梯度算子只能获得图像的一阶边缘信息，导致融合过程中全色图像的高阶空间信息并未得到有效的利用。而图像的结构-纹理分解的理论中，结构分量包含了全局结构信息，而纹理分量中包含着一些梯度图像中不存在的图案细节(即高阶信息)，这对于全锐化过程中的细节信息提取更有利。因此，基于PAN图像和HRMS图像之间结构分量和纹理分量的相似性，我们提出一种卡通-纹理分解与二阶导变分模型结合的变分全锐化方法。我们的模型分别使用梯度稀疏性和低轶性进行先验约束来保留空间信息，通过数据保真项保留MS图像的光谱信息，通过交替方向乘子法(alternative direction multiplier method，ADMM)制定并求解模型，从而获得具有高分辨率的多光谱图像。

附图说明

图1是本发明对于全色图像的梯度图像和结构-纹理分量的比较。

图2是本发明对全色图像和高分辨率多光谱图像的结构和纹理分量的相似性描述，(a)其结构分量的梯度比较(采用来自同一行的数据进行演示)，(b)为纹理分量的比较，其中方框中的块具有相似的结构。

具体实施方式

为了便于本领域普通技术人员理解和实施本发明，下面结合附图及实施例对本发明做进一步的详细描述，应当理解的是，此处所描述的实施例仅用于解释本发明，并不用于限定本发明。

我们所使用的数据来自高分一号卫星数据集，该数据集的MS和PAN图像的空间分辨率分别为8m和2m，MS图像均具有4个波段，即蓝色，绿色，红色，和近红外波段。为了定量分析融合的结果，基于Wald合成协议[1]，我们将经过下采样后得到的MS(75×75×4)和PAN(300×300)数据进行全锐化，并将初始MS图像作为参考图像(300×300×4)。本实施提供一种基于结构-纹理分解的多光谱与全色图像变分融合方法，包括以下步骤：

步骤1：选取同一时间同一画面的经过配准等预处理的多光谱图像和全色图像。低分辨率多光谱图像的记为M，全色图像记为P，融合得到的高分辨率多光谱图像记为F。

在本实施例中，选用的高分一号卫星数据，全色图像与多光谱图像之间的分辨率之比为4。

步骤2：构建数据保真项J(F)来保留低分辨率多光谱图像M的光谱信息，构建先验约束项K(F)以确保融合图像F与全色图像P之间的空间结构一致性。从而获得目标能量函数：

argmin_FJ(F)+K(F) (1)

步骤2的具体实现包括以下子步骤：

步骤2.1：构建数据保真项J(F)，来保留低分辨率多光谱图像M的光谱信息。低分辨率多光谱图像M通常被视为高分辨率多光谱图像F的下采样模糊版本，即M＝ΨF+n，其中Ψ表示下采样运算符，n表示均值为0的高斯噪声。因此，为了保留LR MS图像的光谱信息，我们的保真项为：

其中

表示Frobenius范数。

步骤2.2：基于全色图像和相应的高分辨率多光谱图像之间结构分量和纹理分量的相似性，在结构-纹理空间构建先验约束项K(F)。具体原理如下：

融合图像F可以表示为纹理分量R和结构分量S的组合：

F＝R+S (3)

其中结构分量S代表的是图像中的全局结构信息，即大面积平滑区域和锐利的边缘。纹理分量R代表的是图像中的小尺寸细节内容。

图2展示了我们选取的一组PAN图像与相同分辨率的HRMS图像的结构-纹理分量的相似性关系。首先我们在PAN图像的结构分量图像和HRMS图像的四个波段的结构分量图像上随机扫描相同位置的一条线，并在图2(a)中显示出该条线上的梯度值，该图表明，PAN图像和HRMS图像的结构分量的梯度差异较小，因此，我们选择使用梯度算子来制定结构分量规范。图2(b)中展示了一组HRMS图像的三个波段的纹理分量与PAN图像纹理分量的残差图像，由于遥感图像中存在许多自相似结构，如建筑物，道路等，因此这些残差图像中存在许多重复的纹理，从而具有低轶性，我们选择使用核范数来对纹理分量进行约束。因此，先验约束项K(F)：

其中，

和

表示将U和V分别复制到b个波段。U和V为全色图像P的结构分量和纹理分量，由图像分解算法获得[2]。R和S分别为F的结构和纹理分量。

为梯度算子。||·||_2,1表示群体稀疏[3]；||·||_*表示和范数，通常用于低轶近似。ρ和τ为用于平衡各项的参数。

为块映射函数，用于将矩阵

无重叠的分割为r×r的块并将每个块重新排列成长度为r²的列向量，其中m和n为矩阵T的行数和列数。在本实施例中，b＝4，ρ＝0.05，τ＝1200，m＝300，n＝300，r＝30。

步骤2.3：基于子步骤2.1和2.2，我们所提出的目标能量函数模型为：

步骤3：基于步骤2的目标能量函数，通过交替方向乘子法(alternativedirection multiplier method，ADMM)来求解优化问题，从而得到融合图像。步骤3的具体实现包括以下子步骤：

步骤3.1：对于步骤2得到的模型，我们通过交替方向乘子法来求解该优化问题。通过引入辅助变量Y，(5)的增广拉格朗日函数表示为：

其中λ为拉格朗日参数，在本实施例中，λ＝2.0。(6)式的优化问题可以被分解成4个子问题，我们将分别进行求解。

步骤3.2：求解F子问题：

其中下标t表示迭代次数。上述问题相当于:

I是单位矩阵。方程(2)可以用许多有效的数值方法求解，如高斯消去法。

步骤3.3：求解R子问题：

等式(9)类似于矢量全变分问题[4]，可以通过快速迭代收缩阈值算法来解决[5]。

步骤3.4：求解S子问题：

等式(10)是低轶优化问题，可以通过奇异值阈值算法解决[6]。

步骤3.5：求解Y子问题。该子问题可由下式直接更新：

Y_t+1＝Y_t+λ(F_t-R_t-S_t) (11)

优化求解的流程如表1所示。

表1

步骤4：基于上述步骤得到高分辨率多光谱图像，为了对融合图像进行定量评价，我们使用了PCA[7]，GSA[8]，NIHS[9]，Indusion[10]，HPF[10]，DLP[12]和LGC[13]六种方法进行对比，使用了ERGSA[14]，RASE[15]，RMSE，UIQI[16]和SAM[17]和五种质量评价指标。定量分析结果如表2所示，粗体表示最佳结果。

表2不同方法的融合图像的定量指标对比

可以看出，由于我们的融合方法采用了二阶导来捕获更高阶的边缘，因此可以更精确地保留空间信息并同时减少因空间信息不准确而导致的光谱失真，因此我们的融合图像的各项指标均最接近于理想值。

本发明主要针对由低分辨率多光谱图像和高分辨率全色图像融合得到高分辨率多光谱图像的应用需求。考虑到，我们的模型将图像的cartoon+texture分解理论与图像全锐化相结合，以有效的确保PAN图像和pan-sharpened MS图像的空间一致性。此外，我们的模型中基于加权核范数的低轶的纹理分量先验项和基于一阶梯度的结构分量先验项能够进一步保留图像空间结构和光谱信息。最后，所提出的模型使用ADMM算法求解，从而获得较小空间失真及光谱失真的高分辨率多光谱图像。

应当理解的是，本说明未详细阐述的部分均属于现有技术。

应当理解的是，上述针对实施例的描述较为详细，并不能因此而认为是对本发明专利保护范围的限制，本领域的普通技术人员在本发明的启示下，在不脱离本发明权利要求所保护的范围情况下，还可以做出替换或变形，均落入本发明的保护范围之内，本发明的请求保护范围应以所附权利要求为准。

参考文献

[1]Lucien Wald,Thierry Ranchin,and Marc Mangolini,“Fusion ofsatellite images of different spatial resolutions:Assessing the quality ofresulting images,”Photogrammetric Engineering and Remote Sensing,vol.63,no.6,pp.691–699,1997.

[2]S.Ono,T.Miyata,and I.Yamada,“Cartoon-texture image decompositionusing blockwise low-rank texture characterization,”IEEE Transactions on ImageProcessing,vol.23,no.3,pp.1128–1142,March 2014.

[3]C.Chen,Y.Li,W.Liu,and J.Huang,“Sirf:Simultaneous satellite imageregistration and fusion in a unified framework,”IEEE Transactions on ImageProcessing,vol.24,no.11,pp.4213–4224,Nov 2015.

[4]Xavier Bresson,“Fast dual minimization of the vectorial totalvariation norm and applications to color image processing,”2008.

[5]A.Beck and M.Teboulle,“A fast iterative shrinkagethresholdingalgorithm for linear inverse problems,”SIAM Journal on Imaging Sciences,vol.2,no.1,pp.183–202,2009.

[6]Jian-Feng Cai,Emmanuel Candes,and Zuowei Shen,“A`singular valuethresholding algorithm for matrix completion,”SIAM Journal on Optimization,vol.20,pp.1956–1982,03 2010.

[7]Zhijun Wang,D.Ziou,C.Armenakis,D.Li,and Qingquan Li,“A comparativeanalysis of image fusion methods,”IEEE Transactions on Geoscience and RemoteSensing,vol.43,no.6,pp.1391–1402,June 2005.

[8]B.Aiazzi,S.Baronti,and M.Selva,“Improving component substitutionpansharpening through multivariate regression of ms+pan data,”IEEETransactions on Geoscience and Remote Sensing,vol.45,no.10,pp.3230–3239,Oct2007.

[9]M.Ghahremani and H.Ghassemian,“Nonlinear ihs:A promising methodfor pan-sharpening,”IEEE Geoscience and Remote Sensing Letters,vol.13,no.11,pp.1606–1610,Nov 2016.

[10]M.M.Khan,J.Chanussot,L.Condat,and A.Montanvert,“Indusion:Fusionof multispectral and panchromatic images using the induction scalingtechnique,”IEEE Geoscience and Remote Sensing Letters,vol.5,no.1,pp.98–102,Jan 2008.

[11]PSJr Chavez,SC Sides,and JA Anderson,“Comparison of threedifferent methods to merge multiresolution and multispectral data:Landsat tmand spot panchromatic,”Photogrammetric Engineering and Remote Sensing,vol.57,no.3,pp.295—303,March 1991.

[12]Giuseppe Masi,Davide Cozzolino,Luisa Verdoliva,and GiuseppeScarpa,“Pansharpening by convolutional neural networks,”Remote Sensing,vol.8,no.7,pp.594,2016.

[13]Xueyang Fu,Zihuang Lin,Yue Huang,and Xinghao Ding,“A variationalpan-sharpening with local gradient constraints,”in 2019 IEEE Conference onComputer Vision and Pattern Recognition(CVPR),2019.

[14]Israa Amro,Javier Mateos,Miguel Vega,Rafael Molina,and Aggelos KKatsaggelos,“A survey of classical methods and new trends in pansharpening ofmultispectral images,”Eurasip Journal on Advances in Signal Processing,vol.2011,no.1,pp.79,2011.

[15]Myungjin Choi,“A new intensity-hue-saturation fusion approach toimage fusion with a tradeoff parameter,”IEEE Transactions on Geoscience andRemote Sensing,vol.44,no.6,pp.1672–1682,June 2006.

[16]Zhou Wang and A.C.Bovik,“A universal image quality index,”IEEESignal Processing Letters,vol.9,no.3,pp.81–84,2002-03.

[17]X.Kang,S.Li,and J.A.Benediktsson,“Pansharpening with mattingmodel,”IEEE Transactions on Geoscience and Remote Sensing,vol.52,no.8,pp.5088–5099,Aug 2014.

Claims (2)

1.一种基于结构-纹理分解的遥感图像变分融合方法，其特征在于，包括以下步骤：

步骤1，选取同一时间同一画面的经过配准处理的多光谱图像和全色图像，低分辨率多光谱图像的记为M，全色图像记为P，融合得到的高分辨率多光谱图像记为F，全色图像与多光谱图像之间的分辨率之比为c；

步骤2，构建数据保真项J(F)来保留低分辨率多光谱图像M的光谱信息，构建先验约束项K(F)以确保融合图像F与全色图像P之间的空间结构一致性，从而获得目标能量函数：

arg min_FJ(F)+K(F) (1)

步骤2的具体实现包括以下子步骤，

步骤2.1，构建数据保真项J(F)，来保留低分辨率多光谱图像M的光谱信息，低分辨率多光谱图像M被视为高分辨率多光谱图像F的下采样模糊版本，即M＝ΨF+n，其中Ψ表示下采样运算符，n表示均值为0的高斯噪声；因此，为了保留LR MS图像的光谱信息，构造保真项为：

其中

表示Frobenius范数；

步骤2.2，基于全色图像和相应的高分辨率多光谱图像之间结构分量和纹理分量的相似性，在结构-纹理空间构建先验约束项K(F)；具体原理如下：

融合图像F表示为纹理分量R和结构分量S的组合：

F＝R+S (3)

其中结构分量S代表的是图像中的全局结构信息，即大面积平滑区域和锐利的边缘，纹理分量R代表的是图像中的小尺寸细节内容；

选择使用核范数来对纹理分量进行约束，因此，先验约束项K(F)：

其中，

和

表示将U和V分别复制到b个波段，U和V为全色图像P的结构分量和纹理分量，由图像分解算法获得；R和S分别为F的结构和纹理分量，

为梯度算子，‖·‖_2,1表示群体稀疏；‖·‖_*表示核范数，ρ和τ为平衡各项的参数；

为块映射函数，用于将矩阵

无重叠的分割为r×r的块并将每个块重新排列成长度为r²的列向量；

步骤2.3：基于子步骤2.1和2.2，提出的目标能量函数模型为：

步骤3：基于步骤2的目标能量函数，通过交替方向乘子法来求解优化问题，从而得到融合图像。

2.如权利要求1所述的一种基于结构-纹理分解的遥感图像变分融合方法，其特征在于：步骤3的具体实现包括以下子步骤，

步骤3.1：对于步骤2得到的目标能量函数模型，通过交替方向乘子法来求解该优化问题，通过引入辅助变量Y，(5)的增广拉格朗日函数表示为：

其中λ为拉格朗日参数，(6)式的优化问题被分解成4个子问题，分别进行求解；

步骤3.2：求解F子问题：

其中下标t表示迭代次数，上述问题相当于：

其中I是单位矩阵，等式(8)通过高斯消去法求解,Ψ′表示Ψ的转置；

步骤3.3：求解R子问题：

等式(9)通过快速迭代收缩阈值算法来解决；

步骤3.4：求解S子问题：

等式(10)是低轶优化问题，通过奇异值阈值算法解决；

步骤3.5：求解Y子问题，该子问题由下式直接更新：

Y_t+1＝Y_t+λ(F_t-R_t-S_t) (11)。

Images