瞬变电磁数据L-PSO反演方法

王书明， 底青云， 夏彤， 任子乾， 宋江涛， 邹贵安

1 中国地质大学(武汉)地球物理与空间信息学院， 地球内部多尺度成像湖北省重点实验室， 武汉 430074 2 中国科学院地质与地球物理研究所， 北京 100029 3 中国科学技术大学地球和空间科学学院， 合肥 230026

0 引言

反演问题是地球物理学中的核心问题，是从观测数据求解地质模型的过程.由于实际地层理论上属于无限维的空间，然而观测数据往往是离散且有限的，因此，反演问题的解必定是非唯一的(杨文采，1997).除了观测资料的局限性之外，场的等效性和观测数据的误差也会使得解的非唯一性更加严重.

基本的反演方法可分为线性反演和非线性反演两大类方法.由于大多数地球物理反问题都是非线性的，用传统的线性反演方法很容易陷入到局部最优解当中，并且线性反演方法非常依赖初始设置的模型.因此，20世纪中期之后，出现了地球物理非线性反演方法，PSO(Particle Swarm Optimization)算法就是其中运用较为广泛的一种反演方法.

PSO算法在地球物理领域的发展可追溯到20世纪，由印度的Shaw和Srivastava(2007)首先应用于大地电磁(MT)数据的反演中，并做了可行性分析，结果表明该方法可以解决电磁数据的反演问题.PSO方法被引进国内后，主要应用于地震和磁法勘探反演解释工作中，如在地震波阻抗反演(易远元等，2007)、磁资料井地联合反演(张大莲等，2009)、板状体磁异常数据反演(吴招才和刘天佑，2009)；在电磁法及瞬变电磁法方面，Olalekan和Di(2018)、易远元和王家映(2009)、李明星等(2014)、肖敏(2010)、师学明等(2009)、程久龙等(2014)、底青云等(Di et al., 2018)学者发表了一些相关文章.PSO优化算法公式简单实用，经过多年发展，衍生出很多改进型的算法.但是以上多数工作都是针对于参数c1，c2以及惯性权重ω的改进，并没有从根本上改善算法前期容易陷入局部最优解的境况.

Levy飞行轨迹是Levy提出，Benoist-Madelbrot进行描述的一种稳定分布(Mandelbrot et al., 1983).飞行是一种随机游走过程，它的步长是一种连续的重尾分布.Yang和Deb(2009)首次提出了一种称为布谷鸟算法(CS)的优化算法，在该算法中引入了Lévy flight策略，以此提高长距离搜索能力，避免搜索陷入局部解的情况，取得了良好的效果.侯慧超(2014)针对粒子群算法容易陷入局部最优解的问题，引入莱维飞行的策略，增强了种群的多样性，也提高了粒子群算法的收敛速度和求解能力.牛海帆等(2016)完成了莱维飞行与粒子群的混合搜索算法.李焰华等(2018)将Lévy flight与乌鸦算法(CSA)相结合，实验结果表明LCSA算法(Levy Crow Search Algorithm)有效地降低了“乌鸦群”在没有“领导者”情况下搜索的盲目性.王学武等(2017)利用莱维飞行粒子群算法求解焊接机器人路径优化问题.因此本文将此策略引入PSO算法中，通过与目前经典的PSO算法及其他优化算法做对比，探究Lévy flight 对PSO算法的影响.

1 方法原理

1.1 经典PSO优化算法及其改进型算法

粒子群优化算法，是受自然中群体行为的启发发展而来的优化算法，在一个种群中，每个个体都是一个“粒子”，种群(如鸟群)在寻找食物过程中，每个个体都在自己周围的一个范围内寻找，通过行为的信息交换，群体会整体向着食物源的方向移动，最后寻找到食物源.在PSO算法中，每个粒子都有个体的位置和速度，对于第i个粒子，定义自己的位置：xi=(xi1,xi2,…,xin)和速度：vi=(vi1,vi2,…,vin).通过三个量确定和更新自身的位置和速度，分别是自身惯性速度vi=(vi1,vi2,…,vin)、个体最优位置pbesti=(pbesti1,pbesti2,…,pbestin)和群体最优位置gbest=(gbest1,gbest2,…,gbestn).并且每个粒子通过(1) (2)式更新自身的速度和位置：

(1)

(2)

其中，ω是惯性权重；rand()是(0,1)区间内的随机数；c1,c2为学习因子，取值在区间[0,4]，一般取c1=c2=2；t为迭代次数.

由(1)式可以看出，惯性权重系数ω决定了个体的惯性速度vi，较大的惯性权重系数使得迭代时的速度变大，使得粒子“飞”的更远，有利于算法跳出局部极小，但是较大的ω会不利于算法的收敛，因此为了算法能够在早期迅速收敛到全局极值附近，应选取较大的ω，在迭代次数增加之后，为了算法能够迅速收敛，应选取较小的ω.对于参数选取策略，有很多学者做了改进工作.

林卫星和陈炎海(2011)提出的改进粒子群(Modified Particle Swarm Optimization, MPSO)算法，针对惯性权重系数早期较大、晚期较小的特点，根据最佳阻尼比的原理，对标准PSO算法的惯性权重系数做出修改：

(3)

(4)

其中，φ1=c1rand()，φ2=c2rand()，rand()为符合均匀分布的区间(0,1)上的随机数.

赵玉静(2011)提出的MPSO，是在公式(3)、(4)的基础上对速度更新方程做出了调整，将第二项当前迭代的位置用当前迭代的局部最优解代替，公式如下：

(5)

其中，惯性权重系数ω的取值策略与(4)式相同.该方法使得全局最优解与局部最优解影响更大，加快算法收敛速度.

师学明等(2009)的离散粒子群(Discrete Particle Swarm Optimization,DPSO)算法，对惯性因子采取非线性波动递减的取值策略.使得惯性权重系数ω的振荡递减呈现出与模拟退火法类似的退火过程，其速度更新方程等同于(3)式，惯性权重系数取值策略如下：

ω=0.99k·rand()/2+α,

(6)

其中α为常数，取值为[0.1,0.5]区间内，α常取0.1.

1.2 Lévy flight基本理论

Lévy flight是一种非高斯随机过程(李煜和马良，2012)，其步长服从Lévy分布.Lévy分布是一种“重尾”的稳定分布，用以描述自然界中多数种群的觅食游走行为特征.同时，Lévy flight也是以偶尔发生长步长跳跃为特点的一类马尔科夫性质的随机游走.其位置的更新方程为：

(7)

其中，α为步长因子，常取α=0.01，⊕算子表示点对点的乘积，Lévy(λ)表示服从参数为λ的连续的Lévy分布，其中1<λ<3：

Lévy(λ)～μ=t-λ,

(8)

地球物理反演模型参数通常是有限的，即xi=(xi1,xi2,…,xin)为有限且离散的，所以一般采用Mantegna(1994)提出的离散化公式模拟Lévy分布函数：

(9)

L(λ)为(8)式的Lévy(λ)，参数μ和v都服从正态分布，其中，v服从标准正态分布：N(0,1)；μ服从：N(0,σ2)；σ由(10)式确定：

(10)

其中，Γ(λ)是参数为λ的标准Gamma函数.λ一般取值为1.5.

由以上可得Lévy flight的位置更新方程为：

(11)

1.3 L-PSO算法基本理论

PSO算法在迭代初期由于粒子的多样性，使得在早期具有较强的全局搜索能力，但当迭代进行，粒子会向着当前最优解聚集，使得算法出现早熟现象，尤其是当粒子数不足时，存在着忽略全局最优的小概率事件，一旦迭代进行，容易陷入局部最优.而Lévy flight重尾分布的特点，使得粒子即使在局部最优的情况下，也能够有机会进行长距离跳跃，正好弥补了PSO算法容易早熟的不足.

粒子群算法根据不同的参数策略形成了不同变种的算法，其统称为PSO算法，Lévy flight与具体的算法并无关系，因此可以引入任何一种变体的PSO算法之中，本文将引入Lévy flight之后的任何PSO算法统称为L-PSO算法，涉及具体算法时，例如：MPSO引入Lévy flight之后称之为LMPSO；经典PSO算法引入Lévy flight之后简称为LPSO算法.

L-PSO算法的具体步骤，实际上是在粒子群进行完一次位置更新后，再进行一步Lévy flight随机游走，因为Lévy flight偶尔长距离的跳跃的特点，使得部分粒子能够跳出局部进行搜索.从而增强粒子群搜索全局极值的能力.具体实施如下：

(1) 初始化粒子群：设定粒子总数，并对每个粒子初始位置xi=(xi1,xi2,…,xin)和初始的速度vi=(vi1,vi2,…,vin)随机选择.设定最高迭代次数tmax以及截断误差RMS，每次迭代时速度v的取值范围，c1、c2的取值以及惯性权重系数ω的取值策略；

(2)计算粒子群的适应度值：根据适应度函数(即目标函数)，计算每个粒子的适应度值，并根据适应度值更新每个粒子的个体最优位置pbesti=(pbesti1,pbesti2,…,pbestin)和群体最优位置gbest=(gbest1,gbest2,…,gbestn)，在第一次迭代时，个体最优位置可设为初始位置，群体最优位置是适应度最好的粒子位置；

(4)判断：是否达到误差要求范围内或达到最大迭代次数，如否，进行步骤(3)，依此进行，直到满足终止条件，停止迭代，输出结果并保存.

图1展示了L-PSO优化算法的基本流程.

图1 L-PSO优化算法流程图

2 数值模拟

为了验证L-PSO算法的效果，本文在Matlab2014a平台上编程实现了不同的惯性权重的取值策略的算法，并设计了数值模拟实验，取得了理想的结果.

2.1 常用多极值测试函数

全局优化算法主要有两个指标决定其优化能力，一个是在多极值函数中搜索到全局最优即全局搜索能力；第二是收敛能力，即收敛于全局最优解的精度和速度；因此我们选取了两个典型的函数：Aekley函数与Schaffer函数.

Aekley函数的表达式如下：

xi∈[-5.12,+5.12].

(12)

Aekley函数在定义域内具有极密集的极值点，经常被用来测试优化算法的全局搜索能力.该函数只存在一个全局极值点，即(x,y)=(0,0)时，f(0,0)=0(如图2).

图2 Aekley函数图像

为对比引入Lévy flight后的PSO算法收敛效果，记加入Lévy flight后的经典PSO算法、MPSO(林卫星和陈炎海，2011)以及DPSO算法(师学明等，2009)分别为：LPSO、LMPSO和LDPSO.因为PSO算法也是一种随机算法，因此我们对每种(引入Lévy flight前和引入Lévy flight后，共6种)算法都进行多次试验，不设精度，以迭代50次结果平均值作为参考值.

经过10次重复试验，6种算法都搜索到了全局极值点，说明PSO算法的优良特性.由图3a对比来看，凡是引入Lévy flight游走的之后的算法，如：PSO与LPSO、MPSO与LMPSO以及DPSO与LDPSO都比未经过Lévy flight的算法得到了更好的结果.其中，师学明(2009)的惯性因子阻尼振荡递减取值策略，在加入Lévy flight之后取得了最为理想的收敛效果，其迭代误差最小，并且也是最快达到收敛的算法，仅用了不到10次就搜索到全局极值点.原因在于，在引入Lévy flight之后，算法在每次迭代都有机会进行“大步长”的搜索，如果搜索效果不好则待在原位置，如果“有幸”得到更好的解，则会加快算法的收敛速度，因此可以看出，Lévy flight具有提高算法收敛速度的特性.

图3 (a)表示迭代50次，绝对误差随迭代次数的变化； (b)可以看出6种算法最终都收敛到全局极值

Schaffer函数是一种常见的用来检测优化算法全局收敛能力的一类特殊函数，其公式如下：

xi∈[-5.12,+5.12]

(13)

其图像如图4所示.该函数特点就是在中心点(x1=x2=0处)存在全局最大值f(x)=1，而在其周围存在一圈次极值，并且次极值约为0.99，与极值相差甚小，且全局最优解周围斜率较大，反演结果很容易越过最优极值，落在次级值上.因此该函数适合测试算法的收敛能力以及全局搜索能力.依旧以三种PSO算法对比加入Lévy flight之后的效果.由于其随机性，进行了大量的测试.

图4 Schaffer函数图像

图5a展示的是6种算法的迭代误差的递减，可以看到LDPSO最终误差最小，对应图5b，前5种算法最终都收敛到了次极值的圈脊上(两个绿色的点代表PSO与LPSO，两个蓝色的点代表MPSO与LMPSO，黄色的点代表DPSO，黑色的点代表LDPSO)，只有LDPSO收敛到了全局最优点，精度达到了足以跳出次极值点的精度.这归功于加入Lévy flight之后，使得算法在搜索到局部解之后，仍然具备一定全局搜索的能力，从而最终找到最优的全局解.可见LDPSO算法具有更好的全局搜索能力.

图5 6种PSO算法的计算结果

2.2 瞬变电磁L-PSO算法反演

为探究引入Lévy flight策略之后的粒子群算法在瞬变电磁数据反演方面的效果，分别对三种不同惯性参数策略的粒子群算法：PSO、MPSO、DPSO以及其引入Lévy flight之后的改进算法：LPSO、LMPSO、LDPSO，共6种算法进行同类对比.并且与另一种常用的随机算法——模拟退火算法(SA)进行对比，比较L-PSO算法与经典随机优化算法的区别.

反演数据由正演模拟得到，模拟装置采用40 m×40 m的单匝方形中心回线装置，发射电流为10 A，根据一维层状模型的解析解(Nabighian, 1992)，首先计算得到频率域的解(牛之琏, 2007)，然后通过Euler方法进行时-频转化，得到时间域瞬变电磁响应(Li et al., 2016; 王鹏飞等, 2019).在Matlab2014a平台下编程实现一维正演与本文提及的反演方法，处理器为Inter Core i5-6300HQ，主频为2.3 Hz(双核四线程).

(1)PSO及其改进型算法的同类对比

设置含有低阻与高阻夹层的7层复杂地层作为理论模型，其地层参数如下：

ρ1=50 Ωm，ρ2=10 Ωm，ρ3=100 Ωm，ρ4=200 Ωm，ρ5=100 Ωm，ρ6=600 Ωm，ρ7=200 Ωm；h1=50 m，h2=5 m，h3=50 m，h4=50 m，h5=100 m，h6=10 m.

分别采用：PSO、LPSO、MPSO、LMPSO、DPSO、LDPSO六种粒子群算法进行反演.其中DPSO采用ω=0.99k·rand()/2+0.1的策略，其他参数统一设置为：c1=c2=2，粒子数W=30.迭代次数上限50次，误差RMS设置为3%.为了与真实模型进行更准确的对比、方便计算拟合误差，反演模型层数与真实模型层数相同，并且其初始模型范围相同，经过5次试验反演之后，平均结果如表1及图6所示.

图6第一行的三张图为采用不同惯性参数的PSO算法，第二行分别对应其引入Lévy flight之后的改进型算法.从表1的统计结果来看，所有方法最后都得到了收敛，误差在3%以下.首先从惯性参数的取值策略来看，DPSO的惯性参数策略好于PSO和MPSO.根据图6的结果，在6种算法中，明显可以看到LDPSO算法的拟合结果最优，并且相对于3种未采用Lévy flight的PSO算法，在其引入Lévy flight之后，拟合结果都得到了相应的提升，拟合误差进一步减小.

图6 六种PSO算法反演得出的模型拟合图

表1 六种PSO算法反演对比结果

对5次结果的统计表明LDPSO算法对理论模拟地层的反演精度明显优于其他PSO算法，尤其是对薄层的反演上具有更好的反演效果.值得提及的是，在这10次反演中，LDPSO算法都只用数次迭代就达到了误差精度内，展现了LDPSO算法高效的特点.

(2)与模拟退火算法对比

粒子群算法的同类算法对比表明LDPSO算法具有更优秀的性能，因此将LDPSO算法与经典模拟退火(SA)算法作对比，以四层地层模型：典型KH型和低阻薄层地层模型为例.设置参数如下：

KH型：ρ1=10 Ωm，ρ2=100 Ωm，ρ3=20 Ωm，ρ4=400 Ωm，h1=50 m，h2=100 m，h3=100 m；

低阻薄层：ρ1=50 Ωm，ρ2=10 Ωm，ρ3=200 Ωm，ρ4=500 Ωm，h1=50 m，h2=10 m，h3=100 m.

两种算法的反演结果如图7.

图7 (a) LDPSO算法反演KH型地层模型的拟合图，(b) SA算法反演KH型地层模型的拟合图

对于KH模型，为更合理的评估两种算法和真实模型的拟合误差，反演模型层数与真实模型采用相同层数，两种算法采用相同的最大迭代次数，并且初始模型相同，反演模型参数ρ1为0～100 Ωm，ρ2为10～300 Ωm，ρ3为0～100 Ωm，ρ4为100～1000 Ωm；h1为0～100 Ωm，h2为10～200 Ωm，h3为10～200 Ωm.由上述反演拟合结果来看，在相同迭代次数的情况下，LDPSO算法效果明显优于经典SA算法.

图8是低阻薄层的反演结果.两种算法同样具有相同的最大迭代次数，迭代次数上限为100次，反演模型参数设置如下：ρ1为0～150 Ωm，ρ2为0～100 Ωm，ρ3为50～600 Ωm，ρ4为100～1000 Ωm；h1为0～100 Ωm，h2为0～100 Ωm，h3为10～200 Ωm.由图8可以看出，虽然两种算法对低阻薄层都有所反映, 但LDPSO算法反演的薄层厚度比SA算法效果更为精确.

图8 (a) LDPSO算法反演低阻薄层地层模型的拟合图，(b) SA算法反演低阻薄层地层模型的拟合图

3 实测数据试验

为验证L-PSO算法的实际应用效果，结合陕西韩城某煤矿井采空区瞬变电磁调查数据进行测试，并与常用的瞬变电磁反演方法进行对比，分析反演效果.

3.1 地质背景

矿井位于陕西省韩城市东北部，矿井已经开采完毕，主采1#煤层，煤层平均厚度6.1 m，采煤高度2.3～2.5 m，放煤高度3.6～3.8 m，主采煤层海拔约为+230至+240 m.矿区地表为典型黄土高原地貌，沟壑纵横，平均海拔约为+450 m，地表平均高差约为150 m.浅层以黄土层覆盖，经受强烈的剥蚀和地表水的长期冲蚀切割，下切很深.中深部因构造应力较弱，地质构造相对比较简单，除宽缓的褶曲外，岩层倾角平缓，一般10°以下.矿井工作面处于地层走向方向上的宽缓褶曲构造的背斜区域，同时距离边浅部复杂地质构造带较远.在该井边界内及其附近煤矿，发育中、小型断层和褶曲构造9条，落差3～9 m不等，造成该地区采空区容易充水，为探明该矿井积水范围和采空区分布特征，在地面进行了瞬变电磁法探测.

3.2 数据采集与预处理

现场数据采集仪器使用重庆璀陆探测技术有限公司的FCTEM-60瞬变电磁探测系统，该套采集系统采用重叠小线圈收发装置，高密度时序序列数据采样方式，横向分辨率高.现场施工严格按照地面瞬变电磁施工规范执行.瞬变电磁发射线框为半径为1.55 m的多匝圆形线框，发射电压为21 V，点距2 m.本地区实际完成45条测线，实际探测瞬变电磁物理点1866个，测线分布见图9.

图9 工区测线分布图.由于地表沟壑纵横，山体起伏较大且陡峭,很难布设贯穿测区的南-北走向的较长测线

3.3 反演

预处理之后，根据瞬变电磁法的电磁场扩散规律，用“烟圈法”进行半定性的解释.根据“烟圈法”和已掌握的地质资料设置反演初始模型，在反演过程中设置30个粒子数，最大迭代次数为50次，反演采用10层模型，为显示反演效果，选取了常用的阻尼最小二乘算法进行对比，图10将阻尼最小二乘算法和LDPSO算法反演的模型参数响应与实测数据的拟合情况同时展示.

图10 (a) 中左图为阻尼最小二乘算法反演拟合结果，右图是LDPSO优化算法反演拟合结果； (b) 是两种算法在50次迭代中的误差变化曲线

由拟合结果图10a和误差变化折线图10b的对比可以看出，对于瞬变电磁法早期和晚期数据差异较大的数据，LDPSO算法的反演效果明显优于阻尼最小二乘算法，并且在50次的迭代过程中，尽管LDPSO在初始位置的解相较于最小二乘方法更差，但其在初始的几步迭代中得到迅速收敛，到迭代晚期搜寻到了更好的解，再次印证了该方法在瞬变电磁实测资料的一维反演中的适用性.

3.4 资料解释

我们采用LDPSO算法，对所有测线数据进行了反演，根据测线的实际坐标和海拔标高，绘制了几条主要的测线剖面在测区的实际位置(如图11)；并根据反演电阻率的相对值，估计采空区电阻率范围勾勒出测区的三维异常图(如图12)，以及反演目标层范围内的异常体平面分布(如图13).

图11 左图表示工区内穿过异常体的几条测线反演的电阻率剖面在空间的位置，右图为这几条测线的俯视图

图13 测区异常体分布与塌陷规律对比图，黄色线段表示推测异常区在空间上的分布联系

从反演电阻率剖面切片图11可以看出，通过异常体的几条主要测线，反演的异常体均在煤层所在深度，并且各个测线之间的反演结果在平面上能够相互对应，具有一致性.从图12的异常体分布范围可以看出，低阻异常体呈现条带状分布，符合采空区的空间分布，且异常体的高度约为15 m，根据煤矿采空区一般的垮落规律，采空区冒落带高度一般为采煤高度的3倍左右.因此推测低阻异常体是由于采空区积水造成的，而反演目标深度范围内的高阻异常体，推断为由塌陷引起的主要空洞、裂隙、破碎带非积水区.

通过掌握的已有地质资料结合瞬变电磁反演结果，我们认为，由于该工作面已经停产了相当的一段时间，地表沟壑纵横，且土质疏松，由于采空区塌陷影响和原本的裂隙发育，形成渗水通道，易于大气降水及上部含水层中水的补给和储存，造成了采空区积水.其次，根据图13的异常体分布范围，可以看出，工区东西两侧的采空区范围明显比中间区域的采空区范围大，这是由于单一工作面塌陷时，靠近中心位置的采空区越塌实，不易形成积水空间，而越靠近两侧时采空区塌陷情况复杂，不能完全塌实，因此形成了可能的积水空间.通过本次探测的结果，该区的塌实为“畸变”的OX型结构，从整体上符合采空区塌陷的一般规律.

4 结论与建议

数值模拟及实测数据的反演案例表明常用的粒子群类算法在引入Lévy flight之后，在全局搜索能力上具有更出色的能力.其具有两个优点：首先，在保留了经典PSO算法收敛的快速性的同时，提高了算法游走全局的能力，大大增加了算法搜索到全局最优的概率；其次，由于Lévy flight重尾的特点，在早期加快了PSO算法的收敛速度，而在极值附近保持了PSO算法的性能.相比蒙特卡罗类算法和模拟退火等优化算法，L-PSO算法调用参数较少.但值得注意，虽然Lévy flight增加了PSO算法在迭代的任何时刻都具备探索全局的能力，但面对单峰函数，Lévy flight不仅不会提高PSO算法的性能，反而会因为多一步搜索而略微降低收敛速度.因此，我们认为Lévy flight更适用于多极值的非线性问题.

L-PSO算法相较于其他改进型的PSO算法，其多出了一步游走的步骤，因此多调用了一次正演，使得速度方面有所降低，对于瞬变一维反演来说，时间增加并不明显，但是对于二、三维反演，瞬变正演过程耗时巨大，每多一次正演对反演的效率都造成很大影响，因此，如何将L-PSO算法应用于高维瞬变电磁反演是下一步需要研究的问题.

致谢感谢重庆璀陆探测技术有限公司提供的仪器，感谢审稿人给出的建设性意见和建议.