基于二进制小波变换的MODIS多波段数据融合

许 晨,康 雪,吕达仁

(1.成都市气象局，成都 610000; 2.中国科学院大气物理研究所，北京 100029)

针对多个传感器获得的同一地区卫星数据，遥感数据融合算法通过空间配准等方法将各源数据中所包含的信息互补有机地结合起来，改善待观测对象的精确性和定量可靠性，从而有效提高数据利用率[1]。本文应用小波变换针对Terra MODIS多波段数据进行融合，采用最近邻采样、双线性内插、三次卷积法将数据进行重采样，并评价融合效果。

1 资料来源与方法

1.1 资料来源

本文选取MODIS上午星Terra于5月29日在中国西部地区的数据(Latitude：22.1150—59.6326°N，Longitude： 90.0406—130.9920°E)进行分析，MODIS Level-1B数据包含250m、500m和1000m三种空间分辨率，数据格式为HDF格式。读取Terra MODIS数据的波段1、4、3，其波段宽度分别为0.62—0.67m、0.545—0.5657m、0.459—0.4797m，近似于可见光的红、绿、蓝波段(红色：656.37nm，绿：546.1nm，蓝：486.1nm)，因1、4、3波段组合比较接近真彩色，故常选用这三个波段来表示MODIS影像，使得视觉效果更加真实。

1.2 二进制小波变换

本文基于Harr小波将待融合数据分解为3层小波系数，按照不同融合策略分别对每一级小波系数上的高、低频信息进行处理，综合高、低频信息的特征，再通过小波逆变换重构融合数据。融合结果不仅能保留源数据光谱特性，还能增强数据的空间细节层次信息。

图1 空间V0

Haar函数φ(x)的定义为[2]：

(1)

φ(x-k)为φ(x)向右(左)平移k(正数或负数)个单位后的结果。将所有线性组合

∑akφ(x-k),k∈Z,ak∈R

(2)

所构成的空间用V0表示，如图1所示[2]。

用Haar函数分析信号的高、低频，由函数

{…,φ(2nx+1),φ(2nx),φ(2nx-1),φ(2nx-2),…}

(3)

生成的函数Vn描述的是时间区间长度为1/2n内取常数的信号。n越大，信号变化越剧烈；n越小，表示信号变化越缓慢。

具体实现过程与算法为：

(1)配准源数据，采用最近邻内插、双线性内插法、三次卷积内插法三种方法分别将源数据重采样为大小一致的数据。其中，最近邻内插法基于待求点的四个邻域像素灰度，选取距离最近的像素灰度值赋给待求点。双线性内插法基于待求点的四个邻域像素灰度，通过从两个方向线性内插，获得待求点的灰度值。三次卷积内插法基于待求点的十六个邻域坐标，通过三个方向线性内插，得到待求点灰度值。

(2)设计小波变换的滤波器系数数列，设计分解滤波器，逐层抽样分解数据。

(3)按照不同融合策略分别对每一级小波系数上的高、低频信息进行处理，并合并系数。本文选取的融合规则是：

图2 融合规则

加权平均法的基本原理是[5]：

CF(m,n)=α*CA(m,n)+(1-α)*CB(m,n)

(4)

其中，CA,CB,CF分别表示源数据A、B及融合后数据F的在最后一个分解层的子数据，m,n为最后一个分解层子数据的像素位置，α为权值，值在0-1。

(5)

最后使用重构滤波器逆向逐层二插值重构数据，获取融合后数据。

1.3 算法实现

图3 算法整体实现框图

图4 小波分解实现框图

图5 小波重构框图

2 结果与分析

2.1 读取数据

读取Terra MODIS数据的波段1、4、3。如图6所示。

观察图6可见光单通道数据，图中白色区域为云区，灰色区域为地表，黑色区域为水体。250m分辨率的数据(a)中的细节纹理比500m数据(b)、(c)、(d)更加丰富。对比(b)、(c)、(d)，发现随着波段增加，影像颜色越白亮。原因有二：一是不同通道对不同地物的反射率不同；二是这与仪器本身的光学性能和电子性能有关，比如放大器倍数不同等。

(a) 250m Band 1 (b) 500m Band 1 (c) 500m Band 3 (d) 500m Band 4图6 源数据

2.2 数据配准

在数据配准时，需要对500m数据进行重采样为250m数据，采用三种方法进行：最近邻内插法、双线性内插法以及三次卷积内插法。

其结果如图7-图9所示。

(a) Band 1 (b) Band 3 (c) Band 4图7 最近邻内插法对500m数据重采样

(a) Band 1 (b) Band 3 (c) Band 4图8 双线性内插法对500m数据重采样

(a) Band 1 (b) Band 3 (c) Band 4 图9 三次卷积内插法对500m数据重采样

2.3 小波融合

采用二进制小波融合算法对以下数据进行融合。

(1)对原500m 分辨率的Band 1、4、3进行融合，如图10 (a)所示；将分辨率为500m数据重采样为250m后，对Band 1、4、3进行融合，10(b)、(c)、(d)分别对应用最近邻内插、双线性内插以及三次卷积内插法采样后的融合结果。

(2)将空间分辨率为500m的数据重采样为250m后， Band 3、4与原250m分辨率的Band 1 融合，如图11，(a)、(b)、(c)分别对应最近邻内插、双线性内插以及三次卷积内插法采样后的融合结果。

图10 小波融合的数据结果Ⅰ

图11 小波融合的数据结果Ⅱ

从肉眼上看，图10和11的结果相似，故采用融合评价标准对其进行评价。

3 融合评价

融合评价通过分析和计算评估融合算法的性能、效果，可根据评价结果改进和优化算法，方法主要有主观和客观评价两种。

3.1 主观评价

主观评价基于人的视觉器官和经验知识，通过目视判别来分析评价数据融合的质量，该结果受人为因素影响较大。

分析融合结果中不同颜色代表的不同对象：亮色区域(白色和粉红色)代表云区；黑色区域是水体；绿色区域代表地表，有植被等信息。因本文选取的是中国西部地区的数据，故结合地理信息进行分析，图像中下部分偏左位置的小块黑色区域为青海湖，图像上方月牙形的黑色区域为贝加尔湖。

从视觉效果来看，三波段的融合结果增强了观测对象的空间细节表达能力，整幅图的纹理特征更加清晰，地表、水体等地物之间的分界线都能清晰辨认，所以融合结果图10和图11比图6的源数据细节信息更丰富，色彩细节信息越丰富对特征的提取和分类越有利，数据质量明显变好。但用肉眼观察图10和图11结果相似，判别不出采样方法的优劣。

故先针对三种采样方式的算法本身进行分析：最近邻内插值法计算过程最简单，但插值结果可能存在灰度不连续性，甚至在图像中表现为锯齿状；双线性内插法比最近邻插值法计算量大，克服了灰度不连续性的缺点，但算法本身有低通滤波的性质，导致高频分量受损，会使图像的轮廓有一定模糊；三次卷积内插法计算量最大、计算精度较高，能够克服以上两种算法的缺陷，但运算速度相对较慢。由于本文处理的数据像素很高，像素之间的梯度较小，所以从肉眼分辨不出三种采样方式的优劣。

3.2 客观评价

客观评价对数据的特性和统计参数进行分析[3-4]，一般综合两类统计参数，如表1，计算结果如表2。

表1 客观评价标准

表2 融合后数据的性能参数结果

从表2的统计参数结果中可看出，

(1)两组融合后数据——图10和11的方差、信息熵和清晰度均大于融合前的源数据，说明融合后数据包含的信息量增大，像素组合方式种类更多，纹理变化特征更加丰富。

(2)比较源数据——图6(a)、(b)，由表可知，(a)的三项指标均优于(b)。说明选择同一通道时，数据的分辨率越高，其方差、信息熵、清晰度越大，表达的信息量越大。

(3)分析融合结果Ⅰ——图10，(a)是原500m 分辨率的Band 1、4、3融合后结果，其各项指标最差，说明三波段融合后，融合结果的分辨率越高，指标越好；比较采用三种采样方式后的融合结果(b)、(c)、(d)，发现采用三次卷积内插后的结果(d)最优，其方差、信息熵、清晰度和相关系数最大，扭曲程度和偏差指数最小。双线性内插后的结果(c)次之，最近邻内插后的结果(a)最差。这和3.1中对采样算法分析的结论一致。

(4)分析融合结果Ⅱ，结论同(3)，三次卷积内插后的融合效果最优。

(5)比较融合结果Ⅰ和Ⅱ，Ⅰ中(b)、(c)、(d)是将重采样后的Band 1、4、3融合，Ⅱ是将采样后的Band3、4与源Band 1数据融合。由表可见，当选用的采样方式相同时，对应的结果Ⅱ优于结果Ⅰ，说明采样丢失了细节信息。

(6)观察对比所有融合结果，发现各参数之间的变化并不大，这是因为本文处理的数据像素很高，像素之间的梯度较小，所以差别不大。综合来看，融合结果Ⅱ中(c)的各项指标最优，其方差、信息熵、清晰度和相关系数最大，扭曲程度和偏差指数最小。证明采用三次卷积内插法重采样Band 3、4后与源250mBand 1融合具有最佳效果。

4 结论

本文提出了一种快速有效的多波段遥感数据融合方法，并通过融合评价比较了三次卷积内插法、双线性内插法、邻内插法的优劣，结果表示采用三次卷积内插法重采样后的融合效果最好。本文选用MODIS Band 1、4、3融合彩色图像，下阶段应进一步分析其他波段组合的物理意义，理解不同组合所突出的地物特征有何不同；分析融合后数据的物理特性，理解不同融合结果的差异在何处、具体表现为何特征、对地物分类有何意义等。