积雪作为最活跃的地面覆盖类型之一，通过影响地表辐射平衡和地气之间的水热交换来反馈气候变化，是冰冻圈内最为敏感的环境变化响应因子^[1-4]。同时，积雪覆盖也影响区域生态地理环境，其异常变化常导致如雪崩、干旱、洪灾等自然灾害，对当地自然生态环境、社会经济发展以及人类活动等造成严重影响^[4-5]。黑龙江流域位于东北亚地区，是全球重要的季节性积雪分布地区^[6]。黑龙江流域生态地理环境复杂，人类活动以及社会经济区域间差异明显，流域内河系发达，支流众多，融雪径流约占径流补给的15%~20%^[7]，积雪在黑龙江流域维持生态环境平衡、社会经济发展以及防治自然灾害等方面具有十分显著的现实作用。因此，在全球气候改变的大背景下，掌握黑龙江流域积雪时空动态变化与气候因子之间的关系是非常必要且有意义的。

针对大尺度范围内的积雪时空变化，国内外学者开展了一系列的相关研究。目前针对积雪时空变化研究的手段主要有常规的地面气象站点监测和卫星遥感监测^[8-9]。MODIS积雪产品因其具有约500 m的空间分辨率、每日重访一到两次的时间分辨率和良好的产品开发团队等优势，成为进行大尺度范围内积雪时空动态变化研究的重要数据源^[10-11]。于灵雪等^[4]基于MODIS8d合成的积雪产品对黑龙江流域的积雪覆盖面积进行提取并分析其时空变化，Wang等^[6]则基于MODIS逐日积雪覆盖产品在黑龙江流域开展积雪时空动态变化研究，并提取积雪日数、积雪初日以及积雪终日等积雪评价参数，结合土地覆盖数据探讨森林对积雪的影响。综上所述，MODIS积雪产品在黑龙江流域进行积雪时空动态研究应用已有一定的研究基础。雷小春等^[12]在黑龙江流域就MODIS积雪产品的精度进行分析，表明云是影响该产品精度的主要原因。因此在使用MODIS积雪产品进行分析应用之前，还需要对云进行处理，减少云对积雪产品精度的影响。目前，针对这一科学问题，已有一系列的成熟算法^[13-22]，但在黑龙江流域的应用较为薄弱，在实际应用过程中，需综合考虑黑龙江流域地理环境特征，做进一步的调整。黑龙江流域作为典型的季节性积雪区域，对全球气候变化响应十分敏感，但现阶段对黑龙江流域的积雪日数、积雪面积变化的综合研究相对薄弱，针对黑龙江流域的积雪时空动态变化与气候要素间的相关性研究仍有不足，特别是针对不同时期的积雪面积与气候因子之间的相关性有待进一步探讨。

本文以MODIS逐日积雪产品MOD10A1和MYD10A1为数据源，参考现有的去云算法并加以调整，对数据源进行预处理，得到逐日无云积雪覆盖数据集，基于该数据集分析黑龙江流域2002年9月—2016年8月共14 a的积雪日数和积雪面积的变化特征，并将积雪面积与气候因子相结合，分析不同时间段的积雪面积对气候变化的响应关系。研究结果旨在提高对积雪时空变化及其与气候因子之间相关关系的认识，更好地应对全球气候变化。

1 研究区概况

黑龙江流域(图 1，图中所含的行政数据来源于http://www.naturalearthdata.com)是一个流经中国、俄罗斯、蒙古以及朝鲜小部分地区的跨国流域，流域面积约2.08×10⁶ km²，是东北亚最大的流域之一。黑龙江流域地形特征复杂，海拔跨度为0~2 585 m，地貌整体上呈现西高东低的趋势，流域的西边主要以丘陵山地和平原为主，而中东部则主要以平原为主。流域属于寒温带与温带季风气候，受内陆及海岸气压与季风气候的影响，年平均气温由北纬的-7 ℃至南纬6 ℃不等，不同区域温差较大。受区域气候和地形地貌等因素的影响，流域降水较为集中，主要发生在6—9月，而旱季相应较长。降水量年际差异较大，且因地形条件和海洋影响，流域内部降水分布差异明显，年降水均值约500 mm^{[4, 7, 23-24]}。

2 数据及预处理 2.1 雪深数据和气象数据

本文收集到中国东北黑龙江省78个站点实测的2002—2010年共8年的雪深数据，用于对无云积雪覆盖数据集的精度评估，数据来源于黑龙江气候数据中心。文中涉及的中国气象数据来源于国家气象信息中心(http://data.cma.cn)，国外气象数据来源于美国国家海洋和大气管理局(https://www.climate.gov/)。雪深数据和气候数据等气象站点分布如图 1所示。

2.2 MODIS逐日积雪覆盖数据

美国国家雪冰数据中心免费提供的每日积雪覆盖产品MOD10A1 /MYD10A1为源数据，投影格式为正弦曲线投影SIN(sinusoidal projection)，分辨率约为500 m^[25]。MODIS系列遥感影像数据在全球范围内共按36个条带进行存储，黑龙江流域占其中6个条带(h24v03, h25v03, h25v04, h26v3, h26v04, h27v04)，时间序列选择为2002年9月1日—2016年8月31日。基于IDL编程环境，调用MRT软件将获得的MODIS积雪数据进行批量预处理，预处理主要包括拼接、投影转换、裁剪以及重编码，其中重编码规则见表 1。

2.3 数字高程模型数据

数字高程模型数据(digital elevation model，DEM)用于去云算法辅助研究，采用USDS EROS数据中心发布的90 m分辨率的SRTM数字高程模型数据(http://srtm.csi.cgiar.org)。为了匹配MODIS积雪产品分辨率，将DEM数据重采样为MOD10A1/MYD1A1同等分辨率。

2.4 逐日无云积雪覆盖数据集生成及精度验证

由于光学遥感无法穿透云层，无法获得云下的地表覆盖信息，因此MODIS逐日积雪覆盖产品存在大量的云像元干扰。据统计，黑龙江流域2002—2016年MOD10A1和MYD10A1影像的平均云含量百分比达到53%和60%，云含量超过50%的天数比例高达59%和74%左右。大量的云像元造成地表覆盖信息在空间上不连续，直接使用该数据做后续的研究，明显是不合理的，因此，使用该数据进行积雪研究，去云是必不可少的一部分。

参考文献[19]的去云算法，依次通过包括基于Terra星和Auqa星过境时间不同的双星数据合成^{[13, 15, 18-19]}、基于时间组合的连续3 d临近日分析^{[3, 17, 19, 21, 26]}、基于空间连续性的邻近4像元法^{[13-14, 17, 26]}、基于数字高程模型的Snowl方法^{[3, 13-14, 17, 19, 21-22]}以及高程滤波算法^{[3, 14, 19, 26]}、基于长时间序列分析的季节性过滤^{[13-14, 19, 21, 26]}等6步去云算法进行去云处理，其算法流程如图 2所示。

针对黑龙江流域的地理环境特点，主要对区域Snowl算法和季节性过滤算法作了调整。区域Snowl算法的使用对云像元含量有限制，并且当积雪空间异质性及破碎化情况比较严重时，也会受到一定的限制^[26]。常规的处理方式是进行海拔分带处理，黑龙江流域海拔集中在1 000 m以内，并且纬度跨度较大，流域范围广。采用海拔分带处理，易造成海拔相近纬度不同的区域的积雪误判。因此参考Paudel和Andersen^[13]的处理方法，基于Arc-swat模型将黑龙江流域细化为多个子流域，在子流域范围内运用Snowl算法，降低积雪误判的影响。季节性过滤是基于黑龙江流域为典型季节性积雪区域这一事实，对每个像元以每年的9月1日至次年的8月31日为时间序列进行分析。假定积雪稳定积累后不出现偶然融雪事件，研究区的每个像元主要以陆地或云为开始，直到出现积雪(即积雪开始累积)，随后春季积雪融化再次出现陆地或云。利用连续3 d滑动窗口过滤积雪初期和融雪末期的偶然积雪事件，将积雪稳定积累的时间作为积雪积累的阈值，将积雪完全融化后的时间为积雪融化的阈值，不考虑积雪期内的偶然降雪和融雪事件，将在阈值范围内的云像元判为积雪，阈值范围外的云像元判为陆地。

使用调整后的算法对MODIS积雪数据进行去云处理，得到逐日无云积雪覆盖数据集。利用气象站点的实测雪深数据对该数据集进行精度评估。基于气象站点的位置提取影像像素值，与站点实测值成对作为验证样本数据，分别对去云前的MODIS原始逐日积雪覆盖数据和去云后的逐日无云积雪覆盖数据进行精度验证。

验证结果如表 2所示，在晴空条件下(排除被云覆盖的站点)原始MOD10A1影像的积雪分类精度为83.26%，陆地分类精度为98.56%，整体分类精度为90.91%；移除云阻碍的影响后(即所有站点)逐日无云积雪覆盖数据的积雪分类精度为83.69%，陆地分类精度为94.41%，整体分类精度为89.05%。以MOD10A1影像为基准，云像元被分类为积雪和陆地准确性结果如表 3所示，整体精度为89.39%，积雪分类精度为83.27%，陆地分类精度为95.51%。结合两表数据分析，可以看出云像元被重分类为积雪的精度与原始数据的分类精度基本一致，陆地分类精度损失约3%。去云后的无云积雪覆盖数据集的精度与原始数据在晴空条件下的精度相比，积雪分类精度轻微提升，陆地分类稍有损失，但增加了大量的可用样本。可见该去云算法可以有效地减少云像元的影响，并保证精度。无云积雪覆盖数据集弥补了数据空间不连续的不足，提供更多的可用信息，更加适用于后续的积雪研究。

3 研究方法 3.1 积雪参数提取 3.1.1 积雪日数

积雪日数(snow-cover days, SCD)^{[21, 27]}表示一个像元在指定时间序列上观测到积雪的次数，数学描述为

${\rm{SCD}} = \sum\limits_{i = 1}^n {({s_i})} , $

(1)

式中：SCD为积雪日数，取值范围为0~366；n为总观测数，本文中为一整年的天数(当年9月1日至次年8月31日，365或366)；s_i表示逐日无云积雪数据集，将该数据集重编码为1代表积雪，0代表非雪。

3.1.2 积雪面积

积雪面积(snow cover area, SCA)^{[5, 28]}表示指定区域内的积雪覆盖范围，可以直观描述研究区的积雪覆盖量，通过时间序列的分析，简单明了地掌握研究区积雪的年际变化。SCA是通过统计逐日影像中积雪像元的个数，乘以栅格像元的大小来表示当日的积雪面积，基于不同时间尺度的需求，求取均值用以后续的研究。

3.2 趋势分析

采用一元线性回归最小二乘拟合模型模拟每个像元j的积雪日数的变化趋势^{[2, 5]}，数学描述为

$\begin{array}{l} {\rm{Slop}}{{\rm{e}}_j} = \frac{{n\sum\limits_{i = 1}^n {i \cdot {\rm{SC}}{{\rm{D}}_{ij}}} - \sum\limits_{i = 1}^n i \cdot \sum\limits_{i = 1}^n {{\rm{SC}}{{\rm{D}}_{ij}}} }}{{n \cdot \sum\limits_{i = 1}^n {{i^2}} - {{(\sum\limits_{i = 1}^n i )}^2}}}, \end{array} $

(2)

式中: Slope_j表示第j个像元在研究期间积雪日数的变化倾向率；i表示年序号；n表示研究时间长度，本文中为14 a。若Slope_j>0表示该像元在研究期间变化趋势是增加的，反之是减少的。并通过F检验判断其变化趋势的显著性水平。

3.3 相关性分析

基于相关性分析探究积雪特征参数与温度、降水的关系^[2]，通过计算各变量间的相关系数以及显著性。数学描述为

$\begin{array}{l} R = \frac{{\sum\limits_{i = 1}^{14} {({X_i} - \bar X)({Y_i} - \bar Y)} }}{{\sqrt {\sum\limits_{i = 1}^{14} {{{({X_i} - \bar X)}^2}{{({Y_i} - \bar Y)}^2}} } }}, \end{array} $

(3)

式中：X、Y分别表示相关分析的变量，X、Y分别表示2个变量的均值。相关系数R值的范围为[-1, 1]，绝对值越大，相关性越强，反之越弱。

4 结果与讨论 4.1 SCD变化分析

统计黑龙江流域在2002—2016年间的逐年SCD，得到14 a平均SCD的空间分布如图 3(a)所示。可以看出，黑龙江流域积雪空间分布极不均匀。整体来看，SCD与纬度的相关性较好，随纬度增加呈逐渐增大的趋势。且海拔相对较高的地区的SCD也较大。统计表明，黑龙江流域14 a的平均SCD范围为0 ~ 270 d，主要集中在90~150 d，年平均SCD为123 d。其中小于90 d的区域主要在流域中南部地区，占比13%，大于180 d的比例仅占5%，主要分布在海拔较高的山区，如流域西部边缘的肯特山脉区域，东北部的雅布洛夫山脉、斯塔洛夫山脉等。另外，在森林密度较大的地区如大兴安岭地区，由于光学遥感手段自身对积雪的低估，加之不同的森林密度对积雪的积累、消融以及再分布等过程有着复杂的影响，导致森林密度较大的地区的SCD比周围纬度和海拔相似的非森林区较短，该现象与Wang等^[6]研究结果一致。

利用式(2)得出研究区14 a SCD的变化情况以及显著性分布如图 3(b)所示，经过统计发现，约45%的区域SCD呈减少的趋势，分布在流域的西部、北部以及东部的边缘地区，南部也有小部分地区出现减少的趋势。约1.5%的区域呈显著减少(P < 0.05), 主要分布在海拔较高的山区。流域的中南部并向东西两侧延伸约54%的区域SCD呈增加趋势，仅2%的区域通过0.05显著性水平检验。整体上来看，SCD呈增加趋势的占比稍大于减少趋势的比例。

4.2 SCA变化分析 4.2.1 SCA年际变化

统计14 a的年均SCA如图 4(a)所示，可以看出年际间有较大的波动性，整体呈微弱的增加趋势，这与全球积雪逐渐减少的趋势相反，2007、2012年为特殊年份。黑龙江流域2002—2016年间平均温度和年降水量分布如图 4(b)所示，可以看出其年际差异较大，均表现为增加趋势，这表明在研究期间黑龙江流域的气候在变暖、变湿。综合三者来看，2007年SCA最小，平均温度最高，降水量较低。而2012年SCA最大，平均温度最低，降水量最高，表明SCA的变化是气候变化的具体表现。

为进一步探索黑龙江流域的SCA变化，本研究按季度(9—11月为秋季，12—2月为冬季，3—5月为春季，6—8月为夏季)进行分组讨论(图 5)。与年SCA变化相比，月SCA的异常变化更加明显。从图 5可以看出研究期间月SCA变化波动最大的为春季的3月和秋季的11月，其次是冬季的12月和2月以及春季的4月。黑龙江流域的积雪在11月份大范围积累，年际间SCA呈减少趋势，且冬季整体上表现为增加趋势，表明部分地区可能有积雪积累期后推的迹象。冬季黑龙江流域几乎整个都被积雪覆盖，年内均值SCA最大在1月，根据趋势线可知，12月和1月的SCA正向增加，2月的SCA却反向增加。这意味着2月的SCA有衰退迹象。另外，对比冬季2月和春季3、4月的SCA，可以看出，SCA在3、4月份骤减，说明黑龙江流域融雪期主要集中在3、4月份，但3、4月份的年际线性变化趋势不甚明显。黑龙江流域夏季积雪覆盖区域较小，但各月份储存的SCA年际间波动性较大。夏季SCA整体上呈增加趋势，反映了夏季气候的异常变化。

4.2.2 SCA与气候因子的关系

气温和降水是影响积雪面积的重要因素，也是直接影响因素^[29]。研究不同时段的气候因子与SCA的关系可以更好地揭示气候因子对积雪的影响。结合气象资料，分析研究区14 a的月平均SCA及季平均SCA与同期气象因子之间的相关关系如图 6和表 4所示。结果表明在黑龙江流域气温和降水对积雪面积的影响是复杂的，不同时间段的SCA与气象因子相关性有所差异，不同温度指标的影响程度也不尽相同。其中春季与温度呈明显负相关(P < 0.01)，与降水相关性较差，在3、4月份表现尤为明显。夏季SCA变化主要受降水的影响(P < 0.01)，主要表现在6月份。秋季的SCA与气候因子的相关性十分复杂，9月仅与最高温度有较好的负相关关系(P < 0.05)，10月与最高、最低气温和降水都有较好的相关性(P < 0.05)，与平均温度无明显相关性。整个季节来看，秋季仅与最高温度有较好的负相关关系(P < 0.05)。冬季12、2月与降水有明显正相关关系(P < 0.01)。1月、5月、7月和11月与气候因子各指标的相关性均未通过0.05显著性水平检验。

为了更加全面地探究气候因子对SCA的影响，对季节性的SCA以及同期气候因子的年际变化做了分析。图 7表明14年间温度和降水在不同季节的变化趋势不同，平均温度在春、夏、秋季略有升高，在冬季轻微降低，降水的变化在各个季节均表现为略微增加。综合前文相关性的分析，可以看出黑龙江流域2002—2016年14年间春季SCA变化主要受温度的影响，对比温度和SCA的变化曲线，并没有发现较为明显的负相关关系。冬、夏季的SCA与降水的年际变化表现出一致的增加趋势。

尽管得到了上述研究结果，但也存在一些不确定性，本研究选择的时间序列仅14 a，各因子的线性变化趋势都不十分显著。比如春季，SCA与温度呈明显的负相关关系，但线性变化趋势没有表现出这一特征，这可能就是由于研究的时间序列过短造成。另一方面是黑龙江流域地域面积较大，拥有多种小气候，常有流域北部涌入寒流而南部地区仍处于酷暑的情况，不同地区间的气候差异可能造成某些流域均值失去代表性。此外，积雪变化是多种气候要素共同作用的结果^[30]，比如秋季的SCA与降水、温度都没有十分显著的相关关系，想要探索其面积减少的原因，还需进一步增加其他气候因子进行讨论。后续将补充更多研究年份，针对不同气候的区域进行分区分析，降低流域过大的影响，并引入更多的气候因子开展深入研究。

5 结论

本研究利用去云后的MODIS逐日积雪覆盖产品MOD10A1/MYD10A1数据，开展黑龙江流域2002—2016年间的积雪时空变化及其与气候因子之间的相关关系的研究，得到以下结论：

1) 黑龙江流域的平均SCD主要集中在90~150 d，14 a平均SCD为123 d。空间分布上受地理环境影响较大，年均SCD的时空变化表现为增加的比例稍大于减少的比例。其中增加的区域主要集中在流域的中南部。

2) 积雪面积(SCA)能够直观的表达流域积雪在年际间的变化。14年来，黑龙江流域年平均SCA呈微弱的增加趋势，这与全球的积雪变化趋势不同。季节上表现为冬、夏季SCA呈微弱增加，秋季SCA有轻微的减少趋势，春季无明显的线性变化趋势。

3) 通过对SCA与降水、温度进行相关性分析可知，气温和降水与积雪的关系密切且复杂。整体上表现为，春季SCA与不同的温度指标均有明显负相关关系(P < 0.01)，夏季、冬季则与降水有明显的正相关关系(P < 0.01)，而秋季仅与最高气温有较好的负相关关系(P < 0.05)。