分布式模糊聚类微动法铁路路基岩溶地球物理探测:以皖赣铁路宁国改线工程为例

王其合, 苏本玉, 王国林

(1.中铁上海设计院集团有限公司, 上海 200070; 2.中国矿业大学资源与地球科学学院, 徐州 221116)

随着中国既有铁路升级改造和邻近既有铁路新建铁路的推进,地表浅层的地质问题愈发突出,尤其是在路基段浅地表岩溶发育地区。岩溶发育区常伴有溶沟、溶隙、土洞等,绝大部分岩体完整性遭到破坏,使得岩石强度、稳定性大大降低。岩溶地区在基础施工过程中,可能出现严重漏浆、地下水位突然显著变化、溶洞顶板坍塌等现象,造成周围岩土体稳定受到严重影响,从而引发地面沉降甚至塌陷等病害,形成次生灾害。微动勘探法是一种高精度非侵入式地球物理探测新方法[1],它利用天然源获得地下介质结构,具有较强的抗干扰能力且不影响既有线的运营安全。该方法被广泛应用于城市轨道交通、地下空间开发、水利工程、地下河道、城市塌陷等领域的岩溶探测。刘铁华[2]利用综合微动技术从覆盖层厚度和异常规模两个角度对昆明地铁5号线地质勘探中的成果精度进行量化分析,提高了探测精度。张陈等[3]采用微动探测方法对成都某批发市场地表建筑密集区进行地下浅表2～40 m精细分层,有效支撑了地铁建设。何军等[4]采用综合物探方法对武汉市隐伏岩溶区地下空间进行探测,提出微动对岩溶发育程度、土层厚度及分层均有较好效果。汤克轩等[5]采用瞬变电磁法和微动勘探法对云南某水利工程进行区域岩溶发育程度探测,依据物性参数分布特征总结可溶岩地层地质解译判定依据。梁东辉等[6]在岩溶区成功使用微动HVSR(horizontal-to-vertical spectral ratio)法对地下管线及岩溶进行探测。范长丽等[7]利用微动勘探法在怀宁县月山镇探测岩溶塌陷区,表明该方法对地下岩溶甚至充填型溶洞有很好的反映。针对反演成果图中异常边界模糊、异常划定受人为影响大、存在误判、漏判等问题,Ward等[8]对三维反演电阻率数据采用模糊识别算法进行分类,实现自动划分基岩面。Chambers等[9]对采石场动态监测地下水位变化的图像进行模糊聚类处理,自动圈定水位变化界面。王国林[10]采用分布式模糊聚类算法对煤矿底板电阻层析成像(electrical resistance tomography, ERT)监测数据进行聚类分析,实现自动圈定异常范围。但是目前模糊聚类算法在微动勘探识别异常的研究还很少见,对提高微动异常识别的准确度值得研究。

现以皖赣铁路宁国改线K115+000～K115+500区间既有铁路路基岩溶勘察为例,开展分布式模糊聚类微动勘探,反演得到横波图像,并基于分布式模糊聚类法进行地层分层、溶洞自动划分,消除人为因素引起的异常判定,准确地探测到路基溶洞位置。

1 微动勘探法关键技术

20世纪50年代,Aki[11]、Capon[12]最早利用空间自相关法和频率-波数滤波法在微动信号中成功提取面波并求得频散曲线,从此微动走向应用阶段。早期的微动主要用于研究深部大地构造,近10年来随着中国基础设施建设如火如荼,微动勘探法已在浅地表工程探测中飞速发展。

微动勘探法是一种利用S波(横波)速度差异进行岩土分层及溶(土)洞分析的一种新兴地球物理勘探方法,以平稳随机过程理论为依据,从天然微动信号中提取瑞雷波频散曲线,再通过对频散曲线的反演获得地下不同介质结构的S波速度值[13-14]。故频散曲线的提取和横波深度反演是微动勘探技术的两项关键技术。

频散曲线提取主要包括SPAC法(spatial auto-correlation method)和F-K法(frequency wave number method)。SPAC法原理如下:设a(0, 0)、b(r,θ)两点为地表微动观测信号,两点在时间t内的空间自相关函数[15]为

(1)

(2)

式(1)中:irk为贝塞尔函数的展开形式,r为两点距离,k为波数;h(ω,φ)为频率方位密度函数;ω为角频率;θ为波的入射角;φ为观测角。

对式(2)取方位平均,即

(3)

定义a、b两点间的空间自相关系数ρ(ω,r,θ)为

(4)

F-K法原理如下,首先利用最大似然法求出不同频率的功率谱,由于采集的原始信号是由体波和面波组成的,如果某一信号优势明显,就会在功率谱上存在对应的最大值。根据这个最大值可以求出相应的传播速度和传播方向,故求出的不同相速度就可得到实测相速度的频散曲线[16]。由于本次探测场地有限,难以形成圆形观测系统,故本次研究使用F-K法提取频散曲线。

微动勘探应用的另一项关键技术即横波速度反演,目前相关方法主要有最小二乘法、遗传算法等人工智能算法。这些方法各有优缺点,其中最小二乘法比较依赖初始速度模型,而遗传算法可以在较大范围内进行全局搜索,对初始模型依赖较小,但精度相对较低。由于本项目前期对探测区域地质条件有一定了解,故采用最小二乘法反演横波速度。针对反演成果中土、岩体分界面模糊不清,异常位置及边界不准确等问题,本文中采用分布式模糊聚类算法分析以上反演成果数据。

2 分布式模糊聚类算法

利用分布式模糊聚类算法分析反演横波数据,根据概率密度曲线峰值确定质心位置和聚类数,实现物探异常自动判定,从而解决反演成果存在界面模糊、异常范围受人为影响、误判漏判等问题。

2.1 模糊聚类C均值

模糊聚类分析是以数据对象之间的相似性为基础,对数据对象进行分类的聚类方法。本文采用模糊C均值聚类算法,根据反演横波值特征进行分析。该算法的核心思路是将数据分为若干类,类与类之间数据差异较大,同类数据之间相似度较高,利用模糊隶属度值表征分类相关性,以欧式距离定义非相似度,从而构建目标函数,通过反复迭代数据归类分析,直至目标函数最小时达到最优隶属度矩阵和聚类中心,迭代停止[17]。

设定反演横波值为集合s={s1,s2,…sj,…,sN},sj为第j个数据,共有N个数据。建立集合S的隶属度矩阵U=C×[uij],C表示集合S被分为c类,uij表示sj隶属于第i个分类的隶属度,需满足条件:

(2)0

在上述准则基础上随机生成模糊矩阵U(0),进行初始化迭代。隶属度uij求解公式为

(5)

式(5)中:m∈[1,∞)是模糊指数或加权指数;dik=xk-vi表示数据sk与聚类质心vi的欧式距离;聚类质心集合V={v1,v2,…,vc},其中vi表示第i类聚类质心,vi的计算公式为

(6)

2.2 分布式模糊聚类

模糊聚类算法的初始矩阵具有随机性,存在无法保证每次聚类成果均达到全局最优解,并且存在反演数据量大,迭代耗时等问题。故提出分布式模糊聚类算法。该算法首先将全部数据进行核密度估计算出最优全局聚类质心,依据最优质心进行模糊聚类C均值分类,进而克服原算法对初始值的依赖和易陷入局部最优解,并跳出迭代寻优过程,减少计算时间。

核密度估计是将一个对称加权函数求和(简称“核函数”)按概率密度的差异应用到每一个数据中,而不是简单地将每点分配到每个间隔中,实质是近似地创建一个表示数据分布的直方图[18]。在核密度估计中,集合s={s1,s2,…sj,…,sN}的拟估计概率密度函数形式一般为

(7)

根据经验,核函数的选取对估计结果影响较小,带宽选择的成败直接影响核密度估计结果[19]。因此建立寻找最优带宽模型显得尤为重要,然而实际数据的核密度函数f是未知的,直接建立最优带宽模型较困难,本文以正态分布准则为未知函数f的参考分布函数建立最优带宽模型,进而求出最优带宽σ[20]为

(8)

采用Gauss函数作为核函数,概率密度函数为

(9)

3 数据采集与处理

3.1 地质情况

皖赣铁路宁国改线工程线路全长39.034 km,设港口镇站、罗田站及宁国南站3个车站。改建后,铁路线将外绕至宁国市区西侧,城市空间规划将更加合理,大大提高了区域间交通便利度,对宁国市市域集中、经济增长、社会发展、改善城市面貌具有重要而深远的意义。戴家塔大桥控制性工程测区位于新华夏系第二隆起带西翼区域的既有铁路线,该线分布有宁国组钙质页岩、西阳山组钙质页岩夹灰岩、新岭组页岩夹砂岩等多种岩性,地质构造复杂,经过多次构造运动,演变历史悠久,褶皱断裂发育,岩浆活动频繁,线路工程地质图见图1。其中K114+500～K115+035区段为探测区段,属寒武系西阳山组钙质页岩夹灰岩地层。场地范围内岩层面倾角30°～60°,平均43°,层面倾角较大,有利于溶蚀作用发育;区域内局部地层受隐伏断层构造作用,岩性破碎,岩层节理裂隙发育,地下水沿岩层面向下侵蚀,局部岩溶发育,见图2。

图1 线路工程地质图Fig.1 Line engineering geological map

图2 K114+500～K115+035段右侧山体(寒武系炭质页岩夹灰岩、灰岩溶蚀)Fig.2 K114+500～K115+035 section right mountain (Cambrian carbonaceous shale interbedded with limestone, limestone dissolved)

既有路基K114+500～K115+035段地层参考邻近新建铁路地层,该地层表层为第四系覆薄层种植土、粗圆砾土,下伏基岩为寒武系西阳山群(∈3x)钙质页岩夹灰岩组成,分为全风化、强风化和弱风化三个亚层。灰岩地层中遇溶洞,溶洞局部半填充状态,填充物为碎石土;局部溶洞无填充物,为空洞。桩基施工中可见地表塌陷,如图3所示。

图3 新建25#墩钻孔附近地面塌陷Fig.3 Ground subsidence near the borehole of newly built 25# pier

图4 微动法测点布置平面图Fig.4 Plan layout of measuring points for microtremor survey method

3.2 数据采集

新建铁路特大桥中心线距离既有路基K114+500～K115+035 段中心线仅22～35 m,新建25#墩钻探时已发生地面塌陷事故,既有路基受桩基施工扰动威胁显著,故推测 K114+500～K115+035 段既有路基塌陷风险较高,列车安全运营风险显著升高。查明该区段路基下方岩溶发育情况尤为重要,用以分析评估桩基施工对铁路临近路基稳定性的影响,确保列车安全运行。

采用微动勘探法详细查明既有路基段岩溶分布平面及深度范围。首先根据实际情况设置观测系统,如图5所示,本次工作测线位于现有铁路路肩上,周边较为狭窄,因此测量工作采用直线型布置方式进行微动勘探,共布2条测线,分别于路肩左右线,总长度1 070 m,点间距2 m,测点数536个,测点平面布置图见4。数据采集使用合肥国为电子有限公司的GN309微动探测系统,连接主频为2 Hz的垂直分量宽频带拾震器进行数据采集,每个测点采集至少30 min,采用4G远程通信对现场数据实时质量控制。

图5 野外物探测线布置图Fig.5 Layout of field geophysical detection line

3.3 数据处理

对采集的数据进行分布式模糊聚类反演分析,数据处理流程如图6所示。图7(a)为某测点原始微动信号,经过数据预处理得到S波速度随深度变化曲线如图7(b)所示,然后采用F-K法得到相速度频散曲线如图7(c)所示,最后对频散曲线进行反演获得该点S波速度结构。图8(a)为既有路基K114+500～K115+035左线微动勘探反演成果图,根据横波波速Vs大小,分为上中下三层,上层Vs呈低速,解译该层主要为第四系粗圆砾土;中层Vs呈中低值,解译该层主要为全风化～强风化钙质页岩夹灰岩;下层Vs相对较高,解译该层主要为弱风化钙质页岩夹灰岩。通过对反演横波数据进行分布式模糊聚类分析,分析得出最优带宽,聚类类型为4类,如图9所示。根据分布式模糊聚类分析成果,绘制微动勘探分布式模糊聚类图,如图8(b)所示。

图6 分布式模糊聚类微动法数据处理流程图Fig.6 Data processing flow chart of distributed fuzzy clustering of microtremor survey method

图8 K114+500～K115+035综合解译图Fig.8 Comprehensive interpretation chart of K114+500～K115+035

图9 横波速度概率密度曲线Fig.9 Probability density curve of s S-wave velocity

4 结果分析

微动法经过反演得出地下各层介质横波速度,依据《建筑抗震设计规范》(GB 50011—2010)[21]表4.1.3岩土的类型划分原则,可将反演横波速度大致将测区地层划分4层,分别为粗圆砾石层(0～400 m/s)、钙质页岩夹灰岩(W4～W3:400～600 m/s)、钙质页岩夹灰岩(W2:600～2 200 m/s)、溶洞(600～1 000 m/s),地层分层详见图8(a);采用分布式模糊聚类分析方法对反演横波速度值进行分析,地层分层情况为粗圆砾石层(0～298.2 m/s)、钙质页岩夹灰岩(W4～W3:298.2～594.9 m/s)、钙质页岩夹灰岩(W2:1 030.6～2 200 m/s)、溶洞(594.9～1 030.6 m/s),地层分层详见图8(b)。根据两种方法地层分层情况与钻孔验证情况分析对比,如表1所示。

表1 钻孔验证情况下微动反演与分布式模糊聚类分析地层分层情况的对比

从表1得出,反演横波顶板深度平均误差为22.7%,分布式模糊聚类分析顶板深度平均误差为8.9%,并且从钻孔21BZ戴家塔L03揭露地层可知,反演横波的岩体完整性深度分层误差偏大,明显与钻探结果不吻合,分析主要是由于全风化钙质页岩夹灰岩层深度厚,该层与粗圆砾土层横波速度差异较小,难以肉眼辨别区分。分布式模糊聚类分析的岩土层界面、岩体完整性分界面明显好于反演横波速度界面,且与钻探揭露深度吻合度较高。

根据微动反演横波速度大小划定溶洞异常,图8中分别划定了8个溶洞异常。其中两种方法判定的物探异常范围均较大的编号是:A-1和B-1、A-4和B-4、A-6和B-6、A-7和B-7、A-8和B-8。根据上述异常位置,布置5个钻探验证孔,分别为21BZ戴家塔L01、21BZ戴家塔L03、21BZ戴家塔L04、21BZ戴家塔L05、21BZ戴家塔L06,分别验证上述岩溶异常情况。

从图8可以看出,微动反演法和分布式模糊聚类法划定异常范围处钻探均揭露有岩溶,岩溶大小不一、以全充填-半充填为主,局部含空洞,存在串珠状溶洞,局部掉钻;两种方法划定的物探异常范围总体大于实际钻探揭露溶洞范围,如K114+839处反演横波物探异常A-6为9.3 m,分布式模糊聚类分析物探异常B-6为6 m,钻探揭露异常1.4 m;K114+976处反演横波物探异常A-8为12.0 m,分布式模糊聚类分析物探异常B-8为5 m,钻探揭露异常3.1 m。

从表2得出,分布式模糊聚类分析圈定的异常深度位置与反演横波圈定的异常深度位置相比,前者与钻探揭露深度吻合度更高,甚至可圈定局部微小异常。K114+762 处反演横波异常A-4范围-15.0～-33.0 m,分布式模糊聚类分析物探异常B-4范围-21.5～-41.5 m,钻孔21BZ 戴家塔L03揭露溶洞范围-24.1～-42.8 m,分布式模糊聚类分析圈定的异常层底深度与实际深度仅相差1.3 m,另外该位置分布式模糊聚类分析多判定异常B-5,钻孔也揭露该异常,说明该方法可更加准确判定微小异常。

表2 物探异常与钻探揭露情况对比

验证孔21BZ戴家塔L01、21BZ 戴家塔L03、21BZ 戴家塔L05中均存在串珠状溶洞,从图8所示,两种异常判定方法均难以区分串珠状溶洞。

5 结论

本文系统回顾了微动勘探法和分布式模糊聚类算法基本原理,以皖赣铁路宁国改线某区间路基岩溶勘察为例,开展分布式模糊聚类微动勘探进行地层分层、溶洞自动划分。得出如下主要结论。

(1)针对环境复杂与敏感地段(电磁干扰严重、接地不良、场地受限等因素)可采用微动勘探法进行岩溶探测,微动法可以较为准确的圈定岩溶发育区。

(2)本研究地层分界反演横波顶板深度平均误差为22.7%,分布式模糊聚类分析顶板深度平均误差为仅为8.9%。说明分布式模糊聚类算法可显著提高岩土、岩体完整性分界面划定的准确率,大大降低人为判定误差。

(3)利用反演横波速度和分布式模糊聚类分析方法划定的物探异常范围总体大于实际钻探揭露溶洞范围;分布式模糊聚类分析圈定的异常深度位置比反演横波圈定的异常深度位置更加精确,与钻探揭露深度吻合度更高,可圈定局部微小异常。

(4)地层中当存在串珠状溶洞时,利用反演横波速度和分布式模糊聚类分析法均难以区分串珠状溶洞,后续可考虑加强微动反演算法方面研究改进成像结果。