不同空间分辨率的多光谱遥感图像融合方法

摘要

本发明属于遥感图像处理领域，具体涉及一种不同空间分辨率的多光谱遥感图像融合方法。该方法首先将高空间分辨率多光谱图像分成n类地物，计算高空间分辨率多光谱图像中各类地物的光谱矢量均值，然后对高空间分辨率多光谱图像进行混合像素分解，利用高空间分辨率多光谱图像中每一个像素中各类地物所占的百分此表示低空间分辨率多光谱图像，求解低空间分辨率多光谱图像中各类地物的光谱矢量均值，最后计算融合图像的每个像素的像素值，生成融合图像。本发明能够将高空间分辨率多光谱图像和低空间分辨率的多光谱图像进行融合，使融合后的图像更加清晰。

说明书

技术领域

本发明属于遥感图像处理领域，具体涉及一种不同空间分辨率的多光谱遥感图像融合方法。

背景技术

目前各种融合算法一般是通过建立适当的融合模型来解决融合过程中已知条件不足所导致的方程欠定性问题，如IHS变换方法，首先将多光谱图像通过插值算法采样到与全色图像相同的分辨率，然后将图像的谱信息与强度信息相分离，用全色图像代替强度信息与谱信息反变换成融合图像。在融合过程中融合图像的谱信息(在图像中表现为图像的色彩)是通过插值算法对多光谱图像进行插值得到的，因此融合图像谱信息的准确度在一定程度上由插值算法决定。再如高通滤波或小波融合方法，它们将高空间分辨率图像的高频信息从图像中分离出来，迭加在低空间分辨率的图像上，但是事实上，不同波段图像的高频信息是不相同的，因此这种迭加方法必然会对融合结果的准确性造成一定程式度的影响。事实上这些融合的方法都同样的忽略掉了一个问题，即图像上每一个像元与其周围的临近像元有着高度的相关性，而只是单独的考虑每一个像元，因此在理论上这些算法是不完备的，也就是说这些方法即使在理论上讲也只能产生与真实图像相近似的融合图像。而混合像元分解技术能够将像元分解成主要特征地物在其中的百分含量，表达了图像各点与主要特征地物之间的关系，因此如果能将这种关系运用到图像融合中，就可以在一定程度上解决融合方程的欠定性问题。

目前IKONOS、QuickBird等高空间分辨率遥感卫星所携带的传感器可以采集较高空间分辨率的多光谱图像，但是只具有三到四个波段的多光谱信息，而TM，Aster、Hyperion的数据具有较多的光谱信息，但其分辨率却很低。因此可以利用IKONOS、QuickBird的中分辨率多光谱图像对具有更多光谱信息的低空间分辨率图像进行锐化。

发明内容

本发明的目的在于针对现有技术中存在的问题，提供一种基于混合像元分解技术的不同空间分辨率的多光谱遥感图像融合方法，从而能够将高空间分辨率多光谱图像和低空间分辨率的多光谱图像进行融合，使融合后的图像更加清晰。

为实现上述目的，本发明的技术方案如下：一种不同空间分辨率的多光谱遥感图像融合方法，包括如下步骤：

(1)对同一地区的高空间分辨率的多光谱图像以及低空间分辨率的多光谱图像进行预处理去噪、配准；

(2)对高空间分辨率多光谱图像进行非监督分类，将其分成n类地物；

(3)计算高空间分辨率多光谱图像中各类地物的光谱矢量均值；

(4)将步骤(3)中得到的高空间分辨率多光谱图像中各类地物的光谱均值作为混合像元分解的端元，对高空间分辨率多光谱图像进行混合像素分解，计算高空间分辨率多光谱图像中每一个像素中各类地物所占的百分比；

(5)利用高空间分辨率多光谱图像中每一个像素中各类地物所占的百分比表示低空间分辨率多光谱图像，求解低空间分辨率多光谱图像中各类地物的光谱矢量均值；

(6)利用步骤(4)中得到的高空间分辨率多光谱图像中每一个像素中各类地物所占的百分比，以及步骤(5)中得到的低空间分辨率多光谱图像中各类地物的光谱矢量均值，计算融合图像的每个像素的像素值，生成融合图像。

进一步，如上所述的不同空间分辨率的多光谱遥感图像融合方法，步骤(2)中所述的非监督分类采用K-means方法，地物的类数n大于200个。

进一步，如上所述的不同空间分辨率的多光谱遥感图像融合方法，步骤(3)中计算高空间分辨率多光谱图像中各类地物的光谱矢量均值Z^A_i的公式如下：

$\overline{Z_i^A}=\frac{1}{{\mathrm{NUM}}_i}·\underset{(x,y)\in {\mathrm{class}}_i}{\Sigma }\mathrm{AI}(x,y)$

其中，i＝1，2…n ，NUM_i为高空间分辨率多光谱图像中第i类地物的像素总数，class_i为第i类地物，AI(x，y)为高空间分辨率多光谱图像中第x行第y列处的光谱矢量。

进一步，如上所述的不同空间分辨率的多光谱遥感图像融合方法，步骤(4)中计算高空间分辨率多光谱图像中每一个像素中各类地物所占的百分比e_i(x，y)的公式如下：

$\mathrm{AI}(x,y)=\underset{i=1}{\overset{n}{\Sigma }}\{e_i(x,y)·\overline{Z_i^A}\}$

其中，AI(x，y)为高空间分辨率多光谱图像中第x行第y列处的光谱矢量，Z^A_i为高空间分辨率多光谱图像中各类地物的光谱矢量均值。

进一步，如上所述的不同空间分辨率的多光谱遥感图像融合方法，步骤(5)中利用高空间分辨率多光谱图像中每一个像素中各类地物所占的百分比表示低空间分辨率多光谱图像的公式如下：

$\mathrm{BI}(x,y)=\frac{1}{{\mathrm{NUM}}_{(x^{\prime },y^{\prime })\in D}}·\underset{i=1,(x^{\prime },y^{\prime })\in D}{\overset{n}{\Sigma }}{e}_i(x^{\prime },y^{\prime })·\overline{Z_i^B}$

其中，BI(x，y)为低空间分辨率多光谱图像中第x行第y列处的光谱矢量，D为低空间分辨率多光谱图像中x，y像素的对应区域，x′，y′为高空间分辨率多光谱图像中的像素坐标，NUM_{(x′，y′)∈D}为D区域中高空间分辨率多光谱图像的像素总数，Z^B_i为低空间分辨率多光谱图像中各类地物的光谱矢量均值。

进一步，如上所述的不同空间分辨率的多光谱遥感图像融合方法，步骤(6)中计算融合图像的每个像素的像素值CI(x，y)的公式如下：

$\mathrm{CI}(x,y)=\underset{i=1}{\overset{n}{\Sigma }}{e}_i(x,y)·\overline{Z_i^B}$

其中，e_i(x，y)为高空间分辨率多光谱图像中每一个像素中各类地物所占的百分比，Z^B_i为低空间分辨率多光谱图像中各类地物的光谱矢量均值。

本发明的有益效果如下：利用本发明提出的方法可以实现不同空间分辨率多光谱(或高光谱)图像的融合，尤其是高空间分辨率多光谱图像和低空间分辨率高光谱图像，进行融合的两景影像即使空间分辨率差异巨大，也不会产生类似其它融合方法的生成的色斑现象，生成的融合图像仍然纹理清晰。

附图说明

图1为本发明具体实施例的方法流程图。

具体实施方式

下面结合附图和实施例对本发明进行详细的描述。

为了描述方便，令A代表某地区的高空间分辨率的多光谱图像，AI(x，y)为第x行第y列处的光谱矢量，则AI(x，y)＝{A₁(x，y)，A₂(x，y)…A_l(x，y)}；B代表同一地区的低空间分辨率的多光谱图像，BI(x，y)为第x行第y列处的光谱矢量，则BI(x，y)＝{B₁(x，y)，B₂(x，y)…B_p(x，y)}。

如图1所示，本发明所提供的不同空间分辨率的多光谱遥感图像融合方法，包括如下步骤：

步骤1)对同一地区的高空间分辨率的多光谱图像以及低空间分辨率的多光谱图像进行预处理去噪、配准。这一步骤使用图像处理技术中惯用的去噪、配准方法即可，此处不再过多描述。

步骤2)对高空间分辨率多光谱图像进行非监督分类。

选取参加图像融合中空间分辨率较高的一幅图像作为A图像，对其进行非监督分类。分类的方法为K-means方法，设置的分类数可以尽可能的高，一般应在200个以上。例如共分成的n类地物：class₁ class₂…class_i…class_n。

步骤3)计算各类地物光谱均值。

计算A图像每一类别地物的光谱矢量均值，即：

$\overline{Z_i^A}=\frac{1}{{\mathrm{NUM}}_i}·\underset{(x,y)\in {\mathrm{class}}_i}{\Sigma }\mathrm{AI}(x,y)$

其中：i＝1，2…n  ，NUM_i为A图中第i类地物的像素总数。

步骤4)对A图像进行混合像素分解。

以计算得到的光谱均值Z^A_i(i＝1…n)作为混合像素分解的端元，对A图像进行混合像素分解。计算A图像中每一像素的光谱组成，即计算A图像中每一个像素中各类地物所占的百分比。则A图像中的每一像素可以表示为：

$\mathrm{AI}(x,y)=\underset{i=1}{\overset{n}{\Sigma }}\{e_i(x,y)·\overline{Z_i^A}\}$

其中e_i(x，y)表示在A图像x行y列处，i类别地物所占百分比。

步骤5)利用e_i(x，y)表示低空间分辨率多光谱图像。

假设B图像的空间分辨率是A图像空间分辨率的m倍，则B图像某像素对应的区域D中包含m²个A图像的像素。由于B图像与A图像象素位置存在这样的对应关系，在求解出A图像每一象素中各类地物的组份后，则可以将B图像进行如下表示：

$\mathrm{BI}(x,y)=\frac{1}{{\mathrm{NUM}}_{(x^{\prime },y^{\prime })\in D}}·\underset{i=1,(x^{\prime },y^{\prime })\in D}{\overset{n}{\Sigma }}{e}_i(x^{\prime },y^{\prime })·\overline{Z_i^B}$

式中D为B图像中x，y像素的对应区域，x′，y′为A图像中的像素坐标，NUM_{(x′，y′)∈D}为D区域中A图像的像素总数，Z^B_i为B图像中各类地物的光谱矢量均值。

步骤6)迭代求解B图像中各地物的光谱矢量。

由于BI(x，y)已知，e_i(x′，y′)在上面的步骤中已经求出，可以利用建立的关于Z^B_i的方程式对各类地物在对应B图像波段的光谱进行求解。每一个B图像的像素可建立一个方程式，因此方程数量非常庞大，且为超定方程(一般情况下，B图像素的个数远远大于地物的种类)，利用迭代方法求算最佳解。

步骤7)生成融合图像。

根据获得的B图像各类地物的光谱矢量均值Z^B_i，将融合图像每一像素可以表示为：

$\mathrm{CI}(x,y)=\underset{i=1}{\overset{n}{\Sigma }}{e}_i(x,y)·\overline{Z_i^B}$

其中，e_i(x，y)表示在A图像×行y列处，i类别地物所占百分比。最后计算得到的CI即为融合图像。

表1四种方法的融合结果的客观评价数据

表1中列出了本发明方法以及IKONOS多光谱法、ASTER多光谱法、HIS融合法这四种方法的融合结果的客观评价数据，从融合结果和客观评价指标进行分析，本发明方法虽在保持光谱特性方面即相关系数与HIS融合法相比略有差距，但是在信息量的融入程度方面远远超过IHS融合法，客观评价指数很高。