2023, Vol. 43
2. 重庆市地震局，重庆市红黄路339号，401147;
3. 青海省地震局，西宁市兴海路1号，810001
早在20世纪80年代，地震相关分析和地质建模研究中就引入了数值模拟方法[1]，特别是有限元数值模拟方法，可用于三维非线性多场耦合的复杂建模分析，目前已成为主流建模分析方法之一[2-7]。本文利用有限元分析方法，基于地质构造资料构建五指山台台址几何模型，结合GPS和重力观测结果，分析五指山台定点形变观测数据趋势变化的物理意义，为该区域地壳应力与应变状态的长期变化研究提供资料。
1 五指山台概况五指山台位于五指山市冲山镇阿里山度假村旁太平山半山腰(109.53°E，18.79°N)，台站正南方25 km处发育有尖峰-吊罗深大断裂，距台站10 km处则有一条隐伏平移断层，该隐伏断裂极有可能是尖峰-吊罗深大断裂的延长发育断裂，对台站观测有一定影响。研究区域广泛出露二叠纪和三叠系花岗岩[8]。台站洞体基岩坚硬完整、致密均匀，测点周边观测环境较为稳定，为台站提供了较好的观测背景场。台站为无人值守台站，“十五”期间架设，观测山洞由部队防空洞改建。山洞内观测仪器有DSQ型水管倾斜仪、VP型宽频带垂直摆倾斜仪和SS-Y型铟瓦棒伸缩仪，其中DSQ型水管倾斜仪和SS-Y型铟瓦棒伸缩仪布设在主洞室内，两者相同分量的仪器部件平行布设且共用同一槽体，分别监测该区域的应变和倾斜变化；VP型宽频带倾斜仪布设在主洞室南侧的侧洞室。
2 五指山台地质环境五指山地区在大地构造上属于西太平洋地壳构造不同发展阶段的大陆边缘区，成土母岩为花岗岩和流纹岩，地貌主要为陡坡火山地貌和侵入岩体构造地貌。区内出露地层有长城系抱板群峨文岭组，是海南岛结晶基底的主要组成部分，侵入岩主要形成于二叠纪、三叠纪、侏罗纪和白垩纪。研究区位于昌江-琼海构造带以南、尖峰-吊罗断裂带以北，主要发育NW向构造[9]。
2.1 地层特征以横贯海南岛的EW向王五-文教大断裂及九所-陵水大断裂为界，2条断裂之间为五指山地层分区。区内发育前震旦系、震旦系、古生界地层，前震旦系地层以具复理石韵律结构的砂泥质岩及火山岩为主，少量碳酸盐岩；震旦系、下古生界地层以具复理石韵律结构的粉砂泥质岩石为主，少量砂岩、碳酸盐岩及酸性、基性、超基性火山岩，具地槽沉积特点。石炭纪-早三叠世为相对稳定的准地台发展阶段[10]。
观测洞室的台基为角闪黑云二长花岗岩，厚度较大，约为40 m，下部为薄层岩层碳质粉砂岩，有时与呈条带状、条带纹白云岩互层，其下与碳质板岩呈过渡状态，或夹碳质板岩，厚约50 m，其下与石碌群上亚群为假整合接触。
2.2 岩性特征由于不同岩性的物质成分和结构不相同，所以其物理特性和岩石力学参数具有很大的差异。海南岛岩浆岩分布广泛，具有多期次活动特征，其中侵入岩以海西-印支期斑状花岗岩、燕山期花岗岩及花岗闪长岩为主[11]。海南岛内花岗岩出露非常广泛，乐东-五指山-万宁地区出露面积约800 km2中二叠世同碰撞强变形花岗岩。结合五指山区域的地层分布特征，该区域出露的岩性主要为二叠纪花岗岩，以黑云母二长花岗岩和角闪黑云二长花岗岩为主。
2.3 构造应力场特征长期以来，由于植被覆盖严重，对海南岛的地质研究较少。统计雷琼地区MS≥4.5地震和局域小震的震源机制解可知，琼北至雷州半岛地区平均主压应力方向为NNW向，琼南地区平均主压应力方向为近NS向；断层面走向主要为NW向和NE向，断层面倾角大。全区主压应力轴(与海岸线正交)具有明显的扇形分布特征，反映了东南沿海构造应力场的基本状况，即来源于太平洋板块向西的挤压和印度洋板块向北和东的推挤[12]。
3 GPS及重力场分析 3.1 GPS应力场从单个GPS站点数据(图 1)来看，海南三亚(HISY)站2015-11以来呈SW向偏移，海南琼中(QION)站2010-06以来呈NW向偏移，海南海口(HIHK)站则一直较为平稳。从GPS站点间基线数据(图 1)来看，HIHK-QION虽表现为缓慢缩短挤压状态，但HIHK-HISY、QION-HISY均表现为快速拉伸状态，标志着整个海南岛主要处于地壳拉伸状态。地壳拉伸可能会引起地壳下降(沉)，而不均匀下降可引起地壳S倾，从而使山体出现S倾。
|
图 1 GPS基线变化曲线 Fig. 1 GPS baseline change curve |
从海南重力场变化(图 2)来看，2015-09~2016-09海南呈现东北部正变化、南部和西部负变化的过程，零值线穿琼东南而过，形成重力场正负异常变化高梯度区，五指山台周边地区重力变化约为－10 μGal。2016-09~2017-10海南重力场变化比较平稳，基本在30 μGal以内，南部为正变化，北部为负变化，五指山台周边地区重力变化约为10 μGal。故2015-09以来，五指山台周边地区的重力场变化经历了先负后正的变化过程。
|
图 2 海南重力场动态变化 Fig. 2 Dynamic variation of gravity field in Hainan |
研究结果表明，重力值减小说明地表隆升或地下物质减少或两者叠加；重力值增大说明地表沉降或地下物质增多或两者叠加。比较琼中台重力非潮汐变化(图 3)和GPS站基线变化(图 1)可知，2017-01以来琼中台重力非潮汐变化处于小幅上升变化趋势，琼中GPS基线垂直分量自2017年中期开始处于缓慢下降过程，两者的变化均表明2017年初以来琼中台地区地表小幅沉降，而2017-01以来GPS基线琼中-三亚呈伸长趋势(图 2)。结合这些变化特征认为，2017年以来海南南部地区均处于地表沉降的过程，且北部的沉降幅度小于南部。
|
图 3 琼中台重力非潮汐变化和五指山台水管倾斜仪NS分量变化 Fig. 3 Non tidal variation of gravity at Qiongzhong station and NS component variation of water pipe tiltmeter at Wuzhishan station |
选用ANSYS Workblench 16.0软件进行建模分析。
4.1 问题描述与条件假设自2015年观测以来，五指山台3套仪器工作正常，积累了丰富的观测资料。经对比分析发现：1)在基线布设完全一致的情况下，倾斜类仪器的映震能力远高于应变类仪器；2)2套倾斜仪能记录到明显的固体潮，而应变仪记录的固体潮幅度较小，特别是EW分量几乎记录不到同震响应；3)2套倾斜仪的长期趋势不一致，垂直摆趋于N倾，而水管仪出现S倾现象。
为便于建模分析，进行以下假设：1)用简单的纵向拉伸地质剖面图模拟实际的山体模型；2)由于模型是针对山体整体的特性进行分析，故对岩性和构造要素作简单的归一化处理；3)由于大部分研究结果都是关于应力方向和GPS速率的，为便于显示，假设应力大小与杨氏模量量级相同。
4.2 模型建立与分析 4.2.1 几何模型建立为便于模型分析、提高运算效率，选用应力45°方向AB剖面(图 4)进行建模，并结合实际地质资料，参考岩石力学参数表设定模型的材料属性(表 1)。
|
图 4 分析模型 Fig. 4 Analysis model |
|
表 1 五指山台地壳介质岩石力学参数汇总 Tab. 1 Summary of rock mechanical parameters of crustal media at Wuzhishan station |
1) 网格划分。选用Workblench中的自动网格划分工具，为保证模型网格密度和计算质量，综合考虑计算精度和计算量，选择网格尺寸为1.5 m，其他参数均采用默认值，满足模型需求共生成95 400个单元、409 128个节点。
2) 接触关系。Workblench提供多种接触类型，但由于不同岩层间的摩擦系数和粘合系数很难确定，本文采用简单的绑定接触以达到模型需求。
3) 边界条件。载荷和约束是ANSYS软件求解的边界条件，它们是以所选单元的自由度形式定义的。载荷主要选择力载荷工具的压力载荷，使用固定底部、两边拉张及两边压缩模型计算的位移和应力结果一致。本文在此只根据区域构造应力方向计算总体变形结果和等效应力分布结果，大小设定为1 010 Pa。
4.3 模型结果解释本模型主要研究AB剖面侧面45°受压时台址山体的变形特征和等效应力分布特征。但由于区域应力大小的不确定性，本文同时模拟了两边压缩和两边拉张工况(图 5和表 2)。可以发现，两边压缩和两边拉张工况下，除了变形大小和应力影响范围变大外，其他结果与侧面45°受压时趋于一致。等效应力分布结果显示，最大应力为78.288 MPa，两边受压或者两边受拉最大总体变形结果一致，受压应力45°方向后，在剖面A岩层上部应力的方向发生明显变化。
|
图 5 两边压缩、两边拉张、侧面45°受压情况下有限元模型分析结果 Fig. 5 Finite element model analysis results under compression on both sides, tension on both sides and 45° compression on one side |
|
|
表 2 对剖面等效应力的模型结果 Tab. 2 Model results of equivalent stress of profile |
总体变形结果显示，最大总体变形为19.116 mm，水平最大变形为17.881 mm，垂向最大变形为6.803 mm。黑云母二长花岗岩和断层碎裂岩更容易发生变形，而角闪黑云二长花岗岩变形幅度较小，这可能与研究区的岩石占比有关。受岩层之间摩擦力和其他接触力的影响，模型上部的变形幅度明显大于下部，当下部持续缩短，便容易形成上拱的背斜，这与实际的地貌相符合。由于岩层的倾向趋于SE，造成变形趋势也趋于SE。
4.4 实测数据验证结合实际观测资料(图 6)发现，2015-01-01~2019-12-01水管倾斜仪数据曲线整体呈SE倾(NS分量趋势下降，EW分量趋势上升)，反映观测岩层呈SE倾，与数值模拟结果相似；垂直摆倾斜仪数据曲线长期趋势不明显(自观测以来，NS分量趋势上升，EW分量数据趋势无规律)，部分时间段出现较大的波动，之后趋于恢复，其物理意义需进一步讨论；伸缩仪NS分量呈长时间张性(上升)，EW分量2018-08-06急速压性变化(下降)，2018-11-26趋势转折缓慢张性(上升)，整体呈张性趋势。
|
图 6 五指山台形变仪器观测数据曲线 Fig. 6 Observation data curve of deformation instrument at Wuzhishan station |
在建台初期，水管倾斜仪与伸缩仪位于同一仪器墩、同一水平面上，由于构造地质作用岩层呈小角度S倾。当受到构造应力挤压及南部地表沉降作用后，地层发生塑性缩短变形并挤压上层岩石，造成上覆岩层上拱形成一定弧度，使对应基线拉长，这是伸缩仪两分量观测结果为张性变化的原因。而水管倾斜仪利用连通器原理观测基线两端水位变化，由于基线较长，两侧的高差容易发生变化，因此2017年起其NS分量呈缓速S倾长期变化趋势。垂直摆倾斜仪通过摆体的倾角记录基点岩层的倾斜变化，NS分量长期趋势为N倾，与仪器受零漂影响有关；EW分量趋势不规律，部分时段出现较大波动，之后趋于恢复，其物理意义需进一步讨论。
5 结语1) 五指山地区构造较复杂，因地质资料的稀缺，仅以出露地质图作为山体模型，基底岩性差异未知，且长期受NNW向挤压应力影响，在一定程度上影响着区域应力场的大小和方向。
2) 有限元模型分析结果显示，在NNW向压力作用下，模型上部的变形幅度明显大于下部，下部持续缩短便容易形成上拱，这与GPS及重力场结果相符。由于岩层的倾向趋于SE，造成变形的趋势也趋于SE倾。
3) 结合五指山台各仪器观测数据进行分析可知，水管倾斜仪利用连通器原理观测基线两端水位变化，由于基线较长，基线两侧的高差容易发生变化，这是2017年起水管仪NS分量出现缓速S倾长期变化趋势的原因。数据曲线整体呈SE倾(NS分量趋势下降，EW分量趋势上升)，反映观测岩层呈SE倾，这与数值模拟结果相似。垂直摆通过摆体的倾角记录基点岩层的倾斜变化，其长期趋势不明显(自观测以来NS分量趋势上升，EW分量无规律)，部分时间段出现较大波动，之后趋于恢复，其物理意义需进一步讨论。伸缩仪NS分量呈长时间张性(上升)，EW分量2018-08-06急速压性变化(下降)，2018-11-26趋势转折呈缓慢张性(上升)趋势，整体呈张性。当受构造应力挤压及南部地表沉降作用后，地壳发生塑性缩短变形并挤压上层岩石，使上覆岩层上拱形成一定弧度，对应基线拉长，很好地解释了伸缩仪两分量观测结果整体为张性变化的原因。
| [1] |
England P, McKenzie D. A Thin Viscous Sheet Model for Continental Deformation[J]. Geophysical Journal International, 1982, 70(2): 295-321 DOI:10.1111/j.1365-246X.1982.tb04969.x
(  0) |
| [2] |
Parsons T. Post-1906 Stress Recovery of the San Andreas Fault System Calculated from Three-Dimensional Finite Element Analysis[J]. Journal of Geophysical Research, 2002, 107(B8): 2 162 DOI:10.1029/2001JB001051
( 0) |
| [3] |
邓志辉, 宋键, 孙君秀, 等. 数值模拟方法在地震预测研究中应用的初步探讨(Ⅰ)[J]. 地震地质, 2011, 33(3): 660-669 (Deng Zhihui, Song Jian, Sun Junxiu, et al. Preliminary Study on Application of Numerical Simulation Methods to Earthquake Prediction Research(Ⅰ)[J]. Seismology and Geology, 2011, 33(3): 660-669 DOI:10.3969/j.issn.0253-4967.2011.03.015)
( 0) |
| [4] |
申通, 王运生, 吴龙科. 重庆小南海滑坡形成机制离散元模拟分析[J]. 岩土力学, 2014, 35(增2): 667-675 (Shen Tong, Wang Yunsheng, Wu Longke. Discrete Element Simulation Analysis of Formation Mechanism of Xiaonanhai Landslide in Chongqing City[J]. Rock and Soil Mechanics, 2014, 35(S2): 667-675)
( 0) |
| [5] |
叶正仁, 王建. 中国大陆现今地壳运动的动力学机制[J]. 地球物理学报, 2004, 47(3): 456-461 (Ye Zhengren, Wang Jian. Dynamics of Present-Day Crustal Movement in the China Mainland[J]. Chinese Journal of Geophysics, 2004, 47(3): 456-461 DOI:10.3321/j.issn:0001-5733.2004.03.014)
( 0) |
| [6] |
龚丽文, 邓志辉, 陈丽娟, 等. 基于台址构造环境的有限元建模分析: 以黔江台为例[J]. 地震学报, 2019, 41(1): 80-91 (Gong Liwen, Deng Zhihui, Chen Lijuan, et al. Analyses of Finite Element Model Based on Station's Tectonic Environment: Taking Qianjiang Station for Example[J]. Acta Seismologica Sinica, 2019, 41(1): 80-91)
( 0) |
| [7] |
曹建玲, 张晶, 王辉, 等. 首都圈现今断层活动方式的数值模拟[J]. 地震, 2013, 33(3): 116-123 (Cao Jianling, Zhang Jing, Wang Hui, et al. Numrical Simulation of Fault Deformation in the Capital Region of China[J]. Earthquake, 2013, 33(3): 116-123)
( 0) |
| [8] |
王超, 魏昌欣, 云平, 等. 海南岛五指山地区顺作花岗岩锆石U-Pb年龄、地球化学特征及其地质意义[J]. 地质通报, 2019, 38(8): 1 352-1 361 (Wang Chao, Wei Changxin, Yun Ping, et al. Zircon U-Pb Age, Geochemistry and Geological Significance of Shunzuo Granite in Wuzhishan Area, Hainan Island[J]. Geological Bulletin of China, 2019, 38(8): 1 352-1 361)
( 0) |
| [9] |
陈哲培. 海南省岩石地层[M]. 武汉: 中国地质大学出版社, 1997 (Chen Zhepei. Lithostratigraphy of Hainan Province[M]. Wuhan: China University of Geosciences Press, 1997)
( 0) |
| [10] |
许德如, 范蔚茗, 梁新权, 等. 海南岛元古宙变质基底性质和地壳增生的Nd、Pb同位素制约[J]. 高校地质学报, 2001, 7(2): 146-157 (Xu Deru, Fan Weiming, Liang Xinquan, et al. Characteristics of Proterozoic Metamorphic Basement in Hainan Island and Its Implications for Crustal Growth: Nd and Pb Isotope Constraints[J]. Geological Journal of China Universities, 2001, 7(2): 146-157)
( 0) |
| [11] |
陈新跃, 王岳军, 张玉芝, 等. 海南晨星安山质火山岩地球化学、年代学特征及其构造意义[J]. 大地构造与成矿学, 2013, 37(1): 99-108 (Chen Xinyue, Wang Yuejun, Zhang Yuzhi, et al. Geochemical and Geochronological Characteristics and Its Tectonic Significance of Andesitic Volcanic Rocks in Chenxing Area, Hainan[J]. Geotectonica et Metallogenia, 2013, 37(1): 99-108)
( 0) |
| [12] |
李辉, 申重阳, 孙少安, 等. 中国大陆近期重力场动态变化图像[J]. 大地测量与地球动力学, 2009, 29(3): 1-10 (Li Hui, Shen Chongyang, Sun Shaoan, et al. Dynamic Gravity Change in Recent Years in China Continent[J]. Journal of Geodesy and Geodynamics, 2009, 29(3): 1-10)
( 0) |
2. Chongqing Earthquake Agency, 339 Honghuang Road, Chongqing 401147, China;
3. Qinghai Earthquake Agency, 1 Xinghai Road, Xining 810001, China