基于低频相关性分析的全色图像与多光谱图像融合方法

摘要

本发明是一种基于低频相关性分析的全色图像与多光谱图像融合方法，利用匹配更新后的全色图像的梯度场作为图像空间几何信息表达，形成对高分辨率融合图像的空间几何信息约束；同时引入光谱信息约束，使高分辨率融合图像的低频部分与原多光谱图像吻合；通过对低频全色图像与多光谱图像的相关性分析，自动确定两种约束之间的比重关系，从而保证空间几何信息与光谱信息在不同的图像区域合理融合；最后，在两种约束间加入一个可调参数，可对融合结果中的几何与光谱信息相对含量进行调节。本发明的融合效果优于当前流行的IHS，PCA，DWT等算法，并且几何信息与光谱信息相对含量可自由调节，提高了本方法的实用性。

说明书

技术领域

本发明涉及一种基于低频相关性分析的全色与多光谱图像融合方法，属于图像融合领域。

背景技术

遥感影像数据种类越来越多，不同传感器、不同波段与不同时相数据如何进行优势互补和有机结合是广泛关注的问题。遥感图像融合很好地解决了这个问题，它是指不同图像在经过空间配准之后，按照一定的原则与算法进行整合，结合互补信息形成新的影像数据。全色图像具有较高的空间分辨率，多光谱图像保存着不同波段的光谱信息。大量研究工作表明，把全色图像和多光谱图像融合起来，在保持原光谱信息的前提下提高多光谱图像的空间分辨率，有利于提高多光谱遥感影像的信息分析和提取能力。

多数学者认为现有全色与多光谱图像融合算法主要分为两大类：分量替换法(Component-substitution Method)和多尺度分析法(Multi-resolution Analysis)。如IHS、PCA、Gram-Schmidt变换法等属于分量替换法，它们的共同特点是：利用某种变换方法将多光谱影像变换成不同分量，然后将全色图像替代某一分量，再进行反变换得到融合结果。由于全色图像具有较高的空间分辨率，反变换得到的多光谱融合图像也具有相同的空间分辨率。如，IHS算法是将多光谱图像进行IHS变换，得到I、H和S三个分量，然后用高分辨率全色图像代替I分量，最后进行IHS反变换得到融合图像。但是这种简单的“变换-替换-反变换”方法容易破坏源多光谱图像的光谱信息，特别是当全色图像波段范围与多光谱图像波段范围相差较大时尤为严重。

多尺度分析法是指对待融合的影像数据进行多尺度分解，然后再重构的融合技术，如离散小波变换法(DWT)、高通滤波法(HPF)等。对于多光谱图像，一般认为低频部分保存着图像光谱信息，高频部分体现的是图像的纹理细节信息。基于这种原理，多尺度分析法的目的是融合多光谱波段低频光谱信息和全色波段高频几何信息，从而生成高空间分辨率的多光谱图像。多尺度分析法的研究主要集中在两方面：多尺度分解方法和高低频信息融合策略。经常采用小波变换进行多尺度分解，要求其具有平移不变性，以避免微小的平移变化造成小波系数突变。常用的是具有近似平移不变性，较为容易实现的“átrous”小波(ATW)。其他的分解方法还有：广义拉普拉斯金字塔(GLP)、Curvelet变换、Contourlet变换等。

多尺度分析法往往在保持光谱信息上表现较好，但是在增强空间分辨率上效果不佳，主要原因在于没有找到合适的高低频信息融合策略。国外学者较早地提出了高低频融合的ARSIS思想，即通过有选择性地注入高频几何结构信息来增强多光谱图像的空间分辨率。这种“高频信息注入”思想的目的是，使多光谱图像合理地接受和融合全色波段的高频信息，从而使高空间分辨率融合图像在物理上具有较高的真实性，相当于高空间分辨率多光谱成像传感器的实际拍摄效果。基于这种思想，提出了一些高频信息注入模型和方法，如RWM、SDM、CBD模型及其各种改进版本，以及基于MMSE、遗传算法、Kalman滤波的模型参数估计和优化方法。它们的基本思想是，在低尺度上构造高频信息注入模型并估计模型参数，然后利用这些估计的模型在高频段执行信息注入。这些模型在一定程度上提升了多尺度分析法的融合效果。

尽管如此，多尺度分析法还是很难在“保持光谱信息”和“合理提高空间分辨率”上都达到理想效果。近些年来，基于分量替换和多尺度分析优缺点互补的原则也提出了不少混合型的融合算法，如DWT-IHS、FFT-IHS、DWT-PCA等。但是，融合过程中缺乏对数据的详细分析，在很大程度上限制了这些方法的融合性能。

发明内容

现有全色图像与多光谱图像融合技术难以在“保持光谱信息”和“合理提高空间分辨率上”都达到很好的效果。为此，本发明基于低频相关性分析的全色图像与多光谱图像融合方法的目的是通过利用低频全色图像与多光谱图像的相关性分析，使全色图像的空间几何信息与多光谱图像的光谱信息合理融合，在保持光谱信息的基础上最大限度地提高空间分辨率。

为实现所述的目的，本发明提供一种基于低频相关性分析的全色图像与多光谱图像融合方法，所采用的技术方案包含以下步骤：

步骤S1：利用低通滤波器对输入的全色图像P进行低通滤波，得到低频全色图像其中k为滤波核；

步骤S2：使用计算机中的相关性分析模块对低频全色图像和输入的多光谱图像进行相关性分析，计算并得到低频全色图像与多光谱图像每一区域处的相关系数；

步骤S3：根据多光谱图像与低频全色图像的亮度比值关系，对全色图像进行匹配更新得到新全色图像，并对匹配更新后的新全色图像进行计算，得到新全色图像的梯度场；

步骤S4：以新全色图像梯度场作为待求融合图像的空间几何信息约束，同时使所述融合图像的低频部分充分接近于多光谱图像，得到光谱信息约束，利用这两种约束并结合使用相关性分析模块得到的相关系数构造融合能量泛函E(I)；

步骤S5：使用梯度下降法最小化融合能量泛函，求解能量泛函最小化时的融合图像I。

本发明方法与现有融合方法相比较，具有以下突出特点和优点：

(1)以匹配更新后的全色图像梯度场来表示空间几何信息，具有简单、准确、直接的优点，便于在融合过程中引入空间几何信息约束。

(2)充分考虑和分析了融合数据的特点。对低频全色图像与多光谱图像数据进行了相关性分析，根据不同的相关性来调整空间几何信息与光谱信息的约束关系，能够在保持光谱信息的前提下，最大限度地提高空间分辨率。

(3)保持光谱信息和提高空间分辨率是全色与多光谱图像融合的两大目标，本发明方法在融合能量泛函中显式地表达出了这两个融合目标，并给出了解决目标间矛盾关系的方法，最后通过引入一个总的调节参数，调整融合结果中空间几何信息和光谱信息的相对含量，可以满足不同的融合需求，便于人机交互，具有很高的实用价值。

附图说明

图1为本发明方法的融合流程图

图2(a)为输入的全色图像

图2(b)为输入的多光谱图像

图2(c)为IHS方法融合图像

图2(d)为DWT方法融合图像

图2(e)为本发明方法融合图像

图2(f)为实际的高分辨率多光谱图像

图3为不同融合方法的客观评价指标比较

具体实施方式

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

本发明首先在融合时引入了空间几何信息约束和光谱信息约束；以匹配更新后的全色图像梯度场表示图像的空间几何信息，使融合图像具有类似的几何结构；基于“图像低频部分体现光谱信息”这一思想，引入了光谱信息约束，使融合图像的低频部分与原多光谱图像吻合。使用这两种约束构造融合能量泛函，能量泛函值越小表示越接近既定的融合目标。

为了更好地整合两种约束关系，使全色图像的空间几何信息与多光谱图像的光谱信息合理融合，首先对低频全色图像和多光谱图像进行相关性分析，得出每一图像区域的相关系数。相关系数高，表明可以加大该区域空间几何信息约束力度，进一步提高空间分辨率，同时不会引起光谱失真；相关系数低时，表明该区域不宜从全色图像引入过多的几何结构，否则容易引起光谱畸变，此时应相对地加大光谱信息约束力度。为了便于对图像整体的空间几何信息与光谱信息相对含量进行统一调整，加入了一个调节参数，可以满足不同的融合需求，也便于人机交互。

图1是本发明的融合流程图。图2(a)与图2(b)是本实例输入的全色和多光谱图像。全色图像空间分辨率为0.7m，多光谱图像具有3个波段，空间分辨率为2.8m。上述图像已经过精确的配准预处理。

本实例将多光谱图像中的各波段分别与全色图像进行融合处理，利用计算机实现低频全色图像与多光谱图像的相关性分析，具体步骤如下：

(1)利用低通滤波器对输入的全色图像P进行低通滤波，得到低频全色图像其中k为滤波核。

(2)使用计算机中的相关性分析模块对低频全色图像与多光谱图像S进行相关性分析，计算并得到低频全色图像与多光谱图像每一区域处的相关系数ρ为：

$\rho (x,y)=\frac{{\tilde{P}}_1{S}_1+{\tilde{P}}_2{S}_2+···+{\tilde{P}}_n{S}_n}{\sqrt{{{\tilde{P}}_1}^2+{{\tilde{P}}_2}^2+···+{{\tilde{P}}_n}^2}\sqrt{{S_1}^2+{S_2}^2+···+{S_n}^2}},$

式中，表示低频全色图像中，以(x，y)为中心的区域内的像元值，[S₁，S₂，…，S_n]表示在多光谱图像S中，以(x，y)为中心的同样大小区域内的像元值。

(3)根据多光谱图像S与低频全色图像的亮度比值关系，对全色图像P进行匹配更新得到新全色图像，得到全色图像P与多光谱图像S相匹配的新全色图像P′：

$P^{\prime }=\frac{S}{\tilde{P}}·P,$

对新全色图像P′进行计算，得到新全色图像P′的梯度场：

$▿P^{\prime }=\left(\begin{array}{c}\partial P^{\prime }/\partial x\\ \partial P^{\prime }/\partial y\end{array}\right),$

(4)以新全色图像的梯度场作为待求融合图像I的空间几何信息约束，同时使所述融合图像的低频部分k*I与多光谱图像S尽量保持一致，得到光谱信息约束，利用这两种约束并结合得到的相关系数ρ构造能量泛函，所述构造融合能量泛函的步骤包括：

定义相关系数ρ的单调下降函数f(·)，利用相关系数ρ构造函数C(x，y)＝f(ρ(x，y))作为光谱信息约束项权重系数，构造融合能量泛函E(I)如下：

$E\left(I\right)=\underset{\Omega }{\int }\underset{\Omega }{\int }{|▿I-▿P^{\prime }|}^2\mathrm{dxdy}+\beta \underset{\Omega }{\int }\underset{\Omega }{\int }C·{(k*I-S)}^2\mathrm{dxdy},$

式中，C为光谱信息约束项的权重系数，由式C(x，y)＝f(ρ(x，y))定义，其中函数f(·)单调下降；β为可调参数，用来调节融合结果中空间几何信息与光谱信息相对含量，通过利用几何光谱信息调节模块实现对β的调节，为待求融合图像I的梯度场，为新全色图像P′的梯度场。

(5)使用梯度下降法最小化融合能量泛函，求解得到能量泛函最小化时的融合图像I。能量泛函E(I)的一阶变分为：

δE(I)＝2(ΔP′-ΔI)+2β[C(k*I-S)]*k′，

式中，k′(x，y)＝k(-x，-y)。第i步迭代过程为：

I_i+1＝I_i+α·δE，

其中，α为固定迭代步长。取初始迭代图像I₀＝S+(P-k*P)，当两次迭代间的变化足够小时，即停止迭代计算，得到融合结果。

下面对该实例中的其他技术细节问题作详细说明：

(1)所述的相关性分析模块采用的是数字信号相关性分析技术，在计算机中构建专用的相关性分析模块，计算输入信号间的局部相关系数。在对输入的低频全色图像与多光谱图像进行相关性分析时，需要设置相关性分析区域的窗口大小，适宜的窗口大小为3～15个像元，本实例中取为7。

(2)在对全色图像P进行匹配更新中，为防止过大，使新全色图像P′超出图像动态范围，需要对的最大值进行限制。本实例中，其值被限制为不能大于3，即当大于3时，其值被置为3。

(3)在构造融合能量泛函的步骤中，β为可调参数，C为根据相关性分析得到的光谱信息约束的相对权重系数，具体表示为：C(x，y)＝f(ρ(x，y))。其中，f(·)应为单调下降函数，表示随着相关系数ρ(x，y)的提高，空间几何约束相对加强，光谱约束相对减弱。本实例中，取f(·)为线性函数：f(ρ)＝2.01-2ρ，同时令β的取值区间为[0，100]，具体实施时β取值定为80。

(4)为减少求取融合图像时的迭代次数，加快计算速度，取初始迭代图像I₀＝S+(P-k*P)。I₀相当于利用高通滤波融合法(HPF)得到的融合结果，以它作为初始迭代图像可以有效减少迭代次数，提高计算速度。本实例中，将迭代次数设置为30即可。

(5)低通滤波核k大多数情况下可以根据传感器特性给定，也可以通过估计得到。本实例中采用高斯函数近似，对于全色与多光谱图像空间分辨率相差4倍的情况，高斯函数的标准偏差通常可取为2.5。

图2(a)示出输入的全色图像；图2(b)示出输入的多光谱图像；图2(c)、图2(d)分别是IHS，DWT方法的融合结果，图2(e)是本发明方法的融合结果，图2(f)是真实的高分辨率多光谱图像。从目视效果来看，本发明方法的融合结果与图2(f)的实际高分辨率多光谱图像最为接近。图2(c)的融合图像在色彩上发生了较大的改变，表明IHS方法光谱失真比较严重；图2(d)的融合图像光谱扭曲要小些，但在很多局部区域也产生了明显的色彩变化，且不如图2(e)和实际图像清晰。其中，图2(f)中用圆圈框出的屋顶、小车等几个局部区域对比明显。图3列出了与实际图像比较得到的几个客观统计数据，用来评价各算法的融合效果。其中，指标CC表示与真实数据的相关系数，RMSE表示与真实数据的均方根误差。从客观评价指标来看，本发明方法的融合效果也远优于其他方法，在保持光谱信息和提高空间分辨率方面都表现出良好的性能。

以上所述，仅为本发明中的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉该技术的人在本发明所揭露的技术范围内，可理解想到的变换或替换，都应涵盖在本发明的权利要求书的保护范围之内。