基于多源数据的松嫩平原黑土区亚像元雪盖率算法研究

王子龙 胡石涛 付 强 姜秋香 印玉明

(东北农业大学水利与土木工程学院， 哈尔滨 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)