遥感图像双线性插值小波变换融合的优化方法

摘要

一种遥感图像双线性插值小波变换融合的优化方法，采用双线性插值算式对高光谱图像进行重采样优化计算，使其空间分辨率等于全色图像，根据小波分解理论对全色图像与重采样后的高光谱图像进行小波分解，获取各自的近似低频和细节高频分量。遥感图像的小波系数融合时，将高光谱图像的小波分解的近似分量作为新图像的近似分量，而细节分量融合，则利用线性计算来获取，根据小波系数融合所得新图像数据块的小波系数，进行小波重建和数据块拼接，生成所需的遥感融合图像。本发明在保留高光谱图像的光谱信息的同时，能够提高其空间分辨率，同时显著提高了计算机处理速度，获得的融合图像质量更加可靠。

说明书

技术领域：

本发明涉及一种遥感图像双线性插值小波变换融合的优化方法，是高光谱遥感图像经双线性重采样后，与高空间分辨率全色图像融合时进行高精度及快速计算的融合优化方法，是多源遥感影像后处理阶段的一项核心技术，在各类军用或民用的遥感信息处理系统、数字城市空间信息系统等领域中均可有广泛的应用。

背景技术：

二十多年来，遥感技术获得了迅猛发展，各种面向复杂应用背景的多平台、多时相、多光谱和高分辨率卫星遥感系统也随之大量涌现，它广泛应用于地表植被的分类和环境观测等领域。为了能进行更准确地分类决策或增加定位精度，对高空间与光谱分辨率图像的需求非常迫切。欲实现这个目标，有两条途径：一方面，采用高空间与光谱分辨率的遥感器，例如法国SPOT卫星上遥感器的XS1、XS3波段的地面分辨率为10m×10m；其IPS的PAN图像地面分辨率为5.8m×5.8m。IPS的1CLISS III则为23.5×23.5m。美国IKONOS卫星的全色(PAN)图像地面分辨率为1m×1m、波长0.45-0.90μm；多光谱图像则为4m×4m、其波长范围为蓝0.45-0.52μm、绿0.52-0.60μm、红0.63-0.69μm、近红外0.69-0.90μm，这也与陆地卫星LAND SAT的1-4波段的光谱相同。而另一方面，则利用图像处理技术，即进行不同光谱与空间分辨率遥感图像的融合，或称多光谱图像空间分辨率增强(Enhancement of the Spatial resolution ofMulti-spectral Image)、多光谱图像的锐化(Sharpening)，可节省大量的研究经费。因此，基于多传感器的信息融合理论，多源遥感影像信息融合(Remote Sensing ImageInformation Fusion)研究便于20世纪80年代应运而生。

一般地说，遥感图像融合就是图像合成技术，将不同平台(卫星与机载)上的同一或不同传感器获取的不同空间与光谱分辨率图像按特定的算法进行处理，以使所产生的新图像同时具有原来图像的多光谱特性以及高地面分辨率，来实现不同的应用需求。目前，应用较多的遥感图像融合是不同分辨率的多光谱(Multi-Spectral)与PAN(Panchromatic)图像融合或高光谱(Hyper-Spectral)图像与PAN图像的融合，主要方法有IHS变换、主成分分析(PCA)、BROVERY变换以及小波变换合等。

对这些融合算法进行归纳，它们的一个共同特征是，融合前均将高光谱图像重采样为与全色图像空间分辨率相等的新图像，以供进一步处理。为此，面临两个方面的问题需要解决：

1.高光谱图像重采样算法对融合图像的光谱及空间分辨率影响；

2.由于高光谱图像的大小一般为几百兆的数量级，若其重采样算法变得复杂时，须考虑计算机的运算速度。

目前，进行上述遥感图像融合时，在精度许可的条件下，主要应用近邻插值方法，以满足算法的处理速度要求。这里采用的近邻插值算法的原理是，像素点插值搜索过程中，若高光谱图像中某素点离重采样算法所得的新图像中某像素点的空间位置最近，则将其像素值作为新图像对应点的像素值。该方法通常情况下可以满足图像处理要求。但若要求光谱信息损失较小时，则它难以达到所需的融合精度，究其原因，是由于近邻算法本身的局限性造成的，所以需研究新的算法。

发明内容：

本发明的目的在于针对现有技术的上述不足，提供一种遥感图像双线性插值小波变换融合的优化方法，不仅能够最有效的保持原始高光谱图像的光谱信息，而且还能在提高光谱遥感图像空间分辨率的同时，使计算机运算速度加快，满足应用需求。

为实现这样的目的，本发明的技术方案中，利用遥感图像处理与开发平台提供的开发语言，针对高光谱图像的文件大小及计算机的处理速度，采用双线性插值算式对高光谱图像进行重采样优化计算，使其空间分辨率等于全色图像。然后，根据小波分解理论对全色图像与重采样后的高光谱图像进行小波分解，获取各自的近似低频和细节高频分量。遥感图像的小波系数融合时，将高光谱图像的小波分解的近似分量作为新图像的近似分量，而细节分量融合，则利用线性计算来获取，这样有利于保留高光谱图像的光谱信息。最后，根据小波系数融合所得新图像数据块的小波系数，进行小波重建和数据块拼接，生成所需的遥感融合图像。

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

1.数据特性的提取：根据应用的需求，获得配准后待融合的高光谱和全色遥感

  图像文件，通过遥感图像处理与开发平台提供的接口函数，取得它们的特性

  数据，包括文件大小，空间分辨率以及图像高光谱遥感图像的波段数等。

2.高光谱图像的双线性重采样：以全色图像的空间分辨率为标准，获取计算机

  的内存与CPU处理速度等指标后，根据高光谱文件对应的数据量的大小将其

  分割为几个连续的数据块(数目为计算机内存量/高光谱图像数据量)，按双线

  性插值规则，对各数据块分别逐点进行搜索、插值，生成具有与全色图像空

  间分辨率相同的新高光谱图像数据块。

3.全色图像的数据分块，按照第二步对应的数据分块数，将全色图像分成与经

  双线性插值后所得的新高光谱图像数据块数目相等的数据块。

4.循环调入对应的全色与高光谱图像数据块，进行各数据块的二维小波分解。

  就是在水平和垂直方向对全色与高光谱图像数据块各自独立地进行一次一维

  小波变换。相当于先在水平方向对图像分别进行低通和高通滤波，并作下采

  样，然后在垂直方向上实施同样的操作。最终获得各个尺度下的低频近似分

  量和三个方向(水平、垂直、斜线)的高频细节分量，它们构成对应数据块的小

  波系数。

5.高光谱与全色图像各数据块的小波系数融合：将高光谱图像数据块的低频近

  似分量作为融合图像数据块的低频分量，而融合图像数据块的高频分量，则

  使用线性加权算法由全色与高光谱图像数据块的高频细节分量取得，由此构

  建新图像的小波系数(低频与高频分量)。

6.融合图像数据块的小波变换重建及数据块的拼接：将获得的融合图像数据块

  的小波系数按小波变换的重构原理，进行图像重构，获得所需的高空间分辨

  率和光谱分辨率的数据块，然后将这些数据块按前述高光谱图像的分割顺序

    进行拼接，组成新图像。

本发明的特色在于，由于同一地区高光谱与全色遥感图像的空间分辨率及光谱分辨率等特性上相差很大，按双线性规则对高光谱图像进行重采样后，利用小波变换进行融合，在保留高光谱图像的光谱信息的同时，能够提高其空间分辨率，使地物空间结构和纹理特性能更加明确的显示。同时，由于进行了数据分块，将各数据块依次调入计算机内存，能显著地提高计算机处理速度。总之，本发明提出的方法可广泛地用于遥感图像数据的融合处理，获得更加可靠的融合图像质量。

附图说明：

图1为本发明所提出处理方法的总体框图。

如图1所示，首先利用基于Windows的ENVI/IDL遥感图像处理与开发平台提供的接口函数来获取高光谱与全色图像的特性数据。将高光谱图像按一定的比例进行数据分块后对各块再实施双线性重采样，使各数据块的空间分辨率等于全色图像的空间分辨率。此后，按相同的比例对全色图像进行数据分块。循环调入全色与重采样后高光谱图像数据块，对各对应的数据块进行小波分解，获得小波分解的低频近似分量与高频细节分量。循环调入小波分解后的数据块进行小波系数融合，构造新图像数据块的小波系数。最后，对新图像数据块进行小波变换重建和数据块的拼接，得到所需的融合图像。

图2为双线性重采样算法所使用的插值模型。

图3是全色和高光谱图像对应数据块小波变换融合的计算过程示意图。其中C_a^p、C_h^p、C_v^p、C_d^p分别为全色图像数据块小波分解后的低频分量、高频水平、高频垂直和高频斜线分量；

C_a^m、C_h^m、C_v^m、C_d^m分别为高光谱图像重采样后数据块小波分解后的低频分量、高频水平、高频垂直和高频斜线分量。

图4为本发明采用的方法与其它方法的融合效果对比。

其中，图4(a)的Land TM图像为高光谱图像，本实施例取其波段4、3、2，文件大小1007×560，空间分辨率28m。图4(b)的SPOT遥感影像是全色图像，文件大小2820×1569，窨分辨率10m。图4(c)是近邻插值小波融合效果，图4(d)是双线性插值小波融合效果，从中可以看出，双线性插值小波变换融合有利于保留高光谱图像的光谱信息。

具体实施方式：

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

图1为本发明提出的一种用于高光谱与全色图像双线性插值小波变换融合优化方法的总体框图，数据文件的输入是图4(a)的Land TM图像为高光谱图像，取其波段4、3、2，文件大小1007×560，空间分辨率28m。图4(b)的SPOT遥感影像是全色图像，文件大小2820×1569，空间分辨率10m。各部分具体实施细节如下：

1.数据特性的提取，获得高光谱与全色图像的特性数据利用ENVI/IDL开发平台的ENVI DISPLAY BANDS函数，可以显示出文件的大小、波段数及空间分分辨率等。

2.高光谱图像的双线性重采样，利用ENVI/IDLIDL的C语言接口调用WINDOWS2000操作系统提供的API函数，获取计算机的硬件特性，如CPU处理速度、内存大小后，对高光谱图像按图2所示的方法进行双线性重采样以及按内存/高光谱文件数据量的比例进行高光谱数据文件的分块。

对于高光谱图像中某像素点(x，y)及其周围的点(0，0)、(0，1)、(1，0)、(1，1)，它们对应的图像灰度值分别为f(x，y)、f(0，0)、f(0，1)、f(1，0)和f(1，1)，则插值过程是：

首先，通过一阶线性插值得出f(x，0)：

      f(x，0)＝f(0，0)+x[f(1，0)-f(0，0)]                    (1)

类似地，对f(x，1)进行一阶线性插值：

      f(x，1)＝f(0，1)+x[f(1，1)-f(0，1)]                     (2)

最后，对垂直方向进行一阶线性插值，以确定f(x，y)

      f(x，y)＝f(x，0)+x[f(x，1)-f(x，0)]                     (3)

3.全色图像的数据分块

按第二步中获得的计算机内存量/高光谱图像数据量的数值，对全色图像分为与高光谱图像数据块数目相等的数据块。

4.循环调入对应的全色与高光谱图像数据块，依次进行各数据块的的小波分解。

图3提供了一个简洁的各数据块小波融合层次框架，其中全色与高光谱数据块被分成两个部分：近似分量(Approximate)和细节分量(Detail)，而后者又可分为水平、垂直和斜线分量，计算式是：

设ψ(t)∈L²(R)(L²(R)表示平方可积的实数空间，即能量有限的信号空间)，其傅里叶变换为ψ^*(ω)，当ψ^*(ω)满足允许条件： _a，b(t)： $>>>\psi >>a>,>b>>>>(>t>)>>=>>1>>|>a>|>>>\psi >>(>>>t>->b>>a>>)>>>s>a，b\in R；a\ne 0 (5)式中a为尺度因子，b为平移因子。$

对于任意函数f(t)∈L²(R)的连续小波变换W_f(a，b)有 $>>>W>f>>>(>a>,>b>)>>=>>>|>a>|>>>->1>/>2>>>>\int >R>>f>>(>t>)>ver>>>>\psi >*>>>(>>>t>->b>>a>>)>>>‾>>dt>->->->>(>6>)>>>s>式中，带符号“-”表示复数共轭。$

目前，已经形成了多种实用的小波族，它们各有长处，在本方案中，选用Symlet小波族。

5.高光谱与全色图像各数据块小波系数融合

将上述依次调入计算机内存的各数据块小波分解后，进行小波系数融合。设融合图像数据块的小波系数的低频分量为C_a^*，高频分量分别为C_h^*、C_v^*、C_d^*，则，新图像数据块小波系数采用的计算式为：式中a_x、b_x分别是在数据块的像素点中，从C_x^P到C_x^M直线的斜率与截距，而x则可取h、v、d，m表示高光谱图像，p表示全色图像。

6.融合图像数据块的小波变换重建及数据块拼接

融合图像数据块的像素值看作为函数f(x，y)，则f(x，y)∈V_j²(j∈Z)，它在V_j²空间的投影为A_jf(x，y)，这里V_j²是L²(R)的一个多尺度空间剖分，于是有

A_jf(x，y)＝A_j+1f(x，y)+D_j+1¹f(x，y)+D_j+1²f(x，y)+D_j+1³f(x，y)       (8)式中 $>>>A>>j>+>1>>>f>=>>\Sigma >>>m>1>>=>->\infty >>\infty >>>\Sigma >>>m>2>>=>->\infty >>\infty >>>C>>j>+>1>>>>(>>m>1>>,>>m>2>>)>>>\psi >>j>+>1>>>>(>>m>1>>,>>m>2>>)>>>s>$ >sup>>D>>j>+>1>>isup>>f>=>>\Sigma >>>m>1>>=>->\infty >>\infty >>>\Sigma >>>m>2>>=>->\infty >>\infty >sup>>D>>j>+>1>>isup>>>(>>m>1>>,>>m>2>>)>>>\psi >>j>+>1>>>>(>>m>1>>,>>m>2>>)>>>s>$$

i＝1，2，3

C_j+1(m₁，m₂)——融合图像数据块在二进制网格2^j+1尺度下的离散逼近(低频轮廓)

D_j+1ⁱ(m₁，m₂)——融合图像数据块在二进制网格2^j+1尺度下的离散细节(高频纹理)，分别对应于水平、垂直和斜线分量。

(m₁，m₂)——表示融合图像数据块像素的位置

用H_r、H_c和G_r、G_c分别表示镜像共轭滤波器，H、G作用于阵列{C_j(m₁，m₂)，(m₁，m₂)∈Z²}的行和列，则本方案使用的二阶Mallat分解计算式为相应地，其重构公式是： $>>>C>j>>=ver>>>H>r>>‾>ver>>>H>c>>‾>>>C>>j>+>1>>ver>>>H>r>>‾>ver>>>G>c>>‾>sup>>D>>j>+>1>>1sup>>+ver>>>G>r>>‾>ver>>>H>c>>‾>sup>>D>>j>+>1>>2sup>>+ver>>>G>r>>‾>ver>>>G>c>>‾>sup>>D>>j>+>1>>3sup>>->->->>(>10>)>>>s>式中 H、G分别为H、G的共轭。$

融合图像数据块小波重建后，按照原先的高光谱图像分割的顺序进行数据拼接，即可得到所需的融合图像。