利用二维动态数值模型评估风化岩石对局部地震危险性的影响：以新加坡武吉蒂马花岗岩为例

《Geomatics, Natural Hazards and Risk》：Assessment of weathered rock effects on local seismic hazard using 2D dynamic numerical modeling: a case study of the Bukit Timah granite, Singapore

本研究针对新加坡Bukit Timah花岗岩区，通过整合HVSR、MASW、MAM及2.5D地质建模，构建了五条二维剖面并开展动态数值模拟，揭示了风化岩层厚度与地震动放大效应（PGA达0.31g）的密切相关性，验证了数值模型在地震微分区中的应用可靠性，为热带岩区地震风险评估提供了新方法。

　　在新加坡的Bukit Timah地区，对地震影响的评估是一项重要的研究工作，尤其关注花岗岩风化层对局部场地效应的影响。本研究通过整合地层、地质图和地质物理数据，构建了一个2.5D地质模型，并从中提取了五个2D剖面，用于动态数值模拟。模拟过程涵盖了共振频率、放大模式、峰值地面加速度（PGA）和剪切波速度（Vs）等关键参数的分析。研究结果表明，主导频率在1.8到3.8 Hz之间，放大系数最高可达4.5，这一结果与高频谱比（HVSR）数据高度一致，突显了风化花岗岩层对场地响应的显著影响。通过将HVSR数据与合成曲线进行反演，我们还构建了一个2D剪切波速度模型，进一步揭示了地层和地形对地震响应的复杂作用，特别是风化层的横向异质性。

### 地质与地震背景

新加坡位于苏门答腊板块上，属于地质稳定性较高的地区，但仍然受到来自苏门答腊活动构造带，如苏门答腊俯冲带和苏门答腊断层带的地震波影响。尽管该地区的地震活动相对较低，但由于其地质构造的复杂性，远场地震如2004年苏门答腊-安达曼地震（Mw 9.1）和2005年尼亚斯-西米卢埃地震（Mw 8.6）在新加坡仍然可以产生可感知的地面运动。这种现象表明，即使在地震活动较低的地区，也必须考虑地震波的局部放大效应，特别是在城市化程度较高的区域。

在新加坡的Bukit Timah地区，花岗岩和沉积岩构成了地表的主要地质单元，其中花岗岩占据了约三分之二的面积。由于热带气候、高降雨量和生物活动的影响，花岗岩表现出广泛的风化特征，从新鲜花岗岩（Grade I）到残余土壤和高度风化的材料（Grade VI）的过渡层，厚度可达70米，平均约为30米。这些风化层富含高岭石、氧化铁和水铝石，呈现出独特的红色土壤特征。风化层的厚度和性质对地震波的传播和放大具有重要影响，因此对这些地层的准确识别和建模对于地震灾害评估至关重要。

为了更精确地刻画地层结构，本研究结合了多种非侵入性地质物理方法，包括多道表面波分析（MASW）、微震阵列测量（MAM）和高频谱比（HVSR）分析。这些方法为估算剪切波速度（Vs）和识别风化层提供了有效的工具。例如，HVSR分析揭示了共振频率与层厚之间的强相关性，有助于风化花岗岩深度的确定。此外，研究表明，阵列几何和场地条件对MAM结果有显著影响，而剪切波速度约为500 m/s时，通常表示中度风化的花岗岩。MASW和MAM的结果也与HVSR数据一致，说明这些方法在刻画地层特征方面具有良好的一致性。

在Bukit Timah地区，地质结构和风化程度对地震放大效应具有决定性影响。该地区的地震活动和构造特征，如Bukit Timah断层，影响了研究区域内共振频率和剪切波速度的分布。由于风化层的存在，不同区域的地震响应特征存在显著差异，这使得传统的场地效应分析方法在缺乏详细地层信息的情况下显得不足。因此，本研究采用了一种结合多种方法的综合策略，以提高对花岗岩风化层影响的刻画精度。

### 方法论

为了评估Bukit Timah地区的地震地面运动和场地放大效应，本研究采用了2D动态数值模拟方法。该方法基于离散元法（DEM），在UDEC（版本6.37）软件中实现。DEM方法结合了有限差分法（FDM）的基本原理，能够模拟静态和动态行为。模型中的介质由离散、绑定的块体组成，这些块体被划分为具有线性或非线性机械特性的有限差分区域。为了减少边界效应并模拟真实的波传播，采用了Rayleigh阻尼，同时设置了自由场边界条件。

在静态建模阶段，模型的构建基于数字高程模型（DEM）数据，使用Shuttle Radar Topography Mission（SRTM）提供的地形信息，并在ArcGIS中绘制剖面图。剖面采样间隔为25米，以确保数据的准确性和一致性。地层几何结构基于地质图、钻孔数据和2.5D地质模型，将材料分为残余土壤（Mat 1）、中度风化花岗岩（Mat 2）和新鲜花岗岩（Mat 3）。岩体属性，如剪切波速度（Vs）、密度（ρ）、泊松比（ν）等，均根据地层类型进行赋值。

在动态建模阶段，为了确保地震波模拟的准确性，必须选择合适的有限差分（FD）区域大小。根据Kuhlemeyer和Lysmer（1973）的研究，FD区域长度（Δl）应小于最短波长的十分之一，以避免数值色散。例如，对于剪切波速度为230 m/s和频率为8 Hz的软材料，FD区域长度约为5米。在深度方向上，采用较粗的网格划分以优化计算效率，同时保持模型的准确性。相邻FD区域和接触面的匹配对于数值稳定性至关重要，能够减少反射伪影。

输入的地震波包括两个Ricker波形（Ricker, 1953），中心频率分别为1.4 Hz和3.5 Hz。这些波形用于模拟波传播和块体位移，是地震响应分析的关键参数。为了减少边界反射，采用了静边界和自由场边界条件，以模拟更真实的波传播过程。网格划分遵循Kuhlemeyer和Lysmer的建议，以确保数值稳定性并准确模拟波传播。

### 数据处理与结果

在Bukit Timah地区，通过HVSR测量获得了丰富的场地响应数据。这些测量在五个剖面上进行，共包括10个站点（BT1–BT10），每个站点使用三轴4.5 Hz地震仪进行30分钟的记录。大部分站点靠近钻孔位置，以确保校准和一致性。数据处理使用Geopsy（版本3.5.2），通过Konno–Omachi算法（Konno和Omachi, 1998）对环境信号进行滤波和平滑，平滑常数为40。基本频率（f?）和其标准差通过在时间窗口内对频谱峰值进行平均计算，而峰值振幅（A?）则作为阻抗对比度和数据可靠性的代理指标，根据SESAME（2004）指南进行估算。分析范围为0.5–15 Hz，大多数站点的主导峰值位于1.8–5.3 Hz之间。

基本频率为场地响应提供了重要的线索，反映了阻抗对比度和层厚。在西北和中央剖面（BT1, BT4–BT6）中，较低的基本频率（1.8–2.5 Hz）表明深层风化基岩，伴随较低的放大系数（A? < 3）。相比之下，东南剖面（BT7, BT9, BT10）记录了较高的频率（3.7–5.7 Hz），表明较薄的沉积层和较高的放大系数（A?达到6–8.5）。值得注意的是，BT9站点虽然显示出较低的基本频率（f? = 2.5 Hz），但其A?值约为8.5，表明在深层阻抗对比度下，仍存在显著的局部场地效应。类似的差异也出现在BT3和BT8站点，进一步突出了风化层和阻抗的横向异质性。

通过三种主题图，我们展示了HVSR参数的空间变化：(i) 基本频率（f?；图4a），(ii) 峰值振幅（A?；图4b），以及 (iii) 基本频率的极化与方位趋势（图4c）。在这些图中，颜色代表基本频率的类别，符号大小与峰值振幅相关，对角线指示波的方向性。这些可视化结果揭示了所有剖面中的显著横向变化，为后续剪切波速度建模提供了重要的参考信息。

在HVSR和2D SSR分析的椭圆率曲线反演过程中，我们采用了61个合成数据集，通过在五个剖面中进行2D SSR动态数值模拟，以解决反演问题并比较2D SSR结果与HVSR数据之间的差异。在NW–SE方向的剖面1中，软沉积物的厚度呈现出一致的趋势：中央部分（BT4–BT6）显示最大厚度，向西北（BT1–BT3）和东南（BT7, BT9）逐渐减小。这一模式得到了2D SSR反演和钻孔数据的验证。图5比较了Borehole B14、HVSR数据（BT7）和2D SSR模拟（P1.24）的结果，显示了不同方法之间的良好一致性。

通过整合MASW和MAM结果与SPT-N基础数据，我们估计了平均剪切波速度（Vs）。在0.5–15 Hz的频率范围内，HVSR和2D SSR曲线均进行了分析。初始模型参数基于Abdialim等（2021）的研究，使用了每剖面的已发表密度和平均速度。在本研究中，这些参数进行了优化以提高分辨率，如表2所示。之前的研究表明，联合反演瑞利波色散曲线和HVSR曲线可以减少仅基于HVSR的反演中常见的模糊性。

在剖面1上，BT7和P1.24的HVSR反演结果显示出三个层状结构，其下为半空间。对于BT7，半空间边界位于约70米深度，反演模型误差为0.23，使用了192,603个潜在模型。对于P1.24，最佳拟合模型同样显示三个层状结构，半空间边界位于约70米深度，反演模型误差为0.34，使用了171,432个潜在模型。值得注意的是，HVSR结果通常将半空间边界置于2D SSR结果稍浅的位置。Borehole B14数据（图5a）显示了四个层状结构（表1），而反演将其简化为三个主要单元：软沉积物、风化花岗岩和花岗岩。

HVSR曲线显示出明显的峰值，表明基岩的共振效应。对于BT7，峰值出现在1.9–2.5 Hz，而对于P1.24，峰值出现在2.3–2.8 Hz。表面层厚度（h?）的估算采用公式：h? = Vs / (4 × f?)，其中f?为共振频率，Vs为平均剪切波速度。计算结果显示，BT7的表面层厚度约为30米，P1.24约为26米，而Borehole B14约为28米。这些结果在剖面1上显示出一定的空间差异，表明不同位置的地质特征存在显著变化。

在某些情况下，HVSR和2D SSR结果之间存在微小差异，这主要归因于两个因素。首先，地表地形和地形变化引入了局部效应，这些效应在2.5D地质模型中未被完全捕捉。其次，HVSR测量点与2D SSR合成接收器之间存在空间偏移，这可能导致模型结果的不一致。然而，尽管存在这些差异，比较结果（图5c和图5e）表明两种方法在共振频率和振幅方面高度一致，这些差异更多地反映了局部地质变化而非方法论问题。

### 地质模型与动态模拟

为了进一步刻画Bukit Timah地区的地下结构，我们采用2.5D地质建模方法，结合了10个HVSR测量点、4个MASW剖面、3个MAM阵列和17个钻孔数据（图1b）。表面层厚度通过公式（5）进行计算，该公式仅适用于具有明确共振峰的HVSR记录，如SESAME（2004）建议的那样。剪切波速度（Vs）来自Moon等（2019）和Ku等（2021）的MASW和MAM数据，而HVSR峰值频率则与附近的钻孔地质剖面进行了交叉验证。

地质模型构建在一个8×10公里的框架内，整合了2.5D数据。地层边界和结构特征从地质图中数字化，并在剖面之间进行横向插值，以提供比传统2D模型更广泛的空间表示，同时避免了全3D建模所需的大量数据。最终模型包括三个主要地质单元：软沉积物（黄色）、风化花岗岩（浅灰色）和花岗岩（深灰色），如图6所示，展示了五个解释的地质剖面。尽管输入数据密度有限，该模型仍能有效捕捉研究区域的关键地下特征，并为后续的2D动态数值模拟提供了坚实的基础。

在某些区域，由于测量点稀疏或分布不均，模型存在一定的不确定性。横向插值可能造成边界定位不准确，因此在建模过程中，优先考虑具有明确HVSR峰值和一致钻孔标记的区域。然而，这些外推边界应被视为指示性，随着新数据的出现，应进行进一步的修正。

### PGA分布与动态模拟结果

本研究进一步分析了Bukit Timah地区五个剖面的PGA分布，以评估地震风险。通过2D动态数值模拟，我们发现不同剖面的地震波传播特性存在显著差异，这与地层和地形条件有关。模拟数据与现场测量结果相吻合，确保了模型与实际观测的一致性。

在动态模拟过程中，使用了61个合成接收器和五个参考接收器（位于模型底部），以评估模型的频谱特性。通过比较地表加速度记录与参考接收器的相应数据，我们获得了更精确的地震响应信息。这些记录随后被加倍并进行滤波（图S3，补充材料）。对于每个地表接收器，创建了一个包含两个水平分量和一个垂直分量的合成数据集（图7），用于HVSR分析。通过计算地表数据与参考数据的频谱比，我们得到了清晰的放大和衰减模式。例如，图7显示在0.6 Hz附近存在明显的放大效应，而2.5 Hz附近则表现出显著的衰减。这些频谱模式与HVSR结果进行了比较，以评估两种方法的一致性。

通过将模拟结果与HVSR数据进行对比，我们发现PGA值在某些区域显著增加，特别是在软沉积物层较厚的区域。例如，在剖面1的P1.7和P1.24站点，PGA值分别达到2.5和2.8 g，而在剖面3的P3.3和P3.6站点，PGA值进一步增加至3.0和3.1 g。剖面5的P5.3和P5.6站点的PGA值则相对较低，约为2.2和2.4 g。这些结果表明，软沉积物层的存在对地震波的传播和放大具有重要影响。

在Bukit Timah地区的西北部，P1.2站点的PGA值为0.25 g，这表明该区域的地震响应相对较弱。而在西北与剖面2的交界处，P1.3站点的PGA值为0.24 g，低于剖面2的P2.3站点的0.26 g。在剖面1和剖面3的交界处，P1.4和P2.9站点的PGA值几乎相同，均为0.24–0.25 g，这可能与它们位于相似厚度的软沉积物层有关。进一步沿剖面1分析，P1.7和P3.3站点的PGA值从0.27 g增加至0.29 g，表明随着剖面向东南延伸，地震响应逐渐增强。这种趋势在剖面1的东南部分更加明显，P1.22、P1.24、P4.4、P4.5、P5.3和P5.6站点的PGA值范围为0.26–0.29 g，显示出明显的地震响应变化。

在西南部，剖面1与剖面4和剖面5的交界处，P4.4、P4.6、P5.3和P5.6站点的PGA值相对较低，表明该区域的地震响应较弱。然而，P5.6站点的PGA值为0.28 g，与东南部的其他站点一致，这表明在东南部区域，软沉积物层的存在导致了较高的地震响应。这些结果表明，软沉积物层的厚度和性质对地震波的传播和放大具有显著影响。

### 结论

本研究通过整合被动地震技术（HVSR分析和椭圆率曲线反演）、钻孔数据和2D动态数值模拟，展示了评估风化花岗岩地区地震场地效应的有效方法。在新加坡的Bukit Timah花岗岩体中，我们准确刻画了共振频率、剪切波速度（Vs）剖面、剪切波速度（Vs30）分布和局部PGA，突出了风化岩和地形放大对地面运动的重要影响。研究的创新点在于同时使用HVSR导出的参数和2D动态模拟，提供了详细的地下模型和地震响应预测，无需大规模钻探即可实现微区划分和抗震土地利用规划。

尽管本研究采用的2D数值模型在某些方面存在局限性，例如无法完全模拟复杂的3D波传播，以及由于现场数据分布稀疏导致的插值不确定性，但这些方法为提高地震风险评估的精度提供了可行的途径。未来的研究应致力于发展全3D地质模型，增加现场测量密度，并结合真实地震记录以实现更精确的校准。同时，应特别关注地层变化对地震放大的影响，以进一步提升地震风险评估的准确性。

总之，本研究为提高新加坡及类似地质条件下城市地区的地震风险评估和抗震规划提供了重要的数据支持和方法论框架。通过结合多种地质物理方法，我们能够更全面地理解场地效应和地震响应特性，从而为地震灾害防控提供科学依据。