专利名称：基于贝叶斯线性估计的遥感图像融合方法
技术领域：
本发明属于遥感图像处理技术领域，具体涉及一种基于贝叶斯线性估计的遥感图像融合方法。
背景技术：
由于遥感传感器设计的限制，遥感图像一般在空间分辨率和光谱分辨率之间进行折衷，具有高的光谱分辨率的图像一般不具备高的空间分辨率，反之亦然。例如，Landsat ETM+传感器提供的就是6幅30m空间分辨率的多光谱波段图像和一幅15m空间分辨率的全色波段图像。在实际应用中，那些同时具有高的空间分辨率和高的光谱分辨率的图像能够有效地提高解译和分类的精度，因此不同分辨率的遥感图像的融合成为研究的热点，尤其是低分辨率的多光谱图像和高分辨率的全色波段图像的融合。一般来说，融合后的图像既要求保留多光谱图像的光谱特性，又要求融入全色图像的空间信息。
目前常见的融合方法有HIS方法[1]、PCA方法[2，3]以及小波变换方法[4，5]。HIS方法和PCA方法常常显著地改变原多光谱图像的光谱特性，而小波变换方法对于分解层次和小波基的选择比较敏感，并且会因操作人员的不同，而有不同的融合效果。
近年来，基于统计参数估计的方法[6-9]开始应用于遥感图像的融合中。Nishii等[8]假设高分辨率多光谱图像和全色图像的概率分布服从联合高斯分布，用条件均值作为估计。Hardie等[9]在上述联合高斯分布假设的基础上引入了高分辨率多光谱图像和低分辨率多光谱图像之间的观测模型，得到了最大后验概率(MAP)意义上的估计。Nishii方法和Hardie方法在低分辨率多光谱图像和全色图像相关性不强的时候，难以融入全色图像的空间信息。
针对以上问题，在遥感图像融合的研究中，在增强空间细节的同时很好地保持光谱特性，以及保证算法的鲁棒性成为目前研究的热点。
发明内容
本发明的目的是提出一种基于贝叶斯线性估计的遥感图像融合方法，以解决传统的基于统计参数估计方法依赖于多光谱图像和全色图像相关系数的问题，增强空间细节，保持光谱特性。
本发明提出的基于贝叶斯线性估计的遥感图像融合方法，具体步骤如下引入高分辨率多光谱图像和低分辨率多光谱图像之间的观测模型，以及高分辨率多光谱图像和全色图像之间的观测模型，并将上述两个观测模型联立成一个贝叶斯线性模型；应用贝叶斯高斯-马尔科夫定理，计算得到线性最小均方误差(LMMSE)意义上的高分辨率多光谱图像的估计。
下面对各步骤作进一步具体描述。
1、引入观测模型假设同一地区分别被低分辨率的多光谱波段传感器和高分辨率的全色波段传感器拍摄到，所谓图像分辨率是指图像中每个像素覆盖地面的范围，图像分辨率的高低是个相对的概念，根据需要可具体划分。全色图像按照如下方式排列成一维的列矢量x＝[x1，x2，…，x1，…，xN]T(1)其中xi表示全色图像在空间位置i上的像素值，而N是全色图像的像素数目。低分辨率的多光谱图像也按照类似的方式排列成一维的列矢量y=[y1T,y2T,···,yjT,···,yMT]T---(2)]]>其中yj表示低分辨率的多光谱图像在空间位置j上的像素值(有K个波段，yj＝[yj，1，yj，2，…，yj，k]T)，而M是低分辨率多光谱图像的像素数目。
假设高分辨率的多光谱图像存在，那么它应该既包含多光谱图像的光谱信息，又具有和全色图像一样的空间分辨率，这里用如下的一维列矢量来表示z=[z1T,z2T,···,ziT,···,zNT]T---(3)]]>其中zi表示高分辨率的多光谱图像在空间位置i上的像素值(有K个波段，zi＝[zi，1，zi，2，…，zi，k]T)，N是高分辨率多光谱图像的像素数目。
一般，低分辨率的多光谱图像可以认为是由高分辨率的多光谱图像(假如存在)通过低通滤波和降采样过程得到的[10]。本发明引入高分辨率的多光谱图像与低分辨率多光谱图像之间的观测模型，具体如下所示y＝Hz+u(4)其中u是随机噪声，它的均值是零，协方差矩阵是Cu，且与z是不相关的；H矩阵表示低通滤波和降采样过程。
另外，我们在全色图像和高分辨率多光谱图像之间，引入如下所示的观测模型x＝GTz+v (5)其中v是随机噪声，它的均值为零，协方差矩阵是Cv，且与z是不相关的；G矩阵表示对高分辨率多光谱图像的K个波段做加权平均，权重因子如下所示
gl=ccl/Σl=1K|ccl|---(6)]]>其中ccl表示全色图像与低分辨率多光谱图像的第l波段的相关系数。
由于式(4)和(5)提供的方程数目M×K+N小于待估计的未知量N×K(一般M＜N)，因此不能直接解出z。本发明将从贝叶斯线性模型的角度对上述两个观测模型进行推导，从而得到高分辨率多光谱图像z的估计量。
在估计高分辨率的多光谱图像时，观测模型(4)和(5)可以按照如下的形式联立成一个贝叶斯线性模型yx=HGTz+uv---(11)]]>2、应用贝叶斯高斯-马尔科夫定理，进行线性估计假设数据由贝叶斯线性模型来描述＝Aθ+w(7)其中是L×1数据矢量，A是已知的L×p观测矩阵，θ是p×1的随机矢量。θ的现实是要估计的，它的均值是E(θ)，协方差矩阵是Cθ。w是L×1的随机矢量，它的均值是零，协方差矩阵是Cw，且与θ是不相关的。
首先假定θ的估计量 可通过数据集按下式求得， 选择加权系数aj和bi来使贝叶斯均方误差(Bayesian mean square error，BMSE)BMSE(θ^)=E[(θ-θ^)2]---(9)]]>最小，得到的估计量称为线性最小均方误差(Linear minimum mean square error，LMMSE)估计量(贝叶斯高斯-马尔可夫定理[11])，如下所示 此时BMSE(θ^l)=[(C0+ATCw-1A)-1]ii.]]>由于通常情况下θ不能完美地表示为j的线性组合，所以LMMSE估计量不是最佳的，但是在实际中它是相当有用的，因为它具有闭合形式的解，并且只与均值和协方差有关。
对于贝叶斯线性模型(11)，应用(10)式中的LMMSE估计量可以得到高分辨率的多光谱图像的估计(这里，令n＝[uT，vT]T)
z^=E(z)+CzHGTT(HGTGzHGTT+Gn)-1y-E(y)x-E(x)---(12)]]>其中的统计参数估计如下在式(12)中要得到估计量 需要知道高分辨率多光谱图像的均值E(z)和协方差Cz。这里假设高分辨率多光谱图像的像素之间是相互独立的，因此整个图像的均值和协方差可以由各个像素的均值和协方差组成[12]，具体如下所示E(z)=[E(z1)T,E(z2)T,···,E(zi)T,···,E(zN)T]---(13)]]> 其中E(z)是使用双线性插值(B)后的多光谱图像，具体为E(z)＝B(y) (15)为了估计Cz，我们用矢量量化算法对上述矢量E(z1)根据欧式距离进行分类，计算每一类矢量集的协方差矩阵，并把它作为类中各矢量对应的协方差矩阵。
另外，式(12)中的E(y)的估计是通过对插值图像低通滤波和降采样(H)得到的，E(y)＝H(E(z))＝H(B(y)) (16)而E(x)的估计是通过对全色波段图像进行低通滤波和降采样，然后再双线性插值得到的，E(x)＝B(H(x))(17)在实际的计算中，如果协方差矩阵维数比较大，(12)式中的求逆运算将会产生困难，因此，在这种情况下，本发明将图像分割成许多图像小块再进行估计。
上述的贝叶斯线性估计方法对于融合前图像的波段数目没有限制，低分辨率的多光谱图像y的波段数目可以大于3，而高分辨率的全色图像既可以是单波段也可以是多波段的。因此，贝叶斯线性估计方法既适用于多波段的多光谱图像和单波段的全色图像的融合，也适用于多波段的超光谱图像和多波段的多光谱图像的融合。
图1为基于贝叶斯线性估计的遥感图像融合方法的几何解释。
图2为Landsat ETM+的全色图像和多光谱图像融合实验结果。其中，图2(a)是30m空间分辨率的多光谱图像，图2(b)是30m空间分辨率的全色图像，图2(c)是120m空间分辨率的多光谱图像，图2(d)是HIS方法，图2(e)是PCA方法，图2(f)是小波变换方法，图2(g)是Nishii方法，图2(h)是Hardie方法，图2(i)是贝叶斯方法图3为SPOT的全色图像和TM的多光谱图像融合实验结果。其中，图3(a)是30m空间分辨率的多光谱图像，图3(b)是30m空间分辨率的全色图像，图3(c)是150m空间分辨率的多光谱图像，图3(d)是HIS方法，图3(e)是PCA方法，图3(f)是小波变换方法，图3(g)是Nishii方法，图3(h)是Hardie方法，图3(i)是贝叶斯方法具体实施方式
以下通过实施例，进一步对发明中的各个组成加以描述。
1设置观测模型低分辨率的多光谱图像可以认为是由高分辨率的多光谱图像(假如存在)通过低通滤波和降采样过程得到的[10]，该观测模型如下所示y＝Hz+u(4)其中u是随机噪声，它的均值是零，协方差矩阵是Cu，且与z是不相关的；H矩阵表示低通滤波和降采样过程。
在全色图像和高分辨率多光谱图像之间，存在着如下所示的观测模型x＝GTz+v (5)其中v是随机噪声，它的均值为零，协方差矩阵是Cv，且与z是不相关的；G矩阵表示对高分辨率多光谱图像的K个波段做加权平均，权重因子如下所示gl=ccl/Σl=1K|ccl|---(6)]]>其中ccl表示全色图像与低分辨率多光谱图像的第l波段的相关系数。
2将观测模型联系成贝叶斯线性模型在估计高分辨率的多光谱图像时，观测模型(4)和(5)可以按照如下的形式联立成一个贝叶斯线性模型yx=HGTz+uv---(11)]]>3应用贝叶斯线性估计对于一个贝叶斯线性模型(11)，应用(10)式中的LMMSE估计量可以得到高分辨率的多光谱图像的估计(这里，令n＝[uT，vT]T)z^=E(z)+CzHGTT(HGTGzHGTT+Gn)-1y-E(y)x-E(x)---(12)]]>4估计统计参数假设高分辨率多光谱图像的像素之间是相互独立的，因此整个图像的均值和协方差可以由各个像素的均值和协方差组成[12]，如下所示E(z)＝[E(z1)T，E(z2)T，…，E(zi)T，…，E(zN)T] (13) 其中E(z)是使用双线性插值(B)后的多光谱图像，E(z)＝B(y)(15)为了估计Cz，我们用矢量量化算法对上述矢量E(zi)根据欧式距离进行分类，计算每一类矢量集的协方差矩阵，并把它作为类中各矢量对应的协方差矩阵。
另外，式(12)中的E(y)的估计是通过对插值图像低通滤波和降采样(H)得到的，E(y)＝H(E(z))＝H(B(y)) (16)而E(x)的估计是通过对全色波段图像进行低通滤波和降采样，然后再双线性插值得到的，如下所示E(x)＝B(H(x)) (17)对本发明的上述方法，进行了仿真计算。具体的仿真条件如下(1)Landsat 7 ETM+的全色图像和多光谱图像；(2)SPOT的全色图像和TM的多光谱图像。
图2显示的是Landsat 7 ETM+传感器在2000年6月14日拍摄的上海地区的多光谱图像和全色图像。其中，全色图像具有15m的空间分辨率，而多光谱图像具有30m的空间分辨率，这里以第3，2和1波段分别作为RGB通道。
由于Landsat 7 ETM+不提供15m分辨率的真实多光谱图像用来比较，因此我们很难评价各种方法得到的融合效果。为了解决这个问题，我们将全色图像和多光谱图像的空间分辨率分别退化到30m和120m。通过对上述退化图像进行融合，并将融合结果与原30m分辨率的多光谱图像进行比较。
图3显示了在1995年10月26日拍摄的Hanoi地区的SPOT卫星的全色图像和TM的多光谱图像(http//earth.esa.int/mcities/images/cases)。本发明将SPOT的全色图像从10m的分辨率退化到30m，TM的多光谱图像从30m的分辨率退化到150m。这里以TM的第3，2，和1波段分别作为RGB通道。
本发明所提出的方法将和以下方法做比较(1)传统的HIS方法、PCA方法；
(2)小波变换方法，选用4阶Daubechies小波基，分解的层次是3层。
(3)Nishii方法和Hardie方法。
实验结果如下1、Landsat 7 ETM+的全色图像和多光谱图像从视觉效果上分析，图2中基于HIS方法和PCA方法的融合图像显著地改变了真实多光谱图像的光谱特性，其融合效果明显劣于其它图像融合方法，因此在以下的基于统计参数的定量比较中，不再考虑HIS方法和PCA方法。
小波变换方法导致右边的水体中出现了一些斑点。Nishii方法和Hardie方法的融合图像很模糊。这里需要指出的是低分辨率多光谱图像和全色图像的相关系数分别是0.57、0.18和-0.20。贝叶斯方法的融合图像比较接近于真实的多光谱图像，不仅增强了空间细节，而且很好地保持了原多光谱图像的光谱特性。
以下部分我们定量地比较小波变换方法、Nishii方法和Hardie方法以及贝叶斯方法在增强空间细节和保持光谱特性上的统计参数变化情况。
为了衡量上述方法增强空间细节的能力，我们计算融合图像减去120m分辨率的多光谱图像得到的差值图像在各个波段上的标准差(SDD)，如表1所示。其中第一列显示了真实的多光谱图像(30m分辨率)的SDD参数。在表1中，Nishii方法和Hardie方法在各通道上的SDD参数都远远小于真实图像的值，说明Nishii方法和Hardie方法在本实验中不能有效地增强空间细节。小波变换方法在各通道上的SDD参数是相同的，这说明小波变换方法主要取决于分解层次，而对各通道上的具体情况不作针对性的处理。贝叶斯方法的SDD参数在各通道上比较接近于真实图像的值，这说明通过高分辨率多光谱图像的协方差估计来决定融合图像在各通道上的幅度的做法是合理的。
表1 各种方法增强空间细节的统计参数(Landsat ETM+)
为了衡量上述方法保持光谱特性的能力，我们采用如下的统计参数(1)峰值信噪比[9](PSNR)用来衡量图像灰度的峰值和两幅图像的误差之间的比率，定义如下PSNR＝20×log10(b/rms)(20)其中b表示图像灰度的峰值，rms是两幅图像的误差均方根，峰值信噪比的单位是db。一般峰值信噪比越大，两幅图像之间的差异越小。
(2)相关系数(CC)的定义如下C(f,g)=Σ[(f(i,j)-f‾)×(g(i,j)-g‾)]Σ[(f(i,j)-f‾)2]×Σ[(g(i,j-g‾))2]---(21)]]>其中f(i，j)和g(i，j)表示图像的灰度，f和g是图像的均值。一般相关系数越高，两幅图像越相似。峰值信噪比和相关系数在融合图像和真实的多光谱图像的各波段上分别计算。
(3)相对全局误差[5](ERGAS)如下所示，RASE参数越低，光谱的扭曲程度越小。
ERGAS=100hl1KΣk=1Krms(K)/mz(K)---(22)]]>其中，l和h是融合前和融合后多光谱图像的分辨率(这里分别取120m和30m)，rms(K)表示融合图像和真实的多光谱图像在各波段上的误差均方根，mz(K)是真实的多光谱图像在各波段上的均值。一般相对全局误差越小，说明融合图像的光谱特性保持的越好。
上述统计参数如表2所示。贝叶斯方法在各通道上都具有最高的峰值信噪比，和最大的相关系数，而相对全局误差是最小的，说明贝叶斯方法的融合图像和真实的多光谱图像之间的差异最小，因此贝叶斯方法在光谱特性的保持上优于小波变换方法、Nishii方法和Hardie方法。
表2 各种方法保持光谱特性的统计参数(Landsat ETM+)
2、SPOT的全色图像和TM的多光谱图像从视觉效果上分析，图3中的HIS方法和PCA方法显著地改变了多光谱图像的光谱特性，尤其是PCA方法将河流的色彩转换。Nishii方法和Hardie方法的融合图像比较模糊。这里，低分辨率多光谱图像和全色图像的相关系数分别是0.52、0.35和0.29。小波变换方法和贝叶斯方法的融合效果相对比较好。
在表3中，Nishii方法和Hardie方法的SDD参数都很小，说明它们对于空间分辨率的提高确实很有限。小波变换方法在各通道上具有相同的SDD参数，仍然没有对不同通道上的具体情况作有针对性的处理。贝叶斯方法在各通道上的SDD参数和真实多光谱图像的值几乎相同。
另外，在表4中，贝叶斯方法在各通道上都具有最高的峰值信噪比和最大的相关系数，而相对全局误差是最小的，说明贝叶斯方法的融合图像和真实的多光谱图像之间的差异较小，因此，可以说贝叶斯方法的性能优于小波变换方法、Nishii方法和Hardie方法。
表3 各种方法增强空间细节的统计参数(SPOT和TM)
表4 各种方法保持光谱特性的统计参数(SPOT和TM)
综上所述，贝叶斯方法解决了Nishii方法和Hardie方法依赖于多光谱图像和全色图像的相关系数的问题，同时实验结果证明了贝叶斯方法的融合性能明显优于HIS方法、PCA方法和小波变换方法。
最后，特别值得指出的是，一些方法如小波变换方法，由于参数设置和操作人员的不同，会导致融合结果不同，而本发明所提出的方法将自动设置参数，在不需要人为干预的情况下仍能取得较好的融合效果。
参考文献[1]J W Carper，T M Lillesand，R W Kiefer.The use of intensity-hue-saturation transformationfor merging SPOT panchromatic and multispectral image data[J].PhotogrammetricEngineering and Remote Sensing，1990，56459-467. V K Shettigara.A generalized component substitution technique for spatial enhancement ofmultispectral images using a higher resolution data set[J].Photogrammetric Engineering andRemote Sensing，1992，58561-567. P S Chavez，S C Sides，J A Anderson.Comparison of three different methods to mergemulti-resolution and multispectral daraLandsat TM and SPOT panchromatic[J].Photogrammetric Engineering and Remote Sensing，1991，57(3)295-303. J X Otazu，O Fors，et al.Multiresolution-based image fusion with additive waveletdecomposition[J].IEEE Transactions on Geoscience and Remote Sensing，1999，371204-1211. M A González，J L Saleta，R G Catalán，et al.Fusion of mnltispectral and panchromaticimages using improved IHS and PCA mergers based on wavelet decomposition[J].IEEETransactions on Geoscience and Remote Sensing，2004，23(18)1291-1299. J C Price.Combining panchromatic and multispectral imagery from dual resolution satelliteinstruments[J].Remote Sensing of Environment，1987，21119-128. C K Munechika，J S Warnick，C Salvaggio，et al.Resolution enhancement of multispectralimage data to improve classification accuracy[J].Photogrammetric Engineering and RemoteSensing，1993，5967-72. R Nishii，S Kusanobu，S Tanaka.Enhancement of low resolution image based on highresolution bands[J].IEEE Transactions on Geoscience and Remote Sensing，1996，341151-1158. R C Hardie，M T Eismann，G L Wilson.MAP estimation for hyperspectral image resolutionenhancement using an auxiliary sensor.IEEE Transactions on image processing，2004，13(9)1174-1184. M T Eismann，R C Hardie.Application of stochastic mixing model to hyperspectralresolution enhancement[J].IEEE Transactions on Geoscience and Remote Sensing，2004，42(9)1924-1933. S M Kay，Fundamentals of statistical signal processingEstimation theory[M].EnglewoodsCliffs，NJPrentice-Hall，1993391-392. H Eves.Elementary matrix theory[M].New YorkDover，1966107. J Zhou，D L Civco，J A Silander.A wavelet transform method to merge Landsat TM andSPOT panchromatic data[J].International Journal of Remote Sens.，1998，19(4)743-757.
权利要求
1.一种基于贝叶斯线性估计的遥感图像融合方法，其特征在于引入高分辨率多光谱图像和低分辨率多光谱图像之间的观测模型，以及高分辨率多光谱图像和全色图像之间的观测模型，并将上述两个观测模型联立成一个贝叶斯线性模型；应用贝叶斯高斯-马尔科夫定理，计算得到线性最小均方误差LMMSE意义上的高分辨率多光谱图像的估计。
2.根据权利要求1所述的基于贝叶斯线性估计的遥感图像融合方法，其特征在于所述的高分辨率多光谱图像和低分辨率多光谱图像之间的观测模型如下y＝Hz+u (4)其中u是随机噪声，它的均值是零，协方差矩阵是Cu，且与z是不相关的；H矩阵表示低通滤波和降采样过程；y为低分辨率的多光谱图像排成一维的列矢量y=[y1T,y2T,...,yjT,...,yMT]T---(2)]]>其中yj表示低分辨率的多光谱图像在空间位置j上的像素值(有K个波段，yj＝[yj，1，yj，2，…，yj，K]T)，而M是低分辨率多光谱图像的像素数目；z为高分辨率的多光谱图像排成一维的列矢量z=[z1T,z2T,...,ziT,...,zNT]T---(3)]]>其中zi表示高分辨率的多光谱图像在空间位置i上的像素值(有K个波段，zi＝[zi，1，zi，2，…，zi，K]T)，N是高分辨率多光谱图像的像素数目；所述高分辨率多光谱图像和全色图像之间的观测模型如下x＝GTz+v(5)其中v是随机噪声，它的均值为零，协方差矩阵是Cv，且与z是不相关的；G矩阵表示对高分辨率多光谱图像的K个波段做加权平均，权重因子如下所示gl=ccl/Σl=1K|ccl|---(6)]]>其中ccl表示全色图像与低分辨率多光谱图像的第l波段的相关系数；x为全色图像排成的一维矢量x＝[x1，x2，…，xi，…，xN]T(1)其中xi表示全色图像在空间位置i上的像素值，而N是全色图像的像素数目；观测模型(4)和(5)按照如下的形式联立成一个贝叶斯线性模型yx=HGTz+uv---(11)]]>
3.根据权利要求2所述的基于贝叶斯线性估计的遥感图像融合方法，其特征在于高分辨率多光谱图像的估计式如下z^=E(z)+CzHGTTHGTCzHGTT+Cn-1y-E(y)x-E(x)---(12)]]>其中，E(z)和Cz分别为高分辨率多光谱图像的均值的协方差，它们由各个像素的均值和协方差组成，具体如下所示E(z)＝[E(z1)T，E(z2)T，…，E(zi)T，…，E(zN)T](13) 其中E(z)是使用双线性插值(B)后的多光谱图像，具体为E(z)＝B(y)(15)为了估计Cz，用矢量量化算法对上述矢量E(zi)根据欧式距离进行分类，计算每一类矢量集的协方差矩阵，并把它作为类中各矢量对应的协方差矩阵；式(12)中的E(y)的估计是通过对插值图像低通滤波和降采样(H)得到的E(y)＝H(E(z))＝H(B(y)) (16)而E(x)的估计是通过对全色波段图像进行低通滤波和降采样，然后再双线性插值得到的E(x)＝B(H(x)) (17)。
全文摘要
本发明属于遥感图像处理技术领域，具体为一种基于统计参数估计的遥感图像融合方法。该方法引入高分辨率多光谱图像和低分辨率多光谱图像之间的观测模型，以及高分辨率多光谱图像和全色图像之间的观测模型，并将上述两个观测模型联立成一个贝叶斯线性模型。通过应用贝叶斯高斯—马尔科夫定理，得到线性最小均方误差意义上的高分辨率多光谱图像的估计。本发明不仅能够增强空间细节，而且很好地保持了光谱特性，其性能优于传统的HIS方法、PCA方法和小波变换方法以及现存的基于统计参数估计的Nishii方法和Hardie方法。新方法可为改善遥感图像的目视判读精度，提高信息清晰度和可靠性上提供新的有效的技术支持。
文档编号G01S7/48GK1808181SQ20061002411
公开日2006年7月26日 申请日期2006年2月23日 优先权日2006年2月23日
发明者葛志荣, 王斌, 张立明 申请人:复旦大学