CN102609930B - 一种基于多方向梯度场的图像融合方法 - Google Patents

Abstract

本发明公开了一种高质量快速融合多光谱图像与全色图像的方法。该方法通过多方向梯度算子保持全色图像的细节信息，并利用下采样算子保持多光谱的颜色信息，从而得到高分辨率的多光谱图像。此外，本发明通过一等价转换将复杂的空间域计算转换到频域上，不仅能减少计算机的内存占用，而且可以减小运算复杂度，从而能够实时计算。本发明解决了现有图像融合方法无法同时保持全色图像的细节信息和多光谱图像的细节信息的难题，通过本发明的方法，能够使融合后的图像在空间分辨率和光谱分辨率上同时达到最优。

Description

一种基于多方向梯度场的图像融合方法

技术领域

本发明属于遥感图像处理领域，特别涉及一种基于多方向梯度场的全色图像与多光谱图像的融合方法。

背景技术

在遥感图像处理领域，对于光学传感器，图像的空间分辨率和光谱分辨率是一对矛盾的因素。要获得高空间分辨率的图像就只能以单光谱工作，而要获得更多光谱就必须降低空间分辨率。为了同时提高图像的空间分辨率和光谱分辨率，人们提出了利用融合全色图像(高空间分辨率单光谱图像)和多光谱图像(低空间分辨率)来得到高空间分辨率多光谱图像的方法。总结已有的图像融合方法，可以将其分为三类：

第一类是基于颜色空间的转换方法。这些方法先将多光谱从原颜色空间转换到另外一种颜色空间，这样可以将亮度信息分量和颜色信息分量分开，再将亮度信息分量用单光谱的全色图像替换得到替换后的颜色空间表示，最后再对新的颜色空间进行逆变换，就得到融合后的图像。如，IHS(Intensity，Hue，Saturation)变换法先将原多光谱信息看作是RGB空间的颜色图像，对它进行IHS变换分别得到亮度通道I、色调通道H和饱和度通道S。之后用全色图像替代亮度通道I，再进行IHS逆变换从IHS空间变换到RGB空间，从而得到融合后的图像。基于颜色空间转换的方法实现原理简单，这种基于颜色空间转换的方法存在着两个方面的问题：(1)只适合于三通道图像融合；(2)虽然用全色图像直接替代了亮度通道，但是色调通道和饱和通道仍是原多光谱的信息的简单上采样，这必然会导致空间细节上的损失。

第二类是基于统计的主成份方法，这类方法与颜色空间转换类似，不同的是它可以适用于任何波段的图像融合，这类方法先将多光谱的每个像素看成是一个多维向量，对其进行统计成份分析(如主成分分析PCA)，得到投影向量，再将多光谱图像在这些投影向量上进行投影，在主成份投影的值保留了图像大部分的信息。再用全色图像替换主成分图像，之后进行相应的反变换就可以得到融合图像。基于统计的主成份方法虽然可以进行任意波段图像的融合，但融合后的图像很难保持原多光谱图像的颜色信息(即发生颜色失真)，而且在空间细节上与原全色图像相比也存在一定的损失。

第三类方法是基于小波分解法，这类方法的基本思想是利用小波变换提取全色图像的高频信息(细节信息)和多光谱图像的低频信息(轮廓信息)。将通过小波变换提取的全色图像的高频信息和多光谱图像的低频信息组合成一组新的小波系数，再对这些小波系数进行小波逆变换，从而得到融合图像。但是基于小波变换的图像融合方法得到的融合图像容易产生人造痕迹，从而降低了融合图像的质量。

综上所述，虽然现有的全色图像与多光谱图像融合方法很多，但这些方法都不能在空间分辨率和光谱分辨率上同时达到最优。而对于遥感图像的应用而言，得到高质量的高分辨率多光谱图像至关重要。如何由已有的全色图像和多光谱图像得到高质量的高分辨率多光谱图像在遥感图像处理领域中仍是一个颇具挑战性的难题。

发明内容

本发明的目的旨在使融合后的图像尽可能的保持全色图像的细节信息和多光谱图像的颜色信息，为此本发明提供了一种高质量的图像融合方法。

为达成所述目的，本发明提供一种基于多方向梯度场的图像融合方法，其特征在于，该方法包括以下步骤：

步骤S1，为了使融合图像保持全色图像的细节信息，基于多方向梯度算子，根据融合图像的多方向梯度场与全色图像的多方向梯度场的差值，得到融合图像在多方向梯度场上的目标函数L₁(f)，其中，f为融合图像；

步骤S2，为了使融合图像保持多光谱图像的颜色信息，根据融合图像的下采样图像与多光谱图像之间的差别，得到融合图像在颜色信息上的目标函数L₂(f)；

步骤S3，基于所述目标函数L₁(f)和所述目标函数L₂(f)，得到融合图像的整体目标函数L(f)；

步骤S4，在频域最小化所述整体目标函数L(f)，得到融合图像的频谱

步骤S5，对频谱

进行FFT逆变换，得到融合图像。

本发明的所述方法针对全色图像和多光谱图像融合问题，可实现空间分辨率和光谱分辨率的最优融合，得到高质量的高空间分辨率多光谱图像。其优点体现在：1)该方法利用多方向梯度算子，使融合后的图像尽可能地保持全色的细节信息；2)该方法利用下采样算子使融合后的图像下采样后与多光谱图像保持一致，这最大可能的保持了多光谱图像的颜色信息；3)通过对下采样目标函数进行近似转换后，图像融合的目标函数可以转换到频率域求解，避免在空间域上求解超大规模的稀疏线性方程组问题，这不仅解决了内存占用问题，同时降低了运算复杂度，即使对于大规模遥感图像的融合，不仅能得到高质量的融合结果同时也可达到实时计算。该方法对于遥感图像处理有着重要的应用价值。

附图说明

图1是本发明提出的基于多方向梯度场的图像融合方法的总体流程图。

图2是融合图像保持高分辨全色图像细节信息的示意图。

图3是融合图像保持多光谱图像颜色信息的示意图。

图4是经过近似转换后保持多光谱图像颜色信息的示意图。

图5是根据本发明的2米多光谱图像与1米全色图像融合的效果图。

图6是根据本发明的2米多光谱图像与1米全色图像融合的另一效果图。

图7是根据本发明的2米多光谱图像与0.5米全色图像融合的效果图

图8是根据本发明的2米多光谱图像与0.5米全色图像融合的另一效果图。

具体实施方式

为使本发明的目的、技术方案和优点更加清楚明白，以下结合具体实施例，并参照附图，对本发明进一步详细说明。

本发明的输入图像为配准好的高分辨率全色图像和低分辨率多光谱图像，输出为高分辨率多光谱图像。

如图1所示，本发明提供了一种基于多方向梯度场的图像融合方法，该方法具体包括以下步骤：

步骤S1，为了使融合图像保持全色图像的细节信息，基于多方向梯度算子，根据融合图像的多方向梯度场与全色图像的多方向梯度场的差值，得到融合图像在多方向梯度场上的目标函数；

利用多方向梯度场的目的是为了保持全色图像的细节信息，为此，本发明定义了4个多方向梯度算子：

$g_1=\left(\begin{array}{ccc}-1& 0& 1\\ -2& 0& 2\\ -1& 0& 1\end{array}\right),$ $g_2=\left(\begin{array}{ccc}-1& -2& -1\\ 0& 0& 0\\ 1& 2& 1\end{array}\right),$

$g_3=\left(\begin{array}{ccc}-2& -1& 0\\ -1& 0& 1\\ 0& 1& 2\end{array}\right),$ $g_4=\left(\begin{array}{ccc}0& -1& -2\\ 1& 0& -1\\ 2& 1& 0\end{array}\right),$

其中，g₁，g₂分别为x，y方向的梯度算子，g₃，g₄分别为xy，yx方向的梯度算子，利用多方向的梯度算子提取全色图像的多方向梯度场，就可以尽可能的提取全色图像的细节信息。

本发明的融合目标之一是使融合后的图像保持全色图像的多方向梯度场，这就意味着融合后的图像能尽可能的保持全色图像的细节信息。图2给出了多方向梯度场保持的示意图。因此，根据融合图像的多方向梯度场与全色图像的多方向梯度场的差值，得到融合图像在多方向梯度场上的目标函数为：

$L_1\left(f\right)=\underset{i=1}{\overset{4}{\mathrm{\Σ}}}{|f\⊗g_i-I\⊗g_i|}^2,$

其中，f为融合图像，I为全色图像，

为卷积操作。

最小化上面的目标函数就可以使融合后的图像尽可能保持全色图像的多方向梯度场，利用上式对f求导可得：

$\frac{{\mathrm{dL}}_1}{\mathrm{df}}=\underset{i=1}{\overset{4}{\mathrm{\Σ}}}((f\⊗g_i-I\⊗g_i)\⊗{\stackrel{\‾}{g}}_i),$

其中，

表示为将梯度算子g_i上下、左右翻转后的翻转梯度算子。

上式中的卷积操作可以转换为频率域的乘积，对上式等号的右边进行FFT变换可得：

其中，

分别为f、g_i、I、

的傅里叶变换。

步骤S2，为了使融合图像保持多光谱图像的颜色信息，根据融合图像的下采样图像与多光谱图像之间的差别，得到融合图像在颜色信息上的目标函数L₂(f)；

本发明的第二个目标是保持多光谱图像的光谱信息：即融合图像的下采样图像应与多光谱图像相等，图3给出了这个目标的示意图。为了保持多光谱图像的颜色信息，根据融合图像的下采样图像与多光谱图像之间的差别，得到融合图像在颜色信息上的目标函数：

L₂(f)＝|f↓-C|²，

其中，↓表示下采样算子，f↓表示融合图像f的下采样图像，C为多光谱图像。

除了保持全色图像的细节信息，本发明的融合方法的另外一个目标是最小化上述目标函数来保持多光谱图像的颜色信息。要最小化上述目标函数，需要对上式中的f进行求导，但是上式中存在的下采样算子使得对f的求导不是一个显然的过程，实际上由于下采样的操作一般是对图像进行高斯平滑，然后用最近邻采样得到下采样图像，这些操作都是线性操作，那么如果将融合图像f按列拉成列向量，那么下采样操作就可以写成矩阵-向量相乘的形式即：

其中，

表示融合图像f按列拉成的列向量，D为一个大规模的稀疏矩阵，其作用等价于下采样。

比如，如果全色图像I的空间分辨率为M×N，多光谱图像C的空间分辨率为m×n，则图像融合的目是得到空间分辨率为M×N的多光谱融合图像f。根据上述描述，目标函数L₂(f)中对图像f的下采样操作f↓可写为矩阵-向量相乘的形式D

此时矩阵D为一个mn×MN的稀疏矩阵，D中每一行的非零元素由所需采样图像f中的对应像素及采样算法确定。具体而言，当多光谱图像C的空间分辨率在x，y方向上的分辨率为全色图像I的1/2，即m×n等于M/2×N/2时，如果采用双线性插值算法来进行下采样，由于多光谱图像C的分辨率在x，y方向上的分辨率为全色图像I的1/2，那么下采样图像的每个像素由图像f中所对应的四个相邻像素双线性插值得到，这四个像素的加权系数为1/4。此时稀疏矩阵D为如下形式：

$D=\left[\begin{array}{cccccccccc}\frac{1}{4}& \frac{1}{4}& 0& \·\·\·& \frac{1}{4}& \frac{1}{4}& 0& \·\·\·& 0& 0\\ 0& \frac{1}{4}& \frac{1}{4}& \·\·\·& 0& \frac{1}{4}& \frac{1}{4}& \·\·\·& 0& 0\\ & & & & & \·& & & & \\ & & & & & \·& & & & \\ & & & & & \·& & & & \\ & \·\·\·& & \·\·\·& & \·\·\·& & \·\·\·& \frac{1}{4}& \frac{1}{4}\end{array}\right],$

此时D为一个MN/4×MN的稀疏矩阵。虽然矩阵D是一个稀疏矩阵，但是对于遥感图像来说，由于图像的空间分辨率都比较大，即使采用这种稀疏矩阵的表示方式，对计算机的内存来说仍然是一个很大的需求，而且求解的计算复杂度也很大，不适用于高分辨率遥感图像的处理。

为了解决上述问题，我们将融合目标进行一个近似的转换，如图4所示，原来的目标是使融合后的图像其下采样图像等于多光谱图像，现在将目标转化为使融合后的图像经过模糊后等于多光谱图像的上采样图像。所以，为了保持多光谱图像颜色信息，定义如下转换后的目标函数：

$L_2^{\′}\left(f\right)={|f\⊗k-\tilde{C}|}^2,$

其中，k为高斯模糊核，为多光谱图像C采用双三次插值的上采样图像。

最小化上式的目标函数就可以使融合后的图像尽可能地保持多光谱图像C的颜色信息，利用上式对f求导可得：

$\frac{{\mathrm{dL}}_2^{\′}}{\mathrm{df}}=(f\⊗k-\tilde{C})\⊗\stackrel{\‾}{k},$

其中，

表示为将高斯模糊核k上下、左右翻转后的翻转模糊算子。

经过上述近似转换后，上式中的卷积操作同样也可以转换到频率域上来计算，对上式等号的右边进行FFT变换可得到：

其中，

分别为f、k、

的傅里叶变换。

步骤S3，基于所述目标函数L₁(f)和所述目标函数L′₂(f)，得到融合图像的整体目标函数L(f)；

本发明提出的图像融合方法的目的是在空间分辨率和光谱分辨率上同时达到最优。为了保持全色图像的细节信息，利用多方向的梯度算子来使融合图像的多方向梯度场等于全色图像的多方向梯度场，其目标优化可以转换到频率域上计算；为了保持多光谱图像的颜色信息，利用下采样操作使融合图像的下采样图像等于多光谱图像，并对这一目标进行近似变换即：融合图像经过模糊后与多光谱图像的上采样图像相等。本发明的图像融合方法就是同时优化这两个目标以到达空间分辨率和光谱分辨率的同时最优。

综上，基于所述目标函数L₁(f)和所述目标函数L′₂(f)，本发明所提出的图像融合的整体目标函数可以定义为：

$L\left(f\right)=L_1\left(f\right)+{\mathrm{\λL}}_2^{\′}\left(f\right)$

$=\underset{i=1}{\overset{4}{\mathrm{\Σ}}}{|f\⊗g_i-I\⊗g_i|}^2+\mathrm{\λ}{|f\⊗k-\tilde{C}|}^2,$

其中，λ为权值系数。由于保持全色图像的空间细节和保持多光谱图像的颜色信息同等重要，一般情况下λ的取值为1，实际应用中可根据需求调整λ的取值，建议λ的取值范围为[0.8，1.2]。

步骤S4，在频域最小化所述整体目标函数L(f)，得到融合图像的频谱

最小化上述整体目标函数L(f)就可以使融合图像既保持全色图像的细节，同时又能保持多光谱图像的细节。

步骤S4进一步包括以下步骤：

步骤S41，为了最小化L(f)，对f进行求导并令其为0，得到：

$\begin{array}{c}\frac{\mathrm{dL}}{\mathrm{df}}=\frac{{\mathrm{dL}}_1}{\mathrm{df}}+\mathrm{\λ}\frac{{\mathrm{dL}}_2^{\′}}{\mathrm{df}}=0\\ \⇓\end{array},$

$\underset{i=1}{\overset{4}{\mathrm{\Σ}}}((f\⊗g_i-I\⊗g_i)\⊗{\stackrel{\‾}{g}}_i)+\mathrm{\λ}(f\⊗k-\tilde{C})\⊗\stackrel{\‾}{k}=0.$

步骤S42，将上式转换到频域求解，得到融合图像f的频谱

如果在空间域直接计算上式，需要将卷积先转换为矩阵-向量乘积的形式，这样不仅会占用大量的计算机内存，而且目标函数的优化也会很慢，不利用实际的应用。因此，本发明将上式转换到频率域求解，得到融合图像f的频谱

为：

从上式可以看出计算

的最耗时部分是FFT变换，而FFT变换的复杂度为

(NlogN)，这大大加快了求解速度。

步骤S5，对频谱进行FFT逆变换，得到融合后的图像。

由上面的求解得到了融合图像的频谱

，再对频谱

进行傅里叶反变换就可以得到融合图像f：

其中，ifft(·)表示傅里叶反变换。

图像融合结果：为验证本发明的技术，进行一系列真实多光谱图像和全色图像的融合实验。图5至图8所示是真实全色图像和多光谱图像融合效果图。图5和图6为2米多光谱图像与1米全色图像的融合实验效果图。图5中从左往右分别为2米多光谱图像、1米全色图像和融合图像，可以看出利用本发明的图像融合方法不仅能很好的保持全色图像的细节，同时也保持了多光谱图像的细节。图6中从上到下分别为2米多光谱图像、0.5米全色图像和融合图像，从图中可以看出融合后的图像质量既具有全色图像的清晰度同时又具有多光谱图像的颜色。图7和图8为2米多光谱与0.5米全色图像的融合实验效果图。图7中从左往右分别为2米多光谱图像、0.5米全色图像和融合图像，从中可以看出在边缘细节上融合后的图像非常好的保持了全色图像的细节(如实线框区域)，同时在颜色上也保持了多光谱的颜色信息(如虚线框区域)。图8中从左往右分别为2米多光谱图像、0.5米全色图像和融合图像。从融合后得到的图像质量有效地说明了本发明所述方法能得到高质量的高分辨率多光谱图像。上述两个实验的图像分辨率都在1000×1000左右，在一台2.4Ghz CPU/2G内存的机器上所需的运行时间不超过0.5秒。

以上所述的具体实施例，对本发明的目的、技术方案和有益效果进行了进一步详细说明，所应理解的是，以上所述仅为本发明的具体实施例而已，并不用于限制本发明，凡在本发明的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (5)

1.一种基于多方向梯度场的图像融合方法，其特征在于，该方法包括以下步骤：

步骤S1，为了使融合图像保持全色图像的细节信息，所述全色图像是高空间分辨率的单光谱图像，基于多方向梯度算子，根据融合图像的多方向梯度场与全色图像的多方向梯度场的差值，得到融合图像在多方向梯度场上的目标函数L₁(f)，所述目标函数L₁(f)为：

$L_1\left(f\right)=\underset{i=1}{\overset{4}{\mathrm{\Σ}}}{|f\⊗g_i-I\⊗g_i|}^2,$

其中，f为融合图像，I为全色图像，g_i为多方向梯度算子，

为卷积操作，所述多方向梯度算子分别为：

$g_1=\left(\begin{array}{ccc}-1& 0& 1\\ -2& 0& 2\\ -1& 0& 1\end{array}\right),g_2=\left(\begin{array}{ccc}-1& -2& -1\\ 0& 0& 0\\ 1& 2& 1\end{array}\right),$

$g_3=\left(\begin{array}{ccc}-2& -1& 0\\ -1& 0& 1\\ 0& 1& 2\end{array}\right),g_4=\left(\begin{array}{ccc}0& -1& -2\\ 1& 0& -1\\ 2& 1& 0\end{array}\right);$

步骤S2，为了使融合图像保持多光谱图像的颜色信息，所述多光谱图像是低空间分辨率的多光谱图像，根据融合图像的下采样图像与多光谱图像之间的差别，得到融合图像在颜色信息上的目标函数L₂(f)，所述目标函数L₂(f)为：

L₂(f)=|f↓-C|²，

其中，↓表示下采样算子，f↓表示融合图像f的下采样图像，C为多光谱图像；

步骤S3，基于所述目标函数L₁(f)和所述目标函数L₂(f)，得到图像融合的整体目标函数L(f)，所述整体目标函数L(f)为：

L(f)=L₁(f)+λL′₂(f)，

其中，λ为权值系数，k为高斯模糊核，

为多光谱图像C的上采样图像；

步骤S4，在频域最小化所述整体目标函数L(f)，得到融合图像的频谱

步骤S5，对频谱

进行FFT逆变换，得到融合图像。

2.根据权利要求1所述的方法，其特征在于，所述步骤S4进一步包括以下步骤：

步骤S41，对L(f)中的f进行求导并令其为0，得到：

$\frac{\mathrm{dL}}{\mathrm{df}}=\frac{{\mathrm{dL}}_1}{\mathrm{df}}+\mathrm{\λ}\frac{{\mathrm{dL}}_2^{\′}}{\mathrm{df}}=0;$

步骤S42，将上式转换到频域求解，得到融合图像f的频谱

3.根据权利要求2所述的方法，其特征在于，所述步骤S41中：

$\frac{{\mathrm{dL}}_1}{\mathrm{df}}=\underset{i=1}{\overset{4}{\mathrm{\Σ}}}((f\⊗g_i-I\⊗g_i)\⊗{\stackrel{\‾}{g}}_i),\frac{{\mathrm{dL}}_2^{\′}}{\mathrm{df}}=(f\⊗k-\tilde{C})\⊗\stackrel{\‾}{k},$

其中，为将梯度算子g_i上下、左右翻转后得到的翻转梯度算子；

为将高斯模糊核k上下、左右翻转后得到的翻转模糊算子。

4.根据权利要求3所述的方法，其特征在于，所述步骤S42中

的频域表示为：

其中，

分别为f、g_i、I、

的傅里叶变换；

的频域表示为：

其中，

分别为k、的傅里叶变换。

5.根据权利要求4所述的方法，其特征在于，所述融合图像f的频谱

为：

Images