基于Shearlet变换的多通道卫星云图融合方法

摘要

本发明涉及一种基于Shearlet变换的多通道卫星云图融合方法，属于气象预测领域。基于已配准的两幅卫星云图，首先对卫星云图进行Shearlet变换，得到低频系数和高频系数；然后，对Shearlet域低频部分，利用拉普拉斯金字塔进行再次分解，对其顶层取均值，其他层取灰度绝对值大的部分后重构；在Shearlet域高频部分，先计算每个高频子图的信息熵、平均梯度和标准差这三个值，然后分别进行归一化处理，再对处理后的三个值做乘积，最后取乘积大的子图为融合子图；再对融合后的高频子图用非线性算子增强图像细的细节；最后，通过Shearlet逆变换得到最终的融合图像。本方法可以推广到三幅及以上的卫星云图的融合，实现多通道卫星云图融合，得到较高精确度的台风中心定位结果。

说明书

技术领域

本发明属于气象预测领域。具体来说，涉及一种以提高台风中心定位精度为目的的基于Shearlet变换的多通道卫星云图融合方法。 

背景技术

气象卫星云图在天气监测预报和大气环境检测中发挥了极其重要的作用，特别是对一些极端气象灾害监测起到了关键作用。因此，对卫星云图进行后续分析处理，能更好地获取大气、陆地、海洋、云层等信息，为监测预报提供可靠的数据支持，而且能提高预报的自动化和准确性，具有重要的现实意义。 

我国风云2号C星（FY-2C）地球静止轨道气象卫星通过一个可见光通道、三个红外通道和一个水气通道的扫描辐射仪接收来自地球的可见光、红外和水气辐射，五个通道每半小时各获取一幅覆盖1/3地球的全景云图，多时次的频繁观测特别适合于检测生命史较短而危害又极大的暴雨、台风、沙尘暴等灾害性天气的发生和发展。但是，各通道的图像成像原理不同，得到的数据信息也有所不同，从单一通道的卫星云图中获取的信息有一定的局限性，不利于反映观测目标的特点。图像融合方法结合了不同通道的卫星云图信息，能够提供更全面的云图信息，利于获取更可靠的数据，提高预报和监测的精度。因此，国内外学者对多通道卫星云图的融合技术都进行了不断地探索。 

小波分析理论经过长足的发展，在图像融合领域已经得到了极为广泛的应用。A.Abd-Elrahman等人提出了一种在利用小波变换进行卫星云图融合时提高其云相关阴影区域的增强方法，并保留了细节信息，其结果有效地提高了云图质量。Lee,Y.等人提出一种新的小波域卫星图像融合算法，考虑了每幅源图像的强度和光谱范围以及相关光谱响应。用高斯函数之和表示每个通道的光谱响应，然后用高斯函数建模调整图像的空间和光 谱分辨率。该方法融合结果的PSNR（Peak Signal-Noise Ratio）值、均方根误差和相关系数比传统方法的效果更优。Yang W.等人将压缩感知（CS）引入到卫星图像的融合算法中，提出一个基于对称部分的b样条小波的融合算法CS-FWT-PCA，该方法使用Hama Da矩阵作为测量矩阵，稀疏度自适应匹配追踪算法（SAMP）作为重建算法，采用一种改进的基于局部方差的融合规则，得到优于传统融合方法的融合效果。 

然而，小波变换存在着一定的局限性，由于只有较少的方向数，即水平、垂直及对角三个方向，因而只能捕获到有限方向上的信息，容易造成信息的丢失。针对小波分离的局限性和不足，近几年来多分辨分析理论得到了进一步的发展，多尺度几何分析工具应运而生。多尺度几何分析工具不仅具有小波的多分辨率，同时还具有多尺度、良好的时频局部特性以及高度的方向性和各向异性。常见的多尺度几何分析工具包括：Bandelet、Ridgelet、Curvelet、Contourlet、NSCT（NonSubsampled Contourlet Transform）和Shearlet等。其中，Curvelet变换比小波更加适合分析二维图像中曲线或直线的边缘特征，具有更高的逼近精度和稀疏表达能力，将Curvelet变换引入图像融合能够更好地提取原始图像的特征，为融合图像提供更多的信息。Shutao Li等人提出了一种结合Curvelet和小波变换的多聚焦图像融合算法，该图像融合结果优于其他任何单独的多尺度融合方法。但是Curvelet变换算法的冗余度较高，没有基于临界采样的滤波器组，而且有吉布斯现象。此外，Contourlet变换是由严格采用的扇形滤波器和重采样实现的，通常由拉普拉斯金字塔变换和一个方向滤波器构成。能用不同尺度、不同频率的子带更精确地捕获图像中的分段二次连续曲线，从而使图像的边缘能量更加集中。Miao Qiguang等人提出了一种基于Contourlet变换的图像融合方法，对高频部分进行区域能量比较，并结合一致性检验，得到在边缘保持和纹理信息方面优于小波变换方法和拉普拉斯金字塔的图像融合方法的结果。Juan Lu等人提出了一种基于NSCT和能量熵的图像融合算法，该算法的融合结果具有更丰富的方向信息和很强 的噪声鲁棒性。但Contourlet变换应用于图像融合时容易引入伪吉布斯现象，克服该现象所需时间过长且数据量大。 

而Shearlet变换引入Shear滤波器，该滤波器没有方向数的约束，还可以用一个窗函数表示。Shearlet变换不仅可以检测到所有的奇异点，而且可以自适应跟踪奇异曲线的方向，反变换时也不需要进行逆扇形滤波器组变换。Shearlet变换既克服了小波变换在图像融合应用中信息丢失的不足，还成功摆脱了滤波时方向数的限制。因此，随着图像处理技术的发展，Shearlet变换越来越受到广大研究人员的关注和重视，已经成为研究热点之一。Qi-guang Miao等人利用Shearlet变换的方向性、定位性、各向异性、多尺度等优势，将Shearlet变换用于图像融合，其融合结果比其他方法包含更多的细节和较小的失真信息。Cheng S.等人提出一种基于Shearlet变换和PCNN（Pulse Coupled Neural Network）的图像融合算法。提取每个方向Shear矩阵的梯度特征，用小波进行多尺度分解，并用PCNN融合高频系数，取得了良好的融合效果。Guorong G.等人提出一种新的基于NSST（Non-Subsampled Shearlet Transform）变换的多聚焦图像融合算法，针对多聚焦图像的小方误差特点提出集中区域信息的融合规则，实验结果显示该方法的视觉质量和客观评价明显优于离散小波变换的融合结果。 

发明内容

本发明提供一种基于Shearlet变换的多通道卫星云图融合方法，目的是提供一种综合考虑图像各评价参数以保留云图高信息度的融合方法，实现较好地台风中心定位，并可推广到多通道卫星云图的融合。 

本发明采取的技术方案是，包括下列步骤： 

步骤1、对配准后的源图像A和B，图像大小为M×N，分别进行Shearlet变换，分解层数为W，分解方向个数为T，T=2^r,r∈Z^*，得到高频系数SH_A和SH_B、低频系数SL_A和SL_B； 

步骤2、分别对低频系数SL_A和SL_B做拉普拉斯金字塔分解，分解层数为Q，得到分解图像LA和LB，第q层子图分别为LA_q和LB_q，1≤q≤Q； 

步骤3、对拉普拉斯金字塔顶层子图LA_Q和LB_Q用均值法进行融合，得到融合结果LF_Q为： 

$>{\mathrm{LF}}_Q(i,j)=\frac{{\mathrm{LA}}_Q(i,j)+{\mathrm{LB}}_Q(i,j)}{2}$>

其中，1≤i≤CL_Q，1≤j≤RL_Q，CL_Q是分解子图第Q层图像的行数，RL_Q是分解子图第Q层图像的列数； 

步骤4、对拉普拉斯金字塔其他层子图LA_q和LB_q用灰度绝对值取大的融合规则融合，1≤q≤Q-1，则融合结果LF_q为： 

$>{\mathrm{LF}}_q(i,j)=\left(\begin{array}{cc}{\mathrm{LA}}_q(i,j),& \left|{\mathrm{LA}}_q(i,j)\right|\ge \left|{\mathrm{LB}}_q(i,j)\right|\\ {\mathrm{LB}}_q(i,j),& \left|{\mathrm{LA}}_q(i,j)\right|<\left|{\mathrm{LB}}_q(i,j)\right|\end{array}\right)$>

步骤5、将融合后得到的拉普拉斯金字塔LF进行重构，得到低频部分的融合结果TL_F； 

步骤6、在Shearlet变换域的高频部分，分别求每层每个方向子图的信息熵、平均梯度和标准差，记第w层第t个方向的高频子图分别为和1≤w≤W，1≤t≤T，其大小和原图大小一样是M×N，则其信息熵E为： 

$>E=-\underset{i=0}{\overset{L-1}{\Sigma }}{P}_i{\log }_2{P}_i$>

其中，P_i表示子图中像素灰度值为i的概率,L表示图像中灰度级的个数，高频子图的平均梯度表示为 

$>\bar{G}=\frac{\underset{\mathrm{ii}=1}{\overset{M-1}{\Sigma }}\underset{\mathrm{jj}=1}{\overset{N-1}{\Sigma }}\sqrt{\frac{{\left(\frac{{\partial \mathrm{SH}}_w^t(x_{\mathrm{ii}},y_{\mathrm{jj}})}{{\partial x}_{\mathrm{ii}}}\right)}^2+{\left(\frac{{\partial \mathrm{SH}}_w^t(x_{\mathrm{ii}},y_{\mathrm{jj}})}{{\partial y}_{\mathrm{jj}}}\right)}^2}{2}}}{(M-1)(N-1)}$>

其中，表示高频子图或中x_ii行y_jj列的像素点，1≤ii≤M,1≤jj≤N，高频子图的标准差σ表示为： 

$>\sigma =\sqrt{\underset{\mathrm{ii}=1}{\overset{M}{\Sigma }}\underset{\mathrm{jj}=1}{\overset{N}{\Sigma }}\frac{{({\mathrm{SH}}_w^t(x_{\mathrm{ii}},y_{\mathrm{jj}})-\bar{h})}^2}{M\times N}}$>

其中，表示高频子图的灰度均值； 

步骤7、对高频子图和相对应的信息熵E、平均梯度标准差σ分别进行归一化处理，得到归一化后的信息熵E_g、平均梯度标准差σ_g，选择三者乘积值大的作为融合后的高频子图，即 

$>{\left({\mathrm{SH}}_F\right)}_w^t=\left(\begin{array}{cc}{\left({\mathrm{SH}}_A\right)}_w^t,& {\left(E_g\right)}_A\times {\left({\bar{G}}_g\right)}_A\times {\left({\sigma }_g\right)}_A\ge {\left(E_g\right)}_B\times {\left({\bar{G}}_g\right)}_B\times {\left({\sigma }_g\right)}_B\\ {\left({\mathrm{SH}}_B\right)}_w^t,& {\left(E_g\right)}_A\times {\left({\bar{G}}_g\right)}_A\times {\left({\sigma }_g\right)}_A<{\left(E_g\right)}_B\times {\left({\bar{G}}_g\right)}_B\times {\left({\sigma }_g\right)}_B\end{array}\right)$>

步骤8、对融合后的高频子图进行基于离散平稳小波变换的非线性增 强处理，设高频子图的所有像素灰度中绝对值最大的为maxh，则可以得到增强后的高频子图为： 

$>{\left({E\_\mathrm{SH}}_F\right)}_w^t(\mathrm{ii},\mathrm{jj})=a·\max h\left\{\mathrm{sigm}\right[c({\mathrm{Sh}}_w^t(\mathrm{ii},\mathrm{jj})-b)]-\mathrm{sigm}[-c({\mathrm{Sh}}_w^t(\mathrm{ii},\mathrm{jj})+b)\left]\right\}$>

其中，b=0.35，c=20；a=1/(d₁-d₂)，d₁=sigm(c×(1+b))，d₂=sigm(-c×(1-b))； 

$>{\mathrm{Sh}}_w^t(\mathrm{ii},\mathrm{jj})={\left({\mathrm{SH}}_f\right)}_w^t(\mathrm{ii},\mathrm{jj})/\max h;$>

步骤9、对融合处理后的Shearlet系数值进行Shearlet逆变换，得到最终的融合图像F。 

本发明所述的步骤1中，所述配准后的源图像A和B选用我国气象卫星FY-2C返回的卫星云图图像，该卫星有5个通道：红外1通道、红外2通道、水气通道、红外4通道、可见光通道，任选其中两个通道的图像进行配准。 

本发明还包括：在前面的融合结果上再加上其他通道的云图，实现三幅及以上的云图融合，达到多通道卫星云图融合的效果。 

本发明将Shearlet变换应用于多通道卫星云图的融合，基于已配准的两幅云图，结合拉普拉斯金字塔分解的融合方法，提出一种综合考虑图像各评价参数以保留云图高信息度的融合规则，实现较好地台风中心定位，并可推广到多通道卫星云图的融合。 

发明所述的技术方案能够充分融合各个通道的有用信息，很好地实现多通道台风云图融合，最大限度的保留各个通道云图的细节，保持融合后图像的清晰度。通过利用融合后的云图对有眼和无眼台风进行了中心定位，得到较高精确度的台风中心定位结果，可以说明融合效果有良好的实用价值。 

附图说明

图1为本发明基于Shearlet变换的多通道卫星云图融合方法的流程图； 

图2（a）为我国气象卫星FY-2C返回的5通道卫星云图中的红外1通道云图； 

图2（b）为我国气象卫星FY-2C返回的5通道卫星云图中的红外2通道云图; 

图2（c）为我国气象卫星FY-2C返回的5通道卫星云图中的水气通道云图; 

图2（d）为我国气象卫星FY-2C返回的5通道卫星云图中的红外4通道云图; 

图2（e）为我国气象卫星FY-2C返回的5通道卫星云图中的可见光通道云图; 

图3(a)为2005年8月31日12时00分台风“泰利”的红外2通道和水气通道云图（有眼台风）的多通道卫星云图融合实验结果中的红外2通道云图； 

图3(b)为2005年8月31日12时00分台风“泰利”的红外2通道和水气通道云图（有眼台风）的多通道卫星云图融合实验中的水气通道云图； 

图3(c)为2005年8月31日12时00分台风“泰利”的红外2通道和水气通道云图（有眼台风）的多通道卫星云图融合实验中的拉普拉斯金字塔的融合结果； 

图3(d)为2005年8月31日12时00分台风“泰利”的红外2通道和水气通道云图（有眼台风）的多通道卫星云图融合实验结果中的经典离散正交小波的融合结果； 

图3(e)为2005年8月31日12时00分台风“泰利”的红外2通道和水气通道云图（有眼台风）的多通道卫星云图融合实验结果中的Curvelet的融合结果； 

图3(f)为2005年8月31日12时00分台风“泰利”的红外2通道和水气通道云图（有眼台风）的多通道卫星云图融合实验结果中的Contourlet的融合结果； 

图3(g)为2005年8月31日12时00分台风“泰利”的红外2通道和 水气通道云图（有眼台风）的多通道卫星云图融合实验结果中的NSCT的融合结果； 

图3(h)为2005年8月31日12时00分台风“泰利”的红外2通道和水气通道云图（有眼台风）的多通道卫星云图融合实验结果中的本发明算法的融合结果； 

图4（a）为图3(c)中融合结果的部分放大图像； 

图4(b)为图3(d)中融合结果的部分放大图像； 

图4(c)为图3(e)中融合结果的部分放大图像； 

图4(d)为图3(f)中融合结果的部分放大图像； 

图4(e)为图3(g)中融合结果的部分放大图像； 

图4(f)为图3(h)中融合结果的部分放大图像； 

图5（a）为图3(c)中“泰利”云图的融合结果截取云图，便于进行台风中心定位； 

图5(b)为图3(d)中“泰利”云图的融合结果截取云图，便于进行台风中心定位； 

图5(c)为图3(e)中“泰利”云图的融合结果截取云图，便于进行台风中心定位； 

图5(d)为图3(f)中“泰利”云图的融合结果截取云图，便于进行台风中心定位； 

图5(e)为图3(g)中“泰利”云图的融合结果截取云图，便于进行台风中心定位； 

图5(f)为图3(h)中“泰利”云图的融合结果截取云图，便于进行台风中心定位； 

图6(a)为利用图3(a)中“泰利”云图的红外2通道云图的台风中心定位结果示意图； 

图6(b)为利用图3(b)中“泰利”云图的水气通道云图的台风中心定位结果示意图； 

图6(c)为利用图3(c)中“泰利”云图的融合结果图像的台风中心定位结果示意图； 

图6(d)为利用图3(d)中“泰利”云图的融合结果图像的台风中心定位结果示意图； 

图6(e)为利用图3(e)中“泰利”云图的融合结果图像的台风中心定位结果示意图； 

图6(f)为利用图3(f)中“泰利”云图的融合结果图像的台风中心定位结果示意图； 

图6(g)为利用图3(g)中“泰利”云图的融合结果图像的台风中心定位结果示意图； 

图6(h)为利用图3(h)中“泰利”云图的融合结果图像的台风中心定位结果示意图； 

图7(a)为2006年5月11日00时00分台风“珍珠”的红外1通道云图（无眼台风）； 

图7(b)为2006年5月11日00时00分台风“珍珠”的水气通道云图（无眼台风）； 

图7(c)为2006年5月11日00时00分台风“珍珠”的红外1通道和水气通道云图（无眼台风）的多通道卫星云图的拉普拉斯金字塔的融合实验结果； 

图7(d)为2006年5月11日00时00分台风“珍珠”的红外1通道和水气通道云图（无眼台风）的多通道卫星云图的经典离散正交小波的融合实验结果； 

图7(e)为2006年5月11日00时00分台风“珍珠”的红外1通道和水气通道云图（无眼台风）的多通道卫星云图的Curvelet的融合实验结果； 

图7(f)为2006年5月11日00时00分台风“珍珠”的红外1通道和水气通道云图（无眼台风）的多通道卫星云图的Contourlet的融合实验结果； 

图7(g)为2006年5月11日00时00分台风“珍珠”的红外1通道和水气通道云图（无眼台风）的多通道卫星云图的NSCT的融合实验结果； 

图7(h)为2006年5月11日00时00分台风“珍珠”的红外1通道和水气通道云图（无眼台风）的多通道卫星云图的本发明的融合实验结果； 

图8（a）为图7（c）中融合结果的部分放大图像； 

图8（b）为图7（d）中融合结果的部分放大图像； 

图8（c）为图7（e）中融合结果的部分放大图像； 

图8（d）为图7（f）中融合结果的部分放大图像； 

图8（e）为图7(g)中融合结果的部分放大图像； 

图8（f）为图7(h)中融合结果的部分放大图像； 

图9(a)为图7(c)中“珍珠”云图的融合结果截取云图，便于进行台风中心定位； 

图9（b）为图7（d）中“珍珠”云图的融合结果截取云图，便于进行台风中心定位； 

图9（c）为图7（e）中“珍珠”云图的融合结果截取云图，便于进行台风中心定位； 

图9（d）为图7（f）中“珍珠”云图的融合结果截取云图，便于进行台风中心定位； 

图9（e）为图7(g)中“珍珠”云图的融合结果截取云图，便于进行台风中心定位； 

图9（f）为图7(h)中“珍珠”云图的融合结果截取云图，便于进行台风中心定位； 

图10(a)为利用图7(a)中“珍珠”的红外1通道云图的台风中心定位结果示意图； 

图10（b）为利用图7（b）中“珍珠”的水气通道云图的台风中心定位结果示意图； 

图10（c）为利用图7（c）中“珍珠”云图的拉普拉斯金字塔的融合 结果图像的台风中心定位结果示意图； 

图10（d）为利用图7（d）中“珍珠”云图的经典离散正交小波的融合结果图像的台风中心定位结果示意图； 

图10（e）为利用图7（e）中“珍珠”云图的Curvelet的融合结果图像的台风中心定位结果示意图； 

图10（f）为利用图7（f）中“珍珠”云图的Contourlet的融合结果图像的台风中心定位结果示意图； 

图10(g)为利用图7(g)中“珍珠”云图的NSCT的融合结果图像的台风中心定位结果示意图； 

图10(h)为利用图7(h)中“珍珠”云图的本发明的融合结果图像的台风中心定位结果示意图。 

具体实施方式

本发明提出一种基于Shearlet变换的多通道卫星云图融合方法。基于已配准的两幅卫星云图，首先对卫星云图进行Shearlet变换，得到低频系数和高频系数；然后，对Shearlet域低频部分，利用拉普拉斯金字塔进行再次分解，对其顶层取均值，其他层取灰度绝对值大的部分后重构；在Shearlet域高频部分，先计算每个高频子图的信息熵、平均梯度和标准差这三个值，然后分别进行归一化处理，再对处理后的三个值做乘积，最后取乘积大的子图为融合子图；再对融合后的高频子图用非线性算子增强图像细的细节；最后，通过Shearlet逆变换得到最终的融合图像。本方法可以推广到三幅及以上的卫星云图的融合，实现多通道卫星云图融合。 

如图1所示为本发明基于Shearlet变换的多通道卫星云图融合方法的流程示意图。基于Shearlet变换的多通道卫星云图融合方法，包括下列步骤： 

步骤1、对配准后的源图像A和B，图像大小为M×N，分别进行Shearlet变换，分解层数为W，分解方向个数为T（T=2^r,r∈Z^*），得到高频系数SH_A和SH_B、低频系数SL_A和SL_B； 

步骤2、分别对低频系数SL_A和SL_B做拉普拉斯金字塔分解，分解层数为Q，得到分解图像LA和LB，第q（1≤q≤Q）层子图分别为LA_q和LB_q； 

步骤3、对拉普拉斯金字塔顶层子图LA_Q和LB_Q用均值法进行融合，得到融合结果LF_Q为 

$>{\mathrm{LF}}_Q(i,j)=\frac{{\mathrm{LA}}_Q(i,j)+{\mathrm{LB}}_Q(i,j)}{2}$>

其中，1≤i≤CL_Q，1≤j≤RL_Q，CL_Q是分解子图第Q层图像的行数，RL_Q是分解子图第Q层图像的列数； 

步骤4、对拉普拉斯金字塔其他层子图LA_q和LB_q（1≤q≤Q-1）用灰度绝对值取大的融合规则融合，则融合结果LF_q为 

$>{\mathrm{LF}}_q(i,j)=\left(\begin{array}{cc}{\mathrm{LA}}_q(i,j),& \left|{\mathrm{LA}}_q(i,j)\right|\ge \left|{\mathrm{LB}}_q(i,j)\right|\\ {\mathrm{LB}}_q(i,j),& \left|{\mathrm{LA}}_q(i,j)\right|<\left|{\mathrm{LB}}_q(i,j)\right|\end{array}\right)$>

步骤5、将融合后得到的拉普拉斯金字塔LF进行重构，得到低频部分的融合结果TL_F； 

步骤6、在Shearlet变换域的高频部分，分别求每层每个方向子图的信息熵、平均梯度和标准差，记第w（1≤w≤W）层第t（1≤t≤T）个方向的高频子图分别为和其大小和原图大小一样是M×N，则其信息熵E为 

$>E=-\underset{i=0}{\overset{L-1}{\Sigma }}{P}_i{\log }_2{P}_i$>

其中，P_i表示子图中像素灰度值为i的概率,L表示图像中的灰度级的个数；高频子图的平均梯度表示为 

$>\bar{G}=\frac{\underset{\mathrm{ii}=1}{\overset{M-1}{\Sigma }}\underset{\mathrm{jj}=1}{\overset{N-1}{\Sigma }}\sqrt{\frac{{\left(\frac{{\partial \mathrm{SH}}_w^t(x_{\mathrm{ii}},y_{\mathrm{jj}})}{{\partial x}_{\mathrm{ii}}}\right)}^2+{\left(\frac{{\partial \mathrm{SH}}_w^t(x_{\mathrm{ii}},y_{\mathrm{jj}})}{{\partial y}_{\mathrm{jj}}}\right)}^2}{2}}}{(M-1)(N-1)}$>

其中，表示高频子图或中x_ii行y_jj列的像素点，1≤ii≤M,1≤jj≤N。高频子图的标准差σ表示为 

$>\sigma =\sqrt{\underset{\mathrm{ii}=1}{\overset{M}{\Sigma }}\underset{\mathrm{jj}=1}{\overset{N}{\Sigma }}\frac{{({\mathrm{SH}}_w^t(x_{\mathrm{ii}},y_{\mathrm{jj}})-\bar{h})}^2}{M\times N}}$>

其中，表示高频子图的灰度均值； 

步骤7、对高频子图和相对应的信息熵E、平均梯度标准差σ分别进行归一化处理，得到归一化后的信息熵E_g、平均梯度标准差σ_g，选择三者乘积值大的作为融合后的高频子图，即 

$>{\left({\mathrm{SH}}_F\right)}_w^t=\left(\begin{array}{cc}{\left({\mathrm{SH}}_A\right)}_w^t,& {\left(E_g\right)}_A\times {\left({\bar{G}}_g\right)}_A\times {\left({\sigma }_g\right)}_A\ge {\left(E_g\right)}_B\times {\left({\bar{G}}_g\right)}_B\times {\left({\sigma }_g\right)}_B\\ {\left({\mathrm{SH}}_B\right)}_w^t,& {\left(E_g\right)}_A\times {\left({\bar{G}}_g\right)}_A\times {\left({\sigma }_g\right)}_A<{\left(E_g\right)}_B\times {\left({\bar{G}}_g\right)}_B\times {\left({\sigma }_g\right)}_B\end{array}\right)$>

步骤8、对融合后的高频子图进行基于离散平稳小波变换的非线性增强处理，设高频子图的所有像素灰度中绝对值最大的为maxh，则可以得到增强后的高频子图为 

$>{\left({E\_\mathrm{SH}}_F\right)}_w^t(\mathrm{ii},\mathrm{jj})=a·\max h\left\{\mathrm{sigm}\right[c({\mathrm{Sh}}_w^t(\mathrm{ii},\mathrm{jj})-b)]-\mathrm{sigm}[-c({\mathrm{Sh}}_w^t(\mathrm{ii},\mathrm{jj})+b)\left]\right\}$>

其中，b=0.35，c=20；a=1/(d₁-d₂)，d₁=sigm(c×(1+b))，d₂=sigm(-c×(1-b))； 

$>{\mathrm{Sh}}_w^t(\mathrm{ii},\mathrm{jj})={\left({\mathrm{SH}}_f\right)}_w^t(\mathrm{ii},\mathrm{jj})/\max h.$>

步骤9、对融合处理后的Shearlet系数值进行Shearlet逆变换，得到最终的融合图像F。 

本发明还包括：在前面的融合结果上再加上其他通道的云图，实现三幅及以上的云图融合，达到多通道卫星云图融合的效果。 

其中所述的对Shearlet分解系数的处理，是对图像的Shearlet分解系数分别按照各自的融合规则融合。Shearlet域低频部分，利用拉普拉斯金字塔进行再次分解，对其顶层取均值，其他层取灰度绝对值大的部分后重构；Shearlet域高频部分，求得高频子图的信息熵、平均梯度和标准差，取三者乘积大的部分作为融合后的高频子图。 

Shearlet低频系数的融合规则是首先做拉普拉斯金字塔分解，然后对拉普拉斯金字塔顶层子图用均值法进行融合，对拉普拉斯金字塔其它层子图用灰度绝对值取大的融合规则进行融合，最后对融合后的拉普拉斯金字塔进行重构，得到新的Shearlet低频系数。 

下边通知实验例进一步说明本发明的效果。 

实验例1： 

如图3（a）～(h)所示，我们选取来源于2005年8月31日12时00分台风“泰利”的红外2通道和水气通道云图图像作为原图进行融合处理。通过MATLAB7.0将其处理为灰度图像，均是从图2所示的2288×2288大小这种类型的卫星云图中截取512×512像素的融合实验图像。其中每个像素点用数字表示其亮度。数字越大说明该点越亮。 

我们分别对两幅待融合云图图像进行Shearlet变换，所述步骤的分解过程中分解层数均为一层，八个方向。为了验证本文提出的融合算法的有效性，将本文方法的融合结果与拉普拉斯金字塔图像融合方法、经典离散正交小波图像融合方法、Contourlet图像融合方法（融合规则为低频系数取均值，高频系数取区域能量大的部分，分解方向设定为[0,2]）、Curvelet图像融合方法（单独的Curvelet图像融合方法，融合规则低频系数取均值，高频系数取窗口区域内能量大的部分，其中窗口大小为3×3）和NSCT图像融合方法（NSCT结合能量融合的算法，其中NSCT的分解方向设定为[3,3]）这5种方法的融合结果进行对比。其中，拉普拉斯金字塔图像融合方法和经典离散正交小波图像融合方法的融合规则相同，均采用低频部分取均值，高频部分取灰度绝对值较大的部分的方法。 

如图3(a)和图3(b)所示为2005年8月31日12时00分台风“泰利”的红外2通道和水气通道云图（512×512）。图3(c)拉普拉斯金字塔的融合结果，图3(d)是经典离散正交小波的融合结果，图3(e)是Curvelet图像融合方法的融合结果，图3(f)是Contourlet图像融合方法的融合结果，图3(g)是NSCT图像融合方法的融合结果，图3(h)是本文融合算法的融合结果。 

从图3(a)～(h)中可以看到，图3(c)拉普拉斯金字塔融合算法的融合图像比图3(d)经典离散正交小波融合算法的融合结果更清晰些；图3(e)Curvelet图像融合算法的融合图像和图3(g)NSCT图像融合算法的融合图像比较接近于图3(b)水气通道的融合源图，图像灰度值稍偏大，台风眼和周围的云区别较小；图3(f)Contourlet图像融合算法的融合图像有细小的网格现象。图3(h)本文图像融合算法的融合图像和图3(c)拉普拉斯金字塔融合算法的融合图像比较相似，图像清晰，风眼附近信息突出。为了更清晰的对比细节部分，我们截取上述融合结果的部分图像，如图4所示。 

从图4(a)～(f)中可以看到图4(e)Curvelet图像融合算法的融合结果和图4(f)NSCT图像融合算法的融合结果中云团和风眼灰度差别不是很大，相 较于其他组的结果区分度更差些。其他几组的融合结果中台风旋的图像效果比较相近。本文算法的融合结果能有效地突出台风眼信息，台风主体云系整体比较平滑，有利于提高基于卫星云图的台风中心定位的精度。 

为了客观地评价以上图像的融合效果，本发明分别计算了以上融合图像的信息熵E、平均梯度G、标准差σ和平均相关系数Average_Corr，并计算了这四个评价参数的乘积。由于本融合算法是针对多通道卫星云图的，目的是提高台风中心定位的精度，所以十分注重融合图像的信息量、空间分辨率、清晰度等指标，保证融合图像具有较优的细节信息和纹理特征。图像的信息量可以用信息熵来评价，信息熵可以客观地评价图像在融合前后信息量的多少，信息熵E越大表示融合图像所包含的平均信息量有所增加，信息越丰富融合效果就越好。图像的空间分辨率可以用相关系数、标准差来评价。相关系数可以用来衡量两幅图像之间的相关程度。如果融合结果和融合原图的相关系数越接近1，相关程度就越大，也就是融合效果越好。对于融合源图像A和融合图像F，其相关系数为Corr(A,F)；对于融合源图像B和融合图像F，其相关系数为Corr(B,F)；则平均相关系数为 

$>\mathrm{Average}\_\mathrm{Corr}=\frac{\mathrm{Corr}(A,F)+\mathrm{Corr}(B,F)}{2}$>

平均相关系数Average_Corr越接近于1，说明融合结果越好。标准差σ反映了图像灰度值相对于图像灰度平均值的离散程度。此外，标准差σ越大表示融合图像的信息反差越大，越容易体现信息。反之，标准差σ越小，表示图像的灰度级分布比较集中，对比度不明显，不易体现融合图像的细节信息。图像的清晰度可以用平均梯度来评价，平均梯度可敏感地反映图像对微小细节反差表达的能力。一般来说，平均梯度越大表示图像灰度变化率越大，图像越清晰。综合考虑这四个评价参数，本发明提出用这四个参数的乘积综合评价融合图像的效果，乘积值越大，说明融合图像的效果越好，信息量越丰富，图像越清晰，越有利于台风中心定位。 

台风“泰利”的红外2通道和水气通道云图融合结果的性能指标如表1所示。 

表1  图3(c)～(h)中台风“泰利”的红外2通道和水气通道云图的各种融合结果的性能参数比较 

从表1中可得，本文算法融合结果的平均梯度和标准差比其他融合算法的结果都好，四个评价参数的乘积也是最优的，说明本算法的融合效果的综合性能是最优的。而信息熵和平均相关系数不是最优，但跟其他融合算法的结果相差不大，信息熵最大相差为0.008，平均相关系数最大相差为0.003，所以可以认为这两个参数与其他融合方法效果相当。 

如图5（a）～(f)所示，是对图3(c)～(h)中的各种融合结果截图39×39大小的图像，然后用台风中心定位算法定位台风中心。台风中心定位算法是首先定出台风密闭云区域，然后基于密闭云内台风中心区域梯度信息最丰富的特点，用9×9大小窗口遍历密闭云区，选出密闭云区中纹理线交点最多的窗口定位为台风中心区域，再取中心区域的几何中心为台风中心。找到台风中心后，用“+”符号在512×512的融合结果图中标记出中心位置，如图6(a)～(h)所示。 

从图6(c)～(h)中可以看到，图6(c)拉普拉斯金字塔融合结果的台风中心定位结果、图6(d)经典离散正交小波融合结果的台风中心定位结果、图6(f)Contourlet融合结果的台风中心定位结果以及图6(h)本文算法融合结果的台风中心定位结果都比较接近台风中心，细微差别难以用肉眼观察， 所以我们根据台风中心定位的经纬度误差计算台风中心的距离误差，2005年8月31日12时00分台风“泰利”红外2通道和水气通道云图融合结果的台风中心定位误差如表2所示。 

表2  2005年8月31日12时00分台风“泰利”红外2通道和水气通道云图各种融合方法结果的中心定位误差比较 

从表2中可以看到，本文算法台风中心误差为39.37km，中心定位结果误差最小，优于单独红外2通道、水气通道以及其他融合方法的中心定位结果。 

实验例2： 

如图7（a）～（h）所示，我们选取来源于2006年5月11日00时00分台风“珍珠”的红外1通道和水气通道云图图像作为原图进行融合处理。其红外1通道和水气通道云图如图7(a)和图7(b)所示。图7(c)拉普拉斯金字塔的融合结果，图7(d)是经典离散正交小波的融合结果，图7(e)是Curvelet图像融合方法的融合结果图7(f)是Contourlet图像融合方法的融合结果，图7(g)是NSCT图像融合方法的融合结果，图7(h)是本图像融合算法的融合结果。由于图7中是无眼台风，从气旋外围的云图细节来看， 图7(e)Curvelet图像融合算法的融合结果和图7(g)NSCT图像融合算法的融合结果中灰度值大的成分比较多，细节较模糊；其他组融合结果中图7(c)拉普拉斯金字塔融合算法的融合结果稍优于图7(d)经典正交小波融合算法的融合结果，图7(h)本文算法融合算法的融合结果和图7(c)拉普拉斯金字塔融合算法的融合结果比较接近，图7(f)Contourlet融合算法的融合结果也较好。为了更清晰的对比细节部分，我们截取上述融合结果的部分图像，如图8(a)～(f)所示。 

图8中可以看到，这组无眼台风云图的风旋中心云图仍然较亮，对比来看差异不是很明显，几乎都相差不大。从外围云图信息来看，图8(a)拉普拉斯金字塔融合结果、图8(b)经典正交小波融合结果和图8(f)本发明算法融合结果稍优，图8(c)Curvelet融合结果和图8(e)NSCT融合结果优于过亮，云图不是很清晰，不能突出云团信息，图8(d)Contourlet融合结果的细节、边缘部分不够清晰。 

各种融合算法对图7中台风“珍珠”的红外1通道和水气通道云图进行融合，其融合结果的性能指标如表3所示。 

表3  图7（c）～(h)中台风“珍珠”的红外1通道和水气通道云图的各种融合结果的性能参数比较 

从表3中可得，本文算法融合结果的平均梯度和标准差比其他融合算 法的结果都好，四个评价参数的乘积也是最优的，说明本算法的融合效果是最优的。而信息熵和平均相关系数不是最优，但跟其他融合算法的结果相差不大，信息熵最大相差为0.025，平均相关系数最大相差为0.005，所以可以认为这两个参数与其他融合方法的性能基本相当。 

然后对各种方法的融合结果截取39×39大小的图像,如图9（a）～(f)所示，用台风中心定位算法对图像进行台风中心定位，以验证本融合算法的有效性。由于这组台风云图是无眼的，所以图像灰度值比较大，但几幅融合结果的截图看起来差异也不是很大。台风中心定位的结果用“+”符号在512×512的融合结果图中标记中心位置，如图10(a)～(h)所示。图10(a)～(h)中各种融合方法的台风中心定位结果都有所不同，融合源图中图10(a)红外1通道的台风中心定位偏离中心较远，图10(b)水气通道的台风中心定位比较接近中心。图10(c)拉普拉斯金字塔融合结果的台风中心定位结果和图10(g)NSCT融合结果的台风中心定位结果偏下方，图10(d)经典离散正交小波融合结果的台风中心定位结果、图10(e)Curvelet融合结果的台风中心定位结果和图10(f)Contourlet融合结果的台风中心定位结果都稍稍偏右些，比较接近中心。图10(f)本发明融合结果的台风中心定位结果比较靠近中心，效果较优。根据台风中心定位的经纬度误差计算台风中心的距离误差，2006年5月11日00时00分台风“珍珠”红外1通道和水气通道云图融合结果的台风中心定位误差如表4所示。 

表4  2006年5月11日00时00分台风“珍珠”的红外1通道和水气通道云图各种融合方法结果的中心定位误差比较 

从表4中可以看出，本文算法台风中心误差为76.21km，明显优于单独红外1通道以及其他融合方法的中心定位结果，效果最优。 

实验例3： 

为了进一步说明本文所提出融合算法的有效性，下面分析一下本发明所提出方法的计算复杂度。本文图像融合算法是在MatLab R2009a软件运行的，软件运行在处理器为英特尔酷睿2四核Q94002.66GHz，内存为2GB（金士顿DDR31333MHz），操作系统为Windows XP专业版32位SP3（DirectX9.0c）的戴尔OptiPlex780台式电脑上。此处对各类融合方法的运行时间进行度量，用第二组实验图像做测试，各种融合算法的运行时间如表5所示。 

表5  各种融合算法的运行时间 

从表5中可以看出，除了基于拉普拉斯金字塔的图像融合算法和经典正交离散小波的图像融合算法运行时间比较短外，本发明提出的图像融合算法其他图像融合算法的所用时间都少。因此，本发明提出的融合算法的计算复杂度低，且能得到较好的融合效果。 

通过上述三组实验说明本发明算法能够很好地实现两两卫星云图的融合，通过与拉普拉斯金字塔图像融合方法、经典离散正交小波图像融合方法、Contourlet图像融合方法、Curvelet图像融合方法和NSCT图像融合方法这5种方法的融合结果进行对比，证明本文算法具有最优的标准差和平均梯度，与其他方向相当的信息熵和平均相关系数，且综合评价指标最好，融合图像视觉效果好，能清晰地保留台风眼和云系细节信息，且利用融合结果进行台风中心定位的精度较高，适用于有眼和无眼台风，其卫星云图融合结果的综合效果是最好的，按本方法进行三幅及以上的卫星云图的融合，实现多通道卫星云图融合，有利于结合云图的更多信息，提高台风中心定位的精度。