通海台垂直摆受建筑载荷影响的数值模拟研究

岳 冲 李利波 闫 伟 洪 敏

1 中国地震台网中心，北京市三里河南横街5号，100045 2 中国地震局地质研究所地震动力学国家重点实验室，北京市华严里甲1号，100029 3 云南省地震局，昆明市北辰大道148号，650224

为精确测量地倾斜潮汐变化，目前常用的仪器有水管倾斜仪、水平摆倾斜仪和垂直摆倾斜仪，上述仪器的观测精度可以达到10-10～10-9。以往研究结果显示，地震之前往往存在趋势转折、破年变、短临加速等异常变化[1-3]。柯昌安等[4]研究发现，2008年宁强5.7级地震前7 d，安康、宝鸡、宁陕和西安台垂直摆、水管仪和伸缩仪均出现短临加速的异常现象；马震等[5]研究发现，2016年门源6.4级地震前，周边门源、湟源等形变台站倾斜类观测仪器出现短期趋势转折异常变化，各测项的异常时间较为同步，且各仪器出现异常的分量方位与震中方位较为吻合。地形变观测仪器能够直观地反映区域构造应力场的调整变化，进而捕捉地震前兆异常，但由于测量精度较高，观测仪器易受到各种干扰因素的影响，如气象因素、人为干扰、环境干扰、仪器更换或故障等[6-9]。山东泰安台应变及倾斜类观测仪器受8 km外水库抽水、补水等产生的载荷变化影响，出现2～3 d持续的短临加速变化[10]；2014年鲁甸6.5级地震前，距离震中不到50 km的昭通地震台出现多个测项的短临异常变化，研究发现该异常与降雨相关性较高[11]；Yue等[12]利用有限元模型分析发现，南北向长、东西向短的双阳湖水位增加更易在湖周边地表引起沿南北向形变，造成双阳台水管倾斜仪出现NS及EW方向差异性变形特征。因此，形变台站观测仪器出现变化后，除分析与地震相关的区域构造应力调整引起的变化外，也需要结合仪器不同测项的变化特征，分析其受到干扰因素的影响，进而对仪器变化性质做进一步判定。

通海台垂直摆NS和EW测项自2019-04起出现同步趋势转折变化，类似的趋势转折变化在2014年景谷6.6级地震前也出现过，但2019-04以来的趋势转折过程中存在建筑施工的影响。因此，本文针对通海台垂直摆NS与EW测项趋势转折变化，构建通海台建筑载荷影响三维有限元模型，定量分析建筑载荷对通海台垂直摆NS与EW测项影响的差异性特征，进一步探讨通海台垂直摆趋势转折变化机理，为进一步分析此次趋势转折变化对未来该区域应力调整、地震危险性评价影响等提供参考依据。

1 通海台垂直摆

1.1 通海台地质概况

通海台位于通海盆地南缘，距离杞麓湖南侧约3 km，通海盆地为断陷盆地，呈长条形EW走向。通海台始建于1971-07，位于宝秀-建水断裂北侧约1 km处[13]，地理坐标为102.78°E、24.10°N，海拔高程为1 821 m，观测山洞深38 m，全长54 m，观测硐室覆盖层厚23 m，硐室岩性主要为白云岩、石灰岩，山洞内部岩体部分破碎。

通海台地处川滇菱形块体南部，该区域受NW向石屏-建水断裂带、曲江断裂以及近NS-NNE向小江断裂带南段等控制和影响(图1)，上述区域历史地震活动强烈，通海、峨山、石屏等地发生过10多次6.0级以上破坏性地震，震级最大的一次为1970-01-05通海M7.8 地震。自1970年以来(截至2021-01-01)，该地区发生多次5级左右地震，尤其是2018-08-13与2018-08-14两次MS5.0地震均发生在通海盆地北缘，距离台站仅约12 km。

1.2 通海台垂直摆变化特征

目前通海台硐室内共布设有VS型和VP型垂直摆倾斜仪，仪器观测墩类型为石灰岩块墩。其中，VS型垂直摆倾斜仪自2008-11-08开始观测，截至2021-06，除NS测项于2010-01～02间存在部分缺数外，其他时段内数据完整性较好，可为本文研究提供较长的背景数据资料。VP型垂直摆倾斜仪于2019-08才开始观测，观测时间较短，因此数据缺乏可靠性。本文重点结合VS型垂直摆2019年前数据变化的背景趋势，对2019-04后NS和EW测项出现的同步变化特征进行分析。

图1 通海台位置Fig.1 Location of Tonghai station

VS型垂直摆2009年以来的观测数据资料如图2(a)、(b)所示，图中已对通海台周边历史地震进行标注(标注规则：震中距100 km内MS5.0～5.5地震，200 km内MS5.5～6.0地震，300 km内MS6.0以上地震)。由原始观测曲线可知，垂直摆NS测项2009～2014年一直持续N倾变化，2014年初N倾趋势出现转平后逐渐转为S倾变化，2019-04后再次转为N倾变化；垂直摆EW测项自观测以来一直持续E倾变化，但E倾趋势在2013年出现速率变缓，2019-04后E倾趋势再次加速。

图2 通海台垂直摆NS和EW测项原始曲线、拟合曲线、去趋势时序曲线Fig.2 Time series of original data, fitting data and de-trended data of vertical pendulum inclinometer NS and EW components of Tonghai station

单从垂直摆两个测项时序曲线来看，两个测线变化较为同步，均在2013年前后和2019-04出现趋势转折变化，而震例统计结果显示，2013年出现趋势转折后，距离台站259 km的景谷于2014-10-07发生MS6.6地震，2013年前后NS和EW测项的趋势转折变化是否与此次地震有关，还需进一步深入研究。自2019-04起，距离硐室北东侧约50 m处陆续开始建设医院等建筑物(图2中红色阴影时段)，垂直摆NS和EW测项数据变化时间与建筑物载荷影响时间较为吻合，但该测项变化是否均由建筑物施工造成以及NS和EW测项变化量为何相差较大等问题仍需作进一步分析研究。

进一步对垂直摆NS和EW测项作趋势拟合处理，以消除仪器零漂等产生的影响，其中NS测项采用幂函数进行拟合，相关系数为0.923；EW测项采用对数函数进行拟合，相关系数为0.97，结果如图2(c)、(d)所示。去趋势后的NS和EW测项在2019-04后趋势转折更加突出，且NS测项趋势变化量显著大于EW测项，2019-04后NS测项N倾变化量达到约2 750 ms，而EW测项E倾变化量达到约580 ms，NS测项约为EW测项的4.7倍。

为排除观测仪器更换、格值变化等对仪器观测的影响[9]，计算2017～2020年NS和EW测项M2波潮汐因子(图3)。由图可知，2017～2020年NS测项M2波潮汐因子在0.7左右波动，其中2018-08～09出现较为突出的M2波潮汐因子增大现象，可能是因为受到通海两次MS5.0以及墨江MS5.9地震的影响；EW测项M2波潮汐因子4 a间较为稳定，均在0.52左右波动。NS与EW测项在2019年的变化趋势(图中蓝色曲线)与其他年份相差不大，表明2019-04至今仪器观测稳定，未发生仪器标定等格值变化现象，而数据的转折变化未能在潮汐频段引起同步变化。

图3 通海台NS和EW测项M2波潮汐因子Fig.3 M2 wave tidal factor of vertical pendulum inclinometer NS and EW components of Tonghai station

2 模型设计

2.1 模型及边界条件设置

2019年观测硐室北东侧约50 m的县医院住院部开始施工，该工程2019-09开始搭建主体建筑，整个楼层共设计12层，总面积为12 643.1 m2。整个模型尺寸设计为600 m×600 m，为避免模型边界效应等产生的影响，将通海台设置于模型中央，距离模型左右边界分别为300 m，模型深度设计为800 m。医院位于通海台东北部，其西南角距通海台50 m，医院模型尺寸设计为50 m×50 m。对医院工程量载荷进一步简化，每层为0.5 m钢混结构(ρ=2 500 kg/m3)，共计12层，计算医院对地面压强为150 000 Pa(重力加速度取值10 m/s2)。地下介质结构参考距台站30 km的玉溪市防震减灾局地震监测井(DZKC1)，结合钻井及通海台山洞岩性特征，对模型沿深度方向进行剖分，其中0～100 m为砂岩，100～200 m为粘土层，200～800 m为石灰岩(图4(a))。整个模型底部沿垂直方向进行固定，模型表面自由，仅考虑建筑物产生的载荷压强作用，不考虑重力作用，采用弹性体模型进行模拟计算。整个垂直摆布设于石灰岩块墩上，折合摆长为10 cm，因此在实际模拟计算中将垂直摆简化为模型中心点位。为保证计算结果的准确性，对台站及医院区域网格进行加密，其中模型边界网格大小为50 m，台站及医院周边网格间隔为1 m，网格类型为C3D8R，共计获得466 760个网格单元(图4(b))。

图4 模型设置及网格剖分Fig.4 Model setting and grid segmentation

2.2 模型参数

为对比模型参数对模拟结果的影响，构建两类模型参数进行模拟计算：模型1中岩石杨氏模量参数取值为岩石参数下限，砂岩为2.5×107Pa，粘土层为1.0×108Pa，石灰岩为5.0×1010Pa[14]；模型2中砂岩和粘土层杨氏模量取值为模型1的4～5倍，具体参数设置见表1。

2.3 模拟结果

通过施加建筑物载荷作用，模拟计算两种模型地表位移结果(图5)。由图5(a)、(b)可知，沿NS向位移(N向位移为正，S向位移为负)结果显示，医院(图5中红色虚线框)南部形成N向位移，医院北部形成S向位移，通海台位于N向位移影响范围内；沿EW向位移(E向位移为正，W向位移为负)结果显示，医院西部形成E向位移，东部形成W向位移，通海台位于E向位移影响范围内。为进一步研究医院建筑载荷对通海台的影响，在模型表面作穿过通海台沿NS向剖面1和沿EW向剖面2(剖面1由N向S，剖面2由W向E)，剖面位移场结果如图5(c)、(d)所示。由图可知，剖面1三维形变场以垂向沉降为主，最大沉降位移在剖面线225 m距离处，最大沉降量达40.45 mm，通海台站处沉降量为-3.67 mm；剖面1沿EW向位移均为正值，表明该剖面均呈现E向运动特征，在225 m距离处E向位移量达到最大10.09 mm，通海台站处E向位移量为0.72 mm；剖面1沿NS向位移首先呈现S向位移运动特征，NS向位移在225 m距离处变为0后逐渐转为N向位移，通海台站处N向位移量为2.21 mm。由于剖面2未穿过医院建筑载荷区，该剖面三维形变场量级低于剖面1，三维形变场仍以垂向位移变化为主，在325 m距离处达到最大沉降量4.25 mm，EW向位移呈现E向运动转W向运动特征，而NS向位移则均为N向运动。

表1 模型参数设置及模拟结果

图5 模拟位移结果及剖面1和剖面2三维形变场Fig.5 Simulated displacement results and 3D deformation field of profile 1 and profile 2

由此可知，位于通海台东北部的医院建筑物载荷在台站附近产生N向以及E向位移加载，这与垂直摆NS和EW测项在2019-04出现的N倾和E倾趋势转折变化较为一致。模拟结果显示，在不考虑地形等因素影响的条件下，模型1和模型2的岩石参数均显示，建筑载荷对通海台产生的N向位移量大于E向位移量，模拟得到的N向位移量约为E向位移量的3倍。

3 讨 论

“十五”数字化改造后，由前兆数据观测原理和数据产出过程可知，绝大多数地壳倾斜、应变的初始产出物理量为电压，通过与仪器“格值”相乘后可得到相应的物理观测量，如应变量、倾斜量等[9,15]。由此可知，仪器观测结果并非地壳“真实”形变量，这会进一步增大本文模拟结果与实际观测结果直接对比的不确定性。通海台垂直摆摆体折合摆长为10 cm，垂直摆NS向自2019-04开始趋势转折，N倾转折变化量达到约2 750 ms。假设2019-04趋势转折前垂直摆仪器处于绝对垂直状态，N倾转折变化量转换为N向位移仅为0.001 3 mm，而该变形量远远小于模型1与模型2中NS向位移量。由于仪器硐室岩性主要为白云岩、石灰岩，而上述模型结果为建筑载荷对地表较软砂岩和粘土层的影响。为探究上述差异是否由模型参数设置所引起，本文对模型参数进一步调整，将模型中砂岩、粘土层参数均设置为石灰岩参数，模拟建筑物直接作用在石灰岩地层所产生的影响，得到通海台N向位移量为0.031 mm，E向位移量为0.004 4 mm，变化量仍大于仪器观测的结果。

建筑载荷作用到沉积层后在传递至仪器所在岩体的过程仍需更加深入的研究，但受到该地区精细的地下波速结果、钻孔信息等资料限制，本文对地下介质连续性、各向异性、物性参数、地表起伏变化等已进行必要的简化，且未考虑地下流体等因素的影响。通过上述简化模型模拟不同岩石参数下建筑载荷对通海台NS和EW测项的影响，结果均显示建筑载荷对通海台产生的形变场存在不均匀性，建筑载荷在NS向位移变化显著大于EW向位移变化，这与实际仪器中NS测项变化量大于EW测项变化量的结论一致。

4 结 语

针对通海台垂直摆NS和EW测项自2019-04出现同步趋势转折变化，利用三维有限元模型模拟建筑载荷对通海台NS向和EW向位移影响的差异特征，主要结论如下：

1)对垂直摆NS、EW测项作趋势拟合处理，以消除仪器零漂等产生的影响，结果表明，2019-04后NS测项N倾变化量达到约2 750 ms，而EW测项E倾变化量达到约580 ms，NS测项趋势变化量显著大于EW测项，NS测项约为EW测项的4.7倍。

2)计算2017～2020年NS和EW测项M2波潮汐因子，NS测项M2波潮汐因子在0.7左右波动，EW测项M2波潮汐因子4 a间较为稳定，均在0.52左右波动；2019-04至今仪器观测稳定，未发生仪器标定等格值变化现象，而数据的转折变化未能在潮汐频段引起同步变化。

3)位于通海台东北部的医院建筑载荷在台站附近产生N向及E向位移加载，与垂直摆NS和EW测项在2019-04后出现的N倾和E倾趋势转折变化较为一致。在不考虑地形等因素的影响下，模型1和模型2结果均显示，建筑载荷对通海台产生的N向位移量均大于E向位移量，N向位移量约为E向位移量的3倍。

致谢：感谢中国电建集团北京勘测设计研究院有限公司李善飞高级工程师在软件使用过程中提供的帮助。