积雪是冰冻圈的重要组成部分，是气候系统中比较活跃的因子，对气候环境变化较为敏感，其年内和年际变化影响全球气候和水文的变化^［1-4］。积雪还是冰川积累的物质基础，高山冰雪融水是内陆河重要的补给来源^［5］。中国西北地区属于干旱半干旱地区，积雪融水对当地农牧业灌溉有重要意义，积雪面积的变化会直接影响当地的环境演变及人类生存。然而，随着每天太阳辐射能量的不断变化，积雪可能会发生剧烈变化。因此，快速、实时并精确地进行积雪面积变化监测对于气候演变模拟、水资源利用管理、灾害分析评估具有重要意义^［6-8］。

近年来，随着遥感理论和技术日新月异的发展，国内外学者针对一定区域内长时间序列的积雪面积变化已经完成了大量研究^{［10，13-15］}。利用AVHRR传感器获得的遥感数据资料对加拿大每日积雪深度数据研究，发现1955~1989年加拿大积雪面积持续减少，并发现面积减少的原因与季节有关系^［11］；利用AVHRR传感器获取的数据对欧洲阿尔卑斯山脉积雪时空变化研究，发现27 a间阿尔卑斯山脉高山地区降雪量的下降趋势有所减缓，并提出了一套新的无云卫星积雪产品^［12］；使用MOD10A2数据对2001~ 2017年祁连山区的积雪面积变化进行研究，发现积雪面积年际波动趋势较大，总体呈减小趋势，多年平均积雪面积约为5×10⁴ km²，占祁连山总面积的25.9%^［7］；使用FY-3A/MERSI数据建立NDSI积雪指数模型，并对2010~2011年祁连山区积雪进行研究，发现模型能够客观精细地反映祁连山区积雪的空间分布与动态发展过程^［6］。

随着中国风云系列气象卫星的蓬勃发展及数据量的剧增^［16］，有必要结合FY-3C MULSS（多仪器融合数据）积雪覆盖日产品较高空间分辨率的特点和FY-4A AGRI（多通道扫描成像辐射计）数据高时间分辨率的特点，对祁连山区积雪面积变化和时空演变进行研究。一方面有利于发挥我国风云系列气象卫星数据的优势，另一方面弥补了目前大多数研究中仅用单一数据研究积雪面积变化的不足。

在以往学者对祁连山区积雪面积和时空分布特征研究的基础上，提出了一套适合于FY-4A AGRI传感器数据的动态多阈值多时相云隙间积雪判识方法，以及基于像元级和决策级的FY-3C MULSS积雪日产品和FY-4A积雪判识结果的融合方法。利用Landsat 8 OLI数据和MOD10A2积雪产品对方法判识结果的精度进行验证，并对祁连山区2018年3月~2019年3月的积雪面积及其时空演变展开研究，以期为政府部门决策提供参考指标。

祁连山区（94°10'~103°04' E，35°50'~39°19' N）地处甘肃、青海两省交界处，东起乌鞘岭的松山，西到当金山口，北临河西走廊，南靠柴达木盆地^［8］。山峰海拔多为4 000~5 500 m，最高峰疏勒南山的团结峰海拔5 808 m。祁连山区地处北温带，受海拔及地理位置的影响，年平均气温低于4 ℃。降水受海拔、地理位置、坡度和坡向等因素的影响，季节和年际变化较大，降水量总体表现为从东向西递减的变化趋势。祁连山区是永久积雪区、稳定季节积雪区、年周期性不稳定季节积雪区、非年周期性不稳定季节积雪区的综合积雪区域^［7］。且祁连山区具有特殊的地理位置和地形，下垫面类型十分丰富，在我国地理位置和自然环境中都占有重要的地位。

采用的数据如表1所示，包括FY-3C MULSS（多仪器融合数据）积雪覆盖日产品数据、FY-4A AGRI（多通道扫描成像辐射计）成像仪全圆盘2KM L1数据、数字高程模型数据（SRTM）、MOD10A2积雪覆盖产品数据和Landsat 8 OLI卫星遥感数据。其中，FY-4A AGRI波段特征及其主要用途如表2所示。

（1）FY-3C MULSS积雪覆盖日产品的预处理。FY-3C MULSS积雪覆盖日产品包含地理位置信息，通过HDFView-3.0.0软件读出影像的位置信息和数据编码信息，覆盖研究区的两景数据的位置信息分别为3090（坐标系：WGS84左上角地理坐标：90°E、40°N）和30A0（坐标系：WGS84左上角地理坐标：100°E、40°N）。利用ENVI5.3完成影像投影、几何校正、拼接和裁剪等预处理工作。

（2）FY-4A AGRI成像仪全圆盘2KM L1数据的预处理。FY-4A AGRI成像仪全圆盘2KM L1数据是具有7个波段的一级原始数据。首先，使用查找表获取了可见光近红外波段的反射率。其次，使用官方提供的经纬度查找表建立研究区的GLT文件完成了影像的几何校正。最后，利用余弦校正法进行了太阳高度角校正的工作。

（3）MOD10A2积雪覆盖产品的预处理。借助EOS/MODIS 数据中心提供的软件MODIS Reprojection Tool（MRT）对MOD10A2 数据进行格式转换、镶嵌、投影变换等3个步骤的处理。

（4）Landsat 8 OLI数据的预处理。使用ENVI 5.3对Landsat 8 OLI数据进行了波段融合、影像镶嵌和影像裁剪的预处理工作，由于Landsat 8 OLI数据预处理方法发展成熟，所以此处不再赘述相关过程。

如图2所示，将FY-3C号和FY-4A号数据的融合方法分为两部分展开研究。第一部分，构建适合于FY-4A号卫星数据的动态多阈值多时相云隙间积雪识别方法。首先，对每一景FY-4A数据使用多阈值积雪识别方法获得每一景的积雪识别结果。其次，FY-4A数据具有高时间分辨率的优势。因此，在云流动的间隙中判识积雪像元，可以大大降低云对光学影像判识积雪的影响程度。即每日可利用的10~12景数据分别使用积雪识别方法，判识出每天不同时间的积雪覆盖情况，再以大量数据试验后的概率值为基准，融合每日的10~12景影像并得到除去云后的FY-4A日积雪覆盖数据。

第二部分，将FY-3C MULSS积雪覆盖日产品和FY-4A数据的积雪判识结果按一定方法融合以提高判识精度，从而达到综合判识祁连山区积雪覆盖情况的目的。融合后的数据将同时具有FY-3C较高空间分辨率和FY-4A高时间分辨率的优势。最后利用Landsat 8 OLI数据和MOD10A2数据对融合后的数据进行精度评价并提取祁连山区积雪面积。

通过分析地物反射率波谱曲线，发现仅仅依靠单一波段的反射率特征无法很好地区分积雪和其他地物。为了建立适合于FY-4A AGRI数据的积雪识别方法，引入目前光学遥感常用的归一化差分积雪指数（NDSI）进行积雪识别^［17-20］，具体计算公式为：

其中：Green和SWIR1分别代表FY-4A号卫星数据的第二波段（0.55~0.75 μm）和第五波段（1.58~1.64 μm）。然而，大量实验表明归一化差分积雪指数可以区分雪和土壤、植被和薄云等地物，但是无法区分雪、卷云和水体，如表3所示。

根据官方提供的FY-4A号卫星数据的性能参数如表2所示，发现其第4波段（1.36~1.39 μm）设计目的是分离卷云，此属性也和地物反射率特征一致。NDSI对区分水体和积雪也不敏感，通过密度分割及样本采样，发现利用FY-4A号卫星数据的第5波段（1.58~1.64 μm）可以区分积雪和水体。综上所述，本文选择NDSI、第4波段（1.36~1.39 μm）和第5波段（1.58~1.64 μm）为积雪识别因子，并使用如图3所示的FY-4A动态多阈值积雪识别方法。

由于静止轨道卫星成像方式为全圆盘扫描、成像范围广且时间分辨率高。因此，同一观测时刻不同地区的太阳方位角和卫星观测方位角差异较大，造成同一地物在一天内的反射率不断变化，该现象对可见光波段的影响较为显著，这就对完成积雪识别的阈值选择产生了较大的影响。为了减小这种影响，根据研究区样本数据的实验结果，按照影像成像的季节和时间间隔设置了如表4和表5所示的动态阈值。通过影像获取的时间即可确定唯一对应的决策树阈值。

为了充分利用FY-4A数据高时间分辨率的优势，从而达到减少云量的目的。本文提出了如图4所示的FY-4A多时相云隙间积雪识别方法。首先，通过人工判别将数据分为两类，一类命名为晴间云数据（通过云流动可以判识云间隙中的积雪像元），另一类命名为云数据（研究区整天被卷云覆盖即无法判识云下或云隙间的积雪像元）。其次，针对晴间云数据进行了大量实验，发现如果一天内10~12景数据的同一像元被判识为积雪的次数为2~3次以上时，则认为该像元真实为积雪像元。最后，针对云数据，选择与其时间最接近的前后两景晴间云数据的面积的平均值代替该时间的积雪面积。

图像融合的目的是结合FY-4A积雪判识结果和FY-3C MULSS积雪覆盖日产品各自的优势，以更准确地反映祁连山区积雪面积变化和时空分布特征。融合后积雪判识结果以FY-4A号数据为准，每日融合图像FY3C4数据的第i行j列像元F_ij赋值方法如图5所示。①FY-3C数据中编码为0、1、11、254、255定义为无意义并且定义为0，若FY-3C数据像元编码为0时，融合后图像像元值F_ij⁴赋值为FY-4A动态多阈值多时相云隙间积雪识别方法的结果值；②当FY-3C数据像元值为云和雪（即像元值编码为50和200），融合后图像像元值F_ij⁴赋值为FY-4A动态多阈值多时相云隙间积雪判识方法的结果值，同时保留FY-3C的积雪像元；③当FY-3C值为陆地、水体、海洋（编码为37、39）等时，融合后图像F_ij⁴像元赋值为FY-3C MULSS积雪覆盖日产品的值。

使用Landsat 8 OLI数据和MOD10A2积雪产品，对积雪判识结果进行综合验证。Landsat 8 OLI数据具有较高的空间分辨率，并且在监测陆地资源时有出色的表现^［21］。因此，本文利用Landsat 8 OLI数据作为检验融合方法精度的一项指标。目前，研究积雪面积或时空演变时，国内外常采用的数据是MODIS系列的积雪产品，许多学者已经对该数据进行了精度分析，因此本文不再赘述其精度。

Landsat 8 OLI数据的时间分辨率为16 d，本文选择每月云量较小的数据对融合方法的精度加以验证。通过与Landsat 8 OLI数据叠加分析发现，融合后的FY3C4数据相较于单一的FY-3C MULSS积雪覆盖日产品能更加准确地实现积雪像元的判识如图6所示。从Landsat8 OLI数据中可以明显地看到，FY-3C MULSS积雪日产品对于山体边缘的积雪像元存在有错判和漏判的情况，而融合后的FY3C4数据可以较好地弥补这一缺陷，对于山体边缘的积雪像元也可以较为准确地判识。

MOD10A2积雪产品是由美国国家冰雪产品数据中心（NSIDC）提供的，空间分辨率为500 m^［22］。融合后的FY3C4数据空间分辨率为1 000 m。因此，需要将MOD10A2数据的空间分辨率重采样为融合后影像的空间分辨率，以避免出现由于空间位置错位而造成的错误验证。将MOD10A2积雪产品记录的数值作为FY3C4数据精度评价的地面真实值进行随机精度验证。总体精度（OA）是指在MOD10A2积雪产品判识为积雪同时融合后影像判识为积雪、MOD10A2积雪产品判识没有积雪和影像判识没有雪的像元数总和除以所有像元数。积雪非一致性检验包括多分误差（EU）和漏分误差（EO），EU表示积雪判识结果将无雪判识为有雪的误差比例，EO表示积雪判识结果将有雪判识为无雪的误差比例^［21-24］。

其中：SS为MOD10A2积雪产品和融合后影像均为积雪的像元数，DD为MOD10A2积雪产品和融合后影像均不为积雪的像元数，DS为MOD10A2积雪产品不是积雪而融合后影像判识为积雪的像元数量，SD为MOD10A2积雪产品是积雪而融合后影像判识不是积雪的像元数量，AC为除云以外的所有像元数量。经过计算，无云情况下FY3C4数据的平均总体精度、多分误差和漏分误差如表6所示。

通过使用Landsat 8 OLI数据和MOD10A2数据对积雪判识精度进行评价，发现融合FY-3C 数据和FY-4A数据可以大大降低光学数据中云对积雪判识造成的影响。如图7所示为2019年1月挑选一景为例的FY-3C MULSS积雪日产品、MOD10A2数据和FY3C4数据的积雪判识情况。积雪识别方法基本可以对于同一像元进行实时观测，充分利用了FY-4A气象卫星数据的高时间分辨率的优势，很大程度上减少了云对积雪判识产生的影响。从对比图中可以看出，相比于MOD10A2数据由于FY-4A数据低空间分辨率和波段少的问题，导致FY3C4对细碎的积雪识别能力不足存在有错判和多判的现象，但是与单一的FY-3C积雪产品相比，能更真实地反映地面上积雪覆盖的情况。

通过精度验证后，为了方便使用栅格数计算面积，对融合后的影像重新定义为等积投影，融合后影像的空间分辨率为1 000 m。因此，融合后影像的一个像元在等积投影下代表1 km²，像元统计提取出祁连山区日积雪面积如图8所示。

为方便探寻积雪面积变化的规律，以月为基本单位进行研究。祁连山区积雪面积的最大值出现在11月其面积值为11×10⁴ km²，积雪面积的最小值出现在7月其面积值为2×10³ km²，一年内积雪面积的最大值和最小值差距悬殊，极值相差10.8×10⁴ km²。积雪积累季节时段的面积值在年内面积平均值以上的时间高达96%。祁连山区年内积雪面积值分布大体上呈现“驼峰”形状，即在年内有两个波峰分别为11月和1月，一个波谷为7月。根据祁连山区多年气象资料显示，一般在每年9月中下旬开始为积雪积累期，经过45~50 d的积累期，在11月中下旬积雪面积值达到峰值，出现了第一个波峰。12月天气现象较为稳定。因此，相比较于11月的积雪面积值略微有所下降。而1月由于寒潮等极端天气的影响又出现了第二个波峰，2月份气温转暖、降水减少，导致祁连山区积雪面积逐渐减小，一直持续到7月中旬积雪面积达到最低值。

祁连山区气温随季节变化明显，总体呈现倒三角形特征。即每年1月中下旬或2月上旬至7月中旬气温逐渐回升、降水量逐渐下降。7月中旬至8月下旬祁连山区气温变化复杂，但是从总体上看，气温呈现降低的趋势。每年9月至次年1月中旬，气温逐渐下降、降水量逐渐增加。祁连山区积雪面积变化与气温呈显著相关性，气温升高时为积雪消融期，积雪面积持续下降直到最小值，气温下降时为积雪积累期，积雪面积达到最大值。

祁连山区积雪面积变化，从空间上看差异显著，总体上呈现出东部积雪面积大于西部积雪面积的特征。海拔和坡向是造成积雪空间差异的主要原因。对于海拔，海拔不同的区域气温和降水量也不同，根据4.2节中的分析，祁连山区积雪面积变化与气温和降水有显著相关性。因此，本文将祁连山区的海拔高度划分为4个高程带，分别为第一高程带2 000~3 000 m，第二高程带3 000~4 000 m，第三高程带4 000~5 000 m和第四高程带5 000 m以上，4个高程带分别占祁连山区总面积的12%、52.5%、34.8%和0.7%。总体来看，祁连山区积雪面积较大的区域集中第二高程带和第三高程带。原因在于，第一高程带海拔较低气温高，所以积雪面积值相对较小。虽然第四高程带海拔高气温低，但是此类地区占地面积较小，从而积雪面积较小。尽管如此，积雪覆盖率整体随着海拔的增加呈现增加的趋势。

对于坡度，迎风坡和背风坡、阴坡和阳坡所受到气流类型和太阳辐射量不同，导致其所在区域的气温和降水也不同。因此，不同坡度的积雪面积值不同。积雪面积值一般表现为迎风坡高于背风坡，阳坡低于阴坡。将研究区按照方位角重分类为8个坡向，分别为北坡（0°~22.5°及337.5°~360°）、东北坡（22.5°~67.5°）、东坡（67.5°~112.5°）、东南坡（112.5°~157.5°）、南坡（157.5°~202.5°）、西南坡（202.4°~247.5°）、西坡（247.5°~292．5°）以及西北坡（292.5°~337.5°），8个坡向分别占祁连山总面积的12.8%、14.8%、9.7%、18.7%、8.7%、16.7%、9.9%、8.7%。提取各坡向的积雪像元，进行积雪覆盖率计算并统计，发现祁连山区不同坡向积雪面积差异较大，西坡、西北坡、北坡以及东北坡积雪面积较大，远远多于东坡、东南坡、南坡以及西南坡。

我国风云系列气象卫星自1977年开始研发至今已经积累了很大的数据量。本文基于FY-3C数据和FY-4A数据的特性，提出了FY-4A动态多阈值多时相云间隙积雪识别方法，并且将该方法的判识结果与FY-3C积雪日产品融合达到了较高的积雪判识精度，这表明FY3C4积雪日数据在积雪监测中有着一定的应用价值。由于祁连山海拔高地势起伏较大，有时甚至导致研究区内云的分布高达90%以上，而利用光学数据研究积雪大多都受制于云。目前，MODIS系列的积雪产品和FY-3C MULSS积雪产品以其高空间分辨率的优势为积雪研究提供着服务，但是MODIS系列的积雪产品和FY-3C MULSS积雪产品重访周期长易受到云的影响，这种现象尤其体现在高海拔区域。在云层分布较多的区域，积雪判识精度就有所降低。FY-4A数据在其高时间分辨率的背景下，体现出对积雪判识的巨大优势，但是FY-4A用于构造积雪指数的波段空间分辨率相对较低，也制约了积雪判识的精度。祁连山区的云量相对低海拔区域的云量更大，经过影像融合后部分云像元被替换，在无云层覆盖的情况下总体精度能达到85.25%。

通过对祁连山区2018年3月~2019年3月积雪面积的监测，发现祁连山区积雪分布与其气候有显著相关性，研究发现积雪面积与气温呈负相关，与降水量呈正相关。在时间上，祁连山区积雪面积变化大体上呈现“驼峰”型，积雪面积变化季节性较强，具体表现为春季至夏季，温度上升降水量下降导致积雪面积持续下降，秋季至冬季温度下降降水量上升，积雪面积逐渐积累达到年内峰值。在空间上，祁连山区积雪主要分布在海拔3 000~5 000 m之间的范围内，积雪覆盖率随着海拔高度的增加呈显著正相关，对于不同的坡度，积雪覆盖分布也不同，其中西坡、西北坡、北坡以及东北坡积雪面积较大。

在对FY3C4数据积雪识别精度进行分析时，发现融合后数据的精度具有不确定性。主要体现在以下4个方面：①在研究区和数据方面，FY-4A数据空间分辨率低、波段数少且研究区内长时间被云层覆盖的频率较高，因此对小区域进行研究时混合像元分类精度低；②在样本选择方面，本文选择样本时仅考虑了研究区域内的样本，导致该方法对其他区域的积雪识别能力弱；③在积雪识别方法方面，积雪识别过程中需要人工对数据类型进行区分，从而引入了相应的误差；④在精度验证方面，MOD10A2数据是8 d合成产品与FY3C4数据的时间分辨率不一致也会增加精度验证的不确定性，在以后的研究中会针对这些问题对积雪识别方法做进一步改进。