地热资源作为新型接替能源, 因其具有可再生、节能及环保等特点, 越来越受到广泛重视^[1]。据中国大陆地区大地热流数据汇编显示, 我国地热资源十分丰富, 具有东高、中低、南高以及西北低等特征^[2-3]。其中, 深层(3~10 km埋深)干热岩地热资源总量为20.9×10⁶ EJ, 相当于714.9×10¹² t标准煤, 是极具开发潜力的战略接替能源^[4]。

应用于地热资源勘探的地球物理方法以间接勘探方法为主, 主要依据岩石物性参数与温度的相关性开展研究。目前, 常用的方法主要有重力法、磁法、电磁法、以及地震方法等^[5-9]。大地电磁法(MT)因具有勘探深度大、工作效率高和成本低廉等优点, 已成为深层地热勘探最常用的方法之一^[10-13]。常规MT资料处理解释主要基于电磁感应理论, 然而, 地下岩石在低频电磁场作用下, 会同时发生电磁感应和激发极化现象, 此时的岩石电阻率是一个随频率变化的复电阻率^[14-18]。多年来, 国内外学者在考虑激电效应的MT正演方面做了很多研究。1978年, PELTON等^[19]基于实验结果, 提出由激电效应引起的复电阻率随频率的变化规律能用Cole-Cole模型表示。随后, 罗延钟等^[20]基于电子导体激发极化的电化学机理, 通过研究过电位充、放电特性和频谱特征, 建立了面极化和体极化的等效电路和数学模型, 从而在理论上证明了采用Cole-Cole模型描述激电复电阻率频谱的合理性。2011年, 朱占升等^[21]采用有限单元法, 通过MT二维正演计算非极化均匀空间中存在的不同极化体, 发现高阻极化体对视电阻率值影响强于低阻极化体。2016年, 符超等^[22]利用视电阻率比值与相位比值, 详细讨论了大地电磁场激发极化效应与极化体埋深和厚度的关系。2018年, 张志勇等^[23]采用超松弛迭代双共轭梯度法, 开展了同时考虑激电效应和电磁效应的复电阻率法二维数值模拟研究。2020年, 熊治涛等^[24]采用有限元数值计算方法, 利用Galerkin加权余量法, 实现了二维MT各向异性极化介质正演, 该研究成果对频率域电磁法野外勘探具有指导意义。

随着基于激电效应的MT正演理论的不断完善, 野外MT实测资料激电参数提取取得了突破性进展。2006年, 陈清礼等^[25]基于岩矿石激电响应的Dias模型, 开展MT资料真谱参数反演方法研究, 成功从MT资料中提取激电参数, 为地质解释提供更多依据。曹中林等^[26]在油气检测中进行MT激电效应的模拟研究, 该方法利用岩石电阻率、极化率、时间常数以及频率相关系数等参数, 同时进行MT资料处理解释, 提高了MT方法的分辨能力和可靠性。HE等^[27]和罗卫峰等^[28]在柴达木盆地开展MT实测资料激发极化信息提取研究, 试验结果表明, 气田呈高电阻率、高极化率异常模式, 该成果对研究区含油气性检测评价起到了指示性作用。董莉等^[29]通过改进自适应策略, 提出一种基于非均匀统计分布的自适应差分进化两阶段最小构造反演方法, 提高了对微弱激电信息提取的精度, 降低了MT资料解释的多解性。截止目前, MT激电勘探主要应用于矿产和油气勘探领域, 针对深层地热资源勘探的MT方法仅限于常规方法, 为提高MT方法在深层地热勘探中的应用效果, 有必要对激电条件下MT勘探在深层热储中的响应特征进行深入分析。

本文采用有限差分数值计算方法, 通过引入Cole-Cole复电阻率模型, 完成了不同激电模型的二维MT正演计算。通过分析不同激电参数对深层热储视电阻率与相位响应特征, 结合热储岩石温度与激电参数关系, 为深层地热电磁勘探提供新方法。

1 正演理论 1.1 电磁场基本理论

当平面入射的电磁波在地球介质中传播时, 假定近地表大地电磁场为稳态场, 取时谐因子为e^－iwt, 基于本构方程, 频率域Maxwell方程组的旋度方程可以表示为:

$ \left\{\begin{array}{l} \nabla \times \boldsymbol{E}=\mathrm{i} \omega \mu \boldsymbol{H} \\ \nabla \times \boldsymbol{H}=(\sigma-\mathrm{i} \omega \varepsilon) \boldsymbol{E} \end{array}\right. $

(1)

式中: E, H分别为电场强度和磁场强度; ω为角频率; μ为磁导率; ε为介电常数; σ为电导率。

设定二维模型构造走向沿x方向, 倾向沿y方向, 即∂/∂x=0, 在笛卡尔坐标系中, 旋度方程能解耦得到如公式(2)所示的TE、TM极化模式电磁场方程组。

$ \left\{\begin{array}{l} \frac{\partial^{2} E_{x}}{\partial y^{2}}+\frac{\partial^{2} E_{x}}{\partial z^{2}}+\mathrm{i} \omega \mu(\sigma-\mathrm{i} \omega \varepsilon) E_{x}=0 \\ \frac{1}{\sigma-\mathrm{i} \omega \varepsilon} \cdot\left(\frac{\partial^{2} H_{x}}{\partial y^{2}}+\frac{\partial^{2} H_{x}}{\partial z^{2}}\right)+\mathrm{i} \omega \mu H_{x}=0 \end{array}\right. $

(2)

式中: E_x, H_x分别为电场强度和磁场强度沿x方向的分量。

1.2 有限差分数值模拟

本文采用有限差分数值计算方法进行大地电磁二维正演计算, 利用差分方法可以将公式(2)中的亥姆赫兹方程组分别转化为TE和TM极化模式的二阶偏微分方程, 因此, 不需要采用交错网格, 而是分别求取E_x, H_x对应的二阶偏微分方程。构建如图 1所示的矩形剖分网格, 其中, y方向的节点编号用j表示, z方向的节点编号用k表示。

在忽略位移电流项的条件下, 采用三点差分公式, 将公式(2)中的微分方程组改写为差分方程组, 则获得节点P(j, k)处的差分公式为:

$ \begin{aligned} &\frac{2}{\left(\mathrm{~d} z_{k}+\mathrm{d} z_{k-1}\right) \mathrm{d} z_{k-1}} E_{x}(j, k-1)+ \\ &\qquad \frac{2}{\left(\mathrm{~d} y_{j}+\mathrm{d} y_{j-1}\right) \mathrm{d} y_{j-1}} E_{x}(j-1, k)+ \\ &\qquad {\left[\mathrm{i} \omega \mu \bar{\sigma}(j, k)-\frac{2}{\mathrm{~d} y_{j} \mathrm{~d} y_{j-1}}-\frac{2}{\mathrm{~d} z_{k} \mathrm{~d} z_{k-1}}\right]}\cdot \\ &\qquad E_{x}(j, k)+\frac{2}{\left(\mathrm{~d} y_{j}+\mathrm{d} y_{j-1}\right) \mathrm{d} y_{j}} E_{x}(j+1, k)+ \\ &\qquad \frac{2}{\left(\mathrm{~d} z_{k}+\mathrm{d} z_{k-1}\right) \mathrm{d} z_{k}} E_{x}(j, k+1)=0 \end{aligned} $

(3)

$ \begin{aligned} &\frac{2 \bar{\rho}_{z}(j, k-1)}{\left(\mathrm{d} z_{k}+\mathrm{d} z_{k-1}\right) \mathrm{d} z_{k-1}} H_{x}(j, k-1)+ \\ &\qquad \frac{2 \bar{\rho}_{y}(j-1, k)}{\left(\mathrm{d} y_{j}+\mathrm{d} y_{j-1}\right) \mathrm{d} y_{j-1}} H_{x}(j-1, k)+ \\ &\qquad {\left[\mathrm{i} \omega \mu-\frac{2 \bar{\rho}_{y}(j-1, k)}{\left(\mathrm{d} y_{j}+\mathrm{d} y_{j-1}\right) \mathrm{d} y_{j-1}}-\frac{2 \bar{\rho}_{y}(j, k)}{\left(\mathrm{d} y_{j}+\mathrm{d} y_{j-1}\right) \mathrm{d} y_{j}}-\right.} \\ &\qquad \left.\frac{2 \bar{\rho}_{z}(j, k-1)}{\left(\mathrm{d} z_{k}+\mathrm{d} z_{k-1}\right) \mathrm{d} z_{k-1}}-\frac{2 \bar{\rho}_{z}(j, k)}{\left(\mathrm{d} z_{k}+\mathrm{d} z_{k-1}\right) \mathrm{d} z_{k}}\right]\cdot \\ &\qquad H_{x}(j, k)+\frac{2 \bar{\rho}_{y}(j, k)}{\left(\mathrm{d} y_{j}+\mathrm{d} y_{j-1}\right) \mathrm{d} y_{j}} H_{x}(j+1, k)+ \\ &\qquad \frac{2 \bar{\rho}_{z}(j, k)}{\left(\mathrm{d} z_{k}+\mathrm{d} z_{k-1}\right) \mathrm{d} z_{k}} H_{x}(j, k+1)=0 \end{aligned} $

(4)

式中: σ(j, k)为节点P(j, k)处的平均电导率; ρ_y(j, k), ρ_z(j, k)分别为节点P(j, k)在y方向和z方向的平均电阻率。其计算公式为:

$ \begin{aligned} \bar{\sigma}&(j, k)=\left(\mathrm{d} y_{j}+\mathrm{d} y_{j-1}\right)\left(\mathrm{d} z_{k}+\mathrm{d} z_{k-1}\right) /\left[\rho(j, k) \mathrm{d} y_{j}\cdot\right. \\ &\mathrm{d} z_{k}+\rho(j, k-1) \mathrm{d} y_{j} \mathrm{~d} z_{k-1}+\rho(j-1, k)\cdot \\ &\left.\mathrm{d} y_{j-1} \mathrm{~d} z_{k}+\rho(j-1, k-1) \mathrm{d} y_{j-1} \mathrm{~d} z_{k-1}\right] \end{aligned} $

(5a)

$ \begin{aligned} \bar{\rho}_{y}&(j, k)= {\left[\rho(j, k) \mathrm{d} y_{j} \mathrm{~d} z_{k}+\rho(j, k-1) \mathrm{d} y_{j} \mathrm{~d} z_{k-1}\right] /} \\ & \left(\mathrm{d} y_{j} \mathrm{~d} z_{k}+\mathrm{d} y_{j} \mathrm{~d} z_{k-1}\right) \end{aligned} $

(5b)

$ \begin{aligned} &\bar{\rho}_{z}(j, k)=\left[\rho(j, k) \mathrm{d} y_{j} \mathrm{~d} z_{k}+\rho(j-1, k) \mathrm{d} y_{j-1} \mathrm{~d} z_{k}\right] / \\ &\quad\left(\mathrm{d} y_{j} \mathrm{~d} z_{k}+\mathrm{d} y_{j-1} \mathrm{~d} z_{k}\right) \end{aligned} $

(5c)

本文进行正演计算时统一采用第一类边界条件, 4个边界都距离探测目标区足够远, 不受地下不均匀体的影响, 上边界在距离地表足够远(一般大于100 km)的高空中, 取场值为常数(通常取1), 下边界在地下足够深的地方(探测深度的5倍以上), 取值为0, 左、右边界条件分别为一维正演结果, 则本文差分方程求解的边界条件为:

$ \left\{\begin{array}{l} E_{x}(j, 0)=1, H_{x}(j, 0)=1 \\ E_{x}(j, \infty)=0, H_{x}(j, \infty)=0 \\ E_{x}(1, k)=E_{x \mathrm{~L}}(k), H_{x}(1, k)=H_{x \mathrm{L}}(k) \\ E_{x}\left(N_{y}+1, k\right)=E_{x \mathrm{R}}(k), H_{x}\left(N_{y}+1, k\right)=H_{x \mathrm{R}}(k) \end{array}\right. $

(6)

式中: E_xL(k), H_xL(k), E_xR(k), H_xR(k)分别为一维正演得到的左、右边界电、磁场值; N_y为沿y方向的第N个节点编号。

联立公式(3)、公式(4)和公式(6), 则可得到大型线性方程组:

$ \boldsymbol{K} \boldsymbol{x}=\boldsymbol{b} $

(7)

式中: K为大型稀疏矩阵; x为待求解的电磁场值; b为矢量边界条件。

求解公式(7)可得到电磁场分量, 结合公式(8)所示的卡尼亚视电阻率公式, 即可计算出地表观测点对应的视电阻率与阻抗相位。

$ \left\{\begin{array}{l} \rho_{x y}=\frac{1}{\omega \mu}\left|\frac{E_{x}}{H_{y}}\right|^{2}, \varphi_{x y}=\arctan \frac{\operatorname{Im}\left(\frac{E_{x}}{H_{y}}\right)}{\operatorname{Re}\left(\frac{E_{x}}{H_{y}}\right)} \\ \rho_{y x}=\frac{1}{\omega \mu}\left|\frac{E_{y}}{H_{x}}\right|^{2}, \varphi_{y x}=\arctan \frac{\operatorname{Im}\left(\frac{E_{y}}{H_{x}}\right)}{\operatorname{Re}\left(\frac{E_{y}}{H_{x}}\right)} \end{array}\right. $

(8)

1.3 基于Cole-Cole模型的二维MT正演

本文采用COLE等^[30]提出的Cole-Cole复电阻率模型, 开展基于激电效应的二维MT正演研究, 其中Cole-Cole复电阻率数学表达式如下:

$ \rho(\mathrm{i} \omega)=\frac{\rho}{1-m}\left\{1-m\left[1-\frac{1}{1+(\mathrm{i} \omega \tau)^{c}}\right]\right\} $

(9)

式中: ρ(iω)为不同频率条件下岩石的复电阻率; ρ为不存在极化时的直流电阻率; m为极化率; τ为时间常数; c为频率相关系数。

基于激电模型的二维MT正演计算时, 由于激电效应的存在, 网格单元的电阻率应为复电阻率。因此, 将公式(3)、公式(4)中的平均电导率和平均电阻率, 替换成公式(9)中的平均复电导率和平均复电阻率。通过改变激电参数即可实现不同激电模型的二维MT正演。

2 深层热储激电响应特征分析 2.1 基于激电参数的深层热储MT勘探可行性分析

在针对岩石电阻率与温度关系的研究中, 国内外很多学者认为: 以离子导电方式为主的地下岩石, 随着温度的升高, 离子活性增强, 岩石内部电离程度增加, 使得导电离子总数增多, 从而电阻率降低。大量实验研究发现, 岩石电阻率与温度的经验公式^[31-33]为:

$ \rho=\frac{\rho_{i}}{1+\alpha\left(T-T_{i}\right)} $

(10)

式中: ρ为温度T时的电阻率; ρ_i为温度T_i时的电阻率; α为温度系数, 其值与其岩性和地下水溶液矿化度相关, 一般取0.02。

将公式(10)变型, 即可得到电阻率变化率与温度变化的关系式:

$ \varepsilon_{\rho}=\left|\frac{\rho-\rho_{i}}{\rho_{i}}\right|=\left|\frac{-\alpha \Delta T}{1+\alpha \Delta T}\right| $

(11)

式中: ε_ρ为电阻率相对变化率; ΔT=T－T_i为温度差。

计算不同温度差条件下岩石电阻率相对变化率, 结果如表 1所示。从表 1中可以看出, 当温度差达到20℃时, 能够引起28.57%的电阻率相对变化率, 说明电磁方法进行地热勘探具有可行性。

由于地下岩石具有频散效应, 其电阻率是一个随频率变化的值, 即复电阻率。本文引入Cole-Cole复电阻率模型来表征岩石复电阻率, 联立公式(9)和公式(10)得到温度与岩石激电参数的关系式:

$ \begin{gathered} T=T_{i}+\frac{1}{\alpha}\left\{\frac { \rho _ { i } } { {\bar{\rho}} ( 1 - m ) ^ { - 1 } } \cdot \left\{1-m \cdot\right.\right. \\ \left.\left.\left[1-\frac{1}{1+(\mathrm{i} \omega \tau)^{c}}\right]\right\}^{-1}-1\right\} \end{gathered} $

(12)

公式(12)说明, 电阻率、极化率、时间常数及频率相关系数能够反映温度变化, 为基于激电效应的MT方法进行深层地热勘探提供了理论依据。

2.2 算法验证

本文采用3层层状介质模型(K型)进行二维MT有限差分正演算法验证, 层状介质模型参数如表 2所示。

对比二维MT有限差分正演结果与一维MT解析解结果发现, 两种方法的计算结果基本一致(图 2), 从而验证了二维MT有限差分方法的正确性。

2.3 深层热储极化异常响应特征研究

为研究深层地热热储MT激电响应特征, 本文设置如图 3、表 3所示的激电模型。如图 3所示, 在电阻率为100 Ω·m的均匀半空间中, 存在一个电阻率为10 Ω·m, 埋深为3 km, 大小为1 km×2 km的热储极化异常体, 其激电参数如表 3所示。本文正演计算频率范围为10^-3~10³ Hz, 按对数等间距取61个频点, 地表测点范围为-5~5 km, 每0.1 km一个测点, 共101个测点。

本文设置3种不同激电参数模型进行深层热储极化异常体响应特征研究。

2.3.1 不同极化率的MT响应

如表 3中异常体模型一所示, 在100 Ω·m的均匀半空间中, 存在一个电阻率为10 Ω·m、时间常数为10^-1、频率相关系数为0.4、极化率分别为0, 0.2, 0.4, 0.6, 0.8等5种情况的极化异常体。基于Cole-Cole复电阻率模型, 开展不同极化率条件下的二维MT正演, 得到如图 4所示的视电阻率和相位拟断面图。

图 4a至图 4d自上而下分别显示了极化率为0, 0.2, 0.4, 0.6, 0.8时, TE、TM极化模型下二维MT正演视电阻率和相位响应。由图 4可知, 地下介质的视电阻率和相位随极化率参数的增加而变大, 其中TM极化模式响应明显强于TE极化模式, 且TM极化模式的视电阻率响应大于相位响应。基于此, 开展TM极化模式下视电阻率相对变化率计算, 结果如图 5所示, 研究发现, 当极化率差大于0.2时, 视电阻率相对变化率最大值高于20%, 表明因极化率参数改变而产生的激电响应较强, 说明该参数可作为深层地热勘探的有效参数。

2.3.2 不同时间常数的MT响应

设置如图 3、表 3中异常体模型二所示的二维MT正演模型, 在100 Ω·m的均匀半空间中, 存在一个电阻率为10 Ω·m、极化率为0.4、频率相关系数为0.4、时间常数分别为10^-2, 10^-1, 10⁰, 10¹, 10²等5种情况的极化异常体。基于Cole-Cole复电阻率模型, 采用有限差分算法, 开展不同时间常数条件下二维MT正演, 得到如图 6所示的视电阻率和相位拟断面图。

图 6清晰反映了时间常数对大地电磁视电阻率和相位有较大影响, 随着时间常数的增大, 视电阻率和相位有明显的减小趋势, 其中视电阻率的响应大于相位的响应, 且TM极化模式对时间常数的变化相对敏感。计算TM极化模式下视电阻率相对变化率, 结果如图 7所示, 可以发现, 当时间常数差大于(10²-10¹)时, 视电阻率相对变化率最大值大于20%, 表明时间常数改变引起的激电响应较弱。

2.3.3 不同频率相关系数的MT响应

在笛卡尔直角坐标系下, 设计100 Ω·m的均匀半空间中, 存在一个电阻率为10 Ω·m、极化率为0.4、时间常数为10^-1、频率相关系数分别为0.2, 0.4, 0.6, 0.8等4种情况的极化异常体, 如图 3、表 3中异常体模型三所示。基于MT正演理论, 采用有限差分算法, 通过引入Cole-Cole复电阻率模型, 得到不同频率相关系数条件下, 二维MT正演视电阻率和相位拟断面图, 结果如图 8所示。

由图 8可知, 大地电磁视电阻率和相位随着频率相关系数的增大而增大, 其中TM极化模式响应特征比TE极化模式的响应特征明显, 且TM极化模式条件下, 视电阻率的变化更敏感。计算TM极化模式下视电阻率相对变化率, 结果如图 9所示, 可以发现, 频率相关系数差大于0.6时, 视电阻率相对变化率最大值为20%, 表明频率相关系数改变引起的激电响应也相对较弱。

3 结论

基于有限差分方法, 同时引入Cole-Cole模型, 实现了二维激电介质的MT正演。通过设置不同激电参数(极化率、时间常数、频率相关系数)模型, 模拟计算并分析深层热储模型的MT视电阻率和相位响应特征, 得到如下结论与建议:

1) 研究不同激电模型二维MT响应发现, 视电阻率和相位随极化率、频率相关系数的增大而增大, 随时间常数的增大而减小;

2) 二维极化介质条件下, TM极化模式的响应特征明显强于TE极化模式的响应特征, 且视电阻率与相位相比, 对激电参数的变化更敏感;

3) 计算TM极化模式下视电阻率相对变化率, 结合温度与电阻率变化率理论关系, 发现由极化率参数引起的激电响应较强, 而由时间常数和频率相关系数引起的激电响应相对较弱。

综上所述, 在实际深层地热二维MT勘探时, 为提高资料处理解释的精度与速度, 建议采用TM极化模式进行资料反演, 并采用电阻率和极化率多参数联合解释方式, 以提高MT方法在深层地热勘探中的应用能力。