基于MODIS反演雪深的融雪徑流模擬

陳智梁，王 娟，李春紅，王冉旋，王 奕，馬志貴

（1.國(guó)家能源集團(tuán)新疆吉林臺(tái)水電開(kāi)發(fā)有限公司，新疆伊犁835000；2.南瑞集團(tuán)有限公司/國(guó)網(wǎng)電力科學(xué)研究院，南京211106）

新疆地區(qū)河流大多以冰雪融水為主要補(bǔ)給水源，其春、秋季來(lái)水以冰雪融水為主，夏季洪水多為融雪與暴雨的疊加。對(duì)于此類地區(qū)而言，融雪徑流是河流的基礎(chǔ)性來(lái)水，準(zhǔn)確的融雪徑流預(yù)報(bào)可以提高水資源利用效率，有效預(yù)防洪澇災(zāi)害［1-3］。

SRM 融雪徑流模型因結(jié)構(gòu)簡(jiǎn)單、應(yīng)用效果較好已廣泛應(yīng)用于融雪徑流預(yù)報(bào)［4-18］，它以積雪覆蓋面積作為模型驅(qū)動(dòng)之一，計(jì)算逐日融雪徑流過(guò)程。隨著遙感影像時(shí)空分辨率和光譜分辨率的逐漸提高，遙感已成為有效的積雪監(jiān)測(cè)手段［19］，其中MODIS 遙感數(shù)據(jù)因具有較高的時(shí)空分辨率，其積雪覆蓋產(chǎn)品通常被作為SRM 模型的輸入［4-15］。NSIDC（National Snow and Ice Data Center）發(fā)布的基于MODIS 傳感器的積雪覆蓋產(chǎn)品有兩類：全球逐日積雪覆蓋產(chǎn)品MOD10A1/MYD10A1 和全球8日合成的積雪覆蓋產(chǎn)品MOD10A2/MYD10A2，空間分辨率均為500 m［16-18］。逐日積雪覆蓋產(chǎn)品的積雪分類精度受云層影響嚴(yán)重，在多云或陰天時(shí)，平均積雪識(shí)別率僅為17.8%［18，19］，因此融雪徑流預(yù)報(bào)通常使用合成產(chǎn)品MOD10A2。在春秋積雪不穩(wěn)定季節(jié)，積雪覆蓋面積在短時(shí)間內(nèi)就會(huì)發(fā)生明顯變化，譬如一場(chǎng)薄薄的新雪會(huì)帶來(lái)較大的積雪覆蓋面積，又會(huì)在短短幾天內(nèi)融化殆盡，因此8日合成的遙感積雪產(chǎn)品無(wú)法反映積雪覆蓋信息的逐日變化過(guò)程，逐日融雪徑流計(jì)算也將產(chǎn)生較大的偏差。

SRM 模型假定雪蓋區(qū)在單位計(jì)算時(shí)段內(nèi)有充足的積雪可融化，雪深不設(shè)限，僅根據(jù)流域的積雪覆蓋率、溫度和度日因子計(jì)算逐日融雪徑流過(guò)程。若每日的積雪覆蓋率準(zhǔn)確，則該假定對(duì)融雪徑流計(jì)算結(jié)果的影響不大；但如前所述，可用的遙感積雪覆蓋產(chǎn)品為8日合成，無(wú)準(zhǔn)確的逐日積雪覆蓋面積，因此要保證融雪徑流計(jì)算精度，需要獲取相對(duì)準(zhǔn)確的逐日積雪覆蓋信息。積雪覆蓋面積是積雪的二維空間體現(xiàn)，精細(xì)化空間分布的雪深既可計(jì)算雪蓋面積也反映了積雪量，同時(shí)也是融雪徑流計(jì)算的中間變量，因此將空間分布的雪深信息作為SRM 模型的輸入，進(jìn)行雪蓋面積、積雪量、融雪量間的轉(zhuǎn)換計(jì)算是一種可行的解決途徑。首先由空間分布的雪深信息計(jì)算雪蓋面積、雪蓋率，并將其代入SRM 模型計(jì)算融雪徑流；而后用融雪徑流計(jì)算所得的融雪深度計(jì)算雪深變化，進(jìn)而得到雪深變化后的雪蓋率，此雪蓋率即為SRM 模型下一時(shí)段的輸入。依次循環(huán)，即可計(jì)算逐日的融雪量、融雪深度和積雪覆蓋面積。實(shí)測(cè)融雪徑流和雪深空間分布信息是計(jì)算過(guò)程中模型參數(shù)調(diào)整校核的依據(jù)。

綜上，為了減少因積雪信息不足造成的融雪徑流預(yù)報(bào)誤差，本文采用MODIS 遙感信息反演積雪深度，用雪深信息作為SRM 模型的輸入，并在模型中增加逐日雪深和雪蓋率計(jì)算模塊，最終以實(shí)測(cè)徑流和雪深數(shù)據(jù)作為融雪徑流模擬的校核依據(jù)，得到雪深、積雪覆蓋和徑流的逐日變化過(guò)程。該方法在新疆某河A水庫(kù)以上流域（以下簡(jiǎn)稱A流域）進(jìn)行了試驗(yàn)應(yīng)用。

1 研究區(qū)概況

A 流域?yàn)楠M長(zhǎng)形，集水面積2 576 km2，高程范圍1 891～4 614 m，流域上游高海拔區(qū)分布有永久性積雪和冰川；洪水成因主要是河源區(qū)高山帶冰川及永久性積雪的融冰雪水，年內(nèi)最大洪峰主要集中在6-7月，洪峰過(guò)程一般為5～10 d。夏季降雨多集中在上游河段，洪水過(guò)程可分為消融型洪水和消融與暴雨疊加型的混合洪水；全年過(guò)程多為一個(gè)融雪型主峰，在此基礎(chǔ)上疊加若干短時(shí)段的暴雨洪水。

圖1 A流域DEMFig.1 DEM of A basin

流域內(nèi)布設(shè)有水文自動(dòng)測(cè)報(bào)系統(tǒng)，包括5 個(gè)雨量氣溫站（A1～A5）、2個(gè)雨量氣溫雪深站（S1、S2）和1個(gè)水庫(kù)站，雨量氣溫站大多布設(shè)在河谷地區(qū)，測(cè)站分布見(jiàn)圖2，圖中測(cè)站旁標(biāo)注數(shù)字表示測(cè)站位置的高程。采集的水文氣象數(shù)據(jù)為2016-2020年，數(shù)據(jù)時(shí)段長(zhǎng)為日；其中僅2019、2020年有實(shí)測(cè)雪深數(shù)據(jù)。

圖2 A流域測(cè)站分布圖Fig.2 Station distribution in A basin

2 原理方法

2.1 MODIS反演雪深

已發(fā)布的MODIS 逐日和合成積雪資料均為積雪覆蓋信息［16-19］，無(wú)雪深數(shù)據(jù)。為獲取空間分布的雪深信息，本文基于傅華等［20］研制的MODIS積雪遙感監(jiān)測(cè)系統(tǒng)，集合地表不同覆蓋物的光譜特征，結(jié)合季節(jié)、地形、下墊面等積雪深度的影響因素，利用MODIS 高光譜、多波段反射率和測(cè)站觀測(cè)雪深，采用數(shù)學(xué)統(tǒng)計(jì)方法，建立MODIS積雪深度回歸模型。從而在計(jì)算積雪面積、覆蓋率的基礎(chǔ)上，獲取空間分布的雪深、雪水當(dāng)量信息，制作500 m 分辨率的積雪深度產(chǎn)品。因受云層影響，雪深為晴空日或10日合成產(chǎn)品［20，21］。

積雪深度的反演主要依賴于積雪深度與可見(jiàn)光波段的反射率間存在較好的相關(guān)性［21-24］：當(dāng)積雪深度小于20 cm 時(shí)，雪面反射率隨積雪深度的增加而增加；當(dāng)積雪深度大于20 cm 時(shí)，雪面反射率隨積雪深度增加的趨勢(shì)減緩，當(dāng)積雪達(dá)到一定深度后，雪面反射率趨于飽和。

當(dāng)積雪深度相同時(shí)，下墊面類型不同，反射率也存在差別，因此進(jìn)行積雪深度反演時(shí)，必須要考慮下墊面的差異，根據(jù)積雪深度反演原理，用數(shù)理統(tǒng)計(jì)方法（逐步回歸判別、貝葉斯判別等），建立積雪深度與影響因子之間的線性關(guān)系［21］：

式中：S為積雪深度；x1，...，xn為影響因子，根據(jù)雪深反演原理定義為MODIS 多個(gè)特定通道反射率的函數(shù)；a1，...，an+1為回歸系數(shù)，通過(guò)歷史資料的統(tǒng)計(jì)分析確定。因影響因子較多，本文僅針對(duì)A 流域山地地形，依據(jù)土地覆蓋類型（森林、灌木、草原、荒漠等）將其進(jìn)行分類處理，每一類單獨(dú)進(jìn)行回歸計(jì)算；計(jì)算結(jié)果結(jié)合坡度、坡向等對(duì)積雪分布及反射率的影響進(jìn)行訂正，從而提高反演精度［25，26］。

2.2 基于雪深的SRM模型

因不同垂向高程帶的融雪徑流特征差異明顯，SRM 模型計(jì)算中通常將流域劃分為若干區(qū)（高程帶），每區(qū)單獨(dú)計(jì)算融雪徑流過(guò)程，再匯流至流域出口斷面進(jìn)行疊加。SRM 模型的大致原理為度日法［26］，即基于積雪覆蓋面積和度日因子逐日計(jì)算融雪徑流和降水徑流，并將它們疊加到當(dāng)日的退水流量上，得到每日的日徑流量［4-14］。本文將積雪覆蓋率替換為雪深，則基于雪深的SRM融雪徑流模型的計(jì)算公式如下：

式中：Q為日平均徑流量，m3/s；C是徑流系數(shù)，其中Cs、Cr分別為融雪和降雨的徑流系數(shù)；α為度日因子，表示單位時(shí)間、溫度的融雪徑流深，cm/℃·d；T為度日數(shù)，℃·d；△T為計(jì)算區(qū)域平均高程與氣溫站高程不同而產(chǎn)生的氣溫調(diào)整值，℃·d；P為日降水量，cm；K是退水系數(shù)，表示在沒(méi)有融雪或降水時(shí)期逐日徑流的下降率；n為徑流計(jì)算的日數(shù)序列；常數(shù)項(xiàng)是徑流深到徑流量的轉(zhuǎn)換系數(shù)［25］。遙感反演雪深產(chǎn)品的空間分辨率為500 m 網(wǎng)格，As為計(jì)算區(qū)域內(nèi)雪深>0 的所有網(wǎng)格面積，A為計(jì)算區(qū)域總面積，As與A的比值為無(wú)量綱數(shù)，代表了區(qū)域的積雪覆蓋率。

對(duì)于10日合成雪深產(chǎn)品而言，假定某產(chǎn)品為第n日的雪深，則下一個(gè)合成產(chǎn)品為第n+9日雪深，兩個(gè)產(chǎn)品的間隔期內(nèi)無(wú)雪深數(shù)據(jù)，因此需要以合成雪深為初始值，依據(jù)降水、氣溫等信息，采用公式計(jì)算得到間隔期內(nèi)逐日雪深。

雪深計(jì)算首先需判別當(dāng)日（第n日）氣溫是否超過(guò)臨界溫度：低于臨界溫度時(shí)，降水通過(guò)雪水當(dāng)量轉(zhuǎn)化為雪深，無(wú)融雪，雪深增加；高于臨界溫度時(shí)，積雪消融，雪深減小。積雪深度計(jì)算公式為：

式中：Hn為第n日雪深，cm；Hm，n為第n日積雪消融的深度，cm；Hs，n為第n日降雪深度，cm。

降雪深度和消融的積雪深度均通過(guò)相應(yīng)的水量換算得到，計(jì)算公式為：

式中：Rm，n為降水量或徑流量，m3，為SRM 模型計(jì)算量；ρw、ρs分別為水、雪密度，g/cm3。

通常，新下的松軟的雪密度為0.04～0.1 g/cm3，融雪期的雪密度為0.6～0.7 g/cm3，雪密度的影響因素較多，新疆地區(qū)雪密度時(shí)空分布可參見(jiàn)黃慰軍等研究成果［28］。因融雪徑流量、積雪深度互為相關(guān)計(jì)算值，因此需借助合成的反演積雪深和實(shí)測(cè)徑流量校核，由融雪徑流模型迭代計(jì)算得到逐日的雪深和徑流過(guò)程。

3 結(jié)果分析

3.1 流域分區(qū)

因氣溫隨高程變化，而氣溫是融雪徑流的重要影響因素，因此需要對(duì)垂直高程大于500 m的流域進(jìn)行高程帶劃分?；贏STER GDEM-Global Elevation Data 的DEM 數(shù)據(jù)（空間分辨率30 m），采用ArcGIS 軟件的Arc Hydro Tools 工具包生成數(shù)字流域，提取相應(yīng)的地形信息；在此基礎(chǔ)上采用ArcGIS 的空間分析工具條Contour 進(jìn)行高程帶劃分，每500 m 左右分為一帶。A 流域依據(jù)500 m等高線進(jìn)行整合，分為4個(gè)高程帶，見(jiàn)表1、圖3。

圖3 A流域高程分帶圖Fig.3 Elevation distribution in A basin

表1 A流域高程帶信息表Tab.1 Elevation information in A basin

3.2 MODIS反演雪深分析

MODIS 反演的積雪深度產(chǎn)品為二進(jìn)制的無(wú)符號(hào)短整型格式（unsigned short），頭文件后綴為*.HDR；數(shù)據(jù)文件采用等經(jīng)緯度投影；A 流域反演數(shù)據(jù)行列數(shù)為：columns=1 518，lines=548；數(shù)據(jù)范圍為0～255，其中0～49 代表雪深值，cm；250～252 代表水體（雪深為0），253、254代表云層（無(wú)效信息）。將雪深產(chǎn)品與研究區(qū)進(jìn)行空間疊置，可得流域內(nèi)500 m分辨率的雪深分布信息，圖4為2019-01-08 A流域雪深空間分布圖。

圖4 遙感反演雪深分布圖（2019-01-08）Fig.4 Distribution of retrieved snow depth（2019-01-08）

在具備實(shí)測(cè)雪深數(shù)據(jù)的2019、2020年，選取積雪期內(nèi)的MODIS 反演雪深產(chǎn)品，從中提取雪深站位置的數(shù)據(jù)，依據(jù)積雪等級(jí)和判識(shí)標(biāo)準(zhǔn)劃分表［24］，將其與圖2 中S1、S2 雪深站的實(shí)測(cè)雪深進(jìn)行對(duì)比分析，結(jié)果見(jiàn)表2。

表2 可見(jiàn)，2019、2020年遙感反演雪深數(shù)據(jù)的準(zhǔn)確率均高于85%，遙感反演的年積雪日平均雪深與實(shí)測(cè)平均值的誤差均小于10%，2020年誤差在2%以內(nèi)，表明雪深反演精度較好。

表2 遙感反演與實(shí)測(cè)雪深對(duì)比結(jié)果統(tǒng)計(jì)Tab.2 Comparison of the snow depth between RS retrieved and measured

在遙感反演雪深的基礎(chǔ)上，對(duì)雪深的空間分布信息進(jìn)行處理，將空間網(wǎng)格積雪數(shù)據(jù)與流域分區(qū)疊置，計(jì)算各分區(qū)的積雪相關(guān)信息，包括雪深面積分布、平均雪深、積雪覆蓋率等。此處的平均雪深指有積雪覆蓋區(qū)域的平均雪深，無(wú)積雪區(qū)不統(tǒng)計(jì)在內(nèi)，因此雪量應(yīng)表示為平均雪深和雪蓋率的函數(shù)，對(duì)于非穩(wěn)定積雪期的3-6月，僅依據(jù)平均雪深無(wú)法確定分區(qū)的雪量和積雪融水量。處理后流域各分區(qū)2016-2020年平均積雪深度和雪蓋率變化過(guò)程見(jiàn)圖5、圖6，由圖可看出，①雪深和雪蓋率呈現(xiàn)明顯的年際變化，其中雪蓋率的年際變化規(guī)律更明顯；②雪深、雪蓋率與高程均呈現(xiàn)一定的相關(guān)性，其中雪蓋率隨高程變化的相關(guān)性更強(qiáng)；③融雪期的雪蓋率為與高程相關(guān)的緩慢變化的過(guò)程，而積雪期的雪蓋率變化迅猛，且與高程關(guān)系不大；④較雪蓋率而言，積雪區(qū)內(nèi)不同高程帶的平均雪深增加了隨時(shí)間跳變得無(wú)序性。

圖5 各區(qū)2016-2020年平均雪深變化過(guò)程圖Fig.5 Variation process of average snow depth in 2016-2020

圖6 各區(qū)2016-2020年雪蓋率變化過(guò)程圖Fig.6 Variation process of snow cover in 2016-2020

3.3 模型參數(shù)確定

基于雪深的SRM 模型參數(shù)與原基于雪蓋的SRM 模型基本相同，參數(shù)確定方式也類似，部分參數(shù)采用RS、GIS 手段或直接通過(guò)歷史資料推求得到［1，4-9］。

氣溫直減率可表示為氣溫隨高程變化的函數(shù)。在對(duì)應(yīng)的高程區(qū)域內(nèi)，根據(jù)不同高程測(cè)站的歷史氣溫?cái)?shù)據(jù)，采用γ =公式推求。其中γ 為氣溫直減率，表示每100 m高程氣溫下降的度數(shù)，單位為℃/100 m；h1、h2和T1、T2分別為兩個(gè)測(cè)站的高程和歷史平均氣溫。依據(jù)2016-2019年實(shí)測(cè)氣溫?cái)?shù)據(jù)計(jì)算各高程帶的氣溫直減率，當(dāng)1 個(gè)高程帶包含2 個(gè)以上測(cè)站時(shí)，首先計(jì)算每相鄰高程兩測(cè)站間的氣溫直減率，再依據(jù)高程差進(jìn)行加權(quán)計(jì)算，最終計(jì)算所得各高程帶氣溫直減率見(jiàn)表3。

表3 不同高程帶氣溫直減率Tab.3 The lapse rate of air temperature in elevation bands

無(wú)降雨和融雪徑流情況下，退水系數(shù)可表示為徑流量隨時(shí)間變化的函數(shù)。退水系數(shù)為第n天的實(shí)測(cè)徑流量，x、y為需要確定的常數(shù)，可根據(jù)退水段的歷史實(shí)測(cè)流量資料，構(gòu)建基于Qn和Qn+1的雙對(duì)數(shù)散點(diǎn)圖推求得到［4，5］。依據(jù)A站2016-2019年退水期日徑流繪制退水過(guò)程散點(diǎn)圖（圖7），根據(jù)擬合公式選取兩點(diǎn)，Q1=100 m3/s，Q2=250 m3/s 代入上式計(jì)算得x=906，y=-0.003 8。

圖7 A站2016-2019年退水過(guò)程散點(diǎn)圖Fig.7 Scatter plot of recession flow for A Station，2016-2019

基于雪深的SRM 模型中的融雪徑流系數(shù)CS、降雨徑流系數(shù)CR和度日因子α無(wú)法通過(guò)歷史資料直接確定，需要將降水、氣溫、蒸發(fā)和雪深的時(shí)間序列數(shù)據(jù)作為模型輸入，在設(shè)定初始參數(shù)的情況下，通過(guò)模型計(jì)算得到各高程帶逐日雪深和流域出口的流量過(guò)程，經(jīng)過(guò)與實(shí)測(cè)雪深和徑流對(duì)比，進(jìn)行參數(shù)的不斷調(diào)試，使計(jì)算和實(shí)測(cè)流量過(guò)程擬合最優(yōu)，從而確定最終參數(shù)［29，30］。

3.4 模擬精度分析

以MODIS 反演雪深為初始值，采用基于雪深的SRM 模型計(jì)算逐日雪深和徑流過(guò)程。以徑流誤差、Nash系數(shù)R2作為計(jì)算徑流與實(shí)測(cè)徑流吻合程度的評(píng)價(jià)指標(biāo)，以反演雪深作為計(jì)算雪深的校核標(biāo)準(zhǔn)，對(duì)A流域2016-2020年徑流過(guò)程進(jìn)行模擬，其中2016-2019年為模型率定期，2020年為檢驗(yàn)期。2016-2020年A 流域徑流模擬結(jié)果見(jiàn)表4（2020年數(shù)據(jù)截至12月10日），2019、2020年徑流模擬過(guò)程見(jiàn)圖8。結(jié)果表明，基于雪深的融雪徑流模擬精度較高，基本反映了流域徑流過(guò)程，其中2016-2020年模擬年徑流誤差均小于7%，各年的徑流擬合Nash系數(shù)都大于0.85。

表4 徑流模擬結(jié)果Tab.4 Simulation results of snowmelt runoff

圖8 模擬徑流與實(shí)測(cè)徑流對(duì)比圖Fig.8 Comparison of simulated runoff and measured runoff

融雪徑流計(jì)算中雪深和雪蓋率的準(zhǔn)確性直接影響融雪徑流的模擬精度，因雪深空間分布的不均勻性難以刻畫(huà)，本文依據(jù)雪深的空間分布計(jì)算區(qū)域雪蓋率，并在融雪徑流模擬中依據(jù)遙感反演信息對(duì)計(jì)算雪蓋率進(jìn)行校核。圖9為A流域不同高程帶2019、2020年3-6月融雪期的遙感反演和模擬計(jì)算雪蓋率變化過(guò)程，該圖清晰反映了高程對(duì)積雪覆蓋率的影響，同時(shí)可看出經(jīng)校核的逐日計(jì)算雪深較遙感雪深更精細(xì)地反映雪蓋率隨時(shí)間的變化。

圖9 各區(qū)計(jì)算雪蓋率與反演雪蓋率對(duì)比圖Fig.9 Comparison of calculated snow cover and retrieved snow cover

4 結(jié) 論

為了減少M(fèi)ODIS 合成雪蓋產(chǎn)品因時(shí)間分辨率不足造成的融雪徑流預(yù)報(bào)誤差，本文以新疆A流域?yàn)檠芯繉?duì)象，在建立MO?DIS積雪深度回歸模型、反演積雪深度的基礎(chǔ)上，將空間分布的雪深信息作為SRM 模型的輸入，并在模型中增加逐日雪深和雪蓋率計(jì)算模塊，最終以實(shí)測(cè)徑流和雪深數(shù)據(jù)作為融雪徑流模擬校核的依據(jù)，模擬計(jì)算A流域2016-2020年雪深、積雪覆蓋和徑流的逐日變化過(guò)程，確定SRM模型參數(shù)。

A流域MODIS反演雪深與2個(gè)實(shí)測(cè)站雪深相比，2019、2020年準(zhǔn)確率分別為85.29%和88.00%，遙感反演的日平均雪深與實(shí)測(cè)雪深的誤差小于10%?；谘┥畹娜谘搅髂M在2016-2020年模擬年徑流誤差均小于7%，最高的為2020年檢驗(yàn)期，誤差小于1%；同時(shí)每年的徑流擬合Nash系數(shù)都大于0.85。A 流域試驗(yàn)應(yīng)用結(jié)果表明，MODIS 反演雪深和基于雪深的融雪徑流模擬精度較高，MODIS 反演雪深有助于準(zhǔn)確刻畫(huà)流域積雪的時(shí)間、空間分布情況，采用基于雪深的SRM 模型計(jì)算，可更精細(xì)地反映積雪隨時(shí)間的變化過(guò)程，同時(shí)促進(jìn)融雪徑流預(yù)報(bào)精度的提升。

積雪信息的時(shí)空分辨率和準(zhǔn)確性是影響融雪徑流預(yù)報(bào)的關(guān)鍵，因受云層影響，基于MODIS 信息反演的積雪產(chǎn)品時(shí)間分辨率較低，而基于被動(dòng)微波遙感的積雪產(chǎn)品空間分辨率低，將兩者進(jìn)行時(shí)空尺度上的融合可為融雪徑流預(yù)報(bào)提供更精確的積雪信息數(shù)據(jù)源，將是遙感和融雪徑流預(yù)報(bào)的發(fā)展方向。 □