CN110133733B - 一种基于粒子群优化算法的电导-极化率多参数成像方法 - Google Patents

Abstract

本发明涉及一种基于粒子群优化算法的电导‑极化率多参数成像方法，目的在于同时获取地下极化介质的导电和极化信息，提高成像精度。本发明基于柯尔‑柯尔模型，获得了梯形波激励下的感应‑极化磁场响应表达式，根据理论模型的磁场响应与实测磁场响应构建目标函数，利用超导时域电磁测量的早期、晚期数据分别获得零频电导率和时间常数，二者对粒子群搜索范围进行约束，缩小搜索范围，最终采用粒子群优化算法提取零频电导率、时间常数、极化率、频散系数、深度。本发明相比传统电阻率成像方法，更符合实际极化介质频散传播规律，提高了电导率、极化率的成像精度。

Description

一种基于粒子群优化算法的电导-极化率多参数成像方法

技术领域

本发明涉及地球物理勘探领域，具体地来讲为一种基于粒子群优化算法的 电导-极化率多参数成像方法。

背景技术

在电法勘探领域，电阻率解释是一个非常重要的内容，通常将其用来估算 地下介质的分布。受激发极化效应的影响，岩(矿)石电阻率是一个与频率有 关的复数，因此在向地下供入稳定电流时，不仅能产生以二次涡流场为主的感 应场，也能产生极化场。极化场的存在，使得电磁响应产生“反号”现象，该“反 号”包含大量的极化信息。因此，研究极化现象，准确提取地下介质中的极化参 数，进行极化率、电导率-深度成像是十分必要的。

目前已经应用到实际电法勘探极化参数提取的方法中，最为广泛的是 Constable(1987)提出的奥克姆反演方法，但在存储反演时面临着灵敏度矩阵需 要多次正演的问题。闫国祥(2015)应用最小二乘法对极化参数进行了分步提 取，但结果过于依赖第一步反演结果，而且仅是理论模型研究，没有应用于实 际例子。

CN109376446A公开了一种基于粒子群优化的极化曲线拟合方法，特别是 涉及在电化学领域腐蚀动力学参数的提取，待拟合的动力学电流曲线呈e指数 非线性变化，但是其使用的是基本粒子群算法对整条曲线进行拟合，并没有对 特征区间进行划分，不足以描述极化特征，而且随着种群数目、待求参数量的 增大，时间成本也在累计增加，为此将粒子群算法进行优化是十分必要的。

CN101706587A公开了一种电法勘探激电模型参数的提取方法，特别是涉 及在柯尔-柯尔模型中的未知参数。该方法结合随机统计算法与最小二乘法，对 单个柯尔-柯尔模型相对相位谱和振幅联合提取极化效应参数，既避免了随机统 计算法中过于依赖初值的问题，又保留了最小二乘法的提取速度较快的特点。 但以上参数提取方法针对的仅是复电导率，在实际野外测量中，不能直接测得 复电导率值，更多的是得到电磁响应数据，因此严重影响参数提取在实际野外 中的应用，所以研究实际极化地质的极化参数提取方法具有重要意义。

发明内容

本发明的目的在于针对现有极化参数提取方法的不足，根据实际地质中存 在的极化问题，建立柯尔-柯尔模型，提供一种基于粒子群优化算法的电导-极 化率多参数成像方法。

本发明是这样实现的，

一种基于粒子群优化算法的电导-极化率多参数成像方法，该方法包括：

1)基于地下极化介质的柯尔-柯尔复电导率模型，获得梯形波激励下的时 域感应-极化磁场响应；

2)采用超导量子传感器时域电磁探测系统进行野外磁场数据观测，并对实 测数据进行叠加取样、滤波、以及一次场剔除处理；

3)利用步骤2)数据的早期数据计算零频电导率值、晚期数据e指数多项 式拟合得到时间常数，将获得的零频电导率值和时间常数对粒子群搜索范围进 行约束；

4)根据实测磁场数据与步骤1的时域感应-极化磁场响应构建目标函数， 利用粒子群优化算法，根据步骤3约束的范围，提取极化参数，包括零频电导 率σ₀、时间常数τ、极化率η、频散系数c、深度d；

5)利用步骤4提取到实测数据的极化参数，带入广义趋肤深度计算公式， 进行极化介质的电导率、极化率-深度成像。

进一步地，步骤3中，发射电流关断后，早期响应以感应场为主，利用早 期数据计算零频电导率值σ₁；晚期数据以极化场为主，对其晚期数据进行e指 数拟合，得到时间常数值τ₁；零频电导率值σ₁和时间常数值τ₁的±20％分别作为 零频电导率和时间常数值的粒子群搜索范围上下限。

进一步地，步骤4中具体包含以下步骤：

4Ⅰ、输入相关参数，随机初始化粒子的位置、速度和个体最优值P_best、种 群最优值G_best；

4Ⅱ、构建目标函数：minφ(X)＝φ₁(X)+λφ₂(X)，X＝σ₀、η、τ、c，

其中

φ₂(X)＝f_t(X)，f_t(X)为步骤1的时域感应 -极化磁场响应，B_t为实测磁场数据，λ为收敛因子；

4Ⅲ、计算每个粒子当前适应度，与P_best、G_best进行比较，若当前适应度值 较小，则令其替代P_best、G_best，反之继续保留P_best、G_best值；

4Ⅳ、实测曲线依据斜率划分为两段，据此分别按照式(1)和式(2)更新 粒子的速度和位置：

x_ij(k+1)＝x_ij(k)+v_ij(k+1)·t (2)

式中，

t₀为曲线斜率由负变为正的时刻，p_k为最大迭代次数，iter为当前迭代次数；

迭代过程中，当粒子的速度或位置超出搜索范围，产生变异时，根据二者 移动方向分为两种情况，若二者变异方向相反，则不做任何处理，反之，将由 下式对位置和速度重新定义，当超出上边界时：

v_ij(k+1)＝v_max-2·r₃·v_max (3)

x_ij(k+1)＝x_max-0.5·r₄·x_max (4)

当粒子的位置或速度超出下边界时：

v_ij(k+1)＝v_min-2·r₃·v_min (5)

x_ij(k+1)＝x_min-0.5·r₄·x_min (6)

4Ⅴ、判断是否达到最大迭代次数，若是，获得P_best、G_best的位置信息，输 出σ₀、τ、η、c、d；否则返回步骤4Ⅲ，继续执行。

进一步地，步骤5中包含以下步骤：

5Ⅰ、应用matlab编程，输入步骤4提取出的零频电导率σ₀、时间常数τ、 极化率η、频散系数c到程序中，根据

公式，获得复电 导率值；

5Ⅱ、根据广义趋肤深度公式，计算该复电导率σ(ω)的趋肤深度；

5Ⅲ、利用步骤5Ⅰ中的零频电导率、极化率和步骤5Ⅱ获得的趋肤深度， 分别进行电导率-深度、极化率-深度成像。

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

本发明的极化介质电导、极化率参数成像方法，相比传统电阻率成像方法， 更符合实际地下极化介质的传播扩散规律，从而提高了电导率探测深度的解释 精度。本方法为开展电磁探测寻找矿产资源提供新的技术保障，有利于瞬变电 磁勘探方法的实用化。

附图说明

图1是基于地下极化介质模型的电导率、极化率多参数成像方法流程图；

图2是粒子群优化算法提取极化参数的流程图；

图3是根据地质资料建立的三层层状大地模型；

图4是本发明一个实施例的层状模型提取结果成像对比图。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合实施例， 对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以 解释本发明，并不用于限定本发明。

本发明在电流关断后100ms内，采用超导量子传感器时域电磁探测系统观 测电磁数据，晚期数据采用e指数拟合确定时间常数τ1；根据发射电流计算一 次磁场，从实测数据中减去一次磁场后，计算零频电导率σ1；二者大小的±20％ 作为搜索范围输入粒子群中，根据理论模型的磁场响应与实测磁场响应构建目 标函数，最后采用粒子群优化算法进行零频电导率σ0、时间常数τ、极化率η、 频散系数c、深度d的提取。

本发明是这样实现的，一种基于粒子群优化算法的电导-极化率多参数成像 方法包括：

1)基于地下极化介质的柯尔-柯尔复电导率模型，获得梯形波激励下的时 域感应-极化磁场响应表达式；

先获取层状模型下负阶跃波的磁场响应，将梯形波发射电流一阶导数与磁 场响应卷积，可得到梯形波激励下的感应磁场响应；再带入柯尔-柯尔模型频域 表达式

其中σ₀为零频电导率、τ为时间常数、η为极 化率、c为频散系数，最后采用基于贝塞尔函数的快速Hankel变换、余弦变换 算法，将结果从频域转换到时域，从而获得层状大地模型下时域感应-极化双场 电磁响应，即理论模型磁场响应。

2)使用超导量子传感器时域电磁探测系统接收磁场，对电磁实测数据进行 处理，包括叠加取样、抽道、滤波、一次场剔除等处理；

3)处理后的数据，利用早期数据(一般为100微秒～1毫秒)计算零频电 导率值，利用晚期数据(一般二十个毫秒以后)e指数拟合得到时间常数，二者对 粒子群搜索范围进行约束；

电流关断后，晚期数据以极化场为主，对其进行e指数拟合，得到时间常 数值τ₁；早期响应以感应场为主，根据发射电流计算一次磁场，从实测数据中 减去一次磁场后，计算零频电导率σ₁。

4)根据实测磁场数据与理论模型磁场响应构建目标函数，利用粒子群优化 算法进行参数提取，包括零频电导率σ₀、时间常数τ、极化率η、频散系数c、 深度d；

基于粒子群优化算法的电导-极化率多参数提取，具体步骤如下：

4Ⅰ、输入粒子群优化算法相关参数D、p_k、b_k、

x_limit、v_limit，其中 D为问题维数，b_k为种群规模，p_k为最大迭代次数，

为加速度因子，x_limit为粒子的位置限制，v_limit为速度限制；随机初始化粒子的位置、速度、个体最 优值P_best和种群最优值G_best；

4Ⅱ、构建目标函数：

minφ(X)＝φ₁(X)+λφ₂(X)，X＝σ₀、η、τ、c，

其中

φ₂(X)＝f_t(X)，B_t为实测磁场数据，f_t(X)为理论模型磁场响应，λ为收敛因子；

4Ⅲ、计算每个粒子当前适应度，与P_best、G_best进行比较，若当前适应度值 较小，则令其替代P_best、G_best，反之继续保留P_best、G_best值；

4Ⅳ、实测曲线依据斜率划分为两段，据此更新粒子的速度(式(1))和 位置(式(2))：

x_ij(k+1)＝x_ij(k)+v_ij(k+1)·t (2)

式中，

t₀为斜率由负变为正的时刻，p_k为最大迭代次数，iter为当前迭代次数；下标i表示第i个粒 子(i＝1,2,3......M)，下标j表示粒子的第j维(j＝1,2,3......D)，i表示当前的迭代次 数，假设粒子群搜索空间为D维，粒子群有M个粒子(通常将粒子数设置为问 题维数的5～10倍)，r₁、r₂是(0，1)之间的随机数。

为了保证粒子在搜索空间内搜索，通常将粒子的速度限制在[-v_max v_max]区间 中，将步骤3得到的搜索范围作为粒子的位置限制[-x_min x_max]，则粒子的最大 速度可以取x_max的10％～20％。

由于粒子的速度和位置均有方向性，迭代过程中，当其到达边界，即产生 变异时，大致可分为两种情况，若二者变异方向相反，则不做任何处理，反之， 将由下式对位置和速度重新定义，当粒子的位置或速度超出上边界时：

v_ij(k+1)＝v_max-2·r₃·v_max (3)

x_ij(k+1)＝x_max-0.5·r₄·x_max (4)

当粒子的位置或速度超出下边界时：

v_ij(k+1)＝v_min-2·r₃·v_min (5)

x_ij(k+1)＝x_min-0.5·r₄·x_min (6)

式中r₃、r₄是(0，1)之间的随机数。

Ⅴ、判断是否达到最大迭代次数，若是，获得P_best、G_best的位置信息，输出 σ₀、η、τ、c、d；否则返回步骤Ⅲ，继续执行；

5)利用提取到的实测数据的极化参数，带入深度计算公式中，绘制极化介 质的电导率、极化率-深度图。包括

5Ⅰ、应用matlab编程，输入步骤4提取出的零频电导率、时间常数、充电 率、频散系数到程序中，带入

中获得复电导率值；

5Ⅱ、根据广义趋肤深度公式，计算该复电导率σ(ω)的趋肤深度；

5Ⅲ、利用步骤5Ⅰ中的零频电导率、极化率和步骤5Ⅱ获得的趋肤深度， 分别进行电导率-深度、极化率-深度成像。

实施例

结合图1所示，一种基于粒子群优化算法的电导-极化率多参数成像方法， 包括：

1)基于地下极化介质的柯尔-柯尔复电导率模型，获得了梯形波激励下的 时域感应-极化磁场响应表达式；

如图1所示的步骤中，先根据图4(a)中的地质资料，建立一个三层大地模 型，初始化各层电性参数，获取该层状大地模型的负阶跃波磁场响应，将梯形 波发射电流与磁场响应卷积，得到梯形波下的感应磁场响应；再带入柯尔-柯尔 模型表达式

最后采用基于贝塞尔函数的快速Hankel 变换、余弦变换算法，将结果从频域转换到时域，从而获得层状大地模型下时 域感应-极化双场电磁响应。

2)使用超导量子传感器时域电磁探测系统接收磁场，对电磁实测数据进行 处理，包括叠加取样、抽道、滤波、一次场剔除等处理；

进一步地，超导量子传感器时域电磁探测系统接收到的原始数据需要进行 转换及归一化，使原始数据变为有意义的电磁响应量。由于电磁勘探系统在实 际测量过程中，有很多因素使得数据产生噪声，为了提高数据质量，要对噪声 特点分析，并对实测数据进行噪声滤波，采用多次叠加消除漂移量。根据发射 电流计算一次场

其中R表达收发距，实测数据减去一次场。

3)利用早期数据计算零频电导率值，晚期数据e指数多项式拟合得到时间 常数，二者对粒子群搜索范围进行约束；

进一步地，晚期以极化放电为主，因此利用晚期数据采用e指数拟合确定 时间常数τ1＝0.0055ms，利用早期数据计算零频电导率σ1＝0.058S/m。

4)根据实测磁场数据与理论模型磁场响应构建目标函数，利用粒子群优化 算法进行极化参数提取，包括零频电导率、时间常数、充电率、频散系数、深 度五个参数；

进一步地，基于粒子群优化算法的柯尔-柯尔模型参数提取，如图2所示， 具体步骤如下：

Ⅰ、输入粒子群优化算法相关参数D、p_k、b_k、

x_limit、v_limit，其中D 为5，b_k为20，p_k为300，

均相等，大小均设定为2，

x_limit为[0.0464 0.0696；0.0044 0.0066；0 1；0 1；-50 -300]；v_limit为[-0.001,0.001；-0.001,0.001；-0.0001,0.0001；-0.002,0.002；-0.002,0.002]，依次为零频 电导率、时间常数、充电率、频散系数、深度的搜索范围；随机初始化粒子的 位置、速度和个体最优值P_best、种群最优值G_best；

Ⅱ、构建目标函数：minφ(X)＝φ₁(X)+λφ₂(X)，X＝σ₀、η、τ、c，其中

φ₂(X)＝f_t(X)，B_t为实测数据，f_t(X)为理论模型磁 场响应，λ为收敛因子；

Ⅲ、计算每个粒子的当前适应度，并与P_best、G_best进行比较，若当前适应 度值较小，则令当前适应值替代P_best、G_best；

Ⅳ、实测曲线依据斜率划分为两段，据此更新粒子的速度和位置：

x_ij(k+1)＝x_ij(k)+v_ij(k+1)·t (2)

式中，

t₀为斜率由负变为正的时刻，p_k为最大迭代次数，iter为当前迭代次数，下标i表示第i个粒 子(i＝1,2,3......M)，下标j表示粒子的第j维(j＝1,2,3......D)，假设粒子群搜索空 间为D维，粒子群有M个粒子(通常粒子数设置为问题维数的5～10倍)，r₁、 r₂是(0，1)之间的随机数。

为了保证粒子在搜索空间内搜索，通常将粒子的速度限制在[-v_max v_max]区间 中，将步骤3得到的搜索范围做为粒子的位置限制[-x_min x_max]，则粒子的最大 速度可以取x_max的10％～20％。

由于粒子的速度和位置均有方向性，迭代过程中，当其到达边界，即产生 变异时，大致可分为两种情况，若二者变异方向相反，则不做任何处理，反之， 将由下式对位置和速度重新定义，当粒子的位置或速度超出上边界时：

v_ij(k+1)＝v_max-2·r₃·v_max (3)

x_ij(k+1)＝x_max-0.5·r₄·x_max (4)

当粒子的位置或速度超出下边界时：

v_ij(k+1)＝v_min-2·r₃·v_min (5)

x_ij(k+1)＝x_min-0.5·r₄·x_min (6)

式中r₃、r₄是(0，1)之间的随机数。

Ⅴ、判断是否达到最大迭代次数，若是，获得P_best、G_best的位置信息，输出 σ₀＝0.051S/m，τ＝0.0050ms，η＝0.791，c＝1，d＝-107m；否则返回步骤Ⅲ，继续 执行；

5)进行极化介质的电导率、极化率-深度成像。

参见图1所示，利用提取到的极化参数，带入深度计算公式，绘制极化介 质的电导率、极化率-深度图像。

Ⅰ、应用matlab编程，输入步骤4提取出的零频电导率、时间常数、充电 率、频散系数到程序中，根据

获得复电导率值；

Ⅱ、将复电导率带入广义趋肤深度公式

中，可得

令i^c＝a+ib，

可得

其中ε为 介电常数，ω为角频率，σ_∞为高频交流电导率，μ₀是磁导率，计算该复电导率 σ(ω)的趋肤深度；

Ⅲ、利用步骤Ⅰ中的零频电导率、极化率和步骤Ⅱ获得的趋肤深度，分别 进行电导率-深度、极化率-深度成像。

图2是粒子群优化算法提取参数的流程图；

图3是根据地质资料建立的三层层状大地模型；

图4为采用图1所示的本发明一个实施例的层状模型成像对比图，其中a 是提取结果的零频电导率-深度成像图，b是提取结果的极化率-深度成像图。结 果表明，本方法的电导-极化率深度的成像结果，与实际大地一致，相比传统电 阻率计算方法，更符合实际地下极化介质的传播扩散规律，提高了电导率探测 深度的解释精度，为电磁勘探方法野外高精度解释提供了新的思路和方法。

以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发 明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明 的保护范围之内。

Claims (2)

1.一种基于粒子群优化算法的电导-极化率多参数成像方法，其特征在于，该方法包括：

1)基于地下极化介质的柯尔-柯尔复电导率模型，获得梯形波激励下的时域感应-极化磁场响应；

2)采用超导量子传感器时域电磁探测系统进行野外磁场数据观测，并对实测数据进行叠加取样、滤波、以及一次场剔除处理；

3)利用步骤2)数据的早期数据计算零频电导率值、晚期数据e指数多项式拟合得到时间常数，将获得的零频电导率值和时间常数对粒子群搜索范围进行约束；

4)根据实测磁场数据与步骤1的时域感应-极化磁场响应构建目标函数，利用粒子群优化算法，根据步骤3约束的范围，提取极化参数，包括零频电导率σ₀、时间常数τ、极化率η、频散系数c、深度d；

5)利用步骤4提取到实测数据的极化参数，带入广义趋肤深度计算公式，进行极化介质的电导率、极化率-深度成像；

步骤4中具体包含以下步骤：

4Ⅰ、输入相关参数，随机初始化粒子的位置、速度和个体最优值P_best、种群最优值G_best；

4Ⅱ、构建目标函数：minφ(X)＝φ₁(X)+λφ₂(X)，X＝σ₀、η、τ、c，

其中

φ₂(X)＝f_t(X)，f_t(X)为步骤1的时域感应-极化磁场响应，B_t为实测磁场数据，λ为收敛因子；

4Ⅲ、计算每个粒子当前适应度，与P_best、G_best进行比较，若当前适应度值较小，则令其替代P_best、G_best，反之继续保留P_best、G_best值；

4Ⅳ、实测曲线依据斜率划分为两段，据此分别按照式(1)和式(2)更新粒子的速度和位置：

x_ij(k+1)＝x_ij(k)+v_ij(k+1)·t (2)

式中，

t₀为曲线斜率由负变为正的时刻，p_k为最大迭代次数，iter为当前迭代次数；

迭代过程中，当粒子的速度或位置超出搜索范围，产生变异时，根据二者移动方向分为两种情况，若二者变异方向相反，则不做任何处理，反之，将由下式对位置和速度重新定义，当超出上边界时：

v_ij(k+1)＝v_max-2·r₃·v_max (3)

x_ij(k+1)＝x_max-0.5·r₄·x_max (4)

当粒子的位置或速度超出下边界时：

v_ij(k+1)＝v_min-2·r₃·v_min (5)

x_ij(k+1)＝x_min-0.5·r₄·x_min (6)

4Ⅴ、判断是否达到最大迭代次数，若是，获得P_best、G_best的位置信息，输出σ₀、τ、η、c、d；否则返回步骤4Ⅲ，继续执行；

步骤5中包含以下步骤：

5Ⅰ、应用matlab编程，输入步骤4提取出的零频电导率σ₀、时间常数τ、极化率η、频散系数c到程序中，根据

公式，获得复电导率值；

5Ⅱ、根据广义趋肤深度公式，计算该复电导率σ(ω)的趋肤深度；

5Ⅲ、利用步骤5Ⅰ中的零频电导率、极化率和步骤5Ⅱ获得的趋肤深度，分别进行电导率-深度、极化率-深度成像。

2.按照权利要求1所述的成像方法，其特征在于，步骤3中，发射电流关断后，早期响应以感应场为主，利用早期数据计算零频电导率值σ₁；晚期数据以极化场为主，对其晚期数据进行e指数拟合，得到时间常数值τ₁；零频电导率值σ₁和时间常数值τ₁的±20％分别作为零频电导率和时间常数值的粒子群搜索范围上下限。

Images