刘彦涛,陈伟,张文峰,程莎莎,骆燕
(1.核工业航测遥感中心,河北 石家庄 050002;2.中核集团铀资源地球物理勘查技术(重点实验室),河北 石家庄 050002;3.河北省航空探测与遥感技术重点实验室,河北 石家庄 050002)
航空大地电磁法(AFMAG)是一种被动源频率域电磁勘探方法,具有探测深度大、工作效率高、数据采集受地形影响小的优点,近几年得到了快速的发展。与地面大地电磁原理相同,由远方雷电、太阳风等产生的电磁信号在一定范围内可视为均匀平面电磁波[1-2],该变化的平面电磁波垂直入射到地下,与地下不同岩性、构造产生不同电磁信号,通过测量这些微弱的电磁信号,从而达到勘查地下地质体的目的。
1972 年Vozoff Keeva 对大地电磁倾子在二维条件下响应特征及原理进行了系统的阐述,但受限于当时设备及数据噪声处理技术,野外实测数据中磁场Hz分量受噪声影响严重,不能进行倾子参数换算用于实际数据解释[3]。1985年Labson 等人提出远参考噪声压制技术[4],才使得倾子在实际数据解释中得到成功应用。
航空大地电磁方法通过飞行载体在空中进行数据采集,相对于传统地面大地电磁法,不能采集电场强度分量,通过采集磁感应强度x、y、z3 个分量,换算倾子参数[5-6],倾子是航空大地电磁数据解释中的重要参数。以往对倾子参数响应特点研究较多,但多数为数据模拟研究,且算法多为二维算法[7-9]。本文通过COMSOL 软件进行数值模拟,求取倾子总散度,分析倾子响应特征,并通过航空大地电磁ZTEM 系统实测数据进一步对倾子响应特征进行分析,探讨倾子资料解释的优缺点。
自然界产生的电磁场具有水平极化的特点,当地下地质体电阻率存在横向变化时,就会激发出垂直磁场分量。地下地质体所激发的磁场垂直分量Hz与两个水平分量Hx、Hy存在复线性关系[3],这种复线性关系可以用两个复数系数在任意频率表示:
式(1)中Hx(ω)、Hy(ω)、Hz(ω)分别为频率域中磁场的x、y、z分量,Tx(ω)、Ty(ω)为倾子系数,倾子是一复数值且随频率变化。
倾子参数与地下地质体电性分布相关而不与激发电磁波极化方向相关,所以倾子可以作为探测地质体的参数。
式(1)中Hx(ω)、Hy(ω)、Hz(ω)为野外实测时间域磁场经傅里叶变换后的值。为求取Tx、Ty值,方程式(1)分别乘以Hx(ω)、Hy(ω)的共轭复数值[10],构建方程组:
式(2)中,[Hz H*x]表示Hz与Hx共轭复数乘积的累加和,其他表示类似。
通过方程组(2)可求取倾子值Tx、Ty:
ZTEM 系统由加拿大Geotech 公司研发,是目前较为成熟的航空大地电磁测量系统,自2007 年投入使用以来,在全球范围内进行了大量的试验与找矿探测,取得了较好的探测效果[11-13]。ZTEM 系统由空中接收系统和地面基站接收系统两部分组成,如图1 所示。
图1 ZTEM 系统实拍Fig.1 Photos of ZTEM system
空中接收系统,接收线圈直径7.4 m,有效面积1 200 000 m2,吊挂在直升机下方90 m 处,线圈平面平行于地面,接收磁场垂直分量,设计线圈离地高度100 m。线圈上安装有3 个GPS 定位天线,用于监控线圈姿态。实际飞行数据采集过程中,由于风速等外界不可控因素影响,线圈不是完全水平,接收到的Hz分量耦合有Hx、Hy分量,可通过3 个GPS 天线校正线圈姿态误差。
地面基站接收系统主要由两个磁探头组成,磁探头有效面积660 000 m2。同一测区范围内,磁场水平分量Hx、Hy变化很小,测区内选定干扰较小位置作为地面基站放置点,两个磁探头水平放置,相互垂直,分别采集Hx、Hy分量。
对采集得到的Hz分量与Hx、Hy分量磁场数据进行傅里叶变换,提取25、37、75、150、300、600 Hz 频率域信息,通过式(3)求取这6 个频率沿测线倾子参数Tx和垂直测线倾子参数Ty,通过分析解释倾子参数及相关衍生值从而达到探测地下地质体的目的。
使用基于有限元的三维正演程序对倾子响应进行数值模拟,分析倾子响应特征规律,为倾子资料解释提供依据。
为求解倾子值需要求解两线性相关极化的磁场值。将大地电磁场源分解为两正交场源S1和S2的等效作用,为便于分析,假设这两个场源极化方向分别沿主轴方向,源1沿y轴方向,设为S1,源2沿x轴方向,设为S2。在源S1激励下,求取磁场3个分量在源S2激励下,求取磁场3个分量,利用式(4)求取倾子参数[2]。
模拟中综合使用磁绝缘和理想磁导体边界条件,使场源平行于磁绝缘边界,垂直于理想磁导体边界,从而达到模拟无限域的效果[14-15]。使用自由四面体对模型进行网格剖分,对异常体网格加密,提高求解精度,剖分效果如图2 所示。
图2 模型网格剖分Fig.2 Diagram of model meshing division
数值模拟中使用25、37、75、150、300、600 Hz共6个频率,模型背景电阻率100 Ω·m,充分考虑模拟边界效应以及趋肤深度的影响,构建如图3所示“回”型异常体,模型长2 000 m,宽2 000 m,高1 000 m。外层异常体电阻率10 Ω·m,内层异常体电阻率1 000 Ω·m,x轴沿异常体长轴方向。本文主要讨论倾子响应平面特征,为突出异常响应特点,异常体上方不设置覆盖层。
图3 模型尺寸图示Fig.3 Model size diagram
研究表明[16-21],倾子资料实部和虚部分量响应特征类似,且倾子实部分量对异常体反应更好。实际生产应用中只对实部倾子数据成图,只有在反演过程中才使用到虚部倾子资料,所以文中只对倾子实部资料进行讨论分析。
图4为沿y轴方向模型中线电阻率曲线图,图4a 为ρxy视电阻率曲线,图4b 为ρyx视电阻率曲线,图中虚线为不同电阻率介质界线。从图中可以看出,两个视电阻率曲线对于电阻率突变界面都有较好的反应,其中ρyx视电阻率曲线对电阻率突变界面反应更好。ρxy和ρyx视电阻率参数对背景围岩和低阻体电阻率值反应较为准确,但对模型中间高阻体电阻率值反应存在偏差,且ρxy视电阻率曲线对中间高阻体的反应极不明显,反映了大地电磁法视电阻率参数对低阻体响应优于对高阻体响应的特点。
图4 y 轴方向视电阻率曲线Fig.4 Apparent resistivity curve along the y axis
图5和图6 分别为不同频率倾子Tx实部和倾子Ty实部平面图。其中倾子Tx只对x轴方向电阻率变化有响应,倾子Ty只对y轴方向电阻率变化有响应。倾子实部在横向电阻率交界位置处形成异常中心,当电阻率值从高阻到低阻的横向分界面处,倾子形成正异常中心;当电阻率值从低阻到高阻的横向分界面处,倾子形成负异常中心。
图5 不同频率倾子Tx实部平面图Fig.5 Plan map of real part of tipper Tx at different frequencies
图6 不同频率倾子Ty实部平面图Fig.6 Plan map of real part of tipper Ty at different frequencies
不同频率倾子正演结果对于电阻率横向变化边界皆有响应。通过对比不同频率倾子响应平面图可以看出,频率越高,倾子值正负差异越大,对比越大,对横向电性分界面的分辨能力增强。探测频率从高频变化到低频时,探测深度增加,体积效应增强,异常中心范围变大,探测频率较高时,体积效应小,对横向电性界面反应更为精确。
从上述结果可以看出,倾子对于横向电性分界面响应较好,并在边界处形成相应的极值异常区,实际解释中希望通过结合倾子Tx和倾子Ty参数圈定异常体,并在异常体上方形成相应的极值异常区,为此引入倾子总散度参数TD(Total Derivative),TD是通过计算沿测线方向倾子Tx对x的导数与垂直测线方向倾子Ty对y的导数之和求得[22],具体计算公式如下:
图7为不同频率倾子总散度实部平面图,从图中可以看出倾子总散度对异常体圈定效果较好,能有效地显示高阻与低阻异常体位置,对于模型中间的高阻体边界圈定也相当准确,但总散度平面图不能反应异常体实际电阻率值,只能显示相对差异。
图7 不同频率倾子总散度实部平面图Fig.7 Plan map of real part of total divergence of tipper at different frequencies
对比不同频率总散度平面图可以看出,当频率为25 Hz 时,图中数据范围为-0.005~0.005;当频率为600 Hz 时,图中数据范围为-0.022~0.022。即高频数据范围大于低频数据范围,相应高频数据对比更大,异常体反应更加明显。实际上DT参数反应沿相应方向的变化率,对于低频倾子响应,体积效应大,倾子极值对横向电阻率突变界面圈定范围大,变化率小,相应总散度值偏低,而高频倾子响应恰恰相反,对于高频倾子响应,体积效应小,倾子极值对横向电阻率突变界面圈定范围小,变化率大,相应总散度值偏高。
2015 年核工业航测遥感中心在某试验区开展了航空电磁测量,利用航空瞬变电磁VTEM 系统和航空大地电磁ZTEM 系统对测区进行探测,查明地质构造,推断成矿靶区。截取一部分测区数据,将航空大地电磁倾子总散度成像效果与航空瞬变电磁和航空磁法探测效果对比,说明倾子总散度在实际探测中的应用效果。
图8为不同频率倾子总散度实部平面图,6幅分图中均能反映出两条南西-北东向低值异常条带,推断为岩石破碎带充水形成的低阻断裂带。6 个频率对两条低值异常带的反应特征类似,对比各图数据值范围,高频相对低频数据范围更大,数据差异性更明显,相应探测效果更好,与正演模拟的规律相同。
图9为试验区地质示意图。试验区北部岩性主要为奥长花岗岩,南部岩性主要为太古宇龙岗群杨家店组下段和太古宇龙岗群四道砬子河组上段。区域构造方向为北东-南西向,小型断裂发育,东双丫倒转背斜东北端位于试验区中部。将图8 中航空大地电磁推测断裂投影到地质图中,分别命名为F1 断裂,F2 断裂。对比推测断裂位置与实际地质图,发现实际地质图中并未显示F1 与F2 位置存在断裂构造,推断F1 与F2 为本次测量发现的隐伏断裂。
图8 实测不同频率倾子总散度实部平面图Fig.8 Plan map of real part of total divergence of field tipper at different frequencies
图9 研究区地质示意图Fig.9 Geological sketch of the study area
为验证上述推断,对比分析航空瞬变电磁和航磁探测结果,对上述推断进行验证说明。图10 为航空大地电磁倾子总散度与航空瞬变电磁和航空磁法探测效果对比图。图10a 为75 Hz 总散度图,根据图中低阻条带延伸标记两条断裂F1,F2。图10b 为航空瞬变电磁时间常数平面图,图中F2 断裂延伸方向时间常数表现为高值条带,反应低阻异常,与倾子总散度断裂F2 异常特征相吻合;图10b 中F1 断裂延伸效果不明显,实际航空瞬变电磁探测深度较浅(一般地质条件下小于600 m),推断F1 断裂埋深超过航空瞬变电磁探测深度。图10c 与图10d 皆为航空磁法探测结果,图中对两条断层的位置与延伸皆有反应,其中图10 d 为航磁三维反演-1 200 m 磁化率等值平面图,图中断层F1 反应较好,表现为低磁化率条带,而对断层F2 的反应效果较差,实际航磁勘探中对深部异常探测效果较好,而对浅部异常反应效果较弱,推测断裂F2 埋深较小,断层F1 埋深较大,该推测与航空瞬变电磁探测结果推测一致。
图10 不同勘探方法探测效果对比Fig.10 Comparison of detection effects by different geophysical exploration methods
3 种探测方法对断裂F1 与断裂F2 皆有反应,3 种方法相互印证,发现研究区两条隐伏断裂,断裂F2 埋深较小,断层F1 埋深较大。由上述分析可以看出,航空大地电磁倾子总散度对于横向电性界面变化的探测行之有效且效果较好,对于深部断层F1 和浅部断层F2 皆有反应,即该方法具有兼顾浅部与深部的探测能力。6 个频率的总散度平面图中对于两条断层反应特征相似,单纯通过该方法不能很好的推断断层的埋深,即航空大地电磁倾子总散度参数的纵向分辨率较低,为到达更好的勘探效果,需结合其他物探方法对资料进行解释。
航空大地电磁法利用飞行器作为探测载体,具有探测效率高,勘探深度大,受地形影响小等优点,拥有良好的推广前景。本文通过数据模拟和实例分析得出以下结论:
1)航空大地电磁勘探以倾子作为解释参数,对横向电性突变界面反应灵敏,可用于岩性分界面、断裂等构造解释。倾子资料Tx只对x轴方向电阻率变化有响应,倾子Ty只对y轴方向电阻率变化有响应,并在电性交界面上形成异常中心。
2)倾子总散度可以将电性交界面的极值异常转换为异常体上方的极值异常,并且倾子总散度对于低阻异常和高阻异常边界圈定皆较为精确。同等地电条件下,倾子总散度高频响应相对于低频响应数据范围大,对比度高,探测效果更好。
3)实际应用表明航空大地电磁对于断层破碎带具有很好的探测效果,推断研究区中两条隐伏断裂,并通过航空瞬变电磁和航磁探测结果进行印证。对比航空瞬变电磁和航空磁法勘探,航空大地电磁倾子总散度参数用于地质解释行之有效,对浅部和深部探测效果皆较好。但倾子总散度参数的纵向分辨能力较差,为到达更好的勘探效果,需结合其他物探方法对资料进行综合解释。