贝叶斯同化重力反演方法构建龙门山地壳密度模型
发布时间：2021-05-15
摘要重力异常对地壳横向密度变化敏感，而无约束重力反演得到的密度模型其垂向分辨能力往往不理想.为了改善反演结果的垂向分辨率，本文参考已有先验分层模型，基于贝叶斯原理，提出了一种重震联合反演的新策略，可实现多种参考模型和复杂加权参数条件下的最大后验概率估计.理论模型测试结果表明，对于深度加权、多参考模型约束等多种问题，本文提出的新方法都可以稳健地获得最优化的模型参数.本文同时以中国地震科学台阵在龙门山地区及周边的一维接收函数分层模型和地震层析成像结果为参考，通过此方法对该区的重力异常进行反演，获得了该区的高精度三维密度结构，其水平分辨率优于10km，垂直分辨率优于5km.结合四条通过汶川和芦山地震震中的剖面进行分析后发现，反演得到的密度结构模型在过强震震源区位置横向变形显著，其揭示的分层地壳结构和变形模式与地表已知断裂构造具有相关性.本文提出的重震联合反演新策略，可为研究潜在强震风险源区的地壳结构和物性特征提供有效的科技方法支撑.
关键词重力反演;贝叶斯推断;三维密度结构;龙门山构造;数据同化
0引言
地球物理反演技术是研究地壳内部结构和性质的重要工具.由于地球物理反演问题先天具有多解性，通常反演结果在数学上适定，但可能与实际地质情况并不相符(LiandOldenburg，1998;Williamsetal.，2004;Howe，2009;Dirianetal.，2016).现今在很多研究热点地区，利用多种手段获得的地壳结构模型和认识越来越多.如何充分融合已有的模型，发挥各种地球物理手段的优势，减小反演结果的不确定性，无疑是地球物理联合反演研究要解决的核心问题(Welfordetal.，2018;Yangetal.，2012).
如果将所有已有模型解释结果看作先验认识，当今地学家面临的首要挑战是如何从已有模型数据中提取尽可能多的有用信息以获取新的见解.与常规技术相比，数据同化提供了一套新模式，用于发现常规技术不容易揭示的数据和模型结构之间的关系(Bergenetal.，2019).数据同化将观测数据、计算模型和先验推断集成到一个系统中，通过对系统中不确定性的估计和控制，来提高后验估计模型的精度.在当今这个大数据时代中，数据同化方法已成功地应用到了多个科学领域的建模、数据分析和预测中(Riggelsenetal.，2007;deWitetal.，2013).
重力反演是探测地球深部构造的有效手段之一，具有对地下物质密度变化敏感、水平分辨能力强、成本低等优势(Keareyetal.，2002;Hinzeetal.，2013).然而，由于观测数据误差和核函数算子本身的性质，重力反演具有多解性(Green，1975;FediandRapolla，1999).除了通过增加观测数据和减小观测误差，在一定程度上降低反演的多解性.之外，在反演中加入先验信息进行约束是更有效的方法(LiandOldenburg，1998;陈石等，2010;李红蕾等，2016).近年来，重力约束条件越来越受到重视，许多地球物理学者都提出了各种各样的约束条件以及相应的反演方法.约束总体上可以分为两类，即数学结构约束和参考模型约束.前者包括最小构造约束(LastandKubik，1983)、深度加权约束(LiandOldenburg，1998)、平坦度和光滑度约束(BoulangerandChouteau，2001)等，后者主要包括地震波速转化密度约束(王新胜等，2013)、地质约束(Nabighianetal.，2005)、岩性约束(Gengetal.，2019)等.在传统的三维重力反演算法中，上述约束和参考都是通过阻尼因子被引入到反演方程中，并通过广义交叉验证(GCV)(Golubetal.，1979)或L曲线准则(Hansen，2001)获取的单个阻尼因子实现对观测数据和先验信息权重的平衡.然而，当在先验假设(光滑先验、深度加权)和参考权重设置增多的情况下，传统算法很难依据数据特征实现对多个约束信息权重参数的优选，这很容易造成反演结果偏离实际(LiandOldenburg，1998，2000;Williams，2008).
Akaike(1977)最早利用贝叶斯方法来解决气象领域的数据同化问题.这种方法已成为数据同化的重要手段，其是以数据驱动模式量化反演参数，通过最小化Akaike贝叶斯信息量(ABIC)实现对观测数据及先验约束不确定性的估计与控制，从而显著提高模型估计精度，减少模型不确定性.同样方法随后成功应用到时-空传染型余震序列模型(ETAS模型)的参数估计(OgataandKatsura，1993)、三维地震层析成像、应力图像反演(TerakawaandMatsu′ura，2008)、GPS断层数据反演(Fukahataetal.，2004;FukahataandWright，2008)、重力平差参数优化(Chenatal.，2019)等多个地球物理学领域.
本文探索通过数据同化手段，将贝叶斯同化策略引入到传统的重力约束反演中，设计一种可以融合多种先验推断的重力异常贝叶斯同化反演算法，以期得到更加准确的最大后验概率意义下的模型参数估计.文中第2部分给出了三维重力贝叶斯同化反演的求解方程和ABIC目标函数的构建;第3部分设计了两个综合实验模型，测试了该方法在解决多种已知先验参考约束和深度加权约束中超参数和后验模型估计问题的优化效果，验证了该方法的可行性和有效性;第4部分和第5部分是应用该方法和龙门山及周边地区2′×2′的高精度WGM2012布格重力异常数据反演了该区的地壳密度模型，并对龙门山地区地壳内低密度分布、隆升变形机制、强震震源区环境等进行了深入分析;最后，第6部分对本文研究方法和认识进行了总结和讨论.
2实验测试
为测试上述方法的可行性和有效性，我们分别设计了两组仿真测试模型，来检验贝叶斯同化反演策略的实际效果.
2.1参考模型约束
引入适当的参考模型约束能在垂向分层上改善重力异常反演结果.通常认为参考模型越准确，反演结果越接近真实情况.但在实际问题中，可参考的模型不止一个，其误差、不确定性差异都不尽相同，同时不同参考模型之间也往往不一致.对于多参考模型约束问题，特别是在先验参考模型误差未知情况下，如何分配每个参考模型的权重问题一直是研究热点(Hansen，1994;Vogel，1996).下面就本文提出的方法如何分配权重和改善反演结果进行测试.
这个数值实验中，真实模型在y方向没有变化，在y=0处的垂直切片，模型起伏界面上方蓝色部分密度为0.1g·cm-3，下方红色部分密度为0.2g·cm-3，起伏界面上下密度差为0.1g·cm-3(如图1a).图1b是仿真的重力异常数据，在理论真实模型正演异常的基础上，加入了5%的高斯误差.
在测试中，设计了两个简单的参考模型，每个模型与真实模型都存在一定的误差.其中，参考模型M1密度起伏界分界面形态与真实模型一致，但界面两侧的密度差比真实模型小10%，如图1c所示;参考模型M2中的密度起伏界面位置与真实模型不一致，界面两侧的密度差真实模型大10%，如图1d所示.此外，参考模型1和2中还分别加入了2%和1%高斯误差.
在模型测试过程中，我们选择了四种不同的反演方案，方案1是仅有参考模型M1优化约束反演;方案2是仅有参考约束模型M2优化约束反演;方案3是参考模型M1和M2联合非优化约束反演;方案4是参考模型M1和M2联合优化约束反演.每种反演都实现了模型异常与数据异常的拟合，但反演的模型结果差异性显著，详见图2所示.
从图2中的结果可以看出，图2a的反演模型相比图2b结果界面更清晰，但上下密度差偏差较大，图2b的界面凸起部分效果比图2a差.图2c的联合反演结果，相比前两个无论在界面起伏和上下密度差方面都有一定程度改善，而与图2d的贝叶斯同化联合反演结果相比，在上下界面位置和密度差方面图2d明显优于图2c.
在该模型测试中，参考模型1和2的归一化权重参数分别为λ1和λ2.图2a—d的反演结果，对应的参数详见表1.从表1中的结果可以看出，在图2d的最优化过程中，参考模型1的权重更大，对应的ABIC值也最小.
相关期刊推荐：《地球物理学报》创刊于1948年，主要刊载固体地球物理、应用地球物理、地磁和空间物理、大气和海洋地球物理，以及与地球物理密切相关的交叉学科研究成果的高质量论文。作者和读者对象主要为从事地球物理学、地球科学及其他相关学科的国内外科技工作者和大专院校师生。
结合表1中的反演统计参数结果来看，参考模型1约束反演优化得到的超参数λ1为129.055，以此计算得到的ABIC值为-1355.069;参考模型2约束反演优化得到的超参数λ2为160.062，计算得到的ABIC值为-1406.770;将上述两个模型单独优化得到的超参数直接代入多模型约束，计算得到的ABIC值为-3217.861;将两个模型同时优化反演得到的最优λ1和λ2分别为1798.5和1094.8，最小化ABIC值为-3700.195.
将表1反演统计参数与图2中的反演结果结合来看，与单个参考模型和非优化权重的多个参考模型约束反演相比，ABIC最小化(ABIC最小值为-3700.195)给出的最优超参数可有效降低反演的卡方值，提高反演的准确性和精度.与手动指定的模型约束权重进行反演计算得到的反演结果(图2c)相比，基于优化ABIC的反演算法可依据观测数据来选择提取不同参考之间的有用信息，并实现不同参考模型与反演结果之间的数据同化.
2.2深度加权约束
在式(1)的核函数矩阵G中，每个元素数值大小都与观测点和模型单元之间的距离平方成正比，在不同深度位置上，模型单元的核函数大小数值差异明显，浅层单元的数值远大于深层.体现在反演结果中，常常会出现异常的变化仅集中在浅部模型单元上，通常称之为重力位场反演的趋肤效应(LiandOldenburg，1998).
3龙门山地壳模型
3.1龙门山构造背景
龙门山地区位于南北地震带中南部位，拥有复杂的构造边界条件、活动断层系统.其西部是活跃的青藏高原边界，东部是稳定的华南地台，是青藏高原主体向东挤出构造边界，地壳变形强烈，结构复杂.地块内部褶皱断裂广泛发育，其中包括了多个大地构造区边界断裂和控制强震活动的活断层：鲜水河断裂带、龙门山断裂带、东昆仑断裂带、龙泉山断裂带、龙日坝断裂带、毛尔盖分支断层等(徐锡伟等，2008).此外，沿着龙门山断裂带，还存在着约5km的强烈高程梯度带(如图4所示).
同时，该地区也属于中国地震科学实验场区，过去几十年中以中国地震科学台阵项目为代表的大规模地球物理观测相继开展，已积累了大量的地球物理探测工作成果(如Yaoetal.，2008;Lietal.，2011;Zhangetal.，2011;王绪本等，2018).这些深部成果对该区的壳幔结构及其相关动力学问题达成了部分共识，但在一些基本问题上依旧存在争议，如在该区的地壳上地幔变形机制的解释上，就有壳幔连续变形机制(EnglandandMolnar，1997)、块体挤出机制(Tapponnieretal.，1982，2001)、下地壳流机制(ClarkandRoyden，2000;Roydenetal.，1997)等多种模式
.此外，龙门山地区地震多发，如2008年汶川8.0级和2013年芦山7.0级强震.虽然国内外研究学者对该区大震震源区的深部孕震结构等进行了大量的研究，但对于地震孕育深部地壳结构特征及相关动力学机制还存在分歧.如：Zhang等(2011)通过地震层析结果，认为龙门山断裂带两侧高低速异常的边界是该区强震的凹凸区;而房立华等(2013)和王夫运等(2015)的地震测深剖面显示断裂带下方中滑脱层是产生研究区地震的主要原因;张培震等(2008)、杨文采等(2015)和王绪本等(2018)的多种研究结果表明研究区强震的发生与地壳内部存在的低速高导层有关
.密度作为地球内部构造最重要的参数之一，能够很好地反应地下物质的运动和变化，高精度的地壳三维密度结构对该区构造形成及演化、地震孕育等深部动力学过程的深入认识具有重要作用.虽然前人已经在该研究区内进行了一定的重力密度探测研究工作，这些探测成果为我们理解和认识研究区地幔变形及强震风险源区的地壳结构和物性特征研究提供了重要的深部资料(姜文亮和张景发，2011;申重阳等，2015)，但这些成果主要集中在二维.已有三维地壳密度数据分辨率和精度较低(方剑和许厚泽，1997;杨文采等;2015;李红蕾等，2017;Lietal.，2017).不足以识别地壳和上地壳深度的细节特征，也不能为解决该区壳幔变形、地震孕育及相关动力学过程争议提供很好的论据(王椿镛等，2016).
3.2研究区已有地壳结构
众所周知，重力数据水平分辨率高，但垂向分辨率低，在实际反演过程中，要想得到有地质意义的结果，还需要添加深部参考约束.相对重力方法来讲，地震成像方法具有较好的垂向分辨率，重力反演中将地震成像结果作为参考约束，可以集两种数据之长，提高重力深部结构的探测能力.本文选择的研究区，在地震学研究方面，已有体波成像(如Zhangetal.，2011;Wangetal.，2017)、面波成像(如Yaoetal.，2008;Lietal.，2011)、接收函数反演(如Baoetal.，2015;Liuetal.，2014)、地震测深成像(如张先康等，2007;王夫运等，2008)等多个结果.然而，由于不同学者使用的研究数据及方法不同，获得的参考地震模型结果在数据分布、分辨率、误差及物性等方面常存在较大差异，如表3所示.——论文作者：李红蕾1，2，陈石1，2*，庄建仓3，张贝1，2，石磊1，2
