壹、方法原理
Monte Carlo反演方法的基本原理是：采用偽隨機數發生器產生壹系列先驗(a priori)模型，構成先驗模型空間；再計算各個先驗模型所對應的正演解；而後將正演解與觀測值相比，根據壹定的定量判別標準，決定該先驗模型的取舍，從而得到後驗(a posteriori)模型構成的後驗模型空間(Press，1968；Tarantola，1987)。
Monte Carlo反演方法的優點在於該方法的簡潔性。該方法不要求數據和反演參數之間存在線性關系，同時先驗模型和後驗模型的概率密度分布並不局限於Gauss分布或其他特定的概率分布。Monte Carlo反演方法唯壹所需的是存在解正演問題的算法，而且對於先驗和後驗模型空間中的模型采樣(sampling)在統計上是獨立的且具有代表性。只要Monte Carlo反演中先驗和後驗模型空間的采樣數量足夠大，反演結果將不會陷入局部最小。同時，Monte Carlo反演方法可以處理多峰態分布的參數反演問題。因此，Monte Carlo反演方法特別使用於地熱學反演，因為在地熱學中模型的非唯壹性問題非常突出。
Monte Carlo反演方法的主要缺點是計算量大。為達到統計上的代表性，需要在先驗模型空間中至少采樣幾千個正演模型，同時完成各個模型的正演計算。在後驗模型空間中新模型是否被接受的速率取決於所采用的隨機算法。因此，在先驗模型空間采樣過程中采用高效采樣算法，同時在後驗模型空間采樣中適當的采樣算法是十分重要的。
在先驗模型空間中采樣的最簡單的算法是搜索空間內的每壹個模型，但這樣做的代價是耗費機時，其低效是計算設備所無法接受的。因此，根據模型各參數的概率密度分布進行采樣是非常有效的算法。根據先驗模型計算的正演結果與觀測值之間的擬合差(misfit)，可以計算其最大似然估計(maximum likelihood estimate)。而在後驗模型空間中的采樣可以擬合差的最大似然估計值為準，采用Metropolis算法或模擬退火法進行采樣。
在我們的研究中，采用Mosegaard和Tarantola(1995)文章中的原理進行地熱學的Monte Carlo反演。對於壹個模型m而言，其含義是在模型空間M中的壹系列參數的組合，這些參數包括：熱導率和生熱率及其在深度上的變化。每個模型對應的先驗信息的概率密度為ρ(m)，而相應的後驗概率密度為σ(m)。後驗信息依賴於相應信息和似然函數L(m)：
σ(m)=kρ(m)L(m) (6)
在此，k是歸壹化系數，L(m)是正演計算值與觀測值之間的擬合差的度量函數。對於任壹給定的模型，其各個參數值均按同樣的先驗概率得到，而溫度和熱流值則根據熱傳導方程正演計算得出。這樣得到的大量參數組合及其相應的溫度和熱流值即代表了給定參數範圍內可能的巖石圈熱結構的統計結果。
在我們的研究中，設定觀測值的誤差擾動(noise)為Gauss分布。因此，可以采用觀測值和正演計算值的差值的平方和作為擬合參量S(mi)：
S(mi)=0.5[q(model)-q(measured)]2 (7)
式中：q(model)是正演計算得到的熱流值；q(measured)是觀測的熱流值(在本研究中僅處理壹維問題)。
模型空間的采樣采用Metropolis算法。該算法在采樣時是根據由先驗模型分布轉化為後驗模型分布的判據準則實現由壹個采樣點向下壹個采樣點的遷躍。在後驗模型空間中壹個模型的接受或拒絕則是根據其似然函數值的大小。我們定義先驗模型空間中兩個相鄰的模型mi和mj的似然函數分別為L(mi)和L(mj)。對應的擬合參量分別為S(mi)和S(mj)：
L(mi)=exp(- S(mi)/s2) (8)
L(mj)=exp(-S(mj)/s2) (9)
式中s為觀測值的標準偏差(即其誤差擾動)。當模型L(mj)≥L(mi)時接受模型mj為新的後驗模型。若L(mj)＜L(mi)，則按 P=L(mj)/L(mi)的概率接受模型mj為新的後驗模型。
如果模型mj被拒絕，則其前壹個模型(即模型mi)作為新的壹個模型加入後驗模型空間中。就這樣，上述算法將先驗模型空間中的各采樣模型轉化為後驗模型空間中的各個模型(Mosegaard和Tarantola，1995；Jokinen，2000)。
在本研究中，我們采用FORTRAN語言編寫程序，在PC機上實現上述算法。偽隨機發生器的算法據Press等(1988)的程序。溫度場的正演計算是通過求解溫度邊界條件的壹維穩態熱傳導方程實現。其中熱導率隨溫壓條件的變化通過叠代算法實現。采用地表平均熱流值做為擬合目標，觀測值的標準偏差(即其誤差擾動s)取3μW·m-3。
二、數值實驗
Jokinen(2000)采用上述Monte Carlo反演方法對北歐地盾區的巖石圈地熱結構進行了研究。當地的平均熱流值為46mW·m-2，地殼厚度為50km。反演結果表明，反演出的參數值與設計的正演模型之間不存在偏差，即可以應用Monte Carlo方法反演穩定地盾區的巖石圈熱結構及其相應的地熱學參數。
在本研究中，針對中國大陸東部地區熱流值明顯高於地盾區以及地殼厚度明顯薄於地盾區的實際情況，我們設計了地表熱流值為70mW·m-2，地殼厚度為30km的正演模型。該模型中上、中、下地殼的厚度各為10km，相應的生熱率分別為1.5μW·m-3，0.9μW·m-3和0.6μW·m-3，地殼平均生熱率為1.0μW·m-3，相應的地幔熱流值為40mW·m-2。常溫常壓下上、中、下地殼的熱導率λ0分別為3.0 Wm-1K-1，2.6 Wm-1K-1和2.6 Wm-1K-1(Chapman和Furlong，1992)；地殼熱導率隨溫壓變化的公式和參數也據Chapman和Furlong(1992)。地幔的生熱率為0.03μW·m-3，熱導率按Doin和Fleitout(1997)的公式取值；地表溫度取0℃。按這些參數進行的壹維穩態熱傳導正演計算求出64km處溫度為1100℃，78.5km處溫度為1300℃。
根據前述Monte Carlo反演程序，以地表0℃和78.5km處1300℃(模型1)或64km處1100℃(模型2)為邊界條件，采用地表平均熱流值作為擬合目標進行反演，求取地殼各層的生熱率和熱導率(地幔生熱率和熱導率與正演模型相同)。地殼各層的反演參數取值範圍見表4-9，其在反演時的抽樣采用隨機分布(evenly distribution)。反演中的先驗和後驗模型集合的數量均為10 000，反演結果見圖4-2和表4-10。
表4-9 反演實驗中熱導率和生熱率的取值範圍 Table4-9 The applied a priori values and ranges of thermal parameters used in the inversion experiments
註：熱導率λ0單位為Wm-1K-1；生熱率單位為μW·m-3。
表4-10 反演的地殼各層生熱率 Table4-10 The heat production rates of each layer from inversion experiments
註：生熱率單位μW·m-3。
圖4-2(a) 反演模型1的結果
Fig.4-2(a)The inversion results of Model 1
從圖4-2和表4-10可以看出，反演得到的Moho界面處的溫度值高於正演模型值，而地幔熱流值低於正演模型值，相應的地殼平均生熱率高於正演模型值；但反演結果中各參數的算術平均值與正演結果的差值均在反演值的標準偏差範圍內。究其原因，作者認為：在穩態熱傳導假設和取溫度邊界條件的情況下，基於變分原理，反演過程不可避免地導致地殼生熱率在允許的範圍內取最大值，以使穩態地溫線滿足極值條件，即在給定深度上的最大值。所以，按Monte Carlo方法反演得到的地殼平均生熱率代表地殼平均生熱率和Moho界面溫度值的上限估計值；相應地，反演得到的地幔熱流值代表穩態地幔熱流值的下限值。
圖4-2(b) 反演模型2的結果
Fig.4-2(b)The inversion results of Model 2
三、華北和華南地區地殼平均生熱率的Monte Carlo反演約束
根據以上研究結果，我們對華北和華南地區的主要構造單元，采用其地表平均熱流值作為擬合目標，正演計算的底界取大地電磁測深方法得到的巖石圈/軟流圈界面深度，對各構造單元的地殼生熱率進行Monte Carlo反演。
對各個構造單元而言，反演時上、中、下地殼的生熱率和常溫常壓下的熱導率的取值範圍及其中值與上述數值反演實驗的取值完全相同(參見表4-9)；而地殼內各層的厚度則根據人工地震剖面的資料取值(表4-11)。各構造單元的巖石圈/軟流圈界面深度按各單元的大地電磁測深結果的平均值取值，溫度壹律取1300℃；地表溫度取0℃。地幔的生熱率和熱導率取值亦與反演實驗相同。反演出的各構造單元的地殼平均生熱率和地幔熱流值見表4-11，反演的先驗和後驗模型的統計分布直方圖見圖4-3～圖4-14。
表4-11 華北及華南地區主要構造單元地殼平均生熱率的反演結果 Table4-11 The crustal heat production of tectonic domains in North and South China by Monte Carlo inversion
註：熱流值單位為mW·m-2；生熱率單位為μW·m-3。
圖4-3 鄂爾多斯盆地地殼生熱率、莫霍面溫度和地幔熱流值的Monte Carlo反演結果
Fig.4-3 Crustal heat production，temperature at Moho and mantle heat flow ofOrdos Basin by Monte Carlo inversion
圖4-4 華北盆地地殼生熱率、莫霍面溫度和地幔熱流值的Monte Carlo反演結果
Fig.4-4 Crustal heat production，temperature at Moho and mantle heat flow ofNorth China Basin by Monte Carlo inversion
圖4-5 河淮盆地地殼生熱率、莫霍面溫度和地幔熱流值的Monte Carlo反演結果
Fig.4-5 Crustal heat production，temperature at Moho and mantle heat flow of Hehuai Basin by Monte Carlo inversion
圖4-6 南陽盆地地殼生熱率、莫霍面溫度和地幔熱流值的Monte Carlo反演結果
Fig.4-6 Crustal heat production，temperature at Moho and mantle heat flow ofNanyang Basin by Monte Carlo inversion
圖4-7 蘇北盆地地殼生熱率、莫霍面溫度和地幔熱流值的Monte Carlo反演結果
Fig.4-7 Crustal heat production，temperature at Moho and mantle heat flow of Subei Basin by Monte Carlo inversion
圖4-8 東南沿海造山帶地殼生熱率、莫霍面溫度和地幔熱流值的Monte Carlo反演結果
Fig.4-8 Crustal heat production，temperature at Moho and mantle heat flow of Southeast Coast Orogenic Belt by Monte Carlo inversion
圖4-9 武夷山地區地殼生熱率、莫霍面溫度和地幔熱流值的Monte Carlo反演結果
Fig.4-9 Crustal heat production，temperature at Moho and mantle heat flow of Wuyi Area by Monte Carlo inversion
圖4-10 湘中地區地殼生熱率、莫霍面溫度和地幔熱流值的Monte Carlo反演結果
Fig.4-10 Crustal heat production，temperature at Moho and mantle heat flow ofCentral Yangtze by Monte Carlo inversion
圖4-11 四川盆地地殼生熱率、莫霍面溫度和地幔熱流值的Monte Carlo反演結果
Fig.4-11 Crustal heat production，temperature at Moho and mantle heat flow of Sichuan Basin by Monte Carlo inversion
圖4-12 康滇構造帶地殼生熱率、莫霍面溫度和地幔熱流值的Monte Carlo反演結果
Fig.4-12 Crustal heat production，temperature at Moho and mantle heat flow ofKangdian Tectonic Belt by Monte Carlo inversion
圖4-13 楚雄盆地地殼生熱率、莫霍面溫度和地幔熱流值的Monte Carlo反演結果
Fig.4-13 Crustal heat production，temperature at Moho and mantle heat flow ofChuxiong Basin by Monte Carlo inversion
圖4-14 蘭坪-思茅盆地地殼生熱率、莫霍面溫度和地幔熱流值的Monte Carlo反演結果
Fig.4-14 Crustal heat production，temperature at Moho and mantle heat flow ofLanping Simao Basin by Monte Carlo inversion
從表4-11可以看出華北和華南地區的地區平均生熱率的上限值為1.35μW·m-3。這同本章第三節得出的中國大陸主要構造單元地區平均生熱率不超過1.3μW·m-3的結論基本相符。
采用兩種方法得到的中國若幹主要構造單元地殼平均生熱率列於表4-12。從中可以看出，除少數外，采用Monte Carlo法反演得到的地殼平均生熱率在絕大多數構造單元都高於應用地下流體氦同位素比值與殼幔熱流比值關系得到的地殼平均生熱率數值。這說明前述Monte Carlo法反演結果至少代表當地地殼平均生熱率上限值的結論是可靠的。相反的情況僅在鄂爾多斯、四川和南陽盆地出現。其中，在鄂爾多斯盆地兩種方法的差值僅0.01μW·m-3，四川盆地僅差0.05μW·m-3。考慮到熱流數據和其他地球物理數據本身的誤差，可以認為此種差異並不影響我們對當地地殼平均生熱率所做估計的合理性。
值得指出，無論是本章第四節應用殼幔熱流比值和地下流體氦同位素的相關關系求取地殼平均生熱率，還是本節采用Monte Carlo反演方法計算地殼生熱率，或是本章第二節根據地盾區地幔熱流值估算的中國構造單元地殼生熱率的上限值，雖然都是利用大地熱流值資料進行，但是具體的方法卻各有不同。然而幾種方法得到的結果卻基本壹致，沒有嚴重的內在矛盾，說明大地熱流數據在涉及與放射性生熱元素有關的大陸地殼成分研究方面確實能夠起到其獨特的作用。由於在本章中我們所用的方法相對於傳統的區域地球化學方法研究深部地殼成分而言，對地震波速資料及其解釋的依賴性很小。因此，我們在本章中發展的應用大地熱流數據約束大陸地殼平均生熱率的研究方法，與基於地震波速的巖性解釋的傳統區域地球化學方法在技術方法層次上具有相當大的獨立性。所以，大地熱流資料和深部地熱學研究可以為地殼整體成分研究提供重要的信息，應當在地殼成分研究中加以充分關註。
表4-12 采用兩種方法得到的中國若幹主要構造單元地殼生熱率(A)的比較 Table4-12 Comparison of crustal heat production rates of somemajor tectonic units in China by twomethods