雅魯藏布江流域極端降水模擬及預(yù)估

高佳佳， 杜 軍

（1.中國(guó)氣象局成都高原氣象研究所拉薩分部，西藏 拉薩 850000； 2.西藏高原大氣環(huán)境科學(xué)研究所，西藏 拉薩 850000；3.西藏高原大氣環(huán)境研究重點(diǎn)實(shí)驗(yàn)室，西藏 拉薩 850000）

0 引言

極端降水是全球最受關(guān)注、影響最大的自然災(zāi)害之一，是短期氣候預(yù)測(cè)研究的重點(diǎn)［1］。IPCC［2］曾指出，隨著全球氣候變暖，內(nèi)陸地區(qū)的極端降水事件頻率呈現(xiàn)出增加趨勢(shì)。近年來，區(qū)域性洪水、干旱、高溫、雨雪冰凍等極端事件頻發(fā)，尤其是20世紀(jì)80 年代以后，頻繁的極端事件給生態(tài)環(huán)境、經(jīng)濟(jì)發(fā)展和人民生活造成了嚴(yán)重影響［3-4］。因此，研究區(qū)域性極端降水事件對(duì)科學(xué)認(rèn)識(shí)氣候變化背景下水循環(huán)的時(shí)空演變，把握氣候異常對(duì)極端降水的影響規(guī)律，評(píng)估水資源管理及區(qū)域水資源安全具有十分重要的理論和實(shí)踐意義。

極端氣候事件歸根到底是氣候極值問題，氣候極值是極端事件產(chǎn)生的必要條件，極端事件發(fā)生發(fā)展的預(yù)測(cè)首先要考慮極值的分布規(guī)律。國(guó)內(nèi)外部分學(xué)者研究表明，中國(guó)西部的極端降水天數(shù)呈增加趨勢(shì)，尤其是西北地區(qū)［5-8］。Fischer 等［9］分析了珠江的極端降水分布特征，并估算了極端降水指數(shù)。Hong 等［10］認(rèn)為海河流域的極端降水主要發(fā)生在38°N，大部分站點(diǎn)的降水呈現(xiàn)出減少趨勢(shì)。程炳巖等［11］、江志紅等［12］研究認(rèn)為廣義帕累托分布函數(shù)（generalized Pareto distribution，GPD）在重慶、中國(guó)東部的日降水模擬中具有更高的擬合度；李占玲等［13］基于GPD 函數(shù)分析了黑河流域的極端降水頻率特征，得出該流域20 世紀(jì)60 年代發(fā)生的極端降水次數(shù)最多，90 年代以后次數(shù)較少。Eylon 等［14］運(yùn)用極值理論分析了巴拿馬運(yùn)河的極端降水分布特征，并估算了極端降水的重現(xiàn)期及相應(yīng)的置信區(qū)間。劉彩紅等［15］運(yùn)用CMIP5 模式指出，青海高原的降水有極端化的趨勢(shì)，極端降水頻次增加，強(qiáng)度增大。韓國(guó)軍［16］、游慶龍等［17］運(yùn)用統(tǒng)計(jì)方法計(jì)算出青藏高原極端降水大部分呈增加趨勢(shì)，且逐年平均降水強(qiáng)度和逐年連續(xù)降水天數(shù)均有所增加，90 年代以來增加明顯。

雅魯藏布江發(fā)源并流經(jīng)西藏高原，地理位置特殊，是世界上海拔最高的大河之一，平均海拔4 000 m 以上，是全球氣候主要變化區(qū)與敏感區(qū)。流域沿岸為西藏主要農(nóng)、牧業(yè)生產(chǎn)區(qū)，其洪澇和干旱的頻繁發(fā)生導(dǎo)致了水資源分布不均，進(jìn)而影響了流域的用水矛盾和生態(tài)環(huán)境的惡化，而極端事件的發(fā)生是對(duì)區(qū)域氣候、環(huán)境變化的重要響應(yīng)。目前對(duì)氣候極值進(jìn)行定量評(píng)估的方法以氣候動(dòng)力模式為主，從概率論角度對(duì)極端氣候事件及可預(yù)測(cè)性研究并不多，尤其是預(yù)測(cè)方法。IPCC報(bào)告中特別強(qiáng)調(diào)統(tǒng)計(jì)方法對(duì)極端氣候事件的重要性［2］，因此，本文基于廣義帕累托方法（GPD）的分布參數(shù)模型，針對(duì)超出閾值的數(shù)據(jù)作為樣本數(shù)據(jù)來建模，從氣候極值的分布規(guī)律出發(fā)，揭示極端降水的發(fā)生發(fā)展規(guī)律，探索極端降水的可預(yù)測(cè)性，從而更好地預(yù)估極端事件，為提高防災(zāi)減災(zāi)能力提供科學(xué)依據(jù)。

1 研究區(qū)與數(shù)據(jù)

1.1 研究區(qū)域

雅魯藏布江（簡(jiǎn)稱雅江）全部在中國(guó)境內(nèi)，橫貫西藏高原南部，干流全長(zhǎng)約2.1×103km，流域面積2.4×105km2，雅魯藏布江干流河谷沿東西向的斷裂帶發(fā)育，流域呈東西向的狹長(zhǎng)帶，支流多而短小，較大支流有拉薩河、雅魯藏布江帕隆藏布、易貢藏布、拉喀藏布、尼澤曲、年楚河等。干流在仲巴縣里孜以上為上游，是高寒河谷地帶。里孜到米林縣派區(qū)為中游，支流眾多，流量增大，河谷展寬，氣候溫和，水利條件較好，是西藏農(nóng)業(yè)最發(fā)達(dá)的地區(qū)。下游位于林芝一帶。截止國(guó)境線，年徑流總量為1.1×1010m3，洪水由強(qiáng)降水形成，持續(xù)時(shí)間較長(zhǎng)。流域東部地區(qū)降水量充足，年平均降水量超過900 mm，達(dá)到半濕潤(rùn)地區(qū)水平，西部地區(qū)降水量少，年均降水量不足100 mm，為干旱地區(qū)，整個(gè)流域的降水量從東至西呈現(xiàn)遞減的趨勢(shì)。

圖1 研究區(qū)域位置示意Fig.1 Location of study area

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

選取流域內(nèi)9 個(gè)氣象站近50 年（1959—2017年）5—9 月逐日降水資料作為研究對(duì)象。部分站點(diǎn)開始于1978年和1979年。

2 研究方法

2.1 GPD分布函數(shù)和密度函數(shù)

GPD 可以直接利用歷年的原始數(shù)據(jù)，人為設(shè)置閾值，在設(shè)置好閾值后，以此為標(biāo)準(zhǔn)來抽取每一年超過此閾值的極大或極小值，即“超門限峰值POT”（peaks over threshold），可以提高估算精度［18］。具體為：

GPD的分布函數(shù)為：

其相應(yīng)的密度函數(shù)為：

式中：ξ為門限值；a為尺度參數(shù)；k為形狀參數(shù)，自變量x的取值，取決于k的值，當(dāng)k< 0時(shí)，ξ 0 時(shí)，GPD 為指數(shù)分布；當(dāng)k=1 時(shí)，GPD 在區(qū)間[ξ，ξ+a]上均勻分布。

表1 雅魯藏布江流域內(nèi)各站點(diǎn)資料長(zhǎng)度Table 1 Basic information of each stations in Yarlung Zangbo River basin

2.1.1 閾值的選取

GPD 模型的核心在于確定閾值，它是正確估計(jì)參數(shù)的前提。如果閾值選取的過高，會(huì)使得超額數(shù)據(jù)量太少，導(dǎo)致估計(jì)出來的參數(shù)方差很大；如果閾值選取的過低，則不能保證函數(shù)的收斂性，所估計(jì)參數(shù)有較大偏差。本文主要使用Hill 圖估計(jì)、百分位法和年交叉率法來確定各站點(diǎn)日降水量閾值。

Hill 圖法是基于Hill 估計(jì)量的一種閾值圖形法［19］，由點(diǎn)（k，1/H（k，n））構(gòu)成的曲線，通過觀察圖中尾部指數(shù)穩(wěn)定的區(qū)域來選擇閾值。其定義為：

百分位法是將該站點(diǎn)的日降水量從小到大排序，并計(jì)算相應(yīng)的累計(jì)百分位，某一百分位所對(duì)應(yīng)數(shù)據(jù)的值就稱為這一百分位的百分位數(shù)，文中分別計(jì)算了90、95、97、98、99百分位。

年交叉率法是假定超閾值降水極值出現(xiàn)次數(shù)服從泊松分布，以一年為時(shí)段所量度，極值超過閾值的次數(shù)，λ=n/T，其中，λ為年交叉率，T為資料總年數(shù)。

2.1.2 參數(shù)估計(jì)

對(duì)GPD 分布進(jìn)行參數(shù)估計(jì)的方法有很多，本文主要通過極大似然估計(jì)法對(duì)參數(shù)進(jìn)行估計(jì)。該方法具有很強(qiáng)的靈活性，可以適應(yīng)不同模型的需求，且統(tǒng)計(jì)特性良好，能夠綜合各種有關(guān)信息到統(tǒng)計(jì)推斷中去。模型估計(jì)時(shí)，樣本變異可能會(huì)導(dǎo)致模型參數(shù)的不確定性，由于極大似然方法具有漸進(jìn)正態(tài)性，容易給出估計(jì)值及其標(biāo)準(zhǔn)誤差（標(biāo)準(zhǔn)誤差是參數(shù)不確定性或變異性的度量之一）。該方法唯一的缺點(diǎn)是計(jì)算時(shí)迭代繁瑣。具體方法見文獻(xiàn)［20］。

2.1.3 重現(xiàn)期計(jì)算

某一指定重現(xiàn)期T時(shí)間的降水量分位計(jì)算公式：

2.2 獨(dú)立性和穩(wěn)定性檢驗(yàn)

GPD 擬合需要超閾值數(shù)據(jù)序列滿足平穩(wěn)性的條件，因此，擬合之前需對(duì)超閾值序列進(jìn)行穩(wěn)定性檢驗(yàn)。文章中使用Mann-Kendall（M-K 檢驗(yàn)）對(duì)序列的變化趨勢(shì)和突變點(diǎn)進(jìn)行檢驗(yàn)。M-K 檢驗(yàn)［21］是氣象學(xué)、氣候?qū)W中經(jīng)常用來進(jìn)行突變檢驗(yàn)的一種非參數(shù)檢驗(yàn)方法，它不要求樣本符合一定的分布。即給定顯著性水平α=0.05，則統(tǒng)計(jì)量的臨界值為±1.96。統(tǒng)計(jì)量大于0，表示序列呈上升趨勢(shì)；反之，表明呈下降趨勢(shì)，大于或小于±1.96，表示上升或下降趨勢(shì)明顯。該方法能有效區(qū)分某一自然過程是處于自然波動(dòng)還是存在確定的變化趨勢(shì)，常用于氣候變化影響下的降水、干旱頻次趨勢(shì)檢測(cè)。

2.3 擬合度檢驗(yàn)

通過對(duì)超閾值序列進(jìn)行平穩(wěn)性檢驗(yàn)后，即可對(duì)GPD 進(jìn)行擬合，擬合結(jié)果需通過Kolmogorov-Simimov（K-S 檢驗(yàn)）。一般在K-S 檢驗(yàn)中，先計(jì)算需要做比較的兩組觀察數(shù)據(jù)的累積分布函數(shù)，然后求這兩個(gè)累積分布函數(shù)的差的絕對(duì)值中的最大值D。最后通過查表以確定D值是否落在所要求對(duì)應(yīng)的置信區(qū)間內(nèi)。若D值落在了對(duì)應(yīng)的置信區(qū)間內(nèi)，說明被檢測(cè)的數(shù)據(jù)滿足要求。

3 結(jié)果與分析

3.1 閾值的選取

通過Hill 圖尾部特征的穩(wěn)定性來選取閾值（圖2），由于Hill圖具有較強(qiáng)的主觀性，從而會(huì)導(dǎo)致選取的閾值不同。以中游地區(qū)的拉孜、拉薩和墨竹工卡站為例，拉孜站點(diǎn)的Hill 圖尾部指數(shù)趨于穩(wěn)定大致位于65～70 附近，對(duì)應(yīng)的降水量為7.9～15.4 mm；拉薩和墨竹工卡站的尾部指數(shù)均位于160左右時(shí)趨于穩(wěn)定，對(duì)應(yīng)的降水量分別為8.7～18.5 mm 和9.9～20.3 mm。同樣的，下游地區(qū)的林芝站在80～120位置時(shí)趨于穩(wěn)定，對(duì)應(yīng)的降水量為14.4～17.6 mm。為更好的確定閾值，在此基礎(chǔ)上，我們結(jié)合百分位法和年交叉率法進(jìn)行閾值確定（表2）。

圖2 各站點(diǎn)Hill圖指數(shù)分布Fig.2 The distribution of index of Hill plot for each stations

表2 各站點(diǎn)的百分位閾值（mm）和年交叉率Table 2 Precipitation threshold selection in GPD and the average annual occurrence number for each stations

與Hill 圖相比，百分位法可以更精確的確定閾值。以拉孜站為例，Hill 圖顯示的閾值為7.9～15.4 mm，而對(duì)應(yīng)的是該站點(diǎn)的95 百分位；同樣拉薩站、墨竹工卡站的Hill 圖確定的閾值均位于93 百分位和95 百分位?？傮w而言，Hill 圖確定的閾值要小于百分位確定的閾值。根據(jù)前人研究結(jié)果得出，年交叉率為1～2 時(shí)對(duì)應(yīng)的閾值可作為GPD 擬合分析時(shí)的參考閾值。結(jié)合表1 得出，當(dāng)各站點(diǎn)日降水量達(dá)到99 百分位時(shí)，年交叉率均穩(wěn)定在1.5 附近，因此我們確定99 百分位時(shí)的閾值為最佳閾值。

3.2 超閾值序列的平穩(wěn)性檢驗(yàn)

根據(jù)閾值的平穩(wěn)性和穩(wěn)定性要求，我們用M-K法在顯著性水平α=0.05 條件下檢測(cè)超閾值降水序列的變化趨勢(shì)和突變。結(jié)果顯示，下游地區(qū)的林芝和米林站的統(tǒng)計(jì)量略大于1.96，分別為2.3 和2.8，其余大部分站點(diǎn)均位于臨界區(qū)域內(nèi)，通過了顯著性水平檢驗(yàn)。通過時(shí)間序列曲線（UFk曲線）可以看出（圖3），米林、墨竹工卡、南木林站呈顯著增長(zhǎng)趨勢(shì)，說明三個(gè)站點(diǎn)的超閾值序列的日降水量呈逐漸增加趨勢(shì)；其余站點(diǎn)呈下降趨勢(shì)（此部分只給出拉孜、墨竹工卡、日喀則、米林站的趨勢(shì)變化圖，其余圖表省略）。各站點(diǎn)超閾值序列沒有明顯的突變性。

圖3 流域內(nèi)四個(gè)站點(diǎn)的時(shí)間序列統(tǒng)計(jì)量變化圖Fig.3 Time series statistics of four sites in the basin

3.3 模型擬合效果檢驗(yàn)

通過對(duì)流域內(nèi)各站點(diǎn)進(jìn)行GPD 擬合（圖4），并使用K-S檢驗(yàn)看其是否符合已知理論分布函數(shù)。大部分站點(diǎn)的統(tǒng)計(jì)值小于0.01顯著性水平，接受原假設(shè)。說明雅江流域各站點(diǎn)之間雖存在差異，但由GPD 擬合曲線可知，理論頻數(shù)和實(shí)測(cè)頻數(shù)基本相符。另外，雖然西藏地區(qū)大部分臺(tái)站的觀測(cè)記錄起始年代不一致，如拉薩站、拉孜站，資料開始時(shí)間分別為1955 年、1977 年，但從圖4 中可看出，數(shù)據(jù)滿足方程需求，擬合結(jié)果表明資料長(zhǎng)度并不影響降水極值的統(tǒng)計(jì)推斷，且資料年限越長(zhǎng)擬合結(jié)果越好。

圖4 各站點(diǎn)累積頻率和實(shí)際頻率分布曲線對(duì)比Fig.4 The distribution of cumulative frequency and empirical frequency over the observation stations

3.4 雅江流域極端降水的分布特征

從各站點(diǎn)的極端降水閾值分布來看，林芝地區(qū)的閾值最大，閾值最小的站點(diǎn)為拉孜站。這與我們的觀測(cè)事實(shí)一致，林芝地區(qū)日降水量大，連續(xù)降雨日數(shù)長(zhǎng)，得到的閾值就大。

尺度參數(shù)主要是描述極值分布的變率，尺度參數(shù)越大，極值波動(dòng)范圍越大，表明打破極端降水的記錄值也越大。整體而言，雅江流域的尺度參數(shù)由下游向上游是逐漸減小的，平均值為5.95。由表3可知，下游地區(qū)的尺度參數(shù)最大，約為7.00，表明這一區(qū)域的極端降水變化幅度很大。從氣候背景來看，該地區(qū)位于高原季風(fēng)區(qū)，受印緬槽和西風(fēng)帶影響，季節(jié)性降水較大［22-23］，5—9 月的降水總量可達(dá)600 mm，是西藏地區(qū)夏季降水量最大的區(qū)域，因此可能出現(xiàn)的破極端降水記錄值要高于其他地區(qū)。尺度參數(shù)最小的區(qū)域位于流域中上游地區(qū)，受地理位置和大氣環(huán)流影響，雨期短，降水量少，且連續(xù)降水日數(shù)也少，降水極值的范圍比較小，區(qū)域打破極端降水的記錄值要比下游地區(qū)低。這與前人研究結(jié)果一致［24］，我國(guó)的干旱地區(qū)大部分位于非季風(fēng)區(qū)，降水極值范圍較小，破紀(jì)錄的可能性較季風(fēng)區(qū)小。

表3 流域內(nèi)各站點(diǎn)的GPD模型參數(shù)估計(jì)及檢驗(yàn)Table 3 Estimation and validation of parameters in GPD model

形狀參數(shù)作為模型的第二個(gè)重要參數(shù)，不同的形狀有不同的尾部分布特征，它表示該區(qū)域極端降水的破紀(jì)錄率。由表3 看出，形狀參數(shù)正值區(qū)主要位于拉孜地區(qū)，這些地區(qū)發(fā)生破紀(jì)錄降水事件的可能性比其他地區(qū)大。正是由于該區(qū)域降水日數(shù)少，所以一旦有降水過程，就可能會(huì)打破降水極值。而下游地區(qū)的形狀參數(shù)為負(fù)值，說明這些區(qū)域的降水發(fā)生破紀(jì)錄的概率偏小。因?yàn)橄募荆涤暝筋l繁的區(qū)域，極值變率大，較均值離散程度大，則破紀(jì)錄的概率較小。有降水的地區(qū)，不是有較大的形狀參數(shù)就是有較大的尺度參數(shù)，不可能兩個(gè)參數(shù)都大。形狀參數(shù)較大的地區(qū)，稱為“形狀參數(shù)主導(dǎo)區(qū)”，該區(qū)域多持續(xù)降水，極端氣候事件較少；尺度參數(shù)較大的區(qū)域稱為“尺度參數(shù)主導(dǎo)區(qū)”，該區(qū)域降水較少，且多變，極端氣候事件較多［23］。

3.5 雅江流域極值重現(xiàn)期預(yù)估

極值模型建立最重要的目的之一就是預(yù)測(cè)極端事件的重現(xiàn)期或重現(xiàn)水平。如表4 所示，根據(jù)重現(xiàn)期公式可以得出，從5 年一遇和10 年一遇的極端降水值來看，雅江流域除拉孜站外，其他地區(qū)降雨極值均超過30 mm，日喀則地區(qū)的降水極值達(dá)50 mm，其中拉薩、澤當(dāng)、墨竹工卡和南木林站5 年一遇和10年一遇的極端降水量分別在40 mm 左右。相關(guān)研究指出［26-27］，當(dāng)溫度上升為2 ℃時(shí)，青藏高原的強(qiáng)降水距平百分率平均增多44.5%～59.5%，大值區(qū)出現(xiàn)在山南附近，這與我們的研究結(jié)論相符合。各站點(diǎn)在15年一遇的極端降水值之后，極值水平的增長(zhǎng)變得非常緩慢，其中林芝地區(qū)的增長(zhǎng)最緩慢，以0.32 mm·a-1的速率增長(zhǎng)；日喀則地區(qū)的降雨極值增長(zhǎng)率最快，約0.72 mm·a-1。西藏地區(qū)極端降水頻率一般為每年4.3 次，強(qiáng)度在20 mm·d-1以上，林芝地區(qū)為極端降雨的高值區(qū)，且沿雅江一線極端降水的頻次呈增加趨勢(shì)［28-29］。

表4 雅江流域各站點(diǎn)日降水量極大值重現(xiàn)水平（單位：mm）Table 4 The maximum daily precipitation of flood season reappeared in Yarlung Zangbo River（unit：mm）

為檢驗(yàn)各站點(diǎn)不同重現(xiàn)期水平的合理性，將其帶入各站點(diǎn)逐日降水序列中進(jìn)行驗(yàn)證。以拉薩和拉孜站為例，拉薩站5 年重現(xiàn)期水平值為39 mm，在1967—2017 年所有逐日降水中，共有10 次大于39 mm的降水過程，平均每5年一次。拉孜站5年重現(xiàn)期水平值為35.1 mm，在1977—2017 年期間，共有9 次大于35.1 mm 的降水過程，平均每4.4 年一次。由GPD 擬合計(jì)算出的極端降水重現(xiàn)期水平基本符合實(shí)際。值得注意的是這里的“重現(xiàn)期”并不意味著經(jīng)過T年之后一定會(huì)出現(xiàn)的“周期”，它是概率意義上的“統(tǒng)計(jì)周期”。

4 結(jié)論

雅魯藏布江作為高原河流，由于強(qiáng)降水的時(shí)空分布不均而引起洪澇和干旱，并對(duì)流域內(nèi)的農(nóng)牧業(yè)產(chǎn)生重要威脅。因此在氣候日益增暖趨勢(shì)下，評(píng)估極端降水規(guī)律及其發(fā)生概率十分必要。本文通過引進(jìn)GPD 概率分布模型，對(duì)西藏地區(qū)汛期強(qiáng)降水規(guī)律進(jìn)行模擬。結(jié)果表明：

（1）通過Hill 圖法選取的流域內(nèi)各站點(diǎn)的閾值序列小于百分位法選取的閾值序列，綜合考慮Hill圖法、百分位法及年交叉率法最終確定99百分位時(shí)的閾值為最佳閾值。

（2）各站點(diǎn)閾值序列在M-K 顯著性水平檢驗(yàn)下，無明顯突變。擬合效果通過K-S檢驗(yàn)，各站點(diǎn)擬合的理論頻數(shù)和實(shí)測(cè)頻數(shù)基本相符，且資料長(zhǎng)度并不影響降水極值的統(tǒng)計(jì)推斷。

（3）通過分析流域內(nèi)各站點(diǎn)擬合的極端降水特征可知，尺度參數(shù)的大值區(qū)位于流域下游，即林芝、米林地區(qū)，表明該地區(qū)的極值波動(dòng)大；相反地，小值區(qū)位于流域中上游的拉孜站附近，表明極值波動(dòng)小。形狀參數(shù)正值區(qū)位于流域中上游地區(qū)，說明發(fā)生破紀(jì)錄的降水事件概率較大，擬合結(jié)果與實(shí)際觀測(cè)一致。

（4）從5 年一遇和10 年一遇的極值水平看，雅江流域除拉孜站外，其他地區(qū)降雨極值均超過30 mm，日喀則地區(qū)的降水極值達(dá)50 mm；各站點(diǎn)在15 年一遇的極端降水值之后，極值水平的增長(zhǎng)變得非常緩慢，由GPD 擬合計(jì)算出的降水極值具有一定的合理性。