基于概率Budyko方程的流域未來可用水資源比例預估

陳澤峰，王衛光，李長妮，丁一民，傅健宇

(1.河海大學 水文水資源與水利工程科學國家重點實驗室，南京 210098；2.河海大學水文水資源學院，南京 210098；3.陜西省河流工程技術研究中心，西安 710018)

近年來，伴隨著全球變暖和人類活動的影響，降水、陸面蒸散發、流域下墊面、大氣化學組分等均發生顯著變化[1,2]，進而影響水文循環過程[3]，使得流域徑流過程發生改變[4]、極端水文事件頻繁發生，也加劇了我國水資源短缺、供需矛盾突出、水環境惡化等的水資源狀況。因此，精確的評估陸面可利用水資源量變得至關重要。

Budyko水熱耦合平衡方程[5,6]由于物理機制明確、輸入資料少，因此在大時間尺度上可利用水資源量的評估和預測工作中得到了廣泛的應用[7,8]。Budyko理論最初是表示多年平均蒸發率(實際蒸散發AET/降水P)和干旱指數(潛在蒸散發PET/降水P)關系的無參數方程。隨后，為了刻畫不同下墊面特征的流域Budyko曲線形狀，含有參數ω的Budyko方程被逐漸發展起來。含參數的Budyko方程對Budyko曲線的擬合精度較高，也因此應用較為廣泛。雖然，眾多學者研究過參數ω與單個或多個流域要素間的關系[7,9]，但由于參數ω與觀測變量間的關系復雜，對于不同氣候狀況、不同流域的ω和各因素間關系的系統變化仍然沒有一個確定的結論。對此，Greve等[10]用美國流域的觀測數據，建立了從理論上評估不同流域特征可用水資源量的概率Budyko方程。隨后，基于上述理論，Wang和Fu[11]結合流域下墊面特征，改進了概率Budyko方程，提出了雙變量概率性Budyko框架，提高了可利用水資源量的準確性。概率Budyko的發展，為評估相同氣候條件下，不同下墊面特征的流域可用水資源及其不確定性提供了新的視角。將基于下墊面因素的參數ω的概率分布代入Budyko方程，可以定量評估未來水資源預估的不確定性，從而為水資源管理者和決策者提供更精確的依據。

本文結合IPCC 第五次報告 CMIP5 的氣候模式數據(CanESM2)，運用統計降尺度方法，預估我國不同干濕特征的4個典型流域(黑河流域、灤河山區、漢江丹江口以上流域和贛江棟背以上流域)未來水文氣象變量，分析其變化特征；同時，基于構建的雙變量概率性Budyko方程，預估典型流域不同排放情景下未來可用水資源比率的分布，并評估其不確定性。

1 流域概況及資料

本文選取4個典型流域(黑河流域、灤河山區、漢江丹江口以上流域和贛江棟背以上流域)為研究區域。4個流域具有不同的干濕特征：黑河流域代表干旱流域、灤河山區代表半干旱流域、漢江丹江口以上流域代表半濕潤流域、贛江棟背以上流域代表濕潤流域。研究區位置及氣象站點分布情況如圖1所示。

圖1 4個典型流域及其氣象站點的分布Fig.1 Spatial distribution of four typical catchments and locations of meteorological stations

研究區歷史觀測數據包括從國家氣象局(http:∥cdc.cma.gov.cn/)搜集的4個流域內共21個氣象站點1961-2000年間逐日降水、最低氣溫、最高氣溫、日照時數、平均相對濕度和平均風速資料。潛在蒸散發通過具有明確物理機制的Penman-Monteith方法[12]計算。各流域上降水、潛在蒸散發通過泰森多邊形法[13]計算流域內部氣象站點的面積權重來得到。

未來數據本文選擇IPCC 第五次報告 CMIP5 的CanESM2(the second generation of Canadian Earth System Model)氣候模式網格數據，CanESM2氣候模式由加拿大環境氣候變化部門的氣候模型和分析中心提供，其空間分辨率為2.812 5°×2.812 5°，時間序列為2006-2100年。CanESM2氣候模式包括3個排放情景，分別是RCP 26、RCP 45和RCP 85，依次描述的是低排放情景，中等排放情景和高排放情景。

再分析資料(NCEP，National Centers for Environmental prediction)同樣來自加拿大環境氣候變化部門的氣候模型和分析中心。該NCEP數據集的時間長度是1961-2005年，空間尺度與CanESM2氣候模式相同(2.812 5°×2.812 5°)。NCEP數據包括26個大氣因子，主要用來率定和驗證統計降尺度模型，以預測未來不同排放情景下的氣候變量。

2 研究方法

2.1 概率Budyko方程

1981年，傅抱璞[14]在Budyko假設的基礎上，根據流域水文氣候的物理機制，經過量綱分析和數學推導，得到了Budyko假設的首個解析表達式：

(1)

式中：ω為表征流域自然地理特性的無量綱參數，且ω≥1。

由于傅抱璞公式的推導過程為Budyko假設提供了數理基礎，具有物理和數學意義，所以與其他經驗的Budyko方程相比，傅抱璞公式的認可程度較高并且應用較為廣泛。

相比于傳統的Budyko方程給定一個干旱指數值只能得到一個蒸發率，概率Budyko方程通過建立參數ω的概率分布得到的是蒸發率的一個分布，繼而可以定量評估相同氣候條件下，不同下墊面特征流域可利用水資源量的不確定性。

本文以傅抱璞的Budyko水熱耦合平衡方程為依據，通過Copula函數，構造Budyko方程參數ω與流域下墊面因素(植被覆蓋指數、土壤可用含水量、土壤容重、流域高程、坡度、復合地形指數)的邊緣分布函數。經過參數計算、擬合優度評價等選取最優的Copula函數，分別建立ω和各因素的聯合分布函數，以植被覆蓋指數(NDVI)為例，聯合概率分布見圖2。隨后，計算ω在流域下墊面條件下的條件分布函數，反推出在下墊面各因素限制后的參數ω，進而構建基于流域下墊面特征的概率Budyko方程。具體方法見Wang和Fu[11]。

圖2 參數ω與植被覆蓋指數的聯合概率分布圖Fig.2 Joint probability distribution of ω and NDVI

2.2 統計降尺度模型SDSM

統計降尺度(Statistical Downscaling Model)的核心是建立大尺度氣候因子與站點觀測變量之間的關系。由于運行速度快且成本低，統計降尺度相比動力降尺度成為眾多學者首選的方法[15]，也因此被廣泛應用于評估氣候變化背景下氣象、水文、環境等的響應[16,17]。

統計降尺度的主要過程可分為以下4步：預測因子的篩選；構建模型；模型的率定和驗證；未來氣象變量的預測。

3 結果與分析

3.1 統計降尺度模型的率定與驗證

預測因子的篩選是統計降尺度的核心環節，主要依據觀測變量和多個再分析大氣因子間的相關性和偏相關性，從而篩選出最相關的預測因子。本文計算了4個典型流域日尺度的潛在蒸散發與降水和26個氣候因子間的相關系數(corr)和偏相關系數(p-corr)，并從中選取最相關的預報因子。

因子篩選之后，進行模型計算和率定驗證。本文以1961-1990年作為模型率定期，1991-2000年為模型驗證期，用相關系數R、模擬值標準差與觀測值標準差的比值RS、確定性系數Ens和偏差Bias4個指標來評價模型，以灤河山區的6個氣象站點為例，評價結果見表1。

可以看出，模型對潛在蒸散發的模擬效果較好，R和RS在率定期和驗證期均超過了0.92，多倫站和青龍站的RS甚至等于1；確定性系數Ens也都大于0.80；偏差Bias大都接近于0。降水的模擬結果雖然不如潛在蒸散發好，但R都在0.43以上，RS都大于0.57，也能用來預估未來的氣象變量。

流域尺度上，比較了驗證期(1991-2000年)實測月潛在蒸散發和模擬月潛在蒸散發。黑河流域、灤河山區、漢江丹江口以上流域和贛江棟背以上流域的R2分別達到0.993、0.992、0.977和0.975，說明整體模擬效果好。同時比較了驗證期(1991-2000年)實測月降水和模擬月降水。黑河流域、灤河山區、漢江丹江口以上流域和贛江棟背以上流域的R2分別達到0.800、0.912、0.811和0.812，4個流域的R2都大于0.80，說明模擬效果較好。總之，SDSM在構建降水和潛在蒸散發時具有良好的模擬效果，可用于對未來數據進行降尺度。

3.2 流域未來水文氣象變量的預估

根據建立好的SDSM，輸入不同排放情景的預報因子，便可得到未來不同情景下的日潛在蒸散發量和日降水量，進而累加計算得到年潛在蒸散發量和降水量，其結果如圖3所示。4個流域的潛在蒸散發量在未來RCP 26、RCP 45和RCP 85情景下均呈上升趨勢，且RCP 85情景下的增加趨勢更為明顯，尤其在2050年之后。降水的整體增加趨勢相對較小，尤其在濕潤的贛江棟背以上流域。

4個流域未來時段2021-2050年3個排放情景的多年平均降水量(P)和潛在蒸散發量(PET)，繼而計算得到相應的干旱指數(AI)。結果表明，在RCP26情景下，黑河流域、灤河山區、漢江丹江口以上流域以及贛江棟背以上流域的AI分別為8.944、1.651、1.144和0.884；在RCP45情景下，黑河流域、灤河山區、漢江丹江口以上流域以及贛江棟背以上流域的AI分別為8.797、1.727、1.145和0.867；在RCP85情景下，黑河流域、灤河山區、漢江丹江口以上流域以及贛江棟背以上流域的AI分別為8.745、1.710、1.131和0.879。

表1 統計降尺度模型在構建灤河山區各氣象站點日潛在蒸散發和日降水時的模擬結果評價Tab.1 Evaluation of daily potential evapotranspiration and precipitation at each at each meteorological station of Luanhe Mountain that simulated by statistical downscaling model

圖3 流域歷史實測潛在蒸散發、降水和CanESM2模式不同氣候情景下的潛在蒸散發、降水Fig.3 Measured potential evapotranspiration and the potential evapotranspiration under the different emission scenarios based on CanESM2 for four catchments

3.3 流域2021-2050年可用水資源比率的預估

在構建概率Budyko方程的過程中，通過對比下墊面各因素限制后的參數ω′和ω的概率密度函數，得到了各下墊面因素對參數ω的影響程度。結果表明，黑河流域和贛江棟背以上流域的最大影響因素為植被覆蓋指數，漢江丹江口以上流域的最大影響因素為土壤可用含水量，而灤河山區的最大影響因素為高程。因此，根據建立的概率Budyko方程，假定流域2021-2050年參數ω分布與歷史時期(1982-2008年)ω的分布保持一致，即黑河流域和贛江棟背以上流域的參數是植被覆蓋指數限定后的ω′，灤河山區和漢江丹江口以上流域的參數分別是由高程和土壤可用含水量限定后的ω′。將4個典型流域對應的ω′及得到的未來不同排放情景下干旱指數(AI)分別代入Budyko方程，便可得到未來可用水資源比率的分布，從而評估其不確定性。

圖4展示了4個流域不同排放情景下2021-2050年多年平均可用水資源比率的分布。可發現，黑河流域可用水資源比率值分布范圍較窄，多在0～0.06之間，不同排放情景的可用水資源比率值較為接近[圖4(a)]。灤河山區可用水資源比率范圍相對較寬，多集中在0～0.5之間，且有雙峰現象，說明有兩段值出現的頻率較大，0.1～0.25、0.25～0.4為兩個雙峰的大概位置；RCP 26 情景與另外兩個排放情景下的可用水資源比率略有差異[圖4(b)]。漢江丹江口以上流域可用水資源比率大都集中在0.2～0.6間，雙峰現象略不明顯，不同排放情景的可用水資源比率值較為接近[圖4(c)]。贛江棟背以上流域可用水資源比率多分布在0.3～0.8的范圍內，雙峰現象明顯，且在峰值以前的曲線斜率相對峰值以后的曲線斜率較小，不同排放情景的可用水資源比率值差異較小[圖4(d)]。

圖4 流域2021-2050年不同排放情景下可用水資源比率的分布Fig.4 Water availability ratio under the different emission scenarios during 2021-2050 for four catchments

從可用水資源比率的箱線圖(圖5)可看出，各流域可用水資源比率在不同排放情景下的分布極為相似，表明不同的排放情景對可用水資源比率的影響不大；與常理一致，可用水資源比率值隨著流域濕潤程度的增加而增加，黑河流域、灤河山區、漢江丹江口以上流域和贛江棟背以上流域可用水資源比率的中位數依次為0.10、0.20、0.38和0.57；箱子的高度(即四分位距IQR)在黑河流域最矮(0.022)，漢江丹江口以上流域次之(0.144)，接著是贛江棟背以上流域(0.151)，灤河山區最高(0.194)，說明可用水資源比率的不確定性在灤河山區最大，黑河流域最小，其余兩個流域介于之間。此外，黑河流域和贛江棟背以上流域的下墊面最大影響因素同為植被覆蓋指數，兩個流域的Budyko參數均為由植被覆蓋指數限定之后的參數ω，但其可用水資源比率的大小及不確定性卻相差甚大，說明流域氣候特征對流域水資源狀況影響強烈，且加入植被覆蓋指數后可用水資源比率的不確定性在干旱區域(黑河流域)比濕潤區域(贛江棟背以上流域)小。

箱子最底端和最頂端的短線分別表示可用水資源比率的最小值和最大值，箱子從下往上的三條線依次表示25%、50%、75%處的值，實心圓代表可用水資源比率的平均值圖5 流域2021-2050年不同排放情景下可用水資源比率的箱線圖Fig. 5 Boxplot of water availability ratio under the different emission scenarios during 2021-2050 for four catchments

總之，不同排放情景對流域可用水資源比率的影響較小，流域可用水資源比率值隨著流域濕潤程度的增加而增加，可用水資源比率的不確定性從高到低依次是灤河山區、贛江棟背以上流域、漢江丹江口以上流域和黑河流域。此外，濕潤區贛江棟背以上流域的可用水資源比率的不確定性略高于半濕潤區漢江丹江口以上流域的這一結果也反映出研究不同下墊面因素對參數ω影響程度的必要性。

4 結 語

本文運用統計降尺度模型(SDSM)，結合CanESM2氣候模式、實測降水、由Penman-Monteith公式計算的潛在蒸散發，預估未來3種不同輻射強迫和大氣排放標準(RCP 26、RCP 45、RCP 85)下的降水量和潛在蒸散發量；基于構建的雙變量概率Budyko方程，對我國不同干旱特征的典型流域(黑河流域、灤河山區、漢江丹江口以上流域和贛江棟背以上流域)未來可用水資源比率進行預估，并分析其不確定性。主要結果如下：

(1)4個典型流域的潛在蒸散發量在未來RCP 26、RCP 45和RCP 85情景下均呈上升趨勢，且在2050年之后，RCP 85情景下的潛在蒸散發量最大，RCP 45情景下的潛在蒸散發量次之，RCP 26情景下的潛在蒸散發量最小；降水的整體增加趨勢則相對較小，尤其在濕潤的贛江棟背以上流域。

(2)黑河流域、灤河山區、漢江丹江口以上流域和贛江棟背以上流域可用水資源比率的中位數依次為0.10、0.20、0.38和0.57，即4個流域2021-2050年多年平均可用水資源比率值隨著濕潤程度的增加而增加；不同排放情景下的可用水資源比率分布的差異較小；灤河山區和贛江棟背以上流域的可用水資源比率分布有雙峰現象。

(3)不同排放情景對流域可用水資源比率的不確定性影響較小；4個流域可用水資源比率的不確定性從高到低依次是灤河山區(IQR為0.194)、贛江棟背以上流域(IQR為0.151)、漢江丹江口以上流域(IQR為0.144)和黑河流域(IQR為0.022)。