基于多源數據的松嫩平原黑土區亞像元雪蓋率算法研究

王子龍 胡石濤 付 強 姜秋香 印玉明

(東北農業大學水利與土木工程學院， 哈爾濱 150030)

0 引言

積雪是重要的地表覆蓋物之一，北半球冬季的雪蓋面積可達5×107km2，積雪覆蓋率達到30%以上[1-2]。積雪的時空分布及其參數特性對高緯度、高海拔地區的氣象、水文狀況有顯著影響，并且對全球水資源利用、大氣循環、農作物墑情預報和氣候變暖等研究具有重要意義[3]，同時積雪還是一種至關重要的淡水資源[4]。因此，完善積雪監測系統，提高積雪面積制圖精度，實現對積雪覆蓋信息的量化研究，對我國寒區農業種植及經濟社會發展具有獨特作用[5]。

目前，國內外學者針對不同的地表環境開展了一系列雪蓋制圖研究，通用的算法主要是二值分類法和亞像元制圖法[6]。二值分類法是基于積雪與其他地物不同的光譜特性，借助設定NDSI閾值、決策樹分類等方法，最終將影像中地物區分為積雪或非雪。HALL等[7]最先提出SNOMAP算法，利用Landsat/TM數據將地物進行分類并完成雪蓋制圖。但是該方法只能將像元內地物分為積雪或非雪像元，導致雪蓋反演中存在明顯的高估或低估現象，因此很難滿足高精度反演研究的需求。

為了解決普遍存在的像元內地物混合問題，研究者提出了混合像元分解法[6,8-10]和統計回歸法[11-13]，實現了亞像元雪蓋制圖。通過混合像元分解完成雪蓋反演的算法雖然極大地提高了積雪制圖精度，但其要求提取“純凈端元”并應用最小二乘法計算端元組分比例，計算過程復雜，很難在大范圍內進行實踐，從而限制了此算法的推廣[14]。統計回歸法通過分析NDSI與積雪覆蓋率之間的關系，建立線性回歸模型反演逐像元雪蓋率[1]，借助區域內通用模型，利用MODIS數據的NDSI值反演像元雪蓋率，提高大范圍區域內像元雪蓋率的計算效率，有利于推廣應用。SALOMOSON等[13]將美國、俄羅斯、加拿大境內的3個地區作為研究區，基于NDSI值與由Landsat/ETM+數據提取的真實雪蓋率之間的線性回歸關系建立反演模型。國內學者也開展了不同地區的相關研究，張穎等[15]針對MOD10A1數據在青藏高原雪蓋率反演中存在精度較低等問題，利用MODIS地表反射率數據建立分段模型反演雪蓋率，提高了反演精度；曹云剛等[3]基于MODIS數據建立雪蓋率與NDSI、歸一化植被指數(Normalized differential vegetation index，NDVI)等因子間的回歸模型，改進了反演模型。

鑒于不同的地理環境會對雪蓋率與NDSI之間關系的穩定性產生影響，為了進一步提高局部區域內雪蓋率估算的精度，本文改進SALOMOSON模型，建立逐像元雪蓋率估算模型，以松嫩平原黑土區為研究區，將Landsat8 OLI影像提取的積雪像元假設為地表真實雪蓋，進而逐像元建立MODIS像元雪蓋率與相應NDSI值間的回歸關系模型，并與FSC數據進行對比分析，驗證本研究反演方法提取積雪信息的優越性。結果將作為計算雪水當量的輸入數據，以期為實現寒區農田土壤春熵預報提供數據支持。

1 研究區概況與數據集

1.1 研究區概況

研究區為松嫩平原黑土區，松嫩平原土壤肥沃、地勢平坦，是黑龍江省糧食主產區和我國商品糧生產基地；松嫩平原黑土區是我國東北黑土區的主要組成部分[16]，作為松嫩平原內集中連片的黑土區域，具有良好的地理典型性和區域代表性。圖1為研究區高程圖與地表反射率數據MOD09GA 6、4、2波段合成影像，其中藍色部分為雪蓋范圍。研究區所在的東北地區是我國三大季節性積雪區域之一，是典型的高緯積雪帶，范圍為44°～52°N、124°～130°E。

圖1 研究區高程圖與MOD09GA 6、4、2波段合成影像Fig.1 Elevation map and MOD09GA 6, 4 and 2 band synthesis images of study region

1.2 數據集

1.2.1MODIS數據

MODIS是搭載于美國對地觀測系統(Earth observation system，EOS)Terra和Aqua衛星上的主要傳感器之一，光譜范圍0.4～14.4 μm，共有36個光譜波段，具有較高的時空分辨率。在綜合考慮穩定雪蓋期、避免云層影響、傳感器周期及數據代表性等因素的基礎上，本文收集了2016年12月17日的4景無云MODIS/Terra地表反射率數據(MOD09GA)反演積雪覆蓋率，此時距初雪日已有一個月間隔，地表雪蓋對反映同時期研究區內的整體降雪情況具有較好的代表性。另外選用了4景同時相的MOD10A1 FSC數據與反演結果進行對比驗證，兩種數據均來源于美國國家雪冰數據中心(National snow and ice data center，NSIDC)，版本為V005，分辨率為500 m，數據格式為HDF，一般投影為正弦曲線投影[17]。

1.2.2OLI數據

OLI數據波段信息如表1所示，相對于MODIS數據，具有更高的空間分辨率，可以獲取更加詳細的地表積雪參數信息，還包括了對積雪信息識別極其重要的短波紅外(Shortwave infrared，SWIR)波段，有利于中小范圍的積雪反演研究[18]，同時，具有較高分辨率的OLI影像往往被當做驗證低分辨率影像(如MODIS數據)的“真值”數據。數據可以從美國地質勘探局(United states geological survey，USGS)官網下載。本文主要收集了11景與MOD10A1 FSC、MOD09GA數據同時相的Landsat8 OLI數據。

表1 Landsat8 OLI數據波段信息Tab.1 Band information of Landsat8 OLI data

2 研究方法

2.1 數據預處理

2.1.1MODIS數據

MOD09GA原始數據已經過輻射定標、大氣校正等處理[19]，利用MODIS數據處理軟件MRT(MODIS reprojection tools)對獲得的原始數據進行格式轉換、坐標變換，文件格式輸出為Geo-TIFF格式，將原投影轉換成基準面為WGS84的通用墨卡托(Universal transverse mercator projection，UTM)投影系統，空間分辨率為500 m，重采樣選用最鄰近法，提取影像的反射率值，利用研究區“松嫩平原黑土區”的行政邊界矢量數據進行裁剪得到研究區范圍內的影像數據，最后通過ENVI軟件計算得到NDSI；針對MOD10A1數據，預處理手段與MOD09GA數據一致，最終輸出FSC數據[20]，然后根據FSC數據集的像元編碼意義(表2)，提取積雪覆蓋率信息，由于積雪面積比例數據中包含多類無意義編碼，文中只提取(0,100]范圍內的像元值，其余值(無效值)不參與計算。NDSI計算公式為

NDSI=(b4-b6)/(b4+b6)

(1)

式中b4——MODIS數據第4波段反射率

b6——MODIS數據第6波段反射率

表2 MOD10A1 FSC像元編碼及含意Tab.2 Pixel coding and meaning of MOD10A1 FSC

2.1.2OLI數據

OLI數據與其他Landsat系列數據類似，標示為L1T級，原始影像格式默認為Geo-TIFF格式，投影系統默認是基準面為WGS84的UTM投影系統，數據已經過地形參與的幾何精校正[21]。由于本文定量反演積雪信息需要利用地表反射率數據，因此需要對OLI數據進行大氣校正處理。針對OLI數據的預處理主要包括大氣校正、輻射校正、影像拼接、裁剪等。通過輻射校正將原始影像的像元亮度(Digital number，DN)轉換為大氣頂層表觀輻射亮度，之后需要對影像進行大氣校正處理，將輻射校正結果轉換為地表反射率，去除大氣因素對影像的影響。

L=ad+b

(2)

式中L——轉換得到的大氣頂層表觀輻射亮度

d——影像DN值

a——增益系數b——偏移系數

2.2 逐像元估算模型建立

2.2.1OLI數據積雪信息提取

針對Landsat系列數據，一般利用SNOMAP算法識別積雪信息，該算法的關鍵是NDSI閾值法。對于OLI數據，首先計算NDSI

NDSI=(RVIS-RSWIR)/(RVIS+RSWIR)

(3)

式中RVIS——可見光波段反射率

RSWIR——SWIR波段反射率

通過設定NDSI閾值并利用決策樹分類法完成積雪與非積雪地物的區分，OLI數據中的第3、6波段分別對應式(3)中的可見光和短波紅外波段。

根據HALL等[7]研究，SNOMAP算法中將NDSI閾值設為0.4，能夠有效地區分云霧并識別積雪信息，即滿足像元NDSI大于等于0.4，且同時滿足R2>0.10、R4>0.11(R2、R4分別為MODIS數據第2、4波段的反射率)時，將該像元識別為積雪像元。該閾值在NSIDC發布的MODIS全球雪蓋產品中作為通用閾值，并且在國內不同研究區中得到了驗證[15,17,22]，同時考慮到本研究區內沒有相應結論可作參考的現實情況，本文也將此閾值作為識別積雪像元信息的標準。本文采用經驗的決策樹分類方法分別設定R3>0.10、R5>0.11(R3、R5分別為OLI數據第3、5波段的反射率)兩個附加條件以避免暗物質及結冰水體被誤識為積雪[15]，生成分辨率為30 m的二值積雪分類圖。

2.2.2雪蓋率“真值”計算

本文提出的真實雪蓋率是指任一MODIS像元內二值圖中積雪像元所占的比例。以高分辨率為30 m的OLI數據提取得到的積雪數據作為積雪像元“真值”，MODIS數據經重采樣后分辨率為480 m，從而每個MODIS像元內包含對應的256個OLI像元。對每個MODIS像元內包含的二值積雪像元數進行統計，計算相應MODIS像元的雪蓋率(Snow cover fraction，FRA)，針對二值圖與對應MODIS像元疊合時邊緣處存在的零碎(非整個)像元，最后通過所占面積比例統計，計算公式為

(4)

式中nsnow——MODIS像元內包含二值圖中積雪像元的數目

N——MODIS像元內二值圖包含的像元個數，為256

利用Create Fishnet工具建立與像元大小一致網格，進一步計算出每個MODIS像元對應網格內包含的OLI數據值為1(即雪像元)的個數，最后利用式(4)統計得到每個MODIS像元對應的FRA值。

2.2.3模型建立

積雪混合像元通常是指低分辨率影像識別的積雪像元內并非完全是積雪，而是包括積雪和非積雪兩類地物信息的現象。SALOMOSON等[13]基于統計回歸方法建立反映MODIS影像的NDSI與對應像元內雪蓋率之間關系的線性回歸模型。模型分為2種，分別為：

模型1

NDSI=a1FRA+b1

(5)

模型2

FRA=a2NDSI+b2

(6)

式中a1、b1、a2、b2——模型參數

由于受到“Errors-in-variables”問題[23-24]的影響，以上2個模型最后構建的線性方程不能重合。因此，需要探討選用哪個模型更為合理。在線性回歸模型的構建中往往是假設自變量值固定，因變量值隨自變量變化并受到隨機因素影響而存在相應誤差[24]。因為MODIS影像的分辨率較低且計算得到的NDSI值容易受到大氣、傳感器本身等因素的影響，利用OLI數據提取的MODIS雪蓋率真值比利用MODIS數據通過波段運算得到的NDSI值更為固定，因此，本文選用模型1建立回歸模型。最后對比MOD10A1 FSC數據與反演模型結果，通過計算與基于OLI影像獲得的MODIS雪蓋率真值的均方根誤差及平均絕對誤差，分析反演模型與標準產品的精度。

3 反演算法驗證與精度分析

對比OLI數據提取的雪蓋率“真值”與FSC數據之間的相關關系，表明FSC數據僅能反映出松嫩平原黑土區的地表積雪的基本分布狀況，在山區及城市邊緣地區兩者積雪分布存在較大差別(圖2)，同時其易受到自然因素(云、霧霾等)的影響導致部分地區的積雪數據異常，從而不能準確反映研究區地表真實的積雪面積等參數信息，很難滿足當前積雪參數反演研究的高精度要求，因此，需要建立能夠更好適應特定地區雪蓋率反演研究的算法及模型。

圖2 MOD10A1 FSC數據與OLI影像提取的積雪分類圖Fig.2 Snow classification maps of MOD10A1 FSC data and OLI image

圖3為MODIS各像元的NDSI值與雪蓋率“真值”之間的散點圖。由圖可見雪蓋率值主要集中在2個區間(右上角高值區和左下角低值區)，高值區(雪蓋率大于0.8)數據集聚是由于區域內存在連續的大面積雪蓋，同時在季節性積雪區往往會將結冰后河流湖泊誤識為地表雪蓋，從而高估了雪蓋率；而低值區(雪蓋率小于0.1)出現集中主要是由于大量非積雪像元的存在[25]，大面積的山區、林地也導致地表積雪覆蓋不均勻，此外還有光學傳感器極易受到自然因素影響的原因。本文采用模型1建立線性模擬方程，可以在一定程度上改善模型2存在的在極值區高估、低估雪蓋率的問題，同樣也存在低估高值區積雪系數的問題[16]。研究認為：在雪蓋率高值區，當雪蓋率值超過一定限度后，對NDSI值影響最大的因素是雪粒徑等積雪性質，而不再是積雪的覆蓋率；對于低值區，受復雜地表覆蓋物等因素的影響，NDSI與FRA之間的相關度呈現下降趨勢。因此，通過NDSI值與像元雪蓋率的相關關系建立回歸模型反演得到的積雪覆蓋率“估值”會存在一定的誤差，之后的研究中需要針對此問題不斷進行完善。

圖3 NDSI值與雪蓋率“真值”散點圖Fig.3 Scatter diagram of NDSI and FRA

為了評估反演模型的精度并對其進行對比驗證，本文選擇將反演估算結果與同時相的FSC數據進行對比，將得到的反演模型估算結果、FSC數據分別與基于OLI數據提取的雪蓋率真值進行誤差統計分析(表3)并對模型估算值的精度進行驗證，結果表明：①FSC數據在松嫩平原黑土區精度較低，平均雪蓋率為80.21%，與同時相OLI影像的平均雪蓋率(87.71%)相差較大，兩者之間的相關系數僅為0.58。②利用MOD09GA數據建立的亞像元雪蓋率反演模型得到的平均雪蓋率為85.28%，與同時相的OLI影像較為接近，且兩者相關系數為0.66，相對于FSC數據來說，雪蓋率與相關度有了明顯提高。③本文反演模型的估算結果與FSC數據相比，誤差統計結果(均方根誤差、平均絕對誤差)均有所降低。其中，對亞像元反演模型得到的估算結果進行統計時，雪蓋率大于1的則將值重新賦值為1，而雪蓋率小于0的則賦值為0[19]。

表3 反演模型與FSC數據的誤差統計Tab.3 Error statistics of inversion model and FSC data

4 結論

(1)借鑒SALOMOSON模型對松嫩平原黑土區建立逐像元的雪蓋率反演模型，利用MOD09GA、OLI數據參與計算和分析，對模型反演結果進行了精度驗證和不足分析。研究發現，利用OLI數據作為數據源進行亞像元雪蓋率研究是可行的，提高了研究區的雪蓋率反演精度，解決了二值分類法存在的雪蓋率估算誤差較大等問題。

(2)SNOMAP算法利用NDSI閾值法進行積雪像元識別時，將閾值設為0.4作為判識標準在研究區內是可行的，結果基本符合地表雪蓋現狀，可為之后的積雪參數反演研究提供參考。

(3)基于反演算法得到的雪蓋率與雪蓋率“真值”數據相比，相對于MOD10A1 FSC數據在精度上有了明顯提高，誤差相對減小，在一定程度上滿足當前大范圍雪蓋率反演精度的要求，對該地區雪蓋監測提供數據支持。

1 劉良明,徐琪,胡玥,等. 利用非線性NDSI模型進行積雪覆蓋率反演研究[J]. 武漢大學學報：信息科學版,2012,37(5):534-536.

LIU Liangming,XU Qi,HU Yue,et al. Estimating fractional snow cover based on nonlinear NDSI model[J]. Geomatics and Information Science of Wuhan University,2012,37(5):534-536.(in Chinese)

2 ROBINSON D A, DEWEY K F, HEIM J R R. Global snow cover monitoring: an update[J]. Bulletin of the American Meteorological Society, 1993, 74(9): 1689-1696.

3 曹云剛,劉闖.一種簡化的MODIS亞像元積雪信息提取方法[J].冰川凍土,2006,28(4):562-567.

CAO Yungang,LIU Chuang. A simplified algorithm for extracting subpixel snow cover information from MODIS data[J]. Journal of Glaciology &Geocryology, 2006, 28(4):562-567.(in Chinese)

4 劉海,陳曉玲,宋珍,等. MODIS影像雪深遙感反演特征參數選擇與模型研究[J]. 武漢大學學報：信息科學版,2011,36(1):113-116,121.

LIU Hai,CHEN Xiaoling,SONG Zhen,et al. Study of characteristic parametric selection and model construction for snow depth retrieval from MODIS image[J].Geomaticsand Information Science of Wuhan University,2011, 36(1):113-116,121.(in Chinese)

5 何詠琪,黃曉東,方金,等.基于HJ-1B衛星數據的積雪面積制圖算法研究[J].冰川凍土,2013,35(1): 65-73.

HE Yongqi,HUANG Xiaodong,FANG Jin,et al. Snow cover mapping algorithm based on HJ-1B satellite data[J]. Journal of Glaciology & Geocryology, 2013, 35(1):65-73.(in Chinese)

6 施建成. MODIS亞像元積雪覆蓋反演算法研究[J]. 第四紀研究,2012,32(1):6-15.

SHI Jiancheng.An automatic algorithm on estimating subpixel snow cover from MODIS[J].Quaternary Sciences,2012,32(1):6-15.(in Chinese)

7 HALL D K,RIGGS G A, SALOMOSON V V. Development of methods for mapping global snow cover using moderate resolution imaging spectroradiometer data[J]. Remote Sensing of Environment, 1995, 54(2): 127-140.

8 ROSENTHAL W,DOZIER J. Automated mapping of montane snow cover at subpixel resolution from the Landsat Thematic Mapper[J]. Water Resources Research, 1996, 32(1): 115-130.

9 PAINTER T H,DOZIER J,ROBERTS D A,et al. Retrieval of subpixel snow-covered area and grain size from imaging spectrometer data[J]. Remote Sensing of Environment, 2003, 85(1): 64-77.

10 PAINTER T H,RITTGER K,MCKENZIE C,et al. Retrieval of subpixel snow covered area, grain size, and albedo from MODIS[J]. Remote Sensing of Environment, 2009, 113(4): 868-879.

11 KAUFMAN Y J,KLEIDMAN R G,HALL D K,et al. Remote sensing of subpixel snow cover using 0.66 and 2.1 μm channels[J]. Geophysical Research Letters, 2002, 29(16):28-1-28-4.

12 BARTON J S,HALL D K,RIGGS G A. Remote sensing of fractional snow cover using moderate resolution imaging spectroradiometer (MODIS) data[C]∥Proceedings of the 57th Eastern Snow Conference, 2000: 171-183.

13 SALOMOSON V V,APPEL I. Estimating fractional snow cover from MODIS using the normalized difference snow index[J]. Remote Sensing of Environment, 2004, 89(3): 351-360.

14 唐志光,王建,彥立利,等.基于MODIS的青藏高原亞像元積雪覆蓋反演[J].干旱區資源與環境,2013,27(11):33-38.

TANG Zhigang,WANG Jian,YAN Lili,et al. Estimating sub-pixel snow cover from MODIS in Qinghai-Tibet Plateau[J]. Journal of Arid Land Resources & Environment, 2013, 27(11):33-38.(in Chinese)

15 張穎,黃曉東,王瑋,等. MODIS逐日積雪覆蓋率產品驗證及算法重建[J]. 干旱區研究,2013,30(5):808-814.

ZHANG Ying,HUANG Xiaodong,WANG Wei,et al. Validation and algorithm redevelopment of MODIS daily fractional snow cover products[J]. Arid Zone Research, 2013, 30(5):808-814.(in Chinese)

16 劉繼龍,任高奇,付強,等. 松嫩平原黑土區玉米穗質量構成要素的空間變異性研究[J/OL]. 農業機械學報,2016,47(12):178-184,222. http:∥www.j-csam.org/jcsam/ch/reader/view_abstract.aspx?flag=1&file_no=20161222&journal_id=jcsam.DOI:10.6041/j.issn.1000-1298.2016.12.022.

LIU Jilong,REN Gaoqi,FU Qiang,et al.Spatial variability of components of corn ear weight in black soil region of Songnen Plain[J/OL].Transactions of the Chinese Society for Agricultural Machinery,2016,47(12):178-184,222.(in Chinese)

17 李云,馮學智,肖鵬峰,等. 巴音布魯克典型區MODIS亞像元積雪覆蓋率估算[J]. 南京大學學報：自然科學,2015,51(5):1022-1029.

LI Yun,FENG Xuezhi,XIAO Pengfeng,et al.Estimating per-pixel snow cover fraction from MODIS in typical area of Bayanbulak[J].Journal of Nanjing University:Natural Sciences,2015,51(5):1022-1029.(in Chinese)

18 王曉艷,王建,李弘毅,等. NDSI與NDFSI結合的山區林地積雪制圖方法[J]. 遙感學報,2017,21(2):310-317.

WANG Xiaoyan,WANG Jian,LI Hongyi,et al.Combination of NDSI and NDFSI for snow cover mapping in Mountainous-Forested[J].Journal of Remote Sensing,2017,21(2):310-317.(in Chinese)

19 王杰,黃春林,郝曉華.一種考慮雪粒徑變化的積雪面積反演算法[J].地球信息科學學報,2017,19(1):101-109.

WANG Jie,HUANG Chunlin,HAO Xiaohua.An algorithm of snow cover fraction retrieval considering the variability of snow particle size[J].Journal of Geo-information Science,2017,19(1):101-109.(in Chinese)

20 梁天剛,高新華,黃曉東,等.新疆北部MODIS積雪制圖算法的分類精度[J].干旱區研究,2007,24(4):446-452.

LIANG Tiangang,GAO Xinhua,HUANG Xiaodong,et al. Study on the accuracy of MODIS snow cover mapping algorithm in northern Xinjiang[J]. Arid Zone Research, 2007, 24(4):446-452.(in Chinese)

21 樊曉兵. 復雜山區MOD10A1積雪面積數據精度變化研究[D].成都：西南交通大學,2016.

FAN Xiaobing.The study of accuracy variation of MOD10A1 snow cover datas in Rugged Terrain[D].Chengdu:Southwest Jiaotong University,2016.(in Chinese)

22 鄧婕. 基于多源遙感資料的中國積雪制圖及其時空變化研究[D].蘭州：蘭州大學,2016.

DENG Jie.Snow mapping and its dynamics using multi-source remote sensing data in China[D].Lanzhou: Lanzhou University,2016.(in Chinese)

23 FULLER W A．Measurement error models[M]．New York: Wiley，1987: 30-59.

24 時正華,袁永生. Errors-in-variables模型的參數估計[J].曲阜師范大學學報：自然科學版,2005,31(1):35-37.

SHI Zhenghua,YUAN Yongsheng. The parameter estimation of errors-in-variables model[J]. Journal of Qufu Normal University:Natural Sciences, 2005,31(1):35-37.(in Chinese)

25 周強,王世新,周藝,等.MODIS亞像元積雪覆蓋率提取方法[J].中國科學院研究生院學報,2009,26(3):383-388.

ZHOU Qiang,WANG Shixin,ZHOU Yi,et al. Algorithm for MODIS subpixel snow fraction[J]. Journal of the Graduate School of the Chinese Academy of Sciences, 2009,26(3):383-388.(in Chinese)