基于两步变异差分进化算法的激电测深一维反演

doi:10.11720/wtyht.2023.0058

基于两步变异差分进化算法的激电测深一维反演

丁志军^,¹, 罗维斌^,², 连伟章¹, 张星¹, 何海颦²

1.甘肃省有色地质调查院,甘肃兰州 730000

2.兰州资源环境职业技术大学地质与珠宝学院,甘肃兰州 730021

One-dimensional inversion of induced polarization sounding data based on the differential evolution algorithm with two-step mutation

DING Zhi-Jun^,¹, LUO Wei-Bin^,², LIAN Wei-Zhang¹, ZHANG Xing¹, HE Hai-Pin²

1. Gansu Nonferrous Geological Survey Institute, Lanzhou 730000, China

2. Lanzhou Resources and Environment Vocational and Technical University, College of Geology and Jewelry, Lanzhou 730021, China

通讯作者: 罗维斌(1972-),男,地球探测与信息技术工学博士学位,正高级工程师,主要从事地球物理电磁法勘探教学和应用研究工作。Email:lwb210521@lzre.edu.cn

第一作者: 丁志军(1987-),男,地球物理学学士学位,工程师,主要从事电磁法勘探应用研究工作。Email:609415517@qq.com

责任编辑: 王萌

收稿日期: 2023-02-15 修回日期: 2023-05-30

基金资助:

甘肃省教育厅产业支撑计划项目(2021CYZC-67)
兰州资源环境职业技术大学科技项目(Y2021B-01)

Received: 2023-02-15 Revised: 2023-05-30

摘要

激电测深一维反演是一个多参数非线性优化问题。本文基于改进的两步变异差分进化全局最优化算法实现了激电测深的一维反演。传统的差分进化算法包含变异、交叉和选择操作,变异为单步变异。本文提出的两步变异法分步将最优个体与随机选取的两个个体经变异后产生新个体。加强了最优个体的影响度,提高了全局寻优能力。通过模型试算结果表明,两步变异法比传统方法寻优能力更强。利用等效电阻率法加载极化率参数,通过数字滤波算法可快速正演计算层状模型表面激电测深电阻率曲线,在此基础上应用两步变异差分进化算法不断变异产生新个体,正演计算电阻率与观测值进行拟合,选择适应度值趋近于最大适应度值的个体作为反演结果。本文反演方法操作简便,计算速度快。通过对H型和KH型地电模型进行计算,得出本反演方法有较高的拟合精度。

关键词： 激电测深; 非线性一维反演; 两步变异差分进化优化算法; 电阻率; 极化率

Abstract

The one-dimensional inversion of induced polarization (IP) sounding data involves multi-parameter nonlinear optimization. This study achieved the one-dimensional (1D) inversion of IP sounding data based on the improved global optimization algorithm of differential evolution (DE) with two-step mutation. The conventional DE algorithm includes mutation (single-step), crossover, and selection operations. The two-step mutation method proposed in this study can produce new individuals through the mutation of the optimal individual and two randomly selected individuals in steps, thus enhancing the influence of the optimal individual and the global optimization ability. The model test results show that the two-step mutation method has a higher optimization ability than the conventional method. Specifically, the polarizability parameters were loaded using the equivalent resistivity method, and the surface IP sounding resistivity curves of a layered model can be quickly calculated through forward modeling using the digital filtering algorithm. Based on this, the DE algorithm with two-step mutation was employed to produce new individuals through continuous mutation. Then, the resistivity obtained through forward modeling was fitted with the observed values, and the individuals whose fitness approached the maximum fitness were selected as the inversion results. The inversion method proposed in this study features simple operations and fast calculations. As verified through the calculations of H- and KH-type geoelectric models, the inversion method enjoys high fitting accuracy.

Keywords： IP sounding; nonlinear one-dimensional inversion; optimization algorithm of differential evolution with two-step mutation; resistivity; polarizability

PDF (1299KB) 元数据多维度评价相关文章导出 EndNote| Ris| Bibtex 收藏本文

本文引用格式

丁志军, 罗维斌, 连伟章, 张星, 何海颦. 基于两步变异差分进化算法的激电测深一维反演[J]. 物探与化探, 2023, 47(4): 1033-1039 doi:10.11720/wtyht.2023.0058

DING Zhi-Jun, LUO Wei-Bin, LIAN Wei-Zhang, ZHANG Xing, HE Hai-Pin. One-dimensional inversion of induced polarization sounding data based on the differential evolution algorithm with two-step mutation[J]. Geophysical and Geochemical Exploration, 2023, 47(4): 1033-1039 doi:10.11720/wtyht.2023.0058

0 引言

激发极化法是一种应用广泛的传导类电法,利用岩矿石在人工电流场作用下产生激发极化效应的差异性进行地质勘查,可获得视电阻率和视极化率等多个参数来表征岩矿石物理属性,可在时间域进行观测,也可在频率域进行观测^[1-2]。由于激电异常对地下水、硫化物矿体等地质体分布指示性强,因此在浅地表800 m深度范围内的各类地质勘查中激电法发挥了重要作用^[3]。因观测方式不同,时间域激发极化法利用关断电流后测量的二次电位与关断电流前的总电位的比值,即极化率,来表征岩矿石的激发极化效应;在频率域则形成了基于双频或多频响应的相对相位或百分比频率效应参数表征岩矿石激发极化效应的差异^[4]。

在近一个世纪的发展过程中,激发极化法在方法原理、观测设备、数据处理及正反演解释等方面均有长足的进步^[4]。在时间域需要测量断电后的二次电场,信号微弱,且易受人文或工业噪声干扰,因此频率域激发极化法表现出大深度探测能力。从变频法到双频法,再到多频激电法,频率域激电法在深地勘探中发挥着重要作用。中南大学何继善院士基于2ⁿ序列伪随机信号理论,创新发展了双频激电法和伪随机多频激电法,有力推动了频率激电法的应用发展。陈儒军等^[5]将扩频信号(PRBS)引入频谱激电法中,研制了精密扩频激电仪,找矿案例显示出大深度的探测能力,且有较强的抗干扰能力。

由于传导类电法对地探测的体积效应明显,直接利用观测结果进行地质解释分辨率不高,精细化反演是提高解释水平的必由之路^[6]。激发极化效应是地质体在人工电流场作用下产生的电化学反应,是典型的非线性问题。目前常规反演方法均采用将非线性问题线性化后基于最小二乘法或共轭梯度法拟合获得反演结果。近些年来随着计算力的快速提升,一些智能优化算法也应用到了电阻率和极化率的联合反演当中,包括模拟退火算法、遗传算法、神经网络算法和群智能算法等^[7⇓-9],这些算法的核心思想都是进行全局寻优,并且已经在一维和二维反演中得到了较好的应用。

差分进化(differential evolution,DE)最优化算法是一种新兴的随机启发式进化计算技术,因其特有的记忆能力,且不需要借助问题的特征信息,简单易用、稳健性好且全局搜索能力强大,已在多个领域得到应用^[10]。高颖等^[11]提出了一种差分进化协同定位算法,通过差分进化算法对物联网定位问题进行优化,改善了定位性能。闫秀英等^[12]提出了一种具有差分进化算子的粒子群优化算法(改进 DEPSO 算法)对中央空调冷水机组负荷分配进行优化,降低了多冷水机组系统在不同典型负荷率下的能耗。丁学振等^[13]提出了基于改进差分进化算法的误差参数估计方法显著降低了磁梯度张量系统的测量误差,使测量精度得到明显改善。王天意等^[9]基于改进的差分进化算法实现了大地电磁的一维及二维完全非线性反演,有效提高了大地电磁反演的精度。廖作文^[14]研究了基于差分进化算法在一次运行中实现非线性方程组的多根联解。

本文通过两步变异法对差分进化算法进行了改进,较原型算法优化性能有显著改善。在层状模型直流电阻率测深一维反演基础上^[15],通过等效电阻率法引入极化率参数,运用改进的差分进化算法优化这一多目标函数,实现了激电测深非线性一维反演。

1 激电测深正演计算方法

正演是反演的基础,在一维层状模型垂向直流电阻率测深数字滤波算法基础上,采用等效电阻率法加载极化率参数实现激电测深一维正演。地表观测的电阻率值是各层层参数的非线性综合响应:

(1)

ρ_{s} = f (ρ_{1}, d_{1}, η_{1}, ρ_{2}, d_{2}, η_{2}, \dots, ρ_{N}, η_{N}, A B / 2)

式中: $ρ_{1}, ρ_{2}, \dots, ρ_{N}$ 为各层电阻率; $d_{1}, d_{2}, \dots, d_{N - 1}$ 为层厚度,最底层为厚度无限大的均匀半空间; $η_{1}, η_{2}, \dots, η_{N}$ 为各层极化率参数; $A B / 2$ 为观测电极距。

根据地表四电极地电阻率测量视电阻率的微分形式,复杂地质体视电阻率是MN测量电极实测的电流密度 $J_{M N}$ 或电场强度 $E_{M N}$ 与水平地形均匀半空间模型同一装置的电流密度 $J_{0}$ ,或电场强度 $E_{0}$ 的比值与 $ρ_{1}$ 的乘积:

(2)

ρ_{s} = ρ_{1} \frac{J_{M N}}{J_{0}} = ρ_{1} \frac{E_{M N}}{E_{0}}

对于层状地电模型,有:

(3)

E_{M N} = \frac{I ρ_{1}}{2 π} \int_{0}^{\infty} R_{N}^{0} (λ) J_{1} (λ r) λ d λ

而均匀半空间地电模型, $E_{0} = \frac{I ρ_{1}}{2 π r^{2}}$ ,故,视电阻率,由式(2),有:

(4)

ρ_{s} = ρ_{1} r^{2} \int_{0}^{\infty} R_{N}^{0} (λ) J_{1} (λ r) λ d λ

在层参数已知情况下,式(4)可由Hankel数值积分数字滤波算法快速计算。

式中电极距 $r = A B / 2$ ,对于对称四极测深装置,极距影响不容忽略^[15],可用 $r = [{(\frac{A B}{2})}^{2} - {(\frac{M N}{2})}^{2}] / [4 (\frac{M N}{2})]$ 替换; $J_{1}$ 为第一类Bessel函数; $R_{N}^{0}$ 为地表的地层校正系数,第N层均匀半空间顶面的地层校正系数 $R_{N}^{N} = 1$ ,由底层到地表面的地层校正系数由下式递推可得:

R_{N}^{i} = \frac{1 - K e^{- 2 d_{i} λ}}{1 + K e^{- 2 d_{i} λ}}

(5)

K = \frac{1 - (ρ_{i + 1} / ρ_{i}) R_{N}^{i}}{1 + (ρ_{i + 1} / ρ_{i}) R_{N}^{i}}; i = 1,2, \dots, N - 1

。

式(5)中的层电阻率是加载了极化率的等效电阻率:

(6)

ρ_{i}^{*} = ρ_{i} / (1 - η_{i})

。

式(4)的计算采用了数字滤波算法^[16-17]。

2 改进的两步变异差分进化算法

差分进化(differential evolution,DE)算法是模拟自然界生物种群以“优胜劣汰,适者生存”的进化发展规律而形成的一种随机启发式搜索算法,是一种新兴的进化计算技术。它于1995年由Rainer Storn和Kenneth Price提出。方法简单易用、稳健型好,有强大的全局搜索能力^[18-19]。

差分进化算法保留了基于种群的全局搜索策略,采用实数编码、基于差分的简单变异操作和一对一的竞争生存策略,降低了遗传操作的复杂性。同时,差分进化算法特有的记忆能力使其可以动态跟踪当前的搜索情况,以调整其搜索策略,具有较强的全局收敛能力和鲁棒性,且不需要借助问题的特征信息,适于求解一些利用常规的数学规划方法所无法求解的复杂环境中的优化问题。采用差分进化算法可实现复杂系统的参数辨识^[18]。

实验结果表明,差分进化算法的性能优于粒子群算法和其他进化算法,该算法已成为一种求解非线性、不可微、多极值和高维复杂函数的有效方法。是目前最强大的随机实参数优化算法之一。差分进化算法中差分变异算子是最重要的算子,在种群中随机选择两个个体向量进行差分与缩放,同时与第三个个体向量相加得到一个变异的子个体,然后与父个体向量进行杂交得到新个体向量,最后将其与父个体比较适应值的大小,选择较优的用于子代的进化中。差分进化算法是根据个体之间的差异来实现个体进化的。开始先根据搜索空间的均匀分布,随机生成初始种群。获得初始种群之后,该算法通过变异、交叉和选择操作来更新个体的种群。这3个操作构成了差分进化算法的主循环,直到算法结束(图1)。

图1

新窗口打开| 下载原图ZIP| 生成PPT

图1 差分进化流程

Fig.1 Flow chart of differential evolution

2.1 种群初始化

设置种群规模为NP,模型参量的空间维数为D,及每个参量的上、下界范围,在模型的可行解空间内随机产生实数编码的初始种群pop为

(7)

\begin{array}{l} p o p (N P, j) = X V_{j m i n} + r a n d (N P, 1) \times \\ (X V_{j m a x} - X V_{j m i n}), j = 1,2, \dots, D \end{array}

式中: $X V_{j m i n}$ 为第j个模型参量的取值下界; $X V_{j m a x}$ 为第j个模型参量的取值上界;D为模型参量维数;rand为产生的(0,1)间的随机数。

2.2 差分变异操作

现有差化进化算法有以下几种差分变异方法:

1)DE/best/1: $V_{i}^{} = b m + F \times (p m_{1} - p m_{2})$ ;

2)DE/rand/1: $V_{i}^{} = p m_{1} + F \times (p m_{2} - p m_{3})$ ;

3)DE/rand-to-best/1: $V_{i}^{} = X V_{i} + F \times (b m - X V_{i}) + F \times (p m_{1} - p m_{2})$ ;

4)DE/best/2: $V_{i}^{} = b m + F \times (p m_{1} - p m_{2} + p m_{3} - p m_{4})$ ;

5)DE/rand/2: $V_{i}^{} = p m_{1} + F \times (p m_{2} - p m_{3} + p m_{4} - p m_{5})$ ;

其中:bm为种群中的最优个体;XV_i为父种群旧个体;F为变异系数;pm_n为从种群NP中随机选取的个体;V_i为变异产生的新个体。产生的变异新个体还需要检查各参量是否在参数范围内。若小于参量范围,以参量范围下界值代替,若大于参量范围,以参量范围上界值代替。

2.3 交叉操作

交叉操作的目的是提高种群个体的多样性,由交叉系数CR控制,方式如下:

(8)

V m_{i, j} = \{\begin{array}{l} V_{i, j}, r a n d \leq C R, j = 1,2, \dots, D; \\ X V_{i, j}, r a n d > C R, j = 1,2, \dots, D 。 \end{array}

式中:V为变异产生的新个体;XV为种群旧个体;根据产生随机数,分别从V和XV个体中选取第j个参数,组成新的个体Vm;rand为产生(0,1)范围内的随机数;CR为交叉控制系数。

2.4 选择操作

计算当前操作个体XV_i和新变异、交叉产生新个体Vm_i的个体适应度值,保留适应度值高的个体进入下一轮变异:

(9)

X V_{i} = \{\begin{array}{l} V m_{i}, f (V m_{i}) > f (X V_{i}); \\ X V_{i}, f (V m_{i}) \leq f (X V_{i}) 。 \end{array}

其中:f()表示适应度计算。

2.5 改进的差分进化方法

本文对差分进化算法做两步变异法改进,每步均对种群所有个体遍历完成变异、交叉、选择操作,流程如图2所示:

图2

新窗口打开| 下载原图ZIP| 生成PPT

图2 两步变异法流程

Fig.2 Flow chart of two-step variation method

第一步变异法: $V_{i}^{} = b m + F \times (p m_{1} - p m_{2})$ ;

第二步变异法: $V_{i}^{} = p m_{1} + F \times (b m - p m_{2})$ ;

其中:bm为最大适应度值个体;pm_n为在种群NP中随机选取的2个个体,F取值1.21×rand()。

2.6 性能测试

模型参数如表1所示。先正演计算出180个电极距(极距范围1~4 000 m)的地表视电阻率,再用传统的单步变异差分算法和本文提出的两步变异法反演进行比较。单步变异循环一次按一次统计,两步变异法循环一次按两次统计,在变异次数相同的情况下进行比较。从图3可以看出,两步变异法(Modified DE)较传统单步变异法(DE)在550次变异中有明显的性能提升。计算参数设置:种群规模NP=45,模型参数维度D=8,模型参数取值及反演取值范围见表1第2、3列。个体的适应度值MF通过设置最大适应度值FML=1 500,和计算的拟合均方差Nerr通过式(8)转换为趋近于1的值,当Nerr趋近于0,MF趋近于1,当MF大于FML或达到设定最大迭代次数则终止计算。经550次迭代后,传统单步变异法的拟合均方差Nerr=4.39,MF=0.994 164,两步迭代法的Nerr=0.061 7,MF=0.999 996 263。算法改进效果明显。

表1 H型地电模型参数及反演结果

Table 1 Parameters and inversion results of H-type earth resistivity model

模型参数	真值	取值范围	反演次数					均值	误差/%
模型参数	真值	取值范围	1	2	3	4	5	均值	误差/%
ρ₁/(Ω·m)	100	[10,300]	99.72	99.50	99.78	99.50	100.32	99.76	-0.24
ρ₂/(Ω·m)	30	[10,150]	31.15	29.93	31.80	30.09	30.31	30.66	2.19
ρ₃/(Ω·m)	110	[10,300]	110.56	110.56	107.63	107.27	107.34	108.67	-1.21
d₁/m	180	[10,300]	180.00	180.00	180.00	180.00	180.00	180.00	0.00
d₂/m	60	[10,300]	59.99	60.01	59.99	60.01	60.00	60.00	0.00
η₁/m	0.5	[0.1, 1]	0.78	1.00	0.72	1.00	0.18	0.74	0.24
η₂/m	8.5	[3, 11]	4.97	8.72	3.00	8.24	7.56	6.50	-2.00
η₃/m	1	[0.5, 3.5]	0.50	0.50	3.13	3.46	3.40	2.20	1.20

新窗口打开| 下载CSV

图3

新窗口打开| 下载原图ZIP| 生成PPT

图3 传统单步变异差分进化与两步变异差分进化性能比较

Fig.3 Performance comparison between traditional differential evolution and two-step mutation differential evolution

(10)

M F = (F M L - N e r r) / F M L

(11)

N e r r = \sqrt{\overset{n}{\sum_{i = 1}} {(o b d a t - m d a t)}^{2}}

式中:obdat是地表电阻率观测值,算例中是按第2列模型参数正演计算值;mdat是种群个体的模型正演计算电阻率值;n为AB/2的电极距个数。以拟合均方差最小的个体参数作为反演结果。

3 基于改进差分进化的激电测深一维反演

3.1 H型模型

模型参数如表1所示。运用改进的两步变异差分进化进行反演计算。按模型层参数,给每个参数设置取值范围(表1第3列)。设置FML=9,MF=0.999 999 85,种群规模NP=45,最大迭代次数为210。

由于种群初始化的随机性,因此多次反演结果会有差异。给定参数不变,计算10次,取拟合误差最小的5次反演结果进行统计,得到反演结果的均值(表1第9列)和误差(表1第10列),其中极化率值为百分数。层电阻率和层深度采用均方误差统计,层极化率参数采用绝对误差统计。从表1反演结果可以看出,反演结果与模型真值很接近。第二层电阻率均方误差最大,为2.19,层厚的均方误差最小,几乎为0。反演结果的变化也是一维层状模型等值现象的响应。图4是模型参数电阻率正演值与反演结果模型的电阻率正演值,视电阻率曲线有较高的拟合度(图中两个曲线几乎重合)。最小拟合均方差为1.08×10^-4。

图4

新窗口打开| 下载原图ZIP| 生成PPT

图4 H型模型真实模型和反演结果的电阻率测深结果对比

Fig.4 Comparison of resistivity sounding results between real model and inversion results of H-type model

3.2 KH型模型

四层模型参数如表2所示。运用改进的两步变异差分进化进行反演计算。按模型层参数,给每个参数设置取值范围(表2第3列)。设置FML=6,MF=0.999 999 85,模型参数维度D为11,种群规模NP=60,最大迭代次数为500。

表2 KH型地电模型参数及反演结果

Table 2 Parameters and inversion results of KH-type earth resistivity model

模型参数	真值	取值范围	反演次数					均值	误差/%
模型参数	真值	取值范围	1	2	3	4	5	均值	误差/%
ρ₁/(Ω·m)	120	[10,300]	119.52	120.24	120.39	120.48	120.47	120.22	0.18
ρ₂/(Ω·m)	560	[100,800]	557.21	557.26	559.88	559.71	557.59	558.33	-0.30
ρ₃/(Ω·m)	25	[10,100]	26.04	26.85	26.10	24.37	24.32	25.53	2.14
ρ₄/(Ω·m)	310	[50,500]	310.58	306.77	303.68	304.13	310.54	307.14	-0.92
d₁/m	110	[50,300]	110.00	110.00	110.04	110.00	110.00	110.01	0.01
d₂/m	120	[50,300]	120.00	119.98	119.20	119.98	120.02	119.83	-0.14
d₃/m	150	[50,300]	150.03	150.29	160.59	150.31	149.68	152.18	1.45
η₁/m	0.6	[0.2,1]	1.00	0.40	0.28	0.20	0.21	0.42	-0.18
η₂/m	0.3	[0.1,0.8]	0.80	0.80	0.65	0.36	0.72	0.67	0.37
η₃/m	9.5	[3,12]	5.76	3.00	11.59	11.95	11.80	8.82	-0.68
η₄/m	1	[0.8,3.5]	0.81	2.03	3.08	2.88	0.82	1.93	0.93

新窗口打开| 下载CSV

同样,多次反演结果因种群初始化的随机性而有差异。给定参数不变,计算10次,取拟合误差最小的5次反演结果进行统计,得到反演结果的均值(表2第9列)和误差(表2第10列),表中极化率值为百分数,层电阻率和层深度采用均方误差统计,层极化率参数采用绝对误差统计。从表2反演结果可以看出,反演结果与模型真值很接近。第三层电阻率均方误差最大,为2.14,第一层厚d₁的均方误差最小,均方误差为0.01。图5是KH模型参数电阻率正演值与反演结果模型的电阻率正演值,视电阻率曲线同样有较高的拟合度(2个电阻率测深曲线重合在一起)。最小拟合均方差Nerr为2.73×10^-4。

图5

新窗口打开| 下载原图ZIP| 生成PPT

图5 KH型模型真实模型和反演结果的电阻率测深结果对比

Fig.5 Comparison of resistivity sounding results between real model and inversion results of KH-type model

4 结果讨论

两步变异差分进化算法中均利用了最优个体参与变异,增加了最优个体的影响力。通过模型试算,两步变异法能使适应度值快速增加,趋近于最优解。

差分进化算法适合于多参数目标函数的反演。尽管层状模型地表电阻率可以通过数值计算快速实现,但是考虑极化率参数的反演并不容易实现。通过激电测深极化率数据辅助确定极化率层参数变化范围,等效电阻率法拟合电阻率数据的反演方式,基于两步变异差分进化全局优化算法实现了激电测深的一维反演,验证了方法的可行性。

差分进化算法的种群规模NP对反演效果有很大影响,需要合理选择。一般取参数维度D的5~10倍为宜,本文试算中NP取值45~70。

5 结论

非线性最优化算法是地球物理勘探领域一种重要的反演方法。差分进化全局优化算法应用于反演问题表现出一定的优良性能。本文在常规差分进化单步变异算法基础上,用两步变异法加以改进,在每一步变异中都有最优个体,增加了最优个体的影响度。从验证计算结果看,两步变异法较常规一步变异法有明显的寻优能力。

基于层状模型直流电测深电阻率数字滤波算法,采用等效电阻率法加载极化率层参数,应用改进的两步变异差分进化算法实现了层状模型电阻率、层厚度和极化率的多参数反演。反演仅依赖观测数据和正演计算方法。初始种群根据种群规模、参数维度及其取值范围随机产生,建议多次试算确定种群规模和参数取值范围。

本文计算是以含极化率信息的正演电阻率结果作为观测值实现反演的,验证了方法的有效性。下一步实现含噪声实测数据的反演,改进算法的健壮性。

参考文献

原文顺序

文献年度倒序

文中引用次数倒序

被引期刊影响因子

[1]

李金铭

. 地电场与电法勘探[M]. 北京: 地质出版社, 2005.