差分进化算法在MT信号激电信息提取中的应用研究
&&&&2015, Vol. 30 Issue (4):
董莉, 李帝铨, 江沸菠. 2015. 差分进化算法在MT信号激电信息提取中的应用研究[J]. 地球物理学进展, 30(4): , doi: 10.6038/pg&&
DONG Li, LI Di-quan, JIANG Fei-bo. 2015. The DE algorithm and its application for IP effect extraction of magnetotelluric sounding data. Progress in Geophysics, 30(4): , doi: 10.6038/pg &&
差分进化算法在MT信号激电信息提取中的应用研究
董莉1,2, 李帝铨1 , 江沸菠1,3&&&&
1. 中南大学地球科学与信息物理学院, 长沙 410083;2. 湖南涉外经济学院信息科学与工程学院, 长沙 410205;3. 湖南师范大学物理与信息科学学院, 长沙 410081
基金项目: 国家自然科学基金资助项目()、中国博士后科学基金资助项目()和湖南科技大学页岩气资源利用省重点实验室开发基金资助项目(E21423)联合资助.
作者简介: 董莉,女,1982年生,湖南长沙人,博士研究生,主要从事电磁法研究工作.(E-mail:)
通讯作者: 李帝铨,男,1982年生,广西玉林人,博士,副教授,主要从事电磁法研究工作.(E-mail:lidiquan@)
摘要: 从MT信号中提取激电信息的研究对深部矿产资源和油气资源的开发具有重要的现实意义.在激电信息提取方法的研究中,过去的理论多以线性反演方法为主,主要存在三个问题:(1)依赖初始模型,(2)容易陷入局部极值,(3)多解性严重.本文通过以下方法解决上述问题:(1)采用差分进化算法对MT信号中的激电信息以反演的方式进行提取,该方法受初始模型影响小,不易陷入局部极值.(2)在差分进化算法的适应度目标函数中引入最小构造约束,并提出对极化率和电阻率分别约束来提高反演结果的稳定性.对极化层处于不同位置的不同地电类型的反演仿真结果表明,本文算法能有效获得电阻率和激电参数的光滑模型,反演结果具有较高的稳定一致性和准确性.加入噪声测试后的实验结果表明,本文算法对高斯白噪声有较强的鲁棒性.与微粒群优化算法(PSO)比较结果表明,本文算法具有更为优越的收敛速度,能够获得更好的反演效果.上述结果表明,用非线性方法提取MT信号激电信息具有积极的意义.
差分进化算法&&&&
大地电磁测深&&&&
非线性反演&&&&
激电效应提取&&&&
The DE algorithm and its application for IP effect extraction of magnetotelluric sounding data
DONG Li1,2, LI Di-quan1 , JIANG Fei-bo1,3&&&&
1. School of Geosciences and Info-Physics, Central South University, Changsha 410083, C2. School of Information Science and Engineering, Hunan International Economics University, Changsha 410205, C3. College of Physics and Information Science, Hunan Normal University, Changsha 410081, China
Abstract: The study of the IP information extraction from magnetotelluric sounding data is of great and practical significance to the exploitation of deep mineral resources and oil and gas resources.The linear inversion method is given priority to on previous research of IP information extraction method, which has three main problems: (1) Depends on the initial model, (2) Easily falls into the local minimum, (3) Serious non-uniqueness of solutions. The following methods are used to solve the above problems in this paper: (1) Differential Evolution algorithm is used to extract the IP information in the form of inversion. The inversion results are less affected by the initial model and are not easy to get into the local extremum. (2)Minimum structure is employed to constrain the fitness function of DE algorithm so as to solve the problem of multi-solution in inversion. And we propose that restriction should be imposed respectively on polarizability and resistivity in order to improve the stability of the inversion results. The inversion simulation results of polarization layer in different stratum of various geoelectric model show that the smooth models of resistivity and IP parameters can be obtained by the proposed algorithm, the results of which have higher stability and accuracy. The experiment results added with noise indicate that this method is robust to Gaussian white noise. And when compared with PSO algorithm, this algorithm has better convergence rate and inversion effects. The above results show that nonlinear method plays a positive role in the IP information extraction of MT signal.
Key words:
magnetotelluric sounding&&&&
nonlinear inversion&&&&
IP effect extraction&&&&
激发极化(Induced Polarization，IP)是发生在地质介质中因外电流激发而引起介质内部出现电荷分离，产生一个附加的“过电位”(over voltage)的一种物理化学现象().基于岩矿石激发极化属性的激发极化法作为一种有效的地球物理勘探方法，已广泛地应用于矿产勘探、地下水勘查、地震监测、油气勘探等领域.
早在20世纪70年代，吴汉荣等就开展了对天然场源作激电法测量的可行性研究().2003年，罗延钟通过理论分析和数值模拟详细论证了当时已有的天然场源激电法方案的利弊，认为具有普适性的天然场源及人工源长偏移距激电法方案尚待探索().随着近年来对深部矿产资源和油气资源需求的不断增加，越来越多的研究人员致力于从探测深度较大的大地电磁测深法(MT)资料中提取激发极化参数的研究工作().2006，陈清礼等基于Dias模型，对三层地电模型中间层有无极化情况的正演进行了模拟对比，给出了从MT资料中反演地层真谱参数的算法，但未给出模拟仿真结果().同年，曹中林等使用Cole-Cole模型模拟激发极化效应，采用广义逆方法对MT理论模型和实际资料进行反演，获得了较好的效果().2007年，岳安平等尝试从可控源音频大地电磁法(CSAMT)信号中提取IP信息()，给出了提取的基本思路和步骤.2009年，岳安平等在MT正演研究的基础上，基于Dias模型，对含激电参数的CSAMT信号展开了正演研究，为提取CSAMT信号中所含激电信息提供理论基础().
2013年，冯兵等分析了Cole-Cole模型参数对CSAMT信号的Ex、Hy的影响特征，采用最小二乘反演方法，分别利用Ex和Hy分量反演电阻率、厚度和激电参数().总的来说，目前对于电磁法信号中激电信息的提取研究多以激电信息提取方法的可行性探讨及验证为主，重点讨论岩矿石激电响应模型参数对不同地电类型的激电响应的影响，通过对比有无激电效应的正演结果，总结由极化模型参数变化所引起的极化异常和真谱参数之间的对应关系的相关规律.
对极化率提取的具体实现基本上以传统线性反演方法为主.将非线性的地球物理问题线性化或局部线性化的局部搜索算法，如最速下降法、共轭梯度法、最小二乘算法、广义逆算法等较为简单，能较好地满足实际工程的需要，但往往依赖于接近真实模型的初始值才能成功收敛，容易陷入局部极值().差分进化算法(Differential Evolution，简称DE)是一种高效的全局优化算法.该方法采用实数编码，具有控制参数简单，全局搜索能力和鲁棒性强、收敛速度快并且结果精准等特点，非常适合求解许多不连续、不可微，或者没有函数解析式的全局优化问题(; ).潘克家等采用改进的混合差分算法在不同噪声水平下对四种正则化方法的反演结果进行了详细比较，验证了HDE算法的有效性和稳定性，为求解地球物理反问题提供了一种新方法().宋维琪等研究了微地震贝叶斯差分进化反演方法，对比分析了网格搜索方法反演结果，通过理论和实际验证了该研究方法的正确合理性.()
本文基于Dias复电阻率模型，在含激电效应的MT信号正演研究的基础上，采用差分进化算法对激电信息以反演的方式进行提取.为解决由于反演的多解性所带来的反演模型过于复杂的问题，消除多余构造，提高反演结果的稳定性，同时增强反演算法的抗噪声性能，本文将最小构造约束引入MT 1-D差分进化反演算法中.鉴于模型参数中极化率m的取值范围与ρ电阻率的取值范围有较大差距，本文提出对两者分别采用不同最小构造约束共同构成DE算法的评估函数.实验结果显示，模型参数反演结果稳定一致，边界分辨清楚，数值较为准确，较好地重现了真实的地电结构和激电信息，上述结果表明本文算法在解决MT信号激电信息提取问题时具有较佳的有效性和稳定性，并对高斯白噪声具有较强的鲁棒性.说明用非线性方法提取MT信号激电信息具有积极的意义.
1 方法原理 1.1 激电响应Dias模型
在含有可极化岩石的地球物理反演问题中，前人已提出Dias模型、Cole-Cole模型、复Cole-Cole模型等十几种激电响应模型，其中Dias模型从提出到改进求精，经历了三十多年的历史，无论是在对复杂实验数据的拟合精度还是参数物理含义的明确性方面，Dias模型都极具优势.本文亦采用Dias模型来模拟地电介质的激电效应.Dias模型等效电路如图 1所示.该模型中复电阻率ρ的表达式为()
其中，μ=iωτ+(iωτ″)1/2，τ=rCdl，τ'=(R+Rs)Cdl，τ″=(aCdl)2，(ρ0为不考虑激电效应的直流电阻率)，.时间常数τ为正实数，与产生极化效应的岩土颗粒的平均大小有关；电化学参数η与产生极化效应的电化学环境的特征相关，表征了感应电流部分和扩散系数部分的相对重要性，具有S-1/2的量纲.极化电阻率系数δ是定义在[0，1)区间的无量纲参数，只与孔隙路径有关，孔隙路径代表了极化作用产生的位置.极化率m也是定义在［0，1)区间的无量纲参数，其取值通常与极化效应的强度有关，m越大，极化效应越强.
1.2 正演理论
在此以一维水平层状地层模型为研究对象，设有n层地电断面，各层厚度分别为h1，h2，…，hn(hn→∞)，各层电阻率分别为ρ1，ρ2，…，ρn，各层极化率分别为m1，m2，…，mn.在大地电磁测深方法中，地面波阻抗是根据各层厚度及电阻率从下往上递推得到的，具体计算公式为()
其中.Zi为第i层地层顶面的复波阻抗；Z0i为第i层的本征阻抗；ki为第i层的复波数；hi为第i层的厚度；ρs为复视电阻率；Z为地面波阻抗；ρi为第i层的电阻率，与第i层的电导率σi互为倒数；ω为角频率；μ0为大地导磁率，通常取μ0=4π×10-7 H/m.当考虑激电效应时，相应层的实电阻率用(1)式所示的复电阻率取代，从而获得激电效应和电磁效应的综合响应.
图 1Fig. 1图 1 Dias模型等效电路Fig. 1 The equivalent circuit of Dias mode
1.3 差分进化算法
差分进化算法由美国学者Storn和Price于1995年提出()，因其形式简洁，控制参数少，收敛速度快等特点成为最优化问题研究方向的一个热点，在许多领域都得到了应用().作为一种进化计算类算法，差分进化算法同样基于自然界生物体物竞天择、适者生存的本质规律.与早期的遗传算法()类似，差分进化算法的基本操作也包括初始化种群、变异、交叉、选择.
DE/r and /1/bin，，简称DE1，是被广泛采用的差分策略()，具有实现简单，有利于保持种群多样性等特点.本文亦采用该种策略来实现MT信号激电信息的提取.
(1)初始种群
设种群规模为NP，可行解空间的维数为D，用X(t)来表示进化到第t代时的种群.在问题的可行解空间内随机产生初始种群为
其中xi0=［xi，10，xi，20，…xi，D0]，表示第i个个体解.个体的各个分量可按下式产生：
式中，xj，max和xj，min分别为解空间第维的上下界.
对于第G代种群中任意一个目标向量xi，G，DE1根据下列公式生成变异向量viG+1为
式中，{xr1G，xr2G，xr3G}是在第G代种群中随机选择的三个不同个体，且r1≠r2≠r3≠i.F被称为放缩因子，是介于[0，2]间的实常量，用于控制差分向量(xr2G-xr3G)的影响.
交叉操作的目的旨在通过变异向量viG+1和目标向量xiG各维分量的随机重组以提高种群个体的多样性.新的交叉向量uiG+1通过以下公式生成：
式中，r and b是[0，1]间的随机数；CR被称为交叉因子，是介于[0，1]间的常数，取值越大发生交叉的概率越高，同时影响着收敛速度.r and j是在[1，D]间随机选择的整数，它确保交叉向量uiG+1至少从变异向量viG+1中获得一个元素.
差分进化算法的选择操作是一种“贪婪”选择模式，当且仅当新的向量个体uiG+1的适应度函数值比目标向量个体xiG的适应度值更好时，uiG+1才会被种群接受.该过程可由下式描述：
式中，f(uiG+1)和f(xiG)分别表示uiG+1和xiG的适应度函数.选择的结果是使得子代个体总是优于或等于父代个体，从而使种群始终向着最优解的方向进化().
实验中，为平衡反演效果和计算效率之间的矛盾，我们选取的参数分别为：种群规模NP=30，迭代次数1000次，交叉概率系数CR=0.9，放缩因子F=0.6.
2 反演方案
(1)不考虑最小构造约束时的目标函数
在含激电效应的层状地电模型的应用背景下，MT信号DE算法反演的模型参数通常包括各层的厚度hi、电阻率值ρi以及极化率mi(h1ρ1m1，h2ρ2m2，…hiρimi，…hnρnmn，n为地电模型层数，i∈[1，n]，hn→∞.)故此激电效应反演问题的可行解空间的维数为D=3n-1.
差分进化算法中的“贪婪”选择模式保证了在目标函数的约束下，种群的更新始终向着保留精英个体的方向进行.在上述激电效应反演问题的D维可行解空间中，根据卡尼亚视电阻率的正演计算公式，我们采用观测数据和正演视电阻率之间误差向量 e 的L2范数来构建目标函数E(m，f)为()：
其中，ρiobs表示观测数据，ρipre表示预测数据，F(·)为MT正演算子，miinv为反演得到的模型参数，i为频点序号.
(2)最小构造约束时的目标函数
寻求能够拟合观测数据的唯一模型一直是地球物理反演所关心的问题，然而众所周知，反演的结果存在着多解性，如果仅以数据拟合度作为评价标准，那么往往会得到过于复杂的模型，引入多余的构造，造成反演结果的不稳定.Constable等为提高反演过程的稳定性，采用正则化思想，提 出了最小构造的概念，通过增加对模型参数变化过程的约束函数来消除伪构造，进一步得到简单光滑的地电模型().该方法简单有效，李帝铨()、汤井田()等在CSAMT反演研究中均采用了这一思想.本文采用的最小构造目标函数为
其中，miinv为反演得到的模型参数，包括各层的电阻率值ρi以及极化率mi.ρ1m1，ρ2m2，…ρimi，…ρnmn，n为地电模型层数，i∈[1，n]，在此除半空间外，各层层厚取固定值H，即h1=h2=…=hn-1=H，hn→∞.故该激电效应反演问题的可行解空间的维数为D=2n.
(3)最小构造差分进化反演算法的适应度函数
在差分进化算法中，获得具有较高精度和稳定性的模型参数的关键在于施加适当的约束条件，具体对于考虑激电效应的MT反演问题而言，便是找到观测和预测数据的数据拟合度与地电构造粗糙程度之间的最佳平衡点.因此，综合(7)式和(9)式，构成求解条件极值的辅助函数，即最小构造差分进化算法反演的适应度函数为
其中λ为某一常数，称为拉格朗日乘数.改变λ的大小，则调整了反演模型粗糙度在差分进化算法中的影响权重.λ值越大，则反演模型光滑度好，但观测数据和预测数据的数据拟合度较差；λ值越小，反演模型粗糙度大，而观测数据和预测数据的数据拟合度较好.
最小构造约束的引入，有效增强了差分进化算法的稳定性.但从实验仿真中我们发现，调整λ值对电阻率模型参数反演结果的改进效果要大大优于对极化率反演结果的改进效果，反演极化率的准确性和稳定性均差于电阻率的反演结果.作者个人认为，产生上述结果的主要原因在于极化率的取值空间(m∈[0，1])较另一模型参数电阻率的取值空间有较大差异(一般可认为ρ>>m)，由于相邻层极化率的差值相对于相邻层电阻率的差值要小得多，因此(9)式给出的最小构造约束方式对极化率的控制未能起到明显作用，而对电阻率反演效果的控制作用明显.
(4)改进后的最小构造差分进化反演算法适应度函数
因此，我们尝试对(10)式所示的适应度函数进行改进，将对极化率的最小构造约束独立于对电阻率的最小构造约束.从而得到改进后的适应度函数为
其中，R1(m，f)、R2(m，f)分别为对电阻率和极化率的最小构造约束函数，在此均采用(9)式.λ1、λ2分别为R1(m，f)、R2(m，f)对应的拉格朗日系数.我们对三层K型地层(中间层极化率m=0.3，λ1=1×10-4，λ2=5×10-4)随机重复五次实验，实验结果表明，采用新的适应度函数后，极化率的反演结果的稳定性有了明显改善，电阻率和极化率的反演结果如图 1所示.极化率的五次随机实验结果均不存在虚假激电信息，极化层的极化率无振荡，极化层的相对深度、范围，盖层与中间极化层的边界分辨效果改善作用明显.观测数据和预测数据的数据拟合度均较好，如图 2所示，计算“随机1”至“随机5”各组电阻率振幅实验曲线与正演曲线的决定系数分别为0.8，0.6，0.9998；电阻率相位实验曲线与正演曲线的决定系数分别为0.4，0.0，0.9978.差分进化算法适应度曲线如图 3所示，五次随机实验过程中，本文算法均显示出较优的适应度下降速 度，在100代之内适应度曲线基本上均呈现出线性下降状态，直至1000代终结，适应度值始终在减小，体现了DE的“贪婪”选择属性，显示了该算法高速高效的收敛性能.
图 1Fig. 1图 1 电阻率和极化率反演结果曲线(λ1=1×10-4，λ2=5×10-4)Fig. 1 Resistivity & polarizability inversed result curves(λ1=1×10-4，λ2=5×10-4)图 2Fig. 2图 2 正演曲线与反演响应曲线对比图Fig. 2 Forward and inversed curves图 3Fig. 3图 3 差分进化算法适应度曲线Fig. 3 Fitness curves of Differential Evolution algorithm
3 含激电效应的MT 1D 反演数值模拟
3.1 考虑激电效应时三层介质模型
我们研究了A、H、K、Q四种不同地电类型中极化层分别位于首层、中间层和底层时采用最小构造约束的差分进化反演方法的有效性.给出了模型参数(包括电阻率和极化率)反演结果的稳定性和准确性评估结果.从以上实验结果的比较可知：
(1)对于三层A、H、K、Q地电模型中极化层分别处于首层、中间层和底层的情况，本算法实现了对模型参数电阻率和极化率的反演提取，获得了良好的效果；
(2)极化率的提取效果如图 4所示，各地电模型中当中间层为极化层时的提取效果相对最好，边界分辨清楚，数值反演较为准确，不同地电类型对参数提取效果的影响较小；首层极化时，K型地电模型的提取效果相对最好，极化率搜索空间最大.A型稳定性较好，收敛速度较慢.H型收敛速度较快，稳定性相对较差.实验过程中，Q型控制参数调整难度较大，极化率搜索空间受限较大；底层极化时，A、H、K型提取效果较为接近，在合适的控制参数作用下，均获得了较为稳定一致的极化率模型参数.Q型的提取效果相对次于以上三种地电类型，曲线爬升相对缓慢，但仍然可以较为准确地反映底层极化的激电特性.
图 4Fig. 4图 4 极化率模型与反演结果对比图(λ1=1×10-4)(A)A型;(B)H型;(C)K型;(D)Q型.
(a)首层极化m1，2，3=[0.3，0，0]；(b)中间层极化m1，2，3=[0，0.3，0]；(c)底层极化m1，2，3=[0，0，0.3]Fig. 4 Polarizability model and inversed result(λ1=1×10-4)(A)Type A;(B)Type H;(C)Type K;(D)Type Q.
(a)Polarized first layer m1，2，3=[0.3，0，0];(b)Polarized middle layer m1，2，3=[0，0.3，0];(c)Polarized bottom layer m1，2，3=[0，0，0.3].图 5Fig. 5图 5 电阻率模型与反演结果对比图(λ1=1×10-4)(A)A型;(B)H型;(C)K型;(D)Q型.
(a)首层极化m1，2，3=[0.3，0，0]；(b)中间层极化m1，2，3=[0，0.3，0，0]；(c)底层极化m1，2，3=[0，0，0.3].Fig. 5 Resistivity model and inversed result(λ1=1×10-4)(A)Type A;(B)Type H;(C)Type K;(D)Type Q.
(a)Polarized first layer m1，2，3=[0.3，0，0];(b)Polarized middle layer m1，2，3=[0，0.3，0];(c)Polarized bottom layer m1，2，3=[0，0，0.3]
(3)在电阻率反演的准确性方面，如图 5所示，无论极化层位置如何，H、K型反演结果相对而言数值最为接近真谱参数，边界分辨最为清楚，稳定性最好，在此我们用归一化方差来衡量其稳定性.方差越小，反演的模型参数与其均值的偏离程度越小，则解越稳定，计算结果如表 1所示.
表 1(Table 1)表 1 参数拟合度评估表Table 1 Parameters fitting evaluation
地电类型ρ1(Ωm)ρ2(Ωm)ρ3(Ωm)m1m2m3模型参数稳定性评估(方差)电阻率拟合决定系数相位拟合决定系数最佳适应度值
ρ(Ωm)m平均最佳平均最佳
A型1005009000.3000.05700.05580.13330.03410.29360.13870.17460.17460.23440.05491.00001.00000.99480.997713.7561
A型10050090000.300.09860.32610.04720.09040.00890.16261.00000.25830.27030.04300.99991.00000.99380.997640.6687
A型100500900000.30.97761.00000.7250 0.41520.27780.11470.12670.01900.04300.00710.99991.00000.99580.998631.8958
H型4002008000.3000.01210.02030.01020.00780.00340.10280.29420.07890.09080.04300.99991.00000.99570.998623.8970
H型40020080000.300.00900.01020.04040.04550.02040.28220.16260.24630.16260.17460.99950.99970.98990.994045.7520
H型400200800000.30.02270.02040.01090.00910.02780.00710.06690.03100.00710.13870.99970.99990.98650.994644.2375
K型2005003000.3000.01220.00730.00400.03630.00400.29420.11470.01900.44970.01900.99940.99950.99440.996312.8574
K型20050030000.300.00020.02970.00830.00850.00590.04300.41380.13870.06690.07890.99970.99990.99700.99837.6109
K型200500300000.30..00220.01650.00140.13870.06690.01900.12670.03100.99810.99860.99520.996014.6711
Q型1000500500.3000.03000.02630.00650.02170.01220.04030.05460.00000.00740.01130.99991.00000.99950.999927.9151
Q型10005005000.300.08270.13780.15770.06450.07050.76070.77270.66500.18650.47360.99991.00000.99900.999541.2563
Q型100050050000.30.03570.01980.05110.05330.02370.04300.13870.24630.23440.11470.99950.99990.99870.999466.6681
表 1 参数拟合度评估表Table 1 Parameters fitting evaluation
(4)本文算法中对于极化层位于不同位置的不同地电模型情况均显示出较优的适应度下降速度.在相同实验条件下(如图 4(B)(b))，我们比较了本文算 法和微粒群优化算法(PSO)的适应度下降曲线，如图 7所示，PSO算法当迭代次数超过240代后，收敛速度明显变慢，直至1000代终止时适应度值也下降不多.实验结果直观显示了差分进化算法具有对模型参数较强的寻优能力，比PSO算法收敛速度更快，适应度值更低，算法结果稳定无振荡，体现了其在解决MT 1-D激电信息提取问题时的可行性及优越性.
图 7Fig. 7图 7 DE算法和PSO算法适应度比较曲线Fig. 7 Comparison curves of the fitness of DE and PSO algorithm
(5)图 6显示了观测数据和预测数据的拟合情况，图中的电阻率曲线跟观测数据均具有较高的拟合度，而相位曲线实际也具有尚佳的拟合效果.在此我们使用决定系数(DeterminationCoefficient，R2)来衡量数据拟合程度，其定义如(12)式所示为()
式中yi为预测数据，Yi为观测数据，n为频点数目.R2表征了预测值与测量值之间的相关程度，其值越大，表示两组数据之间存在着越明显的线性相关性.如表 1所示，在无噪的情况下，观测数据和预测数据具有较高的拟合系数.
图 6Fig. 6图 6 正演曲线与反演响应曲线对比图(A)A型;(B)H型;(C)K型;(D)Q型.
(a)首层极化m1，2，3=[0.3，0，0]；(b)中间层极化m1，2，3=[0，0.3，0，0]；(c)底层极化m1，2，3=[0，0，0.3].Fig. 6 Forward and inversed curves(A)Type A;(B)Type B;(C)Type C;(D)Type Q.
(a)Polarized first layer m1，2，3=[0.3，0，0];(b)Polarized middle layer m1，2，3=[0，0.3，0];(c)Polarized bottom layer m1，2，3=[0，0，0.3] 3.2 加入噪声以后的模型反演
实际工程中所采集的勘探数据通常都含有噪声干扰，我们对三层A、H、K、Q地电模型中极化层分别处于首层、中间层和底层的情况进行了算法的抗噪声性能分析，对4.1中的正演数据分别加入10%和20%的高斯白噪声以模拟含噪的勘探数据再进行反演，结果如图 8~图 11及表 2所示.
图 8Fig. 8图 8 模型参数加入10%噪声后的极化率反演结果(A)A型；(B)H型；(C)K型；(D)Q型.
(a)首层极化m1，2，3=[0.3，0，0]；(b)中间层极化m1，2，3=[0，0.3，0]；(c)底层极化m1，2，3=[0，0，0.3].Fig. 8 The inversion result of polarizability added 10% noise(A)Type A;(B)Type H;(C)Type K;(D)Type Q.
(a)Polarized first layer m1，2，3=[0.3，0，0];(b)Polarized middle layer m1，2，3=[0，0.3，0];(c)Polarized bottom layer m1，2，3=[0，0，0.3].图 9Fig. 9图 9 模型参数加入10%噪声后的电阻率反演结果(A)A型;(B)H型;(C)K型;(D)Q型.
(a)首层极化m1，2，3=[0.3，0，0]；(b)中间层极化m1，2，3=[0，0.3，0]；(c)底层极化m1，2，3=[0，0，0.3].Fig. 9 The inversion result of resistivity added 10% noise(A)Type A;(B)Type H;(C)Type K;(D)Type Q.
(a)Polarized first layer m1，2，3=[0.3，0，0];(b)Polarized middle layer m1，2，3=[0，0.3，0];(c)Polarized bottom layer m1，2，3=[0，0，0.3].图 10Fig. 10图 10 模型参数加入20%噪声后的极化率反演结果(A)A型;(B)H型;(C)K型;(D)Q型.
(a)首层极化m1，2，3=[0.3，0，0]；(b)中间层极化m1，2，3=[0，0.3，0]；(c)底层极化m1，2，3=[0，0，0.3].Fig. 10 The inversion result of polarizability added 20% noise(A)Type A;(B)Type H;(C)Type K;(D)Type Q.
(a)Polarized first layer m1，2，3=[0.3，0，0];(b)Polarized middle layer m1，2，3=[0，0.3，0];(c)Polarized bottom layer m1，2，3=[0，0，0.3].图 11Fig. 11图 11 模型参数加入20%噪声后的电阻率反演结果(A)A型(B)H型(C)K型(D)Q型
(a)首层极化m1，2，3=[0.3，0，0]；(b)中间层极化m1，2，3=[0，0.3，0]；(c)底层极化m1，2，3=[0，0，0.3].Fig. 11 The inversion result of resistivity added 20% noise(A)Type A(B)Type H(C)Type K(D)Type Q
(a)Polarized first layer m1，2，3=[0.3，0，0];(b)Polarized middle layer m1，2，3=[0，0.3，0];(c)Polarized bottom layer m1，2，3=[0，0，0.3].
从以上实验结果的比较可知：
(1)对勘探数据添加10%的高斯白噪声干扰后，本算法能实现对电阻率和极化率模型参数的反演提取，获得了与真实模型较为接近的光滑模型.与不考虑噪声的情况相比较，电阻率模型参数无论是在准确性还是稳定性方面仍然能显示出较佳的反演效果.实验结果显示，H型和K型反演效果最佳，光滑模型曲线与真谱参数曲线拟合度好，边界分辨清楚，数值反演较准确.Q型和A型相对次于前面两种类型，主要体现在中间层分界线相对模糊，但基本不影响对地电类型判断的辨识.极化层位于首层的Q型和H型的极化率稳定性低于无噪声的情况，其他模型参数反演质量在边界分辨，曲线形态，数值准确性方面与无噪情况相比无明显下降；
(2)对勘探数据添加20%的高斯白噪声干扰后，电阻率的反演结果仍然保持了较高的稳定一致性和准确性.极化率反演稳定性较前两种情况有所下降.A型和K型中极化层位于首层时，极化率下降速度较加噪10%的情况有所减小.A型和H型极化层位于底层时，极化率数值误差增大，勘探数据与预测数据的数据拟合度明显下降.Q型当极化层位于首层时，受噪声影响，极化率搜索空间受限增大.对于极化层处于不同位置的A、H、K、Q四种地电模型，算法的适应度值较加噪10%情况均有较大增幅.
(3)添加噪声后的视电阻率拟合决定系数较无噪时有所下降.噪声系数越大，决定系数下降越多.总体来说，视电阻率振幅拟合决定系数受噪声影响的降幅低于视电阻率相位拟合决定系数.其主要原因如前面所述，考虑到稳定可靠的高精度相位参数在实际工程中通常难以获得，故在本文实验环境下，未将相位列入模型参数参与反演，也没有与之相应的最小构造约束，故其对噪声的鲁棒性相对较弱.
(4)以上实验结果验证了差分进化算法在解决MT 1-D激电信息提取问题时的抗噪声性能.表 2所给实验结果直观显示了本文方法在不同的噪声环境中对模型参数进行反演和激电信息进行提取的优越性能.
表 2(Table 2)表 2 加噪数据拟合度评估表Table 2 Noise data fitting evaluation
地电类型ρ1(Ωm)ρ2(Ωm)ρ3(Ωm)m1m2m3噪声系数视电阻率拟合决定系数视相位拟合决定系数适应度
A型1000.310%0.98770.8705343.6770
900020%0.94580.5192785.9816
A型100010%0.97990.6926428.8146
900020%0.97660.5032487.4979
A型100010%0.98490.7532388.3384
9000.320%0.96200.5574603.1002
H型4000.310%0.96530.7653389.6014
800020%0.89370.3146632.6293
H型400010%0.97970.8533300.6401
800020%0.85090.4044807.9118
H型400010%0.97680.6857305.2360
8000.320%0.84210.3585745.7184
K型2000.310%0.95870.5748126.3131
300020%0.74590.5590331.3120
K型200010%0.80110.2441186.1597
300020%0.53550.1669367.4753
K型200010%0.84860.7346179.6518
3000.320%0.42250.2533385.8220
Q型10000.310%0.98900.9105255.8369
50020%0.96960.8012406.0653
Q型1000010%0.99890.8900116.1279
50020%0.98080.6042477.8089
Q型1000010%0.98580.6683416.2237
500.320%0.98330.8246472.4030
表 2 加噪数据拟合度评估表Table 2 Noise data fitting evaluation
4.1&&&&本文采用差分进化算法用反演的方式来提取含激电效应的MT 1-D信号中的激电信息并重现地电结构.结合了最小构造思想的差分进化算法在解决此类反演问题时显示出了良好的寻优能力，获得较为理想的实验结果.在适应度目标函数中加入对模型参数光滑度的约束，并将对极化率的最小构造约束独立于对电阻率的最小构造约束，选取合适的拉格朗日系数后，同时获得了极化率及电阻率反演的光滑模型，提高了模型参数反演结果的稳定性和准确性.
4.2&&&&在不考虑噪声的情况下，对三层A、H、K、Q地电类型的实验仿真结果较好地反映了地层的激电特性，证明了该方法的有效性；为模拟实际工程中难以避免的噪声和干扰影响，我们在正演数据中分别加入10%和20%的高斯白噪声，算法的抗噪声性能显示，本方法对噪声干扰具有较强的鲁棒性，获得的反演参数值具有较高的精度和准确性.表明了该方法在解决此类问题时具有一定的优越性.
4.3&&&&由于实际的地电结构非常复杂，实际工程环境中存在的噪声组成也十分复杂.因此上述方法用于实测数据反演时，还存在以下问题：(1)算法中的拉格朗日系数对不同地电类型以及极化层所在不同位置的反演效果有较大影响，只有合理地设置该系数才能取得较好的反演结果.在本实验中是采用凑试法确定该参数的，下一步希望通过自适应方法来确定该参数以提高算法的适用性；(2)尽管对两个参数分别采用独立最小构造约束，但在实际调试的过程中，电阻率反演的稳定性和准确性要优于极化率.并且两个约束系数λ1、λ2之间是否存在对反演结果的相互影响以及产生怎样的影响，尚需进一步定量分析研究；(3)DE算法本身的主要参数F、CR、NP的设置与DE算法的计算效率、有效性、鲁棒性均紧密相关().寻求DE算法参数的合理设置以有效改善激电信息提取效果，加快寻优收敛速度的研究也极具现实意义.
致 谢　感谢审稿专家提出的修改意见和编辑部的大力支持！
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..
Chen Q L,Hu W B,Li J M.2007.Method for Inversing Real Spectral Parameters Based on Magnetotelluric Method(MT)..
Cohn H, Fielding M. 1999. Simulated annealing:searching for an optimal temperature schedule. .
Constable S C, Parker R L, Constable C G. 1987. Occam's inversion:A practical algorithm for generating smooth models from electromagnetic sounding data. .
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):..
Dias C A. 2000. Developments in a model to describe low-frequency electrical polarization of rocks. Geophysics.65(2):437-451.
Duan H B,Zhang X Y,Xu CH F.2011.Bio-inspired Computing.Beijing:Science Press.
Feng B,Wang J L,Wang Y,et al.2013.The preliminary exploration to extract the IP effect which using CSAMT electromagnetic field respose method.PROGRESS IN GEOPHYSICS. 28(4):..
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 SH.2006.Dual frequency induced polarization method.Beijing:Higher education press.
Kim H-K, Chong J-K, Park K-Y, et al. 2007. Differential evolution strategy for constrained global optimization and application to practical engineering problems. .
Li D Q,Wang G J,Di Q Y,Wang M Y,Wang R.2008.The application of Genetic Algorithm to CSAMT inversion for minimum structure..
Liu G D,Shi Y J.1985.Magnetotelluric Sounding Method Tutorials.Beijing:Earthquake press.
Luo Y Zh,Zhang Sh Ye,Xiong B.2003.Feasibility of natural source induced polarization..
Pan K J,Wang W J,Tan Y J,et al.2009.Geophysical linear inversion based on hybrid differential evolution algorithm..
Shaw R, Srivastava S. 2007. Particle swarm optimization:A new tool to invert geophysical data. .
Storn R, Price K. 1995..
Song W Q,Gao Y K,Zhu H W.2013.The differential evolution inversion method based on Bayesian theory for micro-seismic data. CHINESE JOURNAL OF GEOPHYSICS.56(4):..
Tang J T,Zhang L C,Xiao X.2011.One dimension CSAMT mimimum structure inversion based on the frequency..
Wang P Ch,Qian X,Wang Y,et al.2009.Overview of differential evolution algorithm..
Wu H R,Wang S M.1978.The feasibility of nature source induced polarization..
Yang Q W,Cai L,Xue Y C.2008.A Survey of Differential Evolution Algorithms..
Yang SH Y,Zhang H.2012.The Realization of Swarm Intelligence and Bio-inspired Computing Using Matlab Technology.Beijing:PUBLISHING HOUSE OF ELECTRONICS INDUSTRY.
Yao Y.2002.The Basic Geophysics Inversion Theory & Application Method.Wuhan:China University of Geosciences press.
Yue A P,Di Q Y,Shi K F.2008.The discussion on the IP information extraction from CSAMT signal..
Yue A P,Di Q Y,Wang M Y.2009b.1-D forward modeling of the CSAMT signal incorporating IP effect..
Yue A P,Di Q Y,Wang M Y.2009a.1D forward simulation of MT induced polarization(IP) effect of reservoir..
曹中林,何展翔,昌彦君. 2007. MT激电效应的模拟研究及在油气检测中的应用..
陈清礼,胡文宝,李金铭. 2007.由 MT资料反演真谱参数的基本原理..
戴前伟,江沸菠,董莉.2014.基于汉南-奎因信息准则的电阻率层析成像径向基神经网络反演.地球物理学报.57(4):..
段海滨,张祥银,徐春芳. 2011.仿生智能计算.北京:科学出版社.
冯兵,王珺璐,王玉,等. 2013.利用 CSAMT电磁场响应提取激电效应的方法初探.地球物理学进展.28(4):..
何继善. 2006.双频激电法.北京:高等教育出版社.
李帝铨,王光杰,底青云,等. 2008.基于遗传算法的 CSAMT最小构造反演..
刘国栋,石应骏. 1985.大地电磁测深法教程.北京:地震出版社.
罗延钟,张胜业. 2003.天然场源激电法的可行性..
潘克家,王文娟,谭永基,等. 2009.基于混合差分进化算法的地球物理线性反演..
宋维琪,高艳珂,朱海伟,等.2013.微地震资料贝叶斯理论差分进化反演方法.地球物理学报. 56(4):..
汤井田,张林成,肖晓. 2011.基于频点 CSAMT-维最小构造反演..
王培崇,钱旭,王月,等. 2009.差分进化计算研究综述..
吴汉荣,王式铭. 1978.利用天然电磁场进行激发极化法测量的可能性.物探与化探.(1):62-64.
杨启文,蔡亮,薛云灿. 2008.差分进化算法综述.模式识别与人工智能.21(4):506-513.
杨淑莹,张桦. 2012.群体智能与仿生计算 Matlab技术实现.北京:电子工业出版社.
姚姚. 2002.地球物理反演基本理论与应用方法:武汉:中国地质大学出版社.
岳安平,底青云,石昆法. 2008.从 CSAMT信号中提取 IP信息探讨..
岳安平,底青云,王妙月,等.2009b.含激电效应的 CSAMT一维正演研究..
岳安平,底青云,王妙月,等. 2009a.油气藏 MT激电效应一维正演研究..