第46卷 第2期2024年3月
物探化探计算技术
C O M P U T I N G T E C H N I Q U E S F O R G E O P H Y S I C A L A N
D G
E O C H E M I C A L E X P L O R A T I O N
V o l .46 N o .2
M a r .2024
文章编号:1001-1749(2024)02-0197-09起伏地形的2.5维主轴各向异性海洋C S E M 反演研究
覃金生
(东华理工大学地球物理与测控技术学院，南昌 330013
)摘 要：笔者利用O C C A M 方法实现了M P I 并行的2.5维海洋C S E M 反演算法㊂正演采用基于非结
构网格的有限元方法，在精确模拟海底起伏地形和构造模型的同时，对核心区域和扩展区域进行不同尺度的剖分以减少反演参数㊂水平海床和起伏海床的理论模型算例表明，沿测线发射的2.5维剖面数据对垂直各向异性模型的恢复效果优于水平各向异性模型㊂复杂海底地形模型的试算结果表明，带地形沉积层的存在会对反演结果产生较大影响，在实际资料解释中应该引起重视㊂
关键词：海洋可控源电磁；起伏地形；反演；各向异性
中图分类号:P 631  文献标志码:A  D O I :10.3969/j
.i s s n .1001-1749.2024.02.09收稿日期:2023-08-11
作者简介：覃金生(1999-)，男，硕士，主要研究方向为地球物理电磁勘探,E -m a i l :2116205112@q q
.c o m ㊂0 引言
近二十年是海洋勘探技术发展的黄金时期，作为海洋地震勘探的有效辅助手段，海洋可控源电磁方法(M a r i n e C o n t r o l l e d S o u r c e E l e c t r o m a g
n e t i c M e t h o d ,M C S E M )
被广泛应用于海底油气资源与深部地质构造的探测[1-3
]㊂我国对M C S E M 的研究
及其应用起步较晚，数值模拟方面仍以二维和三维各向同性为主㊂但由于海洋油气成藏环境的复杂性，使得储层在宏观上往往呈现电导率各向异性特
征[4-5
]，此时，在二维地质构设条件下，各向同
性的2.5维反演算法并不能适应这种模型，往往难以准确反映储层特征㊂因此，从物理模型的适应性出发，在现阶段2.5维反演技术作为主流的条件下，对反演算法进行拓展以支持各向异性介质反演是极其必要的，可综合对比同性和异性反演结果，有效提升解释的可靠性，并有望通过各向异性参数的反演恢复对储层电性解释提供更多参考㊂
近年来，海洋可控源电磁各向异性介质影响的
研究受到了广大学者的关注㊂罗鸣等[6-7
]详细分析
了海洋覆盖层和高阻储层电阻率各向异性对电磁场响应的影响，并在此基础上实现了一维垂直各向异性的频率域海洋可控源数据的反演;M a s n a g
h e t t i 等[8]对2.5维海洋可控源电磁垂直各向异性介质反
演进行分析研究，结果表明必须将o f f l i n e 数据考虑
在内可以充分恢复水平电阻率分量;K e y 等[9]和
R a m a n a n j
a o n a 等[10
]实现了2.5维垂直各向异性介质条件下的M C S E M 反演研究;W i i k 等[11
]年采用基于电磁场的积分方程的对比源方法实现了垂直各向异性介质的M C S E M 三维反演,
指出垂直各向异性介质反演对弱异常的识别更加有效;N e w m a n 等[12]㊁彭荣华等[13]㊁赵宁等[14]和M o r t e n 等[1
5
]采用不同的正演求解方法分别实现了V T I 介质的M C -S E M 三维反演㊂在各向异性条件下,
海床地形和沉积层往往会对反演结果产生影响，而该方面的研究工作还相对较少，需要进一步的探索㊂
笔者对起伏地形2.5维海洋可控源电磁各向异
性问题进行反演研究，为精确模拟海底起伏地形和地质构造模型，采用非结构三角网格进行剖分，同时
对核心区域和扩展区域分别采用不同网格尺寸进行剖分，确保计算精度的条件下减小计算工作量㊂采
用基于M P I 并行的O c c a m 反演方法[1
6-17
]，对不同情况下的各向异性模型分别进行试算㊂
1 主轴各向异性正演理论
海洋电磁勘探利用的频率范围通常为0.01H z
~10H z ，满足准静态条件(σ≫ωε)
，
因此求解电磁场时可忽略位移电流的影响㊂此时可建立以下频率
域电磁场控制方程(取时间因子e
-i ωt
)：浣ˑE -i ωμ
H =0(1
)浣ˑH -σE =J s
(2
)其中,E 为电场强度;H 为磁场强度;J s 为电性激励源的电流密度矢量;σ为模型电导率;μ取真空磁导
率μ0(
不考虑介质的磁导率变化)㊂在主轴各向异性反演研究中，介质的电导率只有三个沿坐标轴方向的分量，其张量形式可表示为:
σ=σx
0 00 σy
00 0 σz
(3
)假设地质体走向为x 方向,σ张量只沿坐标轴
y ㊁
z 方向变化，沿走向不变，此时由于场源为三维响应，所以海洋可控源电磁主轴各向异性问题是一个
2.5维问题㊂考虑到海洋可控源电磁勘探多场源的特性以及复杂地形的影响，采用直接求解总场的方
法，在直角坐标系中，对式(1)和式(2)沿走向进行傅里叶变换并展开得到如下方程:
∂E ^z ∂y -∂E ^y
∂
z =i ωμ0H ^x (4)∂E ^x
∂z
-i k x E ^z =i ωμ0H ^y (5)i k x E ^y -∂E ^x
∂y
=i ωμ
0H ^z (6)∂H ^z
∂y
-
∂H ^y
∂z
-σx
E ^x =J ^x (7)∂H ^x
∂z
-i k x H ^z -σy
E ^y =J ^y (8)i k x H ^y -∂H ^x
∂y
-σz E ^z =J ^z
(9
)其中,k x 是x 方向上的空间波数,J ^x ㊁J ^y 和J ^z 代表
波数域中场源电流密度沿直角坐标系的坐标轴分
量㊂
结合方程(4)~(9
)，可得到以下微分方程:∂∂y (σy λy ∂E ^x ∂y )+∂∂z (σz λz ∂E ^x ∂z )-∂∂y (i k x λy ∂H ^x ∂
z )+∂∂z (i k x λz ∂H ^x ∂y )-σx E ^x =-∂∂y (i k x λy J ^y )-∂∂z (i k x λz J ^z )
+J ^x
(10
)∂∂y (C λz ∂H ^x ∂y )+∂∂z (C λy
∂H ^x ∂z )-∂∂y (i k x λz ∂E ^x ∂z )+
∂∂z (i k x λy ∂E ^x ∂y )-C H ^x =-∂∂y (C λz J ^z )+
∂∂z (C λy
J ^y )(11
)其中,C =i ωμ0,λy =k 2x -i ωμ0σy ,λz =k 2
x -i ωμ
0σz ㊂为尽量准确模拟海床起伏，在正演计算区域Ω
内使用非结构有限单元法数值求解以上控制方程的边值问题，式(10)和(11
)的加权余量积分如下:ʏΩA 1㊃δE ^x d Ω=ʏΩ
f 1㊃δE ^x d Ω(12
)ʏΩA 2㊃δH ^x d Ω=ʏΩ
f 2㊃δH ^x d Ω(13)其中,A 1㊁A 2分别为方程(10)㊁(11)的左端项,f 1㊁
f 2分别为方程(10)㊁(11)的右端场源项㊂假设n y
㊁n z 表示单元边界外法线与坐标轴之间夹角的三角函数，利用格林公式[18
]，结合方程(12)和(13
)，方程(10)㊁(11
)变为:-ʏΩ(
σy λy ∂
E ^x ∂y ∂δE ^x ∂y +σz λz ∂E ^x ∂z ∂δE ^x ∂z -i k x λy ∂H ^x ∂
z ∂δE ^x ∂y +i k x λz ∂H ^x ∂y ∂δE ^x ∂z )d Ω+ɥΓ(σy λy ∂E ^x ∂
y δE ^x n y +σz λz ∂E ^x ∂z δE ^x n z -i k x λy
∂H ^x ∂z δE ^x n y +i k x λz ∂
H ^x ∂y δE ^x n z )d Γ-ʏΩσx E ^x δE ^x d Ω=ʏΩ(i k x λy
J ^y ∂δE ^x ∂y
+i k x λz
J ^z
∂δE ^x ∂z
+J ^x δE ^x )d Ω-ɥΓ(i k x λy
J ^y δE ^x n y +
i k x λz
J ^z δE ^x n z )d Γ(14
)-ʏΩ(
C λz ∂H ^x ∂y ∂δH ^x ∂y +C λY ∂H ^x ∂z ∂δH ^x ∂z -i k x λz ∂E ^x ∂
z ∂δH ^x ∂y +i k x λy ∂E ^x ∂y ∂δH ^x ∂z )d Ω+ɥΓ(C λz ∂H ^x ∂y δH ^x n y +C λy
891    物探化探计算技术
46卷
∂H ^x
∂z δH ^x n z -i k x λz ∂E ^x ∂z δH ^x n y +i k x λy ∂
E ^x ∂
y δH ^x n z )d Γ-ʏΩC H ^x δH ^x d Ω=ʏΩ(C λz
J ^z
∂δH ^x ∂y
-C λy
J ^y ∂δH ^x ∂z
)d Ω-ɥΓ(C λz J ^z δH ^x n y -C λy
J ^y
δH ^x n z )d Γ(15)采用开源代码T r i a n g l e 程序[19
]构建模型的非结构有限单元网格，将区域Ω离散为多个三角单
元，对每个单元的式(14)㊁(15
)分别积分并联立为线性方程组㊂采用P a r d i s o 直接求解器求解该方程组即可得到波数域的电磁场解E ^x 和H ^x ，再通过傅里叶反变换即可得到空间域电磁场值㊂
2 O c c a m 反演
O c c a m 反演算法是基于正则化思想的最小二
乘方法[20-22
]，它对以下无约束最优化问题求极小:
U =μ
R m  2+ W (d -F (m )) 2(16)其中,m 是n 维模型参数向量;R 为粗糙度算子矩阵;μ是正则化算子，用于平衡模型粗糙度与数据拟合差，若μ取值较大，反演结果倾向于结果的光滑，反之则更倾向于拟合数据;W 是对拟合差的对角加权矩阵;d 为观测数据矢量;F (m )为模型m 对应的正演响应㊂
对于这里所研究的主轴各向异性模型，模型向量可表示为:
m =m x m y
m z
(17
)其中,m x =(σx 1,σx 2, ,σx i )
，表示第i 个剖分单元内电导率的x 分量㊂则模型粗糙度可表示为:
R m  2= R m x  2+ R m y  2
z  2'2(
)其中，标量λ为各向异性惩罚项,m '=m y m z m x  T
㊂对于第k 次迭代模型m k ，使目标函数(16
)充分下降的迭代式取以下形式:
m k +1=m k +μR T R +(W J k )T W J k  -1
ˑ(W J k )T W d ^-μR T
R
m k  (19
)其中，拟合残差向量d ^=d -F (m k )
㊂电磁反演工作会消耗大量的计算成本，尤其是在各向异性情况下，二维海洋C S E M 各向异性反演
是一个非常耗时㊁耗内存的计算过程，需要根据计算
平台应用合理的并行加速技术[
23
]㊂由于海洋可控源电磁正演各频率的方程相互独立，具有良好的并行性质，因此目前的并行加速方案多采用O p e n M P 或M P I 按频率并行的方案[24
]㊂其中,M P I 协议
(M e s s a g e P a s s i n g I
n t e r f a c e P a r a l l e l P r o t o c o l )为进程级并行编程协议，每个进程都有独立的存储空间,
可适应于集环境[25-26
]，能够实现大规模反演问题的加速计算㊂本文将不同场源㊁不同频率和不同接
收位置的计算任务进行合理分配，在个人计算机上实现了M P I 并行加速计算㊂
3 模型算例
3.1 海床水平地形模型
如图1(a )所示，模型y 方向水平向右,z 方向垂
直向下，顶部空气层电阻率设定为1ˑ1012Ω㊃m ,
中间海水层深度为1k m ，电阻率为0.3Ω㊃m ,
海底沉积层中含有一个4k mˑ1k m 的高阻油气储层，其顶部埋深2.5k m ㊂通过水平电偶源沿测线进行
激发，发射频率采用0.1H z 和1.0H z
，
在海底上方50m 处-6k m<y <6k m 范围内设置了11个发射场源，在海底上方0.1m 处-8k m<y <8k m 范围内布设了个接收器㊂
正则化反演图1 水平地形模型
F i g .1 H o r i z o n t a l t o p o g r a p h y m
o d e l 9
912期覃金生：起伏地形的2.5维主轴各向异性海洋C S E M 反演研究
该模型储层取两种各向异性参数:①储层为
V T I 介质,
水平方向和垂直方向的电阻率分别为20Ω㊃m 和100Ω㊃m ;②储层为H T I 介质,x 方向和y 方向的电阻率分别为20Ω㊃m 和100Ω㊃m ，沉积层电阻率均为1.0Ω㊃m ㊂
对于沉积层，同样取两种各向异性参数:①沉积层为V T I 介质，水平方向和垂直方向的电阻率分别为1.0Ω㊃m 和5.0Ω㊃m ;②沉积层为H T I 介质,x 方向和y 方向的电阻率分别为1.0Ω㊃m 和
5.0Ω㊃m ，储层电阻率均为100Ω㊃m ㊂
正演计算得到的模型电磁场分量的幅值与相位，在添加4%的随机高斯噪声之后,
作为反演计算所需的观测数据㊂反演初始模型如图1(b
)所示，仅对-6.5k m ੭y ੭
6.5k m ,0k m ੭z ੭6.5k m 的目标区域进行精细网格剖分，而远离求解目标的扩展区域则选择粗网格进行剖分，且空气层和海水层作为固定参数不参与反演计算，共生成4927个三角单元和2556个节点，反演初始模型电阻率设定为1.0Ω㊃m ，并将反演过程的电阻率值限制在
㊃㊃范围内以获得更为准确的结果㊂
上述模型的反演结果如图2所示，图中的白圆点和三角为布设的发射器和接收器㊂由图2可见，储层为各向异性介质时，反演结果清晰恢复了储层的电阻率㊁位置和形态特征，且垂直各向异性模型的恢复效果比水平各向异性模型的恢复效果要好，更接近真实模型㊂沉积层为垂直各向异性介质时，储层会变薄，并向水平方向延伸，同时在
其周围出现
伪影
㊂沉积层为水平各向异性介质时，储层会向垂直方向拉伸，且在下方出现了虚假低阻异常㊂
图3为海床水平地形模型的反演均方根拟合差(R M S )
和粗糙度(R o u g h n e s s )随迭代过程的变化曲线㊂当迭代进行到第10次时，四种模型已开始收敛到了真实模型附近，表明该反演方法真实有效㊂同时，从图可知，沉积层为各向异性介质时，模型收敛速度较慢㊂在粗糙度变化曲线图中，随着迭代次数的增加，模型的粗糙度在总体上逐渐增大，与
R M S 呈现出相反的变化趋势
㊂
图2 水平地形模型反演结果
图3 R M S 和R o u g h n e s s 随迭代次数的变化F i g .3 V a r i a t i o n o f R M S a n d R o u g
h n e s s w i t h i t e r a t i o n t i m e s 002    物探化探计算技术46卷
3.2海床起伏地形模型
设计如图4(a)的海洋模型以分析起伏地形下的各向异性反演效果㊂海水层深度为455.59m~ 1167.69m，电阻率为0.3Ω㊃m，储层的大小㊁位置以及场源的激发㊁接收装置的布设均与上述的水平地形模型一致㊂同样构建储层与沉积层各向异性介质模型，参数设置与上述的海床水平地形模型一致，此处不再重复叙述㊂
反演初始模型如下图4(b)所示，对-6.5k m< y<6.5k m,0k m<z<6.5k m的目标区域进行精细三角网格剖分，共生成5063个三角单元和2629个节点，反演初始模型电阻率设定为1.0Ω㊃m，反演结果如图5所示㊂
由图5(a)和图5(b)可以看出，在起伏地形条件下，同样可以较为准确地恢复储层的电阻率值㊁位置和形态特征，且垂直各向异性模型的恢复效果要优于水平各向异性模型㊂由图5(c)可见，沉积层为垂直各向异性介质时，储层范围在垂直方向上有所缩减，在水平方向上的延伸更为明显，同时在其周围也出现了 伪影 ㊂图5(d)中，储层位置相比于真实模型有所偏移，且向垂直方向拉伸，在其下方出现了虚假低阻异常和分层现象㊂由此可以判断，起伏地形对储层各向异性模型的影响较小，而对沉积层各向异性模型的影响较大㊂
图6为海底起伏地形模型的R M S和R o u g h-n e s s随迭代过程的变化曲线㊂对比这四种模型的R M S变化曲线可以看出，沉积层为各向异性介质时，其收敛速度较慢，需要更多次的迭代计算㊂粗糙度变化曲线图中，相比于海床水平地形模型，沉积层为各向异性介质时对应的粗糙度曲线有较大幅度的起伏变化
㊂
图4起伏地形模型
F i g.4 U n d u l a t i n
g t o p o g r a p h y m o d e
l
图5起伏地形模型反演结果
F i g.5I n v e r s i o n r e s u l t s o f u n d u l a t i n g t o p o g r a p h y m o d e l 102
2期覃金生：起伏地形的2.5维主轴各向异性海洋C S E M反演研究

版权声明：本站内容均来自互联网，仅供演示用，请勿用于商业和其他非法用途。如果侵犯了您的权益请与我们联系QQ:729038198，我们将在24小时内删除。