芦山余震震源参数及震源区品质因子反演

温瑞智，王宏伟，任叶飞，冀 昆

(中国地震局工程力学研究所，150080哈尔滨)

2013年4月20日四川省芦山县发生MS7.0级地震，芦山主震震中位于龙门山断裂带西南端的彭县-灌县断裂上，主震过后余震活动十分频繁，余震沿发震断层向主震两侧延伸，主要分布在长约32 km、宽约15～20 km、深度为5～24 km的范围内［1］，余震破裂类型以逆冲型为主［2］.芦山地震发震区域位于龙门山断裂及鲜水河-小江断裂的“Y”形交界处，在过去的40年里，距离芦山地震震中200 km范围内发生过多次MS6.0级以上地震，其中包括2008年MS8.0级汶川地震及其余震，这一区域地震多发人口密集，地震危险性很高.

我国测震资料相对丰富，国内研究人员大多基于Atkinson方法与遗传算法采用测震记录反演我国地震多发地区的品质因子及震源参数，比如甘肃、云南和四川地区［3-5］.随着我国高质量强震动记录的增加，采用强震动记录反演中小地震震源参数及品质因子已经逐渐得到推广［6-8］.我国数字强震动台网在龙门山地区密度相对较大，芦山余震序列中台站多次触发，记录到2.0～5.4级余震176次，共收集到超过1 000组加速度记录，这为反演芦山余震震源参数及品质因子提供了基础数据.本文基于参考事件的广义反演方法估计了发震区域剪切波品质因子及震源参数，包括地震矩、拐角频率、应力降、破裂半径及地震波能量，进一步确定了震源参数的定标关系.

1 广义反演方法

自由表面水平地震动剪切波加速度傅里叶幅值谱可表示为

式中:Oij(f)表示第j个台站观测到第i个地震的记录剪切波水平方向加速度傅里叶幅值谱;Si(f)表示第i个地震的加速度震源谱;GS(Rij)表示几何扩散;Rij表示第i个地震到第j个台站的震源距;exp(-πfRij/Q(f)Vs)表示地震波非弹性衰减;Q(f)表示与频率f相关的品质因子;Vs表示震源处剪切波速，取为3.6 km/s;Gj(f)表示第j个台站的场地反应.

本文几何扩散采用Atkinson等［9］给出的三段线性衰减曲线，表示为

式中:R01=1.5D，R02=2.5D，b1=1.0，b2=0.5，D表示地壳厚度.芦山地震序列震中附近区域地壳结构复杂，地壳厚度从西北侧的52.5 km减小至东南侧的 41.5 km［1］，近似以平均地壳厚度47 km作为地壳厚度，则R01=70.5 km，R02=117.5 km.

式(1)两边同时取对数，得到线性叠加形式:

式(3)矩阵形式为

式中:A是每一行包含3个非零项(两个1和-πfRij/Vs)的稀疏矩阵;X表示式(3)右边所有未知量的向量;b表示式(3)左边所有已知量的向量.该矩阵方程可详细表示为式(5)所示形式.采用奇异值分解方法求解式(4)，在每个频率上确定Ⅰ(地震个数)+J(台站个数)+1(Q值)个未知数，由于存在一个未加约束的自由度，需要考虑震源与场地间的权衡问题，通常采用参考场地法或参考事件法来解决.芦山余震序列中6个触发基岩台站水平/垂直谱比法估计的场地反应均不能满足参考场地应在整个频段上不存在场地放大的要求，因此本文采用参考事件法.

2 数据及数据处理

采用以下原则选取合适的强震记录用于广义反演:

1)记录震源距范围为25～100 km.芦山余震震源深度主要分布于5～24 km，为避免震源深度不确定性对震源距计算的影响和减少截取的剪切波中面波引起的干扰分别设定最小和最大震源距限制［10］.

2)对原始数据进行零线校正及0.1～30 Hz巴特沃斯带通滤波处理后，为减小噪声对记录的影响和避免场地出现非线性反应选取三分量峰值地面加速度均满足5～100 cm/s2的记录［11］.

3)为减小反演结果的离散性，在满足1)和2)的基础上，选择同时满足不少于4次地震中触发的台站以及不少于4个台站记录到的地震所对应的记录.

根据上述3个原则选取了25个台站在34次地震中记录到的262组记录，台站记录震源距分布见图1，图中阴影区域表示5 km间隔范围内记录数量分布直方图围成的区域，64表示最大频数.由图1可见，纵坐标51DXY及以上台站(A组)，记录震源距大多不超过70 km，主要分布于25～45 km范围内;纵坐标 51DXY以下台站(B组)，震源距主要分布于80～90 km范围内.A、B两组台站记录震源距差异决定了传播路径不同，由于品质因子对地壳介质横向和垂直方向的不均匀性特别敏感，A、B两组台站记录品质因子差异较大［12］.

图1 台站记录震源距分布

本文在A组台站数据中选取了30次地震中15个台站记录到的179组记录用于广义反演.台站及余震震中分布见图2.所有传播路径分布均匀且交织在一起，方位角对震源影响及传播介质不均匀性对品质因子影响可忽略.A组记录震源距不超过R01(即70.5 km)，几何扩散GS(R)取为R-1.

图2 台站、地震震中位置及震源与台站间传播路径分布

广义反演方法需截取记录的剪切波部分，剪切波抵达时间和结束时间分别定义为Husid函数中地震波能量开始急剧增加的点和累加均方根函数开始下降的点对应的时间［13］.在截取的剪切波前后各乘以10%的剪切波持时的余弦边瓣窗口以消除截断误差，计算剪切波两个水平方向的傅里叶谱，采用b=20的Konno和 Ohmachi窗口函数平滑傅氏谱［14］，并以两个水平方向谱的矢量和作为记录水平傅氏谱.

3 参考事件震源谱

Ren等［11］以62WIX台站作为参考场地基于广义反演方法估计了96个汶川余震加速度震源谱.本文选取2013年4月20日09:37:29发生的4.9级地震作为参考事件，从96个汶川余震中选取与参考事件震级相差不超过 0.1(4.9±0.1)的余震13个.采用下文4.1节提到的网格搜索方法确定选取的13个余震的平均位移震源谱的理论震源参数，地震矩M0=2.427×1023dyne·cm、拐角频率fc=0.89 Hz、高频滚降系数 γ=2.0，以该理论震源谱作为参考事件的实际震源谱用于广义反演.

4 结果分析及讨论

4.1 震源参数

基于参考事件的广义反演方法分离了上文所选15个台站的场地反应、30个地震的加速度震源谱以及区域品质因子.反演得到的余震位移震源谱见图3，图中可见，位移震源谱较符合Aki［15］提出的ω2震源谱模型.采用震源参数(M0、fc)及γ表示的实际位移震源谱的理论形式:

图3 反演得到的30次余震位移震源谱

式中:RΘΦ为点源辐射图型因子，随方位角呈对称性变化，在此取平均值0.55;V表示地震波水平分量含剪切波能量的比例，取为0.707;F为半空间表面放大，本文取为1;ρs为震源处介质密度，取为2 700 kg/m3;R0为参考距离，通常取为1 km.

其中矩震级MW与地震矩有如下关系［16］:

采用网格搜索方法在0.5～20 Hz频段内搜索满足理论位移震源谱与实际位移震源谱相对面积差最小的MW、fc及γ，即

式中:n表示0.5～20 Hz内所有频率点数;fi表示第i个点的频率;SS(fi)与SO(fi)分别表示搜索的理论位移震源谱及基于广义反演方法确定的实际位移震源谱.理论震源谱的搜索过程中MW变化范围为余震震级±0.5，变化步长0.01;fc变化范围为0.01～5.0 Hz，变化步长为 0.01 Hz;γ 变化范围为2±0.3，变化步长 0.1.

最小二乘拟合地震矩与拐角频率可得:

Aki［15］认为M0fc3是一个与应力降有关的常数，故将式(9)斜率固定为-3.00后拟合得:

拟合结果见图4，M0fc3=1.73×1023dyne·cm·s-3，略低于Dutta等［17］研究阿拉斯加中南部地区中小地震得到的M0fc3=2.09×1023dyne·cm·s-3和Hassani等［13］研究伊朗中东部地区中小地震得到的M0fc3=2.48×1023dyne·cm·s-3.芦山余震平均应力降为 3.14 MPa，略低于喻畑等［8］研究13次汶川余震 (MW＞5.0)震源参数得到的3.80 MPa平均等效应力降.

图4 地震矩与拐角频率关系(阴影区域表示加减一倍标准差范围)

根据Brune［18］提出的圆盘形应力脉冲震源模式，计算了余震震源半径r及应力降Δσ.地震矩与震源半径r的关系见图5，图中同时给出了余震的应力降分布范围(1.0～10.0 MPa).应力降与震级没有明显的相关关系，与Moya等［19］研究阪神地震余震得到的结论一致.

图5 地震矩与震源半径散点分布图

基于Izutani-kanamori理论谱方法同时考虑地震波能量补偿将积分上限提高至最大拐角频率(拐角频率上限为 5 Hz)的 10 倍，即 50 Hz［20］.本文分别计算了地震波能量Es和视应力σA.

地震波能量与地震矩的关系见图6，双对数坐标下线性拟合得到 logEs=0.82logM0-7.7，相关系数R=0.91，拟合直线斜率0.82 接近1.0，近似有Es∝M0，固定斜率为1拟合得

Es/M0平均值为2.19×10-12，与 Hassani等［13］得到的Es/M0=2.50×10-12较为接近，但明显大于Dutta 等［17］得到的Es/M0=1.20×10-12.

图6 地震波能量与地震矩关系

本文确定的视应力为 0.17～2.21 MPa，低于程万正等［21］研究2000—2004年四川地区中小地震得到的视应力变化范围 0.1～10.0 MPa，推断可能是发生于主震破裂面的余震初始应力降低所致.根据应力降与视应力的计算结果，线性拟合得Δσ=6.78σA-1.05，R=0.91，表明应力降越高视应力越大，单位地震矩辐射的地震波能量越大，本文余震σA/Δσ的平均值为0.258.

不考虑面波震级MS和地方震级ML之间的换算，通过最小二乘法拟合可得到地震面波震级与矩震级的关系:MW=0.876 2MS+0.563 9(R=0.83)，如图7所示，本文确定的矩震级与面波震级关系与Ren等［11］反演汶川余震得到的结果较为一致，整体上矩震级拟合结果与面波震级一致，对于小震级余震，矩震级略高于面波震级;对于中等震级余震，矩震级略低于面波震级.

图7 矩震级与面波震级关系

4.2 品质因子

品质因子Q(f)通常表示为Q0fn的形式，本文反演确定的芦山余震区域剪切波品质因子见图8，拟合得到0.5～20 Hz频段品质因子为Q(f)=31.867f1.0375，芦山地震近场区地震波衰减表现出吸收快且与频率依赖性强的特点(低Q0高n).芦山发震区与汶川地区临近，其频率相关的品质因子对频率的依赖性相近，即n值较一致，一定程度上可推断本文反演的品质因子较为可靠.与四川其他地区相比［4，5，11］，本文反演确定的芦山余震区域品质因子结果较小，主要原因有两方面:1)广义反演方法确定的品质因子反映了所有传播路径的平均非弹性衰减，本文选取的余震震源深度主要分布于10～25 km范围内，该研究区域地壳平均厚度约为47 km，震源距集中在25～45 km范围，地震波传播路径主要集中在上地壳，而其他研究成果采用较远震源距的记录，地震波传播路径主要集中在下地壳，众所周知地壳深处Q值相对较大;2)本文研究区域集中于芦山地震发震区，该区域地壳活动强烈，一般认为，在地壳活动相对稳定地区Q值相对较高，在地壳活动强烈的地区Q值相对较小.

图8 频率相关的品质因子

5 结论

1)采用芦山强余震记录基于广义反演方法估了余震震源参数 (M0、fc、r、Δσ、Es、σA)及该区域剪切波品质因子Q(f)，并给出了震源参数间的定标关系.

2)芦山余震M0fc3=1.73×1023dyne·cm·s-3，相应于3.14 MPa的平均应力降，略低于汶川地震余震平均等效应力降;应力降主要在0.1～10.0 MPa范围内变化，与震级没有明显的相关性;视应力与应力降正相关，由芦山余震视应力低于2000—2004年四川余震结果推断发生于主震断层面上的余震初始应力可能较低.芦山余震震源参数M0fc

3及Es/M0与 Hassani等［13］研究伊朗中东部中小震级地震的结果十分接近，初步推断两地区余震震源破裂过程以及区域构造环境可能存在相似性.

3)0.5～20 Hz频段内本文估计的剪切波品质因子Q(f)=31.867f1.0375，分析表明地震波衰减吸收较快且与频率的依赖性较强;本文获得的芦山地震区域品质因子低于四川其他地区的结果，主要是由于所用记录震源距相对较小且该区域地壳活动较为活跃.另外，需要说明的是，本文研究区域位于青藏高原与四川盆地过渡地带，存在地壳介质横向不均匀性.广义反演确定的品质因子仅体现该区域的平均水平，该方法无法考虑单一区域内品质因子的横向变化性.考虑到本文研究区域相对较小(半径约70 km)，暂不考虑地壳介质横向不均匀性对本文反演结果的影响.

［1］陈晨，胥颐.芦山MS7.0级地震余震序列重新定位及构造意义［J］.地球物理学报，2013，56(12):4028-4036.

［2］林向东，葛洪魁，徐平，等.近场全波形反演:芦山7.0级地震及余震矩张量解［J］.地球物理学报，2013，56(12):4037-4047.

［3］华卫，陈章立，郑斯华，等.三峡水库地区震源参数特征研究［J］.地震地质，2010，32(4):533-542.

［4］乔慧珍，张永久，程万正.川西北地区介质衰减特性研究［J］.地震地磁观测与研究，2006，27(4):1-7.

［5］张永久，乔慧珍，程万正.四川盆地地区介质衰减特性研究［J］.地震研究，2007，30(1):43-48.

［6］章文波，谢礼立，郭明珠.利用强震记录分析场地的地震反应［J］.地震学报，2001，23(6):604-614.

［7］刘杰，郑斯华，黄玉龙.利用遗传算法反演非弹性衰减系数、震源参数和场地响应［J］.地震学报，2003，25(2):211-218.

［8］喻畑，李小军.汶川地震余震震源参数及地震动衰减与场地影响参数反演分析［J］.地震学报，2012，34(5):621-632.

［9］ATKINSON G M，MEREU R F.The shape of ground motion attenuation curves in southeastern Canada［J］.Bulletin of the Seismological Society of America，1992，82(5):2014-2031.

［10］MCNAMARA D，MEREMONTE M，MAHARREY J Z，et al.Frequency-dependent seismic attenuation within the Hispaniola island region of the Caribbean sea ［J］.Bulletin of the Seismological Society of America，2012，102(2):773-782.

［11］REN Yefei，WEN Ruizhi，HIROAKI Yamanaka，et al.Site effects by generalized inversion technique using strong motion recordings ofthe 2008 Wenchuan earthuake［J］.Earthquake Engineering and Engineering Vibration，2013，12(2):165-184.

［12］PETUKHIN A，IRIKURA K.A method for the separation of source and site effects and the apparent Q structure from strong motion data ［J］.Geophysical Research Letters，2000，27(20):3429-3432.

［13］HASSANI B，ZAFARANI H，FARJOODI J，et al.Estimation of site amplification，attenuation and source spectra of S-waves in East-Central Iran ［J］.Soil Dynamic and Earthquake Engineering，2011，31(2011):1397-1413.

［14］KONNO K，OHMACHI T.Ground-motion characteristics estimated from spectral ratio between horizontal and vertical components ofmicrotremor ［J］.Bulletin of the Seismological Society of America，1998，88(1):228-241.

［15］AKI K.Scaling law of seismic spectrum［J］.Journal of Geophysical Research，1967，72(4):1217-1231.

［16］HANKS T C，KANAMORI H.A moment magnitude scale［J］.Journal of Geophysical Research，1979，84(B5):2348-2350.

［17］DUTTA U，BISWAS N，MARTIROSYAN A，et al.Estimation of earthquake source parameters and site response in Anchorage，Alaska from strong motion network data using generalized inversion method ［J］.Physics of the Earth and Planetary Interiors，2003，137:13-29.

［18］BRUNE J N.Tectonic stress and the spectra of seismic shear-waves from earthquake ［J］. Journal of Geophysical Research，1970，75(26):4997-5009.

［19］MOYA A，IRIKURA K.Estimation of site effects and Q factor using a reference event［J］.Bulletin of the Seismological Society of America，2003，93(4):1730-1745.

［20］华卫，陈章立，郑斯华.2010年4月14日青海玉树7.1级地震序列中小地震辐射能量的估计［J］.地球物理学进展，2012，27(1):8-17.

［21］程万正，陈学忠，乔慧珍.四川地震辐射能量和视应力的研究［J］.地球物理学进展，2006，21(3):692-699.