張梅潔， 呂海深， 劉 娣， 朱永華， 孫銘悅

（1.河海大學(xué)水文水資源與水利工程科學(xué)國家重點實驗室，江蘇 南京 210098；2.河海大學(xué)水文水資源學(xué)院，江蘇 南京 210098）

中國西北地區(qū)的河流多發(fā)源于山區(qū)，冰雪融水是河流春季的重要補給來源［1］。西營河位于甘肅省境內(nèi)，其徑流主要來源于祁連山東段融雪和山區(qū)降水。在西營河流域，多數(shù)研究側(cè)重于西營灌區(qū)的生態(tài)及水利建設(shè)，進行水文模擬極少，迫切需要探索適用于該區(qū)域的水文模型。在該流域范圍內(nèi)使用水文模型進行模擬時，需要考慮春季融雪的影響。Wu［2］等通過研究長序列融雪徑流過程，完成了融雪組分的評估，指出春夏季融雪徑流比例隨氣候變暖而有所改變，融雪過程會影響高寒山區(qū)水文過程。李志龍［3］在研究資料缺乏的寒區(qū)流域水文模擬中，將改進的新安江模型作為混合模型的一部分，進行了融雪徑流的模擬。目前，關(guān)于融雪徑流預(yù)測已經(jīng)有應(yīng)用較廣的模型，如SRM模型［4］；此外還有包含融雪模塊的模型，如SWAT 模型［5］、MIKE SHE 模型［6］等。然而，SRM 模型主要對于春季融雪模擬較好，多數(shù)只能應(yīng)用于年內(nèi)短期的徑流預(yù)報［7］。SWAT 模型與SHE模型又較為復(fù)雜，輸入的數(shù)據(jù)需要包含大量信息（包括數(shù)字高程數(shù)據(jù)、土地利用數(shù)據(jù)、積雪覆蓋特征、土壤特征、降水氣象數(shù)據(jù)等）［8］，在西營河這種數(shù)據(jù)較為缺乏的流域，數(shù)據(jù)搜集困難且數(shù)量較少，代表性不強。故而需要探索更加適用于西營河流域的融雪徑流模型。

新安江模型［9］作為國內(nèi)應(yīng)用較廣的模型之一，其模型結(jié)構(gòu)清晰、層次分明，在實際應(yīng)用中得到不斷地完善，可以通過改進的模型進行數(shù)據(jù)缺乏地區(qū)的水文模擬。劉金濤等［10］總結(jié)了新安江模型發(fā)展的特點，主要是架構(gòu)形式趨于多樣性、產(chǎn)匯流模型不斷物理化、應(yīng)用范圍不斷延伸。Lyu 等［11］采用粒子群算法和集合卡爾曼濾波數(shù)據(jù)同化的3 種變化，對羅江和資料匱乏地區(qū)進行模型參數(shù)批量估計。韓元元等［12］將變動態(tài)存儲系數(shù)法結(jié)合新安江模型進行河道洪水演算。以上研究說明在模型驅(qū)動數(shù)據(jù)及方法上均可與新安江模型進行配合，從而提高模型的應(yīng)用范圍及模擬效果。Alazzy［13］采用不確定性估計（GLUE）方法，對中國新安江模型估算月徑流的不確定性進行檢驗和預(yù)測，發(fā)現(xiàn)采用Nash效率系數(shù)評估的新安江模型的參數(shù)不確定性更小。據(jù)此，本文率定參數(shù)時，主要以Nash 效率系數(shù)作為率定參數(shù)時的目標(biāo)函數(shù)。鄧元倩等［14］、邵成國等［15］分別將新安江模型應(yīng)用于灃河流域、新疆烏魯木齊河，模擬效果較好，說明新安江模型適用流域逐漸擴大到西北地區(qū)，不再僅僅局限于濕潤地區(qū)。趙蘭蘭等［16］根據(jù)干旱地區(qū)實際情況，將超滲產(chǎn)流模式加入新安江模型在拒馬河流域進行了暴雨洪水預(yù)報及對比分析，在新安江模型各層次內(nèi)部使用的方法可根據(jù)流域的地理位置等條件進行適當(dāng)改進。至此，超滲產(chǎn)流與新安江模型的融合，使新安江模型突破了原有的固化模式，從最根本的產(chǎn)流層次進行改進，加強了新安江模型在全國多地的適用性。姜卉芳［17］最早將融雪徑流模擬加入新安江模型應(yīng)用于切德克流域，隨后田龍［18］等人在該版本基礎(chǔ)上考慮能量概念進行高寒山區(qū)融雪徑流模擬。然而，該模型仍需對太陽輻射、地形影響、積雪面積等進行大量分析，在資料缺乏的山區(qū)應(yīng)用較難。本文選取應(yīng)用廣泛、層次結(jié)構(gòu)化清晰的新安江模型，在其中加入僅以氣溫為主要閾值的融雪模塊，從而擴大新安江模型的適用范圍。

在資料匱乏地區(qū)，如西營河九條嶺站以上流域無降水觀測，且流域內(nèi)地形差異大，九條嶺水文站降水無法代表全流域情況。較為準(zhǔn)確的降水驅(qū)動數(shù)據(jù)對水文模擬具有重要意義，直接影響水文模擬精度。在模擬時，需要選取最適合的驅(qū)動數(shù)據(jù)。孫銘悅［19］等在進行水文模型模擬時，比較了站點數(shù)據(jù)和格點降水?dāng)?shù)據(jù)的模擬效果，得出格點降水比站點平均的降水?dāng)?shù)據(jù)更好。本文首先進行格點降水?dāng)?shù)據(jù)修正，隨后將融雪模塊與新安江模型耦合，利用修正后的降水?dāng)?shù)據(jù)驅(qū)動耦合模型進行徑流模擬。

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

1.1 研究區(qū)概況

西營河位于武威市西南部，是石羊河流域的一條支流，上游由寧昌河和水關(guān)河組成，發(fā)源于祁連山東段冷龍嶺。西營河屬于高寒半干旱半濕潤氣候區(qū)。徑流主要由高山融雪以及祁連山區(qū)降水組成。西營河多年平均出山徑流達3.184×108m3，位列石羊河上游八條河流之首，占比22.3%。出山口附近設(shè)有九條嶺水文站，也是該流域出山口以上唯一的水文測站。九條嶺站以上區(qū)域植被良好，水土流失輕微，此處與上游由于地勢的差異，降水差別較大［20］。上游年降水量可達500～800 mm，九條嶺站年降水量僅為200～400 mm。流域4—5月多為春季融雪徑流，6—9月徑流為山區(qū)降水產(chǎn)生。本文主要研究區(qū)為西營河九條嶺站以上流域（圖1），集水面積為1077 km2。該站是石羊河流域的八大河流中不受人工設(shè)施影響的水文測站之一，積累了長序列水文實測資料。

圖1 研究區(qū)域示意圖Fig.1 Schematic diagram of the study area

1.2 數(shù)據(jù)來源

本文改進的新安江模型輸入的降水和氣溫數(shù)據(jù)采用國家氣象信息中心發(fā)布的數(shù)據(jù)集。降水采用中國地面降水日值0.5°×0.5°格點數(shù)據(jù)集（V 2.0），代碼為SURF_CLI_CHN_PRE_DAY_GRID_0.5。氣溫采用中國地面氣溫日值0.5°×0.5°格點數(shù)據(jù)集（V2.0），代碼為SURF_CLI_CHN_TEM_DAY_GRID_0.5。此外，九條嶺水文站2011—2018 年的逐日蒸發(fā)和逐日流量資料來源于中華人民共和國水文年鑒第10卷第5冊甘肅河西地區(qū)內(nèi)陸河水文資料。

1.3 研究方法

1.3.1 格點降水?dāng)?shù)據(jù)修正方法 格點數(shù)據(jù)集主要由中國地面高密度臺站的降水結(jié)合薄盤樣條法制作，依托于站點數(shù)據(jù)，插值時存在一定誤差。數(shù)據(jù)集（V 2.0）制作人趙煜飛［21］對該數(shù)據(jù)集進行精度評估時指出，夏季誤差較其他季節(jié)較為明顯。將格點降水?dāng)?shù)據(jù)用于降水（融雪）徑流模擬時發(fā)現(xiàn)在夏季會出現(xiàn)少量徑流異常值，需要進行降水?dāng)?shù)據(jù)修正。經(jīng)試驗，未經(jīng)修正的格點降水，在研究區(qū)每年都會出現(xiàn)少數(shù)的異常情況，徑流未出現(xiàn)峰值，然而作為徑流主要成因的降水卻出現(xiàn)較大值，且與序列前后日突變較大，與臨近站點降水及相關(guān)記錄比較發(fā)現(xiàn)此日并非暴雨日且前后日并無較大降雨，從而考慮應(yīng)是數(shù)據(jù)集制作產(chǎn)生的誤差導(dǎo)致。常用的修正方式是與附近站點序列進行比較。然而，在西營河流域，附近的氣象站點距離流域較遠（圖2），且分布多偏東北、西北方向，作為降水產(chǎn)流主要來源的西南部山區(qū)沒有代表性測站，無法全面展現(xiàn)整個流域降水、氣溫等情況。

圖2 西營河流域鄰近氣象站點分布Fig.2 Distribution of meteorological stations near Xiying River Basin

本文探索性嘗試通過年最大降水和徑流系數(shù)進行約束，篩選出每年的異常格點降水值進行修正，從而提高徑流模擬的精度。經(jīng)過多年格點降水徑流模擬試驗，歸納得出以下降水量和徑流系數(shù)兩個約束能檢測出較明顯的格點降水異常數(shù)據(jù)。同時違背降水量約束（公式1a）和徑流系數(shù)約束（公式1b）兩個約束條件判斷為異常點：

式中：R表示出口九條嶺水文站日尺度實測徑流深（mm）；P為逐日格點降水值（mm）；α為徑流系數(shù)，徑流系數(shù)0.3 的約束條件主要是采用3σ原則，通過正態(tài)分布置信區(qū)間算出，即平均年徑流系數(shù)為0.503，標(biāo)準(zhǔn)差σ為0.067，在小于μ-3σ即0.3 時，概率僅為0.27%，認為該測量值為壞值。

盡管下墊面條件和氣候條件的復(fù)雜性會對不同自然地理條件的降水-徑流關(guān)系產(chǎn)生一定影響，降雨和徑流仍具有較強的成因關(guān)系，可以作為兩個序列的初步評定參考。本研究識別出異常點后，利用降水和徑流相關(guān)關(guān)系進行修正，將識別出的異常點修正到每年繪制的剔除異常點的降水-徑流相關(guān)關(guān)系圖上，具體流程如圖3所示。

圖3 格點降水異常點識別與修正流程Fig.3 Grid precipitation anomaly identification and correction process

某些年份單日格點降水最大值達到40 mm甚至更大時，若此點為異常點，會影響降水量約束的異常點挑選，檢測時遇此情況，需要以年第二大降水值作為降水量約束的基礎(chǔ)值，再乘以試驗率定出的比例系數(shù)0.7。

1.3.2 融雪模塊與新安江模型耦合 本文在新安江模型基礎(chǔ)上，增加了融雪模塊，主要包括融雪模塊、蒸散發(fā)計算、蓄滿產(chǎn)流計算、三水源劃分和匯流計算。

（1）融雪模塊

改進后的新安江模型在第一層次前，加入融雪模塊，主要采用度日因子法進行融雪計算。度日法假設(shè)融雪量與氣溫為顯著線性關(guān)系，其計算公式如下：

式中：M為融雪水量（mm·d-1）；d為度日因子（mm·℃-1·d-1）；T0為融雪臨界溫度（℃）；T為氣溫（℃）。

在融雪模塊中，主要輸入數(shù)據(jù)為氣溫和降水。積雪消融的能量主要來源于凈輻射和感熱通量，其中，凈輻射包含長波輻射和短波輻射。冰雪消融以近地表層的長波輻射和感熱通量原因為主，氣溫又是地表長波輻射的主要影響因素，從而本方法將能量的影響直接由較易獲取的氣溫作為模型的驅(qū)動。加入融雪模塊后，將融雪水量和降雨量綜合的液態(tài)水作為后續(xù)模型的輸入。

（2）蒸散發(fā)計算

模型采用三層蒸發(fā)模式，其中，蒸散發(fā)能力使用蒸發(fā)折算系數(shù)與九條嶺水文站觀測值的乘積。

式中：EP為蒸散發(fā)能力（mm）；KC為蒸散發(fā)折算系數(shù)；Epan為蒸發(fā)皿觀測蒸發(fā)（mm）。運用的主要思想是土壤蒸發(fā)能力與相同氣象條件下的水面蒸發(fā)一般呈線性關(guān)系。

（3）蓄滿產(chǎn)流

蓄滿產(chǎn)流模式為包氣帶達到田間持水量的產(chǎn)流模式。概化為達到田間持水量后，才開始產(chǎn)流。由于流域各處蓄水能力不同，降水與初始包氣帶蓄水量之和小于流域包氣帶蓄水容量最大值時，為流域部分面積產(chǎn)流。同理，降水與初始包氣帶蓄水量之和大于流域包氣帶蓄水容量最大值時為全流域產(chǎn)流。具體計算公式如下：

式中：A為流域初始包氣帶蓄水量所對應(yīng)的縱坐標(biāo)值（mm）；WMM為流域單點最大蓄水容量（mm）；WM為流域任一點的包氣帶蓄水容量（mm）；B反映流域包氣帶蓄水容量分布的不均勻性，B越小表示越均勻；R為產(chǎn)流量（mm）；W為流域蓄水量（mm）。值得注意的是，PE中的P指的是融雪計算后的液態(tài)水量（mm）。

（4）三水源劃分

產(chǎn)流量的劃分延續(xù)傳統(tǒng)新安江模型的思路，分為地表徑流、壤中流和地下徑流。考慮產(chǎn)流面積的影響和下墊面條件，采用與上一模塊包氣帶蓄水容量分布不均類似的處理，仍將產(chǎn)流面積上的自由水容量分布用拋物線近似。其中，自由水是毛管水和重力水的合稱，是在土壤中受力的作用可以自由流動的水。

（5）匯流計算

地面徑流直接流入河網(wǎng)，計算公式如下：

式中：QS為地面徑流匯流（m3·s-1）；RS為地面徑流（mm）；U為單位轉(zhuǎn)換系數(shù)，為常數(shù)；i指日數(shù)，時間尺度為日尺度。

壤中流匯流采用線性水庫法，計算公式如下：

式中：QI為壤中流匯流（m3·s-1），CI為壤中流消退系數(shù)，RI為壤中流徑流（mm），U 為單位轉(zhuǎn)換系數(shù)，為常數(shù)。

地下徑流匯流仍采用線性水庫法，采用公式如下：

式中：QG為地下徑流匯流（m3·s-1）；CG為地下水消退系數(shù)；RG為地下水徑流（mm）；U 為單位轉(zhuǎn)換系數(shù)，為常數(shù)。

1.3.3 評估準(zhǔn)則 根據(jù)水文預(yù)報評價標(biāo)準(zhǔn)，采用相對誤差RE和確定性系數(shù)DC進行評價。二者公式如下：

2 結(jié)果與分析

2.1 格點降水?dāng)?shù)據(jù)修正結(jié)果分析

根據(jù)流域在降水網(wǎng)格的位置，提取流域內(nèi)格點降水。經(jīng)過格點降水異常點識別與修正方法進行修正（圖3）。在約束條件所構(gòu)成的區(qū)域內(nèi)尋找異常點，每年需要修正的格點并不多，但這些格點屬于偏離降水徑流成因關(guān)系曲線較大的降水?dāng)?shù)據(jù)點。由2018 年的格點降水的修正可以看出（圖4），2018年僅有兩日的格點降水需要修正，即紅色標(biāo)注部分（圖4a），將異常點處的降水?dāng)?shù)據(jù)修正到剔除異常點后的降水-徑流關(guān)系曲線上（圖4b），這兩個異常點的降水值分別由24.3 mm、21.8 mm修正為0.6 mm和4.3 mm。修正后，率定期的模擬效果確定性系數(shù)從0.59 提升到0.63，提升了0.04。驗證期的模擬效果確定性系數(shù)從0.69 提升到0.72，提升了0.03（表1）。2013 年的格點降水修正對確定性系數(shù)的提升較為明顯，從0.23 提升到了0.69。格點降水修正對于75%的年份均表現(xiàn)為模擬效果提升。

圖4 2018年格點降水的修正Fig.4 Grid precipitation correction in 2018

2.2 融雪徑流模擬分析

分別用融雪模型融合前后及降水修正后的新安江模型進行模擬，模擬效果見表1。由表1 可知，融雪模塊耦合后的模擬結(jié)果通常比傳統(tǒng)的新安江模型模擬得更好?？梢宰⒁獾剑倪M過程中，確定性系數(shù)的提高和相對誤差的降低并不同步，這是由于在參數(shù)優(yōu)化時，選擇了確定性系數(shù)作為目標(biāo)函數(shù)，其更能反應(yīng)模擬過程，體現(xiàn)流量模擬值的變化在總體變化中的比例。而相對誤差更多地評價模擬總流量與實測值之間的關(guān)系，重點在總量的控制上。除2013 年由于以確定性系數(shù)為目標(biāo)函數(shù)的參數(shù)優(yōu)化使得確定性系數(shù)大幅提升，相對誤差為23%，超過了20%外，其余年份的各種情況模擬，相對誤差均在20%以下，在后續(xù)研究中若能繼續(xù)進行多目標(biāo)參數(shù)優(yōu)化，有望將2013年的相對誤差進一步降低。

表1 不同情況下新安江模型模擬結(jié)果Tab.1 Simulation results of Xin’anjiang model under different conditions

選擇2011—2016 年作為模型率定期，2017—2018年作為模型驗證期。率定期參數(shù)優(yōu)化后，確定性系數(shù)達到0.63，此套參數(shù)用于驗證期，確定性系數(shù)為0.72。說明此模型在西營河流域具有一定適用性。經(jīng)過融雪模塊融合、格點降水修正，87%以上的模擬年確定性系數(shù)在0.6以上、62%以上的年份模擬達0.69以上，模擬結(jié)果較好。新安江模型平均確定性系數(shù)為0.62，融雪耦合和格點降水修正后平均確定性系數(shù)為0.67，提高0.05。將變動態(tài)存儲系數(shù)法用于新安江模型河道洪水演算的研究在印江河流域只提高了0.0182［12］；改進新安江模型在烏江獨木河流域模擬時確定性系數(shù)僅提高0.02［22］。與類似研究相比，本文運用改進的新安江模型在西營河流域的平均確定性系數(shù)提高0.05，效果顯著。

采用2011—2016年進行參數(shù)優(yōu)化模型率定，利用2017—2018年進行模型的驗證，參數(shù)率定時采用遺傳算法優(yōu)化參數(shù)，其中以確定性系數(shù)作為算法的目標(biāo)函數(shù)進行率定，表2顯示了優(yōu)化后的模型參數(shù)。

表2 含融雪的新安江模型參數(shù)及物理意義Tab.2 The parameters of snowmelt-Xin’anjiang model and its physical significance

耦合融雪的新安江模型率定期和驗證期模擬結(jié)果如圖5所示，除第二年洪峰模擬偏低外，其余年份洪峰模擬效果均很好。率定期的多年模擬結(jié)果中，確定性系數(shù)達到0.63，比未進行格點降水修正的0.59 提高了0.04。在驗證期，可以看到模型已模擬出春季融雪徑流，夏季的洪水部分前期模擬較好，后期模擬效果較低，可能與當(dāng)年的格點降水?dāng)?shù)據(jù)仍存在誤差有關(guān)。驗證期耦合融雪的新安江模型經(jīng)降水修正后比新安江模型模擬的確定性系數(shù)提升0.05。

圖5 模型率定期和驗證期模擬結(jié)果Fig.5 Simulation results during model calibration period and validation period

在融雪徑流的模擬過程中，率定期內(nèi)，2011年、2013 年、2014 年、2015 年的春季徑流模擬很好。在未出現(xiàn)較大降水時，也會出現(xiàn)春季小洪峰，此時主要是由于流域上游祁連山東段的山間積雪融化導(dǎo)致。如2011 年4 月20 日，此前4 d 及當(dāng)天未出現(xiàn)降水，出現(xiàn)了13.8 m3·s-1的流量峰值。由于當(dāng)時氣溫逐漸升高，山上的積雪逐漸融化，形成了此次春季非降雨成因的融雪徑流。改進的新安江模型在夏季洪峰的模擬上效果很好，這可能是因為在西營河流域，產(chǎn)流層較薄，多數(shù)情況下當(dāng)日降水會產(chǎn)生當(dāng)日出流，流域匯流延時較短。在驗證期，總體模擬效果普遍不如率定期，但春季徑流仍然能夠模擬較好，說明模型加入的融雪模塊運用較好，以氣溫為基礎(chǔ)足以計算出研究區(qū)的春季徑流。2017 年的夏季洪水模擬較好，而2018 年7 月、8 月的徑流，雖然趨勢一致，但多出現(xiàn)模擬值較高的情況。對照降水?dāng)?shù)據(jù)可以發(fā)現(xiàn)，在2018 年修正的降水僅為9 月降水。這是由于在降水修正時，沒有將非較大值的降水考慮在修正范圍內(nèi)，修正的兩個約束并不能識別較小降水的趨勢性問題，這也是在后續(xù)降水?dāng)?shù)據(jù)使用中，需要深入探討的問題。

3 結(jié)論

本文考慮到研究區(qū)實際情況將新安江模型加入融雪模塊，考慮融雪對徑流的影響。另一方面，由于研究區(qū)水文站點少，采用格點降水和氣溫數(shù)據(jù)集作為改進的新安江模型的驅(qū)動數(shù)據(jù)。運用過程中發(fā)現(xiàn)異常點并進行了識別方法及修正方法的探索，將格點降水進行修正作為模型驅(qū)動，進行水文模擬，進而分析新安江模型改進前后、格點降水?dāng)?shù)據(jù)修正前后的模擬效果，得到如下結(jié)論：

（1）在站點資料缺乏的干旱半干旱區(qū)，利用全國范圍內(nèi)制作的氣象數(shù)據(jù)集解決由于資料缺乏，數(shù)據(jù)集存在一定誤差的問題，運用降水-徑流關(guān)系曲線修正少數(shù)偏差較大的格點降水后，能夠提高數(shù)據(jù)可用性，模擬的確定性系數(shù)在8 a 年中有6 a 均有提高，其余2 a差異不大。

（2）將融雪模塊耦合新安江模型，使用的模型驅(qū)動數(shù)據(jù)較少，運用簡單。新安江模型在西營河流域模擬流量過程的年確定性系數(shù)均在0.5 以上，經(jīng)過融雪耦合和格點降水修正后87%以上的模擬年確定性系數(shù)0.6以上、62%以上的模擬年確定性系數(shù)達0.69 以上。新安江模型平均確定性系數(shù)為0.62，融雪耦合和格點降水修正后平均確定性系數(shù)為0.67，提高了8%。總體模擬效果較好，在西營河流域運用良好。