1 引言

近年来，电阻率探测技术日益成熟和完善，其应用领域和服务对象不断扩展，促使电阻率反演方法研究得到了长足发展(Zohdy, 1989; Li and Oldenburg, 1992; Sasaki, 1994; Zhang et al., 1995; 毛先进等, 2000; 吴小平和徐果明, 2000；底青云和王妙月, 2001)，主要分为线性方法和非线性方法.对于三维电阻率探测而言，线性反演方法是最主要的手段，三维电阻率线性反演方法较为成熟，局部搜索能力强、收敛速度快.在线性反演中，正则化、先验约束等被用来压制多解性和降低假异常，如：光滑约束(阮百尧等, 1999; 黄俊革等, 2004; 宛新林等, 2005)、不等式约束(Kim et al., 1999; 刘斌等, 2012a, b)等，起到了积极作用.由于场的等效性、外界干扰大、初始模型选取不当等原因，多解性问题有时仍较严重，特别是在初始模型选取不当的情况下，容易导致陷入局部极小，使得线性反演结果与实际情况相差较大，出现假异常、分辨率低等问题，影响探测效果.近年来，随着非线性反演方法的发展，人们开始研究电阻率非线性反演方法，如神经网络法(Stephen et al., 2004; Singh et al., 2005, 2010; 徐海浪和吴小平, 2006)、模拟退火法(Dittmer and Szymanski, 1995; 卢元林等, 1999)、遗传算法(王兴泰等, 1996; Jha et al., 2008; 李帝铨等, 2008; 闫永利等, 2009；Liu et al., 2012)等.从已有研究来看，非线性反演适合于反演参数较少的一维和二维反演情况，因此，目前非线性反演方法基本都集中在一维或二维反演方面.非线性反演方法的全局搜索能力强，对初始模型的依赖程度很低.但是对于同样的问题，非线性方法的耗时是线性方法的十几倍甚至几十倍，搜索效率低、耗时长是非线性方法不可回避的缺点.同时由于三维反演中反演参数数量多达成百上千，远远多于一维和二维情况，使得非线性反演方法的全局搜索的优化方向难以控制，随机性和不确定性大大增加，容易导致出现早熟的问题，多解性问题有时也很严重.总体而言，目前尚未见到关于三维电阻率非线性反演的文献.

通过上述分析可见三维电阻率反演主要存在着以下几方面问题：(1)从方法原理和研究实践中来讲，电阻率的线性反演方法与非线性反演方法之间存在着较强的互补特性，如何将二者相互结合来解决各自的弊端是一个重要问题；(2)多解性问题是电阻率反演的固有难题，如何较好的利用先验信息来尽量压制多解性问题值得思考；(3)对于非线性方法而言，三维电阻率反演参数数量较大，非线性反演方法全局搜索方向难以优化控制，搜索效率极低，如何合理控制非线性反演搜索方向也是需要思考的.

针对上述问题，本文将介质电阻率变化范围作为不等式约束引入到三维电阻率反演中，定义了较宽松不等式约束和局部较严格不等式约束两种类型，分别提出了基于不等式约束的最小二乘线性反演方法和遗传算法反演方法(属非线性方法)，发现线性和非线性反演方法具有很强的互补性，融合二者的优势对于提高三维电阻率探测效果是十分有益的.因此，本文提出了混合反演方法的概念，将基于不等式约束的改进线性方法与改进非线性方法相互嵌套，旨在压制多解性、摆脱初始模型依赖、提高反演效果.最后结合合成数据与实际数据反演算例，验证了混合反演方法的可行性与有效性.

2 基于不等式约束的非线性与线性反演方法 2.1 电阻率反演中的不等式约束先验信息获取及定义

在反演中引入已知信息，对于降低多解性有直接作用.在传统的电阻率反演中，光滑约束经常被用来改善反演方程的病态程度，使得相邻网格间“光滑过渡”.但是，光滑约束是一种较为宽松的约束，基于介质物性连续变化的“天然特性”来对反演区域施加约束.虽然光滑约束改善了反演方程的病态程度，但对多解性的改善作用往往不理想，反演结果中往往会出现较多假异常，与实际情况相差较大.实际上，在电阻率探测中，有一种很重要的已知信息未得到很好的利用，即介质电阻率变化的上下限范围.如果能够把这种先验信息施加到反演方程中，可以将反演中电阻率的搜索范围限定在一个封闭范围内，缩小了反演求解的可行空间，可起到压制多解性的效果.很显然，我们可以用不等式约束来表征这种已知的电阻率变化范围信息.

在实际的探测工作中，这种介质电阻率的变化范围可以通过如下方式获取和定义：

(1)半开放式的不等式约束：根据电阻率的非负特性，将电阻率的变化下限定义为0，这点已经在最小二乘反演中通过反演“对数化”实现了，这是一种最简单的不等式约束.

(2)较宽松的不等式约束：根据探测区域岩性和含水性等地质情况，依据常识推断区域电阻率的变化范围.或者根据视电阻率数据及其分布特征，基于正演经验，大致推断反演区域电阻率的变化范围，当然，通过这两种方式获取的不等式约束都是较为“宽松”的.而且由于这种方式难以对探测区域的局部获取较详细的上下限范围，也就使得这种不等式约束通常是对整个反演区域施加一个相同的较宽松的不等式约束.

(3)局部较严格的不等式约束：在探测区域进行钻探，并测量岩芯电阻率，进而推断出更加确切的介质电阻率变化范围，并通过多个钻孔取芯测试结果推断，可推断出较严格不等式约束在反演区域内的分布.从而在反演中施加上这种局部严格不等式约束，有助于较大程度的改善多解性问题.

可见，按照上述不等式约束的定义，其对多解性问题的改善作用由强到弱排序应是：局部严格不等式约束、较宽松不等式约束、半开放不等式约束.实质上，在实际应用中往往是上述三种不等式约束的综合应用.

2.2 携带不等式约束的三维电阻率改进遗传算法非线性反演

在本节中三维电阻率非线性反演选择遗传算法(Genetic Algorithm，GA).传统的遗传算法是无法用于三维电阻率反演的，必须对其进行改进和创新.三维电阻率反演是一种大型反演问题(反演参数多达成千甚至上万)，对于这种大型反问题，传统遗传算法无法适用的主要原因是遗传过程中变异方向的选取是随机的，“变好”与“变差”的概率是相同的，而不是沿着最优的方向进化，不利于保护优良个体，导致遗传搜索效率和速度低下，反演不稳定或反复震荡，更不利于改善多解性和探测可靠性及稳定性，甚至会导致反演失败.针对上述问题，本文提出了一种变异方向控制方法和参数范围不等式约束，并对传统遗传算法进行了一些必要的改进.

假设N为观测数据的数量，M为模型参数的数目，假定N_p为遗传群体的大小(即群体中包含的个体的数量)，m_i(i=1, 2, 3, …N_p)所表示的模型参数向量为群体中的第i个个体.

2.2.1 基于不等式约束的三维电阻率反演目标函数

本文将不等式约束施加到三维电阻率遗传算法反演的目标函数中，式(1)中包含了表征模型参数变化范围的不等式约束，使得搜索空间缩小，有利于减小多解性.公式(1)为

(1)

式中，Φ_i为对应于第i个个体的目标函数值，d为实际观测数据，d_mi表示对应于模型参数m_i的理论观测数据，C为光滑度矩阵，m_ik为第i个个体中的第k个参数，ρmin_k和ρmax_k分别为第k个模型参数取值范围的下限与上限，λ为光滑约束的权重系数，经多次验证，认为λ的合理取值在0.02~0.1，本文取λ=0.05.

对于第i个网格而言，光滑约束可表示为

(2)

式中，m_i为m_i群体中第i个网格的模型参数修正量，m_i^F, m_i^B, m_i^L, m_i^R, m_i^U, m_i^D分别表示第i个网格的前、后、左、右、上、下的相邻网格的电阻率修正量.将整个模型的光滑约束用矩阵形式表示为

(3)

式中，C为光滑度矩阵.

目标函数的计算是以数值正演为基础的，通过数值计算可求得理论观测数据d_mi.本文采用三维有限单元法进行正演模拟，利用Cholesky分解法求解三维有限元大型线性方程组，得到各节点的电场响应，进而得到d_mi.

2.2.2 初始群体的近似均布产生

传统遗传算法采用随机方法来产生初始群体，由于其随机性导致产生相似个体的几率增加，或者个体分布不均匀.不均匀的个体分布使得初始群体接近“最优解”的几率大大降低，或者接近“最优解”的机会不稳定，使得后续反演搜索效率降低甚至是不收敛.对此，有学者提出了严格均布产生的方式，可有效的控制初始群体的多样性和均匀性.但这种严格均布产生方法不适用于反演参数较多的情况.对于反演参数较多的三维电阻率反演问题，严格均布产生方式需要足够多(大约上百个)的个体才能得到接近最优解的初始个体，导致计算量剧增，现有PC机或计算工作站无法承担如此大量的计算任务.本文则采用了一种用海明距离控制的近似均布产生方式(刘斌, 2010)，主要流程是：(1)首先，通过分析观测数据和经验判断，在给定的电阻率变化范围内(即不等式约束)给出第一个初始个体；(2)然后，在给定的电阻率变化范围内(即不等式约束)随机产生下一个个体，使得该个体与已产生的所有个体之间的距离都大于海明距离, 使得种群中的个体分布趋于均布化；(3)最后，将所有已产生的个体放入初始群体中，构成近似均匀分布的初始群体.

个体i与个体j之间的海明距离d_ij表示为

(4)

通过定义遗传种群中任意两个个体之间的海明距离并使之满足d_ij>d_min(d_min为最小体海明距离，一般略小于).

2.2.3 适应度值计算与选择操作

个体的优良程度与其适应度的大小呈正比关系，根据经典的遗传算法的适应度值计算方法，本文取个体i的适应度值F_i为：

(5)

式中，C_max为一个相对较大的实数.由于针对不同的实测数据和地电模型，目标函数的具体计算值是有很大区别的，因此，本文采用如下方法来确定C_max的取值：

(1)针对某一具体的算例，进行试算，得到一组个体对应的目标函数值Φ_i.

(2)找出上述Φ_i中的最大值，将C_max设定为最大值的10~20倍.

本文采用比例选择的策略和最优保存策略，以使得优良个体获得更多的遗传优势.

2.2.4 交叉操作与隐含的不等式约束引入

遗传算法中交叉和变异是产生下一代的主要手段，交叉主要体现由父代到子代的传承性，而变异主要体现由父代到子代的变化性.本文采用自适应交叉概率和算术交叉策略(王小平和曹立明, 2002; 姚磊华和李竞生, 2004)，算术交叉是将群体中的个体进行两两随机配对，将选中配对的两个个体做线性组合就可以产生对应的子代新个体.假设在个体m_i^(k)与m_j^(k)之间进行交叉运算，则有：

(6)

式中，m_i^(k+1)与m_j^(k+1)为新一代个体，μ为[0, 1]范围之内的一个常数.

实质上，在执行算术交叉操作的同时也实现了不等式约束的引入.分析式(4)可看出，假如a与b分别为m_i^(k)与m_j^(k)中的第n个网格的电阻率值，且a>b，那么下一代中第n个网格的电阻率值只能限定在[a, b]区间内，不会超越这个范围.上述规律具有普遍性，可推广到任意交叉操作中.因此，只要是初始群体产生是控制在给定的电阻率变化范围内的(本文第2.2.2节就采用了这个限定)，那么后续的遗传进化中的每一代的每个个体都不会超越给定的不等式约束.所以，在遗传算法中的算术交叉策略隐含的实现了不等式约束的引入，原理简单，易于理解.

2.2.5 变异操作

变异操作决定了遗传算法的局部搜索能力.从根本上讲，传统的遗传算法无法用于三维电阻率反演的主要原因就是其遗传搜索方向无法优化控制.在变异操作中常用的非均匀随机算术变异算法采用了随机变异策略，个体“变好”和“变差”的几率相同，也就使得变异方向无法控制，导致反演效率低下，甚至由于个体“变差”的程度过高或者“变差”的个数数量过多而导致反演失败.

最小二乘法局部搜索能力较强，其线性反演方程能够对下一代优化方向进行较好的控制和引导，这为遗传变异提供了很好借鉴.本文引入了最小二乘法，对变异方向进行控制(刘斌, 2010).对于某一参与变异操作的个体i，假设其变异增量为Δm，求解公式为

(7)

式中，Δd=d-d_m为实际观测数据与正演理论观测数据的差向量，Δm为模型参数的增量向量，A为偏导数矩阵，又称敏感度矩阵，表示模型的理论观测数据对模型参数的偏导数，C为光滑度矩阵，λ₂为拉格朗日常数，一般取值为0.1~0.5之间，本文取值0.25.

利用方程(5)求得可控变异增量Δm来产生下一代新个体，取代了传统遗传算法中变异增量随机扰动的方法，从理论上来讲，个体“变差”的可能性大大降低了，可有效的避免“优良基因”的丢失，对优化遗传方向和进化效率有积极作用.生物学和生物技术领域早已出现了“控制变异”的技术，即通过一定的生物技术手段，对基因变异进行控制，得到想要的生物新品种，本文借鉴这种技术思路，基于线性反演的方法对变异操作进行控制，优化变异方向和幅度，这应是对遗传算法的一种有益的探索.

2.3 携带不等式约束的三维电阻率最小二乘线性反演

对于线性反演方法，如何将这种不等式约束施加到线性反演方程中是一个难题.以下反演目标函数同时包含了光滑约束和不等式约束，公式为

(8)

式中，d_obs为实际观测数据，d_m为正演得到的理论观测数据，m为模型参数，m_i为第i个网格的电阻率，ρmin_i和ρmax_i分别为第i个网格的电阻率的下限和上限.

刘斌等(2012a)采用引入相关参数向量的方法来引入不等式约束，起到了较好效果.但是这种不等式约束的施加算法对初始模型依赖较强，若模型参数初值选取不当，反演过程中模型参数接近约束边界时，导致新的模型参数出现模型电阻率值在上下边界两端来回振荡的情况，无法收敛，难以搜索到全局最优解.针对上述问题，本文采用障碍函数法将不等式约束的信息嵌入到目标函数中，在式(8)的基础上构造增广目标函数，公式为(刘斌等, 2012b; Li and Oldenburg, 2000)：

(9)

式中，μ恒大于零，为障碍因子.

求解上述增广目标函数(9)的最优化问题，得到公式为(刘斌等, 2012b)：

(10)

式中，A为偏导数矩阵，Δm为模型参数增量向量，Δd为观测数据，Δd=d_obs-d_m，e=(1, 1, …, 1)^T，X和Y都是对角矩阵，X矩阵的对角线元素为m_i-ρmin_i(i=1, 2, M)，Y矩阵的对角线元素为ρmax_i-m_i(i=1, 2, M).需要说明的是，A为n×m阶矩阵(其中m为参与反演的网格数量，n为观测数据的数量)，表示模型的理论观测数据对模型参数的偏导数，设A_ij为矩阵A的元素，则

(11)

对偏导数矩阵的计算，本文采用常用的电极互换定律法，阮百尧等介绍了用互换定理来计算偏导数矩阵方法(阮百尧, 2001)，黄俊革提出了在异常电位法情况下利用互换定理求取偏导数矩阵的方法(黄俊革, 2003)，具体原理本文不再赘述.

线性方程组(8)是含有不等式约束的三维电阻率线性反演方程，在线性反演方程中已经嵌套了不等式约束，这种正则化处理包含了重要的先验信息，使得反演参数变化被限定在较小的封闭范围内，有利于压制反演多解性、提高反演可靠性.

上述线性反演收敛的标准是理论观测数据与实际观测数据之间的均方差满足收敛判据rus < 为反演收敛的容许值).

3 三维电阻率探测混合反演方法

本文提出了线性方法与非线性方法相互嵌套的混合反演方法，融合了最小二乘法反演方法和改进遗传算法反演方法(其流程如图 1，刘斌, 2010; 刘斌等, 2012c)，主要如下：

(1)按照2.1节提出的不等式约束的确定方法，通过对观测数据分析、现场岩样电阻率测试、钻孔岩芯测试等方式，给出不等式约束的上下限范围.

(2)在全局反演阶段(即混合反演的第一阶段)采用全局搜索能力较强的改进遗传算法，它同时具有对初始模型的依赖程度低的优势，为了控制总体耗时，一般第一阶段的改进遗传算法的遗传进化代数上限t_max为4~7次.待达到规定的遗传代数之后(即t < t_max，其中t是遗传代数，t_max是遗传代数上限)，输出当前最优个体，作为最小二乘反演的初始模型.在混合反演中，改进遗传算法的任务是通过少数几次遗传进化，搜索到最优解附近，为后续的线性反演提供一个较优的初始模型.

(3)在局部搜索阶段(即混合反演的第二个阶段)，采用局部搜索能力较强的最小二乘法，将遗传算法得到的最优个体作为初始模型，在接近最优解的局部空间内执行高效率搜索和反演，直到满足收敛标准后输出反演结果，最后对反演结果进行三维可视化，并对反演结果的质量进行评价.由于最小二乘法计算效率较高，不对其反演次数设置上限.

可见，混合反演算法既发挥了非线性方法不依赖初始模型的优势，又融合了线性方法局部搜索能力强的特点，以其兼顾全局搜索与局部搜索的优势，在摆脱初始模型依赖、降低多解性、改善反演质量等方面起到作用.

4 合成算例验证

为了验证本文提出的携带不等式约束的三维电阻率混合反演方法，设计了一个合成算例.三维地电模型见图 2，该模型由两层电阻率不同的介质构成，上一层为高阻(1000 Ωm)，下一层为相对低阻(800 Ωm)，在下一层中赋存一个倾斜低阻带(50 Ωm).设计了两条平行测线，采用施伦贝谢尔装置形式，每条测线的电极数为35个.

首先建立三维有限元模型，其反演目标区采用等边长八节点六面体单元剖分，在边界区域采用逐渐增加的八节点六面体单元，网格数量为8976个，其中反演目标区域(见图 2a)的网格数量为2856个(注：本算例中的反演模型都采用上述同一模型).采用三维有限单元法计算得到合成观测数据，见图 3.为更贴近实际情况，合成观测数据中含有3%的随机噪声.

4.1 基于光滑约束的传统最小二乘反演结果

为了对比携带不等式约束的混合反演方法的有效性，首先采用传统的最小二乘法对合成数据进行反演.通过分析合成视电阻率剖面图，推断探测区域的背景电阻率在1000 Ωm左右，因此将最小二乘反演的初始模型设定为均一模型，初始电阻率设定为1000 Ωm.此次反演共迭代6次，耗时约25 min (计算机配置：酷睿四代i7处理器，速度3.6 GHz，下同).反演结果见图 4，反演结果中大体反映了倾斜低阻带的存在，但存在着较大问题：

(1)反演结果中在浅部出现了两处高阻区，推断是浅部高阻层的反映.但两处高阻区彼此独立，未连接成层状，其最高电阻率幅值为2571.24 Ωm，与原型中电阻率幅值(1000 Ωm)相差较大.

(2)蓝色低阻区是原型中倾斜低阻带的反映，但其形态近似团状，与原型的条带形状相差较大.其电阻率最低为217.49 Ωm，与原型中的电阻率(50 Ωm)也有较大差异.

总体而言，传统最小二乘反演结果与原型在异常体形态、位置和电阻率幅值方面有很大区别，给地质解释带来困难.

4.2 携带较宽松不等式约束的混合反演结果

按照本文中对不等式约束的定义，首先对混合反演施加较为宽松的不等式约束.对合成视电阻率剖面(图 3)进行分析，发现视电阻率剖面中底部存在低阻区，推断探测区域内存在低阻异常，所以不等式约束的下限定义的要“宽松”一些；而视电阻率剖面中除了在浅部存在层状高阻区外，没有发现单个高阻体的反映，所以不等式约束的上限可以稍严格一些，比视电阻率的最高值稍高即可.综合分析，本算例给出了一个较为宽松的不等式约束上下限：0~1600 Ωm.

在定义好不等式约束的基础上，利用混合反演方法对合成数据进行反演，反演过程中rus随反演代数的变化情况如图 5，其中改进遗传算法中群体数量设为20个，共进化5代，耗时约210 min，最小二乘法反演阶段迭代5代，耗时约20 min.

由图 5可见，整个反演过程中rus由初始的136.22 Ωm降至4.56 Ωm.在遗传算法阶段(前五代)，rus由136.22 Ωm降至43.2 Ωm，未出现波动震荡，呈现较快的单调递减的趋势，表明其搜索方向和遗传效率得到了优化，第3~5次遗传操作中rus的变化比较平稳，出现了一个平缓台阶，表明非线性反演呈现收敛趋势，应采取线性反演来进一步深化反演；而最小二乘反演阶段(后五代) rus由43.26 Ωm降至4.21 Ωm，rus在平稳台阶之后出现了又一个较大降幅，反演进一步收敛，表明线性反演方法局部搜索能力较强.

基于上述不等式约束的混合反演结果如图 6.同时，为了对比验证，将上述不等式约束施加到了传统最小二乘反演方程中，对于合成观测数据单独利用方程(7)开展了不等式约束最小二乘反演，反演结果见图 7.

对比分析图 4、图 6和图 7，可见：

(1)由于引入了不等式约束，最小二乘反演结果(图 7)比传统最小二乘的反演结果(图 4)有了一定改善，浅部高阻体连接呈短条带状，浅部高阻最高电阻率为1567.33 Ωm，与原型浅部高阻层电阻率值较接近；蓝色低阻区域虽仍呈团状，但其形态趋于扁平，呈现出的倾斜角度与原型在趋势上较一致，其“宽度”有所减小.这表明不等式约束对于提高最小二乘线性反演的反演质量、降低多解性方面有积极作用.但总体上评价，携带不等式约束的最小二乘反演结果仍与原型有较大差异，蓝色低阻区的电阻率幅值与原型有较大不同(最低电阻率值为189.51 Ωm，而原型中为50 Ωm).

(2)对最小二乘法反演结果相比(图 4和图 7)，基于不等式约束的混合反演结果(图 6)有较大改进，浅部高阻层出现，虽然高阻层在z方向上“起伏变化”，左右两侧的高阻层形态、厚度均与原型有一定差别，但能够较确切的反映高阻层的存在，比最小二乘法有较大改善.同时，蓝色低阻区的形态呈近似“条带状”，“宽度”进一步减小，其倾斜程度与原型也较为一致，最低电阻率幅值为113.69 Ωm.从实际工程的探测需求来衡量，图 6所示反演结果基本满足实际要求.

可见，虽然采用了相同的不等式约束，混合反演结果要明显优于最小二乘反演结果.由于混合反演将遗传算法作为第一阶段反演方法，其初始群体采用了近似均布产生技术，对摆脱初始模型依赖、搜索全局最优有显著作用.可见，混合反演方法以牺牲计算时间为代价(是最小二乘法耗时的约9倍)，显著提高了异常体的成像效果.对于电阻率探测而言，探测效果是人们最为关心的因素，以牺牲计算时间来改善探测效果，是被认可的，具有重要的实际意义.

另外，需要说明的是，由于将遗传算法作为第一阶段反演方法，其初始群体采用了近似均布产生技术，对摆脱初始模型依赖有显著作用.笔者针对同一算例，开展过多次不同近似均布初始模型(通过控制个体之间的最小海明距离来产生不同初始模型)条件下的反演验证，反演结果基本一致，较为稳定，也就是说最终反演结果对于初始模型不依赖，这一点要明显优于最小二乘方法.

4.3 携带局部较严格不等式约束的混合反演结果

按照本文中对不等式约束的定义，本小节模拟实际探测工作中经常采用的钻探方法来获取更加严格的不等式约束信息，钻孔分布见图 8a.显然，通过对岩芯进行电阻率测试，可得到沿钻孔深度的介质导电性分布，将4个钻孔的测试结果进行综合对比，可得到较为严格的不等式约束信息的分布，见图 8b.考虑贴近工程实际情况，为使得不等式约束推断过程较为合理，本小节提出了由钻孔测试结果推断不等式约束的原则：

(1)对于几个钻孔出现相同或类似的电阻率分布的情况，将它们的电阻率分布做“横向连接”和推断解释是合理的，按照图 8a给出的钻孔分布，4个钻孔中均可揭露浅部高阻层，可连接起来做统一推断，经分析，将该区域的不等式约束上下限定为800~1200 Ωm；而钻孔2和钻孔3同时揭露了倾斜低阻带，也可连接起来做统一推断，将该区域的不等式约束范围定为0~200 Ωm.

(2)对于仅在单个钻孔中出现的电阻率结构，不宜做横向连接推断，可以做“局部延伸推断”，例如由于倾斜低阻带的存在，4个钻孔对于第二层相对低阻层的揭露很难做横向推断，因此对这种情况可将钻孔的已知信息向钻孔左右各延伸推断1个网格，在图 8b中给出了4个独立区域，这四个独立区域的不等式约束范围均定为600~1000 Ωm.

(3)对于未能通过钻孔揭露推断的剩余区域，则经过综合分析钻孔结果，给出一个较为宽松的不等式约束范围，在本节中将剩余区域的不等式约束范围定为0~1200 Ωm.

实质上，上述不等式约束是“局部较严格不等式约束”与“较宽松不等式约束”的集成.本小节中混合反演中初始群体的产生是在上述不等式约束范围内以近似均布的方式产生的，图 9是在上述较严格不等式约束下的混合反演结果，可见图 9所示反演结果比宽松约束下的反演结果(图 6)有进一步改善：(1)浅部高阻层起伏程度降低，平均厚度与原型接近(图 9中推断出了浅部高阻层与下部相对低阻层的界面)，电阻率在900~1200 Ωm之间；(2)倾斜低阻带的倾斜程度、“宽度”、电阻率幅值等与原型较为接近，其“宽度”进一步缩小，与原型的宽度基本一致，电阻率幅值在21.33~184.16 Ωm之间(最低为21.33 Ωm)，接近原型中倾斜低阻带的电阻率(50 Ωm).

从实际探测工作的技术要求来分析，图 9所示反演结果的质量是较高的，满足实际工程的需要.可见由于采用了局部较严格不等式约束，使得关键区域的反演效果大大提高，多解性问题得到了较好的控制，可为后续地质解释和推断提供可靠的基础.这表明，在实际电阻率探测中应充分利用可靠的先验信息，一些较确切和严格的约束信息将对反演结果的改善起到较大作用.

5 实际工程应用

南方某山区高速公路隧道施工中发生较大涌水，掌子面初始涌水量约800~1200 m³/h，而后隧道内出现了衬砌开裂、位移增加等问题.由于该段隧道埋深较浅(100 m左右)，隧道掌子面正上方出现塌坑和冒顶，同时，在隧道西侧约800 m处的村庄中相继出现地面塌陷、房屋开裂、水塘/泉水干涸、河水倒灌等问题，给当地居民的生产生活造成影响.图 10显示了整个测区地面塌坑、干涸泉水、隧道轴线位置等分布.结合区域地质、地形、涌水后的各种地质现象等综合分析，推断村庄至隧道涌水点之间应存在地下水导水通道.为探明是否存在导水通道及其深度位置，根据场地条件，在一条山路上开展了三维电阻率探测和钻探工作，测线见图 10a(刘征宇, 2014)，钻孔分布见图 10b，钻孔编号为ZK1-ZK5.沿着小路走向布置了三条平行测线, 测线间距15 m，测线长度为325 m，电极间距5 m，采用温纳观测装置.

在钻探工作中有两个钻孔(ZK2和ZK4)是在电阻率探测之前完成的，在实施电阻率反演前对ZK2和ZK4的岩芯进行了电阻率测试，可为电阻率反演提供可靠的已知信息.图 11是两个钻孔的柱状图及电阻率变化范围的推断结果.按照2.1节提出了不等式约束定义方法和4.3节提出的“由钻孔信息推断不等式约束”的方法，由于仅有两个钻孔信息，本案例中采用“较宽松不等式约束”+“局部较严格不等式约束”的方式来确定探测区域介质电阻率的变化范围.对于局部严格不等式约束，将图 11给出的钻孔已知信息向左右各延伸1个网格即可；对于剩余区域的宽松不等式约束，通过岩芯(样)测试确定为0~2700 Ωm.三维反演模型的单元数量为26000个，混合反演的初始群体在上述不等式约束范围内近似均布产生，混合反演中遗传算法共反演4代，最小二乘法反演7代，总耗时约为370 min.

基于三维混合反演方法的三维电阻率探测结果见图 12，结合已有两个钻孔揭露情况对三维电阻率图像及其二维切片进行分析发现：

(1)在左下部位存在着一个明显的低阻区域(电阻率在100 Ωm左右)，其具体位置是：x方向50~100 m之间，z方向80~100 m之间.由于该低阻区域在水平方向和深度方向均与隧道洞内涌水点接近，因此推断该低阻区域为地下水导通通道.

(2)在x方向200~230 m，z方向55~75 m之间的区域内也存在着一个低阻异常体，电阻率在500 Ωm左右，推断该区域是岩体局部破碎含水或充泥.

(3)同时，根据粘土层比基岩电阻率低的特性，推断出了土岩界面(见图 12b).

图 13是五个钻探的柱状图(包括电阻率探测完成后施工的3个钻孔)，其中ZK2和ZK3号钻孔在底部出现钻液、护壁泥浆全漏失的现象，揭露的岩体较破碎，局部有溶蚀痕迹，该部位与掌子面涌水位置接近，推断为地下水通道.将钻孔剖面图与三维反演结果切片对比，发现三维混合反演中推断导水通道的范围比钻孔揭露的通道范围稍大，但深度和位置基本一致.而根据混合反演结果推断出的土岩界面大体与钻孔揭露情况一致.另外，在反演结果中揭露了另一处低阻异常区，这个异常未在钻孔范围内，未能验证.

为了对比验证，图 14给出了传统最小二乘法(未携带不等式约束)的反演结果，与钻探结果存在着较多不一致：(1)反演结果中电阻率幅值最高达到7919.73 Ωm，与现场实测的介质电阻率范围明显不符；(2)反演结果中在钻孔ZK2和ZK3附近存在两个低阻异常区，位于左下角位置，推断应是导水通道的反映，但其位置、深度、范围均与钻孔揭露情况有较大差距，导致地下水导通通道这一关键目标体的判断和解释不准确；(3)反演结果中在钻孔ZK4和ZK5底部存在高阻异常区，而这与钻孔揭露情况不一致.

可见，本文提出的三维电阻率探测混合反演方法的可靠性和有效性在现场应用中得到了较好的验证，其探测精度亦满足工程需要，不足之处是反演耗时较长.高效率反演和数据处理是未来需重点解决的问题，特别对于抢险性的任务，对于反演计算速度有更高的要求.

6 结论

(1)本文针对三维电阻率反演中存在的问题，鉴于线性方法和非线性方法的互补优势，提出了二者相互嵌套的混合反演方法，从本质上来讲，混合反演方法将全局搜索和局部搜索相结合，利用近似均布初始群体产生方法使初始群体近似均匀分布在解空间上，利用非线性方法全局搜索能力强的特点，搜索到最优解的附近，再适时发挥线性方法局部搜索能强的特点，在最优解附近空间开展高效率的“小步长”搜索，有利于摆脱对初始模型的依赖，从理论上来讲，混合反演方法有其合理性和先进性.

(2)本文采用表征反演区域电阻率变化范围的不等式先验约束分别对最小二乘反演方法和遗传算法进行正则化处理，将原来只有非零约束的开放搜索空间改变成一个范围相对较小的搜索域，显然对压制多解性有积极作用.特别是在利用钻孔等手段获得较严格的不等式约束的情况下，反演效果和质量会大大提高，这一点在合成算例和实际工程案例中得到了验证.

(3)合成算例和工程实例中将携带不等式约束的混合反演算法分别与传统最小二乘算法、携带不等式约束的最小二乘算法的反演结果进行了对比分析，发现混合反演方法可较好的反映原模型的情况，其探测效果、定位精度、分辨能力、电阻率幅值均明显优于最小二乘法，可实现异常体的较精确三维成像.特别对于能够施加较严格不等式约束的情况，对多解性的改善作用十分明显.这表明，混合反演方法与不等式约束二者相结合，在提高反演质量、结果可靠性和工程实用性方面起到了较好作用.混合反演算法的主要问题是耗时要远高于线性方法，下一步的研究重点是采用并行计算策略，研究并行-混合反演方法，大幅缩短混合反演耗时，进一步提高其工程适应性.