2. 中国矿业大学(北京)地球科学与测绘工程学院, 北京 100083;
3. 中国科学院矿产资源研究重点实验室, 中国科学院地质与地球物理研究所, 北京 100029;
4. 煤炭科学研究总院矿山安全技术研究分院, 北京 100013
2. College of Geoscience and Surveying Engineering, China University of Mining and Technology, Beijing 100083, China;
3. Key Laboratory of Mineral Resources, Institute of Geology and Geophysics, Chinese Academy of Sciences, Beijing 100029, China;
4. Mine Safety Technology Branch of China Coal Research Institute, Beijing 100013, China
近几年来,矿井瞬变电磁法在矿井工作面顶板、底板岩层富水性探测和掘进巷道超前探测中发挥了重要作用(程久龙等,2014).由于矿井瞬变电磁法勘探是体积勘探,受体积效应的制约,其分辨率难以满足生产实际要求(Cheng et al.,2015),并且受井下施工环境和金属物体的影响,采集到的数据中往往存在一定的噪声,致使信噪比降低,而且,矿井瞬变电磁法探测成果大多是以电阻率等值线断面图的形式表达,不能准确地分辨实际地质情况中最为重要的地层电性界面.基于波场变换的瞬变电磁拟地震成像技术因能够相对突出异常体的物性和几何参数而越来越受到人们的重视(李貅和薛国强,2013),然而由于巷道超前探测方向不可能很密,加上采样点相对稀疏等问题,所 以波场变换成像的信噪比和分辨率有待进一步提 高,合成孔径雷达及航空瞬变电磁中合成孔径成像的思想为矿井瞬变电磁法合成孔径成像提供了新的思路.
合成孔径成像技术最早是应用于雷达方面,国内外学者做了一系列的研究,研究内容主要包括算法、应用及成像技术的对比等(Cetin and Karl,2001;王治乐等,2004;Samadi et al.,2009),应用领域较为广泛.在超声波探测方面,Jensen等(2006)利用合成孔径成像技术来提高医学超声波检测图的成像分辨率,Ganguli等(2012)研究了混凝土超声波探伤的合成孔径聚焦成像技术.在电磁法探测方面,Mason等(2001)研究了利用宽频井中雷达合成孔径干涉法来绘制3D矿体图,张春城和周正欧(2004)研究了基于stolt偏移的探地雷达合成孔径成像.在全空间瞬变电磁数值模拟方面,李展辉和黄清华(2014)研究了将复频率参数完全匹配层吸收边界应用到均匀全空间模型和包含巷道及异常体的全空间模型.上述研究没有涉及异常体界面的定位.
在瞬变电磁勘探方面,国内学者鉴于合成孔径雷达和航空瞬变电磁之间的相似性首先把合成孔径成像的思想应用于航空瞬变电磁,研究了利用合成孔径方法对航空瞬变电磁数据进行成像,以提高分辨率,增加勘探深度(李貅等,2010).之所以可以把合成孔径雷达的思想引入到地面瞬变电磁数据处理中,是因为在不同测点所接收到的由同一个发射源所激励的二次场数据等效于合成孔径雷达单点发射源时的阵列数据,而地面瞬变电磁所有的发射源则构成了动源观测.由于井下瞬变电磁法探测装置及数据采集方式与地面基本上相似,所以,像地面瞬变电磁合成孔径成像方法一样,对井下瞬变电磁数据做类似的处理是可行的.本文利用合成孔径成像技术对矿井瞬变电磁波场变换后的拟地震信号进行处理,结合超前探测的工程应用实例对方法效果进行对比分析和验证,实现突出弱异常且提高界面分辨率的目的.
2 波场变换原理瞬变电磁场与其相对应的虚拟波场之间的数学积分形式为(Lee et al.,1989)
其中P(r,t)为瞬变电磁场,U(r,q)为待求拟地震 波场,q为与时间t相对应的类时间变量,单位为
将公式(1)改写为数值积分形式为(李貅和薛国强,2013)
其中m为采样个数,n为数值积分时积分区间剖分的个数,hj为积分步长.
将公式(2)写成矩阵形式为
其中 U =(u1,…,uj,…,un)T为离散化的虚拟子波, P =(p1,…,pi,…,pm)T为离散化的瞬变电磁场,为m×n阶矩阵并且m≥n,求解方程组(3)就可得出对应瞬变电磁场 P 的拟地震波场 U .
需要注意的是方程组(3)中的系数矩阵 A 常呈病态.为了获得稳定的解,可采用阻尼法最小二乘法(杨文采,1997),建立目标函数为
其中ε2为阻尼因子,令公式(4)的梯度为0可得
进而有
其中(A T A +ε2 I)+为(A T A +ε2 I)的Moore-Penrose(M-P)广义逆.对于广义逆的求解,可采用IMqrginv 算法(Ataei,2014).阻尼因子需在解估计的分辨率和方差间折衷并通过实验选取(杨文采,1997; 陈本池等,1999).矿井瞬变电磁超前探测数据是在同一掘进工作面处不同方向进行采集的,在数据采集 段内,误差水平是一定的,所以对同一个掘进工作面处的超前探测数据,其波场变换可采用同一个阻尼因子.
波场变换时,瞬变电磁场中的误差会给变换后的波场造成较大的影响,特别是当解估计的单位协方差矩阵(杨文采,1997)较大时,这种影响会更大.为此,在波场变换前需对所采集到的数据进行滤波处理,根据实验选取小波变换法(To et al.,2009)进行滤波效果较好.
由于矿井瞬变电磁超前探测位置是在矿井巷道掘进工作面处,故采集的信号是前后双向回波,可借助于半空间和全空间瞬变响应关系将全空间瞬变响应校正为半空间瞬变响应,也即相当于仅有超前探测时的前方回波信号.半空间和全空间瞬变响应磁场间的关系(Kaufman and Eaton,2001)为
其中BzW为全空间瞬变磁场响应,BzH为半空间瞬变磁场响应,k为全空间瞬变响应磁场分量和半空间 瞬变响应磁场分量间的比例系数.Kaufman和Eaton(2001)研究了在均匀全空间介质中,对于中心回线或重叠回线装置,当τ/r>15时,公式(7)中系数k为5/2,其中r为发射回线的半径,,t为采样延迟时间,σ为介质电导率,μ为介质磁导率,可近似取为真空中的磁导率μ0.由于矿井瞬变电磁的施工环境位于煤矿巷道中,围岩介质为煤岩地层,其电阻率一般大于20 Ωm,因此,对于边长为2 m×2 m的发射线圈,在矿井瞬变电磁观测时间段内,τ/r>15,所以公式(7)中k取值为5/2. 3 矿井瞬变电磁合成孔径成像方法 3.1 雷达和地面瞬变电磁合成孔径成像原理
合成孔径雷达(SAR)就是利用雷达-目标的相对运动,以单元阵在不同相对空间位置上接收的时间采样信号替代大尺寸源的空间采样信号.由雷达-目标相对运动形成的轨迹构成一个合成孔径来代替庞大的真实孔径,把尺寸较小的真实天线孔径用数据处理的方法合成一较大的等效天线孔径的雷达.合成孔径雷达的特点是分辨率高,能有效地识别目标物(Chan and Koo,2008),合成孔径成像原理如图 1a所示.
地面瞬变电磁法是在地表布设一个矩形回线发射框,在发射框的内部和外部接收瞬变电磁信号.地面瞬变电磁法也是在移动方式中观测瞬变电磁信号,并且不同发射源间的信号满足类相干性.因此地面瞬变电磁法也可以使用合成孔径成像技术(Guo et al.,2012;李貅等,2012).地面瞬变电磁合成孔径原理如图 1b所示.
地面瞬变电磁合成孔径成像公式为(李貅等,2012)
其中U(ri,t)表示i点的波场值,ri为i点到-N,…,N 内某点的距离,τ为相对时移量,2N+1为合成孔径的长度.τkm为对应于最大归一化互相关系数ρmax(rk,τkm)时的时移量,称为最佳时移量,U′(ri,tj)表示合成孔径成像后的值.
从以上雷达和地面瞬变电磁合成孔径成像技术中可以看出,合成孔径成像方法的本质是不同测点信号间的相关叠加.
3.2 矿井瞬变电磁超前探测合成孔径成像原理和方法矿井瞬变电磁法超前探测的地点是在巷道掘进工作面处,如图 2a所示,图中矩形小回线源紧靠掘进工作面掌子面上,虚线圆圈代表掘进工作面前方瞬变电磁一次场不同时刻的“烟圈”.常采用扇形布置,即在掘进工作面位置放置小线框天线,以巷道正前方为中心轴,按一定的角度(10°或15°)往两帮分别进行不同方向探测,最终形成扇形探测断面,如图 2b所示,图中箭头表示同一断面上不同的探测方向(发射线框的法线方向),虚线框表示掘进工作面前方的异常体,设掘进工作面正前方为z轴,z轴所指方向为0°,用与z轴所夹的角度θ来定义探测方向(z轴左侧的角度为负,右侧的角度为正),即:U(-θi)表示z轴左侧与z轴成θi角度方向上的实测信号值,U(+θi)表示z轴右侧与z轴成θi角度方向上的实测信号值,U(θ0)表示z轴方向上的实测信号值.
从图 2矿井瞬变电磁超前探测的观测方式中可以看出:同一个掘进工作面处其周围全空间范围内的地质条件在观测期间是一定的,不同探测方向(也即不同发射源,相当于地面瞬变电磁法的不同测点)观测到的瞬变电磁信号间具有类相干性.因此,矿井瞬变电磁数据具备相干叠加的条件,也即可以采用合成孔径成像技术.
在矿井瞬变电磁拟地震波场变换、反褶积和时深转换等的基础上(程久龙等,2013),根据合成孔径成像的思想,本文提出合成孔径成像公式为
其中U′(θm+N,tk)为合成孔径成像后θm+N探测方向上的信号值,θm是指合成孔径成像的起始探测方向,而θm+2N为合成孔径成像的结束探测方向,m和m+2N分别为起始探测方向编号和结束探测方向编号,2N+1为合成孔径的长度,M是指同一探测方向上重复测量次数,Wij是权重函数,U(θi,tj)是 矿井瞬变电磁波场变换后θi探测方向上第j次测 量的拟地震波振幅值,τijk为权函数Wij取极值时的时移量.
权函数Wij为相关系数,表达式为
其中i依次取m、m+1、…、m+2N,j取1、2、…M,τijk为公式(10)取极大值时的相对时移量,K为波场变换后某一探测方向上的类时间门数.矿井瞬变电磁超前探测合成孔径成像原理如图 3所示.
相关系数表征了两个探测方向上波场变换后信号值间的相关程度.相关系数越接近于1或-1,相关性越强;反之,相关系数越接近于0,相关性越弱.在通常情况下,通过以下取值范围判断信号间的相关性(表 1).
某矿巷道掘进至G22+5 m位置,需要查清掘 进迎头前方是否存在含水地质体,对迎头前方100 m 范围内砂岩富水情况、断层带富水情况进行总体控制,对含水异常区进行地质解释.选择矿井瞬变电磁法进行超前探测,装置方式为中心回线,回线边长为2 m×2 m,电流1.8 A,关断时间约为110 μs,实测数据的时间区间为[116.8 μs,7088 μs],均在关断时间之外.巷道围岩为砂岩,电阻率大于20 Ωm.
图 4是G22+5 m探测点水平断面11个方向(以巷道掘进方向为0°,往左帮每隔10°一个探测方向,到左帮50°,往右帮每隔10°一个探测方向,到右帮50°)的实测感应电动势曲线图,图中曲线表示不同探测方向上相同时间门间的连线.
对图 4实测数据分别进行滤波和全空间校正后,进行波场转换和时深转换等数据处理(程久龙等,2013),可以得到该测点水平断面上11个探测方向(从左帮-50°至右帮+50°,每隔10°依次探测)的拟地震信号,横轴表示深度,纵轴代表探测方向 角,如图 5所示.从图中可以看出,拟地震信号有明显的两个同相轴,但总体偏弱.第一个相位在30~40 m范围内,能量相对较强,反映的是高阻到低阻的界面特征;第二个相位在65~80 m范围内,能量相对较弱,反映的是低阻到高阻的界面特征(程久龙等,2013).
根据公式(10),对水平断面11个探测方向的拟地震数据分别计算相邻两探测方向之间的相关系数,判断相关性,见表 2.从表中可以看出,水平断面上相邻两测点数据之间最小的相关系数为0.12,最大的相关系数为0.9726.
孔径长度的选择,即中心探测方向确定以后,具 体应该选取几个探测方向上的数据进行叠加而且叠加后会有好的效果,一般情况下应通过试验来选取.根据矿井瞬变电磁超前探测的特点,水平断面上只有11个方向的探测数据,探测方向数较少,选择了3个方向和4个方向上的数据分别进行合成孔径试验,最终确定选择3个方向进行合成孔径成像.
图 6为巷道迎头G22+5m位置处矿井瞬变电磁超前探测合成孔径成像后拟地震信号,可以看出,波场变换曲线的极值位置没有改变,深度方向未发生偏移,即合成孔径处理不影响波场曲线极值点与深度的对应关系.
对比该断面合成孔径前(图 5)与合成孔径后(图 6)拟地震信号波形图,合成孔径之后的波形图在相关性较好的探测方向较合成孔径之前的信号得到了明显加强,如在相关性较好的-20°探测方向到+20°探测方向间的探测数据,合成孔径前的两处曲线幅值较为平缓,但合成孔径后,波形幅值增大,信号得到了显著增强,并且波形宽度未发生变化;而在相关性相对较差的-50°探测方向以及50°探测方向,合成孔径前曲线幅值较为平缓,在合成孔径后,波形幅值被削弱;其他探测方向合成孔径后波场曲线幅值得到了不同程度的加强.整体而言,各个探测方向上有意义的弱异常得到了加强,信噪比得到了提高,合成孔径效果比较显著.
由以上分析可以看出:合成孔径处理不影响拟地震波场曲线极值点与深度的对应关系.相关性较强的两个探测方向上,其拟地震信号极值点在合成孔径后得到了明显增强;而相关性较弱的两个探测方向上,其拟地震信号极值点在合成孔径后被削弱.这说明合成孔径成像在矿井瞬变电磁超前探测中突出弱异常、压制假异常进而提高分辨率的有效性.
将图 6的拟地震波场曲线进行归一化处理,并辅以实际断面扇形坐标,标定电性界面位置,以利于解释.图 7是G22+5 m位置瞬变电磁超前探测合成孔径成像断面图与常规视电阻率断面图对比,可以看出:常规视电阻率断面图上有一明显的低阻异常,位于巷道前方约20~60 m范围,中心在34 m附近;合成孔径成像断面图有两个电性界面,第一电性界面位于前方35 m处,信号能量较强,应该是高 阻到低阻的突变界面;第二个电性界面位于前方73 m处,信号能量相对较弱,应该是低阻到高阻的突变界面.两处界面大体上为直线,并与探测方向近垂直.巷道掘进实际揭露了探测位置为前方33~70 m范围砂岩含水,为弱富水,与合成孔径成像解释成果吻合较好,而常规视电阻率断面图所反映的异常与实际揭露存在一定的偏差.这进一步说明,仅依据断面图视电阻率的变化来确定电性界面位置存在较大误差,而合成孔径成像确定电性界面具有明显优势.
论文在将矿井瞬变电磁波场转换为拟地震波场的基础之上,实现了矿井瞬变电磁合成孔径成像,并结合工程应用实例对方法的有效性和实用性进行了探讨和验证,主要结论如下:
(1)通过对比分析合成孔径前后波形图可知,合成孔径之后的波形图在相关性较好的探测方向上较合成之前的信号得到了加强,且波形宽度未发生变化;而在相关性相对较差的方向上合成孔径后较合成孔径前由噪声引起的小幅值异常进一步减小,降低了随机噪声,提高了信噪比.
(2)合成孔径成像的效果与各探测方向之间的相关系数或相关性密切相关,较大的相关系数会获得较好的合成孔径成像效果.
(3)对实际资料的处理和分析可知,低阻异常体在常规的视电阻率等值线图中表现为一个放大的低阻区域,而在合成孔径成像后的成果图上则表现为更为精确的两个电性界面.合成孔径成像处理后电性界面的分辨率显著提高.一般来说,在高阻到低阻的突变界面处,波场变换后信号较强,对应界面合成孔径成像后的能量较强;而在低阻到高阻的突变界面处,波场变换后信号较弱,对应界面合成孔径成像后的能量也相对较弱.对于后者,增强波场变换后信号有待进一步研究.
(4)合成孔径成像方法能够突出瞬变电磁数据中所包含的电性界面信息,弥补了矿井瞬变电磁法在异常体深度确定方面的不足,提高信噪比,突出弱异常进而提高分辨率及勘探精度,发展了矿井瞬变电磁法对低阻异常电性界面的精细探测.
[1] | Ataei A. 2014. Improved Qrginv algorithm for computing moore-penrose inverse matrices. ISRN Applied Mathematics, 2014:Article ID 641706. |
[2] | Cetin M, Karl W C. 2001. Feature-enhanced synthetic aperture radar image formation based on nonquadratic regularization. IEEE Transactionson Image Processing, 10(4):623-631. |
[3] | Chan Y K, Koo V C. 2008. An introduction to synthetic aperture Radar(SAR). Progress in Electromagnetics Research B, 2:27-60. |
[4] | Chen B C, Li J M, Zhou F T. 1999. Wave-field conversion method for transient electromagnetic field. OGP (in Chinese), 34(5):539-545. |
[5] | Cheng J L, Li F, Peng S P, et al. 2014. Research progress and development direction on advanced detection in mine roadway working face using geophysical methods. Journal of China Coal Society (in Chinese), 39(8):1742-1750. |
[6] | Cheng J L, Li F, Peng S P, et al. 2015. Joint inversion of TEM and DC in roadway advanced detection based on particle swarm optimization. Journal of Applied geophysics,123:30-35. |
[7] | Cheng J L, Qiu H, Ye Y T, et al, 2013. Research on wave-field transformation and data processing of the mine transient electromagnetic method. Journal of China Coal Society (in Chinese), 38(9):1646-1650. |
[8] | Ganguli A, Rappaport C M, Abramo M, et al. 2012. Synthetic aperture imaging for flaw detection in a concrete medium. NDT&E International, 45(1):79-90. |
[9] | Guo W B, Xue G Q, Li X, et al. 2012. Correlation analysis and imaging technique of TEM data. Exploration Geophysics, 43(3):137-148. |
[10] | Jensen J A, Nikolov S I, Gammelmark K L, et al. 2006. Synthetic aperture ultrasound imaging. Ultrasonics, 44:e5-e16. |
[11] | Kaufman A A, Eaton P A. 2001. The Theory of Inductive Prospecting. New York:Elsevier,365-373. |
[12] | Lee K H, Liu G, Morrison H F. 1989. A new approach to modeling the electromagnetic response of conductive media. Geophysics, 54(9):1180-1192. |
[13] | Li X, Qi Z P, Liu Y A, et al. 2010. A research of synthetic aperture rapid imaging of ATEM. //26rd CGS, Chinese Geophysics 2010 (in Chinese). 639. |
[14] | Li X, Xue G Q. 2013. Study on Pseudo-Seismic Migration Imaging of Transient Electromagnetic Method. Beijing:Science Press. |
[15] | Li X, Xue G Q, Liu Y A, et al. 2012. A research on TEM imaging method based on synthetic-aperture technology. Chinese J. Geophys. (in Chinese), 55(1):333-340, doi:10.6038/j.issn.0001-5733.2012.01.034. |
[16] | 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. |
[17] | Mason I, Osman N, Liu Q, et al. 2001. Broadband synthetic aperture borehole radar interferometry. Journal of Applied Geophysics,47(3-4):299-308. |
[18] | Qi Z P. 2011. 3D continuation imaging technology with a synthetic aperature based on the transient electromagnetic method (in Chinese). Xi'an:Chang'an University. |
[19] | Samadi S, Çetin M, Masnadi-Shirazi M A. 2009. Multiple feature-enhanced synthetic aperture radar imaging.//Proc. SPIE Algorithms Synthetic Aperture Radar Imagery, XVI.Orlando, Florida, USA. |
[20] | To A C, Moore J R, Glaser S D. 2009. Wavelet denoising techniques with applications to experimental geophysical data. Signal Processing, 89(1):144-160. |
[21] | Wang Z L, Zhou Y P, Zhang W, et al. 2004. Synthetic aperture imaging techniques. Journal of Harbin Institute of Technology (in Chinese), 36(3):384-387. |
[22] | Yang H Y, Deng J Z, Zhang H, et al. 2010. Research on full-space apparent resistivity interpretation technique in mine transient electromagnetic method. Chinese J. Geophys.(in Chinese), 53(3):651-656, doi:10.3969/j.issn.0001-5733.2010.03.020. |
[23] | Yang W C. 1997. Theory and Methods of Geophysical Inversion (in Chinese). Beijing:Geological Publishing House. 45-67. |
[24] | Zhang C C, Zhou Z O. 2004. Ground penetrating radar synthetic aperture imaging based on Stolt migration. Chinese Journal of Radio Science (in Chinese), 19(3):316-320. |
[25] | 陈本池, 李金铭, 周凤桐. 1999. 瞬变电磁场的波场变换算法. 石油地球物理勘探, 34(5):539-545. |
[26] | 程久龙, 李飞, 彭苏萍等. 2014. 矿井巷道地球物理方法超前探测研究进展与展望. 煤炭学报, 39(8):1742-1750. |
[27] | 程久龙, 邱浩, 叶云涛等. 2013. 矿井瞬变电磁法波场变换与数据处理方法研究. 煤炭学报, 38(9):1646-1650. |
[28] | 李貅, 戚志鹏, 刘银爱等. 2010. 航空瞬变电磁合成孔径快速成像方法研究.//中国地球物理2010-中国地球物理学会第二十六届年会论文集, 639. |
[29] | 李貅, 薛国强. 2013. 瞬变电磁法拟地震偏移成像研究. 北京:科学出版社. |
[30] | 李貅, 薛国强, 刘银爱等. 2012. 瞬变电磁合成孔径成像方法研究. 地球物理学报, 55(1):333-340, doi:10.6038/j.issn.0001-5733.2012.01.034. |
[31] | 李展辉, 黄清华. 2014. 复频率参数完全匹配层吸收边界在瞬变电磁法正演中的应用. 地球物理学报, 57(4):1292-1299, doi:10.6038/cjg20140426. |
[32] | 戚志鹏. 2011. 瞬变电磁三维合成孔径延拓成像方法研究. 西安:长安大学. |
[33] | 王治乐, 周彦平, 张伟等. 2004. 合成孔径成像技术对比研究. 哈尔滨工业大学学报, 36(3):384-387. |
[34] | 杨海燕, 邓居智, 张华等. 2010.矿井瞬变电磁法全空间视电阻率解释方法研究. 地球物理学报, 53(3):651-656, doi:10.3969/j.issn.0001-5733.2010.03.020. |
[35] | 杨文采. 1997. 地球物理反演的理论与方法. 北京:地质出版社,45-67. |
[36] | 张春城, 周正欧. 2004. 基于Stolt偏移的探地雷达合成孔径成像研究. 电波科学学报, 19(3):316-320. |