1979年—2019年兴凯湖湖冰物候变化的被动微波遥感监测

柯长青，蔡宇，肖瑶

南京大学 地理与海洋科学学院,南京210023

1 引 言

湖冰有很强的季节性，气温变化能够直观地反映在湖冰冻结和消融的时间上，因此湖冰覆盖状况、冻结消融时间等既是良好的气候变化指示器（汪关信等，2020），也是冬季气温变化的指示器，湖冰有时甚至比短期内的气温变化更能反映气候的变化趋势（Kouraev 等，2007）。同时湖冰冻结与消融速率、湖冰厚度、湖冰温度梯度等也是区域气候变化的重要指标（Gould 和Jeffries，2005）。过去几十年来，北半球的湖冰在不断地消失（Sharma 等，2019），而且在未来升温的情境下这种情况会更加严峻（Sharma 等，2020）。传统的野外观测活动和气象、水文站的记录虽然可以获取湖冰相关参数，但是野外观测耗时耗力，尤其是位置偏远的大型湖泊。气象站、水文站的分布并不均匀，而且这些观测站一般位于湖泊岸边，由于视线的限制，仅能够记录观测站视野范围内的湖泊冻结与消融状态，很难获得大范围的湖冰记录。因而对于面积比较大的湖泊，地面观测记录无法准确反映整个湖泊的湖冰冻结与消融变化。

利用遥感技术监测湖冰是比较快捷方便的方法，它能获取长时间序列、大范围的湖冰观测数据，弥补了人工观测、以及水文与气象观测站点分布不均匀的缺点，为湖冰监测研究提供了非常有效的技术方法（魏秋方和叶庆华，2010）。陈贤章等（1995）利用高级甚高分辨率辐射计AVHRR（Advanced Very High Resolution Radiometer） 数 据分析过青海湖湖冰的冻结消融过程，并用观测站点数据对湖冰厚度与气温之间的关系进行过探讨。殷青军和杨英莲（2005）基于中等分辨率成像光谱仪MODIS（Moderate-resolution Imaging Spectroradiometer）数据通过设定阈值来自动提取湖冰覆盖范围，并用气象观测数据进行了验证。Cai 等（2019）利用2000年—2017年的MODIS 数据分析了青藏高原50 多个湖泊冻融日期和湖冰覆盖时间的变化，并进一步探讨了湖冰物候变化的时空差异和原因。尽管AVHRR、MODIS 等中等分辨率的光学遥感数据有高时间分辨率的优势，可以提取每日的湖冰覆盖范围并进行制图（Weber 等，2016；祁苗苗等，2018；Wu 等，2021），进而可以得到湖泊冻融日期，但是会受到云的影响，使最终的日期结果产生一定的误差。

与光学遥感技术方法不同，微波遥感不受云雨天气条件的影响，而且被动微波遥感有较高的重访周期，可以获取每日的湖泊表面亮温信息，根据湖冰和湖水显著的亮温差异，即可提取湖泊的冻融日期（车涛等，2009）。虽然被动微波遥感数据空间分辨率通常较低，但在监测大型湖泊的湖冰覆盖范围方面有很大的优势（汪关信等，2021）。Cai 等（2017）利用1979年—2016年被动微波遥感数据分析了青海湖湖冰冻结消融的长时间序列变化，而Howell 等（2009）利用SeaWinds/QuickSCAT微波遥感数据监测过加拿大两个大型湖泊Great Bear Lake 和Great Slave Lake 的湖冰物候变化。

国内利用遥感技术开展湖冰物候变化的监测研究主要集中在青藏高原湖区（Ke 等，2013；姚晓军 等，2015）和新疆湖区（Cai 等，2020），而东北湖区的湖冰遥感监测研究较少（Yang 等，2019）。兴凯湖（本文指大兴凯湖）面积大，湿地分布广，是三江平原地区生物多样性最为丰富的自然生态系统，是国家级自然保护区，同时也是世界生物圈保护区。兴凯湖是以湖泊、沼泽、草地和森林组成的一个完整湿地生态系统，是东北亚重要的鸟类停歇繁殖基地和迁徙通道，2001年入选拉姆尔国际湿地公约组织的国际重要湿地名录。因此开展兴凯湖的湖冰物候变化对全球气候变暖的响应研究具有重要的科学意义。Yang 等（2019）利用2003年—2016年MODIS 数据分析了兴凯湖的湖冰物候变化，然而MODIS 数据的时间序列短，而且还会受到云的干扰。因而本文利用1979年—2019年的被动微波遥感数据分析兴凯湖湖冰物候的长时间序列变化及其原因。

2 兴凯湖概况

兴凯湖位于黑龙江省东南部的乌苏里江流域，是中国、俄罗斯之间的界湖，介于44°32′N—45°21′N，131°58′E—132°51′E。兴凯湖是由造山运动导致地壳陷落而形成的构造湖，呈椭圆形，南北长（约90 km），东西窄（约50 km）。湖面海拔69 m，正常蓄水位的湖泊面积约为4010 km2，水深为2—10 m。兴凯湖是河流直接补给的湖泊，入湖河流有10 余条，而出湖河流仅有唯一的松阿察河。兴凯湖属于湿润半湿润的温带大陆性季风气候，年平均气温2.9—3.1 ℃，年平均降水量750 mm，夏季降水量约占全年的70%。兴凯湖冬季极其寒冷，11 月开始结冰，15 d 左右湖面全部冻结，湖冰厚度0.8—1.5 m，积雪深度可达20 cm，4 月中下旬解冻，湖冰覆盖时间150 d 左右。本文设定一个湖冰年周期从09-01 开始至下一年的08-31，例如，2018-09-01—2019-08-31 称为2019年的湖冰周期。

3 数 据

3.1 被动微波遥感数据

用美国雪冰中心（https：//nsidc.org/data/［2021-04-29］）CETB（Calibrated Enhanced Resolution Passive Microwave）产品中的SMMR（Scanning Multi-channel Microwave Radiometer）、SSM/I（Special Sensor Microwave Image）和SSMIS（Special Sensor Microwave Imager/Sounder）数据来提取兴凯湖的湖冰物候。其中美国雨云Nimbus-7 卫星搭载的SMMR 传感器，共有6.6 GHz、10.7 GHz、18.0 GHz、21.0 GHz 和37.0 GHz 5个频段，时间分辨率是1 d，空间分辨率是25 km。SSM/I 是搭载在DMSP（Defense Meteorological Satellite Program）系列卫星F8、F10、F11、F13、F14、F15 上的传感器，SSMIS 搭载在同系列F16、F17、F18、F19等卫星上，其中F15、F16、F17、F18仍在运行，共有19.35 GHz、22.2 GHz、37.0 GHz、85.5 GHz（SMMIS 上用91.655 GHz 通道代替85.5 GHz）4 个频段。各频段时间分辨率均为1 d。

CETB 是NASA（National Aeronautics and Space Administration）MEaSUREs（Making Earth System Data Records for Use in Research Environments）计划的一部分，是一种改进的格网地球系统数据记录，用于冰冻圈和水圈动态变化的监测。该数据集根据目前最成熟的星载被动微波二级数据生成，而且经过了跨传感器校准和质量检查，拥有更好的质量控制和改进的投影网格，并进行了本地时间的改进处理。数据集包含所有频段的25 km 分辨率产品，同时根据通道频段不同，提供最高可达3.125 km的增强分辨率产品，是目前空间分辨率最高的被动微波亮温数据。本文使用3.125 km 分辨率的CETB 产品，时间序列为1979-10-26—2019-08-31。由于1986-04-02—06-23 的Nimbus-7 卫星数据缺失，因此无法获得该时间段兴凯湖的开始消融和完全消融日期。

3.2 MODIS数据

MODIS 每日积雪产品来自美国雪冰中心（https：//nsidc.org/data/［2021-04-29］）。美国NASA地球观测系统的Terra 和Aqua 卫星上均搭载有MODIS 传感器，共有36 个波段，空间分辨率包括250 m、500 m和1 000 m。基于归一化积雪指数NDSI（Normalized Difference Snow Index）方法生成了MODIS 每日积雪产品MOD10A1/Terra 和MYD10A1/Aqua，提供空间分辨率为500 m的分类后地表覆盖类型数据，包括陆地、云层、积雪等类型。其中MOD10A1 产 品自2000-02-24 起，MYD10A1 产品自2002-07-04 起。NASA 于2017年发布了第6 个版本的MODIS 积雪产品，提供像元百分比的积雪覆盖。2000-09—2019-08 第六版MODIS 每日积雪产品也被用于提取兴凯湖的湖冰物候，主要目的是验证CETB数据提取的湖冰物候结果。

3.3 气候数据

使用兴凯湖附近鸡西气象站的数据（http：//data.cma.cn/［2021-04-29］）分析湖冰物候变化的原因。该站点位于130.57°E，45.17°N，获取了1979-09—2019-08的气温、降水和风速数据。与湖冰物候周期一致，采用9 月至次年8 月的数据计算年均气温、年均风速和年降水量。

4 方 法

4.1 被动微波数据湖泊冻融日期提取

冰和水的介电特性有显著差异，湖冰冻结和融化过程中，微波亮温在多个频段都会发生急剧变化，冰的亮温高于水的亮温，因此可以利用湖泊在有冰和无冰期间显著的亮温差异来获取湖冰日期。滑动t检验MTT（Moving T Test）阈值法在湖冰物候提取的过程中被证明非常有效（Du 等，2017）。MTT 方法对时间序列上的每个点依次进行t检验，检查时间点前后一定长度的子样本均值是否具有显著差异，进而获得突变点，常用于气候因素的突变检测（Xiao 和Li，2007）。获取湖冰冻融日期序列上的突变点时，判断某日前后20 d 的平均亮温是否有显著差异，进而判断该日是否发生了湖冰覆盖状况的变化。

湖泊的冻结和融化是一个变化的过程，在一个年周期内，至少会有两组亮温突变点（冻结和融化过程）。此外，由于湖冰的反复冻结、湖泊表面粗糙度变化等噪声存在，通常存在两组以上的连续突变点。对于一组突变点，第一个时间点前20 d 的平均亮温，和最后一个时间点后20 d 的平均亮温，分别可以代表变化过程前后的稳定状态。通常大型湖泊在形成湖冰覆盖前，会经历几天至几十天的冷却过程，这段时间湖泊的亮温处于全年最低时期。因此，选取所有突变组中最低的变化前亮温为湖水的亮温参考值，同组的变化后亮温为湖冰的亮温参考值。计算湖冰和湖水亮温参考值的均值作为划分湖冰与湖水状态的阈值。考虑到短期天气变化或传感器观测周期变化都可能引起湖泊亮温发生短暂变化，采用前后10 d 的亮温均值代替原始亮温，并与阈值比较，进而判断当日的湖冰/湖水状态。

由于CETB 产品原始观测数据的足迹点为25 km，导致靠近湖泊边界的像元很可能包含陆地混合信息，不能真实表现湖泊亮温在水/冰转换时的变化。为了避免陆地的干扰，并尽可能多地获取可用的、纯净像元，设置两个像元大小（6.25 km）的缓冲区，排除严重受陆地影响的像元，进而可以保证湖冰物候的提取精度。对于经过筛选的像元，统计每天湖冰和湖水的像元数量，设置年内湖水像元最大值为湖泊的总像元量，计算5%和95%的阈值，用于提取湖泊冻融日期（Kropáček 等，2013）。当湖冰像元数量大于5%且保持冻结状态，认为湖泊开始冻结；当湖冰像元数量大于95%且保持冻结状态，认为湖泊完全冻结；反之，当湖冰像元数量小于95%，认为湖冰开始消融；当湖冰像元数量小于5%，认为湖冰完全消融。

以2018年—2019年兴凯湖的湖冰冻融过程为例，6.25 km 的缓冲区可以筛选得到良好的湖泊冻融过程曲线（图1、图2（a））。兴凯湖共有246个像元参与湖泊冻融日期提取，在观测到超过12 个湖冰像元（5%）时，就认为该年开始冻结；湖冰像元达到234 个（95%）时，认为完全冻结。反之，春季湖水像元超过12个，就认为开始消融；湖水像元超过234个时，就认为完全消融（图1）。2018年—2019年兴凯湖冻融日期分别为：2018-11-21 开始冻结（图2（b）），2018-12-03 完全冻结（图2（c））；2019-03-27 开始消融（图2（d）），2019-04-09 完全消融（图2（e））。

图1 2018年—2019年兴凯湖冻融日期的提取过程Fig.1 The extraction process of freeze-thaw dates of Khanka Lake from 2018 to 2019

4.2 MODIS数据去云处理与湖泊冻融日期提取

由于Terra 是上午星，Aqua 是下午星，可以通过双星合成法和临近日五天合成法两个步骤完成去云（Gafurov 和Bárdossy，2009）。由于两颗卫星的发射时间不同，只能对2002-07-04 以后的数据进行去云处理。以2018-09—2019-08 MODIS 数据的去云处理为例，处理前Terra 卫星的年均云覆盖率为49.78%，第1 步双星合成法去云后，平均云量降低到了40.33%，第2 步临近日五天合成法去云后，平均云量只剩10.24%。

由于湖泊面积变化和不可避免的像元分类错误，首先对兴凯湖范围内每年的湖水像元数量取组中值，将其划分为大于组中值和小于组中值的两组，分别计算均值作为开阔水域状态和湖冰覆盖状态的像元数量参考值，利用两个参考值的5%和95% 作为阈值提取每年的湖泊冻融日期（Kropáček 等，2013）。为了避免湖面状况突变影响提取过程，去云后云量仍大于50%的影像不参与湖泊冻融日期的提取。

4.3 误差及不确定性分析

被动微波遥感数据提取湖冰物候的误差来源于两个方面：卫星运行模式导致的周期性数据缺失，以及湖泊冻融日期提取过程产生的误差。理论上被动微波数据的时间分辨率为1 d，但由于存在轨道条带间隔，实际上并非如此。采样间隔在1979年—1987年的数据中产生4—6 d 的误差，1988年—2019年的数据中产生1—3 d的误差（Cai等，2017）。总体上，随着传感器更迭，被动微波数据的采样频率在逐渐提高，数据缺失情况在逐渐改善。在湖冰物候提取过程中，如果遇到数据缺失，会向后顺延采用有数据的一日作为当年的结果。因此，被动微波获取的日期可能会比实际情况稍晚。

湖泊冻融日期提取过程中的误差主要由混合像元造成。以SSM/I为例，虽然增强分辨率产品的空间分辨率为3.125 km，但原始数据的足迹点约为25 km。因此，单个像元可以获取多种不同地表特征信息（Bennartz，1999），混合像元效应严重，尤其是靠近湖岸的像元，不可避免会混入陆地信息。纯净水体像元在冻结/融化的变化过程中会导致高达80—100 K 的亮温变化，但陆地成分的存在会导致湖泊未结冰时的亮温上升、湖泊结冰时的亮温下降，无冰/有冰的亮温差异显著减小。而MTT 算法主要依靠湖泊冻融时产生的亮温突变，因此混合像元会直接影响算法判断的结果。文中设置了两个像元大小（6.25 km）的缓冲区剔除靠近湖岸的混合像元，但是缓冲区的设置会导致被动微波数据的结果只包含了湖泊中心部分的信息。由于湖泊的冻融通常从湖岸开始，与实际冰情相比，被动微波对湖冰变化开始的信号响应可能有延迟，同时会提前认为湖泊完成变化。反映在湖泊冻融日期的结果中，被动微波的开始冻结和开始消融日期可能比实际情况晚，而完全冻结和完全消融日期可能比实际情况早。

4.4 数据分析方法

采用相关系数（r）、平均偏差（bias）和平均绝对误差MAE（Mean Absolute Error）评价被动微波和MODIS 数据获取的湖冰物候结果的一致性和差异性。用Mann-Kendall 检验评估湖冰物候变化趋势的方向和显著性，用Sen’s 斜率法计算湖冰冻融日期和部分气候要素时间序列的变化速率。同时相关系数也被用于评价湖冰物候与气候要素（气温、降水和风速）之间的相关性。

5 结 果

5.1 湖泊冻融日期提取结果的验证

无论是开始冻结日期、完全冻结日期、开始消融日期还是完全消融日期，被动微波和MODIS提取的结果都有较大的相关系数（图3）。整体上，被动微波获取的冻结日期比MODIS 获取的结果稍早，而消融日期更晚。冻结日期之间的偏差比消融日期小，但消融日期之间的相关性更高。其中，开始冻结日之间的偏差为3.95 d，MAE 为5.95 d（图3（a））。兴凯湖在冻结过程中容易受到云雨天气影响，持续的云覆盖可能导致了MODIS 获取的开始冻结日比实际情况有一定延后。完全冻结日之间的偏差为1.21 d，MAE 为3.42 d，是4 个日期中偏差和MAE 最小的（图3（b）），被动微波由于空间分辨率较低，容易忽略细节信息，因此对完全冻结日的判断可能比实际情况稍早。开始消融日之间的偏差最大，为-7.53 d，MAE 为7.53 d（图3（c）），表明被动微波获取的开始消融日每年都比MODIS晚，对于被动微波数据，像元内只要有部分的湖冰覆盖，都可通过高亮温识别出来，但对小范围的消融可能并不敏感，此外，被动微波的轨道条带间隔也可能进一步导致了开始消融日的延后。完全消融日之间的偏差为-5.89 d，MAE 为6.00 d（图3（d）），这同样可能是由于MODIS 易受云雨天气影响，宽松的阈值排除了噪声像元，也排除了小范围未消融的湖冰，导致判断的完全消融日提前，而被动微波对小范围的湖冰更敏感，轨道条带间隔也可能导致判断的完全消融日延后。

由于MODIS 和被动微波数据类型不同、分辨率不同、提取湖泊冻融日期的方法不同，两者获取的日期结果有不可避免的差异。MODIS 的空间分辨率更高，比被动微波更能反映实际情况，但云的存在增加了数据结果的不确定性；被动微波容易忽略湖面冻融的细节信息，且有周期性的数据缺失，但提取湖泊冻融日期的过程没有额外的干扰信息。尽管有一定的差异，总体上，被动微波和MODIS 获取的湖泊冻融日期都呈现显著的相关性，表明基于被动微波数据提取湖泊冻融日期的方法是可行的，其结果也是可靠的。

5.2 兴凯湖的湖冰物候变化

1979年—2019年，平均而言，兴凯湖大约每年11-13开始冻结，11-23完全冻结，湖冰冻结持续时间9.80 d（图4）。兴凯湖最早的开始冻结时间发生在1985年，是10-31，最晚在1979年，是11-23。最早的完全冻结时间发生在1980年，是11-10，最晚在2004年和2018年，都是12-03。一般而言，次年04-23 兴凯湖湖冰开始消融，04-30 湖冰完全消融，消融持续时间8.03 d。最早的开始消融日期在2019年，是03-27，最晚开始消融日期是1984年和2013年的05-05。最早的湖冰完全消融时间发生在2019年，04-09，最晚在1984年，05-15。

图4 1979年—2019年兴凯湖湖冰物候的变化Fig.4 The ice phenology variations of Khanka Lake from 1979 to 2019

过去41年，平均而言，兴凯湖湖冰完全封冻时间150.50 d，湖冰覆盖时间168.03 d。最短的完全封冻时间发生在2018年/2019年，具体日期为2018-12-03—2019-03-27，完全封冻时间为114 d。最长的完全封冻时间发生在1980年/1981年，也即1980-11-10—1981-04-27，完全封冻时间长达168 d。最短的湖冰覆盖时间也发生在2018年—2019年，具体日期为2018-11-21—2019-04-09，湖冰覆盖时长为139 d。最长的湖冰覆盖时间发生在1984年—1985年，也即1984-11-07—1985-05-14，湖冰覆盖时间超过半年，长达188 d。

1979年—2019年，兴凯湖开始冻结日期没有明显的变化，完全冻结日期平均推后了0.19 d/a，总共推迟了7.79 d（图4（a））。开始消融日期提前了0.16 d/a，完全消融日期提前了0.13 d/a，过去41 a开始消融日期和完全消融日期分别提前了6.56 d和5.33 d（图4（b））。兴凯湖湖冰完全封冻时间平均缩短0.31 d/a，过去41 a 总共缩短了12.71 d，同样，湖冰覆盖时间也是缩短的，速率为0.07 d/a，总共缩短了2.87 d（图4（c））。尽管本文使用的遥感数据及其时间序列不同，但是兴凯湖湖冰持续时间缩短的趋势与Yang 等（2019）的结果是一致的。

6 讨 论

将1979年—2019年兴凯湖的湖冰冻融日期变化与鸡西气象站的年总降水、年均风速和气温变化进行比较（表1）。可以发现湖冰冻融日期、湖冰持续时间与降水都没有显著的相关关系，表明降水变化对兴凯湖的湖冰物候变化几乎没有影响。风速与完全冻结日期呈显著的正相关关系，表明风速增大，湖冰扩散破碎，不易聚集，导致湖泊完全冻结日期的延后。除完全冻结日期外，风速与湖冰封冻持续时间也有负相关关系。冻结过程中，风的动力作用会机械地破坏较薄的湖冰，从而推迟了完整湖冰覆盖的形成，进而导致封冻持续时间的缩短。

表1 1979年—2019年兴凯湖湖冰物候变化与气候因素的相关关系Table 1 Correlations among ice phenology of Khanka Lake and climate factors from 1979 to 2019

鸡西气象站的数据表明1979年—2019年兴凯湖流域气温明显升高，升高速率为0.024 ℃/a。气温与冻结日期有一定的正相关关系，与消融日期和湖冰持续时间呈显著的负相关关系（表1，图5），是影响兴凯湖湖冰物候变化的主要原因。与消融日期相比，湖泊冻结受到的影响因素更多，除温度、风力作用外，还取决于湖泊形态和大小等特征。湖冰的消融主要取决于吸收的热量，升高的温度不仅会直接影响湖冰持续时间的长短，秋冬季节温和的气候也会造成冰厚减薄、湖面积雪减少，从而导致湖冰更易消融。在平均气温低的年份湖冰持续时间相应较长，如1981年、1982年、1984年、1985年、2010年，湖冰完全封冻时间和湖冰覆盖时间分别都超过了160 d和180 d。而平均气温高的年份湖冰持续时间相对较短，如2014年、2019年，湖冰完全封冻时间和湖冰覆盖时间分别都少于130 d和150 d。

图5 1979年—2019年兴凯湖湖冰覆盖时间与气温之间的关系Fig.5 The relationship between ice cover duration of Khanka Lake and air temperature from 1979 to 2019

7 结 论

利用1979年—2019年的被动微波遥感数据获取了兴凯湖的湖冰物候，与2000年—2019年MODIS 遥感数据的结果具有较好的一致性，这表明用被动微波亮温数据获取湖冰物候的方法是可行的，结果也是可靠的。过去41 a，兴凯湖湖冰开始冻结日期没有明显变化，完全冻结日期推迟了7.79 d，开始消融日期和完全消融日期分别提前了6.56 d 和5.33 d。湖冰完全封冻时间和湖冰覆盖时间分别缩短了12.71 d 和2.87 d。气温升高是湖冰冻结日期推后、消融日期提前以及湖冰持续时间缩短的主要原因。因此湖冰冻结、消融日期以及湖冰持续时间可以作为湖区气候变化的指示器。

湖冰作为气候变化的一种指示器，仅靠地面观测无法满足全方位、多视角湖冰研究的需求。对于兴凯湖这样偏远、跨国界、大面积的湖泊来说，用遥感方法来获取湖泊冰情信息进而分析其对区域气候变化的响应非常有必要。长时间序列、高时间分辨率、免费的被动微波遥感数据在湖泊冰情对气候变暖的响应研究方面发挥了举足轻重的作用。然而被动微波遥感数据空间分辨率低，大大限制了数据的应用范围。在今后的研究中，混合像元的分解可以作为新的研究方向，如能解决这个问题则可以在一定程度上扩大数据的应用范围，分析一些中小湖泊的湖冰物候变化。对亮温曲线不明显的年份，湖冰物候的提取方法也需要进一步优化和调整，寻找一种能够降噪的算法以提高湖冰物候提取结果的精度。此外，湖泊冰情变化的影响因素众多，除了气温之外，还有地形、风向风速、盐度、水深等因素，将高分辨率的多源遥感与气象、水文等观测或模式数据相结合，综合分析湖泊冻结消融状况及其对气候变化的响应也是未来重要的发展方向。