基于改进粒子群算法的大地电磁反演

李丽丽,, 李长伟,2*, 程勃,2, 陈汉波, 吕玉增,2, 熊彬,2, 张媛, 黄杨

(1. 桂林理工大学地球科学学院, 桂林 541000; 2. 广西隐伏金属矿产勘查重点实验室, 桂林 541000)

大地电磁(magnetotelluric sound, MT)反演方法大体可分为线性和非线性反演[1]。线性反演迭代收敛快,主要有高斯牛顿法[2]、共轭梯度法[3]、马夸特法[4]等,但是这些方法存在反演结果依赖于初始模型,易陷入局部最优解的缺点。非线性反演方法既能克服线性反演的缺点,又能有效地避免算法陷入局部最优[5]。因此,非线性反演方法一直是学者们研究的热点。常见的非线性反演方法包括粒子群算法(particle swarm optimization, PSO)、蒙特卡洛、遗传算法,神经网络算法等。由于粒子群算法原理简明、不依赖于初始模型、具有全局寻优、收敛速度快等优点,因此广泛应用于电力、地球物理、机器学习等各种领域中[6-9]。

PSO算法是由Shi等[10]在1995年通过模拟鸟觅食行为而发展起来的全空间搜索的算法,由Shaw等[11]最先应用在MT数据的反演中,从可行性分析来看,该方法可以有效解决电磁数据的反演问题。PSO算法目前主要应用于磁法、地震、测井、电法等方面。

赵平起等[12]在效差分算子的基础上,将BFP算法中的趋化、复制和扩散3个步骤加入PSO算法中进行反演,弹性波在不同模型上的数值模拟结果表明,改进算法不仅具有更高的精度,而且能够有效压制数值频散。陈杰等[13]将GA算法和PSO算法同时用于铁路隧道实测数据进行反演,结果表明两种方法对理论模型模拟的观测数据均具有较高的反演精度,且PSO算法的迭代收敛速度更快,但在加入噪声后的反演中GA算法却具有更强的抗干扰能力。王书明等[14]采用Levy flight 飞行策略结合了短距离搜索与偶尔长距离游走的特点,用于求解瞬变电磁二维反演问题,反演结果表明改进算法在保证算法收敛效率的同时提高了算法收敛全局的能力。李曦等[15]在优化支持向量机的基础上采用PSO算法进行岩性识别,不仅正确识别了岩性,而且提高了测井储层解释精度。Abril等[16]将带有遗传算子的PSO算法(EMPSO)用来解决二维电阻率成像问题,并通过将纯MPI和混合MPI-OpenMPI并行来解决时间问题,真实数据反演效果良好,与没有并行的相比,提高了计算时间。PSO算法经过多年地发展,已经衍生出很多优化的算法。但是大多数学者都是针对算法参数学习因子, 和惯性权重的优化,并未从根本上改善算法前期容易陷入局部极值的情况。

针对PSO算法维度越高,越容易陷入局部极值,大地电磁传统反演方法依赖初始模型的缺点,现提出一种改进的粒子群算法(cooperation adaptivity particle swarm optimization, CA-PSO),来进行大地电磁反演模拟,并在几个测试函数和物理模型中进行了测试,与其他改进的粒子群算法进行了对比,讨论该算法的性能。CA-PSO算法一维反演的目标函数是利用观测视电阻率和正演得到的视电阻率构建,二维反演的目标函数在一维的基础上添加核函数以及先验信息的限制。为了改善PSO算法容易陷入局部极值,增加粒子的活性,与标准粒子群算法相比,改进算法的迭代方式多了局部进化,将静态惯性权重替换为动态惯性权重,并且增加了收缩因子参数。

1 改进的粒子群算法(CA-PSO)

1.1 粒子群算法

粒子群算法起源于对鸟群运动行为的研究。每只鸟在某处位置能够找到食物的可能通过适应度来刻画,每只鸟都能记着它的觅食地点,并找到最好的(局部最优,相当于极值点),鸟群中所有个体的最佳位置就可以看作整个鸟群的最佳觅食点(全局最优,相当于最值点)。整个鸟类群体的觅食活动必然会移动到整个搜索范围的最佳觅食区,通过不断改变鸟群觅食位置,即不断迭代更新位置,速度也不断更新,鸟群向最优位置靠拢。

原理如下:在一个D维解空间中,有N个粒子构成一个群体,每个粒子均以一定的速度移动。其中第i个粒子表示为一个D维向量Xi=[Xi1,Xi2,…,Xid]T(i=1,2,…,N,1≤d≤D),即第i个粒子在D维空间中的位置为Xi,将Xi代入目标函数就可以计算出粒子i的适应度,依靠适应度值的大小来衡量粒子Xi的好坏;第i个粒子的速度表示为Vi=[Vi1,Vi2,…,Vid]T;第i个粒子自身的最优值为Pbesti=[Pbesti1,Pbesti2,…,Pbestid]T;种群粒子最优值为Gbesti=[Gbesti1,Gbesti2,…,Gbestid]T。粒子速度和位置每次更新公式为

(1)

(2)

式中:w为惯性权重,是用来均衡全局以及局部搜索的能力;c1和c2为学习因子,代表粒子向自己和群体中有优秀粒子学习的能力;然后朝着历史最优点Pbest靠拢,一般取c1=c2=2;r1和r2为[0.1]的伪随机数。粒子通过多次迭代更新来贴近最优解,最终得到的Gbest就是算法找到的最优解。

由式(1)可知粒子下一步运动速度由三部分决定:惯性权重和上一步的速度项,自我学习因子和个体最佳位置项,群体学习因子和群体最佳位置项。由式(2)可知下一步粒子,移动的位置由上一步的位置和上一步的速度来决定,具体流程如图1所示。

在CA-PSO算法中,用局部进化替代流程图中红框圈出来的步骤

1.2 改进粒子群算法

本文提出的CA-PSO优化算法的基本原理是:在原标准PSO算法的基础上,增加局部进化策略来提高粒子迭代的效率,引入线性递变的惯性权重来确保粒子速度过大跳出搜索区域,过小陷入局部极值。同时也引入收缩因子[17]来均衡算法的收敛性,同时在算法迭代时,采用轮换法则进行计算。

算法局部进化策略:假设初始种群的个数为N,小组个数为Nm,小组内的粒子数为Nn=N/Nm,将种群粒子代入适应度函数获得适应度,这里为了编程方便,按适应度大小进行降序排序,按照转轮法则分为Nmm。例如:分为10组,每组粒子数就有N/10个,则Nn=N/10,所有粒子就被表示为一个Nn行、10列的矩阵,为了便于计算,第一行数编号为(1, 1),(1, 2),……,(1,Nn);第二行编号(2, 1),(2, 2),……,(2,Nn), 其他行依次类推,这里采用轮换法则是为了便于编程计算。

收缩因子计算公式为

(3)

式(3)中:Nn为小组内粒子个数;i为当前粒子数。

惯性权重的计算公式为

(4)

式(4)中:wmax为最大惯性权重,一般取0.9;wmin为最小惯性权重,一般取0.4;T为最大迭代次数;t为当前迭代次数。

每组中有一个最优粒子,组内迭代Nger次,适应度越好收缩因子S越小,适应度越差,收缩因子S越大;当某个粒子连续Nger次都是适应度比较差的,则将这个粒子初始化;组内迭代完后再次按适应度分组进行组内迭代,直到迭代次数足够或者最佳适应度满足所需条件时迭代停止。

小组内更新公式为

(5)

(6)

1.3 算法流程

基于以上理论得到本文算法(CA-PSO),算法的实现步骤分为全局进化步骤和局部进化步骤两部分。

全局进化步骤如下:①初始化种群(粒子群共有N个粒子)中的粒子数量,并给每个粒子一个随机的初始位置和初始速度;②根据个体粒子的适应度,将它们按降序排列,得到种群的全局最优粒子Pbest;③按照转轮法则将粒子划分为若干个小组;④在每个小组内部实行CA-PSO的局部进化策略;⑤判断组内局部进化步骤是否结束,是转入步骤⑥,否转入步骤④;⑥判断种群迭代有没有结束,是转入步骤⑦,若没有,则转入步骤②。⑦最佳适应度值所对应的全局最优粒子,就是最终结果,算法结束。

局部进化步骤如下:①初始化小组计数器ik和小组内粒子迭代次数Nger;②根据小组计数器的值选择局部进化的小组;③计算出小组内的最好个体xmi;④依据式(5)和式(6)更新粒子位置和速度;⑤判断粒子位置是否越过边界,若是,则按无形/吸收的方式重新生成粒子;⑥计算组内粒子的适应度,并排序,按适应度更新收缩因子;⑦确定组内的迭代次数是否大于最大迭代次数,是继续步骤⑧,否继续步骤③;⑧令迭代计数器ik=ik+1,判断ik计数器是否大于种群数,否转入步骤②,是判断组内每个粒子的历史最优位置是否发生变化,若没有,且适应度排序靠后,否则重找此粒子,跳出局部进化。

2 数值实验

2.1 测试函数

基于上述分析,利用MATLAB实现了CA-PSO算法,选取6个经典测试函数[18]来检验CA-PSO算法的综合性能,函数详细信息如表1所示,最佳值是测试函数在搜索空间中找到的最佳结果。

表1 测试函数

2.2 对比算法

为进一步验证CA-PSO算法的求解精度和收敛速度,将CA-PSO与其他PSO改进算法进行对比,分别为:粒子群算法PSO,包含变异的惯性和经验相互影响的算法(BEPSO)[19],重构了模型收缩因子的算法(CIPSO)[20]、自适应粒子群算法(IPSO)[21]。实验在相同的硬件环境下进行,种群大小N为500;维数D为500;最大的迭代次数为1 000;算法运行次数是20,上述对比算法的详细参数取值如表2所示,每种算法运行20次后取其平均最优值,结果如表3所示。对比表3中的数据可以看出,与其他4种算法相比,CA-PSO在函数f1～f6上的均值都为0,这说明算法的寻优精度比其他4种算法都高。同时BEPSO算法在f1偶次多项式函数和f4具有较高寻优难度上的均值不为0,说明该算法收敛度和收敛精度较低,陷入了局部误差。

表2 对比算法的参数列表

表3 算法结果对比

2.3 收敛曲线分析

图2和图3显示了4种算法在6个测试函数上的适应度平均值的变化曲线。可以很直观地看出,随着迭代次数的增多,CA-PSO算法在6种测试函数上都能快速收敛到最优解。f2、f3、f6中CIPSO是最慢找到最优解的,f1、f4中BEPSO是最慢收敛到最优解的。结果表明,CA-PSO算法比其他4种对比算法有更好的全局搜索能力和更快的收敛速度。

图2 单峰测试函数收敛图

图3 多峰测试函数收敛图

3 MT反演

3.1 MT一维反演

一维定义理论模型和实际模型的相对均方差作为目标函数φ,即

(7)

构建了2个三层模型:K和H型,模型参数如表4所示。取初始种群数50,空间维度为5(3个电阻率,2个层厚),最大迭代次数1 000,最大误差界限常数ε=1×10-4,反演结果如表5所示。

表4 模型参数

表5 反演结果

一维大地电磁改进粒子群算法反演具体步骤如下。

(1)初始化:根据模型层数和反演参数的多少及问题复杂程度确定相关参数,粒子维度D与地层参数(2n-1)个数一致。

(3)择优:根据评价的情况,选择当前最优适应度值对应的粒子,并判断是否满足中止条件,如满足,则停止搜素;否则,进入下一步种群更新。

(4)更新:对粒子群按照式(5)和式(6)更新当前种群(小组内)个体对应的参数值,再更新个体粒子适应度,从而完成当前整个种群的更新。

(5)进入下一代循环:算法转至步骤(2)继续执行,直到:①不满足迭代次数t>2且φ<ε(最大误差界限常数)时;②最大迭代次数为100;③φ<1×10-5中的任何一个条件时迭代中止。

CA-PSO算法反演与PSO算法反演流程区别在于步骤(4)中多了一步添加惯性权重和收缩因子的小组迭代操作。

由表5可以看出,H型曲线的反演电阻率均比理论模型电阻率小,最大60 Ω·m,最小0.1 Ω·m,层厚但也在40 m以内有略微差异;K型曲线反演电阻率均比理论模型大,范围在0.02～42 Ω·m,层厚更接近理论模型,差距在10 m以内,K型反演效果更好,总体来说该方法可用。

一般来说,H型测深曲线模型在反演中最难解决,所以这里只给出H型曲线模型的反演结果图,如图4所示。从图4可以看出,LS与模型的一致性最高,而粒子群算法的一致性最差。CA-PSO方法与理论模型在650 m左右的深度范围内曲线基本一致,表明CA-PSO方法适用于一维MT是可行的。

图4 H型模型结果图

3.2 MT二维反演

CA-PSO用在MT二维反演上,由于维度增加,计算量大且计算时间长。为了改善计算效果,二维反演的目标函数跟一维相比加入核函数[22-23]。例如:一组输入样本{xi,yi}i=1,2,…,m,输入量为xi∈Rn,输出量为yi∈R,通过核函数k(·,·)将在低维映射到更高的维度。目标函数为

φ=k(xi,yi)w+b

(8)

式(8)中:w为惯性权重;b为偏置量,增加平移的能力,一般取值为2;{xi,yi}i=1,2,…,m类比为

式(8)的求解转化为规划问题的最优解。

核函数取为Gaussian核函数,即

(9)

式(9)中:σ为输入参数,取值为1×10-3。

(10)

为验证算法在二维反演的适用性,设计了从简单到复杂的地电模型进行测试,地表布设一条长6 km的测线,13个测点,频率取0.01、1、100 Hz,计算这些测点处的TE(横电波)和TM(横磁波)极化模式下的视电阻率,来进行拟合反演寻找与所设异常体模型相对应的在背景场中反演异常体的位置。设置初始种群数N=100, 搜索空间维度D=10,最大迭代次Nger=1 000,学习因子c1=c3=2,c2=0.8,搜索空间的其他参数跟一维类似。背景电阻率均为100 Ω·m。

(1)模型一:该模型(图5)反映的是1个低阻异常体,埋深1 km,规模为2 km×1 km,异常体电阻率10 Ω·m。由图6可以清楚地看出,异常体位置与实际模型相符合,表明该方法反演单个规则异常体效果很好。

图5 单异常低阻模型

图6 单异常低阻反演结果

(2)模型二:该模型(图7)反映的是2个一样大小的低阻异常体,规模为2 km×2 km,埋深 1 km,异常体电阻率均为30 Ω·m。由图8可以明显地看出,异常体与实际模型的位置非常相似,表明该方法对于反演2个同样大小的低阻异常体是可行的。

图7 双异常低阻模型

图8 双异常低阻反演结果

(3)模型三:该模型(图9)反映的是一个电阻率为30 Ω·m,规模3 km×0.5 km,埋深 0.5 km的低阻异体,和一个电阻率为200 Ω·m,规模3 km×3 km,埋深 2 km的高阻异常体。 由图10可以看出,图10(b)比图10(a)能更好地反映实际模型,这表明在目标函数中只引入核函数和又加入模型修改量相比,后者反演效果更好,而且可以很清楚地区分高低阻。对于图10(a)反演结果较差,是因为少了模型修改量和一些限制。综上,改进算法CA-PSO应用在MT二维反演上可以很好地定位异常体位置,而且对于多个异常体也可以明显地反演出来,该改进方法用在二维反演是可行的。

图9 高低阻双异常模型

图10 高低阻双异常反演结果

4 结论

基于以上内容,得到以下结论。

(1)提出了一种基于增加局部进化策略和加入收缩因子的CA-PSO优化算法。算法加入了含有递变的惯性权重和收缩因子的局部进化步骤,将原来需要每个粒子之间相互对比的方式,变成了只需每个组的最优粒子进行对比,不仅兼顾了粒子集的多样性和收敛性而且提高了粒子进化的效率,也较好地避免了算法过早收敛和陷入局部误差的问题;在测试函数上的对比结果表明,该算法不依赖于初始模型,计算简单,在一定程度上改善了PSO算法容易陷入局部极值的优点。

(2)通过设置目标函数,在目标函数中加入先验信息限制,利用CA-PSO算法进行一、二维MT反演,能够准确地进行异常体的定位,反演结果与实际模型吻合,表明了算法应用于MT反演的可行性。

(3)利用CA-PSO算法实现二维MT反演,在目标函数中引入核函数,将传统拟合函数所使用的二范数改为一范数,能够减少计算量,得到准确的反演结果。

(4)由于粒子群算法是全局搜索,维度越高,计算量越大,计算时间也越长,所以目前很少有人用改进的PSO算法来做MT的二维反演,大都做的是一维反演。所以本文研究做二维反演只是一种尝试,而且每计算一次时间较长,计算结果也是理想化的,没有考虑实际的因素。对于减少二维反演计算时间和实际问题的研究,是下一步的研究计划。