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

孫 歡， 王寧練， 張 泉

（1. 陜西省地表系統(tǒng)與環(huán)境承載力重點實驗室，陜西 西安 710127； 2. 西北大學(xué) 城市與環(huán)境學(xué)院地表系統(tǒng)與災(zāi)害研究院，陜西 西安 710127）

0 引言

對過去冰面溫度變化過程的研究是氣候重建和預(yù)測冰川對未來氣候變化響應(yīng)的重要問題之一［1］。由于熱傳導(dǎo)的作用，冰面溫度變化信號向冰面以下傳播，對穩(wěn)態(tài)冰溫垂直分布造成了相應(yīng)的擾動。通過測量冰川鉆孔溫度可知擾動后的冰溫-深度剖面，該剖面能夠保存過去幾個世紀(jì)甚至更長時間的冰面溫度變化信息［2-3］。通過冰川鉆孔溫度方法重建氣候變化歷史就是建立在冰面溫度變化和冰溫垂直分布的直接聯(lián)系之上；考慮到冰川運動的影響，可建立耦合的熱傳導(dǎo)-冰流物理模型來描述冰面溫度變化信號向下傳播的過程［4］。模型以穩(wěn)態(tài)冰溫垂直分布和冰面溫度變化為初始條件和邊界條件，其中穩(wěn)態(tài)冰溫垂直分布由冰川底部地?zé)崃骱烷L期平均冰面溫度共同確定。與求解模型中任意時刻的冰溫-深度剖面這一正問題不同，通過測量冰川鉆孔溫度重建冰面溫度變化是求解模型的邊界條件，是一個典型的反問題。由于熱擴散的影響，利用冰川鉆孔溫度重建氣候變化歷史具有較低的分辨率；雖然無法體現(xiàn)短期氣候事件，卻可以在一定時間尺度得到完整的冰面溫度變化過程［5］。

近三十年來，冰川鉆孔溫度方法在兩極和高緯度冷冰川地區(qū)的氣候重建研究方面得到了廣泛的應(yīng)用［6-15］。這部分地區(qū)缺少早期的氣象觀測數(shù)據(jù)，冰川鉆孔溫度方法無疑提供了一個獨立完善的研究工具，具有十分重要的意義。在實際應(yīng)用中，受冰川鉆孔深度、實地環(huán)境等方面的影響，應(yīng)用冰川鉆孔溫度方法重建氣候有著十分嚴格的適用條件。要求冰川處于相對穩(wěn)定的狀態(tài)，底部凍結(jié)在基巖上無滑動過程；冰面消融較少，不涉及徑流、滲浸作用等熱交換過程，即冰川內(nèi)部相變對鉆孔溫度的影響較小；冰川鉆孔位于分冰嶺附近，水平方向的運動速度可忽略等。熱傳導(dǎo)-冰流模型的反演算法是該研究的基礎(chǔ)和關(guān)鍵。以往對于熱傳導(dǎo)-冰流模型的反演算法研究存在不足，如：沒有詳細討論該模型解的穩(wěn)定性問題；對各類反演算法的優(yōu)劣性和適用性也沒有具體的比較和說明。本文首先通過構(gòu)造冰面溫度變化和模型參數(shù)，進行理想條件下的數(shù)值實驗，分析了最小二乘法、Tikhonov 正則化方法和蒙特卡羅算法這三種反演算法重建百年尺度冰面溫度變化的結(jié)果。然后，以我國藏北高原腹地的馬蘭冰川鉆孔資料為例設(shè)定模型參數(shù)，結(jié)合氣象觀測數(shù)據(jù)，詳細比較了上述反演算法在10 a 和40 a 尺度的重建結(jié)果。最后，對構(gòu)造理想條件下和應(yīng)用真實鉆孔數(shù)據(jù)的兩組數(shù)值實驗結(jié)果進行了比較，并進一步討論該反問題中解的穩(wěn)定性問題和不同算法的優(yōu)劣性，指出算法的適用條件。本項研究可為今后利用冰溫-深度剖面進行氣候重建研究提供借鑒。

1 模型描述

早在1955 年，Robin［16］研究了冰流運動對冰川鉆孔溫度的影響。Dansgaard 等［17］1969 年建立了格陵蘭冰蓋的運動模式。這些研究為利用冰川鉆孔溫度重建冰面溫度變化奠定了基礎(chǔ)。耦合的熱傳導(dǎo)-冰流模型如下［4］：

式中：T為冰川冰在某深度任意時刻的溫度（℃）；t為時間（s）；v→為冰流運動速度（m·s-1）。ρ、c和λ分別為冰的密度（kg·m-3）、比熱容（J·kg-1·℃-1）和導(dǎo)熱率（W·m-1·℃-1）。為溫度隨時間的變化率；?T為溫度在空間各方向上的全微分，即溫度的空間變化率。f為冰川內(nèi)部與熱傳導(dǎo)過程無關(guān)的熱源項（J·s-1·m-3），如：冰川內(nèi)部由于復(fù)雜相變產(chǎn)生的熱量等。

在冰川內(nèi)部與熱傳導(dǎo)過程無關(guān)的熱源項、冰流運動水平方向速度可忽略的條件下，式（1）可簡化為如下的一維方程［18］：

且滿足：

式中：z為豎直方向（m），向下為正；v(z)為冰川垂直于地表面方向的運動速度（m·s-1）。κ、H和tf分別為冰的熱擴散率（m2·s-1）、冰川厚度（m）和重建時間（s）。T(z，t)是在t時刻深度z處的瞬時溫度（℃），即由冰面溫度變化μ(t)導(dǎo)致的溫度擾動。

模型以穩(wěn)態(tài)冰溫垂直分布函數(shù)U(z)為初始條件。在初始時刻t= 0 的冰溫-深度剖面擾動（瞬時溫度）為零，即T(z，0) = 0。θ(z)對應(yīng)測量時刻t=tf的冰溫-深度剖面擾動，如圖1 所示?？梢钥闯?，μ(t)和θ(z)大于零代表冰面溫度較穩(wěn)態(tài)冰面溫度U(0)高，擾動為正；小于零則代表冰面溫度較U(0)低，擾動為負。

圖1 冰面溫度變化對穩(wěn)態(tài)冰溫垂直分布的影響Fig. 1 Influence of ice surface temperature change on the vertical distribution of steady-state ice temperature

穩(wěn)態(tài)冰溫垂直分布函數(shù)U(z)由冰川底部地?zé)崃髅芏萹（W·m-2）和長期平均冰面溫度Us（℃）共同確定：

式中：z為豎直方向（m），向下為正；v(z)為冰川垂直于地表面方向的運動速度（m·s-1）；U(0)為穩(wěn)態(tài)冰面溫度。κ、λ和H分別為冰的熱擴散率（m2·s-1）、冰的導(dǎo)熱率（W·m-1·℃-1）和冰川厚度（m）。

冰川底部地?zé)崃髅芏萹可看作是短期內(nèi)恒定不變的。在假設(shè)冰川穩(wěn)定的條件下，冰川厚度H不變；冰川垂直于地表面方向的運動速度v(z)與平均積累率相等；v(z)隨深度的增加呈線性減小，至冰川底部減小為零［13］?；跓醾鲗?dǎo)-冰流模型的氣候重建研究就是應(yīng)用相關(guān)的反演算法解決以下問題：已知式（2）在某一時刻的解T(z，tf) =θ(z)，求出邊界條件μ(t)。

需要注意的是，冰面溫度變化信號向下傳播的過程中，由于熱擴散的影響，其振幅隨著深度的增加而指數(shù)衰減、位相延遲。要利用冰川鉆孔溫度重建氣候，在較大深度的冰溫測量就需要很高的精度，一般不低于0.01 ℃［4］。實際中要達到這一精度和可靠測定難度較大。通常冰川鉆孔溫度由熱敏電阻溫度計測量，即通過測定冰層電阻，利用電阻-溫度轉(zhuǎn)換關(guān)系確定待測深度的溫度［19］。具體做法是：將按照一定長度間隔分布的若干個熱敏電阻探頭放置待測深度；待系統(tǒng)達到熱平衡狀態(tài)后，利用數(shù)字萬用表逐個測量對應(yīng)深度電阻傳感器的電阻值；電阻值經(jīng)計算程序轉(zhuǎn)換為對應(yīng)的溫度值。事實證明，已有的這種測溫系統(tǒng)在負溫條件下的分辨率可達0.01 ℃，取得了很好的結(jié)果［19-20］。

另外，在實際中應(yīng)用熱傳導(dǎo)-冰流模型有非常嚴格的適用條件，通常需要滿足以下幾點要求：（1）冰川底部呈負溫，與基巖凍結(jié)在一起，如極大陸型冰川；（2）冰面消融較少，內(nèi)部相變對鉆孔溫度的影響可忽略；（3）冰川鉆孔位于分冰嶺附近，水平方向的運動速度可忽略。由此可知，兩極和高緯度冷冰川分冰嶺附近的鉆孔溫度可以較好地用來重建氣候變化過程。

2 反演算法

2.1 最小二乘法

最小二乘法（最小平方法）通過最小化誤差的平方和尋找未知數(shù)據(jù)的最佳匹配，使得計算值與測量值之間的誤差達到最小，在多個學(xué)科領(lǐng)域具有廣泛的應(yīng)用［21-22］。假設(shè)Tc(z，tf)是由某一冰面溫度變化μ(t)計算的冰溫-深度剖面擾動，最小二乘法構(gòu)造了如下式的誤差函數(shù)J作為目標(biāo)函數(shù)。當(dāng)J取最小值時對應(yīng)的邊界條件μ(t)為最優(yōu)解［7］。

為了表示和計算方便，將μ(t)寫成如下傅里葉級數(shù)形式［15，18］：

任意給定估計系數(shù)a0、am和bm的初值，將μ(t)作為邊界條件代入式（2）計算Tc(z，tf)；并由式（5）求出J；a0、am和bm的迭代公式通過梯度下降法給出，如下式：

式中：n為迭代次數(shù)，γn為步長。直至J的取值達到給定的誤差水平停止上述迭代過程，找到最優(yōu)解。

不難計算，式（7）中的各項偏導(dǎo)值如下式：

式中：Wa0(z)、Wam(z)和Wbm(z)分別是當(dāng)μ(t)= 0.5，和時在tf時刻的冰溫-深度剖面擾動。將迭代過程中的式（8）代入式（7）最終可求出μ(t)。

2.2 Tikhonov正則化方法

已知熱傳導(dǎo)-冰流模型中的冰溫-深度剖面擾動θ(z)求解邊界條件μ(t)是一個典型的不適定問題。Tikhonov正則化方法常用來解決不適定問題。該方法是對式（5）的誤差函數(shù)J加一個約束條件；這樣的約束條件可以解釋為先驗信息。在該問題中，通常長期的冰面溫度變化μ(t)不會存在很大的波動，而由最小二乘法構(gòu)造的目標(biāo)函數(shù)對μ(t)沒有任何限制。因此，需要對μ(t)施加一個獨立的約束。Tikhonov正則化方法構(gòu)造了如下式的目標(biāo)函數(shù)Ψ；當(dāng)Ψ達到最小時，對應(yīng)的邊界條件μ(t)為最優(yōu)解［10，15］。

式中：α為正則化參數(shù)（α> 0）；Ω{μ(t)}稱為穩(wěn)定函數(shù)，其取值與μ(t)的光滑性有關(guān)。定義如下式：

求Ψ最小值的方法與最小二乘法相同，式（6）中估計系數(shù)a0、am和bm的迭代公式也通過梯度下降法給出，如下式：

式中各項偏導(dǎo)值如下式：

不難計算，穩(wěn)定函數(shù)Ω{μ(t)}可近似表示如下：

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

式（12）中各式最右邊一項分別為：

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

式（9）中的正則化參數(shù)α，相當(dāng)于對μ(t)添加了約束：α越大，約束越嚴格，μ(t)波動越??；α越小，約束越寬松，μ(t)波動越大。α調(diào)節(jié)著誤差函數(shù)J和穩(wěn)定函數(shù)Ω{μ(t)}之間的平衡。Tikhonov正則化方法不僅考慮了計算值與測量值之間的誤差水平；也考慮了μ(t)的波動大小，即μ(t)的穩(wěn)定性。可根據(jù)L-曲線準(zhǔn)則求出合適的α［23］。具體做法是在平面直角坐標(biāo)系中以(xi，yi) =(J，Ω{μ(t)})(i= 1，2，…，p)為坐標(biāo)點做出曲線（p為正則化參數(shù)的取值個數(shù)）。曲線通常呈L 形；該L 形曲線的拐角處對應(yīng)最佳的正則化參數(shù)。另外，考慮到J和Ω{μ(t)}的數(shù)量級問題，可對xi和yi兩者或其中之一取對數(shù)后再做出L-曲線；例如，在平面直角坐標(biāo)系中，坐標(biāo)點也可設(shè)置為(xi，yi) =(log(J)，log(Ω{μ(t)}))［23］。

2.3 蒙特卡羅算法

蒙特卡羅算法將冰面溫度變化μ(t)在不同時刻的取值μ→=(μ0，μ1，…，μtf)作為參數(shù)組合向量，隨機選擇μ→作為模型的邊界條件，采用隨機步的方法選取滿足給定誤差水平的參數(shù)進行統(tǒng)計分析［8，18］。具體的步驟如下：

步驟1：任意給定參數(shù)組合向量μ→n=作為模型邊界條件，計算誤差函數(shù)

式中：Tμ→n(z，tf)為由μ→n計算的冰溫-深度剖面擾動。

步驟2：在附近找到一個新的向量計算為了高效地尋找所有可能的參數(shù)組合，每一步都是隨機選擇中的某一個參數(shù)，對該參數(shù)添加微小擾動得到。

步驟4：以概率Pn接受新的向量：若接受，則；若不接受，則重復(fù)步驟1 到4。

通過以上步驟可以看出，蒙特卡羅算法是以產(chǎn)生一系列逐步改善的參數(shù)組合為基礎(chǔ)的。最后，選取滿足給定誤差水平的參數(shù)組合，運用統(tǒng)計學(xué)方法獲取各個參數(shù)取值μi(i= 0，1，…，tf)的頻數(shù)分布。這種分布通常存在區(qū)域極大值，對應(yīng)參數(shù)μi最有可能（最大概率）的取值。需要說明的是，通過統(tǒng)計方法得到的μi的頻數(shù)分布極大值可能不止一個；這反映了在對應(yīng)時刻存在多個溫度取值能夠與測量值吻合較好［4，18］。在這種情況下，由于實際中冰面溫度隨時間變化是平滑連續(xù)的，短期內(nèi)的溫度變化幅度不會太大。因此，可通過μi鄰近點的取值，找到合理的分布區(qū)域來確定μi。同時，為避免最終選取的參數(shù)組合之間存在相關(guān)性，也可對滿足誤差水平的μi進行等間距取值來統(tǒng)計頻數(shù)分布［4］。

3 反演算法的應(yīng)用和比較

3.1 構(gòu)造數(shù)據(jù)模擬實驗

為比較不同反演算法的重建結(jié)果，我們構(gòu)造了式（3）中的模型參數(shù)和冰面溫度變化μ(t)，進行理想條件下的數(shù)值實驗。假設(shè)式（3）中的各個參數(shù)取值為：冰川厚度H= 500 m，年均凈積累量0.3 m冰當(dāng)量厚度；重建時間tf= 400 a；冰的熱擴散率κ=1.3×10-6m2·s-1。冰面速度v(0) = 0.3 m·a-1，冰川垂直于地表面方向的運動速度v(z)隨深度的增加呈線性減小，至冰川底部減小為零。模型采用Crank-Nicolson 有限差分法求解［24］，它在時間方向上是隱式的二階方法，具有數(shù)值穩(wěn)定的特點。

首先，分別以冰面溫度變化μ1(t) = 2sin(2πt/tf)和μ2(t) = 2sin(πt/tf)作為邊界條件代入式（2），計算對應(yīng)的冰溫-深度剖面擾動θ1(z)和θ2(z)，如圖2 所示。然后，以θ1(z)和θ2(z)作為測量值重建冰面溫度變化。圖3 給出了在測量值為θ1(z)時不同反演算法的部分相關(guān)結(jié)果。最后，為驗證解的穩(wěn)定性，分別對θ1(z)和θ2(z)添加±0.02 ℃的隨機擾動并比較結(jié)果，對應(yīng)的重建結(jié)果如圖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)時不同反演算法的相關(guān)結(jié)果［a0、am和bm初值為0（a）和1（b）時最小二乘法重建的冰面溫度變化；Tikhonov正則化方法對應(yīng)的L-曲線（c）和重建的冰面溫度變化（0 < α < 0.01）（d）；蒙特卡羅算法不同時刻冰面溫度變化值的頻數(shù)分布（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 ℃隨機擾動前后的重建結(jié)果（a）、（b）；對θ2(z)添加±0.02 ℃隨機擾動前后的重建結(jié)果（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）由最小二乘法求得，分別是當(dāng)式（6）中的估計系數(shù)a0、am和bm初值為0［圖3（a）］和初值為1 時［圖3（b）］的重建結(jié)果。當(dāng)a0、am和bm初值為0時，重建結(jié)果與真實溫度變化基本一致；而當(dāng)a0、am和bm初值為1時，重建結(jié)果與真實溫度變化存在較大差異：近期冰面溫度變化波動明顯且振蕩幅度非常大。由此可知，應(yīng)用最小二乘法時a0、am和bm初值對重建結(jié)果有顯著影響。除此之外，梯度步長γ的選擇對重建結(jié)果影響也很大：步長太大會使估計系數(shù)在迭代過程中不斷增大，導(dǎo)致結(jié)果不收斂；而步長太小會使收斂速度太慢，可能需要迭代很多次，導(dǎo)致時間成本過高。因此，實際中往往需要設(shè)置一個合適的步長γ。實踐表明，本例中當(dāng)γ取0.0001 時，重建結(jié)果與真實溫度變化最接近［圖3（a）中紅色曲線］。圖3（b）是當(dāng)a0、am和bm初值為1，γ分別取0.001 和0.00001 時的重建結(jié)果。事實上，當(dāng)γ取0.01 時，無論a0、am和bm初值如何選擇，重建結(jié)果都趨于無窮大。

圖3（c）和圖3（d）是由Tikhonov正則化方法得到的L-曲線［其中坐標(biāo)點為（xi，yi）=（J，log（Ω｛μ（t）｝））］和正則化參數(shù)取值0＜α＜0.01的重建結(jié)果。其中，估計系數(shù)a0、am和bm初值為1。由式（9）可知，當(dāng)α取值為0 時，Tikhonov 正則化方法與最小二乘法的結(jié)果相同。與最小二乘法不同，應(yīng)用Tikhonov 正則化方法只需要給定合適的梯度步長（0.0001），可設(shè)定任意值作為a0、am和bm的初值，均不影響最終結(jié)果。本例中L-曲線拐角處x的取值為0.0012，對應(yīng)的α取值為0.00014。這里將α限定在［0，0.01］范圍內(nèi)，是因為在這個范圍內(nèi)的L-曲線上能夠找到合適的拐角［21］。當(dāng)α＞0.01 時，從結(jié)果上看，會使冰面溫度的變化幅度非常小；從L-曲線形狀上看，拐角右側(cè)的曲線會向著幾乎平行于x軸的方向繼續(xù)向右延伸，不存在其他的拐角。由圖3（d），冰面溫度變化整體上呈正弦函數(shù)形式，只是由于正則化參數(shù)取值不同，溫度變化振幅有所不同。隨著α取值的增大，冰面溫度變化幅度逐漸減小并趨于平滑和穩(wěn)定，即α對邊界條件起著施加約束的作用。可以看出，與最小二乘法相比，Tikhonov 正則化方法考慮了反問題中解的穩(wěn)定性問題：由正則化參數(shù)α調(diào)節(jié)著誤差函數(shù)和冰面溫度波動大小之間的平衡。這也在一定程度上限制了模型的復(fù)雜度。

圖3（e）和圖3（f）是由蒙特卡羅算法這一統(tǒng)計學(xué)方法獲取的不同時刻冰面溫度變化值的頻數(shù)分布，最終選取了滿足給定誤差水平的2 200 個參數(shù)組合進行統(tǒng)計分析。由于蒙特卡羅算法對隨機選擇的邊界條件沒有限制，即對模型的解沒有任何約束；在滿足誤差函數(shù)很小的條件下，往往會出現(xiàn)冰面溫度在短時間內(nèi)波動非常大的結(jié)果，失去了實際意義。為避免這種情況，本例中將參數(shù)組合向量=(μ0，μ1，…，μtf)各個分量的取值控制在±2 ℃的范圍內(nèi)進行模擬。由圖3（e），在t= 20 a和t= 140 a時刻，頻數(shù)分布極大值對應(yīng)的溫度變化值為0.6 ℃和1.5 ℃，是參數(shù)μi在對應(yīng)時刻最有可能的取值。事實上，這兩個時刻真實的溫度變化值為0.62 ℃和1.62 ℃，與計算值十分接近。結(jié)果表明，大多數(shù)時刻的溫度變化值都有著類似圖3（e）中的分布：能夠較容易地確定極大值，即選取的μi與最大頻數(shù)相對應(yīng)。而圖3（f）中，在t= 160 a時刻，頻數(shù)分布極大值對應(yīng)的溫度變化值分別為-1.8 ℃和1.1 ℃；其中，最大頻數(shù)對應(yīng)的取值是-1.8 ℃。一般情況下，根據(jù)冰面溫度變化的連續(xù)性，可通過鄰近時刻的頻數(shù)分布找到合理的分布區(qū)域，由該區(qū)域來確定極大值。在本例中，由于重建時間為百年尺度（400 a），長期平均冰面溫度變化的振幅不會太大。事實上，在t=160 a 鄰近時刻t= 120 a 至t= 180 a，頻數(shù)分布均集中在0～2 ℃這一區(qū)域。因此，t= 160 a 時刻選取1.1 ℃較-1.8 ℃更為合理。該時刻的真實溫度變化值為1.17 ℃，與1.1 ℃這一取值接近。需要指出的是，應(yīng)用蒙特卡羅算法有時很難找到合理的分布區(qū)域，重建結(jié)果會在短期內(nèi)存在明顯振蕩。在這種情況下，應(yīng)該選用其他的反演算法進行求解。

由圖4可以看出，最小二乘法、Tikhonov 正則化方法和蒙特卡羅算法的重建結(jié)果與真實溫度變化趨勢基本一致。其中，蒙特卡羅算法的結(jié)果不僅相對誤差較大，而且也不平滑。對θ1(z)和θ2(z)添加±0.02 ℃的隨機擾動后，Tikhonov 正則化方法的重建結(jié)果比最小二乘法和蒙特卡羅算法的結(jié)果更接近真實的溫度變化過程，解的穩(wěn)定性也最好。

3.2 真實鉆孔數(shù)據(jù)模擬實驗

真實鉆孔資料受實地環(huán)境影響。應(yīng)用熱傳導(dǎo)-冰流模型重建冰面溫度變化需要考慮鉆孔深度、年內(nèi)冰面溫度變化幅度、季節(jié)變化影響深度和式（4）中的穩(wěn)態(tài)冰溫垂直分布等實際問題。為與構(gòu)造的數(shù)據(jù)模擬實驗進行比較，我們以我國藏北高原腹地的馬蘭冰川鉆孔（位置點35°48.4′ N，90°45.3′ E；海拔5 680 m；鉆孔深度102 m）資料為例設(shè)定模型參數(shù)［25］。結(jié)合已有的、距離鉆孔位置最近的五道梁氣象站（35.13°N，93.05°E；海拔4 612 m；距離鉆孔位置約220 km）的觀測數(shù)據(jù)假定邊界條件，位置如圖5所示。最后，分別將假定的10 a 和40 a 冰面溫度變化代入式（2），計算對應(yīng)的冰溫-深度剖面擾動作為測量值反演邊界條件并比較結(jié)果。另外，為驗證解的穩(wěn)定性，對由40 a 冰面溫度變化計算的冰溫-深度剖面擾動添加±0.01 ℃的隨機擾動并對比結(jié)果。

圖5 藏北高原腹地馬蘭冰川鉆孔位置圖Fig. 5 Location map： the Malan Glacier on the northern Qinghai-Tibet Plateau showing the borehole site

馬蘭冰川位于我國藏北高原腹地的可可西里地區(qū)，最高峰海拔為6 056 m，是一個很大的典型山地冰帽冰川群；其南部平緩，冰面開闊平坦，屬于極大陸型冰川類型［26］。馬蘭冰川發(fā)育于寒冷干燥的氣候環(huán)境下，其嚴寒程度與東南極冰蓋邊緣地區(qū)相近。有關(guān)研究表明，位于藏北高原腹地的冰川變化幅度較邊緣山區(qū)要小，而自小冰期以來馬蘭冰川的變化非常小，表明了它處于比較穩(wěn)定的狀態(tài)［27］。

由于研究區(qū)域是地形復(fù)雜的山區(qū)，而五道梁氣象站距離馬蘭冰川鉆孔位置約220 km 且海拔相差較大。因此，氣溫觀測值與冰面溫度值相差較大。盡管如此，式（3）中的μ(t)并非冰面溫度的取值（絕對值），而是冰面溫度相對于穩(wěn)態(tài)冰面溫度U(0)的變化值（相對值）。因此，我們利用五道梁氣象站記錄的10 a 和40 a 氣溫變化幅度（趨勢）假定邊界條件。為驗證這一假定的準(zhǔn)確性，我們同時獲取了MODIS 地表溫度產(chǎn)品（MOD11A1）中鉆孔處像元在鉆取后一年（2000 年）的月平均地表溫度數(shù)據(jù)［數(shù)據(jù)來源：美國航空航天局（NASA）的地球科學(xué)數(shù)據(jù)系統(tǒng)Earthdata，https：//www.earthdata.nasa.gov/］與五道梁氣象站記錄的月平均氣溫變化幅度進行對比。

3.2.1 季節(jié)變化影響深度

馬蘭冰川鉆孔于1999 年5 月鉆取并測溫［25］。鉆孔底部未到冰床；20 m 深度以下冰溫變化幅度很小，僅為1.1 ℃。經(jīng)計算，冰川厚度H為130 m，年均凈積累量0.23 m 冰當(dāng)量厚度；冰床溫度-1.7 ℃。冰面豎直速度v(0) = 0.23 m·a-1。假設(shè)冰川垂直于地表面方向的運動速度v(z)隨深度的增加呈線性減小，至冰川底部減小為0。根據(jù)式（1）中各個物理參數(shù)取值的經(jīng)驗公式［13，26］：導(dǎo)熱率λ=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。經(jīng)驗證，本例中的導(dǎo)熱率和熱擴散率取值在上述范圍內(nèi)對實驗結(jié)果影響較小。因此，為便于計算，模型中取λ= 2 W·m-1·℃-1；κ= 1.4 × 10-6m2·s-1。由于各地的地?zé)崃髅芏萹差異很大，我們由馬蘭冰川鉆孔近底部的溫度梯度計算出地?zé)崃髅芏燃s為0.04 W·m-2。另外，根據(jù)《中國大陸地區(qū)大地?zé)崃鲾?shù)據(jù)匯編（第三版）》［30］，現(xiàn)有距馬蘭冰川鉆孔最近、具有地?zé)崃鲾?shù)據(jù)實測值的位置點是柴達木（38°14′ N，90°47′ E），如圖5所示。該位置點的校正熱流為0.043 W·m-2，與計算值接近。由上述冰床溫度和地?zé)崃髅芏韧扑闶剑?）中的長期平均冰面溫度Us為-4.2 ℃，由此計算穩(wěn)態(tài)的冰溫-深度剖面U(z)。

通常冰面以下十幾米（冰溫年變化深度）的鉆孔溫度反映的是年內(nèi)的季節(jié)變化。因此，重建長時間尺度的冰面溫度變化過程需要選取季節(jié)影響深度以下的冰溫-深度剖面。計算該深度需要估算一年內(nèi)的冰面溫度變化幅度。由于馬蘭冰川鉆孔的鉆取時間是1999年，我們參考五道梁氣象站前一年（1998 年）月平均氣溫的變化幅度［數(shù)據(jù)來源：國家氣象科學(xué)數(shù)據(jù)中心（CMDSC），http：//data. cma. cn］來假定一年內(nèi)冰面溫度變化這一邊界條件μ3(t)，由此計算季節(jié)影響深度，如圖6（a）所示。氣象資料表明，一年的氣溫變化幅度約為20 ℃（-15～5 ℃）。因為冰面溫度不會高于0 ℃，即夏季冰面溫度最高升至0 ℃。因此，我們假定一年的冰面溫度變化是取值范圍在-20～0 ℃的正弦函數(shù)，如圖6（b）所示。為驗證該冰面溫度變化幅度的準(zhǔn)確性，我們將鉆孔處像元2000 年的月平均地表溫度數(shù)據(jù)與其對比。結(jié)果表明，2000 年鉆孔處的最低月平均地表溫度約為-23 ℃，與假定的年內(nèi)冰面溫度變化（-20～0 ℃）基本一致，說明了該假設(shè)可靠。將其作為邊界條件代入式（2），得到冰溫-深度剖面擾動θ3(z)和擾動后的冰溫-深度剖面如圖6（c）和圖6（d）所示。

圖6 五道梁氣象站1998年的月平均氣溫（a），假定的一年內(nèi)冰面溫度變化（b），由（b）計算的冰溫-深度剖面擾動θ3(z)（c）以及穩(wěn)態(tài)與季節(jié)影響剖面（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)呈明顯的負向擾動。當(dāng)溫度剖面擾動小于0.1 ℃時，可認為在該深度上，年內(nèi)溫度變化的振幅消失［15］。經(jīng)計算，θ3(16)=-0.14 ℃，θ3(17)=-0.08 ℃，即季節(jié)影響深度為16 m（冰溫年變化深度）。

3.2.2 10 a尺度重建結(jié)果比較

在冰面消融微弱的條件下，可通過重建季節(jié)影響深度處的冰溫變化代替年均冰面溫度變化。由3.2.1 節(jié)可知，季節(jié)影響深度為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 左右。由于冰面溫度變化信號緩慢向下傳播，鉆孔越深處的溫度擾動對應(yīng)著更早期的冰面溫度變化。從深度100 m 向上至30 m，負向擾動逐漸增大；而30 m 以上的負向擾動逐漸減??；這說明了冰面溫度變化大致經(jīng)歷了先降溫，再升溫的過程。事實上，將θ4(z)由下而上觀察，可以看到與圖7（a）的冰面溫度變化相似的趨勢，即由θ4(z)可大致了解冰面溫度變化趨勢。但是，冰面溫度變化的具體時間和振幅需要通過相關(guān)的反演算法求解。

我們以θ4(z)為已知條件，分別應(yīng)用最小二乘法、Tikhonov 正則化方法和蒙特卡羅算法反演邊界條件，分析各個算法的重建結(jié)果并討論這些算法的優(yōu)劣性。圖7（c）給出了上述三種反演算法10 a 尺度的重建結(jié)果。整體結(jié)果與圖4 類似：最小二乘法和Tikhonov 正則化方法的重建結(jié)果較好；而蒙特卡羅算法的結(jié)果不連續(xù)，振蕩幅度大，準(zhǔn)確性也較低。

3.2.3 40 a尺度重建結(jié)果比較與穩(wěn)定性檢驗

與3.2.2 節(jié)的研究方法類似，這里我們結(jié)合五道梁氣象站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 尺度的重建結(jié)果。對θ5(z)添加±0.01 ℃的隨機擾動，重建結(jié)果如圖8（d）所示。由圖8 可知，最小二乘法和Tikhonov 正則化方法的重建結(jié)果與真實的溫度變化接近；Tikhonov 正則化方法的結(jié)果最平滑和穩(wěn)定；而蒙特卡羅算法結(jié)果的準(zhǔn)確性和穩(wěn)定性都較差。

將3.1 節(jié)構(gòu)造數(shù)據(jù)模擬實驗結(jié)果與3.2 節(jié)利用馬蘭冰川真實鉆孔數(shù)據(jù)模擬實驗結(jié)果進行對比，可得到以下一致的結(jié)論：最小二乘法、Tikhonov 正則化方法和蒙特卡羅算法這三種反演算法的重建結(jié)果與真實的冰面溫度變化基本一致；蒙特卡羅算法的誤差較大，且結(jié)果不平滑；Tikhonov 正則化方法是目前求解該反問題的最優(yōu)方法，能夠有效降低噪聲干擾，使得到的解穩(wěn)定。

需要說明的是，實際中馬蘭冰川的表面積雪消融和融水的再凍結(jié)等依然存在，對冰體和冰面水熱平衡有一定影響。而應(yīng)用熱傳導(dǎo)-冰流模型無法考慮這些因素的影響，這可能導(dǎo)致重建結(jié)果中近期的升溫幅度偏大［10］。盡管如此，為了更好地比較反演算法，在利用馬蘭冰川鉆孔數(shù)據(jù)模擬實驗中，冰面溫度變化是已知的；用于重建冰面溫度變化的冰溫-深度剖面擾動也是由給定的邊界條件代入模型計算得到的。因此，重建結(jié)果與真實溫度變化趨勢基本一致。此外，由于實際中無論測量精度多高，總會存在一定的測量誤差，由測量得到的冰溫-深度剖面相當(dāng)于添加了隨機擾動。為了能夠得到準(zhǔn)確的結(jié)果，可對測量的冰溫-深度剖面平滑后再進行反演。

4 結(jié)論

本文基于冰面溫度變化信號在冰層中的傳播特性和熱傳導(dǎo)-冰流模型，研究并比較了利用冰川鉆孔溫度重建冰面溫度變化這一反問題的反演算法，包括了最小二乘法、Tikhonov 正則化方法和蒙特卡羅算法。由于該反問題存在解的唯一性和穩(wěn)定性問題，不論是哪種反演算法，目標(biāo)都是找到相對穩(wěn)定的邊界條件使得冰溫-深度剖面的計算值最大程度地接近測量值，實質(zhì)是最優(yōu)化問題的求解。

由于蒙特卡羅算法隨機選擇邊界條件，再輸入模型進行計算和統(tǒng)計，因此，它可看作是正問題的求解方法。主要存在兩個問題：一是冰面溫度變化值的頻數(shù)分布會出現(xiàn)多個極值；二是求解的冰面溫度變化不平滑。盡管如此，因為蒙特卡羅算法是將模型的未知參數(shù)組合作為多維向量尋找最優(yōu)解，所以，在實際求解過程中，在部分參數(shù)（如初始條件等）未知的情況下，可選用蒙特卡羅算法這一求解方法。即考慮到能夠獲取的地學(xué)要素的限制，蒙特卡羅算法是最優(yōu)的。最小二乘法中估計系數(shù)a0、am和bm初值對重建結(jié)果影響很大，且僅考慮了誤差函數(shù)這一目標(biāo)函數(shù)，因此，它的穩(wěn)定性差。Tikhonov正則化方法將誤差函數(shù)和穩(wěn)定函數(shù)相加作為目標(biāo)函數(shù)，綜合考慮了解的準(zhǔn)確性和穩(wěn)定性問題，是目前求解該反問題的最優(yōu)方法。與其他反演算法相比，Tikhonov 正則化方法能夠有效降低噪聲干擾，使得到的解穩(wěn)定。