20世纪90年代以来，冰下湖曾被认为是彼此孤立、相对稳定的封闭系统，分布于分冰岭附近^[1-3]。2005年，Gray等的研究表明，这种观点需要修正^[4]，且实际上冰下湖水会发生流动，并使冰面高程产生变化^{[5, 6]}。卫星测高技术在冰下湖的活动监测中发挥了重要的作用。2008年，Smith等利用卫星测高发现124个活动的冰下湖，且主要分布于冰川流域的上游以及快速流动的冰流分支等区域^[7]。众多学者利用卫星测高数据研究表明，冰下湖水会向其他冰下湖流动，形成相互联系的冰下水文系统^[8-10]。冰下湖的活动对冰下水文、冰流速、冰盖物质平衡以及海平面变化等有着重要的影响^[11]。

由于人们对冰下湖活动的了解尚少，为了分析其对冰流速和冰盖物质平衡等的具体影响，需要确定冰下湖的位置、水量变化、活动周期以及冰下水文关系等。早期的雷达测高卫星地面脚斑直径达千米级，在冰下湖的应用中受到限制^[12]。2003年发射的冰、云和陆地高程卫星(ICESat)激光高度计有较小的地面脚斑直径(约为70 m)，提供的高精度高程数据集被广泛应用于冰下湖的监测。2010年发射的CryoSat-2卫星，其全新的SARIn模式测高数据在冰下湖的监测中同样发挥着重要作用^[13-15]。

Totten冰川流域是东南极最大的冰川流域，冰盖底部消融导致Totten冰川长期处于负物质平衡状态，在全南极物质平衡中占有重要的作用^{[16, 17]}。冰下湖是冰盖底部独特的组成部分，由于Totten冰川流域不稳定的冰盖底部状态，研究该区域冰下湖活动变得尤为重要。本文选取Totten冰川流域的冰下湖，利用ICESat激光测高数据，采用重复轨迹最小二乘平面拟合的方法，计算冰下湖的高程和水量变化，并结合BEDMAP 2数据生成的冰下水流路径图分析冰下湖的水文关系。

1 实验区域和数据

实验区域为位于东南极Totten冰川流域的3个冰下湖，即Totten₁、Totten₂和Wilkes₁冰下湖(见图 1)。图 1中湖的边界是通过内插获取的冰面高程变化值为零的等值线，底图是2004年MODIS影像镶嵌图^[18]，灰色直线表示经过该区域全部的ICESat重复轨迹，方框表示本文中冰下湖的研究范围，每个冰下湖表面均有2~3条ICESat轨迹覆盖。

本文使用的研究数据为美国国家冰雪数据中心(NSIDC)提供的ICESat/GLA12 L2级数据产品，时间为2003年10月~2009年4月。为了保证数据质量，本文首先对由大气、地表反射和仪器等造成的误差根据GLA12的数据质量标准做了数据筛选，去除不合格的数据，再对数据进行饱和度改正^[19]，并利用文献[20]中提供的偏差值对ICESat任务间的高程偏差进行改正，最后将测高数据统一到WGS84坐标系上。为了分析冰下湖水的流动方向，利用BEDMAP 2数据集中的冰底地形和冰盖表面地形生成冰下水流路径图。

2 数据处理方法 2.1 ICESat测高数据监测冰下湖

冰下湖的活动会引起冰盖表面高程发生变化，一般将冰盖剖面高程变化量具有“抛物线”型的特征视为冰下湖的活动。卫星测高重复轨道算法常用于冰下湖的活动监测。本文将每组ICESat重复轨迹测高数据分割成沿轨迹方向长700 m、垂直于轨迹方向宽350 m的矩形条带。为了保证计算的连续性，每个矩形条带沿轨迹方向有200 m的重叠。将每个矩形条带内的测高数据利用最小二乘方法拟合成一个随时间变化的平面(式(1))，然后计算出冰盖的高程变化率k和消除地形影响的高程随时间的变化量dh(式(2))。

$ \begin{array}{l} H = D\left( {x, y} \right) + \Delta h\left( t \right) = \\ {m_1}\left( {x - \overline x } \right) + {m_2}\left( {y - \overline y } \right) + {m_3} + kt \end{array} $

(1)

$ \begin{array}{l} {\rm{d}}h = {H_p} - D\left( {x,y} \right) = \\ {H_p} - {m_1}\left( {x - \overline x } \right) + {m_2}\left( {y - \overline y } \right) + {m_3} \end{array} $

(2)

式中，H表示拟合的高程；D(x, y)表示某一时刻的基础地形项；Δh(t)表示高程随时间变化项；(x, y)表示脚点坐标；(x, y)为矩形条带内点坐标的平均值；t是测高数据的获取时间；m₁、m₂、m₃表示待求系数；k是高程变化率，且通常认为是恒定的；H_p表示ICESat各任务期间的测高数据脚点高程。为了保证观测的准确性，至少保证有2条重复轨迹穿过该冰下湖，每条重复轨迹不少于6个，矩形条带内的点不少于10个。

通过上述方法可以得到每条重复轨迹随时间的高程变化量。假设有n条重复轨迹穿过冰下湖，首先将重复轨迹k中高程变化量大于0.1 m的连续段作为冰下湖内的点，计算该段连续轨迹的长度l_k和每个ICESat任务期间的测高数据高程变化平均值${\rm{d}}\overline h \left( t \right)$，将${\rm{d}}\overline h \left( t \right)$减去冰下湖外高程变化平均值，以消除积雪的影响，最后得到由于冰下湖活动引起的平均高程变化Δ ${\overline h _k}\left( t \right)$，利用轨迹长度加权平均的方法计算出冰下湖平均高程变化Δh(t)。对高程变化值内插，将高程变化值为零的等值线作为冰下湖的边界^[7]，从而确定冰下湖的面积A，即可得到冰下湖近似的平均水量变化ΔV(t)：

$ \Delta V\left( t \right) = \Delta \overline h \left( t \right)A = \frac{{\sum\limits_{k = 1}^n {\Delta {{\overline h }_k}\left( t \right){l_k}} }}{{\sum\limits_{k = 1}^n {{l_k}} }} $

(3)

2.2 水流路径获取

冰盖底部的水压势能是冰盖底部水的重力势能与冰下水压的和，冰盖底部水流方向受到水压势能梯度的影响。由于冰下水压近似等于冰层的负重压力，因此水压势能方程近似为^[21]：

$ \mathit{\Phi} = {\rho _i}g{z_s} + \left( {{\rho _w}g - {p_i}g} \right)h $

(4)

式中，Φ表示水压势能；ρ_i和ρ_w分别为冰和水的密度，取值分别为917 kg/m³和1 000 kg/m³；g为重力加速度；z_s为冰盖表面高程；h表示冰下基岩高程。分析发现，由于冰盖巨大压力的影响，使得冰盖的影响是冰底地形的10倍^[22](ρ_i/(ρ_w-ρ_i)≈10)。基于上述方程，利用BEDMAP 2数据集中冰盖表面及冰底地形数据，获取冰下水压势能表面栅格，最后利用水文分析模型得到冰下水流路径图。

3 实验结果和分析 3.1 冰下湖高程变化率分析

本文首先计算了实验区域所有ICESat重复轨迹的高程变化率，分析发现冰下湖区域的高程变化率具有一致性，通过此规律可确定3个冰下湖的基本位置和经过每个冰下湖的ICESat测高数据的重复轨迹编号。如图 2所示，覆盖Totten₁的轨迹有Track107和Track145，覆盖Totten₂的轨迹有Track11、Track211和Track130，覆盖Wilkes₁的轨迹有Track145和Track360。

3.2 冰下湖活动分析

确定各冰下湖的大致位置和经过表面的重复轨迹后，对每条重复轨迹利用§2.1中的方法进行高程变化量解算，并内插得到重复轨迹在湖内的长度，进而估算冰下湖水量变化。

图 3(a)是经过Totten₁冰下湖的Track107和Track145高程变化结果, 分别表示两条重复轨迹在ICESat不同任务期间的高程变化量，竖直的红色虚线表示冰下湖边界的位置，l₁和l₂为湖内各重复轨迹的长度，经内插计算得到的长度约为28 km和42 km。由图 3(a)可见，Track107和Track145两条轨迹随时间的高程变化趋势基本是一致的。在2003年10月到2005年6月期间，由于水量减少导致冰下湖表面产生最大约3.8 m的下降；随后水量逐渐增加，使冰盖表面在2005年6月到2006年11月缓慢上升，最大约为0.5 m；在2006年11月到2007年10月，冰下湖水大量增加，导致冰盖表面产生约2 m的上升；之后湖水继续缓慢增加，导致冰盖表面最大上升约0.5 m。

覆盖Totten₂冰下湖的Track11、Track211和Track130在冰下湖内的轨迹长度分别约为21 km、36.5 km和29 km。分析3条轨迹的高程变化量(见图 3(b)，其中Track11未列出)发现，Totten₂冰下湖的活动较频繁，其中在2005年6月，由于水量减少导致冰盖表面最大下降约1.6 m。此外，3条轨迹的高程变化趋势在时间上并不一致，原因是水向冰下湖不同区域流动时，由于冰下地势及阻滞等作用，导致各区域高程变化时间存在偏差^[10]。

对于经过Wilkes₁冰下湖的Track145和Track360，其湖内的轨迹长度分别为31 km和28 km。分析图 3(c)发现，冰下湖在该段时间水量持续增加，导致冰盖表面最大上升约2.4 m。此外，图 3(c)中冰下湖外一侧的高程变化也比较明显，表明该区域可能存在较明显的冰盖底部消融现象^[17]。

考虑到冰下地势差异和阻滞作用导致不同ICESat重复轨迹的高程变化趋势在时间上存在偏差，本文采用加权平均计算每个冰下湖的平均高程变化(图 4(a))，再根据3个冰下湖的面积计算水量变化(图 4(b))。经计算，Totten₁、Totten₂和Wilkes₁的面积分别为567 km²、710 km²和596 km²。图 4中3个冰下湖高程和体积相对变化的起始时间是2003年10月。Totten₁冰下湖在2003年10月到2005年6月，水量减少约为1 km³；2005年6月到2006年11月，水量基本维持在稳定状态；在2006年11月后，约有0.6 km³的水量增加。Totten₂冰下湖的水量变化分别有两个最大值和最小值，在2005年3月达到最大，减少量约为0.69 km³，到2006年11月又有约等量水的增加，从后续变化推测Totten₂冰下湖活动呈周期性的变化。对于Wilkes₁冰下湖，在观测时间段内约有0.54 km³的水量增加。

3.3 冰下湖水文关系分析

为确定3个冰下湖的冰下水文关系，本文利用BEDMAP 2数据集获取该区域冰下水流路径图。如图 5所示，底图表示水压势能，灰色线表示水压势能等值线，蓝色线表示水流路径。由图 5可以看出，冰下水的主要流动方向是由西北方流向东南方。3个冰下湖在水流路径图中的水文分布关系为Totten₂位于Totten₁的上游，Totten₁位于Wilkes₁的上游，它们之间具有一定的关联性。由于冰下湖的活动受到上游水流、湖的容水量及下游水流等因素的影响，从而形成了图 4中的变化特征。Totten₁和Totten₂冰下湖水流入和流出较频繁，其中Totten₂冰下湖位于最上游，其变化特征表现出一定的周期性；Totten₁冰下湖位于Totten₂冰下湖的下游，湖水在该时段内整体表现出先减少后增加的现象；而3个湖中Wilkes₁冰下湖位于最下游，在本观测时间内呈现水量持续增加的现象。

4 结束语

本文利用ICESat卫星激光测高数据，采用重复轨道最小二乘平面拟合的方法，分别获取了东南极Totten冰川流域Totten₁、Totten₂和Wilkes₁三个冰下湖2003-2009年的高程和水量变化，并结合BEDMAP 2数据集生成的冰下水流路径图分析这3个冰下湖的水文关系。实验结果表明，Totten₁冰下湖在观测时间内水量先逐渐减少，经过一段稳定期后水量逐渐增加；Totten₂冰下湖的活动存在一定的周期性；Wilkes₁冰下湖在观测时间内水量持续增加。水流路径图显示3个冰下湖具有一定的关联性。

ICESat卫星测高数据能有效应用于冰下湖的活动监测，并且能辅助分析冰下水文关系，然而其时间和空间覆盖范围均有限。今后将利用CryoSat-2及其后续卫星测高数据^[23]对冰下湖进行监测，分析冰下湖的活动规律及其对冰下水文和冰盖运动等的影响。

致谢: 感谢美国国家冰雪数据中心(NSIDC)提供的MODIS影像镶嵌图和ICESat测高数据，以及英国南极局提供的BEDMAP 2数据。