基于耦合热传导-冰流模型的冰川钻孔温度重建的气候反演算法比较

孙 欢， 王宁练， 张 泉

（1. 陕西省地表系统与环境承载力重点实验室，陕西 西安 710127； 2. 西北大学 城市与环境学院地表系统与灾害研究院，陕西 西安 710127）

0 引言

对过去冰面温度变化过程的研究是气候重建和预测冰川对未来气候变化响应的重要问题之一［1］。由于热传导的作用，冰面温度变化信号向冰面以下传播，对稳态冰温垂直分布造成了相应的扰动。通过测量冰川钻孔温度可知扰动后的冰温-深度剖面，该剖面能够保存过去几个世纪甚至更长时间的冰面温度变化信息［2-3］。通过冰川钻孔温度方法重建气候变化历史就是建立在冰面温度变化和冰温垂直分布的直接联系之上；考虑到冰川运动的影响，可建立耦合的热传导-冰流物理模型来描述冰面温度变化信号向下传播的过程［4］。模型以稳态冰温垂直分布和冰面温度变化为初始条件和边界条件，其中稳态冰温垂直分布由冰川底部地热流和长期平均冰面温度共同确定。与求解模型中任意时刻的冰温-深度剖面这一正问题不同，通过测量冰川钻孔温度重建冰面温度变化是求解模型的边界条件，是一个典型的反问题。由于热扩散的影响，利用冰川钻孔温度重建气候变化历史具有较低的分辨率；虽然无法体现短期气候事件，却可以在一定时间尺度得到完整的冰面温度变化过程［5］。

近三十年来，冰川钻孔温度方法在两极和高纬度冷冰川地区的气候重建研究方面得到了广泛的应用［6-15］。这部分地区缺少早期的气象观测数据，冰川钻孔温度方法无疑提供了一个独立完善的研究工具，具有十分重要的意义。在实际应用中，受冰川钻孔深度、实地环境等方面的影响，应用冰川钻孔温度方法重建气候有着十分严格的适用条件。要求冰川处于相对稳定的状态，底部冻结在基岩上无滑动过程；冰面消融较少，不涉及径流、渗浸作用等热交换过程，即冰川内部相变对钻孔温度的影响较小；冰川钻孔位于分冰岭附近，水平方向的运动速度可忽略等。热传导-冰流模型的反演算法是该研究的基础和关键。以往对于热传导-冰流模型的反演算法研究存在不足，如：没有详细讨论该模型解的稳定性问题；对各类反演算法的优劣性和适用性也没有具体的比较和说明。本文首先通过构造冰面温度变化和模型参数，进行理想条件下的数值实验，分析了最小二乘法、Tikhonov 正则化方法和蒙特卡罗算法这三种反演算法重建百年尺度冰面温度变化的结果。然后，以我国藏北高原腹地的马兰冰川钻孔资料为例设定模型参数，结合气象观测数据，详细比较了上述反演算法在10 a 和40 a 尺度的重建结果。最后，对构造理想条件下和应用真实钻孔数据的两组数值实验结果进行了比较，并进一步讨论该反问题中解的稳定性问题和不同算法的优劣性，指出算法的适用条件。本项研究可为今后利用冰温-深度剖面进行气候重建研究提供借鉴。

1 模型描述

早在1955 年，Robin［16］研究了冰流运动对冰川钻孔温度的影响。Dansgaard 等［17］1969 年建立了格陵兰冰盖的运动模式。这些研究为利用冰川钻孔温度重建冰面温度变化奠定了基础。耦合的热传导-冰流模型如下［4］：

式中：T为冰川冰在某深度任意时刻的温度（℃）；t为时间（s）；v→为冰流运动速度（m·s-1）。ρ、c和λ分别为冰的密度（kg·m-3）、比热容（J·kg-1·℃-1）和导热率（W·m-1·℃-1）。为温度随时间的变化率；∇T为温度在空间各方向上的全微分，即温度的空间变化率。f为冰川内部与热传导过程无关的热源项（J·s-1·m-3），如：冰川内部由于复杂相变产生的热量等。

在冰川内部与热传导过程无关的热源项、冰流运动水平方向速度可忽略的条件下，式（1）可简化为如下的一维方程［18］：

且满足：

式中：z为竖直方向（m），向下为正；v(z)为冰川垂直于地表面方向的运动速度（m·s-1）。κ、H和tf分别为冰的热扩散率（m2·s-1）、冰川厚度（m）和重建时间（s）。T(z，t)是在t时刻深度z处的瞬时温度（℃），即由冰面温度变化μ(t)导致的温度扰动。

模型以稳态冰温垂直分布函数U(z)为初始条件。在初始时刻t= 0 的冰温-深度剖面扰动（瞬时温度）为零，即T(z，0) = 0。θ(z)对应测量时刻t=tf的冰温-深度剖面扰动，如图1 所示。可以看出，μ(t)和θ(z)大于零代表冰面温度较稳态冰面温度U(0)高，扰动为正；小于零则代表冰面温度较U(0)低，扰动为负。

图1 冰面温度变化对稳态冰温垂直分布的影响Fig. 1 Influence of ice surface temperature change on the vertical distribution of steady-state ice temperature

稳态冰温垂直分布函数U(z)由冰川底部地热流密度q（W·m-2）和长期平均冰面温度Us（℃）共同确定：

式中：z为竖直方向（m），向下为正；v(z)为冰川垂直于地表面方向的运动速度（m·s-1）；U(0)为稳态冰面温度。κ、λ和H分别为冰的热扩散率（m2·s-1）、冰的导热率（W·m-1·℃-1）和冰川厚度（m）。

冰川底部地热流密度q可看作是短期内恒定不变的。在假设冰川稳定的条件下，冰川厚度H不变；冰川垂直于地表面方向的运动速度v(z)与平均积累率相等；v(z)随深度的增加呈线性减小，至冰川底部减小为零［13］。基于热传导-冰流模型的气候重建研究就是应用相关的反演算法解决以下问题：已知式（2）在某一时刻的解T(z，tf) =θ(z)，求出边界条件μ(t)。

需要注意的是，冰面温度变化信号向下传播的过程中，由于热扩散的影响，其振幅随着深度的增加而指数衰减、位相延迟。要利用冰川钻孔温度重建气候，在较大深度的冰温测量就需要很高的精度，一般不低于0.01 ℃［4］。实际中要达到这一精度和可靠测定难度较大。通常冰川钻孔温度由热敏电阻温度计测量，即通过测定冰层电阻，利用电阻-温度转换关系确定待测深度的温度［19］。具体做法是：将按照一定长度间隔分布的若干个热敏电阻探头放置待测深度；待系统达到热平衡状态后，利用数字万用表逐个测量对应深度电阻传感器的电阻值；电阻值经计算程序转换为对应的温度值。事实证明，已有的这种测温系统在负温条件下的分辨率可达0.01 ℃，取得了很好的结果［19-20］。

另外，在实际中应用热传导-冰流模型有非常严格的适用条件，通常需要满足以下几点要求：（1）冰川底部呈负温，与基岩冻结在一起，如极大陆型冰川；（2）冰面消融较少，内部相变对钻孔温度的影响可忽略；（3）冰川钻孔位于分冰岭附近，水平方向的运动速度可忽略。由此可知，两极和高纬度冷冰川分冰岭附近的钻孔温度可以较好地用来重建气候变化过程。

2 反演算法

2.1 最小二乘法

最小二乘法（最小平方法）通过最小化误差的平方和寻找未知数据的最佳匹配，使得计算值与测量值之间的误差达到最小，在多个学科领域具有广泛的应用［21-22］。假设Tc(z，tf)是由某一冰面温度变化μ(t)计算的冰温-深度剖面扰动，最小二乘法构造了如下式的误差函数J作为目标函数。当J取最小值时对应的边界条件μ(t)为最优解［7］。

为了表示和计算方便，将μ(t)写成如下傅里叶级数形式［15，18］：

任意给定估计系数a0、am和bm的初值，将μ(t)作为边界条件代入式（2）计算Tc(z，tf)；并由式（5）求出J；a0、am和bm的迭代公式通过梯度下降法给出，如下式：

式中：n为迭代次数，γn为步长。直至J的取值达到给定的误差水平停止上述迭代过程，找到最优解。

不难计算，式（7）中的各项偏导值如下式：

式中：Wa0(z)、Wam(z)和Wbm(z)分别是当μ(t)= 0.5，和时在tf时刻的冰温-深度剖面扰动。将迭代过程中的式（8）代入式（7）最终可求出μ(t)。

2.2 Tikhonov正则化方法

已知热传导-冰流模型中的冰温-深度剖面扰动θ(z)求解边界条件μ(t)是一个典型的不适定问题。Tikhonov正则化方法常用来解决不适定问题。该方法是对式（5）的误差函数J加一个约束条件；这样的约束条件可以解释为先验信息。在该问题中，通常长期的冰面温度变化μ(t)不会存在很大的波动，而由最小二乘法构造的目标函数对μ(t)没有任何限制。因此，需要对μ(t)施加一个独立的约束。Tikhonov正则化方法构造了如下式的目标函数Ψ；当Ψ达到最小时，对应的边界条件μ(t)为最优解［10，15］。

式中：α为正则化参数（α> 0）；Ω{μ(t)}称为稳定函数，其取值与μ(t)的光滑性有关。定义如下式：

求Ψ最小值的方法与最小二乘法相同，式（6）中估计系数a0、am和bm的迭代公式也通过梯度下降法给出，如下式：

式中各项偏导值如下式：

不难计算，稳定函数Ω{μ(t)}可近似表示如下：

式中：ξm=α0+α1(2πm/tf)2，αi=αtf/2，i= 0，1。

式（12）中各式最右边一项分别为：

将式（14）代入式（12）可得出式（11）中的各项a0、am和bm取值，进而求出μ(t)。

式（9）中的正则化参数α，相当于对μ(t)添加了约束：α越大，约束越严格，μ(t)波动越小；α越小，约束越宽松，μ(t)波动越大。α调节着误差函数J和稳定函数Ω{μ(t)}之间的平衡。Tikhonov正则化方法不仅考虑了计算值与测量值之间的误差水平；也考虑了μ(t)的波动大小，即μ(t)的稳定性。可根据L-曲线准则求出合适的α［23］。具体做法是在平面直角坐标系中以(xi，yi) =(J，Ω{μ(t)})(i= 1，2，…，p)为坐标点做出曲线（p为正则化参数的取值个数）。曲线通常呈L 形；该L 形曲线的拐角处对应最佳的正则化参数。另外，考虑到J和Ω{μ(t)}的数量级问题，可对xi和yi两者或其中之一取对数后再做出L-曲线；例如，在平面直角坐标系中，坐标点也可设置为(xi，yi) =(log(J)，log(Ω{μ(t)}))［23］。

2.3 蒙特卡罗算法

蒙特卡罗算法将冰面温度变化μ(t)在不同时刻的取值μ→=(μ0，μ1，…，μtf)作为参数组合向量，随机选择μ→作为模型的边界条件，采用随机步的方法选取满足给定误差水平的参数进行统计分析［8，18］。具体的步骤如下：

步骤1：任意给定参数组合向量μ→n=作为模型边界条件，计算误差函数

式中：Tμ→n(z，tf)为由μ→n计算的冰温-深度剖面扰动。

步骤2：在附近找到一个新的向量计算为了高效地寻找所有可能的参数组合，每一步都是随机选择中的某一个参数，对该参数添加微小扰动得到。

步骤4：以概率Pn接受新的向量：若接受，则；若不接受，则重复步骤1 到4。

通过以上步骤可以看出，蒙特卡罗算法是以产生一系列逐步改善的参数组合为基础的。最后，选取满足给定误差水平的参数组合，运用统计学方法获取各个参数取值μi(i= 0，1，…，tf)的频数分布。这种分布通常存在区域极大值，对应参数μi最有可能（最大概率）的取值。需要说明的是，通过统计方法得到的μi的频数分布极大值可能不止一个；这反映了在对应时刻存在多个温度取值能够与测量值吻合较好［4，18］。在这种情况下，由于实际中冰面温度随时间变化是平滑连续的，短期内的温度变化幅度不会太大。因此，可通过μi邻近点的取值，找到合理的分布区域来确定μi。同时，为避免最终选取的参数组合之间存在相关性，也可对满足误差水平的μi进行等间距取值来统计频数分布［4］。

3 反演算法的应用和比较

3.1 构造数据模拟实验

为比较不同反演算法的重建结果，我们构造了式（3）中的模型参数和冰面温度变化μ(t)，进行理想条件下的数值实验。假设式（3）中的各个参数取值为：冰川厚度H= 500 m，年均净积累量0.3 m冰当量厚度；重建时间tf= 400 a；冰的热扩散率κ=1.3×10-6m2·s-1。冰面速度v(0) = 0.3 m·a-1，冰川垂直于地表面方向的运动速度v(z)随深度的增加呈线性减小，至冰川底部减小为零。模型采用Crank-Nicolson 有限差分法求解［24］，它在时间方向上是隐式的二阶方法，具有数值稳定的特点。

首先，分别以冰面温度变化μ1(t) = 2sin(2πt/tf)和μ2(t) = 2sin(πt/tf)作为边界条件代入式（2），计算对应的冰温-深度剖面扰动θ1(z)和θ2(z)，如图2 所示。然后，以θ1(z)和θ2(z)作为测量值重建冰面温度变化。图3 给出了在测量值为θ1(z)时不同反演算法的部分相关结果。最后，为验证解的稳定性，分别对θ1(z)和θ2(z)添加±0.02 ℃的随机扰动并比较结果，对应的重建结果如图4所示。

图2 分别由μ1(t) = 2sin(2πt/tf)（a）计算的冰温-深度剖面扰动θ1(z)（b）和μ2(t) = 2sin(πt/tf)（c）计算的冰温-深度剖面扰动θ2(z)（d）［（b）和（d）中的红色曲线是对θ1(z)和θ2(z)添加了±0.02 ℃的随机扰动］Fig. 2 The glacier surface temperature functions μ1(t) = 2sin (2πt/tf) （a） used to generate the temperature-depth profile θ1(z) （b） and μ2(t) = 2sin (πt/tf) （c） used to generate θ2(z) （d） ［red lines in （b） and （d）：the calculated temperature-depth profiles with the random errors ±0.02 ℃］

图3 测量值为θ1(z)时不同反演算法的相关结果［a0、am和bm初值为0（a）和1（b）时最小二乘法重建的冰面温度变化；Tikhonov正则化方法对应的L-曲线（c）和重建的冰面温度变化（0 < α < 0.01）（d）；蒙特卡罗算法不同时刻冰面温度变化值的频数分布（e）、（f）］Fig. 3 The reconstructed solutions be the least square method （a） and （b）， Tikhonov regularization method （c） and （d），and Monte Carlo algorithm （e） and （f） ［（c） and （d）： the L-curve and the reconstructed glacier surface temperature change with 0 < α < 0.01； the probability distributions of the past glacier surface temperatures at selected times using Monte Carlo algorithm， （e） and （f）］

图4 边界条件分别为μ1(t) = 2sin (2πt/tf)和μ2(t) = 2sin (πt/tf)时三种反演算法重建的冰面温度变化［对θ1(z)添加±0.02 ℃随机扰动前后的重建结果（a）、（b）；对θ2(z)添加±0.02 ℃随机扰动前后的重建结果（c）、（d）］Fig. 4 Reconstructed changes in glacier surface temperature under the boundary conditions of μ1(t) = 2sin(2πt/tf) and μ2(t) = 2sin(πt/tf) （black line） using the least square method （blue line）， Tikhonov regularization method （red line）and Monte Carlo algorithm （green line） ［the input data θ1(z) and θ2(z) without errors， （a） and （c），the input data θ1(z) and θ2(z) with the random errors ±0.02 ℃， （b） and （d）］

图3（a）和图3（b）由最小二乘法求得，分别是当式（6）中的估计系数a0、am和bm初值为0［图3（a）］和初值为1 时［图3（b）］的重建结果。当a0、am和bm初值为0时，重建结果与真实温度变化基本一致；而当a0、am和bm初值为1时，重建结果与真实温度变化存在较大差异：近期冰面温度变化波动明显且振荡幅度非常大。由此可知，应用最小二乘法时a0、am和bm初值对重建结果有显著影响。除此之外，梯度步长γ的选择对重建结果影响也很大：步长太大会使估计系数在迭代过程中不断增大，导致结果不收敛；而步长太小会使收敛速度太慢，可能需要迭代很多次，导致时间成本过高。因此，实际中往往需要设置一个合适的步长γ。实践表明，本例中当γ取0.0001 时，重建结果与真实温度变化最接近［图3（a）中红色曲线］。图3（b）是当a0、am和bm初值为1，γ分别取0.001 和0.00001 时的重建结果。事实上，当γ取0.01 时，无论a0、am和bm初值如何选择，重建结果都趋于无穷大。

图3（c）和图3（d）是由Tikhonov正则化方法得到的L-曲线［其中坐标点为（xi，yi）=（J，log（Ω｛μ（t）｝））］和正则化参数取值0＜α＜0.01的重建结果。其中，估计系数a0、am和bm初值为1。由式（9）可知，当α取值为0 时，Tikhonov 正则化方法与最小二乘法的结果相同。与最小二乘法不同，应用Tikhonov 正则化方法只需要给定合适的梯度步长（0.0001），可设定任意值作为a0、am和bm的初值，均不影响最终结果。本例中L-曲线拐角处x的取值为0.0012，对应的α取值为0.00014。这里将α限定在［0，0.01］范围内，是因为在这个范围内的L-曲线上能够找到合适的拐角［21］。当α＞0.01 时，从结果上看，会使冰面温度的变化幅度非常小；从L-曲线形状上看，拐角右侧的曲线会向着几乎平行于x轴的方向继续向右延伸，不存在其他的拐角。由图3（d），冰面温度变化整体上呈正弦函数形式，只是由于正则化参数取值不同，温度变化振幅有所不同。随着α取值的增大，冰面温度变化幅度逐渐减小并趋于平滑和稳定，即α对边界条件起着施加约束的作用。可以看出，与最小二乘法相比，Tikhonov 正则化方法考虑了反问题中解的稳定性问题：由正则化参数α调节着误差函数和冰面温度波动大小之间的平衡。这也在一定程度上限制了模型的复杂度。

图3（e）和图3（f）是由蒙特卡罗算法这一统计学方法获取的不同时刻冰面温度变化值的频数分布，最终选取了满足给定误差水平的2 200 个参数组合进行统计分析。由于蒙特卡罗算法对随机选择的边界条件没有限制，即对模型的解没有任何约束；在满足误差函数很小的条件下，往往会出现冰面温度在短时间内波动非常大的结果，失去了实际意义。为避免这种情况，本例中将参数组合向量=(μ0，μ1，…，μtf)各个分量的取值控制在±2 ℃的范围内进行模拟。由图3（e），在t= 20 a和t= 140 a时刻，频数分布极大值对应的温度变化值为0.6 ℃和1.5 ℃，是参数μi在对应时刻最有可能的取值。事实上，这两个时刻真实的温度变化值为0.62 ℃和1.62 ℃，与计算值十分接近。结果表明，大多数时刻的温度变化值都有着类似图3（e）中的分布：能够较容易地确定极大值，即选取的μi与最大频数相对应。而图3（f）中，在t= 160 a时刻，频数分布极大值对应的温度变化值分别为-1.8 ℃和1.1 ℃；其中，最大频数对应的取值是-1.8 ℃。一般情况下，根据冰面温度变化的连续性，可通过邻近时刻的频数分布找到合理的分布区域，由该区域来确定极大值。在本例中，由于重建时间为百年尺度（400 a），长期平均冰面温度变化的振幅不会太大。事实上，在t=160 a 邻近时刻t= 120 a 至t= 180 a，频数分布均集中在0～2 ℃这一区域。因此，t= 160 a 时刻选取1.1 ℃较-1.8 ℃更为合理。该时刻的真实温度变化值为1.17 ℃，与1.1 ℃这一取值接近。需要指出的是，应用蒙特卡罗算法有时很难找到合理的分布区域，重建结果会在短期内存在明显振荡。在这种情况下，应该选用其他的反演算法进行求解。

由图4可以看出，最小二乘法、Tikhonov 正则化方法和蒙特卡罗算法的重建结果与真实温度变化趋势基本一致。其中，蒙特卡罗算法的结果不仅相对误差较大，而且也不平滑。对θ1(z)和θ2(z)添加±0.02 ℃的随机扰动后，Tikhonov 正则化方法的重建结果比最小二乘法和蒙特卡罗算法的结果更接近真实的温度变化过程，解的稳定性也最好。

3.2 真实钻孔数据模拟实验

真实钻孔资料受实地环境影响。应用热传导-冰流模型重建冰面温度变化需要考虑钻孔深度、年内冰面温度变化幅度、季节变化影响深度和式（4）中的稳态冰温垂直分布等实际问题。为与构造的数据模拟实验进行比较，我们以我国藏北高原腹地的马兰冰川钻孔（位置点35°48.4′ N，90°45.3′ E；海拔5 680 m；钻孔深度102 m）资料为例设定模型参数［25］。结合已有的、距离钻孔位置最近的五道梁气象站（35.13°N，93.05°E；海拔4 612 m；距离钻孔位置约220 km）的观测数据假定边界条件，位置如图5所示。最后，分别将假定的10 a 和40 a 冰面温度变化代入式（2），计算对应的冰温-深度剖面扰动作为测量值反演边界条件并比较结果。另外，为验证解的稳定性，对由40 a 冰面温度变化计算的冰温-深度剖面扰动添加±0.01 ℃的随机扰动并对比结果。

图5 藏北高原腹地马兰冰川钻孔位置图Fig. 5 Location map： the Malan Glacier on the northern Qinghai-Tibet Plateau showing the borehole site

马兰冰川位于我国藏北高原腹地的可可西里地区，最高峰海拔为6 056 m，是一个很大的典型山地冰帽冰川群；其南部平缓，冰面开阔平坦，属于极大陆型冰川类型［26］。马兰冰川发育于寒冷干燥的气候环境下，其严寒程度与东南极冰盖边缘地区相近。有关研究表明，位于藏北高原腹地的冰川变化幅度较边缘山区要小，而自小冰期以来马兰冰川的变化非常小，表明了它处于比较稳定的状态［27］。

由于研究区域是地形复杂的山区，而五道梁气象站距离马兰冰川钻孔位置约220 km 且海拔相差较大。因此，气温观测值与冰面温度值相差较大。尽管如此，式（3）中的μ(t)并非冰面温度的取值（绝对值），而是冰面温度相对于稳态冰面温度U(0)的变化值（相对值）。因此，我们利用五道梁气象站记录的10 a 和40 a 气温变化幅度（趋势）假定边界条件。为验证这一假定的准确性，我们同时获取了MODIS 地表温度产品（MOD11A1）中钻孔处像元在钻取后一年（2000 年）的月平均地表温度数据［数据来源：美国航空航天局（NASA）的地球科学数据系统Earthdata，https：//www.earthdata.nasa.gov/］与五道梁气象站记录的月平均气温变化幅度进行对比。

3.2.1 季节变化影响深度

马兰冰川钻孔于1999 年5 月钻取并测温［25］。钻孔底部未到冰床；20 m 深度以下冰温变化幅度很小，仅为1.1 ℃。经计算，冰川厚度H为130 m，年均净积累量0.23 m 冰当量厚度；冰床温度-1.7 ℃。冰面竖直速度v(0) = 0.23 m·a-1。假设冰川垂直于地表面方向的运动速度v(z)随深度的增加呈线性减小，至冰川底部减小为0。根据式（1）中各个物理参数取值的经验公式［13，26］：导热率λ=9.828e-0.0057(273.15+T)W·m-1·℃-1，比热容c=(2098 +7.122T) J·kg-1·℃-1。代入本例计算，得λ为2～2.1W·m-1·℃-1。冰川冰的密度ρ介于830～917 kg·m-3之间，而山地冰川冰ρ为（850 ± 60） kg·m-3具有广泛的适用性［28-29］。代入热扩散率公式，得κ为1.2 × 10-6～1.4 × 10-6m2·s-1。经验证，本例中的导热率和热扩散率取值在上述范围内对实验结果影响较小。因此，为便于计算，模型中取λ= 2 W·m-1·℃-1；κ= 1.4 × 10-6m2·s-1。由于各地的地热流密度q差异很大，我们由马兰冰川钻孔近底部的温度梯度计算出地热流密度约为0.04 W·m-2。另外，根据《中国大陆地区大地热流数据汇编（第三版）》［30］，现有距马兰冰川钻孔最近、具有地热流数据实测值的位置点是柴达木（38°14′ N，90°47′ E），如图5所示。该位置点的校正热流为0.043 W·m-2，与计算值接近。由上述冰床温度和地热流密度推算式（4）中的长期平均冰面温度Us为-4.2 ℃，由此计算稳态的冰温-深度剖面U(z)。

通常冰面以下十几米（冰温年变化深度）的钻孔温度反映的是年内的季节变化。因此，重建长时间尺度的冰面温度变化过程需要选取季节影响深度以下的冰温-深度剖面。计算该深度需要估算一年内的冰面温度变化幅度。由于马兰冰川钻孔的钻取时间是1999年，我们参考五道梁气象站前一年（1998 年）月平均气温的变化幅度［数据来源：国家气象科学数据中心（CMDSC），http：//data. cma. cn］来假定一年内冰面温度变化这一边界条件μ3(t)，由此计算季节影响深度，如图6（a）所示。气象资料表明，一年的气温变化幅度约为20 ℃（-15～5 ℃）。因为冰面温度不会高于0 ℃，即夏季冰面温度最高升至0 ℃。因此，我们假定一年的冰面温度变化是取值范围在-20～0 ℃的正弦函数，如图6（b）所示。为验证该冰面温度变化幅度的准确性，我们将钻孔处像元2000 年的月平均地表温度数据与其对比。结果表明，2000 年钻孔处的最低月平均地表温度约为-23 ℃，与假定的年内冰面温度变化（-20～0 ℃）基本一致，说明了该假设可靠。将其作为边界条件代入式（2），得到冰温-深度剖面扰动θ3(z)和扰动后的冰温-深度剖面如图6（c）和图6（d）所示。

图6 五道梁气象站1998年的月平均气温（a），假定的一年内冰面温度变化（b），由（b）计算的冰温-深度剖面扰动θ3(z)（c）以及稳态与季节影响剖面（d）Fig. 6 Monthly mean air temperature at the Wudaoliang meteorological station in 1998 （a）， the yearly glacier surface temperature oscillations （b）， the temperature-depth profile θ3(z) generated from （b） （c）， and the temperature log constructed by adding θ3(z) to the steady-state temperature-depth profile U(z) （d）

由图6，由于冰面温度的平均值（-10 ℃）比长期平均冰面温度Us（-4.2 ℃）低很多，所以θ3(z)呈明显的负向扰动。当温度剖面扰动小于0.1 ℃时，可认为在该深度上，年内温度变化的振幅消失［15］。经计算，θ3(16)=-0.14 ℃，θ3(17)=-0.08 ℃，即季节影响深度为16 m（冰温年变化深度）。

3.2.2 10 a尺度重建结果比较

在冰面消融微弱的条件下，可通过重建季节影响深度处的冰温变化代替年均冰面温度变化。由3.2.1 节可知，季节影响深度为16 m。我们以16 m的冰温变化为边界条件，用该深度以下的冰温-深度剖面重建1988—1998 年的冰面温度变化过程。图7（a）是五道梁气象站1988—1998年的年均气温。用该时期的气温变化趋势作为冰面温度变化这一边界条件μ4(t)，计算的冰温-深度剖面扰动θ4(z)如图7（b）所示。

图7 五道梁气象站1988—1998年的年均气温及变化趋势（a），由（a）计算的冰温-深度剖面扰动θ4(z)（b）以及三种反演算法重建的冰面温度变化（c）Fig. 7 Annual mean air temperature from 1988 to 1998 collected by the Wudaoliang meteorological station （a），the temperature-depth profile θ4(z) generated from （a） （b）， and the inversion results with the least square method（blue line）， Tikhonov regularization method （red line） and Monte Carlo algorithm （green line） （c）

由图7（b），θ4(z)在深度约50 m 处小于0.1 ℃，60 m深度以下小于0.05 ℃，即10 a冰面温度变化信号向下传播至50 m 左右。由于冰面温度变化信号缓慢向下传播，钻孔越深处的温度扰动对应着更早期的冰面温度变化。从深度100 m 向上至30 m，负向扰动逐渐增大；而30 m 以上的负向扰动逐渐减小；这说明了冰面温度变化大致经历了先降温，再升温的过程。事实上，将θ4(z)由下而上观察，可以看到与图7（a）的冰面温度变化相似的趋势，即由θ4(z)可大致了解冰面温度变化趋势。但是，冰面温度变化的具体时间和振幅需要通过相关的反演算法求解。

我们以θ4(z)为已知条件，分别应用最小二乘法、Tikhonov 正则化方法和蒙特卡罗算法反演边界条件，分析各个算法的重建结果并讨论这些算法的优劣性。图7（c）给出了上述三种反演算法10 a 尺度的重建结果。整体结果与图4 类似：最小二乘法和Tikhonov 正则化方法的重建结果较好；而蒙特卡罗算法的结果不连续，振荡幅度大，准确性也较低。

3.2.3 40 a尺度重建结果比较与稳定性检验

与3.2.2 节的研究方法类似，这里我们结合五道梁气象站1959—1998 年这40 a 的年均气温资料，如图8（a）所示。用该时期的气温变化趋势作为边界条件μ5(t)，计算的冰温-深度剖面扰动θ5(z)如图8（b）所示。

图8 五道梁气象站1959—1998年的年均气温及变化趋势（a），由（a）计算的冰温-深度剖面扰动θ5(z)（b）以及分别对θ5(z)添加±0.01 ℃随机扰动前后［（c）和（d）］三种反演算法重建的冰面温度变化Fig. 8 Annual mean air temperature from 1959 to 1998 collected by the Wudaoliang meteorological station （a）， the temperaturedepth profile θ5(z) generated from （a） （b）， the inversion results by the input data θ5(z) with the least square method（blue line）， Tikhonov regularization method （red line） and Monte Carlo algorithm （green line） （c），and the inversion results by the input data θ5(z) with the random errors ±0.01 ℃ （d）

由图8（b）可知，40 a 冰面温度变化信号向下传播至78 m 左右［θ5(78) < 0.1 ℃］。θ5(z)在任意深度均为正向扰动，且扰动自下而上逐渐增大，说明了冰面温度变化的升温趋势。图8（c）给出了三种反演算法40 a 尺度的重建结果。对θ5(z)添加±0.01 ℃的随机扰动，重建结果如图8（d）所示。由图8 可知，最小二乘法和Tikhonov 正则化方法的重建结果与真实的温度变化接近；Tikhonov 正则化方法的结果最平滑和稳定；而蒙特卡罗算法结果的准确性和稳定性都较差。

将3.1 节构造数据模拟实验结果与3.2 节利用马兰冰川真实钻孔数据模拟实验结果进行对比，可得到以下一致的结论：最小二乘法、Tikhonov 正则化方法和蒙特卡罗算法这三种反演算法的重建结果与真实的冰面温度变化基本一致；蒙特卡罗算法的误差较大，且结果不平滑；Tikhonov 正则化方法是目前求解该反问题的最优方法，能够有效降低噪声干扰，使得到的解稳定。

需要说明的是，实际中马兰冰川的表面积雪消融和融水的再冻结等依然存在，对冰体和冰面水热平衡有一定影响。而应用热传导-冰流模型无法考虑这些因素的影响，这可能导致重建结果中近期的升温幅度偏大［10］。尽管如此，为了更好地比较反演算法，在利用马兰冰川钻孔数据模拟实验中，冰面温度变化是已知的；用于重建冰面温度变化的冰温-深度剖面扰动也是由给定的边界条件代入模型计算得到的。因此，重建结果与真实温度变化趋势基本一致。此外，由于实际中无论测量精度多高，总会存在一定的测量误差，由测量得到的冰温-深度剖面相当于添加了随机扰动。为了能够得到准确的结果，可对测量的冰温-深度剖面平滑后再进行反演。

4 结论

本文基于冰面温度变化信号在冰层中的传播特性和热传导-冰流模型，研究并比较了利用冰川钻孔温度重建冰面温度变化这一反问题的反演算法，包括了最小二乘法、Tikhonov 正则化方法和蒙特卡罗算法。由于该反问题存在解的唯一性和稳定性问题，不论是哪种反演算法，目标都是找到相对稳定的边界条件使得冰温-深度剖面的计算值最大程度地接近测量值，实质是最优化问题的求解。

由于蒙特卡罗算法随机选择边界条件，再输入模型进行计算和统计，因此，它可看作是正问题的求解方法。主要存在两个问题：一是冰面温度变化值的频数分布会出现多个极值；二是求解的冰面温度变化不平滑。尽管如此，因为蒙特卡罗算法是将模型的未知参数组合作为多维向量寻找最优解，所以，在实际求解过程中，在部分参数（如初始条件等）未知的情况下，可选用蒙特卡罗算法这一求解方法。即考虑到能够获取的地学要素的限制，蒙特卡罗算法是最优的。最小二乘法中估计系数a0、am和bm初值对重建结果影响很大，且仅考虑了误差函数这一目标函数，因此，它的稳定性差。Tikhonov正则化方法将误差函数和稳定函数相加作为目标函数，综合考虑了解的准确性和稳定性问题，是目前求解该反问题的最优方法。与其他反演算法相比，Tikhonov 正则化方法能够有效降低噪声干扰，使得到的解稳定。