基于ASSA方法的等值性电测深数据反演及数据误差相关性分析

2017-11-01付代光廖锦芳周黎明肖国强王法刚

石油地球物理勘探 2017年5期

关键词：单纯形模拟退火等值

付代光廖锦芳周黎明肖国强王法刚罗荣

(①长江科学院水利部岩土力学与工程重点实验室，湖北武汉 430010; ②湖北省页岩气开发有限公司，湖北武汉 430071)

·非地震·

基于ASSA方法的等值性电测深数据反演及数据误差相关性分析

付代光*①廖锦芳②周黎明①肖国强①王法刚①罗荣①

(①长江科学院水利部岩土力学与工程重点实验室，湖北武汉 430010; ②湖北省页岩气开发有限公司，湖北武汉 430071)

付代光，廖锦芳，周黎明，肖国强，王法刚，罗荣.基于ASSA方法的等值性电测深数据反演及数据误差相关性分析.石油地球物理勘探，2017,52(5):1077-1084.

针对直流电测深数据等值性反演的病态性(等值性层参数及层参数间的狭长平坦的强相关性)，尝试采用自适应单纯形模拟退火(ASSA)方法反演具等值性电测深模型，通过与Lamarckian-marked-constraint(LMC)算法对比，发现ASSA方法能较好地解决电测深等值性反演问题，且反演精度较高。目前地球物理反演多基于数据误差为高斯分布和空间不相关(非对角元素为零)的假设，对贝叶斯反演而言，后一假设如果运用不当往往会低估反演结果的不确定性，从而导致对反演结果可靠性的误判。为了合理评价反演结果的不确定性，文中采用Runs-test非参数化判别方法，通过对等值性模型的归一化数据残差的判别，确定数据误差是空间相关的。由此得出的协方差矩阵应用于贝叶斯反演，进而对反演结果的不确定性、相关性等进行更合理、准确的评价。

电测深等值性贝叶斯反演 ASSA 数据误差相关性不确定性相关性

1 引言

直流电测深方法是直流电阻率法的一种，通过控制极距达到探测不同深度目标体的目的。该方法应用范围广，不但可以应用于寻找金属和非金属矿产、圈定油藏边界、探测天然气水合物，而且在水文与工程地质、城市环境调查、建筑基础勘查等方面也发挥着重要的作用[1-5]。

直流电测深具有等值性，而具有等值性的电测深数据反演是病态的[6]。电测深等值性是由其方法本身的特性所决定的，无法从理论上消除，因而常采用改进的反演算法提高反演精度，但相关的研究并不多见。Basokur等[7]提出了混合算法Lamarckian-marked-constraint(LMC)，该算法对传统的遗传算法进行了改进，将自然选择方式变为人为选择，即增加约束条件，对不满足条件的模型进行淘汰，同时将最小二乘与标记约束遗传算法相结合，利用最小二乘的梯度信息，提高算法对强相关性参数反演的能力。但该算法只能应用于有限个参数的反演中，另外对不符合实际情况的约束条件，反演算法将会失效；同时，该算法的反演过程消耗大量的计算时间，尤其二维反演更是如此，因为二维正演过程需要大量的计算时间，网格越小所需时间就越多，因而限制了在实际中的应用。为了实现基于LMC方法的二维电阻率成像， Akça等[8]采用Delaunay 三角剖分法极大地减少了地电模型的网格数，从而基于反演结果实现了岩性单元的划分。Fallat等[9]和Dosso等[10]提出的单纯形模拟退火(SSA)和自适应单纯形模拟退火是一种局部与全局、线性与非线性相结合的确定性搜索算法，该方法对解决强参数相关性问题非常有效，而且无需求导计算。针对ASSA方法的特点，本文尝试采用该方法对等值性电测深数据进行反演。

贝叶斯反演是一种把观测数据和模型参数都看成随机变量的反演方法。传统的反演过程是一个寻找最优解的过程，它只提供一个确定的反演结果，但并未定量描述反演结果的可靠度，也未将含有误差和噪声的数据加入到反演过程中；而贝叶斯反演能对上述问题提供较满意的解决方法[11-13]。贝叶斯反演是关于协方差矩阵最大似然函数的求解过程，如果假设协方差矩阵中非对角元素为零，那么将得到贝叶斯反演中常规的目标函数公式(非对角元素为零)。目前很多地球物理反演都是基于此假设的[14-18]，但该假设往往造成对反演结果不确定性的低估。为了对反演结果的不确定性进行合理的评价，Dosso等[19，20]提出了一种非参数化方法Runs-test(亦可称为Runs检验)用于检验协方差矩阵元素是否空间相关。通过Runs检验得到的协方差矩阵应用于贝叶斯反演中，能对反演结果的不确定性、相关性等给出合理的评价。

本文首先从正演角度对等值性模型的一维和二维参数变化引起的误差变化程度曲线的特征进行了分析，并应用ASSA算法对等值性电测深模型进行反演，通过与LMC方法反演结果的对比验证了ASSA算法的有效性。为了分析数据误差相关性对反演结果的不确定性的影响，本文对比了数据误差空间相关与不相关反演结果的差异，并采用非参数化方法Runs检验数据误差相关与否。将Runs检验判定出的数据误差相关性结果应用到贝叶斯反演中，从而对反演结果的不确定性、参数间相关性等进行更加合理的评价。

2 贝叶斯反演

2.1 贝叶斯定理

假定d表示实测数据，m为模型参数，则对于任意的d和m，满足如下的贝叶斯关系

(1)

式中：p(m|d)表示已知d时，所确定的m的函数，称为后验概率密度函数(probability density function，简称为后验PDF)；同理，p(d|m)表示在给定m时，所得d的概率密度函数；p(m)称为m的先验PDF；p(d)是d的先验概率密度，通常d是已知的，因此该项往往都看作常数，此时式(1)可写为

p(m|d)∝p(d|m)p(m)

(2)

p(d|m)与似然函数L(m)是对等的，即L(m)∝p(d|m)，因此式(2)可进一步写为

p(m|d)∝L(m)p(m)

(3)

式中L(m)∝exp[-E(m)]，E(m)为目标函数。广义目标函数包含了数据误差p(m)和先验信息，即

φ(m)=E(m)-lnp(m)

(4)

联立式(1)～式(4)，并将向量m的PDF归一化得到

(5)

式中V表示模型的积分域。式(5)即所求问题的解。为了对反演结果进行更合理的解释，需要对反演参数进行估计、分析不确定性并计算参数相关性，即对后验概率进行积分求解，具体求解公式参见文献[21]。

2.2 自适应单纯形模拟退火

ASSA方法是将模拟退火与单纯形方法相结合的一种非线性方法，该方法具备全局与局部搜索能力，与传统的全局优化方法(如模拟退火和遗传算法等)和局部优化方法(如拟牛顿法、最小二乘法等)相比具有以下特点：

(1)反演结果与初始值无关，克服了局部算法对初始估计值的较高要求；

(2)能够保留搜索过程中得到的最优解，克服了模拟退火不能保存最优解的缺点；

(3)由于单纯形具有梯度搜索能力(但不需要计算导数或灵敏度矩阵)，因此该方法反演效率高，尤其对相关性强或在收敛值附近的问题能够高效解决。

自适应单纯形模拟退火的自适应体现在对参数扰动的自动更新，其模型扰动公式为

(6)

ASSA方法将单纯形嵌入到模拟退火中，因此每一次的降温过程不再是对单一模型进行操作，而是对所有模型的单纯形进行操作。通过模拟退火扰动产生的初始模型，经由单纯形算法进行优化，再由模拟退火对优化参数进行随机扰动，然后利用Metropolis准则判断新生成的模型是否被接受。在上述过程中需要注意的是，单纯形会产生M+1个模型(M表示模型m的元素个数)，且模型每次都会被保留，因此避免了模拟退火中好的模型被丢弃；同时，在模拟退火后期，随着温度降低(温度为0°C时)，会变成纯粹的单纯形的搜索模式[8]。

3 协方差计算

对贝叶斯反演中数据协方差矩阵的计算采用Dosso[19]提出的计算方式。数据协方差矩阵可定义为

Cd=〈(d-〈d〉)T(d-〈d〉)〉

(7)

式中〈·〉表示总体平均。由于采集数据的误差不能代表整个样本的误差，因此该公式很难应用于实际。

如果将协方差矩阵近似为最常用的形式——对角协方差矩阵，即非对角元素为零的矩阵Cd=σ2I(I是单位矩阵)，那么此时的非线性贝叶斯公式的似然函数(目标函数)可简化为常用的二范数形式。对角协方差矩阵其实忽略了非对角元素的相关性，因而容易导致对参数不确定性的低估。为了合理估计参数的不确定性，本文采用优化模型的数据残差计算全元素协方差矩阵

(8)

i=1,…,N;j=1,…，N

(9)

(10)

4 ASSA反演

为验证ASSA算法的有效性，本文选择具有等值性特征的四层地电模型[7](表1)。表中第一行给出了模型各层的电阻率和厚度。在验证ASSA方法之前，先对该模型进行正演分析。

表1 模型参数值及LMC与ASSA反演结果

4.1 电测深正演

参数敏感性与参数相关性是影响反演结果的关键因素。单个参数变化引起的适值变化程度，称为敏感性[10]。如果随着参数变化，适值发生较明显改变(曲线的曲率较大)，则说明参数敏感，较窄的参数敏感性有利于反演；反之，则不利于参数的反演。图1给出了第二、第三层电阻率和厚度的参数适值范围。从图中可以看出各参数相对敏感，在不考虑参数相关性、等值性等其他因素的情况下，各参数的反演并不会太复杂。

从图2部分参数间的误差能量等值线图可知，参数间存在强相关性，且这种强相关性表现为狭长型。其中，ρ2-h2(图2g)、ρ3-h3(图2h)表现得尤为明显，因此对其反演时更难捕捉到有效的参数值。因此，单纯采用线性反演方法需要有较好的初始猜测值，否则对这种狭长相关性的参数反演可能会失败；而对非线性全局优化反演方法而言，则要求大量空间搜索(大量的计算时间)，或许可以得到较为合理的解，有时也可能得到局部最优解(对等值性电测深而言，可能为等值性参数对)。因此该模型的反演具有一定的挑战性。

图1 一维参数误差能量断面图

图2 误差能量等值线图(黑色区域为低适值区)

4.2 电测深反演

为验证ASSA算法的有效性，本文选择Basokur等[7]提出的LMC反演方法进行对比。LMC法是一种将线性最小二乘法与遗传算法相结合的混合算法。为了方便对比，设定与文献[7]相同的参数搜索范围(见表1)。ASSA反演参数为：初始温度为0.03°C，迭代1500次，每个温度扰动5次。图3为ASSA反演的适值变化与迭代次数的关系，从图中可以看出适值降低速度较快，即算法的收敛速度较快。根据图4所示的反演结果可知：视电阻率—极距和电阻率—深度曲线的反演结果精度较高，具体的反演值见表1。由表1可见，ASSA的反演精度与LMC方法相当甚至更高。但需要注意的是，单纯形采用近似梯度搜索方法，而对于等值性强且窄的相关性参数而言，有时单次反演也未必一定能取得很好的结果，一般取几次反演结果的均值为宜。

图3 ASSA反演的适值与迭代次数的关系曲线

5 数据误差分析

数据误差的相关与否，决定了目标函数的形式，而目标函数的形式是否会对反演结果产生影响，为此，采用文献[7]的模型，ASSA反演参数设定为：迭代次数为500次，初始温度为0.3°C，每个温度扰动20次，均取2次反演结果的平均值。图5为电阻率反演结果。从图5b看出，空间相关协方差矩阵对反演结果的拟合精度影响较大，其中对第三层电阻率和厚度的影响最为明显；但从图5a发现，视电阻率—极距的反演拟合精度较好，究其原因，可能是由数据误差的空间相关性及模型的等值性引起电阻率—深度剖面的反演拟合效果较差。在对电阻率剖面进行反演时，如何判断数据误差是否空间相关，至关重要。

本文采用非参数化Runs方法判断数据误差是否相关。从图6a可以看出，对角协方差矩阵的自相关图形主峰较宽，而且跳跃幅度较大，因此初步判定数据误差相关，通过Runs定量计算可知p=0.0001，再次证明数据误差相关。由图6b也可看出，自相关的波峰较窄，且p=0.94，说明数据误差是空间相关的。需要注意的是，对于不同轮次的反演结果，数据残差自相关图形和p值大小会有所变化，但这并不影响最终的结论，因此不同轮次的反演对残差自相关结果的估计不敏感。

前文已证明数据误差空间相关，因此在贝叶斯反演和后验概率积分时都应将其考虑进去。图7～图9为数据误差相关时所得的一维边缘概率分布。由图7可见，ρ3和h3的不确定性较大，呈双峰状，且峰值均不在理论值附近，ρ3和h3在其上、下边界附近的采样值较多；h1和ρ4主峰较窄，h1基本在9～10间变化，ρ4基本在30～50间变化，且与理论值基本一致，不确定性小；ρ2和h2的峰值较缓，ρ2在20～80间变化，h2在5～20间变化，置信区间相对较小，不确定性相对较大。根据图8可知，ρ3-ρ2、ρ3-h2、h3-ρ2、h3-h2的二维边缘概率分布的采样点极其分散，这是由于ρ3和h3为多峰多模态，ρ2和h2的不确定性较大，置信区间较小，且ρ2、h2、ρ3和h3彼此间存在不同程度的相关性(图9)，影响了各参数的采样值；对于h1和ρ4这样主峰较窄的参数，ρ3-h1、ρ3-ρ4、h3-h1、h3-ρ4的二维边缘概率分布的采样点相对集中。结合图8和图9综合分析发现，ρ2-h2呈强线性正相关，而ρ3-h3呈强非线性负相关，此结果与图2中的结果也是吻合的。这些相关性将会影响反演结果的精度，同时也会影响参数的不确定性。强的相关性会导致反演精度降低，增加反演结果的不确定性，这也是ρ3和h3在反演过程中很难得到较好结果的原因。如何有效降低参数间的相关性、提高反演精度，是亟待解决的问题。

图4 采用文献[7]参数的ASSA反演结果

图5 电阻率反演结果

图6 Runs-test法反演结果

图7 部分参数一维边缘概率分布

图8 部分参数二维边缘概率分布白色虚线表示参数真值

图9 模型参数间相关系数

6 结论

(1)通过对等值性电测深数据一维和二维正演分析发现，单个参数敏感性较窄，利于参数反演。但从二维正演曲线发现，大多数参数间存在强相关性，且在真值附近存在许多等值性参数对，因此要获得好的反演结果较困难。本文根据ASSA方法的近似梯度搜索能力和收敛速度快等特点，对等值性电测深数据进行了反演，其反演结果精度较高，验证了ASSA算法的优越性。

(2)通过对比数据误差空间相关与不相关的反演结果发现，数据误差相关与否对反演精度具有较大的影响。为了对数据误差的相关性进行判断，本文采用Runs检验法定性和定量地判定了电测深模型数据残差的相关性，最终证明其数据误差是空间相关的。

(3)通过对后验概率数值积分，获得反演结果的一维和二维边缘概率分布、相关系数等特征量，由此得到本文所述模型的第二层与第三层参数的不确定性大、且第三层的厚度和电阻率呈多模态性(双峰状)的结论。参数间存在不同程度的相关性，其中以第二和第三层相关性最为明显，第三层厚度和电阻率呈强非线性负相关，第二层参数间呈强正相关。相关性影响了反演结果精度并增大了参数的不确定性。

[1] 王桥，万汉平，王闻文等.综合物探方法在铝土矿勘查中的应用.地球物理学进展,2012,27(2):709-714. Wang Qiao，Wan Hanping，Wang Wenwen et al.The application of integrated geophysical exploration in Bauxite.Progress in Geophysics,2012,27(2):709-714.

[2] 翁爱华，祝嗣安.天然气水合物的近海底直流电测深响应.石油地球物理勘探,2010,45(3)：458-461. Weng Aihua,Zhu Si’an.Near sea-bottom DC sounding response for gas hydrate.OGP,2010,45(3)：458-461.

[3] 张天伦.用直流电阻率法确定油气藏边界的初步实验.石油地球物理勘探,1990，25(5)：584-593. Zhang Tianlun.An experiment in locating reservoir boundary using direct-current resistivity method.OGP，1990，25(5)：584-593.

[4] 王玉，郭秀军，贾永刚等.土壤—地下水系统非水相液体污染区的电性特征及电阻率法探测.地球物理学进展,2009,24(6):2316-2323. Wang Yu，Guo Xiujun，Jia Yonggang et al.Overview of the electrical characters and the electrical survey methods for the soil-groundwater system in organic contaminated areas.Progress in Geophysics,2009,24(6):2316-2323.

[5] 贾志宽，安西峰，李振峰等.非常规电测深法在工程地质勘查中的应用.地球物理学进展,2006,21(1):296-299. Jia Zhikuan,An Xifeng,Li Zhenfeng et al.Application of unconventional electric detecting method in engineering geology survey.Progress in Geophysics,2006,21(1):296-299.

[6] Alvarez J P F,Martinez J L F,Perez C O M.Feasibility analysis of the use of binary genetic algorithms as importance samplers application to a 1-D DC resistivity inverse problem.International Association for Mathematical Geology,2008,40(4):375-408.

[7] Basokur A T,Akca I,Siyam N W A.Hybrid genetic algorithms in view of the evolution theories with application for the electrical sounding method.Geophy-sical Prospecting,2007，55(3)：393-406.

[8] Akça I,Basokur A T.Extraction of structure-based geoelectric models by hybrid genetic algorithms.Geophysics,2010，75(1):F15-F22.

[9] Fallat M R,Dosso S E.Geoacoustic inversion via local,global,and hybrid algorithms.Acoustical Society of America.1999，105 (6):3219- 3230.

[10] Dosso S E,Wilmut M J,Lapinski A L S.An adaptive-hybrid algorithm for geo-acoustic inversion.Acous-tical Society of America，2001，26(3):324-336.

[11] Guo R W,Dosso S E,Liu J X.Non-linearity in Baye-sian 1-D magnetotelluric inversion.Geophysical Journal International，2011，185(2):663-675.

[12] Dong H F,Dosso S E.Bayesian inversion of interface-wave dispersion for seabed shear-wave speed profiles.IEEE Journal of Oceanic Engineering,2011,36(1):1-11.

[13] 尹彬,胡祥云.基于并行回火技术优化的大地电磁数据贝叶斯变维反演.石油地球物理勘探,2016,51(6):1212-1218. Yin Bin,Hu Xiangyun.Magnetotelluric Bayesian trans-dimensional inversion improved by parallel tempering algorithm.OGP,2016,51(6):1212-1218.

[14] 宗兆云,印兴耀,张繁昌.基于弹性阻抗贝叶斯反演的拉梅参数提取方法研究.石油地球物理勘探,2011,46(4):598-604,609. Zong Zhaoyun,Yin Xingyao and Zhang Fanchang.Elastic impedance Bayesian inversion for Lame parameters extracting.OGP,2011,46(4):598-604,609.

[15] 张丰麒，金之钧，盛秀杰等．一种稳健的弹性阻抗反演方法．石油地球物理勘探，2017，52(2)：294-303. Zhang Fengqi,Jin Zhijun,Sheng Xiujie et al.A stable elastic impedance inversion approach.OGP,2017，52(2)：294-303

[16] 张广智，潘新朋，孙昌路等．纵横波联合叠前自适应MCMC反演方法．石油地球物理勘探，2016，51(5)：938-946. Zhang Guangzhi,Pan Xinpeng,Sun Changlu et al.PP- & PS-wave prestack nonlinear inversion based on adaptive MCMC algorithm.OGP,2016，51(5)：938-946.

[17] 陈建江,印兴耀.基于贝叶斯理论的AVO三参数波形反演.地球物理学报,2007,50(4):1251-1260. Chen Jianjiang,Yin Xingyao.Three parameter AVO waveform inversion based on Bayesian theorem.Chinese Journal of Geophysics,2007,50(4):1251-1260

[18] 侯栋甲，刘洋，胡国庆等．基于贝叶斯理论的叠前多波联合反演弹性模量方法．地球物理学报，2014,57(4)：1251-1264. Hou Dongjia,Liu Yang,Hu Guoqing et al.Prestack multiwave joint inversion for elastic moduli based on Bayesian theory.Chinese Journal of Geophysics,2014,57(4):1251-1264.

[19] Dosso S E,Nielsen P L,Wilmut M J.Data error covariance in matched-field geoacoustic inversion.A-coustical Society of America，2006，119(1):208-219.

[20] Dosso S E,Holland C W.Geoacoustic uncertainties from viscoelastic inversion of seabed reflection data.IEEE Journal of Oceanic Engineering,2006,31(3):657-671.

[21] 付代光，刘江平，周黎明等.基于贝叶斯理论的软夹层多模式瑞雷波频散曲线反演研究.岩土工程学报，2015,37(2):321-329. Fu Daiguang,Liu Jiangping,Zhou Liming et al.Inversion of multimode Rayleigh-wave dispersion curves of soft interlayer based on Bayesian theory.Chinese Journal of Geotechnical Engineering,2015,37(2):321-329.

(本文编辑：刘海樱)