基于全区视电阻率瞬变电磁一维正演中滤波系数的选取
2021-02-24郭嵩巍刘小畔
郭嵩巍, 郑 凯, 刘小畔, 王 成, 王 淼
(1.中国地质大学(北京)地球科学与资源学院, 北京 100083; 2.长江大学地球物理与石油资源学院, 武汉 430100;3.浙江省地矿勘探院, 杭州 310000; 4.成都理工大学地球物理学院, 成都 610059; 5.内蒙古国土资源勘查开发有限责任公司, 呼和浩特 010020)
瞬变电磁法(transient electromagnetic method,TEM)是一种较为成熟的地球物理电法勘探方法,被广泛应用于矿产、水文、工程等领域[1-2]。在工程应用领域,一维反演依然是目前瞬变电磁法资料解释的重要依据之一[3-4]。然而,瞬变电磁正演计算效率低下是制约反演应用的重要因素,为提高正演的计算效率,选取合适的滤波系数是关键。
瞬变电磁法装置多样,发射场源包括磁性源、电性源[5-8],接收装置普遍采用回线线框,采集的数据大多为感应电动势,可解析得到晚期视电阻率。由于晚期视电阻率在早期段,即测深曲线的首支存在明显的高阻假象,而全区视电阻率很好地解决了这个问题,对地质体电性表达更为真实,因此转换全区视电阻率逐渐成为处理解释的主流。全区视电阻率一般可通过数值计算得出,无论是磁性源还是电性源,很多学者都将感应电动势转换成垂直磁场,进而计算出全区视电阻率[9-12]。因此以全区视电阻率作为滤波系数选取研究的正演参数。
通过对麦克斯韦方程组的推导,可得出瞬变电磁一维正演响应的解析表达式。其中核函数的贝塞尔函数积分问题,通常采用汉克尔变换[13-15];频率域到时间域的转换方面,有G-S变换[16-17]、余弦变换[18]、逆Laplace变换[19]等,本文采用余弦变换。汉克尔变换和余弦变换数值算法均采用数值滤波算法,选择合适的滤波系数决定着瞬变电磁正演的计算效率、精度以及有效范围。分别选取正演效果较好的三组滤波系数,对比分析102Ω·m均匀半空间模型正演响应,得出全区视电阻的有效范围,并以该模型为参照模型,利用瞬变电磁响应的平移伸缩特性[20],提出参照模型法,以为工程应用的正反演计算中选取滤波系数,提供一个简单可行的估算方法。
1 选取滤波系数问题的由来
以最为典型的中心回线装置为例,一维正演的数值计算可简化为核函数的两层积分,内层积分是贝塞尔函数积分问题,采用汉克尔变换,用Ht{}表示;外层积分是电磁场响应从频率域到时间域的转换问题,本文采用余弦变换,用Ct{}表示;核函数用fc表示,瞬变响应的垂直磁场用hz表示,于是中心回线瞬变电磁[21]一维正演响应可简化表示为
(1)
(2)
式中:I0为发射电流;a为发射线框半径;t为正演自变量时间;u(1)为顶层波阻抗;λ为核函数汉克尔变换无穷积分的自变量;ω为余弦变换的自变量频率。
中心回线瞬变电磁一维正演响应进一步简化为
(3)
对于均匀半空间地电模型,可以推导出垂直磁场hz的解析表达式为
(4)
式(4)中:erf(u)为误差函数。
归一化hz得到
(5)
于是就将瞬变响应转换为基于垂直磁场的全区视电阻率,即
(6)
式(6)中:μ0为真空中磁导率。
由于函数Z(u)为隐函数,无法解析导出反函数,但可利用函数在0~1单调递增的特点得到结果通过二分查找的办法实现了该全区视电阻率的数值计算,并与其他方法做了对比分析,得出了该全区视电阻反演效果更好的结论[9]。
于是瞬变电磁的正演数值算法的核心问题,就转变为汉克尔变换和余弦变换的数值滤波计算问题。一般来说,滤波系数采样点越多,正演的计算精度越高,正演响应的有效范围就越大,但计算效率就越低。若滤波系数不足或不符合正演响应有效范围要求,将不能模拟准确完整的正演曲线,反演时会造成拟合失真,再做解释将是灾难性的。因此,有必要择优选取滤波系数。
2 数值滤波算法
数值滤波算法,通常是为了解决无穷积分而设计的,汉克尔变换解决的是含有贝塞尔函数的无穷积分,余弦变换是对傅里叶变换的化简,实现了从频率域到时间域转换的无穷积分。
滤波系数的得出是通过变量替换,将无穷积分式变成卷积的形式。在满足抽样定理的条件下,通过求解卷积方程得到指数间隔采样的滤波系数,即反算滤波器,方法有傅里叶变换法、维纳-霍普夫极小化法[22]等。
数值滤波算法通常形式为
(7)
(8)
式中:B和W为一组滤波系数变换对,数组长度为n,B一般采用指数间隔采样,λ为数值滤波算法中卷积形式的采样点,在汉克尔变换中对应贝塞尔函数的采样点,在余弦变换中对应频率ω的采样点。
推荐汉克尔变换为Guptasarma47点算法[13],即
B=10[a+(i-1)s],i=1,2,…,n
(9)
式(9)中:n=47,a=-5.08 25,s=1.166 383 038 62×10-1。
余弦变换为Kong21点算法[23],即
B=exp(S+iT),i=-L~L
(10)
式(10)中:L=11,S=1.40,T=0.35。
3 参照模型法
由于实测数据的发射装置不同、对应地电模型电阻率不同、采集数据响应时间的上下限不同,要选取合适的滤波系数,就必须将“不同”变成“相同”,因此提出参照模型法。该方法是一种简单的估算方法,可将发射装置不同、模型电阻率不同、响应时间上下限不同的实测数据,转换到对应参照模型的时间范围,再与参照模型不同滤波系数组合正演模拟的有效范围做比对,从而做到对滤波系数的选取。
这里需要说明的是,本文参照模型的有效范围是中心回线装置的计算结果,其他装置的有效范围可通过相应正演模拟计算得出。
现有地质地球物理资料显示,地壳内大多数地质体的电阻率范围在1~104Ω·m,取对数中间值102Ω·m作为正演参照模型的电阻率(ρr)。发射装置方面,煤田首采区勘查多采用边长3×102~103m的大回线框扫面;水文勘查多采用102~3×102m中心回线装置;工程勘查常采用多扎小线框,边长n~n×10 m;取对数中间值102m作为参照模型的发射线框边长,参照发射线框面积(Sr)为104m2。
根据瞬变电磁的平移特性
(11)
(12)
得到实测数据时间与其参照模型时间的关系:
(13)
对式(13)取对数,得
lgtr=lgts+lg(ρs/ρr)+lg(Sr/Ss)
(14)
式中:Sr为参照模型发射线框的面积,104m2;ρr为参照模型电阻率,102Ω·m;tr为参照模型时间;Ss为实测数据发射线框的面积,104m2;ρs为实测数据模型电阻率,102Ω·m,如假设为反演模型的电阻率上下限范围等;ts为实测数据的时间,s。
由于实测数据并非对应均匀半空间模型,可根据全区视电阻率的范围推断实测数据的模型电阻率范围,或者采用反演模型约束电阻率的上下限范围估算,总之需要对实测数据对应的模型电阻率估算一个上下限范围。令
(15)
c2=lg(Sr/Ss)
(16)
最终得到
(17)
4 不同滤波系数组合的有效范围
经过对多组汉克尔变换滤波系数、余弦变换滤波系数的对比试验,择优选出3组汉克尔变换系数和3组余弦变换系数,计算参照模型的正演响应,并转换为全区视电阻率,设定相对误差在10%以内为有效数据,得到有效范围(表1)。余弦变换固定81点滤波系数,汉克尔变换分别对Guptasama 47点、201点、401点[24]滤波系数做正演,对比结果见图1、表1(组合编号1、2、3);汉克尔变换固定Guptasarma47点滤波系数,采用3组不同的余弦变换分别采用21点、81点、601点滤波系数[24],对比结果见图2、表1(组合编号4、1、5)。
表1 参照模型有效时间范围统计表Table 1 Statistical table of effective time range of reference model
图1 3组汉克尔变换滤波系数对比Fig.1 Comparison of the Hankel transform filters coefficient
图2 3组余弦变换滤波系数对比Fig.2 Comparison of the 3 groups cosine transform filter coefficient
5 实测数据的参照模型时间范围
根据前面推导的参照模型法,分别列出了内蒙古自治区不同地区的不同工区、不同发射装置、不同地电模型的实测数据(全部采集自V8电法工作站),估算出对应参照模型的时间范围。由于变化范围大,用对数值lgt表示,见表2。
与(表1)比对可以发现,通过参照模型法的估算,本文提供的滤波系数组合的有效时间范围,全部满足实测数据的计算需求。
6 数据处理软件模块的开发
择优选取滤波系数,目的是为了更好地实现正演建模和对实测数据的反演计算,因此基于VB2015开发了瞬变电磁数据处理软件模块,实现了数据读取、显示、编辑,全区视电阻率转换、烟圈反演、一维反演等功能(图3)。作为正演、反演计算的设置参数,程序中可根据实际情况需要,选择滤波系数的任意组合,其中组合4计算效率最高,“参照模型法”估算其可以满足大部分实测数据的计算需求(表1、表2),因此程序将组合4设置为缺省值。
图3 数据处理模块软件界面Fig.3 Data processing software module interface
7 实测剖面算例
为证明“参照模型法”估算选取滤波系数的可行性,选择内蒙古自治区凉城县某地热勘查项目的一条剖面的实测数据作为测试数据,反演算法选择正则化拟二维反演(ARIA)[25-26],该项目施工条件良好,数据采集质量高,最终在北东向和北西向构造交汇处布置地热井(DR2),取得了不错的找水效果[27]。
首先,实测数据剖面位于地热钻孔(DR2)附近,地层电性分层明确,基底埋深确定,可确保反演结果的正确性。其次,该数据发射线框600 m×400 m,采集时间8.9×10-5~2.1×10-2s,按照“参照模型法”估算的时间范围(表2),对数值在-3.67~0.40,范围较大,有利于测试各个滤波系数组合的有效范围是否满足实际需求,若有效范围偏小,对曲线尾支的拟合影响较大,在反演断面图上表现为基底偏差较大,不同滤波系数组合呈现出不同的反演结果,因此该数据具备对滤波系数选取的测试条件。
表2 实测数据转换参照模型的时间范围Table 2 The transform of reference model time range from measured data
7.1 地质概况
7.1.1 地层
工作区地层主要有太古界桑干群、中生界侏罗系及新生界古近系、新近系、第四系。由老至新分述如下:
太古界桑干群(Ar1Sg),主要分布于蛮汉山区及平原区深部,为一套黑云矽线榴石斜长片麻岩、钾长片麻岩与浅色麻粒岩。
中生界侏罗系上侏罗统(J3),零星分布于蛮汉山区大西沟一带,为一套巨厚层状粉白色、灰紫色、灰白色凝灰岩、凝灰角砾岩、凝灰集块岩及安山岩组成,倾向北西,不整合覆于桑干群片麻岩之上。
新生界可划分为古近系与新近系:
古近系(E):见于钻孔,为一套深棕红色泥岩、泥质砂岩、砂质泥岩夹灰白色、深棕红色砂砾岩,不整合于太古界地层之上。
新近系(N):见于钻孔,上部为浅棕红色泥岩、泥质砂岩、砂质泥岩夹灰黄色砂砾岩,下部为土黄色、灰白色砂砾岩,与古近系地层呈假整合接触。
新生界第四系(Qal+pl):主要分布于蛮汉山山前倾斜平原,岩性为浅黄棕色、黄褐色粉土、粉质黏土与浅黄褐色砂砾卵石、含卵砾中粗砂、中粗砂等。粉质黏土或粉砂土中含有砂与砾石;砂砾卵石分选差,多棱角状,成分以片麻岩、凝灰岩为主。向平原南部岩性变细,为黄褐色粉土、粉质黏土与砂砾石。
7.1.2 构造
工作区区域地质构造(图4)属阴山东西向构造带的南缘。因受早期构造运动的影响,呈北东30°~40°延伸,系高角度压性断层,而北西向断裂据浅层地震资料所反映为一条北西325°的破碎带,宽度为100~200 m。
图4 工作区地质图Fig.4 Geological map of work area
7.1.3 岩浆岩
主要分布在蛮汉山区及平原区深部,为太古代早期侵入岩,岩性以苏长岩和似斑状花岗岩为主。
7.2 瞬变电磁反演效果对比
测试剖面800线,点距10 m,测点共51个,剖面总长500 m,位于山前倾斜平原。地热钻孔(DR2)位于剖面附近,钻孔显示地层分层明确、岩性差异明显,存在明显电性差异:浅部第四系(Q)砂岩黏土总体呈高阻,电阻率为30~80 Ω·m;中间层为古近系(E)、新近系(N)泥岩地层,为低阻层,电阻率为10~30 Ω·m,顶板埋深约150 m;基底花岗岩高阻,电阻率大于100 Ω·m,顶板埋深大于300 m。
反演结果(图5)显示,不同滤波系数组合的反演结果基本一致。电性分层与DR2钻孔基本吻合,说明决定反演效果主要因素是正演响应的有效范围,正演响应的精度(相对误差)对反演结果的影响相对较小。因此认为正演的有效范围和计算效率更具实际意义(表3),对测点较多、计算量较大的工程应用可影响反演的效率和效果。
图5 不同滤波系数组合方案正则化反演结果对比Fig.5 Comparison of ARIA inversion results of different filter coefficient selection schemes
表3 实测剖面反演效率统计表Table 3 Statistics of inversion efficiency of measured profile
8 结论
(1)瞬变电磁一维正演响应的计算效果(有效范围)与滤波系数的取值相关,计算效率与滤波系数的数量相关,成正相关关系。共推荐并测试了正演效果良好的三组汉克尔变换和余弦变换的滤波系数。
(2)通常情况下,计算效果好的滤波系数组合普遍积分区间较长,计算效率不高,可通过参照模型法估算出实测数据反演的计算要求,做到合理选取滤波系数组合。
(3)对实测剖面的反演测试显示,有效范围并非越大越好,正演响应(全区视电阻率)有效范围大牺牲的是计算效率,而计算精度的提升对于反演结果的影响较小,因此对于实际工程应用而言,选取滤波系数的原则是满足参照模型有效时间范围的情况下,尽可能选择滤波系数少,计算效率高的滤波系数组合。
综上所述,推荐滤波系数组合4,即汉克尔变换Gaptarsma47点、余弦变换Kong21点,为目前瞬变电磁实测数据反演计算的最佳滤波系数组合。