赵明杰， 王宁练,3， 石晨烈， 侯靖琪
（1.陕西省地表系统与环境承载力重点实验室，陕西 西安 710127；
2.西北大学城市与环境学院地表系统与灾害研究院，陕西 西安 710127；
3.中国科学院青藏高原研究所青藏高原地球系统与资源环境全国重点实验室，北京 100101）
气候变化是世界范围内一个重要的共同课题［1］。据IPCC 第六次报告称，相对于1850 到1900年，全球地表温度已经升高1.09 ℃，并呈现出加速上升趋势［2］。气候变化对全球陆地和海洋生态系统产生了巨大而深远的影响，同时给人们生产生活带来了诸多不利影响，因此明确气候变化响应机制迫在眉睫［3］。自然界中有很多介质可以用来指示气候变化，比如海平面、冰川、冻土、湖泊等，湖泊在全球水循环中发挥着至关重要的作用，其湖泊面积、水位以及湖冰物候等都可以作为区域气候变化的指标［4］，其中湖冰物候对区域气候变化响应异常敏感［5-8］。高纬度地区的湖泊会随着天气变冷而产生湖冰［9-10］，湖冰物候变化对气候、环境和社会经济活动会产生重要的影响［11-13］。此外，湖冰物候与湖水能量平衡变化有直接的关系［14］，可以很好地反映区域乃至全球气候变化［15-16］。过去几十年来，北半球的湖冰正在持续减少［8］，而且随着未来气温升高，这种情况将会更加严重［10］。因此，研究湖冰物候特征及相关因素，对于进一步认识湖泊对全球气候变化的响应规律和反馈机理，具有十分重要的科学意义。
湖冰物候研究早期主要采用野外实地观测的方法来获取湖冰物候数据。但这种方法容易受到天气、环境、技术的影响，无法获取准确、长时间序列、大范围的湖冰物候数据。近年来，随着遥感技术的蓬勃发展，利用高时空分辨率的遥感影像来获取湖冰物候信息已经成为一种研究趋势［17］。基于遥感的湖冰物候监测主要有微波遥感和光学遥感。微波遥感又有被动微波和主动微波之分，其中被动微波遥感具有高时间分辨率（每天2 次或更多），但空间分辨率相对较低（>10 km），因此被动微波适合监测大型湖泊的湖冰变化，例如青藏高原上的青海湖［18-20］，以及加拿大的大熊湖和大奴湖［21］和中国、俄罗斯之间的界湖兴凯湖［22］。而主动微波遥感具有较高的空间分辨率，ERS（方位方向<30 m，距离方向<26.3 m）被用于监测湖泊结冰过程和冰厚度［23-24］；
Radarsat（1～100 m）被用于监测湖泊的冻结和消融过程［25-26］。相比被动微波遥感，主动微波遥感的时间分辨率较低（ERS 为3 d，Radarsat 为24 d）,不能满足对湖冰进行高频次检测的需求［27-28］。光学遥感如MODIS影像，其同时具有高时间和高空间分辨率（1 d，空间分辨率为250 m），被广泛应用于监测湖冰物候变化。如：邰雪楠等［29］使用MODIS 数据，分析了2000—2020 年色林错湖冰物候特征及其影响因素。Yao 等［30］利用MODIS 和Landsat 数据提取了2000—2011 年可可西里地区22 个湖泊的湖冰物候，并对相关的影响因子进行了讨论。吴艳红等［31］使用MODIS数据，提取了纳木错湖2000—2015年间的湖冰物候数据，基于相关模型重建了1963—2018年纳木错湖的湖冰物候序列，分析了60多年来纳木错湖冰物候的变化特点。
先前湖冰物候的研究多集中在高海拔高纬度地区，对于中亚地区湖冰物候的研究较少。中亚位于欧亚大陆中部，75%的地区属于干旱半干旱的大陆性气候［32］，其中包括数千个湖泊，这些湖泊对中亚地区生态系统的可持续性和人类福祉有着重要的意义［33］。中亚湖泊的相关研究主要集中在湖泊面积以及水位与水量的变化。如Che 等［33］利用terraPulse™月度Landsat衍生的地表水域范围数据集和HydroLAKES 数据集，对2000—2015 年中亚湖泊面积的时空变化进行了研究。Huang 等［34］利用SRTM和Landsat 影像构建了面积-水位-水量变化的经验模型，并估计9619 个湖泊中最大水量和最小水量。Hu 等［35］利用3 个GRACE 卫星数据集和5 个全球水文模型对2003—2014 年中亚干旱区陆地储水量变化进行研究。对于中亚地区湖泊湖冰物候的研究相对较少，而且精度较差。因此，本文基于MODIS地表反射率数据，提取了2000—2020年中亚地区湖泊的湖冰物候信息，分析了湖冰物候变化特征，并利用气象数据和现有湖泊资料，讨论了气温、降水、湖泊面积以及海拔对湖冰物候的影响。
中亚一般意义上包括土库曼斯坦、哈萨克斯坦、吉尔吉斯斯坦、塔吉克斯坦和乌兹别克斯坦。中亚的地势总体上东南高、西北低，气候干燥，降水稀少，昼夜温差较大［36-37］，占世界干旱半干旱地区总面积的34%［38］，是全球最大的非地带性干旱区［39］（图1）。中亚干旱地区缺水严重，可供开发的淡水资源很少［40］，因此对于中亚地区，湖泊是非常重要的水资源［41］。中亚湖泊总面积超过88000 km2，其中约3000 个湖泊面积>1 km2，40 多个湖泊面积超过100 km2［42-43］。由于独特的气候，中亚大多数湖泊的水源来自冰川融化、山区降水和河流径流［44］。
图1 研究区湖泊分布Fig.1 Distribution of lakes in the study area
其中选择7 个湖冰存在期稳定的并且面积>100 km2湖泊（云量较少，像素误差分类较少，湖泊冻结消融明显）作为研究湖泊，它们分别是卡拉库尔湖、巴尔喀什湖、咸海、阿拉湖、斋桑泊、查蒂尔-科尔湖以及马卡科尔湖。表1 为7 个湖泊的详细信息，包括海拔高度、面积大小以及湖泊类型。
表1 研究中使用的湖泊信息Tab.1 Lake information used in the study
2.1 数据
2.1.2 Landsat 数据Landsat 是美国航空航天局与美国地质调查局的联合计划，该计划提供了现有的世界上最长的连续空基观测记录，从计划开始到现在，共发射了9颗卫星。本研究用到的Landsat数据（时间分辨率为16 d，空间分辨率为30 m）主要包括：Landsat7、Landsat8（表2），通过获取MODIS 和Landsat同时过境的数据，来检验MODIS数据提取湖冰物侯属性的精度。
表2 研究中使用的Landsat数据Tab.2 Landsat data used in the study
2.1.3 再分析气象数据由于中亚地区气象站稀少且连续观测数据周期短，气象站数据不足以支撑研究湖冰物候变化的因果分析。因此，本研究使用再分析气象数据（CRUTS v4.05）来分析湖冰物候变化的影响因素。CRUTS v4.05数据（时间分辨率为1个月，空间分辨率为0.5°），包括1901—2020年逐月的平均气温与降水量，在中亚气候研究中表现出良好的适用性［45-46］。本研究获取了2000—2020 年的气温以及降水数据来对湖冰物候变化进行相关性分析。
2.2 方法
2.2.1 湖冰信息提取冰在可见光与近红外波段的反射率高，而水在可见光和近红外波段的反射率较低［47］，因此根据冰和水在红光与近红外波段的反射率差异，利用直方图和目视解译，找到合适阈值来区分冰和水，该方法为阈值法［19］。具体公式为：
式中：Band1、Band2 分别为MODIS 影像250 m 的红光、近红外波段；
a、b 为阈值，当Band1-Band2>a、Band1>b时，则可以认为是湖冰。
NDSI（Normalized difference snow index）方法也可以用于提取湖冰信息［17］。本文使用NDSI 方法提取Landsat数据的湖冰特征信息，其公式如下：
式中：Band2 为绿光波段；
Band5 为短波红外波段。依据计算出来的NDSI 结合反射率直方图以及目视解译最终确定阈值，进而区分出湖冰和湖水。
2.2.2 湖冰时间属性定义湖冰时间属性具体分为：湖泊开始冻结日期（FUS）、完全冻结日期（FUE）、开始消融日期（BUS）、完全消融日期（BUE）。湖冰冻结期（FUD）指湖泊从开始冻结直至完全冻结所用的时间；
湖冰消融期（BUD）指湖泊从开始消融直至完全消融所用的时间；
湖泊湖冰存在期（ICD）指从湖泊开始冻结到完全消融这段时间；
湖泊完全冻结期（CFD）指从湖泊完全冻结至开始消融的这段时期［48］。一般情况下，在湖泊冻结期间，湖泊边缘附近会出现部分未冻结的区域，在湖冰消融过程中，湖岸处也可能会存在零星的湖冰［30,49］。此外，由于天气变化、湖泊边界不匹配、像素分类错误以及不可避免的噪声（尤其是云层）导致的反复冻结和消融，会降低湖冰信息提取的准确性。因此，本研究根据Kropáek 等［48］的方法，将湖冰占湖泊面积的10%和90%作为阈值进行湖冰物候的提取。具体公式如下：
(ⅱ) 对任意F∈CIrr(X),若Fδ∩(∪i∈IUi)=∪i∈I(Fδ∩Ui)≠Ø,则存在i∈I使得Fδ∩Ui≠Ø,由Ui∈τCSI及F∈CIrr(X),F∩Ui≠Ø,于是F∩(∪i∈IUi)≠Ø,从而∪i∈IUi∈τCSI。
3.1 MODIS 与Landsat 提取湖冰物候信息的交叉验证
为了评估阈值法提取MODIS湖冰信息的精度，本文使用了Landsat数据进行了验证。在本研究中，选用了测量平均绝对误差（MAE）、测定校正（R2）和偏差（Bias）等指标来评估了MODIS数据提取湖冰物候的精度。首先，选取了20幅没有云层干扰以及湖泊处于冻结与消融过程中的Landsat影像，然后分别计算出每一幅影像中湖冰面积占整个湖泊面积的百分比。最后将这些比例与MODIS 数据提取的湖冰覆盖百分比进行对比。
图2 为Landsat 数据与MODIS 数据提取湖冰面积比例的对比，其结果为：R2为0.999，MAE为0.92%，Bias 为1.15%，这表示MODIS 数据根据阈值法提取中亚湖泊的湖冰面积精度较高，可用于提取中亚地区大型湖泊的湖冰物候信息。图3 为基于NDSI 与阈值法提取2017 年12 月22 日卡拉库尔湖的湖冰信息，其中Landsat 影像使用NDSI 方法提取的湖冰面积为297.69 km2，MODIS 数据基于阈值法提取的湖冰面积为286.15 km2，2 种提取方法的误差约为0.58%。
图2 Landsat数据与MODIS数据湖冰面积比例对比Fig.2 Comparison of lake ice area ratio between Landsat data and MODIS data
图3 NDSI与阈值法提取湖冰信息Fig.3 Extraction of lake ice information by NDSI and threshold method
3.2 2000—2020年中亚7湖的湖冰物候特征分析
基于MODIS影像，使用阈值法结合人工目视解译，得到2000—2020年中亚地区所选湖泊的湖冰物候平均特征（表3、图4）。
表3 2000—2020年中亚地区所选湖泊平均湖冰物候统计Tab.3 Average lake ice phenology of selected lakes in Central Asia from 2000 to 2020/d
图4 2000—2020年中亚地区所选湖泊平均湖冰物候特征Fig.4 Average lake ice phenological characteritics of selected lakes in Central Asia from 2000 to 2020
从图4 和表3 可以看出，7 个湖泊从每年9—11月开始结冰，其中阿拉湖和查蒂尔-科尔湖开始冻结日期是早于7个湖泊的平均冻结日期。阿拉湖开始冻结日期最早，咸海开始冻结日期最晚。7 个湖泊的平均完全冻结时间是在当年的11月底到12月底，平均冻结期为35 d。
7 个湖泊的湖冰在次年的3—4 月开始消融，其中巴尔喀什湖和阿拉湖日期是早于7个湖泊的平均消融日期。阿拉湖日期开始消融时间最早，而查蒂尔-科尔湖开始消融的时间最晚。7个湖泊的平均完全消融时间是在次年的4—6 月，平均消融期为18 d。7个湖泊中，有6个湖泊的冻结期比消融期长，说明中亚地区所选湖泊冻结的速度比消融的速度慢。
7个湖泊在湖冰存在期和完全冻结期上差异很大。湖泊平均湖冰存在期为171 d，其中咸海湖冰的湖冰存在期最短为127 d，查蒂尔-科尔湖冰存在期最长为237 d；
湖泊平均完全冻结期为126 d，其中查蒂尔-科尔湖完全冻结期最长为170 d，巴尔喀什湖完全冻结期最短为94 d。
3.3 2000—2020年中亚7湖的湖冰物候变化趋势
2000—2020年7个湖泊湖冰物候变化趋势存在较大的差异（表4、图5～6）。从表4 和图5 中可以看出，2000—2020 年除了巴尔喀什湖［1.44 d·(10a)-1］和卡拉库尔湖［0.30 d·(10a)-1］的FUS 呈缓慢提前趋势外，剩下5个湖的冻结日期都是呈延后趋势，平均延后率为4.86 d·(10a)-1，其中查蒂尔-科尔湖推迟最为明显，推迟率为18.00 d·(10a)-1。BUE除了卡拉库尔湖外缓慢延后外，其余的呈现提前的趋势，平均提前率为2.90 d·(10a)-1，其中阿拉湖提前最明显，提前率为5.70 d·(10a)-1。
表4 2000—2020中亚地区所选湖泊湖冰物候变化趋势Tab.4 Change trend of lake ice phenology in selected lakes in Central Asia from 2000 to 2020/d·(10a)-1
图5 2000—2020中亚地区所选湖泊湖冰物候变化趋势空间分布Fig.5 Spatial distribution of lake ice phenology change trend in selected lakes in Central Asia from 2000 to 2020
图6 2000—2020年中亚所选湖泊湖冰物候年际变化Fig.6 Interannual variation of lake ice phenology in selected lakes in Central Asia from 2000 to 2020
2000—2020年查蒂尔-科尔湖、马卡科尔湖、阿拉湖和斋桑泊的湖冰存在期呈现出缩短趋势。例如，阿拉湖湖冰存在期缩短趋势为8.7 d·(10a)-1；
马卡科尔湖、阿拉湖以及斋桑泊的完全冻结期和湖冰存在期都呈现缩短趋势。在7 个湖泊中，除了阿拉湖和咸海的冻结期与消融期的趋势相反，其余5 个湖泊在冻结期和消融期方面表现出相同的趋势。
4.1 湖泊冻结-消融空间模式特征
通过对中亚地区所选6 个湖泊的冻结-消融过程统计分析，湖泊冻结-消融空间格局分为2 种：（1）冻结自湖岸一侧延伸至另一侧，越先冻结的湖区越先消融。符合这种空间格局的湖泊包括：卡拉库尔湖、巴尔喀什湖、阿拉湖、斋桑泊以及马卡科尔湖。为了直观展示这种模式的湖冰物候演变过程，本研究选用巴尔喀什湖为例，如图7 和图8［50］所示，由于巴尔喀什湖西部淡水区的盐度通常低于0.5 g·L-1，而东部半咸水区的盐度可以高达3 g·L-1以上，所以巴尔喀什湖于11 月下旬左右先在湖西南方向开始冻结，湖冰从西南向东北缓慢覆盖，到12 月下旬完全被冰覆盖，整个湖冰冻结期比较漫长。湖冰消融时，主要从湖的西南岸先消融，然后向东北逐渐消融，到4 月中旬，湖冰完全消融。（2）湖水从两岸向湖心逐渐冻结，消融时从一侧到另一侧，这种模式的湖泊有查蒂尔-科尔湖。图9 为查蒂尔-科尔湖湖冰物候演化过程，查蒂尔-科尔湖于10 月上旬左右先在湖南北两岸形成岸冰，湖冰由两岸向湖心逐渐蔓延，到12 月下旬完全被冰覆盖，整个湖冰冻结期比较长。湖冰消融时，主要从湖的南岸先消融，由南向北迅速消融，到5 月中旬，湖冰完全消融。
图7 2010—2011年巴尔喀什湖冻结-消融过程Fig.7 Freezing-thawing process of Balkhash Lake from 2010 to 2011
图8 巴尔喀什湖盐度分布Fig.8 Salinity distribution map of Balkhash Lake
图9 2010—2011年查蒂尔-科尔湖冻结-消融过程Fig.9 Freezing-thawing process of Chatir Kol from 2010 to 2011
4.2 气象因素对湖冰物候变化的影响
湖冰物候的形成受到气候条件和湖泊理化特征的影响［18,51］，其中对湖冰物候产生影响的气象因素有：气温、降水、太阳辐射、风速、风向等，以往在其他地区的研究发现，气温是影响湖冰物候的关键因素［52］。为了研究中亚地区选定湖泊的湖冰物候变化与气候变化之间的响应规律，本研究分析了气温和降水量与湖冰物候之间的相关性。
为了验证再分析气象数据的准确性，本研究使用2013年6月到2020年12月卡拉库尔湖气象站的月平均气温与CRUTS v4.05对应的数据进行了交叉验证，如图10 所示R2为0.961，MAE 为0.211%，Bias为1.807%，这表明CRUTS v4.05 数据集可用于讨论中亚地区的气候变化状况。
图10 CRUTS v4.05数据集精度验证Fig.10 Precision validation of CRUTS v4.05 dataset
图11为湖冰物候与气温、降水的相关性。由图11可知，湖冰的存在期、完全冻结期、开始和完全融化期与年平均气温有较强的相关性。年平均气温较低时，湖冰的存在期和完全冻结期较长，年平均气温较高时，湖冰的存在期较短。此外，中亚湖泊开始和完全融化的日期与年平均气温呈较强的负相关。这表明气温是中亚湖泊湖冰物候变化的决定因素。此外，分析结果还表明，降水对中亚地区湖泊湖冰物候的变化影响不大。
图11 湖冰变化物候与气象因子变化的相关性热图Fig.11 Correlation heat map of lake ice change phenology and changes in meteorological factors
4.3 湖水理化性质对湖冰物候变化的影响
4.3.1 湖冰物候特征与湖泊面积之间的关系湖泊面积会影响湖冰物候的变化［53］。图12为近20 a 7个湖泊的平均面积与FUS、BUE、ICD 的相关性。从图中可以看出，中亚地区7 个湖泊面积与FUS 呈正相关，与BUE呈负相关，这表明湖泊面积越大，湖水开始结冰的时间愈迟，融化愈快，湖冰的存在期就越短。
图12 中亚地区湖泊湖冰物候与湖泊面积关系图Fig.12 Relationship between lake ice phenology and lake area in Central Asia
图13 为2000—2020 年中亚湖泊面积与FUS 变化率、BUE 变化率、ICD 变化率的相关性，从图中可以看出面积与FUS变化率、BUE变化率、ICD变化率的相关系数分别为-0.42、0.24、-0.45，湖泊面积与FUS 和ICD 变化率呈现负相关，这说明湖泊面积越大，湖泊的FUS和ICD的变化速率越慢；
湖泊面积与BUE 变化率呈现正相关，表明湖泊面积越大，BUE变化速率越快。
图13 中亚地区湖泊湖冰物候平均变化率与湖泊面积关系图Fig.13 Relationship between the average change rate of lake ice phenology and lake area in Central Asia
4.3.2 湖冰物候特征与湖泊海拔之间的关系湖面
海拔会影响湖冰物候的变化［50］。图14 为湖泊海拔与FUS、BUE、ICD 的相关性。从图中可以看出海拔与FUS 呈负相关，与BUE 和ICD 呈显著正相关。这说明湖面海拔越高，湖水开始结冰的时间越早，消融的时间越晚，湖冰存在期越长。
图14 中亚地区湖泊湖冰物候与湖泊海拔关系图Fig.14 Relationship between lake ice phenology and lake altitude in Central Asia
图15 为2000—2020 年中亚湖泊海拔与FUS 变化率、BUE变化率和ICD变化率的相关性，从图中可以看出海拔与FUS、BUE、ICD 的相关系数分别为0.49、0.63、0.43，都呈现出正相关趋势，这说明随着海拔的升高，湖泊的FUS、BUE以及ICD的变化速率都在加快。
图15 中亚地区湖泊湖冰物候平均变化率与湖泊海拔关系图Fig.15 Relationship between the average change rate of lake ice phenology and lake altitude in Central Asia
4.4 无量纲化分析
无量纲化是一种数据分析方法，旨在将数据集中的多个变量转换为一组无单位的变量，以便更好地理解和分析数据。这个过程可以减少数据的复杂性，提高数据的处理效率，并提高数据分析的精度和准确性。本文采用标准化进行无量纲分析，具体公式如下：
式中：xj为n个训练样本中第j个特征值组成的向量；
μj为训练样本中的均值；
σj为训练样本的标准差。
4.4.1 单个湖泊的无量纲化分析以巴尔喀什湖为例，将湖冰物候信息、冬半年气温、降水量、面积采用标准化方式进行无量纲化处理，然后进行多元回归分析，进而分析出冬半年气温、降水量、面积对湖冰物候的影响程度。具体分析结果如下：
式中：X为冬半年气温；
Y为降水量；
Z为面积；
*、**、***表示通过P值为0.1、0.05、0.01 的显著性检验。结果表明，气温是影响湖冰物候的关键因素，面积主要影响湖冰物候的FUE和CFD，降水量对湖冰物候的影响相对较小。
4.4.2 不同湖泊的无量纲化分析以6 个湖泊2000—2020年的平均湖冰物候、平均湖泊面积以及湖泊海拔进行无量纲处理，然后进行多元回归分析，进而分析出面积和海拔对湖冰物候的影响程度。具体分析结果如下：
式中：X为海拔；
Y为面积。结果表明，除了FUE主要受湖泊面积影响外，其他湖冰物候主要受海拔影响。
4.5 现有研究结果对比
通过文献搜索查阅发现，研究中亚地区湖冰物候变化的文章非常少，我们统计了现有研究结果和本研究的结果的差异（表5）。Hou 等［54］使用了MODIS 地表温度数据，通过设置一定阈值来提取湖冰的范围；
而本研究采用MODIS 地表反射率数据，通过设置动态阈值来提取湖冰。前者虽然可以一天获取多景数据来减少云的影响，但由于不同湖泊自身理化性质（如盐度等）的差异，使用固定阈值带提取湖冰物候属性会存在一定的误差。而本研究采用250 m的反射率数据，具有更高的空间分辨率，且动态的设置阈值，其精度更高。
表5 本研究与现有研究结果对比Tab.5 Comparison of the results of this study with existing research
本文基于MODIS 与Landsat 数据并结合相关的气象数据，分析了中亚地区7 个湖泊的长期湖冰物候特征及其影响因素，主要结论如下：
（1）2000—2020 年中亚的湖泊从每年9—11 月开始结冰，湖冰在3—6 月逐渐完全消融。7 个湖整体平均冻结期和平均消融期为35 d 和18 d，湖冰平均完全冻结期和平均存在期为126 d和171 d。
（2）2000—2020 年中亚7 个湖泊中有5 个湖泊开始冻结日期呈现延后的趋势，平均延后率为4.86 d·(10a)-1。总体上，完全消融日期有提前的趋势，平均提前率为2.90 d·(10a)-1。其中有4个湖泊湖冰存在期呈缩短趋势，平均缩短率为9.90 d·(10a)-1。完全冻结期呈整体缩短趋势，平均缩短率为5.00 d·(10a)-1。
（3）中亚7 个湖泊湖冰的冻结-消融空间模式主要分为两类：湖水自湖岸冻结至对岸，越先冻结的湖区越先消融；
湖水从两岸向湖心逐渐冻结，消融时从湖岸到对岸。
（4）中亚地区湖冰物候变化受到多种因素的综合影响，包括湖泊自身的特征（如海拔、面积等）和气候变化（如气温和降水量等）。气温是影响湖冰物候的关键因素，气温越高，湖冰存在期越短；
面积越大，湖泊的湖冰存在期越短；
随着海拔的升高，湖泊的湖冰存在期越长。
