黃河源區(qū)氣象干旱與水文干旱關(guān)聯(lián)性分析

石 朋，詹慧婕，瞿思敏，馮 進(jìn)，管曉祥

(河海大學(xué)水文水資源學(xué)院，江蘇 南京 210098)

全球氣候變暖背景下的干旱問(wèn)題日益突出，干旱導(dǎo)致的水資源短缺進(jìn)一步加劇了人類(lèi)的用水緊張程度[1]，同時(shí)導(dǎo)致土地退化以及極端天氣和復(fù)合型災(zāi)害的發(fā)生[2-3]。1972—1997年黃河經(jīng)歷了漫長(zhǎng)的干涸時(shí)期，并在20世紀(jì)90年代初出現(xiàn)干旱歷時(shí)延長(zhǎng)、干旱頻率增加的現(xiàn)象[4]。1997年在中國(guó)北方發(fā)生的嚴(yán)重干旱導(dǎo)致黃河出現(xiàn)為期226 d的零徑流，造成了嚴(yán)重的社會(huì)經(jīng)濟(jì)損失[5]。準(zhǔn)確表征黃河流域干旱特征及變化規(guī)律，對(duì)確保黃河流域生態(tài)保護(hù)與高質(zhì)量發(fā)展具有極為重要的意義。

干旱分為4種類(lèi)型，氣象干旱、水文干旱、農(nóng)業(yè)干旱和社會(huì)經(jīng)濟(jì)干旱[4]。降水的虧缺會(huì)導(dǎo)致氣象干旱，地表水及地下水的虧缺會(huì)引發(fā)水文干旱，這兩類(lèi)干旱對(duì)人類(lèi)社會(huì)產(chǎn)生了廣泛的影響。一般使用干旱指數(shù)來(lái)表征氣象和水文干旱，標(biāo)準(zhǔn)化降水指數(shù)(standardized precipitation index，SPI)只需要降水序列作為輸入，并且能將降水分布轉(zhuǎn)化為正態(tài)分布，從而克服非標(biāo)準(zhǔn)化分布帶來(lái)的差異，因此被廣泛應(yīng)用于氣象干旱分析[6]。同樣，標(biāo)準(zhǔn)化徑流指數(shù)(standardized runoff index，SRI)只需流量序列，且計(jì)算簡(jiǎn)便，因而被應(yīng)用于水文干旱分析[7]。

一般情況下首先發(fā)生氣象干旱，氣象干旱會(huì)向水文干旱和農(nóng)業(yè)干旱傳播[8]，因而越來(lái)越多的學(xué)者研究氣象干旱與水文干旱之間的聯(lián)系。準(zhǔn)確描述氣象干旱和水文干旱的關(guān)系可以在氣象干旱的初期及時(shí)向決策管理者預(yù)警，減輕水文干旱帶來(lái)的社會(huì)經(jīng)濟(jì)損失。Li等[9]采用小波分析法發(fā)現(xiàn)，沙潁河上游氣象干旱和水文干旱之間的傳播時(shí)間隨季節(jié)變化明顯；Dehghani等[10]基于阿基米德Copula氣象干旱指數(shù)對(duì)水文干旱進(jìn)行概率性預(yù)報(bào)；Yu等[11]研究了淮河上游嵌套流域的氣象干旱到水文干旱的傳播時(shí)間，發(fā)現(xiàn)流域越大，氣象干旱引發(fā)的水文干旱持續(xù)時(shí)間越長(zhǎng)。但是，由于不同流域降水條件及調(diào)蓄能力的不同，氣象干旱到水文干旱的傳播存在較大差異，目前對(duì)黃河源區(qū)開(kāi)展的研究局限于氣象及水文干旱特征、干旱的時(shí)空變化規(guī)律等[12-13]，黃河源區(qū)氣象干旱與水文干旱的關(guān)聯(lián)性及干旱傳播時(shí)間尚未探明。

劉永佳等[14]利用皮爾遜相關(guān)分析方法分析了黃土高原氣象干旱到水文干旱的季節(jié)性傳播規(guī)律；楊銘珂等[15]利用皮爾遜相關(guān)分析方法，發(fā)現(xiàn)黔中地區(qū)不同季節(jié)氣象和水文干旱的最佳時(shí)間尺度和滯后時(shí)間不同；Wang等[16]利用交叉小波及小波相干方法探究了黃河流域氣象干旱到水文干旱的傳播特征。目前研究干旱的傳播規(guī)律多使用皮爾遜相關(guān)分析方法，結(jié)合1～2種小波方法進(jìn)行干旱周期變化規(guī)律的研究，而利用多種小波方法探明區(qū)域化干旱傳播特征及氣象干旱與水文干旱的內(nèi)在聯(lián)系性的研究相對(duì)較少。

黃河源區(qū)是重要的水資源保護(hù)地，本研究聚焦黃河源區(qū)區(qū)域化氣象干旱與水文干旱的關(guān)聯(lián)性及水文干旱的滯后效應(yīng)，側(cè)重對(duì)黃河源區(qū)區(qū)域化規(guī)律進(jìn)行研究，利用最優(yōu)的SPI-n(n為月時(shí)間尺度，n=1,2,…,24)和SRI擬合結(jié)果，探究黃河源區(qū)氣象干旱向水文干旱的最佳傳播時(shí)間，基于連續(xù)小波變換、交叉小波變換和小波相干來(lái)分析黃河源區(qū)氣象干旱和水文干旱的周期性變化規(guī)律及其關(guān)聯(lián)性，以期為流域干旱災(zāi)害預(yù)警及防治提供參考。

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

黃河源區(qū)位于青藏高原東北部(95°50′E～103°30′E，32°30′N(xiāo)～36°00′N(xiāo))，流域面積約12.19萬(wàn) km2，主要控制水文站為唐乃亥水文站。黃河源區(qū)地勢(shì)起伏較大，從西到東地勢(shì)逐漸下降，平均高程約為4 425 m，多年平均降水約510 mm[17]，汛期集中在夏季，多年平均徑流量約占黃河流域的38%。黃河源區(qū)下墊面類(lèi)型主要以高寒草甸、季節(jié)性凍土和多年凍土區(qū)為主[18]，流域的水文站、氣象站及水系分布如圖1所示。

圖1 黃河源區(qū)水系分布

本文采用的90 m分辨率SRTM DEM原始高程數(shù)據(jù)來(lái)自地理空間數(shù)據(jù)云(http://www.gscloud.cn/)。黃河源區(qū)內(nèi)部7個(gè)氣象站(瑪多、興海、達(dá)日、河南、瑪曲、若爾蓋、紅原)1967—2020年的日降水?dāng)?shù)據(jù)來(lái)自中國(guó)氣象數(shù)據(jù)網(wǎng)(htttp://data.cma.cn/)，經(jīng)過(guò)日累加處理得到月降水?dāng)?shù)據(jù)。唐乃亥水文站1967—2020年的月流量資料來(lái)源于水文年鑒，數(shù)據(jù)經(jīng)過(guò)嚴(yán)格質(zhì)量控制和檢查。

2 研究方法

2.1 Kolmogorov-Smirnov檢驗(yàn)

Kolmogorov-Smirnov(KS)檢驗(yàn)是一種擬合優(yōu)度檢驗(yàn)方法。根據(jù)一個(gè)給定的概率分布，利用經(jīng)驗(yàn)累積分布函數(shù)(empirical cumulative distribution function，ECDF)和累積分布函數(shù)(cumulative distribution function，CDF)之間的差異計(jì)算擬合優(yōu)度。KS檢驗(yàn)包含一個(gè)統(tǒng)計(jì)量D，若D值大于某一臨界值，則拒絕假設(shè)的概率分布[19-20]。

2.2 氣象干旱與水文干旱指標(biāo)

SPI是一種常見(jiàn)的氣象干旱指標(biāo)[21]，一般選用Gamma分布來(lái)擬合降水時(shí)間序列[22]。SRI被用來(lái)刻畫(huà)水文干旱[7]，計(jì)算SRI的步驟和計(jì)算SPI的類(lèi)似，需要將降水時(shí)間序列替換成流量時(shí)間序列，并且分布函數(shù)因區(qū)域而異，流量序列可能適用于Gamma分布、皮爾遜Ⅲ型(P-Ⅲ型)分布、廣義極值(GEV)分布、對(duì)數(shù)正態(tài)(Log-Normal)分布等。

2.3 皮爾遜相關(guān)系數(shù)

皮爾遜相關(guān)系數(shù)可以表示兩個(gè)變量之間的相關(guān)關(guān)系[23]，相關(guān)系數(shù)越大，表示兩個(gè)變量線性相關(guān)越強(qiáng)，范圍在[-1,1]之間。

2.4 面平均降水量計(jì)算方法

由于降水空間分布不均勻，不同的面平均降水量計(jì)算方法可能會(huì)影響計(jì)算結(jié)果，本文挑選3種常用的降水空間插值方法來(lái)計(jì)算面平均降水量，分別為算術(shù)平均法、泰森多邊形法和距離平方倒數(shù)法[24]。

2.5 小波分析方法

連續(xù)小波變換可以用來(lái)分析時(shí)間序列中的局部間歇性振蕩，同時(shí)可以分析出時(shí)間序列的周期性變化規(guī)律[25]；交叉小波變換和小波相干可以用來(lái)檢驗(yàn)兩個(gè)時(shí)間序列之間的相關(guān)關(guān)系，其中，交叉小波變換用來(lái)檢測(cè)兩個(gè)序列高能量區(qū)域的相關(guān)性[16]，小波相干則用來(lái)檢測(cè)兩個(gè)序列間低能量區(qū)域的相關(guān)性[26]。

3 結(jié)果與分析

3.1 最優(yōu)流量分布篩選

由于氣候變化、人類(lèi)活動(dòng)等因素，流量的概率分布存在區(qū)域性差異[16]，需要依據(jù)實(shí)測(cè)數(shù)據(jù)篩選出最適宜的流量分布函數(shù)。本研究選擇6種常見(jiàn)的分布函數(shù)[19]來(lái)擬合唐乃亥站1967—2020年的月流量，包括GEV分布、威布爾(Weibull)分布、廣義帕累托(廣義Pareto)分布、Log-Normal分布、對(duì)數(shù)Logistic(Log-Logistic)分布、P-Ⅲ型分布。

利用KS檢驗(yàn)來(lái)驗(yàn)證6種分布函數(shù)與經(jīng)驗(yàn)分布的逐月差異，差異程度由假設(shè)檢驗(yàn)結(jié)果的判定參數(shù)p值來(lái)表征。p越接近1，表示分布越接近真實(shí)情況，在p<0.05的閾值下[9,23]，除了廣義Pareto分布被拒絕，其他5個(gè)分布均通過(guò)KS檢驗(yàn)，如圖2所示，其中紅空心點(diǎn)代表異常值，黑實(shí)心點(diǎn)代表均值，箱內(nèi)黑線代表中位數(shù)。

分析發(fā)現(xiàn)，GEV分布p最小值(0.321)和平均值(0.796)都高于其他分布；Weibull分布雖然也通過(guò)了KS檢驗(yàn)，但在唐乃亥水文站的應(yīng)用上并不理想；Log-Normal和Log-Logistic分布擬合結(jié)果較為接近，P-Ⅲ型分布擬合結(jié)果略次于Log-Normal和Log-Logistic分布。

根據(jù)p值選擇每個(gè)月最優(yōu)的流量CDF如圖3所示，GEV分布在6個(gè)月份中表現(xiàn)出最優(yōu)的性能，綜合考慮選擇GEV分布作為最優(yōu)流量CDF進(jìn)行后續(xù)計(jì)算。

圖2 5種分布經(jīng)過(guò)KS檢驗(yàn)后的p值

3.2 SPI和SRI的計(jì)算

由于降水具有空間不均勻性[27]，為了使計(jì)算更切合實(shí)際，減少降水輸入的不確定性，利用算術(shù)平均法、泰森多邊形法和距離平方倒數(shù)法分別處理1967—2020年的月降水量，得到流域月尺度面降水量。

根據(jù)3種降水插值方法得到的月尺度流域面降水量，使用Gamma分布計(jì)算氣象干旱指數(shù)SPI-n(n=1,2,…,24)。以往的研究結(jié)果表明[14]，1月尺度的SRI指數(shù)最具有代表性，所以最優(yōu)流量分布利用GEV分布函數(shù)，計(jì)算水文干旱指數(shù)SRI-1。利用皮爾遜相關(guān)系數(shù)，分析由3種降水插值方法計(jì)算的SPI-n與SRI-1之間的相關(guān)關(guān)系。結(jié)果表明，3種方法得出的SPI-9與SRI-1擬合效果最優(yōu)，皮爾遜相關(guān)系數(shù)分別為0.68、0.67、0.68，3種方法計(jì)算結(jié)果較接近，黃河源區(qū)氣象干旱向水文干旱的傳播時(shí)間為9月。基于以上結(jié)果，本次研究選擇距離平方倒數(shù)法進(jìn)行后續(xù)研究。

基于距離平方倒數(shù)法計(jì)算的SPI-9和基于GEV分布函數(shù)計(jì)算的SRI-1如圖4所示。由于SPI-9采用了累計(jì)9月的降水進(jìn)行計(jì)算，為了保證流量序列和降水序列的一致性，降水和流量序列統(tǒng)一取1968—2020年。根據(jù)SPI和SRI的值劃分干濕等級(jí)[14]，結(jié)果表明，SPI-9和SRI-1表征的干濕程度基本一致，氣象干旱與水文干旱發(fā)生年份幾乎重疊。黃河源區(qū)在1969—1971年和1973年發(fā)生了嚴(yán)重氣象干旱和中度水文干旱事件，在1977—1978年發(fā)生了嚴(yán)重氣象干旱和水文干旱事件，在1985—2005年間歇性地發(fā)生了氣象和水文干旱事件，每次發(fā)生干旱的持續(xù)時(shí)間較長(zhǎng)，在2003年發(fā)生了極端氣象和水文干旱事件，在2005—2020年發(fā)生干旱的嚴(yán)重程度減小，干旱發(fā)生頻率也較1985—2005年減少。周帥等[5]指出2005年以后由于降水減少幅度緩解的原因，黃河流域干旱嚴(yán)重程度減小，這與本文的研究結(jié)論一致。

(a)1月 (b)2月 (c)3月 (d)4月

(e)5月 (f)6月 (g)7月 (h)8月

(i)9月 (j)10月 (k)11月 (l)12月

(a)SPI-9

(b)SRI-1

3.3 氣象干旱和水文干旱關(guān)聯(lián)性

先用連續(xù)小波變換分析氣象干旱和水文干旱的周期性變化規(guī)律，再用交叉小波變換分析高能量區(qū)域的氣象干旱和水文干旱的聯(lián)系，對(duì)于低能量區(qū)域，則選擇小波相干法分析兩者的聯(lián)系。

圖5給出了1968—2020年SPI-9和SRI-1的連續(xù)小波功率譜，圖中影響錐表示不受小波光譜邊緣效應(yīng)影響的彩色區(qū)域，粗輪廓線表示相對(duì)于紅色噪聲的95%置信水平。結(jié)果表明，對(duì)于SPI-9序列，在95%的置信水平下，1972年5月至1978年2月、1988年3月至1992年2月、1997年12月至2002年10月、2002年10月至2008年12月的小波高功率分別在37～50月、29～42月、12～16月、18～26月之間，持續(xù)時(shí)間約為4～6 a。1971年9月至1973年4月、1988年7月至1990年2月、2015年5月至2017年9月的小波高功率都在17～22月之間，持續(xù)時(shí)間較短，約為2 a。對(duì)于SRI-1序列，在1978—2016年間，2～8月周期的小波功率值較大，持續(xù)時(shí)間較短，表現(xiàn)為準(zhǔn)周期的間歇性振蕩。SRI-1最明顯的周期是78～105月，持續(xù)時(shí)間是1977年6月至1986年7月。SRI-1在1988年9月至1992年2月的小波高功率在37～47月之間，與SPI-9在相同時(shí)間段的周期性較接近。

(a)SPI-9

(b)SRI-1

在37～42月的周期，SPI-9和SRI-1序列都在1988年9月至1992年2月有較高的功率值，說(shuō)明黃河源區(qū)的氣象干旱和水文干旱具有相似的周期特征。

圖6(a)為SPI-9和SRI-1的交叉小波功率譜圖，向右箭頭表示SPI-9和SRI-1序列相位變化相同，向左箭頭表示兩序列相位變化相反。SPI-9和SRI-1序列在1973年4月至1979年4月和在1978年9月至1984年8月分別呈現(xiàn)出37～53月和83～105月的顯著正相關(guān)，在1987年8月至1994年3月間呈現(xiàn)出31～50月的顯著正相關(guān)，在2002年3月至2009年6月間呈現(xiàn)出14～25月的顯著正相關(guān)。SPI-9和SRI-1的小波相干圖如圖6(b)所示，向右的黑色箭頭表明SPI-9和SRI-1序列具有相同的相位，且存在比較連續(xù)的小波顯著相干，具體表現(xiàn)在1971—2016年。SRI-1的振蕩可以較好地用SPI-9來(lái)表述，周期為16～176月不等。

連續(xù)小波分析、交叉小波分析和小波相干分析結(jié)果表明，在不同的時(shí)間段，SPI-9和SRI-1均存在明顯的相關(guān)性，SPI-9序列可以用于黃河源區(qū)的水文干旱監(jiān)測(cè)。

4 結(jié) 論

a.選擇的6種分布中，廣義Pareto分布函數(shù)被KS檢驗(yàn)拒絕，其他5種分布均通過(guò)檢驗(yàn)，GEV分布表現(xiàn)最優(yōu)，然后依次是Log-Normal分布、Log-Logistic分布、P-Ⅲ型分布，Weibull分布表現(xiàn)最差。

(a)交叉小波

(b)小波相干

b.3種降水插值方法中，算術(shù)平均法、泰森多邊形法、距離平方倒數(shù)法計(jì)算結(jié)果接近，在黃河源區(qū)適用效果均良好。由距離平方倒數(shù)法計(jì)算的SPI-9和GEV分布計(jì)算的SRI-1擬合效果最優(yōu)，皮爾遜相關(guān)系數(shù)為0.68，氣象干旱和水文干旱發(fā)生的年份基本一致，兩者存在密切聯(lián)系。

c.連續(xù)小波變換分析表明SPI-9和SRI-1具有類(lèi)似的周期特征。交叉小波變換表明SPI-9和SRI-1在1973年4月至1979年4月、1978年9月至1984年8月、1987年8月至1994年3月、2002年3月至2009年6月分別呈現(xiàn)出37～53月、83～105月、31～50月、14～25月的顯著正相關(guān)。小波相干表明SPI-9能較好地表述SPI-1的振蕩性，在1971—2016年呈現(xiàn)出16～176月的周期。