2019, Vol. 62
2. 南方科技大学地球与空间科学系, 广东深圳 518055
2. Department of Earth and Space Sciences, Southern University of Science and Technology, ShenzhenGuangdong 518055, China
背景噪声信号是由许多随机产生的不同类型的源激发的包括体波、面波、衍射波以及散射波组成的复杂波场.背景噪声信号中大部分能量是由面波构成的(Aki, 1957; Okada, 2006)，而面波的频散特性是与地下速度结构密切相关的，因而我们常常利用背景噪声的面波信息研究地下的速度结构.而面波勘探的重要一步是提取相速度的频散曲线.提取相速度频散曲线的方法主要有：1.空间自相关法(SPAC)(Aki, 1957)；2.频率-波数谱法(f-k法)(Lacoss et al., 1969; Capon, 1969)；3.τ-p变换法(McMechan and Yedlin, 1981)；4.多道面波分析方法(MASW)(Park et al., 1998, 2007)；5.倾斜叠加算法(Xia et al., 2007)；6.高分辨率线性拉东变换法(LRT)(罗银河等, 2008, 2009)；7.窄带滤波分析方法(姚华建等, 2004, 2006).上述方法对高阶面波的反应并不灵敏，虽然2—7种方法在理论上可以得到含有高阶的频散曲线，但很难在实际背景噪声中提取出清晰的频散曲线.然而，高阶面波对地下结构的反应较敏感，尤其是对含有软弱层的结构.因此，使用基阶和高阶面波联合反演可以降低反演的非唯一性(Luo et al., 2008)，不仅可以使反演收敛更快(Xia et al., 2000, 2003; 何耀锋, 2005)，还可以提高反演结果的精度.而本文将使用Wang等(2019)提出频率-贝塞尔变换法，使我们能从背景噪声中提取出包含高阶振型的频散曲线.
关东盆地是日本最大都市区所在地.关东盆地是大型的冲积盆地，面积约为17000 km2，最大的沉积厚度约为4 km.关东盆地和它的环绕区域地下结构非常复杂，菲律宾板块在此俯冲到大陆板块之下，而太平洋板块又俯冲到菲律宾板块之下，造成此处为地震多发区域.盆地的沉积层会放大地面运动，并且使地震波的传播过程复杂化.1923年发生在日本都市区域的关东地震，1985年发生在墨西哥城附近的米却肯州地震以及1989年发生在旧金山附近的洛马·普雷塔大地震都证明了地震对人口密集的盆地地区的巨大危害.因而，建立一个精确的关东盆地地下结构模拟强地面运动对此处的危害影响是相当有必要的，Yamanaka和Yamada(2002)曾建立了东京地区浅沉积层的三维速度结构；Japan Seismic Hazard Information Station(J-SHIS)和Headquarters for Earthquake Research Promotion也各自建立并多次更新了关东地区的三维速度模型.Koketsu等(2008, 2009)利用地质填图、地震反折射、重力调查、地震尾波的H/V反演、钻井剖面以及微震调查等多种方法得出全日本综合速度结构模型JIVSM，这是用于评价日本大地震长周期强地面运动使用最为广泛的速度结构模型.本文中，我们将频率-贝塞尔变换法应用于MeSO-net背景噪声数据，反演得到S波速度结构，对关东盆地的速度结构作出评价以及研究方法上的补充.
1 频率-贝塞尔变换法(F-J method)地震噪声互相关函数的频率-贝塞尔变换的定义为(Wang et al., 2019)：
|
(1) |
其中
|
(2) |
其中，A为常数，代表I于Im之间的比例系数.gz是垂直位移场k-积分放入核函数.可以证明：满足频散关系的频散点[ω, cn(ω)](n=0, 1, 2, …)是核函数gz的极点，即当速度c趋近相速度cn(ω)，核函数gz会趋于无穷大(Kennett, 1985; Chen, 1999; Wang et al., 2019).由于实际资料有限，只能近似计算F-J变换的积分(1), 因此公式(2)的关系只能近似成立，I(ω, c)在频散点虽然不会达到无穷大但仍具有突出的极大值.利用I(ω, c)这一特点，我们可以提取频散曲线(Wang et al., 2019).在数值计算中，我们采用如下近似公式计算噪声相关函数
|
(3) |
该公式的详细推导请参考Wang等(2019).
2 反演算法从频率-相速度能量谱或频率-波数能量谱提取频散曲线后，我们需要利用频散曲线反演地下S波速度结构.在本研究中，我们采用BFGS拟牛顿法(何耀锋, 2005; 田玥和陈晓非, 2006; Pan et al., 2019).
为了利用多阶频散曲线数据进行反演，我们的目标函数如下
|
(4) |
其中VS为待求S波速度模型，i为频率序号，j为频散曲线的阶数，cijs为实测相速度，cij为反演过程中当前模型的理论相速度，我们通常需要对不同阶的频散曲线设置不同的权重weightj.理论频散曲线的计算采用Chen在1993年提出的利用广义反透射系数计算面波问题的算法.相速度是S波速度等模型参数的函数.计算目标函数对S波速度的偏导数，我们可以使用链式求导法则
所使用的背景噪声数据是日本National Research Institute for Earth Science and Disaster Resilience(NIED)在日本首都圈地区布设的地震台网MeSO-net的垂直分量数据.MeSO-net台网一共包含298个地震台站，台站分布如图 1所示，图中的黑色三角形代表地震台站.MeSO-net地震台网于2017年4月1日开始工作.选取其中部分台站2018年1月1日至2018年6月30日一共6个月的垂直分量连续波形数据进行处理.首先将数据以天为单位进行划分，对每天的数据进行频率-贝塞尔变换，然后将得到的能量谱进行叠加.在数据预处理阶段，我们需要对数据进行去均值、去势、重采样、滤波、时域归一化和谱白化等处理，根据区域大小和台网的分布，我们的重采样频率是5 Hz，滤波器是0.01~2 Hz的带通滤波器，值得注意的是MeSO-net台网数据的有效频段是DC-80 Hz，这个频段远远宽于我们研究的频段，因此不需要做仪器响应.
|
图 1 MeSO-net台网台站空间分布图 Fig. 1 Map showing distribution of MeSO-net stations |
频率-贝塞尔变换法是以一维水平层状结构为基础的，如果研究区域过大，一维水平层状结构的近似将不成立(Wang et al., 2019).综合考虑区域的大小和台站数量等因素，我们将关东盆地划分为A、B、C、D、E和F一共6个区域分别进行研究.区域的划分如图 2所示，不同颜色的矩形框是对应颜色字母对应区域范围，对应颜色的五角星位置代表该区域的台站中心位置.各个区域的台站中心的经纬度和台站数目如表 1所示.
|
图 2 研究区域划分图 Fig. 2 Map showing block subdivision of the study area |
|
表 1 区域台站信息 Table 1 Information of regional stations |
将频率-贝塞尔变换法的算法分别应用于6个区域，得到6个区域的F-J能量谱.6个区域的F-J能量谱如图 3所示.图 3(a—f)分别是区域A、B、C、D、E和F的频率-速度能量谱.6个子图均可以看到，频率-贝塞尔变换算法提取的F-J能量谱有不止一条连续的能量峰值连线，其中两条明显的峰值连线对应基阶和一阶频散曲线，频散曲线的有效范围为0.15~0.6 Hz，在F区域可以高达0.8 Hz.图 3(b—f)的能量谱的频散曲线均比较连续，能量峰值较窄，很容易并且能够比较准确地提取频散曲线，但图 3a的能量谱相对粗糙，这增加了提取频散曲线的难度，造成此现象的主要原因是A区域的台站数目太少，对此Wang等(2019)已进行过相应的讨论.图 3b和3e的能量谱都出现基阶频散曲线和一阶频散曲线在某一频率处无限接近的现象，图 3b出现在约0.25 Hz处，图 3e出现在约0.2 Hz处，这造成了我们无法区分低于该频率处的能量峰值连线是基阶频散曲线还是一阶频散曲线，这一现象是由地下结构决定的.由于反演需要明确频散曲线的阶数，错误地提取频散曲线，把基阶频散曲线当作一阶频散曲线或把一阶频散曲线当作基阶频散曲线都会对反演的结果造成很大的影响，因此我们需要在提取频散曲线时仔细甄别，多做尝试.
|
图 3 用频率-贝塞变换法得到的区域A-F对应的F-J能量谱 Fig. 3 F-J spectrograms using F-J method corresponding to blocks A-F |
得到F-J能量谱后，我们从中手动拾取频散曲线，得到一系列的频率-相速度对，然后利用拟牛顿法反演S波速度结构.
由于拟牛顿法反演可能陷入局部极小，因此选取合适的初始模型变得至关重要.Koketsu等(2008, 2009)使用地质填图、地震反折射、重力调查、地震尾波的H/V反演、钻井剖面以及微震调查等方法得出全日本的综合三维地下结构，命名为JIVSM.我们以JIVSM为参考，对区域内JIVSM的三维速度结构做一定的平均和近似，每个区域得到40个一维模型，将这些模型当作反演的初始模型，反演得到40个S波速度结构模型，将这些模型取平均当作最终的反演结果.
根据JIVSM模型，关东盆地的沉积层和地震基岩层可分为四层，每一层的参数如表 2所示.
|
|
表 2 JIVSM关东盆地模型参数表 Table 2 Parameters of JIVSM in the Kanto Basin |
将模型以0.4 km为间隔，将模型划分为13层，最后一层为4.8 km~∞，在反演过程中，层数和层厚不变.
当我们仅使用基阶频散曲线反演的时候，如图 4所示，黑色的点线代表真实数据，红色实线代表反演模型模拟得到的频散曲线，可以看到基阶频散曲线得到很好的拟合，但是一阶频散曲线的拟合效果不好.反演得到的S波速度结构模型如图 5所示，灰色虚线表示初始模型的取值范围，红色实线表示20个误差最小的反演模型，黑色实线表示反演的平均模型.可以看到，当仅利用基阶频散曲线反演的时候，得到的结果不稳定，多解性问题严重.
|
图 4 基阶频散曲线反演结果的频散曲线拟合情况 Fig. 4 Dispersion curves fitting using fundamental-mode dispersion curve inversion |
|
图 5 基阶频散曲线反演结果的S波速度结构模型 Fig. 5 S-wave velocity models from fundamental-mode dispersion curve inversion |
使用基阶频散曲线和高阶频散曲线联合反演时，频散曲线的拟合情况如图 6所示.可以看到基阶频散曲线和一阶频散曲线均能得到较好拟合，区域A的频散曲线拟合情况较差，这与A区域的台站数目较少从而导致提取的频散曲线较为粗糙有关.而利用基阶频散曲线和高阶频散曲线联合反演得到S波速度结构模型如图 7所示.可以看到，相较于仅利用基阶频散曲线得到的S波速度结构，基阶面波和高阶面波联合反演得到的速度结构更为收敛，不确定性更少.
|
图 6 基阶面波联合高阶面波反演结果的频散曲线拟合情况 Fig. 6 Dispersion curves fitting using fundamental-mode and higher mode dispersion curves joint-inversion |
|
图 7 基阶面波联合高阶频散曲线反演结果的S波速度结构模型 Fig. 7 S-wave velocity models from fundamental-mode and higher-mode dispersion curves joint-inversion |
取40次反演的S波速度模型的平均值作为最终S波速度结果，A-F区域的反演结果与JIVSM的S波速度结构如图 8所示，红线代表反演的S波速度模型，绿线代表JIVSM的S波速度模型.
|
图 8 平均反演模型和JIVSM对比图 Fig. 8 Average inversion S-wave models estimated by JIVSM |
Suzuki(1996)根据地质调查和钻井资料，提出关东盆地的沉积层可分为三层，最表层为Shimosa层，该沉积层为中后更新世沉积，S波速度约0.4~0.6 km·s-1；Shimosa层下是Kazusa层，该层是上新世到早更新世沉积，S波速度约1.0 km·s-1；最下层是Miura层，S波速度约1.7 km·s-1；Miura层下是地震基岩层，S波速度约为3.0 km·s-1.
从图 7中可以看到，反演得到的地震基岩层界面深度与JIVSM的地震基岩层差别不大.区域A、B、C、E的地震基岩层上界面较深，约为3 km，即这几个区域的沉积层厚度约3 km；区域D和区域F的地震基岩较浅，沉积层厚度约1.8 km.这与前人所做的多个模型的结论一致，即关东盆地东北地区的沉积层厚度较小；Miura层沉积厚度从关东盆地的西南往东北方向逐渐减少，A区域处的Miura层厚达2 km.
5 结语我们将频率-贝塞尔变换法应用于关东盆地的MeSO-net台网的背景噪声数据中，从背景噪声中提取基阶和高阶频散曲线.然后利用基阶频散曲线联合高阶频散曲线反演关东盆地沉积层以及基岩层的S波速度结构.最终证明频率-贝塞尔变换法能有效地从背景噪声中提取出基阶和高阶频散曲线，在基阶频散曲线的基础上，高阶频散曲线能有效地减少反演中的不确定性，得到更为准确稳定的速度模型.但是频率-贝塞尔变换是基于一维水平层状结构的，当研究区域过大，地下结构不满足一维水平层状结构的假设时，并不能提取到清晰的频散曲线，因此在保证台站数目的情况下，进行分块处理会取得更好的效果.
致谢 感谢日本National Research Institute for Earth Science and Disaster Resilience (NIED)提供的数据.
Aki K. 1957. Space and time spectra of stationary stochasticwaves, with special reference to microtremors. Bulletin of the Earthquake Research Institute, 35: 415-416. |
Aki K, Richards P G. 1980. Quantitative Seismology: Theory and Methods. San Francisco: W. H. Freeman and Company.
|
Capon J. 1969. High-resolution frequency-wavenumber spectrum analysis. Proceedings of the IEEE, 57(8): 1408-1418. DOI:10.1109/PROC.1969.7278 |
Chen X. 1999. Seismogram synthesis in multi-layered half-space Part Ⅰ. theoretical formulations. Earthquake Research in China, 13(2): 53-78. |
Chen X F. 1993. A systematic and efficient method of computing normal modes for multilayered half-space. Geophysical Journal International, 115(2): 391-409. DOI:10.1111/j.1365-246X.1993.tb01194.x |
He Y F. 2005. A preliminary research on inverse problems of surface wave dispersion[Master's thesis](in Chinese). Beijing: Peking University.
|
Kennett B L N. 1985. Seismic wave propagation in stratified media. Cambridge: Cambridge University Press.
|
KoketsuK, Miyake H, Fujiwara H, Hashimoto T. 2008. Progress towards a Japan integrated velocity structure model and long-period ground motion hazard map.//Proceedings of the 14th World Conference on Earthquake Engineering, S10-038, Beijing, October 12-17.
|
Koketsu K, Miyake H, Afnimar, et al. 2009. A proposal for a standard procedure of modeling 3-D velocity structures and its application to the Tokyo metropolitan area, Japan. Tectonophysics, 472(1-4): 290-300. DOI:10.1016/j.tecto.2008.05.037 |
Lacoss R T, Kelly E J, Toksöz M N. 1969. Estimation of seismic noise structure using arrays. Geophysics, 34(1): 21-38. DOI:10.1190/1.1439995 |
Luo Y H, Xia J H, Liu J P, et al. 2008. Joint inversion of fundamental and higher mode Rayleigh waves. Chinese Journal of Geophysics (in Chinese), 51(1): 242-249. |
Luo Y H, Xu Y X, Liu Q S, et al. 2008. Rayleigh-wave dispersive energy imaging and mode separating by high-resolution linear Radon transform. The Leading Edge, 27(11): 1536-1542. DOI:10.1190/1.3011026 |
Luo Y H, Xia J H, Miller R D, et al. 2009. Rayleigh-wave mode separation by high-resolution linear Radon transform. Geophysical Journal International, 179(1): 254-264. DOI:10.1111/j.1365-246X.2009.04277.x |
McMechan G A, Yedlin M J. 1981. Analysis of dispersive waves by wave field transformation. Geophysics, 46(6): 869-874. DOI:10.1190/1.1441225 |
Okada H. 2006. Theory of efficient array observations of microtremors with special reference to the SPAC method. Exploration Geophysics, 37(1): 73-85. DOI:10.1071/EG06073 |
Pan L, Chen X F, Wang J N, et al. 2019. Sensitivity analysis of dispersion curves of rayleigh waves with fundamental and higher modes. Geophysical Journal International, 216(2): 1276-1303. DOI:10.1093/gji/ggy479 |
Park C B, Miller R D, Xia J H. 1998. Imaging dispersion curves of surface waves on multi-channel record.//68th Ann. Internat Mtg., Soc. Expi. Geophys.. Expanded Abstracts, 1337-1380.
|
Park C B, Miller R D, Xia J H. 2000. Advantages of calculating shear-wave velocity from surface waves with higher modes.//70th Ann. Internat Mtg., Soc. Expi. Geophys.. Expanded Abstracts.
|
Park C B, Miller R D, Xia J H, et al. 2007. Multichannel analysis of surface waves (MASW)-active and passive methods. The Leading Edge, 26(1): 60-64. DOI:10.1190/1.2431832 |
Sanchez-Sesma F J, Campillo M. 2006. Retrieval of the Green's function from cross correlation:The canonical elastic problem. Bulletin of the Seismological Society of America, 96(3): 1182-1191. DOI:10.1785/0120050181 |
Sato H, Fehler M C, Maeda T. 2012. Green's function retrieval from the cross-correlation function of random waves.//Seismic Wave Propagation and Scattering in the Heterogeneous Earth: Second Edition. Berlin Heidelberg: Springer.
|
Suzuki H. 1996. Geology of the Koto deep borehole observatory and geological structure beneath the Metropolitan Area, Japan (In Japanese). Report of the National Research Institute for Earth Science and Disaster Prevention, 56, 77-123.
|
Tian Y, Chen X F. 2006. Simultaneous inversion of hypocenters and velocity using the quasi-Newton method and trust region method. Chinese Journal of Geophysics (in Chinese), 49(3): 845-854. |
Wang J N, Wu G X, Chen X F. 2019. Frequency-Bessel transform method for effective imaging of higher-mode Rayleigh dispersion curves from ambient seismic noise data. Journal of Geophysical Research:Solid Earth, 124(4): 3708-3723. DOI:10.1029/2018JB016595 |
Xia J H, Miller R D, Park C B, et al. 2003. Inversion of high frequency surface waves with fundamental and higher modes. Journal of Applied Geophysics, 52(1): 45-57. DOI:10.1016/S0926-9851(02)00239-2 |
Xia J H, Xu Y X, Miller R D. 2007. Generating an image of dispersive energy by frequency decomposition and slant stacking. Pure and Applied Geophysics, 164(5): 941-956. DOI:10.1007/s00024-007-0204-9 |
Yamanaka H, Yamada N. 2002. Estimation of 3D S-wave velocity model of deep sedimentary layers in Kanto plain, Japan, using microtremor array measurements. Geophysical Exploration (Butsuri-Tansa), 55: 53-65. |
Yao H J, Van DerHilst R D, DeHoop M V. 2006. Surface-wave array tomography in SE Tibet from ambient seismic noise and two-station analysis-Ⅰ. Phase velocity maps. Geophysical Journal International, 166(2): 732-744. DOI:10.1111/j.1365-246X.2006.03028.x |
Yao H J, Xu G, Xiao X, et al. 2004. A quick tracing method based on image analysis technique for the determination of dual stations phase velocities dispersion curve of surface wave. Seismological and Geomagnetic Observation and Research, 25(1): 1-8. |
何耀锋. 2005.面波频散反演问题的初步研究[硕士论文].北京: 北京大学. http://www.wanfangdata.com.cn/details/detail.do?_type=degree&id=Y814623
|
罗银河, 夏江海, 刘江平, 等. 2008. 基阶与高阶瑞利波联合反演研究. 地球物理学报, 51(1): 242-249. DOI:10.3321/j.issn:0001-5733.2008.01.030 |
田玥, 陈晓非. 2006. 利用拟牛顿法和信赖域法联合反演震中分布与一维速度结构. 地球物理学报, 49(3): 845-854. DOI:10.3321/j.issn:0001-5733.2006.03.029 |
姚华建, 徐果明, 肖翔, 等. 2004. 基于图像分析的双台面波相速度频散曲线快速提取方法. 地震地磁观测与研究, 25(1): 1-8. DOI:10.3969/j.issn.1003-3246.2004.01.001 |