基于優化聚類的人工源電磁法數據信噪分離方法

胡艷芳, 劉子杰, 李帝銓*, 張賢, 索光運

1 湖南工商大學微電子與物理學院, 長沙 410205 2 有色金屬成礦預測與地質環境監測教育部重點實驗室(中南大學), 長沙 410083 3 中南大學地球科學與信息物理學院, 長沙 410083

0 引言

電磁法是地球物理勘探的重要分支,根據場源不同,主要分為天然源和人工源電磁法(Nabighian,1988;何繼善,1990, 2010).人工源電磁法(CSEM)是一類利用可控源作為激勵信號的電磁勘探方法,其中包括了可控源音頻大地電磁法(Controlled Source Audio-frequency Magnetotellurics, CSAMT)、廣域電磁法(Wide Field Electromagnetic Method, WFEM)、瞬變電磁法(Transient Electromagnetic Method, TEM)以及時頻電磁法(Time-frequency Electromagnetic Method, TFEM)等,該類方法采用人工源取代天然場源作為激勵源,克服了天然場源非平穩性、信號微弱以及極化方向隨機的缺點.相比于天然場其能量有很大的提升,在一定程度上提高了方法的抗干擾能力.然而隨著人類活動范圍的擴大,噪聲的影響也在不斷加強,尤其在強干擾的礦集區以及人類活動頻繁的城區及周邊,各種大尺度的非周期噪聲和工頻及其諧波干擾,使得電磁法勘探無法獲得高信噪比的數據(湯井田等,2012, 2018; Li G et al., 2017; 李晉等, 2018, 2019),因此針對CSEM信號的去噪研究也越來越多.雷達等(2010)聯合數據信息熵和有理函數濾波的方法對CSAMT數據進行了處理,有效地提高了信噪比,改善了實測數據質量;Streich等(2013)、Imamura等(2018)利用Robust估計的方法對可控源電磁數據進行了處理;張必明等(2015)提出了一種自適應雙向均方差閾值法的有源電磁法信噪分離技術,在少量頻點信息被噪聲淹沒時,能有效剔除頻譜數據中的粗大誤差,但僅采用單一的閾值控制可能造成有效數據的過處理及部分噪聲不明顯的數據未被處理的情況,同時也沒有考慮頻點之間的連續性;Mo等(2017)提出了一種基于灰色建模系統結合閾值法對異常進行剔除,然后利用穩健M估計對剩余樣本數據進行估計,得到對應頻率測量值;陳超健等(2019)提出一種分步處理方法,首先基于灰色判別理論對數據中明顯的離群點進行處理,恢復曲線的基本形態,再通過有理函數濾波剔除數據中的殘余噪聲;Yang等(2018)提出了一種基于最小二乘反演的CSEM信號處理方法,利用不含噪或含隨機噪聲的時間序列片段,求解噪聲頻譜,獲得了高質量數據;Yang Z 等(2021)提出一種基于消除法、相關辨識與維納濾波反褶積求解大地脈沖響應的可控源電磁法周期性工頻干擾消除方法,其結果能有效消除電力線干擾;楊洋等(2022a)通過傅里葉級數構造超定方程,基于最小二乘反演方法,有效識別并剔除了半航空瞬變電磁法(SA-TEM)數據中的運動噪聲.

除此之外,隨著機器學習和人工智能方法的迅速發展,不少學者將其應用到信號處理方面.如基于人工神經網絡的時間序列編輯,識別并刪除具有噪聲特征的序列片段以及對缺失或不連續時間序列進行預測等(韓盈等, 2021).由于人文干擾通常具有持續性、規律性的特點,Li等(2021)提出了一種基于字典學習和移不變稀疏編碼與互補集合經驗模態分解結合的CSEM去噪方法,該方法利用周期信號可以進行稀疏表示的特性,引入移不變稀疏編碼(Shift Invariant Sparse Coding, SISC)字典學習方法進行稀疏分解去噪;Xue 等(2020)利用字典學習的方法對時間域航空電磁數據中的運動噪聲進行了處理,取得了較好的應用效果;許滔滔等(2020)提出了一種基于機器學習的工頻干擾去除方法,該方法首先建立長短時記憶的循環神經網絡模型,然后構建數據集訓練模型,最后將含噪數據輸入到訓練好的模型中,提取工頻干擾信號,從而獲取有效信號;李晉等(2017)采用K-means和模糊C均值聚類的方法進行信噪辨識,并結合數學形態濾波對強干擾的時間序列進行噪聲壓制,取得了比較好的效果;Zhang等(2021)根據信號時間域特征,提出了一種基于精細復合多尺度散布熵結合正交匹配追蹤的模糊C均值聚類方法,通過聚類識別噪聲,再將其進行去除,對于在時間域特征比較明顯的噪聲信號,該方法識別噪聲的效果明顯.綜上所述,不論是時間域還是頻率域數據處理方法,應用條件都具有局限性,對單一類型的噪聲或者特征明顯的噪聲,處理效果比較好,對于復雜噪聲,處理效果欠佳.

基于上述分析,在實際應用中,噪聲的復雜性決定了單一的去噪方法不可能對所有類型的干擾噪聲都能取得很好的效果.因此,本文針對電磁噪聲的特點,提出了一種基于加權自適應帶寬估計的均值漂移聚類方法,對CSEM數據在復平面內進行信噪分離.首先在傳統的MSC算法基礎上,引入權重函數,降低MSC處理結果對帶寬選擇的敏感度,提高算法穩健性;其次針對CSEM實測數據分布特征提出了一種基于局部密度梯度的帶寬估計方法,使帶寬選取自適應化;最后通過仿真數據和實測廣域電磁法數據對本文方法進行驗證,結果表明本文方法能有效降低強電磁干擾影響,改善數據質量.

1 CSEM信號與噪聲特點

CSEM信號大多是基于偽隨機編碼發送的,其最大的特點是具有周期性,頻率已知,如周期方波信號就是最簡單的偽隨機信號.在實際應用中,通常假設大地是一個線性時不變系統,CSEM有效信號經過大地濾波之后,其振幅和相位會發生變化.在未受到電磁干擾時,這種變化不會因為采集時間的不同而發生改變,因此在復平面內,同一測點同一頻率不同時間對應的值是固定的,多周期采集的數據會集中在一點;當受到非周期噪聲干擾時,其振幅和相位會發生改變,不同周期采集的數據由于受到的干擾類型和強度不同,振幅和相位的變化不盡相同,因此會以真值為中心出現不同程度的偏離,偏離程度與受噪聲影響的程度相關,影響程度越大,與真值的距離越遠;反之,影響程度越小,與真值的距離越近.在實際應用中,最理想的狀態是只受高斯噪聲影響,此時呈現出高斯分布特征.圖1展示了在2n序列偽隨機信號(7-2頻組)(何繼善, 2010; Yang Y et al.,2021; Zhang et al.,2022)中加入幾種典型噪聲之后的時間域波形與相應頻譜,有效信號的頻率為64 Hz、32 Hz、16 Hz、8 Hz、4 Hz、2 Hz、1 Hz.受噪聲影響,不同有效頻率受到影響的尺度不同,且影響范圍向高頻擴展.圖2為不同頻率信號受到典型噪聲影響時在復平面內的分布特征,噪聲和信號在時間域的疊加,在復平面內表現出來的就是矢量相加,矢量的相加改變的不只是振幅,還有相位.噪聲的幅度、尺度以及引入時間不同,對有效頻率振幅和相位的影響也不同,而不同類型噪聲疊加,頻譜會變得更加復雜.在實際應用中,噪聲的影響是復雜多樣的,除周期噪聲,長時段實測非周期噪聲的影響多數可近似為服從高斯分布.圖3為實測不同時長噪聲在有效頻率位置的分布特征,由圖可知,噪聲在復平面內呈現近似高斯分布的特征,且采用不同時長對信號進行分段,噪聲影響范圍會縮小,說明采用冗余觀測能有效壓制高斯噪聲的影響.

圖1 不同噪聲類型對有效信號的影響規律

圖2 復平面內不同噪聲類型對有效信號的影響

圖3 實測噪聲監測點不同時長數據在復平面內的分布特征

2 算法原理

2.1 傳統MSC算法原理

均值漂移(Mean-shift,簡稱MS)的概念在Fukunaga和Hostetler(1975)的論文中出現之后,許多學者對該算法進行了改進并將其應用在目標跟蹤、圖像平滑及分割等領域.MS算法最初被用來進行密度估計,而均值漂移聚類(Mean-shift Clustering,簡稱MSC)實際上也是一種基于概率密度的聚類方法,通過計算偏移向量從低密度區域向高密度區域移動的一種滑動算法.

給定一個n維的數據空間S={s1,s2,…,sn},對于空間內的任意樣本點sj,以樣本sj為圓心,帶寬bandwidth為半徑形成的圓形區域稱為樣本點sj的局部空間,則樣本點sj對應的偏移向量為以該樣本點為起點到其局部空間內所有樣本點的向量的均值.即

(1)

式中M(sj)表示樣本點sj對應的局部空間的偏移向量,m為sj局部空間內的樣本個數.MSC算法的核心就是計算偏移向量,局部樣本空間向偏移向量所指的方向移動,移動量為偏移向量的模,而偏移向量的終點代表的是新的質心(即圓心)所在位置,如圖4.當偏移向量的模達到一定的閾值或者說等于0時,即質心的位置不再改變時,聚類結束.

圖4 MSC算法原理示意圖

2.2 加權偏移向量計算

對于傳統的MSC算法,局部樣本空間中所有樣本對偏移向量的計算具有同等權重,當存在噪聲點時,偏移向量的計算就會受到影響,因此在傳統算法中增加核函數,降低噪聲點對偏移向量的影響.在公式(1)中引入核函數,此時偏移向量的計算公式如下:

(2)

式中κ表示核函數.以高斯核函數為例,對應的加權偏移向量如下:

(3)

2.3 Bandwidth估計

MSC算法最重要的參數是帶寬的選擇,算法輸出結果對于帶寬的選擇是敏感的,本文針對CSEM實測數據分布特征提出了一種基于局部密度梯度的bandwidth估計方法.

圖5所示為自適應帶寬估計的過程.對實測CSEM時間序列進行等時長分段,并將所有分段序列進行快速傅里葉變換(FFT),提取不同時間段的頻率信號形成一個m維的數據空間S={s1,s2,…,sm},其中m代表CSEM有效頻率數,任一點sj=(sj1,sj2,…,sji,…,sjn)T,代表第j個頻率的所有觀測數據,n代表樣本點數,其頻率域分布如圖5a所示(36 Hz).Bandwidth估計的具體流程如下:

圖5 帶寬估計過程

(1)計算第j個頻率的觀測樣本集sj中所有樣本點兩兩之間的距離,得到該頻率所有樣本之間的距離矩陣dj如下:

dj=(d1,2,d1,3,d1,4,…,d1,n,d2,3,d2,4,…,d2,n,…,

di,i+1,di,i+2,…,di,j,…,di,n,…,dn-1,n)T

(2)將矩陣dj中的元素按從小到大的順序排列,得到新的距離樣本矩陣d′j,如圖5b,矩陣d′j中的任意樣本所在的位置用R表示:

(4)

式中N(d′0,d′i)表示任意樣本d′i到最小樣本d′0之間的樣本長度;N(d′j)表示距離樣本矩陣d′j的長度.則m個頻率對應的所有距離樣本矩陣d如下:

d=(d′1,d′2,…,d′j,…,d′m).

(3)將L=max(d′j)-min(d′j)進行分段處理,計算每一小段對應的樣本局部密度δi,i=1,2,…,N′,其中N′為分段個數,

(5)

(4)對密度距離分布進行處理,計算局部密度隨距離下降最快的位置.

首先對上述密度分布進行多項式擬合,如圖5d,然后求解擬合結果的梯度,得到min(dδi).

(5)估計帶寬.本文處理的實測數據在未受到電磁干擾或受到電磁干擾較小時,其在頻率域呈現高斯分布特征,而當受到強電磁干擾時,部分樣本會大尺度偏離樣本分布中心.大量實測數據分析結果表明帶寬的選擇應該以多數數據分布范圍為宜,即求取密度下降最快位置對應的距離作為聚類的帶寬,如圖5e所示.

h=distance{min(dδi)}.

(6)

(6)以h為帶寬進行聚類,結果如圖5f所示.

2.4 算法流程

綜上所述,基于2.3節最優帶寬進行MSC聚類的步驟如下:

(1)計算偏移均值向量.

給定樣本數據集S={s1,s2,…,sn},選取任意點sj,計算其偏移均值向量:

(7)

其中Sh表示以sj為圓心,以h為半徑的數據區域,m為區域內樣本點的個數,M(sj)為樣本點sj為起點到其數據空間Sh內所有點的向量均值,稱為偏移均值向量.其中Sh中的數據s′滿足以下條件:

Sh(s′)={s′|(s′-sj)T(s′-sj)

(8)

(2)更新中心點.

以M(sj)向量的終點作為新的圓心,重新構建數據空間計算新的均值向量,此時的圓心sj+1由第j次的圓心sj移動到第j次的偏移均值向量的終點位置,移動距離為偏移均值向量的模.

sj+1=M(sj)+sj.

(9)

(3)重復步驟(1)和(2)直到偏移均值向量的模滿足設定的閾值或者不變為止,本文設定閾值為10-6,當均值向量的模小于10-6時,停止聚類.

(4)樣本點分類.根據每個劃分類對每個樣本點的訪問次數對樣本點進行分類,最終輸出劃分簇C={C1,C2,…,Ck},其中k為聚類簇數,需要注意的是與常規的K-means方法不同,此時的k值是在聚類過程中根據樣本數據的分布情況以及所選定的帶寬確定的,并不是事先預設的.

MSC算法是通過從低密度區域向高密度樣本區域滑動對樣本進行劃分,因此算法本身具有自適應優化的特點,在處理的過程中會對離群點進行標注,對于大尺度非周期噪聲的影響具有很好的處理效果.

3 仿真實驗分析

3.1 仿真信號

通過對海量實測數據分析可知,CSEM實測數據受到的典型非周期噪聲有方波、脈沖、三角波、衰減信號等多種類型,本文以WFEM為例,通過在偽隨機時間域信號中疊加典型噪聲,構造了一組仿真信號,其中噪聲類型有方波、脈沖、衰減、振蕩衰減和高斯噪聲等,噪聲信號的幅度是有效信號的幾十倍,甚至上百倍,如圖6所示.圖6a為加噪前后的時間域波形,藍色為偽隨機信號,黑色為加噪后的偽隨機信號.有效信號包含1 Hz、2 Hz、4 Hz、8 Hz、16 Hz、32 Hz以及64 Hz共計7個主要頻率.信號時長100 s,采樣率4096 Hz.加噪后,每個周期的信號均受到了影響.將加入噪聲的仿真信號整體進行傅里葉變換之后,對應的頻率-振幅曲線如圖6b所示,藍色為加噪前,黑色為加噪后.加入噪聲后,所有頻率均受到了影響,尤其是1 Hz信號受干擾程度最嚴重.圖6c為仿真數據加噪后不同頻率的100個樣本數據的振幅分布.由圖可知,100個數據的振幅均受到不同程度的影響.

圖6 仿真信號

3.2 不同核函數處理效果

核函數的種類有很多,理論上任意一個具有半正定核矩陣的對稱函數均可作為核函數.本文引入了Gaussian核函數、Triangle核函數和Triweight核函數,見公式(10)—(12),上述三種核函數給予中心點的權重值各不相同.圖7為三種核函數在示值區間內的曲線對比,其中高斯核函數的中心為0,帶寬為h=1.對比圖中的三種核函數可知,核函數對于越接近中心點的樣本給予的權重越高,在相同帶寬的條件下,樣本越靠近中心點,Triweight核給予樣本的權重最高,其次是Triangle核,Gaussian核的權重值最小.

Gaussian核函數表達式:

(10)

Triangle核函數表達式:

y=1-|x|, (|x|≤1).

(11)

Triweight核函數表達式:

(12)

為分析不同核函數對于仿真數據的處理結果,本文固定帶寬選擇位置分別為R=0.7、0.5和0.3,表1為不同頻率在不同R時對應的帶寬.圖8—10為不同R時,未加核函數和加入不同核函數的處理效果.

表1 不同頻率在不同R時對應的帶寬

圖8 R=0.7時不同核函數處理后的頻率-電場曲線

圖9 R=0.5時不同核函數處理后的頻率-電場曲線

圖10 R=0.3時不同核函數處理后的頻率-電場曲線

分析圖8—10可知,在加入不同的核函數時,處理結果均獲得不同程度的改善.表2為加入核函數后處理結果與真值之間的均方誤差,對比幾種核函數處理均方誤差可知,Triweight核函數的處理結果均方誤差最小,說明對于仿真數據而言,Triweight核函數的處理效果最好;其次是Triangle和Gaussian核函數.帶寬選擇越大,加入核函數的處理效果越明顯,說明在傳統MSC算法中引入核函數能有效降低帶寬選擇的敏感度,提高算法的穩健性.

表2 加入不同核函數前后處理結果均方誤差

3.3 帶寬估計處理結果

對上述仿真信號,采用本文提出的基于局部密度梯度的方法進行帶寬估計,并利用估計帶寬進行聚類計算得到的結果如圖11所示.

圖11 本文方法處理后的頻率-電場曲線

分析圖11可知,在未加入核函數時,本文帶寬估計方法的處理結果誤差為0.002986,在加入Gaussian核函數時,處理結果有所改善,但在加入其他兩種核函數時,處理結果比加入之前的要差,說明本文方法獲得的帶寬已基本達到最優,且此帶寬獲得的局部空間數據的分布接近于高斯分布(只受高斯噪聲的影響),因此給予中心位置更高權重的核函數處理效果反而不好,而Gaussian核函數對于服從高斯分布的數據體處理效果會比較好,如圖12所示為處理前后仿真信號的頻率域分布.

圖12 仿真數據處理前后的頻率域分布

4 實例分析

為驗證本文方法對CSEM實測信號的處理效果,選擇山東某城區長時段采集的WFEM信號進行處理,該信號利用高階偽隨機信號作為激勵源,相同采集時間,數據量得到了成倍的提升,尤其是低頻段(山東大學等, 2021;Yang et al., 2022b).城區由于強電磁干擾影響,原始數據信噪比低.圖13為研究區內某一實測點WF-4的原始數據,其中包括時間域波形、頻譜和時頻譜.從其時間域波形(圖13a)中可以看出在7500～8000 s之間該測點受到了大尺度的方波干擾,結合其時頻譜(圖13c)可以看出,此時段0.1～500 Hz之間能量均很強;尤其是在低頻段(<50 Hz),有效信號幾乎完全被淹沒,如果采用全時段數據進行有效信號提取,無疑會引入較大的觀測誤差.圖14a為不同時長分段數據進行FFT變換后通過振幅平均獲得的歸一化電場;圖14b為剔除7500～8000 s大尺度方波干擾后,不同時長分段數據通過振幅平均獲得的歸一化電場.從結果中可以看出,采用不同時長進行數據分段,處理結果差異較大.受大尺度干擾的影響,原始數據全時段獲得的頻率-歸一化電場曲線畸變嚴重,而剔除強干擾后的全時段曲線除低頻段外整體較為光滑連續;隨著分段時長減小,噪聲影響范圍逐漸縮小,曲線形態趨于穩定,但經過振幅平均的方法獲得的歸一化電場是將強噪聲影響進行平均,獲得的結果依然存在很大誤差.

圖13 實測點WF-4的原始數據

圖14 實測點WF-4不同時長分段數據取平均獲得的歸一化振幅

圖15為本文方法處理后部分頻率的聚類結果,圖中不同顏色代表不同的簇.由圖可知,受強噪聲影響較大的數據樣本單獨成簇,通過本文方法能將噪聲點進行有效地標注、剔除,降低了強噪聲對有效信號的影響.表3為WF-4測點處理前后部分頻率的樣本數與均方誤差對照表.由表可知,處理前樣本數均為2048個,處理后樣本數占比在91%以上,均方誤差均小于5%,說明本文帶寬估計方法能盡可能保留受干擾較小的數據樣本,且處理后的數據質量得到了有效的改善.圖16為WF-4測點處理前后頻率域分布圖,從圖中可以看出,處理前如圖16a,受大尺度噪聲影響,WF-4測點數據頻率域分布出現大量離群點,尤其是低頻段,結合圖15中部分頻率分布可知,0.25 Hz與0.5 Hz部分時段振幅已大于100 mV,超出實測信號強度上千倍,說明WF-4測點受噪聲影響嚴重.采用本文方法處理后的頻率域分布如圖16b所示,信號主體分布在0.05 mV范圍內,并呈現高斯分布特征.

表3 處理前后不同頻率對應的樣本數量及均方誤差

圖15 本文方法處理后部分頻率的聚類結果

圖16 本文方法處理前后測點WF-4頻率域分布

圖17為測點WF-4去除大尺度干擾時段前后,采用本文方法與Robust估計方法進行處理的頻率-歸一化電場曲線對比圖.從圖中可以看出,原始數據受噪聲影響,多數頻點產生畸變,曲線整體形態不清;采用本文方法處理后,頻點連續性好,曲線形態光滑連續,說明本文方法在處理強電磁噪聲方面具有良好的效果.

圖17 本文方法與Robust估計方法處理后測點WF-4歸一化電場對比曲線

圖18為采用本文方法對部分實測WFEM測點處理前后的頻率-視電阻率曲線對比圖,其中黑色為原始數據,紅色為處理后,從圖中可以看出,原始數據受噪聲影響較大,部分頻率畸變嚴重,曲線形態不清晰.采用本文方法處理之后數據得到了很好的改善,曲線整體光滑連續.

圖18 部分測點采用本文方法處理前后頻率-視電阻率曲線對比

5 結論

(1) 本文針對CSEM信號周期性以及噪聲影響的復雜性、隨機性等特點,將不同時間的觀測數據在頻率域進行量化,用樣本間的相對距離作為衡量噪聲對有效信號影響程度的指標,利用均值漂移聚類(MSC)算法實現了復平面內的CSEM噪聲分離.

(2) 傳統MSC算法對帶寬選擇具有較高的敏感性,因此本文在傳統MSC算法基礎上進行了改進,通過增加核函數有效降低聚類結果對帶寬選擇的敏感度,提高了算法的穩健性,并通過不同的樣本比例R驗證了核函數算法的處理效果.

(3) 對于帶寬參數選擇的復雜性與不可控性,本文提出了基于樣本局部密度梯度的自適應帶寬估計方法,該方法實現了帶寬的自適應最優化估計,在一定程度上提高了算法的處理效果與計算效率.

(4) 通過本文方法對仿真數據與實測數據進行了處理,仿真數據處理后均方誤差小于0.3%,實測數據處理后均方誤差小于4%.說明本文提出的加權自適應帶寬均值漂移算法能夠有效降低噪聲對有效信號的影響程度,提高數據信噪比與處理效率,是一種較有效的CSEM信噪分離方法.