基于SWAT模型與Copula修正的融雪徑流模擬

邢貞相，金超群，紀 毅，付 強，劉 東

（東北農(nóng)業(yè)大學(xué)水利與土木工程學(xué)院，哈爾濱 150030）

受溫度季節(jié)性變化影響，中高緯度地區(qū)積雪融化為典型自然現(xiàn)象，其產(chǎn)生的融雪徑流是融雪季節(jié)徑流重要組成部分，也是河流主要補給來源。在全球變暖背景下，全球范圍內(nèi)積雪正在逐漸消融[1]，積雪融水變化對當?shù)毓┧?、防洪、灌溉和雪崩等生態(tài)環(huán)境和生產(chǎn)生活具有重要影響[2]。因此，模擬融雪徑流對流域內(nèi)水資源持續(xù)利用具有現(xiàn)實意義[3]。

融雪過程受溫度、輻射、土地覆蓋等自然因素影響[4]。融雪徑流模擬采用兩類方法，一類方法基于溫度指數(shù)法，假設(shè)溫度是融雪主要影響因子，另一類是能量平衡法，考慮所在流域內(nèi)能量交換[5-6]。目前，全球已開發(fā)多個帶有融雪模塊分布式水文模型，如VIC、UEB、DHSVM和SWAT等，可對融雪徑流作精準模擬[7-8]。

溫度指數(shù)方法由于輸入?yún)?shù)較少、模型計算簡單等特點應(yīng)用廣泛[9-10]。Marks等研究表明，雨夾雪天氣產(chǎn)生融雪對空氣與積雪間能量交換最敏感[11]，溫度指數(shù)方法忽略風(fēng)速等自然因素作用，低估融雪前和融雪季節(jié)融雪。能量平衡方法在融雪模擬中綜合考慮多因素條件影響。當一定溫度與高濕度和風(fēng)速共同作用時，顯熱和潛熱成為能量主要來源，導(dǎo)致積雪開始融化[12]。Fontaine等認為能量平衡方法需額外氣象和地形參數(shù)輸入[13]。孟現(xiàn)勇等研究表明，在無額外數(shù)據(jù)要求情況下，能量平衡方法在融雪預(yù)測中表現(xiàn)優(yōu)于溫度指數(shù)方法[14]。Fuka等將不同形式能量平衡方法嵌入SWAT模型中[15]。Zhang等使用3種融雪模型應(yīng)用于大型山區(qū)流域（中國黃河上游水域）發(fā)現(xiàn)，能量平衡方法性能優(yōu)于溫度指數(shù)方法[16]。Debele等發(fā)現(xiàn)模擬3個不同流域（美國兩個小流域和中國黃河一個大流域）日徑流量方面[12]，能量平衡方法不如溫度指數(shù)方法。因模擬3個流域時，能量平衡方法中校準參數(shù)始終保持不變。針對不同流域環(huán)境應(yīng)分別校準。為提高融雪徑流模擬精度，Piani等提出利用數(shù)理統(tǒng)計方法預(yù)測模擬后誤差[17]，尋找模擬值與其誤差間相關(guān)性并修正模擬值。

為保證融雪徑流模擬精度，本文利用分布式SWAT模型模擬我國東北地區(qū)三江平原撓力河流域融雪期（3～5月）內(nèi)徑流，采用Copula函數(shù)預(yù)報徑流模擬誤差，修正徑流模擬值提高模擬精度，為流域春季農(nóng)業(yè)灌溉用水配置及春季防汛抗旱提供參考。

1 研究方法

1.1 SWAT融雪模塊

SWAT融雪模塊采用溫度指數(shù)法模擬融雪過程[13]。當日平均氣溫低于SWAT參數(shù)（SFTMP），將降水視為雪。當日平均氣溫高于SWAT參數(shù)（SMTMP）時，積雪開始融化。則雪水當量平衡計算方程為：

式中，SNOi和SNOi-1分別是第i天和第i-1天雪水當量；Rsf,i是第i天降雪量；Esub,i是蒸發(fā)量；SNOmlt,i是融雪量。

積雪溫度隨著氣溫變化而變化，主要受前1 d平均積雪溫度影響。前日積雪溫度對今日積雪溫度影響可用滯后系數(shù)（TIMP）描述。該參數(shù)用于解釋積雪密度、濕度和陽光。積雪溫度計算如下:

式中，Tsp,i和Tsp,i-1分別為第i天和第i-1天積雪溫度；Ta,i平均氣溫。則融雪量計算如下:

其中，Tmax為最高氣溫；bmlt,i為第i天融雪系數(shù)其計算如下:

式（4）中，SMFMX和SMFMN分別為最大值和最小值融雪因子。

1.2 二維Copula函數(shù)表達式

Copula函數(shù)克服傳統(tǒng)多變量聯(lián)合分布要求邊際分布為同一類型分布不足，通常水文模型模擬值與其誤差間具有相關(guān)關(guān)系，提高水文模型模擬精度，本文在率定期采用Copula函數(shù)描述模擬值和誤差值聯(lián)合概率分布，探討兩者關(guān)系，預(yù)報誤差修正模擬值。

在眾多Copula家族中，阿基米德Copula在水文研究應(yīng)用廣泛。阿基米德Copula中以Clayton Copula、Gumbel Copula和Frank Copula 3個函數(shù)應(yīng)用最廣。其中，F(xiàn)rank Copula結(jié)構(gòu)簡單，描述對稱相關(guān)結(jié)構(gòu)，對正負相關(guān)性均可適用且不限制相關(guān)性程度，而其他兩種Copula函數(shù)僅適用于正相關(guān)水文序列[18]。

本文研究融雪徑流及其模擬誤差序列關(guān)系不一定為正相關(guān)，因此選擇Frank Copula作為聯(lián)合分布函數(shù)，構(gòu)建融雪徑流及其模擬誤差聯(lián)合分布。

在Copula理論中，具有不同邊緣分布變量X和Y聯(lián)合分布函數(shù)表示為式（5）:

式中，H（x,y）為聯(lián)合分布函數(shù)；C（·，·）為連接2個邊緣概率分布Copula函數(shù)；F（x）和G（y）分別是x，y邊緣概率分布，若F、G連續(xù)，則C唯一。

Frank Copula函數(shù)作為Copula成員之一，最早由Frank于1979年提出，表達式為式（6）：

式中，u和v分別為變量邊緣累積概率；θ為參數(shù)，其取值范圍為-∞＜θ＜+∞，可由Kendall秩相關(guān)系數(shù)τ求得，而τ與θ表達式為式（7）和（8）為：

式中，（xi,yi）為樣本數(shù)據(jù)，sign為signum函數(shù)。

1.3 誤差修正方案

在傳統(tǒng)水文模擬中，利用多年歷史水文數(shù)據(jù)率定一組固定模型參數(shù)，代入水文模型，用于模擬驗證期水文過程。由于模型輸入誤差、結(jié)構(gòu)不完善及參數(shù)不確定性，導(dǎo)致模擬值存在不同程度誤差。

1.3.1 Copula修正

為降低模擬值誤差，本文采用Copula函數(shù)描述模擬徑流值與誤差值聯(lián)合概率分布，揭示兩者關(guān)系，并通過預(yù)報誤差以修正模擬值，減小模擬誤差，達到提高模擬精度目的，此方法在文中稱為Copula修正。

本文將2004～2008年作為率定期（2004年為模型預(yù)熱期），將該時期數(shù)據(jù)作為此率定模型參數(shù)值，并構(gòu)建模擬值及其誤差相關(guān)關(guān)系。具體步驟為：首先，將DEM、土地利用、土壤數(shù)據(jù)、氣象數(shù)據(jù)輸入SWAT模型中，模擬融雪期（3～5月）日徑流。計算徑流模擬值與實測值誤差，并確定誤差值與模擬值概率分布、計算兩變量間Kendall相關(guān)系數(shù)（τ），確定 Frank Copula參數(shù)（θ）值，建立Frank Copula聯(lián)合概率分布函數(shù)。將2009～2014年作為驗證期（2009年為模型預(yù)熱期），利用已率定模型參數(shù)，模擬融雪期徑流，采用Coupla修正方法修正模擬值。

1.3.2 優(yōu)化處理

采用滑動平均法平滑修正后得到最終模擬值，本文將該方法稱為優(yōu)化處理。優(yōu)化處理基于統(tǒng)計規(guī)律，將連續(xù)數(shù)據(jù)定義為一個長度固定為N序列，在新數(shù)據(jù)加入時，去掉該序列首個數(shù)據(jù)，其余N-1個數(shù)據(jù)依次前移，在第N個位置插入新數(shù)據(jù)，形成新序列，并對新序列作均值運算處理，將運算結(jié)果作為本次測量值，依次類推[19]。該方法在本文中主要用于修正后徑流過程線平滑。

2 應(yīng)用實例

2.1 研究區(qū)域

撓力河流域（如圖1）位于中國黑龍江省三江平原東北部（東經(jīng)45°43′～47°35′，北緯131°21′～134°10′），流域以東南部完達山為界，東接烏蘇里江，流域總面積2.42×104km2；撓力河流域地處寒溫帶大陸性濕潤季風(fēng)氣候區(qū)，年平均氣溫2.3～3.4℃。流域內(nèi)河流一般在11月中旬結(jié)冰，次年3月初解凍。年降雪量30～150 mm，其融雪徑流占全年徑流總量13.3%～24.9%。

圖1 撓力河流域概況Fig.1 General map of Naoli River Basin

2.2 數(shù)據(jù)

SWAT模型主要輸入數(shù)據(jù)包括數(shù)字高程模型（DEM）、土地利用、土壤屬性和日氣象數(shù)據(jù)（降水、最高和最低氣溫、相對濕度、風(fēng)速和太陽輻射）。本研究數(shù)據(jù)來源于網(wǎng)上公開數(shù)據(jù)，部分數(shù)據(jù)來自黑龍江省水文年鑒，來源如表1所示。

2009年頒布的《中華人民共和國郵政法》中為快件損壞賠償問題提供了法律依據(jù)。根據(jù)新郵政法的相關(guān)規(guī)定，普遍郵政服務(wù)業(yè)務(wù)范圍以外的郵件的損失賠償，適用有關(guān)民事法律的規(guī)定，并明確要求快遞企業(yè)應(yīng)當在其營業(yè)場所公示或者以其他方式公布其損失賠償辦法，以及用戶對其服務(wù)質(zhì)量的投訴辦法，以保障公眾的知情權(quán)?！翱爝f企業(yè)要加強對消費者合法權(quán)益的保護，今后可以采取保價、保險的方法解決賠償問題，即用戶根據(jù)寄遞物的貴重程度，選擇辦理相應(yīng)的保價業(yè)務(wù)，如發(fā)生丟失或損毀，由快遞企業(yè)依據(jù)保價金額進行賠償?！?/p>

表1 研究數(shù)據(jù)來源Table 1 Research data source

2.2.1 數(shù)字高程模型

數(shù)字高程模型（DEM）在流域和子流域所有地形屬性數(shù)據(jù)集推導(dǎo)過程中發(fā)揮重要作用，其屬性包括面積、坡度、坡度長度等，本文采用分辨率為30 m DEM作為輸入數(shù)據(jù)集。

SWAT輸入包括1個最小值-80 m和最大值831 m數(shù)字高程模型（DEM）。流域主要地形為平原，主要分布于流域中北部，坡度在0～5°。山區(qū)分布在流域西部和南部，大部分山區(qū)坡度在9°以上。

2.2.2 土地使用

采用2012年土地利用地圖，其原始空間分辨率為30 m。輕微修改原數(shù)據(jù)與SWAT數(shù)據(jù)庫一致，主要包括農(nóng)用地（80.17%）、林地（17.37%）和其他土地用途包括草、水、城市（2.46%）。

2.2.3 土壤類型與特征

根據(jù)土壤科學(xué)數(shù)據(jù)庫定義不同土壤類別，可以查詢我國主要土壤類型、面積、分類名稱和典型剖面分布情況，確定研究流域14種土壤類型，主要包括普通淋溶土（22.24%）、薄層黑土（18.41%）、松軟潛育土（17.89%），而其余11種類型土壤平均占總流域面積3.8%。

2.2.4 氣象和徑流數(shù)據(jù)

2.3 模型設(shè)置與配置

運用SWAT模型模擬研究流域內(nèi)水量平衡，DEM對流域劃分是模擬過程中第一步。在流域劃分后，定義水文響應(yīng)單元（HRUs），HRUs是地形、土地利用和土壤類型相似區(qū)域，下墊面特征相對單一和均勻，具有相似水文特征。因此，3個空間數(shù)據(jù)集（坡度、土地利用和土壤屬性）對于HRUs定義具有重要意義。將所需氣象變量輸入至模型，包括降雨量、最低和最高溫度、相對濕度、平均風(fēng)速和太陽輻射數(shù)據(jù)。應(yīng)用SWAT模型模擬。在模型啟動過程中，將2004～2008年作為率定期，確定適用于此流域參數(shù)值。而2009～2014年作為驗證期，評估模型適用性，其中2004、2009年則分別作為模型在相應(yīng)時段的熱身期，在模型模擬過程中僅模擬融雪期（3～5月）日徑流。

3 結(jié)果與分析

3.1 模型率定

SWAT-CUP廣泛用于SWAT校準計算機程序，作為SWAT操作重要補充，SWAT-CUP工具提供5種SWAT驗證方法，即GLUE、Para Sol、SUFI2、MCMC和PSO算法。李倩楠等比較5種不確定性評價方法，認為SUFI-2方法運行次數(shù)少、運算效率高[21]，故本研究采用SUFI2標定方法，并根據(jù)SWAT輸入/輸出文檔中參數(shù)特征、物理特征和參數(shù)范圍確定撓力河流域SWAT模型參數(shù)。選用參數(shù)及其描述、調(diào)整下限和上限見表2。

選用納什效率系數(shù)（NSE）、確定性系數(shù)（R2）、平均絕對誤差（MAE）及均方根誤差（RMSE）作為模型評價指標。利用NSE及R2作為本文主要評價方法，估計模型參數(shù)適用性，而MAE與RMSE對比相同時段不同模擬方法結(jié)果。NSE＞0.5，R2＞0.6，且MAE與RMSE均小于30，模擬效果顯著。本文率定期（2004～2008年）內(nèi)每年融雪期日徑流模擬值與實測值對比見圖2。模型評價結(jié)果見表3?？梢姡氏禂?shù)僅0.26，2008年融雪期徑流峰值小于其他年份，這是SWAT模擬徑流不準原因，但整體模擬效果良好，故SWAT模型在撓力河流域融雪徑流模擬中具有良好適用性。

表2 撓力河流域SWAT模型參數(shù)率定值Table 2 Calibration of SWAT model parameters in the Naoli River basin

圖2 率定期（2005～2008）徑流過程線Fig.2 Runoff process line in the calibration(2005-2008)

表3 率定期（2005～2008）模擬結(jié)果Table 3 Result of simulation in the calibration(2005-2008)

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

本文采用SWAT-CUP中全局敏感性分析[22]，利用t值和P值評價模型參數(shù)敏感性。參數(shù)敏感性與t值呈正相關(guān)、與P值呈負相關(guān)，結(jié)果如圖3所示，CN2和ALPHA_BF其t值均在4～5，其中CN2的t值大于ALPHA_BF，P值均接近于0，故在研究區(qū)域中參數(shù)敏感性CN2排名第一，ALPHA_BF排名第二。而SFTMP、SMTMP和TIMP其t值均在1.5～2.5，P值均在0～0.4。而其他參數(shù)t值均小于1.5，P值則在0.4以上。綜上所述，在12個參數(shù)中，CN2和ALPHA_BF是最敏感參數(shù)，其次為SFTMP、SMTMP和TIMP。

圖3 參數(shù)敏感性關(guān)系Fig.3 Relation diagram of parameter sensitivity

3.3 基于Copula驗證期模擬修正

由表3可知，率定期參數(shù)值（即由參數(shù)確定模型）對2006年徑流模擬精度較高，2008年較差。但整體上，率定期模型模擬精度良好。為提高個別年份及整體模擬精度，本文利用Frank Copula對模擬徑流值與誤差值建立聯(lián)合概率分布，用于融雪期模擬徑流Copula修正。由圖2可知，前期融雪期徑值流較低，模擬精度較高，修正必要性不大，因此本文僅對自融雪期開始第40天后徑流模擬值作Copula修正。利用非線性擬合技術(shù)確定模擬徑流值及誤差值概率分布（見圖4），其中圖4a、b分別為模擬值和誤差值概率分布。本文對模擬值分別采用韋伯分布、正態(tài)分布和伽馬分布作擬合嘗試，而對誤差值采用正態(tài)分布、極值分布和Logistic分布作擬合嘗試。擬合結(jié)果發(fā)現(xiàn)，正態(tài)分布對模擬值和誤差值擬合效果良好，其中模擬徑流值主要分布在30～60，而誤差主要分布在5～25。利用Cop?ula函數(shù)對兩者聯(lián)合概率分布描述見圖5。

圖4 模擬徑流值與誤差值概率密度分布Fig.4 Probability density distribution of simulated runoff values and error values

圖5 誤差V和模擬值U聯(lián)合概率分布Fig.5 Joint probability distribution of errors and simulation values

構(gòu)建SWAT模擬驗證期（2010～2014年）融雪期徑流，采用Copula修正并經(jīng)優(yōu)化處理，融雪徑流模擬結(jié)果如表4所示，模擬徑流過程如圖6所示。由表4可知，驗證期SWAT模擬納什效率系數(shù)值為0.75，采用Copula修正方法納什效率系數(shù)為0.92，利用移動平均法優(yōu)化處理納什效率系數(shù)為0.93。對于確定性系數(shù)、平均絕對誤差和均方根誤差，與傳統(tǒng)SWAT模型相比，經(jīng)Copula修正和優(yōu)化處理后驗證期整體模擬結(jié)果均有明顯提高，單一年份模擬結(jié)果也有所提高（2014年NSE提高139%）。然而，2012年Copula修正后較SWAT模擬結(jié)果有所降低（NSE降低22.5%），但優(yōu)化處理后模擬結(jié)果提高且與SWAT模擬相當。由2012年SWAT模擬徑流過程可知，SWAT模擬徑流值遠大于真實值。但在率定期，SWAT模擬值低于真實值，Copula修正后無理想結(jié)果，但優(yōu)化處理后可彌補。綜上所述，經(jīng)修正處理后，融雪期日徑流模擬精度有明顯改善。

表4 驗證期模型模擬結(jié)果、修正結(jié)果及優(yōu)化處理結(jié)果Table 4 Results of the simulation,correction and optimization in the validation

圖6 驗證期模型模擬、修正及優(yōu)化處理徑流過程Fig.6 Runoff process of the simulation,correction and optimization in the validation

4 討論與結(jié)論

通過敏感性分析發(fā)現(xiàn)，CN2（徑流曲線數(shù)）和ALPHA_BF（基流系數(shù)）為主要敏感參數(shù)（敏感程度分列1～2位），并與控制地表徑流和地下徑流密切相關(guān)。CN2取決于土地利用類型和土壤含水量，而ALPHA_BF主要控制地下徑流量，因SWAT模型模擬融雪徑流在計算過程中主要通過影響地表徑流與地下水補給控制河道徑流量。與融雪有關(guān)參數(shù)（SMTMP、SFTMP和TIMP）分別排在 3～5位。其中，SMTMP是融雪溫度，為模擬融雪徑流重要參數(shù)。在SWAT融雪模塊中，當平均氣溫達到SMT?MP值時，雪開始融化。SFTMP是降雪溫度。在融雪期前幾個月，由于氣溫較低，大量降水以雪形式儲存在流域中。評價融雪期前以雪形式儲存水量在徑流模擬中具有重要作用。因此，SFTMP在模型模擬中與SMTMP有重要作用。積雪溫度受前1 d積雪平均溫度影響極大，隨氣溫升降而變化。前日積雪溫度對今日積雪溫度影響可用滯后因子描述。因素由TIMP指定，TIMP是集總參數(shù)，綜合考慮積雪密度、含水量和日照等因素，靈敏度較高。在SWAT模型中，由于模擬融雪徑流采用溫度指數(shù)法，溫度相關(guān)參數(shù)靈敏度處于中等水平。

模擬結(jié)果分析發(fā)現(xiàn)，無論是率定期還是驗證期均采用同一套模型參數(shù)，對有些年份模擬精度較好，但個別年份模擬精度不理想，如2008和2014年，其納什效率系數(shù)分別為0.26和0.33。原因是在模擬過程中模型參數(shù)存在不確定性。

本文提出Copula修正和優(yōu)化處理并未改變SWAT模型融雪徑流計算方法，僅修正和處理模擬結(jié)果，提高模型模擬精度。通過圖2和6發(fā)現(xiàn)，SWAT模擬徑流值低于實測值，因SWAT融雪徑流模擬采用是溫度指數(shù)法，未考慮融雪期風(fēng)速對融雪影響。Copula修正和優(yōu)化處理方法可彌補SWAT模型不足。通過驗證期模擬徑流過程（圖6）發(fā)現(xiàn)，修正后徑流過程線并非平滑曲線，利用滑動平均法對Copula修正后徑流作平滑優(yōu)化處理。由表4發(fā)現(xiàn)，2010～2014年Copula修正后納什效率系數(shù)為0.92，優(yōu)化處理后納什效率系數(shù)為0.93，模擬精度提高。優(yōu)化處理方法可平滑Copula修正過程線，提高徑流模擬精度。