王志成

(新疆塔里木河流域管理局， 新疆 庫爾勒 841000)

基于GEV模型的開都河上游降水模擬

王志成

(新疆塔里木河流域管理局， 新疆 庫爾勒 841000)

基于開都河流域上游降水?dāng)?shù)據(jù)，采用GEV模型進行極值降水的模擬和分析。結(jié)果表明：AMS1序列的形狀參數(shù)大于0，服從Fréchet分布；在較高重現(xiàn)期下由輪廓似然方法估計的置信區(qū)間比Delta方法更準確；重現(xiàn)水平的輪廓似然函數(shù)曲線在較高重現(xiàn)期之下呈較顯著不對稱性；基于非平穩(wěn)GEV模型得到極值降水設(shè)計值，其在1958年的100年一遇設(shè)計值到2010年下降為接近50年一遇，預(yù)示著未來發(fā)生極值降水和洪災(zāi)的風(fēng)險加大。

GEV模型； 極值降水； 開都河上游； 非平穩(wěn)

1 引 言

全球氣候變化背景下，水文氣象極值事件增加的可能性已經(jīng)引起全世界各國政府和學(xué)術(shù)界的關(guān)注，因為相比于平均值的變化，極值事件頻率和強度的變化對自然和人類社會將產(chǎn)生更加意義深遠的影響。極值水文氣象事件對于全球氣候變化的響應(yīng)十分敏感，已有研究表明，當(dāng)前世界上許多國家和地區(qū)的極端水文氣象事件發(fā)生概率在增加。此外，在氣候變化和人類活動影響越來越顯著的條件下，極值水文氣象數(shù)據(jù)常常表現(xiàn)出非平穩(wěn)性，如存在顯著趨勢或發(fā)生突變等[1]，針對具有非平穩(wěn)特征的水文氣象序列如何進行模擬，已經(jīng)引起了廣泛的關(guān)注。

概化極值分布(Generalized Extreme Value，GEV)將Weibull分布、Gumbel分布和Fréchet分布統(tǒng)一起來，不必考慮原始分布類型，可以避免單獨采用某一種分布的不足，因而被廣泛用于極值事件的模擬。值得注意的是，在極值模擬中，GEV分布所依據(jù)的極值理論假定序列具有平穩(wěn)性，然而水文氣象時間序列平穩(wěn)性的假定一直受到質(zhì)疑[2]。葉長青等[3]研究表明，當(dāng)洪水序列存在非一致性時，如果仍采用傳統(tǒng)方法，會高估設(shè)計洪水的量級。針對此問題，謝平等[4]提出了基于時間序列分析的非一致性年徑流序列水文頻率計算方法；宋松柏等[5]推導(dǎo)了具有跳躍變異的非一致分布水文序列頻率計算公式。

當(dāng)前基于GEV模型進行開都河流域上游極值降水模擬方面的研究尚不多見，尤其是針對該流域極值降水序列可能存在的非平穩(wěn)特征進行模擬的研究尚未開展。本文選擇開都河上游巴音布魯克氣象站數(shù)據(jù)，基于GEV模型進行該流域1958—2010年極值降水演變的模擬和分析，以期有助于該區(qū)域未來洪災(zāi)風(fēng)險的評估。

2 數(shù)據(jù)來源和研究方法

2.1 數(shù)據(jù)來源

數(shù)據(jù)來源于中國氣象局國家氣象信息中心(中國氣象數(shù)據(jù)網(wǎng)http://cdc.cma.gov.cn)。選取流域內(nèi)的巴音布魯克氣象站為典型站，逐日降水?dāng)?shù)據(jù)序列跨度為1958—2010年，提取得到逐年最大1d降水量序列(annual maximum series 1-day,AMS1)，進一步進行極值降水序列的模擬。

2.2 非平穩(wěn)GEV模型

平穩(wěn)GEV模型，其位置參數(shù)、尺度和形狀參數(shù)是常數(shù)，不隨時間而變化；而非平穩(wěn)GEV模型的參數(shù)可以隨時間而變化，其關(guān)系可以是線性函數(shù)，也可以是二次函數(shù)等，本文采用線性函數(shù)。位置參數(shù)和尺度參數(shù)對極值降水的變化趨勢起主要作用，而形狀參數(shù)的作用則較小，另外，模擬中如允許形狀參數(shù)進行變化，可能造成數(shù)值計算方面的問題。因此本文將形狀參數(shù)作為常數(shù)。此外，模擬時一般先考慮位置參數(shù)隨時間變化，在此基礎(chǔ)上再考慮尺度參數(shù)隨時間的變化。

2.3 Delta方法與輪廓似然函數(shù)方法

置信區(qū)間的估計方法有Delta方法與輪廓似然函數(shù)方法。其中Delta方法假定參數(shù)的估計結(jié)果是對稱的，但是對于形狀參數(shù)和重現(xiàn)水平來說往往并不如此，該方法所得置信區(qū)間關(guān)于重現(xiàn)水平的估計值也是對稱的，這在極值統(tǒng)計的應(yīng)用時存在不足。輪廓似然函數(shù)方法是在僅對極值分布中的一部分參數(shù)感興趣的情況下，構(gòu)造感興趣參數(shù)的置信區(qū)間。研究表明，輪廓似然函數(shù)曲線的不對稱性，正好與極值變量取大值時的不確定性大于取小值時的不確定性這一特點相吻合，因此該方法的區(qū)間估計結(jié)果可以得到很大改善，優(yōu)于Delta方法。

3 結(jié)果與分析

3.1 平穩(wěn)序列GEV模型

首先進行序列平穩(wěn)性特征的檢驗，主要考察其是否存在趨勢或突變。計算巴音布魯克站最大1d降水量序列的1至10階自相關(guān)系數(shù)，由于各序列的1階自相關(guān)系數(shù)均不顯著，其他各階自相關(guān)系數(shù)中絕大部分也不顯著，因此在趨勢檢驗中沒有進行自相關(guān)系數(shù)的校正。對該站AMS1序列進行MK(Mann-Kendall)和線性趨勢檢驗，發(fā)現(xiàn)其極值降水序列存在較顯著的上升趨勢，其MK趨勢Z值和線性趨勢t值分別為2.78和1.26，達到0.05顯著性水平。

根據(jù)趨勢和突變檢驗的結(jié)果，巴音布魯克站極值降水序列存在非平穩(wěn)性，適合于采用非平穩(wěn)GEV模型進行模擬。對該站AMS1序列進行平穩(wěn)GEV模擬，模擬診斷結(jié)果見圖1。

圖1 巴音布魯克站AMS1序列GEV模擬結(jié)果診斷

圖1(a)為Q-Q圖，縱坐標與橫坐標分別為經(jīng)驗和模型擬合分位數(shù)。在模擬結(jié)果非常好的情況下，數(shù)據(jù)將沿著對角線排列?？梢?，絕大部分點據(jù)位于對角線上，僅僅在高分位數(shù)部分有個別奇異點。圖1(b)與Q-Q圖很相似，但它是從擬合分布中隨機取一個樣本，將其分位數(shù)和實際數(shù)據(jù)的經(jīng)驗分位數(shù)繪制而得。其中1∶1線是當(dāng)x=y時的對角線，回歸直線是數(shù)據(jù)點的擬合直線，可見點據(jù)擬合線與1∶1線非常接近。圖中的置信限提供了關(guān)于其尾部不確定性的直觀了解，可見隨著分位數(shù)增加，不確定性顯著增大。分析這種不確定性可能與降水序列在1975年和2005年存在兩個超大值有關(guān)。概率密度圖[圖1(c)]中經(jīng)驗和模擬的密度曲線較一致，也說明參數(shù)估計結(jié)果是合理的。從重現(xiàn)水平圖[圖1(d)]來看，所有的點據(jù)都位于95%置信區(qū)間內(nèi)，且大多集中在置信區(qū)間中央的凹曲線附近，表明所構(gòu)建模型是穩(wěn)定可靠的；其曲線下凹也說明形狀參數(shù)大于0，服從Fréchet分布。

3.2 兩種方法計算的置信限

分析基于Delta方法和輪廓似然函數(shù)方法分別得到的置信區(qū)間，結(jié)果見圖2。在重現(xiàn)期低于5年左右的情況下，Delta方法得到的置信區(qū)間在置信下限和置信上限均大于輪廓似然函數(shù)方法；重現(xiàn)期大于5年左右的情況下則相反，Delta方法的置信下限和置信上限均低于輪廓似然函數(shù)方法。在重現(xiàn)期大于100年之后，兩種方法得到的置信限均較寬，反映對遠超過序列長度的重現(xiàn)期進行推斷時所產(chǎn)生的固有不確定性很大。但在重現(xiàn)期較大的情況下，輪廓似然函數(shù)方法給出的結(jié)果更加準確，因為該方法考慮了參數(shù)分布的偏態(tài)特征。此外，由圖2可見，在重現(xiàn)期較小的情況下，Delta方法的結(jié)果與輪廓似然函數(shù)方法比較接近，也可以得到較合理的置信區(qū)間。

圖2 AMS1序列重現(xiàn)水平

3.3 不同重現(xiàn)期下的極值降水

在非平穩(wěn)GEV模型中，其位置參數(shù)μ是隨時間變化的，因此計算得到的極值降水重現(xiàn)水平(極值降水設(shè)計值)也將隨著時間的推移而變化，稱為等效重現(xiàn)水平(effective return level)，它是一組數(shù)值而不是唯一確定的值，因此等效重現(xiàn)水平擴展了重現(xiàn)水平的概念。由圖3可見，基于非平穩(wěn)GEV模型計算得到的極值降水設(shè)計值，是隨著時間的推移而不斷上升的一組值，在1958年的100年一遇設(shè)計值到2010年下降為接近50年一遇，預(yù)示著未來發(fā)生極值降水和洪水的風(fēng)險會加大。

圖3 基于非平穩(wěn)GEV模型的AMS1重現(xiàn)水平序列

等效重現(xiàn)水平的概念，與時變矩方法計算的設(shè)計值很類似。在當(dāng)前水文時間序列常常表現(xiàn)出明顯非平穩(wěn)特征的情況下，傳統(tǒng)頻率分析方法得到的設(shè)計值結(jié)果常常難以滿足需求。例如，當(dāng)徑流序列的均值存在上升趨勢時，其洪水設(shè)計值將被低估，從而造成水工建筑物的失事風(fēng)險加大；反之，若存在下降趨勢，設(shè)計值將被高估，從而增加建筑物的造價。

4 結(jié) 論

a.巴音布魯克站年最大1d降水量序列(AMS1)存在顯著上升趨勢，該站以時間作為位置參數(shù)的協(xié)變量進行非平穩(wěn)GEV模擬，似然比檢驗結(jié)果表明其模擬效果優(yōu)于平穩(wěn)GEV模型。

b.與Delta方法相比較，在較高重現(xiàn)期下輪廓似然函數(shù)方法給出的置信區(qū)間更準確；同時，重現(xiàn)水平的輪廓對數(shù)似然函數(shù)曲線在較高重現(xiàn)期的情況下呈較顯著的不對稱性。

c.采用等效重現(xiàn)水平的概念，基于非平穩(wěn)GEV模型計算得到了不同重現(xiàn)期之下的巴音布魯克站極值降水的設(shè)計值；它們分別是隨著時間推移而不斷上升的一組值，在1958年的100年一遇設(shè)計值到2010年下降為接近50年一遇。該結(jié)果對于評估當(dāng)前氣候變化下的區(qū)域極值降水與洪災(zāi)風(fēng)險具有重要的啟示。

[1] 顧西輝,張強,王宗志.1951—2010年珠江流域洪水極值序列平穩(wěn)性特征研究[J].自然資源學(xué)報,2015,30(5):824-835.

[2] MILLY P C D,BETANCOURT J.Stationarity is dead:Whither water management[J].Science,2008(319):573-574.

[3] 葉長青,陳曉宏,張家鳴.具有趨勢變異的非一致性東江流域洪水序列頻率計算研究[J].自然資源學(xué)報,2013,28(12):2105-2116.

[4] 謝平,陳廣才,夏軍.變化環(huán)境下非一致性年徑流序列的水文頻率計算原理[J].武漢大學(xué)學(xué)報(工學(xué)版),2005,38(6):6-15.

[5] 宋松柏,李揚.具有跳躍變異的非一致分布水文序列頻率計算方法[J].水利學(xué)報,2012,43(6):734-739,748.

Simulation of upstream precipitation of Kaidu River based on GEV model

WANG Zhicheng

(XinjiangTarimRiverBasinAdministration,Korla841000,China)

GEV model is adopted for simulation and analysis of extreme precipitation on the basis of upstream precipitation data in Kaidu River Basin. Results show that the shape parameter of AMS1 sequence is greater than 0, which is consistent with Frechet distribution. The confidence interval estimated by the profile likelihood method is more accurate than Delta method under higher return period. The profile likelihood function curve of the return level shows more prominent asymmetry under higher return period. The extreme precipitation design value is obtained on the basis of non-stationary GEV model. The once-in-a-century design value in 1958 is reduced to once-every-50 years in 2010. It indicates that the extreme precipitation and flooding risks are increased in the future.

GEV model; extreme precipitation; Kaidu River upstream area; non-stationary

10.16616/j.cnki.10-46/TV.2017.03.016

TV125

A

2096-0131(2017)03- 0052- 04