基于Budyko假設(shè)和微分方程的淮河上中游徑流變化歸因分析

葉 婷，石 朋，鐘 華，瞿思敏，吳洪石，沈玲怡

(1.河海大學(xué)水文水資源學(xué)院，江蘇 南京 210098； 2.南京水利科學(xué)研究院水文水資源研究院，江蘇 南京 210098)

徑流特征變化受到多方面因素的綜合影響，包括氣候變化、土地利用/覆蓋變化和人類取用水等。徑流變化的主要驅(qū)動(dòng)因素被廣泛歸因于氣候變化和人類活動(dòng)[1-4]。然而，由于氣候變化和人類活動(dòng)之間復(fù)雜的相互作用，水文界在解釋其對(duì)徑流變化驅(qū)動(dòng)機(jī)制方面仍面臨巨大挑戰(zhàn)[5-6]。目前，評(píng)估徑流變化的主要方法有微觀角度的牛頓法和宏觀角度的達(dá)爾文法[7]。多數(shù)水文模型(如SWAT模型、VIC模型、TOPMODEL等)被歸類為牛頓法，并被廣泛應(yīng)用于國內(nèi)眾多流域。白琪階等[8]運(yùn)用SWAT模型研究發(fā)現(xiàn)漳衛(wèi)南運(yùn)河流域人類活動(dòng)引起的徑流變化是自然因素影響的4倍；鄧鵬等[9]通過VIC模型研究發(fā)現(xiàn)，未來氣候變化使得鄱陽湖汛期徑流增加，枯季徑流減小。雖然水文模型方法能較好地定量評(píng)估徑流變化，但建模過程復(fù)雜且需要大量數(shù)據(jù)。而達(dá)爾文法通過系統(tǒng)思維解釋水文行為，計(jì)算過程簡潔，在徑流變化歸因分析研究中也有較為廣泛的應(yīng)用，基于Budyko假設(shè)的方法是其代表性方法之一，該方法來源于Budyko假設(shè)和水平衡理論，已在人類活動(dòng)劇烈的地區(qū)得到廣泛應(yīng)用[10-11]。例如：張麗梅等[12]利用Budyko假設(shè)研究發(fā)現(xiàn)，劇烈的人類活動(dòng)如工業(yè)用水以及退耕還林等引起的下墊面變化是渭河流域徑流量減少的主要原因；王衛(wèi)光等[13]證實(shí)了Budyko水熱耦合平衡方程在黃河流域的適用性；Xin等[14]利用Budyko假設(shè)研究發(fā)現(xiàn)氣候變化和直接人類活動(dòng)都會(huì)導(dǎo)致輝發(fā)河年尺度流量的減少，其中68%的變化是由人類活動(dòng)引起的。

淮河流經(jīng)豫、皖、蘇三省，是我國重要的商品糧基地之一。當(dāng)前基于Budyko假設(shè)的淮河徑流變化歸因研究多集中于徑流年際間的變化，忽略了年內(nèi)其他關(guān)鍵徑流特征[15-17]。基于此，本文選取淮河上中游為研究區(qū)域，利用多重對(duì)數(shù)轉(zhuǎn)換線性回歸方法擬合各徑流特征指標(biāo)與年徑流深、氣候指標(biāo)的相關(guān)關(guān)系，結(jié)合基于Budyko假設(shè)的彈性系數(shù)法建立評(píng)估框架，探究氣候變化和人類活動(dòng)對(duì)不同時(shí)間尺度徑流特征的影響，以期能為淮河水資源配置及河流治理提供參考。

1 研究區(qū)概況

淮河干流自西向東，經(jīng)河南省南部、安徽省中部，在江蘇省中部注入洪澤湖，經(jīng)洪澤湖調(diào)蓄后，主流經(jīng)江水道至揚(yáng)州三江營注入長江?；春由嫌蝺砂渡角鹌鸱?，支流眾多；中游地勢(shì)平緩，多湖泊洼地。上中游流域內(nèi)息縣、王家壩、魯臺(tái)子和蚌埠4個(gè)水文站的控制流域面積分別為10 190 km2、30 630 km2、88 630 km2和121 330 km2?；春由现杏瘟饔虻靥幈本?0°55′～34°56′、東經(jīng)111°55′～117°30′之間，其中息縣以上主河道長250 km，建有南灣和石山口2座大型水庫，多年平均年徑流量為37.34億m3。王家壩站位于淮河上中游分界處，干流右岸多為山區(qū)河流，急流迅速，左岸為平原區(qū)，地勢(shì)平坦，多年平均年徑流量為84.75億m3。魯臺(tái)子以上主河道長511 km，多年平均年徑流量為211.5億m3。魯蚌區(qū)間主要支流為茨淮新河和渦河，存在1個(gè)蓄洪區(qū)和7個(gè)行洪區(qū)，蚌埠站多年平均徑流量為305億m3。

水文氣象數(shù)據(jù)包括淮河上中游息縣、王家壩、魯臺(tái)子和蚌埠4個(gè)水文站1955—2015年的長系列日徑流數(shù)據(jù)，信陽等14個(gè)氣象站(圖1)的逐日氣象數(shù)據(jù)。DEM數(shù)據(jù)采用地理空間數(shù)據(jù)云(http://www.gscloud.cn/)提供的SRTMDEM 90 m分辨率數(shù)據(jù)產(chǎn)品。土地利用數(shù)據(jù)采用中國科學(xué)院資源環(huán)境科學(xué)與數(shù)據(jù)中心(resdc.cn)提供的1 km空間分辨率產(chǎn)品，并根據(jù)中國現(xiàn)行一級(jí)分類標(biāo)準(zhǔn)對(duì)其進(jìn)行了重分類。

2 研 究 方 法

利用Mann-Kendall趨勢(shì)檢驗(yàn)(M-K檢驗(yàn))[18]、滑動(dòng)T檢驗(yàn)[19]識(shí)別徑流特征突變年份，并依此劃分基準(zhǔn)期和變化期。應(yīng)用多重對(duì)數(shù)線性回歸方法(MLLR)[20-21]研究各徑流特征指標(biāo)與年徑流深、氣候指標(biāo)之間的相關(guān)關(guān)系[22]。采用基于Budyko假設(shè)的經(jīng)驗(yàn)公式將年徑流深變化歸因于氣候變化和人類活動(dòng)。結(jié)合基于Budyko假設(shè)的彈性系數(shù)和徑流特征指標(biāo)微分方程，求解由于年徑流深變化和氣候指標(biāo)變化引起的徑流特征指標(biāo)變化，其余變化則被看作人類活動(dòng)的影響。最后建立土地利用轉(zhuǎn)移矩陣，分析人類活動(dòng)對(duì)徑流特征變化的驅(qū)動(dòng)機(jī)制。

2.1 徑流特征和氣候指標(biāo)

徑流特征簡明地描述了長期觀測(cè)到的徑流狀態(tài)，包括反映徑流分布以及徑流動(dòng)態(tài)的指數(shù)。研究選取了7個(gè)流量指標(biāo)表征淮河上中游的徑流特征，分別為3個(gè)平均流量指標(biāo)(年徑流深MGAF、汛期徑流深MFF和非汛期徑流深MNF)、2個(gè)低流量指標(biāo)(年最小日徑流深A(yù)min和90%低流量徑流深L90)、2個(gè)高流量指標(biāo)(年最大日徑流深A(yù)max和10%高流量徑流深H10)。采用11個(gè)氣候指標(biāo)表征研究流域的氣候變化，分別是6個(gè)降水指標(biāo)(年降水量P、最大7 d降水量P7d、汛期降水量Pf、非汛期降水量Pnf、無降雨天數(shù)Pz和非汛期與汛期降水量比值Psea)、5個(gè)潛在蒸散發(fā)量指標(biāo)(年潛在蒸散發(fā)量E0、7 d最大潛在蒸散發(fā)量E7d、汛期潛在蒸散發(fā)量Ef、非汛期潛在蒸散發(fā)量Enf和非汛期與汛期潛在蒸散發(fā)量比值Esea)。采用逐步回歸方法選取每個(gè)徑流特征的3個(gè)最顯著氣候指標(biāo)。流域的潛在蒸散發(fā)根據(jù)FAO修正的Penman-Monteith公式[23]進(jìn)行計(jì)算。

2.2 氣候變化和人類活動(dòng)的影響評(píng)估

2.2.1 自然條件下的徑流特征

通過M-K檢驗(yàn)和滑動(dòng)T檢驗(yàn)，將所有徑流特征時(shí)間序列分為基準(zhǔn)期和變化期。為了保持同一站點(diǎn)前后周期的一致性，將各指標(biāo)序列的平均突變年份作為站點(diǎn)徑流變化的突變年份[22]。選擇年徑流深和氣候指標(biāo)作為回歸建模的自變量，各個(gè)徑流特征作為因變量，采用MLLR方法擬合前期徑流特征：

(1)

式中：b0、b1、…、bN為MLLR的系數(shù)；MAF為年徑流深；bR為年徑流深的MLLR系數(shù)；S為徑流特征；C為氣候指標(biāo)。

根據(jù)MLLR方法進(jìn)行轉(zhuǎn)換，式(1)變?yōu)?/p>

(2)

由于淮河上中游的站點(diǎn)可能出現(xiàn)零徑流日，而式(2)并不適用此情景，因此在考慮零徑流日的低流量指標(biāo)時(shí)，使用3年平均值建立方程。

2.2.2 基于Budyko假設(shè)的方法

Budyko假設(shè)揭示了年徑流深與其影響因子(降雨、潛在蒸散發(fā)和流域下墊面特征)之間的非線性關(guān)系。傅抱璞[24]在Budyko假設(shè)的基礎(chǔ)上推導(dǎo)出具有普適性的年尺度流域水熱耦合平衡方程，隨后Yang等[25]改進(jìn)該方程，提出水熱耦合平衡方程的經(jīng)驗(yàn)公式。基于此經(jīng)驗(yàn)公式，忽略流域長期蓄水量變化，則流域多年水量平衡方程為

(3)

其中n=ω-0.72

式中：E為年實(shí)際蒸散發(fā)量；n為無量綱的下墊面參數(shù)；ω為反映流域下墊面特征的無量綱數(shù)[26]，如土壤屬性、坡度、土地利用和植被覆蓋和氣候季節(jié)性，取值范圍為[1,+∞]；E0/P為干燥度指數(shù)。

研究表明，ω可以很好地表征人類活動(dòng)[12,17,20]，ω越大，流域平均蒸散發(fā)越大且徑流深越小?；鶞?zhǔn)期的n可以根據(jù)基準(zhǔn)期多年平均年徑流深、降水量和蒸散發(fā)量求得。假設(shè)P、E0和n是獨(dú)立變量，由氣候變化和人類活動(dòng)導(dǎo)致的徑流深變化可表示為徑流深對(duì)該變量的偏導(dǎo)數(shù)與其本身變化量的乘積，根據(jù)彈性系數(shù)法[12]可得

(4)

設(shè)E0/P為φ，可得：

(5)

2.2.3 微分方程的推導(dǎo)與求解

流域?qū)傩栽诤艽蟪潭壬蠜Q定了徑流特征，通過建立基準(zhǔn)期流域?qū)傩院蛷搅魈卣髦g的多重對(duì)數(shù)線性回歸關(guān)系，可以定量分析人類活動(dòng)和氣候變化對(duì)徑流變化的貢獻(xiàn)率。該方法已被廣泛應(yīng)用于徑流特征的預(yù)測(cè)分析[20]。根據(jù)式(2)可知，徑流特征變化的新微分方程可以推導(dǎo)為

(6)

根據(jù)式(4)和式(6)可知，氣候變化對(duì)徑流特征變化的影響(ΔSC)可以表示為

(7)

式中?MAF/?P和?MAF/?E0可以根據(jù)式(5)計(jì)算。

假設(shè)人類活動(dòng)的影響與氣候變化無關(guān)，人類活動(dòng)對(duì)徑流特征變化的影響(ΔSH)可以表示為

ΔSH=ΔS-ΔSC

(8)

因此，在徑流成因分析中，氣候變化以及人類活動(dòng)對(duì)徑流特征變化影響的貢獻(xiàn)率PC和PH分別為為

(9)

3 結(jié) 果 分 析

3.1 徑流突變分析

采用M-K檢驗(yàn)、滑動(dòng)T檢驗(yàn)方法確定淮河上中游4個(gè)水文站徑流突變年份，得到息縣、王家壩、魯臺(tái)子和蚌埠站的徑流突變年份分別為1999年、2008年、1988年和1991年。由表1可知，4個(gè)水文站的多年平均年徑流深在突變年份后均呈下降趨勢(shì)，這與孫鵬等[27]的研究結(jié)論一致。其中息縣站變化期多年平均年徑流深相較于基準(zhǔn)期下降了6.25%，而王家壩、魯臺(tái)子和蚌埠站分別下降了39.67%、13.54%和20.16%，中游多年平均年徑流深下降趨勢(shì)較上游明顯，主要是因?yàn)橹杏蔚貐^(qū)受到的人類活動(dòng)影響更多。

表1 4個(gè)水文站徑流特征、年降水量和潛在蒸散量的平均值及其變化率

除息縣站外，其下游水文站的MFF下降幅度均大于MNF，而王家壩站的下降幅度更是達(dá)到了44.25%。王家壩站位于淮河上中游交界處，其下游的濛洼蓄洪區(qū)1954—2015年累計(jì)12個(gè)年份15次滯蓄洪水，蓄洪總量達(dá)75億m3。重大水利工程汛期防洪調(diào)控不僅大幅度減少了下游汛期徑流，同時(shí)確保了沿淮大中城市的防洪安全。

此外，變化期的高、低流量指標(biāo)與基準(zhǔn)期相比，息縣站和王家壩站大部分都呈上升趨勢(shì)(除王家壩站H10下降47.59%)，而魯臺(tái)子站和蚌埠站大部分呈下降趨勢(shì)(除魯臺(tái)子站L90無變化)。這也是上游水利工程對(duì)徑流的調(diào)控所致。

3.2 徑流特征的估計(jì)

由表2可知，高流量指標(biāo)、汛期徑流深與年徑流深以及氣候指標(biāo)的相關(guān)性最高，相關(guān)系數(shù)范圍為0.859～0.987；而低流量指標(biāo)、非汛期徑流深用年徑流和氣候指標(biāo)估計(jì)較為困難，相關(guān)性偏低，大部分指標(biāo)方程相關(guān)系數(shù)小于0.8(除魯臺(tái)子站L90徑流特征方程相關(guān)系數(shù)為0.843)。

表2 基準(zhǔn)期徑流特征方程以及相關(guān)系數(shù)

息縣、王家壩、魯臺(tái)子和蚌埠4個(gè)水文站徑流特征多重相關(guān)系數(shù)變化范圍分別為0.467Amin～0.967MFF、0.612Amin～0.967MFF、0.757Amin～0.987H10和0.541MNF～0.984MFF?？傮w而言，各站點(diǎn)的平均相關(guān)系數(shù)均達(dá)到0.75以上，徑流特征與氣候變量呈現(xiàn)出明顯的非線性關(guān)系。表明本文方法能夠很好地應(yīng)用于徑流特征估計(jì)。

3.3 氣候變化和人類活動(dòng)的影響分析

根據(jù)Budyko曲線擬合可知，息縣、王家壩、魯臺(tái)子和蚌埠站的n值分別為1.883、2.118、1.827和2.141，干燥度指數(shù)分別為0.928、1.116、0.997和1.105，由n和干燥度指數(shù)的數(shù)值變化可知，隨著干旱程度的增加，4個(gè)水文站的徑流逐漸減小。在息縣、王家壩、魯臺(tái)子和蚌埠4個(gè)水文站控制的流域中，氣候變化對(duì)徑流特征的影響貢獻(xiàn)率變化范圍分別為1.54%～12.45%、11.50%～19.09%、7.58%～19.58%和0.01%～16.77%。人類活動(dòng)對(duì)徑流特征的影響貢獻(xiàn)率變化范圍分別為87.55%～98.46%、80.91%～88.50%、80.42%～92.42%和83.23%～99.99%，貢獻(xiàn)率均超過80%。人類活動(dòng)在淮河上中游的徑流特征變化中起主導(dǎo)作用。王國慶等[28]對(duì)中國主要江河徑流變化成因進(jìn)行定量分析，發(fā)現(xiàn)人類活動(dòng)對(duì)包括淮河在內(nèi)的大多數(shù)河流徑流變化的影響率超過50%，是徑流變化的主要驅(qū)動(dòng)因素。城市化進(jìn)程的推進(jìn)與“三生”用水需求的不斷擴(kuò)大，都可能使徑流呈現(xiàn)下降趨勢(shì)。

對(duì)于淮河流域而言，土地利用變化是徑流變化的關(guān)鍵影響因素，本文將其納入人類活動(dòng)的影響[16-17]。息縣、王家壩、魯臺(tái)子和蚌埠站的徑流特征平均突變年份分別為1999年、2008年、1988年和1991年，依照各站點(diǎn)突變年份，構(gòu)建1980—1990年(魯臺(tái)子和蚌埠站)、1980—1995年(息縣站)和1980—2005年(王家壩站)3個(gè)時(shí)期土地利用轉(zhuǎn)移矩陣(表3)，從而分析不同站點(diǎn)基準(zhǔn)期和變化期土地利用類型空間分布變化。

表3 淮河流域上中游1980—1990、1980—1995和1980—2005年土地利用轉(zhuǎn)移矩陣 單位：km2

1980—1990年土地利用類型變化較大的有：耕地向水域及水利設(shè)施用地、住宅用地轉(zhuǎn)變，變化面積分別為218 km2和208 km2。1980—1995年土地利用類型變化較大的有：耕地向住宅用地轉(zhuǎn)變，變化面積為865 km2；其次，耕地向水域及水利設(shè)施用地轉(zhuǎn)變、林地和草地向耕地轉(zhuǎn)變的面積都大于200 km2；而水域及水利設(shè)施用地、住宅用地向耕地轉(zhuǎn)變的面積均大于100 km2。1980—2005年土地利用類型變化較大的有：耕地向住宅用地轉(zhuǎn)變，轉(zhuǎn)變面積高達(dá)1 391 km2；其次，林地和草地向耕地、草地向林地轉(zhuǎn)變的面積均大于200 km2；水域及水利設(shè)施用地向耕地轉(zhuǎn)變面積也達(dá)到了100 km2。20世紀(jì)末以及21世紀(jì)初淮河上中游住宅用地大幅度增加，城市化進(jìn)程加快，農(nóng)業(yè)、工業(yè)及生活用水增加使得河流徑流減少。劉曉麗等[17]研究下墊面變化對(duì)淮河流域中上游徑流變化的影響中也得到了相同的結(jié)論。

4 結(jié) 論

a.大多數(shù)徑流特征序列都檢測(cè)到了基準(zhǔn)期與變化期，突變年份主要出現(xiàn)在20世紀(jì)90年代與21世紀(jì)初。與基準(zhǔn)期相比，4個(gè)水文站變化期的大部分徑流特征都呈現(xiàn)下降趨勢(shì)，尤其是平均流量指標(biāo)，下降幅度從2.25%到44.25%不等。

b.氣候指標(biāo)對(duì)徑流特征有著重要影響。其中，高流量指標(biāo)、汛期徑流深與年徑流深以及氣候指標(biāo)的相關(guān)性最高，相關(guān)系數(shù)范圍為0.859～0.987；低流量指標(biāo)、非汛期徑流深與年徑流深以及氣候指標(biāo)的相關(guān)性相對(duì)較低，相關(guān)系數(shù)小于0.8(除魯臺(tái)子站L90相關(guān)系數(shù)為0.843)。

c.淮河上中游氣候變化和人類活動(dòng)對(duì)大部分徑流特征起到負(fù)面作用，其中人類活動(dòng)為主導(dǎo)因素。在息縣、王家壩、魯臺(tái)子和蚌埠4個(gè)水文站控制流域，人類活動(dòng)貢獻(xiàn)率范圍分別為87.55%～98.46%、80.91%～88.50%、80.42%～92.42%和83.23%～99.99%，4個(gè)站點(diǎn)控制流域貢獻(xiàn)率都達(dá)到80%以上。