2018, Vol. 61
2. 吉林大学仪器科学与电气工程学院, 长春 130026;
3. 吉林大学地球探测科学与技术学院, 长春 130026
2. College of Instrumentation and Electrical Engineering, Jilin University, Changchun 130026, China;
3. College of Geo-Exploration Science and Technology, Jilin University, Changchun 130026, China
固定翼时间域电磁(Fixed-wing Time-Domain Electromagnetic，简称FTEM)系统，采用固定翼飞机作为飞行平台(Fountain, 2008)，具有激发磁矩大、探测效率高等优点，被广泛应用于高山、沙漠、湖泊、沼泽、海陆交互带等复杂区域的快速勘查，在油气、金属矿藏、地下水探测与地质成图等方面发挥越来越重要的作用(Witherly and Sattel, 2012；殷长春等，2015).近年来，地下复杂结构的精细化勘探受到广泛关注，断裂带作为一种复杂的地质结构逐渐成为研究热点.通常，学者们认为断裂带伴随油气产生，并且在油气勘探中起到传输与阻塞的双重作用(沈忠山等，2016).其中活跃的断裂带可以进行油气传输，而静止的断裂带将对油气起到封闭作用(Hooper, 1991；樊计昌和刘明军，2007).
随机介质模型是多种成分随机混合分布的地质结构，自然界的多数物质都符合这一分布规律.典型地质结构如：溶洞、缝洞型油气藏、金属矿、断裂带等，它们的物性参数呈随机分布，因此适用于随机介质模型理论(奚先，2005；Jiang et al., 2013；Ge et al., 2015; Liu et al., 2016).目前，已经证明可以采用统计学方法描述随机介质的电导率分布(Bonnet et al., 2001; Klimeš, 2002).在探地雷达、地震和可控源电磁勘探等领域，学者们已经开展了随机介质的数值模拟研究.Ge等(2012, 2015)基于频率域差分方法，采用Gauss自相关函数、Von Kármán函数建立了任意方向随机裂隙并且模拟了二维、三维接地长导线的可控源电磁响应.Zeng等(2015)采用Gauss自相关函数建立了随机介质模型并且研究了探地雷达的波场特点.郭士礼等(2015)基于随机介质模型理论研究了探地雷达波的随机扰动特征与多相离散随机介质模型参数之间的关系.Sato(2016)在二维大地随机介质中采用差分方法讨论了包络的合成问题.
目前，时间域三维电磁数值模拟主要针对均匀异常体地质结构开展研究.采用有限差分(Finite-Difference Time-Domain，简称FDTD)方法进行数值模拟近年来发展较快.许洋铖等(2012)采用三维有限差分方法实现了全波形时域航空电磁响应的数值计算.李建慧等(2013)采用有限差分方法和有限体积法实现了回线源瞬变电磁法的三维数值模拟.李展辉和黄清华(2014)在时间域数值模拟中采用了复频移PML吸收边界条件.殷长春等(2015)采用积分方程法计算了频率域航空电磁响应，并且经过Hankel变换将其从频率域变换到了时间域.Commer等(2015)采用了一种新的迭代步长算法，提高了有限差分方法的时间域电磁响应计算速度.Chang等(2016)采用有限差分方法模拟了煤炭充水采空区的全波形时间域电磁响应.杨海燕等(2016)采用时间域有限差分方法模拟了覆盖层影响下的地-井瞬变电磁响应.Hu等(2017)结合虚拟波动场电磁传播特点，采用CFS-PML吸收边界条件提高了高阻模型的计算精度.Ji等(2017)采用有限差分方法求解了虚拟波动场扩散方程，实现了起伏地形条件下的时间域电磁响应数值模拟.Commer等(2017)采用有限差分方法模拟了时间域电磁法中的感应-极化效应.余翔等(2017)直接求解了有源Maxwell差分方程，可计算地形影响下的时间域电磁响应.学者们针对三维模型的数值模拟虽然已开展大量研究，但是基于三维断裂带模型的时间域数值模拟与特征分析等研究尚处于起步阶段.
本文采用Von Kármán函数建立随机介质模型，通过一维傅里叶变换得到条状断裂带模型、三维傅里叶变换生成背景围岩电导率模型；并讨论了Hurst指数与自相关长度对电导率的产生范围以及条状断裂带数量及宽度的影响，通过建立三组断裂带模型，分析了均匀断裂带(断裂带的电导率值为一个常数)与随机断裂带(断裂带的电导率随空间位置变化)的响应差异，并且探讨了基于断裂带模型的航空时间域电磁系统的响应特征.
1 三维随机介质建模为了计算随机断裂带地质结构的FTEM响应，需要采用统计学的方法离散地下介质的电导率，一般采用自相关函数描述断裂地质结构(Klimeš, 2002).在建立电导率模型时，通常做法是将地壳模拟为均匀层，而其中的三维异常体的电导率值同样设置为定值.然而在随机介质中，电导率由两部分组成：大尺度变量和小尺度变量.大尺度变量σc为传统意义上的地质模型均匀层电导率，小尺度变量σd(s)为在均匀层电导率上的一个随机扰动，可以采用统计学的方法进行构造，因此最终的介质电导率可以表示为
|
(1) |
其中x=(x, y, z).
1.1 三维随机介质建模方法结合Klimeš(2002)和Ge等(2015)的自相关函数算法，建立三维随机介质电导率模型，具体步骤如下：
(1) 首先采用Yee网格将大地进行空间剖分，确定出与网格步长相关的自相关函数，并且将其变换到波数域
|
(2) |
其中C[r(x)]是自相关函数，可以表示为
|
(3) |
Γ(v)是伽马函数，
(2) 计算白噪声的傅里叶变换N(k).
(3) 将C(k)开平方并和N(k)相乘
|
(4) |
(4) 将P(k)从波数域变换到空间域可以得到随机扰动电导率σd(x).
(5) 最后，将扰动电导率σd(x)与均匀电导率σc求和，得到随机介质电导率σ(x).
1.2 自相关长度和Hurst指数对随机介质模型的影响本节从条状断裂带电导率模型和背景围岩电导率模型两个方面分析了自相关长度和Hurst指数对随机介质建模的影响.结合油气藏断裂带的特征，当空间域到波数域采用一维傅里叶变换时，将生成三维条状断裂带异常体，当采用三维傅里叶变换时，将产生在三个方向随机变化的背景围岩电导率模型.
1.2.1 条状断裂带电导率模型在1.1节中，第1步、第2步和第4步包含傅里叶变换或逆变换，如果对其中的三维函数仅做一维变换，可得到
|
(5) |
|
(6) |
将生成如图 1、图 2所示的条状断裂带模型.图 1为不同自相关长度时的条状随机断裂带电导率分布图，图 1a、1b、1c、1d的随机断裂带模型中自相关长度分别为1、2、3、9.图 2为不同Hurst指数时的条状随机断裂带电导率分布图，图 2a、2b、2c、2d的随机模型中Hurst指数分别为0.01、0.05、0.3、1.
|
图 1 不同自相关长度时随机条状断裂带电导率分布 (a) a=b=c=1; (b) a=b=c=2; (c) a=b=c=3; (d) a=b=c=9. Fig. 1 Conductivity distribution of random streak fractures for different autocorrelation lengths |
|
图 2 不同Hurst指数时随机条状断裂带电导率分布 (a) v=0.01; (b) v=0.05; (c) v=0.3; (d) v=1. Fig. 2 Conductivity distributions of random streak fractures for different Hurst numbers |
如图 1所示，当Hurst指数不变，自相关长度由小变大时，一维傅里叶变换算法可以产生不同数量的平行条状断裂带.当自相关长度增加时，断裂带的数量减少，但宽度增加，自相关长度等于9时，仅有一条断裂带存在.虽然当自相关长度增加时，生成模型的电导率发生变化，但是变化不大，仅从0.01~0.07 S·m-1变化到0.01~0.11 S·m-1，最大值增大了不到一倍.
如图 2所示，当自相关长度不变，Hurst指数由小变大时，条状断裂带的数量保持不变，但是对随机模型电导率的变化范围影响较大，从0.01~0.02 S·m-1变化到0.01~0.54 S·m-1，最大值增大了26倍.
1.2.2 背景围岩电导率模型在1.1节中的傅里叶变换与逆变换采用三维变换时，可得到
|
(7) |
|
(8) |
可生成随机分布的背景围岩电导率模型.图 3为不同自相关长度时的背景围岩分布图，图 3a、3b、3c、3d的随机围岩模型中自相关长度分别为3、10、20、30.图 4为不同Hurst指数时的背景围岩分布图，图 4a、4b、4c、4d的随机模型中Hurst指数分别为0.01、0.1、1、10.
|
图 3 不同自相关长度时背景围岩分布 (a) a=b=c=3; (b) a=b=c=10; (c) a=b=c=20; (d) a=b=c=30. Fig. 3 Conductivity distribution of surrounding rock for different autocorrelation lengths |
|
图 4 不同Hurst指数时背景围岩分布 (a) v=0.01; (b) v=0.1; v=0.1; (c) v=1; (d) v=10. Fig. 4 Conductivity distribution of surrounding rock for different Hurst numbers |
从图 3、图 4中可以看出，当采用三维傅里叶变换算法时，只能产生三维方向的随机变化背景围岩，无法产生条状断裂带模型.如图 3所示，当Hurst指数不变，自相关长度增加时，电导率变化范围同样很小，仅从0.01~0.11 S·m-1变化到0.01~0.17 S·m-1.如图 4所示，对于三维背景围岩模型，当自相关长度保持不变，Hurst指数从0.01变化到10时，电导率从0.01~0.06 S·m-1变化到0.01~1.02 S·m-1，最大值增大了16倍.
分析上述结果，只有当空间域到波数域的变换采用一维傅里叶变换时，才能生成条状断裂带模型，并且自相关长度是确定条状断裂带数量与宽度的重要参数；与自相关长度相比，Hurst指数是决定电导率变化的重要参数.
2 基于FDTD的麦克斯韦方程组离散 2.1 控制方程和模型离散Wang和Hohmann(1993)给出了修正后的麦克斯韦方程组时域形式
|
(9) |
|
(10) |
其中e(x, y, z, t)为电场, h(x, y, z, t)为磁场，μ为磁导率, 
电场在Yee元胞的棱边上采样，磁场在Yee元胞的面中心处采样(Yee, 1966).根据图 5，电场周围环绕了四个磁场，而磁场周围环绕着四个电场.
|
图 5 Yee元胞中的电场和磁场 Fig. 5 Electric field and magnetic field at Yee′s cell |
将麦克斯韦方程组写成标量形式并且用中心差分进行离散，推导后可得到三个电场分量ex、ey、ez与三个磁场分量hx、hy、hz的显式迭代方程，以x方向电场为例，有公式
|
(11) |
其中exn+1(i+1/2, j, k)表示在tn+1时刻，对应于位置(i+1/2, j, k)时x方向的电场.Δy与Δz为空间步长.其它场分量方程可以采用同样的方法得到.
2.2 初始条件FTEM系统一般可近似为偶极-偶极装置，装置示意图如图 6所示.发射线圈与接收线圈距离地面的高度分别为120 m和60 m，发射线圈与接收线圈的水平距离为130 m.大地模型被剖分为175×175×78个网格，计算区域包括吸收边界条件的空间范围为：6870 m×6870 m×3620 m.
|
图 6 FTEM系统装置示意图 Fig. 6 FTEM system sketch map |
数值模拟中利用初始条件替代实际探测中的激励源，在t=t0时刻x与y方向初始电场可表示为(Nabighian，1992)
|
(12) |
|
(13) |
|
(14) |
其中m为磁偶极矩, 
采用F.N.Kong241点的滤波系数(Kong, 2007)计算式(12)和(13)中的Bessel函数积分.频率域到时间域的变换则采用Gaver-Stehfest算法(罗延钟和昌彦君, 2000).在计算Ex, Ey, Ez后，磁场分量可通过麦克斯韦方程组计算得出.
为了保证初始条件的计算精度，初始时刻一定要取的非常小.根据数值经验，初始时刻应为
|
(15) |
其中Δ1和σ1分别为顶层网格步长与电导率.
2.3 稳定性条件在时间域有限差分方法中，为保证计算的稳定，需要遵循CFL(Courant-Friedrichs-Lewy，简称CFL)稳定性条件.在实际的差分离散中，最大时间步长可以表示为(Wang and Hohmann, 1993)
|
(16) |
其中系数α的范围从0.1到0.2.
2.4 吸收边界条件本文采用C-PML(Convolutional Perfect Matched Layer，简称C-PML)作为吸收边界条件.同样以x方向电场为例，吸收边界层中的Maxwell时域方程组标量形式可表示为(Elsherbeni and Demir, 2015)
|
(17) |
其中
仔细研究式(17)可发现，它与传统的控制方程相比，差别在于：磁场的空间偏导多了两个系数以及两项卷积项，在迭代程序中，仅需对C-PML层中的磁场项系数进行修改，最后将两个卷积项计算结果加在ex的表达式中.其它各方向的电场和磁场计算方法与其相似.
2.5 精度验证文中采用均匀半空间模型验证FTEM算法的正确性.均匀半空间模型的垂直磁场频率域半解析解可表示为(Nabighian，1992)
|
(18) |
对磁场而言，基于时间域电磁法观测的是脉冲响应，即负阶跃响应对时间的微分，可表示为
采用以下参数验证算法的精度，大地电导率为0.01 S·m-1，磁偶极矩为7×106 A·m2，程序的运行环境为：Intel(R) Core(TM) i7-5820K CPU@3.30Ghz，NVIDIA GeForce GTX970，基于GPU的并行运算，程序的运行时间为3327 s.图 7a为FDTD方法与半解析解计算结果对比图.图 7b为两者之间的误差曲线.从图中可以看出，FDTD方法求解结果与半解析解吻合很好，最大相对误差为5.754%，晚期误差大主要是由于吸收边界条件引起的.
|
图 7 (a) FDTD(黑色虚线)与半解析解(黑色实线)计算结果对比图，(b)误差曲线 Fig. 7 (a) Comparison of the time derivative of magnetic field between FDTD results (black dotted line) and semi-analytical solution (black solid line), (b) Relative error |
为了分析随机断裂带的时间域电磁响应，我们建立了三种典型断裂带模型，第一种为三条平行走向的断裂带模型；第二种保持断裂带的走向不变，在xoz平面改变它的倾斜角度与倾向；第三种在xoy平面中建立不同倾角的随机断裂带模型.下面详细分析不同三维断裂带的电磁响应特征.
3.1 均匀和随机断裂带的电磁响应特征分析为了分析随机断裂带与均匀断裂带之间的响应差异，在xoz平面建立了产状NE90°，∠0°的随机断裂带和均匀断裂带模型.大地模型被分为两层，第一层电导率是0.01 S·m-1，厚度为150 m，第二层随机电导率为0.01~0.012 S·m-1.当Hurst指数v=0.35，自相关长度a=b=c=2时，共产生8条断裂带，选择其中一条作为基准断裂带，该条断裂带位于第二层模型中间位置，距离地表深度为150 m，长度为4950 m，宽度为170 m，高度为230 m，电导率范围为0.01~0.12 S·m-1.在同一位置，同时建立了均匀断裂带模型，围岩电导率为0.01 S·m-1，为体现出与围岩电导率之间的差异，断裂带电导率取值0.1 S·m-1.
为了更好地分析随机断裂带与均匀断裂带的电磁响应特征，将背景围岩响应剔除，提取了随机和均匀断裂带的纯异常响应.图 8a、8b、8c为不同时刻的沿x方向产状NE90°，∠0°随机和均匀断裂带的电磁响应特征曲线；图 8d、8e、8f为不同时刻沿y方向产状NE90°，∠0°随机和均匀断裂带的电磁响应特征曲线.图 8中黑色实线为均匀断裂带响应特征曲线，黑色虚线为随机断裂带响应特征曲线.从图中可见，均匀断裂带和随机断裂带的特征曲线存在差异，该差异在沿x方向的响应特征曲线上表现的更加明显，最大相对差异率为33.7%，平均相对差异率为4.76%.
|
图 8 产状NE90°，∠0°的均匀断裂带(黑色实线)和随机断裂带(黑色虚线)的电磁响应特征曲线 (a) t=0.5 ms，沿x方向响应曲线; (b) t=3 ms，沿x方向响应曲线; (c) t=8 ms，沿x方向响应曲线; (d) t=0.5 ms，沿y方向响应曲线; (e) t=3 ms，沿y方向响应曲线; (f) t=8 ms，沿y方向响应曲线. Fig. 8 Comparisons of EM responses between the homogeneous fractures(black solid lines) and random fractures(black dotted lines) for the model NE90°, ∠0° (a) EM response along x-axis at t=0.5 ms; (b) EM response along x-axis at t=3 ms; (c) EM response along x-axis at t=8 ms; (d) EM response along y-axis at t=0.5 ms; (e) EM response along y-axis at t=3 ms; (f) EM response along y-axis at t=8 ms. |
图 9a、9c分别为产状NE90°，∠0°随机和均匀断裂带沿x方向的电磁响应二维切片图，切片方向与断裂带方向垂直.图 9b、9d分别为产状NE90°，∠0°随机和均匀断裂带沿y方向的电磁响应二维切片图，切片方向与断裂带方向平行.从图 9中可见，中心区域的响应能够清晰的反应出电磁场扩散特点.对应于均匀断裂带产生的电磁响应，从图 9a中可以清晰的观察到电磁场扩散边界与断裂带的形状相似，但是图 9c中随机断裂带模型所产生的响应扩散边界已经变得模糊不清，无法识别出断裂带的形状.从图 9a、图 9c和图 9b、图 9d中均可以看出，均匀断裂带与随机断裂带之间的响应同样存在差异，这意味着在实际工程探测时，不能将随机断裂带近似为均匀断裂带处理，否则会对数据的解释带来误差.
|
图 9 产状NE90°，∠0°的电磁响应在t=3 ms时的切片图 (a)在y=3000 m, 3500 m, 4000 m时的均匀断裂带响应切片图；(b)在x=3000 m, 3500 m, 4000 m时的均匀断裂带响应切片图；(c)在y=3000 m, 3500 m, 4000 m时的随机断裂带响应切片图；(d)在x=3000 m, 3500 m, 4000 m时的随机断裂带响应切片图 Fig. 9 EM response in slice figure for NE90°, ∠0° at t=3 ms (a) EM response slice figure of homogeneous fracture at y=3000 m, 3500 m, 4000 m; (b) EM response slice figure of homogeneous fracture at x=3000 m, 3500 m, 4000 m; (c) EM response slice figure of random fracture at y=3000 m, 3500 m, 4000 m; (d) EM response slice figure of random fracture at x=3000 m, 3500 m, 4000 m. |
为了分析平行走向随机断裂带的时间域电磁响应特征，选取与断裂带走向垂直的x测线和与断裂带走向平行的y测线电磁响应进行分析.在图 10中建立了三组平行走向的断裂带模型.第一组断裂带位于x=2130 m至x=2300 m之间(简称x=2130 m断裂带)，对应于计算区域的第48~64个空间网格；第二组断裂带位于x=2540 m至x=2710 m之间(简称x=2540 m断裂带)，对应计算区域的第80~96个空间网格；第三组断裂带位于x=2950 m至x=3120 m之间(简称x=2950 m断裂带)，对应计算区域的第112~128个空间网格，随机断裂带的物性参数与3.1节中相同.图 11a、11b、11c分别为不同时刻沿x测线方向的平行走向断裂带电磁响应特征.图 12a、12b、12c分别为不同时刻沿y测线方向的平行走向的断裂带电磁响应特征.
|
图 10 平行条状断裂带的三维模型 (a)平行断裂带主视图; (b) x=2130 m的断裂带俯视图; (c) x=2540 m的断裂带俯视图; (d) x=2950 m的断裂带俯视图. Fig. 10 3D models for parallel streak fractures (a) Section view for parallel fracture; (b) Plane view of random fracture at x=2130 m; (c) Plane view of random fracture at x=2540 m; (d) Plane view of random fracture at x=2950 m. |
|
图 11 沿x测线方向的平行随机断裂带电磁响应 (a) t=0.5 ms时，x=2130 m、x=2540 m、x=2950 m断裂带电磁响应；(b) t=3 ms时，x=2130 m、x=2540 m、x=2950 m断裂带电磁响应；(c) t=8 ms时，x=2130 m、x=2540 m、x=2950 m断裂带电磁响应. Fig. 11 EM responses for parallel random fractures along the x-axis (a) EM responses for parallel random fractures at x=2130, x=2540 and x=2950 when t=0.5 ms; (b) EM responses for parallel random fractures at x=2130 m, x=2540 m and x=2950 m when t=3 ms; (c) EM responses for parallel random fractures at x=2130 m, x=2540 m and x=2950 m when t=8 ms. |
|
图 12 沿y测线方向的平行随机断裂带电磁响应 (a) t=0.5 ms时，x=2130 m、x=2540 m、x=2950 m断裂带电磁响应; (b) t=3 ms时，x=2130 m、x=2540 m、x=2950 m断裂带电磁响应; (c) t=8 ms时，x=2130 m、x=2540 m、x=2950 m断裂带电磁响应. Fig. 12 EM responses for parallel random fractures along the y-axis (a) EM responses for parallel random fractures at x=2130 m, x=2540 m and x=2950 when t=0.5 ms; (b) EM responses for parallel random fractures at x=2130 m, x=2540 m and x=2950 m when t=3 ms; (c) EM responses for parallel random fractures at x=2130 m, x=2540 m and x=2950 m when t=8 ms. |
从图 11a中可以看出，当剖面曲线与断裂带方向垂直，当t=0.5 ms时，位于x=2130 m的断裂带连续出现两正峰一负谷的异常形态，两正峰值相差2.64×10-7 V·m-2，如黑色实线所示；当断裂带位于x=2540 m时，出现单峰正异常为主，响应曲线是近似对称的，如黑色虚线所示；当断裂带位于x=2950 m时，与x=2130 m的模型响应曲线具有对偶关系，如黑色点划线所示.从图 11b、11c中可以看出随着观测时间的加长，响应幅值逐渐减少，并且x=2130 m和x=2950 m的断裂带响应由三个峰值变为仅具有正、负两个峰值.
从图 12a中可以看出，当剖面曲线与断裂带方向平行时，取位于计算区域中心x=2620 m的剖面曲线，此时x=2540 m断裂带在t=0.5 ms时响应最大，左右两侧断裂带响应相对较小.随着时间的增大，平行断裂带的响应幅值减小.比较x=2130 m与x=2950 m断裂带的响应，两者形态相似，幅值相差很少，不易分辨.图 13为平行断裂带的电磁响应切片图，从中可以清晰的看出不同位置平行断裂带的扩散场分布情况.因此，应参考沿x测线方向和y测线方向的电磁响应特征曲线，并结合二维切片图(断裂带两侧的电磁响应发生明显变化)，进行断裂带走向的判断.
|
图 13 t=3 ms，z=305 m时电磁响应切片图 (a) x=2130 m断裂带电磁响应; (b) x=2540 m断裂带电磁响应; (c) x=2950 m断裂带电磁响应. Fig. 13 Slice figures of EM responses at z=305 m, when t=3 ms (a) EM response of random fracture at x=2130 m; (b) EM response of random fracture at x=2540 m; (c) EM response of random fracture at x=2950 m. |
为分析倾角与倾向变化的断裂带，共建立了三个模型.产生断裂带的Hurst指数v=0.5，自相关长度a=b=c=2.5，电导率范围为0.02~0.2 S·m-1.图 14为倾角与倾向变化的随机断裂带模型.图 14a为三个模型的俯视图.图 14b、图 14c、图 14d是在xoz平面产状分别为SW225°, ∠45°、NE90°, ∠0°、SE135°, ∠45°的主视图.图 15a、15b、15c分别为不同产状模型在不同时刻沿x测线方向的断裂带电磁响应曲线.图 16a、16b、16c分别为不同产状模型在不同时刻沿y测线方向的断裂带电磁响应曲线.
|
图 14 倾角与倾向变化的随机断裂带模型 (a)产状为SW225°, ∠45°、NE90°, ∠0°、SE135°, ∠45°的断裂带俯视图; (b)产状为SW225°, ∠45°的断裂带主视图; (c)产状为NE90°且∠0°的断裂带主视图; (d)产状为SE 135°且∠45°的断裂带主视图. Fig. 14 Random fracture changed with dip angle and inclination (a) Plane view of random fracture with SW225°, ∠45°、NE90°, ∠0°、SE 135°, ∠45°; (b) Section view of random fracture with SW 225°, ∠45°; (c) Section view of random fracture with NE90°, ∠0°; (d) Section view of random fracture with SE 135°, ∠45°. |
|
图 15 沿x方向测线的产状为SW225°，∠45°、NE90°且∠0°、SE135°且∠45°的断裂带电磁响应曲线 Fig. 15 Comparisons of EM responses between the models of random fracture SW225°, ∠45°、NE90°, ∠0° and SE135°, ∠45° along the x-axis (a) t=0.5 ms; (b) t=3 ms; (c) t=8 ms. |
|
图 16 沿y方向测线的产状为SW225°，∠45°、NE90°, ∠0°、SE135°, ∠45°的断裂带电磁响应曲线 Fig. 16 Comparisons of EM responses between models of random fracture SW225°, ∠45°, NE90°, ∠0° and SE135°, ∠45° along the y-axis (a) t=0.5 ms; (b) t=3 ms; (c) t=8 ms. |
从图 15a中可以看出，当t=0.5 ms时，产状为SW225°, ∠45°模型的电磁响应连续出现两正峰两负谷的异常形态，且正峰一大一小，两者幅值差为8.29×10-7V·m-2，与平行走向断裂带的大、小两正峰幅值差相比大了近3倍；断裂带产状为NE90°, ∠0°时，单峰正异常为主，响应曲线不完全对称；产状为SE135°, ∠45°与SW225°, ∠45°断裂带模型的响应具有对偶关系.随着观测时间增大，随机断裂带的响应幅值变小，正峰异常向断裂带顶板方向倾斜，负峰异常幅值变大而且变宽.
从图 16可以看出，当断裂带的倾角与倾向变化时，在同一时间取样道沿y测线方向的特征曲线形态相似，仅幅值发生一些变化，因此，仅通过y方向响应较难判断随机断裂带的倾向.
通过分析倾角与倾向变化时的断裂带响应特征，可以得出如下结论：结合不同时刻，从电磁响应曲线上的正峰和负峰幅值、倾斜方向可以判断出断裂带的倾向；与断裂带垂直时测线的电磁响应特征能够更清晰地描述断裂带的倾向，但是无法判断出具体的倾斜角度.
3.4 xoy平面不同倾角的断裂带电磁响应特征分析为了进一步分析随机断裂带的复杂电磁响应特征，在xoy平面中建立了不同倾角的随机断裂带模型.图 17分别为倾角45°、90°和135°的随机断裂带模型图.图 18a、18b、18c分别为沿x轴测线不同时刻时，倾角为45°、90°和135°的电磁响应特征曲线.图 19a、19b、19c分别为沿y轴测线不同时刻时，倾角为45°、90°和135°的电磁响应特征曲线.
|
图 17 xoy平面不同倾角时三维断裂带模型 (a)主视图; (b) 45°倾斜断裂带俯视图; (c) 90°倾斜断裂带俯视图; (d) 135°倾斜断裂带俯视图. Fig. 17 3D random fracture models for different dip angles on the xoy plane (a) Section view; (b) Plane view of 45° slant fracture; (c) Plane view of 90° slant fracture; (d) Plane view of 135° slant fracture. |
|
图 18 xoy平面不同倾角的随机断裂带沿x轴测线的电磁响应 (a) t=0.5 ms时，45°、90°、135°倾斜断裂带电磁响应; (b) t=3 ms时，45°、90°、135°倾斜断裂带电磁响应; (c) t=8 ms时，45°、90°、135°倾斜断裂带电磁响应. Fig. 18 EM responses of random fractures with different dip angle on xoy plane along the x-axis (a) EM responses of random fractures with the dip angles of 45°, 90°, 135° when t=0.5 ms; (b) EM responses of random fractures with the dip angle of 45°, 90°, 135° when t=3 ms; (c) EM responses of random fractures with the dip angles of 45°, 90°, 135° when t=8 ms. |
|
图 19 xoy平面不同倾角的随机断裂带沿y轴测线的电磁响应 (a) t=0.5 ms时，45°、90°、135°倾斜断裂带电磁响应; (b) t=3 ms时，45°、90°、135°倾斜断裂带电磁响应; (c) t=8 ms时，45°、90°、135°倾斜断裂带电磁响应. Fig. 19 EM responses of random fractures with different dip angle on xoy plane along the y-axis (a) EM responses of random fractures with the dip angles of 45°, 90°, 135° when t=0.5 ms; (b) EM responses of random fractures with the dip angle of 45°, 90°, 135° when t=3 ms; (c) EM responses of random fractures with the dip angles of 45°, 90°, 135° when t=8 ms. |
从图 18a中可以看出，在t=0.5 ms时，45°、90°和135°倾斜断裂带的电磁响应形态基本一致，倾角为90°的断裂带电磁响应幅值最大.早期电磁响应受断裂带倾角影响较小，当测线方向垂直于断裂带时，可以获得最大电磁响应；从图 18b中可以看出，在t=3 ms，90°断裂带的响应仍然较大，45°断裂带先出现峰值，其次是135°；t=8 ms时的响应差别较小，假设测量系统的最小分辨率为5 nV，此时的电磁响应已经无法区分.从图 19中可以看出，当断裂带倾角为45°和135°时，沿y轴测线的电磁响应基本重合，因此很难分辨；在t=0.5 ms和t=3 ms时，90°断裂带的响应特征形态与45°和135°倾斜断裂带形态相似，但响应幅值要更大.
在判别xoy平面倾斜断裂带时，可以通过分析早期响应幅值大小、观察峰值出现的早晚并需结合响应的二维切片图进行综合判别.由于沿x测线方向所成剖面曲线与断裂带的走向不再垂直，因此图 18、图 19中的曲线不如3.2与3.3中的响应曲线特征明显，所以在实际测量中，采用与断裂带垂直方向所接收到的响应值对断裂带的判别最有效.
4 结论采用Von Kármán自相关函数建立了随机断裂带模型, 基于有限差分方法直接实现了随机断裂带模型的时间域电磁响应数值求解，基于均匀半空间模型验证了算法的有效性与精度，通过分析三维随机断裂带的建模及电磁响应特征得出以下结论：
(1) 在随机断裂带模型建立过程中，采用一维傅里叶算法实现空间域到波数域的变换时，将产生条状断裂带异常体，而采用三维傅里叶变换将产生随机分布的背景围岩.
(2) 在建立断裂带模型时，Hurst指数与自相关长度两个参数起到非常重要的作用.当建立条状断裂带时，自相关长度是确定断裂带数量与宽度的重要参数.Hurst指数对电导率变化范围的影响更大.
(3) 通过数值模拟分析，随机断裂带和均匀断裂带产生的电磁响应存在差异，因此不能简单的将随机断裂带近似为均匀断裂带，否则会对数据解释带来误差.
(4) 通过分析随机断裂带的电磁响应特征，随机断裂带模型可以准确描述地下介质的分布特征，而且垂直于断裂带测线方向的电磁响应特征清晰地描述了断裂带的倾向、走向与位置，为随机断裂带结构探测提供了理论依据和技术指导.
随机断裂带模型的建立以及电磁响应特征分析尚处于起步阶段，断裂带的电磁探测识别方法，以及空间分数阶差分方程的求解将是下一步研究的重点.
致谢 感谢编辑与审稿人对论文提出的宝贵意见.感谢李静老师在论文修改过程中的帮助.本文受国家自然科学基金项目(41604084，41674109)以及吉林省优秀青年人才基金(20170520104JH)项目的资助，深表感谢.
Bonnet E, Boun O, Odling N E, et al. 2001. Scaling of fracture systems in geological media. Reviews of Geophysics, 39(3): 347-383. DOI:10.1029/1999RG000074 |
Chang J H, Yu J C, Liu Z X. 2016. Three-dimensional numerical modeling of full-space transient electromagnetic responses of water in goaf. Applied Geophysics, 13(3): 539-552. DOI:10.1007/s11770-016-0572-y |
Commer M, Hoversten G M, Um E S. 2015. Transient-electromagnetic finite-difference time-domain earth modeling over steel infrastructure. Geophysics, 80(2): E147-E162. DOI:10.1190/geo2014-0324.1 |
Commer M, Petrov P V, Newman G A. 2017. FDTD modelling of induced polarization phenomena in transient electromagnetics. Geophysical Journal International, 209(1): 387-405. |
Elsherbeni A Z, Demir V. 2015. The Finite-Difference Time-Domain Method for Electromagnetics with MATLAB Simulations. 2nd ed. Edison, NJ: SciTech Publishing, an Imprint of the IET, 229-275.
|
Fan J C, Liu M J. 2007. An approach determining internal structure and physical parameters of fractural zone. Oil Geophysical Prospecting (in Chinese), 42(2): 164-169. |
Fountain D. 2008. 60 Years of airborne EM-focus on the last decade.//AEM2008-5th International Conference on Airborne Electromagnetics. Finland: Haikko Manor.
|
Ge J C, Everett M E, Weiss C J. 2012. Fractional diffusion analysis of the electromagnetic field in fractured media part I:2D approach. Geophysics, 77(4): WB213-WB218. DOI:10.1190/geo2012-0072.1 |
Ge J C, Everett M E, Weiss C J. 2015. Fractional diffusion analysis of the electromagnetic field in fractured media——part 2:3D approach. Geophysics, 80(3): E175-E185. DOI:10.1190/geo2014-0333.1 |
Guo S L, Ji M E, Zhu P M, et al. 2015. Study on multiphase discrete random medium model and its GPR wave field characteristics. Chinese J. Geophys. (in Chinese), 58(8): 2779-2791. DOI:10.6038/cjg20150813 |
Hooper E C D. 1991. Fluid migration along growth faults in compacting sediments. Journal of Petroleum Geology, 14(S1): 161-180. |
Hu Y P, Egbert G, Ji Y J, et al. 2017. A novel CFS-PML boundary condition for transient electromagnetic simulation using a fictitious wave domain method. Radio Science, 52(1): 118-131. DOI:10.1002/2016RS006160 |
Ji Y J, Hu Y P, Imamura N. 2017. Three-dimensional transient electromagnetic modeling based on fictitious wave domain methods. Pure and Applied Geophysics, 174(5): 2077-2088. DOI:10.1007/s00024-017-1528-8 |
Jiang Z M, Zeng Z F, Li J, et al. 2013. Simulation and analysis of GPR signal based on stochastic media model with an ellipsoidal autocorrelation function. Journal of Applied Geophysics, 99: 91-97. DOI:10.1016/j.jappgeo.2013.08.005 |
Klimeš L. 2002. Correlation functions of random media. Pure and Applied Geophysics, 159(7-8): 1811-1831. DOI:10.1007/s00024-002-8710-2 |
Kong F N. 2007. Hankel transform filters for dipole antenna radiation in a conductive medium. Geophysical Prospecting, 55(1): 83-89. DOI:10.1111/gpr.2007.55.issue-1 |
Li J H, Hu X Y, Zeng S H, et al. 2013. Three-dimensional forward calculation for loop source transient electromagnetic method based on electric field Helmholtz equation. Chinese J. Geophys. (in Chinese), 56(12): 4256-4267. DOI:10.6038/cjg20131228 |
Li Z H, Huang Q H. 2014. Application of the complex frequency shifted perfectly matched layer absorbing boundary conditions in transient electromagnetic method modeling. Chinese J. Geophys. (in Chinese), 57(4): 1292-1299. DOI:10.6038/cjg20140426 |
Liu D Y, Han L G, Zhang P, et al. 2016. Comparison of finite difference and pseudo-spectral methods in forward modelling based on metal ore model of random media. Global Geology, 19(2): 102-108. |
Luo Y Z, Chang Y J. 2000. A rapid algorithm for G-S transform. Chinese J. Geophys. (in Chinese), 43(5): 684-690. |
Nabighian M N. 1992. Electromagnetic Methods in Applied Geophysics (in Chinese). Zhao J X, Wang Y J, trans. Beijing: Geological Publishing House, 191-231.
|
Sato H. 2016. Envelope broadening and scattering attenuation of a scalar wavelet in random media having power-law spectra. Geophysical Journal International, 204(1): 386-398. DOI:10.1093/gji/ggv442 |
Shen Z S, He X, Mei M, et al. 2016. Geological characteristics of the faulted zone in North Xingshugang oilfield of Songliao Basin. Petroleum Geology and Oilfield Development in Daqing (in Chinese), 35(3): 22-25. |
Wang T, Hohmann G W. 1993. A finite-difference, time-domain solution for three-dimensional electromagnetic modeling. Geophysics, 58(6): 797-809. DOI:10.1190/1.1443465 |
Witherly K, Sattel D. 2012. The application of ZTEM to porphyry copper-gold exploration.//82nd Ann. Internat Mtg., Soc. Expi. Geophys.. Expanded Abstracts, 26-29. http://www.publish.csiro.au/EX/ASEG2012ab329
|
Xi X, Yao Y, Gu H M. 2005. Simulations of random cave medium model. Journal of Huazhong University of Science and Technology (Nature Science Edition) (in Chinese), 33(9): 105-108. |
Xu Y C, Lin J, Li S Y, et al. 2012. Calculation of full-waveform airborne electromagnetic response with three-dimension finite-difference solution in time-domain. Chinese J. Geophys. (in Chinese), 55(6): 2105-2114. DOI:10.6038/j.issn.0001-5733.2012.06.032 |
Yang H Y, Yue J H, Xu Z Y, et al. 2016. Transient electromagnetic method modeling in ground-borehole model with overburden influence. Journal of Jilin University (Earth Science Edition) (in Chinese), 46(5): 1527-1537. |
Yee K S. 1966. Numerical solution of initial boundary value problems involving Maxwell's equations in isotropic media. IEEE Transactions on Antennas and Propagation, 14(3): 302-307. DOI:10.1109/TAP.1966.1138693 |
Yin C C, Zhang B, Liu Y H, et al. 2015. Review on airborne EM technology and developments. Chinese J. Geophys. (in Chinese), 58(8): 2637-2653. DOI:10.6038/cjg20150804 |
Yu X, Wang X B, Li X J, et al. 2017. Three-dimensional finite difference forward modeling of the transient electromagnetic method in the time domain. Chinese J. Geophys. (in Chinese), 60(2): 810-819. DOI:10.6038/cjg20170231 |
Zeng Z F, Chen X, Li J, et al. 2015. Recursive impedance inversion of ground-penetrating radar data in stochastic media. Applied Geophysics, 12(4): 615-625. DOI:10.1007/s11770-015-0514-0 |
樊计昌, 刘明军. 2007. 确定断裂带内部结构和物性参数的一种方法. 石油地球物理勘探, 42(2): 164-169. DOI:10.3321/j.issn:1000-7210.2007.02.008 |
郭士礼, 冀孟恩, 朱培民, 等. 2015. 多相离散随机介质模型及其探地雷达波场特征研究. 地球物理学报, 58(8): 2779-2791. DOI:10.6038/cjg20150813 |
李建慧, 胡祥云, 曾思红, 等. 2013. 基于电场Helmholtz方程的回线源瞬变电磁法三维正演. 地球物理学报, 56(12): 4256-4267. DOI:10.6038/cjg20131228 |
李展辉, 黄清华. 2014. 复频率参数完全匹配层吸收边界在瞬变电磁法正演中的应用. 地球物理学报, 57(4): 1292-1299. DOI:10.6038/cjg20140426 |
罗延钟, 昌彦君. 2000. G-S变换的快速算法. 地球物理学报, 43(5): 684-690. DOI:10.3321/j.issn:0001-5733.2000.05.012 |
Nabighian M N.1992.勘查地球物理电磁法.赵经祥, 王艳君, 译.北京: 地质出版社, 191-231.
|
沈忠山, 何欣, 梅梅, 等. 2016. 松辽盆地杏树岗北部油田断裂带地质特征. 大庆石油地质与开发, 35(3): 22-25. DOI:10.3969/J.ISSN.1000-3754.2016.03.004 |
奚先, 姚姚, 顾汉民. 2005. 随机溶洞介质模型的构造. 华中科技大学学报(自然科学版), 33(9): 105-108. DOI:10.3321/j.issn:1671-4512.2005.09.033 |
许洋铖, 林君, 李肃义, 等. 2012. 全波形时间域航空电磁响应三维有限差分数值计算. 地球物理学报, 55(6): 2105-2114. DOI:10.6038/j.isn.0001-5733.2012.06.032 |
杨海燕, 岳建华, 徐正玉, 等. 2016. 覆盖层影响下典型地-井模型瞬变电磁法正演. 吉林大学学报(地球科学版), 46(5): 1527-1537. |
殷长春, 张博, 刘云鹤, 等. 2015. 航空电磁勘查技术发展现状及展望. 地球物理学报, 58(8): 2637-2653. DOI:10.6038/cjg20150804 |
余翔, 王绪本, 李新均, 等. 2017. 时域瞬变电磁法三维有限差分正演技术研究. 地球物理学报, 60(2): 810-819. DOI:10.6038/cjg20170231 |