CN107084709A - 一种多弹性对径流变化驱动因素的定量分割方法 - Google Patents

Abstract

本发明公开了一种多弹性对径流变化驱动因素的定量分割方法，通过收集流域内长序列水文气象资料，计算潜在蒸发量和实际蒸发量，再采用单参数Budyko方程计算代表流域特性的参数，分析流域内各水文气象变量相对于多年平均值的变化程度，计算径流对七个驱动因素的弹性，判断年径流过程的突变点，突变点之前为基准期，之后为变化期，利用各水文气象变量的弹性系数以及在变化期相对于基准期的差异，分析七个驱动因素对径流变化的贡献率。本发明结合流域内水文气象要素变化，将气候弹性法推广到径流变化对降水，最高、最低气温，风速，日照，相对湿度和流域特性的多弹性分析，定量辨识各驱动因素对径流变化的影响。该方法适用性广，计算成果更全面，更科学。

Description

一种多弹性对径流变化驱动因素的定量分割方法

技术领域

本发明涉及一种多弹性对径流变化驱动因素的定量分割方法，属于水文水资源应用技术领域。

背景技术

径流是地表水文循环的重要环节，作为人类最早可利用的水资源形式之一，其变化过程很早以前就受到人们的关注。近年来，以气候变化和强人类活动影响为特征的环境变化导致中国乃至全球很多流域的河川径流都发生了显著的变化。特别是对于我国而言，在社会经济迅猛发展和水资源短缺、洪旱灾害问题并存的背景下，中国北方的海河、黄河和辽河等河流的径流量的显著下降对我国水资源可持续利用和管理提出了严峻的挑战。因此，精确解读气候变化和强人类活动下水资源转化过程的响应规律和机理，定量辨识气候变化和人类活动的各因子对地表径流的影响，是减少水资源评估不确定性，从而保障区域水资源安全的关键科学基础。

如何准确的分离出气候变化的各变量和人类活动对地表径流的影响十分具有挑战性，目前在气候变化水文响应评估中被广泛采用的方法主要有两类：基于气候弹性系数的方法和基于水文模拟的方法。水文模拟方法采用水文模型来量化气候变化及人类活动引起的下垫面改变对径流的影响。该方法具有较好的物理基础，但模型结构和参数等存在一定的不确定性，可能导致评价结果产生较大偏差。另一种基于气候弹性系数的方法，已有研究证明该方法是一种简单有效的方法。随着流域水文技术的发展，气候弹性法已由最初的只考虑降水对径流影响的单参数弹性法发展到考虑降水，温度，辐射，风速以及相对湿度的多弹性法。然而，辐射数据测站少且不易获取，辐射与温度和湿度之间的相关性不可忽略。如何消除变量间的相关性影响，更精确地定量评估相对独立的各气候因子对径流变化的贡献？此外，如何利用弹性法计算人类活动对径流变化的影响？这些都是需要解决的关键问题。

发明内容

本发明所要解决的技术问题是：提供一种多弹性对径流变化驱动因素的定量分割方法，可准确定量地评估降水、最高气温、最低气温、风速、日照、相对湿度和流域特性对径流变化的影响。

本发明为解决上述技术问题采用以下技术方案：

一种多弹性对径流变化驱动因素的定量分割方法，包括如下步骤：

步骤1，收集流域内至少40年的年径流量资料以及日气象要素资料，其中，气象要素资料包括降水量、最高气温、最低气温、平均相对湿度、平均风速、日照时数；

步骤2，结合最高气温、最低气温、平均相对湿度、平均风速、日照时数，运用Penman-Montieth公式计算流域日潜在蒸发量，进而得到流域多年平均潜在蒸发量；

步骤3，根据流域内的水量平衡方程，结合流域多年平均径流量与多年平均降水量，计算流域多年平均实际蒸发量；

步骤4，采用单参数的Budyko方程，结合多年平均降水量、多年平均潜在蒸发量、多年平均实际蒸发量，计算多年平均的流域特性参数；

步骤5，计算径流对降水量、最高气温、最低气温、平均相对湿度、平均风速、日照时数以及流域特性这七个驱动因素的弹性系数；

步骤6，分析流域内年径流量的变化趋势，判断年径流过程的突变点，将突变点之前作为基准期，突变点之后作为变化期；

步骤7，根据弹性系数、多年平均径流量、各个驱动因素的多年平均值、变化期径流量均值与基准期径流量均值的差值、各个驱动因素变化期均值与基准期均值的差值，计算各个驱动因素对径流变化的贡献率。

作为本发明的一种优选方案，步骤2所述日潜在蒸发量的计算公式为：

其中，E_{0_d}为日潜在蒸发量，Δ为饱和水汽压斜率，R_n为净辐射，G为土壤热通量，γ为干湿常数，T为平均温度，U为平均风速，VPD为饱和水汽压差。

作为本发明的一种优选方案，步骤3所述多年平均实际蒸发量的计算公式为：

E＝P-Q

其中，E为多年平均实际蒸发量，P为多年平均降水量，Q为多年平均径流量。

作为本发明的一种优选方案，步骤4所述多年平均的流域特性参数的计算公式为：

E＝P/[1+(P/E₀)ⁿ]^1/n

其中，E为多年平均实际蒸发量，P为多年平均降水量，E₀为多年平均潜在蒸发量，n为多年平均的流域特性参数。

作为本发明的一种优选方案，步骤6所述判断年径流过程的突变点方法为：Pettitt检验方法。

作为本发明的一种优选方案，步骤7所述各个驱动因素对径流变化的贡献率计算公式为：

其中，P，TA，TI，U，N，RH，n分别为多年平均降水量、多年平均最高气温、多年平均最低气温、多年平均风速、多年平均日照时数、多年平均相对湿度、多年平均的流域特性参数；η_X为贡献率，ε_X为弹性系数，Q为多年平均径流量，ΔQ为变化期径流量均值与基准期径流量均值的差值，ΔX为各个驱动因素变化期均值与基准期均值的差值。

本发明采用以上技术方案与现有技术相比，具有以下技术效果：

1、本发明基于Budyko方程和Penman-Montieth公式，结合流域内的水文气象要素变化，将气候弹性法推广到径流变化对降水，最高、最低气温，风速，日照，相对湿度和流域特性的多弹性分析，从而定量辨识各驱动因素对径流变化的影响。

2、本发明方法具有很好的适用性，所需资料极易获取，大大降低了计算成本，计算成果也更全面、更科学，既消除了各气象变量间的相关性，考虑了各气象因子对径流的影响，也利用表征流域特性的参数考虑了人类活动对径流的影响。

附图说明

图1是本发明多弹性对径流变化驱动因素的定量分割方法的流程框图。

具体实施方式

下面详细描述本发明的实施方式，所述实施方式的示例在附图中示出。下面通过参考附图描述的实施方式是示例性的，仅用于解释本发明，而不能解释为对本发明的限制。

如图1所示，为本发明一种多弹性对径流变化驱动因素的定量分割方法的流程框图，包括以下步骤：

(1)基础资料收集：收集流域内的长序列(一般为40年以上)水文气象资料(降水量、最高气温、最低气温、平均相对湿度、平均风速、日照时数均以日为单位；径流量以年为单位)。

(2)结合日最高气温、最低气温、平均相对湿度、平均风速、日照时数，运用Penman-Montieth公式计算流域长序列的日潜在蒸发量，进而可求得流域多年平均潜在蒸发量。

日潜在蒸发量E_{0_d}：

其中，Δ为饱和水汽压斜率(kPa℃-1)，R_n为净辐射(MJ m-2d-1)，G为土壤热通量(MJ m-2d-1)，γ为干湿常数(kPa℃-1)，T为平均温度(℃)，U为平均风速，VPD为饱和水汽压差，e_s为平均饱和水汽压(kPa)，e_a为实际水汽压(kPa)。

多年平均潜在蒸发量E₀：

其中，M为序列的年份，m为一年中的天数。多年平均降水量、多年平均最高气温、多年平均最低气温、多年平均风速、多年平均日照时数、多年平均相对湿度均可参照公式(2)。

(3)运用流域长期平均的水量平衡方程(即P＝E+Q)结合流域多年平均径流量与降水量来计算流域多年平均实际蒸发量E：

E＝P-Q (3)

多年平均径流量由各年的年径流量总和除以长序列的长度得到，长序列以40年为例，则多年平均径流量由1到40年各年的年径流量总和除以40得到。

(4)选择一种形式的单参数Budyko方程，目前常用的有四种形式，包括：

Fu：

Choudhury：E＝P/[1+(P/E₀)ⁿ]^1/n

Zhang：E＝P[1+n(E₀/P)]/[1+n(E₀/P)+P/E₀]

Wang-Tang：

以上各种形式的Budyko方程均可表示为E＝f(E₀,P,n)，根据流域多年平均实际蒸发量，降水量和潜在蒸发量来计算出代表流域特性的参数值。如选用Choudhury的Budyko方程，计算代表流域特性的参数值n。

(5)分析流域内各水文气象变量相对于多年平均值的变化程度，计算径流对降水、日最高气温、最低气温、平均相对湿度、平均风速、日照时数以及流域特性共七个驱动因子的弹性。

径流对各驱动因素的弹性定义为单位气候要素或流域特性参数的变化导致的流域径流量的变化程度，如降水量增加1％导致流域年径流量相对于多年平均值的变化百分比，其中，温度的弹性习惯上定义为温度增加1℃导致流域年径流量相对于多年平均值的变化百分比。其具体推导过程如下：

根据流域长期平均的水量平衡方程，多年平均降水量P的变化可以表示为以下全微分形式：

dP＝dE+dQ (4)

同样地，根据Budyko方程，多年平均实际蒸发量E的变化可以表示为以下全微分形式：

将公式(5)带入公式(4)可得以下形式：

等式两边同时除以Q＝P-E，年径流量的相对变化可表示为：

令和上式可化简为：

由定义可知，ε_P、和ε_n分别为径流对降水、潜在蒸发和流域特性的弹性。以上P、E₀、Q均表示各变量的多年平均值，n表示多年平均的流域特性参数。

相似地，根据Penman-Montieth公式，日潜在蒸发量的变化可表示为以下全微分形式：

公式(9)是建立在日尺度基础上的，因为Penman-Montieth公式是利用日尺度的气象变量来计算的日潜在蒸发量，但由于弹性的定义是建立在多年平均尺度上的，所以为了更清晰地表达尺度转换问题，这里引入下标_d表示日变量，_y表示年变量。则公式(9)可表示为：

则年潜在蒸发量的变化可表示为：

式中，m为一年中的天数，进一步地，年潜在蒸发量的相对变化可表示为：

假设在长序列(一般大于40年)的估算中，以下等式能成立：

式中，表示各变量的多年平均值，将公式(13)带入公式(12)，可得以下形式：

为了与公式(8)统一表达形式，公式(14)可进一步整理化简为：

式中，ε_TA′、ε_TI′、ε_U′、ε_N′和ε_RH′分别表示潜在蒸发对最高温度、最低温度、风速、日照和相对湿度的弹性。将公式(15)带入公式(8)，可得以下形式：

令和则上式可写为：

式中，ε_P、ε_TA、ε_TI、ε_U、ε_N、ε_RH和ε_n分别为径流对降水、最高温度、最低温度、风速、日照、相对湿度和流域特性的弹性。

(6)分析流域内径流量的历史演变趋势，判断年径流过程的突变点，突变点之前为基准期，突变点之后为变化期。运用Pettitt检验方法检验年径流系列的突变点，具体步骤如下：

①使用统计量U_t,M来检验同一个总体x(t)的两个样本，统计量U_t,M的公式：

式中，t表示年，M为序列长度，本实施例中为40年，若x_t-x_i>0，则sgn(x_t-x_i)＝1；若x_t-x_i＝0，则sgn(x_t-x_i)＝0；若x_t-x_i<0，则sgn(x_t-x_i)＝-1。

②Pettitt检验的零假设为序列无变异点，其统计量K_M和相关概率的显著性检验公式如下：

K_t,M＝max|U_t,M|,(1≤t≤M)

式中，若p≤0.05，则认为t点为显著变异点。由此将年径流系列分为基准期和变化期，t点包括在变化期中。

(7)利用径流对降水、日最高气温、最低气温、平均相对湿度、平均风速、日照时数以及流域特性共七个驱动因子的弹性，全时段径流量与七个驱动因子的均值，以及变化期径流量和七个驱动因子与基准期的差值，分析七个驱动因素对径流变化的贡献率。

变化期相对于基准期径流量的变化可以表示为七个驱动因子的变化导致的径流变化之和，即：

计算各驱动因子对径流变化的贡献率，计算公式如下：

式中，X表示各变量的多年平均值，ΔX表示各变量在变化期均值相对于基准期均值的差值。

以上实施例仅为说明本发明的技术思想，不能以此限定本发明的保护范围，凡是按照本发明提出的技术思想，在技术方案基础上所做的任何改动，均落入本发明保护范围之内。

Claims (6)

1.一种多弹性对径流变化驱动因素的定量分割方法，其特征在于，包括如下步骤：

步骤1，收集流域内至少40年的年径流量资料以及日气象要素资料，其中，气象要素资料包括降水量、最高气温、最低气温、平均相对湿度、平均风速、日照时数；

步骤2，结合最高气温、最低气温、平均相对湿度、平均风速、日照时数，运用Penman-Montieth公式计算流域日潜在蒸发量，进而得到流域多年平均潜在蒸发量；

步骤3，根据流域内的水量平衡方程，结合流域多年平均径流量与多年平均降水量，计算流域多年平均实际蒸发量；

步骤4，采用单参数的Budyko方程，结合多年平均降水量、多年平均潜在蒸发量、多年平均实际蒸发量，计算多年平均的流域特性参数；

步骤5，计算径流对降水量、最高气温、最低气温、平均相对湿度、平均风速、日照时数以及流域特性这七个驱动因素的弹性系数；

步骤6，分析流域内年径流量的变化趋势，判断年径流过程的突变点，将突变点之前作为基准期，突变点之后作为变化期；

步骤7，根据弹性系数、多年平均径流量、各个驱动因素的多年平均值、变化期径流量均值与基准期径流量均值的差值、各个驱动因素变化期均值与基准期均值的差值，计算各个驱动因素对径流变化的贡献率。

2.根据权利要求1所述多弹性对径流变化驱动因素的定量分割方法，其特征在于，步骤2所述日潜在蒸发量的计算公式为：

<mrow> <msub> <mi>E</mi> <mrow> <mn>0</mn> <mo>_</mo> <mi>d</mi> </mrow> </msub> <mo>=</mo> <mfrac> <mrow> <mn>0.408</mn> <mi>&amp;Delta;</mi> <mrow> <mo>(</mo> <msub> <mi>R</mi> <mi>n</mi> </msub> <mo>-</mo> <mi>G</mi> <mo>)</mo> </mrow> <mo>+</mo> <mi>&amp;gamma;</mi> <mrow> <mo>(</mo> <mn>900</mn> <mo>/</mo> <mo>(</mo> <mi>T</mi> <mo>+</mo> <mn>273</mn> <mo>)</mo> </mrow> <mo>)</mo> <mi>U</mi> <mi>V</mi> <mi>P</mi> <mi>D</mi> </mrow> <mrow> <mi>&amp;Delta;</mi> <mo>+</mo> <mi>&amp;gamma;</mi> <mrow> <mo>(</mo> <mn>1</mn> <mo>+</mo> <mn>0.34</mn> <mi>U</mi> <mo>)</mo> </mrow> </mrow> </mfrac> </mrow>

其中，E_{0_d}为日潜在蒸发量，Δ为饱和水汽压斜率，R_n为净辐射，G为土壤热通量，γ为干湿常数，T为平均温度，U为平均风速，VPD为饱和水汽压差。

3.根据权利要求1所述多弹性对径流变化驱动因素的定量分割方法，其特征在于，步骤3所述多年平均实际蒸发量的计算公式为：

E＝P-Q

其中，E为多年平均实际蒸发量，P为多年平均降水量，Q为多年平均径流量。

4.根据权利要求1所述多弹性对径流变化驱动因素的定量分割方法，其特征在于，步骤4所述多年平均的流域特性参数的计算公式为：

E＝P/[1+(P/E₀)ⁿ]^1/n

其中，E为多年平均实际蒸发量，P为多年平均降水量，E₀为多年平均潜在蒸发量，n为多年平均的流域特性参数。

5.根据权利要求1所述多弹性对径流变化驱动因素的定量分割方法，其特征在于，步骤6所述判断年径流过程的突变点方法为：Pettitt检验方法。

6.根据权利要求1所述多弹性对径流变化驱动因素的定量分割方法，其特征在于，步骤7所述各个驱动因素对径流变化的贡献率计算公式为：

<mfenced open = "" close = ""> <mtable> <mtr> <mtd> <mrow> <msub> <mi>&amp;eta;</mi> <mi>X</mi> </msub> <mo>=</mo> <msub> <mi>&amp;epsiv;</mi> <mi>x</mi> </msub> <mfrac> <mi>Q</mi> <mi>X</mi> </mfrac> <mfrac> <mrow> <mi>&amp;Delta;</mi> <mi>X</mi> </mrow> <mrow> <mi>&amp;Delta;</mi> <mi>Q</mi> </mrow> </mfrac> </mrow> </mtd> <mtd> <mrow> <mi>X</mi> <mo>=</mo> <mi>P</mi> <mo>,</mo> <mi>T</mi> <mi>A</mi> <mo>,</mo> <mi>T</mi> <mi>I</mi> <mo>,</mo> <mi>U</mi> <mo>,</mo> <mi>N</mi> <mo>,</mo> <mi>R</mi> <mi>H</mi> <mo>,</mo> <mi>n</mi> </mrow> </mtd> </mtr> </mtable> </mfenced>

其中，P，TA，TI，U，N，RH，n分别为多年平均降水量、多年平均最高气温、多年平均最低气温、多年平均风速、多年平均日照时数、多年平均相对湿度、多年平均的流域特性参数；η_X为贡献率，ε_X为弹性系数，Q为多年平均径流量，ΔQ为变化期径流量均值与基准期径流量均值的差值，ΔX为各个驱动因素变化期均值与基准期均值的差值。