一、图像融合简介

1 图像融合的具体步骤
(1)对全色图像和多光谱图像进行图像预处理,包括图像滤波、重采样、图像配准。
(2)将预处理后的多光谱图像fmul进行IHS变换,分别得到fmul-i(亮度)、fmul-h(色调)、fmul-s(饱和度)3个分量。
(3)对全色图像fpan和多光谱图像亮度分量fmul-i进行J层Contourlet变换,得到以下分量:
其中,AJ表示低频分量;Djk表示各尺度上各方向的高频分量;j=1,…,J,表示分解层次;k=1,…,2lj,表示第j层的各个方向。
(4)对Contourlet变换分解后的低频系数和高频系数分别进行融合,对低频系数采用基于粒子群优化的自适应加权融合规则,对高频系数采用基于粒子群优化的区域结构相似度融合规则,得到满足融合要求的低频和高频分量。
(5)通过Contourlet逆变换得到新的亮度分量fmul-i′,最后对新的亮度分量fmul-i′以及多光谱图像分量fmul-h与fmul-s进行IHS逆变换得到融合图像。

2 低频系数融合规则
Contourlet变换分解后的低频系数反映了源图像的近似信息,它包含了源图像的平均特征、光谱信息和绝大多数的能量信息,决定了融合图像的近似轮廓。低频系数融合的目的是在有效保持多光谱图像的光谱能量信息的基础上适当地融入全色图像的特征信息。目前较为常用的低频系数融合规则有加权平均法、绝对值取大(小)法以及标准差取大法,这些方法直观、简单且速度快,但是未能准确地在多光谱图像的光谱信息和全色图像的空间信息中进行优化选择。因此,本文提出了一种结合粒子群优化算法的低频系数融合规则,以线性加权值作为需要优化的决策变量,以信息熵和相对偏差之间的差值作为目标适应度函数,通过粒子群优化算法的进化迭代自适应地寻找最优加权融合系数。具体融合规则如下:

其中,AJfnew-i表示融合后的低频分量;w1,w2表示加权系数;AJfpan表示全色图像fpan的低频分量;AJfmul-i表示多光谱图像亮度分量fmul-i的低频分量。
基于粒子群优化算法的低频系数融合,以信息熵和相对偏差之间的差值作为目标适应度函数,通过粒子群优化算法迭代寻找信息熵与相对偏差的差值最大时的最优加权系数w1和w2。目标适应度函数=信息熵(E)-相对偏差(RD),其中信息熵为融合图像的信息熵,是基于信息量的评价;相对偏差为融合图像与原多光谱图像之间的相对偏差,是基于光谱性能的评价。为搜索最优的加权系数,粒子群优化算法的参数设置如下:c1为0.01,c2为0.02,惯性权重ω为0.02,种群大小为20,最大迭代次数为600。

3 高频系数融合规则
高频系数反映源图像的细节信息,主要包含边缘、纹理、亮线和区域等特征信息,而单一的图像像元不能很好地表征这些特征信息,因此需要结合多源图像对应像元之间的关系进行综合考虑与分析,通过这一区域特征的多个像元来进行表征和体现。高频系数融合的目的是在保持较优的光谱特性的基础上更多地保留原全色图像的空间细节信息。本文参考文献[19] 对高频系数采用基于区域结构相似度的融合规则进行融合的思想,在此基础上利用粒子群优化算法的全局寻优能力寻找出区域结构相似度的最优阈值p来进行高频系数的融合,融合规则如下。

(1)利用粒子群优化算法寻找高频系数间结构相似度的最优阈值p,以确定高频系数融合的方式。在本文提出的算法中,以融合图像和全色图像之间的结构相似度(SS)作为目标适应度函数,通过粒子群优化算法迭代寻找结构相似度最大时的阈值p,再采用基于区域结构相似度的融合规则进行融合。粒子群优化算法的参数设置如下:c1为0.01,c2为0.02,惯性权重ω为0.6,种群大小为20,最大迭代次数为300。
(2)对高频方向子带进行窗口运算,计算它们对应区域的结构相似度,并记录相似度的值。
(3)若相似度小于p,则采用标准差最大原则进行融合,公式为:
若相似度大于p,采用以下加权融合规则,公式为:
其中,std表示邻域窗口内系数标准差;j表示分解层次,j=1,…,J;k表示第j层的各个方向,k=1,…,2lj;SSIMjk(fpan,fmul-i)为fpan和fmul-i对应区域的结构相似度;E1jk和E2jk分别为fpan、fmul对应区域的权值。

当两幅源图像对应区域的结构相似度小于最优阈值p时,说明图像相关性较小,采用区域方差最大的融合方式可以尽可能多地增加融合图像的细节信息;当两幅源图像的对应区域的结构相似度大于阈值最优p时,说明两幅图像相关性较大且比较相似,采用加权平均融合方式可以更多地保留源图像所共有的区域结构特征。两幅图像X、Y的结构相似度定义如下:

其中,mX和mY分别表示图像X、Y的均值;σX2和σY2分别表示图像X、Y的方差;βXY表示图像X、Y的协方差;L(X,Y)、C(X,Y)和S(X,Y)分别表示图像X、Y的亮度、对比度、结构比较;C1、C2、C3为小的常数,以避免式中分母为零时出现不稳定现象。

二、部分源代码

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clc;
clear;
close all;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% 加载图像和预处理步骤
addpath QuickBird_Data  %% 数据集路径
load PAN;               %% 加载 PAN 图像
load  MS;               %% 加载 MS 图像

%% 使 PAN 和 MS 数据准备好进行处理
MSWV_db  = double(MS);
PANWV_db = double(PAN);
MS_ORG   = double(MS);

%% 调整大小，将 MS 数据上采样到 PAN 的大小

MSWV_US  = imresize(MSWV_db,  1/4, 'bicubic');
MSWV_US  = imresize(MSWV_US,  4,   'bicubic');
MSWV_DG  = MSWV_US;
PANWV_DS = imresize(PANWV_db, 1/4, 'bicubic');
PANWV_US = imresize(PANWV_DS, 4,   'bicubic');

%% 分离光谱带

R   = MSWV_US(:,:,1); 
G   = MSWV_US(:,:,2); 
B   = MSWV_US(:,:,3); 
NIR = MSWV_US(:,:,4); 

%%数据规范化

for i=1:size(MSWV_US,3)
    bandCoeffs(i)      = max(max(MSWV_US(:,:,i)));
    MSWV_US(:,:,i)     = MSWV_US(:,:,i)/bandCoeffs(i);
end

P = PANWV_DS;
panCoeff = max(max(P));
P = P/panCoeff;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% 问题定义
CostFunction=@(x) ERGAS_Index(x);          % 成本函数

nVar = 4;                                  % 决策变量的数量

VarSize = [1 nVar];                        % 决策变量矩阵的大小

VarMin = 0;                                % 变量的下界
VarMax = 1;                                % 变量的上限

三、运行结果

【图像融合】基于matlab粒子群优化自适应多光谱图像融合【含Matlab源码 004期】_matlab
【图像融合】基于matlab粒子群优化自适应多光谱图像融合【含Matlab源码 004期】_相似度_02
【图像融合】基于matlab粒子群优化自适应多光谱图像融合【含Matlab源码 004期】_图像处理_03

四、matlab版本及参考文献

1 matlab版本
2014a

2 参考文献
[1] 蔡利梅.MATLAB图像处理——理论、算法与实例分析[M].清华大学出版社，2020.
[2]杨丹,赵海滨,龙哲.MATLAB图像处理实例详解[M].清华大学出版社，2013.
[3]周品.MATLAB图像处理与图形用户界面设计[M].清华大学出版社，2013.
[4]刘成龙.精通MATLAB图像处理[M].清华大学出版社，2015.
[5]谷志鹏,贺新光.Contourlet变换与粒子群优化相耦合的遥感图像融合方法[J].计算机科学. 2016,43(S2)

3 备注
简介此部分摘自互联网，仅供参考，若侵权，联系删除