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

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

（1.國家能源集團新疆吉林臺水電開發有限公司，新疆伊犁835000；2.南瑞集團有限公司/國網電力科學研究院，南京211106）

新疆地區河流大多以冰雪融水為主要補給水源，其春、秋季來水以冰雪融水為主，夏季洪水多為融雪與暴雨的疊加。對于此類地區而言，融雪徑流是河流的基礎性來水，準確的融雪徑流預報可以提高水資源利用效率，有效預防洪澇災害［1-3］。

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

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

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

1 研究區概況

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

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

流域內布設有水文自動測報系統，包括5 個雨量氣溫站（A1～A5）、2個雨量氣溫雪深站（S1、S2）和1個水庫站，雨量氣溫站大多布設在河谷地區，測站分布見圖2，圖中測站旁標注數字表示測站位置的高程。采集的水文氣象數據為2016-2020年，數據時段長為日；其中僅2019、2020年有實測雪深數據。

圖2 A流域測站分布圖Fig.2 Station distribution in A basin

2 原理方法

2.1 MODIS反演雪深

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

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

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

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

2.2 基于雪深的SRM模型

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

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

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

雪深計算首先需判別當日（第n日）氣溫是否超過臨界溫度：低于臨界溫度時，降水通過雪水當量轉化為雪深，無融雪，雪深增加；高于臨界溫度時，積雪消融，雪深減小。積雪深度計算公式為：

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

降雪深度和消融的積雪深度均通過相應的水量換算得到，計算公式為：

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

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

3 結果分析

3.1 流域分區

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

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

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

3.2 MODIS反演雪深分析

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

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

在具備實測雪深數據的2019、2020年，選取積雪期內的MODIS 反演雪深產品，從中提取雪深站位置的數據，依據積雪等級和判識標準劃分表［24］，將其與圖2 中S1、S2 雪深站的實測雪深進行對比分析，結果見表2。

表2 可見，2019、2020年遙感反演雪深數據的準確率均高于85%，遙感反演的年積雪日平均雪深與實測平均值的誤差均小于10%，2020年誤差在2%以內，表明雪深反演精度較好。

表2 遙感反演與實測雪深對比結果統計Tab.2 Comparison of the snow depth between RS retrieved and measured

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

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

圖6 各區2016-2020年雪蓋率變化過程圖Fig.6 Variation process of snow cover in 2016-2020

3.3 模型參數確定

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

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

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

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

圖7 A站2016-2019年退水過程散點圖Fig.7 Scatter plot of recession flow for A Station，2016-2019

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

3.4 模擬精度分析

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

表4 徑流模擬結果Tab.4 Simulation results of snowmelt runoff

圖8 模擬徑流與實測徑流對比圖Fig.8 Comparison of simulated runoff and measured runoff

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

圖9 各區計算雪蓋率與反演雪蓋率對比圖Fig.9 Comparison of calculated snow cover and retrieved snow cover

4 結 論

為了減少MODIS 合成雪蓋產品因時間分辨率不足造成的融雪徑流預報誤差，本文以新疆A流域為研究對象，在建立MO?DIS積雪深度回歸模型、反演積雪深度的基礎上，將空間分布的雪深信息作為SRM 模型的輸入，并在模型中增加逐日雪深和雪蓋率計算模塊，最終以實測徑流和雪深數據作為融雪徑流模擬校核的依據，模擬計算A流域2016-2020年雪深、積雪覆蓋和徑流的逐日變化過程，確定SRM模型參數。

A流域MODIS反演雪深與2個實測站雪深相比，2019、2020年準確率分別為85.29%和88.00%，遙感反演的日平均雪深與實測雪深的誤差小于10%。基于雪深的融雪徑流模擬在2016-2020年模擬年徑流誤差均小于7%，最高的為2020年檢驗期，誤差小于1%；同時每年的徑流擬合Nash系數都大于0.85。A 流域試驗應用結果表明，MODIS 反演雪深和基于雪深的融雪徑流模擬精度較高，MODIS 反演雪深有助于準確刻畫流域積雪的時間、空間分布情況，采用基于雪深的SRM 模型計算，可更精細地反映積雪隨時間的變化過程，同時促進融雪徑流預報精度的提升。

積雪信息的時空分辨率和準確性是影響融雪徑流預報的關鍵，因受云層影響，基于MODIS 信息反演的積雪產品時間分辨率較低，而基于被動微波遙感的積雪產品空間分辨率低，將兩者進行時空尺度上的融合可為融雪徑流預報提供更精確的積雪信息數據源，將是遙感和融雪徑流預報的發展方向。 □