2. 中国科学技术大学地球和空间科学学院,合肥 230026;
3. 西安交通大学电子与信息工程学院,西安 710049
2. School of Earth and Space Science, University of Science and Technology of China, Hefei 230026, China;
3. School of Electronic and Information Engineering, Xi'an Jiao Tong University, Xi'an 710049, China
由于自由表面效应与随深度而改变的覆盖层自重压力,地下核爆炸过程不能简单地用爆炸源来表示,而是增加一个补偿线性矢量偶极(CLVD)辅助源来最低阶地近似源区的复杂过程[1],用以描述核爆炸过程伴随的逆倾滑断层形式的构造应力释放[2].CLVD 源的物理意义可以简单地理解为纯张裂位错矩张量的偏量部分,利用其激发的Rg 波频谱中的低谷点是估算爆炸源深度的重要方法[3, 4].Langston依据模式理论从P 波和SV 波干涉抵消的角度演示了CLVD 源下Rg波低谷点的存在[5].Patton等人给出了CLVD 源深度与Rg波低谷点的解析关系[6],而CLVD 源的深度通常约为爆炸源深度的1/3[3, 4],据此可得到爆炸源深度.但是,CLVD源深度与Rg波低谷点的解析关系理论上只严格适用于半空间泊松体.何永峰、陈晓非等人对CLVD源Rg波低谷点产生的机理做了阐述并考察了东哈萨克斯坦地区一维平行分层模型下该解析关系的适用情况[7],并对区域Lg震相中的低谷点形成与Rg波的关系做了研究[8].
实际情况下,地表存在一定的地形起伏,此时,CLVD 源深度与Rg波低谷点的解析关系是否适用及其计算结果是否具有物理含义是值得考察的问题.实际工作中,地下核爆炸又常常在靠近山体的位置进行,利用地震学方法得到的爆炸源深度是否涵盖了地形因素是有争议的焦点问题.因此,利用数值模拟的方法研究地表起伏对CLVD 源Rg波低谷点的影响具有重要的学术意义和应用价值.
数值计算中,地表起伏属于复杂介质结构情况,基于一维分层模型的传统理论地震图算法不再适用.地球物理学家针对复杂介质结构问题发展了一系列数值算法,主要有纯数值法、高频射线法、边界积分方程法等.Zhou Hong和ChenXiaofei提出的局域离散波数边界积分法[9],朱德瀚等人提出的基于任意曲线网格的有限差分法[10]都是对原有算法的重大发展,能够较高效率地处理不规则介质结构问题.谱元法是最近20年新发展起来的地震波数值模拟技术,它兼具了有限元法处理复杂结构灵活性与伪谱法较高精度的优点.Komatitsch等人编写了从一维到三维的谱元计算程序,并将源代码公布在互联网上供全球科研人员下载研究使用.网格划分上,谱元法首先根据设定的分层面坐标及每层的网格数生成物理空间网格,计算过程中采用坐标变换映射到规则的自然坐标网格空间.由于采取了高阶的特殊的插值函数,使得质量矩阵恰好为对角矩阵,既避免了低阶有限元数值频散的问题,又利于并行计算,缩短了计算时间[11~14],并在实际问题的数值模拟中得到成功运用[15].本文利用二维谱元程序对地表起伏情况下CLVD 源Rg波低谷点的变化特征加以研究考察.
2 数值模拟研究 2.1 半空间泊松体情况考察介质设置为泊松体,密度为2600kg/m3,纵波速度为6.8km/s, 横波速度为3.9km/s.计算模型的尺度选取为宽度500km, 高度50km.震源坐标设在(50km, -1km),3个位于地表的接收点坐标为(250 km, 0 km),(350 km, 0 km),(450 km, 0km),从左至右编号为1~3.我们首先考察水平地表的情况,计算模型如图 1所示.
计算得到的各接收点垂向地震图如图 2a所示,由图看到,CLVD 源激发出的主要是Rg面波,前面的P波幅值相对很小.对Rg波做频谱分析,如图 2b所示.这里需要指出的是,本文考察的是地表起伏对Rg波频谱低谷点的影响,频谱的幅值没有什么含义,为了图形的显示清楚,我们将不同接收点的频谱在幅值上有个错距.由频谱图可以清楚地看到低谷点存在,各接收点Rg波频谱低谷点位置均为0.43Hz, 即低谷点位置不随震中距改变.根据Patton等人给出的CLVD 源深度与低谷点频率关系[6]:
(1) |
其中hCLVD为CLVD源深度,Vp为源区介质纵波速度,fNull为低谷点频率.计算得到震源深度为1km.可见,计算结果与理论相吻合.
2.2 源区山体起伏斜度对低谷点的影响在震源处的地表叠加一个山体,震源位于山体的正下方,距水平地表的距离为1km.为了作对比研究,这里分别考察两种情况:
(1) 山体的底部水平宽度为10km, 高度分别为0.5km, 1km;
(2) 山体的高度为1km, 底部水平宽度分别为10km, 20km.
计算模型如图 3所示.为了图形的显示清楚,山体的宽度比实际扩大了4倍.
图 4为山体水平宽度10km, 高度分别为0.5km、1km 情况的Rg 波位移谱对比.可以看到,两种情况下3个接收点的Rg波位移谱的低谷点位置均基本一致,没有随震中距发生变化,体现了波的传播特性.山体高度0.5km 时低谷点位置约为0.35Hz.根据(1)式,震源深度为1.23km.Rg波低谷点位置反映了源区的地形起伏特征变化情况.山体高度1km时低谷点变得不明显,可见如果地形起伏过于陡峭,低谷点趋于消失.此时低谷点对应的频率约为0.29Hz, 计算得到震源深度为1.48km.
图 5为山体水平宽度20km、高度1km 情况的Rg波位移谱.
对比图 5和图 4b,可以看出,在同样高度情况下,山体水平宽度为20km 时的低谷点效应较水平宽度为10km 时的明显.山体水平宽度20km 时,低谷点位置对应0.25Hz, 据此计算震源深度为1.72km, 而山体宽度10km 时,计算得到的震源深度为1.48km.可见,山体宽度增加后,计算得到的震源深度加深,更接近震源距山顶的距离.这个结果是合理的,如果山体的高度为1km, 宽度为无穷大,则相当于在地表覆盖了1km 厚的介质层,震源深度则变为2km.因此,当山体宽度增加时,震源深度应有接近2km 的趋势.
2.3 源区山体个数对低谷点的影响设置两个等高山体,其水平宽度为10km, 高度分别为0.5km、1km, 震源距水平地表的距离为1km, 计算模型如图 6所示.
图 7为两种情况的Rg 波位移谱对比.两个山体高度为0.5km 时低谷点对应的频率为0.33Hz, 两个山体高度1km 时低谷点对应的频率为0.29Hz.与一个山体计算结果相比,低谷点位置基本不变,低谷点更加明显了.可见,低谷点受源区山体起伏的影响十分复杂,可能与山体对自由表面施加的面力分布有关,需要进一步研究.
计算模型如图 8所示,山体水平宽度为10km, 高度1km.图 8a、8b、8b分别为震源距水平地表 0km、1km 情况,每种情况5个震源距左端山脚的水平距离分别为0km, 3km, 5km, 7km, 10km, 从左至右编号为1~5,3号震源位于山顶的正下方.接收点坐标为(450km, 0km).
图 9为震源相对于山体不同位置的Rg波位移谱图.图 9a为震源距水平地表 0km 情况的Rg波位移谱.从图中看到,位于山脚的1,5号震源情况的Rg波位移谱没有低谷点.这是合理的,根据式(1),震源深度为0,低谷点频率则为无穷大;山体内部的2~4号震源,其Rg波位移谱均有低谷点存在,距左端山脚3km 的2 号震源情况低谷点最明显,要大大显著于位于山顶正下方的3号震源情况,这与2.3节得到的两个山体情况下,即震源与山体位置不对称时,低谷点效应显著的结论相同.但是,与2 号震源处于对称位置的4号震源情况低谷点却不明显,可见,震源与山体的相对位置对低谷点有着复杂的影响,其原因为面波传播路径的地形起伏特征不同.本质上依旧是与山体对自由表面施加的面力分布有关.2~4号震源,根据其Rg波低谷点频率计算得到震源深度分别为0.72km, 0.88km, 0.63km.图 9b为震源距水平地表 1km 情况的Rg 波位移谱.此时,5个震源的Rg波位移谱均有低谷点,2 号震源的低谷点依旧最明显,1,4,5号震源情况的Rg波低谷点频率位置基本相同.1~5号震源,根据其Rg波位移谱低谷点频率计算得到震源深度分别为1.30km, 1.65km, 1.48km, 1.38km, 1.30km.
我们注意到,在两种计算模型情况下,利用同一接收点的Rg波低谷点计算得到的处于对称位置的2号震源与4 号震源深度结果却不相同,如果将接收点做相对于震源的180°平移,则可以想象2号震源与4号震源的深度结果会有对换.这说明利用Rg波低谷点计算震源深度结果与接收点的方位角有关.我们认为这是合理的,Rg波为沿自由表面传播的面波,其传播特性受地表起伏影响.不同方位角传播路径上的地表起伏特征不同,利用Rg 波中的低谷点计算得到的震源深度理应有所不同.
地形起伏情况下,震源的深度很难定义,如果我们定义震源到地表的竖向距离为理论值,则理论值与实际计算值的对比如图 10 所示.可以看到,计算值与理论值大小变化趋势整体上一致,但有一定偏差.
计算模型如图 11所示,在远区的传播路径上叠加一个山体,山体的水平宽度为10km, 高度分别为0.5km, 1km.接收点1位于山体的左侧,接收点2,3位于山体的右侧.
图 12为山体两种高度情况Rg波位移谱.由图看到,同样的地表起伏,如果放在远区位置,Rg波低谷点效应基本与半空间模型一样,低谷点位置为0.43Hz, 没有发生改变.可见,Rg波低谷点受远区地表起伏的影响较小.
本文利用谱元算法对山体起伏情况下利用Rg波低谷点计算震源深度的方法做了数值模拟研究.结果表明,在源区地表存在一定起伏的情况下,Rg波位移谱中的低谷点依然存在,不随震中距改变,低谷点位置反映了源区地表起伏的变化,但与理论值仍有一定偏差.山体越陡峭,低谷点越不明显,趋于消失;山体高度不变,宽度增加时,低谷点变得明显,并向低频移动,计算得到的震源深度增大.计算中发现,震源与山体位置不对称时,低谷点比较明显,其原因可能与山体对自由表面施加的面力分布有关,需要进一步研究.值得注意的是,利用Rg波低谷点计算得到的震源深度依赖于接收点的方位角.相似地形起伏放在远区位置,Rg波低谷点效应基本与水平地表时相同,没有受到影响.这些结论加深了我们对地表起伏情况下Rg波低谷点变化特征的认识,为在实际工作中约束爆炸震源深度提供了一定的理论参考.
由于低谷点频率处于分母位置,当低谷点频率较小时,深度对低谷点变化很敏感.经计算,低谷点频率低于0.6Hz时,低谷点位置0.01Hz变化就会带来深度10m 以上的改变,在实际深度定位中,应利用多种方法综合考虑以提高精度.本文是基于二维简单山体模型的模拟,至于山体的三维复杂形态对Rg波低谷点的影响需要进一步研究.
[1] | Stevens J L, Xu H M, Baker G E. An upper bound on Rg to Lg scattering using modal energy conservation. Bull. Seism. Soc. Am. , 2009, 99(2A): 906-913. DOI:10.1785/0120080213 |
[2] | 何永锋, 陈晓非, 张海明. 层裂对区域震相Lg波的影响. 地震学报 , 2005, 27(3): 309–316. He Y F, Chen X F, Zhang H M. The effect of spall on Lg waves. Acta Seismologica Sinica (in Chinese) (in Chinese) , 2005, 27(3): 309-316. |
[3] | Gupta I N, Chan W W, Zhang T R, et al. A comparison of regional phases from underground nuclear explosions at East Kazakh and Nevada Test Sites. Bull. Seism. Soc. Am. , 1992, 82(1): 352-382. |
[4] | Gupta I N, Zhang T R, Wagner R A. Low frequency Lg from NTS and Kazakh nuclear explosions-observations and interpretation. Bull. Seism. Soc. Am. , 1997, 87(5): 1115-1125. |
[5] | Langston C A. A note on spectral nulls in Rayleigh waves. Bull. Seism. Soc. Am. , 1980, 70(4): 1409-1414. |
[6] | Patton H J, Taylor S R. Analysis of Lg spectral ratios from NTS explosions: Implications for the source mechanisms of spall and the generation of Lg waves. Bull. Seism. Soc. Am. , 1995, 85(1): 220-236. |
[7] | 何永锋, 陈晓非, 何耀峰, 等. 地下爆炸Rg波低谷点激发机理. 地球物理学报 , 2005, 48(3): 643–648. He Y F, Chen X F, He Y F, et al. Generation of null in Rg wave by underground explosions. Chinese J. Geophys. (in Chinese) , 2005, 48(3): 643-648. DOI:10.1002/cjg2.697 |
[8] | 何永锋, 陈晓非, 张海明. 地下核爆炸Lg波的激发机制. 地球物理学报 , 2005, 48(2): 367–372. He Y F, Chen X F, Zhang H M. The excitation of Lg wave by underground nuclear explosions. Chinese J. Geophys. (in Chinese) , 2005, 48(2): 367-372. |
[9] | Zhou H, Chen X F. The localized boundary integral equation-discrete wavenumber method for simulating P-SV wave scattering by an irregular topography. Bull. Seism. Soc. Am. , 2008, 98(1): 265-279. DOI:10.1785/0120060249 |
[10] | 朱德瀚, 张伟, 祝贺君, 等. 利用曲线网格有限差分方法研究三维倾斜断层的破裂动力学. 地震学报 , 2010, 32(4): 401–411. Zhu D H, Zhang W, Zhu H J, et al. 3-D rupture dynamics of dipping faults with curvilinear-grid finite difference method. Acta Seismologica Sinica (in Chinese) (in Chinese) , 2010, 32(4): 401-411. |
[11] | Komatitsch D, Ritsema J, Tromp J. The spectral-element method, beowulf computing, and global seismology. Science , 2002, 298(5599): 1737-1742. DOI:10.1126/science.1076024 |
[12] | Komatitsch D, Tromp J. Introduction to the spectral element method for three-dimensional seismic wave propagation. Geophys. J. Int. , 1999, 139(3): 806-822. DOI:10.1046/j.1365-246x.1999.00967.x |
[13] | Komatitsch D, Tromp J. Specral-element simulations of global seismic wave propagation-I: Validation. Geophys. J. Int. , 2002, 149(2): 390-412. DOI:10.1046/j.1365-246X.2002.01653.x |
[14] | Komatitsch D, Tromp J. Spectral element simulations of global seismic wave propagation-II. Three-dimensional models, oceans, rotation and self gravitation. Geophys. J. Int. , 2002, 150(1): 303-318. DOI:10.1046/j.1365-246X.2002.01716.x |
[15] | Lee S J, Komatitsch D, Huang B S, et al. Effects of topography on seismic wave propagation: An example from northern Taiwan. Bull. Seism. Soc. Am. , 2009, 99(1): 314-325. DOI:10.1785/0120080020 |