CN1254770C - 基于期望值最大与离散小波框架的图像融合方法 - Google Patents

Abstract

一种基于期望值最大与离散小波框架的图像融合方法，首先采用小波框架将待融合图像进行多尺度分解，得到小波框架的一系列高频分量和一个最低频分量，然后对图像的高频部分采用基于视觉特性的特征提取方法进行融合，对图像的低频部分采用基于期望值最大方法进行融合，最后将高频部分的融合结果和低频部分的融合结果，通过小波反变换得到最终的融合图像。本发明使融合后的图像质量得到大幅度地提高，对于各种应用系统的后续处理和图像显示具有重要意义和实用价值。

Description

基于期望值最大与离散小波框架的图像融合方法

技术领域

本发明涉及一种基于期望值最大与离散小波框架图像融合方法，是信息融合领域中一项多尺度和统计方法结合的图像融合方法，在光学成像、目标监视、安全检查等系统中均可有广泛应用。

背景技术

图像融合技术是多传感器信息融合中可视信息的融合，它利用各种成像传感器不同的成像方式，为不同的影像提供互补信息，增加图像信息量，减少原始图像数据量，提高对环境的适应性，以获得更可靠、更准确的有用信息供观察或进一步处理。图像融合技术是一门综合了传感器、信号处理、图像处理和人工智能等的新兴技术，近年来已成为一种十分重要和有用的图像分析与计算机视觉技术，在自动目标识别、计算机视觉、遥感、机器人、医学图像处理以及军事应用等领域有着广泛的应用前景。

目前具有代表性的方法有塔形变换的方法、子波变换等基于多尺度的方法。多尺度图像融合的过程是首先将配准后图像经多尺度分解，分解方法包括拉普拉斯、梯度金字塔以及小波分解等方法；分解后图像的每层看作为图像在这一尺度或波段的特征，根据这些特征所反应的能量测度进行加权平均或选择，以达到融合的目的。这类方法的主要缺点是没有考虑图像的统计特性，从而导致了其对噪声的敏感性。根据国内外文献，目前已有的基于统计理论的图像融合方法只有贝叶斯融合方法和基于期望值最大的方法。现存的基于贝叶斯的融合方法是在假设图像和噪声都服从高斯分布基础上进行的，在实际融合过程中，需要有先验图像的参与，这种方法近似于局部主成分分析方法。现存的期望值最大的方法对上述贝叶斯融合方法的假设条件有所拓宽，在多尺度域假定图像偏移或噪声服从高斯混合模型，应用期望值最大算法估计出融合图像。但在多尺度域的统计特性正是近似服从高斯混合模型，由于将图像高频或通频带部分的信号模拟成噪声，因此容易造成融合结果高频信息的损失。因此不能简单的将这种图像的形成统计模型应用到多尺度域。

发明内容

本发明的目的在于针对现有多尺度图像融合技术存在的不足，提供一种基于期望值最大与离散小波框架的图像融合方法，能够提高融合后的图像质量，达到理想的实用效果。

为实现这样的目的，本发明技术方案的创新点在于将图像在多尺度分析中的小波框架中进行分析，分别采用不同的方法对图像的高频部分和低频部分进行融合。首先采用小波框架将待融合图像进行多尺度分解，得到小波框架的一系列高频分量和一个最低频分量，然后对图像的高频部分采用基于视觉特性的特征提取方法进行融合，对图像的低频部分采用基于期望值最大方法进行融合，最后将高频部分的融合结果和低频部分的融合结果，通过小波反变换得到最终的融合图像。

本发明的方法具体包括如下步骤：

1、采用小波框架将待融合图像进行多尺度分解：将原始图像信号输入小波基函数构造的两个滤波器进行分解，得到一个高频分量和一个低频分量，然后将这两个滤波器进行升采样处理，将分解后的低频分量作为输入信号送入经升采样处理后的滤波器中进行小波框架的两层分解，得到又一组高频分量和低频分量，同时滤波器再一次进行升采样处理，将两层分解后得到的低频分量再次作为滤波器的输入信号，可进一步得到更高的分解层数，最终得到小波框架的一系列高频分量和一个最低频分量。

2、对图像的高频部分融合处理。首先定义一个大小适当的扫描窗口，在窗口内计算边缘特征幅值和边缘连接概率。以边缘特征幅值和边缘连接概率的乘积作为多尺度空间(小波域)的高频部分的视觉特征。然后选取具有最大视觉特征的那一部分高频信息作为窗口内高频分量的融合结果。扫描窗口对所有高频分量进行融合处理，最终可以得到所有高频分量的融合结果。

3、对图像的最低频分量融合处理。首先建立图像的形成模型。图像的形成模型反映了图像与场景、成像传感器偏移量以及噪声三者之间的关系。根据期望值最大算法和图像形成模型推导出图像的叠代公式。通过这种叠代运算得到场景图像。最后对这一场景图像进行校正，保留成像传感器的对比度极性，即在估计真实场景时不对成像传感器的负对比度极性做反转。从而得到最低频分量的融合结果。

4、由以上步骤得到图像高频部分和低频部分的融合结果，通过小波框架反变换可以得到最终的融合结果。小波框架的反变换是将最后得到的低频分量和高频分量分别送入相应的构造滤波器中进行滤波，这构造滤波器也是通过小波基函数得到的。对滤波后的两组图像信号进行求和，这样就将小波框架变换层数降了一层。在作下一层反变换时，将构造滤波器进行一次降采样处理。以此方式进行叠代，就可以完成整个小波框架的反变换从而得到最终的融合图像。

本发明的图像融合方法具有如下有益效果：

在对图像基于小波框架域(多尺度域)进行分析，由于这种分析方法的空间平移不变性，有利于改善图像融合结果。在高频部分进行基于视觉特征的图像融合方法，使得在融合过程中可以尽量提取具有边缘特征的信息，而避免由于噪声干扰情况下产生的错误决策。对低频部分的融合过程考虑了传感器成像的偏移，构造了一种新的图像形成模型，从而实现了图像融合结果更加接近真实场景。基于以上三点原因，本发明的基于期望值最大与离散小波框架的图像融合能够充分提高图像的融合性能，使融合后的图像质量得到大幅度地提高，对于各种应用系统的后续处理和图像显示具有重要意义和实用价值。

附图说明

图1为本发明的基于期望值最大与离散小波框架的图像融合方法流程示意图。

图1中，待融合图像为A和B。首先分别对图像A、B进行小波框架变换，将图像信号以多尺度的形式表示出来。信号的多尺度表示有两部分，一部分是反映信号突变的高频部分，也就是信号的细节部分，另一部分是反映信号概貌的低频部分。对高频部分采用基于视觉特性的选择方法进行融合，对低频部分采用基于期望值最大的估计算法进行融合。最后将得到的高频和低频部分经小波框架反变换就可以得到最终的融合图像。

图2为本发明方法实施例的融合结果与其它不同融合方法的结果对照。

图2中，a为红外图像，b为可见光图像，c-e分别为采用DB3小波，方向金字塔和拉普拉斯金字塔的融合方法得到的图像，f为采用本发明融合方法得到的图像。

图3为受噪声干扰情况下不同融合方法的结果对照。

图3中，a为与图2中a对应的受噪声干扰的红外图像，b为与图2中b对应的受噪声干扰的可见光图像，c-e分别为采用DB3小波，方向金字塔和拉普拉斯金字塔的融合方法得到的图像，f为采用本发明融合方法得到的图像。

具体实施方式

以下结合附图和实施例对本发明的技术方案作进一步描述。

本发明的实施例中，分别对图2(a)表示的红外图像和图2(b)表示的可见光图像，以及对受噪声干扰情况下图3(a)表示的红外图像和图3(b)表示的可见光图像进行融合，采用的基于期望值最大与离散小波框架的图像融合方法按图1所示流程进行。各部分具体实施细节如下：

1采用小波框架将待融合图像进行多尺度分解

离散小波框架是对离散小波变换的修正，在对图像进行离散小波变换的时，一般需要将通过小波变换的高频和低频滤波器的输出进行降采样。小波框架的不同之处就在于没有对这一部分作降采样处理，而是将小波的滤波器系数进行修正。具体的叠代变换过程如公式(1)和公式(2)：

$w_{i+1}\left(n\right)=\underset{k}{\mathrm{\Σ}}g(2^i\·k)\·s_i(n-k)---\left(1\right)$

$s_{i+1}\left(n\right)=\underset{k}{\mathrm{\Σ}}h(2^i\·k)\·s_i(n-k)---\left(2\right)$

g(k)是小波框架的分解高频滤波器，h(k)是小波框架的分解低频滤波器。通过公式(1)和公式(2)得到小波框架的高频分量w₀，w₁，...，w_N和最低频分量s_N；N为小波框架变换的层数。

2采用基于视觉特性的特征提取方法对图像的高频部分进行融合

首先定义一个大小适当的扫描窗口，如3*3，就是在一个像素位置包含其相邻的8个像素一起进行处理。

视觉特性主要考虑了两部分：

            PI_X(m，n)＝C(m，n)·I(m，n)               (3)

PI_X(m，n)是指视觉特性，下标X表示待融合图像，C(m，n)是信号变化的幅度，也就是图像信号的高频部分的绝对值，本发明利用小波框架的高频部绝对值作为信号变化的幅度；I(m，n)是图像信号变化的拓扑性，通过如下公式计算：

            sign＝sign(C_X(m，n))                      (4)

            I(m，n)＝p_X(m，n)·(1-p_X(m，n))          (5)

如果小波框架的高频部分为大于等于零，则sign为1；如果小于零，则sign为0。p_X(m，n)为与中心位置的符号sign值相同的周围的像素点个数的概率值。

融合策略如下：

3采用基于期望值最大算法对图像的低频部分进行融合

定义图像的形成模型：

图像成像过程中，系统对真实场景的影响分为三个部分。第一部分是传感器的敏感系数，也就是真实场景的目标是否能够在传感器中感应出来；第二部分是噪声的影响，这里假定为非高斯分布的(高斯混合分布)；第三部分为成像传感器的偏移，不同的传感器对场景的偏移不同。

通过这三部分，可以得到图像的形成模型：

I_i(j)＝α_i(j)S(j)+β_i(j)+ε_i(j)   i＝1，...，q            (7)

i＝1，...，q表示不同的传感器；j表示系数位置，是一个矢量，如在x，y位置的系数可以表示为j≡(x，y)；I_i(j)表示第i个传感器的观测图像；S(j)为真实场景；α_i(j)＝±1或0表示传感器的敏感系数，当α_i(j)＝1时，表示传感器能够反映场景，当α_i(j)＝0时，表示传感器不能够反映场景，当α_i(j)＝-1时，表示传感器能够以负的对比度反映真实场景；ε_i(j)是传感器噪声，这里用K项混合高斯概率密度函数表示：

$f_{{\mathrm{\ϵ}}_i\left(j\right)}\left({\mathrm{\ϵ}}_i\left(j\right)\right)=\underset{k=1}{\overset{K}{\mathrm{\Σ}}}{\mathrm{\λ}}_{k,i}\left(j\right)\frac{1}{\sqrt{2\mathrm{\π}{\mathrm{\σ}}_{k,i}^2\left(j\right)}}\exp (-\frac{{\mathrm{\ϵ}}_i{\left(j\right)}^2}{2{\mathrm{\σ}}_{k,i}^2\left(j\right)})---\left(8\right)$

β_i(j)是指图像对场景的偏移程度。

多源图像的融合过程实际上就是对模型(7)中真实场景S(j)的估计过程，这个估计过程是在灰度系数的位置j的邻域中进行的，邻域大小的选择需要慎重，根据具体实验发现选取5×5的邻域大小是比较合适的。假设图像在小区域中，模型参数β(j)、高斯混合分布噪声的参数(λ_k，l(j)、σ_k，l ²(j))以及传感器的敏感系数α_i(j)是恒定的。

叠代过程的初始化，初始化过程和传统的融合过程非常类似，采用平均成像传感器图像的方法，假设真实场景为S。初始化S为

$\left(l\right)=\underset{i=1}{\overset{q}{\mathrm{\Σ}}}{w}_i{I}_i\left(l\right)---$ $\left(9\right)$

其中w_i＝l/q为各个成像传感器的权系数。从而可以得到各个传感器局部窗口L＝h×h的总方差系数：

${\mathrm{\σ}}_i^2=\underset{l=1}{\overset{L}{\mathrm{\Σ}}}[I_i\left(l\right)-S\left(l\right){]}^2/L---\left(10\right)$

${\mathrm{\σ}}_i^2=\underset{k=1}{\overset{K}{\mathrm{\Σ}}}{\mathrm{\λ}}_{k,i}{\mathrm{\σ}}_{k,i}^2---\left(11\right)$

其中λ_k，i＝l/K从而可以通过下式初始化混合高斯模型的方差：

${\mathrm{\σ}}_{k,i}^2=10{\mathrm{\σ}}_{k-1,i}^2---\left(12\right)$

传感器的偏移量为： ${\mathrm{\β}}_i=\frac{\underset{l=1}{\overset{L}{\mathrm{\Σ}}}{I}_i\left(l\right)}{L}---\left(13\right)$

在叠代过程开始之前，需要对图像数据进行标准化和归一化：

I′_i(j)＝(I_i(j)-μ)/H

其中I′_i和I_i分别是标准归一化后的图像和原图像，μ为整幅图像的均值，H为图像的灰度级别。

叠代过程是根据期望值最大算法中的SAGE算法推导出来的。

用S′(l)表示对S(j)的更新，类似的，α′_i，λ′_k，i和σ′_k，i是相应系数的更新。迭代过程如下：

第一步计算条件概率密度：

$g_{k,i,l}\left[I_i\left(l\right)\right]=\frac{\frac{{\mathrm{\λ}}_{k,i}}{\sqrt{2\mathrm{\π}{\mathrm{\σ}}_{k,i}^2}}\exp (-\frac{{(I_i\left(l\right)-{\mathrm{\α}}_{i}s\left(l\right)-{\mathrm{\β}}_i)}^2}{2{\mathrm{\σ}}_{k,i}^2})}{\underset{p=1}{\overset{K}{\mathrm{\Σ}}}\frac{{\mathrm{\λ}}_{p,i}}{\sqrt{2\mathrm{\π}{\mathrm{\σ}}_{p,i}^2}}\exp (-\frac{{(I_i\left(l\right)-{\mathrm{\α}}_{i}s\left(l\right)-{\mathrm{\β}}_i)}^2}{2{\mathrm{\σ}}_{p,i}^2})}---\left(14\right)$

在期望值最大算法中，这一步是计算在参数以及原始数据(图像)给定情况下，混合高斯分布第k项的条件概率密度。

第二步更新参数α_i，α_i的值在{-1，0，1}中选择以使(5)式的值最大

$Q=-\frac{1}{2}\underset{i=1}{\overset{q}{\mathrm{\Σ}}}\underset{l=1}{\overset{L}{\mathrm{\Σ}}}\underset{k=1}{\overset{K}{\mathrm{\Σ}}}[\mathrm{In}\left({\mathrm{\σ}}_{k,i}^2\right)+\frac{{(I_i\left(l\right)-{{\mathrm{\α}}^{\′}}_{i}S\left(l\right)-{\mathrm{\β}}_i)}^2}{2{\mathrm{\σ}}_{k,i}^2}]\·g_{k,i,l}\left[I_i\left(l\right)\right]---\left(15\right)$

第三步重新计算条件概率密度分布g_k，i，l，更新真实场景S(l)：

$S^{\′}\left(l\right)=\frac{\underset{i=1}{\overset{q}{\mathrm{\Σ}}}\underset{k=1}{\overset{K}{\mathrm{\Σ}}}(I_i\left(l\right)-{\mathrm{\β}}_i){{\mathrm{\α}}^{\′}}_i\frac{g_{\mathrm{kil}}\left(I_i\left(l\right)\right)}{{\mathrm{\σ}}_{\mathrm{ki}}^2}}{\underset{i=1}{\overset{q}{\mathrm{\Σ}}}\underset{k=1}{\overset{K}{\mathrm{\Σ}}}{\mathrm{\α}}_i^{{\′}^2}\frac{g_{\mathrm{ki},l}\left(I_i\left(l\right)\right)}{{\mathrm{\σ}}_{\mathrm{ki}}^2}}---\left(16\right)$

${\mathrm{\β}}_i=\frac{\underset{l=1}{\overset{L}{\mathrm{\Σ}}}\underset{k=1}{\overset{K}{\mathrm{\Σ}}}(I_i\left(l\right)-{\mathrm{\α}}_i^{\′}{S}^{\′}\left(l\right))\frac{g_{k.\mathrm{il}}\left(I_i\left(l\right)\right)}{{\mathrm{\σ}}_{k.i}^2}}{\underset{l=1}{\overset{L}{\mathrm{\Σ}}}\underset{k=1}{\overset{K}{\mathrm{\Σ}}}\frac{g_{\mathrm{ki}.l}\left(I_i\left(l\right)\right)}{{\mathrm{\σ}}_{\mathrm{ki}}^2}}---\left(17\right)$

第四步重新计算条件概率密度分布g_k，i，l、真实场景S(l)和传感器偏移量β_i，更新噪声的模型参数λ_k，i、σ_k，i ²和β_i：

${{\mathrm{\λ}}^{\′}}_{k,i}=\frac{1}{L}\underset{l=1}{\overset{L}{\mathrm{\Σ}}}{g}_{\mathrm{kil}}\left(I_i\left(l\right)\right),k=1,...,K;i=1,...,q---\left(18\right)$

${{\mathrm{\σ}}^{{\′}^2}}_{\mathrm{ki}}=\frac{\underset{l=1}{\overset{L}{\mathrm{\Σ}}}{(I_i\left(l\right)-{{\mathrm{\α}}^{\′}}_{i}s\left(l\right)-{\mathrm{\β}}_i)}^2{g}_{k.\mathrm{il}}\left(I_i\left(l\right)\right)}{\underset{l=1}{\overset{L}{\mathrm{\Σ}}}{g}_{\mathrm{kil}\left(I_i\left(l\right)\right)}},k=1,...,K;i=1,...,q---\left(19\right)$

第五步，使用新参数重复第一到第四步。

上述叠代过程中，当参数收敛到某一个特定的范围，确定融合图像为：

$S^{\′}\left(l\right)=\frac{\underset{i=1}{\overset{q}{\mathrm{\Σ}}}\underset{k=1}{\overset{K}{\mathrm{\Σ}}}{I}_i\left(l\right){\mathrm{\α}}_i^{{\′}^2}\frac{g_{\mathrm{kil}}\left(I_i\left(l\right)\right)\·C_i\left(l\right)}{{\mathrm{\σ}}_{k.l}^2}}{\underset{i=1}{\overset{q}{\mathrm{\Σ}}}\underset{k=1}{\overset{K}{\mathrm{\Σ}}}{\mathrm{\α}}_i^{{\′}^2}\frac{g_{k.\mathrm{il}}\left(I_i\left(l\right)\right)}{{\mathrm{\σ}}_{\mathrm{ki}}^2}}---\left(20\right)$

其中C_i是边缘特征影响因子， $\underset{i=1}{\overset{q}{\mathrm{\Σ}}}{C}_i=1.$

求取真实场景S的估计过程就是图像的融合过程。从图像融合的角度考虑，为了保证融合后图像有较高的对比度，融合过程应该尽量保持成像传感器的特性，保留成像传感器的对比度极性，即在估计真实场景S时不对成像传感器的负对比度极性做反转。这样做并没有减弱真实场景的信息，而且从视觉角度看，增大了对比度信息，更加有利于人眼的观察。

4采用小波框架反变换得到最终的融合结果。

$s_i\left(n\right)=\underset{k}{\mathrm{\Σ}}\tilde{h}(2^i\·n-k)\·s_{i+1}\left(n\right)+\underset{k}{\mathrm{\Σ}}\tilde{g}(2^i\·n-k)\·w_{i+1}\left(n\right)---\left(21\right)$

公式(21)表示的是小波框架的反变换过程(信号重构)。

是小波框架的重构高频滤波器，

是小波框架的重构低频滤波器。最终得到融合图像，如图2(f)和图3(f)。

表1、表2分别为可见光/红外图像以及受噪声干扰的可见光/红外图像的融合结果指标评价指标。从表中可以看出，当采用本发明提出的方法时，融合性能已经超过了DB3小波方法，方向金字塔方法以及拉普拉斯金字塔方法。

      表1可见光和红外图像的融合结果评价测度

评价测度	熵值	像素互信息	边缘互信息
				DB3小波	4.6220	0.1985	0.3719
方向金字塔	4.6071	0.2097	0.4232
				拉普拉斯金字塔	4.6974	0.2543	0.4880
本发明的方法	4.9835	0.2898	0.4848

表2可见光和红外图像在受噪声干扰情况下，融合结果的评价测度

评价测度	熵值		像素互信息		边缘互信息
								测度	相对差	测度	相对差	测度	相对差
DB3小波	5.0671	0.0963	0.3580	0.8035	0.2739	0.2635
							方向金字塔	5.0620	0.0987	0.3801	0.8126	0.3310	0.2179
拉普拉斯金字塔	5.0853	0.0826	0.4205	0.6536	0.3541	0.2744
							本发明的方法	5.0631	0.0157	0.3798	0.2370	0.3941	0.2301

Claims (1)

1、一种基于期望值最大与离散小波框架的图像融合方法，其特征在于包括如下具体步骤：

1)采用小波框架将待融合图像进行多尺度分解：将原始图像信号输入小波基函数构造的两个滤波器进行分解，得到一个高频分量和一个低频分量，然后将这两个滤波器进行升采样处理，将分解后的低频分量作为输入信号送入经升采样处理后的滤波器中进行小波框架的两层分解，得到又一组高频分量和低频分量，同时滤波器再一次进行升采样处理，将两层分解后得到的低频分量再次作为滤波器的输入信号，可进一步得到更高的分解层数，最终得到小波框架的一系列高频分量和一个最低频分量；

2)对图像的高频部分融合处理：首先定义一个扫描窗口，在窗口内计算边缘特征幅值和边缘连接概率，以边缘特征幅值和边缘连接概率的乘积作为多尺度空间的高频部分的视觉特征，然后选取具有最大视觉特征的那一部分高频信息作为窗口内高频分量的融合结果，扫描窗口对所有高频分量进行上述融合处理，最终得到所有高频分量的融合结果；

3)对图像的最低频分量融合处理：首先建立反映图像与场景、成像传感器偏移量以及噪声三者之间关系的图像形成模型，根据期望值最大算法和图像形成模型得到图像的叠代公式，通过这种叠代运算得到场景图像，然后对这一场景图像进行校正，保留成像传感器的对比度极性，即在估计真实场景时不对成像传感器的负对比度极性做反转，从而得到最低频分量的融合结果；

4)小波框架反变换得到最后融合结果：将最后得到的图像低频分量和高频分量分别送入一个同样由小波基函数构造的滤波器中滤波并对滤波后的图像信号进行求和，将小波框架变换层数降一层，将构造滤波器进行一次降采样处理后再作下一层反变换，以此类推，完成整个小波框架的反变换从而得到最终的融合图像。

Images