2017, Vol. 32
2. 云南省建设投资控股集团有限公司, 昆明 650000
3. 浙江水利水电学院, 杭州 310018
2. Yunnan Construction Investment holding group Co., LTD, Kunming 650000, China
3. Zhejiang University of Water Resources and Electric Power, Hangzhou 310018, China
随着城市地下工程和水利工程建设事业的日益发展，出现了大量作为围护和防渗功能的地下结构连续墙，其施工可能且难以保证不出现缺陷隐患.以软土地区为例，由于基坑深度的愈来愈大，作为地下空间围护结构的地下连续墙的施工也愈加困难，会由于施工技术或者周边地质条件的影响，使得地连墙接缝出现错位或地连墙渗漏的问题，若不能及时发现并加以处理，会造成地连墙坍塌，危及周边环境的安全，造成人民财产损失.然而，由于地下连续墙宽度小且缺陷隐患深度往往又较深，在连续墙顶部即开展地面物探难以满足工程需要.
针对此上述难题，本文借助地下连续墙的测斜孔，开展跨孔电阻率层析成像，对地下连续墙病害进行探测研究.
1987年，日本学者Shima首次提出电阻率层析成像(ERT)这一概念，同时给出了相应的反演解释.跨孔电阻率成像勘探方式是在两个测孔中放置一定数量的电极，其中一部分电极发射电流进行供电，另一部分电极进行电位观测，利用观测电位值转换的视电阻率值进行反演成像计算，最终得到两测孔之间区域电阻率的分布图像.目前电阻率成像的反演以线性反演为主(Daniels and Dyck, 1984; Shima，1992; Bing and Greenhalgh, 2000)，反演方法本身也相当成熟.但线性化反演方法对初始模型的依赖性很强，且实际工程中的反演问题是复杂而非线性，如果初值选取不当，反演结果很可能会陷入局部极值，从而得不到全局最优解.相比较而言，非线性反演方法不依赖于初始模型，具有较强的全局搜索和非线性映射能力.研究表明，非线性反演比线性反演得到的结果更加接近实际(吴小平和汪彤彤，2002；杨文采, 2002；王家映，2007；杨磊等，2016).
在电阻率层析成像的非线性反演中，Martínez等(2010)将粒子群优化算法引入到电阻率反演，并取得了不错的结果；张倩等(2015)利用充分混沌振荡粒子群优化方法提升了成像质量；徐海浪和吴小平(2006)使用人工神经网络进行了一维和二维的电阻率反演研究；卢元林等(1999)将模拟退火算法运用到二维的电法资料反演；Sen等(1993)实现了结合模拟退火的贝叶斯一维电阻率测深反演；Yang和Torres-Verdín(2011)利用多重网格和滑动时窗的方法优化改进了测井电阻率层析成像的贝叶斯反演；Bouchedda等(2015)将Matérn协方差矩阵引入到电阻率层析成像的贝叶斯反演当中，且取得不错结果.但是，目前的这些随机反演的结果在精度和成像质量上有所欠缺，且反演时间过长，工作效率较低.在实际工程探测中，葛双成等(2006, 2013)都做了许多的研究工作，并取得了不错的效果.
在本文所提出并实现的二维跨孔电阻率反演成像策略中，电阻率的计算借鉴Pidlisecky和Knight(2008)提出的正演算法，运用奇异/边界校正算法，消除了近电源处的奇异性提高了剖分网格边界的精度；再将Xu等(1998)提出的最优化波数运用到正演算法中，以提高正演数值模拟精度.在反演算法中，基于贝叶斯反演框架，使用快速傅里叶变换移动平均算法简化先验信息约束初始模型生成的过程，结合扩展的Metropolis算法使用连续吉布斯采样生成后验概率密度，以达到提高反演精度和成像质量的结果.
1 方法原理假设地下空间的离散模型参数为m，通过电极观测的数据为d，地下模型参数通过某种适当的物理映射关系影响着观测数据的大小.因此数据和模型参数的正演关系可以表示为
|
(1) |
式中，g是线性或非线性的映射算子，其通常遵循着某种物理定律.本文中，模型参数表示的是地下地电空间的电阻率或电导率，观测数据通常为电位值或电位差值.本文使用有限单元法(FEM)，结合最优化波数和奇异/边界校正来求解式(1) 的正演问题，数值模拟结果如图 1所示.
|
图 1 正演数值模拟结果曲线图 (a)电位曲线图；(b)相对误差图. Figure 1 Forward modeling numerical simulation results (a) Potential curve; (b)Relative error. |
从图 1中上可以看出，采用了最优化波数和奇异/边界校正进行数值模拟使得近电源处的数值模拟精度上升，且在边界处的误差远远低于传统算法(徐世浙，1988)得到的结果，这也佐证了使用该正演算法的正确性和可行性.
1.1 贝叶斯反演理论贝叶斯理论是由Thomas Bayes在1763年提出来的，他首次将归纳推理方法应用于概率统计理论中，并由此创立了贝叶斯理论.在贝叶斯反演框架中，反演目标不再仅仅寻找最优解，而是以概率论为基础，使用模型参数的均值、方差和概率密度函数(Probability Density Function，PDF)，来描述反演结果不确定的分析，通常用后验概率密度来描述反演结果，其表达式为
|
(2) |
式中，c正则化常数，ρM(m)表示先验概率密度，L(m)为似然函数.ρM(m)描述的是多种先验信息对模型参数的约束.L(m)表示在考虑数据噪声和模型不确定性的统计描述下，模拟数据和观测数据的拟合程度.所以，σM(m)描述的是数值模拟数据给定的信息与模型参数给定的先验信息相互调整的结果状态.因此，充足的数据不确定信息和较好的先验信息是获取准确后验信息的关键.
1.2 采样后验概率密度在高度非线性反演问题中，由于模型参数与观测数据的非线性正演关系，假设先验概率密度为高斯分布，使用快速傅里叶变换移动平均算法(Ravalec et al., 2000)进行先验概率的生成，加之似然函数也假设为高斯分布.根据公式(2) 可知，后验概率密度为两个高斯分布的概率密度的积也是高斯分布的.
在本文中，选用扩展的Metropolis算法(Mosegaard and Tarantola, 1995)对后验概率密度进行采样，该算法非常方便，因为在采样过程中不需要先验概率密度精确的表达式.扩展的Metropolis算法主要由以下两个步骤构成：
1) 探索阶段：随机生成一个当前模型mcurrent，然后随机扰动生成一个新的候选模型mpropose，与此同时生成了先验概率密度.
2) 开发阶段：决定候选模型是否被接受.候选模型满足Metropolis接受概率则(参照Metropolis准则)被接受.公式为
|
(3) |
式中，L(mcurrent)/L(mpropose)为当前模型的似然函数与候选模型的似然函数的比值.如果候选模型被接受了，则候选模型变为当前模型，相反，则候选模型被舍去，当前模型继续代入上述步骤中再次计算.因此，每次迭代过程中，后验概率密度的样本量在增加.
在探索阶段中，候选模型从先验概率密度中随机生成.与较大探索步长相比，较小的探索步长，候选模型与当前模型的相关度较高.因此，根据公式(3)，小的探索步长会导致较高的接收概率，反之亦然.由于每算一次似然函数就要进行一次正演计算，这会增加计算成本，所以对加快探索后验概率密度来说，选择一个合适的步长是非常重要的.但在非线性反演中要想减少太多的由于计算似然函数产生的计算量是非常困难的，所以探索策略在处理多维先验信息时就变得愈加重要.
1.3 先验信息采样由于传统的采样策略有时会非常的低效率，Hansen等(2012)提出了一种灵活的、可变采样步长的采样算法(连续吉布斯法)来对高维先验信息进行采样.因此，本文利用连续吉布斯采样来对先验概率密度进行采样，即连续吉布斯采样算法充当扩展的Metropolis算法中的“黑箱子”算法，用来对先验概率密度P(m)进行采样.连续吉布斯采样算法的流程如下：
1) 定义一个二维矩形网格下的模型物性参数(一般来说为均值和方差)，基于先验概率密度随机生成一个无约束的初始的模型mcurrent.
2) 在当前模型mcurrent中，边长为查找步长Estep的正方形区域内，物性随机扰动生成一个与mcurrent相关的新的一个模型msubarea.
3) 利用连续模拟获取条件概率密度
4) 候选模型mpropose模型就变成步骤2) 中的当前模型mcurrent，重复步骤2) 和步骤3).
重复上述步骤，把连续吉布斯采样运用到扩展的Metropolis算法中，当候选模型mpropose被Metropolis准则剔除时，当前模型mcurrent要重复用来计算.此外，当前模型和候选模型之间的关联性是由查找步长(也就是重复模拟方形区域的边长)来控制的.
1.4 似然函数由式(2) 可知，似然函数L(m)表达的是模型参数m与观测数据d之间的相对概率，一般用它的大小来反映模型响应与观测数据的失配程度.一般来说，实际观测数据会存在一些噪声，反演结果也会有相当的误差.假设这些噪声和误差的分布为高斯分布，则似然函数表达为
|
(4) |
式中，g(m)和dobs分别表示电流激发的模拟数据和实测数据.CD为正则化常数，它是定义数据不确定度的方差和协方差的协方差矩阵，它包含了观测数据和模拟误差，若不知道每个模型的观测误差和模拟误差的大小，可选取为单位矩阵.
1.5 贝叶斯反演算法流程本文采用基于贝叶斯反演算法为基础的反演理论，其关键是在反演过程中进行大量的抽样模型样本，以便计算后验概率密度和进行样本的概率统计来获得模型均值和其对应的方差.贝叶斯反演方法与传统的反演方法差别主要在产生新模型后，采用的贝叶斯方法需要判断新模型是否是全模型空间中的一个样本，“是”则该模型被储存，“否”则舍弃，而这个判断标准是新模型应该满足Metropolis接受准则.基于贝叶斯反演方法的二维跨孔电阻率成像演整体上可以分为三大部分：先验信息的收集与转化、反演计算与样本模型收集、模型参数统计，具体流程如图 2所示.
|
图 2 贝叶斯反演算法流程图 Figure 2 Flow chart of Bayesian inversion algorithm |
在上述的反演理论基础上，首先考虑建立简单模型异常体模型来验证反演的可靠性.跨孔电阻率的观测方式有多种，比如单极-单极、单极-偶极(偶极-单极)和偶极-偶极排列方式，本文采用稳定性较好、数据较为丰富的AM-BN偶极-偶极装置(Bing and Greenhalgh, 2000).建立如图 3所示模型，在均匀半空间中存在一个电阻率2 Ω·m为的圆形低阻体(模拟地下渗流通道)，赋存环境电阻率为10 Ω·m.模型中异常体位于探测的两个孔位之间，为半径1 m的圆形，圆心位于地下5 m处，电极间距分别为0.5 m、1 m和1.5 m，选取的孔间距为3 m，反演迭代为20000次.
|
图 3 单个圆形低阻异常体示意图 Figure 3 Diagram of a single circular low resistivity anomaly |
由图 4a、图 4b、图 4c可看出，当反演迭代次数足够多，最终迭代解、最大后验概率解和均值解均能很好地刻画出模型中的低阻异常体，而且均值解的方差在异常体空间位置处的值是很小的，说明了模型参数(电阻率)在该区域的扰动是微小的和稳定的，综上反演结果验证了本文所使用的反演算法的可靠性.由于贝叶斯反演是基于概率统计原理和Metropolis准则，当采样次数足够多的时候，接受的当前样本较之前样本的扰动会很小，这就会使得反演结果中的最终迭代解、最大后验概率解和均值解相差不大.
|
图 4 单个圆形低阻异常体反演结果 (a)最终迭代解；(b)最大后验概率解；(c)均值解；(d)均值解的方差. Figure 4 Inversion results of single circular low resistivity anomaly (a) The final iteration solution; (b)Maximum posterior probability solution; (c)Mean solution; (d)Variance of mean solution. |
接下来考虑跨孔电阻率层析成像对条带状异常体的识别特点，按照图 5建立了条带状异常体三种不同的空间位置模型进行数值模拟，图中存在一低阻条带状异常体，电阻率为2 Ω·m，周围赋存环境电阻率为10 Ω·m，水平与竖直条带状异常体尺寸为1 m×5 m，中心位置位于地下4.5 m处，倾斜条带状异常体倾斜角度为18°，厚度为0.36 m，长度为6 m，中心位置位于地下4.5 m处，网格剖分大小为0.1 m×0.1 m，电极间距分别为0.5 m、1 m和1.5 m，水平条带状模型选取的孔间距为6 m，其他两个模型孔间距为3 m，采用AM-BN阵列方式进行数据采集，反演迭代为20000次.
|
图 5 不同空间位置条带状模型示意图 (a)水平条带状模型；(b)竖直条带状模型；(c)倾斜条带状模型. Figure 5 Patterns of stripes with different spatial distribution (a)Horizontal strip pattern; (b)Vertical strip pattern; (c)Oblique strip pattern. |
从图 6中可以看出，跨孔电阻率层析成像对水平条带状异常体水平方向的分辨效果较为理想，而在竖直方向的刻画并不是很精确，这个情况很类似在井间地震勘测中的情况，井间地震勘测在竖直方向的分辨率不高，解决的办法一般是在地面加检波器，增加竖直方向的数据，跨孔电阻率层析成像可以借鉴该思路，扩展一下跨孔电阻率勘测，在地表加一系列的电极，扩展为井地电阻率勘测手段，这样可以到达提高地下空间探测的的分辨率.
|
图 6 水平条带状异常体反演结果 (a)最终迭代解；(b)最大后验概率解；(c)均值解；(d)均值解的方差. Figure 6 Inversion results of horizontal strip (a) The final iteration solution; (b)Maximum posterior probability solution; (c)Mean solution; (d)Variance of mean solution. |
从图 7中可看出，竖直条带状异常体反演的结果很好，对异常体位置、形态刻画比较准确，只是在异常区域周边会出现一片较高电阻的假象区域，这样虽然低阻异常体可以勘测出来，但在实际勘测中可能会对数据的解释带来一定的扰乱.
|
图 7 竖直条带状异常体反演结果 (a)最终迭代解；(b)最大后验概率解；(c)均值解；(d)均值解的方差. Figure 7 Inversion results of vertical strip (a) The final iteration solution; (b)Maximum posterior probability solution; (c)Mean solution; (d)Variance of mean solution. |
从图 8看出，对于倾斜的条带状异常体的成像效果一般，只能看到有一条很细小的低阻带区域，但在其周边会伴生出现相对应的两个高阻带，这与原模型的背景场值相差甚远，是由于异常体是倾斜状态，会与电流线产生一定程度上不规则的干扰，所以反演结果会产生一定程度的畸变.
|
图 8 倾斜条带状异常体反演结果 (a)最终迭代解；(b)最大后验概率解；(c)均值解；(d)均值解的方差. Figure 8 Inversion results of Oblique strip (a) The final iteration solution; (b)Maximum posterior probability solution; (c)Mean solution; (d)Variance of mean solution. |
不同地连墙内部隐患的理论计算结果表明：
(1) 基于贝叶斯反演方法的二维跨孔电阻率成像方法，使用了扩展的Metropolis算法，避免求取准确的先验概率密度所带来的计算耗时，使用快速傅里叶变换移动平均算法简化了先验信息约束初始模型的过程.
(2) 在传统采样算法的基础上，使用了自适应可变采样步长的连续吉布斯采样算法，使得能够更精确的逼近全空间的最优解，避免陷入局部极值的问题，通过与理论数值解的对比验证了此方法的有效性和准确性.
(3) 基于贝叶斯反演的二维跨孔电阻率成像结果可得到最大后验概率解、均值解还有评价均值解的方差，可以更加充分地认识反演结果的不确定性，丰富了反演结果评价，最终的反演结果可以准确的反演地下地点空间的物性特征，说明了该算法是稳定的、可信的.为后期实际工程中的地下连续墙隐患探测奠定了理论基础.
3.2在贝叶斯反演框架中，影响最终反演结果的因素有许多，比如先验信息对初始模型的约束影响着反演的收敛速度，跨孔间距和电极间距的大小对异常体的分辨率也有一定的影响.另外，由于模拟数据量很大，如果使用全部数据进行反演会花费较长的时间，如果考虑将跨孔雷达层析成像中的“覆盖角度”的思想引入到对数据的筛选中，提取对异常体较为敏感的数据进行反演，可有望提高反演速度及准确度，以上因素具有重要研究价值.这些问题，将在后续工作中展开进一步研究.
致谢 感谢审稿专家提出的修改意见和编辑部的大力支持！| [] | Bing Z, Greenhalgh S A. 2000. Cross-hole resistivity tomography using different electrode configurations[J]. Geophysical Prospecting, 48(5): 887–912. DOI:10.1046/j.1365-2478.2000.00220.x |
| [] | Bing Z, Greenhalgh S A. 2001. Finite element three-dimensional direct current resistivity modelling:Accuracy and efficiency considerations[J]. Geophysical Journal International, 145(3): 679–688. DOI:10.1046/j.0956-540x.2001.01412.x |
| [] | Bouchedda A, Giroux B, Gloaguen E. 2015. Bayesian ERT inversion using non-stationary inverse Matérn covariance matrix[C].//SEG Technical Program. Expanded Abstracts, 2246-2251, doi:10.1190/segam2015-5884510.1. |
| [] | Daniels J J, Dyck A V. 1984. Borehole resistivity and electromagnetic methods applied to mineral exploration[J]. IEEE Transactions on Geoscience and Remote Sensing, GE-2(1): 80–87. DOI:10.1109/TGRS.1984.350582 |
| [] | Ge S C, Liang G Q, Lu Y X, et al. 2006. The study and test of geophysical exploration technique applied to harmless detection in seepage control wall[J]. Journal of Geotechnical Investigation & Surveying(10): 72–76. |
| [] | Ge S C, Liu X L, Zhao Y H, et al. 2013. Application research of detecting technology for seawall hidden defects and grouting effect[J]. Chinese Journal of Underground Space and Engineering, 9(1): 42–47. |
| [] | Hansen T M, Cordua K S, Mosegaard K. 2012. Inverse problems with non-trivial priors:Efficient solution through sequential Gibbs sampling[J]. Computational Geosciences, 16(3): 593–611. DOI:10.1007/s10596-011-9271-1 |
| [] | Lu Y L, Wang X T, Wang R, et al. 1999. The simulated annealing method of electrical resistivity tomography[J]. Chinese Journal of Geophysics, 42(S): 225–233. |
| [] | Martínez J L F, Gonzalo E G, Álvarez J P F, et al. 2010. PSO:A powerful algorithm to solve geophysical inverse problems:Application to a 1D-DC resistivity case[J]. Journal of Applied Geophysics, 71(1): 13–25. DOI:10.1016/j.jappgeo.2010.02.001 |
| [] | Mosegaard K, Tarantola A. 1995. Monte Carlo sampling of solutions to inverse problems[J]. Journal of Geophysical Research, 100(B7): 12431–12447. DOI:10.1029/94JB03097 |
| [] | Pidlisecky A, Knight R. 2008. FW2_5D:A MATLAB 2.5-D electrical resistivity modeling code[J]. Computers & Geosciences, 34(12): 1645–1654. DOI:10.1016/j.cageo.2008.04.001 |
| [] | Ravalec M L, Noetinger B, Hu L Y. 2000. The FFT moving average (FFT-MA) generator:An efficient numerical method for generating and conditioning gaussian simulations[J]. Mathematical Geology, 32(6): 701–723. DOI:10.1023/A:1007542406333 |
| [] | Sen M K, Bhattacharya B B, Stoffa P L. 1993. Nonlinear inversion of resistivity sounding data[J]. Geophysics, 58(4): 496–507. DOI:10.1190/1.1443432 |
| [] | Shima H. 1992. 2-D and 3-D resistivity image reconstruction using crosshole data[J]. Geophysics, 57(10): 1270–1281. DOI:10.1190/1.1443195 |
| [] | Shima H, Sakayama T. 1987. Resistivity tomography:An approach to 2-D resistivity inverse problems[C].//Expanded Abstracts of 57th SEG Annual Meeting. New Orleans:Society of Exploration Geophysicists, 59-61, doi:10.1190/1.1892038. |
| [] | Wang J Y. 2007. Lecture on non-linear inverse methods in geophysics (一):Introduction to geophysical inverse problems[J]. Chinese Journal of Engineering Geophysics, 4(1): 1–3. DOI:10.1088/1742-2132/4/1/001 |
| [] | Wu X P, Wang T T. 2002. Progress of study on 3D resistivity inversion methods[J]. Progress in Geophysics, 17(1): 156–162. DOI:10.3969/j.issn.1004-2903.2002.01.023 |
| [] | Xu H L, Wu X P. 2006. 2-D resistivity inversion using the neural network method[J]. Chinese Journal of Geophysics, 49(2): 584–589. DOI:10.3321/j.issn:0001-5733.2006.02.035 |
| [] | Xu S Z. 1988. Selection of wave number k in Fourier inverse transformation for point source and 2-D electric field problems[J]. Computing Techniques for Geophysical and Geochemical Exploration, 10(3): 235–239. |
| [] | Xu S Z, Duan B C, Zhang D H. 1998. Selection of the wavenumbers k using an optimization method for the inverse Fourier transform in 2.5D electrical modelling[J]. Geophysical Prospecting, 48(5): 789–796. DOI:10.1046/j.1365-2478.2000.00210.x |
| [] | Yang L, Zhang Z Y, Li M, et al. 2016. 2D joint inversion of direct current resistivity and magnetotelluric sounding data[J]. Progress in Geophysics, 31(2): 851–855. DOI:10.6038/pg20160247 |
| [] | Yang Q S, Torres-Verdín C. 2011. Efficient 2D Bayesian inversion of borehole resistivity measurements[C].//SEG Technical Program. Expanded Abstracts, 427-431, doi:10.1190/1.3628101. |
| [] | Yang W C. 2002. Non-linear geophysical inversion methods:Review and perspective[J]. Progress in Geophysics, 17(2): 255–261. |
| [] | Zhang Q, Wang L, Jiang F B. 2015. 2-D Improved particle swarm optimization algorithm for electrical resistance tomography inversion[J]. Geophysical and Geochemical Exploration, 39(5): 1047–1052. DOI:10.11720/wtyht.2015.5.27 |
| [] | 葛双成, 梁国钱, 陆云祥, 等. 2006. 物探技术在防渗墙质量无损检测中的应用试验研究[J]. 工程勘察(10): 72–76. |
| [] | 葛双成, 刘喜亮, 赵永辉, 等. 2013. 海塘隐患及防渗注浆效果检测技术应用研究[J]. 地下空间与工程学报, 9(1): 42–47. |
| [] | 卢元林, 王兴泰, 王若, 等. 1999. 电阻率成像反演中的模拟退火方法[J]. 地球物理学报, 42(S): 225–233. |
| [] | 王家映. 2007. 地球物理资料非线性反演方法讲座(一)地球物理反演问题概述[J]. 工程地球物理学报, 4(1): 1–3. |
| [] | 吴小平, 汪彤彤. 2002. 电阻率三维反演方法研究进展[J]. 地球物理学进展, 17(1): 156–162. DOI:10.3969/j.issn.1004-2903.2002.01.023 |
| [] | 徐海浪, 吴小平. 2006. 电阻率二维神经网络反演[J]. 地球物理学报, 49(2): 584–589. DOI:10.3321/j.issn:0001-5733.2006.02.035 |
| [] | 徐世浙. 1988. 点电源二维电场问题中付氏反变换的波数k的选择[J]. 物探化探计算技术, 10(3): 235–239. |
| [] | 杨磊, 张志勇, 李曼, 等. 2016. 直流电阻率法与大地电磁法的二维联合反演[J]. 地球物理学进展, 31(2): 851–855. DOI:10.6038/pg20160247 |
| [] | 杨文采. 2002. 非线性地球物理反演方法:回顾与展望[J]. 地球物理学进展, 17(2): 255–261. |
| [] | 张倩, 王玲, 江沸菠. 2015. 电阻率层析成像的二维改进粒子群优化算法反演[J]. 物探与化探, 39(5): 1047–1052. DOI:10.11720/wtyht.2015.5.27 |