CN102855616B - 基于多尺度字典学习的图像融合方法 - Google Patents

Abstract

本发明提供了一种基于多尺度字典学习的图像融合方法，首先进行多尺度学习，对每一幅训练图像分解出S个子带，每个子带对应学习一个子字典；然后对源图像进行小波变换，得到所有源图像的子带，采用SOMP算法求解子带的稀疏表示系数并进行融合，最后进行逆小波变换，得融合图像。本发明提高了图像表示系数的稀疏度、拟合度，增强了融合图像的细节表现能力，具有更优的融合效果和更好的抑制噪声能力，字典的泛化能力也更强。

Description

基于多尺度字典学习的图像融合方法

技术领域

本发明涉及一种图像融合方法。

背景技术

图像融合是指把同一对象经多个传感器的不同成像或经单一传感器的多次成像通过特定方法加以综合，以获得更全面、更准确描述的技术，该技术在医学成像、遥感成像、机器视觉、安全监控等领域有广阔的应用前景。

由于基于稀疏表示的图像融合方法表现出了比经典的基于小波、曲波、非下采样轮廓波等融合方法更好的融合效果，从而成为当前图像融合领域的一个非常活跃的研究方向。

文献1“Yang Guang,Xu Xing-zhong,and Man Hong,Optimum image fusion viasparse representation[C],WOCC 2011-20th Annual Wireless and Optical CommunicationsConference,2011.”和文献2“Yu Nan-nan,Qiu Tian-shuang,and Bi Feng,Image featuresextraction and fusion based on joint sparse representation[J],IEEE Journal of SelectedTopics in Signal Processing,2011,5(5):1074-1082.”重点研究了基于稀疏表示的图像融合中系数融合规则的选取问题。其中文献1利用字典求得源图像对的稀疏表示向量，再通过优化源图像和融合图像的欧式距离，并用梯度对距离赋予权重融合系数。而文献2则利用联合的稀疏表示提取源图像的共有和特有特征，再按照特有特征的平均绝对值加权求得最终融合图像的稀疏系数。这两种方法均取得了比离散小波或独立成分分析等方法更优的融合效果。文献1和文献3“Yang Bin and Li Shu-tao,Pixel-level imagefusion with simultaneous orthogonal matching pursuit[J],Information Fusion,2012,1(13):10-19.”比较研究了分别采用DCT过完备字典、混合字典、学习型字典的稀疏表示系数表征源图像，再将这些系数按照选择最大的准则融合图像的三种基于稀疏表示的图像融合方法，结果表明基于学习型字典稀疏表示的融合效果最优，其原因在于基于学习型冗余字典的信号稀疏表示能够更准确有效地提取出源图像中的信息，具有更优秀的稀疏性与特征保持性。

上述基于稀疏表示的图像融合方法都包含四个基本步骤：①构造字典，②求解源图像的稀疏表示，③通过特定规则融合系数，④重构融合图像。其中，步骤①和③是决定融合效果的最关键步骤，是值得重点研究的两个方面的内容。

如前所述，目前基于字典学习的图像融合方法都是在图像域单一尺度下进行字典学习，这种学习型字典虽然可以精细地拟合数据，但是不能多尺度地分析数据。

发明内容

为了克服现有技术的不足，本发明提供一种基于多尺度字典学习的融合方法，一方面通过在小波域中学习字典，从而具有小波分析的多尺度、多方向特点；另一方面，在小波系数的基础上进一步稀疏，从而达到双重稀疏的效果。因此，基于多尺度字典学习方法既能够较小波分析方法具有更好的稀疏表示源图像的能力，又能够较单尺度学习型字典具有多尺度分析的能力，从而使得该方法能够更稀疏更有效地表示图像信息。

本发明解决其技术问题所采用的技术方案包括以下步骤：

(1)对于已经经过配准的K幅大小为M×N的源图像I₁,...,I_K，用2D小波变换分解每一幅训练图像，每一幅训练图像分解出S个子带，其中，训练图像是源图像本身或与源图像采集方式相同的图像；

(2)对所有子带，以步长为1，大小为的滑动窗口按照从左上到右下的顺序提取块，再将块拉直并依次排列组成矩阵，每一子带排列出一个矩阵；

(3)对每一个矩阵用K-SVD算法学习一个子字典D_b；

(4)保存所有学习到的子字典，即子字典D_b∈R^n×m,b＝1,2,...,S，n为字典原子的大小，m为每个子字典的原子数；

(5)对源图像进行小波变换，得到所有源图像的子带

(6)将按照步骤(2)的方法排列成矩阵得到大小为 $(\sqrt{n}\×\sqrt{n})\×[(M-\sqrt{n}+1)\×(N-\sqrt{n}+1)]$ 矩阵，从而矩阵中的第j列对应第k幅源图像经小波分解在第b个子带中第j块；

(7)采用SOMP算法并结合子字典求解下式：

$\underset{D_b,{\mathrm{ca}}_{b_j}^k}{\arg \min }{\left|\right|{\mathrm{cV}}_{b_j}^k-D_b{\mathrm{ca}}_{b_j}^k\left|\right|}_F^2\mathrm{subjectto}{\left|\right|{\mathrm{ca}}_{b_j}^k\left|\right|}_0^0\≤T,k=\{\mathrm{1,2},...,K\}$

从而得到的稀疏表示系数

（8）按照以下步骤融合源图像各个子带的稀疏表示系数：

①求解的活动水平 $A_{b_j}^k={\left|\right|{\mathrm{c\α}}_{b_j}^k\left|\right|}_1;$

②选取活动水平最大的系数作为融合图像在b子带的稀疏表示系数 的第j列 ${\mathrm{c\α}}_{b_j}^F={\mathrm{c\α}}_{b_j}^{k^{*}},$ $k^{*}=\underset{k=\mathrm{1,2},...,K}{\arg \max }\left(A_{b_j}^k\right).;$

③历所有源图像各个子带每一列的稀疏表示系数，重复步骤①和②，得到融合图像的所有子带的稀疏表示系数

④分别计算融合图像每一子带的融合矩阵得到融合图像每一子带的融合矩阵

(9)分别遍历每一子带的融合矩阵将矩阵中的每一列排列成大小的块，按照提取的顺序放到的对应位置，对同一位置的子带系数累加求平均，从而得到b＝1,2，...,S；

(10)最后将进行逆小波变换，得融合图像I_F。

本发明的有益效果是：本发明较单尺度字典学习和多尺度分析的融合方法都具有优势，在稀疏表示源图像的同时对其进行多尺度分析，进一步提高了图像表示系数的稀疏度、拟合度，增强了融合图像的细节表现能力，所以具有更优的融合效果。

本发明方法采用基于稀疏表示的多尺度字典学习策略，之所以能有较优的抑制噪声的能力，其根本原因是源图像的有用信息与噪声在稀疏表示上具有可分性，即源图像在字典上具有稀疏表示，而噪声则没有。在融合时，我们只提取在字典上有稀疏表示的源图像信息，所以达到了在融合的同时抑制噪声的效果。进一步的，本发明方法较其它基于稀疏的融合方法（OMP方法和SOMP方法）效果更优，主要因为本方法的字典具有多尺度性，能在各个尺度提取源图像的信息，抑制噪声，所以较传统单尺度学习字典的OMP、SOMP图像融合方法有更好的抑制噪声能力。

由于字典的每个原子代表源图像的一个结构原型，所以对一类图像训练好的字典可以直接应用于同类型的另一组图像，这样在融合时不必重新训练字典，这是字典在图像融合中的泛化能力。本发明方法不仅使字典中的每一个原子表征源图像的一个结构原型，还能在不同方向和尺度下表征某一结构原型，因此特征保持性更好，字典的泛化能力也更强。

附图说明

图1是红外可视图像及几种方法的融合结果，其中，(a)是可视图像，(b)是红外图像，(c)是DWT方法，(d)是OMP方法，(e)是SOMP方法，(f)是本发明方法；

图2是遥感图像及几种方法的融合结果，其中，(a)是遥感图像A，(b)是遥感图像B，(c)是DWT方法，(d)是OMP方法，(e)是SOMP方法，(f)本发明方法；

图3是几种方法对噪声图像的融合，其中，(a)是近焦距图像，(b)是远焦距图像，(c)是标准融合图像，(d)是近焦距含噪声图像，(e)是远焦距含噪声图像，(f)是DWT方法，(g)是OMP方法，(h)是SOMP方法，(i)是本发明方法；

图4是第1组遥感图像融合结果比较，其中，(a)是遥感图像A，(b)是遥感图像B，(c)是DWT方法，(d)是OMP方法，(e)是SOMP方法，(f)是本发明方法；

图5是第2组遥感图像融合结果比较，其中，(a)是遥感图像A，(b)是遥感图像B，(c)是DWT方法，(d)是OMP方法，(e)是SOMP方法，(f)是本发明方法；

图6是本发明方法框架，其中，(a)是多尺度字典学习示意图，(b)是融合过程的示意图。

具体实施方式

由于图像在不同尺度、不同方向下常常包含不同的特征，这些特征往往是图像融合需要区分和保留的突出信息，并且鉴于以往多尺度分析思想成功应用于图像融合领域，因此，我们认为在小波域多个尺度上训练的字典，将既能保持稀疏表示数据的能力，又具备多尺度分析的性能，进而在图像融合问题中能够在多尺度下用少量的系数提取显著的特征，达到比单尺度下基于字典学习的图像融合方法或者小波多尺度融合方法更优的融合效果。因此，本发明提出一种基于多尺度字典学习的图像融合新方法。

基于多尺度字典学习的融合方法具体分为两大部分：①多尺度字典学习、②融合 过程，以下分述：

多尺度字典学习如下述5步所示：

假设源图像已经经过配准，有K幅大小为M×N的源图像，并分别记作I₁,...,I_K。

(1)用2D小波变换分解每一幅训练图像，每一幅训练图像分解出S个子带，其中，训练图像可以是源图像本身，也可以是与源图像采集方式相同的图像；

(2)初始化所有子字典D_b∈R^n×m,b＝1,2,...,S，其中，S表示子字典个数，也等于源图像或训练图像经小波分解后子带个数，n为字典原子的大小，m为每个子字典的原子数；

(3)对所有子带，以步长为1，大小为的滑动窗口按照从左上到右下的顺序提取块，再将块拉直并依次排列组成矩阵，每一子带排列出一个矩阵；

(4)对每一个矩阵用K-SVD算法学习一个子字典D_b；

(5)保存所有学习到的子字典，即子字典D_b∈R^n×m,b＝1,2,...,S。

融合过程如下述6步所示：

(1)对K幅大小为M×N的源图像进行小波变换，得到所有源图像的子带其中，小波变换类型与分解级数与学习部分一致。

(2)将按照“多尺度字典学习”中步骤(3)方法排列成矩阵得到大小为 $(\sqrt{n}\×\sqrt{n})\×[(M-\sqrt{n}+1)\×(N-\sqrt{n}+1)]$ 矩阵，从而矩阵中的第j列对应第k幅源图像经小波分解在第b个子带中第j块。

(3)采用SOMP算法并结合在“多尺度字典学习”中学习到的子字典求解下式：

$\underset{D_b,{\mathrm{ca}}_{b_j}^k}{\arg \min }{\left|\right|{\mathrm{cV}}_{b_j}^k-D_b{\mathrm{ca}}_{b_j}^k\left|\right|}_F^2\mathrm{subjectto}{\left|\right|{\mathrm{ca}}_{b_j}^k\left|\right|}_0^0\≤T,k=\{\mathrm{1,2},...,K\}---\left(1\right)$

从而得到的稀疏表示系数

(4)按照以下4个步骤，融合源图像各个子带的稀疏表示系数：

①求解的活动水平，它反映了第k幅源图像在b子带第j列（块）的特征：

$A_{b_j}^k={\left|\right|{\mathrm{c\α}}_{b_j}^k\left|\right|}_1---\left(2\right)$

②选取活动水平最大的系数作为融合图像在b子带的系数 的第j列由下式获得：

${\mathrm{c\α}}_{b_j}^F={\mathrm{c\α}}_{b_j}^{k^{*}},$ $k^{*}=\underset{k=\mathrm{1,2},...,K}{\arg \max }\left(A_{b_j}^k\right).---\left(3\right)$

③遍历所有源图像各个子带每一列的稀疏表示系数，按照步骤①和②得到融合图像的所有子带的稀疏表示系数

④分别计算融合图像每一子带的融合矩阵

${\mathrm{cV}}_b^F=D_b{\mathrm{c\α}}_b^F---\left(4\right)$

最终，得到融合图像每一子带的融合矩阵

(5)分别遍历每一子带的融合矩阵将矩阵中的每一列排列成大小的块，再将这些块按照提取的顺序放到的对应位置，然而同一子带系数位置可能会出现在多个块中，因此对这些子带系数采用平均的方法处理，即对同一位置的子带系数进行累加并除以累加的次数，从而得到b＝1,2，...,S。

(6)最后将b＝1,2,...,S进行逆小波变换，即得融合图像I_F。

下面结合附图和实施例对本发明进一步说明。

本实例的图像取自图像融合数据库（http://www.Imagefusion.org/），并将本发明方法与基于单尺度字典学习的OMP、SOMP图像融合方法以及基于DWT的多尺度融合方法比较。其中，OMP，SOMP图像融合方法以及本发明方法的滑动窗口大小均为8×8像素，字典大小均为64×256，DWT和本发明方法的小波类型均选用3级db4小波。

实例的客观评价指标选择如下：①没有标准融合结果的情况，采用图像信息熵E和Q₀、Q_W、Q_E进行评价，其中，指标E用来评价融合图像中所含的信息量，越大越好；Q₀是源图像与融合图像结构化相似度的一种度量，Q_W是源图像与融合图像窗口加权的融合质量评价，Q_E反映了融合图像融合源图像边缘的情况，Q₀、Q_W、Q_E的值均在[0，1]之间，越接近1表明融合质量越好。②有标准融合结果的情况，采用RMSE和PSNR评价，以衡量融合结果与标准结果的差异，RMSE值越小越好，PSNR越大越好。

实例1.无污损图像的融合实例

本实例分别对红外可视图像对和遥感图像对进行融合，没有标准的融合结果。本实例中多尺度字典学习与融合过程的实现步骤如下：

多尺度字典学习如下述所示：

(1)用3级db4小波变换分解每一幅训练图像，训练数据采用源融合图像本身；

(2)初始化所有子字典D_b∈R^64×256,b＝1,2,...,10(3×3+1)；

(3)对所有子带，以步长为1，大小为8×8的滑动窗口按照从左上到右下的顺序提取块，再将块拉直并依次排列组成矩阵，每一子带排列出一个矩阵；

(4)对每一个矩阵用K-SVD算法学习一个子字典D_b，K-SVD算法允许误差ε取常用数值0.01；

(5)保存学习到的子字典，即子带字典D_b∈R^64×256,b＝1,2,...,10。

融合过程如下述6步所示：

(1)对源图像进行3级db4小波变换，得到所有源图像的子带 $\{{\mathrm{cI}}_b^k,k=\mathrm{1,2};b=\mathrm{1,2},...10\};$

(2)将按照“多尺度字典学习”部分步骤(3)方法排列成矩阵本实例中红外可视图像对和遥感图像对大小分别240×320和256×256，故大小分别为64×72929和64×62001；

(3)采用SOMP算法结合已经学习到的子字典求解式(1)，从而得到的稀疏表示系数 其中SOMP算法允许误差ε取常用数值0.01；

(4)按照以下4个步骤，得到所有子带的融合矩阵

①按照式(2)求解的活动水平

②选取活动水平最大的系数作为融合图像在b子带的系数 的第j列由求解式(3)获得；

③遍历所有源图像各个子带每一列的稀疏表示系数，按照步骤①和②得到融合图像的所有子带的稀疏表示系数

④按照式(4)分别计算融合图像每一子带的融合矩阵

(5)分别遍历每一子带的融合矩阵将矩阵中的每一列排列成8×8大小的块，再将这些块按照提取的顺序放到的对应位置，并取平均，即对同一位置的子带系数进行累加并除以累加的次数，从而得到b＝1,2，...,10。

(6)最后将b＝1,2,...,10进行3级db4小波逆变换，即得融合图像I_F。

图1中是对红外与可视图像的融合效果。由图1可以看出，DWT融合结果的房屋与周围景物对比度减弱，SOMP较OMP较优，但二者的清晰度较差，本发明方法则将栏杆、房屋等地物融合得更为清晰，连贯性好，视觉效果最佳。图2是对不同波段采集到的遥感图像融合。可以看出DWT的融合结果几乎无法区分建筑、陆地和树木。OMP和SOMP方法在建筑物的融合上中间亮、两边暗，存在失真现象，而本发明方法较好的体现了建筑物和树木等目标。

表1几种融合方法对无污损图像融合的客观评价

从表1的客观评价指标上看（其中粗体表示最优的指标值），对于红外可视图像的融合，本发明方法均为最优。对于遥感图像的融合，本发明方法只在评价指标E上略低于SOMP方法，其它指标优势明显。其主要原因是本发明方法为了检验多尺度学习字典这一工具在图像融合领域的有效性与准确性问题，所以在融合规则上仅是选取系数取大的传统最简单策略。

实例2.加噪图像的融合实例

在实际应用中，经常获取到包含噪声的源图像，因此如何在融合图像的过程中抑制噪声对融合结果的影响也是现在融合方法需要考虑的一个重要方面。以往大多数图像融合方法的去噪和融合处理是分开进行的，但是近几年，基于稀疏表示的图像融合方法表现出在融合同时完成去噪的优良性能，所以在本实例中进一步阐述本发明融合方法对噪声的抑制能力。

本实例的图像为有标准融合结果的多焦距图像，图4(a)(b)是无噪声的近焦和远焦图像，(c)是标准融合结果，大小均为256×256。在图4(a)(b)上分别加上均值为0，方差为σ(σ＝5,10,15,20,25)的高斯白噪声，图4(d)(e)显示的是当σ＝15时的噪声图像，(f)到(i)依次是DWT方法、OMP、SOMP和本发明方法对(d)(e)图像的融合结果。

本实例的实现步骤除在多尺度字典学习中第4步K-SVD算法允许误差ε取值和在融合过程的第3步SOMP算法允许误差ε取值与实例1不一致外，其他步骤均相同。其中，两处允许误差ε＝Cσ，C＝1.15。

由图3可以看出，DWT方法的融合结果中明显含有噪声，抗噪性能较差；基于稀疏表示的OMP、SOMP和本发明方法对噪声均有一定的抑制能力，SOMP略优于OMP方法，本发明方法效果最好，图像整体最为平滑，近焦和远焦的细节纹理也更清晰。

表2几种融合方法对噪声图像融合的客观评价

表2给出了几种融合方法的客观评价结果，其中，本发明方法在RMSE和PSNR两项指标上均为最优，显然本发明方法对噪声的抑制能力最强。

实例3.字典的泛化实例

之前两小节实例分别讨论了几种方法在无污损图像和噪声图像中的融合结果，值得注意的是字典的训练数据均是来自于源图像本身。事实上，由于字典的每个原子代表源图像的一个结构原型，所以对一类图像训练好的字典可以直接应用于同类型的另一组图像，这样在融合时不必重新训练字典，基于这一点，本小节讨论字典在图像融合中的泛化能力。

本实例的多尺度字典学习部分没有重新学习多尺度字典，而是直接将实例1中已经训练好的遥感图像多尺度字典直接用于其它两组遥感图像的融合；融合过程部分的步骤及参数设置同实例1。

融合结果如图4和图5所示。可以看出，在图4中本发明方法较其它三种方法水域与陆地更易区分，对陆地部分融合的也更为均匀；而OMP和SOMP的陆地中间亮边缘暗；DWT方法对比度较弱不利于观察。在图5中，DWT方法的融合结果显得较为杂乱，没有体现源图像的显著特点；OMP和SOMP方法在平坦区域的边界处融合的较为生硬，出现边界效应；而本发明方法的视觉效果较好。本发明方法不仅使字典中的每一个原子表征源图像的一个结构原型，还能在不同方向和尺度下表征某一结构原型，因此特征保持性更好，字典的泛化能力也更强。从表3的客观评价结果也能看出本发明提出的方法效果最好。

表3融合结果客观评价

Claims (1)

1.一种基于多尺度字典学习的图像融合方法，其特征在于包括下述步骤：

(1)对于已经经过配准的K幅大小为M×N的源图像I₁,...,I_K，用2D小波变换分解每一幅训练图像，每一幅训练图像分解出S个子带，其中，训练图像是源图像本身或与源图像采集方式相同的图像；

(2)对所有子带，以步长为1，大小为的滑动窗口按照从左上到右下的顺序提取块，再将块拉直并依次排列组成矩阵，每一子带排列出一个矩阵；

(3)对每一个矩阵用K-SVD算法学习一个子字典D_b；

(4)保存所有学习到的子字典，即子字典D_b∈R^n×m,b＝1,2,...,S，n为字典原子的大小，m为每个子字典的原子数；

(5)对源图像进行小波变换，得到所有源图像的子带

(6)将按照步骤(2)的方法排列成矩阵得到大小为矩阵，从而矩阵中的第j列对应第k幅源图像经小波分解在第b个子带中第j块；

(7)采用SOMP算法并结合子字典求解下式：

$\begin{array}{ccc}\underset{D_b,{\mathrm{ca}}_{b_j}^k}{\arg \min }{\left|\right|c{V}_{b_j}^k-D_b{\mathrm{ca}}_{b_j}^k\left|\right|}_F^2& \mathrm{subject\; to}& {\left|\right|{\mathrm{ca}}_{b_j}^k\left|\right|}_0^0\≤T,k=\{\mathrm{1,2},...,K\}\end{array}$

从而得到的稀疏表示系数

(8)按照以下步骤融合源图像各个子带的稀疏表示系数：

①求解的活动水平

②选取活动水平最大的系数作为融合图像在b子带的稀疏表示系数的第j列 $c{\mathrm{\α}}_{b_j}^F={\mathrm{c\α}}_{b_j}^{k^{*}}=\underset{k=\mathrm{1,2},...,K}{\arg \max }\left(A_{b_j}^k\right).;$

③遍历所有源图像各个子带每一列的稀疏表示系数，重复步骤①和②，得到融合图像的所有子带的稀疏表示系数

④分别计算融合图像每一子带的融合矩阵得到融合图像每一子带的融合矩阵

(9)分别遍历每一子带的融合矩阵将矩阵中的每一列排列成大小的块，按照提取的顺序放到的对应位置，对同一位置的子带系数累加求平均，从而得到 ${\mathrm{cI}}_b^F,b=\mathrm{1,2},...,S;$

(10)最后将进行逆小波变换，得融合图像I_F。