CN1313972C

CN1313972C - 基于滤波器组的图像融合方法

Info

Publication number: CN1313972C
Application number: CNB031417884A
Authority: CN
Inventors: 敬忠良; 王宏; 李建勋
Original assignee: Shanghai Jiaotong University
Current assignee: Shanghai Jiaotong University
Priority date: 2003-07-24
Filing date: 2003-07-24
Publication date: 2007-05-02
Anticipated expiration: 2023-07-24
Also published as: CN1484040A

Abstract

本发明涉及一种基于滤波器组的图像融合方法，在对原图像进行配准的基础上，对低空间分辨率多光谱图像进行插值重采样，以实现与高空间分辨率全色遥感图像具有相同的尺寸大小，在设计出满足融合要求的滤波器组的基础上，对高空间分辨率全色遥感图像和低空间分辨率多光谱图像分别进行滤波器组分解，分别得到多个频带的子图像，采用低空间分辨率多光谱图像的低频部分替换相应高空间分辨率全色遥感图像的低频部分的方式，融合不同分辨率的遥感图像，最后，进行滤波器组重构以获得融合后的图像。本发明在保留融合图像的光谱信息的同时，提高图像的空间信息，能够使空间细节信息和光谱信息二者之间达到更好的均衡，使得融合后图像的质量得以较大的提高。

Description

基于滤波器组的图像融合方法

技术领域：

本发明涉及一种基于滤波器组的图像融合方法，用以融合高空间分辨率全色遥感图像和低空间分辨率的多光谱图像。在各类民用或军用的遥感信息处理系统、数字城市空间信息系统等领域中均有广泛的应用。

背景技术：

随着成像传感器技术的飞速发展，使得遥感图像处理从单平台，单一传感器的处理方式向多平台，多传感器的多图像处理方式发展。由于不同成像传感器的使用，不同的遥感数据具有不同的空间分辨率和光谱分辨率。对于多光谱传感器来说，其能够获得光谱分辨率高而空间分辨率低的多光谱图像；对于全色波段的传感器，可得到空间分辨率高的遥感图像。利用图像融合技术能将它们各自的优势结合起来，得到空间分辨率高的多光谱图像，弥补了单一图像上信息的不足，不仅扩大了应用范围，而且大大提高了遥感影像分析的精度。图像融合的宗旨就是利用各种成像传感器不同的成像方式，为不同的影像提供互补信息，增加图像的信息量。

对遥感图像而言，多光谱图像包含着丰富的光谱信息而全色图像其空间分辨率较高，这样其互补信息分布有着明显的不同，多光谱图像其光谱信息集中于图像的低频部分；而全色图像其高分辨率空间信息分布于高频部分。融合的目的就是在融合图像中低频部分具有多光谱图像的光谱信息而高频部分具有全色图像的高分辨率空间信息。因此，在遥感图像融合中，一般采用将全色图像的低频部分替换为多光谱图像的低频部分的方式或者将全色图像的高频部分叠加到多光谱图像中的方式进行融合处理。

至今为止，人们已发展了多种方法用以融合不同分辨率的遥感图像，其中包括：IHS(强度、色度、饱和度)变换、主成分分析、高通滤波法等经典方法。随着小波理论的发展，基于离散小波变换(DWT)以及小波框架(DWF)等融合方法已广泛应用于遥感领域。近年来，随着研究的深入，人们已将多速率滤波器组的方法引入到图像融合领域。Blanc提出了基于两通道分数滤波器组的融合方法；Argenti设计了用于图像融合的余弦调制滤波器组。小波变换和滤波器组这两种融合方法的融合思想是相同的，即将高空间分辨率全色遥感图像的高频部分替换低空间分辨率多光谱图像的高频部分或将高空间分辨率全色遥感图像的高频部分叠加到低空间分辨率的全色图像。对于二进制小波变换，它只能实现两通道信号分解，不能作到多个通道的信号分解。

发明内容：

本发明的目的在于针对现有技术存在的不足，提供一种基于滤波器组的图像融合方法，通过滤波器组的创新设计以获得较满意的融合信号，设计所得的滤波器组可用于不同分辨率遥感图像融合领域。

为实现这样的目的，本发明在对原图像进行配准的基础上，对低空间分辨率多光谱图像进行插值重采样，以实现与高空间分辨率全色遥感图像具有相同的尺寸大小。在设计出满足融合要求的滤波器组的基础上，对高空间分辨率全色遥感图像和低空间分辨率多光谱图像分别进行滤波器组分解，分别得到多个频带的子图像。采用低空间分辨率多光谱图像的低频部分替换相应高空间分辨率全色遥感图像的低频部分的方式，融合不同分辨率的遥感图像；最后，进行滤波器组重构以获得融合后的图像。

本发明的方法包括如下具体步骤：

1、在对同一地貌的低空间分辨率多光谱图像和高空间分辨率全色遥感图像进行配准的基础上，通过对低空间分辨率多光谱图像进行插值重采样，以实现二者具有相同的尺寸大小。

对于图像的插值重采样，可采用常规的最近邻法、双线性插值法以及三次卷积插值等方法。

2、考虑高空间分辨率全色遥感图像和低空间分辨率多光谱图像之间的光谱差别以及不同气候和照明条件的影响，对高空间分辨率全色遥感图像和低空间分辨率多光谱图像进行直方图匹配，使它们的灰度均值和方差相一致。

3、设计满足融合要求的滤波器组。融合的目的就是产生一个融合信号，使该信号的低频部分具有低分辨率信号的低频信息；而高频部分具有高分辨率信号的高频信息。在融合过程中，使频率信息的失真最小。本发明采用对原型滤波器进行余弦调制的方法来产生分析和综合滤波器组，当原型滤波器满足：

\{\begin{matrix} G (ω) = H (ω) & 0 \leq ω < π / M_{1} \\ 0 & π / M_{1} \leq ω \leq π \end{matrix} - - - (1)

将使得融合信号的频率信号失真最小。其中，G(ω)为用于分解高分辨率信号的原型滤波器；H(ω)为用于分解低分辨率信号的原型滤波器；M₁为通道数。

4、在得到满足融合要求的滤波器组的基础上，将高空间分辨率全色遥感图像和低空间分辨率多光谱图像分别进行滤波器组分解，得到多个频带的子图像。在分解时，可分为保留降采样过程和忽略降采样过程两种情况。其中，原型滤波器满足(1)式的要求。

5、将得到的低空间分辨率多光谱图像低频部分替换高空间分辨率全色遥感图像的低频部分，高空间分辨率全色遥感图像的高频部分保持不变。

6、对替代后得到的子图像进行滤波器组重构以获得融合后的图像。

本发明提出了一种基于滤波器组的图像融合方法，具体的有益效果为：在对输入原图像进行基于滤波器组分解的基础上，采用将低空间分辨率多光谱图像的低频部分替换高空间分辨率全色遥感图像的低频部分的方式进行不同分辨率的图像融合，在保留融合图像的光谱信息的同时，提高了图像的空间信息，能够使空间细节信息和光谱信息二者之间达到更好地的均衡，从而使得融合后图像的质量得以较大的提高，对于应用系统的后续处理具有重要意义和实用价值。

附图说明：

图1为本发明的基于滤波器组图像融合方法流程示意图。

如图1所示，本发明在对原图像配准的基础上，对低空间分辨率多光谱图像MS进行插值重采样，然后对高空间分辨率全色遥感图像PS和低空间分辨率多光谱图像MS进行滤波器组分解，分别得到多个频带的子图像P₀₀～P₂₂和Q₀₀～Q₂₂，通过用低空间分辨率多光谱图像的低频部分Q₀₀替换高空间分辨率全色遥感图像的低频部分P₀₀的方式，提高融合后图像的光谱信息和空间信息。最后，进行滤波器组重构以获得融合后的图像MS‘。

图2为用于融合的滤波器组设计示意图。

图3为本发明遥感影像原图以及与小波变换方法的融合结果比较。

其中，(a)为高空间分辨率全色遥感图像；(b)低空间分辨率多光谱图像；(c)～(f)3～6通道滤波器组融合图像；(g)小波变换融合图像(两层分解)；(h)小波变换融合图像(三层分解)；(i)～(l)3～6通道滤波器组融合图像(忽略降采样)；(m)小波框架融合图像(两层分解)；(n)小波框架融合图像(三层分解)；

具体实施方式：

为了更好地理解本发明的技术方案，以下结合附图对本发明的实施方式作进一步描述。本发明具体实施如下：

1、对同一地貌的低空间分辨率多光谱图像和高空间分辨率全色遥感图像进行配准，在此基础上，再对低空间分辨率多光谱图像进行插值重采样，以实现二者具有相同的尺寸大小。

采用三次卷积插值方法进行图像重采样的实施例所得的低空间分辨率多光谱图像和高空间分辨率全色遥感图像如图-3(a)，(b)所示。

2、对高空间分辨率全色遥感图像和低空间分辨率多光谱图像进行直方图匹配，使它们的灰度均值和方差相一致。

3、设计满足融合要求的滤波器组。在对滤波器组进行设计时，以一维情况为例，对于二维图像可采用可分离方式直接将一维滤波器组推广到二维情况。

具体如下：设x₁(n)和x₂(n)分别为具有不同分辨率的信号，对低分辨率的信号

x₁(n)进行插值重采样使得信号x₁′(n)与x₂(n)具有相同的长度。如图2所示，融合的目的就是产生一个融合信号使该信号的低频部分具有信号x₁′(n)的低频信息；而高频部分具有信号x₂(n)的高频信息。从图-2中，可得输出融合信号x′(n)表示为：

X^{'} (z) = \frac{1}{M_{1}} Σ_{l = 0}^{M_{1} - 1} X_{1} (z W^{l}) Σ_{k = 0}^{M_{2} - 1} H_{k} (z W^{l}) F_{k} (z) + \frac{1}{M_{1}} Σ_{l = 0}^{M_{1} - 1} X_{2} (z W^{l}) Σ_{k = M_{2}}^{M_{1}} G_{k} (z W^{l}) V_{k} (z) - - - (2)

在(2)式中，对于l＝0所得的信号为所需要的输出融合信号；而对于l≠0的那些信号为混叠信号。对于l＝0，可得，

X^{'} (z) = \frac{1}{M_{1}} Σ_{k = 0}^{M_{2} - 1} H_{k} (z) F_{k} (z) \cdot X_{1} (z) + \frac{1}{M_{1}} Σ_{k = M_{2}}^{M_{1} - 1} G_{k} (z) V_{k} (z) \cdot X_{2} (z) - - - (3)

从(3)式，可看到融合信号的低频部分与信号x₁(n)有关而其高频部分与信号x₂(n)有关。那么，若有：

\frac{1}{M_{1}} Σ_{k = 0}^{M_{2} - 1} H_{k} (z) F_{k} (z) = c \cdot e^{- jαw}, 0 \leq w \leq \frac{M_{2}}{M_{1}} π - ϵ - - - (4)

\frac{1}{M_{1}} Σ_{k = M_{2}}^{M_{1} - 1} G_{k} (z) V_{k} (z) = c \cdot e^{- jαw}, \frac{M_{2}}{M_{1}} π + ϵ < w \leq π

成立就可得到需要的融合结果。其中c，α为常数，对于ε是由于过渡带可能是理想的，频率信息的替换也不可能是完全的，为此在ω＝M₂·π/M₁处存在频率信息失真。

在此，采用余弦调制滤波器组进行信号的分解和重构，分析滤波器和重构滤波器分别为：

\{\begin{matrix} h_{k} (n) = 2 h (n) \cos ((2 k + 1) \frac{π}{2 M_{1}} (n - \frac{N_{H} - 1}{2}) + {(- 1)}^{k} \frac{π}{4}) \\ f_{k} (n) = 2 h (n) \cos ((2 k + 1) \frac{π}{2 M_{1}} (n - \frac{N_{H} - 1}{2}) - {(- 1)}^{k} \frac{π}{4}) \end{matrix} \{\begin{matrix} 0 \leq n \leq N - 1 \\ 0 \leq N \leq M_{1} - 1 \end{matrix} - - - (5)

\{\begin{matrix} g_{k} (n) = 2 g (n) \cos ((2 k + 1) \frac{π}{2 M_{1}} (n - \frac{N_{G} - 1}{2}) + {(- 1)}^{k} \frac{π}{4}) \\ v_{k} (n) = 2 g (n) \cos ((2 k + 1) \frac{π}{2 M_{1}} (n - \frac{N_{G} - 1}{2}) - {(- 1)}^{k} \frac{π}{4}) \end{matrix} \{\begin{matrix} 0 \leq n \leq N - 1 \\ 0 \leq N \leq M_{1} - 1 \end{matrix} - - - (6)

其中，h(n)和g(n)为具有线性相位，功率互补特性的原型滤波器。相应的频率响应为：

\{\begin{matrix} H_{k} (z) = a_{k} c_{k} H (z W_{2 M}^{k + 0.5}) + a_{k}^{*} c_{k}^{*} H (z W_{2 M}^{- (k + 0.5)}) \\ F_{k} (z) = a_{k}^{*} c_{k} H (2 W_{2 M}^{k + 0.5}) + a_{k} c_{k}^{*} H (z W_{2 M}^{- (k + 0.5)}) \end{matrix} 0 \leq k \leq M_{1} - 1 - - - (7)

\{\begin{matrix} G_{k} (z) = b_{k} d_{k} G (z W_{2 M}^{k + 0.5}) + b_{k}^{*} d_{k}^{*} G (z W_{2 M}^{- (k + 0.5)}) \\ V_{k} (z) = b_{k}^{*} d_{k} G (2 W_{2 M}^{k + 0.5}) + b_{k} d_{k}^{*} G (z W_{2 M}^{- (k + 0.5)}) \end{matrix} 0 \leq k \leq M_{1} - 1 - - - (8)

其中，上标“*”为复共轭；W_2M＝e^-jπ/M；

c_{k} = W^{(k + 0.5) N_{H} / 2}; d_{k} = W^{(k + 0.5) N_{G} / 2};

a_{k} = b_{k} = e^{j θ_{k}},

θ_k＝(-1)^kπ/4。

将(7)和(8)式代入(3)中，则有：

X^{'} (ω) = X_{1} (ω) \cdot \frac{1}{M_{1}} Σ_{k = 0}^{M_{2} - 1} [G^{2} (ω - \frac{π}{M_{1}} (k + \frac{1}{2})) + G^{2} (ω + \frac{π}{M_{1}} (k + \frac{1}{2}))] +

X_{2} (ω) \cdot \frac{1}{M_{1}} Σ_{k = M_{2}}^{M_{1} - 1} [H^{2} (ω - \frac{π}{M_{1}} (k + \frac{1}{2})) + H^{2} (ω + \frac{π}{M_{1}} (k + \frac{1}{2}))]

(9)

当

\{\begin{matrix} G (ω) = H (ω) & 0 \leq ω < π / M_{1} \\ 0 & π / M_{1} \leq ω \leq π \end{matrix} - - - (10)

那么，由于原型滤波器h(n)的功率互补性质，可得(4)式，其中c＝1/M₁及α＝N_H-1＝N_G-1。

考虑l≠0情况，混叠失真主要位于M₂·π/M₁处，在其他部分由于余弦调制滤波器组本身的设计已将其混叠相互抵消，而在M₂·π/M₁处混叠的产生是由不同信号造成的，因此它们不能相互抵消。

如果忽略降采样过程，那么混叠失真将被避免。输出融合信号x′(n)能够被简化为：

X^{'} (z) = Σ_{k = 0}^{M_{2} - 1} H_{k} (z) F_{k} (z) \cdot X_{1}^{'} (z) + Σ_{k = M_{2}}^{M_{1} - 1} G_{k} (z) V_{k} (z) \cdot X_{2} (z) - - - (11)

所得的原型滤波器仍满足(10)式。

将设计所得的滤波器组应用到不同分辨率遥感图像融合中，如图1所示。

4、将高空间分辨率全色遥感图像PS和低空间分辨率多光谱图像MS通过分析滤波器组，分别得到M*M频带的子图像P_ij(z₁，z₂)和Q_ij(z₁，z₂)：

P_{ij} (z_{1}, z_{2}) = \frac{1}{M^{2}} Σ_{l 1 = 0}^{M - 1} Σ_{l 2 = 0}^{M - 1} H_{i} ({z_{1}}^{1 / M} W^{l 1}) H_{j} ({z_{2}}^{1 / M} W^{l 2}) PS ({z_{1}}^{1 / M} W^{l 1}, {z_{2}}^{1 / M} W^{l 2}) - - - (12)

Q_{ij} (z_{1}, z_{2}) = \frac{1}{M^{2}} Σ_{l 1 = 0}^{M - 1} Σ_{l 2 = 0}^{M - 1} G_{i} ({z_{1}}^{1 / M} W^{l 1}) G_{j} ({z_{2}}^{1 / M} W^{l 2}) MS ({z_{1}}^{1 / M} W^{l 1}, {z_{2}}^{1 / M} W^{l 2}) - - - (13)

其中，0≤i，j≤M；W＝e^-jπ/M；H_i和G_i分别为分析滤波器组；P_ij(z₁，z₂)和Q_ij(z₁，z₂)为原图像在各频带上的子图像。

如果忽略降采样过程，那么(12)式和(13)式将简化为：

P_ij(z₁，z₂)＝H_i(z₁)Hj(z₂)PS(z₁，z₂)0≤i，j ≤M (14)

Q_ij(z₁，z₂)＝G_i(z₁)G_j(z₂)MS(z₁，z₂)0≤i，j ≤M (15)

5、用低空间分辨率多光谱图像的低频部分替换高空间分辨率全色遥感图像的低频部分。融合后低频子图像为：{Q₀₀}；融合后高(带)频子图像为：{P_ij，0≤i，j≤M且i+j＞0}；

6、对替代后的图像进行滤波器组重构以获得融合后的图像。

FZ (z_{1}, z_{2}) = Σ_{i = 0}^{M} Σ_{j = 0}^{M} F_{i} (z_{1}) F_{j} (z_{2}) X_{ij} (z_{1}, z_{2}), 0 \leq i, j \leq M - - - (16)

其中，

将本发明所得的融合结果，与其他融合方法所得的融合结果进行了对比，评价结果对照如表1～4所示，图3为融合所得图像。其中，(c)～(f)3～6通道滤波器组融合图像；(g)小波变换融合图像(两层分解)；(h)小波变换融合图像(三层分解)；(i)～(l)3～6通道滤波器组融合图像(忽略降采样)；(m)小波框架融合图像(两层分解)；(n)小波框架融合图像(三层分解)；结果表明，本发明在保留融合图像的光谱信息的同时，提高了图像的空间信息，能够使空间细节信息和光谱信息二者之间达到更好地的均衡，使得融合后图像的质量得以较大的提高，对于应用系统的后续处理具有重要意义和实用价值。

表1融合图像与多光谱图像间的平均偏差

	波段1	波段2	波段3	波段4
	波段1	波段2	波段3	波段4	Filterbanks3通道Filterbanks4通道Filterbanks5通道Filterbanks6通道DWF(2层)DWF(3层)	7.994810.510312.108813.327010.212214.8207	8.267710.832612.432213.581810.545214.9313	8.086810.600212.164513.345310.305814.7579	8.181810.739312.406813.668910.453715.2672

表2融合图像与全色图像间高频部分的相关系数

	波段1	波段2	波段3	波段4
	波段1	波段2	波段3	波段4	Filterbanks3通道Filterbanks4通道Filterbanks5通道Filterbanks6通道DWT(2层)DWT(3层)	8.758011.226212.747513.920911.106115.5235	9.044311.564613.025714.066411.387315.5279	8.855311.313712.767913.984711.176915.4085	8.950511.471913.056814.283011.329315.9421

表3融合图像与多光谱图像间的平均偏差(滤波器组忽略降采样过程)

	波段1	波段2	波段3	波段4
	波段1	波段2	波段3	波段4	Filterbanks3通道Filterbanks4通道Filterbanks5通道Filterbanks6通道DWF(2层)DWF(3层)	0.93810.98900.99120.99600.98190.9989	0.93940.98610.99210.99650.98290.9991	0.93950.98980.99180.99630.98260.9990	0.93750.98810.99080.99580.98120.9989

表4融合图像与全色图像间高频部分的相关系数(滤波器组忽略降采样过程)

	波段1	波段2	波段3	波段4
	波段1	波段2	波段3	波段4	Filterbanks3通道Filterbanks4通道Filterbanks5通道Filterbanks6通道DWT(2层)DWT(3层)	0.91270.97100.98750.99400.96690.9981	0.91380.97260.98850.99470.96800.9985	0.91410.97220.98810.99440.96800.9983	0.91200.96970.98690.99370.96580.9981

Claims

1、一种基于滤波器组的图像融合方法，其特征在于包括如下具体步骤：

1)在对同一地貌的低空间分辨率多光谱图像和高空间分辨率全色遥感图像进行配准的基础上，对低空间分辨率多光谱图像进行插值重采样，以实现二者具有相同的尺寸大小；

2)对高空间分辨率全色遥感图像和低空间分辨率多光谱图像进行直方图匹配，使它们的灰度均值和方差相一致；

3)设计满足融合要求的滤波器组，使产生的融合信号的低频部分具有低分辨率信号的低频信息，而高频部分具有高分辨率信号的高频信息，采用对原型滤波器进行余弦调制的方法来获得分析和综合滤波器组，当原型滤波器满足：

\{\begin{matrix} G (ω) = H (ω) & 0 \leq ω < π / M_{1} \\ 0 & π / M_{1} \leq ω \leq π \end{matrix}

使融合信号的频率信号失真最小，其中，G(ω)为用于分解高分辨率信号的原型滤波器；H(ω)为用于分解低分辨率信号的原型滤波器；M₁为通道数；

4)在得到满足融合要求的滤波器组的基础上，将高空间分辨率全色遥感图像和低空间分辨率多光谱图像分别进行滤波器组分解，得到多个频带的子图像，分解时保留降采样过程或忽略降采样过程；

5)用低空间分辨率多光谱图像的低频部分替换高空间分辨率全色遥感图像的低频部分，高空间分辨率全色遥感图像的高频部分保持不变；

6)对替代后的图像进行相应的滤波器组重构以获得融合后的图像。