CN111369483A - 一种融合多源遥感数据生成高时空分辨率遥感数据的方法 - Google Patents

Abstract

本发明提供一种融合多源遥感数据生成高时空分辨率遥感数据的方法，包括：对输入的不同时空分辨率遥感数据进行预处理；对高时间、低空间分辨率遥感数据，提取像元的变化趋势信息；分别在高时间分辨率和高空间分辨率遥感数据的像元上减去趋势值，使残差值满足一阶和二阶平稳性条件；根据高时间、低空间分辨率像元的残差值构建时间协方差模型；基于该时间协方差模型，对不同时间上观测的高空间、低时间分辨率的遥感数据的残差项做时间Kriging插值，估计其任意未观测时刻的残差值；将得到的趋势值与的残差值相加，得到任意高空间分辨率遥感未观测时刻的像元值。本发明的计算效率、数据利用率和精确度都明显提升，可灵活用于不同场景。

Description

一种融合多源遥感数据生成高时空分辨率遥感数据的方法

技术领域

本发明涉及遥感影像融合算法技术领域，特别涉及一种融合多源遥感数据生成高时空分辨率遥感数据的方法。

背景技术

近年来，遥感应用如土地覆盖变化、农情监测、资源环境调查等等领域的不断发展，推动了遥感数据的需求，遥感数据的空间分辨率、光谱分辨率、时间分辨率不断提高。但是单一的传感器获取的影像，仍然难以满足应用的需求。一方面，由于遥感传感器的硬件技术瓶颈，难以直接获得具有高空间和高时间分辨率的遥感影像；另一方面，由于云污染等影响，获得的影像往往存在像元值的缺失。在此背景下，遥感影像融合技术是解决该问题的有效方法。

随着遥感数据融合相关研究的增多，涌现出了大量的融合算法。其中应用最为广泛且的算法为STARFM以及ESTARFM算法。ESTARFM算法相较于STARFM算法，在异质性较强的地区更为适用。

然而，ESTARFM算法具有以下缺陷：首先，需要通过准确的寻找相似像元来保证融合的精度，在变化比较复杂的地区往往难以找到对的相似像元导致精度的下降；其次，它通过求两个不同时刻相似像元的交集来提升找到的相似像元的精度，但是在多云的地区，这种方式反而造成融合结果中像元的缺失；最后，寻找相似像元是个耗时的过程。总之，目前尚没有算法能够构建规则且完整的遥感影像时间序列，这也是本发明要解决的问题。

发明内容

本发明的目的在于提供一种融合多源遥感数据生成高时空分辨率遥感数据的方法，实现构建规则且完整的高空间、低时间分辨率遥感影像时间序列，提高数据利用率、计算效率和精确度。

为解决上述技术问题，本发明的实施例提供如下方案：

一种融合多源遥感数据生成高时空分辨率遥感数据的方法，包括以下步骤：

S1、对输入的高时间、低空间分辨率遥感影像和高空间、低时间分辨率遥感影像进行预处理，包括：使两者投影信息以及观测范围相同；根据数据的质量控制文件判断数据的质量情况，若存在被质量控制文件标识为低质量的像元，则除去该像元值将其视为未观测的像元；在高时间、低空间分辨率遥感影像中，若存在被除去像元值的低质量的像元，用多年观测数据的均值作为背景值，填充被除去的像元值，并对每个像元在时间上做线性插值，插值出背景值不能填充的像元值；

S2、对高时间、低空间分辨率遥感数据，提取像元的变化趋势信息，通过对得到的趋势值做重采样，得到高空间分辨率尺度下像元的趋势值；

S3、在高时间、低空间分辨率遥感影像上除去像元的趋势值，得到满足一阶和二阶平稳性条件的残差值，并根据高时间、低空间分辨率遥感影像上像元的残差值构建时间协方差模型，同样在高空间、低时间分辨率遥感影像上除去像元的趋势值得到残差值；

S4、利用构建的时间协方差模型，通过Kriging插值求得任意未观测时刻的高空间、低时间分辨率遥感影像上像元的残差值；

S5、对于高空间、低时间分辨率遥感影像，根据步骤S2中得到的趋势值以及步骤S4中得到的残差值，两者求和得出任意未观测时刻的像元值，逐像元预测直到得到整个高空间、低时间分辨率遥感影像的像元值，完成影像预测。

优选地，所述步骤S1还包括：

对线性插值后的高时间、低空间分辨率遥感影像的像元时间序列进行滤波处理，滤除噪声。

优选地，所述步骤S2包括：

S201、根据数据特点，选择相应的趋势提取方式，提取高时间、低空间分辨率遥感影像上像元的趋势值，所述趋势提取方式包括移动窗口法、模型拟合法；

S202、将得到的高时间、低空间分辨率遥感影像像元的趋势值重采样到高空间、低时间分辨率遥感影像尺度，得到高空间、低时间分辨率遥感影像像元的趋势值。

优选地，所述步骤S3包括：

S301、分别在高时间、低空间分辨率遥感影像和高空间、低时间分辨率遥感影像的像元上除去对应的趋势值C_{μ(x，y，t)}与F_{μ(x，y，t)}，公式如下：

C_{δ(x，y，t)}＝C_Z(x，y，t)-C_{μ(x，y，t)}

F_{δ(x，y，t)}＝F_Z(x，y，t)-F_{μ(x，y，t)}

其中，C_Z(x，y，t)与F_Z(x，y，t)分别表示高时间、低空间分辨率遥感影像和高空间、低时间分辨率遥感影像的像元值，C_{δ(x，y，t)}与F_{δ(x，y，t)}分别表示高时间、低空间分辨率遥感影像和高空间、低时间分辨率遥感影像像元的残差值；

S302、根据得到高时间、低空间分辨率遥感影像像元的残差值C_{δ(x，y，t)}构建时间协方差模型，选用的模型包括kriging适用的模型，根据拟合情况进行选择。

优选地，所述步骤S4包括：

通过Kriging插值求得任意未观测时刻t_p的残差值

公式如下：

其中，λ_i表示插值的权重，通过以下公式求得：

其中，

表示一个高时间、低空间分辨率遥感影像的像元在t_i与t_j时刻残差值的相关性，由时间协方差模型得到，

表示拉格朗日乘子。

优选地，所述步骤S5包括：

对于高空间、低时间分辨率遥感影像，通过以下公式得到任意未观测时刻t_p的像元值

其中，

为任意未观测时刻t_p高空间、低时间分辨率遥感影像像元的趋势值，

为t_p时刻高空间、低时间分辨率遥感影像像元的残差值。

本发明的上述方案至少包括以下有益效果：

1)不用寻找相似像元，主要运算仅一个插值过程，计算效率明显提升；

2)能够使用时间变程内所有有效的像元值，数据利用率大大提升；

3)避免寻找相似像元带来的误差，精度提高；

4)不需要过多调参，可灵活应用于不同场景。

附图说明

图1是本发明实施例提供的融合多源遥感数据生成高时空分辨率遥感数据的方法的流程图；

图2是本发明实施例中融合多源遥感数据生成高时空分辨率遥感数据的方法的具体实施过程图。

具体实施方式

为使本发明要解决的技术问题、技术方案和优点更加清楚，下面将结合附图及具体实施例进行详细描述。

本发明的实施例提供了一种融合多源遥感数据生成高时空分辨率遥感数据的方法，如图1所示，该方法包括以下步骤：

S1、对输入的高时间、低空间分辨率遥感影像和高空间、低时间分辨率遥感影像进行预处理，包括：使两者投影信息以及观测范围相同；根据数据的质量控制文件判断数据的质量情况，若存在被质量控制文件标识为低质量的像元，则除去该像元值将其视为未观测的像元；在高时间、低空间分辨率遥感影像中，若存在被除去像元值的低质量的像元，用多年观测数据的均值作为背景值，填充被除去的像元值，并对每个像元在时间上做线性插值，插值出背景值不能填充的像元值；

S2、对高时间、低空间分辨率遥感数据，提取像元的变化趋势信息，通过对得到的趋势值做重采样，得到高空间分辨率尺度下像元的趋势值；

S3、在高时间、低空间分辨率遥感影像上除去像元的趋势值，得到满足一阶和二阶平稳性条件的残差值，并根据高时间、低空间分辨率遥感影像上像元的残差值构建时间协方差模型，同样在高空间、低时间分辨率遥感影像上除去像元的趋势值得到残差值；

S4、利用构建的时间协方差模型，通过Kriging插值求得任意未观测时刻的高空间、低时间分辨率遥感影像上像元的残差值；

S5、对于高空间、低时间分辨率遥感影像，根据步骤S2中得到的趋势值以及步骤S4中得到的残差值，两者求和得出任意未观测时刻的像元值，逐像元预测直到得到整个高空间、低时间分辨率遥感影像的像元值，完成影像预测。

本发明方法对输入的观测高时间、低空间分辨率遥感影像以及高空间、低时间分辨率遥感影像进行预处理；在具有高频次观测的高时间、低空间分辨率遥感影像上提取像元的趋势信息；分别在观测的高时间、低空间分辨率遥感影像以及高空间、低时间分辨率遥感影像的像元上除去得到的趋势值，得到满足平稳性假设的残差值，并根据高频次的高时间、低空间分辨率像元的残差值构建时间协方差模型，得到残差值在不同时间间隔上的相关性；利用该相关性对高空间、低时间分辨率遥感影像的残差值做时间上的Kriging插值，求得未观测时刻的高空间、低时间分辨率残差值；在求得的高空间、低时间分辨率残差值上返回趋势信息，得到未观测时刻的像元值，完成影像预测。

本发明方法具有以下优点：无需寻找相似像元，主要运算仅一个插值过程，计算效率明显提升；能够使用时间变程内所有有效的像元值，数据利用率大大提升；有效避免了寻找相似像元带来的误差，精度提高；不需要过多调参，可灵活应用于不同场景。

进一步地，步骤S1还包括：

根据数据特点，对线性插值后的高时间、低空间分辨率遥感影像的像元时间序列可进行适当的滤波处理，滤除噪声。

该步骤在保证高时间、低空间分辨率遥感影像时间序列是完整的前提下，避免噪声干扰，为之后的趋势计算提供充分的信息。

进一步地，步骤S2包括：

S201、根据数据特点，选择相应的趋势提取方式，提取高时间、低空间分辨率遥感影像上像元的趋势值，所述趋势提取方式包括移动窗口法、模型拟合法；

以移动窗口法为例，通过选用时空移动窗口提取每个窗口内高时间、低空间分辨率遥感影像上像元值的中位数作为像元的趋势值，计算公式如下：

其中，C_{μ(x，y，t)}表示高时间、低空间分辨率遥感影像像元的趋势值，x、y和t分别表示像元所在的空间位置以及时刻，a、b和c表示移动窗口的大小为a*b*c。

S202、将得到的高时间、低空间分辨率遥感影像像元的趋势值重采样到高空间、低时间分辨率遥感影像尺度，得到高空间、低时间分辨率遥感影像像元的趋势值F_{μ(x，y，t)}。

进一步地，步骤S3包括：

S301、分别在高时间、低空间分辨率遥感影像和高空间、低时间分辨率遥感影像的像元上除去对应的趋势值C_{μ(x，y，t)}与F_{μ(x，y，t)}，公式如下：

C_{δ(x，y，t)}＝C_Z(x，y，t)-C_{μ(x，y，t)}

F_{δ(x，y，t)}＝F_Z(x，y，t)-F_{μ(x，y，t)}

其中，C_Z(x，y，t)与F_Z(x，y，t)分别表示高时间、低空间分辨率遥感影像和高空间、低时间分辨率遥感影像的像元值，C_{δ(x，y，t)}与F_{δ(x，y，t)}分别表示高时间、低空间分辨率遥感影像和高空间、低时间分辨率遥感影像像元的残差值；

S302、根据得到高时间、低空间分辨率遥感影像像元的残差值C_{δ(x，y，t)}构建时间协方差模型，选用的模型包括kriging适用的模型，根据拟合情况进行选择。

作为本发明的一种实施方式，选用的模型包括kriging插值适用的模型。

以球面模型为例，时间协方差模型的公式为：

其中，c₁表示时间协方差模型的基台值，D表示像元间的时间间隔，a表示时间协方差模型的变程，c(t)表示协方差值。

进一步地，步骤S4包括：

通过Kriging插值求得任意未观测时刻t_p的残差值

公式如下：

其中，λ_i表示插值的权重，通过以下公式求得：

其中，

表示一个高时间、低空间分辨率遥感影像的像元在t_i与t_j时刻残差值的相关性，由时间协方差模型得到，

表示拉格朗日乘子。

进一步地，步骤S5包括：

对于高空间、低时间分辨率遥感影像，通过以下公式得到任意未观测时刻t_p的像元值

其中，

为任意未观测时刻t_p高空间、低时间分辨率遥感影像像元的趋势值，

为t_p时刻高空间、低时间分辨率遥感影像像元的残差值。

本发明方法的具体应用过程如图2所示。

(1)预处理

输入数据为高时间、低空间分辨率遥感影像以及高空间、低时间分辨率遥感影像时间序列，经过投影转换、影像配准、影像裁剪等处理，使高空间、低时间分辨率遥感影像以及高时间、低空间分辨率遥感影像的投影信息与观测范围相同。

根据数据的质量控制文件判断数据的质量情况，若存在被质量控制文件标识为低质量的像元，则除去该像元值将其视为未观测的像元，避免云污染等造成的干扰。

之后，在高时间、低空间分辨率遥感影像中，用多年观测数据的均值作为背景值，填充被除去的低质量的像元值，然后对每个像元在时间上做线性插值，插值出背景值不能填上的像元值，保证时间序列的完整。

根据数据的特点，这里可做适当的滤波处理，使得时间序列变化更加平滑，避免噪声干扰。

本步骤的目的是使高时间、低空间分辨率遥感影像整个时间序列是完整的，排除云污染以及噪声干扰，为之后的趋势计算提供充分的信息。

(2)趋势计算

根据数据特点，选择适当的趋势提取方式，提取高时间、低空间分辨率遥感影像上像元的趋势值，可使用移动窗口法、模型拟合法等；

以移动窗口法为例，通过选用时空移动窗口提取每个窗口内高时间、低空间分辨率遥感影像上像元值的中位数作为像元的趋势值，计算公式如下：

其中，C_{μ(x，y，t)}表示高时间、低空间分辨率遥感影像像元的趋势值，x、y和t分别表示像元所在的空间位置以及时刻，a、b和c表示移动窗口的大小为a*b*c。

将得到的高时间、低空间分辨率遥感影像像元的趋势值重采样到高空间、低时间分辨率遥感影像尺度，得到高空间、低时间分辨率遥感影像像元的趋势值F_{μ(x，y，t)}。

(3)时间协方差模型建立

分别在高时间、低空间以及高空间、低时间分辨率遥感影像像元上除去步骤(2)中得到的对应趋势C_{μ(x，y，t)}与F_{μ(x，y，t)}，公式如下：

C_{δ(x，y，t)}＝C_Z(x，y，t)-G_{μ(x，y，t)}

F_{δ(x，y，t)}＝F_Z(x，y，t)-F_{μ(x，y，t)}

式中C_Z(x，y，t)与F_Z(x，y，t)分别表示高时间、低空间以及高空间、低时间分辨率遥感影像像元值，C_{δ(x，y，，t)}与F_{δ(x，y，t)}分别表示高时间、低空间以及高空间、低时间分辨率遥感影像像元的残差值。根据地统计学可知，此时的残差值是满足一阶与二阶平稳性条件的，这里能够得到满足平稳性假设的残差值C_{δ(x，y，t)}与F_{δ(x，y，t)}。

根据得到的C_{δ(x，y，t)}构建时间协方差模型，该时间协方差模型反应了同一位置的像元在不同时间间隔下相关性强弱，通常时刻越接近相关性越高。选用的模型可以是任何kriging适用的模型，根据拟合情况灵活选择。

(4)Kriging插值

在步骤(3)中求得的残差值F_{δ(x，y，t)}满足平稳性的情况下，可以通过Kriging插值求得高空间、低时间分辨率影像任意未观测时刻t_p的残差值

公式如下：

式中的λ_i表示插值的权重，在Kriging插值中它可以通过以下公式求得：

式中的

表示一个高时间、低空间分辨率像元在t_i与t_j时刻残差值的相关性，由构建的时间协方差模型得到，

表示拉格朗日乘子。

通过Kriging插值可得出所需任意未观测时刻的高空间、低时间分辨率像元残差值。

(5)影像预测

对于任意未观测时刻t_p，根据步骤(2)中得到的趋势

以及步骤(4)中得到的残差

通过以下公式可以得到该时刻的高空间、低时间分辨率像元值

逐像元预测进而得到整个高时间、低空间分辨率遥感影像的像元值，完成影像的预测。

与现有技术相比，本发明方法不同于传统融合算法的基于增量再分配的方式，而是将像元值当作随机变量，把地统计学引入像元值预测当中，将像元值的预测转换为一个插值过程。

本发明方法中的Kriging插值不同于传统的空间与时空Kriging插值，提出了仅在时间上进行的Kriging插值方式，避免了引入空间维度带来的平滑效应。

本发明方法用高时间、低空间分辨率遥感影像提供趋势变化信息并构建时间协方差模型，为预测影像提供了像元的趋势和时间相关性，用高空间、低时间分辨率遥感影像插值出残差项，为预测影像提供了空间细节的信息。

本发明的计算效率、数据利用率和精确度都明显提升，可灵活应用于不同场景。

以上所述是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明所述原理的前提下，还可以做出若干改进和润饰，这些改进和润饰也应视为本发明的保护范围。

Claims (6)

1.一种融合多源遥感数据生成高时空分辨率遥感数据的方法，其特征在于，包括以下步骤：

S1、对输入的高时间、低空间分辨率遥感影像和高空间、低时间分辨率遥感影像进行预处理，包括：使两者投影信息以及观测范围相同；根据数据的质量控制文件判断数据的质量情况，若存在被质量控制文件标识为低质量的像元，则除去该像元值将其视为未观测的像元；在高时间、低空间分辨率遥感影像中，若存在被除去像元值的低质量的像元，用多年观测数据的均值作为背景值，填充被除去的像元值，并对每个像元在时间上做线性插值，插值出背景值不能填充的像元值；

S2、对高时间、低空间分辨率遥感数据，提取像元的变化趋势信息，通过对得到的趋势值做重采样，得到高空间分辨率尺度下像元的趋势值；

S3、在高时间、低空间分辨率遥感影像上除去像元的趋势值，得到满足一阶和二阶平稳性条件的残差值，并根据高时间、低空间分辨率遥感影像上像元的残差值构建时间协方差模型，同样在高空间、低时间分辨率遥感影像上除去像元的趋势值得到残差值；

S4、利用构建的时间协方差模型，通过Kriging插值求得任意未观测时刻的高空间、低时间分辨率遥感影像上像元的残差值；

S5、对于高空间、低时间分辨率遥感影像，根据步骤S2中得到的趋势值以及步骤S4中得到的残差值，两者求和得出任意未观测时刻的像元值，逐像元预测直到得到整个高空间、低时间分辨率遥感影像的像元值，完成影像预测。

2.根据权利要求1所述的融合多源遥感数据生成高时空分辨率遥感数据的方法，其特征在于，所述步骤S1还包括：

对线性插值后的高时间、低空间分辨率遥感影像的像元时间序列进行滤波处理，滤除噪声。

3.根据权利要求1所述的融合多源遥感数据生成高时空分辨率遥感数据的方法，其特征在于，所述步骤S2包括：

S201、根据数据特点，选择相应的趋势提取方式，提取高时间、低空间分辨率遥感影像上像元的趋势值，所述趋势提取方式包括移动窗口法、模型拟合法；

S202、将得到的高时间、低空间分辨率遥感影像像元的趋势值重采样到高空间、低时间分辨率遥感影像尺度，得到高空间、低时间分辨率遥感影像像元的趋势值。

4.根据权利要求3所述的融合多源遥感数据生成高时空分辨率遥感数据的方法，其特征在于，所述步骤S3包括：

S301、分别在高时间、低空间分辨率遥感影像和高空间、低时间分辨率遥感影像的像元上除去对应的趋势值C_{μ(x，y，t)}与F_{μ(x，y，t)}，公式如下：

其中，

与

分别表示高时间、低空间分辨率遥感影像和高空间、低时间分辨率遥感影像的像元值，

与

分别表示高时间、低空间分辨率遥感影像和高空间、低时间分辨率遥感影像像元的残差值；

S302、根据得到高时间、低空间分辨率遥感影像像元的残差值

构建时间协方差模型，选用的模型包括kriging插值适用的模型，根据拟合情况进行选择。

5.根据权利要求4所述的融合多源遥感数据生成高时空分辨率遥感数据的方法，其特征在于，所述步骤S4包括：

通过Kriging插值求得任意未观测时刻t_p的残差值

公式如下：

其中，λ_i表示插值的权重，通过以下公式求得：

其中，

表示一个高时间、低空间分辨率遥感影像的像元在t_i与t_j时刻残差值的相关性，由时间协方差模型得到，

表示拉格朗日乘子。

6.根据权利要求5所述的融合多源遥感数据生成高时空分辨率遥感数据的方法，其特征在于，所述步骤S5包括：

对于高空间、低时间分辨率遥感影像，通过以下公式得到任意未观测时刻t_p的像元值

其中，

为任意未观测时刻t_p高空间、低时间分辨率遥感影像像元的趋势值，

为t_p时刻高空间、低时间分辨率遥感影像像元的残差值。

Images