基于PSO-PPE模型的寧夏揚(yáng)水灌區(qū)水土資源承載力綜合評價

鄒全程， 伍昊洋， 葉 威， 王春懿， 趙 璐，2， 崔寧博， 梁 川

(1.四川大學(xué) 水利水電學(xué)院， 四川 成都 610065； 2.農(nóng)業(yè)部旱作節(jié)水農(nóng)業(yè)重點開放實驗室， 北京 100081)

1 研究背景

土地資源是人類賴以生存和發(fā)展的基礎(chǔ)性自然資源，是農(nóng)業(yè)生產(chǎn)發(fā)展中的核心要素。通過對農(nóng)業(yè)水土資源匹配特征和承載力現(xiàn)狀的評價，可為土地規(guī)劃、管理及資源優(yōu)化配置提供理論依據(jù)，對區(qū)域的可持續(xù)發(fā)展具有重要作用[1]。

目前，國內(nèi)外有很多學(xué)者對水土資源承載力進(jìn)行了研究[2]，且有大量學(xué)者從不同角度，采用不同方法對水土資源承載力進(jìn)行了評價。拓學(xué)森等[3]綜合考慮了民勤縣水土資源與生態(tài)環(huán)境的實際特點，建立了民勤縣水土資源承載力系統(tǒng)動力學(xué)模型，用于預(yù)測未來水土資源承載力的變化趨勢；封志明等[4]綜合人糧、人水關(guān)系定量評價分析了雄安新區(qū)的水土資源承載力，提出了若干解決建議；張晶等[5]應(yīng)用GIS技術(shù)，研究了寧夏平原在天然和人工不同條件下的水土資源特征；李天霄等[6]采用實碼加速遺傳算法優(yōu)化投影指標(biāo)函數(shù)，建立水土資源承載力DPSIR模型指標(biāo)，分析了齊齊哈爾市市區(qū)的水土資源承載力等級情況。而投影尋蹤模型應(yīng)用到水土資源承載力評價的研究越來越多[7-9]，且隨著算法不斷改進(jìn)，粒子群優(yōu)化算法已經(jīng)逐漸被更多學(xué)者用于不同地區(qū)的水土資源承載力評價[10-11]，基于粒子群優(yōu)化算法的投影尋蹤模型已經(jīng)是目前應(yīng)用較成熟的評價方法之一。

自20個世紀(jì)70年代以來，為逐步解決寧夏回族自治區(qū)中部干旱帶日趨嚴(yán)重的干旱缺水狀況，該地區(qū)已經(jīng)興建了紅寺堡、固海、固海擴(kuò)灌和鹽環(huán)定等揚(yáng)水工程。為實現(xiàn)揚(yáng)水灌區(qū)水土資源科學(xué)的可持續(xù)利用和管理，對灌區(qū)進(jìn)行農(nóng)業(yè)水土資源匹配，承載力研究和農(nóng)業(yè)干旱風(fēng)險的評價十分重要?；诖?，本文構(gòu)建了農(nóng)業(yè)水土資源承載力綜合評價指標(biāo)體系，采用莊妍[8]研究的粒子群優(yōu)化投影尋蹤模型對寧夏揚(yáng)水灌區(qū)各個分區(qū)農(nóng)業(yè)水土資源承載力進(jìn)行了評價，以期為寧夏揚(yáng)水灌區(qū)水土資源可持續(xù)利用提供科學(xué)依據(jù)。

2 材料與方法

2.1 研究區(qū)域概況

寧夏揚(yáng)水灌區(qū)是自20世紀(jì)70年代發(fā)展的大型農(nóng)業(yè)灌溉區(qū)，該區(qū)域?qū)儆谥袦貛Ц珊祬^(qū)，太陽輻射強(qiáng)，大陸性氣候明顯，干旱少雨且強(qiáng)烈蒸發(fā)，致使區(qū)域土壤較干旱。

為研究灌區(qū)農(nóng)業(yè)水土資源的可持續(xù)利用情況，本文選擇比較完整的紅寺堡揚(yáng)水灌區(qū)(以下簡稱為紅寺堡灌區(qū))、固海揚(yáng)水灌區(qū)和固海擴(kuò)灌灌區(qū)(合并簡稱為固海固擴(kuò)灌區(qū))作為研究區(qū)域。為方便數(shù)據(jù)收集與研究，將灌區(qū)基本按所涉及縣域進(jìn)行分區(qū)，其中紅寺堡灌區(qū)分為H1、H2、H3、H4四個分區(qū)，所涉及縣域為利通區(qū)、紅寺堡區(qū)、同心縣、中寧縣；固海固擴(kuò)灌區(qū)分為G1、G2、G3、G4、G5、G6六個分區(qū)，所涉及縣域為紅寺堡區(qū)、同心縣、原州區(qū)、沙坡頭區(qū)、中寧縣、海原縣，圖1為研究區(qū)域分區(qū)概況圖。

圖1 揚(yáng)水灌區(qū)分區(qū)概況圖

2.2 評價指標(biāo)體系建立

水土資源存在相互制約的特點，本文在對寧夏揚(yáng)水灌區(qū)土地資源及各影響因素綜合分析的基礎(chǔ)上，參考國內(nèi)學(xué)者相關(guān)研究成果[13-16]，以灌區(qū)水土資源承載力評價為總目標(biāo)A，水土資源利用評價B1、社會發(fā)展評價B2、經(jīng)濟(jì)技術(shù)評價B3、生態(tài)環(huán)境評價B4為子目標(biāo)，構(gòu)建了包括20個指標(biāo)的分層評價體系，如圖2所示。

其中，指標(biāo)數(shù)據(jù)來源于《2015寧夏水利統(tǒng)計公報》、《2015寧夏水資源公報》、2016年《寧夏統(tǒng)計年鑒》、《寧夏大型揚(yáng)黃灌區(qū)用水現(xiàn)狀調(diào)查報告》等。根據(jù)指標(biāo)性質(zhì)，若指標(biāo)xi值越大，所反映的農(nóng)業(yè)水土資源承載力越高，稱為正向指標(biāo)，用(+)標(biāo)識；反之，則稱為負(fù)向指標(biāo)，用(-)標(biāo)識。

2.3 指標(biāo)分級標(biāo)準(zhǔn)確立

考慮紅寺堡灌區(qū)、固海固擴(kuò)灌區(qū)的特點，結(jié)合文獻(xiàn)[8,17-18]，將灌區(qū)農(nóng)業(yè)水土資源承載力評價結(jié)果劃分為5個等級，分別是高(Ⅰ)、較高 (Ⅱ)、中等 (Ⅲ)、 較低(Ⅳ)和低(Ⅴ)。Ⅰ～Ⅲ級代表區(qū)域生態(tài)功能完善，水土資源與社會經(jīng)濟(jì)發(fā)展系統(tǒng)很穩(wěn)定，處于可持續(xù)狀態(tài)；Ⅳ級代表生態(tài)功能相對完善，水土資源與社會經(jīng)濟(jì)平衡發(fā)展的狀態(tài)被破壞，系統(tǒng)不穩(wěn)定；Ⅴ級代表區(qū)域生態(tài)環(huán)境受到破壞，水土資源承載力已接近飽和。具體指標(biāo)見表1。

注：A、B1～B4、C1～C20分別為評價指標(biāo)體系內(nèi)的總目標(biāo)、4個子目標(biāo)、20個指標(biāo)。

指標(biāo)類型Ⅰ高Ⅱ較高Ⅲ中等Ⅳ較低Ⅴ低C1/%正(+)>55.141.4～55.127.6～41.413.8～27.60.0～13.8C2負(fù)(-)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)負(fù)(-)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)負(fù)(-)0.0～37.137.1～74.174.1～111.2111.2～148.0>148.3C10/‰負(fù)(-)0.0～4.54.5～8.98.9～13.413.4～17.8>17.8C11/104元負(fù)(-)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負(fù)(-)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)負(fù)(-)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)負(fù)(-)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：評價指標(biāo)歸一化。設(shè)樣本數(shù)據(jù)x*(i,j)(i=1,2,…,n;j=1,2,…,p)為P維數(shù)據(jù)，其中n為樣本總個數(shù)，p為評價指標(biāo)總數(shù)目。歸一化指標(biāo)值按公式(1)計算。

(1)

式中：xmax(j)、xmin(j)分別為i分區(qū)j指標(biāo)的最大值和最小值；x(i,j)為i分區(qū)j指標(biāo)歸一化后的值。

(2)步驟2：構(gòu)造投影指標(biāo)函數(shù)Q(a)。投影尋蹤模型是將P維數(shù)據(jù){x(i,j)|j=1,2,…,p}綜合成以a={a(1),a(2),…,a(p)}為投影方向，并將歸一化的指標(biāo)簡化為一維，從而可以構(gòu)造投影指標(biāo)函數(shù)。其中：

(2)

(i=1,2,…,n)

Q(a)=SZ·DZ

(3)

式中：a為投影方向，是一單位長度向量；Q(a)為投影指標(biāo)函數(shù)，是Z(i)的函數(shù)，其值越大越好，所以投影尋蹤就是尋求單位向量a(i)使Q(a)取最大值；SZ為類間距離，是投影特征值的標(biāo)準(zhǔn)差，SZ越大，散布越開；DZ為類內(nèi)密度，是局部點的密集度，SZ及DZ的計算公式如下：

(4)

(5)

式中：E(z)為投影指標(biāo)值Z(i)的平均值；R為局部密度的窗口半徑，通?？梢匀≈?.1SZ；r(i,j)為樣本投影值之間的距離,即|z(i)-z(j)|；μ(t)函數(shù)為單位越階函數(shù)，其中變量t為R-r(i,j),當(dāng)t≥0時，其值取1，當(dāng)t<0時，其值取0。

(3)步驟3：利用粒子群算法，優(yōu)化指標(biāo)函數(shù)。投影尋蹤問題之關(guān)鍵在于尋找最大可能暴露高維數(shù)據(jù)某特征結(jié)構(gòu)的最佳投影方向a，當(dāng)評價指標(biāo)、分級標(biāo)準(zhǔn)的樣本集給定時，Q(a)只與投影方向a有關(guān)。因此，可以通過尋求Q(a)最大值，實現(xiàn)估計最佳投影方向。本文利用PSO算法求得最佳投影方向a*與Z*(i)。這里可以借助MATLAB編程實現(xiàn)，也可通過其他編程軟件實現(xiàn)，其數(shù)學(xué)模型如下：

(6)

(4)步驟4：粒子群算法理論及步驟。設(shè)種群規(guī)模為N，xi為第i(i=1,2,…,N)個粒子的位置，vi為速度，fi(t)為t時刻的適應(yīng)值。在隨機(jī)產(chǎn)生初始位置和粒子速度之后的每一次迭代中，顆粒在t時刻通過跟蹤個體極值pbesti(t)(粒子飛行過程中經(jīng)歷的最好位置)和全局極值gbest(t)(粒子群所經(jīng)歷過的最好位置)來更新自己在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為慣性權(quán)重；c1、c2為學(xué)習(xí)因子；r1、r2為(0，1)之間的均勻分布隨機(jī)數(shù)。每個個體的局部極值和全部個體的全局極值更新公式如下：

pbesti(t+1)=

(9)

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

(10)

式中：xmax(t+1)為t+1時刻f(pbesti(t))取得最大值對應(yīng)的粒子位置。其中的具體過程如下：

①規(guī)模為N的粒子群體，確定起始參數(shù)，如最大允許迭代次數(shù)為Gmax，計算各個粒子初始適應(yīng)值fi(0)；

②根據(jù)公式(7)和(8)進(jìn)行多次迭代，結(jié)合實際要求計算每個粒子新的速度、位置與適應(yīng)值；

③逐個比較各粒子的最新適應(yīng)值與前一時刻最優(yōu)適應(yīng)值的大小，若前者高，則以最新適應(yīng)值替代前一時刻最優(yōu)值；

④迭代過程中，若適應(yīng)值誤差或迭代次數(shù)超過相應(yīng)的極限，則停止迭代，并輸出結(jié)果。否則，繼續(xù)返回第②步進(jìn)行迭代。

(5)步驟5：將最佳投影方向代入公式(2)，得到指標(biāo)及指標(biāo)標(biāo)準(zhǔn)的一維投影值Z*(i)、Zr*(i)。將兩者進(jìn)行比較可得到灌區(qū)各分區(qū)的農(nóng)業(yè)水土資源承載力評價等級。

3 結(jié)果分析與討論

3.1 水土資源承載力綜合分析結(jié)果

本文根據(jù)紅寺堡灌區(qū)、固海固擴(kuò)灌區(qū)的20個指標(biāo)數(shù)據(jù)及分級指標(biāo)中分級標(biāo)準(zhǔn)分界點按公式(1)進(jìn)行歸一化處理，采用Matlab R2017a及Python3.7.0編制程序計算，通過調(diào)整粒子群優(yōu)化算法經(jīng)驗參數(shù)的數(shù)據(jù)，以盡量更好地找到函數(shù)最優(yōu)值為目標(biāo)，得到兩個灌區(qū)各個分區(qū)的農(nóng)業(yè)水土資源承載力評價結(jié)果。最佳投影方向a*及最佳投影值Z*(i)、Zr*(i)結(jié)果見表2～4。其中，為對結(jié)果進(jìn)一步分析，再分別將4個子目標(biāo)層B1、B2、B3、B4所對應(yīng)的評價指標(biāo)及分級標(biāo)準(zhǔn)進(jìn)行相同步驟計算，從而得到紅寺堡灌區(qū)、固海固擴(kuò)灌區(qū)的4個系統(tǒng)層的最佳投影值。上述所得判定結(jié)果數(shù)據(jù)的粒子群優(yōu)化算法參數(shù)設(shè)置如下：種群規(guī)模N=400；學(xué)習(xí)因子c1、c2分別取0.5和1.25；慣性權(quán)重W=0.9；最大迭代次數(shù)Gmax=200。

表2 指標(biāo)分級標(biāo)準(zhǔn)分界點歸一化及最佳投影值結(jié)果

表3 紅寺堡灌區(qū)各分區(qū)指標(biāo)歸一化及水土資源評價結(jié)果

表4 固海、固擴(kuò)區(qū)各分區(qū)指標(biāo)歸一化及水土資源評價結(jié)果

為了對結(jié)果進(jìn)一步分析比較，分別求得各分區(qū)B1、B2、B3、B4層指標(biāo)最佳投影方向a*和最佳投影值Z*(i)，再求出對應(yīng)分級指標(biāo)最佳投影值Zr*(i)，從而可以得到各分區(qū)子目標(biāo)的水土資源承載力評價結(jié)果。將各評價值及等級判別結(jié)果分別繪制成坐標(biāo)圖，見圖3。

表5 灌區(qū)各分區(qū)總目標(biāo)層A水土資源承載力等級判別

圖3 灌區(qū)各分區(qū)子目標(biāo)層B1～B4水土資源承載力等級判別

3.2 討論

由表5知，灌區(qū)各分區(qū)評價等級基本處于Ⅲ級和Ⅳ級，其中紅寺堡灌區(qū)H2、H3分區(qū)和固海固擴(kuò)灌區(qū)的G2、G6分區(qū)為Ⅲ級，水土資源承載力為中等。但更多分區(qū)都為Ⅳ級，其中紅寺堡區(qū)以及固海固擴(kuò)區(qū)的多個分區(qū)的承載力評價值小于2，且處于Ⅲ級的4個分區(qū)最佳投影指標(biāo)值也較低，最高僅為2.809?？梢姡瑢幭膿P(yáng)水灌區(qū)水土資源承載力總體水平較低。

總體上看，灌區(qū)水土資源承載力指標(biāo)評價處于Ⅲ級的較少。由圖3可看出，在紅寺堡灌區(qū)中H2、H3雖然在A總目標(biāo)層中都為Ⅲ級，但H2在B1、B2、B3系統(tǒng)層中都為Ⅳ級，只有在B4中為Ⅲ級；H3在B1中接近Ⅱ級，在B3中甚至只有Ⅴ級；可以看出B1、B4系統(tǒng)層往往對總體指標(biāo)評價影響較大。在固海固擴(kuò)灌區(qū)的G2、G6分區(qū)同樣也是在B1、B4的系統(tǒng)層評價等級與總體指標(biāo)評價相等或接近。因此從B1、B4系統(tǒng)層的指標(biāo)分析，總指標(biāo)評價等級為Ⅲ的分區(qū)其水資源利用的評價等級和生態(tài)環(huán)境的評價等級基本上處于中等，其分區(qū)都有一定量的供水保證，高耗水作物播種面積比(C2)較小，人均供水量(C7)較高，但單位耕地面積農(nóng)業(yè)用水量(C5)較小，水資源需求量依然隨著經(jīng)濟(jì)快速發(fā)展和人口數(shù)量增多而增加。在今后的水土資源開發(fā)利用中，應(yīng)當(dāng)充分保證調(diào)用揚(yáng)水水資源的同時，還應(yīng)提高農(nóng)業(yè)灌溉效率，采用高效節(jié)水灌溉措施，以確保水資源的供給能力。

另外，從總體上看，灌區(qū)更多分區(qū)的總承載力評價結(jié)果較低，紅寺堡灌區(qū)的H1、H4以及固海固擴(kuò)管區(qū)的G1、G3、G4、G5分區(qū)處于Ⅳ級。從子系統(tǒng)層看，對總體指標(biāo)評價影響較大主要是B1、B3、B4系統(tǒng)層的指標(biāo)。這些分區(qū)在B1、B2指標(biāo)中的評價相差較大，且針對B4系統(tǒng)層指標(biāo)的評價等級主要為Ⅳ級和Ⅴ級?？梢姡@些分區(qū)的水土資源與社會經(jīng)濟(jì)發(fā)展不平衡，水資源的開發(fā)利用潛力較小，且其生態(tài)環(huán)境已經(jīng)遭到極大的破壞。具體的，其B1系統(tǒng)層的指標(biāo)高耗水作物面積比重(C2)大，農(nóng)業(yè)用水量大；B3系統(tǒng)層的指標(biāo)水利工程投資(C15)?。籅4系統(tǒng)層的指標(biāo)水土流失治理面積(C20)低，生態(tài)環(huán)境用水率(C19)也很低。而這些指標(biāo)對這些分區(qū)的較低承載力等級有著主要影響，因此，這些分區(qū)應(yīng)當(dāng)從水資源利用率入手，減少水土流失，增加水資源供給；在生態(tài)環(huán)境上要留足生態(tài)用水量，保證環(huán)境自凈能力；同時在作物結(jié)構(gòu)上，要減少高耗水作物的種植面積，充分開發(fā)利用旱地，并采用高效節(jié)水灌溉措施等。從外部社會經(jīng)濟(jì)層面考慮，還應(yīng)當(dāng)加大對水資源調(diào)配的工程建設(shè)等水利工程的投資，以此綜合改善灌區(qū)水土資源承載能力，避免情況惡化。

4 結(jié) 論

(1)本文通過建立粒子群優(yōu)化投影模型對寧夏揚(yáng)水灌區(qū)進(jìn)行了多目標(biāo)、多因素的水土資源承載力綜合評價。從評價方法上看，評價結(jié)果與原始數(shù)據(jù)所直觀表征的特點相符合，評價方法是科學(xué)而切實可行的。

(2)從評價結(jié)果上看，寧夏揚(yáng)水灌區(qū)水土資源承載力總體水平較低，其中紅寺堡區(qū)以及固海固擴(kuò)區(qū)的多個分區(qū)的承載力評級等級只有Ⅳ級，且處于Ⅲ級的4個分區(qū)最佳投影指標(biāo)值也較低(最高為2.809)，表明其地區(qū)水土資源與社會經(jīng)濟(jì)平衡發(fā)展的狀態(tài)被破壞，系統(tǒng)不穩(wěn)定。

(3)從4個子系統(tǒng)評價結(jié)果分析來看，各個分區(qū)的水土資源承載力的主要影響因素有差異，因此可以有針對性地對總目標(biāo)層、子系統(tǒng)層的結(jié)果進(jìn)行分析，從水資源利用、生態(tài)環(huán)境維護(hù)、種植作物類型以及水利工程投資等多角度對寧夏揚(yáng)水灌區(qū)各分區(qū)實施綜合改善措施。以此針對不同分區(qū)不同等級特點采取必要改善手段，以實現(xiàn)更為經(jīng)濟(jì)有效的揚(yáng)水灌區(qū)水土資源可持續(xù)利用與管理。