【Hydro】部分基流分割方法及程序代码说明

说明
一、数字滤波法
- 单参数数字滤波
- - Lyne-Hollick滤波法
  - Chapman滤波法
  - Chapman-Maxwell滤波法
  - Boughton-Chapman滤波法
- 双参数滤波法
- - Eckhardt滤波法
二、其他基流分割方法
- 基流指数（BFI）法
- 时间步长（HYSEP）法
- PART法
- 加里宁-阿里巴扬基流估算改进方法
三、计算程序及代码
- SWAT Baseflow Program基流计算程序
- - 基本原理
  - Baseflow输入输出文件
- Lyne-Hollick滤波法Python代码
- grwat中的基流分割方法
参考文献

说明

径流分为直接径流和基流，直接径流退水快而陡，基流退水慢而平缓。基流是指来源于地下水的径流，在枯水季节是径流的主要组成部分。基流相对稳定，对于维持河流生态、保障流域生产生活用水以及水环境保护具有重要的意义。
基流分割方法根据基流计算的原理，分为图解法、物理化学法、数学物理法、水文模型法和数值模拟法。常用的数值模拟法有数字滤波法、 BFI法、 HYSEP法、 PART法和加里宁法。

一、数字滤波法

数字滤波法利用数字滤波器将径流过程划分为直接径流和基流两个部分。直接径流代表快速响应特征，与数字信号中的高频部分相似，而基流代表慢速响应特征，与数字信号中的低频部分相似。该方法包括单参数滤波法和双参数滤波法（又称递归数字滤波法）。

单参数数字滤波

Lyne-Hollick滤波法

Lyne-Hollick滤波法由Lyne和Hollick于1979年首次提出，并由Nathan和Mcmahon于1990年引入到水文中进行基流分割，其基流分割方程为：
在这里插入图片描述

Chapman滤波法

Chapman 于1991年对Lyne-Hollick滤波法进行了改进，提出了Chapman滤波法，改进后的基流分割方程为：
在这里插入图片描述

Chapman-Maxwell滤波法

Chapman和Maxwell于1996把基流假定为同时刻的地表径流和前一时刻基流的加权平均，提出了Chapman-Maxwell滤波法，基流分割方程如下：
在这里插入图片描述

Boughton-Chapman滤波法

为使基流分割更为平滑， Boughton于1993年提出Boughton-Chapman滤波法，基流分割方程为：
在这里插入图片描述

双参数滤波法

Eckhardt滤波法

Eckhardt在单参数的基础上提出了一种双参数数字滤波方法，包含2个滤波参数（退水常数α和最大基流指数BFImax），方法更加灵活准确，其计算如式所示：
在这里插入图片描述 MRC方法可用BASEFLOW程序实现，双参数方法相对于单参数方法的优势在于可以通过调节 $BFI_{max}$ 来修改基流过程。 $BFI_{max}$ 不是一个可观测的量，Eckhardt给出了推荐值：对于常年性排水流域为0.8，对于季节性排水流域为0.5，对于常年性排水但是硬质基岩的流域为0.25。需要指出的是，当 $BFI_{max}$ 为0.5时，双参数和单参数数字滤波方法等价。

二、其他基流分割方法

基流指数（BFI）法

基流指数（BFI）法又称最小滑动法，该方法以基流指数为权系数来计算基流量，在区域枯季径流的研究中，该值可反映区域地下水资源量的参数。通常情况下，主要由地下水补给的河流BFI接近1，而季节性河流的BFI值趋近0。

BFI法是英国水文所 (Institute of Hydrology) 提出的一种基流分割方法，又叫最小滑动法。其原理为：首先将年流量过程按照给定的时间间隔N分割为m(m=365/N) 段，确定每段的最小流量序列 (q1，q2，…，qm)；然后从中确定拐点 (如果满足kqt小于qt-1和qt+1，则qt为拐点)。最后将所有的拐点连接得到基流过程，拐点之间的基流则通过线性插值得到。如果某时刻基流超过总径流，则将该时刻基流修改为总径流量。在实际分割中，通常k取值0.9或者0.979能得到较好的效果。此外，在地下水分析中基流指数BFI (一般以年或者多年为统计时期) 是反映流域特性的一个重要特征量，其定义为统计时期内基流所占总径流的比。

时间步长（HYSEP）法

Petty John和Henning于1979年提出了HYSEP程序，此方法首次实现了基流计算的程序。HYSEP法是美国地质调查局（USGS）使用的主要基流计算方法。HYSEP法有固定步长法、滑动步长法和局部最小值法3种方式。固定步长法（Fixed interval），在上述所选取的时间间隔内，将该时间间隔内的最小流量作为该时段内任意一天的基流。滑动步长法（Sliding interval），将某天前后（2N-1）/2d内的最小流量作为该天的基流，由此计算出每一天的基流。局部最小值法（Local minimum），选择时间步长内中心点前后（2N-1）/2d内的最小流量作为相邻时间步长内中心点的基流值，然后通过线性内插得到步长中心点之外时段的基流。3种方法均以本次计算的终点作为下次时间的起点，重复以上过程便可计算出基流分割过程。

PART法

PART是美国地质调查局（USGS）提出的用于分割基流的另一种计算机程序， PART法适用于以谷地泻流方式排泄或测站位于下游终点的河流，所需径流资料不少于一个水文年。该方法以日流量为数据源，基于前期径流消退进行分割。这种方法首先将日平均流量数据排列成单维数组，然后在数组中选择符合前期衰退要求的日值，在符合条件的这些天数中，如果日衰退量小于0.1个对数周期，日径流量值即作为基流量，在其余的日期中，基流值则通过线性插值的方法来获得。

加里宁-阿里巴扬基流估算改进方法

基于水量平衡原理，加里宁提出了基流的分割方法，杨远东和陈利群等对其进行改进。实际资料表明，改进的方法避免了地下水流量过程估算的任意性，且减少了计算工作量，所得成果是合理的，方法是简易可行的。其估算方程为：
在这里插入图片描述

三、计算程序及代码

SWAT Baseflow Program基流计算程序

基本原理

该程序使用递归滤波技术来分离基本流量，并计算流量衰退常数（α），SWAT官网包含该程序的下载链接。SWAT Bflow程序的算法参照Arnold和Allen等人提出的单参数数字滤波方法,其滤波方程为：
在这里插入图片描述
滤波器可以在流量数据上通过三次（正向、反向和正向），这取决于用户从流量数据的试点研究中选择的基本流量估计值。通常，每次通过都会导致基本流量占总流量的百分比减少。对于所分析的数据，显示基流的最高、最低和平均减少量。此选项为用户提供了更多的灵活性，可以调整分离，以更准确地接近现场条件。可以参考相关论文，以及实际情况选择合适的通道。

根据base flow源代码，Baseflow 中3个通道（Bflow Pass1、2、3），的基流计算代码如下：

!! perform passes to calculate base flow
f1 = 0.925
f2 = 0.
f2 = (1. + f1) / 2.
surfq(1) = strflow(1) * .5
baseq(1,1) = strflow(1) - surfq(1)
baseq(2,1) = baseq(1,1)
baseq(3,1) = baseq(1,1)
!! make the first pass (forward)    
do i = 2, isumd
  surfq(i) = f1 * surfq(i-1) + f2 * (strflow(i) - strflow(i-1))
  if (surfq(i) < 0.) surfq(i) = 0.
  baseq(1,i) = strflow(i) - surfq(i)
  if (baseq(1,i) < 0.) baseq(1,i) = 0.
  if (baseq(1,i) > strflow(i)) baseq(1,i) = strflow(i)
end do

!! make the second pass (backward)     
baseq(2,isumd - 1) = baseq(1,isumd - 1)
do i = isumd - 2, 1, -1
  surfq(i) = f1 * surfq(i+1) + f2 * (baseq(1,i) - baseq(1,i+1))
  if (surfq(i) < 0.) surfq(i) = 0.
  baseq(2,i) = baseq(1,i) - surfq(i)
  if (baseq(2,i) < 0.) baseq(2,i) = 0.
  if (baseq(2,i) > baseq(1,i)) baseq(2,i) = baseq(1,i)
end do

!! make the third pass (forward)    
baseq(3,isumd - 1) = baseq(1, isumd - 1)
do i = 2, isumd
  surfq(i) = f1 * surfq(i-1) + f2 * (baseq(2,i)- baseq(2,i-1))
  if (surfq(i) < 0.) surfq(i) = 0.
  baseq(3,i) = baseq(2,i) - surfq(i)
  if (baseq(3,i) < 0.) baseq(3,i) = 0.
  if (baseq(3,i) > baseq(2,i)) baseq(3,i) = baseq(2,i)
end do

Baseflow输入输出文件

输入文件主要是参数文件file.lst、流量序列文件xxxxxxxx.prn。

对于file.lst，文件内容如下，对于第2-4行的变量，程序不要求数字位于该行的某个位置。但是，数字必须是行上的第一个数字，后面必须有一个空格，以将其与行上写的任何其他内容隔开。
在这里插入图片描述
文件变量格式为

NDMIN，用于计算地下水衰退方程中 alpha 值的最小天数。NDMIN 的默认设置为 10。
NDMAX，用于计算地下水衰退方程中 alpha 值的最大天数。NDMAX 的默认设置为 300。
IPRINT，该变量控制每日数据的打印。0 不打印每日值，1 打印每日数值如果 IPRINT 设置为 1，则必须定义输出文件名 FLWFILE_OUT。 IPRINT 的默认设置为 0。
FLWFILE，包含每日流量数据的文件的名称。文件名最多可包含 15 个字符。可以同时处理多个的流量数据文件。
FLWFILE_OUT，每日过滤器结果写入的文件的名称。仅当 IPRINT = 1 时才需要。文件名最多可包含 15 个字符。必须为每个流数据文件指定不同的每日输出文件名。

对于xxxxxxxx.prn，
baseflow 程序以空格分隔数据，第一列为时间YYYYMMDD，第二列为流量，
除了转换制表符之外，还必须删除 USGS 放置在数据文件顶部的所有标题信息。基流程序假定文件中的第一行用于标头信息。所有其他行应包含日期和每日流量值。

输出文件主要为baseflow.dat，xxxxxxxx.out（当IPRINT为1时），baseflow.dat各指标为
Baseflow Fr1，由第一通道估计的基流贡献系数。
Baseflow Fr2 由第二通道估计的基流贡献系数。
Baseflow Fr3 由第三通道估计的基流贡献系数。
NPR 用于计算主衰退曲线的单个基流衰退数量。
Alpha Factor 基流衰退常数
Baseflow days 基流衰退在一个对数周期内下降的天数。

xxxxxxxx.out中各列为：
Streamflow 流量数据文件中的逐日流量， $m^3/s$ ；
Bflow Pass1 第一通道估计的逐日基流， $m^3/s$ ；
Bflow Pass2 第二通道估计的逐日基流， $m^3/s$ ；
Bflow Pass3 第三通道估计的逐日基流， $m^3/s$ 。

Lyne-Hollick滤波法Python代码

#! /usr/bin/env python
# -*- coding: utf-8 -*-
import numpy as np

def hollickLyneFilter(data, beta=0.925, invert=True):
    """
    Calculate baseflow from a given discharge timeseries using the 
    digital filter method proposed by Hollick and Lyne (1979).


    Definition
    ----------
    def hollickLyneFilter(data,beta=0.925,invert=True)

    Input
    -----
    data   : numpy.ndarray 1D

    Optional Input
    --------------
    beta   : float -> dimensionless filter parameter. The default value 
                      was proposed by Nathan and McMahon (1990)
    invert : True  -> add a backward filter pass
           : False -> no backward filterpass, will result in a phase shift

    Output
    ------
    numpy.ndarray 1D

    Literature
    ----------
    Lyne, V. & Hollick, M. 1979, “Stochastic time-variable rainfall-runoff modelling”, 
    Proceedings of the Hydrology and Water Resources Symposium, Perth, 10-12 September, 
    Institution of Engineers National Conference Publication, No. 79/10, pp. 89-92. 

    Nathan, R. J. & McMahon, T. A. 1990a, “Evaluation of automated techniques for base 
    flow and recession analysis”, Water Resources Research , Vol. 26, pp. 1465-1473. 
    """
    out = data.copy()

    for i in range(1, len(data)):
        bflow = beta * out[i-1] + (1 - beta) * .5 * (data[i]+data[i-1])
        out[i] = min(bflow, data[i])

    if invert:
        out = hollickLyneFilter(out[::-1], beta, invert=False)[::-1]

    return out[1:]


if __name__ == "__main__":
    data = np.array([368.500, 408.100, 493.600, 502.500, 453.100, 439.600, 408.100,
                     430.600, 417.100, 412.600, 390.100, 381.300, 360.300, 340.000,
                     293.100, 276.600, 283.100, 279.800, 270.200, 263.900, 239.800])
    bflow = hollickLyneFilter(data)
    print("", bflow)
    blow = hollickLyneFilter(data, invert=False)
    print("", bflow)
    blow = hollickLyneFilter(data, .8)
    print("", bflow)

grwat中的基流分割方法

grwat实现了几种基流滤波方法，包括Lyne和Hollick（1979）、Chapman（1991）、Boughton（1993）、Jakeman和Hornberger（1993）以及Chapman and Maxwell（1996）的方法。

hdata = spas %>%
  mutate(lynehollick = gr_baseflow(Q, method = 'lynehollick', a = 0.9),
         boughton = gr_baseflow(Q, method = 'boughton', k = 0.9),
         jakeman = gr_baseflow(Q, method = 'jakeman', k = 0.9),
         maxwell = gr_baseflow(Q, method = 'maxwell', k = 0.9)) %>% 
  pivot_longer(lynehollick:maxwell, names_to = 'Method', values_to = 'Qbase')

ggplot(hdata) +
  geom_area(aes(Date, Q), fill = 'steelblue', color = 'black') +
  geom_area(aes(Date, Qbase), fill = 'orangered', color = 'black') +
  scale_x_date(limits = c(ymd(19810101), ymd(19811231))) +
  facet_wrap(~Method)
#> Warning: Removed 93508 rows containing non-finite values (`stat_align()`).
#> Removed 93508 rows containing non-finite values (`stat_align()`).

在这里插入图片描述