基于全區(qū)視電阻率瞬變電磁一維正演中濾波系數(shù)的選取

郭嵩巍, 鄭 凱, 劉小畔, 王 成, 王 淼

(1.中國地質(zhì)大學(xué)(北京)地球科學(xué)與資源學(xué)院， 北京 100083； 2.長江大學(xué)地球物理與石油資源學(xué)院， 武漢 430100；3.浙江省地礦勘探院， 杭州 310000； 4.成都理工大學(xué)地球物理學(xué)院， 成都 610059； 5.內(nèi)蒙古國土資源勘查開發(fā)有限責(zé)任公司, 呼和浩特 010020)

瞬變電磁法(transient electromagnetic method，TEM)是一種較為成熟的地球物理電法勘探方法，被廣泛應(yīng)用于礦產(chǎn)、水文、工程等領(lǐng)域[1-2]。在工程應(yīng)用領(lǐng)域，一維反演依然是目前瞬變電磁法資料解釋的重要依據(jù)之一[3-4]。然而，瞬變電磁正演計算效率低下是制約反演應(yīng)用的重要因素，為提高正演的計算效率，選取合適的濾波系數(shù)是關(guān)鍵。

瞬變電磁法裝置多樣，發(fā)射場源包括磁性源、電性源[5-8]，接收裝置普遍采用回線線框，采集的數(shù)據(jù)大多為感應(yīng)電動勢，可解析得到晚期視電阻率。由于晚期視電阻率在早期段，即測深曲線的首支存在明顯的高阻假象，而全區(qū)視電阻率很好地解決了這個問題，對地質(zhì)體電性表達(dá)更為真實，因此轉(zhuǎn)換全區(qū)視電阻率逐漸成為處理解釋的主流。全區(qū)視電阻率一般可通過數(shù)值計算得出，無論是磁性源還是電性源，很多學(xué)者都將感應(yīng)電動勢轉(zhuǎn)換成垂直磁場，進(jìn)而計算出全區(qū)視電阻率[9-12]。因此以全區(qū)視電阻率作為濾波系數(shù)選取研究的正演參數(shù)。

通過對麥克斯韋方程組的推導(dǎo)，可得出瞬變電磁一維正演響應(yīng)的解析表達(dá)式。其中核函數(shù)的貝塞爾函數(shù)積分問題，通常采用漢克爾變換[13-15]；頻率域到時間域的轉(zhuǎn)換方面，有G-S變換[16-17]、余弦變換[18]、逆Laplace變換[19]等，本文采用余弦變換。漢克爾變換和余弦變換數(shù)值算法均采用數(shù)值濾波算法，選擇合適的濾波系數(shù)決定著瞬變電磁正演的計算效率、精度以及有效范圍。分別選取正演效果較好的三組濾波系數(shù)，對比分析102Ω·m均勻半空間模型正演響應(yīng)，得出全區(qū)視電阻的有效范圍，并以該模型為參照模型，利用瞬變電磁響應(yīng)的平移伸縮特性[20]，提出參照模型法，以為工程應(yīng)用的正反演計算中選取濾波系數(shù)，提供一個簡單可行的估算方法。

1 選取濾波系數(shù)問題的由來

以最為典型的中心回線裝置為例，一維正演的數(shù)值計算可簡化為核函數(shù)的兩層積分，內(nèi)層積分是貝塞爾函數(shù)積分問題，采用漢克爾變換，用Ht{}表示；外層積分是電磁場響應(yīng)從頻率域到時間域的轉(zhuǎn)換問題，本文采用余弦變換，用Ct{}表示；核函數(shù)用fc表示，瞬變響應(yīng)的垂直磁場用hz表示，于是中心回線瞬變電磁[21]一維正演響應(yīng)可簡化表示為

(1)

(2)

式中：I0為發(fā)射電流；a為發(fā)射線框半徑；t為正演自變量時間；u(1)為頂層波阻抗；λ為核函數(shù)漢克爾變換無窮積分的自變量；ω為余弦變換的自變量頻率。

中心回線瞬變電磁一維正演響應(yīng)進(jìn)一步簡化為

(3)

對于均勻半空間地電模型，可以推導(dǎo)出垂直磁場hz的解析表達(dá)式為

(4)

式(4)中：erf(u)為誤差函數(shù)。

歸一化hz得到

(5)

于是就將瞬變響應(yīng)轉(zhuǎn)換為基于垂直磁場的全區(qū)視電阻率，即

(6)

式(6)中：μ0為真空中磁導(dǎo)率。

由于函數(shù)Z(u)為隱函數(shù)，無法解析導(dǎo)出反函數(shù)，但可利用函數(shù)在0～1單調(diào)遞增的特點得到結(jié)果通過二分查找的辦法實現(xiàn)了該全區(qū)視電阻率的數(shù)值計算，并與其他方法做了對比分析，得出了該全區(qū)視電阻反演效果更好的結(jié)論[9]。

于是瞬變電磁的正演數(shù)值算法的核心問題，就轉(zhuǎn)變?yōu)闈h克爾變換和余弦變換的數(shù)值濾波計算問題。一般來說，濾波系數(shù)采樣點越多，正演的計算精度越高，正演響應(yīng)的有效范圍就越大，但計算效率就越低。若濾波系數(shù)不足或不符合正演響應(yīng)有效范圍要求，將不能模擬準(zhǔn)確完整的正演曲線，反演時會造成擬合失真，再做解釋將是災(zāi)難性的。因此，有必要擇優(yōu)選取濾波系數(shù)。

2 數(shù)值濾波算法

數(shù)值濾波算法，通常是為了解決無窮積分而設(shè)計的，漢克爾變換解決的是含有貝塞爾函數(shù)的無窮積分，余弦變換是對傅里葉變換的化簡，實現(xiàn)了從頻率域到時間域轉(zhuǎn)換的無窮積分。

濾波系數(shù)的得出是通過變量替換，將無窮積分式變成卷積的形式。在滿足抽樣定理的條件下，通過求解卷積方程得到指數(shù)間隔采樣的濾波系數(shù)，即反算濾波器，方法有傅里葉變換法、維納-霍普夫極小化法[22]等。

數(shù)值濾波算法通常形式為

(7)

(8)

式中：B和W為一組濾波系數(shù)變換對，數(shù)組長度為n，B一般采用指數(shù)間隔采樣，λ為數(shù)值濾波算法中卷積形式的采樣點，在漢克爾變換中對應(yīng)貝塞爾函數(shù)的采樣點，在余弦變換中對應(yīng)頻率ω的采樣點。

推薦漢克爾變換為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 參照模型法

由于實測數(shù)據(jù)的發(fā)射裝置不同、對應(yīng)地電模型電阻率不同、采集數(shù)據(jù)響應(yīng)時間的上下限不同，要選取合適的濾波系數(shù)，就必須將“不同”變成“相同”，因此提出參照模型法。該方法是一種簡單的估算方法，可將發(fā)射裝置不同、模型電阻率不同、響應(yīng)時間上下限不同的實測數(shù)據(jù)，轉(zhuǎn)換到對應(yīng)參照模型的時間范圍，再與參照模型不同濾波系數(shù)組合正演模擬的有效范圍做比對，從而做到對濾波系數(shù)的選取。

這里需要說明的是，本文參照模型的有效范圍是中心回線裝置的計算結(jié)果，其他裝置的有效范圍可通過相應(yīng)正演模擬計算得出。

現(xiàn)有地質(zhì)地球物理資料顯示，地殼內(nèi)大多數(shù)地質(zhì)體的電阻率范圍在1～104Ω·m，取對數(shù)中間值102Ω·m作為正演參照模型的電阻率(ρr)。發(fā)射裝置方面，煤田首采區(qū)勘查多采用邊長3×102～103m的大回線框掃面；水文勘查多采用102～3×102m中心回線裝置；工程勘查常采用多扎小線框，邊長n～n×10 m；取對數(shù)中間值102m作為參照模型的發(fā)射線框邊長，參照發(fā)射線框面積(Sr)為104m2。

根據(jù)瞬變電磁的平移特性

(11)

(12)

得到實測數(shù)據(jù)時間與其參照模型時間的關(guān)系：

(13)

對式(13)取對數(shù)，得

lgtr=lgts+lg(ρs/ρr)+lg(Sr/Ss)

(14)

式中：Sr為參照模型發(fā)射線框的面積，104m2；ρr為參照模型電阻率，102Ω·m；tr為參照模型時間；Ss為實測數(shù)據(jù)發(fā)射線框的面積，104m2；ρs為實測數(shù)據(jù)模型電阻率，102Ω·m，如假設(shè)為反演模型的電阻率上下限范圍等；ts為實測數(shù)據(jù)的時間，s。

由于實測數(shù)據(jù)并非對應(yīng)均勻半空間模型，可根據(jù)全區(qū)視電阻率的范圍推斷實測數(shù)據(jù)的模型電阻率范圍，或者采用反演模型約束電阻率的上下限范圍估算，總之需要對實測數(shù)據(jù)對應(yīng)的模型電阻率估算一個上下限范圍。令

(15)

c2=lg(Sr/Ss)

(16)

最終得到

(17)

4 不同濾波系數(shù)組合的有效范圍

經(jīng)過對多組漢克爾變換濾波系數(shù)、余弦變換濾波系數(shù)的對比試驗，擇優(yōu)選出3組漢克爾變換系數(shù)和3組余弦變換系數(shù)，計算參照模型的正演響應(yīng)，并轉(zhuǎn)換為全區(qū)視電阻率，設(shè)定相對誤差在10%以內(nèi)為有效數(shù)據(jù)，得到有效范圍(表1)。余弦變換固定81點濾波系數(shù)，漢克爾變換分別對Guptasama 47點、201點、401點[24]濾波系數(shù)做正演，對比結(jié)果見圖1、表1(組合編號1、2、3)；漢克爾變換固定Guptasarma47點濾波系數(shù)，采用3組不同的余弦變換分別采用21點、81點、601點濾波系數(shù)[24]，對比結(jié)果見圖2、表1(組合編號4、1、5)。

表1 參照模型有效時間范圍統(tǒng)計表Table 1 Statistical table of effective time range of reference model

圖1 3組漢克爾變換濾波系數(shù)對比Fig.1 Comparison of the Hankel transform filters coefficient

圖2 3組余弦變換濾波系數(shù)對比Fig.2 Comparison of the 3 groups cosine transform filter coefficient

5 實測數(shù)據(jù)的參照模型時間范圍

根據(jù)前面推導(dǎo)的參照模型法，分別列出了內(nèi)蒙古自治區(qū)不同地區(qū)的不同工區(qū)、不同發(fā)射裝置、不同地電模型的實測數(shù)據(jù)(全部采集自V8電法工作站)，估算出對應(yīng)參照模型的時間范圍。由于變化范圍大，用對數(shù)值lgt表示,見表2。

與(表1)比對可以發(fā)現(xiàn)，通過參照模型法的估算，本文提供的濾波系數(shù)組合的有效時間范圍，全部滿足實測數(shù)據(jù)的計算需求。

6 數(shù)據(jù)處理軟件模塊的開發(fā)

擇優(yōu)選取濾波系數(shù)，目的是為了更好地實現(xiàn)正演建模和對實測數(shù)據(jù)的反演計算，因此基于VB2015開發(fā)了瞬變電磁數(shù)據(jù)處理軟件模塊，實現(xiàn)了數(shù)據(jù)讀取、顯示、編輯，全區(qū)視電阻率轉(zhuǎn)換、煙圈反演、一維反演等功能(圖3)。作為正演、反演計算的設(shè)置參數(shù)，程序中可根據(jù)實際情況需要，選擇濾波系數(shù)的任意組合，其中組合4計算效率最高，“參照模型法”估算其可以滿足大部分實測數(shù)據(jù)的計算需求(表1、表2)，因此程序?qū)⒔M合4設(shè)置為缺省值。

圖3 數(shù)據(jù)處理模塊軟件界面Fig.3 Data processing software module interface

7 實測剖面算例

為證明“參照模型法”估算選取濾波系數(shù)的可行性，選擇內(nèi)蒙古自治區(qū)涼城縣某地?zé)峥辈轫椖康囊粭l剖面的實測數(shù)據(jù)作為測試數(shù)據(jù)，反演算法選擇正則化擬二維反演(ARIA)[25-26]，該項目施工條件良好，數(shù)據(jù)采集質(zhì)量高，最終在北東向和北西向構(gòu)造交匯處布置地?zé)峋?DR2)，取得了不錯的找水效果[27]。

首先，實測數(shù)據(jù)剖面位于地?zé)徙@孔(DR2)附近，地層電性分層明確，基底埋深確定，可確保反演結(jié)果的正確性。其次，該數(shù)據(jù)發(fā)射線框600 m×400 m，采集時間8.9×10-5～2.1×10-2s，按照“參照模型法”估算的時間范圍(表2)，對數(shù)值在-3.67～0.40，范圍較大，有利于測試各個濾波系數(shù)組合的有效范圍是否滿足實際需求，若有效范圍偏小，對曲線尾支的擬合影響較大，在反演斷面圖上表現(xiàn)為基底偏差較大，不同濾波系數(shù)組合呈現(xiàn)出不同的反演結(jié)果，因此該數(shù)據(jù)具備對濾波系數(shù)選取的測試條件。

表2 實測數(shù)據(jù)轉(zhuǎn)換參照模型的時間范圍Table 2 The transform of reference model time range from measured data

7.1 地質(zhì)概況

7.1.1 地層

工作區(qū)地層主要有太古界桑干群、中生界侏羅系及新生界古近系、新近系、第四系。由老至新分述如下：

太古界桑干群(Ar1Sg)，主要分布于蠻漢山區(qū)及平原區(qū)深部，為一套黑云矽線榴石斜長片麻巖、鉀長片麻巖與淺色麻粒巖。

中生界侏羅系上侏羅統(tǒng)(J3)，零星分布于蠻漢山區(qū)大西溝一帶，為一套巨厚層狀粉白色、灰紫色、灰白色凝灰?guī)r、凝灰角礫巖、凝灰集塊巖及安山巖組成，傾向北西，不整合覆于桑干群片麻巖之上。

新生界可劃分為古近系與新近系：

古近系(E)：見于鉆孔，為一套深棕紅色泥巖、泥質(zhì)砂巖、砂質(zhì)泥巖夾灰白色、深棕紅色砂礫巖，不整合于太古界地層之上。

新近系(N)：見于鉆孔，上部為淺棕紅色泥巖、泥質(zhì)砂巖、砂質(zhì)泥巖夾灰黃色砂礫巖，下部為土黃色、灰白色砂礫巖，與古近系地層呈假整合接觸。

新生界第四系(Qal+pl)：主要分布于蠻漢山山前傾斜平原，巖性為淺黃棕色、黃褐色粉土、粉質(zhì)黏土與淺黃褐色砂礫卵石、含卵礫中粗砂、中粗砂等。粉質(zhì)黏土或粉砂土中含有砂與礫石；砂礫卵石分選差，多棱角狀，成分以片麻巖、凝灰?guī)r為主。向平原南部巖性變細(xì)，為黃褐色粉土、粉質(zhì)黏土與砂礫石。

7.1.2 構(gòu)造

工作區(qū)區(qū)域地質(zhì)構(gòu)造(圖4)屬陰山東西向構(gòu)造帶的南緣。因受早期構(gòu)造運動的影響，呈北東30°～40°延伸，系高角度壓性斷層，而北西向斷裂據(jù)淺層地震資料所反映為一條北西325°的破碎帶，寬度為100～200 m。

圖4 工作區(qū)地質(zhì)圖Fig.4 Geological map of work area

7.1.3 巖漿巖

主要分布在蠻漢山區(qū)及平原區(qū)深部，為太古代早期侵入巖,巖性以蘇長巖和似斑狀花崗巖為主。

7.2 瞬變電磁反演效果對比

測試剖面800線,點距10 m，測點共51個，剖面總長500 m，位于山前傾斜平原。地?zé)徙@孔(DR2)位于剖面附近，鉆孔顯示地層分層明確、巖性差異明顯，存在明顯電性差異：淺部第四系(Q)砂巖黏土總體呈高阻,電阻率為30～80 Ω·m；中間層為古近系(E)、新近系(N)泥巖地層，為低阻層，電阻率為10～30 Ω·m，頂板埋深約150 m；基底花崗巖高阻，電阻率大于100 Ω·m，頂板埋深大于300 m。

反演結(jié)果(圖5)顯示，不同濾波系數(shù)組合的反演結(jié)果基本一致。電性分層與DR2鉆孔基本吻合，說明決定反演效果主要因素是正演響應(yīng)的有效范圍，正演響應(yīng)的精度(相對誤差)對反演結(jié)果的影響相對較小。因此認(rèn)為正演的有效范圍和計算效率更具實際意義(表3)，對測點較多、計算量較大的工程應(yīng)用可影響反演的效率和效果。

圖5 不同濾波系數(shù)組合方案正則化反演結(jié)果對比Fig.5 Comparison of ARIA inversion results of different filter coefficient selection schemes

表3 實測剖面反演效率統(tǒng)計表Table 3 Statistics of inversion efficiency of measured profile

8 結(jié)論

(1)瞬變電磁一維正演響應(yīng)的計算效果(有效范圍)與濾波系數(shù)的取值相關(guān)，計算效率與濾波系數(shù)的數(shù)量相關(guān)，成正相關(guān)關(guān)系。共推薦并測試了正演效果良好的三組漢克爾變換和余弦變換的濾波系數(shù)。

(2)通常情況下，計算效果好的濾波系數(shù)組合普遍積分區(qū)間較長，計算效率不高，可通過參照模型法估算出實測數(shù)據(jù)反演的計算要求，做到合理選取濾波系數(shù)組合。

(3)對實測剖面的反演測試顯示，有效范圍并非越大越好，正演響應(yīng)(全區(qū)視電阻率)有效范圍大犧牲的是計算效率，而計算精度的提升對于反演結(jié)果的影響較小，因此對于實際工程應(yīng)用而言，選取濾波系數(shù)的原則是滿足參照模型有效時間范圍的情況下，盡可能選擇濾波系數(shù)少，計算效率高的濾波系數(shù)組合。

綜上所述，推薦濾波系數(shù)組合4，即漢克爾變換Gaptarsma47點、余弦變換Kong21點，為目前瞬變電磁實測數(shù)據(jù)反演計算的最佳濾波系數(shù)組合。