模拟退火法与阻尼最小二乘法联合反演岩石激电谱参数

陈 恒， 严良俊，代小威，李鹏飞

(1.长江大学 湖北非常规油气协同创新中心，湖北 武汉 430100；2.中国石油塔里木油田分公司 勘探开发研究院，新疆 库尔勒 841000)



模拟退火法与阻尼最小二乘法联合反演岩石激电谱参数

陈 恒1， 严良俊1，代小威1，李鹏飞2

(1.长江大学 湖北非常规油气协同创新中心，湖北 武汉 430100；2.中国石油塔里木油田分公司 勘探开发研究院，新疆 库尔勒 841000)

Cole—Cole模型是一种定量描述激电响应谱特性的方法，即可以用4个参数来定量描述岩石激电响应特性,这4个参数是岩矿石激电响应特征的基本参数。由于不同岩矿石具有不同的激电谱特征，因而存在着利用这4个Cole—Cole参数来识别它们的可能性。阻尼最小二乘法是应用地球物理学中常用的一种最优化反演技术，广泛地用在谱激电反演之中，具有很快的收敛速度。然而阻尼最小二乘法对初始值很依赖，非线性模拟退火反演方法能够很好地消除这种影响。因此联合这两种方法反演岩石激电谱参数，经理论和实测数据计算表明，反演能较好地拟合实测曲线，取得了一定的效果。

模拟退火法；阻尼最小二乘法；激电谱；反演； Cole—Cole模型

1 引 言

岩矿石的复电阻率频谱，一般都满足Cole—Cole模型，可以用四个Cole—Cole参数描述。一般，岩矿石四个Cole—Cole参数是通过对岩矿石的复电阻率频谱测量数据进行最优化反演来确定的。当其频谱能很好拟合实测的岩矿石的复电阻率频谱时，拟合所用的Cole—Cole参数就是岩矿石的Cole—Cole参数。频谱激电资料的谱激电反演，不论是常规的还是求真参数的，一般都采用阻尼最小二乘法最优化技术[1]。

阻尼最小二乘法又称马奎特迭代法，是应用地球物理学中常用的一种最优化反演技术。它通过自适应调整阻尼因子来达到收敛特性，具有更高的迭代收敛速度，但在阻尼最小二乘法的计算过程中，初值是一个很重要的因素，若选择的初值接近真解时，收敛速度很快且能得到收敛，但如果初值远离真解时，优化结果往往过早地陷入局部最优解，从而得到的结果完全背离真解。要解决该问题，方法一是要通过大量的原始信息对其初值有一个大概的估计，但这在现实工作中是很难做到的，另外一种方法就是选择一种合适的方法与其相结合，消除阻尼最小二乘法对初始值的依赖，并且还应具有很快的收敛速度[2,4]。

非线性模拟退火反演方法很好地克服了初始值设定这一缺点。每修改一次模型，将会随机产生一个新的状态模型，如果新模型能够被系统接受，则进行‘降温’处理，产生下一个新模型。这种接受新模型的方式使其成为一种全局最优算法[2]。该方法能跳出局部最优的“陷阱”，避免了在反演迭代过程中陷入目标函数的局部最小值，从而得到全局最优解，并得到理论证明和实际应用验证[3]。模拟退火方法虽然可以达到全局最小，但计算量太大，尤其是解的个数太多和解空间范围太大时，该方法很难收敛[5]。

利用阻尼最小二乘法和模拟退火方法各自的优缺点，国内外利用这两种方法进行反演的例子很多。刘建军等利用阻尼最小二乘法与模拟退火法结合实现非线性模型参数估计，并将此算法运用于植物光合作用光响应新模型对水稻的光响应曲线的拟合，得到的拟合参数与DPS拟合值极为接近[3]。王伟男等讨论了对于观测数据的非线性拟合，目前广为采用的是基于蒙特卡罗法的模拟退火方法和阻尼最小二乘法[6]。王芳等基于模拟退火算法中的全局收敛思想，建立了一种解非线性反演问题的新策略[7]。

2 基本原理与方法

2.1 复电阻率的Cole—Cole表达式

大多数岩，矿石的激电谱可以用4个Cole—Cole参数来很好的定量描述，一般都满足Cole—Cole模型，可以根据岩矿石的这4个参数来评价激电异常。其复电阻率的Cole—Cole表达式为：

(1)

式中，ρ0是零频率时的电阻率，单位是Ω·m；m是极化率，无量纲；c为复电阻率随频率变化程度的参数，称为频率相关系数，无量纲；τ为时间常数[1]，单位为s。

2.2 阻尼最小二乘法基本原理

阻尼最小二乘法是在最小二乘法和最速下降法的基础上发展起来的，它结合了最小二乘法与最速下降法的优点，并在两种方法之间取某种折衷，力图以最大的步长，同时又靠近最速下降方向，以保证稳定收敛，并加快收敛速度。这种方法又称马奎特法[8,9]。

假设要反演的一组实测地球物理数据向量用D表示，设其有n个数据；设正演模型的已知模型参数向量为X(频点的频率)，它也有n个数；设待求参数向量为P(Cole—Cole参数)，有m个数[1]。

正演模型是P和X的函数，设为G(P,X)，显然应当有D=G(P,X)。

对于正演模型G(P,X)，按照台劳级数展开并略去所有二阶以上的项可以得到：

G(P,X) =G(P0,X)+

(2)

上式可以写成矩阵形式：ΔG=AΔP

(3)

其中，ΔG是数据差列向量；A是雅克比矩阵；ΔP是参数差列向量；最小二乘法的拟合优度用x2值衡量，根据上式即可得到目标函数：

(4)

上式是Δpj的多元函数，取得极小值的充分条件是ε(P0)对所有的Δpj的偏导数为零。即：

=0

(5)

进一步化简可以得到:

(ATA)ΔP=ATΔG

(6)

由式(6)可以得到：

ΔP=(ATA)-1ATΔG

(7)

加入阻尼因子k,其中I是单位矩阵，即得

ΔP=(ATA+kI)-1ATΔG

(8)

ΔP即是使目标函数ε往极小方向减小的最佳修正步长，对于绝大多数地球物理问题来说，G(P,X)不是P的线性函数，因此要经过若干次类似迭代，才能使目标函数满足某一精度要求，这里采用一边迭代一边调整阻尼因子[10]。

2.3 模拟退火法基本原理

近年来，模拟退火法已在地球物理资料的反演中被广泛地采用，并取得了一些令人瞩目的成果，它已成为目前非线性反演的重要方法之一。它的基本思想是，把物理系统的能量作为最优化问题的目标函数，其能量达到最小，系统处于最优状态。整个最优化过程中，温度是一个至关重要的参数[8]。据统计热力学，物体中每个分子的状态服从吉布斯概率分布，即为：

(9)

式中：E(ri)为第i个分子的能量函数；ri为第i个分子所处的状态；k为玻尔兹曼常数；T为温度；而ρ(ri)为第i个分子的概率密度。其中温度T只是一个控制参数，玻尔兹曼常数k起一个尺度因子的作用，实际应用取为1。显然，T较小时，对ΔE起着一种放大作用，使低温时不易接受模型修改。这相当于物体在低温时分子被束缚在平衡位置附近[8,9]。

2.4 模拟退火法与阻尼最小二乘法联合反演流程

图1 模拟退火法与阻尼最小二乘法联合反演流程Fig.1 Flow chart of the joint inversion of simulated annealing method and damped least square method

3 反演结果分析

3.1 理论模型计算

选择Cole—Cole模型参数ρ0=25 Ω·m，m=0.6，τ=1 s，c=0.5，频率范围为0.01Hz～1 kHz，用所给的参数带入Cole—Cole模型进行理论计算[11，12]，然后把计算出来的数据作为初始理论数据，用模拟退火法和阻尼最小二乘法联合反演进行复电阻率的振幅反演。本文先随机给一个初始值，用模拟退火法进行反演，然后把模拟退火法反演出来的参数当作阻尼最小二乘法初始值进行反演，反演出来的参数见表1，联合反演结果见图2。为了使本方法更具有说服力，笔者在计算出来的理论数据中加入±5%以内的随机相对噪声，把加入噪声的数据再来进行联合反演，反演出来的参数见表1，联合反演结果见图3。

表1 理论模型联合反演参数

图2 理论模型联合反演结果Fig.2 Joint inversion result of theoretical model

图3 加入噪声理论模型联合反演结果Fig.3 Joint inversion result of theoretical model after adding noise

3.2 阻尼最小二乘法单独反演结果分析

本次反演通过实验实测某地区岩石1号点、2号点和3号点的复电阻率数据，研究复电阻率随频率的变化关系。实测频率范围是0.01～10 kHz，考虑到电磁耦合效应，我们只研究0.01～1 kHz的实验数据[13,14]。本次反演用单Cole—Cole模型进行反演,并只对幅值进行反演。先编程用阻尼最小二乘法进行反演Cole—Cole模型参数。在对1号点反演过程中，取初值ρ0=200 Ω·m，m=0.2，τ=100 s,c=0.1无法收敛，然后通过不断变化取初值，最后当取初值ρ0=100 Ω·m，m=0.6,τ=0.5 s,c=0.25时反演得到了岩石频谱参数，但也不是岩石真频谱参数。同样2号点和3号点也是通过不断变化取初值才最终反演出Cole—Cole模型参数。1号点、2号点和3号点的反演参数见表2，反演结果见图4～图6。这说明阻尼最小二乘法反演对初值严重依赖，如果初值选取不当，反演可能无法收敛或者不能收敛到全局最小[12]。

表2 各个点实测曲线反演参数

图4 1号点的实测曲线联合反演结果 Fig.4 Joint inversion result of the measured curve of the 1st point

图5 2号点的实测曲线联合反演结果Fig.5 Joint inversion result of the measured curve of the 2nd point

图6 3号点的实测曲线联合反演结果Fig.6 Joint inversion result of the measured curve of the 3rd point

3.3 阻尼最小二乘法和模拟退火法联合反演分析

由于阻尼最小二乘法对初始值严重依赖，因此本文通过阻尼最小二乘法和模拟退火法联合反演的思想来消除这种影响。先编程用非线性模拟退火法反演得到Cole—Cole模型参数，并作为阻尼最小二乘法反演的初始参数，这样阻尼最小二乘法就会快速反演得到岩石频谱参数。1号点、2号点和3号点的联合反演参数见表2。

图4～图6分别给出了1号点、2号点和3号点实测曲线，阻尼最小二乘法单独反演结果，模拟退火法单独反演结果，以及阻尼最小二乘法和模拟退火法联合反演结果。将实测曲线与反演得到的电图阻率进行对比，结果如图4～图6所示。

4 结 论

本文充分利用模拟退火法不依赖初值以及能收敛到全局最小的优点，编程实现了阻尼最小二乘法和模拟退火法联合反演的程序，并反演得到了岩石频谱参数，取得了一定的效果。经过对理论和实测的数据计算表明，这种联合反演能较好地拟合实测曲线，降低了阻尼最小二乘法对初始值的依赖程度，并具有一定的抗噪能力。

[1]刘崧.谱激电法[M].武汉：中国地质大学出版社，1997.

[2]罗延钟,张桂青.频率域激电法原理[M].北京:地质出版社,1988.

[3]刘建军，陈明峰.阻尼最小二乘法与模拟退火法结合实现非线性模型参数估计[J].井冈山大学学报,2010,31(6):10-14.

[4]武菊，任鹏.基于LM和SA的混合优化算法[J].内江师范学院学报，2010，25 (8):32-34.

[5]严良俊，胡文宝.大地电磁测深资料的二次函数逼近非线性反演[J].地球物理学报，2004,47(5):935-940.

[6]王伟男，童茂松，陈国华.泥质砂岩的物理性质及其测井应用[M].北京：石油工业出版社,2004 .

[7]王芳，陈生昌.解非线性反演问题的新策略[J].石油物探,1999,38(3):76-85.

[8]王家映.地球物理反演理论[M].北京：高等教育出版社，2002.

[9]姚姚.地球物理反演基本理论与应用方法[M].武汉：中国地质大学出版社,2002.

[10]薛嘉庆.最优化原理与方法[M].北京：冶金工业出版社，1992.

[11]李鹏飞，严良俊，谢兴兵,等.基于最小二乘的岩石频谱参数反演初始值选取的讨论[J].地震工程学报,2014，36(3):569-574.

[12]向葵，胡文宝，严良俊,等.川黔地区页岩复电阻率的频散特性[J].石油地球物理勘探，2014，49(5):1 013-1 019.

[13]周新鹏，邹长春，乔祎,等.岩矿石电性参数测定实验研究[J].工程地球物理学报，2015,12(2):205-213.

[14]李鹏飞，严良俊，余刚.南方页岩岩芯的复电阻率测试与分析[J].工程地球物理学报，2014,11(3):383-386.

The Joint Inversion of Simulated Annealing Method and Damped Least Square Method for IP Spectrum Parameters of Rock

Chen Heng1,Yan Liangjun1,Dai Xiaowei1,Li Pengfei2

(1.HubeiCooperativeInnovationCenterofUnconventionalOilandGas,YangtzeuUniversity,WuhanHubei430100,China；2.ResearchInstituteofExplorationandDevelopment,TarimOilfieldCompany,PetroChina,KorlaXinjiang841000,China)

The Cole—Cole model is a quantitative method for describing the characteristics of the excitation response spectrum, that is, the spectral characteristics of the IP response can be quantitatively described with 4 parameters, and these 4 parameters are the basic parameters of the characteristics of rock IP response. Because different rocks have different IP spectral characteristics, it is possible to identify them by using the four Cole—Cole parameters. Damped least square method is one of the commonly used optimization inversion techniques in the application of geophysics, which is widely used in the spectral inversion and has a fast convergence rate. However, the damped least squares method is very dependent on the initial value, and the nonlinear simulated annealing inversion method can eliminate the effect of this method. So the two methods were combined to inverse the parameters of the rock. The calculation of the theory and the measured data showed that the joint inversion can better fit the measured curves, and good results had been obtained.

simulated annealing method; damped least squares method; IP spectrum; inversion; Cole—Cole model

1672—7940(2016)02—0170—05

10.3969/j.issn.1672-7940.2016.02.006

国家自然科学基金项目(编号：U1562109,41274082)；国家“973”项目(编号：2013CB228605)

陈 恒(1990-)，男，硕士研究生，主要研究方向为电磁法勘探。E-mail：1033256984@qq.com

P631.3

A

2015-11-10