🧄Python时间序列模型使用河流-地下水-泵水数据计算地下水位下降

关键词

Python | 时间序列 | 地下水 | 水文 | 泵水 | 河流水位 | 传递函数噪声建模 | 物理脉冲响应函数 | 降水过剩 | 河流水位变化 | 自回归噪音模型 | Makkink参考蒸发 | Hantush井函数 | 贝塞尔函数 | 赤池信息准则 | 畜水层 | 稳态回撤

梗概

方法

时间序列模型表示对观察井中不同应力对水头的作用的独立估计， 应用多模型方法来确定哪些水文应力与描述观测井中的水头动力学相关。

降水过剩、河流水位和地下水抽取被列为潜在相关的水文压力。 每个头部时间序列测试了八种不同的模型结构。 最简单的模型仅考虑降水过剩，由降水和潜在蒸发量计算得出。 下一个模型将河流阶段添加为压力。 在接下来的三个模型中，最多添加三个井场作为潜在应力，从最近的井场开始向最远的井场移动。 最后三个模型重复这最后一步，忽略河流作为压力。

最好的模型结构是从每个观测井的可靠模型集中选出的。 拆分样本测试，其中一部分时间序列保持独立，用于测试校准模型。

时间序列模型

时间序列建模方法，也称为传递函数噪声建模，使用基于物理的脉冲响应函数来描述头部对不同压力的响应。 降水过剩、河流水位变化和噪声建模的影响模拟基于 Von Asmuth 等人的标准方法。 时间序列模型写为

$h(t)=\sum_{m=1}^M h_m(t)+d+r(t)$

模型校准

数据集准备

import json
import logging
import sys

sys.path.insert(1, "..")
from util import get_pastastore

logging.basicConfig(level=logging.INFO)

with open("../tsa_settings.json", "r") as f:
    settings = json.load(f)

name = f"frontiers_{settings['freq']}"
pstore = get_pastastore(name, settings["connector_type"])
conn = pstore.conn

prec = hpd.ObsCollection.from_knmi(xmid=oseries.x, ymid=oseries.y,
                                   meteo_vars=["RH"], start=["1980-01-01"],
                                   ObsClass=[hpd.PrecipitationObs])

evap = hpd.ObsCollection.from_knmi(xmid=oseries.x, ymid=oseries.y,
                                   meteo_vars=["EV24"], start=["1980-01-01"],
                                   ObsClass=[hpd.EvaporationObs])

# %% write to store

for o in tqdm(prec.loc[:, "obs"].values, desc="Add precipitation series"):

    meta = _get_metadata_from_obs(o)
    # convert to mm/day
    conn.add_stress(o["RH"] * 1e3, o.name, "prec",
                    metadata=meta, overwrite=True)

for o in tqdm(evap.loc[:, "obs"].values, desc="Add evaporation series"):

    meta = _get_metadata_from_obs(o)
    # convert to mm/day
    conn.add_stress(o["EV24"] * 1e3, o.name, "evap",
                    metadata=meta, overwrite=True)

读取数据

建模和模型选择

数据结果可视化

源代码

🏈page指点迷津 | Brief