2. 中国科学院地质与地球物理研究所, 北京 100029
2. Institute of Geology and Geophysics, Chinese Academy of Sciences, Beijing 100029, China
瞬变电磁拟地震延拓成像技术是一个较前沿的研究热门[1~4].瞬变电磁场满足的微分方程实际上是一个扩散方程,不利于传统意义上的成像.国内、外学者一直在寻找一种数学上的解决方法,即通过数学积分变换,将满足扩散方程的瞬变电磁场转换为满足波动方程的虚拟波场[5~8],实现这种转换后,可以借助于地震勘探中发展起来的一些比较成熟的成像方法技术,解决复杂情况下瞬变电磁场的成像问题.
1989年Lee[5]建立了满足时域扩散方程的电磁场与虚拟波场的关系式,通过这个变换式,可以将对应地电模型的波动方程模拟结果变换成时域电磁响应,但该变换的意义并不在于对时域电磁场的正演模拟,而在于其逆变换过程.这种逆变换过程是病态的,若不采用优化算法和正则化处理,将会使干扰扩大而破坏结果.文献[9, 10]提出了从瞬变电磁场到波场的优化算法,在保证瞬变电磁场的计算精度的同时,应用两步最优化和正则化算法成功地解决了波场反变换问题,大量的理论计算表明:这种变换得到的虚拟波场,不仅满足波动方程,而且还类似于地震子波,具有传播、反射、透射等特征.这就为瞬变电磁场的偏移成像创造了条件.
本文在以往瞬变电磁场的波场变换计算的基础上,用克希霍夫积分法对瞬变电磁虚拟波场进行延拓成像处理,实现了地表以为曲面的向下延拓成像.通过对理论模型计算和实际资料处理,证明了该方法可以增强瞬变电磁法识别地下目标体的能力,使瞬变电磁法对地的三维精细解释成为可能.
2 基本原理 2.1 瞬变电磁场的波场变换Lee[5]以电性源为基础,建立了满足时域扩散方程的电场强度与虚拟波场的关系式,依据这一思路,本文从麦克思韦方程出发,建立以大回线源为基础的时域瞬变响应与虚拟波场的对应关系.
在导电介质中,如果忽略位移电流,磁场满足如下扩散方程:
(1) |
引入一函数U(r,τ)使其满足如下波动方程:
(2) |
分别对扩散方程(1)中的Hm(r,t)从t到s和波动方程(2)中的U(r,τ)从τ到p进行Laplace变换,得:
(3) |
(4) |
其中
如果令s=p2,则方程(3)变为
(5) |
对比(4)和(5)式,可证得
成立,去掉空间变量r,并写成标量形式:
由于
将上式两端进行Laplace反变换,得到从波场到时域场的波场正变换式[5]:
(6) |
波场正变换(6)式是典型的第一类Fredholm型算子方程,它的反问题具有“不适定性”[9].目前野外实测数据的时间范围一般为几十微妙到几十毫秒,如果直接用(6)式的反变换求波场,不适定性是相当严重的,为了保证反变换式中系数矩阵的条件数不至于过大,可将该时间范围分为7个时段来计算,即:第1时段:32.5~80μs;第2时段:80~325μs;第3时段:325~800μs;第4时段:800μs~2.4ms;第5时段:2.4~8.7ms;第6时段:8.7~27 ms;第7时段:27~81ms.利用文献[9]中的方法,采用最优化技术和正则化算法首先解决了(6)式的反变换问题,实现了从时域场到波场的变换.
通过波场变换获得的图像,能够显示不同电性层位的位置,但由于反射面的“强弱”不同,使得结果波形能量相差悬殊.对于一些能量过弱的波形甚至无法识别,以致遗漏层位.为改善剖面的成像效果,需采用相关的处理手段如:地震中常用的能量均衡处理手段[11~13],该方法在一定程度上减少了遗漏分层的风险.另外,在瞬变电磁场转化为虚拟波场的过程中,所得虚拟波场波形会随着偏移距和地层电导率的增大而不可避免地展宽,严重影响成像的分辨率,这一现象与地震勘探中的“大地滤波作用”非常相似,从滤波角度看,要提高分辨率抵消大地滤波器的作用,需求取一反滤波因子,令其抵消大地滤波的作用[13],即可大幅提高分辨率.
将瞬变场转换成虚拟波场,就是提取电磁响应中与传播有关的特征,而压制或去除电磁波传播过程中与频散、衰减有关的特性.这种变换得到的虚拟波场,完全可以借助于地震勘探解释中的方法,实现瞬变电磁场的拟地震解释.
2.2 瞬变电磁场拟波动方程三维曲面延拓成像 2.2.1 瞬变电磁虚拟波场克希霍夫积分(曲面延拓)将瞬变电磁扩散场变换为虚拟波场后,就可以将地震波场中的偏移成像方法用于对瞬变电磁场的解释[14~17],也就是用克希霍夫积分法进行瞬变电磁虚拟波场偏移成像处理,就可以实现虚拟电磁波场的延拓计算.
波场在地下传播可以用波动方程描述:
(7) |
上述方程的克希霍夫积分解为[13]
(8) |
其中Q为包围波场分布空间的外边界,Q=Q0 + Q1,Q0为地面,Q1为无限大半球面,v为波场传播速度,F=μ0δ(t-o),如图 1所示.考虑到波场是随着距离衰减的,故Q1上的积分可以不考虑.
对于u(x,y,z,t),当把t改变成-t时,既是w(x,y,z,t)=u(x,y,z,-t)时,w(x,y,z,t)仍可满足同样的波动方程.对于w是时间向前的问题,对于u就是时间“倒退”的问题.可以把反射界面的各点看为同时激发上行波的源点.这样,可以将地面上的接收点接收到的信号视为由反射界面的“二次发射源”产生的,将这些信息值时间“倒退”到“二次发射源”,寻找反射界面的波场函数,以确定反射界面.
设自激自收的上行波为G(x,y,z0,t),它是地下反射界面作为源点发射的波场g(x,y,z,t)在地面上z=z0上的值,由(8)式可得:
(9) |
式中,
将三维边界元技术引入到克希霍夫积分计算中,其优点是能够实现三维曲面延拓计算,这在地震勘探中是很难实现的,由于公式(9)中存在∂G/∂n项,地震勘探中无法进行地表垂向梯度测量,相反瞬变电磁法利用梯度探头可以很容易地实现地表垂向梯度测量.通过对时间域的瞬变电磁场进行逆时偏移成像的深入研究,可将瞬变电磁虚拟波场从地面向地下反向外推进行偏移成像,建立有瞬变电磁法特点的偏移成像方法.
2.2.2 波场延拓的边界单元法(1)克希霍夫积分的离散化
用三角形单元对边界Q0进行剖分,共剖分成ne个单元,n个节点.将(9)式中的边界积分分解为诸单元积分的积分之和,当p点位于Ω中时
(10) |
其中
(2)单元分析
设三角形单元顶点编号为j,k,m;各节点坐标分别为(xj,yj,zj),(xk,yk,zk),(xm,ym,zm),如图 2所示,首先做如下坐标变换:
用形函数ξj,ξk,ξm表示单元中任意点的坐标.
(11) |
其中ξj,ξk,ξm是0→1的函数,且ξj+ξk+ξm=1,
(12) |
可见,ξj,ξk,ξm是x,y,z的线性函数.
由于单元Qe一般取得很小,故可假定波场G在各单元上是线性变化的,即单元上的G可表示为
(13) |
式中,Gj,Gk,Gm是地面Q0节点上的G值.
考虑单元积分:
其中
(14) |
其中,
同理,单元积分
(15) |
其中,
(16) |
单元积分:
(17) |
其中,
(18) |
(3)合成矩阵
将(14)、(15)、(17)式代入(10)式中,i点处的波场为
由n个节点得一方程组:
(19) |
其中,F=[Fij],D=[Dij],C=[Cij],Fij为j节点上相邻单元fij之和;Dij为j节点上相邻单元dij之和,Cij为j节点上相邻单元cij之和.若已知地表上的波场值
在以上的瞬变电磁虚拟波场延拓成像计算中,直接影响延拓成像质量的连续速度建模技术显得十分重要.(9)式中的
由等效导电平面法可获得地电断面总纵向电导:
(20) |
其中H为地层埋深,h为地层厚度.利用相邻层的纵向电导可以导出第i层的电导率值
(21) |
某一时刻地下瞬变电磁虚拟波场延拓点的瞬时速度为
(22) |
利用纵向电导与虚拟波速的对应关系,将纵向电导的数据转换为虚拟速度数据体,从实际测量的数据中得到的纵向电导数据体不可能很大,在空间分布上也会表现得较为稀少,为了解决数据不够丰富的问题,使用近点线性三维空间插值计算的办法,来使该数据量在保证正确的基础上得以扩大,实现虚拟波场连续速度建模.
3 模型计算在地表布置面积工作如图 3所示,共布置11条剖面,每条剖面布置11个测点,取发射回线边长为200m,时窗采样宽度为32.5μs-8.7 ms,对全部121个测点进行正演计算,对所得电磁响应进行波场变换,得到各测点的虚拟波场值,以数据文件的形式输入计算机,计算地电断面上的波场延拓值.
给出一个瞬变电磁中心回线装置下的三层A型地电模型(ρ1=50 Ωm,ρ2=500 Ωm,ρ3=5000 Ωm,h1=100m,h2=200m).在这个模型条件下,对每个测点进行正演计算,得到电磁响应值,之后进行波场变换,为了清楚起见,只对变换后的波峰进行显示(如图 4所示).
从图中可以看出,随着虚拟时间τ的增大,虚拟波场出现了两个正的峰值变化,这两个峰值就标示出了这个A型地电模型的两个地下分界面的位置.
利用每个测点的电磁响应值,采用式(22)计算虚拟速度值.速度计算的结果如图 5所示.可以看到,虚拟速度的分布明显呈现三层分布,与模型的电阻率值分布是对应的,这是因为虚拟波场的速度值取决于地下电导率的分布,因此,虚拟速度表现出与电阻率一致的三层分布规律.
最后,将三维边界元技术应用克希于霍夫积分计算中,利用速度分析结果实现三维延拓成像,得到图 6所示的结果,两个峰值显示明显表示出了A型地电模型的两个地下分界面的位置.这两个位置所对应的深度也恰好在100m与300m的附近,与模型的参数基本吻合.
在煤矿开采过程中,积水、采空区、陷落柱等水体病害一直困扰着生产.在山西某煤矿区开展采空区调查,目的是要搞清煤田地下充水采空区的采空深度、采空高度、分布范围,避免在采掘过程中发生突水灾害.
勘查区表层主要为第四系黄土覆盖,厚度一般不超过20m,下伏地层主要为第三系砂砾岩层、二叠系砂泥岩夹煤系地层、石炭系砂泥岩夹煤系地层及奥陶系灰岩地层,该区分布的主要煤层有两层,分别位于地下150m和260m处.各地层岩矿石电阻率差异较小,难以引起明显的电阻率异常,当地下煤层局部被采出后,在岩体内形成一个有一定规模的空间,使周围的应力平稳状态遭受破坏,产生局部的应力集中,采空区顶板在上覆岩层压力的作用下,发生变形、断裂、位移、冒落,形成的冒落带、断裂带、变形弯曲带,在地下水的充填及地表水沿裂缝向采空区渗漏,其电阻率将明显发生变化,形成一个低阻电性体,也与围岩电性形成较明显的差异,这样就为采空区的勘查提供了前提条件.
根据该工区附近的采空区的埋藏深度和厚度情况,结合瞬变电磁特点,选择大定源回线装置,发射回线边长为200m,测线间距20m,测点间距20m,布置了面积工作.
4.1 波场变换结果图 7是经波场变换,并做均衡处理、反褶积、多点平滑后得到最终结果的三维成像全图,波场值颜色深的部分为虚拟波能量高的区域,也是反射界面所在处.
为了更好地展示采空区的空间分布,将波场延拓成像的结果采用两种显示方式,一种是去掉低能反射(低速或者低阻分布)部分如图 8a所示,可以清楚地显示采空区在上下围岩中的分布位置.另一种是去掉高能反射(高速或者高阻分布)部分如图 8b所示,展示出了采空区的空间分布形态.通过对成像结果进行解释,认为有两层采空区分布,上层采空区位于地下140m处,充水较少,冒落带分部范围小;下层采空区位于地下260 m处,充水较多,冒落带分部范围大,由于充水的影响会使异常产生放大作用,采空区的分布范围也会有一定的放大.经验证解释结果与实际地质情况吻合较好,这也证明了本文算法的有效性.
满足扩散方程的瞬变电磁场经积分逆变换后,转换为满足波动方程的虚拟波场.这一变换的实质是提取瞬变电磁响应中与传播有关的特征,压制或去除电磁波传播过程中与频散和衰减有关的特征.将瞬变电磁场的求解问题转化为波动方程的求解问题,这样就能将地震偏移成像技术用于瞬变电磁场的反演解释中.
用本文提出的方法对各种模型进行了数值计算,同时对野外实际资料进行了处理,处理结果表明,瞬变电磁虚拟波场三维曲面延拓成像法理论及方法是有效的,实际例子的成像结果与地质资料吻合,表明该法在实际中的应用是可行的.
[1] | 何继善. 电法勘探的发展和展望. 地球物理学报 , 1997, 40(Suppl): 308–316. He J S. Development and prospect of electrical prospecting method. Chinese J. Geophys. (in Chinese) , 1997, 40(Suppl): 308-316. |
[2] | 于鹏, 王家林, 吴健生. 大地电磁场成像方法综述与新进展. 地球物理学进展 , 2003, 18(1): 53–58. Yu P, Wang J Y, Wu J S. Summary and advance looking of magnetotelluric field imaging method. Progress in Geophysics (in Chinese) , 2003, 18(1): 53-58. |
[3] | 李貅, 薛国强, 郭文波. 瞬变电磁法拟地震成像研究进展. 地球物理学进展 , 2007, 22(3): 811–816. Li X, Xue G Q, Guo W B. Research progress in TEM pseudo-seismic imaging. Progress in Geophysics (in Chinese) , 2007, 22(3): 811-816. |
[4] | 薛国强, 李貅, 底青云. 瞬变电磁法正反演研究进展. 地球物理学进展 , 2008, 23(4): 1165–1172. Xue G Q, Li X, Di Q Y. Research progress in TEM forward modeling and inversion calculation. Progress in Geophysics (in Chinese) , 2008, 23(4): 1165-1172. |
[5] | Lee K H. Liu G, Morrison H F. A new approach to modeling the electromagnetic response of conductive media. Geophysics , 1989, 54(6): 1180-1192. |
[6] | Lee S, Memechan G A. Phase-field imaging:The electromagnetic equivalent of seismic migration. Geophysics , 1987, 52(5): 678-693. DOI:10.1190/1.1442335 |
[7] | Zhdannov M S, Traynin P N. Resistivity imaging by time domain electromagnetic migration. Exploration Geophysics , 1995, 26: 186-194. DOI:10.1071/EG995186 |
[8] | Zhdanov M S, Portniaguine O. Time-domain electro-magnetic migration in the solution of inverse problems. Geophysics , 1997, 62(2): 293-309. |
[9] | 李貅, 薛国强, 宋建平, 等. 从瞬变电磁场到波场的优化算法. 地球物理学报 , 2005, 48(5): 1185–1190. Li X, Xue G Q, Song J P, et al. An optimize method for transient electromagnetic field-wave field conversion. Chinese J. Geophys. (in Chinese) , 2005, 48(5): 1185-1190. |
[10] | Xue G Q, Yan Y J, Li X. Pseudo-Seismic Wavelet Transformation of transient electromagnetic Response in geophysical exploration. Journal of Geophysical Research Letters , 2007, 34: L16405. DOI:10.1029/2007GL031116, 2007 |
[11] | 李蓉, 胡天跃. 时移地震资料处理中的互均衡技术. 石油地球物理勘探 , 2004, 39(4): 424–427. Li R, Hu T Y. Cross-equilibration technique in time-lapse seismic data processing. Oil Geophysical Prospecting (in Chinese) , 2004, 39(4): 424-427. |
[12] | 金龙, 陈小宏, 刘其成. 基于奇异值分解的时移地震互均衡方法. 地球物理学进展 , 2005, 20(2): 294–297. Jin L, Chen X H, Liu Q C. New method for time-lapse seismic matching based on SVD. Progress in Geophysics (in Chinese) , 2005, 20(2): 294-297. |
[13] | 何樵登. 地震勘探原理和方法. 北京: 地质出版社, 1985 : 13 -176. He Q D. Principles and Methods of Seismic Exploration (in Chinese). Beijing: Geological Publishing House, 1985 : 13 -176. |
[14] | 吕国印. 瞬变电磁法二维逆时电磁偏移. 物探与化探 , 1998, 22(2): 139–142. Lü G Y. Two-dimensional counter-time shift in transient electromagnetic method. Geophysical & Geochemical Exploration (in Chinese) , 1998, 22(2): 139-142. |
[15] | Zhdanov M S, Booker J R. Underground imaging by electromagnetic migration.In:63rd Ann.Internat.Mtg. Expl.Geophys., Expanded Abstracts, 1993.355~357 |
[16] | Zhdanov M S, Li W D. 2D finite-difference time domain electromagnetic migration.In:67th SEG EM2.1.1997.370~373 |
[17] | Zhdanov M S, Portniaguine O. Time-domain electro-magnetic migration in the solution of inverse problems. Geophysics , 1997, 62(2): 293-309. |