CN110764154A - 基于改进粒子群算法的时间域航空电磁一维反演方法 - Google Patents

Abstract

本发明涉及一种基于改进粒子群算法的时间域航空电磁一维反演方法，包括以下步骤：a、输入含激发极化效应的时间域航空电磁观测数据；b、设置种群个数，阈值，最大迭代次数，随机产生种群各粒子初始位置和速度；c、计算各个粒子的正演响应值，以及每个粒子的RMS值；d、计算各个粒子的适应度值fitness，并按照从大到小顺序排序；e、判断是否满足反演终止条件，若满足则输出最优粒子，迭代结束，否则进行步骤f；f、对不同fitness值的粒子采用相应的更新粒子方式，产生新的粒子；g、返回步骤c。本发明采用的改进粒子群算法不仅具有不易陷入局部极值，不依赖初始模型的优点，而且大大改善了经典粒子群算法存在的收敛速度慢，易陷入“早熟”的问题，将其应用到含激发极化效应的时间域航空电磁法中，可得到良好的反演效果。

Description

基于改进粒子群算法的时间域航空电磁一维反演方法

技术领域

本发明属于地球物理勘探技术领域，具体涉及一种时间域航空电磁反演方法，特别涉及一种基于改进粒子群算法的含激发极化效应的时间域航空电磁一维反演方法。

背景技术

航空电磁法有速度快、成本低和无需地面人员操作等优点，特别适合于地形条件复杂的无人区，如高山、沙漠、沼泽、森林、湖泊等进行大面积勘查工作，在矿产探测、地质普查等方面具有广泛的应用。近些年，时间域航空电磁法勘探在观测数据中常常出现符号反转现象，给时间域航空电磁法的解释带来了一定的困扰。这种符号反转现象通常认为是由于地下介质存在激发极化效应(Induced Polarization Effect，简称IP效应)引起的。针对时间域航空电磁法观测数据中出现符号反转现象，对含有IP效应的观测数据进行考虑激发极化效应的反演，能够解决符号反转带来的解释困扰，从而更真实的反映地下地电参数的分布。

含有IP效应的时间域航空电磁数据是由地下各地层的电阻率，极化率，时间常数，频率相关系数共同决定的，由于其反演参数较多，而且各个参数灵敏度差异较大，往往具有严重的多解性，给反演带来了巨大挑战。传统的梯度类反演方法，如Occam，GN，NLCG等方法以收敛速度快，反演精度高等优势广泛应用于地球物理反演中，然而，这类方法依赖于选择的初始模型，容易陷入局部极值，针对含有IP效应的时间域航空电磁数据，由于具有严重的多解性，很难得到有效的反演结果。

粒子群优化算法(Particle Swarm Optimization，PSO)是近年来发展起来的一种全局性的群体智能随机搜索算法，最早由Kennedy and Eberhart于1995年提出，通过对简单社会系统的模拟，在多维解空间中构造“粒子群”，把待求解问题潜在的可行解看作群体中的每个粒子，并赋予与目标函数相关的适应度值。群体中每个粒子通过跟踪自己和群体所发现的最优值，修正粒子飞行的速度，从而带领整个粒子群体寻优，实现对问题最优解的搜索。虽然，经典粒子具有全局搜索，不受初始模型影响的优势，但该方法收敛速度慢，且易出现“早熟”现象，应用于含有IP效应的时间域航空电磁数据反演仍十分困难。

发明内容

本发明的目的在于针对上述现有技术的不足，提供一种基于改进粒子群算法的含激发极化效应的时间域航空电磁一维反演方法，在经典粒子群算法的基础上引进Occam的搜索方式，在保证全局搜索的基础上加快局部搜索速度，大大提高了反演速度；此外，在每次迭代过程中，选出一部分适应度值较差的粒子在全局最优粒子附近扰动，生成的一部分新的粒子，对其余的适应度值较差的粒子随机变异生成新的粒子，从而保持群体的多样性，避免算法出现“早熟”现象。

本发明的目的是通过以下技术方案实现的：

一种基于改进粒子群算法的时间域航空电磁一维反演方法，包括以下步骤：

a、输入含激发极化效应的时间域航空电磁观测数据；

b、设置种群个数np，阈值eps，最大迭代次数t_max，设置模型参数的上下限，在模型参数上下限范围内随机生成初始粒子位置x_i(0)和速度v_i(0)；

其中，x_i(0)和速度v_i(0)代表第0次迭代的第i个粒子位置和速度；

c、正演计算，对种群中每个粒子计算其含激发极化效应的时间域电磁响应值，并计算每个粒子的RMS值；

d、计算各个粒子的适应度值fitness，这里，fitness＝1/φ_i，φ_i为各个粒子建立的反演目标函数，并对fitness按从大到小排序；

e、判断适应度值RMS是否小于设定的阈值eps，以及迭代次数t是否等于t_max，若RMS<eps或者t＝t_max，则输出最优粒子，迭代结束，否则进行步骤f；

f、对fitness排在前50％的粒子，以Occam的搜索方式更新粒子位置，对fitness排在前50％-60％的粒子，以经典粒子群算法的搜索方式更新粒子位置，对fitness排在前60％-70％，70％-80％，80％-90％的粒子，在历史全局最优值附近分别以20％，50％，80％的扰动幅度随机变异生成新的粒子位置，对fitness排在90％-100％的粒子，随机变异生成新的粒子位置；

g、返回步骤c。

与现有技术相比，本发明的有益效果在于：

本发明采用的改进粒子群算法不仅具有不易陷入局部极值，不依赖初始模型的优点，而且大大改善了经典粒子群算法存在的收敛速度慢，易陷入“早熟”的问题。本发明对部分粒子采用Occam的搜索方式，相较于传统的粒子群算法，大大提高了收敛速度；每次迭代过程中选出一部分适应度值较差的粒子在全局最优粒子附近扰动，生成的一部分新的粒子，同时随机变异一部分适应度值较差的粒子，从而保持群体的多样性，避免粒子陷入“早熟”。改进的粒子群算法大大提高了经典粒子群算法的优化能力，将其应用到含IP效应的时间域航空电磁反演中，可得到良好的反演效果。

附图说明

图1基于改进粒子群算法的含激发极化效应的时间域航空电磁一维反演方法流程图；

图2航空电磁中心回线装置示意图；

图3a三层地质模型电阻率反演结果图；

图3b三层地质模型极化率反演结果图；

图4a四层地质模型电阻率反演结果图；

图4b四层地质模型极化率反演结果图。

具体实施方式

下面结合附图和实施例对本发明作进一步的详细说明。

本发明是在经典粒子群算法的基础上进行改进，将其应用于含激发极化效应的时间域航空电磁数据一维反演中，方法流程如图1所示，包括以下步骤：

a、输入含激发极化效应的时间域航空电磁观测数据。

b、设置种群个数np，阈值eps，最大迭代次数t_max，设置模型参数的上下限，在模型参数上下限范围内随机生成初始粒子位置x_i(0)和速度v_i(0)。

经典粒子群可以用以下方程来描述：

v_ij(t+1)＝ω(t)·v_ij(t)+c₁·r_1j(t)·(p_ij(t)-x_ij(t))+c₂·r_2j(t)·(p_gj(t)-x_ij(t))， (1)

x_ij(t+1)＝x_ij(t)+v_ij(t+1)， (2)

ω(t)＝ω_max-(ω_max-ω_min)·t/t_max， (3)

其中，下标i表示第i个粒子，j表示第j个参数；x_i＝[x_i1,x_i2,...,x_iN]表示第i个粒子的位置，这里N代表模型参数的总个数，将地下分为n层，每层地层均含有Cole-Cole模型的四种参数：零频电阻率ρ_∞，极化率η，时间常数τ，频率相关系数c，则粒子x_i＝[ρ_∞1,ρ_∞2,…,ρ_∞n；η₁,η₂,…,η_n；τ₁,τ₂,…,τ_n；c₁,c₂,…,c_n]^T，共有4n个模型参数；v_i＝[v_i1,v_i2,...,v_iN]表示粒子更新速度；t表示第t次粒子群迭代；p_i＝[p_i1,p_i2,...,p_iN]表示第i个粒子的历史个体最优值；p_g＝[p_g1,p_g2,...,p_gN]表示群体的历史全局最优值；r₁、r₂为[0，1]间相互独立的随机函数值；ω为惯性权重系数，决定上次迭代速度所占的权重，ω_max取值为0.9，ω_min取值为0.4；c₁、c₂为学习因子，分别影响“自我学习”能力和“社会学习”能力，一般取值为2。

模型参数的上下限是指四种模型参数的范围，即零频电阻率的下限ρ_∞min，上限ρ_∞max，极化率的下限η_min，上限η_max，时间常数的下限τ_min，上限τ_max，频率相关系数的下限c_min，上限c_max。

随机生成初始粒子位置x_i(0)和速度v_i(0)，即随机生成第0次迭代的第i个粒子的位置和速度，这里，一共生成np个粒子位置和粒子速度。

每个粒子x_i均可构建相应的反演目标函数，其形式可表示为：

φ_i＝(Rx_i)^TRx_i+μ^-1[W(d-F(x_i))]^TW(d-F(x_i))， (4)

式中，右端第一项是模型粗糙度项，反映模型的不光滑程度，R是模型约束矩阵；第二项是数据拟合项，反映模型参数正演响应与观测数据的拟合程度，d＝[d₁,d₂,...,d_M]表示观测数据，M代表观测数据的个数，本文d取各个时间道的含激发极化效应的时间域电磁响应值，W＝[1/e₁,1/e₂,...,1/e_M]表示对角加权矩阵，这里e代表观测误差，F是正演算子；μ是拉格朗日乘子，用来权衡模型粗糙度项和数据拟合项的权重。

数据拟合差可以表示为：

RMS可以判断数据拟合的程度，其理论最佳值为1.0，阈值eps是设定的终止条件。

c、正演计算，对种群中每个粒子计算其含激发极化效应的时间域电磁响应值，包括以下步骤：

频率域中准静态条件下，水平圆形回线源产生的磁场强度z分量可表示为：

这里，x、y轴位于水平地表面，取z轴垂直向下为正，水平圆形回线源位于直角坐标系原点正上方h处，圆形回线发射电流强度为I，回线半径为a，r_TE为反射系数，J₁(λa)为一阶贝塞尔函数。式(6)的结果可通过汉克尔积分计算得到。

频率域响应经过正余弦变换，并根据正余弦函数和贝塞尔函数间的转换关系式，可以得到含1/2阶贝塞尔函数的积分。通过汉克尔积分计算，便可得到时间域响应值，再通过殷长春等(2013.时间域航空电磁系统瞬变全时响应正演模拟.地球物理学报，2013，56(9)：3153-3162)提出的波形褶积算法便可得到任意波形下的时间域航空电磁响应。

要想计算含激发极化效应的时间域航空电磁响应，需要使用复电阻率。Cole–Cole模型下的复电阻率表达式如下：

其中，ρ₀表示零频电阻率，η为极化率，τ为时间常数，c表示频率相关系数，i为虚数单位，ω为角频率。零频电阻率ρ₀与无穷频电阻率ρ_∞存在如下关系式：

ρ₀＝ρ_∞(1-η)， (8)

因此，式(7)可以转化为以下形式：

公式(9)中得到的是含极化率的复电阻率，将该复电阻率代替时间域航空电磁响应计算公式中的实电阻率，即可得到含有激发极化效应的时间域航空电磁响应值。

依照上述过程，可计算每个粒子含激发极化效应的时间域电磁响应值，之后，利用公式(5)可计算每个粒子的RMS值。

d、计算各个粒子的适应度值fitness，并按从大到小排序。适应度值fitness主要用来描述粒子的好坏，这里，用反演目标函数的倒数来描述：

fitness＝1/φ_i (10)

e、判断适应度值RMS是否小于设定的阈值eps，以及迭代次数t是否等于t_max，若RMS<eps或者t＝t_max，则输出最优粒子，迭代结束；否则进行步骤f。

f、对fitness排在前50％的粒子，以Occam的搜索方式更新粒子位置，对fitness排在前50％-60％的粒子，以经典粒子群算法的搜索方式更新粒子位置，对fitness排在前60％-70％，70％-80％，80％-90％的粒子，在历史全局最优值附近分别以20％，50％，80％的扰动幅度随机变异生成新的粒子位置，对fitness排在90％-100％的粒子，随机变异生成新的粒子位置；

经典PSO算法尽管具有全局搜索，不受初始模型影响的优势，然而也有一定的局限性：

1)算法后期收敛速度较慢，计算量较大；

2)算法容易出现“早熟”现象，容易陷入局部极值；

针对经典PSO算法的局限性，本申请针对经典PSO算法进行适当的改进.：

1)引入Occam反演算法中的搜索方式提高PSO算法每次迭代的收敛速度；

2)选出一部分适应度值较差的粒子在全局最优粒子附近以一定幅度扰动，生成新的粒子；

3)变异一部分适应度值较差的粒子，生成新的粒子，从而保持群体的多样性。

由于经典粒子群算法后期收敛速度较慢，计算量较大，本申请将Occam的搜索方式引进粒子群算法中，以Occam的搜索方式更新粒子位置，其迭代方程表达形式如下：

其中，

k代表第k次粒子群迭代，μ的取值按照Zhdanov(2009.GeophysicalElectromagnetic Theory and Methods.New York:Elsevier.)自适应算法求取，即

μ_k＝μ₀q^k-1，q>0，k＝1，2，3，…，K (13)

初始值μ0由数据项与模型项的比值给定，定义相邻两次迭代数据拟合项相对变化量为

Δφ＝(||W(d-F(x_i(k)))||²-||W(d-F(x_i(k+1)))||²)/||W(d-F(x_i(k))||²， (14)

式中i为迭代次数，衰减因子q的大小按照Δφ的大小去调整：

对fitness排在前50％的粒子，以Occam的搜索方式更新粒子位置，即采用公式(11)-(15)的形式更新粒子位置。

对fitness排在前50％-60％的粒子，以经典粒子群算法的搜索方式更新粒子位置，即采用公式(1)-(3)的形式更新粒子位置。

经典粒子群算法搜索全局最优值的时，由于所有粒子都要向历史全局最优值学习，从而不断降低粒子群体的多样性，容易陷入局部极值，出现“早熟”现象。针对这个问题，本申请提出在历史全局最优粒子的每个参数上进行一定幅度的扰动生成部分新粒子，其表达式如下：

x_ij(t+1)＝(1-percent×r_ij(t))p_gj(t)， (16)

其中，r_ij(t)为[-1，1]间相互独立的随机值，percent为扰动幅度。

通过扰动生成新的粒子，能够增加粒子群的群体多样性，一定程度上避免出现“早熟”现象。

对fitness排在前60％-70％，70％-80％，80％-90％的粒子，在历史全局最优值附近分别以20％，50％，80％的扰动幅度随机变异生成新的粒子位置，即为了增加群体多样性，在历史全局最优值附近分别以20％，50％，80％的扰动幅度随机生成新的粒子，更新适应度值排在60％-70％，70％-80％，80％-90％的粒子。

为进一步增加粒子群的多样性，避免群体出现“早熟”现象，在每次迭代过程中，变异一部分适应度值较差的粒子，从而增加粒子群的多样性。对于选取变异的个体数量，选择太多个体可能影响算法的性能，选择太少个体可能无法明显增加群体的多样性，经试验，本申请选择总个数10％的粒子较为合适。

对fitness排在前90％-100％的粒子，随机变异生成新的粒子位置，即选取竞争力较差的，适应度值排在90％-100％的粒子进行随机变异，重新生成新的粒子。

g、返回步骤c。

实施例1

设计理论模型验证本发明的反演效果。

理论模型试算采用时间域航空电磁系统装置如图2所示，发射线圈半径为15m，发射电流峰值为870A，发射磁矩峰值为615,000Am²，发射波形为负阶跃波，接收机位于发射线圈的中心处，发射线圈和接收机距离地面的高度均为30m。首先设计理论模型，通过正演计算得到其响应数据，再对响应数据添加2％的随机噪声合成理论数据d，下面将理论数据作为观测数据进行反演试算。

由于时间域航空电磁响应对四个Cole-Cole模型参数的灵敏度相差较大，且部分参数之间还存在一定的相关性，导致完全反演四个Cole-Cole模型参数十分困难。由于时间域航空电磁响应对电阻率和充电率的灵敏度较大，对时间常数和频率相关系数的灵敏度相对较小，故本申请在反演过程中，将时间域航空电磁响应对时间常数和频率相关系数的导数强制为零，从而使电阻率和充电率两个参数的反演结果更加稳定，也大大降低了反演结果的多解性。

首先，讨论三层模型的反演结果。设计模型Model1：第一层电阻率为200Ω.m，极化率为0.1，厚度为60m；第二层为电阻率为20Ω.m，极化率为0.5的低阻高极化层，厚度为30m；第三层基底层电阻率为500Ω.m，极化率为0.2。各层的时间常数均为0.01s，频率相关系数均为0.3。将地下分成35层，首层厚度为2m，其他层厚按照1.05倍向下递增。

实施基于改进粒子群算法的含激发极化效应的时间域航空电磁数据一维反演，包括以下步骤：

a、输入含激发极化效应的时间域航空电磁观测数据。

b、设置种群个数np＝100；阈值eps＝1.02；最大迭代次数t_max＝200；设置模型参数的上下限，即ρ_∞min＝1，ρ_∞max＝2000，η_min＝0.01，η_max＝0.99，τ_min＝0.0001，τ_max＝1000，c_min＝0.1，c_max＝0.8；在模型参数上下限范围内随机生成初始粒子位置x_i(0)和速度v_i(0)。

c、正演计算，对种群中每个粒子计算其含激发极化效应的时间域电磁响应值，并利用公式(5)计算每个粒子的RMS值。

d、按照公式(10)计算各个粒子的适应度值fitness，并按从大到小排序。

e、判断适应度值RMS是否小于设定的阈值eps，以及迭代次数t是否等于t_max，若RMS<eps或者t＝t_max，则输出最优粒子，迭代结束，否则进行步骤f。

f、对fitness排在前50％的粒子，依据公式(11)-(15)，以Occam的搜索方式更新粒子位置；对fitness排在前50％-60％的粒子，依据公式(1)-(3)，以经典粒子群算法的搜索方式更新粒子位置；对fitness排在前60％-70％，70％-80％，80％-90％的粒子，依据公式(16)，在历史全局最优值附近分别以20％，50％，80％的扰动幅度随机变异生成新的粒子位置，对fitness排在前90％-100％的粒子，随机变异生成新的粒子位置。

g、返回步骤c。

Model1合成理论数据的反演结果如图3所示。从图3中可以看出：电阻率的反演结果与真实值十分接近，极化率的反演结果略差于电阻率的反演结果，这主要是因为极化率的灵敏度相较于电阻率小很多，从而影响极化率的反演结果；综合来看，电阻率和极化率的反演结果均能反映地下模型参数的变化趋势，可以得到较好的反演效果。

实施例2

为进一步验证算法的效果，设计一个四层模型Model2。Model2第一层电阻率为50Ω.m，极化率为0.1，厚度为22m；第二层为电阻率为300Ω.m，极化率为0.3的高阻高极化层，厚度为61m；第三层为电阻率为10Ω.m，极化率为0.6的低阻高极化层，厚度为42m；第四层基底层电阻率为100Ω.m，极化率为0.1；各层的时间常数均为0.01s，频率相关系数均为0.3。将地下分成35层，首层厚度为2m，其他层厚按照1.05倍向下递增。

实施基于改进粒子群算法的含激发极化效应的时间域航空电磁数据一维反演，包括以下步骤：

a、输入含激发极化效应的时间域航空电磁观测数据。

b、设置种群个数np＝100；阈值eps＝1.02；最大迭代次数t_max＝200；设置模型参数的上下限，即ρ_∞min＝1，ρ_∞max＝2000，η_min＝0.01，η_max＝0.99，τ_min＝0.0001，τ_max＝1000，c_min＝0.1，c_max＝0.8；在模型参数上下限范围内随机生成初始粒子位置x_i(0)和速度v_i(0)。

c、正演计算，对种群中每个粒子计算其含激发极化效应的时间域电磁响应值，并利用公式(5)计算RMS值。

d、按照公式(10)计算各个粒子的适应度值fitness，并按从大到小排序。

e、判断适应度值RMS是否小于设定的阈值eps，以及迭代次数t是否等于t_max，若RMS<eps或者t＝t_max，则输出最优粒子，迭代结束，否则进行步骤f。

f、对fitness排在前50％的粒子，依据公式(11)-(15)，以Occam的搜索方式更新粒子位置；对fitness排在前50％-60％的粒子，依据公式(1)-(3)，以经典粒子群算法的搜索方式更新粒子位置；对fitness排在前60％-70％，70％-80％，80％-90％的粒子，依据公式(16)，在历史全局最优值附近分别以20％，50％，80％的扰动幅度随机变异生成新的粒子位置，对fitness排在前90％-100％的粒子，随机变异生成新的粒子位置。

g、返回步骤c。

Model2合成理论数据的反演结果如图4所示。从图4中可以看出：与Model1反演结果类似，电阻率的反演结果与真实值十分接近，极化率的反演结果略差于电阻率的反演结果，这主要还是因为极化率的灵敏度相较于电阻率小很多，从而对极化率的反演结果有一定影响；综合来看，电阻率和极化率的反演结果均能反映地下模型参数的变化趋势，可以得到较好的反演效果。

Claims (1)

1.一种基于改进粒子群算法的时间域航空电磁一维反演方法，包括以下步骤：

a、输入含激发极化效应的时间域航空电磁观测数据；

b、设置种群个数np，阈值eps，最大迭代次数t_max，设置模型参数的上下限，在模型参数上下限范围内随机生成初始粒子位置x_i(0)和速度v_i(0)；

其中，x_i(0)和速度v_i(0)代表第0次迭代的第i个粒子位置和速度；

c、正演计算，对种群中每个粒子计算其含激发极化效应的时间域电磁响应值，并计算每个粒子的RMS值；

d、计算各个粒子的适应度值fitness，这里，fitness＝1/φ_i，φ_i为各个粒子建立的反演目标函数，并对fitness按从大到小排序；

e、判断适应度值RMS是否小于设定的阈值eps，以及迭代次数t是否等于t_max，若RMS<eps或者t＝t_max，则输出最优粒子，迭代结束，否则进行步骤f；

f、对fitness排在前50％的粒子，以Occam的搜索方式更新粒子位置，对fitness排在前50％-60％的粒子，以经典粒子群算法的搜索方式更新粒子位置，对fitness排在前60％-70％，70％-80％，80％-90％的粒子，在历史全局最优值附近分别以20％，50％，80％的扰动幅度随机变异生成新的粒子位置，对fitness排在90％-100％的粒子，随机变异生成新的粒子位置；

g、返回步骤c。

Images