任志远李永华强正阳石磊
1） 中国北京 100081 中国地震局地球物理研究所
2） 中国北京 100081 中国地震局震源物理重点实验室
地壳厚度和波速比是描述地壳结构与物质组成的重要参数，其中：地壳厚度不仅可为地震波成像研究中地壳改正、重力补偿与模拟等提供基本参数，还可为地壳形成演化及其动力学过程研究提供重要约束；
地壳波速比则可为地壳物质组成的确定提供约束（Christensen，Fountain，1975）.
确定地壳厚度的方法较多，包括接收函数、面波、重力等方法.接收函数是由远震水平分量与垂直分量之间的反褶积产生的时间序列，由于该方法消除了震源和传播路径的影响，更适于获取地壳和上地幔不连续界面的信息，是获取地壳和地幔结构最常用的工具之一（Burdick，Langston，1977；
Vinnik，1977；
吴庆举，曾融生，1998）.Zhu和 Kanamori （2000）提出了接收函数H-κ叠加技术来估计地壳厚度H和地壳平均波速比κ，该方法根据地壳参数对H值和κ值求取理论到时，并依据该理论到时所对应接收函数的转换波和多次波的振幅大小来确定最优地壳厚度和波速比.该方法不需要人工拾取地震震相到时，且操作简单，被广泛应用于各种构造区域的地壳厚度和波速比κ的计算中（Julià，Mejía，2004；
Wanget al，2010）.然而由于地震台站下方复杂地壳结构的影响，当多次波震相无法被清晰识别时，则很难得到可靠的地壳结构（Zhu，Kanamori，2000；
Liet al，2014）.
重力法也是研究地壳结构和成分的重要工具，该方法具有较高的水平分辨率（Guo，Gao，2018；
Zhaoet al，2018）.完整的布格重力异常包括莫霍界面起伏引起的莫霍重力异常和地壳内密度分布不均匀引起的地壳重力异常（Blakely，1995；
王谦身，2003），该异常不仅与地壳厚度和密度有关，还与地壳波速比κ有关.为此，Lowry和Pérez-Gussinyé （2011）提出了利用接收函数与重力联合反演地壳厚度和波速比的方法，并将其用于美国西部地区地壳厚度和波速比的确定.Shi等（2018）对 Lowry 和 Pérez-Gussinyé （2011）提出的联合估计技术进行了简化和改进，使其易于操作，并具有较高的效率和精度.在不考虑地热影响的情况下，利用接收函数与完全布格重力异常的联合约束，确定了地壳厚度和vP/vS.该方法的优点在于，即使接收函数多次波不清晰，也可以很好地确定地壳厚度和波速比（Christensen，1996；
Lowry，Pérez-Gussinyé，2011）.
上述的接收函数H-κ叠加方法和重力与接收函数联合约束方法，都需要给定一个初始的地壳平均P波速度.当地壳平均P波速度变化时，地壳厚度和波速比也会发生变化（Ammonet al，1990）.为了减少给定 P 波速度所产生的误差，Ma和 Zhou （2007）将面波频散引入接收函数H-κ叠加方法中，通过同时拟合瑞雷波群速度频散和接收函数转换波、多次波震相到时来估计地壳厚度、vP/vS和地壳平均vP.不足之处在于，如果接收函数的多次波不清晰，该方法也无法给出可信的地壳厚度和vP/vS.
本研究试图借鉴上述接收函数与重力、接收函数与面波频散联合方法各自的优势，利用接收函数、重力和面波频散信息联合约束地壳厚度、地壳波速比和平均P波速度，并采用理论和真实观测资料进行测试，以证明该方法的可行性和有效性.
1.1 接收函数H-κ叠加方法
接收函数波形记录中，在直达P波之后是莫霍面的Ps转换波和多次波（PpPs和PsPs＋PpSs），它们相对于直达 P 波的到时可以表达为 （Zhu，Kanamori，2000）：
式中H为地壳厚度，vP和vS分别为P波和S波速度，p为入射波射线参数.
假设vP为常数，则波速比κ＝vP/vS随vS的变化而变化.对于给定的H和κ，相应的转换波和多次波振幅叠加谱为（Zhu，Kanamori，2000）：
式中S（H，κ）为H-κ叠加谱，w1，w2和w3为加权系数，r（ti） （i＝1，2，3）为转换波和多次波的波形振幅.
本研究中设置H和κ的范围分别为20——60 km和1.50——2.00，步长分别为1 km和0.01.根据式（1）——（4）执行扫描和叠加，得到H-κ叠加谱，再得到地壳厚度H和vP/vS的最优估计.
1.2 完全布格重力异常的似然估计
Lowry和Pérez-Gussinyé （2011）认为布格重力异常ΔgB由莫霍面起伏引起的重力异常∆gMolo、地壳密度分布不均匀引起的地壳重力异常Δgcrust和热流引起的重力异常组成.在重力似然反演方法中引入滑动窗口技术，在如此小的区域内，地温的变化通常较小，热流引起的重力影响可以忽略不计（Shiet al，2018；
Guoet al，2019），因此 ΔgB的表达式为
莫霍面重力异常与地壳厚度的关系式、地壳重力异常与地壳波速比的关系式分别如下（Lowry，Pérez-Gussinyé，2011）：
式中：ΔρMoho为莫霍面上下密度差，可以描述地壳密度相对于波速比的变化率；
D为莫霍面深度，D＝H－E，其中H为地壳厚度，E为高程；
D为研究区莫霍面深度的平均值， κ为研究区波速比κ的平均值；
F{·}和F−1{·}分别表示傅里叶变换和反傅里叶变换；
G为万有引力常数，f为波数，c为地壳厚度变化的校正因子.
若已知若干点的地壳厚度、波速比和实测布格重力异常，通过线性回归方法可以得到式（6）和（7）中的 ΔρMoho和 ∂ ρcrust/∂κ；
将其重新代入式（5）——（7），则可正演得到理论重力异常；
再将实测布格重力异常与理论布格重力异常作差求得噪声异常.
其似然函数为
将其取对数，并令 µ和 σ2的一阶导数为0，则
从而可得
这样，采用似然估计法 ［ 式（12）和（13） ］ 分别计算出μ和σ2，再根据式（9）计算似然函数值.以接收函数H-κ叠加相同的范围和步长选取H和κ值，根据式（5）——（9）进行扫描并计算似然函数值，则得到重力H-κ似然谱.
1.3 面波频散似然估计
观测频散和合成频散之间的最小二乘拟合m（H，κ）由理论面波频散与观测频散差的均方根和得到.Ma和Zhou （2007）定义了一个适度函数来得到面波H-κ似然谱：
式中mmax和mmin是H-κ平面内m（H，κ）的最大值和最小值.在计算m（H，κ）时，H值和κ值的选择与接收函数H-κ叠加方法相同，进行网格搜索时固定vP不变，vS随κ值变化.
1.4 联合约束
利用接收函数、重力和面波频散数据联合约束的流程如下：
1） 计算接收函数H-κ叠加谱，将最优地壳厚度和vP/vS估计作为下一步重力估计的初始值.对于采用接收函数H-κ叠加方法不能得到地壳厚度和vP/vS的地震台站，采用频域滤波技术（Guoet al，2013）和密度界面反演技术（Oldenburg，1974）估计莫霍面深度.将高程值加上莫霍面深度，即可得地壳厚度的初始值，然后由接收函数H-κ叠加谱导出对应的vP/vS初始值.
2） 计算重力估计的第i个台站的ΔρMoho和∂ρcrust/∂κ参数.设定滑动窗口的大小，搜索所有地壳厚度和vP/vS以及分布在该观测站中心窗口内的重力异常，然后经网格化后，在一个规则的网格中获取地壳厚度、vP/vS和观测的重力异常数据，将这三个网格数据带入式（5）——（7），利用线性回归算法求解 ΔρMoho和∂ρcrust/∂κ.
3） 将求得的 ΔρMoho和∂ρcrust/∂κ重新代入式（5）——（7），求出模型理论重力异常.
4） 将观测重力异常与模型理论重力异常作差求得噪声异常，采用似然估计法的式（12）和（13）分别计算出μ和σ2，再根据式（9）计算似然函数值.
5） 以与接收函数H-κ叠加相同的范围和步长选取H和κ值，重复步骤3）和4），形成重力反演的H-κ似然谱.
6） 将接收函数H-κ叠加谱与重力反演的H-κ似然谱点乘，得到接收函数和重力联合谱.
7） 采用面波频散似然估计方法（Ma，Zhou，2007），以接收函数H-κ叠加的范围和步长选取H和κ值，通过改变初始κ值改变vS，得到理论面波频散；
求取其与观测频散差的均方根，得到面波拟合的H-κ似然谱.
8） 设定vP范围为 6.0——6.5 km/s，通过改变初始vP（步长为 0.02 km/s），根据步骤 7）得到不同vP所对应的面波H-κ似然谱，将步骤6）得到的联合谱最优估计值与面波H-κ似然谱相比较.当依据前述两个谱分别确定的最优κ值最接近时，所对应的vP为最优vP.
9） 设定最终vP，分别求取接收函数H-κ叠加谱、面波H-κ似然谱和重力H-κ似然谱，依据
求取联合谱，并拾取最终的H值和κ值.
1.5 不确定性估计
本研究将最大振幅标准误差范围内的解定义为联合估计解，参照Eaton等（2006）评估接收函数叠加结果误差的公式，拟采用下式来确定解的不确定性：
式中ηR和ηG分别为接收函数高斯滤波器和重力似然滤波器的宽度，NR和NG分别为用于叠加的接收函数和重力观测的数量.当所有解在H-κ网格空间上趋于高斯分布时，产生最大振幅的解不一定位于高斯分布的中心.因此解的相关不确定性不能用解的总体均值和标准差来描述.为了计算联合约束方法中的不确定度，我们取高斯正态分布函数下15.9%——84.1%范围内即68.2%的范围用于计算平均值，类似于高斯总体分布的一个标准偏差，而最大振幅代表最优解（Delphet al，2019）.
2.1 理论模型及地震、重力数据合成
为检测该方法的可行性，本文首先建立了两种不同的地壳速度模型（图1），并基于这两种模型采用Hrftn96和Surf96正演模拟程序（Herrmann，2013）合成理论接收函数和瑞雷波群速度频散（图2）.
图1 用于正演测试的两个速度模型（a） 简单速度模型（模型Ⅰ）；
（b） 含低速层的速度模型（模型Ⅱ）Fig.1 Two velocity models used for forward testing（a） Simple velocity model；
（b） Velocity model including a low velocity layer
图2 基于图1 中的模型Ⅰ（a）和Ⅱ（b）合成的接收函数（上）和瑞雷波群速度频散 U （下）Fig.2 Synthetic receiver function （upper panels） and Rayleigh wave group velocity dispersion U（lower panels） based on the models Ⅰ （a） and Ⅱ （b） shown in Fig.1
为了对理论重力异常进行正演，我们以待反演台站为中心构建一个半窗口为150 km的地壳厚度和波速比模型（图3a，b），其中数据网格的间距为50 km，地壳厚度为36.2——43.1 km，vP/vS范围为1.59——2.05.本文所用的联合估计方法假设布格重力异常是由莫霍面起伏和壳内密度不均匀共同引起.假设其观测面的高程为0 km，莫霍面深度等于地壳厚度.地壳地幔密度差ΔρMoho为0.5 g/cm3，地壳密度与地壳波速比之比的偏导数∂ρcrust/∂κ为0.25.图4为基于理论模型计算得到的布格重力异常.
图3 用于理论重力异常正演而构建的地壳厚度（a）和波速比（b）分布模型Fig.3 The crustal thickness （a） and vP/vS ratio （b） used for synthetic gravity anomaly
图4 基于模型（图3）正演得到的布格重力异常（a） 莫霍面起伏引起的重力异常；
（b） 地壳内密度不均匀引起的重力异常；
（c） 理论合成的布格重力异常Fig.4 Synthetic Bouguer gravity anomalies for the model in Fig.3（a） The modeled gravity anomalies associated with undulation Moho interface；
（b） The modeled gravity anomalies associated with crustal heterogeneous density distribution；
（c） The modeled Bouguer gravity anomalies with 5% Gaussian noise
2.2 简单地壳模型数据测试
假定由地壳和上地幔组成的简单模型（图1a），中心台站下方的地壳厚度和vP/vS分别为40 km 和 1.75，在莫霍面上方和下方的vP分别为 6.10 km/s和 8.15 km/s，图2a展示了由该模型合成的中心台站的接收函数.在这个简单模型中，接收函数的多次波震相清晰.然而，在实际应用中，多次波的震相往往难以识别，Ps波震相是唯一比较明显的震相.
为了模拟多次波不易识别这一现实情况，我们假设只存在Ps波震相，即在式（4）中设置加权系数w1为1，w2和w3均为0.通过展开接收函数H-κ叠加，得到单极值条带的H-κ叠加谱（图5a）；
对重力异常进行似然估计算法，得到重力H-κ似然谱（图5d）；
将重力H-κ似然谱与接收函数H-κ叠加谱相乘，得到重力与接收函数联合约束H-κ反演谱（图5e）；
根据面波频散似然估计方法得到面波H-κ似然谱（图5b），根据初始P波速度的变化范围6.0——6.5 km/s和步长0.02 km/s对比联合谱与面波频散H-κ似然谱，两图极值的最大值κ的坐标最接近时，取此时的P波速度为H-κ叠加方法的初始P波速度；
将H-κ反演谱与面波频散H-κ似然谱点乘，得到H-κ联合谱（图5g）.由此可以得出，中心台站的地壳厚度、vP/vS比值和初始P波速度的最优估计值分别为（40±1.62） km，（1.75±0.032）和 6.1 km/s，与理论地壳模型相同.
图5 基于模型Ⅰ（图1a）得到的地壳厚度H-波速比κ约束图（图中白线代表最优值68%的置信区间）（a） 接收函数 H-κ 叠加谱；
（b，c） 初始 P 波速度为 6.1 和 6.2 km/s时的面波 H-κ 似然谱；
（d） 重力 H-κ 似然谱；
（e，f） 初始P波速度为6.1和6.2 km/s时，接收函数与重力联合约束谱；
（g） 接收函数、面波和重力联合约束谱Fig.5 H-κ stacking map based on model Ⅰ（Fig.1a） where white lines delineate 68% confidence interval（a） The receiver function H-κ stacking map；
（b，c） The surface wave H-κ likelihood map with initial vP＝6.1 and 6.2 km/s；
（d） The gravity H-κ likelihood map；
（e，f） Normalized SRSG with vP＝6.1 km/s and vP＝6.2 km/s；
（g） Joint H-κ stacking map
2.3 含低速层地壳模型数据测试
假设简单地壳模型中存在一个10 km厚的壳内低速层（图1b），各层的密度分别为2.65，2.7，2.8和 3.3 g/cm3，P 波速度分别为 5.60，5.20，6.80和 8.15 km/s，图2b所示为由该模型合成的中心台站的接收函数（高斯系数为2.5）以及由该模型合成的面波群速度频散.
在式（4）中设置w1＝0.6，w2＝0.3，w3＝0.1，中心台站所产生的接收函数H-κ叠加谱如图6a所示.由于地壳内部存在低速层，接收函数H-κ叠加结果存在两个极值条带，难以估计地壳厚度和vP/vS.本研究利用重力、接收函数和面波频散联合约束生成联合H-κ叠加谱（图6d），据此给出的地壳厚度和vP/vS最优估计值为H＝（40±0.78） km，vP/vS＝1.76±0.028，vP＝6.1 km/s，这与本文构建的速度结构模型相符.
图6 基于模型Ⅱ（图1b）得到的H-κ约束图（图中白线代表最优值68%的置信区间）（a） 接收函数 H-κ 叠加谱；
（b） 面波 H-κ 似然谱；
（c） 重力 H-κ 似然谱；
（d） 初始 P 波速度为 6.1 km/s时，接收函数与重力联合约束谱；
（e） 接收函数、面波和重力联合约束谱Fig.6 H-κ stacking map based on model Ⅱ（Fig.1b） where white lines represent 68% confidence interval（a） Receiver function H-κ stacking map；
（b） Surface wave H-κ likelihood map with vP＝6.1 km/s；
（c） Gravity H-κ likelihood map；
（d） Normalized SRSG with vP＝6.1 km/s；
（e） Joint H-κ stacking map
我国华南地区位于秦岭——大别造山带以南、青藏高原以东，该地区的地质构造演变过程复杂，研究该地区的地壳厚度和波速比可为研究区域地壳变形和动力学机制提供重要依据.这方面已经受到关注并开展了一系列研究（邓阳凡等，2011；
Zhouet al，2012；
Tenget al，2013；
Huanget al，2015；
Shenet al，2016；
Guoet al，2019）.
本文使用杨晓瑜和李永华（2021）的接收函数资料用于实际数据分析.该研究从国家测震台网数据备份中心（郑秀芬等，2009）下载了中国国家地震数字台网在华南地区的336个固定台站在2009年1月至2018年2月期间记录的远震波形数据，并从中挑选震中距在30°——90°之间、M＞5.5且具有清晰P波初至和高信噪比的远震资料用于P波接收函数计算.面波频散数据源自Li等（2013），该研究利用中国国家地震数字台网及GEOSCOPE，K-NET和KZ-NET等计划在东亚大陆周边地区布设的184个台站所记录的区域地震波形数据，采用单台法提取了面波群速度频散，并构建了周期为10——145 s的瑞雷波群速度频散分布图，检测板测试结果显示其多数周期的频散图分辨率为2°.重力数据选取世界重力数据库WGM2012中华南地区的完整布格重力异常数据（Balminoet al，2012），其网格间距为 2″ （图7）.
图7 研究区布格重力异常Fig.7 The complete Bouguer gravity anomalies in study area
为证实本文研究方法的有效性，通过对比前人研究结果（Huanget al，2015；
Guoet al，2019；
Luoet al，2019；
杨晓瑜，李永华，2021），挑选前人接收函数H-κ叠加结果差异较大的两个台站（HB_NZH和HB_YDU）进行深入分析.
3.1 HB_NZH台站结果
图8a为本文计算得到的HB_NZH台站的226个接收函数波形记录，图8b为其观测频散图（Liet al，2013）.
假定地壳厚度H和κ分别在20——60 km和1.5——2.0之间变化，依据前述方法分别计算接收函数H-κ叠加谱（图8c）、面波H-κ似然谱（图8d）和重力H-κ似然谱（图8e）.在接收函数H-κ叠加过程中，当设定w1＝0.6，w2＝0.3，w3＝0.1时，台站接收函数H-κ叠加谱呈单极值条带（图8c）.由此可见，仅依据接收函数得到的地壳厚度和波速比具有较大的误差.
图8f展示了本文利用改进的联合约束方法生成的联合H-κ叠加谱.从该结果来看，该台站下方的地壳厚度和波速比最优估计分别为（36±1.05） km和1.75±0.036，平均P波速度为 6.34 km/s，这与前人（Heet al，2013；
Guoet al，2019）利用接收函数H-κ叠加方法以及接收函数与重力联合反演方法得到的地壳厚度和波速结果（Guo，Gao，2018）较为一致，但是本文利用联合方法给出的误差估计明显要小.
图8 台站HB_NZH数据资料及计算得到的H-κ约束图（图中白线代表最优值68%的置信区间）（a） HB_NZH 台站观测接收函数；
（b） 实测瑞雷波群速度频散 U （Li et al，2013）；
（c） 接收函数 H-κ 叠加谱；
（d） 面波 H-κ 似然谱；
（e） 重力 H-κ 似然谱；
（f） 联合约束谱Fig.8 Data and H-κ stacking map of station HB_NZH where white lines in Figs.（c）–（f）represent the 68 percent confidence interval（a） Observed receiver functions of station HB_NZH；
（b） Observed dispersions U （Li et al，2013）；
（c） Receiver function H-κ stacking map；
（d） Surface wave H-κ likelihood map；
（e） Gravity H-κ likelihood map；
（f） Joint H-κ stacking map
3.2 HB_YDU台站结果
图9a给出了从HB_YDU台站计算的接收函数中选取的133个接收函数，图9b为HB_YDU 台站的观测频散图（Liet al，2013）.
与HB_NZH台站相同，假定地壳厚度H和κ分别在20——60 km和1.5——2.0之间变化，依据前述方法分别计算了接收函数H-κ叠加谱（图9c）、面波频散似然谱（图9d）和重力似然谱（图9e）.在接收函数H-κ叠加过程中，当设定w1＝0.6，w2＝0.3，w3＝0.1时，台站接收函数H-κ叠加谱呈单极值条带（图9c）.由此可见，仅仅依据接收函数得到的地壳厚度和波速比具有较大的误差.
图9f展示了我们利用改进的联合约束方法生成的联合H-κ叠加谱.从该结果来看，该台站下方的地壳厚度和波速比最优估计分别为（40±1.34） km和1.67±0.027，平均P波速度为6.38 km/s，这与Luo等（2019）利用接收函数与PmP走时共同约束给出的结果（H＝40.4 km，κ＝1.67，vP＝6.34 km/s）一致，但与前人利用接收函数、接收函数与重力联合反演方法得到的地壳厚度与波速比结果相差较大.例如：Huang等（2015）利用接收函数H-κ叠加方法得到的结果约为33.5 km和1.73，杨晓瑜和李永华（2021）利用接收函数H-κ叠加方法得到的结果为32.5 km和1.75，而Guo等（2019）利用接收函数与重力联合反演方法得到的结果为32.5 km和1.75.事实上，该台站的接收函数多次波并不清晰，因此，依靠接收函数或接收函数和重力联合地约束很难准确约束地壳厚度和波速比.
图9 台站HB_YDU数据资料及计算得到的H-κ约束图（图中白线代表最优值68%的置信区间）（a） HB_YDU 台站实测接收函数；
（b） 实测瑞雷波群速度频散 U （Li et al，2013）；
（c） 接收函数 H-κ叠加谱；
（d） 面波频散 H-κ 似然谱；
（e） 重力 H-κ 似然谱；
（f） 联合约束谱Fig.9 Data and H-κ stacking map of station HB_YDU where white lines in Figs.（c）−（f）represent 68 percent confidence interval（a） Observed receiver functions；
（b） Observed dispersions U （Li et al，2013）；
（c） Receiver function H-κ stacking map；
（d） Surface wave H-κ likelihood map；
（e） Gravity H-κ likelihood map；
（f） Joint H-κ stacking map
由于接收函数对P波速度敏感，因此，利用接收函数H-κ叠加方法（Zhu，Kanamori，2000）、重力和接收函数联合约束方法（Shiet al，2018；
Lowry，Pérez-Gussinyé，2011）计算地壳厚度和波速度比时，都需要给定一个初始的地壳平均P波速度.
接收函数和面波频散联合约束方法（Ma，Zhou，2007）可以用于同时估计地壳厚度、vP/vS和地壳平均vP，但其缺点是当接收函数多次波不清晰时，估计的地壳结构参数存在较大误差.本文提出了一种利用接收函数、重力数据及面波频散联合约束地壳厚度、波速比和平均P波速度的新方法.通过对合成数据和实际数据的测试，结合以往研究，验证了该方法的可行性，证明了重力和面波约束的介入可以减少H-κ叠加结果的不确定性，提高估计的精度，并得到可靠的平均P波速度.精准的地壳参数可为研究地壳物质的组成、地壳形成演化及其动力学过程提供有利条件.此外，当研究区存在较厚的低速沉积层时，会导致接收函数多次波和转换波的延时效应，该方法计算得到的地壳厚度较真实值大.因此，该方法仍需进一步优化.
感谢审稿专家和编辑提出的宝贵修改意见.
猜你喜欢面波波速台站2013-12-16巴东MS5.1地震前后波速比异常特征大地测量与地球动力学(2022年10期)2022-10-10中国科学院野外台站档案工作回顾中国科学院院刊(2022年8期)2022-09-02土层剪切波速与埋深间的统计关系研究四川工商学院学术新视野(2021年1期)2021-07-22基于实测波速探讨地震反射波法超前预报解译标志铁道建筑技术(2021年4期)2021-07-21地震台站基础信息完善及应用分析科技研究·理论版(2021年20期)2021-04-20gPhone重力仪的面波频段响应实测研究地震研究(2021年1期)2021-04-13一种适用于高铁沿线的多台站快速地震预警方法中小企业管理与科技·下旬刊(2021年1期)2021-03-02自适应相减和Curvelet变换组合压制面波石油地球物理勘探(2020年5期)2020-10-17灰岩声波波速和力学参数之间的关系研究湖南理工学院学报（自然科学版）(2020年1期)2020-04-17铁路无线电干扰监测和台站数据管理系统应用研究铁道通信信号(2020年7期)2020-02-06