CN108154193B - 一种长时间序列降水数据降尺度方法 - Google Patents

Abstract

本发明公开了一种长时间序列降水数据降尺度方法，它包括以下步骤：(一)基于降水的空间差异特征构建水汽运移成本距离因子；(二)由MODIS EVI数据和步骤(一)得出的水汽运移成本距离因子得出2001年至2013年内分辨率为1km的降水数据。本发明具有如下有益效果：本发明可以有效解决卫星降水反演产品空间分辨率较低的问题，使在有限的数据支持下获取长时序高分辨率空间化降水成为现实；本发明算法结构简单，避免了物理模型结构复杂、输入参数多等问题。

Description

一种长时间序列降水数据降尺度方法

技术领域

本发明涉及一种降水数据降尺度方法，特别是一种耦合地表植被与地形信息等辅助参数的长时序降水数据降尺度方法，属于遥感技术领域。

背景技术

降水是全球水循环的重要组成部分，对地表物质循环和能量交换具有重要作用，降水是生态学、水文学和气象学的重要参数，特别是高分辨率的空间化降水资料对于区域水文水资源分析、旱涝灾害预报和生态环境治理有着重要的基础作用。近年来，遥感(尤其是微波遥感)以其空间连续性等特点在流域降水检测中得到广泛应用。然而，由于微波遥感技术发展历史较短，加之降水发生的高随机性与异质性，使得已有微波遥感降水数据分辨率较低，难以反映降水的空间细节特征。降尺度方法是获取高分辨率降水数据集的有效途径，不少学者在此方面开展了广泛研究，从方法上大致可分为动力降尺度方法和统计降尺度方法两种。

前者主要根据降水产生的物理环境进行模型模拟，通过将区域气候模式内嵌至全球气候模式中，或根据全球气候模式的输出作为边界条件进行阈值限定，生成高分辨率的区域尺度降水资料。由于模型复杂且通常难以获取模型所需的多种高空间分辨率的输入变量，因此其使用受到一定程度的限制。

第二种方法主要基于降水量与其它高空间分辨率地表参数(如海拔高度、植被覆盖程度等)之间的经验性统计关系进行尺度转换。它通过引入先验知识抓住降水时空分布规律性与随机性这一本质，方法灵活多样，简便易行，因此被广泛使用，例如申请号为201610305772.5的中国发明专利申请“一种基于数据挖掘的干旱监测方法”公开了如下技术方案：首先把1km的环境变量因子NDVI、DEM聚合到25km，作为自变量，对应25km分辨率的TRMM数据为因变量，建立地理加权回归模型，并将建立的模型应用到相应地理区域的1km环境变量因子上，最终得出1km的高精度降水预测数据。类似地，申请号为201611078012.1的发明专利申请中将所用到的环境变量增加到包括DNVI、地表温度、高程、坡度、坡向、经度、纬度在内的7个数据；申请号为201610307333.8的发明专利申请中使用的环境变量包括植被指数、数字高程模型、白天地表温度、晚上地表温度、地形湿度指数、坡度、坡向、坡长8个数据；申请号为201610305772.5的发明专利申请中使用的环境变量因子包括植被指数、数字高程模型、白天地表温、晚上地表温、地形湿度指数、坡度、地表粗糙度、地表反射率和谷底平坦指数9个数据。

但是，由于降水的高时空异质性，不同区域间影响降水因素复杂多样，各个区域降水模式间差异通常十分巨大，传统的环境因子(如NDVI与海陆位置)在模拟降水时存在一定缺陷，例如：NDVI在高植被覆盖区易于饱和；而海陆位置(通常用经度、纬度和目标区域与水汽来源地的欧几里得距离表示)往往不能有效模拟出地形起伏明显区域降水量的空间差异情况。同时，不断增加的预测因子一方面增加极大地增加数据准备的难度，比如地表温湿度和反射率数据很难保证数据的空间上的完整性和时间上的连续性；另一方面，多个环境变量间通常存在多重共线性，可能导致数据的冗余同时也增加了拟合模型的复杂性。

综上所述，现有技术在构建高分辨率长时序降水量数据集是仍然面临许多挑战，特别是提出有效的降尺度环境因子，成为流域降水的长期精准检测亟待解决的问题。

发明内容

针对现有技术存在的上述不足，本发明的目的是提供一种耦合降水与地表辅助参数间关系以实现流域降水长期精准检的长时序降水数据降尺度方法。

为了解决上述问题，本发明采用了以下的技术方案。

一种长时间序列降水数据降尺度方法，其特征在于，它包括以下步骤：

(一)基于降水的空间差异特征构建水汽运移成本距离因子：

子步骤1-1：获取90m空间分辨率的STRM DEM数据，并将数据重采样至空间分辨率为1km；

子步骤1-2：获取全球海岸线矢量数据，并将矢量数据转换为1km分辨率的栅格数据；

子步骤1-3：将研究区划分为1km大小的网格，以海岸线作为起点，以DEM数据为路径成本，计算来自海洋的水汽运移到研究区内各网格的成本距离；各网格成本距离计算方法为：

COST_m＝(DEM_i-DEM_j)/2 (1)

其中，DEM_i和DEM_j分别为两个相邻网格的海拔高程；COST_m、COST_n分别为以DEM数据为阻抗的邻接像元成本距离和对角像元成本距离；

研究区内各网格的水汽运移成本距离定义为从起始网格到各终止网格的累积成本距离；累积成本距离为起始像元移动到终止像元的成本距离之和；

(二)由MODIS EVI数据和步骤(一)得出的水汽运移成本距离因子得出2001年至2013年内分辨率为1km的降水数据：

子步骤2-1：接收2001年至2013年内分辨率为1km的MODIS增强型植被指数EVI数据；

子步骤2-2：对水汽运移成本距离和EVI进行归一化处理，消除数量级上的差异性：

MMCD^*＝(MMCD-MMCD_min)/(MMCD_max-MMCD_min) (3)

EVI^*＝(EVI-EVI_min)/(EVI_max-EVI_min) (4)

其中，MMCD_max为水汽运移成本距离最大值，MMCD_min为水汽运移成本距离最小值，EVI_max为增强型植被指数最大值，EVI_min为增强型植被指数最小值，MMCD^*为归一化后的水汽运移成本距离，EVI^*为归一化后的增强型植被指数；

子步骤2-3：对归一化后的水汽运移成本距离数据MMCD^*和归一化后的增强型植被指数数据EVI^*进行重采样至空间分辨率为25km；

基于重采样得到的2001至2013年内的归一化后的水汽运移成本距离数据MMCD^*和归一化后的增强型植被指数数据EVI^*建立与TRMM降水量的统计关系：

其中a_hv表示权重系数，h，v分别代表自变量的维数，TRMM为2001至2013年内的降水数据；通过最小二乘拟合得到权重系数a_hv；

子步骤2-4：将子步骤2-2中获取的分辨率为1km的归一化后的水汽运移成本距离和分辨率为1km的增强型植被指数数据代入公式(5)中，获取1km分辨率降水模拟值pre_1km；

子步骤2-5：将子步骤2-3中获取的分辨率为25km的归一化后的水汽运移成本距离和分辨率为25km的增强型植被指数数据代入公式(5)中，获取25km分辨率降水模拟值pre_25km；

子步骤2-6：将原始分辨率为25km的原始TRMM降水数据减去pre_25km获取空间分辨率为25km的残差residual_25km，并将所述残差residual_25km内插至1km分辨率残差residual_1km；

子步骤2-7：将子步骤2-4得到的1km分辨率降水模拟值pre_1km与子步骤2-6得到的1km分辨率残差residual_1km相加，得出空间分辨率为1km的降水数据。

进一步的，所述水汽运移成本距离MMCD和MODIS增强型植被指数EVI数据为降水数据辅助环境因子，时间尺度为2001年至2013年。

更进一步的，所述起始网格为海岸线。

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

受观测条件和与技术等因素限制，通常难以获取高空间分辨率的降水资料，由于降水的高时空异质性，导致现有技术中通常使用的环境因子在降水降尺度时依然面临许多挑战。本发明根据降水与环境因子的相关关系，提出了基于水汽运移成本距离和增强型植被指数的长时序遥感降水数据降尺度方法，相比现有技术，发明具有以下有益效果：(1)本发明可以有效解决降水数据降尺度对高分辨率数据的依赖性，在有限数据支持下获取长时间高分辨率降水数据成为现实；(2)本发明算法结构简单，具有很好的可操作性，避免了物理模型结构复杂、输入参数多等问题，具有较强的实用性与推广性。

附图说明

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

图2为基于本算法构建的水汽运移成本距离因子空间分布图。

图3为基于原始TRMM降水数据的河潼区间2000～2013年年均降水分布图。

图4为基于本算法的降水空间特分布图。

图5为基于本算法的降尺度结果与实测降水量散点图。

具体实施方式

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

本发明以美国国家航空航天局(NASA)和日本宇航局(JAXA)联合发布的2001至2013年TRMM卫星逐月降水数据产品(分辨率25km)为例，公开了耦合地表辅助参数的长时序卫星降水数据降尺度方法。本发明所采用降水数据辅助环境因子分别为水汽运移成本距离MMCD与MODIS EVI，时间尺度为2001年至2013年。

一、核心方法和过程

如图1所示，本发明主要包括两大步骤：(一)根据DEM和海岸线数据构建水汽运移成本距离因子，将其与MODIS增强型植被指数一同作为指示降水空间差异的辅助环境变量；(二)根据水汽运移成本距离、增强型植被指数两个地表环境变量与TRMM降水量的统计关系，由高分辨率环境参数得出其对应时间的高分辨率降水数据。

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

(一)利用DEM数据和海岸线数据构建“水汽运移成本距离”因子。

基本思路和原理是：我国主体大部分位于东亚季风区，夏季东南季风从西太平洋带来的暖湿气流是我国大部分区域水汽的主要来源。因此认为降水量大小的空间分布与水汽运移的距离存在相关性，已有研究将由西太平洋上空出发的含水大气通过大气环流直接进入到研究区内各个地区的欧氏(Euclidean)距离作为降水预测因子之一。

考虑到地形起伏对水汽运移造成的影响(阻抗)，假设水汽在运移过程中越过山区时含水量损失大于平原地区。因此，为更合理地表示水汽从海洋中运移到研究区各部分后的含水量损失状况，本发明将水汽运移路径的海拔高程作为水汽的运移成本，计算水汽从海洋通过大气环流运移到研究区各部分的成本距离(Moisture Migration CostDistance)，以此作为表征研究区降水差异的环境因子，在本发明中我们将该环境因子称为水汽运移成本距离(MMCD)。简言之，也即是根据降水的空间变异特征提出了基于海拔高程的“水汽运移成本距离”降水预测因子，该降水预测因子同时考虑了影响降水空间变异性的纬度位置、海陆位置、地形对降水的影响在表现降水的空间差异特征方面与已有的环境因子相比具有一定优势。

在具体操作中，将西太平洋海岸线设为起点，把研究区划分为1km大小的网格，以DEM为路径成本，计算来自海洋的水汽运移到研究区内各网格的成本距离。相邻像元成本距离计算方法为：

COST_m＝(DEM_i-DEM_j)/2 (1)

其中，COST_m、COST_n分别为以DEM为阻抗的邻接像元、对角像元成本距离，累积成本为起始像元移动到终止像元的成本距离之和，该累积成本即是水汽运移成本距离MMCD。成本距离计算中，要分配至输出栅格，像元必须具有到达某个源像元的次最低成本路径，因此这是起始于源像元的迭代过程，每个像元的目标是分配到输出成本距离栅格中，为所有像元执行此过程，直到遇到栅格的边、窗口的边界或达到最大距离为止。

(二)由水汽运移成本距离因子和MODIS增强型植被指数得出2001年至2013年内分辨率为1km的逐月分降水数据。

基本思路和原理在于：本领域公知的，降水量与植被覆盖、海陆位置等因素密切相关。基于该原理，本发明采用分辨率为1km的MODIS增强型植被指数(EVI)产品及步骤(一)构建的水汽运移成本距离(MMCD)因子，建立与降水量的经验方程，实现降水数据的降尺度处理。简言之，也即是根据MODIS增强型植被指数产品和“水汽运移成本距离”两个环境变量与降水量在空间上的相关关系建立降尺度转换函数，获取高分辨率的空间化降水数据。

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

MMCD^*＝(MMCD-MMCD_min)/(MMCD_max-MMCD_min) (3)

EVI^*＝(EVI-EVI_min)/(EVI_max-EVI_min) (4)

其中，MMCD_max，MMCD_min分别为最大与最小水汽运移成本距离，EVI_max，EVI_min分别为最大与最小增强型植被指数，MMCD^*和EVI^*为归一化后数据。

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

其中a_ij表示权重系数，i，j分别代表自变量的维数，一次至三次回归方程即可取得较好的结果。由于TRMM、MMCD和EVI三个变量均已知，通过最小二乘拟合即可得到权重系数a_ij，进而可利用该系数将原始的1km分辨率的MMCD及EVI数据代入公式(5)，获取1km分辨率的降水模拟值(pre_1km)。同时，将分辨率为25km的归一化后的MMCD与EVI数据代入公式(5)中，得到分辨率为25km的降水模拟值(pre_25km)。

然后将TRMM原始数据(分辨率为25km)减去25km的降水模拟值(pre_25km)，得出分辨率为25km的降水残差(residual_25km)，并通过空间内插方法将残差内插为1km分辨率(residual_1km)。最后将1km分辨率的降水模拟值(pre_1km)和1km分辨率的降水量残差(residual_1km)相加，获取最终的分辨率为1km的降水数据降尺度结果。

这里所依据的基本原理是：由于TRMM降水数据与EVI、MMCD数据存在较为稳定的相关性，本发明假设这种相关性在0.25°空间尺度和1km空间尺度上相同，并且可以通过某种函数模型来表示(指数模型、对数模型、线性模型等)。因此可以利用这种相关性对TRMM降水数据进行降尺度。

此外，需要说明的是，由于降水的高时空异质性，必然存在不能被EVI和MMCD有效指示的那部分降水量，即残差。本发明通过将这部分残差进行降尺度并与1km分辨率降水模拟值相加的方式来最小化模型误差。由于TRMM原始数据为0.25°，因此可通过TRMM数据值减去回归模型拟合值计算得出0.25°分辨率下的残差。由于残差产生的随机性，因此选择通过空间内插的方式对残差进行降尺度，在众多的空间内插方法中，推荐使用薄板样条函数插值法(Thin-plate Spline)对残差进行插值处理，来获取最优的降尺度结果。

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

本发明最后从空间分布趋势准确性和数值精度两个方面评价降尺度结果的有效性。在空间分布趋势准确性方面，本方法采用与原始25km分辨率的TRMM数据产品直接对比的方式来评估降尺度结果在保持原有降水空间差异及空间趋势特征上的有效性；在数值精度方面，通过计算国家基准和基本气象站的实测降水量与其所在降尺度结果像元的相关系数(CC)、相对偏差(BIAS)、均方根误差(RMSE)来评估降尺度结果在数值上的准确性，其计算方法为：

式中，y_i为实测降水量，均值为

是气象站点位置对应的降尺度结果像元值，均值为

i是观测站标识，n是观测站个数。CC反映了TRMM降水量与实测值的一致性，值域为0～1，数值越大反映出两者一致性越高；BIAS能表征降尺度结果对实测降水的偏离程度，BIAS大于0为正偏差，反之为负偏差；RMSE反映了数据在数值上差异大小，对离群值较为敏感，越接近0表明数据精度越高。

三、应用实例

以黄河中游河口镇至潼关区间为例，验证了本算法的效果。图2为利用本算法获得的基于高程的水汽运移成本距离降水辅助参数示意图，对比该因子与图3研究区的降水空间分布趋势图可以发现，水汽运移成本距离因子与降水量有明显的负相关关系，能够较好地反映出了研究区降水的空间差异与变化趋势特征。

图3和图4显示了2001～2013年原始TRMM年均降水数据与利用地表辅助参数尺度转换方法降水数据降尺度结果对比图。在降尺度后，TRMM降水数据空间分辨率得到明显提高，降水的空间分布变异趋势得到平滑，细节更为突出，降水随海拔高程、经纬度的变化趋势明显。可以看出，研究区内年降水量由东南向西北递减，最高值出现在南部的秦岭山麓，年降水量超过900mm，降水量最低值出现在研究区西北部的毛乌素沙地，年均降水量不足100mm，降尺度后的降水数据较好地保存了原始TRMM数据在反映降水空间差异及变化趋势。

为验证降尺度结果的数值准确性，以研究区内国家基准和基本气象站实测降水量为自变量，以对应像元的降尺度结果为因变量进行一元线性回归分析，结果如图5。检验结果表明，降尺度结果与气象站点实测值的决定系数R²＝0.85，算法结果与实测值具有较高的一致性。同时，计算降尺度结果与站点实测值的各评价指标可得，相关系数CC＝0.92，相对偏差BIAS＝5.00％，均方根误差RMSE＝43.97mm，说明算法结果与实测值存在极高的相关性。

表1

为进一步验证降尺度结果在研究区内误差的空间分布，计算各站点实测值与对应1km像元值的偏差，结果如表1(基于本算法的降尺度结果与实测降水量差异表)所示。根据表1可以知，降尺度结果与实测值较为接近，在研究区内的23个站点中，有18个站点误差小于50mm，表明降水数据降尺度方法取得较好效果。

此外，需要特别指出的是：本领域技术人员也可以利用2001年至2013年的MODISEVI数据及水汽运移成本距离数据将TRMM数据进行空间降尺度。

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

Claims (3)

1.一种长时间序列降水数据降尺度方法，其特征在于，它包括以下步骤：

(一)基于降水的空间差异特征构建水汽运移成本距离因子：

子步骤1-1：获取90m空间分辨率的STRM DEM数据，并将数据重采样至空间分辨率为1km；

子步骤1-2：获取全球海岸线矢量数据，并将矢量数据转换为1km分辨率的栅格数据；

子步骤1-3：将研究区划分为1km大小的网格，以海岸线作为起点，以DEM数据为路径成本，计算来自海洋的水汽运移到研究区内各网格的成本距离；各网格成本距离计算方法为：

COST_m＝(DEM_i-DEM_j)/2 (1)

其中，DEM_i和DEM_j分别为两个相邻网格的海拔高程；COST_m、COST_n分别为以DEM数据为阻抗的邻接像元成本距离和对角像元成本距离；

研究区内各网格的水汽运移成本距离定义为从起始网格到各终止网格的累积成本距离；累积成本距离为起始像元移动到终止像元的成本距离之和；

(二)由MODIS EVI数据和步骤(一)得出的水汽运移成本距离因子得出2001年至2013年内分辨率为1km的降水数据：

子步骤2-1：接收2001年至2013年内分辨率为1km的MODIS增强型植被指数EVI数据；

子步骤2-2：对水汽运移成本距离和EVI进行归一化处理，消除数量级上的差异性：

MMCD^*＝(MMCD-MMCD_min)/(MMCD_max-MMCD_min) (3)

EVI^*＝(EVI-EVI_min)/(EVI_max-EVI_min) (4)

其中，MMCD_max为水汽运移成本距离最大值，MMCD_min为水汽运移成本距离最小值，EVI_max为增强型植被指数最大值，EVI_min为增强型植被指数最小值，MMCD^*为归一化后的水汽运移成本距离，EVI^*为归一化后的增强型植被指数；

子步骤2-3：对归一化后的水汽运移成本距离数据MMCD^*和归一化后的增强型植被指数数据EVI^*进行重采样至空间分辨率为25km；

基于重采样得到的2001至2013年内的归一化后的水汽运移成本距离数据MMCD^*和归一化后的增强型植被指数数据EVI^*建立与TRMM降水量的统计关系：

其中a_hv表示权重系数，h，v分别代表自变量的维数，TRMM为2001至2013年内的降水数据；通过最小二乘拟合得到权重系数a_hv；

子步骤2-4：将子步骤2-2中获取的分辨率为1km的归一化后的水汽运移成本距离和分辨率为1km的增强型植被指数数据代入公式(5)中，获取1km分辨率降水模拟值pre_1km；

子步骤2-5：将子步骤2-3中获取的分辨率为25km的归一化后的水汽运移成本距离和分辨率为25km的增强型植被指数数据代入公式(5)中，获取25km分辨率降水模拟值pre_25km；

子步骤2-6：将原始分辨率为25km的原始TRMM降水数据减去pre_25km获取空间分辨率为25km的残差residual_25km，并将所述残差residual_25km内插至1km分辨率残差residual_1km；

子步骤2-7：将子步骤2-4得到的1km分辨率降水模拟值pre_1km与子步骤2-6得到的1km分辨率残差residual_1km相加，得出空间分辨率为1km的降水数据。

2.根据权利要求1所述的一种长时间序列降水数据降尺度方法，其特征在于，所述水汽运移成本距离MMCD和MODIS增强型植被指数EVI数据为降水数据辅助环境因子，时间尺度为2001年至2013年。

3.根据权利要求1或者2所述的一种长时间序列降水数据降尺度方法，其特征在于，所述起始网格为海岸线。

Images