2017, Vol. 60
, ZHANG Xin-Xin, LIU Bin
, LIU Zheng-Yu, WANG Chuan-Wu, GUO Qian, LIU Hai-Dong, WANG Hou-Tong
隧道施工过程中经常遭遇突涌水等地质灾害，给隧道施工带来了严重的经济损失、人员伤亡以及环境破坏(Taherian，2015).在隧道施工期开展突涌水灾害源的超前探测，提前探明含水构造的位置具有重要意义.对于隧道突涌水灾害源超前探测而言，三维电阻率探测方法因其对含水体响应敏感，被越来越多的用于隧道含水体的超前探测中.隧道三维电阻率观测方式可以分为定点源三极测深类和聚焦探测类(李术才等，2014).定点源三极测深方法首先被引入到隧道超前探测中，对低阻异常体的电阻率响应特征进行了研究，采用反演方法实现了对含水构造的定位，并在实际工程中得到验证(程久龙等，2000；黄俊革等，2007；刘斌等，2012a).德国GD(geohydraulic data)公司基于聚焦电法的原理，研发了用于隧道超前探测的BEAM(bore-tunneling electrical ahead monitoring)技术，利用不同里程的连续探测结果来定性推断掘进面前方的含水情况(Kaus and Boening, 2008).李术才等(2015)基于聚焦与测深的思路，提出了一种多同性源阵列电阻率的超前探测方法，实现了对隧道掌子面前方含水构造的反演成像.
三维电阻率在隧道超前探测中的发展，对三维电阻率探测的快速化解译与成像提出了更高的要求，而三维电阻率反演是实现隧道含导水构造成像的关键.三维电阻率反演方法主要包括线性反演方法与非线性反演方法.线性反演方法具有局部搜索能力强，收敛速度快的优势，在电阻率反演中发挥了重要作用(吴小平和徐果明，2000；宛新林等，2005；刘斌等，2012b；Khalil and Santos, 2013), 但其对初始模型的依赖性较强，使反演结果容易陷入局部最优.非线性反演方法因其全局搜索能力强，不依赖初始模型，可较好的弥补线性反演方法的缺点，但往往收敛很慢，计算时间消耗大(徐海浪和吴小平，2006；Singh et al., 2010；Liu et al., 2012；戴前伟等，2014).
为了发挥线性反演与非线性反演的优势，克服两者的缺点，线性与非线性相结合的混合反演方法已成为一种新的发展趋势.张凌云和刘鸿福(2011)用神经网络优化的方法与蚁群算法实现二维电阻率的非线性联合反演，相比较单一反演方法，联合反演有效减少了迭代次数，提高了反演精度，获得更好的反演结果.李术才等(2015)提出了一种将光滑约束的最小二乘方法与蚁群算法相结合的混合反演算法，通过线性反演与非线性反演相互提供初始模型，实现了三维电阻率“嵌套式”反演.刘斌等(2017)提出了一种携带不等式约束的最小二乘法与改进遗传算法相结合的混合反演方法，实现地电结构的三维成像并成功应用于实际工程中，压制了反演的多解性、提高了反演效果.高明亮等(2016)利用免疫遗传算法对神经网络的反演参数进行同步优化，设计了一种将神经网络和免疫遗传算法进行有机结合的全局优化反演策略，提高了电阻率反演的精度.
以上研究成果为混合反演方法进行了有益的尝试与探索，为解决三维电阻率超前探测快速化解译与成像提供了可借鉴的思路.目前，混合反演主要集中在二维模型，三维电阻率混合反演相对于二维反演而言，增加了大量的网格数，计算速度问题一直未得到有效解决.同时，为了实现超前探测的较准确成像，必然在三维反演的基础上继续增加网格的数量，这将大大增加三维混合反演的计算耗时.混合反演的计算速度已经是制约混合反演发展的瓶颈问题，也成为超前探测快速化解译的必然要求.随着计算机技术的飞速发展，GPU(Graphics Processing Unit)并行计算的出现为解决混合反演计算速度的问题提供了一种解决方案.本文在笔者团队提出的三维混合反演基础上，改进蚁群算法与最小二乘线性反演的混合策略，提出了基于GPU并行的蚁群算法与最小二乘相结合的混合反演方法.利用蚁群算法的全局搜索能力强的优势为最小二乘线性反演提供较优的初始模型，克服最小二乘算法易陷入局部最优的缺点.基于GPU并行计算能力，增加反演模型的网格数量并提高蚁群算法中初始群体的数量，可以提高混合反演的搜索速度和反演效果.这种混合反演方法解决了单一反演方法反演效果差与成像速度慢的问题，提高了蚁群算法的搜索速度与全局搜索能力，确保三维反演收敛于全局最优解，改善反演结果的精度.同时提高了混合反演的搜索速度与计算速度，大大节约了计算耗时，可以实现隧道三维电阻率超前探测的快速化解译与成像.
2 基于GPU并行的蚁群算法三维电阻率非线性反演结合线性反演与非线性反演的优势，将蚁群算法与最小二乘线性反演进行联合，提出了一种三维电阻率超前探测的混合反演方法，可以压制三维反演的多解性并提高反演成像效果.然而由于快速化解译与反演成像的要求，需要增大计算模型的网格且提高混合反演的计算速度与搜索效率.混合反演耗时部分主要集中在蚁群算法部分，而蚁群算法具有天然的可并行性，基于CUDA(Compute Unified Device Architecture)运行平台实现蚁群算法的GPU加速算法，可以提高反演效率，减少三维电阻率反演的耗时.
2.1 基于蚁群算法的三维电阻率非线性反演的并行性分析蚁群算法的灵感来源于大自然界蚂蚁的觅食行为，蚂蚁在寻求食物的过程中，通过信息素的诱导作用，总是能够找到蚁巢与食物之间的最短路径(Colorni et al., 1992；Conti et al., 2013).因此将蚁群算法用于三维电阻率反演中，通过对蚂蚁群体选择的多个模型进行正演，并利用每个模型正演数据与实测数据的拟合差作为目标，引导蚂蚁在电阻率节点移动选择，最终得到全局最优路径(蚁群算法反演的全局最优解).蚁群算法三维电阻率反演目标函数如式(1)所示(张凌云，2011)：
|
(1) |
式中，Φmn是数据拟合差，N为蚂蚁个数，ρdi为第i个实测数据，ρsi为第i个正演数据.Φmn越小则表明蚂蚁选择的模型与实际模型越接近.当拟合差满足收敛条件时，此时蚂蚁选择的模型即为三维非线性反演的结果.假定蚂蚁群体数量为N，三维反演区域网格数量为Q，每个网格上可供蚂蚁选择的电阻率节点数为k个，蚁群算法的过程及并行性的分析如下：
(1) 蚂蚁群体(数量为N)根据选择概率选择三维反演区域网格的电阻率节点值，可得到N个模型(每个模型包含Q个电阻率值)，其中选择概率如式(2)所示(Dorigo，1992):
|
(2) |
式中，Pijk为蚂蚁k选择第i单元上第j个节点电阻率值的概率，τij为第i个网格单元上节点j的信息素，α为信息素强度的权重因子.
蚂蚁群体在选择模型时，每只蚂蚁的移动与群体中其他蚂蚁的移动相互独立，蚂蚁群体选择N个模型可以并行执行.同时对于单只蚂蚁而言，每个网格单元上电阻率节点的选择之间也是相互独立的，蚂蚁选择不同网格单元的电阻率节点也可以并行执行.
(2) 蚂蚁群体选择N个模型之后，需要对选择的N个模型进行有限元正演，分别求解N个模型的正演数据ρs，以求解每一个模型对应的数据拟合差Φmn.其中有限元正演如式(3)所示:
|
(3) |
式中K和K′分别为异常电位和正常电位向量的总体系数矩阵，d为理论观测数据，u0为正常电位.
由于N个蚂蚁选择的N个模型进行有限元正演过程中，N个模型是相互独立的，可以采用并行算法对N个模型进行正演计算，从而得到N个模型对应的正演数据ρs.
(3) 分别计算每个选择模型的观测数据与实测数据的拟合差，找出本次选择最小拟合差对应的模型，即本次选择中最优的解.然后按照选择的模型对相应网格上的节点的信息素强度进行更新，如式(4)(张凌云，2011)所示:
|
(4) |
式中，τijnew为信息素更新后的强度，τijold为当前信息素强度，λ为信息素挥发因子，Δτijk为蚂蚁k所增加的信息素的值.
在网格单元不同节点上信息素更新过程中，每个网格单元之间是独立的，且同一个网格单元不同电阻率节点也是独立的，因此可以采用并行计算的方法对信息素的强度进行更新.
(4) 网格上不同电阻率节点信息素更新后，蚂蚁在选择电阻率节点时的选择概率发生了变化，对应于目标函数最小的路径上的信息素被增大，其他路径上的信息素被减少，正是由于目标函数控制下信息素的引导与反馈作用，使得一定次数之后，蚂蚁群体可以找到一条全局最优路径，可以得到三维电阻率蚁群算法的最优反演结果.
2.2 基于蚁群算法的电阻率非线性反演并行策略与实现从2.1节的分析可知，蚁群算法具有天然的可并行性，利用并行算法可以加速蚁群算法的计算速度.GPU作为一种多线程的、高度并行的、具有强大浮点运算能力的众核处理器，拥有超多计算核心，远超传统的CPU多核计算(Delévacq et al., 2013；Angelo，2013).因此，可以利用GPU并行计算技术开展三维电阻率混合反演.
GPU以CUDA作为新的并行计算设备平台，其编程模型以CPU作为主机(host)，主要负责逻辑性强的串行计算；GPU作为设备(device)，主要负责执行高度并行的线程化任务.在GPU上运行的CUDA并行计算函数定义为内核函数.内核函数以线程格(Grid)的形式组织，每个线程格中包含若干个线程块(Block)，每个线程块中又包含若干线程(Thread).如图 1所示.因此，设计了三个CUDA的内核函数(AntsMoveKernel，Choleskykernel，Update Kernel)分别对蚂蚁选择模型、求解模型的正演观测数据以及信息素强度的更新进行并行加速，如图 2所示.具体方法如下：
|
图 1 CUDA的线程组织(陈召羲等，2012) Fig. 1 Thread organization on CUDA(Chen et al., 2012) |
|
图 2 基于GPU并行算法的蚁群算法反演流程 Fig. 2 Process of Improved ACO algorithm Based on GPU Parallel Algorithm inversion |
(1) 蚂蚁选择模型的并行化
设计了AntsMoveKernel核函数用来实现蚂蚁群体选择模型(蚂蚁选择路径)，如图 2中Grid1所示.根据CUDA线程组织模式，将每个蚂蚁个体转化成一个线程块(Block)，N只蚂蚁则有N个线程块，一个线程块完成一只蚂蚁的路径选择.每个线程块里分配若干线程(Thread)，每个线程负责每只蚂蚁的一步移动，根据蚂蚁所走路径中移动步数的总数决定线程块中分配线程的个数.N个线程块中线程同时并发，完成蚂蚁的路径选择，产生N组模型mN.
(2) 模型的正演观测数据的并行化
Choleskykernel核函数负责N个模型进行有限元正演获取观测数据，如图 2中Grid2所示.将N个选择的模型的正演计算映射到N个线程块上.而模型的三维有限元正演的耗时主要集中在总体系数矩阵K的Cholesky分解中(Rennich et al., 2014；李术才等，2016).根据CUDA线程组织模式，将每一个总体系数矩阵K映射到一个线程块上.而在矩阵K的Cholesky分解中(K=LTL，其中L是对角阵)，对于L中的每一列元素lij，在求得该列主对角线元素lii(i=1, 2，…，F)(F为矩阵K的行数)后，该列的其他元素在求解时相互之间没有影响，因此，在每个线程块中分配若干线程，用以计算L每一列中各元素的值.因此，利用Choleskykernel核函数可以并行执行蚂蚁选择的N个模型的正演计算，从而可以快速得到N个模型的正演观测数据.
(3) 信息素强度更新的并行化
设计了UpdateKernel核函数负责对蚂蚁信息素强度τij的更新，也就是对拟合差最小的蚂蚁选择节点的信息素强度进行更新.每个线程块中分配与网格的电阻率节点数目相同的线程.每一个线程负责计算网格的电阻率节点上蚂蚁留下的信息素.从而可以实现信息素强度更新的并行计算.
3 基于GPU并行加速的三维电阻率混合反演算法基于蚁群算法与最小二乘线性反演的混合反演方法，利用ACO的全局搜索能力为最小二乘提供较优的初始模型，以克服最小二乘算法容易陷入局部最优的缺点.笔者团队曾对蚁群算法与最小二乘混合反演进行过探讨，提出了最小二乘算法与蚁群算法相互提供初始模型的“嵌套式”反演方法(李术才等，2015).首先进行最小二乘线性反演为蚁群算法提供初始模型，用所提供的初始模型进行蚁群算法非线性反演，根据非线性反演的结果，重新设置线性反演初始模型进行线性反演，直到满足收敛条件.这种“嵌套式”反演方法在一定程度上提高了反演效果.然而，由于受计算速度的影响，反演时将掌子面前方的反演区域以层状网格剖分，大量减少了网格数量，影响反演成像的精度.因此，在文献(李术才等，2015)的基础上，引入了GPU并行算法对混合反演进行加速，实现掌子面前方三维网格剖分的混合反演成像.同时，将“嵌套式”反演方法进行优化，提出了一种基于GPU并行的蚁群算法与最小二乘混合反演策略.利用蚁群算法的全局搜索能力为最小二乘提供较优的初始模型，然后进行线性反演，理论上具有改善反演成像效果并提高计算速度的优势，如图 3所示，其具体流程如下：
|
图 3 三维电阻率混合反演流程图 Fig. 3 Process of 3D resistivity joint inversion |
(1) 建立3D反演模型，设定模型初始参数.
(2) 设置蚁群算法的初始信息素强度以及蚂蚁数目N等参数.
(3) 计算蚂蚁选择每个网格单元中每个电阻率节点的选择概率.
(4) 根据蚂蚁的选择概率，利用GPU并行策略，同时得到N只蚂蚁选择的N组模型参数.
(5) 利用GPU并行策略，对获取的N个模型进行基于Cholesky分解的有限元正演计算，得到N个模型的观测数据，并求解每一个模型对应的观测数据拟合差.
(6) 找出本次最小拟合差对应的模型，即本次选择中最优的解.判断是否满足收敛条件，如不满足收敛条件，则按照本次最优模型对相应网格单元的电阻率节点进行信息素强度更新，重复以上过程.
(7) 如满足收敛条件，则进入最小二乘线性反演.
3.2 最小二乘线性反演阶段(1) 将非线性反演所得的模型作为线性反演的初始模型，利用有限元正演得到理论观测数据.
(2) 计算正演观测数据与实际观测数据的拟合差，并判断是否满足收敛条件.
(3) 如不满足收敛条件，则求解反演方程，得到新的模型参数，进行下一次迭代循环，直至满足收敛条件.
(4) 若观测数据与理论数据的拟合差小于反演收敛的容许误差，则将此时得到的模型作为最终结果输出，混合反演算法结束.
4 数值反演算例为了评价基于GPU并行的蚁群算法与最小二乘方法相结合的混合反演计算速度，并检验基于GPU混合反演在隧道超前探测中的成像效果，开展了2个隧道超前探测的数值算例，其中1个算例对比GPU并行计算与串行顺序计算的计算效率，并探讨混合反演与最小二乘线性反演的成像结果.另外1个算例用于评价基于GPU混合反演方法在探测复杂倾斜含水构造的成像效果.
4.1 算例一设隧道断面宽12 m×高8 m，隧道空腔电阻率ρ1=1×107 Ωm，围岩电阻率ρ2=1000 Ωm.含水溶洞的大小为8 m×8 m×4 m，其电阻率为ρ3=10 Ωm.溶洞的几何形心位于隧道掌子面前方22 m处，如图 4所示，以掌子面中心点为坐标原点建立坐标系，高度方向为X轴，宽度方向为Y轴，开挖方向为Z轴.
|
图 4 含水溶洞超前探测示意图 Fig. 4 Model for ahead prospecting of karst cave filled with water |
(1) GPU加速蚁群算法的计算速度对比
由于混合反演中大量的计算耗时花费在非线性反演中，为了对比混合反演的GPU并行计算效果，对蚁群非线性计算过程进行GPU并行计算与串行计算计算耗时对比.
建立三维有限元计算模型，计算网格数量为23800，由于获得的理论观测数据中视电阻率值基本在1000 Ωm，因此设定计算区域中模型网格电阻率初值为1000 Ωm.反演网格的数量为3520个，即X、Y与Z方向上网格剖分的数量分别是8、22与20个.分别采用基于GPU并行蚁群算法与串行的蚁群算法进行非线性反演.蚁群算法中蚂蚁的数量为50个，初始信息素强度为1.0.计算测试在一台有4个CPU和2个GPU的服务器上.CPU主频为2.50 GHz，10核心；GPU显卡型号为NVidia Tesla K20X，流处理单元为2688个，显存容量为6 GB.
对蚁群算法进行了8次反演迭代，GPU并行计算总耗时约416 min，CPU单线程串行蚁群算法总耗时约3320 min，表 1给出了两种方法在一次反演迭代中的计算时间以及主要耗时单元的时间消耗.与CPU单线程串行蚁群算法反演的计算效率对比可发现，采用并行算法后，蚂蚁选择模型的数值正演耗时相比CPU单线程串行计算时间大幅度降低，计算效果得到显著的提高.在一次迭代中，50个模型同时进行有限元正演耗时约45 min，相比串行时所用的400 min，计算效率提高了约7.89倍.在蚂蚁选择模型阶段，采用GPU并行策略后，时间由串行时的9 min缩短到5 min，计算效率提高了约0.8倍.在对信息素强度更新采取并行策略后，时间由串行时的5 min提速到2 min，计算效率提高1.5倍.由此可见，基于cholesky分解的有限元正演是蚁群算法中最耗时的部分，对正演部分进行并行改造是十分有必要的.GPU并行改进的蚁群算法总耗时是传统串行蚁群算法的约1/8，在提高计算效率方面的效果是十分显著的.
|
表 1 一次迭代中两种计算模式各部分计算耗时对比 Table 1 Computing time using two inversion model in one iteration |
(2) 基于GPU混合反演的超前探测成像结果
为了检验基于GPU混合反演超前探测的成像效果，采用与4.1节中相同的模型，进行GPU混合反演超前探测的数值模拟，其他参数与4.1节相同.基于GPU并行的混合反演结果如图 5所示，在掌子面前方存在条带状低阻区域，低阻区域在Y轴方向宽度约为9 m，这与实际模型在Y轴的宽度8 m接近.低阻区域在Z轴方向厚度约为4.5 m，位于掌子面前方约21~25 m，这与实际模型在Z轴方向的厚度4 m且位于掌子面前方20~24 m较为接近.且能够反映出宽度与厚度约2:1的条带状低阻体，与实际模型的宽厚比一致.
|
图 5 混合反演结果 (a)三维反演结果; (b)异常体提取三维图; (c) x=0 m处切片，图中红色虚框表示反演结果的低阻异常区，红色竖线表示真实异常体在Z轴方向的中心位置. Fig. 5 Joint inversion results using GPU parallel on ACO algorithm and the least-squares method (a) 3D view; (b) 3D graph of anomaly-body; (c) Slice at x=0 m. The red square represents low resistivity region in the inversion result, the red dotline represents the center of the abnormal body in the Z-axis direction. |
(3) 基于最小二乘线性反演的超前探测成像结果
为了与最小二乘线性反演的成像结果进行对比，仍采用与4.1节中相同的模型与计算参数，进行最小二乘线性反演超前探测数值模拟.反演结果见图 6，反演图像中有较大范围的低阻区域，能够反映出掌子面前方低阻体的存在，但存在较大问题，具体表现为：1)反演结果中除设定的低阻模型外，出现了两处明显的低阻异常，判定为假异常，给解译工作带来了干扰.2)低阻区域是原型中含水溶洞的反映，但反演出的形态近似于圆形，与原型中宽厚比为2:1的长方形条带状存在较大差异.3)低阻区域在Y轴方向的宽度约为12 m，与实际模型在Y轴方向的宽度8 m相比差异较大.低阻区域在Z轴方向的厚度约为10 m，与实际模型在Z轴的厚度4 m相差较大.
|
图 6 最小二乘反演结果 (a)三维反演结果; (b)异常体提取三维图; (c) x=0 m处切片. Fig. 6 Inversion results using the least-squares method (a) 3D view; (b) 3D graph of anomaly-body; (c) Slice at x=0 m. |
对比两种反演结果发现，最小二乘算法的反演结果(图 6)能够大致反映含水溶洞的存在，但异常体的区域较大，定位不准确，其规模与位置，与原模型相比存在较大差异.而GPU并行的混合反演相比最小二乘反演，结果有了较大改进，低阻区呈近似长方形条带状，其规模与位置与实际模型较为一致.另外，混合反演的反演结果中没有假异常的出现.这表明混合反演在提高三维电阻率的反演质量、降低多解性方面有积极作用.
4.2 算例二为了探讨基于GPU并行混合反演探测复杂异常体的成像效果，设计一个倾斜含水构造的探测模型，如图 7所示.位于隧道掌子面前方8 m处设计了一个倾斜含水断层，其在Y轴方向的宽度为12 m，在Z轴方向的长度为8 m，在X轴方向的高度为8 m.背景电阻率为1000 Ωm，含水断层的电阻率值为10 Ωm.
|
图 7 含水断层超前探测示意图 Fig. 7 Model for ahead prospecting of fault filled with water |
(1) 基于GPU混合反演的超前探测成像结果
同样，分别采用基于GPU并行蚁群算法与串行的蚁群算法分别对合成观测数据进行反演.蚁群算法中，蚂蚁数目为50只，反演迭代8次，计算参数与算例一均相同.表 2给出了一次迭代过程中两种计算模式的主要耗时单元耗时对比以及八次迭代总耗时对比.同样发现，蚁群算法串行计算总耗时约3440 min，基于GPU并行化蚁群算法总耗时约440 min，计算时间大幅度降低.与传统CPU单线程串行计算相比，在一次迭代过程中，蚁群算法各主要单元的耗时经过GPU并行化处理之后，计算效率显著提升.
|
|
表 2 两种计算模式耗时情况对比 Table 2 Computing time using two inversion methods |
基于GPU并行的混合反演结果如图 8所示，在三维电阻率反演图像上出现明显的低电阻区域，且其形状为倾斜条带状.倾斜条带状低阻体位于掌子面前方约9~16 m范围(Z轴方向)，这与实际异常体在Z轴位置8~16 m接近.倾斜条带状低阻体在Y轴的宽度约为10 m，这与实际模型宽度12 m相差不大.需要指出的是反演出的倾斜条带状低阻体的倾斜角度与实际模型接近.可见，基于GPU并行的混合反演方法能够反映出倾斜含水断层，且能够反映出倾斜含水断层的位置、规模.
|
图 8 混合反演结果 (a)三维反演结果; (b)异常体提取三维图; (c) x=0 m处切片. Fig. 8 Joint inversion results using GPU parallel on ACO algorithm and the least-squares method (a) 3D view; (b) 3D graph of anomaly-body; (c) Slice at x=0 m. |
(2) 基于最小二乘线性反演的超前探测成像结果
同时开展了最小二乘线性反演，其结果见图 9.反演图像中有较大范围的低阻区域，能够反映倾斜含水断层的存在.但反演出的异常体形态的倾斜角度与原型不相符，将倾斜断层反演成了近似圆形的团状低阻体.另外，反映出的低阻体位于掌子面左侧，与实际模型的位置相差较大.
|
图 9 最小二乘反演结果 (a)三维反演结果; (b)异常体提取三维图; (c) x=0 m处切片. Fig. 9 Inversion results using the least-squares method (a) 3D view; (b) 3D graph of anomaly-body; (c) Slice at x=0 m. |
由此可见，在掌子面前方存在较复杂含水构造(如倾斜含水断层)时，基于GPU并行的混合反演结果可以反演出含水构造的位置、规模与形态，而仅采用基于最小二乘的线性反演方法，可以反映出含水构造的存在，但在含水构造的位置与规模方面分辨能力较差.因此，基于GPU的混合反演对于提高三维电阻率的反演精度与降低多解性有积极作用.
5 三维电阻率超前探测模型试验研究为了进一步验证基于GPU并行的蚁群算法与最小二乘方法相结合的混合反演的实际应用效果，开展了隧道超前探测的物理模型试验.
5.1 模型试验装置与观测系统电阻率物理模型试验需满足几何相似及电阻率相似.几何相似指模型包括尺寸、形状、电极极距等参数与实际几何相似；电阻率相似指模型围岩与异常体的电阻率，与实际工程现场围岩与异常体的电阻率相似.本模型试验的几何相似比是1:5，如图 10所示.模型的尺寸为17.0 m×8.4 m×6.7 m，隧道腔体采用玻璃钢材料制作，尺寸为2.0 m×1.7 m×6.0 m(聂利超等，2015).含水构造的尺寸为1.2 m×1.0 m×0.4 m.含水构造中充填一定比例的泥沙和水，控制其电阻率为10 Ωm，位于隧道掌子面前方1.6 m处.围岩采用已配置好的相似材料，分层压实，对试验区进行填充，电阻率为200 Ωm.
|
图 10 物理模型试验装置图 Fig. 10 Device of physical model test |
采用多同性源阵列电阻率超前探测的观测方式，如图 11所示.在隧道掌子面上布置2条测线，测线间距为0.4 m，每条测线上布置10个测量电极(M1—M10，M11—M20)，其中一个电极作为供电电极，电极间距为0.15 m.在测量电极的外围掌子面周边布置4个供电电极(A1—A4)，在测量过程中4个供电电极逐次向后移动.
|
图 11 物理模型试验观测方式示意图 Fig. 11 Observation system in physical model test |
利用本文提出的并行蚁群算法与最小二乘算法相结合的混合反演算法对实验观测数据进行反演，反演网格数量在X，Y，Z方向上分别为8、21与16，共2688个.蚁群算法中，蚂蚁数量为50个，初始信息素强度为1.0，反演迭代10次，共耗时约500 min(传统非并行蚁群算法的耗时为4000 min，是并行蚁群算法的8倍).基于GPU并行的混合反演方法的三维电阻率反演结果见图 12，反演图像中存在明显的低阻区域，在Z轴方向位于约1.5~2.2 m范围，与实际模型在Z轴方向的位置1.6~2.0 m较接近.低阻区域在Y轴方向宽度约为1.4 m，这与实际模型在Y轴宽度1.2 m接近.可见，基于GPU混合反演可以反演出低阻体的规模与位置，且与实际模型较为一致.
|
图 12 混合反演结果 (a)三维反演结果; (b)异常体提取三维图; (c) x=0 m处切片. Fig. 12 Joint inversion results using GPU parallel on ACO algorithm and the least-squares method (a) 3D view; (b) 3D graph of anomaly-body; (c) Slice at x=0 m. |
同时采用模型试验数据进行了传统最小二乘线性反演，反演结果如图 13所示.反演结果能够反映含水构造的存在，其位置与真实模型相差较大.低阻区在Y轴方向的位置位于约为-1.8 m至0.2 m范围内，与实际模型在Y轴方向的位置-0.6 m至0.6 m差异较大，其在Z轴方向位置约为0.6~2.2 m也与实际模型在Z轴方向的位置1.6~2.0 m偏差较大.
|
图 13 最小二乘反演结果 (a)三维反演结果; (b)异常体提取三维图; (c) x=0 m处切片. Fig. 13 Inversion results (a) 3D view; (b) 3D graph of anomaly-body; (c) Slice at x=0 m. |
综合对比图 12，图 13可知，基于GPU并行的混合反演结果相比最小二乘的反演结果有了明显改善.在模型试验的实际应用中，混合反演结果中低阻区域能够反映含水构造的大体位置，且位置与模型较为一致.
6 结论(1) 针对隧道三维电阻率超前探测成像与快速化解译的需求，提出了基于GPU并行的蚁群算法与最小二乘相结合的混合反演方法.利用蚁群算法的全局搜索能力强的优势为最小二乘线性反演提供较优的初始模型，克服最小二乘算法易陷入局部最优的缺点.同时，基于GPU并行加速混合反演算法，提高了隧道三维电阻率反演成像的精度，实现了三维电阻率的快速化成像.
(2) 数值算例将基于GPU混合反演方法与最小二乘线性反演方法进行了对比，发现基于GPU并行的混合反演方法可以更好的反映出含水构造的规模、位置与形态，且与实际模型较为接近.相比较最小二乘线性反演结果，基于GPU混合反演在异常体定位、识别与快速化成像方面具有显著优势.
(3) 基于GPU并行的混合反演方法在隧道超前探测物理模型试验中进行了应用，混合反演方法能够较准确的反映低阻体的位置与规模，表明基于GPU并行的混合反演方法在模型试验的实际应用中得到较好的效果，具有在工程实际中应用的前景，为解决隧道超前探测快速化解译与成像提供了一条可行途径.
Angelo J S. 2013. Strategies for parallel ant colony optimization on graphics processing units.//Ant Colony Optimization-Techniques and Applications. InTech, 63-83.
|
|
Chen Z X, Meng X H, Guo L H, Liu G F, et al.
Three-dimensional fast forward modeling and the inversion strategy for large scale gravity and gravimetry data base on GPU. Chinese J. Geophys, 55(12): 4069-4077.
|
|
Cheng J L, Wang Y H, Yu S J, et al.
2000. The principle and application of advance surveying in roadway excavation by resistivity method. Coal Geology and Exploration, 28(4): 60-62.
|
|
Colorni A, Dorigo M, Maniezzo V. 1992. An investigation of some properties of an "ant aigorithm".//Proc. of the Parallel Problem Solving from Nature Conference (PPSN92). Blegium:Elsevier Publishing, 509-520.
|
|
Conti C R, Roisenberg M, Neto G S, et al.
2013. Fast Seismic inversion methods using ant colony optimization algorithm. IEEE Geoscience & Remote Sensing Letters, 10(5): 1119-1123.
|
|
Dai Q W, Jiang F B, Dong L.
2014. RBFNN inversion for electrical resistivity tomography based on Hannan-Quinn criterion. Chinese Journal of Geophysics, 57(4): 1335-1344.
DOI:10.6038/cjg20140430 |
|
Delévacq A, Delisle P, Gravel M, et al.
2013. Parallel ant colony optimization on graphics processing units. Journal of Parallel and Distributed Computing, 73(1): 52-61.
DOI:10.1016/j.jpdc.2012.01.003 |
|
Dorigo M. 1992. Optimization, learning and natural algorithms[Ph. D. thesis]. Italy:Politecnico Di Milano.
|
|
Gao M L, Yu B S, Zheng J B, et al.
2016. Research of resistivity imaging using neural network based on immune genetic algorithm. Chinese J. Geophys., 59(11): 4372-4382.
DOI:10.6038/cjg20161136 |
|
Huang J G, Ruan B Y, Wang J L.
2007. The fast inversion for advanced detection using DC resistivity in tunnel. Chinese Journal of Geophysics, 50(2): 619-624.
|
|
Kaus A, Boening W.
2008. BEAM-Geoelectrical ahead monitoring for TBM-drives. Geomechanics and Tunneling, 1(5): 442-450.
DOI:10.1002/geot.v1:5 |
|
Khalil M A, Santos F A M.
2013. 2D and 3D resistivity inversion of Schlumberger vertical electrical soundings in Wadi El Natrun, Egypt:A case study. Journal of Applied Geophysics, 89: 116-124.
DOI:10.1016/j.jappgeo.2012.11.014 |
|
Li S C, Liu B, Sun H F, et al.
2014. State of ART and trends of advanced geological prediction in tunnel construction. Chinese Journal of Rock Mechanics and Engineering, 33(6): 1090-1113.
|
|
Li S C, Nie L C, Liu B, et al.
2015. Advanced detection and physical model test based on multi-electrode sources array resistivity method in tunnel. Chinese Journal of Geophysics, 58(4): 1434-1446.
DOI:10.6038/cjg20150429 |
|
Li S C, Wang C W, Nie L C, et al.
2016. Parallel 3D electrical resistivity inversion method with inequality constraint based on slack variables. Chinese Journal of Rock Mechanics and Engineering, 35(6): 1122-1132.
|
|
Liu B, Li S C, Li S C, et al.
2012a. 3D electrical resistivity inversion with least-squares method based on inequality constraint and its computation efficiency optimization. Chinese J. Geophys., 55(1): 260-268.
DOI:10.6038/j.issn.0001-5733.2012.01.025 |
|
Liu B, Li S C, Nie L C, et al.
2012b. Advanced detection of water-bearing geological structures in tunnels using 3D DC resistivity inversion tomography method. Chinese Journal of Geotechnical Engineering, 34(10): 1866-1876.
|
|
Liu B, Liu Z Y, Song J, et al.
2017. Joint inversion method of 3D electrical resistivity detection based on inequality constraints. Chinese J. Geophys., 60(2): 820-832.
DOI:10.6038/cjg20170232 |
|
Liu B, Li S C, Nie L C, et al.
2012. 3D resistivity inversion using an improved Genetic Algorithm based on control method of mutation direction. Journal of Applied Geophysics, 87: 1-8.
DOI:10.1016/j.jappgeo.2012.08.002 |
|
Nie L C, Li S C, Liu B, et al.
2015. Development of comprehensive model similitude material for multiple geophysical detection. Chinese Journal of Rock Mechanics and Engineering, 34(S1): 3087-3096.
|
|
Rennich S C, Stosic D, Davis T A. 2014. Accelerating sparse cholesky factorization on GPUs.//Proceedings of the 4th Workshop on Irregular Applications:Architectures and Algorithms. New Orleans, Louisiana, 9-16.
|
|
Singh U K, Tiwari R K, Singh S B.
2010. Inversion of 2-D DC resistivity data using rapid optimization and minimal complexity neural network. Nonlinear Processes in Geophysics, 17(1): 65-76.
DOI:10.5194/npg-17-65-2010 |
|
Taherian A R.
2015. Experiences of TBM operation in gas bearing water condition-A case study in Iran. Tunnelling & Underground Space Technology, 47: 1-9.
|
|
Wan X L, Xi D Y, Gao E G, et al.
2005. 3-D resistivity inversion by the least-squares QR factorization method under improved smoothness constraint condition. Chinese J. Geophys., 48(1): 439-444.
|
|
Wu X P, Xu G M.
2000. Study on 3-D resistivity inversion using conjugate gradient method. Chinese J. Geophys., 43(3): 420-427.
|
|
Xu H L, Wu X P.
2006. 2D resistivity inversion using the neural network method. Chinese J. Geophys., 49(2): 584-589.
|
|
Zhang L Y. 2011. The study of nonlinear inversion method in high-density resistivity method inversion[Ph. D. thesis] (in Chinese). Taiyuan University of Technology.
|
|
Zhang L Y, Liu H F.
2011. The application of ABP method in high-density resistivity method inversion. Chinese J. Geophys., 54(1): 227-233.
DOI:10.3969/j.issn.0001-5733.2011.01.024 |
|
陈召曦, 孟小红, 郭良辉, 等.
2012. 基于GPU并行的重力、重力梯度三维正演快速计算及反演策略. 地球物理学报, 55(12): 4069–4077.
DOI:10.6038/j.issn.0001-5733.2012.12.019 |
|
程久龙, 王玉和, 于师建, 等.
2000. 巷道掘进中电阻率法超前探测原理与应用. 煤田地质与勘探, 28(4): 60–62.
|
|
戴前伟, 江沸菠, 董莉.
2014. 基于汉南-奎因信息准则的电阻率层析成像径向基神经网络反演. 地球物理学报, 57(4): 1335–1344.
DOI:10.6038/cjg20140430 |
|
高明亮, 于生宝, 郑建波, 等.
2016. 基于IGA算法的电阻率神经网络反演成像研究. 地球物理学报, 59(11): 4372–4382.
DOI:10.6038/cjg20161136 |
|
黄俊革, 阮百尧, 王家林.
2007. 坑道直流电阻率法超前探测的快速反演. 地球物理学报, 50(2): 619–624.
|
|
李术才, 刘斌, 孙怀凤, 等.
2014. 隧道施工超前地质预报研究现状及发展趋势. 岩石力学与工程学报, 33(6): 1090–1113.
|
|
李术才, 聂利超, 刘斌, 等.
2015. 多同性源阵列电阻率法隧道超前探测方法与物理模拟试验研究. 地球物理学报, 58(4): 1434–1446.
DOI:10.6038/cjg20150429 |
|
李术才, 王传武, 聂利超, 等.
2016. 基于松弛变量的不等式约束三维电阻率并行反演方法研究. 岩石力学与工程学报, 35(6): 1122–1132.
|
|
刘斌, 李术才, 李树忱, 等.
2012a. 基于不等式约束的最小二乘法三维电阻率反演及其算法优化. 地球物理学报, 55(1): 260–268.
DOI:10.6038/j.issn.0001-5733.2012.01.025 |
|
刘斌, 李术才, 聂利超, 等.
2012b. 隧道含水构造直流电阻率法超前探测三维反演成像. 岩土工程学报, 34(10): 1866–1876.
|
|
刘斌, 刘征宇, 宋杰, 等.
2017. 基于不等式约束的三维电阻率探测混合反演方法. 地球物理学报, 60(2): 820–832.
DOI:10.6038/cjg20170232 |
|
聂利超, 李术才, 刘斌, 等.
2015. 多元地球物理综合探测模型试验相似材料的研制. 岩石力学与工程学报, 34(S1): 3087–3096.
|
|
宛新林, 席道瑛, 高尔根, 等.
2005. 用改进的光滑约束最小二乘正交分解法实现电阻率三维反演. 地球物理学报, 48(1): 439–444.
|
|
吴小平, 徐果明.
2000. 利用共轭梯度法的电阻率三维反演研究. 地球物理学报, 43(3): 420–427.
|
|
徐海浪, 吴小平.
2006. 电阻率二维神经网络反演. 地球物理学报, 49(2): 584–589.
|
|
张凌云. 2011. 高密度电阻率勘探反演的非线性方法研究[博士论文]. 太原: 太原理工大学.
http://cdmd.cnki.com.cn/Article/CDMD-10112-1011082080.htm |
|
张凌云, 刘鸿福.
2011. ABP法在高密度电阻率法反演中的应用. 地球物理学报, 54(1): 227–233.
DOI:10.3969/j.issn.0001-5733.2011.01.024 |
|