李 茉，郭 萍，郭珊珊

(中國農(nóng)業(yè)大學(xué)中國農(nóng)業(yè)水問題研究中心，北京 100083)

水資源的合理利用在農(nóng)業(yè)生產(chǎn)中發(fā)揮關(guān)鍵的作用[1]。農(nóng)業(yè)是我國用水大戶，用水比例達(dá)63.5%[2]。社會經(jīng)濟的快速發(fā)展使工業(yè)、生活和生態(tài)用水與農(nóng)業(yè)用水競爭激烈，導(dǎo)致灌溉用水供需矛盾日益突出，且我國灌溉用水管理不當(dāng)，效率低。如何在考慮區(qū)域社會經(jīng)濟和生態(tài)可持續(xù)發(fā)展的前提下合理且高效的配置灌溉水資源，成為區(qū)域水資源決策者面臨的緊迫任務(wù)。

國內(nèi)外關(guān)于灌溉水資源優(yōu)化配置已有大量研究，其中普通線性規(guī)劃模型由于計算簡便被廣泛應(yīng)用[3]。針對我國灌溉用水效率低的特點，優(yōu)化配置有限的灌溉水資源以提高區(qū)域灌溉水生產(chǎn)率比單純提高產(chǎn)量或產(chǎn)值更有實際意義，特別是在干旱和半干旱地區(qū)。因此，基于分式線性規(guī)劃的灌溉水資源優(yōu)化配置可被采用[4,5]?？紤]可持續(xù)發(fā)展的灌溉水資源優(yōu)化配置需建立在灌溉用水安全閾值的基礎(chǔ)上。水資源安全閾值的研究主要是確定區(qū)域水資源系統(tǒng)安全區(qū)間的上下限值[6]，在此區(qū)間內(nèi)，水資源處于安全狀態(tài)，社會經(jīng)濟與生態(tài)可持續(xù)發(fā)展。對于水資源緊缺的干旱半干旱地區(qū)，尋找水資源安全閾值的下限值即不發(fā)生干旱的最小用水指標(biāo)顯得尤為重要[7]。近些年，國內(nèi)關(guān)于水資源安全閾值的研究成為熱點[7-9]，但多數(shù)集中于區(qū)域尺度上的研究，對側(cè)重考慮灌溉用水安全閾值的研究較少，尤其是基于灌溉用水安全閾值的灌溉水資源高效配置的研究鮮有報道。

本文以黑河中游的甘州區(qū)、臨澤縣、高臺縣為研究區(qū)域，首先采用模糊識別模型計算區(qū)域灌溉水資源與社會經(jīng)濟生態(tài)協(xié)調(diào)發(fā)展的協(xié)調(diào)度，并以此協(xié)調(diào)度作為狀態(tài)變量，與由因子分析法確定的控制變量一起擬合基于突變理論的尖點突變模型，并以此確定灌溉水資源安全閾值。在此基礎(chǔ)上，構(gòu)建灌溉用水高效配置優(yōu)化模型并求解，得到的優(yōu)化配水方案對研究區(qū)域水資源管理具有實際指導(dǎo)意義。

1 可變模糊識別模型與熵值法

可變模糊識別模型理論與方法是建立在相對隸屬度與相對隸屬函數(shù)概念的基礎(chǔ)上的。該方法可以反映水資源系統(tǒng)的復(fù)雜性、不確定性、一定時段內(nèi)具有的動態(tài)性等特點，并且能夠客觀分析水資源與社會經(jīng)濟生態(tài)協(xié)調(diào)發(fā)展的狀況。可變模糊識別模型可表示為[10]：

(1)

式中：vj為方案集綜合相對優(yōu)屬度；i=1,2,…,I為指標(biāo)種類；J=1,2,…,J為時段數(shù)；uij為指標(biāo)i時段j的因子值；ωi為指標(biāo)i的權(quán)重；μij(uij)為指標(biāo)i時段j的相對隸屬度；α為優(yōu)化準(zhǔn)則參數(shù)，常取α=1(最小一乘方準(zhǔn)則)和α=2最小二乘方準(zhǔn)則)；p為距離參數(shù)，常取p=1(海明距離)和p=2(歐氏距離)。

令dbj={∑Ii=1[ωi(1-μij(uij)]p}1/p表示時段j的最優(yōu)距離，dwj={∑Ii=1[ωiμij(uij)]p}1/p表示時段j的最劣距離，公式(1)中不同α和p取值的隨機組合可形成以下4種模型：①當(dāng)α=1,p=1時，vj=∑Ii=1ωiμij(uij)，相當(dāng)于模糊綜合評判模型；②當(dāng)α=1,p=2時，vj=dwj/(dbj+dwj)，相當(dāng)于理想點模型；③當(dāng)α=2,p=1時，vj=1/{1+[(1-dwj)/dwj]2}，此時為sigmoid函數(shù)，是一個良好的閾值函數(shù)；④當(dāng)α=2,p=2時，vj=1/[1+(dbj/dwj)2]，此時為模糊優(yōu)選模型。

在上述各模糊識別模型中，一個重要的參數(shù)就是指標(biāo)的權(quán)重ωi。關(guān)于指標(biāo)權(quán)重的確定方法有很多，其中熵值法確定權(quán)重由于其客觀合理性在水資源評價領(lǐng)域中被廣泛應(yīng)用[11]，其確定過程可總結(jié)如下：①歸一化非負(fù)處理評價因子集u′ij=[(uij-min(uij))/(max(uij)-min(uij))]+1，其中“+1”項為避免求熵時取對數(shù)的無意義；②計算指標(biāo)i時段j占總體比例pij=u′ij/∑Jju′ij；③計算各因子權(quán)重ωi=di/∑Ii=1di，其中di=(1-ei)/(1-∑Ii=1ei)，ei=-[1/(ln(J)∑Jj=1pijln(pij))]。

選取代表灌溉水資源、經(jīng)濟、社會及生態(tài)各維度的典型指標(biāo)集合，記為B，在此基礎(chǔ)上，根據(jù)上述各可變模糊識別模型可以計算研究區(qū)域不同年份的灌溉水資源與社會經(jīng)濟生態(tài)協(xié)調(diào)發(fā)展的協(xié)調(diào)度(以下簡稱“協(xié)調(diào)度”)。由于該“協(xié)調(diào)度”代表灌溉水資源的可持續(xù)利用狀態(tài)，在一定程度上反映灌溉水資源的安全程度，因此本文選用此“協(xié)調(diào)度”作為后述尖點突變模型的狀態(tài)變量。

2 尖點突變模型

突變理論是根據(jù)勢函數(shù)來研究對象的變化過程和突變線性，系統(tǒng)的勢函數(shù)可表示系統(tǒng)任一狀態(tài)值，這個值是控制變量與狀態(tài)變量的統(tǒng)一。水資源系統(tǒng)滿足尖點突變的兩個特性，突跳性和發(fā)散性，因此可采用尖點突變模型對水資源安全系統(tǒng)進行突變分析[7]。尖點突變模型由兩個控制變量和一個狀態(tài)變量構(gòu)成，其勢函數(shù)可表示成：

V(x)=x4+qx2+rx

(2)

式中：x為狀態(tài)變量；q為主控制變量；r為次控制變量。

對勢函數(shù)求導(dǎo)可得到其臨界點，其表達(dá)式如下：

(3)

式(3)為尖點突變模型的平衡曲面，聯(lián)立式(2)和式(3)，得到分歧點方程為

8q3+27r2=0

(4)

令Δ=8q3+27r2代表突變判別式，當(dāng)Δ>0時，系統(tǒng)是穩(wěn)定安全的，當(dāng)Δ<0時，系統(tǒng)會發(fā)生突變，當(dāng) 時是發(fā)生突變的臨界值。

實際應(yīng)用中的具體步驟可表述如下：①選取勢函數(shù)中的控制變量和狀態(tài)變量。其中狀態(tài)變量為根據(jù)模糊識別模型求得的“協(xié)調(diào)度”，對于主控制變量和次控制變量的選擇，若指標(biāo)集合B中指標(biāo)較多，可先采用因子分析方法對指標(biāo)進行篩選，根據(jù)旋轉(zhuǎn)后的因子負(fù)荷矩陣[12]排序選取主、次控制變量；②變量(包括控制變量和狀態(tài)變量)歸一化處理，其中“越大越優(yōu)”指標(biāo)采用公式u″ij=(u″ij-uijmin)/(uijmax-uijmin)計算，而“越小越優(yōu)”指標(biāo)采用公式u″ij=1-(u″ij-uijmin)/(uijmax-uijmin)計算；③熵值法計算各控制變量權(quán)重；④模型擬合，將尖點突變模型的平衡曲面寫成4x3=-2qx-r的形式，令y=4x3，則尖點突變平衡曲面擬合式可表述為y=k1(-2q′x′)+k2(-r′)+k3，其中q′，r′和x′分別為主控制變量、次控制變量和狀態(tài)變量的歸一化值，采用多元回歸分析法求得k1，k2，k3值；④根據(jù)擬合結(jié)果計算判別式Δ從而獲得灌溉用水安全閾值。

上述步驟中涉及用因子分析法對初始指標(biāo)體系進行篩選得到主、次控制變量，其中因子分析法是從研究變量內(nèi)部相關(guān)的依賴關(guān)系出發(fā)，把一些具有錯綜復(fù)雜關(guān)系的變量歸納為少數(shù)綜合因子的一種多變量統(tǒng)計方法[13]，計算步驟可歸納如下：①將原始指標(biāo)體系標(biāo)準(zhǔn)化；②求標(biāo)準(zhǔn)化矩陣的相關(guān)矩陣；③求相關(guān)矩陣的特征值和特征向量；④計算各因子方差貢獻率和累計方差貢獻率；⑤根據(jù)方差累計貢獻率程度(不低于85%)確定公共因子；⑥通過正交旋轉(zhuǎn)法進行因子旋轉(zhuǎn)；⑦根據(jù)旋轉(zhuǎn)后因子負(fù)荷矩陣篩選主要影響因子。

3 灌溉水資源優(yōu)化配置模型

以區(qū)域灌溉水生產(chǎn)率最大為目標(biāo)函數(shù)，在有限水資源和灌溉用水安全閾值等約束下，對研究區(qū)域的地表水和地下水進行綜合配置。配水模型為線性分式規(guī)劃模型，分子為灌溉所獲得的產(chǎn)量，分子為區(qū)域分配的水量，分子分母同時優(yōu)化，以期獲得區(qū)域最小灌水量下的最大產(chǎn)量。其數(shù)學(xué)表達(dá)式如下：

(5)

地表水可利用量約束：

(6)

地下水可利用量約束：

GWi≤Gi?i

(7)

灌溉用水安全閾值約束：

SWi+GWi≥Mi?i

(8)

非負(fù)約束：

SWi≥0,GWi≥0 ?i

(9)

式中：Z為目標(biāo)函數(shù)，kg/m3；i為研究區(qū)域；SWi、GWi分別為區(qū)域i的地表配水量和地下配水量(決策變量)，m3；Yi為區(qū)域i的單方水產(chǎn)量，kg/m3；p為灌溉用水比例；以由各區(qū)域i組成的整體為閉合系統(tǒng)，Qu，Qs，Qd分別為河流上游來水，自產(chǎn)水與承諾的下泄水量，m3；Gi為區(qū)域i可開采的被農(nóng)業(yè)利用的地下水量，m3；Mi為灌溉用水安全閾值下限，m3。

4 實例應(yīng)用與結(jié)果分析

研究區(qū)域選定為黑河中游的甘州區(qū)、臨澤縣和高臺縣。黑河流域中游段即從鶯落峽至正義峽中間段，多年平均降雨125 mm，年均蒸發(fā)量1 200 mm，蒸發(fā)量遠(yuǎn)大于降雨量，水資源短缺嚴(yán)重。中游農(nóng)業(yè)用水占總用水量比例高達(dá)90%，灌溉水資源主要引黑河干流水和抽取地下水，還有少量小河流產(chǎn)水。根據(jù)國務(wù)院關(guān)于黑河分水目標(biāo)，當(dāng)鶯落峽多年平均來水為15.8億m3時，正義峽斷面的泄水量需達(dá)到9.5億m3，以保證黑河下游生態(tài)健康。中游灌溉可用水量因此被大量縮減，如何在保證中游灌溉用水安全的前提下高效的配置有限水資源，使得黑河中游段的水資源與社會經(jīng)濟生態(tài)協(xié)調(diào)發(fā)展是水資源管理者需要考慮的問題。

根據(jù)多年統(tǒng)計與調(diào)研資料，獲得能夠代表三區(qū)(縣)灌溉水資源、經(jīng)濟、社會和生態(tài)相關(guān)的典型指標(biāo)體系如表1所示(以甘州區(qū)為例，各研究區(qū)域指標(biāo)一致)。首先根據(jù)確定的指標(biāo)體系計算三區(qū)(縣)的“協(xié)調(diào)度”，本文將四種可變模糊識別模型均做計算，最終取四種模型的平均值作為各研究區(qū)域的“協(xié)調(diào)度”，如圖1所示，其中由熵值法確定的三區(qū)(縣)各指標(biāo)權(quán)重見表2，“協(xié)調(diào)度”的判斷標(biāo)準(zhǔn)為見表3。結(jié)合圖1和表3得到甘州區(qū)的水資源與社會經(jīng)濟生態(tài)發(fā)展的協(xié)調(diào)程度在2000-2013年間呈好轉(zhuǎn)趨勢，其“協(xié)調(diào)度”從“明顯失調(diào)-動態(tài)平衡”狀態(tài)逐漸轉(zhuǎn)變成“動態(tài)平衡-基本協(xié)調(diào)”狀態(tài)，水資源基本呈可持續(xù)利用狀態(tài)。而臨澤縣和高臺縣的“協(xié)調(diào)度”整體呈下降趨勢，其中臨澤縣2000-2009年間的“協(xié)調(diào)度”呈波動狀態(tài)，2009-2012年間的“協(xié)調(diào)度”大幅度下降，2013年有所好轉(zhuǎn)，高臺縣2000-2002年間的“協(xié)調(diào)度”呈上升趨勢，之后保持三年平穩(wěn)狀態(tài)，然后呈逐年下降趨勢。整體來講，臨澤縣和高臺縣水資源呈不可持續(xù)利用狀態(tài)，其原因可能是由于“分水”計劃導(dǎo)致中游水資源短缺嚴(yán)重，而甘州區(qū)作為張掖市的政治、經(jīng)濟、文化中心，在有限的水資源情況下，優(yōu)先保證甘州區(qū)的水資源可持續(xù)利用。將計算得到的三區(qū)(縣)4種模型“協(xié)調(diào)度”的歸一化平均值作為尖點突變模型的狀態(tài)變量。

表1 甘州區(qū)指標(biāo)體系Tab.1 Index system of Ganzhou district

圖1 甘州區(qū)、臨澤縣和高臺縣的“協(xié)調(diào)度”Fig.1 Coordination degree of Ganzhou district, Linze County and Gaotai County

區(qū)域水資源指標(biāo)用水總量/億m3灌溉定額/(m3·hm-2)灌溉水利用系數(shù)/%地下水量/億m3萬元GDP用水量/(m3·萬元-1)社會經(jīng)濟指標(biāo)GDP總量/億元人口/(m3·萬元-1)人均GDP/(萬元·人-1)糧食總產(chǎn)/億kg有效灌溉面積/萬hm2城鎮(zhèn)人均收入/(萬元·人-1)農(nóng)民人均收入/(萬元·人-1)水費收入/萬元單方水產(chǎn)量/(kg·m-3)糧食單位面積產(chǎn)量/(kg·hm-2)生態(tài)指標(biāo)年降雨量/mm林草灌溉面積/萬hm2甘州區(qū)0．060．050．050．040．07950．060．040．060．060．110．0550．0630．0550．070．070．040．03臨澤縣0．040．040．040．050．06670．070．070．050．060．140．0650．0690．0560．0520．070．040．03高臺縣0．050．040．070．050．04610．050．060．060．070．080．0530．0550．0710．0340．050．050．12

表3 “協(xié)調(diào)度”的判斷標(biāo)準(zhǔn)Tab.3 Criteria of coordination degree

各區(qū)(縣)指標(biāo)體系均包含17個指標(biāo)，指標(biāo)相對較多，本文通過因子分析方法對各區(qū)(縣)指標(biāo)進行篩選，根據(jù)旋轉(zhuǎn)后的因子負(fù)荷矩陣篩選主、次控制變量。表4為三區(qū)(縣)旋轉(zhuǎn)后因子負(fù)荷矩陣，根據(jù)各區(qū)(縣)各因子所占負(fù)荷值并協(xié)調(diào)各綜合因子，同時考慮與灌溉用水密切相關(guān)的指標(biāo)，最終確定甘州區(qū)的主控制指標(biāo)為用水總量、年降雨量、灌溉定額，次控制指標(biāo)為農(nóng)民人均純收入、糧食單方水產(chǎn)量、水費收入；臨澤縣的主控制指標(biāo)為用水總量、灌溉水利用系數(shù)，次控制指標(biāo)為萬元GDP用水量、有效灌溉面積、林草灌溉面積；高臺縣的主控制指標(biāo)為用水總量、灌溉水利用系數(shù)、年降雨量，次控制指標(biāo)為萬元GDP用水量、GDP總量、單方水產(chǎn)量。對選好的控制變量進行歸一化處理，并根據(jù)熵值法對選好的指標(biāo)進行權(quán)重計算，擬合尖點突變模型，并擬合判別式與灌溉毛用水量之間的函數(shù)關(guān)系，如圖2所示。由于黑河中游處于水資源短缺地區(qū)，因此尋找的灌溉用水安全閾值為對應(yīng)的下限值，即不發(fā)生干旱災(zāi)害的最小灌溉用水量。從圖2中可以看出，甘州區(qū)和高臺縣歷年灌溉用水情況安全，只有個別年份判別式小于零，而臨澤縣的灌溉用水因其歷年判別式均小于零呈現(xiàn)不安全狀態(tài)，臨澤縣的灌溉用水在很大程度上不能滿足需求，極易發(fā)生干旱。根據(jù)圖2擬合的結(jié)果，可以得到甘州區(qū)、臨澤縣和高臺縣判別式與毛灌溉水量之間的函數(shù)關(guān)系分別為Δ甘=49x2-754.76x+2 901.4，Δ臨=-6.29x2+77.93x-247.35，Δ高=276.69x2-2 563.2x+5 941.4。通過對判別式求極值得到甘州區(qū)、臨澤縣和高臺縣的毛用水安全閾值分別為7.70，6.20和4.63億m3。根據(jù)三區(qū)(縣)多年渠系水利用系數(shù)平均值(甘州區(qū)：68.79%，臨澤縣：63.23%，高臺縣60.80%)和三區(qū)(縣)灌溉用水比例，得到甘州區(qū)、臨澤縣和高臺縣的凈灌溉用水安全閾值(Mi)分別為4.88，3.02和2.69億m3。

將得到的凈灌溉用水安全閾值輸入到高效配水優(yōu)化模型中并對模型進行求解得到三區(qū)(縣)地表水和地下水優(yōu)化配水方案。地表總可利用水量為：鶯落峽徑流(Qu=15.8億m3)+境內(nèi)河流產(chǎn)流(梨園河2.37億m3+其他小河流1億m3，即Qs= 3.37億m3)-下泄到正義峽的水量(Qd= 9.5億m3)。甘州區(qū)、臨澤縣和高臺縣的地下水可開采量(Gi)分別為2.01，1.3和1.5億m3。甘州區(qū)、臨澤縣和高臺縣的單方水產(chǎn)量(Yi)分別為1.60，1.79和1.23 kg/m3。模型優(yōu)化結(jié)果見圖3。優(yōu)化的灌溉水生產(chǎn)率為1.58 kg/m3，比多年平均情況增長1.3%；優(yōu)化的總灌水量為10.28 m3，比多年平均情況降低3%。優(yōu)化的結(jié)果在保證三區(qū)(縣)灌溉用水安全的情況下，節(jié)約了水資源，并使灌溉水生產(chǎn)率有小幅度提高。從優(yōu)化結(jié)果中可以得到，分配給三區(qū)(縣)的總地下水達(dá)到4.7億m3，雖然在總地下水可開采能力(4.81 m3)范圍內(nèi)，但地下水用量比多年平均多了1.32億m3。由于配水是在滿足下泄到正義峽水量達(dá)到9.5億m3和灌溉用水安全閾值的基礎(chǔ)上的，此時地表徑流全部利用，若使地下用水量維持多年平均水平，則會導(dǎo)致下泄到正義峽的水量相應(yīng)減少或不能完全滿足灌溉用水安全閾值要求，決策者可根據(jù)實際情況與需求對配水方案進行相應(yīng)調(diào)整。

5 結(jié) 語

(1)利用模糊可變識別模型計算甘州區(qū)、臨澤縣和高臺縣水資源與社會經(jīng)濟生態(tài)協(xié)調(diào)發(fā)展的協(xié)調(diào)度，其中甘州區(qū)的“協(xié)調(diào)度”呈上升趨勢，臨澤縣和高臺縣的“協(xié)調(diào)度”呈下降趨勢。

表4 因子分析旋轉(zhuǎn)后因子負(fù)荷矩陣Tab.4 The rotating factor loading matrix of factor analysis method

圖2 甘州區(qū)、臨澤縣和高臺縣的判別式值與用水量擬合Fig.2 Fitting results of the discriminant and irrigation water use of the three study areas

圖3 三區(qū)(縣)的優(yōu)化配水結(jié)果Fig.3 Optimal water allocation of the three study areas

(2)以四種模糊可變識別模型計算的“協(xié)調(diào)度”均值作為狀態(tài)變量，結(jié)合利用因子分析法篩選出的主、次控制變量擬合尖點突變模型，根據(jù)判別式確定甘州區(qū)、臨澤縣和高臺縣的灌溉用水安全閾值下限分別為4.88，3.02和2.69億m3。

(3)以灌溉用水安全閾值為約束條件，以最大化灌溉水生產(chǎn)率為目標(biāo)函數(shù)，構(gòu)建灌溉水資源優(yōu)化配置模型，優(yōu)化得到灌溉水生產(chǎn)率為1.58 kg/m3，比多年平均增長1.3%；優(yōu)化的總灌水量為10.28 m3，比多年平均降低3%。

□

[1] 張展羽, 司 涵, 馮寶平, 等. 缺水灌區(qū)農(nóng)業(yè)水土資源優(yōu)化配置模型[J]. 水利學(xué)報, 2014,45(4):403-409.

[2] 中華人民共和國水利部. 中國水資源公報[R]. 2014.

[3] Singh A. Irrigation planning and management through optimization modelling [J]. Water Resources Management, 2014,28(1):1-14.

[4] Li M, Guo P, Singh V P, et al. An efficient irrigation water allocation model under uncertainty [J]. Agricultural Systems, 2016, 144:46-57.

[5] Guo P, Chen X H, Li M, et al. Fuzzy chance-constrained linear fractional programming approach for optimal water allocation [J]. Stochastic Environmental Research and Risk Assessment, 2014,28:1 601-1 612.

[6] 王靄景. 天津市水資源優(yōu)化配置及安全閾值研究[D]. 北京：華北電力大學(xué)，2013.

[7] 張玉山, 李繼清, 梅艷艷, 等. 基于突變理論的天津市水資源安全閾值分析模型[J]. 遼寧工程技術(shù)大學(xué)學(xué)報(自然科學(xué)版), 2013,32(4):562-567.

[8] 暢明琦. 水資源安全理論與方法研究[D]. 西安：西安理工大學(xué), 2006.

[9] 秦長海, 甘 泓, 汪 林, 等. 海河流域水資源開發(fā)利用閾值研究[J]. 水科學(xué)進展, 2013,24(2):220-227.

[10] 蓋 美, 李偉紅. 基于可變模糊識別模型的大連市水資源與社會經(jīng)濟協(xié)調(diào)發(fā)展[J]. 資源科學(xué), 2008,30(8):1 141-1 145.

[11] 羅軍剛, 解建倉, 阮本清. 基于熵權(quán)的水資源短缺風(fēng)險模糊綜合評價模型及應(yīng)用[J]. 水利學(xué)報, 2008,39(9):1 092-1 104.

[12] 劉 迅, 耿進強, 畢遠(yuǎn)志. 基于因子分析法的水利工程標(biāo)段劃分影響因素研究[J]. 中國農(nóng)村水利水電, 2014,(12):91-95.

[13] 楊 娜, 李慧明. 基于因子分析與熵值法的水資源承載力研究----以天津市為例[J]. 軟科學(xué), 2010,24(6):66-70.