基于SWAT模型和降水隨機(jī)模擬的徑流預(yù)測

張金萍，王宇昊

（1.鄭州大學(xué)水利科學(xué)與工程學(xué)院，鄭州450001；2.鄭州大學(xué)黃河生態(tài)保護(hù)與區(qū)域協(xié)調(diào)發(fā)展研究院，鄭州450001）

SWAT 模型基于物理過程建立，且能夠充分反映流域內(nèi)降水和下墊面要素空間變化對水文過程的影響，是研究流域水文循環(huán)的有效工具，是探究水文規(guī)律和機(jī)理的有效手段［1，2］。SWAT 模型在國內(nèi)外應(yīng)用廣泛，包括對水質(zhì)、徑流、泥沙等方面的研究，其中最為常見的就是徑流模擬［3］。例如王博威等依據(jù)四期土地利用數(shù)據(jù)構(gòu)建SWAT 模型，從而分析土地利用變化對流域徑流及其空間分布的影響［4］；廖亞一等將蒙特卡洛隨機(jī)采樣引入SWAT 模型，探討氣象數(shù)據(jù)輸入不確定性對徑流模擬結(jié)果的影響［5］；這些研究主要集中于通過構(gòu)建SWAT 模型探究土地利用變化、人類活動、氣候變化等因素對徑流的影響［6，7］，但降水是氣象數(shù)據(jù)中對徑流模擬結(jié)果影響最大的因素，降水的總量和年內(nèi)分布均會對徑流過程造成較大的影響［8］，而對不同量級和年內(nèi)分布降水情景相應(yīng)徑流過程的研究較少。

降水模擬包括降水概率的模擬和降水量的模擬，目前國內(nèi)外運(yùn)用比較多的降水概率模型主要有一階Markov 鏈、高階Markov鏈、間隔長度和半經(jīng)驗(yàn)分布法4種，而分布模型廣泛用于降水量模擬［9］。其中耦合Markov 鏈-Gamma 分布的隨機(jī)模擬模型既能利用Markov鏈將降水日和非降水日相聯(lián)系，也能在考慮降水年內(nèi)分布不均勻的前提下利用降水的概率分布較好的模擬降水日的降水量［10］。例如王斌等應(yīng)用模擬的降水序列，結(jié)合缺水率指標(biāo)劃分農(nóng)業(yè)旱情等級［11］；秦道清等將比例調(diào)整算法引入隨機(jī)模擬模型，生成日降水模擬序列［12］。但這些研究并沒有考慮到不同降水量級下降水年內(nèi)分布規(guī)律和模擬參數(shù)的差異性，因此根據(jù)年降水量級的不同分別構(gòu)建隨機(jī)模擬模型能有效提高模擬的代表性和準(zhǔn)確性。

本文的創(chuàng)新點(diǎn)在于將降水隨機(jī)模擬和SWAT 模型相結(jié)合，基于不同降水量級下降水年內(nèi)分布規(guī)律和模擬參數(shù)的差異性，通過構(gòu)建不同量級降水下的隨機(jī)模擬模型，生成未來年降水量為800 mm（豐水）、550 mm（平水）、300 mm（枯水）相對應(yīng)的日降水過程，并將其輸入SWAT模型預(yù)測徑流，從而進(jìn)一步探究降水和徑流的關(guān)系，具體流程見圖1。

1 數(shù)據(jù)與方法

1.1 研究區(qū)域概況及數(shù)據(jù)來源

沙河屬于海河流域大清河水系，發(fā)源于山西省靈丘縣，流經(jīng)河北省阜平境內(nèi)，最后注入王快水庫。王快水庫以上流域控制面積為3 770 km2，主河道長166 km，占沙河流域6 420 km2的59%。該流域多為山區(qū)，地表植被覆蓋較差，水土流失較為嚴(yán)重。王快水庫以上流域地處中溫帶半濕潤半干旱大陸季風(fēng)氣候區(qū)，多年平均氣溫12.5 ℃，最高氣溫41.7 ℃，多發(fā)生在7月份，最低氣溫-18.9 ℃，多發(fā)生在1月份［13］。

本文氣象數(shù)據(jù)選用CMADS 數(shù)據(jù)集［14，15］在研究區(qū)域內(nèi)及周邊的15 個站點(diǎn)和阜平站的氣象資料（最高、最低氣溫，平均風(fēng)速，日降水量，相對濕度，日照時數(shù)），徑流資料選用阜平水文站2008-2015年實(shí)測徑流資料，研究區(qū)域地理位置和站點(diǎn)分布見圖2。

圖2 研究區(qū)域和站點(diǎn)分布Fig.2 Study area and distribution of stations

1.2 SWAT模型

SWAT 是一種基于物理機(jī)制的半分布式水文模型，該模型以日為時間步長，能夠在不同的時間尺度（日、月和年）上進(jìn)行連續(xù)的長時段的模擬，通常設(shè)計用于利用流域和流域尺度的輸入數(shù)據(jù)對連續(xù)的長期徑流、泥沙和農(nóng)藥產(chǎn)量進(jìn)行建模和模擬［16］。SWAT模型將流域水文過程模擬分為水循環(huán)陸面過程和水循環(huán)水面過程，其中前者綜合了氣候、冠層截留、蒸騰蒸發(fā)、入滲、植被、地下水等產(chǎn)流和坡面匯流的影響因素，后者綜合了子河道傳輸過程、主河道匯流、洼地滯蓄過程等河網(wǎng)匯流的影響因素。因此，相比集總式模型，SWAT模型對于流域下墊面非均質(zhì)的描述更符合實(shí)際情況［17］，SWAT 模型的具體運(yùn)行流程見圖3。

圖3 SWAT模型運(yùn)行流程Fig.3 Operating process of SWAT model

1.3 耦合Markov鏈和Gamma分布的降水模擬

一階Markov 鏈模型將自然日分為降水日（P≥0.1 mm）和非降水日（P＜0.1 mm），模型的兩個基本參數(shù)是非降水日轉(zhuǎn)移到降水日的概率P（W/D）和降水日轉(zhuǎn)移到降水日的概率P（W/W）。研究表明當(dāng)月的兩個狀態(tài)轉(zhuǎn)移概率值均與同月降水日出現(xiàn)頻率q密切相關(guān)且存在線性關(guān)系，且轉(zhuǎn)移概率可用下式估算［18］：

若連續(xù)型隨機(jī)變量X服從Gamma 分布，則X的密度函數(shù)為：

式中：α為形狀參數(shù)，決定了Gamma 分布的形狀；β為尺度參數(shù)，決定了Gamma 分布的分散度。一階Markov 鏈模型的兩個轉(zhuǎn)移概率值和Gamma分布的兩個參數(shù)合稱為降水模擬參數(shù)。

隨機(jī)模擬步驟如下：

（1）以1月為例，假設(shè)1月1日為非降水日，在確定模擬參數(shù)的基礎(chǔ)上，將計算機(jī)產(chǎn)生的［0，1］之間的隨機(jī)數(shù)和P（W/D），P（W/W）進(jìn)行比較。

（2）當(dāng)前一日為非降水日時，若隨機(jī)數(shù)大于當(dāng)月的P（W/D），則當(dāng)日仍為非降水日，反之則為降水日；前一日為降水日時，若隨機(jī)數(shù)小于當(dāng)月的P（W/W），則當(dāng)日仍為降水日，反之則為非降水日；

（3）當(dāng)日為降水日時，將當(dāng)月降水的Gamma 分布結(jié)合隨機(jī)數(shù)產(chǎn)生當(dāng)日降水量。以此類推，即可產(chǎn)生若干年的逐日降水序列［19］。

2 SWAT模型的構(gòu)建

2.1 基礎(chǔ)數(shù)據(jù)庫的構(gòu)建

除氣象數(shù)據(jù)和徑流數(shù)據(jù)以外，SWAT 模型的構(gòu)建還需要流域的地形高程數(shù)據(jù)、土地利用數(shù)據(jù)、土壤分類數(shù)據(jù)等資料。其中地形數(shù)字高程圖（DEM）來源于地理空間數(shù)據(jù)云，分辨率為30 m；土地利用數(shù)據(jù)選取2015年全國土地利用數(shù)據(jù)，分辨率為30 m；土壤數(shù)據(jù)來自HWSD 全球土壤數(shù)據(jù)庫1∶100 萬中國土壤數(shù)據(jù)集。此外，由于水庫位于流域出口位置且研究區(qū)域多為山區(qū)，人類活動對研究區(qū)域內(nèi)水文循環(huán)過程的影響較弱，因此在SWAT模型中可忽略不計。

根據(jù)DEM 圖提取坡向、坡度數(shù)據(jù)，生成王快水庫以上流域河網(wǎng)，并根據(jù)設(shè)定的臨界集水面積將王快水庫以上流域劃分為14 個子流域（圖4），其中阜平水文站位于5 號子流域的出口；對流域內(nèi)的土地利用/覆蓋數(shù)據(jù)進(jìn)行重分類，共分為耕地、林地、草地、水域、城鄉(xiāng)工礦居民用地、未利用地等6個種類；以HWSD土壤數(shù)據(jù)集為基礎(chǔ)，結(jié)合SPAW軟件計算土壤層有效持水量、土壤濕密度等參數(shù)，完成土壤數(shù)據(jù)庫的制作并對土壤進(jìn)行重分類，共分為26 個種類，其中土地利用分類和土壤分類見圖5。根據(jù)土地利用、土壤、坡向設(shè)置，共劃分水文響應(yīng)單元（HRU）119個。

圖4 子流域劃分Fig.4 Watershed Subdivision of study area

圖5 土地利用分類和土壤分類Fig.5 The classification of land use and soil

2.2 參數(shù)敏感性分析

在模型初步構(gòu)建完成后，應(yīng)當(dāng)在其運(yùn)行結(jié)果的基礎(chǔ)上進(jìn)行參數(shù)敏感性分析，從而確定該模型中對徑流模擬影響較大的參數(shù)，以此為依據(jù)求解使模型模擬效果最佳的參數(shù)組合。本文依托SWAT-CUP 軟件的SUFI2算法進(jìn)行參數(shù)敏感性分析，參數(shù)的相對顯著性可由t-Stat檢驗(yàn)值確定，其中t-Stat檢驗(yàn)值的絕對值越大，參數(shù)敏感性越高［20］。選取影響徑流的28 個參數(shù)［21］，將模型輸出結(jié)果和阜平水文站的實(shí)測徑流數(shù)據(jù)結(jié)合，進(jìn)行1 000 次迭代運(yùn)算，挑選敏感性排名前17 的參數(shù)進(jìn)行下一步的率定，敏感性分析見表1。

表1 參數(shù)敏感性分析Tab.1 Parameter sensitivity analysis

由表1 可得：各參數(shù)對徑流模擬結(jié)果均有不同程度的相關(guān)性，其中最為敏感的參數(shù)是濕潤條件下的初始SCS 徑流曲線數(shù)，該參數(shù)的大小代表著下墊面不透水性的強(qiáng)弱，即CN2 值越大，徑流量越大；第二敏感的參數(shù)為地下水延遲時間，反映的是基流快慢，對水文過程有重要影響；此外，土壤飽和容重、最大冠層截留量等參數(shù)對模型輸出結(jié)果的影響也較為顯著。

2.3 參數(shù)率定與驗(yàn)證

將2008，2009 作為預(yù)熱年份，選用阜平水文站2010.1.1-2012.12.31 的實(shí)測流量對模型進(jìn)行參數(shù)率定，本文選取決定系數(shù)（R2）、納什效率系數(shù)（NSE）、相對誤差（Re）3 個指標(biāo)用于評價模型的模擬精度，R2和NSE越接近1，Re越接近0，模擬效果越好。一般來說，當(dāng)模擬結(jié)果滿足R2＞0.6，NSE＞0.5，|Re|≤20%時視為模擬結(jié)果可信；當(dāng)模擬結(jié)果滿足R2＞0.7，NSE＞0.65 時視為模擬結(jié)果良好［22］。計算公式分別如下：

式中：Qk為第k個實(shí)測流量；Sk為第k個模擬流量為實(shí)測流量的均值為模擬流量的均值。

在參數(shù)敏感性分析的基礎(chǔ)上，進(jìn)行多次迭代，使模擬結(jié)果與實(shí)測結(jié)果較為擬合，此時月徑流R2為0.92，NSE為0.90，Re為-9.1%。選取2013.1.1-2015.12.31 的實(shí)測流量數(shù)據(jù)作為驗(yàn)證資料，將參數(shù)組代入SWAT模型進(jìn)行模擬驗(yàn)證，模擬結(jié)果與實(shí)測徑流擬合較好，此時月徑流R2為0.85，NSE為0.81，Re為8.6%，精度達(dá)到模型的良好標(biāo)準(zhǔn)，率定期和驗(yàn)證期的模擬結(jié)果如圖6。

圖6 阜平站實(shí)測月徑流與模擬徑流對比Fig.6 Comparison of measured monthly runoff and simulated runoff at Fuping Station

3 基于降水隨機(jī)模擬的徑流預(yù)測

3.1年際-年內(nèi)降水統(tǒng)計分析

由圖7 可得：1958-2017年間，阜平站年降水量在均值線附近上下波動，21世紀(jì)以來，年降水量的波動幅度逐漸趨于平緩。汛期降水量占全年降水量比重大多為80%～90%，且年降水量過大或過小時，汛期降水占比相比均值存在較大的偏差，即年降水量的量級和降水年內(nèi)分布的規(guī)律存在一定的相關(guān)性。因此針對不同量級年降水量進(jìn)行年內(nèi)分布模擬更準(zhǔn)確，根據(jù)多年降水量均值和均方差將年降水量劃分為3 個狀態(tài)，劃分標(biāo)準(zhǔn)見表2。

圖7 1958-2017年阜平站年降水量Fig.7 Annual precipitation at Fuping Station from 1958 to 2017

為合理預(yù)測未來不同量級年降水量下的徑流量，設(shè)置300 mm 為某一枯水年的年降水量，550 mm 為某一平水年的年降水量，800 mm 為某一豐水年的年降水量。根據(jù)狀態(tài)I內(nèi)22年的逐日降水量進(jìn)行300 mm年降水量的年內(nèi)分布模擬，同理根據(jù)狀態(tài)II、狀態(tài)III 內(nèi)逐日降水量分別進(jìn)行550、800 mm年降水量的年內(nèi)分布模擬。

表2年降水量狀態(tài)劃分標(biāo)準(zhǔn)Tab.2 Standard for dividing annual precipitation status

3.2 降水隨機(jī)模擬

根據(jù)阜平站1958-2017年的逐日降水?dāng)?shù)據(jù)，分別計算各狀態(tài)下各月降水日的出現(xiàn)頻率，結(jié)合式（1）、（2）估算降水狀態(tài)轉(zhuǎn)移概率，其中各月P（W/D）見圖8。

圖8 各月降水狀態(tài)轉(zhuǎn)移概率P（W/D）Fig.8 Probability of monthly precipitation state transition P（W/D）

從圖8 中可以看出P（W/D）總體上均呈現(xiàn)先增大后減小的趨勢，且在汛期（5-9月）較大，在非汛期（10-4月）較小。當(dāng)年降水量為狀態(tài)III 時，P（W/D）最大值出現(xiàn)在8月；當(dāng)年降水量為狀態(tài)I 或II 時，P（W/D）最大值出現(xiàn)在7月。以各狀態(tài)1月、7月、12月降水量為例，其Gamma 分布概率Q-Q 圖見圖9，同理可得不同狀態(tài)年降水量的各月降水基本符合Gamma分布，滿足構(gòu)建耦合Markov 鏈-Gamma 分布隨機(jī)模擬模型的要求，并計算相應(yīng)的分布參數(shù)見圖10。

圖9 各狀態(tài)1、7、12月降水量Gamma分布Q-Q圖Fig.9 Q-Q chart of Gamma distribution of precipitation in January，July and December under different states

圖10 不同狀態(tài)年降水量各月Gamma分布參數(shù)Fig.10 Gamma distribution parameters of annual precipitation in different states

以800 mm年降水量的年內(nèi)分布模擬為例，耦合狀態(tài)III 中各月的降水轉(zhuǎn)移概率和Gamma 分布參數(shù)，根據(jù)生成的隨機(jī)數(shù)，連續(xù)進(jìn)行100年的模擬，綜合考慮年降水總量接近800 mm、汛期降水量占比和各月降水天數(shù)與樣本均值相近等條件選取模擬年，對其進(jìn)行適當(dāng)修正后即可得800 mm年降水量在日尺度上的分配。同理，求得550 mm 和300 mm年降水量的年內(nèi)分配見圖11。

圖11 不同量級年降水量年內(nèi)分配Fig.11 Annual distribution of precipitation in different magnitudes

3.3 不同量級降水下的徑流預(yù)測

將模擬的日降水輸入到SWAT 模型中，結(jié)合已驗(yàn)證的參數(shù)進(jìn)行輸出，即可得出3種降水情景下的月徑流過程線見圖12。

由圖12 可得，年降水量為800、550、300 mm 的隨機(jī)降水模擬結(jié)果輸出的阜平水文站月徑流過程中，三者均呈現(xiàn)減小-增大-減小的變化趨勢，且1 至3月的徑流量基本相等，這是由于汛期降水量占了年降水量較大的比重，不同量級的年降水量在1 至3月份的降水分配相差不大；年降水量為800 mm 相對應(yīng)的月徑流過程中，月徑流量在8月份達(dá)到最大，為17.81 m3/s，全年平均徑流量為5.97 m3/s；年降水量為550 mm 相對應(yīng)的月徑流過程中，月徑流量在9月份達(dá)到最大，為7.01 m3/s，全年平均徑流量為2.79 m3/s；年降水量為300 mm 相對應(yīng)的月徑流過程中，月徑流量在9月份達(dá)到最大，為4.08 m3/s，全年平均徑流量為1.51 m3/s。

圖12 不同量級年降水量下的徑流模擬Fig.12 Runoff simulation under annual precipitation of different magnitudes

4 結(jié) 論

基于2015年土地利用數(shù)據(jù)構(gòu)建相應(yīng)的SWAT 模型，可得CN2、GW_DELAY、SOL_BD等參數(shù)對徑流模擬結(jié)果的影響較為顯著；根據(jù)流域內(nèi)阜平水文站2008-2015年的實(shí)測徑流數(shù)據(jù)對模型進(jìn)行率定和驗(yàn)證，率定期和驗(yàn)證期的R2和NSE均高于0.8，|Re|均小于10%，表明該SWAT 模型對于研究區(qū)域徑流的模擬結(jié)果良好。

根據(jù)阜平站1958-2017年逐日降水?dāng)?shù)據(jù)可得，不同降水量級對應(yīng)的年內(nèi)分布規(guī)律存在一定差異。在此基礎(chǔ)上，分別構(gòu)建了耦合Markov-Gamma 分布的降水隨機(jī)模擬模型，并生成了年降水量為800、550、300 mm 對應(yīng)的日降水過程，輸入模型后得三種降水情景下的年徑流量分別為5.97、2.79、1.51 m3/s，最大月徑流量分別為17.81、7.01、4.08 m3/s。這一結(jié)論體現(xiàn)了不同量級降水對徑流過程的影響，可為后續(xù)降雨-徑流關(guān)系的研究提供理論基礎(chǔ)。此外，耦合Markov-Gamma 分布的降水隨機(jī)模擬模型沒有充分考慮特大暴雨等極端降水事件，從而在一定程度上影響徑流的模擬，因此對極端降水的發(fā)生過程進(jìn)行研究并和降水模擬相結(jié)合是一個研究方向。