张志斌 梁晓峰 周贝贝 刘代芹 唐明帅

1)新疆维吾尔自治区地震局， 乌鲁木齐 830011 2)中国地震局地球物理研究所， 北京 100081 3)中国科学院地质与地球物理研究所， 北京 100029

0 引言

印度板块持续向N挤压被认为是东亚地区一系列变形的主要驱动机制之一(Zhangetal.， 2012； Chatterjeeetal.， 2013)； 而天山造山带作为世界上最年轻(年龄<20Ma)与最宏伟(海拔>7000m)的陆内造山带之一， 经历了显著的以褶皱逆冲构造发育为特征的地壳缩短变形(Tapponnieretal.， 2001； Leietal.， 2013)， 是中国内陆地震活动最为强烈的地区之一。北天山中段地处天山与准噶尔盆地的交会区， 是新疆人口和经济分布最为集中的地区， 同时也是丝绸之路经济带的重要窗口。在印度板块和欧亚板块碰撞远程效应的影响下， 该区山前发育一系列活动断裂， 从南向北依次为亚马特断裂、 准噶尔南缘断裂、 博格达弧形断裂带和霍尔果斯-玛纳斯-吐谷鲁断裂(图 1)。这些褶皱逆冲构造带上曾发生数次大地震(Molnaretal.， 1975； Avouacetal.， 1993； 邓起东等， 1999)。自1600年以来， 该区域共记录到5级以上地震20次， 6级以上强震5次， 最大的为1906年玛纳斯7.7级地震； 最近一次为2016年呼图壁MS6.2 地震(下文简称呼图壁地震)， 该地震发生在北天山山前逆断裂带上， 发震断层表现为逆冲特征， 在震中附近未发现地表破裂。由于目前对大陆内部地震发震机制的认识不足， 且对该地区深部结构特征缺少针对性研究， 导致在了解该地区地震孕震机制方面存在一定困难， 至今对呼图壁地震发震断层也未有统一的认识(Lietal.， 2018； 刘建明等， 2018； Luetal.， 2018； 杨文等， 2018)。

图 1 研究区主要断裂和所使用台站及地震震中分布图Fig. 1 Distribution of major faults， stations and earthquake epicenters in the study area.图中红色实线为断层： F1亚马特断裂； F2准噶尔南缘断裂； F3博格达弧形断裂带； F4霍尔果斯-玛纳斯-吐谷鲁断裂。紫色圆圈为与储气库的距离≤30km的范围， 震源机制解为研究区1976—2018年MS>5.0地震， 数据来自GCMT。右上小图为2016年1月28日阜康 MS3.7 地震在50km震中距范围内台站所记录的地震波形， 震源深度为30km

地球物理探测是获取研究区深部构造和孕震环境的有效手段。近年来， 不同学者在新疆天山中段开展了接收函数、 体波层析成像和人工震源剖面(莎雅—布尔津和乌鲁木齐—库尔勒测线)等大量地球物理探测， 获取了该区地壳上地幔地震波速度结构的基本特征， 初步探明了天山与准噶尔盆地及塔里木盆地在较大尺度上的盆山耦合关系， 探讨了天山的隆升机制和深部变形特征， 发现地幔物质上涌对天山隆升所发挥的重要作用， 并提出了 “岩石圈拆沉”和 “层间插入消减”等地球动力学模型(胥颐等， 2000a； Zhaoetal.， 2003； 郭飚等， 2006； Leietal.， 2007； Heetal.， 2014， 2018)。

前人的研究成果为认识天山的隆升机制提供了重要深部约束， 但这些研究重点关注岩石圈结构， 北天山—准噶尔盆地过渡区较为精细的地壳三维结构较为缺乏， 这限制了对天山北麓深浅耦合关系及破坏性地震发震机制的认识。近年来， 新疆测震台网数字网络化的建设及震后流动台显著改善了天山中段的地震台站分布， 为进行地壳结构成像提供了可靠的数据保证。本文利用区域内台站记录的近震走时数据， 采用近震走时层析成像(Simul2000)(Linetal.， 2012， 2014； Douillyetal.， 2016)反演了新疆北天山中段地壳三维速度结构， 进而探讨该区域地震活动性与深部结构的关联性， 为新疆天山地区的孕震机制及构造演化动力学研究提供更多约束。

1 数据与方法

1.1 数据资料的选取

新疆测震台网在经过数字化网络改造后， 对乌鲁木齐周边地区的监测能力得到了一定提升。此外， 为监测呼图壁地震的余震和乌鲁木齐周边的地震， 新疆测震台网在这2个区域分别架设了2套流动台(图 1)， 区内共计18个台站为本研究提供了丰富的近震到时数据。本研究选取了新疆测震台网2009年1月—2018年11月期间区内所记录的917个MS≥1.5的地震事件。为保证反演结果的准确性， 对917个地震事件的P波和S波到时进行了人工检查和拾取， 舍弃了记录不清晰的震相， 并对重新拾取的震相采用Hyposat定位方法进行了重定位， 定位速度模型采用 “3400走时表”解算出来的模型(陈向军等， 2014)。由于采用Pg和Sg的走时作为反演地壳三维速度结构， 为避免Pn震相干扰， 只选择震中距≤200km的到时数据。同时， 按照每个地震事件至少被5个台站记录到且方位角覆盖>180°的要求进行筛选， 最终得到629个地震事件， 其中包括5238条P波和2144条S波到时数据(图2a)。S波的射线条数明显较P波射线条数少。从定位结果分析， 该区域存在一定数量震源深度在中下地壳范围内的地震事件(图2c)， 图 1 展示了研究区内震中距≤50km的台站所记录的1次震源深度为30km的地震波形。分析走时数据可以发现(图2a)， 震中距较近的走时分布较为零散。为此， 采用 “3400走时表”解算的速度模型分别计算了震源深度为10km、 20km、 30km和40km的理论走时曲线， 发现随着震源深度的加深， 震中距较近时走时曲线之间的差异越大， 即本研究中使用的走时数据在近震中距端分布范围较宽的现象反映了研究区内存在部分震源深度较深的地震事件(图2a)。

图 2 走时数据分布及初始速度模型Fig. 2 Travel time data distribution and initial velocity model.a 研究中选用的地震震相数据及不同震源深度的理论走时曲线， 其中红线、 黄线、 紫线和蓝线分别代表震源深度在10km、 20km、 30km和40km情况下的理论走时曲线； b 评估波速比平均变化的和达直线图； c 研究中使用的地震震源深度分布直方图； d本研究使用的初始速度模型

地震层析成像所选用初始速度模型综合了邵学钟等(1996)沿乌鲁木齐—冰达板—库尔勒布设的测深剖面和Zhao等(2003)沿阿勒泰—奎屯—库车沿线人工爆炸宽角剖面获取的P波速度结构(图2d)。初始波速比选择为1.72(图2b)， 该值是基于和达直线法拟合走时数据得出的。该计算波速比的方法在前人的研究中已被广泛验证(Wuetal.， 1987； Langinetal.， 2003； Zhangetal.， 2005； Zhangetal.， 2011)。反演速度模型中的0km是相对于海平面的， 因此本文后续讨论所涉及的震源深度均为相对于海平面的深度。

表1 “3400走时表”速度模型Table1 Crustal velocity model of the“3400 travel time table”

1.2 层析成像方法的原理

Simul2000是目前被广泛使用的近震层析成像反演程序(Thurber， 1983， 1993； Eberhart-Phillips， 1990； Thurberetal.， 1999)， 本研究采用该程序同时反演VP和VP/VS速度模型， 反演中应用一种全矩阵阻尼最小二乘法， 适用于近震和可控源地震事件。在反演迭代时， P波和S-P的走时残差可表达为VP、VP/VS和地震位置改变量的方程：

(1)

(2)

(3)

2 建立反演模型

2.1 反演参数的选择

采用Simul2000反演时需对网格间距、 初始速度模型和阻尼值等参数进行选择， 其中网格间距根据研究区域和台站分布进行设置。经过多次测试后， 选择在水平方向由间隔为40km的均匀网格构成， 在垂直方向为-5km、 0km、 5km、 10km、 15km、 20km、 25km、 30km和40km网格节点组成的三维网格节点。笛卡尔坐标系的中心选择为(44.5°N， 86.5°E)， 为适应研究区的构造走向， 将模型坐标进行了旋转，x轴和y轴的正向分别指向 W16°N 和 N16°E。

由于层析成像反演问题本身的混定特性， 阻尼因子的选择对反演结果的稳定性影响较大。为使反演结果稳定， 通过一系列大范围阻尼因子的单次迭代反演， 并分析反演后数据拟合残差和模型方差的折中曲线来选择合适的阻尼因子， 以避免只反演部分折衷关系曲线而影响获取真实阻尼因子(1)User’s manual for SIMULPS12 for imaging VP and VP/VS： A derivative of the ‘Thurber’ tomographic inversion SIMUL3 for local earthquakes and explosions。。具体做法为： 利用P波和S波数据， 首先将VP/VS的阻尼值固定为一个较大值(9999)， 以将S波的数据影响控制为最小； 根据折中曲线选择VP的阻尼值， 由图3a可知， 此时选择700作为VP的阻尼值； 之后， 将VP的阻尼值固定为700， 再进行系列反演得到折衷曲线， 选择VP/VS的阻尼值为600(图3b)。为验证700是否为合适的VP阻尼值， 再固定VP/VS的阻尼值为600， 进行一系列大范围的VP阻尼值单次迭代测试， 发现将VP的阻尼值设定为700较为合适。

图 3 确定反演中VP及VP/VS阻尼因子的折衷曲线Fig. 3 Determination of the trade-off curves of VP and VP/VS damping factors in inversion.a 确定VP反演的阻尼因子； b 确定VP/VS反演的阻尼因子

2.2 模型分辨率

为评估反演模型的可靠性及空间分辨率， 本文利用 “棋盘格”方法测试层析成像反演结果能达到的分辨率。“棋盘格”分辨率测试方法(Linetal.， 2010)将不同反演节点的VP和VP/VS设置为空间 “正负相间的棋盘格”， 其扰动值分别为初始值的±5%。在分辨率测试中， 震源位置、 台站位置、 初始模型和反演参数与实际数据反演保持一致。为检测计算结果的稳定性， 在计算理论走时时加入标准差为0.2s的随机走时误差。图 4 和图 5 分别为研究区内初始和反演后的VP和VP/VS速度模型及其深度剖面。DWS(Distribution of Derivative Weight Sum)描述了射线长度及射线到每个节点的距离加权， 反映了相对射线密度(Evansetal.， 1994)。测试显示，DWS>100的区域能得到较好恢复， 因此最终结果中只展示DWS>100的区域。此外， 当减小阻尼值时， 走时数据的拟合会更好， 可显示更多细节。但为了保证与反演结果保持一致， 仍使用与最终反演结果一致的阻尼因子来评估模型的分辨率。

图 4 P波及VP/VS速度模型 “棋盘式”分辨率测试及结果Fig. 4 Results of checkerboard test for VP and VP/VS ratio.

图 5 4条剖面VP及VP/VS速度模型 “棋盘式”分辨率测试及结果Fig. 5 Result of checkerboard test for VP and VP/VS ratio at 4 sections.AA′、 BB′、 CC′和DD′的具体位置见图1

“棋盘格”测试显示， 研究区内中心部位的P波速度模型在5～30km深度之间存在较好的分辨率(图 4)； 同时由于浅层地震数据较少(图2c)， 导致5km以浅区域的分辨率受到了限制。分辨率最好的区域主要集中在北天山中段东侧山前区域， 即天山与准噶尔盆地的接触部位， 这与地震分布特征一致。而北天山中段西侧由于地震个数较少， 该区域的输入异常基本无法恢复。而VP/VS速度模型则不同， 由于研究区浅层结构复杂， S波到时拾取存在困难， 导致反演过程中使用的S-P走时差数据较少， 因此VP/VS速度模型的分辨率有限。在5～20km的深度范围内， 由于地震活动性较高，VP/VS速度模型有最佳的分辨率(图 4)。在深度剖面上， P波速度模型在10～20km的深度范围内分辨率最好， 同时也在30km能得到较好的分辨率， 这与地震深度分布相关， 射线可集中穿过这些区域。而与P波速度结构相比，VP/VS速度结构在AA′剖面的分辨率较好， 而在BB′剖面的分辨率较差， 这主要是受限于地震震中分布及台站方位(图 5)。

3 地震重定位结果

采用近震走时反演三维速度模型时， 速度结构和震中位置之间总是存在一定的折中， 可通过多次迭代修正速度模型和地震位置的偏差(Shearer， 2009)。首先， 根据初始的速度模型和Pg及Sg-Pg的走时反演初步的三维速度模型； 然后， 根据初步的三维速度模型和走时进行重定位， 得到重定位的震源位置并重新计算走时。循环上述步骤， 直到震源位置不再改变， 且走时残差均方根达到最小。在每次迭代前和迭代后计算震中位置的改变量， 迭代最后震源位置的改变量平均约为0.5km(即在水平和垂直方向的改变量约为0.3km)(图6a)。在经历25次迭代后， 走时残差稳定在0.52s左右(图6b)， 重定位震源位置的理论偏差约为0.5km。

图 6 迭代前、 后震源位置的改变量及走时残差值的均方根Fig. 6 The root-mean-square of the source position change and travel time residual value before and after iteration.

图 7 给出了层析成像反演前、 后的走时残差均方根直方分布图。图中， 灰色直方分布图显示在层析成像反演前残差值绝大多数分布在±1s的范围内， 透明矩形显示层析成像反演后残差值集中分布在±0.25s范围内， 到时残差均方根由0.8s减少至0.5s， 减少幅值为40%， 地震定位精度较第1次迭代时得到明显提升。

图 7 层析成像反演前、 后走时残差分布直方图Fig. 7 Histogram of travel time residual distribution before and after tomography inversion.灰色矩形为反演前走时残差分布直方图； 透明矩形为反演后走时残差分布直方图

图 8 给出了重定位前后与地震目录震源位置对比图， 重定位后震源深度分布相对零散， 但局部有明显的成簇特征， 震源深度分布的范围拓宽到整个地壳范围内， 剖面AB揭示了天山向博格达弧形构造的地震震源深度下界在变深。历史地震分布表明， 天山地区的地震主要集中发生于天山与塔里木盆地、 准噶尔盆地和哈萨克地台的接触带上， 在山间也分布着一些与断层活动有关的地震， 但天山边缘的震源深度下界明显更深。地震的发生往往被认为是岩石发生脆性破裂所导致的， 但下地壳和上地幔被认为是塑性的， 现阶段该区域深度范围内的发震机理还未形成统一认识， 部分研究者认为高度变形的天山与地壳薄且内部无变形的块体接触， 准噶尔盆地基底的俯冲和现今的碰撞构造可能涉及到天山造山带的下地壳， 导致发生天山前陆深地震事件(Zhaoetal.， 2003)。

图 8 地震目录与重定位后的震中位置分布图Fig. 8 Earthquake catalog， and distribution of epicenters after relocation.a 重定位前震中位置分布图； b 重定位后震中位置分布图

4 反演结果与讨论

图 9 展示了5～30km处每层P波相对于平均速度的扰动量， 由于地表附近地震活动较为稀疏， 且无人工地震事件参与反演， 致使分辨率在浅层较差， 故未展示地表附近的速度模型。在各层速度结果上分别投影了其深度±2.5km范围内的地震。

图 9 5km、 10km、 15km、 20km、 25km和30km深度处P波相对于每层平均速度的扰动量Fig. 9 Map of P wave perturbation relative to the average velocity at depth of 5km， 10km， 15km， 20km， 25km and 30km.

在5km深处， 由于天山造山带的沉积厚度小， 大部分地区基岩出露地表， 准噶尔盆地沉积厚度大， 盆地内沿昌吉一直到沙湾都表现为低速异常。VP速度模型中最显著的特征是沿着天山造山带的高速异常， 比该层速度的平均值(5.801km/s)高约5%～10%(6.09～6.38km/s)， 显示出较好的盆山分界线。随着深度的变化， 在分辨率较好的10km处高速异常区逐渐变宽， 向准噶尔盆地扩散， 以呼图壁为界， 以东表现为高速异常， 以西为低速异常， 呈现明显的速度分界面。在15km的深处， 呼图壁以东表现为相对低速异常， 并逐渐向乌鲁木齐方向扩散； 以西表现为高速异常， 在此范围内以呼图壁为界， 区域两侧的P波速度在横向上存在较大梯度。

在20km深处， 受分辨率限制， 成像区域缩小， 区内主要以一条SN向相对高速异常为主， 同时在呼图壁地震震中位置即准噶尔南缘断裂带上表现为高速异常。高速异常往往与地震活动有关， 通常P波高速异常表明介质具有相对较高的强度， 不易发生破裂， 承载力更强的岩石能够维持地震应力并在地震事件中释放(吴建平等， 2009)， 这为强震的孕育提供了载体。该深度处沿准噶尔南缘断裂有较强的地震活动性， 地震集中分布在该断裂的两侧。在准噶尔南缘断裂10～20km深度处地震活动最为集中， 同时该处也是高、 低速过渡带， 区域介质变化比较剧烈， 地震活动密集。

在25～30km深度范围内， 区内以低速异常为主， 分布于准噶尔南缘断裂及博格达弧形断裂附近， Xin等(2019)和蔡妍等(2019)的研究成果显示在该区域中下地壳存在低速特征。地质资料表明， 天山造山带的热流值是准噶尔盆地和塔里木盆地的2倍以上(Zhaoetal.， 2003； Leietal.， 2007； Pirajno， 2010)， 且中下地壳存在广泛分布的低密度体， 电阻率变化较大。结合构造成因分析， 大型韧性剪切系统往往表现为低速体， 且在浅层脆性破裂机制控制下的大型断裂向下延伸， 一般与中—下地壳的低速体相连(杨晓平等， 2002)， 该区域展布的低速带分布在准噶尔南缘断裂两侧。准噶尔南缘断裂作为天山和准噶尔盆地的构造边界， 其埋深约为25km， 现今构造活动强烈。此结果表明天山造山带中下广泛分布地壳低速体， 可能为天山内部脆韧性过渡带的底部。

受限于S波观测数量少，VP/VS成像结果只在5～20km之间呈现较好的分辨率， 故只展示5～20km范围内的成像结果。研究区内VP/VS在深度为5km时以乌鲁木齐附近的高波速比异常为主， 其值约达2.0， 在20km深度变为1.6～1.7(图 10)。其中从15km深度开始， 以博格达弧形构造带附近的低波速比异常为主， 且该特征一直延伸到20km深度处都比较清楚。而博格达弧形构造因博格达峰位于其附近而得名， 是研究区东部天山与准噶尔盆地的分界线(罗福忠等， 2006)， 其中上地壳呈现出相对低波速、 低波速比等特征， 指示该区域地温梯度较低， 构造变形强烈， 介质比较破碎。著名的博格达峰是天山的第三高山峰， 博格达弧形构造中上地壳的低波速和低波速比异常及中下地壳地震多发的特征是否与博格达山峰的隆起相关， 值得进一步研究。

图 10 5km、 10km、 15km和20km波速比(VP/VS)水平切片Fig. 10 Map of VP/VS ratios at depth of 5km， 10km， 15km and 20km.

为分析研究区内中小地震的活动性与速度结构特征的联系及呼图壁地震发震构造背景， 在呼图壁地震周边分别沿着断层走向和垂直于断层走向绘制了4条P波速度剖面(图 11， 位置如图 1 所示)， 并对分辨率较好的区域进行成像， 同时将发生在剖面两侧20km范围内的地震垂直投影在相应剖面上。

图 11 沿图1中4条深度剖面的P波速度结构图像Fig. 11 P wave velocity structure beneath the 4 sections shown in Fig. 1.

AA′剖面横穿整个研究区域， 与准噶尔南缘断裂走向一致， 该区域范围内的地震主要发生在速度梯度变化较大的区域内， 在低速体和高速体内部地震数量较少。重定位后地震事件主要集中在约15km深度附近， 其中呼图壁MS6.2 地震发生在高、 低速体过渡带的上方。天山造山带在SN向水平挤压和地幔形变力的共同作用下， 刚性的天山上地壳可能沿着中地壳的韧性剪切带发生滑脱， 在相邻块体的阻挡下不得不向两侧盆地发生推覆逆冲， 造成脆性岩体破裂(胥颐等， 2000b)， 使得大部分地震， 尤其是5级以上中强震主要发生在中地壳介质速度变化较大的区域。

BB′剖面由博格达弧形构造穿越进入准噶尔盆地南缘， 经过呼图壁储气库。在乌鲁木齐—玛纳斯一带下方约5km表现为明显的低速异常； 石河子与玛纳斯下方15km的位置呈现局部高速体； 同时， 在博格达弧形构造下方的上地壳存在局部高速的特征。其中可明显看出在呼图壁储气库周边存在自地表到中上地壳的地震孕震带， 该条带在延伸到地下约10km深处。Tang等(2018)的研究表明， 在呼图壁储气库30km的范围内， 存在因注气和抽气导致的一系列诱发地震， 诱发地震的震源最深约为10km， 因此该区的地震活动性与呼图壁储气库的注气及抽气有重要关联。

CC′剖面为SN走向， 横穿了呼图壁储气库和准噶尔南缘断裂。在43.7°N、 20km深的位置， 高、 低速体的交界带上地震活跃， 该地震活跃带即处于准噶尔南缘断裂带上。

DD′剖面横穿呼图壁MS6.2 地震震中。已有相关研究表明， 呼图壁MS6.2 地震为一次逆冲型地震事件， 其发震断层面的走向为292°， 倾角为62°， 滑动角为80°； 双差定位的结果显示， 本次地震为一次N倾的地震， 倾角约为45°～60°(刘建明等， 2018)。DD′剖面显示， 呼图壁MS6.2 地震及其余震发生在S倾的断面上， 且本次地震的余震基本都位于主震上方， 初始破裂点位于最下部， 自下而上产生破裂。2018年5月在呼图壁震中区域发生了一系列震群， 初步判定为呼图壁地震的余震。冯先岳(1997)的研究表明， 准噶尔南缘断裂发育在北天山山前， 长约280km， 具有逆冲性质， 总体走向近EW， 断面S倾， 倾角为45°～75°。而根据本研究在水平和深度方向的投影综合分析认为， 呼图壁地震的发震断层应为准噶尔南缘断裂。本文的定位结果与双差定位结果不一致， 主要可能是受到重定位使用的三维速度模型影响； 同时本研究对震相到时进行了重新拾取， 舍弃了记录不清晰的地震事件， 以保证反演结果的可靠性。这2点因素可能是造成余震定位结果差异的主要原因。

5 结论

本研究利用新疆测震台网固定与流动台站近10a记录的917个地震事件的到时数据， 采用近震体波走时层析成像方法反演了新疆天山中段较为精细的三维速度结构， 给出了利用三维速度结构进行重定位后的震源位置。反演后， 研究区内地震走时残差均方根由重定位前的0.8下降到0.5， 地震震源位置的偏差减小为0.5km， 且重定位后地震的分布更加集中。研究结果显示， 在乌鲁木齐NE侧， 即博格达弧形断裂带区域中下地壳有一定数量的地震发生， 该区域在速度结构上显示出相对较低的P波速度， 同时中上地壳的波速比也相对较低。这些观测表明， 该地区岩石圈的局部温度可能比较低且介质变形强烈， 这有待利用更加精确的下地壳波速比数据进行验证。而获得区域内下地壳地震活动的证据， 有助于我们今后进一步认识该区域的深部构造特征。

此外， 新疆天山中段的速度结构呈现出明显的纵向不均匀性， 浅层受天山造山带快速隆升的影响， 在准噶尔南缘形成较厚的沉积层， 盆山界线明显。在10km深度和15km深度上分别以呼图壁为界存在高、 低速转换带。在中地壳以下， 准噶尔南缘及博格达弧形断裂表现为明显的低速异常， 推测该低速带是天山地壳内部活动的大型韧性剪切带， 该区域历史上发生的如玛纳斯7.7级强震等地震恰好位于该韧性剪切带的上缘(杨晓平等， 2002)， 即新疆北天山中段的中地壳很有可能是天山浅层的脆性变形到地壳深部韧性剪切变形的分界线。

地震活动性与地壳速度结构具有较强的对应关系。研究区内地震主要发生在高速异常体， 高、 低速异常体过渡带及介质结构变化剧烈的区域。呼图壁地震的震源区附近存在局部P波高速异常， 邻近波速比明显变化区域， 该地震可能发生在S倾的准噶尔南缘断裂上， 倾角约为50°， 且余震主要沿着准噶尔南缘断裂带展布。此外， 呼图壁地震震源区周边都表现为低速异常， 而震中附近的高速异常代表介质易产生应力集中和释放， 这为该区域强震的孕育和发生创造了条件。

致谢本研究所使用的地震到时数据来自新疆维吾尔自治区地震局监测中心； 本文的图件使用GMT绘制； Simul2000程序来自Cliff Thurber教授(2)http: ∥www.geology.wisc.edu/～thurber/simul2000/。； 审稿专家对本文提出了宝贵意见。在此一并表示感谢!