CN104504670A - 多尺度梯度域图像融合算法 - Google Patents

Abstract

本发明提供一种多尺度梯度域图像融合算法，将红外热像图和可见光图像相融合，其特征在于，包括以下步骤：构造梯度场步骤；修正梯度融合权重步骤；重建梯度场步骤；以及透明度法图像融合步骤。本发明所提供的多尺度梯度域图像融合算法能够同时保留红外热像图和可见光图像的温升区域和细节信息，去除图像之间相互作用产生的光晕伪影现象和冗余信息，并避免了红外热像图温升区域的颜色失真。

Description

多尺度梯度域图像融合算法

技术领域

本发明涉及计算机算法，特别涉及一种基于结构张量的多尺度梯度域图像融合算法。

背景技术

红外热成像技术能够有效地检测设备的热隐患，将电气设备的维护提升到了对设备状态过程的检测，保证了供电、用电的可靠性，因而广泛地应用于电力设备的实时监测中。红外热成像技术通过红外热像仪采集红外辐射，转换为温度场，然后表现为物体的热像分布图。热像分布图可以将温度场信息瞬间可视化，清晰直观地判别高温区域的温度情况，但是会导致非高温区域的图像细节信息缺失。而可见光图像具有完整的视觉信息，将热红外的伪彩色图像与可见光彩色图像融合可以得到更加本质、客观的互补信息，有助于检测人员或系统对故障点快速准确地定位。

近年来，关于将红外热图像与可见光图像融合的研究日益增多。该类研究大多将图像源转换到灰度域进行融合，由此得到的图像丢失了红外热像图中直观的温度信息，无法判断热点是否超出常规温度，因而为故障区域的定位增加了难度。

“Application of visible image mixing function for thermography”公开了通过混合函数将红外热像图与可见光图像融合的方法，通过二者的融合，温升区域与设备可见光图像同时呈现在一张图片中，根据该图片可以分辨出热点大概位置，但是可见光图像中的信息被严重覆盖，对热点的定位造成了干扰。

“Thermal distribution monitoring of the container data center by afast infrared image fusion technique”提出了快速红外图像融合法(FIIF)，将伪色彩热点图像与彩色可见光图像融合。该方法把源图像中的色彩信息很好地保留，使热点定位的研究取得了重要进展，然而该方法会导致部分热点信息丢失，并导致融合图像存在一定程度的冗余信息。

“A color daytime and nighttime image fusion algorithm based onIHS and Multi-Wavelet transform”提出将红外热像图与可见光图像转换到IHS空间域融合，将彩色源图像的RGB转换为色度、饱和度与亮度分量，运用加权或小波变换将亮度分量融合。由于该方法获得的融合图像改变了图像的亮度元素，且在图像恢复过程时忽略了可见光图像的色度元素，导致融合图像色彩失真，严重情况下可能会致使巡检人员对温度的误判。

发明内容

本发明是针对上述问题进行的，目的在于提供一种基于结构张量的多尺度梯度域图像融合算法，尽可能地保留红外热像图的温升区域、背景信息以及可见光图像的细节信息。

本发明为实现上述目的，采用了以下的技术方案：

本发明提供一种多尺度梯度域图像融合算法，将红外热像图和可见光图像相融合，其特征在于，包括以下步骤：构造梯度场步骤，设灰度图像为I(x,y)，该灰度图像的特征信息由梯度表示为：(,x)y是灰度图像的像素坐标，|▽I|是灰度图像的灰度的变化大小，▽I/|▽I|是梯度的方向，▽I^⊥/|▽I|表示垂直于梯度的方向，将红外热像图和可见光图像都作为源图像，设多通道图像I_m(x₁,x₂):Ω→[1,M]^N的各个标量通道来自于多个源图像I_n(x₁,x₂),n＝1,...,N，式中(x₁,x₂)表示多通道图像某一像素点的灰度值，设多通道图像中存在点a＝(x₁,x₂)和点b＝(y₁,y₂)，点a和点b的差分为I_m(x₁,x₂)-I_m(y₁,y₂)，当点a和点b的距离趋近于无穷小，差分可以用微分代替，表示为：微分dI_m的平方范数为： ${\left|{\mathrm{dI}}_m\right|}^2={\left(\frac{{\∂I}_m}{{\∂x}_1}\right)}^2{\mathrm{dx}}_1{\mathrm{dx}}_1+2(\frac{{\∂I}_m}{{\∂x}_1}\·\frac{{\∂I}_m}{{\∂x}_2}){\mathrm{dx}}_1{\mathrm{dx}}_2+{\left(\frac{{\∂I}_m}{{\∂x}_2}\right)}^2{\mathrm{dx}}_2{\mathrm{dx}}_2---\left(3\right),$ 定义矩阵G，该矩阵G的各个元素的结构形式为：则平方范数等价于：公式(5)在微分几何曲面中为第一基本形式，矩阵G为结构张量，将该结构张量G转换到梯度场表示为：

$G=\underset{n=1}{\overset{N}{\mathrm{\Σ}}}{\▿I}_n\·{{\▿I}_n}^T=\left(\begin{array}{cc}\underset{n=1}{\overset{N}{\mathrm{\Σ}}}{I_{\mathrm{nx}}}^2& \underset{n=1}{\overset{N}{\mathrm{\Σ}}}{I}_{\mathrm{nx}}\·I_{\mathrm{ny}}\\ \underset{n=1}{\overset{N}{\mathrm{\Σ}}}{I}_{\mathrm{nx}}\·I_{\mathrm{ny}}& \underset{n=1}{\overset{N}{\mathrm{\Σ}}}{I_{\mathrm{ny}}}^2\end{array}\right)=\left(\begin{array}{cc}\underset{n=1}{\overset{N}{\mathrm{\Σ}}}{\left(\frac{{\∂I}_n}{\∂x}\right)}^2& \underset{n=1}{\overset{N}{\mathrm{\Σ}}}\frac{{\∂I}_n}{\∂x}\·\frac{{\∂I}_n}{\∂y}\\ \underset{n=1}{\overset{N}{\mathrm{\Σ}}}\frac{{\∂I}_n}{\∂x}\·\frac{{\∂I}_n}{\∂y}& \underset{n=1}{\overset{N}{\mathrm{\Σ}}}{\left(\frac{{\∂I}_n}{\∂y}\right)}^2\end{array}\right)---\left(6\right),$

修正梯度融合权重步骤，采用均值法提取红外热像图中的温升区域，设源图像的特征图为S_n，梯度融合权重为：将高斯滤波G_σ(x₁,x₂,σ)的多尺度变换与红外热像图的图像梯度协方差C相结合，对图像梯度协方差C进行多尺度变换，有： $G_{\mathrm{\σ}}=\left(\begin{array}{cc}{I_{{\mathrm{nx}}_2}}^2*G_{\mathrm{\σ}}(x_1,x_2,\mathrm{\σ})& \left(I_{{\mathrm{nx}}_1}{I}_{{\mathrm{nx}}_2}\right)*G_{\mathrm{\σ}}(x_1,x_2,\mathrm{\σ})\\ \left(I_{{\mathrm{nx}}_1}{I}_{{\mathrm{nx}}_2}\right)*G_{\mathrm{\σ}}(x_1,x_2,\mathrm{\σ})& {I_{{\mathrm{nx}}_2}}^2*G_{\mathrm{\σ}}(x_1,x_2,\mathrm{\σ})\end{array}\right)---\left(8\right),$ 式中，*代表卷积运算，σ是尺度因子，将图像衔接处的结构特征设计为： $S^2={(\sqrt{{\mathrm{\λ}}_1}+\sqrt{{\mathrm{\λ}}_2})}^2+0.5{(\sqrt{{\mathrm{\λ}}_1}-\sqrt{{\mathrm{\λ}}_2})}^2---\left(9\right),$ 式中，λ₁和λ₂是经过多尺度变换后的红外热像图的图像梯度协方差的特征值，采用soft-max函数修正特征图边缘区域的梯度权重：式中，S_n是特征图边缘区域的多尺度结构特征，且最终的带有梯度权重的结构张量G为：

$G=\left(\begin{array}{cc}\underset{n=1}{\overset{N}{\mathrm{\Σ}}}{\left(w_n\frac{{\∂I}_n}{{\∂x}_1}\right)}^2& \underset{n=1}{\overset{N}{\mathrm{\Σ}}}{w_n}^2\frac{{\∂I}_n}{{\∂x}_1}\·\frac{{\∂I}_n}{{\∂x}_2}\\ \underset{n=1}{\overset{N}{\mathrm{\Σ}}}{w_n}^2\frac{{\∂I}_n}{{\∂x}_1}\·\frac{{\∂I}_n}{{\∂x}_2}& \underset{n=1}{\overset{N}{\mathrm{\Σ}}}{\left(w_n\frac{{\∂I}_n}{{\∂x}_2}\right)}^2\end{array}\right)---\left(11\right);$

重建梯度场步骤，结构张量G是半正定矩阵，特征值为λ₁和λ₂，其中最大特征值λ₁表示源图像在某一点的最大变化率，最小特征值λ₂代表该点的最小变化率，对应最大变化率和最小变化率的方向分别用特征向量e₁和e₂表示，将结构张量G对角化为： $G={\mathrm{Q\ΛQ}}^T=Q\left[\begin{array}{cc}{\mathrm{\λ}}_1& 0\\ 0& {\mathrm{\λ}}_2\end{array}\right]Q^T$ (12)，式中，Q为正交矩阵，结合公式(11)和公式(12)得到：Det(G)＝λ₁·λ₂  (13)， $\mathrm{Trace}\left(G\right)={\mathrm{\λ}}_1+{\mathrm{\λ}}_2=\underset{n=1}{\overset{N}{\mathrm{\Σ}}}{\left|{\▿I}_n\right|}^2---\left(14\right),$ 为保持源图像的基本几何性质，重建图像I_re的二阶矩阵应尽可能接近源图像的结构张量G＝QΛQ^T，于是有 $\widetilde{\mathrm{\Λ}}=\left[\begin{array}{cc}{\mathrm{\λ}}_1& 0\\ 0& 0\end{array}\right],$ 则重建图像的结构张量为 $\tilde{G}={\mathrm{\λ}}_1{e}_1{e_1}^T,$ 重建梯度场为 $\▿\tilde{I}=\sqrt{{\mathrm{\λ}}_1}{e}_1,$ 特征向量e₁由 $\mathrm{sign}(e_1,\frac{1}{N}{\mathrm{\Σ}}_{n=1}^N{\▿I}_n)$ 决定，重建图像I_re的梯度场▽I_re要尽可能地接近重建梯度场使目标函数h(I_re)最小： $h\left(I_{\mathrm{re}}\right)={\∫\∫}_{\mathrm{\Ω}}H({\▿I}_{\mathrm{re}},\▿\tilde{I})\mathrm{dxdy}---\left(15\right),$ 式中 $H({\▿I}_{\mathrm{re}},\▿\tilde{I})={|\▿I_{\mathrm{re}}-\▿\tilde{I}|}^2={(I_{\mathrm{rex}}-\▿\tilde{I})}^2+{(I_{\mathrm{rey}}-\▿\tilde{I})}^2---\left(16\right),$ 目标函数h(I_re)取极值时有 $H_{I_{\mathrm{re}}}-\frac{\∂}{\∂x}\left(\frac{\∂H}{\∂I_{\mathrm{rex}}}\right)-\frac{\∂}{\∂y}\left(\frac{\∂H}{\∂I_{\mathrm{rey}}}\right)=0---\left(17\right),$ 公式(12)的拉格朗日方程为 $2(\frac{{\∂}^2{I}_{\mathrm{re}}}{{\∂x}^2}-\frac{\∂\▿\tilde{I}}{\∂x})+2(\frac{{\∂}^2{I}_{\mathrm{re}}}{{\∂y}^2}-\frac{\∂\▿\tilde{I}}{\∂y})=0---\left(18\right),$ 简化为 ${\mathrm{\ΔI}}_{\mathrm{re}}=\mathrm{div}(\▿\tilde{I})---\left(19\right),$ 式中Δ是拉普拉斯算子，求该解泊松方程即可得到重建图像I_re；

透明度法图像融合步骤，采用R、G、B三通道方式处理重建图像和可见光图像，然后采用透明度法将重建图像和可见光图像融合，用公示表示为：I_F＝f(I_re,I_v,τ)＝(1-τ)I_re+τI_v(20)，式中，I_v是可见光图像，I_F是融合图像，τ是浑浊因子，且τ∈(0,1)，计算公式(20)得到融合图像I_F。

本发明所提供的多尺度梯度域图像融合算法，还可以具有这样的特征：其中，尺度因子σ＝0.5。

本发明所提供的多尺度梯度域图像融合算法，还可以具有这样的特征：其中，浑浊因子τ＝0.3。

发明的作用与效果

根据本发明所提供的多尺度梯度域图像融合算法，因为采用多尺度变换与温升区域特征构造梯度权重，结合带有梯度权重的结构张量与变分技术重建图像，然后采用透明度法进一步融入可见光图像的细节信息，因此能够同时保留红外热像图和可见光图像的温升区域和细节信息，去除图像之间相互作用产生的光晕伪影现象和冗余信息，并避免了红外热像图温升区域的颜色失真，解决了人眼无法识别设备热隐患的问题。

附图说明

图1是多尺度梯度域图像融合算法的流程图；

图2(a)为红外热像图，图2(b)为σ＝4时的特征图，图2(c)为σ＝2时的特征图，图2(d)为σ＝0.5时的特征图，图2(e)为σ＝0.1时的特征图；

图3(a)是人物图像的红外热像图，图3(b)是可见光图像，图3(c)是混合函数法得到的融合图像，图3(d)是采用小波变换得到的融合图像，图3(e)是采用IHS变换得到的融合图像，图3(f)采用快速红外图像融合法得到的融合图像，图3(g)是采用本实施例的多尺度梯度域融合算法得到的融合图像，图3(h)是图3(f)中I部分的放大图，图3(i)是图3(g)中II部分的放大图；

图4是发动机的红外热像图，图4(b)是可见光图像，图4(c)是混合函数法得到的融合图像，图4(d)是采用小波变换得到的融合图像，图4(e)是采用IHS变换得到的融合图像，图4(f)采用快速红外图像融合法得到的融合图像，图4(g)是采用本实施例的多尺度梯度域融合算法得到的融合图像，图4(h)是图4(f)中III部分的放大图，图4(i)是图4(g)中IV部分的放大图；以及

图5(a)是高密度数据中心开关组件的红外热像图，图5(b)是可见光图像，图5(c)是混合函数法得到的融合图像，图5(d)是采用小波变换得到的融合图像，图5(e)是采用IHS变换得到的融合图像，图5(f)采用快速红外图像融合法得到的融合图像，图5(g)是采用本实施例的多尺度梯度域融合算法得到的融合图像，图5(h)是图5(f)中V部分的放大图，图5(i)是图5(g)中VI部分的放大图。

具体实施方式

以下结合附图，对本发明所提供的多尺度梯度域图像融合算法作详细阐述。

<实施例>

图1是多尺度梯度域图像融合算法的流程图。

如图1所示，多尺度梯度域图像融合算法包括以下步骤：

步骤S-1，构造梯度场。图像的梯度值代表了图像的纹理特征，梯度值越大，图像的细节特征就越明显。梯度的矢量场由灰度图像的每个像素点的梯度构成。设灰度图像为I(x,y)，该灰度图像的特征信息由梯度表示为： $\&dtri;I={[\frac{\&PartialD;I}{\&PartialD;x},\frac{\&PartialD;I}{\&PartialD;y}]}^T---\left(1\right),$

(x,y)是灰度图像的像素坐标，|▽I|是灰度图像的灰度的变化大小，▽I/|▽I|是梯度的方向，▽I^⊥/|▽I|表示垂直于梯度的方向。

将红外热像图和可见光图像都作为源图像，在多源图像结构张量的梯度域融合中，假设多通道图像I_m(x₁,x₂):Ω→[1,M]^N的各个标量通道来自于多个源图像I_n(x₁,x₂),n＝1,...,N，式中，坐标(x₁,x₂)表示多通道图像某一像素点的灰度值。为描述多通道图像的梯度场信息，设多通道图像中存在点a＝(x₁,x₂)和点b＝(y₁,y₂)，点a和点b的差分为I_m(x₁,x₂)-I_m(y₁,y₂)，当点a和点b的距离趋近于无穷小，差分可以用微分代替，表示为：

${\mathrm{dI}}_m=\frac{\&PartialD;I_m}{{\&PartialD;x}_1}{\mathrm{dx}}_1-\frac{{\&PartialD;I}_m}{{\&PartialD;x}_2}{\mathrm{dx}}_2---\left(2\right),$

微分dI_m的平方范数为：

${\left|{\mathrm{dI}}_m\right|}^2={\left(\frac{{\&PartialD;I}_m}{{\&PartialD;x}_1}\right)}^2{\mathrm{dx}}_1{\mathrm{dx}}_1+2(\frac{{\&PartialD;I}_m}{{\&PartialD;x}_1}\&CenterDot;\frac{{\&PartialD;I}_m}{{\&PartialD;x}_2}){\mathrm{dx}}_1{\mathrm{dx}}_2+{\left(\frac{{\&PartialD;I}_m}{{\&PartialD;x}_2}\right)}^2{\mathrm{dx}}_2{\mathrm{dx}}_2---\left(3\right),$

定义矩阵G，该矩阵G的各个元素的结构形式为：

$G_{\mathrm{ij}}=\frac{{\&PartialD;I}_m}{{\&PartialD;x}_i}\&CenterDot;\frac{{\&PartialD;I}_m}{{\&PartialD;x}_j},(i,j=\mathrm{1,2})---\left(4\right),$

则矩阵G的平方范数等价于： ${\left|{\mathrm{dI}}_m\right|}^2=\left(\begin{array}{cc}{\mathrm{dx}}_1& {\mathrm{dx}}_2\end{array}\right)G_{\mathrm{ij}}\left(\begin{array}{c}{\mathrm{dx}}_1\\ {\mathrm{dx}}_2\end{array}\right)---\left(5\right),$

公式(5)在微分几何曲面中为第一基本形式，矩阵G即为结构张量，将该结构张量G转换到梯度场表示为：

$G=\underset{n=1}{\overset{N}{\mathrm{\&Sigma;}}}{\&dtri;I}_n\&CenterDot;{{\&dtri;I}_n}^T=\left(\begin{array}{cc}\underset{n=1}{\overset{N}{\mathrm{\&Sigma;}}}{I_{\mathrm{nx}}}^2& \underset{n=1}{\overset{N}{\mathrm{\&Sigma;}}}{I}_{\mathrm{nx}}\&CenterDot;I_{\mathrm{ny}}\\ \underset{n=1}{\overset{N}{\mathrm{\&Sigma;}}}{I}_{\mathrm{nx}}\&CenterDot;I_{\mathrm{ny}}& \underset{n=1}{\overset{N}{\mathrm{\&Sigma;}}}{I_{\mathrm{ny}}}^2\end{array}\right)=\left(\begin{array}{cc}\underset{n=1}{\overset{N}{\mathrm{\&Sigma;}}}{\left(\frac{{\&PartialD;I}_n}{\&PartialD;x}\right)}^2& \underset{n=1}{\overset{N}{\mathrm{\&Sigma;}}}\frac{{\&PartialD;I}_n}{\&PartialD;x}\&CenterDot;\frac{{\&PartialD;I}_n}{\&PartialD;y}\\ \underset{n=1}{\overset{N}{\mathrm{\&Sigma;}}}\frac{{\&PartialD;I}_n}{\&PartialD;x}\&CenterDot;\frac{{\&PartialD;I}_n}{\&PartialD;y}& \underset{n=1}{\overset{N}{\mathrm{\&Sigma;}}}{\left(\frac{{\&PartialD;I}_n}{\&PartialD;y}\right)}^2\end{array}\right)---\left(6\right).$

步骤S-2，修正梯度融合权重。为了保留红外热像图中的温度场信息，需要对其进行特征提取，进而结合图像特征设计梯度权重。由于红外热像图的特殊性，其伪色彩图像高温与低温区域色差明显，因此采用均值法提取红外热像图中的温升区域，不但降低了算法的复杂度，而且提高了运算效率。

假设源图像的特征图为S_n，梯度融合权重为：

$w_n=\frac{s_n}{{\mathrm{\&Sigma;}}_{k=1}^N{s}_k}---\left(7\right),$

由于红外热像图与可见光图像的色差明显，单一地使用特征图来设计梯度权重，融合衔接处的误差会导致融合图像出现光晕伪影现象，为消除光晕伪影现象，采用多尺度结构特点对图像融合边缘信息的梯度权重作进一步修正。

由于图像结构与图像梯度协方差密切相关，所以将高斯滤波G_σ(x₁,x₂,σ)的多尺度变换与红外热像图的图像梯度协方差C相结合，从而对图像梯度协方差C进行多尺度变换，得到：

$G_{\mathrm{\&sigma;}}=\left(\begin{array}{cc}{I_{{\mathrm{nx}}_2}}^2*G_{\mathrm{\&sigma;}}(x_1,x_2,\mathrm{\&sigma;})& \left(I_{{\mathrm{nx}}_1}{I}_{{\mathrm{nx}}_2}\right)*G_{\mathrm{\&sigma;}}(x_1,x_2,\mathrm{\&sigma;})\\ \left(I_{{\mathrm{nx}}_1}{I}_{{\mathrm{nx}}_2}\right)*G_{\mathrm{\&sigma;}}(x_1,x_2,\mathrm{\&sigma;})& {I_{{\mathrm{nx}}_2}}^2*G_{\mathrm{\&sigma;}}(x_1,x_2,\mathrm{\&sigma;})\end{array}\right)---\left(8\right),$

式中，*代表卷积运算，σ是尺度因子，其大小体现了图像的平滑程度。根据结构张量行列式与迹的关系，即不同区域模式中图像的平滑、边缘和角点特征，可以将图像衔接处的结构特征设计为：

$S^2={(\sqrt{{\mathrm{\&lambda;}}_1}+\sqrt{{\mathrm{\&lambda;}}_2})}^2+0.5{(\sqrt{{\mathrm{\&lambda;}}_1}-\sqrt{{\mathrm{\&lambda;}}_2})}^2---\left(9\right),$

式中，λ₁和λ₂是经过多尺度变换后的红外热像图的图像梯度协方差的特征值。

图2(a)为红外热像图，图2(b)为σ＝4时的特征图，图2(c)为σ＝2时的特征图，图2(d)为σ＝0.5时的特征图，图2(e)为σ＝0.1时的特征图。

如图2所示，为选取合适的尺度因子，本实施例对红外热像图进行了多组实验比较分析。图2(a)为红外热像图原图；如图2(b)所示，当σ取较大值4时，特征图只表现了红外热像图的概貌特征；如图2(c)所示，当σ为2时，图像细节信息逐渐明显；如图2(d)所示，当σ取较小值0.5时，图像融合边缘细节特征能够得到很好的体现；如图(e)所示，当σ取0.1时，虽然边缘细节特征清晰，但特征图的内部信息部分丢失。故本实施例选用较小的尺度因子σ＝0.5。

由于soft-max函数可以柔滑梯度场的权重，并进一步消除光晕伪影现象，故本文结合soft-max函数修正特征图像边缘区域的梯度权重，表示为：

$w_n=\frac{1}{1+{\exp }^{(S_n^{\&prime;}-S_n)/(S_n^{\&prime;}+S_n)}}---\left(10\right),$

式中，S_n是特征图像边缘区域的多尺度结构特征，且最终的带有梯度权重的结构张量G为：

$G=\left(\begin{array}{cc}\underset{n=1}{\overset{N}{\mathrm{\&Sigma;}}}{\left(w_n\frac{{\&PartialD;I}_n}{{\&PartialD;x}_1}\right)}^2& \underset{n=1}{\overset{N}{\mathrm{\&Sigma;}}}{w_n}^2\frac{{\&PartialD;I}_n}{{\&PartialD;x}_1}\&CenterDot;\frac{{\&PartialD;I}_n}{{\&PartialD;x}_2}\\ \underset{n=1}{\overset{N}{\mathrm{\&Sigma;}}}{w_n}^2\frac{{\&PartialD;I}_n}{{\&PartialD;x}_1}\&CenterDot;\frac{{\&PartialD;I}_n}{{\&PartialD;x}_2}& \underset{n=1}{\overset{N}{\mathrm{\&Sigma;}}}{\left(w_n\frac{{\&PartialD;I}_n}{{\&PartialD;x}_2}\right)}^2\end{array}\right)---\left(11\right).$

步骤S-3，重建梯度场。结构张量G是半正定矩阵，其特征值为λ₁和λ₂，代表多值图像在某一点的变化率。其中结构张量G的最大特征值λ₁表示源图像在某一点的最大变化率，最小特征值λ₂代表该点的最小变化率，对应最大变化率和最小变化率的方向分别用特征向量e₁和e₂表示。

因为G为半正定矩阵，则G可以对角化。即：

$G={\mathrm{Q\&Lambda;Q}}^T=Q\left[\begin{array}{cc}{\mathrm{\&lambda;}}_1& 0\\ 0& {\mathrm{\&lambda;}}_2\end{array}\right]Q^T---\left(12\right),$

其中，Q为正交矩阵，结合公式(11)和公式(12)得到：

Det(G)＝λ₁·λ₂              (13)，

$\mathrm{Trace}\left(G\right)={\mathrm{\&lambda;}}_1+{\mathrm{\&lambda;}}_2=\underset{n=1}{\overset{N}{\mathrm{\&Sigma;}}}{\left|{\&dtri;I}_n\right|}^2---\left(14\right),$

行列式和迹的关系在结构张量中代表了图像的区域模式，当Trace(G)＝0时，图像处于平坦区域；当Trace(G)＞0且Det(G)＝0时，图像处于边缘区域；当Trace(G)＞0且Det(G)＞0时，图像处于角点区域。

故对于一幅灰度图像(即N＝1时)，可以得知其变化率和特征向量分别为：λ₁＝|▽I|²，λ₂＝0，e₁＝▽I/|▽I|，e₂＝▽I^⊥/|▽I|。为保持源图像的基本几何性质，重建图像I_re的二阶矩阵应尽可能接近源图像的结构张量G＝QΛQ^T，于是有 $\widetilde{\mathrm{\&Lambda;}}=\left[\begin{array}{cc}{\mathrm{\&lambda;}}_1& 0\\ 0& 0\end{array}\right].$ 则重建图像的结构张量为 $\tilde{G}={\mathrm{\&lambda;}}_1{e}_1{e_1}^T,$ 重建梯度场为 $\&dtri;\tilde{I}=\sqrt{{\mathrm{\&lambda;}}_1}{e}_1,$ 特征向量e₁由 $\mathrm{sign}(e_1,\frac{1}{N}{\mathrm{\&Sigma;}}_{n=1}^N{\&dtri;I}_n)$ 决定。

重建图像I_re的梯度场▽I_re要尽可能地接近重建梯度场使目标函数h(I_re)最小： $h\left(I_{\mathrm{re}}\right)={\&Integral;\&Integral;}_{\mathrm{\&Omega;}}H({\&dtri;I}_{\mathrm{re}},\&dtri;\tilde{I})\mathrm{dxdy}---\left(15\right),$

式中 $H({\&dtri;I}_{\mathrm{re}},\&dtri;\tilde{I})={|\&dtri;I_{\mathrm{re}}-\&dtri;\tilde{I}|}^2={(I_{\mathrm{rex}}-\&dtri;\tilde{I})}^2+{(I_{\mathrm{rey}}-\&dtri;\tilde{I})}^2---\left(16\right),$

由变分技术可知，目标函数h(I_re)取极值时有

$H_{I_{\mathrm{re}}}-\frac{\&PartialD;}{\&PartialD;x}\left(\frac{\&PartialD;H}{\&PartialD;I_{\mathrm{rex}}}\right)-\frac{\&PartialD;}{\&PartialD;y}\left(\frac{\&PartialD;H}{\&PartialD;I_{\mathrm{rey}}}\right)=0---\left(17\right),$

公式(12)的拉格朗日方程为

$2(\frac{{\&PartialD;}^2{I}_{\mathrm{re}}}{{\&PartialD;x}^2}-\frac{\&PartialD;\&dtri;\tilde{I}}{\&PartialD;x})+2(\frac{{\&PartialD;}^2{I}_{\mathrm{re}}}{{\&PartialD;y}^2}-\frac{\&PartialD;\&dtri;\tilde{I}}{\&PartialD;y})=0---\left(18\right),$

简化为 $\mathrm{\&Delta;}{I}_{\mathrm{re}}=\mathrm{div}(\&dtri;\tilde{I})---\left(19\right),$

式中Δ是拉普拉斯算子，求该解泊松方程即可得到重建图像I_re。

步骤S-4，透明度法图像融合。为了使容和图像包含更均衡的可见光元素，采用透明度法将重建图像与可见光图像进一步融合，可调节式保留高温区域中可见光图像的细节信息。

采用R、G、B三通道方式处理重建图像和可见光图像，避免不同通道色彩相互影响而产生色度偏差。用公示表示为：

I_F＝f(I_re,I_v,τ)＝(1-τ)I_re+τI_v             (20)，

式中，I_v是可见光图像，I_F是融合图像，τ是浑浊因子，且τ∈(0,1)，代表重建图像与可见光图像融合的调节因子。浑浊因子可以根据可见光图像与红外热像图信息量的不同需求进行选取，若浑浊因子为0.5即为平均混合。在本实施例中取τ＝0.3，以得到更完整的热红外色彩信息。

计算公式(20)，得到最终的融合图像I_F。

图3(a)是人物图像的红外热像图，图3(b)是可见光图像，图3(c)是混合函数法得到的融合图像，图3(d)是采用小波变换得到的融合图像，图3(e)是采用IHS变换得到的融合图像，图3(f)采用快速红外图像融合法得到的融合图像，图3(g)是采用本实施例的多尺度梯度域融合算法得到的融合图像，图3(h)是图3(f)中I部分的放大图，图3(i)是图3(g)中II部分的放大图。

如图3(a)和(b)所示，从主观视觉角度出发，经过红外热图像与可见光图像的融合，能够轻易分辨出被黑色袋子挡住部分手臂的位置和手势。如图3(c)～(g)所示，对于融合效果来说，混合图像法过多的保留了红外热像图中黑色背景元素；小波变换融合效果图的颜色信息产生了严重失真；IHS变换的融合图像相对前两种融合效果具有明显的提高，然而可见光图像的颜色信息仍然存在失真现象；快速红外图像融合图像将背景部分颜色得以恢复，但是人物身体的温度信息大量丢失；采用本实施例所提供的算法，得到的融合图像尽可能大的避免了光晕伪影现象的产生，准确还原出可见光背景部分颜色信息，并且完整的保留了人物温度信息。对比图3(h)和图3(i)可以看出，本实施例的算法将温度信息保留的更加标准。

图4是发动机的红外热像图，图4(b)是可见光图像，图4(c)是混合函数法得到的融合图像，图4(d)是采用小波变换得到的融合图像，图4(e)是采用IHS变换得到的融合图像，图4(f)采用快速红外图像融合法得到的融合图像，图4(g)是采用本实施例的多尺度梯度域融合算法得到的融合图像，图4(h)是图4(f)中III部分的放大图，图4(i)是图4(g)中IV部分的放大图。

如图4所示，采用混合函数法、小波变换和HIS变换得到的融合图像都过多的融入了红外热像图中多余的背景颜色信息。快速红外图像融合的效果较前三种有所改善，但是根据融合图像的局部细节放大图(如图4(h)所示)可以明显看出，结果中热点信息严重丢失，在实际应用中可能导致设备热故障监测发生误检、漏检的现象，此类安全隐患造的后果成不堪设想。如图4(i)所示，采用本实施例的图像融合算法得到的融合图像，其热点信息与源图像最接近，目视效果较好，融合效果与其他融合方法相比较具有显著提高。

图5(a)是高密度数据中心开关组件的红外热像图，图5(b)是可见光图像，图5(c)是混合函数法得到的融合图像，图5(d)是采用小波变换得到的融合图像，图5(e)是采用IHS变换得到的融合图像，图5(f)采用快速红外图像融合法得到的融合图像，图5(g)是采用本实施例的多尺度梯度域融合算法得到的融合图像，图5(h)是图5(f)中V部分的放大图，图5(i)是图5(g)中VI部分的放大图。

如图5所示，采用混合函数法、小波变换和HIS变换得到的融合图像都过多的保留了红外热图像中蓝色部分多余的背景信息。快速红外图像融合技术虽然有了一定程度的改善，但存在于两部分热点信息之间的冗余信息并不能准确的去除，如图5(h)所示，冗余信息对热点位置的识别造成了很大的干扰。对比图5(h)和图5(i)可见，本实施例的图像融合算法尽可能多地剔除了不相关信息，提高了融合图像质量，更有利于热点定位，给监测工作提供了极大的便利，为电力设备的热点监测提供了有效的方法。

为了更有效将现有的图像融合算法与本实施例的图像融合算法相比较，对各个算法所得的的融合图像作了定量比较分析。分别采用平均梯度、标准差以及信息熵来对融合图像进行评价，平均梯度(Average gradient,AG)反映了图像对细节信息反差的表达能力，平均梯度越大，图像越清晰；标准差(Stand deviation,STD)代表了融合图像的灰度相对于灰度平均值的离散状况，标准差越大，说明图像的灰度级分布越散，对比度越高；信息熵(Information entropy,IE)用来度量融合图像的平均信息量，如果图像的信息熵越大，则融合后的信息量就越丰富，融合效果越好，一般信息熵表示为：

$H\left(x\right)=-{\mathrm{\&Sigma;}}_{i=0}^{255}p\left(i\right){\log }_{2}p\left(i\right)---\left(21\right)$

式中x为输入的图像变量，p(i)是像素灰度值i的概率。

各个图像融合算法得到的融合图像的AG、STD和IE值计算结果如表1所示。

表1

如表1所示，采用本实施例的多尺度梯度域图像融合算法得到的融合图像的AG、STD和IE值均最大，说明由多尺度梯度域图像融合算法得到的融合图像对比度高，相对其他融合图像更清晰，图中包含的信息量最大，融合效果最好。

实施例的作用与效果

根据本实施例所提供的多尺度梯度域图像融合算法，因为采用多尺度变换与温升区域特征构造梯度权重，结合带有梯度权重的结构张量与变分技术重建图像，然后采用透明度法进一步融入可见光图像的细节信息，并选择合适的尺度因子和混沌因子，因此能够同时保留红外热像图和可见光图像的温升区域和细节信息，去除图像之间相互作用产生的光晕伪影现象和冗余信息，并避免了红外热像图温升区域的颜色失真，解决了人眼无法识别设备热隐患的问题。

当然，本发明所提供的多尺度梯度域图像融合算法并不仅仅限定于以上实施例所述的内容，以上仅为本发明构思下的基本说明，而依据本发明的技术方案所作的任何等效变换，均应属于本发明的保护范围。

Claims (3)

1.一种多尺度梯度域图像融合算法，将红外热像图和可见光图像相融合，其特征在于，包括以下步骤：

构造梯度场步骤，设灰度图像为I(x,y)，该灰度图像的特征信息由梯度表示为： $\&dtri;I={[\frac{\&PartialD;I}{\&PartialD;x},\frac{\&PartialD;I}{\&PartialD;y}]}^T---\left(1\right),$

(x,y)是所述灰度图像的像素坐标，|▽I|是所述灰度图像的灰度的变化大小，▽I/|▽I|是所述梯度的方向，▽I^⊥/|▽I|表示垂直于所述梯度的方向，

将所述红外热像图和所述可见光图像都作为源图像，设多通道图像I_m(x₁,x₂):Ω→[1,M]^N的各个标量通道来自于多个所述源图像I_n(x₁,x₂),n＝1,...,N，式中(x₁,x₂)表示所述多通道图像某一像素点的灰度值，设所述多通道图像中存在点a＝(x₁,x₂)和点b＝(y₁,y₂)，所述点a和所述点b的差分为I_m(x₁,x₂)-I_m(y₁,y₂)，当所述点a和所述点b的距离趋近于无穷小，所述差分可以用微分代替，表示为：

${\mathrm{dI}}_m=\frac{\&PartialD;I_m}{\&PartialD;x_1}d{x}_1-\frac{\&PartialD;I_m}{\&PartialD;x_2}{\mathrm{dx}}_2---\left(2\right),$

所述微分dI_m的平方范数为：

${\left|{\mathrm{dI}}_m\right|}^2={\left(\frac{\&PartialD;I_m}{\&PartialD;x_1}\right)}^2{\mathrm{dx}}_1{\mathrm{dx}}_1+2(\frac{\&PartialD;I_m}{\&PartialD;x_1}\&CenterDot;\frac{\&PartialD;I_m}{\&PartialD;x_2}){\mathrm{dx}}_1{\mathrm{dx}}_2+{\left(\frac{\&PartialD;I_m}{\&PartialD;x_1}\right)}^2{\mathrm{dx}}_2{\mathrm{dx}}_2---\left(3\right),$

定义矩阵G，该矩阵G的各个元素的结构形式为：

$G_{\mathrm{ij}}=\frac{\&PartialD;I_m}{\&PartialD;x_i}\&CenterDot;\frac{\&PartialD;I_m}{\&PartialD;x_j},(i,j=\mathrm{1,2})---\left(4\right),$

则所述平方范数等价于： ${\left|d{I}_m\right|}^2=\left(\begin{array}{cc}{\mathrm{dx}}_1& {\mathrm{dx}}_2\end{array}\right)G_{\mathrm{ij}}\left(\begin{array}{c}{\mathrm{dx}}_1\\ {\mathrm{dx}}_2\end{array}\right)---\left(5\right),$

所述公式(5)在微分几何曲面中为第一基本形式，所述矩阵G为结构张量，将该结构张量G转换到梯度场表示为：

$G=\underset{n=1}{\overset{N}{\mathrm{\&Sigma;}}}\&dtri;I_n\&CenterDot;\&dtri;{I_n}^T=\left(\begin{array}{cc}\underset{n=1}{\overset{N}{\mathrm{\&Sigma;}}}{I_{\mathrm{nx}}}^2& \underset{n=1}{\overset{N}{\mathrm{\&Sigma;}}}{I}_{\mathrm{nx}}\&CenterDot;I_{\mathrm{ny}}\\ \underset{n=1}{\overset{N}{\mathrm{\&Sigma;}}}{I}_{\mathrm{nx}}\&CenterDot;I_{\mathrm{ny}}& \underset{n=1}{\overset{N}{\mathrm{\&Sigma;}}}{I_{\mathrm{ny}}}^2\end{array}\right)=\left(\begin{array}{cc}\underset{n=1}{\overset{N}{\mathrm{\&Sigma;}}}{\left(\frac{\&PartialD;I_n}{\&PartialD;x}\right)}^2& \underset{n=1}{\overset{N}{\mathrm{\&Sigma;}}}\frac{\&PartialD;I_n}{\&PartialD;x}\&CenterDot;\frac{\&PartialD;I_n}{\&PartialD;y}\\ \underset{n=1}{\overset{N}{\mathrm{\&Sigma;}}}\frac{\&PartialD;I_n}{\&PartialD;x}\&CenterDot;\frac{\&PartialD;I_n}{\&PartialD;y}& \underset{n=1}{\overset{N}{\mathrm{\&Sigma;}}}{\left(\frac{\&PartialD;I_n}{\&PartialD;y}\right)}^2\end{array}\right)---\left(6\right),$

修正梯度融合权重步骤，采用均值法提取所述红外热像图中的温升区域，设所述源图像的特征图为S_n，梯度融合权重为：

$w_n=\frac{s_n}{{\mathrm{\&Sigma;}}_{k=1}^N{s}_k}---\left(7\right),$

将高斯滤波G_σ(x₁,x₂,σ)的多尺度变换与所述红外热像图的图像梯度协方差C相结合，对所述图像梯度协方差C进行多尺度变换，有：

$C_{\mathrm{\&sigma;}}=\left(\begin{array}{cc}{I_{{\mathrm{nx}}_1}}^2*G_{\mathrm{\&sigma;}}(x_1,x_2,\mathrm{\&sigma;})& \left(I_{{\mathrm{nx}}_1}{I}_{{\mathrm{nx}}_2}\right)*G_{\mathrm{\&sigma;}}(x_1,x_2,\mathrm{\&sigma;})\\ \left(I_{{\mathrm{nx}}_1}{I}_{{\mathrm{nx}}_2}\right)*G_{\mathrm{\&sigma;}}(x_1,x_2,\mathrm{\&sigma;})& {I_{{\mathrm{nx}}_2}}^2*G_{\mathrm{\&sigma;}}(x_1,x_2,\mathrm{\&sigma;})\end{array}\right)---\left(8\right),$

式中，*代表卷积运算，σ是尺度因子，将图像衔接处的结构特征设计为： $S^2={(\sqrt{{\mathrm{\&lambda;}}_1}+\sqrt{{\mathrm{\&lambda;}}_2})}^2+0.5{(\sqrt{{\mathrm{\&lambda;}}_1}-\sqrt{{\mathrm{\&lambda;}}_2})}^2---\left(9\right),$

式中，λ₁和λ₂是经过所述多尺度变换后的所述红外热像图的图像梯度协方差的特征值，

采用soft-max函数修正所述特征图边缘区域的梯度权重：

$w_n=\frac{1}{1+{\exp }^{(S_n^{\&prime;}-S_n)/(S_n^{\&prime;}+S_n)}}---\left(10\right),$

式中，S_n是所述特征图边缘区域的多尺度结构特征，且最终的带有梯度权重的结构张量G为：

$G=\left(\begin{array}{cc}\underset{n=1}{\overset{N}{\mathrm{\&Sigma;}}}{\left(w_n\frac{\&PartialD;I_n}{\&PartialD;x_1}\right)}^2& \underset{n=1}{\overset{N}{\mathrm{\&Sigma;}}}{w_n}^2\frac{\&PartialD;I_n}{\&PartialD;x_1}\&CenterDot;\frac{\&PartialD;I_n}{\&PartialD;x_2}\\ \underset{n=1}{\overset{N}{\mathrm{\&Sigma;}}}{w_n}^2\frac{\&PartialD;I_n}{\&PartialD;x_1}\&CenterDot;\frac{\&PartialD;I_n}{\&PartialD;x_2}& \underset{n=1}{\overset{N}{\mathrm{\&Sigma;}}}{\left(w_n\frac{\&PartialD;I_n}{\&PartialD;x_2}\right)}^2\end{array}\right)---\left(11\right);$

重建梯度场步骤，所述结构张量G是半正定矩阵，特征值为λ₁和λ₂，其中最大特征值λ₁表示所述源图像在某一点的最大变化率，最小特征值λ₂代表该点的最小变化率，对应所述最大变化率和所述最小变化率的方向分别用特征向量e₁和e₂表示，将所述结构张量G对角化为：

$G={\mathrm{Q\&Lambda;Q}}^T=Q\left[\begin{array}{cc}{\mathrm{\&lambda;}}_1& 0\\ 0& {\mathrm{\&lambda;}}_2\end{array}\right]Q^T---\left(12\right),$

式中，Q为正交矩阵，结合所述公式(11)和所述公式(12)得到：

Det(G)＝λ₁·λ₂   (13)，

$\mathrm{Trace}\left(G\right)={\mathrm{\&lambda;}}_1+{\mathrm{\&lambda;}}_2=\underset{n=1}{\overset{N}{\mathrm{\&Sigma;}}}{|\&dtri;I_n|}^2---\left(14\right),$

为保持所述源图像的基本几何性质，所述重建图像I_re的二阶矩阵应尽可能接近所述源图像的结构张量G＝QΛQ^T，于是有 $\widetilde{\mathrm{\&Lambda;}}=\left[\begin{array}{cc}{\mathrm{\&lambda;}}_1& 0\\ 0& 0\end{array}\right],$ 则所述重建图像的结构张量为重建梯度场为所述特征向量e₁由决定，

所述重建图像I_re的梯度场▽I_re要尽可能地接近所述重建梯度场使目标函数h(I_re)最小： $h\left(I_{\mathrm{re}}\right)={\&Integral;\&Integral;}_{\mathrm{\&Omega;}}H(\&dtri;I_{\mathrm{re}},\&dtri;\tilde{I})\mathrm{dxdy}---\left(15\right),$

式中 $H(\&dtri;I_{\mathrm{re}},\&dtri;\tilde{I})={|\&dtri;I_{\mathrm{re}}-\&dtri;\tilde{I}|}^2={(I_{\mathrm{rex}}-\&dtri;\tilde{I})}^2+{(I_{\mathrm{rey}}-\&dtri;\tilde{I})}^2---\left(16\right),$

所述目标函数h(I_re)取极值时有 $H_{I_{\mathrm{re}}}-\frac{\&PartialD;}{\&PartialD;x}\left(\frac{\&PartialD;H}{\&PartialD;I_{\mathrm{rex}}}\right)-\frac{\&PartialD;}{\&PartialD;y}\left(\frac{\&PartialD;H}{\&PartialD;I_{\mathrm{rey}}}\right)=0---\left(17\right),$

所述公式(12)的拉格朗日方程为

$2(\frac{{\&PartialD;}^2{I}_{\mathrm{re}}}{\&PartialD;x^2}-\frac{\&PartialD;\&dtri;\tilde{I}}{\&PartialD;x})+2\left(\frac{{\&PartialD;}^2{I}_{\mathrm{re}}}{\&PartialD;y^2}\frac{\&PartialD;\&dtri;\tilde{I}}{\&PartialD;y}\right)---\left(18\right),$

简化为 $\mathrm{\&Delta;}{I}_{\mathrm{re}}=\mathrm{div}(\&dtri;\tilde{I})---\left(19\right),$

式中Δ是拉普拉斯算子，求该解泊松方程即可得到所述重建图像I_re；

透明度法图像融合步骤，采用R、G、B三通道方式处理所述重建图像和所述可见光图像，然后采用透明度法将所述重建图像和所述可见光图像融合，用公示表示为：

I_F＝f(I_re,I_v,τ)＝(1-τ)I_re+τI_v   (20)，

式中，I_v是所述可见光图像，I_F是融合图像，τ是浑浊因子，且τ∈(0,1)，计算所述公式(20)得到所述融合图像I_F。

2.根据权利要求1所述的多尺度梯度域图像融合算法，其特征在于：其中，所述尺度因子σ＝0.5。

3.根据权利要求1所述的多尺度梯度域图像融合算法，其特征在于：其中，所述浑浊因子τ＝0.3。