CN1770201A - 基于小波变换的可调节遥感影像融合方法 - Google Patents

Abstract

本发明涉及一种基于小波变换的可调节遥感影像融合方法，通过将高分辨率影像和多光谱影像进行小波多层分解，得到不同分解水平上各自的低频近似子图像和高频细节子图像，在每一分解水平上，采用一个滑动窗口，计算高分辨率影像和多光谱影像的子图像在滑动窗口内的方差之比，以该比值为依据，使用可调节参数对待融合子图像进行线性计算，依次得到各分解水平上对应子图像融合后的小波系数，最后进行小波逆变换得到融合影像。本发明在传统的小波分解与重构算法的基础上，通过引入可调节参数，很方便地在融合后影像的细节保留和光谱保持之间进行合理的取舍，可以得到单独使用某种传统的融合方法所无法实现的效果，从而可以满足不同的应用需求。

Description

基于小波变换的可调节遥感影像融合方法

技术领域

本发明属于图像分析领域，尤其是涉及多源遥感影像后处理阶段的一种基于小波变换的可调节遥感影像融合方法。

背景技术

遥感影像融合是一种高级影像处理技术，其目的是将单一传感器的多波段信息或不同类别传感器所提供的信息加以综合，消除多传感器信息之间可能存在的冗余和矛盾，加以互补，改善遥感信息提取的及时性和可靠性，提高数据的使用效率。将高空间分辨率的影像和多光谱影像进行融合的主要技术有：IHS变换法、PCA(主成分分析)变换法、HPF(高通滤波)方法以及小波变换法等。

然而，对上述融合算法进行归纳，其共同的缺点是：对一个特定的方法来说，只能得到一种特定的融合效果，即：用户无法根据影像的特点和实际应用的需要在高分辨率的保留和光谱信息保持之间进行动态的调整。但是，不同的应用对这两个面有不同的要求。对于特征提取、纹理提取、地图制图等应用来说，影像的细节信息较为重要；对于遥感影像分类和城市制图来说，保持光谱信息的准确、完整，则是至关重要的。所以，提出能够为不同的应用提供不同的融合效果的方法，就显得尤为重要。

发明内容

本发明针对上述现有影像融合技术存在的不足，提供一种基于小波变换的可调节遥感影像融合方法，其可以在融合影像的光谱保持能力和高空间分辨率保留能力之间进行动态调整。

为实现这样的目的，本发明在将高分辨率影像和多光谱影像进行小波多层分解之后，得到不同分解水平上各自的低频近似子图像和高频细节子图像。为了综合考虑两幅影像的细节特征，在小波分解的每一水平上，采用一个滑动窗口，计算高分辨率影像的子图像和多光谱影像的子图像在滑动窗口内的方差之比，以该比值为依据，使用可调节参数对待融合的子图像进行线性计算，依次得到各个分解水平上多光谱影像和高分辨率影像对应子图像融合后的小波系数，最后进行小波逆变换得到融合影像。

本发明提供的技术方案具体如下：

一种基于小波变换的可调节遥感影像融合方法，其包含有如下步骤：

A、将高分辨率影像A进行线性拉伸并与多光谱影像B进行直方图匹配，然后进行小波分解，得到低频近似子图像S_A(2^j；x，y)和高频细节子图像W_kA(2^j；x，y)，其中k＝1，2，3表示3个方向；j＝1，2，...，J表示不同的分辨率；

将影像B进行小波分解，得到低频近似子图像S_B(2^j；x，y)和高频细节子图像W_kB(2^j；x，y)；

B、确定权系数q_k值：

在第j分解水平上，分别对影像A和B的子图像进行统计，构造R_k(2^j；x，y)，用以反映两幅影像对应子图像的细节信息的相对大小，公式(1)中：D_kA(2^j；x，y)与D_kB(2^j；x，y)分别是W_kA(2^j；x，y)与W_kB(2^j；x，y)中以(x，y)为中心像元的n×n，n＝3，5，...，K窗口内的方差；

$R_k(2^j;x,y)=\frac{D_{\mathrm{kA}}(2^j;x,y)}{D_{\mathrm{kB}}(2^j;x,y)}$ 公式(1)

将R_k(2^j；x，y)进行归一化，公式(2)中：R_{k_normal}(2^j；x，y)是归一化的R_k(2^j；x，y)，R_kmin(2^j)是R_k(2^j；x，y)的最小值，R_kmax(2^j)是R_k(2^j；x，y)的最大值；

$R_{k\_\mathrm{normal}}(2^j;x,y)=\frac{R_k(2^j;x,y)-R_{k\min }\left(2^j\right)}{R_{k\max }\left(2^j\right)-R_{k\min }\left(2^j\right)}$ 公式(2)

计算影像融合权系数q_k的值，公式(3)中：a和b是可调节参数(0≤a≤b≤1)；

$q_k=\left\{\begin{array}{cc}0& R_{k-\mathrm{normal}}\&le;a\\ \frac{R_{k\_\mathrm{normal}}-a}{b-a}& a<R_{k\_\mathrm{normal}}<b\\ 1& R_{k\_\mathrm{normal}}\&GreaterEqual;b\end{array}\right.$ 公式(3)

C、确定融合值：

根据权系数q_k的值，将两幅影像对应子图像W_kA(2^j；x，y)和W_kB(2^j；x，y)在(x，y)处的小波系数进行融合，得到(x，y)处的融合值，其中：F_k(2^j；x，y)是融合后的子图像的滑动窗口中心点(x，y)的值；

  F_k(2^j；x，y)＝q_k×W_kA(2^j；x，y)+(1-q_k)×W_kB(2^j；x，y)

D、重复步骤B和C，依次对每一组低频近似子图像S_A(2^j；x，y)和S_B(2^j；x，y)，以及所有的分解水平上进行融合操作，得到融合影像在各分解水平上的低频近似子图像和高频细节子图像；

对得到的小波系数进行小波逆变换，得到高光谱影像的一个波段与高分辨率影像的融合影像；再依次将多光谱影像B的其他波段也与高分辨率影像A按照前述步骤进行融合，从而完成所有波段的融合。

本发明是在传统的小波分解与重构算法的基础上，通过引入两个可调节参数，来灵活控制两幅不同分辨率影像的融合效果，很方便地在融合后影像的细节保留和光谱保持之间进行合理的取舍。由于引入了可调节参数，可以获得单独使用某种传统的融合方法所无法实现的效果。在合理的参数组合下，融合图像在细节保留和光谱保持方面均要超过某些传统的融合方法，从而可以满足不同的应用需要，增强遥感影像的信息清晰度、改善解译的精度以及提高信息的可靠性和使用效率。

附图说明

图1是本发明所述方法与习用方法的融合效果比较图；

图2是本发明所述方法的流程图。

具体实现方式

为了更加清楚本发明的技术方案，兹配合附图和具体实施例详细说明如下：

本发明选取SPOT全色波段影像作为高分辨率影像A(256×256)如图1(a)，选取Landsat TM影像的4、3、2波段作为多光谱影像B(256×256)，如图1(b)；A、B严格配准后，具体步骤如下，如图2所示：

1、将高分辨率影像A进行线性拉伸并与B进行直方图匹配，然后进行小波分解，选用的小波基为db4，分解层数为2。该小波分解既可以用DB4，也可用其他的小波基，如：DB1，DB2，Harr，...，同时分解层数也可以是1，2，3，4，...选用的小波基和小波分解的层数可以是任意的。本实施例仅以小波基为db4，分解层数为2进行说明。得到第j分解水平的低频近似子图像S_A(2^j；x，y)和3个高频细节子图像W_1A(2^j；x，y)，W_2A(2^j；x，y)，W_3A(2^j；x，y)；

2、将多光谱影像B的第2波段进行小波分解，同样选用db4作为小波基，分解层数为2。得到第j分解水平的低频近似子图像S_B(2^j；x，y)和3个高频细节子图像W_1B(2^j；x，y)，W_2B(2^j；x，y)，W_3B(2^j；x，y)；

3、确定一个3×3的窗口，在第j分解水平上，分别对影像A和B的子图像进行统计，求得小窗口内的方差D_1A(2^j；x，y)和D_1B(2^j；x，y)，按照公式(1)在每一分解水平上构造R_k(2^j；x，y)。在第j分解水平上：

$R_1(2^j;x,y)=\frac{D_{1A}(2^j;x,y)}{D_{1B}(2^j;x,y)}$

$R_2(2^j;x,y)=\frac{D_{2A}(2^j;x,y)}{D_{2B}(2^j;x,y)}$ 公式(1)

$R_3(2^j;x,y)=\frac{D_{3A}(2^j;x,y)}{D_{3B}(2^j;x,y)}$

4、在第j分解水平上，将R₁(2^j；x，y)、R₂(2^j；x，y)和R₃(2^j；x，y)进行归一化处理，得到：

$R_{1\_\mathrm{normal}}(2^j;x,y)=\frac{R_1(2^j;x,y)-R_{k\min }\left(2^j\right)}{R_{k\max }\left(2^j\right)-R_{k\min }\left(2^j\right)}$

$R_{2\_\mathrm{normal}}(2^j;x,y)=\frac{R_2(2^j;x,y)-R_{k\min }\left(2^j\right)}{R_{k\max }\left(2^j\right)-R_{k\min }\left(2^j\right)}$ 公式(2)

$R_{3\_\mathrm{normal}}(2^j;x,y)=\frac{R_3(2^j;x,y)-R_{k\min }\left(2^j\right)}{R_{k\max }\left(2^j\right)-R_{k\min }\left(2^j\right)}$

根据公式(3)计算在第j分解水平上子图像间融合的权系数q_k，对于W_1A(2^j；x，y)和W_1B(2^j；x，y)来说：

$q_1=\left\{\begin{array}{cc}0& R_{1-\mathrm{normal}}\&le;a\\ \frac{R_{1\_\mathrm{normal}}-a}{b-a}& a<R_{1\_\mathrm{normal}}<b\\ 1& R_{1\_\mathrm{normal}}\&GreaterEqual;b\end{array}\right.$ 公式(3)

通过调节参数a和b，就可以达到在细节保留和光谱保持两方面进行调整的目的：若R_{k_normal}大于或等于b，则R_{k_normal}等于1，中心点的值将完全由高分辨率影像所决定；若R_{k_normal}小于或等于a，则R_{k_normal}等于0，融合影像中心点的值将取决于多光谱影像；否则，就使用高分辨率影像和多光谱影像的加权值作为融合结果。R_{k_normal}将决定最终的融合效果：R_{k_normal}越大(≤1)，融合影像的细节保留能力越强，而光谱保持能力下降；反之，R_{k_normal}越小(≥0)，融合影像的光谱保持能力越强，而细节信息的保留则相对减弱。

根据权系数q_k的值，将两幅影像对应子图像W_kA(2^j；x，y)和W_kB(2^j；x，y)进行融合，得到(x，y)处经过融合后的小波系数的融合值F_k(2^j；x，y)：

   F₁(2^j；x，y)＝q₁×W_1A(2^j；x，y)+(1-q₁)×W_1B(2^j；x，y)

   F₂(2^j；x，y)＝q₂×W_2A(2^j；x，y)+(1-q₂)×W_2B(2^j；x，y)

   F₃(2^j；x，y)＝q₃×W_3A(2^j；x，y)+(1-q₃)×W_3B(2^j；x，y)

7、将窗口进行滑动，依次计算第j分解水平上的子图像内每个小波系数融合后的融合值；

8、对每一组低频近似子图像S_A(2^j；x，y)和S_B(2^j；x，y)，也按照步骤3～7进行融合；

9、在所有的分解水平上按照步骤3～8进行操作，得到融合影像在各分解水平上的低频近似子图像和高频细节子图像；

10、对得到的小波系数(子图像)进行小波逆变换，得到高光谱影像的一个波段与高分辨率影像的融合影像；

11、将多光谱影像B的第3和第4波段也与高分辨率影像A按照步骤2～10进行融合，完成所有波段的融合。

如图1所示，为本发明采用的方法与其他方法的融合效果比较。所有影像的大小均为256×256。

其中，图1(a)的SPOT全色影像是高分辨率影像，其分辨率为10米；图1(b)的Landsat TM图像为多光谱遥感影像，本实例取其波段4、3、2，空间分辨率为30米；图1(c)为IHS方法的融合结果；图1(d)为PCA方法的融合结果；图1(e)为小波系数替换法的融合结果；图1(f)为基于特征选择的小波系数替换法的融合结果；图1(g)为本发明a＝0.001b＝0.1时的融合结果，图1(h)为本发明a＝0.001 b＝0.3时的融合结果，图1(i)是本发明a＝0.001 b＝0.7时的融合结果，图1(j)是本发明a＝0.01 b＝0.2时的融合结果，图1(k)为本发明a＝0.1 b＝0.2时的融合结果，图1(l)是本发明a＝0.7 b＝0.9时的融合结果。从中可以看出，不同的参数组合可以得到不同的融合效果，某些效果在光谱保持和分辨率保留上均优于其他方法。

性能分析：表1列出了选取不同的融合参数时，融合结果的第2波段的均值、信息熵、平均梯度、光谱扭曲以及相关系数的定量值，并与其他方法进行了对比。其融合后是有3个波段，但是如果在此将这3个波段一一列出则比较凌乱。因此为了表述清楚，在此只列出第2波段的统计参数值。

参数a和b的改变，对融合结果有很大的影响。当a和b接近于0时，融合结果的清晰度很高，影像的特征较接近于SPOT影像，而光谱保持能力相对有所下降；a和b增大，融合影像的光谱保持能力增强，图像特征越来越接近于TM影像，而细节和纹理的保留能力则不断下降。这说明，通过调节参数，可以在细节保留和光谱保持之间达到不同程度的平衡。而且，某些参数组合的结果要优于传统的融合方法。比如，从均值来看，新方法的结果比其他方法都更接近于原来的TM影像的均值，而光谱扭曲都低于其他方法，相关系数则均在0.9以上，说明新方法具有较强的光谱保持能力；从信息熵来看，除a＝0.7 b＝0.9以外，其他参数组合的信息熵都接近或高于其他的融合方法。平均梯度则可以从2.205变换到4.867，基本可以达到其他几种融合方法的融合效果。

表1 融合结果的性能评价指标比较

图像	均值	信息熵	平均梯度	光谱扭曲	相关系数
						SPOT	102.832	5.8888	5.026
TM	93.771	5.3640	2.185
						IHS	104.610	5.7012	5.107	13.428	0.61288
PCA	110.530	5.7558	4.767	16.858	0.89350
						WR	121.842	5.9328	5.825	28.089	0.92808
WS	125.013	5.9758	3.843	31.245	0.96708
						新方法(a＝0.001 b＝0.1)(a＝0.001 b＝0.3)(a＝0.001 b＝0.7)(a＝0.01  b＝0.2)(a＝0.1   b＝0.2)(a＝0.7   b＝0.9)	102.69299.82496.510101.694100.81093.785	5.9015.7475.5375.8565.8495.365	4.8673.8342.8514.3664.1762.205	9.5206.4702.9528.4137.5980.039	0.9040.9300.9690.9170.9170.999

Claims (4)

1、一种基于小波变换的可调节遥感影像融合方法，其特征在于，该方法包含有如下步骤：

A、将高分辨率影像A进行线性拉伸并与多光谱影像B进行直方图匹配，然后进行小波分解，得到低频近似子图像S_A(2^j；x，y)和高频细节子图像W_kA(2^j；x，y)，其中k＝1，2，3表示3个方向；j＝1，2，…，J表示不同的分辨率；

将影像B进行小波分解，得到低频近似子图像S_B(2^j；x，y)和高频细节子图像W_kB(2^j；x，y)；

B、确定权系数q_k值：

在第j分解水平上，分别对影像A和B的子图像进行统计，构造R_k(2^j；x，y)，用以反映两幅影像对应子图像的细节信息的相对大小，公式(1)中：D_kA(2^j；x，y)与D_kB(2^j；x，y)分别是W_kA(2^j；x，y)与W_kB(2^j；x，y)中以(x，y)为中心像元的n×n，n＝3，5，…，K窗口内的方差；

$R_k(2^j;x,y)=\frac{D_{\mathrm{kA}}(2^j;x,y)}{D_{\mathrm{kB}}(2^j;x,y)}$ 公式(1)

将R_k(2^j；x，y)进行归一化，公式(2)中：R_{k_normal}(2^j；x，y)是归一化的R_k(2^j；x，y)，R_kmin(2^j)是R_k(2^j；x，y)的最小值，R_kmax(2^j)是R_k(2^j；x，y)的最大值；

$R_{k\_\mathrm{normal}}(2^j;x,y)=\frac{R_k(2^j;x,y)-R_{k\min }\left(2^j\right)}{R_{k\max }\left(2^j\right)-R_{k\min }\left(2^j\right)}$ 公式(2)

计算影像融合权系数q_k的值，公式(3)中：a和b是可调节参数(0≤a≤b≤1)；

$q_k=\left\{\begin{array}{cc}0& R_{k-\mathrm{normal}}\&le;a\\ \frac{R_{k\_\mathrm{normal}}-a}{b-a}& a<R_{k\_\mathrm{normal}<b}\\ 1& R_{k\_\mathrm{normal}}\&GreaterEqual;b\end{array}\right.$ 公式(3)

C、确定融合值：

根据权系数q_k的值，将两幅影像对应子图像W_kA(2^j；x，y)和W_kB(2^j；x，y)在(x，y)处的小波系数进行融合，得到(x，y)处的融合值，其中：F_k(2^j；x，y)是融合后的子图像的滑动窗口中心点(x，y)的值；

       F_k(2^j；x，y)＝q_k×W_kA(2^j；x，y)+(1-q_k)×W_kB(2^j；x，y)

D、重复步骤B和C，依次对每一组低频近似子图像S_A(2^j；x，y)和S_B(2^j；x，y)，以及所有的分解水平上进行融合操作，得到融合影像在各分解水平上的低频近似子图像和高频细节子图像；

对得到的小波系数进行小波逆变换，得到高光谱影像的一个波段与高分辨率影像的融合影像；再依次将多光谱影像B的其他波段也与高分辨率影像A按照前述步骤进行融合，从而完成所有波段的融合。

2、如权利要求1所述的基于小波变换的可调节遥感影像融合方法，其特征在于：所述的步骤B更具体为：

确定一个3×3的窗口，在第j分解水平上，分别对影像A和B的子图像进行统计，求得小窗口内的方差D_1A(2^j；x，y)和D_1B(2^j；x，y)，按照公式(1)在每一分解水平上构造R_k(2^j；x，y)。在第j分解水平上：

$R_1(2^j;x,y)=\frac{D_{1A}(2^j;x,y)}{D_{1B}(2^j;x,y)}$

$R_2(2^j;x,y)=\frac{D_{2A}(2^j;x,y)}{D_{2B}(2^j;x,y)}$

$R_3(2^j;x,y)=\frac{D_{3A}(2^j;x,y)}{D_{3B}(2^j;x,y)}$ 公式(1)

在第j分解水平上，将R₁(2^j；x，y)、R₂(2^j；x，y)和R₃(2^j；x，y)进行归一化处理，得到：

$R_{1\_\mathrm{normal}}(2^j;x,y)=\frac{R_1(2^j;x,y)-R_{k\min }\left(2^j\right)}{R_{k\max }\left(2^j\right)-R_{k\min }\left(2^j\right)}$

$R_{2\_\mathrm{normal}}(2^j;x,y)=\frac{R_2(2^j;x,y)-R_{k\min }\left(2^j\right)}{R_{k\max }\left(2^j\right)-R_{k\min }\left(2^j\right)}$

$R_{3\_\mathrm{normal}}(2^j;x,y)=\frac{R_3(2^j;x,y)-R_{k\min }\left(2^j\right)}{R_{k\max }\left(2^j\right)-R_{k\min }\left(2^j\right)}$ 公式(2)

根据公式(3)计算在第j分解水平上子图像间融合的权系数q_k：

$q_1=\left\{\begin{array}{cc}0& R_{1-\mathrm{normal}}\&le;a\\ \frac{R_{1\_\mathrm{normal}}-a}{b-a}& a<R_{1\_\mathrm{normal}<b}\\ 1& R_{1\_\mathrm{normal}}\&GreaterEqual;b\end{array}\right.$ 公式(3)

通过调节参数a和b，可以达到在细节保留和光谱保持两方面进行调整的目的，即：若R_{k_normal}大于或等于b，则R_{k_normal}等于1，中心点的值将完全由高分辨率影像所决定；若R_{k_normal}小于或等于a，则R_{k_normal}等于0，融合影像中心点的值将取决于多光谱影像；否则，就使用高分辨率影像和多光谱影像的加权值作为融合结果。

3、如权利要求1所述的基于小波变换的可调节遥感影像融合方法，其特征在于：所述的步骤C更具体为：

根据权系数q_k的值，将两幅影像对应子图像W_kA(2^j；x，y)和W_kB(2^j；x，y)进行融合，得到(x，y)处经过融合后的小波系数的融合值F_k(2^j；x，y)：

          F₁(2^j；x，y)＝q₁×W_1A(2^j；x，y)+(1-q₁)×W_1B(2^j；x，y)

          F₂(2^j；x，y)＝q₂×W_2A(2^j；x，y)+(1-q₂)×W_2B(2^j；x，y)

          F₃(2^j；x，y)＝q₃×W_3A(2^j；x，y)+(1-q₃)×W_3B(2^j；x，y)

4、如权利要求1所述的基于小波变换的可调节遥感影像融合方法，其特征在于：所述的步骤A中，将高分辨率影像A进行线性拉伸并与B进行直方图匹配，进行小波分解时，选用的小波基为db4，分解层数为2；将多光谱影像B进行小波分解，同样选用db4作为小波基，分解层数为2。

Images