传统的瑞利波成像方法主要基于主动源,常采用炸药或重锤激发,炸药具有较强的破坏性,在城镇人口密集地区无法应用。锤击震源能量有限,主要集中在几到几十赫兹的频带范围,对深部信息探测不足,且极易受城市活动干扰。近年来,迅速发展起来的被动源瑞利波勘探技术将环境噪声作为信号源,具有抗干扰能力强、施工条件受限少等优点,更适合于城市周边等地区的勘探工作。因此,被动源瑞利波勘探技术受到了越来越多的关注。

1957年,Aki^[1]提出利用背景噪声研究地下结构,从台阵信号中提取出瑞利波频散信息,这种方法被称为空间自相关法(SPAC)。1968年,Claerbout^[2]从理论上指出对经过两个台站的透射波场进行互相关计算,可以获得台站间介质的反射响应函数。2003年,Campillo与Paul^[3]根据Claerbout的理论,对不同台站的地震尾波做互相关计算,提取出台站间的格林函数。在此之后,2004年,Shapiro等人^[4]成功利用背景噪声提取出了格林函数,并得到了美国加州地区的瑞利波速度结构图像。2007年,Bensen等人^[5]对被动源噪声成像方法的流程进行了总结,探讨了从经验格林函数中提取面波速度频散曲线的问题。近年来,被动源瑞利波勘探方法在实际应用方面取得了越来越多的研究成果。2008年,Brenguier等人^[6]研究了Foumaise火山地区内的三维速度结构。2009年,房立华等人^[7]在华北地区运用噪声互相关技术得到了该地区的短周期面波群速度。2011年,Ozalaybey S等人^[8]利用微动数据横—纵分量谱比法(HVSR)调查了土耳其Izmit Bay区域的盆地结构。2012年,郑现^[9]通过地震噪声信号获取了中国中东部地区的群速度结构。2013年,徐佩芬等人^[10]利用空间自相关法从地震台站微动信号的垂直分量中提取瑞雷波相速度频散曲线,反演得到了首都圈地区30 km深度内的S波速度变化特征。2016年,张宝龙等人^[11]在五大连池火山区布设微动台阵,利用扩展空间自相关法(ESPAC)提取频散曲线结合面波层析成像技术获得了700 m以浅地层的S波速度结构信息。上述研究通常被用于大尺度地下结构探测,观测台站之间的距离可达上千米,数据采集通常需要持续几个月甚至一年。由于被动源信号的信噪比较低,信号中有效波的能量较弱,频散曲线提取时,需要对不同时段的噪声数据进行叠加,以便更好地突显出有效信号。

近几年,被动源探勘方法在小尺度结构探测的应用也逐渐受到关注,如2009~2012年,徐佩芬等^[12-14]将被动源勘探方法分别应用到煤矿陷落柱、地铁孤石、地热探测中。2016年,吴学明等人^[15]将被动源空间自相关法和折射微动法(ReMi法)应用到浅层水电工程勘探中。2017年,翟法智等人^[16]利用折射微动探测方法进行了轨道交通暗浜勘查。与大尺度深部结构探测相比,对于小尺度浅层实际工程探测,人们总希望数据观测周期越短越好(如一天或几个小时)。但短时间的背景噪声数据会大大限制互相关计算过程中去除干扰的效果,在利用Aki法拾取频谱实部曲线零点时,互相关函数频谱会出现零点增多或缺失的现象,给频散曲线的提取带来了困难。因此,亟需采用适当的质量控制措施来确保频散曲线的提取精度。

针对这一问题,笔者拟在数据预处理、时间域和频率域归一化、时窗长度选择、互相关信号滤波以及互相关函数频谱零点选取等各个环节进行改进,形成一整套频散曲线质量控制流程。然后,通过理论模型被动源数值模拟结果以及陕西凤翔县野外实际噪声数据处理,进一步验证本文频散曲线质量控制方法的可行性和有效性。

2003年,Campillo与Paul^[3]对不同台站的地震尾波记录进行互相关并叠加,发现得到的波形与理论模型合成的格林函数一致。2008年,Gouedard等人^[17]总结前人的研究成果,得出两点间的格林函数可由两点间互相关函数的导数近似表示。研究表明,由于噪声源不均匀随机分布,实际计算时所选取的检波器距离至少应大于3倍波长。当地层横波速度越大时,所应选取的检波器距离也应越大^[5]。

Aki^[1]指出在噪声场均匀分布的情况下,空间自相关系数可由时间域的互相关谱代替。Tasi和Moschetti等人^[18]证明了空间自相关和互相关函数之间的等价关系。2009年,Ekstrom等人^[19]通过实际数据证明了可以将Aki公式与互相关法相结合,来提取被动源瑞利波的频散曲线。

Aki^[1]指出,空间自相关系数与第一类零阶贝塞尔函数的关系为:

式中,ω表示频率,ρ(r,ω)定义为空间自相关系数。J₀表示第一类零阶贝塞尔函数,r为台站间的距离,c(ω₀)表示频率ω₀处的相速度值。张宝龙等^[11]指出在长时间观测条件下,两台站的相关结果可等效于台站间的方位平均。因此式(1)表明如果得到了(r,ω₀)对应的相关系数ρ(r,ω₀),那么就可以求得对应ω₀频率处的相速度值c(ω₀)。

在计算得到的频谱中,有时会因为环境的干扰造成频谱上零点数量的变化,在某一小段频率上出现多个零点,使频谱实部曲线上的零点无法很好地与贝塞尔函数的零点对应,同时由于在低频处出现了较多个多余的零点,导致后续有效波的频散零点无法完全表现出来。李欣欣^[20]的研究表明,可以加入一个调控参数m,来改变频谱实部曲线零点与贝塞尔函数零点的对应关系,以此来调控频散曲线的形态,得到相速度计算公式为:

式中:f_n表示互相关谱上第n个零点,Z_n表示贝塞尔函数的第n个根,r表示两检波器的距离,v(ω_n)表示面波相速度,m=±1,2,…表示增加或缺失的零点个数。在计算相速度时,首先对互相关结果做傅里叶变换,求得实部频谱曲线,再利用三次多项式插值对实部频谱曲线的零点进行拾取,最后将互相关频谱的零点与第一类零阶贝塞尔函数的零点代入式(2),通过不断改变两种零点的对应关系来对频散曲线的形态进行调整,然后根据地质资料确定相速度值的变化范围选出拟合效果最好的频散曲线。

文中分别采用一组被动源数值模拟数据和一组野外实际观测数据来验证被动源瑞利波两道法提取频散曲线质量控制措施的有效性,野外实际数据采集时分别进行了主动源观测和被动源观测。

Lawrence等人^[21]的研究表明,由噪声源激发的高频面波抵达检波器时已经可近似为稳定的平面波,因此被动源信号可以看作由多个不同位置、不同深度、不同频率以及不同激发时间的主动源信号的叠加。文中采用一个两层模型进行被动源波场模拟,模型参数如表1所示。地层模型大小为200 m×100 m,波场模拟采用交错网格有限差分算法^[22],网格间距为1 m,时间步长为0.1 ms。模拟过程中,需保证每个时段内有足够数量的噪声源(通常不小于16个)。文中设置了多组不同位置、不同深度以及不同频率的人工震源,在不同时刻进行激发,在相同时窗内分别计算它们的地震记录,最后对相同检波器处的地震记录进行叠加,就可以看作该台站(检波点)位置处的被动源信号。震源分布情况如图1所示,采用雷克子波,依次激发(激发时设置不同的延迟参数来满足时间上的随机性)。地震记录的采样率为0.1 ms,采样点数为5 000个,检波器起始位置为0 m,终止位置为200 m,道间距为2 m,共101道。合成地震记录如图2所示,图中显示的是每隔10道的地震记录。

野外背景噪声实测试验在陕西省凤翔县西南庞家务村附近进行,采用14台DZS-1数字地震仪进行背景噪声数据采集。测线排列方式如图3,台站间距为50 m,采样间隔为10 ms,每条测线连续记录10个小时的背景噪声。同时,在原有的测线上进行了一次主动源瑞利波探测,以便将主动源频散曲线与被动源频散曲线进行对比。图4a为测线6上的野外背景噪声信号(0~300 s),图4b为测线6上第1个台站位置处对应的主动源面波记录。

两道法提取频散曲线的流程如图5所示,其中每一环节所采用的处理方法和计算参数都会影响频散曲线的提取效果。所有处理步骤的最终目的是尽量避免互相关函数频谱上出现零点增多或缺失的问题,得到与实际情况尽量符合的频散曲线。为评价实测频散曲线与理论频散曲线的接近程度,文中引入均方误差和相关系数来评价频散曲线的提取效果。均方误差是观测值与真实值之差平方的期望值,均方误差的值越小,频散曲线的提取精度越好。相关系数是研究变量之间线性相关程度的量,相关系数的绝对值越接近1,说明两变量的相关性越好。

不同型号的仪器在出厂时都会提供准确的响应参数,但是随着仪器的反复使用,可能会出现老化或者其他偏差。当仪器的响应参数不同时,采集到的数据将不具有可比性,因此需要先将这种影响消除。去除仪器响应的方法通常采用SAC软件中的transfer命令,常用的仪器响应文件有两种格式,即 RESP 文件和 PZ 文件。从易用性来看,处理RESP 文件过程简单,只需要采用一个evalresp命令语句,即可同时对所有波形数据进行处理。而处理PZ文件则需要人为的从数据文件的文件名或头段中提取信息,指定对应的响应文件,这种方法操作起来比较繁琐,但数据处理速度更快。在本次野外测量试验之前通过对仪器进行检验,仪器工作状态良好,因此这里可以略去此步骤。

去趋势处理可以消除由于检波器在采集数据时出现零漂现象而对计算产生的影响。通常去除趋势的方法为趋势拟合和希尔伯特黄变换,原理是将数据减去一条最优(最小二乘)拟合直线,使去趋势后的数据均值为零。本实验采用MATLAB中的detrend命令进行处理,设置语句参数“linear”、“constant”即可去除地震数据中的趋势和均值。

被动源噪声信号的频率通常都比较低,而一些突发事件(如汽车经过、工业生产活动)会产生一些高频信号,导致在信号频谱上出现高频峰值,影响对同时间段内有效信号的提取工作,因此需要窄通滤波器对数据进行滤波,消去高频噪声。在本实验中设置一个巴特沃斯带通滤波器,通带频率为0.2~35 Hz。图6所示为野外实测背景噪声信号的滤波结果,图6a为野外实测数据测线6第一个台站采集的地震信号频谱,在40~50 Hz之间存在有部分高频噪声干扰,而被动源噪声信号主要频率分布在0~20 Hz范围内,其中在5、12 Hz处出现峰值;滤波后的频谱如图6b所示,可以看出40~50 Hz范围内的高频噪声得到了压制,同时保持原有频率分布趋势不变。

时间域归一化在数据预处理中起到十分关键的作用,其目的是消除天然地震以及采集过程中出现的瞬时信号的干扰。常用的方法有“one-bit”法和滑动绝对平均法。“one-bit”法的计算原理是将原始信号中的正值用1来代替,负值用-1代替。滑动绝对平均方法的计算原理是选取一个时间窗口,计算该窗口内所有振幅的平均值,将平均值的倒数作为权系数与该段数据的中心点相乘得到新的地震数据。之后逐段移动时窗,直至整条地震序列的最后一位,便可得到归一化之后的地震数据。滑动绝对平均方法可以用公式表示为:

式(3)中,d_j为经过前几步处理后的地震信号,n表示时窗函数的中心位置,W_n表示每个数据点的权重,d_n和 d˙n分别表示滑动绝对平均处理前后的地震记录。“one-bit”法计算方便,操作简单,但是由于改变了原有振幅值,可能会减弱频谱中部分有效信息。而滑动绝对平均法则可以很好地保留原始噪声信号的频谱信息,缺点是计算过程繁琐复杂,花费时间长。此外对于不同的时窗长度,信号的归一化程度也会不同。

仍选用野外实测数据测线6第一个检波器采集的地震信号,对其分别运用两种方法进行处理,计算结果如图7所示,其中图7c和图7d是采用“one-bit”法处理得到的地震信号及其频谱,具体做法是将振幅大于零的值设置为1,将振幅小于零点值设置为-1。时间域归一化处理对信号的频率成分和频谱分布趋势不会造成大的影响,同时该方法对信号的振幅大小进行了调整,可以最大程度地压制信号的瞬时“突跳”。图7e和图7f是采用滑动绝对平均法处理得到的地震信号及其频谱,虽然滑动绝对平均法也可以保留原有地震信号的频谱特征,但是对比图7a和图7e可以看出,处理后的地震信号中仍然存信号“突跳”现象(11 s附近),对时间域信号的压制效果不如“one-bit”法, 即在进行时间域归一化时,“one-bit”法处理效果更好。

频率域归一化也叫做频谱的白化处理。在一般的频谱曲线上,通常会有一些频率高峰的存在,而通过白化处理,可以降低某一单独频率的幅值,使得整个频谱曲线更为均匀,避免在之后的计算中出现“谱孔”现象。频率域归一化可以拓宽背景噪声信号的频带,对能量较弱的某些频段内的波进行增强,使得处理后的信号包含有各个频率成分。

具体做法是对经过时间域归一化后的野外实测噪声信号(图7c)进行傅里叶变换得到其频谱数据(图7d),然后分别采用“one-bit”法和滑动绝对平均法对频谱数据进行处理,结果如图8所示。其中图8a和图8b分别是运用“one-bit”法处理后得到的地震信号和频谱,由于“one-bit”法是将频谱大于0的值统一设定为1,而信号的频谱值均大于0,因此信号频谱变为一条数值为1的直线,从而保证各个频段内波的能量均衡,达到拓宽频带的目的。而滑动绝对平均法可以通过改变窗口长度来控制归一化效果,在保留原有频率分布情况下拓宽信号的频带宽度,如图8d在5 Hz和10~15 Hz处的频谱变化趋势得以保留。运用滑动绝对平均法时,窗口长度需要根据噪声信号的频率和振幅特征进行判定,假设出现幅值高峰的频率范围内有N个数据点(如图8d中红框所示),则窗口可以选择为包含2N或3N个点的数据段,这是为了保证在滑动过程中,有足够数量的非高峰幅值与高峰幅值做平均,达到频率归一化的目的。

两种方法的处理效果可通过考察互相关结果的信噪比来衡量,信噪比计算方法为:

式(4)表示将互相关函数中面波最大到时之后的部分作为噪声信号,其标准差作为分母,将互相关函数中最大振幅值作为分子,两者相除即为互相关函数的信噪比。将归一化后的数据通过傅里叶反变换转换为时间域信号(即图8a和图8c),采用间隔5道的两地震记录(即第1道与第6道、第2道与第7道等)通过互相关运算提取瑞利面波信号,得到多个互相关结果的信噪比,如表2所示。可以看出采用滑动绝对平均法处理后的信噪比均高于one-bit法处理后的信噪比,这是由于在采用one-bit法进行归一化时,与滑动平均法相比过多地抬高了高频信号的幅值,导致了信噪比降低。因此,我们将滑动绝对平均法作为谱白化过程中的首选方法。

在进行互相关计算时,首先需进行数据分段,然后将不同时段的互相关结果进行叠加计算。当记录时长一定时,分段时窗长度越长,面波发育完善,单次互相关效果越好,但叠加次数也会越少。而分段时窗长度越短,面波信息可能无法被完全包含在时窗内,单次互相关效果越差,但叠加次数越多。因此如何平衡好时窗长度与叠加次数就显得尤为重要。

时窗长度的选取原则是在确保面波信息完整的情况下尽可能多地增加叠加次数,从而提高互相关结果的可靠性。根据噪声源均分理论,两台站间距大于三倍波长时,互相关计算才能可靠地恢复对应波长的面波信号,定义一个地层速度窗口(v_min,v_max),根据最远两台站之间的距离L计算面波的到时窗口(t_min=L/v_max,t_max=L/v_min),时窗长度尽量选择为最大到时或略大于最大到时,这样无论是最大速度的波还是最小速度的波,都可以最大程度地保证在该时段内传播三个波长的距离。对于本次野外试验,根据当地的地质资料,地层横波速度范围是200~3 000 m/s,检波器阵列长度约为600 m,根据上述原则选定的时窗长度约为3 s,可恢复最大波长为200 m,频率范围为1~15 Hz的面波,与图6a上显示的主要频率范围相同,满足计算要求。为了验证窗口长度选取的合理性,分别将时窗长度设为1、3、6、12、30、60 s,依次进行互相关叠加运算,并通过对互相关结果进行滤波、求取信噪比来衡量互相关计算结果的信号质量。滤波时采用的高斯滤波器如下:

式(5)中,α是影响滤波器带宽的因子,如果选择过小,导致滤波带宽过大,滤波后会包含很多其他频率的信号,如果α选择过大,则导致频率成分少。

对野外实测中第六条测线上的数据进行处理,采用不同时窗长度依次计算第一道信号与后续13道信号的互相关函数。结果如图9所示,图9a中面波发育不全,无法清晰的展现出来,同时由于互相关叠加次数过多,严重影响了数据处理的效率。图9c和图9d中时窗长度逐渐变长,迭代次数减少,信噪比也逐渐降低。不同时窗长度的信噪比结果如表3所示,信噪比随着时窗长度的增加,出现先增后减的现象,当时窗长度为3 s时信噪比最高。图10为经过高斯滤波后的互相关信号,可以看出经过滤波处理后的互相关记录更为平滑,低频范围内的频谱信息更为突出。

对于理论模型模拟信号(即2.1节得到的模拟数据,如图2所示),得出互相关结果后,便可根据式(2)提取频散曲线。提取结果如图11所示。图11中从上至下m的取值依次为-2、-1、0、1,通过Knopoff ^[23]提出的Schwab-Knopoff快速计算方法得出模拟记录的理论频散曲线,如图11中的黑色实线所示。通过对比可以得知,当m取-1时,频散曲线与理论曲线拟合度最高。

图12中的频散曲线分别为:只经历时间域归一化处理的数据的频散曲线(蓝色实心圆点)、只经历频率域归一化处理的数据的频散曲线(绿色实心圆点)以及经历完整处理流程的频散曲线(红色实心三角)。可以看出只经历部分处理流程的频散曲线与理论频散曲线的吻合度较差。表4给出了图12中三条频散曲线的评价结果,通过均方误差与相关系数可以看出,经历过完整的处理流程后得到的频散曲线精度最高。

在处理野外实测背景噪声信号时,信噪比高的信号,零点数量变化的情况不严重,m通常只需要取到±1,而信噪比低的信号,就需要设置更大绝对值的m。用测线6上的1号和6号检波器数据进行互相关提取频散曲线,提取结果如图13所示。

根据当地的地质资料判断横波速度大致范围,并将不同m值的被动源频散曲线与主动源方法提取的频散曲线进行比较,挑选出具有良好一致性的被动源频散曲线。对比结果如图14所示,蓝色线条为运用两道法得到的被动源频散曲线,青色线条为通过主动源方法提取得到的基阶频散曲线,灰色线条为主动源方法提取得到的高阶频散曲线。可以看出当m=-2时,被动源频散曲线与主动源基阶频散曲线可以很好的拟合。之后依次对每一道互相关信号的频谱曲线进行筛选,通过改变m的取值,将所有频散曲线拟合在一起,从而得到一条测线上所有的频散曲线。

利用理论模型被动源数值模拟数据和野外实际观测数据,采用两道法计算格林函数,运用Aki公式的原理提取瑞利面波的频散曲线,所得结果与主动源方法计算得到的频散曲线具有很好的一致性。总结归纳出一套被动源瑞利波频散曲线质量控制流程,通过理论模型和实测数据试算,得到如下结论:

1)文中给出的背景噪声数据处理流程有效地改善了被动源瑞利波频散曲线的提取质量,验证了方法的可行性和有效性。

2)对数据进行时间域归一化时,采用“one-bit”法会取得更好的效果,而对数据进行频率域归一化时,滑动绝对平均法表现更好。

3)互相关运算时,时窗长度的选取原则应在确保面波信息完整的情况下尽可能多地增加叠加次数。实际应用中,可根据数据采集阵列长度与地层的最小横波速度的比值计算得到。

4)改变计算瑞利波相速度公式中的m,结合当地地质资料判断横波速度大致范围,挑选与主动源基阶频散曲线具有一致性的被动源频散曲线,可以消除因频谱实部零点的缺失或增加造成的干扰。