基于EVT法和MC模型的極限水文干旱歷時(shí)計(jì)算研究

趙 緯，宋松柏，張更喜

(西北農(nóng)林科技大學(xué) 水利與建筑工程學(xué)院，陜西 楊凌712100)

干旱是世界上最主要的自然災(zāi)害之一，一直受到環(huán)境學(xué)、生態(tài)學(xué)、水文學(xué)和氣象學(xué)等諸多領(lǐng)域的密切關(guān)注。由于干旱形成條件的復(fù)雜性和對(duì)社會(huì)影響的廣泛性，干旱一般可分為水文干旱、氣象干旱、農(nóng)業(yè)干旱和社會(huì)經(jīng)濟(jì)干旱。其中，水文干旱的影響因素較氣象、農(nóng)業(yè)等干旱更加復(fù)雜，且更容易受到人為因素的影響。因而，其能較全面地反映區(qū)域的整體干旱狀況[1]。水文干旱通常是指在一定時(shí)期內(nèi)河流徑流量低于臨界水平而造成的水分虧缺現(xiàn)象。其中，虧缺量和持續(xù)時(shí)間可以分別用干旱烈度和干旱歷時(shí)表示[2]。當(dāng)水文干旱歷時(shí)持續(xù)的時(shí)間超過河流、水庫(kù)等的承載能力時(shí)，就會(huì)對(duì)農(nóng)業(yè)生產(chǎn)和社會(huì)發(fā)展造成嚴(yán)重影響[3]。由于現(xiàn)有的實(shí)測(cè)水文序列較短，從而使得在有限實(shí)測(cè)序列中提取出的干旱序列更短。因此，上述問題給區(qū)域?qū)嶋H旱情和抗旱年限的確定帶來了困難。在實(shí)際應(yīng)用中，人們常通過嚴(yán)密的數(shù)學(xué)推導(dǎo)獲取表征極限水文干旱歷時(shí)的解析表達(dá)式，從而為區(qū)域水資源規(guī)劃與利用、抗旱年限確定以及區(qū)域社會(huì)與經(jīng)濟(jì)發(fā)展計(jì)劃制定等提供依據(jù)[4]。

極限干旱歷時(shí)是指在一定時(shí)期內(nèi)可能出現(xiàn)的水文干旱的最大干旱歷時(shí)[5]，其計(jì)算方法是水文學(xué)家和水資源管理者共同感興趣的問題之一。20世紀(jì)60年代末，F(xiàn)eller[6]應(yīng)用游程理論首次推導(dǎo)出最大游程長(zhǎng)度的概率分布函數(shù)。之后，關(guān)于極限干旱歷時(shí)概率計(jì)算方法的研究得到了國(guó)內(nèi)外學(xué)者的廣泛關(guān)注。Moye等[7-8]以游程理論為基礎(chǔ)，運(yùn)用極值順序統(tǒng)計(jì)量概念，經(jīng)過嚴(yán)密的數(shù)學(xué)推導(dǎo)得到了計(jì)算極限干旱歷時(shí)期望值和方差值的解析表達(dá)式。馬秀峰等[9-10]從簡(jiǎn)單的伯努利試驗(yàn)出發(fā)，結(jié)合隨機(jī)模擬試驗(yàn)，發(fā)現(xiàn)平穩(wěn)隨機(jī)序列的輪次長(zhǎng)度服從二階線性齊次差分方程。馮國(guó)章[11-13]用枚舉法推導(dǎo)出獨(dú)立同分布隨機(jī)變量序列最大負(fù)游程長(zhǎng)度的概率分布函數(shù)，并將其作為極限干旱歷時(shí)的概率分布函數(shù)，與相同結(jié)構(gòu)的年徑流序列最大水文干旱歷時(shí)的經(jīng)驗(yàn)頻率擬合，取得了較好的擬合效果。Sharma等[14-15]研究了不同截?cái)嗨綏l件下，當(dāng)徑流序列服從正態(tài)分布、gamma分布和對(duì)數(shù)正態(tài)分布的極限干旱歷時(shí)和強(qiáng)度的概率分布函數(shù)。丁晶等[16]以年徑流數(shù)據(jù)為基礎(chǔ)，分析水文干旱中負(fù)輪長(zhǎng)的變化特征，發(fā)現(xiàn)P-Ⅲ型分布可以用來表征負(fù)輪長(zhǎng)的統(tǒng)計(jì)特性。Cancelliere等[17]假設(shè)水文數(shù)據(jù)服從周期性的簡(jiǎn)單馬爾科夫鏈，推導(dǎo)出了極限干旱歷時(shí)概率分布和重現(xiàn)期的解析表達(dá)式。袁超等[4]根據(jù)游程理論和隨機(jī)變量極值理論，采用解析法和模擬法進(jìn)行了極限水文干旱歷時(shí)概率分布和徑流序列參數(shù)的影響研究。Sharma等[18]和Panu等[19]以月、年時(shí)間尺度上的徑流序列為例，運(yùn)用極值理論法和馬爾科夫模型研究了極限干旱歷時(shí)、最大干旱強(qiáng)度和最大干旱烈度之間的關(guān)系。綜上可知，關(guān)于水文干旱歷時(shí)概率的研究雖已經(jīng)取得了一些成果，但對(duì)于極值理論(Extreme value theory,EVT)法和馬爾科夫鏈(Markov chain,MC)模型在計(jì)算極限干旱歷時(shí)中的對(duì)比研究尚比較少，且缺少應(yīng)用案例。為此，本研究應(yīng)用游程理論(Run theory)，以黃河流域3個(gè)水文站的年徑流資料為研究對(duì)象，應(yīng)用EVT法和MC模型計(jì)算極限干旱歷時(shí)，將計(jì)算值與實(shí)際觀測(cè)結(jié)果進(jìn)行對(duì)比，并采用納什效率系數(shù)(NS)和平均絕對(duì)百分誤差(MAPE)對(duì)計(jì)算精度進(jìn)行評(píng)價(jià)，分析EVT法和MC模型在該流域的適用性，以期為區(qū)域干旱歷時(shí)的預(yù)報(bào)提供參考和依據(jù)。

1 極限干旱歷時(shí)的計(jì)算方法

1.1 EVT法

隨機(jī)變量極值理論是Todorovic等[20]在1975年提出的，該方法有以下3點(diǎn)假設(shè)：

1) 在重現(xiàn)期T內(nèi)，干旱發(fā)生次數(shù)的概率服從泊松分布，其概率密度函數(shù)為：

(1)

式中：P為離散事件的概率值；n為干旱發(fā)生次數(shù)；j為干旱發(fā)生次數(shù)值，取離散值；q為干旱概率；qq為前一時(shí)段發(fā)生干旱條件下下一時(shí)段發(fā)生干旱的條件概率。

2) 對(duì)干旱歷時(shí)進(jìn)行分析時(shí)，需要用到條件概率qq。qq可由求解徑流序列的相關(guān)計(jì)算參數(shù)得到，計(jì)算公式[21]為：

(2)

式中：ρ1為一階自相關(guān)系數(shù)，z0為標(biāo)準(zhǔn)正態(tài)變量值，ρ為積分變量。

3) 在某一游程中，干旱歷時(shí)發(fā)生的概率一般服從幾何分布，即：

P(L=l)=ql-1·p,l=1,2,3,…。

(3)

式中：L為干旱歷時(shí)；l為干旱歷時(shí)值變量，取離散值；p為不發(fā)生干旱的概率，p=1-q。

基于以上3點(diǎn)假設(shè)，可以得出在重現(xiàn)期T內(nèi)，極限干旱歷時(shí)的概率分布、期望及方差，即：

P(LT=i)=exp[-T·q(1-qq)qqi-1]·
[exp(T·q(1-qq)2qqi-1)-1],i=1,2,3,…。

(4)

(5)

(6)

式中：LT為極限干旱歷時(shí)。對(duì)于年徑流序列，i的最大取值為25，因?yàn)楫?dāng)i的取值大于25時(shí)，式(4)中的概率值很小可以忽略不計(jì)。因此，式(5)又可以表示為：

E(LT)=1P(LT=1)+2P(LT=2)+3P(LT=3)+
…+25P(LT=25)。

(7)

1.2 MC模型

根據(jù)Sharma等[22]對(duì)于極限水文干旱歷時(shí)的定義，馬爾科夫鏈中有關(guān)干旱事件的概率可約定如下：

q=P(d);p=P(w),

(8)

qp=P(d/w);pqq=P(w/d,d),

(9)

qqp=P(d|d,w);qqq=P(d|d,d)。

(10)

式中：d為干旱事件；w為不干旱事件；qp為前一時(shí)段不發(fā)生干旱條件下下一時(shí)段發(fā)生干旱的條件概率；pqq為前兩時(shí)段發(fā)生干旱條件下下一時(shí)段不發(fā)生干旱的條件概率；qqp為前兩時(shí)段中第一時(shí)段發(fā)生干旱而第二時(shí)段不發(fā)生干旱條件下下一時(shí)段不發(fā)生干旱的條件概率；qqq為前兩時(shí)段發(fā)生干旱的條件下下一時(shí)段仍然發(fā)生干旱的條件概率。

根據(jù)干旱的定義，可以將干旱歷時(shí)的發(fā)生概率表示為：

P(L=l)=P(w)P(d|w)P(d|d,w)·

[P(d|d,d)]l-2P(w|d,d),l=2,3,4,…。

(11)

在重現(xiàn)期T內(nèi)，干旱發(fā)生的超越概率為1/T。如果干旱歷時(shí)的發(fā)生概率≥l，則在滿足二階馬爾科夫鏈(two-order Markov chain,MC-2)模型的條件下，公式(11)可進(jìn)一步表示為：

(12)

將式(8)～(10)代入式(12)，可得重現(xiàn)期為T年的超越概率為：

(13)

由式(13)可得，滿足MC-2模型的極限干旱歷時(shí)的表達(dá)式為：

(14)

本研究中極限干旱歷時(shí)的計(jì)算選用Sharma等[15]修正后的公式。因此，滿足MC-2模型的極限水文干旱歷時(shí)的表達(dá)式為：

(15)

式中：F為計(jì)算因子，可以通過經(jīng)驗(yàn)繪點(diǎn)位置公式得到，例如F=1，1.6，2分別代表weibull、blom、hazen繪點(diǎn)位置公式。

當(dāng)qqp=qq且qqq=qq時(shí)，MC-2模型可簡(jiǎn)化為一階馬爾科夫鏈(one-order Markov chain,MC-1)模型。因此，滿足MC-1模型的極限水文干旱歷時(shí)的表達(dá)式為：

(16)

當(dāng)qp=q和qq=q時(shí)，MC-1模型可簡(jiǎn)化為零階馬爾科夫鏈(zero-order Markov chain,MC-0)模型。因此，滿足MC-0模型的極限水文干旱歷時(shí)的表達(dá)式為：

(17)

由于年徑流序列具有弱相關(guān)性,因此本研究選取MC-1模型計(jì)算極限干旱歷時(shí)。其中條件概率qp可由式(18)和式(19)得到，即：

(18)

qp+pp=1。

(19)

式中：pp為前一時(shí)段不發(fā)生干旱條件下下一時(shí)段也不發(fā)生干旱的條件概率。

2 徑流序列的標(biāo)準(zhǔn)化處理

對(duì)于水文序列xi，經(jīng)過標(biāo)準(zhǔn)化處理后得到標(biāo)準(zhǔn)水文序列ei，其轉(zhuǎn)化公式為：

(20)

式中：μx和σx分別為水文序列的均值和標(biāo)準(zhǔn)差。

如果標(biāo)準(zhǔn)水文序列ei服從正態(tài)分布，則有ei=zn，那么干旱概率q可以通過標(biāo)準(zhǔn)化正態(tài)分布函數(shù)表或者下面的標(biāo)準(zhǔn)化正態(tài)概率積分得到,即：

(21)

式中：e0為標(biāo)準(zhǔn)截?cái)嗨?，e0=(x0-μ0)/σx，其中x0為截?cái)嗨?；zn為標(biāo)準(zhǔn)正態(tài)變量值；z為積分變量。

當(dāng)標(biāo)準(zhǔn)水文序列ei服從對(duì)數(shù)正態(tài)分布或兩參數(shù)gamma分布時(shí)，可以通過公式(22)或(23)變換為近似標(biāo)準(zhǔn)化的正態(tài)序列zl或zg，然后再將zl=zn或zg=zn代入式(21)中，從而得到干旱概率q。其中：

(22)

(23)

式中：Cv為水文序列的變差系數(shù)。

在本研究中，可以利用公式(24)中的Wilson-Hilferty反轉(zhuǎn)換法[23]，將標(biāo)準(zhǔn)水文序列ei變換為近似標(biāo)準(zhǔn)化的正態(tài)序列z0，然后再將z0=zn代入式(21)中，從而得到干旱概率q。

Z0=(6/Cs)[(0.5Csei+1)0.333-1]+Cs/6。

(24)

式中：Cs為水文序列的偏態(tài)系數(shù)。

取不同的截?cái)嗨絜0和變差系數(shù)Cv，由式(2)、(21)、(22)、(23)分別計(jì)算正態(tài)分布、對(duì)數(shù)正態(tài)分布和兩參數(shù)gamma分布的標(biāo)準(zhǔn)正態(tài)值和干旱概率值，結(jié)果見表1。

表1 不同截?cái)嗨较碌臉?biāo)準(zhǔn)正態(tài)值和干旱概率值Table 1 Standard normal values drought probability in different threshold levels

3 實(shí)例應(yīng)用

3.1 黃河流域極限水文干旱歷時(shí)的確定

以黃河流域蘭州、龍門、白馬寺3個(gè)水文站的年徑流序列為例，應(yīng)用上述公式，計(jì)算基于EVT法和MC-1模型的極限干旱歷時(shí)。經(jīng)審查，各站資料均符合可靠性、一致性和代表性的要求，3個(gè)水文站徑流序列的基本統(tǒng)計(jì)參數(shù)，包括均值μx、標(biāo)準(zhǔn)差σx、變差系數(shù)Cv、偏態(tài)系數(shù)Cs和一階自相關(guān)系數(shù)ρ1均列于表2。

表2 黃河流域3個(gè)水文站年徑流序列的統(tǒng)計(jì)參數(shù)Table 2 Statistical parameters of 3 hydrological stations in the Yellow River Basin

應(yīng)用EVT法和MC-1模型，計(jì)算3個(gè)水文站在10個(gè)概率水平下的極限干旱歷時(shí)，并將E(LT)與LT-ob進(jìn)行對(duì)比。各測(cè)站年徑流序列E(LT)和LT-ob的計(jì)算結(jié)果見表3～5。從表3～5可以看出：當(dāng)干旱概率q以0.05為步長(zhǎng)均勻減小時(shí)，標(biāo)準(zhǔn)化截?cái)嗨絜0、干旱概率q和條件概率qq、qp均呈減小趨勢(shì)。與LT-ob相比，EVT法計(jì)算的E(LT)的偏差值在1年之內(nèi)；而MC-1模型結(jié)合weibull、hazen、和blom繪點(diǎn)位置公式計(jì)算的E(LT)，在有些概率水平下的偏差值達(dá)到了2年。

為了更加直觀地對(duì)比分析2種計(jì)算方法的優(yōu)劣，將E(LT)和LT-ob繪制在斜率為1的直線兩側(cè)。由于MC-1模型結(jié)合3種繪點(diǎn)位置公式計(jì)算得到的E(LT)和LT-ob的對(duì)比圖趨勢(shì)相似，因此本研究只將MC-1模型結(jié)合hazen繪點(diǎn)位置公式計(jì)算的E(LT)和LT-ob對(duì)比繪制在圖1中。

由圖1可以看出，各測(cè)站年徑流序列的E(LT)和LT-ob點(diǎn)值分散在斜率為1的直線周圍，這是因?yàn)閷?duì)計(jì)算值進(jìn)行了取整。例如：在年尺度上計(jì)算的E(LT)是4.49和4.51，分別被近似為4和5年。當(dāng)計(jì)算值的范圍在1～11時(shí)，即便是180個(gè)點(diǎn)被計(jì)算出來，但大量的點(diǎn)都是彼此重疊的，并且只有有限的點(diǎn)被繪制在圖形中。在不同的截?cái)嗨较拢珽VT法計(jì)算的E(LT)與LT-ob的偏差在2年之內(nèi)，而MC-1模型結(jié)合hazen繪點(diǎn)位置公式計(jì)算的E(LT)與LT-ob的偏差值在有些截?cái)嗨较聞t達(dá)到了3年。

表3 蘭州站極限干旱歷時(shí)計(jì)算值和對(duì)應(yīng)觀測(cè)值的比較Table 3 Comparison of calculated and observed values of critical drought duration at Lanzhou station

注：第6～9列括號(hào)內(nèi)數(shù)值為計(jì)算值，括號(hào)外數(shù)值為取整值；表4，5同。

Note:The values in the parentheses are calculated values and outside are the integer values from sixth to ninth columns.The same for Table 4 and 5.

表4 龍門站極限干旱歷時(shí)計(jì)算值和對(duì)應(yīng)觀測(cè)值的比較Table 4 Comparison of calculated and observed values of critical drought duration at Longmen station

表5 白馬寺站極限干旱歷時(shí)計(jì)算值和對(duì)應(yīng)觀測(cè)值的比較Table 5 Comparison of calculated values and observed values of critical drought duration for Baimasi station

圖1 3個(gè)水文站極限干旱歷時(shí)計(jì)算值與觀測(cè)值的對(duì)比Fig.1 Comparison of calculated and observed values of critical drought duration at 3 hydrological stations

3.2 計(jì)算結(jié)果的精度評(píng)價(jià)

基于年時(shí)間尺度,通過MC-1模型和EVT 2種方法計(jì)算的極限干旱歷時(shí)，很難從數(shù)值上判斷哪種方法更加優(yōu)越，故需要用NS和MAPE 2個(gè)統(tǒng)計(jì)量對(duì)計(jì)算結(jié)果進(jìn)行評(píng)價(jià)。NS值越接近1，或MAPE值越接近0，表明預(yù)測(cè)效果越準(zhǔn)確。

(1)納什效率系數(shù)(NS)。計(jì)算公式為：

(2)平均絕對(duì)百分誤差(MAPE)。計(jì)算公式為：

(26)

式中：Pi為計(jì)算值；Oi為實(shí)測(cè)值；Omean為實(shí)測(cè)值的平均值；N為實(shí)測(cè)干旱歷時(shí)長(zhǎng)度。

本研究將黃河流域蘭州、龍門、白馬寺3個(gè)水文站極限干旱歷時(shí)計(jì)算精度的評(píng)價(jià)結(jié)果列于表6。

表6 3個(gè)水文站極限干旱歷時(shí)計(jì)算結(jié)果的精度評(píng)價(jià)Table 6 Accuracy evaluation of predictive results for critical drought duration at 3 stations

從表6可以看出：對(duì)于蘭州站而言，在年時(shí)間尺度上,EVT法和MC-1模型計(jì)算結(jié)果的NS值均超過0.824 7，并且兩者很接近；其中EVT法計(jì)算的NS值達(dá)到了0.921 0，可認(rèn)為其計(jì)算精度最高。同時(shí)2種方法計(jì)算結(jié)果的MAPE值也以EVT法最小。對(duì)于龍門站而言，MC-1模型結(jié)合3種繪點(diǎn)位置公式所得結(jié)果的NS值和MAPE值差異均比較大，其中MC-1模型結(jié)合weibull繪點(diǎn)位置公式時(shí)的精度評(píng)價(jià)結(jié)果與EVT法最為接近，只是EVT法的MAPE值略小于MC-1模型結(jié)合weibull繪點(diǎn)位置公式的值。對(duì)于白馬寺站而言，EVT法和MC-1模型結(jié)合blom繪點(diǎn)位置公式所得結(jié)果的NS值和MAPE值較為接近，二者無明顯差異。因此，從以上3個(gè)水文站的精度評(píng)價(jià)可以看出：在年徑流時(shí)間尺度上，與MC-1模型相比，EVT法能夠更好地用于極限干旱歷時(shí)的計(jì)算。

4 結(jié) 論

通過EVT法和MC-1模型來計(jì)算黃河流域的極限水文干旱歷時(shí)，主要得到以下幾點(diǎn)結(jié)論：

1)黃河流域蘭州、龍門、白馬寺3個(gè)水文站年徑流序列的計(jì)算表明，利用MC-1模型結(jié)合weibull、hazen和bolm 3種繪點(diǎn)位置公式計(jì)算得到的極限干旱歷時(shí)結(jié)果較為接近，在有些截?cái)嗨较掠?jì)算得到的E(LT)與LT-ob的偏差達(dá)到了3年，而利用EVT法計(jì)算的E(LT)與LT-ob的偏差在2年之內(nèi)。

2)用NS和MAPE對(duì)MC-1模型和EVT法的計(jì)算結(jié)果進(jìn)行精度評(píng)價(jià)，結(jié)果表明，在年徑流時(shí)間尺度上，對(duì)于黃河流域極限干旱歷時(shí)的計(jì)算，EVT法優(yōu)于MC-1模型。

3)EVT法是干旱歷時(shí)發(fā)生的泊松概率法則和干旱長(zhǎng)度的幾何概率法則的整合，并且公式中只涉及3個(gè)參數(shù)，即干旱概率q、條件概率qq和重現(xiàn)期T。由于EVT法只需要獲得q和T就可以求得極限干旱歷時(shí)，且計(jì)算精度較高，因此可以作為極限干旱歷時(shí)的求解方法。