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

孔繁祥, 谭捍东, 刘建勋

(1. 中国地质科学院 地球物理地球化学勘查研究所,河北 廊坊 065000;2. 中国地质大学(北京) 地球物理与信息技术学院,北京 100083)

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]中关于一次场计算公式的详细推导过程,推导了任意方向的电偶源作用下水平层状介质中各个电磁场分量的计算公式。

(1)

(2)

(3)

式中: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),

(4)

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

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

(5)

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

图1 反演网格离散化示意Fig.1 Inversion grid discrete schemes

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

(6)

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

为验证上述交叉梯度算法的正确性,建立了如图2所示的2种模型Model 1和Model 2。其中,Model 1的模型参数从中心开始,以放射状向周围递增,Model 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 联合反演网格示意Fig.3 The sketch of the Cross-Gradient mesh

2.3 联合反演理论

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

(7)

其中:

(8)

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

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

(9)

其中:

(10)

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

(11)

(12)

ΔSi,j=[Δzi+(Δzi+1+Δzi-1)/2][Δxj+(Δxj+1+Δxj-1)/2]。

(13)

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

t=Wr,crossmr=Wv,crossmv,

(14)

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

(15)

(16)

(17)

(18)

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

(19)

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

(20)

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

2.4 联合反演流程

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

图4 联合反演流程Fig.4 Joint-inversion flow

3 算例分析

3.1 单一板状体模型算例

建立如图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 单一板状体模型反演结果及拟合差迭代曲线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 复杂板状体理论模型示意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 复杂板状体模型反演结果及拟合差迭代曲线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 凹槽体理论模型示意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。

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