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

图1 冰面温度变化对稳态冰温垂直分布的影响

Fig. 1 Influence of ice surface temperature change on the vertical distribution of steady-state ice temperature

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

\{\begin{array}{l} v (z) \frac{d U}{d z} = κ \frac{d}{d z} (\frac{d U}{d z}) \\ U (0) = U_{s} \\ \frac{d U}{d z} (H) = - \frac{q}{λ} \end{array}

（4）

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

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

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

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

2 反演算法

2.1　最小二乘法

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

J = \frac{1}{2} \int_{0}^{H} {[T_{c} (z, t_{f}) - θ (z)]}^{2} d z

（5）

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

μ (t) = \frac{a_{0}}{2} + \sum_{m = 1}^{M} a_{m} c o s (\frac{2 π m t}{t_{f}}) + b_{m} s i n (\frac{2 π m t}{t_{f}})

（6）

任意给定估计系数 $a_{0}$ 、 $a_{m}$ 和 $b_{m}$ 的初值，将 $μ (t)$ 作为边界条件代入式（2）计算 $T_{c} (z, t_{f})$ ；并由式（5）求出 $J$ ； $a_{0}$ 、 $a_{m}$ 和 $b_{m}$ 的迭代公式通过梯度下降法给出，如下式：

a_{0}^{n + 1} = a_{0}^{n} - γ^{n} \frac{\partial J^{n}}{\partial a_{0}^{n}}

，

a_{m}^{n + 1} = a_{m}^{n} - γ^{n} \frac{\partial J^{n}}{\partial a_{m}^{n}}

，

b_{m}^{n + 1} = b_{m}^{n} - γ^{n} \frac{\partial J^{n}}{\partial b_{m}^{n}}

，

m = 1,2, \dots, M,

（7）

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

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

\{\begin{array}{l} \frac{\partial J^{n}}{\partial a_{0}^{n}} = \int_{0}^{H} W_{a_{0}} (z) [T_{c}^{n} (z, t_{f}) - θ (z)] d z \\ \frac{\partial J^{n}}{\partial a_{m}^{n}} = \int_{0}^{H} W_{a_{m}} (z) [T_{c}^{n} (z, t_{f}) - θ (z)] d z \\ \frac{\partial J^{n}}{\partial b_{m}^{n}} = \int_{0}^{H} W_{b_{m}} (z) [T_{c}^{n} (z, t_{f}) - θ (z)] d z \end{array}

（8）

式中： $W_{a_{0}} (z)$ 、 $W_{a_{m}} (z)$ 和 $W_{b_{m}} (z)$ 分别是当 $μ (t) = 0.5$ ， $μ (t) = c o s (\frac{2 π m t}{t_{f}})$ 和 $μ (t) = s i n (\frac{2 π m t}{t_{f}})$ 时在 $t_{f}$ 时刻的冰温-深度剖面扰动。将迭代过程中的式（8）代入式（7）最终可求出 $μ (t)$ 。

2.2　Tikhonov正则化方法

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

Ψ = J + α Ω \{μ (t)\}

（9）

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

Ω \{μ (t)\} = \int_{0}^{t_{f}} [μ^{2} (t) + {(\frac{d μ (t)}{d t})}^{2}] d t

（10）

求 $Ψ$ 最小值的方法与最小二乘法相同，式（6）中估计系数 $a_{0}$ 、 $a_{m}$ 和 $b_{m}$ 的迭代公式也通过梯度下降法给出，如下式：

a_{0}^{n + 1} = a_{0}^{n} - γ^{n} \frac{\partial Ψ^{n}}{\partial a_{0}^{n}}, a_{m}^{n + 1} = a_{m}^{n} - γ^{n} \frac{\partial Ψ^{n}}{\partial a_{m}^{n}}, b_{m}^{n + 1} = b_{m}^{n} - γ^{n} \frac{\partial Ψ^{n}}{\partial b_{m}^{n}}

（11）

式中各项偏导值如下式：

\{\begin{array}{l} \frac{\partial Ψ^{n}}{\partial a_{0}^{n}} = \frac{\partial J^{n}}{\partial a_{0}^{n}} + α \frac{\partial Ω^{n}}{\partial a_{0}^{n}} \\ \frac{\partial Ψ^{n}}{\partial a_{m}^{n}} = \frac{\partial J^{n}}{\partial a_{m}^{n}} + α \frac{\partial Ω^{n}}{\partial a_{m}^{n}} \\ \frac{\partial Ψ^{n}}{\partial b_{m}^{n}} = \frac{\partial J^{n}}{\partial b_{m}^{n}} + α \frac{\partial Ω^{n}}{\partial b_{m}^{n}} \end{array}

（12）

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

α Ω \{μ (t)\} = α_{0} \frac{a_{0}^{2}}{2} + \sum_{m = 1}^{M} (a_{m}^{2} + b_{m}^{2}) ξ_{m}

（13）

式中： $ξ_{m} = α_{0} + α_{1} {(2 π m / t_{f})}^{2}$ ， $α_{i} = α t_{f} / 2, i = 0,1$ 。

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

α \frac{\partial Ω^{n}}{\partial a_{0}^{n}} = α_{0} a_{0}^{n}, α \frac{\partial Ω^{n}}{\partial a_{m}^{n}} = 2 a_{m}^{n} ξ_{m}, α \frac{\partial Ω^{n}}{\partial b_{m}^{n}} = 2 b_{m}^{n} ξ_{m}

（14）

将式（14）代入式（12）可得出式（11）中的各项 $a_{0}$ 、 $a_{m}$ 和 $b_{m}$ 取值，进而求出 $μ (t)$ 。

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

2.3　蒙特卡罗算法

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

步骤1：任意给定参数组合向量 $\vec{μ}$ $^{n} = (μ_{0}^{n}, μ_{1}^{n}, \dots, μ_{t_{f}}^{n}$ ）作为模型边界条件，计算误差函数

S ({\vec{μ}}^{n}) = \int_{0}^{H} {[T_{{\vec{μ}}^{n}} (z, t_{f}) - θ (z)]}^{2} d z

，（15）

式中： $T_{{\vec{μ}}^{n}} (z, t_{f})$ 为由 ${\vec{μ}}^{n}$ 计算的冰温-深度剖面扰动。

步骤2：在 ${\vec{μ}}^{n}$ 附近找到一个新的向量 ${\vec{μ}}^{n + 1}$ $= (μ_{0}^{n + 1}, μ_{1}^{n + 1}, \dots, μ_{t_{f}}^{n + 1}$ ），计算 $S ({\vec{μ}}^{n + 1})$ 。为了高效地寻找所有可能的参数组合，每一步都是随机选择 ${\vec{μ}}^{n}$ 中的某一个参数，对该参数添加微小扰动得到 ${\vec{μ}}^{n + 1}$ 。

步骤3：定义转移概率 $P_{n} = m i n (1, e^{- (S ({\vec{μ}}^{n + 1}) - S ({\vec{μ}}^{n}))})$ 。

步骤4：以概率 $P_{n}$ 接受新的向量 ${\vec{μ}}^{n + 1}$ ：若接受 ${\vec{μ}}^{n + 1}$ ，则 ${\vec{μ}}^{n}$ $= {\vec{μ}}^{n + 1}$ ；若不接受 ${\vec{μ}}^{n + 1}$ ，则 ${\vec{μ}}^{n}$ $= {\vec{μ}}^{n}$ 。重复步骤1到4。

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

3 反演算法的应用和比较

3.1　构造数据模拟实验

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

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

图2

图2 分别由 $μ_{1} (t) = 2 s i n (2 π t / t_{f})$ （a）计算的冰温-深度剖面扰动 $θ_{1} (z)$ （b）和 $μ_{2} (t) = 2 s i n (π t / t_{f})$ （c）计算的冰温-深度剖面扰动 $θ_{2} (z)$ （d）［（b）和（d）中的红色曲线是对 $θ_{1} (z)$ 和 $θ_{2} (z)$ 添加了 $\pm 0.02$ ℃的随机扰动］

Fig. 2 The glacier surface temperature functions $μ_{1} (t) = 2 s i n (2 π t / t_{f})$ （a） used to generate the temperature-depth profile $θ_{1} (z)$ （b） and $μ_{2} (t) = 2 s i n (π t / t_{f})$ （c） used to generate $θ_{2} (z)$ （d）［red lines in （b） and （d）： the calculated temperature-depth profiles with the random errors $\pm 0.02$ ℃］

图3

图3 测量值为 $θ_{1} (z)$ 时不同反演算法的相关结果［ $a_{0}$ 、 $a_{m}$ 和 $b_{m}$ 初值为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

图4 边界条件分别为 $μ_{1} (t) = 2 s i n (2 π t / t_{f})$ 和 $μ_{2} (t) = 2 s i n (π t / t_{f})$ 时三种反演算法重建的冰面温度变化［对 $θ_{1} (z)$ 添加 $\pm 0.02$ ℃随机扰动前后的重建结果（a）、（b）；对 $θ_{2} (z)$ 添加 $\pm 0.02$ ℃随机扰动前后的重建结果（c）、（d）］

Fig. 4 Reconstructed changes in glacier surface temperature under the boundary conditions of $μ_{1} (t) = 2 s i n (2 π t / t_{f})$ and $μ_{2} (t) = 2 s i n (π t / t_{f})$ （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 $\pm 0.02$ ℃，（b） and （d）］

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

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

图3（e）和图3（f）是由蒙特卡罗算法这一统计学方法获取的不同时刻冰面温度变化值的频数分布，最终选取了满足给定误差水平的2 200个参数组合进行统计分析。由于蒙特卡罗算法对随机选择的边界条件没有限制，即对模型的解没有任何约束；在满足误差函数很小的条件下，往往会出现冰面温度在短时间内波动非常大的结果，失去了实际意义。为避免这种情况，本例中将参数组合向量 $\vec{μ} = (μ_{0}, μ_{1}, \dots, μ_{t_{f}})$ 各个分量的取值控制在 $\pm 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)$ 添加 $\pm 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冰面温度变化计算的冰温-深度剖面扰动添加 $\pm 0.01$ ℃的随机扰动并对比结果。

图5

图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.828 e^{- 0.0057 (273.15 + T)}$ W·m^-1·℃^-1，比热容 $c = (2098 + 7.122 T)$ J·kg^-1·℃^-1。代入本例计算，得 $λ$ 为2～2.1 W·m^-1·℃^-1。冰川冰的密度 $ρ$ 介于830～917 kg·m^-3之间，而山地冰川冰 $ρ$ 为（ $850 \pm 60$ ） kg·m^-3具有广泛的适用性^［28-29］。代入热扩散率公式 $κ = \frac{λ}{c ρ}$ ，得 $κ$ 为 $1.2 \times 10^{- 6}$ ～ $1.4 \times 10^{- 6}$ m²·s^-1。经验证，本例中的导热率和热扩散率取值在上述范围内对实验结果影响较小。因此，为便于计算，模型中取 $λ = 2$ W·m^-1·℃^-1； $κ = 1.4 \times 10^{- 6}$ m²·s^-1。由于各地的地热流密度 $q$ 差异很大，我们由马兰冰川钻孔近底部的温度梯度计算出地热流密度约为0.04 W·m^-2。另外，根据《中国大陆地区大地热流数据汇编（第三版）》^［30］，现有距马兰冰川钻孔最近、具有地热流数据实测值的位置点是柴达木（38°14′ N，90°47′ E），如图5所示。该位置点的校正热流为0.043 W·m^-2，与计算值接近。由上述冰床温度和地热流密度推算式（4）中的长期平均冰面温度 $U_{s}$ 为-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

图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 ℃）比长期平均冰面温度 $U_{s}$ （-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

图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

图8 五道梁气象站1959—1998年的年均气温及变化趋势（a），由（a）计算的冰温-深度剖面扰动 $θ_{5} (z)$ （b）以及分别对 $θ_{5} (z)$ 添加 $\pm 0.01$ ℃随机扰动前后［（c）和（d）］三种反演算法重建的冰面温度变化

Fig. 8 Annual mean air temperature from 1959 to 1998 collected by the Wudaoliang meteorological station （a）， the temperature-depth 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 $\pm 0.01$ ℃ （d）

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

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

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

4 结论

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

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

参考文献

原文顺序

文献年度倒序

文中引用次数倒序

被引期刊影响因子

[1]

Van de Wal

R S W

， Mulvaney

， Isaksson

， et al.

Reconstruction of the historical temperature trend from measurements in a medium-length borehole on the Lomonosovfonna plateau， Svalbard

［J］. Annals of Glaciology， 2002， 35： 371-378.

[2]

Bodri

， Cermak

. Borehole Climatology： a new method on how to reconstruct climate［M］. Amsterdam： Elsevier， 2007.

[3]

Zagorodnov

， Nagornov

， Thompson

L G

Influence of air temperature on a glacier’s active-layer temperature

［J］. Annals of Glaciology， 2006， 43： 285-291.

[4]

Dahl-Jensen

， Mosegaard

， Gundestrup

， et al.

Past temperatures directly from the Greenland ice sheet

［J］. Science， 1998， 282 （5387）： 268-271.

[本文引用: 5]

[5]

Liu

Jia

Reconstruction of the past ground surface temperature changes using borehole paleothermometry

［D］. Lanzhou： Lanzhou University， 2015.

刘佳

利用钻孔温度梯度重建过去地表温度变化研究

［D］. 兰州：兰州大学， 2015.

[6]

Dahl-Jensen

， Johnsen

S J

Palaeotemperatures still exist in the Greenland ice sheet

［J］. Nature， 1986， 320（6059）： 250-252.

[7]

Macayeal

D R

Paleothermometry by control methods

［J］. Journal of Glaciology， 1991， 37（127）： 326-338.

[8]

Dahl-Jensen

， Morgan

V I

， Elcheikh

Monte Carlo inverse modelling of the Law Dome （Antarctica） temperature profile

［J］. Annals of Glaciology， 1999， 29（1）： 145-150.

[9]

Kotlyakov

V M

， Arkhipov

S M

， Henderson

K A

， et al.

Deep drilling of glaciers in Eurasian Arctic as a source of paleoclimatic records

［J］. Quaternary Science Reviews， 2004， 23： 1371-1390.

[10]

Nagornov

O V

， Konovalov

Y V

， Tchijov

Temperature reconstruction for Arctic glaciers

［J］. Palaeogeography， Palaeoclimatology， Palaeoecology， 2006， 236（1/2）： 125-134.

[11]

Barrett

B E

， Nicholls

K W

， Murray

， et al.

Rapid recent warming on Rutford Ice Stream， West Antarctica， from borehole thermometry

［J］. Geophysical Research Letters， 2009， 36 （2）： 349-363.

[12]

Muto

Multi-decadal surface temperature trends in East Antarctica inferred from borehole firn temperature measurements and geophysical inverse methods

［D］. Boulder， America： University of Colorado， 2010.

[13]

Zagorodnov

， Nagornov

O V

， Thompson

L G

， et al.

Borehole temperatures reveal details of 20th century warming at Bruce Plateau， Antarctic Peninsula

［J］. The Cryosphere， 2012， 6（3）： 675-686.

[14]

Yang

J W

， Han

， Orsi

A J

， et al.

Surface temperature in twentieth century at the Styx Glacier， northern Victoria Land， Antarctica， from borehole thermometry

［J］. Geophysical Research Letters， 2018， 45： 9834-9842.

[15]

Sun

， Wang

， Hou

Twentieth century warming reflected by the Malan Glacier borehole temperatures， northern Tibetan Plateau

［J］. Arctic， Antarctic， and Alpine Research， 2021， 53（1）： 227-236.

[本文引用: 4]

[16]

Robin G de

Ice movement and temperature distribution in glaciers and ice sheets

［J］. Journal of Glaciology， 1955， 2（18）： 523-532.

[17]

Dansgaard

， Johnsen

S J

A flow model and a time scale for the ice core from Camp Century， Greenland

［J］. Journal of Glaciology， 1969， 8（53）： 215-223.

[18]

Nagornov

O V

， Tyuflin

S A

， Nikitaev

V G

， et al.

Comparative analysis of the past glacier surface temperatures based on various probe-functions

［J］. Journal of Physics Conference Series， 2017， 788（1）： 012055.

[本文引用: 4]

[19]

Zhen

， Yao

Tandong

， Tian

Lide

， et al.

Borehole temperature at the ice-core drilling site in the Muztag Ata Glacier， east Pamirs

［J］. Journal of Glaciology and Geocryology， 2004， 26（3）： 284-288.

李真，姚檀栋，田立德，等.

慕士塔格冰川海拔7 000 m处冰芯钻孔温度

［J］. 冰川冻土， 2004， 26（3）： 284-288.

[20]

Guangjian

， Yao

Tandong

， Xu

Baiqing

， et al.

Ice-core borehole temperature in the Muztag Ata， east Pamirs

［J］. Journal of Glaciology and Geocryology， 2003， 25（6）： 676-679.

邬光剑，姚檀栋，徐柏青，等.

慕士塔格冰芯钻孔温度测量结果

［J］. 冰川冻土， 2003， 25（6）： 676-679.

[21]

Jia

Xiaofeng

， Hu

Tianyue

Solving seismic wave equation by moving least squares （MLS） method

［J］. Progress in Geophysics， 2005， 20（4）： 920-924.

贾晓峰，胡天跃.

滑动最小二乘法求解地震波波动方程

［J］. 地球物理学进展， 2005， 20（4）： 920-924.

[22]

Shen

P Y

， Beck

A E

Least squares inversion of borehole temperature measurements in functional space

［J］. Journal of Geophysical Research， 1991， 96（B12）： 19965-19979.

[23]

Calvetti

， Morigi

， Reichel

， et al.

Tikhonov regularization and the L-curve for large discrete ill-posed problems

［J］. Journal of Computational and Applied Mathematics， 2000， 123（1/2）： 423-446.

[24]

Sweilam

N H

， Khader

M M

， Mahdy

A M S

Crank-Nicolson finite difference method for solving time-fractional diffusion equation

［J］. Journal of Fractional Calculus and Applications， 2012， 2（2）： 1-9.

[25]

Wang

Ninglian

， Yao

Tandong

， Pu

Jianchen

， et al.

Climatic and environmental changes over the last millennium recorded in the Malan ice core from the northern Tibetan Plateau

［J］. Science in China Series D： Earth Sciences， 2006， 36（8）： 723-732.

王宁练，姚檀栋，蒲建辰，等.

青藏高原北部马兰冰芯记录的近千年来气候环境变化

［J］. 中国科学（D辑）， 2006， 36（8）： 723-732.

[26]

Jianchen

， Yao

Tandong

， Wang

Ninglian

， et al.

Recently variation of Malan Glacier in Hoh Xil region， center of Tibetan Plateau

［J］. Journal of Glaciology and Geocryology， 2001， 23（2）： 189-192.

蒲健辰，姚檀栋，王宁练，等.

可可西里马兰山冰川的近期变化

［J］. 冰川冻土， 2001， 23（2）： 189-192.

[27]

Xie

Zichu

， Han

Jiankang

， Feng

Qinghua

， et al.

Primary study on the glaciers of Mountain Malan， Hoh Xil region， Qinghai-Xizang Plateau

［J］. Journal of Natural Science of Hunan Normal University， 2000， 23（1）： 83-88.

谢自楚，韩健康，冯清华，等.

青藏高原可可西里地区马兰山冰川的初步研究

［J］. 湖南师范大学学报（自然科学学报）， 2000， 23（1）： 83-88.

[28]

Paterson

W S B

. The physics of glaciers［M］. Oxford： Pergamon， 1994.

[29]

Huss

Density assumptions for converting geodetic glacier volume change to mass change

［J］. The Cryosphere， 2013， 7（3）： 877-887.

[30]

Shengbiao

， He

Lijuan

， Wang

Jiyang

Compilation of heat flow data in the China continental area （3rd edition）

［J］. Chinese Journal of Geophysics， 2001， 44（5）： 611-626.

胡圣标，何丽娟，汪集旸.

中国大陆地区大地热流数据汇编（第三版）

［J］. 地球物理学报， 2001， 44（5）： 611-626.