E-mail Alert Rss
 

物探与化探, 2024, 48(1): 67-76 doi: 10.11720/wtyht.2024.2583

方法研究·信息处理·仪器研制

海洋可控源电磁法与地震全波形二维联合反演研究

孔繁祥,1, 谭捍东2, 刘建勋1

1.中国地质科学院 地球物理地球化学勘查研究所,河北 廊坊 065000

2.中国地质大学(北京) 地球物理与信息技术学院,北京 100083

Two-dimensional joint inversion based on the marine controlled-source electromagnetic method and seismic full-waveform

KONG Fan-Xiang,1, TAN Han-Dong2, LIU Jian-Xun1

1. Institute of Geophysical and Geochemical Exploration, Chinese Academy of Geological Sciences, Langfang 065000, China

2. School of Geophysics and Information Technology, China University of Geosciences, Beijing 100083, China

责任编辑: 叶佩, 沈效群

收稿日期: 2022-11-22   修回日期: 2023-03-27  

基金资助: 中国地质调查局地质调查项目“柴达木盆地盐湖区物探综合调查”(DD20230298)
中国地质科学院地球物理地球化学勘查研究所所长基金基本科研业务费项目“现代地质勘查工程技术集成与创新”(A2022P02)

Received: 2022-11-22   Revised: 2023-03-27  

作者简介 About authors

孔繁祥(1995-),男,硕士,助理工程师,从事地震数据处理以及算法开发工作。Email:kfanxiang@mail.cgs.gov.cn

摘要

为降低单一地球物理反演方法的局限性以及反演中所存在的多解性等问题,开展了海洋可控源电磁法(MCSEM)与地震全波形的二维联合反演研究。MCSEM采用数据空间OCCAM反演算法,地震全波形采用梯度法反演算法,引入交叉梯度函数实现两种物性参数结果的相互耦合;开发出一套二维联合反演算法,并通过3组理论模型算例验证了算法的准确性。研究结果表明:MCSEM联合反演结果相对于单方法反演结果有显著改善和提升,主要体现在异常体的形态刻画、结构构造以及物性数值的恢复,这说明全波形的反演方法能够提高MCSEM反演结果的可靠性。

关键词: 海洋可控源电磁法; 地震全波形; 交叉梯度; 联合反演

Abstract

To reduce the limitations and the multiplicity of solutions of a single geophysical inversion method, this study investigated the two-dimensional joint inversion based on the marine controlled-source electromagnetic (MCSEM) method and seismic full-waveform inversion. The MCSEM method employs the data-space Occam’s algorithm, while the seismic full-waveform inversion utilizes the gradient algorithm. By incorporating a cross-gradient function for the mutual coupling of the two types of physical property parameters, this study developed a two-dimensional joint inversion method, whose accuracy was verified using three different models. As indicated by the results, compared to a single inversion method, the MCSEM-based joint inversion yielded significantly improved inversion results, predominantly in terms of the morphology characterization of anomalous bodies, as well as the reconstruction of their structure and textures and their physical property values. Therefore, the full-waveform inversion can enhance the reliability of the MCSEM inversion results.

Keywords: marine controlled-source electromagnetic (MCSEM); seismic full-waveform; cross-gradient; joint inversion

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

本文引用格式

孔繁祥, 谭捍东, 刘建勋. 海洋可控源电磁法与地震全波形二维联合反演研究[J]. 物探与化探, 2024, 48(1): 67-76 doi:10.11720/wtyht.2024.2583

KONG Fan-Xiang, TAN Han-Dong, LIU Jian-Xun. Two-dimensional joint inversion based on the marine controlled-source electromagnetic method and seismic full-waveform[J]. Geophysical and Geochemical Exploration, 2024, 48(1): 67-76 doi:10.11720/wtyht.2024.2583

0 引言

地球物理联合反演是利用多种方式采集所取得的数据来实现多种物性参数反演的方法[1]。在综合地球物理勘探中,联合反演是一种重要的定量解释手段,近些年的大量研究结果表明,联合反演方法能够较好地压制干扰信息(如背景噪声),减小地球物理反演问题中的多解性,获取更为可靠、精准的反演结果,可进一步提高反演解释的精度[2-4]。国外学者对于海洋可控源电磁法(marine controllable source electromagnetic method,MCSEM)和地震资料的联合反演有一定的研究:Hu等[5]成功实现了基于交叉梯度的MCSEM资源和地震走时资料的交互式联合反演(alternating joint inversion),成像精度有显著提升;Brown等[6]实现了地震全波形资料约束MCSEM资料的二维Occam反演,利用澳大利亚Pluto气田数据反演结果验证了该方法的有效性;Da Silva等[7]实现了使用地震全波形资料约束的MCSEM三维反演,并取得了较好的成果;Mac-Gregor等[8]应用Archie公式来估算西非某海域海底底层的电阻率和地震波速度之间的经验关系。然而,这些研究主要还是基于岩石物理关系的同步联合反演,这种方法过于依赖电阻率和地震波波速之间的岩石物理关系[9]

而国内对此研究较少,近年来只有杜润林[10]、李刚[11]在博士论文中进行了相关研究:杜润林提出基于储层构造约束的MCSEM和地震全波形联合反演方法,利用非线性共轭梯度方法推导出地震和电磁梯度;李刚采用构造JTV(joint total variation)耦合泛函成功开发了MCSEM与地震资料构造一维和二维联合反演算法。

本文在前人提出的交叉梯度的基础上,建立了基于交叉梯度函数约束的MCSEM和地震全波形反演的目标函数,详细推导了梯度表达式,实现地下电阻率和地震波速度这两种物性参数之间相互约束的联合反演算法,并通过理论模型的测试,探讨了联合反演算法的可行性及适用性。

1 地震全波形和海洋可控源电磁二维正反演方法

1.1 地震全波形二维正反演方法

本文使用Hustedt等[12]提出的九点混合网格有限差分方法来完成地震波场的数值模拟计算,并且利用了Jo等[13]提出的方法来构造混合网格条件下的差分方程:对正交坐标系与45°旋转坐标系的二阶空间差分算子进行优化结合,同时,采取了完美匹配层吸收边界条件,又可称之为PML(perfectly matched layer)。PML边界条件最早是由Berenger[14]所提出,该方法是利用设置在目标网格边界周围的PML层数,在PML层数中加入衰减函数可吸收传至人工边界的反射波,进而消除边界反射波的干扰,在实际计算中需要划分若干网格,定义阻尼函数,在PML内部能量衰减。

全波形反演是以最小二乘理论为基础实现数据拟合,反演过程可概括为以下3个主要步骤[15-16]:第一步是获取正演响应(即理论地震波场值),根据所构建的初始速度模型,通过数值模拟方法计算而得到;第二步是建立目标函数,基于正演所计算的波场与观测波场之间的波场残差值来实现;第三步是求解反演问题,利用梯度法并结合Pratt[17]所提出的频率域全波形反演方法,可节省内存空间并提升反演速度。

1.2 海洋可控源电磁二维正反演方法

MCSEM利用总场=一次场(Ep,Hp)+二次场(Es,Hs)的方法求解正演响应。一次场的计算采用Ward和Hohmann[18]提出的方法,并且参照文献[19]中关于一次场计算公式的详细推导过程,推导了任意方向的电偶源作用下水平层状介质中各个电磁场分量的计算公式。

对于MCSEM中二次场的计算,本文使用基于四边形网格剖分的有限单元法[20]:在二维空气—海水—海底岩石的层状模型中,设走向方向为y,水平电偶源位于海底上方的海水中,产生如下关于 EysHys 的偏微分方程:

xy︿ke2E︿ysx+ zy︿ke2E︿ysz- y︿E︿ys+ xikyke2H︿ysz- zikyke2H︿ysx=- xikyσske2E︿xp- zikyσske2E︿zpsE︿yp- xz︿ke2H︿ysx- zz︿ke2H︿ysz+ z︿H︿ys+ xikyke2E︿ysz- zikyke2E︿ysx=- xz︿σske2E︿zp+ zz︿σske2E︿xp,

式中:z︿=iμ0ω为阻抗率;y︿=σ是忽略位移电流情况下的导纳率;ke2=ky2-k2,k2=-iμ0ωσ=-z︿·y︿;σ为介质电导率;σs为异常体电导率。

利用加权余量法对有限元方程进行了推导。参照Mitsuhata[20],本文同样使用伽辽金法,将格林公式加入到式(1)中,对所有网格单元进行叠加,并考虑到如下边界条件:在内边界上,电场和磁场的切向分量(E︿tH︿t)是连续的,而求取积分实质上是沿着每个边界都进行了两次方向相反的积分,两者会相互抵消,这使得内边界的积分之和为0;而对于外边界,采用了Dirichlet边界条件。经上述推导,最终得到正演所要求解的有限元方程组:

e=1NeDeNixy︿ke2E︿ysx+Nizy︿ke2E︿ysz+Niy︿E︿ys+Nixikyke2H︿ysz-Nizikyke2H︿ysxdxdz=- e=1NeDeNiσsE︿yp+Nixikyσske2E︿xp+Nizikyσske2E︿zpdxdz,
e=1NeDeNixz︿ke2H︿ysx+Nizz︿ke2H︿ysz+Niz︿H︿ys-Nixikyke2E︿ysz+Nizikyke2E︿ysxdxdz= e=1NeDeNixz︿σske2E︿zp-Nizz︿σske2E︿xpdxdz,

式中:Ne表示网格单元的总数。

数据空间Occam反演理论是由Siropunvarapor W与Egbert G[21]在传统Occam反演方法的基础上提出来的,该方法改进了模型空间Occam反演方法,是通过数学方法将模型迭代更新公式的计算和储存从模型空间转换到数据空间,这种数学转换不仅降低了计算量和存储空间,同时也提高了计算效率。

数据空间Occam反演流程如下:首先,输入观测参数以及初始模型(ρ0),计算雅可比矩阵(Jk),这里使用自动搜索的方式获取了拉格朗日算子(λ),并计算当前模型响应的数据拟合差(RMS),同时选择RMS最小时所对应的λ值;然后,根据当前的λ值更新电阻率模型。为提高计算效率,选择当反演迭代次数满足所规定的最大值时终止迭代进程,并输出反演模型。

2 基于交叉梯度耦合的二维联合反演方法

2.1 二维交叉梯度函数

交叉梯度函数是由Meju与Gallardo等[22]首次提出的:

t(x,y,z)=Ñmr(x,y,z)×Ñmv(x,y,z),

式中:mrmv表示2种不同的参与联合反演的物性参数,他们在空间上的分布是已知的,文中这2种物性参数分别为取对数的电阻率和取对数的速度;Ñ表示对模型参数取梯度。

对于二维模型,设地质体走向方向为沿y轴正方向,模型在y方向上的偏导数为0,那么txz方向上的分量均为0,只有y方向上有值,即

tx,z= mrx,zzmvx,zx- mr(x,z)xmv(x,z)z

使用中心差分格式将方程(5)进行如图1所示的离散化处理,mrmv共用一套相同的网格,其中的网格可以是等间距剖分也可以是不等间距剖分。

图1

图1   反演网格离散化示意

Fig.1   Inversion grid discrete schemes


用差商替代微商,得到交叉梯度函数的数值表达式:

ti,j= mr(i+1,j)-mr(i-1,j)Δzi+(Δzi+1+Δzi-1)/2· mv(i,j+1)-mv(i,j-1)Δxj+(Δxj+1+Δxj-1)/2- mr(i,j+1)-mr(i,j-1)Δxj+(Δxj+1+Δxj-1)/2· mv(i+1,j)-mv(i-1,j)Δzi+(Δzi+1+Δzi-1)/2,

式中:ij分别表示z方向和x方向的网格数目。

为验证上述交叉梯度算法的正确性,建立了如图2所示的2种模型Model 1和Model 2。其中,Model 1的模型参数从中心开始,以放射状向周围递增,Model 2的模型参数呈水平状分布,由浅至深逐增;2组模型具有相同的网格剖分。

图2

图2   两种模型参数的交叉梯度计算结果

Fig.2   The cross-gradient calculation results of two models


对这2个模型进行交叉梯度计算,结果如图2c所示:在中心位置的垂直方向上,即2组模型参数的梯度方向相互平行时,其交叉梯度值为零;当2种模型参数的梯度方向存在一定夹角时,其交叉梯度值不为零,且随着角度的增大而增大;当在深度为50 m处即梯度方向相互垂直时,交叉梯度值达到最大。

2.2 联合反演网格

全波形反演使用等间距网格进行剖分,边界条件采用的是PML吸收边界条件。海洋可控源电磁法的单方法反演,其异常区域内的网格是等间距的,垂直向下以及水平两端方向上的网格是不等间距的,且离异常体的距离越远间距越大,边界条件采用的是Dirichlet外边界条件。2种方法在进行交叉梯度计算时选取公共网格,对于其他部分,令交叉梯度的权重因子为0。

在进行MCSEM与地震全波形二维联合反演时,采用如图3所示的反演网格:红色线框内为计算交叉梯度的网格,就是地震全波形所剖分的网格与MCSEM的等间距网格的重合部分,即联合反演区域,剩余部分的网格为MCSEM反演扩充区域。

图3

图3   联合反演网格示意

Fig.3   The sketch of the Cross-Gradient mesh


2.3 联合反演理论

基于海洋可控源电磁法与地震全波形的联合反演,并结合两者反演的目标函数,构建如下基于交叉梯度约束的联合反演目标函数:

Φ(mr,mv)= Φmcsem+ λcsem-1Φdcsem+ λv-1Φdseis+γΦcross,

其中:

Φmcsem=(mr-m0)TCm-1(mr-m0),Φdcsem=[F(mr)-dobscsem]TCd-1[F(mr)-dobscsem],Φdseis=[P(mv)-dobsseis]Wd[P(mv)-dobsseis],Φcross=t(mr,mv)Tt(mr,mv)

式中:Φmcsem为MCSEM的模型光滑项;Φdcsem为MCSEM的数据拟合差项;Φdseis为地震的数据拟合差项;Φcross为交叉梯度项;λcsem-1λv-1以及γ分别为调节各项的权重因子;F(mr)与dobscsem分别表示MCSEM的正演响应和观测数据向量;F(mv)与dobsseis分别表示地震全波形的正演响应和观测数据向量;Cd-1为数据协方差矩阵的逆矩阵;m0表示初始模型的电阻率;P(v)和dobsseis分别为地震的正演响应及观测数据向量;t(mr,mv)表示交叉梯度函数;T表示向量的转置矩阵,†表示共轭转置矩阵。

由此可见,通过交叉梯度函数耦合两种物性,在单方法反演过程中受到另外一种方法反演结果的影响,增强了两种方法之间反演结果在结构上的一致性。

采用Gallardo和Meju[22-23]提出的泰勒级数展开法,对t(mr,mv)求取偏导数并取其一阶项,可近似表示为

t(mr,mv)≈t(mr,0,mv,0)+(BcsemBseis) mr-mr,0mv-mv,0,

其中:

Bcsem=t(mr,mv)mr,Bseis=t(mr,mv)mv

对于求解交叉梯度的问题,本文采用文献[24-25]中的方法,将式(9)进行整理变形,可得到:

ti,j= 00-mr(i+1,j)-mr(i-1,j)ΔSi,jmr(i,j+1)-mr(i,j-1)ΔSi,j0-mr(i,j+1)-mr(i,j-1)ΔSi,jmr(i+1,j)-mr(i-1,j)ΔSi,j00· (0mv(i,j-1)mv(i-1,j)mv(i+1,j)mv(i,j+1)0)T,
ti,j= 00-mv(i+1,j)-mv(i-1,j)ΔSi,jmv(i,j+1)-mv(i,j-1)ΔSi,j0-mv(i,j+1)-mv(i,j-1)ΔSi,jmv(i+1,j)-mv(i-1,j)ΔSi,j00· (0mr(i,j-1)mr(i-1,j)mr(i+1,j)mr(i,j+1)0)T ;
ΔSi,j=[Δzi+(Δzi+1zi-1)/2][Δxj+(Δxj+1xj-1)/2] 。

求取每个网格单元的ti,j并将对应位置相加,式(11)、(12)可简化为

t=Wr,crossmr=Wv,crossmv,

式中:Wr,crossWv,cross为交叉梯度的系数矩阵,Wr,cross与速度模型mv的参数信息相关,Wv,cross与电阻率模型mr的参数信息相关。在式(14)的基础上,将交叉梯度函数对模型向量求偏导,得到:

ΦcrossmrT=2Wr,crossTWr,crossmr,ΦcrossmvT=2Wv,crossTWv,crossmv

对目标函数式(7)求取关于mrT的偏导数,得

Φ(mr,mv)mrT=2 Cm-1(mr-m0)-2 λcsem-1JTCd-1· [dobscsem-F(mr)]+2γ Wr,crossTWr,crossmr,

Cm'=Cm-1+γWr,crossTWr,cross,可推导出电阻率模型的最终迭代公式:

mr,k+1-mr,0=Cm'JkTΒk+1,Βk+1=(λcsemCd+JkCmJkT)-1Χk,Χk=[dobscsem-F(mr,k)+Jk(mr,k-mr,0)]

对目标函数式(7)求取关于mvT的偏导数,得到关于梯度向量的表达式:

g=Real{FT[S-1]TWdδ d︿*}+2γvWv,crossTWv,crossmv

其关于mvT的近似Hessian矩阵的对角矩阵可表示为

Ha=diag{Real(JTWdJ*)+2γvWv,crossTWv,cross},

进而可得到速度模型的最终迭代公式:

mv,k+1-mv,0=mv,mv=-α(diagHa+εI)-1Rmg,F=-AvP

在式(17)、(20)的基础上完成基于交叉梯度约束的MCSEM与地震全波形的二维联合反演研究。

2.4 联合反演流程

联合反演流程如图4所示。首先,输入观测参数和初始模型:对参数进行初始化,MCSEM和地震全波形分别输入电阻率和速度的海底空间模型,然后进行一次地震频率域全波形反演,得到速度模型,利用该速度模型去计算交叉梯度,以此来约束MCSEM的反演;在完成一次模型的更新后,将更新了的电阻率模型加至地震全波形反演中,并计算交叉梯度达到更新速度模型的效果;两种反演方法交替迭代进行,相互约束,直至反演迭代满足本文中所设定的终止要求(即最大迭代次数)。

图4

图4   联合反演流程

Fig.4   Joint-inversion flow


3 算例分析

3.1 单一板状体模型算例

建立如图5所示的高阻低速板状体理论模型,模拟海底地层中的油气储层构造。

图5

图5   单一板状体理论模型示意

Fig.5   Synthetic model of the single tabular body


设模型海水层深500 m,电阻率为0.3 Ω·m,速度为1 500 m/s;在海底之下800 m处有一个厚200 m、横向宽度2 000 m、电阻率100 Ω·m、速度1 800 m/s的储层;背景电阻率为1 Ω·m,速度为4 500 m/s。接收器布置在海底,沿测线方向每隔100 m放置一个,共计41个(-2~2 km);水平电偶源放置在距离海底50 m的高度上,长度为1 m,电流为1 A,沿测线方向设置21个发射源,间隔200 m(-2~2 km),发射频率分别为0.25、0.5、1.5 Hz。地震全波形使用海底放炮、海底接收的观测方式,检波器同样设置在海底沿y轴-2~2 km范围内,道间距100 m,共计41个;炮点沿海底面测线方向(-2~2 km) 设置41个,炮间距100 m。

理论模型正演数据(Ey的实部和虚部)加入2%的噪声作为观测数据进行反演;取1、2、3、5、6、8和10 Hz共7个频率进行反演。地震全波形使用51×81的均匀网格进行剖分,大小为100 m×100 m;MCSEM使用的网格大小为49×69,为不等间距网格,与地震全波形网格剖分的异常区域是相同的。MCSEM和地震全波形的反演初始模型均为空气—海水—岩层模型,反演设置在25次迭代后停止。MCSEM和地震全波形的反演结果见图6

图6

图6   单一板状体模型反演结果及拟合差迭代曲线

Fig.6   The profile of the inversion as well as the iterative curve of the RMS misfit of the single tabular body model


图中可见,MCSEM单独反演的结果只能恢复异常体的基本位置,形状呈倒三角状,底部收敛性较差,电阻率为25~60 Ω·m,周围存在部分的冗余构造(图6a、b);而联合反演的结果明显改善了异常体的位置信息以及模型的整体形态特征(图6c、d),其电阻率异常值为60~85 Ω·m,对真值的恢复能力有明显的提高,压制了周围的虚假异常。地震全波形单方法反演对异常体的速度和边界的恢复效果较好,但速度联合反演结果相比于单方法没有明显的改进。

图6e、f给出了MCSEM、地震全波形反演的数据拟合差收敛曲线。反演过程中,两种方法的RMS值均稳定下降,在后期趋于稳定,其中MCSEM单独反演的RMS值降到1.88,联合反演的RMS值降到约1.42;地震全波形单独反演与联合反演的RMS曲线基本重合,迭代22次后RMS值稳定在1.35。

3.2 复杂板状体模型算例

建立如图7所示的复杂板状体理论模型,模拟海底地层中存在上覆异常体的油气储层。

图7

图7   复杂板状体理论模型示意

Fig.7   Synthetic model of the complex tabular body


海水层深500 m,电阻率为0.3 Ω·m,速度为1 500 m/s;海底300 m处有一个厚200 m、横向宽度2 000 m、电阻率10 Ω·m、速度5 000 m/s的低阻高速板状体,在其下方500 m处有一个厚200 m、横向宽度2 000 m、电阻率100 Ω·m、速度1 800 m/s的高阻低速板状体;背景电阻率为1 Ω·m,速度为3 200 m/s。地震与MCSEM的观测系统、反演网格和迭代次数的设置与3.1节相同。

MCSEM和地震的反演结果见图8。MCSEM单独反演结果中的异常体整体形态与理论模型偏差较大,收敛性很差,上层低阻板状体没有显示出来,下层高阻板状体反演后的位置在1 200 m处(图8a),这主要是由于上层低阻板状体对下层高阻板状体的影响使得高阻板状体的位置与原始模型有较大的偏差,电阻率值在15~60 Ω·m,周围存在较大的冗余构造。联合反演的结果明显改善了异常体的位置信息以及模型的形态,虽然对于上层的低阻板状体反演效果还是较差,仅仅是电阻率值稍有改善,但是对于下层高阻板状体的位置,形状都基本恢复了(图8c),其电阻率异常值为50~80 Ω·m,对真值的恢复能力有了明显的提高,并且明显压制了周围的虚假异常。而地震的单方法反演对异常体的速度和边界的恢复效果较好,只是异常体周围存在少许的假异常(图8b),联合反演速度结果相比于单方法没有明显的提高(图8d)。

图8

图8   复杂板状体模型反演结果及拟合差迭代曲线

Fig.8   The profile of the inversion as well as the iterative curve of the RMS misfit of the complex tabular body model


图8e、f给出了MCSEM与地震全波形反演的数据拟合差收敛曲线,可以发现在反演过程中,MCSEM的 RMS值前期下降速度较快,在后期趋于稳定,其单独反演的RMS值降到2.85,联合反演的RMS值下降到1.78;地震全波形单独反演与联合反演的RMS曲线基本重合,迭代24次后RMS值稳定在1.63。

3.3 凹槽体模型算例

建立如图9所示的低阻低速凹槽体理论模型,模拟海底地层中的背斜构造。

图9

图9   凹槽体理论模型示意

Fig.9   Synthetic model of the groove body


海水层深为500 m,电阻率为0.3 Ω·m,速度为1 500 m/s;在海底300 m处有一个厚1 000 m、横向宽1 500 m、电阻率10 Ω·m、速度1 800 m/s的凹槽体;地震与MCSEM采用与3.1节相同的观测系统、反演网格和迭代次数设置,反演结果见图10

图10

图10   凹槽体模型反演结果及拟合差迭代曲线

Fig.10   The profile of the inversion as well as the iterative curve of the RMS misfit of the groove body model


MCSEM单独反演的结果只能恢复异常体的基本位置,形状与理论模型偏差较大,没有刻画出异常体的边界,整体的收敛性较差(图10a),电阻率为15~60 Ω·m;联合反演的结果明显改善了异常体的位置信息以及模型的形态,基本恢复到了理论模型的位置与形状(图10c),其电阻率异常值为10~30 Ω·m,对真值的恢复能力有了明显的提高。全波形单方法反演对异常体的速度和边界的恢复效果较好,联合反演速度结果相比于单方法没有明显的改进(图10b、d)。

图10e、f显示在反演过程中两者的RMS值均稳定下降,在后期趋于稳定,其中MCSEM单独反演的RMS值降到1.92,联合反演的RMS值降到1.41左右;地震全波形单独反演与联合反演的RMS曲线基本重合,迭代23次后RMS值稳定在1.51。

4 结论与讨论

1)对MCSEM和地震全波形的二维联合反演理论公式进行了推导,推导出含有交叉梯度项的MCSEM的电阻率迭代公式以及地震全波形的速度迭代公式;建立了基于交叉梯度约束的MCSEM与地震全波形反演的目标函数,开发了一套二维联合反演算法,并验证了该算法的可行性与适用性。

2)建立多组不同的理论模型,进行联合反演计算,并与单一方法反演的结果进行了对比分析。电阻率模型的联合反演结果有显著改善和提升,而速度模型的联合反演结果基本与单方法的结果相同,说明全波形的反演方法能够提高MCSEM反演结果的可靠性。

本文所研究的MCSEM与地震全波形二维联合反演算法并没有加入实际数据的应用,这就使得对于算法精确性的验证具有一定的局限性,下一步需结合实测资料进行正、反演的计算。

参考文献

Haber E, Oldenburg D.

Joint inversion:A structural approach

[J]. Inverse Problems, 1997, 13(1):63-77.

DOI:10.1088/0266-5611/13/1/006      URL     [本文引用: 1]

Pak Y C, Li T L, Kim G S.

2D data-space cross-gradient joint inversion of MT,gravity and magnetic data

[J]. Journal of Applied Geophysics, 2017, 143:212-222.

DOI:10.1016/j.jappgeo.2017.05.013      URL     [本文引用: 1]

Coutant O, Bernard M L, Beauducel F, et al.

Joint inversion of P-wave velocity and density,application to La Soufrière of Guadeloupe hydrothermal system

[J]. Geophysical Journal International, 2012, 191(2):723-742.

DOI:10.1111/gji.2012.191.issue-2      URL     [本文引用: 1]

敬荣中, 鲍光淑, 陈绍裘.

地球物理联合反演研究综述

[J]. 地球物理学进展, 2003, 18(3):535-540.

[本文引用: 1]

Jing R Z, Bao G S, Chen S Q.

A review of the researches for geophysical combinative inversion

[J]. Progress in Geophysics, 2003, 18(3):535-540.

[本文引用: 1]

Hu W Y, Abubakar A, Habashy T M.

Joint electromagnetic and seismic inversion using structural constraints

[J]. Geophysics, 2009, 74(6):R99-R109.

DOI:10.1190/1.3246586      URL     [本文引用: 1]

We have developed a frequency-domain joint electromagnetic (EM) and seismic inversion algorithm for reservoir evaluation and exploration applications. EM and seismic data are jointly inverted using a cross-gradient constraint that enforces structural similarity between the conductivity image and the compressional wave (P-wave) velocity image. The inversion algorithm is based on a Gauss-Newton optimization approach. Because of the ill-posed nature of the inverse problem, regularization is used to constrain the solution. The multiplicative regularization technique selects the regularization parameters automatically, improving the robustness of the algorithm. A multifrequency data-weighting scheme prevents the high-frequency data from dominating the inversion process. When the joint-inversion algorithm is applied in integrating marine controlled-source electromagnetic data with surface seismic data for subsea reservoir exploration applications and in integrating crosswell EM and sonic data for reservoir monitoring and evaluation applications, results improve significantly over those obtained from separate EM or seismic inversions.

Brown V, Singh S C, Key K.

Using seismic full waveform inversion to constrain controlled-source electromagnetic inversion

[J]. SEG Expanded Abstracts, 2010, 29:619-623.

[本文引用: 1]

da Silva N V, Morgan J, Warner M, et al.

3D constrained inversion of CSEM data with acoustic velocity using full waveform inversion

[C]// SEG Technical Program Expanded Abstracts 2012.Society of Exploration Geophysicists, 2012.

[本文引用: 1]

Mac-Gregor L, Bouchrara S, Tomlinson J, et al.

Integrated analysis of CSEM,seismic and well log data for prospect appraisal:A case study from West Africa

[J]. First Break, 2012, 30(4):43-49.

[本文引用: 1]

Saunders J H, Herwanger J V, Pain C C, et al.

Constrained resistivity inversion using seismic data

[J]. Geophysical Journal International, 2005, 160(3):785-796.

DOI:10.1111/gji.2005.160.issue-3      URL     [本文引用: 1]

杜润林. 海洋可控源电磁场和地震波场联合反演方法研究[D]. 东营: 中国石油大学(华东), 2015.

[本文引用: 1]

Du R L. The research of joint inversion of marine controlled source electromagnetic and seismic[D]. Dongying: China University of Petroleum (Huadong), 2015.

[本文引用: 1]

李刚. 海洋可控源电磁与地震资料构造联合反演方法研究[D]. 青岛: 中国海洋大学, 2015.

[本文引用: 1]

Li G. Jont inversion of marine controlled-source electromagnetic and seismic data using the structural constraints[D]. Qingdao: Ocean University of China, 2015.

[本文引用: 1]

Hustedt B, Operto S, Virieux J.

Mixed-grid and staggered-grid finite-difference methods for frequency-domain acoustic wave modelling

[J]. Geophysical Journal International, 2004, 157(3):1269-1296.

DOI:10.1111/gji.2004.157.issue-3      URL     [本文引用: 1]

Jo C H, Shin C, Suh J H.

An optimal 9-point,finite-difference,frequency-space,2-D scalar wave extrapolator

[J]. Geophysics, 1996, 61(2):529-537.

DOI:10.1190/1.1443979      URL     [本文引用: 1]

In this study, a new finite‐difference technique is designed to reduce the number of grid points needed in frequency‐space domain modeling. The new algorithm uses optimal nine‐point operators for the approximation of the Laplacian and the mass acceleration terms. The coefficients can be found by using the steepest descent method so that the best normalized phase curves can be obtained. This method reduces the number of grid points per wavelength to 4 or less, with consequent reductions of computer memory and CPU time that are factors of tens less than those involved in the conventional second‐order approximation formula when a band type solver is used on a scalar machine.

Berenger J P.

A perfectly matched layer for the absorption of electromagnetic waves

[J]. Journal of Computational Physics, 1994, 114(2):185-200.

DOI:10.1006/jcph.1994.1159      URL     [本文引用: 1]

Virieux J, Operto S.

An overview of full-waveform inversion in exploration geophysics

[J]. Geophysics, 2009, 74(6):WCC1-WCC26.

DOI:10.1190/1.3238367      URL     [本文引用: 1]

Full-waveform inversion (FWI) is a challenging data-fitting procedure based on full-wavefield modeling to extract quantitative information from seismograms. High-resolution imaging at half the propagated wavelength is expected. Recent advances in high-performance computing and multifold/multicomponent wide-aperture and wide-azimuth acquisitions make 3D acoustic FWI feasible today. Key ingredients of FWI are an efficient forward-modeling engine and a local differential approach, in which the gradient and the Hessian operators are efficiently estimated. Local optimization does not, however, prevent convergence of the misfit function toward local minima because of the limited accuracy of the starting model, the lack of low frequencies, the presence of noise, and the approximate modeling of thewave-physics complexity. Different hierarchical multiscale strategies are designed to mitigate the nonlinearity and ill-posedness of FWI by incorporating progressively shorter wavelengths in the parameter space. Synthetic and real-data case studies address reconstructing various parameters, from [Formula: see text] and [Formula: see text] velocities to density, anisotropy, and attenuation. This review attempts to illuminate the state of the art of FWI. Crucial jumps, however, remain necessary to make it as popular as migration techniques. The challenges can be categorized as (1) building accurate starting models with automatic procedures and/or recording low frequencies, (2) defining new minimization criteria to mitigate the sensitivity of FWI to amplitude errors and increasing the robustness of FWI when multiple parameter classes are estimated, and (3) improving computational efficiency by data-compression techniques to make 3D elastic FWI feasible.

周丽芬. 大地电磁与地震数据二维联合反演研究[D]. 北京: 中国地质大学(北京), 2012.

[本文引用: 1]

Zhou L F. Two dimensional joint inversion of MT and seismic data[D]. Beijing: China University of Geosciences, 2012.

[本文引用: 1]

Pratt R G.

Seismic waveform inversion in the frequency domain,Part 1:Theory and verification in a physical scale model

[J]. Geophysics, 1999, 64(3):888-901.

DOI:10.1190/1.1444597      URL     [本文引用: 1]

Seismic waveforms contain much information that is ignored under standard processing schemes; seismic waveform inversion seeks to use the full information content of the recorded wavefield. In this paper I present, apply, and evaluate a frequency‐space domain approach to waveform inversion. The method is a local descent algorithm that proceeds from a starting model to refine the model in order to reduce the waveform misfit between observed and model data. The model data are computed using a full‐wave equation, viscoacoustic, frequency‐domain, finite‐difference method. Ray asymptotics are avoided, and higher‐order effects such as diffractions and multiple scattering are accounted for automatically. The theory of frequency‐domain waveform/wavefield inversion can be expressed compactly using a matrix formalism that uses finite‐difference/finite‐element frequency‐domain modeling equations. Expressions for fast, local descent inversion using back‐propagation techniques then follow naturally. Implementation of these methods depends on efficient frequency‐domain forward‐modeling solutions; these are provided by recent developments in numerical forward modeling. The inversion approach resembles prestack, reverse‐time migration but differs in that the problem is formulated in terms of velocity (not reflectivity), and the method is fully iterative. I illustrate the practical application of the frequency‐domain waveform inversion approach using tomographic seismic data from a physical scale model. This allows a full evaluation and verification of the method; results with field data are presented in an accompanying paper. Several critical processes contribute to the success of the method: the estimation of a source signature, the matching of amplitudes between real and synthetic data, the selection of a time window, and the selection of suitable sequence of frequencies in the inversion. An initial model for the inversion of the scale model data is provided using standard traveltime tomographic methods, which provide a robust but low‐resolution image. Twenty‐five iterations of wavefield inversion are applied, using five discrete frequencies at each iteration, moving from low to high frequencies. The final results exhibit the features of the true model at subwavelength scale and account for many of the details of the observed arrivals in the data.

Ward S H, Hohmann G W.

Electromagnetic theory for geophysical applications

[M]//Electromagnetic Methods in Applied Geophysics. Society of Exploration Geophysicists, 1988:130-311.

[本文引用: 1]

刘颖. 海洋可控源电磁法二维有限元正演及反演[D]. 青岛: 中国海洋大学, 2014.

[本文引用: 1]

Liu Y. 2D finite element modeling and inversion for marine controlled-source electromagnetic fields[D]. Qingdao: Ocean University of China, 2014.

[本文引用: 1]

Mitsuhata Y.

2-D electromagnetic modeling by finite-element method with a dipole source and topography

[J]. Geophysics, 2000, 65(2):465.

DOI:10.1190/1.1444740      URL     [本文引用: 2]

I present a method for calculating frequency‐domain electromagnetic responses caused by a dipole source over a 2-D structure. In modeling controlled‐source electromagnetic data, it is usual to separate the electromagnetic field into a primary (background) and a secondary (scattered) field to avoid a source singularity, and only the secondary field caused by anomalous bodies is computed numerically. However, this conventional scheme is not effective for complex structures lacking a simple background structure. The present modeling method uses a pseudo‐delta function to distribute the dipole source current, and does not need the separation of the primary and the secondary field. In addition, the method employs an isoparametric finite‐element technique to represent realistic topography. Numerical experiments are used to validate the code. Finally, a simulation of a source overprint effect and the response of topography for the long‐offset transient electromagnetic and the controlled‐source magnetotelluric measurements is presented.

Siropunvarapor W, Egbert G.

An efficient data-subspace inversion method for 2-D magnetotelluric data

[J]. Geophysics, 2000, 65(3):791-803.

DOI:10.1190/1.1444778      URL     [本文引用: 1]

There are currently three types of algorithms in use for regularized 2-D inversion of magnetotelluric (MT) data. All seek to minimize some functional which penalizes data misfit and model structure. With the most straight‐forward approach (exemplified by OCCAM), the minimization is accomplished using some variant on a linearized Gauss‐Newton approach. A second approach is to use a descent method [e.g., nonlinear conjugate gradients (NLCG)] to avoid the expense of constructing large matrices (e.g., the sensitivity matrix). Finally, approximate methods [e.g., rapid relaxation inversion (RRI)] have been developed which use cheaply computed approximations to the sensitivity matrix to search for a minimum of the penalty functional. Approximate approaches can be very fast, but in practice often fail to converge without significant expert user intervention. On the other hand, the more straightforward methods can be prohibitively expensive to use for even moderate‐size data sets. Here, we present a new and much more efficient variant on the OCCAM scheme. By expressing the solution as a linear combination of rows of the sensitivity matrix smoothed by the model covariance (the “representers”), we transform the linearized inverse problem from the M-dimensional model space to the N-dimensional data space. This method is referred to as DASOCC, the data space OCCAM’s inversion. Since generally N ≪ M, this transformation by itself can result in significant computational saving. More importantly the data space formulation suggests a simple approximate method for constructing the inverse solution. Since MT data are smooth and “redundant,” a subset of the representers is typically sufficient to form the model without significant loss of detail. Computations required for constructing sensitivities and the size of matrices to be inverted can be significantly reduced by this approximation. We refer to this inversion as REBOCC, the reduced basis OCCAM’s inversion. Numerical experiments on synthetic and real data sets with REBOCC, DASOCC, NLCG, RRI, and OCCAM show that REBOCC is faster than both DASOCC and NLCG, which are comparable in speed. All of these methods are significantly faster than OCCAM, but are not competitive with RRI. However, even with a simple synthetic data set, we could not always get RRI to converge to a reasonable solution. The basic idea behind REBOCC should be more broadly applicable, in particular to 3-D MT inversion.

Meju M A, Gallardo L A, Mohamed A K.

Evidence for correlation of electrical resistivity and seismic velocity in heterogeneous near-surface materials

[J]. Geophysical Research Letters, 2003, 30(7):1373-1376.

[本文引用: 2]

Gallardo L A, Meju M A.

Joint two-dimensional DC resistivity and seismic travel time inversion with cross-gradients constraints

[J]. Journal of Geophysical Research:Solid Earth, 2004, 109(B3):B03311.

[本文引用: 1]

王涛. 磁法与可控源音频大地电磁法二维联合反演研究[D]. 北京: 中国地质大学(北京), 2016.

[本文引用: 1]

Wang T. Two-dimensional joint inversion of magnetic and CSAMT methods[D]. Beijing: China University of Geosciences(Beijing), 2016.

[本文引用: 1]

郭长安. 电阻率法和地震全波形二维联合反演研究[D]. 北京: 中国地质大学(北京), 2018.

[本文引用: 1]

Guo C A. Two-dimensional joint inversion of resistivity and seismic full-waveform method[D]. Beijing: China University of Geosciences(Beijing), 2018.

[本文引用: 1]

/

京ICP备05055290号-3
版权所有 © 2021《物探与化探》编辑部
通讯地址:北京市学院路29号航遥中心 邮编:100083
电话:010-62060192;62060193 E-mail:whtbjb@sina.com , whtbjb@163.com