考虑水资源短缺风险的不确定性水资源优化配置

付 强，王珊珊，李天霄

（东北农业大学水利与土木工程学院，哈尔滨 150030）

全球水资源需求以每年1%速度增长，未来20年还将继续加快[1]，关注水可获得性、水质和涉水风险，改善水资源管理和利用可持续性。我国农业灌溉用水短缺问题严峻，亩均水资源量仅占世界平均水平1/2[2-3]。水资源短缺突出，洪涝、干旱灾害与全球气候变化变暖情况下，可利用水资源总量无法满足农业生产需求，农业水资源供需矛盾加剧[4-5]。灌区作为农业用水主体，优化配置灌区农业灌溉水资源，提高水资源综合效益，保持水资源可持续利用具有重要意义。

水资源优化配置方法很多，灌区水资源管理系统受经济、资源等因素影响，为复杂系统动态配水过程。现实系统中水源降雨量时空分布不均匀性、可供水量、需水量、经济参数涨落不确定性，导致传统确定性优化方法无法准确求解[6-7]。以往研究将水资源配置系统简化为确定性系统，忽略农业水资源系统中不确定性，在水资源优化配置上存在局限，决策方案使用价值有限[8-9]。学者针对系统中不确定性和复杂性，提出有效的系统分析方法。Chen等将随机规划框架中引入区间参数，解决农业用水系统中区间和概率分布不确定性[10]；Li等将区间线性多目标规划应用于农业灌溉水资源配置[11]；Xie等将鲁棒性区间随机规划方法应用于区域水资源管理[12]；李晨洋等将区间两阶段模糊随机规划模型用于红兴隆灌区水资源管理[13]；Yang等将区间参数线性规划模型应用于联合使用地表水和地下水优化调度[14]，区间规划不确定性理论方法可有效处理农业系统复杂性问题。

本文在前人研究基础上，构建区间线性模型，解决灌区水资源系统中存在的不确定性。降雨量和参考作物蒸腾量是灌溉系统中两个具有相关关系的重要水文随机变量，可表征灌区供需水变化关系[15-16]。由于二者变化规律存在相关性和随机性，若仅分析单一变量变化规律和分布特征，不能反映灌区降雨量和参考作物蒸腾量的组合情况，无法科学指导灌溉规划。以往研究很少将不同情景下水资源供需短缺风险引入灌溉用水分配方案。采用Copula函数[17-18]对降雨量和参考作物蒸腾量二维联合概率分布模型量化，分析灌区水资源供需短缺风险组合。本文针对灌区水资源系统中不确定性，考虑供需关系，将区间线性规划模型和水资源风险短缺结合在各作物间优化配水。模型中以灌区总收益最大为目标函数，采用Copula函数、离散区间解决不确定性问题，更适于灌区内作物间优化配置，合理优化有限灌区农业水资源，对提高黑龙江省灌区水资源利用效率具有重要意义。

1 模型建立及求解

1.1 Copula函数-联合概率分布模拟

假设（X1，X2，…，Xd）为d个随机变量，其边缘分布函数分别为FX1，FX2，…，FXd。基于Sklar定理Copula函数是多个随机变量边缘分布联合分布函数，定义域为[0，1]，公式如下：

式中，C()为Copula函数;θ为Copula函数参数。

将灌区降雨量和参考作物蒸腾量统计变化概率和概率分布特征结合，构建水资源供需联合概率分布模型。基于此，引入Copula函数模拟其联合分布。设X和Y分别为参考作物蒸腾量、降雨量随机变量，其联合分布函数为H(x,y)，边缘分布函数分别为Fx(x)和Fy(y)，则存在联合函数Copula C，满足H(x,y)＝C(Fx(x)，Fy(y))。边缘分布模拟采用领域应用较广泛分布模型，即广义极值（Generalized extreme value，GEV）和广义正态（Generalized normal，GNO）。

Copula函数各变量间相依结构和边缘分布可独立考虑，应用较优优势。其中变量相关性函数有Frank、Clayton、Gumbel-Hougaard Copula 函 数 。通过采用非参数Kolmogorov-Simirnov（K-S）拟合度检验，离差平方和（OLS）最小准则对Copula函数进一步优选，构建和平灌区1958～2013年降雨量和参考作物蒸腾量联合概率分布模型。

采用频率法划分降雨量和参考作物蒸腾量丰枯。其中，

1.2 区间线性规划模型构建

灌区水资源配置目标是将水资源合理配置到各种作物中，确保满足作物生长所需最优水量，获得最大系统收益。但实际灌区研究中所获得数据或信息具不确定性，将不确定相关参数值表达为一个区间范围，可解决优化模型中不确定性目标函数和约束区间问题。通过区间规划（Interval parameter programming，IPP）和线性模型（Linear programming，LP），构建区间线性规划模型（Interval parameter linear programming，IPLP）表达和求解。

子模型下限：

子模型上限：

求解上限模型时考虑下限模型求解结果的约束，得到上限目标函数决策目标值和上限决策

1.3 模型求解

在优化配置中，需考虑水源可供水量、灌溉需求、作物需水量等约束，考虑不同量级降雨量和参考作物蒸腾量遭遇组合，构建区间线性优化配置模型。

模型目标函数和约束表达如下：

约束条件：

①水源可供水量约束

②灌溉需求约束

③作物需水量约束

④非负约束

式中， f±为灌溉净效益（元）； j为不同作物，其中 j=1，2，3，分别代表水稻、玉米、大豆；为作物市场单价（元·kg-1）；为灌区作物j种植面积（hm2）；YWj为灌区作物 j水分生产率（kg·m-3）；CW±为单方水价（元·m-3）；为供水目标（m3）；为灌区初期蓄水量（m3）；为灌区可利用来水量（m3）； D±i为灌区蒸发、渗漏等损失水量（m3）； EPj为灌区有效降雨量，m3·hm-2； Ri为灌区不同情景下灌溉水量（m3·hm-2）；为灌区内作物 j充分灌溉条件下最小需水量（m3）；灌区内作物 j充分灌溉条件下最大需水量（m3）；η为灌溉水利用效率。

1.4 作物灌溉水量优化

水资源有限情况下，一旦供水量未达到供水目标值，需减少灌溉水量，但影响作物产量。因此，根据作物不同生育阶段缺水对产量影响程度，判断各生育阶段需水程度，是否处于需水关键期，在生育阶段应用灌溉优化模式[19-20]采用水分生产函数每个生育阶段作物实际蒸腾量与作物最大蒸腾量比值，使产量最大化[21]，保证灌溉水量在作物各生育阶段合理配置，其中约束条件为每个生育阶段

2 案例研究

2.1 研究区域概况

黑龙江省和平灌区位于黑龙江省中部绥化市庆安县境内，东经 127°20'～127°49'，北纬 46°41'～47°04'。地区气候属温带半干旱半湿润大陆性季风气候，年平均气温1.7℃，多年平均蒸发量664.5 mm，多年平均降水量545.3 mm，主要降雨量集中在6～9月份，无霜期120～140 d。土壤以黑土和水稻土为主，灌区位居呼兰河流域上游左岸一级阶梯，灌区范围由东向西呈带状分布，主要作物为水稻、玉米、大豆。和平灌区主要取水工程包含和平渠首和安邦河渠首。

2.2 数据来源

根据《呼兰河灌区工程初期设计报告》及当地水务局提供调研数据，得到不同作物（水稻、玉米、大豆）种植面积。通过黑龙江省用水定额标准及相关文献田间试验数据[22-24]，得到3种作物单位面积最大和最小灌溉用水量（见表1）。

基本气象数据（降水量、平均风速、最高气温、最低气温、平均相对湿度、日照时间等）来自中国气象科学数据共享服务网CMDC（http：//data.cma.cn），利用FAO推荐Penman-Monteith公式[25]计算ET0，采用作物系数法[26]计算作物需水量。本文参照文献[15]确定水稻、玉米、大豆水分生产率分别为1.61、0.82、0.72；灌溉水利用效率为0.4408。

表1 和平灌区作物参数Table 1 Crops economic parameters in Heping irrigation district

3 结果与分析

3.1 边缘分布函数模拟

由图1可知，年参考作物蒸腾量变化较小，而年降雨量与ET0间变化差异显著。因此，对于灌溉水资源确定不同组合遭遇情景条件下ET0和降水有意义。应用水文频率分析方法，选用相应随机变量概率分布函数，和平灌区降雨量服从正态分布，参考作物蒸腾量服从广义极值分布。通过拟合结果判断ET0和降水量理论分布函数与经验分布函数大致相同（见图2），可假定和平灌区降雨量与ET0服从累积概率密度分布模型合理。

借助变量相关性指标分析水资源供需随机变量间相关关系，采用Pearson相关系数γ、Kendall秩相关系数τ、Spearman秩相关系数ρ度量和平灌区1958～2013年降水量和参考作物蒸腾量间相关性，得出γ=-0.108、τ=-0.204、ρ=-0.146，模拟分析得到降雨量及参考作物蒸腾量边缘分布函数，其边缘分布函数参数值见表3。可知和平灌区降雨量与ET0间存在负相关。因此，可采用Copula函数构建水资源供需联合分布模型描述和平灌区水资源供需相关特征。由Copula函数参数和Kendall秩相关系数间关系式，基于此运用Copula函数建立二维联合分布模型，分别采用常见二维对称Archimedean型（Frank、Clayton、Gumbel-Hougaard Copula函数）构建和平灌区1958～2013年降雨量和参考作物蒸腾量联合概率分布模型。采用Kolmogorov-Simirnov（K-S）检验统计量D值均小于临界值0.20056，因此3种Copula函数均可构建联合概率分布函数模型。通过离差平方和（OLS最小准则）检验Copula函数拟合度，分别为0.1723、0.1931、0.1835。其中Frank Copula函数建立联合分布模型下年降水量和参考作物蒸腾量联合分布情况拟合较好，较好模拟降水量与ET0间相依关系。考虑灌区降雨量和参考作物蒸腾量不同组合条件下联合分布概率：丰枯划分为降雨量超过某一特定值663.5 mm、不超过某一特定值571.3 mm，其边缘分布概率值分别为0.6209、0.3728；ET0超过某一特定值914 mm，不超过某一特定值852 mm，其边缘分布概率值分别为0.6395、0.3615。采用频率法分析和平灌区降雨量和ET0丰枯划分情况频率，丰枯遭遇情况分为9种。其中丰枯同步类型：P-ET0同丰，P-ET0同平，P-ET0同枯；分别为8.88%，7.54%，8.31%。丰枯异步类型：P-ET0丰平，PET0丰枯，P-ET0平枯，P-ET0平丰，P-ET0枯丰，P-ET0枯平；分别为9.13%，19.74%，8.33%，8.64%，19.92%，9.51%。降雨量和ET0丰枯异步频率远大于其丰枯同步频率，分别为75.27%，24.73%。P-ET0丰枯和P-ET0枯丰情况是9种遭遇组合中最不利于灌溉调度情况，其遭遇频率约达20%。因此在农业生产中，灌溉活动、科学调配水资源尤为必要。

图1 和平灌区1958～2013年降雨量、参考作物蒸腾量数据Fig.1 Long serial changes of precipitation and ET0for 1958-2013 years

图2 和平灌区降雨量、参考作物蒸腾量分布函数拟合Fig.2 Estimation of distribution functions of precipitation and ET0

表2 水文变量基本统计参数Table 2 Basic statistical parameters of hydrological variables (mm)

表3 边缘分布函数参数值Table 3 Parameter values of marginal distribution function (mm)

3.2 最优配水方案

由图3可知，丰枯同步情境下，最优总配水量变化较明显，在P-ET0同枯8.31%，P-ET0同平7.54%，P-ET0同丰8.88%频率下分别为[2 980，3 083]×104m3，[3 676，3 920]×104m3，[4 875，5 140]×104m3。其中，在P-ET0同丰可利用水量较大情景下，模型通过增大每种作物单位面积配水获得目标最优化，总水量分配（净灌溉水量和有效降雨量之和）达到作物最大需水量上限。而P-ET0同枯情景下为满足作物最低需水量，确保作物基本生理需求，选择加大经济效益获得目标最优化，存在补给灌溉水量比例高达70.5%。

丰枯异步情境下，频率最高两种遭遇组合P-ET0丰枯19.74%，P-ET0枯丰19.92%最优总配水量分别为[4 470，4 462]×104m3，[3 420，3 785]×104m3。P-ET0平枯8.33%，P-ET0枯平9.51%频率下最优配水量分别为[4 602，4 754]×104m3，[3 874，3 990]×104m3。P-ET0平丰 8.64%，P-ET0丰平9.13%频率下最优配水量分别为[3 318，3 625]×104m3，[3 586，3 860]×104m3。表明不同情境下有效利用有限水资源分配方案重要性。

整体看，P-ET0同枯，P-ET0枯丰，P-ET0枯平情境下，受最小需求量约束，模型将水量较为均匀分配。但P-ET0同丰，P-ET0丰平，P-ET0丰枯，总配水量相差有较大幅度差距，各类作物水量分配分化较明显，实现模型目标有利于作物获得水量较多。因为降雨量具有明显干湿特征，对作物水量配置贡献影响较大，在降雨量较大湿润时作物产量高于正常情况下，正常降雨下产量高于干旱条件。结果表明，作物对降水较敏感，通过可利用有效水资源应对不同降雨量和参考作物蒸腾量联合分布水资源短缺情况。

图3 不同情境下最优水量分配及降雨量Fig.3 Optimal water allocation schemes and precipitation under different scenarios

由图4可知，由于市场价格和水分生产率相对较高，作物种植面积及配水量为水稻＞玉米＞大豆。分析所有遭遇组合频率最大P-ET0丰枯和PET0枯丰情境对各作物水量分配，水量充足P-ET0丰枯19.74%情景下水稻、玉米、大豆配水量分别为 [3 702， 3 821]× 104m3， [359， 396] × 104m3，[409，443×104m3]，对于模型结果较为敏感，水量区间分配变化幅度较大。

各作物分配水量平均值比水量不足条件下PET0枯丰19.92%情景下高794×104m3，110×104m3，63×104m3。由于水稻种植面积过大，水稻配水必须满足最低要求，而面积占比较小的玉米和大豆，少量水量即可满足目标优化。可适当调整农作物种植结构，增加旱作物玉米、大豆种植面积，减少水稻种植面积，有效减少缺水量，提高水资源利用率。

图4 不同情境下作物最优水量分配Fig.4 Optimal water allocation for each crop under different scenarios

3.3 作物配水量优化

以耗水量最大水稻为典型作物，研究生育阶段水量优化情况，为便于定量分析水稻生育期降水量和需水量，将和平灌区水稻生育期确定为5月17日～9月20日，共127 d。其中返青期（5月17日～5月30日）、分蘖期（5月31日～7月7日）、拔节孕穗期（7月8日～7月25日）、抽穗开花期（7月26日～8月4日）、乳熟期（8月5日～8月24日）、黄熟期（8月25日～9月20日）。由于返青期稻苗较小、黄熟期自然落干，水稻蒸腾量较小，因此仅选取分蘖期、拔节期、抽穗期、乳熟期。采用相关文献[21]分析各生育阶段作物系数，本文对所有情境中3种极端遭遇组合P-ET0同丰，P-ET0同平，P-ET0同枯计算水稻分蘖期、拔节期、抽穗期、乳熟期4个生育阶段配水。在3种典型极端情境下满足各生育期净灌溉需水量时即达到ETa/ETm=1，作物从开始播种到成熟未遭受水分胁迫，产量达最高。如果作物生长发育阶段，当灌溉水量可用性达到限制时，产量下降。为避免水资源限制，考虑4个总体缺水目标ETa/ETm（0.9，0.8，0.7，0.6），得到不同极端情境水分亏缺状态时每个生育阶段作物实际蒸腾量（ETa）和最大蒸腾量（ETm）间比值。

由图5可知，分蘖期、拔节期亏缺灌溉水平比例均低于抽穗期、乳熟期，因两生育阶段水分敏感指数高于其他，处于需水关键期。为获得更高产量，分蘖期、拔节期分配水量显著高于抽穗期、乳熟期。

在P-ET0同丰情境下，ETa/ETm总体目标为0.9、0.8时，分蘖期ETa/ETm为1.0，避免缺水灌溉。另外两种目标中，分蘖期、拔节期与其他生育阶段ETa/ETm相比较大，抽穗期、乳熟期生育阶段经历亏缺灌溉。在P-ET0同平和P-ET0同枯条件下，4个总体ETa/ETm目标情况下，分蘖期无法达到充分灌溉。当处于极度缺水条件ETa/ETm总体目标为0.7、0.6时，由于3种情境下来水量年内分配不均，在4个生育阶段间水量分配间余水量供下一阶段来水量较小生育阶段，保证生育阶段实际蒸腾量与最大蒸腾量比值（ETa/ETm）高于0.5。抽穗期、乳熟期ETa/ETm接近0.5，说明抽穗期、乳熟期是应用亏缺灌溉主要生育阶段，缺水水平较高。上述优化结果表明，对于不同水分亏缺情景，分蘖期、拔节期属于作物需水关键期，作物产量较大，因此为达高产目的，选择较充分满足作物水分基本需求，而在需水非关键期（抽穗期、乳熟期），为提高灌溉水利用效率，节约水资源，则选择减少灌溉水量。P-ET0同枯，P-ET0同平，P-ET0同丰情景下配水量见图6。

图5 3个典型情景4个缺水水平下水稻生育阶段ETa/ETmFig.5 Optimized stage deficit targets（ETa/ETm）for four overall deficit targets for rice under different growth stages and scenarios

图6 3个典型情景4个缺水水平下水稻生育阶段配水量变化情况Fig.6 Water amount change for four overall deficit targets for rice under different growth stages and scenarios

3种情境下，P-ET0同丰配水量最大，P-ET0同枯反之；通过4个缺水水平对比，总体缺水目标ETa/ETm为0.9时配水量最大。P-ET0同丰情境下供水目标未达到初始目标情况下缺水水平为0.9、0.8时，其中配水量分别为2 640.5×104m3、2 422.5×104m3，略下降，分蘖期、拔节期水量稳定，抽穗期、乳熟期水量变化显著。说明在保证产量同时，通过非工程措施应用调亏灌溉优化模式，管理者可根据不同缺水情境下判断各生育阶段实际蒸腾量与最大蒸腾量比值（ETa/ETm），使灌溉水资源在作物各生育阶段合理配置，有效避免水资源浪费。而在其他情境缺水水平为0.7、0.6时配水量有显著变化，产量均大幅下降。表明对于缺水水平较严重、模型稳定性较差情况应尽量避免。不同缺水条件下管理者应优先满足作物分蘖期、拔节期用水需求，保证产量。

4 讨论

本文通过建立区间不确定性模型解决水资源系统中的不确定性，考虑不同组合联合概率，研究和平灌区作物水资源配置。以往水资源配置重点研究3种水平情境下配水方案，本文重点讨论降雨量和ET09种量级组合下，应用Copula函数定量分析灌区天然供水条件下水资源供需特征变量发生各种组合风险，为灌区水资源短缺风险分析提供技术指导。肖圆圆等仅研究高中低来水流量下水稻用水优化[27]，本文选择3种作物不同条件展开研究，可更精确制定和调整灌溉规划。通过将Stewart模型与水稻不同生育阶段结合，减少灌溉水量达到较高产量。刘银凤研究认为，水稻需水敏感期，应充分满足分蘖期、拔节期生育阶段需水量，而需水非敏感期对产量影响较小，可选择性减少灌溉水量，与本研究结果一致[28]。本研究假设水源同时向3种作物配水，与生产实际存在差异，对模型预测精度有一定影响。今后应将环境效益、灌溉方式引入灌区作物水资源配置方案研究，使水资源优化配置更符合实际，更具代表性。

5 结论

a.以黑龙江省和平灌区1958～2013年降水量P和ET0为研究对象，借助相关性指标分析水资源随机变量间相关性，两者间存在负相关。根据边缘分布函数极大似然法估计参数，使用Copula函数构建灌区自然供水条件下水资源联合分布模型，应用该模型分析相应短缺风险遭遇组合联合概率。其中Frank Copula函数可较好拟合和平灌区降水量P和ET0理论联合分布与经验联合分布，运用联合概率分布模型可较好描述不同遭遇组合事件，定量分析灌区水资源供需特征变量发生各种遭遇组合风险情境。其中P-ET0丰枯和P-ET0枯丰情况是9种遭遇组合中最不利于灌溉调度情况，其遭遇频率约达20%。

b.用Copula函数及离散区间表示水资源系统不确定性，采用区间线性规划模型优化配置各种作物水量分配，确保满足作物生长所需最优水量，获得最大系统收益。P-ET0同丰，P-ET0丰平，P-ET0丰枯总配水量较大，各类作物水量分配分化较明显。从各种作物分配水量来看，耗水作物种植面积最大，水稻配水量最高，水量区间分配变化幅度较大。

c.以水稻为典型作物优化研究各生育阶段用水量，考虑4个总体缺水目标ETa/ETm（0.9，0.8，0.7，0.6），计算3种极端遭遇组合（P-ET0同丰，P-ET0同平，P-ET0同枯）下水稻分蘖期、拔节期、抽穗期、乳熟期阶段配水，其中分蘖期、拔节期亏缺灌溉水平比例低于抽穗期、乳熟期，表明分蘖期、拔节期为需水关键期。在一定缺水条件，可据此模型在满足产量同时对各生育阶段灌溉水量做出调整，以高效利用水资源。