Matlab在CSAMT數(shù)據(jù)“去噪”及近場(chǎng)校正中的應(yīng)用

胡慶輝

(1.山東省地質(zhì)礦產(chǎn)勘查開(kāi)發(fā)局 第四地質(zhì)大隊(duì)，山東 濰坊 261021; 2.山東省地礦局 海岸帶地質(zhì)環(huán)境保護(hù)重點(diǎn)實(shí)驗(yàn)室，山東 濰坊 261021)

0 引言

可控源音頻大地電磁測(cè)深法(CSAMT),是在大地地磁測(cè)深法(MT)和音頻大地電磁測(cè)深法(AMT)的基礎(chǔ)上發(fā)展起來(lái)的一種人工源頻率域測(cè)深方法[1]。自上世紀(jì)70年代引入國(guó)內(nèi)，在油氣勘查、金屬礦普查、水文及環(huán)境地質(zhì)等多個(gè)領(lǐng)域顯示出廣闊的應(yīng)用前景。隨著經(jīng)濟(jì)社會(huì)發(fā)展，CSAMT野外數(shù)據(jù)采集所面臨的“電磁環(huán)境”越發(fā)復(fù)雜，勘探信號(hào)不可避免地受到“污染”，導(dǎo)致最終反演數(shù)據(jù)質(zhì)量不高，加大了資料解譯的難度甚至造成失誤。此外，受限于CSAMT方法的原理及假設(shè)前提，在中、低頻率上電磁場(chǎng)進(jìn)入“近區(qū)”或者“過(guò)渡區(qū)”，卡尼亞視電阻率將發(fā)生畸變，無(wú)法實(shí)現(xiàn)勘查目的[2]。底青云等[3]認(rèn)為CSAMT資料處理的質(zhì)量依賴于近場(chǎng)校正的質(zhì)量。

國(guó)內(nèi)學(xué)者針對(duì)CSAMT方法的電磁干擾壓制及近場(chǎng)校正技術(shù)進(jìn)行了相關(guān)研究。偽隨機(jī)編碼技術(shù)在供電環(huán)節(jié)的使用能夠有效提高野外觀測(cè)數(shù)據(jù)信噪比[4-5]，而雷達(dá)等采用信息熵和有理函數(shù)濾波相結(jié)合的方式對(duì)干擾數(shù)據(jù)進(jìn)行處理也取得到比較好的效果[6]。通常的近場(chǎng)校正的方法是“三角形校正法”[3，7-9]，而采用“全頻域視電阻率法”進(jìn)行近場(chǎng)校正得到越來(lái)越廣泛的關(guān)注和應(yīng)用[10-11]。

用戶界面(或接口)是指人與機(jī)器(或程序)之間交互作用的工具和方法[12]。筆者采用Matlab GUI編程，實(shí)現(xiàn)友好的可視化交互式操作界面，剔除“畸變”頻點(diǎn)數(shù)據(jù)，實(shí)現(xiàn)CSAMT數(shù)據(jù)的“去噪”處理。采用數(shù)值計(jì)算功能，編寫(xiě)了基于全頻域視電阻率計(jì)算的程序以實(shí)現(xiàn)近場(chǎng)校正。程序的應(yīng)用提高了CSAMT數(shù)據(jù)的處理效率和效果，保證后續(xù)處理的數(shù)據(jù)質(zhì)量，在實(shí)際生產(chǎn)中得到應(yīng)用并取得了很好的效果。

圖1 CSAMT野外實(shí)測(cè)干擾數(shù)據(jù)典型特征Fig.1 Typical feature of interfering survey data of CSAMT(a)卡尼亞視電阻率曲線；(b)阻抗相位曲線；(c)電場(chǎng)分量幅值曲線；(d)磁場(chǎng)分量幅值曲線

1 CSAMT干擾及“近場(chǎng)”數(shù)據(jù)特征

CSAMT野外數(shù)據(jù)環(huán)節(jié)的電磁干擾主要來(lái)自天然電磁噪聲和人文電磁噪聲[13]，其中天然電磁噪聲主要來(lái)自雷電活動(dòng)，其引起的電磁場(chǎng)在CSAMT常用工作頻段內(nèi)能量較高。這類(lèi)噪聲一般可采用多次疊加進(jìn)行壓制。人文電磁噪聲主要包括工業(yè)用電、大功率發(fā)射臺(tái)以及隨機(jī)的電磁脈沖干擾(如動(dòng)力電開(kāi)關(guān)轉(zhuǎn)換、鐵磁性物體靠近磁探頭等)。固定頻率的電磁干擾可以通過(guò)設(shè)置陷波器進(jìn)行壓制，而能量較高、隨機(jī)的強(qiáng)干擾還不能完全克服，需要野外施工人員正確判斷電磁干擾特點(diǎn)，采取必要的規(guī)避措施，盡量避開(kāi)干擾源或者錯(cuò)開(kāi)干擾時(shí)段，以便采集到相對(duì)高質(zhì)量的野外數(shù)據(jù)。

隨機(jī)電磁干擾在CSAMT實(shí)測(cè)數(shù)據(jù)中主要表現(xiàn)為電場(chǎng)及磁場(chǎng)分量復(fù)測(cè)離差增大，單頻點(diǎn)數(shù)據(jù)離差加大，數(shù)據(jù)一致性較差，阻抗相位偏差較大，甚至出現(xiàn)反相，計(jì)算得出的卡尼亞電阻率出現(xiàn)“飛點(diǎn)”，曲線形態(tài)不能反映真實(shí)大地電阻率的變化特征(圖1)。

基于上述兩個(gè)方面的問(wèn)題，筆者利用Matlab數(shù)值計(jì)算及GUI功能，設(shè)計(jì)了針對(duì)野外數(shù)據(jù)預(yù)處理和近場(chǎng)校正的相關(guān)程序，在CSAMT數(shù)據(jù)“去噪”預(yù)處理及近場(chǎng)校正等方面進(jìn)行了探索，取得較好的應(yīng)用效果。

2 CSAMT“去噪”及“近場(chǎng)校正”方法

現(xiàn)階段，可控源音頻大地電磁方法(CSAMT)“去噪”多采用濾波方法，在干擾比較小的情況下，一般能夠得到比較好的效果。但隨著干擾的加強(qiáng)，濾波方法往往難以取得理想的電磁噪聲壓制效果。

本次“去噪”處理思路比較直觀，即從多次觀測(cè)數(shù)據(jù)中挑選“合理”數(shù)據(jù)，舍棄畸變數(shù)據(jù)，從而達(dá)到“去偽存真”的目的。所謂“合理”，即從電極布設(shè)及采集時(shí)就保證觀測(cè)質(zhì)量，合理避開(kāi)干擾源或干擾地段，錯(cuò)開(kāi)干擾時(shí)段進(jìn)行數(shù)據(jù)采集，使多數(shù)觀測(cè)數(shù)據(jù)為真實(shí)地電反映，其數(shù)據(jù)分布應(yīng)具有“集中分布”特征，重復(fù)性較好。而少數(shù)觀測(cè)值由于隨機(jī)干擾產(chǎn)生“畸變”，相關(guān)參數(shù)如相位、電場(chǎng)分量幅值、磁場(chǎng)分量幅值等偏離合理區(qū)間，出現(xiàn)明顯的“飛點(diǎn)”特征，通過(guò)友好界面對(duì)畸變數(shù)據(jù)予以剔除，從而實(shí)現(xiàn)“去噪”處理。

圖2 CSAMT數(shù)據(jù)預(yù)處理GUI界面Fig.2 Matlab GUI interface for data pre-processing of CSAMT

數(shù)據(jù)處理通過(guò)鼠標(biāo)在界面上框選畸變數(shù)據(jù)實(shí)現(xiàn)便捷操作(圖2)。用戶界面分為三部分：①工具欄；②繪圖區(qū)；③信息提示區(qū)。工具區(qū)包括數(shù)據(jù)加載、重置、顯示參數(shù)選項(xiàng)、數(shù)據(jù)操作、數(shù)據(jù)保存等功能按鍵或下拉菜單。繪圖區(qū)則顯示相應(yīng)選項(xiàng)的散點(diǎn)圖、擬合曲線圖或者擬斷面等值線圖。界面下方為信息提示區(qū)，顯示部分文件信息以及當(dāng)前點(diǎn)號(hào)等。

加載CSAMT數(shù)據(jù)完成后，繪圖區(qū)默認(rèn)顯示首個(gè)測(cè)點(diǎn)的卡尼亞電阻率數(shù)據(jù)，綠色圓點(diǎn)為某個(gè)頻率重復(fù)觀測(cè)數(shù)據(jù)點(diǎn)，藍(lán)色曲線為其樣條插值擬合曲線?？梢钥闯觯陔姶鸥蓴_條件下采集的數(shù)據(jù)存在明顯的離散特征，即重復(fù)觀測(cè)數(shù)據(jù)一致性較差，“飛點(diǎn)”較多。同時(shí)，每個(gè)頻率的觀測(cè)數(shù)據(jù)點(diǎn)有重復(fù)性較好或相對(duì)較為集中的特征，我們認(rèn)為這樣的數(shù)據(jù)為“合理”數(shù)據(jù)。用戶使用鼠標(biāo)選擇工具欄“刪除”按鍵，在繪圖區(qū)框選需要剔除的數(shù)據(jù)點(diǎn)(圖2中紅色框及紅色點(diǎn))，它們將不參與后續(xù)的數(shù)據(jù)處理。利用下拉菜單選擇“阻抗相位”、“電場(chǎng)分量幅值”及“磁場(chǎng)分量幅值”三個(gè)選項(xiàng)在相應(yīng)數(shù)據(jù)散點(diǎn)上完成類(lèi)似操作，實(shí)現(xiàn)多個(gè)參數(shù)對(duì)同一個(gè)測(cè)點(diǎn)的“去噪”處理。當(dāng)該測(cè)點(diǎn)數(shù)據(jù)選擇完畢后，通過(guò)“保存”按鍵實(shí)現(xiàn)AVG數(shù)據(jù)保存，為后續(xù)處理提供高信噪比數(shù)據(jù)。

CSAMT數(shù)據(jù)近場(chǎng)校正是在上述AVG數(shù)據(jù)的基礎(chǔ)上，利用“等效電阻率全頻域視電阻率法”，采用迭代擬合方式計(jì)算等效電阻率。

通常采用赤道偶極方式進(jìn)行CSAMT數(shù)據(jù)采集時(shí)，采集參數(shù)為電場(chǎng)和磁場(chǎng)的水平分量Ex、Hy，卡尼亞視電阻率可以表示為式(1)。

(1)

(2)

利用ρs是P的單調(diào)下降函數(shù)，通過(guò)實(shí)測(cè)ρs確定P值，并通過(guò)式(3)計(jì)算大地電阻率。

(3)

Matlab數(shù)值計(jì)算的實(shí)現(xiàn)過(guò)程中，采用迭代方式求取F(P)=0的根。具體思路為：給定取值范圍為[a,b] ，給P賦初值c=(a+b)/2，判斷F(a).F(c)符號(hào)，根據(jù)函數(shù)的單調(diào)下降性質(zhì)，采用二分法逐步逼近零點(diǎn)，從而求取P值進(jìn)而得到大地電阻率ρ。

近場(chǎng)校正處理通過(guò)Matlab GUI編程實(shí)現(xiàn)友好用戶界面進(jìn)行操作。通過(guò)主界面“近場(chǎng)及靜校正”按鍵，可以調(diào)用NFC(近場(chǎng)校正)模塊(圖3)。NFC界面包含工具欄及繪圖區(qū)兩部分，工具欄主要包括數(shù)據(jù)加載、數(shù)據(jù)列表、數(shù)據(jù)編輯、近場(chǎng)校正及數(shù)據(jù)保存按鍵。在加載數(shù)據(jù)環(huán)節(jié)需要同時(shí)加載CSAMT實(shí)測(cè)數(shù)據(jù)(*.AVG)、觀測(cè)點(diǎn)坐標(biāo)文件(*.stn)，用于計(jì)算每個(gè)測(cè)點(diǎn)的收發(fā)距及相對(duì)角度，計(jì)算結(jié)果將在相應(yīng)繪圖區(qū)顯示。繪圖區(qū)顯示測(cè)點(diǎn)卡尼亞視電阻率曲線和阻抗相位曲線，用戶可以通過(guò)進(jìn)一步的編輯功能對(duì)曲線進(jìn)行處理，調(diào)整曲線形態(tài)，使其更符合“非波區(qū)”曲線特征。計(jì)算完成后則顯示近場(chǎng)校正前后曲線對(duì)比情況，使處理人員能夠初步了解校正效果。

圖3 CSAMT數(shù)據(jù)近場(chǎng)校正GUI界面Fig.3 Matlab GUI interface for near-field correction of CSAMT

3 近場(chǎng)數(shù)據(jù)正演模擬計(jì)算

3.1 均勻半空間地電模型

均勻半空間地電模型電阻率為1 000 Ω·m，正演計(jì)算收發(fā)距r=6 km及r=14 km兩種情況下的卡尼亞視電阻率特征(圖4)?？梢钥闯?，高頻段視電阻率能夠正確反映模型電阻率特征，當(dāng)進(jìn)入到中低頻率時(shí)，曲線發(fā)生畸變，近45°上揚(yáng)，嚴(yán)重偏離模型電阻率。采用“等效電阻率全頻域視電阻率法”校正后，全頻域視電阻率曲線已經(jīng)基本水平，較準(zhǔn)確反映了模型電阻率數(shù)值，說(shuō)明利用全頻域視電阻率法進(jìn)行CSAMT近場(chǎng)校正效果非常理想。

圖4 均勻半空間CSAMT正演數(shù)據(jù)近場(chǎng)校正處理曲線對(duì)比Fig.4 Near-field correction curve contrast of CSAMT forward data in the homogeneous half-space model

圖5 層狀介質(zhì)CSAMT正演數(shù)據(jù)近場(chǎng)校正處理曲線Fig.5 Near-field correction curve of CSAMT forward data in the layered model(a)H型測(cè)深曲線；(b)A型測(cè)深曲線；(c)K型測(cè)深曲線；(d)Q型測(cè)深曲線

圖6 某測(cè)線CSAMT 1900點(diǎn)去噪及近場(chǎng)校正處理效果Fig.6 Data denoise and near-field correction process results of survey point 1900(a)實(shí)測(cè)卡尼亞電阻率數(shù)據(jù)散點(diǎn)及平均曲線；(b)實(shí)測(cè)阻抗相位數(shù)據(jù)散點(diǎn)及平均曲線；(c)“去噪”及近場(chǎng)校正后曲線對(duì)比；(d)“去噪”后阻抗相位數(shù)據(jù)點(diǎn)及平均曲線

圖7 CSAMT測(cè)量L4線視電阻率斷面圖(預(yù)處理前)Fig.7 Apparent resistivity section of CSAMT survey of line 4(before preprocessing)

3.2 層狀地電模型

層狀地電模型采用三層模型，分別設(shè)置了H型、A型、K型及Q型四種測(cè)深曲線類(lèi)型(圖5)。在中低頻段，各地電模型正演曲線進(jìn)入“近區(qū)”，曲線形態(tài)發(fā)生畸變。采用“全頻域視電阻率法”校正后，測(cè)深曲線形態(tài)基本反映了真實(shí)的地電模型特征。

同時(shí)也應(yīng)該注意到，在均勻半空間模擬中，過(guò)渡區(qū)出現(xiàn)局部“下凹”，校正值偏離模型電阻率。對(duì)于層狀模型，尤其是K型測(cè)深曲線，中間高阻層雖然趨勢(shì)較明顯，但其值在校正后明顯偏低。對(duì)于過(guò)渡區(qū)“校正不足”現(xiàn)象，筆者也做了探索性工作，主要思路是在中低頻段數(shù)據(jù)近場(chǎng)校正時(shí)增加修正系數(shù)。根據(jù)大量正演模型試算數(shù)據(jù)，對(duì)近區(qū)和過(guò)渡區(qū)近場(chǎng)校正修正系數(shù)采用多項(xiàng)式擬合。修正系數(shù)公式見(jiàn)式(4)及式(5)。

近區(qū)修正系數(shù)：

K(f)=-0.0071f3+0.0773f2-

0.2017f+1.0634

(4)

過(guò)渡區(qū)修正系數(shù)：

K(f)=-0.0018f3+0.0111f2-

0.0164f+1.6322

(5)

其中，f為校正點(diǎn)頻率。在進(jìn)行近場(chǎng)校正時(shí)，合理選擇近區(qū)及過(guò)渡區(qū)頻率，并在校正后數(shù)值基礎(chǔ)上進(jìn)行修正。該方法在正演數(shù)據(jù)中得到了很好的校正效果，過(guò)渡區(qū)“下凹”現(xiàn)象得到改善，近區(qū)數(shù)據(jù)也得到適當(dāng)調(diào)整，在野外實(shí)測(cè)數(shù)據(jù)中的應(yīng)用效果有待進(jìn)一步驗(yàn)證。

圖8 CSAMT測(cè)量L4線視電阻率斷面圖(預(yù)處理后)Fig.8 Apparent resistivity section of CSAMT survey of line 4(after preprocessing)

4 處理實(shí)例及效果分析

某高速公路連接線擬建路線穿越煤炭采空區(qū)，為確保擬建項(xiàng)目的安全，在工作區(qū)內(nèi)開(kāi)展CSAMT測(cè)量工作，查清采空區(qū)位置、分布范圍等特征。觀測(cè)裝置采用赤道偶極標(biāo)量方式，采集參數(shù)為Ex和Hy，發(fā)射偶極AB=1 200 m，測(cè)量偶極MN=50 m，最大收發(fā)距為6.4 km，測(cè)量頻段為1 Hz～8 192 Hz。由于工作場(chǎng)地為城市近郊，村莊廠房密集且有高速公路和多條高壓線穿過(guò)，電磁干擾非常嚴(yán)重。在施工過(guò)程中采取多種措施保證野外采集數(shù)據(jù)質(zhì)量，如加大供電電流、增加重復(fù)觀測(cè)次數(shù)、合理避開(kāi)干擾時(shí)段等。

以某測(cè)線1900點(diǎn)的數(shù)據(jù)散點(diǎn)及卡尼亞電阻率和阻抗相位曲線為例(圖6)。圖6(a)、圖6(b)分別為卡尼亞電阻率及阻抗相位實(shí)測(cè)數(shù)據(jù)散點(diǎn)及平均值曲線，由圖6(a)、圖6(b)可以看出，各個(gè)頻段的數(shù)據(jù)點(diǎn)均存在重復(fù)性差、離散度高等特點(diǎn)，說(shuō)明該測(cè)點(diǎn)在采集過(guò)程中受到較強(qiáng)烈電磁干擾，造成數(shù)據(jù)信噪比下降。中低頻段進(jìn)入“近區(qū)”，視電阻率曲線以近45°上揚(yáng)，且在32 Hz處出現(xiàn)局部極小值，應(yīng)為過(guò)渡區(qū)反映。阻抗相位曲線也顯示了近區(qū)及過(guò)渡區(qū)信息。

圖6(c)、圖6(d)為數(shù)據(jù)處理后數(shù)據(jù)散點(diǎn)及曲線特征。利用預(yù)處理程序，剔除離散度大的點(diǎn)，保留重復(fù)性好、相對(duì)集中分布的點(diǎn)作為“合理”的實(shí)測(cè)數(shù)據(jù)，曲線形態(tài)得到適當(dāng)調(diào)整，數(shù)據(jù)信噪比得到提高。采用處理后的數(shù)據(jù)進(jìn)行近場(chǎng)校正(圖6(c)紅色曲線)，選取1 Hz～8 Hz為近區(qū)頻段，16 Hz及32 Hz為過(guò)渡區(qū)頻段，增加相應(yīng)的修正系數(shù)，采用“全頻域視電阻率法”進(jìn)行計(jì)算，可以看出曲線形態(tài)發(fā)生較大變化，低頻段呈45°上揚(yáng)段斜率變小，過(guò)渡區(qū)電阻率抬升，高頻段視電阻率沒(méi)有發(fā)生變化，曲線整體形態(tài)更加合理，便于反演解釋。

選取L4線進(jìn)行平滑濾波二維反演，對(duì)處理前后的數(shù)據(jù)進(jìn)行對(duì)比，以檢驗(yàn)該程序的處理效果。L4線東西向布設(shè)，穿過(guò)煤炭采區(qū)，沿線多為廠區(qū)，并與高壓輸電線和鐵路相交，電磁干擾較強(qiáng)烈。采用原始數(shù)據(jù)進(jìn)行反演(圖7)，視電阻率等值線“下拉”現(xiàn)象明顯，剖面后段約-400 m深度近水平展布的低阻等值線出現(xiàn)波狀起伏且埋深偏淺，新太古代泰山巖群片麻巖地層起伏形態(tài)也發(fā)生改變。尤為明顯的是石24井位置埋深-200 m處出現(xiàn)團(tuán)塊狀高阻異常區(qū)，與鉆井揭示巖性不符，出現(xiàn)“假異?！?。

經(jīng)過(guò)“去噪”預(yù)處理及近場(chǎng)校正，斷面等值線形態(tài)更加合理，基本反映了真實(shí)的地質(zhì)信息，異常形態(tài)符合地質(zhì)規(guī)律，地電層劃分較明顯，其形態(tài)與已知地質(zhì)斷面基本一致。淺部為不連續(xù)的高阻特征，推測(cè)為青山群八畝地組凝灰質(zhì)砂礫巖電性反映，局部低阻異常推測(cè)為斷裂構(gòu)造發(fā)育部位。中部存在較連續(xù)水平展布的低阻異常帶，推測(cè)為淄博群坊子組含煤地層，低阻異常為采空區(qū)充水所致。深部則以高阻為主要特征，推測(cè)為新太古代泰山巖群片麻巖地層，局部低阻凹槽推測(cè)為斷裂構(gòu)造發(fā)育部位。后經(jīng)鉆探驗(yàn)證，煤炭采空區(qū)實(shí)際位置及埋深與物探解譯結(jié)果基本一致，L4線驗(yàn)證孔鉆遇采空區(qū)頂板深度與解譯深度誤差小于10%。

5 結(jié)論

多次的項(xiàng)目實(shí)踐表明，利用Matlab編程實(shí)現(xiàn)的CSAMT數(shù)據(jù)處理程序可以有效提高處理效率和效果，已經(jīng)成為GDP32多功能電法工作站CSAMT數(shù)據(jù)處理的有力工具。

1)隨著數(shù)據(jù)采集環(huán)境的變化，尤其當(dāng)前城市地質(zhì)勘查活動(dòng)越來(lái)越多，CSAMT方法遇到的“電磁干擾”也會(huì)越來(lái)越嚴(yán)重，甚至造成數(shù)據(jù)無(wú)法利用。所以，為了實(shí)現(xiàn)勘查目的和任務(wù)，既要完善野外采集環(huán)節(jié)的“抗噪”措施，也要通過(guò)室內(nèi)適當(dāng)?shù)摹叭ピ搿碧幚?，以提高?shù)據(jù)信噪比，克服“假異?！?，保證后續(xù)處理質(zhì)量。

2)近場(chǎng)校正可以有效提高CSAMT低頻段數(shù)據(jù)的利用程度，加大勘查深度，改善地電斷面形態(tài)，提高物探解譯的準(zhǔn)確性和精度。正演地電模型通過(guò)“全頻域視電阻率法”校正，得到了很好的效果。過(guò)渡區(qū)有“校正不足”現(xiàn)象，通過(guò)數(shù)據(jù)擬合得到的修正公式與頻率選擇有關(guān)系，在實(shí)際應(yīng)用的效果有待進(jìn)一步檢驗(yàn)。

3)Matlab作為強(qiáng)大的數(shù)值計(jì)算工具，可以在多個(gè)領(lǐng)域發(fā)揮其顯著的作用，其豐富的函數(shù)庫(kù)為數(shù)據(jù)處理功能的實(shí)現(xiàn)提供了有效支撐，并可以實(shí)現(xiàn)模塊擴(kuò)展和多語(yǔ)言聯(lián)合編程。