基于Pareto解集的治沙区非线性分式多目标水土资源优化配置模型

郭文贤 郭 萍 潘 琦 罗 彪

(中国农业大学 水利与土木工程学院，北京 100083)

土地沙漠化是干旱、半干旱及部分半湿润地区由于人地关系不协调所造成的土地退化，以风沙活动为主要标志[1]。全球土地沙漠化面积为3 600×104km2，对全球24%的陆地面积和14%的人口存在潜在影响[2]。我国是世界上土地沙漠化危害最严重的国家之一[3]。截至2014年，我国土地沙漠化总面积达172.12×104km2，占国土总面积的17.93%[4]。沙漠化地区往往具有水资源匮乏、土壤保水能力差、地表植被稀疏等特征，其脆弱的生态环境严重制约了区域经济发展和人民生活质量。沙漠化治理是指由水(文)、土(壤)、气(候)和生(物)4个主要要素组成的生态系统退化之后在人工干预下的整体恢复过程[5]。常见的沙漠化治理措施主要有机械固沙、化学固沙和生物固沙，其中，机械固沙和化学固沙分别通过设置障碍物和改变沙地土壤孔隙结构来稳固流沙，生物固沙则通过土壤-植物-大气连续体改善沙漠化地区生态环境，从而保证沙漠化治理的可持续性，在实际工程中普遍采用。国内外对于生物固沙的研究多从固沙植物特性[6-9]、固沙林效益[10-11]、固沙林体系建设[12]、固沙过程中水资源的利用[13-14]等方面展开探索，取得了一定成果。但现有的研究大多只聚焦于治沙过程中的单一方面，较少从生态系统整体上对防沙治沙进行综合分析和治理。因此，如何充分考虑治沙过程中水、土地、生物等要素之间的关系和作用，从生态系统可持续发展的角度进行沙漠化防治，是当前亟待解决的问题。

依据决策主体的意愿，农业水土资源优化配置模型可分为提高灌溉效率的模型[15-16]、水土资源决策系统开发模型[17]和调整农业结构的模型[18]等三类。Ren等[19]构建的考虑地下水和地表水联合运用的模糊分式规划模型在优化灌水的同时提高了用水效率；Guo等[20]构建了考虑节水和质量保证的不确定规划模型来提高河西走廊制种玉米的灌溉效率；Zhang等[21]开发了集决策支持、信息管理和检索于一体的灌区水资源优化配置与灌溉决策支持系统；张帆等[22]构建了分式两阶段随机规划模型来对甘肃省民勤县作物进行种植结构优化。模型的求解方法也从最小偏差法[23]、模糊协调法[24]、理想点法[25]等传统方法逐渐扩展到智能算法[26-27]。对于复杂的多目标问题，智能算法具有更加强大的全局寻优能力。随着对智能算法的不断改进，蚁群算法[28]、粒子群算法[29]、遗传算法[30]等智能算法在求解多目标优化问题方面得到广泛应用，Guo等[31]将遗传算法和逼近理想解排序法相结合求解了非线性多层次多目标优化模型；Wang等[32]采用并列选择遗传算法求解了分布式AquaCrop—非线性多目标模糊相关机会规划模型。沙漠化治理问题与农业水土资源优化配置问题有着极高的相似性，都包含了种植结构、系统效益、用水效率和植物耗水量等指标的优化问题，然而当前水土资源优化配置在治沙过程中的研究还存在一定的空白。

基于以上问题，本研究构建了非线性分式多目标水土资源优化配置模型，模型考虑了治沙过程中水资源、土地资源和生物资源之间的作用与制约关系，将治沙过程中水、土地、生物要素联合配置，以期在提高沙漠化治理效益的同时兼顾水资源利用效率，采用多目标遗传算法求解模型，以讨论不同目标之间的博弈过程，最终得到一系列优化配置方案，为沙漠化治理的管理者和决策者提供丰富有效的理论依据和科学指导。

1 模型建立与求解

1.1 非线性分式规划原理

管理者进行决策时，为衡量系统的效率，通常会进行比率优化，当构建的模型中目标为物理或经济数量的分数，我们将这种模型形式称为分式规划问题，分式规划可以很好的反映研究对象的效率问题，同时也体现了分子分母2个模型之间的均衡[22]，当分式规划的子函数存在非线性关系时，该分式规划为非线性分式规划，典型的非线性分式规划表达式如下：

(1)

式中：f(x)为目标函数；g(x)、h(x)为关于x的函数，g(x)和h(x)至少存在一个是关于x的非线性函数；A为m×n的矩阵；x和b分别是n维和m维的列向量。

1.2 非线性分式多目标水土资源优化配置模型

人工治沙林的生态效益体现在固碳制氧、改善小气候、保育土壤、降低风速、保护生物多样性、改变沙地水分状况等方面。治沙初期，管理和决策者要考虑治沙林种植植物存活的问题，植物存活数量占种植总数量的比例即为植物存活率。谢菲等[33]通过对沙生植物不同梯度的水分处理，得到植物存活率与植物用水量呈二次函数关系的结论，本研究引用此二次函数模型来表征植物存活率。植物存活率和累积生物量结合体现为治沙区固碳制氧总价值，这一指标在治沙初期能够很好地描述治沙林生态效益，故本研究选择治沙初期固碳制氧总价值反映治沙林的生态效益。然而，水资源是制约沙漠化地区生态发展的主要因素，干旱区沙漠治理不仅要考虑生态效益的最大化，还要考虑有限水资源的利用效率，本研究以植物存活数量和总用水量的比值表征治沙区水资源利用效率。国内外关于水土资源优化配置的研究成果丰富，但大多将水资源和土地资源分开研究，水土资源联合优化配置较少，治沙区植物配水量和种植面积相互影响，如梭梭种植比例的增大会导致配水总量的增大，配水量减小时管理者会倾向于种植耗水量更小的花棒，因此本研究将植物配水量和种植面积同时设定为决策变量。

为对治沙区水土资源进行优化配置，采用多目标非线性分式规划方法构建水土资源优化模型，目标函数为最大化固碳制氧总价值和最大化用水效率，约束条件有可供水量约束、最小需水量约束、面积约束和非负约束。模型涉及到的相关参数和变量意义见表1。

本模型中的目标函数如下：

1)最大化固碳制氧总价值

(2)

Vi=1.63×Bi×0.272 9×10-6×250+ 1.2×Bi×10-6×369.7

(3)

2)最大化治沙区用水效率

(4)

约束条件：

1)可供水量约束

(5)

2)最小需水量约束

0.1mij+P0j≥Wij,min

(6)

3)面积约束

(7)

Ai,min≤Ai

(8)

4)非负约束

处女作《回到明朝当王爷》横扫网络，囊括多项年终大奖；《大争之世》《步步生莲》入围台湾图书借阅榜前二十名；《锦衣夜行》荣获台湾金石堂销量第一名。

mij≥0

(9)

Ai≥0

(10)

(11)

1.3 模型求解

本研究建立的非线性分式多目标优化模型的2个目标相互冲突，不存在能使2个目标同时达到最优的解，但存在这样的一个解：进一步优化一个目标时，必然存在另一个目标的劣化，这样的解称为非劣最优解。对于一个多目标优化问题可能存在多个非劣最优解，即多目标优化问题的非劣解集。遗传算法作为一种鲁棒的全局优化算法，提供了一种求解非线性、多模型、多目标等复杂优化问题的通用框架，并且已应用于许多水量和水质管理问题[34]。模型采用带精英策略的并列选择法进行求解，采用实数编码的遗传算法，相较于二进制编码的遗传算法，实数编码的遗传算法具有收敛速度更快、精度更高等优势。求解过程中的每一步算子结束后都进行验证，确保个体在可行域内，加快了模型的收敛速度。模型求解思路如图1所示，求解步骤的具体描述如下：

表1 本研究中相关参数意义说明Table 1 Description of related parameters in this studynni

图1 遗传算法流程Fig.1 Flow chart of genetic algorithm

1)产生初始种群。初始化运行参数，采用英国谢菲尔德大学遗传算法工具箱在解空间内生成随机分布的种群，筛选出满足约束条件的个体，种群个体数为100，将其作为初始种群。

2)根据适应度函数选择。本研究选择并列选择法[35]求解算法，构建的水土资源优化配置模型有2个目标，故将父代种群中的个体均等的划分为2个种群，种群1以固碳制氧总价值为目标函数，种群2以用水效率为目标函数，2个目标函数在各自种群独立的进行运算，以目标函数值作为适应度，根据得到的2个种群的适应度，以0.9为代沟选择出适应度高的个体，合并2个种群并去除不满足约束条件的个体，得到选择后的种群。

3)交叉和变异。将选择后的种群内的个体随机搭配成对作为2个父代，对2个父代进行整体算术杂交，得到重组种群，将交叉重组得到的种群以0.05的概率进行变异，本研究采用实值变异算子，保证了种群的多样性，将交叉变异之后得到的种群作为新一代的种群。

4)终止条件。设置最大迭代次数为100，当迭代次数超过最大迭代次数时，终止运行，记录得到的最终种群。判断最终种群是否使Pareto前沿面稳定，若前沿面稳定，结束程序，反之重新进入循环。

2 案例研究

2.1 研究区概况

治沙区S共有面积500 hm2，地处干旱半干旱过渡区，属于温带荒漠干旱气候区，干燥少雨，蒸发强烈，年蒸发量为1 800～2 800 mm，年蒸发量为年降雨量的15倍以上，有效降雨为0[36]，水资源十分匮乏，存在较为严重的荒漠化现象，风沙流动强烈，生态环境极其脆弱，仅仅依靠天然沙生植物远不足以缓解生态危机，人工种植沙生植物是遏制沙化面积扩大、改善当地生态环境的重要途径。区域内治沙林建设主要存在以下问题：1)水资源极度紧缺，区域内恶劣的自然环境导致水资源总量不足，在满足生活生产用水后可用于治沙林建设的水资源进一步减少；2)需要人工灌水保证治沙植物的存活，区域内土壤为风沙土，养分特征为平氮、少磷、缺钾、有机质含量低，土壤贫瘠，持水能力差，不利于植物生长，无人工灌水时治沙植物无法存活，无法达到理想的生态恢复效果。面对以上问题，该区域治沙管理者开始重点关注如何兼顾治沙林生态效益和水资源利用效率，因此，上节中提出的模型在治沙区水土资源优化配置研究中意义重大。研究表明花棒，梭梭，沙拐枣均具有根系发达、耐旱保水、生长迅速的特点，为治沙优势植物[37]，故治沙区选择花棒，梭梭，沙拐枣为治沙植物。

2.2 相关参数的确定

花棒、梭梭、沙拐枣根据相关试验，获得其灌水量与存活率之间的函数关系[33]。3种植物的主要生长季为4—10月，种植密度为2 000株/hm2，生长季内的每个月均需灌溉，生长季内可供水量为3 500 m3/hm2。3种植物的腾发量根据相关试验获得[38-39]，考虑到沙漠区植物的抗旱性，植物最小需水量由植物腾发量乘以折减系数0.5确定，模型所需的植物生物量数据、最小种植面积、存活率相关关系等参数已有相关研究[40-43]，具体数据见表2和表3。

3 结果与分析

3.1 水土资源优化结果

图3为非劣集的植物配水量。可见：随着目标倾向于提高固碳制氧总价值(方案编号减小)，花棒、梭梭、沙拐枣全生长季配水量逐渐增大，这是因为配水量直接影响植物的存活率，从而影响植物的固碳制氧价值。图4为3种植物固碳制氧价值对配水量的响应曲线。可见：3种植物的单位面积固碳制氧价值都随着配水量的增加而增大，但固碳制氧价值大小差异明显，从大到小依次为梭梭、沙拐枣、花棒；曲线的斜率代表响应速率，3种植物的响应速率都随着配水量的增加而减小,模型中3种植物灌水量与存活率之间都为二次函数关系，因此3种植物对配水量的响应存在相似的变化规律；植物单株生物量只影响植物固碳制氧价值的大小，对响应趋势和响应速率没有影响，3种植物的生物量关系为梭梭>沙拐枣>>花棒，说明生物量会使同一配水量下的植物固碳制氧价值存在差异。因此当决策者倾向于追求固碳制氧总价值时，建议增加系统配水量，同时可以通过调整植物种植比例来提高梭梭的配水比例，限制花棒配水比例，以保证治沙林生态效益的提升。

表2 最小种植面积、单株生物量及灌溉水量和存活率相关关系的参数Table 2 Minimum planting area, biomass and correlation coefficient between irrigation amount and survival rate

表3 4—10月植物最小需水量Table 3 Minimum water requirement of plants from April to October mm

图2 用水效率和固碳制氧总价值博弈下初始种群、最终种群与Pareto解集的表现Fig.2 Initial population, final population and Pareto solution set under the trade-off between water use efficiency and total value of carbon sequestration and oxygen production

图3 非劣方案集配水量Fig.3 Water distribution of noninferior scheme set

图4 植物配水量与固碳制氧价值相关性Fig.4 Correlation between plant water distribution and income from carbon sequestration and oxygen

本研究将用水效率定义为植物总成活数量和总配水量的比值，即灌水量与存活率函数关系的导数。可见：随着单位面积配水量的增加，3种植物的用水效率逐渐减小。花棒用水效率从0.84株/m3下降至0.54株/m3，梭梭用水效率从0.69株/m3下降至0.45株/m3，沙拐枣用水效率从0.60株/m3下降至0.42株/m3；同时，在配水量相同的情况下，花棒的用水效率显著高于梭梭和沙拐枣。随着配水量增加，函数导数逐渐减小，用水效率逐渐降低，模型若想获得较大的用水效率，将倾向于减少植物单位面积配水量。因此，当决策者倾向于追求用水效率时，应减小配水量，同时提高花棒的配水比例，以保障对水资源的高度节约。

图5 植物配水量与用水效率相关性Fig.5 Correlation between plant water distribution and water use efficiency

非劣集植物种植面积配置结果见图6。可见：当模型倾向固碳制氧总价值最大时，花棒和沙拐枣配置面积接近最小配置面积，这是因为梭梭的单株生物量远大于花棒和沙枣拐，生物量越大，固碳制氧总价值越大；当模型倾向用水效率最大时，花棒为模型面积配置中的优势植物，由谢菲等[33]的研究可知，在相同灌水量处理下，花棒的存活率大于梭梭和沙拐枣，3种植物种植密度相同，相同面积下花棒存活数量最多，存活数量和配水量的比值最大，即花棒的用水效率最大；为保证用水效率最大化，模型更倾向于减小梭梭、沙拐枣种植面积，增加花棒种植面积；无论优先考虑固碳制氧总价值目标还是优先考虑用水效率目标，沙拐枣在面积配置上均不占优势，这是因为沙枣拐的最小需水量远高于花棒和梭梭，分别是花棒的1.40倍和梭梭的1.46倍，且与梭梭相比，沙拐枣单株生物量不占优势。在对有限的土地资源进行配置时，模型更倾向于降低沙拐枣的配置面积，以保证花棒、梭梭的种植面积处于较高水平，从而获得更大的固碳制氧总价值以及更高的用水效率。然而，治沙林植物越单一，抗病虫害能力越弱，越容易退化。梭梭、花棒和沙拐枣从形态特征来看分别属于小乔木、半灌木和灌木，种植多种形态的植物增加了治沙林物种丰富度和组成多样性，同时为其他动植物提供了多元的生存空间，有利于增强生态系统的稳定性，降低治沙林退化风险，能更好的满足治沙区不同阶段的植被演替，廖空太等[44]通过对治沙林适应性和抗逆性的研究也表明梭梭、花棒和沙拐枣混交种植优于梭梭和花棒混交种植。因此建议决策者在对植物种植面积进行规划配置时仅将沙拐枣作为丰富物种多样性的植物，给予其允许的最小配置面积。

图6 非劣方案集种植面积配置Fig.6 Planting area configuration of non-inferior scheme

3.2 Pareto解集和目标博弈过程

本研究将水资源和土地资源同时优化，进一步扩大了寻优范围，单一优化水资源或者土地资源时方案选择范围大大缩小，得到的非劣解数量也低于水土资源同时优化时的非劣解数量，同时优化水土资源能够保证治沙区“地尽其力，水尽其用”。治沙区固碳制氧总价值和用水效率2个目标均为越大越优，在Pareto前沿面上，非劣配置方案的一个目标值增大必然会引起另一个目标值的减小，从第1号方案至111号方案，模型优化结果由倾向固碳制氧总价值最大逐步过渡到倾向用水效率最大，灌溉用水效率从0.49株/m3增至0.80株/m3，而固碳制氧价值从4.98万元下降至1.38万元。对比发现，不存在一种方案使得优化的2个目标值同时达到最大，即2个目标之间存在较为激烈的博弈关系。为量化固碳制氧总价值目标和用水效率目标之间的博弈关系，本研究引入变化率CR。

(12)

式中：CR为用水效率和固碳制氧总价值博弈的变化率;f2，m+1和f2，m分别为为第m+1和第m个方案时用水效率目标值;f1，m+1和f1，m分别为为第m+1和第m个方案时固碳制氧总价值目标值。

图7为变化率实际值和变化率的拟合曲线。可见：方案1～21中，变化率拟合值大于1，边际效应大，即固碳制氧总价值的减小对于用水效率的提升贡献很大，但此时用水效率小，不适合水资源匮乏的治沙区；方案104～111中，变化率拟合值小于0.1，边际效应小，即固碳制氧总价值的减小对于用水效率的提升贡献很小，若选择该范围内的方案，治沙区生态恢复速度无法保障；方案22～103在保障生态恢复速度的同时兼顾了用水效率，是较为理想的选择范围。

图7 用水效率和固碳制氧总价值博弈的变化率及其拟合曲线Fig.7 Change rate and fitting curve of the trade-off between water use efficiency and total value of carbon sequestration and oxygen production

3.3 典型方案的对比

为具体分析不同目标倾向下植物的水土资源配置情况，选取第1、55、111号方案为代表进行分析，3组方案分别代表固碳制氧总价值目标权重最大时、2个目标权重处于平衡时和用水效率目标权重最大时的最优解，3个典型方案的固碳制氧总价值依次为4.93、4.68和1.38万元，用水效率依次为0.49、0.63 和0.79株/m3，图8为典型方案配水量和植物种植面积。当追求固碳制氧总价值最大时(1号方案)，花棒、梭梭和沙拐枣的配水量得到充分满足，面积配置优先满足单位面积固碳制氧价值更高的梭梭，花棒和沙拐枣种植面积均为最小种植面积；当追求用水效率最大时(111号方案)，可供水量不足，面积配置优先满足用水效率更高的花棒，梭梭和沙拐枣种植面积接近各自最小种植面积，总配水量相较于55号方案减少了27.08%，总配水量和梭梭种植面积同时减小导致固碳制氧总价值减小了70.29%，用水效率增加了26.11%。对比以上3个典型方案，方案1虽然能得到最大的固碳制氧价值，但配水量远高于另外2种方案，对水资源造成一定程度的浪费，难以适用于水资源匮乏、取水不便的治沙区；方案111用水效率最大、配水量最小，但固碳制氧总价值远小于方案1和方案55的固碳制氧总价值，很有可能造成前期投入巨大但生态恢复成果不佳的情况，难以完成治沙区生态建设的目标；55号方案的固碳制氧总价值略小于1号方案，用水效率在3种方案中居中，即在少量减少固碳制氧总价值的情况下用水效率提升明显；考虑到治沙区需要保证一定的治沙效益，又存在水资源极度匮乏的问题，建议决策者在55号方案附近选择决策结果。

图8 不同典型方案下的配水量和种植面积Fig.8 Water distribution and acreage under different typical schemes

3种典型方案下植物对固碳制氧总价值的贡献和用水效率的贡献见表4。可以看出，3种方案下植物对固碳制氧总价值的贡献从大到小排序均为梭梭>>沙拐枣>花棒，这是因为植物单株生物量关系为梭梭>沙拐枣>>花棒，方案111虽然增加了花棒的种植面积，但梭梭最小种植面积产生的固碳制氧价值仍大于花棒最大种植面积产生的固碳制氧价值；3种方案下植物对用水效率的贡献为花棒略大于梭梭和沙拐枣，灌水量和成活率的函数关系决定了植物用水效率，从2个目标的角度看，治沙优势物种分别为梭梭和花棒，与3.1小节的结论一致。因此建议决策者优先保证梭梭的种植面积，以获得更大的固碳制氧总价值，提升花棒的成活率，当配水量对花棒成活率提升不明显时可以考虑对花棒增加施肥、病虫害防治等农业措施，以保证治沙区用水效率。

3.4 敏感性分析

沙漠化地区水资源问题突出，水是沙漠化地区修复植被、恢复生态的关键因素。为测试模型对水资源的敏感性，选择9种供水量情景对模型进行求解，求解得到的Pareto前沿面如图9所示，可以看出，供水量对Pareto前沿面存在明显的影响，供水量越小，Pareto前沿面越接近原点，这意味着随着可供水量减小，固碳制氧总价值和用水效率都相应减小；当供水量超过3 500 m3/hm2时，Pareto前沿面不再向远离原点方向移动，即超过3 500 m3/hm2的供水量存在水资源浪费的问题，因为本研究基于灌水量和存活率的二次函数关系，当灌水量增加到一定程度时，植物存活率不再增大，表现在Pareto前沿面上即为前沿面不再移动。

表4 不同典型方案下植物对固碳制氧总价值(f1)和用水效率(f2)的贡献Table 4 Contribution of plants to the total value of carbon sequestration and oxygen production (f1) and water efficiency (f2) under different typical scheme

图9 不同供水量情形下Pareto前沿面Fig.9 Pareto fronts under different water supplies

不同用水量情景下最大固碳制氧总价值存在较大差异，而最大用水效率相差较小，都在0.78～0.80株/m3之间，说明可供水量对固碳制氧总价值的影响大于对用水效率的影响。为量化供水量的变化对固碳制氧总价值的影响，选取上文55号典型方案(目标处于平衡倾向)的用水效率附近不同供水量的最优解，得到目标倾向处于平衡时固碳制氧总价值对不同供水量的响应关系，当供水量小于3 000 m3/hm2时，单位面积上(1 hm2)增加1 m3供水量，固碳制氧总价值增加20.76元，供水量的增大对固碳制氧总价值提升明显，供水量超过3 000 m3/hm2时，单位面积上(1 hm2)增加1 m3供水量，固碳制氧总价值只增加0.84元，供水量的增大对固碳制氧总价值的提升影响微弱，即当目标倾向处于平衡时，最佳供水量为3 000 m3/hm2，最佳供水量能在获得理想的固碳制氧总价值的同时充分利用水资源。

4 结 论

本研究将水土资源优化配置方法应用于防沙治沙，构建了治沙区非线性分式多目标水土资源优化配置模型，采用遗传算法求解该模型，主要结论如下：

1)当决策偏好倾向于固碳制氧总价值目标时，植物配水量会向灌水上限靠近，当决策偏好从固碳制氧价值最大向用水效率最大转变时，植物单位面积配水量呈减少趋势。无论目标倾向如何变化，沙拐枣面积配置均不占优势。

2)水资源和土地资源同时优化扩大了寻优范围，可为决策者提供更加丰富的Pareto解集，Pareto解集可以体现目标之间丰富的博弈信息，通过量化2个目标之间的博弈信息可以筛选出能同时兼顾生态效益和用水效率的方案选择范围。

3)55号方案同时兼顾了治沙区的生态恢复和水资源利用效率，是更合理的选择。从固碳制氧总价值角度看，梭梭为治沙优势物种，从用水效率角度分析，花棒为优势物种。

4)超过3 500 m3/hm2的供水量存在水资源浪费现象；当目标处于平衡倾向时，治沙区最佳供水量为3 000 m3/hm2。

综上所述，本研究的优点在于将水土资源优化配置与防沙治沙相结合，可以为治沙区水资源和土地资源配置提供决策依据。但在实际治沙过程中，基础数据可能存在较大的不确定性，不同生长季的植物对水分的敏感程度也存在差异，在今后的研究中将会继续对模型进行改进，增加基础数据的不确定性和不同生长季配水的优先序，使之适用于不确定情形下治沙区水土资源优化配置。