image
Home Journals Progress in Geophysics
Progress in Geophysics

Abbreviation (ISO4): Prog Geophy      Editor in chief:

About  /  Aim & scope  /  Editorial board  /  Indexed  /  Contact  / 
Online first  /  Current issue  /  All issues  /  Top access  /  RSS

Progress in Geophysics >

2024 , Vol. 39 >Issue 4: 1565 - 1585

DOI: https://doi.org/10.6038/pg2024HH0474

Forward modeling study of three-dimensional time-lapse magnetotelluric method considering induced polarization effect

  • BoShuai DAI ,
  • Xin HUANG ,
  • XiaoYue CAO , * ,
  • LiangJun YAN ,
  • XingBing XIE
Expand

Received date: 2023-10-25

  Online published: 2024-09-29

Copyright

Copyright ©2024 Progress in Geophysics. All rights reserved.
Fold

Abstract

The Magnetotelluric (MT) method is widely used in the exploration of oil and gas resources. The time-lapse MT method can monitor reservoir dynamic distribution and interface changes by observing the time-lapse MT response caused by the change of underground electrical structures and interfaces. Traditional time-lapse MT methods monitoring simulation are based on isotropic theory and regular grids. However, the induced polarization (IP) effect is widely present in reservoirs, ignoring which can lead to interpretation errors and difficulties in time-lapse monitoring of reservoirs, especially for unconventional reservoirs with complex terrain and high water content in the later phases of development. First of all, in order to solve these problems, the Cole-Cole complex resistivity model is used in the vector Helmholtz equation to characterize the IP effect of the reservoir medium. Then, the vector Helmholtz equation discretized by the Galerkin finite element method with unstructured tetrahedral grids, has been implemented to a 3D time-lapse electromagnetic algorithm that considers the IP effect. Thirdly, the analytical solution of the ID model considering IP effect is used to verify the correctness of the 3D algorithm in this paper. Further, different IP parameter variations are set and analyzed to determine the characteristic responses of the reservoirs and discuss the resolution ability of the time-lapse electromagnetic monitoring. At last, a realistic reservoir model in the Fuling area is set up to analyze the time-lapse electromagnetic anomalies generated before and after the displacement of the reservoir. The results show that the significant differences in the time-lapse MT responses caused by changes in the IP parameters of the reservoir can be used to infer the time-lapse change process of the reservoir. In the process of reservoir displacement enhancement, the response difference is distinguished by the sensitivity to the boundary of the displacement swept region, where △ρxya can effectively respond to the change of the interface in the x direction, while △ρyxa targets the y-direction.

Cite this article

BoShuai DAI , Xin HUANG , XiaoYue CAO , LiangJun YAN , XingBing XIE . Forward modeling study of three-dimensional time-lapse magnetotelluric method considering induced polarization effect[J]. Progress in Geophysics, 2024 , 39(4) : 1565 -1585 . DOI: 10.6038/pg2024HH0474

0 引言

为满足油气资源的需求，对剩余油气的时移监测已成为提高油气田采收率和经济效益的有效方法之一(刘文岭和韩大匡，2022).时移地球物理方法主要被应用于监测地下介质随时间的变化情况，目前最为成熟的时移地震技术，由于水力压裂过程的现场噪声，微地震往往难以获得较高的信噪比(石玉梅等, 2004; He et al., 2012)，而且微地震既不直接反映流体包裹体和裂缝连通性，也不能提供响应信息(Peacock et al., 2012; Wang et al., 2016; Zhao et al., 2019).油藏的驱替和改造过程中会形成低阻高极化激电异常，因此，时移电磁法监测相比微地震类方法在流体特征识别上具有独特优势，已成为国内外非常规储层监测的研究热点(屈文璋和安志国, 2020).Wirianto等(2010)设置一系列复杂三维电性模型对开采过程中的油气藏变化进行数值模拟，证明在陆地上CSEM方法对油气藏变化的监测是可行的.谢兴兵等(2016)利用长偏移距瞬变电磁法开展剩余油监测实验，实验结果初步表明时移长偏移距瞬变电磁法在剩余油边界的探测和预测上具有良好的应用前景.Yan等(2018)在涪陵地区开展了页岩压裂过程的连续时域电磁法动态监测实验，结果表明该方法能有效观测到压裂液引起的电磁信号变化.Liu等(2020)用可控源电磁法(CSEM)对涪陵页岩气田的水力压裂过程进行了数值模拟和现场试验，模拟结果和现场数据表明，CSEM在页岩气生产中，监测水力压裂过程具有良好的潜力.相比以上人工源监测方法，天然场源的大地电磁法(Magnetotelluric, MT)探测深度大、效率高，对流体导电性和裂缝连通性变化的敏感度能有效反映流体驱替过程中的电性时移变化，已在油藏监测方面取得良好的效果(Rosas-Carbajal et al., 2015; Thiel, 2017; Bai et al., 2022).彭荣华等(2011)用三维有限差分法计算了时移MT模型的正演响应，数值模拟结果表明时移MT能有效地反映地下介质电性随时间的动态变化.胡祖志等(2014)在柴达木盆地涩北气田开展了时移大地电磁测深实验，监测结果证明了时移MT监测技术的可行性.Thiel (2017)数值正演模拟研究强调了裂缝网络的渗透阀值对电阻率和渗透率的重要性.论证了在水力压裂过程中，这可能解释地表电磁测量的时间变化原因.流体泵入地层后，时移MT监测会产生明显的变化，且这些变化具有方向性，注入前的时移MT正演研究对储层流体监测识别、改造评估可行性和设计施工意义重大.Rosas-Carbajal等(2015)在南澳大利亚帕拉拉纳的增强型地热系统的注入实验中，与三维确定性时移MT反演数据进行比较，他们认为从非均质性可以推断出时间推移的变化.Bai等(2022)基于有限差分法研究了注水压裂中MT响应，利用残差相位张量法计算裂缝走向和倾角特征，分析了不同角度和电阻率条件下在地热储层调查中的可行性.研究表明时移MT监测可以为裂缝的分布方向和发育方向提供可靠的科学依据.胡琪璇等(2023)基于有限差分法，在三维AMT反演理论上，通过加入相邻时刻模型约束，实现了三维时移AMT反演，理论模型反演结果表明时移反演可以更好地聚焦异常体的位置.
但以上油气藏时移大地电磁监测的研究很少考虑激电效应(Induced Polarization, IP)，然而油气动态改造过程(例如驱替和压裂过程)中会形成明显的低阻富激化特征，忽略激电效应会对实际油藏时移监测解释准确性造成影响，甚至导致油气动态监测解释错误.
当前应用于油田监测的时移MT算法，多采用有限差分和规则网格理论(彭荣华等，2011；张弘宇等，2019).然而，这些算法都没有考虑激电效应对油藏动态改造过程的影响.MT属于频率域激发极化法的一个分支，有激电效应和电磁感应双重影响.由于天然源激发极化法复杂的电化学过程，许多学者研究并提出了近十多种地电模型(刘明等，2022).其中，Pelton等(1978)提出随频率变化的Cole-Cole复电阻率模型能较好的解释激电现象，此时的岩石电阻率是一个随频率变化的复电阻率.Burtman和Zhdaov (2015)提出了一个极化团簇的概念模型来解释激电现象.研究结果表明含烃砂岩、砂岩样品和碳酸盐岩具有明显的激电响应特征，实验观测结果与基于GEMTIP方法的理论建模结果进行了对比，验证了之前应用激电法进行油气勘探的地球物理实验结果.Tong等(2020)对Burtman基于GEMTIP理论建模结果进行了改进，提出了MGEOTIP模型，实验结果表明该模型可以更准确地预测非均质介质的零频电阻率和极化率.然而，以上研究均应用于岩石物理测量和分析上，并未能深入应用于三维大地电磁测深法正反演理论和应用研究.万伟等(2019)在ModEM的基础上研究了陆地可控源考虑激电效应的正演.张鑫崇和殷长春(2023)基于Debye模型研究了时间域航空电磁的正演响应特征.同时，传统的MT算法基于规则化网格剖分的不同维度模型进行正反演模型算法研究(Newman and Alumbaugh, 2000; Egbert and Kelbert, 2012)，对于复杂地形，规则化网格难以灵活构造，更加难以精确揭示真实的三维地下复杂电性响应(赵晓博等，2011；Cao et al., 2018)，基于非结构化网格的有限元法能够用更少的单元，更准确地表示任意地形结构(Jahandari and Farquharson, 2017).
综上所述，本文采用基于非结构四面体单元的有限元法实现时移MT正演算法，首先推导了基于Cole-Cole模型的三维MT有限元法基本理论，并通过对比考虑激电效应的一维解析解，验证了本文算法的正确性.其次，通过理论模型研究了不同激电参数下，异常体的响应特征.最后，结合理论模型和涪陵页岩气田区焦页30井地区的实际模型，研究了不同驱替油藏时段的三维时移MT响应特征.研究结果证明了三维时移大地电磁法在油气藏动态监测中的可行性.

1 三维时移MT有限元正演算法

1.1 控制方程

针对三维MT正演，假设时谐因子为e－iωt, 并忽略位移电流，则电场满足的矢量亥姆霍兹方程为：
式中：E为电场强度；ω为角频率；μ为磁导率.其中，代表包含激电的复数电导率，本文采用Cole-Cole模型定义复电阻率(Pelton et al.，1978)：
其中，，ρ为真电阻率；ρ0为零频电阻率；m为极化率；τ为时间常数；c为频率相关系数.

1.2 非结构有限元法

方程式(1)利用伽辽金法获得有限元法弱形式为：
其中，Ω为离散区域，Φ为非结构四面体矢量基函数.
采用四面体网格和对应的矢量基函数进行离散，将离散区域Ω剖分为有限个不重叠的四面体单元.图 1为四面体单元的局部节点和棱边编号.
图1 四面体单元节点和棱边示意图

Fig 1 Schematic diagram of tetrahedral element nodes and edges

将有限维函数代入式(3)可得：
其中Eh为有限元法的近似解，Ei为单元棱边上的电场值，M为计算域Ω内棱边的个数. 在强加狄利克雷边界条件后，可获得大型复线性方程组：
其中，K为有限元单位总体矩阵，Sb为与边界条件相关的右端项.

1.3 时移大地电磁响应函数

通过求解方程组(5)，可以得到任意单元棱边上的电场值，则区域内的任意电场值可通过上述有限维函数插值得到.对应磁场H通过法拉第定理获得：
MT的响应函数视电阻率、相位及其对应的时移响应可通过以下方程获得：
其中，i和j分别为x或y方向，Zij为阻抗，()t0和()t1分别为原始阶段和t1阶段. 式(7)的阻抗分量可以通过求解两个极化方向的电磁场得到：
其中，1和2表示不同的极化方向.

2 数值算例

为验证本文算法的正确性，设计一个如图 2所示的3层层状介质模型(K型)，通过对比一维MT解析解与三维有限元正演数值解，验证本文三维考虑激电效应MT有限元法正演的正确性，层状介质模型参数如表 1所示.如图 3a所示，对比本文三维有限元正演结果与一维解析解视电阻率ρxya和相位φxy结果发现，两者吻合较好.图 3b表明，本文算法视电阻率ρxya和相位φxy相对解析解的最大误差均小于5%, 验证了考虑激电效应的三维MT非结构有限元算法的精确性.
图2 三层激电介质模型图

Fig 2 Diagram of the three-layer IP model

表1 层状介质模型参数

Table 1 Parameters of the layered medium model

地层序号 真电阻率
ρ/(Ω·m)		 厚度/m 极化率m 时间常数τ/s 频率相关系数c
1 100 2000 0.0 0.0 0.0
2 10 1000 0.4 100.0 0.4
3 100 — 0.0 0.0 0.0
图3 (a) 本文算法与解析解对比; (b)相对误差

Fig 3 (a) Comparison between the proposed algorithm and analytical solutions; (b) Relative errors

2.1 油藏模型的激电响应分析

为研究深层非常规油气储层MT激电响应特征，本文设置如图 4所示模型，在电阻率为100 Ω·m的均匀半空间中，存在一个有含激电效应的异常体，中心埋深为2.5 km，大小为2 km×2 km×1 km的油藏，设置3组可变激电参数模型如表 2所示.本文正演的计算频率范围为10－4到10－1 Hz，按对数等间距取4个频率，地表的测点范围为-4 km到4 km，点距和线距均间隔为200 m，共计1681个测点.
图4 深层非常规油气藏激电模型

Fig 4 IP models for the deep unconventional reservoir

表2 油气藏模型激电参数

Table 2 IP parameters of the reservoir models

模型 真电阻率ρ
/(Ω·m)		 极化率m 时间常数τ/s 频率相关系数c
模型一 10.0 0.8, 0.6, 0.4, 0.2 100.0 0.4
模型二 10.0 0.4 102, 101, 10-1, 10-2 0.4
模型三 10.0 0.4 100.0 0.6, 0.4, 0.2, 0.1
本文针对如表 2中的3组不同变化的激电参数正演模型，对含激电效应的非常规油藏异常体三维MT特征响应进行研究.

2.1.1 不同极化率m的有限元MT响应

设置如表 2中模型一，在地表下电阻率为100 Ω·m的均匀半空间中，存在真电阻率ρ为10 Ω·m、时间常数τ为100.0 s、频率相关系数c为0.4、极化率m分别为0.8、0.6、0.4和0.2时4种情况的激电异常体.
图 5a—d依次显示了模型一，在4种频率10－1、10－2、10－3和10－4 Hz下，极化率m为0.8、0.6、0.4和0.2时，视电阻率ρxya和ρyxa在主侧线y=0时曲线图.图 5结果显示，在所有频率中异常体中心上随着极化率m的减小，视电阻率ρxya的幅值也减小；不同m时的ρxya幅值变化随着远离异常中心而增大，在异常体边界外出现交叉现象.而在所有频率中，视电阻率ρyxa数据随着极化率m的减小，视电阻率的幅值均减小.在相同频率时，ρxya和ρyxa在异常体的正上方不同m之间的变化值最大，且它们的最大变化值在每个频率上大致相等，变化量随着测点位置远离异常中心而变小；不同频率时，ρxya和ρyxa的变化值随频率减小而增大.
图5 视电阻率ρxya和ρyxa曲线图

(a)10－1 Hz; (b)10－2 Hz; (c)10－3 Hz; (d)10－4 Hz.

Fig 5 Graphs of the apparent resistivityand

从异常中心到远离异常体的位置，视电阻率ρxya的变化大于ρyxa，其中在10－2 Hz和m=0.2时，ρxya从83.63 Ω·m变化到108.89 Ω·m，而视电阻率ρyxa从83.80 Ω·m变化到100.64 Ω·m.随着频率降低，极化率m的变化产生的视电阻率异常特征更显著.
图 6a—d和图 7a—d为10－2 Hz时不同极化率m的视电阻率ρxya和ρyxa对应的等值线图，其中白色方框为异常体位置.图 6和图 7显示，在异常体上方随着m变大，ρxya和ρyxa均变大.图 6显示在异常体周边，ρxya在x方向变化梯级大，在y方向具有拉伸的趋势.而图 7显示，在异常体周边ρyxa在y方向变化梯级大，在x方向具有拉伸的趋势.ρxya和ρyxa变化梯级最大的地方恰好分别对应异常体的边界位置.这种现象很容易理解.事实上，ρxya和ρyxa的响应分别由x和y方向的极化电场决定(曹晓月等，2018).
图6 频率10－2 Hz视电阻率ρxya等值线图

Fig 6 Contour maps of apparent resistivity ρxya at frequency 10－2 Hz

(a)m=0.8; (b)m=0.6; (c)m=0.4; (d)m=0.2.

图7 频率10－2 Hz视电阻率ρyxa等值线图

Fig 7 Contour maps of apparent resistivity ρyxa at frequency 10－2 Hz

(a)m=0.8; (b)m=0.6; (c)m=0.4; (d)m=0.2.

x方向极化电场时，极化电流从x方向穿过异常体，由于电流密度在界面的垂向连续性，异常体电阻率10 Ω·m远小于围岩的100 Ω·m，由欧姆定律导致异常体内的电场远小于围岩内，从而异常体界面处极化电场变化最大，因此异常体上方ρxya出现明显的低阻异常，其变化梯级最大处对应异常体的x界面.同理y方向极化电场时，在异常体上方ρyxa出现明显的低阻异常，而变化梯级最大处对应异常体的y界面.由于ρxya、ρyxa对x和y方向界面不同敏感性的特点，可以由不同视电阻率特征提取极化异常体不同方向的边界信息.

2.1.2 不同时间常数τ的有限元MT响应

设置如表 2中的模型二，在100 Ω·m的均匀半空间中，存在一个真电阻率ρ为10 Ω·m，极化率m，频率相关系数c均为0.4，时间常数τ分别取102、101、10－1和10－2 s时的4种情况的激电异常体.
图 8a—d依次显示了模型二，在4种频率10－1、10－2、10－3和10－4 Hz下，时间常数为102、101、10－1和10－2 s时，视电阻率ρxya和视电阻率ρyxa在主侧线y=0时曲线图.图 8结果显示，在所有频率中异常体中心上随着时间常数τ的增大，视电阻率ρxya的幅值减小；不同τ时的ρxya幅值随着远离异常中心而变大，在异常体边界外出现交叉现象.而在所有频率中，视电阻率ρyxa数据随着τ的增大，视电阻率的幅值均逐渐变小.在相同频率时，ρxya和ρyxa在异常体的正上方不同τ之间的变化值最大；变化量随着测点位置远离异常中心而变小；不同频率时，ρxya和ρyxa的变化值随频率减小而减小.从异常中心到远离异常体的位置，视电阻率ρxya的变化大于ρyxa，其中在10－2 Hz和τ=102 s时，ρxya从84.69 Ω·m变化到108.64 Ω·m，而视电阻率ρyxa从84.88 Ω·m变化到100.86 Ω·m.随着频率的降低，时间常数τ的变化产生的视电阻率异常特征越小.
图8 视电阻率ρxya和ρyxa曲线图

Fig 8 Graphs of apparent resistivity ρxya and ρyxa

(a)10－1 Hz; (b)10－2 Hz; (c)10－3 Hz; (d)10－4 Hz.

图 9a—d和图 10a—d为10－2 Hz时不同时间常数τ的视电阻率ρxya和ρyxa对应的等值线图，其中白色方框为异常体位置.图 9和图 10显示，在异常体上方随着τ的减少，ρxya和ρyxa均有所变大，然而变化幅值较小.图 9同图 6的特征一致，在异常体周边ρxya在x方向变化梯级大，在y方向具有拉伸趋势.而图 10同图 7的特征一致，在异常体周边ρyxa在y方向变化梯级大，在x方向具有拉伸趋势.ρxya和ρyxa变化梯级最大的地方恰好也分别对应异常体的边界位置.结合本文正演算法ρxya和ρyxa的等值线图上x和y异常梯级带，τ类比m变化时，提取的异常体边界信息均与实际基本一致.本节算例进一步说明，可以用ρxya和ρyxa的等值线梯级带去分别定性判断极化异常体的x和y界面.对比模型一和模型二的视电阻率曲线和等值线图，可知在对有激电效应的异常体研究时，τ的参数变化对实际的激电特征响应差异小于m，说明该条件下τ的影响小于m.
图9 频率10－2 Hz视电阻率ρxya等值线图

Fig 9 Contour maps of apparent resistivity ρxya at frequency 10－2 Hz

(a)τ=102 s; (b)τ=101 s; (c)τ=10－1 s; (d)τ=10－2 s.

图10 频率10－2 Hz视电阻率ρyxa等值线图

Fig 10 Contour maps of apparent resistivity ρyxa at frequency 10－2 Hz

(a)τ=102 s; (b)τ=101 s; (c)τ=10－1 s; (d)τ=10－2 s.

2.1.3 不同频率相关系数c的有限元MT响应

设置如表 2中模型三，在地表下电阻率为100 Ω·m的均匀半空间中，存在一个真电阻率ρ为10 Ω·m、时间常数τ为102 s、极化率m为0.4、频率相关系数c分别为0.6、0.4、0.2和0.1时4种情况的激电异常体.
图 11a—d依次显示了模型三，在4种频率10－1、10－2、10－3和10－4 Hz下，频率相关系数c为0.6、0.4、0.2和0.1时，视电阻率ρxya和ρyxa在主侧线y=0时曲线图.图 11结果显示，在频率10－1 Hz和10－2 Hz中，异常体中心上视电阻率ρxya和ρyxa随着频率相关系数c的减小，其幅值增大；频率10－3 Hz异常体中心上视电阻率ρxya和ρyxa随着频率相关系数c的减小，其幅值变化很小；而在频率10－4 Hz中，异常体中心上视电阻率ρxya和ρyxa随着频率相关系数c的减小，其幅值减小.在相同频率时，ρxya和ρyxa在异常体的正上方不同频率相关系数c之间的变化值最大，变化量随着测点位置远离异常中心而变小.从异常中心到远离异常体的位置，视电阻率ρxya的变化大于ρyxa，其中在10－2 Hz和c=0.1时，ρxya从85.44 Ω·m变化到108.45 Ω·m，而视电阻率ρyxa从85.64 Ω·m变化到101.01 Ω·m.随着频率的降低，频率相关系数c的变化产生的视电阻率异常特征区别不明显.
图11 视电阻率ρxya和ρyxa曲线图

Fig 11 Graphs of apparent resistivity ρxya and ρyxa

(a)10－1 Hz; (b)10－2 Hz; (c)10－3 Hz; (d)10－4 Hz.

图 12a—d和图 13a—d为10－2 Hz时不同频率相关系数c的视电阻率ρxya和ρyxa对应的等值线图，其中白色方框为异常体位置.图 12和图 13显示，在异常体上方出现低阻异常，异常体周边趋近于均匀半空间电阻率.然而，随着c的减少，ρxya和ρyxa变化幅值微小.同模型一和模型二一致，本节算例可以用ρxya和ρyxa的等值线梯级带获得极化异常体的x和y界面.类似于模型二，对比模型一和模型三的视电阻率曲线和等值线图，可知在对有激电效应的异常体研究时，c的参数变化对实际的激电特征响应差异也小于m，说明该条件下c的影响也小于m.
图12 频率10－2 Hz视电阻率ρxya等值线图

Fig 12 Contour maps of apparent resistivity ρxya at frequency 10－2 Hz

(a)c=0.6; (b)c=0.4; (c)c=0.2; (d)c=0.1.

图13 频率10－2 Hz视电阻率ρyxa等值线图

Fig 13 Contour maps of apparent resistivity ρyxa at frequency 10－2 Hz

(a)c=0.6; (b)c=0.4; (c)c=0.2; (d)c=0.1.

2.2 油藏动态监测模拟理论算例

油气勘探领域中目标体埋深更深，极化效应更强，且地下介质随着开采不断变化，因此有必要研究石油介质的激电效应时移变化，通过组合模型模拟出石油驱替开采的过程，分析驱替液在油藏驱替动态过程的大地电磁法时移响应，设置理论的驱替组合模型模拟油气开发过程，如图 14所示，在地下电阻率为100 Ω·m的均匀半空间中，存在一个长方体形状具有激电效应的非常规储层，大小为2 km×2 km×1 km，顶面埋深为2 km，储层的真电阻率为10 Ω·m，其激电参数为m=0.2，τ=0.1 s，c=0.1.其驱替过程的四个阶段如图 14a、b所示，驱替液的电阻率为0.1 Ω·m，激电参数为m=0.6，τ=150 s，c=0.6.考虑到2.1节研究的带激电效应异常体的响应特征，本文的时移监测模拟针对10－2 Hz进行研究.地表时移监测网范围为8 km×8 km，点距和线距均间隔320 m，共676个测点.
图14 深层非常规油气藏激电模型

(a)xy方向断剖面图; (b)xz方向断剖面图.

Fig 14 IP models for the deep unconventional reservoir

(a) Section diagram in xy direction; (b) Section diagram in xz direction.

本文设置5种不同驱替状态的组合模型去模拟油气的开采阶段，即整个储层驱替过程的0%(0 km×2 km×1 km)、25%(0.5 km×2 km×1 km)、50%(1 km×2 km×1 km)、75%(1.5 km×2 km×1 km)和100%(2 km×2 km×1 km)，来分析4个阶段的时移视电阻率Δρxya和Δρyxa响应.图 15为四个驱替阶段上y=0测线上的Δρxya和Δρyxa响应曲线图.
图15 驱替阶段时移视电阻率响应曲线图

(a)Δρxya曲线图; (b)Δρyxa曲线图.

Fig 15 Time-lapse apparent resistivity response graphs of displacement stage

(a) Δρxya graphs; (b) Δρyxa graphs.

图 15a的曲线图表明，不管有无激电效应的油藏层驱替模型，随着驱替阶段的增进，无驱替与不同驱替阶段的时移电阻率差值Δρxya不断增大，考虑激电效应后的驱替模型视电阻率差值响应峰值更加拟合目标异常体的中心位置.同时随着驱替波及体积变大，Δρxya异常在x轴上沿着驱替界面不断推进，考虑激电效应后的驱替趋势状态更加明显，更能反映驱替方向.考虑激电效应后，Δρxya曲线的四个峰值位置恰好和驱替模型四个阶段的中心位置重合，比不含激电效应的中心拟合度更好.图 15b的曲线图表明，与图 15a类似，随着驱替阶段的增进，Δρyxa也不断增大，同时随着驱替波及区推进，沿着x轴不断地扩展.在考虑激电效应的油藏层驱替模型后，Δρyxa的四个峰值位置也恰好和驱替模型四个阶段的中心位置吻合.说明随着不同时间段的驱替液推进，驱替波及区的中心位置和走向可以由Δρxya和Δρyxa异常的峰值位置和推进趋势来确定.图 15a和图 15b不同点在于未完全驱替的三个阶段，Δρyxa始终大于Δρxya，在驱替100%阶段二者大致相等.这种现象很容易理解.由2.1节的结论，ρxya和ρyxa的响应分别对x和y方向的识别界面最敏感.在驱替小于100%阶段，驱替波及区y方向(长度2 km)始终大于x方向(长度<=75%×2 km)，由此ρyxa响应会大于ρxya.图 15曲线图表明，通过分析4种驱替阶段的Δρxya和Δρyxa，可以发现，随着驱替油藏的进行，考虑激电效应的时移响应异常能准确地反映驱替的每段过程，对异常中心位置拟合更加准确.Δρxya对x方向增进的识别能力强，Δρyxa对y方向增进的识别能力强.两者结合可以确定驱替变化在水平方向的中心位置和边界变化信息，通过观察异常体响应的大小，判断压裂是否充分.通过观察每阶段的过程，分析驱替走向信息，减少施工成本，增加施工效率，对指导油气工程中方向性施工和开采具有一定的意义.
图 16a—d和图 17a—d分别为10－2 Hz时，有无激电效应的驱替阶段时移视电阻率Δρxya和Δρyxa对应的等值线图，其中白色方框对应不同驱替阶段(50%和100%)的波及区域.有激电效应的视电阻率差值比没有激电效应的视电阻率差值小，说明激电效应会影响驱替过程电性变化.图 16显示，驱替不同阶段Δρxya在驱替波及区上方最大，并随着驱替推进，两个Δρxya的异常中心均对应了驱替波及区的中心位置，并随着驱替的推进Δρxya的极大值也不断增加.图 16a—d显示Δρxya异常的梯级带与驱替过程的波及区x方向的边界具有良好的对应，进一步说明了Δρxya对x方向增进的识别能力强.图 17显示，驱替不同阶段Δρyxa在驱替波及区上方出现异常，随着驱替推进，两个Δρyxa的异常中心也沿着x方向推进，并随着驱替的推进Δρyxa的极大值也不断增加.然而，图 17a—d显示异常的梯级带在y方向，且在驱替过程中大致不变，这从反面进一步说明由于驱替只在x方向推进，Δρyxa对y方向增进的识别能力强，所以Δρyxa异常的梯级带在y方向位置基本保持不变.结合图 16和图 17的结果，进一步证明了考虑激电效应的时移视电阻率响应可以监测驱替的走向和边界.
图16 驱替阶段Δρxya等值线图

(a)驱替50%; (b)驱替100%; (c)有IP驱替50%; (d)有IP驱替100%.

Fig 16 Contour maps of Δρxyaduring displacement stages

(a) Displacement 50%; (b) Displacement 100%; (c) Displacement 50% with IP; (d) Displacement 100% with IP.

图17 驱替阶段Δρyxa等值线图

(a)驱替50%; (b)驱替100%; (c)有IP驱替50%; (d)有IP驱替100%.

Fig 17 Contour maps of Δρyxaduring displacement stages

(a) Displacement 50%; (b) Displacement 100%; (c) Displacement 50% with IP; (d) Displacement 100% with IP.

3 实际应用

理论模型的驱替响应，为简单条件下的正演模拟，结果相比实际应用侧重不同，为更加有效的证明本算法在实际监测压裂工程中的应用，需要考虑激电效应的理论和算法在复杂环境下的驱替响应，应对不同的地质环境下油气压裂施工.在石油勘探开发的电磁监测中，射孔压裂后，压裂液往往会沿着最大主应力方向推进，垂直为水力压裂的水平井方向，随着压裂的进行，压裂区域沿着水平井方向推进.本节参考Liu等(2020)参考涪陵页岩气田区焦页30井测井资料，本文设计了如图 18所示的该地区页岩气储层模型.本文基于Liu等(2020)建立的模型，采用考虑激电效应的三维时移大地电磁法对该地区的页岩气压裂模型，进行动态监测模拟分析.为研究压裂液在页岩气储层压裂过程中的时移响应，本文设置岩气储层激电模型规模为2300 m×840 m×300 m，参考前文驱替过程将其分为四个压裂阶段.在页岩气储层压裂前油藏的真电阻率为42 Ω·m，其激电参数为m=0.2、τ=0.1 s和c=0.1，压裂后真电阻率为5 Ω·m，激电参数为m=0.6、τ=150 s和c=0.6.本模型的正演计算频率为0.01 Hz.地表处测点覆盖范围为8 km×8 km，点距和线距均间隔320 m，共676个测点.
图18 涪陵页岩气田区焦页30井模型剖面图

Fig 18 Cross-section of Jiaoye 30 well model in Fuling shale gas field

图 19为压裂阶段时移视电阻率响应曲线图，异常响应比等值线图简单直接，反映驱替的中心位置走向，清晰地反映压裂液在油藏区不同百分段驱替的过程响应.图 19a的曲线图表明，随着压裂阶段的增进，压裂前与压裂后各阶段的时移电阻率差值Δρxya不断增大，从第一段的8.85 Ω·m增大到第四段的29.81 Ω·m.同时随着压裂的推进，Δρxya异常在x方向上沿着不断推进.而且Δρxya曲线的四个峰值位置恰好和压裂四个阶段的中心位置比较吻合.图 19b的曲线图表明，随着压裂的推进，Δρyxa也不断增大，从第一段的8.46 Ω·m增大到第四段的21.72 Ω·m，同时也沿着x方向不断地推进.Δρyxa的四个峰值位置也和压裂四个阶段的中心位置比较吻合.图 19说明随着在实际的页岩气储层压裂动态监测过程中，压裂区的中心位置和走向可以由Δρxya和Δρyxa异常的峰值位置和推进趋势来确定.图 19a和图 19b不同点在压裂改造区上方的Δρyxa除了峰值位置都小于Δρxya，改造区域外的四个阶段，Δρxya都小于Δρyxa.这可以由上文得到的结论来解释.ρxya和ρyxa的响应分别对x和y方向的识别界面最敏感，而压裂的第一个阶段压裂区y方向大(800 m>25%(2300)，而后三个压裂阶段压裂区x方向大(800 m < 50%(2300).图 19曲线图表明，通过分析4种压裂阶段的Δρxya和Δρyxa，在实际的涪陵页岩气田大地电磁时移监测中，时移响应异常明显，能通过曲线图峰值的位置变化趋势，发现压裂走势是否和水平井的走向一致，判断压裂是否充分，与预期的施工参数是否相同，有助于调整实际油气压裂施工.
图19 压裂阶段时移视电阻率响应曲线图

(a)有IP Δρxya曲线图; (b)有IP Δρyxa曲线图.

Fig 19 Time-lapse apparent resistivity response graphs of fracturing stages

(a) Δρxya graphs with IP; (b) Δρyxa graphs with IP.

由图 20和图 21的压裂阶段时移等值线图可以看出.图 20反映的是Δρxya的响应，由于对x方向电阻率敏感，y方向会产生偏大于实际大小的异常响应，可以通过x方向的梯级带，判断压裂阶段x方向的边界.同理在图 21中，可根据Δρyxa的响应y方向的梯级带，判断压裂阶段y方向的边界.图 20和图 21，都表明随着压裂过程的不断推进，压裂区时移异常逐渐变大.这是由于压裂区储层的视电阻率，不断被富极化且低阻的压裂液灌进，压裂区相对于围岩的激化电性差异不断变大.图 20和图 21的结果，进一步证明了时移视电阻率响应可以用于实际储层条件下时移监测压裂施工，通过等值线图的直观响应，可以分析出压裂液驱替前后的位置和方向差异，判断不同阶段的驱替程度，通过边界的变化综合分析每段压裂效果是否充分，以及压裂方向是否符合预期.
图20 压裂阶段Δρxya等值线图

(a)有IP压裂25%; (b)有IP压裂50%; (c)有IP压裂75%; (d)有IP压裂100%.

Fig 20 Contour maps of Δρxya during fracturing stages

(a) Fracturing 25% with IP; (b) Fracturing 50% with IP; (c) Fracturing 75% with IP; (d) Fracturing 100% with IP.

图21 压裂阶段Δρyxa等值线图

(a)有IP压裂25%; (b)有IP压裂50%; (c)有IP压裂75%; (d)有IP压裂100%.

Fig 21 Contour maps of Δρyxa during fracturing stages

(a) Fracturing 25% with IP; (b) Fracturing 50% with IP; (c) Fracturing 75% with IP; (d) Fracturing 100% with IP.

4 结论

本文通过引入Cole-Cole复电阻率模型并结合非结构四面体单元的有限元法，开发了一种考虑激电效应的三维时移大地电磁法正演方法.首先，通过一维MT解析解，验证了本文算法的正确性.其次针对典型三维油藏激电模型，研究了不同激电参数的响应特征，并分析了油藏时移监测异常和识别方法.最后，将本文算法应用到涪陵页岩气田实际储层条件的激电模型三维时移模拟和分析.通过理论和实际油藏模型算例，得出以下结论.
(1) 激电参数影响三维大地电磁测深法视电阻率响应，在油气藏大地电磁时移监测中极化率m的影响最为明显，而时间常数τ和频率相关系数c较m影响较小.
(2) 在油气藏驱替和压裂动态变化中，时移视电阻率响应可以有效的表征储层激电参数变化，进而来监测储层动态变化的空间展布.
(3)Δρxya对油藏x方向上的变化更为敏感，而Δρyxa对油藏y方向上的变化更为敏感.在储层动态改造中，可以综合分析二者的特征，对储层的时移变化在x和y方向的展布进行有效地分析.
本文研究成果可为三维时移大地电磁测深法在油气藏动态监测中的应用，提供理论基础和参考.

感谢长江大学EMLAB研究团队成员在文章准备过程中提供的帮助，感谢审稿专家对本文提出的修改建议.

References

Publishing order | Descend order by publishing year | Descend order by cited within
Bai L G , Li J , Zeng Z F . Numerical simulation of hot dry rock fracture monitoring by time-lapse magnetotelluric method. Energies, 2022, 15(19): 7203

DOI

Burtman V , Zhdaov M S . Induced polarization effect in reservoir rocks and its modeling based on generalized effective-medium theory. Resource-Efficient Technologies, 2015, 1(1): 34 48

DOI

Cao X Y , Yin C C , Zhang B . A goal-oriented adaptive finite-element method for 3D MT anisotropic modeling with topography. Chinese Journal of Geophysics, 2018, 61(6): 2618-2628

DOI

Cao X Y , Yin C C , Zhang B . 3D magnetotelluric inversions with unstructured finite-element and limited-memory quasi-Newton methods. Applied Geophysics, 2018, 15(3-4): 556-565

DOI

Egbert G D , Kelbert A . Computational recipes for electromagnetic inverse problems. Geophysical Journal International, 2012, 189(1): 251-267

DOI

He L F , Hu X M , Xu L G . Feasibility of monitoring hydraulic fracturing using time-lapse audio-magnetotellurics. Geophysics, 2012, 77(4): WB119-WB126

DOI

Hu Q X , Tan H D , Yu C . 3D Inversion of time-lapse controlled source audio-frequency magnetotellurics. Geoscience, 2023, 37(1): 90-98

DOI

Hu Z Z , He Z X , Li D C . Reservoir monitoring feasibility study with time lapse magnetotelluric survey in Sebei Gas Field. Oil Geophysical Prospecting, 2014, 49(5): 997-1005

DOI

Jahandari H , Farquharson C G . 3-D minimum-structure inversion of magnetotelluric data using the finite-element method and tetrahedral grids. Geophysical Journal International, 2017, 211(2): 1189-1205

DOI

Liu M , Li W B , Zhu Z L . Effective parameters of natural field induced polarization method based on Dias model. Science Technology and Engineering, 2022, 22(3): 939-944

DOI

Liu R , Liu J X , Wang J X . A time-lapse CSEM monitoring study for hydraulic fracturing in shale gas reservoir. Marine and Petroleum Geology, 2020, 120 104545

DOI

Liu W L , Han D K . Digital twin system of oil and gas reservoirs: a new direction for smart oil and gas field construction. Acta Petrolei Sinica, 2022, 43(10): 1450-1461

DOI

Newman G A , Alumbaugh D L . Three-dimensional magnetotelluric inversion using non-linear conjugate gradients. Geophysical Journal International, 2000, 140(2): 410-424

DOI

Peacock J R , Thiel S , Reid P . Magnetotelluric monitoring of a fluid injection: example from an enhanced geothermal system. Geophysical Research Letters, 2012, 39(18): L18403

DOI

Pelton W H , Ward S H , Hallof P G . Mineral discrimination and removal of inductive coupling with multifrequency IP. Geophysics, 1978, 43(3): 588-609

DOI

Qu W Z , An Z G . Numerical simulation of time-lapse audio magnetotelluric monitoring. Progress in Geophysics, 2020, 35(4): 1595-1604

DOI

Rosas-Carbajal M , Linde N , Peacock J . Probabilistic 3-D time-lapse inversion of magnetotelluric data: application to an enhanced geothermal system. Geophysical Journal International, 2015, 203(3): 1946-1960

DOI

Shi Y M , Liu W L , Yao F C . Effects of signal-to-noise ratio on seismic reservoir monitoring. Petroleum Exploration and Development, 2004, 31(S1): 113-116

Thiel S . Electromagnetic monitoring of hydraulic fracturing: relationship to permeability, seismicity, and stress. Surveys in Geophysics, 2017, 38(5): 1133-1169

DOI

Tong X L , Yan L J , Xiang K . Modifying the generalized effective-medium theory of induced polarization model in compacted rocks. Geophysics, 2020, 85(4): MR245-MR255

DOI

Wan W , Tang X G , Huang Q H . Influence of induced polarization effects on 3D land CSEM sounding. Progress in Geophysics, 2019, 34(6): 2328-2335

DOI

Wang X B , Zhang B , He Z X . Electrical properties of Longmaxi organic-rich shale and its potential applications to shale gas exploration and exploitation. Journal of Natural Gas Science and Engineering, 2016, 36 573-585

DOI

Wirianto M , Mulder W A , Slob E C . A feasibility study of land CSEM reservoir monitoring in a complex 3-D model. Geophysical Journal International, 2010, 181(2): 741-755

Xie X B , Zhou L , Yan L J . Remaining oil detection with time-lapse long offset & window transient electromagnetic sounding. Oil Geophysical Prospecting, 2016, 51(3): 605-612

DOI

Yan L J , Chen X X , Tang H . Continuous TDEM for monitoring shale hydraulic fracturing. Applied Geophysics, 2018, 15(1): 26-34

DOI

Zhang H Y , Wen J L , Zhang F M . Finite difference method two-dimensional geoelectric structure model forward modeling of magnetotelluric method. Mineral Exploration, 2019, 10(8): 1988-1992

DOI

Zhang X C , Yin C C . Forward modeling of airborne electromagnetic induced polarization effect in time-domain based on Debye model. Journal of Jilin University (Earth Science Edition), 2023, 53(5): 1573-1581

DOI

Zhao X B , Zhu Z Q , Li J H . Finite element modeling of 2.5D TEM using unstructured meshes. Computing Techniques for Geophysical and Geochemical Exploration, 2011, 33(5): 517-521

DOI

Zhao Y L , Li N Y , Zhang L H . Productivity analysis of a fractured horizontal well in a shale gas reservoir based on discrete fracture network model. Journal of Hydrodynamics, 2019, 31(3): 552-561

DOI

曹 晓月 , 殷 长春 , 张 博 . 面向目标自适应有限元法的带地形三维大地电磁各向异性正演模拟. 地球物理学报, 2018, 61(6): 2618-2628

DOI

胡 琪璇 , 谭 捍东 , 于 翠 . 时移可控源音频大地电磁法三维反演研究. 现代地质, 2023, 37(1): 90-98

DOI

胡 祖志 , 何 展翔 , 李 德春 . 涩北气藏时移大地电磁监测技术可行性研究. 石油地球物理勘探, 2014, 49(5): 997-1005

DOI

刘 明 , 李 文奔 , 朱 占龙 . 基于Dias模型的天然场激发极化法有效参数. 科学技术与工程, 2022, 22(3): 939-944

DOI

刘 文岭 , 韩 大匡 . 数字孪生油气藏: 智慧油气田建设的新方向. 石油学报, 2022, 43(10): 1450-1461

DOI

屈 文璋 , 安 志国 . 时移音频大地电磁监测数值模拟研究. 地球物理学进展, 2020, 35(4): 1595-1604

DOI

石 玉梅 , 刘 雯林 , 姚 逢昌 . 地震资料信噪比对油藏地震监测的影响. 石油勘探与开发, 2004, 31(S1): 113-116

万 伟 , 唐 新功 , 黄 清华 . 陆地可控源电磁法三维勘探的激电效应影响研究. 地球物理学进展, 2019, 34(6): 2328-2335

DOI

谢 兴兵 , 周 磊 , 严 良俊 . 时移长偏移距瞬变电磁法剩余油监测方法及应用. 石油地球物理勘探, 2016, 51(3): 605-612

DOI

张 弘宇 , 温 建亮 , 张 富明 . 有限差分法二维地电构造模型大地电磁法正演模拟. 矿产勘查, 2019, 10(8): 1988-1992

DOI

张 鑫崇 , 殷 长春 . 基于Debye模型的时间域航空电磁激电效应正演模拟. 吉林大学学报(地球科学版), 2023, 53(5): 1573-1581

DOI

赵 晓博 , 朱 自强 , 李 建慧 . 基于非结构化网格的瞬变电磁2.5维有限元正演模拟. 物探化探计算技术, 2011, 33(5): 517-521

DOI

Options
Outlines

/

〈 〉