余鹏飞,熊维,陈威,乔学军*,王迪晋,刘刚,赵斌,聂兆生,李瑜,赵利江,张怀
1 中国地震局地震研究所,武汉 430071 2 中国地震局地震大地测量重点实验室,武汉 430071 3 湖北省地震局,武汉 430071 4 中国地震局地球物理研究所,北京 100081 5 中国地震台网中心,北京 100045 6 青海省基础测绘院,西宁 810001
据中国地震台网(CENC)测定,北京时间2021年5月22日2时4分,在青海省果洛藏族自治州玛多县(98.34°E,34.59°N)发生MS7.4地震,震源深度17 km.据野外地质调查,此次地震造成了约150~160 km的地表破裂及大面积的砂土液化,多处桥梁、房屋等损毁(李智敏等,2021;潘家伟等,2021).美国地质调查局(USGS)、全球矩心矩张量目录(GCMT)等给出了此次地震的震源机制解(表1、图1),均显示此次地震为一次走滑型事件.
玛多地震位于青藏高原中部巴颜喀拉次级块体内,距巴颜喀拉块体北部边界东昆仑断裂约70 km.
表1 玛多地震的震源机制解Table 1 The focal mechanism solutions for the Madoi earthquake
图1 玛多地震震中及周边地区地形、断层、历史强震、震间GNSS水平速度场及InSAR数据覆盖地区红色五角星为玛多地震震中,红色震源球为USGS及GCMT发布的玛多地震震源机制解,黄色圆圈为重定位的余震(王未来等,2021),黑色曲线为断层(邓起东,2007),黑色震源球为来源于GCMT的历史强震,黑色箭头为陆态网络站点1999—2020年间相对欧亚板块的震间GNSS水平速度场,蓝色方框为ALOS-2卫星升降轨影像覆盖地区,绿色为Sentinel-1卫星升降轨影像覆盖地区.Fig.1 Distribution map of topography, fault, major historical earthquake in the Madoi areasThe red star marks the epicenter of the Madoi earthquake from CENC. The red beach balls are focal mechanism solutions of the Madoi earthquake from USGS and GCMT. The yellow circles represent the relocated aftershocks (Wang et al.,2021). The black lines are faults (Deng,2007). The black beach balls are focal mechanism solutions of the historical earthquake from the GCMT catalog. The black arrows represent the horizontal velocity field relative to the Eurasian plate (from CMONOC, between 1999 and 2020) provided by the GNSS Data, Products and Services Platform (http:∥www.cgps.ac.cn/) of the China Earthquake Administration. Blue box is the coverage area of ALOS-2 satellite image and the green box is the coverage area of Sentinel-1 satellite image.
青藏高原是我国现代构造活动最强烈的地区,而巴颜喀拉块体是青藏高原内部近20年来构造活动最强烈的次级块体(邓起东等,2002,2014)(图1).自1997年以来,中国大陆所有7级以上地震都发生在巴颜喀拉块体边界断裂上,包括1997年玛尼MS7.5地震、2001年昆仑山MS8.1地震、2008年于田MS7.3地震、2008年汶川MS8.0地震、2010年玉树MS7.1地震、2013年芦山MS7.0地震、2014年于田MS7.3地震和2017年九寨沟MS7.0地震(图1).此次玛多地震的发生表明巴颜喀拉块体仍是中国大陆强震活动的主体地区.深入研究此次地震的形变及震源特征,对理解巴颜喀拉块体的变形机制和构造活动特征以及区域地震危险性评估都有重要科学意义.
大地测量技术获取的同震形变是强震发生的直接表现,也是解剖地震的重要资料.以GNSS和InSAR为代表的现代大地测量技术在震源机制等研究中已获得了广泛的应用(Xu et al.,2010;乔学军等,2014;单新建等,2015;李永生等,2016;Du et al.,2018).玛多地震发生后,不同学者利用连续GNSS或InSAR观测等大地测量资料开展了相关研究.李志才等(2021)利用21个连续GNSS测站获取了玛多地震的滑动分布,由于观测数据太少且大部分为远场信息,滑动模型相对粗糙.华俊等(2021)基于Sentinel-1卫星升降轨数据计算了滑动分布及应力扰动,但由于影像失相关在近场缺少约束,近场模拟残差相对较大.Wang等(2022)利用DInSAR技术,结合像素偏移追踪补全近场形变信息,反演的同震滑动最大高达7 m,远大于其他学者的结果.因此本文利用震后及时获取的近场流动GNSS观测,结合GNSS连续观测、Sentinel-1和ALOS-2 InSAR观测,获取玛多地震精细的同震形变场,反演震源参数及同震滑动分布模型,并对区域地震活动性的影响进行初步探讨.
新生代以来,在印度板块以~50 mm·a-1的速度向欧亚板块碰撞和持续汇聚作用下(Molnar and Tapponnier,1975),青藏高原地壳南北向缩短,垂向增厚.在这一动力过程中,青藏高原若干次级构造块体分别沿大型走滑断裂带以不同速率向东或东南方向滑动(Tapponnier et al.,1982,2001;Peltzer et al.,1989;闻学泽等,2011).巴颜喀拉块体是这些块体中最活跃的块体之一.
被一系列强烈活动的构造带所分割和围陷,青藏高原可划分为祁连山、东昆仑—柴达木、巴颜喀拉等6个次级块体.巴颜喀拉次级块体位于青藏高原中部,为一走向北西的长条状块体,主要受其南北两侧的走滑断裂所控制(邓起东等,2010).其北边界为东昆仑断裂,南边界为玛尔盖茶卡-若拉岗日断裂、甘孜—玉树断裂和鲜水河断裂,东边界为龙门山断裂.东昆仑断裂由一系列断裂按一定的几何结构组合而成,各断裂多呈左行左阶羽列组合形式,其间夹有拉分盆地或挤压隆起.东昆仑断裂东起若尔盖以东,与塔藏断裂相接,向西经玛曲、玛沁、托索湖、阿拉克湖、纳赤台、至库赛湖以西,继续西延与阿尔金断裂相交,全长超过1600 km,是一条巨型走滑断裂(青海省地震局和中国地震局地壳应力研究所,1999;邓起东等,2010).走向北西西,各段倾向不一,总体倾向北,倾角55°~85°.阿拉克湖—托索湖—玛沁—玛曲段,断层走向由290°转为310°再转为280°,形成反“S”型构造.在断裂走滑运动的情况下,断裂走向的急剧变化导致挤压隆起形成阿尼玛卿山断隆.沿东昆仑断裂带的滑动速率,在中西段(西大滩—玛沁)为10 ~ 12.5 mm·a-1,向东递减,在玛沁—玛曲段衰减到5~7 mm·a-1,至塔藏段降至3 mm·a-1以下(Van der Woerd et al.,2002;Lin,2008;李陈侠,2009;Kirby and Harkins,2013;Ren et al.,2013).自1900年以来,该断裂发生了6次7.0级以上地震(青海省地震局和中国地震局地壳应力研究所,1999;中国地震局震害防御司,1999;邓起东等,2002),形成了长达1200 km的地震破裂带(李建军等,2017).
现代GNSS观测结果表明,在印度板块的碰撞挤压作用下,高原地壳在南北向缩短的同时还存在东向挤出,并绕东喜马拉雅构造结旋转(Wang et al.,2001; Zhang et al.,2004).其中,巴颜喀拉地块的运动方向约为NE61.45°,速率约为21 mm·a-1(张培震等,2003).王伟和王琪(2008)利用GNSS资料建立的中国大陆活动地块运动学模型中巴颜喀拉块体的运动速率则为16.5 mm·a-1.由GNSS观测获得的垂直运动速度场表明巴颜喀拉块体109国道以东地区均表现为持续的隆升,速率在1~3 mm·a-1,其中龙门山、贡嘎山隆升速率达2~3 mm·a-1(Liang et al., 2013).40余年的精密水准资料也显示龙门山断裂以西的高原地区隆升速率达3 mm·a-1以上(王双绪等,2013).
地震发生后,中国地震局地震研究所迅速组织野外测量与科学考察队于2021年5月23日前往地震现场,开展了二十多天的震后流动GNSS观测.流动观测站点包括陆态网络区域站、测绘局B级网点等.此外,我们还收集了震中400 km范围陆态网络基准站,测绘局CORS站等连续GNSS测站的震前与震后数据.所有观测数据基于GAMIT软件,采用常规精密处理方法,获得了玛多地震水平同震形变场(图2).GPS解算结果,北方向、东方向中误差分别为0.003 m、0.003 m.GNSS形变场显示玛多地震同震形变非常显著,距震中200 km仍有1 cm左右同震形变,断层北侧4454测站水平形变达1.208 m,向西运动;断层南侧2836测站水平形变达0.683 m,向南东东向运动,同震形变场在南北两侧分别表现为顺时针和逆时针旋转的特点,表明此次地震是左旋走滑错动的结果.
欧空局(ESA)的C波段Sentinel-1系列卫星在IW(Interferometric Wide)观测模式下采用TOPS技术进行观测,不仅能够提高影像质量,还能够提高影像的干涉性能.工作于L波段的ALOS-2卫星,由于拥有更长的波长,能够监测到更大的形变梯度,因而在强震的同震形变监测中应用广泛.Sentinel-1和ALOS-2卫星在震后分别对本次震中区进行了及时的应急观测,使我们获得了很好的升降轨InSAR同震形变场(表2,图3).
我们使用瑞士GAMMA软件对所有SAR数据进行处理(Werner et al.,2000).其中,利用美国宇航局发布的30 m分辨率的SRTM数字高程模型进行地形相位的模拟和消除,采用加权功率谱技术对干涉图进行滤波以提高信噪比(Goldstein and Werner,1998),使用最小费用流法进行相位解缠(Pepe and Lanari,2006),采用多项式拟合以消除轨道误差,经过地理编码最终获取玛多地震视线向(LOS)的同震形变场(图3).Sentinel-1升、降轨及ALOS2升、降轨LOS向形变中误差分别为0.012 m、0.015 m、0.034 m、0.046 m.
图2 玛多地震GNSS水平同震形变场Fig.2 The GNSS horizontal coseismic deformation field of the Madoi earthquake
表2 InSAR同震形变数据信息Table 2 InSAR coseismic deformation data information
图3 玛多地震同震形变场(a—c)、(d—f)、(g—i)和(j—l)分别为Sentinel-1升轨、Sentinel-1降轨、ALOS2升轨和ALOS2降轨的干涉条纹图、LOS向形变场及降采样的结果.黑色实线为玛多地震的地表破裂迹线.Fig.3 Coseismic deformation field of the Madoi earthquake(a—c), (d—f), (g—i), (j—l) represent the interferograms, the LOS deformation fields, and the down-sampled results for Sentinel-1 ascend track, Sentinel-1 descend track, ALOS2 ascend track and ALOS2 descend track, respectively.The black lines are the ruptured fault traces of the 2021 Madoi earthquake.
Sentinel-1升降轨影像完整地捕捉了玛多地震的同震形变场,整个形变场呈椭圆状(图3a,3d),形变场长轴方向为105°(或285°).此次地震造成了十分显著的地表形变,最大视线向形变达0.9 m.由于升降轨观测的视线向不一致,升轨和降轨图像上发震断层两侧的形变符号相反(图3b,3e).地震造成的地表形变范围非常大,长宽约220 km×160 km.除靠近地表破裂的位置因失相关导致部分条纹缺失,升、降轨干涉条纹整体上连续光滑.干涉条纹以发震断层为中心,两侧总体呈对称分布.靠近发震断层位置有多个不连续的同心圆,说明位错量沿断层的不均匀分布,推测沿断层走向可能存在多个凹凸体.条纹十分密集,说明具有较大的形变梯度.采用Stripmap模式的ALOS-2升轨影像由于刈幅相对较窄仅捕捉到东侧的同震形变.SacnSAR模式的ALOS-2降轨影像则完整捕捉了同震形变.由于波长更长,升降轨ALOS-2干涉图条纹数量相对较少.整体来看,ALOS2获取的同震形变场同Sentinel-1的结果相似,但ALOS-2降轨影像视线向最大形变达1.2 m.图3显示此次地震造成了明显的地表破裂,西起鄂陵湖,东达昌麻河,长约160 km.各形变场图像上地表破裂带的位置十分吻合,并没有因为卫星视线向不一致而导致偏移.
InSAR同震形变场数据量大,且相邻观测值之间具有很高的空间相关性,将全部InSAR观测值用于反演计算没有必要且会影响反演效率.我们利用四叉树采样法(Jónsson et al.,2002)对InSAR同震形变场进行降采样处理,InSAR观测值数量大幅减少,且原始形变场的主要特征在降采样过程中得到较好的保留.
Sentinel-1和ALOS-2卫星升降轨形变场显示此次地震造成了非常显著的地表破裂,且各形变场图像上地表破裂带的位置几乎重合,这使我们可以勾勒出一条走向NWW-SEE的发震断层地表迹线(图3).余震精定位结果显示在余震区东端出现马尾状分叉特征,升降轨干涉图上也显示了与此相符的形变特征,因此,我们在东端主断层南侧相应位置设置了一条分支断层(图3).据此,我们可以确定主断层和分支断层的空间位置,其走向由断层迹线确定并随迹线变化.不同于主断层,余震精定位的结果显示发生于分支断层的余震位于地表破裂迹线的南侧,表明其可能南倾,因此我们也考虑了其倾向南的可能性.结合InSAR形变特征及余震精定位结果,主断层和分支断层长度分别设定为180 km和30 km,宽度均设定为30 km.断层倾角则参考各机构震源机制解、余震精定位结果、震中区域地质构造特征限定在65°~90°范围.通过网格搜索,最终确定了主断层倾角85°,倾向北,分支断层倾角68°,倾向南.
基于弹性半空间位错模型,以GNSS和InSAR同震形变场为约束,采用基于约束条件的最小二乘原理及最速下降法进行发震断层的同震滑动分布反演(Wang et al.,2013).泊松比设为0.25,剪切模量设为30 GPa.为更好的获取断层面滑动分布细节,将主断层和分支断层按2 km×2 km划分为若干个子断层.在以大地测量数据反演同震滑动分布研究中,对多个数据集的相对权重关系,有经验定权(Feng et al., 2010),模型残差最小定权(申文豪等,2019),观测数据误差定权(He et al., 2021)等多种定权方法,本文根据观测数据误差定权,将GNSS和Sentinel1升、降轨,ALOS2升、降轨观测值相对权重分别设为1、0.25、0.20、0.09、0.07,既反映了各观测数据的精度水平,也是大地测量相关计算的通行做法.最佳光滑因子则由模型粗糙度和相对拟合残差的折中曲线确定为0.07.最后获取了基于均匀介质模型的玛多地震断层滑动分布(图4).结果显示,主断层的滑动以走滑为主,平均滑动角为-4.36°.地震引起的同震位错长约170 km,破裂主要集中在0~15 km深度范围内.最大滑动量为4.4 m,对应深度为6.97 km.分支断层最大滑动量2.9 m,平滑滑动角为-11.84°.由滑动分布导出的矩震量为1.61×1020N·m,相应的矩震级为MW7.4.
图4 GNSS与InSAR联合反演的玛多地震滑动分布(a)、(b)分别为主断层和分支断层上的滑动分布.Fig.4 Coseismic slip model of the Madoi earthquake inverted by GNSS and InSAR(a) and (b) represent the coseismic slip on the main fault and the secondary fault.
我们采用棋盘测试进行GNSS与InSAR联合反演的滑动分布分辨率检验(图5).保持断层模型的几何参数不变,仅对断层面上滑动矢量进行更改,但保持断层面上平均滑动量不变.结果表明GNSS与InSAR联合反演模型对15 km以上的滑动具有较好的分辨能力,在浅部能够支持的分辨尺寸为10 km×10 km.前述玛多地震的同震滑动位于0~15 km深度,对此区间内滑动分布特征有足够的分辨力.
GNSS水平同震位移的观测与模拟表明(图6),除4499测站外,震中附近GNSS测站的模拟结果均较好.远离震中区域的测站,断层南侧的测站拟合结果较好,断层北侧的测站拟合残差相对较大.从InSAR获取的同震形变模拟及残差分布来看(图7),模拟得到的形变场与观测值具有较好的一致性,ALOS2降轨影像鄂陵湖东南残差相对较大,可能与其包含了更多的震后余滑有关.整体残差,GNSS数据东方向为2.4 cm,北方向为2.4 cm,Sentinel-1升轨为4.0 cm,降轨为6.7 cm,ALOS2升轨为4.6 cm,降轨为5.5 cm.观测值与模拟值的相关系数为99.03%.
地震发生时,发震断层通常会产生相当大的静态同震滑移,并在近场和远场引起静态同震应力变化.计算和分析同震库仑应力变化是评价地震对周围断层影响的重要手段.根据库仑破裂准则,库仑应力变化定义为
Δδf=Δτ+μ′Δσn,
其中,Δτ代表沿断层面滑动方向的剪切应力变化;Δσn代表断层面法向正应力变化;μ′代表有效摩擦系数,其取值范围为0~1(Scholz,1990).库仑应力为正时表示加载,会促进断层的破裂,为负时表示卸载,会抑制断层的破裂.
基于反演的断层滑动模型,我们利用PSCMP/PSGRN软件(Wang et al.,2006)计算玛多地震引起的同震应力变化.有效摩擦系数设为0.4(King et al.,1994),剪切模量设为30 GPa,根据反演的滑动分布将深度设为7 km.结果表明,玛多地震使东昆仑断裂玛沁段应力明显上升,加载达0.05 MPa,超过了0.01 MPa的触发阈值,而东昆仑断裂托索湖段、玛多—甘德断裂、达日断裂及久治断裂上的应力均被卸载,后续发生地震的可能性降低.
图5 棋盘测试(a)、(b)分别为对主断层和分支断层的棋盘测试.Fig.5 Checkboard tests(a) and (b) represent the checkboard tests on the main fault and the secondary fault.
图6 GNSS同震位移观测值及模拟值Fig.6 The observed and simulated horizontal GNSS coseismic deformation
图7 InSAR同震形变场的观测、模拟及残差分布(a—c)、(d—f)、(g—i)、(j—l)分别为Sentinel-1升轨、Sentinel-1降轨、ALOS2升轨和ALOS2降轨的观测值、模拟值和残差.Fig.7 Observation, simulation, and residual for InSAR coseismic deformation(a—c), (d—f), (g—i), (j—l) represent the observation, simulation, and residual for Sentinel-1 ascend track, Sentinel-1 descend track, ALOS2 ascend track and ALOS2 descend track, respectively.
同震破裂模型显示主发震断层上的滑动以走滑为主,断层中部约100 km范围近乎纯走滑,东西两端则带有拉张分量.同震位错长达170 km,并上至地表,造成约160 km长的地表破裂.破裂主要集中在15 km以浅,最大滑动量为4.4 m,与由统计得出的最大地震滑动量与震级关系式(Wells and Coppersmith,1994)估算的最大滑动量一致.最大滑动量对应深度为6.97 km,与重定位后的震源深度(7.6 km)(王未来等,2021)十分接近.由滑动分布得出的矩震级为MW7.4,比USGS公布的结果略大,可能是因为同震形变场中包含了震后形变的信息.断层上的滑动呈不均匀分布,沿断层走向可识别出4个凹凸体,分别对应野外地质考察的昌麻河乡段、冬草阿龙湖段、黄河乡段和鄂陵湖南段(李智敏等,2021;潘家伟等, 2021).从东至西,4个凹凸体上最大滑动量分别为3.6 m、4.4 m、3.2 m和2.9 m.震源处的滑动约2 m,相对较弱,而后向东西两侧传播,且东侧的滑动强于西侧,表明玛多地震是一次不对称的双侧破裂事件,与其西侧的2001昆仑山口西8.1级地震类似(许力生和陈运泰,2004).
滑动分布模型表明玛多地震是一个左旋走滑型地震.InSAR同震形变场显示玛多地震造成了西起鄂陵湖,东至昌麻河,长约160 km的地表破裂,地表破裂带与昆仑山口—江错断裂重合.因此,我们认为,第四纪左行走滑的昆仑山口—江错断裂为2021年玛多MS7.4地震的发震断裂.该断裂与东昆仑断裂斜交,并平行于东昆仑断裂秀沟—托索湖段,应是东昆仑断裂南侧的分支断裂.
东昆仑断裂自托索湖往东,断裂走向急剧变化导致挤压隆起形成阿尼玛卿山.由于挤压弯曲的存在,阻止了断裂带上的滑动,因此会在其边缘形成次级断裂以释放积累的应变.震前卫星影像上,发震断裂的鄂陵湖段和昌麻河段地貌迹线不清晰,说明这2段断层较为年轻,正处于形成和贯通过程中(李智敏等,2021).潘家伟等(2021)认为该挤压弯曲阻止了玛多地震破裂沿发震断裂向东的继续传播,而InSAR同震形变场和反演的滑动分布表明,玛多—甘德断裂在此次同震破裂中已被贯穿,江错段和昌麻河段可能已实现贯通.江错段和昌麻河段若已完全贯通,该地区的应变分配可能会出现新的模式.
东昆仑断裂带玛沁—玛曲段为无历史地表破裂的地震空段(Van der Woerd et al.,2002;徐锡伟等,2002;Wen et al.,2007),一直备受国内外学者关注.通过探槽开挖,玛沁段早全新世以来古地震复发周期为500~1000 a,2000 a以来的地震复发周期则为600±100 a.玛曲段古地震复发周期则为1000~2000 a,2000 a以来的地震复发间隔为1000 a(李陈侠,2009).玛沁段最晚一次古地震事件距今已有514~534 a,玛曲段为1055~1524 a,地震离逝时间均已接近或超过其复发周期.且汶川地震对东昆仑断裂玛沁—玛曲段有应力加载作用(Toda et al.,2008;单斌等,2012).同震库仑应力结果(图8)显示玛多地震使东昆仑断裂玛沁段库仑应力加载达0.05 MPa,已经超过了触发阈值,该段发生强震的风险进一步升高.
图8 玛多地震的同震库伦应力变化Fig.8 Coseismic Coulomb stress changes of the Madoi earthquake
中国大陆几乎所有8级和80%~90%的7级以上强震都发生在活动地块边界带上(张培震等,2003),而此次玛多地震却发生在活动地块内部,这对活动地块理论的丰富和完善提出了新的要求.事实上,玛多地震并不是孤例,其南侧的达日断裂就曾于1947年发生M73/4地震,同震位错达2~4 m(梁明剑等,2020).巴颜喀拉块体内部发生7级以上强震或许与其东部的变形模式有关.东昆仑断裂自托索湖往东,走向急剧变化导致挤压隆起并在其南侧发育一系列与之近乎平行的断裂,沿东昆仑断裂的滑动速率也从10 mm·a-1降至3 ~ 5 mm·a-1.东昆仑南侧的次级断裂第四纪以来都可能有过强烈的活动并至今活跃,如达日断裂、玛多—甘德断裂和昆仑山口—江错断裂都有过历史地震活动(张裕明等,1996;中国地震局震害防御司,1999;熊仁伟等,2010).这些次级断裂同东昆仑断裂,甘孜—玉树—鲜水河断裂一起,共同调节着巴颜喀拉块体的向东挤出.同东昆仑断裂相比,这些次级断裂滑动速率相对较低,在以往的研究中对其关注有限.此次玛多地震的发生表明,尽管这些断裂滑动速率低,应变积累慢,但仍然可以孕育7级以上强震.同震破裂显示玛多—甘德断裂可能已在此次地震中被昆仑山口—江错断裂所贯穿,这意味着该地区的应变分配可能会出现新的模式,块体内部的次级断裂活动可能会增强,地震危险性也会进一步升高.
(1)利用GNSS和多源InSAR监测,获取了2021年5月22日玛多MS7.4地震的同震形变场.最大水平形变达1.2 m,最大LOS向形变约0.9 m.基于同震形变场反演的断层滑动模型显示,玛多地震为一次左旋走滑事件.同震位错达170 km,并上至地表,造成了约160 km的地表破裂.最大滑动量为4.4 m,对应深度6.97 km.玛多地震释放的矩震量为1.61×1020N·m,对应矩震级MW7.4.
(2)破裂模型表明沿主发震断层存在4个凹凸体,其中东部的两个凹凸体破裂规模较大,最大滑动量分别为4.4 m和3.6 m.西部的两个凹凸体规模相对较小,最大滑动量分别为3.2 m和2.9 m.玛多地震是一次不对称的双侧破裂事件.次断裂倾向与主发震断裂相反,倾向南,与余震精定位的结果相符.
(3)玛多地震的发震断层为昆仑山口—江错断裂.玛多—甘德断裂在此次地震中被昆仑山口—江错断裂所贯穿.同震库伦应力结果显示,玛多地震使东昆仑断裂玛沁—玛曲段库仑应力加载超过了触发阈值,该段发生强震的风险进一步升高.
致谢感谢审稿专家为本文提供宝贵的修改意见.感谢汪荣江教授提供的SDM及PSCMP/PSGRN软件.文中图件均采用GMT软件包绘制,在此一并感谢.