2017, Vol. 60
2. 湖南涉外经济学院信息科学与工程学院, 长沙 410205;
3. 湖南师范大学物理与信息科学学院, 长沙 410081
2. Department of Information Science and Engineering, Hunan International Economics University, Changsha 410205, China;
3. College of Physics and Information Science, Hunan Normal University, Changsha 410081, China
可控源音频大地电磁法(Controlled source audio frequency magnetotellurics, CSAMT)和激发极化法(Induced Polarization，IP)均是地球物理勘探中应用较为广泛的电法勘探技术.CSAMT研究岩矿石间导电性质的差异，IP主要研究地下岩矿石间激发极化效应的差异.两者均有成熟的理论方法和实际应用.地电体对电磁波的响应为电磁感应和激电效应的综合反映，由于激电效应，地下可极化体的电阻率是一个与频率有关的复数.将两者结合进行研究，不但更加接近真实的地电情况，而且可以获得更多的物性参数，有助于提高频率域电磁法的反演与解释精度.
目前在频率域电磁法中提取激电信息的研究主要集中在MT和CSAMT领域，近年来取得了一些重要的进展：吴汉荣(1978)早在20世纪70年代就开展了对天然场源作激电法测量的可行性研究；Murali等(1980)认为利用低频天然电场作为频率域激电研究的场源是可行的；Gasperikova(2001)对天然场源的IP效应检测进行了深入的研究，并得到了肯定的结论，验证了天然场源IP提取的可行性；陈清礼等(2006)介绍了由MT资料反演真谱参数的基本原理，并结合Dias模型，给出了含激电效应的MT正演和基于梯度法反演的计算方法；曹中林等(2007)使用Cole-Cole模型模拟激发极化效应，采用广义逆方法对MT理论模型和实际资料进行反演，获得了较好的效果；岳安平等(2007, 2009a, 2009b)进行了含激电效应的CSAMT一维正演研究，分析了极化层位置和极化率变化对正演结果的影响；Ghorbani等(2009)利用Cole-Cole模型研究了含电磁效应的频谱域极化数据的正演和反演方法，其中反演采用的是高斯-牛顿法，同时还给出了相应正反演算法的开源软件；He等(2010)利用3D MT技术在柴达木盆地的三湖地区进行油气藏的勘探，并指出高电阻率和高极化率是寻找油气藏的两个关键指标；冯兵等(2013)分析了Cole-Cole模型参数对CSAMT信号的Ex、Hy的影响特征，采用最小二乘法反演电阻率、厚度和激电参数；Yu等(2013)采用均匀半空间的Cole-Cole模型, 研究了CSAMT电磁场各分量、视电阻率和相位与相关参数的关系，得到了充电率对正演结果的影响最大，频率相关系数次之，且时间指数影响最小的结论；徐汶东(2011)研究仿真了CSAMT中的IP效应影响，从参数变化的角度分析了Cole-Cole模型各参数对视电阻率幅值的影响水平；Tang等(2014)比较了Dias模型、Cole-Cole模型和Multi-Cole-Cole模型对MT激电响应信号的影响，分析了不同极化模型在含激电信息的MT信号中的适用情况.
上述研究表明从频率域电磁法信号中提取激电信息具有较为乐观的前景，在理论上切实可行，但是也存在一些局限性：(1) 考虑激电效应的正演研究能够分析激电参数对频率域电磁法勘探数据的影响，但无法提取激电参数本身；(2) 用反演来提取激电信息的研究多数采用传统的线性反演方法，其反演结果依赖于初始模型的选择，并容易陷入局部极值.
完全非线性反演方法为处理和解释电磁法资料提供了新的研究思路，受到地球物理学研究人员的广泛关注.Bäumer(1996)较早开展了使用遗传算法和模拟退火算法解释大地电磁资料的研究；Spichak(2000)利用BP神经网络进行了三维MT反演的研究，给出了神经网络建模时结构和参数的调整方法；Moorkamp等(2007)利用遗传算法进行了大地电磁资料和地震资料的联合反演；Shaw(2007)实现了粒子群优化算法对DC数据、IP数据和MT数据的反演，结果表明，粒子群优化算法与遗传算法在反演中效果相近，性能相当；李帝铨等(2008)采用遗传算法实现了CSAMT的最小构造反演；师学明等(2009)提出了一种阻尼粒子群优化算法，该方法采用惯性权重的振荡递减策略，能够快速有效地对大地电磁测深的理论模型和实测数据进行反演；胡祖志等(2015)利用人工鱼群算法对MT资料进行了最优化约束反演，验证了该算法的有效性；Jiang等(2016)利用神经网络对视电阻率数据进行了快速反演，并对神经网络的结构和训练过程进行了优化.
以上非线性反演方法虽然在电磁法资料的反演中已经有了较为广泛的应用，但是直接用在解释含激电信息的CSAMT资料时存在着以下不足：(1) 当考虑激电效应时，地下可极化体的电阻率是一个与频率有关的复数，在进行正演迭代时运算量有所增加，采用Monte Carlo类搜索算法如模拟退火SA、粒子群优化PSO和遗传算法GA，运算效率偏低；(2) 考虑激电效应后，反演模型参数在原有电阻率和厚度的基础上增加了激电参数，使得非线性反演过程变得更加难以收敛；(3) 由于地下导电性不均匀所引起的异常强度远大于激电效应引起的异常，在进行极化率搜索时极易陷入局部极值，加大了激电信息提取的难度.因此寻找一种高效、并能够兼顾全局和局部搜索的非线性反演方法，成为含激电信息的CSAMT资料IP提取的关键问题.
混合蛙跳算法(Shuffled frog leaping algorithm，SFLA)是一种启发式群体智能搜索算法，通过分组算子和模因搜索将全局信息交互和局部搜索有机地结合起来，同时在进化过程中只对群体中最劣个体的位置进行调整，相较于GA，PSO等需要对所有个体进行调整的群智能搜索算法具有更高的执行效率.因此，在CSAMT信号激电信息提取的应用背景下，我们提出了一种基于柯西分布和惯性权重的混合蛙跳二阶段最小构造激电信息提取算法.该算法在标准混合蛙跳算法非线性反演的基础上，进行了以下改进：(1) 采用柯西分布随机算子取代均匀分布随机算子来提高模因组内进化的种群多样性，增强算法的全局搜索能力；(2) 通过引入Logistic振荡惯性权重来均衡算法的前期搜索和后期收敛能力；(3) 通过引入二阶段最小构造反演加强后期极化率参数对观测数据的影响并抑制虚假构造，提高激电信息提取的准确性.同时，结合SFLA模因分组的特点，给出了本算法的多核CPU并行加速方案.理论模型的仿真结果表明，该方法反演结果稳定，能够较好地重构地电结构和提取激电信息，同其他非线性反演算法相比，具有优越的收敛速度和计算效率.
2 含激电信息的CSAMT正演理论 2.1 Dias模型Dias模型是应用较为广泛的激电响应模型之一，在1968年由Dias提出并改进，其参数的物理含义明确，对复杂实验数据的拟合精度高(Dias，2000).Dias模型的复电阻率计算公式为：
|
(1) |
式中
Cole-Cole模型是目前采用较多的另一种激电响应模型，在1941年由Cole兄弟提出，Pelton等(1978)对其进行了改进，其数学表达式比较简洁，能够较准确地描述大多数均匀岩矿石的激电效应.Cole-Cole模型的复电阻率计算公式为：
|
(2) |
式中c为频率相关系数.
2.3 含激电信息的CSAMT正演步骤含激电效应的CSAMT正演过程如下：(1) 设置地电模型的模型参数；(2) 给定激电响应模型(Dias模型或Cole-Cole模型)；(3) 计算与频率相关的复电阻率ρ(ω)取代地电模型参数中的实电阻率ρ；(4) 由各层的复电阻率，计算相邻层间的变换阻抗；(5) 由递推公式计算地面理论阻抗、场量和卡尼亚视电阻率，进而得到含激电响应的CSAMT正演结果(岳安平等, 2009a, 2009b).
3 改进的混合蛙跳算法 3.1 混合蛙跳算法的基本理论混合蛙跳算法数学模型描述如下：设蛙群规模为P，初始蛙群X(0)={X1, X2, …, XP}，可行解空间的维数为d，则解空间中第i只青蛙表示为Xi=[xi1, xi2, …, xid].混合蛙跳算法的初始蛙群随机产生，并均匀地分布于解空间中.生成初始蛙群之后，将蛙群内青蛙按照个体适应度降序排列，然后将整个蛙群分成m个模因组，每个模因组含n只蛙，满足关系P=m×n.具体的分组方式如下：第一只青蛙放入第一个模因组，第二只青蛙放入第二个模因组，直到第m只青蛙放入第m个模因组；然后第m+1只青蛙又重新放入第一个模因组，第m+2只青蛙放入第二个模因组，以此类推，直到所有青蛙分组完毕.具体的模因分组公式为：
|
(3) |
式中Mi为第i个模因组.每个模因组中具有最优和最差适应度的青蛙分别记为Xb和Xw，整个蛙群中的最优适应度青蛙记为Xg.蛙群分组后，分别对各个模因组内最差青蛙Xw进行局部搜索操作(局部搜索次数为iterl)，其更新公式为：
|
(4) |
|
(5) |
式中rand为[0, 1]间的独立随机数，t为迭代次数，Rmin和Rmax为青蛙允许移动的距离范围.经过更新后，如果得到的青蛙Xw(t+1) 适应度优于Xw(t)，则取代模因组内的Xw(t)，否则使用(6) 式取代式(4) 并重新执行局部搜索过程：
|
(6) |
如果仍然没有改进，则随机产生一个新解Xw(t+1) 取代Xw(t).当完成局部搜索后, 将所有模因组内的青蛙进行全局混洗，即重新混合排序.然后再次划分模因组，进行局部搜索，如此反复，直到达到定义的全局搜索次数iterg为止.
3.2 基于惯性权重-柯西分布的混合蛙跳算法对CSAMT信号中IP信息提取而言，标准混合蛙跳算法较好地均衡了算法的全局和局部搜索性能，有利于全局最优解的有效搜索以及整体执行效率的提高.但单一的随机算子rand使得各模因组内的局部搜索范围受限于Xb-Xw或Xg-Xw，在算法的后期易陷入局部极值；且青蛙的位置更新量Δw(t)未考虑自身的更新经验，不利于算法前期扩展搜索空间.因此本文对标准蛙跳算法的公式(4) 和公式(6) 进行改进，提出一种基于惯性权重-柯西分布的蛙跳算法(WCSFLA).具体如下：
|
(7) |
|
(8) |
式中w为惯性权重，Δw(t)为上一代的位置更新量，初始为0，Cauchy为柯西分布算子.其中Cauchy柯西分布算子的位置参数为对应的Xb-Xw或Xg-Xw，尺度参数为0.1.
由式(7) 和式(8) 可知，该算法的主要改进包括以下几个方面：(1) 采用柯西分布算子代替均匀分布算子，这样当组内的最差青蛙距离组内最优青蛙或全局最优青蛙较远时，意味着该青蛙需要更大的搜索空间，而柯西分布算子能够根据该距离生成更大的随机数，提高蛙群中个体的多样性；(2) 在位置更新中加入了最差青蛙自身的更新经验，即上一次的位置更新量Δw(t)，这样有利于蛙群在进化过程中引入个体经验，保持青蛙的搜索惯性，增强算法的全局搜索能力，但同时也可能影响算法后期的收敛；(3) 为解决算法的收敛问题，引入了粒子群算法中常用的惯性权重来均衡进化过程中的个体经验和群体经验，通常将w设定为线性递减，这样在算法的初期可以实现广泛的搜索，而在算法的后期则将搜索范围逐渐缩小，以实现稳定收敛，从而获得更好的算法性能.戴前伟(2013)提出了一种基于Logistic振荡的惯性权重模型，并验证了惯性权重采用振荡递减的下降形态能够帮助群体搜索算法加快收敛速度，跳出局部极值，优于传统的线性递减方法.因此本文采用该惯性权重模型来实现提出的惯性权重-柯西分布蛙跳算法，其计算公式为：
|
(9) |
|
(10) |
式(10) 为混沌Logistic方程，μ=4是控制参数；ws和we分别为惯性权重的初值和终值；t为迭代次数.
4 反演方法及流程 4.1 二阶段最小构造反演方法由于地下导电性不均匀所引起的异常强度远大于激电效应引起的异常(罗延钟，2003)，在Monte Carlo类非线性反演(如遗传算法、粒子群算法、差分进化算法、混合蛙跳算法)过程中，种群中个体电阻率参数对适应度函数的影响远大于极化率参数对适应度函数的影响.因此采用混合蛙跳算法进行激电信息的搜索时，蛙群的进化过程将明显地分为两个阶段：在第一阶段，由于电阻率参数对适应度函数曲线的主要影响，个体在解空间内将快速收敛至正确的电阻率参数附近；第二阶段，个体开始在电阻率参数附近进行微调，此时电阻率参数对适应度函数的影响趋于稳定，极化率参数的优化将成为适应度曲线下降的主要原因，但由于极化率参数在数值上远小于电阻率参数，同时对适应度函数的影响也远小于电阻率参数，算法在此时极易陷入局部极值，得到错误的极化率参数，从而加大了激电信息提取的难度.
针对以上激电信息提取过程中的反演特征，很多学者均认为应该首先对电阻率数据进行反演，然后以反演得到的电阻率数据为初始模型进行极化率参数的提取(He et al., 2010；罗卫锋等，2011；冯兵等，2013).
何继善(2010)指出CSAMT中电场分量Ex和磁场分量Hy中都含有介质的电阻率因素，只是不同的分量采用的公式互不相同，对电阻率变化的敏感程度有所差异.经过比较，电场分量Ex更加实用. 而冯兵等(2013)也分析了不同极化率参数对Ex和Hy的影响，并认为电场分量Ex相较磁场分量Hy对激电效应的反应更加灵敏.
在上述研究基础上，我们采用电场分量Ex作为误差项，给出一种应用于CSAMT信号激电信息提取的二阶段最小构造反演方法，该方法构造的适应度函数如下：
|
(11) |
式中R(ρ)和R(m)分别为对电阻率和极化率的最小构造约束函数；λ1、λ2分别为R(ρ)和R(m)对应的正则化因子，采用两个独立正则化因子的原因是极化率的取值空间(m∈[0, 1])较电阻率的取值空间有较大差异(一般可认为ρ≫m)，同时两者对适应度函数的影响程度也不一致，如果采用统一的正则化因子将无法有效约束相对较小的极化率参数；R(ρ)和R(m)在此均采用(12) 式进行计算(Constable et al., 1987)：
|
(12) |
式中，Miinv为反演得到的模型参数，包括各层的电阻率ρiinv以及极化率miinv.
E(e)为目标误差函数，在反演时为电场分量Ex的拟合误差，采用(13) 式进行计算.
|
(13) |
其中，Exobs表示观测数据矢量；Expre表示预测数据矢量，其计算公式如下：
|
(14) |
其中，F(·)为含激电信息的CSAMT正演算子；ρinv和minv分别为反演得到的模型电阻率参数和极化率参数；hinv为层厚，其初始化为对数等间距；R2为决定系数，是评价反演参数拟合程度的指标，其具体的定义参见5.1节；th为设定阈值.
由式(14) 可知，本文所提出的激电信息提取方法的反演过程分为两个阶段，在第一阶段中，电阻率参数ρ和极化率参数m同时参与反演，进行地电模型参数的搜索，此时由于ρ对观测数据的影响远大于m，混合蛙跳算法中的个体将快速地收敛于正确的ρ附近，R2的值迅速增加，当R2的值达到设定的阈值th时，电阻率参数的搜索结束，反演算法进入第二阶段；在第二阶段中，将搜索到的ρ和h值固定不变，反演算法仅针对个体中的m参数进行搜索，在此阶段中，一方面强化了极化率m对观测数据的影响，能够有效地克服算法后期的局部收敛问题，当达到设定的搜索代数时，算法将收敛于正确的m参数附近；另一方面由于搜索参数的减少，提高了反演算法在第二阶段的执行效率.
4.2 CPU并行加速模因组搜索采用Monte Carlo类搜索机制的非线性反演算法存在运算效率低，计算时间长的缺点.本文所提出的混合蛙跳反演算法只在初始化蛙群时需要调用正演算法计算每个青蛙的适应度，以后在进行组内模因搜索时仅需要对新生产的青蛙调用正演算法进行适应度评估，无需重复计算每个青蛙的适应度，提高了IP提取的效率.同时各模因组内的模因搜索是独立进行的，在全局混洗前不需要进行组间通信，因此便于利用Matlab的并行计算工具箱(Parallel Computing Toolbox)来进行并行优化.Matlab的并行计算工具箱可以使用多核处理器来解决计算问题和数据密集型问题，仅需要对Matlab程序的for循环进行改进就可以实现Matlab应用程序的并行化.本文的计算环境为：CPU为Core i7-4710HQ，内存为8 GB，操作系统为Windows 8.1.由于i7-4710HQ为4核CPU，因此可以通过本地运行多个本地的worker(独立的Matlab计算引擎)来充分利用多核处理器的计算资源.一般在配置本地Matlab并行计算池时，worker的数目应与多核CPU的核心数目一致，所以本文启动4个worker实现模因组搜索的并行计算.为了减少各worker之间的通信开销以及降低Matlab在协调和管理各worker时的管理开销，模因组的数目应该设置为4的倍数，本文经过凑试法设置为4，这样将有效提高混合蛙跳算法的并行计算效率.
4.3 反演步骤及流程本文提出了一种基于惯性权重-柯西分布的混合蛙跳算法用于激电信息的二阶段最小构造反演提取，并采用CPU并行计算来加速模因组的组内模因搜索过程.具体实现步骤如下：
(1) 初始化算法参数(迭代次数，种群规模，模因组数，惯性权重参数，柯西算子参数等)和反演参数(λ1、λ2和th等)；
(2) 根据地电模型参数ρ、m和h的搜索范围产生初始种群，开始第一阶段反演；
(3) 根据种群中青蛙的电阻率参数和极化率参数构造最小构造约束函数，同时调用含激电信息的CSAMT正演算法计算青蛙的误差函数E(e)和决定系数R2.由最小构造约束函数和误差函数共同评估种群中青蛙的适应度；
(4) 按照适应度进行排序，并更新全局最优青蛙；
(5) 根据决定系数与阈值的关系判断是否进入第二阶段反演；如果为否，则执行混合蛙跳算法的模因分组、并行模因组内搜索、全局混洗等操作，产生新一代种群，并转至步骤(3)；如果为是，则保持ρ和h参数不变，并转至步骤(6) 并进行第二阶段反演；
(6) 在当前最优解中m参数基础上进一步执行混合蛙跳算法的模因分组、并行模因组内搜索、全局混洗等操作，产生新一代种群，并更新全局最优解；
(7) 判断是否满足算法的终止条件，如果为否，则转至步骤(6)；如果为是，保存当前的全局最优解；
(8) 根据全局最优解输出并评估反演结果.反演算法的流程图如图 1.
|
图 1 反演算法流程图 Fig. 1 Flow chart of the inversion algorithm |
本文采用均方误差(Mean Square Error，MSE)、绝对百分比误差(Absolute Percent Error，APE)和决定系数(Determination coefficient，R2)来衡量算法的性能(Liu et al., 2011)，其相关定义如下：
|
(16) |
|
(17) |
式中yi为反演结果代入含激电信息的CSAMT正演计算后的第i个数据，Yi为对应的第i个观测数据，n为观测数据的数量.以上指标中MSE代表预测误差，APE代表误差的相对百分比，其值越小，表示反演模型的预测误差越小；R2代表预测值与观测值之间的相关度，其值越大，表示两者间存在着越明显的线性相关性.
5.2 三层K型地电模型实验以三层K型地电模型为例，研究极化层位于不同位置时，本文算法在不同阶段的反演性能.其中地电模型参数(Dias模型)、反演算法参数和二阶段最小构造参数的设置如表 1所示.
|
表 1 三层K型地电模型的仿真参数设置 Table 1 Parameter setting for K-type three-layer geoelectric model |
表 2给出了当极化层分别位于首层、中间层和底层的三种情况下，本文算法在不同阶段时的反演结果.由表 2可知，当极化层位于不同位置时，本文算法中第二阶段反演均能够获得较第一阶段反演更低的MSE、APE和更高的R2.这是由于在第二阶段中电阻率参数ρ和层厚参数h保持不变，反演性能的提升主要取决于极化率参数m，因此证明了二阶段最小构造反演在提取激电信息应用中的有效性.同时由表 2可知，在本文实验中，当首层极化时，反演结果的拟合性能最佳，当底层极化时，反演结果的拟合性能最差，但总的来说，当极化层位于K型模型的不同位置时，反演均能够获得较高的R2和较低的MSE、APE，说明本文所提出的基于惯性权重-柯西分布的混合蛙跳非线性反演能够有效地避免局部极值，提高算法的全局搜索能力.图 2进一步给出了算法在首层极化时模型的反演结果和Ex及Hy分量的拟合情况.
|
|
表 2 三层K型地电模型激电信息提取反演结果 Table 2 Inversion results of IP information extraction for K-type three-layer geoelectric model |
|
图 2 三层K型地电模型首层极化时反演结果 (a)电阻率；(b)极化率；(c) Ex拟合曲线；(d) Hy拟合曲线. Fig. 2 Inversion results of K-type three-layer geoelectric model with polarized upper layer (a) Resistivity; (b) Polarizability; (c) Ex fitting curve; (d) Hy fitting curve. |
以三层A型地电模型为例，研究不同极化模型下，本文算法的反演性能.其中地电模型参数、反演算法参数和二阶段最小构造参数的设置如表 3所示.
|
|
表 3 三层A型地电模型的仿真参数设置 Table 3 Parameter setting for A-type three-layer geoelectric model |
表 4给出了当极化层位于中间层时，本文算法采用Dias模型和Cole-Cole模型时的反演结果.由表 4可知，在本文实验中，Dias模型的反演结果优于Cole-Cole模型的反演结果.Tang等(2014)基于MT方法的研究中指出，采用Dias模型的MT正演曲线的激电响应主要体现在高频段，而采用Cole-Cole模型的MT正演曲线的激电响应主要体现在低频段.因此当极化层位于不同深度时，两种极化模型各有优势.图 3进一步给出了算法在Dias模型下的反演结果和Ex及Hy分量的拟合情况.
|
|
表 4 三层A型地电模型激电信息提取反演结果 Table 4 Inversion results of IP information extraction for A-type three-layer geoelectric model |
|
图 3 三层A型地电模型中间层极化时反演结果 (a)电阻率；(b)极化率；(c) Ex拟合曲线；(d) Hy拟合曲线. Fig. 3 Inversion results of A-type three-layer geoelectric model with polarized middle layer (a) Resistivity; (b) Polarizability; (c) Ex fitting curve; (d) Hy fitting curve. |
以三层H型地电模型为例，我们研究了分别采用电场分量Ex或视电阻率ρa作为适应度函数中的拟合误差时，算法对反演结果的影响.其中地电模型参数(Dias模型)、反演算法参数和二阶段最小构造参数的设置如表 5所示.
|
|
表 5 三层H型地电模型的仿真参数设置 Table 5 Parameter setting for H-type three-layer geoelectric model |
表 6给出了中间层极化时，本文算法适应度函数中拟合误差取电场分量Ex或视电阻率ρa时的反演结果.由表 6可知，首先，由于电场分量Ex和视电阻率ρa的单位和数量级均不同，所以MSE在比较时没有太大的参考价值；其次，适应度函数中拟合误差取电场分量Ex或视电阻率ρa均能够得到较好的反演结果.相较而言，采用电场分量Ex作为适应度函数中的拟合误差，反演得到的APE更小，R2更大.最后从计算时间可知，采用电场分量Ex作为适应度函数中的拟合误差，能够减少正演时求解Hy和ρa时的部分计算内容，计算效率更高，计算时间更少.图 4进一步给出了算法在Dias模型下的反演结果和Ex及Hy分量的拟合情况.
|
|
表 6 三层H型地电模型激电信息提取反演结果 Table 6 Inversion results of IP information extraction for H-type three-layer geoelectric model |
|
图 4 三层H型地电模型中间层极化时反演结果 (a)电阻率；(b)极化率；(c) Ex拟合曲线；(d) Hy拟合曲线. Fig. 4 Inversion results of H-type three-layer geoelectric model with polarized middle layer (a) Resistivity; (b) Polarizability; (c) Ex fitting curve; (d) Hy fitting curve. |
以三层Q型地电模型为例，研究本文算法与PSO算法(Eberhart，1995)、DE算法(Storn and Price, 1995)、标准SFLA算法(Eusuff et al., 2006)的反演性能.为便于比较，不同非线性反演算法的参数尽量保持一致，具体的地电模型参数(Dias模型)、反演算法参数和二阶段最小构造参数的设置如表 7所示.
|
|
表 7 三层Q型地电模型的仿真参数设置 Table 7 Parameter setting for Q-type three-layer geoelectric model |
由于非线性反演算法的计算时间主要来自于反复调用正演算法评估适应度的计算时间，所以为了较为公平地评估不同算法的计算性能，所有算法的种群规模均设置为20，总的迭代次数均设为600次(其中蛙跳算法的局部搜索次数为10，全局搜索次数为60，各模因组的总迭代次数为600次)，同时当适应度值连续100次迭代不变时认为算法收敛，停止进化.表 8给出了采用Dias模型，中间层极化时，不同非线性反演算法的性能和计算时间.由表 8可知，PSO算法和DE算法的性能相当，均优于标准的SFLA算法，这是由于激电信息提取的反演过程对非线性反演算法的全局搜索能力提出了更高的要求，而PSO算法和DE算法基于个体的进化在全局搜索性能上优于标准SFLA算法基于模因组的进化(因为PSO算法和DE算法的每个个体均参与了进化，而SFLA算法仅模因组内的最差个体参与了进化).而改进后的WCSFLA算法的反演性能，基本与PSO算法和DE算法相当.这是因为柯西随机算子和个体进化经验的引入提高了算法的搜索能力.由计算时间可知，SFLA算法和WCSFLA算法相较于PSO算法和DE算法最大的优势在于模因组内的局部搜索极大地提高了算法的收敛速度，同时最差青蛙的更新策略有效减少了算法的计算时间，这对于将非线性反演算法推广至实际应用具有非常重要的意义.通过采用Matlab的并行计算工具箱，我们改写了以上四种非线性反演算法的个体进化部分(即个体进化时采用并行计算)，并行计算时间表明，SFLA算法和WCSFLA算法的并行加速效果更加明显，这是因为模因组的数量与多核CPU的个数一致，同时模因组内的局部搜索过程是相互独立的，尽可能地减少了worker间的通信开销.虽然PSO算法和DE算法通过并行计算也能够减少一部分运行时间，但是由于个体间的进化需要更多的交互，worker间的通信开销较大，因此其减少的运行时间并不显著.这同时也说明了CPU并行计算在处理大规模种群进化时并不具有明显的计算优势，CPU+GPU的并行框架是目前处理此类问题的主要发展方向.但是总的来说，对于本文中蛙跳算法的模因搜索过程，多核CPU并行计算还是能够取得较好的加速效果.图 5进一步给出了本文算法在Dias模型下的反演结果和Ex及Hy分量的拟合情况.
|
|
表 8 不同非线性反演方法的反演性能 Table 8 Performances of different nonlinear inversion methods |
|
图 5 三层Q型地电模型中间层极化时反演结果 (a)电阻率；(b)极化率；(c) Ex拟合曲线；(d) Hy拟合曲线. Fig. 5 Inversion results of Q-type three-layer geoelectric model with polarized middle layer (a) Resistivity; (b) Polarizability; (c) Ex fitting curve; (d) Hy fitting curve. |
由于野外采集的数据均含有一定的噪声，本节以五层KHA型地电模型为例，研究本文算法在加入不同噪声情况下的反演稳定性和鲁棒性.其中地电模型参数(Dias模型)、反演算法参数和二阶段最小构造参数的设置如表 9所示.考虑到噪声对决定系数的影响，设定阈值th=0.96.
|
|
表 9 五层KHA型地电模型的仿真参数设置 Table 9 Parameter setting for KHA-type five-layer geoelectric model |
表 10给出了采用Dias模型，中间层极化时，加入不同程度噪声时本文算法的反演性能.图 6—8给出了本文算法在不同程度噪声下的反演结果和Ex及Hy分量的拟合情况.由图 6可知，在无噪声的情况下，本文算法较为准确地重构了KHA型地电模型的电阻率参数和极化率参数.同时由图 7和图 8可知，反演算法在10%和20%随机噪声的干扰下也能够较好地反演出模型的电阻率参数，同时也能够基本反演出极化率的变化形态.从数据拟合的结果来看，反演模型的Ex及Hy分量均位于加噪后的曲线之间，体现了反演算法的稳定性.由表 10可知，随着噪声的增加，反演的误差MSE和APE逐渐增大，决定系数R2逐渐下降，但总体来说本文算法在加噪条件下保证了较低的反演误差和较高的决定系数，这是因为最小构造约束函数和正则化因子增强了算法的鲁棒性.最后，五层模型反演的计算时间较三层模型反演的计算时间有所增加，这主要因为针对更加复杂的地电模型结构反演，增加了算法的种群规模和进化代数.
|
|
表 10 五层加噪地电模型激电信息提取反演结果 Table 10 Inversion results of IP information extraction for five-layer noisy geoelectric model |
|
图 6 五层地电模型不加噪声时反演结果 (a)电阻率；(b)极化率；(c) Ex拟合曲线；(d) Hy拟合曲线. Fig. 6 Inversion results of five-layer noiseless geoelectric model (a) Resistivity; (b) Polarizability; (c) Ex fitting curve; (d) Hy fitting curve. |
|
图 7 五层地电模型加10%噪声时反演结果 (a)电阻率；(b)极化率；(c) Ex拟合曲线；(d) Hy拟合曲线. Fig. 7 Inversion results of five-layer geoelectric model with 10% noise (a) Resistivity; (b) Polarizability; (c) Ex fitting curve; (d) Hy fitting curve. |
|
图 8 五层地电模型加20%噪声时反演结果 (a)电阻率；(b)极化率；(c) Ex拟合曲线；(d) Hy拟合曲线. Fig. 8 Inversion results of five-layer geoelectric model with 20% noise (a) Resistivity; (b) Polarizability; (c) Ex fitting curve; (d) Hy fitting curve. |
由于地下导电性不均匀所引起的异常强度远大于激电效应引起的异常，提取被强大导电和电磁效应异常“淹没”的弱小激电异常是一件非常困难的工作(罗延钟和张胜业，2003).本文对CSAMT信号中激电信息提取的非线性反演方法进行了深入研究.从反演的角度来看，从CSAMT信号中提取弱小的激电异常是一个复杂的多极值和不适定问题，需要采用非线性反演方法来提高算法的全局搜索能力，以便精确地拟合弱小的激电异常所引起的观测数据的微弱变化；同时需要突出激电异常对适应度函数的影响并提高算法的稳定性.本文从以上两个方面入手，提出了一种基于惯性权重-柯西分布的混合蛙跳算法，并结合二阶段最小构造反演的思想对含激电效应的CSAMT信号进行激电信息提取，在标准混合蛙跳算法的基础上，采用非均匀的柯西分布算子来取代均匀分布的随机算子并引入个体进化的经验来增强算法的全局搜索能力，通过采用混沌振荡的惯性权重来均衡进化过程中的个体经验和群体经验，提高算法后期的收敛速度和稳定性.在激电信息提取的反演过程中，通过引入二阶段的最小构造反演方法来进一步突出激电异常对适应度函数的影响，并将极化率参数的最小构造约束独立于对电阻率参数的最小构造约束，提高了激电信息提取的准确性.最后考虑到模因组间搜索的独立性，采用Matlab的并行计算工具箱实现了混合蛙跳算法模因组间搜索的CPU并行计算，提高了算法的计算效率.
对不同层状模型反演的结果表明：(1) 对于极化层位于不同位置时的层状模型，本文算法均能够较为准确地提取出极化率参数并重构真实的地电结构；(2) 采用不同的极化模型和不同的误差函数指标对反演的结果均存在一定的影响，在本文的实验条件下，采用Dias模型和电场分量Ex作为误差函数指标取得了更好的反演结果；(3) 相较于同类型的标准SFLA，DE和PSO等非线性反演算法，本文算法在适应度、反演结果和计算速度方面体现出了优越性，这是因为：首先基于柯西分布的随机算子令算法在前期获得了比均匀分布的随机算子更大的搜索范围，有利于保持前期种群的多样性，获得更优的全局搜索能力；其次混沌振荡惯性权重的引入控制了个体经验和群体经验之间的比重，保证了算法后期的收敛速度和稳定性；最后，基于CPU并行计算的模因组间搜索提高了算法的执行效率，其计算时间相较于其他非线性反演算法有明显的减少；(4) 为模拟实际工程中难以避免的噪声和干扰的影响，在正演数据中分别加入10%和20%的随机噪声，在噪声的干扰下，观测数据与预测数据的拟合度存在一定的下降.但电阻率和极化率模型参数的反演结果仍然能较为真实准确地显示实际地电构造和所含激电信息.这是因为最小构造约束在一定程度上降低了因为噪声所引起的虚假构造，而第二阶段的反演过程能够保证弱小的激电参数可以通过适应度函数对个体的搜索过程产生更好的引导，促使个体尽快地聚集在全局最优解附近，提高了激电信息提取的准确性和鲁棒性.
| Bäumer O. 1996. Inverting magnetotelluric data using genetic algorithms and simulated annealing. Inverse Methods: 223-230. | |
| Cao Z L, He Z X, Chang Y J. 2007. A simulation study of induced polarization effect of magnetotelluric and its application in oil and gas detection. Progress in Geophysics (in Chinese), 21(4): 1252-1257. | |
| Chen Q L, Hu W B, Li J M. 2007. Method for inversing real spectral parameters based on magnetotelluric method(MT). Journal of Oil and Gas Technology (in Chinese), 28(6): 61-64. | |
| Constable S C, Parker R L, Constable C G. 1987. Occam's inversion: A practical algorithm for generating smooth models from electromagnetic sounding data. Geophysics, 52(3): 289-300. DOI:10.1190/1.1442303 | |
| Dai Q W, Jiang F B. 2013. Nonlinear inversion for electrical resistivity tomography based on chaotic oscillation PSO-BP algorithm. The Chinese Journal of Nonferrous Metals (in Chinese), 23(10): 2897-2904. | |
| Dias C A. 2000. Developments in a model to describe low-frequency electrical polarization of rocks. Geophysics, 65(2): 437-451. DOI:10.1190/1.1444738 | |
| Eberhart R C, Kennedy J. 1995. Particle swarm optimization. Proceeding of IEEE International Conference on Neural Network. Perth, Australia, 1942-1948. | |
| Eusuff M, Lansey K, Pasha F. 2006. Shuffled frog-leaping algorithm: a memetic meta-heuristic for discrete optimization. Engineering Optimization, 38(2): 129-154. DOI:10.1080/03052150500384759 | |
| Feng B, Wang J L, Wang Y. 2013. The preliminary exploration to extract the IP effect which using CSAMT electromagnetic field response method. Progress in Geophysics (in Chinese), 28(4): 2116-2122. | |
| Gasperikova E, Morrison H F. 2001. Mapping of induced polarization using natural fields. Geophysics, 66(1): 137-147. DOI:10.1190/1.1444888 | |
| Ghorbani A, Camerlynck C, Florsch N. 2009. CR1Dinv: a matlab program to invert 1D spectral induced polarization data for the Cole-Cole model including electromagnetic effects. Computers & Geosciences, 35(2): 255-266. | |
| He J S. 2010. Widefield Electromagnetic Method and Pseudo Random Signal Electrical Prospecting (in Chinese). Beijing: Higher Education Press. | |
| He Z X, Hu Z Z, Luo W F. 2010. Mapping reservoirs based on resistivity and induced polarization derived from continuous 3D magnetotelluric profiling: Case study from Qaidam basin, China. Geophysics, 75(1): B25-B33. DOI:10.1190/1.3279125 | |
| Hu Z Z, He Z X, Yang W C. 2015. Constrained inversion of magnetotelluric data with the artificial fish swarm optimization method. Chinese Journal of Geophysics (in Chinese), 58(7): 2578-2587. | |
| Jiang F B, Dai Q W, Dong L. 2016. An ICPSO-RBFNN nonlinear inversion for electrical resistivity imaging. Journal of Central South University, 23(8): 2129-2138. DOI:10.1007/s11771-016-3269-8 | |
| Jiang F B, Dai Q W, Dong L. 2016. Nonlinear inversion of electrical resistivity imaging using pruning Bayesian neural networks. Applied Geophysics, 13(2): 267-278. DOI:10.1007/s11770-016-0561-1 | |
| Li D Q, Wang G J, Di Q Y. 2008. The application of genetic algorithm to CSAMT inversion for minimum structure. Chinese Journal of Geophysics (in Chinese), 51(4): 1234-1245. | |
| Liu M, Liu X, Wu M. 2011. Integrating spectral indices with environmental parameters for estimating heavy metal concentrations in rice using a dynamic fuzzy neural-network model. Computers & Geosciences, 37(10): 1642-1652. | |
| Luo W F, He Z X, Wang C F. 2011. A test study of induced polarization effects of magnetotelluric in oil and gas detection. Oil Geophysical Prospecting (in Chinese), 46(6): 978-983. | |
| Moorkamp M, Jones A, Eaton D. 2007. Joint inversion of teleseismic receiver functions and magnetotelluric data using a genetic algorithm: Are seismic velocities and electrical conductivities compatible?. Geophysical Research Letters, 34(16). DOI:10.1029/2007GL030519 | |
| Murali S R, Bhimasankaram V. 1980. Comparison of anomalous effects determined using telluric fields and time domain IP technique (test results). Exploration Geophysics, 11(12): 45-46. | |
| Pelton W, Ward S, Hallof P. 1978. Mineral discrimination and removal of inductive coupling with multifrequency IP. Geophysics, 43(3): 588-609. DOI:10.1190/1.1440839 | |
| Shaw R, Srivastava S. 2007. Particle swarm optimization: A new tool to invert geophysical data. Geophysics, 72(2): F75-F83. DOI:10.1190/1.2432481 | |
| Shi X M, Xiao M, Fan J K. 2009. The damped PSO algorithm and its application for magnetotelluric sounding data inversion. Chinese Journal of Geophysics (in Chinese), 52(4): 1114-1120. | |
| Spichak V, Popova I. 2000. Artificial neural network inversion of magnetotelluric data in terms of three-dimensional earth macroparameters. Geophysical Journal International, 142(1): 15-26. DOI:10.1046/j.1365-246x.2000.00065.x | |
| Storn R, Price K. 1997. Differential evolution—a simple and efficient heuristic for global optimization over continuous spaces. Journal of Global Optimization, 11(4): 341-359. DOI:10.1023/A:1008202821328 | |
| Tang R, Yu P, Xiang Y. 2014. The sensitivity analysis of different induced polarization models used in magnetotelluric method. Acta Geodaetica et Geophysica, 49(2): 225-233. DOI:10.1007/s40328-014-0050-z | |
| Wu H R, Wang S M. 1978. The feasibility of nature source induced polarization. Geophysical and Geochemical Exploration(1): 62-64. | |
| Xu W D. 2011. Study on effection and application of IP from CSAMT (in Chinese). Jilin: Jilin University. | |
| Yu C T, Liu H F, Zhang X J. 2013. The analysis on IP signals in TEM response based on SVD. Applied Geophysics, 10(1): 79-87. DOI:10.1007/s11770-013-0366-4 | |
| Yue A P, Di Q Y, Shi K F. 2007. The discussion on the IP information extraction from CSAMT signal. Progress in Geophysics (in Chinese), 22(6): 1925-1930. | |
| Yue A P, Di Q Y, Wang M Y. 2009a. 1-D forward modeling of the CSAMT signal incorporating IP effect. Chinese Journal of Geophysics (in Chinese), 52(7): 1937-1946. | |
| Yue A P, Di Q Y, Wang M Y. 2009b. 1D forward simulation of MT induced polarization (IP) effect of reservoir. Oil Geophysical Prospecting (in Chinese), 44(3): 364-370. | |
| 曹中林, 何展翔, 昌彦君. 2007. MT激电效应的模拟研究及在油气检测中的应用. 地球物理学进展, 21(4): 1252–1257. | |
| 陈清礼, 胡文宝, 李金铭. 2007. 由MT资料反演真谱参数的基本原理. 石油天然气学报, 28(6): 61–64. | |
| 戴前伟, 江沸菠. 2013. 基于混沌振荡PSO-BP算法的电阻率层析成像非线性反演. 中国有色金属学报, 23(10): 2897–2904. | |
| 冯兵, 王珺璐, 王玉. 2013. 利用CSAMT电磁场响应提取激电效应的方法初探. 地球物理学进展, 28(4): 2116–2122. DOI:10.6038/pg20130456 | |
| 何继善. 2010. 广域电磁法和伪随机信号电法. 北京: 高等教育出版社. | |
| 胡祖志, 何展翔, 杨文采. 2015. 大地电磁的人工鱼群最优化约束反演. 地球物理学报, 58(7): 2578–2587. DOI:10.6038/cjg20150732 | |
| 李帝铨, 王光杰, 底青云. 2008. 基于遗传算法的CSAMT最小构造反演. 地球物理学报, 51(4): 1234–1245. | |
| 罗卫锋, 何展翔, 王财富. 2011. 大地电磁激电效应油气检测试验. 石油地球物理勘探, 46(6): 978–983. | |
| 罗延钟, 张胜业, 熊彬. 2003. 天然场源激电法的可行性. 地球物理学报, 46(1): 125–130. | |
| 师学明, 肖敏, 范建柯. 2009. 大地电磁阻尼粒子群优化反演法研究. 地球物理学报, 52(4): 1114–1120. | |
| 吴汉荣, 王式铭. 1978. 利用天然电磁场进行激发极化法测量的可能性. 物探与化探(1): 62–64. | |
| 徐汶东. 2011. CSAMT中的IP效应影响及应用研究. 长春: 吉林大学. http://cdmd.cnki.com.cn/Article/CDMD-10183-1011103107.htm | |
| 岳安平, 底青云, 石昆法. 2007. 从CSAMT信号中提取IP信息探讨. 地球物理学进展, 22(6): 1925–1930. | |
| 岳安平, 底青云, 王妙月. 2009a. 含激电效应的CSAMT一维正演研究. 地球物理学报, 52(7): 1937–1946. | |
| 岳安平, 底青云, 王妙月. 2009b. 油气藏M激电效应一维正演研究. 石油地球物理勘探, 44(3): 364–370. | |