基于WEP的橫斷山垂直帶水文過程模擬研究

蘇輝東，賈仰文，牛存穩(wěn)，倪廣恒，劉佳嘉，劉 歡，杜軍凱

（1.中國水利水電科學(xué)研究院 流域水循環(huán)模擬與調(diào)控國家重點(diǎn)實(shí)驗(yàn)室，北京 100038；2.清華大學(xué) 水利水電工程系，北京 100084）

1 研究背景

我國是一個山地多而人均水資源少的國家。我國70%的陸域面積為山地，40%多的人口生活在山區(qū)。而我國又是一個人均水資源短缺的國家，人均水資源量僅為全球1/4的水平，且絕大部分水資源分布在山區(qū)[1]。橫斷山是長江、怒江、瀾滄江等江河重要的水源地，僅長江上游流域（宜昌以上）多年平均徑流量約為4515億m3，約占長江平均總水量的46%。以橫斷山為代表的亞洲高山區(qū)為“亞洲的水塔”，對數(shù)億人的持續(xù)供水做出了重大貢獻(xiàn)[2-3]。開展山地水文過程模擬及水資源演變規(guī)律研究，可為我國水資源安全保障提供支撐，具有重要科學(xué)意義。

山地水文過程及演變規(guī)律研究需要用到水文模型。水文模型可以分為集總式模型（如：TANK模型[4]、新安江模型[5]）和分布式水文模型（如：SHE[6]、TOPMODEL[7]、SWAT[8]、VIC[9]、WEP（Water and Energy transfer Processes）等）。山地水文過程模擬面臨很多復(fù)雜問題，主要包括時空尺度問題和不均介質(zhì)問題。時空多尺度一直以來是山地水文過程研究的重點(diǎn)問題，也是難點(diǎn)問題。在不同水文尺度上，其水循環(huán)和水文要素表現(xiàn)不一。針對這個問題，國內(nèi)外學(xué)者進(jìn)行了許多研究。Samuel和Sivapalan[10]在澳洲的三個不同流域分析了場次降雨、季節(jié)降雨、年內(nèi)降雨、年際降雨等不同尺度對洪水頻率曲線的影響規(guī)律。Liu等[11]使用EEMD法從3年、5年、12年和46年等周期分析了雅魯藏布江流域徑流波動和它對氣候變化的響應(yīng)，結(jié)果表明年際振蕩對雅魯藏布江的徑流變化起著重要作用。徐飛[12]從年季月多時間尺度模擬分析了五臺山清水河流域、崇陵小流域、拒馬河流域多空間尺度的水循環(huán)過程，結(jié)果揭示了不同時空尺度下流域水文過程對氣候變化和土地利用的響應(yīng)不同。不均介質(zhì)問題研究包括大孔隙優(yōu)先流、喀斯特溶洞、巖石裂隙、基巖凹凸面、土石介質(zhì)等方面研究。另外動植物影響[13]、干濕/凍融交替作用[14-15]、人為擾動（如耕地等）等因素會增加介質(zhì)的不均質(zhì)性。土石二元介質(zhì)對入滲產(chǎn)流影響具有不確定性，碎石一方面會減少入滲，而碎石與土壤接觸部分又會增加孔隙度，從而增加入滲[16]。例如：片麻巖介質(zhì)和石灰?guī)r介質(zhì)透水性和持水性各異，對入滲產(chǎn)流會造成不同的影響[17]。受巖溶裂隙[18]影響，喀斯特山地[19]土層薄、蓄水能力弱，補(bǔ)徑排迅速，水位流量季節(jié)變化劇烈。此外，土壤斥水性[20]也影響土壤入滲產(chǎn)流，斥水性越大，入滲率就會減小。

國內(nèi)外已在這方面開展了一些相關(guān)實(shí)驗(yàn)和分析研究，集中在對野外小流域典型觀測[21]，通過分析土壤含水率、蒸散發(fā)量等因素，研究山地流域不同山坡結(jié)構(gòu)、地質(zhì)和地形對徑流的影響。就山坡水文過程與產(chǎn)匯流機(jī)制研究總體而言，需要綜合考慮氣候、地形和地質(zhì)條件約束。目前，盡管國內(nèi)外學(xué)者提出了眾多分布式水文模型，但針對山區(qū)尤其是山坡產(chǎn)匯流、無資料地區(qū)水文模擬等問題的研究仍較薄弱。分布式流域水文模型對山地水文過程垂直地帶性與演變特征的研究不夠，尤其是受山地復(fù)雜地形、氣候、下墊面條件的限制，給參數(shù)方案和模擬結(jié)果帶來了很大不確定性，相關(guān)研究尚不足。

橫斷山海拔落差大，導(dǎo)致水-熱-植被分布具有顯著的垂直地帶性，有“一山四季”和“一山四帶”的現(xiàn)象。氣象、植被、土壤、水分具有顯著垂直地帶性特點(diǎn)，山地垂直地帶性導(dǎo)致參數(shù)不確定性很大，氣候條件惡劣，地形復(fù)雜，位置偏僻，導(dǎo)致數(shù)據(jù)缺乏（1萬km2不足一個氣象觀測站點(diǎn)）。因此，研究橫斷山垂直地帶性對水文過程的影響具有挑戰(zhàn)性與科學(xué)價(jià)值，同時又有很重要的實(shí)踐意義。本文以橫斷山垂直帶水文過程模擬研究為目的，采用分布式流域水文模型WEP，提出垂直帶參數(shù)化方案以更好地刻畫垂直帶水文過程，提高山地垂直帶的水文模擬精度。

2 研究區(qū)概況與研究方法

2.1 研究區(qū)概況橫斷山（24°35′—33°34′N，96°56′—104°30′E）位于我國西南部，青藏高原東南部，面積約50萬km2，其中90%以上為山地。橫斷山位于我國第一、第二階梯的交界地帶，新構(gòu)造運(yùn)動活躍。橫斷山海拔高程差值大，導(dǎo)致水分、溫度、植被垂直地帶性顯著。橫斷山區(qū)自西向東有怒江、瀾滄江、金沙江、雅礱江、大渡河、岷江等河流，其分布如圖1所示。橫斷山地區(qū)山高谷深，水流豐沛，主要河流均發(fā)源于青藏高原地區(qū)，大江大河主要為自北向南流。

圖1 橫斷山流域邊界及水系分布

2.2 研究方法WEP分布式流域水文模型[22]，能夠用于模擬流域內(nèi)水、熱（通量）的運(yùn)移過程。采用子流域套等高帶作為基本計(jì)算單元[23]、并開發(fā)了流域二元水循環(huán)模型[24]、嵌套了跨區(qū)域水資源配置耦合模型[25]，使得模型日趨完善。WEP模型先后在日本、韓國、中國等各種尺度、復(fù)雜地質(zhì)氣候條件下的流域應(yīng)用，取得了良好模擬效果。WEP模型結(jié)構(gòu)如圖2所示，水循環(huán)、能量過程的模擬方法詳見參考文獻(xiàn)[20]。

2.3 數(shù)據(jù)來源橫斷山分布式流域水文模型的基礎(chǔ)數(shù)據(jù)主要包括以下五大類：（1）水文氣象數(shù)據(jù)主要來自觀測站點(diǎn)；（2）地理高程及地形數(shù)據(jù)；（3）河網(wǎng)水系數(shù)據(jù)；（4）土壤信息及水文地質(zhì)等數(shù)據(jù)；（5）土地利用及植被覆蓋數(shù)據(jù)。其數(shù)據(jù)來源匯總?cè)绫?所示。

圖2 WEP模型的水平結(jié)構(gòu)與垂直結(jié)構(gòu)

表1 數(shù)據(jù)資料描述及其來源匯總

3 WEP模型構(gòu)建與垂直帶參數(shù)化方案

3.1 橫斷山WEP模型的構(gòu)建

3.1.1 流域邊界確定 橫斷山區(qū)只是一個地理單元，并不是一個完整的流域。要構(gòu)建分布式流域水文模型，需要構(gòu)建一個或者若干個閉合的流域，使得這些流域能夠全部或者大體上囊括橫斷山區(qū)。根據(jù)橫斷山的三級水資源分區(qū)，需要適當(dāng)擴(kuò)大模擬域的邊界。橫斷山區(qū)東北角的小部分區(qū)域?yàn)榘缀铀担鋵儆邳S河流域，南部有小部分區(qū)域?qū)儆谠谒樱@兩部分沒有進(jìn)入本次的研究區(qū)范圍。本次構(gòu)建的橫斷山流域面積約85萬km2，主要包含“橫斷六江”和“長江源區(qū)”。橫斷山流域分布如下圖1所示。

3.1.2 河網(wǎng)提取 WEP河網(wǎng)包括實(shí)際河網(wǎng)和虛擬河網(wǎng)。實(shí)際河網(wǎng)由測量所得，包括五級河流（如圖3所示）。虛擬河網(wǎng)提取主要包括以下步驟：對DEM進(jìn)行裁剪、重分類、填洼等處理；根據(jù)實(shí)測的河網(wǎng)水系，對DEM進(jìn)行降高程處理；計(jì)算“流向”、“匯”、確定積水閾值，確定流域水庫、湖泊、流域出口等位置，制取分割點(diǎn)；提取計(jì)算得到虛擬河網(wǎng)。橫斷山流域的河網(wǎng)如圖3所示。

圖3 基于WEP提取的實(shí)際——虛擬水系

3.1.3 計(jì)算單元劃分 為既滿足垂直帶的模擬分析要求，也不造成過重的計(jì)算負(fù)擔(dān)，將橫斷山流域85萬km2的面積劃分成了具有拓?fù)潢P(guān)系的3089個子流域，再按照WEP的等高帶劃分規(guī)則，細(xì)分成14 995個基本計(jì)算單元。

3.1.4 輸入數(shù)據(jù)制備 WEP輸入的數(shù)據(jù)主要包括DEM、徑流數(shù)據(jù)、氣象數(shù)據(jù)、土地利用和土壤數(shù)據(jù)，其中DEM和氣象站點(diǎn)分布如圖4所示。

圖4 橫斷山DEM和氣象站點(diǎn)分布

3.2 垂直帶參數(shù)化方案的改進(jìn)垂直帶參數(shù)化方案主要是針對山地氣象、植被、土壤參數(shù)的垂直帶的空間分異而進(jìn)行處理的，旨在加強(qiáng)對山地垂直帶的刻畫和模擬，達(dá)到提高模擬精度和提升水文過程刻畫的準(zhǔn)確度。本文提出的垂直帶參數(shù)化的步驟包括：劃分獲得山地垂直帶譜；根據(jù)WEP及三級水資源區(qū)劃對參數(shù)進(jìn)行分區(qū)，通過ArcGIS將參數(shù)分區(qū)和垂直帶嵌套，獲得參數(shù)編號；通過衛(wèi)星遙感、站點(diǎn)觀測、野外試驗(yàn)、模型經(jīng)驗(yàn)等方法確定垂直帶參數(shù)值。

3.2.1 橫斷山垂直帶劃分 根據(jù)橫斷山的子流域形心點(diǎn)的高程，對橫斷山區(qū)進(jìn)行垂直帶劃分。橫斷山（291～7556 m）一共劃分成6個垂直帶，分別為：0～1000 m帶，1000～2000 m帶，2000～3000 m帶，3000～4000 m帶，4000～5000 m帶及＞5000 m帶。

3.2.2 垂直帶參數(shù)分區(qū) 根據(jù)橫斷山所屬的二、三級水資源分區(qū)和橫斷山6個垂直帶對橫斷山流域進(jìn)行參數(shù)化分區(qū)，構(gòu)建WEP模型，將橫斷山區(qū)分為9個基本參數(shù)區(qū)，其分區(qū)結(jié)果如圖5所示。

圖5 橫斷山流域海拔高程（DEM）及子流域劃分

3.2.3 垂直帶參數(shù)值的確定

（1）氣溫的垂直帶修正。WEP采用中國地面氣候資料數(shù)據(jù)集中的降水、氣溫、日照、風(fēng)速、相對濕度等5類數(shù)據(jù)進(jìn)行空間插值。受高程差大的影響，簡單的不考慮海拔變化的插值，不再適合橫斷流域氣象數(shù)據(jù)空間插值展布，因此需要進(jìn)行高程的修正。降水的修正比較復(fù)雜，本文主要針對溫度進(jìn)行修正。溫度的垂直帶修正步驟如下：

將氣象站點(diǎn)的溫度通過氣溫～高程的關(guān)系（海拔每上升100 m下降0.6℃）將氣溫進(jìn)行換算到參考的高程點(diǎn)溫度。本文將海平面作為參照點(diǎn)時，溫度修正為：

根據(jù)RDS反距離平均法計(jì)算插值點(diǎn)的參考溫度。結(jié)合RDS法計(jì)算公式如式（2）和式（3）：

將得到的氣溫?fù)Q算到子流域形心點(diǎn)的氣溫，將該氣溫作為輸入溫度，展布到整個子流域：

（2）植被參數(shù)的垂直帶方案。為了得到植被垂直帶的參數(shù)方案值，本文先用遙感數(shù)據(jù)進(jìn)行分析，然后結(jié)合實(shí)測數(shù)據(jù)進(jìn)行校正、規(guī)律分析總結(jié)，得到概化后的垂直帶植被參數(shù)值。

本文基于遙感（Modis）解譯的橫斷山流域多年平均月葉面積指數(shù)LAI和NDVI進(jìn)行分析，以100 m為等高帶，分別統(tǒng)計(jì)該等高帶內(nèi)的平均LAI和葉面積指數(shù)veg。其結(jié)果分別如圖6和圖7所示。其中Modis遙感數(shù)據(jù)沒有植被覆蓋度（veg）產(chǎn)品，需要根據(jù)NDVI數(shù)據(jù)估算區(qū)域植被覆蓋度[26]。植被覆蓋度veg的計(jì)算公式：

圖6 基于Modis的葉面積指數(shù)（LAI）隨高程變化結(jié)果

圖7 基于Modis的植被覆蓋度（veg）隨高程變化結(jié)果

式中：veg為像元內(nèi)的植被覆蓋度；N為像元NDVI數(shù)值；Nveg為像元全由植被所覆蓋下的 NDVI；Nsoil表示像元全由土壤覆蓋下的NDVI，Nsoil會隨時間變化，一般在-0.1～0.2之間[27]。根據(jù)NDVI制成的橫斷山流域多年平均月植被覆蓋度（veg）分布如圖7所示。這里需指出，概化后的垂直帶參數(shù)化表僅供參考，在率參時還需適當(dāng)調(diào)整參數(shù)。

最終，得到了橫斷山垂直帶植被葉面積指數(shù)（LAI）參數(shù)化的結(jié)果（表2）和垂直帶植被覆蓋度（veg）的參數(shù)化結(jié)果（表3）。其按照高程分成了6個參數(shù)化帶，且每個帶都概化了參考的1—12月際變化值。

表2 橫斷山垂直帶植被葉面積指數(shù)（LAI）參數(shù)

表3 橫斷山垂直帶植被覆蓋度（veg）參數(shù)

（3）土壤參數(shù)的垂直帶方案。土壤具有垂直地帶性，土壤參數(shù)對水文模型的模擬也影響很大。土壤層厚度和飽和含水率是高敏感性影響因子，飽和導(dǎo)水系數(shù)和濕潤鋒土壤吸力是中敏感性影響因子。橫斷山流域土壤類型分布見圖8。

圖8 橫斷山流域土壤類型分布

垂直帶土壤參數(shù)確定的基本步驟是：確定垂直帶上土壤的主要質(zhì)地組成，其方法如式（6）所示；建立土壤質(zhì)地和土壤水力特征參數(shù)（如土壤層厚度和飽和含水率等）之間的聯(lián)系或者一一對應(yīng)關(guān)系；再結(jié)合野外實(shí)測的土壤電導(dǎo)率、土壤溫度、濕度數(shù)據(jù)進(jìn)行校驗(yàn)分析；按照本文提出的垂直帶參數(shù)分區(qū)方法，確定每個垂直帶上的土壤參數(shù)的平均值：

微生物接種劑在農(nóng)業(yè)領(lǐng)域應(yīng)用技術(shù)已逐步成熟，廣泛應(yīng)用于多種牧草、糧食作物、蔬菜和中草藥等領(lǐng)域。研究表明，大豆植株接種根瘤菌和假單胞菌，大豆磷含量、鐵離子含量分別提高了88.9%、115.7%，且其全碳、全氮和豆血紅蛋白的含量以及根瘤數(shù)目均有增加[1]。Abbasi等[2]從小麥根際分離出具有分泌IAA功能的細(xì)菌并接種于小麥根際，發(fā)現(xiàn)該菌對小麥促生效果明顯，給辣椒接種芽孢桿菌(Bacillus)后辣椒產(chǎn)量和次生代謝產(chǎn)物(包括吲哚乙酸、鐵載體和幾丁質(zhì)酶)含量均可增加[3]。

式中：M，N分別為垂直帶上土壤柵格的行和列數(shù)；Sandij、Siltij和Clayij分別為砂粒、粉粒和黏粒的百分比含量；θsand、θsilt和θclay砂土、粉土和壤土的各類土壤參數(shù)。

再結(jié)合野外實(shí)測的土壤電導(dǎo)率、土壤溫度、濕度數(shù)據(jù)進(jìn)行分析。綜合上述分析，得到了垂直帶土壤的參數(shù)化結(jié)果，其結(jié)果如表4所示。

表4 垂直帶土壤水力特性參數(shù)初值

3.3 方案驗(yàn)證設(shè)計(jì)本文選擇的驗(yàn)證時間段為1960—2015年逐月天然徑流過程（其中參數(shù)率定期為1960—1990年，驗(yàn)證期為1991—2015年），選擇橫斷山流域的14個水文站點(diǎn)進(jìn)行驗(yàn)證對比。分別運(yùn)用三種參數(shù)方案運(yùn)行WEP模型，進(jìn)行模擬結(jié)果對比分析。14個站點(diǎn)分別為巴塘、直門達(dá)、奔子欄、石鼓、瀘寧、雅江、甘孜、屏山、夾江、彭山、三頭壩、洼里、瀘定、華彈，分布如圖9所示。

圖9 垂直帶參數(shù)化化方案選擇驗(yàn)證的站點(diǎn)分布

為驗(yàn)證垂直帶參數(shù)化方案的模擬效果，本文設(shè)置三種不同的參數(shù)化方案進(jìn)行對照。三種參數(shù)化方案分別為：

“垂直帶參數(shù)化方案”。如前文所述，將流域先劃分成若干個參數(shù)區(qū)，根據(jù)高程——空間插值、遙感——野外實(shí)測結(jié)合確定氣象、植被、土壤的概化垂直帶參數(shù)值，即為垂直帶參數(shù)化方案。在垂直帶參數(shù)化方案中，植被變量被當(dāng)成隨高程變化的概化參數(shù)，在模型中可以進(jìn)行調(diào)節(jié)。

“普通參數(shù)化方案”。氣象只采用RDS空間插值、土地利用、土壤、植被都是變量，不作參數(shù)處理，在模型中不能調(diào)節(jié)。其獲取方法為直接采用Modis數(shù)據(jù)等進(jìn)行統(tǒng)計(jì)獲得、不做參數(shù)處理。每個等高帶 （計(jì)算單元） 具有一個單獨(dú)的植被、土壤、水文參數(shù)。

“均一參數(shù)方案”。整個流域不進(jìn)行水資源分區(qū)，全流域用一個相同的參數(shù)。對整個流域Modis數(shù)據(jù)根據(jù)土地利用統(tǒng)計(jì)獲取森林、草地、作物三種類型的LAI和veg，所有等高帶具有相同的植被、土壤參數(shù)。

4 結(jié)果與分析

4.1 橫斷山垂直帶水文過程模擬驗(yàn)證采用水文站控制斷面月平均徑流量對整個模型進(jìn)行參數(shù)率定和驗(yàn)證，時間段為1960—2015年共56年的連續(xù)模擬（部分站徑流數(shù)據(jù)有缺失）。垂直帶參數(shù)化方案驗(yàn)證站點(diǎn)、驗(yàn)證時段設(shè)定和驗(yàn)證效果，見表5和圖10。

表5 基于垂直帶參數(shù)化方案的橫斷山流域徑流模擬結(jié)果

圖10 基于垂直帶參數(shù)方案的徑流過程模擬結(jié)果

在此基礎(chǔ)上，利用改進(jìn)后的WEP模型對橫斷山區(qū)域各水文站點(diǎn)徑流過程進(jìn)行模擬。14個水文站點(diǎn)的模擬月徑流與還原徑流在率定期的NS效率系數(shù)介于0.73～0.84，相對偏差RE范圍介于-12.1%～11.3%；驗(yàn)證期的NS效率系數(shù)介于0.71～0.85，相對偏差RE范圍介于-3.1%～8.7%。因此，總得來說，模型精度尚可，垂直帶參數(shù)化方案較好地模擬了橫斷山區(qū)垂直帶的水文過程。

4.2 垂直帶參數(shù)化方案對模擬改進(jìn)效果分析本研究通過設(shè)置三種不同的植被參數(shù)化方案進(jìn)行對比，分析植被參數(shù)隨高程變化對模型模擬結(jié)果的影響。各方案分別按月份設(shè)置不同的參數(shù)集，將這些參數(shù)集帶入模型中分別進(jìn)行模擬，不同方案下不同水文站模擬的NS效率系數(shù)和RE結(jié)果如圖11和圖12所示。

圖11 三種方案的NS效率變化對比

圖12 三種方案的相對偏差RE變化對比

從圖11中可以看出，“普通參數(shù)化”方案和“垂直帶參數(shù)化”方案NS效率系數(shù)基本相同，但前者相對誤差略高于后者；“均一參數(shù)化”方案的效果要低于其他兩種方案。由此可知，通過考慮植被參數(shù)的高程影響，能夠提高模型模擬效率。

從圖11及圖12中可以看出，垂直帶參數(shù)化方案比普通方案的NS效率系數(shù)和RE偏差都有所改進(jìn)，NS提高0.03，RE減少1%?？傮w來說改進(jìn)的并不明顯，但是也有部分改進(jìn)比較大的站點(diǎn)。而相比均一化方案，垂直帶參數(shù)化方案NS提高0.3、RE減少19.6%。垂直帶參數(shù)化方案的優(yōu)點(diǎn)：（1）相對于固定的計(jì)算，垂直帶參數(shù)化方案具有更多的靈活性，提高精度；（2）在沒有連續(xù)年份Modis數(shù)據(jù)的情況下，可通過高程分組獲取森林、草地、作物三種類型的植被參數(shù)；（3）參數(shù)簡化，易于獲??；（4）更好地反應(yīng)了氣象、植被、土壤、水分等信息的垂直地帶性，彌補(bǔ)了以往水文模型對參數(shù)的垂直變化描述不足。

高原山地區(qū)域在氣候、地形以及垂直梯度變化和地質(zhì)等方面存在著差異，植被垂直帶狀分布明顯，而原有WEP模型對植被參數(shù)的垂直變化描述不足。通過考慮氣象、植被、土壤參數(shù)隨高程變化的空間變異性，可提高山區(qū)水文模擬精度。

5 結(jié)論

本文針對橫斷山海拔落差大的特點(diǎn)，構(gòu)建了分布式流域水文模型WEP，根據(jù)橫斷山氣象—植被—土壤的垂直帶特征提出了參數(shù)化方案，評估了徑流模擬精度。主要結(jié)論如下：

（1）14個水文站點(diǎn)的模擬月徑流與還原徑流在率定期NS效率系數(shù)介于0.73～0.85，相對偏差RE范圍介于-12.1%～11.3%；驗(yàn)證期的NS效率系數(shù)介于0.73～0.88，相對偏差RE范圍介于-9.3%～11.3%?？偟脕碚f，垂直帶參數(shù)化方案較好地模擬了橫斷山區(qū)的徑流過程。

（2）“垂直帶參數(shù)化方案”比“普通的參數(shù)方案”的NS提高0.03、RE降低1%，比“均一化參數(shù)方案”的NS提高0.3、RE降低19.6%。通過考慮參數(shù)的垂直變化影響，能夠提高模型模擬效果。

（3）垂直帶參數(shù)方案具有以下特點(diǎn)：相對于固定普通參數(shù)化方案的計(jì)算，垂直帶參數(shù)化方案具有更多的靈活性，可提高精度；在沒有連續(xù)年份Modis數(shù)據(jù)的情況下，可通過高程分組獲取森林、草地、作物三種類型的植被參數(shù)；參數(shù)簡化，易于獲??；更好地反映了氣象、植被、土壤、水分等信息的垂直地帶性，彌補(bǔ)以往水文模型對參數(shù)的垂直變化描述不足。