The Extraction Method of Swath Topographic Profile in Spherical Coordinate SystemA Study from the Himalayan Arc orogenic Belt
-
摘要: 带状地形剖面是以高程为纵向参考,以距剖面起始端点的距离为横向参考进行制图表达。由于其可真实地反映地表形态,是新构造和活动构造地貌研究中的基本参考和研究对象。随着数字化基础地理信息成果的普及,很多GIS软件可实现地形剖面的自动化提取。但以往方法的位置信息均是建立在剖面投影平面的基础上,当研究区域范围过大,会引入非构造变形,为消除该影响,提出基于球面坐标系统的带状地形剖面图制作方法。该方法不仅适用于小型地质构造,也适用于喜马拉雅弧形造山带等大型地质构造,可提取出更接近真实地形的地形剖面。Abstract: The swath topographic profile is represented by elevation as the vertical axis and the distance from the beginning end of the profile as the horizontal axis. Because it can truly reflect geomorphological features, it is a basic reference and research object in the study of neotectonics and active tectonics. With the popularization of basic digital geographic information, many GIS software can realize the automatic extraction of topographic profile. However, the generation of topographic profile is almost based on projection.Therefore, no matter what kind of map projection is chosen, non-structural deformation will be introduced when the study area is too large. In order to eliminate this artificial distortion, a new method of making swath topographic profiles based on spherical coordinate system is proposed in this paper. It is not only suitable for small geological structures, but also for huge crust-scale structures, such as the Himalayan orogenic belt. The extracted topographic profiles could be much closer to the real terrain.
-
引言
地壳浅部的速度结构对震后地面的震动强度有很大影响,低速、松散沉积层对地震波具有强烈的放大作用,地震发生后,易造成建筑物倒塌,引发严重地震灾害。同时,浅层的速度结构影响了地震定位结果,其可靠性直接影响了震源位置及深度精度(谢军等,2012)。因此,明确地震高发区的浅层速度结构,对于震源研究、减轻地震灾害具有重要意义(沈伟森等,2010)。
目前,针对浅层地质体结构的研究较少,主动源探测及钻孔方法探测深度有限,且成本相对较高,因此使用受限(罗磊等,2019)。广大学者在对天然地震波形研究的过程中发现,地震面波具有频散特征,即面波的传播速度随频率的变化而改变,不同周期的面波具有不同的群速度,群速度的变化主要受一定深度范围内的地震波速度与密度影响,尤其对剪切波(S波)速度较敏感(刘丽等,2012),其变化对应不同深度的地下介质结构的横向差异,因此,可利用不同周期面波的群速度变化反映该周期对应深度的横向不均匀性,同时其采样深度一般随着周期而增加,因此可通过面波频散信息有效约束地下的速度结构(刘成林等,2018),该方法已取得丰硕成果(刘庆华等,2016;杨雅慧,2016;郑晨等,2016;郭希等,2017;冯策等,2018)。面波通常在远震中最为发育,但远震面波周期通常大于20 s,其对浅部地壳结构的约束不足,Shapiro等(2005)发现可从背景噪声中提取短周期面波信号,有效地提高了面波成像对地壳结构的分辨率,极大地推动了面波成像的发展(刘成林等,2018)。
为得到浅层的速度结构,李细兵等(2014)基于三角形台阵微动观测数据,提取场地小范围内的频散曲线,反演得到场地内数百米深度的纵、横波速度及深度信息等。张维等(2012)分别利用人工源与天然源面波探测浅部横波速度结构,结果表明,人工源探测深度较浅,锤击震源可得到数十米深度速度结构,而天然源探测深度可达几百米至上千米。罗艳等(2011)在四川遂宁-重庆潼南地震波形中,观察到部分台站短周期(1~2 s)瑞雷波特别发育。谢军等(2012)基于多重滤波法提取8个台站的面波频散曲线,并反演得到5个台站深度为10 km的S波速度结构。2015年,天津地区发生爆炸,此次事件具有强烈的面波信息,冀战波等(2018)基于同样原理提取面波频散曲线,反演得到天津地区浅层S波速度结构。
河北省唐山市位于玉田东南莫氏面隆起的斜坡上,断裂发育受地质构造等原因影响,一直为地震高发区。研究该地区浅层速度结构,对于地震精定位、活动性分析、灾害评估等均具有重要意义。2019年8月2日12时24分,河北省唐山市路南区某矿井区发生巷道塌陷(ML2.4),附近台站完整记录了此次塌陷的地震动波形,部分台站短周期瑞雷波发育,据此建立唐山及周边地区浅层速度结构,不仅对该地区地面振动及地震动灾害研究具有重要意义,且对于其他地区利用非天然地震事件短周期瑞雷波探测浅层速度结构具有借鉴意义。
1. 研究数据与模型
1.1 地震数据
唐山市某矿井区巷道塌陷形成的震中及部分台站波形如图1所示。震中附近CLI台站记录波形如图2所示,由上至下分别为径向(r)、切向(t)、垂直向(z)波形。根据面波周期与振幅均大于体波,传播速度小于体波,可以判断图2中红框内波形为面波,z向与r向分量中短周期瑞雷波较发育,且有明显的90°相移,周期较短,长周期面波较短周期面波先到,有较明显的频散现象。已有研究表明,面波同一周期群速度对1/3波长深度的剪切波速度结构最敏感(刘丽等,2012;谢军等,2012),探测深度约为半个波长(何正勤等,2007),因此本研究的短周期瑞雷波数据对浅层结构具有较好的约束。波形中t向分量面波发育一般,z向分量面波最为发育,因此仅利用z向分量短周期瑞雷波反演速度结构。
地震波形中一般含有多种频率成分,包括背景噪声、P波、S波、面波等。本研究需利用面波,提取其频散曲线,反演速度结构,因此在数据处理过程中,需准确提取面波信号,据此利用快速Fourier变换(FFT)对面波发育台站波形进行频谱分析,确定面波发育频率,用于后期滤波。LUX与CLI台站面波频谱分析结果如图3所示,由图3可知,本次地震面波主要的发育频率为0.1~0.4 Hz,因此,完成数据预处理后,选择0.1~0.4 Hz对所有台站进行滤波。LUX与CLI台站滤波前、后的对比如图4所示,由图4可知,滤波后体波被压制,面波得到了放大,有助于后续频散曲线的提取。
1.2 初始速度模型
初始速度模型对最终反演结果的影响较大,因此为得到较好的速度结构,通过查阅文献,最终选择冀战波等(2018)利用2015年天津爆炸事件得到的速度模型,并结合万永革等(2008)计算得到的首都圈地区速度模型及河北台网目前在用的华南速度模型,建立本研究使用的初始模型(图5)。
2. 处理过程
根据面波研究理论可知,不同频率的面波受地质结构的速度及密度影响,传播速度不同,且对S波更敏感,即可利用面波频散曲线反演得到S波速度变化。因此,本研究数据处理过程主要包括提取单台波形中短周期瑞雷波的群速度频散曲线,并在给定的初始模型基础上不断修改速度参数,使其反演频散曲线与提取的频散曲线之间的误差最小,即可反演得到震中至台站范围内S波速度结构。
2.1 提取频散曲线
目前提取面波频散曲线的方法主要有多重滤波法、移动窗分析法和连续小波变换法等,其中,多重滤波法在提取面波群速度频散曲线中的应用最广。多重滤波法又分为单台法与双台法,二者的主要区别在于震中是否位于研究区域内,单台法更适用于震中位于研究区域内的情形(蒋婵君等,2019),因此本研究利用单台法提取面波频散曲线,基本原理如下(周青云等,2006):
对原始地震图在频率域用中心频率为
$ {\omega _0} $ 的高斯无相移带通滤波器进行滤波:$$ {H_{\rm{n}}}\left( \omega \right) = \left\{ {\begin{array}{*{20}{l}} {0,}&{\omega < {\omega _{{\rm{ln}}}}}\\ {{\rm{exp}}\left[ { - \alpha {{\left( {\dfrac{{\omega - {\omega _{\rm{n}}}}}{{{\omega _{\rm{n}}}}}} \right)}^2}} \right],}&{{\omega _{{\rm{ln}}}}\leqslant\omega \leqslant{\omega _{{\rm{un}}}}}\\ {0,}&{\omega > {\omega _{{\rm{un}}}}} \end{array}} \right.$$ (1) 式中,
$ \alpha $ 为可调和的折衷因子(刘丽等,2012),用于平衡解在时间域、频率域与震中距之间的关系,取值范围一般为6.25~50,谢军等(2012)研究中表明震中距为0~200 km时,$ \alpha $ 取6.25,本研究震中距约为100 km,因此$ \alpha $ 取6.25;$ {\omega _{\text{n}}} $ 为中心频率;$ {\omega _{\ln }} $ 为中心频率$ {\omega _{\text{n}}} $ 的下限频率;$ {\omega _{{\text{un}}}} $ 为中心频率$ {\omega _{\text{n}}} $ 的上限频率。用以上滤波器滤除通带以外的所有分量,再用傅里叶反变换回到时间域,得到最大振幅对应的时间,认为其是以
$ {\omega _{\text{n}}} $ 为中心频率的群速度波包的到达时间。假设面波沿大圆路径传播,根据震中距与群速度波包的到达时间,可计算出以$ {\omega _{\text{n}}} $ 为中心频率的群速度(蒋婵君等,2019)。本次塌陷能量较小,震中距>120 km后,面波发育程度下降,所以选择震中距120 km以内的台站波形按照上述流程处理,并提取频散曲线。此范围内共有39个台站,其中宽频带台站12个,其余多为井下短周期台站,部分台站记录的面波不发育或淹没在背景噪声中,最终共有9个台站提取频散曲线(图6)。由于唐山市南部及天津地区多为井下短周期台站,所以本次提取频散曲线的台站主要分布在唐山市北部至秦皇岛地区(图7)。分析各台站提取的频散曲线可知,除距震中最近的DOH台站外,其余台站频散曲线一致性相对较好,均呈现出速度随周期增大而增大的趋势,即S波速度随着深度的增加而增大,符合地质规律。
2.2 反演速度结构
瑞雷波反演速度结构方法主要有遗传算法、最小二乘法、人工神经网络算法、蚁群算法及拟牛顿算法等(翟佳羽等,2010)。由于本研究反演一维速度结构,较简单,因此采用效率较高且结果较精确的带权重的多次线性迭代反演法(谢军等,2012),每次线性迭代反演均采用阻尼最小二乘法求解(吴腾飞,2018),当频散曲线拟合的标准误差小于可容忍的误差(0.002 km/s)时,迭代反演自动终止,此时获得最终的反演S波速度结构,如图8所示,图中蓝线表示初始速度模型,红线表示反演速度结构,黑色三角形表示提取频散曲线,绿线表示模拟频散曲线。由图8可知,模拟频散曲线与提取频散曲线拟合一致性较好,说明反演结果准确性较高。与初始模型相比,9个台站的反演结果均发生较大变化,地表S波速度由0.9 km/s增至2.4 km/s左右,5 km深度范围内,速度随着深度的增大变化明显,5~9 km深度范围内增大趋势变缓,部分台站附近出现低速层。
3. 结果分析
3.1 速度结构分析
图8中所有台站反演结果表明,表层S波速度约为2.46 km/s,与张学民等(2001a)和刘丽等(2012)的研究结果基本一致。闫静茹等(2017)对唐山响嘡台阵2号测井附近(DOH与LUX台站周边)开展钻孔测试时,发现该地区地下1.1~76 m岩层由土层变为花岗岩,S波速度由189.1 m/s增至2.17 km/s,本研究结果表明,DOH台站与LUX台站表层S波速度分别为2.25、2.24 km/s,结果较接近。
深度为2~3 km时,距震中位置最近的DOH、LUX台站速度最小,约为2.57 km/s,其他台站S波速度约为2.86 km/s,即震中附近存在1个低速区,其他台站构成1个相对高速区,这与王峻等(2009)的研究结果相符。
深度为4 km时,S波速度达3.47 km/s,与Crust 2.0结果基本一致。对比LUX、CLI、BDH、TLK、QIX台站反演结果,深度达4 km后,S波速度增长缓慢。其中BDH、TLK、QIX台站反演结果表明,当深度为5~9 km时,S波出现低速层,QIX台站反演结果尤为明显,速度最低,为3.02 km/s。相对于震中,KUC与QIX台站方位基本一致,但二者的反演结果差异较大,KUC台站反演结果中未表明深度达5 km后出现低速层,经分析认为,KUC台站提取的频散曲线周期较短,仅为6 s左右,根据面波频散对于S波的敏感程度可知,深度达6 km后的结果可信度下降,应以QIX台站反演结果为准。受传播路径的影响,JIX、ZUH、KUC台站提取的频散曲线周期最短,为5~6 s,此范围内S波速度随着深度的增加而增大。
根据唐山地区沉积特征可知,第四纪沉积层整体表现为东南部厚、西北部薄的特点,所以浅表S波速度也应表现为东南部速度略小于西北部。根据台站分布图可知,DOH、LUX、CLI、BDH、TLK台站位于东南部,浅表S波速度分别为2.25、2.24、2.54、2.69、2.57 km/s,平均速度为2.45 km/s;JIX、ZUH、QIX、KUC台站位于西北部,浅表S波速度分别为2.52、2.69、2.35、2.47 km/s,平均速度为2.50 km/s,西北部速度略高于东南部。同时对比DOH与QIX台站的反演结果可知,DOH台站浅表区域速度增加较慢,当深度为5 km时,S波速度仅达3.00 km/s,说明此范围内沉积层相对较厚,而QIX台站浅表S波速度相对较低,但在深度为1 km左右时,S波速度增大,由2.34 km/s增至2.87 km/s,结合地质情况可知,该台站位于沉积盆地的边界处,所以此区域范围内沉积层相对较薄,即受断裂沉积影响,沉积层厚度由DOH台站向QIX台站方向变薄。
根据唐山地区断裂带分布可知,DOH与LUX台站所处区域,南部为宁河-昌黎断裂,西部为陡河断裂,北部为丰台-野鸡坨断裂,东部为滦县-乐亭断裂(刘保金等,2011),该区域被4个断裂与周边割裂开来,形成小的凹陷,沉积层相对较厚。与其他台站相比,DOH台站与LUX台站周边整体处于低速区。根据张学民等(2001b)的研究结果,S波速度越低,上部脆性岩石所受剪切应力越大,越易产生破裂,形成易震层,通过查看河北历年地震目录可知,1976年唐山大地震位于此低速区内,相较于周边地区,DOH与LUX台站附近确为地震多发区(图9)。
3.2 精定位结果分析
拟合9个台站提取的频散曲线,并反演得到整个研究区域深度为10 km的速度结构。与单台结果相比(图10),5 km深度范围内结果一致性较好,波速随着深度的增加而增大,5~7 km深度范围内出现低速层,较单台反演结果范围小,应是各台站结果综合导致的。
为验证浅层速度结构在精定位中的作用,将本研究反演的模型结果与河北台网现用的华南速度模型相补充得到新模型(其中,P波波速利用经验公式 VP=1.73× VS计算得到),需注意的是,新模型仅丰富并细化了深度10 km以内的速度结构,其他并无变化。
选择2008年至今唐山地区ML1.5以上地震,利用双差定位方法,基于新模型与华南速度模型,进行精定位结果对比(图9)。相对于华南速度模型,新模型定位结果的震中位置更集中。2种模型深度定位结果均显示本地区震源深度由西南向东北变小的趋势,这与王晓山(2017)的研究结果一致。对比2种模型深度定位结果可知,未增加浅层结构的深度定位结果主要集中在5~15 km,深度为5 km以内的地震数量较少;增加浅层结构后,深度范围集中在0~5 km与5~25 km,深度集中性虽稍差,但深度<5 km的地震数量明显增加,由原来的15条增至220条。整体定位结果显示,增加浅层速度结构后,震源位置条带性更明显,主要集中于DOH与LUX台站下方。据前文可知,DOH台站与LUX台站在深度为10 km范围内为低速区,是易震区,与定位结果一致。另外,根据震源深度的条带性分布,DOH台站与LUX台站对应的低速区有近20 km的向下延伸。同时,新模型定位结果显示,在5 km深度处存在分界线,经分析认为,新模型是2个模型的简单叠加,仅丰富了浅层结构,并不能较好地约束深层地震定位结果,为进一步提高定位精度,可利用其他波形数据或技术手段丰富深层的速度结构。
综上所述,本研究所得模型不仅细化了当地的速度结构,且能较好地适用于精定位研究。
4. 结论与讨论
(1)本研究利用ML2.4塌陷产生的面波反演得到唐山北部地区S波速度结构,结果表明,该区域浅表S波速度约为2.46 km/s,深度为4 km时,S波速度约为3.47 km/s,塌陷周边沉积层相对较厚,沉积盆地边界处沉积层变薄。横向上,深度为2 km时,塌陷位置向东延伸至LUX台站范围内,存在1个低速区,S波速度约为2.57 km/s,且向下延伸近20 km,其他台站S波速度约为2.86 km/s,表现出较明显的横向变化。纵向上,整个研究区域内S波速度随着深度的增加而增大,但在唐山东部沉积盆地内5~9 km深度范围内存在S波低速层。
(2)9个台站中,JIX台站与ZUH台站震中距均为80 km,BDH台站与KUC台站震中距均为113 km,其中JIX台站与ZUH台站均位于盆地边缘,BDH台站位于沉积盆地内部,KUC台站位于沉积盆地外部。面波由震源传至台站过程中,JIX台站与ZUH台站路径基本相同,二者提取的频散曲线虽形态有所差异,但周期长度大致相同,均约为6 s;BDH台站与KUC台站路径相差较大,前者相对后者简单,提取的频散曲线周期长度相差较大,前者为8 s,后者仅为5 s(最短周期长度),说明瑞雷面波在相对简单的路径中传播衰减慢于复杂路径。
(3)本研究所得浅层速度结构在精定位中有较好的应用效果,说明唐山部分地区震源深度相对较浅,当震级较大时易造成较严重的灾害。本研究塌陷事件震级仅为ML2.4,能量较低,面波发育程度相对较低,且震中附近短周期台站较多,提取的频散曲线较少,周期多小于10 s,因此本研究初步得到小区域、浅层的速度结构,不能表征大区域、深层的地壳信息,后续可借助背景噪声或远震波形,反演大区域、深层的速度结构,为地震活动性的进一步研究提供数据。
致谢 感谢审稿专家提出的宝贵修改意见!
-
-
[1] 丁林, 钟大赉, 2013. 印度与欧亚板块碰撞以来东喜马拉雅构造结的演化. 地质科学, 48(2): 317—333Ding L. , Zhong D. L. , 2013. The tectonic evolution of the eastern Himalaya syntaxis since the collision of the Indian and Eurasian plates. Chinese Journal of Geology, 48(2): 317—333. (in Chinese) [2] 王春范, 2018. 基于ArcGIS环境的地形剖面图制作方法. 城市勘测, (6): 124—127 doi: 10.3969/j.issn.1672-8262.2018.06.031Wang C. F. , 2018. The method of making profile of topography based on ArcGIS environment. Urban Geotechnical Investigation & Surveying, (6): 124—127. (in Chinese) doi: 10.3969/j.issn.1672-8262.2018.06.031 [3] 张会平, 刘少峰, 2004. 利用DEM进行地形高程剖面分析的新方法. 地学前缘, 11(3): 226. doi: 10.3321/j.issn:1005-2321.2004.03.036 [4] Bendick R., Bilham R., 2001. How perfect is the Himalayan arc? Geology, 29(9): 791—794. [5] Bian S. , Gong J. F. , Zuza A. V. , et al. , 2020. Late Pliocene onset of the Cona rift, eastern Himalaya, confirms eastward propagation of extension in Himalayan-Tibetan orogen. Earth and Planetary Science Letters, 544: 116383. doi: 10.1016/j.jpgl.2020.116383 [6] Burbank D. W. , Anderson R. S. , 2013. Tectonic geomorphology, second edition. Environmental and Engineering Geoscience, 19(2): 198—200. doi: 10.2113/gseegeosci.19.2.198 [7] Ding L. , Zhong D. L. , Yin A. , et al. , 2001. Cenozoic structural and metamorphic evolution of the eastern Himalayan syntaxis (Namche Barwa). Earth and Planetary Science Letters, 192(3): 423—438. doi: 10.1016/S0012-821X(01)00463-0 [8] Fielding E., Isacks B., Barazangi M., et al., 1994. How flat is Tibet? Geology, 22(2): 163—167. [9] Grohmann C. H. , 2004. Morphometric analysis in geographic information systems: applications of free software GRASS and R. Computers & Geosciences, 30(9—10): 1055—1067. [10] Molnar P. , Lyon-Caent H. , 1989. Fault plane solutions of earthquakes and active tectonics of the Tibetan Plateau and its margins. Geophysical Journal International, 99(1): 123—153. doi: 10.1111/j.1365-246X.1989.tb02020.x -