CN111078678B - 一种基于多源信息融合与降尺度的卫星降水数据校正方法 - Google Patents

Abstract

本发明公开了一种基于多源信息融合与降尺度的卫星降水数据校正方法，包括：S1)获取实况降水数据及卫星反演降水产品；S2)根据实况降水数据和卫星反演降水产品进行最优插值；S3)获得低空间分辨率NDVI数据；S4)选取环境变量以及低空间分辨率NDVI数据，进行自变量筛选；S5)建立GWR降尺度模型；S6)利用所述GWR模型获得OI‑GWR校正的研究区降水格点数据；S7)采用10折交叉验证对OI‑GWR校正的研究区降水格点数据进行检验。本发明引入最优插值对原始卫星降水数据进行质量改进，获得更好的初始场；在初始场得到改进的基础上运用GWR方法进行降尺度，得到的降尺度数据在空间分辨率和精度两方面均得到提高。

Description

一种基于多源信息融合与降尺度的卫星降水数据校正方法

技术领域

本发明涉及水文及气象数据分析领域，尤其涉及一种基于多源信息融合与降尺度的卫星降水数据校正方法。

背景技术

降水是地球水循环的主要驱动因子。准确获取降水时空分布信息是提高地球系统功能以及更好地预测天气、气候和自然灾害事件的根本。当前获取降水信息的方式主要有3种：台站实况观测、雷达估测和卫星遥感反演。站点观测具有较高的精度，但台站大多位于低海拔地带，且分布不均，不能有效地反映降水的空间变化特征，空间代表性不足；基于气象雷达回波推算的降水数据具有一定区域估测的优点，然而Z-I关系在不同的降水系统、不同的季节、不同的区域都有很大差异，影响了雷达估测降水的准确性，此外，雷达探测易受到地物阻挡、波束效应以及速度模糊等限制而使降水估测精度不高，实际应用中仍然存在较大限制；雨量计和雷达在陆地上分布不均，且在海洋上分布更加稀少，很难准确地获得大区域和全球性的降水分布。与之相比，卫星反演降水产品具有全天候、全球覆盖以及较好反映降水“空间分布”的独特优势。然而精度低与空间分辨率不高是卫星降水产品的两个主要缺陷，因此有必要在应用之前对其进行校正及降尺度研究。

在现有的卫星降水降尺度研究中，一般是引入环境变量建立与降水的函数映射实现卫星降水产品空间分辨率的降解，由于这些都是以卫星降水作为初始场，因此卫星降水产品的质量成为影响最终降尺度结果的主要因素。然而现有的卫星降水降尺度研究中，大多是在原始卫星降水的基础上运用各种统计模型达到降尺度的目的。比如在中国专利申请公开说明书CN106776481A中公开了一种作用于卫星降水数据的降尺度校正方法，该发明首先基于混合地理加权回归法，构建低分辨率卫星降水数据与多个高分辨率变量之间的混合地理加权回归模型；其次，基于混合地理加权回归模型、多个高分辨率变量，计算具有高分辨率的降尺度降水数据；再次，将具有高分辨率的降尺度降水数据作为背景场，计算得到高分辨率的差值分布图；最后，根据差值分布图，对具有高分辨率的降尺度降水数据进行计算，得到卫星降尺度降水校正数据。再比如在中国专利申请公开说明书CN108154193A中公开了一种长时间序列降水数据降尺度方法，包括以下步骤：(一)基于降水的空间差异特征构建水汽运移成本距离因子；(二)由MODIS EVI数据和步骤(一)得出的水汽运移成本距离因子得出2001年至2013年内分辨率为1km的降水数据。然而这些方法均是在原始卫星降水产品基础上运用各种统计降尺度方法进行的研究，虽然得到的降尺度结果空间分辨率得到了改进，然而由于其以原始卫星降水数据为初始场，在降尺度前并未将原始卫星降水产品进行质量改进，精度提高有限。因此需要提供一种基于多源信息融合与降尺度相结合的卫星降水数据校正方法，在改进降尺度结果空间分辨率的同时，提高降水数据的精度。

发明内容

本发明的目的在于为了解决雨量计和雷达很难准确地获得大区域和全球性的降水分布，并且现有卫星降水产品精度低与空间分辨率不高的问题，提供一种基于多源信息融合与降尺度的卫星降水数据校正方法。本发明通过首先引入最优插值(OI,OptimumInterpolation)方法对原始卫星降水数据进行质量改进，以获得精度更高的初始场；在初始场的基础上运用地理加权回归(GWR，Geographically weighted regression)方法进行降尺度，最终得到更高精度、更高空间分辨率的格点降水数据集。

为了实现上述目的，本发明采用的技术方案如下：

一种基于多源信息融合与降尺度的卫星降水数据校正方法，包括步骤：

S1)获取由气象站点测得的实况降水数据以及与实况降水数据对应的卫星反演降水产品；

S2)将实况降水数据作为观测场，将卫星反演降水产品作为初估场，根据观测场和初估场进行最优插值，获得低空间分辨率OI融合降水数据；

S3)获取高空间分辨率NDVI数据，对高空间分辨率NDVI数据进行预处理，获得低空间分辨率NDVI数据；

S4)选取与降水密切相关的环境变量以及低空间分辨率NDVI数据，进行自变量的筛选；

S5)根据低空间分辨率OI融合降水数据和低空间分辨率NDVI数据建立GWR模型；

S6)利用GWR模型获得OI-GWR校正的研究区降水格点数据；

S7)采用10折交叉验证对所述OI-GWR校正的研究区降水格点数据进行检验。

本发明首先将初始的卫星反演降水产品与实况降水数据(OBS)进行OI最优插值，引入OI方法对原始卫星降水数据进行质量改进，以获得更好的初始场，得到校正后的低空间分辨率OI融合降水数据；其次将筛选得到的解释因子，即低空间分辨率NDVI数据作为自变量，以校正后的低空间分辨率OI融合降水数据为因变量，运用地理加权回归(GWR)模型对低空间分辨率OI融合降水数据进行空间降尺度，获得OI-GWR校正降水数据，即降尺度产品。本发明在OI方法的基础上运用GWR方法进行降尺度，最终得到更高精度、更高空间分辨率的格点降水数据集。

进一步的，步骤S2)中根据观测场和初估场进行最优插值，包括步骤：

S21)将分析半径设为X，以待分析卫星数据栅格格点k为中心，在分析半径内获得距所述卫星数据栅格格点最近的n个有效格点(卫星降水格点中包含至少一个气象站点)；

S22)利用n个有效格点参与最优插值，获得第k个卫星数据栅格格点最终降水分析值

i为有效格点序数，F_k为第k个卫星降水栅格格点初估值，

为第k个卫星降水栅格格点的第i个有效格点实况观测值，

为第k个卫星降水栅格格点的第i个有效格点卫星降水初估值，

为第k个卫星降水栅格格点的第i个有效格点上实况观测值与初估值的偏差在估计时分配的权重；

S23)求解线性方程组

得到权重

i＝1,2,...,n,

表示第i个有效格点与第j个有效格点之间的初估场误差协相关，

表示第i个有效格点与第j个有效格点之间的实况观测误差协相关，λ_i为第i个有效格点的实况观测误差的标准差

与初估场误差标准差

之比，λ_j为第j个有效格点的实况观测误差的标准差

与初估场误差标准差

之比。

在站网稀疏区域需不断调整分析半径，以保证在该区域能搜索到一定数量的有效格点，有效格点即气象站点，再从中选取距分析格点最近的几个有效格点参与最优插值。本发明采用的最优插值方法以卫星反演降水产品作为初估场，以气象站点实况降水为真值，每个格点上的最终降水分析值A_k等于该点的初估值F_k加上该格点上实况观测值与初估值的偏差

而这个偏差

由一定范围内n个格点上已知的实况观测值

与初估值F_k的偏差

加权估计得到，W_i为权重函数，表示第k个卫星降水栅格格点的第i个有效格点上实况观测值与初估值的偏差在估计时分配的权重。通过求解线性方程组定权重回归系数

然后计算出最终降水分析值A_k。

进一步的，步骤S3)中对高空间分辨率NDVI数据进行预处理，包括步骤：

S31)剔除积雪和水体区域的NDVI异常值，获得剔除异常值后的NDVI数据；

S32)采用预设网格对步骤S31)中得到的NDVI数据进行滑动平均滤波，得到平滑的NDVI数据；

S33)对平滑的NDVI数据进行重采样，获得低空间分辨率NDVI数据。

针对剔除异常值后的NDVI数据用网格滑动平均进行低通滤波，以消除人类活动和各种自然因素造成的影响，然后得到平滑的NDVI数据。

进一步的，步骤S32)中设置阈值范围[-0.15,0.15]，根据阈值范围并采用预设大小的网格对步骤S31)中得到的NDVI数据进行滑动平均滤波，获得在阈值范围内平滑的NDVI数据。

滑动平均滤波是将在阈值范围外的NDVI数据作为异常值，首先选定某一尺寸的窗口，即预设网格的尺寸，然后将窗格内的NDVI数据做算术平均，将所求的平均值作为当前窗口中心点的NDVI数据。最后以[-0.15,0.15]为阈值进行筛选，获得平滑的NDVI数据。

进一步的，步骤S4)中选取与降水密切相关的环境变量以及低空间分辨率NDVI数据，进行自变量的筛选，包括步骤：

S41)与降水密切相关的环境变量包括高空间分辨率DEM数据，利用像素平均方法将高空间分辨率DEM数据重采样至低空间分辨率DEM数据，低空间分辨率DEM数据与低空间分辨率NDVI数据的空间分辨率值相同；

S42)从低空间分辨率DEM数据中提取环境变量数据，环境变量数据包括坡度、坡向、遮蔽度、曲率和/或地形湿度指数，将环境变量数据以及所述低空间分辨率NDVI数据作为解释因子；

S43)将解释因子进行不同变量组合，采用VIF方法对不同变量组合进行筛选，获得若干个防多重共线性处理后的变量组合；

S44)建立进行自变量筛选过程中的GWR模型，将若干防多重共线性处理后的变量组合分别带入所述自变量筛选过程中的GWR模型，得到若干个降尺度结果；

S45)对若干个降尺度结果进行对比评估，选取获得最佳降尺度效果的变量组合为最终的解释变量。

本发明选取了与降水密切相关的环境变量以及低空间分辨率NDVI数据，将环境变量数据以及低空间分辨率NDVI数据作为解释因子，对解释因子中的变量进行筛选，选出降尺度效果最好的变量组合，该变量组合可以是一个变量，也可以是多个变量。进行筛选时，先利用方差膨胀系数(Variance inflation factor，VIF)进行变量组合的筛选，方差膨胀系数值越接近于1，多重共线性越轻，反之越重。比如进行实验时从解释因子中逐一加入变量作为一个变量组合，计算其方差膨胀系数，若增加一个变量后导致多重共线，则把该变量删除。本发明通过逐一实验以及不同变量组合，通过VIF方法防止变量的多重共线性。获得若干个防多重共线性处理后的变量组合后，以卫星降水为因变量，建立进行自变量筛选过程中的GWR模型，将每个防多重共线性处理后的变量组合分别带入进行自变量筛选过程中的GWR模型中，获得与变量组合相应的降尺度结果，然后进行降尺度结果的比较，选取获得最佳降尺度效果的变量组合为最终的解释变量。

进一步的，最终的解释变量为低空间分辨率NDVI数据。通过大量实验测试表明当防多重共线性处理后的变量组合仅为低空间分辨率NDVI数据时，取得的降尺度效果最佳。

进一步的，步骤S5)中根据低空间分辨率OI融合降水数据和低空间分辨率NDVI数据建立GWR模型，包括步骤：

S51)将低空间分辨率OI融合降水数据作为因变量，将低空间分辨率NDVI数据作为自变量，根据因变量和自变量建立GWR模型；

S52)获得研究区高空间分辨率格点经纬度坐标，将研究区高空间分辨率格点经纬度坐标带入步骤S51)中GWR模型，获得研究区高空间分辨率格点的截距和NDVI回归系数；

S53)将高空间分辨率NDVI数据应用于所述高空间分辨率格点的截距和NDVI回归系数，获得高空间分辨率的基于NDVI的模型预测降水值；

S54)利用Spline插值获得NDVI异常值所对应的高空间分辨率降水估测值，得到研究区最终的高空间分辨率格点降水数据。

本发明将低空间分辨率OI融合降水数据作为因变量，将低空间分辨率NDVI数据作为自变量，根据因变量和自变量建立GWR模型

为第j个卫星数据栅格格点的低空间分辨率OI校正降水数据；

是围绕第j个卫星数据栅格格点的第i个有效格点的低空间分辨率NDVI数据；β₀(u_j,v_j)和β_i(u_j,v_j)分别为第j个卫星数据栅格格点的截距和斜率，(u_j,v_j)为第j个卫星数据栅格格点在二维空间的坐标，ε_j为残差。然后获取第j个卫星数据栅格格点的斜率β_i(u_j,v_j)＝(X^T(W(u_j,v_j))X)^-1X^TW(u_j,v_j)Y，X为所述自变量，Y为所述因变量，W(u_j,v_j)为权重矩阵。将研究区高空间分辨率格点经纬度坐标带入GWR模型，获得研究区高空间分辨率格点的截距和NDVI回归系数。接着将所述高空间分辨率NDVI数据应用于所述高空间分辨率格点的截距和NDVI回归系数，获得高分辨率的基于NDVI的模型预测降水值

最后将三次样条插值方法，即Spline插值应用于步骤S53)得到的高空间分辨率降水估测值以填充异常值所在空白，得到研究区最终的高空间分辨格点降水数据。计算权重矩阵

d_ij为第j个卫星数据栅格格点与第i有效格点的距离，b为带宽阈值。由权重矩阵能够确保第i个有效格点距离第j个卫星数据栅格格点越近则权重越大。

进一步的，步骤S7)中采用10折交叉验证将气象站点随机分为10组，每次选取每组中的9份用于与卫星反演降水产品的OI融合，剩余每组中的1份构成独立数据集进行低空间分辨率OI校正降水数据以及降尺度产品的精度验证。

将地面气象站点的降水观测数据融合到卫星降水数据的降尺度校正过程中，通过交叉验证，增加数值实验的次数来保障分析结果的稳定性，最大程度地利用了实测数据的优势。

进一步的，采用统计指标对降尺度产品进行精度验证，统计指标为平均绝对误差(MAE)、均方根误差(RMSE)和/或相关系数(CC)。

平均绝对误差

均方根误差

相关系数

式中：

n为样本容量，x_i和y_i分别为降水数据估测值和气象站观测值。

本发明的有益效果是：通过首先引入最优插值方法对原始卫星降水数据进行质量改进，以获得更佳的初始场；在初始场的基础上运用地理加权回归方法进行降尺度，最终得到更高精度、更高空间分辨率的格点降水数据集。本发明为区域降水数据集的开发提供一种更为有效的方法，即提出的OI-GWR两步校正方法先提高原始卫星降水产品精度再进行降尺度，获得了更高精度、更高空间分辨率的格点降水数据集，相关成果可为系统掌握区域降水时空分布特征提供有力的支撑。

附图说明

图1是本发明提出的OI-GWR方法获得天山山区格点月降水数据集DS_CIMERG的整体流程图；

图2是本发明实施例一获得降尺度数据DS_OIMERG的流程图；

图3是本发明实施例一获得降尺度数据DS_Spline的流程图；

图4是本发明实施例一2014-2018年暖季逐月OI处理前后各项指标对比图；

图5是本发明实施例一的三种不同方法逐月各项指标对比示意图；

图6是本发明实施例一的五种降水估测产品与实况降水整体及各月的散点图；

图7是本发明实施例一的五种降水估测产品的各项指标箱线图分布；

图8是本发明实施例一2016年7月五种降水估测数据与实况降水空间分布对比图；

图9是本发明实施例一伊犁河谷地区2016年7月五种降水估测数据与实况降水空间分布对比图；

图10是本发明提出的OI-GWR方法校正得到的研究区5月平均降水量的空间分布；

图11是本发明提出的OI-GWR方法校正得到的研究区6月平均降水量的空间分布；

图12是本发明提出的OI-GWR方法校正得到的研究区7月平均降水量的空间分布；

图13是本发明提出的OI-GWR方法校正得到的研究区8月平均降水量的空间分布；

图14是本发明提出的OI-GWR方法校正得到的研究区9月平均降水量的空间分布；

图15是本发明提出的OI-GWR方法校正得到的研究区年平均降水量的空间分布。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施方式仅仅用以解释本发明，并不用于限定本发明。

实施例一，如图1所示，一种基于多源信息融合与降尺度的卫星降水数据校正方法，本实施例一以新疆天山山区为研究区域。天山是世界上距离海洋最远的山系，世界七大山系之一，面积约为57万平方公里。作为南疆和北疆的地理分界线，天山山区平均海拔4,000m，由一系列山脉、山间盆地、谷地及山前平原等地貌单元组成。该区域属于我国典型的高寒山区，降水分布不均，时间和空间差异悬殊，夏季降水量大，而冬季降水量较小。本实施一例针对GPM(Global Precipitation Measurement)的3级月降水产品IMERG(IntegratedMulti-satellit E Retrievals for GPM)开展校正研究。首先将初始的3级月降水产品(OIMERG)与实况降水(OBS)进行最优插值(OI)，得到校正后的10km校正数据CIMERG；其次将1km NDVI重采样至10km作为自变量，运用地理加权回归(GWR)对10km校正数据CIMERG进行空间降尺度，得到1×1km的天山山区格点月降水数据集DS_CIMERG。包括以下步骤：

S1)获取由气象站点测得的实况降水数据以及与实况降水数据对应的卫星反演降水产品；

所用卫星降水产品来自全球降水测量计划(即GPM，Global PrecipitationMeasurement Mission)，是由美国国家宇航局(NASA)开展的建立在热带降水测量(Tropical Rainfall Measurement Mission，TRMM)基础上的卫星遥感降水测量计划，其目的是提供精度和分辨率更高的新一代全球卫星遥感数据产品。选取GPM的3级产品IMERG为实验数据，时间段为2014-2018年暖季(5-9月)。

S2)将实况降水数据作为观测场，将卫星反演降水产品作为初估场，根据观测场和初估场进行最优插值，在站网稀疏区域需不断调整分析半径，以保证在该区域能搜索到一定数量的有效格点，有效格点即气象站点，再从中选取距分析格点最近的几个有效格点参与最优插值。包括步骤：

S21)本实施例将分析半径设为100km，以待分析卫星数据栅格格点k为中心，在分析半径内获得距卫星数据栅格格点最近的9个有效格点，卫星降水格点中包含至少一个气象站点；

S22)利用n个有效格点参与最优插值，n＝9，获得第k个卫星数据栅格格点最终降水分析值

i为有效格点序数，F_k为第k个卫星降水栅格格点初估值，

为第k个卫星降水栅格格点的第i个有效格点实况观测值，

为第k个卫星降水栅格格点的第i个有效格点卫星降水初估值，

为第k个卫星降水栅格格点的第i个有效格点上实况观测值与初估值的偏差在估计时分配的权重；

S23)求解线性方程组

得到权重

表示第i个有效格点与第j个有效格点之间的初估场误差协相关，

表示第i个有效格点与第j个有效格点之间的实况观测误差协相关，λ_i为第i个有效格点的实况观测误差的标准差

与初估场误差标准差

之比，λ_j为第j个有效格点的实况观测误差的标准差

与初估场误差标准差

之比。

依次获得所有卫星数据栅格格点最终降水分析值，即获得10km分辨率OI融合降水数据。

S3)获取1km分辨率NDVI数据，对1km分辨率NDVI数据进行预处理，获得10km分辨率NDVI数据，包括步骤：

S31)剔除积雪和水体区域的NDVI异常值，获得剔除异常值后的NDVI数据；

S32)设置阈值范围[-0.15,0.15]，根据阈值范围并采用11×11网格对步骤S31)中得到的NDVI数据进行滑动平均滤波，得到平滑的NDVI数据；

S33)对平滑的NDVI数据进行重采样，获得10km分辨率NDVI数据。

对1km分辨率NDVI数据进行预处理，剔除积雪和水体区域的NDVI异常值，提高后续数据在使用过程中的精度。针对剔除异常值后的NDVI数据用网格滑动平均进行低通滤波，以消除人类活动和各种自然因素造成的影响，然后得到平滑的NDVI数据。

S4)选取与降水密切相关的环境变量以及10km分辨率NDVI数据，进行自变量的筛选，包括步骤：

S41)与降水密切相关的环境变量提取自30m分辨率DEM数据，利用像素平均方法将30m分辨率DEM数据重采样至10km分辨率DEM数据；

S42)从10km分辨率DEM数据中提取环境变量数据，所述环境变量数据包括坡度、坡向、遮蔽度、曲率和地形湿度指数，将环境变量数据以及10km分辨率NDVI数据作为解释因子；

S43)将解释因子进行不同变量组合，采用VIF方法对不同变量组合进行筛选，获得若干个防多重共线性处理后的变量组合；

S44)建立进行自变量筛选过程中的GWR模型，将若干防多重共线性处理后的变量组合分别带入所述自变量筛选过程中的GWR模型，得到若干个降尺度结果；

S45)对若干个降尺度结果进行对比评估，选取获得最佳降尺度效果的变量组合为最终的解释变量。最终的解释变量为10km分辨率NDVI数据。通过大量实验测试表明当防多重共线性处理后的变量组合仅有10km分辨率NDVI数据时，取得的降尺度效果最佳。

S5)根据10km分辨率OI融合降水数据和10km分辨率NDVI数据建立GWR模型，包括步骤：

S51)将10km分辨率OI融合降水数据作为因变量，将10km分辨率NDVI数据作为自变量，根据因变量和自变量建立GWR模型；

S52)获得研究区高空间分辨率格点经纬度坐标，将研究区高空间分辨率格点经纬度坐标带入步骤S51)中GWR模型，获得研究区高空间分辨率格点的截距和NDVI回归系数；

S53)将高空间分辨率NDVI数据应用于步骤S51)得到的高空间分辨率格点的截距和NDVI回归系数，获得1km分辨率的基于NDVI的模型预测降水值；

S54)利用Spline插值获得NDVI异常值所对应的高空间分辨率降水估测值，得到研究区最终的1km分辨率格点降水数据。

S6)利用GWR模型获得OI-GWR校正的研究区降水格点数据。

S7)采用10折交叉验证对所述OI-GWR校正的研究区降水格点数据进行检验。

将气象站点随机分为10组，每次选取每组中的9份用于与IMERG降水的OI融合，剩余每组中的1份构成独立数据集进行OI融合降水数据精度验证以及降尺度产品的精度验证。将地面气象站点的降水观测数据融合到卫星降水数据的降尺度校正过程中，通过交叉验证，增加数值实验的次数来保障分析结果的稳定性，最大程度地利用了实测数据的优势。采用统计指标对降尺度产品进行精度验证，统计指标为平均绝对误差(MAE)、均方根误差(RMSE)和相关系数(CC)。

平均绝对误差

均方根误差

相关系数

式中：

n为样本容量，x_i和y_i分别为降水数据估测值和气象站观测值。

为验证本发明提出方法的有效性，实施例一将其与传统的方法进行对比，传统的方法包括以往的GWR(Ori_GWR)降尺度方法和样条插值。如图4所示2014-2018年暖季各月OI处理前后各项指标对比，可以看出，经过最优插值处理后降水数据精度明显比最优插值处理前降水数据精度高。

在第一步经OI处理后取得卫星降水产品质量改进的基础上，进一步验证GWR方法的降尺度效果。将本发明与以往方法进行对比，Ori_GWR降尺度方法以校正前初始的IMERG降水为因变量进行GWR降尺度得到降尺度数据DS_OIMERG(见图2)；运用Spline插值方法对校正前初始的IMERG降水进行降尺度得到降尺度数据DS_Spline(见图3)。

图5显示了三种方法得到的天山山区2014-2018年5-9月降尺度数据DS_CIMERG、DS_OIMERG和DS_Spline的各项统计指标对比图。图6为2014-2018年暖季OIMERG、CIMERG、DS_Spline、DS_OIMERG以及DS_CIMERG与实况降水的散点图。图7为2014-2018年暖季OIMERG、CIMERG、DS_Spline、DS_OIMERG以及DS_CIMERG的各项指标箱线图分布图。

从整体统计(图6)也可明显地看出GWR降尺度后降水数据取得了进一步改进，10km校正数据CIMERG的相关系数CC为0.616,均方根误差RMSE为26.56mm,平均绝对误差MAE为17.59mm；利用本发明得到1×1km的天山山区格点月降水数据集DS_CIMERG相关系数CC为0.635，均方根误差RMSE为25.93mm，平均绝对误差MAE为17.15mm，整体评估结果反映出对于经过本发明提出的改进后的降尺度方法所获得的降尺度数据，大大优于其他两种降尺度产品，同时也显著优于初始OIMERG降水。

图8以天山地区降水最多的2016年7月为例，分析了实况降水、OIMERG、CIMERG以及降尺度处理后的DS_OIMERG、DS_Spline和DS_CIMERG的月降水空间分布。8a-8f分别为天山地区2016年7月实况降水、原始卫星降水估测数据OIMERG，OI融合降水估测数据CIMERG、Spline插值降尺度数据DS_Spline、OI-GWR降尺度数据DS_CIMERG、Ori_GWR降尺度数据DS_OIMERG，图8中8b至8f这5种格点降水估测数据均表现出相同的空间分布趋势，显示伊犁河谷为降水较为集中且降水最大区域。从细节来看，天山地区2016年7月格点降水估测数据OIMERG值明显低估了实况降水，尤其是伊犁河谷地区，而经过OI后的格点降水估测数据CIMERG则将实况数据很好地融合进了OIMERG降水。通过本发明改进的降尺度方法(OI-GWR)得到的DS_CIMERG,不但显著提高了初始的IMERG降水的精度，同时空间分辨率也从初始的约10km提高到1km，反映出更为精细且与实况更为一致的降水分布。图9与图8数据相同，以图像放大的方式给出伊犁河谷地区在2016年7月各降水数据集的详细分布，9a-9f分别为伊犁河谷地区2016年7月实况降水、原始卫星降水估测数据OIMERG，OI融合降水估测数据CIMERG、Spline插值降尺度数据DS_Spline、OI-GWR降尺度数据DS_CIMERG、Ori_GWR降尺度数据DS_OIMERG。

最后给出了本发明方法得到的2014-2018年暖季的各月平均降水量的空间分布见图10～图14，2014-2018年暖季的年平均降水量的空间分布见图15。

通过采用本发明公开的上述技术方案，得到了如下有益的效果：

本发明首先将初始的卫星反演降水产品与实况降水数据进行OI最优插值，通过引入OI方法对原始卫星降水数据进行质量改进，以获得更好的初始场，得到校正后的低空间分辨率OI融合降水数据；其次将筛选得到的低空间分辨率NDVI数据作为自变量，运用地理加权回归模型对低空间分辨率OI融合降水数据进行空间降尺度，获得OI-GWR校正降水数据，即降尺度产品。本发明在OI方法的基础上运用GWR方法进行降尺度，最终得到更高精度、更高空间分辨率的格点降水数据集。最后运用交叉验证通过与传统的降尺度方法Ori_GWR以及Sline插值对比，显示本发明提出的OI_GWR方法具有明显的改进效果，从而为区域降水数据集的开发提供一种更为有效的方法。

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

Claims (8)

1.一种基于多源信息融合与降尺度的卫星降水数据校正方法，其特征在于，包括步骤：

S1)获取由气象站点测得的实况降水数据以及与所述实况降水数据对应的卫星反演降水产品；

S2)将所述实况降水数据作为观测场，将卫星反演降水产品作为初估场，根据所述观测场和所述初估场进行最优插值OI，获得低空间分辨率OI融合降水数据；

S3)获取高空间分辨率NDVI数据，对所述高空间分辨率NDVI数据进行预处理，获得低空间分辨率NDVI数据；

S4)选取与降水密切相关的环境变量以及所述低空间分辨率NDVI数据，进行自变量的筛选，获得最终的解释变量；

S5)根据所述低空间分辨率OI融合降水数据和所述最终的解释变量建立GWR模型，所述最终的解释变量为所述低空间分辨率NDVI数据；

S6)利用所述GWR模型获得OI-GWR校正的研究区降水格点数据；

S7)采用10折交叉验证对所述OI-GWR校正的研究区降水格点数据进行检验。

2.根据权利要求1所述的基于多源信息融合与降尺度的卫星降水数据校正方法，其特征在于，步骤S2)中根据观测场和初估场进行最优插值OI，包括步骤：

S21)将分析半径设为X，以待分析卫星数据栅格格点k为中心，在分析半径内获得距所述卫星数据栅格格点最近的n个有效格点；

S22)利用所述n个有效格点参与最优插值，获得第k个卫星数据栅格格点最终降水分析值

i为有效格点序数，F_k为第k个卫星降水栅格格点初估值，

为第k个卫星降水栅格格点的第i个有效格点实况观测值，F_i ^k为第k个卫星降水栅格格点的第i个有效格点卫星降水初估值，W_i ^k为第k个卫星降水栅格格点的第i个有效格点上实况观测值与初估值的偏差在估计时分配的权重；

S23)求解线性方程组

得到权重W_i ^k，i＝1,2,...,n,

表示第i个有效格点与第j个有效格点之间的初估场误差协相关，

表示第i个有效格点与第j个有效格点之间的实况观测误差协相关，λ_i为第i个有效格点的实况观测误差的标准差

与初估场误差标准差

之比，λ_j为第j个有效格点的实况观测误差的标准差

与初估场误差标准差

之比。

3.根据权利要求1所述的基于多源信息融合与降尺度的卫星降水数据校正方法，其特征在于，步骤S3)中对所述高空间分辨率NDVI数据进行预处理，包括步骤：

S31)剔除积雪和水体区域的NDVI异常值，获得剔除异常值后的NDVI数据；

S32)采用预设大小网格对步骤S31)中得到的NDVI数据进行滑动平均滤波，得到平滑的NDVI数据；

S33)对所述平滑的NDVI数据进行重采样，获得低空间分辨率NDVI数据。

4.根据权利要求3所述的基于多源信息融合与降尺度的卫星降水数据校正方法，其特征在于，步骤S32)中设置阈值范围[-0.15,0.15]，根据所述阈值范围并采用预设大小的网格对步骤S31)中得到的NDVI数据进行滑动平均滤波，获得在所述阈值范围内平滑的NDVI数据。

5.根据权利要求1或3所述的基于多源信息融合与降尺度的卫星降水数据校正方法，其特征在于，步骤S4)中选取与降水密切相关的环境变量以及低空间分辨率NDVI数据，进行自变量的筛选，包括步骤：

S41)所述与降水密切相关的环境变量包括高空间分辨率DEM数据，利用像素平均方法将所述高空间分辨率DEM数据重采样至低空间分辨率DEM数据,所述低空间分辨率DEM数据与所述低空间分辨率NDVI数据的空间分辨率值相同；

S42)从所述低空间分辨率DEM数据中提取环境变量数据，所述环境变量数据包括坡度、坡向、遮蔽度、曲率和/或地形湿度指数，将所述环境变量数据以及所述低空间分辨率NDVI数据作为解释因子；

S43)将所述解释因子进行不同变量组合，对所述不同变量组合采用VIF方法进行筛选，获得若干个防多重共线性处理后的变量组合；

S44)以卫星降水为因变量，建立进行自变量筛选过程中的GWR模型，将所述若干个防多重共线性处理后的变量组合分别带入所述自变量筛选过程中的GWR模型，得到若干个降尺度结果；

S45)对所述若干个降尺度结果进行对比评估，选取获得最佳降尺度效果的变量组合为最终的解释变量。

6.根据权利要求1所述的基于多源信息融合与降尺度的卫星降水数据校正方法，其特征在于，步骤S5)中根据所述低空间分辨率OI融合降水数据和所述最终的解释变量建立GWR模型，包括步骤：

S51)将所述低空间分辨率OI融合降水数据作为因变量，将所述低空间分辨率NDVI数据作为自变量，根据所述因变量和所述自变量建立GWR模型；

S52)获得研究区高空间分辨率格点经纬度坐标，将研究区高空间分辨率格点经纬度坐标带入步骤S51)中所述GWR模型，获得研究区高空间分辨率格点的截距和NDVI回归系数；

S53)将所述高空间分辨率NDVI数据应用于所述高空间分辨率格点的截距和NDVI回归系数，获得高空间分辨率的基于NDVI的模型预测降水值；

S54)利用Spline插值获得NDVI异常值所对应的高空间分辨率降水估测值，得到研究区最终的高空间分辨率格点降水数据。

7.根据权利要求1所述的基于多源信息融合与降尺度的卫星降水数据校正方法，其特征在于，步骤S7)中采用10折交叉验证将气象站点随机分为10组，每次选取每组中的9份用于与卫星反演降水产品的OI融合，剩余每组中的1份构成独立数据集进行低空间分辨率OI校正降水数据以及降尺度产品的精度验证。

8.根据权利要求7所述的基于多源信息融合与降尺度的卫星降水数据校正方法，其特征在于，采用统计指标进行低空间分辨率OI融合降水数据以及降尺度产品的精度验证，所述统计指标为平均绝对误差、均方根误差和/或相关系数。

Images