膠東半島典型流域地表水文過程變化及其歸因分析

王宏宙，孫善磊，王 潔，于志國，閆桂霞，周舒佳，畢早瑩

（1.南京信息工程大學(xué)氣象災(zāi)害教育部重點(diǎn)實(shí)驗(yàn)室，南京210044；2.常州市金壇區(qū)氣象局，江蘇常州213200）

0 引 言

隨著社會(huì)經(jīng)濟(jì)的快速發(fā)展和人口的不斷增長(zhǎng)，人類對(duì)水資源的需求增大，水資源短缺問題日益加劇，已成為影響社會(huì)經(jīng)濟(jì)可持續(xù)發(fā)展的重要約束條件［1］。作為水文循環(huán)的重要環(huán)節(jié)，地表水文過程始終與水資源的合理開發(fā)利用密切相關(guān)；因此，了解和認(rèn)識(shí)地表水文過程對(duì)全球變化的響應(yīng)已成為水文、氣候等科學(xué)研究的熱點(diǎn)科學(xué)問題［2-4］。

地表水文過程始終與氣候變化和人類活動(dòng)密切相關(guān)［5-8］。如降水直接決定了地表水分供給；溫度、風(fēng)速、相對(duì)濕度和輻射等氣候要素通過改變蒸散發(fā)（ET）過程，影響地表水分支出［6］。而人類活動(dòng)通過改變流域下墊面特征（如水利工程設(shè)施的修建、生態(tài)恢復(fù)、城市擴(kuò)張等）影響產(chǎn)匯流機(jī)制及地表水文過程［6-8］。在全球氣候變化加劇的背景下，人類改造自然環(huán)境的能力不斷提高，使得區(qū)域乃至全球水文循環(huán)變化變得愈加復(fù)雜。為厘清相關(guān)機(jī)制，學(xué)者們采用各種手段和方法，剝離了不同氣候要素變化和人類活動(dòng)對(duì)徑流、ET的影響，取得了大量研究成果［9-13］。常用的區(qū)分氣候變化和人類活動(dòng)影響的手段有：流域?qū)Ρ确治龇ǎ?］、敏感性分析［10］、陸面模式法［11］、水文模型法［12］及彈性系數(shù)法［13］等；但這些方法存在一定的局限性。比如，流域?qū)Ρ确治龇ê茈y找到兩個(gè)完全一致的流域，且費(fèi)時(shí)費(fèi)力［9］；而水文模型和陸面模式法對(duì)驅(qū)動(dòng)場(chǎng)要求較高，耗費(fèi)計(jì)算機(jī)時(shí)，尤其對(duì)大流域［11，12］；彈性系數(shù)法未考慮影響因子間相互作用對(duì)水文過程的影響，存在不確定性［10］。近來，Sun等［14］原創(chuàng)性地提出了一種基于敏感性試驗(yàn)的多控制因子聯(lián)立求解方法，可以有效而準(zhǔn)確地剝離各因子對(duì)水文分量變化的影響，已被成功應(yīng)用于鄱陽湖流域和中國西南地區(qū)。另外，自Budyko 框架理論提出以來，因其需要資料少、結(jié)構(gòu)簡(jiǎn)單、物理機(jī)制明確、模型參數(shù)少（僅1 個(gè)流域特征參數(shù)），已被廣泛應(yīng)用于水文循環(huán)變化研究［15，16］；更重要的是，Budyko 參數(shù)反映了流域特征，為綜合研究氣候變化和人類活動(dòng)對(duì)水文循環(huán)的影響提供了重要的模型工具。

膠東半島包括青島市、煙臺(tái)市、威海市及濰坊市的一部分，總面積42 680 km2，為我國北方主要的經(jīng)濟(jì)技術(shù)開發(fā)區(qū)。膠東半島區(qū)域水系多發(fā)源于中部山地，河川徑流洪枯期懸殊，徑流集中在汛期；由于特殊的地形地貌和持續(xù)性少雨，該區(qū)域已出現(xiàn)嚴(yán)重的水資源短缺問題，嚴(yán)重制約了膠東半島經(jīng)濟(jì)發(fā)展［17-23］。近年來，膠東半島地區(qū)氣候變化加劇，加之人類活動(dòng)規(guī)模增大、強(qiáng)度增強(qiáng)，造成了該地區(qū)水土流失嚴(yán)重、涵養(yǎng)水分能力降低等一系列環(huán)境問題［19，23］。目前對(duì)膠東半島流域水文循環(huán)變化的關(guān)注較少，尤其是ET變化尚不清楚；雖然少數(shù)學(xué)者對(duì)該地區(qū)徑流變化特征進(jìn)行了研究［17-23］，但依然缺乏關(guān)于水文循環(huán)變化機(jī)理的定量化歸因分析。因此，系統(tǒng)分析膠東半島氣候和水文過程變化規(guī)律，并揭示氣候變化和人類活動(dòng)在其中的作用，已成為理解本地區(qū)水文循環(huán)變化機(jī)制、改善水資源管理措施及合理開發(fā)利用的當(dāng)務(wù)之急。

綜上，擬選取膠東半島4個(gè)典型流域，詳細(xì)分析氣候和主要陸面水文分量（如ET和徑流）的變化特征，并采用Budyko 方程結(jié)合多控制因子聯(lián)立求解方法，定量化估算氣候變化和人類活動(dòng)（以Budyko 參數(shù)表示）對(duì)ET和徑流的影響，以此探討膠東半島地表水文循環(huán)變化機(jī)理。本研究有望為應(yīng)對(duì)膠東半島水資源短缺，建立科學(xué)的水資源管理制度，及維持社會(huì)經(jīng)濟(jì)健康發(fā)展提供重要的定量化信息。

1 資料與方法

1.1 資料及處理

氣象要素采用國家氣候中心提供的膠東半島1961-2010年間16個(gè)常規(guī)氣象站（圖1）的逐日降水（mm）、日照時(shí)數(shù)（h）、10 m風(fēng)速（m/s）、相對(duì)濕度（%），及最低（oC）、最高（oC）和平均溫度（oC）觀測(cè)資料。由于缺少輻射觀測(cè)，利用日照時(shí)數(shù)和Allen等［30］推薦的公式計(jì)算凈輻射。水文要素為1961-2010年4 個(gè)水文站（圖1）的逐年徑流數(shù)據(jù)（m3/s），來自煙臺(tái)市水文水資源勘測(cè)局和青島市水文水資源勘測(cè)局。采用1 km×1 km SRTM 數(shù)字高程模型（http：//westdc.westgis.ac.cn）和ArcGIS 10.2 水文分析工具，提取流域邊界和面積；然后，根據(jù)流域面積估算地表徑流（mm）?？紤]世界糧農(nóng)組織推薦的FAO56 Penman-Monteith 公式具有明確的物理含義，擬采用該公式估算各氣象站點(diǎn)PET，詳見Allen等［24］。采用流域水量平衡方法估算流域不同年代ET；首先，將1961-2010年分為1961-1980年（基準(zhǔn)期）、1981-1990年（1980s）、1991-2000年（1990s）和2001-2010年（2000s）4 個(gè)時(shí)段，采用反距離權(quán)重法對(duì)年降水（P）進(jìn)行空間插值，獲得各流域不同時(shí)段年降水；以P減去徑流獲得ET。氣候和水文分量年代變化，以各時(shí)段減去基準(zhǔn)期表示。

圖1 流域、氣象站和水文站分布Fig.1 Locations of the selected 4 basins，weather and hydrological sites

1.2 Budyko方程

由Budyko 假設(shè)，ET受水分供給（即降水）和能量供給（即PET）的共同影響；流域尺度上，降水和PET之間存在耦合平衡關(guān)系。在此基礎(chǔ)上，傅抱璞［25］和Yang 等［26］基于無量綱分析和數(shù)理方程，提出了包括流域特征影響的Budyko模型。本研究采用Yang等［26］的方程，

一般，n越大，ET消耗降水越多，反之相反。利用降水、PET和ET，分別擬合基準(zhǔn)期、1980s、1990s 和2000s 參數(shù)n，即n基準(zhǔn)期、n1980s、n1990s和n2000s。關(guān)于徑流（Q）的方程可寫為，

1.3 歸因方法

為減少影響因素間相互作用帶來的不確定性，擬參考Sun等［14］提出的基于敏感性試驗(yàn)的多控制因子聯(lián)立求解方法，設(shè)計(jì)試驗(yàn)，分離各因子對(duì)ET和徑流變化的貢獻(xiàn)。具體如下：首先，由影響ET或徑流的氣候因子（溫度、凈輻射、相對(duì)濕度、風(fēng)速和降水）和參數(shù)n，設(shè)計(jì)7 組試驗(yàn)，其中1 組基準(zhǔn)試驗(yàn)（EXP-BASE）和6 組敏感性試驗(yàn)（EXP_non-Y，Y為影響因子），具體見表1。以1980s的溫度試驗(yàn)（EXP_non-Tave）為例，溫度設(shè)置為基準(zhǔn)期，凈輻射、相對(duì)濕度、風(fēng)速維持在1980s，采用FAO56 Penman-Monteith 公式計(jì)算PET，然后，將計(jì)算的PET、1980s 降水和n1980s代入方程（1）或（2），估算EXP_non-Tave 的ET或徑流。最后，聯(lián)立方程，獲得各因素對(duì)水文分量變化的貢獻(xiàn)。以EXP_non-Y為例，其對(duì)應(yīng)的ET或徑流（相對(duì)于基準(zhǔn)期）變化可認(rèn)為是除Y以外的其余因子共同引起的，可表示為：

表1 試驗(yàn)設(shè)計(jì)Tab.1 Configurations of experiments

需要指出的，PET對(duì)ET或徑流變化的貢獻(xiàn)為溫度、輻射、相對(duì)濕度和風(fēng)速貢獻(xiàn)之和。為確定ET、徑流變化的主控因子，首先，判斷參數(shù)n的貢獻(xiàn)是否大于氣候變化（即降水和PET貢獻(xiàn)之和）；如果大于，則ET或徑流變化的主控因子為參數(shù)n；相反，氣候變化影響較大且降水貢獻(xiàn)大于（小于）PET，則主控因子為降水（PET）。

2 結(jié)果分析

2.1 降水、PET及參數(shù)n變化特征

由圖2（a），各流域基準(zhǔn)期降水基本在700 mm 左右，以團(tuán)旺和臧格莊較大（＞760 mm）；PET均在1 000 mm以上，且臧格莊最大（1 118.40 mm）；產(chǎn)芝、團(tuán)旺、尹府和臧格莊的基準(zhǔn)期參數(shù)n分別為1.40、1.34、1.11 和1.42，說明產(chǎn)芝和臧格莊有較多的降水通過ET的形式被耗散，而尹府有較少的降水被ET耗散。圖2（b）為降水較基準(zhǔn)期的變化，總體而言，各流域不同時(shí)段降水均減小，特別是1980s 和1990s（除尹府）減小較大（＜-100 mm）；2000s 降水依然偏少，但減小幅度基本在50 mm 以內(nèi)。由圖2（c），PET變化表現(xiàn)出一定的區(qū)域性差異；在產(chǎn)芝、團(tuán)旺和尹府，1980s 和1990s 的PET增加，但2000s 下降，變化幅度均小于30 mm；臧格莊各時(shí)段PET均減小，且1990s 和2000s 減小明顯（＜-55 mm）。較基準(zhǔn)期，各流域參數(shù)n在各時(shí)段一致性地增加［圖2（d）］，說明由于人類活動(dòng)的增強(qiáng)，ET分配的降水增加；除個(gè)別時(shí)段外，參數(shù)n的增幅均大于0.5，尤以產(chǎn)芝1980s 和2000s增幅最大，在1.0以上。

圖2 不同時(shí)期各流域主要水文氣候分量和參數(shù)n及降水、PET、參數(shù)n較基準(zhǔn)期的變化Fig.2 Values of the major hydro-meteorological elements and parameter n over the 4 basins during each period and changes in precipitation，PET and parameter n relative to the baseline period

2.2 歸因方法評(píng)估

在進(jìn)行歸因分析之前有必要對(duì)剝離因子貢獻(xiàn)的方法進(jìn)行評(píng)估，本研究采用各因子貢獻(xiàn)之和與觀測(cè)的ET、徑流變化散點(diǎn)圖，以及相關(guān)系數(shù)（R）、趨勢(shì)系數(shù)和均方根誤差（RMSE）等指標(biāo)，結(jié)果如圖3。很明顯，該方法剝離的各因子貢獻(xiàn)之和與觀測(cè)的ET、徑流變化幾乎處在1∶1 線上。由定量化指標(biāo)，各因子貢獻(xiàn)之和與觀測(cè)的ET、徑流變化的R幾乎為1，RMSE為5.12 mm，擬合趨勢(shì)線的趨勢(shì)系數(shù)也近乎為1；均說明該剝離因子貢獻(xiàn)方法可以較為準(zhǔn)確地估算各影響因子對(duì)ET、徑流變化的貢獻(xiàn)。

圖3 歸因方法評(píng)估Fig.3 Validation for the selected attribution method

2.3 蒸散發(fā)變化及其歸因

如圖2（a），尹府流域基準(zhǔn)期ET為454.98 mm，其他流域均在500 mm 左右，以臧格莊最大，為550.53 mm。由圖4（a），就產(chǎn)芝、團(tuán)旺和尹府，各時(shí)段ET均增加，且產(chǎn)芝和尹府2000s的增加最明顯（～100 mm）；相反，臧格莊不同時(shí)段ET均減小，且1990s最大，為-35.60 mm。圖4（b）～（d）為各流域降水和參數(shù)n對(duì)不同年代ET變化的貢獻(xiàn)。1980s，由于降水減少，產(chǎn)芝和尹府ET減小了～95 mm，而團(tuán)旺和臧格莊減小了～80 mm；隨著PET增大，產(chǎn)芝、團(tuán)旺和尹府ET增大（＞10 mm），其中尹府增加最大（20.51 mm），需要注意的，盡管臧格莊PET減小，但其貢獻(xiàn)為正，這可能與剝離貢獻(xiàn)方法的誤差有關(guān)；就參數(shù)n，其對(duì)ET變化的貢獻(xiàn)均為正，且大于50 mm，尤以在尹府最大（111.55 mm）。1990s，降水貢獻(xiàn)均為負(fù)，除尹府外，其他3個(gè)流域的貢獻(xiàn)率均大于65 mm，且臧格莊最大（-99.94 mm）；由于產(chǎn)芝、團(tuán)旺和尹府PET增加，ET增大（＜15 mm），而臧格莊PET的減小使得ET略微減小。2000s，降水對(duì)ET變化的貢獻(xiàn)依然為負(fù)（＞-25 mm）；PET的減小使得各流域ET減?。ǎ?12 mm）；隨參數(shù)n的增大，各流域ET不同程度的增大，其中產(chǎn)芝和尹府的貢獻(xiàn)大于105 mm，而臧格莊的貢獻(xiàn)最?。?6.46 mm）。由圖4（b）～（d），臧格莊各時(shí)段PET、參數(shù)n對(duì)ET變化的貢獻(xiàn)均小于降水，故臧格莊1980s、1990s和2000s的ET變化主要?dú)w因于降水的減少；與之相反，其他流域各時(shí)段均顯示參數(shù)n的貢獻(xiàn)最大，故產(chǎn)芝、團(tuán)旺和尹府各時(shí)段ET變化的主控因子為參數(shù)n。另外，需要指出的是，臧格莊ET變化主控因子不同于其他流域，其原因可能為降水［圖2（b）］、參數(shù)n［圖2（d）］的變化幅度以及ET對(duì)它們的敏感性（即單位降水或參數(shù)n變化引起的ET變化）存在一定區(qū)域差異。

圖4 不同時(shí)段ET變化及降水、PET和參數(shù)n對(duì)各時(shí)段ET變化的貢獻(xiàn)Fig.4 Changes in ET during different periods，accompanied with contributions of precipitation，PET and parameter n

2.4 徑流變化及其歸因

由圖2（a），除產(chǎn)芝基準(zhǔn)期徑流小于200 mm，其他流域均大于200 mm，且尹府最大，為242.65 mm。就各時(shí)期徑流變化，各流域均不同程度減?。蹐D5（a）］，1980s 徑流下降最大（＜-130 mm），而2000s 下降較?。ǎ?120 mm）。如圖5（b）～（d），1980s，降水減小使各流域徑流減小，除產(chǎn)芝外，其貢獻(xiàn)率均大于50 mm；由于PET增大，產(chǎn)芝、團(tuán)旺和尹府徑流減小（＞-20 mm）；參數(shù)n的增加使各流域徑流減小，以產(chǎn)芝和尹府減小最大（＜-90 mm）。1990s，降水對(duì)徑流變化的貢獻(xiàn)為負(fù)，其中團(tuán)旺和臧格莊的貢獻(xiàn)較大（＜-60 mm），且以臧格莊最大（-70.11 mm）；隨著PET的增加，產(chǎn)芝、團(tuán)旺和尹府徑流減?。ǎ?15 mm），而臧格莊PET使得徑流略微增大。2000s，降水對(duì)各流域徑流變化的貢獻(xiàn)為負(fù)（＞-30 mm）；PET對(duì)各流域徑流變化貢獻(xiàn)為正（＜12 mm）；隨著參數(shù)n的增大，各流域徑流不同程度減小，其中產(chǎn)芝和尹府貢獻(xiàn)較大（＜-105 mm），而團(tuán)旺和臧格莊貢獻(xiàn)較?。ǎ?50 mm）。對(duì)比各因子對(duì)徑流變化的貢獻(xiàn)［圖5（b）～（d）］，可以發(fā)現(xiàn)，1980s，產(chǎn)芝和尹府流域參數(shù)n對(duì)徑流下降的貢獻(xiàn)最大，而團(tuán)旺和臧格莊流域降水對(duì)徑流減小的貢獻(xiàn)最大，故產(chǎn)芝和尹府、團(tuán)旺和臧格莊徑流變化的主控因子分別為參數(shù)n、降水；1990s 和2000s，除臧格莊主控因子為降水（降水貢獻(xiàn)最大），而其他流域均顯示參數(shù)n對(duì)徑流變化的貢獻(xiàn)最大，故其他3 個(gè)流域徑流的變化可歸因于參數(shù)n的增大。值得注意的是，流域間徑流變化主控因子存在差異，這可能與降水［圖2（b）］、參數(shù)n［圖2（d）］的變化幅度以及徑流變化對(duì)它們的敏感性（即單位降水或參數(shù)n變化引起的徑流變化）的區(qū)域性差異有關(guān)。

圖5 不同時(shí)段徑流變化及降水、PET和參數(shù)n對(duì)各時(shí)段徑流變化的貢獻(xiàn)Fig.5 Changes in runoff during different periods，accompanied with contributions of precipitation，PET and parameter n

3 結(jié)論與討論

以膠東半島4 個(gè)典型流域?yàn)槔?，本研究詳?xì)分析了各流域降水、PET、ET和徑流變化特征，并基于Budyko方程和多控制因子聯(lián)立求解方法，估算了氣候和參數(shù)n對(duì)ET和徑流變化的貢獻(xiàn)，開展了歸因分析，得到以下結(jié)論。

（1）較基準(zhǔn)期1961-1980，各流域各時(shí)段降水均減小，以1980s 和1990s 減小最大（＜-100 mm）。PET變化表現(xiàn)出明顯的區(qū)域性差異，產(chǎn)芝、團(tuán)旺和尹府在1980s 和1990s 增大，而在2000s減小，臧格莊各時(shí)段PET均減小，且1990s和2000s減小幅度較大（＜-55 mm）。各流域不同時(shí)段參數(shù)n均增大，尤其產(chǎn)芝1980s和2000s增加最大（＞1.0）。

（2）定性和定量化指標(biāo)均顯示多控制因子聯(lián)立求解方法可以較為準(zhǔn)確地估算各影響因子對(duì)ET、徑流變化的貢獻(xiàn)。

（3）產(chǎn)芝、團(tuán)旺和尹府各時(shí)段ET均增加，其中產(chǎn)芝和尹府2000s 的ET增加最明顯（＞100 mm）；相反，臧格莊不同時(shí)段的ET較基準(zhǔn)期均減小，且1990s 減小最大（-35.60 mm）。不同時(shí)段ET的變化，除臧格莊主控因子為降水，其余流域ET增加均歸因于參數(shù)n的增大。

（4）各流域徑流均減小，1980s 下降最大（＜-130 mm），而2000s 下降較?。ǎ?120 mm）。產(chǎn)芝和尹府1980s 徑流的下降歸咎于參數(shù)n的增大，而團(tuán)旺和臧格莊徑流減小的主控因子為降水；1990s 和2000s，除臧格莊主控因子為降水，其余流域徑流變化的主控因子均為參數(shù)n。

通過分析膠東半島4 個(gè)流域不同時(shí)段ET和徑流（較1961-1980年）的變化，揭示了在全球變化背景下該地區(qū)水資源量明顯下降。采用Budyko方程和多控制因子聯(lián)立求解方法對(duì)ET和徑流變化進(jìn)行了歸因分析，甄別了控制各流域地表水文分量變化的主要因子，有助于各流域制定合適的應(yīng)對(duì)水資源短缺的措施。對(duì)于參數(shù)n（即人類活動(dòng)）主導(dǎo)的流域，政策制定者可通過改變?nèi)祟惢顒?dòng)等，減小ET耗散，增大徑流量，提高水資源儲(chǔ)備；而對(duì)于降水變化主導(dǎo)的流域，政策制定者可考慮修建水利工程設(shè)施，儲(chǔ)備降水，提高水資源利用效率，以減少水資源短缺帶來的負(fù)面影響。

本研究存在一定的不確定性。如資料的處理：未考慮地表反照率（設(shè)置為常數(shù)0.23）的時(shí)空變化［27，28］估算凈輻射，可能會(huì)對(duì)PET估算造成偏差。氣象觀測(cè)（尤其是降水）經(jīng)常受儀器和觀測(cè)場(chǎng)條件的影響［29］，進(jìn)而給本研究帶來不確定性。FAO56 Penman-Monteith 公式：FAO56 Penman-Monteith 公式并未考慮CO2濃度升高對(duì)植被生理特征（如氣孔導(dǎo)度和植被結(jié)構(gòu)特征）的影響，進(jìn)而造成PET估算存在偏差［30］。剝離各因子貢獻(xiàn)的方法：該方法的核心假設(shè)為影響因子的貢獻(xiàn)是線性疊加的，但事實(shí)并非如此，進(jìn)而造成一定的不確定性。 □