基于PSO-PPE模型的宁夏扬水灌区水土资源承载力综合评价

邹全程， 伍昊洋， 叶 威， 王春懿， 赵 璐，2， 崔宁博， 梁 川

(1.四川大学 水利水电学院， 四川 成都 610065； 2.农业部旱作节水农业重点开放实验室， 北京 100081)

1 研究背景

土地资源是人类赖以生存和发展的基础性自然资源，是农业生产发展中的核心要素。通过对农业水土资源匹配特征和承载力现状的评价，可为土地规划、管理及资源优化配置提供理论依据，对区域的可持续发展具有重要作用[1]。

目前，国内外有很多学者对水土资源承载力进行了研究[2]，且有大量学者从不同角度，采用不同方法对水土资源承载力进行了评价。拓学森等[3]综合考虑了民勤县水土资源与生态环境的实际特点，建立了民勤县水土资源承载力系统动力学模型，用于预测未来水土资源承载力的变化趋势；封志明等[4]综合人粮、人水关系定量评价分析了雄安新区的水土资源承载力，提出了若干解决建议；张晶等[5]应用GIS技术，研究了宁夏平原在天然和人工不同条件下的水土资源特征；李天霄等[6]采用实码加速遗传算法优化投影指标函数，建立水土资源承载力DPSIR模型指标，分析了齐齐哈尔市市区的水土资源承载力等级情况。而投影寻踪模型应用到水土资源承载力评价的研究越来越多[7-9]，且随着算法不断改进，粒子群优化算法已经逐渐被更多学者用于不同地区的水土资源承载力评价[10-11]，基于粒子群优化算法的投影寻踪模型已经是目前应用较成熟的评价方法之一。

自20个世纪70年代以来，为逐步解决宁夏回族自治区中部干旱带日趋严重的干旱缺水状况，该地区已经兴建了红寺堡、固海、固海扩灌和盐环定等扬水工程。为实现扬水灌区水土资源科学的可持续利用和管理，对灌区进行农业水土资源匹配，承载力研究和农业干旱风险的评价十分重要。基于此，本文构建了农业水土资源承载力综合评价指标体系，采用庄妍[8]研究的粒子群优化投影寻踪模型对宁夏扬水灌区各个分区农业水土资源承载力进行了评价，以期为宁夏扬水灌区水土资源可持续利用提供科学依据。

2 材料与方法

2.1 研究区域概况

宁夏扬水灌区是自20世纪70年代发展的大型农业灌溉区，该区域属于中温带干旱区，太阳辐射强，大陆性气候明显，干旱少雨且强烈蒸发，致使区域土壤较干旱。

为研究灌区农业水土资源的可持续利用情况，本文选择比较完整的红寺堡扬水灌区(以下简称为红寺堡灌区)、固海扬水灌区和固海扩灌灌区(合并简称为固海固扩灌区)作为研究区域。为方便数据收集与研究，将灌区基本按所涉及县域进行分区，其中红寺堡灌区分为H1、H2、H3、H4四个分区，所涉及县域为利通区、红寺堡区、同心县、中宁县；固海固扩灌区分为G1、G2、G3、G4、G5、G6六个分区，所涉及县域为红寺堡区、同心县、原州区、沙坡头区、中宁县、海原县，图1为研究区域分区概况图。

图1 扬水灌区分区概况图

2.2 评价指标体系建立

水土资源存在相互制约的特点，本文在对宁夏扬水灌区土地资源及各影响因素综合分析的基础上，参考国内学者相关研究成果[13-16]，以灌区水土资源承载力评价为总目标A，水土资源利用评价B1、社会发展评价B2、经济技术评价B3、生态环境评价B4为子目标，构建了包括20个指标的分层评价体系，如图2所示。

其中，指标数据来源于《2015宁夏水利统计公报》、《2015宁夏水资源公报》、2016年《宁夏统计年鉴》、《宁夏大型扬黄灌区用水现状调查报告》等。根据指标性质，若指标xi值越大，所反映的农业水土资源承载力越高，称为正向指标，用(+)标识；反之，则称为负向指标，用(-)标识。

2.3 指标分级标准确立

考虑红寺堡灌区、固海固扩灌区的特点，结合文献[8,17-18]，将灌区农业水土资源承载力评价结果划分为5个等级，分别是高(Ⅰ)、较高 (Ⅱ)、中等 (Ⅲ)、 较低(Ⅳ)和低(Ⅴ)。Ⅰ～Ⅲ级代表区域生态功能完善，水土资源与社会经济发展系统很稳定，处于可持续状态；Ⅳ级代表生态功能相对完善，水土资源与社会经济平衡发展的状态被破坏，系统不稳定；Ⅴ级代表区域生态环境受到破坏，水土资源承载力已接近饱和。具体指标见表1。

注：A、B1～B4、C1～C20分别为评价指标体系内的总目标、4个子目标、20个指标。

指标类型Ⅰ高Ⅱ较高Ⅲ中等Ⅳ较低Ⅴ低C1/%正(+)>55.141.4～55.127.6～41.413.8～27.60.0～13.8C2负(-)0.0～0.10.1～0.30.3～0.40.4～0.6>0.6C3正(+)>0.90.6～0.90.4～0.60.2～0.40.0～0.2C4正(+)>6.04.5～6.03.0～4.51.5～3.00.0～1.5C5/104(m3·hm-2)负(-)0.0～0.10.1～0.30.3～0.40.4～0.5>0.5C6/104 t正(+)>7.65.7～7.63.8～5.71.9～3.80.0～1.9C7/104 m3正(+)>0.210.16～0.210.10～0.160.05～0.100.00～0.05C8/hm2正(+)>0.410.31～0.410.21～0.310.11～0.210.0～0.11C9/(人·km-2)负(-)0.0～37.137.1～74.174.1～111.2111.2～148.0>148.3C10/‰负(-)0.0～4.54.5～8.98.9～13.413.4～17.8>17.8C11/104元负(-)0.00～0.980.98～1.961.96～2.942.94～3.92>3.92C12/%正(+)>7.25.4～7.23.6～5.41.8～3.60.0～1.8C13负(-)0.0～0.10.1～0.20.2～0.30.3～0.4>0.4C14/104(元·hm-2)正(+)>3.62.7～3.61.8～2.70.9～1.80.0～0.9C15/104元正(+)>2.72.0～2.71.3～2.00.7～1.30.0～0.7C16/(m3·hm-2)负(-)0～40654065～81158115～1218012180～16230>16230C17/108 m3正(+)>0.240.18～0.240.12～0.180.06～0.120.00～0.06C18/(t·hm-2)负(-)0.00～0.60.6～1.351.35～1.951.95～2.55>2.55C19/%正(+)>150113～15075～11338～750～38C20/km2正(+)>2116～2110～165～100～5

2.4 具体建模步骤[8,11-12]

(1)步骤1：评价指标归一化。设样本数据x*(i,j)(i=1,2,…,n;j=1,2,…,p)为P维数据，其中n为样本总个数，p为评价指标总数目。归一化指标值按公式(1)计算。

(1)

式中：xmax(j)、xmin(j)分别为i分区j指标的最大值和最小值；x(i,j)为i分区j指标归一化后的值。

(2)步骤2：构造投影指标函数Q(a)。投影寻踪模型是将P维数据{x(i,j)|j=1,2,…,p}综合成以a={a(1),a(2),…,a(p)}为投影方向，并将归一化的指标简化为一维，从而可以构造投影指标函数。其中：

(2)

(i=1,2,…,n)

Q(a)=SZ·DZ

(3)

式中：a为投影方向，是一单位长度向量；Q(a)为投影指标函数，是Z(i)的函数，其值越大越好，所以投影寻踪就是寻求单位向量a(i)使Q(a)取最大值；SZ为类间距离，是投影特征值的标准差，SZ越大，散布越开；DZ为类内密度，是局部点的密集度，SZ及DZ的计算公式如下：

(4)

(5)

式中：E(z)为投影指标值Z(i)的平均值；R为局部密度的窗口半径，通常可以取值0.1SZ；r(i,j)为样本投影值之间的距离,即|z(i)-z(j)|；μ(t)函数为单位越阶函数，其中变量t为R-r(i,j),当t≥0时，其值取1，当t<0时，其值取0。

(3)步骤3：利用粒子群算法，优化指标函数。投影寻踪问题之关键在于寻找最大可能暴露高维数据某特征结构的最佳投影方向a，当评价指标、分级标准的样本集给定时，Q(a)只与投影方向a有关。因此，可以通过寻求Q(a)最大值，实现估计最佳投影方向。本文利用PSO算法求得最佳投影方向a*与Z*(i)。这里可以借助MATLAB编程实现，也可通过其他编程软件实现，其数学模型如下：

(6)

(4)步骤4：粒子群算法理论及步骤。设种群规模为N，xi为第i(i=1,2,…,N)个粒子的位置，vi为速度，fi(t)为t时刻的适应值。在随机产生初始位置和粒子速度之后的每一次迭代中，颗粒在t时刻通过跟踪个体极值pbesti(t)(粒子飞行过程中经历的最好位置)和全局极值gbest(t)(粒子群所经历过的最好位置)来更新自己在t+1时刻的位置，粒子群更新公式为：

vi(t+1)=Wvi(i)+c1r1(t)(pbest1(t)-xi(t))+

c2r2(t)(gbest(t)-xi(t))

(7)

xi(t+1)=xi(t)+vi(t+1)

(8)

式中：W为惯性权重；c1、c2为学习因子；r1、r2为(0，1)之间的均匀分布随机数。每个个体的局部极值和全部个体的全局极值更新公式如下：

pbesti(t+1)=

(9)

gbest(t+1)=xmax(t+1)

(10)

式中：xmax(t+1)为t+1时刻f(pbesti(t))取得最大值对应的粒子位置。其中的具体过程如下：

①规模为N的粒子群体，确定起始参数，如最大允许迭代次数为Gmax，计算各个粒子初始适应值fi(0)；

②根据公式(7)和(8)进行多次迭代，结合实际要求计算每个粒子新的速度、位置与适应值；

③逐个比较各粒子的最新适应值与前一时刻最优适应值的大小，若前者高，则以最新适应值替代前一时刻最优值；

④迭代过程中，若适应值误差或迭代次数超过相应的极限，则停止迭代，并输出结果。否则，继续返回第②步进行迭代。

(5)步骤5：将最佳投影方向代入公式(2)，得到指标及指标标准的一维投影值Z*(i)、Zr*(i)。将两者进行比较可得到灌区各分区的农业水土资源承载力评价等级。

3 结果分析与讨论

3.1 水土资源承载力综合分析结果

本文根据红寺堡灌区、固海固扩灌区的20个指标数据及分级指标中分级标准分界点按公式(1)进行归一化处理，采用Matlab R2017a及Python3.7.0编制程序计算，通过调整粒子群优化算法经验参数的数据，以尽量更好地找到函数最优值为目标，得到两个灌区各个分区的农业水土资源承载力评价结果。最佳投影方向a*及最佳投影值Z*(i)、Zr*(i)结果见表2～4。其中，为对结果进一步分析，再分别将4个子目标层B1、B2、B3、B4所对应的评价指标及分级标准进行相同步骤计算，从而得到红寺堡灌区、固海固扩灌区的4个系统层的最佳投影值。上述所得判定结果数据的粒子群优化算法参数设置如下：种群规模N=400；学习因子c1、c2分别取0.5和1.25；惯性权重W=0.9；最大迭代次数Gmax=200。

表2 指标分级标准分界点归一化及最佳投影值结果

表3 红寺堡灌区各分区指标归一化及水土资源评价结果

表4 固海、固扩区各分区指标归一化及水土资源评价结果

为了对结果进一步分析比较，分别求得各分区B1、B2、B3、B4层指标最佳投影方向a*和最佳投影值Z*(i)，再求出对应分级指标最佳投影值Zr*(i)，从而可以得到各分区子目标的水土资源承载力评价结果。将各评价值及等级判别结果分别绘制成坐标图，见图3。

表5 灌区各分区总目标层A水土资源承载力等级判别

图3 灌区各分区子目标层B1～B4水土资源承载力等级判别

3.2 讨论

由表5知，灌区各分区评价等级基本处于Ⅲ级和Ⅳ级，其中红寺堡灌区H2、H3分区和固海固扩灌区的G2、G6分区为Ⅲ级，水土资源承载力为中等。但更多分区都为Ⅳ级，其中红寺堡区以及固海固扩区的多个分区的承载力评价值小于2，且处于Ⅲ级的4个分区最佳投影指标值也较低，最高仅为2.809。可见，宁夏扬水灌区水土资源承载力总体水平较低。

总体上看，灌区水土资源承载力指标评价处于Ⅲ级的较少。由图3可看出，在红寺堡灌区中H2、H3虽然在A总目标层中都为Ⅲ级，但H2在B1、B2、B3系统层中都为Ⅳ级，只有在B4中为Ⅲ级；H3在B1中接近Ⅱ级，在B3中甚至只有Ⅴ级；可以看出B1、B4系统层往往对总体指标评价影响较大。在固海固扩灌区的G2、G6分区同样也是在B1、B4的系统层评价等级与总体指标评价相等或接近。因此从B1、B4系统层的指标分析，总指标评价等级为Ⅲ的分区其水资源利用的评价等级和生态环境的评价等级基本上处于中等，其分区都有一定量的供水保证，高耗水作物播种面积比(C2)较小，人均供水量(C7)较高，但单位耕地面积农业用水量(C5)较小，水资源需求量依然随着经济快速发展和人口数量增多而增加。在今后的水土资源开发利用中，应当充分保证调用扬水水资源的同时，还应提高农业灌溉效率，采用高效节水灌溉措施，以确保水资源的供给能力。

另外，从总体上看，灌区更多分区的总承载力评价结果较低，红寺堡灌区的H1、H4以及固海固扩管区的G1、G3、G4、G5分区处于Ⅳ级。从子系统层看，对总体指标评价影响较大主要是B1、B3、B4系统层的指标。这些分区在B1、B2指标中的评价相差较大，且针对B4系统层指标的评价等级主要为Ⅳ级和Ⅴ级。可见，这些分区的水土资源与社会经济发展不平衡，水资源的开发利用潜力较小，且其生态环境已经遭到极大的破坏。具体的，其B1系统层的指标高耗水作物面积比重(C2)大，农业用水量大；B3系统层的指标水利工程投资(C15)小；B4系统层的指标水土流失治理面积(C20)低，生态环境用水率(C19)也很低。而这些指标对这些分区的较低承载力等级有着主要影响，因此，这些分区应当从水资源利用率入手，减少水土流失，增加水资源供给；在生态环境上要留足生态用水量，保证环境自净能力；同时在作物结构上，要减少高耗水作物的种植面积，充分开发利用旱地，并采用高效节水灌溉措施等。从外部社会经济层面考虑，还应当加大对水资源调配的工程建设等水利工程的投资，以此综合改善灌区水土资源承载能力，避免情况恶化。

4 结 论

(1)本文通过建立粒子群优化投影模型对宁夏扬水灌区进行了多目标、多因素的水土资源承载力综合评价。从评价方法上看，评价结果与原始数据所直观表征的特点相符合，评价方法是科学而切实可行的。

(2)从评价结果上看，宁夏扬水灌区水土资源承载力总体水平较低，其中红寺堡区以及固海固扩区的多个分区的承载力评级等级只有Ⅳ级，且处于Ⅲ级的4个分区最佳投影指标值也较低(最高为2.809)，表明其地区水土资源与社会经济平衡发展的状态被破坏，系统不稳定。

(3)从4个子系统评价结果分析来看，各个分区的水土资源承载力的主要影响因素有差异，因此可以有针对性地对总目标层、子系统层的结果进行分析，从水资源利用、生态环境维护、种植作物类型以及水利工程投资等多角度对宁夏扬水灌区各分区实施综合改善措施。以此针对不同分区不同等级特点采取必要改善手段，以实现更为经济有效的扬水灌区水土资源可持续利用与管理。