CN103700075B - 基于Tetrolet变换的多通道卫星云图融合方法 - Google Patents

Abstract

本发明的目的是一种基于Tetrolet变换的多通道卫星云图融合方法。首先对待融合的多通道卫星云图分别进行图像直方图均衡化处理，再分别进行Tetrolet变换，得到低频系数、高频系数及其相对应的拼板覆盖值。然后，在Tetrolet域低频部分，利用拉普拉斯金字塔进行再次分解，对其顶层取均值，其他层取灰度绝对值大的部分后重构；在Tetrolet域高频部分，取每个图像拼板中高频系数标准差较大的部分，拼板覆盖值也取其相对应的值。最后，通过Tetrolet逆变换得到最终的融合图像。实验结果表明本文方法能够很好地实现多通道卫星云图融合，融合图像视觉效果好，能清晰地保留台风眼和云系细节信息，且利用融合结果的台风中心定位精度高，适用于有眼和无眼台风。

Description

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

技术领域

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

背景技术

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

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

A.Abd-Elrahman等人提出了一种在利用小波变换进行卫星云图融合时提高其云相关阴影区域的增强方法，并保留了细节信息，其结果有效地提高了云图质量。Lee,Y.等人提出一种新的小波域卫星图像融合算法，考虑了每一幅源图像的强度和光谱范围以及方面的相关光谱响应。用高斯函数之和表示每个通道的光谱响应，然后用高斯函数建模调整图像的空间和光谱分辨率。该方法融合结果的PSNR（Peak Signal-Noise Ratio）值、均方根误差和相关系数比传统方法的效果更优。V.Harikumar等人提出一种新的基于压缩感知和图像切割法的多分辨率图像融合方法，并对MS（multispectral）和PAN（panchromatic）图像进行融合。该方法用基于正则化的方法解决MS图像不适定的反转问题，用缩短的二次平滑法保留融合图像的间断点，对卫星图像的实验表明该方法优于其他几种同类图像融合方法。N.H.Kaplan等人针对多光谱和全色卫星图像提出了一种新的基于双边金字塔图像融合方法，MS和PAN图像用双边金字塔分解，按照替换法和附加法得到融合图像，与比较广泛使用的HIS（Intensity-Hue-Saturation）变换、àtrous小波变换（àtrous Wavelet Transform，ATWT）融合方法对比，其提出方法的融合结果的评价参数较优。姜波利用自适应的PCNN（Pulse CoupledNeural Network）和NSCT（NonSubsampled Contourlet Transform）提出一种将红外云图与可见光云图进行融合的图像融合方法，该方法既能获得较好的方向信息，去除频谱混叠效应，又具有很好的抗噪性。Amr M.Ragheb等人提出一种由IHS变换和离散小波变换（Discrete Wavelet Frame Transform，DWFT）结合的融合方法对卫星图像进行融合，这种方法是为了在图像不变形的情况下得到尽可能高的光谱和空间分辨率，其实验结果表明该方法优于传统融合方法。从上述来看，目前对多通道卫星云图的融合方法主要有基于小波变换的图像融合方法、基于压缩感知的图像融合方法、基于金字塔分解的图像融合方法、基于多尺度变换的图像融合方法、基于IHS变换的对全彩色图像进行融合处理的方法等，并基于这些方法进行一些图像细节的优化，获得较优的融合结果。

多源图像融合经过几十年的发展，其技术的理论依据越来越完备，应用场合也越来越广泛。目前，多源图像融合方法主要是基于空间域的图像融合方法和基于变换域的图像融合方法。基于空间域的图像融合方法主要有线性加权图像融合、假彩色图像融合、基于调制的图像融合、基于统计的图像融合以及基于神经网络的图像融合等。基于变换域的融合方法包括基于离散余弦变换（Discrete Cosine Transform,DCT）的图像融合算法、基于快速傅里叶变换（FastFourier Transformation,FFT）的图像融合算法、基于多尺度分解的图像融合算法等。其中，多尺度分解的图像融合算法应用广泛且得到了越来越多研究者的重视。多尺度分解的图像融合方法始于Burt P.J.在1984年提出的基于拉普拉斯金字塔（Laplacian Pyramid）分解的图像融合算法。随后，Ranchin和Wald在1993年提出了基于离散小波变换的图像融合算法应用于遥感图像。随着多尺度几何分析的发展，研究者们对基于多尺度分解的图像融合技术进行了不断地探索。

Miao Qiguang等人提出了一种基于Contourlet变换的图像融合方法，对高频部分进行区域能量比较，并结合一致性检验，得到在边缘保持和纹理信息方面优于小波变换方法和拉普拉斯金字塔的图像融合方法的结果。Shutao Li等人提出了一种结合Curvelet和小波变换的多聚焦图像融合算法，该图像融合结果优于其他任何单独的多尺度融合方法。Juan Lu等人提出了一种基于NSCT和能量熵的图像融合算法，该算法的融合结果具有更丰富的方向信息和很强的噪声鲁棒性。Qi-guang Miao等人利用Shearlet变换的方向性、定位性、各向异性、多尺度等优势，将Shearlet变换用于图像融合，其融合结果比其他方法包含更多的细节和较小的失真信息。

Jens Krommweh于2009年提出了Tetrolet变换，它兼顾了Curvelet等方法的优点，是一种新的稀疏图像表示的自适应Haar小波变换，能很好地表示图像几何结构特性，稀疏化程度非常高。至今，Tetrolet变换已应用于图像去噪、图像稀疏表示和图像复原领域，在图像融合方向还没有相关的研究成果，本发明尝试将Tetrolet变换引入到图像融合领域，并将其结合拉普拉斯金字塔分解的融合方法实现多通道卫星云图融合。

发明内容

本发明的目的是提供一种基于Tetrolet变换的多通道卫星云图融合方法。先对待融合的多通道卫星云图分别进行图像直方图均衡化处理，再分别进行Tetrolet变换，得到低频系数、高频系数及其相对应的拼板覆盖值。然后，在Tetrolet域低频部分，利用拉普拉斯金字塔进行再次分解，对其顶层取均值，其他层取灰度绝对值大的部分后重构；在Tetrolet域高频部分，取每个图像拼板中高频系数标准差较大的部分，拼板覆盖值也取其相对应的值。最后，通过Tetrolet逆变换得到最终的融合图像。

为了达到所述目的，本发明基于Tetrolet变换的多通道卫星云图融合方法，具体步骤如下：

步骤1对配准后的源图像A和B分别进行直方图均衡化处理，得到图像A'和B'；

步骤2对图像A'和B'进行Tetrolet变换，分解层数为M，得到高频系数TH_A和TH_B、低频系数TL_A和TL_B及相应的拼板覆盖值TC_A和TC_B；

步骤3分别对低频系数TL_A和TL_B做拉普拉斯金字塔分解，分解层数为N，得到分解图像LA和LB，第n（0＜n≤N）层子图分别为LA_n和LB_n；

步骤4对拉普拉斯金字塔顶层子图LA_N和LB_N用均值法进行融合，得到融合结果LF_N为

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

其中，0＜i≤CL_N，0＜j≤RL_N，CL_N是分解子图第N层图像的行数，RL_N是分解子图第N层图像的列数；

步骤5对拉普拉斯金字塔其他层子图LA_n和LB_n用灰度绝对值取大的融合规则融合，则融合结果LF_n为

${\mathrm{LF}}_n(i,j)=\left\{\begin{array}{cc}{\mathrm{LA}}_n(i,j),& \left|{\mathrm{LA}}_n(i,j)\right|\&GreaterEqual;\left|{\mathrm{LB}}_n(i,j)\right|\\ {\mathrm{LB}}_n(i,j),& \left|{\mathrm{LA}}_n(i,j)\right|<\left|{\mathrm{LB}}_n(i,j)\right|\end{array}\right.;$

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

步骤7在Tetrolet变换域高频系数部分，对每个分解子块按照标准差取大的融合规则融合，记第m（0＜m≤M）层p行q列的高频系数分别为(TH_A)_m(p,q)和(TH_B)_m(p,q)，其大小为12×1的矩阵，其中0＜p≤CH_m，0＜q≤RH_m，CH_m是分解子图第m层高频系数的行数，RH_m是分解子图第m层高频系数的列数，则

${\left({\mathrm{TH}}_F\right)}_m(p,q)=\left\{\begin{array}{cc}{\left({\mathrm{TH}}_A\right)}_m(p,q),& \mathrm{std}\left({\left({\mathrm{TH}}_A\right)}_m(p,q)\right)\&GreaterEqual;\mathrm{std}\left({\left({\mathrm{TH}}_B\right)}_m(p,q)\right)\\ {\left({\mathrm{TH}}_B\right)}_m(p,q),& \mathrm{std}\left({\left({\mathrm{TH}}_A\right)}_m(p,q)\right)<\mathrm{std}\left({\left({\mathrm{TH}}_B\right)}_m(p,q)\right)\end{array}\right.$

其中std是标准差，表示为

$\mathrm{std}=\sqrt{\frac{\underset{i=1}{\overset{K}{\mathrm{\&Sigma;}}}\underset{j=1}{\overset{G}{\mathrm{\&Sigma;}}}{(x(i,j)-\stackrel{\&OverBar;}{x})}^2}{K\&times;G-1}}$

这里，是x的均值，1≤i≤K，1≤j≤G。在Tetrolet域的高频系数矩阵中K＝12，G＝1。

步骤8Tetrolet变换中的拼板覆盖值的选取根据高频部分的选择而定，即

${\left({\mathrm{TC}}_F\right)}_m(p,q)=\left\{\begin{array}{cc}{\left({\mathrm{TC}}_A\right)}_m(p,q),& {\mathrm{TH}}_m(p,q)={\left({\mathrm{TH}}_A\right)}_m(p,q)\\ {\left({\mathrm{TC}}_B\right)}_m(p,q),& {\mathrm{TH}}_m(p,q)={\left({\mathrm{TH}}_B\right)}_m(p,q)\end{array}\right.$

其中，(TC_A)_m(p,q)和(TC_B)_m(p,q)分别表示图像A和B经Tetrolet变换后第m层p行q列的拼板覆盖值。

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

优选的，所述的步骤2中基于Tetrolet变换的多通道卫星云图融合方法，其主要步骤包括：

步骤1把图像分成4×4子块；

步骤2考虑117种四格拼板分割方案，应用Haar小波变换获得每一种方案的高频和低频系数，选择Tetrolet系数小的方案，得到每个子块中的稀疏Tetrolet表示；

步骤3重新排列每个子块中的高频和低频部分的系数，把它们排列成2×2的子块；

步骤4存储Tetrolet系数的高频部分；

步骤5对低频部分重复进行步骤1-4的操作。即对图像A'和B'进行Tetrolet变换，分解层数为M，把图像分成4×4大小的子块，根据每个子块的几何特征，用117种四格拼板的分割方案分割子块，得到每个子块的稀疏Tetrolet表示，即高频系数TH_A和TH_B、低频系数TL_A和TL_B及相应的拼板覆盖值TC_A和TC_B。

优选的，所述的步骤3-6中对融合源图的Tetrolet低频系数按照融合规则进行处理。首先对低频系数做拉普拉斯金字塔分解，然后对拉普拉斯金字塔顶层子图用均值法进行融合，Tetrolet低频系数的融合步骤为：

步骤1别对低频系数TL_A和TL_B做拉普拉斯金字塔分解，分解层数为N，得到分解图像LA和LB，第n（0＜n≤N）层子图分别为LA_n和LB_n；

步骤2对拉普拉斯金字塔顶层子图LA_N和LB_N用均值法进行融合，得到融合结果LF_N为

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

其中，0＜i≤CL_N，0＜j≤RL_N，CL_N是分解子图第N层图像的行数，RL_N是分解子图第N层图像的列数；

步骤3对拉普拉斯金字塔其他层子图LA_n和LB_n用灰度绝对值取大的融合规则融合，则融合结果LF_n为

${\mathrm{LF}}_n(i,j)=\left\{\begin{array}{cc}{\mathrm{LA}}_n(i,j),& \left|{\mathrm{LA}}_n(i,j)\right|\&GreaterEqual;\left|{\mathrm{LB}}_n(i,j)\right|\\ {\mathrm{LB}}_n(i,j),& \left|{\mathrm{LA}}_n(i,j)\right|<\left|{\mathrm{LB}}_n(i,j)\right|\end{array}\right.;$

步骤4将融合后得到的拉普拉斯金字塔LF进行重构，得到低频部分的融合结果TL_F。对拉普拉斯金字塔其它层子图用灰度绝对值取大的融合规则进行融合，最后对融合后的拉普拉斯金字塔进行重构，得到新的Tetrolet低频系数。

优选的，所述的步骤7中对融合源图的Tetrolet高频系数按照融合规则进行处理。在Tetrolet变换域高频系数部分，对每个分解子块按照标准差取大的融合规则融合，记第m（0＜m≤M）层p行q列的高频系数分别为(TH_A)_m(p,q)和(TH_B)_m(p,q)，其大小为12×1的矩阵，其中0＜p≤CH_m，0＜q≤RH_m，CH_m是分解子图第m层高频系数的行数，RH_m是分解子图第m层高频系数的列数，则

${\left({\mathrm{TH}}_F\right)}_m(p,q)=\left\{\begin{array}{cc}{\left({\mathrm{TH}}_A\right)}_m(p,q),& \mathrm{std}\left({\left({\mathrm{TH}}_A\right)}_m(p,q)\right)\&GreaterEqual;\mathrm{std}\left({\left({\mathrm{TH}}_B\right)}_m(p,q)\right)\\ {\left({\mathrm{TH}}_B\right)}_m(p,q),& \mathrm{std}\left({\left({\mathrm{TH}}_A\right)}_m(p,q)\right)<\mathrm{std}\left({\left({\mathrm{TH}}_B\right)}_m(p,q)\right)\end{array}\right.$

其中std是标准差，表示为

$\mathrm{std}=\sqrt{\frac{\underset{i=1}{\overset{K}{\mathrm{\&Sigma;}}}\underset{j=1}{\overset{G}{\mathrm{\&Sigma;}}}{(x(i,j)-\stackrel{\&OverBar;}{x})}^2}{K\&times;G-1}}$

这里，是x的均值，1≤i≤K，1≤j≤G。在Tetrolet域的高频系数矩阵中K＝12，G＝1。

优选的，所述的步骤8中对融合源图的Tetrolet拼板覆盖值按照融合规则进行处理。Tetrolet变换中的拼板覆盖值的选取根据高频部分的选择而定，即

${\left({\mathrm{TC}}_F\right)}_m(p,q)=\left\{\begin{array}{cc}{\left({\mathrm{TC}}_A\right)}_m(p,q),& {\mathrm{TH}}_m(p,q)={\left({\mathrm{TH}}_A\right)}_m(p,q)\\ {\left({\mathrm{TC}}_B\right)}_m(p,q),& {\mathrm{TH}}_m(p,q)={\left({\mathrm{TH}}_B\right)}_m(p,q)\end{array}\right.$

其中，(TC_A)_m(p,q)和(TC_B)_m(p,q)分别表示图像A和B经Tetrolet变换后第m层p行q列的拼板覆盖值。

优选的，所述的步骤9中对融合处理后的Tetrolet系数值（即低频系数TL_F、高频系数TH_F和拼板覆盖值TC_F）进行Tetrolet逆变换，得到最终的融合图像F。

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

附图说明

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

图2为我国气象卫星FY-2C返回的5通道卫星云图；

图3为2007年9月17日12时00分0713号台风“韦帕”的红外2通道和水气通道云图（有眼台风）的多通道卫星云图融合实验结果；

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

图5为图3中“韦帕”云图的融合结果截取云图，便于进行台风中心定位；

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

图7为2008年9月19日06时00分0814号台风“黑格比”的红外1通道和水气通道云图（无眼台风）的多通道卫星云图融合实验结果；

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

图9为图7中“黑格比”云图的融合结果截取云图，；

图10为利用图7中“黑格比”云图的融合结果图像的台风中心定位结果示意图。

具体实施方式

本发明提出一种基于Tetrolet变换的多通道卫星云图融合方法。首先对待融合的多通道卫星云图分别进行图像直方图均衡化处理，再分别进行Tetrolet变换，得到低频系数、高频系数及其相对应的拼板覆盖值。然后，在Tetrolet域低频部分，利用拉普拉斯金字塔进行再次分解，对其顶层取均值，其他层取灰度绝对值大的部分后重构；在Tetrolet域高频部分，取每个图像拼板中高频系数标准差较大的部分，拼板覆盖值也取其相对应的值。最后，通过Tetrolet逆变换得到最终的融合图像。

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

步骤1对配准后的源图像A和B分别进行直方图均衡化处理，得到图像A'和B'；

步骤2对图像A'和B'进行Tetrolet变换，分解层数为M，得到高频系数TH_A和TH_B、低频系数TL_A和TL_B及相应的拼板覆盖值TC_A和TC_B；

步骤3分别对低频系数TL_A和TL_B做拉普拉斯金字塔分解，分解层数为N，得到分解图像LA和LB，第n（0＜n≤N）层子图分别为LA_n和LB_n；

步骤4对拉普拉斯金字塔顶层子图LA_N和LB_N用均值法进行融合，得到融合结果LF_N为

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

其中，0＜i≤CL_N，0＜j≤RL_N，CL_N是分解子图第N层图像的行数，RL_N是分解子图第N层图像的列数；

步骤5对拉普拉斯金字塔其他层子图LA_n和LB_n用灰度绝对值取大的融合规则融合，则融合结果LF_n为

${\mathrm{LF}}_n(i,j)=\left\{\begin{array}{cc}{\mathrm{LA}}_n(i,j),& \left|{\mathrm{LA}}_n(i,j)\right|\&GreaterEqual;\left|{\mathrm{LB}}_n(i,j)\right|\\ {\mathrm{LB}}_n(i,j),& \left|{\mathrm{LA}}_n(i,j)\right|<\left|{\mathrm{LB}}_n(i,j)\right|\end{array}\right.;$

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

步骤7在Tetrolet变换域高频系数部分，对每个分解子块按照标准差取大的融合规则融合，记第m（0＜m≤M）层p行q列的高频系数分别为(TH_A)_m(p,q)和(TH_B)_m(p,q)，其大小为12×1的矩阵，其中0＜p≤CH_m，0＜q≤RH_m，CH_m是分解子图第m层高频系数的行数，RH_m是分解子图第m层高频系数的列数，则

${\left({\mathrm{TH}}_F\right)}_m(p,q)=\left\{\begin{array}{cc}{\left({\mathrm{TH}}_A\right)}_m(p,q),& \mathrm{std}\left({\left({\mathrm{TH}}_A\right)}_m(p,q)\right)\&GreaterEqual;\mathrm{std}\left({\left({\mathrm{TH}}_B\right)}_m(p,q)\right)\\ {\left({\mathrm{TH}}_B\right)}_m(p,q),& \mathrm{std}\left({\left({\mathrm{TH}}_A\right)}_m(p,q)\right)<\mathrm{std}\left({\left({\mathrm{TH}}_B\right)}_m(p,q)\right)\end{array}\right.$

其中std是标准差，表示为

$\mathrm{std}=\sqrt{\frac{\underset{i=1}{\overset{K}{\mathrm{\&Sigma;}}}\underset{j=1}{\overset{G}{\mathrm{\&Sigma;}}}{(x(i,j)-\stackrel{\&OverBar;}{x})}^2}{K\&times;G-1}}$

这里，是x的均值，1≤i≤K，1≤j≤G。在Tetrolet域的高频系数矩阵中K＝12，G＝1。

步骤8Tetrolet变换中的拼板覆盖值的选取根据高频部分的选择而定，即

${\left({\mathrm{TC}}_F\right)}_m(p,q)=\left\{\begin{array}{cc}{\left({\mathrm{TC}}_A\right)}_m(p,q),& {\mathrm{TH}}_m(p,q)={\left({\mathrm{TH}}_A\right)}_m(p,q)\\ {\left({\mathrm{TC}}_B\right)}_m(p,q),& {\mathrm{TH}}_m(p,q)={\left({\mathrm{TH}}_B\right)}_m(p,q)\end{array}\right.$

其中，(TC_A)_m(p,q)和(TC_B)_m(p,q)分别表示图像A和B经Tetrolet变换后第m层p行q列的拼板覆盖值。

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

其中所述的基于Tetrolet变换的图像分解算法步骤包括：

步骤1把图像分成4×4子块；

步骤2考虑117种四格拼板分割方案，应用Haar小波变换获得每一种方案的高频和低频系数，选择Tetrolet系数小的方案，得到每个子块中的稀疏Tetrolet表示；

步骤3重新排列每个子块中的高频和低频部分的系数，把它们排列成2×2的子块；

步骤4存储Tetrolet系数的高频部分；

步骤5对低频部分重复进行步骤1-4的操作。

其中所述的对Tetrolet分解系数的处理，是对图像的Tetrolet分解系数分别按照各自的融合规则融合。Tetrolet域低频部分，利用拉普拉斯金字塔进行再次分解，对其顶层取均值，其他层取灰度绝对值大的部分后重构；Tetrolet域高频部分，取每个图像拼板中高频系数标准差较大的部分，拼板覆盖值也取其相对应的值。

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

步骤1别对低频系数TL_A和TL_B做拉普拉斯金字塔分解，分解层数为N，得到分解图像LA和LB，第n（0＜n≤N）层子图分别为LA_n和LB_n；

步骤2对拉普拉斯金字塔顶层子图LA_N和LB_N用均值法进行融合，得到融合结果LF_N为

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

其中，0＜i≤CL_N，0＜j≤RL_N，CL_N是分解子图第N层图像的行数，RL_N是分解子图第N层图像的列数；

步骤3对拉普拉斯金字塔其他层子图LA_n和LB_n用灰度绝对值取大的融合规则融合，则融合结果LF_n为

${\mathrm{LF}}_n(i,j)=\left\{\begin{array}{cc}{\mathrm{LA}}_n(i,j),& \left|{\mathrm{LA}}_n(i,j)\right|\&GreaterEqual;\left|{\mathrm{LB}}_n(i,j)\right|\\ {\mathrm{LB}}_n(i,j),& \left|{\mathrm{LA}}_n(i,j)\right|<\left|{\mathrm{LB}}_n(i,j)\right|\end{array}\right.;$

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

对于Tetrolet变换域高频系数部分，对每个分解子块按照标准差取大的融合规则融合，记第m（0＜m≤M）层p行q列的高频系数分别为(TH_A)_m(p,q)和(TH_B)_m(p,q)，其大小为12×1的矩阵，其中0＜p≤CH_m，0＜q≤RH_m，CH_m是分解子图第m层高频系数的行数，RH_m是分解子图第m层高频系数的列数，则

${\left({\mathrm{TH}}_F\right)}_m(p,q)=\left\{\begin{array}{cc}{\left({\mathrm{TH}}_A\right)}_m(p,q),& \mathrm{std}\left({\left({\mathrm{TH}}_A\right)}_m(p,q)\right)\&GreaterEqual;\mathrm{std}\left({\left({\mathrm{TH}}_B\right)}_m(p,q)\right)\\ {\left({\mathrm{TH}}_B\right)}_m(p,q),& \mathrm{std}\left({\left({\mathrm{TH}}_A\right)}_m(p,q)\right)<\mathrm{std}\left({\left({\mathrm{TH}}_B\right)}_m(p,q)\right)\end{array}\right.$

其中std是标准差，表示为

$\mathrm{std}=\sqrt{\frac{\underset{i=1}{\overset{K}{\mathrm{\&Sigma;}}}\underset{j=1}{\overset{G}{\mathrm{\&Sigma;}}}{(x(i,j)-\stackrel{\&OverBar;}{x})}^2}{K\&times;G-1}}$

这里，是x的均值，1≤i≤K，1≤j≤G。在Tetrolet域的高频系数矩阵中K＝12，G＝1。

Tetrolet变换中的拼板覆盖值的选取根据高频部分的选择而定，即

${\left({\mathrm{TC}}_F\right)}_m(p,q)=\left\{\begin{array}{cc}{\left({\mathrm{TC}}_A\right)}_m(p,q),& {\mathrm{TH}}_m(p,q)={\left({\mathrm{TH}}_A\right)}_m(p,q)\\ {\left({\mathrm{TC}}_B\right)}_m(p,q),& {\mathrm{TH}}_m(p,q)={\left({\mathrm{TH}}_B\right)}_m(p,q)\end{array}\right.$

其中，(TC_A)_m(p,q)和(TC_B)_m(p,q)分别表示图像A和B经Tetrolet变换后第m层p行q列的拼板覆盖值。

实施例1：

如图3所示，我们选取来源于2007年台风“韦帕”的红外2通道和水气通道云图图像作为原图进行融合处理。通过MATLAB7.0将其处理为灰度图像，均是从图2所示的2288×2288大小这种类型的卫星云图中截取512×512像素的融合实验图像。其中每个像素点用数字表示其亮度。数字越大说明该点越亮，则越显示白色。

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

如图3(a)和图3(b)所示为2005年8月30日09时00分台风“泰利”的红外2通道和水气通道云图（512×512）。对图3(a)和图3(b)分别进行直方图均衡化处理，得到图3(c)和图3(d)。图3(e)拉普拉斯金字塔的融合结果，图3(f)是经典离散正交小波的融合结果，图3(g)是Contourlet图像融合方法的融合结果，图3(h)是Curvelet图像融合方法的融合结果，图3(i)是NSCT图像融合方法的融合结果，图3(j)是Shearlet图像融合方法的融合结果，图3(k)是本发明融合算法的融合结果。

从图3中可以看到，图3(e)拉普拉斯金字塔融合算法的融合图像和图3(f)经典离散正交小波融合算法的融合结果比较相近，图3(g)Contourlet图像融合算法的融合图像有细小的网格现象，图3(h)Curvelet图像融合算法的融合图像比较接近于图3(d)水气通道的融合源图，图像灰度值稍偏大，台风眼和周围的云区别较小。图3(j)Shearlet图像融合算法的融合图像相较于图3(i)NSCT图像融合算法的融合图像有些模糊，细节不够突出。图3(k)本发明算法的融合图像在风眼外围的云细节上稍弱于图3(i)NSCT图像融合算法的融合图像，部分信息轮廓较粗，但主要的信息特点都有所突出。为了更清晰的对比细节部分，我们截取上述融合结果的部分图像，如图4所示。

从图4中可以看到图4(d)Curvelet图像融合算法的融合结果和图4(e)NSCT图像融合算法的融合结果中台风旋过亮，边缘信息有所丢失。其他几组的融合结果中台风旋的图像效果比较相近。本发明算法的融合结果能有效地突出台风眼信息，台风主体云系整体比较平滑，有利于提高基于卫星云图的台风中心定位的精度。

为了客观地评价以上图像的融合效果，本发明分别计算了以上融合图像的信息熵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}$

平均相关系数越接近于1，说明融合结果越好。标准差反映了图像灰度值相对于图像灰度平均值的离散程度。标准差越大，表示融合图像的信息反差越大，越容易体现信息。反之，标准差越小，表示图像的灰度级分布比较集中，对比度不明显，不易体现融合图像的细节信息。联合熵也是信息论里的一个重要概念，可以作为两幅图像相关性的量度，表示两幅图像之间综合信息量的大小。对于融合源图像A和融合图像F，两者的联合熵越大，说明融合效果越好。融合源图像A、B和融合图像F三者的联合熵为

$\mathrm{Joint}\_E(F,A,B)=-\underset{k=1}{\overset{L}{\mathrm{\&Sigma;}}}\underset{i=1}{\overset{L}{\mathrm{\&Sigma;}}}\underset{j=1}{\overset{L}{\mathrm{\&Sigma;}}}{p}_{\mathrm{FAB}}(k,i,j){\log }_2\left[p_{\mathrm{FAB}}(k,i,j)\right]$

其中，p_FAB(k,i,j)表示图像A、B和F的归一化联合直方图。

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

表1图3中台风“韦帕”的红外2通道和水气通道云图的各种融合结果的性能参数比较

从表1中可得，本发明算法融合结果的信息熵、标准差和联合熵比其他融合算法的结果都好，仅仅有平均相关系数不是最优，但跟其他融合算法的结果相差不大，最大相差为0.002，所以认为效果相当。

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

从图6中可以看到，各种融合方法的台风中心定位结果比较接近，细微差别难以用肉眼观察，所以我们根据台风中心定位的经纬度误差计算台风中心的距离误差，2007年9月17日12时00分“韦帕”红外2通道和水气通道云图融合结果的台风中心定位误差如表2所示。

表22007年9月17日12时00分“韦帕”红外2通道和水气通道云图各种融合方法结果的中心定位误差比较

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

实施例2：

如图7所示，我们选取来源于2008年9月19日06时00分台风“黑格比”的红外1通道和水气通道云图图像作为原图进行融合处理。其红外1通道和水气通道云图如图7(a)和图7(b)所示。对图7(a)和图7(b)分别进行直方图均衡化处理，得到图7(c)和图7(d)。图7(e)拉普拉斯金字塔的融合结果，图7(f)是经典离散正交小波的融合结果，图7(g)是Contourlet图像融合方法的融合结果，图7(h)是Curvelet图像融合方法的融合结果，图7(i)是NSCT图像融合方法的融合结果，图7(j)是Shearlet图像融合方法的融合结果，图7(k)是本发明融合算法的融合结果。由于图7中是无眼台风，从风旋外围的云图细节来看，图7(h)Curvelet图像融合算法的灰度值大的成分比较多，细节较模糊；其他组融合结果中图7(e)拉普拉斯金字塔融合结果和图7(f)经典正交小波融合结果比较接近，图7(i)NSCT融合结果和图7(k)本发明算法融合结果比较相似。为了更清晰的对比细节部分，我们截取上述融合结果的部分图像，如图8所示。

图8中可以看到，这组无眼台风云图的风旋中心云图仍然较亮，对比来看差异不是很明显，几乎都相差不大。从外围云图信息来看，图8(a)拉普拉斯金字塔融合结果和图8(b)经典正交小波融合结果稍优，其他融合结果的效果相当。

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

表3图7中台风“黑格比”的红外1通道和水气通道云图的各种融合结果的性能参数比较

从表3中可得，本发明算法融合结果的信息熵和Contourlet融合结果的信息熵一样，优于其他融合结果的信息熵值；本发明算法融合结果的标准差和联合熵明显优于其他融合结果；平均相关系数跟其他融合算法的结果相差不大，最大相差为0.003，所以效果几乎相当。

然后对各种方法的融合结果截取39×39大小的图像（如图9所示），用台风中心定位算法对图像进行台风中心定位，以验证本融合算法的有效性。由于这组台风云图是无眼的，所以图像灰度值比较大，但几幅融合结果的截图看起来差异也不是很大。台风中心定位的结果用“+”符号在512×512的融合结果图中标记中心位置，如图10所示。图10中各种融合方法的台风中心定位结果都有所不同，融合源图中图10(a)红外1通道的台风中心定位比较接近中心，图10(b)水气通道的台风中心定位偏离中心较远。图10(c)～(h)融合结果的台风中心定位都有点稍稍偏左侧，图10(i)本发明算法融合结果的台风中心定位比较靠近中心。根据台风中心定位的经纬度误差计算台风中心的距离误差，2008年9月19日06时00分台风“黑格比”红外1通道和水气通道云图融合结果的台风中心定位误差如表4所示。

表42008年9月19日06时00分台风“黑格比”的红外1通道和水气通道云图各种融合方法结果的中心定位误差比较

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

实施例3：

为了进一步说明本发明所提出融合算法的有效性，下面分析一下本发明所提出方法的计算复杂度。此处对各类融合方法的运行时间进行度量，用第二组实验图像做测试，各种融合算法的运行时间如表5所示。

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

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

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

Claims (2)

1.基于Tetrolet变换的多通道卫星云图融合方法，其特征在于：首先对待融合的多通道卫星云图分别进行图像直方图均衡化处理，再分别进行Tetrolet变换，得到低频系数、高频系数及其相对应的拼板覆盖值；然后，在Tetrolet域低频部分，利用拉普拉斯金字塔进行再次分解，对其顶层取均值，其他层取灰度绝对值大的部分后重构；在Tetrolet域高频部分，取每个图像拼板中高频系数标准差较大的部分，拼板覆盖值也取其相对应的值，最后，通过Tetrolet逆变换得到最终的融合图像,

步骤1对配准后的源图像A和B分别进行直方图均衡化处理，得到图像A'和B'；

步骤2对图像A'和B'进行Tetrolet变换，分解层数为M，得到高频系数TH_A和TH_B、低频系数TL_A和TL_B及相应的拼板覆盖值TC_A和TC_B；

步骤3分别对低频系数TL_A和TL_B做拉普拉斯金字塔分解，分解层数为N，得到分解图像LA和LB，第n层子图分别为LA_n和LB_n，其中0＜n≤N；

步骤4对拉普拉斯金字塔顶层子图LA_N和LB_N用均值法进行融合，得到融合结果LF_N为

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

其中，0＜i≤CL_N，0＜j≤RL_N，CL_N是分解子图第N层图像的行数，RL_N是分解子图第N层图像的列数；

步骤5对拉普拉斯金字塔其他层子图LA_n和LB_n用灰度绝对值取大的融合规则融合，则融合结果LF_n为

${\mathrm{LF}}_n\left(i,j\right)=\left\{\begin{array}{cc}{\mathrm{LA}}_n\left(i,j\right),& \left|{\mathrm{LA}}_n\left(i,j\right)\right|\&GreaterEqual;\left|{\mathrm{LB}}_n\left(i,j\right)\right|\\ {\mathrm{LB}}_n\left(i,j\right),& \left|{\mathrm{LA}}_n\left(i,j\right)\right|<\left|{\mathrm{LB}}_n\left(i,j\right)\right|\end{array}\right.;$

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

步骤7在Tetrolet变换域高频系数部分，对每个分解子块按照标准差取大的融合规则融合，记第m层p行q列的高频系数分别为(TH_A)_m(p,q)和(TH_B)_m(p,q)，其中0＜m≤M，其大小为12×1的矩阵，其中0＜p≤CH_m，0＜q≤RH_m，CH_m是分解子图第m层高频系数的行数，RH_m是分解子图第m层高频系数的列数，则

${\left({\mathrm{TH}}_F\right)}_m\left(p,q\right)=\left\{\begin{array}{cc}{\left({\mathrm{TH}}_A\right)}_m\left(p,q\right),& std\left({\left({\mathrm{TH}}_A\right)}_m\left(p,q\right)\right)\&GreaterEqual;std\left({\left({\mathrm{TH}}_B\right)}_m\left(p,q\right)\right)\\ {\left({\mathrm{TH}}_B\right)}_m\left(p,q\right),& std\left({\left({\mathrm{TH}}_A\right)}_m\left(p,q\right)\right)<std\left({\left({\mathrm{TH}}_B\right)}_m\left(p,q\right)\right)\end{array}\right.$

其中std是标准差，表示为

$std=\sqrt{\frac{\underset{i=1}{\overset{K}{\&Sigma;}}\underset{j=1}{\overset{G}{\&Sigma;}}{\left(x\right(i,j)-\stackrel{\&OverBar;}{x})}^2}{K\&times;G-1}}$

这里，是x的均值，1≤i≤K，1≤j≤G，在Tetrolet域的高频系数矩阵中K＝12，G＝1；

步骤8Tetrolet变换中的拼板覆盖值的选取根据高频部分的选择而定，即

${\left({\mathrm{TC}}_F\right)}_m(p,q)=\left\{\begin{array}{cc}{\left({\mathrm{TC}}_A\right)}_m(p,q),& {\mathrm{TH}}_m(p,q)={\left({\mathrm{TH}}_A\right)}_m(p,q)\\ {\left({\mathrm{TC}}_B\right)}_m(p,q),& {\mathrm{TH}}_m(p,q)={\left({\mathrm{TH}}_B\right)}_m(p,q)\end{array}\right.$

其中，(TC_A)_m(p,q)和(TC_B)_m(p,q)分别表示图像A和B经Tetrolet变换后第m层p行q列的拼板覆盖值；

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

2.根据权利要求1所述的基于Tetrolet变换的多通道卫星云图融合方法，其特征在于：对图像A'和B'进行Tetrolet变换，基于Tetrolet变换的图像分解算法主要步骤包括：

步骤1把图像分成4×4子块；

步骤2考虑117种四格拼板分割方案，应用Haar小波变换获得每一种方案的高频和低频系数，选择Tetrolet系数小的方案，得到每个子块中的稀疏Tetrolet表示；

步骤3重新排列每个子块中的高频和低频部分的系数，把它们排列成2×2的子块；

步骤4存储Tetrolet系数的高频部分；

步骤5对低频部分重复进行步骤1-4的操作。