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

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

（1.南京信息工程大學氣象災害教育部重點實驗室，南京210044；2.常州市金壇區氣象局，江蘇常州213200）

0 引 言

隨著社會經濟的快速發展和人口的不斷增長，人類對水資源的需求增大，水資源短缺問題日益加劇，已成為影響社會經濟可持續發展的重要約束條件［1］。作為水文循環的重要環節，地表水文過程始終與水資源的合理開發利用密切相關；因此，了解和認識地表水文過程對全球變化的響應已成為水文、氣候等科學研究的熱點科學問題［2-4］。

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

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

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

1 資料與方法

1.1 資料及處理

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

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

1.2 Budyko方程

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

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

1.3 歸因方法

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

表1 試驗設計Tab.1 Configurations of experiments

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

2 結果分析

2.1 降水、PET及參數n變化特征

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

圖2 不同時期各流域主要水文氣候分量和參數n及降水、PET、參數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 歸因方法評估

在進行歸因分析之前有必要對剝離因子貢獻的方法進行評估，本研究采用各因子貢獻之和與觀測的ET、徑流變化散點圖，以及相關系數（R）、趨勢系數和均方根誤差（RMSE）等指標，結果如圖3。很明顯，該方法剝離的各因子貢獻之和與觀測的ET、徑流變化幾乎處在1∶1 線上。由定量化指標，各因子貢獻之和與觀測的ET、徑流變化的R幾乎為1，RMSE為5.12 mm，擬合趨勢線的趨勢系數也近乎為1；均說明該剝離因子貢獻方法可以較為準確地估算各影響因子對ET、徑流變化的貢獻。

圖3 歸因方法評估Fig.3 Validation for the selected attribution method

2.3 蒸散發變化及其歸因

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

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

2.4 徑流變化及其歸因

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

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

3 結論與討論

以膠東半島4 個典型流域為例，本研究詳細分析了各流域降水、PET、ET和徑流變化特征，并基于Budyko方程和多控制因子聯立求解方法，估算了氣候和參數n對ET和徑流變化的貢獻，開展了歸因分析，得到以下結論。

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

（2）定性和定量化指標均顯示多控制因子聯立求解方法可以較為準確地估算各影響因子對ET、徑流變化的貢獻。

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

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

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

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