唐晓宇，何 英，彭 亮，王 俊

(1.新疆农业大学 水利与土木工程学院，乌鲁木齐 830052；2.新疆水利水电勘测设计研究院，乌鲁木齐 830000)

0 引 言

在内陆干旱地区，水资源十分短缺，降雨稀少蒸发强烈，灌溉是实现作物稳产、高产的重要措施[1]。但由于水资源短缺及不同产业的用水竞争，灌溉可用水量日益减少[2]，难以保证作物充分灌溉以实现高产，人们开始从根本上寻求灌溉用水最合理的利用方式，提高灌溉用水利用率[3]。在这种情况下，学者们提出非充分灌溉，其实质是允许作物一定程度的减产，从而使作物灌溉用水量和产量达到新的平衡，是应对水资源短缺的优化措施。

采用传统算法对作物非充分灌溉制度进行优化计算时，所得到的解易陷入局部最优而不能保证达到全局最优[4]。因此，学者们提出了解决此缺点的优化算法。付强等[5]建立遗传动态规划模型(RA-GA-DP)对三江平原井灌水稻的灌溉制度进行优化计算，取得了满意的结果。王斌等[6]在对灌溉制度优化研究中，分别采用自由搜索优化算法(Free Search，FS)与传统算法对优化模型进行求解，其优化结果表明FS算法的性能更好，且提高了寻优精度。杨习清等[7]利用混合粒子群算法 (SELPSO)对景泰川一期灌区春小麦的灌溉制度进行优化研究，结果表明该算法相较于遗传算法寻优结果更佳。以上3种算法都没有涉及将灌溉日期考虑在决策变量中，难以指导农田灌溉实践。尚松浩[8]采用单纯形搜索法，以灌水日期为决策变量，很好地解决了农田灌溉实践指导问题，但单纯形搜索法作为一种局部搜索方法，计算所得解为局部最优解，存在一定的缺陷性。霍军军等[9]以灌水日期作为决策变量，将模拟技术及遗传算法应用于冬小麦灌溉制度优化中，结果表明相较于传统算法，遗传算法能得到全局最优解，且结果稳定。郄志红等[10]基于改进的NSGA-Ⅱ算法，以灌水日期和灌溉水量作为决策变量，建立了能够同时对灌水日期和灌溉水量进行优化的多目标模型，优化目标由以往单目标优化模型转变成多目标优化模型，但其考虑的约束条件过少，易导致灌水日期的选择对作物生长产生不利影响。

本文以新疆阿瓦提丰收灌区棉花灌溉制度为研究对象，以农田水量平衡模型和作物水分生产函数为理论依据，以日为时间步长，以作物的灌水时间和灌水定额作为自变量，将作物的相对产量和全生育期总灌水量作为目标函数，构建内陆干旱区作物灌溉制度多目标优化模型，采用CNSGA-Ⅱ算法对优化模型求解，最终确定了该地区棉花最优的灌水方案。

1 材料与方法

1.1 研究区概况

丰收灌区位于新疆阿克苏地区阿瓦提县的西南部，北纬40°20′～40°38′，东经80°10′～80°20′，处于亚欧大陆腹地的塔里木盆地西北边缘[11]。该地区属暖温带大陆性荒漠气候，夏季炎热，冬季寒冷，昼夜温差较大，光照充足，降雨稀少，蒸发强烈。年蒸发量1 905.2 mm，年平均气温10.5 ℃，年平均降水量47.8 mm，主要集中在6-8月，约占全年降水量的60%[12]。灌区主要种植作物为棉花、玉米、小麦等，灌区棉花生育期的灌溉制度如表1所示。

表1 丰收灌区棉花的灌溉制度Tab.1 Irrigation system of cotton in Fengshou area

1.2 农田水量平衡模型

农田水量平衡模型是用于描述旱作物在生育期任一时段内，根系层(土壤计划湿润层)内的水量平衡关系[13]，其基本表达式如下：

Wt-W0=WT+P0+M+K-E-T

(1)

式中：Wt、W0分别为根系层在时段末与时段初的储水量；WT为随计划湿润层增加而增加的水量；P0为根系层内保存的有效降水量；K为时段t内的地下水补给量；M为时段t内的灌溉水量；E、T为时段t内作物实际的蒸发量和蒸腾量。

农田中主要的水分消耗项是作物需水量，其主要受大气蒸发能力(通常以参考作物蒸发蒸腾量ET0表示)、生长状态(以作物系数Kc表示)及土壤水分条件(由土壤水分胁迫系数Ks表示)等影响。作物蒸发蒸腾量(作物需水量)可以采用作物系数法计算，即：

ET=KsETm

(2)

ETm=KcET0

(3)

式中：ETm为供水充足时的田间最大蒸发蒸腾量；K2为土壤水分胁迫系数；ET0为参考作物蒸发蒸腾量；Kc为作物系数。

Kc与作物种类、生育期和作物的群体叶面积指数等因素有关，可以将其概化为生育期内时间t的函数[14]，即：

Kc=Kcmexp{-[(t-tm)2/c2]}

(4)

式中：Kcm为生育期内最大作物系数；tm为最大作物系数所对应的时间；c为形状参数。

采用FAO推荐的方法计算Ks，根区土壤含水量用根区耗水量Dr表示[15]。Ks计算式如下：

(5)

式中：TAW为土壤总有效水量；Dr为作物根区消耗水量；RAW为土壤易利用有效水分。

TAW采用下式计算：

TAW=1 000 (θFC-θWP)Zr

(6)

式中：θFC为田间持水量；θWP为永久凋萎点含水量；Zr为根区深度。

RAW的计算式如下：

RAW=ρTAW

(7)

式中：ρ为水分胁迫发生前土壤易利用有效水分与土壤总有效水量的比值。

1.3 水分生产函数模型

1.3.1 模型的选取

作物水分生产函数表示田间水分消耗与作物产量之间的关系[16]。本研究采用Jensen乘法模型，具体表达式如下：

(8)

式中：Y为作物实际产量；Ym为作物最高产量；EM、TM为作物潜在蒸发量和蒸腾量；n为作物生育期的总阶数；i为作物生育阶段划分序号；λi为作物在第i生育阶段的缺水敏感指数。

1.3.2 参数求解

传统的回归分析法在求解水分敏感指数时，所得解不稳定，且经常出现负值[17]，为避免出现负值，影响计算的精确性，本文采取水分敏感指数累积函数进行求解[18]：

(9)

λ(Δti)=z(ti)-z(ti-1)

(10)

式中：ti为作物从某一指定日期算起的生长天数；λ(Δti)为时段Δt=ti-ti-1的水分敏感指数值；z(t)为水分敏感指数累积曲线，a、b、c为待求系数。根据试验资料[19]，拟合得阿克苏棉花a=5.0072、b=0.0497、c=1.23。

1.3.3 水分敏感指数随时间的变化过程

根据上述所求得的a、b、c系数值，可以推算阿克苏地区棉花水分敏感指数随时间的变化过程，如图1和图2所示。根据当地棉花的生育期划分，得出各阶段的缺水敏感指数：λ1=0.186 4、λ2=0.202 9、λ3=0.678、λ4=0.135 8。

图1 棉花水分敏感指数累计函数变化过程Fig.1 The change process of the cumulative function of cotton moisture sensitivity index

图2 棉花时段水分敏感指数变化过程Fig.2 Change process of water sensitive index in cotton

2 优化模型与优化算法

2.1 优化模型

在非充分灌溉条件下，以作物相对产量最大和全生育期总灌水量最小为优化目标，建立以下优化模型。

(1)目标函数。

(11)

(12)

当为作物根系提供吸水作用的土层含水量达到凋萎系数含水量时，作物相对产量为零。

(2)约束条件。

灌水日期约束：

dmin≤dj≤dmax

(13)

土壤计划湿润层深度的约束：

Hmini≤Hi≤Hmaxi

(14)

土壤含水率约束：

θwp≤θi≤θf

(15)

式中：f1为作物生育期总灌水量，mm；f2为作物相对产量；Ei、Ti为作物在第i时段的蒸发量和蒸腾量，mm；EMi、TMi为作物在第i时段潜在蒸发量和蒸腾量，mm；xj为作物在第j次灌水量，mm；dmin、dmax为作物生育期的最小和最大天数，d；dj为作物从生育期开始到第j次灌水时的天数，mm；Hmini、Hmaxi为土壤在作物第i生长阶段内最小和最大的湿润层深度，mm；θi为土壤第i阶段含水率；θwp和θf为凋萎系数和田间持水率。

2.2 优化算法

1994年Deb和Srinivas提出非支配排序遗传算法(Non-dominated Sorting Genetic Algorithm，NSGA)[20]。该算法采用非支配分层的思想，利用适应度共享策略使得准Pareto面上的个体均匀分布，解决了算法早熟收敛的问题。2002年，Deb等[21]提出了快速非支配排序算法(Non-dominated Sorting Genetic Algorithm-II，NSGA-Ⅱ)。该算法主要在原NSGA算法的基础上引入精英策略，使用拥挤度概念对非支配层级中同一级的个体进行排序，提高了算法速度。

NSGA-Ⅱ算法主要用于求解无约束多目标优化问题，而本研究是以总灌水量最小及相对产量大为目标函数的多约束多目标优化模型，所以将原算法进行改进，将Deb约束准则引入到NSGA-Ⅱ算法中，记作CNSGA-II算法。

Deb约束支配准则中每个个体根据Pareto支配、约束违反度两种信息来选择较优个体。CNSGA-Ⅱ算法中约束处理的方法：如对于一个解x，若其满足约束条件，则称该解为可行解，若不满足，则称之为不可行解。对于不可行解，一般使用约束违反值(Constraint Violation value，CV)，该值用来定量描述一个解违反约束条件的程度。对于一个解x，其值可如下表达：

(16)

其中当g(x)<0时，g(x)= -g(x),否则g(x)=0，当x满足约束条件中任何条件，即x在可行域内，则CV(x)=0；当x不完全满足约束条件时，即x不在可行域内，则CV(x)≠0，当CV值越小，x越靠近可行域。

3 模型求解及结果分析

3.1 模型参数设定

本研究以丰收灌区棉花为例，该灌区棉花的生育期从4月15日到8月20日，共138 d。已收集阿瓦提县气象站(东经80°24′，北纬40°39′，海拔1045.8m)1960年到2001年共42年的降水资料。模型参数，如：土壤初始含水量W、田间持水量θFC、作物缺水敏感指数λi，生育期潜在腾发量ETM、生育期实际腾发量ET。其中ET、ETM分别由公式(2)、(3)计算。

本研究假定棉花在生育期总共灌5次水，将每次灌水量xi与灌水间隔ti作为优化参数，灌水时间间隔设置为整数。算法中参数选取：种群规模P=100，最大进化代数g=500，运行次数20，变异概率Pm=0.1，交叉概率Pc=0.9，交叉分布指数etac=20，突变分布指数etam=80。

3.2 优化结果与分析

采用CNSGA-II算法对灌溉制度多目标优化模型进行求解。为提高模型的适用性，本文选取不同典型年的降水优化相应灌溉制度。典型年的选取方法：收集当地1960-2001年共42 a的降水资料，按照总降水量在作物全生育期的大小进行频率计算，选取1996年(生育期降水量109.5 mm)作为丰水年(频率P=10%)；选取1992年(生育期降水量52.7 mm)作为平水年(频率P=50%)；选取1973年(生育期降水量13.2 mm)作为枯水年(频率P=90%)。典型年灌水总量与相对产量的优化结果如图3～图5所示。

图3 枯水年相对产量关系优化结果Fig.3 Optimal results of the relative output relationship in dry years

图4 平水年相对产量关系优化结果Fig.4 Optimal results of the relative output relationship in the average year

图5 丰水年相对产量关系优化结果Fig.5 Optimal results of the relationship between the relative output in high water years

由图3～图5可知，丰水年、平水年和枯水年图形中优化解的变化趋势大致相同，在一定范围内，棉花的产量都随着灌水量的增加而显著增长。但在平水年，当棉花的总灌水量达到385 mm，其相对产量增加的幅度明显减小且随着灌水量的增加而趋近于直线，即灌水所带来的效益微乎其微，而在丰水年需要总灌水量352 mm、枯水年则需总灌水量416 mm才能达到相同的效果。

在3个典型年的相对产量稳定不变时，他们的相对产量分别为：丰水年0.954、平水年0.929、枯水年0.926，说明降水对平水年与枯水年的相对产量影响不大。在选取典型年为平水年时，该模型相较于当地以往所制定的灌溉制度：在灌溉定额一致时，作物相对产量提高了2.5%～3.1%；在相对产量一致时，每公顷土地的灌溉用水节约了495～705 m3，表明对灌水日期和灌溉水量进行合理的安排能在保证作物产量稳定的情况下节约灌溉用水，对于作物生长依靠灌溉的缺水地区，该模型的合理利用能使农户在作物产量和总灌水量不变的情况下，增加作物种植面积，提高当地农户的经济收入。从平水年的优化结果中选取部分Pareto解如表2所示。

4 结 论

本研究利用农田水量平衡模型模拟作物根系层水量之间的关系，采用作物水分生产函数模型计算其相对产量，以日为步长，以作物的灌水时间和灌水定额作为自变量，将作物的相对产量和全生育期总灌水量作为目标函数，构建内陆干旱区作物灌溉制度多目标优化模型；为处理含约束的多目标优化问题，将Deb约束准则引入到NSGA-Ⅱ算法中，记作CNSGA-Ⅱ算法，并利用CNSGA-Ⅱ算法对优化模型进行求解，得到了实验地区最优的灌溉制度。

表2 棉花在灌溉制度优化下的部分Pareto解Tab.2 Part of Pareto solution of cotton under the optimization of irrigation system

通过对棉花在3种典型年的灌溉制度优化表明：在内陆干旱区作物的生长主要依靠农业灌溉，降水对其相对产量的影响不大。在典型年为平水年时，相较于当地以往的灌溉制度：在灌溉定额一致时，作物产量提高了2.5%～3.1%；在相对产量一致时，每公顷土地的灌溉用水节约了495～705 m3。

棉花灌溉制度多目标优化模型在优化灌水量的同时也将灌水日期一并优化，相较于以往只考虑单目标的优化模型，其优化结果更加合理，而灌水时间间隔设置为整数，便于指导农业灌溉实践。目前模型输入的降水数据为历史数据，具有一定的局限性，为制定适用性更强的灌溉制度，在进一步研究中将预测未来气候变化的模型与灌溉制度优化模型相结合，既能预测未来极端天气，又能动态调整灌水量与灌水日期。