高密度电法反演数据预处理方法研究

尚新磊，于 悦，王天宇，王浩宇，王晓光

(1.吉林大学 仪器科学与电气工程学院，吉林 长春 130012；2.吉林大学 公共计算机教学与研究中心,吉林 长春 130012)

1 引 言

高密度电法因其较高的施工效率和勘探分辨率被广泛应用在地质调查中[1],起源于20世纪70年代末期的阵列电法探测思想[2]，高密度电法装置具有高效率、可以一次性获得更多地电信息的优点而被广泛应用于堤防隐患探测[3]以及水文[4]、工程[5]、环境等地质勘探中[6]。随着高密度电法的发展，高密度电法反演技术也得到迅速发展。2.5维(2.5D)反演是在2D反演的基础上发展起来的，对多个平行的剖面测量得到的视电阻率数据进行反演计算，得到多个剖面电阻率信息，高密度电法的2.5D反演能够更多地解释地质构造和地下环境。反演技术是正演技术的逆过程，是通过高密度电法仪器测量的电阻率数据获得地质信息的重要过程。自动反演法是电法反演的重要手段，随着计算机技术的发展，该方法也得到了广泛的应用，其过程为：通过给定初始模型，计算初始模型的视电阻率数据，与实际测得的视电阻率数据进行对比获得误差值，利用误差值获得模型修正参数对模型进行修正，利用修正后的模型再进行正演计算，不断重复修正过程，直到满足精度[7]，最后获得的模型参数就是实测曲线的解释结果

但是在反演中，由于仪器测量过程中会引入噪声，主要包括白噪声、工频噪声等，反演会出现不收敛的情况。在反演算法不收敛时，会重复改变初始模型但不能满足精度要求，这会使得反演时间增加，且反演结果不可靠。因此，对仪器测量得到的电阻率数据进行预处理，尽可能保证反演过程中的收敛性，对缩短自动反演时间，并使反演结果可靠是十分重要的。

对测得的视电阻率进行预处理，主要是对测得的视电阻率中的噪声进行滤除，现在对噪声数字滤波手段较多，包括中值滤波[8]、平滑滤波[9]、限幅平均滤波[10]等。

这些滤波手段都针对特定噪声效果明显，但是高密度电法工作环境特殊，往往伴随特殊噪声，例如尖峰噪声等，单一的滤波方法往往不能够应对高密度电法仪器测量的视电阻率中存在的误差。自适应滤波技术对各种“不确定性”、不知道来源的噪声都有着很好的滤除效果[11]，现在在各个领域得到广泛应用，并且有较好的表现。2016年，李菲等[12]应用自适应滤波法于电子天平抗震方面，通过实验验证，利用该方法研制的200 g/0.001 g应变片电子天平能够在木制试验台准确称重，并且不受榔头敲击实验台等冲击性震动的影响。2019年，吴作鹏等[13]基于归一化LMS滤波校正的核辐射法液位计，并且实验证明，该方法能够消除测量系统噪声干扰,提高液位测量精度。2020年，李杨等[14]将自适应滤波方法应用于多相电机故障信号提取方面，通过自适应滤波算法滤除含有噪声污染的故障信号中电机产生的电磁噪声，有效地提取出了故障信号。除此外，在图像信号处理[15]，超声信号[16]，语音端点检测[17]等方面等领域，自适应滤波方法都有应用。

而为尽可能防止数据反演发散的情况发生，考虑利用自适应滤波方法对视电阻率进行滤波处理，达到提高反演质量的目的。

2 自适应滤波及自动反演基本原理

为防止反演算法不收敛导致的反演结果不可靠和反演时间增加,先对视电阻率进行自适应滤波处理。再通过自动反演算法对预处理后的视电阻率数据进行反演。

2.1 用于视电阻率滤波的自适应滤波算法基本原理

针对视电阻率无采集时间先后差异的特点，避免数据抽头损失，采用逆向滤波弥补抽头的方式减少数据损失。

首先通过正向自适应滤波算法对视电阻率数据进行滤波处理，正向自适应滤波算法可以表示为：

e1(n)=E1(n)-RT(n)W1(n)

(n=K∶N)

(1)

W1(n+1)=W1(n)+2u1R(n)e1(n)

(n=K∶N)

(2)

(3)

其中，R1(n)为n点及之后测量的视电阻率，单位为Ω·m；W1(n)为n点及之后N阶自适应滤波器的权重系数；E1(n)表示期望信号；e1(n)表示误差值，单位为Ω·m；u1是滤波过程中的步长因子；K1为滤波选用的抽头长度，用来预测的数据长度。可以看出，获得的滤波电阻率数据长度为N-K，会损失抽头长度，也就是y1只具有除抽头外的滤波视电阻率数据，为弥补损失的视电阻率数据，再进行反向自适应滤波。

反向自适应滤波算法可以表示为：

e2(m)=E2(m)-RT(m)W2(m)

(m=N-K∶1)

(4)

W2(m+1)=W2(m)+2u2R(m)e2(m)

(m=N-K∶1)

(5)

(6)

其中，R2(m)表示m点及之前测量的视电阻率值，单位为Ω·m；W2(m)表示m点及之前测量的N阶自适应滤波器的权重系数；E2(m)表示当前长度值的期望信号；e2(m)表示误差值，单位为Ω·m；u2是步长因子；K2为抽头长度。反向抽头滤波会损失尾部视电阻率数据，可通过互补的方式获得完整的视电阻率滤波数据，取Y作为最后的数据。

Y=[y1(K+1∶N),y2(1∶K+1)

](7)

2.2 自动反演基本原理

采用自适应滤波算法对具有不确定噪声的视电阻率进行滤波后，采用自动反演方法对预处理后的视电阻率进行反演运算。自动反演的基本原理可归结为寻找一种地下模型，使地下模型的视电阻率与测得的真实视电阻率尽可能接近。

设置待寻找的模型为Mj(j=1∶MN)，j为模型参数，包括每层厚度Th(单位为m)，和每层的电阻率ρi(单位为Ω·m)。为寻找模型Mj，首先设置初始模型，初始模型主要包括厚度Th，单位为m，和每层的电阻率ρi,单位为Ω·m。

使寻找的模型的视电阻率与真视电阻率满足最小二乘法的精度要求:

(8)

为寻找模型Mj，将模型M在预测模型处展开，只保留二次项:

(9)

其中,ΔMj为修正量。通过式(9)来获得模型的修正量，从而获得接近真实模型的反演模型。在式(9)中要获得修正量，则需要计算视电阻率的偏导值。对于数值反演可通过差分方法获得偏导值:

(10)

则修正量可表示为：

(11)

获得模型修正参数后，对修正后的模型参数进行正演计算，与测量参数进行比较，可以得到误差值。通常情况下，一次修正量是不能满足精度要求的。为了保证反演过程是收敛的，每次修正都限制模型参数的修正量在一定的范围之内。因此一般需要多次迭代获得多次修正量，逐步对模型进行修正，直到满足参数误差精度要求为止，最后获得的反演参数即为模型参数。

3 仿真验证

为验证视电阻率进行自适应滤波对反演过程收敛的有效性，通过仿真分别对存在白噪声、工频噪声的多模型视电阻率数据进行数据预处理前和预处理后反演结果的对比实验。

实验设备：自制的2.5D电法正反演仿真平台，该仿真平台具有如下功能：①能够根据需要对不同形状地下异常体进行建模，并能够对模型进行正演，获得视电阻率曲线；②为观测视电阻率曲线中存在噪声情况对反演结果的影响，该仿真平台能够在视电阻率曲线中加入白噪声、工频噪声等，并能够进行视电阻率自适应滤波处理；③对正演生成的视电阻率进行反演计算。

3.1 实验步骤

1)设置地下电阻率异常体形状为球体具体参数如图1所示。

图1 地下异常体示意Fig.1 Schematic diagram of underground abnormal body

其中，球体异常区域与总区域同心，设置总区域边长d0为3 000 cm，球体半径r为1 000 cm。

2)正演模型的建立：对异常体模型进行正演计算，正演方式采用温纳四极高密度电法探测方式，其探测示意如图2所示。

其中共铺设20道电极，设置电极间距为150 cm，被探测区域为边长3 000 cm的正方形，为能够将被探测区域探测完全，电极分布距离须为探测深度的两倍以上，电极分布在6 000 cm×6 000 cm的区域内。

3)进行正演计算，获得20条视电阻率曲线R1∶20。

4)同时向视电阻率曲线中增加白噪声和工频噪声获得两组存在不同噪声的视电阻率曲线RWF1∶20。

图2 电法探测示意Fig.2 Schematic diagram of electrical detection

对存在白噪声和工频噪声的视电阻率和利用自适应滤波算法对存在噪声的视电阻率进行滤波后的视电阻率进行反演，并记录每次反演时间。

3.2 数据分析

根据上述实验步骤，其球体异常区域的正演模型如图3所示，球体半径为1 000 cm。因为设置20道测量线，得到20组剖面图，其中每组剖面图间距为150 cm,实现2.5 D反演。设置异常体电阻为35 Ω·m，区域内土壤电阻率为200 Ω·m。对所建模型进行正演计算，获得20条视电阻率曲线。在视电阻率曲线中同时加入一个固定幅值白噪声和工频噪声，其中白噪声对电阻率测量结果的影响幅值设置为0.1 Ω·m，工频噪声波动幅值为0.1 Ω·m。

图3 地下异常体正演模型Fig.3 Forward modeling of underground anomalous body

图4 20道视电阻率预处理曲线Fig.4 20 apparent resistivity preprocessing channel

图5 滤波前反演结果Fig.5 Inversion results before filtering

图6 滤波后反演结果Fig.6 Inversion results after filtering

预处理后的20条视电阻率曲线如图4所示。通过对20条预处理前存在噪声的视电阻率曲线和预处理后的视电阻率曲线进行反演计算，得到2.5 D反演结果，其中存在噪声的视电阻率反演结果如图5所示，对噪声进行自适应滤波预处理后反演结果如图6所示。

从图5可以看出，对不进行预处理、存在白噪声和工频噪声的球体异常区域的视电阻率曲线反演结果异常区域中心电阻率约为65 Ω·m，边缘电阻率约为150 Ω·m，异常区域半径约为1 400 cm，电阻率测量相对误差约为10 %，异常区域半径相对误差约为40 %，反演时间为27 min。由图6可知，采用自适应滤波算法对存在噪声的视电阻率曲线进行滤波处理后，反演结果中异常体中心的视电阻率约为35 Ω·m，异常体边缘视电阻率约为120 Ω·m，异常体半径约为1 100 cm，反演时间为19 min。通过比较，对存在白噪声和工频噪声的进行自适应滤波处理能够明显改善反演结果，并且缩短反演时间约30 %。

4 总 结

在进行高密度电法探测地下异常体时，测量获得的视电阻率曲线中往往存在各种不确定噪声，带来反演过程不收敛，反演结果不准确，反演时间增加等问题。针对存在的各种“不确定性”噪声，提出采用自适应滤波方法对存在噪声的视电阻率进行预处理，再进行反演计算，通过仿真验证该方法的有效性。首先构建地下球体半径为1 000 cm，电阻率为35 Ω·m的异常区域，在假设铺设20道测量线的情况下，进行正演计算，获得20组视电阻率曲线，在视电阻率曲线中加入白噪声和工频噪声，再通过自适应滤波算法对存在噪声的视电阻率进行滤波预处理。对存在噪声的视电阻率和滤波后的视电阻率进行2.5 D反演计算，不进行预处理，存在白噪声和工频噪声的球体异常区域的视电阻率曲线反演结果，异常区域中心电阻率约为50 Ω·m，边缘电阻率约为150 Ω·m，异常区域半径约为1 300 cm，电阻率测量相对误差约为10 %，异常区域半径相对误差约为40 %，反演时间为27 min。采用自适应滤波算法对存在噪声的视电阻率曲线进行滤波处理后，反演结果中异常体中心的视电阻率约为35 Ω·m，异常体边缘电阻率约为120 Ω·m，异常体半径约为1 100 cm，反演时间为19 min。通过比较，对存在白噪声和工频噪声的进行自适应滤波处理能够明显改善反演结果，并且缩短反演时间约30 %。