基于滤波算法的隧道全时域激发极化超前探测正演方法

第 55 卷第 2 期2024 年 2 月
中南大学学报(自然科学版)
Journal of Central South University (Science and Technology)
V ol.55 No.2Feb. 2024
基于滤波算法的隧道全时域激发极化超前探测正演方法
李志强1，褚国安1，聂利超1，裴文兵1，郭一凡1，周维1，胡军2，闫冰1，邓朝阳1
(1. 山东大学 岩土与结构工程研究中心，山东 济南，250061；2. 云南省滇中引水二期工程有限公司，云南 昆明，650205)
摘要：为进一步提高隧道全时域激发极化多参数超前探测的效果，采用滤波算法进行计算时间域等效电阻率，采用有限单元法进行计算，提出一种基于滤波算法的隧道全时域激发极化超前探测正演方法。

该方法针对时间域Cole-Cole 表达式计算收敛性差的问题，引入滤波算法进行时间域激发极化等效电阻率的计算，提高算法的收敛性；针对正演过程中大型矩阵求解影响计算效率的问题，引入OpenMP 并行算法，实现隧道全时域激发极化超前探测正演的加速运算。

最后，开展数值算例和物理模型试验。

研究结果表明：该正演方法对异常体有敏感响应，且在加快正演速度具有显著优势。

通过基于滤波算法的隧道全时域激发极化超前探测获得的观测数据与实际探测数据较为相符，这为基于全时域激发极化隧道超前探测提供了一种可行的正演算法。

关键词：隧道超前地质预报；全时域激发极化；滤波算法；快速三维正演模拟；并行算法中图分类号：P631 文献标志码：A 文章编号：1672-7207（2024）02-0638-11
Forward method for full time domain induced polarization
advance detection of tunnel based on filtering algorithm
LI Zhiqiang 1, CHU Guoan 1, NIE Lichao 1, PEI Wenbing 1, GUO Yifan 1, ZHOU Wei 1,
HU Jun 2, YAN Bing 1, DENG Zhaoyang 1
(1. Geotechnical & Structural Engineering Research Center of Shandong University, Jinan 250061, China;
2. Yunnan Water Diversion Phase II Project Co. Ltd., Kunming 650205, China)
收稿日期： 2023 −06 −14； 修回日期： 2023 −10 −18
基金项目(Foundation item)：国家重点研发专项课题(2021YFC3100805)；国家自然科学基金资助项目(52079076)；山东省自然
科学基金资助项目(ZR2022YQ46，ZR2022QD014)；云南省科技厅重大科技专项计划项目(202002AF080003) (Project (2021YFC3100805) supported by the National Key Research and Development Program of China; Project(52079076) supported by the National Natural Science Foundation of China; Projects(ZR2022YQ46, ZR2022QD014) supported by the Natural Science Foundation of Shandong Province; Project(202002AF080003) supported by the Key Science and Technology Special Program of Yunnan Province)
通信作者：聂利超，博士，教授、博士生导师，从事地球物理正反演理论与隧道不良地质超前探测的研究；E-mail ：
*****************
DOI: 10.11817/j.issn.1672-7207.2024.02.017
引用格式： 李志强, 褚国安, 聂利超, 等. 基于滤波算法的隧道全时域激发极化超前探测正演方法[J].中南大学学报(自然科学版), 2024, 55(2): 638−648.
Citation: LI Zhiqiang, CHU Guoan, NIE Lichao, et al. Forward method for full time domain induced polarization advance detection of tunnel based on filtering algorithm[J]. Journal of Central South University(Science and Technology), 2024, 55(2): 638−
648.
第 2 期李志强，等：基于滤波算法的隧道全时域激发极化超前探测正演方法
Abstract:In order to meet the requirements of high detection efficiency in actual tunnel engineering, the filtering algorithm was used to calculate the equivalent resistivity in the time domain, and the finite element method was used for calculation. A forward modeling method of full-time induced polarization advance tunnel detection based on the filtering algorithm was proposed. To address the problem of poor convergence of Cole-Cole expression in time domain, the filtering algorithm was introduced to calculate the equivalent resistivity in time domain induced polarization, improving the convergence of algorithm. To solve the problem that the large-scale matrix solution affects the computational efficiency in the forward modeling process, the OpenMP parallel algorithm was introduced to realize the accelerated operation of the full-time domain induced polarization advanced detection forward modeling of the tunnel. Finally, numerical examples and physical model tests were carried out. The results show that the forward modeling method has a sensitive response to abnormal bodies. The full-time domain induced polarization forward acceleration algorithm has significant advantages in accelerating forward modeling speed. The forward data obtained by the full-time domain induced polarization forward detection method based on the filtering algorithm is in good agreement with the actual detection data, which lays a foundation for the application of the full-time domain induced polarization inversion method in the actual project of tunnel geological prediction.
Key words: tunnel advanced prediction; full-time domain induced polarization; filtering algorithm; fast three-dimensional forward modeling; parallel algorithm
近年来，随着中国基础设施建设和西部大开发进程的加快，水利水电、市政工程等领域所修建的隧道数量大幅增长。

在隧道建设过程中，常遭遇裂隙发育带、断层破碎带、软弱地层、溶洞、暗河等不良地质情况，易诱发突水突泥、塌方、大变形等地质灾害[1−3]。

给隧道的安全高效建设带来极大挑战。

因此，提前探明掌子面前方不良地质至关重要。

当前隧道地质超前预报技术主要有地震波类、电磁类、电法类等[4−6]。

激发极化法因其对水体的响应敏感，且受地形影响小的特点，常被用来探测地下水赋存情况。

CHIRINDJA等[7]采用激发极化法与电阻率法成功探明了莫桑比克楠普拉省沿海含水层的构成。

LI等[8]综合运用地质分析、激发极化与地震超前地质预报方法，对TBM施工隧洞掌子面前方的溶洞、断层破碎带等不良地质进行了准确预报，避免了隧洞突水突泥、塌方、卡机灾害的发生。

刘春生等[9]采用巷道钻孔激发极化超前探水方法，解决了煤矿突涌水灾害问题。

LING等[10]采用激发极化法，结合水文地质资料和钻孔信息，成功划分含水层，为激发极化法在城市地下水调查的应用提供了理论支持。

目前激发极化法多利用电阻率、极化率、衰减时差等激发极化充放电曲线的部分信息，而全时域激发极化法是当前的研究热点，利用了整个充放电过程的信息，可实现描述充放电过程的激发极化多参数成像，可从多角度对目标体进行较准确刻画[11−13]。

分析全时域激发极化信息时，需要一个能够准确描述时间域激发极化场时变效应、且携带激发极化谱参数的模型，目前常用的模型是时间域Cole-Cole模型。

然而，由于时间域Cole-Cole模型表达中存在无穷级数，当t/τ(t为衰减时间，τ为弛豫时间)在2π附近时，级数的收敛性差，当频率相关系数c较小时，收敛性更差，而常见的岩土体及矿物的频率相关系数一般不大，这导致采用时间域Cole-Cole模型的正演算法在实际应用时常存在难以收敛的情况。

针对上述问题，相关学者对时频转换进行了研究。

KEMNA等[14]将观测数据转化至频率域进行求解，实现将复电阻率成像应用于任意电极排列；FIANDACA等[15]采用时频转换的方法，将频率域复电阻率的正演方程作为正演核，然后将频率域数据转化到时间域，与全时域的观测数据进行拟合，有效提高了全时域激发极化正演算法的收敛性，但其需要经过复杂的傅里叶变换，过程繁琐，影响正演效率[14−16]。

因此如何获得高精度的快速三维正演模拟结果仍然是全时域激化极化法的一个技术难点，在一定程度上制约了全时域激发极化技术的应用和发展。

GUPTASARMA[17]提出了基于线性数字滤波求取时间域Cole-Cole模型解的方法，其计算精度高，能够满足正演的要求，因而也得到了广泛应用[18−19]，
639
第 55 卷
中南大学学报(自然科学版)MENG 等[20]提出了基于Weibull 增长模型的激发极化正演算法，这为在宏观尺度上进一步阐释激发极化效应方面的研究提供了新的途径。

以上学者的研究，为研究隧道超前预报的全时域激发极化正演算法提供了参考。

针对隧道激发极化正演的收敛性问题，本文作者提出了一种基于滤波算法的隧道全时域激发极化超前探测正演方法。

引入滤波算法，解决了时间域Cole-Cole 计算过程中存在无穷级数的问题，有效提高了算法的收敛性；针对全时域多时窗正演耗时较多的难题，采用数组映射方法以此实现对于某节点周围节点的快速查找。

同时，引入OpenMP 并行算法解决了正演计算速度的问题，实现了隧道全时域激发极化超前探测正演加速运算。

最后，通过数值与物理模型试验验证了该全时域激发极化正演方法的可行性。

1 隧道全时域激发极化正演方法
全时域激发极化正演计算包含多时窗，每个时窗均需进行一次完整电阻率正演，本文采用滤波算法计算了时间域等效电阻率，然后对等效电阻率进行正演模拟，在等效电阻率基础上求取了电位曲线。

提出了基于滤波算法的隧道全时域激发极化超前探测正演方法。

1.1　基于滤波算法的等效电阻率求解
PELTON 等[21]将Cole-Cole 模型转至时间域进行计算，Cole-Cole 模型能够准确表达激发极化参数之间的相关性，Cole-Cole 模型时间域的表达式为：
η(t )=η0éëêê∑k =1
¥(-1)k -1
x kc Γ()1+kc ùûú
ú(1)ρl (t )=ρ0(1-η(t ))
(2)
式中：Γ(x )为Gamma 函数；η0为固有极化率；c 为频率相关系数；k =1，2，3，…；x=t /τ，ρl (t )为在t 时刻的等效电阻率；ρ0为零频电阻率。

本文采用滤波算法计算了时间域等效电阻率，依据滤波算法公式，涉及频域数据到时域数据的传递函数为[17]：
ρ(ω)=ρ0(1-η0éë
êêê1-11+()i ωτc ùûúúú)
(3)
式中：ω为角频率。

由传递函数得到的时间域等效电阻率的计算公式为[17]：
ρl (t )=∑r =1
21ξr Â[]
ρl (ωr )(4)
式中：ξr 为线性滤波系数，Â[ρl (ωr )]表示对ρl (ωr )取实部；ωr =10(αr -lg t )，其中，αr 为位移系数。

对于全时域激发极化正演的研究多采用时频转换的方式，过程繁琐。

另外，由于Cole-Cole 模型表达式中包含无穷级数，当时间域弛豫时间的比值在2π附近时，级数的收敛性较差。

而本文提出的滤波算法将频域数据转至时域数据，可以很好地解决收敛问题。

为了验证本文方法求取等效电阻率的准确性，分别通过时间域Cole-Cole 模型和滤波算法计算等效电阻率，设置频率相关系数为0.5，弛豫时间为0.3 s ，固有极化率为0.2，零频电阻率为100 Ω·m ，通过MATLAB 编程计算等效电阻率，结果如图1所示。

由图1可得时间域Cole-Cole 模型与滤波算法计算结果均差为10−5，两者精度相近。

但由于Cole-Cole 模型等效电阻率公式中含有无穷级数，因此，收敛性较差，且收敛耗时较长，本文进行了时间域Cole-Cole 模型和滤波算法计算等效电阻率的耗时对比。

在相同计算条件下，时间域Cole-Cole 模型耗时为1 s ，而滤波算法计算耗时0.1 s ，另外，滤波算法有效避免了时间域Cole-Cole 模型无穷级数不收敛问题，提升了计算效率。

1.2　隧道全时域激发极化正演模拟
HORDT 等[22]
提出时域激电正演模拟可以忽略
图1　时间域Cole-Cole 模型与滤波算法数值结果对比Fig. 1　Comparison of numerical results between Cole-Cole model of time domain and filtering algorithm
640
第 2 期李志强，等：基于滤波算法的隧道全时域激发极化超前探测正演方法随时间变化的电导率与电场的卷积，可以采用直
流电阻率法分别计算各时窗电位数据，这种近似
方法极大地简化了时域激发极化数据的正演数值
模拟过程。

本文采用直流电阻率正演算法计算各
时窗数据，将使用滤波算法得到的等效电阻率进
行正演模拟。

在三维电场中，当点电源位于隧道内时，则
可近似为全空间条件，此时均匀全空间下的电位
应满足的微分方程为[23]：
Ñ·(σÑv)=-Iδ(A)(5)
式中：Ñv为全空间下电位v的梯度；σ为介质的电
导率，I为供电电流强度，δ(A)是以A为中心的δ
函数。

当地下为均匀介质时，由式(5)可解得全空间
正常电位u
为[23]：
u 0=
I
4πrσ
(6)
式中：r为测点至电源的距离；σ
为电源处的介质
电导率。

假设在电导率为σ
1
的均匀介质中，存在一个
电导率为σ
2
的异常体。

通过滤波算法计算等效电
导率σ
1l 和σ
2l
：
σ
1l
=
1
ρ
1l
(t)
(7)
σ
2l
=
1
ρ
2l
(t)
(8)
电源处介质等效电导率为σ
0=σ
1l
；v
1
和v
2
分别
表示σ
1l 和σ
2l
所占区域Ω
1
和Ω
2
中的总电位；u
1
和u
2
分别表示Ω
1和Ω
2
区域中的异常电位，有[23]：
ì
í
î
v
1
=u
+u
1
v
2
=u
+u
2
(9)
异常电导率σ′=σ-σ
，令[23]
Ñ·(σ0Ñu0)=-Iδ(A)(10)可得异常电位u满足的微分方程为[23]：
Ñ·(σÑu)=-Ñ·(σ′Ñu0)(11)为简化计算，用八节点六面体的规则六面体单元进行建模，对求解区域的网格进行加密。

将求解区域剖分后，异常电位的变分问题中的积分被离散为各单元中的积分，则可得线性方程组[23]：
KU=--K′U
(12)式中：K为求解总电位的总体系数矩阵，K′为求解异常电位的总体系数矩阵；U为全部节点的异常电
位向量，U
0为全部节点的正常电位向量。

求解方程(12)可得到各节点的异常电位，从而
得到各节点的总电位为[23]：
V=U+U
(13)
式中，V为全部节点的总电位向量。

采用多元非线性最小二乘拟合，由Cole-Cole
模型作为数据的拟合模型，通过测量参数与Cole-
Cole模型的拟合结果，计算视零频电阻率、视固
有极化率、视弛豫时间等激发极化参数。

2 全时域激发极化正演加速算法
三维全时域激发极化正演算法中需进行多个
时刻的电阻率正演，而三维电阻率正演通常存在
耗时多、效率低的问题，严重影响算法在隧道探
测时的实用性，因此，有必要进行正演算法的加
速。

在进行三维电阻率正演计算时，需要调用三
维模型中各个网格节点及其周围节点，且每次正
演需要遍历全部节点。

然而，常规的存储方式会
造成正演过程中包含一部分冗余节点信息，造成
较大的内存占用；并且正演计算过程涉及大型矩
阵求解，造成耗时多、效率低。

针对上述问题，
本文使用数组映射的存储方式，同时引入并行算
法，实现对正演过程的加速运算。

对于本文使用的六面体八节点单元(图2[23])，
冗余节点存在的原因是，对于每个节点默认其相
邻节点为27个，然而例如边界上的节点其相邻节
点数目并不足27个；同时，由于之前的存储方式
未对相邻节点分组存储，每次调用都需要从头开
始。

为了解决这两个问题，需要对相邻节点分组
存储，且不再存储无用点，故将其存储结构设置
为：首先按顺序对单元节点进行存储，存至ID数
组中；然后按照单元编号顺序分组储存其相邻节
点并从0开始重新储存在First数组中；之后将相邻
节点本身在ID中的编号按照在数组First中的顺序
存储在Second数组中；最后设置一个数组IDID用
以存相邻节点分组的起始编号。

通过这种数组映
射的存储方式，实现缩短正演程序耗时的目标。

为解决正演过程中涉及大型矩阵求解等问题，
减少计算时间、提升计算效率，本文将引入并行
算法以此实现对正演过程的加速运算。

OpenMP是
基于CPU多核多线程的共享式内存结构并行计算，
具有对源程序改变较小等优点，因而被广泛应用
641
第 55 卷
中南大学学报(自然科学版)于数据处理计算中[24]。

若将完成一个过程中所有元素的求解称为“任务”，则求解其中一个元素称为“子任务”，当对每个任务来说其中每个子任务相对独立，各个子任务的求解之间互不影响时，可采用OpenMP 算法将子任务均分至多个线程同时计算，实现对这一过程的求解并行，在同一个程序中可嵌套多个并行区间。

由式(12)，令−K ′U 0=B ，则可得方程：
KU =B
(14)将对称正定矩阵K 进行分解[25]：
K =LL T
(15)
这里可引入中间向量Y ，求解方程(14)等价于以下方程组：
L T U =Y (16)LY =B
(17)
在求解下三角矩阵L 的元素时，首先求对角线元素L jj ，再利用L jj 求解L ij (j <i )，求得的矩阵L 的元素为[25]：
L ij =ì
íî
ïïïïïï
()
K ij -∑p =i j -1L ip L jp /L jj j <i (
)
K ii -∑p =i i -1L 2
ip 1
2 j =i (18)对于方程组(16)和(17)，首先求解式(17)，可
得到中间向量Y ，将Y 代入式(16)，可得到异常电位向量U 。

由于对向量Y 和U 元素逐个求解耗时较多，因此当对一个求解过程来说其中每个元素相对独立且互不影响时，采用OpenMP 算法将任务均分至多个线程同时计算，在同一个程序中嵌套多个并行区间，实现对于不同电源点供电的求解并行。

本文选择15个线程同时运算，在行数小于15
时按照所需线程并行运算，详细计算流程如图3所示。

其中，在求解U 的过程中，若供电方式为单点供电，则B 为列向量。

因本文采用多同性源阵列供电，故右端向量B 为矩阵：
B =--K ′U 0=--K ′[U 01 U 0i U 0n ]
(19)
式中，U 0i 为第i 个电源点供电时的正常电位向量，共n 个供电点。

3 全时域激发极化超前探测数值
算例
为了研究基于滤波算法的全时域激发极化正演算法对含水目标体的响应特征，并验证该方法
的可行性，通过数值模拟分别测试了隧道掌子面
图3　多点源并行计算流程示意图
Fig. 3　Multi-point source parallel computing process
diagram
图 2　单元网格示意图[23]Fig.2　Unit grid diagram [23]
642
第 2 期李志强，等：基于滤波算法的隧道全时域激发极化超前探测正演方法
前方不同距离异常体的响应以及不同供电极距探测方法对异常体的响应。

本文采用多同性源阵列激发极化观测模式[26]，供电电流为1 A ，从上到下布置3条测线，依次为测线1、2、3。

多同性源阵列激发极化观测模式三维示意图如图4所示。

每条测线布置个9测点(测点1号~9号)，如图4(b)所示。

在隧道掌子面后方布置6圈供电电极，每圈设置4个供电点(A1~A4)，如图4(a)所示，每圈供电电极距离掌子面分别为0、5、15、30、50、75 m ，多同性源供电电极按照每圈设定位置依次向掌子面后方移动探测。

3.1　算例一
时间域激发极化各参数对异常体的反映规律尚不清晰，因此，有必要在隧道掌子面前方设置不同距离的低阻异常体模型，研究隧道全时域激发极 化正演响应特征，设隧道掌子面长×宽为10 m× 10 m ，低阻体前界与掌子面距离d 分别为4、8、12、16、20、24、28 m ，异常体的长×宽×高为 10 m×10 m×6 m ，模型示意图如图5所示。

选用六面体为网格剖分单元，采用有限单元法进行正演
计算。

异常体和围岩参数设置如表1所示。

采用基于滤波算法求解等效电阻率后，为验证全时域激发极化正演加速算法的效果，还进行了无并行加速算法的隧道全时域激发极化正演计算，采集了30个时刻的正演数据，即将整个正演过程分为30个时窗。

采用正演加速算法将计算时长由30 min 减少至3 min ，实现了隧道全时域激发极化的快速正演。

分别采集了充电过程中30个时刻的电位，由Cole-Cole 模型作为数据的拟合模型，采用多元非
线性最小二乘拟合，计算出视零频电阻率、视固有极化率、视弛豫时间等激发极化参数。

最后，绘制了不同距离异常体的各参数正演结果曲线，如图6所示。

由图6可见：全时域激发极化各参数随异常体位置改变而变化，其曲线形态的变化趋势存在明显特征。

视零频电阻率随着异常体与掌子面距离的增大而增大，数值趋向于985 Ω·m ，逐渐接近围岩的零频电阻率；视固有极化率随异常体距离掌子面距离的增大而减小，数值趋向于0.016 5，且逐渐接近围岩的极化率；视弛豫时间随异常体距离掌子面距离的增大而减小，并且数值趋向于0.3，即逐渐接近围岩的弛豫时间；视频率相关系数随异常体距离掌子面距离的增大而减小，数值趋向于0.6
，且逐渐接近围岩的频率相关系数。

全时域
(a) 观测方式三维示意图；(b) 掌子面阵列电极布置示意图
图4　多同性源电极激发极化超前探测示意图Fig. 4　Schematic diagrams of advanced prediction with
induced polarization method
图5　正演计算模型示意图
Fig. 5　Forward calculation model diagram
表1　数值模型中使用的参数
Table 1　Parameters used in numerical model 介质围岩异常体
频率相关系数c 0.50.8
弛豫时间τ/s 0.32
固有极化率η00.0050.2
零频电阻率ρ0/(Ω·m)
1 000100
643
第 55 卷
中南大学学报(自然科学版)激发极化参数曲线变化特征反映了掌子面前方异常体距离变化产生的影响，表明了该正演方法可反映掌子面前方含水目标体的存在。

3.2　算例二
算例二研究了不同供电极距下全时域激发极化对含水目标体的正演响应特征。

异常体前界面与掌子面距离d 分别为4 m 和28 m ，其他参数设置与算例一保持一致，采集30个时刻的正演数据后，通过多元非线性最小二乘拟合，得到了各参数正演结果曲线，分别如图7和图8所示。

由图7和图8可见：当异常体距离掌子面距离d =4 m 时，随着供电测量极距的增大，视零频电阻率逐渐增大，视固有极化率、视弛豫时间和视频率相关系数则逐渐减小；当异常体距离掌子面距离d =28 m 时，随着供电测量极距的增大，视零频电阻率先下降而后增大，当供电测量极距为15 m 时，视零频电阻率正演结果出现极小值；视固有极化率、视弛豫时间和视频率相关系数则先增大而后减小，且当供电测量极距为15 m 时，这3个参数均出现极大值。

通过对比图7和图8可以看出：当掌子面前方存在含水体时，视零频电阻率
减小，视频率相关系数、视弛豫时间、视固有极化率数值增大，且视零频电阻率、视固有极化率、视弛豫时间和视频率相关系数均在同样的供电测量电极距上同步产生异常，因此，全时域激发极化4个参数对于含水目标体的影响较为一致。

4 全时域激发极化超前探测物理
模型试验
为了进一步验证全时域激发极化超前探测正演方法的有效性，并初步探讨全时域激发极化超前探测的实际应用效果，开展了隧道全时域激发极化法超前探测的物理模型试验。

本模型试验是在已建的隧道超前预报物理模型上开展的[27]，物理模型试验示意图如图9所示。

本模型试验的几何相似比为1꞉5，模型的长×宽×高为17.0 m×8.4 m×6.7 m ，隧道腔体采用玻璃钢材料制作，隧道腔体长×宽为2.0 m×1.7 m 。

含水目标体中充填一定比例的泥沙和水，控制其电阻率为10 Ω·m
，含水目标
(a) 视零频电阻率曲线；(b) 视固有极化率曲线；(c) 视弛豫时间曲线；(d) 视频率相关系数曲线
图6　距掌子面不同距离异常体的激发极化各参数的正演结果
Fig.6　Forward modeling results of induced polarization parameters of anomalous bodies at different distances from the
tunnel face
644
第 2 期李志强，等：基于滤波算法的隧道全时域激发极化超前探测正演方法
体的长×宽×高为1 m×3 m×3 m ，位于掌子面前方5 m 处。

围岩采用已配置好的相似材料，分层压实，
对试验区进行填充，电阻率为856 Ω·m 。

采用多同性源阵列超前探测的观测方式，在掌子面布置2
条
(a) 视零频电阻率曲线；(b) 视固有极化率曲线；(c) 视弛豫时间曲线；(d) 视频率相关系数曲线
图7　异常体距离掌子面4 m 时的正演结果
Fig.7　
Forward result with distance between abnormal body and working face of 4 m
(a) 视零频电阻率曲线；(b) 视固有极化率曲线；(c) 视弛豫时间曲线；(d) 视频率相关系数曲线
图8　异常体距离掌子面28 m 时的正演结果
Fig. 8　Forward result with distance between abnormal body and working face of 28 m
645
第 55 卷
中南大学学报(自然科学版)间距为0.4 m 的测线，每条测线上共布置10个测量电极组成阵列测量电极系(M1~M10，M11~M20)。

在测量电极系的外围，掌子面的周边布置由4个供电电极组成的供电电极系(A1~A4)，每圈供电电极
距掌子面分别为0、0.8、1.8、2.6、3.8 m ，共5圈。

通过试验，测得激发极化4个参数的结果如图10所示。

通过分析实测数据曲线发现全时域激发极化
(a) 物理模型示意图；(b) 物理模型试验观测方式；(c) 掌子面阵列电极布置示意图
图9　物理模型试验示意图Fig. 9　
Diagram of Physical model test
(a) 零频电阻率曲线；(b) 固有极化率曲线；(c) 弛豫时间曲线；(d) 频率相关系数曲线
图10　各参数实测结果
Fig. 10　Measured results of parameter s
646
第 2 期李志强，等：基于滤波算法的隧道全时域激发极化超前探测正演方法
各参数随供电测量极距的变化亦存在明显特征。

由图10可知：随着供电测量极距的增大，零频电阻率在d=2.6 m时出现极小值，固有极化率、弛豫时间在供电测量极距d=2.6 m时出现极大值，这一现象可以反映掌子面前方含水目标体的存在，为定性判断异常体位置提供有效依据。

通过与数值算例结果进行对比分析发现，数值模拟中视弛豫时间、视频率相关系数均出现极大值，全时域激发极化4个参数对目标体的响应特征与数值正演结果类似，这也验证了全时域激发极化正演方法的有效性，同时初步验证了隧道全时域激发极化超前探测的可行性。

5 结论
1) 本文建立了基于滤波算法的隧道全时域激发极化正演方法，利用滤波算法进行时间域激发极化中等效电阻率的计算，有效解决了Cole-Cole 模型在计算过程中收敛性较差的问题。

引入OpenMP并行算法，实现了隧道全时域激发极化加速正演。

2) 隧道全时域激发极化各参数曲线呈现极值响应特征，且均在相同位置附近产生异常，表明了基于滤波算法的隧道全时域激发极化对于含水构造等目标体具有一致性反应。

3) 全时域激发极化超前探测对目标体的响应结果与数值模拟结果类似，这为后续隧道全时域激发极化超前探测工程现场试验奠定了基础。

参考文献：
[1]钱七虎. 隧道工程建设地质预报及信息化技术的主要进展
及发展方向[J]. 隧道建设, 2017, 37(3): 251−263.
QIAN Qihu. Main developments and directions of geological prediction and information technology of tunnel construction
[J]. Tunnel Construction, 2017, 37(3): 251−263.
[2]洪开荣. 近2年我国隧道及地下工程发展与思考(2017—
2018年)[J]. 隧道建设(中英文), 2019, 39(5): 710−723.
HONG Kairong. Development and thinking of tunnels and underground engineering in China in recent 2 years (From 2017 to 2018)[J]. Tunnel Construction, 2019, 39(5): 710−723.
[3]田四明, 王伟, 巩江峰. 中国铁路隧道发展与展望(含截至
2020年底中国铁路隧道统计数据)[J]. 隧道建设(中英文),
2021, 41(2): 308−325.
TIAN Siming, WANG Wei, GONG Jiangfeng. Development and prospect of railway tunnels in China(including statistics of railway tunnels in China by the end of 2020)[J]. Tunnel Construction, 2021, 41(2): 308−325.
[4]OTTO R, BUTTON E, BRETTEREBNER H, et al. The
application of TRT-true reflection tomography at the Unterwald tunnel[J]. Felsbau, 2002, 20(2): 51−56.
[5]ALIMORADI A, MORADZADEH A, NADERI R, et al.
Prediction of geological hazardous zones in front of a tunnel face using TSP-203 and artificial neural networks[J].
Tunnelling and Underground Space Technology, 2008, 23(6): 711−717.
[6]李术才, 刘斌, 李树忱, 等. 基于激发极化法的隧道含水地
质构造超前探测研究[J]. 岩石力学与工程学报, 2011, 30(7): 1297−1309.
LI Shucai, LIU Bin, LI Shuchen, et al. Study of advanced detection for tunnel water-bearing geological structures with induced polarization method[J]. Chinese Journal of Rock Mechanics and Engineering, 2011, 30(7): 1297−1309.
[7]CHIRINDJA F J, DAHLIN T, JUIZO D, et al. Reconstructing
the formation of a costal aquifer in Nampula Province, Mozambique, from ERT and IP methods for water prospection [J]. Environmental Earth Sciences, 2016, 76(1): 36.
[8]LI Shucai, NIE Lichao, LIU Bin. The practice of forward
prospecting of adverse geology applied to hard rock TBM tunnel construction: the case of the Songhua River water conveyance project in the middle of Jilin Province[J].
Engineering, 2018, 4(1): 281−294.
[9]刘春生, 武强, 曾一凡. 基于钻孔激发极化法的井下超前探
水技术研究[J]. 煤炭工程, 2018, 50(3): 99−103.
LIU Chunsheng, WU Qiang, ZENG Yifan. Study on advance underground water forecasting based on borehole induced polarization method[J]. Coal Engineering, 2018, 50(3): 99−103.
[10]凌丹丹, 薛胜利, 周少伟. 激电测深在城市浅层地下水
调查评价中的应用[J]. 地球物理学进展, 2019, 34(5): 2079−2086.
LING Dandan, XUE Shengli, ZHOU Shaowei. Application of induced polarization sounding method to exploration shallow groundwater in the city[J]. Progress in Geophysics, 2019, 34(5): 2079−2086.
[11]LESMES D P, FRYE K M. Influence of pore fluid chemistry
on the complex conductivity and induced polarization responses of Berea sandstone[J]. Journal of Geophysical Research Solid Earth, 2001, 106(B3): 4079−4090.
[12]KEMNA A, BINLEY A, CASSIANI G, et al. An overview
of the spectral induced polarization method for near-surface
647。