位于青藏高原东北缘的海原断裂带是中国大陆重要的活动地块边界断裂带之一. 该断裂带由冷龙岭断裂、 毛毛山断裂、 老虎山断裂、 海原断裂(狭义)等断裂组成，在西边与祁连山主干断裂相接. 在该断裂带上曾发生1920年8.5级大地震(图 1). 研究海原断裂带的活动特征和应力积累状态对该断裂带及其附近地震预测有很重要的意义，而应力的积累则主要与断层滑动速率和闭锁深度有关(国家地震局地质研究所，宁夏回族自治区地震局，1990; Smith，S and well， 2003). 海原断裂带的分段活动特点，前人已进行了大量的研究(冉勇康，邓起东，1998; 徐锡伟等，2007; 张培震等，2003; 刘静等，2007; Gan et al，2007). 一种是通过研究地形地貌的特征和测年确定断层滑动速率. 国家地震局地质研究所和宁夏回族自治区地震局(1990)、 袁道阳等(1997)、 何文贵等(1994，1996，2000)用此方法研究表明，海原断裂带的左旋走滑速率一般在2-6 mm/a之间，小于Lasserre等(1999，2002)用同样方法确定的(8±4)mm/a-(12±4)mm/a的滑动速率量级. 另一种是用大地测量方法来确定. 中国地壳运动观测网络的GPS观测结果表明，海原断裂带南北两侧表现为比较显著的左旋运动特征，其远场位移所揭示的海原断裂左旋活动速率只有5-6 mm/a(王敏等，2003; Thatcher，2007; Zhang et al，2004)，与袁道阳等(1997)、 何文贵等(1994，1996，2000)利用地质地貌方法确定的断层活动速率比较一致. 张希等(2005)用负位错模型反演的滑动速率为1-2 mm/a，与Lasserre等(1999，2002)的结果同样存在很大差异. 产生这种差异的原因主要有以下几个方面:对于用地形地貌方法来说，差异的主要来源为探槽位置差别和测年误差; 对于用大地测量方法来说主要由测量点位稀疏、 观测时间太短而导致的约束不足，如GPS观测站间距在50 km以上，所用数据仅是1999年和2001年的观测结果; 另一种可能的原因是所使用的模型差别. 因此，使用最新观测数据，用新的模型、 方法反演断裂滑动速率是非常必要的. 对于断层闭锁深度的确定更为困难，用地震波成像方法来反演海原断裂带地壳结构(李松林等，2001)，由于在该断裂带上只有一条地震剖面，在现阶段用此方法很难确定断层闭锁深度分段特点. 因此，用大地测量数据反演断层闭锁深度成为确定断层闭锁深度主要途径.

图 1(Fig.1)

图 1 海原断裂带断层和GPS测站分布图. 白色圆为1900年以来MS ≥ 5地震 Fig.1 Map of faults and GPS sites in Haiyuan fault zone White circles indicate MS ≥ 5.0 earthquakes after 1900

Smith和S and well(2003)提出了三维弹性半空间形变分析方法. 该方法与Okada(1985)位错模型的主要区别是:前者将断层远场块体位移转换为地面下一定深度上的体力矢量，而Okada模型则认为地面位移是由断层位错引起的; Smith的模型不是解析解，而是在波数域解弹性方程，然后进行傅里叶逆变换得到空间域的解，该方案在保持计算精度与数值解一致的前提下大大地提高了计算效率，特别适合于断层比较复杂的断裂带; Okada模型的断层面为有限宽度，Smith模型的断层下界可以为无限深，它更适合由非震断层滑动引起的形变模拟. Smith模型已经用于研究美国圣安德列斯断裂系形变分析，取得了满意的结果(Smith，S and well， 2003，2004). 本文应用多个项目的多次GPS观测数据作为反演约束条件，采用Smith三维弹性半空间形变分析方法反演海原断裂带断层滑动速率和闭锁深度.

本文使用了以下3个项目64个观测站的观测数据:① 徐锡伟、 王庆良等依托科技部的重点平台项目、 国家自然科学基金课题、 中法合作项目建立的古浪-永登GPS剖面、 景泰-白银GPS观测剖面1999年、 2004年和2005年3期观测数据，测点均为基岩标志，观测采用三角架、 光学对中，连续观测3天; ② 中国地震局第二监测中心承担的"'十五'数字地震网络建设项目"建立的华藏寺、 沙沟河GPS、 水准、 重力综合观测剖面2005年、 2006年和2007年的GPS观测数据; ③ "中国地壳运动观测网络"区域站1999年、 2001年、 2004年和2007年在该地区的GPS观测数据. ②和③项目中测站均为强制对中钢筋混凝土观测墩，连续观测4天. 所有项目的观测采用Ashtech Z12接收机，扼流圈天线.

如上所述，我们使用了多个项目的多次观测数据. 为了避免由于约束条件及处理方法不同而产生的速度基准及差别，我们用统一的数据处理方案对所有数据进行了重新处理，这样可以避免采用不同的固体潮、 极潮等地球物理模型，以及卫星和接收机的天线相位中心模型所导致的系统偏差(王敏，2008). 具体方案如下:

首先采用GAMIT软件进行单日解计算(King，Bock，2006)，用GRED计算时间序列，以检查基线及坐标重复性，最后用GLOBK计算站点速度(Herring，2006). 其中，GAMIT数据处理时把中国大陆及周边约10-15个IGS站的原始数据与区域站观测数据一起处理. IGS站坐标(x，y，z)约束取(0.05 m，0.05 m，0.05 m)，区域站坐标(x，y，z)统一给定约束(1.0 m，1.0 m，1.0 m)，卫星轨道约束取10^-8.

用相似变换法计算速度，即通过对坐标和速度的平移和旋转实现ITRF参考框架. 这种方法的优点是不会引起网扭曲变形，基准的稳定性依赖于所选取的IGS站的稳定性，如果选取的IGS稳定点的速度线性比较好，基准就比较稳定. 具体做法是:① 用GLOBK把全球H文件和区域H文件组合起来; ② 用GLORG求速度场，选择位于中国周边和欧亚板块的10个速度稳定的IGS站，以ITRF2000提供的站坐标和速度作为先验值，以其误差的2倍作为约束，区域站的坐标约束为1.0 m，速度约束为0.1 m/a; ③ 把由②求得的速度中扣除欧亚板块旋转影响后的值作为区域地壳运动的速度场，速度的误差一般在0.5-3.0 mm/a.

Smith和S and well(2003，2004)提出的半无限空间3-D体力模型见图 2. 在均匀各项同性的弹性介质中，得到直立断层体力矢量(F_x，F_y，F_z)作用于断层下界d₁和上界d₂引起的位移公式，任何深度z处的位移为源部分、 像部分和Boussinesq改正3项之和

图 2(Fig.2)

图 2 3-D 半空间断层模型示意图 断层由上界d2深度延伸至下界深度d1,位移的非连续性通过置于网格中有限宽度的力偶F来模拟 Fig.2 Sketch of a 3-D fault in elastic space

这里

Boussinesq改正为

上式中，为拉梅常数，i为虚数单位.

通过求点源在断层法方向倒数来设置一个力偶，在实际应用中用高斯函数的导数来近似计算由跨断层引起的应力不连续，由此产生一个有限厚度的断层. 断层埋置在经傅里叶变换后的二维网格中，与上述转换函数相乘，再进行傅里叶逆变换就得到了网格每个结点上的速度矢量.

在该模型中有两种模式:一种是浅部模型，即断层滑动发生在浅部矩形断裂面内，这种模型与Okada(1985)模型一致; 第二种是深部模型，这种模型假设断层下界面为无穷深，它与负位错模型相似(图 3).

图 3(Fig.3)

图 3 断层走滑位移随断层距离分布 浅部模型的位移达到最大值后随距断层距离增大而变小,深部模型的位移达到最大值后随距断层距离增大变化而很小 Fig.3 Variation of displacement with distance away from a strike slip fault

在进行断层参数反演时，采用遗传算法. 遗传算法是基于生物进化论的原理发展起来的搜索优化算法，该算法有利于避免最终解陷入适应度局部最大处，而能搜索到全局最大处. 遗传算法不需要对模型函数求导，因而被广泛应用于地学问题反演的多个领域. 我们在使用遗传算法计算时初始种群数为100，交叉概率为0.6，变异概率为0.01，目标函数为模型值与观测值残差加权平方和，具体如下式:

式中，N为观测站数的2倍，vⁱ_GPS和vⁱ_m分别为第i号点的GPS观测速度分量和模型速度分量.

适应度函数为

根据已有的研究结果(邓起东，1982; 国家地震局地质研究所，宁夏回族自治区地震局，1990)，把海原断裂段分为5段，由西向东依次为毛毛山断裂、 老虎山断裂、 海原断裂西段、 海原断裂中段和海原断裂东段(图 1). 断层的位置走向等参数由MapSIS 2.0中的中国活动断裂数据库中得到，在反演中把它们作为已知参数，断层下界深度d₁取无穷大; 把每一段的断层闭锁深度d₂、 走滑速率和垂直于断层的速率(缩短速率)作为模型参数. 从海原断裂带GPS水平运动图象看(图 4)，存在明显的顺时针旋转. 因此，我们把旋转轴位置及旋转角速度也作为待定参数，这样反演的参数总共为18个. 计算中，剪切模量均为4.12×10¹⁰ Pa，泊松比取0.25. 为了克服闭锁深度与滑动速率的相关性，我们先把根据地质方法确定的断层滑动速率作为固定值，只反演闭锁深度，然后将反演的闭锁深度值加50%的变化作为闭锁深度的初始种群范围，走滑速率初始种群范围为0-20 mm/a，拉张速率为0-5 mm/a. 通过20 000次的迭代，残差均方根为1.1 mm/a，小于速率分量的观测误差.

图 4(Fig.4)

图 4 海原断裂带GPS水平运动观测值与模拟值比较 Fig.4 Comparison between the simulated (void arrow) and observed (solid arrow) rates in Haiyuan fault zone

初步反演结果显示，老虎山断裂走滑速率为10.5 mm/a. 这一结果与何文贵等(1994)用地质方法的结果(4.1-5.5 mm/a)相差太大，这主要是由于位于老虎山断裂以北的GPS测站太少导致约束太弱. 对此我们把老虎山断裂走滑速率的初始范围设为0-7 mm/a，反演结果为6.5 mm/a. 反演最终结果见表 1和图 4. 由观测速度场与模型速度场对比来看，总体大小和方向比较一致，而东部速度在NS方向上差异稍大，这可能是由于断层模型与实际存在一定差异，实际断层面并非绝对直立，而所用模型假定断层为直立走滑断层.

用遗传算法反演，不能直接给出反演参数的精度. 由于我们使用的3-D体力模型的傅里叶解在2维、 直立、 走滑断层假设下，断层中部位置剖面上的位移与Weertman(1964)的代数模型结果一致(Smith，S and well， 2003)，因此，我们从Weertman的代数模型出发，来近似估计反演参数的精度.

Weertman模型的位移表达式为

式中，v为地表观测速度，v₀为断层水平走滑速率，x为观测点距断层的水平距离，w_u为 断层上界面深度，w_l为断层下界面深度. 我们使用的模型中假设断层下界面深度w_l为无穷深，所以式(11)变为

上式两边对反演参数v₀和w_u微分

根据误差传播定律得(武汉大学测绘学院，1982)

式中，m_v，m_v0，m_wu分别为v，v₀和w_u的误差. 由于

所以

同理可得

我们取v₀=30 mm/a，距离x取中等距离50 km，w_u取中等闭锁深度5 km，m_v取拟合残差均方根1.1 mm/a，得断裂滑动速率误差m_v0=2.7 mm/a，闭锁深度误差m_wu=3.2 km. 因此反演结果是可信的.

毛毛山断裂. 毛毛山断裂西部与北祁连山活动断裂带相连，东端与老虎山断裂相接，长约40 km，总体走向N70°W. 大约自中更新世以来其断裂力学性质由压性转变为以左旋走滑为主的压扭性. 反演结果表明(表 1)，毛毛山断裂左旋走滑运动速率为3.5 mm/a，为5个断裂中最小的. 这与何文贵等(1996)利用地质地貌方法得到的全新世滑动速率3.69 mm/a，袁道阳等(1997)利用黄土剖面古土壤年龄方法得到的结果(2.3-3.9 mm/a)非常一致. Lasserre等(1999)同样用地质地貌方法计算的滑动速率为(12±4)mm/a，我们的结果与其相差较大. 这可能有两个原因:首先Lasserre的结果是毛毛山-老虎山断裂平均结果，并未分段; 其次是由于研究的地点和方法有所不同. 该段有跨断裂、 站间距为5 km左右的华藏寺GPS剖面作为约束，反演结果比较可靠.

老虎山断裂. 老虎山断裂带位于青藏高原隆起区的东北缘，该断裂主要形成于加里东期，后经多次构造变动，在中更新世断裂的力学性质由挤压逆冲转化为左旋走滑. 反演结果显示，老虎山断裂左旋走滑速率为6.5 mm/a，与Lasserre等(1999)根据地形地貌计算的毛毛山-老虎山断裂的滑动速度(12±4)mm/a差异较大，与甘卫军等(2005)根据GPS资料正演的结果(8.3 mm/a)比较接近. 何文贵等(1994)根据最新测年资料，求得老虎山断裂中更新世中期以来的水平滑动速率为2.4-2.8mm/a，中更新世晚期以来为3.65-4.17 mm/a，晚更新世早期以来为4.1-4.8 mm/a，最大为5.5 mm/a，该断裂全新世以来滑动速率具有明显加快的趋势. 与何文贵等结果相比我们的结果稍偏大，这种差异一方面可能是由于方法不同，反映的断裂滑动的时间尺度也不同. 后者反映的是全新世以来1万年尺度的平均速度，而用GPS观测数据反演得到的是断裂几年至几十年的滑动速度. 另一方面，在该断裂附近，GPS观测站较少，反演结果的可靠性有所降低.

海原断裂带(狭义). 海原断裂带是1920年海原8.5级地震产生过破裂和错动的断裂段. 其西段西自景泰东至大营水，长约80 km; 中段西自大营水东到干盐池盆地，位于南华山与西华山之间，长约40 km; 东段西起干盐池经南华山东至固原以西长约110 km. 其晚第四纪活动以左旋走滑为主要特征，1920年海原8.5级大地震产生了10.5 m的最大左旋位移(张培震等，2003). 反演结果显示，海原断裂带西段、 中段和东段滑动速率依次为4.0 mm/a、 5.6 mm/a和5.5 mm/a. 西段最小，中段和东段的滑动速率相当，这与国家地 震局地质研究所和宁夏回族自治区地震局(1990)、 甘卫军等(2005)的结果(2.9-3.3 mm/a，4.1-5.4 mm/a，5.8-7.7 mm/a)比较一致. 中段和东段与Ding等(2004)的结果(6 mm/a)很接近，也接近Lasserre等(1999)结果(8±2)mm/a的下边界值. 但与张希等(2005)的结果(1.1-1.7 mm/a)相差很大. 他们仅用了网络工程区域站1999-2001年的GPS观测数据，观测站间距平均在50 km以上，也未使用GPS剖面观测数据. 综上所述，在海原断裂带(狭义)上，大多数的研究结果比较一致.

反演结果还显示(表 1)，毛毛山断裂闭锁深度为22.0 km，为5段断裂中最深的一段; 老虎山断裂闭锁深度为10.3 km; 而海原断裂带(狭义)西段、 中段东段闭锁深度依次为8.4 km、 3.6 km和4.3 km. 这与该地区地壳厚度沿该断裂带由西向东逐渐变小的变化趋势基本一致. 老虎山断裂20世纪以来发生的最大一次地震为1990年天祝-景泰6.1级地震，其深度为12 km. 我们反演的闭锁深度与其基本一致. 1920年海原8.5级地震的震源深度一般认为在17-20 km，而我们反演的闭锁深度与之相比要小得多. 这可能是由于海原8.5级地震使断裂贯通比较彻底，在地壳下部还没有完全形成闭锁. 由于闭锁深度浅，不能形成大量应变能积累，这可能是海原断裂带(狭义)自1920年海原8.5级地震后88 a来未发生6.0级以上地震的主要原因.

表 1(Table 1)

表 1 海原断裂带断层速率与闭锁深度反演结果 Table 1 Inversion results of fault slip rates and locking depths 序号 断裂名称 起点坐标 终点坐标 长度/km 闭锁深度/km 走滑速率/mm·a-1 拉张速率/mm·a-1 λE/° φN/° λE/° φN/° 注:表中右边一列数值前的负号表示压性. 1 毛毛山断裂 102.1 37.6 103.3 37.2 37.37 22.0 3.5 -0.4 2 老虎山断裂 103.3 37.2 104.1 37.1 68.01 10.3 6.5 -0.2 3 海原断裂西段 103.9 37.1 104.8 36.9 80.56 8.4 4.0 -0.3 4 海原断裂中段 104.8 36.9 105.2 36.7 37.60 3.6 5.6 -0.1 5 海原断裂东段 105.1 36.8 106.1 36.3 101.9 4.3 5.5 -0.6

表 1 海原断裂带断层速率与闭锁深度反演结果 Table 1 Inversion results of fault slip rates and locking depths

1)用基于力偶的3-D体力模型，通过傅里叶变换模拟任意深度的位移场、 断层数目和断层迹的复杂性对计算量影响很小，在保持计算精度与数值方法一致的前提下大大地提高了计算效率. 其缺陷是只适用于直立型断裂.

2)毛毛山断裂左旋走滑运动速率为3.5 mm/a，闭锁深度为22.0 km; 老虎山断裂左旋走滑速率为6.5 mm/a，闭锁深度为10.3 km; 海原断裂带西段、 中段和东段的滑动速率依次为4.0 mm/a、 5.6 mm/a和5.5 mm/a，闭锁深度依次为8.4 km、 3.6 km和4.3 km.

3)毛毛山断裂左旋走滑运动速率小，但闭锁深度大，有利于积累应变能. 该断裂也是6级以上地震的空区. 因此该断裂及附近地区具有发生强震背景条件. 老虎山断裂左旋走滑速率为5段中最大的，其闭锁深度为10.3 km，有利于中强地震的孕育. 海原断裂带中西段，断裂闭锁深度小，不利于能量积累.