一种基于主成分回归的多源降水融合方法与流程

1.本发明涉及降水空间估计技术领域，具体涉及一种基于主成分回归的多源降水融合方法。


背景技术：

2.降水是水分和能量循环的重要驱动因素，具有时空非连续性和差异性，是较难准确获取的一种气象变量。准确可靠的降水数据不仅是研究降水时空变化规律的关键，也是提高水文模拟精度的重要输入条件参数，对区域灾害监测、防洪减灾和水资源管理有着至关重要的作用。
3.长期以来，估计降水空间分布是基于地面雨量站点的观测数据来进行的，该方式的特点是测量点处精确度高且空间插值方法多样，然而地面站点观测的降水仅仅能够反映离散点降水信息，空间分布不均匀，受地形、环境、技术等因素影响。
4.随着科技的发展，基于卫星遥感观测降水成为了第二种途径，其具有覆盖面积广、时空分辨率高以及能够实现全球尺度降水观测的特点弥补了地面站点观测的短板，但是卫星降水产品精度差异大，适用性有待验证。
5.可以将空间连续的卫星降水数据与精确的地面站点观测数据相结合，为获取高质量高时空分辨率降水数据提供了新的思路。目前主流的融合方法有概率匹配、客观分析、贝叶斯技术、地理加权回归等，这些方法中大多数都是基于传统的多元回归模型，在降水空间估计的过程中会造成自变量间多重线性自相关的问题，导致计算的回归系数的符号及其数值与理论不一致，从而影响了结果的准确度。


技术实现要素：

6.发明目的：本发明的目的是提供一种能够解决融合过程中自变量之间存在的多重共线性问题，提高降水空间估计准确性的基于主成分回归的多源降水融合方法。
7.技术方案：本发明所述的一种基于主成分回归的多源降水融合方法，包括以下步骤：
8.(1)收集整理多源降水信息、匹配降水数据空间尺度、获取融合所需的dem数据、空间信息、地形因子数据，其中dem数据指数据高程模型，空间信息和地形因子数据为栅格的经度、维度、高程、坡度和坡向；
9.(2)通过步骤(1)中获取的信息和数据构建降水观测场p0和降水背景场pb；
10.(3)根据步骤(2)中获得的降水观测场、背景场信息，获取降水背景值误差e＝p0(i)-pb(i)，i为地面站点观测站，i＝1,2,3,...,n，基于主成分回归模型对背景值误差e拟合估计背景残差场拟合公式为：
[0011][0012]
在上式中，x
1i
,x
2i
,x
3i
,x
4i
,x
5i
为第i个网格处的经度、维度、高程、坡度和坡向，
为各变量的回归系数；
[0013]
(4)基于加法模型获取研究区域的多源融合降水空间分布。
[0014]
在上述技术方案中，考虑了宏观地形因子(经纬度、高程)与微观地形因子(坡度、坡向)的影响，解决了在融合了卫星降水产品与地面观测降水过程中自变量之间存在的多重共线性问题，有效提高了降水空间估计的准确性。
[0015]
优选的，在步骤(1)中，获取的多源降水信息包括地面站点观测降水数据和遥感反演降水数据。
[0016]
优选的，在获取了地面站点观测降水数据和遥感反演降水数据之后，确定研究区域降水估计的栅格空间分辨率，将不同遥感反演降水的空间分辨率统一。
[0017]
优选的，步骤(1)中，采用主成分回归法拟合背景残差场的步骤包括：(3.1)对自变量数据的标准化处理；(3.2)指标之间的相关性判定；(3.3)指标之间的相关矩阵特征值及对应的特征向量；(3.4)确定主成分的表达式；(3.5)成分与考察因变量的回归分析
[0018]
优选的，在步骤(4)中，未知点处的降水估计值通过背景残差场的估计值与背景场叠加得到，栅格降水估计值计算公式为：
[0019][0020]
在上式中，j为研究区域内任一处栅格。
[0021]
有益效果：本发明与现有技术相比，其具有的优点：本方案中采用了能够消除变量共线性问题的主成分回归法，构建了可同时集成空间连续的卫星降水产品数据与地面离散的雨量站点数据，并考虑局部流域区域化特征与地形、气象因子的多源降水信息融合模型，有效提高了降水空间估计的准确性。
附图说明
[0022]
图1为本发明的方法流程图；
[0023]
图2为实施例中研究区地理位置及地面水文气象站分布图；
[0024]
图3为黄河源区卫星降水产品精度评估折线图；
[0025]
图4为黄河源区融合降水结果精度评估折线图；
[0026]
图5为黄河源区卫星将水产品与站点观测的混淆矩阵；
[0027]
图6为黄河源区融合降水结果与站点观测的混淆矩阵。
具体实施方式
[0028]
下面结合附图对本发明的技术方案作进一步说明。
[0029]
实施例：在该实施例中，如图1所示，以我国黄河源区作为范例研究区，具体包括如下步骤：
[0030]
(1)收集整理多源降水信息、匹配降水数据空间尺度、获取融合所需的dem数据、空间信息、地形因子数据，其中，多源降水信息具体为：
[0031]
如图2所示，收集黄河源区2008.01.01-2013.12.31逐日气象站日降水观测资料，作为地面站点降水数据，其中黄河源区共有13个站点；收集覆盖黄河源区的四种卫星降水产品(tmpa 3b42v7、tmpa 3b42rt、cmorph_crt、imerg_final)，作为遥感反演降水数据；随
后采用双线性插值方法将获取的0.1
°×
0.1
°
imerg_final数据重采样到0.25
°×
0.25
°
空间分辨率，统一四种遥感反演数据的空间分辨率，最后利用srtm v4.1的90mdem数据，采用地理信息系统软件裁剪研究区域边界，提取研究个气象站与0.25
°×
0.25
°
栅格的经度、维度、高程、坡度、坡向信息；
[0032]
(2)通过地面站点观测降水数据和遥感反演降水数据构建黄河源区的降水观测场p0和降水背景场pb，具体为：
[0033]
构建2008.01.01-2013.12.31黄河源区13个气象站逐日降水观测场p0，tmpa 3b42v7、tmpa 3b42rt、cmorph_crt、imerg_final在250个栅格处的逐日降水背景场pb；
[0034]
(3)根据步骤(2)中获取的降水观测场和背景场信息，基于主成分回归模型构建空间范围内的任一没有地面雨量观测处的背景残差场具体为：
[0035]
首先基于降水观测场p0和降水背景场pb，计算降水背景值误差e：e＝p0(i)-pb(i)；
[0036]
在上式中，i为地面站点观测站，i＝1,2,3,...,n，黄河源区对应n为13；
[0037]
随后对于空间范围内的任一没有地面雨量观测处的位置j，降水背景值残差可以利用实测点降水数据与背景值之差e进行拟合估计，即：
[0038][0039]
在上式中，j为研究区域内任一处栅格，j＝1,2,...,k，黄河源区对应j为250；
[0040]
最后，采用主成分回归的方法来拟合背景残差场回归计算中使用的背景误差e作为因变量，气象站的经度、维度、高度、坡度和坡向作为输入变量，得到背景残差场：
[0041][0042]
在上式中，x
1i
,x
2i
,x
3i
,x
4i
,x
5i
为第i个网格处的经度、维度、高程、坡度和坡向，为各变量的回归系数。
[0043]
上述的主成分回归法的具体步骤为：对于n个降水观测数据，每个数据对应有p个变量x1,x2,...,x
p
，第k个观测数据可表示为(x
k-1
,x
k-2
,...,x
k-p
)，zi(i＝1,2,...,p)上的载荷为a
ij
(i＝1,2,...,p,j＝1,2,...,p)，步骤(3)中采用主成分回归方法拟合背景残差场的步骤包括：
[0044]
(3.1)对自变量数据的标准化处理：先得出样本标准差sj和样本均值xj，记sj为xj的样本标准差，即的样本标准差，即为xj的样本均值，即原始数据的标准化变换为x是标准化后的数据矩阵，
[0045]
(3.2)指标之间的相关性判定：x的相关系数即为x的协方差阵，
[0046]
在上式中，r为半正定矩阵；
[0047]
(3.3)指标之间的相关矩阵特征值及对应的特征向量：第一主成分z1反映最多的原p个变量的信息，即var(z1)方差最大，求第一主成分z1＝a1'x，即为求a1＝(a
11
,a
21
,...,a
p1
)'，使得在条件a1'a1＝1下var(z1)达到最大值，已知var(z1)＝var(a1'x)＝a1'ra1，r为x的相关系数阵，且的相关系数阵，且
[0048]
由上式可得，即ra1＝λ1a1，λ1为r的特征值，a1为λ1对应的单位正交特征向量，且根据var(z1)＝var(a1’
x)＝a1’
ra1＝a1’
λ1a1＝λ1可得，第一主成分的方差var(z1)达到最大值，为r的特征值的最大值λ1，而z1＝a1'x中的a1＝(a
11
,a
21
,...,a
p1
)'为λ1对应的正交特征向量，第二主成分方差var(z2)达到小于var(z1)的最大值，即为r的特征值λ2，且λ2＜λ1；
[0049]
设r的特征值λ1≥λ2≥
…
≥λ
p
≥0,a1,a2,...,a
p
是对应单位正交特征向量，记a＝(a1,a2,...a
p
)，满足aa'＝i，主成分z＝(z1,z2,...,z
p
)'，其中zi＝ai'x，根据var(z)＝a'xx'a＝a'ra＝λ，λ＝diag(λ1,λ2,...,λ
p
)，ra＝aλ，通过矩阵形式表示为变换后得
[0050]
(3.4)确定主成分的表达式：通过主成分回归可得p个主成分，第i个主成分的贡献率用表示，采用以下两个准则确定主成分个数：
[0051]
一是以累计贡献率确定，前m个主成分累计贡献率为该贡献率达到70％～85％时，选取前m个主成分；二是根据p个特征值的大小，假设p-m个特征值有λ
m 1
，λ
m 2
，
…
，λ
p
≈0，则p-m个特征值对应的主成分z
m 1
，z
m 2
，
…
，z
p
的贡献率小，舍去p-m个主成分，选取m个对应特征值偏大的主成分z1，z2，
…
，zm；
[0052]
(3.5)主成分与考察因变量的回归分析：对于p个主成分z1，z2，
…
，z
p
按照模型需求
及选取前m个主成分z1，z2，
…
，zm，该模型假设为y
t
＝b0 b1z
t1
b2z
t2
... bmz
tm
ε
t
,t＝1,2,...,n，e(ε
t
)＝0，var(ε
t
)＝σ2，cov(εi,εj)＝0，i≠j，记则回归模型的矩阵为y＝cb ε,ε～nn(0,σ2in)，在模型中，通过最小二乘法得到参数的估计值使得误差平方和达到最小即其中可得为b的最小二乘估计，由此得到y与z1，z2，
…
，zm的线性回归方程，根据的线性回归方程，根据将前m个主成分z1，z2，
…
，zm与x1，x2，
…
，x
p
之间的关系带入y＝cb ε得到y与自变量x1，x2，
…
，x
p
的线性回归方程。
[0053]
(4)基于加法模型，获取研究区域内的多源融合降水空间分布，具体为：
[0054]
采用加法模型叠加背景残差场和背景场，得到黄河源区2008～2013年逐日0.25
°×
0.25
°
栅格降水估计值pa(j)：
[0055][0056]
基于主成分回归的多源降水融合方法，得到四种地面观测-卫星产品融合降水，分别为tmpa 3b42v7_m，tmpa 3b42rt_m，cmorph_crt_m和imerf_final_m。采用黄河源区2008～2013年地面站点逐日降水观测作为基准数据，评估四种融合降水的精度。采用的精度评价指标包括：相关系数(correlation coefficient，cc)、偏差(relative bias，bias)、均方根误差(root mean squared error，rmse)三种定量指标，命中率(probability of detection，pod)、误报率(false alarm ratio，far)、临界成功指数(critical success index，csi)三种分类指标，此外，利用混淆矩阵(confusion matrix)评估融合结果对不同量级降水的探测能力。
[0057]
下表为本实施例的评价指标表：
[0058][0059][0060]
如图3和图4所示，分别给出了黄河源区遥感反演降水和融合降水的精度统计折线图。在定量指标方面，tmpa 3b42rt_m的cc、bias和rmse的范围提升至0.25～0.39/27.8％～168％、18mm～30.1mm，其余三种融合结果的范围提升至0.39～0.58、-33.1％～31.2％、8.7mm～15.6mm；对探测能力而言，融合后四种产品的pod提升至0.6～0.8，far降至0.3～0.56，csi最高值为0.57。总体而言，imerg_final_m的各统计指标均表现出较好的结果，对应cc、bias和rmse的均值为0.52、-1.6％、10.8mm；探测能力也明显高于其他融合结果，pod、
far和csi的均值为0.77、0.39、0.47。
[0061]
图5和图6分别给出了黄河源区遥感反演降水和融合降水重现站网观测能力的混淆矩阵。融合四种卫星降水产品与地面站点的降水数据能够更准确地探测各不同等级下的降水，对角线处网格的值也有所增加。融合后的各降水序列减轻了原始卫星数据在黄河源区降水量集中于0.1mm～1mm之间的现象，提高了与实测降水的相关性。具体而言，tmpa 3b42v7经融合处理后，较小等级降水准确度由0.43提升至0.77，此外中高类别降水(10mm～20mm、20mm～30mm、30mm～40mm)的准确度显著提高。tmpa 3b42rt_m的混淆矩阵相对于tmpa 3b42rt在无降水日和大雨等级(＞40mm)的准确度提升较明显，但和其余融合结果相较对角线处的统计指标仍然有很大的差距。cmorph_crt_m和imerg_final_m的混淆矩阵显示出较好的结果，尤其对于中高类别的降水(10mm～20mm、20mm～30mm)识别有所提升。
[0062]
本方案针对卫星降水空间连续与地面站点精确度高的特点，采用主成分回归方法融合了卫星降水产品与地面观测数据，有效地避免了自变量间多重共线性的问题。四种融合降水结果在黄河源区整体精度均有所提升，提高了对降水的探测率。本方案构建的融合方法也可以推广到其他流域，为利用卫星-地面降水数据融合进行水文气象研究提供了新的思路。
[0063]
以上显示描述了本发明的基本原理、主要特征、操作步骤。本领域技术人员应该了解，本发明不受上述优选实施例的限制，上述实例和说明书中描述的只是说明本发明的原理，凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。本发明要求保护范围由所附的权利要求书及其等效物界定。