CN114462518B - 考虑多要素空间依赖性的区域蒸散发变化归因分析方法 - Google Patents

Abstract

本发明公开了一种考虑多要素空间依赖性的区域蒸散发变化归因分析方法，所述方法包括：步骤一、搜集研究区的气象要素、地理要素、人类活动要素数据集及研究区内水文站点径流数据；步骤二、数据前处理；步骤三、利用空间聚类方法对研究区内不同流域进行聚类分区，在每个分区内采用多元线性回归方法筛选Budyko水热耦合平衡方程中流域特征综合参数的主控因子；步骤四、基于层次贝叶斯模型构建考虑多要素空间依赖性的流域水热耦合平衡模型；步骤五、利用弹性系数法定量区分各要素对区域蒸散发变化的影响。本发明方法可以考虑多要素空间依赖性对于区域蒸散发变化的影响，提高了流域水热耦合平衡模型的模拟精度，进而有助于准确定量区分各要素对蒸散发变化的贡献，对于区域水资源调控和科学管理、农业生产和生态建设具有重要意义。

Description

考虑多要素空间依赖性的区域蒸散发变化归因分析方法

技术领域

本发明涉及流域蒸散发演变及归因领域，具体涉及一种考虑多要素空间依赖性的区域蒸散发变化归因分析方法。

背景技术

地表蒸散发是陆地水热平衡的重要组成部分，其在地表与大气的物质循环与能量交换中扮演着重要角色。地表蒸散发将约60％的降水量输送回大气，进而对降水过程产生影响；同时，蒸散发过程中吸收的潜热通量又保证了地表气温的稳定。因此区域蒸散发变化对区域水文循环及气候系统变化具有重要意义。由于蒸散发量难以直接测定，需要寻找有效途径或方法对其进行准确估算；同时，区域蒸散发过程受当地气候、下垫面及人类活动等多种因素的影响，如何对其进行归因及量化是当下的研究重点之一。

在蒸散发研究中，重点是对实际蒸散发的获取。实际蒸散发的计算过程以气候学计算方法为主，但是该方法计算繁琐，并且需要大量基础数据，实际使用较为繁琐。此外，一些理论被广泛使用，如Penman正比假设理论、Bouchet互补相关理论、Budyko水热耦合平衡理论等。Penman正比假设认为实际蒸散量与可能蒸散量成正比，但是这一假设过于粗糙，不能充分反映蒸散发过程中各环节的影响。Bouchet互补相关理论认为区域实际蒸散量与潜在蒸散量之间为互补关系，该理论与正比假设相反，这说明在特定地区，实际蒸散量与潜在蒸散量只能符合某一种关系，因此他们的使用具有一定局限性。

Budyko水热耦合平衡理论认为区域的水分供给(降水)与大气蒸散发需求(潜在蒸散发)控制区域长期水量平衡，这一假设阐明了流域蒸散发和降水及能量间的关系。该理论结构简单，物理机制明确，因此被广泛用于区域蒸散发或产流估算。在Budyko方程的实际使用过程中，学者们发现除气候特征外，流域植被、土壤、地形等特征也会对水文循环过程产生影响，因此流域特征综合参数被纳入到Budyko框架中。同时，包含流域特征综合参数的Budyko方程的微分形式极易推导出径流或蒸散发对气候变量及流域特征变量的敏感性解析式，因而可用于定量归因区域径流或蒸散发变化。

基于Budyko水热耦合平衡方程的区域蒸散发研究多基于流域独立性假设，而忽视了不同区域多要素空间依赖性对蒸散发的影响，即某一流域的蒸散发变化不仅受区域内气候、自然地理和社会经济等多要素相互作用关系的影响，也同时受到邻近流域多要素的协同作用影响，对空间依赖性的考虑不足将影响蒸散发变化归因分析结果的准确性，导致区域蒸散发变化响应机制区域差异化解析不充分。因此，亟需构建考虑空间依赖性的流域水热平衡模型，进而准确反映蒸散发变化的区域差异，提高区域蒸散发变化归因分析的准确性，这对水资源区域化管理的科学决策具有重要意义。

发明内容

本发明的目的在于提供一种考虑多要素空间依赖性的区域蒸散发变化归因分析方法，旨在解决定量确定多要素及其空间依赖性对蒸散发变化影响的问题。

为实现上述目的，本发明提供如下技术方案：

一种考虑多要素空间依赖性的区域蒸散发变化归因分析方法，包括以下步骤：

步骤一、数据收集：收集目标研究区在研究时段内的气象要素数据、地理要素数据、人类活动要素数据以及径流数据；利用气象数据及FAO Penman–Monteith方程计算目标研究区在研究时段内的潜在蒸散发；

步骤二、数据前处理：包括缺失数据的插补延长、数据空间插值和重采样、数据转换中的一项或几项；计算各要素的面平均值；

步骤三、多要素空间聚类及流域特征综合参数主控因子识别：根据降水数据、径流数据及潜在蒸散发数据求解各流域的Budyko水热耦合平衡方程中的流域特征综合参数；根据已有研究及研究区域实际情况，挑选可以反映流域特征的因子；基于多要素空间聚类方法，利用挑选的反映流域特征的因子对研究区内的多个子流域进行空间聚类分区，并对聚类结果进行分析；对于聚类后的各分区，基于多元线性回归方法，识别影响流域特征综合参数的主控因子；

步骤四、构建考虑多要素空间依赖性的流域水热平衡模型：根据先验知识确定流域特征综合参数的先验分布，先验分布中的未知参数为超参数；根据流域相邻关系确定因子间的空间依赖关系；将因子空间依赖关系作为变量纳入到超参数和主控因子的逐步回归模型中，该模型即为流域特征综合参数和主控因子间的层次贝叶斯模型；将基于层次贝叶斯模型估计的流域特征综合参数经验关系式带入Budyko水热平衡关系式中，此时的Budyko方程即为考虑多要素空间依赖性的流域水热耦合平衡模型；分区利用层次贝叶斯模型估计Budyko方程中的流域特征综合参数是本发明的一个重要改进点；

步骤五、定量归因区域蒸散发变化：基于弹性系数法，定量分析目标研究区域内各分区各影响因素对蒸散发变化的影响。

进一步的，步骤一中所述气象要素数据包括降水、气温、气压、蒸发、相对湿度、风速风向、日照及辐射数据；地理要素数据包括高程、坡度、坡向、植被覆盖、土壤饱和水力传导度数据；人类活动要素数据包括单位面积农业灌溉用水、单位面积建成区比例、单位面积GDP数据。

步骤一、FAO Penman–Monteith方程的具体计算步骤为：

式中，PE为蒸散发量，mm/d；Δ为饱和水汽压斜率，kPa/℃；R_n为地表净辐射，MJ/(m²*d)；G为土壤热通量，MJ/(m²*d)；γ为干湿表常数，kPa/℃；T_mean为日平均温度，℃；u₂为2米高处风速，m/s；e_s为饱和水汽压，kPa；e_a为实际水汽压，kPa。

进一步的，步骤二所述的各要素面平均值的计算方法为泰森多边形法，其具体步骤为：利用Gis软件将研究区内的各待研究对象值提取到气象站点；将研究区域内所有相邻气象站连成三角形；依次做这些三角形各边的垂直平分线；将每个气象站点周围的若干垂直平分线围城一个多边形，该图形称作泰森多边形；用泰森多边形内包含的唯一气象站点处的要素值表示该多边形区域内的平均值；区域内所有的泰森多边形构成多边形网，将每个多边形在多边形网的权重系数与要素区域平均值相乘后累加可获得要素在整个研究区域的面平均值：

X＝f₁X₁+f₂X₂+…f_mX_m

其中，m为站点个数，X₁,X₂,…,X_m为各站点处的要素值，f₁,f₂,…,f_m为各多边形的权重系数。

进一步的，步骤三所述的流域水量平衡方程是指，在多年平均尺度下，流域蓄水变化(ΔS)可以忽略不计，流域径流可视作降水与蒸散发的差值，即：

E＝P-R

其中，E为多年平均实际蒸散量，P为多年平均降水量；R为多年平均径流量。

步骤三所述Budyko水热耦合平衡方程有多种形式，本发明采用的是Choudhury-Yang Budyko方程，其具体形式如下：

其中，ET₀为潜在蒸散量，n为Budyko参数，与流域属性有关。

进一步的，步骤三中反映流域特征的因子包括：气候特征因子、下垫面特征因子、地形特征因子、人类活动特征因子；所述气候特征因子包括年平均降水量、逐日平均气温、气候季节性；所述下垫面特征因子包括土壤饱和水力传导度、植被覆盖；所述地形特征因子包括高程、坡度、坡向、地形湿度指数；所述人类活动特征因子包括单位面积灌溉用水量、单位面积建成区比例、单位面积GDP。

进一步的，步骤三所述的多要素空间聚类方法为自组织适应地图(SOM)神经网络聚类方法。

进一步的，步骤四所述的层次贝叶斯模型具体计算步骤为：

根据先验信息得出流域特征综合参数n的先验分布，以正态分布为例，构建流域特征综合参数n的层次贝叶斯模型：

n_i～N(μ_i,σ²)，μ_i＝α₀+α₁x₁+…α_kx_k+b_i

其中，μ_i为随机变量n的均值，σ²为随机变量n的方差，且σ²～gamma(0.001,0.001)，μ_i为该分布的超参数，即μ_i为待估参数；x₁,x₂,…,x_k为对n有影响的主控因子，α₀,α₁,…,α_k为系数，并且α₀,α₁,…,α_k服从正态分布；b_i是对各自变量空间依赖性的描述，其先验分布由CAR模型给出，即：

b_i～NormalCAR(adj_i,w_i,num_i,τ)

其中adj_i为第i个区域的邻近地区集合；w_i为空间权重因子，邻近地区为1，不相邻地区为0；num_i为第i个区域的邻近区域个数；τ～gamma(0.5,0.0005)，表示一种精度。

进一步的，步骤五所述的弹性系数法具体计算过程如下：

其中，E、P、ET₀和n的意义同前文所述；S_{E_P}、

和S_{E_n}分别表示蒸散发对降水、蒸散发和Budyko流域特征综合参数的弹性；v_i表示影响流域特征综合参数的因子。

本发明的优点和有益效果是：利用层次贝叶斯模型考虑了因子空间依赖性的作用，提高了流域特征综合参数的模拟精度；将利用层次贝叶斯模型估计的流域特征综合参数经验公式纳入到Budyko框架中，有助于准确归因和定量分解蒸散发变化。

附图说明

下面结合附图和实施例对本发明作进一步说明。

图1为本发明实施例的流程示意图；

图2为研究区域Budyko流域特征综合参数n的分布图；

图3为研究区域要素演变特征图；

图4为逐步回归法与层次贝叶斯模型对比；

图5为研究区域蒸散发演变归因。

具体实施方式

一种考虑多要素空间依赖性的区域蒸散发变化归因分析方法，包括以下步骤：

步骤一、数据收集：收集目标研究区在研究时段内的气象要素数据、地理要素数据、人类活动要素数据以及径流数据；利用气象数据及FAO Penman–Monteith方程计算目标研究区在研究时段内的潜在蒸散发。

气象要素数据包括降水、气温、气压、蒸发、相对湿度、风速风向、日照及辐射数据；地理要素数据包括高程、坡度、坡向、植被覆盖、土壤饱和水力传导度数据；人类活动要素数据包括单位面积农业灌溉用水、单位面积建成区比例、单位面积GDP数据。

FAO Penman–Monteith方程的具体计算步骤为：

假设作物植株高度为0.12m，固定的作物表面阻力为70m/s，反射率为0.23，即认为非常类似于表面开阔、高度一致、生长旺盛、完全遮盖地面而水分充分适宜的绿色草地的蒸散量为潜在蒸散量：

式中，PE为蒸散发量(mm/d)；Δ为饱和水汽压斜率(kPa/℃)；R_n为地表净辐射(MJ/(m²*d))；G为土壤热通量(MJ/(m²*d))；γ为干湿表常数(kPa/℃)；T_mean为日平均温度(℃)；u₂为2米高处风速(m/s)；e_s为饱和水汽压(kPa)；e_a为实际水汽压(kPa)。

式中关键参数的计算如下：

步骤1-1，计算Δ：

其中，T选用日平均气温；

步骤1-2，计算R_n：

R_n＝R_ns-R_nl (3)

即净辐射R_n为收入的短波辐射R_ns和支出的净长波辐射R_nl之差；

R_ns＝(1-α)R_s (4)

其中，α＝0.23，此处取绿色草地参考作物的反照率；

R_s为太阳辐射；a_s和b_s表示回归系数，推荐使用a_s＝0.25，b_s＝0.50；n为实际日照时数(h)，N为最大可能日照时数(h)；

w_s为日出时角；

表示纬度(rad)；

δ为太阳磁偏角；J为日旭，取值范围为1到365或366，1月1日为1；

G_sc为太阳常数(MJ/(m²*min)),G_sc＝0.0820，

和δ意义同上，R_a为日地球外辐射；

d_r为日地平均距离，J的含义同上；

其中σ为斯蒂芬-玻尔兹曼常数(MK·K^-4·m^-2·d^-1)，σ＝4.903×10^-9；

T_max,K＝T_max+272.15 (12)

T_min,K＝T_min+272.15 (13)

T_max,K、T_min,K为最高、最低绝对温度(K)，T_max、T_min为日最高、最低温度(℃)；

e_a＝RH×e_s (14)

e_a为实际水汽压，RH为相对湿度，e_s为饱和水汽压；

T_max、T_min意义同上；

R_so＝(a_s+b_s)R_a (18)

R_a、a_s、b_s意义同上；

步骤1-3，计算G：

G≈0 (19)

步骤1-4，计算γ：

γ＝0.665×10³P (20)

其中，P为大气压；

步骤1-5，计算T_mean：

其中，T_max、T_min意义同上；

步骤1-6，计算u₂：

其中，u_z为10m高度处的风速，z＝10。

步骤二、数据前处理：数据前处理包括缺失数据的插补延长、数据空间插值和重采样、数据转换；缺失数据利用线性内插或是线性外延方法进行插补延长；数据的空间插值采用反距离加权平均方法，并利用GIS软件中的Create Fishnet功能按1km分辨率在空间上生成均匀分布的采样点，以此获取采样点位置处的数据序列；数据转换是指根据研究需要对数据进行如对数化、无量纲化或比例化等处理；数据面平均值计算采用泰森多边形法；

数据转换的方法包括对数化、无量纲化及比例化，其计算公式分别为：

对数化：

以10为底，对数据取对数：

X_i＝logx_i (23)

以e为底，对数据取对数：

X_i＝lnx_i (24)

其中x为变量变换前的数值，X为变换后的数值；i表示第i个变量，i＝1,2,…,k，k为因子个数；

无量纲化：

其中i、x、X的意义同前文所述；j表示变量的第j个值，j＝1,2,…,l，l为各自变量的序列长度；

比例化：

X_i＝αx_i (26)

其中i、x、X的意义同前文所述,α为比例系数。

各要素面平均值的计算方法为泰森多边形法，其具体步骤为：

步骤2-1：利用Gis软件将研究区内的各待研究对象值提取到气象站点；

步骤2-2：将研究区域内所有相邻气象站连成三角形；

步骤2-3：依次做这些三角形各边的垂直平分线；

步骤2-4：将每个气象站点周围的若干垂直平分线围城一个多边形，该图形称作泰森多边形；

步骤2-5：用泰森多边形内包含的唯一气象站点处的要素值表示该多边形区域内的平均值；

步骤2-6：区域内所有的泰森多边形构成多边形网，将每个多边形在多边形网的权重系数与要素区域值相乘后累加可获得要素在整个研究区域的面平均值：

X＝f₁X₁+f₂X₂+…f_mX_m (27)

其中，m为站点个数，X₁,X₂,…,X_m为各站点处的要素值，f₁,f₂,…,f_m为各多边形的权重系数。

步骤三、多要素空间聚类及流域特征综合参数主控因子识别：

步骤3-1：根据降水数据、径流数据及潜在蒸散发数据求解各流域的Budyko水热耦合平衡方程中的流域特征综合参数；

所述的求解Budyko水热耦合平衡方程中的流域特征综合参数是指利用降水、径流(或蒸散发)及潜在蒸散发数据，根据Choudhury-Yang Budyko水热平衡方程进行求解，Budyko方程的具体形式如下：

在多年平均尺度下，流域蓄水变化(ΔS)可以忽略不计，根据水量平衡公式R＝P-E，E可用P-R替代。其中，P为降水量，E为实际蒸散发量，ET₀为潜在蒸散发量，R为径流量，n为Budyko参数，与流域气候、土壤、植被覆盖等有关。

步骤3-2：根据已有研究及研究区域实际情况，挑选可以反映流域特征的因子；

其中反映流域特征的因子包括：气候特征因子、下垫面特征因子、地形特征因子、人类活动特征因子；所述气候特征因子包括年平均降水量、逐日平均气温、气候季节性；所述下垫面特征因子包括土壤饱和水力传导度、植被覆盖；所述地形特征因子包括高程、坡度、坡向、地形湿度指数；所述人类活动特征因子包括单位面积灌溉用水量、单位面积建成区比例、单位面积GDP。

步骤3-3：基于多要素空间聚类方法，利用挑选的反映流域特征的因子对研究区内的多个子流域进行空间聚类分区，并对聚类结果进行分析；

本实施例的多要素空间聚类方法为自组织适应地图(SOM)神经网络，SOM聚类方法可以使用MATLAB的Neural Net Clustering工具包完成。其具体步骤如下：

3-3-1：以向量X＝(x₁,x₂,…,x_k)作为神经网络的输入向量；

3-3-2：确定神经网络结构，其中输入向量的元素数量为k；

3-3-3：对连接权重向量W_i＝(w_i1,w_i2,…,w_ik)随机赋值，其取值范围为[0,1]，其中1≤i≤j，表示在输出层中共有j个权重向量及j个输出向量；

3-3-4：对输入向量做归一化处理并计算欧式距离d_i：

3-3-5：对连接权重向量w_ik进行迭代：

其中，t表示时间，η为学习效率；

3-3-6：重复以上步骤，直至学习效率小于规定的最小值，即η≤η_min，此时权重向量位于向量聚类的中心，因此性质相似的向量的聚集分布在输出层中，以此实现了多要素空间聚类的目的。

3-4对于聚类后的各分区，基于多元线性回归方法，识别影响流域特征综合参数的主控因子；所述的多元线性回归方法具体步骤为：

3-4-1：对选取的k个主控因子X₁,X₂,…,X_k作为回归自变量，分别同流域特征综合参数n建立一元回归模型：

n＝β₀+β_iX_i+ε,i＝1,2,…,k (32)

计算各回归系数的F检验统计量值，记为F₁ ⁽¹⁾,F₂ ⁽¹⁾,…,

记最大值为

取显著性水平α＝0.05，则当F_i ⁽¹⁾＞F_0.05时，将自变量X_i引入回归模型；

3-4-2：建立因变量n与自变量子集{X_i,X₁},…,{X_i,X_j},…,{X_i,X_k}，的二元回归模型，其中j＝1,…,i-1,i+1,…,k，并分别计算回归系数的F统计量值，记最大值为

若

则引入二元变量{X_i,X_j}；

3-4-3：重复上述步骤，分别建立因变量n与自变量子集的多元回归模型，直至无变量引入为止。

步骤四、构建考虑多要素空间依赖性的流域水热平衡模型：基于层次贝叶斯模型，进一步量化各分区主控因子及其空间依赖性对流域特征综合参数的影响，构建考虑多要素空间依赖性的流域水热耦合平衡模型；根据先验知识确定流域特征综合参数的先验分布，先验分布中的未知参数为超参数；根据流域相邻关系确定因子间的空间依赖关系；将因子空间依赖关系作为变量纳入到超参数和主控因子的逐步回归模型中，该模型即为流域特征综合参数和主控因子间的层次贝叶斯模型；将基于层次贝叶斯模型估计的流域特征综合参数经验关系式带入Budyko水热平衡关系式中，此时的Budyko方程即为考虑多要素空间依赖性的流域水热耦合平衡模型；

步骤四所述的层次贝叶斯模型具体计算步骤为：

4-1：检验流域特征综合参数n的先验分布，以正态分布为例：

n_i～N(μ_i,σ²) (33)

其中，μ_i为随机变量n的均值，σ²为随机变量n的方差，且σ²～gamma(0.001,0.001)，μ_i为该分布的超参数，即μ_i为待估参数；

4-2：构建层次贝叶斯模型：

μ_i＝α₀+α₁x₁+…α_kx_k+b_i (34)

其中，x₁,x₂,…,x_k为经逐步多元回归识别出的影响流域特征综合参数n的主控因子，α₀,α₁,…,α_k为系数，并且α₀,α₁,…,α_k服从正态分布；b_i是对各自变量空间依赖性的描述，其先验分布由CAR模型给出，即：

b_i～NormalCAR(adj_i,w_i,num_i,τ) (35)

其中adj_i为第i个区域的邻近地区集合；w_i为空间权重因子，邻近地区为1，不相邻地区为0；num_i为第i个区域的邻近区域个数；τ～gamma(0.5,0.0005)，表示一种精度。

步骤五、定量归因区域蒸散发变化：基于弹性系数法，定量分析目标研究区域内各分区各影响因素对蒸散发变化的影响。

弹性系数法具体计算过程如下：

其中，E表示蒸散发，P表示降水，ET₀表示潜在蒸散发，n表示流域特征综合参数；S_{E_P}、

和S_{E_n}分别表示蒸散发对降水、蒸散发和Budyko流域特征综合参数的弹性；v_i表示影响流域特征综合参数的因子。

实施例1：

以长江流域下游某一区域连续30年的数据为案例进行分析，对本发明技术方法的具体应用作进一步说明，具体应包括以下步骤：

步骤一、数据收集：

挑选研究区域内共计15个水文站点控制的子流域为研究流域。研究区内的气象站点数据从中国气象局国家气象信息中心(NMIC)(http://cdc.cma.gov.cn)下载。30m分辨率的DEM数据从USGS/NASA航天飞机雷达地形测绘任务(SRTM)数字高程数据库(http://srtm.csi)下载。土壤饱和水力传导度Kc数据从国家冰川冻土沙漠科学数据中心(https://www.crensed.ac.cn/portal/)获取，分辨率为1km。1982-2013年的GIMMS NDVI数据从NASA官网(https://www.nasa.gov/nex)获取，2014-2018年的SPOT/VEGETATION NDVI数据从中科院中国科学院资源与环境数据云平台(http://www.resdc.cn)获取，空间分辨率为1km。建成区面积和灌溉用水量由各省及各市统计年鉴获取。利用气象数据根据FAO Penman–Monteith方程计算潜在蒸散发量，主要是利用降水、气温、辐射、气压、相对湿度、风速及日照时数等数据。水文站点的径流数据从水文年鉴中获取。

表1长江流域蒸散发变化影响因素数据准备

序号	数据类型	时间分辨率	空间分辨率	数据类型
					1	气象数据	逐日	-	站点
2	DEM	-	30m	栅格
					3	土壤Kc数据	-	1km	栅格
4	NDVI	逐年	1km	栅格
					5	建成区面积	逐年	-	矢量
6	灌溉用水	逐年	-	矢量
					7	径流数据	逐年	-	站点

步骤二、数据前处理：

采用反距离加权平均法对栅格数据进行空间插值，使各要素的空间分辨率一致，并计算栅格数据的面平均值；对于站点数据，利用泰森多边形计算面平均值；

步骤三、多要素空间聚类及流域特征综合参数主控因子识别：

利用径流、降水及潜在蒸散发数据构建各水文站点控制流域的Budyko水热耦合平衡方程，并计算Choudhury-Yang Budyko模型的流域特征综合参数n(图2)；根据已有研究及研究区域实际情况挑选可以充分反映流域特征的因子，本实施例挑选的要素为降水(P/mm)、气温(T/℃)、降水季节性(SI)、高程(ELEV/m)、坡度(SLP/°)、地形湿度指数(TWI)、土壤饱和水力传导度(KC/cm·d-1)、植被覆盖度(M)、单位面积灌溉水量(IRR/mm)及单位面积建成区面积比例(B/％)，其中降水季节性由降水数据计算而来，地形湿度指数由地形数据计算而来，具体属性见表2。

表2流域属性

部分具有时变特征的要素的演变特征图3所示。降水在流域上游部分区域及下游区域呈现上升趋势，在流域中游呈现下降趋势；气温在流域上游和下游呈现下降区域，在流域中游呈现上升趋势；潜在蒸散发在流域上游为上升趋势，在流域下游为下降趋势；植被覆盖仅在流域西南和北部部分地区有所下降，在其余区域均呈现上升趋势；降水季节性的变化趋势在空间上分布不均匀，上游主要为下降趋势，下游主要呈现上升趋势；单位面积灌溉水量的变化趋势在空间上分布不均匀，流域上游主要为下降趋势，下游地区主要为上升趋势。

根据所选的要素，依据SOM方法对研究区域进行空间聚类，被分到一类的流域具有相同的属性，15个小流域共被分为两类。对流域特征综合参数n分区进行逐步多元回归，可知影响区域1流域特征综合参数n的主要影响因子为气温、土壤饱和水力传导度、高程、地形湿度和单位面积灌溉水量；影响区域2流域特征综合参数n的主要影响因子为降水、土壤饱和水力传导度、坡度、地形湿度指数和单位面积灌溉用水量。

步骤四、构建考虑多要素空间依赖性的流域水热平衡模型：

首先对流域特征综合参数的先验分布进行检验，经检验，流域特征综合参数符合正态分布，其中流域特征综合参数先验分布中的期望均值为超参数，即待估参数；对所选的子流域进行编号，并构建空间相邻矩阵，相邻区域记为1，否则记为0；利用CAR模型给出空间依赖性的描述；将描述空间依赖性的参数作为变量，纳入到待估参数和主控因子的线性回归模型中，该模型即为流域特征综合参数和主控因子的层次贝叶斯模型。分别用构建好的逐步多元回归模型和层次贝叶斯模型对流域特征综合参数进行模拟，并将其与用Budyko方程直接计算的结果对比。选用决定系数(R²)、相对误差(RB)、均方根误差(RMSE)作为评价模型模拟效果的指标，其中决定系数越接近1，相对误差和均方根误差越小表明模型效果越好，其结果如图4所示。层次贝叶斯法的决定系数为0.9814，接近1，而逐步回归法的决定系数仅为0.7787；层次贝叶斯法的相对误差为1.68％，均方根误差为0.03，小于逐步回归法的6.67％和0.30。结果表明考虑空间依赖性的层次贝叶斯模型可以改进对流域特征综合参数n的模拟效果。

步骤五、蒸散发演变归因分析：

利用弹性系数法，识别并量化研究区域30年来各要素对蒸散发变化的影响(图5)。如图所示，对于区域1而言，造成蒸散发变化的主要原因为降水、气温、潜在蒸散发和灌溉用水变化；其中气温变化是最主要原因，气温升高导致区域1潜在蒸散发在研究时段内增加了约12.7％。造成区域2蒸散发变化的主要原因为降水、潜在蒸散发和灌溉用水变化；其中降水是主要因素，降水增多导致区域2蒸散发在研究时段内增加了约7.2％。

最后应说明的是，以上仅用以说明本发明的技术方案而非限制，尽管参照较佳布置方案对本发明进行了详细说明，本领域的普通技术人员应当理解，可以对本发明的技术方案进行修改或者等同替换，而不脱离本发明技术方案的精神和范围。

Claims (9)

1.一种考虑多要素空间依赖性的区域蒸散发变化归因分析方法，其特征在于，所述方法包括以下步骤：

步骤一、数据收集：收集目标研究区在研究时段内的气象要素数据、地理要素数据、人类活动要素数据以及径流数据；利用气象数据及FAO Penman–Monteith方程计算目标研究区在研究时段内的潜在蒸散发；

步骤二、数据前处理：包括缺失数据的插补延长、数据空间插值和重采样、数据转换中的一项或几项；计算各要素的面平均值；

步骤三、多要素空间聚类及流域特征综合参数主控因子识别：

3-1根据降水数据、径流数据及潜在蒸散发数据求解各流域的Budyko水热耦合平衡方程中的流域特征综合参数；

3-2根据已有研究及研究区域实际情况，挑选可以反映流域特征的因子；

3-3基于多要素空间聚类方法，利用挑选的反映流域特征的因子对研究区内的多个子流域进行空间聚类分区，并对聚类结果进行分析；

3-4对于聚类后的各分区，基于多元线性回归方法，识别影响流域特征综合参数的主控因子；

步骤四、构建考虑多要素空间依赖性的流域水热平衡模型：根据先验知识确定流域特征综合参数的先验分布，先验分布中的未知参数为超参数；根据流域相邻关系确定因子间的空间依赖关系；将因子间的空间依赖关系作为变量纳入到超参数和主控因子的逐步回归模型中，该模型即为流域特征综合参数和主控因子间的层次贝叶斯模型；将基于层次贝叶斯模型估计的流域特征综合参数经验关系式带入Budyko水热平衡关系式中，此时的Budyko方程即为考虑多要素空间依赖性的流域水热耦合平衡模型；

步骤五、定量归因区域蒸散发变化：基于弹性系数法，定量分析目标研究区域内各分区各影响因素对蒸散发变化的影响。

2.根据权利要求1所述的一种考虑多要素空间依赖性的区域蒸散发变化归因分析方法，其特征在于，步骤一中所述气象要素数据包括降水、气温、气压、蒸发、相对湿度、风速风向、日照及辐射数据；地理要素数据包括高程、坡度、坡向、植被覆盖、土壤饱和水力传导度数据；人类活动要素数据包括单位面积农业灌溉用水、单位面积建成区比例、单位面积GDP数据。

3.根据权利要求1所述的一种考虑多要素空间依赖性的区域蒸散发变化归因分析方法，其特征在于，步骤一FAO Penman–Monteith方程的具体计算步骤为：

式中，PE为蒸散发量，mm/d；△为饱和水汽压斜率，kPa/℃；R_n为地表净辐射，MJ/(m²*d)；G为土壤热通量，MJ/(m²*d)；γ为干湿表常数，kPa/℃；T_mean为日平均温度，℃；u₂为2米高处风速，m/s；e_s为饱和水汽压，kPa；e_a为实际水汽压，kPa。

4.根据权利要求1所述的一种考虑多要素空间依赖性的区域蒸散发变化归因分析方法，其特征在于，步骤二所述的各要素的面平均值的计算方法为泰森多边形法，其具体步骤为：

步骤2-1：利用Gis软件将研究区内的各待研究对象值提取到气象站点；

步骤2-2：将研究区域内所有相邻气象站连成三角形；

步骤2-3：依次做这些三角形各边的垂直平分线；

步骤2-4：将每个气象站点周围的若干垂直平分线围城一个多边形，该多边形称作泰森多边形；

步骤2-5：用泰森多边形内包含的唯一气象站点处的要素值表示该多边形区域内的平均值；

步骤2-6：区域内所有的泰森多边形构成多边形网，将每个多边形在多边形网的权重系数与要素区域平均值相乘后累加可获得要素在整个研究区域的面平均值：

X＝f₁X₁+f₂X₂+…f_mX_m

其中，m为站点个数，X₁,X₂,…,X_m为各站点处的要素值，f₁,f₂,…,f_m为各多边形的权重系数。

5.根据权利要求1所述的一种考虑多要素空间依赖性的区域蒸散发变化归因分析方法，其特征在于，步骤三3-1所述的Budyko方程的具体形式如下：

在多年平均尺度下，根据水量平衡公式R＝P-E，E用P-R替代，其中，P为降水量，E为实际蒸散发量，ET₀为潜在蒸散发量，R为径流量，n为Budyko流域特征综合参数。

6.根据权利要求1所述的一种考虑多要素空间依赖性的区域蒸散发变化归因分析方法，其特征在于，步骤三3-2中反映流域特征的因子包括：气候特征因子、下垫面特征因子、地形特征因子、人类活动特征因子；所述气候特征因子包括降水量、气温、气候季节性；所述下垫面特征因子包括土壤饱和水力传导度、植被覆盖；所述地形特征因子包括高程、坡度、坡向、地形湿度指数；所述人类活动特征因子包括单位面积灌溉用水量、单位面积建成区比例、单位面积GDP。

7.根据权利要求1所述的一种考虑多要素空间依赖性的区域蒸散发变化归因分析方法，其特征在于，步骤三3-3所述的多要素空间聚类方法为自组织适应地图神经网络聚类方法。

8.根据权利要求1所述的一种考虑多要素空间依赖性的区域蒸散发变化归因分析方法，其特征在于，步骤四所述的层次贝叶斯模型具体计算步骤如下：

根据已有的流域特征综合参数n的样本信息和先验信息，得出n的先验分布；以n的先验分布为正态分布，构建n的层次贝叶斯模型：

n_i～N(μ_i,σ²)，μ_i＝α₀+α₁x₁+…α_kx_k+b_i

其中，μ_i为随机变量n的均值，σ²为随机变量n的方差，且σ²～gamma(0.001,0.001)，μ_i为该分布的超参数，即μ_i为待估参数；x₁,x₂,…,x_k为可能影响n的主控因子，α₀,α₁,…,α_k为系数，并且α₀,α₁,…,α_k服从正态分布；b_i是对各自变量空间依赖性的描述，其先验分布由CAR模型给出，即：

b_i～NormalCAR(adj_i,w_i,num_i,τ)

其中adj_i为第i个区域的邻近地区集合；w_i为空间权重因子，邻近地区为1，不相邻地区为0；num_i为第i个区域的邻近区域个数；τ～gamma(0.5,0.0005)，表示一种精度。

9.根据权利要求1所述的一种考虑多要素空间依赖性的区域蒸散发变化归因分析方法，其特征在于，步骤五所述的弹性系数法的具体计算过程如下：

其中，S_{E_P}、

和S_{E_n}分别表示蒸散发对降水、潜在蒸散发和Budyko流域特征综合参数的弹性；v_i表示影响流域特征综合参数的因子。

Images