任梅芳,徐宗學(xué),龐 博

(1. 北京師范大學(xué)水科學(xué)研究院,北京 100875；2. 城市水循環(huán)與海綿城市技術(shù)北京市重點(diǎn)實(shí)驗(yàn)室,北京 100875)

近年來,變化環(huán)境下河川徑流發(fā)生了明顯的變化,在氣候變化和人類活動(dòng)雙重影響下,干旱地區(qū)流域徑流顯著減少[1- 2],濕潤地區(qū)流域徑流增加[3],而城市流域的洪澇事件增多增強(qiáng)[4]。變化環(huán)境下流域徑流演變及歸因是國際水文科學(xué)協(xié)會(huì)(International Association of Hydrological Sciences,IAHS)最新10年科學(xué)研究計(jì)劃的核心研究主題之一[5],研究水文過程變化及其驅(qū)動(dòng)機(jī)理對(duì)評(píng)價(jià)氣候變化和人類活動(dòng)對(duì)全球和區(qū)域水循環(huán)演變過程的影響十分重要[6]。通常情況下,降水是流域洪水變化最主要的驅(qū)動(dòng)要素[7],在全球氣候變暖及城市化高度發(fā)展的背景下,全球及區(qū)域降水發(fā)生了很大的改變,同時(shí),流域洪水的變化也受到其他因素的影響,如氣溫、地下水變化及人類活動(dòng)[8]。近年來,有關(guān)學(xué)者對(duì)全球及區(qū)域徑流的驅(qū)動(dòng)機(jī)制進(jìn)行了分析。Zhang等[9]對(duì)中國降水和徑流的關(guān)系研究發(fā)現(xiàn),不同流域的人類活動(dòng)伴隨著降水的變化對(duì)徑流表現(xiàn)出不同的影響；Wang和Hejazi[10]研究氣候變化和人類活動(dòng)對(duì)美國413個(gè)流域徑流的影響發(fā)現(xiàn),氣候變化導(dǎo)致流域徑流量增加,而由人類活動(dòng)導(dǎo)致的徑流變化在不同流域具有不同的特性；王國慶等[11]研究表明,人類活動(dòng)對(duì)中國北方河川徑流的影響大于氣候變化帶來的影響,而在中國淮河以南,由氣候變化引起的徑流變化占主導(dǎo)地位；Li等[12]研究表明,降水、氣溫、農(nóng)業(yè)用地和水庫的修建對(duì)黃河流域流量的變化有顯著影響,由水庫和農(nóng)業(yè)用地相關(guān)因素變化引起流量的變化是錯(cuò)綜復(fù)雜的,量化每個(gè)驅(qū)動(dòng)要素對(duì)流量變化的影響仍然面臨著較大的挑戰(zhàn)。在城市流域中,城市化的發(fā)展對(duì)城市流域流量的變化起到至關(guān)重要的作用。城市化進(jìn)程使得城市區(qū)域人口急劇上升,城市區(qū)域不透水面積顯著增加,對(duì)城市流域的水循環(huán)過程造成一定的影響,而在高度城市化地區(qū),多種人類活動(dòng)共同作用使得城市洪水的發(fā)生機(jī)理更為復(fù)雜。

大多數(shù)傳統(tǒng)的水文變量統(tǒng)計(jì)分析是基于一致性假設(shè)基礎(chǔ)之上的,換言之,假設(shè)水文變量在時(shí)間序列上沒有發(fā)生突變[13- 14]。然而,由于受到氣候變化和人類活動(dòng)的影響,無法保證水文資料在過去、現(xiàn)在和未來都服從同一分布,這就導(dǎo)致了水文變量在時(shí)間序列上的非一致性[15]。大量研究表明一致性假設(shè)在水資源規(guī)劃及危險(xiǎn)性評(píng)估中不再適用,非一致性洪水頻率分析已成為前沿水文科學(xué)問題[16- 17]。目前,國內(nèi)常采用還原法或還現(xiàn)法來考慮非一致性水文頻率計(jì)算[18]。還原法是指將非一致性時(shí)間序列修正為滿足一致性條件的時(shí)間序列,還現(xiàn)法是將水文資料修正為在現(xiàn)狀條件下滿足一致性要求的時(shí)間序列[16],而大量研究表明傳統(tǒng)的還原法或還現(xiàn)法存在“還原失真”和“還原失效”的現(xiàn)象[18- 19]。對(duì)于非一致性水文資料的分析常采用3種方法,分別為基于混合分布、基于時(shí)變矩和基于條件概率分布的非一致性水文頻率分析方法。由于時(shí)變矩法可以靈活地選取與洪水變化相關(guān)的解釋變量來描述洪水的變化,并且能夠?qū)Ψ且恢滦赃M(jìn)行成因分析,因此,時(shí)變矩法成為近年來國內(nèi)外研究非一致性洪水頻率分析最常采用的方法之一。廣義加和模型(Generalized Additive Models for Location,Scale and Shape,GAMLSS)是一種引入位置、尺度和形狀的廣義參數(shù)可加模型,該模型能夠靈活描述統(tǒng)計(jì)參數(shù)和解釋變量之間的線性或非線性關(guān)系,為時(shí)變矩法提供了強(qiáng)大和便捷的工具[20- 23]。

北京市作為中國的政治、文化和經(jīng)濟(jì)中心,是城鎮(zhèn)化程度最高、經(jīng)濟(jì)最發(fā)達(dá)的地區(qū)之一[24]。溫榆河流域位于北京市的核心區(qū)域,該流域流經(jīng)北京市昌平、海淀、順義、朝陽和通州5個(gè)行政區(qū)。近年來,北京市洪澇事件頻繁發(fā)生,全面科學(xué)識(shí)別不同等級(jí)洪水發(fā)生的關(guān)鍵要素,研究氣候變化和城市化如何影響城市流域流量的變化對(duì)城市防洪工作具有重要的科學(xué)意義。

本文以高度城市化的溫榆河流域?yàn)檠芯繉?duì)象,借助GAMLSS模型,對(duì)溫榆河流域夏季全概率洪水演變的關(guān)鍵驅(qū)動(dòng)要素進(jìn)行識(shí)別,揭示變化環(huán)境下城市流域不同頻率洪水演變的主要驅(qū)動(dòng)機(jī)制,以期對(duì)中國海綿城市建設(shè)和城市洪水管理提供科學(xué)依據(jù)。

1 研究區(qū)概況

溫榆河流域位于北京市中部,是北京市城區(qū)主要的防洪和排水河道,全長47.5 km,流域面積為2 478 km2,被譽(yù)為北京市的“母親河”,流域?qū)儆跍貛Т箨懶约撅L(fēng)氣候,夏季炎熱多雨,冬季寒冷干燥。多年平均降雨量為619.0 mm,大多數(shù)雨量集中在6—9月,約占全年的84%。流域位于北京市核心區(qū)域,自20世紀(jì)80年代以來隨著北京市城市化的高度發(fā)展,溫榆河流域的土地利用發(fā)生了顯著變化,不透水面積顯著增加,城市化率從1985年的約4%增長到2016年的約42%[25- 26]。

2 數(shù)據(jù)與方法

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

本研究采用溫榆河流域出口北關(guān)閘1985—2016年夏季日尺度流量數(shù)據(jù)和流域內(nèi)10個(gè)氣象站日尺度降水?dāng)?shù)據(jù),水文、氣象站地理位置如圖1所示。將流域內(nèi)10個(gè)氣象站日降雨量的平均值作為流域平均降雨量,并計(jì)算夏季累計(jì)日降雨量作為降雨量驅(qū)動(dòng)因子,水文數(shù)據(jù)來自北京市水文總站。氣溫?cái)?shù)據(jù)來源于中國區(qū)域高時(shí)空分辨率地面氣象要素驅(qū)動(dòng)數(shù)據(jù)集(China Meteorological Forcing Dataset,CMFD),該數(shù)據(jù)集是中國科學(xué)院青藏高原研究所研發(fā)的近地面氣象與環(huán)境要素再分析數(shù)據(jù)集,CMFD數(shù)據(jù)集包括降水、近地表氣溫、地表氣壓、近地表空氣比濕、近地面全風(fēng)速、向下短波輻射和向下長波輻射等7個(gè)氣象變量,時(shí)間分辨率為3 h,空間分辨率為0.1°×0.1°,數(shù)據(jù)下載地址為http:∥westdc.westgis.ac.cn/data/7a35329c- c53f- 4267- aa07- e0037d913a21[27]；將CMFD數(shù)據(jù)疊加轉(zhuǎn)換為日尺度氣溫?cái)?shù)據(jù),計(jì)算研究區(qū)域季節(jié)氣溫平均值作為流域氣溫。本研究采用北京市逐年地下水埋深反映溫榆河流域地下水埋深狀況,該數(shù)據(jù)來源于北京市水文總站。采用全球逐年人工不透水面積地圖來分析溫榆河在過去幾十年間人工不透水面積的變化過程[28],該數(shù)據(jù)集是采用Landsat遙感影像數(shù)據(jù),借助于Google Earth Engine平臺(tái)建立的1985—2018年全球不透水面積地圖,數(shù)據(jù)空間分辨率為30 m×30 m,數(shù)據(jù)來源于http:∥data.ess.tsinghua.edu.cn/。圖2為溫榆河流域在1985—2016年不透水面積空間變化過程示意圖(紅色區(qū)域?yàn)椴煌杆?。

圖1 溫榆河流域位置及水文氣象站點(diǎn)分布Fig.1 Location of the Wenyu River basin and the spatial distribution of hydrological and meteorological stations

圖2 研究期溫榆河不透水面積空間變化過程Fig.2 Spatial changes of impervious areas in the Wenyu River basin during the study periods

2.2 研究方法

為了對(duì)溫榆河夏季不同概率洪水的非一致性進(jìn)行診斷,本文選用Pettitt突變檢驗(yàn)法對(duì)溫榆河夏季不同概率流量數(shù)據(jù)進(jìn)行突變檢驗(yàn),Pettitt突變檢驗(yàn)?zāi)軌蜉^好地識(shí)別出數(shù)據(jù)序列的突變點(diǎn),其優(yōu)點(diǎn)為計(jì)算簡便,可以明確突變的時(shí)間[29]；采用GAMLSS模型對(duì)研究區(qū)不同概率洪水建立模型；最后采用皮爾遜相關(guān)系數(shù)法評(píng)估GAMLSS模型模擬效果及檢驗(yàn)優(yōu)選模型的魯棒性。

2.3 GAMLSS模型簡介

GAMLSS是一種引入位置、尺度和形狀的廣義參數(shù)可加模型,該模型于2005年由Rigby和Stasinopoulos提出[30],GAMLSS模型是假設(shè)因變量服從一種參數(shù)分布的廣義回歸模型,這種分布的所有參數(shù)可以通過解釋變量的函數(shù)來估計(jì),因此,GAMLSS模型的主要特點(diǎn)是根據(jù)解釋變量的值允許不同因變量的分布形式[30]。近年來,國內(nèi)外學(xué)者將GAMLSS模型應(yīng)用到水文分析中,該模型不僅能夠?qū)Ψ且恢滦缘乃臅r(shí)間序列進(jìn)行模擬分析,并且能夠?qū)Ψ且恢滦赃M(jìn)行成因分析。

(1)

式中:θk和ηk為長度為n的向量；βk=(β1k,…,βJkk)T為長度為Jk的回歸參數(shù)向量；Xk為一個(gè)已知的n×Jk的解釋變量矩陣；hjk為解釋變量的函數(shù),可以較靈活地采用協(xié)變量通過所選的分布函數(shù)來描述。

2.4 GAMLSS模型建立與優(yōu)選

表1 GAMLSS分布函數(shù)表達(dá)式及連接函數(shù)形式

表2 GAMLSS模型表達(dá)式及協(xié)變量因子

為了選取各概率流量的最優(yōu)擬合模型,引入廣義AIC(Akaike Information Criterion)準(zhǔn)則GAIC(GAIC,Generalized Akaike Information Criterion)作為評(píng)判標(biāo)準(zhǔn)[34]:

GAIC=DG+ωdf

(2)

式中:DG為模型的全局偏差；df為模型自由度；ω為懲罰因子。當(dāng)ω=2時(shí)稱為AIC準(zhǔn)則;當(dāng)ω=lgk時(shí)稱為SBC(Schwarz Bayesian Criterion)準(zhǔn)則,k為解釋變量的個(gè)數(shù)；AIC和SBC準(zhǔn)則是GAIC準(zhǔn)則的2種特例,取GAIC最小的模型作為最優(yōu)模型,本文選用SBC準(zhǔn)則作為評(píng)判準(zhǔn)則。

由于SBC準(zhǔn)則不能較直觀地反映出優(yōu)選模型的擬合效果,GAMLSS模型中worm圖可反映出優(yōu)選模型的殘差分布,可作為模型的評(píng)判標(biāo)準(zhǔn)。為了評(píng)估優(yōu)選模型的魯棒性,本研究采用留一法對(duì)最優(yōu)模型的預(yù)報(bào)流量值進(jìn)行交叉驗(yàn)證,即為將各個(gè)概率各個(gè)年份的流量逐一作為預(yù)報(bào)項(xiàng),將其余年份的驅(qū)動(dòng)因子作為預(yù)報(bào)因子,采用最優(yōu)模型來預(yù)報(bào)去除年份的流量值,逐一進(jìn)行這一過程,直到得到所有的流量值,采用Pearson相關(guān)系數(shù)法,分別計(jì)算觀測流量值與模型交叉驗(yàn)證預(yù)報(bào)值之間的相關(guān)系數(shù),從而驗(yàn)證模型對(duì)不同概率流量值的擬合效果。

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

皮爾遜(Pearson)相關(guān)系數(shù)是描述2個(gè)隨機(jī)變量線性相關(guān)的統(tǒng)計(jì)量[35],取值范圍為(-1,1),假設(shè)2個(gè)變量x(x1,x2,…,xn)和y(y1,y2,…,yn),2個(gè)變量之間的Pearson相關(guān)系數(shù)(r)計(jì)算公式式(3)所示:

(3)

當(dāng)r>0時(shí),表明2個(gè)變量呈正相關(guān),當(dāng)r越接近1表明正相關(guān)越顯著；當(dāng)r<0時(shí),表明2個(gè)變量呈負(fù)相關(guān),當(dāng)r越接近-1表明負(fù)相關(guān)越顯著；當(dāng)r=0時(shí),表明2變量相互獨(dú)立。

3 結(jié)果與討論

3.1 Pettitt test徑流突變點(diǎn)分析

采用Pettitt法對(duì)溫榆河夏季不同概率流量數(shù)據(jù)進(jìn)行突變檢驗(yàn),由于篇幅所限,本文僅列出95%、75%、50%、25%、10%和5%等6個(gè)概率流量的Pettitt突變點(diǎn)檢驗(yàn)結(jié)果(如圖3所示)。從圖中可以看出,不同概率的流量均存在突變點(diǎn),95%和75%概率的流量突變點(diǎn)發(fā)生在2011年,50%概率流量突變點(diǎn)發(fā)生在1996年,25%、10%和5%概率流量序列突變點(diǎn)均發(fā)生在1998年。

圖3 溫榆河夏季不同概率洪水Pettitt突變點(diǎn)檢驗(yàn)結(jié)果Fig.3 Pettitt test results of different probability discharges in summer in the Wenyu River basin

3.2 GAMLSS模型優(yōu)選

本研究對(duì)溫榆河夏季不同概率(僅列出部分概率模擬圖)的流量數(shù)據(jù)建立GAMLSS模型進(jìn)行擬合分析,GAMLSS模型可對(duì)預(yù)測變量的概率分布進(jìn)行擬合,換言之,GAMLSS模型擬合的不是1個(gè)單一值,而是全概率分布,圖4為各概率流量最優(yōu)模型擬合的全概率分布(主要計(jì)算95%、75%、50%、35%、20%和5%分位數(shù))。從圖4中可以看出,溫榆河夏季高頻流量呈現(xiàn)出逐漸上升的趨勢(shì),而中高頻流量比較穩(wěn)定,從圖中還可以看出,絕大多數(shù)觀測流量點(diǎn)據(jù)位于模擬值范圍內(nèi),說明優(yōu)選模型能夠較好地捕捉到觀測流量的變化特征。

圖4 不同概率洪水最優(yōu)模型擬合全概率分布Fig.4 Probabilistic distribution of selected best models for different probability of floods

根據(jù)SBC準(zhǔn)則,將SBC值最小的模型作為最優(yōu)模型,表3給出了各個(gè)概率優(yōu)選模型殘差分布的平均值、方差、偏態(tài)系數(shù)、峰度系數(shù)和Filliben相關(guān)系數(shù),全局偏差和SBC值,殘差分布的平均值越接近0,方差越接近1,偏度系數(shù)越接近0,峰度系數(shù)越接近3,Filliben相關(guān)系數(shù)都大于0.95,證明所選模型的殘差分布符合正態(tài)分布[35]。由于SBC準(zhǔn)則無法直觀反映出優(yōu)選模型的擬合效果,采用worm圖來反映GAMLSS模型的擬合效果,圖5給出了各個(gè)分位數(shù)優(yōu)選模型的殘差分布worm圖,圖中模型擬合殘差值位于兩條黑色曲線包圍的區(qū)域,代表模型擬合效果較好,從worm圖可以看出,各個(gè)概率的流量擬合殘差滿足要求。為了驗(yàn)證優(yōu)選模型的魯棒性,采用留一法對(duì)模型擬合結(jié)果進(jìn)行驗(yàn)證,計(jì)算觀測值與預(yù)報(bào)值序列的皮爾遜相關(guān)系數(shù)(如表3所示),經(jīng)計(jì)算,預(yù)測值與觀測值之間的r平均值為0.74,說明所選最優(yōu)模型的模擬效果較好。

表3 優(yōu)選模型擬合效果相關(guān)參數(shù)及相關(guān)系數(shù)

圖5 不同概率洪水模擬殘差worm圖Fig.5 Worm plots for different probability of floods

3.3 協(xié)變量分析

由于GAMLSS模型位置參數(shù)(μ)為反映預(yù)測變量量級(jí)變化的參數(shù),為了識(shí)別溫榆河流域夏季不同概率流量非一致性的關(guān)鍵驅(qū)動(dòng)要素,分析不同概率洪水非一致性的成因,探討該流域全概率洪水演變的驅(qū)動(dòng)機(jī)制,本文對(duì)μ表達(dá)式的系數(shù)進(jìn)行討論分析。圖6為不同概率洪水優(yōu)選模型μ表達(dá)式的系數(shù)變化圖,從圖中可以看出,μ的截距(如圖6(a))隨著洪水發(fā)生概率的降低而增高,這與前人研究結(jié)果一致[12,33]；從圖6(b)可以看出,發(fā)生概率高于70%的夏季小洪水主要受到城市不透水面積比變化的影響,其不透水面積變化比的系數(shù)要遠(yuǎn)大于降水的系數(shù)；而發(fā)生概率為45%到70%的夏季洪水受到降水和不透水面積比變化的雙重影響,但降水的系數(shù)要大于不透水面積比的系數(shù),說明這一發(fā)生概率階段的洪水,相比于降水,受到不透水面積比變化的影響較小；對(duì)于低于概率45%的夏季洪水,不透水面積比變化不再選為關(guān)鍵驅(qū)動(dòng)要素,降雨的系數(shù)越來越大,說明夏季低于45%概率洪水演變的主要驅(qū)動(dòng)因素為降水的變化。

圖6 夏季不同概率洪水優(yōu)選模型μ截距及協(xié)變量系數(shù)變化Fig.6 Intercept and the coefficients of covariates of location parameters μ of selected models for different probability floods in summer

綜上所述,溫榆河夏季中小洪水(中高頻洪水)的變化主要受到城市化的影響,而溫榆河在過去幾十年間城市化發(fā)展十分迅速,這也就較好地詮釋了溫榆河夏季中小洪水呈現(xiàn)出逐漸上升的趨勢(shì)(如圖4)；而夏季大洪水(低頻流量)變化的主要驅(qū)動(dòng)因素為降水的變化。

4 結(jié) 論

由于受到氣候變化和高強(qiáng)度人類活動(dòng)的影響,城市洪水在長時(shí)間序列上表現(xiàn)出非一致性特性,本文采用GAMLSS模型對(duì)溫榆河夏季全概率洪水變化的主要驅(qū)動(dòng)機(jī)制進(jìn)行分析探討。主要結(jié)論如下：

(1) 溫榆河流域夏季不同概率的洪水在研究期呈現(xiàn)出非一致性,并且不同等級(jí)的洪水出現(xiàn)突變點(diǎn)的年份不一致；GAMLSS模型能夠較好地應(yīng)用到城市洪水分析中,該模型可以較好地捕捉到城市洪水的變化特征,優(yōu)選模型預(yù)測的中位數(shù)值與觀測流量值之間的Pearson相關(guān)系數(shù)平均值為0.74,說明模型擬合效果較好。

(2) 降水和不透水面積比的變化是溫榆河夏季洪水演變的主要驅(qū)動(dòng)要素,不同等級(jí)洪水呈現(xiàn)出的非一致性特性具有不同的驅(qū)動(dòng)成因。溫榆河夏季高于概率70%的小洪水變化的主要驅(qū)動(dòng)要素為不透水面積比變化；概率為65%到50%之間的中小洪水變化的主要驅(qū)動(dòng)要素為降水和不透水面積比變化,其中受到降水的影響要顯著大于城市不透水面積變化的影響；而對(duì)于低于45%概率的大洪水,城市不透水面積變化比不再作為主要驅(qū)動(dòng)要素,而其主要受到降水變化的影響。

(3) 下墊面類型改變的影響主要體現(xiàn)在中小洪水上,而大洪水及極端洪水的演變主要受到降水變化的影響,這一結(jié)論可對(duì)中國海綿城市建設(shè)提供一定的科學(xué)借鑒價(jià)值。

(4) GAMLSS模型可以靈活地選取與洪水變化相關(guān)的解釋變量來描述洪水序列的非一致性,并且能夠?qū)Ψ且恢滦赃M(jìn)行歸因分析,但GAMLSS模型更擅長于評(píng)估水文資料在時(shí)間系列上的非一致性成因,而在對(duì)水文序列的預(yù)報(bào)方面可能存在一定的缺陷。