劉江濤 ，徐宗學(xué) ，趙 煥 ，彭定志

（1.北京師范大學(xué) 水科學(xué)研究院，北京 100875；2.城市水循環(huán)與海綿城市技術(shù)北京市重點(diǎn)實(shí)驗(yàn)室，北京 100875）

1 研究背景

積雪在地球的物質(zhì)循環(huán)和水量平衡中扮演著十分重要的角色，積雪能夠增加地表反照率，影響能量的收支平衡［1-2］。積雪融化產(chǎn)生的融雪徑流是包括黃河、長(zhǎng)江和雅魯藏布江在內(nèi)的眾多河流的主要補(bǔ)給源，為下游數(shù)以億計(jì)人民的生產(chǎn)和生活提供十分重要的支撐作用。在全球變暖背景下，全球范圍內(nèi)積雪正在逐漸融化［3］。尤其是在北半球的中高緯度地區(qū)，積雪融水的變化更易給當(dāng)?shù)氐纳鷳B(tài)環(huán)境和生產(chǎn)生活帶來影響，雪消融會(huì)引發(fā)洪澇及泥石流等自然災(zāi)害［4］。因此，研究融雪徑流動(dòng)態(tài)變化尤其是在融雪期的變化對(duì)高原河流流經(jīng)地區(qū)的水資源持續(xù)利用和防災(zāi)減災(zāi)具有十分重要的現(xiàn)實(shí)意義［5］。

融雪徑流的模擬主要有基于物理機(jī)制的能量平衡法和度日因子法，基于物理機(jī)制的模型能夠反映流域中積雪在消融過程中能量的動(dòng)態(tài)變化特征，比如積雪表面所接受到的長(zhǎng)波輻射、短波輻射以及積雪內(nèi)部的熱通量交換等能量變化細(xì)節(jié)［6］?；谖锢頇C(jī)制的模型有SHE模型［7］、PRMS模型［8］、VIC模型［9］等，它們?cè)诒姸嗟谋└采w區(qū)得到了廣泛的應(yīng)用［10-12］。然而，基于物理機(jī)制的模型所需的數(shù)據(jù)較多，我國(guó)的西部高寒地區(qū)，氣象站點(diǎn)分布極度稀缺，水文氣象數(shù)據(jù)無論是在數(shù)量上還是時(shí)間序列的長(zhǎng)度上都較難滿足模型的輸入需求。度日因子模型假設(shè)融雪量與氣溫存在線性關(guān)系，模型的結(jié)構(gòu)簡(jiǎn)單，所需的數(shù)據(jù)容易獲得，并且模擬的結(jié)果較為理想。常見的度日因子模型有SRM模型［13］、HBV模型［14］、WetSpa模型［15］等，國(guó)內(nèi)外研究學(xué)者利用度日因子模型或改進(jìn)的度日因子模型在融雪過程模擬中開展了廣泛的研究。李蘭海等［16］利用APHRODITE驅(qū)動(dòng)SRM模型在開都河流域進(jìn)行研究，認(rèn)為有效活動(dòng)積溫改進(jìn)的度日數(shù)能夠提高模型中融雪速率和徑流系數(shù)的計(jì)算精度。Immerzeel等［5］在融雪模型中增加冰川模塊，以喜馬拉雅山流域?yàn)檠芯繀^(qū)，分析了流域內(nèi)的積雪變化，并對(duì)冰川融水進(jìn)行了模擬，認(rèn)為在模擬期冰川融化逐漸加速，對(duì)印度河流域的水文循環(huán)產(chǎn)生了影響。綜合眾多研究，度日因子模型能夠在高寒山區(qū)具有很好的適用性，且模型計(jì)算易于實(shí)現(xiàn)，結(jié)構(gòu)易于與其他模塊耦合，被廣泛的應(yīng)用于冰雪消融及冰雪融水徑流模擬研究中［6］。

在青藏高原地區(qū)，氣象站點(diǎn)往往設(shè)立在低海拔便利地區(qū)，空間分布十分不均，很難從地面站點(diǎn)中獲得流域氣象要素時(shí)空分布的真實(shí)情況。氣象衛(wèi)星能夠提供高時(shí)空分辨率的氣象要素反演產(chǎn)品，對(duì)氣象站點(diǎn)覆蓋空白地區(qū)能夠起到很好的補(bǔ)充作用，在過去的十幾年中得到了較為廣泛的應(yīng)用［17-20］。然而，氣象衛(wèi)星數(shù)據(jù)在影像獲取和數(shù)據(jù)處理過程中存在較多不確定性，在某些區(qū)域的適用性可能相對(duì)較差，因此在區(qū)域中使用氣象衛(wèi)星數(shù)據(jù)前，應(yīng)對(duì)衛(wèi)星產(chǎn)品的精度進(jìn)行評(píng)估，并對(duì)數(shù)據(jù)精度進(jìn)行校正，可以有效降低衛(wèi)星產(chǎn)品的系統(tǒng)誤差，提高衛(wèi)星數(shù)據(jù)產(chǎn)品的質(zhì)量。降水衛(wèi)星精度校正過程中常使用的方法包括逐步訂正［27］、最優(yōu)內(nèi)插［28］、概率密度匹配法［21］和線性模型［22］等，孫樂強(qiáng)等［23］采用逐步訂正和最優(yōu)插值兩種方法提高TMPA衛(wèi)星數(shù)據(jù)對(duì)站點(diǎn)降水量的反演精度，但是該校正方法只使用了最簡(jiǎn)單的降水徑流系數(shù)法進(jìn)行模擬驗(yàn)證，沒有與較為成熟的融雪徑流模型進(jìn)行耦合，且該方法在校正的權(quán)重函數(shù)選擇中只考慮衛(wèi)星降水網(wǎng)格與站點(diǎn)之間的距離關(guān)系，沒有考慮地形因素對(duì)權(quán)重函數(shù)的影響，使其在海拔變化劇烈的半干旱的高寒地區(qū)的應(yīng)用可能存在一定的局限性。

本文通過改進(jìn)逐步訂正法，構(gòu)建降水輸入模塊，并將其與較為廣泛應(yīng)用的融雪徑流模型（SRM）進(jìn)行耦合，使其具備校正高寒地區(qū)降水?dāng)?shù)據(jù)精度并形成半網(wǎng)格半站點(diǎn)的降水輸入特點(diǎn)，以期提高融雪模型中降水產(chǎn)流項(xiàng)的模擬效果，使模型更加準(zhǔn)確地模擬融雪徑流的變化情況。并以拉薩河流域?yàn)檠芯繀^(qū)，構(gòu)建改進(jìn)后的融雪模型，通過與原模型的精度進(jìn)行對(duì)比，驗(yàn)證改進(jìn)后模型的適用性，為半干旱高寒區(qū)水資源分配及管理提供有用的信息。

2 模型構(gòu)建

2.1 度日因子模型介紹 對(duì)高寒地區(qū)而言，選擇合適的模型對(duì)融雪過程進(jìn)行模擬是影響融雪徑流模擬效果的關(guān)鍵。度日因子模型是基于冰雪消融與正積溫之間的線性關(guān)系所建立的［6］，Snowmelt Runoff Model（SRM）是一種概念性度日因子模型，已在29個(gè)國(guó)家、100多個(gè)流域得到了廣泛應(yīng)用［13，24］。SRM模型計(jì)算日融雪徑流的原理是將當(dāng)日融雪量和降水量疊加到前日徑流退水量上。SRM模型在各個(gè)高程分區(qū)上降水?dāng)?shù)值是基于參考?xì)庀笳军c(diǎn)的數(shù)據(jù)，利用降水遞減率進(jìn)行推求，而流域地形特征對(duì)降水時(shí)空分布產(chǎn)生的差異不予考慮［25］。因此，本文在采用SRM模型中產(chǎn)流模塊的基礎(chǔ)上，耦合考慮地形因子的逐步訂正法的降水模塊。度日因子模型需要輸入的變量包括日氣溫、日降水以及日雪蓋數(shù)據(jù)等，產(chǎn)流計(jì)算方法如下［26］：

式中：Qn+1為第n+1天的流量，m3/s；Csn為融雪徑流系數(shù)；αn為第n天的度日因子，cm·℃-1·d-1； Tn為度日數(shù)，℃·d；ΔTn為由氣溫直減率得出的度日修正值，℃·d；Sn為積雪覆蓋率；Crn為降雨徑流系數(shù)；Pn為第n天的降雨量，cm；A為流域面積，km2；kn+1為第n+1天的徑流衰退系數(shù)。

按照融雪模型中徑流產(chǎn)生的來源，式（1）可以被簡(jiǎn)化為如下形式：

式中：Sn為每天新產(chǎn)生融雪產(chǎn)生的流量，m3/s；Rn為每天新產(chǎn)生的降水產(chǎn)生的流量，m3/s；QSn為第n天實(shí)測(cè)徑流中由融雪產(chǎn)生的部分，m3/s，QRn為第n天實(shí)測(cè)徑流中由降水產(chǎn)生的部分，m3/s。

當(dāng)流域海拔高程超過500 m時(shí)，需要對(duì)流域進(jìn)行分區(qū)處理，此時(shí)降水產(chǎn)流項(xiàng)公式為：

式中：R為流域中降雨的產(chǎn)流量，m3/s；Ri為第i個(gè)分區(qū)降雨的產(chǎn)流量，m3/s；Crni為第i個(gè)分區(qū)的退水系數(shù)；Pi為第i個(gè)分區(qū)的降雨量，cm；Ai為第i個(gè)分區(qū)的面積，km2；N表示分區(qū)總數(shù)。

2.2 降水輸入模塊原理 在青藏高原的高寒地區(qū)，山區(qū)的海拔變化和地形起伏變化幅度均比較大，氣象站點(diǎn)一般建設(shè)在低海拔的平原地區(qū)，而在高海拔地區(qū)無氣象站點(diǎn)覆蓋。各高程帶上的降水量大小通常使用降水遞減率進(jìn)行推算［25］，這樣會(huì)忽略區(qū)域中地形及水汽分布對(duì)降水的影響。氣象衛(wèi)星數(shù)據(jù)是無資料或缺資料地區(qū)氣象數(shù)據(jù)的重要補(bǔ)充，可以有效克服氣象站點(diǎn)時(shí)空分布極度不均問題，但是在半干旱高寒區(qū)域使用時(shí)應(yīng)該對(duì)其精度進(jìn)行校正，使之更加符合區(qū)域降水的實(shí)際情況。

逐步訂正法在降雨衛(wèi)星數(shù)據(jù)校正過程中使用權(quán)重函數(shù)來確定降水場(chǎng)各個(gè)格點(diǎn)所分配的殘差系數(shù)，但當(dāng)前研究中，大多僅考慮網(wǎng)格各個(gè)位置處的數(shù)據(jù)與地面站點(diǎn)數(shù)據(jù)的空間關(guān)系，而沒有考慮地形因素對(duì)校正過程所帶來的影響，因而在海拔變化幅度巨大的青藏高原地區(qū)的適用性可能相對(duì)較差。本文在降水輸入模塊中考慮地形因子影響，改進(jìn)了逐步訂正法［23，27-30］對(duì)降水衛(wèi)星數(shù)據(jù)的精度進(jìn)行校正，在選用距離權(quán)重時(shí)考慮地形因素對(duì)校正結(jié)果的影響，并與廣泛應(yīng)用的融雪模型進(jìn)行耦合。本文對(duì)降水輸入模塊進(jìn)行改進(jìn)的過程為：首先，以原始降水衛(wèi)星網(wǎng)格數(shù)據(jù)作為模塊的初始降水場(chǎng)，將初始場(chǎng)中實(shí)測(cè)站點(diǎn)位置處的衛(wèi)星數(shù)據(jù)降水量與實(shí)測(cè)降水量之間的差值通過權(quán)重函數(shù)插值到初始場(chǎng)網(wǎng)格的各個(gè)格點(diǎn)位置處，從而得到第一猜測(cè)場(chǎng)；然后，將第一猜測(cè)場(chǎng)作為初始場(chǎng)，重復(fù)上述步驟，最終得到較符合地面真實(shí)降水情況的降水網(wǎng)格數(shù)據(jù)；最后，將網(wǎng)格數(shù)據(jù)與地面實(shí)測(cè)數(shù)據(jù)按照高程分帶的不同進(jìn)行組合（當(dāng)高程分帶的海拔高度大于氣象站點(diǎn)平均海拔高度時(shí)，高程帶使用校正后的降水衛(wèi)星數(shù)據(jù)，反之則使用地面站點(diǎn)的降水?dāng)?shù)據(jù)），輸入到融雪模型中，模擬研究區(qū)的融雪徑流情況。具體計(jì)算步驟如下：

（1）待校正的初始場(chǎng)確定。待校正的初始場(chǎng)降水?dāng)?shù)據(jù)是從降水衛(wèi)星數(shù)據(jù)中提取，本文首先確定研究區(qū)域的矩形范圍［L，W］，此范圍內(nèi)的降水衛(wèi)星網(wǎng)格數(shù)據(jù)即為待校正的初始場(chǎng)Γ0。在對(duì)降水場(chǎng)進(jìn)行校正時(shí)，地面氣象站點(diǎn)的位置與降水網(wǎng)格節(jié)點(diǎn)位置可能不相匹配，本文使用雙線性插值法推求站點(diǎn)位置處的降水衛(wèi)星數(shù)據(jù)。

（2）降水輸入模塊中降水場(chǎng)的校正。降水輸入模塊得到降水初始場(chǎng)Γ0，其中第k個(gè)地面站點(diǎn)的實(shí)測(cè)值Mk與在初始場(chǎng)中的該點(diǎn)位置處的數(shù)值Γ0k的差值，通過權(quán)重函數(shù)對(duì)初始場(chǎng)進(jìn)行訂正得到第一猜測(cè)場(chǎng)Γ1，并將第一猜測(cè)場(chǎng)作為初始場(chǎng)重復(fù)上述步驟，直到訂正后的場(chǎng)的數(shù)值逼近觀測(cè)值：

式中：Γab1為在第一猜測(cè)場(chǎng)的格點(diǎn)坐標(biāo)（a，b）位置處的第一猜測(cè)值，mm；Γab0為在初始場(chǎng)的格點(diǎn)坐標(biāo)（a，b）位置處上的初始值，mm；Mk為站點(diǎn)k上的實(shí)測(cè)值，mm；Γab0k為初始場(chǎng)中在站點(diǎn)k位置處的衛(wèi)星數(shù)據(jù)初始值，mm；Wabk為權(quán)重因子，它的數(shù)值范圍為0～1；K為搜索半徑中的站點(diǎn)數(shù)（圖1）。

圖1 降水輸入模塊中逐步訂正法

（3）降水輸入模塊中權(quán)重函數(shù)的確定。權(quán)重函數(shù)Wabk有圓形、橢圓形和曲率橢圓形等不同的形式。本文選用比較常用的圓形權(quán)重函數(shù)［23］。

式中，R為搜索長(zhǎng)度，降水輸入模塊的搜索長(zhǎng)度的初始值為0.5°，在模型中會(huì)以0.5°為迭代步長(zhǎng)，增加至區(qū)域邊界的范圍最大值，降水輸入模塊最后選擇校正效果最好的步長(zhǎng)作為模塊最終使用的搜索長(zhǎng)度；dabk為格點(diǎn)（a，b）到地面站點(diǎn)k的距離，m。

（4）地形因子修正。青藏高原地區(qū)地形復(fù)雜，海拔的高度變化幅度大，因此在訂正過程中應(yīng)該考慮地形因素對(duì)訂正效果的影響。降水輸入模塊中地形因子為：

式中：Vabk為地形因子，mm；Dab為格點(diǎn)（a，b）位置處的海拔高度數(shù)值，m；hk為站點(diǎn)k位置處的海拔高度，m；α為地形修正系數(shù)，mm/m，其數(shù)值與地形、植被覆蓋、水汽梯度等因素有關(guān)，本文根據(jù)參考文獻(xiàn)［31］的研究結(jié)果選擇α的值為0.005 mm/m。

此時(shí)考慮地形因素影響的逐步訂正公式（4）可以修改為：

2.3 降水輸入模塊與融雪徑流模型耦合 高寒區(qū)徑流的來源一般可以分為兩部分，一部分是累積在地表面的積雪在溫度變暖的時(shí)候由消融產(chǎn)生的融水，另一部分是降水經(jīng)過下滲蒸發(fā)攔截等過程形成的降水徑流。本文積雪信息是通過MODIS衛(wèi)星數(shù)據(jù)獲得，降水信息通過TRMM衛(wèi)星數(shù)據(jù)和地面降水?dāng)?shù)據(jù)獲得。采用改進(jìn)的逐步訂正法對(duì)降水衛(wèi)星數(shù)據(jù)進(jìn)行校正后，降水衛(wèi)星數(shù)據(jù)的反演精度在一定程度得到提高，但與地面站點(diǎn)的“真實(shí)”值相比仍具有一定的差距。因此在降水輸入模塊中會(huì)對(duì)校正后的降水衛(wèi)星數(shù)據(jù)與地面站點(diǎn)數(shù)據(jù)進(jìn)行組合，被地面站點(diǎn)所覆蓋到的高程帶分區(qū)應(yīng)使用地面站點(diǎn)數(shù)據(jù)，無法被地面站點(diǎn)覆蓋到的分區(qū)使用校正后的衛(wèi)星數(shù)據(jù)。模型耦合具體流程如圖2所示。

圖2 融雪徑流模型與降水輸入模塊耦合

圖3 累積面積與高程曲線

首先，融雪模型將流域每500 m劃分一個(gè)高程帶單元，通過累積面積與高程曲線圖（圖3）得到各個(gè)分帶的平均高程（單位m），同時(shí)計(jì)算流域內(nèi)所有氣象站點(diǎn)的平均高程（單位m）。具體而言，流域內(nèi)逐點(diǎn)的高程值集合其中N表示高程點(diǎn)個(gè)數(shù)，流域中最小高程值為hmim，最大高程值為hmax，按照每500 m一個(gè)高程分區(qū)對(duì)流域進(jìn)行分帶處理，得到N個(gè)高程分區(qū)，其中（INTup表示向上取整）。對(duì)于第i個(gè)高程分區(qū)Zonei，其高程的最低點(diǎn)為hi_min，高程的最高點(diǎn)為hi_max。則此時(shí)的Elvi應(yīng)該滿足：式中 f()x表示累積面積高程曲線的擬合方程，Ai表示第i個(gè)分區(qū)的面積變量，Ai1和Ai2分別表示第i個(gè)高程分區(qū)最低點(diǎn)高程所對(duì)應(yīng)的累積面積和最高點(diǎn)高程對(duì)應(yīng)的累積面積。Elvi表示第i個(gè)高程分區(qū)的平均高程。

融雪模型計(jì)算融雪徑流的原理是將每日的融雪量與降水量疊加到前日的徑流退水量中。每個(gè)分區(qū)的融雪量由度日因子進(jìn)行計(jì)算，每個(gè)分區(qū)的降水量由地面站點(diǎn)數(shù)據(jù)和降水衛(wèi)星數(shù)據(jù)進(jìn)行計(jì)算。對(duì)于第i個(gè)分區(qū)Zonei的降水輸入項(xiàng)Pi，其數(shù)值為：

降水模塊判斷每個(gè)分區(qū)的平均高程Elvi與氣象站點(diǎn)平均高程的大小關(guān)系。分別選用地面站點(diǎn)數(shù)據(jù)或降水衛(wèi)星數(shù)據(jù)（圖4）。若≤Elvi，則該分區(qū)的降水?dāng)?shù)據(jù)使用校正后的降水衛(wèi)星數(shù)據(jù)PTrmm，若＞Elvi，則該分區(qū)的降水?dāng)?shù)據(jù)使用地面站點(diǎn)數(shù)據(jù)PSta。

3 應(yīng)用案例

圖4 不同高程帶對(duì)應(yīng)不同的降雨數(shù)據(jù)源

3.1 研究區(qū)概況 雅魯藏布江位于青藏高原的西南部，海拔范圍大，氣象站點(diǎn)稀少。本文選取雅魯藏布江流域中的拉薩河流域作為研究區(qū)（圖5），拉薩河流域面積為3.18×105km2，拉薩河流域的氣候類型是高原溫帶半干旱氣候類型，多年平均降水量為445 mm，時(shí)空分布十分不均勻，表現(xiàn)為從東南到西北地區(qū)逐漸遞減［32-33］。徑流的補(bǔ)給主要來源于雨水、融雪融水和地下水，其中融水補(bǔ)給占到徑流總量25%［34］，是雅魯藏布江徑流的重要組成部分。

圖5 研究區(qū)位置以及氣象站點(diǎn)和水文站點(diǎn)分布

3.2 數(shù)據(jù)來源和處理 本文使用的數(shù)據(jù)主要包括數(shù)字高程（DEM）、積雪覆蓋、降水?dāng)?shù)據(jù)等。

3.2.1 數(shù)字高程模型（DEM）DEM來源于NASA的“EARTHDATA”平臺(tái)，數(shù)據(jù)版本為ASTER Global Digital Elevation Model V002，空間分辨率為30 m，重采樣為500 m。拉薩河流域的高程范圍為3581～7112m，將DEM按照500 m進(jìn)行高程帶劃分，拉薩河流域可以分為7個(gè)高程區(qū)（如表1）。

表1 拉薩河流域高程分區(qū)

3.2.2 MODIS積雪數(shù)據(jù) 融雪是高寒地區(qū)徑流的重要來源，但是在人跡罕至的高寒地區(qū)，傳統(tǒng)的地面觀測(cè)很難獲得區(qū)域內(nèi)積雪的時(shí)空分布。MODIS（Moderate-Resolution Imaging Spectroradiometer）產(chǎn)品能夠提供多時(shí)段不同分辨率的積雪覆蓋產(chǎn)品，在積雪研究中發(fā)揮著十分重要的作用。本文使用MODIS系列數(shù)據(jù)中MOD10A1積雪二元值數(shù)據(jù)，其空間分辨率為500 m，時(shí)間分辨率是1 d。從MODIS中獲取積雪數(shù)據(jù)的算法是SNOMAP算法［35］，利用積雪在MODIS的可見光波段具有較強(qiáng)的反射率，而在近紅外波段具有較強(qiáng)的吸收能力，通過計(jì)算歸一化差值積雪指數(shù)NDSI來識(shí)別積雪：

式中：band4為第四波段（綠光）（0.545～0.565 μm）；band6為第六波段（近紅外波段）（1.652～1.680 μm）。

在非林積雪區(qū)域，當(dāng)像元的NDSI≥0.4時(shí)，該像元會(huì)被定義為雪。由于水體的反射特性和積雪相似，因此利用NDSI識(shí)別出的積雪中會(huì)存在水體，本文利用MODIS的第2波段band2（0.841～0.876 μm），設(shè)置波段閾值去除水體的影響，即設(shè)置band2≥0.11的像元為無水體的雪［36-37］。通過時(shí)間序列和空間相關(guān)性對(duì)積雪識(shí)別結(jié)果進(jìn)行去云處理［38］，最后得到拉薩河流域多年逐日積雪覆蓋率變化曲線（3—12月）（圖6）。

為了驗(yàn)證NDSI和去云處理法得到的積雪數(shù)據(jù)的可靠性，本文將計(jì)算得到的積雪覆蓋率與“寒區(qū)旱區(qū)科學(xué)數(shù)據(jù)中心（http：//westdc.westgis.ac.cn/）”的“青藏高原地區(qū)MODIS逐日無云積雪產(chǎn)品（2002—2010）”［39］產(chǎn)品在日尺度上進(jìn)行對(duì)比，計(jì)算各個(gè)分區(qū)和流域整體上的相關(guān)系數(shù)、均方根誤差和標(biāo)準(zhǔn)差，結(jié)果如表2所示。在流域整體上，本文計(jì)算得到的積雪數(shù)據(jù)與寒旱所的積雪產(chǎn)品之間的相關(guān)系數(shù)為0.816，均方根誤差為0.076，標(biāo)準(zhǔn)差為0.131，因此兩套數(shù)據(jù)在流域整體上的相關(guān)性強(qiáng)，偏差小，精度相當(dāng)。從各個(gè)分區(qū)來看，相關(guān)系數(shù)數(shù)值最小的分區(qū)為G分區(qū)，相關(guān)系數(shù)僅有0.451，但是G分區(qū)的面積只占整個(gè)流域總面積的0.02%，對(duì)積雪產(chǎn)品整體精度影響較小。C、D、E分區(qū)的相關(guān)系數(shù)都在0.76以上，相關(guān)性較強(qiáng)，且均方根誤差和標(biāo)準(zhǔn)差數(shù)值較小，偏差較小，它們的面積占流域總面積的81.10%。因此本文通過NDSI和去云處理所得到的積雪數(shù)據(jù)與寒區(qū)旱區(qū)科學(xué)數(shù)據(jù)中心所發(fā)布的逐日無積雪產(chǎn)品在精度上相當(dāng)，數(shù)據(jù)的可靠性較強(qiáng)，可以進(jìn)一步進(jìn)行融雪模擬研究。

表2 積雪數(shù)據(jù)的可靠性驗(yàn)證

從積雪覆蓋率變化曲線（圖6）中可以看出，拉薩河流域不同高程分區(qū)的積雪覆蓋率多年逐日變化趨勢(shì)是相同的，從積雪覆蓋率的逐月變化情況上看，積雪覆蓋率總是在4月份達(dá)到最大值然后逐月消退，在8月份多年平均積雪覆蓋率達(dá)到了全年的最低值。因此本文選擇4—11月份作為拉薩河流域融雪季進(jìn)行融雪模擬。從積雪覆蓋率的分帶情況上看，從A分帶到G分帶的積雪覆蓋率逐漸增加，積雪覆蓋率與各個(gè)分區(qū)的高程值呈現(xiàn)出正相關(guān)的關(guān)系，即高程分區(qū)的高程越高，分帶的積雪覆蓋率數(shù)值越大。

3.2.3 降水?dāng)?shù)據(jù) 降水?dāng)?shù)據(jù)是SRM模型的重要輸入變量，降水在各個(gè)分區(qū)帶的精度影響了模型模擬融雪徑流的效果。拉薩河流域氣象站點(diǎn)的平均海拔高度是3942.80 m，而流域的海拔高度范圍是3581～7112 m，在傳統(tǒng)的SRM模型中只有A分區(qū)具有實(shí)測(cè)的降水?dāng)?shù)據(jù)，在其他分區(qū)中的降水?dāng)?shù)據(jù)一般使用降水遞減率推求。降水?dāng)?shù)值隨海拔高度線性變換的處理，忽略了流域中地形、海拔高度以及水汽等因素的影響，使得降水?dāng)?shù)據(jù)出現(xiàn)均坦化的問題［25］。融雪模型與降水輸入模塊耦合的模型，使用降水衛(wèi)星數(shù)據(jù)對(duì)氣象站點(diǎn)覆蓋不到的高程帶進(jìn)行補(bǔ)充。TRMM（Tropical Rainfall Measuring Mission）3B42數(shù)據(jù)產(chǎn)品的空間分辨率是 0.25°×0.25°，時(shí)間分辨率是日，空間范圍是50°S～ 50°N，180°W～180°E，TRMM數(shù)據(jù)相對(duì)較高的時(shí)間和空間分辨率使其在眾多區(qū)域上得到了廣泛的應(yīng)用，降水輸入模塊調(diào)用的降水?dāng)?shù)據(jù)包括地面降水站點(diǎn)數(shù)據(jù)（圖5）和TRMM數(shù)據(jù)，降水?dāng)?shù)據(jù)信息見表3。

圖6 拉薩河流域不同高程分區(qū)多年平均逐日積雪覆蓋率（3—12月）

表3 降水?dāng)?shù)據(jù)的時(shí)空分辨率和使用的時(shí)間段

為了保證降水輸入模塊對(duì)降水衛(wèi)星數(shù)據(jù)校正的可靠性，本文使用的氣象站點(diǎn)數(shù)據(jù)向拉薩河流域外進(jìn)行了擴(kuò)展，氣象站點(diǎn)使用流域內(nèi)及流域外的11個(gè)站點(diǎn)，站點(diǎn)分布如圖5所示。

3.3 模型率定和驗(yàn)證 在模擬中，利用2001—2007年的日流量數(shù)據(jù)進(jìn)行率定，利用2008—2014年的日流量數(shù)據(jù)進(jìn)行驗(yàn)證，并利用Nash-Sutcliffe系數(shù)NSE和百分偏差PBIAS指標(biāo)對(duì)改進(jìn)的融雪模型精度進(jìn)行評(píng)價(jià)［40］：

式中：Qi為實(shí)測(cè)的日流量，m3/s；Q′i為模擬的日流量，m3/s；為實(shí)測(cè)的流量平均值，m3/s；n為模擬的天數(shù)。

NSE范圍為-∞～1，其數(shù)值愈接近1，說明模型模擬的精度愈高。當(dāng)NSE的數(shù)值介于0～1之間時(shí)，說明模型具有一定的模擬能力，而當(dāng)NSE的數(shù)值小于0時(shí)，說明模型對(duì)融雪徑流不具備模擬能力。

百分偏差的公式如下：

式中，PBIAS的數(shù)值愈接近0說明模型的精度愈高。當(dāng)PBIAS的數(shù)值為正時(shí)，說明模型低估了融雪流量，當(dāng)PBIAS的數(shù)值為負(fù)時(shí)，說明模型高估了融雪流量。

4 結(jié)果與分析

本文基于積雪數(shù)據(jù)以及站點(diǎn)和TRMM降水?dāng)?shù)據(jù)，對(duì)2001—2014年的拉薩河流域的融雪徑流過程進(jìn)行模擬，其中2001—2007年為率定期，2008—2014年為驗(yàn)證期，得到耦合降水輸入模塊的融雪模型，以及由此估算得到的融雪流量過程和精度評(píng)價(jià)結(jié)果，并與改進(jìn)前的融雪模型得出的結(jié)果進(jìn)行對(duì)比。

4.1 降水衛(wèi)星數(shù)據(jù)校正結(jié)果 氣象衛(wèi)星數(shù)據(jù)具有觀測(cè)范圍廣，觀測(cè)范圍連續(xù)等特點(diǎn)，能夠?yàn)闅庀笳军c(diǎn)稀缺的高寒地區(qū)有效地提供數(shù)據(jù)上的補(bǔ)充，但是在區(qū)域應(yīng)用前應(yīng)先對(duì)其進(jìn)行評(píng)估和校正。降水輸入模塊對(duì)降水衛(wèi)星的校正結(jié)果如表4所示，從表4可以看出，原始的TRMM降水?dāng)?shù)據(jù)的精度較低，地面站點(diǎn)與站點(diǎn)位置處的TRMM數(shù)據(jù)的相關(guān)系數(shù)（CC）為0.307，百分偏差為49.0%，說明在絕大部分時(shí)間內(nèi)TRMM日降水?dāng)?shù)據(jù)反演拉薩河流域的降水效果相對(duì)較差。在降水輸入模塊中，一共進(jìn)行了5次迭代，即校正步數(shù)為5，每經(jīng)過一次迭代降水輸入模塊會(huì)計(jì)算地面站點(diǎn)數(shù)據(jù)與在站點(diǎn)位置處的衛(wèi)星數(shù)據(jù)的相關(guān)系數(shù)和百分偏差，并將生成的猜測(cè)場(chǎng)作為初始場(chǎng)進(jìn)行重新校正。經(jīng)過降水輸入模塊修正的逐步訂正法校正后，TRMM降水衛(wèi)星的精度改善較為明顯，相關(guān)系數(shù)從0.307上升到0.613～0.839，而百分偏差從49.0%下降到0.8%～1.5%。尤其是當(dāng)校正步數(shù)等于4的時(shí)候降水衛(wèi)星對(duì)地面降水的反演精度達(dá)到最大，此時(shí)TRMM降水衛(wèi)星數(shù)據(jù)與地面站點(diǎn)數(shù)據(jù)之間的相關(guān)系數(shù)為0.839，百分偏差為0.7%。可以看出改進(jìn)后逐步訂正法對(duì)衛(wèi)星的精度提升較大，經(jīng)過校正后的衛(wèi)星降水?dāng)?shù)據(jù)通過降水輸入模塊與地面站點(diǎn)進(jìn)行組合驅(qū)動(dòng)融雪徑流進(jìn)行模擬。

表4 拉薩河流域TRMM降水衛(wèi)星校正結(jié)果

圖7 拉薩河流域融雪徑流模擬值和實(shí)測(cè)值

4.2 模型模擬效果評(píng)價(jià) 本文基于改進(jìn)的逐步訂正法構(gòu)建了降水輸入模塊，并將其與融雪模型進(jìn)行耦合，利用高程分帶組合輸入半網(wǎng)格半站點(diǎn)的降水輸入數(shù)據(jù)，以拉薩河為例對(duì)半干旱高寒區(qū)的融雪徑流進(jìn)行了模擬，從圖7和表5中可以看出，傳統(tǒng)的融雪徑流模型和采用耦合降水輸入模塊的融雪模型模擬值與實(shí)測(cè)值有很好的擬合效果，能夠大體上模擬出融雪徑流的實(shí)際變化趨勢(shì)。在率定期2001—2007年，原模型的納什效率系數(shù)NSE和PBIAS分別為0.695和30.3%，改進(jìn)后的融雪模型的NSE和PBIAS分別為0.741和14.5%。在驗(yàn)證期2008—2014年，原模型的納什效率系數(shù)NSE和PBIAS分別為0.744和22.8%，改進(jìn)后的融雪模型的NSE和PBIAS分別為0.770和15.4%。改進(jìn)后的融雪模型提高了NSE的數(shù)值和降低了PBIAS數(shù)值，減小了相對(duì)誤差，提高了模擬精度，說明融雪模型降水項(xiàng)改進(jìn)后的高寒地區(qū)的融雪徑流過程模擬要優(yōu)于改進(jìn)前的模擬效果，更加接近流域的實(shí)際產(chǎn)流情況。

在拉薩河流域應(yīng)用SRM模型和耦合降水輸入模塊的模型的模擬值和實(shí)測(cè)值擬合程度較好，它們的變化趨勢(shì)一致。但無論是原融雪模型還是改進(jìn)后的模型，對(duì)退水的效果模擬都很差，總是在模擬年開始時(shí)降到“無水區(qū)”，這主要是因?yàn)樵诎敫珊蹈吆貐^(qū)，冬季和春季的降水相對(duì)較少，由降水產(chǎn)生的降水徑流相對(duì)較少，而此時(shí)的溫度較低融雪量也相對(duì)較少，因此融雪徑流總量相對(duì)于其他降水充沛和溫度偏暖時(shí)期徑流總量要少的多，更重要的是SRM融雪徑流模型中不包含地下水模擬的模塊，無法對(duì)研究區(qū)的基流進(jìn)行模擬。

5 討論

度日因子模型是半干旱高寒地區(qū)融雪徑流模擬的重要工具，模型結(jié)構(gòu)簡(jiǎn)單，所需的數(shù)據(jù)較少且易得，因此在融雪過程模擬中得到了廣泛的應(yīng)用。隨著遙感產(chǎn)品種類越來越多，精度越來越高，融雪過程的模擬效果也不斷增強(qiáng)。但是在融雪徑流模擬中仍然存在一系列的不確定性。模型不確定性包括水文模型參數(shù)的不確定性和水文模型輸入變量的不確定性等內(nèi)容。對(duì)于模型參數(shù)而言，在度日模型中會(huì)將風(fēng)吹雪和雪的積累與消融等動(dòng)態(tài)過程進(jìn)行簡(jiǎn)化，只使用簡(jiǎn)單的經(jīng)驗(yàn)公式或者經(jīng)驗(yàn)參數(shù)代替。對(duì)于模型的輸入變量，降水和積雪是模型的重要輸入變量，青藏高原地區(qū)海拔較高，下墊面變化復(fù)雜，植被、坡度和坡向等都會(huì)對(duì)模型的輸入變量帶來一定的誤差。如何減小輸入變量的誤差，提高融雪徑流過程模擬的精度是未來研究的重要方向［43］。

5.1 降水類型對(duì)模擬結(jié)果影響 降水類型（雨、雪和雨夾雪）會(huì)影響融雪徑流的模擬精度，不同的降水類型形成徑流的過程不同，降雪在冬季累積到地面上，在溫度升高時(shí)，通過雪的累積和消融過程形成徑流。降水直接降到河流中，或經(jīng)過土壤下滲，蒸散等過程后匯聚到河流中［44］。因此不同的降水類型對(duì)融雪徑流模擬的影響也是不同的，可以通過分離降水類型來提高青藏高原地區(qū)融雪徑流過程模擬的精度。在水文模型中，常利用單溫度閾值法對(duì)降水進(jìn)行分割，單溫度閾值法常設(shè)定0℃為臨界溫度閾值，當(dāng)溫度高于臨界溫度時(shí)，降水的形態(tài)為雨；當(dāng)溫度低于臨界溫度時(shí)，降水的形態(tài)為雪［45］。單溫度閾值法會(huì)導(dǎo)致融雪模擬過程存在一定的不確定性，首先，不同降水形態(tài)會(huì)在-2℃～4℃之間進(jìn)行轉(zhuǎn)換，即在此溫度區(qū)間內(nèi)，雨、雪和雨夾雪均有可能發(fā)生，單溫度閾值法設(shè)定固定的閾值會(huì)導(dǎo)致水文模型高估或低估降水量和降雪量［46］。其次，在流域復(fù)雜的地形條件下，使用單溫度閾值法忽略了下墊面的異質(zhì)性。雙溫度閾值法，設(shè)定臨界最低溫度和臨界最高溫度，當(dāng)溫度高于臨界最高溫度時(shí)，降水的形態(tài)為雨；當(dāng)溫度低于臨界最低溫度時(shí)，降水的形態(tài)與雪；當(dāng)溫度介于臨界最低溫度和臨界最高溫度之間時(shí)，降水的形態(tài)為雨夾雪［47］。Ding等［48］在青藏高原地區(qū)通過研究不同的降水類型與溫度之間的關(guān)系，使用一種新的指數(shù)曲線法對(duì)降水類型進(jìn)行分割，他認(rèn)為雙溫度閾值法對(duì)青藏高原地區(qū)降水類型的分割十分準(zhǔn)確，但是該方法并沒有在融雪徑流模型中進(jìn)行驗(yàn)證。

本文所使用的地面站點(diǎn)數(shù)據(jù)中只有1979年以前的數(shù)據(jù)標(biāo)記了降水類型。本文基于拉薩河流域內(nèi)外4個(gè)站點(diǎn)（拉薩站、當(dāng)雄站、那曲站和嘉黎站），計(jì)算了1957—1979年拉薩河流域不同降水類型的比例。從圖8中可以看出，不同的降水形態(tài)轉(zhuǎn)化主要發(fā)生在1℃～9℃之間。當(dāng)溫度小于1℃時(shí)，降水的主要形態(tài)為降雪；隨著溫度的升高，降雪發(fā)生的頻率逐漸減小，但降雪發(fā)生的頻率要大于降雨的頻率；當(dāng)溫度升高到2.93℃時(shí)，降雪的頻率等于雨夾雪的頻率，此時(shí)降水與降雪開始轉(zhuǎn)換；當(dāng)溫度達(dá)到4.52℃時(shí)，降雨的頻率等于雨夾雪的頻率，此時(shí)雨夾雪的頻率也在峰值附近，此后隨著溫度的升高降雪的頻率逐漸降低，直至頻率為0。

圖8 不同降水類型在降水中頻率隨溫度變化

圖9 使用單溫度閾值法和雙溫度閾值法在拉薩河流域融雪徑流模擬（2011年）

根據(jù)Liu和Ding等［47-48］的指數(shù)曲線法，本文計(jì)算了拉薩河流域降水類型分割的雙溫度閾值函數(shù)（式（12））。在拉薩河流域臨界最低溫度和臨界最高溫度分別為2.93℃和4.52℃，也即當(dāng)溫度小于2.93℃時(shí)，降水中100%的降水量均為降雪量；當(dāng)溫度高于4.52℃時(shí)，降水中100%的降水量均為降雨量；當(dāng)溫度介于兩者之間時(shí)，降水事件中的降雪量所占的比例根據(jù)式（12）計(jì)算，而降雨量等于所占的比例等于1-降雪量比例。

式中：Sp為降水量中降雪量所占的比例；T為溫度，℃。

以2011年站點(diǎn)為例，利用雙溫度閾值法對(duì)降水進(jìn)行分離，分別計(jì)算降水量中降雨和降雪所占的比例，并將分離后的降水變量驅(qū)動(dòng)融雪徑流模型進(jìn)行模擬。使用單閾值法分離降水類型的NSE數(shù)值為0.896，使用雙溫度閾值法分離降水類型的NSE數(shù)值為0.908。使用兩種方法模擬結(jié)果的不同主要體現(xiàn)在9月中旬到12月初，特別是在10月份左右（圖9），雙溫度閾值法在一定程度上降低了模擬的高估值。

5.2 季節(jié)尺度對(duì)模擬結(jié)果影響 時(shí)間尺度同樣影響融雪徑流的模擬精度，一般而言，模型在融雪季的模擬精度要高于非融雪季的精度。對(duì)融雪徑流模型在春、夏、秋、冬尺度上分別進(jìn)行模擬，使用單閾值法的NSE數(shù)值分別為-2.321、0.625、0.748和-0.817。從模擬結(jié)果中可以看出，融雪徑流模型在夏季和秋季的模擬精度遠(yuǎn)高于春季和冬季的模擬精度，也就是融雪季的精度要遠(yuǎn)高于非融雪季的精度。利用雙溫度閾值法在四季尺度上模擬，則夏季和秋季的NSE數(shù)值分別為0.621與0.750，春季和冬季的NSE數(shù)值分別為-1.037和-0.782。從結(jié)果中可以看出，模型在9月中旬到2月初的精度得到一定程度的提高，雙溫度閾值法的融雪徑流模擬主要提高了非融雪期的模擬精度。因此使用雙溫度閾值法對(duì)降水類型進(jìn)行分割，可以在一定程度上提高融雪徑流的模擬精度，這也是融雪徑流模擬接下來進(jìn)行改進(jìn)的一個(gè)方向。

6 結(jié)論

本文基于改進(jìn)的逐步訂正法構(gòu)建了降水輸入模塊，并與融雪徑流模型SRM相耦合，以拉薩河作為研究區(qū)評(píng)估改進(jìn)后的耦合模型在半干旱高寒地區(qū)對(duì)融雪徑流的模擬效果。相比于傳統(tǒng)的SRM模型，耦合降水輸入模塊的融雪模型能夠有效提高融雪徑流模擬效果。本文得出以下主要結(jié)論：（1）原始的TRMM降水衛(wèi)星反演地面站點(diǎn)降水的相關(guān)系數(shù)、百分偏差分別為0.307和49.0%，精度相對(duì)較低。經(jīng)過降水輸入模塊校正的降水衛(wèi)星反演地面降水的精度得到顯著的提高，其相關(guān)系數(shù)、百分偏差分別達(dá)到0.839和0.7%。（2）改進(jìn)后的融雪模型的NSE和PBIAS在率定期分別為0.741和14.5%，比原始的SRM模型分別高出了0.046和減小了15.8%。NSE和PBIAS在驗(yàn)證期分別為0.770和15.4%，比原始的SRM模型分別高出了0.026和減小了7.4%。經(jīng)過改進(jìn)的融雪模型能夠提高半干旱高寒區(qū)融雪徑流模擬的精度。

本文利用基于地形因子改進(jìn)的逐步訂正法構(gòu)建降水輸入模塊并將其與融雪模型進(jìn)行耦合，在半干旱高寒區(qū)的拉薩河流域進(jìn)行了試驗(yàn)，驗(yàn)證了模型的合理性和有效性，為缺資料地區(qū)的半干旱高寒區(qū)融雪徑流模擬提供了一種有效的研究思路。