Joint inversion of multi-station receiver functions and gravity data for imaging Moho variations and average crustal vP/vS ratios
-
摘要: 地壳厚度和波速比是研究地壳结构和组分的两个重要参数，可为区域构造研究提供重要约束。接收函数被广泛地用于确定地壳厚度和波速比，例如H-κ方法或H-κ-c方法，但是该方法只能确定台站下方的地壳厚度和速度比，当地震台站分布稀疏时，很难约束台站间的横向不均匀性。另一方面，重力数据也可用于莫霍面的起伏变化研究，它在横向上覆盖很好，有较高的分辨率，但在纵向上分辨率相对较低。为此，本研究提出了一种联合反演算法求解莫霍面深度和地壳波速比参数。联合反演算法综合考虑了接收函数在纵向上的较高分辨率和重力数据在横向上的较高分辨率，同时拟合区域内所有台站上的接收函数和区域重力数据。模型测试表明联合反演算法较单一的接收函数反演更精确，特别是对于地壳厚度的确定。Abstract: Crustal thickness and vP/vS ratio are two important parameters for understanding crustal structure and composition, which can help to study regional tectonics. Receiver function analysis has been widely used for determining crustal thickness and vP/vS ratio by the H-κ method or the H-κ-c method. However, it can only acquire average crustal thickness and vP/vS ratio beneath each seismic station, but cannot constrain their lateral variations among seisimic stations due to their sparse and irregular distribution. On the other hand, the gravity data has been widely used to derive the Moho variaitons, which has a good coverage and resolution laterally but poor resolution vertically. Therefore, in this study we have developed a new joint inversion method of receiver functions and gravity data to simultaneously invert for variations of Moho depths and average crustal vP/vS ratios in a region. The method takes advantage of complementary strengths of receiver functions and gravity data, and can simultaneously fit all receiver functions and gravity data in the region. The synthetic tests show that the proposed joint inversion method produces more reliable results than only receiver function analysis, especially for the crustal thickness.
-
Keywords:
- joint inversion /
- receiver function /
- gravity /
- Moho /
- vP/vS ratio
-
引言
地震地表破裂带指震源断层错动在地表产生的破裂和形变的总称，由地震断层、地震鼓包、地震裂缝、地震沟槽等组成（中华人民共和国国家质量监督检验检疫总局，中国国家标准化管理委员会，2012），中强地震产生的地震地表破裂带的分布范围、运动特征、位移分布、古地震特征等方面的研究对于中强地震特别是缺乏仪器测定地震参数的历史地震震中、震级及发震构造等的确定具有重要意义（聂宗笙等，2010；李文巧等，2011；李西等，2018）。
我国活动构造研究始于20世纪二三十年代（邓起东，闻学泽，2008），其中包括对广西灵山地震的考察与研究。关于1936年灵山M6¾ 地震的研究始于陈国达院士，他通过对震区地表破坏、房屋破坏、极震区房屋破坏及同震运动方向等的调查，认为极震区及发震断层呈NE向展布于罗阳山西北坡及山麓的高塘、鸦山岭、六鹅、夏塘和山鸡麓一带（陈国达，1939）。70年代以后，相关单位按照自身的任务对该地震再度进行调查和研究。陈恩民和黄咏茵（1984），李伟琦（1992）以及任镇寰等（1996）重新修正了极震区的烈度分布，增加了NNW方向的极震区长轴；潘建雄和黄日恒（1995），黄河生等（1990），周本刚等（2008），张沛全和李冰溯（2012）以及何军等（2012）对震区的断裂系统及断裂活动性进行了研究，其结果表明震区NE−ENE向的灵山断裂、NW向的友僚—蕉根坪断裂、NNW向的泗州断裂具有较强的活动性。
总体上看，前人对于灵山地震的研究主要集中于烈度分布和断裂活动性这两方面，而对于地震地表破裂带的研究多集中于地裂缝、崩塌、陷落等地表破坏现象调查（陈国达，1939；任镇寰等，1996）。关于灵山地震的地震参数及发震构造，前人也多从极震区烈度分布或断裂活动性来探讨，而从地震地表破裂带的空间展布及位移数据出发探讨此次地震的参数和发震构造则较少涉及。目前对于该地震的微观震中、宏观震中、震级等地震参数及发震构造均存在不确定性甚至是不同的认识（国家地震局全国地震烈度区划编图组，1979；陈恩民，黄咏茵，1984；黄河生等，1990；李伟琦，1992；潘建雄，黄日恒，1995；任镇寰等，1996；周本刚等，2008；张沛全，李冰溯，2012）。通过研究古地震地表破裂带讨论无仪器地震记录的历史地震的震中、震级、烈度等地震参数和发震构造已经在我国西部历史地震遗迹保留较好的地区有较多报道（聂宗笙等，2010；李文巧等，2011），但在华南中强地震构造区的陆域部分，由于地震和构造活动强度较弱、气候湿热以及地表人工改造较多等原因难以保留历史地震所产生的地震地表破裂带等遗迹，所以对于该地区缺乏仪器地震记录的历史中强地震的地震参数，多通过调查震区建筑物或地表破坏程度等所圈定的等震线图间接地获取，自然存在一定的不确定性。
鉴于此，本文拟在详细总结前人研究成果的基础上，从研究1936年灵山M6¾ 地震极震区的地表破裂特征入手，综合分析该地表破裂带的分布特征、运动性质、同震位移量等参数，并进一步讨论灵山地震的宏观震中、微观震中、震级、地震烈度，以期确定该地震的地震参数和发震构造。
1. 地震概况及区域地震构造背景
1936年4月1日9时31分，灵山县（震时属广东省）平山圩东罗阳山发生一次M6¾强烈地震，震中位于灵山县东北约20 km的罗阳山西北麓一带，震中烈度达Ⅸ度或Ⅸ度强。该地震造成101人死亡，263人受伤，8 000余间房屋倒塌破坏，是华南沿海地震带内陆地区自有地震记载以来发生的最大地震（国家地震局全国地震烈度区划编图组，1979；陈恩民，黄咏茵，1984；李伟琦，1992；任镇寰等，1996）。
灵山地震震区的大地构造位于扬子板块与华夏板块的结合部位，此次地震发生在区域性NW向巴马—博白断裂带与NE向防城—灵山断裂带的交会部位，罗阳山前的灵山断裂为震区断裂活动性最强的段落（黄河生等，1990；任镇寰等，1996；周本刚等，2008；何军等，2012） （图1）。灵山断裂沿罗阳山北西坡山麓展布，断裂南东侧罗阳山山地主要由印支期花岗岩组成，断裂北西侧钦江谷地主要由古生代、中生代沉积岩组成。
震区新构造位置处于广西新构造运动和地震活动均较强的桂东南断块区（李伟琦，1989）。新构造运动主要为间隙性抬升运动及断块差异运动。由于NE向防城—灵山断裂带的活动形成了NE向的断块隆起和断陷，造就了“两隆夹一陷”的新构造格局，断裂带南北两侧分别为罗阳山断块隆起和东山断块隆起，中间为钦江断陷谷地。由于间歇性抬升运动，在钦江谷地形成三至四级河流阶地，在罗阳山山前形成三至六级台地。在罗阳山前至少发育了五级台地或洪积残留台地（任镇寰，1983），灵山断裂切割了第二级洪积台地，在其上形成断层槽地或陡坎并右旋错移山前水系，周本刚等（2008）的研究表明其在晚更新世以来有明显活动，李细光等（2017a）的最新调查表明沿该断裂在高塘—夏塘—六蒙、蕉根坪—合口一带断续出露了1936年灵山地震地表破裂带（图1）。
![]()
图 1 1936年灵山M6¾地震震区的地震构造图左下角为区域构造简图。F1-2：灵山断裂；F2：蕉根坪—友僚断裂；F3：石塘断裂；F4：六银—寨圩断裂；F5：寨圩—浦北断裂；F6：佛子圩断裂；F7：泗州断裂Figure 1. Local seismic structure of the M6¾ Lingshan earthquake in 1936The inset in the bottom-left corner shows the regional geological structure. F1-2：Lingshan fault；F2：Jiaogenping-Youliao fault；F3：Shitang fault；F4：Liuyin-Zhaixu fault；F5：Zhaixu-Pubei fault；F6：Fozi fault；F7：Sizhou fault2. 地表破裂带与地震参数讨论
根据我们的最新调查与核实，1936年4月1日灵山M6¾地震的地表破裂带分布在罗阳山西北麓，沿着灵山断裂的高塘—六蒙、蕉根坪—合口一带形成了断续的地震裂缝、地震陡坎、断层槽地（图2a）、地震滑坡、冲沟错移等同震地表破裂现象，全长约12.5 km （图1），总体走为NE−ENE，同震右旋位移介于36—290 cm之间，同震垂直位移介于15—102 cm之间，同震最大右旋水平、垂直位移分别可达2.9 m和1.02 m （李细光等，2017b）。
整个地震地表破裂带的分布并不连续，有的段落上地震陡坎、地裂缝不发育，仅发育地震滑坡或断层槽地等，随着观察尺度的缩小，每一个地震陡坎均是由更次一级的阶梯状断层陡坎、地震裂缝等组成，总体运动形式表现为正-右旋走滑性质（图2b）。
![]()
图 2 高垌南东断层槽地和老陡坎（a）以及校椅麓南东的地震陡坎、地震裂缝和同震运动方向（b）Figure 2. Fault trough and old fault scarp in southeastern Gaodong （a） as well as earthquake scarp，fissures，old fault scarp and co-seismic motion direction in southeastern Xiaoyilu （b）2.1 震中位置
地震震中分为宏观震中和微观震中，其中宏观震中一般指极震区的几何中心（鄢家全等，2010），随着震级增大宏观震中与微观震中的差距可能愈加明显。特别是近年来发生的一系列大震如汶川地震改变了传统的宏观震中是一个点的认识，例如李志强等（2008）认为汶川地震的宏观震中是一条狭长的中间断开的线或窄带。
根据陈国达（1939）的考察结果，1936年灵山地震烈度最高、面积最小、同震地表破裂现象最明显的范围集中在罗阳山西北坡及山麓的高塘、鸦山岭、六鹅、夏塘、山鸡麓一带，并据此推断宏观震中位于罗阳山西北麓高塘—夏塘一线，呈一狭窄带状。陈恩民和黄咏茵（1984）以及李伟琦（1992）通过调查震区房屋破坏情况修正了极震区等震线的形状（图3c），陈恩民和黄咏茵（1984）认为本次地震的震中位于校椅麓附近。任镇寰等（1996）认为存在平山—蕉根坪和龙湾—高架岭这两个长轴走向分别为ENE和NNW的极震区，而宏观震中位于NNW向极震区长轴延长线与ENE向极震区长轴交点的高塘附近。已有研究显示地震地表破裂带与极震区分布范围具有很强的一致性（李志强等，2008；徐锡伟等，2008；马寅生等，2010），根据位错理论，断层每一段落所释放的能量与其错距平方成正比，地震断层上的最大位错点为初始破裂点，即宏观震中位置（张四昌，1989）。本研究调查显示，友僚—蕉根坪断裂以西的最大位移带位于夏塘水库东北至鸭子塘—蕉根坪一线附近，友僚—蕉根坪断裂以东发育两条地震地表破裂带（图3c），因此可以认为同震位错量为这两条地震地表破裂带位移量之和，即最大位移带位于蕉根坪至镇安一带，以此最大同震位移带为中心向NE和SW两个方向位移量呈递减的趋势。为此，我们在此最大位移带开挖了一系列探槽，其结果也揭示了地震断层、地震陡坎及地震崩积楔等丰富的地震地表破裂现象，据此推断本次地震的宏观震中极有可能位于灵山断裂北段与友僚—蕉根坪断裂交会处附近。
此前由于观测资料的缺失，本次地震没有可靠的微观仪器震中数据。据郭培兰等（2017）的最新研究成果，微观震中测定结果如图1所示，两种不同方法测定的震中位置分别位于灵家南东侧约30 km和灵家北西侧约15 km。由于灵山地震发生的时间较早，台站记录缺乏，且在灵山主震前不足5分钟发生了印度尼西亚卡拉克隆岛MW7.7地震，其记录强烈影响或覆盖了灵山地震的记录（郭培兰等，2017），依据微观震中位置在误差范围内与宏观震中位置基本保持一致，结合发震断层倾向SE，本文推断微观震中位置应在灵山断裂南东罗阳山一侧。
2.2 震级
由于灵山地震长期以来缺乏仪器记录，目前多数文献中均根据震中烈度Ⅸ度推算出灵山地震震级为M6¾ （李善邦，1960；顾功叙，1983；莫敬业，1990；中国地震局震害防御司，1999），郭培兰等（2017）根据仅有的上海徐家汇台记录图纸测定该地震的震级为MS7.0。但由于当时监测手段落后且记录单一，所测地震震级无法达到现代多台仪器所测定的地震震级精度。
鉴于上述原因，本文利用我们最新获得的地震地表破裂带参数，对比前人研究成果来推算1936年灵山地震震级，各种方法计算出的地震震级列于表1。
表 1 根据多种经验关系式推算的1936年灵山地震震级Table 1. Magnitude of 1936 Lingshan earthquake estimated according to severval empirical relationships between earthquake parameters and magnitude of historical earthquakes计算方法 公式 地震地表破裂
长度L/km最大同震地表
位移D/m震中烈度
Ie推算震级
M邓起东等（1992） M＝6.25＋0.8lgL 12.5 7.1 Wells和Coppersmith （1994） M＝5.16＋1.12lgL 12.5 6.4 陈达生（1984） M＝6.636 2＋0.565 1lgL 12.5 7.3 Wells和Coppersmith （1994） M＝6.81＋0.78lgD 2.9 7.2 李善邦（1960） M＝0.58Ie＋1.5 9.5 7.0 刘昌森（1989） M＝0.67Ie＋0.66 9.5 7.0 许卫晓等（2016） M＝0.549Ie＋1.859 9.5 7.1 根据陈恩民和黄咏茵（1984）、任镇寰等（1996）及李细光等（2017a）可知，极震区发育了长约12.5 km断续分布的地震地表破裂带以及丰富的滑坡、崩塌、砂土液化、地陷、地裂缝等同震地表破裂现象，所以我们认为震中烈度为Ⅸ度强并将其用于相关计算中。对比表1中结果可知，通过地表破裂带长度推算出的震级与通过震中烈度推算出的震级有较好的一致性，大部分均为M7左右，Wells和Coppersmith （1994）的计算结果较小，可能与我国采用的面波震级普遍大于国际上采用的矩震级且二者差值平均为0.3有关（戴志阳等，2008）。
由于1936年灵山M6¾ 地震的Ⅸ度强区域面积非常狭小，不到1平方千米（图3c），且灵山地区湿热多雨的环境和较为频繁的人类活动可能导致同震位移测量值偏大，因此，结合华南沿海地震带内陆地区的历史地震特征，我们推测1936年灵山地震的震级应在M6.8左右，这一结果与地表破裂带参数及极震区地表破坏现象吻合得较好。但值得注意的是，灵山M6¾ 地震的烈度衰减速度远大于其它地区发生的类似震级地震（陈国达，1939），所以通过烈度区的长轴半径来计算震级会出现计算结果偏小的现象。
2.3 震中烈度
前人对于灵山地震的烈度分布特别是震中（极震区）的烈度分布已作过很多研究（陈国达，1939；李伟琦，1992；任镇寰等，1996），由于评定烈度的标准不同，掌握的资料有异，所给出的地震烈度图和极震区形态也不尽相同（图3）。陈国达（1939）根据梅卡里（Mercalli）烈度表评定标准，将灵山地震震中的烈度评定为Ⅹ级；国家地震局全国地震烈度区划编图组（1979）汇编的 《中国地震等烈度线图集》 中灵山地震的震中烈度为Ⅸ度（图3a）， 《中国地震目录》 （顾功叙，1983）和 《中国近代地震目录》 （中国地震局震害防御司，1999）均采用了此图。此后陈恩民和黄咏茵（1984）、李伟琦（1992）以及任镇寰等（1996）均对震中区烈度进行了重新评定，震中烈度达Ⅸ度强，分别如图3b和图3c所示。
野外调查显示1936年灵山地震地表破裂带延伸至蕉根坪以东约3.5 km的合口和六蒙一带，与上述前人所绘制的极震区烈度图对比可见，地震地表破裂带分布范围与NE向Ⅸ度区的长轴方向基本一致，据此我们将灵山地震震中Ⅸ度区范围向东延伸至蕉根坪以东约4 km，沿地表破裂带两侧约1—2 km范围分布。
2.4 发震构造
对于灵山地震的发震构造，前人的认识不尽相同。陈国达（1939）认为，灵山地震震中为一长短轴比约5:1的狭长形椭圆，从“罗阳山脉西北麓一带水源断绝及水井干涸之原因推测，此次地震之发生，似即为该处原有断层继续活动之结果”，再结合烈度向南东衰减相对慢的特点推断此次地震的发震断层为罗阳山西北麓山前的SE倾向的NE向断裂。陈恩民和黄咏茵（1984）根据极震区的长轴方向认为，本次灵山地震的震源断裂面以ENE走向为主，NNW走向为辅；李伟琦（1992）根据极震区的等震线形状和低烈度区的长轴方向推测灵山地震可能是ENE向断裂与NNW向断裂共轭破裂的结果，其中ENE向构造是控震构造；潘建雄（1994）认为，在现代WNW−ESE向区域构造应力场作用下，NE向的防城—灵山断裂带和NW向的巴马—博白断裂带拟合为活动性较高的共轭构造，1936年灵山地震可能是这组共轭构造同时活动的结果，主破裂面以ENE向断裂为主；任镇寰等（1996）认为“ENE向、NNW向断裂均参与了本次地震的孕育过程，ENE向是主破裂”。
![]()
图 3 1936年灵山M6¾地震等烈度线图（a，b）以及极震区等烈度线与地震构造叠合图（c）Figure 3. Isoseismal contours of the 1936 Lingshan M6¾ earthquake （a，b） and its superposition with seismic structure （c）Fig. （a） is from Mapping Group of Seismic Intensity Zoning in China，State Seismological Bureau （1979），Fig. （b） is from Chen and Huang （1984）. In Fig. （c），the isoseismal contours are from Li （1992） and Ren et al （1996），and fault number is the same as Fig. 1罗阳山西北麓的高塘—合口、六蒙一带沿灵山断裂发育了长约12.5 km的地震地表破裂带，并且该断裂在地貌上表现为断层槽地、断层陡坎等；而NW向的蕉根坪断裂未有错断地表表现（图4），不是此次地震的发震构造。罗阳山南麓的泗州、根竹水、龙湾等地也有地裂缝，房屋破坏较为严重（陈恩民，黄咏茵，1984；李伟琦，1992；任震寰等，1996）。本研究调查发现罗阳山南麓的破坏以小规模滑坡为主，地裂缝等其它地表破坏的规模较罗阳山西北麓偏小，发育程度较其偏低，所以罗阳山南麓的地表破坏应为沿泗州断裂等NNW向小断裂的同震感应震动造成。综上认为，1936年灵山地震的发震构造为罗阳山西北麓NE−ENE走向的灵山断裂，南麓的NNW向、NW向断裂在主断层破裂影响下发生了感应震动，造成了局部烈度增强。
![]()
图 4 蕉根坪—友僚断裂野外照片及构造剖面图Figure 4. Outcrop and structural section of Jiaogenping−Youliao fault3. 讨论与结论
基于最新发现的1936年灵山地震地表破裂带的特征，本文分析探讨了此次地震的参数及发震构造，主要结论如下：
1) 在罗阳山西北麓，沿着灵山断裂的高塘—六蒙、蕉根坪—合口一带发育断续的以地震裂缝、地震陡坎、地震滑坡、冲沟错移等为标志的地震地表破裂带，全长约12.5 km，总体走向为NE−ENE，沿此地震地表破裂带同震右旋位移介于36—290 cm之间，同震垂直位移介于15—102 cm之间，最大同震水平、垂直位移分别为2.9 m和1.02 m，运动性质为正-右旋走滑。
2) 根据地震地表破裂带分布特征及极震区烈度分布判断灵山地震的宏观震中位于灵山断裂北段与友僚—蕉根坪断裂交会处附近，结合发震断层倾向SE，微观震中位置应在灵山断裂南东罗阳山一侧；综合多种方法认为1936年灵山地震震级为M6.8左右，目前仅有的单一仪器记录所测震级为MS7.0；前人调查结果显示震中烈度达Ⅸ度强，根据地震地表破裂带展布，将Ⅸ度区范围向东延伸至蕉根坪以东约4 km，沿地表破裂带两侧约1—2 km范围分布；罗阳山前灵山断裂构造地貌特征及受其控制的地表破裂带的发育说明灵山地震的发震构造为晚更新世以来活动的灵山断裂北段。
本文获得的同震位移量略大于其它地区相似震级地震所产生的地表位移（如1936年甘肃康乐M6¾ 地震），这一方面与历史地震参数的不确定性有关，另一方面1936年灵山地震的地表破裂带多发育在山麓地带，垂直位移可能受后期重力作用的影响偏大，而南方湿热多雨的环境可能会导致水平位移测量值偏大（如流水侧蚀作用较强）。此外，1936年灵山地震烈度的衰减速度远大于其它地区发生的类似震级的地震（陈国达，1939），这可能是导致地震地表破裂长度偏小的原因之一。
在本研究实施过程中，中国地震局地震研究所徐锡伟、冉永康、杨晓平、汪一鹏、宋方敏、田勤俭、何正勤、李伟琦等研究员在探槽选址、探槽编录解译、年龄样品采集等方面给予了悉心指导和帮助，桂林理工大学的陈磊、蒙南忠、李金峰参与了野外调查，审稿专家提出了重要的修改意见，作者在此一并表示衷心的感谢！
-
![]()
图 5 （a，b） 仅采用接收函数数据反演和联合反演获取的莫霍面；（c，d） 采用接收函数反演和联合反演所获取的莫霍面结果与理论模型的残差分布；（e，f） 采用接收函数和联合反演得到的地壳平均波速比；（g，h） 反演的平均地壳波速比结果与理论模型之间的残差
Figure 5. （a，b） The Moho results determined by receiver function analysis and joint inversion，respectively； （c，d） The deviations of inverted Moho models in Figs.（a） and （b） from theoretical Moho model in Fig. 1a，respectively；（e，f） The average crustal vP/vS ratios by receiver function analysis and joint inversion，respectively； （g，h） The deviations of inverted vP/vS models in Figs. （e） and （f） from the theoretical vP/vS model in Fig. 1b，respectively
![]()
图 1 莫霍面起伏模型（a）和地壳平均波速比模型（b），其中三角形表示虚拟台站
Figure 1. Synthetic Moho model （a） and the average crustal vP/vS model （b） where the black triangles indicate virtual seismic stations
![]()
图 2 虚拟台站（东向坐标x＝350 km，北向坐标y＝250 km）下方的简单P波速度结构（a）和根据该模型正演出的对应不同射线参数的理论接收函数（b），以及图1a中莫霍面起伏所引起的重力异常（c）
Figure 2. The simple crustal P-wave velocity structure beneath one virtual seismic station at x＝350 km in the east direction and y＝250 km in the north direction （a），and the theoretical receiver functions for different ray parameters （b），and gravity anomalies caused by Moho variations in Fig.1a （c）
![]()
图 3 联合反演L-曲线分析
（a，b，c）通过接收函数确定的不同平滑参数和阻尼参数下的归一化模型与数据残差的关系，最优参数$ {\varpi }_{H}＝300 $，$ {\lambda }_{H}＝300 $，$ {\lambda }_{k}＝8\;000 $；（d） 接收函数与重力之间权重关系曲线，重力的最优参数$ \gamma ＝25 $
Figure 3. L-curve analysis for the joint inversion
（a，b，c） Trade-off between the normalized model residuals and data residuals for different smoothing or damping parameters used in receiver function inversion （$ {\varpi }_{H}＝300 $，$ {\lambda }_{H}＝300 $，${\lambda }_{\kappa}＝8\;000$）； （d） Tradeoff analysis between the normalized model residuals and data residuals for different weights between receiver function and gravity data，and the optimal weight $ \gamma ＝25 $
![]()
图 4 （a） 接收函数的RMS迭代收敛曲线；（b） 重力异常的RMS迭代收敛曲线
图中黑色圆点为只采用接收函数到时数据，红色方块为采用接收函数和重力异常两种数据联合反演的结果
Figure 4. （a） The RMS residuals of receiver function；（b） The RMS residuals of gravity data
The black dots denote the results only by receiver function data，and the red diamonds denote the results by joint inversion
![]()
图 6 虚拟台站（x＝350 km，y＝250 km）下的地壳P波速度结构（a）和正演的接收函数（b）
Figure 6. The crustal P-wave velocity structure （a） at one virtual seismic station （x＝350 km，y＝250 km） and the theoretical receiver functions （b）
![]()
图 7 复杂速度模型情况下接收函数（a）和重力异常（b）的RMS迭代收敛曲线
图中黑色圆点为只采用接收函数到时数据，红色方块为采用接收函数和重力异常两种数据联合反演的结果
Figure 7. The RMS residuals of receiver function data （a） and gravity data （b） with iterations for the complex velocity model
The black dots denote the results only by receiver function data，and the red diamonds denote the results by joint inversion
![]()
图 8 只采用接收函数反演（a）和联合反演（b）得到的莫霍面结果与原始模型的残差以及只采用接收函数反演（c）和联合反演（d）获取的速度比结果与原始模型的残差
Figure 8. Deviations between theoretical Moho model in Fig. 1a and the inverted Moho models from only receiver functions （a） and joint inversion （b），and deviations between theoretical vP/vS model in Fig. 1b and the inverted vP/vS models from only receiver functions （c） and joint inversion （d），respectively
表 1 不同P波速度对联合反演的影响
Table 1 The effect of different P-wave velocities on joint inversion
vP
/（km·s−1）联合反演RMS拟合 联合反演残差 莫霍面深度/km vP/vS 接收函数/s 重力异常/（10−5 m·s−2） 最大残差 标准差 最大残差 标准差 6.0 0.362 2.905 0.664 0.285 0.177 0.041 6.1 0.361 2.842 0.607 0.263 0.164 0.033 6.2 0.363 2.794 0.554 0.251 0.152 0.033 6.3 0.363 2.743 0.526 0.249 0.139 0.041 6.4 0.363 2.695 0.569 0.256 0.126 0.054 6.5 0.362 2.649 0.612 0.268 0.123 0.062 6.6 0.364 2.602 0.653 0.286 0.142 0.068 
下载: 导出CSV
表 2 不同剩余密度参数对联合反演结果的影响
Table 2 Effect of different contrast densities on joint inversion
剩余密度
/（kg·m−3）联合反演RMS拟合 联合反演残差 莫霍面深度/km vP/vS 接收函数/s 重力异常/（10−5 m·s−2） 最大残差 标准差 最大残差 标准差 350 0.375 6.550 1.410 0.794 0.168 0.055 400 0.371 3.560 1.150 0.638 0.156 0.048 450 0.364 3.070 0.940 0.372 0.147 0.044 500 0.363 2.740 0.529 0.249 0.139 0.041 550 0.365 2.510 0.710 0.318 0.132 0.040 600 0.368 2.330 0.990 0.456 0.126 0.040 650 0.372 2.190 1.250 0.595 0.121 0.041
下载: 导出CSV
表 3 不同参考界面深度对联合反演影响
Table 3 Effect of different reference interfaces on joint inversion
参考深度/km 联合反演RMS拟合 联合反演残差 莫霍面深度/km vP/vS 接收函数/s 重力异常/（10−5 m·s−2） 最大残差 标准差 最大残差 标准差 38.5 0.363 2.745 1.908 1.410 0.136 0.076 39.0 0.363 2.745 1.439 0.953 0.121 0.063 39.5 0.363 2.745 0.971 0.515 0.130 0.051 40.0 0.363 2.743 0.529 0.249 0.139 0.041 40.5 0.363 2.742 0.994 0.545 0.148 0.034 41.0 0.363 2.742 1.463 0.986 0.156 0.032 41.5 0.363 2.743 1.931 1.444 0.165 0.034
下载: 导出CSV
-
冯锐. 1986. 三维物性分布的位场计算[J]. 地球物理学报,29(4):399–406. doi: 10.3321/j.issn:0001-5733.1986.04.010 Feng R. 1986. A computation method of potential field with three-dimensional density and magnetization distributions[J]. Acta Geophysica Sinica,29(4):399–406 (in Chinese).
郭良辉,孟小红,石磊,陈召曦. 2012. 优化滤波方法及其在中国大陆布格重力异常数据处理中的应用[J]. 地球物理学报,55(12):4078–4088. doi: 10.6038/j.issn.0001-5733.2012.12.020 Guo L H,Meng X H,Shi L,Chen Z X. 2012. Preferential filtering method and its application to Bouguer gravity anomaly of Chinese continent[J]. Chinese Journal of Geophysics,55(12):4078–4088 (in Chinese).
宋婷,沈旭章,梅秀苹. 2020. 利用接收函数频率特征研究莫霍面形态及应用[J]. 地震学报,42(2):135–150. doi: 10.11939/jass.20190149 Song T,Shen X Z,Mei X P. 2020. Constraining Moho characteristics with frequency-dependence of receiver function and its application[J]. Acta Seismologica Sinica,42(2):135–150 (in Chinese).
张洪双,田小波,滕吉文. 2009. 接收函数方法估计Moho倾斜地区的地壳速度比[J]. 地球物理学报,52(5):1243–1252. doi: 10.3969/j.issn.0001-5733.2009.05.013 Zhang H S,Tian X B,Teng J W. 2009. Estimation of crustal vP/vS with dipping Moho from receiver functions[J]. Chinese Journal of Geophysics,52(5):1243–1252 (in Chinese).
张盛,孟小红. 2013. 约束变密度界面反演方法[J]. 地球物理学进展,28(4):1714–1720. doi: 10.6038/pg20130411 Zhang S,Meng X H. 2013. Constraint interface inversion with variable density model[J]. Progress in Geophysics,28(4):1714–1720 (in Chinese).
张雪敏,付丽华,张海江,彭佳明. 2019. 基于正交秩-1矩阵追踪的天然地震数据重建研究:以加州San Jacinto断层密集地震台阵为例[J]. 地球物理学报,62(4):1427–1439. doi: 10.6038/cjg2019M0352 Zhang X M,Fu L H,Zhang H J,Peng J M. 2019. Reconstruction of natural earthquake data based on Orthogonal Rank-one Matrix Pursuit and its application to dense seismic array around the San Jacinto fault zone in California[J]. Chinese Journal of Geophysics,62(4):1427–1439 (in Chinese).
Aitken A R A. 2010. Moho geometry gravity inversion experiment (MoGGIE):A refined model of the Australian Moho,and its tectonic and isostatic implications[J]. Earth Planet Sci Lett,297(1/2):71–83.
Aitken A R A,Salmon M L,Kennett B L N. 2013. Australia’s Moho:A test of the usefulness of gravity modelling for the determination of Moho depth[J]. Tectonophysics,609:468–479. doi: 10.1016/j.tecto.2012.06.049
Aster R C, Borchers B, Thurber C H. 2013. Parameter Estimation and Inverse Problems[M]. 2nd ed. Burlington, MA, USA: Academic Press: 93–95.
Barbosa V C F,Silva J B C,Medeiros W E. 1999. Stable inversion of gravity anomalies of sedimentary basins with nonsmooth basement reliefs and arbitrary density contrast variations[J]. Geophysics,64(3):754–764. doi: 10.1190/1.1444585
Chai C P,Ammon C J,Maceira M,Herrmann R B. 2015. Inverting interpolated receiver functions with surface wave dispersion and gravity:Application to the western U.S. and adjacent Canada and Mexico[J]. Geophys Res Lett,42(11):4359–4366. doi: 10.1002/2015GL063733
Chakravarthi V,Sundararajan N. 2007. 3D gravity inversion of basement relief:A depth-dependent density approach[J]. Geophysics,72(2):I23–I32. doi: 10.1190/1.2431634
Cordell L,Henderson R G. 1968. Iterative three-dimensional solution of gravity anomaly data using a digital computer[J]. Geophysics,33(4):596–601. doi: 10.1190/1.1439955
Feng J,Meng X H,Chen Z X,Zhang S. 2014. Three-dimensional density interface inversion of gravity anomalies in the spectral domain[J]. J Geophys Eng,11(3):035001. doi: 10.1088/1742-2132/11/3/035001
Guo L H,Gao R,Shi L,Huang Z R,Ma Y W. 2019. Crustal thickness and Poisson’s ratios of South China revealed from joint inversion of receiver function and gravity data[J]. Earth Planet Sci Lett,510:142–152. doi: 10.1016/j.jpgl.2018.12.039
Han S C,Zhang H J,Xin H L,Shen W S,Yao H J. 2021. USTClitho2.0:Updated unified seismic tomography models for continental China lithosphere from joint inversion of body-wave arrival times and surface-wave dispersion data[J]. Seismol Res Lett,93(1):201–215.
He R Z,Shang X F,Yu C Q,Zhang H J,van der Hilst R D. 2014. Detailed Moho depth mapping of continental China by receiver function analysis[J]. Geophys J Int,199:1910–1918. doi: 10.1093/gji/ggu365
Hu S Q,Jiang X H,Zhu L P,Yao H J. 2019. Wavefield reconstruction of teleseismic receiver function with the stretching-and-squeezing interpolation method[J]. Seismol Res Lett,90(2A):716–726. doi: 10.1785/0220180197
Julià J,Mejía J. 2004. Thickness and vP/vS ratio variation in the Iberian crust[J]. Geophys J Int,156(1):59–72. doi: 10.1111/j.1365-246X.2004.02127.x
Li J T,Song X D,Wang P,Zhu L P. 2019. A generalized H-κ method with harmonic corrections on Ps and its crustal multiples in receiver functions[J]. J Geophys Res:Solid Earth,124(4):3782–3801. doi: 10.1029/2018JB016356
Li Y G,Oldenburg D W. 1998. 3-D inversion of gravity data[J]. Geophysics,63(1):109–119. doi: 10.1190/1.1444302
Li Y H,Gao M T,Wu Q J. 2014. Crustal thickness map of the Chinese mainland from teleseismic receiver functions[J]. Tectonophysics,611:51–60. doi: 10.1016/j.tecto.2013.11.019
Lowry A R,Pérez-Gussinyé M. 2011. The role of crustal quartz in controlling Cordilleran deformation[J]. Nature,471(7338):353–357. doi: 10.1038/nature09912
Maceira M,Ammon C J. 2009. Joint inversion of surface wave velocity and gravity observations and its application to central Asian basins shear velocity structure[J]. J Geophys Res:Solid Earth,114(B2):B02314.
Niu F L,Bravo T,Pavlis G,Vernon F,Rendon H,Bezada M,Levander A. 2007. Receiver function study of the crustal structure of the southeastern Caribbean plate boundary and Venezuela[J]. J Geophys Res:Solid Earth,112(B11):B11308. doi: 10.1029/2006JB004802
Oldenburg D W. 1974. The inversion and interpretation of gravity anomalies[J]. Geophysics,39(4):526–536. doi: 10.1190/1.1440444
Owens T J,Zandt G. 1997. Implications of crustal property variations for models of Tibetan Plateau evolution[J]. Nature,387(6628):37–43. doi: 10.1038/387037a0
Paige C C,Saunders M A. 1982. LSQR:An algorithm for sparse linear equations and sparse least squares[J]. CAM Trans Math Softw,8(1):43–71. doi: 10.1145/355984.355989
Parker R L. 1973. The rapid calculation of potential anomalies[J]. Geophys J R astr Soc,31(4):447–455. doi: 10.1111/j.1365-246X.1973.tb06513.x
Parker R L. 1974. Best bounds on density and depth from gravity data[J]. Geophysics,39(5):644–649. doi: 10.1190/1.1440454
Pavlis G L. 2011. Three-dimensional,wavefield imaging of broadband seismic array data[J]. Comput Geosci,37(8):1054–1066. doi: 10.1016/j.cageo.2010.11.015
Shi L,Guo L H,Ma Y W,Li Y H,Wang W L. 2018. Estimating crustal thickness and vP/vS ratio with joint constraints of receiver function and gravity data[J]. Geophys J Int,213(2):1334–1344. doi: 10.1093/gji/ggy062
Silva J B,Costa D C,Barbosa V C. 2006. Gravity inversion of basement relief and estimation of density contrast variation with depth[J]. Geophysics,71(5):J51–J58. doi: 10.1190/1.2236383
Song P H,Zhang X M,Liu Y S,Teng J W. 2017. Moho imaging based on receiver function analysis with teleseismic wavefield reconstruction:Application to South China[J]. Tectonophysics,718:118–131. doi: 10.1016/j.tecto.2017.05.031
Syracuse E M,Maceira M,Prieto G A,Zhang H J,Ammon C J. 2016. Multiple plates subducting beneath Colombia,as illuminated by seismicity and velocity from the joint inversion of seismic and gravity data[J]. Earth Planet Sci Lett,444:139–149. doi: 10.1016/j.jpgl.2016.03.050
Syracuse E M,Zhang H J,Maceira M. 2017. Joint inversion of seismic and gravity data for imaging seismic velocity structure of the crust and upper mantle beneath Utah,United States[J]. Tectonophysics,718:105–117. doi: 10.1016/j.tecto.2017.07.005
Tao K,Liu T Z,Ning J Y,Niu F L. 2014. Estimating sedimentary and crustal structure using wavefield continuation:Theory,techniques and applications[J]. Geophys J Int,197(1):443–457. doi: 10.1093/gji/ggt515
Tian X B,Zhang Z J. 2013. Bulk crustal properties in NE Tibet and their implications for deformation model[J]. Gondwana Res,24(2):548–559. doi: 10.1016/j.gr.2012.12.024
Zhang H J,Maceira M,Roux P,Thurber C. 2014. Joint inversion of body-wave arrival times and surface-wave dispersion for three-dimensional seismic structure around SAFOD[J]. Pure Appl Geophys,171(11):3013–3022. doi: 10.1007/s00024-014-0806-y
Zhang J H,Zheng T Y. 2015. Receiver function imaging with reconstructed wavefields from sparsely scattered stations[J]. Seismol Res Lett,86(1):165–172. doi: 10.1785/0220140028
Zhao Y,Guo L H,Guo Z,Chen Y J,Shi L,Li Y H. 2020. High resolution crustal model of SE Tibet from joint inversion of seismic P-wave travel-times and Bouguer gravity anomalies and its implication for the crustal channel flow[J]. Tectonophysics,792:228580. doi: 10.1016/j.tecto.2020.228580
Zhu L P,Kanamori H. 2000. Moho depth variation in southern California from teleseismic receiver functions[J]. J Geophys Res:Solid Earth,105(B2):2969–2980. doi: 10.1029/1999JB900322
-
期刊类型引用(3)
1. 郑礼洋，万永革. 广西及邻区构造应力场分区研究. 科学技术与工程. 2025(06): 2227-2236 . 
百度学术
2. 陈刚，李茂峰，李克华，唐勇，张清，马桂芳，申文豪，姜文亮. 基于随机振动模型的强地面运动数值模拟及其在广西灵山M6(3/4)地震中的应用. 中国地震. 2023(04): 832-844 . 百度学术
3. 潘黎黎，李细光. 广西灵山1936年6■级地震震中第四系地层研究. 地质科技情报. 2019(05): 147-154 . 百度学术
其他类型引用(0)
计量
- 文章访问数: 884
- HTML全文浏览量: 363
- PDF下载量: 267
- 被引次数: 3
