复杂地质构造空间分布评价的多源信息融合方法及应用

2020-05-31李培锋李斯涛袁从华

水利与建筑工程学报 2020年2期

李培锋，李斯涛，李春，严松，袁从华

(1.云南玉临高速公路建设有限责任公司，云南临沧 677000；2.中国科学院武汉岩土力学研究所岩土力学与工程国家重点实验室，湖北武汉 430071；3.中国科学院大学，北京 100000)

近年来，计算机技术迅速发展，新信息技术在岩土工程中逐渐得到应用，图像技术、数据库及GIS等新技术的研发推动着岩土工程的数字化及信息化进程。岩土工程与地理、空间位置信息紧密相关，并且地质信息通常具有多变性、复杂性，而信息技术能将采集到的信息进行关联及整合，因此在工程地质调查领域的应用前景十分广阔[1]。在岩土工程地质勘察中，测绘、踏勘、钻探及地球物理勘探等方法提供了不同尺度、不同层面的信息[2-6]，综合全面地处理这些信息是复杂地质构造勘查的关键，因此研究信息融合方法有利于获取更准确地地质构造信息。

在地质勘查领域，数据分析及图像融合等信息融合方法都已得到应用。2006年吴庆鸣等[7]提出了应用于超前地质预报的信息融合系统，可对采集到的数据进行分析从而得出前方掌子面的围岩类别。在遥测领域，图像融合方法的应用已十分广泛，由于地质调查最后也需得到图像信息，因此可以借此提高勘查结果图像的精度。2007年，赵珍梅等[8]将多光谱图像与雷达图像进行融合，对比分析后发现的图像融合后色域信息完整，同时兼有雷达图构造信息丰富的特性。2008年，印兴耀等[9]分析了传统的线型变换主成分分析法的劣势，提出了基于核主成分的PCA方法，并应用于对弹性波属性的降维处理中，提高了该方法的提取特征能力。2015年郭彪等[10]将小波变换应用于地震波多特性信息融合，对数据进行小波分解及重构，充分利用了各属性的主要特征信息，从而为地质勘查提供准确可信的资料。2017年，余洋等[11]针对探地雷达法，研究了信息融合算法，将不同频率探地雷达的信息进行融合，综合分析了不同频率探地雷达采集的信息，提高了结果的准确性并提升了图像的分辨率。

本文将信息融合技术应用于复杂地质构造勘探中，提出了应用于复杂地质勘探的数据相关性分析及图像融合方法，以某采空区为对象，采用三维激光扫描技术及高密度电法等手段对地下采空区进行综合勘测，并将数据结果进行融合，评价采空区的空间分布状态。

1 工程勘探中的信息融合方法

信息融合作为一种数据综合和处理技术，实际上是许多学科和新技术的集成和应用[12]，其基本原理是充分利用多物理量的信息，把多渠道、多尺度所获得的不完整的信息加以综合，消除多物理量数据间在空间大小及位置上可能存在的矛盾或多余，将各自的优点互补，降低不确定性，以形成对勘探目标的相对完整的一致性描述[13]。

本文针对采空区地下洞室调查的信息融合，提出了形态包络的概念，并采用了统计学分析方法进行数据相关性分析，采用累计标记法进行信息收敛从而实现图像融合。

(1) 形态包络。地下空间的地质体常常成群体分布，例如矿脉可长达几千米，厚度从不足一毫米以及到几米，形状可以比较简单成单脉出现，也可能成平行的、交错的脉群出现，矿脉包含很多分支，但总体上是沿着主矿脉延伸的。溶洞通常也是成群分布的，溶洞的发展经历了溶沟、孤立洞穴、洞室群、连通洞室到大型溶洞等一系列演化过程[14]。实际工程勘探的结果通常以剖面展示，在探测到成群分布的地质体时，剖面上会出现多个离散的异常体，因此需对其进行包络处理。

平面中包络线方程的求解如下所示：包络线与一族曲线系中的每一条都相切，如图1所示。假设某曲线族的每条曲线Cs可表示为t：x(s,t),y(s,t)其中s是曲线族的参数，t是特定曲线的参数。若包络线存在，则由s:[x(s,h(s)),y(s,h(s))] 得出，其中h(s)以下的方程求得[15]：

(1)

若曲线族以隐函数形式F(x,y,s)=0表示，其包络线的隐方程由下式消去参数s得出：

(2)

图1 形态包络

在地质勘探中常用到的高密度电法，其反演结果为等值线分布，因此结果与勘探目标实际形态不匹配，将同位置包络后的形态叠加到高密度电法剖面结果上，采用临界值原则，如图2所示，将临界值ρc的等值线包络的范围判定为勘探目标，并以等高线包络的形态作为推测的真实形态。

(2) 统计学分析方法。在进行数据分析时，通常需对多个信息源的数据进行相关性分析，其最基本的方法主要是统计学分析方法。常见的有数据分布直方图，回归分析，协方差矩阵及主成分分析等，本文采用了分布直方图分析高密度电法临界值的取值范围。

图2 高密度电法临界值

(3) 累计标记法。对于x×y×z的离散数据点阵，累计标记法能快速准确地确定三维数据点阵中的聚集簇，该方法沿三轴向依次对数据点进行扫描，在访问每一个数据点时，通过前邻点的标记状态确定当前数据点的标记，并根据数据点的标记状态进行依次分组，从而完成聚集簇中孤立数据点的识别。

对于数据点(x,y,z)，其与相邻6个数据点的位置关系如图3所示，根据扫描次序，将当前数据点的6个邻点，分为3个前邻和3个后邻，当前数据点和三个前邻点的值D(x,y,z)、D(x-1,y,z)、D(x,y-1,z)和D(x,y,z-1)分别记为Ds、D1、D2和D3。

图3 邻点示意图

采用依次沿x、y、z正方向的扫描方式对数据点阵进行顺序标记，若首先以D(x,y,z)=1的数据点为需要标记的目标，则D(x,y,z)=0 的数据点暂无需标记。综合所有聚集簇数据点整合情况，当前数据点为目标点时，当前元素的8种前邻状态及聚集簇整合情况列于表1中。其中m1、m2、m3分别为三个前邻点D1、D2和D3所在原聚集簇中数据点数目。

完成上述累计标记扫描后，以D(x,y,z)=0的数据点为需要标记的目标，再次进行扫描。所有数据点根据其前邻状态被分配至不同的聚集簇，将数据点数m小于一定值的聚集簇定义为孤立点簇，即完成孤立点簇过滤。

表1 8种前邻状态及聚集簇整合情况

2 工程背景及数据采集

本文采用GPS定位、三维激光扫描及高密度电法进行综合勘测，用形态包络及主成分分析的信息融合方法应用于某地下采空区的空间分布调查。

2.1 工程背景

项目区域基岩裂隙发育，水系也较发育，地下水丰富。河谷和盆地地区地形平坦，容易引起降水汇集，形成地表溪流、孔隙水及地下裂隙水。地表沿线水系较发育，河流呈树枝状，主要为山间溪流，呈西北-东南走向，河流形态多呈V字型，且河曲发育。地表水流量受降雨量影响，水位受季节影响大，在雨季水位暴涨，旱季地表水水位较低，工作期间应考虑到洪水的影响。降雨是该地区地下水供应的主要来源，地下水水位受降雨的影响大。

已发现一处较大采空区地下洞室，由于采空区洞室内状况、延伸与走向等不详，加之山坡上到处都有洞口露出，直接施工存在很大安全隐患，因此需要通过各种勘探手段把握内部状况、分布。测区位于华南活动带之东南沿海中生代火山断陷带中段，闽东火山断拗带之北段，中生代以来，区内岩浆活动频繁、强烈，其中火山活动尤为突出，形成了分布广泛、厚度巨大的一套陆相沉积-火山岩系，洞室内岩体状态如图4所示。

图4 测区地形地貌及洞内状态

洞室空间窄小，洞内蜿蜒曲折，可进入洞长约74.6 m，洞室高度约2 m～3 m，均布与山顶相通，宽度约40 cm～60 cm。与主洞室平行在洞室南部存在另外一条洞室。洞室周边火山岩体完整，节理裂隙不发育。洞内湿润，下部存在大量积水，水流不断，局部存在大面积积水。

2.2 现场数据采集

针对采空区这类空间形态复杂的构造，本文提出了一种信息融合技术。采用GPS定位，三维激光扫描及高密度电法进行勘测，用形态包络及统计学方法进行数据相关性分析，最终通过图像融合得出采空区的空间分布，其勘探流程如图5所示。

图5 采空区勘探流程

(1) GPS定位。采用GPS定位确定采空区所有露头洞口的位置信息A (x,y,z)，得出其在地形图上的分布，如图6所示，标注出所有露头洞口的位置，共12个，形态及尺寸如表2所示。另外在电法测线布设时，利用GPS对测线首尾及地形突变点定位，获取其大地坐标及高程，可将电法结果准确地与地形图信息结合。

图6 露头洞口位置

GPS设备型号为灵锐S82 GPS RTK，测量系统包括基准站和移动站，基准站由主机、数传电台、发射天线与电瓶组成。水平精度为±1 cm+1 ppm，垂直精度为±2 cm+1 ppm。

表2 露头洞口形态及尺寸

(2) 三维激光扫描。对已知的可进入区域可采用三维激光扫描技术进行精确的扫描，获取其详细的空间形态，为信息融合提供准确的信息，同时还可为之后本项目开发的设计和施工提供依据。

设备采用瑞士徕卡三维脉冲式激光扫描仪ScanStation C10，具有高达50 000点/秒的扫描速度，可以进行360°×270°全视场角扫描，获取目标顶部，水平方向以及垂直方向的数据。测量距离精度为2 mm，角度精度为12″。测量时首先在洞内初步勘探地形，根据需要确定扫描站数、中转基站的位置和控制标靶(用来匹配每站扫描的点)的个数和位置。然后安置三维激光扫描仪，调整好方向和水平。扫描仪预热后，设置好扫描参数(分辨率、曝光度和全景扫描等)，扫描仪自动进行扫描。

由于洞内曲折，需设置多个中转站，三维激光扫描在主洞共测量10站，左支洞共测量6站。大部分分支洞室尺寸较小，错综复杂，总体呈洞室群分布。

(3) 高密度电法。采空区往往存在顶板塌陷导致内部通路堵塞等情况，因此可能存在很多难以进入或未知的地下空洞。由于洞室沿矿脉分布并富含水，可期望利用洞室与周围岩体的电阻率差异推测出洞室分布范围，因此决定采用高密度电法进行勘探。根据前期高密度电法勘探的结果，表层电阻率在100 Ω·m以下，山坡的地表电阻率随含水程度在10 Ω·m～数100 Ω·m范围内变化。火山碎屑岩的电阻率随风化程度而变化，中—强风化火山碎屑岩电阻率约在1 000 Ω·m以下，弱风化—未风化火山碎屑岩的电阻率约数1 000 Ω·m；因此采用高密度电法对大范围内采空区进行勘测，总体上把握洞室的分布范围及走向。

高密度电法数据采集系统由主机、升压器、电极系三部分组成，野外工作一般由外接计算机、主机、电池组、主电缆等组成，如图7所示。主机包括发送控制命令、电流和电压测量、接收信号等部分，主电缆由多芯电线组成，主要作用是信号传输。

图7 高密度电法仪器

洞室主要分布在地表以下30 m以内，综合考虑勘探深度以及地表条件，高密度电法勘探采用2极法(Pole-Pole法)，远电位电极布设在测线2 km之外，每条测线布设60个电极，间距2 m，总长118 m，供电电压360 V。在采空区上方共布设54条测线，测线间距5 m，测线布设如图8所示，定义水平面上沿测线为x方向，垂直测线为y方向，竖直为z方向。

图8 高密度电法测线布设图

根据当地的地形、地质与地球物理条件，结合以往工作经验，勘探时采取了以下应对措施：在遇到松散土层或植被较发育地表，挖坑埋置电极，使电极与大地耦合良好，测区内大部分为山坡地貌，地形起伏大且坡度陡，为保证成果的解释精度，根据地形条件每隔15 m～30 m采用GPS放桩，并记录地形的突变的情况，以用作地形校正。高密度电法测试中采用360 V电压供电以加大电流密度，数据采集过程中加强监控，一旦发现仪器出现非正常反应或者数据出现非正常变化现象，立即停止观测，直到排除故障后再重新开始观测。通过采取以上若干措施，保证了物探现场采集数据的质量。同时前期详细记录了地形、地貌、地质变化，特别是断层露头、基岩出露等关键部位，为后期勘探结果的解释提供辅助和验证依据。

3 数据处理与融合

3.1 数据处理

(1) 三维激光扫描数据处理。三维激光扫描数据处理是通过Cyclone提供的坐标匹配功能，将各站的点云数据拼合成一个完整的测量目标点云模型。删除干扰的云点，通过分布框选点云，最终完成了对洞室的建模，并可对模型进行加工，生成剖面图，计算高程等，根据应用需求以AutoCAD和GIS的格式输出数据。原始点云包括了部分误差、错误和无关信息。这些误差信息的删减也是扫描数据处理的一个关键。影响三维激光扫描仪精度的要素较多，主要包括环境因素、激光信号的信噪比、激光信号的反射率、激光脉冲接受器的灵敏度、仪器和被测点间的距离、仪器和被测目标间的角度等。这些因素直接导致了误差的产生。通常情况下，三维激光扫描仪数据的模型精度显著高于单点的精度，可通过设立球靶、球状目标的坐标来改正点云的坐标，提高扫描仪模型的精度。

洞内宽窄不一，长52.7 m，最窄处约50 cm，最矮处约2 m，前段较水平，后半段坡度较大，有一长度约为40 m的上坡，洞内潮湿。将扫描数据整合后，可得到这一部分洞室的三维扫描图像，如图9(a)所示，将扫描图像的坐标进行降低精度的数据提取，得到了洞室的三维坐标，如图9(b)所示。

图9 三维激光扫描结果

(2) 高密度电法数据处理。高密度电法数据分析采用美国GEOTOMO公司的商用反演分析软件Res2dinv进行处理，对原始数据进行格式转换、噪音剔除、平滑与插值等预处理后，通过格式转换把数据重排成该软件处理格式，然后进行反演。高密度电法结果剖面共54个，前四条测线剖面结果如图10所示。将54个切面整合到同一空间坐标中，如图11所示，得到电法三维效果图。

图10 高密度电法剖面图

图11 电法结果三维效果图

3.2 数据分析

本次信息融合数据有3个来源，分别是GPS定位、三维激光扫描及高密度电法。区域GPS定位采集到地表的地形点的高程信息以及12个洞口的坐标信息。三维激光扫描对可进入区域洞室的轮廓进行高精度的扫描，得到洞室内壁上大量点的坐标，数据量庞大，可根据需要对其进一步处理。由GPS定位的露头洞口位置信息及三维激光扫描定义出的可进入区域洞室位置信息的可靠性及精度较高，可作为采空区勘探的边界条件。高密度电法能获取大范围内地下空间岩土体电阻率，其结果展示了观测剖面的地下电阻率分布，下一步需要根据电阻率信息解释出洞室在剖面上分布的位置信息。在观测区域总共布设了54条测线，因此共得出电阻率剖面图54张。信息评估如图12及表3所示。

表3 信息评估

图12 信息评估示意图

三维激光扫描能对仪器可进入区域进行精准的扫描，获取这部分区域地下洞室的空间分布与形态，数据形式以洞室轮廓上所有点的空间坐标展现。截取结果中的竖直方向截面，可以发现截面上包含多个局部的小洞室，但这些洞室呈群分布。采用形态包络的原则，描绘出洞室群的包络线，以此作为洞室群的位置，如图13所示。

图13 洞室群形态包络

电法勘探区域覆盖了可进入区域，可将此范围内电法剖面结果与三维激光扫描结果进行对比分析，取同一位置处的剖面，将二者结果叠加，如图14所示。由图可见，可进入的洞室与高密度电法结果在主洞室的位置对应性良好，洞室位置对应于电法结果中的高阻区域，因此可以将三维激光扫描结果融合到电法剖面结果中。对图中二者结果重合的部分，取此范围内电阻率最小值ρc作为临界值，将电阻率大于ρc的区域推定为洞室。

采空区的可进入区域上共布设了电法测线20条，将每个剖面上的电法结果与激光扫描的剖面结果进行对照，取得所有剖面上的电阻率临界值，绘制电阻率临界值分布直方图，如图15所示。可以发现临界值处于12 000 ohm·m～17 000 ohm·m，期望值约为14 500 ohm·m，因此以14 500 ohm·m作为临界值对所有电法剖面结果进行解释，即推定剖面上电阻率ρ>ρc的区域为洞室。

图14 剖面洞室分布

图15 电阻率临界值分布

3.3 图像融合

根据数据分析得到的临界值，可以对高密度电法的结果进行解释，得到电法结果解释空间图，如图16所示。可以发现解释后的结果信息繁多，且由于电法测量时存在一定误差，导致结果中包含离散的孤立点。采用累计标记法对孤立点进行过滤剔除，得到洞室的空间展布，如图17所示。

可以考虑以洞室区域的中心点确定其位置，并忽略其形态大小。结合地形图与GPS定位，可得出各电法剖面中洞室中心点在地下空间中的坐标信息C(x,y,z)。最终在地下空间中，根据矿脉的连续性，结合洞口的定位信息及地下剖面的洞室中心点的坐标信息进行连接，就可得出整个采空区地下洞室的分布图，如图18所示。

根据结果可发现勘探区域内除可进入区域洞室外，还存在多条主要贯通性洞室，与可进入洞室相连，区域北侧形成错综复杂的多条相连通道。项目结束后在区域北侧的路面开口和洞口进行了详细踏勘，现场情况与信息融合后的结果进行比对验证，发现二者在空间位置上基本吻合。另外，此结果表明南北两端洞室是不连通的，在中间的不连通区域进行了多次钻探，同样未能探测出空洞。