一种基于小波域统计信号的图像融合处理方法

摘要

本发明涉及一种基于小波域统计信号处理的图像融合方法，属于图像处理的技术领域。首先对一组初始图像进行小波分解；对低频尺度系数采用加权平均的方法得到融合后的低频尺度系数；对于高频小波系数，在三个方向上建立小波域隐马尔可夫树模型，使用期望最大算法估计融合后的高频小波系数；对融合后低频尺度系数和高频小波系数进行小波逆变换，得到融合图像。本发明在处理高频小波系数时，设各初始图像的高频小波系数为“模糊”的真实图像的高频小波系数与非高斯噪声之和，并使用了连续的“模糊因子”，这简化了图像处理方法的运算，同时能够达到较好的融合效果。

说明书

技术领域

本发明涉及一种基于小波域统计信号的图像融合处理方法，属于图像处理的技术领域。

背景技术

图像融合处理方法利用从同一场景获得的多幅初始图像之间信息的冗余性和互补性，最终得到一幅综合之前多幅初始图像优点的内容丰富的融合图像。其中，数码成像时对多聚焦图像的融合是融合处理方法的主要应用之一。

在图像融合处理方法中，基于小波变换的融合处理方法是其中很重要的一类，并已经得到广泛应用。由于人眼视觉系统是在不同尺度上处理接收到的视觉信号，而且对细节信息非常敏感，基于小波变换的图像融合处理方法就是针对上述特性，将图像分解到多个尺度上，并对低频能量和高频细节分别进行处理，从而达到比其他融合处理方法更好的融合效果。其中，基于小波域统计信号的融合处理方法，将待融合的一组初始图像建模为被“模糊”的真实图像与非高斯噪声之和，再利用估计理论估计出真实图像(即融合图像)。

尽管小波变换通常被认为是“去相关”变换，但Crouse证实了小波系数之间存在着一定的相关性，它表现为尺度内的聚集性和尺度间的持续性。小波域隐马尔可夫树模型能够充分反映小波系数在尺度内和尺度间的相关性，利用这个模型的树状结构，使用期望最大的迭代算法，可以更准确估计出融合图像。

Blum提出初始图像的高频小波系数形成模型的“模糊因子”为离散值-1，0，1，这使得期望最大算法中迭代时更新“模糊因子”的步骤略为复杂。同时，在某些场景下，“模糊因子”的离散取值显得较为单一。

发明目的

本发明的目的是提出一种基于小波域统计信号的图像融合处理方法，根据待融合的初始图像的成像特点，使用连续的“模糊因子”，并且在期望最大算法更新“模糊因子”时使用显式更新方法，以此减少图像处理的运算量。

本发明提出的基于小波域统计信号的图像融合处理方法，包括以下步骤：

(1-1)对待融合的一组初始图像中的每个图像分别进行N层小波分解，分别得到每个初始图像的低频尺度系数和每个初始图像在水平方向、垂直方向和对角方向上的N层高频小波系数；

(1-2)对初始图像的上述低频尺度系数进行加权融合，对显著性高的图像尺度系数赋予大的权值，得到融合的低频尺度系数；

(1-3)设各初始图像的高频小波系数为“模糊”的真实图像的高频小波系数与非高斯噪声之和，对上述初始图像在水平方向、垂直方向和对角方向上的N层高频小波系数分别建立隐马尔可夫树模型，使用期望最大的迭代算法对真实图像的高频小波系数进行估计，得到融合的高频小波系数；

(1-4)将上述融合的低频尺度系数和融合的高频小波系数进行小波逆变换，得到融合图像。

上述方法中，步骤(1-2)对初始图像的低频尺度系数进行加权融合的方法，包括以下各步骤：

(2-1)设z_i为第i幅初始图像的低频子带位于坐标(u，v)处的尺度系数，则初始图像的尺度系数的显著性为：${\Omega }_i=\underset{s=-2}{\overset{2}{\Sigma }}\underset{t=-2}{\overset{2}{\Sigma }}p(s+3,t+3){\left[z_i(u+s,v+t)\right]}^2,$其中，

$p=\frac{1}{48}\left(\begin{array}{ccccc}1& 1& 1& 1& 1\\ 1& 2& 2& 2& 1\\ 1& 2& 16& 2& 1\\ 1& 2& 2& 2& 1\\ 1& 1& 1& 1& 1\end{array}\right),i=1,...,q,s=-2,...,2,t=-2,...,2;$

(2-2)根据上述尺度系数的显著性，得到坐标(u，v)处各幅初始图像的加权值为：

${\lambda }_i={\Omega }_i/\underset{i=1}{\overset{q}{\Sigma }}{\Omega }_i;$

(2-3)根据上述各幅初始图像的加权值，得到坐标(u，v)处的融合低频尺度系数为：

$w=\underset{i=1}{\overset{q}{\Sigma }}{\lambda }_i{z}_i.$

上述方法中，步骤(1-3)得到融合的高频小波系数的方法，包括以下各步骤：

(3-1)设各初始图像在X方向上的高频小波系数为$z_{i,k}\left(j\right)={\beta }_{i,k}\left(j\right)w_k^{-}\left(j\right)+{ϵ}_{i,k}\left(j\right),$其中，X＝HL、LH或HH，z_i，k(j)表示第i幅初始图像在X方向上的第k棵小波树上第j个节点处的小波系数，k＝1，…，K，K为初始图像在X方向上的小波树总数，j＝1，…，P，P为每棵小波树上的节点总数，一般P＝(4^N-1)/3，w_k(j)表示真实图像在X方向上的第k棵小波树上第j个节点处的小波系数，即待估计的融合后小波系数，β_i，k(j)为“模糊因子”，-1≤β_i，k(j)≤1，ε_i，k(j)为加性非高斯噪声，由零均值和M态高斯混合模型拟合得到；

(3-2)使用期望最大算法推导出迭代公式，迭代估计上述参数，得到真实图像在X方向上的小波系数w_k(j)，即融合后的小波系数。

本发明提出的基于小波域统计信号的图像融合处理方法，具有以下优点：

(1)本发明方法在参数估计时，使用了小波域隐马尔可夫树模型，考虑了小波系数尺度内和尺度间的相关性，因此估计的参数更准确，得到的融合图像效果更好。

(2)本发明方法中的“模糊因子”采用连续值，更符合待融合的初始图像成像的实际情况，这使得对初始图像的高频小波系数建模更加准确，因此能够得到较好的融合图像。

(3)本发明方法中的期望最大算法在更新“模糊因子”时，有显式表达，使得计算简化，减少了图像处理的运算量。

附图说明

图1是本发明的图像融合处理方法的原理框图。

图2是本发明方法中建立的高频小波系数的小波树树状模型的平面图。其中，S1是一棵小波树的根节点，ρ(i)是节点i的父节点。

图3是单棵小波树的树状结构图。其中，S1是这棵小波树的根节点，ρ(i)是节点i的父节点。

图4是本发明方法中对高频小波系数进行融合的方法流程图。

图5和图6分别是两幅待融合的初始图像。其中图5在右边大表面清晰聚焦，左边小表面失焦。图6相反。

图7是利用本发明方法对图5和图6的图像进行融合后的融合图像。

具体实施方式：

本发明提出的基于小波域统计信号的图像融合处理方法，其原理框图如图1所示，包括以下步骤：

(1-1)对待融合的一组初始图像中的每个图像分别进行N层小波分解，分别得到每个初始图像的低频尺度系数和每个初始图像在水平方向、垂直方向和对角方向上的N层高频小波系数，一般N＝3或4，图1所示，N＝3；

(1-2)对初始图像的上述低频尺度系数进行加权融合，对显著性高的图像尺度系数赋予大的权值，得到融合的低频尺度系数；

(1-3)设各初始图像的高频小波系数为“模糊”的真实图像的高频小波系数与非高斯噪声之和，对上述初始图像在水平方向、垂直方向和对角方向上的N层高频小波系数分别建立隐马尔可夫树模型，使用期望最大的迭代算法对真实图像的高频小波系数进行估计，得到融合的高频小波系数；

(1-4)将上述融合的低频尺度系数和融合的高频小波系数进行小波逆变换，得到融合图像。

本发明方法中，对所有低频尺度系数进行加权融合的方法，包括以下各步骤：

(2-1)设z_i为第i幅初始图像的低频子带位于坐标(u，v)处的尺度系数，则初始图像的尺度系数的显著性为：${\Omega }_i=\underset{s=-2}{\overset{2}{\Sigma }}\underset{t=-2}{\overset{2}{\Sigma }}p(s+3,t+3){\left[z_i(u+s,v+t)\right]}^2,$其中，

$p=\frac{1}{48}\left(\begin{array}{ccccc}1& 1& 1& 1& 1\\ 1& 2& 2& 2& 1\\ 1& 2& 16& 2& 1\\ 1& 2& 2& 2& 1\\ 1& 1& 1& 1& 1\end{array}\right),i=1,...,q,s=-2,...,2,t=-2,...,2;$

(2-2)根据上述尺度系数的显著性，得到坐标(u，v)处各幅初始图像的加权值为：

${\lambda }_i={\Omega }_i/\underset{i=1}{\overset{q}{\Sigma }}{\Omega }_i;$

(2-3)根据上述各幅初始图像的加权值，得到坐标(u，v)处的融合低频尺度系数为：

$w=\underset{i=1}{\overset{q}{\Sigma }}{\lambda }_i{z}_i.$

本发明方法中，得到融合的高频小波系数的方法，包括以下各步骤：

(3-1)设各初始图像在X方向上的高频小波系数为z_i，k(j)＝β_i，k(j)w_k(j)+ε_i，k(j)，其中，X＝HL、LH或HH，HL为水平方向，LH为垂直方向，HH为对角方向，z_i，k(j)表示第i幅初始图像在X方向上的第k棵小波树上第j个节点处的小波系数，k＝1，…，K，K为初始图像在X方向上的小波树总数，j＝1，…，P，P为每棵小波树上的节点总数，一般P＝(4^N-1)/3，w_k(j)表示真实图像在X方向上的第k棵小波树上第j个节点处的小波系数，即待估计的融合后小波系数，β_i，k(j)为“模糊因子”，由于初始图像可以认为由真实图像部分模糊形成，变换到小波域则表现为对真实图像的小波系数的正负加权，所以，这里β_i，k(j)可以在[-1，1]之间连续取值，ε_i，k(j)为加性非高斯噪声，由零均值和M态高斯混合模型拟合得到；

(3-2)使用期望最大算法推导出迭代公式，迭代估计上述参数，得到真实图像在X方向上的小波系数w_k(j)，即融合后的小波系数。

上述得到融合的高频小波系数的方法中，对水平方向、垂直方向和对角方向上的小波系数分别建立模型，如图2所示。其中，S1是一棵小波树的根节点，一般选取最细尺度(J＝1)上的小波系数作为小波树的根节点；节点ρ(i)与节点i互为父、子节点。在这个树状模型中，尺度J(1≤J＜N)上的每一个节点都在尺度J+1上有四个子节点。设尺度J上的节点在子带中的坐标为(u₀，v₀)，则它在尺度J+1上的子节点坐标分别为(2u₀-1，2v₀-1)，(2u₀，2v₀-1)，(2u₀-1，2v₀)和(2u₀，2v₀)。

上述期望最大的迭代算法在估计真实图像在X方向上的小波系数时，还要涉及到一系列的参数，包括：状态概率$P_{S_{i,k}\left(j\right)}\left(m\right)=P(S_{i,k}\left(j\right)=m),$表示噪声ε_i，k(j)为状态m的概率，同时由于ε_i，k(j)w_k(j)是确定值，所以也表示z_i，k(j)为状态m的概率，这里认为在X方向上对任意的i，k，是一样的，因此可简写为P_s(j)(m)，其中m＝1，…，M；状态转移概率$a_{j,\rho \left(j\right)}^{m,n}=P\{S\left(j\right)=m|S\left(\rho \left(j\right)\right)=n\},$表示父节点ρ(j)处状态为n时其子节点j处状态为m的概率，其中n＝1，…，M；σ_m²(j)表示ε_i，k(j)为状态m时的高斯方差。上述所有的未知参数构成了待估计的模型参数集：$\Phi =\{P_{s\left(j\right)}\left(m\right),a_{j,\rho \left(j\right)}^{m,n},{\sigma }_m^2\left(j\right),{\beta }_{i,k}\left(j\right),w_k\left(j\right)\}.$在期望最大算法中，上述模型参数集被迭代更新，从而得到真实图像在X方向上的小波系数w_k(j)，即融合后的小波系数。

上述期望最大算法中，对模型参数进行估计的具体步骤包括：

(1)初始化模型参数集Φ；

(2)利用模型参数集Φ中的参数，使用前向后向算法计算P(S_i，k(j)＝m|z，Φ)和P(S_i，k(j)＝m，S_i，k(ρ(j))＝n，|z，Φ)，其中，z表示所有初始图像在X方向上的所有小波系数；

(3)根据(2)步得到的概率值，逐步更新模型参数集，并记为Φ′；

(4)若||Φ′-Φ||充分小，则循环停止；否则，Φ＝Φ′，回到(2)。

上述初始化模型参数集，包括以下步骤：

(1)状态概率P_s(j)(m)＝1/M，其中m＝1，…，M；

(2)状态转移概率$a_{j,\rho \left(j\right)}^{m,n}=1/M,$其中m＝1，…，M，n＝1，…，M；

(3)将初始图像的小波系数的加权平均值作为真实图像小波系数w_k(j)的初始值，即$w_k\left(j\right)=\underset{i=1}{\overset{q}{\Sigma }}{\lambda }_{i,k}\left(j\right)z_{i,k}\left(j\right).$其中，权值λ_i，k(j)由初始图像小波系数的显著性来决定。小波系数的显著性为${\Omega }_{i,k}\left(j\right)=\underset{s=-2}{\overset{2}{\Sigma }}\underset{t=-2}{\overset{2}{\Sigma }}p(s+3,t+3){\left[z_{i,k}(u+s,v+t)\right]}^2$。其中$p=\frac{1}{48}\left(\begin{array}{ccccc}1& 1& 1& 1& 1\\ 1& 2& 2& 2& 1\\ 1& 2& 16& 2& 1\\ 1& 2& 2& 2& 1\\ 1& 1& 1& 1& 1\end{array}\right),$(u，v)表示z_i，k(j)在子带中的坐标。这时，权值${\lambda }_{i,k}\left(j\right)={\Omega }_{i,k}\left(j\right)/\underset{i=1}{\overset{q}{\Sigma }}{\Omega }_{i,k}\left(j\right);$

(4)“模糊因子”β_i，k(j)＝1，其中i＝1，…，q，k＝1，…，K，j＝1，…，P，即每个初始图像的小波系数都可以完全反映真实图像的小波系数；

(5)高斯方差σ_m²(j)的选取。令${\sigma }_m^2\left(j\right)=\gamma {\sigma }_{m-1}^2\left(j\right),$这里取γ＝10。同时，各状态方差要满足${\sigma }^2\left(j\right)=\frac{1}{M}\underset{m=1}{\overset{M}{\Sigma }}{\sigma }_m^2\left(j\right)=\underset{i=1}{\overset{q}{\Sigma }}\underset{k=1}{\overset{K}{\Sigma }}{[z_{i,k}\left(j\right)-w_k\left(j\right)]}^2/\mathrm{qK}.$

上述使用前向后向算法计算P(S_i，k(j)＝m|z，Φ)和P(s_i，k(j)＝m，S_i，k(ρ(j))＝n，|z，Φ)时，假设了小波树之间彼此独立，所以需要对每个初始图像在X方向上的每棵小波树分别计算概率值。对于单棵小波树的计算，描述时可省略z_i，k(j)，β_i，k(j)和w_k(j)的下标，同时设小波分解的层数为N，J＝1为最细尺度，J＝N为最粗尺度，ρ(j)是节点j的父节点，c(j)是节点j的子节点的集合。前向后向算法包括以下步骤：

(1)设T_j是以j节点为根的子树。若T_i是T_j的子树，则T_j\i是T_j除去T_i。由此定义：

B_j(m)≡f(T_j|S(j)＝m，Φ)

B_j，ρ(j)(m)≡f(T_j|S(ρ(j))＝m，Φ)

B_ρ(j)\j(m)≡f(T_ρ(j)\j|S(ρ(j))＝m，Φ)

A_j(m)≡p(S(j)＝m，T_1\j|Φ)

对于M态混合高斯模型，设z(j)为状态m的概率密度函数为：

$g(z;\beta ,w,{\sigma }_m^2)=\frac{1}{\sqrt{2\pi {\sigma }_m^2}}\exp \{-\frac{{(z\left(j\right)-\beta \left(j\right)w\left(j\right))}^2}{2{\sigma }_m^2}\}$

(2)前向算法中，初始化最细尺度(J＝1)上节点j处的似然函数：

$B_j\left(m\right)=g(z\left(j\right);\beta \left(j\right),w\left(j\right),{\sigma }_m^2\left(j\right))$

由此向粗尺度循环计算，对尺度J上的节点有：

$B_{j,\rho \left(j\right)}\left(m\right)=\underset{n=1}{\overset{M}{\Sigma }}{a}_{i,\rho \left(i\right)}^{n,m}{B}_j\left(n\right)$

$B_{\rho \left(j\right)}\left(m\right)=g(z\left(\rho \left(j\right)\right);\beta \left(\rho \left(j\right)\right),w\left(\rho \left(j\right)\right),{\sigma }_m^2\left(\rho \left(j\right)\right))·\underset{j\in c\left(\rho \left(j\right)\right)}{\Pi }{B}_{j,\rho \left(j\right)}\left(m\right)$

$B_{\rho \left(j\right)\backslash j}\left(m\right)=\frac{B_{\rho \left(j\right)}\left(m\right)}{B_{j,\rho \left(j\right)}\left(m\right)}$

J＝J+1

直至J＝N时停止；

(3)后向算法中，初始化最粗尺度(J＝N)根节点处：

A₁(m)＝P_S(1)(m)

由此向细尺度循环计算，对尺度J上的节点有：

J＝J-1

$A_j\left(m\right)=\underset{n=1}{\overset{M}{\Sigma }}{a}_{j,\rho \left(j\right)}^{m,n}{A}_{\rho \left(j\right)}\left(n\right)B_{\rho \left(j\right)\backslash j}\left(n\right)$

直至J＝1时停止；

(4)由上述求出的A和B值，最终得到条件概率：

$P(S\left(j\right)=m|z,\Phi )=\frac{A_j\left(m\right)B_j\left(m\right)}{\underset{n=1}{\overset{M}{\Sigma }}{A}_j\left(n\right)B_j\left(n\right)}$

$P(S\left(j\right)=m,S\left(\rho \left(j\right)\right)=n|z,\Phi )\text{=}\frac{B_j\left(m\right)a_{j,\rho \left(j\right)}^{m,n}{A}_{\rho \left(j\right)}\left(n\right)B_{\rho \left(j\right)\backslash j}\left(n\right)}{\underset{n=1}{\overset{M}{\Sigma }}{A}_j\left(n\right)B_j\left(n\right)}$

上述更新模型参数集的过程，包括以下步骤：

(1)更新状态概率P_S(j)(m)和状态转移概率a_j，ρ(j)^m，n

${P^{\prime }}_{S\left(j\right)}\left(m\right)=\frac{\underset{i=1}{\overset{q}{\Sigma }}\underset{k=1}{\overset{K}{\Sigma }}P(S_{i,k}\left(j\right)=m|z,\Phi )}{\mathrm{Kq}}$

$a_{j,\rho \left(j\right)}^{\prime m,n}\frac{\underset{i=1}{\overset{q}{\Sigma }}\underset{k=1}{\overset{K}{\Sigma }}P(S_{i,k}\left(j\right)=m,S_{i,k}\left(\rho \left(j\right)\right)=n|z,\Phi )}{\mathrm{Kq}{P}_{S\left(\rho \left(j\right)\right)}\left(n\right)}$

(2)更新高斯方差σ_m²(j)

${\sigma }_m^{\prime 2}\left(j\right)=\frac{1}{\mathrm{Kq}{P}_{S\left(j\right)}\left(m\right)}\underset{i=1}{\overset{q}{\Sigma }}\underset{k=1}{\overset{K}{\Sigma }}{(z_{i,k}\left(j\right)-{\beta }_{i,k}\left(j\right)w_k\left(j\right))}^2·P(S_{i,k}\left(j\right)=m|z,\Phi )$

(3)更新“模糊因子”β_i，k(j)

设${\beta }_j=\frac{z_{i,k}\left(j\right)}{w_k\left(j\right)},$则${\beta }_{i,k}^{\prime }\left(j\right)=\left(\begin{array}{cc}-1& {\beta }_j<-1\\ {\beta }_j& {\beta }_j\in [-\mathrm{1,1}]\\ 1& {\beta }_j>1\end{array}\right)$

(4)更新w_k(j)值

${w^{\prime }}_k\left(j\right)=\frac{\underset{i=1}{\overset{q}{\Sigma }}\underset{m=1}{\overset{M}{\Sigma }}\frac{{{\beta }^{\prime }}_{i,k}\left(j\right)z_{i,k}\left(j\right)}{{{\sigma }^{\prime }}_m^2\left(j\right)}P(S_{i,k}\left(j\right)=m|z,\Phi )}{\underset{i=1}{\overset{q}{\Sigma }}\underset{m=1}{\overset{M}{\Sigma }}\frac{{{\beta }^{\prime }}_{i,k}\left(j\right)}{{{\sigma }^{\prime }}_m^2\left(j\right)}P(S_{i,k}\left(j\right)=m|z,\Phi )}$

以图5和图6两幅初始图像为例，使用本发明方法进行融合，得到图7的融合结果。在两幅初始图像中，两个表面分别聚焦，而得到的融合图像在两个表面都清晰聚焦，得到了比较好的融合结果。