积雪作为地球上重要的水资源之一，融化的冰川积雪是许多淡水湖的源头（张佳华等，2008）。南北半球的中高纬度地区有着大量的积雪和冰川，终年不化的积雪和季节性积雪广泛存在于高海拔地区（孙燕华等，2014）。相比于终年积雪，季节性积雪对天气和气候的影响更为显著，季节性积雪对地热状况有着重要的反馈作用，积雪开始和结束时间、积雪覆盖范围与深度都会对区域内的土壤含水量造成影响（杨林等，2019），进而对地气系统的辐射平衡、水文循环产生显著的影响（唐志光等，2017），中国黄河中下游和长江中下游的旱涝灾害又受到二者的影响。全世界98%的季节性积雪位于北半球，北半球每年12月到次年1月的最大积雪面积占整个陆地面积的45%以上，欧亚大陆作为北半球的主要积雪区，其积雪面积更是达到北半球积雪面积的60%以上（Lemake，2007; 杨林等，2019）。在全球变暖的气候变化背景下，北半球的积雪面积正逐年下降，中国的积雪面积也处在缓慢下降的趋势中（吴彦，2005）。大量的积雪融化，造成淡水湖泊和河流水面上升，空气中水分充足，降水降雪频发，易造成暴雪等极端天气气候事件，对中国的农业生产和自然生态平衡造成严重的后果，引发一系列问题。因此掌握中国的积雪时空动态变化是十分重要的。

对于大范围的积雪时空变化，国内外近几十年已经开展了大量的研究工作。当前比较主流的获取积雪信息的方法有气象站点监测和卫星遥感数据获取。Huang et al.（2019）使用气象站点数据，得出1952—2012年中国年积雪持续时间在减少的结论。郝祥云等（2018）通过分析锡林河流域周围1981—2016年的逐日气象数据，发现近35 a内雪深和积雪日数呈上升趋势，且两者显著相关。气象站点数据虽然准确度较高，方便记录，但对高海拔和气候恶劣的地区难以覆盖，卫星遥感监测解决了此类难题。当前卫星遥感对积雪的监测分为主动微波遥感和被动微波遥感。主动微波遥感借助积雪与其他覆盖物在不同波段的反射系数不同来判别，被动微波遥感可以全天候穿过云层，提取地表的亮温，通过雪深反演公式得到积雪深度。基于MODIS传感器得到的MOD10A1、MYD10A1积雪产品数据集，其空间分辨率500 m，时间分辨率1 d，Huang et al.（2016）利用此数据集对中国积雪面积、积雪天数进行了探讨，结果表明短暂积雪分布区为中国的东部和东南部以及新疆和内蒙古的部分地区，内蒙古、新疆北部和青藏高原为稳定积雪区，新疆西天山和青藏高原山区是永久性积雪的主要地区。众多研究表明，中国积雪面积在近几十年内逐步减少，并且减少的速度也在加快（Déry and Brown，2007; Brown and Robinson，2011）。中国积雪分布具有明显的空间差异性，随着纬度的升高，积雪逐步增多，西部积雪明显多于东部; 山地积雪多于盆地和平原，且山地积雪随着海拔的升高而增多（李栋梁和王春学，2011）。青藏高原作为世界第三极（除多等，2017; 姚檀栋等，2017），地形复杂，山脉、峡谷、盆地纵横交错，近20 a川西高原（青藏高原东南缘和横断山脉一部分）的积雪覆盖率整体呈缓慢减少的趋势，积雪覆盖时间的空间分布差异明显，随海拔的增高而增加，其他地区变化不显著（钟鼎杰等，2021）。高原积雪由于其海拔高，持续时间长等特点，持续作用于大气环流，从而影响中国夏季降水量（丁峰等，2009）。

总体上，基于国产风云三号气象卫星对中国区域近10 a积雪状况的研究较少。为了更加详细的对中国区域内近10 a的积雪变化进行分析，本文采用FY-3气象卫星的积雪覆盖产品数据，对2010—2019年中国积雪覆盖频率（Snow Cover Frequency，SCF）进行逐月演变分析和线性回归提取长期趋势，提出基于经验模态分解（Empirical Mode Decomposition，EMD）（Huang et al.，1998）的EMDTrend指标来表现积雪覆盖率（Snow Cover Rate，SCR）的长期趋势，使用Mann-Kendall检验（符淙斌和王强，1992; Wu et al.，2011）方法来衡量SCR的变化情况。

1 数据和方法

1.1 研究区

中国国土幅员辽阔，南北跨越纬度约50°，东西跨越经度约62°，占欧亚大陆面积的17.5%，自南向北包含热带、亚热带、暖温带、中温带、寒温带5个温度带，以及处于高原气候带的青藏高原周边地区（郑景云等，2010）。中国地势东低西高，拥有高原、山岭、平原、丘陵、盆地等地形以及昆仑山脉、横断山脉等大小山脉（方如康，1995），复杂的地形、走向多样的山脉以及沿海与内陆城市的差别形成了多种多样的气候，降雪的季节性特征显著，每年冬季为降雪频发时期，而冬季降雪在地理位置上尤为不同，积雪分布具有明显的空间异质性，其中新疆、东北、青藏高原是积雪分布的重要区域（Huang et al.，2016）。中国两大河流黄河、长江分别起源于青藏高原巴颜喀拉山脉和唐古拉山，青藏高原的冰川积雪融化状况对黄河、长江中下游乃至中国大部分河流流域影响显著。因此本研究选择中国陆地领土为研究区（图1）。

1.2 研究数据

FY-3气象卫星是中国继FY-1的第二代极轨气象卫星，携带了中分辨率成像光谱仪、地球辐射探测仪、扫描辐射计、微波温度计等11台设备，可以做到大气三维探测，实现全球高分辨观测和全天候全天时工作（谷松岩等，2022），为天气预报、自然灾害和生态环境监测、军事服务等提供有力帮助，搭载于FY-3上的中分辨率成像光谱仪相比于MODIS可活动角度更大，扫描范围更广（范天锡，2002）。通过FY气象卫星得到的遥感数据已发展出覆盖大气、陆表、海洋等多类型产品，并且应用到数值模拟和预测等各种业务模式中（杨军等，2018）。本研究使用的数据来自风云卫星遥感数据网（http://satellite.nsmc.org.cn）提供的FY-3的多仪器融合数据（MULSS）产生的积雪覆盖产品（SNC），该产品基于可见光红外扫描辐射计和中分辨率成像光谱仪融合而成，通过归一化积雪指数对云和雪判别，对像元逐个进行分类（杨军和董超华，2011）。SNC产品代码共有11个分类，定义如表1所示，分类实例如图1。

2010年9月—2016年5月的数据使用FY-3A，2016年9月—2020年1月的数据使用FY-3C。该数据空间分辨率为1 km，时间分辨率为10 d，定义每年9月至次年5月为一个研究年份。分别选取30张FY-3A、15张FY-3C的栅格数据，共计约7.2×10⁸个像元与作为真值参考的中国雪深长时间序列数据集（SD）比较（Che et al.，2008）。如表2所示，以五个指标对FY-3的积雪产品进行精度评估，FY-3A和FY-3C的积雪判别准确率分别为83.20%、85.92%。张帅等（2018）以MOD10C1（MYD10C1）全球日积雪覆盖数据集作为参考，对2010—2014年风云三号MULSS积雪产品从总精度、偏差等多个方面做了评估，得出两者有雪期积雪分布趋势一致，在空间分布和时间演变上保持了较好的一致性。闵文彬等（2021）以2018年10月—2019年4月青藏高原的气象站点积雪观测数据做真值，得出MULSS积雪产品判别准确率为87.18%，误判率为12.81%。综上所述，FY-3的MULSS积雪产品质量稳定并拥有较高的精度，满足应用需求。

积雪覆盖频率（Snow Cover Frequency，SCF）为每个像元在一段时间内有积雪覆盖的天数除以总天数所得到的结果，反映积雪覆盖持续时间的总体特征（秦艳等，2018），具体定义如下:

其中:F _SCF-k为n年内累计某月像元k的积雪覆盖频率; m为该月的天数; f_snow（k）=1表示像元被积雪覆盖，f_snow（k）=0表示像元未被积雪覆盖（钟鼎杰等，2021）。

积雪覆盖率（Snow Cover Rate，SCR）为积雪覆盖的像元面积除以研究区总面积所得的结果，具体定义如下:

注:将SNC和SD同时有雪的像元数记为a; SNC和SD同时无雪的像元数记为b; SNC有雪而SD无雪的像元数记为c; SNC无雪而SD有雪的像元数记为d; 准确率为: $\frac{a}{a+c}\times 100\mathrm{\%}$; 召回率为: $\frac{a}{a+d}\times 100\mathrm{\%}$; 误判率为: $\frac{c}{a+c}\times 100\mathrm{\%}$; 漏判率为: $\frac{d}{a+d}\times 100\mathrm{\%}$; 总体精度为: $\frac{a+b}{a+b+c+d}\times 100\mathrm{\%}$.

其中:P_i（1）为研究区内第i个为积雪的像元; n为总像元数。

1.3 方法

1.3.1 经验模态分解

积雪覆盖具有季节性特征和明显的周期，这种拥有周期性或季节性的长时间序列是非平稳的，经验模态分解（Empirical Mode Decomposition，EMD）是一种处理非平稳数据的方法（Huang et al.，1998）。EMD在处理线性和非线性数据时均有良好的效果，相比于小波分析和傅立叶变换，EMD基于数据自身驱动，无需基函数，每次分解根据约束条件（极值点的数量和过零点的数量之差小于等于1; 上下包络线的平均值为0）判断h（t）是否为一个内涵模态分量（Intrinsic Mode Functions，IMF）。EMD分解步骤如下:

1）对原始序列去趋势和标准化处理，得到处理后的数据x（t）。

2）找到数据x（t）所有的上下极值点，用三次样条插值分别拟合上下极值点，得到极大值包络线和极小值包络线。

3）求出上下包络线的平均值m（t），然后令:

判断h（t）是否满足约束条件，如果不满足则重复步骤2）和3）。

4）每得到一阶IMF就从原序列中减去它，直到最后剩余部分r_n是单调序列或者常值序列，此时r_n为周期很长的波，因此可以将此看作是数据的长期趋势（Wu et al.，2011; 刘权和姚凤梅，2021）。

这样经过EMD就将原始序列x（t）分解成N阶IMF（用F_i表示）和一个代表趋势的序列:

1.3.2 基于EMD的趋势提取

EMD分解的r_n可以被看作数据随时间变化的趋势，但经过标准化处理和EMD分解的序列无单位、无物理意义。在量化不同研究区的数据时，无法对分解后的数据类比分析。为了解决这个问题，提出一个变量E _Trend（t）。该变量的求解过程如下:

该变量代表了数据经分解后，以初始日期数据为基准随时间的增量，可以表明数据的长期变化趋势。在对不同研究区对比时，以百分比的形式展现趋势变化。

1.3.3 Mann-Kendall检验法

Mann-Kendall检验是一种非参数统计方法，常用于气象方面的数据变化检验，其特点是不要求数据服从正态分布、不受少数异常值和缺失值的影响（魏凤英，2007）。

对于具有n个样本量的时间序列$x_1，x_2，\cdots ，x_n$，构造一个秩序列:

其中:秩序列S_k是第i时刻数值大于j时刻数值个数的累计数。在时间序列随机独立的假定下，计算累计数的期望值E（s_k），以及方差var（s_k），定义统计量:

其中:U _F-k为标准正态分布，它是按时间序列x顺序$x_1，x_2，\cdots ，x_n$计算出的统计量序列，然后按时间序列x的逆序$x_n，x_{n-1}，\cdots ，x_1$计算出逆序列$U_{\mathrm{F}-k}^{\mathrm{\text{'}}}$。

若U _F-k的值大于0，则表明原样本呈上升趋势，小于0则呈下降趋势。当它们超出临界线时，表明上升或下降趋势显著。若U _B-k与U _F-k两条曲线出现交叉点，且交叉在临界线之间，它们交叉点对应的时刻便是突变开始的时间。

2 结果和讨论

2.1 积雪覆盖频率（SCF）

研究表明，中国积雪每年随着冬季的来临和结束呈现“增加-减少”的年际波动（Huang et al.，2016）（下文中提到的年份均指每年9月至次年5月，如2016年指的是2016年9月至2017年5月），图2给出了2010—2019年中国SCF在秋、冬、春季的变化情况。

整体来看，典型积雪覆盖区有东北、新疆部分地区、青藏高原，SCF在空间分布上差异明显（覃郑婕等，2017）。新疆北部和东北地区的SCF季节性差异显著，每年九月该地区积雪偏少，随着冬季的来临，准噶尔盆地和大兴安岭北部SCF显著增加。东北地区、新疆北部以及内蒙古东北部二月份的SCF可以达到80%，华北平原和黄河中下游流域SCF达到40%。进入春季后，黄河流域以南的地区积雪率先消融，SCF减少至0，东北和新疆北部四月底SCF逐渐减少到20%以下。青藏高原除昆仑山脉外也呈现出单峰的变化情况，昆仑山脉平均海拔为5 500~6 000 m，海拔5 600~5 900 m处为雪线（韩芳等，2014），终年不化的积雪覆盖在雪线之上，9月时SCF达到20%。青藏高原的月最大SCF相比其他地区滞后1个月，昆仑山脉西部的SCF在四月达到80%。五月份，中国大部分区域SCF为0，青藏高原仍有大量积雪，SCF平均为40%。青藏高原大部分地区的夏季平均温度在10℃以下，高原腹地的年平均气温不足0℃，降雪落到地面，形成积雪和冰川，气温未达到积雪消融的程度（Leathers and Luff，1997; 李海生等，2017），海拔较高的区域冰川终年不化（Schmidt et al.，2009; Zhou et al.，2017）。尽管如此，青藏高原的积雪相比于世界上其他的冰原地区融化得更快，过去的一百年里，青藏高原的升温速度是全世界的2倍（Monitoring and Programme，2011），升高了近1.5℃，短期内的积雪消融会增加下游江河水量，但随着积雪冰川的减少，部分径流甚至会消失（王顺久，2017）。

逐月SCF变化过程能较好展现多年平均SCF季节性差异，但无法表明长时间SCF趋势在空间上的变化。因此对每个像元使用基于最小二乘法的一元线性回归，对近10 a的中国SCF空间分布进行趋势分析，并进行显著性检验。图3显示了通过置信度为95%显著性检验的2010年9月至2020年1月中国SCF趋势分布。

中国SCF变化比较复杂，具有显著的空间异质性（陈海山等，2019），近10 a的SCF减少出现于大兴安岭、东北平原、内蒙古中东部、华北平原以及长江中下游平原，西部塔里木盆地周边以及青藏高原西南部也处在减少的状态（Zhao et al.，2007; 唐志光等，2017），大部分地区SCF每年减少0.1%至0.2%，内蒙古中东部每年减少0.35%。大兴安岭和长白山脉附近处于显著减少的现状，每年减少0.5%。SCF增加的地区偏少，黄土高原北部、祁连山脉每年增加0.2%，位于高原山地气候的柴达木盆地周边、昆仑山脉、巴颜喀拉山脉以及唐古拉山脉，SCF每年增加0.3%。显著增加的地区有小兴安岭、柴达木盆地与祁连山脉交接的区域，每年增加0.4%，这与Huang et al.（2016）得出的中国东北部小兴安岭的积雪天数显著增加、青藏高原减少的研究结论十分吻合。总体上看，中国SCF趋势由东向西处于“减少-增加-减少”的现状，增加的为盆地、山脉此类地形复杂的区域，大部分区域处于减少的状态（陈文倩等，2018），尤其我国东北部减少的非常明显。东北是中国纬度最高的地区，SCF却逐年减少，具体原因后续需结合气温、降水等因素进一步研究。

2.2 积雪覆盖率（SCR）

逐月SCF和SCF趋势从季节性和空间差异两个方面说明了近10 a中国SCF的演变情况。为了更深层次的对中国领土内积雪面积的变化进行分析，由2010年9月至2020年1月共86个月的风云SNC产品，得出每月平均SCR。鉴于对SCF分析时发现部分地区降雪丰富，如小兴安岭、东北平原、青藏高原、祁连山脉等，因此选取了青海、内蒙古、新疆、西藏4个省（区）及东北地区做对比研究。对中国及中国典型积雪覆盖区的SCR进行线性回归前，使用差分法去除SCR的季节性周期，由表3的SCR线性趋势可以看到只有新疆的SCR趋势通过了显著性检验（P ＜0.05），每年减少0.03%。由此看出线性回归对周期性数据的拟合效果较差，因此使用EMD提取SCR趋势，由图4分解得到的EMDTrend趋势可以看到，内蒙古、西藏与中国整体的趋势均处于缓慢减少的过程，内蒙古和西藏近10 a SCR下降了约1.5%，中国下降了约1.2%。青海的SCR在EMDTrend趋势中有一个“减少-增加”的波动，SCR至2016年下降了0.9%，随后开始增加，与2010年相比，最终下降了0.6%。经过EMD分解的数据最后得到的是常值序列或单调序列，青海地区的SCR在提取四阶IMF后，不满足上下包络线均值为0的条件，因此选取第五阶IMF作为青海的EMDTrend趋势。新疆和东北的EMDTrend趋势中SCR下降缓慢，10 a分别减少了0.18%、0.1%，东北地区2010—2012年的SCR有微小的增加。东北地区SCF下降显著，SCR下降程度微小，但SCF和SCR有没有必然的联系并不能确定。SCR表明的是积雪覆盖面积的变化情况，SCF受地表温度和地形多种因素影响，两者之间的相关性有待进一步考证。对东北近10 a的SCR数据线性回归时，剔除周期性的影响后效果仍不好（P =0.66），东北处于高纬度地区，每年冬季降雪非常丰富，SCR有显著的季节性特征，EMDTrend趋势相比于线性趋势能更好地表示SCR长期变化的情况。

EMDTrend代表了SCR的长期变化趋势，为了进一步探究SCR在10 a内的详细变化情况，对SCR进行Mann-Kendall检验，结果如图5所示。中国、新疆以及西藏的SCR总体趋势都是先上升后下降，中国2010年到2015年期间SCR呈上升的趋势，随后一直下降。新疆在2017年10月由上升转为下降，西藏则在2016年3月开始下降。东北、内蒙古和青海处于“上升-下降-上升”的波动中。东北地区在2011年10月至2012年10月以及2016年3月至2020年有两个下降的时期，在2018年底至2020年1月SCR显著下降。内蒙古在2010年10月至2014年1月SCR处于上升的趋势，随后开始下降，但在2015年11月至2016年6月有一个短暂的上升期，2018年10月至2020年1月显著下降。青海省在2011年前有一个短暂的下降期，随后开始上升，2015年3月继续下降，2017年以及2018年SCR显著下降。综上所述，中国典型积雪覆盖区SCR总体上均是先上升后下降的变化情况，发生转变的年份均在2016年前后。

3 结论与讨论

本研究基于FY-3气象卫星的MULSS积雪产品数据集，在对原始数据集定位、拼接、裁剪、提取积雪信息的基础上，从SCF和SCR两个方面对中国近10 a的积雪时空动态进行了比较深入的研究。对SCF从季节性演变和长时间趋势两个方面在空间上的特征进行了分析。在SCR方面，本研究提出了一个基于EMD的EMDTrend变量来衡量SCR的长期趋势，对积雪显著的五个省份及地区使用Mann-Kendall检验，来探究SCR的详细变化状况。

1）中国的SCF具有显著的季节性特征，每年9月至次年5月，呈现先增加后减少的特点，这与大部分学者研究成果一致。新疆北部、东北地区、内蒙古中东部、青藏高原在每年2、3月份达到当年最大SCF。春季开始，黄河流域积雪率先开始消融，除青藏高原外的其他地区在5月时，积雪基本消融，青藏高原的积雪因气温、地形等因素仍有保留。由东向西，我国SCF趋势处于“减少-增加-减少”的现状，增加的多为盆地、山脉类地形复杂的区域（Zhao et al.，2007），如小兴安岭、柴达木盆地与祁连山脉交接的区域。东北地区、内蒙古北部、华北平原以及长江中下游流域SCF处于减少的状态，尤其我国东北部减少得非常明显。希爽和张志富（2013）也指出东北地区的积雪天数在减少，造成东北SCF减少的原因，待基于FY系列精度更高的再处理资料发布后结合气温、降水等多种因素进一步分析。

2）对典型积雪覆盖区的SCR研究表明，内蒙古、西藏的SCR有长期下降的趋势，近10 a下降了1.5%，下降的趋势逐渐变缓。青海在经历6 a的下降后，有一个增加的趋势，至2020年最终下降了0.6%。新疆和东北的SCR分别减少了0.18%、0.1%，长期趋势比较稳定。中国及中国典型积雪覆盖区的SCR在2016—2017年均发生由增加到减少转变，造成变化的原因可能是副热带高压异常，影响了高空的水汽输送，对气温和降雪的变化造成了一定作用，最终导致SCR改变（余志豪和李海盛，2000; 张人禾等，2016）。全球SCR减少已成为大多数学者的共识（希爽和张志富，2013; 孙燕华等，2014），影响SCR减少的因素包括气温降水、地形、海拔、植被等（杨林等，2019），史建桥（2014）指出风速也会影响积雪面积变化，这对后续研究工作有很大的启示。

本研究分析了中国积雪的时空变化动态，对掌握中国典型积雪覆盖区的时空动态变化有积极作用，后续将结合气象数据对造成SCF和SCR显著变化的原因进行探讨分析。

致谢:中国气象局国家卫星气象中心提供了FY-3气象卫星的积雪产品在线下载服务。

参考文献