高密度電法反演數(shù)據(jù)預(yù)處理方法研究

尚新磊，于 悅，王天宇，王浩宇，王曉光

(1.吉林大學(xué) 儀器科學(xué)與電氣工程學(xué)院，吉林 長春 130012；2.吉林大學(xué) 公共計算機教學(xué)與研究中心,吉林 長春 130012)

1 引 言

高密度電法因其較高的施工效率和勘探分辨率被廣泛應(yīng)用在地質(zhì)調(diào)查中[1],起源于20世紀(jì)70年代末期的陣列電法探測思想[2]，高密度電法裝置具有高效率、可以一次性獲得更多地電信息的優(yōu)點而被廣泛應(yīng)用于堤防隱患探測[3]以及水文[4]、工程[5]、環(huán)境等地質(zhì)勘探中[6]。隨著高密度電法的發(fā)展，高密度電法反演技術(shù)也得到迅速發(fā)展。2.5維(2.5D)反演是在2D反演的基礎(chǔ)上發(fā)展起來的，對多個平行的剖面測量得到的視電阻率數(shù)據(jù)進行反演計算，得到多個剖面電阻率信息，高密度電法的2.5D反演能夠更多地解釋地質(zhì)構(gòu)造和地下環(huán)境。反演技術(shù)是正演技術(shù)的逆過程，是通過高密度電法儀器測量的電阻率數(shù)據(jù)獲得地質(zhì)信息的重要過程。自動反演法是電法反演的重要手段，隨著計算機技術(shù)的發(fā)展，該方法也得到了廣泛的應(yīng)用，其過程為：通過給定初始模型，計算初始模型的視電阻率數(shù)據(jù)，與實際測得的視電阻率數(shù)據(jù)進行對比獲得誤差值，利用誤差值獲得模型修正參數(shù)對模型進行修正，利用修正后的模型再進行正演計算，不斷重復(fù)修正過程，直到滿足精度[7]，最后獲得的模型參數(shù)就是實測曲線的解釋結(jié)果

但是在反演中，由于儀器測量過程中會引入噪聲，主要包括白噪聲、工頻噪聲等，反演會出現(xiàn)不收斂的情況。在反演算法不收斂時，會重復(fù)改變初始模型但不能滿足精度要求，這會使得反演時間增加，且反演結(jié)果不可靠。因此，對儀器測量得到的電阻率數(shù)據(jù)進行預(yù)處理，盡可能保證反演過程中的收斂性，對縮短自動反演時間，并使反演結(jié)果可靠是十分重要的。

對測得的視電阻率進行預(yù)處理，主要是對測得的視電阻率中的噪聲進行濾除，現(xiàn)在對噪聲數(shù)字濾波手段較多，包括中值濾波[8]、平滑濾波[9]、限幅平均濾波[10]等。

這些濾波手段都針對特定噪聲效果明顯，但是高密度電法工作環(huán)境特殊，往往伴隨特殊噪聲，例如尖峰噪聲等，單一的濾波方法往往不能夠應(yīng)對高密度電法儀器測量的視電阻率中存在的誤差。自適應(yīng)濾波技術(shù)對各種“不確定性”、不知道來源的噪聲都有著很好的濾除效果[11]，現(xiàn)在在各個領(lǐng)域得到廣泛應(yīng)用，并且有較好的表現(xiàn)。2016年，李菲等[12]應(yīng)用自適應(yīng)濾波法于電子天平抗震方面，通過實驗驗證，利用該方法研制的200 g/0.001 g應(yīng)變片電子天平能夠在木制試驗臺準(zhǔn)確稱重，并且不受榔頭敲擊實驗臺等沖擊性震動的影響。2019年，吳作鵬等[13]基于歸一化LMS濾波校正的核輻射法液位計，并且實驗證明，該方法能夠消除測量系統(tǒng)噪聲干擾,提高液位測量精度。2020年，李楊等[14]將自適應(yīng)濾波方法應(yīng)用于多相電機故障信號提取方面，通過自適應(yīng)濾波算法濾除含有噪聲污染的故障信號中電機產(chǎn)生的電磁噪聲，有效地提取出了故障信號。除此外，在圖像信號處理[15]，超聲信號[16]，語音端點檢測[17]等方面等領(lǐng)域，自適應(yīng)濾波方法都有應(yīng)用。

而為盡可能防止數(shù)據(jù)反演發(fā)散的情況發(fā)生，考慮利用自適應(yīng)濾波方法對視電阻率進行濾波處理，達到提高反演質(zhì)量的目的。

2 自適應(yīng)濾波及自動反演基本原理

為防止反演算法不收斂導(dǎo)致的反演結(jié)果不可靠和反演時間增加,先對視電阻率進行自適應(yīng)濾波處理。再通過自動反演算法對預(yù)處理后的視電阻率數(shù)據(jù)進行反演。

2.1 用于視電阻率濾波的自適應(yīng)濾波算法基本原理

針對視電阻率無采集時間先后差異的特點，避免數(shù)據(jù)抽頭損失，采用逆向濾波彌補抽頭的方式減少數(shù)據(jù)損失。

首先通過正向自適應(yīng)濾波算法對視電阻率數(shù)據(jù)進行濾波處理，正向自適應(yīng)濾波算法可以表示為：

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階自適應(yīng)濾波器的權(quán)重系數(shù)；E1(n)表示期望信號；e1(n)表示誤差值，單位為Ω·m；u1是濾波過程中的步長因子；K1為濾波選用的抽頭長度，用來預(yù)測的數(shù)據(jù)長度?？梢钥闯觯@得的濾波電阻率數(shù)據(jù)長度為N-K，會損失抽頭長度，也就是y1只具有除抽頭外的濾波視電阻率數(shù)據(jù)，為彌補損失的視電阻率數(shù)據(jù)，再進行反向自適應(yīng)濾波。

反向自適應(yīng)濾波算法可以表示為：

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階自適應(yīng)濾波器的權(quán)重系數(shù)；E2(m)表示當(dāng)前長度值的期望信號；e2(m)表示誤差值，單位為Ω·m；u2是步長因子；K2為抽頭長度。反向抽頭濾波會損失尾部視電阻率數(shù)據(jù)，可通過互補的方式獲得完整的視電阻率濾波數(shù)據(jù)，取Y作為最后的數(shù)據(jù)。

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

](7)

2.2 自動反演基本原理

采用自適應(yīng)濾波算法對具有不確定噪聲的視電阻率進行濾波后，采用自動反演方法對預(yù)處理后的視電阻率進行反演運算。自動反演的基本原理可歸結(jié)為尋找一種地下模型，使地下模型的視電阻率與測得的真實視電阻率盡可能接近。

設(shè)置待尋找的模型為Mj(j=1∶MN)，j為模型參數(shù)，包括每層厚度Th(單位為m)，和每層的電阻率ρi(單位為Ω·m)。為尋找模型Mj，首先設(shè)置初始模型，初始模型主要包括厚度Th，單位為m，和每層的電阻率ρi,單位為Ω·m。

使尋找的模型的視電阻率與真視電阻率滿足最小二乘法的精度要求:

(8)

為尋找模型Mj，將模型M在預(yù)測模型處展開，只保留二次項:

(9)

其中,ΔMj為修正量。通過式(9)來獲得模型的修正量，從而獲得接近真實模型的反演模型。在式(9)中要獲得修正量，則需要計算視電阻率的偏導(dǎo)值。對于數(shù)值反演可通過差分方法獲得偏導(dǎo)值:

(10)

則修正量可表示為：

(11)

獲得模型修正參數(shù)后，對修正后的模型參數(shù)進行正演計算，與測量參數(shù)進行比較，可以得到誤差值。通常情況下，一次修正量是不能滿足精度要求的。為了保證反演過程是收斂的，每次修正都限制模型參數(shù)的修正量在一定的范圍之內(nèi)。因此一般需要多次迭代獲得多次修正量，逐步對模型進行修正，直到滿足參數(shù)誤差精度要求為止，最后獲得的反演參數(shù)即為模型參數(shù)。

3 仿真驗證

為驗證視電阻率進行自適應(yīng)濾波對反演過程收斂的有效性，通過仿真分別對存在白噪聲、工頻噪聲的多模型視電阻率數(shù)據(jù)進行數(shù)據(jù)預(yù)處理前和預(yù)處理后反演結(jié)果的對比實驗。

實驗設(shè)備：自制的2.5D電法正反演仿真平臺，該仿真平臺具有如下功能：①能夠根據(jù)需要對不同形狀地下異常體進行建模，并能夠?qū)δＰ瓦M行正演，獲得視電阻率曲線；②為觀測視電阻率曲線中存在噪聲情況對反演結(jié)果的影響，該仿真平臺能夠在視電阻率曲線中加入白噪聲、工頻噪聲等，并能夠進行視電阻率自適應(yīng)濾波處理；③對正演生成的視電阻率進行反演計算。

3.1 實驗步驟

1)設(shè)置地下電阻率異常體形狀為球體具體參數(shù)如圖1所示。

圖1 地下異常體示意Fig.1 Schematic diagram of underground abnormal body

其中，球體異常區(qū)域與總區(qū)域同心，設(shè)置總區(qū)域邊長d0為3 000 cm，球體半徑r為1 000 cm。

2)正演模型的建立：對異常體模型進行正演計算，正演方式采用溫納四極高密度電法探測方式，其探測示意如圖2所示。

其中共鋪設(shè)20道電極，設(shè)置電極間距為150 cm，被探測區(qū)域為邊長3 000 cm的正方形，為能夠?qū)⒈惶綔y區(qū)域探測完全，電極分布距離須為探測深度的兩倍以上，電極分布在6 000 cm×6 000 cm的區(qū)域內(nèi)。

3)進行正演計算，獲得20條視電阻率曲線R1∶20。

4)同時向視電阻率曲線中增加白噪聲和工頻噪聲獲得兩組存在不同噪聲的視電阻率曲線RWF1∶20。

圖2 電法探測示意Fig.2 Schematic diagram of electrical detection

對存在白噪聲和工頻噪聲的視電阻率和利用自適應(yīng)濾波算法對存在噪聲的視電阻率進行濾波后的視電阻率進行反演，并記錄每次反演時間。

3.2 數(shù)據(jù)分析

根據(jù)上述實驗步驟，其球體異常區(qū)域的正演模型如圖3所示，球體半徑為1 000 cm。因為設(shè)置20道測量線，得到20組剖面圖，其中每組剖面圖間距為150 cm,實現(xiàn)2.5 D反演。設(shè)置異常體電阻為35 Ω·m，區(qū)域內(nèi)土壤電阻率為200 Ω·m。對所建模型進行正演計算，獲得20條視電阻率曲線。在視電阻率曲線中同時加入一個固定幅值白噪聲和工頻噪聲，其中白噪聲對電阻率測量結(jié)果的影響幅值設(shè)置為0.1 Ω·m，工頻噪聲波動幅值為0.1 Ω·m。

圖3 地下異常體正演模型Fig.3 Forward modeling of underground anomalous body

圖4 20道視電阻率預(yù)處理曲線Fig.4 20 apparent resistivity preprocessing channel

圖5 濾波前反演結(jié)果Fig.5 Inversion results before filtering

圖6 濾波后反演結(jié)果Fig.6 Inversion results after filtering

預(yù)處理后的20條視電阻率曲線如圖4所示。通過對20條預(yù)處理前存在噪聲的視電阻率曲線和預(yù)處理后的視電阻率曲線進行反演計算，得到2.5 D反演結(jié)果，其中存在噪聲的視電阻率反演結(jié)果如圖5所示，對噪聲進行自適應(yīng)濾波預(yù)處理后反演結(jié)果如圖6所示。

從圖5可以看出，對不進行預(yù)處理、存在白噪聲和工頻噪聲的球體異常區(qū)域的視電阻率曲線反演結(jié)果異常區(qū)域中心電阻率約為65 Ω·m，邊緣電阻率約為150 Ω·m，異常區(qū)域半徑約為1 400 cm，電阻率測量相對誤差約為10 %，異常區(qū)域半徑相對誤差約為40 %，反演時間為27 min。由圖6可知，采用自適應(yīng)濾波算法對存在噪聲的視電阻率曲線進行濾波處理后，反演結(jié)果中異常體中心的視電阻率約為35 Ω·m，異常體邊緣視電阻率約為120 Ω·m，異常體半徑約為1 100 cm，反演時間為19 min。通過比較，對存在白噪聲和工頻噪聲的進行自適應(yīng)濾波處理能夠明顯改善反演結(jié)果，并且縮短反演時間約30 %。

4 總 結(jié)

在進行高密度電法探測地下異常體時，測量獲得的視電阻率曲線中往往存在各種不確定噪聲，帶來反演過程不收斂，反演結(jié)果不準(zhǔn)確，反演時間增加等問題。針對存在的各種“不確定性”噪聲，提出采用自適應(yīng)濾波方法對存在噪聲的視電阻率進行預(yù)處理，再進行反演計算，通過仿真驗證該方法的有效性。首先構(gòu)建地下球體半徑為1 000 cm，電阻率為35 Ω·m的異常區(qū)域，在假設(shè)鋪設(shè)20道測量線的情況下，進行正演計算，獲得20組視電阻率曲線，在視電阻率曲線中加入白噪聲和工頻噪聲，再通過自適應(yīng)濾波算法對存在噪聲的視電阻率進行濾波預(yù)處理。對存在噪聲的視電阻率和濾波后的視電阻率進行2.5 D反演計算，不進行預(yù)處理，存在白噪聲和工頻噪聲的球體異常區(qū)域的視電阻率曲線反演結(jié)果，異常區(qū)域中心電阻率約為50 Ω·m，邊緣電阻率約為150 Ω·m，異常區(qū)域半徑約為1 300 cm，電阻率測量相對誤差約為10 %，異常區(qū)域半徑相對誤差約為40 %，反演時間為27 min。采用自適應(yīng)濾波算法對存在噪聲的視電阻率曲線進行濾波處理后，反演結(jié)果中異常體中心的視電阻率約為35 Ω·m，異常體邊緣電阻率約為120 Ω·m，異常體半徑約為1 100 cm，反演時間為19 min。通過比較，對存在白噪聲和工頻噪聲的進行自適應(yīng)濾波處理能夠明顯改善反演結(jié)果，并且縮短反演時間約30 %。