2017, Vol. 32
2. 中国建筑科学研究院, 北京 100013
3. 中铁西北科学研究院有限公司, 兰州 730000
4. 河南亚新投资集团, 郑州 450052
2. China Academy of Building Research, Beijing 100013, China
3. Northwest Research Institute Co., Ltd of C. R. E. C, Lanzhou 730000, China
4. Henan Yasin Holdings, Zhengzhou 450052, China
21世纪以来，瑞利面波勘探方法 (Park et al., 1999; Aki and Richards, 2002) 以其可在复杂干扰波场条件下能够有效的提取地质信息的能力，从而在浅层工程地质勘探和工程质量检测等领域中得到了广泛应用.实验室物理模拟研究是认识和发展瑞利面波理论和方法的主要手段之一.波场实验相似准则是实验室设计物理模型的尺度的确立依据，而超声地震模型实验技术解决了在实验室小尺度模型上研究勘探地震学问题 (赵鸿儒等，1986；彭一民和孙进忠，1987；孙进忠等，1997).在瑞利面波的理论和模拟研究中，常将三维地质模型简化为二维模型来处理，以此降低分析问题的难度和节约建模成本.但常规的实验室二维物理模型的制作需要在走向方向上有足够长的延伸，以近似走向方向为无限延伸的二维理论要求.实验研究表明，在Rayleigh面波实验中，走向方向的延伸应大于4~5倍的实验波长，才能满足这一近似条件 (韩超等，2016).这导致建模难度和成本的急剧增高.但基于地震波场理论，地震波的横向分辨率或水平分辨率是以第一菲涅尔带半径R为度量的 (于振清等，1990；王绪松等，2006).据此，若实验物理模型的走向方向上延伸或模型厚度小于第一菲涅尔带半径R，则实验超声波应该不能识别模型的厚度，即走向方向的维度影响可以忽略.因此从理论上可以推断：满足该条件的物理模型可视为二维物理模型.显然，若在实验中的三维模型厚度小于第一菲涅尔带半径R, 则可近似为二维模型.这样，一方面可以极大的削减建模难度和成本，另一方面可实现二维物理模拟的结果更好的与二维的理论或数值模拟结果匹配分析.因此，该方法的实际可行性研究对于促进Rayleigh面波物理实验技术和理论发展都具有和重要的意义.
本文首先基于数值模拟方法，模拟研究了具有不同菲涅尔带半径R厚度的三维均匀半空间和横向非均匀模型的Rayleigh面波波场和频散特征，并与对应的二维数值模拟结果进行对比分析.然后，以数值模拟研究确定的合理的菲涅尔带半径R厚度制作三维模型，对比分析该模型和二维数值模拟的Rayleigh面波波场和频散特征.
1 理论基础 1.1 第一菲涅尔带半径点震源发出的球面波到达界面时的波前面，与前面相距1/4波长先期到达的另一波前面在界面上形成的圆称第一菲涅尔带，该圆的半径称为第一菲涅尔带半径，记为R(姚姚，2006).如图 1所示.目前，理论上把第一菲涅尔带半径R作为地震资料横向分辨力 (王绪松等，2006；于振清等，1990).也就是说第一菲涅尔带表达的含义是一个反射体的长度a超过第一菲涅尔带半径R(a>R) 时才能够被地震波识别.第一菲涅尔带半径的计算公式为
|
图 1 第一菲涅尔带半径示意图 Figure 1 The first Fresnel zone radius |
|
(1) |
式中f为地震波频率，λ为地震波波长，h为点震源到界面的深度.
根据第一菲涅尔带原理，设计的实验物理模型的厚度a若小于第一菲涅尔带半径R，则可推断超声波应该不能识别模型的厚度，即满足该条件的三维物理模型可近视为二维物理模型.
1.2 数值模拟方法本文采用商用Abaqus有限元软件进行数值模拟分析.Abaqus是目前较为先进的大型有限元软件之一，其具有强大、广泛的计算模拟能力，并且能够驾驭高度非线性问题，适合瑞利面波波场模拟.数值模拟过程中采用Abaqus/Explicit分析模块，该模块采用显式动力有限元列式，适用于像冲击和爆炸这类短暂，瞬时的动态事件 (庄茁等，2005；石亦平和周玉蓉，2006).
数值模拟中采用雷克子波作为震源，雷克子波的加载通过编写用户子程序uamp予以实现.雷克子波表达式为
|
(2) |
其中，t0表示激励震源的雷克子波主峰时刻，f0表示激励震源的主频率.在以下的数值模拟中均采用和实验室物理模型实验中的激发超声声源相近主频f0=48 kHz的雷克子波作为波场激发源.
研究中希望对应力波的传播过程能够正确的模拟，综合考虑计算量、计算时间等因素后，选用线性减缩积分单元CPS4R (4节点四边形平面应力线性减缩积分单元)、C3D8R (8节点六面体线性减缩积分单元) 单元作为二维、三维模型的单元类型.
单元网格划分的大小对瑞利面波波场模拟的结果精度和收敛性都有很大的影响，所以单元网格划分大小一般采用主波长的1/20，网格达到此要求时，模拟精度较高，计算量在可控范围内.
2 数值模拟分析 2.1 均匀介质数值模型 2.1.1 模型和观测参数基于实验室中常用的模型材料，本文主要采用了有机玻璃的物性参数进行数值模拟建模，具体材料参数见表 1.
|
表 1 模型材料参数 Table 1 Model material parameters |
为较系统的分析模型厚度的影响，本节以具有不同菲涅尔带半径R厚度的三维均匀半空间和对应的横向非均匀半空间的数值模型，研究了二维物理模拟的模型厚度和影响因素.
依据表 1中的有机玻璃模型的物性参数，建立了5个均匀介质数值模型 (模型1、2、3、4、5)，其长度和高度 (或深度) 相同，但模型的厚度不同.具体模型尺度参数见表 2.表 2中模型1为二维模型，故其厚度以0表示.由有机玻璃材料的模型参数和主频为48 kHz的雷克子波可计算，激发的面波主波长为2.5 cm左右，远小于上述均匀介质模型的41 cm高度.因此，对于瑞利面波的传播，建立的数值模型在深度尺度上是可视为均匀介质模型的.在综合考虑近场效应、空间假频影响和三维模拟的有限元计算量问题后，选用了如图 2所示的瞬态多道瑞利面波实验观测排列，以及1.024 ms采样窗长，0.5 μs采样间隔的观测参数.为更清楚的从理论上分析模型厚度的影响, 依据图 1中观测排列的12 cm最小偏移距和模拟的48 kHz的激发源主波长，由 (1) 式计算出观测排列的最小第一菲涅尔带半径R为3.9 cm.表 2中显示了以最小第一菲涅尔带半径R归一后的各模型厚度.
|
|
表 2 数值模型尺度列表 Table 2 The size parameters of homogeneous medium models |
|
图 2 有机玻璃数值模型观测排列示意图 Figure 2 The observation layout on the top surface of the lucite model |
图 3a，c，e，g，i分别显示了模型1，2，3，4，5的数值模拟单炮观测记录的垂直分量剖面.由图 3c，e，g可看到，厚度远小于最小第一菲涅尔带半径R的模型2，3，4，其单炮记录上可清楚的识别出直达P波、直达面波，其特征与图 3a中的二维模型-模型1记录非常相近.而厚度大于最小第一菲涅尔带半径R的模型5，图 3i中显示的其直达面波记录与二维模型1的记录在具有很大差异.根据之前研究结果 (韩超等，2016)，这是由于在瑞利面波的物理模拟实验中，模型厚度较大时侧边界的反射面波的叠加干扰所致.图 3b，d，f，h，j显示了各模型单炮记录的频率-速度 (f-v) 域的面波能量分布图或称为面波频散谱 (张碧星等，2002；张碧星和鲁来玉，2005).在面波频散谱提取前, 采用了常规的视速度开窗滤波 (各单炮记录剖面中的阴影区) 对各模型的单炮记录进行干扰波切除处理.从瑞利面波频散谱中可以看到，图 3b，d，f，j面波能量频散分布特征极为相似，均表现出均匀半无限空间无频散特征，其拾取的无频散速度为1200 m/s左右，非常接近于表 1中有机玻璃的理论面波速度.但是对于着模型厚度的大于最小第一菲涅尔带半径R的模型5，由于侧边界的反射面波的叠加干扰，其瑞利面波频散谱 (图 3i) 呈现出强的频散特征，已经不能反映出均匀半无限空间介质的信息了.
|
图 3 均匀介质模型数值模拟记录和面波频散分析图 (a) 模型1的数值模拟单炮记录和干扰波切除视速度窗口; (b) 干扰波切除后的f-v域图像; (c) 模型2的数值模拟单炮记录和干扰波切除视速度窗口; (d) 干扰波切除后的f-v域图像; (e) 模型3的数值模拟单炮记录和干扰波切除视速度窗口; (f) 干扰波切除后的f-v域图像; (g) 模型4的数值模拟单炮记录和干扰波切除视速度窗口; (h) 干扰波切除后的f-v域图像; (i) 模型5的数值模拟单炮记录和干扰波切除视速度窗口; (j) 干扰波切除后的f-v域图像. Figure 3 The shot records and the dispersion analysis data of the Rayleigh waves of the homogeneous medium model with numerical data (a) The common shot records of model 1; (b) The image of direct arrival Rayleigh waves in the f-v domain; (c) The common shot records of model 2; (d) The image of direct arrival Rayleigh waves in the f-v domain; (e) The common shot records of model 3; (f) The image of direct arrival Rayleigh waves in the f-v domain; (g) The common shot records of model 4; (h) The image of direct arrival Rayleigh waves in the f-v domain; (i) The common shot records of model 5; (j) The image of direct arrival Rayleigh waves in the f-v domain. |
以下，我们以含矩形空洞的横向非均匀模型的数值模拟结果来进一步分析在复杂条件下上述结论是否成立.
2.2.1 模型和观测参数为类比，在上节设计的5个具有不同最小第一菲涅尔带半径R的均匀半空间模型均插入了一相同尺度和位置的矩形空洞而构成5个横向非均匀数值模型 (见表 2中模型6、7、8、9、10)，矩形空洞具体几何尺度和位置见图 2.在数值模拟中，采用的观测排列参数和激发震源参数与上节中均匀介质模型完全一致.
2.2.2 数值模拟结果分析图 4a，c，e，g，i分别为模型6，7，8，9，10的数值模拟的单炮观测记录的垂直分量剖面.与前述的均匀介质模型类似，由图 4c，e，g可看到，厚度远小于最小第一菲涅尔带半径R的模型7，8，9的单炮记录上可清楚的识别出直达P波，直达面波，其特征与图 4a中的二维模型-模型6记录非常相近.而图 4i中显示的厚度大于最小第一菲涅尔带半径R的模型10的单炮记录中, 由于存在模型侧边界的反射面波干扰，其直达面波记录与二维模型-模型6的记录存在很大的差异.图 4b，d，f，h，j为经简单的视速度开窗滤波 (单炮记录图中的阴影区) 进行干扰波切除处理后，各模型单炮记录的频率-速度 (f-v) 域的面波能量分布图.图 4b，d，h显示，厚度远小于最小第一菲涅尔带半径R的模型7，8所示的面波能量分布特征与图 4a中的二维模型-模型6的极为相似，都出现了明显的“之”字形曲线，且“之”字形拐点处频率大致为60 kHz，对应面波波长为2 cm，与模型实际空洞埋深范围较一致.而对于模型厚度的大于最小第一菲涅尔带半径R的模型10, 由于存在强烈的侧边界反射面波干扰，图 4j所示的面波能量分布的“之”字形曲线与二维模型-模型6的明显不同.
|
图 4 横向非均匀介质模型数值模拟记录和面波频散分析图 (a) 模型6的数值模拟单炮记录和干扰波切除视速度窗口; (b) 干扰波切除后的f-v域图像; (c) 模型7的数值模拟单炮记录和干扰波切除视速度窗口; (d) 干扰波切除后的f-v域图像; (e) 模型8的数值模拟单炮记录和干扰波切除视速度窗口; (f) 干扰波切除后的f-v域图像; (g) 模型9的数值模拟单炮记录和干扰波切除视速度窗口; (h) 干扰波切除后的f-v域图像; (i) 模型10的数值模拟单炮记录和干扰波切除视速度窗口; (j) 干扰波切除后的f-v域图像. Figure 4 The shot records and the dispersion analysis data of the Rayleigh waves of the inhomogeneous medium model with numerical data (a) The common shot records of model 6; (b) The image of direct arrival Rayleigh waves in the f-v domain; (c) The common shot records of model 7; (d) The image of direct arrival Rayleigh waves in the f-v domain; (e) The common shot records of model 8; (f) The image of direct arrival Rayleigh waves in the f-v domain; (g) The common shot records of model 9; (h) The image of direct arrival Rayleigh waves in the f-v domain; (i) The common shot records of model 10; (j) The image of direct arrival Rayleigh waves in the f-v domain. |
基于上述数值模拟研究结果和观测的可行性，本节以厚度小于最小第一菲涅尔带半径R的有机玻璃材料分别制作了均匀和横横向非均匀三维物理模型，在实验室中进行了超声物理模拟实验，并将观测到的Rayleigh面波波场和频散特征与二维数值模拟的进行了对比分析.
3.1 均匀介质物理模型 3.1.1 模型和观测参数参考前述的数值模拟结果和物理模拟的实验参数，制作了表 3中所示尺度的完整有机玻璃 (模型11) 的均匀介质物理模型.对比表 2中数据, 可知该物理模型与数值模型2具有完全相同的几何参数和物性参数.
|
|
表 3 物理模型尺度列表 Table 3 The size parameters of homogeneous medium models |
物理模拟的观测排列布置参数与数值模拟中近似 (见图 5).物理模拟中采用一48 kHz的超声波换能器作为声激发源，一宽带超声换能器为接收检波器.记录窗长和采样间隔与数值模拟中参数相同.考虑到实验中的随机噪声干扰，每道采样记录均以1024次重复激发-采样的记录叠加作为最终记录，该重复激发-采样记录和叠加过程由计算机控制的虚拟仪器自动完成 (Xiong et al., 2010, 2013).
|
图 5 有机玻璃物理模型观测排列示意图 Figure 5 The observation layout on the top surface of the lucite model |
图 6显示了物理模型-模型11的实验观测得到的单炮记录和按数值模拟研究中的相同参数及数据处理流程进行数据处理的Rayleigh面波的能量谱图像.图 6a的单炮记录清楚的显示了直达P波，直达面波，其与二维数值模拟的均匀介质模型-模型1的数值模拟结果 (图 3a) 具有相似的特征.图 6b中的物理模拟的面波能量谱分布特征与图 3b中的二维数值模型的相似，均表现出相速度为1200 m/s的无频散均匀半无限空间特征.
|
图 6 均匀介质物理模拟记录和面波频散分析图 (a) 数值模拟单炮记录和干扰波切除的视速度窗口; (b) 干扰波切除后的f-v域图像. Figure 6 The records and the dispersion analysis data of the Rayleigh waves of the homogeneous medium model with physical modeling (a) The common shot records; (b) The image of direct arrival Rayleigh wavesin the f-v domain. |
在上节的均匀半空间模型基础上, 设计了一圆形空洞制作为横向非均匀介质物理模型 (模型12).圆形空洞半径2 cm，位置如图 5所示.该物理模型的观测排列，模拟中采用的声激发源和接收检波器、记录窗长、采样间隔、数据采集和处理方法等均与上节均匀介质物理模型一致 (见图 5).
3.2.2 物理模拟结果分析图 7显示了物理模型-模型12的实验观测得到的单炮记录和按数值模拟研究中的相同参数及数据处理流程进行数据处理的Rayleigh面波的能量谱图像.图 7b中，面波能量分布图中出现了明显的“之”字形曲线，且“之”字形曲线拐点处频率大致为30 kHz，对应面波波长约为4 cm，与空洞实际大小位置一致，正确的反映了真实的物理模型特征.
|
图 7 横向非均匀介质物理模拟记录和面波频散分析图 (a) 数值模拟单炮记录和干扰波切除的视速度窗口; (b) 干扰波切除后的f-v域图像. Figure 7 The records and the dispersion analysis data of the Rayleigh waves of the inhomogeneous medium model with physical modeling (a) The common shot records; (b) The image of direct arrival Rayleigh wavesin the f-v domain. |
本文以第一菲涅尔带半径R的波场理论为指导, 基于均匀半空间和横向非均匀半空间模型的数值模拟和物理模拟，研究了瑞利面波物理模拟实验中模型厚度对瑞利面波波场和频散曲线提取的影响.研究结果表明，在模型厚度小于最小第一菲涅尔带半径R条件下, 三维物理模型实验结果与二维数值模拟模型实验结果基本一致.因此，在实际物理模拟实验中, 可用厚度小于最小第一菲涅尔带半径R的三维物理模型很好的近似二维模型，既进行二维物理模型实验.这不仅极大的节约物理建模成本、降低建模难度本，也可实现二维物理模拟的结果更好的与二维的理论或数值模拟结果的匹配分析.研究的结果还表明，物理模型的厚度选择还受限于物理实验中超声换能器的半径尺度，尤其是厚度小于接受超声换能器半径的过薄模型会导致观测难以工作的问题.
致谢 本项目由国家自然科学基金—山区深部岩体结构的瑞利面波椭圆极化特征研究 (41274107) 资助.本文的部分观点由麻省理工大学朱正亚教授提出，在此深表谢意.在数值模型的设计、研究方法和论文书写方面，得到了于文福师兄和熊自英师姐的热心帮助，在此深表谢意.对于文福师兄前期工作和平台的搭建所作的贡献，衷心感谢.| [] | Aki K, Richards P G. 2002. Quantitative Seismology[M]. Sausalito, CA: University Science Books. |
| [] | Han C, Liu Z P, Xiong Z Y, et al. 2016. An approach to the model scale influence in the experiments of multi-channel Rayleigh waves[J]. Progress in Geophysics (in Chinese), 31(3): 1313–1319. DOI:10.6038/pg20160353 |
| [] | Park C B, Miller R D, Xia J H. 1999. Multichannel analysis of surface waves[J]. Geophysics, 64(3): 800–808. DOI:10.1190/1.1444590 |
| [] | Sun J Z, Guo T S, Tang W B, et al. 1997. Theory, practice and application of ulltrasonic seismic modeling experiment in China[J]. Chinese Journal of Geophysics (in Chinese), 40(S): 266–274. |
| [] | Xiong Z Y, Liu Z P, Song B, et al. 2010. Physical experiments for Rayleigh surface wave field[C].//Near-Surface Geophysics and Geohazards-Proceedings of the 4th International Conference on Environmental and Engineering Geophysics (Volume 1). Chengdu, China. |
| [] | Xiong Z Y, Liu Z P, Zhang K. 2013. An experimental research of Rayleigh wave based on the acoustic-electric conversions effect[C].//Near Surface Geophysics Asia Pacific Conference. Beijing, China:SEG, 258-261. |
| [] | Zhang B X, Lu L Y, Bao G S. 2002. A study on zigzag dispersion curves in Rayleigh wave exploration[J]. Chinese Journal of Geophysics (in Chinese), 45(2): 263–275. DOI:10.3321/j.issn:0001-5733.2002.02.013 |
| [] | Zhang B X, Lu L Y. 2005. Investigation on the dispersion curves of Rayleigh wave by frequency-wavenumber analysis method[J]. Chinese Journal of Engineering Geophysics (in Chinese), 2(4): 245–255. |
| [] | 韩超, 刘争平, 熊自英, 等. 2016. 瞬态多道瑞利面波实验中的模型尺度影响研究[J]. 地球物理学进展, 31(3): 1313–1319. DOI:10.6038/pg20160353 |
| [] | 彭一民, 孙进忠. 1987. 超声地震模型实验相似准则[J]. 勘察科学技术, (5):58-61, 52(5): 58–61, 52. |
| [] | 石亦平, 周玉蓉. 2006. ABAQUS有限元分析实例详解[M]. 北京: 机械工业出版社. |
| [] | 孙进忠, 郭铁栓, 唐文榜, 等. 1997. 我国超声地震模型实验的理论研究与实践[J]. 地球物理学报, 40(S): 266–274. |
| [] | 王绪松, 钱雪峰, 李波涛. 2006. 第一菲涅尔带半径与地震资料的横向分辨力[J]. 勘探地球物理进展, 29(4): 244–248. |
| [] | 姚姚. 2006. 地震波场与地震勘探[M]. 北京: 地质出版社. |
| [] | 于振清, 张玉芬, 候慎健, 等. 1990. 用菲涅尔理论解释反射地质体边界的理论模型研究[J]. 石油物探, 29(1): 63–73. |
| [] | 张碧星, 鲁来玉, 鲍光淑. 2002. 瑞利波勘探中"之"字形频散曲线研究[J]. 地球物理学报, 45(2): 263–275. DOI:10.3321/j.issn:0001-5733.2002.02.013 |
| [] | 张碧星, 鲁来玉. 2005. 用频率-波数法分析瑞利波频散曲线[J]. 工程地球物理学报, 2(4): 245–255. |
| [] | 赵鸿儒, 唐文榜, 郭铁栓. 1986. 超声地震模型试验技术及应用[M]. 北京: 石油工业出版社. |
| [] | 庄茁. 2005. ABAQUS非线性有限元分析与实例[M]. 北京: 科学出版社. |