基于模拟退火算法的大地电磁和地震数据贝叶斯联合反演

陈　晓，　周　磊,　于　鹏，　邓居智，　易寒婷

(1.东华理工大学放射性地质与勘探技术国防重点学科实验室，江西 南昌　330013；2.中国地质大学(武汉)地球内部多尺度成像湖北省重点实验室，湖北 武汉　430074；3.长江大学油气资源与勘探技术教育部重点实验室，湖北 武汉　430100；4.同济大学海洋地质国家重点实验室，上海　200092；5胜利油田测井公司，山东 东营　257000)



基于模拟退火算法的大地电磁和地震数据贝叶斯联合反演

陈晓1，2，周磊3,于鹏4，邓居智1，易寒婷5

(1.东华理工大学放射性地质与勘探技术国防重点学科实验室，江西 南昌330013；2.中国地质大学(武汉)地球内部多尺度成像湖北省重点实验室，湖北 武汉430074；3.长江大学油气资源与勘探技术教育部重点实验室，湖北 武汉430100；4.同济大学海洋地质国家重点实验室，上海200092；5胜利油田测井公司，山东 东营257000)

摘要：基于模拟退火算法实现了大地电磁和地震数据贝叶斯同步联合反演。贝叶斯反演是反演领域的研究热点。目前，虽已有其他方法之间的贝叶斯联合反演文献发表，但鲜见大地电磁和地震数据贝叶斯联合反演文献。充分考虑数据噪音不确定性的影响，将联合反演的解表达为后验概率密度分布，并提取了电阻率和速度的最大后验概率解、均值解和最终迭代解。

关键词：大地电磁；地震；联合反演；贝叶斯反演；模拟退火；噪音

陈晓，周磊,于鹏，等.2016. 基于模拟退火算法的大地电磁和地震数据贝叶斯联合反演[J].东华理工大学学报：自然科学版，39(1)：59-66.

Chen Xiao, Zhou Lei,Yu Peng, et al.2016. Bayesian joint inversion of magnetotelluric and seismic data based on simulated annealing method[J].Journal of East China University of Technology (Natural Science), 39(1):59-66.

Tarantola(1987)是地球物理贝叶斯反演的先驱，将数理统计学的思想引入地球物理反演，此后贝叶斯反演理论被广泛应用于各种地球物理反问题的求解中。Kane等(2004)将贝叶斯理论应用于三维慢度反演，在空间上按照深度拾取了慢度场的不确定度的分布，进而提出了预测井位。Buland等(2006)将贝叶斯理论应用于时移反演，将油气藏弹性参数的解表达为后验概率密度分布(PPD)，所得结果与实际相吻合。Hoversten等(2006)在贝叶斯框架下实现了海洋电磁法和地震数据联合反演，并在有井资料的地区将联合反演结果与单独反演进行了对比，认为联合反演确定的含气饱和度、含油饱和度以及孔隙度更符合井资料。Chen等(2009)实现了海洋电磁法和地震数据贝叶斯反演，并应用于含气饱和度的估计中，指出贝叶斯反演可以定量的评价参数的不确定度，此外尽管地震方法有着较高的分辨率，但是由于地震方法本身对深层气含量不敏感，因此很难区分深层气含量的状况，电磁方法的加入可以降低含气饱和度估计的不确定度，增强对深层含气饱和度的分辨能力。王文涛(2009)将贝叶斯反演应用于地震储层评价，利用模拟退火算法进行采样，将解表述为最优解和均值解，并指出均值解结果优于最优解。赵峦啸(2010)在贝叶斯框架下综合了地震数据、地震信息、测井信息和岩芯信息来对储层进行定量的评价。杨迪坤等(2011)在大地电磁反演中，将数据空间的不确定度转移到模型空间，使用概率分布反映了模型参数取值的可能性。郭荣文(2011)采用改进的自适应纯形模拟退火算法(ASSA)实现了一维大地电磁贝叶斯反演，使用Metroplois-Hastings采样，获取地电介质的不确定度信息。刘文劼(2012)将噪声相关性引入一维大地电磁反演中，并构建非对角协方差矩阵。

通过以上分析，可以看出贝叶斯反演是目前反演领域研究热点，目前已有基于贝叶斯概率思想的联合反演文献发表，但在大地电磁和地震数据联合反演领域中却鲜见文献涉及。作者曾对大地电磁和地震数据联合反演的研究现状做过详细分析(陈晓等，2010，2011；陈晓，2013)，不再赘述。

1基于模拟退火的贝叶斯联合反演

1.1传统的正规化联合反演目标函数

依据正则化反演理论(Zhdanov等，2002)，大地电磁和地震数据同步联合反演的目标函数包括了两大部分，分别为数据拟合泛函和模型约束泛函，具体为：

P(m,d)α=φ(m)+αS(m)=WMT(WdAMT(m)-WddMT)T(WdAMT(m)-WddMT)+Ws(WdAS(m)-WddS)T(WdAs(m)-Wdds)+α(WeWmm-WeWmmapr)T(WeWmm-WeWmmapr)

(1)

其中，d是观测数据，m是模型参数，Pα(m,d)为参数泛函，φ(m)为数据拟合泛函, S(m)为模型稳定泛函,α为正则化因子，WMT、WS分别是大地电磁和地震方法的权系数，AMT、AS分别是大地电磁和地震方法正演算子，dMT为视电阻率log(ρa)或相位，ds为地震旅行时，m为电阻率和速度，mapr为大地电磁和地震的先验信息，Wd、Wm为数据、模型加权矩阵,We为稳定泛函矩阵，具体形式参见文献(陈晓，2013)。

1.2贝叶斯基本理论

在地球物理反演领域，Tarantola(1987)提出了结合数据信息和模型信息的模型后验概率的公式：

(2)

其中，P(d|m)是模型m下的条件概率，也称似然函数；P(m)是模型m的先验概率，它包含了多种先验信息(地质信息、钻井信息等)；P(d)是观测数据的全空间概率，它与模型参数无关，在仅关心后验概率密度分布形态的条件下，P(d)可被忽略(陈建江等，2006)；σ(m|d)是在观测数据d下模型m的后验概率，可表述为：

σ(m|d)∞P(d|m)P(m)

(3)

假定观测数据的噪音和计算误差为高斯分布时，似然函数P(d|m)等价于函数exp[-Sφ(m)],S为一个比例因子，φ(m)为数据拟合泛函。

在获得后验概率分布之后，便可以根据概率理论，求取均值模型(王文涛，2009；刘文劼，2012)，公式如下：

I(m)=∫mσ(m|d)dm

(4)

1.3非常快速的模拟退火算法

本文的大地电磁和地震数据贝叶斯联合反演，是结合全局寻优的非常快速的模拟退火算法(杨辉，2000；于鹏等，2009)实现的。简而言之，模拟退火算法主要有以下三部分组成：模型扰动，接受准则，退火计划。

(1) 模型扰动

采用Ingber(1989)提出的与温度相关的似Cauchy分布随机产生新模型：

mi′=mi+yi[Bi-Ai]

(5)

yi=Tsgn(u-0.5)[(1+1/T)|2u-1|-1]

(6)

其中，mi为当前模型的第i个参数，u为(0，1)区间内的随机数，[Ai,Bi]为mi的扰动区间，mi′为新模型中第i个参数，sgn为符号参数,T是当前温度。

(2) 退火计划

采用Ingber(1989)提出的退火计划来缓慢降低温度。

T(C)=T0Dk

(7)

其中，T0为初始温度，k为迭代次数，D为衰减因子。

(3) 接受准则

采用张霖斌等(1997)提出的基于广义Boltzmann-Gibbs分布的接受准则：

P=[1-(1-q)ΔE/T]1/(1-q)

(8)

其中，ΔE=E(m＇)-E(m)，E(m＇)和E(m)分别为m＇,m状态下的能量，q为一实数。

1.4基于模拟退火算法的贝叶斯联合反演流程

结合贝叶斯反演的基本理论和模拟退火算法的基本原理，作者梳理了基于模拟退火算法的贝叶斯联合反演的具体流程：

(1)提取先验概率分布。整理收集各种地质和地球物理的先验信息，并将这些信息转换为先验概率分布P(m)。

(2)积累样本。进入模拟退火算法流程，利用公式(5)和(6)产生新模型，利用公式(7)控制温度，根据公式(8)判断模型是否接受，若当前模型被接受，则把接受的模型作为样本存储下来。

(3)计算后验概率分布。计算所积累的样本的似然函数，进而利用公式(3)计算相应的后验概率分布。

(4)缓慢降温，重复步骤(2)和(3)直至联合反演结束。

(5)提取解的表达。将后验概率密度均一化，提取出最大后验概率解，利用公式(4)提取均值解。

可以看出，完成一次贝叶斯联合反演便可以得到三种解：最终迭代解，最大后验概率解和均值解。

2联合反演的实现

图1　理论模型、共界面反演结果以及初始扰动结果Fig.1　Theoretical model, joint inversed results in common interface and initial perturbed results (color scale is the same as figure 2 and 5)a.电阻率理论模型(色标同图2和图5)；b.速度理论模型(色标同图2和图5)；c.电阻率共界面反演结果；d.速度共界面反演结果；e.电阻率初始扰动模型；f.速度初始扰动模型

本文采用于鹏等(2008,2009)提出的物性参数随机分布共网格建模技术以满足物性参数横向变化剧烈的复杂模型，该技术的本质为不同物性参数同界面分布这一联合反演基本条件。每个矩形网格都有独立的物性参数分布，按网格随机扰动的本质是将联合反演的基本条件推向网格单元的情况。该技术利用网格组成灵活地描述地质界面和块体，并可以根据正反演精度要求来调整网格间距。

二维大地电磁和地震方法的数值模拟分别利用通用的有限单元法和速度随机分布的射线追踪方法来实现．电阻率和速度同步联合反演具体的建模方法、正演公式、正则化因子选取，模拟退火算法等，在于鹏等(2009)和陈晓(2013)中有详细描述，本文不再赘述。

3模型试验

笔者设计了电阻率和速度分布非共界面的断块模型(图1a,b)。自上而下共五套电性层，电阻率分别为20.0,100.0,2000.0,5.0,200.0 Ω·m，四套速度层，速度分别为3.0，4.0，5.0，6.0 km/s，即最后两套地层的速度分布相同。大地电磁计算频点在320～0.000 55 Hz之间取40个，测点间距为1 km。

为了降低联合反演初始建模的人为性，选取均匀半空间为初始模型，依次经过单独反演、共网格联合反演和共界面联合反演，最后选取大地电磁和地震共界面联合反演结果(图1c,d)为模拟退火贝叶斯反演的先验信息。与以往贝叶斯反演不同，即往往仅考虑单一噪音条件下的不确定性，为了充分考虑数据噪音即数据空间不确定性的影响，本文共进行了10次随机高斯噪音条件下的贝叶斯联合反演。高斯噪音随机生成，标准差均为5%。

将电阻率和速度共界面联合反演结果分别扰动±20 %和±5 %作为模拟退火物性参数扰动空间(表1)。将模型中所有的物性参数，以及四套地层的底界面予以松约束，即放开反演。MT反演以TE模式为例。电阻率和速度贝叶斯联合反演的初始扰动结果见图1e,f。

表1　模拟退火物性参数扰动空间

图2　贝叶斯联合反演结果Fig.2　Inversion results of Bayesian joint inversiona.电阻率最大后验概率反演结果(标号1和2代表两种不同的随机高斯噪音条件)；b.速度最大后验概率反演结果；c.电阻率均值反演结果；d.速度均值反演结果；e.电阻率最终迭代结果；f.速度最终迭代结果

由于篇幅所限，仅列2次标准差均为5%的随机高斯噪音下的电阻率和速度反演结果。电阻率、速度最大后验概率反演结果见图2a，b，电阻率、速度均值反演结果见2c，d，电阻率、速度最终迭代结果见图2e，f。每一层的物性参数后验概率密度分布结果见图3。

最终迭代解是一般意义上的解，均值解和最大后验概率解是概率意义上的解。结合图2和图3，可以看出三种反演结果都可以接受。随着数据噪音的变化，最终迭代解差别较大，而最大后验概率解和均值解的一致性优于最终迭代解。从物性后验概率密度分布可以清楚的看出每层物性的概率分布情况，峰值处代表最大的可能，对解的不确定性表达更清晰。

图2和图3说明了不同噪音条件下不确定性反演解的差异，作者以第1和第50个测点为例，列出了标准差为5 %的十种不同的随机高斯噪音下的最大后验概率解对应的物性分布图(图4中黑色线条)。可以看出在不同噪音条件下，最大后验概率解不是某个确定值，而是一个变化的区间，这也是本文提出需要充分考虑数据不确定性空间对反演结果影响的原因，为了在不同的数据噪音条件下获取更可靠的解。将标准差为5 %的十种不同的随机高斯噪音下的所有可能解，进行了统一整理，并提取相应的最大后验概率解和均值解，列出第1和第50个测点在标准差为5 %的十种不同的随机高斯噪音下的最大后验概率解对应的物性分布图(图4中蓝色线条)。可看出蓝色线条是10种随机高斯噪音条件下的综合反映，体现了不同噪音的影响，体现了解的可能分布。综合了标准差为5 %的十种不同的随机高斯噪音的最大后验概率解和均值解见图5(色标同图1)，物性参数的后验概率分布见图6。

对比图2和图5的最大后验概率解和均值解，可看出图5的最大后验概率解表现的界面有一定“跳跃”，体现了数据不确定性对反演结果的影响。均值解为所有样本的综合反映，界面横向比较连续，解比较稳定。对比图6和图3的物性参数后验概率密度，可看出图6后验概率密度函数的峰值更加集中，说明样本采样更充分，体现了不同数据噪音的影响。

4结论

本文充分考虑数据噪音的不确定性，实现了大地电磁和地震数据同步贝叶斯联合反演，表达了解的可能分布，提取了最大后验概率解和均值解，增强了对反问题解的表达和认识。通过本文论述，可以看出：

(1)基于模拟退火算法实现大地电磁和地震数据贝叶斯联合反演是可行的且是有优势的。其一，模拟退火本身借鉴了概率的思想，模型接受的准则本身就是和概率相关的；其二，模拟退火是全局寻优算法，需要进行大量的正演计算和模型判断，为样本的收集提供了便利。

(2)在不同数据噪音条件下，最大后验概率解，均值解以及最终迭代解都有较好地表达反演结果的能力，任何一种解都可以从不同角度作为一种地球物理认识为地质解释提供模型，为地质学家提供更广阔的视野。数据空间的不确定性对反演结果影响比较大，在实际资料处建议重视数据噪音的评估和利用。

参考文献

陈建江, 印兴耀, 张广智.2006.层状介质AVO叠前反演[J]. 石油地球物理勘探,41( 6):656-662.

陈晓，于鹏，张罗磊，等. 2010.大地电磁与地震正则化同步联合反演[J].地震地质,32(3):402-408.

陈晓，于鹏，张罗磊，等.2011.地震与大地电磁测深数据的自适应正则化同步联合反演[J].地球物理学报,54(10):2678-2681.

陈晓.2013.大地电磁与地震正则化联合反演研究[D].上海:同济大学.

郭荣文.2011.贝叶斯MT反演的非线性和不确定度分析[D].湖南：中南大学.

刘文劼.2012.频率和空间相关性观测数据条件下的大地电磁贝叶斯反演[D].湖南：中南大学.

王文涛.2008.地震储层评价与预测的贝叶斯反演方法研究[D].武汉：中国地质大学.

杨迪琨,胡祥云.2008.含噪声数据反演的概率描述[J].地球物理学报,51(3):901-907.

杨辉.2000.大地电磁、地震资料的模拟退火约束联合反演[D].上海:同济大学.

于鹏,戴明刚,王家林，等.2008.密度和速度随机分布共网格模型的重力与地震联合反演[J]. 地球物理学报,51(3):845-852.

图3　各层物性参数后验概率密度Fig.3　The posterior distribution of parameters in every layera.第一层电阻率和速度概率密度分布(标号1-4代表两种不同的随机高斯噪音条件下电阻率和速度后验概率密度)；b.第二层电阻率和速度概率密度分布；c.第三层电阻率和速度概率密度分布；d.第四层电阻率概率密度分布；e.第五层电阻率概率密度分布

图4　十种不同的随机高斯噪音下反演所得最大后验概率解的物性分布Fig.4　The posterior distribution of parameters for ten different kind of random Gaussian noisea.第1个测点电阻率分布；b.第1个测点速度分布；c.第50个测点电阻率分布；d.第50个测点速度分布

图5　十种不同的随机高斯噪音下最大后验概率解和均值解Fig.5　The maximum posterior and the mean solutions for ten different kind of random Gaussian noisea.电阻率最大后验概率解；b.速度最大后验概率解；c.电阻率均值解；d.速度均值解

图6　十种不同的随机高斯噪音下物性参数后验概率密度Fig.6　The posterior distribution of parameters in every layer for ten different kind of random Gaussian noisea.第一层电阻率后验概率密度；b.第一层速度后验概率密度；c.第二层电阻率后验概率密度；d.第二层速度后验概率密度；e.第三层电阻率后验概率密度；f.第三层速度后验概率密度；g.第四层电阻率后验概率密度；h.第五层电阻率后验概率密度

于鹏,戴明刚,王家林，等.2009.电阻率和速度随机分布的MT与地震联合反演[J]. 地球物理学报,52(4):1089-1097.

张霖斌,姚振兴,纪晨，等.1997,快速模拟退火算法及应用[J],石油地球物理勘探,32(5):654-660.

赵峦啸.2010.贝叶斯框架下基于叠前地震数据的岩性(流体)预测[D].上海：同济大学.

Buland A, ElOuair Y.2006. Bayesian time-lapse inversion[J]. Geophysics, 71(3):R43-R48.

Chen J, Hoversten G M.2009.Joint inversion of marine seismic AVA and CSEM data using statistical rock-physics models and Markov random fields[C]. SEG, 714-718.

Hoversten G M, Cassassuce F, Gasperikova E, et al. 2006.Direct reservoir parameter estimation using joint inversion of marine seismic AVA and CSEM data[J]. Geophysics, 71(3):C1-C13.

Ingber L.1989. Very fast simulated reannealing[J]. Math Conput Modeling,967-973.

Kane J A, Rodi W, Kennth P, et al..2004. Structural uncertainty and Bayesian inversion[C]. SEG：1511-1514.

Tarantola A.1987.Inverse Problem Theory and Methods for Model Parameter Estimation[M].Elsevier Science

Zhdanov M S. 2002. Geophysical inverse theory and regularization problems[M]. Netherlands:Elsevier Science,36-45.

Bayesian Joint Inversion of Magnetotelluric and Seismic Data Based on Simulated Annealing Method

CHEN Xiao1,2,ZHOU Lei3,YU Peng4,DENG JU-zhi1，YI Han-ting5

(1. Fundamental Science on Radioactive Geology and Exploration Technology Laboratory, East China University of Technology, Nanchang, JX 330013,China; 2. Hubei Subsurface Multi-scale Imaging Key Laboratory (SMIL), China University of Geosciences, Wuhan, HB 430074,China;3.Key Laboratory of Exploration Technologies for Oil and Gas Resources, Ministry of Education, Yangtze University, Wuhan, HB 430100,China; 4.State Key Laboratory of Marine Geology,Tongji University, Shanghai 200092,China;5.Shengli oilfield well logging company, Dongying, SD 257000,China)

Abstract:We are successful in fulfilling the synchronous Bayesian joint inversion of MT and seismic based on simulated annealing method. Bayesian inversion is the hotspot of inversion research. By now, lots of papers about Bayesian joint inversions of other geophysical methods are published, but less work in the joint inversion of MT and seismic. This thesis considered the effect of noise uncertainty, expressed the solution of joint inversion as posterior probability density distribution, and extracted the maximum posterior solution, the mean solution and final iterative solution of resistivity and velocity.

Key Words:magnetotelluric; seismic; joint inversion; bayesian inversion; simulated annealing; noises

中图分类号：P631

文献标识码：A

文章编号：1674-3504(2016)01-0059-08

doi:10.3969/j.issn.1674-3504.2016.01.010

作者简介：陈晓(1986—)，男，博士，研究方向为综合地球物理，特别是多方法间的联合反演，E-mail:dwjhtj@hotmail.com通讯作者：周磊(1985—)，男，博士，研究方向为电磁法理论及应用技术，Email:zhoulei-85-7-8@163.com

基金项目：同济大学海洋地质国家重点实验室开放基金(MGK1506), 中国地质大学(武汉)地球内部多尺度成像湖北省重点实验室开放基金(SMIL-2015-09)，东华理工大学放射性地质与勘探技术国防重点学科实验室开放基金(RGET1405)和东华理工大学博士科研启动基金(DHBK201403)联合资助。

收稿日期：2014-09-05