趙 玄,嚴(yán)家斌,2*,皇祥宇,2,胡 濤
(1.中南大學(xué) 地球科學(xué)與信息物理學(xué)院,湖南 長(zhǎng)沙 410083; 2.中南大學(xué) 有色資源與地質(zhì)災(zāi)害探查湖南省重點(diǎn)實(shí)驗(yàn)室,湖南 長(zhǎng)沙 410083)
大地電磁法(MT)誕生于20世紀(jì)50年代初,該方法具有工作效率高、成本低廉、勘探深度大、高阻層的屏蔽作用小等特點(diǎn)[1-2],廣泛應(yīng)用于礦產(chǎn)勘查、地下水與地?zé)峥碧健⒂蜌馄詹榈确矫?,產(chǎn)生了巨大的經(jīng)濟(jì)效益[3-5]。隨著城市的擴(kuò)展和人類工業(yè)活動(dòng)的加劇,大地電磁場(chǎng)信號(hào)的觀測(cè)越來(lái)越困難,干擾越來(lái)越嚴(yán)重,從受干擾的觀測(cè)數(shù)據(jù)中估計(jì)穩(wěn)定的阻抗或視電阻率是當(dāng)前地球物理電磁信號(hào)處理的難點(diǎn)之一。Sims等較早提出依據(jù)最小二乘法的張量阻抗估計(jì),該方法要求各數(shù)據(jù)道之間的噪聲相互獨(dú)立,但當(dāng)電磁資料存在相關(guān)噪聲以及觀測(cè)誤差不服從高斯分布規(guī)律時(shí),由最小二乘法估算的阻抗會(huì)產(chǎn)生嚴(yán)重偏差[6]。為了克服上述問(wèn)題,Egbert等將穩(wěn)健性估計(jì)(Robust Estimation)引入到地磁轉(zhuǎn)換函數(shù)中,并詳細(xì)介紹了加權(quán)最小二乘法和改進(jìn)算法[7]。為了進(jìn)一步進(jìn)行穩(wěn)健性估計(jì),Chave等詳細(xì)介紹了地球物理中幾個(gè)重要的物理量:功率譜、相干度和轉(zhuǎn)換函數(shù),并綜述了Robust估計(jì)理論,從統(tǒng)計(jì)學(xué)原理出發(fā)梳理了M估計(jì)和極大似然估計(jì)法[8]。針對(duì)海洋環(huán)境下的大地電磁法,柳建新等提出了基于相關(guān)歸一Robust法,根據(jù)相關(guān)系數(shù)的變化改進(jìn)Robust法的權(quán)系數(shù)[9]。張剛等將基于重復(fù)中位數(shù)估計(jì)的Robust法應(yīng)用于長(zhǎng)周期大地電磁阻抗張量估算中[10]。Chave等提出了基于Stable分布的最大似然阻抗估計(jì)法,根據(jù)數(shù)據(jù)自身的噪聲分布進(jìn)行加權(quán),能更合理地給出阻抗估計(jì)值[11-12]。常規(guī)的阻抗估計(jì)法對(duì)噪聲辨識(shí)度很低,未利用大地電磁阻抗的物理特性分辨高質(zhì)量數(shù)據(jù),導(dǎo)致噪聲嚴(yán)重干擾信號(hào)數(shù)據(jù),使得阻抗結(jié)果發(fā)生嚴(yán)重偏差。由于大地電磁數(shù)據(jù)中高質(zhì)量信號(hào)數(shù)據(jù)不能有效利用,所以聚類分析算法被引入電磁數(shù)據(jù)處理中。李晉等介紹了基于遞歸分析和聚類分析的大地電磁信噪辨識(shí)及分離方法[13],該方法雖然保留了大地電磁信號(hào)低頻段的緩慢變化信息,改善了低頻段大地電磁數(shù)據(jù)質(zhì)量,但對(duì)電場(chǎng)和磁場(chǎng)的處理是獨(dú)立進(jìn)行的,對(duì)電磁場(chǎng)的分量間相關(guān)性考慮不足。
本文從大地電磁阻抗的實(shí)虛分量特性出發(fā),通過(guò)定義阻抗的歐式距離構(gòu)建阻抗相似性度量,利用其相似性對(duì)受干擾阻抗、輕微受干擾阻抗、未受干擾阻抗進(jìn)行K中心點(diǎn)聚類分析,識(shí)別出高質(zhì)量的信號(hào),并基于仿真實(shí)驗(yàn)和實(shí)例分析驗(yàn)證該方法的有效性。
聚類分析(Clustering Analysis)是數(shù)據(jù)挖掘的主要方法之一[14-17],主要研究數(shù)據(jù)之間的相關(guān)性,這種相關(guān)性大小是通過(guò)一定的相似性準(zhǔn)則來(lái)判斷的,然后依此將對(duì)象分簇或聚類。在頻率域內(nèi)針對(duì)某一地質(zhì)構(gòu)造,高質(zhì)量的大地電磁信號(hào)數(shù)據(jù)的阻抗實(shí)虛分量必然在真實(shí)阻抗附近浮動(dòng),其阻抗間存在相似性,而與噪聲數(shù)據(jù)阻抗間存在差異性[18]。這一特性十分適合采用聚類分析進(jìn)行阻抗數(shù)據(jù)的優(yōu)選,篩選出高質(zhì)量的數(shù)據(jù)。采用聚類分析的結(jié)果使得信號(hào)數(shù)據(jù)盡可能歸為一類,信號(hào)數(shù)據(jù)和噪聲歸為不同類。
一個(gè)聚類分析過(guò)程的質(zhì)量取決于對(duì)度量標(biāo)準(zhǔn)的選擇。為了度量數(shù)據(jù)對(duì)象之間的相似或者接近程度,需要定義一些相似性度量標(biāo)準(zhǔn)。聚類分析中,為了表示兩個(gè)數(shù)據(jù)對(duì)象之間的相似度,一般采用特征空間中的距離作為度量標(biāo)準(zhǔn)來(lái)計(jì)算兩個(gè)樣本間的相異度[16-17]。為了描述阻抗數(shù)據(jù)之間的相似或者接近程度,提出了阻抗(Z)的歐式距離。其表達(dá)式為
(1)
由式(1)可知,歐式距離越小表示高質(zhì)量阻抗之間的相似或接近程度越大。相對(duì)于受噪聲干擾的數(shù)據(jù),受干擾較小的信號(hào)數(shù)據(jù)的阻抗歐式距離小。
K中心點(diǎn)(K-medoids)聚類分析是為了降低K-means算法[19]對(duì)噪聲孤立點(diǎn)的敏感度而提出的新算法。K中心點(diǎn)聚類分析選擇類或簇中最靠近均值中心點(diǎn)的一個(gè)對(duì)象來(lái)代表該類或簇的中心點(diǎn),可以有效消除孤立點(diǎn)對(duì)聚類效果的影響[20]。在聚類分析中只需要計(jì)算類或簇的中心點(diǎn),并選擇該類或簇內(nèi)靠近中心點(diǎn)的對(duì)象作為新中心點(diǎn),依次循環(huán)迭代,直到得到穩(wěn)定的類中心[21]。
以頻率域阻抗為例,對(duì)于包含N個(gè)阻抗數(shù)據(jù)對(duì)象的集合X,隨機(jī)選取K個(gè)數(shù)據(jù)對(duì)象作為初始聚類中心點(diǎn),并把K個(gè)初始聚類中心點(diǎn)表示成K個(gè)初始類。計(jì)算集合內(nèi)所有數(shù)據(jù)對(duì)象與各個(gè)初始聚類中心點(diǎn)的阻抗歐式距離,依據(jù)就近分配原則,將數(shù)據(jù)對(duì)象劃分到最靠近的類內(nèi),重新分配每個(gè)數(shù)據(jù)對(duì)象直至所有數(shù)據(jù)對(duì)象分配完畢,計(jì)算新的聚類中心點(diǎn)。K中心點(diǎn)聚類分析步驟為:①輸入K,X=[ReZ,ImZ];②選擇K個(gè)初始聚類中心點(diǎn),C[1]=X[1],C[2]=X[2],…,C[K]=X[K];③計(jì)算第i個(gè)類X[i]與第②步中K個(gè)初始聚類中心點(diǎn)的阻抗歐式距離,找到離它最近中心點(diǎn)C[m],將其分配到C[m]所在的類內(nèi),并對(duì)該點(diǎn)作標(biāo)記;④計(jì)算出m個(gè)類內(nèi)的所有數(shù)據(jù)對(duì)象均值,選擇靠近均值的數(shù)據(jù)對(duì)象作為m個(gè)類的新中心點(diǎn);⑤重復(fù)第③、④步,直至所有數(shù)據(jù)分配到類內(nèi)以及達(dá)到設(shè)定迭代次數(shù)。
經(jīng)過(guò)K中心點(diǎn)聚類分析后,每個(gè)頻點(diǎn)下的阻抗數(shù)據(jù)被劃分為K個(gè)類。把每個(gè)類的中心點(diǎn)作為該類阻抗的代表值,每個(gè)頻點(diǎn)下將會(huì)得到K個(gè)阻抗值和對(duì)應(yīng)的視電阻率。從K個(gè)結(jié)果中挑選出最佳阻抗值,依據(jù)電磁法理論定義類的選取準(zhǔn)則。
(1)相干度準(zhǔn)則。相干度可以衡量輸入信號(hào)和輸出信號(hào)的相關(guān)性,相干度越大,數(shù)據(jù)相關(guān)性越好,信號(hào)質(zhì)量越好。因此,設(shè)定相干度閾值,篩選出滿足閾值對(duì)應(yīng)的類。
(2)緊湊性準(zhǔn)則。為了描述數(shù)據(jù)聚集程度,提出了緊湊性概念。以某一頻點(diǎn)下的阻抗為例,對(duì)于包含N個(gè)阻抗數(shù)據(jù)對(duì)象的集合X={X1,X2,…,XN},X=[ReZ,ImZ],假設(shè)數(shù)據(jù)集合X被劃分為K個(gè)類,每個(gè)類中數(shù)據(jù)對(duì)象個(gè)數(shù)為{n1,n2,…,nk},且n1+n2+…+nk=N。以第j類為例,該類的緊湊性Cj(f)表達(dá)式為
緊湊性Cj(f)值越小,說(shuō)明類內(nèi)數(shù)據(jù)對(duì)象聚集的越緊湊,數(shù)據(jù)之間相似程度越高,更符合高質(zhì)量信號(hào)數(shù)據(jù)的阻抗高相似性特征,緊湊性好的類被認(rèn)為是高質(zhì)量信號(hào)數(shù)據(jù)所在的類。滿足相干度準(zhǔn)則的條件下,選擇緊湊性好的類[22]。
算法步驟為:①分別對(duì)M段時(shí)間域數(shù)據(jù)采用最小二乘法求取初始阻抗值,得到每個(gè)頻點(diǎn)下M個(gè)阻抗值;②輸入K,X=[ReZ,ImZ],對(duì)第①步中阻抗組成的集合X采用K中心點(diǎn)聚類分析,每個(gè)頻點(diǎn)下的阻抗被劃分到K個(gè)類內(nèi);③計(jì)算每個(gè)頻點(diǎn)下各類的相干度和緊湊性;④依據(jù)相干度準(zhǔn)則和緊湊性準(zhǔn)則,篩選出符合要求的類以及相對(duì)應(yīng)的阻抗值;⑤根據(jù)第④步中的阻抗值計(jì)算出相對(duì)應(yīng)的視電阻率。
大地電磁資料中存在著多種噪聲,使得數(shù)據(jù)處理異常困難,常規(guī)方法估算的結(jié)果會(huì)發(fā)生偏移[23-24]。通過(guò)在大地電磁仿真數(shù)據(jù)中加入噪聲,模擬實(shí)際情況下噪聲對(duì)信號(hào)的影響,并對(duì)比基于K中心點(diǎn)聚類分析和Robust法的估算結(jié)果。首先利用蒙特卡羅(Monte Carlo)法產(chǎn)生互不相關(guān)的隨機(jī)時(shí)間序列,作為電道或者磁道數(shù)據(jù),通過(guò)傅里葉變換把時(shí)間域序列轉(zhuǎn)化到頻率域,根據(jù)層狀介質(zhì)的視電阻率理論公式,加入設(shè)定的電性模型參數(shù),計(jì)算出磁道或者電道頻率域數(shù)據(jù),再通過(guò)傅里葉反變換轉(zhuǎn)換到時(shí)間域。本文設(shè)定的電性模型參數(shù)為:地下結(jié)構(gòu)為兩層均勻介質(zhì),第一層介質(zhì)電阻率為500 Ω·m,厚度為250 m,第二層介質(zhì)電阻率為100 Ω·m。上述仿真實(shí)驗(yàn)產(chǎn)生了x方向電道分量Ex和y方向電道分量Ey各50段,x方向磁道分量Hx和y方向磁道分量Hy各50段,各段數(shù)據(jù)采樣長(zhǎng)度為2 048,采樣頻率為10 000 Hz。在大地電磁仿真實(shí)驗(yàn)的電道或者磁道中加入類三角波噪聲和脈沖噪聲。
圖1 原始磁場(chǎng)和加噪聲后磁場(chǎng)的時(shí)間序列Fig.1 Time Series of Magnetic Field with and Without Noise
圖2 原始磁場(chǎng)和加噪聲后磁場(chǎng)的頻譜Fig.2 Frequency Spectra of Magnetic Field with and Without Noise
圖3 磁場(chǎng)噪聲阻抗聚類圖(f=24.41 Hz)Fig.3 Clustering Diagram of Impedance of Magnetic Field with Noise (f=24.41 Hz)
(1)隨機(jī)選出15段磁道仿真數(shù)據(jù),在選擇的數(shù)據(jù)段不同位置加入類三角波噪聲、脈沖噪聲。信號(hào)與類三角波噪聲幅值之比為1∶4,噪聲數(shù)據(jù)長(zhǎng)度占每段數(shù)據(jù)總長(zhǎng)度的10%;信號(hào)與脈沖噪聲幅值之比為1∶50。對(duì)50段數(shù)據(jù)分別采用Robust法和K中心點(diǎn)聚類分析估計(jì)視電阻率和相位。由圖1(a)和圖2(a)可知,不含噪聲的磁場(chǎng)頻譜能量是逐漸增大的。在磁場(chǎng)中加入類三角波噪聲和脈沖噪聲[圖1(b)],這兩種噪聲影響所有頻點(diǎn)的磁場(chǎng)頻譜[圖2(b)],使得頻譜能量明顯高于原始信號(hào)能量。對(duì)上述50段數(shù)據(jù)先進(jìn)行K中心點(diǎn)聚類分析,對(duì)分類后的阻抗按照相干度準(zhǔn)則和緊湊性準(zhǔn)則篩選出最佳阻抗值所在的類,計(jì)算最佳阻抗值所對(duì)應(yīng)的視電阻率和相位。以24.41 Hz頻點(diǎn)為例,阻抗數(shù)據(jù)被劃分到3個(gè)區(qū)域內(nèi)(圖3)。表1是該頻點(diǎn)下劃分成3個(gè)類對(duì)應(yīng)的參數(shù),由于第3類的相干度和緊湊性很好,所以選擇第3類為最佳類,并計(jì)算最佳類所對(duì)應(yīng)的視電阻率及相位。對(duì)其他頻點(diǎn)采用相同計(jì)算方法,得到視電阻率曲線(圖4)和相位曲線(圖5)。從圖4、5可以看到,基于K中心點(diǎn)聚類分析估算的結(jié)果更接近理論曲線,當(dāng)頻率高于1 000 Hz時(shí),基于Robust法估算的視電阻率小于理論值。相對(duì)于Robust法,基于K中心點(diǎn)聚類分析可以分離和識(shí)別出高質(zhì)量數(shù)據(jù)對(duì)應(yīng)的阻抗,削弱了噪聲的干擾,使結(jié)果更為可靠。
表1 基于K中心點(diǎn)聚類分析的磁場(chǎng)阻抗識(shí)別評(píng)價(jià)參數(shù)(f=24.41 Hz)Tab.1 Evaluation Parameters of Impedance Recognition for Magnetic Field Based on K-medoids Clustering Analysis (f=24.41 Hz)
注:理論上,視電阻率為194.3 Ω·m,相位為31.0°。
圖4 磁場(chǎng)噪聲視電阻率曲線Fig.4 Apparent Resistivity Curves of Magnetic Field with Noise
圖5 磁場(chǎng)噪聲相位曲線Fig.5 Phase Curves of Magnetic Field with Noise
圖6 原始電場(chǎng)和加噪聲后電場(chǎng)的時(shí)間序列Fig.6 Time Series of Electric Field with and Without Noise
(2)隨機(jī)選出10段電道仿真數(shù)據(jù),在選擇數(shù)據(jù)段的不同位置加入類三角波噪聲、脈沖噪聲。信號(hào)與類三角波噪聲幅值之比為1∶4,噪聲數(shù)據(jù)長(zhǎng)度占每段數(shù)據(jù)總長(zhǎng)度的10%;信號(hào)與脈沖噪聲幅值之比為1∶50。對(duì)50段數(shù)據(jù)分別采用Robust法和K中心點(diǎn)聚類分析估計(jì)視電阻率和相位。由圖6(a)和圖7(a)可知,不含噪聲的電場(chǎng)頻譜能量是逐漸增大的。在電場(chǎng)中加入類三角波噪聲和脈沖噪聲[圖6(b)],這兩種噪聲影響所有頻點(diǎn)的電場(chǎng)頻譜[圖7(b)],使得頻譜能量明顯高于原始信號(hào)能量。對(duì)上述50段數(shù)據(jù)進(jìn)行K中心點(diǎn)聚類分析,對(duì)分類后的阻抗按照相干度準(zhǔn)則和緊湊性準(zhǔn)則,篩選出最佳阻抗值所在的類,并計(jì)算其對(duì)應(yīng)的視電阻率和相位。以43.95 Hz頻點(diǎn)為例,阻抗數(shù)據(jù)被劃分到4個(gè)區(qū)域內(nèi)(圖8)。表2是該頻點(diǎn)下劃分成4個(gè)類分別對(duì)應(yīng)的參數(shù),第2類的相干度和緊湊性很好,選擇第2類作為最佳類,并計(jì)算其對(duì)應(yīng)視電阻率和相位。對(duì)其他頻點(diǎn)采用相同計(jì)算方法,得到視電阻率曲線(圖9)和相位曲線(圖10),基于K中心點(diǎn)聚類分析的估算結(jié)果更接近理論值,而基于Robust法估算的視電阻率曲線和相位曲線都出現(xiàn)了波動(dòng),結(jié)果偏離了理論值,由此說(shuō)明基于K中心點(diǎn)聚類分析可以識(shí)別和篩選出高質(zhì)量信號(hào)數(shù)據(jù)的阻抗,進(jìn)而估算出更為可靠的結(jié)果。
圖7 原始電場(chǎng)和加噪聲后電場(chǎng)的頻譜Fig.7 Frequency Spectra of Electric Field with and Without Noise
圖8 電場(chǎng)噪聲阻抗聚類圖(f=43.95 Hz)Fig.8 Clustering Diagram of Impedance of Electric Field with Noise (f=43.95 Hz)
類序號(hào)相干度緊湊性視電阻率/(Ω·m)相位/(°)點(diǎn)數(shù)10.081.60440 000.067.9420.990.02186.130.93930.092.70480 000.076.5540.168.408 100 000.0-32.02
注:理論上,視電阻率為186.1 Ω·m,相位為30.7°。
圖9 電場(chǎng)噪聲視電阻率曲線Fig.9 Apparent Resistivity Curves of Electric Field with Noise
圖10 電場(chǎng)噪聲相位曲線Fig.10 Phase Curves of Electric Field with Noise
圖11 實(shí)測(cè)磁場(chǎng)和實(shí)測(cè)電場(chǎng)時(shí)間序列Fig.11 Time Series of Measured Magnetic Field and Electric Field
為驗(yàn)證基于K中心點(diǎn)聚類分析的實(shí)際應(yīng)用效果,選取云南省昭通市牛欄江天花板水電站田壩村堆積體處電磁成像系統(tǒng)觀測(cè)數(shù)據(jù)進(jìn)行驗(yàn)證。觀測(cè)區(qū)周圍無(wú)高壓線,但居民較多。圖11是采集的原始時(shí)間序列數(shù)據(jù),磁場(chǎng)在水平軸260~280 ms之間存在似類三角波噪聲,電場(chǎng)在水平軸0~20 ms之間存在明顯的似類三角波噪聲和脈沖噪聲。該測(cè)點(diǎn)電道和磁道信號(hào)中存在明顯的類三角波噪聲和脈沖噪聲(圖11),同時(shí)還存在幅度較小的類方波等噪聲。由仿真實(shí)驗(yàn)可知,不含噪聲的電磁數(shù)據(jù)頻譜能量逐漸增大,而測(cè)點(diǎn)的頻譜能量逐漸減小,說(shuō)明噪聲影響了所有頻點(diǎn)的頻譜(圖12)。對(duì)測(cè)點(diǎn)數(shù)據(jù)采用和仿真實(shí)驗(yàn)一樣的K中心點(diǎn)聚類分析,對(duì)阻抗數(shù)據(jù)進(jìn)行聚類分析。以頻點(diǎn)82.5 Hz為例,該頻點(diǎn)下阻抗經(jīng)過(guò)K中心點(diǎn)聚類分析,劃分成6個(gè)區(qū)域(圖13),每個(gè)區(qū)域代表一個(gè)類。由圖13可知,第2、3、6類內(nèi)數(shù)據(jù)點(diǎn)比較分散,說(shuō)明類內(nèi)的阻抗受噪聲影響很大,而第1、4、5類內(nèi)的數(shù)據(jù)點(diǎn)分布集中,說(shuō)明高質(zhì)量數(shù)據(jù)比較多。由表3可知,第1、4、5類所對(duì)應(yīng)的相干度和緊湊性很好,而第5類所對(duì)應(yīng)的緊湊性最好,因此,選擇第5類為最佳類,把其對(duì)應(yīng)的視電阻率作為最終的視電阻率。對(duì)其他頻點(diǎn)采用相同計(jì)算步驟,得到視電阻率曲線(圖14)。由圖14可知,相對(duì)于基于Robust法的估算結(jié)果,基于K中心點(diǎn)聚類分析的估算結(jié)果更為光滑連續(xù),不會(huì)出現(xiàn)異常跳動(dòng),更貼近真實(shí)地下結(jié)構(gòu),說(shuō)明該方法可以識(shí)別出高質(zhì)量數(shù)據(jù)的阻抗,降低了噪聲干擾,使結(jié)果更為可信。
表3 基于K中心點(diǎn)聚類分析的實(shí)測(cè)阻抗識(shí)別評(píng)價(jià)參數(shù)(f=82.5 Hz)Tab.3 Evaluation Parameters of Measured Impedance Recognition Based on K-medoids Clustering Analysis (f=82.5 Hz)
圖12 實(shí)測(cè)磁場(chǎng)和電場(chǎng)頻譜Fig.12 Frequency Spectra of Measured Magnetic Field and Electric Field
圖13 實(shí)測(cè)數(shù)據(jù)阻抗聚類圖(f=82.5 Hz)Fig.13 Clustering Diagram of Impedance of Measured Data (f=82.5 Hz)
圖14 實(shí)測(cè)數(shù)據(jù)視電阻率曲線Fig.14 Apparent Resistivity Curves of Measured Data
(1)從大地電磁阻抗的實(shí)虛分量特性出發(fā),定義了阻抗歐氏距離,描述阻抗數(shù)據(jù)之間的相似性,提出了基于K中心點(diǎn)聚類分析的電磁信號(hào)識(shí)別及阻抗提取方法,并依據(jù)類的選取準(zhǔn)則,識(shí)別出可靠的阻抗值。通過(guò)仿真實(shí)驗(yàn)和實(shí)例分析,基于K中心點(diǎn)聚類分析可以識(shí)別并提取出高質(zhì)量數(shù)據(jù)的阻抗,得到穩(wěn)定阻抗值和視電阻率,提高估算的可靠性和穩(wěn)定性。
(2)Robust法更適用于磁場(chǎng)受干擾小時(shí)的情況,當(dāng)電場(chǎng)或者磁場(chǎng)含有異常幅度大的噪聲時(shí),Robust法估算的結(jié)果不再穩(wěn)定,會(huì)使結(jié)果發(fā)生偏差;而K中心點(diǎn)聚類分析是利用阻抗的物理特性,高質(zhì)量數(shù)據(jù)的阻抗會(huì)集中分布且具有高度的相似性,因此,受影響較小,得到的阻抗更合理有效。由于個(gè)別異常點(diǎn)的存在和初始聚類數(shù)目的選擇都會(huì)影響聚類算法效率,所以本文的方法還需要進(jìn)一步的優(yōu)化。