基于岩石物性和地下水矿化度约束的三维电性模型地下水分布特征研究

辛会翠,徐志敏,2,肖 晓

(1.承德石油高等专科学校,河北 承德 067000;2.河北省仪器仪表工程技术研究中心,河北 承德 067000;3.中南大学 地球科学与信息物理学院,长沙 410083)

引言

柴达木盆地地下水分布不均是困扰该地区工业发展的主要障碍之一。因此加强柴达木盆地地下水资源调查,查明地下水分布规律,增加水资源储备,对盆地水资源规划与利用具有重要意义。适用于地下水资源探测的地球物理方法较多,如直流电测深(DC)法通过视电阻率参数可确定含水层结构;激电测深法的激电参数可用于了解地层富水性［1］,但受限于对接地条件的较高要求。本文涉及的柴达木盆地研究区,浅部分布了厚度较大的第四系松散沉积层,且在大、小柴旦湖及周边一定范围内分布广泛的盐碱地,形成了浅部电阻率相对较低的电性低阻层,常规直流电法难以克服浅部低阻地层的影响,探测深度和探测精度均受到较大影响。音频大地电磁法(AMT)是基于电磁波的趋肤效应,利用天然电磁场作为场源的一种频率域电磁测深法［2-4］。该方法具有施工简便,探测效率较高、探测精度高以及探测深度大等特点,在水文及环境地质调查等方面得到了广泛的应用。国内外学者利用音频大地电磁法进行了大量的生产科研实践,其在划分深部地质构造、寻找低阻异常体空间分布特征中有着良好的探测效果。例如,在地热水资源探测中AMT法可以有效识别热储盖层,揭示地热空间分布特征及地热分布区域地质构造特征,为查明地热成因,圈定地热靶区提供重要科学依据［5-8］;在第四系含水层水文地质地球物理探测中有着良好的探测效果［9-12］。大量实践证明,AMT法可以有效探测地下水资源空间分布特征。然而众多探测案例中,缺乏对其他物性参数的水文地质解释的讨论,采用多种物性参数进行水文地质综合解释显然可以减少物探结果的多解性,另外水文地质解释需要建立正确的地质-地球物理特征关系,物性资料的研究同样重要,因此加强物性特征以及多种物性参数的水文地质综合解释研究,对于提高水文地质解释精度具有重要意义。

岩石的物性特征是地球物理资料解释的基础,本文首先开展了研究区岩石物性特征分析,采用地层岩石物性测定方法,获得不同岩性区域岩石电性特征。采用AMT电磁数据一维反演方法,获得水文钻孔处岩石电性特征与AMT电磁数据之间的关系;然后,采用三维非线性共轭梯度反演算法(ModEM)进行AMT数据三维反演工作［13］,获得了可靠的柴达木盆地研究区地下3 000 m范围内的三维电导率模型。最后,结合三维电性模型分析结果,对第四系地层、含水层分布特征进行了进一步分析,并采用电阻率参数及地下水矿化度参数对区内不同深度地下水分布特征进行了进一步讨论,结合研究区地下水矿化度对比分析,验证了本文地质解释的准确性。本文研究成果为研究区水资源开发及利用提供了重要科学依据,具有重要的实际应用价值。

1 研究区岩石物性特征分析

1.1 岩石电阻率影响因素

1.1.1 岩石岩性

研究区第四系地层岩性以冲洪积卵石、砾石、粗砂、中砂和粉土等为主,其中卵石、砾石、粗砂、中砂电阻率值较大,粉土电阻率值较小。第四系地层以下为新近系地层,该层上部岩性以泥质粉砂岩及钙质泥岩为主,砾岩成分较少,电性以低阻特征为主;下部岩性以砂泥岩为主,受断裂构造影响,裂隙发育,充填泥质及黏土等良导矿物,电阻率值也较低。第三系地层以下为侏罗系地层,岩性主要为砂砾岩、炭质泥岩,受燕山期构造运动影响,地层结构破坏,沉积厚度及岩性变化较大,高低阻同时存在,电阻率值相对较高。基底岩性为震旦系花岗片麻岩,电性为高阻特征。

1.1.2 孔隙结构

研究区主要为孔隙结构的沉积岩,其中第四系砂砾卵石层粒度较大、分选性差,泥质含量不高,因而电阻率相对较高;而随着沉积深度的增加,岩石电阻率随粒度的减小、分选性变好、泥质含量增高、胶结程度变差和孔隙中水含量的增大,电阻率逐渐降低,因此第三系地层岩石电阻率值较第四系地层小而显示相对低阻特征。侏罗系地层粒度进一步减小,但受构造控制,局部岩体较破碎,电阻率值并未显著升高,至基底震旦系花岗片麻岩,岩石孔隙致密,电阻率值较高。

1.1.3 岩石良导矿物成分

研究区岩石良导矿物主要有黏土及泥岩,其中黏土矿物因其表面能够吸附阳离子产生扩散层产生附加导电作用;泥质矿物在外电场作用下,泥质颗粒表面吸附的离子沿表面移动形成附加导电作用,同时泥质充填在岩石颗粒之间或者附着在岩石颗粒表面,会对地层岩石的孔隙进行改造使得地层孔隙直径变小,孔隙结构变得复杂,地层束缚水含量增大,进而导致岩石电阻率降低［14］。研究区第四系砂砾卵石层在泥质含量较低时,地层赋水性较差,电阻率较高;在泥质含量较高时,由于泥质的存在,地层赋水性较好,电阻率较低。具体来讲,粉土及黏土电阻率最低;砂岩随着泥质含量降低,电阻率变大;砾岩由于泥质含量相对最低,故具有比砂岩高的电阻率。第三系地层泥质粉砂岩、砂泥岩随着地层泥质及黏土矿物的增多而显示低阻特征;侏罗系地层砂砾岩泥质含量较少而显示相对高阻特征;基底震旦系花岗片麻岩泥质含量最少,电阻率最大。

1.2 研究区岩石电阻率测量

1.2.1 均匀大地电阻率测量

基于前文分析,研究区岩石电阻率测定采用均匀大地电阻率测量法,针对不同地质单元布设了数量不等的电测深测点(图1)。本文采用对称四极电测深法［4］,在获得各测点原始测深数据资料后,对电测深数据进行了一维反演,最终获得了区域岩石电阻率特征,成果见表1。

图1 研究区AMT实测点位及物性测试点位分布图Fig.1 Distribution map of AMT sites and physical property test sites in the study area

1.2.2 水文钻孔岩性电性分析

本文针对研究区9个水文钻孔(图2,W1～W9)采用Occam(一维)反演和Bostick反演方法,对位于水文孔处AMT数据进行一维反演,2种反演方法相互约束,共同揭示了水文钻孔处地层岩石电性特征。鉴于篇幅有限,本文以水文钻孔ZK1和ZK2处的AMT测点W7-03和W7-13的一维反演解译成果为例展示钻孔岩性电性分析。

图2 AMT测点视电阻率—相位曲线、AMT测点一维反演拟合曲线及水文钻孔岩性柱状图Fig.2 AMT apparent resistivity phase curve,one-dimensional inversion fitting curve ofAMT measuring point and lithologic histogram of hydrological borehole

研究区地层岩石电阻率统计成果如表2所示。如前文所述,研究区含水层分布在第四系卵、砾石松散堆积层中,特别是第四系中、上更新统的卵、砾石层构成了本区地下水主要含水层。由表2可知,卵、砾石电阻率相对较大,特别当地层含有漂石时,电阻率值较大,因此可以推断的是,本区含水层电性特征显示为相对高阻特征,而第四系以下的第三系地层电阻率值相对较小,为相对低阻特征。

表2 研究区地层岩石电阻率统计表Table 2 Statistics of formation and rock resistivity in study area

通过研究区岩石物性特征分析本文取得如下成果:

第一,通过区域地质背景、地层岩性及水文地质特征分析,初步揭示了研究区水文地质特征,包括第四系地层、含水层分布特征;

第二,通过区域岩石电阻率影响因素分析,获得了本区地层岩石电性特征定性认识;

第三,通过区域岩石电阻率测量揭示了本区地层岩石电阻率特征;

第四,通过水文钻孔处AMT测点一维反演,建立了水文钻孔处地层岩石电性特征与AMT数据资料的相互关系;

综合以上研究成果,本文初步建立研究区水文地质先验模型,即研究区电性结构纵向为浅部第四系地层(含水层)为相对高阻特征,第三系地层为相对低阻特征;横向上冲洪积扇后缘地带电阻率值较大,为相对高阻特征;冲洪积平原区电阻率值逐渐减小,至冲洪积扇前缘地带,电阻率值最小。

研究成果为后续AMT三维电性模型分析及水文地质解释提供了可靠的先验信息和电性依据。特别说明的是,由于岩石电阻率影响因素较多,比如相同地下水矿化度情况下,泥质含量越高,电阻率越低,泥质含量相同时,矿化度越高,电阻率越低,还有孔隙度及地层岩性等,因此以上总结的岩石电阻率仅适用于本区物探资料及水文地质解释。

2 AMT数据采集与处理

2.1 AMT数据采集

本次AMT数据采集使用的仪器为加拿大凤凰公司生产的MTU-5A卫星同步大地电磁观测系统,配备MTC-30磁探头和不极化电极,共投入4台仪器,开工前后均进行了仪器的标定及一致性等工作,并在研究区外优选、布设了远参考站。测点采集坐标系为正南北向坐标系,采用张量采集模式,每个测点采集4道水平电磁场分量(Ex、Ey、Hx、Hy)数据,单点采集时间保持在40 min以上,工频滤波频率设为50 Hz,采集频率范围为0.35～10 400 Hz,实测数据剖面分布如图1所示。

2.2 AMT数据处理

研究区地处柴达木盆地北缘,深居内陆、地处高原,平均海拔2 800 m,人烟稀少,人为干扰较少。野外仪器采集到的数据是随时间变化的电场和磁场的时间序列数据。在大地电磁数据处理中,首先将时间序列数据转换为阻抗数据,本文音频大地电磁时间序列数据处理采用凤凰公司的SSMT2000处理软件,获得相应的阻抗张量、视电阻率和相位等大地电磁响应数据。

2.2.1 Robust阻抗估计

Robust估计算法就是根据数据离差的大小对数据进行加权,尤其注重没有受到干扰的数据,通过降低跳变点的权重系数来使得估算的阻抗效果更好［15］。

如图3所示,测点W1-40和W6-1经过Robust处理后获得的视电阻率-相位曲线光滑连续,数据质量较好,分析上述2个测点位于研究区电磁场环境“平静”区,电磁干扰微弱,可以忽略不计,经过Robust处理获得阻抗数据可以真实反映地下电性结构信息,可以直接参与反演计算。

2.2.2 大地电磁远参考法数据处理

当AMT数据受随机噪声干扰,全频段数据离差均较大时,Robust阻抗估计无法通过离差加权来获得数据质量较好的数据,大地电磁远参考法(Remote Reference MT,简称RRMT)假设两观测点相距较远时,两观测点间电磁场分量中的噪声满足不相关性,同时大地电磁场信号在一定范围内变化缓慢,将远参考点处磁场信号作为测点处的磁场分量来估算张量阻抗,进而提高阻抗张量估算精度［16-18］。

如图4-A和图5-A所示,测点W4-19和W4-24经过Robust阻抗估计,获得的视电阻率及相位曲线“飞点”较多,表明数据受到随机干扰噪声影响;全频段数据离差较大,Robust阻抗估计无法获得离差小的无偏阻抗估计。如图4-B和图5-B所示,测点W4-19和W4-24经过大地电磁远参考法处理后,视电阻率和相位曲线均获得的较大改善,飞点显著较少,曲线变得较为光滑连续。

图4 测点W4-19远参考处理前(A)后(B)视电阻率-相位曲线Fig.4 Apparent resistivity phase curve of measuring point W4-19 before and after remote reference processing

图5 测点W4-24远参考处理前(A)后(B)视电阻率-相位曲线Fig.5 Apparent resistivity phase curve of measuring point W4-24 before and after remote reference processing

2.2.3 AMT数据维性分析

获得大地电磁阻抗张量元素后,需对AMT数据进行维性分析,以确定研究区地下电性结构维性特征。判别电性结构维性的方法较多,主要有:阻抗极化椭圆分析、二维偏离度Skew维性指示因子和相位张量倾角因子。相位张量方法［19］与上述几种张量分解或维性判断方法不同,相位张量不需要事先对地下维性做出假设,因此存在局部三维异常体或者地下介质为三维构造时,相位张量方法仍然适用。图6展示了研究区226个测点在典型频率7 200 Hz、460 Hz、27.5 Hz和1.72 Hz的相位张量椭圆图。当相位张量呈圆形或近圆形、椭圆时表示地下为一维结构,当倾角因子(β)的的取值范围为-3°～3°时地下为二维结构,否则为三维结构［19］。由图6可见,从高频到低频,研究区AMT数据呈现由一维/二维到二维/三维的渐变规律,特别在1.72 Hz,AMT数据显示了非常明显的三维特征,从而说明了本文对选取的数据作三维反演的必要性。结合区域地质资料分析,研究区深部在燕山期地质时期形成了一系列断裂构造,三维性较强。

图6 音频大地电磁数据相位张量椭圆极化分布图Fig.6 Phase tensor elliptic polarization distribution of AMT data

3 AMT三维反演

3.1 三维反演数据准备及反演参数设置

前文经过一系列数据处理获得了质量较好的原始数据资料,以及AMT数据维性分析表明研究区深部三维结构发育,因此本文采用三维反演方法进行研究区AMT数据反演工作［13］,并依据数据质量优选了161个测点进行三维反演;大地电磁三维反演依赖于较大的计算资源,因此本文对参与反演的测点数据频点以对数等间隔计算抽取,从1～10 400 Hz频段选取了总计15个频点数据进行非对角线模式反演。三维反演迭代进行了89次,最终RMS(拟合差)为1.36,反演最终拟合差较小表明该反演结果具有较高的可靠性。

3.2 三维反演不同深度反演结果分析

下面通过不同深度三维反演结果分析,揭示研究区水文地质电性特征。

如图7所示,研究区0～600 m深度电性分区明显,在绿梁山以东,大头羊洪积扇以南大范围区域显示了相对高阻(150～250 Ω·m)特征。由前文区域岩石物性特征分析结果可知,本区含水层主要为第四系的卵、砾石层,电性特征在本区为相对高阻特征。同时,由图6所示,表明该区域含水层厚度较大,水量丰富,由于卵、砾石层孔隙较大,地下水径流较快,因此推测地下水类型为淡水;地下水向洪积扇前缘径流过程中,至前缘地带,由于含水层颗粒由粗变细,透水性减弱,地下水径流缓慢,潜水水位变浅,局部区域由于蒸发量极大造成地表盐碱化,地下水水质变差,因此推测为咸水。

图7 三维反演模型在0 m和600 m深度处的电阻率平面截面图Fig.7 Resistivity plane section of 3D inversion model at depth of 0 m and 600 m(图中蓝色虚线为推测咸淡水分界线)

研究区1.2～1.8 km深度(图8)范围电阻率值范围主要集中在10～200 Ω·m。该层为第三系下干柴沟组泥质粉砂岩及钙质泥岩,砾岩成分较少,因而导致地层电性显著降低。其中,在W5线绿梁山北侧一带显示了高阻区域,一方面可能由于该区域地层沉积结构较为致密,地层含水率较低;另一方面,由于靠近绿梁山一侧,地层抬升,第三系地层遭到剥蚀,而下伏地层电阻率值较高引起。

图8 三维反演模型在1 200 m和1 800 m深度处的电阻率平面截面图Fig.8 Resistivity plane section of 3D inversion model at depth of 1 200 m and 1 800 m

研究区2.4～3.0 km深度(图9)范围电阻率值范围主要集中在150～350 Ω·m。该层为第三系路乐河组砂泥岩和下伏基岩,大头羊洪积扇区域仍然存在一定低阻区域,表明区内存在较大断裂构造,使得第三系地层在大头羊洪积扇区域沉积厚度较大。区内基底东高西低,大致倾向北西,表明库尔雷克山一侧构造较为活跃,地层隆起幅度较绿梁山大,同时伴随绿梁山的隆起,地下水径流方向走向由东向西,在绿梁山向北西大柴旦湖转向。

图9 三维反演模型在2 400 m和3 000 m深度处的电阻率平面截面图Fig.9 Resistivity plane section of 3D inversion model at depth of 2 400 m and 3 000 m

4 研究区地下水矿化度特征分析

前文依据三维电性模型进行了区内咸淡水分区并分析了区内地下水分布特征,本节将通过地下水矿化度分析进一步研究区内咸淡水分区,以期提高研究区咸淡水分区解译精度。

4.1 研究区电阻率分布影响因素分析

4.1.1 岩石空隙结构特征

一般而言,岩石的电阻率由岩石内部胶结物和矿物颗粒的电阻率、形状及相对含量决定。由于本区含水层为孔隙构造岩石,因此本文仅讨论沉积孔隙构造岩石电阻率。对于大部分松散沉积岩类,如研究区第四系砂砾卵石层及砂砾石层孔隙结构较大,地下水连通性较好,其地层胶结物实际是地下水,而他们的矿物颗粒电阻率2远远大于水的电阻率［20］,岩石电阻率受岩石体积含水量或湿度影响,含水量或湿度微小变化可引起ρ岩石电阻率的极大变化。因此,G.E.阿尔奇于1942年通过试验发现了著名的阿尔奇公式［21］,并定义了一个参数,叫孔隙参数(或相对电阻率)用F表示:

(1)

式(1)中ρ0是岩石孔隙中100%含水时的地层电阻率,ρw为含水地层电阻率。F与φ的关系,由经验公式给出:

(2)

式中:a为比例系数(与岩性有关);m为胶结系数,对于砂岩类一般a取值1,m取值1.3;φ为岩石孔隙度,此公式对于孔隙度3%～5%至30%～40%的介质广泛适用。

4.1.2 地下水矿化度

一般而言,水的矿化度变化范围是0.01～1 g/L。矿化度的变化可以导致水的电阻率变化,特别在有盐溶液的情况下,电阻率变化范围可达几个数量级。地下水的电阻率与温度的关系在0 ℃～100 ℃范围内为线性关系。地下水的电阻率随着水矿化度的增长呈线性下降,设水的矿化度为C水,其电阻率满足经验公式:

log10ρ0=a+blog10C水

(3)

4.1.3 地下水温度

地下水电阻率由于温度升高离子活性增加,导致电阻率降低。由表4可知, 温度18 ℃时,C=1 g/L 时,ρ=5.6 Ω·m , 那么, 任意温度(0 ℃～200 ℃)a值为:

(4)

故可得由矿化度和温度同时作用时,水溶液的电阻率表达式:

(5)

式中t为温度(℃),α=0.025。

4.2 地下水矿化度分析

4.2.1 孔隙参数的确定

实际观测并反演后的电阻率, 代表着一定范围内地层的电性特征, 由其值进行地下水矿化度的求取可以反映地下水矿化度的空间分布特征。由于收集资料有限,研究区内水文观测井深一般100～300 m,因此本文研究研究区地下水矿化度空间分布特征深度范围在0～300 m。通过砂岩类孔隙度与渗透率关系公式［22］,计算求得研究区砂岩孔隙度为0.41～0.48,本文取孔隙度参数为0.45。通过收集研究区水文井观测资料,获得了区内砂岩渗透系数统计数据(表3)。

续表3

4.2.2 温度的校正

由前文可知,对电阻率的解释要注意温度的因素, 一般深度大于300 m以上的钻孔需要温度校正地下水电阻率值。本区含水层深度主要在100～300 m,据研究区水文观测井地下水温度资料(孔深120～300 m)显示,本区地下水温度变化较小(见表4),在4.5 ℃～9.0 ℃范围之间,由温度对地下水电阻率的影响公式计算可得地下水电阻率变化最大值为1.23 Ω·m,因此可以忽略温度的影响。

表4 研究区水文测井矿化度及水温统计表［23］Table 4 Statistics of mineralization and water temperature of the hydrological logging in the study area

研究区地下水矿化度计算,首先将孔隙度参数设为0.45可得研究区孔隙参数F为2.82;然后将三维反演获得的地层电阻率代入公式(2)计算可得研究区地下水电阻率(ρ0);其次将ρ0代入公式(3)计算获得研究区地下水矿化度。研究区含水层主要分布在100～300 m,因此本文主要研究0～300 m深度范围内地下水矿化度分布特征。

本文将计算获得的地下水矿化度与研究区水文观测井资料矿化度进行误差分析,由表5可知,仅J7观测孔矿化度误差达到10%,其他孔位误差较小,表明本文计算地下水矿化度较为可靠。

表5 三维反演电阻率计算地下随矿化度与研究区水文资料矿化度相对误差统计表Table 5 Statistical table of relative error of salinity of hydrologic data calculated by 3D inverse resistivity in the study area

如前文所述,研究区主要含水层深度为100～300 m,因此图6展示了研究区深度200 m和300 m地下水矿化度平面分布。

在图10-A中研究区矿化度低值(小于0.5 g/L)区主要集中在绿梁山东北及东南区域,表明该区域地下水运移速度较快,推测该区域为研究区地下水主要径流通道;同时研究区西北部地下水矿化度小于东南侧,表明区内主要淡水分布绿梁山北侧,塔塔凌河洪积扇冲积平原区。特别说明的是图中左侧测点未分布区域出现的高矿化度区域受周围测点共同约束形成,但由于未分布测点,约束力较弱,需要结合区域水文地质资料具体分析,鉴于该区域位于本文主要目标区以外,因此未做进一步解析。在图10-B中,位于绿梁山东侧至塔塔凌河出山口区域,地下水矿化度有所升高,表明随着地下水分布深度加大,地层泥质含量逐渐增多,地下水径流速度逐渐趋于静止,造成地下水矿化度逐渐增大。但在塔塔凌河现代河道位置,矿化度变化较小,分析由于在主河道位置地下水径流较快而显示低矿化度特征,在河道两侧不断沉积盐分而显示高矿化度特征。

图10 研究区地下水矿化度不同深度平面分布图Fig.10 Plane distribution of mineralized depth of groundwater in study area

结合前文三维电性模型不同深度电性特征分析结果总结如下:研究区在现代河床及冲洪积扇中缘矿化度多小于0.5 g/L,电阻率较高;至冲洪积扇前缘矿化度有逐渐升高的趋势,电阻率逐渐降低;塔塔凌河冲洪积扇右翼受绿梁山阻挡,形成阻水作用,地下水径流变缓在山前平原一带补给条件较差,矿化度大于0.5 g/L,电阻率较低。绿梁山南侧泄出带附近,受地层岩性和地下水径流相对缓慢的影响,矿化度为大于0.5 g/L,电阻率较低。

5 结 论

针对单一电阻率参数解译存在多解性问题,本文依据区域岩石物性特征分析、三维电性模型不同深度的平面电阻率分布特征和地下水矿化度分布特征对研究区地层地下水分布特征进行了综合解译,完成了研究区地下水空间分布特征精细化研究,取得了以下成果:

a.研究区有利的富水区域为绿梁山以东,大头羊洪积扇以南区域,同时受绿梁山阻隔,研究区地下水径流方向为向由东向西,在绿梁山向北西大柴旦湖转向的特征。

b.研究区主要含水层为第四系卵、砾石层,该层由于孔隙度大,造成电阻率特征为高阻特征,其中位于现代河床及冲洪积扇中缘区域为淡水区域的主要富集区,绿梁山山前平原一带以及绿梁山南侧为咸水区域。

c.本文研究成果为研究区地下水资源开发与利用提供了科学依据,需要说明的是由于本文区域岩石物性分析部分缺乏测井数据支撑,采用本文技术方案显然无法达到测井数据的精度,因此今后在研究水文地质问题时,物性特征分析应加强水文测井工作。