基于MODIS数据的青海省积雪覆盖范围监测算法探索

王雪璐，王玮，冯琦胜，吕志邦，梁天刚

（兰州大学草地农业科技学院 草地农业生态系统国家重点实验室，甘肃 兰州730020）

雪盖的时空动态分布是全球大气模拟的关键因子，对全球水循环和气候变化的研究起到至关重要的作用［1］。青海高原地区是我国主要的积雪分布区之一，同时也是我国六大牧区之一，雪灾是该地区冬春季最主要的自然灾害，严重影响草地畜牧业的可持续发展。而遥感技术以其视野广、时间分辨率高、信息量大、客观真实性强等优点，可实时监测大范围雪灾的成灾过程。因此，提高对该地区积雪范围的监测精度，对牧区防灾减灾具有特别重要的意义［2－6］。

积雪在可见光波段具有较高的反射率，在短波红外波段有较强的吸收特性，而大多数云在可见光波段和短段红外波段均具有较高的反射率。基于这种特殊的光谱特性，研究发现归一化差值雪被指数 （normalized difference snow index，NDSI）不仅可以识别出积雪，而且还是一个理想的云雪识别器，因此可以使用中分辨率成像光谱仪（moderate resolution imaging spectroradiometer，MODIS）第4波段（0.545～0.565μm）和第6波段（1.628～1.652μm）反射率计算的NDSI，作为监测积雪的重要指标［7］。国内外在MODIS积雪产品精度评价方面，已开展了大量的研究工作。在北美、欧洲和我国北疆等地的研究结果表明，在晴空条件下MOD10A1产品的积雪判识精度在80%以上。Liang等［8］研究结果表明北疆地区的该产品积雪分类精度为98.2%；Keith和Amir［9］在亨伯河流域的研究结果表明该产品2008年1－3月雪盖比例在75.08%～93.87%。然而，在我国青藏高原和高山地区的积雪分类精度相对较低，存在少分积雪而多分陆地的现象。张学通［10］研究表明晴空条件下青海省MOD10A1积雪分类精度仅为43.5%。曹云刚等［11］研究表明由于青藏高原山区地表坡度较大、地表粗糙度高、地表覆盖类型复杂等原因，容易出现斑块状积雪，因此严重影响了MOD10A1的积雪监测精度。郝晓华等［12］利用Landsat ETM＋以祁连山中部地区为例对MOD10A1产品的NDSI阈值进行了检验，结果表明该产品采用的NDSI阈值0.4偏高，造成了山区积雪面积低估的现象。本研究以青海省为例，利用美国国家雪冰数据中心（National Snow and Ice Data Center，NSIDC）的 MODIS地表反射产品 MOD09GA，结合青海省气象台站数据，探索研究了适合青海高原积雪监测的NDSI阈值，以期为准确监测青藏高原地区积雪分布范围提供科学依据。

1 材料与方法

1.1 研究区概况

青海省位于青藏高原的东北部，北纬31°39′～36°12′，东经89°45′～102°23′，东西长度约1 200km，南北宽约800km，总面积72.23万km2。东部地区地势起伏多山，但海拔较低，西部地区地势平缓，主要由高原和盆地结构组成。全省高原地区面积占总面积的80%以上，平均海拔在3 000m以上，海拔3 000m以下区域仅占全省总面积的26.3%；海拔3 000～5 000m的地区占地面积最大，约为全省面积的67%；5 000m以上的地区仅占全省面积的5%。全省依地形可分为祁连山地、柴达木盆地和青南高原三区。气候属典型的高原大陆性气候，年平均气温为－5.6～8.6℃、年降水量为15～750mm。天然草地以草甸植被为主，可利用草场占全省总国土面积的43.8%，其中高寒草甸和高寒草原总面积达2.448 16×107hm2，占全省草地面积的80.88%，其余部分为少量荒漠草地植被［13－15］。雪灾是该省境内主要的自然灾害之一，主要分布在青南高原地区，包括玉树、果洛两州及海南、黄南两州的南部，其中玉树、杂多、曲麻莱、清水河、玛多、达日高海拔区域雪灾尤为频繁，且范围广，影响大，对草地畜牧业所造成的损失也最严重［16］。

1.2 研究数据

通过NSIDC网站，订购并下载了2007年10月－2008年3月的每日MOD09GA（Surface Reflectance Daily L2GGlobal 1km and 500m）和 MOD10A1（MODIS／Terra Snow Cover Daily L3Global 500mGrid）产品。覆盖研究区的图像轨道号分别为h25v05和h26v05，共183d，共计732幅。MOD09GA和MOD10A1产品均以HDF－EOS分层数据格式（hierarchical data format，HDF）保存［17］。同时，为了开展积雪监测精度研究，收集了青海省2007年10月－2008年3月共计183d47个台站的气象观测数据，包括台站经纬度、海拔、雪深、平均气温等资料。在积雪观测时，雪深以整数登记，将小于0.5cm的积雪记录为无雪。共计8 601条记录。

1.3 积雪图像的分类方法

积雪分类算法的核心是归一化雪被指数的计算，此外还要排除水体和一些暗物体的影响，因此利用MOD09GA资料判别积雪的公式如下：

式中，Band2、Band4和Band6分别代表MOD09GA地表反射产品第2波段、第4波段和第6波段的反射率值，其空间分辨率为500m。

在NSIDC发布的全球MODIS雪被产品（如MOD10A1等）中，积雪制图算法主要包括以下3个条件：当NDSI≥0.4时，认为可识别出雪盖面积占50%以上的单个像元，为积雪判识的主要条件之一。水体在可见光波段和短波近红外波段的反射率和积雪类似，但可通过其第2波段的反射率小于11%的特性将其排除；同时通过暗物质在第4波段反射率小于10%这一约束，可防止一些可见光波段反射率极低的暗色物体被误判为积雪。

在森林覆盖区，由于植被对雪盖反射率的影响，需要对积雪判别条件进行相应的调整，将满足0.1＜NDSI＜0.4且 NDVI＞0.1的像元也识别为有雪［4］。

式中，Band1和Band2分别代表MOD09GA地表反射产品第1波段和第2波段的反射率值，其空间分辨率为500m。

1.4 研究区NDSI阈值的确定

影像的积雪分类精度是确定NDSI阈值的关键。依据青海省47个气象台站的实际观测值，比较分析不同NDSI阈值条件下的积雪分类图像精度，确定适合研究区积雪分类的NDSI阈值。积雪分类图像的精度主要通过积雪分类精度R（%）、总精度P（%）、漏测误差D（%）和多测误差V（%）这些指标进行综合评价。其计算公式如下：

式中，K代表影像和地面台站均有雪的样本数，E代表地面台站有雪的样本数，S代表影像和地面台站均无雪的样本数，T代表地面台站总样本数，M代表影像数据无雪但地面台站有雪的样本数，N代表影像数据有雪但地面台站无雪的样本数，Z代表地面无雪样本数。

1.5 自定义积雪分类图像的合成

对MOD09GA和MOD10A1数据进行预处理，分别得到MODIS的1，2，4，6波段反射率数据。依据公式（1）计算研究区逐日NDSI图像，采用不同的NDSI阈值分别对积雪进行分类，并结合NDVI计算林区雪盖，使用band2＞0.11和band4＞0.1两个阈值消除水体和暗物质的干扰，利用MOD10A1产品确定云覆盖范围，合成不同NDSI阈值的积雪分类图像（图1）。合成图像的编码方式和MOD10A1积雪分类图像一致，即25为陆地，50为云，200为积雪。

图1 用户自定义积雪分类图像合成流程图Fig.1 Flowchart for composition of user－defined snow classification image

2 结果与分析

2.1 MOD10A1积雪分类精度分析

青海省2007年10月－2008年3月的积雪季共有183个时相的MOD10A1积雪分类产品，结合47个气象台站观测的8 601条实测雪情记录数据，对该产品进行精度分析的结果表明，MOD10A1产品在晴空状况下的积雪分类精度仅为48.97%，总精度为95.56%。当雪深在1～3cm时，积雪漏测误差为72.65%，受其影响积雪分类精度仅为27.35%，总精度为95.85%；当雪深大于3cm时，积雪分类精度为86.01%（表1），总精度为99.21%。由此可见，研究区MOD10A1积雪分类精度明显低于该产品在其他研究区的精度验证结果。由于MOD10A1积雪分类产品是以MODIS L1B数据为基础，通过大气校正、几何纠正，结合云掩膜产品和水体掩膜产品，利用归一化雪被指数NDSI的阈值来判识积雪。在NSIDC发布的MODIS雪盖产品中，通用的NDSI阈值为0.40。该阈值是Hall和Riggs［18］在美国阿拉斯加、加利福尼亚、蒙大拿、明尼苏达地区利用TM雪盖图验证确定的，他指出部分区域的NDSI阈值是否合理仍然需要试验。因此，有必要对青藏高原地区积雪分类图像的NDSI阈值进行调整，以提高该地区的积雪分类精度。

表1 MOD10A1产品的积雪分类精度Table 1 Snow classification accuracies of composite image MOD10A1for different snow depths

2.2 NDSI阈值比较分析

已有研究表明，当雪深介于0.5～3cm时，MODIS传感器的积雪识别精度较低［8］。黄晓东等［19］对北疆牧区MODIS积雪制图MOD10A1精度评价结果表明，积雪深度1～3cm时积雪分类精度为7.5%。张学通［10］对青海省2002年10月－2008年3月的6个积雪季MOD10A1精度评价研究结果表明，1～3cm厚度的积雪分类精度仅为29.93%。由于MODIS图像的空间分辨率较低（500m），部分地区的地表坡度较大、常有裸露的石块和少量植被等地物干扰，所以积雪深度小于3cm时，地表难以被积雪均匀地覆盖，导致MODSI的积雪分类精度偏低。因此，本研究在探讨青海省积雪分类精度时，主要利用雪深大于3cm的台站积雪实测样本数据，分析确定该研究区的最佳NDSI阈值。

为了更精准地监测青海省雪盖变化情况，提高积雪分类图像的精度，本研究分别合成了NDSI≥0.50，0.49，……，0.21和0.20时的雪被分类图像，结合气象台站实测数据，分别统计了台站有雪图像也有雪、台站有雪而图像无雪、台站无雪但图像有雪、台站无雪图像也无雪的样本数，并计算了相关精度评价指标。研究结果表明，晴空下当雪深大于3cm时，随着NDSI阈值的减小，积雪分类精度和多测误差均逐渐增高，漏测误差逐渐减小。即NDSI阈值越大，识别出来的有雪像元越少。当NDSI阈值小于0.37时，较多非雪像元被误判为有雪像元，而当NDSI阈值大于0.44时，许多有雪像元不能被识别出来。因此，该研究区NDSI的理想阈值应该在0.37～0.44。分析总精度在不同NDSI阈值下的变化曲线，发现当NDSI阈值分别为0.44，0.43，0.38，0.37时，雪被分类图像的总精度最高，均为99.51%（表2）。当NDSI阈值为0.38和0.37时，积雪分类精度明显高于当NDSI阈值为0.44和0.43时的积雪分类精度，而相应的漏测误差也降低了近2%，多测误差基本保持不变（只上升了0.02%）。对NDSI≥0.38和NDSI≥0.37的雪被图像进行精度比较发现，当NDSI≥0.37时，其积雪分类精度比NDSI≥0.38的积雪分类精度高出0.7%，且漏测误差低于NDSI≥0.38时的雪被图像，所以最终确定NDSI≥0.37作为青海省积雪分类的阈值。

为了更全面地评价自定义NDSI阈值的合理性，本研究分析了NDSI≥0.37时雪被图像的雪深分段精度。结果表明，用户自定义NDSI阈值的雪被图像积雪分类精度也具有随积雪深度的增大而提高的特点（表3）。当雪深在1～3cm时，雪被图像的积雪分类精度相对较低，仅为43.56%。这主要是由于研究区内雪深较浅时NDSI雪被算法受其他地物的影响所造成的。当雪深大于3cm时，积雪分类精度可达90.37%，总精度为99.51%。由此可见，研究区NDSI阈值确定为0.37时，可以有效地提高积雪分类精度和总精度。

2.3 自定义图像与MOD10A1的精度比较

利用青海省2007年10月—2008年3月的MOD10A1雪被产品及相应的气象台站实测数据，比较分析了自定义合成雪被图像同MODIS每日雪被产品的积雪分类精度。研究结果表明，自定义雪被产品的积雪分类精度为62.31%，总精度为97.18%，而 MOD10A1产品的积雪分类精度为48.97%，总精度为95.56%。同MOD10A1产品相比较，自定义雪被图像的积雪分类精度和总精度都有不同程度的提高，尤其是积雪分类精度可提高13%，多测误差从0.36%降为0.22%，漏测误差从52.12%降至37.69%。因此，改进后的方法无论在漏测与多测误差上还是分类精度上都好于NSIDC的MODIS雪被产品，更适合于青海高原的积雪分类和动态监测。

对比分析在不同雪深情况下MOD10A1雪被产品和用户自定义图像的积雪分类精度表明，当气象台站记录雪深为0cm时，即实测值为陆地的情况下，自定义雪被图像和MOD10A1雪被产品的多测误差分别为0.22%和0.36%。这说明NDSI阈值的适度下调并没有对多测误差和总精度造成不良的影响。当雪深为1～3cm时，自定义雪被图像和MOD10A1产品的积雪漏测样本数分别是114和178，漏测误差分别为56.44%和72.65%。自定义图像积雪分类精度为43.56%，MOD10A1积雪分类精度为27.35%。说明自定义图像在对浅雪（雪深为1～3cm）区域积雪漏测的状况有明显改善。随着积雪深度的增加，自定义雪被图像和MOD10A1雪被图像的积雪分类精度都在提高，当雪深在4～6cm时，自定义雪被图像的积雪分类精度为89.09%，明显高于MOD10A1的积雪分类精度76.92%。当雪深≥7cm时，其积雪分类精度分别为91.25%和91.36%（表4）。

表2 不同NDSI阈值下的积雪判识精度评价结果Table 2 Accuracy and errors for snow cover images under different NDSI thresholds %

表3 NDSI≥0.37时合成图像的积雪分类精度Table 3 Snow classification accuracies of user－defined snow cover images under NDSI≥0.37for different snow depths

表4 不同雪深下自定义产品与MOD10A1产品的积雪分类精度Table 4 Snow classification accuracies of user－defined snow cover and MOD10A1images for different snow depths

3 结论与讨论

本研究通过对青海省范围内MOD10A1雪盖产品精度分析，发现其积雪分类精度并不理想，在雪深大于3 cm时的晴空条件下为80.61%；同时结合青海省气象台站的实测雪情数据，比较了不同NDSI阈值雪被图像的精度，确定出晴空下雪深大于3cm时NDSI阈值为0.37的雪被图像精度最高，其积雪分类精度可达90.37%，总精度99.51%，多测误差0.22%，漏测误差为9.63%；最后在比较了MOD10A1产品和自定义NDSI≥0.37的积雪分类图像精度后发现，雪深分段统计的结果表明，当雪深＜7cm时雪被图像的分类精度较MOD10A1产品有很大的改善作用；当雪深≥7cm时，积雪分类精度无明显变化。

由此可见，通过下调NDSI阈值至0.37能够较好地提高积雪识别的准确率，特别是当雪深＜7cm时，积雪分类精度较MOD10A1产品有明显的提高。这是由于在薄雪地区，积雪容易受地表其他地物的影响，很难完整地覆盖图像中的整个像元，导致用来判识积雪的反射率特征不够明显。因此，下调NDSI域值可以使更多的浅雪像元识别出来，从而提高浅层积雪和山区积雪的识别精度；当雪深≥7cm时，自定义雪被产品的积雪分类精度同MODIS每日雪被产品的精度基本相同，都高于90%。说明本研究雪被算法的改进并未对≥7cm的积雪识别起到明显效果。这可能是由于本研究是基于NDSI阈值的积雪分类算法改进，以调整NDSI阈值的方法来提高积雪识别精度。而这种方法是通过降低或提升积雪判别条件来实现的，因此对雪盖反射率特征显著的深层积雪影响不大。

总体来说，本研究通过改进青海省积雪覆盖度图像算法，使自定义积雪分类图像的积雪精度明显高于MOD10A1产品，尤其是提高了对浅层积雪的识别率，不仅适用于提取该研究区内的雪被图像，而且更加符合我国青藏高原地区的实际情况，对该研究区的积雪监测具有重要意义。

该项研究存在的一些不足和有待进一步解决的问题主要有以下几点。首先，精度评价准确性受气象台站分布位置限制。本研究所利用的积雪实测数据全部来自于青海省内的47个气象台站，其中大部分台站分布在青海省的东北部，导致实测数据在地理上的分布不均，不能够准确而全面地代表整个研究区的情况，从而影响了研究的最终结果。其次，由于MOD09GA数据量较大，从订购下载到保存再到数据的处理的过程比较复杂，因此本研究只选择了2007—2008年这一积雪季对该区域NDSI阈值进行了初步的探索。在以后的研究中，可扩展MODIS资料的时空分布范围，如开展多个积雪季的整个青藏高原地区的研究，检验自定义NDSI阈值的合理性，从而提高研究结果的准确性。

［1］梁继，张新焕，王建.基于NDVI背景场的雪盖制图算法探索［J］.遥感学报，2007，11（1）：85－93.

［2］张娟，冯蜀青，徐维新.MODIS数据在青海省积雪监测中的应用［J］.青海气象，2006，1：55－57.

［3］周陆生，李海红，汪青春.青藏高原东部牧区大—暴雪过程及雪灾分布的基本特征［J］.高原气象，2000，19（4）：480－458.

［4］王丽红，付培建，鲁安新.遥感技术在牧区雪灾监测研究中的应用［J］.遥感技术与应用，1998，13（2）：32－36.

［5］刘兴元，龙瑞军，尚占环.草地生态系统服务功能及其价值评估方法研究［J］.草业学报，2011，20（1）：167－174.

［6］韩立辉，尚占环，任国华，等.青藏高原“黑土滩”退化草地植物和土壤对秃斑面积变化的响应［J］.草业学报，2011，20（1）：1－6.

［7］Hall D K，Riggs G A，Salomonson V V.Algorithm Theoretical Basis Document（ATBD）for the MODIS Snow and Sea Ice－Mapping Algorithms［DB／OL］.［2011－11－25］.http：／／modis.gsfc.nasa.gov／data／atbd／atbd＿mod10.pdf.2001.

［8］Liang G，Huang D，Wu X.An application of MODIS data to snow cover monitoring in a pastoral area：A case study in Northern Xinjiang，China［J］.Remote Sensing of Environment，2008，112：1514－1526.

［9］Keith N A，Amir A K.Investigation of MODIS Snow Mapping in A Boreal Forest Watershed［C］.Baltimore：ASPRS，2009.

［10］张学通.青海省积雪监测与青南牧区雪灾预警研究［D］.兰州：兰州大学，2010.

［11］曹云刚，杨秀春，徐斌，等.MODIS在青藏高原大范围积雪制图中的应用及存在的问题［J］.科技导报，2007，25（21）：51－54.

［12］郝晓华，王建，李弘毅.MODIS雪盖制图中NDSI阈值的检验——以祁连山中部山区为例［J］.冰川冻土，2008，30（1）：132－138.

［13］张国胜，伏洋，杨琼，等.青海省天然草地类型空间分布特征及气候分区［J］.草业科学，2009，26（1）：23－29.

［14］陈钟，冯琦胜，黄晓东，等.青海省主要类型天然草地EVI动态变化［J］.草业科学，2010，27（5）：23－29.

［15］李东，黄耀，吴琴，等.青藏高原高寒草甸生态系统土壤有机碳动态模拟研究［J］.草业学报，2010，19（2）：160－168.

［16］徐裕华.高原牧区雪冻灾害的气候特征［J］.四川气象，1998，（2）：28－32.

［17］Riggs G A，Hall D K，Salomonson V V.MODIS Snow Products User Guide Collection 5［DB／OL］.［2011－01－20］.http：／／modis－snow－ice.gsfc.nasa.gov／uploads／sug＿c5.pdf.2006.

［18］Hall D K，Riggs G A.Accuracy assessment of the MODIS snow－cover products［J］.Hydrological Processes，2007，21：1534－1547.

［19］黄晓东，张学通，李霞，等.北疆牧区MODIS积雪产品MOD10A1和MOD10A2的精度分析与评价［J］.冰川冻土，2007，29（5）：722－729.