耦合多目標(biāo)優(yōu)化和多屬性決策的水庫生態(tài)調(diào)度方法研究

劉 艷， 李曉春， 羅軍剛

(1.西安理工大學(xué) 西北旱區(qū)生態(tài)水利國家重點(diǎn)實(shí)驗(yàn)室， 陜西 西安 710048；2.陜西省江河水庫管理局， 陜西 西安 710018)

水庫生態(tài)調(diào)度是在傳統(tǒng)水庫調(diào)度的基礎(chǔ)上，將水庫和河道的生態(tài)服務(wù)功能考慮進(jìn)來的一種新的調(diào)度模式。目前已成為國內(nèi)外水庫調(diào)度領(lǐng)域的研究熱點(diǎn)。在構(gòu)建生態(tài)調(diào)度模型時，一般是將生態(tài)因子作為約束條件(約束型)或目標(biāo)函數(shù)(目標(biāo)型)。約束型因其目標(biāo)函數(shù)單一，模型求解較為簡單，而被廣泛應(yīng)用。Richter等[1]提出了一個評估水庫生態(tài)調(diào)度的框架，可以用來恢復(fù)河道環(huán)境流量。Castelletti等[2]為保護(hù)河流生態(tài)，將最小生態(tài)流量作為約束方程加入到水庫調(diào)度模型中。Gates等[3]對美國Snake河流上的大壩進(jìn)行生態(tài)調(diào)度研究時，通過增加季節(jié)性枯水流量，進(jìn)而提高下游生物棲息地的利用性。梅超等[4]對不同生態(tài)流量下的調(diào)度方案進(jìn)行研究，分析了水庫群的綜合效益。范驄驤等[5]建立了一種以魚類產(chǎn)卵棲息地的生態(tài)需水為約束條件的生態(tài)調(diào)度方法。目標(biāo)型則是將生態(tài)因子作為模型的目標(biāo)函數(shù)，因其更符合多目標(biāo)問題的本質(zhì)而受到廣泛關(guān)注。Halleraker等[6]通過水力-生境模型、IHA指標(biāo)法、水溫模擬法等多種方法計(jì)算環(huán)境流量，在此基礎(chǔ)上進(jìn)行生態(tài)調(diào)度。尚文繡等[7]通過水文情勢評價指標(biāo)建立了水庫生態(tài)調(diào)度模型。朱金峰等[8]研究了供水期末生態(tài)調(diào)度對生態(tài)用水產(chǎn)生的影響。張珮綸等[9]建立了基于需水調(diào)控策略的最優(yōu)供水模型。張召等[10]以天生橋一級水庫為例，建立了以生態(tài)保證率為生態(tài)目標(biāo)的優(yōu)化調(diào)度模型。賈磊等[11]構(gòu)建了能同時滿足生產(chǎn)生活用水和生態(tài)需水的水庫調(diào)度模型。薛耀東[12]通過對水庫利益主體及其用水特點(diǎn)的分析，構(gòu)建了面向生態(tài)的水庫多利益主體協(xié)調(diào)調(diào)度模型。

目前，水庫生態(tài)調(diào)度研究一直在不斷進(jìn)步和完善。但是大部分將生態(tài)作為約束條件，不能達(dá)到生態(tài)目標(biāo)與其他目標(biāo)之間的均衡，其次在目標(biāo)構(gòu)建時，沒有對目標(biāo)之間進(jìn)行定量分析，直接構(gòu)建模型，對于目標(biāo)的選擇相對不客觀。因此，本文以林家村水庫為研究對象，在構(gòu)建模型之前，對各目標(biāo)之間進(jìn)行相關(guān)分析。為此，本文提出了一種耦合多目標(biāo)相關(guān)分析、多目標(biāo)優(yōu)化和多屬性決策的、貫通“建模-求解-優(yōu)選”全過程的一體化水庫生態(tài)調(diào)度方法，為水庫生態(tài)調(diào)度提供一種新的思路和方法。

1 研究區(qū)概況

林家村水庫位于陜西省寶雞市以西11 km的渭河峽谷出口處，是一座以農(nóng)業(yè)灌溉為主，兼顧防洪發(fā)電和生態(tài)供水的水庫，具體所在位置見圖1。年調(diào)節(jié)水量0.8億m3，為寶雞峽灌區(qū)水庫補(bǔ)水1.48億m3，是陜西省重點(diǎn)水利工程。有四個主要部分：主體工程，壩后電站，防護(hù)工程和庫區(qū)淹沒處理工程。水庫主要特征參數(shù)見表1。

表1 林家村水庫主要特征參數(shù)

圖1 林家村水庫位置Fig.1 Location of Linjiacun Reservoir

渭水被引入渠道灌溉和發(fā)電，有效地緩解了灌區(qū)嚴(yán)重缺水的情況。但是同時也帶來了生態(tài)問題，在枯水期，林家村到魏家堡之間的河道基本處于斷流狀態(tài)，并且生態(tài)需水滿足率在逐年下降。渭河生態(tài)治理進(jìn)一步增加了對生態(tài)用水的需求，使得用水矛盾變得更加突出。因此，為了緩解生態(tài)流量不足、保證生態(tài)用水需求，開展林家村水庫生態(tài)調(diào)度是很有必要的。

2 調(diào)度目標(biāo)提取及關(guān)系分析

2.1 調(diào)度目標(biāo)提取

生態(tài)調(diào)度是在兼顧水庫調(diào)度的社會、經(jīng)濟(jì)效益的基礎(chǔ)上，重點(diǎn)考慮生態(tài)因素的水庫調(diào)度模式[13]。為了合理有效地利用水庫水資源，降低不同生態(tài)流量指標(biāo)對水庫經(jīng)濟(jì)效益產(chǎn)生的影響，本文同時考慮灌溉目標(biāo)、發(fā)電目標(biāo)和生態(tài)目標(biāo)，建立多目標(biāo)優(yōu)化的水庫生態(tài)調(diào)度模型。在構(gòu)建目標(biāo)之前，本文對灌溉和發(fā)電目標(biāo)之間進(jìn)行相關(guān)性分析，如果灌溉和發(fā)電之間是負(fù)相關(guān)關(guān)系，則建立以灌溉、發(fā)電和生態(tài)為目標(biāo)的多目標(biāo)優(yōu)化模型；否則，建立以灌溉效益最大、修正全年流量偏差最小為目標(biāo)函數(shù)的優(yōu)化調(diào)度模型。

目標(biāo)1：灌溉效益最大，可用式(1)表示。

(1)

式中：maxY為最大灌溉效益；Wt為時段灌溉水量；f(Wt)為灌溉水量與灌溉效益的關(guān)系。

本文根據(jù)文獻(xiàn)[14]的研究成果，在不考慮化肥施用折純量、農(nóng)業(yè)從業(yè)人員和農(nóng)用機(jī)械總勞動力變化的情況下，利用面積比擬法，得到寶雞峽塬上灌區(qū)灌溉效益與灌溉水量的效益關(guān)系如下：

f(Wt)=7.16×Wt0.087

(2)

目標(biāo)2：發(fā)電效益最大，可用式(3)表示。

(3)

目標(biāo)3：生態(tài)效益最大，即生態(tài)AAPFD值最小。采用Ladson等提出的修正全年流量偏差函數(shù)(生態(tài)AAPFD值)，相關(guān)研究表明，生態(tài)AAPFD值越小，表示河流生態(tài)越好[15]。

(4)

約束條件主要有水位約束、水量平衡約束、出力約束、灌溉需水量約束、機(jī)組過流能力約束和渠首引水流量約束等。

2.2 調(diào)度目標(biāo)相關(guān)性分析

根據(jù)不同典型年(平水年、較枯水年、特枯水年)來水量和灌溉需水量，以文獻(xiàn)[16]中計(jì)算的生態(tài)流量結(jié)果作為約束條件(見表2)。

表2 林家村生態(tài)流量Tab.2 Ecological flow of Linjiacun

應(yīng)用NSGA-II算法[17]求解以灌溉和發(fā)電為目標(biāo)函數(shù)的雙目標(biāo)調(diào)度模型，探究灌溉效益與發(fā)電效益之間的關(guān)系，結(jié)果見圖2。從圖2中可以看出，平水年，當(dāng)發(fā)電效益小于0.119億元，發(fā)電效益和灌溉效益近似呈線性關(guān)系，y1=7.213x+6.684；當(dāng)發(fā)電效益大于0.119億元，灌溉效益y2=7.725，不再隨著發(fā)電效益的增加而增加。這是因?yàn)楣喔刃杷渴歉鶕?jù)灌溉面積和灌溉定額確定的，不可能是無限大的，當(dāng)灌溉供水量達(dá)到灌溉需水量之后，不會繼續(xù)增加，這個時候發(fā)電引水只會增加發(fā)電效益，對于灌溉效益影響不大。較枯水年和特枯水年發(fā)電效益與灌溉效益近似呈線性正相關(guān)關(guān)系。進(jìn)一步分析發(fā)現(xiàn)，發(fā)電效益與灌溉效益平水年的相關(guān)系數(shù)為0.952，較枯水年為0.991，特枯水年為0.984，三個典型年的發(fā)電效益與灌溉效益均在99%置信水平下顯著正相關(guān)。

圖2 不同典型年灌溉效益與發(fā)電效益關(guān)系Fig.2 Relationship between irrigation benefit and power generation benefit in different typical years

為提高水庫的綜合利用效益，有必要進(jìn)一步分析生態(tài)與灌溉、生態(tài)與發(fā)電目標(biāo)之間的相互關(guān)系。以特枯水年為例，根據(jù)文獻(xiàn)[16]中計(jì)算的生態(tài)流量結(jié)果，分別構(gòu)建灌溉效益最大優(yōu)化模型，模型目標(biāo)函數(shù)見式(1)；以及發(fā)電效益最大優(yōu)化模型，模型目標(biāo)函數(shù)見式(3)。分別分析計(jì)算了生態(tài)流量在不同滿足度的條件下，林家村水庫灌溉效益與林家村斷面生態(tài)流量的響應(yīng)關(guān)系以及發(fā)電效益與生態(tài)流量的響應(yīng)關(guān)系。計(jì)算結(jié)果見圖3。可以看出，隨著生態(tài)流量約束的增加，即年生態(tài)需水量的增加，灌溉效益逐漸減少且呈近似線性遞減的趨勢。林家村水庫灌溉效益的變化趨勢為：

圖3 年生態(tài)需水量與灌溉效益的關(guān)系Fig.3 Relationship between annual ecological water demand and irrigation benefits

y灌=-0.0704q+8.1269

(5)

式中：q為生態(tài)流量(m3/s)；y灌為灌溉效益(億元)?？梢钥闯?，生態(tài)流量每增加1m3/s，灌溉效益就減少0.070 4億元。

同樣地由圖4可以看出，發(fā)電效益也隨著生態(tài)流量約束的增加呈近似線性遞減的趨勢，發(fā)電效益的變化趨勢為：

圖4 年生態(tài)需水量與發(fā)電效益的關(guān)系Fig.4 Relationship between annual ecological water demand and power generation benefits

y電=-0.7685q+9.9415

(6)

式中：q為生態(tài)流量(m3/s)；y電為發(fā)電效益(百萬元)?？梢钥闯?，生態(tài)流量每增加1 m3/s，發(fā)電效益就減少0.768 5百萬元。

通過以上各調(diào)度目標(biāo)之間的相關(guān)性分析表明，林家村水庫發(fā)電與生態(tài)、灌溉與生態(tài)之間均存在負(fù)相關(guān)關(guān)系。林家村水庫的發(fā)電效益和灌溉效益存在正相關(guān)關(guān)系，所以一個目標(biāo)增大的同時，另一個目標(biāo)也同步增大，反之亦然。因此，在構(gòu)建目標(biāo)時，選擇發(fā)電效益和灌溉效益這兩個目標(biāo)其中的一個目標(biāo)即可。林家村水庫主要是以灌溉用水為首要目標(biāo)，發(fā)電居于次要地位，結(jié)合水庫的實(shí)際運(yùn)行情況，由于灌溉放水的不確定性，使得電站無法長期穩(wěn)定的發(fā)電。故本文建立以灌溉效益最大、修正全年流量偏差最小為目標(biāo)函數(shù)，協(xié)調(diào)建立林家村水庫多目標(biāo)生態(tài)調(diào)度模型，并采用平水年、較枯水年、特枯水年3個典型年的徑流資料計(jì)算。

3 多目標(biāo)調(diào)度模型構(gòu)建與求解

3.1 數(shù)學(xué)模型構(gòu)建

依據(jù)林家村水庫的調(diào)度目標(biāo)之間相關(guān)性分析結(jié)果，可構(gòu)建如下的水庫多目標(biāo)調(diào)度模型。

目標(biāo)函數(shù)：

MinimizeF=f(-maxY, minR)，t=1,2,…,T

(7)

約束條件：

1) 水庫水位過程約束:

Zt,min≤Zt≤Zt,max

(8)

2) 時段出力約束:

Nt,min≤Nt≤Nt,max

(9)

3) 機(jī)組過流能力約束:

(10)

4) 水量平衡約束:

Vt+1=Vt+(It-Qt)Δt

(11)

5) 灌溉需水量約束:

(12)

6) 渠首引水流量約束:

(13)

3.2 模型求解算法設(shè)計(jì)

為解決水庫多目標(biāo)優(yōu)化求解-多方案偏好選擇問題，本文將經(jīng)典的多目標(biāo)優(yōu)化算法(NSGA-II)和多屬性決策方法(SEABODE)耦合，建立了一種基于多目標(biāo)優(yōu)化-多屬性決策的水庫生態(tài)調(diào)度模型求解算法(NSGA-II-SEABODE)。該方法的總體思路是：首先，使用NSGA-II算法[17]解決目標(biāo)優(yōu)化問題，以獲得備選方案集A；然后，用SEABODE法在方案集A中優(yōu)選出最終的偏好方案。SEABODE法是一種基于k階p級有效概念的備選方案逐次淘汰法[18]，與傳統(tǒng)方法相比，SEABODE法可直接從決策矩陣中篩選方案，逐步淘汰劣等方案，無需標(biāo)準(zhǔn)化決策矩陣，也無需對屬性權(quán)重進(jìn)行計(jì)算，因此，降低了方案在優(yōu)選過程中受到的主觀影響。NSGA-II-SEABODE算法的求解流程見圖5。

圖5 NSGA-II-SEABODE算法流程圖Fig.5 NSGA-II-SEABODE combined method flow chart

Step1： 利用NSGA-II算法對多目標(biāo)優(yōu)化問題進(jìn)行求解，得到備選方案集A。

Step2： 基于SEABODE法構(gòu)建四維屬性空間C={α，γ，υ，MSI}對方案集A排序。其中，α為可靠性，γ為可恢復(fù)性，ν為脆弱性，MSI為缺水指數(shù)。

1) 從k=4開始，在四維屬性空間C={α，γ，υ，MSI}中辨別4-Pareto最優(yōu)方案，找出所有最低階(記為kmin)。

2) 從選出的四階有效方案集中進(jìn)一步向三階有效方案進(jìn)行優(yōu)選，即k=3，此時三階子空間分別由{α，γ，υ}、{α，γ，MSI}、{γ，υ，MSI}、{α，υ，MSI}構(gòu)成。

3) 若2)找到的所有kmin-Pareto最優(yōu)方案中沒有一個是(kmin-1)-Pareto最優(yōu)的，則選擇的偏好方案為(kmin-1)階最高級(pmax)有效方案。

Step3： 重復(fù)3)，直到找出[k,p]-Pareto最優(yōu)方案為偏好方案。

3.3 生態(tài)調(diào)度方案合理性評價

基于NSGA-II法對水庫多目標(biāo)生態(tài)調(diào)度模型求解，得到含有100個非劣解的非劣解集，為了找尋綜合效益最大的方案，需要對多目標(biāo)方案進(jìn)行優(yōu)選。因此為了進(jìn)一步淘汰部分較差方案，本文選擇可靠性、可恢復(fù)性、脆弱性和缺水指數(shù)四個評價指標(biāo)構(gòu)成水庫調(diào)度合理性評價指標(biāo)體系的四維屬性空間C={α，γ，υ，MSI},其中α，γ為最大化類型屬性指標(biāo)；ν和MSI為最小化類型屬性指標(biāo)。然后，采用SEABODE方法從備選方案集A的100個方案中篩選出最終的偏好方案。

1) 可靠性α。在水庫調(diào)度期內(nèi)，達(dá)到需求目標(biāo)的概率，其計(jì)算公式為：

(14)

2) 可恢復(fù)性γ。在運(yùn)行分析期內(nèi)，水庫從破壞狀態(tài)恢復(fù)到正常供水的平均頻率，其計(jì)算公式為：

(15)

式中：γ為水庫供水可恢復(fù)性。

3) 缺水深度ν。水庫運(yùn)行分析期內(nèi)，單個時段內(nèi)的最大相對缺水量，其計(jì)算公式為：

ν=max{DR1,DR2,…,DRT}

(16)

4) 缺水指數(shù)MSI。反映水庫供水效益的損失程度。

(17)

式中：T為調(diào)度期總時段數(shù)。

4 實(shí)例計(jì)算與分析

4.1 典型年的選擇

選取1950—2018年共69年的水文序列，將其進(jìn)行排頻計(jì)算，然后根據(jù)Pearson-Ⅲ型曲線對水庫年平均流量進(jìn)行適線，選取P=5%為特豐水年，P=25%為較豐水年，P=50%為平水年，P=75%為較枯水年，P=90%為特枯水年。典型年的來水過程見圖6。

圖6 林家村水庫典型年入庫流量Fig.6 Storage water of typical year of Linjiacun Reservoir

4.2 灌溉需水量計(jì)算

寶雞峽水庫的主要供水目標(biāo)包括灌溉、發(fā)電和生態(tài)用水。根據(jù)文獻(xiàn)[19]可知，特豐水年和較豐水年，渠首引水后下泄流量能夠滿足生態(tài)流量；平水年引水后生態(tài)流量不足；而較枯水年和特枯數(shù)年引水前生態(tài)基流已不足，引水后更加嚴(yán)重。故本文主要計(jì)算平水年、較枯水年和特枯水年的灌溉需水量。

根據(jù)寶雞峽灌區(qū)的灌溉制度和灌溉定額，采用定額指標(biāo)法，通過分析計(jì)算可以得到寶雞峽塬上灌區(qū)的月灌溉需水量，結(jié)果見圖7。其中50%典型年灌溉需水量為2.43億m3，75%典型年灌溉需水量為5.55億m3，90%典型年灌溉需水量為5.91億m3。

圖7 不同典型年灌溉需水量Fig.7 Irrigation water requirement in different typical years

4.3 優(yōu)化求解與結(jié)果分析

采用MATLAB編制水庫生態(tài)調(diào)度的NSGA-Ⅱ-SEABODE優(yōu)化求解算法，以月為調(diào)度時段，總調(diào)度時段12個，運(yùn)行期為7月至次年6月。優(yōu)化問題包括13個決策變量，分別是水庫在12個時間段內(nèi)的13個水位值(包括起始水位和結(jié)束水位)。選取種群規(guī)模N=100，最大迭代次數(shù)Gen=1 000，NSGA-II算法使用模擬的二進(jìn)制交叉算子和多項(xiàng)式突變，交叉概率pc=0.9，變異概率pm=1/n，n是決策變量的個數(shù)，交叉分配指數(shù)ηc和突變分布指數(shù)ηm都是20。通過計(jì)算可得到100組Pareto最優(yōu)解。不同典型年的優(yōu)化調(diào)度結(jié)果見圖8，可以看出灌溉效益與生態(tài)AAPFD值之間存在著競爭關(guān)系。為進(jìn)一步分析優(yōu)化結(jié)果，取其中的兩種典型優(yōu)化方案：方案1(生態(tài)AAPFD值最小)、方案2(灌溉效益最大)。當(dāng)以追求灌溉效益最大為目標(biāo)時，可選擇方案2；當(dāng)追求生態(tài)AAPFD值最小為目標(biāo)時，可選擇方案1。從方案1到方案2，水庫的灌溉效益在增加的同時，生態(tài)AAPFD值也在增加?？梢姙榱吮Ｕ虾拥纼?nèi)具有一定的生態(tài)流量，會使水庫的經(jīng)濟(jì)效益受到一定的影響。

圖8 不同典型年多目標(biāo)Pareto散點(diǎn)圖Fig.8 Pareto scatter diagram of multi-objective in different typical years

分別計(jì)算平水年、較枯水年和特枯水年水庫屬性空間C中4個評價指標(biāo)值，統(tǒng)計(jì)結(jié)果見表3。

由表3可以看出，三個不同典型年中，不同方案的4個評價指標(biāo)均存在差異。

表3 決策空間A中不同典型年評價指標(biāo)統(tǒng)計(jì)結(jié)果Tab.3 Statistical results of evaluation indicators of different typical years in decision space A

因此，可采用SEABODE法對方案集進(jìn)行選擇與淘汰，確定最終的偏好方案。

平水年、較枯水年、特枯水年三種典型年下，四維指標(biāo)空間C及其三維子空間中的Pareto最優(yōu)方案數(shù)統(tǒng)計(jì)結(jié)果見表4。

表4 不同典型年的指標(biāo)空間C中Pareto最優(yōu)方案個數(shù)Tab.4 Number of Pareto optimal schemes in indicator space C of different typical years

表4 不同典型年的指標(biāo)空間C中Pareto最優(yōu)方案個數(shù)Tab.4 Number of Pareto optimal schemes in indicator space C of different typical years

典型年(1234)(123)(124)(134)(234)三階有效方案平水年4639301887較枯水年564229201513特枯水年49413613129

注：1-α，2-γ，3-ν，4-MSI。

在k=4的第一輪方案優(yōu)劣排序識別中，不同典型年的方案A1～方案A100，在四維指標(biāo)空間C中分別有46、56、49個方案是Pareto最優(yōu)的，相當(dāng)于偏好方案的優(yōu)選范圍分別縮小了54%、44%、51%。為了選擇出一個具有代表性的偏好方案，需進(jìn)一步完善和縮小決策者的選擇范圍。在第二輪中，將所選的四階有效解進(jìn)一步優(yōu)化為三階有效解，分別獲得平水年的三階有效方案個數(shù)是7個，較枯水年是13個，特枯水年是9個。第三輪k=2的方案優(yōu)劣排序，同樣使用SEABODE法進(jìn)行優(yōu)選，得到了平水年(P=50%)、較枯水年(P=75%)和特枯水年(P=90%)下，多屬性決策的最終偏好方案分別為方案A29、方案A13、方案A80，結(jié)果見表5。

表5 典型年偏好方案結(jié)果Tab.5 Preference schemes for three typical years

從表5可以看出，不同典型年的生態(tài)調(diào)度性能評價指標(biāo)差別很大。平水年的供水可靠性、可恢復(fù)性、脆弱性指標(biāo)均優(yōu)于較枯水年和特枯水年。其中，較枯水年和特枯水年的調(diào)水可靠性α分別下降了12.4 %、24.7%；可恢復(fù)性γ分別下降了25.4 %、50.7%；脆弱性ν分別下降了25.0%、33.3%。較枯水年的缺水指數(shù)MSI最小，為2.21，表明12月中有2個多月無法滿足正常需水量。平水年、較枯水年、特枯水年的經(jīng)濟(jì)效益分別為7.949、8.452、7.927億元，生態(tài)AAPFD值分別為1.624、1.207、1.425。

5 結(jié) 論

為解決水庫生態(tài)調(diào)度模型構(gòu)建和求解中存在的問題，本文提出并構(gòu)建了一種耦合多目標(biāo)相關(guān)分析、多目標(biāo)優(yōu)化和多屬性決策的水庫生態(tài)調(diào)度方法，并通過實(shí)例分析與驗(yàn)證可得出如下結(jié)論。

1) 通過抽取水庫生態(tài)調(diào)度目標(biāo)，并進(jìn)行相關(guān)性分析，從而獲得相互沖突的調(diào)度目標(biāo)，使得水庫生態(tài)調(diào)度模型的構(gòu)建更具科學(xué)性和定量依據(jù)。

2) 將多目標(biāo)優(yōu)化和多屬性決策方法結(jié)合，構(gòu)建的NSGA-II-SEABODE算法，能夠在優(yōu)化求解模型的同時，從Pareto最優(yōu)調(diào)度方案中優(yōu)選出最終偏好調(diào)度方案。

3) 將建模方法、優(yōu)化算法、優(yōu)選方法耦合，從而為水庫生態(tài)調(diào)度提供了一種貫通“模型構(gòu)建-優(yōu)化求解-方案優(yōu)選”全過程的調(diào)度方法。