顾及非局部特性与时空变化的遥感数据时空定量融合方法

摘要

本发明公开了一种顾及非局部特性与时空变化的遥感数据时空定量融合方法，基于非局部滤波，利用移动窗口技术，首先选取相似像元，用基于小窗口运算的经验型公式筛选相似像元，并在此基础上二次筛选相似像元，以相似像元为中心的邻域矩阵的相似度以及其与中心像元的相对距离衡量权重大小，考虑时空变化，针对区域主要特征，选取相应的方式加权，最终融合得到中心像元的反射率值。本发明利用高空间分辨率数据与高时间分辨率数据的互补信息，顾及非局部特性与时空变化，融合得到既具有高空间有具有高时间分辨率的数据，不仅取得更高的精度，而且更具有实际应用潜力。

说明书

技术领域

本发明属于遥感影像融合领域，涉及一种顾及非局部特性与时空变化的遥感 数据时空定量融合方法，充分利用多时相遥感影像中存在的互补信息以及融合区 域的特征，可以融合得到更可靠的既具有高空间又具有高时间分辨率的数据。

背景技术

高分辨率、长时序、高精度遥感监测对于全球变化研究、资源调查与管理、 环境监测等相关领域具有重要意义。然而，由于传感器设计技术上各指标的制约， 使得其空间分辨率与时间分辨率也相互制衡。为了解决这一问题，产生了相应的 方法，其利用高空间分辨率的卫星数据与高时间分辨率的卫星数据的互补信息， 融合生成既具有高空间分辨率也具有高时间分辨率的数据，即时空融合技术。目 前常用的融合方法有经典的时空自适应反射率融合方法以及基于光谱解混理论 的时空自适应反射率融合方法。两种方法均有各自的优势和局限性，经典的方法 在预测主要是时间变化的区域有优势，而对于异质区域空间上的细节变化往往难 以预测出来；增强型算法在预测空间上的细节变化上有优势，但对于时间变化的 区域预测效果比较差一些。

发明内容

本发明的目的在于，针对现有技术的特点，提供一种基于非局部思想的并顾 及时间与空间变化的时空定量融合方法，考虑时间和空间变化，进行时空定量融 合。

为了更好地捕捉到区域特征，本发明针对两影像对预测，结合非局部滤波， 用基于小窗口运算的经验型公式相似像元筛选方法，并进行二次筛选相似像元， 采用一阶的类似非局部均值权重，基于移动窗口运算，对相似像元加权平均，预 测窗口中心像元的反射率值。

本发明所采用的技术方案是：一种顾及非局部特性与时空变化的遥感数据时 空定量融合方法，其特征在于，包括以下步骤：

步骤1：输入Tm时刻高空间与低空间分辨率影像对和Tn时刻高空间与低空间 分辨率影像对，对输入影像对进行预处理，所述的预处理包括重投影、重采样或 裁剪；

步骤2：Tm时刻高空间与低空间分辨率影像对和Tn时刻高空间与低空间分辨率 影像对，各自利用经验型公式并结合时间差异T_ijk和光谱差异S_ijk筛选相似像元；

步骤3：基于非局部滤波，采用以相似像元为中心的邻域矩阵的相似度以及其与 中心像元的相对距离计算各相似像元的权重；

步骤4：考虑时间变化和空间变化，对时间变化的区域采用经典算法方式融合， 对空间变化的区域采用时间加权的方式融合。

作为优选，步骤2中所述的经验型公式为：

$|\bar{F}(x_i,y_j,B)-\bar{F}(x_{w/2},y_{w/2},B)|<=d*2^{\bar{F}(x_{w/2},y_{w/2},B)}---\left(1\right);$

其中，F(x_i,y_j,B)为基础时间高空间分辨率数据邻近像元的反射率值， F(x_w/2,y_w/2,B)为中心像元反射率，B为波段数，d为自由参数。

作为优选，所述的筛选过程都是在小窗口中进行的，即计算均值。

作为优选，步骤3的具体实现过程为，基于非局部滤波，用基于矩阵的方法 计算各相似像元的权重，采用一阶形式进行空间滤波，并在此基础上加入相对距 离，具体实现方式为:

$W=\exp (-\frac{\bar{S}+\bar{T}}{h})*1/D---\left(2\right);$

其中，和分别为邻域矩阵的光谱差异和时间差异的均值，通过调节滤 波参数h和矩阵窗口来进行试验，最优的h的范围为[0.001，0.1]；D为邻近像 元(x_i，y_j)与中心像元的相对距离。

作为优选，步骤4中所述的对时间变化的区域采用经典算法方式融合，其公 式为：

$F(x_i,y_j,t_k)=\underset{i=1}{\overset{w}{\Sigma }}\underset{j=1}{\overset{w}{\Sigma }}\underset{k=1}{\overset{n}{\Sigma }}{W}_{\mathrm{ijk}}\times (F(x_i,y_j,t_0)+C(x_i,y_j,t_k)-C(x_i,y_j,t_0))---\left(3\right);$

其中，F、C分别代表着高空间分辨率与高时间分辨率数据的反射率， (x_i,y_j)为给定高空间分辨率与高时间分辨率数据像元的位置，t₀为数据获取的 时间，t_k为预测时间，w为搜索窗口大小，n为输入的影像对数；

步骤4中所述的对空间变化的区域采用时间加权的方式融合，其公式为：

F(x_w/2,y_w/2,t_p)＝F(x,y,t₀)+∑w(x,y)*a_w/2*(C(x,y,t_p)-C(x,y,t₀))    (4)； 其特点是计算移动窗口内相似像元的线性回归系数a_w/2，进而计算各基础时间的 时相权重，加权得到最终的融合结果。

本发明的特色在于，基于非局部滤波，进行相似像元的筛选和其权重的计算， 并考虑时空变化，针对区域主要特征，选取相应的方式加权，最终融合得到中心 像元的反射率值。总之，本发明提出的方法可有效的用于遥感影像的时空定量融 合，获得更加准确的预测结果。

附图说明

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

具体实施方式

为了便于本领域普通技术人员理解和实施本发明，下面结合附图及实施例对 本发明作进一步的详细描述，应当理解，此处所描述的实施示例仅用于说明和解 释本发明，并不用于限定本发明。

请见图1，本实施例提供的一种顾及非局部特性与时空变化的遥感数据时空 定量融合方法，用于融合得到中心像元的反射率值；包括以下步骤：

步骤1：输入Tm时刻高空间与低空间分辨率影像对和Tn时刻高空间与低 空间分辨率影像对，对输入影像对进行预处理，所述的预处理包括重投影、重采 样或裁剪；

步骤2：Tm时刻高空间与低空间分辨率影像对和Tn时刻高空间与低空间分 辨率影像对，各自利用经验型公式并结合时间差异T_ijk和光谱差异S_ijk筛选相似像 元；

经验型公式为：

$|\bar{F}(x_i,y_j,B)-\bar{F}(x_{w/2},y_{w/2},B)|<=d*2^{\bar{F}(x_{w/2},y_{w/2},B)}---\left(1\right);$其中，F(x_i,y_j,B)为基础时间高空间分辨率数据邻近像元的反射率值(如果是其他 产品，参数d适当改变)，F(x_w/2,y_w/2,B)为中心像元反射率，B为波段数，d为 自由参数，对给定传感器其值在很小的范围内变动。经过初步的筛选，仍然不可 避免的有些低质量的像元或一些不可靠的像元被筛选进来，因此还要基于筛选出 来的像元进行进一步剔除，引入时间差异T_ijk和光谱差异S_ijk，整个筛选相似像元 的方法同样顾及非局部特性，都是在小窗口中进行的，即计算均值。

步骤3：基于非局部滤波，采用以相似像元为中心的邻域矩阵的相似度以及 其与中心像元的相对距离计算各相似像元的权重；

基于非局部滤波(Non-local means)，用基于矩阵的方法计算各相似像元的权 重，考虑到二阶范式的计算可能会造成过渡平滑，使融合结果损失过多的细节信 息。因此，本发明采用用一阶形式进行空间滤波，并在此基础上加入相对距离， 具体实现方式为:

$W=\exp (-\frac{\bar{S}+\bar{T}}{h})*1/D---\left(2\right);$

其中，和分别为邻域矩阵的光谱差异和时间差异的均值，通过调节滤波 参数h和矩阵窗口来进行试验，最优的h的范围为[0.001，0.1]；D为邻近像元(x_i， y_j)与中心像元的相对距离。

步骤4：考虑时间变化和空间变化，对时间变化的区域采用经典算法方式融 合，对空间变化的区域采用时间加权的方式融合，最终融合得到中心像元的反射 率值。

其中对时间变化的区域(表现在地表覆盖(物候)在某时间段发生较大变化) 采用经典算法方式融合，其公式为：

$F(x_i,y_j,t_k)=\underset{i=1}{\overset{w}{\Sigma }}\underset{j=1}{\overset{w}{\Sigma }}\underset{k=1}{\overset{n}{\Sigma }}{W}_{\mathrm{ijk}}\times (F(x_i,y_j,t_0)+C(x_i,y_j,t_k)-C(x_i,y_j,t_0))---\left(3\right);$

其中，F、C分别代表着高空间分辨率与高时间分辨率数据的反射率， (x_i,y_j)为给定高空间分辨率与高时间分辨率数据像元的位置，t₀为数据获取的 时间，t_k为预测时间，w为搜索窗口大小，n为输入的影像对数；

对空间变化的区域(表现在地表覆盖(物候)随时间基本不发生变化，或变 化很小，一般发生在空间异质区域)采用时间加权的方式融合，其公式为：

F(x_w/2,y_w/2,t_p)＝F(x,y,t₀)+∑w(x,y)*a_w/2*(C(x,y,t_p)-C(x,y,t₀))    (4)； 其特点是计算移动窗口内相似像元的线性回归系数a_w/2，进而计算各基础时间的 时相权重，加权得到最终的融合结果。

应当理解的是，本说明书未详细阐述的部分均属于现有技术。

应当理解的是，上述针对较佳实施例的描述较为详细，并不能因此而认为是 对本发明专利保护范围的限制，本领域的普通技术人员在本发明的启示下，在不 脱离本发明权利要求所保护的范围情况下，还可以做出替换或变形，均落入本发 明的保护范围之内，本发明的请求保护范围应以所附权利要求为准。