一种基于天尺度洪涝指数提高海河流域洪涝监测时效性方法

1.本发明涉及卫星重力学、水文学交叉技术领域，尤其涉及一种基于天尺度洪涝指数提高海河流域洪涝监测时效性方法。


背景技术：

2.全球变暖的大背景下，水循环加快，导致洪涝和干旱的频率和严重程度增加，对地球生态和人类社会的发展造成了严重影响。洪涝已经成为影响经济发展的主要气象灾害，因此迫切需要及时、准确地监测洪涝灾害以减少人民的生命财产损失。
3.现有常用的洪涝检测方法主要分为两种：1)空间分布不均的现场站通过水位测量、河流流量监测可以准确记录水位、水量信息，但难以应对大范围洪涝覆盖。2)高分辨率的融合遥感技术能够持续、大规模地评估洪涝灾害，但它所需成本较高，而且不能探测土壤水分的响应。
4.土壤水对于地表水量变化高度敏感，甚至被当作洪水监测的替代数据。陆地水储量(terrestrial water storage，tws)是存在于地上和地下所有形态的水之总和，反映了水量输入(如降水、冰川融雪等)和输出(如人类耗水量、蒸散发等)之间的净效应。陆地水储量异常(terrestrial water storage anomaly，twsa)是相对某些年份tws的变化，如相对2004～2009年的tws变化。重力恢复和气候实验(gravity recovery and climate experiment，grace)卫星能够高精度地测量地球重力场，为监测大尺度空间的质量迁移提供了新手段。至今为止，grace时变重力场数据已被广泛应用于水储量变化、冰川消融、干旱、地壳形变等多领域研究。同时，grace反演的水储量变化在捕获洪水信息中发挥了重要作用。如：reager和famiglietti基于twsa首次提出了洪水潜力指数(flood potential index，fpi)，并通过与观测站数据进行比较，验证了其在全球范围内的适用性。随后，fpi得到广泛应用，如sun等利用fpi重现长江流域的洪水事件，证明了fpi能够以较高的时间和空间分辨率评估极端水文变化；idowu和zhou基于fpi评估了尼日利亚2012年洪水事件，其结果与洪水观测站的报告基本一致；huang等通过fpi确定了珠江流域2002～2019年间的四次洪水事件，并强调大型水库蓄水对长期干旱和洪水的影响。此外，改进的fpi也能够达到令人满意的效果。
5.另一方面，其他基于标准化思想的grace衍生指数也被用来监测极端气象所引起的水量变化，如sinha等提出基于grace twsa的联合气候偏差指数(combined climatologic deviation index，ccdi)探究印度主要流域的极端气候事件，结果表明ccdi在降雨量充沛的地区效率较高；nigatu等利用储水不足指数(storage deficit index，wsdi)评估了尼罗河流域的8次大洪水事件，通过与帕尔默干旱严重指数(pdsi)比较表明wsdi的结果更接近现实；xiong等基于多源grace数据得到集合fpi，并与wsdi进行比较，结果表明二者在黄河流域具有一致响应，在不同阈值下集合的fpi表现更优。
6.由于grace/grace-fo(后续grace或下一代grace)重力场模型粗糙的时间分辨率(周或月)难以捕捉到更精确的水文信号，因此，以上对洪涝的评估受限于月尺度。然而，目
前洪涝事件的形成可能只需要几个小时或几天，类似的洪水事件往往带有突发性特征。不同于前人已有发明，由于grace/grace-fo twsa粗糙的时间分辨率，基于twsa的洪涝指数仅能探测到月级洪涝灾害，并且难以在空间上准确刻画洪涝灾害的破坏范围。
7.因此，亟需一种基于天尺度洪涝指数提高海河流域洪涝监测时效性方法，以解决上述问题。


技术实现要素：

8.为缓解或解决上述问题中的至少一个方面或者至少一点，提出本发明。
9.本发明的一种基于天尺度洪涝指数提高海河流域洪涝监测时效性方法，包括如下步骤：
10.构建天尺度twsa：利用天尺度的降雨、温度数据及grace/grace-fo的尺度twsa数据构建天尺度twsa；
11.构建天尺度洪涝指数：利用去趋势后的重建的天尺度twsa构建天尺度洪涝指数；
12.天尺度洪涝指数包括fpi：
13.fpa
i,j
＝p
i,j-(max(twsa
i,j
)-twsa
i,j-1
)
[0014][0015]
其中，twsa
i,j
为去趋势后的twsa；p
i,j
为降雨；i为年；j表示某年的第j天。
[0016]
优选的，天尺度洪涝指数还包括ccdi、wsdi：
[0017]
wsd
i,j
＝twsa
i,j-mean(twsaj)
[0018]
pa
i,j
＝p
i,j-mean(p
i,j
)
[0019][0020][0021][0022]
pa
i,j
是降雨p相对于时序均值的异常，twsaj和paj分别为不同年份同一天的twsa和pa的集合；为pa
i,j
的残差，wsd
ij
和ccd
ij
分别是水储存赤字和组合气候偏差。
[0023]
优选的，利用天尺度的降雨、温度数据及grace/grace-fo的月尺度twsa数据构建天尺度twsa包括：
[0024]
明天的水储量等于今天经温度消耗剩下的水储量加上明日的降水:
[0025]
twsa(t)＝(twsa(t-1))e-1/τ(t)
+p(t)
[0026]
其中，t和p(t)分别表示日时间向量和每日降雨输入，e-1/τ(t)
为水储量的衰减函数，与温度相关，取值(0，1)；
[0027]
将降雨和温度驱动的结果平均成月，并作去趋势和去周期处理，最后利用同样去趋势和去周期的grace/grace-fo twsa进行约束：
[0028]
anom(grace(tm))＝β
·
anom(twsa(tm))+ε
[0029]
式中，β为约束因子，ε是残差项，tm表示月时间向量，anom表示从原始序列减去趋势项和周期项。
[0030]
优选的，还包括如下步骤：
[0031]
将twsa的时间序列分解为：
[0032]
twsa
total
＝twsa
trend
+twsa
seasonal
+twsa
residual
[0033]
其中，twsa
total
表示原始grace/grace-fo twsa信号，twsa
trend
是线性趋势分量，twsa
seasonal
表示季节性循环分量，包括年和半年循环；twsa
residual
代表残差分量，包括原始信号去趋势和季节性循环之后的剩余量；
[0034]
去趋势后的twsa为twsa
total-twsa
trend
。
[0035]
优选的，还包括如下步骤：
[0036]
利用最小二乘法将残差分量和其他分量分离：
[0037]
twsa＝a+b
·
t+c
·
cos(2θt)+d
·
sin(2θt)+e
·
cos(4θt)+f
·
sin(4θt)+ε
[0038]
其中,a～f代表各分量的参数，ε是残差项，对应twsa
residual
，a为常数项，t为相对2004的时间向量，a+b
·
t是趋势项，对应twsa
trend
，其余的则对应twsa
seasonal
。
[0039]
优选的，还包括如下步骤：将grace/grace-fo twsa分解的趋势和季节分量同样经过最小二乘法插值，得到对应的天尺度分量，直接加上残差分量，以此重建完整的天尺度twsa信号。
[0040]
优选的，还包括：天尺度twsa评价步骤：采用cc、rmse、nse，对构建的天尺度twsa进行评价。
[0041]
优选的，天尺度twsa评价步骤包括：将重建的天尺度twsa平均成月，比较其与grace/grace-fo twsa格网点之间的差异，并计算二者之间的cc和nse。
[0042]
优选的，还包括不同部分的土壤水对同一洪涝事件的响应分析步骤：将天尺度的sma、gldas-twsa和itsg-grace2018以及jpl-era5与重构的twsa进行比较，以评估它们在极端降雨下的水量变化，以上数据均去除了2004～2009的均值。
[0043]
优选的，还包括结果分析步骤：采用多种温度和降雨数据，以连续的grace/grace-fo twsa作为约束条件，重建天尺度twsa，将多种重建结果与jpl-rea5、itsg-grace2018和gldas-twsa比较来共同确定最优组合。
[0044]
本发明基于归一化和标准化的思想，利用重建的天尺度twsa提出了新型天尺度洪涝指数法，至少包括三种新的天尺度洪涝指数之一：洪水潜力指数(fpi)、储水不足指数(wsdi)以及联合气候偏差指数(ccdi)。本发明的天尺度洪涝指数法将洪涝事件的监测从传统的月尺度提升到了天尺度，较大程度提高了洪涝监测的时效性。
[0045]
本发明利用降雨和温度数据重建的陆地水储量异常(twsa)与重力恢复和气候实验(grace)卫星反演的twsa之间的nse达0.96，证明其具有极高的准确性。本发明新的天尺度洪涝指数法不仅保留了twsa对洪水的探测能力，而且在时间尺度上更加精确。
[0046]
本发明与2016年中国水旱灾害公报对比发现，新型天尺度洪涝指数可以有效探测到洪涝灾害的蔓延过程，并且对于有记录的亚月级洪水事件能够有效评估，验证了利用新型天尺度洪涝指数法监测洪涝事件的可靠性。新型天尺度洪涝指数fpi、wsdi和ccdi在时序上都能探测到海河流域2012年、2016年和2021年的洪涝事件，并且对于2016年的亚月级洪涝事件，fpi、wsdi和ccdi探测到7月19～20日两天内海河流域56％、43％和66％的地区遭受
洪水破坏，其中wsdi的结果与政府公开的信息相符。
[0047]
本发明为短时极端洪涝灾害监测提供了新的方法，证明了本发明天尺度洪涝指数法在洪涝监测方面的潜力，同时帮助流域水资源管理和预防及时做出决策。
附图说明
[0048]
图1海河流域概况：(a)位置和高程；(b)气象站点；
[0049]
图2重建月度twsa与原始twsa对比结果；
[0050]
图3不同降雨产品在海河流域的表现；
[0051]
图4不同温度在海河流域的变化；
[0052]
图5重建结果对比：(a)月尺度总twsa；(b)天尺度次季节分量；
[0053]
图6海河流域2016年7月17～7月21日降雨分布(来自gpm)；
[0054]
图7重建数据和其他产品在洪水中心的时序变化；
[0055]
图8重建数据和其他产品水量上的空间变化；
[0056]
图9洪涝指数时序变化；
[0057]
图102016年7月洪水事件中指数的空间变化；
[0058]
图11洪水中心各部分土壤水的响应；
[0059]
图12grace/grace-fo twsa与重建twsa的比较结果。
具体实施方式
[0060]
下述参照附图对本发明实施方式的说明旨在对本发明的总体发明构思进行解释，而不应当理解为对本发明的一种限制。在本发明中，相同的附图标记表示相同或相似的部件。
[0061]
在此描述的特征可以以不同的形式来实现，而不应被解释为限于在此描述的示例。相反，已提供在此描述的示例，以仅示出实现在此描述的方法、设备和/或系统的许多可行方式中的一些可行方式，所述许多可行方式在理解本发明的公开之后将是清楚的。
[0062]
尽管在此可使用诸如“第一”、“第二”和“第三”的术语来描述各种构件、组件、区域、层或部分，但是这些构件、组件、区域、层或部分不应被这些术语所限制。相反，这些术语仅用于将一个构件、组件、区域、层或部分与另一构件、组件、区域、层或部分进行区分。
[0063]
在说明书中，当元件(诸如，层、区域或基底)被描述为“在”另一元件上、“连接到”或“结合到”另一元件时，该元件可直接“在”另一元件上、直接“连接到”或“结合到”另一元件，或者可存在介于其间的一个或多个其他元件。相反，当元件被描述为“直接在”另一元件上、“直接连接到”或“直接结合到”另一元件时，可不存在介于其间的其他元件。
[0064]
在此使用的术语仅用于描述各种示例，并不将用于限制公开。除非上下文另外清楚地指示，否则单数形式也意在包括复数形式。术语“包含”、“包括”和“具有”说明存在叙述的特征、数量、操作、构件、元件和/或它们的组合，但不排除存在或添加一个或多个其他特征、数量、操作、构件、元件和/或它们的组合。
[0065]
为了使得本领域技术人员能够使用本发明的内容，下文中可能结合特定的应用场景、特定的系统、器件和元件的参数以及特定的连接方式，给出以下示例性实施例。然而，对于本领域技术人员来说，这些实施例仅是示例，在不脱离本发明的精神和范围的情况下，可
以将这里定义的一般原理应用于其他实施例和应用场景。
[0066]
根据本发明的一个示例性实施例，本发明提供一种基于天尺度洪涝指数提高海河流域洪涝监测时效性方法，包括如下步骤：
[0067]
构建天尺度twsa：利用天尺度的降雨、温度数据及grace/grace-fo的月尺度twsa数据构建天尺度twsa；构建天尺度洪涝指数：利用去趋势后的重建的天尺度twsa构建天尺度洪涝指数。还可进一步的包括，根据天尺度洪涝指数判断洪水发生的可能性。
[0068]
如图1所示：本发明以海河流域的洪水预测为例，说明本发明的工作过程。
[0069]
海河流域位于中国华北地区，东临渤海，地处112
°
e～120
°
e和35
°
n～43
°
n之间，总面积约32万km2。全流域总地势为西北高东南低，西部为黄土高原和太行山区，北部为蒙古高原和燕山山区，地跨京、津、冀、晋、鲁、豫、辽、内蒙古八省(市)区，东南部为平原区，是农业、城市的集中区域(图1a)。流域内共有234个气象站点，主要集中在华北平原(图1b)。海河流域属温带东亚季风气候区，降水多发生在汛期(6～9月)，在作物需水高峰期的3～5月降水较少，因此，海河流域水资源因农业灌溉而常年处于亏损状态。由于受人类活动和气候变化的影响，海河流域水量变化备受关注，然而，流域遭受的洪涝灾害却少有研究。当汛期发生短时强降雨，降雨强度远大于土壤水的入渗能力或水库达到饱和状态，海河流域发生洪涝灾害的风险将会增大。根据海河流域水资源公报以及水利部发布的信息，海河流域在2012年7月、2016年7月和2021年7月～10月发生了严重的洪涝灾害。
[0070]
根据本发明的一个示例性实施例：本发明采用的数据包括：grace/grace-fo twsa数据，grace卫星以前所未有的精度对tws变化进行监测，由于处理方法不同，产品可以分为球谐解和mascon solutions。与传统的球谐解相比，mascon solutions减少了陆地和海洋之间的信号泄露误差，并且无需消除南北条带误差。twsa是相对于2004～2009年平均基线的异常，通常以等效水高的形式表示，因此，mascon solutions可以直接绘制区域甚至全球的等效水高。
[0071]
本发明选择的csr rl06 mascon和jpl rl06 mascon产品分别为美国德克萨斯大学奥斯汀分校空间研究中心(csr)和nasa喷气动力实验室(jpl)提供的月值解决方案，包含时间范围2004年01月～2020年12月。通过插值近似估计两个产品中缺失的单独月份。针对grace与grace-fo间连续11个月的数据间隔，通过比较已公开的重建数据集对此空隙进行填补，重建方法包括长期降雨驱动(lpd)、贝叶斯卷积神经网络(bcnn)、贝叶斯框架(bf)。可选的，本发明将csr mascon与空间分辨率重采样至0.25
°×
0.25
°
的jpl mascon平均处理，将其作为约束数据grace/grace-fo twsa保证重建twsa的精度。
[0072]
根据本发明的一个示例性实施例：本发明的数据包括：气象数据:天尺度的降雨和温度是构建天尺度twsa的重要驱动变量。
[0073]
本发明集了5种降雨数据集分别为来自中国国家气象局的每日网格降雨数据(cma)和cn05.1，国家海洋局物理科学实验室的基于全球统一测量的日降水量分析(cpc)，全球降雨测量任务早期产品(gpm)、热带降雨测量任务(trmm)。温度数据来自cn05.1以及全球陆面同化系统-集水区地表模型(gldas-clsm)。不同于gpm、trmm以及cpc，cma和cn05.1均是基于中国2400多个气象站点数据计算生成，并且考虑了地理因素。然而，cma和cn05.1仍因处理方法不同而有所差异，前者是利用anusplin软件的薄盘样条法进行空间插值，后者则类似于xu等，数据集cn05.1使用距平逼近法(anomaly app-roach)，由气候场和距平场分
别插值后叠加得到，该数据最早包括日平均气温、最高和最低气温及降水，后增加了湿度、风速和表面蒸发三个变量。此外，考虑到trmm数据停止更新以及cma的2010、2011年数据难以获取，其他降雨数据的均值被用来代替空缺。
[0074]
根据本发明的一个示例性实施例：本发明的数据包括：辅助数据:土壤水对极端气候变化所引起的水量变化比较敏感，然而，土壤水的不同分层可能会因湿润程度表现出不同响应。gldas-clsm没有明确的土壤水分垂直水平，它以地表(0～2cm)、根区(0～100cm)和剖面(随网格变化)土壤水三部分表示。值得注意的是剖面包括曲面和根分区，根分区包括曲面。在一定程度上，剖面土壤水和tws定义上相同。为了区分土壤水和twsa，本发明取三部分土壤水之和(sma)。此外，itsg-grace2018重力场模型是由graz计算的grace重力场模型，提供卡尔曼平滑每日解决方案，包含与官方grace月度产品相同的完整水文信号。
[0075]
本发明将天尺度的sma、gldas-twsa和itsg-grace2018以及jpl-era5与重构的twsa进行比较，以评估它们在极端降雨下的水量变化，以上数据均去除了2004～2009的均值。
[0076]
根据本发明的一个示例性实施例：本发明的基于天尺度洪涝指数提高海河流域洪涝监测时效性方法，包括构建天尺度twsa。
[0077]
从twsa变化的角度而言，grace/grace-fo反演的tws受气候变化和人为活动的综合影响，降雨和温度作为基本的自然变量，可以近似表示全球或区域twsa的自然变化。
[0078]
humphrey和gudmundsson受水文建模基本原理的启发，以降雨和温度为驱动变量，通过建立一种统计模型构建气候驱动的天尺度twsa，并以grace/grace-fo twsa为约束，利用气候驱动的天尺度twsa和grace/grace-fo twsa之间的信号联系进行缩放变换处理，达到重建天尺度总twsa的目的。
[0079]
根据本发明的一个示例性实施例：本发明基本思想如下：假设一线性蓄水模型，水量输出与蓄水量和蓄水停留时间成正比。明天的水储量等于今天经温度消耗剩下的水储量加上明日的降水：
[0080]
twsa(t)＝(twsa(t-1))e-1/τ(t)
+p(t)
ꢀꢀꢀꢀꢀ
(1)
[0081]
其中，t和p(t)分别表示日时间向量和每日降雨输入，e-1/τ(t)
为水储量的衰减函数，与温度相关，取值(0，1)。
[0082]
根据本发明的一个示例性实施例：为保证重建的天尺度twsa的精度，将降雨和温度驱动的结果平均成月，并作去趋势和去周期处理，最后利用同样去趋势和去周期的grace/grace-fo twsa进行约束：
[0083]
anom(grace(tm))＝β
·
anom(twsa(tm))+ε
ꢀꢀꢀꢀ
(2)
[0084]
式中，β为约束因子(也称校正参数)，ε是残差项，tm表示月时间向量，anom表示从原始序列减去趋势项和周期项。本发明选择时间尺度为2004～2021年进行重建和校正，输入的驱动变量包括十种“温度-降雨”的组合，结果与月度grace数据以及天尺度数据进行对比，最终选择最优方案进行逐格网重建。
[0085]
根据本发明的一个示例性实施例：本发明使用时间序列分解来实现anom处理。时间序列可以分解为若干不同频率的时间序列，本发明使用时间序列分解来实现anom处理，最小二乘拟合方法可以将twsa的时间序列分解为：
[0086]
twsa
total
＝twsa
trend
+twsa
seasonal
+twsa
residual
ꢀꢀꢀꢀꢀ
(3)
[0087]
其中，twsa
total
表示原始grace/grace-fo twsa信号，twsa
trend
是线性趋势分量，twsa
seasonal
表示季节性循环分量，包括年和半年循环。twsa
residual
代表残差分量，包括原始信号去趋势和季节性循环之后的剩余量，反映了年际间的次季节信号和噪声。twsa的趋势部分较大程度反映了由人类活动导致的地下水枯竭的长期变化，而且气候驱动的年际间twsa变化通常与中低纬度地区的洪水和干旱有关。因此，本发明利用最小二乘法将残差分量和其他分量分离：
[0088]
twsa＝a+b
·
t+c
·
cos(2θt)+d
·
sin(2θt)+e
·
cos(4θt)+f
·
sin(4θt)+ε
ꢀꢀꢀ
(4)
[0089]
其中,a～f代表各分量的参数，ε与公式(2)中ε具有相同含义，对应公式(3)中的twsa
residual
，a为常数项，t为相对2004的时间向量，a+b
·
t是趋势项，对应式(3)的twsa
trend
，其余的则对应twsa
seasonal
。
[0090]
本发明将grace/grace-fo twsa分解的趋势和季节分量同样经过最小二乘法插值，得到对应的天尺度分量，直接加上残差分量，以此重建完整的天尺度twsa信号。
[0091]
根据本发明的一个示例性实施例：本发明还包括构建新型天尺度洪涝指数的步骤。
[0092]
为了实现对亚月级洪涝灾害的监测，本发明基于重建的天尺度twsa，提出了新型天尺度洪涝指数，包括fpi、wsdi和ccdi(上述指数虽名称上和现有的相同，本发明重新进行了内容定义)。
[0093]
由于grace/grace-fo twsa的趋势分量主要来自人类活动的影响，因此去趋势的总twsa被用来计算相关洪水指数。
[0094]
为了表示流域容纳和消滞洪水的潜力，降雨和twsa被用来计算fpi：
[0095]
fpa
i,j
＝p
i,j-(max(twsa
i,j
)-twsa
i,j-1
)
ꢀꢀꢀꢀꢀꢀ
(5)
[0096][0097]
其中，twsa
i,j
为去趋势后的twsa；p
i,j
为降雨；i为年(包括2004、2005、2006
……
2021)；j表示某年的第j天，从1到365/366。当fpi接近1时，意味着更高的洪水发生率。fpi强调了在陆地饱和水增加的情况下，流域向洪水易发情况转变的可能性。
[0098]
此外，另外两个基于twsa的指数wsdi和ccdi如：
[0099]
wsd
i,j
＝twsa
i,j-mean(twsaj)
ꢀꢀꢀꢀ
(7)
[0100]
pa
i,j
＝p
i,j-mean(p
i,j
)
ꢀꢀꢀꢀꢀ
(8)
[0101][0102][0103][0104]
其中，twsa和p分别代表去趋势的总twsa和降雨，pa
i,j
是降雨p相对于时序均值的异常，twsaj和paj分别为不同年份同一天的twsa和pa的集合，比如当j＝1时，则对应变量的集合包括20040101、20050101、20060101
……
20210101。为pa
i,j
的残差，wsdi和ccdi分
别是特定时间的水储存赤字和组合气候偏差。值得注意的是，wsdi、ccdi均是基于标准化思想所计算，fpi是基于归一化，此外，降雨是否考虑在内也会引起结果的差异。由于缺少天尺度洪涝指数参考，本发明并没有对得到的结果进行分级，仅把同一指数的相对大小当作是洪水的参考值，同一指数范围内，值越大，表示其发生洪水的概率越高。
[0105]
根据本发明的一个示例性实施例：本发明还包括天尺度twsa评价步骤。
[0106]
由于重建的twsa是根据grace/grace-fo twsa进行校准生成，因此相应的grace/grace-fo twsa可作为评估重建twsa效果的参考值。此外，通过比较月度重建数据来填补空窗期，因此，本发明采用以下三个常用的一致性指标进行量化评估：相关系数(cc)、均方根误差(rmse)、纳什效率系数(nse)，如：
[0107][0108][0109][0110]
其中，oi是由grace/grace-fo twsa的原始数据，mi表示重建数据，横线上标代表均值，i为第i月或日，n是总个数。cc和nse分别在[-1,1]和(-∞,1]的范围内越高，意味着重建的效果越好，rmse则相反，取值范围为[0，+∞]。
[0111]
根据本发明的一个示例性实施例：本发明还包括结果验证步骤。可选的，本发明还包括构建数据筛选步骤。grace/grace twsa作为约束数据，在重建过程中不可或缺，然而，grace与grace-fo之间的数据间隙需要被填补。图2从cc和rmse的角度显示了重建数据与原始twsa的关系，可以看到在2004年1月～2017年6月期间，所有参与比较的重建数据表现较好，cc都在0.93以上，rmse分布在20～30mm。考虑到已发布的数据在重建过程中，grace期间的twsa可能被用来获取先验信息，本发明更加重视指标在grace-fo(2018.06～)期间的表现。如图2(b)和(d)，bf代表的红色方块距离原点最近，与csr和jpl的cc和rmse分别为0.50和27.77mm、0.71和26.90mm，意味着bf与grace-fo twsa在海河流域的表现更好。因此，本发明截取基于bf重建的2017年7月～2018年5月的twsa填补grace与grace-fo间的数据间隙，以期在整个研究期内为模型提供连续的精度保障。
[0112]
本发明分别比较了降雨和温度数据，并将它们的十种组合作为输入变量，以连续的grace/grace-fo twsa作为约束条件，重建天尺度twsa，“温度-降雨”组合包括cn05.1-gpm、cn05.1-trmm、cn05.1-cpc、cn05.1-cma、cn05.1-cn05.1和gldas-gpm、gldas-trmm、gldas-cpc、gldas-cma、gldas-cn05.1。十种重建结果与jpl-rea5、itsg-grace2018和gldas-twsa比较来共同确定最优组合。
[0113]
由图3和图4显示，2004～2021年间海河流域的降雨和温度表现出明显的季节性。从整体上看，二者对于极端气候异常具有相同响应，如2012、2016、2021年的强降雨和2016年的低温。各个产品在数值上仍存在差异，具体而言，cpc的多年降雨极值表现突出，trmm和
cn05.1分别在2012年和2016年表现明显，gpm数值与其他产品也存在较大差别。各降雨产品在2011、2012、2020和2021相差较大。图3(b)和图3(c)分别展示了降雨产品之间的差异，来自中国气象局的数据cma和cn05.1以及gpm和trmm之间的拟合度较高，cc分别为0.90和0.92，rmse分别为1.58mm和1.37mm，然而，cpc与其他四个之间的拟合度总体不高。此外，对于温度，cn05.1在振幅上略低于gldas，如图4所示。
[0114]
如图5所示，对整个海河流域而言，基于不同的组合重建的月尺度结果与grace/grace-fo twsa具有相同的下降趋势，它们之间的nse在0.93～0.96之间。同样地，不同的重建天尺度twsa的次季节项表现也令人满意，天尺度重建nse结果处于0.4～0.81。通过比较十种组合在海河流域的重建效果，cn05.1-cn05.1和cldas-cn05.1重建的效果表现最好，二者的月度比较nse均达到最高值0.96，次季节项分别达到0.80和0.81。因此，基于cn05.1-cn05.1重建的天尺度twsa被用来监测海河流域洪涝灾害。此外，本发明发现降雨对模型的贡献率比温度高，因为拥有同一降雨的组合对比结果近似相等，而同一温度不同降雨的组合结果差异较大，这可能与公式(1)中衰减函数的构建密切相关。
[0115]
根据本发明的一个示例性实施例：本发明还包括洪涝监测能力验证步骤。以2016年海河水灾为例，为了验证重建效果，基于cn05.1-cn05.1重建的天尺度twsa被用来监测海河流域2016年的洪水事件。降雨过程大致为：2016年7月17日～18日降雨逐渐增大，2016年7月19日降雨大幅增强，局部多点最大降雨量超过600mm，海河水利委员会连续发布洪水预警，随后宣布发生一号洪水。图6进一步展示了海河流域2016年7月17～21日降雨的空间分布，从19日开始，强降雨从流域的西南方向东北方移动，运动范围几乎覆盖整个流域，同时降雨中心也沿着安阳-邯郸-保定-北京-天津-渤海路线转移。根据官方的水资源公报，安阳、邯郸、邢台及周围县区在此次洪水事件中受灾较为严重。因此，本发明将受灾严重的山前平原区作为洪涝中心，比较重建结果与其他天尺度产品在此次洪水事件中的水量变化。
[0116]
由图7(a)可以看出，在洪水中心区域，其他产品的水量变化并没有重建twsa强烈，为了消除基准的影响，起始点统一被设置为0。通过比较各数据7.18～7.20两天内的增长量，更直观地观察它们对极端洪水的响应。图7(b)展示了cn05.1-cn05.1在这两天的等效水高增长最高(168.83mm)，sma和gldas-twsa增长量分别为154.44mm和82.86mm，jpl-era5和itsg-grace2018较小，意味着基于cn05.1-cn05.1重建的结果能够敏感地注意到水量的增加，这对于洪涝的监测至关重要。图8展示了这五个数据对应的等效水高以及7.18-7.20两天内水量增长的空间分布，需要指出的是，sma中土壤水存在重合部分，因此sma在整个空间上表现明显大于其他数据。从空间特征可以看出，gldas-twsa与重建的天尺度twsa都能够追踪栅格化降雨运动，然而，在洪水中心的水量变化却不如重建的twsa明显。era5-jpl和itsg-grace2018在水量上同样没有做出明显响应。总之，通过与其他产品在时间和空间上的对比可初步得到结论，基于cn05.1-cn05.1重建的twsa能够敏锐地感应到极端降雨的发生，并及时做出反应，这是监测洪水事件必备的特点，证明了基于cn05.1-cn05.1重建的天尺度twsa具备监测洪涝的能力。
[0117]
根据本发明的一个示例性实施例：本发明还包括新型天尺度洪涝指数的时空分析步骤：时间特征：为量化重建天尺度twsa对洪涝灾害的影响，三个新型天尺度洪涝指数被用来从时间和空间两个方面分析海河流域发生洪水的可能性。对整个海河流域而言，有记录的洪水事件分别发生在2012年7月、2016年7月以及2021年7月～10月。因此，本发明将整个
流域作为对象，利用三次有记录的洪涝事件验证新型洪涝指数的可靠性。如图9所示，从2004年1月1日～2021年12月31日，三次有记录的洪水事件在fpi、wsdi和ccdi三个指数中表现突出，如灰色标记所示，证明了新型天尺度洪涝指数在时序上监测洪涝的准确性。相对而言，仅由twsa计算的wsdi在2016洪水事件中振幅变化较小，有降雨参与的fpi和ccdi对洪水的反应更加显著，这可能与洪水形成于降雨的累积密切相关。
[0118]
空间特征：由图10所示，fpi、wsdi和ccdi的空间分布揭示了2016年7月17日～2016年7月21日洪水由流域的西南部向东北部蔓延的过程，这与政府公布的中国水旱灾害公报的记录一致。三个指数都能够呈现洪水的大致范围，fpi、wsdi和ccdi在洪水中心的表现相同，表明它们对于极端降雨产生的洪水事件具备相同的反应。不同的是，由于fpi和ccdi中降雨的参与，对流域内洪水范围的评估要比wsdi大。
[0119]
具体而言，三个指数在7月17日和18日表现相似，仅在流域南部的小范围降雨表现明显。从19日开始，暴雨过境，有降雨参与的fpi和ccdi对流域的描绘比wsdi广，这种差异在20日较为明显。由于各指数计算方法不同，本发明通过比较fpi、wsdi和ccdi从18日～20日的变化来大致观察此次洪水的破坏范围，由(红色网格数/总网格数)
×
100％进一步计算，三个指数的影响范围分别为56％、43％和66％。根据官方资料记载，流域北部的承德并没有发生洪水，仅在城市的东部受到降雨的干扰，再者流域内辽宁西部有洪水发生的记录。因此，ccdi对此次洪水事件存在高估，意味全流域66％的地区都有可能发生洪水，尤其是承德区域的描绘，而fpi与wsdi的不同仅仅发生在北部的承德市。相比之下，wsdi既没有夸大承德的水情，同时也记录了辽河以西的洪水，比较符合公报信息。总之，基于cn05.1-cn05.1重建的天尺度twsa的新型天尺度指数都能够有效监测流域内亚月级洪涝事件，但从空间特征而言，wsdi对洪涝的破坏范围的描绘更加符合政府公开的信息。
[0120]
根据本发明的一个示例性实施例：本发明还包括不同部分的土壤水对同一洪涝事件的响应分析步骤：土壤水分的变化在一定程度上能够预示洪水信息。图8和图9已经展示了总土壤水(sma)对于暴雨的响应，这里分开讨论各部分土壤水在洪水中的变化。图11展示了在洪水中心，去趋势的各部分土壤水在洪水事件中的水量变化。从图12可以看出，在洪水到来之前，三部分土壤水在水量上并没有明显变化，随着暴雨过境，土壤水在19日和20日有较大增长，这两日剖面土壤水增长最大，变化量分别为52.46mm和31.92mm；根区土壤水次之，分别为48.18mm和23.13mm；表层土壤水变化最小，仅有2.43mm和-1.45mm。暴雨之后，各部分土壤水都呈现不同程度的下降，其中下降最快的是根区土壤水分，这意味着极端降雨所引发的洪水的主要影响范围在根区(0～100cm)。此外，土壤水的变化也体现了此次洪水事件短促的特点。
[0121]
根据本发明的一个示例性实施例：本发明还包括格网点长期变化对比步骤：本发明将重建的天尺度twsa平均成月，比较其与grace/grace-fo twsa格网点之间的差异，并计算二者间的cc和nse。如图12(a)和(b)所示，重建的月度twsa的趋势与grace/grace-fo twsa表现相似，二者都揭示了流域的南部twsa亏损的事实，其中西南部的表现尤为明显，流域内仅有北部的部分地区twsa存在较小的增长趋势，这主要是由于西南部超采地下水进行农业灌溉。如图12(c)和(d)，全流域大部分格网的cc和nse值都在0.8以上。具体来看，基于格网的cc在流域的中南部甚至接近1，nse在相同的区域也展示出较好的效果。然而，东部沿海的nse和cc仅约0.5和0.6，出现这种差异的原因可能有两点：第一是海陆交接地带的重力
信号难以被grace卫星准确捕捉，这与信号泄露具有较大关系；第二，重建方法仅仅是基于气候变化的角度，在插值方法补充剩余趋势项和季节性循环时，与grace分离的信号存在差距。但从本发明可以看出，差异的存在对于流域或者更大尺度的研究影响较小。
[0122]
本发明通过提出的天尺度洪涝指数法，克服了twsa在监测洪涝方面的缺陷，在天尺度上对亚月级洪涝进行探测。新型天尺度洪涝指数不仅能够准确探测到亚月级洪涝灾害的发生，同时由于监测的尺度更加精细，天尺度的twsa信号难以受到干扰，因此，新型天尺度洪涝指数在空间上也能有效描绘洪涝灾害的破坏范围。新型天尺度洪涝指数法充分考虑了重建的天尺度twsa感应水量变化的能力，本发明以破坏严重的2016年海河水灾为例验证新型洪涝指数法的可靠性。
[0123]
尽管已经示出和描述了本发明的实施例，对于本领域的普通技术人员而言，可以理解在不脱离本发明的原理和精神的情况下可以对这些实施例进行变化、要素组合，本发明的范围由所附权利要求及其等同物限定。

技术特征：
1.一种基于天尺度洪涝指数提高海河流域洪涝监测时效性方法，其特征在于：包括如下步骤：构建天尺度twsa：利用天尺度的降雨、温度数据及grace/grace-fo月尺度twsa数据构建天尺度twsa；构建天尺度洪涝指数：利用去趋势后的重建的天尺度twsa构建天尺度洪涝指数；天尺度洪涝指数包括fpi：fpa
i,j
＝p
i,j-(max(twsa
i,j
)-twsa
i,j-1
)其中，twsa
i,j
为去趋势后的twsa；p
i,j
为降雨；i为年；j表示某年的第j天。2.根据权利要求1所述的方法，其特征在于：天尺度洪涝指数还包括ccdi、wsdi：wsd
i,j
＝twsa
i,j-mean(twsa
j
)pa
i,j
＝p
i,j-mean(p
i,j
)))pa
i,j
是降雨p相对于时序均值的异常，twsa
j
和pa
j
分别为不同年份同一天的twsa和pa的集合；为pa
i,j
的残差，wsd
ij
和ccd
ij
分别是水储存赤字和组合气候偏差。3.根据权利要求1或2所述的方法，其特征在于：利用天尺度的降雨、温度数据及grace/grace-fo的月尺度twsa数据构建天尺度twsa包括：确定水储量与温度、降雨之间的关系：明天的水储量等于今天经温度消耗剩下的水储量加上明日的降水:twsa(t)＝(twsa(t-1))e-1/τ(t)
+p(t)其中，t和p(t)分别表示日时间向量和每日降雨输入，e-1/τ(t)
为水储量的衰减函数，与温度相关，取值(0，1)；将降雨和温度驱动的twsa结果平均成月，并作去趋势和去周期处理，最后利用同样去趋势和去周期的grace/grace-fotwsa进行约束：anom(grace(t
m
))＝β
·
anom(twsa(t
m
))+ε式中，β为约束因子，ε是残差项，t
m
表示月时间向量，anom表示从原始序列减去趋势项和周期项。4.根据权利要求1或2所述的方法，其特征在于：还包括如下步骤：将twsa的时间序列分解为：twsa
total
＝twsa
trend
+twsa
seasonal
+twsa
residual
其中，twsa
total
表示原始grace/grace-fotwsa信号，twsa
trend
是线性趋势分量，twsa
seasonal
表示季节性循环分量，包括年和半年循环；twsa
residual
代表残差分量，包括原始
信号去趋势和季节性循环之后的剩余量；去趋势后的twsa为twsa
total-twsa
trend
。5.根据权利要求4所述的方法，其特征在于：还包括如下步骤：利用最小二乘法将残差分量和其他分量分离：twsa＝a+b
·
t+c
·
cos(2θt)+d
·
sin(2θt)+e
·
cos(4θt)+f
·
sin(4θt)+ε其中,a～f代表各分量的参数，ε是残差项，对应twsa
residual
，a为常数项，t为相对2004的时间向量，a+b
·
t是趋势项，对应twsa
trend
，其余的则对应twsa
seasonal
。6.根据权利要求5所述的方法，其特征在于：还包括如下步骤：将grace/grace-fotwsa分解的趋势和季节分量同样经过最小二乘法插值，得到对应的天尺度分量，直接加上残差分量，以此重建完整的天尺度twsa信号。7.根据权利要求1或2所述的方法，其特征在于：还包括：天尺度twsa评价步骤：采用cc、rmse、nse，对构建的天尺度twsa进行评价。8.根据权利要求7所述的方法，其特征在于：天尺度twsa评价步骤包括：将重建的天尺度twsa平均成月，比较其与grace/grace-fotwsa格网点之间的差异，并计算二者之间的cc和nse。9.根据权利要求1或2所述的方法，其特征在于：还包括不同部分的土壤水对同一洪涝事件的响应分析步骤：将天尺度的sma、gldas-twsa和itsg-grace2018以及jpl-era5与重构的twsa进行比较，以评估它们在极端降雨下的水量变化，以上数据均去除了2004～2009的均值。10.根据权利要求1或2所述的方法，其特征在于：还包括结果分析步骤：采用多种温度和降雨数据，以连续的grace/grace-fotwsa作为约束条件，重建天尺度twsa，将多种重建结果与jpl-rea5、itsg-grace2018和gldas-twsa比较来共同确定最优组合。

技术总结
本发明涉及一种基于天尺度洪涝指数提高海河流域洪涝监测时效性方法，基于归一化和标准化的思想，利用重建的天尺度TWSA提出了新的天尺度洪涝指数法，包括三种新的天尺度洪涝指数：洪水潜力指数(FPI)、储水不足指数(WSDI)以及联合气候偏差指数(CCDI)。利用降雨和温度数据重建的陆地水储量异常(TWSA)与重力恢复和气候实验(GRACE)卫星反演的TWSA之间的NSE达0.96。新的天尺度洪涝指数法不仅保留了TWSA对洪水的探测能力，而且在时间尺度上更加精确。而且在时间尺度上更加精确。而且在时间尺度上更加精确。


技术研发人员：郑伟 聂圣琨 沈祎凡 祝会忠
受保护的技术使用者：辽宁工程技术大学
技术研发日：2023.02.15
技术公布日：2023/8/14