地球物理学报  2016, Vol. 59 Issue (3): 959-971   PDF    
基于邻域算法的瑞利面波垂直-水平振幅比及频散曲线联合反演及应用
袁艺1,2, 姚华建1,2, 秦岩1,2    
1. 中国科学技术大学地球和空间科学学院, 地震与地球内部物理实验室, 合肥 230026;
2. 蒙城地球物理国家野外科学观测研究站, 安徽 亳州 233527
摘要: 瑞利面波垂直-水平振幅比(或ZH振幅比)是一个随频率变化的函数,对于台站下方浅层地壳结构非常敏感,且具有和频散资料不同的深度敏感核,是传统频散反演方法的一个很好的补充,从而可以将基阶瑞利面波的ZH振幅比和面波频散数据联合起来更好地反演获得观测台站下方的速度结构.本文提出了基于邻域算法的面波频散曲线与ZH振幅比联合反演方法,我们进行了基于理论模型的模拟测试,证明了联合反演是一种更为可靠的反演方法,且能更好地约束浅层地壳结构.相比于频散曲线单独反演,联合反演不仅可以精确反演获得地壳的Vs结构,对分层地壳的Vp/Vs也能很好地约束.然后我们将联合反演算法应用于实际测量数据,获得了中国西南昆明台(KMI)下方更为准确的地壳横波速度结构及Vp/Vs模型.
关键词: 联合反演     邻域算法     瑞利面波     相速度频散     ZH振幅比     地壳结构    
Joint inversion of Rayleigh wave vertical-horizontal amplitude ratios and dispersion based on the Neighborhood Algorithm and its application
YUAN Yi1,2, YAO Hua-Jian1,2, QIN Yan1,2    
1. Laboratory of Seismology and Physics of Earth's Interior, School of Earth and Space Sciences, University of Science and Technology of China, Hefei 230026, China;
2. National Geophysical Observatory at Mengcheng, Anhui Bozhou 233527, China
Abstract: Rayleigh wave ellipticity(or ZH ratio) is a function of frequency and is particularly sensitive to shallow crustal structure beneath the seismograph station. Since depth sensitivity kernels of ZH ratios are different from those of dispersion data, the ZH ratio provides good complementary information for the dispersion-based inversion method.Therefore, we can combine the ZH ratio and dispersion data of Rayleigh wave fundamental mode to better invert for the velocity structure under a specific seismograph station. In this paper, we propose a joint inversion method using the dispersion and ZH ratio data based on the Neighborhood Algorithm. We conduct synthetic tests based on a theoretical model and prove the robustness of the joint inversion method, which can better constrain the shallow crustal structure. Compared to traditional inversion methods that only use dispersion data, the joint inversion can provide a more accurate crustal Vs model as well as Vp/Vs ratios for the layered crust. Finally, we apply the joint inversion technique to real measurements and obtain a more accurate crust shear velocity and Vp/Vs model beneath the station at Kunming(KMI) in southwest China.
Key words: Joint inversion     Neighborhood Algorithm     Rayleigh waves     Phase velocity dispersion     ZH ratio     Crustal structure    
1 引言

利用地震面波和背景噪声数据测量得到的面波频散数据是反演地下结构一种十分重要且日趋成熟的方法(Shapiro et al.,2005Yao et al.,2008),但是这种方法存在一定的缺陷:因为某一周期的面波相速度或群速度对地下一定深度范围的结构比较敏感,周期越长,敏感深度越大,但敏感深度范围越宽(徐果明等,2007),从而导致深部的反演结果的分辨率降低;且由于高质量的短周期频散数据较难获得,所以通常很难准确获得地壳浅部几公里的速度结构.如果能提高对地表几百米到几千米的速度结构的分辨率,我们将能更好地关联地表的地质构造特征和地壳深部的结构异常及演化过程,且对认识和预测地震地质灾害及矿产油气资源勘查都具有重要的意义.

基阶瑞利面波在传播的过程中径向(H)和垂向(Z)的振动分量具有90°的相位差,其垂向和径向分量的振幅比值与地下浅部结构有密切的关系,我们把这一参数称为ZH振幅比.ZH振幅比是一个随频率变化的函数,对于横向各向同性,纵向随深度变化的分层地球模型来说,这一函数只与基阶瑞利面波所到达位置的地下介质力学参数有关,与其他参数(如震源)无关(Tanimoto and Alvizuri,2006).随着宽频带地震台的增加及数据测量、处理技术的提升,关于ZH比的研究成果也越来越多.Tanimoto和Alvizuri(2006)利用基阶瑞利面波径向和垂向分量之间的90°相位差从连续的噪声信号中搜寻基阶瑞利面波,计算其ZH振幅比并反演地下浅层横波速度结构.Tanimoto和Rivera(2008)利用数值方法得出了长周期(20~250 s)面波ZH振幅比的深度敏感核曲线,相比于基阶瑞利面波相速度频散曲线,ZH比对浅层结构更加敏感.Yano等(2009)利用中长周期瑞利面波的ZH振幅比反演了GEOSCOPE台网不同台站下方的地壳和上地幔结构.Tanimoto等(2013)发展了从连续噪声信号中计算ZH比的新方法.

面波频散曲线和ZH振幅比都可以用来反演地下速度结构,但是同样的频率范围两者的敏感深度范围不同,ZH振幅比正好可以弥补频散曲线对地壳浅层结构约束不足的缺陷.Boore和Toksöz(1969)很早就做了频散曲线和ZH振幅比联合反演的试验.Lin等(2012)将中长周期(24~100 s)的ZH振幅比和8~100 s的频散曲线结合起来用于反演美国西部地区地壳上地幔的结构.Lin等(2014)利用三分量背景噪声的互相关数据计算了台站下方中短周期的ZH振幅比(8~24 s),和之前从地震数据中得到的24~100 s ZH振幅比及8~100 s相速度频散结果结合,可以更好地用于联合反演美国西部地区地壳及上地幔的三维结构.上述这些ZH比和频散曲线的联合反演均是基于线性最小二乘的反演方法,需要计算ZH比和频散数据的敏感核.Chong等(2015)提出了基于模拟退火算法的面波频散曲线及ZH振幅比联合反演,直接在模型空间进行随机搜索从而获得全局最优模型,所以不需要计算敏感核.

本文提出一种基于邻域算法(Neighborhood Algorithm,或NA)(Sambridge,1999a1999b)的面波频散曲线和ZH振幅比的联合反演,并通过贝叶斯分析方法获得模型反演的误差、分辨率及参数之间的相关性.我们首先通过模拟测试数据证明基于邻域算法的联合反演的可靠性,接下来将联合反演算法应用于实际地震数据.

2 瑞利面波ZH振幅比及频散的深度敏感核

我们定义ZH振幅比为基阶瑞利面波在地表处(z=0时)垂直向与水平向振幅的比值,即

其中,ω代表频率,|r1(ω|z)|和|r2(ω|z)|为对应于频率ω在不同深度z处的瑞利面波本征函数值(Aki and Richards,1980).

面波频散资料是反演地壳和上地幔三维横波速度结构的重要工具,对于基于线性的反演方法,关键在于计算不同频率的相速度或群速度对横波波速、纵波波速以及密度的深度敏感核.对于ZH振幅比,Boore和Toksöz(1969)的研究认为影响ZH振幅比函数的介质力学参数同样包括纵波波速、横波波速和介质密度,Q值影响可以忽略,其中ZH振幅比函数对横波波速的变化最敏感.某一频率(f)瑞利面波振幅比(或相速度)的深度敏感核公式可以表示为

η为ZH振幅比(或相速度c);ρ为密度,α为纵波波速,β为横波波速,δp(p=ρ,α,β)分别为ZH振幅比(或相速度)、密度、纵波波速、横波波速的扰动,这些参数均为深度z的函数;Kp被称为敏感核,是频率和深度的函数.

图 1给出了一维平层速度结构下(地壳厚度为45 km)3个不同周期瑞利面波ZH振幅比及相速度对结构参数的深度敏感核,可以清楚地观察到ZH振幅比和相速度具有非常不同的深度敏感核,尤其是对横波速度.与相速度相比,ZH振幅比的敏感深度更浅,约为频散曲线的一半(Tanimoto and Rivera,2008),而且ZH振幅比在地表处的敏感性最强,因此在相同频带下能反演比频散曲线更浅深度的横波速度结构.瑞利波相速度频散对地壳的纵波速度及密度有一定的敏感性,ZH比对地壳(浅部)密度也有一定的敏感性,但对纵波速度的敏感性很低.

图 1 不同周期基阶瑞利面波相速度(a)和ZH振幅比(b)对纵横波速和密度的深度敏感核 Fig. 1 Depth sensitivity kernels at different periods for Rayleigh wave fundamental mode phase velocity (a) and ZH ratio (b) to Vs, Vp and density, respectively
3 ZH振幅比及频散的联合反演方法

因为ZH振幅比具有和面波频散不一样的深度敏感核,且它对浅部结构的约束性更强,所以我们可以将其与面波频散曲线结合进行联合反演,以更好地约束地壳的速度结构(Lin et al.,2012).本文中,我们发展了一套基于邻域算法(NA)(Sambridge,1999a1999b)的ZH振幅比及频散曲线联合反演算法.

3.1 邻域算法(NA)

邻域算法是一种基于全局搜索的反演方法,它包括两个阶段:第一阶段应用“Voronoi cells”(维诺单元)的特殊几何分割方法在模型空间中搜索可接受的(“好”的)模型(Sambridge,1999a),第二阶段对第一阶段生成的模型进行贝叶斯统计分析.得到的模型参数的一维边际后验概率分布和高斯函数分布的偏离程度可以作为反演结果可靠与否的判断标准,且以一维概率分布的宽度(即标准差)来表征模型的不确定性.二维边际后验概率密度函数可以用于衡量两个模型参数之间的相互折衷性(trade-off),其分布越集中越接近圆形,说明二者的相互折衷性越小.一维和二维边际概率分布函数都是估计分辨率和模型参数不确定性的有力工具(Sambridge,1999b).

3.2 理论模型反演

为了验证频散曲线与ZH振幅比联合反演相对于频散曲线单独反演地壳结构的优越性,我们先用理论模型正演产生的数据来进行NA反演测试.

3.2.1 理论模型构建及数据准备

我们构建的理论模型包括地壳和上地幔两部分,上地幔310 km深度之下设为半无限空间.地壳厚度设为45 km,分为四层,第一层的纵横波速比(Vp/Vs)为1.9,其余三层均为1.7.我们先设定了地壳中四层的横波速度Vs,根据Vp/Vs值计算相应层的纵波波速Vp,再根据经验公式(Brocher,2005)可得密度值;上地幔部分从45 km到310 km,来源于ak135模型(Kennett et al.,1995),如表 1所示.根据此模型我们采用Herrmann和Ammon(2004)的方法正演得到瑞利面波相速度频散数据31个(周期范围为5~150 s)以及ZH振幅比数据26个(周期范围为5~120 s),这些数据(未加噪声)将用于下面的NA测试反演之中.

表 1 邻域算法反演的10个参数的参考值及其变化范围 Table 1 The reference values and their perturbation ranges for 10 parameters in the NA inversion
3.2.2 NA联合反演

表 1中给出我们采用NA算法拟反演的10个参数,包括4层地壳和4层上地幔的横波速度,以及地壳第一层(0~6 km)和其余三层(6~45 km)的纵横波速比(Vp/Vs).对于横波速度结构,我们以理论模型为参考值,对每个参数给定相应的变化范围(如表 1所示).这里上地幔的各层包含ak135模型中的若干厚度较薄的层,在反演中模型参数值将在参考模型值的基础上整体变化.由于面波数据对间断面约束较差,且和相邻层的速度结构有很大的折衷性(例如Yao et al.,2008),所以我们通常不反演层厚(或间断面的位置).

NA算法中将产生各层的横波速度,我们将根据纵横波速比(如果也作为模型参数的话)或经验公式(Brocher,2005)计算出地壳部分相应层的纵波速度,之后再根据经验公式(Brocher,2005)从纵波速度计算出密度值.但是该经验公式不适用于上地幔部分,所以我们采用如下另一经验关系式来计算纵波速度和密度相对于横波速度的变化(Masters et al.,2000Yao et al.,2008):

反演中误差函数采用了L1范数的形式:
其中N是包括ZH振幅比和相速度频散数据的所有待研究周期范围内的数据点总数,这里N=57;dipreddiobs是第i个周期处的模型计算值及观测值;σi是第i个周期处观测数据的标准差,这里σi取0.01diobs.

考虑到反演结果的可靠性以及收敛速度,经过反复试验,我们在每一次迭代中将产生20个新模型(Ns=20),结合先前生成的模型,每次迭代时模型空间将被分割生成新的维诺单元图,接着从所有的维诺单元中重新选出10个最好的单元(Nr=10)再生成20个新模型,如此反复迭代,直至搜索过程收敛.之后我们对生成的模型进行贝叶斯统计分析获得每个模型参数的后验平均值(posterior mean)及其标准误差(对应于68%的置信区间).

图 2-4显示了我们直接反演10个参数(见表 1)的结果.图 2黑色细实线给出了后验平均模型,而虚线给出了搜索得到的最优模型,最优模型和理论模型符合得非常好,但后验平均模型和理论模型在部分参数上存在一定的偏差.由于Moho面附近的第4层(下地壳)和第5层(上地幔顶部)的横波速度之间存在相互折衷(如图 4的二维后验概率密度函数所示),导致下地壳后验平均横波速度较真实模型偏高约0.12 km·s-1,而上地幔顶部后验横波速度则偏低约0.18 km·s-1.

图 2 一步联合反演得到地壳及上地幔的横波速度结构(左),由统计后验平均模型和最优模型正演得到的相速度频散曲线(中)和ZH振幅比(右)与理论计算值的比较
左图中粗实线、细实线、点化线分别代笔理论模型、后验平均模型和最优模型,灰色区域给出后验平均模型的标准偏差; 中间及右图中三角形代表由后验平均模型正演出的频散数据及ZH振幅比,星号代表最优值,菱形代表理论模型计算值.
Fig. 2 Shear-wave speed model of crust and upper mantle obtained from the one-step joint inversion (Left), comparison of phase velocity dispersion (Middle) and ZH ratios (Right) calculated from the statistical posterior mean, optimum and theoretical models, respectively
In the left plot, the thick, thin and dashed line represents the theoretical, posterior mean, and optimum model, respectively. The gray region gives the standard error of the posterior mean model. In the middle and right plots, the triangles show dispersion and ZH ratio data computed from the posterior mean model, the stars from the optimum model and the diamonds from the theoretical model.

图 3 频散及ZH比数据同时反演地壳及上地幔10个参数的一维后验概率密度函数分布(阴影区域所示)
横轴为横波速度或纵横波速比的变化范围,纵轴代表归一化的后验概率密度,每幅图中的黑色竖直实线分别 对应于后验平均模型.对于横波速度参数(上两行),横坐标0值处对应于真实模型的值.
Fig. 3 1-D marginal posterior probability density functions (PPDFs) shown as the shaded area of the 10 model parameters from the joint inversion of the dispersion and ZH ratio data
The horizontal axis shows the perturbation range of shear velocity or Vp/Vs. The vertical axis gives the normalized posterior probability density. In each plot, the black solid vertical line shows the parameter value of the posterior mean model. For the shear velocity parameters (top two rows), the zero value of the horizontal axis corresponds to the true model value.

图 4 频散及ZH比数据同时反演地壳及上地幔10个参数的二维后验概率密度函数分布
黑线、蓝线和红线分别代表99%、90%和60%置信级的等值线,白色三角形代表后验平均模型.
Fig. 4 2-D marginal posterior probability density functions (PPDFs) of the 10 model parameters from the joint inversion of the dispersion and ZH ratio data
The black, blue and red lines are the contour lines for the 99, 90 and 60 percent confidence level, respectively. The posterior mean model is shown as the white triangle in each plot.

图 3可见,虽然后验概率密度函数呈高斯函数分布,且大部分模型参数的后验平均值与真实模型值也比较相符,但模型后验函数宽度较大,说明反演模型的标准误差偏大.由于NA算法的全局搜索参数较多时,各参数之间的相互折衷性较大,可能会导致反演结果的不稳定,从而使得模型搜索落入到区域极小值附近,使得最终得到的后验平均模型与真实模型存在一定的偏差.

因此,我们尝试采用两步反演的模式,第一步仅反演地壳上地幔若干层的横波速度结构,这里为8个参数,分别为4层地壳和4层上地幔的横波速度(见表 1的前8个参数),以求将各层的横波速度初步确定.获得较为可靠的上地幔模型后,我们进行第二步反演,共反演7个参数,分别为4层地壳(同第一步分层)和1层上地幔(45~310 km)的横波速度Vs,以及地壳第一层(0~6 km)和其余三层(6~45 km)的纵横波速比(Vp/Vs).第二步的参考模型选取第一步反演获得的地壳上地幔横波速度结构的后验平均模型值,并将上地幔的反演参数由4个变成1个,即Moho-310 km深度的上地幔横波速度在第二次反演时仅相对于第一次反演获得的该深度范围的横波模型做整体的扰动变化.通过反演参数的变化,我们希望更好地约束地壳横波速度结构及第一步反演中未涉及的Vp/Vs值.

两步NA联合反演结果如图 5-7所示.图 5给出NA反演所得后验平均模型值与理论模型的对比图,两模型非常接近,且相比于一步反演的结果,最终反演模型的误差很小.另外,面波频散和ZH比的后验平均模型计算值与理论值几乎重合(图 5),这个结果证明了反演所得后验平均模型的正确性,验证了NA联合反演算法的可靠性和精确度,结果也明显优于单步反演法(图 2-4).

图 5 两步联合反演得到地壳及上地幔的横波速度后验平均模型(实线)及理论模型(虚线)(左),由后验平均模型计算得到的相速度频散曲线(中)和ZH振幅比(右)与理论模型计算值的比较 Fig. 5 Shear-wave speed posterior mean model (solid line) and theoretical model (dashed line) of the crust and upper mantle obtained from the two-step joint inversion (Left), comparison of phase velocity dispersion (Middle) and ZH ratios (Right) calculated from the posterior mean and theoretical models, respectively

图 6 频散与ZH振幅比两步法联合反演(黑色实线及阴影部分)和频散数据单独反演(虚线)的一维后验概率密度函数分布对比. 每幅图中的竖直实线或虚线分别对应于联合反演或频散曲线单独反演的后验平均模型 Fig. 6 1-D marginal posterior probability density functions (PPDFs) of the seven model parameters from the dispersion and ZH ratio two-step joint inversion (black solid lines and the shaded area) in comparison with the results from dispersion data only (dashed lines). In each plot, the solid or dashed vertical line shows the parameter value of the posterior mean model for the joint inversion or inversion with dispersion data only, respectively

图 7 基于理论模型合成数据的NA两步联合反演所获得的不同模型参数之间的二维后验概率密度函数分布 Fig. 7 2-D marginal posterior probability density functions (PPDFs) between different model parameters from the NA joint inversion based on the synthetic data from the theoretical model

图 6的阴影部分给出联合反演所得模型的一维边际后验概率密度函数分布,各模型参数的后验平均值如图中竖直黑线所示,最大似然模型可以从这些一维分布图的峰值得到(Sambridge,1999b).相比于一步反演的结果(见图 3),图 6显示各反演参数的一维边际后验概率密度函数分布基本都呈高斯分布,且范围较窄,显示出模型的标准误差小,分辨能力强.除了横波速度Vs,纵横波速比Vp/Vs也与理论设定值相符.6 km-Moho深度的Vp/Vs值的标准误差较0~6 km深度Vp/Vs值的标准误差明显要大,显示出联合反演对地壳浅部的Vp/Vs值有更强的约束.

图 7给出了不同模型参数之间的二维边际后验概率密度函数分布.与图 4中一步反演结果对比,横波速度参数之间的二维边际后验概率密度函数等值线弥散性明显变小,几乎呈圆形分布,且半径很小,显示出地壳中4层的横波速度均能很好地确定,且各层速度之间的相关性也很小.0~6 km这一层的横波速度与Vp/Vs值几乎没有相关性,但0~6 km的横波速度与6 km-Moho深度的Vp/Vs值存在较为明显的负相关性.0~6 km的Vp/Vs值与6 km-Moho深度的Vp/Vs值没有明显的相关性,但0~6 km的Vp/Vs值的反演误差要明显小一些,这与图 6中一维边际后验概率密度函数的分析结果一致.

经过对比,两步反演的结果更好,分辨率更高,所以之后我们采用两步反演法对实际数据进行反演.

3.2.3 瑞利面波频散数据单独反演

使用频散数据来反演地下速度结构是比较常见的方法.这里我们执行和3.2.2中联合反演相同的反演流程,但仅使用31个面波频散数据参与反演.本次反演中得到的后验平均模型和理论值差距较大(图 6垂直虚线所示),反演得到的一维后验概率分布函数分布较宽,尤其是0~6 km和6~15 km这两层的横波速度;且有的参数出现明显的多峰值情况,说明仅用频散曲线反演得到的模型参数的非唯一性较强.经比较,我们舍弃了反演得到的后验平均模型而采用了所有模型中的“最优”模型(数据拟合误差最小),结果如图 8所示.由图可见,单独反演的最优模型虽然能够较完美地拟合频散数据,但反演得出的ZH振幅比和理论值的差距却较大.这说明传统的频散数据反演方法存在一定的局限性,对一些模型参数的约束不强,从而会导致部分参数的反演结果存在较大的误差.

图 8 频散曲线单独反演结果,左图中黑色实线为NA反演中的最优模型,其它同图6 Fig. 8 Inversion results from the dispersion data only. In the left plot, the black solid line gives the best fitting model from the NA inversion. Other symbols are similar as in Fig.6
4 应用4.1 数据

通常无论是地震面波还是背景噪声,ZH振幅比的实际测量结果都存在较大的误差.Tanimoto和Rivera(2008)认为这种误差分布是稳定的,可以用统计方式进行处理.为了减小误差,我们在提取地震面波振幅比的数据处理过程中设定了多个筛选条件.如信噪比(SNR)须大于10,瑞利面波径向和垂向分量之间要近似满足90°相位差,并且删除ZH振幅比明显异常的数据点等.图 9为我们得到的昆明台下方45 s和85 s周期的ZH振幅比分布,可以看到它们大致呈现高斯分布.最终处理所得ZH振幅比曲线与已有的地壳上地幔模型(Yao et al.,2010)产生的ZH振幅比曲线基本吻合,证明本次反演所用ZH振幅比测量数据较为可靠.

图 9 用不同地震面波信号测量得到的不同周期(45 s,85 s)昆明台ZH振幅比分布的直方图(近高斯分布),图中黑色竖直线代表对应于各自周期的ZH振幅比均值 Fig. 9 The histograms of ZH ratios (close to Gaussian distribution) at different periods (45 s and 85 s) of the KMI station measured from earthquake surface waves. The black vertical line in each figure represents the mean value of ZH ratios at the corresponding period, respectively

本次联合反演频散观测数据共39个,来源于Yao等(2010)对青藏高原东南部地区瑞利面波成像在昆明台所处经纬度处的插值结果,周期范围10~150 s.ZH振幅比的观测值共29个,周期范围25~110 s.

4.2 参考模型

NA联合反演参考模型包括地壳及上地幔两部分,上地幔310 km以下设为半无限空间.地壳分为三层(0~15 km,15~30 km,30~45.554 km),从上到下各层横波速度分别为3.3 km·s-1、3.6 km·s-1、3.8 km·s-1,Moho面深度参考Yao等(2010),设为45.554 km,上地幔部分从45.554 km到310 km仍然取自ak135模型.

4.3 NA联合反演步骤及结果

我们采用和理论测试中相同的反演流程,NsNr取值也与理论测试相同.

第一步:反演7个参数,分别为3层地壳和4层上地幔(Moho-80 km,80~120 km,120~210 km,210~310 km)的横波速度.

第二步:反演6个参数,分别为3层地壳和1层上地幔(Moho-310 km)的横波速度Vs,以及地壳第一层(0~15 km)和其余两层(15 km-Moho)的纵横波速比Vp/Vs.第二步的参考模型选取第一步反演获得的后验平均模型值,在第一步的基础上对地壳各层的横波速度进行调整,并将重点放在纵横波速比的反演上.

联合反演结果如图 10-12所示.图 10中,黑色粗实线为最终反演模型,由该模型正演得到的相速度频散曲线及ZH振幅比曲线基本都处于观测值的误差范围内.但当周期在140 s以上时,计算值与观测值差距较大,可能是因为长周期的观测频散数据存在较大的误差,所以上地幔深部的反演结果不太可信.图 11显示各反演参数的一维边际后验概率密度函数分布基本呈高斯型分布,且范围较窄,显示出很好的模型分辨能力.需要指出的是,15 km-Moho范围内反演得到的Vp/Vs值的后验概率分布有两个峰值,其中一个明显高于另一个,此处我们并不像选取其他参数那样选取后验概率平均模型值(位于两峰之间),而是取高的峰(最大似然值)对应的Vp/Vs值作为最终反演结果.另外,我们将地壳分为两层并分别反演它们的Vp/Vs值,结果都在1.72左右,这和已有的利用接收函数方法得到的该地区的泊松比大致符合(黄建平等,2008徐强等,2009),显示出ZH振幅比和频散曲线的联合反演可以用于获得分层地壳的Vp/Vs值(或泊松比).图 12为各反演参数值之间的二维边际后验概率密度函数等值线,和理论测试结果(图 7)对比,实际数据结果的弥散性更强,但最内层的60%等值线(红色)几乎呈圆形且半径很小.

图 10 实际联合反演得到的昆明台下方地壳及上地幔的横波速度结构(左),实际数据联合反演得到的频散数据(中)和ZH振幅比(右)与观测值的对比.左图中黑色实线为最终反演得到的后验平均模型,灰色区域显示其标准误差,虚线和点虚线分别为第一次和第二次NA反演的参考模型;中间和右图中带误差棒的星号代表观测数据,虚线和三角分别对应于第一步和第二步NA反演的参考模型计算得到的理论数据,实线对应最终反演获得的后验平均模型计算得到的理论数据 Fig. 10 Shear-wave speed model of crust and upper mantle beneath the station KMI obtained from real seismic data through joint inversion (Left), comparison of the phase velocity dispersion obtained from inversion and the observed dispersion data (Middle), same as the Middle plot but for ZH ratios (Right). In the left plot, the black solid line gives the final posterior mean model after inversion with the gray area showing its standard error, and the dashed line and the dot dashed line is the reference model for the first stage and second stage NA inversion, respectively. In the middle and right plots, the stars with the error bar are the observed data; the dotted line and the triangles are the synthetic data from the reference model in the first and second stage NA inversion, respectively; the solid line shows the synthetic data computed from the final posterior mean model after inversion

图 11 实际地震数据联合反演得到的模型参数的一维后验概率密度函数分布(与图6类似).除右下角的图外,黑色竖直实线对应于后验平均模型;右下角的图(Vp/Vs [15-Moho])的黑色竖直实线对应于模型的最大似然值 Fig. 11 1-D marginal posterior probability density functions (PPDFs) of the model parameters from the joint inversion of real data (similar as Fig.6). Except the bottom right plot, the vertical line in each plot represents the posterior mean model; for the bottom right plot, the vertical line shows for the model with the most likelihood

图 12 实际数据联合反演得到的二维后验概率密度函数分布(与图7类似) Fig. 12 2-D marginal posterior probability density functions (PPDFs) of the model parameters from the joint inversion of real data (Similar as Fig.7)

值得注意的是,虽然第一步反演的参考模型的选取比较随意,但如果参数范围合适,最终结果并未受影响:反演第一步所得模型(即为第二步反演的参考模型)的正演结果(频散曲线和ZH比)已经和观测值比较接近了,到反演第二步,其反演模型的正演结果已经基本落在了观测值的误差范围以内.由图 10可见,最初的参考模型对应的正演曲线(虚线)和最终反演模型对应的曲线(粗实线)相差较大,但经过两步联合反演,即使是不太准确的初始模型也可以得到比较准确的结果,有力验证了我们所使用的联合反演算法的稳定性.

5 讨论及结论

用相速度频散曲线来反演地下速度结构的方法应用广泛,但在某些深度范围的分辨率不足,尤其对地壳浅部的约束不够.ZH振幅比有独立于相速度的深度敏感核,对于同一周期,其敏感深度范围与相速度的敏感深度范围显著不同,更加接近地壳浅层.二者结合起来,可显著提高对地壳浅部结构的分辨能力.

本文提出了基于邻域算法(NA)的频散曲线和ZH比联合反演方法,相比于一步反演,采用两步反演更稳定,能够更为准确地获得地壳上地幔横波速度结构及分层地壳的Vp/Vs值.通过理论测试我们发现,当一步反演的参数较多时,即使得到的最优模型和理论模型非常接近,统计所得的后验平均模型却可能和理论模型相差较远,模型搜索容易陷入局部极小.另外,联合反演明显优于频散曲线单独反演,能更好地约束模型参数并降低其反演误差.之后我们将两步联合反演应用于实际测量数据,获得了昆明台下方更为可靠的地壳速度结构模型.此外,基于邻域算法的联合反演可以通过后验概率密度分布函数更好地分析模型的误差和各参数间的相关性.

之前的相关研究或者没有给出Vp/Vs的反演值,或者发现其与密度值之间很强的相互折衷性而没能得到确定的结果.在本文中,我们采用了更为合理的纵波速度与密度之间的经验公式(Brocher,2005),并不反演密度,从而更好地约束Vp/Vs值.从昆明台真实数据的反演结果看,联合反演获得的分层地壳的Vp/Vs值与以前通过接收函数获得的地壳的平均泊松比结果比较吻合,说明频散曲线和ZH比的联合反演可以独立地提供分层地壳泊松比的约束,这为认识地壳不同深度介质的物理和成分属性提供新的约束.

另外有一点值得注意:在用实际地震数据进行反演时,我们将地壳分为了三层,其厚度基本相等.而理论测试时,我们却分了四层,相当于将第一层(15 km)又分为两层(0~6 km和6~15 km).这是因为ZH振幅比主要是对浅层的约束比较强,且起始周期越小反演效果越好.理论测试数据的周期从5 s开始,所以对所划分的四层地壳的约束都很好.而实际反演中,相速度的周期从10 s开始,ZH振幅比的周期则从25 s开始,所以很难约束浅部几千米的结构.对于地震面波,我们通常测量得到的ZH振幅比的周期都在20 s以上(Lin et al.,2012),为了更好地约束地壳浅层的速度结构,我们可以通过连续背景噪声信号(Tanimoto et al.,2013)或背景噪声的互相关数据(Lin et al.,2014)来得到短周期的ZH振幅比.

综上,频散曲线和ZH比的联合反演确实要优于频散曲线的单独反演,但联合反演结果的分辨率也受到数据质量、周期范围等因素的影响.此外,考虑到面波频散及ZH振幅比对界面信息的约束不足,本文并没有反演Moho面或壳内其他界面的深度,若在现有基础上再加入对分界面敏感的接收函数数据,可能会得到更好的效果.面波频散和接收函数均对介质的S波速度结构敏感,但二者的侧重点不同.面波频散曲线可以揭示介质的横波速度结构变化,而接收函数对横波的绝对速度约束比较差,但在确定转换界面深度方面具有独特的优势,故联合面波频散和接收函数资料可以获得比单一方法更为可靠的横波速度结构(胡家富等,2005刘启元等,2010).近年来,面波频散和接收函数的联合反演广泛应用于地壳上地幔速度结构研究中(胡家富等,2005刘启元等,2010).开展联合反演的目的是为了实现优势互补,将面波频散曲线、ZH振幅比曲线及接收函数联合来反演地壳上地幔速度结构,可以更为准确地约束模型信息,减小反演结果的不确定性,这也是目前研究的热点方向之一.

致谢  感谢两位审稿人对本论文提出的修改建议,感谢美国University of Colorado,Boulder的沈伟森博士提供自动的面波时频分析程序用于地震面波振幅比的测量工作.
参考文献
[1] Aki K, Richards P G. 1980. Quantitative Seismology, Vol. 1:Theory and Methods. San Francisco:WH Freeman and Company.
[2] Boore D M, Toksöz M N. 1969. Rayleigh wave particle motion and crustal structure. Bulletin of the Seismological Society of America, 59(1):331-346.
[3] Brocher T M. 2005. Empirical relations between elastic wavespeeds and density in the earth's crust. Bulletin of the Seismological Society of America, 95(6):2081-2092.
[4] Chong J J, Ni S D, Zhao L. 2015. Joint inversion of crustal structure with the rayleigh wave phase velocity dispersion and theZHRatio. Pure and Applied Geophysics, 172(10):2585-2600.
[5] Herrmann R B, Ammon C J. 2004. Computers programs in seismology, surface waves, receiver functions and crustal structure, version 3.30. Department of Earth and Atmospheric Sciences. Saint Louis University.
[6] Hu J F, Zhu X G, Xia J Y, et al. 2005. Using surface wave and receiver function to jointly inverse the crust-mantle velocity structure in the West Yunnan area. Chinese J. Geophys.(in Chinese), 48(5):1069-1076.
[7] Huang J P, Chong J J, Ni S D, et al. 2008. Inversing the crustal thickness under the stations of China via H-Kappa method. Journal of University fo Science and Technology of China (in Chinese), 38(1):33-40.
[8] Kennett B L N, Engdahl E R, Buland R. 1995. Constraints on seismic velocities in the earth from traveltimes. Geophysical Journal International, 122(1):108-124.
[9] Levshin A L, Pisarenko V F, PogrebinskyG A. 1972. On a frequency-time analysis of oscillations. Annales Geophysicae, 28(2):211-218.
[10] Lin F C, Schmandt B, Tsai V C. 2012. Joint inversion of rayleigh wave phase velocity and ellipticity using usarray:Constraining velocity and density structure in the upper crust. Geophysical Research Letters, 39:L12303.
[11] Lin F C, Tsai V C, Schmandt B. 2014. 3-D crustal structure of the western United States:Application of Rayleigh-wave ellipticity extracted from noise cross-correlations. Geophysical Journal International, 198(2):656-670.
[12] Liu Q Y, Li Y, Chen J H, et al. 2010. Joint inversion of receiver function and ambient noise based on Bayesian theory. Chinese J. Geophys.(in Chinese), 53(11):2603-2612.
[13] Masters G, Laske G, Bolton H, et al. 2000. The relative behavior of shear velocity, bulk sound speed, and compressional velocity in the mantle:implications for chemical and thermal structure.//Earth's Deep Interior:Mineral Physics and Tomography From the Atomic to the Global Scale.Washington DC, American Geophysical Union, 63-87.
[14] Sambridge M. 1999a. Geophysical inversion with a neighbourhood algorithm-Ⅰ. Searching a parameter space. Geophysical Journal International, 138(2):479-494.
[15] Sambridge M. 1999b. Geophysical inversion with a neighbourhood algorithm-Ⅱ. Appraising the ensemble. Geophysical Journal International, 138(3):727-746.
[16] Shapiro N M, Campillo M, Stehly L, et al. 2005. High-resolution surface-wave tomography from ambient seismic noise. Science, 307(5715):1615-1618.
[17] Tanimoto T, Alvizuri C. 2006. Inversion of the HZ ratio of microseisms forS-wave velocity in the crust. Geophysical Journal International, 165(1):323-335.
[18] Tanimoto T, Rivera L. 2008. The ZH ratio method for long-period seismic data:sensitivity kernels and observational techniques. Geophysical Journal International, 172(1):187-198.
[19] Tanimoto T, Yano T, Hakamata T. 2013. An approach to improve rayleigh-wave ellipticity estimates from seismic noise:Application to the Los Angeles Basin. Geophysical Journal International, 193(1):407-420.
[20] Xu G M, Yao H J, Zhu L B, et al. 2007. Shear wave velocity structure of the crust and upper mantle in western China and its adjacent area. Chinese J. Geophys.(in Chinese), 50(1):193-208.
[21] Xu Q, Zhao J M, Cui Z X, et al. 2009. Structure of the crust and upper mantle beneath the southeastern Tibetan Plateau by P and S receiver functions. Chinese J. Geophys. (in Chinese), 52(12):3001-3008.
[22] Yano T, Tanimoto T, Rivera L. 2009. The ZH ratio method for long-period seismic data:Inversion for S-wave velocity structure. Geophysical Journal International, 179(1):413-424.
[23] Yao H J, Beghein C, Van Der Hilst R D. 2008. Surface wave array tomography in SE Tibet from ambient seismic noise and two-station analysis-Ⅱ. Crustal and upper-mantle structure. Geophysical Journal International, 173(1):205-219.
[24] Yao H J, Van Der Hilst R D, Montagner JP. 2010. Heterogeneity and anisotropy of the lithosphere of Se Tibet from surface wave array tomography. Journal of Geophysical Research, 115(B12), doi:10.1029/2009JB007142.
[25] 胡家富, 朱雄关, 夏静瑜等. 2005. 利用面波和接收函数联合反演滇西地区壳幔速度结构. 地球物理学报, 48(5):1069-1076.
[26] 黄建平, 崇加军, 倪四道. 2008. 利用h-Kappa方法反演中国地区台站下地壳厚度. 中国科学技术大学学报, 38(1):33-40.
[27] 刘启元, 李昱, 陈九辉等. 2010. 基于贝叶斯理论的接收函数与环境噪声联合反演. 地球物理学报, 53(11):2603-2612.
[28] 徐果明, 姚华建, 朱良保等. 2007. 中国西部及其邻域地壳上地幔横波速度结构. 地球物理学报, 50(1):193-208.
[29] 徐强, 赵俊猛, 崔仲雄等. 2009. 利用接收函数研究青藏高原东南缘的地壳上地幔结构. 地球物理学报, 52(12):3001-3008.