目前,国内外众多学者根据不同地震动特性的地震动参数与地震烈度的关系确定了相应的仪器地震烈度算法(关曙渊等,2020; 金星等,2013; Wu et al,2003; 崔建文等,2006)。《中国地震烈度表》(GB/T 17442—2020)给出了我国现行仪器地震烈度的计算流程及其定量化测定方法,在对数据进行基线校正、记录转换、数字滤波后,三分向合成地震动参数计算公式如下:

根据式(1)(2)分别计算I_A、I_V值:

I_A=3.17lg(PGA)+6.59(3)

I_V=3.00lg(PGV)+9.77(4)

最后根据I_A和I_V值计算仪器测定的地震烈度I_I为:

式中:PGA为合成地震动加速度记录的最大值,单位为m/s²; PGV为合成地震动加速度记录的最大值,单位为m/s; a(t_i)为t_i时刻点合成的三分向加速度值; v(t_i)为t_i时刻点合成的三分向速度值; 下标E-W、N-S、U-D分别表示东西方向、南北方向、垂直方向。

空间插值是一种通过离散空间数据计算未知空间数据的方法,一般分为确定性方法和地质统计学方法两类(李海涛,邵泽东,2019)。

确定性插值方法是基于已知样本点之间的相似程度或者整个曲面的光滑性来创建一个拟合曲面。主要有反距离权重法趋势面法、样条函数法等; 其中反距离权重法作为一种精确插值方法,算法简便易行,当样本数据分布均匀、密度较高足以反映局部差异时具有较好的插值精度(朱求安等,2004; 张海平等,2017)。

地质统计学插值方法是利用样本点的统计规律,使样本点之间的空间自相关性定量化,从而在待预测的点周围构建样本点的空间结构模型,比如克立金插值法。克里金插值法作为地质统计学的主要内容是一种以变异函数理论和结构分析为基础,在有限区域内对区域化变量进行无偏最优估计的一种方法,一般分为普通克里金、简单克里金、协克里金、泛克里金、指示克里金等(牟乃夏等,2012; 郑向向,帅向华,2013)。

克里金插值法与反距离权重法类似,两者均通过已知样本点权重来求得未知样本点的值,可统一表示为:

式中:Z(x₀)为未知样本点的值; Z(x_i)为未知样本点周围的已知样本点的值; λ_i为第i个已知样本点对未知样本点的权重; n为已知样本点的个数。

相对于反距离权重法只考虑已知样本点与未知样本点的距离远近,克里金插值法还通过变异函数和结构分析,考虑了已知样本点的空间分布及与未知样本点的空间方位关系,能较好反映客观地质规律。其次,地质统计学方法的插值精度与选择的变异函数模型密切相关,常用的变异函数理论模型有圆形模型、球形模型、指数模型、高斯模型等(汤国安,杨昕,2012)。

本文基于以上两类插值方法基本原理,考虑到此次漾濞地震具有较密集的预警台网数据,为寻找最优空间插值方法,分别选择反距离权重法和普通克里金法对离散仪器烈度值进行空间插值。

通常在插值前需要根据所选的空间插值方法,对样本数据分布特征、内在规律性以及趋势效应等进行分析,验证确定权重系数及相关参数或选择合适的变异函数模型,从而提高空间插值的准确率。在使用地质统计学方法分析时,首先需要对数据是否符合正态分布进行验证,若不符合正态分布,则需对数据进行变换。如图2a为本次地震仪器烈度数据直方图,柱状图整体呈左偏的钟形曲线,显示数据接近于正态分布。正态QQ分布图通常用来评估两个数据集分布的相似程度,图2b显示了该组数据整体延45°方向近似呈直线展布,同样显示数据具有正态分布特征。因而该组数据满足地质统计学要求,可以使用克里金插值方法。

图2 仪器烈度数据直方图(a)与正态QQ分布图(b)
Fig.2 Histogram of instrumental intensity data(a)and Normal QQ distribution map(b)

在最优权重或变异函数理论模型选择方面,使用反距离权重法插值时,通过预览输出和对比验证确定插值最优权重为3。使用普通克里金插值方法时,根据步长及步长数选择标准,一般最大步长取值应该小于研究区域采样点之间最大尺度的一半,并通过确定点与最近的相邻要素之间的平均距离来确定步长大小,取步长值为0.169,同时为了满足步长数乘以步长大小的值是所有点之间最大距离的一半左右的经验关系,确定步长组数为12; 最后,根据数据空间变量的相关性,普通克里金变异函数理论模型分别选择球面模型、指数模型、高斯模型进行空间插值。

为了评价不同方法获得的空间插值结果,本文将空间插值结果按照地震烈度值进行分类显示,以空间插值结果与宏观地震烈度的吻合程度作为评价依据,从而判断插值结果的优劣。本次漾濞地震现场共调查了218个居民点的震害情况,确定了本次地震最高烈度为Ⅷ度,Ⅵ度区及以上面积约为6 600 km²,共涉及大理州6个县(市)。图3为反距离权重法和普通克里金法对应的4种仪器地震烈度空间插值结果与宏观地震烈度的对比图。

如图3a所示,使用反距离权重法获得的仪器烈度影响场与实际宏观烈度边界在Ⅷ、Ⅶ度区的东边界拟合度较好,Ⅵ度区拟合度较差。拟合得到的仪器烈度影响场长轴走向不明显,尤其是Ⅶ度区的长轴走向整体呈北东向,与宏观地震烈度差别较大。在巍山盆地东缘存在一个孤立Ⅶ度异常区,该区域仪器地震烈度值明显高于周围仪器烈度值,与宏观地震烈度有所差别。

对于地质统计学方法,不同变异函数理论模型所得到的插值结果差距较大。如球面模型与高斯模型(图3b、c)对高烈度区的插值效果较差,基本没有拟合出Ⅷ度区的范围,仪器烈度影响场边界与实际宏观烈度边界拟合吻合度较差。但其影响场范围的整体形态与实际宏观调查所获得的长轴走向拟合较好,可以明显看出仪器烈度影响场长轴走向为北西向。

从指数模型的普通克里金插值结果可以看出,在高烈度区该方法依旧可以根据仅有的一个高烈度值拟合出较小的Ⅷ度区,在Ⅶ、Ⅵ度区的东南边界插值结果和实际宏观烈度也比较吻合,在Ⅵ度区西南边界可以看出受预警台站间距的影响,仅在有样本数据的位置上与宏观Ⅵ度区边界拟合度较好,其边界呈花瓣状。其影响场Ⅵ度区与Ⅶ度区长轴方向也较明显,与宏观地震烈度相符(图3d)。

图3 4种空间插值结果与宏观地震烈度对比图
Fig.3 Comparison of four kinds of spatial interpolation results and macroscopic intensity

结合上述讨论的4种插值结果可知,在极震区内台站数量较少的情况下,产出的仪器地震烈度空间分布图更加依赖于插值算法的适用性,不同的插值方法、参数、变异函数理论模型均会对插值结果造成较大影响。本次漾濞地震仪器地震烈度空间插值方法中基于指数模型变异函数的普通克里金方法的插值结果与宏观地震烈度整体范围及特征吻合度较好,反距离权重法插值结果次之。

对仪器烈度空间插值结果与宏观地震烈度空间范围进行叠加分析,如图4所示,在宏观地震烈度Ⅵ度区范围内,存在部分仪器烈度值小于宏观地震烈度的情况。Ⅵ度区内共有预警台14个,其中8个预警台仪器烈度值小于Ⅵ度,共有7个预警台位于宏观地震烈度Ⅵ度区的西北部,仅有1个预警台位于东南部。仪器烈度影响场存在沿长轴方向在Ⅷ度、Ⅶ度、Ⅵ度区东南方向面积普遍较大的情况,其原因可能是在震中东南方向地震动衰

图4 克里金指数模型插值结果与宏观地震烈度空间分布图
Fig.4 The optimal interpolation result and the spatiae distribution of macroscopic intensity

减较慢,导致在主震东南侧记录到的仪器地震烈度普遍比西北侧高,且距震中更远。而实际调查的宏观地震烈度则沿长轴方向在微观震中两侧面积相差不大,其原因可能是仪器烈度数据仅为漾濞M_S6.4主震产生的地震动强度所计算的烈度值,而漾濞地震除主震以外,前震和余震频发,其中5级以上地震就有3次,宏观地震烈度为主震及多次前震、余震叠加所造成的震害,因而在漾濞地震中同一地区宏观烈度值普遍高于仪器烈度值。

在漾濞地震震中附近仅有漾濞台仪器地震烈度大于Ⅷ度,所以导致仅通过仪器烈度插值获得的Ⅷ度区范围较小,与实际宏观烈度Ⅷ度区范围相差较大。其次,在宏观地震烈度Ⅵ度区南部,受预警台站分布不均匀且密度较低的影响,导致在该处空间插值形成的仪器地震烈度影响场Ⅵ度边界呈花瓣状,与实际调查烈度存在一定差异。可见台站密度和分布均匀程度是造成仪器地震烈度与宏观地震烈度差异的另一个影响因素。

由于宏观地震烈度反映的是震害水平的综合指标,影响因素较多,且在综合制图中有平均化的趋势,因而其影响场边界更为均匀光滑,而仪器烈度值则直接由观测的地震动记录计算得到,两者在概念和属性上存在一定差异,因而两者的影响场分布存在差异是不可避免的。

震区属横断山滇西高山峡谷区,地形起伏较大,现场调查发现,震害情况受地形、场地影响较大。通常现场调查是以居民点为单位对该区域的房屋震害整体情况进行评估,而预警台记录的仪器烈度值则反映的是台站所处区域的场地地震动特征。因目前大部分烈度计与通信基站合作架设,其周边地形及场地条件可能往往与居民点所在的环境有所不同,从而造成仪器烈度值与实际调查结果有一定的差异。

在地形起伏缓和、场地覆盖层较厚的地区,场地的放大作用对仪器烈度具有较明显的影响。图4中在Ⅵ度区东南部,巍山盆地东缘大仓台(震中距47 km)仪器烈度值为7.1,而位于其西南侧山区马鞍山台(震中距43 km)仪器烈度值仅为5.9。在地形起伏较大的山区,台站周边场地通常覆盖层较薄,而居民点通常位于平缓的坡地、谷地等覆盖层相对较厚的地方,从而导致仪器烈度值往往低于实际调查烈度值。关坪台位于Ⅵ度区西北部,台站架设在覆盖层较薄的坡地上(图5a),本次记录的仪器烈度值仅为4.1,而在其北侧、东侧约2 km的关坪村旧街小组与自新村实际宏观烈度值均为Ⅵ度(图5b),两者烈度值相差较大。另外,台站与居民点实际所处的地形地貌对震害情况也有着明显的控制作用,往往造成两者烈度值的差异。

图5 关坪台场地环境照片(a)与旧街小组局部震害照片(b)
Fig.5 Surroundings of Guanping station(a)and the earthquake damage to a house in the Jiujie village(b)