2016, Vol. 59
2. 中国科学院大学, 北京 100049;
3. 长安大学地质工程与测绘学院, 西安 710054;
4. 中国地质科学院地球物理地球化学勘查研究所, 河北廊坊 065000
2. University of Chinese Academy of Sciences, Beijing 100049, China;
3. School of Geology and Engineering, Chang'an University, Xi'an 710054, China;
4. Institute of Geophysical and Geochemical Exploration, Chinese Academy of Geological Sciences, Hebei Langfang 065000, China
多道瞬变电磁法(Multi-channel transient electromagnetic，MTEM)是近年来新发展起来的电磁探测新技术(Wright et al., 2001, 2002; Ziolkowski et al., 2007).与地震法勘探类似，MTEM技术采用一发多收、阵列采集的方式进行数据接收，可以得到与地震数据类似的共偏移距道集、共中心点道集数据体(Li et al., 2015).因此，采用瞬变电磁拟地震解释方法对MTEM数据进行处理解释是提高MTEM数据解释精度的研究方向之一(李貅和薛国强, 2013, Xue et al., 2012).虚拟波场提取是瞬变电磁拟地震处理的前提，因此研究MTEM虚拟波场提取方法，具有重要意义.
国内学者已经对常规瞬变电磁法虚拟波场的提取方法进行了研究.陈本池等(1999)采用奇异值分解法进行虚拟波场提取，并对算法中的影响因素进行了系统的研究和计算.李貅等将正则化算法引入到虚拟波场提取过程中，采用时间窗口分段计算的方式，改善了系数矩阵的病态程度，使得虚拟波场解更加稳定和可靠(李貅等, 2005, 2010；薛国强等，2006；Xue et al., 2011).戚志鹏等在李貅研究基础上，将正则化共轭梯度法与预条件共轭梯度法结合，采用预条件正则化共轭梯度法进行虚拟波场提取，提高了提取精度(薛国强等，2011；戚志鹏等，2013).
本文提出一种新的MTEM虚拟波场提取方法--相关叠加法，并将其与奇异值分解法、预条件正则化共轭梯度法的虚拟波场提取结果进行分析对比，认为相关叠加法可以得到更加光滑的虚拟波场波形，且抗干扰能力较强.
2 MTEM基本原理MTEM方法通过接地电极向大地注入阶跃电流，同时采用接地电极接收瞬变电场轴向分量，最后通过解卷积的方式从采集到的信号中获取与地下介质信息相关的大地脉冲响应.
MTEM方法的数据采集方式如图 1所示.它采用电偶极源发射、电偶极子阵列接收的方式采集不同偏移距下的感应电场轴向分量，同时在发射源附近采集源信号.采用接收阵列采集数据，在一定程度上提高了数据采集效率，海上探测时尤为如此.同时这种采集方式能够对地下目标体进行多次覆盖，从而可将采集到的数据整理成共中心点道集、共偏移距道集，进而可借鉴地震勘探中的方法技术进行处理.
|
图 1 MTEM数据采集示意图 Fig. 1 Sketch of data acquisition in MTEM |
若将大地看成线性时不变系统，则多道瞬变电磁系统采集到的响应信号可以表达成如下形式(Wright et al., 2002)：
|
(1) |
其中，ak(xs, xr, t)表示接收到的总响应，i(k, xs, t)表示源电流，r(xr, t)表示记录系统的响应，k代表第k次数据采集，g(xs, xr, t)为大地脉冲响应，n(xr, t)为噪声，xs为发射源位置，xr为接收机位置，t为时间.根据(1)式可知，MTEM采集到的响应信号在经过去噪处理后，通过解卷积的方式可以获得大地脉冲响应.
为提高信噪比和分辨率，MTEM技术采用m序列发射波形(齐彦福等，2015).图 2为MTEM源电流信号、接收信号及大地脉冲响应示意图.
|
图 2 MTEM源电流信号、接收信号及大地脉冲响应示意图 (a)源电流；(b)实测响应；(c)大地脉冲响应. Fig. 2 Diagram of source current (a), receiving response (b), and earth impulse response (c) of MTEM |
虚拟波场提取是MTEM拟地震处理的前提.Lee等(1987)给出了时域扩散场与虚拟波场之间的数学对应关系，其具体表达式为
|
(2) |
式中，E(x, y, z, t)为时域扩散场，U(x, y, z, τ)为虚拟波场，τ为虚拟时间.通过求解公式(2)，即可完成MTEM数据的虚拟波场提取.
由式(2)可知，虚拟波场的求取是典型的反问题，其很重要的特征就是不适定性.对于反问题的求解，普通的数值化方法已不能满足要求(Xue et al., 2007).本文针对几种虚拟波场提取方法进行介绍，并对其结果进行对比，选择出最优方法.
3.2 虚拟波场提取方法 3.2.1 奇异值分解法公式(2)的离散数值积分形式为
|
(3) |
式中，
将公式(3)写成矩阵形式：
|
(4) |
式中，U=(U1, U2, …, Un)，E=(E1, E2, …, En)，Km×n为系数矩阵.将系数矩阵Km×n进行奇异值分解，
|
(5) |
式中，Lm×m，Vn×n为正交矩阵，Sm×n为对角矩阵.此时有
|
(6) |
奇异值分解后的系数矩阵，其病态程度在求解过程中不会增加.同时，奇异值分解将系数矩阵对角化，可以大大减少计算机运算时间.
3.2.2 预条件正则化共轭梯度法共轭梯度法是解大型非线性方程最有效的算法之一.目前常用的共轭梯度求解方法包括预条件共轭梯度法与正则化共轭梯度法.但本文所求方程系数矩阵条件数较大，单纯的预条件共轭梯度法或正则化共轭梯度法已不能满足要求，因此将两种方法结合起来形成预条件正则化共轭梯度法进行求解.
共轭梯度法要求系数矩阵为正定矩阵，为满足这一要求，将公式(4)转化为
|
(7) |
式中，Km×nT为Km×n的转置，只要Km×n为列满秩矩阵，Km×nTKm×n就是正定矩阵.但Km×nTKm×n的条件数较Km×n更大，因此首先使用超松弛预条件法对条件数进行降低.构造预条件子为
|
(8) |
式中，C1、C2、C3分别为Km×nTKm×n的对角元、下三角元和上三角元，ω为(0, 2)内的参数.
预条件子选定后，就可以通过正则化方法对方程组进行迭代求解.构造新的迭代方程如式(9)所示，
|
(9) |
式中，α为正则化参数，xk为第k次迭代的值，x初值选为
|
(10) |
为压制预条件正则化共轭梯度法提取结果中的虚假波峰，得到更加光滑的虚拟波场波形，采用相关叠加法进行虚拟波场提取.具体步骤为，
①对全时段MTEM数据(图 3)进行预条件正则化共轭梯度法虚拟波场提取，得到对应的虚拟波场Uall.
|
图 3 MTEM全时段数据示意图 Fig. 3 Whole time-domain data diagram of MTEM |
②从数据起始位置开始，选定一个时间段1(图 4)，对该时间段内的数据进行预条件正则化共轭梯度法虚拟波场提取，得到对应的虚拟波场U1.
|
图 4 时间段1数据示意图 Fig. 4 Data sketch diagram of time quantum 1 |
③将时间段1整体推移一个时间单位得到时间段2 (图 5)，对该时间段内的MTEM数据进行正则化共轭梯度法虚拟波场提取，得到对应的虚拟波场U2.
|
图 5 时间段2数据示意图 Fig. 5 Data sketch diagram of time quantum 2 |
再次将时间段整体推移一个时间单位得到时间段3，重复步骤③，如此循环下去，直到推移到最后一个时间段，得到各时间段对应的虚拟波场U3, U4, …, Un.
④依次将Un(n=1, 2, …, n-m+1)与Uall进行相关性分析，如果两者相关性大于某一阈值α，则将Un保留，否则舍去.两者的相关性定义如下：
|
⑤将全部保留的Un与Uall进行叠加，叠加结果作为最终的虚拟波场提取结果.
4 MTEM仿真数据虚拟波场提取效果比较 4.1 仿真模型为验证相关叠加法虚拟波场提取效果，采用三种提取方法对MTEM仿真数据进行虚拟波场提取，并对比其效果.图 6为设计的地电模型，为两层水平模型，地层分界面深度为200 m，上层电阻率为100 Ωm，下层电阻率为10 Ωm.发射源长度500 m，发射电流10 A，发射波形为12阶m序列波形，发射频率128 Hz，接收点位于距发射源中点1000 m处.
|
图 6 MTEM仿真数据模型 Fig. 6 Model of MTEM simulation data |
图 7为采用奇异值分解法进行虚拟波场提取的结果.由图 7可知，奇异值分解法得到的虚拟波场波形跳动比较大，存在比较明显的虚假波形.
|
图 7 奇异值分解法虚拟波场提取结果 Fig. 7 Extracting result of singular value decomposition method |
图 8为采用预条件正则化共轭梯度法进行虚拟波场提取的结果.由图 8可知，预条件正则化共轭梯度法得到的虚拟波场波形，波峰明显，仅存在幅值很小的虚假峰值.
|
图 8 预条件正则化共轭梯度法虚拟波场提取结果 Fig. 8 Extracting result of preconditioned regularized conjugate gradient method |
图 9为采用相关叠加法进行虚拟波场提取的结果.由图 9可知，相关叠加法虚拟波场提取的结果较正则化共轭梯度法更加光滑，虚假波峰很小.
|
图 9 相关叠加法虚拟波场提取结果 Fig. 9 Extracting result of correlation stacking method |
为了对比三种方法的抗干扰性，在仿真数据中分别加入5%、10%、20%、40%的随机噪声，采用三种方法分别进行虚拟波场提取，分析不同噪声背景下三种方法的虚拟波场提取效果.
图 10-13分别为5%、10%、20%、40%噪声环境下三种方法的虚拟波场提取效果.其中，图 10a、图 11a、图 12a、图 13a为采用奇异值分解法求得的虚拟波场；图 10b、图 11b、图 12b、图 13b为采用预条件正则化共轭梯度法求得的虚拟波场；图 10c、图 11c、图 12c、图 13c为采用相关叠加法求得的虚拟波场.由图 10-13可知，当噪声存在时，采用奇异值分解法提取的虚拟波场跳动非常剧烈，已不能得到稳定的提取结果；采用预条件正则化共轭梯度法提取的虚拟波场相对稳定，但是存在幅值较大的虚假波峰；采用相关叠加法提取的虚拟波场，虚假波峰幅值得到了明显压制，可以得到最优提取效果.
|
图 10 5%噪声背景虚拟波场提取结果 Fig. 10 Pseudo-wavefield extracting results with 5% noise |
|
图 11 10%噪声背景虚拟波场提取结果 Fig. 11 Pseudo-wavefield extracting results with 10% noise |
|
图 12 20%噪声背景虚拟波场提取结果 Fig. 12 Pseudo-wavefield extracting results with 20% noise |
|
图 13 40%噪声背景虚拟波场提取结果 Fig. 13 Pseudo-wavefield extracting results with 40% noise |
为对比不同虚拟波场提取方法的实际应用效果，采用三种方法分别对MTEM实测数据进行虚拟波场提取，并对其效果进行对比.
实验区位于内蒙古自治区乌兰察布市兴和县曹四夭村.图 14为MTEM测线图.测线全长4.8 km，角度为北偏东57°.发射极距240 m，接收极距40 m.发射电流30 A，频率128 Hz，采用10台接收机同时采集30道电场数据.采集时，接收机位置固定，发射机沿测线跑极.
|
图 14 MTEM测线图 Fig. 14 Survey line of MTEM |
本次分析所用数据采集点位于测线2900m处，发射机位于3480m处.采集数据经过预处理后得到的大地脉冲响应如图 15所示.
|
图 15 MTEM大地脉冲响应 Fig. 15 Earth impulse response of MTEM |
图 16为大地脉冲响应虚拟波场提取结果.其中，图 16a为奇异值分解法提取结果，图 16b为预条件正则化共轭梯度法提取结果，图 16c为相关叠加法提取结果.对比三种不同方法的虚拟波场提取结果可知，采用奇异值分解法提取的虚拟波场波形剧烈震荡，无法分辨波峰位置；采用预条件正则化共轭梯度法提取的虚拟波场波形较奇异值分解法有很大改善，但波形仍然呈现锯齿状，波峰不明显；采用相关叠加法提取的虚拟波场波形非常光滑，可以明显分辨出波峰位置.这表明，在对MTEM实测数据进行虚拟波场提取时，相关叠加法是最佳方法.
|
图 16 MTEM实测数据虚拟波场提取结果 Fig. 16 Pseudo-wavefield extracting result of MTEM measured data |
本文主要对MTEM大地脉冲响应虚拟波场提取方法进行了研究.首先采用两种已有的虚拟波场提取方法进行模拟数据虚拟波场提取，结果表明，奇异值分解法抗干扰能力较差，在噪声存在时无法得到有用的虚拟波场信息；预条件正则化共轭梯度法得到的虚拟波场波形比奇异值分解法更加光滑，但当噪声存在时得到的虚拟波场波形仍然存在比较明显的跳动.本文在预条件正则化共轭梯度法的基础上，提出采用相关叠加法进行虚拟波场提取，得到了更加光滑的虚拟波场波形，且抗干扰能力较强.实测MTEM数据虚拟波场提取结果表明，采用相关叠加法进行MTEM数据虚拟波场提取可以得到光滑且波峰明显的虚拟波场.
| Chen B C, Li J M, Zhou F T. 1999. Wave-field conversion method for transient electromagnetic field. Oil Geophysical Prospecting (in Chinese), 34(5): 539-545. | |
| Lee S, McMechan G A, Aiken C L V. 1987. Phase-field imaging: The electromagnetic equivalent of seismic migration. Geophysics, 52(5): 678-693. DOI:10.1190/1.1442335 | |
| Li H, Xue G Q, Zhou N N, et al. 2015. Appraisal of an array TEM method in detecting a mined-out area beneath a conductive layer. Pure and Applied Geophysics, 172(10): 2917-2929. DOI:10.1007/s00024-015-1075-0 | |
| Li X, Xue G Q, Song J P, et al. 2005. An optimize method for transient electromagnetic field-wave field conversion. Chinese J. Geophys. (in Chinese), 48(5): 1985-1990. | |
| Li X, Qi Z P, Xue G Q, et al. 2010. Three dimensional curve surface continuation image based on TEM pseudo wave-field. Chinese J. Geophys. (in Chinese), 53(12): 3005-3011. DOI:10.3969/j.issn.0001-5733.2010.12.025 | |
| Li X, Xue G Q. The Research of the Pseudo Seismic of the TEM Data (in Chinese).Beijing: Science Press, 2013. | |
| Qi Y F, Yin C C, Wang R, et al. 2015. Multi-transient EM full-time forward modeling and inversion of m-sequences. Chinese J. Geophys. (in Chinese), 58(7): 2566-2577. DOI:10.6038/cjg20150731 | |
| Qi Z P, Li X, Wu Q, et al. 2013. A new algorithm for full-time-domain wave-field transformation based on transient electromagnetic method. Chinese J. Geophys. (in Chinese), 56(10): 3581-3595. DOI:10.6038/cjg20131033 | |
| Wright D A, Ziolkowski A, Hobbs B A. 2001. Hydrocarbon detection with a multi-channel transient electromagnetic survey.//SEG Technical Program Expanded Abstracts 2001. Tulsa: SEG: 1435-1438. http://library.seg.org/doi/abs/10.1190/1.1816373 | |
| Wright D A, Ziolkowski A, Hobbs B A. 2002. Hydrocarbon detection and monitoring with a multicomponent transient electromagnetic (MTEM) survey. The Leading Edge, 21(9): 852-864. DOI:10.1190/1.1508954 | |
| Xue G Q, Li X, Guo W B, et al. 2006. Equivalent transformation from TEM field sounding data to plane-wave electromagnetic sounding data. Chinese J. Geophys. (in Chinese), 49(5): 1539-1545. | |
| Xue G Q, Yan Y J, Li X. 2007. Pseudo-seismic wavelet transformation of transient electromagnetic response in engineering geology exploration. Geophysical Research Letters, 34(16): 377-396. | |
| Xue G Q, Li X, Qi Z P, et al. 2011. Study of sharpening the TEM pseudo-seismic wave-form. Chinese J. Geophys. (in Chinese), 54(5): 1384-1390. DOI:10.3969/j.issn.0001-5733.2011.05.028 | |
| Xue G Q, Yan Y J, Li X. 2011. Control of the waveform dispersion effect and applications in a TEM imaging technique for identifying underground objects. Journal of Geophysics and Engineering, 8(2): 195-201. DOI:10.1088/1742-2132/8/2/007 | |
| Xue G Q, Bai C Y, Li X. 2012. Extracting the virtual reflected wavelet from TEM data based on regularizing method. Pure and Applied Geophysics, 169(7): 1269-1282. DOI:10.1007/s00024-011-0392-1 | |
| Ziolkowski A, Hobbs B A, Wright D A. 2007. Multitransient electromagnetic demonstration survey in France. Geophysics, 72(4): F197-F209. DOI:10.1190/1.2735802 | |
| 陈本池, 李金铭, 周凤桐. 1999. 瞬变电磁场的波场变换算法. 石油地球物理勘探, 34(5): 539–545. | |
| 李貅, 薛国强, 宋建平, 等. 2005. 从瞬变电磁场到波场的优化算法. 地球物理学报, 48(5): 1185–1190. | |
| 李貅, 戚志鹏, 薛国强, 等. 2010. 瞬变电磁虚拟波场的三维曲面延拓成像. 地球物理学报, 53(12): 3005–3011. DOI:10.3969/j.issn.0001-5733.2010.12.025 | |
| 李貅, 薛国强. 瞬变电磁法拟地震偏移成像研究.北京: 科学出版社, 2013. | |
| 齐彦福, 殷长春, 王若, 等. 2015. 多通道瞬变电磁m序列全时正演模拟与反演. 地球物理学报, 58(7): 2566–2577. DOI:10.6038/cjg20150731 | |
| 戚志鹏, 李貅, 吴琼, 等. 2013. 从瞬变电磁扩散场到拟地震波场的全时域反变换算法. 地球物理学报, 56(10): 3581–3595. DOI:10.6038/cjg20131033 | |
| 薛国强, 李貅, 郭文波, 等. 2006. 从瞬变电磁测深数据到平面电磁波场数据的等效转换. 地球物理学报, 49(5): 1539–1545. | |
| 薛国强, 李貅, 戚志鹏, 等. 2011. 瞬变电磁拟地震子波宽度压缩研究. 地球物理学报, 54(5): 1384–1390. DOI:10.3969/j.issn.0001-5733.2011.05.028 | |