1956—2018年青海省黃河流域降水變化特征分析及預(yù)測?

冷 雪， 龐 亮??， 姜欣彤， 董 勝

(1. 中國海洋大學(xué)工程學(xué)院， 山東 青島 266100； 2. 中水珠江規(guī)劃勘測設(shè)計(jì)有限公司， 廣東 廣州 510610)

黃河流域水資源狀況關(guān)系著社會發(fā)展及生態(tài)安全。降水是黃河流域水資源的主要來源，其時空分布與變化直接影響水資源的合理開發(fā)與利用，也是進(jìn)一步研究水循環(huán)的基礎(chǔ)[1-2]。研究降水量的分布特征及趨勢預(yù)測，可為抗旱防澇提供依據(jù)，對了解降水趨勢變化規(guī)律、保證農(nóng)業(yè)、水利部門決策以及保障當(dāng)?shù)厣鐣?jīng)濟(jì)發(fā)展具有重要意義。

青海省作為黃河流域的發(fā)源省份，降水變化一直備受關(guān)注。戴虎德等[3]分析了黃河源頭地區(qū)近50年降水及氣溫變化特點(diǎn)，并指出黃河源頭年和四季平均氣溫呈上升趨勢，地表蒸發(fā)量增大，而降水量無明顯變化；白淑英等[4]研究了青海省近60年降水?dāng)?shù)據(jù)的空間分布特征及變化規(guī)律，指出降水分布呈現(xiàn)東多西少、南多北少的格局；王文等[5]利用黃河中上游24站夏季降水資料，提出基于EEMD的統(tǒng)計(jì)預(yù)報模型，提高了降水預(yù)報的精度。徐慧等[6]分析了青海省月降水異質(zhì)性指數(shù)和降水集中性指數(shù)的區(qū)域分布及變化趨勢，并指出季節(jié)性降水的增加有助于改善月降水分配的均勻性；王明常等[7]的研究發(fā)現(xiàn)，青海省極端降水事件及極端降水指數(shù)均呈上升趨勢，且極端降水事件空間差異顯著；劉曉瓊等[8]的研究發(fā)現(xiàn)，在過去55年三江源地區(qū)總體降水量變化更趨向于青藏高原的降水變化，呈現(xiàn)弱增趨勢。

青海省降水分析現(xiàn)有研究多采用降水距平分析、小波分析、Arcgis空間插值等方法進(jìn)行降水空間分析，而對于時間尺度上的降水預(yù)測研究較少；在黃河流域降水研究方面，以上、中、下游分區(qū)研究居多，在單個省份的研究較少。為探討青海省降水預(yù)測模型合理性以及豐富黃河流域研究內(nèi)容，本文基于M-K檢驗(yàn)、Pettitt檢驗(yàn)、Morlet小波分析以及ARIMA模型等方法，采用青海省黃河流域31個氣象站1956—2018年月值降水資料，重點(diǎn)探討年、季降水變化特性以及降水預(yù)測變化的合理性，從而對當(dāng)?shù)厮Y源規(guī)劃和管理做出指導(dǎo)。

1 區(qū)域概況與數(shù)據(jù)來源

1.1 區(qū)域概況

青海省黃河流域位于青海省東部，地處95°39′E—103°06′E，32°43′N—39°08′N，包括龍羊峽以上的黃河源頭和龍羊峽至寺溝峽的高山峽谷，地形復(fù)雜，海拔高度呈現(xiàn)出由西南向東北遞減的特點(diǎn)，大部分海拔高程在3 000 m以上，平均氣溫為-6～9.2 ℃。氣候?qū)俑咴箨懶愿吆畾夂?，寒冷干燥，雨熱同期。多年平均降水量?30～650 mm之間，區(qū)域地貌復(fù)雜多樣。

1.2 數(shù)據(jù)來源與處理

本研究選取1956—2018年青海省黃河流域內(nèi)31個氣象站點(diǎn)的逐日降水?dāng)?shù)據(jù)。數(shù)據(jù)由中國氣象科學(xué)數(shù)據(jù)共享服務(wù)網(wǎng)平臺提供(http://data.cma.cn/)，質(zhì)量可靠，站點(diǎn)分布見圖1。首先利用Excel軟件對青海省黃河流域31個氣象站點(diǎn)近63 年逐日降水?dāng)?shù)據(jù)進(jìn)行初步處理，對于霜、雪、霧以及微量降水進(jìn)行特殊處理、對于少量不觀測或者無數(shù)據(jù)情況采用鄰近日數(shù)據(jù)線性回歸插值處理，從而計(jì)算各站點(diǎn)逐月降水?dāng)?shù)據(jù)；其次利用Arcgis軟件采用修正反距離加權(quán)平方法[9]對31個氣象站點(diǎn)月值降水?dāng)?shù)據(jù)進(jìn)行空間差值處理，得到青海省黃河流域區(qū)域近63年逐月降水?dāng)?shù)據(jù)。

圖1 青海省黃河流域氣象站點(diǎn)分布圖Fig.1 Distribution map of meteorological stations in the Yellow River Basin of Qinghai Province

2 研究方法

2.1 Mann-Kendall檢驗(yàn)法

Mann-Kendall檢驗(yàn)法(簡稱M-K檢驗(yàn)法)是世界氣象組織(World Meteorological Organization, WMO)推薦并廣泛使用的非參數(shù)檢驗(yàn)方法，目前被廣泛應(yīng)用于水文序列趨勢變化及突變點(diǎn)檢驗(yàn)[10-13]。在趨勢變化方面，M-K檢驗(yàn)?zāi)軌蛱蕹贁?shù)異常值，且數(shù)據(jù)不需要服從一定的分布，適用于一般數(shù)據(jù)序列趨勢檢驗(yàn)。在突變分析方面，M-K檢驗(yàn)?zāi)軌虿檎覕?shù)據(jù)異常點(diǎn)，但有時會產(chǎn)生虛假突變點(diǎn)的情況，從而難以判斷真正的突變位置[14]。本文采用M-K檢驗(yàn)和Pettitt突變檢驗(yàn)[15]兩種方法來確定符合實(shí)際的突變時間及突變區(qū)域。

2.2 小波分析方法

連續(xù)小波變換是應(yīng)用最為廣泛的變換方法之一，適合分析各評價因子的周期性變化，對原始數(shù)據(jù)要求較低，適用性好[16-18]。目前小波分析技術(shù)在水文學(xué)中的研究主要集中在水文序列多時間尺度分析和序列模擬預(yù)測兩大方面[19-21]。本研究采用Morlet小波為母小波，進(jìn)行多時間尺度的周期變換。

2.3 ARIMA模型

差分自回歸移動平均(Auto Regressive Integrated Moving Average , ARIMA)模型是由Box和Jenkins在20世紀(jì)70年代初提出的時間序列檢測方法，簡稱ARIMA模型。ARIMA模型由3個重要參數(shù)(p,d,q)決定:p為自回歸系數(shù)，表示序列滯后p階；d為時間序列變?yōu)槠椒€(wěn)時間序列需要進(jìn)行的差分次數(shù)；q為滑動平均系數(shù)，表示誤差項(xiàng)滯后q階。ARIMA(p,d,q)模型已較為成熟，通過選擇模型的參數(shù)和辨識模型的系數(shù)即可實(shí)現(xiàn)對時間序列的擬合，模型評價機(jī)制完善，預(yù)測效果較好。此外，ARIMA(p,d,q)模型具有建模簡單，擬合效果好，短期預(yù)測誤差小等優(yōu)點(diǎn)，在水文時間序列的預(yù)測領(lǐng)域有廣泛的應(yīng)用[22-23]。

3 青海省黃河流域降水變化特性分析及預(yù)測

3.1 降水年內(nèi)分配及年代變化特征分析

將青海省黃河流域1956—2018年降水量按照月份為單位進(jìn)行累加統(tǒng)計(jì)可得到表1，季節(jié)劃分按12月至次年2月為冬季，3—5月為春季，6—8月為夏季，9—11月為秋季。統(tǒng)計(jì)顯示1956—2018年累計(jì)降水量達(dá)28 362 mm，其中，春季降水量占全年累計(jì)降水量的18.4%，夏季降水量占全年累計(jì)降水量的57.6%，秋季降水量占全年降水量的21.9%，冬季降水量占全年累計(jì)降水量的2.1%。此外，青海省黃河流域降水年內(nèi)分配極不均勻，進(jìn)入4月后降水逐漸增多，5—9月降水量占年降水量的84.3%左右?？傮w而言，青海省黃河流域具有年降水集中程度高，降水年內(nèi)分配不均，豐枯季分明的特點(diǎn)。

通過計(jì)算青海省黃河流域每個年代的月均降水量，繪制不同年代(20世紀(jì)50年代—21世紀(jì)20年代)降水變化(見圖2)。從圖中可以看出，無論哪個年代，降水峰值均出現(xiàn)在7—8月，降水主要集中在4—10月。經(jīng)統(tǒng)計(jì)分析發(fā)現(xiàn)，青海省黃河流域21世紀(jì)00年代與20世紀(jì)90年代相比，降水量下降了8%；而21世紀(jì)10年代 與21世紀(jì)00年代相比，降水量上升了7%。整體來看，21世紀(jì)10年代降水量達(dá)到流域的頂峰。

3.2 降水時間序列趨勢分析

采用線性回歸方法以及滑動平均方法對青海省黃河流域近63年降水時間序列進(jìn)行線性趨勢分析，線性方程的斜率為降水傾向率。同時采用M-K趨勢分析方法對線性回歸得到的趨勢分析結(jié)果進(jìn)行驗(yàn)證，以證實(shí)趨勢預(yù)測結(jié)論的合理性。

圖3顯示了青海省黃河流域1956—2018年線性趨勢、年降水量變化以及5年滑動平均過程，可以看出：(1)在95%置信區(qū)間下降水趨勢函數(shù)為f(x)=0.907x-1 352，降水增長率為(9.01 mm/10 a)，年降水量整體呈現(xiàn)出上升的波動趨勢。(2)年降水量波動性變化較為明顯，波動范圍為343.30～602.36 mm，在1967、1989、2018年達(dá)到波峰，在1956、1962、2002年達(dá)到波谷，經(jīng)歷了偏枯-偏豐-正常-偏枯-偏豐5個階段，年降水量最大值出現(xiàn)在2018年(602.36 mm)，年降水量最小值出現(xiàn)在1956年(343.30 mm)，最大值和最小值之比為1.75。(3)2014—2018年5年的平均降水量最大(503.27 mm)，1956—1960年5年的平均降水量最小(405.90 mm)，兩者相差97.37 mm。

圖2 青海省黃河流域降水量年代變化特征Fig.2 Time-varying characteristics of precipitation in the Yellow River Basin of Qinghai Province

圖3 青海省黃河流域年降水量變化及5年滑動平均過程Fig.3 The process of annual precipitation change and 5-year moving average in the Yellow River Basin of Qinghai Province

采用M-K突變檢驗(yàn)方法對青海省黃河流域1956—2018年全年以及各季節(jié)分別計(jì)算統(tǒng)計(jì)量Z值(見表2)，可以看出，近63年以來青海省黃河流域年降水量呈現(xiàn)上升的變化趨勢，這與線性回歸方法分析結(jié)果一致。對于各季節(jié)降水量的變化趨勢而言，春季和冬季降水統(tǒng)計(jì)量Z值均滿足|z|≥1.96，說明上升趨勢顯著；夏季和秋季降水統(tǒng)計(jì)量Z值不滿足|z|≤1.96，說明沒有顯著的上升趨勢。

簡言之，通過線性回歸方法和M-K方法分析青海省黃河流域1956—2018年降水時間序列可得出：青海省黃河流域年降水量有上升的變化趨勢，即未來年降水量有增加的趨勢；春季和冬季降水量有顯著上升的趨勢，而夏季和秋季降水量則無顯著的上升趨勢。

表2 青海省黃河流域年及各季節(jié)降水M-K統(tǒng)計(jì)量Z值Table 2 M-K statistic Z of annual and seasonal precipitation in the the Yellow River Basin of Qinghai Province

3.3 降水時間序列突變分析

青海省黃河流域年及各季平均降水量M-K突變檢驗(yàn)結(jié)果如圖4所示。結(jié)果表明，流域年及各季節(jié)平均降水的突變檢驗(yàn)曲線均存在交匯點(diǎn)，依據(jù)圖4(a)可知，青海省黃河流域近63年降水序列M-K突變檢驗(yàn)在0.05顯著性水平下，UF統(tǒng)計(jì)量和UB統(tǒng)計(jì)量在2004年存在交點(diǎn)，說明2004年為年降水序列顯著突變年份；依據(jù)圖4(b)、(c)、(d)、(e)可知，各季節(jié)多年平均降水M-K突變檢驗(yàn)在0.05顯著性水平下，UF統(tǒng)計(jì)量和UB統(tǒng)計(jì)量均存在多個交點(diǎn)情況，因此未明確真正的突變位置。

圖4 青海省黃河流域年及各季平均降水M-K檢驗(yàn)結(jié)果Fig.4 M-K statistic of annual and seasonal mean precipitation in the Yellow River Basin of Qinghai Province

結(jié)合Pettitt檢驗(yàn)法對青海省黃河流域年及各季平均降水量的突變年份作進(jìn)一步分析，其檢驗(yàn)結(jié)果見表3。其中，青海省黃河流域年降水統(tǒng)計(jì)量Umax=428，對應(yīng)一級突變年為2004年，突變顯著；春季和冬季平均降水量存在突變情況，通過了0.05水平的顯著性檢驗(yàn)，說明突變顯著；夏季和秋季平均降水量未通過0.05水平的顯著性檢驗(yàn)，因此突變情況不顯著；這與M-K突變檢驗(yàn)結(jié)果相同。兩種方法綜合分析表明，青海省黃河流域歷史年及各季節(jié)平均降水存在突變情況，且2004年為該流域大的降水轉(zhuǎn)折年。

表3 青海省黃河流域年及各季平均降水量Pettitt突變檢驗(yàn)結(jié)果Table 3 Inspection results of Pettitt abrupt change of annual and seasonal mean precipitation in the Yellow River Basin of Qinghai Province

3.4 降水時間序列周期分析

本文選用Morlet小波作為母小波對青海省黃河流域年降水時間序列進(jìn)行周期分析，由于M-K檢驗(yàn)得到2004年為降水序列突變年份，因此去除2004年的原始降水?dāng)?shù)據(jù)，進(jìn)行插值處理得到新降水?dāng)?shù)據(jù)，并對時間序列進(jìn)行距平處理和延拓處理。Morlet小波變換系數(shù)的實(shí)部反映了不同時間尺度下小波域內(nèi)數(shù)據(jù)的波動特征，實(shí)部的正負(fù)取值交替反映了時間域內(nèi)降水序列豐枯特征，小波方差則為小波系數(shù)的平方值在時間域上的積分，可以確定降水序列存在的多級主周期。據(jù)此，做出年降水序列小波系數(shù)等值線圖(見圖5)和小波方差圖(見圖6)。

圖5 青海省黃河流域年平均降水距平序列 小波變換系數(shù)實(shí)部等值線圖Fig.5 Isoline for the real components of wavelet transform coefficient chart for the annual mean precipitation anomaly in the Yellow River Basin of Qinghai Province

從圖5可以看出，青海省黃河流域年降水量在不同時間尺度表現(xiàn)出不同的周期變化和豐枯交替的變化過程，形成各種尺度正負(fù)相間的振蕩中心，存在明顯的年際變化和年代際變化。從上至下分析得出青海省黃河流域63年(1956—2018年)降水存在23～32 a，15～20 a，9～13 a以及4～6 a的周期變化規(guī)律。從23～32 a較大尺度分析，青海省黃河流域降水存在3個準(zhǔn)振蕩周期，1959年以前，1971—1977年和1996—2003均表現(xiàn)為枯水年，而1960—1970年，1978—1995年和2003—2015年表現(xiàn)為豐水年。青海省黃河流域降水在15～20 a尺度周期集中表現(xiàn)在1986年以后，存在3個準(zhǔn)振蕩周期，1988—1990年，2000—2005年和2013—2016年以后均表現(xiàn)為枯水年，而1991—1999年，2006—2012年和2016年以后均表現(xiàn)為豐水年。9～13 a和4～6 a的年際變化尺度具有全域性，豐枯交替頻繁。不同時間尺度對應(yīng)的降水結(jié)構(gòu)有所不同，較大時間尺度的降水周期變化嵌套小時間尺度的降水變化，同一段時域內(nèi)受多級周期影響。進(jìn)一步結(jié)合圖6小波方差圖確定降水序列存在的多級主周期，可以確定青海省黃河流域的降水序列存在四級主周期，即30、15、11和6 a時間尺度。

圖6 青海省黃河流域年平均降水距平序列小波方差圖

3.5 降水時間序列預(yù)測

采用ARIMA模型對青海省黃河流域1956—2018年63年降水時間序列進(jìn)行模型選擇及降水預(yù)測，主要分為5個步驟：(1)序列的平穩(wěn)化處理。根據(jù)ADF檢驗(yàn)結(jié)果可知，此時間序列為非平穩(wěn)序列，進(jìn)行一階差分后ADF檢驗(yàn)結(jié)果顯著，P-value=0.01<0.05，說明一階差分處理后的序列為平穩(wěn)序列，即ARIMA模型中的參數(shù)d取1。(2)繪制自相關(guān)圖和偏自相關(guān)圖。繪制一階差分后平穩(wěn)序列的自相關(guān)圖(ACF圖)和偏自相關(guān)圖(PACF圖)(見圖7)。(3)模型定階。從PACF圖中可以看出在5%的顯著性水平下，偏自相關(guān)系數(shù)顯著不為零的個數(shù)為3，之后趨于零，并呈現(xiàn)出拖尾現(xiàn)象，可以判斷參數(shù)P值取3或者2(取2是因?yàn)樽詈笠粋€顯著超過0的階數(shù)是由于偶然出現(xiàn)的)；從ACF圖中可看出，在5%的顯著水平下，自相關(guān)系數(shù)顯著為0的個數(shù)為5，可以判斷參數(shù)q值取5或4(取4是因?yàn)樽詈笠粋€顯著超過0的階數(shù)是由于偶然出現(xiàn)的)。因此，模型可初步確定為ARIMA(2,1,4)、ARIMA(2,1,5)、ARIMA(3,1,4)、ARIMA(3,1,5)。其中，ARIMA(3,1,4)模型的AIC值為672.38、ARIMA(2,1,5)模型的AIC值為674.15、ARIMA(3,1,4)模型的AIC值為674.22、ARIMA(3,1,5)模型的AIC值為676.15。由赤池信息準(zhǔn)則(Akaike Information Criterion, AIC)可知，在一組可供選擇的模型中選擇最佳模型時，通常選擇AIC值最小的模型。因此，最終選擇ARIMA(2,1,4)為青海省黃河流域降水預(yù)測最佳模型。(4)模型的檢驗(yàn)。對ARIMA(2,1,4)模型進(jìn)行殘差白噪聲檢驗(yàn)，殘差的Ljung-Box檢驗(yàn)結(jié)果中P-value=0.88>0.05,說明模型的殘差沒用通過顯著性檢驗(yàn)，可以認(rèn)為殘差的自相關(guān)系數(shù)為零，ARIMA(2,1,4)模型能較好地?cái)M合本序列。(5)模型的預(yù)測。本文采用ARIMA(2,1,4)模型對青海省黃河流域2019—2023年降水量進(jìn)行預(yù)測，在95%置信區(qū)間下的2019—2023年降水量預(yù)測結(jié)果如表4所示。

3.6 預(yù)測結(jié)果驗(yàn)證

為驗(yàn)證ARIMA(2,1,4)模型對青海省黃河流域2019—2023年降水序列預(yù)測結(jié)果的合理性，本研究依舊采用線性回歸方法及M-K趨勢分析方法對青海省黃河流域1956—2023年降水序列進(jìn)行趨勢分析，其中，線性回歸方法分析結(jié)果如圖8所示，可以看出，在95%置信區(qū)間下降水趨勢函數(shù)為f(x)=0.922x-1 381.2，降水增長率為(0.922 mm/10 a)，年降水量整體呈現(xiàn)出上升的波動趨勢；通過M-K趨勢檢驗(yàn)法獲得統(tǒng)計(jì)量Z值為3.08，滿足|z|≥1.96，說明青海省黃河流域1956—2023年降水序列上升趨勢顯著。綜上所述，可以看出3.5節(jié)中采用ARIMA(2,1,4)模型對青海省黃河流域2019—2023年降水序列預(yù)測結(jié)果符合3.2節(jié)中得出的青海省黃河流域年降水時間序列未來有增加趨勢的結(jié)論。

表4 95%置信區(qū)間下青海省黃河流域 2019—2023年降水量預(yù)測結(jié)果Table 4 Precipitation prediction results of the Yellow River Basin in Qinghai Province in 2019—2023 under 95% confidence interval

圖7 青海省黃河流域年降水量一階差分序列ACF(a)和PACF(b)圖Fig.7 ACF(a) and PACF(b) maps of the first difference series of annual precipitation in the Yellow River Basin of Qinghai Province

圖8 青海省黃河流域1956—2023年降水時間序列趨勢分析Fig.8 Tendency analysis of precipitation time series from 1956 to 2023 in the Yellow River Basin of Qinghai Province

4 結(jié)論

相較于青海省降水時空分布特性相關(guān)研究，本研究運(yùn)用更為系統(tǒng)的包括降水分配、趨勢、突變、周期以及預(yù)測在內(nèi)的分析方法，并將趨勢檢驗(yàn)結(jié)果作為降水預(yù)測是否合理的一個度量標(biāo)準(zhǔn)，起到了雙重驗(yàn)證的作用，豐富了青海省東南部地區(qū)降水特性的研究成果。相較于黃河流域降水特性相關(guān)研究，本研究僅選取青海省區(qū)域，研究范圍更加精細(xì)化，且為黃河源區(qū)降水趨勢預(yù)測提供了ARIMA(2,1,4)模型，可為水資源規(guī)劃及農(nóng)業(yè)發(fā)展做出參考。本研究在降水年內(nèi)分配及趨勢分析方面與白淑英等[4]研究結(jié)果大致相同；在降水周期分析方面與王玉娟等[24]研究略有差異，而降水趨勢變化研究大致相同，這說明，不同的研究區(qū)域、研究時間和研究尺度都會導(dǎo)致降水量周期分析的差異。在降水預(yù)測方面，本研究只針對了單一變量的降水預(yù)測，可以考慮多氣溫、蒸發(fā)量、海拔高程等多因子對降水預(yù)測的影響。本研究利用M-K檢驗(yàn)、Pettitt檢驗(yàn)、小波分析等方法對青海省黃河流域1956—2018年降水量進(jìn)行趨勢、突變、周期分析，運(yùn)用ARIMA(2,1,4)模型對青海省黃河流域2019—2023年降水量進(jìn)行預(yù)測，并采用線性回歸及M-K檢驗(yàn)方法對預(yù)測趨勢進(jìn)行了合理驗(yàn)證。主要結(jié)論如下：

(1)年降水序列表現(xiàn)為集中程度高、年內(nèi)分配不均、豐枯季明顯的特點(diǎn)，降水集中在夏秋兩季， 5—9月是降水多發(fā)期，預(yù)計(jì)未來該流域年降水量將處于偏豐或正常狀態(tài)；1950s—2010s之間的任一年代，該流域月均降水均集中在4—10月，且月均降水量未來有上升的變化趨勢。

(2)年降水序列存在上升趨勢，趨勢函數(shù)為f(x)=0.907x-1 352，春、冬季降水量有顯著上升的趨勢，夏、秋季降水量無顯著的上升趨勢；2004年為該流域大的降水轉(zhuǎn)折年，采用M-K檢驗(yàn)與Pettitt檢驗(yàn)兩種方法對該流域進(jìn)行突變識別效果優(yōu)于單一使用M-K檢驗(yàn)分析。

(3)年降水序列存在23～32 a，15～20 a，9～13 a以及4～6 a的周期變化規(guī)律，且不同尺度的豐枯變化與降水豐枯分析較為吻合；降水序列的四級主周期分別對應(yīng)30、15、11和6 a的時間尺度。

(4)ARIMA(2,1,4)模型能合理地對該流域未來5年降水序列進(jìn)行預(yù)測，從而對當(dāng)?shù)厮Y源規(guī)劃和管理做出指導(dǎo)。