Reliability tests of shear wave velocity structure from joint inversion of multiple types of seismic data
-
摘要: 本文模拟使用青藏高原东南缘区域台网及国家台网的170个宽频台站基于背景噪声、天然地震面波、P波接收函数反演时的实际数据，对青藏高原东南缘假定的初始模型进行恢复，通过计算初始模型台站下方纯路径频散、提取各台站对间的瑞雷波频散曲线、计算理论接收函数以及反演剪切波速度结构来测试使用不同单项数据与联合使用多种数据反演对初始模型的恢复程度。结果表明，同时使用接收函数、基于噪声经验格林函数的群速度、相速度频散以及基于天然地震面波的相速度频散联合反演的剪切波速度结构，充分利用了几种数据的分辨率优势，清晰地分辨出中下地壳及上地幔顶部的低速层。此外，本文也分析了实际数据处理中出现的计算误差、随机噪声干扰对计算结果稳定性的影响。结果显示：对于面波频散，加入1%的误差后，联合反演的结果仍可很好地反映低速层的形态，但是当误差提升至5%后，对最终结果则产生了一定程度的影响；而在接收函数中加入4%的随机噪声时，虽然地幔低速层的上界面和下界面会略微受到随机噪声的影响，但是低速层的深度范围和速度值均得到了较好的恢复。Abstract: Based on the real data from joint inversion of ambient noise, surface wave data, and P wave receiver functions of 170 broad-band seismic stations of national and regional networks of the southeastern margin of Tibetan Plateau and its adjacent areas, we preformed the recovering tests to the presumed initial model of southeastern margin of Tibetan Plateau. We calculated pure path dispersion curves on the basis of the initial model, then retrieved the Rayleigh wave dispersion curves between station pairs and receiver functions beneath each station. Finally, the recovering tests were taken to measure the recovery ability to the initial model based on different seismic data alone and joint inversion of multiple types of seismic data. Our result reveals that joint inversion of receiver function, dispersion cures based on empirical Green’s functions （EGFs） of ambient noise and on teleseismic surface wave data can take full advantage of the resolution of each seismic data, and can resolve the crustal and upper mantle LVZs perfectly. Additionally, we analyzed the resulting reliability on the condition of adding calculation error or random noise. From these tests, surface wave dispersions with 1% error in joint inversion can resolve the low velocity zone （LVZ） commendably, while with 5% error can cause some differences from the initial model. Though receiver functions with 4% random noise in joint inversion can reduce the resolution of the upper and lower boundaries of the mantle LVZ, they can commendably recover the LVZs in terms of occurrence depth and velocities.
-
引言
地震是地球上频繁发生的一种自然灾害，是不可避免的，给人类社会带来了巨大的损失（杨桐等，2013）。近年来的大量研究表明，电磁异常对地震十分敏感（赵国泽等，2017）。国内外学者研究发现电离层观测到的多种数据，例如电场、磁场、等离子体密度、空间粒子通量变化等，在地震前后均出现异常变化，这些电磁异常现象与地震活动都具有相关性（安张辉等，2011；万剑华等，2012；泽仁志玛等，2012；Zhang et al，2012；闫相相等，2014）。
自2004年第一颗专门用于探测地震电离层扰动的电磁卫星（Detection of Electro-Magnetic Emissions Transmitted from Earthquake Regions，缩写为DEMETER）发射以来，科学工作者们在地震电离层方面作了大量研究。由于DEMETER是太阳同步轨道卫星，卫星飞过地面不同区域的地方时是固定的，昼侧和夜侧的轨道数据对应的地方时分别为10:30和22:30左右（杨牧萍等，2018），因此对DEMETER卫星电磁数据的背景场研究主要分昼侧和夜侧。杨牧萍等（2018）利用五年的DEMETER电场功率谱数据，分昼侧和夜侧构建了东北亚地区电离层电场背景场并分析其动态变化。Bertello等（2018）分析了八年的DEMETER卫星电场和磁场数据，构建了拉奎拉地震区昼侧和夜侧的电磁背景场，并针对2009年拉奎拉地震提取异常轨道。
蜂群（Swarm）卫星一天约飞行16个圆轨，每个圆轨的飞行时间大约为1.5个小时。其飞过同一区域的地方时是不断变化的，而不同地方时下磁场数据的背景值相差几十甚至上百倍，因此对Swarm卫星建立昼侧、夜侧背景场并不适用，还会影响提取的异常轨道的准确性。目前针对Swarm卫星磁场数据的研究主要是利用滑动时窗等方法进行单轨道分析，例如：Marchetti和Akhoondzadeh （2018）利用Swarm卫星磁场数据对墨西哥地震进行异常提取，在震前130天出现了异常累计加速；Akhoondzadeh等（2018）利用Swarm卫星磁场数据对2017年伊朗地震进行分析，发现在震前8—11天出现了地震前兆异常，然而单轨道分析提取出的异常也有可能由区域内其它短期局部性非震因素引起。此外，卫星观测的磁场数据受地磁活动的影响，当地磁活动很强时，电离层的扰动可能会掩盖因地震引起的磁场变化（李凯艳，2020）。目前排除地磁活动影响的常规作法是选择低地磁活动日的卫星轨道、舍弃高地磁活动日的卫星轨道（泽仁志玛等，2012；Zhang et al，2012；杨牧萍等，2018；Marchetti et al，2019；Ouyang et al，2020），但这样会造成卫星观测数据应用不充分。
考虑到上述不足，本文拟利用Swarm卫星磁场数据建立长期的高时间分辨率的时变背景场，去除地方时对磁场数据的影响，以及区域性的短期非震因素对结果的影响。此外，通过变分模态分解（variational mode decomposition，缩写为VMD）去除地磁活动对数据的影响，使卫星全部的观测数据可以应用于地震分析中。利用建立的时变背景场对2020年MW7.7牙买加地震影响区域内的磁场Y分量（东向分量）数据进行异常提取，并与传统的昼、夜侧背景场的异常提取结果进行对比分析。此外，对岩石层、大气层以及电离层的多参数进行震前异常分析，并对三个圈层的异常出现时间进行解释。
1. 地震与数据
1.1 牙买加地震
2020年1月28日19:10:24，在古巴南部和牙买加西北部加勒比海地区发生了MW7.7强震（USGS，2020a），震中位于（19.419°N，78.756°W）（图1），震源深度为14.9 km。根据Dobrovolsky公式R＝100.43M计算地震影响区域的半径R （Dobrovolsky et al，1979），其中，M为震级，R为孕震区半径。针对2020年牙买加MW7.7地震计算得到R＝2046.4 km，对应于图1所示区域（0.58°S—39.42°N，58.76°W—98.76°W）。
![]()
图 1 2020年牙买加 MW7.7地震的地理位置和地震影响区（USGS，2020a）Figure 1. Location of the 2020 Jamaica MW7.7 earthquake and related earthquake affected areas （USGS，2020a）1.2 电磁卫星数据
2013年11月22日，欧洲空间局发射了Swarm卫星群（Friis-Christensen et al，2006），星群包括Alpha，Bravo和Charlie三颗卫星。Swarm卫星群能实现对地球磁场高精度、高分辨率的测量，目前已经广泛应用到地震前兆异常提取中。由于中低纬地区地震磁场Y分量受外源磁场扰动小（在高纬场向电流等也会影响Y分量），更有可能受到岩石圈活动的影响（Pulinets，Ouzounov，2011），因此，本文选择Swarm卫星群Alpha卫星的磁场Y分量数据进行研究，数据采样率为1 Hz。考虑到Alpha卫星完成一个轨道大概需要90 min，地方时变化一个小时需要大概11天，大约133天可以完成地方时的重访（李凯艳，2020），本文选用2017年1月1日至2020年6月30日的磁场Y分量数据构建牙买加地震影响区域内的时变背景场，将该区域内震前90天至地震当天的轨道数据作为异常提取对象。
1.3 太阳与地磁活动指数
太阳是地磁场变化的根本来源（丁鉴海等，2011），太阳射电流量F10.7能表征太阳活动性，一般F10.7低于90 sfu的年份是太阳活动低年（李凯艳，2020）。
卫星观测数据受地磁活动的影响，而地磁指数是对空间电流体系引起的地磁扰动的定量描述，可以表征地磁活动水平（王楚钦等，2015；金巍等，2017）。本文选用的地磁活动指数为ap指数，该指数又被称为“等效的行星性3小时幅度”，是描述全球3小时时段内地磁活动的指标（丁鉴海等，2011）。
2. 时变背景场建立方法
时变背景场的建立分为两个步骤：首先利用太阳活动指数F10.7排除太阳活动性高的数据，并基于VMD将地磁活动对磁场数据的影响去除；然后利用区域网格化建立时间分辨率为1 h的时变背景场，去除地方时对磁场数据的影响。
利用VMD去除数据的部分分量，降低地磁活动对卫星磁场数据的影响。VMD是Dragomiretskiy和Zosso （2014）提出的一种自适应算法，该方法可以自适应地将地磁活动信号从卫星磁场数据中有效分离，这样既可以降低地磁活动这一非震因素的干扰，又可以利用更多的卫星磁场数据。
VMD是通过迭代搜寻约束变分模型的最优解，自适应地将信号分解成多个有限带宽的固有模态函数，每个模态函数的中心频率及带宽在迭代求解变分模型的过程中不断更新（Dragomiretskiy，Zosso，2014；池成全，2020；赵亚军等，2020），其约束变分模型表达式为
$$ \left\{ \begin{array}{l} \mathop {\min }\limits_{\left\{ {\mathop u\nolimits_k } \right\}\left\{ {\mathop \omega \nolimits_k } \right\}} \left\{ {\displaystyle\sum\limits_k {\mathop {\left\| {\mathop {\text{∂}} \nolimits_t \left[\left({\delta} （ t ） ＋ \frac{j}{{{\text{π }}t}}\right) * \mathop u\nolimits_k （ t ） \right]\mathop {\text{e}}\nolimits^{ - {\rm{j}}\mathop \omega \nolimits_k t} } \right\|}\nolimits_2^2 } } \right\} \\ s.t.\displaystyle\sum\limits_k {\mathop u\nolimits_k } ＝ f \qquad k {\text{≥}}2， \end{array} \right. \text{，} $$ (1) 式中，{uk}＝{u1， u2， ···， uk}为分解后的模态，{ωk}＝{ω1， ω2， ···， ωk}为各模态对应的中心频率，∂t为随时间的导数，δ（t）为狄拉克方程，$* $为卷积运算，|| ||2为向量的二范数，f为分解出的k个模态之和。
选用2017年1月1日至2020年6月30日期间Swarm星群Alpha卫星的磁场Y分量数据，对数据作预处理，减去CHAOS-7模型（Finlay et al，2020）以去除主磁场、岩石圈磁场以及外源磁场。对预处理后的每一条轨道数据进行VMD，分解成两个模态。通过计算分解出模态的平均能量与ap指数的相关系数，定量分析两个模态与地磁指数的关系。两个模态的平均能量φIMFi的计算公式为
$$ {\mathop \varphi \nolimits_{ {{\rm{IMF}}i} } ＝ \frac{{\displaystyle\sum\limits_{j ＝ 1}^N {\mathop {\left| { {{\rm{IMF}}i} （ j ） } \right|}\nolimits^2 } }}{N}\qquad i ＝ 1， 2 }， $$ (2) 式中，i代表第i个模态，j代表第j个数据点，N代表该轨道的数据点数。
计算每一条轨道分解出的两个模态的平均能量φIMFi与ap指数的相关系数，结果如图2所示。图2a中蓝线和红线分别是每条轨道分解出的第一个模态的平均能量φIMF1和对应的滑动平均值$\overline \varphi $IMF1；图2b中蓝线和红线分别是每条轨道分解出的第二个模态的平均能量φIMF2和对应的滑动平均值$\overline \varphi $IMF2；图2c中蓝线和红线分别是每条轨道对应的ap指数和对应的滑动平均值$\overline a $p。
![]()
图 2 两个模态的平均能量与对应的$a_p $指数（a） 轨道第一个模态的能量；（b） 轨道第二个模态的能量；（c） 轨道对应的$a_p $指数Figure 2. Results of average energy of the two modes and their corresponding $a_p $ index（a） The average energy of the first mode of all the tracks；（b） The average energy of the second mode of all the tracks； （c） The $a_p $ index at the moment corresponding to the time of the track由图2可知，在$\overline a $p取得高值时，$\overline \varphi $IMF1也取得高值，在$\overline a $p取得低值时，$\overline \varphi $IMF1也取得低值，二者的变化趋势很接近，而$\overline \varphi $IMF2与$\overline a $p的变化趋势没有明显的相似之处。为了进一步判断$\overline \varphi $IMF1，$\overline \varphi $IMF2以及$\overline a $p之间的关系，利用滑动窗分别计算$\overline \varphi $IMF1与$\overline a $p的相关系数λ1和$\overline \varphi $IMF2与$\overline a $p的相关系数λ2，分别计算相关系数的分布，结果如图3所示，图中红点是$\overline \varphi $IMF1与$\overline a $p的相关系数，占比大的值集中在0.7—0.9之间；图中蓝点是$\overline \varphi $IMF2与$\overline a $p的相关系数，较均匀地分布在0.1—0.6之间。其中$\overline \varphi $IMF1与$\overline a $p的相关系数中值为0.800 3，而$\overline \varphi $IMF2与$\overline a $p的相关系数中值为0.305 1。
![]()
图 3 模态的平均能量与对应的$\overline a_p $的相关系数Figure 3. The correlation coefficients between $\overline \varphi_{\rm{IMF1}} $ and $\overline a_p $ and those between $\overline \varphi_{\rm{IMF2}} $ and $\overline a_p $$\overline \varphi $IMF1与$\overline a $p的相关性较高，而$\overline \varphi $IMF2与$\overline a $p相关性不高，故第一个模态受到地磁活动的影响较大，第二个模态受到地磁活动的影响较小。因此，剔除第一个模态以降低地磁活动对磁场数据的影响，接下来只对去除第一个模态后的轨道数据进行分析。
对2020年牙买加地震影响区域进行网格化，划分成5°×5° （地理经纬度）的64个小网格，去除该区域内M6.0以上地震（USGS，2020b）分别在其Dobrovolsky范围内（Dobrovolsky et al，1979）震前90天的数据，以排除地震对背景场的影响；接着把剩余的轨道数据按照地方时（local time，缩写为LT）划分成24个地方时范围，以1 h为间隔，分别是LT0 （00:00—00:59），LT1 （01:00—01:59），···，LT22 （22:00—22:59），LT23 （23:00—23:59）。根据地磁季节划分的方法，把3—4月份划分为春季，5—8月份划分为夏季，9—10月份划分为秋季，11—12月份以及来年1—2月划分为冬季（杨牧萍等，2018）。由于2020年牙买加地震震前90天正处于冬季，因此，本文建立的是冬季时变背景场。
计算每个小网格内24个地方时的数据点的能量，计算方法为
$$ \varepsilon （ t ） ＝ \mathop {\left| {f （ t ） － {{\rm{IMF1}}} （ t ） } \right|}\nolimits^2 \text{，} $$ (3) 式中，f （t）是预处理后的数据，IMF1（t）是f （t）分解出的第一个模态。
分别以网格内所有数据点能量的中值作为对应地方时的背景值来建立冬季时变背景场，如图4所示（图中白色区域代表该区域无数据）。从图中可以看出，不同地方时的背景值相差很大，其中，地方时为白天（0—5时，18时—23时）的背景值远大于地方时为夜晚（6时—17时）的背景值。
![]()
图 4 牙买加地震影响区域在不同地方时（LT）的时变背景场Figure 4. The time-varying background field at different local time inside the earthquake affected areas of the Jamaica earthquake选取整个牙买加地震影响区域，分别计算时变背景场里整个研究区域每个地方时范围内的中值并将其作为区域背景值，进一步分析该区域时变背景场在不同地方时的背景值差异。将该区域24个地方时的区域背景值进行拟合，结果如图5所示。当地方时为10时，背景值最大，为1.059 nT2；当地方时为4时，背景值最小，为0.005 9 nT2，二者相差180倍。由此可见，地方时对Swarm卫星磁场数据的影响非常大，建立时变背景场去除地方时的影响是非常有必要的，这会影响后续异常提取的准确性。
![]()
图 5 牙买加地震影响区域内24个地方时的平均背景值Figure 5. The average background values of 24 local times inside the earthquake-affected area of the Jamaica earthquake3. 震前卫星磁场异常提取结果与讨论
3.1 基于时变背景场与昼侧、夜侧背景场的异常提取结果
本文通过相同的方法建立昼侧、夜侧背景场，选取2017年1月1日至2020年6月30日期间牙买加地震影响区域内的磁场数据，首先将数据划分到四季，然后将冬季数据按照地方时分成昼侧和夜侧的数据集，即地方时为0—5时，18时—23时期间的数据划分到夜侧数据集；地方时为6时—17时期间的数据划分到昼侧数据集，利用区域网格化和中值建立冬季昼侧、夜侧背景场。
分别使用时变背景场与昼侧、夜侧背景场提取牙买加地震震前异常轨道。其中基于时变背景场提取异常轨道是以每个网格内对应地方时的时变背景场为背景值，利用滑动四分位（Liu et al，2000）计算网格内对应地方时的所有数据点能量的上四分位Q1、下四分位Q3和四分距Q1—Q3，以2.5倍的四分距作为阈值偏差。滑动四分位的计算方法如式（4）—（6）所示。
将一个网格内一个地方时的所有数据点的能量进行排序，结果为
$$ \varepsilon （ \mathop t\nolimits_1 ） {\text{≤}} \varepsilon （ \mathop t\nolimits_2 ） {\text{≤}} \cdots {\text{≤}} \varepsilon （ \mathop t\nolimits_k ） \cdots {\text{≤}} \varepsilon （ \mathop t\nolimits_n ） \qquad k ＝ 1， 2， \cdots ， n \;\;\text{，} $$ (4) $$ {{\displaystyle Q}}_{3}＝\left\{\begin{array}{l} \dfrac{\varepsilon （ {{\displaystyle t}}_{k} ） ＋\varepsilon （ {{\displaystyle t}}_{k＋1} ） }{2}\qquad n＝4k\; 或 \; n＝4k＋1，\\ \varepsilon （ {{\displaystyle t}}_{k+1} ） \qquad\;\;\quad\qquad n＝4k＋2\; 或 \;n＝4k＋3，\end{array} \right. $$ (5) $$ \mathop Q\nolimits_{\text{1}} ＝ \left\{ \begin{array}{ll} \dfrac{{\varepsilon （ \mathop t\nolimits_{3k} ） ＋ \varepsilon （ \mathop t\nolimits_{3k ＋ 1} ） }}{2}& \quad n ＝ 4k ，\\ \dfrac{{\varepsilon （ \mathop t\nolimits_{3k + 1} ） ＋ \varepsilon （ \mathop t\nolimits_{3k ＋ 2} ） }}{2}& \quad n ＝ 4k ＋ 1 ，\\ \varepsilon （ \mathop t\nolimits_{3k ＋ 2} ） & \quad n ＝ 4k ＋ 2 ，\\ \varepsilon （ \mathop t\nolimits_{3k ＋ 3} ） & \quad n ＝ 4k ＋ 3 {\text{.}} \end{array} \right. $$ (6) 当一条轨道经过某一网格的数据能量的中值与对应地方时的背景值偏差超过2.5倍四分距时，即认为这条轨道为异常轨道。基于昼侧、夜侧背景场提取异常轨道与基于时变背景场提取异常轨道的方法相同。
分别将基于时变背景场与昼侧、夜侧背景场提取出的异常轨道进行累计，当第k天提取出n条异常轨道时，异常轨道累计数目N（k）在N（k－1）的基础上加n。异常轨道数目累计结果如图6所示。由图6a可知：震前90天至震前50天，基于时变背景场提取出的异常轨道累计数目呈线性增长；震前50天至震前43天，异常轨道累计数目加速增长；震前43天至地震当天，缓慢恢复到线性增长。由图6b可见：震前90天至震前54天，基于昼侧、夜侧背景场提取出的异常轨道累计数目呈现线性增长；震前54天至震前43天，异常轨道累计数目加速增长；震前43天至地震当天，缓慢恢复至基本不变化。
![]()
图 6 基于时变背景场（a）与昼侧、夜侧背景场（b）分别提取的异常轨道累计结果Figure 6. The cumulative number of anomalous tracks over time based on the time-varying background field （a） and the day-side and the night-side background field （b）3.2 异常提取结果对比与分析
为了说明时变背景场相较昼侧、夜侧背景场的优势，本文将从背景场和震前异常提取两方面进行比较分析。考虑到地方时为白天的昼侧背景场与地方时为夜晚的夜侧背景场差异很大，将研究区域内24个地方时的时变背景场分成昼侧和夜侧并分别与昼侧、夜侧背景场进行对比，如图7所示。可见：对于地方时为白天的时变背景场，地方时为10时的背景值最大，为1.059 0 nT2，地方时为17时的背景值最小，为0.030 5 nT2，二者相差35倍；对于地方时为夜晚的时变背景场，地方时为20时的背景值最大，为0.018 5 nT2，地方时为0时的背景值最小，为0.005 9 nT2，二者相差3倍；而该区域内昼侧背景值为0.273 7 nT2，夜侧背景值为0.010 3 nT2。由此可见，不管是白天还是夜晚，在不同的地方时，时变背景场的背景值相差均很大，有的会高达几十倍。因此，建立时变背景场可以有效地把不同地方时的背景值差异凸显出来，这对于后续准确提取异常轨道是非常重要的。
![]()
图 7 地方时为白天（a）和夜晚（b）时昼侧、夜侧的两个背景场对比Figure 7. Comparison between dayside and night side background field in the day time （a） and the night time （b）为了分析基于时变背景场与昼侧、夜侧背景场提取出异常轨道差异的原因，对两种背景场提取异常的阈值作对比分析，如图8a所示，以一个小网格为例。图8b是两个背景场阈值差异对比，若一条轨道在某一网格的能量值落在区域A，则该轨道在时变背景场与昼侧、夜侧背景场都被判定为非异常的轨道；若落在区域B，则该轨道是昼侧、夜侧背景场提取异常少于时变背景场的轨道，集中在白天和夜晚低背景值的地方时范围内；若落在区域C，则该轨道是昼侧、夜侧背景场提取异常多于时变背景场的轨道，集中在白天和夜晚高背景值的地方时范围内。
![]()
图 8 时变背景场与昼侧、夜侧背景场提取异常的阈值（a）及其对比图（b）Figure 8. The threshold （a） of the time-varying background field and the day-side and the night-side background field used to extract anomalous tracks and its comparison （b）昼侧、夜侧背景场比时变背景场多提取出的异常轨道都是集中在背景值高的地方时，这是由于昼侧、夜侧背景场的时间分辨率比时变背景场低，不管是白天还是夜晚，高值背景值会被低值背景值拉低，导致部分不是异常的轨道被错误地识别为异常轨道。同理，昼侧、夜侧背景场比时变背景场少提取出的异常轨道均集中在背景值低的地方时，低值背景值会被高值背景值拉高，导致部分异常轨道不能被识别。因此，建立高时间分辨率的时变背景场相较昼侧、夜侧背景场能更加准确地识别并提取异常轨道。
3.3 震前多圈层异常提取结果与分析
本文分析了震前岩石层、大气层和电离层的参量，对地震影响范围内M＞4.6浅源地震（深度＜50 km）进行震前能量累计（Kanamori，Anderson，1975；Mignan et al，2007；Marchetti et al，2020）。图9a是岩石层的震前地震能量累计结果，对岩石层应变采用了贝尼奥夫应变分析方法累计地震前能量，从图中可以看出，在震前35天出现加速现象；图9b是大气层的总水气柱异常累计结果，大气层的总水气柱（total column water vapour，缩写为TCWV）异常提取对象为震中3°×3°范围内震前90天的TCWV数据，根据地震前兆识别的气候分析算法（climatological analysis for seismic precursor identification，缩写为CAPRI）提取大气参数地震前兆异常（Piscini et al，2019），图中显示震前45天出现异常累计加速；图9c是电离层的磁场异常轨道累计结果，震前50天出现异常轨道累计加速
![]()
图 9 岩石层能量（a）、大气层总水气柱（b）和电离层磁场异常轨道（c）累计结果Figure 9. Accumulative number of anomalies of the lithospheric energy （a）， the total column water vapourof atmosphere （b） and the ionosphere magnetic field （c）自震前78天，岩石层能量开始缓慢累计，震前50天，在与牙买加地震相同的断层上发生MW5.0地震（2019年12月9日，19.08°N，80.44°W），这可能是主震前的微破裂，这种微破裂可能产生超低频（ultra-low frequency，缩写为ULF）电磁波（Molchanov，Hayakawa，1995）。由于电磁波在电离层的扰动，我们推测震前50天观测到电离层磁场数据的异常累计加速现象可能与电磁波的扰动有关。Ventura和Di Giovambattista （2013）发现地震前可能存在流体（或气体）的地下运动，我们推测流体的地下运动可能与微破裂同时发生，可能需要几天的时间从断层向上传播到达海洋底部，并在水中释放一些气体（Ouzounov et al，2018），这些气体可能会导致大气参数的异常，因此，震前45天观测到总水气柱异常累计加速现象可能是由气体的释放导致的。马瑾和郭彦双（2014）提出，断层带上存在相对弱和相对强的部位，前者往往表现为断层预滑、慢地震或弱震，后者是应力锁定部分和快速失稳区域（Noda et al，2013）；在断层的黏滑过程中实际存在两次失稳，前者与弱部位的释放有关，后者与强部位的快速释放有关，表现为强震。因此本文推测岩石层能量累积在震前35至震前34天的加速可能与弱部位的释放有关，在震前22天至震前20天的加速可能与强部位的快速释放有关。岩石层、大气层和电离层中不同参数的异常累积结果表明，这些异常可能与牙买加地震的孕育有关。
4. 讨论与结论
本文针对蜂群卫星地方时不断变化的特点，建立了时间分辨率高的卫星磁场数据时变背景场，有效去除了地方时对数据的影响。基于建立的时变背景场和传统的昼、夜侧背景场，分别对2020年牙买加MW7.7地震进行了震前异常提取，两者的对比分析结果显示：利用时变背景场可以准确地提取出异常轨道，而昼、夜侧背景场由于时间分辨率低、背景值不准确，会导致异常轨道的提取结果不准确。因此，本文提出的时变背景场相较传统的昼、夜侧背景场能更准确地提取出异常轨道，提高了地震震前异常提取的准确度，为地震前兆异常的研究提供了新思路。
通过时变背景场提取的异常轨道累计结果在震前50至43天出现加速增长，联合岩石层和大气层的研究，牙买加地震先在电离层出现异常，然后在大气层出现异常，最后在岩石层出现异常，结合Molchanov和Hayakawa （1995）、Ventura和Di Giovambattista （2013）、Ouzounov等（2018）以及马瑾和郭彦双（2014）的研究成果，我们对牙买加地震震前岩石层、大气层、电离层参数的异常出现时间的机理提出了合理的猜测：电离层的异常可能是由主震前微破裂产生的ULF电磁波导致的，而大气层的异常可能是由于流体的地下运动产生的，岩石层的异常则可能是由断层的黏滑过程中存在的两次失稳导致的。
本文仅对牙买加地震进行了分析，为了获取更加令人信服的结果，需要进一步对多个震例进行统计分析，特别是M7.0以上的大地震。此外，虽然本文联合了岩石层、大气层和电离层的异常提取结果进行了简单的分析，但是如何验证提取的异常与震前活动是否相关尚需进一步的探索和研究。
-
![]()
图 1 青藏高原东南缘台站（a）及剪切波速度分布（b）示意图
Figure 1. Images of station locations （a） and shear wave velocity （b） in the southeastern margin of the Tibetan Plateau
![]()
图 2 成像分辨率检测试验中各周期实际使用的频散数目
红色圆点与绿色圆点分别表示基于背景噪声经验格林函数和基于天然地震面波的频散数目；蓝色圆点表示在重叠周期21—60 s内基于两种数据最终提取的频散数目
Figure 2. Number of dispersion measurements used in resolution test at each periods
The red dots and green dots represent the number of dispersion measurements retrieved from empirical Green’s functions based on ambient noise，two-station analysis using teleseismic sur-face wave data，respectively；blue dots represent the number of final dispersion measurements retrieved from the two types of data in the overlapping periods （21−60 s）
![]()
图 3 格点（102°E，25°N）下方的联合反演结果
（a） 初始模型和由面波频散曲线及接收函数联合反演的输出模型；（b） 根据图（a）中模型计算的P波接收函数；（c） 基于初始模型的瑞雷波相速度、群速度频散曲线
Figure 3. Joint inversion results beneath the grid （102°E，25°N）
（a） The initial model and output model after joint inversion of surface wave dispersion measurements and receiver function； （b） The P wave receiver functions calculated according to the initial models in Fig.（a）； （c） Rayleigh wave phase velocity and group velocity dispersion measurements calculated according to the initial model respectively
![]()
图 4 由图1初始模型计算得到的不同周期T瑞雷波的相速度、群速度图像
Figure 4. Rayleigh wave phase velocity and group velocity images at different periods T calculated from the initial model shown in Fig.1
![]()
图 5 使用不同数据反演所得沿25°N剖面的剪切波速度结构
（a） 仅使用基于背景噪声的瑞雷波群速度频散反演；（b） 仅使用基于天然地震面波的瑞雷波相速度频散反演；（c） 同时使用基于背景噪声及天然地震面波的瑞雷波群速度、相速度频散反演；（d） 仅使用接收函数反演；（e） 同时使用3种数据联合反演
Figure 5. Cross sections of shear wave velocity structure along 25°N from inversion by using different seismic data
（a） The output model inverted from Rayleigh wave group velocity dispersions based on ambient noise；（b） The output model inverted from Rayleigh wave phase velocity dispersions based on teleseismic wave；（c） The output model jointly inverted from Rayleigh wave group velocity and phase velocity dispersions based on ambient noise and teleseismic wave；（d） The output model inverted from receiver function；（e） The output model jointly inverted from three types of seismic wave data
![]()
图 6 沿25°N剖面使用不同数据反演后模型与初始模型的速度扰动分布
（a） 仅使用基于背景噪声的瑞雷波群速度频散反演；（b） 仅使用基于天然地震面波的瑞雷波相速度频散反演；（c） 同时使用基于背景噪声及天然地震面波的瑞雷波群速度、相速度频散反演；（d） 同时使用3种数据联合反演
Figure 6. Cross sections of shear wave velocity perturbation along 25°N using different types of seismic data
（a） Cross sections of shear wave velocity perturbation calculated based on Rayleigh wave group velocity dispersion measurements；（b） Similar to Fig.（a） but calculated based on Rayleigh wave phase velocity dispersion measurements from teleseismic wave；（c） Similar to Fig.（a） but calculated based on both Rayleigh wave group velocity and phase velocity dispersion measurements from ambient noise and tele-seismic wave；（d） Similar to Fig.（a） but calculated based on three types of seismic data
![]()
图 7 初始模型纯路径频散加入1%初始误差后反演的瑞雷波相速度、群速度图像
Figure 7. Rayleigh wave phase velocity and group velocity images at different periods by inversion with 1% initial error added to the pure path dispersions based on the the initial model
![]()
图 8 初始模型纯路径频散加入1% （a）和5% （b）初始误差后反演所得到的剪切波速度结构
Figure 8. Cross sections of shear wave velocity along 25°N inverted using pure path dispersion measurements with 1% （a） and 5% （b） initial error added to the pure path dispersions based on the the initial model
![]()
图 9 格点（102°E，25°N）下方100组基于加入随机噪声的接收函数联合反演的结果
（a） 初始模型和联合反演输出模型；（b） 联合反演使用的P波接收函数；（c） 瑞雷波相速度、群速度频散曲线
Figure 9. The 100 joint inversion results beneath the grid （102°E， 25°N） based on receiver functions with different random noise
（a） The initial model and the output models based on joint inversion；（b） The P wave receiver functions used in joint inversions；（c） The Rayleigh wave phase velocity and group velocity dispersion measurements used in joint inversions
-
胡家富,朱雄关,夏静瑜,陈赟. 2005. 利用面波和接收函数联合反演滇西地区壳幔速度结构[J]. 地球物理学报,48(5):1069–1076. doi: 10.3321/j.issn:0001-5733.2005.05.013 Hu J F,Zhu X G,Xia J Y,Chen Y. 2005. Using surface wave and receiver function to jointly inverse the crust-mantle velocity structure in the west Yunnan area[J]. Chinese Journal of Geophysics,48(5):1069–1076 (in Chinese). doi: 10.1002/cjg2.750
刘启元,李昱,陈九辉,van der Hilst R D,郭飚,王峻,齐少华,李顺成. 2010. 基于贝叶斯理论的接收函数与环境噪声联合反演[J]. 地球物理学报,53(11):2603–2612. Liu Q Y,Li Y,Chen J H,van der Hilst R D,Guo B,Wang J,Qi S H,Li S C. 2010. Joint inversion of receiver function and ambient noise based on Bayesian theory[J]. Chinese Journal of Geophysics,53(11):2603–2612 (in Chinese).
Bao X W,Sun X X,Xu M J,Eaton D W,Song X D,Wang L S,Ding Z F,Mi N,Li H,Yu D Y,Huang Z C,Wang P. 2015. Two crustal low-velocity channels beneath SE Tibet revealed by joint inversion of Rayleigh wave dispersion and receiver functions[J]. Earth Planet Sci Lett,415:16–24.
Bensen G B,Ritzwoller M H,Barmin M P,Levshin A L,Lin F,Moschetti M P,Shapiro N M,Yang Y. 2007. Processing seismic ambient noise data to obtain reliable broad-band surface wave dispersion measurements[J]. Geophys J Int,169(3):1239–1260. doi: 10.1111/gji.2007.169.issue-3
Bensen G B,Ritzwoller M H,Shapiro N M. 2008. Broadband ambient noise surface wave tomography across the United States[J]. J Geophys Res,113(B5):B5306. doi: 10.1029/2007JB005248
Bodin T,Sambridge M,TkalčIć H,Arroucau P,Gallagher K,Rawlinson N. 2012. Transdimensional inversion of receiver functions and surface wave dispersion[J]. J Geophys Res,117(B2):B02301.
Campillo M,Paul A. 2003. Long-range correlations in the diffuse seismic coda[J]. Science,299(5606):547–549. doi: 10.1126/science.1078551
Chang S J,Baag C E,Langston C A. 2004. Joint analysis of teleseismic receiver functions and surface wave dispersion using the genetic algorithm[J]. Bull Seismol Soc Am,94(2):691–704. doi: 10.1785/0120030110
Fang H J,Zhang H J,Yao H J,Allam A,Zigone D,Ben-Zion Y,Thurber C,van der Hilst R D. 2016. A new algorithm for three-dimensional joint inversion of body wave and surface wave data and its application to the Southern California Plate boundary region[J]. J Geophys Res,121(5):3557–3569. doi: 10.1002/2015JB012702
Julià J,Ammon C J,Herrmann R B,Correig A M. 2000. Joint inversion of receiver function and surface wave dispersion observations[J]. Geophys J Int,143(1):99–112. doi: 10.1046/j.1365-246x.2000.00217.x
Kang D,Shen W S,Ning J Y,Ritzwoller M H. 2016. Seismic evidence for lithospheric modification associated with intracontinental volcanism in northeastern China[J]. Geophys J Int,204(1):215–235.
Lawrence J F,Wiens D A. 2004. Combined receiver-function and surface wave phase-velocity inversion using a niching genetic algorithm:Application to Patagonia[J]. Bull Seismol Soc Am,94(3):977–987. doi: 10.1785/0120030172
Li Y H,Wu Q J,Zhang R Q,Tian X B,Zeng R S. 2008. The crust and upper mantle structure beneath Yunnan from joint inversion of receiver functions and Rayleigh wave dispersion data[J]. Phy Earth Planet Inter,170(1/2):134–146.
Lin F C,Ritzwoller M H,Townend J,Bannister S,Savage M K. 2007. Ambient noise Rayleigh wave tomography of New Zea-land[J]. Geophys J Int,170(2):649–666. doi: 10.1111/gji.2007.170.issue-2
Lin F C,Moschetti M P,Ritzwoller M H. 2008. Surface wave tomography of the western United States from ambient seismic noise:Rayleigh and Love wave phase velocity maps[J]. Geophys J Int,173(1):281–298. doi: 10.1111/gji.2008.173.issue-1
Liu Q Y,van der Hilst R D,Li Y,Yao H J,Chen J H,Guo B,Qi S H,Wang J,Huang H,Li S C. 2014. Eastward expansion of the Tibetan Plateau by crustal flow and strain partitioning across faults[J]. Nat Geosci,7(5):361–365. doi: 10.1038/ngeo2130
Love A E H. 1911. Some Problems of Geodynamics[M]. New York: Cambridge University Press: 1–210.
Moschetti M P,Ritzwoller M H,Lin F C,Yang Y. 2010. Crustal shear wave velocity structure of the western United States inferred from ambient seismic noise and earthquake data[J]. J Geophs Res,115(B10):B10306. doi: 10.1029/2010JB007448
Obrebski M,Allen R M,Zhang F X,Pan J T,Wu Q J,Hung S H. 2012. Shear wave tomography of China using joint inversion of body and surface wave constraints[J]. J Geophys Res,117(B1):B01311.
Paige C C,Saunders M A. 1982a. LSQR:An algorithm for sparse linear equations and sparse least squares[J]. ACM Trans Math Software,8(1):43–71. doi: 10.1145/355984.355989
Paige C C,Saunders M A. 1982b. LSQR:Sparse linear equations and least squares problems[J]. ACM Trans Math Software,8(2):195–209. doi: 10.1145/355993.356000
Press F. 1956. Determination of crustal structure from phase velocity of Rayleigh waves part I:Southern California[J]. GSA Bull,67(12):1647–1658. doi: 10.1130/0016-7606(1956)67[1647:DOCSFP]2.0.CO;2
Saint Louis University. 2013. Computer programs in seismology[CP/OL]. [2018−05−01]. http://www.eas.slu.edu/eqc/eqccps.html.
Shapiro N M,Ritzwoller M H. 2002. Monte-Carlo inversion for a global shear-velocity model of the crust and upper mantle[J]. Geophys J Int,151(1):88–105. doi: 10.1046/j.1365-246X.2002.01742.x
Shapiro N M,Campillo M. 2004. Emergence of broadband Rayleigh waves from correlations of the ambient seismic noise[J]. Geophys Res Lett,31(7):L07614.
Shapiro N M,Campillo M,Stehly L,Ritzwoller M H. 2005. High-resolution surface-wave tomography from ambient seismic noise[J]. Science,307(5715):1615–1618. doi: 10.1126/science.1108339
Stoneley R. 1926. The effect of the ocean on Rayleigh waves[J]. Geophys Suppl Mon Not Roy Astronom Soc,1(7):349–356.
Sun X L,Song X D,Zheng S H,Yang Y J,Michiael H R. 2010. Three dimensional shear wave velocity structure of the crust and upper mantle beneath China from ambient noise surface wave tomography[J]. Earthquake Science,23(5):449–463. doi: 10.1007/s11589-010-0744-4
Sun X X,Bao X W,Xu M J,Eaton D W,Song X D,Wang L S,Ding Z F,Mi N,Yu D Y,Li H. 2014. Crustal structure beneath SE Tibet from joint analysis of receiver functions and Rayleigh wave dispersion[J]. Geophys Res Lett,402(5):1479–1484.
Tokam A P K,Tabod C T,Nyblade A A,Julià J,Wiens D A,Pasyanos M E. 2010. Structure of the crust beneath Cameroon,West Africa,from the joint inversion of Rayleigh wave group velocities and receiver functions[J]. Geophys J Int,183(2):1061–1076. doi: 10.1111/j.1365-246X.2010.04776.x
Wang W L,Wu J P,Fang L H,Lai G J,Yang T,Cai Y. 2014. S wave velocity structure in southwest China from surface wave tomography and receiver functions[J]. J Geophys Res,119(2):1061–1078. doi: 10.1002/2013JB010317
Wessel P,Smith W H F. 1998. New,improved version of generic mapping tools released[J]. Eos Trans AGU,79(47):579. doi: 10.1029/98EO00426
Yang Y J,Ritzwoller M H,Levshin A L,Shapiro N M. 2007. Ambient noise Rayleigh wave tomography across Europe[J]. Geophys J Int,168(1):259–274. doi: 10.1111/gji.2007.168.issue-1
Yang Y J,Zheng Y,Chen J,Zhou S Y,Celyan S,Sandvol E,Tilmann F,Priestley K,Hearn T M,Ni J F,Brown L D,Ritzwoller M H. 2010. Rayleigh wave phase velocity maps of Tibet and the surrounding regions from ambient seismic noise tomography[J]. Geochem Geophys Geosyst,11(8):Q08010.
Yao H J,van der Hilst R D,de Hoop M V. 2006. Surface-wave array tomography in SE Tibet from ambient seismic noise and two-station analysis: I . Phase velocity maps[J]. Geophys J Int,166(2):732–744. doi: 10.1111/gji.2006.166.issue-2
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: Ⅱ . Crust and upper-mantle structure[J]. Geophys J Int,173(1):205–219. doi: 10.1111/gji.2008.173.issue-1
Zhang P,Yao H J. 2017. Stepwise joint inversion of surface wave dispersion,Rayleigh wave ZH ratio,and receiver function data for 1D crustal shear wave velocity structure[J]. Earthquake Science,30(5/6):229–238. doi: 10.1007/s11589-017-0197-0
Zheng S H,Sun X L,Song X D,Yang Y J,Ritzwoller M H. 2008. Surface wave tomography of China from ambient seismic noise correlation[J]. Geochem Geophys Geosyst,9(5):Q05020.
Zheng X,Zhao C P,Zhou L Q,Zheng S H. 2019. Crustal and upper mantle structure beneath SE Tibetan Plateau from joint inversion of multiple types of seismic data[J]. Geophys J Int,217:331–345. doi: 10.1093/gji/ggz027
Zhou L Q,Xie J Y,Shen W S,Zheng Y,Yang Y J,Shi H X,Ritzwoller M H. 2012. The structure of the crust and uppermost mantle beneath South China from ambient noise and earthquake tomography[J]. Geophys J Int,189(3):1565–1583. doi: 10.1111/gji.2012.189.issue-3
下载:
计量
- 文章访问数: 1894
- HTML全文浏览量: 1182
- PDF下载量: 83
