本发明公开了一种基于改进Budyko公式的径流变化定量分解方法，主要步骤如下：首先获取研究流域的主要水文变量并估算该流域的潜在蒸散发序列；其次采用闭合流域水量平衡公式获取研究流域的陆地水储量变化序列，并将其耦合加入到原始Budyko公式中进而获得改进的Budyko公式；然后依据改进的Budyko公式构建研究流域在水热耦合平衡状态下的数学模型；最后基于样本方差定义进一步量化不同变量对径流变化的贡献作用，并对研究流域内径流变化的主要驱动因子进行分析和归纳总结。本发明方法不仅考虑了陆地水储量变化对流域水文循环的影响作用，还可以对研究流域内的径流变化进行归因量化分析，为变化环境下的径流预测和极端洪涝灾害事件预报提供了新的理论支撑和技术支持。
一种基于改进Budyko公式的径流变化定量分解方法
技术领域
[0001] 本发明属于水文水资源领域的径流变化归因分析领域，特别是一种基于改进Budyko公式的径流变化定量分解方法。
背景技术
[0002] 水资源是人类社会发展和工业进步最根本、最必要的要素。径流是水资源在地球上的最重要的存在形式，也是社会发展和工业进步所需水资源的主要来源。近年来随着气候变化日益显著以及人类活动的逐渐增加，全球不同研究区域内的水文循环过程发生了显著改变，加剧了不同流域间水资源时空间分布严重不均的局面，进而造成更多的极端干旱和洪涝事件，严重威胁区域人类活动的水资源安全。因此，深入分析并探讨变化环境下引起径流变化的主要原因及背后的水文机理，不仅可以帮助水资源管理者制定更为合理的水资源管理措施和调度方案，同时还可以有效预测未来时期内极端洪涝事件发生的可能性，进而保障人民群众的生命财产安全。
[0003] Budyko假设最早由Budyko等人于1974年提出，该假设认为研究流域内的多年平均蒸散发主要取决于大气对陆面的水分供给和大气蒸散发需求量之间的平衡关系。此后，傅抱璞等人根据研究流域内水文气象的物理机制对原始的Budyko假设进行了补充，并得到一组Budyko假设对应的解析值结果，不仅为Budyko假设赋予了相应的物理意义，也为各研究学者提供了新的理解角度，大大加深了该假设在水文学研究中的应用。
[0004] Budyko公式具有计算简单、准确度高以及可靠性强等优点，能够合理描述特定研究流域在水热耦合平衡状态下的主要水文变量之间相关关系，故国内外学者常常基于该公式对变化环境下导致研究流域内径流变化的主要驱动因子进行归纳分析。然而，该假设往往忽略了陆地水储量对流域径流变化的影响作用，且主要考虑长时间尺度下的流域状态变化，从而导致得到的分析结果与实际结果之间可能产生较大差异，尤其是在陆地水储量显著变化的区域内上述差异更加明显。
发明内容
[0005] 本发明同时考虑特定闭合流域的水热耦合平衡条件以及陆地水储量变化对水文循环过程的影响作用，提出一种基于改进Budyko公式的径流变化定量分解方法。该方法可有效避免传统的Budyko公式未能考虑陆地水储量对流域水文循环过程影响作用，以及未能有效量化不同要素对径流变化贡献作用的不足。
[0006] 本发明采用以下技术方案实现：
[0007] 一种基于改进Budyko公式的径流变化定量分解方法，包括以下步骤：
[0008] Step 1：获取研究流域的主要水文变量，包括研究流域的降水量、径流量、最高温度、最低温度、风速、相对湿度、日照时长、蒸散发等，并基于Penman‑Monteith公式计算该流域的潜在蒸散发；
[0009] Step 2：为有效考虑陆地水储量变化对流域水文循环的影响作用，通过闭合流域水量平衡公式估算该流域在研究时间段内的陆地水储量变化值；
[0010] Step 3：将Step 2中所获得的陆地水储量变化值耦合加入到原始Budyko公式中，通过构建改进的Budyko公式来表征研究流域水文变量在水热耦合平衡状态下的交互作用关系；最后，通过最小二乘法求取改进的Budyko公式中所包含的特征参数；改进的Budyko公式由于能够充分考虑陆地水储量变化对流域水文循环过程的影响作用，故而能够更加准确地描述短时间尺度下动态水文变化过程中各水文变量之间的交互作用关系；
[0011] Step 4：对Step 3中所获得的改进的Budyko公式采用皮亚诺型余项的泰勒公式展开，并依据样本方差的定义，对研究流域的径流变化部分进行定量分解。最后，根据分解结果求取不同水文变量对研究流域径流变化的贡献作用。
[0012] 上述技术方案中，进一步地，所述Step 1具体为：
[0013] Step 1‑1，确定位于研究流域的气象站点，收集并整理计算月尺度潜在蒸散发所需要的水文变量，包括月尺度的最高温度、最低温度、风速、相对湿度、日照时长；
[0014] Step 1‑2，基于Penman‑Monteith公式估算研究流域对应的月尺度潜在蒸散发序列，其中，所述的Penman‑Monteith公式为：
[0015]
[0016] 式中：PET为研究流域的潜在蒸散发，mm·month‑1；Δ为饱和水汽压‑温度曲线斜‑1 ‑2 ‑1 ‑1率，kPa·℃ ；n为净辐射，MJ·m ·month ；γ为湿度计常数，kPa·℃ ；G为地表热通量，‑2 ‑1
MJ·m ·month ；T为距离地面2米处的空气平均温度，℃；u为距离地面2米处的平均风速，‑1
m·s ；es为距离地面2米处的饱和水汽压，kPa；ea为距离地面2米的实测水汽压，kPa；
[0017] Step 1‑3，将对应的水文变量代入式(1)进行计算，即可得到不同站点的月尺度潜在蒸散发；然后，基于不同气象站点的潜在蒸散发估算结果，通过反距离插值方法得到整个研究流域的平均潜在蒸散发。
[0018] 进一步地，所述Step 1‑3具体包括以下步骤：
[0019] Step 1‑3‑1，根据研究流域的面积大小，将其空间划分成若干个网格，分别计算所有网格格点到各个气象站点的距离：
[0020]
[0021] 式中：(xi,yi)为研究流域内对应第i个格点的经纬度；(Xj,Yj)为对应第j个气象站点的经纬度；Dij即为第i个格点到第j个气象站点的距离；
[0022] Step 1‑3‑2，进一步求取各个气象站点对研究流域内第i个格点的影响权重系数，即为：
[0023]
[0024] 式中：wij为第j个气象站点对第i个格点的影响权重系数；m为气象站点的总个数；
[0025] Step 1‑3‑3，最后，求取不同水文变量在流域内各个格点的对应数值，进而求取该水文变量在研究流域内的平均值：
[0026]
[0027] 式中：fj为水文变量在第i个格点的对应数值；fj为第j个气象站点对应的气象数据观测值。
[0028] 进一步地，所述Step 2具体包括以下步骤：
[0029] Step 2‑1，将控制水文站的月尺度径流量转化为月尺度径流深，其计算方法为：
[0030]
[0031] 式中：R为研究流域气象站点的月径流深，单位为mm；R′为研究流域控制水文站的3 2
月径流量，m；A为研究流域的总面积，m；
[0032] Step 2‑2，基于各个气象站点的降水量观测数据，通过泰森多边形法求取研究流域的月尺度面平均降水量：
[0033]
[0034] 式中：P为整个研究流域的月尺度面平均降水量，mm；wj为第j个气象站点占整个研究流域的权重值；Pj为第j个气象站点的降水量观测值，mm；m为气象站点的总个数；
[0035] Step 2‑3，联合利用上述各式所求得的月尺度面平均降水量、月径流深以及蒸散发序列，基于闭合流域水量平衡公式，求取研究流域的月尺度陆地水储量变化值：
[0036] ΔS＝P‑R‑ET (7)
[0037] 式中：ΔS为研究流域在对应时间段内的陆地水储量变化值，mm；ET为研究流域实测蒸散发，mm。
[0038] 进一步地，所述Step 3具体包括以下步骤：
[0039] Step 3‑1，基于原始Budyko公式求取多年尺度(一般规定为10年以上)下水文变量在水热耦合平衡状态下的交互作用关系，并基于最小二乘法求取该公式所包含的特征参数[0040]
[0041] 式中： 为该流域的蒸散发多年平均值，mm； 和 分别为该流域的潜在蒸散发多年平均值及降水量多年平均值，其单位均为mm；为干旱指数多年平均值，是一个无量纲参数；为该流域水热耦合平衡状态下对应的特征参数，是一个无量纲参数；
[0042] Step 3‑2，当考虑其他更短的时间尺度(如月尺度、季尺度或者年尺度等)时，则需通过下式求取研究流域在水热耦合平衡状态对应的水文变量交互作用关系：
[0043]
[0044] 式中：ETi为该流域在短时间尺度下的蒸散发值，mm；PETi和Pi分别为该流域短时间尺度下的潜在蒸散发值及降水量，其单位均为mm； 为短时间尺度下的干旱指数，是一个无量纲参数；n为该流域水热耦合平衡状态下对应的特征参数，同样为一个无量纲参数；
[0045] Step 3‑3，考虑陆地水储量变化对研究流域水文过程的影响作用，需按照下式求取研究流域的有效降水量：
[0046] P′i＝Pi‑ΔSi (10)[0047] 式中：P′i为研究流域有效降水量，mm；Pi为基于气象站得到的研究流域实测降水量，mm；ΔSi为在对应时间段内的研究流域陆地水储量变化值，mm；
[0048] Step 3‑4，用有效降水量P′i替代式(9)中的实测降水量Pi，则可获得考虑陆地水储量变化影响作用的改进的Budyko公式，并基于最小二乘法求取该公式所包含的特征参数n：
[0049]
[0050] 式中：P′i为研究流域的有效降水量，mm； 为短时间尺度下考虑陆地水储量变化影响作用的干旱指数，是一个无量纲参数。
[0051] 更进一步地，所述Step 4具体包括以下步骤：
[0052] Step 4‑1，进一步对式(11)采用带有皮亚诺型余项的泰勒公式展开，并忽略二阶以上参数项对展开结果所带来的误差影响，则得到下式：
[0053]
[0054] 式中：为干旱指数的多年平均值； 为短时间尺度下考虑陆地水储量变化影响作用的干旱指数； 为二阶以上参数项对展开结果所带来的误差影响； 为研究时间段内的干旱指数变化值；
[0055] Step 4‑2，根据闭合流域水量平衡公式，研究流域在特定时间段内的径流深可表述为：
[0056] Ri＝Pi‑ETi‑ΔSi (13)[0057] 式中：Ri为研究流域在对应研究时间段内的径流深，mm；
[0058] Step 4‑3，假设研究流域陆地水储量变化值的多年平均值为0，则该流域在特定时间段内的径流变化值表述为：
[0059]
[0060] 式中：ΔRi为研究流域在对应研究时间段内的径流深变化值，mm；Ri为研究流域的径流深序列值，mm；为研究流域在整个研究时间段内的平均径流深，mm；Pi和ETi分别为研究流域的降水量和蒸散发序列值，其单位均为mm；ΔSi为研究流域的陆地水储量变化值，mm；和 分别为研究流域的降水量多年平均值和蒸散发序列多年平均值，其单位均为mm；为陆地水储量变化值的多年平均值，mm，其数值假设为0；ΔPi和ΔETi分别为研究流域的降水量变化值和蒸散发变化值，其单位均为mm；
[0061] Step 4‑4，依据式(9)和式(11)分别可以得到下式：
[0062]
[0063]
[0064] Step 4‑5，将式(15)和式(16)代入式(14)则可以得到：
[0065]
[0066]
[0067] Step 4‑6，将式(12)代入式(17)可以得到：
[0068]
[0069] 式中：ΔP′i为有效降水量变化值，mm；
[0070] Step 4‑7，将式(18)中的 进行转化则可以得到下式：
[0071]
[0072] 式中：PETi为潜在蒸散发变化值，mm；
[0073] Step 4‑8，依据样本方差的定义，可以获得以下结果：
[0074]
[0075]
[0076]
[0077]
[0078]
[0079] 式中： 为径流变化的样本方差，mm‑2； 为蒸散发变化的样本方差，mm‑2； 为‑2 ‑2降水量变化的样本方差，mm ； 为潜在蒸散发变化的样本方差，mm ； 为陆地水储量‑2
变化的样本方差，mm ；N为样本统计个数；
[0080] Step 4‑9，对式(19)左右项分别求取平方和，然后各项分别乘以系数 并将式(20)至式(24)代入式(19)后可得：
[0081]
[0082] 式中：cov(Pi，ETi)为降水量和潜在蒸散发之间的协方差；cov(Pi，ΔSi)为降水量和陆地水储量之间的协方差；cov(PETi，ΔSi)为潜在蒸散发和陆地水储量之间的协方差；
[0083] Step 4‑10，通过将研究流域在短时间尺度下的水文变量序列依次代入式(25)，分别估算不同水文变量方差变化项所对应的比例系数，即可求取不同水文变量对径流变化的贡献作用。
[0084] 进一步地，所述Step 4‑7中， 的转化方法具体包括以下步骤：
[0085] Step 4‑7‑1，对于特定时间段内的干旱指数变化量 有如下定义：
[0086]
[0087] Step 4‑7‑2，将陆地水储量变化值的多年平均值假设为0，即有 则特定时间段内的有效降水量变化值ΔP′i为：
[0088]
[0089] Step 4‑7‑3，对式(26)进行转化，并将式(27)代入到式(26)则可以都得到：
[0090]
[0091]
[0092] 本发明的发明原理为：
[0093] 现阶段水文研究中所采用的原始Budyko公式仅适用于表征多年尺度下(一般为10年以上)水文变量之间的交互作用关系，且未能够提出行之有效的研究方法来对变化环境下的径流变化特征进行归因分析，也无法量化水文变量对径流变化的贡献作用大小。本发明通过有效考虑陆地水储量变化对流域水文循环过程的影响作用，成功构建改进Budyko公式，该公式不仅可以克服原始Budyko公式无法表征短时间尺度下流域水文变量交互作用关系的缺点，还可以有效考虑陆地水储量变化对该流域水文过程的影响作用，并结合样本方差定义提出一系列对应的径流变化分解方法，来量化变化环境下不同水文变量对径流变化的贡献作用。
[0094] 本发明具有以下有益效果：
[0095] (1)本发明能够充分满足变化环境下研究流域径流变化特征归因分析的有关要求，不仅能够有效量化不同水文变量对径流变化的贡献作用，还可以显著提高研究流域在变化环境下的径流预测准确度，为极端环境下的洪涝事件预测和水资源管理提供相应的理论支撑和研究依据；
[0096] (2)相比较于其他方法，Penman‑Monteith公式具有更高的准确度和可靠度，利用该公式可以准确描述研究流域在变化环境下的潜在蒸散发特征。此外，基于闭合流域水量平衡公式估算该流域的陆地水储量变化序列，能够克服原始Budyko公式未考虑研究流域陆地水储量变化对水文循环影响作用的不足，更加合理地描述变化环境下的水文变量演变规律；
[0097] (3)采用改进的Budyko公式可以同时考虑研究流域水量平衡和能量平衡条件，既能够表征水热耦合平衡状态下各主要水文变量之间的交互作用关系，还能够充分考虑陆地水储量变化对研究流域水文循环动态变化过程中的影响作用，确保该分析方法能够在短时间尺度下也得以应用。相比而言，原始Budyko公式仅仅能够应用于表征多年尺度(一般为10年以上)的水热耦合平衡状态下各主要水文变量之间的交互作用关系；
[0098] (4)本发明所采用的分解方法既能够就不同水文变量对径流变化的贡献作用进行充分量化分析，也能够估算不同水文变量之间的协同变化对研究流域径流变化的贡献作用，遵循严格的数理推导过程，有助于深入理解水文循环过程中不同水文变量对研究流域径流变化的影响机制。
附图说明
[0099] 图1为基于改进Budyko公式的径流变化定量分解方法具体流程；
[0100] 图2为某流域在2003年‑2014年间的月尺度陆地水储量变化值序列；
[0101] 图3为某流域主要水文变量在月尺度的修正Budyko公式拟合曲线；
[0102] 图4为某流域主要水文变量对径流变化的贡献作用占比。
具体实施方式
[0103] 下面结合以下附图，详细阐释实施案例。
[0104] 本发明提出了一种基于改进Budyko公式的径流变化定量分解方法，该方法不仅考虑了陆地水储量变化对水文循环过程的影响作用，还能够量化不同水文变量对径流变化的贡献作用。如图1所示为本发明方法流程图。
[0105] 一种基于改进Budyko公式的径流变化定量分解方法，包括以下步骤：
[0106] Step 1：获取研究流域的主要水文变量，包括研究流域的降水量、径流量、最高温度、最低温度、风速、相对湿度、日照时长、蒸散发等，并基于Penman‑Monteith公式计算该流域的潜在蒸散发；
[0107] Step 2：为有效考虑陆地水储量变化对流域水文循环的影响作用，通过闭合流域水量平衡公式估算该流域在研究时间段内的陆地水储量变化值；
[0108] Step 3：将Step 2中所获得的陆地水储量变化值耦合加入到原始Budyko公式中，通过构建改进的Budyko公式来表征研究流域水文变量在水热耦合平衡状态下的交互作用关系；最后，通过最小二乘法求取改进的Budyko公式中所包含的特征参数；改进的Budyko公式由于能够充分考虑陆地水储量变化对流域水文循环过程的影响作用，故而能够更加准确地描述短时间尺度下动态水文变化过程中各水文变量之间的交互作用关系；
[0109] Step 4：对Step 3中所获得的改进的Budyko公式采用皮亚诺型余项的泰勒公式展开，并依据样本方差的定义，对研究流域的径流变化部分进行定量分解。最后，根据分解结果求取不同水文变量对研究流域径流变化的贡献作用。
[0110] 所述Step 1具体为：
[0111] Step 1‑1，确定位于研究流域的气象站点，收集并整理计算月尺度潜在蒸散发所需要的水文变量，包括月尺度的最高温度、最低温度、风速、相对湿度、日照时长；
[0112] Step 1‑2，基于Penman‑Monteith公式估算研究流域对应的月尺度潜在蒸散发序列，其中，所述的Penman‑Monteith公式为：
[0113]
[0114] 式中：PET为研究流域的潜在蒸散发，mm·month‑1；Δ为饱和水汽压‑温度曲线斜‑1 ‑2 ‑1 ‑1率，kPa·℃ ；n为净辐射，MJ·m ·month ；γ为湿度计常数，kPa·℃ ；G为地表热通量，‑2 ‑1
MJ·m ·month ；T为距离地面2米处的空气平均温度，℃；u为距离地面2米处的平均风速，‑1
m·s ；es为距离地面2米处的饱和水汽压，kPa；ea为距离地面2米的实测水汽压，kPa；
[0115] Step 1‑3，将对应的水文变量代入式(1)进行计算，即可得到不同站点的月尺度潜在蒸散发；然后，基于不同气象站点的潜在蒸散发估算结果，通过反距离插值方法得到整个研究流域的平均潜在蒸散发。
[0116] 所述Step 1‑3具体包括以下步骤：
[0117] Step 1‑3‑1，根据研究流域的面积大小，将其空间划分成若干个网格，分别计算所有网格格点到各个气象站点的距离：
[0118]
[0119] 式中：(xi,yi)为研究流域内对应第i个格点的经纬度；(Xj,Yj)为对应第j个气象站点的经纬度；Dij即为第i个格点到第j个气象站点的距离；
[0120] Step 1‑3‑2，进一步求取各个气象站点对研究流域内第i个格点的影响权重系数，即为：
[0121]
[0122] 式中：wij为第j个气象站点对第i个格点的影响权重系数；m为气象站点的总个数；
[0123] Step 1‑3‑3，最后，求取不同水文变量在流域内各个格点的对应数值，进而求取该水文变量在研究流域内的平均值：
[0124]
[0125] 式中：fj为水文变量在第i个格点的对应数值；fj为第j个气象站点对应的气象数据观测值。
[0126] 进一步地，所述Step 2具体包括以下步骤：
[0127] Step 2‑1，将控制水文站的月尺度径流量转化为月尺度径流深，其计算方法为：
[0128]
[0129] 式中：R为研究流域气象站点的月径流深，单位为mm；R′为研究流域控制水文站的3 2
月径流量，m；A为研究流域的总面积，m；
[0130] Step 2‑2，基于各个气象站点的降水量观测数据，通过泰森多边形法求取研究流域的月尺度面平均降水量：
[0131]
[0132] 式中：P为整个研究流域的月尺度面平均降水量，mm；wj为第j个气象站点占整个研究流域的权重值；Pj为第j个气象站点的降水量观测值，mm；m为气象站点的总个数；
[0133] Step 2‑3，联合利用上述各式所求得的月尺度面平均降水量、月径流深以及蒸散发序列，基于闭合流域水量平衡公式，求取研究流域的月尺度陆地水储量变化值：
[0134] ΔS＝P‑R‑ET (7)
[0135] 式中：ΔS为研究流域在对应时间段内的陆地水储量变化值，mm；ET为研究流域实测蒸散发，mm。
[0136] 所述Step 3具体包括以下步骤：
[0137] Step 3‑1，基于原始Budyko公式求取多年尺度(一般规定为10年以上)下水文变量在水热耦合平衡状态下的交互作用关系，并基于最小二乘法求取该公式所包含的特征参数[0138]
[0139] 式中： 为该流域的蒸散发多年平均值，mm； 和 分别为该流域的潜在蒸散发多年平均值及降水量多年平均值，其单位均为mm；为干旱指数多年平均值，是一个无量纲参数；为该流域水热耦合平衡状态下对应的特征参数，是一个无量纲参数；
[0140] Step 3‑2，当考虑其他更短的时间尺度(如月尺度、季尺度或者年尺度等)时，则需通过下式求取研究流域在水热耦合平衡状态对应的水文变量交互作用关系：
[0141]
[0142] 式中：ETi为该流域在短时间尺度下的蒸散发值，mm；PETi和Pi分别为该流域短时间尺度下的潜在蒸散发值及降水量，其单位均为mm； 为短时间尺度下的干旱指数，是一个无量纲参数；n为该流域水热耦合平衡状态下对应的特征参数，同样为一个无量纲参数；
[0143] Step 3‑3，考虑陆地水储量变化对研究流域水文过程的影响作用，需按照下式求取研究流域的有效降水量：
[0144] P′i＝Pi‑ΔSi (10)[0145] 式中：P′i为研究流域有效降水量，mm；Pi为基于气象站得到的研究流域实测降水量，mm；ΔSi为在对应时间段内的研究流域陆地水储量变化值，mm；
[0146] Step 3‑4，用有效降水量P′i替代式(9)中的实测降水量Pi，则可获得考虑陆地水储量变化影响作用的改进的Budyko公式，并基于最小二乘法求取该公式所包含的特征参数n：
[0147]
[0148] 式中：P′i为研究流域的有效降水量，mm； 为短时间尺度下考虑陆地水储量变化影响作用的干旱指数，是一个无量纲参数。
[0149] 所述Step 4具体包括以下步骤：
[0150] Step 4‑1，进一步对式(11)采用带有皮亚诺型余项的泰勒公式展开，并忽略二阶以上参数项对展开结果所带来的误差影响，则得到下式：
[0151]
[0152] 式中：为干旱指数的多年平均值； 为短时间尺度下考虑陆地水储量变化影响作用的干旱指数； 为二阶以上参数项对展开结果所带来的误差影响； 为研究时间段内的干旱指数变化值；
[0153] Step 4‑2，根据闭合流域水量平衡公式，研究流域在特定时间段内的径流深可表述为：
[0154] Ri＝Pi‑ETi‑ΔSi (13)
[0155] 式中：Ri为研究流域在对应研究时间段内的径流深，mm；
[0156] Step 4‑3，假设研究流域陆地水储量变化值的多年平均值为0，则该流域在特定时间段内的径流变化值表述为：
[0157]
[0158] 式中：ΔRi为研究流域在对应研究时间段内的径流深变化值，mm；Ri为研究流域的径流深序列值，mm；为研究流域在整个研究时间段内的平均径流深，mm；Pi和ETi分别为研究流域的降水量和蒸散发序列值，其单位均为mm；ΔSi为研究流域的陆地水储量变化值，mm；和 分别为研究流域的降水量多年平均值和蒸散发序列多年平均值，其单位均为mm；为陆地水储量变化值的多年平均值，mm，其数值假设为0；ΔPi和ΔETi分别为研究流域的降水量变化值和蒸散发变化值，其单位均为mm；
[0159] Step 4‑4，依据式(9)和式(11)分别可以得到下式：
[0160]
[0161]
[0162] Step 4‑5，将式(15)和式(16)代入式(14)则可以得到：
[0163]
[0164] Step 4‑6，将式(12)代入式(17)可以得到：
[0165]
[0166] 式中：ΔP′i为有效降水量变化值，mm；
[0167] Step 4‑7，将式(18)中的 进行转化则可以得到下式：
[0168]
[0169] 式中：PETi为潜在蒸散发变化值，mm；
[0170] Step 4‑8，依据样本方差的定义，可以获得以下结果：
[0171]
[0172]
[0173]
[0174]
[0175]
[0176] 式中： 为径流变化的样本方差，mm‑2； 为蒸散发变化的样本方差，mm‑2； 为降‑2 ‑2水量变化的样本方差，mm ； 为潜在蒸散发变化的样本方差，mm ； 为陆地水储量变‑2
化的样本方差，mm ；N为样本统计个数；
[0177] Step 4‑9，对式(19)左右项分别求取平方和，然后各项分别乘以系数 并将式(20)至式(24)代入式(19)后可得：
[0178]
[0179] 式中：cov(Pi，ETi)为降水量和潜在蒸散发之间的协方差；cov(Pi，ΔSi)为降水量和陆地水储量之间的协方差；cov(PETi，ΔSi)为潜在蒸散发和陆地水储量之间的协方差；
[0180] Step 4‑10，通过将研究流域在短时间尺度下的水文变量序列依次代入式(25)，分别估算不同水文变量方差变化项所对应的比例系数，即可求取不同水文变量对径流变化的贡献作用。
[0181] 所述Step 4‑7中， 的转化方法具体包括以下步骤：
[0182] Step 4‑7‑1，对于特定时间段内的干旱指数变化量 有如下定义：
[0183]
[0184] Step 4‑7‑2，将陆地水储量变化值的多年平均值假设为0，即有 则特定时间段内的有效降水量变化值ΔP′i为：
[0185]
[0186] Step 4‑7‑3，对式(26)进行转化，并将式(27)代入到式(26)则可以都得到：
[0187]
[0188]
[0189] 将本发明应用于某流域的径流变化归因分析研究中，并采用2003年‑2014年间的月尺度径流量、降水量、最高温度、最低温度、相对湿度、风速、日照时长以及蒸散发作为模型的输入数据。
[0190] 基于Penman‑Monteith公式计算某流域各个气象站点的潜在蒸散发，然后通过反距离插值算法求取该流域各个格点的潜在蒸散发值，并进一步求得整个某流域的月尺度潜在蒸散发序列。
[0191] 基于闭合流域水量平衡公式，联合利用降水、径流深以及蒸散发数据求取某流域在研究时间段内的陆地水储量变化值。
[0192] 将所有的潜在蒸散发、降水量、陆地水储量变化值以及蒸散发分别代入改进的Budyko公式中，并基于最小二乘法求取该公式所包含的特征参数值，得到水热耦合平衡状态下水文变量之间的主要拟合关系。
[0193] 基于本发明所构建的径流变化定量分解方法，求取不同水文变量对应的系数作为该水文变量对某流域径流变化的贡献作用，进而量化不同水文变量对流域径流变化的贡献作用。
[0194] 图2显示了某流域在2003年‑2014年间的月尺度陆地水储量变化值序列；图3显示了某流域主要水文变量在月尺度的修正Budyko公式拟合曲线；图4显示了某流域主要水文变量对径流变化的贡献作用占比，其中，P_ΔS：降水量和陆地水储量协同变化值；PET_ΔS：潜在蒸散发和陆地水储量协同变化值；P_PET：降水量和潜在蒸散发协同变化值。由图2可以发现某流域在2003年‑2014年间的陆地水储量呈现较好的季节性变化特征；由图3可以发现某流域主要水文变量在月尺度上能够较好得符合修正Budyko公式拟合曲线；由图4则可以发现该流域的月径流变化主要是由于降水量变化以及降水量与陆地水储量之间的协同作用所共同导致的。
[0195] 本发明提出的方法不仅能够考虑陆地水储量变化对研究流域水文循环过程的影响作用，还能够准确量化不同水文变量对研究流域径流变化的贡献作用。
[0196] 此实施例仅为本发明较佳的具体实施方式，但本发明的保护范围并不局限于此，在阅读了本发明之后，本领域技术人员对本发明的各种等价形式的修改均落于本申请所附权利要求所限定的范围。
