CN1472544A - 遥感影像像素与特征联合最优融合方法 - Google Patents

Abstract

一种遥感影像像素与特征联合最优融合方法，在IHS空间，将多光谱影像的强度分量经小波多层分解得到的低频基带系数，与高分辨率影像经对应多层小波分解得到低频基带系数进行以空间细节信息和光谱信息两项特征指标的最优像素级融合，对经小波分解后的高频系数进行特征级融合，然后将小波系数进行对应的小波逆变换，得到新的强度分量，再进行IHS逆变换后得到融合后的影像。本发明结合了IHS融合方法和小波融合方法的特点，通过分别对小波基带系数的权系数进行像素级求优融合和高频子带系数的特征级融合，使得融合后的影像既达到最高的空间分辨率，同时又最大限度的降低了彩色的畸变，有效的改善了融合影像的光谱信息指标。

Description

遥感影像像素与特征联合最优融合方法

技术领域：

本发明涉及一种遥感影像像素与特征联合最优融合方法，结合小波多分辨率分解具有的时频特性和IHS(强度Intensity-色度Hue-饱和度Saturation)变换融合方法，进行像素级和特征级联合最优融合，是遥感影像融合的核心技术，在各类军用或民用的遥感信息处理***、数字城市空间信息***等领域中均可有广泛的应用。

背景技术：

如何有效的融合高分辨率全色遥感影像和低分辨率的多光谱遥感影像，均衡融合结果中空间细节信息和光谱信息两项特征指标，是多源遥感影像融合技术的研究热点之一。

Haydn等人首先提出的IHS融合方法是经典的实用算法之一。该方法将多光谱影像通过IHS变换从RGB空间变换到IHS空间，同时将高分辨率的全色影像进行线性拉伸，使得拉伸后影像的均值和方差与IHS空间中的强度分量I₀一致。然后，将拉伸后的高分辨率影像作为新的强度分量，与H和S分量一起按照IHS逆变换公式变换到原RGB空间。这样，使得融合后的影像既具有较高的空间分辨率，同时又保持了原低分辨率多光谱影像相同的色度和饱和度。然而，这种经典的IHS融合方法存在着一定的缺陷，由于不同波段的数据具有不同的光谱特性曲线，IHS融合方法扭曲了原始的光谱特性，产生了不同程度的光谱退化现象，因而不利于影像的正确识别和分类，特别是对于不同时相的多传感器遥感影像的影像融合，IHS融合方法无法使得融合影像的色调和原多光谱影像的色调保持一致，这种因为光谱信息的变换，导致了影像不能用于地物识别和反演。Te-Ming等人在IHS空间进行了数学上的证明，论述了IHS融合方法的缺陷，得到的结论是：尽管用于替换强度分量I₀的高分辨率的全色影像I_new在替换前进行了影像的统计特性的匹配，但是匹配误差δ＝I_new-I导致了彩色的畸变。

发明内容：

本发明的目的在于针对上述IHS变换融合技术的不足，提供一种遥感影像像素与特征联合最优融合方法，既能提高融合后影像的空间分辨率，又能降低彩色的畸变，有效改善融合影像的光谱信息指标。

为实现这样的目的，本发明在IHS空间，将多光谱影像的强度分量经小波多层分解得到的低频基带系数，与高分辨率影像经对应多层小波分解得到低频基带系数进行以空间细节信息和光谱信息两项特征指标的最优像素级融合，对经小波分解后的高频系数进行特征级融合，然后将小波系数进行对应的小波逆变换，得到新的强度分量，再进行IHS逆变换后得到融合后的影像。

由于小波变换在变换域具有良好的分频特性，小波系数的统计特性反映了遥感影像的边缘、线和区域等显著特征，本发明将小波变换的多分辨分析(Multi-resolution Analysis)方法引入高分辨率全色遥感影像和低分辨率的多光谱遥感影像的融合中。

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

(1)对多光谱影像B进行IHS变换，分别得到在IHS彩色空间的色度H、饱和度S和强度分量I，然后对I分量进行3层小波分解；

(2)对高分辨率影像A进行线性拉伸和直方图匹配，然后进行小波分解，分解层数也为3层；

(3)确定一3×3的空域窗口，分别对影像A的I分量和影像B的子带小波系数进行统计，求得均值μ(2^j)和方差D(2^j)；

(4)对应分辨率层的高频子带系数按照(1)式进行特征级融合：

(1)式中，2^j为小波分解层数，W^k(2^j，x，y)为2^j分辨率下得到的高频子带系数融合结果；W_A ^k(2^j，x，y)和W_B ^k(2^j，x，y)为影像A和I’的对应高频子带系数，D_A ^k、D_B ^k分别是以(x，y)为中心像元的3×3的空域窗口的方差；

(5)对影像A的I分量和影像B的子带小波系数的低频基带系数，按照(2)式进行最优像素级融合，权系数的k_opt按照光谱信息评价指标和空间分辨率评价指标进行求优：

A(2^j，x，y)＝k₁A_A(2^j，x，y)+k₂A_B(2^j，x，y)              (2)

(2)式中，A_A(2^j，x，y)、A_B(2^j，x，y)分别为影像A和I’的对应2^j分辨率的低频基带数据，k₁、k₂为需要求优的权系数。按照归一化要求，k₁、k₂满足k₁+k₂＝1，即求优权系数的确定可以归结为满足目标函数的k_opt＝k₁＝1-k₂；

(6)对得到的像素级最优融合的小波系数的低频基带系数，以及进行的特征级融合的各高频小波系数进行对应的小波逆变换，得到新的强度分量I’；

(7)将I’、H、S进行IHS逆变换，得到融合后的影像C。

本发明结合了IHS融合方法和小波融合方法的特点，通过分别对小波基带系数的权系数进行了像素级求优融合和高频子带的特征级融合，其有益效果体现为：使得融合后的影像既达到最高的空间分辨率，同时又最大限度的降低了彩色的畸变。均衡了融合结果中空间细节信息和光谱信息两项特征指标，有效的改善了融合影像的光谱信息指标。

附图说明：

图1为本发明—基于小波统计特性的遥感图像像素与特征联合最优特征融合方法流程图。

图2为本发明遥感影像融合结果和IHS方法的对比。其中，图2(a)为多光谱遥感影像(512×512)；图2(b)为高分辨率全色遥感影像(512×512)；图2(c)为IHS方法的融合结果(512×512)，图2(d)为本发明的融合结果(512×512)。

图3为小波分解后的各子带系数的特征级融合。其中，图3(a)～(c)分别为全色高分辨率影像(P-Panchromatic)和多光谱影像(M-Multispectral)以及特征级融合结果影像(F-Fused)的水平高频图像；图3(d)～(f)分别为全色、多光谱垂直高频图像及其特征级融合结果影像；图3(g)～(i)分别为全色、多光谱斜线高频图像及其特征级融合结果影像。

图4为本发明的性能评价指标曲线。

具体实施方式：

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

本发明的详细流程如图1。本发明选取一多光谱遥感影像A(512×512)如图2(a)，高分辨率全色遥感影像B(512×512)如图2(b)，A、B严格配准后，实施如下步骤：

1、对多光谱影像B进行IHS变换，分别得到在IHS彩色空间的色度H、饱和度S和强度分量I，然后对I分量进行3层小波分解；

2、对高分辨率影像A进行线性拉伸和直方图匹配，然后进行小波分解，分解层数也为3层；

3、确定一3×3的空域窗口，分别对影像A的I和影像B的子带小波系数进行统计，求得均值μ(2^j)和方差D(2^j)；

4、按照(1)式确定对应分辨率层的高频子带系数，进行特征级融合；可以得到对应分辨率层的特征级融合结果，如图3所示；图3为小波分解后的各子带系数的特征级融合。图3(a)～(c)分别为全色高分辨率影像(P-Panchromatic)和多光谱影像(M-Multispectral)以及特征级融合结果影像(F-Fused)的水平高频图像；图3(d)～(f)分别为全色、多光谱垂直高频图像及其特征级融合结果影像；图3(g)～(i)分别为全色、多光谱斜线高频图像及其特征级融合结果影像；

5、对影像A的I和影像B的子带小波系数的低频基带系数，按照(2)式进行像素级融合，权系数k_opt按照(3)式所示的目标函数求优：

F(k_opt)＝F(k₁，k₂)

＝Max(E_SP(k₁，k₂)； E_HF(k₁，k₂))OR Max(E_SP(k₁，k₂)；AG(k₁，k₂))  (3)

＝{E_SP(k_opt)；E_HF(k_opt)}OR{E_SP(k_opt)；AG(k_opt)}

   D:g(k_opt):k_opt＝k₁＝1-k₂(0≤k_opt≤1)

             k_opt∈D∈R

其中，融合图像寻优评价指标分别定义如下：

(i)光谱信息评价指标

利用融合图像与多光谱图像的相关程度来定义光谱信息的评价指标。

令f为融合后图像，f₀为多光谱图像。定义光谱信息评价E_SP指标如(4)式。 $E_{\mathrm{SP}}=\mathrm{Corr}(f,f_0)=\frac{{\mathrm{\Σ}}_{j=1}^{\mathrm{npix}}(f_i-\stackrel{\‾}{f})(f_{0j}-{\stackrel{\‾}{f}}_0)}{\sqrt{{\mathrm{\Σ}}_{j=1}^{\mathrm{npix}}(f_i-\stackrel{\‾}{f}{)}^2{(f_{0j}-{\stackrel{\‾}{f}}_0)}^2}}........\left(4\right)$

其中，npix是图像中像素点的个数， f和 f₀表示图像的灰度均值，相关程度Corr(f，f₀)反映了影像f和f₀的相似程度。

(ii)空间分辨率评价指标

利用融合图像对应的灰度的高频分量与高分辨率图像高频分量之间的相关程度来定义空间分辨率指标。令f_H为高分辨率全色影像。首先将融合的图像转化为灰度图像，然后进行小波分解，得到融合图像的四个分量(f^a，f^h，f^v，f^d)，分别表示融合图像的低频分量、水平方向的高频分量、垂直方向的高频分量和对角线方向的高频分量。同样也可以得到高分辨率图像小波分解的四个分量(f_H ^a，f_H ^h，f_H ^v，f_H ^d)。定义空间分辨率评价指标如(5)式所示： $E_{\mathrm{HF}}=\frac{\mathrm{Corr}(f^h,f_H^h)+\mathrm{Corr}(f^v,f_H^v)+\mathrm{Corr}(f^d,f_H^d)}{3}......\left(5\right)$

另外，图像清晰度质量的改进最直接的指标是图像的平均梯度，它反映了影像的清晰程度，同时还反映出图像中微小细节反差和纹理变换特征，定义平均梯度(Average grads)评价指标如(6)式所示： $\mathrm{AG}=\frac{1}{M*N}\underset{i=1}{\overset{m}{\mathrm{\Σ}}}\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{[\mathrm{\Δxf}(i,j)}^2+\mathrm{\Δyf}{(i,j)}^2{]}^{1/2}......\left(6\right)$

(6)式中，Δxf(i，j)、Δyf(i，j)分别为像元(i，j)在x/y方向上的一阶差分；

依据融合图像的寻优评价指标修正融合准则中的基带数据融合的权系数k₁、k₂。随着高分辨率影像低频基带系数的融合权值k₁的增加，空间分辨率评价指标E_HF和AG随之增大，光谱信息评价指标E_SP减小。因此，(3)式中的k_opt是使得目标函数达到最大的权系数，即使得E_SP、E_HF或者E_SP、AG达到同时最大的权系数。k_opt取值范围为[0，1]。

根据权系数k_opt的求优目标函数，在取值区间[0，1]内，随着权系数融合权值k₁的增加，按(7)式得到空间分辨率评价指标E_HF和AG和光谱信息评价指标E_SP的曲线，曲线交点即为最优权系数的k_opt，如图4所示。 $E\left(i\right)=\frac{E\left(i\right)-\mathrm{MinE}\left(i\right)}{\mathrm{Max}[E\left(i\right)-\mathrm{Min}(E\left(i\right)]}......\left(7\right)$

(7)式中，E(i)为评价指标，i为权系数k₁寻优的次数。图4所示为在[0，1]以寻优步长0.001得到的E_SP、E_HF、AG曲线。可以看到，归一化后的E_SP、E_HF、AG都呈非线性，E_HF随融合权值k₁的增加而增大的趋势大于AG的增大的趋势。满足求优目标函数的权系数k_opt是E_SP和E_HF、E_SP和AG的交点k_opt1＝0.474或k_opt2＝0.586(收敛精度δ≤0.001)。因此，k_opt1和k_opt2均衡了融合结果中空间细节信息和光谱信息两项特征指标，达到最优。换言之，在按照(4)～(6)式，进行客观定性评价计算可以得到最优权系数k_opt1或k_opt2，其使得融合后的影像既达到最高的空间分辨率，同时又最大限度的降低了彩色的畸变；

6、对得到的像素级最优融合的小波系数的低频基带系数，以及进行的特征级融合的各高频小波系数进行对应的小波逆变换，得到新的强度分量I’；

7、将I’、H、S进行IHS逆变换，得到融合后的影像C；

融合结果比较：图2(c)为IHS方法的融合结果(512×512)，图2(d)为本发明的融合结果(512×512)。

性能分析：表1为选取最优权系数时，融合结果的相关系数和平均梯度两项指标定量值，并与IHS方法进行了对比。

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

	波段	相关系数		平均梯度
						IHS	R	0.5624		10.6146
G	0.4644		11.3663
					B		0.5060		9.5875
本发明方法	ki	kopt1＝k1＝0.474；k2＝0.526		kopt1＝k1＝0.586；k2＝0.414
						相关系数	平均梯度	相关系数	平均梯度
		R	0.8488	8.9448	0.7950					9.0954
	G					0.7691	9.4102	0.6992	9.6055
		B	0.8247	7.9353	0.7626					8.1007

Claims (3)

1、一种遥感影像像素与特征联合最优融合方法，其特征在于包括如下具体步骤：

(1)对多光谱影像B进行IHS变换，分别得到在IHS彩色空间的色度H、饱和度S和强度分量I，然后对I分量进行3层小波分解；

(2)对高分辨率影像A进行线性拉伸和直方图匹配，然后进行小波分解，分解层数也为3层；

(3)确定一3×3的空域窗口，分别对影像A的I分量和影像B的子带小波系数进行统计，求得均值μ(2^j)和方差D(2^j)；

(4)对应分辨率层的高频子带系数按照公式进行特征级融合，式中，2^j为小波分解层数，W^k(2^j，x，y)为2^j分辨率下得到的高频子带系数融合结果，W_A ^k(2^j，x，y)和W_B ^k(2^j，x，y)为影像A和I’的对应高频子带系数，D_A ^k、D_B ^k分别是以(x，y)为中心像元的3×3的空域窗口的方差；

(5)对影像A的I分量和影像B的子带小波系数的低频基带系数，按照公式

              A(2^j，x，y)＝k₁A_A(2^j，x，y)+k₂A_B(2^j，x，y)进行最优像素级融合，权系数的k_opt按照光谱信息评价指标和空间分辨率评价指标进行求优，式中，A_A(2^j，x，y)、A_B(2^j，x，y)分别为影像A和I’的对应2^j分辨率的低频基带数据，k₁、k₂为需要求优的权系数，k₁、k₂满足k₁+k₂＝1，即求优权系数的确定归结为满足目标函数的k_opt＝k₁＝1-k₂；

(6)对得到的像素级最优融合的小波系数的低频基带系数，以及进行的特征级融合的各高频小波系数进行对应的小波逆变换，得到新的强度分量I’；

(7)将I’、H、S进行IHS逆变换，得到融合后的影像C。

2、如权利要求1所说遥感影像像素与特征联合最优融合方法，其特征在于权系数k_opt进行求优时，求优目标函数按照式

 F(k_opt)＝F(k₁，k₂)

 ＝Max(E_SP(k₁，k₂)；E_HF(k₁，k₂))OR Max(E_SP(k₁，k₂)；AG(k₁，k₂))

 ＝{E_SP(k_opt)；E_HF(k_opt)}OR{E_SP(k_opt)；AG(k_opt)}

   D:g(k_opt):k_opt＝k₁＝1-k₂(0≤k_opt≤1)

             k_opt∈D∈R

确定。

3、如权利要求1所说的遥感影像像素与特征联合最优融合方法，其特征在于权系数k_opt按照光谱信息评价和空间分辨率评价指标进行求优时，利用式 $E\left(i\right)=\frac{E\left(i\right)-\mathrm{MinE}\left(i\right)}{\mathrm{Max}[E\left(i\right)-\mathrm{Min}(E\left(i\right)]}$

得到光谱信息评价和空间分辨率评价指标的曲线，曲线交点即为最优权系数的k_opt，式中，E(i)为评价指标，i为权系数k₁寻优的次数。