Quantitative evaluation of modeling error in Lg-wave attenuation model
-
摘要: 对Lg波衰减模型中建模误差的统计特征进行了系统研究,并建立了地壳二维Lg波衰减模型。由于Lg波振幅可能受到几何扩散函数的强烈影响,合理评估反演过程的误差对于能否使用最小二乘意义下的反演非常重要。通过在川滇及其邻近地区收集的建模误差样本,使用K-S数值检验方法、Q-Q图和正态分布图形检验方法对Lg波衰减层析成像反演的输入数据中建模误差的分布特征进行了统计分析。采用奇异值分解(SVD)和反投影方法,分别获得了川滇地区的QLg模型,定量计算模型的协方差矩阵和分辨率矩阵,定量评估了QLg模型中每个格点的分辨率和误差。结果表明:在一阶近似条件下建模误差服从正态分布;通过开发的数据筛选程序,可以产生一个接近完美正态分布的数据集;与反投影方法相比,利用SVD方法获得的地壳Q值的分辨率更高;在射线覆盖较好的区域,QLg模型的分辨率达到100 km,相对误差小于3%。Abstract: In this study, we performed comprehensive statistical studies for the characteristics of the modeling error in Lg-wave attenuation model and obtained a two-dimensional crustal Lg-wave attenuation model. The amplitudes of Lg-wave can be strongly influenced by the geometrical spreading function and it is essential to reasonably evaluate the error of the inversion process if we use the inverted method under the meaning of least squares. Based on the modeling error samples collected in the Sichuan-Yunnan region and its adjacent areas, we use three statistical test methods, namely the K-S test, the Q-Q plot, and the normal distribution plot, to analyze the distribution characteristics of modeling errors in the input data of Lg- wave attenuation tomography inversion. We used both the SVD-based tomographic method and the back-projection method to invert for the Lg Q model of the Sichuan-Yunnan region and its adjacent area respectively. Then, we calculated the covariance matrix and the resolution matrix of the model and evaluated the resolution and error of each grid point in the Lg Q model quantitatively. Our results indicate that the modeling errors in the input data obeyed a normal distribution under first-order approximation. By using the data screening technique, we generated a new dataset with nearly perfect normally distributed for Q tomography. Compared with the result of back-projection, the SVD-based method could allow for obtaining a finer resolution of crustal Q value in the areas with dense ray path coverages, where the model resolution can be resolved within a range of 100 km with a relative error of less than 3%.
-
引言
地震Lg波对地壳的横向变化非常敏感,在稳定的地壳区域,如克拉通和地盾等,Lg波可以传播至3 000 km以外。但当Lg波穿过大陆地壳较薄区、地壳厚度突变区、巨厚的沉积盆地或存在部分熔融物质的地区时,Lg波会呈现出明显的衰减差异(Nuttli,1980;Kennett,1986;Mitchell,Hwang,1987;Bostock,Kennett,1990;Mitchell et al,1997;Rodgers et al,1997;Schurr et al,2003;Xie et al,2004).例如在青藏高原北部的强烈衰减区,Lg波的传播不超过700 km (McNamara et al,1996;Xie,2002)。地震波振幅衰减的多少可以用地下介质品质因子Q值来衡量,它是了解地下裂隙的数量、孔隙密度与分布,孔隙中流体含量以及地壳温度的重要参数。例如青藏高原的Lg波Q值极低(约为100—200),Chen和Xie (2017)认为这可能是青藏高原地壳内的高温造成的。由于构造因素导致的Q值变化比速度更加敏感,因此借助Lg波对Q值进行研究可以为了解地下介质提供更多的约束和解释(Akinci et al,1995;Rietbrock,2001;Gung,Romanowicz,2004)。
Q值层析成像是一个非常复杂的问题,因为Q值反演所使用的地震波振幅数据不仅受有限的Q值的影响,而且还受三维速度结构中波前扩展所导致的未知的三维几何扩散的影响(Nolet,1987)。在反演Q值时,一般先用简化过的一维几何扩散修正观测到的地面运动振幅,然后再取对数,将其结果作为Q值反演的输入数据,所得到的输入数据称为“折衷振幅数据”(Menke et al,2006;Chen,Xie,2017)。即使通过一维几何扩散修正,目前未知的三维几何扩散影响依旧保留在输入数据中。Chen和Xie (2017)将这些残留的三维几何扩散的影响当作数据误差来处理并称其为建模误差,相较于振幅测量不准导致的误差,建模误差才是Q值数据的主要误差来源(Xie,1998)。本研究使用双台法测量Q值,并将两个台站相对地震的方位角差也引入了建模误差。
在许多地球物理反演研究中,直接运用中心极限定理隐含地假设输入数据的误差呈正态分布,反演问题就能简化至用最小二乘算法求解。中心极限定理指出:如果观测数据的误差来源很多,无论某一种来源的误差呈何种形态分布,所有来源误差的综合效应可以假设呈正态分布(盛骤等,2008)。如果一组数据在对数域呈正态分布,那么数据本身就可能不是正态分布,所以在使用振幅数据的Q值反演中,对于数据的正态性检验不可忽略。Chen和Xie (2017)首次开展了对“折衷振幅数据”建模误差的统计分析工作,使用中国东南地区的Lg波,得出“折衷振幅数据”的建模误差可以用正态分布描述的结论,此结论是仅适用于当前研究区或所使用的数据,还是一个具有普适性的结论,尚待商榷。本文将使用在川滇及其周边地区收集到的更为密集的地震数据,对Lg波的降幅数据误差进行更详尽的统计分析工作。
另一方面,地球物理层析成像方法主要包括同时迭代重构法(simultaneous iterative reconstruction technique)、反投影法(back-projection)和最小平方QR分解法(least-squares QR)等。当Q值层析成像所使用的振幅数据符合最小二乘算法的要求时,这些算法就能够被广泛使用在Q值层析成像中(Pei et al,2006;Xie et al,2006;Hearn et al,2008; Pasyanos et al,2009;Zhao et al,2010,2013;Zhao,Xie,2016)。这些方法都是利用一个简单的初始Q值模型,使用交互式的方式改进模型以使数据得到更好的拟合,从而使反演快速收敛,得到令数据残差最小化的Q值最终模型,然后再通过“棋盘格”测试(Gök et al,2003;Zhao et al,2010,2013;Liang et al,2014)、合成“点扩散函数”(Dai et al,2020)或者反映射线路径的采样密度的区域(Zhao ,Xie,2016)等各种分辨率测试方法来找出这些模型的近似分辨率和近似误差。需要指出的是,这些方法并不能定量地给出模型的分辨率和误差值,但是在实际应用这些模型时,定量的模型分辨率和误差是非常重要的。例如,对模型的分辨率和误差定量分析的结果有助于推断出所得到的Q模型的可靠部分,从而通过可靠的Q模型恢复特征来推断地壳属性。尽管理论上可以量化Q模型的分辨率和误差,并检查它们之间的影响(Backus,Gilbert,1970),但实际应用中它们会受到计算能力的限制。大规模矩阵奇异值(singular value decomposition,缩写为SVD)分解的数值计算需要使用迭代方法,但由于有限的计算精度,所得到的奇异矢量在迭代过程中会逐渐失去正交性,从而使最终的SVD分解是不完整的,即并非所有的奇异矢量都还保持正交。如果在迭代过程中保持所有奇异矢量的正交性,则需要庞大的计算量和内存。Chen和Xie (2017)开发了一种基于PROPACK算法的新的2D层析成像方法,该算法的优点在于它在迭代过程中只针对失去正交性的奇异矢量做调整,由于失去正交性的只是部分奇异矢量,因此计算量和内存需求大量减少.而正因为PROPACK算法得到的所有奇异矢量都保持了正交性,所以它们不仅可以用来计算模型参数,还可以计算模型的分辨率矩阵和协方差矩阵(Larsen,1998)。
本文拟利用双台法得到川滇及其周边地区地震数据的”折衷振幅数据”,并对其建模误差统计样本进行了提取和分析,采用 K-S (kolmogorov-Smirnov),Q-Q (Quartile-Quartile)图和正态分布图等三种统计检验方法得到研究区建模误差的统计特征,以期验证Q值反演数据中建模误差统计特征的普适性。在证明建模误差近似服从正态分布的基础上,使用最小二乘法求解Q值层析成像问题,并采用基于PROPACK 算法的 2D 层析成像方法得到模型的参数、分辨率矩阵与协方差矩阵,以期得到定量的模型分辨率和误差。最后,采用不同的层析成像方法来求解Q值模型,通过所得到模型结果的异同来检验新层析成像方法的有效性和可靠性。
1. 数据和方法
1.1 数据筛选
川滇地区位于华南地块、印度板块和青藏高原的交会处(95°E—108°E,22°N—35°N),区域内深大断裂纵横交错,地震活动强烈。川滇地区也是印度板块和欧亚板块碰撞带东南缘的一部分及特提斯—喜马拉雅造山系的转折点(邓起东等,2002),川滇地区的主要断裂带和地块分布如图1所示。本研究收集了来自国家数字测震台网数据备份中心(郑秀芬等,2009)提供的2007年7月至2013年5月期间262个固定台站(图1b)记录到的地震波形数据。台站分布在川滇及其周边地区(图1a),挑选了震中距范围在200—3 000 km、震级在M3.5—6.5之间且具有清晰的P波初至和Lg波波包发育的地震事件资料用于Lg波Q值的计算。通过手动拾取震相后,再进行去均值、去趋势、滤波等数据处理,选取133次地震的24 001条地震波垂直分量波形数据用来计算Lg波频谱。Wei和Zhao (2019)在川滇地区分别对频率为0.5,1.0和1.5 Hz的Lg波进行反演,得出在该区域内
$ {Q}_{\mathrm{L}\mathrm{g}} $ 没有很强的频率依赖性。所以在本研究中,我们只分析了1 Hz附近的Lg谱幅值数据,使用一个中心频率为1.0 Hz的窄带宽(0.5—1.5 Hz)来获得平均谱。图 1 研究所用地震事件震中(a)和研究区台站(b)分布图图中红色矩形框代表研究区域和台站所在区域;黑线为主要地块分界线,灰线为主要断裂带,下同Figure 1. Distribution of epicenters of events (a) and seismic stations (b) in this studyRed rectangular represents the extent of the study area;the black solid lines indicate the major tectonic block boundaries and the grey dashed lines indicate faults,the same below1.2 Q值的测量方法
在特定频率
$ f $ 处观测到的Lg波振幅可以通过正演建模得到,即:$$ A ( f ) = S ( f ) G ( \varDelta ) \exp \left( { - \frac{{{\pi }f\varDelta }}{{vQ ( f ) }}} \right)R ( f ) \text{,} $$ (1) 式中
$ ,A ( f ) $ 表示振幅项,$ S ( f ) $ 为震源谱(Sereno et al,1988),$ G ( \mathrm{\varDelta } ) $ 是几何扩散项,$\varDelta$ 表示震中距,$ v $ 表示Lg波的群速度,$ Q ( f ) $ 为衰减值,$ R ( f ) $ 为台站的场地响应项(Street et al,1975;Xie,1993;Yang,2002;Chen,Xie,2017)。真实的几何扩散与3D地壳结构有关,但一维几何扩散可以简化为震中距$ \varDelta $ 的函数(Street et al,1975),表示为$$ G ( \varDelta ) = \left\{ {\begin{array}{*{20}{l}} {\dfrac{1}{\varDelta }} \,\;\;\;\quad\qquad\qquad \varDelta < {\varDelta _0},\\ {\left({\dfrac{1}{{{\varDelta _0}}}} \right){{\left( {\dfrac{{{\varDelta _0}}}{\varDelta }} \right)}^m}} \qquad \varDelta {\text{≥}} {\varDelta _0}{\text{.}} \end{array}} \right. $$ (2) 对于Lg波来说,一维几何扩散的幂指数取0.5,临界距离
$ {\varDelta }_{0} $ 设为100 km。Q值是通过指数函数影响振幅的,经过1D几何扩散校正并取对数后得到$$ \ln {{A_i} ( f ) } - \ln {G ( {{\varDelta _i}} ) } = \ln {{S_{ i}} ( f ) } - \frac{{{\pi }f{\varDelta _i}}}{{v{Q_i} ( f ) }} + \ln {{R_i} ( f ) } \text{,} $$ (3) 式(3)等号左边是反演数据,称为“折衷振幅数据”。本研究使用双台法 (two station based method,缩写为TS),将式(1)应用于与震源位于同一大圆弧的两个台站,因为这两个台站记录了同一次地震,而路径衰减的差别主要体现在两个台站之间,所以两个台站的振幅比直接消除了震源谱,并求出两个台站之间的Q值(Xie,Mitchell,1990;Xie et al,2006;Zor et al,2007;Gallegos et al,2014;Chen,Xie,2017)。但传统双台法不能消除场地效应的影响,只是假设两个台站的场地效应相同。为了避免必须进行Q值和震源参数之间的权衡所导致的Q模型的非唯一性 (Menke et al,2006),研究中又使用了逆双台法(reversed two-station method,缩写为RTS)(Chun et al,1987)和逆双事件法(reversed two-event method,缩写为RTE)(Chen,Xie,2017)。逆双台法是在双台法的基础上,从双台路径的对拓方向再引进一次地震,两个台站记录的两个地震的交叉振幅比,可以将两个台站的场地效应也直接消除,从而得到两个台站之间的Q值。逆双事件法是逆双台法的拓展,即将逆双台法中台站和地震的位置互换就是逆双事件法,它可以测量两个地震事件之间的Q值。
2. 数据误差的统计分析
2.1 建模误差的统计方法
双台法的Q值计算表示为
$$ \frac{{\tilde \varDelta }}{{{Q_{\rm{e}}} ( f ) }} = \frac{{v \ln \tilde A}}{{{\pi }f}} \text{,} $$ (4) 式中,
$ {Q}_{\mathrm{e}} ( f ) $ 为进行一维几何扩散校正后包含误差的$ Q ( f ) $ 的估计值,$\tilde{{\varDelta }}$ 为两个台站之间的有效距离,$ \tilde{A} $ 为双台法的振幅项与几何扩散项乘积的比值。本研究中将建模误差的定义范围进行了拓展,主要包括:① 残留的三维几何扩散导致的误差;② 不能完全去除台站响应部分造成的误差;③ 由Sn尾波引起的“噪声”和地震事件前无法模拟的背景噪声造成的误差;④ 两个台站和地震的方位角差别导致的误差。$ {Q}_{{\rm{e}}} ( f ) $ 中包含的误差可以表述为(Chen,Xie,2017):$$ \varDelta \left[ {\frac{{\tilde \varDelta }}{{Q ( f ) }}} \right] = \frac{{\tilde \varDelta }}{{{Q_e} ( f ) }} - \frac{{\tilde \varDelta }}{{Q ( f ) }} = \frac{v}{{{\pi }f}}\varDelta ( {\ln {\tilde A} } ) = \frac{v}{{{\pi }f}} ( 1 + \varepsilon ) {\text{≈}} \frac{v}{{{\pi }f}}\varepsilon \text{,} $$ (5) 式中,
$\mathrm{ln}\tilde{A}$ 为“折衷振幅数据”,假设误差很小,对$\mathrm{ln}\tilde{A}$ 进行一阶近似,${{\tilde \varDelta }}/{{Q}_{{\rm{e}}} ( f ) }$ 的误差$\varDelta [ {\tilde{\varDelta }}/{Q ( f ) } ] $ ,可以用一个随机变量$\varepsilon$ 表示。由于岩石圈在学者们感兴趣的研究频率范围(>0.1 Hz)是高度不均匀的,所以由$ \varepsilon $ 表示${\tilde{\varDelta }}/{Q ( f ) }$ 中的误差在不同的路径测量中是不相关的,即使是从同一个台站对测量收集的数据,它们的方位角(15°以内)也会略有变化.即使不同的震源位于台站相同的方位角上,仅仅只是震中距不同,也会因为局部波场的组成成分的差异导致产生不相关性(Mitchell,1995;Romanowicz,1998;Ji et al,2005;Gibbons et al,2010)。实验证明,使用不同地震与两个台站的组合,重复测量的两个台站之间Lg波的Q值的随机扰动平均值是没有任何规律的(Chen,Xie,2017),所以在研究过程中,对每一次${\tilde{{\varDelta }}}/{{Q}_{{\rm{e}}} ( f ) }$ 测量都提供一个随机变量$ \varepsilon $ 的样本。假设${\tilde{{\varDelta }}}/{{Q}_{{\rm{e}}} ( f ) }$ 的平均值接近真实的${\tilde{\varDelta }}/{Q ( f ) }$ ,可以用平均值减去单个样本值来近似的估计$\varDelta [ {\tilde{\varDelta }}/{Q ( f ) } ] $ ,从而根据式(4)得到单个$ \varepsilon $ 的值。2.2 建模误差样本采样
在计算1 Hz的频谱后,就可以求出“折衷振幅数据”
$ \mathrm{ln}\tilde{A} $ ,用此数据就可以估计两个台站间的$ {\tilde{\varDelta }}/{Q ( f ) } $ ,其中频率$f$ 设为1 Hz。在研究中分别得到了有效的9 633条双台法(TS)路径、3927条逆双台法(RTS)路径、2888条逆双事件(RTE)路径用于$ {\tilde{\varDelta }}/{Q ( f ) } $ 测量,其中有3 356条TS路径,722条RTS路径和302条RTE路径的唯一路径,在这些基于双台法的唯一路径中,分别有1 943条TS路径,547条RTS路径和253条RTE路径被多次采样,从而使得${\tilde{\varDelta }}/{Q ( f ) }$ 被重复测量.由于逆双台法和逆双事件法对台站和事件位置的要求比双台法更严格,使得可使用的数据大大减少,所以逆双台法和逆双事件法的有效路径相比于双台法要少很多。图2展示了三种方法的各个路径上
${\tilde{{\varDelta }}}/{{Q}_{{\rm{e}}} ( f ) }$ 的重复测量的次数分布图,最大重复测量次数达到100次以上.通过${\tilde{\varDelta }}/{Q ( f ) }$ 的重复测量来收集$ \varepsilon $ 样本,分别用$ {\varepsilon }_{\mathrm{T}\mathrm{S}} $ ,$ {\varepsilon }_{\mathrm{R}\mathrm{T}\mathrm{S}} $ 和$ {\varepsilon }_{\mathrm{R}\mathrm{T}\mathrm{E}} $ 来表示用TS,RTS和RTE方法估计的$ \varepsilon $ 值,如表1所示对于三种方法测量${\tilde{\varDelta }}/{Q ( f ) }$ 获得的$ {\varepsilon }_{\mathrm{T}\mathrm{S}} $ ,$ {\varepsilon }_{\mathrm{R}\mathrm{T}\mathrm{S}} $ 和$ {\varepsilon }_{\mathrm{R}\mathrm{T}\mathrm{E}} $ 的样本数量分别为8 220,3 652和2 839。表 1 使用TS,RTS和RTE方法估计ε的详细信息Table 1. Estimated ε information using TS, RTS and RTE methods方法 $ {\varepsilon } $样本的
数量误差在两个标准差
内的$ \varepsilon $占比ECDF和NCDF的
相关系数K-S测试结果 标准差
std$ ( \varepsilon ) $方差Var$ ( \varepsilon ) $ P值 显著性水平$ \alpha $ TS 8 220 87.09% 99.99% 在置信度为95%下接受零假设 0.158 5 0.025 1 0.693 1 0.05 RTS 3 652 88.12% 99.98% 在置信度为95%下接受零假设 0.081 6 0.006 7 0.324 0 0.05 RTE 2 839 93.59% 99.64% 在置信度为95%下接受零假设 0.104 9 0.011 0 0.280 1 0.05 2.3 统计学方法检验结果
原始的
$ {\varepsilon }_{\mathrm{T}\mathrm{S}} $ ,$ {\varepsilon }_{\mathrm{R}\mathrm{T}\mathrm{S}} $ 和$ {\varepsilon }_{\mathrm{R}\mathrm{T}\mathrm{E}} $ 的样本分布具有中心峰值的特征,是显著的正态分布的特征,如图3左侧所示,然而远离峰值的分布在正负方向上都有较长的尾部,这些分布与Nolet (1987)讨论的远震P波走时残差分布非常相似;$ \varepsilon $ 的特征也符合Chen和Xie (2017)的统计特征结果。图 3 建模误差$\varepsilon$ 样本分布(左)和K-S检测图(右)(a) TS方法;(b) RTS方法;(c)RTE方法。左图中两条黑色直线标记了数据筛选准则|ε-median(ε)|>2std(ε)的区间,红色线表示最佳的正态概率密度函数值;右图中灰色线代表着95%置信区间的CDF的上下界限,绿色小线段代表着K-S检验的经验CDF和理论CDF之间的最大垂向差值Figure 3. Distribution of samples of modeling error$\varepsilon$ (left)and results of K-S test (right)(a) TS method;(b) RTS method ;(c) RTS method. The black lines in the left figure mark the interval of |ε-median(ε)|>2std(ε),which we use as criteria for data screening and the best-fit theoretical normal PDF of the screened samples is plotted over the distribution map with the red line. The grey lines in the right figure indicate lower and upper 95 percent confidence bounds for the CDF and the green segment indicates the maximum distance from the empirical CDF to the theoretical normal CDF which is used in the K-S test在对
$ \varepsilon $ 的样本进行了几次筛选试验后,确定了任何$ \varepsilon $ 偏离中位数超过两个$ \varepsilon $ 的标准差的样本都被剔除的筛选准则,即$\left| {\varepsilon - {\text{median}} ( \varepsilon ) } \right| > 2 {{{\rm{std}} ( }}\varepsilon ) $ ,这也符合统计学上的“68-95-99.7”准则(Pukelsheim,1994),68.27% (n=1),95.45% (n=2)和99.73% (n=3)的值是分别距离样本数据的平均值小于一个、两个和三个标准差内的百分比,超过三个标准差的数据可能是异常值,并且如果有许多数据超过三个标准差,则有理由质疑假定的正态性(Wheelan,2013)。两个标准差范围代表着95%的置信水平,可以同时平衡样本选择的精度和准确度。如表1所示,经过数据筛选后,对于所有的样本$ \varepsilon $ ,剔除的数据量不超过总体的12%。方差$ {\rm{Var}} ( {\varepsilon }_{\mathrm{T}\mathrm{S}} ) ,{\rm{Var}} ( {\varepsilon }_{\mathrm{R}\mathrm{T}\mathrm{S}} ) $ 和$ {\rm{Var}} ( {\varepsilon }_{\mathrm{R}\mathrm{T}\mathrm{E}} ) $ 分别为0.025 1,0.006 7和0.011 0,$ {\varepsilon }_{\mathrm{T}\mathrm{S}} $ 的方差远大于$ {\varepsilon }_{\mathrm{R}\mathrm{T}\mathrm{S}} $ 和$ {\varepsilon }_{\mathrm{R}\mathrm{T}\mathrm{E}} $ 的方差,说明$ {\varepsilon }_{\mathrm{T}\mathrm{S}} $ 的分布比$ {\varepsilon }_{\mathrm{R}\mathrm{T}\mathrm{S}} $ 和$ {\varepsilon }_{\mathrm{R}\mathrm{T}\mathrm{E}} $ 更广,$ {\rm{Var}} ( {\varepsilon }_{\mathrm{R}\mathrm{T}\mathrm{S}} ) $ 和$ {\rm{Var}} ( {\varepsilon }_{\mathrm{R}\mathrm{T}\mathrm{E}} ) $ 之间的微小差异很可能是由每个方法收集的重复测量样本数不同造成的,$ {\rm{Var}} ( {\varepsilon }_{\mathrm{T}\mathrm{S}} ) $ 的值最大可以归因于TS方法在测量中不能消除场地效应的影响。K-S (Kolmogorov-Smirnov)检验是一种非参数测试,测试分布的正态性时,将已知的正态概率分布与数据生成的分布(经验分布函数)进行量化比对(Massey,1951;Vrbik,2020),研究中的零假设为数据是正态分布。本研究采取单样本K-S检验, p值为否定零假设的最低显著性水平(Vrbik,2020)。K-S检验通过p值与显著性水平α进行比较,从而来决定是否拒绝零假设.如表1所示,p值都大于显著性水平,所以在95%的置信水平下接受零假设,建模误差是正态分布。剔除异常数据中新的ε样本的K-S检验图如图3,右侧图所示。使用筛选后的
$ {\varepsilon }_{\mathrm{T}\mathrm{S}} $ 样本计算累积分布函数(cumulation distribution function,缩写为CDF)如图3a右图所示,并与正态分布的CDF进行比较,观测和筛选的经验累积分布函数(empirical cumulation distribution function,缩写为ECDF)与理论最佳正态累积分布函数(normal cumulation distribution function,缩写为NCDF)之间的相关系数是99.99% (见表1),通过K-S检验不仅证明了筛选准则的可靠性,还证明了筛选后的$ {\varepsilon }_{{\rm{TS}}} $ 样本非常接近于正态分布函数。$ {\varepsilon }_{\mathrm{R}\mathrm{T}\mathrm{S}} $ 和$ {\varepsilon }_{\mathrm{R}\mathrm{T}\mathrm{E}} $ 的样本的检测结果如图3b,c右图所示,也得出了相同的结论。研究中还使用Q-Q (Quartile-Quartile)分位数图来检测
$ \varepsilon $ 是否服从正态分布。Q-Q分位数图是一种图形统计检验方法,用此方法可以找出两组数据是否来自同一分布。如果在Q-Q分位数图中纵轴表示理论正态分布分位数,横轴上输入验证数据,拟合的数据如果聚类在$y = x$ 的线上,这表明样本数据呈正态分布(Marden,2004;Wheelan,2013)。如图4a所示,$ {\varepsilon }_{\mathrm{T}\mathrm{S}} $ 的Q-Q分位数图显示$ {\varepsilon }_{\mathrm{T}\mathrm{S}} $ 与正态分布基本一致,仅左下角和右上角存在一定偏差,代表着样本中剔除了长尾异常,即尾部未采样,但总体上可以证明$\varepsilon $ 服从正态分布。$ {\varepsilon }_{\mathrm{R}\mathrm{T}\mathrm{S}} $ 和$ {\varepsilon }_{\mathrm{R}\mathrm{T}\mathrm{E}} $ 的样本的检测结果如图4b,c所示,也得出了相同的结论。研究中还使用了其它图形检验方法如正态概率图以及数值检验方法如Shapiro-Wilk测试应用于检验筛选的
$\varepsilon $ 样本,所有的这些方法都是为了检验随机样本在一定置信水平下符合正态分布的假设而设计的,并互相验证测试结果。本研究中的检验结果与Chen和Xie (2017)开展的Q值模型建模误差的统计特征一致,但本次研究中使用K-S检验、Q-Q分位图并没有像Chen和Xie (2017)的结果一样对TS和RTS方法的正态性检测失效,因为这些方法是为相对较小(不超过几千)的样本量设计的,川滇研究区的Q值模型的建模误差样本都在几千的范围内,可以很好地呈现出正态性。所以通过多种统计检验方法,都证明了本文所采用的数据筛选准则是合适的,可以剔除建模误差样本的异常值,产生一个接近于完美正态分布的新的数据集,从而使用最小二乘算法求解Q值模型。3. 川滇地区Q值模型及评价
3.1 Q值层析成像方法
研究区Lg波Q值的横向变化可以用二维衰减层析成像解决,将研究区域划分为M个小单元来参数化Q值的横向变化,假设每个单元中的Q值恒定. 在收集N条唯一路径上的
${\tilde{\varDelta }}/{Q ( f ) }$ 测量值之后,对每条路径都可以构建线性方程(Chen,Xie,2017)$$ \frac{{{\varDelta _n}}}{{{Q_{{n}}} ( f ) }} = \sum\limits_{m = 1}^M {\frac{{{\varDelta _{nm}}}}{{{Q_m} ( f ) }}}, $$ (6) 式中,n=1,2,···,N代表着第n条射线路径,m=1,2,···,M代表着第m个模型的小单元。
$ {Q}_{n} ( f ) $ 是第n条射线路径上的Q($ f $ )的估计值,第n条射线路径长度为$ {{\varDelta }}_{n} $ ,当它穿过研究区域的第m个小单元的时候,相交的小段距离长度为$ {{\varDelta }}_{nm} $ 。式(6)的矩阵形式为$$ {{\boldsymbol{d}}_N} = {{\boldsymbol{A}}_{N {\text{×}} M}}{{\boldsymbol{m}}_M} $$ (7) 式中,数据矢量d中的元素为
$ {{\varDelta }_{n}}/{{Q}_{n} ( f ) } $ ,转换矩阵中的元素为$ {\varDelta }_{nm} $ ,而模型参数矢量m中的元素为$ {1}/{{Q}_{m} ( f ) } $ 。转换矩阵A是一个稀疏矩阵,因为任意第n条射线路径穿过单元格数量比总的单元格数量M少得多。另外,实际反演的模型值是1/Q,只是按照习惯在图像中转换为Q。一旦证明了Q值反演的数据集受正态分布的误差影响,则可以使用最小二乘法来反演1/Q 值模型。数据的方差可以由
$ \varepsilon $ 的方差得出:$$ {\rm{Var}}\left( {\frac{{ \tilde{\varDelta} }}{{{Q_{\rm{e}}} ( f ) }}} \right) = \frac{{{V^2}}}{{{{\pi }^2}{f^2}}}{\rm{Var}} ( \varepsilon ) {\text{.}}$$ (8) 为了更好地反映通过改变
${{\Delta }_{n}}/{{Q}_{n} ( f ) }$ 测量值之间的方差来量化模型的不确定性,通过式(8)可以计算基于双台方法的任何特定路径上的${\rm{Var}} [ {\tilde{{\varDelta }}}/{{Q}_{{\rm{e}}} ( f ) } ] $ 值。在研究中引入一个$ N{\text{×}} N $ 的数据协方差矩阵${\boldsymbol{C}}$ ,它的第n个对角线元素由${\rm{Var}} [ {{\varDelta }_{n}}/{{Q}_{n} ( f ) } ] $ 给出,然后归一化得到下式:$$ {{\boldsymbol{C}}^{ - \frac{1}{2}}}{\boldsymbol{A}}m = {{\boldsymbol{C}}^{\frac{1}{2}}}{\boldsymbol{d}} {\text{.}}$$ (9) 最小二乘意义下求解模型
$ \mathit{m} $ 需要求解矩阵A的逆矩阵,但是在层析成像中这通常需要很大的计算量和内存,取而代之的是各种迭代技术,如反投影法(Xie et al,2004,2006)。反投影方法主要着力于迭代更新一个测试解$ \mathit{m} $ ,直到预测的数据不匹配而停止改善。研究中使用的新的2D层析成像方法不是迭代技术,而是通过PROPACK算法(Larsen,1998)对矩阵A进行精确的SVD分解,从而求得矩阵A的逆矩阵。3.2 川滇地区Q值模型和分辨率
可视化的棋盘格测试几乎是所有的层析成像方法均采用的一种评估反演可靠性的方法。为了测试研究使用的Lg波频谱振幅数据集能否很好地恢复Lg波衰减特征,使用棋盘格测试作为检测模型的可靠性方法之一,如图5所示。测试了不同网格尺寸的模型参数化范围,最终选择1°×1°作为网格反演大小。周连庆等(2008)研究表明,川滇地区最高Q值不超过800,最低Q值不小于100,因此设置了Q值交错为100和800的棋盘格来生成合成数据。合成数据的路径位置和数量与实际数据相同,然后通过合成数据来恢复原始的棋盘图案来评估反演的质量。基于SVD的新的层析成像方法可以完整地恢复出频率为1 Hz的
$ {Q}_{{\rm{Lg}}} $ 模型的棋盘格,并且恢复的棋盘格图案与射线覆盖范围(图6a)一致。尽管棋盘格测试说明反演过程是稳定可靠的,但它对模型分辨率大小的说明是相对模糊的,仅能说明射线覆盖的好坏,但是密集的射线覆盖并不一定意味着模型的高分辨率。一种极端情况是,如果模型地区仅仅被同一方向的射线覆盖,而缺乏交叉方向的射线,那么模型则不会被很好地解析。
新的层析成像方法同时给出了模型的分辨率矩阵,可以通过构建Backus-Gilbert扩散函数(spread function,缩写为SF)来定量评价模型的分辨率(Chen,Xie,2017),通过分辨率矩阵的元素所控制的均方根距离,可以对模型的有限分辨率进行近似度量,从而给出每个反演格点的可分辨范围。图6b显示出本研究中通过Backus-Gilbert扩散函数对川滇地区Q模型的分辨率进行评估,在地震射线覆盖良好的研究区域中,分辨精度可以达到100 km。通过棋盘格测试和Backus-Gilbert扩散函数都证明了研究中得到的川滇地区Lg波的1/Q模型是可靠的。
基于SVD的新的层析成像方法和反投影方法获得的川滇地区的地壳Q值模型如图7所示。由于两种方法原理不同,但都需要做光滑处理,为方便比较,两种方法采用了相同的光滑参数.在相同的光滑处理下,不同方法的
$ {Q}_{{\rm{Lg}}} $ 的区域变化似乎相当一致,最大值大于700,最小值小于80,横向非均匀性较强,模式具有与该地区速度结构相似的特征(如Tian 等(2021))。由图7可以看出,川滇地区西部及其邻近地区Q值普遍低于东部,四川盆地出现明显的东西向变化与Chen和Xie (2017)的结果相似,同时盆地东部的高Q值与Zhao等(2013)估算的整个四川盆地相似。扬子克拉通的Q值很高(低衰减),松潘甘孜褶皱带Q值很低(高衰减)。本研究中得到的$ {Q}_{{\rm{Lg}}} $ 结果与Wei和Zhao (2019)的川滇地区1 Hz的$ {Q}_{{\rm{Lg}}} $ 的区域分布基本一致。在滇中地块西侧和松潘—甘孜褶皱带西侧,反投影法得到的结果似乎有更精细的结果,但是从分辨率结果可以看出,滇中地块的Backus-Gilbert 扩散函数值(SF)约在100 km,而松潘—甘孜褶皱带西侧以及三江断层西侧的SF值在200至300 km。所以在相同的分辨率下进行构造讨论,新方法相比于反投影方法能够体现出更精细的地壳横向变化结果。不同的反演方法涉及不同的参数选择,在大多数情况下参数的设置是经验决定的,两种方法的比较过程中尽可能使用类似的参数,所以细节方面不尽相同。影响层析成像结果的因素很多,例如反投影方法需要设置一个初始模型,不同的初始条件,反演的结果也可能不同。新方法可以通过PROPACK算法直接对转换矩阵进行SVD分解,所以不需要初始模型,其优越性表现在模型分辨率的定量分析,从而决定了模型所反映的结构的可靠性。3.3 模型的不确定性
通过SVD奇异值分解,得到了Lg波1/Q 模型的协方差矩阵。因为证实了数据误差服从正态分布,误差的大小完全可以由误差的方差表示(式(8))。反演过程中,数据包含的误差会传递到模型中,反映在模型的协方差矩阵中,从而由数据误差导致的1/Q模型的误差可以由模型的协方差矩阵完全量化。Q值层析成像模型的不确定度由模型的分辨率和误差共同决定。一阶近似下,用δQ/Q来标定模型的相对误差,这个值可以使用协方差矩阵对角线元素的平方根来近似估计。在川滇地区,1 Hz的Q模型相对误差如图8所示,一个标准差下的模型的相对误差小于9%,在射线路径覆盖良好的区域,相对误差小于3%,这得益于研究中严格的数据筛选准则和基于SVD的新的层析成像方法。完整的Q值层析成像模型的协方差矩阵可以用来定量估计特定路径预测Q值的误差,当使用路径Q值校正振幅时,这些误差估计可用于建立地震事件识别和震源大小估计的置信水平。
4. 讨论与结论
使用川滇地区及其邻近地区2007年7月至2013年5月的地震数据,对Lg波的“折衷振幅数据”进行了严格的正态性检验工作,这项工作在应用最小二乘算法之前十分重要,但经常被忽略,尤其在应用类似“折衷振幅数据”这样被处理过的数据时。用一个随机变量
$ \varepsilon $ 来表示建模误差,将“折衷振幅数据”误差统计分析的复杂任务转化为更容易处理的对$ \varepsilon $ 样本的统计分析工作,获得了如下结论:1) Q层析成像使用“折衷振幅数据”作为输入数据导致的建模误差主要由几何扩散函数、地震事件前和相关相位前尾波的噪声以及场地效应的方位变化影响,通过多种统计检验方法证明,以前Q值成像研究中仅仅被隐含假设为正态分布的建模误差是真实服从正态分布的,尽管那些研究可能使用了不同的数据筛选方法。
$ \varepsilon $ 的分布以一个拟正态分布的峰为主,辅以向峰两侧扩散的弱长尾异常值。2) 通过恰当的数据筛选准则,可以剔除建模误差样本的异常值,从而产生一个接近于完美正态分布的新的数据集,用于在最小二乘意义下求解Q值模型。
3) 研究中所使用的川滇地区地震数据表明,方差
$ {\rm{Var}} ( {\varepsilon }_{\mathrm{T}\mathrm{S}} ) ,{\rm{Var}} ( {\varepsilon }_{\mathrm{R}\mathrm{T}\mathrm{S}} ) $ 和$ {\rm{Var}} ( {\varepsilon }_{\mathrm{R}\mathrm{T}\mathrm{E}} ) $ 分别为0.025 1,0.006 7和0.011 0,$ {\rm{Var}} ( {\varepsilon }_{\mathrm{T}\mathrm{S}} ) $ 值最大,这可以归因于TS方法在测量中不能消除场地效应的影响.4) 研究中使用了基于SVD分解的新的Q层析成像方法来反演Q值模型,相比于反投影等迭代方法,它最大优点是可以得到模型分辨率矩阵以及协方差矩阵,能够把模型的分辨率和误差定量化。
5) 新方法比反投影法可以获得更精细的地壳Q值横向结构,在射线覆盖良好的区域,新方法的QLg模型可分辨达到100 km左右,模型相对误差小于3%。研究区定量的不确定度结果表明,新的Q层析成像方法是稳定可靠的,可以用来解释地下介质的横向变化,同时通过使用模型的协方差矩阵,可以定量估计特定路径上Q预测的误差,这样的定量误差估计可以在未来继续开展在一定置信水平下的地震动预测、震源谱估计、震源类型识别和震源大小的估计。
国家测震台网数据备份中心为本文提供了数据支持,各位审稿专家为本文提出了宝贵的修改意见,作者在此一并表示感谢。
-
图 1 研究所用地震事件震中(a)和研究区台站(b)分布图
图中红色矩形框代表研究区域和台站所在区域;黑线为主要地块分界线,灰线为主要断裂带,下同
Figure 1. Distribution of epicenters of events (a) and seismic stations (b) in this study
Red rectangular represents the extent of the study area;the black solid lines indicate the major tectonic block boundaries and the grey dashed lines indicate faults,the same below
图 3 建模误差
$\varepsilon$ 样本分布(左)和K-S检测图(右)(a) TS方法;(b) RTS方法;(c)RTE方法。左图中两条黑色直线标记了数据筛选准则|ε-median(ε)|>2std(ε)的区间,红色线表示最佳的正态概率密度函数值;右图中灰色线代表着95%置信区间的CDF的上下界限,绿色小线段代表着K-S检验的经验CDF和理论CDF之间的最大垂向差值
Figure 3. Distribution of samples of modeling error
$\varepsilon$ (left)and results of K-S test (right)(a) TS method;(b) RTS method ;(c) RTS method. The black lines in the left figure mark the interval of |ε-median(ε)|>2std(ε),which we use as criteria for data screening and the best-fit theoretical normal PDF of the screened samples is plotted over the distribution map with the red line. The grey lines in the right figure indicate lower and upper 95 percent confidence bounds for the CDF and the green segment indicates the maximum distance from the empirical CDF to the theoretical normal CDF which is used in the K-S test
表 1 使用TS,RTS和RTE方法估计ε的详细信息
Table 1 Estimated ε information using TS, RTS and RTE methods
方法 $ {\varepsilon } $样本的
数量误差在两个标准差
内的$ \varepsilon $占比ECDF和NCDF的
相关系数K-S测试结果 标准差
std$ ( \varepsilon ) $方差Var$ ( \varepsilon ) $ P值 显著性水平$ \alpha $ TS 8 220 87.09% 99.99% 在置信度为95%下接受零假设 0.158 5 0.025 1 0.693 1 0.05 RTS 3 652 88.12% 99.98% 在置信度为95%下接受零假设 0.081 6 0.006 7 0.324 0 0.05 RTE 2 839 93.59% 99.64% 在置信度为95%下接受零假设 0.104 9 0.011 0 0.280 1 0.05 -
邓起东,张培震,冉勇康,杨晓平,闵伟,楚全芝. 2002. 中国活动构造基本特征[J]. 中国科学:D辑,32(12):1020–1030. Deng Q D,Zhang P Z,Ran Y K,Yang X P,Min W,Chu Q Z. 2003. Basic characteristics of active tectonics of China[J]. Science in China:Series D,46(4):356–372.
盛骤, 谢式千, 潘承毅. 2008. 概率论与数理统计[M]. 北京: 高等教育出版社: 119–126. Sheng Z, Xie S Q, Pan C Y. 2008. Probability Theory and Mathematical Statistics[M]. Beijing: Higher Education Press: 119–126 (in Chinese).
郑秀芬,欧阳飚,张东宁,姚志祥,梁建宏,郑洁. 2009. “国家数字测震台网数据备份中心”技术系统建设及其对汶川大地震研究的数据支撑[J]. 地球物理学报,52(5):1412–1417. Zheng X F,Ouyang B,Zhang D N,Yao Z X,Liang J H,Zheng J. 2009. Technical system construction of Data Backup Centre for China Seismograph Network and the data support to researches on the Wenchuan earthquake[J]. Chinese Journal of Geophysics,52(5):1412–1417 (in Chinese).
周连庆,赵翠萍,修济刚,陈章立. 2008. 川滇地区Lg波Q值层析成像[J]. 地球物理学报,51(6):1745–1752. doi: 10.3321/j.issn:0001-5733.2008.06.015 Zhou L Q,Zhao C P,Xiu J G,Chen Z L. 2008. Tomography of QLg in Sichuan-Yunnan zone[J]. Chinese Journal of Geophysics,51(6):1745–1752 (in Chinese).
Akinci A,Del Pezzo E,Ibáñez J M. 1995. Separation of scattering and intrinsic attenuation in southern Spain and western Anatolia (Turkey)[J]. Geophys J Int,121(2):337–353. doi: 10.1111/j.1365-246X.1995.tb05715.x
Backus G,Gilbert F. 1970. Uniqueness in the inversion of inaccurate gross Earth data[J]. Philosophical Transactions of the Royal Society of London. Series A:Mathematical and Physical Sciences,266(1173):123–192.
Bostock M G,Kennett B L N. 1990. The effect of 3-D structure on Lg propagation patterns[J]. Geophys J Int,101(2):355–364. doi: 10.1111/j.1365-246X.1990.tb06574.x
Chen Y L,Xie J K. 2017. Resolution,uncertainty and data predictability of tomographic Lg attenuation models:Application to southeastern China[J]. Geophys J Int,210(1):166–183. doi: 10.1093/gji/ggx147
Chun K Y,West G F,Kokoski R J,Samson C. 1987. A novel technique for measuring Lg attenuation-results from eastern Canada between 1 to 10 Hz[J]. Bull Seismol Soc Am,77(2):398–419. doi: 10.1785/BSSA0770020398
Dai A,Tang C C,Liu L,Xu R. 2020. Seismic attenuation tomography in southwestern China:Insight into the evolution of crustal flow in the Tibetan Plateau[J]. Tectonophysics,792:228589. doi: 10.1016/j.tecto.2020.228589
Gallegos A,Ranasinghe N,Ni J,Sandvol E. 2014. Lg attenuation in the central and eastern United States as revealed by the Earth Scope Transportable Array[J]. Earth Planet Sci Lett,402:187–196. doi: 10.1016/j.jpgl.2014.01.049
Gibbons S J,Kværna T,Ringdal F. 2010. Considerations in phase estimation and event location using small-aperture regional seismic arrays[J]. Pure Appl Geophys,167(4/5):381–399.
Gök R,Sandvol E,Türkelli N,Seber D,Barazangi M. 2003. Sn attenuation in the Anatolian and Iranian plateau and surrounding regions[J]. Geophys Res Lett,30(24):8042.
Gung Y,Romanowicz B. 2004. Q tomography of the upper mantle using three-component long-period waveforms[J]. Geophys J Int,157(2):813–830. doi: 10.1111/j.1365-246X.2004.02265.x
Hearn T M,Wang S Y,Pei S P,Xu Z H,Ni J F,Yu Y X. 2008. Seismic amplitude tomography for crustal attenuation beneath China[J]. Geophys J Int,174(1):223–234. doi: 10.1111/j.1365-246X.2008.03776.x
Ji C,Tsuboi S,Komatitsch D,Tromp J. 2005. Rayleigh-wave multipathing along the West coast of North America[J]. Bull Seismol Soc Am,95(6):2115–2124. doi: 10.1785/0120040180
Kennett B L N. 1986. Lg waves and structural boundaries[J]. Bull Seismol Soc Am,76(4):1133–1141.
Larsen R M, 1998. Lanczos Bidiagonalization With Partial Reorthogonalization[D]. Jutland: Aarhus University: 8–87.
Liang X F,Sandvol E,Kay S,Heit B,Yuan X H,Mulcahy P,Chen C,Brown L,Comte D,Alvarado P. 2014. Delamination of southern Puna lithosphere revealed by body wave attenuation tomography[J]. J Geophys Res:Solid Earth,119(1):549–566. doi: 10.1002/2013JB010309
Marden J I. 2004. Positions and QQ plots[J]. Statist Sci,19(4):606–614.
Massey F J. 1951. The Kolmogorov-Smirnov test for goodness of fit[J]. J Am Statist Associat,46(253):68–78. doi: 10.1080/01621459.1951.10500769
McNamara D E,Owens T J,Walter W R. 1996. Propagation characteristics of L
g across the Tibetan Plateau[J]. Bull Seismol Soc Am,86(2):457–469. doi: 10.1785/BSSA0860020457 Menke W,Holmes R C,Xie J K. 2006. On the nonuniqueness of the coupled origin time-velocity tomography problem[J]. Bull Seismol Soc Am,96(3):1131–1139. doi: 10.1785/0120050192
Mitchell B J. 1995. Anelastic structure and evolution of the continental crust and upper mantle from seismic surface wave attenuation[J]. Rev Geophys,33(4):441–462. doi: 10.1029/95RG02074
Mitchell B J,Hwang H J. 1987. Effect of low Q sediments and crustal Q on L
g attenuation in the United States[J]. Bull Seismol Soc Am,77(4):1197–1210. doi: 10.1785/BSSA0770041197 Mitchell B J,Pan Y,Xie J K,Cong L L. 1997. L
g coda Q variation across Eurasia and its relation to crustal evolution[J]. J Geophys Res:Solid Earth,102(B10):22767–22779. doi: 10.1029/97JB01894 Nolet G. 1987. Seismic Tomography: With Applications in Global Seismology and Exploration Geophysics[M]. Dordrecht: Springer: 1–24.
Nuttli O W. 1980. The excitation and attenuation of seismic crustal phases in Iran[J]. Bull Seismol Soc Am,70(2):469–485. doi: 10.1785/BSSA0700020469
Pasyanos M E,Matzel E M,Walter W R,Rodgers A J. 2009. Broad-band L
g attenuation modelling in the Middle East[J]. Geophys J Int,177(3):1166–1176. doi: 10.1111/j.1365-246X.2009.04128.x Pei S P,Zhao J M,Rowe C A,Wang S Y,Hearn T M,Xu Z H,Liu H B,Sun Y S. 2006. ML amplitude tomography in North China[J]. Bull Seismol Soc Am,96(4A):1560–1566. doi: 10.1785/0120060021
Pukelsheim F. 1994. The three sigma rule[J]. Am Statist,48(2):88–91.
Rietbrock A. 2001. P wave attenuation structure in the fault area of the 1995 Kobe earthquake[J]. J Geophys Res:Solid Earth,106(B3):4141–4154. doi: 10.1029/2000JB900234
Rodgers A J,Ni J F,Hearn T M. 1997. Propagation characteristics of short-period Sn and Lg in the Middle East[J]. Bull Seismol Soc Am,87(2):396–413.
Romanowicz B. 1998. Attenuation tomography of the earth’s mantle:A review of current status[J]. Pure Appl Geophys,153(2):257–272.
Schurr B,Asch G,Rietbrock A,Trumbull R,Haberland C. 2003. Complex patterns of fluid and melt transport in the central Andean subduction zone revealed by attenuation tomography[J]. Earth Planet Sci Lett,215(1/2):105–119.
Sereno T J Jr,Bratt S R,Bache T C. 1988. Simultaneous inversion of regional wave spectra for attenuation and seismic moment in Scandinavia[J]. J Geophys Res:Solid Earth,93(B3):2019–2035. doi: 10.1029/JB093iB03p02019
Street R L,Herrmann R B,Nuttli O W. 1975. Spectral characteristics of the L
g wave generated by central United States earthquakes[J]. Geophys J Int,41(1):51–63. doi: 10.1111/j.1365-246X.1975.tb05484.x Tian H Y,He C S,Santosh M. 2021. S-wave velocity structure of the Sichuan-Yunnan region,China:Implications for extrusion of Tibet Plateau and seismic activities[J]. Earth Space Sci,8(7):e2021EA001640.
Vrbik J. 2020. Deriving CDF of Kolmogorov-Smirnov test statistic[J]. Appl Math,11(3):227–246. doi: 10.4236/am.2020.113018
Wei Z,Zhao L. 2019. Lg Q model and its implication on high-frequency ground motion for earthquakes in the Sichuan and Yunnan region[J]. Earth Planet Phys,3(6):526–536.
Wheelan C. 2013. Naked Statistics[M]. London: W W Norton & Company: 19–292.
Xie J K. 1993. Simultaneous inversion for source spectrum and path Q using L
g with application to three semipalatinsk explosions[J]. Bull Seismol Soc Am,83(5):1547–1562. doi: 10.1785/BSSA0830051547 Xie J K. 1998. Spectral inversion of Lg from earthquakes:A modified method with applications to the 1995,western Texas earthquake sequence[J]. Bull Seismol Soc Am,88(6):1525–1537.
Xie J. 2002. Lg Q in the eastern Tibetan Plateau[J]. Bull Seismol Soc Am,92(2):871–876. doi: 10.1785/0120010154
Xie J,Mitchell B J. 1990. A back-projection method for imaging large-scale lateral variations of L
g coda Q with application to continental Africa[J]. Geophys J Int,100(2):161–181. doi: 10.1111/j.1365-246X.1990.tb02477.x Xie J,Gok R,Ni J,Aoki Y. 2004. Lateral variations of crustal seismic attenuation along the INDEPTH profiles in Tibet from LgQ inversion[J]. J Geophys Res:Solid Earth,109(B10):B10308.
Xie J,Wu Z,Liu R,Schaff D,Liu Y,Liang J. 2006. Tomographic regionalization of crustal Lg Q in eastern Eurasia[J]. Geophys Res Lett,33(3):L03315.
Yang X. 2002. A numerical investigation of Lg geometrical spreading[J]. Bull Seismol Soc Am,92(8):3067–3079. doi: 10.1785/0120020046
Zhao L F,Xie X B,Wang W M,Zhang J H,Yao Z X. 2010. Seismic Lg-wave Q tomography in and around northeast China[J]. J Geophys Res:Solid Earth,115(B8):B08307.
Zhao L F,Xie X B,He J K,Tian X B,Yao Z X. 2013. Crustal flow pattern beneath the Tibetan Plateau constrained by regional Lg-wave Q tomography[J]. Earth Planet Sci Lett,383:113–122. doi: 10.1016/j.jpgl.2013.09.038
Zhao L F,Xie X B. 2016. Strong Lg-wave attenuation in the Middle East continental collision orogenic belt[J]. Tectonophysics,674:135–146. doi: 10.1016/j.tecto.2016.02.025
Zor E,Sandvol E,Xie J,Türkelli N,Mitchell B,Gasanov A H,Yetirmishli G. 2007. Crustal attenuation within the Turkish Plateau and surrounding regions[J]. Bull Seismol Soc Am,97(1B):151–161. doi: 10.1785/0120050227