基于三维网格模型重构的岩溶地质构造稳定性数值分析

王启明 车爱兰 周晗旭

摘要 由于岩溶地质的离散性大，传统建模方法并不能真实反映其空间分布状态，数值分析结果的可靠性受到限制。为精确分析无锡市某超高层项目基础落位区岩溶地质构造稳定性，本文基于三维网格模型重构技术，建立该区域岩溶地质构造的三维有限元分析模型，分别分析了溶洞处治前及注浆加固后的稳定性，结果表明：溶洞未处治时，溶洞周边岩体应力较大且存在应力集中，不均匀沉降现象明显。对溶洞进行注浆加固处治后，溶洞周边岩体应力集中现象减弱，模型顶面最大沉降量下降了54.55%，不均匀沉降趋势大大减小。通过上述分析，不仅对该超高层项目基础的设计和溶洞的处治具有重要指导意义，也为类似工程的数值分析提供参考。

关 键 词 岩溶地质;三维网格模型重构;注浆加固;数值分析

中图分类号 TU47     文献标志码 A

Abstract Because of the large discreteness of unfavorable geology， the traditional modeling method cannot truly reflect its spatial distribution， the reliability of the numerical analysis is limited. In order to accurately analyze the stability of karst geological structure in the foundation area of a super-high-rise project in Wuxi， a 3D finite analysis model of karst geological structure is established based on 3D mesh model reconstruction technology， and the stability of karst cave is analyzed before and after grouting reinforcement. The results show that when the cave is not treated， the rock mass around the cave has large stress concentration and obvious non-uniform settlement. After the grouting reinforcement， the stress concentration of rock mass around karst cave is weakened， the maximum settlement on top of the model is reduced by 54.55%， and the non-uniform settlement trend is greatly reduced. The above analysis not only has an important guiding significance for the foundation design of the super-high-rise project and the treatment of karst caves， but also provides a reference for the numerical analysis of similar projects.

Key words karst caves; 3D mesh model reconstruction; grouting reinforcement; numerical analysis

0 引言

我國幅员辽阔、地质复杂，岩溶地区分布广泛，大量在建大型工程位于岩溶发育地区，在岩溶地区进行基础工程设计与施工时，下伏溶洞对基础工程稳定性有至关重要的影响[1]。溶洞发育程度、溶洞跨度、顶板厚度、围岩特性等与溶洞上方的地基承载力密切相关[2]。地基滑动和坍塌、承载力不足、不均匀沉降等是溶洞存在常常引起的地基变形问题，影响基础的稳定性[3]。因此，施工前必须对岩溶地质的情况做好仔细的研究和认真的准备，充分考虑其对施工造成的各种不利影响以及应当采取的主要技术措施，以确保施工过程顺利、安全。否则，将会埋下严重的安全隐患，极有可能造成工期延误、重大经济损失和人员伤亡[4-5]。

对于岩溶地质构造稳定性进行准确的分析及评价，为工程问题的及时预测和适时决策提供有力支持，已经成为岩土工程界热点研究课题，大量学者从室内试验、数值分析、理论研究等多方面对岩溶地质展开研究，取得了丰富的成果。张永杰等[6]以湖南某高速公路线岩溶地质发育段为基础，分别从定量和定性评价指标对岩溶区公路路基稳定性进行二级模糊综合评判，使公路路基稳定性评价更具合理性。张聪、阳军生[7]等对新研制的环保可控膏状注浆材料（TGM）进行室内试验，测试了TGM材料的基本力学性能，并应用于工程实际，取得了良好的效果。近年来，数值分析方法在地下岩溶研究中得到了越来越广泛的应用，潘东东等[8]通过建立岩溶隧道承压隐伏溶洞突水模型进行数值分析，研究了地下溶洞突水对隧道开挖的影响，赵明阶[9]针对层间空隙岩溶建立相应数值模拟模型，分析了围岩在开挖过程中的应力、位移和涌水情况。上述数值分析方法根据单一岩溶形态建立数值模型，研究特定岩溶地质状态的数学和力学特性，取得了良好的效果。而实际地下岩溶发育极不均匀且复杂多变，因此根据某种特定岩溶状态建立模型进行数值分析，并不能如实反映地下岩溶地质的空间分布状态，计算结果的可靠性受到限制，故难以在工程中应用。

因此，为精确分析无锡市某超高层项目基础落位区的岩溶地质构造稳定性，对该区域进行弹性波CT探测，利用三维网格模型重构技术建立该区域岩溶地质构造三维计算模型，通过数值模拟分析岩溶地质构造的应力和沉降值，对该项目基础的设计与溶洞处治有重要指导意义。

1 弹性波CT探测现场试验

1.1 弹性波CT工作原理

地球物理勘探最先应用于油气资源开发，随着物理学、数学、电子技术、计算机技术的发展，在不良地质探测领域取得明显的研究进展，比较成熟的运用于不良地质探测方面的方法有高密度电阻率法、探地雷达法、弹性波CT法、电磁波CT法等几种，其中弹性波CT法作为一种领先的地球物理勘探技术，具有探测效率高、探测精度高、抗干扰能力强等优势，在岩溶地质探测中取得良好的应用效果。

利用弹性波对被探测对象一定厚度的层面进行扫描，在某一钻孔中的激发点激发时，分别在其他钻孔不同接收点进行接收，扫描方式有水平同步扫描、斜同步扫描、定点发射扇形接收扫描。图1为弹性波CT射线网格示意图。

1.2 工程概况

无锡市拟建某超高层建筑，初步勘察资料显示，该地区地层岩性中上部为第四系全新统、更新统沉积物，主要成分为粘性土及少量粉土，下部基岩为三叠系下青龙（T1x）灰岩，灰-肉红色。现场钻孔取芯，芯样如图2所示，在基岩面以下，标高-78.00～-88.00 m范围内，发育有空溶洞和半填充溶洞。根据《岩土工程勘察规范》对该项目场地进行岩溶地质弹性波CT探测试验。其中，在基础落位区域共布置钻孔25个（ZK1～ZK25），深度约90 m，钻孔孔径不小于70 mm，孔口至孔底均有PVC管以保证钻孔畅通。钻孔布置如图3所示。

1.3 数据采集

于2018年1月25日至2月10日，对基础落位区域进行彈性波CT数据采集。弹性波CT采集设备主要有数据记录系统、井下水声检波器线阵、大功率电火花震源（如图4），各部分工作原理如下。

数据记录系统：由高精度宽频带数据采集模块（Goed数字地震仪）和笔记本电脑构成，负责对换能器线阵接受到的信号进行模数转换和数字记录。

井下水声检波器线阵：由12组换能器（传感器）按0.5 m间隔排列构成12道检波器线状阵列用于在井孔中接受弹性波。水声检波器由宽频带高灵敏度压电陶瓷元器件和放大电路构成，把作用在器件上的压力波转换成电信号。

大功率电火花震源：由电极高压放电形成电弧，电弧瞬间将电极周围的水体加热成蒸汽，蒸汽膨胀并爆炸在地层中激发出弹性波。

具体采集参数设定如下：接收道数为12道，道间距为0.5 m，激发点间距为0.5 m，每组激发点数为12个，完成一组激发后，12道接收道上移3 m。采样间隔20.833 μs，持续时间0.15 s。在25个钻孔之间，采集72组弹性波CT剖面数据。采集到的波形例图5所示。

1.3 数据处理及结果

由于数据采集现场含有大量的施工噪音、杂波、多次波，而且地层介质对弹性波也有不同程度的吸收，所以需要对数据采取一定的处理，最大可能消除干扰波，使波速构造图像能最有效地展现地下特征。具体处理工作如下。

1）数据预处理。现场记录数据中并不包含检波器和激发点（震源）的位置信息。这些信息一般记录在操作记录中，为了便于数据的一体化管理，需要把这些数据添加到数据文件的文件名中，使数据和采集信息合为一个整体。

2）波形处理。首先进行滤波、视窗和叠加处理抑制噪音，提高信噪比。除此之外，利用数字滤波技术再一次抑制噪音，提高数据质量，以至于首波更为明显。

3）首波走时提取。运用自行编制软件自动读取每个波形的首波运行时间，并自动进行对应于位置信息的记录。

4）波速计算与岩溶分布评价。获取首波走时数据后，采用射线追踪算法、最佳路径算法等方法，通过计算机处理，建立井间介质的二维图像，依据特定插值方法将一系列上述二维图像转化为三维空间分布图像，将地下溶洞的分析结果绘制在三维地层空间内。

以剖面ZK6-ZK19为例，对现场测试结果进行地质解释。测试范围标高-77～-88 m，剖面宽度9 m。剖面内测点间平均纵波速度[Vpa=2 760]m/s，最小波速[Vmin=16 266]m/s，最大波速[Vmax=3 399]m/s。剖面岩层中存在低速区，岩溶较为发育，多为黏土或碎石充填。图6为剖面CT测试分析成果图。

2 岩溶地质三维网格模型重构

2.1 建立三维可视化模型

根据弹性波CT探测结果，对图2虚线标识范围内6个钻孔间探测区域进行三维网格模型重构，首先按照厘定的勘测线依次对弹性波CT切片进行三维空间展布，如图7所示。

在弹性波CT法探测中，CT切片相交将探测区域划分成若干个五面体，若直接进行Kriging插值，每个五面体内离散数据点的插值对象不仅包括五面体的3个侧面，还包括其他弹性波CT切片，模型的精确度大大降低。在插值过程中，以弹性波CT切片及其交线为边界进行区块化Kriging插值，获得三维可视化模型。

弹性波CT切片在三维规范化空间中的区域范围是8.00 m×8.00 m×11.00 m的六面体，根据目标区域大小及所要重构模型的网格单元规模，将三轴向上离散数据点的间隔值设置为0.25 m×0.25 m×0.25 m，则生成37×73×45的规则离散数据点阵。对离散数据点阵进行区块化Kriging插值，所有离散数据点统一进行三维可视化成像，结果如图8a）所示，图8b）为沿z轴方向切片。

2.2 三维网格单元构建

理论上，三维网格构建就是以单类型数据作为建模数据来源构建三维有限元网格模型，并建立不同材料属性与网格间的映射关系。三维网格构建的步骤如下：1）材料阈值分割;2）网格划分及材料属性映射。

2.2.1 材料阈值分割

材料阈值分割的过程：1）确定阈值Te;2）三维可视化模型[F（i，j，k）]中的离散数据点，其值对应不同的值域段;3）区分不同材料的分布区域[10]。可用如下公式表示，即

所以，对于材料阈值分割，最关键的一步是选择合适的阈值T。峰值法[11]是一种全局阈值分割方法，假设不同材料对应不同的值域段，用多个正态分布函数表示，相邻两函数的曲线拟合会得到2个分离的峰值，那么两峰之间的谷底对应的值就是分割阈值T。引入直方图h（F），该直方图通过统计不同值域段离散点的数目来确定，即

2.2.2 网格划分及材料属性映射

建立与数据点阵形成映射关系的三维网格模型，三维网格模型由8节点的六面体单元构成，每个数据点对应1个单元，数据点位于单元中心，单元沿三轴向的尺寸与数据点三轴向间距相同。对网格模型的节点编号，如图9所示，模型第n个节点的坐标为

通过以上步骤，完成了三维网格模型重构，所得三维网格模型如图10所示，网格步长设为0.25 m×0.25 m×0.25 m，共划分121 545个节点，114 048个单元。该三维网格构建方法可以控制三维网格模型的单元规模，在保证分析精度的前提下提高计算效率。

3 数值分析

3.1 模型材料参数及荷载

采用商用有限元分析软件Abaqus的有限元静力分析模块进行承载力数值模拟计算，研究岩溶地质构造稳定性。根据工程勘察报告，基岩为灰岩，且发育空溶洞和半填充溶洞，在数值分析中，将溶洞简化为饱和黏土材料，各材料具体参数如表1所示，表中：ρ为材料重度，E为弹性模量，μ为泊松比。

在数值分析中对图9所示模型的的位移边界进行设置：约束模型前、后、左、右4面x轴和y轴方向的位移，约束模型底面x轴、y轴和z轴方向的位移，模型顶面为自由边界。计算时，首先在模型中施加重力荷载，待初始应力平衡后，在模型顶面施加均布荷载，分别计算均布荷载为50 MPa、100 MPa、150 MPa、200 MPa这4种工况。

3.2 稳定性分析

圖11为均布荷载200 MPa工况下模型的应力等值云图，从图中可知，在模型顶面均布荷载作用下，在溶洞周边岩体应力较大且存在应力集中，有塑性破坏风险，从而使溶洞顶板坍塌进而向上延伸，引起地基承载能力不足，影响上部基础稳定性。

图12为均布荷载200 MPa工况下模型的沉降等值云图，从图中可知模型顶面产生不均匀沉降，最大沉降量达110.73 mm，最小沉降量64.39 mm。

3.3 注浆加固效果评价

从以上分析可以得知，上部基础稳定性受下伏溶洞影响很大，如不进行有效处治，存在安全隐患。综合考虑现场施工条件和经济因素，建议采用注浆填充对溶洞进行加固，提高基础稳定性。根据现有的注浆材料的室内试验报告，注浆填充材料的参数如表2所示。

通过数值模拟分析，计算求得注浆加固后岩溶地质构造的力学特性。图13为注浆加固后均布荷载200 MPa工况下模型的应力等值云图，对比图10可知，溶洞周边的应力集中明显减小。

图14为注浆加固后均布荷载200 MPa工况下模型的沉降等值云图，模型顶面最大沉降量由110.73 mm减小至50.33 mm，下降了54.55%，最小沉降量由64.39 mm减小至43.57 mm，不均匀沉降现象明显减弱。

沿模型顶面y轴方向中轴线，提取点A（4.5，4.5，-77）、B（4.5，9，-77）、C（4.5，13.5，-77）注浆前后的应力和沉降值，如图15所示。

在A、B、C 3点，注浆后不同荷载作用下应力和沉降量都有明显降低。分析其原因，当溶洞填充后，溶洞内部填充的注浆结石体对溶洞顶板有一定的支撑作用，不仅能有效减少不均匀沉降，而且还能使岩溶地质构造内应力分布更加均匀，减少应力集中，并能有效防止溶洞顶板破坏及进一步向上延伸破坏。

4 结论

1）本文提出了基于三维网格模型重构的岩溶地质构造稳定性数值分析方法，该方法根据弹性波CT探测结果，重构能精准反映岩溶地质分布状态的三维网格模型，提高了数值分析结果的可靠性。

2）该方法中，对插值对象进行分区块Kriging插值，避免了非相关切片的影响，提高建模精度的同时，实现岩溶地质弹性波速度构造三维可视化。

3）本文基于岩溶地质三维可视化模型数据结构的特点，设计三维网格构建方法，可控制三维网格模型的单元规模，在保证分析精度的前提下提高计算效率。

4）开展无锡市某超高层项目基础落位区弹性波CT探测现场试验，并重构该区域岩溶地质构造三维网格模型，通过数值分析获取了该区域不同荷载作用下的应力和沉降状况，对该项目基础的设计施工有重要指导意义。同时对溶洞注浆加固处置后的应力和沉降状况进行分析，结果表明：溶洞周边岩体应力集中现象减弱，溶洞顶板稳定性提高，模型顶面最大沉降量下降了54.55%，不均匀沉降趋势大大减小。

参考文献：

[1]    付敬，丁秀丽，张练. 岩溶回填改善地下洞室群围岩稳定性的数值分析[J]. 长江科学院院报，2006，23（4）：47-50.

[2]    唐前松. 岩溶区公路路基稳定性及处治效果三维数值分析[J]. 公路工程，2017，42（1）：143-146.

[3]    苏谦，黄俊杰，白皓，等. 铁路站场岩溶地基加固补强技术与效果分析[J]. 岩土力学，2013，34（3）：776-782.

[4]    李尧，李术才，刘斌，等. 钻孔雷达探测地下不良地质体的正演模拟及其复信号分析[J]. 岩土力学，2017，38（1）：300-308.

[5]    付代光，周黎明，肖国强，等. TSP预报隧道不良地质体应用研究[J]. 地球物理学进展，2016，31（1）：417-426.

[6]    张永杰，曹文贵，赵明华，等. 岩溶区公路路基稳定性的区间模糊评判分析方法[J]. 岩土工程学报，2011，33（1）：38-44.

[7]    张聪，阳军生，谢亦朋，等. 富水岩溶越江盾构隧道注浆材料试验与应用研究[J]. 岩石力学与工程学报，2018，37（9）：2120-2130.

[8]    潘东东，李术才，许振浩，等. 岩溶隧道承压隐伏溶洞突水模型试验与数值分析[J]. 岩土工程学报，2018，40（5）：828-836.

[9]    赵明阶，金鹏. 层间孔隙岩溶影响隧洞稳定性的数值模拟研究[J]. 地下空间与工程学报，2017，13（3）：827-832.

[10]  钟江城，周宏伟，任伟光，等. 基于CT图像灰度分布的含杂质煤体三值化方法[J]. 力学与实践，2018，40（2）：140-147.

[11]  彭启伟，罗旺，冯敏，等. 改进二维Otsu法和果蝇算法结合的图像分割方法[J]. 计算机应用，2017，37（S2）：193-197.