0 引言

中国是亚洲乃至世界上地质灾害最为严重的国家之一，其地质环境复杂、板内构造活动活跃(黄润秋，2007)。20世纪以来，我国发生的内陆破坏性地震已占全国的三分之一(李同林等，2024)，且崩塌、滑坡、泥石流等地质灾害频发，对山区群众生命财产及交通安全造成极大威胁(王涛等，2013)。落石是山区常见的自然地质性灾害，随时威胁着当地人民的居住及生活，而确定落石的位置、发生时间、体积和滑动路径等是改善相关危害和更好地预防这些危害的关键。最新研究表明，落石可以根据其产生的地震信号进行高精度的自动检测和定位(Farin et al，2015)。目前，国内外常规的监测手段有卫星图像分析(黄观文等，2018)、三维激光扫描(戴岚欣等，2017；黄海峰等，2017；王明等，2019)，基于点位移监测的GNSS、伸长计和倾斜计等仪器设备，通过监测边坡位移进行滑坡预警(Intrieri et al，2019；Lombardi et al，2017)。然而，这些预警技术大多是针对面或者有限点的位移变化规律，无法满足针对崩塌落石、泥石流等小区域中快速变化移动的地质灾害监测预警的应用。

过去十几年间，研究人员对岩崩等相关微震信号的研究表明，除了事件监测之外，地震信号还可以提供震源的有关信息，例如体积、持续时间和位置(Hibert et al，2017b；Le Roy et al，2019)。在法国阿尔卑斯山软岩(黑色泥灰岩)沟壑中进行的一系列单个落石实验表明，从地震信号可以获取每个石块撞击前的质量和速度，虽然具有很大的不确定性，但是其值接近于真实数值，且微震信号特征与单块岩崩的能量、速度和质量存在相关性(Hibert et al，2017b)。研究学者还发现地震信号的最大振幅与体积之间存在相关性，但仅适用于具有几乎相同跳动路径的事件(Hibert et al，2017a)。有学者在2011年日本台风引起的山体滑坡的微震观测中获得了相似的结果，在距震源1km左右处测量的信号振幅平方积分随滑坡体积的平方而变化(Yamada et al，2012)。

落石物理参数反演计算是由测量数据反推物理模型的过程，即通过微震监测系统采集到的微震信号反算被测区域落石的相关信息(如体积、掉落位置、传播路径)。因此，微震信号可为相关科学研究提供丰富的科学资料(冯甜等，2021)，例如从地震信号计算的参数与破碎岩崩特性之间的经验函数可以推测出岩崩的体积(Lin et al，2022)。有研究发现震源处的体积与地震记录的振幅之间存在线性关系(Norris，1994)。Rousseau(1999)根据数百次岩崩开发了一个数据模型，可从其最大振幅中得出落石的体积。除此之外，可以根据落石的信号推测出分离后撞击的势能和地震能量，反演出撞击时的体积，发现势能和地震能量之间具有良好的相关性(R²=0.98)(Le Roy et al，2019)。相比之下，Saló等(2018)并未发现动力学参数与测量的地震能量之间存在相关性。所有这些研究均表明根据相关地震信号估计的落石特性(体积、速度、位置等)的巨大不确定性。

本研究利用意大利SARA公司的三轴微震监测设备和大疆无人机远程操控设备在广东省大宝山地区进行为期2天的野外人工滚石实验，探讨关于岩崩落石的物理力学过程与相关微震信号之间的联系，以及落石体积反演在不同地质环境下是否可以进行迁移的问题。对13个落石微震信号特征进行研究，获得大量关于落石特征的数据(位置、发生时间、体积、破坏模式等)并对频谱图进行分析。通过计算落石微震能量和势能并进行幂指拟合，对拟合出的系数a、b在上犹县两组实验数据进行验证(即大宝山的数据被用作训练数据集，以确定落石特性与地震信号特征之间的关系，将这些关系应用于上犹县实验)，对上犹县落石体积结果与实际体积进行对比分析，测试大宝山的数据关系是否能够准确适用于较大体积或较小体积的落石量以及不同的地质环境。本研究不仅为崩塌落石等快速冲击型地质灾害远程监控的应用与发展提供了一种新的思路，也为崩塌落石微震监测预警系统提供一定的理论参考。

1 实验场地和实验介绍 1.1 实验场地

本实验的研究区域广东省大宝山露天采矿场是一座正在开采中的富含铅锌的多金属露天矿山，海拔范围300~1068m，相对高差300~700m，地形起伏较大。地形以丘陵盆地、山地为主，山系总体呈南北走向，地势总体上呈现北高南低，矿区位于两个几乎平行的山脊之间的小型向斜盆地中。

具体的研究实验地点为大宝山矿区开挖过的一个边坡，由于矿山开采，导致该岩质边坡非常陡峭，易发生崩塌、滑坡等地质灾害，雨季更容易发生坍塌等现象。由于矿山附近进行岩石的爆破工作，岩爆产生的冲击、该区域降雨量较大均会增加崩塌、滑坡的发生频次，该实验边坡属于岩质边坡，岩性主要为白垩纪燕山期岩浆作用形成的云英闪长岩，岩石抗压强度为81.8GPa，弹性模量为62.09GPa，密度为2.775g/cm³，边坡表面进行了削坡减载处理，由 图 1(b)可以看到明显的台阶状，本次实验选择在第一个台阶处进行，落石的地貌落差约为23.74m，大致为南北走向，倾角为79°，其地理位置和岩石特征如 图 1(a)、1(c)所示。

1.2 实验方案

本次人工投石实验落石样本数量共计27个，岩块取自坡体自身碎裂岩块，保持岩性一致，落石样本有椭圆形、柱状、立方体形、长条形，扁平状，尺寸大小分布在23cm×20cm×6cm至50cm×40cm×43cm之间，体积分布在0.00276m³至0.086m³之间。岩石样本自第一级台阶顶部边缘滑下，最终跌落停止在边坡底部平台，部分岩石样本停留在边坡中间。

选取27个落石(编号N1~N27)进行实验，在实验监测过程中，存在未能推动和推下来但无法记录掉落位置的坐标，最终只采用准确接收到微震信号的13个落石数据做研究样本(表 1)。

本次人工滚石实验时间选取在中午工人停工时间，以获取一个较为安静的实验环境。首先对抛掷的落石进行编号，并记录其原始体积和初始位置。在整个实验过程中，每一次人工投石均采用高分辨率相机拍摄。实验总共布置了9台仪器，仪器布置方向和落石掉落轨迹如 图 2所示。

1.3 实验仪器

本次实验采用意大利SARA电子仪器公司生产的三轴微震监测设备(S45三轴速度计)，用于接收人工投石的振动信号，换能系数为78V/m/s。监测设备的监测方位为东西方向EW、南北方向NS和垂直方向Z，内置全球定位系统(GPS)实时定位授时，确保单个仪器记录的时间准确性。这些仪器的自然频率为4.5Hz，实验中将仪器的采样频率全部设置为200Hz。监测设备的电源来自不同站台间的户外电源(12V-60A电池)。

实验一共布设了9台微震监测设备，其中一台监测设备(TOR3)为传感器和采集器一体式地面型，3台为分离式钻孔型(TOR1、TOR2、TOR4)，另外5台为分离式地面型(TOR5、TOR6、TOR7、TOR8和TOR9)，每台设备配备了一个SARA 24位A/D转换器(SL06)，通过转换器连接到一个24位数字化仪。每个监测设备均通过水准仪调平，使用三角支撑底座放置在地面上，与地面隔离开。为获得准确的落石振动数据，布设方式采用钻孔布设和平面布设相结合的方式。台网的分布为TOR1、TOR2和TOR4放置在钻孔中，TOR3放置在边坡顶端靠近人工投石的位置，其他仪器在边坡底部平面分散布置，仪器的具体安装位置见 表 2和图 2。使用钻机钻取3个钻孔用于安置TOR1、TOR4和TOR2，钻孔直径为15cm，钻孔深度分别为15m、10m和5m(表 2)。需要注意的是，钻孔型微震基站在布设过程中利用砂土填充，使探头在钻孔中保持垂直并与孔壁紧密贴合。

2 微震信号特征

落石事件通常被记录为一系列短脉冲，持续时间为4~11s，而中途撞击次数则受投石位置的地形影响。前人使用震源距离与微震台网相距600m进行实验，发现落石碎屑流的频率分布在2~10Hz，中心频率为5Hz，而单个落石的中心频率为7Hz，最高频率可达40Hz(Allen，1978)。当震源与微震台网距离100m时，信号频率可达到50Hz，但其频率随基站距离的增大，衰减程度也十分明显。本研究中，基站距离与震源位置的距离大约在16~60m，落石微震事件能量大部分分布在10~90Hz。当监测距离60m时信号频率可高达50Hz，不同台站接收到的振动信号也随距离变大而衰减。

基于大宝山落石的微震监测数据分析微震监测落石的频率特征，为展示对整个数据集进行分析的结果，选择了几个具有代表性的落石事件。同时，为保证地震波形不失真，实验选择在中午进行，排除了外界环境引发的振动(如卡车行驶、工人开工等)。通过微震监测设备获取的数据类型为miniseed，将文件导入Geopsy软件(Picotti et al，2017)处理后，再通过Python编程将读取的数据转换为地面速度信号，采用适合分析落石信号的小波变换cagu 8，小波变换的尺度数目100，得到信号的波形如 图 3所示，其中根据落石撞击特点截取N23、N21落石8~14s和1~8s(N23落石前8s与8~9s波形相同，且N23持续时间较长，为了清晰展示后半部分撞击和反弹的波形，故只截取8~14s)、N5落石0~4s这3个具有意义的时间点的波形，再基于快速傅里叶变换将时域内容转换至频域，得到微震信号的频谱图(每个片段中傅里叶变换的数据点数即窗长度为200)。

落石N23是落石发生类型中较为典型的一类，其体积较大、微震记录较为清楚，落石在启动后经历了自由落体，在撞击斜坡中部小平台后，带动坡体上的碎石一起掉落，最终在坡底经过翻滚后停止。落石N21体积较小，形状近似三棱锥，被推动后带动坡上的土屑、碎石滚落，垂直掉落撞击至斜坡中部小平台后破碎成小体积的碎石屑掉落坡底，直至滚落停止，从推动到结束的时间为7s。落石N5为最常见的落石滑落事件，即推动后沿土坡滚动至中间平台停止，之间可能伴随有微小的撞击情况，其波形从开始至结束并无特别之处。从 图 3(b)、3(c)可见落石N21、N5在被推动后一直伴有碎小落石，而 图 3(a)中N23被推动后直接经历自由落体，波动较N21、N5更小。

对微震信号的分析突出了岩崩的复杂性和可变性。通过对地震信号的分析可以确定地震信号的传播方式，例如Vilajosana等(2008)和Bottelin等(2014)研究观察到的能量最高的阶段是由自由落体后的撞击产生的，这些冲击产生频率高达50Hz的脉冲信号，而连续块反弹的能量较小，出现的波形更多，频率高达30Hz。

通过Python环境下的快速傅里叶变换(FFT)和小波变换的分析结果，显示出微震信号三个方向上的频率波动较为一致。由 图 4(a)可以看到，在12~13s内，落石N23的微震信号波速和频率在三个方向上均达到最大值，其中在E、N方向上波速最大值约为1.8m/s，而在Z方向上波速的最大值约为2.5m/s。对微震信号进行小波变换后得到时频分析图(图 5(a))，显示该事件的微震信号能量主要分布在20~80Hz，主峰为17Hz和50Hz，中心频率为50Hz。由于落石N23是最典型的事件，其掉落过程中带动的碎石较少，因此其反演结果也比N11的效果好。

3 参数反演模型构建 3.1 微震能量

Deparis等(2008)、Dammeier等(2011)、Manconi等(2016)和Fuchs等(2018)对落石参数与地震信号性质之间的比例关系进行了分析，在这些研究中落石体积约为10³~10⁷m³，落石距离微震监测站大部分为几十千米，属于大型崩塌落石，落石的体积、落石到基站之间的距离均比本文研究大。在本研究中，大宝山人工投石实验微震监测基站布设密集，且围绕落石发生地点布设，落石发生位置与微震基站距离不超过100m，落石微震信号记录更加准确真实，更有益于进行落石物理参数反演模型构建。

根据大宝山人工滚石获得的微震信号，并根据已有文献(Le Roy et al，2019；Vilajosana et al，2008)计算出落石的能量(微震能量E_s、势能E_p)和体积，假设各向同性均匀传播介质和点力源，以地震表面波的形式耗散的能量可以使用Vilajosana等(2008)假设的表面波主导的地震信号关系式得到，即

$ E_{\mathrm{s}}=\int_{t_1}^{t_2} 2 \pi r \rho h c u_{\text {env }}(t)^2 \mathrm{e}^{\alpha r} \mathrm{~d} t $

(1)

$ u_{\mathrm{env}}(t)=\sqrt{u(t)^2} $

(2)

其中，t₁和t₂分别为微震信号的起始和终止时间，r为落石事件落点与微震监测仪器之间的距离，h为表面波传播层的厚度，ρ为地面的密度，c为地震波的相位速度，u_env(t)是使用希尔伯特变换H_l(t)获得的地震信号的幅值包络线(这里是地表速度)，α为阻尼因子系数，用于解释地震波的非弹性衰减，该阻尼系数与频率有关，公式为

$ \alpha=\frac{\pi f}{Q c} $

(3)

其中，f为该频带中间的频率，因为在该频带的中心可以观察到颗粒流的大部分能量；根据已有文献(Pavlis，2003)的方法，采用双站台法计算出地震波的相速度为1341m/s，计算方法如下：

(1) 布设两台地震仪，间距L已知；

(2) 采集同一事件的表面波记录；

(3) 通过对比两个波形的特征点(如第一个主峰或零交叉点)得到传播时差Δt，公式如下

$ v_{\mathrm{p}}=\frac{L}{\Delta t} $

(4)

Q为考虑地震波衰减的无量纲品质因数(Xia et al，2002)，按照koyanagi等(1992)的研究设置为50。厚度h被视为以所选速度传播的频率峰值的瑞利波的一个波长，根据已有文献(Le Roy et al，2019)选取h=16m。

3.1.1 落石运动中能量的损失

在地震过程中不可避免地会产生能量的损失。地震能量是落石从初始状态到最终状态运动过程中损失的总能量的一部分，流动过程中损失的能量为

$ \Delta E_{\mathrm{p}}=\left(E_{\mathrm{pf}}+E_{\mathrm{cf}}\right)-\left(E_{\mathrm{pi}}+E_{\mathrm{ci}}\right) $

(5)

其中，E_p为势能，E_c为动能，下标i和f表示初始状态和最终状态。由于实验的落石没有最初速度且最终处于静止状态，故

$ E_{\mathrm{ci}}=E_{\mathrm{cf}}=0 $

(6)

因此，在落石下落过程中最终能量的耗散即为势能的损失，式(4)演变为

$ \Delta E_{\mathrm{p}}=m g\left(Z_{\mathrm{Gf}}-Z_{\mathrm{Gi}}\right) $

(7)

假设落石的落点位置为基准面，则坡高h₀可以通过坡长与坡角计算而得，h=0，则有

$ \Delta E_{\mathrm{p}}=m g l \sin \theta $

(8)

其中，l为该实验边坡的坡长，θ为该实验边坡的坡角。

3.1.2 微震势能比

微震能量与势能的关系很难用物理方法验证，其比值Rs/p往往随岩崩的体积大小而变化，研究表明当势能超过10¹⁰J时，Rs/p大致呈斜率变化。对于不同背景下的岩崩，前人分析得出了不同的微震能量与势能比，结果如 表 3所示。

此外，Vilajosana等(2008)对较小的落石(V=70m³)进行了研究，在这种情况下两个台站形成的地震网络距离震源仅有几十米，其发现E_s/ΔE比率为0.25，该值与其他几位学者研究得到的比率值差异较大，表明场地的地震特性以及落石类型的强烈影响。在本文研究中，产生的落石的性质和现场地质条件在此期间均保持不变。

Farin等(2015)和Bachelet等(2018)通过实验发现，Rs/p的变化可能与石块的大小和受冲击介质的性质有关。Farin等(2015)经过一系列实验发现在撞击物质量为0.001~100g时，Rs/p的值在0.001~1之间。

图 6给出9台微震监测设备三种情况下的势能比Rs/p，除最小值波动较小外，最大值和平均值均有大幅度波动。其中，TOR3位于坡体顶部，具有地震信号的放大效应(Kuehnert et al，2020)，剩余8个微震监测基站距离落点的距离较远，每个基站距离落点的距离分别为TOR3：13.65m、TOR6：37.00m、TOR1：37.55m、TOR4：37.97m、TOR8：37.97m、TOR2：38.74m、TOR9：40.08m、TOR5：51.54m、TOR7：66.07m。

3.2 冲击阶段势能与微震能量分析

冲击阶段的微震能量与势能有较强的相关性，势能公式为

$ E_{\mathrm{p}}=\rho g V H $

(9)

其中，ρ为密度，H为自由落体的高度，g为重力加速度，取9.82m/s²。

由式(1)可以计算出冲击阶段的微震能量(图 7)，根据大宝山的落石数据计算其势能，可以得到如下关系式

$ E_{\mathrm{SI}}=\mathrm{a} E_{\mathrm{p}}^b $

(10)

E_SI的体积可以用下式计算

$ V_{E_{\mathrm{SI}}}=\frac{E_{\mathrm{p}}}{g \rho H_{\mathrm{s}}}=\frac{E_{\mathrm{SI}}^{\frac{1}{b}}}{a g \rho H_{\mathrm{s}}} $

(11)

为了得到E_SI和E_p之间的幂律关系，本文仅使用广东省大宝山的数据进行拟合(即以大宝山为数据集)。图 8(a)中蓝色直线为大宝山TOR9的线性拟合结果，红色三角形为上犹县E_SI、E_p的散点图，针对每次落石事件，每个基站拟合得到的系数a范围在(4.18×10^-4±3.79×10^-4)~(0.42±0.397)，系数b范围在(1.64±0.008)~(0.973±0.088)，拟合相关性系数R²在0.85~0.98(由于TOR6基站信号不明显，计算出的能量对体积迁移反演存在影响，故仅采用剩余基站的系数)。Hibert等(2017b)也发现类似的E_S和E_p关系，其研究给出a=4.6×10^-8，b=1.38。依照式(10)将拟合出的a、b系数用于反算广东省大宝山13个人工滚石事件的体积，发现V_ESI(4533~79260cm³)与V_DEM(5472~86000cm³)的结果十分吻合(图 8(b))，两者相关性系数R²=0.98。

4 实例 4.1 实验区域

为验证落石体积反演计算模型的可迁移性，本研究在江西省赣州市上犹县某小区后侧人工切坡和某道路人工边坡分别进行人工滚石微震监测实验(图 9)。

2022年3月2日的实验区域为赣州市上犹县梦想家园城南社区内的一处人工边坡，坡长约7.2m，坡高约5.35m，坡角为48°，坡向230°。此处边坡有一较明显的汇水沟，其成因可能为长期雨水冲刷。边坡由于节理发育导致滑坡变形破坏，且表面几乎完全无植物覆盖，岩石岩性为半风化红棕色泥质粉砂岩(图 9(b))。

2022年3月6日的实验区域为赣州市上犹县迎宾大道旁一处人工边坡，坡长约11.9m，坡高约8.98m，坡角为49°，坡向为228°。此处部分边坡进行了人工格构梁支护，格构内喷射混凝土用于防护。而实验所选取的地址为未支护区域，其边坡表面几乎完全无植物覆盖，坡底修有排水沟，岩石岩性为中风化青灰色细砂岩，且岩体表面较为破碎(图 9(c))。

本次实验分别从两处斜坡处收集了8块体积不同、但岩性与斜坡相同的岩石，其平均密度分别为2.35g/cm³和2.71g/cm³，尺寸范围分别为0.0029~0.031cm³和0.0067~0.045cm³，两次实验的岩石形状均接近于平行四边形且岩性均为细砂岩(图 10(a)、10(b))。

4.2 模型验证

将广东省大宝山人工滚石物理参数反演数学模型用于赣州市上犹县的实验，用于验证不同地质环境下能否准确反演计算上犹县落石的体积。上犹县实验总共选取16个大小形状不一的落石，落石样本均来自实验坡体碎裂岩石，由于天气、地形条件等不利因素，在微震信号监测过程中仅选取了8组有效数据。利用大宝山实验数据计算出E_p和E_s的关系，拟合出8组系数a、b反演上犹县人工落石的体积(即以大宝山为数据集，上犹县为训练集)。仅选取反演得到的最大体积和最小体积进行绘制(图 11)，发现上犹县落石的体积(260~45144cm³)均在大宝山数学模型反演出的最大和最小体积范围内，某些落石体积的低估或高估可能是由于用于计算冲击阶段地震能量的参数所致，如两组地区实验地点与基站距离差异较大，上犹县两组实验基站与落石抛掷点距离较近(约为10m)，而大宝山基站设置较远，频率随距离会产生衰减而高频信号衰减较快，且花岗岩和细砂岩岩石特性、岩体结构也有所不同，这些均有可能对实验结果产生影响。因此，需要根据上犹县两组实验地点的环境进行调整，以考虑波传播介质的差异，以及针对实验地点测量出地震波速等，使模型迁移反演出的体积更具准确性。

5 结论

本研究基于微震监测技术，根据当前微震技术在崩塌落石物理参数反演和监测预警方面面临的问题，通过解译崩塌落石微震监测信号，获得落石的发生时间、体积、速度等物理参数，对落石滚动的几种情形进行研究，分析物理岩崩过程与微震监测信号之间可能的联系，并将其物理参数模型在不同的地质环境下做可迁移性验证，取得如下初步认识：

(1) 结合无人机摄影和三维建模(DSM、DEM)对人工滚石信号进行分析，获得大量关于落石特征的数据(位置、发生时间、体积、破坏模式等)，识别出落石滚落所经历的过程：分离，自由落体，随后是冲击、反弹和碎石滚动。本研究展示了大宝山实验落石启动后坠落的三种形式：出现明显撞击、随坡滚落、撞击加滚落，从大宝山人工滚石微震信号波形图、频谱图中可以发现80%的落石均有明显的冲击、分离、反弹、碎屑流，且冲击阶段产生的频率最大(高达50Hz)，落石反弹阶段频率也高达30Hz。

(2) 根据大宝山落石的微震能量E_s和势能E_p，估算出以地震波形式耗散的落石势能比，与前人研究范围进行比较后发现，大宝山落石的Rs/p范围均在理论范围之中。此外，推导了大宝山数据集微震能量E_s与势能E_p的关系，发现撞击的微震能量E_s与岩崩的势能E_p有较好的相关性(R²=0.98)。

(3) 将微震能量与势能进行幂律拟合得到每组基站a、b系数(作为数据集)，并在上犹县两组人工落石实验上进行了迁移验证。某些落石体积的低估可能由于基站排布的距离(如两者实验地点距离抛掷点相差6~7倍)、岩石的特性以及落石冲击介质差异不同。此外，大宝山落石体积与上犹县选取的落石体积相差较大均有可能对体积反演产生影响。

本文得出的E_p和E_s关系及系数可能仅适用于特定体积的情况，并可能取决于落石传播的基质，且需考虑地形放大效应，对不同的地质环境应有所调整。在考虑上述影响的基础上，本文研究成果可与多源数据相融合，进行崩塌落石物理参数反演，从而准确预测较大或较小体积的落石量。