2. Department of Geographical and Earth Sciences, University of Glasgow, Glasgow G12 8QQ, United Kingdom
2. Department of Geographical and Earth Sciences, University of Glasgow, Glasgow G12 8QQ, United Kingdom
2008年10月6日,位于西藏拉萨市西北约73km的当雄县辖区内(29.8°N,90.3°E,www.csi.ac.cn)发生了MW6.3级中强地震.此次地震至少造成9人死亡,14人重伤,近2000户房屋破损,直接经济损失达4.1亿元人民币.该事件是继1993年尼木MW6.3地震以来(图 1,表 1),当雄地区活动最强烈、破坏最严重的地震事件.根据区域台网监测结果,震后短短一周内有千余次余震发生,沿北北西方向展布(www.csi.ac.cn)(图 1).地震的宏观震中位于亚东-谷露裂谷带北段的羊易地堑附近[1].该地区长期受EW向拉张作用,1~2.5 Ma以来在主边界断层东发育一系列次生正断层,在长期地质历史演化中逐渐形成多个断陷盆地,且盆中之盆、坎中之坎等地质现象也较为普遍[2].该地区地貌特征丰富,构造环境复杂,是我国地震活动最频繁的区域之一. 1411年当雄南MS8.0地震[3]和1952年当雄北MS7.5地震[4]就是该区内SN向断层活动的产物.
吴中海等[1]在震后迅速展开现场考察,所得地裂缝、地震滑坡、崩塌等现场破坏特征与地震学测定结果基本吻合,等震线分布与中国地震局公布烈度分布趋势较为一致,长轴方向为北北东-北东向.这与余震的空间分布略有不同.结合区域活动构造特征,吴中海等[1]判定此次地震的发震构造为羊易盆地西缘主边界断裂,断层倾角在66°~80°之间,属正断层运动.
我们通过InSAR资料的分析和同震位移场的反演则认定,2008年当雄MW6.3地震的发震构造位于羊易盆地的东缘断裂,是次级构造活动的产物.
2 InSAR同震信息及误差特征依托欧洲空间局(ESA)和中国科技部联合资助的“龙计划”国际合作项目,我们收集了覆盖此次地震13景ASAR图像.由于当雄地震震中地区地形起伏在4000~6000 m之间,高低错落,沟谷纵横,干涉像对的数据相干性较差.经比对筛选,资料20070415和20090419两景图像所构成像对的相干性最佳(表 2),因此,用以获取此次地震的同震位移场.
我们应用加州理工大学(Caltech)/JPL联合发布的开源SAR软件ROI_PACv3.01[5],采取2次差分方案对上述两幅图像进行差分干涉处理,以获取此次地震的同震位移场.由于干涉像对存在明显卫星轨道误差,我们利用2次最佳拟合曲面的方法[6]予以消弱.从图 2a可以看出,除西北方向部分地区未能完全解缠外,整幅干涉图像相位连续,条纹光滑清晰.形变区覆盖了整个羊易盆地,呈SN向25km、EW向20km的椭圆状,形变条纹以盆地东边缘为分界,呈东西两瓣.西瓣形变较大,LOS向最大位移量~0.3m;东瓣形变较小,只出现一个色周变化,LOS向最大位移量~0.06m.同震形变主要集中在图 2a虚线框部分,所以,这一区域的同震形变信息将被用于震源参数的反演.
由于DEM误差、卫星轨道误差和大气噪声影响[7],InSAR干涉相位中存有一定量的误差信息. DEM误差与像对的基线相关,垂直基线越长相位误差越大[8].文中所选像对垂直基线只有3.5m,即1m地形误差仅引入~9.8×10-4cm的形变误差.对地形精度优于10 m的SRTMv4.1数据[9]而言,该误差忽略不计.卫星轨道误差通常会造成长波长变化,远场曲面拟合的方法可以较好地解决该问题[6, 10, 11].干涉像对中大气扰动的影响往往不能被直接忽略[12].
大气扰动(尤其是大气水汽部分)的时空变化会直接引起模型参数的不确定性[13].在没有额外水汽数据情况下[14, 15],直接使用1D或2D统计模型[16, 17]来描述水汽噪声分布特征成为重要的技术手段.在本研究中我们采用这种技术以削弱大气扰动的影响.
首先,对干涉位移场中震中附近主要的形变区进行掩模处理,只保留非形变信息的远场噪声.如图 2b所示,噪声在-2~3.5cm范围内摆动.然后,根据(1)式
(1) |
结合2D快速傅里叶统计方法直接定量评价空间噪声的2D变异函数(Variogram)[18].式(1)中,(x,y)为空间任一点p,(x+hx,y+hy)为p平移后的p'位置.通过分析,我们得到如图 2c所示的|p-p'|2的统计信息.从图中可以清楚地看出,在沿横轴和纵轴方向0~80km范围内图像基本由暗红渐变成蓝色,全局最大的变异函数值为90 mm2,20km以外才出现较明显的色调抖动.0~20km的范围内,等值线近圆形分布,平均变异函数值为~25 mm2.由此看来,用于反演的资料范围内的大气扰动可以被视为各向同性.
2.3 大气噪声的解析描述在各向同性的噪声场中,任意点对之间的方差可以方便地使用与距离h(h>0)相关的理论模型加以描述[19].如上文所述,我们选用的干涉图像中的大气噪声分布呈各向同性特征,因此,可采用如式(2)描述的1D变协方差函数替代经验统计[20].
(2) |
式中,Ckl(h)表示干涉像对中第k和第l个点的噪声变差;h为两点的距离(km,h>0);σ2为标准方差;a和b为调整理论变差与统计变差间相似度的参数;J0为零阶贝塞尔函数.
将如图 2b所示的协方差统计数据作为Ckl(h)的观测值,通过使用混合颗粒群非线性算法[21]拟合观测值,使二者的相关性达0.99,由此确定的最优模型参数:σ2=68,a=19.6,b=1000.因此,我们要采用的大气噪声的解析描述或协方差函数被确定为
(3) |
根据同震位移场的干涉条纹的分布特征可以初步认定发震断层的滑动分布较为简单,因此,我们采用均匀采样法[22]降分辨率提取观测信息.在如图 2a所示的框区内最终获取了966个观测点值.由于发震断层没有明显的地表出露,断层位置(Lon,Lat)、断层埋深(Depth)、断层宽(Width)、断层长(Length)、走向(Strike)、倾角(Dip)等7个非线性参数和滑动矢量的2个分量参数同时成为反演求解的对象.为此,构建如下目标函数:
(4) |
其中,G为根据弹性半空间位错理论[23]计算的格林函数,d为观测数据,W为观测资料的权阵,由WTW=C-1直接确定[24],C为观测数据的协方差矩阵.在噪声可以忽略不计的情况下,可将相等的权重赋予所有观测值,因此,C通常可取为单位矩阵[25, 26].如果观测资料中存在不可忽视的非均匀扰动,而继续使用相等的权重,则反演结果会出现不可忽视的偏差.从下面噪声或不确定性实验可以看出,我们选用C=Ckl,反演结果会更好.同上,我们仍采用颗粒群优化算法[21, 27]求解式(4)描述的非线性系统.
我们仿照Parsons等[28]采用的蒙特卡罗方法,首先基于协方差矩阵模拟100组大气噪声,并添加给观测资料d,形成100组融有大气扰动的观测数据集dp(p={1,…,100}),然后由dp反演确定100组模型参数.同时,为了说明C=Ckl的加权方式更好,我们分别在单位阵(Model1)加权和协方差矩阵(Model2)加权两种模式下进行了上述实验.实验结果如图 3所示.可以清楚地看到,蓝色点集(Model2结果)更集中,且其中心与由原始InSAR观测直接反演的结果几乎完全重合(参见表 3);而灰色点集(Model1结果)比较分散;不过,Model1结果中大多数参数在95%置信水平上的标准方差(ξ)变化也不大,呈现出以中值(μ)为中心的高斯分布.还可以看出,滑动角与纬度、埋深与断层宽度、经度与宽度、经度与埋深之间出现一定的线性折中,而其余任意两维参数之间的关系几乎是随机的(如图 3).比较Model1和Model2两种模式下的反演结果与直接使用协方差函数加权(Optimal)情况下的反演结果(表 3),我们注意到,Model2的反演结果与Optimal的结果更接近,其中滑动角偏差最大,但仅为0.7°;而Model1的反演结果与Optimal的结果中偏差最大的为断层宽度,偏差为6km.由此看来,资料中的大气扰动会给反演结果带来一定影响,但是,只要对这种扰动具有足够的认识,并根据这种认识构建权重矩阵,这种大气扰动的影响可以最小化.
为了反演断层面上的非均匀滑动,以Optimal情况下(表 3)反演得到的均匀模型为基础,将断层面扩展到20km×18km,延伸到地表,并划分成1km×1km的离散网格,形成滑动分布断层模型(Distributed Slip Model,DSM).空间滑动与InSAR观测之间为线性关系.考虑先验光滑约束和权阵W,我们构建出如下线性方程:
(5) |
式中,Gs和Gd分别为DSM单位走滑和倾滑分量形成的地表三维扰动在卫星视线向投影的合成矩阵;d为LOS向InSAR观测数据;L为用于先验光滑约束的拉普拉斯方程[29];κ2为调整模型粗糙度的非负因子,本研究中κ2=2.2;m为待求的走滑和倾滑分量构成的矢量.
已有的研究表明,均匀断层模型的倾角不等于DSM的倾角[22, 30].因此,我们首先使用网格迭代法[22]修订断层倾角,然后反演空间滑动分布m.修订后的断层倾角为55°,反演得到的空间滑动分布如图 4d所示.DSM反演结果显示,该事件的滑动在空间上存在2个独立部分,主要滑动分布在沿倾向方向6~14km范围,最大滑动约2 m,平均滑动角为-124°;断面右上2~5km深度范围内还有最大滑动量为0.4m的滑动区,滑动以右旋为主.
为了检验结果的稳定性,我们应用前文所述dp数据集分别获取100组空间滑动模型,以讨论DSM的不确定性[31].由图 4e和图 4f可知,DSM在走向方向上的最大不确定性水平~0.2 m,位于沿倾向方向12~16km范围,平均不确定水平~0.08 m;倾向滑动的最大不确定性不足0.15m,且分布于沿倾向方向上10~16km范围内.可见,较大的不确定性只出现在较深的范围,而主要滑动区域(6~10km)的不确定性水平仅~0.1 m,仅为最大滑动量的5%.这表明DSM反演结果具有相当的稳定性与可靠性.同时还表明,资料对倾向滑动的约束优于走向滑动.这可能是地表南北形变量在SAR卫星LOS向投影有限所致(参见表 2中合成矢量).从图 4c可知,由最终确定的DSM参数模拟得到的InSAR结果与观测值之间的残差均值仅为0.3cm,最大残差仅为0.8cm.最大残差出现在断层与地表交界位置,此处残差很可能与地震造成的滑坡、塌陷[1]等非弹性形变相关.
4 2008年当雄MW6.3级发震构造藏南地区广泛分布的近南北向裂谷带(包括亚东-谷露裂谷带)是第三纪以来青藏高原区发生强烈东西向扩张的有力证据[2, 32],后期沿主边界正断层以东一系列快速断陷活动逐渐造就了现今当雄一带盆山交错的地貌特征[33].由ETM +资料(www.landcover.org)、经band7、4、2假彩色合成影像可知,在念青唐古拉山东麓出现多处的冲沟与水系扭曲样式非常一致,并存在多次折曲现象,由此可推断亚东-谷露断裂带当雄段在此曾经发生过左旋活动,而且后来的某一时期又发生过右旋运动.该构造是由于印度板块向欧亚板块俯冲,以致青藏高原地壳缩短,尔后夷平,并在第四纪以来急剧隆起等一系列事件作用的结果[34].这次当雄地震震中地区,该断裂带控制着自南向北分布的安岗地堑(Angang)、羊易地堑(Yangyi)和吉达果地堑(Geda),这些断陷盆地的西边界为主断层WZBF(图 5).在图像中还可以清晰识别出羊易盆地东缘分布着的梯状断层三角面和陡坎.水系折曲等线性信息也较突出.由此还可直接解译出两条近南北向展布的正断层(EBF-I和EBF-Ⅱ)(图 5).对上述次级断层的认识,与早期应用地貌学分析的方法[35]和野外探勘所得结果相一致[34].
我们将单一均匀断层模型(白框)和滑动分布模型(黑色网格)叠加于当地的地貌与构造背景(图 5).可以看出,断面在地表的投影覆盖了几乎整个羊易盆地,主要滑动区(白框范围)刚好位于盆地中心区内,断层的地表出露线与EBF-Ⅱ几乎重合,与余震的空间分布具有很好的一致性(图 1).考虑到利用InSAR资料直接确定的发震断层的走向和倾角的不确定性仅4°左右、定位精度在亚公里水平(参见图 3和表 3)、用InSAR资料确定的震中位置距WZBF断层近20km以及WZBF断层为东倾正断层[36]等信息,我们认为这次当雄地震的发震构造为EBF-Ⅱ断层,而不是吴中海等[1]判定的羊易盆地西缘主边界断裂WZBF.
5 结论2008年10月6日西藏当雄MW6.3地震的发震断层是位于羊易盆地东缘的次生断层,走向178°,倾角55°,倾向近正西,以正断层运动为主,兼具少量右旋走滑分量.尽管InSAR结果中仍然分布一定水平的误差,但应用蒙特卡罗方法分析结果表明,该误差对结果的影响有限,反演确定的发震断层的模型参数非常稳定.同时,InSAR资料反演模型与ETM+遥感影像解译结果相吻合,说明反演结果的可靠性.这次地震引起羊易盆地下沉约30cm,似乎是亚东-谷露裂谷带北段的羊易地堑下沉的表现.
基于各向同性变异理论对本文InSAR资料中误差统计分析可以较好描述误差分布特征.显然,由误差协方差矩阵约束下(Model2)的解的可信区间更小,且与等权重约束(Model1)结果存在一定偏差,充分显示了评价误差协方差矩阵的可靠性和必要性.数值实验表明,呈高斯分布的噪声会给反演结果带来一定的影响,但是,如果对噪声分布特征具有足够的认识,并根据此认识构建权重矩阵,这种影响几乎可以完全消除.
致谢感谢Marcotte教授关于快速傅里叶方法评价2D经验协方差图像的讨论;文中图 1,2,4(a~c),5皆由GMT4.5.1[37]绘制而成.感谢两位匿名评审专家的宝贵意见,使得行文更加严谨和完整.文中所使用ASAR资料均由中国科技部和欧洲空间局联合资助的Dragon-Ⅱ项目(ID:5305)提供.
[1] | 吴中海, 叶培盛, 吴珍汉. 2008年10月6日西藏当雄Ms6.6级强震的地震烈度控震构造和发震机理. 地质通报 , 2008, 28(6): 713–725. Wu Z H, Ye P S, Wu Z H. The seismic intensity, seismogenic tectonics and mechanism of the Ms6.6 Damxung earthquake happened on October 6, 2008 in southern Tibet, China. Geological Bulletin of China (in Chinese) , 2008, 28(6): 713-725. |
[2] | 吴珍汉, 胡道功, 刘崎胜, 等. 西藏当雄地区构造地貌及形成演化过程. 地球学报 , 2002, 23(5): 423–428. Wu Z H, Hu D G, Liu Q S, et al. The formation and evolution of tectonic landform of Damxung area in central Tibetan Plateau. Acta Geoscientia Sinica (in Chinese) , 2002, 23(5): 423-428. |
[3] | 吴章明, 曹忠权, 申屠炳明, 等. 1411年西藏当雄南8级地震发震构造. 中国地震 , 1992, 8(2): 46–52. Wu Z M, Cao Z Q, Shentu B M, et al. Seismogenic structure of the 1411 southern Damxung earthquake of M=8 in Tibet. Earthquake Research in China (in Chinese) , 1992, 8(2): 46-52. |
[4] | 吴章明, 曹忠权. 1952年西藏当雄7.5级地震烈度的再评定. 地球物理学报 , 1991, 34(1): 64–72. Wu Z M, Cao Z Q. Review on intensity of Damxung earthquake (Ms=7.5) in 1952, Tibet. Chinese J.Geophys. (in Chinese) , 1991, 34(1): 64-72. |
[5] | Rosen P A, Henley S, Peltzer G, et al. Updated repeat orbit interferometry package released. EOS, Trans.Amer.Geophys.Union , 2004, 85(5). DOI:10.1029/2004EO050004 |
[6] | Wright T J, Lu Z, Wicks C. Source model for the Mw6.7, 23 October 2002, Nenana Mountain earthquake (Alaska) from InSAR. Geophys.Res.Lett , 2003, 30(18). DOI:10.1029/2003GL018014 |
[7] | Li Z, Liu Y, Zhou X, et al. Using small baseline interferometric SAR to map nonlinear ground motion: a case study in northern Tibet. Journal of Applied Geodesy , 2009, 3: 163-170. DOI:10.1515/JAG.2009 |
[8] | 李振洪, 刘经南, 许才军. InSAR数据处理中的误差分析. 武汉大学学报(自然科学版) , 2004, 29(1): 72–76. Li Z H, Liu J N, Xu C J. Error analysis in InSAR processing. Geometries and Information Science of Wuchan University (in Chinese) , 2004, 29(1): 72-76. |
[9] | Farr T G, Caro E, Crippen R, et al. The shuttle radar topography mission. Reviews of Geophysics , 2007, 45. |
[10] | Biggs J, Wright T, Lu Z, et al. Multi-interferomgram method for measuring interseismic deformation: Denali Fault, Alaska. Geophys.J.Int. , 2007, 170: 1165-1179. DOI:10.1111/j.1365-246X.2007.03415.x |
[11] | Wright T J, Parsons B, England P C, et al. InSAR observations of low slip rates on the major faults of western Tibet. Science , 2004, 305: 236-239. DOI:10.1126/science.1096388 |
[12] | Goldstein R. Atmospheric limitations to repeat-track radar interferometry. Geophys.Res.Lett. , 1995, 22(18): 2517-2520. DOI:10.1029/95GL02475 |
[13] | Dawson J, Tregoning P. Uncertainty analysis of earthquake source parameters determined from InSAR: a simulation study. J.Geophys.Res. , 2007, 112: B09406. DOI:10.1029/2007JB005209 |
[14] | Li Z, Fielding E J, Cross P, et al. Interferometric synthetic aperture radar atmospheric correction: GPS topography dependent turbulence model. J.Geophys.Res. , 2006, 111: B02404. DOI:10.1029/2005JB003711 |
[15] | Li Z, Fielding E J, Cross P, et al. Interferometric synthetic aperture radar atmospheric correlation: Medium Resolution Imaging Spectrometer and Advanced Synthetic Aperture Radar integration. Geophys.Res.Lett. , 2006, 33: L06816. DOI:10.1029/2005GL025299 |
[16] | Hanssen R.Radar Interferometry: Data Interpretation and Error Analysis.Netherlands: Kluwer Acad., 2001 |
[17] | Lohman R B, Simons M. Some thoughts on the use of InSAR data to constrain models of surface deformation: noise structure and data downsampling. Geochem.Geophys.Geosyst. , 2005, 6: 12. DOI:10.1029/2004gc000841 |
[18] | Marcotte D. Fast variogram computation with FFT. Computers & Geosciences , 1996, 22(10): 1175-1186. |
[19] | Goovaerts P.Geostatistics for Natural Resources Evaluation.New York:Oxford Univ., 1997 |
[20] | Biggs J, Burgmann R, Freymueller J T, et al. The postseismic response to the 2002 M7.9 Denali Fault earthquake: constraints from InSAR 2003 2005. Geophys.J.Int , 2009, 176: 353-367. DOI:10.1111/j.1365-246X.2008.03932.x |
[21] | 冯万鹏, 李振洪. InSAR资料约束下震源参数的PS0混合算法反演策略.地球物理学进展, 已接受 Feng W P, Li Z H.A novel hybrid PSO/Simplex Algorithm for determining earthquake source parameters using InSAR observations.Progress in Geophysics (in Chinese), (Accepted) |
[22] | 冯万鹏, 李振洪, 李春来.利用InSAR确定2009年4月6日Mw6.3拉奎拉(Italy)地震最优震源模型.地球物理学进展, 已接受 Feng W P, Li Z H, Li C L.Optimal source parameters of the 6 Aprll 2009 Mw 6.3 L'Aquila, Italy earthquake from InSAR observations.Progress in Geophysics (in Chinese), (Accepted) |
[23] | Okada Y. Surface deformation due to shear and tensile faults in a half-space. Bull.Seism.Soc.Am. , 1985, 75(4): 1135-1154. |
[24] | Wright T J, Lu Z, Wicks C. Constraining the slip distribution and fault geometry of the Mw7.9, November 2002, Denali Fault earthquake with Interferometric Synthetic Aperture Radar and Global Positioning System data. Bull.Seism.Soc.Am , 2004, 94(6B): S175-S189. DOI:10.1785/0120040623 |
[25] | 孙建宝, 石耀霖, 沈正康, 等. 基于线弹性模型反演1997年西藏玛尼Mw7.5级地震的干涉雷达同震形变场一Ⅰ均匀滑动反演. 地球物理学报 , 2007, 50(4): 1097–1110. Sun J B, Shi Y L, Shen Z K, et al. Parameter inversion of the 1997 Mani earthquake from InSAR co-seismic deformation field based on linear elastic dislocation model-Ⅰ Uniform slip inversion. Chinese J.Geophys. (in Chinese) , 2007, 50(4): 1097-1110. |
[26] | 冯万鹏, 许力生, 许忠淮, 等. 利用InSAR资料反演2008年西藏改则Mw6.4和Mw5.9地震的断层参数. 地球物理学报 , 2009, 52(4): 983–993. Feng W P, Xu L S, Xu Z H, et al. Source parameters of the 2008 Gerze Mw6.4 and Mw5.9 earthquakes from InSAR measurements. Chinese J.Geophys (in Chinese) , 2009, 52(4): 983-993. DOI:10.3969/.issn.0001~5733.2009.04.015 |
[27] | Li Z, Feng W, Xu Z, et al. The 1998 Mw5.7 Zhangbei-Shangyi (China) earthquake revisited: a buried thrust fault revealed with interferomertic synthetic aperture radar. Geochem.Geophys.Geosyst , 2008, 9(4): Q04026. DOI:10.1029/2007GC001910 |
[28] | Parsons B, Wright T, Powe P, et al. The 1994 Sefidabeh (eastern Iran) earthquakes revisited: new evidence from satellite radar interferometry and carbonate dating about the growth of an active fold above a blind thrust fault. Geophys.J.Int. , 2006, 164: 202-217. DOI:10.1111/gji.2006.164.issue-1 |
[29] | Jónsson S, Zebker H, Segall P, et al. Fault slip distribution of the 1999 Mw7.1 Hector Mine, California Earthquake, estimated from satellite radar and GPS measurements. Bull.Seism.Soc.Am , 2002, 92(4): 1377-1389. DOI:10.1785/0120000922 |
[30] | Burgmann R, Ayhan M E, Fielding E J, et al. Deformation during the 12 November 1999 Duzce, Turkey, earthquake, from GPS and InSAR data. Bull.Seism.Soc.Am. , 2002, 92(1): 161-171. DOI:10.1785/0120000834 |
[31] | Funning G J, Parsons B, Wigth T J. Surface displacements and source parameters of the 2003 Bam (Iran) earthquake from Envisat advanced synthetic aperture radar imagery. J.Geophys.Res. , 2005, 110: B09406. DOI:10.1029/2004JB003338 |
[32] | Yin A, Kapp P A, Murphy M A, et al. Significant Late Neogene east-west extension in northern Tibet. Geology , 1999, 27(9): 787-790. DOI:10.1130/0091-7613(1999)027<0787:SLNEWE>2.3.CO;2 |
[33] | 吴中海, 赵希涛, 吴珍汉, 等. 西藏当雄一羊八井盆地的第四纪地质与断裂活动研究. 地质力学学报 , 2006, 12(3): 305–316. Wu Z H, Zhao X T, Wu Z H, et al. Quaternary geology and faulting in the Damxung-Yangbajing basin, southern Tibet. Journal of Geomechanics (in Chinese) , 2006, 12(3): 305-316. |
[34] | Armijo R, Tapponnier P, Mercier J L, et al. Quaternary extension in southern Tibet: field observations and tectonic implications. J.Geophys.Res. , 1986, 91(B14): 13803-13872. DOI:10.1029/JB091iB14p13803 |
[35] | Fielding E, Isacks B, Barazangi M, et al. How flat is Tibet?. Geology , 1994, 22(2): 163-167. DOI:10.1130/00917613(1994) |
[36] | 邓启东, 张培震, 冉永康, 等. 中国活动构造基本特征. 中国科学(D辑) , 2002, 32(12): 1020–1030. Deng Q D, Zhang P Z, Ran Y K, et al. Basic characteristics of active tectonics in China. Science in China (Series D-Earth Science) (in Chinese) , 2002, 32(12): 1020-1030. |
[37] | Wessel P, Smith W H F. New, improved version of generic mapping tools released.EOS, Trans. Amer.Geophys.Union , 1998, 79(47): 579. DOI:10.1029/98EO00426 |