一种用于时空融合的地理加权空间混合分解方法

本发明涉及遥感影像融合技术领域，尤其是涉及一种用于时空融合的地理加权空间混合分解方法。

背景技术：

landsat和terra\aqua卫星是目前广泛用于全球观测的卫星。因受到技术水平和造价成本等条件的限制，其获取的landsat和modis数据的时间和空间分辨率之间相互制约。自landsat系列和terra\aqua卫星发射以来，landsat和modis遥感数据广泛地应用于监测全球地表变化。但是受到技术水平和造价成本等条件的限制，单个卫星获取的数据无法同时达到高时间和高空间分辨率的要求，从而无法满足地表实时精细监测的需求。具体地，landsat获取的数据的空间分辨率为30m，但重访周期约为16天；terra\aqua卫星获取的modis数据的空间分辨率虽为500m，但每天可获得至少一景。为获取满足应用需求的高时间和空间分辨率的遥感数据，时空融合技术应运而生。目前常用的时空融合方法主要分为两大类：基于空间加权的方法和基于空间混合分解的方法。其中，基于空间混合分解的方法主要包括：以混合分解为基础的数据融合(unmixing-baseddatafusion，ubdf)、遥感数据时空融合方法(spatialandtemporaldatafusionapproach，stdfa)和基于虚拟数据对的空间混合分解方法(virtualimagepair-basedspatio-temporalfusionwithspatialunmixing，vipstf-su)。基于空间混合分解的方法因其明确的数学表达和物理含义而受到广泛的关注和研究。

相比于其他方法，基于空间混合分解的方法对已知信息要求较低，在数据缺乏的区域具有很好的应用价值，同时能够最大限度地利用已知的高空间分辨率信息。一般来说，遥感数据具有空间非平稳性的特点。在一定空间范围内，同一类地物反射率值可能存在较大的变化，像元灰度值间的相关性也随空间位置发生变化。依据地学第一定律，空间距离越近的像元相关性越大。因此，与中心像元距离更近的邻域像元在空间混合分解过程中应当发挥更大的作用。然而，现有的空间混合分解方法没有考虑邻域内同一类地物反射率空间变异性的影响，使得周围邻域像元对中心目标像元的混合分解过程施加了同等的影响。该问题的存在阻碍了基于空间混合分解的时空融合方法的预测精度，限制了其在异质性较强的区域的应用。

技术实现要素：

本发明的目的就是为了克服上述现有技术存在的缺陷而提供一种用于时空融合的地理加权空间混合分解方法。

本发明的目的可以通过以下技术方案来实现：

一种用于时空融合的地理加权空间混合分解方法，该方法包括如下步骤：

s1：根据聚类有效性指数，确定已知的邻近时刻高空间分辨率影像的最优分类数，并依该最优分类数进行分类，获取高空间分辨率分类图。

s2：依次访问预测时刻低空间分辨率数据中每一个低分辨率像元，以该低分辨率像元为中心建立一个窗口，计算窗口内各邻域像元与中心像元之间的空间距离，获取空间距离矩阵，以相应的空间距离计算权重，获取地理加权矩阵。

s3：将步骤s1中得到的高空间分辨率分类图退化至与预测时刻低空间分辨率数据相同的分辨率，计算窗口内各类地物占比矩阵，并对各低分辨率像元构造目标函数，最小化该目标函数，获取各低分辨率像元内各类地物反射率值。

s4：根据分类图和各类地物反射率值构建融合影像。

进一步地，步骤s1中，对不同的分类数c计算聚类有效性指数xb(c)的值，最小的聚类有效性指数值所对应的分类数即为已知的邻近时刻高空间分辨率影像的最优分类数。聚类有效性指数的计算式为：

其中，c为待定的分类数，s为已知的邻近时刻高空间分辨率影像的像元的个数，m为模糊指数，yi为第i个像元对应的光谱特征向量，vc为第c类聚类中心的光谱特征向量，vk为第k类聚类中心的光谱特征向量，c≠k，uci为第i个像元中第c类的隶属度。

进一步地，步骤s2中，根据双重平方函数，以相应的空间距离计算权重，获取的地理加权矩阵为：

式中，b为双重平方函数中的带宽参数，dij为窗口内第j个邻域像元与中心像元i之间的空间距离，wij为在低分辨率像元i空间混合分解时，第j个邻域像元对其施加的权重。

进一步地，步骤s3中，依据加权最小二乘方法及空间混合分解的基本原理对各低分辨率像元构造目标函数。所述目标函数设有追加的约束项，用于与其他空间混合分解模型进行耦合。耦合地理加权模型的广义目标函数表达式为：

式中，n为窗口内低分辨率像元个数，wij为第j个邻域像元对中心像元i的混合分解施加的权重，ei为需要求解的中心像元各类地物反射率向量，qj为窗口内第j个邻域像元的反射率，pj为窗口内第j个邻域像元的各类地物占比行向量，l为广义目标函数中的约束项，α为权衡参数。进一步地，窗口内第j个邻域像元的各类地物占比行向量pj可选择通过软分类方法获取的分类图计算得到。

进一步地，对广义目标函数中的约束项包括但不限于采用砖块效应消除技术，构成基于砖块效应消除的地理加权空间混合分解方法。

本发明提供的用于时空融合的地理加权空间混合分解方法，相较于现有技术至少包括如下有益效果：

一、本发明更准确地考虑了空间混合分解过程中邻域像元对中心像元的影响，提高了融合影像的精度：本发明方法扩展了经典的空间混合分解方法，在空间混合分解模型中创新性地加入了利用双重平方函数量化得到的地理加权矩阵，充分顾及了邻域内同一地物反射率空间变异性的影响，能够有效地恢复地物分布空间异质性，提高时空融合精度；

二、本发明方法具有很好的普适性和扩展性：本发明方法对现有的空间混合分解模型进行扩展，未增加原有方法的输入数据，可直接应用于目前任一种空间混合分解方法，且不会增加现有模型的复杂度，对未来可能提出的新的空间混合分解方法也有很高的应用价值。

附图说明

图1为实施例中用于时空融合的地理加权空间混合分解方法的流程示意图；

图2为实施例仿真实验中异质区域的结果图，其中(a1)为采用原始的ubdf的影像融合结果，(b1)为采用本发明的ubdf-gw的影像融合结果，(c1)为采用原始的stdfa的影像融合结果，(d1)为采用本发明的stdfa-gw的影像融合结果，(e1)为采用原始的vipstf-su的影像融合结果，(f1)为采用本发明的vipstf-su-gw的影像融合结果，(g1)为参考影像；(a2)为采用ubdf-fcm的影像融合结果，(b2)为采用本发明的ubdf-fcm-gw的影像融合结果，(c2)为采用stdfa-fcm的影像融合结果，(d2)为采用本发明的stdfa-fcm-gw的影像融合结果，(e2)为采用vipstf-su-fcm的影像融合结果，(f2)为采用本发明的vipstf-su-fcm-gw的影像融合结果，(g2)为参考影像。

图3为实施例仿真实验中变化区域的结果图，其中(a1)为采用原始的ubdf的影像融合结果，(b1)为采用本发明的ubdf-gw的影像融合结果，(c1)为采用原始的stdfa的影像融合结果，(d1)为采用本发明的stdfa-gw的影像融合结果，(e1)为采用原始的vipstf-su的影像融合结果，(f1)为采用本发明的vipstf-su-gw的影像融合结果，(g1)为参考影像；(a2)为采用ubdf-fcm的影像融合结果，(b2)为采用本发明的ubdf-fcm-gw的影像融合结果，(c2)为采用stdfa-fcm的影像融合结果，(d2)为采用本发明的stdfa-fcm-gw的影像融合结果，(e2)为采用vipstf-su-fcm的影像融合结果，(f2)为采用本发明的vipstf-su-fcm-gw的影像融合结果，(g2)为参考影像。

具体实施方式

下面结合附图和具体实施例对本发明进行详细说明。显然，所描述的实施例是本发明的一部分实施例，而不是全部实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动的前提下所获得的所有其他实施例，都应属于本发明保护的范围。

实施例

如图1所示，本发明涉及一种用于时空融合的地理加权空间混合分解方法，该方法考虑了空间混合分解过程中不同距离的邻域像元对中心像元的影响大小，具体包括如下步骤：

步骤一、根据聚类有效性指数(clustervalidityindexofxieandbeni，xb)，确定已知的邻近时刻高空间分辨率遥感影像的最优分类数。xb指数的计算方式为：

其中，c为待定的分类数，s为邻近时刻高空间分辨率遥感影像的像元的个数，m为模糊指数(通常取2)，yi是第i个像元对应的光谱特征向量，vc为第c类聚类中心的光谱特征向量，vk为第k类聚类中心的光谱特征向量，c≠k；uci为第i个像元中第c类的隶属度。对不同的分类数c计算xb(c)的值，最小的值所对应的分类数即为最优分类数。在确定的最优分类数下，对该数据进行非监督分类，获取分类图。

步骤二、依次访问预测时刻低空间分辨率数据中每一个低分辨率像元，并以该像元为中心建立一个窗口，该像元即称为中心像元。计算窗口内各邻域像元与中心像元之间的空间距离，得到空间距离矩阵。根据双重平方函数，以相应的距离计算权重，得到地理加权矩阵：

其中，b为双重平方函数中的带宽参数，可取为空间混合分解窗口对角线长度的一半。dij为窗口内第j个邻域像元与中心像元i之间的空间距离。wij即为在低分辨率像元i空间混合分解时，第j个邻域像元对其施加的权重。

步骤三、将步骤一中得到的高空间分辨率分类图退化至与预测时刻低空间分辨率数据相同的分辨率，计算窗口内各类地物占比矩阵p。在现有的空间混合分解模型基础上，对各低分辨率像元构造加权目标函数，生成新的su-gw模型。最小化该目标函数，获取各低分辨率像元内各类地物反射率值。对各低分辨率像元构造目标函数为：

式中，n为窗口内低分辨率像元个数，pj为窗口内第j个邻域像元的各类地物占比向量，wij为第j个邻域像元对中心像元i的混合分解施加的权重。在ubdf方法中，ei为需要求解的中心像元各类地物反射率，qj为窗口内第j个邻域像元的反射率。在stdfa方法中，ei为需要求解的已知时刻与预测时刻间中心像元各类地物反射率变化量，qj为窗口内第j个邻域像元相应的反射率变化量。vipstf-su方法中，ei为需要求解的虚拟时刻与预测时刻间中心像元各类地物反射率变化量，qj为窗口内第j个邻域像元相应的反射率变化量。

进一步地，目标函数可继续追加约束项，实现与其他扩展模型的耦合。即在已有的其他的空间混合分解方法上耦合该地理加权模型，得到的广义的地理加权空间混合分解方法的目标函数：

其中，pj可通过软分类方法获取的分类图计算得到，构成基于软分类的地理加权空间混合分解方法。l为对目标函数追加的约束项，可采用砖块效应消除技术，构成基于砖块效应消除的地理加权空间混合分解方法。亦可将软分类方法和约束项同时纳入目标函数，得到同时基于软分类且含有约束的地理加权空间混合分解方法。α为权衡参数。通过最小化该目标函数，获得各低分辨率像元内各类地物反射率值。

步骤四、根据分类图和各类地物反射率值构建融合影像。

为了验证本发明方法的有效性，本实施例采用本发明方法预测融合影像。基于空间混合分解的时空融合方法中包括三种常用经典方法ubdf、stdfa、vipstf-su，本实施例将本发明方法分别用于上述三种方法上。此外，本实施例还将采用基于软分类(即fuzzyc-means，fcm)的地理加权空间混合分解方法作为扩展的空间混合分解方法的代表。以下简称表示的含义为：su-gw：基于地理加权空间混合分解方法；su-fcm：基于软分类的空间混合分解方法；su-fcm-gw：基于软分类的地理加权空间混合分解方法。本实施例同时将su-gw和su-fcm-gw的预测结果分别与现有的原始的空间混合分解方法(su)和扩展的空间混合分解方法(su-fcm)进行比较。两个测试区域皆位于澳大利亚的新南威尔士州北部(异质区域与变化区域)。两个区域的融合影像结果分别如图2和图3所示，第一行为su及su-gw方法整个区域的预测结果，第二行为相应的局部放大的子图，第三行为su-fcm及su-fcm-gw方法整个区域的预测结果，第四行为相应的局部放大的子图。

根据图2、3可知，原始的空间混合分解方法(su)预测的结果存在明显的砖块效应且光谱畸变较为严重；基于软分类的空间混合分解方法(su-fcm)借助软分类技术更好地描述了低分辨率像元内的类内光谱差异，预测结果精度较好，但仍存在砖块效应和光谱畸变现象。在本发明方法中，更加准确地考虑了空间混合分解过程中邻域像元对中心像元的影响，充分顾及了邻域内同一地物反射率的空间变异性，在原始方法和基于软分类的方法上都体现出对砖块效应及光谱畸变较好的修复效果。因此，本发明的结果在视觉展示上有很大的提高。

采用均方根误差(rootmeansquareerror，rmse)和相关系数(correlationcoefficient，cc)评价指标对各方法获取的融合影像进行精度评价，如表1所示。其中rmse度量预测影像与参考影像的差异性，其值越大表明预测影像越偏离参考影像；cc反映预测影像与参考影像之间的相关度，其值越大表示预测影像与参考影像越接近。

表1影像融合结果精度评价

从表1的客观评价结果可以看出，本发明方法的精度相较于非地理加权模型有明显的提高，cc和rmse值均表明本发明方法能够得到更接近真实情况的融合影像。综上所述，本发明的用于时空融合的地理加权空间混合分解方法在视觉展示和精度评价上都有明显优势，得到的融合影像能较好地保持地物的光谱和空间信息，是一种可行有效的时空融合方法。

以上所述，仅为本发明的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的工作人员在本发明揭露的技术范围内，可轻易想到各种等效的修改或替换，这些修改或替换都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应以权利要求的保护范围为准。