基于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個子系統評價結果分析來看，各個分區的水土資源承載力的主要影響因素有差異，因此可以有針對性地對總目標層、子系統層的結果進行分析，從水資源利用、生態環境維護、種植作物類型以及水利工程投資等多角度對寧夏揚水灌區各分區實施綜合改善措施。以此針對不同分區不同等級特點采取必要改善手段，以實現更為經濟有效的揚水灌區水土資源可持續利用與管理。