CN101916436B - 一种多尺度空间投影遥感影像融合方法 - Google Patents

Abstract

本发明涉及一种多尺度空间投影遥感影像融合方法，包括以下步骤：利用原始低空间分辨率的多光谱影像计算其亮度影像，通过改进的高斯尺度空间理论，分别生成该亮度影像和原始高空间分辨率全色影像的高斯影像立方体，然后利用高斯影像立方体的层间差异提取亮度影像和全色影像的空间细节特征，最后将空间细节特征按照加权融合策略投影到原始多光谱影像上，得到具有高空间分辨率和高光谱分辨率的融合多光谱影像。本发明在提高融合影像的空间细节表现能力的同时，很好地保留了原多光谱图像的光谱信息，其物理意义明确，实现结构简单，融合效果好。

Description

一种多尺度空间投影遥感影像融合方法

技术领域

本发明属于遥感图像处理数据融合技术领域，涉及一种多尺度空间投影遥感影像融合方法。

背景技术

遥感技术的发展提供了丰富的多源遥感数据。这些来自不同传感器的数据具有不同的时间、空间和光谱分辨率以及不同的极化方式。单一传感器获取的图像信息量有限，往往难以满足应用需要，通过图像融合可以从不同的遥感图像中获得更多的有用信息，补充单一传感器的不足。图像融合是指将多源遥感图像按照一定的算法，在规定的地理坐标系，生成新的图像的过程。全色图像一般具有较高空间分辨率，多光谱图像光谱信息较丰富，为提高SPOT多光谱图像的空间分辨率，可以将全色图像融合进多光谱图像。通过融合既提高多光谱图像空间分辨率，又保留其多光谱特性。

传统的遥感影像融合算法如IHS变换、PCA变换在进行遥感影像融合时，能够获得具有较高空间质量的融合影像，但也有一定的局限性。IHS变换算法简单，能获得较高的空间分辨率，但要求替换的亮度分量两者之间有较大的相关性，只能适用于三个波段，对于多波段影像显得无能为力，而且由于全色影像与亮度影像之间的巨大差异导致色度改变，会产生严重的光谱失真。PCA变换能够对任意多波段影像进行融合，但是会使融合影像失去原有的物理特性，而且与IHS变换类似，融合效果取决于替换的两个主分量之间的相关性。

近年来多尺度理论在图像处理领域得到越来越多的关注。众所周知，人在不同的距离下，观测同一图像获得的感受是不一样的，如远距离看到的是图像轮廓，近距离下看到的是更多细节，这就是尺度效应。而尺度空间技术是从原始图像出发导出一系列越来越平滑、越来越简化的图像，随着尺度的增加，输出的图像越来越模糊，细节信息越来越不明显。多尺度表达的基本思想就是对原始信号嵌入一簇尺度参数变化的信号，大尺度下信号的结构应该是小尺度下信号结构的简化。“尺度空间”概念是著名的图像金字塔概念的新补充，最先在Kelly的图像处理中被采用，后来又以各种形式进行扩展。在现代的尺度空间表示中，每一层是前一层的模糊结果。通常每一层通过一个方差为σ的高斯函数卷积原图像而产生，而每一层的σ都不一样，这个方差就是“尺度参数”，尺度参数越大，图像越模糊，越多的细节被丢弃，这些细节对人眼视觉起关键作用，利用它可以有效地进行目标识别、特征提取等。

发明内容

针对现有技术的缺点，本发明的目的是提供一种既能提高其空间分辨率，又能很好保留遥感图像的光谱信息的多尺度空间投影融合方法。

为实现上述目的，本发明的技术方案为一种多尺度空间投影遥感影像融合方法，包括以下步骤：

步骤a，将相同地域的原始低空间分辨率的多光谱影像采样到原始高空间分辨率全色影像同等像素大小，进行精确几何校正处理，然后利用经过采样和校正后的多光谱影像计算其亮度影像；

步骤b，分别生成原始高空间分辨率全色影像和步骤a所得亮度影像的高斯影像立方体；所述高斯影像立方体由若干层的影像组成，最底层影像即为原始高空间分辨率全色影像和步骤a所得亮度影像，当前层影像是由上一层影像与体现当前层尺度特征的高斯核函数卷积而成，设最底层的尺度参数为σ，从最底层影像到最顶层影像的尺度参数由σ变化到2σ，并且当前层影像的尺度参数与上一层影像的尺度参数满足倍增关系；

步骤c，利用步骤b所得高斯影像立方体的层间差异，提取亮度影像和原始高空间分辨率全色影像的空间细节特征；

步骤d，将步骤c提取出的空间细节特征按照加权融合策略投影到原始低空间分辨率的多光谱影像上，得到具有高空间分辨率和高光谱分辨率的融合多光谱影像。

而且，步骤a中，采样过程采取最邻近像元法、双线性内插法或双三次卷积法方法。

而且，步骤a中，利用经过采样和校正后的多光谱影像计算其亮度影像的具体方式为，对多光谱影像的所有单波段影像求平均亮度影像。

而且，步骤b中，高斯影像立方体中的各层影像满足以下公式(1)：

$\left\{\begin{array}{c}{I}_p=G(x,y;{\mathrm{\σ}}_p)*I_{p-1}\\ G(x,y;{\mathrm{\σ}}_p)=\frac{1}{2\mathrm{\π}{\mathrm{\σ}}^2}{e}^{-\frac{x^2+y^2}{2{{\mathrm{\σ}}_p}^2}}\end{array}\right.---\left(1\right)$

其中I_p表示当前层影像，I_p-1表示上一层影像，σ_p表示当前层的尺度特征参数，G(x，y；σ_p)表示尺度参数为σ_p的高斯核函数，(x，y)表示高斯窗口内相对于中心的坐标偏移量；层间尺度参数的倍增关系和每一层的尺度参数分别满足下面的公式(2)和公式(3)：

k＝2^1/s    (2)

σ_p＝k^pσ    (3)

其中s表示总层数，p表示当前层数，k是层间尺度参数的倍增关系。

而且，步骤c中，利用高斯影像立方体的层间差异提取出的影像空间细节特征D，如以下公式(4)：

$D=\underset{i=1}{\overset{s}{\mathrm{\Σ}}}{D}_i=\underset{i=1}{\overset{s}{\mathrm{\Σ}}}(I_{i-1}-I_i)---\left(4\right)$

其中I_i-1和I_i分别表示高斯影像立方体中第i-1层和第i层的影像，I_i-1-I_i是这两层影像的层间差异，记为D_i。

而且，步骤c中，根据运算级别，选用以下三种提取空间细节特征的方法之一，实现提取亮度影像的空间细节特征D_Ave和原始高空间分辨率全色影像的空间细节特征D_pan，

第一种：逐层计算差异影像，提取空间细节特征，即按以下公式(5)提取空间细节特征：

$\left\{\begin{array}{c}{D}_{\mathrm{pan}}=\underset{i=1}{\overset{s}{\mathrm{\Σ}}}{D}_i=\underset{i=1}{\overset{s}{\mathrm{\Σ}}}(I_{i-1}-I_i)\\ D_{\mathrm{Ave}}=\underset{i=1}{\overset{s}{\mathrm{\Σ}}}{D_i}^{\′}=\underset{i=1}{\overset{s}{\mathrm{\Σ}}}({I_{i-1}}^{\′}-{I_i}^{\′})\end{array}\right.---\left(5\right)$

其中I_i-1和I_i表示原始高空间分辨率全色影像的高斯影像立方体中第i-1层和第i层的影像，I_i-1′和I_i′表示亮度影像的高斯影像立方体中第i-1层和第i层的影像；

第二种：在逐层计算差异影像的基础上，根据高斯函数的性质提取空间细节特征，即按以下公式(6)提取空间细节特征：

$\left\{\begin{array}{c}{D}_{\mathrm{pan}}=I-I_s=I-G(x,y;\sqrt{s}w,\sqrt{\underset{i=0}{\overset{s-1}{\mathrm{\Σ}}}{i}^2}\mathrm{\σ})*I\\ D_{\mathrm{Ave}}=I^{\′}-{I_s}^{\′}=I^{\′}-G(x,y;\sqrt{s}w,\sqrt{\underset{i=0}{\overset{s-1}{\mathrm{\Σ}}}{i}^2}\mathrm{\σ})*I^{\′}\end{array}\right.---\left(6\right)$

其中I和I′分别表示原始高空间分辨率全色影像和亮度影像，I_s、I_s′分别表示原始高空间分辨率全色影像和亮度影像的第s层影像，s表示总层数，w表示高斯核函数的窗口尺寸，σ表示初始设定的高斯核函数的方差，

表示窗口大小扩大了倍、尺度参数扩大了倍的高斯核函数；

第三种：在满足精度要求的前提下，对第二种方法简化，即按以下公式(7)提取空间细节特征：

$\left\{\begin{array}{c}{D}_{\mathrm{pan}}=I-I_s=I-G(x,y;w,\sqrt{\underset{i=0}{\overset{s-1}{\mathrm{\Σ}}}{i}^2}\mathrm{\σ})*I\\ D_{\mathrm{Ave}}=I^{\′}-{I_s}^{\′}=I^{\′}-G(x,y;w,\sqrt{\underset{i=0}{\overset{s-1}{\mathrm{\Σ}}}{i}^2}\mathrm{\σ})*I^{\′}\end{array}\right.---\left(7\right)$

其中I和I′分别表示原始全色影像和亮度影像，I_s、I_s′分别表示原始高空间分辨率全色影像和亮度影像的第s层影像，s表示总层数，w表示高斯核函数的窗口尺寸，σ表示初始设定的高斯核函数的方差，

表示窗口大小不变、尺度参数扩大了倍的高斯核函数。

而且，当多光谱影像包括三个单波段影像时，

步骤a中，根据下面的公式(8)来计算经过采样和校正后的多光谱影像的亮度影像Ave：

Ave＝(R+G+B)/3           (8)

其中，R、G、B分别是经过采样和校正后的多光谱影像的三个单波段影像；

步骤d中，所示采取加权融合的策略将步骤c中提取出的空间细节特征投影到原始低空间分辨率的多光谱影像上，得到的高空间分辨率融合多光谱影像按照以下公式(9)：

$\left\{\begin{array}{c}\mathrm{FR}=(2-c)*R+c*(D_{\mathrm{pan}}-D_{\mathrm{Ave}})\\ \mathrm{FG}=(2-c)*g+c*(D_{\mathrm{pan}}-D_{\mathrm{Ave}})\\ \mathrm{FB}=(2-c)*B+c(D_{\mathrm{pan}}-D_{\mathrm{Ave}})\end{array}\right.---\left(9\right)$

其中R、G、B分别是经过采样和校正后的多光谱影像的三个单波段影像，FR、FG、FB为融合影像的三个单波段影像，c为权重因子。

本发明提供的技术方案的有益效果为：(1)高斯尺度空间类似人类的视觉系统，从不同尺度的影像立方体上提取的影像细节特征能够更好地描述影像的空间特性。(2)采用空间投影的方法指导融合的过程，能够更好地保持影像的空间细节特征，同时良好地保持多光谱影像的光谱特征。

附图说明

图1为本发明实施例的流程图。

具体实施方式

为了更好地理解本发明的技术方案，下面结合附图对本发明做进一步的详细说明。本发明的实施例是对全色影像Pan和多光谱影像Mul进行融合，参照图1，本发明实施例的步骤如下：

步骤a：

首先进行影像预处理。实施例中影像的预处理包括重采样和精确配准(即本领域的精确几何校正处理)。重采样的目的是让原始多光谱影像具有与原始全色影像同等的像素大小，配准是为了让原始多光谱影像(即原始低空间分辨率的多光谱影像)与原始全色影像(即原始高空间分辨率全色影像)在空间位置上保持一致。根据融合精度要求，采样过程可以采取现有最邻近像元法、双线性内插法、双三次卷积法等方法，其中双三次卷积法精度最高，精确配准的精度在0.5个像素以内。

然后生成多光谱影像的亮度影像。先要根据重采样和配准后的多光谱影像生成其亮度影像，用以代表多光谱影像的空间特性，同时以便接下来生成高斯影像立方体。实施例的原始多光谱影像有三个波段，亮度影像的生成按下面的公式(8进行：

Ave＝(R+G+B)/3            (8)

其中，R、G、B分别是经过采样和校正后的多光谱影像的三个单波段影像。若具体实施时多光谱影像有大于三个的波段，则仍然可以在不影响理解的前提下类似命名并求取亮度影像。

步骤b：生成高斯影像立方体

本发明中引入改进的高斯尺度空间理论，分别生成原始全色影像和亮度影像的高斯影像立方体。高斯影像立方体中的各层影像应当满足以下公式(1)：

$\left\{\begin{array}{c}{I}_p=G(x,y;{\mathrm{\σ}}_p)*I_{p-1}\\ G(x,y;{\mathrm{\σ}}_p)=\frac{1}{2\mathrm{\π}{\mathrm{\σ}}^2}{e}^{-\frac{x^2+y^2}{2{{\mathrm{\σ}}_p}^2}}\end{array}\right.---\left(1\right)$

其中I_p表示当前层影像，I_p-1表示上一层影像，σ_p表示当前层的尺度特征参数，G(x，y；σ_p)表示尺度参数为σ_p的高斯核函数，(x，y)表示高斯窗口内相对于中心的坐标偏移量；层间尺度参数的倍增关系和每一层的尺度参数分别满足下面的公式(2)和公式(3)：

k＝2^1/s    (2)

σ_p＝k^pσ        (3)

其中s表示总层数(即高斯影像立方体包括从最底层的第1层影像到最顶层的第s层影像)，p表示当前层数，σ_p表示当前层的尺度参数。k就是层间尺度参数的倍增关系，即当前层影像(第p层)的尺度参数是上一层影像(第p-1层)的尺度参数的k倍。这个倍增关系跟经典高斯尺度空间理论中倍增关系是一样的。

因此具体实施时，可以设定初始参数包括：尺度参数σ，高斯核函数窗口大小w，高斯尺度空间层数s，也即是高斯影像立方体的层数，然后每一层的影像按以下公式生成：

I_p+1＝G(x，y；k^p+1σ)*I_p，p＝1，2，…，s，且k^s＝2

其中p表示当前层数，I_p+1和I_p分别表示高斯影像立方体的第p+1层和第p层影像，k^p+1σ表示第p+1层影像的尺度参数，G(x，y；k^p+1σ)表示第p+1层影像的高斯核函数。

步骤c：提取空间细节特征

本发明提出，可以根据不同的运算级别，采取下面三种等级的提取亮度影像和全色影像空间细节特征的方法，分别按照公式(5)(6)(7)进行。需要说明的是，这三种方法中只需要采用其中任意一种方法即可，而且采取何种方法视具体情况而定。

第一种：逐层计算差异影像，提取空间细节特征，可称为空间投影(SP，Spatial Projection)方法，满足下面的公式(5)：

$\left\{\begin{array}{c}{D}_{\mathrm{pan}}=\underset{i=1}{\overset{s}{\mathrm{\Σ}}}{D}_i=\underset{i=1}{\overset{s}{\mathrm{\Σ}}}(I_{i-1}-I_i)\\ D_{\mathrm{Ave}}=\underset{i=1}{\overset{s}{\mathrm{\Σ}}}{D_i}^{\′}=\underset{i=1}{\overset{s}{\mathrm{\Σ}}}({I_{i-1}}^{\′}-{I_i}^{\′})\end{array}\right.---\left(5\right)$

其中I_i-1和I_i表示原始全色影像的高斯影像立方体中第i-1层和第i层的影像，I_i-1′和I_i′表示亮度影像的高斯影像立方体中第i-1层和第i层的影像。

第二种：根据高斯函数的性质，在第一种的基础上推导出的方法，可称为快速空间投影(FSP，Fast Spatial Projection)方法，满足下面的公式(6)：

$\left\{\begin{array}{c}{D}_{\mathrm{pan}}=I-I_s=I-G(x,y;\sqrt{s}w,\sqrt{\underset{i=0}{\overset{s-1}{\mathrm{\Σ}}}{i}^2}\mathrm{\σ})*I\\ D_{\mathrm{Ave}}=I^{\′}-{I_s}^{\′}=I^{\′}-G(x,y;\sqrt{s}w,\sqrt{\underset{i=0}{\overset{s-1}{\mathrm{\Σ}}}{i}^2}\mathrm{\σ})*I^{\′}\end{array}\right.---\left(6\right)$

其中I和I′分别表示原始全色影像和亮度影像，I_s、I_s′分别表示原始高空间分辨率全色影像和亮度影像的第s层影像(即最顶层影像)，s表示总层数，w表示高斯核函数的窗口尺寸，σ表示初始设定的高斯核函数的方差，

表示窗口大小扩大了

倍、尺度参数扩大了

倍的高斯核函数。

第三种：在满足精度要求的前提下，对第二种方法进行适当的简化，可称为简化快速空间投影(SFSP，Simplified Fast Spatial Projection)方法，满足下面的公式(7)：

$\left\{\begin{array}{c}{D}_{\mathrm{pan}}=I-I_s=I-G(x,y;w,\sqrt{\underset{i=0}{\overset{s-1}{\mathrm{\Σ}}}{i}^2}\mathrm{\σ})*I\\ D_{\mathrm{Ave}}=I^{\′}-{I_s}^{\′}=I^{\′}-G(x,y;w,\sqrt{\underset{i=0}{\overset{s-1}{\mathrm{\Σ}}}{i}^2}\mathrm{\σ})*I^{\′}\end{array}\right.---\left(7\right)$

其中I和I′分别表示原始全色影像和亮度影像，I_s、I_s′分别表示原始高空间分辨率全色影像和亮度影像的第s层影像(即最顶层影像)，s表示总层数，w表示高斯核函数的窗口尺寸，σ表示初始设定的高斯核函数的方差，

表示窗口大小不变、尺度参数扩大了

倍的高斯核函数。

步骤d：空间投影运算

实施例的多光谱影像包括三个单波段影像，采取加权融合的策略将步骤c中提取出的空间细节特征投影到原始低空间分辨率的多光谱影像上时，按照以下公式得到高空间分辨率融合多光谱影像：

$\left\{\begin{array}{c}\mathrm{FR}=(2-c)*R+c*(D_{\mathrm{pan}}-D_{\mathrm{Ave}})\\ \mathrm{FG}=(2-c)*g+c*(D_{\mathrm{pan}}-D_{\mathrm{Ave}})\\ \mathrm{FB}=(2-c)*B+c(D_{\mathrm{pan}}-D_{\mathrm{Ave}})\end{array}\right.---\left(9\right)$

其中R、G、B为原始多光谱影像的三个波段，FR、FG、FB为融合影像的三个波段，c为权重因子。

具体实施时，可根据运算级别，即运算量由大到小，选择上述SP、FSP、SFSP方法实现空间细节特征的提取。还可参考运算所需的存储空间灵活选择，SP实现过程中要保存每一层的影像，需要的存储量是最大的，而FSP和SFSP实现过程中只需要保存最高一层影像，存储量要小很多。

在上述方法中，步骤a的具体操作属于遥感影像处理领域的常规方法，具体流程在此不再赘述。步骤a中，重采样和高精度配准操作的主要目的是使参与融合的全色影像和多光谱影像在空间地域上处于一致，具有同样的像素大小；获取亮度影像的主要目的是用以描述多光谱影像的空间特性。

在上述方法中，步骤b引入了改进的高斯尺度空间理论来生成高斯影像立方体。在经典高斯尺度空间理论中，尺度空间分为若干段，每段分为若干层，每一层影像都是由代表该层尺度特性的高斯核函数与原始影像卷积而成，而本方法中，只使用一段尺度空间，尺度参数由σ变化到2σ，该段中每一层影像都是由代表该层尺度特性的高斯核函数与下一层影像卷积而得到，一方面避免了经典高斯尺度空间理论中建立高层影像的庞大开销，另一方面也符合人眼视觉特性。可知生成的高斯影像立方体中，最底层影像为原始全色影像或亮度影像，最顶层影像满足下面的公式：

I_s＝G(x，y；2σ)*G(x，y；k^s-1σ)*…*G(x，y；σ)*I

在上述方法中，步骤c的三种提取空间特征的方法的推理和操作过程如下：

在高斯影像立方体中，第p层影像可以描述为：

I_p＝G(x，y；σ_p)*I_p-1

其中σ_p为高斯尺度空间中该层的尺度因子，则可以定义该层影像的空间细节特征D_i为：

D_p＝I_p-1-I_p

若初始设定的高斯尺度空间为s层，那么生成一个s层的原始全色影像的高斯影像立方，以及一个s层的原始全色影像的高斯影像立方体，则原始全色影像和亮度影像的空间细节特征定义为公式(5)：

$\left\{\begin{array}{c}{D}_{\mathrm{pan}}=\underset{i=1}{\overset{s}{\mathrm{\Σ}}}{D}_i=\underset{i=1}{\overset{s}{\mathrm{\Σ}}}(I_{i-1}-I_i)\\ D_{\mathrm{Ave}}=\underset{i=1}{\overset{s}{\mathrm{\Σ}}}{D_i}^{\′}=\underset{i=1}{\overset{s}{\mathrm{\Σ}}}({I_{i-1}}^{\′}-{I_i}^{\′})\end{array}\right.---\left(5\right)$

其中I_i-1和I_i表示原始全色影像的高斯影像立方体中第i-1层和第i层的影像，I_i-1′和I_i′表示亮度影像的高斯影像立方体中第i-1层和第i层的影像。此即为第一种方法：空间投影(SP，Spatial Projection)方法。

根据高斯函数的性质，若有两个高斯核函数与图像进行连续卷积，则卷积的结果图像等效于一个新的高斯核函数与图像进行卷积，这个新的高斯核函数的窗口大小的平方和是原来两个高斯核函数窗口大小的平方和，方差的平方和是原来两个高斯核函数方差的平方和。在生成高斯影像立方体的过程中，其高斯核的方差按照下面的公式从σ变化到2σ：

σ_p＝k^pσ

在整个卷积过程中，高斯核的窗口w大小保持不变。根据上述的高斯函数的性质，第p层影像的生成过程可以写成下面的公式：

$I_p=G(x,y;{\mathrm{\σ}}_p)*I_{p-1}=G(x,y;{\mathrm{\σ}}_p)*G(x,y;{\mathrm{\σ}}_{p-1})*I_{p-2}$

$=G(x,y;{\mathrm{\σ}}_p)*G(x,y;{\mathrm{\σ}}_{p-1})*...*G(x,y;\mathrm{\σ})*I$

$=G(x,y;\sqrt{p}w,\sqrt{\underset{i=0}{\overset{p-1}{\mathrm{\Σ}}}{i}^2}\mathrm{\σ})*I$

即第p次的卷积过程，可以等效于一个窗口大小为

方差为

的高斯核与图像进行卷积。

而在空间特征提取过程中，有下面的关系：

D＝(I-I₁)+(I₁-I₂)+……+(I_s-1-I_s)＝I-I_s

可知：

$I_s=G(x,y;\sqrt{s}w,\sqrt{\underset{i=0}{\overset{s-1}{\mathrm{\Σ}}}{i}^2}\mathrm{\σ})*I$

于是，就有：

$D=I-I_s=I-G(x,y;\sqrt{s}w,\sqrt{\underset{i=0}{\overset{s-1}{\mathrm{\Σ}}}{i}^2}\mathrm{\σ})*I$

则可以得到原始全色影像和亮度影像的空间细节特征如公式(6)：

$\left\{\begin{array}{c}{D}_{\mathrm{pan}}=I-I_s=I-G(x,y;\sqrt{s}w,\sqrt{\underset{i=0}{\overset{s-1}{\mathrm{\Σ}}}{i}^2}\mathrm{\σ})*I\\ D_{\mathrm{Ave}}=I^{\′}-{I_s}^{\′}=I^{\′}-G(x,y;\sqrt{s}w,\sqrt{\underset{i=0}{\overset{s-1}{\mathrm{\Σ}}}{i}^2}\mathrm{\σ})*I^{\′}\end{array}\right.---\left(6\right)$

其中I和I′分别表示原始全色影像和亮度影像，s表示总层数，w表示高斯核函数的窗口尺寸，σ表示初始设定的高斯核函数的方差，

表示窗口大小扩大了倍、尺度参数扩大了

倍的高斯核函数。此即为第二种方法：快速空间投影(FSP，FastSpatial Projection)方法。

从式(6)可以看出，快速算法中，用一个窗口为方差为

的新的高斯核卷积去等效高斯影像立方体的特征提取过程，由于卷积窗口增大而导致效率提高不明显，作为简化算法，在不影响融合效果的前提下，可以将新的高斯核的窗口仍然设置为w，方差设置为

由此，新的特征提取算法可以写为：

$\left\{\begin{array}{c}{h}_{\mathrm{pan}}=I-I_s=I-G(x,y;w,\sqrt{\underset{i=0}{\overset{s-1}{\mathrm{\Σ}}}{i}^2}\mathrm{\σ})*I\\ h_{\mathrm{Ave}}=I^{\′}-{I_s}^{\′}=I^{\′}-G(x,y;w,\sqrt{\underset{i=0}{\overset{s-1}{\mathrm{\Σ}}}{i}^2}\mathrm{\σ})*I^{\′}\end{array}\right.---\left(7\right)$

其中I和I′分别表示原始全色影像和亮度影像，s表示总层数，w表示高斯核函数的窗口尺寸，σ表示初始设定的高斯核函数的方差，

表示窗口大小不变、尺度参数扩大了

倍的高斯核函数。此即为第三种方法：简化快速空间投影(SFSP，SimplifiedFast Spatial Projection)方法。

在上述方法中，步骤d的目的是将提取的空间细节特征添加到原始多光谱影像中，以提高多光谱影像的空间细节信息。加权融合时的权重因子c对于融合结果有着重要影响，默认情况下c＝1，若c＞1，则在融合影像中全色影像的空间细节所占比重比较大，否则多光谱影像的光谱特征所占比重比较大。

综上所述，本发明提出的多尺度空间投影融合方法，在改进的高斯尺度空间理论的基础上，利用高斯影像立方体的层间差异特征提取影像的空间细节特征。高斯尺度空间能够模仿人类视觉感知功能，用一系列尺度参数来描述影像的特征，在生成高斯影像立方体的过程中，使用高斯核函数进行连续卷积，逐渐模糊的影像之间的差异即为空间细节特征，因此从高斯影像立方体中提取的空间细节特征更符合人眼视觉感受。

以下通过仿真实验来验证本发明的有效性：

仿真实验一：QuickBird全色(空间分辨率为0.6米)与多光谱(空间分辨率为2.4米)影像融合，原始影像大小为600*600，选择广为使用的标准IHS变换方法作为对。

融合影像评价指标：对于融合影像的质量评价应该从主观视觉和客观指标两方面进行，同时考察融合结果影像对空间信息和光谱信息的保持度。本发明中的主观视觉评价为目视感受，客观评价指标采取以下三种：

(1)ERGAS：一个从全局综合误差方面来评价融合影像质量的指标，ERGAS值越小，说明融合影像质量越好。由下面的公式来计算：

$\mathrm{ERGAS}=100*\frac{h}{l}*\sqrt{\frac{1}{N}\underset{i=1}{\overset{N}{\mathrm{\Σ}}}\frac{{\mathrm{RMSE}}^2\left(B_i\right)}{\mathrm{Mean}{\left(B_i\right)}^2}}$

其中，h和l分别表示全色和多光谱影像的空间分辨率，N为原始多光谱影像的波段数，RMSE(B_i)和Mean(B_i)分别表示第i波段的均方根误差和均值。

(2)通用图像质量评价指标UIQI：从相关信息损失、辐射值扭曲和对比度扭曲三个方面衡量融合前后影像的相似度，通过比较UIQI的值，其值越大，表示融合质量越高，若两幅影像完全一致，其值为1。由下面的公式来计算：

$\mathrm{UIQI}=\frac{4{\mathrm{\δ}}_{\mathrm{xy}}\stackrel{\‾}{x}\stackrel{\‾}{y}}{({\mathrm{\δ}}_x^2+{\mathrm{\δ}}_y^2)+[{\left(\stackrel{\‾}{x}\right)}^2+{\left(\stackrel{\‾}{y}\right)}^2]}$

其中

和

分别表示原始多光谱影像和融合结果影像的均值，δ_x和δ_y分别表示原始多光谱影像和融合结果影像的方差，而δ_xy表示它们的协方差。

(3)光谱偏差指数：反映了融合影像与原始多光谱影像在光谱方面的扭曲度，光谱偏差指数的值越小，说明融合质量越高。由下面的公式来计算：

$d=\frac{1}{N}\underset{i}{\mathrm{\Σ}}\underset{j}{\mathrm{\Σ}}\frac{|I_{i,j}-{\stackrel{\$}{I}}_{i,j}|}{I_{i,j}}$

其中I_i，j和

分别表示融合前后多光谱影像的灰度值。

仿真结果：

用本发明的方法和标准IHS变换方法对仿真内容的仿真结果影像进行对比，包括原始全色影像、经过采样和配准的原始多光谱影像、标准IHS变换的融合结果影像，及本发明所提供SP、FSP、SFSP方法得到的影像。

(1)按照所述仿真内容的仿真结果客观评价指标如表1所示：

表1实验结果比较

从表1的实验数据可知，本发明提出的多尺度空间投影融合方法在各指标上均比传统的IHS变换方法有非常明显的优势。ERGAS的值越大，说明光谱信息保持度越低，标准IHS变换在这一指标上的值远大于本发明的三种方法。UIQI指标反映了融合结果影像在相关信息损失、辐射值扭曲和对比度扭曲，IHS变换在这一指标上的值远小于本发明的三种方法。偏差指数直接反应了光谱扭曲程度，值越大，说明光谱扭曲越严重，而标准IHS变换的值远大于本发明的三种方法。由此可以得出结论，相比于标准IHS变换，本发明提出的方法更能保持影像的空间细节特征和光谱特征。

(2)按照所述仿真内容的仿真结果主观评价如下：

从视觉上看，本发明所得到的融合影像在亮度上与原始多光谱影像非常接近，明暗适度，色彩鲜艳，边缘清晰，对比度好，层次感强，具有很好的视觉感受。在图像细节方面，本发明所得的融合影像与原始全色影像非常相似，比如房屋边缘、树木的结构、道路的形态及道路上的车辆，都能很清晰地呈现。虽然IHS变换在细节呈现方面也有较好的表现，但是存在较严重的色彩偏差和偏色，不太符合人眼的视觉感受。

综上所述，本发明具有如下的优点：

(1)高斯尺度空间类似人类的视觉系统，从不同尺度的影像立方体上提取的影像细节特征能够更好地描述影像的空间特性。

(2)采用空间投影的方法指导融合的过程，能够更好地保持影像的空间细节特征，同时良好地保持多光谱影像的光谱特征。

与传统影像融合方法比较可知，不管是从客观评价指标上，还是从主观视觉上，本发明的方法都具有很明显的优势，融合结果既能较好地保持高分辨率全色影像的空间细节特征，又能较好地保持低分辨率多光谱影像的光谱特征，增强了主观可视表达效果，是一种可行的遥感影像融合方法。

Claims (5)

1.一种多尺度空间投影遥感影像融合方法，其特征在于包括以下步骤：

步骤a，将相同地域的原始低空间分辨率的多光谱影像采样到原始高空间分辨率全色影像同等像素大小，进行精确几何校正处理，然后利用经过采样和校正后的多光谱影像计算其亮度影像；

步骤b，分别生成原始高空间分辨率全色影像和步骤a所得亮度影像的高斯影像立方体；所述高斯影像立方体由若干层的影像组成，当前层影像是由上一层影像与体现当前层尺度特征的高斯核函数卷积而成，设最底层的尺度参数为σ，最底层影像即为原始高空间分辨率全色影像和步骤a所得亮度影像与体现尺度参数σ的高斯核函数卷积而成，从最底层影像到最顶层影像的尺度参数由σ变化到2σ，并且当前层影像的尺度参数与上一层影像的尺度参数满足倍增关系；高斯影像立方体中的各层影像满足以下公式(1)：

$\left\{\begin{array}{c}{I}_p=G(x,y;{\mathrm{\σ}}_p)*I_{p-1}\\ G(x,y;{\mathrm{\σ}}_p)=\frac{1}{2\mathrm{\π}{\mathrm{\σ}}^2}{e}^{-\frac{x^2+y^2}{2{{\mathrm{\σ}}_p}^2}}\end{array}\right.---\left(1\right)$

其中I_p表示当前层影像，I_p-1表示上一层影像，p的取值从2到s，σ_p表示当前层的尺度特征参数，G(x，y；σ_p)表示尺度参数为σ_p的高斯核函数，(x，y)表示高斯窗口内相对于中心的坐标偏移量；层间尺度参数的倍增关系和每一层的尺度参数分别满足下面的公式(2)和公式(3)：

k＝2^1/(s-1)    (2)

σ_p＝k^p-1σ    (3)

其中s表示总层数，p表示当前层数，p的取值从1到s，k是层间尺度参数的倍增关系；

步骤c，利用步骤b所得高斯影像立方体的层间差异，提取亮度影像和原始高空间分辨率全色影像的空间细节特征；利用高斯影像立方体的层间差异提取出的影像空间细节特征D，如以下公式(4)：

$D=\underset{i=1}{\overset{s-1}{\mathrm{\Σ}}}{D}_i=\underset{i=2}{\overset{s}{\mathrm{\Σ}}}(I_{i-1}-I_i)---\left(4\right)$

其中I_i-1和I_i分别表示高斯影像立方体中第i-1层和第i层的影像，I_i-1-I_i是这两层影像的层间差异，记为D_i；

步骤d，将步骤c提取出的空间细节特征按照加权融合策略投影到原始低空间分辨率的多光谱影像上，得到具有高空间分辨率和高光谱分辨率的融合多光谱影像。

2.根据权利要求1所述的多尺度空间投影遥感影像融合方法，其特征在于：步骤a中，采样过程采取最邻近像元法、双线性内插法或双三次卷积法方法。

3.根据权利要求1所述的多尺度空间投影遥感影像融合方法，其特征在于：步骤a中，利用经过采样和校正后的多光谱影像计算其亮度影像的具体方式为，对多光谱影像的所有单波段影像求平均亮度影像。

4.根据权利要求1所述的多尺度空间投影遥感影像融合方法，其特征在于：步骤c中，根据运算级别，选用以下三种提取空间细节特征的方法之一，实现提取亮度影像的空间细节特征D_Ave和原始高空间分辨率全色影像的空间细节特征D_pan，

第一种：逐层计算差异影像，提取空间细节特征，即按以下公式(5)提取空间细节特征：

$\left\{\begin{array}{c}{D}_{\mathrm{pan}}=\underset{i=1}{\overset{s-1}{\mathrm{\Σ}}}{D}_i=\underset{i=2}{\overset{s}{\mathrm{\Σ}}}(I_{i-1}-I_i)\\ D_{\mathrm{Ave}}=\underset{i=1}{\overset{s-1}{\mathrm{\Σ}}}{D_i}^{\′}=\underset{i=2}{\overset{s}{\mathrm{\Σ}}}({I_{i-1}}^{\′}-{I_i}^{\′})\end{array}\right.---\left(5\right)$

其中I_i-1和I_i表示原始高空间分辨率全色影像的高斯影像立方体中第i-1层和第i层的影像，I_i-1′和I_i′表示亮度影像的高斯影像立方体中第i-1层和第i层的影像；

第二种：在逐层计算差异影像的基础上，根据高斯函数的性质提取空间细节特征，即按以下公式(6)提取空间细节特征：

$\left\{\begin{array}{c}{D}_{\mathrm{pan}}=I-I_s=I-G(x,y;\sqrt{s}w,\sqrt{\underset{i=0}{\overset{s-1}{\mathrm{\Σ}}}{i}^2\mathrm{\σ}})*I\\ D_{\mathrm{Ave}}=I^{\′}-{I_s}^{\′}=I^{\′}-G(x,y;\sqrt{s}w,\sqrt{\underset{i=0}{\overset{s-1}{\mathrm{\Σ}}}{i}^2\mathrm{\σ}})*I^{\′}\end{array}\right.---\left(6\right)$

其中I和I′分别表示原始高空间分辨率全色影像和亮度影像，I_s、I_s′分别表示原始高空间分辨率全色影像和亮度影像的第s层影像，s表示总层数，w表示高斯核函数的窗口尺寸，σ表示初始设定的高斯核函数的方差，

表示窗口大小扩大了

倍、尺度参数扩大了倍的高斯核函数；

第三种：在满足精度要求的前提下，对第二种方法简化，即按以下公式(7)提取空间细节特征：

$\left\{\begin{array}{c}{D}_{\mathrm{pan}}=I-I_s=I-G(x,y;w,\sqrt{\underset{i=0}{\overset{s-1}{\mathrm{\Σ}}}{i}^2}\mathrm{\σ})*I\\ D_{\mathrm{Ave}}=I^{\′}-{I_s}^{\′}=I^{\′}-G(x,y;w,\sqrt{\underset{i=0}{\overset{s-1}{\mathrm{\Σ}}}{i}^2}\mathrm{\σ})*I^{\′}\end{array}\right.---\left(7\right)$

其中I和I′分别表示原始全色影像和亮度影像，I_s、I_s′分别表示原始高空间分辨率全色影像和亮度影像的第s层影像，s表示总层数，w表示高斯核函数的窗口尺寸，σ表示初始设定的高斯核函数的方差，

表示窗口大小不变、尺度参数扩大了

倍的高斯核函数。

5.根据权利要求4所述的多尺度空间投影遥感影像融合方法，其特征在于：当多光谱影像包括三个单波段影像时，

步骤a中，根据下面的公式(8)来计算经过采样和校正后的多光谱影像的亮度影像Ave：

Ave＝(R+G+B)/3    (8)

其中，R、G、B分别是经过采样和校正后的多光谱影像的三个单波段影像；

步骤d中，所示采取加权融合的策略将步骤c中提取出的空间细节特征投影到原始低空间分辨率的多光谱影像上，得到的高空间分辨率融合多光谱影像按照以下公式(9)：

$\left\{\begin{array}{c}\mathrm{FR}=(2-c)*R+c*(D_{\mathrm{pan}}-D_{\mathrm{Ave}})\\ \mathrm{FG}=(2-c)*G+c*(D_{\mathrm{pan}}-D_{\mathrm{Ave}})\\ \mathrm{FB}=(2-c)*B+c^{*}(D_{\mathrm{pan}}-D_{\mathrm{Ave}})\end{array}\right.---\left(9\right)$

其中R、G、B分别是经过采样和校正后的多光谱影像的三个单波段影像，FR、FG、FB为融合影像的三个单波段影像，c为权重因子。

Images