2008年5月12日在四川省汶川县发生的8.0级地震是半个世纪以来中国境内最大的一次地震,造成了巨大的经济损失和人员伤亡。据中国地震台网测定,震源位于龙门山断裂带南段,震中位置(31.0°N,103.4°E),震源深度经重新定位后为16 km(黄媛等,2008)。迄今为止,在板内逆冲型地震中,汶川地震地表破裂带结构最为复杂。震后的野外地质考察表明,龙门山断裂带的映秀—北川断裂是这次地震的主发震断裂,形成的地表破裂带长达240多公里,最大垂直位移和水平位移分别达到6.2 m和4.9 m。龙门山断裂带的灌县—江油断裂同时也发生破裂,形成70多公里长的地表破裂带,最大垂直位移达到3.5 m(徐锡伟等,2008)。

陈运泰等(2008)和王卫民等(2008)利用地震波数据反演地震破裂过程,发现这次地震以逆冲为主、兼具少量右旋走滑,但地震波反演缺少对断层破裂宽度等几何参量的约束。Shen等(2009)利用GPS同震水平和垂直形变及Insar视线向的形变数据反演同震断层模型和滑动分布,反演结果表明断层由南向北倾角逐渐变大,东北部倾角近直立,但缺乏近场同震垂直位移约束,并且Insar形变结果受大气影响较大。王敏(2009)利用GPS同震水平和垂直位移作为约束,反演破裂断层的平均宽度在10～18 km,构造了类似铲状的断层模型,但滑动模型较为粗糙,而且也缺乏近场同震垂向形变约束。许才军等(2009)根据GPS同震位移资料采用基于敏感度迭代拟合法反演分析该次地震的精细同震滑动分布,但采用迭代方法不能对滑动量进行误差估计。Wang等(2009)以基于精密水准测量观测的近场同震垂向位移为约束,构造了倾角随指数变化的铲状断层几何模型,计算的地表垂直形变符合水准测量结果的主要形态,但是断层模型单一并缺少水平方向位移的约束。

本文主要利用GPS和高精度近场水准观测资料获取的同震位移作为约束,建立有限断层震源模型反演汶川地震的同震断层几何模型和断层面上滑动分布。

本文采用157个GPS点的同震水平位移和44个GPS点的同震垂直位移数据来源参见国家重大科学工程“中国地壳运动观测网络项目组”(2008)。

1983、1987和1997年中国地震局曾对汶川地震破裂带中段的茂县—北川—绵竹水准路线进行过3期一等水准测量; 国家测绘局在1983年和1986年对位于破裂带北段的北川—桂溪、桂溪—绵阳、桂溪—平武和平武—青川水准路线进行过二等水准测量。汶川地震后,中国地震局第一监测中心和第二监测中心分别于2008年9月和2010年3月对绵竹—北川、江油—桂溪两条水准路线进行了2期一等水准复测,并分别对桂溪—平武段和平武—青川段进行2期一等水准复测,在2009年5月和2010年3月对茂县—北川水准路线进行2期一等复测。由于唐家山堰塞湖的阻挡,震后复测的茂县—北川水准路线和平武—绵竹的长水准路线没有连接起来,从而成为一条孤立的水准路线。

利用水准数据估计同震垂直位移,一般先计算震前各水准点相对于基准点的高程,再计算震后复测的水准点相对于震前同一基准点的高程,两期的高程之差即为同震垂直位移。但这种做法存在2个方面的不足:(1)如果水准路线过长,系统误差就会累积较大,给观测值带来更大的偏差,不能够真实反映各水准点的同震垂直位移;(2)地震发生后,震前连续的水准路线可能会被破坏,从而出现多个基准的问题,如茂县—北川水准路线。本文从相邻水准点之间的高差出发,以相邻水准点在地震前后高差变化为观测值,因此不再受水准基准点和路线长度的限制,有效解决了汶川地震后水准观测资料中存在多个基准的问题,能够更好地反映地震导致的地表垂直位移变化。表1列出了相邻水准点的高差变化。

汶川地震同震引起的相邻水准点高差变化的观测误差为

σ_L=(L·(σ²_pre+σ²_post))^1/2.(1)

式中,L为相邻水准点之间的距离,σ_pre为震前水准测量每公里中误差,σ_post为震后水准测量每公里中误差。根据国家精密水准测量规范,一等水准测量每公里全中误差为1 mm,二等水准测量每公里全中误差为2 mm。茂县—北川段水准路线震前和震后都是一等水准测量,平武—绵竹和青川—平武水准路线震前是二等、震后是一等水准测量,将相应的每公里全中误差和水准路线长度代入(1)式就可得到相邻水准点高差变化的观测误差。

需要注意的是表1中的高差变化并非完全由汶川地震引起,还包括震前长期地壳形变影响以及震后滑移影响。通过分析茂县—北川—绵竹水准路线震前的3期观测数据,得出龙门山逆冲带与四川盆地之间的长期抬升速率大约为1.5 mm/a,因此1983～2008年期间震前形变不会超过38 mm,若相邻两个水准点之间的距离为5 km,则两水准点间的震前垂直形变为7.6 mm,不会对同震位移造成影响; 根据北川—绵竹以及桂溪—江油在2008年9月和2010年3月的2期复测结果,在这一年半的时间内,位于破裂带中段的相邻水准点高差形变量最大不超过20 mm,由此推测位于破裂带北段的白白线(平武—青川)同震垂向位移不会受震后形变较大的影响。

根据余震序列的重新定位(黄媛等,2008)和地表破裂带分布特征(徐锡伟等,2008),本文建立了一个长度为300 km的映秀—北川主断层模型和80 km的灌县—江油辅助断层模型。主破裂模型中,北川以南(称为北川断裂)的破裂带为铲状模型,北川至青川的破裂带(称为青川断裂)为矩形断层。由于在灌县—江油断裂(称为灌县断裂)与北川断裂之间缺乏观测约束,反演确定灌县—江油破裂断层的全部几何参数和滑移参数存在一定困难。地质考察表明灌县—江油破裂带为纯逆冲断层型破裂,故先验设定该断层起始于(31.62°N,104.37°E)(许才军等,2009),走向为224°,倾角为35°(徐锡伟等,2008)。

北川断裂的铲状断层意为断层的倾角随深度逐渐变缓,断层面成凹面向上的弯曲,其剖面似铲形。断层倾角的表达式采用Wang等(2009)的指数表达形式

δ=δ₀ ×(1-exp(-C/h)).(2)

其中,δ为断层的倾角,δ₀为断层在地表的倾角,C为常数用来调节断层面的曲率,h为断层深度。

表1 相邻水准点高差的同震形变
Tab.1 Co-seismic deformation displacements of height difference between two adjacent benchmarks

本文确定铲形断层的方法如下:将断层沿深度方向平均分成n份,每个子断层的倾角为该段层中部对应的倾角。

根据Okada(1985)半无限空间弹性位错模型,矩形断层面上发生的位错所引起的地表位移是断层深度h、倾角δ以及断层面上的错动量(走滑和倾滑)的函数,因此可建立如下函数模型:

v^k_i=∑^m_j=1f(s_j,d_j,h_j,δ_j)-y^k_i.(3)

其中,i表示观测数据的个数,m表示子断层数目,k表示观测数据的类型,即GPS同震水平位移(东西向和南北向两个分量)和垂直位移,以及水准测量得到的相邻测站高差的同震形变。f为Okada的解析表达式,s_j和d_j为子断层面j上的走滑量和倾滑量,y表示观测值,v为残差。

由(2)式和(3)式可以看出,未知参数包括各子断层面上的走滑量s和倾滑量d、铲形断层在地表的倾角δ₀、常数C及深度h和青川断裂倾角及深度。由于f与s、d是线性关系,而与其他参数是非线性关系。为简化计算,只利用最小二乘原理平差求出走滑和倾滑参量及其中误差,而用搜索的方法确定其他待估参数。从搜索结果中选择具有最小拟合后残差平方和χ²的未知参数作为最佳估值,然后通过F检验方法(Shen等,1993,1994),确定其90%的置信区间。

将(3)式线性化可得

y=Gx+e e～N(0,D_x).(4)

其中,G为设计矩阵,x为各子断层面上的走滑参量和倾滑参量,e为误差向量,D_x为观测值的方差矩阵。

为了避免地表观测到的同震形变数据在空间上分布不均匀,对相邻的两个子断层j和j+1施加以下约束:

x_j-x_j+1=Bx=0+δx δx～N(0,C_x).(5)

其中,B为约束矩阵,δx为误差向量,对角矩阵C_x为模型约束。当C_x=0时,各子断层上的滑动量严格遵循连续性条件; 当C_x无穷大时,各子断层上的滑动量是不相关的。可通过C_x来调节相邻断层面上滑动量的连续性。

联合(4)、(5)两式,利用最小二乘原理可得:

x^{^}=(G^TD^-1_xG+B^TC^-1_xB)^-1G^TD^-1_xy,(6)

解的分辨率可表示为

R=(G^TD^-1_xG+B^TC^-1_xB)^-1G^TD^-1_xG.(7)

本文采用R的对角线元素之和作为定量判断子断层面上的滑动量与邻近子断层错动量平滑程度的指标(王阎昭等,2008; 万永革等,2008)。该指标值越大,对解的分辨率越清楚,验后方差也越小,但解的稳定性较差; 该值越小(即较紧的模型约束),虽然模型稳定性较好(分辨率低),但验后方差较大。

利用不同精度约束水平滑移的先验值,将反演结果分辨率(resolution)和验后方差绘成一条折中曲线(图1)。由图可见,较紧的模型约束虽然模型稳定性较好(分辨率较低),但验后方差较大; 较低的模型约束虽然具有较小的验后方差和较高的模型分辨率,但解的稳定性较差。根据万永革等(2008)采用的方法,综合考虑这条折中曲线,模型约束采用0.41最为合适。

图1 模型分辨率和验后方差的折中曲线
Fig.1 Trade-off curve between model resolution and posterior variance

表2 汶川地震断层模型参数
Tab.2 Fault model parameters of Wenchuan M8.0 earthquake

图2和图3分别给出了GPS测站同震位移和相邻水准点高差形变的模型结果。表2和图4分别是反演得到的断层几何结构及滑动分布。图5所示是每个子断层块上滑动量值及一倍中误差。最终反演的整个断层模型中,断层南部以逆冲为主,北部走滑分量逐渐增加。北川断裂滑动量较大,最大错动量达到7.6 m,对应破坏最严重的北川地区,其次是映秀和清平,错动量为5.6 m和5.3 m。青川断裂最大滑动量为3.0 m,对应南坝地区。灌县断裂逆冲量为1.1 m。模型结果表明,GPS观测数据拟合后的χ²为0.31 m²,水准数据拟合后的χ²为0.33 m²。从图2和图3中可以看出,模型结果和观测数据符合较好,说明断层模型较好的显示了汶川地震同震断层几何结构。

图2 GPS同震位移的观测结果和模型结果
(红线为地表破裂迹线,蓝色箭头为GPS测站的同震水平和垂直位移的观测结果,
洋红和青色箭头分别为水平和垂直位移的模型结果)
Fig.2 Co-seismic displacements observed by GPS and inverted through model(red curve is the trace of the surface rupture; blue arrows are the co-seismic horizontal and vertical displacements observed by GPS; magenta and cyan arrows are the horizontal and vertical displacements inverted through model)

图3 相邻水准点高差变化形变的观测结果和模型结果
(蓝色箭头表示相邻水准点高差变化的观测结果,绿色、洋红和青色箭头表示模型结果)
Fig.3 Deformation of height difference between two adjacent benchmarks observed by GPS and inverted through model(blue arrows are deformation of height difference of two adjacent benchmarks observed by GPS; green,magenta and cyan arrows are the deformation inverted through model)

图4 破裂断层的同震滑动分布
Fig.4 Co-seismic slip distribution of rupture fault

图5 破裂断层同震滑动量的中误差
Fig.5 Errors of co-seismic slips of rupture fault

本文通过反演2008年汶川地震前后的GPS观测资料和1983～2010年的精密水准观测资料得到的汶川地震同震断层模型结构比较复杂,由北川断裂、青川断裂和灌县断裂组成。其中,北川断裂是一个铲形断层,长度为180 km,地表倾角为62°,随指数函数形式逐渐变缓,深度18 km,断层面滑动以逆冲为主; 青川断裂长为120 km,倾角为58°,断层面上走滑分量由南向北逐渐增加; 灌县—江油断裂为纯逆冲断层,倾角35°,长度为80 km,与徐锡伟等(2008)调查的地表破裂长度72 km相近。整个断层模型中最大错动量达到7.6 m,对应破坏最严重的北川地区。反演结果揭示汶川M_S8.0地震的性质主要是逆冲兼右旋走滑,在破裂带南段断层面上的滑动角接近90°,错动以逆冲为主,特别是在发震断层的起始段逆冲分量明显大于走滑分量; 北川断裂的滑动量明显大于青川断裂的滑动量; 在北川附近,断层面上的滑动量最大达到7.6 m,导致该地区破坏最严重; 北川以北的断层面上走滑分量逐步加大,右旋走滑成为断层破裂的主要特征; 地表的错动量基本上都大于断层深部的结果与王卫民等(2008),陈运泰等(2008),张希等(2009),王敏(2009),Shen等(2009),许才军等(2009)的研究结果一致。从图4可以看出,断层深度在18 km,与黄媛等(2008)采用双差算法重新定位余震的深度一致。从图4还可看出,北川断裂深部15～18 km范围内的滑动量为2～4 m,这与Wang等(2011)得出的结论一致。由于部分水准资料是在震后1～2年进行观测所得,所以北川断裂深部较高的滑动量包含了震后形变的因素,但大部分应是同震滑动。断层深部同震滑动特征和北川断裂的铲形结构可能为龙门山地区隆升是由于地壳缩短造成(Wang et al,2011)。

若取地壳介质的平均剪切模量为3.0×10¹⁰ Pa,则反演得到汶川地震释放的标量地震矩为7.62×10²⁰ N·m,得到相应的矩震级为M_W7.9,与美国地质调查局 http://earthquake.usga.gov/eqcenter/eqinthenews/2008/us2008ryan.和美国哈佛大学 http://www.globalcmt.org/CMTsearch. html.给出的震级一致。