CN103116881A - 基于PCA与Shearlet变换的遥感图像融合方法 - Google Patents

Abstract

本发明公开了一种基于PCA与Shearlet变换的遥感图像融合方法，主要解决多光谱和全色图像融合过程中光谱信息和空间分辨率难以平衡的问题。其实现步骤是：对上采样后的多光谱图像进行PCA变换，得到各个分量图像；计算每个分量图像与全色图像的相关系数，并计算其与最大相关系数的差值；对差值小于给定阈值的分量图像和全色图像分别进行Shearlet分解，根据分解结果得到融合后的分量图像；将融合后的分量图像和差值不小于给定阈值的分量图像组成数据集，对该数据集进行逆PCA变换，得到融合图像。本发明具有融合图像的光谱保持性和空间分辨率高的优点，可用于到军事目标识别、气象监测、环境监测、城市规划以及防灾减灾。

Description

基于PCA与Shearlet变换的遥感图像融合方法

技术领域

本发明属于智能图像处理领域，涉及图像融合方法，可用于到军事目标识别、气象监测、环境监测、土地利用、城市规划以及防灾减灾等多个领域的技术。

背景技术

随着遥感技术的发展，由各种卫星传感器对地观测获取同一地区的多源遥感图像数据，即多时相、多光谱、多传感器、多平台和多分辨率的遥感图像数据越来越多，为军事目标的识别、环境监测等提供了丰富而又宝贵的资料。但是由于实际应用中所需的遥感图像数据在时间、空间和光谱方面差异很大，而各种传感器提供的遥感图像数据都具有不同的特点，比如通过卫星获取的多光谱和全色图像，多光谱图像具有丰富的光谱信息，WordView-1卫星甚至能够获取到八个波段的多光谱图像，但是多光谱图像的分辨率一般情况下只有全色图像的四分之一，全色图像具有很高的空间分辨率，但是却没有丰富的光谱信息。所以遥感技术应用的主要障碍不是数据源的不足，而是从这些数据源中提取更丰富、更有用和更可靠信息的能力大小。各种单一的遥感手段获取的图像数据在几何、光谱和空间分辨率等方面存在明显的局限性和差异性，导致其应用能力受到限制，所以在实践中仅仅利用一种遥感图像数据是很难满足要求的，这就要求我们对多源图像数据之间的冗余性进行充分利用。冗余信息的应用，可以降低多源图像融合的误差和不确定性，提高识别率和精确度。多源遥感图像融合，尤其是多光谱和全色图像的融合，被认为是现代多源图像处理和分析中非常重要的一步。

目前，市场上使用的多光谱和全色图像融合方法主要有基于空间域的和基于变换域两种。

空间域中常用的融合方法有基于HIS变换、基于PCA变换和基于Gram-Schmidt变换三种。基于HIS变换的融合方法计算复杂度低，常被用于各种软件中，但是HIS变换融合后的图像光谱失真严重，并且只能用于三个波段的多光谱图像，随着传感器的飞速发展，从卫星传回地面的多光谱图像常常是多于三个波段的，所以基于HIS变换的融合方法的应用就受到了限制；基于PCA变换的融合方法得到的融合的图像具有比HIS变换具有更好的光谱特征和空间分辨率，对多光谱图像进行PCA变换分解，分解后的第一主分量图像含有原始多光谱图像的大部分能量，而当第一主分量图像和全色图像具有较高相关系数的时候，能够得到较高的融合图像，相反，当和全色图像相关系数最高的分量图像不是第一主分量图像的时候，如果仍采用传统的PCA方法，就会对融合图像的空间分辨率造成严重的影响，对具有多于三个波段的多光谱图像进行融合的时候常常会出现这种情况；基于Gram-Schmidt变换的融合方法已经被应用于ENVI软件中。在Gram-Schmidt变换方法的实现过程中，需要首先模拟一个低分辨率全色图像作为Gram-Schmidt分解的第一个分量，并基于此分量对多光谱图像进行正交分解，当模拟的第一分量图像与高分辨全色图像具有较高的相关系数的时候，能够得到较好的融合结果，反之，融合图像的光谱和分辨率就会受到较大的影响。目前，还没有一种更加有效的方式来解决Gram-Schmidt变换中的对第一分量的模拟问题。

基于变换域的常用的方法有基于小波的融合方法和基于多尺度几何分析的融合方法。基于小波的融合方法得到的融合图像能够较好的保持光谱信息，但是由于小波变换只能提取图像三个方向的特征，使得融合的图像空间分辨率不高。基于多尺度几何分析方法是目前图像融合领域研究的热点方法，在图像融合中，常用的多尺度几何分析工具主要有Bandlet、Contourlet和Shearlet。其中，Bandlet变换能够对图像进行任意方向的分解，按照实际的需要提取图像任意方向的特征，但是Bandlet变换不具有平移不变性，也就是当多光谱图像和全色图像在配准过程中存在配准误差的时候，融合图像就会出现双重边缘，严重影响融合图像的质量。Contourlet变换具有平移不变性，但是Contourlet变换在提取图像方向特征的时候，分解的方向数受到限制的，也就是只能对图像进行固定个方向的分解，使得图像的一些方向信息丢失，从而影响融合图像的空间分辨率。Shearlet变换不仅具有平移不变性，还能够对图像进行任意尺度和任意方向的分解，能够提取图像任意方向的特征，所以在融合图像中能够得到较高的空间分辨率。但是Shearlet变换，以及其它所有的多尺度几何分析工具，在多光谱和全色图像的融合过程中，需要对多光谱图像的每一个波段分别进行尺度的分解和方向特征的提取，随着多光谱图像的波段数越来越多，计算的时间复杂度也就会越来越大，单独使用多尺度几何分析工具来融合多光谱和全色图像，在计算的时间复杂度上，很难满足市场的要求。

发明内容

本发明的目的在于针对上述已有技术的不足，提出一种基于PCA与Shearlet变换的遥感图像融合方法，以解决现有技术在多光谱和全色图像融合过程中光谱信息和空间分辨率难以平衡，及多光谱和全色图像融合后光谱失真或者空间分辨率不高的问题，提高融合图像的质量。

为实现上述目的，本发明的实现步骤如下：

1)分别输入一副采样后的多光谱图像I₁和高分辨率全色图像I₂；

2)对多光谱图像I₁进行上采样，上采样后的多光谱图像的大小与高分辨率全色图像I₂的大小相同；

3)对上采样后的多光谱图像进行主成分分析PCA变换，得到变换后的各个分量图像PC_i,i＝1,2,...,Nb,其中Nb是多光谱图像的波段数；

4)分别计算每一个分量图像PC_i与全色图像I₂的相关系数值，记为CC_i,i＝1,2,...,N，并且将最大的相关系数值保存，记为CC_Max；

5)计算每一个相关系数值CC_i与最大相关系数值CC_Max的差值，对所有差值小于给定阈值的分量图像PC_S,S＝1,2,...,N₁，分别执行步骤6）到8），并将其余分量图像保存为PC_NS,NS＝1,2,...,N₂，其中N₁为相关系数值CC_i与最大相关系数值CC_Max的差值小于给定阈值的分量图像的个数，N₂为相关系数值CC_i与最大相关系数值CC_Max的差值大于或者等于给定阈值的分量图像的个数，且N₁+N₂＝Nb；

6)对相关系数差值小于给定阈值的分量图像PC_S,S＝1,2,...,N₁和高分辨率全色图像I₂分别进行Shearlet变换分解，分别得到分量图像的一个低频系数L_k和多个方向子带系数H_k，以及全色图像的一个低频系数L_p和多个方向子带系数H_p；

7)对得到的分量图像PC_S,S＝1,2,...,N₁和全色图像I₂的低频系数L_k，L_p和多个方向子带系数H_k，H_p分别采用不同的提取规则来作为融合后的分量图像的低频系数L_f和多个方向子带系数H_f；

8)对融合后分量图像的低频系数L_f和多个方向子带系数H_f进行逆Shearlet变换，得到融合后的分量图像I_f；

9)对所有融合后的分量图像I_f和步骤5）保存的其余分量图像PC_NS,NS＝1,2,...,N₂进行逆主成分分析PCA变换，得到融合图像I_fus。

本发明与现有技术相比，具有如下效果：

(a)本发明由于使用自适应的主成分分析PCA算法对多光谱图像进行分解，分解后的分量图像的相关系数与全色图像的相关系数两个共同作用，来确定参与后续融合处理的分量图像，克服了传统的只采用PCA分解的第一分量图像参与后续图像的融合处理从而易导致融合图像光谱扭曲的问题，这样就提高了融合图像的光谱信息。

(b)本发明由于使用Shearlet变换分解分量图像与全色图像，利用Shearlet变换能够对图像进行任意方向特征提取的特点，克服了传统小波变换只能提取图像三个方向的特征，易导致图像细节信息丢失的问题，从而提高了融合图像的分辨率。

附图说明

图1是本发明的多光谱图像和全色图像融合流程图；

图2是本发明仿真使用的QuickBrid卫星图像的假彩色多光谱和全色图像；

图3是本发明对图2中的多光谱和全色图像融合的结果图像；

图4是本发明仿真使用的QuickBird卫星图像的真彩色多光谱和全色图像；

图5是本发明对图4中的多光谱和全色图像融合的结果图像。

具体实施方式

参照图1，本发明结合自适应PCA变换和Shearlet变换的多光谱和全色图像融合方法的实现步骤如下：

步骤1，对多光谱图像分别进行上采样采样后的多光谱图像的大小与对应的全色图像图的大小相同，将采样后的多光谱图像记作I₁，原始的全色图像记作I₂。

步骤2，对步骤1中采样后的多光谱图像I₁进行主成分分析PCA变换，得到变换后的各个分量图像，记作PC_i,i＝1,2,...,Nb，其中Nb是多光谱图像的波段数，分量图像的大小和采样后的多光谱图像I₁的大小相同，大小记作M×N，M和N分别是分量图像的宽和高。

由于原始多光谱图像的维数较高，处理的时间复杂度高，并且各个波段之间的数据具有极大的相关性，使得噪声信息难以去除，通过对多光谱图像进行主成分分析PCA变换，将多光谱图像分解成多个分量图像，这些分量图像互不相关，从而达到隔离噪声和减少数据集维数的目的。

步骤3，对步骤2中得到的每一个分量图像PC_i，计算其与全色图像I₂的相关系数值，记作CC_i,i＝1,2,...,N，并且将最大的相关系数值记录下来，记作CC_Max,其中相关系数的计算公式如下：

${\mathrm{CC}}_i=\frac{\underset{p=1}{\overset{M}{\mathrm{\Σ}}}\underset{q=1}{\overset{N}{\mathrm{\Σ}}}({\mathrm{PC}}_i(p,q)-{\mathrm{AP}}_c)(I_p(p,q)-{\mathrm{AI}}_p)}{\sqrt{\underset{p=1}{\overset{M}{\mathrm{\Σ}}}\underset{q=1}{\overset{N}{\mathrm{\Σ}}}{({\mathrm{PC}}_i(p,q)-{\mathrm{AP}}_c)}^2\underset{p=1}{\overset{M}{\mathrm{\Σ}}}\underset{q=1}{\overset{N}{\mathrm{\Σ}}}{(I_p(p,q)-{\mathrm{AI}}_p)}^2}}---\left(A\right)$

其中PC_i(p,q)和I_P(p,q)分别表示分量图像PC_i和全色图像I₂在坐标(p,q)点处的灰度值，AP_c和AI_p分别表示分量图像PC_i和全色图像I₂的灰度平均值。

步骤4，计算步骤3中得到的每一个相关系数值CC_i与最大相关系数值CC_Max的差值，即:

Diff_i＝|CC_i-CC_Max|    (B)将所有差值小于阈值T的分量图像分到第一组，记作PC_S,S＝1,2,...,N₁；将其它的分量图像分到第二组，记作PC_NS,NS＝1,2,...,N₂，其中N₁为相关系数值CC_i与最大相关系数值CC_Max的差值小于给定阈值的分量图像的个数，N₂为相关系数值CC_i与最大相关系数值CC_Max的差值大于或者等于给定阈值的分量图像的个数，且N₁+N₂＝Nb，对第一组中的每一个分量图像，执行步骤5到步骤8。

步骤5，对步骤4中得到的第一个分量图像PC_k,k∈S和全色图像I₂分别进行Shearlet变换分解，将分量图像PC_k,k∈S和全色图像I₂分别分解为一个低频图像和多个方向子带图像。将分量图像的低频系数和多个方向子带系数分别记作L_k和H_k，全色图像的低频系数和多个方向子带系数分别记作L_p和H_p。

其中Shearlet分解分量图像和全色图像的步骤如下：

5a）对分量图像PC_k,k∈S进行多尺度分解。本发明采用拉普拉斯金字塔变换对分量图像进行了4个尺度的分解，即通过拉普拉斯金字塔变换，将分量图像分解为一个低频系数和4个高频系数H1_k,H2_k,H3_k,H4_k的尺度图IS_m,m＝1,2,3,4；

5b）对步骤5a）得到的4个高频系数，分别进行多方向分解，以提取图像在每个尺度下的方向特征。本发明采用Shear滤波器分别对每个尺度图像IS_m,m＝1,2,3,4进行方向分解，从粗尺度到细尺度，四个尺度图像IS_m分别依次被分解为6个,6个,10个和10个方向，从而得到分量图像的多个方向子带系数H_k；

5c）重复步骤5a）和5b），将分量图像替换为全色图像I₂，得到全色图像Shearlet分解后的低频系数L_p和多个方向子带系数H_p。

步骤6，计算分量图像PC_k的灰度方差与全色图像I₂的灰度方差，以及它们的比值，将比值记作R_k，比值的计算公式如下：

$R_k=\frac{\mathrm{\σ}\left({\mathrm{PC}}_k\right)}{\mathrm{\σ}\left(I_2\right)}---\left(C\right)$

其中σ表示图像的灰度方差，其计算公式为：

$\mathrm{\σ}\left({\mathrm{PC}}_k\right)=\frac{\underset{p=1}{\overset{M}{\mathrm{\Σ}}}\underset{q=1}{\overset{N}{\mathrm{\Σ}}}{({\mathrm{PC}}_k(p,q)-{\mathrm{AP}}_c)}^2}{M\×N}---\left(D\right)$

$\mathrm{\σ}\left(I_2\right)=\frac{\underset{p=1}{\overset{M}{\mathrm{\Σ}}}\underset{q=1}{\overset{N}{\mathrm{\Σ}}}{(I_p(p,q)-{\mathrm{AI}}_p)}^2}{M\×N}---\left(E\right)$

其中PC_k(p,q)和I_p(p,q)表示图像分量PC_k和全色图像I₂在坐标(p,q)点处的灰度值，AP_c和AI_p分别表示分量图像PC_k和全色图像I₂的灰度平均值。

步骤7，对步骤5得到的全色图像的每一个方向子带系数，均乘以权值w，得到融合后的分量图像的方向子带系数；将分量图像PC_k的低频系数，作为融合后的分量图像的低频系数，其中权值w等于1+R_k。

步骤8，对步骤7得到的融合后的分量图像的低频系数和方向子带系数进行逆Shearlet变换，得到融合后的分量图像。

步骤9，对第一组中分量图像PC_S,S＝1,2,...,N₁的每一个分量图像执行步骤5到步骤8，均得到对应的融合后的分量图像I_f，将所有融合的分量图像I_f和步骤4中保存在第二组分量图像PC_NS,NS＝1,2,...,N₂组成新的图像集，并对新的图像集进行逆PCA变换，得到融合后的图像I_fus。

本发明的效果可以用下列的仿真实验进一步说明:

(1)仿真条件

本发明的仿真的硬件条件为：windows XP,SPI,CPU Pentium(R)4,基本频率2.4GHZ,软件平台为：MatlabR2009a，仿真选用的图片来源是QuickBrid卫星图像的假彩色多光谱和全色图像和QuickBrid卫星图像的真彩色多光谱和全色图像，分别对应图2和图4,其中图2(a)是本发明仿真使用的QuickBrid卫星图像的假彩色多光谱，图2(b)是本发明仿真使用的QuickBrid卫星图像的全色图像；图4(a)是本发明仿真使用的QuickBrid卫星图像的真彩色多光谱图像，图4(b)为本发明仿真使用的QuickBrid卫星图像的是全色图像。

(2)仿真内容及结果

仿真1，将图2(a)所示的QuickBrid卫星图像的假彩色多光谱和图2(b)所示的QuickBrid卫星图像的全色图像分别作为多光谱图像I₁和全色图像I₂，然后按照本发明上述具体实施方案中的步骤进行仿真处理；并且分别应用传统的结合wavelet和PCA的融合方法和传统的结合Contourlet和PCA融合方法对图2(a)、图2(b)进行融合，融合结果参照图3，其中3(a)是结合wavelet和PCA的融合结果图，3(b)是结合Contourlet和PCA的融合结果图，3(c)是本发明融合的结果图，3(d)是3(a)图的局部放大图，3(e)是3(b)图的局部放大图，3(f)是3(c)图的局部放大图，分别对比结果图和局部放大图可见，本发明融合的结果图更清晰，图像质量更高。

仿真2，将图4(a)所示的QuickBrid卫星图像的真彩色多光谱图像和图4(b)所示的QuickBrid卫星图像的全色图像分别作为多光谱图像I₁和全色图像I₂，然后按照上述具体实施方式中的步骤对I₁和I₂进行仿真处理；并分别应用传统的结合wavelet和PCA的融合方法和传统的结合Contourlet和PCA融合方法对图4(a)和4(b)进行融合，融合结果参照图5，其中5(a)是结合wavelet和PCA的融合结果图，5(b)是结合Contourlet和PCA融合结果图，5(c)是本发明融合的结果图，5(d)是5(a)图的局部放大图，5(e)是5(b)图的局部放大图，5(f)是5(c)图的局部放大图，分别对比结果图和局部放大图可见本发明融合的结果图更清晰，图像质量更高。

以上实验结果表明：本发明相比现有技术在解决多光谱和全色图像融合过程中光谱信息和空间分辨率难以平衡，及多光谱和全色图像融合后光谱失真或者空间分辨率不高的问题上，具有融合图像的光谱保持性和空间分辨率高的优点，提高了融合图像的质量。

Claims (3)

1.一种基于PCA与Shearlet变换的遥感图像融合方法，包括如下步骤：

1)分别输入一副多光谱图像I₁和高分辨率全色图像I₂；

2)对多光谱图像I₁进行上采样，上采样后的多光谱图像的大小与高分辨率全色图像I₂的大小相同；

3)对上采样后的多光谱图像进行主成分分析PCA变换，得到变换后的各个分量图像PC_i,i＝1,2,...,Nb,其中Nb是多光谱图像的波段数；

4)分别计算每一个分量图像PC_i与全色图像I₂的相关系数值，记为CC_i,i＝1,2,...,N，并且将最大的相关系数值保存，记为CC_Max；

5)计算每一个相关系数值CC_i与最大相关系数值CC_Max的差值，对所有差值小于给定阈值的分量图像PC_S,S＝1,2,...,N₁，分别执行步骤6）到8），并将其余分量图像保存为PC_NS,NS＝1,2,...,N₂，其中N₁为相关系数值CC_i与最大相关系数值CC_Max的差值小于给定阈值的分量图像的个数，N₂为相关系数值CC_i与最大相关系数值CC_Max的差值大于或者等于给定阈值的分量图像的个数，且N₁+N₂＝Nb；

6)对相关系数差值小于给定阈值的分量图像PC_S,S＝1,2,...,N₁和高分辨率全色图像I₂分别进行Shearlet变换分解，分别得到分量图像的一个低频系数L_k和多个方向子带系数H_k，以及全色图像的一个低频系数L_p和多个方向子带系数H_p；

7)对得到的分量图像PC_S,S＝1,2,...,N₁和全色图像I₂的低频系数L_k，L_p和多个方向子带系数H_k，H_p分别采用不同的提取规则来作为融合后的分量图像的低频系数L_f和多个方向子带系数H_f；

8)对融合后分量图像的低频系数L_f和多个方向子带系数H_f进行逆Shearlet变换，得到融合后的分量图像I_f；

9)对所有融合后的分量图像I_f和步骤5）中保存的其余分量图像PC_NS,NS＝1,2,...,N₂进行逆主成分分析PCA变换，得到融合图像I_fus。

2.根据权利要求1所述的结合自适应PCA和Shearlet变换的多光谱图像和全色图像融合方法，其中步骤6）所述的对相关系数小于给定阈值的分量图像PC_S,S＝1,2,...,N₁和全色图像I₂分别进行Shearlet分解，按如下步骤进行：

对分量图像PC_S,S＝1,2,...,N₁进行4个尺度的分解，即通过拉普拉斯金字塔变换，将分量图像分解为一个低频系数L_k和4个高频系数H1_k,H2_k,H3_k,H4_k的尺度图IS_m,m＝1,2,3,4；

对步骤6a）得到的4个高频系数H1_k,H2_k,H3_k,H4_k的尺度图IS_m分别进行多方向分解，即从粗尺度到细尺度，四个尺度图IS_m分别依次被分解为6个、6个、10个和10个方向，从而得到分量图像PC_S,S＝1,2,...,N₁的多个方向子带系数H_k；

重复步骤6a）和步骤6b），将分量图像替换为全色图像I₂，得到全色图像Shearlet分解后的低频系数L_p和多个方向子带系数H_p。

3.根据权利要求1所述的结合自适应PCA和Shearlet变换的多光谱和全色图像融合方法，其中步骤7）所述的对分量图像PC_S,S＝1,2,...,N₁和全色图像I₂的低频系数和多个方向子带系数，分别采用不同的提取规则来分别作为融合后分量图像的低频系数和方向子带系数，按如下步骤进行：

对于分量图像PC_S,S＝1,2,...,N₁，则提取该分量图像的的低频系数L_k作为融合后的分量图像的低频系数L_f；

对于全色图像I₂，则提取该全色图像的每一个方向子带系数H_p，用该方向子带系数H_p与权值W的乘积作为融合后的分量图像的方向子带系数H_f＝H_p×W，其中权值W根据分量图像与全色图像的灰度方差比R_k确定，R_k是分量图像PC_k的灰度方差与全色图像I₂的灰度方差比。

Images