地震波衰减性质的研究是地震学研究的重要课题。地震仪器记录到的地震波,包含了地震震源效应、地震波的传播路径效应、台站场地响应及仪器响应。在使用地震波资料研究震源性质时,必须要扣除地震波传播路径效应、台站场地响应及仪器响应的影响。地震波传播的路径效应(地震波衰减),除了随距离存在几何衰减外,还有一个重要的影响因素即介质的非弹性衰减,可用介质品质因子Q值来度量。Q值是地球介质的基本物理参数之一,是通过远离震源的观测资料对地震记录进行定量分析和研究震源性质所必需的重要参数。同时,地震波的衰减特征对于地震危险性分析也具有十分重要的意义,是地震工程研究人员最关心的基础资料之一。由于剪切波的振幅通常比P波振幅大,在有的情况下可达5倍左右,因此对于地震工程来说,研究剪切波衰减的特征将具有更重要的意义。

理论研究和实际观测结果表明,场地响应是影响地震破坏的重要因素,对地震破坏起着重要的作用。已有不少震害调查结果证实,大地震时松软地基上的建筑物的破坏率要比坚硬地基上的建筑物的破坏率高得多。Rogers等(1979)采用内华达试验场核爆破记录的波谱,计算了Long Beach的场地响应,结果表明在0.2～6 s波段冲积层相对于岩石的场地放大倍数可高达11倍。

随着我国区域数字地震台网的建立与完善,地震工作者获取了大量的地震波形资料,数字地震波资料的分析与应用研究逐渐成为一项重要工作,而地震波非弹性衰减Q值和场地响应特征研究是深入开展该项工作的一个重要基础。笔者根据甘肃数字地震台网记录到的青藏高原东北缘地区的中小地震的波形资料,采用互相衔接的三段几何衰减模型,通过在频率域内的分析,采用Atkinson方法和遗传算法(Atkinson等,1992; 黄玉龙等,2003; 刘杰等,2003),研究了青藏高原东北缘地区剪切波的非弹性衰减Q值以及场地响应的征。

本研究的资料取自甘肃数字地震台网的波形记录,该台网由21个子台组成。这些台站都位于基岩上,表1列出了各台站所用仪器的型号和台基岩性。每个台站都安装了三分量数字地震仪,大部分地震仪在1～20 Hz之间具有速度平坦的响应。地震信号以50 Hz的采样率通过有线或无线通讯实时地传到设在兰州市内的台网中心。

表1 地震台站、仪器型号和台基岩性一览表

甘肃数字地震台网从2001年开始有正式的地震记录。通过地震目录和波形考察,我们初步选出3.0级以上的地震159个共749条记录作为候选分析资料。然后,通过信噪比的计算,挑选出81个地震的437条记录用于分析(图1)。这些地震的震级范围为M_L3.0～5.1,震中距为10～500 km。

图1 青藏高原东北缘地区21个台站和81个地震的分布及437条地震记录的传播路径

对于S波的两个水平分量(NS、EW),首先进行带通滤波(Butter滤波器,带宽0.1～22 Hz)和水平校正处理,然后取“S窗”和“噪声窗”,进行傅立叶变换信噪比分析。

如图2所示,把从S波开始到包括S波总能量的90%的时间段定义为“S窗”。对于不同地震,由于震源深度和震中距不同,“S窗”内包含的震相也不尽相同。对于同一次地震,由于各台站的震中距不同,“S窗”的持续时间也不同。为了得到具有相同频率间隔的振幅谱,笔者采用了平移窗谱方法(苏有锦等,2006)。具体步骤如下:把“S窗”内的波形信号分成若干个包含有256个采样点的小段,并使相邻信号段有50%的重叠; 对于采样率为50 Hz的地震记录来说,可得到每个信号段的时间长度是5.12 s; 在每一信号段的起始和末尾各加5%的cos边瓣后,通过傅立叶变换得到每个信号段的傅立叶谱,这样对于每一个台站的记录就可以得到相同频率间隔(0.196 Hz)的傅立叶谱; 最后,对每个信号段的傅立叶谱进行仪器校正,并通过下式得到整个“S窗”内信号的速度振幅谱:

V^-(f)={(∑ⁿ_i=1v_i²(f))·T/(nt)}^1/2.(1)

式中v_i(f)是经过仪器校正的第i个信号段的傅立叶谱,T为“S窗”的持续时间。该“S窗”内包含了n个时间长度为t、每个时间长度包含256个采样点的信号段。取P波初动前256个采样点(1个信号段)的噪声信号,通过(2)式得到与信号相同频率间隔的噪声谱:

N(f)={n²(f)(t_s)/t}^1/2.(2)

式中,n(f)是经过仪器校正的256个数据点的噪声傅立叶谱。然后,得到经过噪声校正的速度振幅谱:

V(f)={V^-²(f)-N²(f)}^1/2.(3)

由于甘肃数字地震台网的资料是速度记录,最后还需要把V(f)除以2πf,把速度谱转换成位移谱。对于S波的两个水平分量分别进行上述处理,并通过(4)式得到S波水平分量合成位移谱:

A(f)={u²_EW(f)+u²_NS(f)}^1/2.(4)

用以上处理方法,按每次地震至少被3个台站记录到、每个台站至少有4条记录的原则,挑选波形较好、能经过信噪比检验的地震记录用于本研究。图2给出了其中的一个实例,为文县台记录到的2005年8月12日卓尼3.6级地震的两个水平分量原始波形及其S波合成位移谱。图2a中,虚线间为“S窗”的范围; 图2b中,上部曲线为经过噪声校正的位移谱曲线,下部曲线为噪声谱曲线。

图2 文县台记录到的2005年8月12日卓尼 3.6级地震的两个水平分量原始波形(a) 及其合成位移谱(b)

根据地震波谱的一般表达,经过上述资料预处理后,任一地震观测记录谱与震源谱均有如下关系:

lgO_ij(f)=lgS_i(f)-lgG(R_ij)-C(f)R_ij

-lgG_j(f).(5)

其中f是频率,O_ij(f)是第j个台站观测到的第i个地震的谱振幅,S_i(f)是第i个地震的震源谱振幅,G(R_ij)为几何衰减函数,C(f)为非弹性衰减系数,R_ij是第i个地震至第j个台站的震中距,G_j(f)是第j个台站的场地响应。

非弹性衰减系数C(f)与区域介质的品质因子Q(f)之间的关系为

Q(f)=(lg(e)πf)/(C(f)β).(6)

式中,β为S波速度,取3.5 km/s。

参考Atkinson等(1992)、黄玉龙等(2002)文献,几何衰减函数G(R_ij)采用互相衔接的三段几何衰减函数表示:

G(R)={R^-b₁, R≤R₀₁;

R-b101·R-b201·R-b2, R01<R≤R02;

R-b101·R-b201·R-b202·R-b302·R-b3, R>R02.

上式中R≤R₀₁时,对应于直达波的几何衰减; 当R₀₁<R≤R₀₂时,对应于过渡区,在该震源距范围内,直达波加入了在地壳几个间断面和莫霍面上的反射波; 当R>R₀₂时,对应于多次折射反射波。根据青藏高原东北缘地区平均地壳厚度H约50 km(李永华等,2006),有关模型参数取值为b₁=1.0、b₂=0.0、b₃=0.5、R₀₁=1.5H=75 km、R₀₂=2.5H=125 km。如果把该几何衰减函数G(R_ij)代入式(5),则式(5)变为对参数C(f)(非弹性衰减系数)和G_j(f)(台站场地响应)进行联合反演的问题。

定义残差为

k_ij=[lgS_i(f)]_j-lgS_i(f)^-.(7)

式中lgS_i(f)^-是第i个地震的震源谱振幅对数的平均值,是对记录到该地震的所有台站计算得到的lgS_i(f)求平均。

非弹性衰减系数C(f)的求解采用计算(8)式求极小值得到,即

sum=∑_i∑_j〖JB<2|〗k_ij〖JB>2|〗.(8)

采用Atkinson方法反演品质因子Q(f)值和场地响应的计算步骤可归纳为:

(1)先将所有台站场地响应G_j(f)设为1,选择合适的参数C(f),使(8)式的残差总和极小。

(2)利用得到的参数,可求得场地响应:

lgGj(f)=(∑mjj=1kij)/(mj).(9)

其中m_j为第j个台站记录的地震事件数。

(3)将场地响应的计算结果代入,重新计算C(f)。重复(2)、(3)步,使残差总和达到最小。

通过上述迭代反演可同时求得该地区的非弹性衰减系数C(f)和各台站的场地响应,进而利用(6)式得到区域介质的品质因子Q。在考虑频率依赖的Q模型时,通常用频率的幂函数,即Q(f)=Q₀f ^η的形式来拟合Q与频率的关系。图3给出了其计算结果,由图可见结果拟合较好。

图3 品质因子Q(f)与频率f的关系

Q(f)=564.7f ^0.3.(11)

图4是迭代反演得到的21个子台的场地响应。由图可见,除了武都、合作、湟源、定西和迭部等台的场地响应总体上接近或小于1以外,其余台站的场地响应均存在较明显的放大效应。

本研究得到的青藏高原东北缘地区Q值与频率f的关系,总体上与广东地区的结果(Q(f)=150f ^0.5)和欧洲中部地区的结果(Q(f)=400f ^0.42)较为接近。

采用Atkinson 方法计算非弹性衰减系数C(f)或品质因子Q时,由于各台站的场地响应是通过与平均值比较来确定的,对于场地响应较小的台站(如基岩台站),该方法所得到的场地响应的对数将是负值,即台站响应将小于1(图4),因此得到的场地响应是相对的,不是真实的场地响应。这表明研究得到的21个台站的场地响应,除了武都、合作、湟源、定西、迭部等台的场地响应总体上没有放大效应以外,其余台站的场地响应均存在较明显的放大效应。

图4 由437条记录通过迭代反演得到的21个子台的场地响应