一种确定蒸散发变化主因及判别因素间耦合作用的方法

摘要

本发明公开了一种确定蒸散发变化主因及判别因素间耦合作用的方法，搜集影响实际蒸散发变化的多个影响因子（包括流域下垫面特征、气候季节性特征以及社会经济指标），基于Budyko水热耦合平衡方程计算流域实际蒸散发量，利用旋转经验正交函数（REOF）对流域进行蒸散发气候敏感性分区，在分区的基础上对每个区域影响蒸散发变化的多种因素进行统计分析，利用多元自适应回归样条（MARS）方法建立模型，从而从多个因素中搜寻出最主要的影响因素，为科学研究以及防洪抗旱等生产工作提供科学依据。

说明书

技术领域

本发明涉及一种确定蒸散发变化主因及判别因素间耦合作用的方法，属于水文水资源应用技术领域。

背景技术

在气候变化过程中，实际蒸散发的变化是一个十分重要的衡量指标，实际蒸散发不仅在陆地-大气系统的能量收支中占有重要的地位，在水量平衡的计算中也十分重要。作为联系能量平衡系统和水量平衡系统的一种水文变量，分析实际蒸散发的变化对研究气候变化具有十分重要的作用。

在一个较大的流域内，影响实际蒸散发的变化因素有很多且具有一定的区域差异性。综合前人的研究可将这些因素大致分为4类：①区域近地面的气候特征，如该地的干湿状况、气候季节性指标和大气中CO₂的增加等；②下垫面表面特征，如土地利用、植被覆盖度、土壤类型和地势地貌等；③人类活动，如人口增长速率、社会经济指数和农业灌溉等；④突发气候或人类活动事件，如火山喷发、森林火灾和大规模开荒造田等。这些因素相互交叉又相互作用，仅单独考虑其中几个因素对参数的直接联系不仅可能会忽略到真正较大作用影响Budyko方程参数的变量，同时也会对水热耦合平衡方程中参数的模拟带来不确定性。流域内土壤、地形、植被和气候之间紧密联系，这种因素之间的相互依赖降低了模型响应的范围和变异性。通过归因分析和定量模拟，明确流域内影响实际蒸散发的主要因素，对科学研究具有十分重要的意义。

发明内容

本发明所要解决的技术问题是：提供一种确定蒸散发变化主因及判别因素间耦合作用的方法，通过对较大数量的样本建立空间分布模型，综合考虑各因素之间的交互作用以及各因素影响的区域性差异，准确的分析出在长期变化中实际蒸散发变化的主要影响因素。

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

一种确定蒸散发变化主因及判别因素间耦合作用的方法，包括如下步骤：

步骤1，确定待研究流域，收集流域内的水文气象数据资料、流域下垫面资料、土壤特性资料以及植被覆盖率资料并进行预处理，上述资料均为逐日数据；

步骤2，将步骤1预处理后的逐日数据划分为逐月数据，根据流域水文年划分方法将逐月数据换算为水文年数据，基于Budyko水热耦合平衡方程计算流域逐水文年实际蒸散发量以及14个影响因素，14个影响因素包括植被覆盖率、流域地形比降、土壤可用含水量、平均暴雨深、降雨集中度、降雨变差系数、降雨季节性指数、气候干旱指数、Milly指数、PfE指数、人口状况、生产总值、有效灌溉面积和作物种植面积；

步骤3，将流域内逐水文年实际蒸散发量数据进行归一化处理，并将归一化处理后的数据转化为变量场矩阵，变量场矩阵的行、列分别代表气象站点、水文年；对变量场矩阵进行经验正交函数分解并构造相关矩阵，求解相关矩阵的特征向量矩阵，根据特征向量矩阵，采用旋转经验正交函数法得到流域内实际蒸散发量的时空分布特征，并对流域进行分区；

步骤4，在各个分区内对影响实际蒸散发的14个影响因素进行聚类分析以及相关分析，并利用多元自适应回归样条法获取14个影响因素中的主要影响因素，建立实际蒸散发量与各主要影响因素之间的回归预测模型；

步骤5，根据回归预测模型得出影响实际蒸散发量的主要影响因素，并根据回归预测模型找出存在耦合关系的主要影响因素，对存在耦合关系的主要影响因素间的交互耦合作用进行识别分析。

作为本发明的一种优选方案，步骤2所述流域逐水文年实际蒸散发量计算公式如下：

根据流域总体蓄水量变化最小的原则划分水文年，即从汛期当月开始往后12个月划分为一个水文年，并根据傅抱璞公式计算逐水文年实际蒸散发量：

其中，E为实际蒸散发量，P为降水量，E₀为潜在蒸散发量，为代表流域下垫面特征的参数。

作为本发明的一种优选方案，所述潜在蒸散发量E₀计算公式如下：

其中，R_n为辐射平衡(J/(cm²·d))；G为土壤热通量(J/(cm²·d))；T为日平均温度(℃)；γ为干湿计常数(kPa/℃)；u₂为两米高度处风速(m/s)；e_s为饱和水汽压；e_a为实际水汽压；Δ为水汽压曲线斜率(kPa/℃)。

作为本发明的一种优选方案，步骤3所述对变量场矩阵进行经验正交函数分解并构造相关矩阵，求解相关矩阵的特征向量矩阵，根据特征向量矩阵，采用旋转经验正交函数法得到流域内实际蒸散发量的时空分布特征，并对流域进行分区具体过程如下：

设变量场矩阵为A_p×n，其中，p为气象站点数，n为水文年数，该矩阵可表示为：

A＝LF

其中，L_p×m为因子载荷矩阵，F_m×n为因子矩阵，利用经验正交函数法提取出m个主分量，各主分量均是均值为0、方差为1的独立变量；

因子载荷矩阵表示为：

其中，V_p×m代表空间函数，由各气象站点实际蒸散发量的变量场矩阵的m个较大的特征值对应的特征向量为列向量所构成，为m个较大的特征值开方为对角元素组成的对角阵；

F为特征向量的权重系数T_m×n的标准化形式，利用由特征向量组成的正交阵V及特征值为对角阵Λ求因子矩阵F，公式为：

式中，为矩阵的逆矩阵；故将A_p×n也表示为：A_p×n＝V_p×mT_m×n；

按照方差极大正交转动原则将F、L进行转动，使得L中各列元素平方的方差之和达到最大；

取前m个因子，使新的因子荷载元素方差达到最大值，其中，l_ij为转动后的因子荷载阵L^*的元素；

计算累计方差贡献，将累计方差贡献较多的前M(M小于m)个主分量作为代表流域内水文变量的空间分布型。

作为本发明的一种优选方案，所述步骤4具体过程如下：

(41)在分区内利用聚类分析方法对该分区内的气象站点个数进行筛选；

(42)对影响实际蒸散发量变化的14个影响因素进行相关分析，按照实际情况，从中选取相关性较大的主要影响因素建立回归预测模型；

(43)设定回归预测模型基函数的最大个数，并通过前向逐步回归过程、后向逐步回归过程以及广义交叉验证得到最终的回归预测模型。

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

本发明运用REOF和MARS方法相结合的方式得出影响实际蒸散发变化的主要影响因素，在样本足量的条件下对流域进行再分区，考虑了不同影响因素之间的交互作用，解决了大流域尺度上由于流域特性等原因造成的分析误差，对分析气候变化下的实际蒸散发的响应具有一定的价值。

附图说明

图1是本发明一种确定蒸散发变化主因及判别因素间耦合作用的方法的流程图。

图2是洛伦兹曲线模拟实测降水天数百分比和降水量百分比分布图。

图3是黄河流域REOF结果分区图。

具体实施方式

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

如图1所示，为本发明一种确定蒸散发变化主因及判别因素间耦合作用的方法的流程图，包括以下步骤：

(1)基础资料收集：收集流域内气象站点的水文气象数据资料(降水量、温度、日照时数、相对湿度、平均风速)以及各种空间数据：流域下垫面条件、土壤特性、植被覆盖率等并进行预处理。

(2)将基础数据划分为逐月值，根据流域水文年划分方法(一般将该年度洪水季节起始月份到次年洪水季节起始月份之前划分为一个水文年)，基于Budyko水热耦合平衡方程计算流域实际蒸散发年值序列以及14个影响因素。

计算流域实际蒸散发年值序列的过程如下：

首先根据流域总体蓄水量变化最小的原则划分水文年，即从当年的第一次涨水当月的第一天开始起算12个月，也就是汛期的当月开始算往后12个月，一般地区都是从3月或4月到次年的2月或3月。

其次，根据以下傅抱璞公式计算逐年实际蒸散发量：

式中：E为实际蒸散发量；P为降水量；E₀为潜在蒸散发量；为参数，代表了流域的下垫面特征。

联合国粮农组织(FAO)通过总结各国学者在不同地区的研究提出了改进后的Penman-Moneith潜在蒸散发的计算公式为：

其中，R_n为辐射平衡(J/(cm²·d))；G为土壤热通量(J/(cm²·d))；T为日平均温度(℃)；γ为干湿计常数(kPa/℃)；u₂为两米高度处风速(m/s)；e_s为饱和水汽压；e_a为实际水汽压；Δ为水汽压曲线斜率(kPa/℃)。

14个影响因素计算方式如下：

1、利用植被指数近似估算植被覆盖度，常用的植被指数为NDVI。本发明采用Yang等(2009)提出的公式来进行计算：

式中：M为植被覆盖率；NDVI_min和NDVI_max根据Yang等(2009)在中国区域分别取值为0.05和0.8。

2、流域地形比降(Relief ratio，Rr)。采用以下公式计算：

式中：Rr为流域地形比降；H_max和H_min分别为流域内最高点与最低点高程，m；L_total为流域内水文控制站点以上的最大河流长度，m。

3、土壤可用含水量(Soil profile available water capacity，PAWC，mm)，使用国际地圈生物圈计划(IGBP)在NASA地球数据网站上公开的数据，该数据来源于GlobalGridded Surfaces of Selected Soil Characteristics(IGBP-DIS)，土壤深度为0-100cm，空间分辨率为5°×5°。

4、平均暴雨深(average storm depth，ASD，mm>-1)。中国气象上规定，24小时之内，由空中降落的雨量在50mm以上的为暴雨，由于本发明研究的流域多处位于我国干旱和半干旱地区，多年平均日降水达不到暴雨程度。因此计算时先统计逐流域逐年最大日降水，然后多年平均可得该流域平均暴雨深。即通过下式计算：

式中：ASD为平均暴雨深，mm>-1；P_k,a为第a年第k天的降雨量，mm>-1；n为统计总年数。

5、采用Martin-vide(2004)提出的计算降雨集中度的方法计算研究的流域内的降雨集中度(precipitation concentration index，CI)。具体计算步骤如下：

a)以1mm为间隔，将日降雨量从0.1–0.9，1.0–1.9开始分组，到本时段降水量最大值为止(假定100mm为这一时段的最大降水量，则以100.0–100.9为最后一组)；

b)将日降雨量按其大小分别对应各个组类，统计出各组的降雨天数；

c)将统计的各个组的降雨天数分别乘以各个组类的均值得出总降雨量；

d)分别将各组的降雨天数和降雨量累加，除以其累加总和得出累积降水天数百分比和累积降水量百分比。如图2所示。

求得累积降水量百分比Y及累积降水天数百分比X后，由降雨集中度的定义知，其分布呈负的指数分布，即符合洛伦兹曲线分布：

Y＝aX exp(bX)

式中：a及b均为系数，用最小二乘法可率定。

a与b一旦求出，则曲线下的面积S可由积分所得，即：

则CI＝2S/10000。

6、降雨变差系数Cv为多年降水序列均方差与均值之比，

7、降雨季节性指数(seasonality index of precipitation，SI)。可由Walsh和Lawler(1981)得到：

式中：SI_a为第a年的SI；R_a为第a年的降雨量，mm；X_at为第a年t月份的降雨量，mm。

8、气候干旱指数采用帕尔默干旱指数(Palmer drought severity index，PDSI)代表。

9、伴随潜在蒸发极值的降水极值发生天数的Milly指数(SI_M)。Budyko和Zubenok(1961)指出P和E₀同相变化时E将增加，反之E将减小，即水热同步性对水热平衡关系的影响(杨汉波等，2012)。SI_M按Milly(1994)提出的公式进行计算：

式中：SI_M代表降水和蒸发的Milly指数；δ_P、分别代表了降水波动幅度和平均降水强度之比、蒸发波动幅度和平均蒸发强度之比，取值在-1.0–1.0之间；DI为气候干旱指数，为多年蒸发强度与降水强度之比。

10、基于实测数据选取符合正弦分布的伴随潜在蒸发极值发生的降水极值发生的月份的个数(PfE)。根据Milly(1994)所定义：即受太阳辐射的影响，理论上一年中的降水量和蒸发量均呈正弦分布。

式中：和分别为当年平均月P(mm)和E₀(mm)；δ_P和代表意义同上；ω为变化周期，2π/ω＝1年；t为月序数。

11、各个省份的人口状况(Pop，万人)、生产总值(GDP，亿元)、有效灌溉面积(IA，千公顷)和作物种植面积(CA，千公顷)。人口是一个内容复杂、综合多种社会关系的社会实体，一切社会活动、社会关系、社会现象和社会问题都同人口发展过程相关。GDP是国民经济核算的核心指标，也是衡量一个国家或地区总体经济状况重要指标。有效灌溉面积在一般情况下等于已经配备灌溉工程或设备的、能够进行正常灌溉的水田和水浇地面积之和，反映了当地耕地抗旱能力的指标。数据整理时首先根据我国各省份的行政区域面积计算这4个变量的密度即Pop密度、GDP密度、IA密度和CA密度，然后根据各个流域占各省份的面积比计算各个流域的Pop、GDP、IA和CA。

(3)分析流域内实际蒸散发的序列资料，进行数据标准化处理，将变量场矩阵进行经验正交函数分解并构造相关矩阵，求其特征向量矩阵，利用旋转经验正交函数法(Rotated Empirical Orthogonal Function，REOF)得到流域内实际蒸散发的时空分布特征，并对流域进行分区。

由于在同一个分区内实际蒸散发的变化趋势较为一致，可将区域内的站点视为区域的代表站。在对年实际蒸散发的变量场矩阵进行REOF分解前需进行EOF分解过程。

31)将年总量实际蒸散发变量长矩阵进行EOF分解。利用经验正交函数(EOF)，将空间场用尽量少的正交向量场表示，有利于集中反映资料要素的主要空间分布特征。设实际蒸散发量组成的变量场矩阵为A_p×n，其中，p为样本总数(气象站点数)，n为时间序列长度(水文年数)，该矩阵可分解为空间函数V和时间函数T两部分：

A＝LF＝(VΛ^1/2)(Λ^-1/2T)

其中，L_p×m为因子载荷矩阵，F_m×n为因子矩阵，且L为A与F的相关矩阵。V是空间函数，每一列为矩阵的归一化特征向量。矩阵T_m×n表示时间函数，为特征向量权重系数。为m个较大的特征值开方为对角元素组成的对角阵。

32)利用旋转经验正交函数(rotated empirical orthogonal function，REOF)将各个站点的水文序列资料展开，按照方差极大正交转动原则进行载荷矩阵和因子矩阵旋转，使得L阵中格列元素平方的相对方差之和达到最大。

若取前m个因子，则使

达到最大值。其中l_ij为转动后的因子荷载阵L^*的元素。

REOF方法能将空间相关的地区减少到几个，用来客观识别空间型。空间场上每一个空间点对应的变量只与一个主成分存在高相关，因此可以用REOF较为客观的进行分区，突出局部特征。

表1实际蒸散发年总量前11个主分量旋转后的方差贡献

由于在同一个分区内实际蒸散发的变化趋势较为一致，可将区域内的站点视为区域的代表站。

旋转经验正交函数(REOF)将空间场分解为正交函数的线性组合，构成位数很少的互不相关的典型模态，每个典型模态都含有尽量多的原始场信息，采用EOF和REOF进行客观分区进一步突出了局部特征。

如表1所示，前6个旋转荷载向量的累计方差贡献超过70％，可选区荷载绝对值大于8％作为区划标准进行分区，由此可得到黄河流域年实际蒸散发的六个主要分区，如图3所示。

(4)在各分区内对影响实际蒸散发的各个影响因素进行聚类分析以及相关分析，为接下来的MARS模型选取主要影响因素提供基础。

(5)利用多元自适应回归样条法(Multivariate Adaptive Regression Splines，MARS)来建立实际蒸散发量与各主要影响因素之间的回归关系的预测模型。具体过程如下：

51)在分区内利用聚类分析对样本点个数进行筛选；

52)对影响实际蒸散发变化的主要因素进行相关分析，以选取相关性较大的几个主要影响因子建立模型进行分析；

在MARS方法中，方程f(x)的样条逼近形式可看作是以下基函数的线性组合：

在这里，为输出变量的预测值；a_h为第h个样条函数的系数；b_h(x)为第h个样条函数，每个基函数代表依赖变量给定的区域；h是模型中含有的样条函数的数目。其中基函数b_h(x)(h＝1,2,…,H)包括：(1)常数。该式中含有且仅含有一个截距，值得注意的是这个常数代表了其他基函数的可能截距；(2)铰链函数(或hockey>

MARS算法的主要过程包括一个前向迭代过程和一个后向迭代过程。前向过程无目的地对训练数据进行迭代分割拟合建模，得到的基函数个数也比较多，精度虽然较高但预测能力较低；而后向过程则会选择性的删除一些基函数，使最终得到的模型具有最好的预测效果。

Friedman(1991)提供了以下推演步骤：

a、考虑到基函数个数过多不仅会影响最终模型且会造成过量计算，首先需设定变量个数的2倍作为基函数的最大个数H_max。模型的初始基函数集只包含一个常数项基函数，每次迭代都会产生两个新的基函数。

b、前向逐步回归过程

一旦确定了D个基函数组，对于输入变量来说可形成新的基函数，通过逐次迭代过程逐渐提高模型的精度。最终基于基函数组和新的基函数找到样条估算的最优拟合模型。产生的新的基函数组合加到原来的基函数上构成新组合。前向回归过程迭代计算直到达到最大基函数个数D_max或模型拟合结果相对均方误差小于0.001为止。

c、后向逐步回归过程

由前向回归过程产生过拟合会增加模型的复杂度，因此后向逐步回归过程也叫后向剪枝过程(即删除前向过程中的多余基函数)。最初产生的那些基函数可能对最终模型的贡献较小或无贡献，这些基函数的作用仅仅是产生后继的基函数。因此MARS的前向逐步回归过程允许构造大量的基函数，正因如此，后向逐步回归过程对MARS算法显得尤为重要。

对于前向过程选定的基函数组，在后向建模过程中每次循环删除一个基函数，基函数B₁(x)＝1作为所有基函数的截距部分是彻底保留到最终模型中的。后向过程是通过首先删除基函数组中最不相关的基函数入手，逐个删除至B₁(x)＝1。

d、广义交叉验证(GCV)

各个子模型的预测误差通过广义交叉验证准则(GCV)进行，在子模型集中取出预测误差最小的子模型作为最终的输出模型。

(6)利用建立的预测模型分析对比得出影响实际蒸散发的主要影响因素，并对关键要素间的交互耦合作用进行识别分析，探讨其内在物理机制。

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