相比于其他地球物理场(重力场、磁场等)，地球深部应力场无法进行直接测量，使得构建准确、高时空分辨率的应力场模型十分困难。地幔对流被认为是影响中下部岩石圈应力场的重要机制，其能够在重力观测中得到反映，故基于重力场推导岩石圈应力场具有一定的可行性。Runcorn^[1-2]在只考虑重力对流的条件下，使用2~8阶扰动位系数，对全球大尺度地幔对流应力场进行初步估计；Liu等^[3-5]在此基础上计算太平洋板块、亚洲、非洲等地区的地幔对流应力场，提出13~25阶扰动位系数对应的是小尺度地幔对流应力场。国内学者同样使用低阶扰动位系数确定了区域地幔对流应力场^[6-7]。

解算地幔对流应力场时，扰动位系数阶数的确定十分关键。在以往研究中，对于该参数的确定常基于经验判断，缺乏定量准则，为避免该问题，国内学者使用其他方法计算地幔对流应力场^[8-9]。但对于全球尺度的地幔对流应力场的研究，Runcorn模型仅使用扰动位阶数进行解算，因此在大尺度地幔对流应力场研究中仍具有优越性。本文结合点质量源法确定截止阶数，并基于EGM2008模型2~213阶扰动位系数，解算深度为30 km的全球尺度地幔对流应力场。此外，由于板块运动是地幔对流的表现形式之一，震源较深的地震(大于等于20 km)成因与地幔对流有关。而前人研究^[2]同样认为，深部地震是地幔对流应力场存在的实际现象，因此本文通过比较2~213阶扰动位系数解算的地幔对流应力场与全球范围5.5级以上地震的震中分布和PB2002板块划分模型，讨论全球尺度地幔对流应力场与地震震源及板块划分模型之间的联系。

1 地幔对流应力场计算中截止阶数选择

Runcorn^[2]在假定地幔是不可压缩流体，地幔上边界速度为0、无下边界、粘滞系数为常数、只考虑密度引起的地幔对流的条件下，给出地幔对流应力场方程：

$ \left\{\begin{aligned} \sigma_\theta= & \frac{M g}{4 {\rm{ \mathsf{ π}}} a^2} \sum\limits_{n=2}^{N_{\max }}\left(\frac{R}{a}\right)^{n+1} \frac{2 n+1}{n+1} \frac{\partial}{\partial \theta} \\ & \sum\limits_{m=0}^n\left(C_{n m} \cos m \lambda+S_{n m} \sin m \lambda\right) P_{n m}(\cos \theta) \\ \sigma_\lambda= & \frac{M g}{4 {\rm{ \mathsf{ π}}} a^2} \sum\limits_{n=2}^{N_{\max }}\left(\frac{R}{a}\right)^{n+1} \frac{2 n+1}{n+1} \frac{1}{\sin \theta} \frac{\partial}{\partial \lambda} \\ & \sum\limits_{m=0}^n\left(C_{m m} \cos m \lambda+S_{m m} \sin m \lambda\right) P_{m n}(\cos \theta) \\ \sigma= & \sqrt{\sigma_\theta^2+\sigma_\lambda^2} \end{aligned}\right. $

(1)

式中，σ_θ、σ_λ、σ分别为地幔对流应力沿经、纬度方向的大小及二者合力的大小，M为地球质量，g为地球平均重力，R为地球半径，a为地球中心到计算深度之间的距离，C_nm、S_nm为扰动位系数，P_nm(cosθ)为勒让德函数。由于R＞a导致地幔对流应力在阶数趋近于无穷时，σ_θ、σ_λ、σ会趋近无穷大，地幔对流应力函数发散，因此为式(1)选择一个合适的阶数是必要的。

在Runcorn模型基础上对地幔对流应力场的研究中，前人对于扰动位系数阶数的选择缺乏定量准则，因此如何确定扰动位系数阶数与研究深度之间的关系是非常重要的。点质量源法建立了扰动位系数与场源深度之间的对应关系，本文结合点质量源法与Runcorn模型，通过扰动位系数阶数、地球半径与场源深度之间的对应关系，初步确定了Runcorn模型的截止阶数。

点质量源法原理如下：将球体分为K层，每层厚度为h，将每层内的质量压缩为一个质面，密度异常在卫星重力测定中可以视为点质量源。使用点质量源法，建立重力扰动位系数阶数与异常埋深体之间的关系。位于点质量源的扰动位可以表示为：

$ U=\frac{G m}{Z} $

(2)

式中，Z为场源深度。大地水准面和扰动位之间的关系可以表示为：

$ N=\frac{U}{g_0}=\frac{G m}{g_0 Z} $

(3)

式中，g₀为地表平均重力。点质量元的垂直分量可以表示为：

$ g=\frac{G m}{Z^2} $

(4)

联立式(3)和式(4)，得到重力异常和大地水准面之间的比率为：

$ Z=\frac{g_0 N}{g} $

(5)

由于Bowin等^[10]给出了每一阶的重力和大地水准面高的关系：

$ \frac{g_n}{N_n}=\frac{g_0(n+1)}{g} $

(6)

联立式(5)和式(6)，将g_n和N_n近似于一个点质量元产生的g和N，可以将结果简化，并得到场源深度Z和地球半径R之间的关系：

$ Z=\frac{R}{(n+1)} \frac{g_n}{N_n} \frac{N}{g}=\frac{R}{n-1} $

(7)

使用式(7)限定Runcorn模型中的截止阶次。将式(1)与式(7)联立得：

$ \left\{\begin{aligned} \sigma_\theta= & \frac{M g}{4 {\rm{ \mathsf{ π}}} a^2} \sum\limits_{n=2}^{[R / Z]+1}\left(\frac{R}{a}\right)^{n+1} \frac{2 n+1}{n+1} \frac{\partial}{\partial \theta} \\ & \sum\limits_{m=0}^n\left(C_{n m} \cos m \lambda+S_{n m} \sin m \lambda\right) P_{n m}(\cos \theta) \\ \sigma_\lambda= & \frac{M g}{4 {\rm{ \mathsf{ π}}} a^2} \sum\limits_{n=2}^{[R / Z]+1}\left(\frac{R}{a}\right)^{n+1} \frac{2 n+1}{n+1} \frac{1}{\sin \theta} \frac{\partial}{\partial \lambda} \\ & \sum\limits_{m=0}^n\left(C_{n m} \cos m \lambda+S_{n m} \sin m \lambda\right) P_{m n}(\cos \theta) \\ \sigma= & \sqrt{\sigma_\theta^2+\sigma_\lambda^2} \end{aligned}\right. $

(8)

本文引入点质量源法中场源深度与扰动位系数阶次的关系方程，初步对Runcorn模型截止阶次进行限定(式(8))。后文采用式(8)中[R/Z]+1计算截止阶数，并利用式(8)计算地幔对流应力场。

2 地幔对流应力场结果与讨论 2.1 全球地幔对流应力与地震分布的关系

板块内部的应力场变化，是导致地震发生的重要原因之一，而地幔对流应力场是板块内部大尺度应力场的起源之一，所以研究地震震中位置与地幔对流应力场之间的空间分布关系是有必要的。前人的研究工作中也将地震震中位置、板块边界与地幔对流应力场汇聚、拉张等高应力区域进行了比较^[6]。

本文使用全球地震目录(2015~2023年深度大于30 km)和板块划分模型(PB2002)与解算结果进行比较。地震数据来源于Seismological Facility for the Advancement of Geoscience提供的地震目录(https://ds.iris.edu/wilber3/find_event)。Runcorn^[2]选择30 km作为地幔对流应力场计算深度，本文亦选择30 km作为计算深度，将该深度代入式(8)中[R/Z]+1后，得到截止阶数为213阶。由于EGM2008模型^[11]应用范围广，模型精度高，本文基于EGM2008重力模型使用式(8)对地幔对流应力场进行计算。

现有结果(图 1(a)，紫红色五角星表示地震)表明，在太平洋东部、南美洲西部智利、秘鲁及邻近海域、南美洲东部、非洲板块东北部与欧亚板块西南部边界、澳大利亚板块南部与南极板块东北部交汇区域存在较大应力区域。在北美南缘可可斯-加勒比板块边界(区域①)、南美洲西部南美-纳兹卡板块边界(区域②)及非洲板块东北部与欧亚板块西南部边界部分，现有结果与全球范围深度大于30 km的地震震中位置对应较好。本文解算结果(图 1(b))表明：除现有结果中高应力部分外，在环太平洋板块区域同样存在高应力区域。与此相对应的是，在日本海沟和阿留申海沟(区域③)、欧亚板块-菲律宾板块边界(区域④)、太平洋板块西北缘菲律宾-太平洋板块交汇处马里亚纳海沟(区域⑤)、欧亚板块南缘爪哇海沟(区域⑥)及太平洋板块西南缘汤加海沟与克马德克海沟(区域⑦)区域，与全球范围内深度大于30 km的地震震中位置对应效果较好。

结果表明，二者差异(图 1(c)和1(d))主要集中在区域①、区域②、区域③、区域④、区域⑤、区域⑥及区域⑦，这些区域都存在大量汇聚或拉张的地幔对流矢量，且绝大部分垂直于板块边界，这些区域的地幔对流可能受到小尺度地幔对流的影响。可以注意到，在阿拉伯-欧亚板块交汇处及喜马拉雅山脉(图 1白框)区域，相较于其他区域而言，地震震中分布更加稀疏，且该地区的地幔对流应力场较小，推测可能与该地区地壳厚度较高、地震主要来源于地壳有关。

2.2 典型区域地幔对流分析

为进一步证明地幔对流与地震震中位置的关联性，选择2015~2023年5.5级以上地震数据，将网格划分为1°×1°，对全球地震数据在对应网格的发震频数进行统计，省略频数在5以下的网格后，将对应数据投影至世界海岸线及板块划分模型PB2002上。可以看到，地震频率较高的位置与图 2(a)(绿色三角为地震频数为5~10的网格，黄色三角为地震频数为10~15的网格，绿色三角为地震频数为15以上的网格)中强震位置较为符合，主要包括欧亚板-菲律宾板块边界(A区域)、南美-纳兹卡板块边界(B区域)和澳大利亚-太平洋边界(C区域)。

本文将这些发震较频繁的典型区域与地幔应力矢量分布进行比较。图 2(b)为欧亚-菲律宾板块边界，图中地震分布与沿板块边界区域的高应力区域显示出较强的关联性。该区域北部边界为琉球海沟，其附近菲律宾板块NW向俯冲至亚欧板块，与本文解算结果的地幔应力矢量方向一致；中南部为马尼拉海沟，其附近地幔应力矢量为NE向，与该地区存在浮力高原且板块俯冲作用减弱的构造背景吻合；南部为菲律宾海沟，其附近地幔应力为SW向，表明该地区菲律宾板块向SW向俯冲并与欧亚板块碰撞，与谭皓原等^[12]的结果一致。

图 2(c)为南美-纳兹卡板块边界部分，同样可以看到，地震主要集中于南美板块，且地震分布与沿板块边界的高应力区域较为吻合。该地区南美洲板块西缘地幔应力为W向；纳兹卡板块北部海岸线边界地幔应力矢量为W向，且随海岸线以西至秘鲁海沟应力逐渐减小，秘鲁海沟以西逐渐增大，西部在智利海沟附近地幔对流应力呈汇聚状态，该地区纳兹卡板块俯冲至南美洲板块下的构造环境，与前人使用地震数据模拟的俯冲带构造结果一致^[13]。

澳大利亚-太平洋板块(图 2(d))中西北部俾斯麦海附近及西南部北斐济海盆地震频发区域与高应力部分较为吻合。总体来看，西北部新几内亚海沟附近的地幔对流应力矢量方向为SW向，与板块移动方向相符，而东南部分北斐济海盆西南边缘地幔对流应力矢量方向较为混乱，但海盆中部存在扩张情况^[14]。从细节来看，地幔对流应力矢量方向在西北部新几内亚海沟附近、东南部北斐济海盆西边界接近垂直，而在北斐济海盆东南边界则为近乎平行状态，可能与该边缘为断裂带有关。

3 结语

本文引入点质量源法中场源深度与扰动位球谐系数的对应关系，初步限定了Runcorn模型的截止阶次。基于该思路，使用EGM2008地球重力场2~213阶扰动位系数解算深度为30 km的地幔对流应力场方程，解算结果整体趋势与2~30阶解算结果相同，且能突出更多小尺度高应力区域。对本文结果进行分析，可得到以下认识：

1) 本文将地幔对流应力场结果与2015~2023年全球范围内5.5级以上地震震中位置进行比较，结果表明，板块边界附近的地震可能受地幔对流应力场控制，但在亚洲中部及青藏高原区域对应关系较差，可能与该地区地壳厚度较高、地震主要来源于地壳有关。

2) 欧亚-菲律宾板块边界、南美-纳兹卡板块边界、澳大利亚-太平洋板块边界，近8 a内深度大于30 km的地震(M≥5.5)频数较高。对这3个地区的地幔对流应力与构造环境进行进一步分析，结果与前人研究结论相符，可以认为地幔对流应力场在板块边界部分能够较好地反映该地区的构造状态。

3) 相比于前人2~8阶解算的应力场结果，本文2~213阶的解算结果能够呈现更多的地幔对流应力场细节，与太平洋板块西南缘、非洲中南部、印度洋中部的地震重合度更高，推测这些地区的地震可能受小尺度地幔对流控制。