基于耦合熱傳導-冰流模型的冰川鉆孔溫度重建的氣候反演算法比較

孫 歡， 王寧練， 張 泉

（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)的穩定性?？筛鶕﨤-曲線準則求出合適的α［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 正則化方法能夠有效降低噪聲干擾，使得到的解穩定。