一种河网结构约束下流域径流的随机时空插值方法

1.本发明涉及一种流域径流时空插值方法，具体涉及一种河网结构约束下流域径流随机时空插值方法，属于流域水文分析技术领域。


背景技术：

2.流域完整径流量序列的获取对于水利基础设施规划、水资源管理、洪水或干旱预报等水文研究必不可少。
3.然而，现实中水文观测数据由于意外或人为原因导致记录不全或信息损失，由此造成对水文设计工作可靠性的巨大威胁，因此亟需插值流域径流缺失数据。
4.受径流天然产汇流过程和河网结构下垫面因素的复杂作用，流域径流情势呈现出显著的空间区域差异和时间波动特性，致使传统的经验正交函数 (empirical orthogonal function；eof)方法难以有效地描述部分地区的径流时序演变模态，进而导致传统eof法插值精度不高的情况。
5.一方面，在河网结构影响下，流域产汇流过程在空间近邻点存在水文相似性，即流域径流在不同近邻点存在嵌套叠加效应，导致流域不同区域径流量序列可能存在高度相关性，如何削弱统计自相关性影响、分离多种径流时序演变模态来改进传统eof法的插值效率问题仍未能得到充分解决；另一方面，流域径流量序列插值具有时间和空间上的二维属性。
6.因此，有必要解决从单一时间或空间维度转向时空二维同步插值的问题。


技术实现要素：

7.为解决现有技术的不足，本发明的目的在于提供一种河网结构约束下流域径流随机时空插值方法。
8.为了实现上述目标，本发明采用如下的技术方案：
9.一种河网结构约束下流域径流随机时空插值方法，包括以下步骤：
10.s1、提取流域河网结构约束下已知m个站点的集水区的面积嵌套关系，计算已知站点的去嵌套集水区面积ωi及去嵌套集水区的径流量序列q(t,ωi)；
11.s2、将流域内不同站点的径流量序列q(t,ωi)进行均一化处理，得到均一化流量序列x(t,ωi)；
12.s3、采用eof法，对于均一化流量序列x(t,ωi)，计算不同去嵌套集水区出口i∈(1,...,m)和j∈(1,...,m)之间的空间方差-协方差矩阵covq(ωi,ωj)；
13.分解x(t,ωi)得到均一化流量序列的时间主成分ψk(t)；k＝1,...,m和空间函数βk(ωi)；
14.s4、通过消除每个去嵌套集水区的径流量序列q(t,ωi)对流域总出口的径流量序列q(t,ω)的线性相关性，创建新残余序列q
(r)
(t,ωi)；ω为流域总出口的集水面积；
15.将eof分解应用于q
(r)
(t,ωi)实施步骤s2和s3，以强制第1阶时间主成分恰等于流域总出口的径流量序列为条件，重新定义并计算均一化后相应的时间主成分φk(t)，记为
ceof法；
16.s5、对基于eof和ceof法分解得到的时间主成分进行归一化处理，使得不同阶时间主成分均满足均值为0，方差为1，记作ψ
′k(t)，计算归一化ψ
′k(t)对应的空间函数βk′
(ω(l))，得到河网上相对流域总出口距离l的任意子流域ω(l)的径流随机时空插值序列q(t,ω(l))；
17.s6、计算不同阶时间主成分对插值序列的相对解释方差评估河网结构上的不同空间点的径流插值精度。
18.上述步骤s1中：
19.a1、去嵌套集水区面积ωi的计算：
20.定义流域总出口的集水面积为ω，总出口的径流量序列表示为q(t,ω)；在流域河网结构约束下，提取流域内m个站点的去嵌套集水面积ωi；i＝1,...,m：
[0021][0022]
满足：
[0023][0024]
式中，s表示站点i上游所有站点数，ai表示站点i的实际总集水面积；若在河网结构下站点i上游没有任何已知站点，则ωi恰等于ai；
[0025]
a2、去嵌套集水区径流量序列q(t,ωi)的计算：
[0026]
为实际流量观测值减去该区上游入流的净流量值，
[0027][0028]
若在河网结构下站点i上游没有任何已知站点，则q(t,ωi)＝q(t,ai)。
[0029]
上述步骤s2中对径流量序列q(t,ωi)进行的均一化处理，包括中心化、标准化；
[0030]
即出口以上去嵌套集水区内的各空间点u的瞬时流量积分，可推求为均一化流量序列x(t,ωi)：
[0031]
中心化为：
[0032][0033]
标准化为：
[0034][0035]
式中，mq(ωi)表示去嵌套集水区的多年平均径流值，σq(ωi)表示去嵌套集水区径流的标准差；mq(ωi)、σq(ωi)由该站点已有的径流记录数据估算。
[0036]
上述步骤s3中eof法，包括：
[0037]
b1、x(t,ωi)线性分解为空间和时间函数的双正交序列组：
[0038][0039]
式中，时间主成分ψk(t)；k＝1,...,m与空间函数βk(ωi)都是阈值(-∞,+∞)内的
无量纲值，βk(ωi)代表某k阶时间主成分ψk(t)；k＝1,...,m在任意去嵌套集水区的权重系数集，满足：
[0040][0041]
式中，任意去嵌套集水区的均一化流量序列x(t,ωi)都可以看作是不同阶的波幅函数投影在该区权重向量上的函数值的线性组合；
[0042]
b2、计算均一化流量序列x(t,ωi)的空间方差-协方差矩阵covq(ωi,ωj)：
[0043]
对于中心化：
[0044][0045]
对于标准化：式(8)转化为空间相关系数ρ(ωi,ωj)的矩阵：
[0046]
covq(ωi,ωj)＝σq(ωi)σq(ωj)ρ(ωi,ωj)＝ρ(ωi,ωj)
ꢀꢀꢀ
(9)
[0047]
b3、求解βk(ωi)：
[0048][0049]
式中，标量λk；k＝1,2,...,m为特征值，通常按值的大小依次递减排序为λk≥λ
k+1
；
[0050]
b4、计算时间主成分ψk(t)：
[0051][0052]
式中，特征值λk也表示ψk(t)的方差值，即λk值越大，表示ψk(t)描述原始序列q(t,ωi)的能力越大，径流信息量损失越小。
[0053]
上述步骤s4中的ceof法，包括：
[0054]
首先依据每个去嵌套集水区的径流量序列q(t,ωi)与流域总出口的径流量序列q(t,ω)的线性相关性，生成新残余序列q
(r)
(t,ωi)：
[0055][0056]
将eof法应用于残余序列q
(r)
(t,ωi)，计算q
(r)
(t,ωi)的方差-协方差矩阵，得到时间主成分
[0057]
定义条件：ceof法的时间主成分φk(t)；k＝1,...,m如下：
[0058]
·
第1阶条件波幅函数为流域出口流量序列，即φ1(t)＝q(t,ω)；
[0059]
·
第2至m阶条件波幅函数依次等于由残余序列确定的m-1个波幅函数，即
[0060]
上述步骤s5中径流随机时空插值序列q(t,ω(l))的计算，包括：
[0061]
计算归一化时间主成分ψ
′k(t)对应的空间函数β
′k(ω(l))：
[0062]
对于中心化：
[0063]
β
′k(ω(l))＝ρ[q(t,ω(l)),ψ
′k(t)]σqꢀꢀꢀ
(13)
[0064]
对于标准化：
[0065]
β
′k(ω(l))＝ρ[q(t,ω(l)),ψ
′k(t)]
ꢀꢀꢀ
(14)
[0066]
式中，ρ表示相关系数，σq表示流域ω(l)径流的标准差；
[0067]
将β
′k(ω(l))和ψ
′k(t)代入式(6)，推导得到随机时空插值计算公式，如下：
[0068][0069]
式中，mq表示流域ω(l)的多年平均径流值；mq、σq以及相关系数ρ由该地点已有的径流记录数据估算或采取空间插值法(如克里金法、反距离加权法等)。
[0070]
上述步骤s6的插值精度的评估，包括：
[0071]
c1、计算k个时间主成分，即{ψ
′1(t),ψ
′2(t)，...,ψ
′k(t)}对于某流域ω(l)径流量序列的相对解释方差为：
[0072][0073]
式中，ω(l)为该流域上游与流域出口距离l的子流域的集水区面积，ρ[q(t,ω(l)),ψ
′k(t)]为q(t,ω(l))和ψk(t)之间的相关系数；
[0074]
当k＝m时，该地点的原始径流量序列可以被完全重构，即越大，表示保留原始数据信息越多；通常只使用少数(k《m)的时间主成分，一方面可以避免冗余信息，另一方面达到对流域径流插值精度的快速收敛；
[0075]
c2、采用效率系数指标，计算径流随机时空插值序列的精度；
[0076]
比如，nse效率系数(nash-sutcliffe)：
[0077][0078]
式中，nse∈[-∞,1]通过衡量长度为n的流量系列的估计(q
esi,t
)和观测系列 (q
obs,t
)的一致程度进行评价；
[0079]
对数化的nselog效率系数(nash-shuffle)：
[0080][0081]
式中，ln(
·
)表示自然对数函数；
[0082]
kge效率系数(kling-gupta)：
[0083][0084]
式中：kge代表了nse系数的分解形式，包含的三个部分依次是对估计序列相对于观测序列的pearson相关系数r，均值偏差(μ
esi
/μ
obs
)，和相对变异性(σ
esi
/σ
obs
)的衡量；
[0085]
c3、采用相对偏差指标re(％)，计算径流随机时空插值序列的统计特征参数θ
esi
相较于原始数据序列统计特征参数θ
obs
的偏差：
[0086]
假设无偏基准为100％，则
[0087][0088]
式中，θ表示上述统计特征参数的通用符号。
[0089]
进一步的，上述统计特征参数包括平均值、标准差、中位数、最小值和最大值。
[0090]
本发明的有益之处在于：
[0091]
本发明的一种河网结构约束下流域径流随机时空插值方法，
[0092]
(1)、通过考虑河网结构约束下流域内不同集水区的面积嵌套关系，有效提取既符合流域水文相似性，又保留不同集水区各自水文特征的径流量序列；
[0093]
(2)、通过分解流域多(水文)站点径流时空序列，有效剔除由于流域水文相似性所造成的不同集水区径流相关性的冗余信息，实现对流域径流量序列时间主成分和空间函数变化特征的单独分离和提取，便于分别诊断流域径流量序列的时间和空间变化规律，直观辨识多种时间主成分在流域内各空间点的演化形式；
[0094]
(3)、通过引入河网结构下不同集水区径流量相关性的强弱信息，灵活地反映各阶时间主成分和空间函数在解释沿河网结构不同径流时空变化规律上的相对代表性，随着输入时间主成分和空间函数的阶数增加，在力保原始数据信息损失最少的原则下，能够更好地提高传统eof法在流域径流时空序列插值上的精度，实现对径流量序列时间和空间上的二维随机插值。
[0095]
本发明的流域径流随机时空插值方法，适用于以地表河川径流为主要特征的流域，除了本实施例中展示的应用结果，本发明尤其对缺乏历史水文流量观测资料的流域集水区，提供了由已知站点径流信息移植到无资料区的一种有效径流随机时空插值方法，可根据不同无资料区径流的空间插值结果信息，由单一维度插值灵活拓展到时空二维的同步插值，为水利行业从事流域水文分析与水资源管理决策提供有力的理论及数据依据，具有很强的实用性和广泛的适用性。
附图说明
[0096]
图1为实测径流量序列中心化和标准化两种均一化处理结果。
[0097]
图2为eof+标准化方法确定的第2阶空间函数。
[0098]
图3为eof+中心化和eof+标准化方法求解出的第2阶时间主成分。
[0099]
图4为eof+中心化和eof+标准化方法求解出的第2阶时间主成分(图a) 及eof+中心化方法给出第2阶空间函数(图b)。
[0100]
图5为沿流域总出口距离l处径流与第1阶时间主成分(eof+标准化方法) 之间相关系数ρ的插值结果。
[0101]
图6为ceof+标准化方法插值结果与原始实测序列对比图。
[0102]
图7为不同插值方法的相对解释方差随不同阶数时间主成分的变化情况。
[0103]
图8为不同插值方法效果的评价指标值对比图。
具体实施方式
[0104]
以下结合附图和具体实施例对本发明作具体的介绍。
[0105]
本实施例的一种河网结构约束下流域径流随机时空插值方法，包括以下步骤：
[0106]
s1、提取流域河网结构约束下已知的m个站点的不同集水区的面积嵌套关系，计算已知站点的去嵌套集水区面积ωi及去嵌套集水区径流量序列q(t,ωi)：
[0107]
a1、去嵌套集水区面积ωi的计算：
[0108]
定义流域总出口的集水面积为ω，总出口的径流量序列为q(t,ω)；在流域河网结构约束下，提取流域内m个站点的去嵌套集水面积ωi；i＝1,...,m：
[0109][0110]
满足：
[0111][0112]
式中，s表示站点i上游所有站点数，ai表示站点i的实际总集水面积；若在河网结构下站点i上游没有任何已知站点，则ωi恰等于ai。
[0113]
a2、去嵌套集水区径流量序列q(t,ωi)的计算：
[0114]
即实际流量观测值减去该区上游入流的净流量值，
[0115][0116]
若在河网结构下站点i上游没有任何已知站点，则q(t,ωi)＝q(t,ai)。
[0117]
s2、将流域内不同站点的长序列径流量各时刻t的q(t,ωi)进行均一化处理，具体可选择中心化或者标准化处理，即出口以上去嵌套集水区ωi内各空间点u的瞬时流量积分可推求为均一化流量序列(中心化或标准化流量序列)：
[0118]
中心化为：
[0119][0120]
标准化为：
[0121][0122]
式中，mq(ωi)表示去嵌套集水区的多年平均径流值，σq(ωi)表示去嵌套集水区径流的标准差；mq(ωi)、σq(ωi)由该站点已有的径流记录数据估算。
[0123]
如图1所示，根据我国赣江流域内站点的径流资料情况，选取流域内23个站点划分集水区。以流域总出口为例，展示了中心化和标准化两种序列均一化处理结果。
[0124]
s3、采用经验正交函数(eof)方法，对于中心化(或标准化)流量序列 x(t,ωi)，计算不同去嵌套区出口i∈(1,...,m)和j∈(1,...,m)之间的空间方差-协方差矩阵covq(ωi,ωj)，分解x(t,ωi)得到中心化(或标准化)流量序列的时间主成分ψk(t)；k＝1,...,m和空间函数βk(ωi)。具体为：
[0125]
b1、x(t,ωi)线性分解为空间和时间函数的双正交序列组：
[0126][0127]
式中，时间主成分ψk(t)；k＝1,...,m与空间函数βk(ωi)都是阈值(-∞,+∞)内的无量纲值，βk(ωi)代表某k阶时间主成分ψk(t)；k＝1,...,m在任意去嵌套集水区ωi的权重系数集，满足：
[0128][0129]
式中，任意去嵌套集水区ωi的中心化(或标准化)流量序列x(t,ωi)都可以看作是不同阶的波幅函数投影在该区权重向量上的函数值的线性组合。
[0130]
b2、计算流量序列x(t,ωi)的空间方差-协方差矩阵covq(ωi,ωj)：
[0131]
对于中心化：
[0132][0133]
对于标准化：式(8)转化为空间相关系数ρ(ωi,ωj)的矩阵：
[0134]
covq(ωi,ωj)＝σq(ωi)σq(ωj)ρ(ωi,ωj)＝ρ(ωi,ωj)
ꢀꢀꢀ
(9)。
[0135]
b3、求解βk(ωi)：
[0136][0137]
式中，标量λk；k＝1,2,...,m为特征值，通常按值的大小依次递减排序为λk≥λ
k+1
；
[0138]
如图2所示，以eof(标准化)确定的第2阶空间函数为例，三角形尺寸反映了空间函数绝对值的大小。由图结果可以发现，分解得到的空间函数在流域内存在显著差异，反映了流域径流空间分布的变化情况。
[0139]
b4、计算时间主成分ψk(t)：
[0140][0141]
式中，特征值λk也表示ψk(t)的方差值，即λk值越大，表示ψk(t)描述原始序列q(t,ωi)的能力越大，径流信息量损失越小。
[0142]
如图3所示，以eof+中心化和eof+标准化方法为例，展示了两种方法求解出的第2阶时间主成分。结果代表了流域实测径流量序列的一种时间演变形式。
[0143]
s4、通过消除每个去嵌套集水区ωi；i＝1,...,m的原始观测流量序列对流域总出口的径流量序列q(t,ω)的线性相关性，创建一组新残余序列q
(r)
(t,ωi)取代q(t,ωi)，将eof分解应用于q
(r)
(t,ωi)实施步骤s2和s3，以强制第1阶时间主成分恰等于流域总出口的流量序列为条件，重新定义并计算中心化和标准化两种情况下相应的时间主成分φk(t)，记为条件ceof(complex empirical orthogonal function) 分解法。具体为：
[0144]
在ceof法中，首先依据每个去嵌套集水区ωi；i＝1,...,m的径流量序列q(t,ωi) 与流域总出口的径流量序列q(t,ω)的线性相关性，生成一组新残余序列q
(r)
(t,ωi)：
[0145][0146]
将eof分解应用于残余序列q
(r)
(t,ωi)，计算q
(r)
(t,ωi)的方差-协方差矩阵(针对中心化处理)或者相关系数矩阵(对于标准化处理)，求解得到时间主成分
[0147]
定义条件：ceof法的时间主成分φk(t)；k＝1,...,m如下：
[0148]
·
第1阶条件波幅函数为流域出口流量序列，即φ1(t)＝q(t,ω)；
[0149]
·
第2至m阶条件波幅函数依次等于由残余序列确定的m-1个波幅函数，即
[0150]
如图4(a、b)所示，以ceof+中心化和ceof+标准化方法为例，展示了两种方法求解出的第2阶时间主成分的结果，并给出了ceof+中心化方法得到的第2阶空间函数。
[0151]
s5、对基于eof和ceof法分解得到的时间主成分进行归一化处理，使得不同阶时间主成分均满足均值为0，方差为1，记作ψ
′k(t)，计算归一化ψ
′k(t)对应的空间函数β
′k(ω(l))，得到河网上相对流域总出口距离l的任意子流域ω(l)的径流随机时空插值序列q(t,ω(l))。具体为：
[0152]
计算归一化时间主成分ψ
′k(t)对应的空间函数β
′k(ω(l))：
[0153]
对于中心化：
[0154]
β
′k(ω(l))＝ρ[q(t,ω(l)),ψ
′k(t)]σqꢀꢀꢀ
(13)
[0155]
对于标准化：
[0156]
β
′k(ω(l))＝ρ[q(t,ω(l)),ψ
′k(t)]
ꢀꢀꢀ
(14)
[0157]
式中，ρ表示相关系数，σq表示流域ω(l)径流的标准差；
[0158]
将β
′k(ω(l))和ψ
′k(t)代入式(6)，推导得到随机时空插值计算公式，如下：
[0159][0160]
式中，mq表示流域ω(l)的多年平均径流值；mq、σq以及相关系数ρ由该地点已有的径流记录数据估算或采取空间插值法(如克里金法、反距离加权法等)。
[0161]
以克里金法的统计估计结果为例，如图5所示，展示了沿流域总出口距离l 处径流与第1阶时间主成分(eof+标准化方法)之间相关系数ρ的插值结果。
[0162]
如图6所示，以流域内吉安站为例，展示了该站点径流量序列插值结果。由图可以看出，ceof+标准化方法计算出的插值结果与原始实测序列吻合度较高，插值效果较好。
[0163]
而对于当相对流域总出口距离l的任意子流域ω(l)落在无资料区时，相应的可采用空间插值法计算相关系数ρ。
[0164]
s6、计算不同阶时间主成分对插值序列的相对解释方差评估河网结构上不同空间点的径流插值精度。具体为：
[0165]
c1、计算k个时间主成分(即{ψ
′1(t),ψ
′2(t)，...,ψ
′k(t)})对于某流域ω(l)径流量序列的相对解释方差为：
[0166]
[0167]
式中，ω(l)表示该流域上游与流域出口距离l的子流域的集水区面积；ρ[q(t,ω(l)),ψ
′k(t)]为q(t,ω(l))和ψk(t)之间的相关系数。当k＝m时，该地点的原始径流量序列可以被完全重构，即越大，表示保留原始数据信息越多。通常只使用少数(k《m)的时间主成分，一方面可以避免冗余信息，另一方面达到对流域径流插值精度的快速收敛。
[0168]
如图7所示，为四种方法插值结果的相对解释方差随着时间主成分阶数的增加逐渐增大。结果表明，当时间主成分达到2阶时，80％的径流量序列演变信息可以被有效描述，并验证了ceof法插值效果优于传统eof法。
[0169]
c2、采用效率系数指标，计算径流随机时空插值序列的精度。
[0170]
比如，nash-sutcliffe效率系数(nse)：
[0171][0172]
式中，nse∈[-∞,1]通过衡量长度为n的流量系列的估计(q
esi,t
)和观测系列 (q
obs,t
)的一致程度进行评价。
[0173]
对数化的nash-shuffle效率系数(nselog)：
[0174][0175]
式中，ln(
·
)表示自然对数函数。
[0176]
kling-gupta效率系数(kge)：
[0177][0178]
式中：kge代表了nse系数的分解形式，包含的三个部分依次是对估计序列相对于观测序列的pearson相关系数r，均值偏差(μ
esi
/μ
obs
)，和相对变异性 (σ
esi
/σ
obs
)的衡量。
[0179]
c3、采用相对偏差指标re(％)(假设无偏基准为100％)，计算径流随机时空插值序列的统计特征参数θ
esi
(包括平均值、标准差、中位数、最小值和最大值)相较于原始数据序列统计特征参数θ
obs
的偏差：
[0180][0181]
式中，θ表示上述统计特征参数的通用符号。
[0182]
如图8所示，以效率系数和统计特征参数中的标准差和中位数为例，不同插值结果评价指标值证明了ceof法相对于eof法偏离100％的标准线最少，其中ceof+标准化的插值效果最优。
[0183]
以上显示和描述了本发明的基本原理、主要特征和优点。本行业的技术人员应该了解，上述实施例不以任何形式限制本发明，凡采用等同替换或等效变换的方式所获得的技术方案，均落在本发明的保护范围内。