岩石圈的力学强度反映了岩石圈长期受构造应力作用的响应，是影响岩石圈变形演化的主要因素，表征岩石圈力学强度的最佳方式之一是通过计算岩石圈有效弹性厚度(effective elastic thickness of the lithosphere, T_e)来实现。T_e定义为与岩石圈板块中实际应力分布所产生的弯矩相等的弯曲弹性薄板的厚度^[1]。T_e概念源于区域均衡补偿模型^[2]或挠曲模型、弹性板模型。不同于局部补偿模型^[3-4](如Airy模型和Pratt模型)，区域均衡模型或挠曲模型其载荷部分是由覆盖在软流圈^[5]之上的岩石圈板内的弹性应力支撑并在一定区域范围内发生补偿。在挠曲模型中，岩石圈弹性板的响应往往使用其挠曲刚度(D)来表征。挠曲模型与Airy模型类似，均是通过岩石圈向下(或向上)发生挠曲以及地壳增厚(或减薄)来实现对地表地形负载的补偿。不同的是，作为局部补偿模式的Airy模型对应的挠曲刚度为0，即代表弱地壳；而挠曲模型则代表具有一定横向强度的地壳^[6]。

基于地形和重力异常数据求取T_e的导纳方法最早由Dorman等^[7]提出，但在早期利用导纳法(admittance)估计T_e时均是仅考虑地表的负载，而忽略地下负载加载的情况。实际上岩石圈不仅会受到地表负载的施压发生形变而产生山根，其同样可能受到地下负载的加载而产生上隆，例如洋中脊与地幔热物质的上涌会造成地壳和地幔物质的密度变化，从而形成地下负载^[8]。Forsyth^[1]首次提出基于地形和布格重力异常的相关函数法(coherence)，引入地表和地下负载的概念，并在导纳法中率先加入地下负载导纳的计算方法。Jin等^[9]在Forsyth基础上又提出准相关理论，该理论的关键是上下负载的相关角φ，当φ= 90°时，地表和地下负载正交，相关函数为传统的相关函数；当φ=0°时，地表和地下负载完全相关，相关函数为1；当0° < φ < 90°时，地表和地下负载具有不同的相关度。

由于采用传统傅里叶变换求取地形和重力异常功率谱方法(称为周期图法)会存在频率泄露问题，文献[10-11]又分别将多窗谱法引入导纳法和相关函数法以改善此问题。随后，Kirby等^[12]基于叠加的Morlet小波，提出一种Fan小波相关函数法来反演岩石圈的有效弹性厚度与初始负载比分布，该方法可克服采用固定尺度窗口而造成的数据频域信息不完整的缺陷^[13]。

华北地区位于欧亚大陆东南部，自中生代以来，相对稳定的华北克拉通遭受强烈的活化改造，发生大规模的构造变形、岩石圈减薄、岩浆活动和盆地再造，致使其原有的克拉通结构遭到明显的破坏^[14]，沉积层结构也相当复杂^[15](图 1)，因此华北克拉通被认为是世界上研究古老岩石圈减薄、改造、置换的最佳场所和天然实验室^[16]。近年来，国内外学者围绕华北克拉通破坏先后提出过不同的动力学模式，主要有热模式、化学侵蚀模式和拆沉模式等。然而无论哪种模式，都存在一些难以回答的问题，只能够解释部分地质现象^[17]。华北地区活动断裂丰富，地震频发，历史上曾发生过多次7级以上强震，因此认识华北地区的深部构造特征及岩石圈动力演化对于认识华北岩石圈力学强度和克拉通破坏具有重要的科学意义。部分学者利用不同的数据和方法对华北地区的岩石圈有效弹性厚度进行估算，如郑勇等^[18]采用Fan小波分析法计算华北克拉通T_e值及其各向异性；陈波等^[6]采用多窗谱法和Fan小波分析法计算中国及邻区T_e值，并分析其动力学意义；Lu等^[19]和侍文等^[20]分别基于导纳法和相关函数法，以ETOPO1地形数据、WGM2012布格重力异常数据和Crust1.0地壳结构数据为基础计算中国大陆地区T_e值分布，分别对中国大陆强震构造区力学特征和影响岩石圈强度的因素进行探讨。由于华北地区的地下结构非常复杂，前人使用的数据和方法存在差异，导致对华北地区T_e的计算结果不尽相同，同时部分学者计算时也未考虑沉积层的影响。

因此，本文使用精度更高的重力异常和地形数据，在完成对重力异常数据的沉积层校正后，基于导纳法和相关函数法联合反演并结合贝叶斯最优参数估计方法，得到华北地区的岩石圈有效弹性厚度以及初始负载比，讨论岩石圈有效弹性厚度与地震活动分布以及地表热流分布的关系，以期为认识华北岩石圈力学特征和克拉通破坏机制提供重要的重力学参考。

1 方法及数据 1.1 导纳法与相关函数法联合反演方法

导纳法的优点是物理意义明确、易于理解和实现；而相关函数法对不同均衡响应模式的转换波长更为敏感，计算精度更高^[20]。Audet^[21]提出结合导纳法和相关函数法的联合反演方法，该方法被国内外学者广泛采用并应用于T_e值的估算^[20-22]。

本文采用基于连续平面小波变换计算得到的导纳和相关函数进行联合反演，估计华北地区T_e值，同时为判断结果的准确性，结合贝叶斯推断法来评价特定情况下模型参数的统计分布。导纳和相关函数的表达式分别为^{[1, 21]}：

$ Z(k)=\frac{\Delta g(k)}{H(k)} $

(1)

$ r^2(k)=\frac{\left|\left\langle g(k) H^*(k)\right\rangle\right|^2}{\left\langle g(k) g^*(k)\right\rangle\left\langle H(k) H^*(k)\right\rangle} $

(2)

式中，Z(k)为导纳，Δg(k)为重力异常谱，H(k)为地形谱，k为波数，〈〉表示平均谱，*表示共轭。

为更准确地求取岩石圈有效弹性厚度，首先需要评价和判断负载来源，即首先考虑地表负载和地下负载同时加载的情况(图 2)。图中上部浅色区域代表地壳，深色区域代表地幔，H_i和W_i分别为初始地表和地下负载，H_T和W_T表示在地表负载作用下产生的地表和莫霍面变形，H_B和W_B表示在地下负载作用时产生的地表和莫霍面变形，H和W表示在共同作用下最终的地表和莫霍面变形。

采用发生挠曲变形前地表负载与地下负载的初始负载比F来判断负载的主要来源^[23-24]，初始负载比定义为：

$ F=\frac{f}{1+f} $

(3)

式中，f表示地下负载产生的高程谱H_B与地表地形产生的高程谱H_T之间的负载系数(loading factor)。当F=0时，f为无穷小，代表负载完全来自于地表；F=1时，f为无穷大，代表负载完全来自于地下。

使用连续平面小波变换分别计算重力异常谱和地形谱，得到每个网格节点的导纳和相干曲线。该方法的优点是可以避免将网格划分为一系列重叠窗口，从而大大降低频谱泄漏的概率。设定好一系列初始模型参数并计算其理论导纳和相关性之后，再使用非线性最小二乘法搜索与实测导纳和相关性拟合效果最佳的T_e和F值来实现对正演数据的反演。理论与实测导纳或相关函数误差使用如下卡方函数(Chi-squared function)进行计算^[21]：

$ \begin{gathered} \chi^2\left(F, T_e\right)= \\ \frac{1}{(2 N-2)} \sum\limits_{j=1}^2 \sum\limits_{i=1}^N\left(\frac{d_{i j}-s_{i j}\left(F, T_e\right)}{\sigma_{i j}}\right)^2 \end{gathered} $

(4)

式中，j表示实际导纳(j=1)或实测相关性(j=2)，d表示实测导纳或相关性，σ为理论与实测导纳或相关性函数之间的方差，s为给定模型的正演结果，N为采样网格数量。

1.2 贝叶斯推断

使用贝叶斯推断方法^[25]在每个采样网格对导纳法或相关函数法获得的有效弹性厚度(T_e)与载荷比(F)数据模型M(θ)进行评估，以确定该网格最优的由{T_e, F}组成的模型参数组合^[19-22]。贝叶斯推断方法表达式为：

$ \begin{gathered} P(\theta \mid x)= \\ \frac{P(x \mid \theta) P(\theta)}{P(x)}=\frac{P(x \mid \theta) P(\theta)}{\int_\theta P(x \mid \theta) P(\theta) \mathrm{d} \theta^{\prime}} \end{gathered} $

(5)

式中，x代表采样网格数据点构成的集合{x₁, x₂, …, x_n}；θ为采样网格数据点所对应的模型参数，由{T_e, F}组成；P(θ|x)为给定一组数据x后得到的模型参数θ的概率(即后验)；P(x|θ)为给定参数θ时生成数据x的概率(即似然)；P(θ)为独立于数据x的θ的分布(称为先验)；P(x)为归一化因子，以保证在整个模型参数空间上评估的后验概率为1。

从贝叶斯推断方法中随机选择3处区域来显示模型参数统计分布特征和曲线拟合效果，包括华北克拉通西部(109.2°E, 37.8°N)(图 3和图 4)、中部(114.1°E, 38.4°N)(图 5和图 6)以及东部地区(118.6°E, 37.5°N)(图 7和图 8)，使用贝叶斯推断法得到的模型参数统计分布以及采用不同反演方法得到的导纳和相关性拟合曲线如图 3~8所示，其中每处均使用4个独立的马尔科夫链(MCMC)从后验序列中抽取2 000个样本，并忽略前500个样本。通过对比相关函数法和联合反演法结果，以展示联合反演方法的优势。3个区域分别使用相关函数法和联合反演法得到的平均值、拟合误差和最大后验估计值见表 1。

从表 1以及图 4、图 6和图 8这3组计算结果可以看出，仅使用相关函数法进行反演时，得到的T_e和F值不确定性更大，最大后验估计MAP远大于平均值μ，T_e值高于联合反演法，标准差σ也明显大于联合反演法。从曲线拟合效果来看(图 3、图 5和图 7)，采用联合反演法得到的导纳曲线与相关性曲线拟合效果更佳，高于单独使用相关函数法的反演效果。因此，不论T_e值高低，相比于单独使用相关函数法，联合反演法均能获得更加准确的T_e和F结果。

1.3 重力数据

本文以华北地区(107°~123°E, 30°~42°N)为研究区，以最新的基于卫星、地面和航空等重力数据融合的高精度全球重力模型数据EIGEN-6C4^[26]为基础，采用导纳法和相关函数法进行联合反演，求取华北地区岩石圈有效弹性厚度T_e。EIGEN-6C4模型融合了GOCE卫星重力数据，在长波长部分较EGM2008模型数据具有更高的精度^[27-28]。地形数据来自斯克里普斯海洋研究所(SIO)的V19.1模型，空间分辨率为1′×1′。莫霍面深度数据和地壳密度数据来自CRUST1.0模型^[29]。使用谱方法进行计算时，需要地形和重力数据在长波长部分具备良好的相关性，但短波长的沉积层会降低地形的长波长成分，导致地形数据与重力异常之间的相关性降低，在计算T_e过程中会引入较大噪声，影响T_e值的准确求取。因此本文根据CRUST1.0模型中沉积层厚度数据，采用Parker^[30]正演法计算沉积层产生的重力异常，然后从布格重力异常中减去其结果，从而实现重力异常数据沉积层校正。

2 反演结果

图 9为反演得到的华北地区有效弹性厚度T_e与初始负载比F分布图。由图 9(a)可以看出，华北地区T_e值分布存在明显的横向差异性，整体呈现东西高、中间低，北部高、南部低特征。鄂尔多斯盆地T_e值较高，大约为30~60 km；华北平原T_e值总体较低，约为5~25 km；鲁东-黄海地块T_e值较高，约为20~50 km。初始负载比F(图 9(b))结果显示，华北地区东南部和太行山构造带北部F值较高(大于0.8)，鄂尔多斯地区F介于0.4和0.6之间，其他地区特别是华北平原地区整体F值较低(小于0.4)，表明在华北地区东南部和太行山构造带负载主要来自于地表，而在华北盆地地区负载主要来自于地下，这也与华北盆地岩石圈拉张减薄具有较好的对应关系。这种明显的不均匀性证实了华北克拉通东部发生克拉通破坏和岩石圈减薄，而西部仍保持相对稳定的特征。其原因可能是自中生代以来太平洋板块俯冲和停滞，地幔过渡带俯冲板片的停滞和脱水造成含水量增加，显著降低岩石圈强度，引起上地幔减压熔融或地幔物质向上流动。这些热物质上升后，不断通过热作用和化学作用从底部侵蚀岩石圈，从而逐渐导致华北克拉通岩石圈的减薄和破坏。

为考察图 9结果的准确性，计算有效弹性厚度T_e和初始负载比F的误差，结果如图 10所示。可以看出，两者误差分别介于0~5 km和0~0.12，且大部分区域的误差分别在0~3 km和0~0.06之内，表明本文得到的T_e和F结果均能较好地反映华北地区岩石圈力学强度及影响岩石圈力学强度的动力学来源。侍文等^[20]和Lu等^[19]采用联合反演计算得到华北地区T_e值在5~50 km之间，与本文计算结果基本一致，且高低值分布区也较为吻合，说明本文计算结果的可靠性较高。

3 结果分析 3.1 T_e分布与华北克拉通破坏

从图 9可以看出，华北克拉通西部鄂尔多斯地区表现为T_e高值区(30~60 km)，表明其岩石圈强度较大，构造相对稳定；F值介于0.4~0.6之间，表明鄂尔多斯地区的地表负载和地下负载大致相当。前人研究结果表明，鄂尔多斯地块上地幔呈现为明显的高速异常^[31]，说明鄂尔多斯地块仍然是一个未遭破坏、相对稳定的克拉通单元。南部秦岭-大别造山带T_e值较低，岩石圈活动性较强，秦岭-大别造山带东段F值较高，表明该地区岩石圈强度主要受地下负载影响。推测是扬子板块对华北地块的深向俯冲，强烈挤压隆起和碰撞造山等复杂构造演化，导致秦岭-大别造山带东段内部出现地壳垂向错断，并且秦岭-大别造山带东段与周围地块之间存在明显的深部构造转换带。在发生克拉通破坏时，东西向延伸的秦岭-大别-苏鲁局部增厚岩石圈将来自于地幔深处的热量转流到北部受扬子板块深俯冲作用影响较小的冀中、渤海湾以及松辽地区的岩石圈底部，使得这些区域的岩石圈成为减薄中心^[16]。

华北克拉通中部呈现为沿太行山构造带的T_e低值区(5~25 km)，通常T_e值较低表明岩石圈的力学强度较低，构造活动性相对较强。重力均衡研究结果表明，太行山构造带处于欠均衡状态^[32]。当岩石圈底部受到软流圈侵蚀、拆沉或机械拉张等破坏机制作用时，力学强度较低的地区将优先被破坏或破坏程度较为显著。北部燕山构造带地壳内普遍存在高、低速互层现象，表现出大规模岩浆活动特征^[33]；地震成像研究表明，燕山地区壳幔边界强烈的速度差异可能来自下地壳和岩石圈拆沉诱发的软流圈抬升，从而形成速度显著不同的物质分界面^[34]。古太平洋的西向俯冲作用导致华北克拉通下方产生不稳定的地幔流动体系，并引起上地幔减压熔融或地幔物质向上流动^[35]，强烈的对流、橄榄岩-熔体相互作用及伸展作用也可能导致岩石圈软化和克拉通破坏^[33]。这与该地区低T_e值(高活动性)和高F值(地表与地下负载共同作用)特征相吻合。

华北克拉通东部为T_e中、低值区(5~35 km)。郯庐断裂带起源于华北与华南板块的碰撞造山过程，郯庐断裂带较低的岩石圈强度及有利的岩浆活动通道，使其在克拉通破坏中成为伸展活动带和岩浆活动带，是华北克拉通破坏最严重的区域^{[34, 36]}。郯庐断裂带全线及其次级断层成为软流圈物质上涌的首选通道，这些热物质上升后，不断通过热作用和化学作用从底部侵蚀岩石圈，从而逐渐导致华北克拉通岩石圈的减薄和破坏。

3.2 T_e分布与地震活动性关系

部分学者认为，岩石圈的强度能够控制岩石圈在构造力作用下的局部变形，因此与地震分布具有一定关联性^[37-38]。从图 9中T_e值与地震空间分布来看，华北地区地震主要沿T_e值较低(5~20 km)或高低值过渡带区域分布，如岩石圈强度较低的山西断裂带、太行山构造带和渤海湾盆地北部等地区。而在T_e值较高区，如华北克拉通西部的鄂尔多斯地块内部，地震则较少发生。这表明在T_e值较高区其岩石圈强度也较高，而T_e值较低区岩石圈强度也较低，在构造应力作用下更易于释放积累的应变，从而更易诱发地震。由于岩石圈有效弹性厚度并不完全代表岩石圈的真实厚度，地震活动性与T_e的关系也较为复杂，T_e低值区并不一定能够代表地震频发区，地震活动性主要与其所处的构造环境有关。如太行山造山带南部和秦岭-大别造山带T_e值虽然较小，但很少发生地震。秦岭-大别造山带形成于扬子板块和华北克拉通碰撞造山运动，区域内发育大规模的岩浆和构造变形活动^[39-40]，表明该区域岩石圈强度的降低主要是由于深部热物质上升和流动过程中发生强烈的壳幔改造作用，发生地震所需的能量被释放，导致该区域地震活动较少。在高T_e值的黄海地区，地震活动却很频繁，这是因为江苏陆区和南黄海海域的基底同属于华南下扬子地块^[41]，嘉山-响水断裂和苏州-湖州断裂等大型断裂带从陆地延伸至海域，作为主要的边界断裂控制着该区域的地震活动特征。此外，地震分布与F值大小之间并未发现明显的对应关系，因此初始负载比F对孕震环境的影响作用尚不明确。

3.3 T_e与地表热流关系

地球内部的热作用是地质构造运动和地震活动的主要动力来源之一。大地热流包含丰富的岩石圈热结构和地球动力学信息^[42]。T_e值大小还受岩石圈热结构影响，因此T_e与地表热流分布之间存在特定联系^[6]。图 11为根据中国大陆地区地表热流数据库^[43]采用克里金插值方法绘制的华北地区地表热流分布图，图中圆点代表地表热流观测点。可以看出，华北地区地表热流数据观测点较少，且分布不均匀。华北地区热流值整体较高，一般来说，高地表热流反映岩石圈构造较为活跃，代表岩石圈强度较弱，即T_e与地表热流值呈现负相关关系。例如低地表热流值和高T_e值的鄂尔多斯地区，构造活动较为稳定，代表稳定的克拉通区域；而高地表热流值和低T_e值的山西断裂带与华北平原，构造活动相对活跃。但是在苏北盆地和渤海湾地区却表现出高热流值和高T_e值，这可能是因为受到华北克拉通破坏的影响。在太平洋板块俯冲以及后撤作用下，中国东部构造体制由早期的挤压转换成伸展^[44]，在该背景之下，深部地幔上隆、岩石圈伸展、大量热物质上涌带来巨大的热量，使整个地区产生活化作用，导致苏北盆地和渤海湾地区地表热流值较高。

4 结语

本文使用EIGEN-6C4模型布格重力异常、V19.1地形和CRUST1.0模型地壳结构数据，采用导纳和相关函数联合反演方法，得到华北地区岩石圈有效弹性厚度T_e和初始负载比F分布，讨论T_e与华北克拉通破坏、地震活动和地表热流的关系，得到以下结论：

1) 华北克拉通T_e值分布整体呈现横向不均匀性，这种明显的不均匀性证实了东部发生克拉通破坏和岩石圈减薄，而西部仍保持相对稳定的特征。

2) 从地震空间分布来看，地震活动与T_e值关系复杂。华北地区地震活动主要沿T_e相对低值区发生，表明在岩石圈强度薄弱区容易释放积累的应力与能量，更有利于地震的孕育和发生。而T_e高值区构造则相对稳定，整体具有较强的抵抗岩石圈变形的能力。地震是构造应力长期积累和突然释放的结果，地震活动性主要受其所处的深部构造环境影响。

3) T_e值大小还受岩石圈热结构的影响，T_e值与地表热流值整体对应关系较好，但在苏北盆地和渤海湾地区，T_e值与地表热流值呈现负相关关系。其原因可能是由于热物质上升和流动，浅部地壳被破坏，产生高地表热流值，而深部地壳仍具有克拉通属性。