一种长时序遥感土壤水分降尺度方法与流程

本发明涉及一种土壤水分降尺度方法，特别是一种耦合地表参数辅助转换及不同尺度间统计关系的长时序降尺度方法，属于遥感技术领域。

背景技术：

土壤水分是地球系统重要的状态变量，其时空变化对于地—气系统能量与水气交换有重要作用，对水循环过程及农业生产等均会造成相应影响，是水文生态领域的重要研究对象。近年来，遥感(尤其是微波遥感)以其空间连续等特点在流域土壤水分监测中得到广泛应用。然而，由于遥感技术发展历史较短，加之不同平台的寿命有限，现有遥感产品的时间范围较短。此外，微波遥感土壤水分空间分辨率较低，难以反映土壤水分的空间细节特征。降尺度方法是获取高分辨率土壤水分数据集的有效途径，不少学者此方面开展了广泛研究，从方法上大致可分为地表参数辅助尺度转换和基于不同尺度间土壤水分统计关系两种。

前者主要利用土壤水分与其它高空间分辨率地表参数(如地表覆盖、地表温度等)之间的物理或经验性统计关系进行尺度转换。

第二种方法主要基于土壤水分“时间稳定性(Temporal stability)”原理，即不同尺度间土壤水分在长时间序列上具有一致的变化趋势，通过不同尺度间土壤水分的相互关系，将低尺度数据转换至较高的空间分辨率。

受观测条件与技术等因素的限制，地面辅助参数和高分辨率土壤水分数据产品通常存在数据缺失等情况，导致上述方法在构建长时序土壤水分数据集时依然面临许多挑战，成为流域土壤水分长期精准监测亟待解决的问题。

技术实现要素：

针对现有技术存在的上述不足，本发明的目的是提供一种耦合地表参数辅助转换及不同尺度间统计关系的长时序降尺度方法。

为实现上述目的，本发明采用如下技术手段：

一种长时序遥感土壤水分降尺度方法，其特征在于，它包括以下步骤：

(一)根据MODIS数据的获取情况对CCI长时序遥感土壤水分数据产品进行划分：

将CCI长时序遥感土壤水分数据产品划分为第一部分CCI数据产品和第二部分CCI数据产品；所述第一部分CCI数据产品为MODIS数据可获取时段内的CCI长时序遥感土壤水分数据产品；所述第二部分CCI数据产品为MODIS数据缺失时段内的CCI长时序遥感土壤水分数据产品；

(二)由MODIS数据得出2001年至2013年内分辨率为1Km的土壤水分数据：

(1)接收2001年至2013年内分辨率为1Km的MODIS数据，所述MODIS数据包括：植被指数数据NDVI和地表温度数据T；

(2)对植被指数所述数据NDVI和地表温度数据T进行归一化处理，消除数量级别上的差异性：

T^*＝(T-T_min)/(T_max-T_min) (1)

NDVI^*＝(NDVI-NDVI_min)/(NDVI_max-NDVI_min) (2)

其中，T_max为地表温度数据T的最大值，T_min为地表温度数据T的最小值；NDVI_max为归一化植被指数数据NDVI的最大值，NDVI_min为归一化植被指数数据NDVI的最小值，T^*为归一化后地表温度数据，NDVI^*为归一化后的植被指数数据；

(3)对归一化后的地表温度数据T^*和归一化后的植被指数数据NDVI^*进行重采样至分辨率为25Km；

基于重采样得到的2001年至2013年内的归一化后的地表温度数据T^*和归一化后的植被指数数据NDVI^*建立第一阶段CCI土壤水分数据的统计关系：

$SSM=\underset{i=0}{\overset{n}{\Σ}}\underset{j=0}{\overset{n}{\Σ}}{a}_{ij}\·{\mathrm{NDVI}}^{*\left(i\right)}\·T^{*\left(j\right)}---\left(3\right)$

其中a_ij表示权重系数，i,j分别代表自变量的维数，SSM为2001年至2013年内的土壤水分数据；

(4)将步骤(1)中获取的分辨率为1Km的植被指数数据NDVI和分辨率为1Km的地表温度数据T数据代入公式(1)、(2)和(3)中，实现土壤水分的降尺度处理，得出2001年至2013年内的分辨率为1Km的土壤水分数据；

(三)以步骤(二)获取的数据为参考数据，获取土壤水分尺度转换参数，进而将步骤(一)划分得到的第二部分CCI数据产品转换为分辨率为1Km的土壤水分数据：

(1)确定拟合转换系数：

不同尺度间土壤水分的关系为

${\mathrm{SSM}}_{x,y,t}^F=a_{x,y}\·{\mathrm{SSM}}_{x,y,t}^C+b_{x,y}---\left(4\right)$

其中，和分别表示为t时刻(x,y)位置上高分辨率及其对应低分辨率像元内的土壤水分；a_x,y和b_x,y表示拟合转换系数；

将步骤(二)得到的2001年至2013年内的分辨率为1Km的土壤水分数据代入式(4)中的将步骤(一)划分得到的第一部分CCI数据代入式(4)中的得出拟合转换系数a_x,y和b_x,y的值；

(2)将步骤(一)划分得到的第二部分CCI数据代入入式(4)中的将步骤(1)中得到的拟合转换系数a_x,y和b_x,y的值分别代入式(4)中的a_x,y和b_x,y，得到1979年至2000年内的分辨率为1Km的土壤水分数据。

进一步的，所述MODIS数据可获取时段为2001年至2013年。

相比现有技术，本发明具有如下有益效果：

受观测条件与技术等因素的限制，地面辅助参数和高分辨率土壤水分数据产品通常存在数据缺失等情况，导致现有技术中通常使用的两种方法在构建长时序土壤水分数据集时依然面临许多挑战。本发明综合已有方法的特点及优劣性，提出了耦合地表参数辅助转换及不同尺度间统计关系的长时序遥感土壤水分降尺度方法，相比现有技术，本发明具有以下优势：(1)本发明可有效解决土壤水分降尺度对高分辨率数据的依赖性，使在有限的数据支持下获取长时序高分辨率土壤水分成为现实；(2)本发明算法结构简单，具有很好的可操作性，避免了物理模型结构复杂、输入参数多等问题，具有较强的实用性与推广性。

附图说明

图1为本发明的流程示意图；

图2为基于地表参数辅助尺度转换方法的土壤水分降尺度结果示意图；

图3为基于本算法的1979-2013年多年平均土壤水分空间分布图

图4为降尺度土壤水分空间特征验证区分布图

图5为TVDI与原始CCI数据、降尺度土壤水分的相关性示意图；

图6为降尺度校准与检验结果示意图。

具体实施方式

下面结合附图和实施例，对本发明做进一步详细说明。

一、本发明综合已有方法的特点及优劣性，以欧空局发表的CCI长时序(1979-2013年，25Km)土壤水分为例，公开了耦合地表参数辅助转换及不同尺度间统计关系的长时序降尺度方法。

如图1所示，本发明主要包括三大步骤：(一)根据高分辨率地表参数的获取情况，将长时序低分辨率遥感数据划分为两个阶段；(二)根据土壤水分与地表参数之间的统计关系，由高分辨率地表参数数据得出其对应年份的高分辨率土壤水分数据；(三)以步骤(二)获取的数据为参考数据，获取土壤水分尺度转换参数，进而将步骤(一)划分得到的无高分辨率地表参数时期的土壤水分转换高分辨率数据。

以下对每一步骤进行详细阐述。

(一)对由欧洲空间局发表的CCI长时序遥感土壤水分数据产品进行划分。

基本思路是：根据高分辨率地表参数数据(以2001-2013年MODIS数据为例，空间分辨率1km)的获取情况，对长时序低分辨率遥感土壤水分数据产品(以欧空局1979-2013年CCI土壤水分为例，空间分辨率25km)进行划分：第一部分对应于地表参数可获取时段(2001-2013年)，第二部分为地表参数缺失时段。

MODIS自2000年4月开始正式发布数据，我国目前已建立了数个接收站并分别于2001年3月前后开始接收数据。因此，由于高分辨率MODIS数据的获取时间始于2001年，因而基于地表参数辅助尺度转换方法无法应用于在此之前的土壤水分降尺度处理。

(二)由MODIS数据得出2001年至2013年内分辨率为1Km的土壤水分数据：

基本思路在于：本领域公知的，土壤水分与植被覆盖、地表温度等因素密切相关。基于该原理，本发明采用分辨率为1Km的MODIS归一化植被指数(Normalized Difference Vegetation Index，NDVI)及地表温度(Land surface temperature，LST)数据产品，建立与土壤水分的经验方程，实现土壤水分的降尺度处理。

基本过程首先将分辨率为1Km的NDVI数据及LST数据进行归一化处理，消除数量级别上的差异性，公式为：

T^*＝(T-T_min)/(T_max-T_min) (1)

NDVI^*＝(NDVI-NDVI_min)/(NDVI_max-NDVI_min) (2)

其中，T_max，T_min分别为最大与最小地表温度，NDVI_max，NDVI_min分别为最大与最小NDVI值，T^*和NDVI^*为归一化后数据。

继而将归一化NDVI与LST数据重采样至分辨率为25Km(在遥感中，重采样是从高分辨率遥感影像中提取出低分辨率影像的过程)，并与相应的CCI数据建立统计关系：

$SSM=\underset{i=0}{\overset{n}{\Σ}}\underset{j=0}{\overset{n}{\Σ}}{a}_{ij}\·{\mathrm{NDVI}}^{*\left(i\right)}\·T^{*\left(j\right)}---\left(3\right)$

其中a_ij表示权重系数，i,j分别代表自变量的维数，二次或三次回归方程即可取得较好的结果。由于SSM，NDVI和T三个变量均已知，通过反算即可得到权重系数a_ij。最终，由于公式(3)中的SSM，NDVI和T等参数已知，通过反算即可得到系数a_ij，进而可利用该系数将原始的1Km分辨率的NDVI及LST数据代入公式(3)，实现土壤水分的降尺度处理，得出2001年至2013年内的分辨率为1Km的土壤水分数据。

(三)将步骤(一)划分得到的第二部分CCI数据产品转换为1979年至2000年内的分辨率为1Km的土壤水分数据：

基本思路在于：将步骤(二)得到的分辨率为1Km降尺度土壤水分数据(2001年-2013年内的数据)与对应的原始25Km CCI土壤水数据(2001年-2013年内的数据)分代入公式(4)，利用土壤水分不同尺度间统计关系方法，获取拟合转换系数a_x,y和b_x,y的值。

不同尺度间土壤水分的关系为

${\mathrm{SSM}}_{x,y,t}^F=a_{x,y}\·{\mathrm{SSM}}_{x,y,t}^C+b_{x,y}---\left(4\right)$

其中，和分别表示为t时刻(x,y)位置上高分辨率及其对应低分辨率像元内的土壤水分；a_x,y和b_x,y表示转换系数。

为消除尺度转换可能导致像元边缘出现不连续的“马赛克”效应，采用地理加权回归(Geographically Weighted Regression，GWR)方法首先将25Km分辨率的CCI土壤水分数据重采样至1Km。由于土壤水分在空间上具有渐变特征，GWR方法基于地理学第一定律，充分考虑了不同空间范围内土壤水分之间的相关性，将某像元土壤水分通过其自身及相邻像元值加权平均得到，权1重系数采用倒距离函数方法确定，从而避免了土壤水分在空间上的突变现象。

在此基础上，利用公式(4)将其它年份(地表参数缺失时段，2001年之前)CCI数据转换为空间分辨率为1Km的数据，也即是将步骤(一)划分得到的第二部分CCI数据(分辨率为25Km)代入入式(4)中的将步骤(1)中得到的拟合转换系数a_x,y和b_x,y的值分别代入式(4)中的a_x,y和b_x,y，得到地表参数(MODIS数据)缺失时段内的分辨率为1Km的土壤水分数据。

这里所依据的基本原理是：土壤水分具有时间稳定性(Temporal stability)特征，即局部土壤水分长时序变化与区域总体趋势具有良好的一致性。简单而言，受大范围土壤水分的总体影响，某空间小区域内土壤水分的时间变化总体上与大范围内的土壤水分保持一致。这为土壤水分的尺度转换奠定了基础，根据这一特征，不同尺度间土壤水分的关系在式(4)中已经给出。

此外，需要指出的是：本发明中MODIS数据可获取时段为2001年至2013年，CCI长时序遥感土壤水分数据可获取的时段为1979年至2013年，这是由本领域当前技术现状决定的，但是在本领域技术水平发展到能够获取上述两个时段以外的相应数据的情况下，本发明提供的方法也能够解决本发明所提出的所要解决的技术问题。

二、对上述方法所得数据进行校准与检验的方法

本发明最后从精度和空间准确性两个方面评价降尺度结果的有效性。在精度检验方面，本方法根据土壤水分的时间稳定性特征(Temporal stability)，采用长时序地面观测点的土壤水分数据，通过回归分析和线性拉伸两种方法校准遥感土壤水分。回归分析方法描述了地面观测点与遥感土壤水分数据之间的线性相关关系，表达式为：

SM_cal＝a×SM_RS+b (5)

式中，SM_cal、SM_RS分别表示校准及原始遥感土壤水分数据；a和b为待定参数，一般根据长时序的地面观测点和对应的遥感土壤水分数据，采用最小二乘法进行估计，使两组数据间的误差最小。

线性拉伸方法则考虑地面观测数据变化特征的实际情况，对遥感土壤水分数据进行拉伸处理，使其与地面观测数据具有相同的均值与方差，其表达式为：

${\mathrm{SM}}_{cal}=({\mathrm{SM}}_{RS}-{\stackrel{\‾}{SM}}_{RS})\×(s_{in}/s_{RS})+{\stackrel{\‾}{SM}}_{in}---\left(6\right)$

式中，和分别为地面观测点和对应遥感长时序土壤水分平均值；s_in和s_RS分别两种数据的标准差。

为验证降尺度土壤水分数据空间特征的准确性，本发明采用TVDI指数方法验证与土壤水分的相关性。本算法利用1Km MODIS地表温度和NDVI数据，构建温度—植被指数空间，验证降尺度数据空间特征的准确性。为避免不同数据量纲上的不一致性，两种数据分别归一化至0-1区间，代表土壤干湿状况的递变程度。采用决定系数(Determination coefficient)R²反映降尺度数据空间准确性的总体情况；同时考虑到气象等因素的地理差异性，选择典型剖面线作为验证区域，检验不同纬度、气象等条件下土壤水分降尺度结果。

三、应用实例

以鄱阳湖流域为例，验证了本算法的效果。图2显示了2001-2013年原始CCI数据与基于地表参数辅助尺度转换方法土壤水分降尺度结果对比图(左图为2001-2013年原始CCI多年平均数据，右图为降尺度后多年平均数据)，可以发现，土壤水分降尺度结果与原始CCI数据具有较强的空间一致性，表现为流域多年平均土壤水分自北向南呈下降趋势，降尺度结果更好地反映了土壤水分空间分布的细节特征。原始CCI数据流域多年平均值为0.350cm³/cm³，标准差为0.035cm³/cm³，降尺度土壤水分的流域多年平均和标准差分别为0.346cm³/cm³和0.008cm³/cm³，空间异质性较原始CCI数据有所降低。流域土壤水分最高出现于湖区，原始CCI数据和降尺度结果土壤水分分别为0.402cm³/cm³和0.351cm³/cm³。

根据所获取的尺度转换系数a_x,y和b_x,y，采用公式(4)获取1979-2013年1km土壤水分数据集，其多年平均空间分布如图3所示，流域多年平均为0.350cm³/cm³，标准差为0.010cm³/cm³。

利用TVDI数据，验证降尺度土壤水分空间分布特征的准确性。结果发现两种数据在流域尺度上具有较强的一致性，相关系数为0.45(R²＝0.20)。考虑到气象等因素的地理差异性，本发明分别选择两条跨越研究区东西方向和南北方向的剖面线作为验证区域，检验不同纬度与气象条件下土壤水分降尺度结果(如图4所示)。图5显示了剖面线方向上降尺度土壤水分和TVDI的分布结果，其中(a)部分为东西方向，(b)部分为南北方向。可以发现，相对于原始CCI数据，降尺度后土壤水分更好地吻合土壤水分的空间分布特征，相关系数R在0.3～0.4之间(R²>0.1)。

采用线性回归方法，利用地面观测数据校准降与检验尺度后土壤水分，结果如图6所示。统计结果表明，校准后降尺度土壤水分较好地吻合了地面观测数据的变化趋势，相关系数和RMSE分别为0.65(R²＝0.42)和0.044cm³/cm³，表明土壤水分降尺度与校准处理取得较好效果。

最后说明的是，以上实施例仅用以说明本发明的技术方案而非限制，尽管参照较佳实施例对本发明进行了详细说明，本领域的普通技术人员应当理解，可以对本发明的技术方案进行修改或者等同替换，而不脱离本发明技术方案的宗旨和范围，其均应涵盖在本发明的权利要求范围当中。