地震破裂过程就是能量释放的过程。随着地震震源研究的不断深入,地震学家们不再满足于只对震源时间函数的了解,他们还试图通过对震源时空破裂过程的反演来了解地震破裂期间发生在断层面上的破裂传播行为以及断层面上不同位置的滑动位错大小分布。中外学者们已研究出了多种由地震观测资料反演震源时空过程的方法。这些方法可分为直接反演和间接反演两类,即由观测地震图直接反演震源时空过程(Kikuchi,Kanamori,1982,1991; Hartzell et al,1991,1996)以及由观测地震图提取震源时间函数,再由震源时间函数反演震源时空过程(Mori,Hartzell,1990; Mori,Frankel,1990,许力生,1995,1996; 周云好,2002a,2002b)。Kikuchi和 Kanamori(1982,1991)假设大地震是由在时间、空间上分布的一系列小破裂事件组成,从地震图上将这些小破裂事件识别出来,分别定出每个小破裂事件的时空位置和震源机制,并用迭代方法求出这些小破裂事件的标量地震矩和震源时间函数,从而得出大地震时空发展的总体过程,该方法不需要事先假定震源时间函数的形状,但识别每个小破裂事件依赖于研究者的经验。Hartzell 等(1991,1996)假设大地震过程中震源机制是不变的,将大地震的平面断层分解为许多子断层,由多个地震台的观测地震图直接反演各子断层的滑动量,从而给出大地震断层面上的滑动分布,但这种方法需要设定破裂起始点位置、破裂传播速度、各子断层滑动函数的形状和持续时间。Mori和Hartzell(1990)也假设大地震在破裂过程中震源机制是不变的,将大地震的平面断层分解为很多子断层,与Hartzell 等(1991,1996)的方法不同的是他们用观测地震图与格林函数提取各台依赖于台站方位的大震震源时间函数,再由时间函数反演各子断层的滑动量,从而求出大地震断层面上的滑动分布,该方法同样也需要设定破裂起始点位置、破裂传播速度、各子断层滑动函数的形状和持续时间。

许力生和陈远泰(1996)、周云好等(2002a,2002b)对上述方法进行了改进,在用震源时间函数反演大地震平面断层上的滑动分布时,不需要预先设定破裂的传播速度、各子断层滑动函数的形状和持续时间等,从而减少了人为因素。

2000年1月15日7时37分5秒,云南省姚安县发生M_S6.5地震。笔者利用全球11个台站的长周期记录的P波波形资料,通过直接波形反演方法(张勇等,2008)对此次地震的震源破裂过程进行了反演。

地震时空破裂过程反演的基本思想是将断层面均匀地划分为若干个子断层,然后利用各个台站接收到的波形信息和相应的格林函数,确定各子断层的滑动幅度和时间函数。综合所有子断层的破裂信息,就得到地震的整个破裂在时间和空间上的分布。位移表示为(陈运泰,顾浩鼎,2000)

U(t)=G(t)·s(t).(1)

式中,U(t)为观测记录,G(t)为特定震源机制对应的格林函数,s(t)为各子断层的矩率时间函数。反演过程中,为了使反演保持稳定,我们在空间和时间域加上光滑约束(Yugi et al,2004)以及地震矩最小约束:

[WU

0

0

0]=[WG

λ₁D

λ₂T

λ₃Z][s].(2)

其中,W为对各台站数据施加的权重,D为空间光滑矩阵,T为时间光滑矩阵,Z为约束地震矩最小的矩阵。λ₁和λ₂分别为空间和时间上的光滑权重,λ₃为地震矩最小约束权重。实际反演过程中,通过调整W来确保矩阵WG绝对值的平均数值在1左右,使得约束权重λ₁、λ₂和λ₃的数值大小与1相当。由于整个系数矩阵较为庞大,我们选择一种反演计算效率较高的带约束的梯度法(Ward,Barrientos,1986)来求解(2)式所表示的线性问题。在梯度法迭代的每一步,我们都对解施以非负的约束条件,即在正值空间中搜索最优解s。

笔者从IRIS(Incorporated Research Institution for Seismology)网站上获取了50个震中距离小于90°的远场波形资料,从这些数据中选用了11个台站的震相清晰的垂直向的P波波形资料。选用台站和震中位置如图1所示,台站参数见表1。

图1 选用台站及地震分布图
Fig.1 The chosen stations and earthquakes distribution

表1 11个选用台站的参数
Tab.1 The parameters of 11 chosen stations

王绍晋等(2002)的研究结果表明:姚安M_S6.5地震序列的主要发震断层是在走向为北50°西、倾角陡立、在南南东向接近水平的主压应力作用下具有以右旋走滑为主的错动性质的构造断裂。毛燕等(2006)从本序列余震的震源机制解得出NW—NWW方位的节面占主导地位,再结合余震分布尺度可以判断发震断裂是北西向的断裂,且应力场以南南东向接近水平的主压应力为主。我们根据不同资料来源给出的地震的震源机制解(表2)、震中分布(图2)、震源机制的分析结果以及地震裂缝等,选取哈佛大学矩张量结果的断层A为发震构造。反演中地震参数为:震中在25.61°N、101.06°E,震源深度30 km,矩震级5.9,走向118°,倾角84°,滑动角-168°。

图2 3.0≤ML≤5.0余震震中分布图F1:建水—楚雄断裂; F2:通海—牟定断裂; F3:滥泥旦断裂; F4:马尾箐断裂; F5:花椒园断裂
Fig.2 Epicenter distribution of aftershocks

表2 2000年姚安地震震源机制解
Tab.2 The focal-mechanism solution of Yaoan earthquake in 2000

挑选出的11个台的记录信噪比较高,反演前我们对该记录进行了预处理,去除了直流分量和倾斜分量,并去除了仪器响应,按5 sps重新采样。

为了客观地反演出实际破裂过程,反演时我们将破裂区域范围选得足够大,使之大于实际破裂区域。整个区域设定为沿断层走向取65 km,沿断层倾向取50 km,再将整个区域划分为130个子断层,即走向上进行13等分,倾向上进行10等分,每个子断层的尺度大小为5 km*5 km。震中位于走向第7、倾向第6的子断层上。

首先运用反射—折射率方法计算全球的格林函数(Kennett,1983)。地震的震源时间函数为有限断层面上所有点的震源时间函数的叠加,地震矩随震源破裂时间的释放率函数如图3所示,可以看出地震矩主要在破裂开始0.2 s后开始释放,总的破裂持续时间约为5.6 s,但主要分两次进行,第二次破裂强度大于第一次。

图3 地震矩释放率函数
Fig.3 Earthquake moment rate function

图4 静态滑动位移分布图
Fig.4 Static slip distribution

图6显示了断层滑动速率的时空分布情况。每幅子图的时间间隔为0.2 s,可以看出最大滑动速率达0.137 m/s。首先是震中区域的滑动,然后是震中以东区域先呈现近南北向的滑动,约0.4 s后立即转为近北西向的滑动,但滑动的总体延伸趋势为北西走向。

图5 应力降分布图
Fig.5 Stress drop distribution

图6 滑动速率时空分布图
Fig.6 Temporal and spatial distribution of the slip rate

图7 不同台站观测波形与合成波形的对比(a)HYB;(b)TLY;(c)AKSU;(d)KMNB;(e)HLQI;(f)SSE; (g)TPUB;(h)ULN;(i)WMQ;(j)MAKZ;(k)HARA
Fig.7 The comparison between observed waves and synthesis waves

综上所述,2000年姚安M_S6.5地震的震源破裂过程反演结果显示:地震的破裂时间约为5.6 s,但地震矩释放主要在两个时间段内进行,且第二次的破裂强度大于第一次; 静态滑动位移和应力降均显示存在震中区域及震中以东区两个变化比较明显的区域,这与地震过程中出现两次破裂较为吻合,在震中以东的区域出现静态滑动位移和应力降最大值,也即这个区域是地震破裂的主要区域,在该区域里东南端的值比西北端的大,由此可以看出破裂主要是由北西走向的断裂引发。

这一研究结果与秦嘉政(2001)得出的姚安地震的发震构造是北西向区域的次级构造的结论一致。

余震的震源机制解(毛燕等,2006,2007)显示有少数余震的发震应力场是以北东向的主压应力为主。静态滑动位移、应力降及滑动速率时空分布图均显示震中区域先滑动后,震中以东的区域产生了滑动,这可能与北东向断裂的活动有关,但整个破裂过程还是以北西向的断裂活动为主(毛玉平,万登堡,2001)。

致谢:向给予指导和帮助的中国地震局地球物理研究所许力生研究员、张勇副研究员及青海省地震局张晓青副研究员表示感谢!