水平电偶极源层状模型垂直磁场全区视电阻率计算方法

图1 水平电偶极源层状地电模型

Fig.1 Layered geoelectricity model with horizontal electric dipole source

略去繁冗的数学推导^{[1-2,7,9,14-15]},此处直接写出极坐标系统中地表电磁场各分量的表达式为:

(1a)

$E_{r} = \frac{I d L}{2 π} \cos ϕ [\frac{i ωμ}{r} \int_{0}^{\infty} \frac{1}{λ + \frac{m_{1}}{R^{*}}} J_{1} (λr) d λ - ρ_{1} \int_{0}^{\infty} \frac{λ m_{1}}{R} J_{0} (λr) d λ + \frac{ρ_{1}}{r} \int_{0}^{\infty} \frac{m_{1}}{R} J_{1} (λr) d λ],$

(1b)

$E_{ϕ} = \frac{I d L}{2 π} \sin ϕ [\frac{ρ_{1}}{r} \int_{0}^{\infty} \frac{m_{1}}{R} J_{1} (λr) d λ - i ωμ \int_{0}^{\infty} \frac{λ}{λ + \frac{m_{1}}{R^{*}}} J_{0} (λr) d λ + \frac{i ωμ}{r} \int_{0}^{\infty} \frac{1}{λ + \frac{m_{1}}{R^{*}}} J_{1} (λr) d λ],$

(1c)

$H_{r} = - \frac{I d L}{2 π r} \sin ϕ [\int_{0}^{\infty} \frac{λ}{λ + \frac{m_{1}}{R^{*}}} J_{1} (λr) d λ + r \int_{0}^{\infty} \frac{m_{1}}{R^{*}} \frac{λ}{λ + \frac{m_{1}}{R^{*}}} J_{0} (λr) d λ],$

(1d)

$H_{ϕ} = \frac{I d L}{2 π r} \cos ϕ [\int_{0}^{\infty} \frac{λ}{λ + \frac{m_{1}}{R^{*}}} J_{1} (λr) d λ],$

(1e)

$H_{z} = \frac{I d L}{2 π} \sin ϕ [\int_{0}^{\infty} \frac{λ^{2}}{λ + \frac{m_{1}}{R^{*}}} J_{1} (λr) d λ] 。$

式中:

$\begin{array}{l} R^{*} = \coth [m_{1} d_{1} + \cot h^{- 1} \frac{m_{1}}{m_{2}} \coth (m_{2} d_{2} + L + \cot h^{- 1} \frac{m_{N - 1}}{m_{N}})]; \\ R = \coth [m_{1} d_{1} + \cot h^{- 1} (\frac{m_{1}}{m_{2}} \frac{ρ_{1}}{ρ_{2}}) \coth (m_{2} d_{2} + L + \cot h^{- 1} (\frac{m_{N - 1}}{m_{N}} \frac{ρ_{N - 1}}{ρ_{N}}))]; \end{array}$

m_j= $\sqrt[]{λ^{2} - k_{j}^{2}}$ ; $k_{j}^{2}$ =iωμ₀σ_n,μ₀为自由空间导磁率;λ称为空间频率,它具有距离倒数的量纲。

特别地,当N=1时,可得到均匀半空间表面电磁场各分量的表达式:

(2a)

$E_{r 0} = \frac{I d L}{2 π σ r^{3}} \cos ϕ [1 + e^{i kr} (1 - i kr)],$

(2b)

$E_{ϕ 0} = \frac{I d L}{2 π σ r^{3}} \sin ϕ [2 - e^{i kr} (1 - i kr)],$

(2c)

$\begin{array}{l} H_{r 0} = \frac{- 3 I d L}{2 π r^{2}} \sin ϕ \{I_{1} (\frac{i kr}{2}) K_{1} (\frac{i kr}{2}) + \\ \frac{i kr}{6} [I_{1} (\frac{i kr}{2}) K_{0} (\frac{i kr}{2}) - I_{0} (\frac{i kr}{2}) K_{1} (\frac{i kr}{2})]\}, \end{array}$

(2d)

$H_{ϕ 0} = \frac{I d L}{2 π r^{2}} \cos ϕ I_{1} (\frac{i kr}{2}) K_{1} (\frac{i kr}{2}),$

(2e)

$H_{z 0} = - \frac{3 I d L}{2 π k^{2} r^{4}} \sin ϕ [1 - e^{i kr} (1 - i kr - \frac{1}{3} k^{2} r^{2})] 。$

以上各式中,μ₀为自由空间的导磁率。ω为谐变电流的圆频率。I₀ $(\frac{i kr}{2})$ 、I₁ $(\frac{i kr}{2})$ 和K₀ $(\frac{i kr}{2})$ 、K₁ $(\frac{i kr}{2})$ 分别为以 $\frac{i kr}{2}$ 为宗量的第一和第二类虚宗量贝寒尔函数,下标“0”或“1”表示贝寒尔函数的阶数,k为波数。

1.1 电磁场计算方法

令:水平电偶极矩P_E= $\frac{I d L}{2 π}$ ,式(1)可写为:

(3a)

$E_{r} = (P_{E} / r) \cos ϕ [i ωμ I_{4} - r ρ_{1} I_{5} + ρ_{1} I_{6}] + E_{r 0},$

(3b)

$E_{ϕ} = (P_{E} / r) \sin ϕ [i ωμ I_{4} - i ωμr I_{7} + ρ_{1} I_{6}] + E_{ϕ 0},$

(3c)

$H_{r} = (- P_{E} / r) \sin ϕ (I_{1} + 0.5 I_{9} + r (I_{2} + 0.5 I_{10}))$

(3d)

$H_{ϕ} = \frac{P_{E} \cos ϕ}{r} (I_{1} + 0.5 I_{9}),$

(3e)

$H_{z} = P_{E} \sin ϕ (I_{3} + 0.5 I_{8}) 。$

式中:

(4a)

$I_{1} = \int_{0}^{\infty} (\frac{λ}{λ + \frac{m_{1}}{R^{*}}} - 0.5) J_{1} (λr) d λ,$

(4b)

$I_{2} = \int_{0}^{\infty} (\frac{m_{1}}{R^{*}} \frac{1}{λ + \frac{m_{1}}{R^{*}}} - 0.5) λ J_{0} (λr) d λ,$

(4c)

$I_{3} = \int_{0}^{\infty} (\frac{λ}{λ + m_{1} / R^{*}} - 0.5) λ J_{1} (λr) d λ,$

(4d)

$I_{4} = \int_{0}^{\infty} (\frac{λ}{λ + m_{1} / R^{*}} - \frac{1}{λ + m_{1}}) J_{1} (λr) d λ,$

(4e)

$I_{5} = \int_{0}^{\infty} λ m_{1} (\frac{1}{R} - 1) J_{0} (λr) d λ,$

(4f)

$I_{6} = \int_{0}^{\infty} m_{1} (\frac{1}{R} - 1) J_{1} (λr) d λ$

(4g)

$I_{7} = \int_{0}^{\infty} (\frac{1}{λ + m_{1} / R^{*}} - \frac{1}{λ + m_{1}}) λ J_{0} (λr) d λ,$

(4h)

$I_{8} = \int_{0}^{\infty} λ J_{1} (λr) d λ,$

(4i)

$I_{9} = \int_{0}^{\infty} J_{1} (λr) d λ = \frac{1}{r},$

(4j)

$I_{10} = \int_{0}^{\infty} λ J_{0} (λr) d λ = 0,$

(4k)

$E_{r 0} = \frac{P_{E}}{r^{3}} ρ \cos ϕ [1 + e^{i kr} (1 - i kr)],$

(4l)

$E_{ϕ 0} = \frac{P_{E}}{r^{3}} ρ \sin ϕ [2 - e^{i kr} (1 - i kr)] 。$

在计算层状介质表面的电场分量和磁场垂直分量时,减去了均匀半空间的场,这是为了保证在λ→∞时积分核快速趋于零,从而保证积分收敛,也节省计算时间。

写出式(1c)、(1d)、(1e)中的核函数:

(5a)

$K_{4} (λ) = \frac{λ}{λ + m_{1} / R^{*}},$

(5b)

$K_{5} (λ) = \frac{m_{1}}{R^{*}} \frac{λ}{λ + m_{1} / R^{*}},$

(5c)

$K_{3} (λ) = \frac{m_{1}}{R^{*}} \frac{λ^{2}}{λ + m_{1} / R^{*}} 。$

比较发现上述三式有类似的特点,对其取极限,则有:

(6a)

$\lim_{λ \to \infty} K_{4} (λ) = 0.5,$

(6b)

$\lim_{λ \to \infty} K_{5} (λ) = 0.5 λ,$

(6c)

$\lim_{λ \to \infty} K_{3} (λ) = 0.5 λ^{2} 。$

这三个极限的收敛速度是很快的。据此,将式(1)各式改写为式(3)。

当采用直角坐标时,

(7a)

$E_{x} = E_{r} \cos ϕ - E_{ϕ} \sin ϕ,$

$H_{y} = H_{r} \sin ϕ + H_{ϕ} \cos ϕ 。$ (7b)

层状地质模型的频率域电磁响应(式3,式4)是以Hankel积分形式给出的,而目前的计算方法主要是数字滤波方法^[5-6,9]。但数字滤波本身计算精度不够高,且滤波系数计算精度取决于具体的计算方法,只有在合适的条件下才有足够的精度。本文采用直接数值积分,将水平电偶源电磁场频率响应所满足的Hankel积分用由Bessel函数零点所分割的子区间上的部分积分的和表示出来,并采用Euler方法加快积分求和的收敛^[20,21]。

文献[20]中,Chave提出了汉克尔变换的数值积分算法。Chave算法的基本思路是,采用高斯求积方法对有限区间积分进行计算,进而结合连分式展开方法加快部分积分和的收敛速度。一反数字滤波法的常规做法,Chave算法使用期望误差(绝对误差/相对误差)来控制计算是继续或是终止。为了直接采用数值求积方法计算式(4),可将积分写成如下部分积分之和:

(8)

f \approx S = \overset{N}{\sum_{n = 0}} P_{n},

(9)

P_{n} = \int_{Z_{n}}^{Z_{n + 1}} K (λ, P) J_{i} (λr) d λ 。

式中:J_i(λr)为第1类的i阶Bessel函数;Z_n为 J_i(λr)的第n个零点用距离r归一化后的值;N为参与求和的部分积分项的个数,取决于计算方法;P_n可采用Gauss求积公式计算

(10)

P_{n} \approx \overset{M}{\sum_{j = 1}} h_{j} K (a_{j}, P) J_{i} (a_{j} r) 。

式(10)中M是求积节点个数,a_j为求积节点,h_j为求积系数。

采用连分式加速式(10)积分项求和的收敛

(11)

S = \frac{d_{0}}{1 + \frac{d_{1}}{1 + \frac{d_{2}}{1 + \frac{d_{3}}{1 + \frac{\dots}{d_{n}}}}}},

式中d₀,d₁,d₂,…,d_N可由区间积分项P_n求得。

本文采用文献[21]介绍的数值积分方法,在matlab计算平台中,预先计算了贝赛尔函数J_i(λr),i=0,1的300个零点,最初的几个零点之间的区间积分采用matlab 的自适应Simpsons积分函数(QUAD),之后的积分区间根据计算精度,采用自适应Lobatto积分函数(QUADL),积分项的求和采用Euler法加速其收敛,积分的终止条件是最后加入的区间积分项小于设定的容许误差,或者已计算了设定的所有区间数。实际上计算不超过30个过零点区间就可达到所需要的精度。

由上,可计算出水平电偶极源层状模型电场和磁场各分量频率响应。

1.2 全区视电阻率计算方法

由计算出的水平电偶极源层状模型的各场分量的频率响应,汤井田给出了各场分量全区视电阻率定义^[14]。本文主要以水平电场E_x和垂直磁场H_z分量计算全区视电阻率。为方便论述,此处重写出均匀半空间表面水平电偶源产生的电磁场:

12a

$E_{x} = \frac{2 I d L}{2 π r^{3}} ρ [3 \cos ϕ - 2 + e^{i kr} (1 - i kr)]$

(12b)

$\begin{array}{l} H_{y} = - \frac{I d L}{4 π r^{2}} [si n^{2} ϕ (6 I_{1} K_{1} + i kr (I_{1} K_{1} - I_{0} K_{1})) \\ - 2 co s^{2} ϕ I_{1} K_{1}], \end{array}$

(12c)

$H_{z} = - \frac{I d Lρ}{2 π r^{4} μω} \sin ϕ [3 - e^{i kr} (3 - 3 i kr - k^{2} r^{2})] 。$

式中:E_x、H_y和H_z分别为电场和磁场的水平分量和垂直分量,r为接收点到偶极中心矢径的模,φ为r和x轴的夹角。I是电流强度,dL是偶极子长度。ρ和μ是均匀半空间的电阻率和导磁率,ω为角频率,k为电磁波的波数,k²=iωμ/ρ,i= $\sqrt[]{- 1}$ 。I₀和K₁、K₀分别是第一和第二类以ikr/2为宗量的虚宗量贝塞尔函数,0和1表示阶数。

若令:

(13a)

$C_{E} = 3 co s^{2} ϕ - 2 + e^{i kr} (1 - i kr),$

(13b)

$\begin{array}{l} C_{H} = [si n^{2} ϕ (6 I_{1} K_{1} + i kr (I_{1} K_{0} - I_{0} K_{1})) \\ - 2 co s^{2} ϕ I_{1} K_{1}] kr, \end{array}$

(13c)

$C_{Z} = 3 - e^{i kr} (3 - 3 i kr - k^{2} r^{2});$

则有:

(14a)

$E_{x} = \frac{I d Lρ}{2 π r^{3}} C_{E},$

(14b)

$H_{z} = - \frac{I d Lρ}{2 π r^{4} μω} \sin ϕ \cdot C_{Z} 。$

定义

(15a)

$ρ_{e}^{| E_{x} |} = \frac{2 π r_{3}}{I d L} |\frac{E_{x}}{C_{E}}|,$

(15b)

$Φ_{e}^{| E_{x} |} = \arctan [\frac{Im (E_{x} / C_{E})}{Re (E_{x} / C_{E})}];$

(16a)

$ρ_{e}^{| H_{z} |} = \frac{2 π r^{2}}{I d L \sin ϕ} μω |\frac{H_{z}}{C_{Z}}| (φ \neq 0^{0}),$

(16b)

$Φ_{e}^{| H_{z} |} = \arctan [\frac{Im (H_{z} / C_{Z})}{Re (H_{z} / C_{Z})}] 。$

本文给出采用循环互相关算法计算全区视电阻率的方法^[18]。观测中需要同时记录发送电流和电场或磁场时间序列。当观测场值是水平电场E_x分量时间序列时,观测值是接收电极间的电位差时间序列, $V_{E_{x}}$ =E_x·MN,MN为接收电场的电极间距,采用循环互相关算法可估计出阻抗谱:

(17)

Z_{E_{x}} (ω) = \frac{P_{C V_{E_{x}} s} (ω)}{P_{C I_{AB} s} (ω)},

式中: $P_{C V_{E_{x}} s}$ (ω)为参考信号SS与接收电位差信号 $V_{E_{x}}$ 的互功率谱, $P_{C I_{AB} s}$ (ω)为参考信号SS与水平电偶源AB发射电流信号I_AB(t)的互功率谱。

由均匀大地表面水平电偶极子源产生的电场E_x计算公式:

(18)

E_{x} = \frac{ρI d L}{2 π r^{3}} [(1 - 3 si n^{2} (φ)) + (1 - i k_{1} r) \exp (i k_{1} r)] 。

式中:ρ为均匀大地的电阻率,I是发射电流,r是接收电极MN中点至发射偶极AB中点的距离,φ是接收电极MN中点至发射AB极中点连线与AB极连线的夹角。dL是发射AB极间距。 $k_{1}^{2}$ =iωμ₀/ρ为均匀大地的波数。

式(18)两边除以电流I与式(17)对应相等,可得:

(19)

Z_{E_{x}} (ω) = \frac{ρ (ω) d L \cdot MN}{2 π r^{3}} C_{E_{x}}, (19)

(20)

C_{E_{x}} = (1 - 3 si n^{2} (φ)) + (1 - i k_{1} r) \exp (i k_{1} r) 。

$\begin{array}{l} \end{array}$

$Z_{E_{x}}$ (ω)就是由观测电场E_x电位差值和激发电流I循环互相关辨识得到的大地的阻抗。由此可得大地的全区视电阻率谱:

(21)

式中: $C_{E_{x}}$ 也是大地电阻率的函数,因此,式(21)全区视电阻率需要采用迭代法求解。这个算法已在实测数据处理中得到验证。

当观测电磁场值是垂直磁场H_z分量时间序列时,由发送电流时间序列和垂直磁场时间序列,采用循环互相关算法可估计出阻抗谱:

(22)

Z_{H_{z}} (ω) = \frac{P_{C H_{z} s} (ω)}{P_{C I_{AB} s} (ω)},

式中: $P_{C H_{z} s}$ (ω)为参考信号SS与接收垂直磁场信号H_z的互功率谱, $P_{C I_{AB} s}$ (ω)为参考信号SS与水平电偶极子源AB发射电流信号I_AB(t)的互功率谱。

由均匀大地表面水平电偶极子源产生的垂直磁场H_z计算公式为:

(23)

H_{Z} = - \frac{3 I d L}{2 π k_{1}^{2} r^{4}} \sin φ [1 - (1 - i k_{1} r - \frac{1}{3} (k_{1} {r)}^{2}) \exp (i k_{1} r)],

式中:I是发射电流,r是磁场接收点P至发射极AB中点的距离,φ是磁场接收点P至发射极AB中点连线与AB极延长线的夹角,dL是发射极AB长度, $k_{1}^{2}$ =iωμ₀/ρ为均匀大地的波数,μ₀是磁导率。

式(23)两边除以电流I与式(22)对应相等,可得:

(24)

Z_{H_{z}} (ω) = - \frac{3 d L \cdot ρ (ω)}{i ω μ_{0} 2 π r^{4}} \sin φ C_{H_{z}}

(25)

\begin{array}{l} C_{H_{z}} = 1 - (1 - i k_{1} r - \frac{1}{3} (k_{1} {r)}^{2}) \exp (i k_{1} r) 。 \end{array}

$Z_{H_{z}}$ (ω)就是由观测垂直磁场H_z和发送电流I循环互相关辨识得到的大地的频率域响应。由此可得大地的全区视电阻率谱:

(26)

ρ_{H_{z}} (ω) = \frac{2 π ω μ_{0} r^{4} \cdot Z_{H_{z}} (ω)}{3 i \cdot d L \cdot \sin φ \cdot C_{H_{z}}} 。

式中 $C_{H_{z}}$ 也是大地电阻率的函数,因此,式(26)需要采用迭代法求解。

对于激发极化效应的模拟计算,采用在频率域Cole-Cole模型来描述:

(27)

ρ (ω) = ρ_{0} [1 - m (1 - \frac{1}{1 + {(i ωτ)}^{c}})],

式中:ρ₀为直流电阻率,m为极化率,ω为角频率,τ为时间常数,c为频率系数。

式(21)和式(26)右端均隐含有视电阻率,需要采用迭代算法求解。全区视电阻率迭代计算过程如下:

1) 给出全区视电阻率初值ρ⁽⁰⁾。此初值是任意的,可选取为相应的波区视电阻率;

2) 将ρ⁽⁰⁾代入式(21)和式(26),并进而求得第一次迭代视电阻率值ρ⁽¹⁾;

3) 判断:

(28)

\frac{| ρ^{(1)} - ρ^{(0)} |}{ρ^{(0)}} \leq ε;

式中ε为给定的正的小数。若式(28)成立,计算相应的相位值,并停止迭代。若式(28)不成立,令ρ⁽⁰⁾=ρ⁽¹⁾,返回步骤(2),直到式(28)成立。

2 模型计算

2.1 模型一——D型

地电模型参数如表1所示。计算参数:供电电流I=40 A,电偶极子长度dL=1 800 m,收发距R=12 000 m,赤道偶极测量装置。计算了水平电场E_x全区视电阻率、垂直磁场H_z全区视电阻率,同时与MT卡尼亚视电阻率和CSAMT卡尼亚视电阻率进行比较。从图2的频谱曲线可以看出,大于1 000 Hz高频段全区视电阻率反映出第一层电阻率300 Ω·m,低于0.01 Hz低频段E_x全区视电阻率和垂直磁场H_z全区视电阻率,与MT卡尼亚视电阻率趋势一致,趋于第二层视电阻率20 Ω·m。而CSAMT卡尼亚视电阻率在低于1 Hz低频段不能正确反映地层特征。

表1 D型地电模型参数

Table 1 D type geoelectricity model parameters

层参数	第一层	第二层
电阻率/(Ω·m)	300	20
层厚/m	1100	Inf
m	0	0
τ/s	0	0
c	0	0

图2

图2 水平电偶源D型地电模型全区视电阻率频谱

Fig.2 whole zone apparent resistivity spectrum of D-type geoelectricity model with horizontal electric dipole Source

2.2 模型二——G型

地电模型参数如表2所示。

表2 G型地电模型参数

Table 2 G type geoelectricity model parameters

层参数	第一层	第二层
电阻率/(Ω·m)	300	1000
层厚/m	1100	Inf
m	0	0
τ/s	0	0
c	0	0

计算参数:供电电流I=40 A,电偶极子长度dL=1 800 m,收发距R=9 000 m,赤道偶极测量装置。计算了水平电场E_x全区视电阻率、垂直磁场H_z全区视电阻率,同时与MT卡尼亚视电阻率和CSAMT卡尼亚视电阻率进行比较。从图3的频谱曲线可以看出,大于1 000 Hz高频段全区视电阻率反映出第一层电阻率300 Ω·m,低于1 Hz低频段E_x全区视电阻率和垂直磁场H_z全区视电阻率,与MT卡尼亚视电阻率趋势一致,趋于第二层视电阻率1 000 Ω·m。E_x全区视电阻率在低于0.4 Hz频段变得平直,不再趋于第二层电阻率真值,是由于电场已进入近区,不再反映地电信息;而H_z全区视电阻率在低于0.1 Hz频段还有逼近第二层电阻率真值的趋势。但同样受近区影响,垂直磁场全区视电阻率只响应到0.085 Hz左右,低于这个频率不再反映地电信息;同时也说明垂直磁场进行近区的频率低于水平电场进行近区的频率。而CSAMT卡尼亚视电阻率在低于10 Hz低频段就不能正确反映地层特征。

图3

图3 水平电偶源G型地电模型全区视电阻率频谱

Fig.3 whole zone apparent resistivity spectrum of G-type geoelectricity model with horizontal electric dipole Source

2.3 模型三——H型

地电模型参数如表3所示。

表3 H型地电模型参数

Table 3 H type geoelectricity model parameters

层参数	第一层	第二层	第三层
电阻率/(Ω·m)	300	20	700
层厚/m	1200	300	Inf
m	0	0.35	0
τ/s	0	0.12	0
c	0	0.25	0

计算参数:供电电流I=30 A,电偶极子长度dL=1 800 m,收发距R=7 000 m,赤道偶极测量装置。第二层含有激电信息:m=0.35,τ=0.12 s, c=0.25。计算了水平电场E_x全区视电阻率、垂直磁场H_z全区视电阻率,同时与MT卡尼亚视电阻率和CSAMT卡尼亚视电阻率进行比较。从图4水平电偶源H型型地电模型全区视电阻率频谱曲线可以看出,大于1 000 Hz高频段E_x和H_z全区视电阻率反映出第一层电阻率300 Ω·m,低于1 Hz低频段E_x全区视电阻率和垂直磁场H_z全区视电阻率,与MT卡尼亚视电阻率趋势一致,趋于第三层视电阻率700 Ω·m。E_x全区视电阻率在低于0.1 Hz频段变得平直,受近区影响不再趋于第三层电阻率真值,而垂直磁场H_z全区视电阻率在低于0.1 Hz频段,响应至0.05 Hz,还有逼近第三层电阻率真值的趋势,低于0.05 Hz不再反映地电信息,低于水平电场进入近区的频率。而CSAMT卡尼亚视电阻率在低于5 Hz低频段就不能正确反映地层特征。

图4

图4 水平电偶源H型地电模型全区视电阻率频谱

Fig.4 whole zone apparent resistivity spectrum of H-type geoelectricity model with horizontal electric dipole Source

2.4 模型四——K型

地电模型参数如表4所示。计算参数:供电电流I=30 A,电偶极子长度dL=1 800 m,收发距R=5 000 m,赤道偶极测量装置。计算了水平电场E_x全区视电阻率、垂直磁场H_z全区视电阻率,同时与MT卡尼亚视电阻率和CSAMT卡尼亚视电阻率进行比较。从图5水平电偶源K型地电模型全区视电阻率频谱曲线可以看出,大于800 Hz高频段E_x和H_z全区视电阻率反映出第一层电阻率300 Ω·m,E_x全区视电阻率在低于1 Hz频段变得平直,说明电场已进入近区,不再趋于第三层电阻率真值。低于1 Hz低频段H_z全区视电阻率与MT卡尼亚视电阻率趋势一致,还有逼近第三层电阻率真值200Ω·m的趋势,但同样受近区影响,只响应到0.1 Hz,之后不再反映地电信息。小收发距条件下,垂直磁场进入近区的频率低于水平电场进入近区的频率。而CSAMT卡尼亚视电阻率在低于100 Hz低频段就不能正确反映地层特征。原因是近区H_y场值太小;同时也说明,垂直磁场H_z全区视电阻率,在较小收发距条件下,就能有大探测深度的能力。

表4 K型地电模型参数

Table 4 K type geoelectricity model parameters

层参数	第一层	第二层	第三层
电阻率/(Ω·m)	300	1000	200
层厚/m	1100	400	Inf
m	0	0	0
τ/s	0	0	0
c	0	0	0

图5

图5 水平电偶源K型地电模型全区视电阻率频谱

Fig.5 whole zone apparent resistivity spectrum of K-type geoelectricity model with horizontal electric dipole Source

2.5 模型五——KH型

地电模型参数如表5所示。

表5 KH型地电模型参数

Table 5 KH type geoelectricity model parameters

层参数	第一层	第二层	第三层	第四层
电阻率/(Ω·m)	300	1000	10	500
层厚/m	800	400	300	Inf
m	0	0	0.35	0
τ/s	0	0	0.1	0
c	0	0	0.25	0

计算参数:供电电流I=30 A,电偶极子长度dL=1 800 m,收发距R分别为7 000、12 000、25 000 m,赤道偶极测量装置。第三层含有激电信息:m=0.35, τ=0.1 s, c=0.25。计算了水平电场E_x全区视电阻率、垂直磁场H_z全区视电阻率,同时与MT卡尼亚视电阻率和CSAMT卡尼亚视电阻率进行比较。

从图6水平电偶源KH型地电模型全区视电阻率频谱曲线可以看出,大于1 000 Hz高频段E_x和H_z全区视电阻率与MT卡尼亚视电阻率和CSAMT卡尼亚视电阻率一样,反映出第一层电阻率300 Ω·m,E_x全区视电阻率在低于0.1 Hz频段变得平直,说明已进入近区,不再趋于第三层电阻率真值。由图6c和图6a比较可知,随收发距增大,E_x全区视电阻率也增大,收发距在7 000 m时,0.01 Hz全区视电阻率为166.95 Ω·m;收发距在12 000 m时,0.01 Hz全区视电阻率为216.69 Ω·m;收发距在25 000 m时,0.01 Hz全区视电阻率为320.5 Ω·m;但不能反映地层电阻率真值,进入近区的频率变低。;

图6

图6 水平电偶源KH型地电模型不同收发距全区视电阻率频谱

Fig.6 whole zone apparent resistivity spectrum of KH-type geoelectricity model with horizontal electric dipole source with different transmit-receive distance

H_z全区视电阻率在低于4 Hz低频段与MT卡尼亚视电阻率趋势一致,还有逼近第四层电阻率真值500 Ω·m的趋势,收发距在7 000 m时,只响应到0.03 Hz,低于这个频率进入近区,不再反映地电信息;增大收发距,进入近区的频率变低。随收发距增大,收发距在12 000 m时,低于0.015 Hz进入近区;收发距在25 000 m时,低于0.0045 Hz进入近区。在不同收发距条件下,H_z全区视电阻率进入近区的频率均低于E_x全区视电阻率进入近区的频率。说明H_z全区视电阻率在小收发距条件下更有利于穿透中间层揭示出基底层。而同样受近区影响,CSAMT卡尼亚视电阻率在低于1 Hz低频段就不能正确反映地层特征。

3 结论

通过计算水平电偶源层状介质的电磁场响应,并反算出视电阻率。结果表明,正演计算层状介质的电磁场响应及全区视电阻率结果是正确的,计算程序是可靠的。

计算表明,水平电偶极源电磁场在不同收发距条件下进入近区的频率不同。收发距越大,磁场和水平电场进入近区的频率越低。在近区,水平电场全区视电阻率不随频率变化了,不能进行频率测深;垂直磁场也有最低响应频率,低于这个频率,不能反映地电信息。在相同收发距条件下,垂直磁场H_z进入近区的频率低于水平电场E_x进入近区的频率近一个数量级。

E_x和H_z全区视电阻率以及可控源卡尼亚视电阻率与层状模型大地电磁测深(MT)卡尼亚视电阻率在高频段(>1 000 Hz)一致性较好,而在低频段(<100 Hz)分异较大。收发距越大,相应频率越低。这说明水平电偶极源在高频段、大收发距条件下满足远区条件,类似平面波了;而在低频段、小收发距条件下,受地层电阻率、偶极源长度、收发距等因素影响,电磁场不再满足平面波场条件了。可控源(CSAMT)卡尼亚视电阻率定义进入近区最早,对深层响应能力最低。E_x全区视电阻率在低频段变得平直,也是因为近区影响所致;H_z低于最低响应频率的全区视电阻率也不反映地电信息了,说明已进入近区。

水平电偶极源赤道偶极装置E_x和H_z全区视电阻率均能反映出地层电性特征变化,在大收发距条件下,对高、低中间薄层均能有较好的分辨。在较小收发距条件下,H_z全区视电阻率较E_x全区视电阻率更利于接近深部电阻率真值。对于水平电偶源频率电磁测深,对于大埋深目标体勘查,磁场传感器应具有宽的频率响应,理想的频率响应范围是10⁴~10^-3 Hz。

根据计算结果,在接收条件允许的条件下,大收发距H_z全区视电阻率和E_x全区视电阻率对地层分辨能力好于小收发距。多收发距观测有利于提高对地层的辨识。

致谢

感谢审稿人提出的中肯建议。

参考文献

原文顺序

文献年度倒序

文中引用次数倒序

被引期刊影响因子

[1]

Kaufman A

, Alekseev

, Oristaglio

Principles of electromagnetic methods in surface geophysics

[M]. Elsevier, 2014.

[2]

Zhdanov M

Foundations of geophysical electromagnetic theory and methods

[M]. Elsevier, 2017.

[3]

Chave A

, Jones A

The magnetotelluric method: theory and practice

[M]. Cambridge University Press, 2012.

[4]

Viacheslav

Spichak. Electromagnetic sounding of the Earth’s interior: theory, modeling, practice(2nd ed)

[M]. Elsevier, 2015.

[5]

殷长春

. 航空电磁理论与勘查技术[M]. 北京: 科学出版社, 2018.

Yin C

Aeromagnetic theory and exploration technology [M]. Beijing: Science Press, 2018.

[6]

周海根

多场源地空频率域电磁探测方法研究

[D]. 长春:吉林大学, 2017.

Zhou H

Research on ground-airborne frequency-domain electromagnetic methods with multiple sources

[D]. Changchun: Jilin University, 2018.

[7]

Richard

. Geyer.

Handbook of geophysical exploration at sea(2nd ed)

[M]. CRC Press, 2018.

[8]

王若, 殷长春, 王妙月 , 等.

CSAMT法一维层状介质灵敏度分析

[J]. 地球物理学进展, 2014,29(3):1284-1291,doi: 10.6038/pg20140339.

Magsci [本文引用: 1]

<p>为改善CSAMT一维反演的实用性，本文以电磁场数值模拟为基础，推导出偶极源激发下频率域电磁反演中的偏导数公式.并以三层水平层状地电模型为例，针对同线装置和赤道装置两种形式，利用扰动法对本文的半解析算法结果进行验证.对比结果表明无论同线还是赤道装置，本文半解析结果和扰动法计算结果拟合程度非常好，验证了偏导数公式推导及程序代码的正确性.进而，分析了CSAMT观测数据对中间薄层参数的灵敏度.灵敏度特征分析表明：CSAMT观测资料对中间低阻薄层比较敏感，对高阻薄层的分辨率较低，但对层厚度，即便是高阻薄层，仍存在较高的灵敏度.频率对层参数的灵敏度有重要影响.实际应用中，应同时考虑资料和参数加权，以提高中间层参数的求解能力.本文的算法不仅有助于CSAMT资料反演解释，同时为地下薄目标层的有效识别提供理论依据.</p>

Wang

, Yin C

, Wang M

, et al.

CSAMT sensitivity analysis for 1D models

[J]. Progress in Geophys (in Chinese), 2014,29(3):1284-1291, doi: 10.6038/pg20140339.

Magsci [本文引用: 1]

[9]

李毓茂

, 等. 电磁频率测深方法与电偶源电磁频率测深量板[M]. 徐州: 中国矿业大学出版社, 2012. 7.

[本文引用: 5]

Li Y

, et al. Electromagnetic frequency sounding method and dipole electromagnetic frequency sounding board [M]. Xuzhou: China University of Mining and Technology Press, 2012.

[本文引用: 5]

[10]

Ziolkowski

, Hobbs

, Wright

Multitransient electromagnetic demonstration survey in France

[J]. Geophysics, 2007,72(4):197-207.

[11]

底青云, 朱日祥, 薛国强 , 等.

我国深地资源电磁探测新技术研究进展

[J]. 地球物理学报, 2019,62(6):2128-2138,doi: 10.6038/cjg2019M0633.

Di Q

, Zhu R

, Xue G

, et al.

New development of the Electromagnetic (EM) methods for deep exploration

[J]. Chinese J. Geophys. (in Chinese), 2019,62(6):2128-2138, doi: 10.6038/cjg2019M0633.

[12]

张文伟, 底青云, 雷达 , 等.

物探新方法——多通道瞬变电磁法在金属矿勘探中的应用

[J]. 黄金科学技术, 2018,26(1):1-8.

Zhang W

, Di Q

, Lei

, et al.

Multi-channel transient electromagnetic method:a new geophysical method and its application in exploring metallic ore deposits

[J]. Gold Science and Technology, 2018,26(1):1-8.

[13]

Strack K

Exploration with deep transient electromagnetic

[M]. Elsevier, 1992.

[14]

汤井田, 何继善 . 可控源音频大地电磁法及其应用[M]. 长沙: 中南大学出版社, 2005.

Tang J

, He J

Controlled source audio magnetotelluric method and its application [M]. Changsha, Central South University Press, 2015.

[15]

何继善

. 广域电磁法和伪随机信号电法[M]. 北京: 高等教育出版社, 2010.

He J

Wide-area electromagnetic method and pseudo-random signal electrical method [M]. Beijing: Higher Education Press, 2010.

[16]

罗维斌, 李庆春, 汤井田 .

编码电磁测深

[J]. 地球物理学报, 2012,55(1):341-349, doi: 10.6038/j.issn.0001-5733.2012.01.035.

Magsci [本文引用: 2]

利用逆重复<em>M</em>序列伪随机信号良好的自相关特性,提出相关检测地电传输特性的编码电磁勘探法.分析了编码电磁测深原理及参数提取方法.同时记录发送电流信号和多收发距电磁场响应,对源信号和场信号进行相关运算,解卷积分离接收系统响应后可得到大地的频率特性或时间特性,以此实现地电断面的精细探测.在时间域,大地冲激响应和阶跃响应含有丰富的地电信息.通过冲激响应的峰值时间或阶跃响应的晚期渐近值可估计地电阻率分布.基于层状模型的大地冲激响应和阶跃响应正演计算结果表明,编码电磁测深法对大埋深薄层目标体有精细的分辨能力.可应用于能源、矿产资源、水资源、环境地质及工程地质勘察,有良好的应用前景.

Luo W

, Li Q

, Tang J

Coded source electromagnetic sounding method

[J]. Chinese J. Geopgys. (in Chinese), 2012,55(1):341-349, doi: 10.6038/j.issn.0001-5733.2012.01.035.

Magsci [本文引用: 2]

[17]

罗维斌

伪随机海洋可控源多道电磁测深法研究

[D]. 长沙:中南大学, 2007.

LUO W

Study on pseudorandom marine controlled-source electromagnetic sounding with multi-offsets

[D]. Changsha: Central South University, 2007.

[18]

罗维斌, 汤井田, 韩海涛 , 等.

编码源电磁测深法获取大地全区视电阻率谱的方法及装置,ZL 2016 1 0013962.X

[P/OL]. 2018-03-23.

Luo W

, Tang J

, Han H

, et al.

The method and device of obtaining whole zone apparent resistivity spectrum of the earth by coded source electromagnetic sounding, ZL 2016 1 0013962.X

[P/OL]. 2018-03-23.

[19]

佟铁钢

E-Hz广域电磁方法研究

[D]. 长沙:中南大学, 2010.

Tong T

E-Hz wide area electromagnetic method research

[D]. Changsha: Central South University, 2010.

[20]

Chave A

Numerical integration of related Hankel transforms by quadrature and continued fraction expansion

[J]. Geophysics, 1983,48(12):1671-1677.

[21]

Nielsena T

, Baumgartnera

CR1Dmod: A Matlab program to model 1D complex resistivity effects in electrical and electro-magnetic surveys

[J]. Computers & Geosciences, 2006,32:1411-1419.