青藏高原纳木错湖冰物候变化遥感监测与模拟

吴艳红，郭立男，范兰馨，文梦宣，迟皓婧，张兵

1.中国科学院空天信息创新研究院 数字地球重点实验室,北京100094;

2.中国科学院大学,北京100049

1 引 言

湖冰物候（湖泊的冻结和消融时间）直接受湖泊能量平衡变化的影响，对区域气候变化有很好的指示作用（Kouraev等，2007；Marszelewski和Skowron，2006）。湖冰物候在反映气候变化方面某种程度上可能比气温更可靠（Livingstone，1997）。湖冰物候的变化影响着湖泊生态系统中生物的繁殖以及渔业、运输业等人类生产活动。因此，近年来中高纬度湖冰物候对全球变暖的响应越来越为研究人员所关注（Salonen等，2009；Bengtsson，2011；Fu和Yao，2015；Gou等，2017；Latifovic和Pouliot，2007；Magnuson 等，2000）。比如，Magnuson 等（2000）利用实测数据分析了150年间北半球39 个湖泊的湖冰物候变化趋势，发现北半球湖冰物候表现出结冰日延后和融冰日提前的趋势，并为后续的研究所证实（Bernhardt 等，2012；Brown 和Duguay，2010；Ghanbari等，2009）。

基于AVHRR、MODIS 等光学传感器（如：Latifovic和Pouliot，2007；Kropáček等，2013；Duguay等，2015）和微波遥感器（如：QuickSCAT，SSM/I，SMMR，AMSR-E）等对湖冰物候进行动态监测已成为当前湖冰物候研究的重要手段。其中，光学遥感湖冰物候监测常用的方法包括阈值法、指数法、反射率光谱曲线法等（庞毓雯等，2020；Guo等，2018）。在微波遥感方面，邱玉宝等（2017）利用微波亮温数据分析了青藏高原湖冰的冻融变化。

遥感是获取湖冰物候信息的有效途径。然而，由于遥感平台的服务年限相对较短，难以揭示湖冰物候的长期变化规律。另外，遥感信息反映了湖冰物候变化的事实而无法解释其变化的物理机制。认识湖泊冻融过程的物理机制，进而量化湖冰物候对气候的响应，通常则需借助数值模拟的方法来实现（任晓倩等，2014）。数理模型考虑了多种气候强迫项和湖泊形态特征对湖冰冻融过程的作用，有助于从机理的层面揭示和预测湖冰的生消过程。现有湖泊过程模式如MIMLAKE96（Fang 和Stefan，1996）、LAKEoneD （Joehnk 和Umlauf，2001）、LAKE模式（Stepanenko等，2016）、SIMSTRAT （Goudsmit 等，2002） CLIM （Duguay等，2003）、LIMNOS（Vavrus等，1996）、HIGHTSI（Launiainen 和Cheng，1998）等。这些模型在世界各地区已经得到了应用和验证。中国湖泊研究起步较晚，现有模型多为外国学者针对当地环境开发，在青藏高原的适用性仍然有待验证（宋兴宇等，2020）。方楠等（2017）改进并评估了WRF模型在纳木错的适用性；苏东生等（2018）利用“中国区域高时空分辨率地表气象驱动数据集（CMFD）”驱动Flake 模型模拟青海湖热力状况对气候变化的响应；Huang 等（2019）评价了3 个一维湖泊模式在纳木错湖泊的适用性并通过敏感性试验对模型进行率定。

青藏高原地区受人类活动影响较小，其气温升温速率大于同期全球陆表气温的变化速率（Huang等，2017），是全球变化的敏感区域（杜娟等，2020）。由于历史观测数据的缺乏，青藏高原湖冰物候的变化研究还不够充分。本文利用过程模型（LAKE 2.3）重建青藏高原纳木错湖泊长时序（1962年—2018年）的湖面温度和湖冰物候，并将模拟结果与遥感监测结果进行对比验证。在此基础上，分析近60 a 来纳木错湖冰物候的变化趋势。本文所用的LAKE 2.3 模型（Stepanenko 等，2016）是综合性的一维湖泊模型，涉及湖泊的热力学过程以及生物地球化学过程。该模型已有许多应用，但在青藏高原湖泊研究中的应用还较少。

2 研究区与数据

2.1 研究区概况

纳木错位于30°30′N—30°55′N，90°16′E—91°03′E（图1），海拔4718 m，湖泊面积1981 km2。纳木错是西藏自治区第二大湖，也是世界上海拔最高的湖泊之一。湖泊集水面积10610 km2，2005年测得的最大深度超过90 m（Wang 等，2009）。纳木错所处流域内冰川占比约1.5%（Wu 等，2014）。近几十年来，随着青藏高原气温的升高，纳木错湖泊的面积、水位、水温、冰期以及生物化学过程也发生着深刻的变化（Guo 等，2020；Huang 等，2017；Zhang等，2011）。纳木错属于高原亚寒带季风半干旱气候带，辐射强烈，年日照时数达2900—3200 h，年均温1.9 ℃，年降水量约500 mm，湖面年蒸发量约1200 mm。

图1 纳木错及邻近气象站点地理位置Fig.1 Locations of the Nam Co and the nearest meteorological station

2.2 遥感数据

本文首先基于多源遥感信息应用不同方法提取纳木错的湖冰物候信息。主要的遥感数据源包括：ARC-Lake（V3），MODIS 以及基于微波辐射计（AMSR-E、AMSR2 和FY-3 MWRI）的湖泊亮温信息。

ARC-Lake（V3）是基于沿轨道扫描辐射计（ATSR2/AATSR）获取的湖泊表面温度日尺度数据集（http：//www.laketemp.net/home_ARCLake/［2021-04-30］），包括日间和夜间的湖面温度。本研究以日间温度和夜间温度的均值作为纳木错湖面温度日值，以此获取的湖面温度序列作为模拟湖面温度的验证数据之一，时间范围1995-06-01—2012-04-05。

中分辨率成像光谱仪MODIS（moderate resolution imaging spectroradiometer）搭载于两颗极轨卫星Terra和Aqua上。Terra上午自北向南穿过赤道，Aqua下午穿过赤道。双星系统1—2 d 就可以覆盖全球一次。MODIS 传感器包含36 个波段，已经拥有应用于各个领域的多种产品。本文应用的是其L3级日尺度陆面温度产品（MOD11A1），空间分辨率为1 km。数据集包含日间和夜间的陆表温度及发射率。为计算湖泊面上的平均温度，本文首先将纳木错湖泊矢量边界图向内缩进1 km 以避免湖岸混合像元的影响。然后，以日间温度和夜间温度的均值作为湖表的日均温度，从而获取湖表日均温度序列。为保证数据集之间的可比性，本文选择的MODIS产品的时间跨度与ARC-Lake一致。

此外，本文还采用了邱玉宝等（2017）生产的“2002年—2016年高亚洲地区中大型湖泊微波亮温和冻融数据集”（http：//www.csdata.org/p/49/［2021-04-30］）。该数据集依据微波辐射计（AMSR-E、AMSR2和FY-3 MWRI）观测像元内湖泊和陆表的面积比例，并应用混合像元分解方法获取像元（亚像元级）的湖泊亮温信息，是通过被动微波遥感对像元级湖冰冻融进行监测的一种方法。

2.3 LAKE模型驱动数据

LAKE 模型需要的输入包括气温、比湿、气压、风速、短波辐射、长波辐射和降水等气象要素。短波辐射和向下的长波辐射根据纬度、实际日照时数、气温和水汽压以计算而得。其他气象变量以距纳木错最近的当雄站（30°29′N，91°06′E，海拔4200 m）的气象数据为依据。基础气象数据来自中国气象数据网，包括1963年—2019年近60 a的气温、降水、风速、气压、相对湿度和日照时数。

考虑到气温随高程的变化，本文首先以0.65℃/100 m 的气温大气直减率（胡芩 等，2014）对气象站气温Tstation进行高程校正得到湖泊气温Tlake，即

式中，Hstation为气象站海拔，Hlake为湖泊海拔。其次，利用气压阶公式对气象站气压Pstation进行高程校正，获取湖面气压Plake，即

式中，比湿q的计算公式为

式中，ε= 0.622，e为水汽压，是相对湿度RH 与饱和水汽压ew的乘积，而饱和水汽压可用Bolton公式（Bolton，1980）计算：

3 方 法

3.1 湖冰物候指标

本文基于遥感与模拟提取的湖冰物候指标包括：开始结冰日FUS（湖面由完全是水到开始有像元变为冰对应的日期）、完全融化日BUE（湖面由存在冰到所有像元为水对应的日期），以及湖泊结冰期（ID），即BUE与FUS二者之间的天数。这3个指标与湖气界面的水热交换过程密切相关。

3.2 湖冰物候遥感监测

本文参照Guo等（2018）的方法，基于不同的遥感信息源提取纳木错2000年—2015年的湖冰物候，并分析其变化规律。湖冰物候的遥感监测方法具体包括：（1）利用MOD09GA 产品生成全湖平均的反射率时间序列，通过阈值法在反射率时间序列上提取湖冰物候（M1）；（2）利用MOD10A1冰雪产品，计算湖面水体占湖泊面积的比，以水占 比10% 为 物 候 提 取 阈 值（M2）；（3） 利 用MOD11A1 产品生成全湖平均湖面温度时间序列，通过阈值法提取湖冰物候（M3）。上述3 种方法的详细过程可参考Guo 等（2018）。此外，本文从“2002年—2016年高亚洲地区中大型湖泊微波亮温和冻融数据集”中获取纳木错2002年—2016年的湖冰物候作为对比数据（M4）。

3.3 基于LAKE模型的湖冰物候模拟

本文采用的LAKE模型（Stepanenko等，2016）其热力学过程表述为

式中，T为水温，垂直坐标z以湖面为z=0，垂直向下为正，h为湖泊深度，并定义为ξ=z/h为坐标。cw是水的比热容，ρ是水的密度，λ涡流扩散率，是湖泊表面水量平衡，S是太阳辐射通量，穿透到深度z，M是浮力混合（对流）。穿透到深度z的太阳辐射计算：

式中，αe是消光系数。

湖泊水温的上边界条件是湖泊表面能量平衡：

式中，S是太阳辐射，Ea向下的大气辐射，Es向上的长波辐射，H和LE分别是显热和潜热通量，α是表面反照率。

本研究中用到的是LAKE 2.3 模型。LAKE 模型包含的参数较多，但参数大多具有明确的物理意义。模拟过程中，除了湖泊基本形态参数（湖泊经纬度、湖泊面积、湖泊深度等），本文参考2005年—2007年实地考察的纳木错湖盆形态估算出各个深度的水平截面面积，作为湖泊形态参数输入模型。在模拟中，本文取纳木错的平均深度为40 m，水体的消光系数为0.1，冰的消光系数为100。主要的模型参数设定如表1。这些参数取值与Huang等（2019）的研究中对纳木错消光系数设置相一致。

表1 纳木错LAKE模型主要参数设定Table 1 Model parameter configuration for the Nam Co

本文设定LAKE 模型时预热时间为20 a，并以模型最终输出的表层温度（skin temperature）的年内变化特征获取湖泊结冰和融化的日期。

4 结果分析与讨论

4.1 基于遥感的湖冰物候

图2 给出了纳木错2003年—2014年基于不同遥感监测方法的湖冰物候。纵轴是年积日（DOY）（从当年01-01 起开始计算的天数），大于365 d 的表明物候日期发生在第2年。由图2 中可以看出，纳木错湖泊结冰日多始于12 月底至次年1 月，次年4—5 月湖冰则完全融化。多年平均情况下，开始结冰日（FUS）为01-07，完全化冻日（BUS）为03-10。年均湖泊结冰期（ID）的长度大致为124 d，标准差为10 d。研究时段内，湖冰物候的变化总体趋势并不显著。

图2 遥感监测的湖冰物候（每个柱状图的上下端分别对应湖冰物候中的FUS和BUE；M1—M4分别为对应的4种遥感监测方法）Fig.2 Lake ice phenology from remote sensing（The upper and lower ends of each histogram correspond to FUS and BUE in lake ice phenology respectively；M1—M4 represents the four different methods used）

不同方法之间的比较可以看出，基于微波亮温数据（M4）获取的冰期通常比其他方法监测所得的冰期短，而其他3种遥感方法在湖冰物候监测上的差异没有明显的规律。这一结果表明，受遥感数据时空分辨率的影响（如：数据缺失），湖冰物候的遥感监测还存在较大的不确定性。

4.2 湖冰物候模拟

由于遥感的数据序列相对较短，难以反映纳木错湖冰物候近几十年来的变化规律，因此本文借助LAKE模型以重建纳木错的湖冰物候序列，并揭示湖冰物候的长期变化规律及其对气候变化的响应。

图3 首先给出了LAKE 模型湖面温度的模拟结果。由图3中可以看出，湖泊表面温度的模拟值与遥感监测的湖面温度一致性较高，模拟结果与ARC-Lake 水温数据的纳什效率系数（Nash-Sutcliffe efficiency coefficient）达0.82，均方根误差为1.75℃。基于LAKE 模型的湖冰物候模拟结果如图4 所示。由图4 可以看出，模拟所得的3 个物候指标中，融化日期BUE 与遥感监测的结果一致性最高，均方根误差为8.5 d。模拟的开始结冰日FUS普遍早于遥感监测的结果，模拟的结冰期总体上大于遥感监测的结冰期。

图3 基于模型和遥感的湖面温度对比Fig.3 Comparison of lake surface temperature based on modelling and remote sensing

图4 模拟湖冰物候与遥感湖冰物候（每个柱状图的上下端分别对应湖冰物候）Fig.4 Comparison of lake ice phenology based on modelling and remote sensing（The upper and lower ends of each histogram correspond to FUS and BUE in lake ice phenology respectively）

LAKE 模拟结果与遥感监测结果之间的差异，一方面受模型的不确定性的影响（包括模型结果、参数以及输入的不确定性）；另一方面也受遥感监测的不确定性的影响。特别是遥感监测的湖面温度通常为卫星过境的瞬时值（或日间夜间瞬时值的均值），而不能真正代表一天内的平均温度，且MODIS 监测的表面温度受云等影响波动较大。此外，确定湖冰物候时相关阈值的设定也会影响湖冰物候的判别。因此，虽然遥感和模拟都可以获取湖冰物候信息，但湖冰物候监测与模拟的准确性还需要进一步的地面观测加以验证。

4.3 近60年来湖冰物候变化趋势

LAKE 等有物理基础的数值模型刻画了气候因子与湖冰物候之间的定量关系，有助于揭示湖冰物候对气候变化的响应。基于LAKE模型，本文重建了近60 a来纳木错的湖冰物候如图5所示。

图5 1962年—2018年纳木错湖冰物候变化趋势（FUS、BUE、ID分别为开始结冰日、完全融化日和冰期）Fig.5 Trends in lake ice phenology during the period 1962—2018 of Nam Co（FUS，BUE，and ID are date of freeze-up start，break-up end and ice duration respectively）

基于Mann-Kendall（MK）非参数检验法（Kendall，1955；Mann，1945）的湖冰物候变化趋势分析表明，近60 a 来纳木错湖泊结冰推迟、融化提前、冰期缩短，且所有物候指标的变化趋势均达到95%的显著性水平。结冰日期（FUS）平均每10 a延后3.3 d，融化日期（BUE）每10 a提前3.0 d。冰期（ID）的变化速率最快，以6.4 d/（10 a）的速率缩短。

4.4 湖冰物候对气候变化的敏感性

已有研究表明，至21世纪末，青藏高原气温相对于1980年—1999年将升高2℃（Chen和Frauenfeld，2014）。为了定量分析湖冰物候对气候变化的响应特征，本文基于LAKE模型开展了湖冰物候对气候变化的敏感性分析。在敏感性分析中，假设纳木错气温以0.1℃的间隔持续升高至2℃，而后基于LAKE 模型的模拟结果，计算分析湖冰物候随气温升高相较于历史时期（1963年—2018年）的可能变化，结果见图6。

图6 湖冰物候对气温变化的敏感性（FUS、BUE、ID分别为开始结冰日、完全融化日和冰期）Fig.6 Sensitiveness of lake ice phenology in response to increasing air temperature（FUS，BUE，and ID are date of freeze-up start，break-up end and ice duration respectively）

模拟结果表明，当气温升高达0.7℃时，开始结冰日（FUS）显著推迟。如果气温升高2℃，结冰日期平均推迟7.3 d。湖冰融化过程对气温的响应比结冰过程更显著。当气温升高2℃，融化日期平均提前12.4 d。因此，冰期随着温度升高将显著缩短，平均缩短19.7 d。湖冰物候的变化会影响湖泊水量平衡、能量平衡过程、湖泊及周围的生态系统和人类活动。因此，在气候变暖的背景下，需要对湖冰物候变化情况进行持续的监测和预测。

5 结 论

本研究以纳木错湖泊为研究区，基于4种不同的遥感监测方式获取了2002年—2015年纳木错的湖冰物候指标，结果表明各方法之间的年际变化情况比较一致，但基于微波亮温数据获取的冰期通常比其他方法监测的时间短。

在遥感监测的基础上，利用一维湖泊模型（LAKE）重建了1962年—2018年纳木错湖冰物候序列，揭示了近60 a 来纳木错湖冰物候的变化趋势及其对气候变化的响应特征。研究结果表明，近60 a 来，纳木错结冰推迟、融化提前、冰期缩短（6.7 d/（10 a））。湖冰物候（特别是融冰日期）对气候变化十分敏感。在气温升高2℃的气候变化情景下，未来纳木错冰期将缩短19.7 d。

本研究将遥感与数值模拟相结合，发挥两种手段的优势，为湖冰物候研究提供了新的思路和方法。但是，在青藏高原湖冰物候的研究中，现有的遥感方法和数值模拟方法都存在不同程度的不确定性。在未来的研究中，除了需要进一步提高遥感信息的精度和分辨率，完善模型的结构和参数化方案外，在条件适宜的地方开展持续的湖冰物候地面观测，将有助于提高青藏高原湖冰物候监测与模拟水平。