基于K中心點聚類分析的大地電磁阻抗識別

趙 玄，嚴家斌,2*，皇祥宇,2，胡 濤

(1.中南大學 地球科學與信息物理學院，湖南 長沙 410083； 2.中南大學 有色資源與地質災害探查湖南省重點實驗室，湖南 長沙 410083)

0 引 言

大地電磁法(MT)誕生于20世紀50年代初，該方法具有工作效率高、成本低廉、勘探深度大、高阻層的屏蔽作用小等特點[1-2]，廣泛應用于礦產勘查、地下水與地熱勘探、油氣普查等方面，產生了巨大的經濟效益[3-5]。隨著城市的擴展和人類工業(yè)活動的加劇，大地電磁場信號的觀測越來越困難，干擾越來越嚴重，從受干擾的觀測數據中估計穩(wěn)定的阻抗或視電阻率是當前地球物理電磁信號處理的難點之一。Sims等較早提出依據最小二乘法的張量阻抗估計，該方法要求各數據道之間的噪聲相互獨立，但當電磁資料存在相關噪聲以及觀測誤差不服從高斯分布規(guī)律時，由最小二乘法估算的阻抗會產生嚴重偏差[6]。為了克服上述問題，Egbert等將穩(wěn)健性估計(Robust Estimation)引入到地磁轉換函數中，并詳細介紹了加權最小二乘法和改進算法[7]。為了進一步進行穩(wěn)健性估計，Chave等詳細介紹了地球物理中幾個重要的物理量：功率譜、相干度和轉換函數，并綜述了Robust估計理論，從統(tǒng)計學原理出發(fā)梳理了M估計和極大似然估計法[8]。針對海洋環(huán)境下的大地電磁法，柳建新等提出了基于相關歸一Robust法，根據相關系數的變化改進Robust法的權系數[9]。張剛等將基于重復中位數估計的Robust法應用于長周期大地電磁阻抗張量估算中[10]。Chave等提出了基于Stable分布的最大似然阻抗估計法，根據數據自身的噪聲分布進行加權，能更合理地給出阻抗估計值[11-12]。常規(guī)的阻抗估計法對噪聲辨識度很低，未利用大地電磁阻抗的物理特性分辨高質量數據，導致噪聲嚴重干擾信號數據，使得阻抗結果發(fā)生嚴重偏差。由于大地電磁數據中高質量信號數據不能有效利用，所以聚類分析算法被引入電磁數據處理中。李晉等介紹了基于遞歸分析和聚類分析的大地電磁信噪辨識及分離方法[13]，該方法雖然保留了大地電磁信號低頻段的緩慢變化信息，改善了低頻段大地電磁數據質量，但對電場和磁場的處理是獨立進行的，對電磁場的分量間相關性考慮不足。

本文從大地電磁阻抗的實虛分量特性出發(fā)，通過定義阻抗的歐式距離構建阻抗相似性度量，利用其相似性對受干擾阻抗、輕微受干擾阻抗、未受干擾阻抗進行K中心點聚類分析，識別出高質量的信號，并基于仿真實驗和實例分析驗證該方法的有效性。

1 原 理

聚類分析(Clustering Analysis)是數據挖掘的主要方法之一[14-17]，主要研究數據之間的相關性，這種相關性大小是通過一定的相似性準則來判斷的，然后依此將對象分簇或聚類。在頻率域內針對某一地質構造，高質量的大地電磁信號數據的阻抗實虛分量必然在真實阻抗附近浮動，其阻抗間存在相似性，而與噪聲數據阻抗間存在差異性[18]。這一特性十分適合采用聚類分析進行阻抗數據的優(yōu)選，篩選出高質量的數據。采用聚類分析的結果使得信號數據盡可能歸為一類，信號數據和噪聲歸為不同類。

1.1 相似性度量

一個聚類分析過程的質量取決于對度量標準的選擇。為了度量數據對象之間的相似或者接近程度，需要定義一些相似性度量標準。聚類分析中，為了表示兩個數據對象之間的相似度，一般采用特征空間中的距離作為度量標準來計算兩個樣本間的相異度[16-17]。為了描述阻抗數據之間的相似或者接近程度，提出了阻抗(Z)的歐式距離。其表達式為

(1)

由式(1)可知，歐式距離越小表示高質量阻抗之間的相似或接近程度越大。相對于受噪聲干擾的數據，受干擾較小的信號數據的阻抗歐式距離小。

1.2 K中心點聚類分析

K中心點(K-medoids)聚類分析是為了降低K-means算法[19]對噪聲孤立點的敏感度而提出的新算法。K中心點聚類分析選擇類或簇中最靠近均值中心點的一個對象來代表該類或簇的中心點，可以有效消除孤立點對聚類效果的影響[20]。在聚類分析中只需要計算類或簇的中心點，并選擇該類或簇內靠近中心點的對象作為新中心點，依次循環(huán)迭代，直到得到穩(wěn)定的類中心[21]。

以頻率域阻抗為例，對于包含N個阻抗數據對象的集合X，隨機選取K個數據對象作為初始聚類中心點，并把K個初始聚類中心點表示成K個初始類。計算集合內所有數據對象與各個初始聚類中心點的阻抗歐式距離，依據就近分配原則，將數據對象劃分到最靠近的類內，重新分配每個數據對象直至所有數據對象分配完畢，計算新的聚類中心點。K中心點聚類分析步驟為：①輸入K，X=[ReZ,ImZ]；②選擇K個初始聚類中心點，C[1]=X[1]，C[2]=X[2]，…，C[K]=X[K]；③計算第i個類X[i]與第②步中K個初始聚類中心點的阻抗歐式距離，找到離它最近中心點C[m]，將其分配到C[m]所在的類內，并對該點作標記；④計算出m個類內的所有數據對象均值，選擇靠近均值的數據對象作為m個類的新中心點；⑤重復第③、④步，直至所有數據分配到類內以及達到設定迭代次數。

1.3 類的選取準則

經過K中心點聚類分析后，每個頻點下的阻抗數據被劃分為K個類。把每個類的中心點作為該類阻抗的代表值，每個頻點下將會得到K個阻抗值和對應的視電阻率。從K個結果中挑選出最佳阻抗值，依據電磁法理論定義類的選取準則。

(1)相干度準則。相干度可以衡量輸入信號和輸出信號的相關性，相干度越大，數據相關性越好，信號質量越好。因此，設定相干度閾值，篩選出滿足閾值對應的類。

(2)緊湊性準則。為了描述數據聚集程度，提出了緊湊性概念。以某一頻點下的阻抗為例，對于包含N個阻抗數據對象的集合X={X1,X2,…,XN}，X=[ReZ,ImZ]，假設數據集合X被劃分為K個類，每個類中數據對象個數為{n1,n2,…,nk}，且n1+n2+…+nk=N。以第j類為例，該類的緊湊性Cj(f)表達式為

緊湊性Cj(f)值越小，說明類內數據對象聚集的越緊湊，數據之間相似程度越高，更符合高質量信號數據的阻抗高相似性特征，緊湊性好的類被認為是高質量信號數據所在的類。滿足相干度準則的條件下，選擇緊湊性好的類[22]。

1.4 算法步驟

算法步驟為：①分別對M段時間域數據采用最小二乘法求取初始阻抗值，得到每個頻點下M個阻抗值；②輸入K，X=[ReZ,ImZ]，對第①步中阻抗組成的集合X采用K中心點聚類分析，每個頻點下的阻抗被劃分到K個類內；③計算每個頻點下各類的相干度和緊湊性；④依據相干度準則和緊湊性準則，篩選出符合要求的類以及相對應的阻抗值；⑤根據第④步中的阻抗值計算出相對應的視電阻率。

2 仿真實驗

大地電磁資料中存在著多種噪聲，使得數據處理異常困難，常規(guī)方法估算的結果會發(fā)生偏移[23-24]。通過在大地電磁仿真數據中加入噪聲，模擬實際情況下噪聲對信號的影響，并對比基于K中心點聚類分析和Robust法的估算結果。首先利用蒙特卡羅(Monte Carlo)法產生互不相關的隨機時間序列，作為電道或者磁道數據，通過傅里葉變換把時間域序列轉化到頻率域，根據層狀介質的視電阻率理論公式，加入設定的電性模型參數，計算出磁道或者電道頻率域數據，再通過傅里葉反變換轉換到時間域。本文設定的電性模型參數為：地下結構為兩層均勻介質，第一層介質電阻率為500 Ω·m，厚度為250 m，第二層介質電阻率為100 Ω·m。上述仿真實驗產生了x方向電道分量Ex和y方向電道分量Ey各50段，x方向磁道分量Hx和y方向磁道分量Hy各50段，各段數據采樣長度為2 048，采樣頻率為10 000 Hz。在大地電磁仿真實驗的電道或者磁道中加入類三角波噪聲和脈沖噪聲。

圖1 原始磁場和加噪聲后磁場的時間序列Fig.1 Time Series of Magnetic Field with and Without Noise

圖2 原始磁場和加噪聲后磁場的頻譜Fig.2 Frequency Spectra of Magnetic Field with and Without Noise

圖3 磁場噪聲阻抗聚類圖(f=24.41 Hz)Fig.3 Clustering Diagram of Impedance of Magnetic Field with Noise (f=24.41 Hz)

(1)隨機選出15段磁道仿真數據，在選擇的數據段不同位置加入類三角波噪聲、脈沖噪聲。信號與類三角波噪聲幅值之比為1∶4，噪聲數據長度占每段數據總長度的10%；信號與脈沖噪聲幅值之比為1∶50。對50段數據分別采用Robust法和K中心點聚類分析估計視電阻率和相位。由圖1(a)和圖2(a)可知，不含噪聲的磁場頻譜能量是逐漸增大的。在磁場中加入類三角波噪聲和脈沖噪聲[圖1(b)]，這兩種噪聲影響所有頻點的磁場頻譜[圖2(b)]，使得頻譜能量明顯高于原始信號能量。對上述50段數據先進行K中心點聚類分析，對分類后的阻抗按照相干度準則和緊湊性準則篩選出最佳阻抗值所在的類，計算最佳阻抗值所對應的視電阻率和相位。以24.41 Hz頻點為例，阻抗數據被劃分到3個區(qū)域內(圖3)。表1是該頻點下劃分成3個類對應的參數，由于第3類的相干度和緊湊性很好，所以選擇第3類為最佳類，并計算最佳類所對應的視電阻率及相位。對其他頻點采用相同計算方法，得到視電阻率曲線(圖4)和相位曲線(圖5)。從圖4、5可以看到，基于K中心點聚類分析估算的結果更接近理論曲線，當頻率高于1 000 Hz時，基于Robust法估算的視電阻率小于理論值。相對于Robust法，基于K中心點聚類分析可以分離和識別出高質量數據對應的阻抗，削弱了噪聲的干擾，使結果更為可靠。

表1 基于K中心點聚類分析的磁場阻抗識別評價參數(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 磁場噪聲視電阻率曲線Fig.4 Apparent Resistivity Curves of Magnetic Field with Noise

圖5 磁場噪聲相位曲線Fig.5 Phase Curves of Magnetic Field with Noise

圖6 原始電場和加噪聲后電場的時間序列Fig.6 Time Series of Electric Field with and Without Noise

(2)隨機選出10段電道仿真數據，在選擇數據段的不同位置加入類三角波噪聲、脈沖噪聲。信號與類三角波噪聲幅值之比為1∶4，噪聲數據長度占每段數據總長度的10%；信號與脈沖噪聲幅值之比為1∶50。對50段數據分別采用Robust法和K中心點聚類分析估計視電阻率和相位。由圖6(a)和圖7(a)可知，不含噪聲的電場頻譜能量是逐漸增大的。在電場中加入類三角波噪聲和脈沖噪聲[圖6(b)]，這兩種噪聲影響所有頻點的電場頻譜[圖7(b)]，使得頻譜能量明顯高于原始信號能量。對上述50段數據進行K中心點聚類分析，對分類后的阻抗按照相干度準則和緊湊性準則，篩選出最佳阻抗值所在的類，并計算其對應的視電阻率和相位。以43.95 Hz頻點為例，阻抗數據被劃分到4個區(qū)域內(圖8)。表2是該頻點下劃分成4個類分別對應的參數，第2類的相干度和緊湊性很好，選擇第2類作為最佳類，并計算其對應視電阻率和相位。對其他頻點采用相同計算方法，得到視電阻率曲線(圖9)和相位曲線(圖10)，基于K中心點聚類分析的估算結果更接近理論值，而基于Robust法估算的視電阻率曲線和相位曲線都出現了波動，結果偏離了理論值，由此說明基于K中心點聚類分析可以識別和篩選出高質量信號數據的阻抗，進而估算出更為可靠的結果。

圖7 原始電場和加噪聲后電場的頻譜Fig.7 Frequency Spectra of Electric Field with and Without Noise

圖8 電場噪聲阻抗聚類圖(f=43.95 Hz)Fig.8 Clustering Diagram of Impedance of Electric Field with Noise (f=43.95 Hz)

類序號相干度緊湊性視電阻率/(Ω·m)相位/(°)點數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 電場噪聲視電阻率曲線Fig.9 Apparent Resistivity Curves of Electric Field with Noise

圖10 電場噪聲相位曲線Fig.10 Phase Curves of Electric Field with Noise

3 實例分析

圖11 實測磁場和實測電場時間序列Fig.11 Time Series of Measured Magnetic Field and Electric Field

為驗證基于K中心點聚類分析的實際應用效果，選取云南省昭通市牛欄江天花板水電站田壩村堆積體處電磁成像系統(tǒng)觀測數據進行驗證。觀測區(qū)周圍無高壓線，但居民較多。圖11是采集的原始時間序列數據，磁場在水平軸260～280 ms之間存在似類三角波噪聲，電場在水平軸0～20 ms之間存在明顯的似類三角波噪聲和脈沖噪聲。該測點電道和磁道信號中存在明顯的類三角波噪聲和脈沖噪聲(圖11)，同時還存在幅度較小的類方波等噪聲。由仿真實驗可知，不含噪聲的電磁數據頻譜能量逐漸增大，而測點的頻譜能量逐漸減小，說明噪聲影響了所有頻點的頻譜(圖12)。對測點數據采用和仿真實驗一樣的K中心點聚類分析，對阻抗數據進行聚類分析。以頻點82.5 Hz為例，該頻點下阻抗經過K中心點聚類分析，劃分成6個區(qū)域(圖13)，每個區(qū)域代表一個類。由圖13可知，第2、3、6類內數據點比較分散，說明類內的阻抗受噪聲影響很大，而第1、4、5類內的數據點分布集中，說明高質量數據比較多。由表3可知，第1、4、5類所對應的相干度和緊湊性很好，而第5類所對應的緊湊性最好，因此，選擇第5類為最佳類，把其對應的視電阻率作為最終的視電阻率。對其他頻點采用相同計算步驟，得到視電阻率曲線(圖14)。由圖14可知，相對于基于Robust法的估算結果，基于K中心點聚類分析的估算結果更為光滑連續(xù)，不會出現異常跳動，更貼近真實地下結構，說明該方法可以識別出高質量數據的阻抗，降低了噪聲干擾，使結果更為可信。

表3 基于K中心點聚類分析的實測阻抗識別評價參數(f=82.5 Hz)Tab.3 Evaluation Parameters of Measured Impedance Recognition Based on K-medoids Clustering Analysis (f=82.5 Hz)

圖12 實測磁場和電場頻譜Fig.12 Frequency Spectra of Measured Magnetic Field and Electric Field

圖13 實測數據阻抗聚類圖(f=82.5 Hz)Fig.13 Clustering Diagram of Impedance of Measured Data (f=82.5 Hz)

圖14 實測數據視電阻率曲線Fig.14 Apparent Resistivity Curves of Measured Data

4 結 語

(1)從大地電磁阻抗的實虛分量特性出發(fā)，定義了阻抗歐氏距離，描述阻抗數據之間的相似性，提出了基于K中心點聚類分析的電磁信號識別及阻抗提取方法，并依據類的選取準則，識別出可靠的阻抗值。通過仿真實驗和實例分析，基于K中心點聚類分析可以識別并提取出高質量數據的阻抗，得到穩(wěn)定阻抗值和視電阻率，提高估算的可靠性和穩(wěn)定性。

(2)Robust法更適用于磁場受干擾小時的情況，當電場或者磁場含有異常幅度大的噪聲時，Robust法估算的結果不再穩(wěn)定，會使結果發(fā)生偏差；而K中心點聚類分析是利用阻抗的物理特性，高質量數據的阻抗會集中分布且具有高度的相似性，因此，受影響較小，得到的阻抗更合理有效。由于個別異常點的存在和初始聚類數目的選擇都會影響聚類算法效率，所以本文的方法還需要進一步的優(yōu)化。