1. 引言
据国务院安委会统计,至2015年底全国金属非金属地下矿山共有采空区12.8亿m³,分布于全国28个省(市、区)。采空区作为诱发重特大事故的重要危险源,其精准探测对人民生命财产安全有着重要意义,有着高度的社会价值[1]。然而金矿开采后形成的采空区往往规模小、埋藏深、形状不规则,地质探测识别难度大。地球物理勘探技术凭借普查范围广、作业效率高、经济性突出等特点,已成为采空区调查的主要技术手段。常规地球物理勘探技术主要包括高密度电阻率法、地质雷达法、瞬变电磁法及微动勘探技术等。
针对金矿采空区探测的问题,地质雷达法探测深度过浅难以达成工程需要,常规电磁法易受环境电磁干扰且对于空区含水与否存在地质多解性。而微动勘探技术恰好可以弥补上述缺点,微动是指地表天然存在的微弱震动信号,震源一方面来自于海浪、潮汐变化、气压变化、构造运动等自然现象,频率往往小于1 Hz;另一方面来自于交通运输、工业活动等人类活动,频率往往大于1 Hz [2]。微动勘探技术便是对天然微动信号进行数据处理,从多分量振动数据中提取瑞雷面波的频散特征曲线,通过反演重构地下介质横波速度剖面,分析不同介质的横波波速差异,从而查明或解决相关地质问题的一种物探技术[3]。近年来,微动勘探技术与地震干涉法的有效结合[4]、密集台阵三维成像带来的良好反演效果[5]、在诸多个领域(地热勘查、隐伏构造勘查、采空区勘查等)的广泛应用[6]-[8],微动勘探技术的应用潜力与发展前景日益凸显。
2. 微动勘探原理
微动探测技术的基本原理是基于平稳随机过程理论和弹性波动理论,通过分析微动信号中的瑞雷波频散特性提取频散曲线,进而反演地下横波速度结构,最后进行解译地质现象。
2.1. 提取面波频散曲线
1957年,Aki在研究地震面波时提出了空间自相关法(Spatial Autocorrelation Method,SPAC法) [9],基于圆形阵列观测数据,通过计算空间自相关函数并拟合零阶贝塞尔函数,提取面波的频散曲线;1969年Capon提出了频率–波数谱分析方法(Frequency Wavenumber Method, F-K),该方法通过对信号进行二维傅里叶变换,得到频率–波数谱,从而提取面波的频散特征曲线。2003年Okada在SPAC法的基础上提出了扩展空间自相关法(Extended Spatial Autocorrelation Method, ESPAC) [10],通过计算不同传感器对之间的空间自相关函数,拟合贝塞尔函数以提取频散曲线。ESPAC法扩展了SPAC法的适用范围,实现了对任意阵列布设的支持,显著提升了微动探测法在复杂实际场景中的实用性和可靠性。具体计算过程如下:
假设布设阵列中存在A、B两点,其空间自相关系数可表示为:
。
式中:
是A与B的ESPAC互相关谱,
是A和B间的距离,M是AB观测数据分段数,
和
分别是A和B的第m段数据的自相关谱,
是两点间的第m段数据的互相关谱。
把所有的观测点(N个)的互相关谱进行叠加,得到平均互相关谱
,
。
空间自相关系数与第一类贝塞尔函数的关系式为:
。
式中:
是台站间距为r、角频率为ω时的空间自相关系数,
表示观测点的自相关函数,
表示其他点和观测点的互相关函数,c是相速度,
代表第一类贝塞尔函数。
2.2. 横波速度结构反演
面波相速度是指同频面波传播的相位速度,其随频率变化产生的差异称为频散现象。在层状介质中,频散特性与地层纵波速度(Vp)、横波速度(Vs)、密度及层厚构成复杂非线性关系。Horike通过长周期微动反演发现,S波速度对相速度的敏感性显著高于P波和密度参数。焦健采用Knopoff算法定量研究表明:当横波速度产生10%波动时,面波相速度平均变化率可达或超过10%;相较而言,地层厚度变化仅引起微弱影响(相速度变化约1%),而Vp和密度参数的影响量级可忽略(不足0.5%)。该结论揭示横波速度是主导面波频散的核心控制因素,因此在工程反演计算中,仅需重点构建横波速度模型即可有效表征面波频散特性,显著降低多维参数反演的计算复杂度。为地下介质横波速度结构的精确反演提供了重要理论依据。
3. 应用实例
3.1. 研究区概况
本次研究区为胶东某金矿采空区,位于山东省烟台市蓬莱区南部,地处胶东半岛西北部,是胶东金矿集中区的重要组成部分。构造上隶属华北克拉通东南缘胶北隆起东翼,处于蓬莱东南部金成矿带核心区。研究区出露地层简单,仅见新生界第四系;构造发育,以断裂构造为主;岩浆岩广布,主要为中生代燕山早期郭家岭超单元;区内构造由韧性剪切带和脆性断裂组成,韧性剪切带主要发育于区域西南部虎路线一带的新太古代栖霞超单元回龙夼单元中,总体走向南东,表现为强烈的韧性变形及糜棱岩化。该地区采矿历史悠久,区内采金业发达,中小型金矿山密集,叠加电网设施、人文活动及工业建筑群构成了复合干扰场。
3.2. 测线及参数布设
微动探测进行测量前首先要确定系统的主要参数,具体包括:采样时长、采样间隔、采集模式、台阵半径,测点间距等。为达到探测目的,结合采空区地质条件,本次微动探测采用了线性排列的观测系统,使用国为公司的GN309智能地震仪,布置L1、L2、L3、L4测线4条,测线总长1200 m,台阵半径为50 m,测点间距为5 m,所有测线施工方向由西北至东南布设在采空区上方地表。为了保证数据质量,应尽可能延长采样时间,确保采样时长内包含的信息更加丰富。但工程上需要兼顾工作效率无法延长采样时间,因此为了尽可能保证数据质量的同时增加工作效率,我们在测量前进行了采样时长试验。
试验采用3组仪器对同一点进行同时测量,采样时长分别设置为20 min、40 min、60 min,得到的频散能量谱如图1所示,从中我们不难发现采样时长为20 min时频散能量谱较散乱、不连续;采样时长40 min与采样时长60 min两者频散能量谱相近,无明显差异。因此综合考虑我们最后采样时长定为40 min。
(a) 20 min频散能量谱 (b) 40 min频散能量谱 (c) 60 min频散能量谱
Figure 1. Dispersion energy spectra corresponding to different sampling durations
图1. 不同采样时长对应的频散能量谱
3.3. 数据处理
图2为本次微动勘探采集数据原始波形,为确保数据可靠需要进行数据实时预处理。该流程主要包括数据格式标准化和有效性验证两个方面:首先将原始观测数据转换为ESPAC分析专用格式,同步生成波形时序图用于监测信号完整性;其次通过噪声识别算法筛除强干扰数据段,并基于空间自相关系数矩阵评估台站阵列数据的可靠性,为后续施工方案优化提供决策依据。
Figure 2. Waveform plot of acquired raw data
图2. 采集数据原始波形图
数据处理阶段采用瑞雷波频散特性进行地层解译。首先从时域信号中提取多阶瑞雷波相速度频散曲线,四条测线提取的频散曲线如图3所示,我们可以发现,频散谱能量集中度高,曲线形态较为标准,数据质量较好;之后运用非线性反演算法将频散特征转换为横波速度模型;最后通过空间插值构建二维视S波速度剖面,并绘制相速度等值线图。其中,半波法经验公式作为经典反演手段,通过建立频散曲线形态与地层速度的统计关系,可有效识别地下构造界面。其数学模型可表示为:
,
式中
为地层类型相关系数,
为波长特征参数,
为相位差函数。最终地质解译需综合速度剖面特征与区域地质背景,实现地球物理参数向地质模型的科学转化,以获取相应的地层构造。
(a) L1线频散曲线
(b) L2线频散曲线
(c) L3线频散曲线
(d) L4线频散曲线
Figure 3. Dispersion curve plot
图3. 频散曲线图
3.4. 应用效果分析
不同地层介质中传播地震波的传播速度不同,总体而言介质密度越高,传播速度越高,在基岩中地震波速比含水空区、裂隙发育区、堆积物松散区中要高,利用这个特性可以在微动视横波速度剖图上对采空区位置分布进行推测。
下图是我们根据频散曲线求取的视横波速度剖面图(图4),从微动视横波速度剖面上看,视横波速度曲线整体较为规律,横向连续性较好,视横波曲线对浅部风化层界面和深部岩体刻画较为清晰,根据视横波曲线形态特征我们在圈定了9处低速异常圈闭,结合已有地质资料我们推测出重点异常位置5个、一般异常位置4个。
(a) L1线视横波速度剖面
(b) L2线视横波速度剖面
(c) L3线视横波速度剖面
(d) L4线视横波速度剖面
Figure 4. Apparent S-wave velocity section
图4. 视横波速度剖面图
3.5. 钻探验证
为进一步验证微动勘探反演结果的准确性,我们于L1线540点、L3线560点附近布置了两个验证钻孔ZK-1、钻孔ZK-2。其中钻孔ZK-1在42 m处见空区,钻孔ZK-2在120 m处见空区,与根据微动视横波速度剖面在此处推测的异常吻合度较高。
通过钻探的验证,表明了微动勘探推测结果的可靠性,上述9处异常区域为空区的可能性较高。
4. 结论
(1) 通过微动勘探技术推测出重点异常位置5个、一般异常位置4个,并采用钻探在L1线540点、L3线560点两处异常点进行验证,且于对应深度位置确认了空区;
(2) 微动勘探技术在金矿采空区探测上可以较好地圈定异常,方法可靠,效果较好;
(3) 经试验采样时长宜采用40 min,既可以保障工作效率,又在一定程度上尽可能多地采集有用信息;
(4) 微动勘探技术无需人工震源,绿色环保且工作效率高,更符合现今物探技术发展趋势。