应用Monte Carlo反演方法研究地壳平均生热率
一、方法原理
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