综合利用多道瞬态面波和微动探测分析斜坡内部结构

李巧灵，张辉，雷晓东，李晨，房浩，关伟，韩宇达，赵旭辰

(1.北京市地质勘察技术院，北京 100120； 2.中国地质科学院，北京 100037; 3.中国地质环境监测院，北京 100081)

通讯作者： 雷晓东(1983-)，男，正高级工程师。Email:lei-xiaodong@139.com

0 引言

杭州市山地丘陵发育，属亚热带季风性气候，雨量充沛，岩石风化作用强烈，斜坡表面多覆盖有残坡积层，裂隙发育，在雨季极易发生滑坡等地质灾害。这些潜在的不稳定斜坡严重威胁着当地居民的生命财产安全，因此，查明重要斜坡地层结构，有针对性地建立灾害预警预报系统意义重大[1]。滑坡体与稳定山体之间存在密度、速度和电性等物性差异，为利用地球物理手段调查斜坡稳定性提供了物性基础。常用的地球物理手段包括地震反射和折射法、面波法、电法及电磁法等。其中面波法包括主动源面波法(多道瞬态面波法)和被动源面波法(微动法)[2-3]。

地球物理手段在世界上很多易滑坡地区都取得了成功的应用。李来喜[4]利用高密度电法和地震折射层析成像法得到甘肃某地目的位置的物性断面，为主滑动面的确定提供了有效依据；Guénolé Mainsant 等[5]通过环境噪声相关技术对瑞士阿尔卑斯山附近某易滑坡地区进行了监测，发现该滑坡体滑动前几天S波速度持续性急速降低；Sherrod L等[6]综合利用电阻率成像和探地雷达确定了纽约Little Porter山活动滑坡内部冰川形态和基岩界面等结构；Shan等[7]综合利用无线电大地电磁(RMT)、电阻率层析成像(ERT)和高分辨率反射地震确定了瑞士西南部某滑坡体的几何特征和物理性质；吕擎峰等[8]基于滑坡体和稳定岩体之间的电阻率和速度差异，利用高密度电法和瑞利面波法确定了滑动面的深度和滑动范围，为滑坡治理提供了依据；Petronio L等[9]利用纵波和横波反射地震及面波勘探对意大利东北部的一个滑坡体进行了探测，获得了滑坡体地表至滑动面深度的2D地质模型和速度参数；Xu X Q等[10]利用双源面波法(DSSW，即综合利用主动源多道面波分析(MASW)方法和被动源微动探测(MSM)，以同时保证大的勘探深度和高的勘探精度)得到四川某滑坡内部地层结构；Wang F等[11]综合利用自电位测量(SP)和微动探测评估了日本某滑坡坝的内部结构以及是否存在管道泄漏。除利用速度信息进行滑坡体内部结构探测外，也可以通过面波信号的频率特征指示土层破损程度和场地共振特性的方位变化[12-13]。

终上所述，地球物理探测不仅可以有效探测滑坡体内部结构，对滑坡机理进行评价，还可以对易滑坡地区进行内部结构探测或长期监测，为滑坡灾害预警提供依据。电(磁)法受山区地形影响大，施工效率低，探测浅部地层结构时易受降雨天气影响，故本次选择面波法。在杭州市余杭区鸬鸟镇残坡积层浅覆盖区，选择典型斜坡布设剖面，通过主动源和被动源两者相结合的方法探测斜坡内部结构，并对比分析两种方法的适用性和有效性，为残坡积层覆盖区地质灾害调查技术方法体系建立进行了探索性研究。

1 方法原理

面波是体波与自由界面或分层介质的弹性分界面相互作用产生的一种弹性波，满足物理上的自由边界条件[14]。在层状介质中，瑞利面波的速度是频率的函数，随频率变化而变化。在不同介质分界面处，频散曲线会出现拐点、斜率变化、频散点疏密变化等特征，该特性是瑞利波勘探的理论基础。20世纪50年代，Haskell利用矩阵计算了瑞利面波在层状介质中的频散曲线，成为天然面波在探测地层内部结构和人工面波在工程勘察等方面的基础；后来又陆续发展了Schwab-Knopoff法、δ矩阵法、Abo-Zena法等。20世纪80年代，面波勘探被引入中国，取得了很大的发展。肖柏勋、李庆春、刘庆华等分析了瑞利面波勘探技术应用情况及面波分析常用方法的优缺点[14-16]。主动源面波勘探从稳态到瞬态再到多道面波勘探，在方法应用和仪器研发等方面都有了很大的进步[16-17]，主要用于工程勘察方面。被动源面波勘探也称作微动和环境随机振动，早期主要是利用长周期信号研究地壳结构等深部信息和用于地热资源勘查等方面[18-19]，近年来，许多学者在利用微动探测进行城市地质和环境地质调查等方面做了大量研究[20-24]。

1.1 多道瞬态面波

瞬态面波法通过锤击、落重乃至炸药震源，产生一定频率范围的瑞利波，再通过振幅谱和相位谱分析把记录中不同频率的瑞利波分离出来，从而得到频散曲线。目前常用的是多道瞬态面波法。多道瞬态面波在工程场地勘察和滑坡等地质灾害调查方面有广泛的应用，如：识别溶洞以及回填土和膨胀土等土层分类[25]；划分斜坡结构；探测滑坡体内部S波速度结构，确定滑坡体岩土分层、软弱带和滑动面，为滑坡体稳定性评价提供依据等[24,26-27]。

多道瞬态面波一般采用单边或双边激发、多道拾震器线性排列接收。数据处理多采用频率—波数法[17]，该方法通过傅里叶变换，将时间—空间域记录信号转换到频率—波数域处理。某排列原始记录信号用f(x,t)表示，傅里叶变换后有：

(1)

式中：k为波数，ω为角频率，x为拾震器到震源的距离，t为时间。由于接收的地震信号中，瑞利波能量占比大于60%，根据F(ω,k)能量谱上极大值对应的频率—波数对，就可以根据式(2)求取频率—相速度对：

(2)

1.2 微动探测

微动探测不需要震源，分为单台和阵列2种方式。单台微动又称为地脉动，主要利用单台三分量拾震器接收信号，通过微动信号的频率振幅特征分析进行场地评价等。阵列微动则是利用多台拾震器按照不同的采集阵列同时接收信号，经处理提取频散曲线，反演得到不同深度的地层结构信息。近年来，微动探测在地质灾害调查方面也取得了很好的应用效果，如：探测地铁沿线的孤石；探测煤矿陷落柱等地质异常体[20-21]；探测深部隐伏岩溶覆盖层结构[22]；反演滑坡体3D模型，利用H/V计算场地主频，为滑坡稳定性评价提供参考等[28-29]。

微动探测频散曲线提取方法主要有空间自相关法(SPAC)、拓展空间自相关法(ESPAC)、频率—波数法(F-K)等[22]。SPAC法是用中心点与各半径圆周上各拾震器采集的信号进行相关，再进行方位平均，得到空间平均相关函数，最后经Hankel转换，得到其功率谱表达式[30]，记作：

(3)

其中：c(ω)是相速度；h(ω)是信号垂直分量功率谱；ω是角频率；r是中心拾震器到圆周拾震器之间的距离，即半径；J0是第一类零阶贝塞尔函数。通过中心频率为ω0的窄带滤波器，则有：

(4)

(5)

2 研究区概况

研究区位于杭州市余杭区西北部鸬鸟镇全城坞村附近(图1)。该地区植被发育，以竹林为主，地形坡度多在20°～35°之间；残积和残坡积层广泛分布于山体顶部、斜坡及坡麓表部，岩性主要为含砾石黏性土及碎石土层等，厚0.2～2.5 m。

图1 研究区地形地貌Fig.1 Surface wave field data acqusion location

研究区典型斜坡东侧，坡底紧邻房屋，曾做过地质灾害工程处理。从现场踏勘情况看，该区山坡较陡，与房屋紧邻且残坡积层易滑脱，仍有发生地质灾害的可能。2019年9月6日强降雨过程中，该区域附近发生过2处崩塌、1处滑坡和1处山洪，是重点地质灾害调查区域，故选为本次方法分析的典型斜坡。研究区整个斜坡单元植被茂密，多为竹林，工程地质分层从上到下依次为表层耕植土、全风化层、强风化层和完整基岩。本次在典型斜坡上顺着高程等值线方向布设综合测线AB，同时进行多道瞬态面波和微动探测(见图1)。

3 数据采集

为了方便进行数据对比分析，在典型斜坡上布设了综合测量剖面AB，走向NNE，测线总长约145 m，微动探测11点，多道瞬态面波33点。

多道瞬态面波一般采用单边或双边激发、多道拾震器线性排列接收，图2g为单边激发，12道接收的等间距线性排列示意。采集参数根据地质任务给出，并根据现场试验适当调整。检波器频率选择根据勘察深度要求，利用f=v/λ和h≈λ/2估算，其中f为检波器频率，v为地层面波速度，λ为波长，h为探测地层的深度。检波器固有频率越高对低频部分的压制越强，在满足分辨率的情况下，应尽量采用固有频率较低的拾震器[24]。偏移距的大小需考虑炮检距和震源弹塑带的范围，偏移距过小，干扰波能量强，偏移距过大，远道有效波能量弱；因此，偏移距选择应根据现场试验确定，目的是在数据处理时直达波、反射波和瑞利波能够有效分离[32]。本次主动源瞬态面波采用北京市水电物探研究所生产的SWS-6A工程勘探与工程检测系统，使用人工锤击的震源激发方式，并在激震点铺设垫板，单边激发，拾震器主频为4 Hz，24道采集，道间距1 m，最小偏移距 5 m，全排列滚动。

a～f—微动十字阵列半径组合示意；g—瞬态面波单端放炮多道采集排列示意a～f—radius combination diagram of microtremor cross array；g—transient surface wave array with single shot and multi-trace acquisition arrangement图2 面波数据采集排列示意Fig.2 Surface wave field data acquisition array diagram

4 数据处理分析

瑞利面波频散曲线提取的常用方法有空间自相关法(SPAC)、表面波谱法(SASW)、频率—波数法(F-K)、τ-p变换法、相位展开法等[15-22]。本次多道瞬态面波选用F-K法，微动探测选用SPAC法。剖面资料处理需先按常规面波勘探方法对野外采集的资料进行单点处理。多道瞬态面波单点数据处理流程一般为：切除干扰波—拾取面波—频率波数域谱分析—提取频散曲线—反演[26]；F-K域分析时，可通过二维滤波突出基阶面波能量，确认基阶面波频散曲线。微动探测单点数据处理流程一般为：数据预处理—计算空间相关系数—第一类零阶贝塞尔函数拟合—提取频散曲线—反演；采用最小二乘法进行贝塞尔函数拟合，不能呈现第一类零阶贝塞尔曲线形态的自相关系数曲线应舍弃，以免影响最终频散曲线的求取精度[23]。图3为微动探测数据W02部分频率点的贝塞尔函数拟合，从图中可以看出，频率低于15 Hz时，各半径相关系数与第一类零阶贝塞尔函数具有很好的拟合关系，误差小于0.01；频率高于15 Hz的拟合精度有所降低，推测是因为高频信号随机性强且易受噪声影响。

图4为部分微动测点与同位置多道瞬态面波测点频散曲线对比。从图中可以看出，2条频散曲线总体形态相似，但多道瞬态面波频散曲线相同相速度对应频率整体较微动探测高。如果直接利用半波长经验公式进行深度转换，则会出现相同深度多道瞬态面波探测相速度高于微动探测的情况，故需结合钻孔数据，对深度转换系数进行拟合。具体拟合公式为H=αV/2f；式中，H为勘探深度，V为相速度，f为频率，α为校正系数。本次微动探测深度校正系数为0.5，瞬态面波校正系数为1。

图3 W02微动探测点贝塞尔函数拟合Fig.3 Bessel function fitting of microtremor survey point W02

图4 部分微动测点与对应多道瞬态面波测点频散曲线对比Fig.4 Dispersion curves of microtremor survey and transient surface wave exploration at the same location

根据拟合后的转换系数，可以把频散曲线从频率域转换到深度域，直接插值得到相速度等值线剖面，如图5所示。剖面走向NE约18°，过两个已知钻孔ZK1和ZK2，钻孔数据详见表1。ZK1终孔深度为16 m，终孔岩性为中风化花岗岩，未对全风化和强风化岩进行划分。ZK2终孔深度12.5 m，终孔岩性为中风化凝灰岩。除ZK1外，在该斜坡其他位置搜集的钻孔资料揭露的基岩岩性皆为凝灰岩，故后续解释中基岩岩性统称为凝灰岩。图5a为多道瞬态面波探测瑞利波相速度等值线剖面，1～30号点为等间距排列，间隔4 m，30～32号点之间由于施工环境限制，间隔调整为12 m。从图中可以看出，完整基岩面埋深总体趋势为南深北浅。图5b为微动探测瑞利波相速度等值线剖面，共11个点，尽量以10 m的点距布设。由于微动每个单点都需要多台拾震器按阵列形态和半径组合布设采集，对采集地面环境有一定要求，阵列内场地高程变化应尽量平缓，故部分点距稍有变化。从图中可以看出，多道瞬态面波与微动探测得到的瑞利波相速度变化趋势具有很好的一致性。含角砾粉质黏土和全风化凝灰岩之间的界面对应的瑞利波相速度约为300 m/s，强风化和中风化凝灰岩之间的界面对应的瑞利波相速度约为600 m/s，2个钻孔资料1个划分有强风化层，1个没有区分，仅参考ZK2数据认为强风化层顶界面对应瑞利波相速度约为400 m/s。

表1 研究区钻孔资料统计

a—多道瞬态面波探测瑞利波相速度等值线断面;b—微动探测瑞利波相速度等值线断面a—phase velocity profile of transient surface wave exploration;b—phase velocity profile of microtremor survey图5 AB测线瑞利波相速度等值线断面Fig.5 Phase velocity profile of AB line

综上所述，多道瞬态面波和微动探测具有很好的分层性，能有效探测残坡积层覆盖斜坡的内部结构。

除直接将频散曲线转换到深度域进行分析外，也可以根据拐点、斜率、频散点疏密以及不收敛段的起始拐点等特征确定初始模型，经反演计算S波速度结构。为避免反演过程中初始模型设置等人为因素的影响，更客观地反映地层的岩性及构造变化，也可以直接根据经验公式计算视S波速度。视S波速度剖面可以突出显示具有较大速度差异的界面，通常也对应地质分界面[21,23]。图6为AB测线视S波速度剖面。由于多道瞬态面波点距更小，测点更密，能更好地反映地层横向变化细节，故本次视S波速度等值线剖面是基于多道瞬态面波数据计算的。从图中可以看出，斜坡沿AB测线横向变化较大，含角砾粉质黏土厚度变化范围0～4 m，全风化凝灰岩厚度变化范围1.5～6 m，强风化凝灰岩厚度变化为2～20 m。总的来说，测线北部基岩完整性较好，南部风化层较厚，R09～R13点之间的中风化界面以上岩层厚度最厚，超过20 m。

5 讨论

瑞利波的频散特性和在介质中的传播速度与岩土介质的物理力学性质有关，同一波长的面波特性反映地层在水平方向上的变化，不同波长的面波特性反映地层介质在竖直方向上的变化[24]。当岩土层介质较为均匀时，瑞利波的相速度随深度的加大而线性增加；存在不同介质分界面时，频散曲线会出现拐点、斜率变化、频散点疏密变化等特征，由此可推断地下岩性变化的分界面[17]。主动源与被动源面波联合勘探通常是将主动源的高频信息与被动源的低频信息进行拼接，以兼顾浅层和深部探测效果[10,31]。本次是综合利用主动源(多道瞬态面波)和被动源(微动探测)勘探形成综合剖面，对比分析两者在残坡积层覆盖斜坡内部结构探测中的应用效果，结果表明两者得到的瑞利波相速度在整体趋势上具有很好的一致性，且各岩土层相速度范围与庄师柳等在佛山飞鹅岭探测滑坡体内部结构总结的数据有较好的对应，即：坡积层和全风化层对应面波速度为40～300 m/s；强风化层对应面波速度为300～600 m/s；中风化层对应面波速度为600～1 200 m/s；微风化层对应面波速度高于1 200 m/s[26]。

图6 AB测线视S波速度等值线断面Fig.6 Apparent S-wave velocity profile of AB line

与钻探等方法相比，面波勘探误差主要来源于反演、介质横向变化以及波速分界面与介质分界面的对应性等[9,15]，为了降低误差，可从提高反演精度和复杂地质条件下精确成像等方面进行尝试。在反演S波速度结构时，可根据实际情况选择合适的算法，如瑞利面波反演常用的Thromson-Haskell法，是基于平面波理论导出相邻两界面的传递矩阵公式，在高频计算时会出现数值溢出和精度降低的现象；Schwab-Knopoff法可以避免高频计算时精度降低，具有快速稳定的特点；Abo-Zena法避免了高频数据的不稳定性，且对存在低速夹层的介质分析很有效[15]。除了反演算法的选择外，还可以综合利用高阶面波和H/V谱比特征进行联合反演，以减少反演的多解性[33]。为避免反演带来的多解性，也可以避开反演过程，直接由剖面上各点面波频散曲线计算得到剖面，直观地反映地层弹性结构特征，或者根据经验公式计算视S波速度，为划分地层内部结构提供客观依据[20-24]。多道瞬态面波和微动探测求取的都是采集阵列下方介质的平均速度，当介质物性参数横向变化剧烈时，会影响探测效果。尹晓菲等基于波场分离技术提出一种方法，利用相位扫描互相关获取任意两道间的面波走时，在此基础上求取频散曲线[34]，并认为该方法可提高浅地表面波勘探水平分辨率，但该方法抗干扰能力值得商榷。对于倾斜界面、复杂地形以及物性参数横向变化剧烈等复杂介质情况，瑞利面波的传播特征还有待进一步研究[15]。

主动源面波勘探对浅部分辨率较好，但很难采集到低频信号；锤击震源的勘探深度通常在25 m左右，被动源面波在提取低频信号方面具有很好的优势。传统的空间自相关法(SPAC)对采集阵列和场地有较高的要求，不少学者提出了改进的方法，如扩展空间自相关法和背景噪声互相关技术(NCF)等。本文采用的十字阵列较常用的嵌套三角阵列施工效率更高，适用范围更广，尤其是在沿城市道路或在有茂密植被覆盖的场地施工时，优势更为明显。NCF对拾震器数量和布置没有要求，尽量保证射线分布均匀即可。在相同拾震器数量下，能更好地覆盖探测区域，降低三维面波勘探成本[35]。环境噪声成像技术目前多用于分析地壳上地幔结构和地震灾害调查中[36-37]。除了以面波为有效信号的面波勘探技术外，传统地震勘探也不再只把面波当作噪声去除，而是加以利用，解决复杂地表静校正和近地表结构等问题[16,38]。

综上所述，面波勘探在地层内部结构调查中有很好的应用效果，但也有很多方面需要进一步研究，例如：当介质具有较强的各向异性，如何提高面波探测的精度；如何减少S波速度结构反演的多解性，提高反演精度；对于微动信号微弱的区域，如何有效提取面波频散信息；如何去除体波和其他面波的干扰；如何有效识别复杂地质情况下面波发育特征等。

6 结论

1)面波勘探具有很好的分层性，能有效区分残坡积层覆盖斜坡内部结构。本次选取的典型斜坡内部结构沿AB测线横向变化较大，测线北部基岩完整性较好，南部风化层较厚。含角砾粉质黏土厚度变化范围为0～4 m，全风化凝灰岩厚度变化范围为1.5～6 m，强风化凝灰岩厚度变化范围为2～20 m。

2)面波勘探在深度转换时多采用半波长经验公式，具有一定的误差，需结合研究区钻孔等已有资料进行校正。多道瞬态面波和微动探测的转换系数不一定相同。

3)多道瞬态面波和微动探测计算得到的瑞利波相速度变化趋势具有较好的一致性。表层含角砾粉质黏土和全风化凝灰岩分界面对应的瑞利波相速度约为300 m/s，全风化凝灰岩与强风化凝灰岩分界面对应的瑞利波相速度约为400 m/s，强风化凝灰岩与中风化凝灰岩分界面对应的瑞利波相速度约为600 m/s，与钻孔资料基本吻合。

4)微动探测空间自相关函数与第一类零阶贝塞尔函数拟合关系中，频率低于15 Hz时，各半径相关系数与第一类零阶贝塞尔函数具有很好的拟合关系，误差小于0.01，高于15 Hz时，拟合精度有所降低。推测是因为高频信号随机性强，且易受噪声影响。