高山复杂地形区域MODIS积雪产品的云移除算法
2020-10-20马强谢绍璞张爽
马强 谢绍璞 张爽
摘 要：免费可用的MODIS积雪覆盖区域图（SCA）因其较高的时空分辨率以及积雪/云的辨别能力一直被广泛用于积雪覆盖变化的研究。然而，云遮挡阻碍了MODIS积雪产品的广泛应用，尤其是在高山区域。在高山复杂地形区域，空间局部异质性相对较高并且在原始的MODIS SCA图中，云覆盖通常是连续的和持续的，尤其是在冬季。基于以上的问题，文章通过使用Terra和Aqua日结合、时间移动平均（TMA）推测、极值雪线法以及K最近邻 （KNN）算法对高山复杂地形区域的云覆盖进行移除，从而产生无云的日MODIS积雪覆盖图。在相对验证中，文章所提出的云移除算法的总体精度达到83.12%。因此，文章中提出的除云算法能够有效的去除高山复杂地形区域中原始MODIS积雪产品的云覆盖。
关键词：MODIS积雪数据;云移除;高山区域;复杂地形;时间移动平均推测;瞬时区域
中图分类号：P426.635 文献标志码：A 文章编号：2095-2945（2020）29-0001-07
Abstract： The freely available MODIS Snow Cover Area （SCA） maps have been widely used to investigate snow coverage changes because of their higher spatio-temporal resolution and snow/cloud discrimination capabilities. However， cloud obscuration hinders their extensive application， especially in mountain regions. In alpine areas with complex to pography， spatial local heterogeneity is relatively higher and cloud coverage in original MODIS SCA maps is frequently continuous and persistent， especially in winter. Aimed at problems mentioned above， Terra and Aqua combination， Temporal Moving-Average（TMA） deduction， and extreme snow-line method and K-Nearest Neighbor （KNN） algorithm were used to remove clouds from MODIS SCA data in Alpine areas with complex to pography and yield cloud-free MODIS SCA data. In relative validation， overall accuracy of the cloud removal methodology proposed in this paper was 83.12%. Therefore， the solution proposed can effectively eliminate clouds from original MODIS snow cover data in alpine areas with complex to pography.
Keywords： MODIS snow data; cloud removal; alpine area; complex to pography; Temporal Moving-Average （TMA） deduction; transient zones
1 概述
由于實地测量的局限性，遥感被广泛用于偏远地区或者高山区域从而探测积雪变化;同时遥感提供了多尺度的积雪观测，这对于监测积雪场的演化是非常有用的信息。积雪覆盖是积雪场的一个重要因子。MODIS二值积雪覆盖图因其高的时间分辨率（一天）、中尺度的空间分辨率（500m）以及积雪/云的辨别能力，从而在众多卫星衍生的积雪覆盖产品中具有明显的优势[2-3]。然而，云遮挡阻碍了其广泛的应用，尤其是在高山区域。
在高山复杂地形区域，积雪累积和消融不仅受到海拔的影响，而且因为不同的太阳辐射、能量平衡以及迎风/背风等因素使得坡度和坡向也会对积雪积累和消融产生影响。这些因素会导致积雪在高山区域是非连续性分布的并且造成较大的空间局部异质性。例如，在冬季坡度较陡的高海拔区域其地表覆盖是陆地，但是毗邻于该陆地区域的低洼平坦地区其地表覆盖却是积雪。即使将坡度、坡向和海拔这些地形因素融入到基于临近像素的空间滤波中，该方法的除云精度和效果依然比较差。同时，高山区域的云覆盖通常是持续的，从而使得时间滤波的精度或者除云效果将会降低。此外，经过时空结合方法处理后MODIS 积雪覆盖图剩余的云主要集中在瞬时区域。然而，瞬时区域的积雪场是极其不稳定的，尤其是在冬季以外的其他季节。因此，基于长时间序列分析的积雪循环方法不能够有效的解决这个问题。
针对上述问题，本文的目的是提出一个有效的方法对高山复杂地形区域中MODIS日积雪产品的云覆盖进行去除并且产生无云的MODIS日积雪覆盖数据。
2研究区和数据源
2.1 研究区
台兰河流域位于中国新疆乌鲁木齐自治区温宿县境内（图1）。以台兰河水文站控制的流域面积总共是1324 km2，其空间范围是从东经80°04′E -80°35′E和北纬41°35′N-42°05′N。在30m DEM数据中所提取的流域面积是1361km2，因此，本研究中所认定的流域面积也是1361km2。流域海拔的变化范围是从1577到7421m asl，其平均海拔是3739m asl。在地形上，大约96%的区域高于2000m asl，73%的区域高于2000m asl以及46%的区域高于4000m asl。在山麓平原区，植被类型主要以荒漠草场和灌丛为主;在中海拔区域分布着大片的云杉和亚高山草甸;在高海拔区域普遍发育着永久性积雪覆盖和冰川冰。
流域受到地形抬升、西风环流和北冰洋水汽的影响，造成大气水汽传送从西到东、从北到南逐渐减少。因此该研究区属于典型的干旱性大陆气候，并且春夏两季高山区域的冰雪融水是台兰河的主要补给源。冰雪融水约占到年径流量的60%以上，雨水的贡献所占比例相对较小并且主要降雨发生在暴雨时期。
2.2 MODIS数据
本研究的除云算法是对MODIS日积雪覆盖数据MOD10A1（V005）和MYD10A1（V005）进行处理。通过访问美国国家冰雪信息中心（National Snow and Ice Data Center， NSIDC）网址，可以获取免费的MODIS日积雪覆盖数据，该数据采用的正弦投影以HDF的格式发布。整个研究区需要两个MODIS块（H23V04和H24V04）进行拼接才能对其完全覆盖。我们使用的MODIS积雪覆盖数据其日期范围是从2007年1月1日开始到2011年12月31日截止。1820幅MOD10A1影像和1825幅MYD10A1影像将会用于该研究（表1）。缺失的影像将会用完全被云覆盖的影像补充。MODIS重投影工具（MODIS Reprojection Tool， MRT） 将用于数字拼接并且将HDF格式转化成GeoTIF格式;转化后生成的影像将会被重投影到WGS1984 UTM 43N带，并且采用最近邻插值将重投影的影像重采样成500m。最后，重采样影像的11种像素类别会被重新分为3类，“湖冰”类别将会被纳入到“积雪”类別;“湖泊”和“海洋”类别将会划归到“非积雪”（陆地）类别;“云”和其他种类将会划分到“云”类别。
2.3 Landsat数据
Landsat TM/ETM+有六个空间分辨率为30m的波段影像，通过访问USGS数据服务器（http：//earthexplorer.usgs.gov/）可以进行免费下载。我们选择了13幅Landsat TM/ETM+影像与同日Terra卫星所提供的MOD10A1积雪产品进行多传感器结合，从而生成非填充验证数据对本文中所描述的除云算法进行精度评估（表2）。Landsat TM/ETM+影像的选取标准如下：（1）Landsat TM/ETM+影像是无云的或者云覆盖极少的;（2）同日的原始MODIS积雪产品（MOD10A1和MYD10A1）的云覆盖高于30%。
因为Landsat TM和ETM+传感器具有相似的光谱波段，积雪制图算法生成二值积雪覆盖图的时候会把它们当成相同的传感器。因此，在本文中，Landsat TM和ETM+将会被简称为“Landsats”。Landsats影像预处理的过程如下：（1）首先将Landsats影像重投影到WGS1984 UTM 43N带;（2）为了能够适应积雪制图算法，将Landsats波段影像初始的DN值转化为大气顶部（Top of the Atmosphere ，TOA）反射率。积雪在可见光波段（Landsats band2）有高的反射率并且在近红外波段（Landsats band5）是接近于0的反射率，积雪制图算法利用这些光谱反射特性探测积雪。归一化积雪指数（Normalized Difference Snow Index， NDSI）的计算公式如下：
一个二值的积雪覆盖图通过设定NDSI≥0.40、近红外波段 （Landsat Band 4/MODIS Band 2）的反射率大于0.11和可见光波段（Landsats Band 3） 的反射率大于0.10而得到。然而，积雪制图算法通常会将阴影区域的积雪错分为陆地即产生低估误差。Macander对误差的原因进行了详细的分析。为了避免这种低估误差，我们对Landsats影像中的阴影进行消除，然后再使用积雪制图算法获取Landsats二值积雪覆盖图。最后，将获取的积雪覆盖图通过判定一个500m栅格内其对应的30m积雪像素所占比例，将其重采样成500m。判定条件如下：如果一个500m栅格内其对应的30m积雪像素所占比例超过50%，则将该栅格像素划分为积雪;如果所占比例小于50%，则将该像素划分为陆地;否则，该500m像素将会被指定为云像素。
2.4 数字高程模型 （Digital Elevation Model， DEM）
ASTER GDEM 版本2数据（GDEM 2）由ASTER GDEM团队研发并且自2011年11月17日公布使用。数据来源于中国科学院计算机网络信息中心地理空间数据云平台（http：//www. gscloud.cn）。整个研究区域需要两景GDEM（N41E080和N42E080）才能够被完全覆盖。首先对两景GDEM进行拼接，然后将其重投影到WGS1984 UTM 43N带，最后重投影的数据将会用于台兰河流域的提取。随后使用双线性插值方法将流域DEM图由30m重采样成500m，同时利用重采样的DEM图获取一个坡向图。
3 方法
3.1 除云算法
如概述所说，本文的目的是提出一个有效的除云算法用于去除高山复杂地形区域中MODIS积雪产品的云覆盖。该方法由四个步骤组成：Terra和Aqua的日结合、时间移动平均（Temporal Moving-Average， TMA）推测、极值雪线法和K最近临（K-Nearest Neighbor， KNN）算法。方法按顺序依次执行每个步骤，每个步骤的输出是后续步骤的输入。整个方法通过Python语言实现。
3.1.1 日结合
步骤1是将MOD10A1（V005版本）和MYD10A1（V005版本）提供的MODIS / Terra和MODIS / Aqua日积雪覆盖产品进行日结合。在本研究中，我们使用MOD10A1优先的合成规则，即MOD10A1中的云覆盖像素将会被MYD10A1中的无云像素替代。
3.1.2 时间移动平均推测
步骤2是基于积雪覆盖的周期性变化对云遮挡下真实的地表覆盖进行恢复。一个区域的积雪覆盖通常从冬季至夏季逐渐减少，从夏季至冬季逐渐增加，呈现出区间变化。在一个水文年内，积雪覆盖面积呈现出四个区间变化：单调递增区间、单调递减区间、凸区间和凹区间。在单调递减区间，当日降雪的概率大于右邻域降雪的概率;在单调递增区间，当日降雪的概率大于左邻域降雪的概率;在凸区间，当天降雪的概率大于左邻域和右邻域降雪的概率;然而，在凹区间，当日降雪的概率小于左邻域和右邻域降雪的概率。通过以上总结，我们发现当日降雪的概率大于增区间、减区间和凸区间左邻域或右邻域的降雪概率，从而可以推测出当日云遮挡下真实的地表覆盖。我们使用上述的积雪区间变化特征，可以检测出位于凹区间的云覆盖像素并且对这些云像素不进行处理。
如图2所示，深青色的线为积雪消耗曲线（A），可以发现日积雪覆盖变化较大; 红线为积雪消耗月移动均值线（B），月移动平均积雪覆盖变化趋势比较稳定，从1月至8月其月移动平均积雪覆盖逐渐减少，8月至12月其月移动平均积雪覆盖逐渐增加，在8月初至8月中旬其平均积雪覆盖达到最低并且较为稳定;蓝线是决策线（C）并且通常低于红线，但是蓝线在凹区间要高于红线。决策线是前一个月和后一个月（31天）的月平均积雪覆盖的最小值确定得到的。
基于以上的理论研究，我们对云像素真实的地表覆盖进行推测。首先，对当日云像素前31天和后31天的降雪概率分别进行计算。降雪概率计算的公式如下：
NS是31天内积雪出现的个数; NL是31天内陆地出现的个数。当日云像素前31天和后31天积雪出现的概率分别标记为P-31→0和P0→+31。其次，选择P-31→0和 P0→+31之间的最小值，并且标记为Pmin。第三，计算当日云像素前15日至后15日积雪出现的概率，作为当日积雪出现概率并且标记为P-15→+15。最后，如果P-15→+15
3.1.3 极值雪线法
在本研究中，我们采用了Gafurov和Bárdossy提出的积雪过渡高程法并且对其进行了改进[17，21]。陆地线（Land Line）仍是积雪像素中最低的海拔作为阈值，然而雪线（Snow Line）是陆地像素中最高的海拔作为阈值。通常，雪线高于陆地线;如果陆地线高于雪线，则表明雪线和陆地线之间的像素全部被云遮挡。因此，我们利用如下规则对位于极高海拔和极低海拔位置的云像素进行分类：
（1）如果陆地线低于雪线，则雪线以上的云像素变成积雪，陆地线以下的云像素变为陆地;
（2）如果陆地线高于雪线，则陆地线变为雪线，而雪线变为陆地线，将雪线以上的云像素指定为积雪像素，陆地线以下的云像素指定为陆地像素;
（3）如果只有雪线而没有陆地线，则将雪线变为陆地线，并将陆地线以下的云像素重新分类积雪像素;
（4）如果只有陆地线而没有雪线，则陆地线变为雪线，并且雪线以上的云像素重新划分为积雪像素。
3.1.4 最近邻（K-Nearest Neighbor， KNN）算法
在经过极值雪线 （步骤3） 处理后， 云覆盖主要集中在瞬时区域，在这个区域的积雪场是极其不稳定的。因此，使用极少的云像素或者一个积雪循环去估算瞬时区域的云像素是不适合的。为了解决这个问题，我们引入了机器学习中的K最近邻（KNN）算法。KNN算法是由Cover和Hart提出的，该算法的基本原理如下：在样本空间中，如果一个未知像素周围最邻近的K个已知像素主体是属于一个特定的类别，则这个未知样本将会被指定为这个特定类别。一个已知样本通常是已知样本和未知样本之间的欧式距离确定的。在三维空间中两个点的欧式距离计算公式如下：
其中y是像素行号索引;x是像素列号索引;t是像素日期索引。KNN算法的样本空间大小取决于样本的个数。为了确保参与估算云像素的无云像素个数充足，我们设定K=27。在高山区域，我们假定时间自相关大于空间自相关。这是因为在高山区域由于较大的地形起伏造成空间局部异质性相对较高，然而，积雪覆盖变化相对较低因为相对积雪覆盖位于相对较高的海拔。因此，必须确保参与估算云覆盖像素的无云像素具有较小空间距离。我们设置时间权重为2，然后将原始的时间轴缩小到一半，空间轴保持不变。在确定27个无云像素之后，将会对云像素的真实地表覆盖进行推测。如果积雪像素的个数大于13，则云覆盖像素将会被重新划分为积雪像素;否则，该云像素将会被指定为陆地像素。
3.2 验证和评估
在本研究中，我们引入一种新的验证方法，将其称作非填充验证。该方法使用Landsats影像对除云算法的精度进行评估。使用Landsats影像用作验证是基于以下的考虑：Landsat卫星的过境时间要早于MODIS卫星，并且在白天随着时间的推移云覆盖逐渐增加。这很有可能在同一天MODIS影像的含云量较高而Landsat卫星获取的波段影像是無云的或者是含云量极少的，这将为非填充验证提供了可能。非填充验证的“地表真实”是通过多传感器结合生成的。MOD10A1中的云像素将会被Landsats二值积雪覆盖图中相同位置的非云像素进行替代，从而生成非填充验证的“地表真实”。在非填充验证中，我们采用一致度 （Degree Of Agreement， DA）这个指标对本研究中所提出的算法精度进行评估。一致度（Degree Of Agreement， DA）是指经过除云算法处理正确分类的像素占总的分类像素的比例，该指标均采用百分比形式，具体计算公式如下：
式中，a和d分别是指正确分类为积雪和陆地的像素个数; c是指将积雪错分为陆地的像素个数; b是将陆地错分为积雪的像素个数。总体一致性（Overall Agreement， OA）是通过各个测试日不同表现性能的平均值得出的，对每个测试日对所检查像素正确率的百分比进行加权得到的，其计算公式为：
其中i代表验证日;RNi代表在原始的Terra影像（MOD10A1）中是云像素而在对应的“地表真实”中是非云像素的个数。
4 结果
本文中描述的四个除云步骤能够完全消除原始的MODIS积雪覆盖数据中的云覆盖。图3a显示了MOD10A1、MYD10A1和使用步骤1-4之后月平均云覆盖。MODIS日积雪覆盖产品（MOD10A1和MYD10A1）以及经过各个步骤后的月平均覆盖如图3b所示。如图3a所示，原始的MODIS积雪产品 （MDO10A1和MYD10A1） 中的月平均云覆盖与季节有关，在冬季的云覆盖达到最高并且在其他季节逐渐减少。MYD10A1的云覆盖高于MOD10A1中的云覆盖是因为Aqua卫星的过境时间晚于Terra卫星并且云覆盖随着白天时间的推移在逐渐增加。与原始的MOD10A1产品相比，经过日结合处理（步骤1）后，月平均云覆盖减少了7.06%-11.29%并且月积雪覆盖略有增加。当MOD10A1和MYD10A1的云结构不同时，日结合的除云效果较好。时间移动平均推测（步骤2）是最有效的步骤并且云覆盖又减少了26.58%-46.06%;同时，积雪覆盖有了明显增加，尤其是在雪季 （10月至3月）。极值雪线法 （步骤3） 能够有效的消除位于极低海拔和极高海拔持续性的云覆盖，并且不需要对经过步骤2处理后的影像设定特定的无云比例进行除云。在经过第3步处理后，云覆盖率降低了2.02%-7.41%。与除云效果相比，极值雪线法对积雪覆盖的恢复在大多数月份都是较少的。这主要是因为极值雪线法去除的云的位置相对固定，并且低海拔区域出现积雪的概率比较低。K最近邻算法（步骤4）能够完全去除影像中剩余的云覆盖，并且增加了更多的积雪。这表明该方法能够有效的解决经过先前三个步骤处理后位于瞬时区域剩余的、持续时间较长的云覆盖。图4是2009年3月6日原始的MODIS积雪产品以及经过各个除云步骤处理后剩余的云覆盖图。
5 验证
在非填充验证中，验证样本由13幅多传感器结合影像组成用于评估本文中所提出的除云算法各个步骤的精度以及除云算法的总体精度，这些影像几乎全部来自冬季（表3）。可以看出日结合（step 1）的总体一致性低于时间移动平均推测（step 2）和极值雪线法（step 3），其原因在3.1.1节进行了详细的分析。虽然移动时间平均推测的输入是精度相对较低的日结合影像，但是当日云像素的推测是基于积雪＼陆地出现的概率从而确保经过移动时间平均推测后日结合方法的不确定性误差能够被最大程度的消除并且获取可靠的结果。移动时间平均推测的云移除效果最好并且其精度最高。同时，该步骤的方差仅仅是0.18%，表明该方法是极其稳定的。极值雪线法的总体精度要比日结合的精度高，这是因为该方法消除的云覆盖位置固定的并且类型是确定的，而日结合去除的云覆盖较为离散且有一定的不确定性。然而，极值雪线法的方差是四个步骤中最高的，这可能是因为探测器饱和、云或云覆盖的影响造成生成的Landsat 积雪覆盖图在较高的海拔有较大的低估误差，尤其是在冬季。在经过上述三步处理后，剩余的云覆盖主要位于瞬时区域，这个区域的积雪场是极其不稳定。因此，该方法的精度较低是可预见的并且能够接受的。K最近邻算法（step 4）的精度最低，但是其方差要低于日结合。这表明更多的无云像素参与瞬时区域云像素的估算，能够确保获取稳定的结果。通过对本文所提出的除云算法进行整体评估，其总体精度达到83.12%，同时该算法的方差仅有0.21%并且原始的MODIS SCA图像中的云覆盖都被消除。因此，本文所提出的算法是适用于去除高山复杂地形区域MODIS积雪产品中的云覆盖。
TFC 和AFC分别代表Terra和Aqua影像中的云比例;RF代表经过各个步骤处理后Terra影像中的剩余的云量;step 1-4是指各云移除步骤， CRM是指由各个步骤组成的整个的云移除算法（Cloud Removal Methodology）。
6 讨论与结论
本文的研究目的是为了提出一个适用于高山复杂地形区域的除云算法。为了克服高山复杂区域连续性和持续性云覆盖以及地物分布空间局部异质性较大的问题，我们提出一个由四个步骤组成的除云算法并且原始的MODIS积雪产品（MOD10A1和MYD10A1） 的云覆盖在经过我们的方法处理后能够被完全的消除。第一个步骤是将Terra和Aqua卫星提供的MODIS日积雪覆盖图（MOD10A1和MYD10A1）进行日结合，该方法是利用了MODIS卫星一日两次获取数据的独特优势。第二个步骤是时间移动平均推测，该方法使用一个水文年内积雪覆盖变化去推测非瞬时区域的云覆盖并且在四个云移除步骤中精度最高、效果最好。同时，该方法具有最小的方差、线性时间复杂度并且易于实施。第三个步骤是极值雪线法，该方法使用陆地像素的最高海拔（雪线）和积雪像素的最低海拔（陆地线，Land Line）对位于极高海拔和极地海拔的云覆盖进行去除。这个方法能够有效的避免了Parajka 提出的区域雪线法将低洼平原区的积雪覆盖错分为陆地的低估误差，并且其除云效果相对较好。位于瞬时区域云覆蓋的去除问题一直是困扰除云算法的一个难题，为了解决这个难题我们在最后一步引入K最近邻（K-Nearest Neighbor， KNN）算法。KNN算法是一种相对简单的机器学习方法，它使用K个已知的最近邻样本来推断未知样本的特定类别。虽然KNN算法的精度相对较低，但是剩余的云覆盖能够被完全去除。
通过与其他除云算法相比，我们提出的解决方案具有最高的除云精度和最佳的除云效果。然而，使用本研究所提出的KNN算法对瞬时区域的云覆盖进行处理后所得结果的精度仍然较低，但是我们认为可能会有更好的机器学习算法从而达到更高的精度。
参考文献：
[1]ARSENAULT K R， HOUSER P R， DE LANNOY G J M. Evaluation of the MODIS snow cover fraction product[J]. Hydrological Processes， 2014， 28（3）： 980-998.
[2]HALL D K， RIGGS G A， SALOMONSON V V， et al. MODIS snow-cover products[J]. Remote Sensing of Environment， 2002， 83： 181-194.
[3]PARAJKA J， BLSCHL G. The value of MODIS snow cover data in validating and calibrating conceptual hydrologic models [J]. Journal of Hydrology， 2008， 358（3-4）： 240-258.
[4]AULT T W， CZAJKOWSKI K P， BENKO T， et al. Validation of the MODIS snow product and cloud mask using student and NWS cooperative station observations in the Lower Great Lakes Region[J]. Remote Sensing of Environment， 2006， 105（4）： 341-353.
[5]PARAJKA J， BLSCHL G. Validation of MODIS snow cover images over Austria[J]. Hydrology and Earth System Sciences Discussions， 2006， 10（5）： 1569-1601.
[6]KRAJí P， HOLKO L， PERDIGO R A P， et al. Estimation of regional snowline elevation （RSLE） from MODIS images for seasonally snow covered mountain basins[J]. Journal of Hydrology， 2014， 519：1769-1778.
[7]ZHOU X， XIE H， HENDRICKX J M H. Statistical evaluation of remotely sensed snow-cover products with constraints from streamflow and SNOTEL measurements[J]. Remote Sensing of Environment， 2005， 94（2）： 214-231.
[8]TONG J， DéRY S J， JACKSON P L. Topographic control of snow distribution in an alpine watershed of western Canada inferred from spatially-filtered MODIS snow products [J]. Hydrology and Earth System Sciences Discussions， 2008， 5（4）： 2347-2371.
[9]RIGGS G A， HALL D K， SALOMONSON V V. MODIS snow products user guide to collection 5[J/OL]. Boulder， Colorado USA： NASA National Snow and Ice Data Center Distributed Active Archive Center， 2006， https：//modis-snow-ice.gsfc.nasa.gov/？c=userguides.
[10]MACANDER M J， SWINGLEY C S， JOLY K， et al. Landsat-based snow persistence map for northwest Alaska [J]. Remote Sensing of Environment， 2015， 163： 23-31.
[11]Salomonson V V， Appel I. Development of the Aqua MODIS NDSI fractional snow cover algorithm and validation results [J]. IEEE Transactions on Geoscience and Remote Sensing， 2006，44 （7）： 1747-1756.
[12]Hall D K， Riggs G A. Accuracy assessment of the MODIS snow products[J]. Hydrological Processes，2007，21（12）：1534-1547.
[13]GAFUROV A， BáRDOSSY A. Cloud removal methodology from MODIS snow cover product [J]. Hydrology and Earth System Sciences， 2009， 13（7）： 1361-1373.
[14]COVER T M， HART P E. Nearest neighbor pattern classification [J]. IEEE Transactions on Information Theory， 1967， 13（1）： 21-27.
[15]KRAJí P， HOLKO L， PERDIGO R A P， et al. Estimation of regional snowline elevation （RSLE） from MODIS images for seasonally snow covered mountain basins[J]. Journal of Hydrology，2014，519：1769-1778.
[16]沈永平，劉时银，丁永建，等.天山南坡台兰河流域冰川物质平衡变化及其对径流的影响[J].冰川冻土，2003，25（2）：124-129.
[17]蒋艳，周成虎，程维明.阿克苏河流域径流补给及径流变化特征分析[J].自然资源学报，2005，20（1）：29-36.
[18]田龙.台兰河流域积雪覆盖面积变化特征分析[J].西北水电，2015（02）：13-16.
[19]谢昌卫，丁永建，刘时银，等.托木尔峰南坡冰川水文特征及其对径流的影响分析[J].干旱区地理，2004，27（4）：570-575.
[20]沈志宝，陈有虞.青藏高原平均云量的估算及其分布特征[J].高原气象，1987，6（4）：49-58.
[21]孙美萍，姚晓军，李忠勤，等.21世纪天山南坡台兰河流域径流变化情景预估[J].气候变化研究进展，2012，8（5）：342-349.
