地球物理学报  2018, Vol. 61 Issue (6): 2380-2395   PDF    
一种基于层析成像技术提高浅地表面波勘探水平分辨率的方法
尹晓菲1, 胥鸿睿2, 夏江海3, 孙石达4, 王芃1     
1. 中国地震局地震预测研究所, 北京 100036;
2. 中国地质大学(武汉)地球内部成像与探测实验室, 武汉 430074;
3. 浙江大学地球科学学院, 杭州 310027;
4. 华中科技大学物理学院基本物理量测量教育部重点实验室, 武汉 430074
摘要:在高频面波方法中,水平分辨率是指水平方向上分辨异常体的能力.异常体在水平方向上的长度可用水平方向上横波速度的异常尺度来确定.面波多道分析(MASW)方法被广泛应用于浅地表横波速度结构的探测,然而该方法确定的横波速度是整个检波器排列的平均计算结果,因此水平分辨率较差.另外,采用共中心点(CMP)多次覆盖的方式采集数据亦增加了野外的工作量.我们在MASW方法的基础上,应用面波层析成像方法,提出一套提高面波勘探水平分辨率的完整方法的技术流程.首先,利用波场分离技术获得准确的基阶或高阶模式面波,采用相位扫描的互相关方法测量多道面波记录中任意两道之间的面波走时;然后根据面波层析成像方法,获得高分辨率的各目标网格内的纯路径相速度频散曲线;最后反演所有目标网格内的纯路径相速度频散曲线,得到研究区域的拟二维横波速度结构.这套方法具有一定的抗噪能力,理论上它可以准确地提取相邻两道之间面波的相速度频散曲线;同时由于该方法最少只需要1个排列就可以获得拟二维横波速度结构,因此它显著减小了野外工作量.理论模型和实际资料都证实了这套方法可有效提高面波勘探的水平分辨率.
关键词: 水平分辨率      纯路径相速度频散曲线      面波层析成像方法      拟二维横波速度结构     
A travel-time tomography method for improving horizontal resolution of high-frequency surface-wave exploration
YIN XiaoFei1, XU HongRui2, XIA JiangHai3, SUN ShiDa4, WANG Peng1     
1. Institute of Earthquake Science, China Earthquake Administration, Beijing 100036, China;
2. Subsurface Imaging and Sensing Laboratory, China University of Geosciences, Wuhan 430074, China;
3. School of Earth Sciences, Zhejiang University, Hangzhou 310027, China;
4. MOE Key Laboratory of Fundamental Physical Quantities Measurements, School of Physics, Huazhong University of Science and Technology, Wuhan 430074, China
Abstract: In high-frequency surface-wave methods, horizontal resolution refers to the ability to distinguish anomalous objects that are laterally displaced from each other. The horizontal length of a recognizable geological anomalous body is measured by the scale of shear (S)-wave velocity anomalies. The multichannel analysis of surface wave (MASW) method is an efficient tool to determine near-surface S-wave velocity. However the inverted 1D S-wave velocity model is an averaged geophysical model under the receiver spread length, thus the MASW method has low horizontal resolution. In addition, the common middle point (CMP) roll-along acquisition method can increase the amount of fieldwork by the roll-along acquisition geometry. To solve these problems, on the basis of the MASW method, we propose a complete technique flow to improve horizontal resolution of surface-wave exploration by the travel-time tomography method. Firstly, we use the wave field separation technique to obtain accurate fundamental or higher mode of surface waves. And surface-wave phase-velocity dispersion curves between any two traces are calculated by the combined method of cross-correlation and phase-shift scanning from a multichannel record. Then with the application of travel-time tomography method, we obtain high resolution pure-path phase-velocity dispersion curves at diverse sizes of grids. Finally, the pseudo-2D S-wave velocity structure is reconstructed by inverting pure-path phase-velocity dispersion curves. The proposed method can effectively enhance the ability of random noise immunity, and can extract accurate surface-wave phase-velocity dispersion curves from a record with a short receiver spacing theoretically. Moreover, the pseudo-2D S-wave velocity structure can be obtained by the proposed method with at least only one seismic record, which reduces the amount of fieldwork dramatically. A synthetic test and a real-data example have demonstrated that the proposed method has a great potential in improving the horizontal resolution of surface-wave exploration using multichannel analysis.
Key words: Horizontal resolution    Pure-path phase-velocity dispersion curves    Travel-time tomography method of surface waves    A pseudo-2D S-wave velocity section    
0 引言

在高频面波方法中,水平分辨率是指水平方向上分辨异常体的能力(Yilmaz,1987).横波速度是浅地表地球物理研究中的基本参数,异常体在水平方向上的长度可用水平方向上横波速度的异常尺度来确定.

面波多道分析(Multichannel Analysis of Surface Waves,MASW)方法可用来获得地下横波速度结构,采用共中心点CMP(Common Middle Point)多次覆盖的方式采集面波数据,并且将反演得到的一维横波速度剖面按一定的顺序排列就可得到拟二维横波速度结构(Miller et al., 1999),因其能够准确地分离不同模式的面波以及便捷环保的野外数据采集方式,该方法在浅地表地球物理界引起了广泛的关注(Xia,2014Yin et al., 2014夏江海等,2015).然而如何提高面波勘探的水平分辨率是MASW方法自问世以来就存在的问题,这是因为MASW方法中面波相速度频散曲线反演获得的横波速度反映的是检波器排列(从第一个检波器到最后一个检波器的距离)的平均结果(Luo et al., 2009a),由该方法进行速度成像时“牺牲”了面波勘探的水平分辨率,滚动排列方式亦需要进行大量的野外数据的采集工作(Mi et al., 2017),同时,通过人工识别频散能量图提取频散曲线也不可避免地会产生人为误差.

Park(2005)指出面波勘探的水平分辨率主要受炮检距和检波器排列长度的限制,实际应用中水平分辨率约等于检波器的排列长度(Xia et al., 2005).为了提高面波勘探的水平分辨率,我们需要利用较短的检波器排列提取频散曲线(Xia and Miller, 2010),而较短的排列长度又会产生较大的频散曲线误差既而被引入拟二维横波速度结构的构建(Hayashi and Suzuki, 2004; Lin C P and Lin C H, 2007).假设水平方向的横波速度的变化较小,Luo等(2008a)利用相位扫描的互相关方法计算多道面波记录中间距较小的两道间的面波相速度频散曲线,再反演频散曲线获得拟二维横波速度结构.尽管这个方法在一定程度上可提高面波勘探的水平分辨率,然而由于该方法没有预先对面波记录进行模式分离,提取的面波相速度频散曲线可能受到体波或者其他模式面波的干扰(Boore,1969),因此该方法很难区分(相邻两道之间)尺度较小的异常体.多偏移距的相位分析(Multi-offset Phase Analysis,MOPA)方法利用最大数量的面波记录中两道间的相位信息探测水平方向上面波相速度频散曲线的变化,但是MOPA方法只适用于已知异常体的水平位置的情况(Vignoli et al., 2011).除此之外,Zeng等(2011)Pan等(2016)认为面波波形反演方法可直接获得浅地表的拟二维横波速度结构,该方法适用于含横向变化较大的速度模型,然而该方法的计算有时并不稳定,反演结果严重依赖于初始模型的选择,且计算成本(耗时)较高.

近些年来,面波走时层析成像的方法是探索浅地表地球结构的最有效的方法(Tondi et al., 2012Boué et al., 2014Xu et al., 2016).采用合适的网格剖分方式将研究区域分为不同大小的网格,可有效地提高速度模型的水平分辨率(Ajo-Franklin et al., 2006).

基于高频近似的假设,面波沿直射线传播(Um and Thurber, 1987),在利用高分辨率线性拉东变换(Luo et al., 2008b, 2009b)对面波记录进行准确的模式分离的前提下,本文提出利用面波走时层析成像的方法提高面波勘探的水平分辨率.该方法共包括四个步骤(图 1):(1)对于仅含某一特定模式面波的炮集记录(本文主要研究基阶模式面波),利用相位扫描的互相关方法提取相邻两道之间的面波相速度频散曲线;(2)根据走时、两道间的距离和面波相速度的关系,直接测量任意两道间的面波走时.通过检测板的测试将研究区域剖分为大小相等或不等的网格,然后把测量得到的面波走时作为输入数据,反演面波的走时获得不同频率时每个网格中的相速度分布;(3)既而得到每个网格的纯路径相速度频散曲线,反演纯路径相速度频散曲线得到每个目标网格的一维横波速度剖面;最后将研究区域的所有一维横波速度剖面按顺序排列,生成拟二维横波速度结构.面波走时层析成像的方法具有一定的抗噪能力,它可准确地提取炮集记录中任意相邻两道之间的纯路径相速度频散曲线.理论模型和实际资料都证实了面波走时层析成像的方法可有效提高面波勘探的水平分辨率.同时,最少只需要1个排列就可以获得排列下方的拟二维横波速度结构,相较于传统的CMP滚动排列的数据采集方式,该方法可显著减小野外工作量.另外,直接测量面波频散曲线避免了MASW方法中手动拾取频散曲线的人为误差.

图 1 利用面波走时层析成像方法获得拟二维横波速度结构的技术流程图 Fig. 1 Technique flow chart of a proposed surface-wave travel-time tomography method for generating a pseudo-2D S-wave velocity section
1 方法原理 1.1 相邻两道面波相速度频散曲线的计算

假设R1(t)和R2(t)是时间-空间(t-x)域的两道地震记录,它们的自相关功率谱可表示为

(1)

(2)

其中,R1(f)和R2(f)分别为R1(t)和R2(t)的傅里叶变换后的结果;R1*(f)和R2*(f)表示它们的复共轭函数;两道地震记录的互相关功率谱如下:

(3)

这里,Δϕ(f)表示某一频率f时所对应的两道地震记录的相位差.相速度VR(f)可用两道地震记录间的距离Δx和它们的相位差Δϕ(f)表示,

(4)

为了提高相邻两道的面波相速度频散曲线的计算精度,引入互相关函数评价频率f时相邻两道地震记录的数据质量和噪声的干扰程度(Guo and Liu, 1999),R1(t)和R2(t)的相干函数为

(5)

如果C(f)越大,则表明两道地震记录具有很好的相干性.在实际情况中,设定一个极限值0.8.若C(f)>0.8,则保留该频率所对应的面波相速度,反之,则应剔除这一频率及其对应的相速度.上述互相关方法可在任意两道地震记录之间提取频散曲线,同SASW方法类似,由于该方法仅用一对地震记录道提频散,表现为对噪声数据较为敏感.如果地震记录中包含噪声,将导致相速度频散曲线的计算存在误差(Al-Hunaidi,1993).因此频散计算前有必要对炮集记录进行预处理,仅保留某一特定模式面波的炮集记录,消除噪声和其他模式面波.

为了获得准确的相速度,我们将相位扫描(Park et al., 1998)引入互相关方法中提取相速度频散曲线.根据需要选择频率范围f∈(fmin, fmax),其中fminfmax分别为频率扫描的最小值和最大值.设置相速度扫描的范围VR∈(vmin, vmax),其中vminvmax分别为相速度扫描的最小值和最大值.根据Δϕ(f)=2πfΔx/VR(f),计算相位差Δϕ(f).最后,相速度VR在某一频率f时所对应的能量V(VR, f)可表示为

(6)

绘制频率-速度(f-v)域的频散能量图,找到每个频率f所对应的能量极大值点,自动拾取基阶模式面波的相速度.

1.2 面波走时层析成像的方法 1.2.1 目标网格的建立和走时的计算

地震层析成像中的速度模型通常是用剖分网格的方式确定的(Zhou,2003陈国金等,2006),如果剖分的网格太大会造成数据的“浪费”,如果剖分的网格太小则会加剧反演过程的多解性甚至在反演结果中引入虚假异常.因此,层析成像之前有必要对反演系统进行分辨率测试(Lévěque et al., 1993)以确定剖分网格的大小.由于不同网格中面波射线覆盖的密集程度不同,为了保证每个网格内有足够多的射线覆盖,本文采用变网格的方式(黄光南等,2013)将研究区域剖分为不同大小的网格,并假定每个网格内的速度为常数;经过大量的检测板测试实验,选择使输入的检测板模型恢复最好的剖分网格方式.

假设已知某一频率f时所对应的射线分布,经过检测板测试将研究区域剖分为n个不同大小的网格,共有m条射线穿过该区域,j=1, 2, …, m.每一条射线通过研究区域的总走时,等于该射线经过每个目标网格的走时的叠加.因此对于第j条射线,它通过研究区域的走时可表示为

(7)

式中,tj为第j条射线穿过研究区域的走时(观测数据);i=1, 2, …, n代表研究区域的所有剖分网格,si表示射线穿过第i个网格ai的慢度;dij为第j条射线在第i个网格中的射线路径长度.图 2a表示将研究区域网格化后,穿过研究区域的一条面波射线,从激发点S到接收点R的面波走时为射线在各个网格a1a2a3a4a5中走时的总和.图 2b表示第j条射线穿过第i个网格的细节图.

图 2 目标网格的剖分 (a)穿过研究区域的一条面波射线; (b)第j条射线穿过第i个网格的细节图. Fig. 2 Dividing the research area into different sizes of grids (a) A ray path passing through grid in research area; (b) Detailed description that jth path travels through the ith grid cell.

j条射线穿过研究区域的面波走时(观测数据)tj,等于第j条射线穿过每个网格的截距dij与其对应的网格的慢度si乘积的总和.(7)式可写成矩阵形式:

(8)

1.2.2 迭代反演

利用面波的直射线理论计算面波传播的走时是正问题;通过走时反演获得剖分网格的相速度分布是反问题.本文采用两轮迭代反演的方法计算目标网格的纯路径相速度频散曲线.对于一个层状地球模型而言,研究区域内面波射线的走时tobs是已知的,可表示为tobs=(tobs 1, tobs 2, …, tobs m)T.

以某一频率f时面波的走时反演为例,已知该频率所对应的射线分布,将每一个网格的慢度so作为初始模型,so=(so1, so2, …, son)T;由射线路径和初始模型计算的面波走时to可看作初始模型的响应,to=(to1, to2, …, tom)T,面波的走时与慢度满足线性关系:

(9)

其中,Δt=tobs-to表示一个m维的面波走时的观测值tobs与初始模型的响应to之间的误差向量;Δs代表一个n维的关于初始模型慢度so的修改向量,Δs=(s1-so1, s2-so2, …, sn-son)Td代表mn列的射线路径的矩阵(mn).

由于实际的观测数据中常含有噪声或其他干扰,为了避免较小的数据的改变引起的较大的模型参数的变化,同时考虑到浅地表剧烈的速度变化并不适用于长周期面波层析成像普遍应用的最光滑模型(例如:Yanovskaya and Ditmar, 1990Barmin et al., 2001),本文采用阻尼最小二乘的迭代反演方法(Xia et al., 1999)求解方程(9),通过多次修改慢度模型以获取最小长度解(王家映,2002).

目标函数可以定义为

(10)

这里,‖‖2代表向量的二范数;α为阻尼系数;W为加权矩阵;‖Δs22为模型长度;‖dΔst2WdΔst2代表数据的拟合误差.阻尼因子α控制Δs的变化方向及其收敛速度.如果反演模型的先验信息未知(例如:本文的研究中目标网格的慢度没有先验资料的情况),可利用L曲线法(Xia et al., 2012)寻找合适的阻尼因子使数据拟合差和模型长度最小.

根据奇异值分解技术最小化目标函数(10).将射线路径d进行奇异值分解,d=UΛVT,则慢度的修改向量可写作

(11)

每一次迭代反演计算中,对比反演结果正演计算的面波走时与观测走时,正演计算的走时同观测走时之间的均方根误差RMS满足

(12)

它们的绝对误差为

(13)

如果正演计算的走时与观测走时的残差大于给定的误差值,则修改慢度模型s=(Δs1+so1, Δs2+so2, …, Δsn+son)T,并将修改后的模型响应赋值于to,进行新一次的迭代反演计算;直到残差小于给定误差值或反演达到最大迭代次数时迭代反演终止,输出反演结果s,将慢度转换为相速度.

第一轮迭代反演计算结束时,我们可以获得频率f时对应的每一个目标网格的相速度.为了保证反演结果的精确性,把正演计算的走时与观测走时的绝对误差大于一倍均方根误差的观测走时剔除,然后将筛选出的观测走时进行第二轮迭代反演的计算(Yanovskaya and Ditmar, 1990),以第二轮迭代反演计算的各目标网格的相速度分布作为最终的反演结果.

1.3 拟二维横波速度结构的构建

通过面波走时层析成像的方法可以得到不同目标网格的纯路径相速度频散曲线.将每一个网格的纯路径相速度频散曲线作为输入数据,然后利用引入奇异值分解技术的阻尼最小二乘的反演方法修改模型(横波速度)以拟合测量的相速度频散曲线(Xia et al., 1999),直到反演达到最大迭代次数或频散数据的拟合差小于0.1 m·s-1时,迭代反演终止.最后,输出每个目标网格的一维横波速度剖面,并插值计算获得研究区域的拟二维横波速度结构.

2 纯路径相速度频散曲线精度的验证

本文选择一个四层模型(Zeng et al., 2011)(表 1)验证面波层析成像方法获得纯路径相速度频散曲线的抗噪能力.采用有限差分数值模拟的方法(Xu et al., 2007)合成炮集记录(图 3a),该排列共有50个检波器,道间距1 m,其中为了保证完整的面波记录,选择2 ms的采样率和1000 ms的记录时长,最小偏移距设为6 m(Xu et al., 2006);用高分辨率线性拉东变换将炮集记录从t-x域转换到f-v域(Luo et al., 2008b)(图 3b).本文对原始炮集记录(图 3a)加入N(0, 0.1)的高斯白噪(均值为0,方差为0.1的服从高斯分布的噪声)(图 3c),图 3d为与其对应的频散能量图,仅研究1~55 Hz频率范围的基阶模式面波的相速度频散曲线.由于面波的高频能量很弱,随机噪声对高频信号的干扰较大,因此仿照图 3b圈定出高阶模式面波的能量团范围,利用高分辨率线性拉东变换去除高阶模式面波能量团,得到仅包含基阶模式面波和随机噪声的t-x域的炮集记录(Luo et al., 2009b).对重构后的炮集记录进行相位扫描的互相关计算,得到相邻两道的面波相速度频散曲线(图 3e).

表 1 一个四层模型的参数 Table 1 Parameters of a four-layer earth model
图 3 加入白噪的炮集记录及其频散分析 (a)一个合成的四层模型的炮集记录; (b)利用高分辨率线性拉东变换获得频散能量图,挑选出高阶模式面波的频散能量团,黑色实心圆点表示用Knopoff方法计算的理论频散曲线,用红色实心圆点圈出高阶模式面波的频散能量团范围; (c)在炮集记录(a)中加入N(0, 0.1)的高斯白噪; (d)加噪后的炮集记录对应的频散能量图; (e)利用相位扫描的互相关方法计算任意两道之间的面波相速度频散曲线. Fig. 3 Analysis of noise contaminated shot data and their dispersion measurements (a) Synthetic shot gather of a four-layer model; (b) Selected higher-mode ranges in the dispersive-energy image. Solid black dots are the theoretical dispersion curves calculated by using Knopoff′s method, red dots represent the corresponding selected first-higher mode range; (c) The Gaussian distribution of N(0, 0.1) is added to the synthetic shot gather; (d) A dispersive-energy image of shot data (c) which has selected multimode ranges in the f-v domain. (e) Surface-wave phase-velocity dispersion curves calculated by any two traces with the noise contaminated data.

首先,将研究区域剖分为1 m网格大小的49个网格.由于地震排列的道间距为1 m,可使得每个目标网格与相邻两个地震记录道对应,例如第一个网格位于第一道和第二道之间,第二个网格位于第二道和第三道之间,以此类推.

然后,计算研究区域内任意两道之间的面波走时.利用L曲线法选择阻尼因子,将阻尼因子以500为间隔从500增加到15000,最终选择6000作为最佳阻尼因子(图 4a,由红色的矩形标注).对求得的面波走时按照合适的剖分网格进行层析成像,得到所有目标网格的纯路径相速度频散曲线(图 4b).将计算的相速度(相位扫描的互相关方法计算的相速度频散曲线或纯路径相速度频散曲线)同Knopoff方法计算的理论解析解对比.频率f(j)时,计算的相速度同理论解析解的相对误差为

(14)

图 4 面波层析成像方法的抗噪能力分析 (a)用L曲线法选择合适阻尼因子(红色的实心矩形标注最佳阻尼因子),其中用蓝色的空心圆形表示不同阻尼因子与数据拟合度和模型长度的关系; (b)利用面波层析成像方法获得的纯路径相速度频散曲线; (c)由相位扫描的互相关方法计算的任意两道之间的相速度频散曲线同理论解析解之间的相对误差“Mode-separated data” (用红色的实心圆形加实线标注),纯路径相速度频散曲线同理论解析解之间的相对误差“Inversion of travel-time” (用蓝色的实线标注); (d)在地震排列的不同水平位置上,相邻两道(间隔1 m)计算的相速度与理论解析解之间的相对误差,其中“Mode-separated data”代表用相位扫描的互相关方法计算的相速度频散曲线同理论解析解之间的相对误差,用蓝色的五角星加虚线标注;“Inversion of travel-time”代表纯路径相速度频散曲线同理论解析解之间的相对误差,用红色的矩形加虚线标注;“Residual error between the two”代表两者之间的绝对误差,用黑色的圆形加虚线标注; (e)面波射线的覆盖情况. Fig. 4 Analysis of anti-noise ability of travel-time tomography method (a) L-curve method used to determine the trade-off value of an optimal damping factor (a solid red square). Figures next to the symbols (hollow blue circles) are values of damping factors associated with the solutions; (b) Pure-path phase-velocity dispersion curves are calculated by travel-time tomography method; (c) Relative errors of dispersion curves calculated by consecutive two traces (labeled as "Mode-separated data") and the pure-path phase-velocity dispersion curves of the same two traces (labeled as "Inversion of travel-time"); (d) Relative errors of different lateral positions calculated by all consecutive two traces, where pentagrams with a dash line and rectangles with a dash line represent the relative errors associated with the dispersion curves of any two traces and pure-path phase-velocity dispersion curves; and circles with a dash line display their absolute errors are labeled as "Residual error between the two"; (e) Source-station surface-wave path distribution.

M个频率f(j)(j=1, 2, …, M),Vcal(j)和Vtheo(j)分别代表频率f(j)时计算的相速度和理论解析解.

图 4c中可见:(1)在整个频率范围内,直接利用相位扫描的互相关方法提取的相速度频散曲线的扰动较大,且同理论解析解的平均相对误差为10.30%;(2)层析成像方法获得的纯路径相速度频散曲线与理论解析解吻合较好,频散数据与理论解析解的平均相对误差为3.74%.研究表明面波层析成像方法可提高频散计算的精度,该方法具有一定的抗噪能力.

在检波器排列的不同水平位置上,相位扫描的互相关方法计算得到的相邻两道(间隔1 m)的相速度频散曲线的相对误差远大于相同位置上纯路径相速度频散曲线的相对误差(图 4d),前者的平均相对误差为9.41%,而后者的平均相对误差仅为2.77%.对于排列的中间区域(第5道到第46道之间),射线覆盖的密度大(图 4e),面波走时反演得到的纯路径相速度频散曲线与理论解析解的偏差较小.然而对于排列的边缘区域(第1道到第4道,第46道到第50道),由于射线分布较少(图 4e),尤其在第1道到第2道、第49道到第50道之间,仅有40条射线通过该研究区域,造成面波走时反演得到的纯路径相速度频散曲线与理论解析解的偏差较大.

因此我们认为面波层析成像方法可获得较精确的小台站间距的面波相速度频散曲线(平均相对误差小于3.0%);纯路径相速度频散曲线的精度主要与研究区域的射线覆盖情况有关:表现为中间区域有较高的精度,而边缘区域精度较低.

3 理论模型试验

一个三层的断层模型(图 5a):其中第一层的横波速度为200 m·s-1,层厚5 m;第二层为地表以下5~10 m的地层,包含两部分,第一部分(位于断层模型的左边,水平位置0~24 m)的横波速度为200 m·s-1,第二部分(位于断层模型的右边,水平位置24~49 m)的横波速度为400 m·s-1;半空间的横波速度等于400 m·s-1.根据有限差分数值模拟的方法合成炮集记录(图 5b),震源设置为主频为20 Hz的高斯子波(延迟时间60 ms),共50道地震记录;道间距1 m,最小偏移距为6 m,选择1 ms的采样率和1000 ms的记录时长;断层位于检波器排列的中点.

图 5 断层模型的炮集记录及其频散分析 (a)断层模型; (b)仅包含基阶模式面波的炮集记录; (c)利用相位扫描的互相关方法计算的任意两道之间的相速度频散曲线,由蓝色的空心圆形标注;用红色实线表示从高分辨率线性拉东变换的频散能量图中挑选的基阶模式面波的相速度频散曲线. Fig. 5 Analysis of shot gather and its dispersion measurements in a single step-shaped earth model (a) Illustration of a single step-shaped earth model; (b) Separated fundamental-mode shot data; (c) Phase-velocity dispersion curves calculated by any two traces; Blue hollow circles represent phase-velocity dispersion curves of any two traces and red line represents the fundamental-mode phase-velocity dispersion curve picked from the dispersive energy which is generated by high-resolution LRT method.

利用高分辨率线性拉东变换对原始炮集记录进行模式分离(Luo et al., 2009b),得到仅包含基阶模式面波的炮集记录.用相位扫描的互相关方法提取相速度频散曲线(图 5c),可见由相位扫描互相关方法提取的基阶模式面波的相速度频散曲线均分布在由高分辨率线性拉东变换拾取的相速度频散曲线的两侧.

研究区域的网格剖分是面波走时反演的基础,本文做了一系列的检测板测试.首先建立一个检测板的速度模型,将220 m·s-1的低速异常体和230 m·s-1的高速异常体呈棋盘状分布于总长度为50 m的一维测线上,且速度异常体的尺度大小为2 m.按照实际的射线分布计算面波走时,采用网格尺度1 m和0.5 m的网格剖分进行走时反演,分别研究不同频率时检测板的测试结果.

图 6a6f为频率8 Hz时检测板模型的测试结果.8 Hz时对应448条射线(448个台站对)(图 6d),沿测线1~40 m的范围内射线分布较为密集,边缘区域(40~50 m)的射线覆盖稀疏(仅有50条射线穿过).当网格尺度为1 m时,理论模型恢复较好(图 6b),残差在沿测线40~50 m的范围内最大(最大值为0.03 m·s-1),残差结果(图 6c)与射线分布(图 6d)有明显的对应关系;当网格尺度为0.5 m时,层析成像结果的误差较大,反演得到的速度和理论模型的最大偏差为0.12 m·s-1(图 6e6f).由此可以看出,面波层析成像方法中网格尺度越小即模型分辨率越高,然而层析成像结果并不是越好,实际应用中需要考虑模型分辨率和数据协方差之间的折衷.

图 6 断层模型中8 Hz时检测板的测试结果 (a)理论的检测板的速度模型; (b)以1 m大小的网格剖分得到的反演结果; (c)以1 m大小的网格剖分得到的反演结果与理论模型的绝对误差; (d)研究区域的面波射线的覆盖情况; (e)以0.5 m大小的网格剖分得到的反演结果; (f)以0.5 m大小的网格剖分得到的反演结果与理论模型的绝对误差. Fig. 6 Checkerboard tests of 8 Hz in a single step-shaped earth model (a) Input velocity on checkerboard grid; (b) Output velocity with grid spacing of 1 m, evaluated by all travel times with any two traces; (c) The residual error between the output velocity (b) and the input velocity (a); (d) Source-station surface-wave path distribution; (e) Output velocity with the grid spacing of 0.5 m evaluated by all travel times with any two traces; (f) The residual error between the output velocity (e) and the input velocity (a).

由于研究区域边缘的射线分布较少,通过不同网格大小的检测板测试实验表明,增加研究区域边缘的网格大小,可以避免由于射线分布稀疏造成的反演结果不稳定的问题.对断层模型的走时反演可采用变网格的方式,即将中间部分的网格大小设置为1 m,边缘区域网格大小设置为2 m.共有27个网格,具体的网格剖分方式如表 2.对实际求得的走时按照表 2的网格剖分进行走时反演,得到不同频率时所有目标网格的相速度分布(图 7).由图 7可见,在第25道的位置上存在一个较明显的速度界限,表现为第1道到第25道的范围内相速度较低(相速度的最大值为300 m·s-1),第26道到第50道的范围内相速度较高(最大值为400 m·s-1),相速度的变化很好地反映了断层的位置.

表 2 断层模型的网格剖分 Table 2 Grid description of a single step-shaped earth model
图 7 断层模型中不同网格的纯路径相速度频散曲线 Fig. 7 Pure-path phase-velocity dispersion curves of a single step-shaped model at different grid cells

整个研究区域总共有27个网格,对每一个网格的纯路径相速度频散曲线进行反演,得到每个目标网格的一维横波速度剖面.初始模型中每层厚度均取为1 m,纵波速度和密度同原模型一致.这里以第一个网格的纯路径相速度频散曲线的反演为例,如图 8a8b,对反演结果进行正演计算得到的相速度频散曲线与纯路径相速度频散曲线完全吻合,对应的横波速度结构的拟合效果也较好,表明反演结果稳定可靠.

图 8 断层模型的反演结果 (a)第一个网格的相速度频散曲线的拟合情况,其中“Measured”代表实测的相速度频散曲线,用蓝色的空心圆标注;“Initial”代表初始模型正演计算得到的相速度频散曲线,用粉色的虚线标注;“Final”代表反演得到的最终模型正演计算获得的频散曲线,用红色的实线标注; (b)第一个网格的横波速度的反演结果,其中“True model”代表真实模型,用蓝色的实线标注;“Initial”代表初始模型,用粉色的虚线标注;“Inverted”代表反演得到的横波速度,用红色的实心菱形加虚线标注; (c)断层模型的拟二维横波速度结构. Fig. 8 Inversion results in a single step-shaped earth model (a) The convergence of phase velocities in the first dissected grid cell. Measured phase-velocity dispersion curve labeled "Measured" is represented by blue hollow circles; those labeled "Initial" is calculated based on the initial S-wave velocity model; and those labeled "Final" is calculated based on the inverted S-wave velocity model; (b) The inverted S-wave velocity profile of the first dissected grid cell. The synthetic S-wave velocity model labeled "True model" is represented by a solid blue line. The initial S-wave velocity model labeled "Initial" is represented by a dash pink line; The inverted S-wave velocity model labeled "Inverted" is represented by red solid diamonds with a dash line; (c) The pseudo 2-D S-wave velocity section of a single step-shaped model.

将所有网格的一维横波速度剖面进行插值,得到研究区域的拟二维横波速度结构,如图 8c可见,断层位于第25道的位置上,第一部分(Trace<25)是表层层厚为10 m的两层模型;中间部分(Trace=25)位于垂直断层,上层和下层转换点的位置位于5 m和10 m;第三部分(Trace>25)是表层层厚为5 m的两层模型.反演结果与理论模型能够较好地吻合.

4 应用实例

1999年,美国堪萨斯大学浅地表地球物理团队对Olathe地区的基岩进行成像(Miller et al., 1999).观测系统如图 9,整个观测系统采用两对相互正交的测线,应用CMP多次覆盖的方式采集面波数据,采集过程中沿测线方向以1.2 m的间隔同时挪动震源和检波器排列,共13个地震排列.测线上的每一个炮集记录设置为:48个4.5 Hz的垂直分量的检波器以间隔0.6 m沿测线排列,最小偏移距为2.4 m,采样率为1 ms,记录时长1024 ms,用一个5 kg的锤子敲击边长为0.3 m的铁板作为震源.由于该研究区域的大部分地区被沥青层覆盖,故使用带有金属托盘的检波器采集地震数据.

图 9 Olathe地区的面波数据采集的测线布置情况(Miller et al., 1999) 共有4条测线,测线2(line 2)和测线4(line 4)平行,测线1 (line 1)和测线3(line 3)平行,黑色圆点表示钻井所在位置. Fig. 9 Layout of survey lines in Olathe, Kansas, U.S. (Miller et al., 1999) There are four lines, line 2 and line 4 are parallel to each other, while line 1 and line 3 are parallel to each other. Black dots represent the locations of boreholes.

选择测线1上的某一炮集记录(图 10a)进行研究.由于采集数据的过程中使用低频检波器和面波数据的采集装置,面波的能量较强;然而采集数据中仍含有高阶模式面波、体波和其他噪声干扰,因此有必要消除这些干扰.图 10b为由高分辨率线性拉东变换得到的频散能量图,在频率25~70 Hz的范围内出现两个能量团;且在频率70~100 Hz的高频范围内,两种能量团汇聚为一个能量团;这说明研究区域的表层横波速度基本一致.结合钻井测量的结果和当地的地质露头,圈定出基阶模式面波的相速度和频率范围(如图 10b中的黑色圆点圈出的范围).利用高分辨率线性拉东变换(Luo et al., 2009b)去除高阶模式面波和其他干扰波,得到仅包含基阶模式面波的炮集记录和频散能量图像(图 10c10d).

图 10 Olathe地区采集实例中的炮集记录及其频散分析 (a) Olathe地区采集的48道的炮集记录; (b)原始炮集记录所对应的面波频散能量图像,并用黑色实心圆点圈出基阶模式面波的频散能量团范围; (c)仅包含基阶模式面波的炮集记录; (d)基阶模式面波的炮集记录所对应的面波频散能量图. Fig. 10 Analysis of shot gather and dispersion in example of Olathe, Kansas, U. S. (a) Data acquired in Olathe, Kansas, U.S., with a 48-channel system; (b) A dispersive image in the f-v domain of shot gather (a) and the selected fundamental mode in black dots; (c) The shot gather only containing fundamental mode surface waves; (d) A dispersive image of shot gather (c).

对于只包含基阶模式面波的炮集记录,利用相位扫描的互相关方法计算任意两道之间的相速度频散曲线.为了保证研究区域的射线分布均匀,进行大量的检测板测试,选择最合适的网格剖分,如表 3.根据目标网格的剖分,计算研究区域任意两道之间的面波走时;并结合L曲线法选取合适的阻尼因子,将阻尼因子以50为间隔从100增加到4000,最终选择2000作为最佳阻尼因子(图 11,由红色的矩形标注).然后对求得的走时按照剖分的目标网格进行层析成像,得到每一个目标网格的纯路径相速度频散曲线.

表 3 Olathe地区测线1上的某一个炮集记录的网格剖分 Table 3 Grid description of a shot gather of survey line 1 in Olathe, Kansas, U. S.
图 11 Olathe地区采集实例中的阻尼因子的选取 用L曲线法选择合适阻尼因子(红色的实心矩形标注最佳阻尼因子),其中用蓝色的空心圆形表示不同阻尼因子与数据拟合度和模型长度的关系. Fig. 11 The selection of an optimal damping factor in example of Olathe, Kansas, U.S. The L-curve is used to determine the trade-off value of an optimal damping factor (a solid red square). Figures next to the symbols (hollow blue circles) are values of damping factors associated with the solutions.

最后,反演每一个目标网格的纯路径相速度频散曲线,获得每个网格的一维横波速度剖面.仿照Miller等(1999)在这个研究区域关于MASW方法的初始模型的布设方式:假设共有十层模型,半空间以上九层地层的层厚从表层向下分别为0.3,0.39,0.6,0.99,1.11,1.23,1.38,1.5,1.5 m;由Xia等(1999)的公式计算初始横波速度;密度和纵波速度都是根据钻井资料而定,且在反演过程中保持不变;各层的密度设为2000 kg·m-3;纵波速度从表层到半空间分别为651,731,808,908,1063,1249,1442,1651,2503,2700 m·s-1.图 12a为所有目标网格的一维横波速度剖面插值后得到的拟二维横波速度结构.由图 12a可见,距表层4 m以下地层为基岩(这与钻井的结果一致);从第1098道到第1130道(地下4~8 m)的范围,横波速度较大,说明该区域的基岩坚固且不易断裂,根据地质露头信息推断该区域的基岩以石灰岩为主;从第1134道到第1145道(地下4~8 m)的范围,横波速度较低,先验地质信息表明该研究区域的基岩为页岩(岩石的硬度较低,且易断裂).

图 12 Olathe地区采集实例中面波层析成像方法与MASW方法的结果的对比 (a)仅利用1个排列基于面波层析成像的方法反演纯路径相速度频散曲线得到的拟二维横波速度结构; (b)利用13个地震排列基于MASW方法得到的拟二维横波速度结构; (c)利用面波层析成像得到的横波速度结构反推从第1108道到第1124道的MASW方法的结果. Fig. 12 Comparison between surface-wave travel-time tomography method and the MASW method (a) The pseudo-2D S-wave velocity section generated by inverting a series of pure-path phase-velocity dispersion curves computed by the travel-time tomography method with only one seismic-spread; (b) The pseudo-2D S-wave velocity section generated by the MASW method with 13 seismic-spreads; (c) The MASW result from station 1108 to station 1124 derived from the surface-wave travel-time tomography.

引入关于该区域的MASW的结果(图 12b),通过对比发现:(1)面波层析成像方法获得的拟二维横波速度结构同MASW方法的结果均显示表层4 m以下为基岩;从第1098道到第1130道的基岩对应的横波速度较大,而从第1134道到第1145道的范围对应较小的横波速度.(2)两者的结果存在不同点:从第1126道到第1134道的范围,层析成像的结果表明该区域的基岩为一个高速块体,而MASW方法没有显示相同区域中异常块体的存在;在地下4 m以上的地层,面波层析成像方法获得的横波速度结构比MASW方法的结果可反映更多水平方向上横波速度变化的信息.

为了进一步验证面波层析成像方法计算的可靠性,用面波层析成像方法计算得到的横波速度结构反推MASW方法的结果(图 12c).将从第1098道到第1118道共20个网格的纯路径相速度频散曲线的平均值作为第1108道所在网格的纯路径相速度频散曲线,以此类推,分别得到从第1108道到第1124道所在网格的纯路径相速度频散曲线,并进行反演获得拟二维横波速度结构.结果表明,由层析成像的结果通过平均效应反推得到的拟二维横波速度结构与MASW方法得到的拟二维横波速度结构高度吻合,证明MASW方法得到的拟二维横波速度结构是一个平均化的结果,层析成像方法得到的拟二维横波速度结构比MASW方法的结果的水平分辨率更高.

5 结论

本文提出利用面波层析成像的方法提高面波勘探的水平分辨率,该方法技术的研究包括特定模式面波的相速度频散曲线的提取、面波层析成像方法的研究、纯路径相速度频散曲线精度的测试、理论和实际数据的应用.

由于高频面波勘探的水平分辨率与检波器排列的间隔有关,如何在较小的两道间得到准确的相速度频散曲线是一个挑战.本文在层析成像的方法中使用两轮迭代反演计算,可有效地去除测量不准的走时,提高反演精度.

将高斯随机噪声加入理论模型,验证纯路径相速度频散曲线的精度.研究表明:(1)面波层析成像方法具有一定的抗噪能力;(2)纯路径相速度频散曲线的精度与射线覆盖的密集程度有关,表现为射线覆盖密集区域的纯路径相速度频散曲线计算的精度高,而射线覆盖稀疏区域的精度较低.由于浅层面波勘探基于一维观测系统,所以本文的方法能够很好地分辨尺度大于或等于相邻两道之间的异常.

断层模型为浅地表地球物理研究中常见的地质体的简化,将面波层析成像方法应用于一个断层模型的研究,获得的拟二维横波速度结构可很好地刻画该模型在水平方向的速度改变.在Olathe基岩成像的实例研究中,将该方法同MASW方法的结果进行比较,可以看到:仅利用1个排列基于面波层析成像技术计算得到的拟二维横波速度结构与利用13个地震排列基于MASW方法的结果相似,并且前者刻画的细节更多.所以该方法较之于MASW方法的水平分辨率明显提高,同时野外工作量显著减小.

References
Ajo-Franklin J B, Urban J A, Harris J M. 2006. Using resolution-constrained adaptive meshes for traveltime tomography. Journal of Seismic Exploration, 14: 371-392.
Al-Hunaidi M O. 1993. Insights on the SASW nondestructive testing method. Canadian Journal of Civil Engineering, 20(6): 940-950. DOI:10.1139/l93-126
Barmin M P, Ritzwoller M H, Levshin A L. 2001. A fast and reliable method for surface wave tomography. Pure and Applied Geophysics, 158(8): 1351-1375. DOI:10.1007/PL00001225
Boore D M. 1969. Effect of higher mode contamination on measured love wave phase velocities. Journal of Geophysical Research, 74(27): 6612-6616. DOI:10.1029/JB074i027p06612
Boué P, Roux P, Campillo M, et al. 2014. Phase velocity tomography of surface waves using ambient noise cross correlation and array processing. Journal of Geophysical Research:Solid Earth, 119(1): 519-529. DOI:10.1002/2013JB010446
Chen G J, Cao H, Wu Y S, et al. 2006. Effects of velocity contrast on the quality of crosswell traveltime tomography and an improved method. Chinese Journal of Geophysics, 49(3): 915-922. DOI:10.3321/j.issn:0001-5733.2006.03.038
Guo T S, Liu L B. 1999. Non-intrusive evaluation of submarine tunnel foundation using dynamic high-frequency surface wave prospecting. //Proceedings of the Symposium on the Application of Geophysics to Engineering and Environmental Problems (SAGEEP 1999). Oakland, California, 67-74.
Hayashi K, Suzuki H. 2004. CMP cross-correlation analysis of multi-channel surface-wave data. Exploration Geophysics, 35(1): 7-13. DOI:10.1071/EG04007
Huang G N, Liu Y, Tryggvason A, et al. 2013. Variable grid spacing velocity tomography. Oil Geophysical Prospecting, 48(3): 379-389.
Lévěque J J, Rivera L, Wittlinger G. 1993. On the use of the checker-board test to assess the resolution of tomographic inversions. Geophysical Journal International, 115(1): 313-318. DOI:10.1111/gji.1993.115.issue-1
Lin C P, Lin C H. 2007. Effect of lateral heterogeneity on surface wave testing:Numerical simulations and a countermeasure. Soil Dynamics and Earthquake Engineering, 27(6): 541-552. DOI:10.1016/j.soildyn.2006.10.008
Luo Y H, Xia J H, Liu J P, et al. 2008a. Generation of a pseudo-2D shear-wave velocity section by inversion of a series of 1D dispersion curves. Journal of Applied Geophysics, 64(3-4): 115-124. DOI:10.1016/j.jappgeo.2008.01.003
Luo Y H, Xia J H, Miller R D, et al. 2008b. Rayleigh-wave dispersive energy imaging using a high-resolution linear Radon transform. Pure and Applied Geophysics, 165(5): 903-922. DOI:10.1007/s00024-008-0338-4
Luo Y H, Xia J H, Liu J P, et al. 2009a. Research on the middle-of-receiver-spread assumption of the MASW method. Soil Dynamics and Earthquake Engineering, 29(1): 71-79. DOI:10.1016/j.soildyn.2008.01.009
Luo Y H, Xia J H, Miller R D, et al. 2009b. Rayleigh-wave mode separation by high-resolution linear Radon transform. Geophysical Journal International, 179(1): 254-264. DOI:10.1111/gji.2009.179.issue-1
Mi B B, Xia J H, Shen C, et al. 2017. Horizontal resolution of multichannel analysis of surface waves. Geophysics, 82(3): EN51-EN66. DOI:10.1190/geo2016-0202.1
Miller R D, Xia J H, Park C B, et al. 1999. Multichannel analysis of surface waves to map bedrock. The Leading Edge, 18(12): 1392-1396. DOI:10.1190/1.1438226
Pan Y D, Xia J H, Xu Y X, et al. 2016. Love-wave waveform inversion in time domain for shallow shear-wave velocity. Geophysics, 81(1): R1-R14. DOI:10.1190/geo2014-0225.1
Park C B, Miller R D, Xia J H. 1998. Imaging dispersion curves of surface waves on multi-channel record. //SEG Technical Program Expanded Abstracts 1998. New Orleans: SEG, 1377-1380.
Park C B. 2005. MASW horizontal resolution in 2D shear-velocity (Vs) mapping. Kansas Geological Survey Open-file Report, 4: 2005.
Tondi R, Cavazzoni C, Danecek P, et al. 2012. Parallel, 'large', dense matrix problems:Application to 3D sequential integrated inversion of seismological and gravity data. Computers & Geosciences, 48: 143-156.
Um J, Thurber C. 1987. A fast algorithm for two-point seismic ray tracing. Bulletin of the Seismological Society of America, 77(3): 972-986.
Vignoli G, Strobbia C, Cassiani G, et al. 2011. Statistical multioffset phase analysis for surface-wave processing in laterally varying media. Geophysics, 76(2): U1-U11. DOI:10.1190/1.3542076
Wang J Y. 2002. Inverse Theory in Geophysics. 2nd ed. Beijing: Higher Education Press.
Xia J H, Miller R D, Park C B. 1999. Estimation of near-surface shear-wave velocity by inversion of Rayleigh waves. Geophysics, 64(3): 691-700. DOI:10.1190/1.1444578
Xia J H, Chen C, Tian G, et al. 2005. Resolution of high-frequency Rayleigh-wave data. Journal of Environmental & Engineering Geophysics, 10(2): 99-110.
Xia J H, Miller R D. 2010. Estimation of near-surface shear-wave velocity and quality factor by inversion of high-frequency Rayleigh waves. //Advances in Near-surface Seismology and Ground-Penetrating Radar. Tulsa, OK: Society of Exploration Geophysicists, 17-36.
Xia J H, Xu Y X, Miller R D, et al. 2012. Estimation of near-surface quality factors by constrained inversion of Rayleigh-wave attenuation coefficients. Journal of Applied Geophysics, 82: 137-144. DOI:10.1016/j.jappgeo.2012.03.003
Xia J H. 2014. Estimation of near-surface shear-wave velocities and quality factors using multichannel analysis of surface-wave methods. Journal of Applied Geophysics, 103: 140-151. DOI:10.1016/j.jappgeo.2014.01.016
Xia J H, Gao L L, Pan Y D, et al. 2015. New findings in high-frequency surface wave method. Chinese Journal of Geophysics, 58(8): 2591-2605. DOI:10.6038/cjg20150801
Xu H R, Luo Y H, Chen C, et al. 2016. 3D shallow structures in the Baogutu area, Karamay, determined by eikonal tomography of short-period ambient noise surface waves. Journal of Applied Geophysics, 129: 101-110. DOI:10.1016/j.jappgeo.2016.03.037
Xu Y X, Xia J H, Miller R D. 2006. Quantitative estimation of minimum offset for multichannel surface-wave survey with actively exciting source. Journal of Applied Geophysics, 59(2): 117-125. DOI:10.1016/j.jappgeo.2005.08.002
Xu Y X, Xia J H, Miller R D. 2007. Numerical investigation of implementation of air-earth boundary by acoustic-elastic boundary approach. Geophysics, 72(5): SM147-SM153. DOI:10.1190/1.2753831
Yanovskaya T B, Ditmar P G. 1990. Smoothness criteria in surface wave tomography. Geophysical Journal International, 102(1): 63-72. DOI:10.1111/gji.1990.102.issue-1
Yilmaz Ö. 1987. Seismic Data Processing. Tulsa, OK: Society of Exploration Geophysicists, 526.
Yin X F, Xia J H, Shen C, et al. 2014. Comparative analysis on penetrating depth of high-frequency Rayleigh and Love waves. Journal of Applied Geophysics, 111: 86-94. DOI:10.1016/j.jappgeo.2014.09.022
Zeng C, Xia J H, Miller R D, et al. 2011. Feasibility of waveform inversion of Rayleigh waves for shallow shear-wave velocity using genetic algorithm. Journal of Applied Geophysics, 75(4): 648-685. DOI:10.1016/j.jappgeo.2011.09.028
Zhou H W. 2003. Multiscale traveltime tomography. Geophysics, 68(5): 1639-1649. DOI:10.1190/1.1620638
陈国金, 曹辉, 吴永栓, 等. 2006. 速度差对井间走时层析成像质量的影响及其改进方法. 地球物理学报, 49(3): 915-922. DOI:10.3321/j.issn:0001-5733.2006.03.038
黄光南, 刘洋, TryggvasonA, 等. 2013. 变网格间距速度层析成像方法. 石油地球物理勘探, 48(3): 379-389.
王家映. 2002. 地球物理反演理论. 2版. 北京: 高等敎育出版社.
夏江海, 高伶俐, 潘雨迪, 等. 2015. 高频面波方法的若干新进展. 地球物理学报, 58(8): 2591-2605. DOI:10.6038/cjg20150801