封常青， 李予國,2,3*， 吳云具， 段雙敏,2

1 中國海洋大學(xué)海洋地球科學(xué)學(xué)院， 青島 266100 2 中國海洋大學(xué)海底科學(xué)與探測技術(shù)教育部重點(diǎn)實(shí)驗(yàn)室， 青島 266100 3 青島海洋科學(xué)與技術(shù)試點(diǎn)國家實(shí)驗(yàn)室海洋礦產(chǎn)資源評(píng)價(jià)與探測技術(shù)功能實(shí)驗(yàn)室， 青島 266237

0 引言

大地電磁法(Magnetotelluric, MT)是一種以天然電磁場為場源來研究地球內(nèi)部電性結(jié)構(gòu)的地球物理方法.海洋大地電磁測深(Marine MT)通過坐底式海底電磁采集站(OBEM)測量感應(yīng)電磁場.高導(dǎo)電的海水相當(dāng)于一個(gè)低通濾波器，在一定程度上可以減少人類活動(dòng)所形成的各類人文噪聲，使得海底成為一個(gè)低噪聲的電磁環(huán)境(Jones, 1999).然而，由于引力、地球自轉(zhuǎn)、空氣流動(dòng)等作用，海水處于不停的運(yùn)動(dòng)之中，且其運(yùn)動(dòng)形式多種多樣(如海浪、潮汐和洋流等).海水的運(yùn)動(dòng)會(huì)對(duì)海洋大地電磁探測尤其是淺水區(qū)電磁探測產(chǎn)生嚴(yán)重干擾.一方面，位于海底的電磁接收儀器易受到海水沖刷作用而發(fā)生振動(dòng)(Lezaeta et al., 2005).另一方面，根據(jù)法拉第電磁感應(yīng)定律，運(yùn)動(dòng)的海水如海浪(Weaver, 1965; Fraser, 1966)和潮汐(Larsen, 1968)等因切割地磁場而產(chǎn)生感應(yīng)電流，進(jìn)而產(chǎn)生二次感應(yīng)磁場.它們都會(huì)對(duì)所采集電磁數(shù)據(jù)的質(zhì)量產(chǎn)生影響，降低海洋MT數(shù)據(jù)的信噪比.

在淺水環(huán)境中，海浪運(yùn)動(dòng)產(chǎn)生的感應(yīng)電磁場是影響海洋電磁數(shù)據(jù)質(zhì)量的主要噪聲之一(Lilley et al., 2004).海浪波動(dòng)所產(chǎn)生的磁感應(yīng)強(qiáng)度通常為幾納特到十幾納特，電場強(qiáng)度為幾至幾十微伏每米(Chave, 1983).海浪運(yùn)動(dòng)對(duì)電磁場的影響主要出現(xiàn)在1～30 s的周期范圍內(nèi)，導(dǎo)致此頻段內(nèi)MT數(shù)據(jù)信噪比降低，視電阻率和相位曲線發(fā)生嚴(yán)重畸變(Duan et al., 2020).根據(jù)海浪感應(yīng)電磁場的特性，可以采用信噪分離算法提取海浪感應(yīng)電磁場.于彩霞等(2010)和張寶強(qiáng)(2018)分別將經(jīng)驗(yàn)?zāi)B(tài)分解和小波閾值去噪方法應(yīng)用到淺水區(qū)MT數(shù)據(jù)處理中，視電阻率曲線得到明顯改善，但相位效果不理想.

K-SVD字典學(xué)習(xí)方法作為一種成熟的稀疏表示算法已經(jīng)被成功地應(yīng)用于圖像去噪(Elad and Aharon, 2006)、語音信號(hào)增強(qiáng)(Wang et al., 2016)和地震勘探中隨機(jī)噪聲壓制(Tang et al., 2012)等領(lǐng)域.字典學(xué)習(xí)算法在電磁數(shù)據(jù)去噪中也得到初步應(yīng)用.湯井田等(2018)基于人文噪聲的特點(diǎn)，利用字典學(xué)習(xí)算法構(gòu)建過完備字典，分離音頻大地電磁數(shù)據(jù)中的人文噪聲.Xue等(2020)將K-SVD字典學(xué)習(xí)方法應(yīng)用于時(shí)間域航空電磁數(shù)據(jù)去噪處理中，取得了較好效果.字典學(xué)習(xí)算法在大地電磁信號(hào)(Li et al., 2020, 2021b)和可控源電磁信號(hào)(Li et al., 2021a)噪聲壓制中也得到了應(yīng)用.海浪運(yùn)動(dòng)會(huì)對(duì)淺水區(qū)大地電磁磁場數(shù)據(jù)造成嚴(yán)重干擾(Constable et al., 1998).考慮到海浪感應(yīng)磁場的窄帶譜特性，以及其在時(shí)間域表現(xiàn)出較強(qiáng)的規(guī)律性(張寶強(qiáng), 2018).我們使用K-SVD字典學(xué)習(xí)方法提取海浪感應(yīng)磁噪聲，并結(jié)合視電阻率信息進(jìn)行相位校正，期望提高海洋大地電磁數(shù)據(jù)的質(zhì)量.

在本文中，我們首先闡述K-SVD字典學(xué)習(xí)方法與相位校正的基本原理，其次提出海浪感應(yīng)磁噪聲去噪流程.然后分別將該方法應(yīng)用于仿真數(shù)據(jù)和南黃海實(shí)測大地電磁數(shù)據(jù)噪聲抑制中，并分析其效果.

1 海浪感應(yīng)磁噪聲壓制方法

1.1 K-SVD字典學(xué)習(xí)方法

在信號(hào)處理領(lǐng)域，某些信號(hào)經(jīng)過合適的變換后，可以用一個(gè)僅含有少量非零元素的矩陣進(jìn)行表示，即信號(hào)在變換域中是稀疏的.在稀疏表示算法中，這種變換域被稱為“字典”，其列向量則被稱為“原子”.從綜合角度來看，信號(hào)可以表示為少量字典原子的線性組合，即信號(hào)在字典表示下是稀疏的(練秋生等, 2015).

假設(shè)一維離散信號(hào)為y∈M，并可以表示為

y=Dx+ε,

(1)

式中，x∈K是稀疏的系數(shù)矩陣，D∈M×K稱為字典，這里M和K分別為字典D的行數(shù)和列數(shù).字典D通常是過完備字典，即M

(2)

式中，‖x‖0是x的l0范數(shù)，它表示x中非零元素的數(shù)量，ε為殘差閾值.

在字典D給定的情況下，Mallat和Zhang(1993)采用匹配追蹤算法(Matching Pursuit, MP)求解稀疏系數(shù)x.Pati等(1993)用正交匹配追蹤算法(Orthogonal Matching Pursuit, OMP)求解過完備字典下的稀疏表示問題.傳統(tǒng)的匹配追蹤算法在求解稀疏系數(shù)時(shí)，往往通過事先定義某種數(shù)學(xué)變換或調(diào)和分析方法構(gòu)造過完備字典，例如過完備離散余弦變換(Overcomplete Discrete Cosine Transform, ODCT)字典.雖然該類算法在構(gòu)造字典時(shí)相對(duì)簡單，但是受限于原子固定的形狀，它對(duì)信號(hào)的適應(yīng)性較差.為了提高字典對(duì)信號(hào)的適用性，Aharon等(2006)和Rubinstein等(2008)提出了K-SVD字典學(xué)習(xí)方法.該方法采用數(shù)據(jù)驅(qū)動(dòng)方式直接從訓(xùn)練集中獲得與目標(biāo)信號(hào)特征相符合的原子，使得更新后的字典對(duì)目標(biāo)信號(hào)具有更好的適應(yīng)性.

對(duì)于訓(xùn)練信號(hào)集Y=[y1,y2,…,yN]∈M×N，假設(shè)其每一列信號(hào)均可以表示為字典D中少量原子的線性組合，則稀疏表示問題可以寫成為

(3)

利用K-SVD字典學(xué)習(xí)方法求解(3)式的過程可分為稀疏編碼和字典更新兩個(gè)階段，二者交替進(jìn)行.在稀疏編碼階段，采用正交匹配追蹤算法獲得稀疏系數(shù)矩陣X=[x1,x2,…,xN]∈K×N，在字典更新階段則利用奇異值分解(Singular Value Decomposition, SVD)方法依次對(duì)字典D逐列進(jìn)行更新.我們期望在更新第k列原子時(shí)，(4)式取最小值(Aharon et al., 2006)：

(4)

(5)

定義Ωk∈N×|ωk|，其中在(ωk(i),i)處元素值為1，其余位置元素值為0.式(4)可以等價(jià)為

(6)

下面，舉例說明稀疏表示方法消除信號(hào)噪聲的原理及效果.我們模擬一段多頻信號(hào)：

T=sin(2π×100t+0.1)+3cos(2π×130t+0.9)

+0.7sin(2π×170t+1.5).

(7)

圖1 不同噪聲增益系數(shù)S的重構(gòu)信號(hào)信噪比 當(dāng)噪聲增益系數(shù)S=1.02時(shí)，去噪后的 重構(gòu)信號(hào)信噪比達(dá)到極值13.14 dB.Fig.1 SNR of reconstructed signals with different noise gain coefficient S When S=1.02, the SNR of the reconstructed signal reaches the extreme value of 13.14 dB.

對(duì)模擬得到的信號(hào)加入高斯白噪聲后得到信噪比為0 dB的加噪信號(hào).分別對(duì)無噪信號(hào)、高斯白噪聲和加噪信號(hào)以不同的噪聲增益系數(shù)S進(jìn)行OMP稀疏表示，重構(gòu)信號(hào)的信噪比隨噪聲增益系數(shù)的變化曲線如圖1所示.由圖1可以看出，在重構(gòu)信號(hào)信噪比相同的情況下，無噪信號(hào)所使用的噪聲增益系數(shù)比高斯白噪聲的更大，這說明在相同字典的情況下，噪聲信號(hào)更難以表達(dá).這就是稀疏表示算法可以用于信號(hào)去噪的本質(zhì).而由加噪信號(hào)的重構(gòu)結(jié)果可以看出，過小或過大的噪聲增益系數(shù)均不能實(shí)現(xiàn)最佳去噪效果，存在最優(yōu)的噪聲增益系數(shù)使得重構(gòu)信號(hào)的信噪比取得極值.在實(shí)際去噪實(shí)驗(yàn)中，最優(yōu)的噪聲增益系數(shù)具有未知性.在去噪過程中，通常通過試驗(yàn)不同的增益系數(shù)，并結(jié)合去噪后信號(hào)的時(shí)頻信息(如振幅譜、時(shí)頻圖等)確定最優(yōu)噪聲增益系數(shù).

1.2 基于視電阻率信息的相位校正方法

在利用K-SVD字典學(xué)習(xí)方法實(shí)現(xiàn)信噪分離后，對(duì)去噪后的信號(hào)進(jìn)行Robust阻抗估計(jì)往往可以得到較為光滑的視電阻率曲線，但是所得到的相位曲線仍然顯得不夠理想.Weidelt(1972)、陳樂壽和王天生(1985)論述了大地電磁阻抗視電阻率與相位之間的轉(zhuǎn)換關(guān)系.基于該轉(zhuǎn)換關(guān)系可以實(shí)現(xiàn)相位曲線和視電阻率曲線的相互校正(劉建利, 2013; 楊生等, 2001).因此，我們可以基于校正后的視電阻率信息對(duì)相位進(jìn)行再校正.在一維情況下，大地電磁阻抗視電阻率和相位之間的轉(zhuǎn)換關(guān)系可以表示為(陳樂壽和王天生, 1985)：

(8)

其中，φ(f)為相位，z(g)為阻抗振幅，這里f和g為頻率.由式(8)，可以得到下列近似表達(dá)式(陳樂壽和王天生, 1985)：

(9)

(9)式為視電阻率ρa(bǔ)(f)和相位φ(f)之間的轉(zhuǎn)換關(guān)系式.該關(guān)系式雖然是在一維地電條件下導(dǎo)出的，但它在非一維地電條件下亦基本成立(楊生等, 2001).

用差分近似式(9)中的微分，可得：

(10)

利用上式便可以實(shí)現(xiàn)基于視電阻率信息的相位校正.

1.3 海浪感應(yīng)磁噪聲壓制與相位校正流程

海浪感應(yīng)磁噪聲壓制與相位校正的主要步驟如下：

(1)由于訓(xùn)練字典過程中的輸入信號(hào)為二維訓(xùn)練信號(hào)集，所以需將濾波后的一維時(shí)間序列按一定步長進(jìn)行分段，并將分段信號(hào)作為二維矩陣的列向量進(jìn)行組合；

(2)設(shè)置初始字典D為ODCT字典，利用正交匹配追蹤算法對(duì)二維訓(xùn)練集進(jìn)行稀疏表示，獲得稀疏系數(shù)矩陣X；

(3)利用K-SVD字典學(xué)習(xí)算法更新字典D；

(4)重復(fù)步驟(2)和(3)，在滿足迭代條件后，利用更新的字典DK-SVD和稀疏矩陣XK-SVD重構(gòu)海浪感應(yīng)磁噪聲；

(5)對(duì)去噪后的MT數(shù)據(jù)進(jìn)行Robust阻抗估計(jì)，并結(jié)合視電阻率信息進(jìn)行相位校正.

綜上所述，海浪感應(yīng)磁噪聲壓制與相位校正流程如圖2所示.

圖2 海浪感應(yīng)磁噪聲壓制與相位校正流程圖Fig.2 Flow chart of wave-induced magnetic noise reduction and phase correction

2 模擬數(shù)據(jù)噪聲抑制

采用模擬數(shù)據(jù)檢驗(yàn)本文方法抑制海浪感應(yīng)磁噪聲的效果.為此，首先分別模擬海洋大地電磁信號(hào)時(shí)間序列和海浪感應(yīng)磁場時(shí)間序列.其次，將大地電磁水平磁場分量和海浪感應(yīng)磁噪聲疊加獲得含噪信號(hào)，然后再利用K-SVD方法進(jìn)行噪聲壓制，并與自適應(yīng)噪聲完備集合經(jīng)驗(yàn)?zāi)B(tài)分解(Complete Ensemble Empirical Mode Decomposition with Adaptive Noise, CEEMDAN)(Torres et al., 2011; Colominas et al., 2014)和小波閾值去噪(張寶強(qiáng), 2018)方法的處理結(jié)果進(jìn)行對(duì)比.最后對(duì)去噪后的電磁數(shù)據(jù)進(jìn)行Robust阻抗估計(jì)，并計(jì)算視電阻率和相位.

2.1 模擬海洋MT信號(hào)時(shí)間序列

利用構(gòu)造系統(tǒng)函數(shù)法模擬大地電磁數(shù)據(jù)，構(gòu)造系統(tǒng)函數(shù)可以表示為(Loddo et al., 2002)：

(11)

式中，K1為比例系數(shù)，Ri為零點(diǎn)，Pi為極點(diǎn)(i=1,2,3,4).考慮到海水具有低通濾波的特性，為得到符合海洋大地電磁場特性的時(shí)間序列，經(jīng)過測試實(shí)驗(yàn)，本文選取的系數(shù)如下：

K1=0.1;

R1=0.01,R2=1.05,R3=10.0,R4=15.0;

P1=-0.01,P2=-7.0,P3=-10.0,P4=-2.0.

(12)

利用式(11)計(jì)算得到MT磁場信號(hào)的單邊譜，再將它擴(kuò)展為雙邊譜.然后利用逆傅里葉變換，可得到MT磁信號(hào)時(shí)間序列.在阻抗Z(ω)已知的情況下，再利用阻抗關(guān)系式E(ω)=Z(ω)H(ω)便可以計(jì)算得到電場的頻譜序列.對(duì)電場頻譜序列進(jìn)行逆傅里葉變換得到電場時(shí)間序列.

為了獲得海洋大地電磁模擬數(shù)據(jù)，設(shè)計(jì)如圖3所示兩個(gè)一維地電模型并進(jìn)行正演計(jì)算，將得到的阻抗作為阻抗張量矩陣的非對(duì)角線元素.利用合成的阻抗張量矩陣可以得到四個(gè)電磁場水平分量的時(shí)間序列(湯井田等, 2013; 張寶強(qiáng)等, 2018).假定采樣頻率為10 Hz，采樣時(shí)長為1 h.利用上述方法得到的海洋大地電磁四個(gè)水平分量的時(shí)間序列如圖4所示.

圖3 兩個(gè)一維水平層狀地電模型 (a) 用于模擬Ex與Hy分量； (b) 用于模擬Ey與Hx分量.Fig.3 Two 1-D horizontal layered geoelectric models (a) Used to simulate the Ex and Hy components; (b) Used to simulate the Ey and Hx components.

圖4 模擬得到的四分量海洋大地電磁數(shù)據(jù) (a),(b) 電場水平分量時(shí)間序列； (c),(d) 磁場水平分量時(shí)間序列.Fig.4 Simulated marine MT data (a),(b) Time series of horizontal components of electric field; (c),(d) Time series of horizontal components of magnetic field.

2.2 模擬海浪感應(yīng)磁場時(shí)間序列

根據(jù)法拉第電磁感應(yīng)定律，海浪切割地磁場會(huì)產(chǎn)生感應(yīng)電流，從而產(chǎn)生二次感應(yīng)磁場.海浪運(yùn)動(dòng)感應(yīng)電場E和磁感應(yīng)強(qiáng)度B滿足下列麥克斯韋方程組：

(13)

其中，電流密度J=σmE+σmV×(F+B)，σm為介質(zhì)電導(dǎo)率，V為海浪運(yùn)動(dòng)速度，F(xiàn)為地磁場.由于B遠(yuǎn)小于F，故電流密度可以近似為J=σm(E+V×F).海浪感應(yīng)電磁場的模擬是在Weaver(1965)提出的模型基礎(chǔ)上進(jìn)行的.基于線性疊加理論將同一頻率對(duì)應(yīng)的所有傳播方向組成波產(chǎn)生的感應(yīng)電磁場進(jìn)行疊加，得到不同頻率海浪運(yùn)動(dòng)產(chǎn)生的感應(yīng)電磁場(張寶強(qiáng), 2018).在得到海浪運(yùn)動(dòng)感應(yīng)磁場頻譜后，再對(duì)其進(jìn)行逆傅里葉變換可得到相應(yīng)的時(shí)間序列.

假設(shè)地磁場為50000 nT，風(fēng)速為8 m·s-1，海水深度為50 m，采樣率為10 Hz，采樣時(shí)間為1 h.圖5顯示了模擬得到的海浪感應(yīng)磁場時(shí)間序列和振幅譜(張寶強(qiáng), 2018).模擬海浪感應(yīng)磁噪聲的能量主要集中在0.08～0.2 Hz頻段內(nèi)，表現(xiàn)為窄帶譜特征.在處理過程中，只在上文模擬獲得的MT磁場時(shí)間序列中加入海浪感應(yīng)噪聲.圖6展示了加入海浪感應(yīng)磁噪聲前后的時(shí)間序列和時(shí)頻譜.從圖中可以看到，加噪數(shù)據(jù)時(shí)頻譜在0.1 Hz附近出現(xiàn)了強(qiáng)能量帶.

為了驗(yàn)證K-SVD字典學(xué)習(xí)方法抑制海浪感應(yīng)磁噪聲的效果，將經(jīng)過帶通濾波處理后的數(shù)據(jù)作為輸入數(shù)據(jù)，利用本文提出的方法進(jìn)行去噪處理.K-SVD字典的大小為512×1024，X分量和Y分量數(shù)據(jù)的噪聲增益系數(shù)分別選為2.66和2.25，經(jīng)過10次迭代提取出海浪感應(yīng)磁噪聲.為了對(duì)比起見，分別將CEEMDAN和小波閾值去噪方法應(yīng)用到海浪感應(yīng)磁噪聲壓制中.利用CEEMDAN算法可以將輸入的含噪聲信號(hào)分解成有限個(gè)具有不同頻率成分的固有模態(tài)函數(shù)(Intrinsic Mode Function, IMF)，通過舍棄與海浪感應(yīng)磁噪聲相關(guān)的IMF實(shí)現(xiàn)壓制海浪感應(yīng)磁噪聲的目的(于彩霞等, 2010).在實(shí)驗(yàn)中，假定模態(tài)平均次數(shù)為500.在小波閾值去噪實(shí)驗(yàn)中，選用dmey小波作為基函數(shù)，采用折衷閾值規(guī)則對(duì)經(jīng)過小波分析得到的第5層和第6層細(xì)節(jié)系數(shù)進(jìn)行閾值處理.最終得到的去噪結(jié)果如圖7所示.K-SVD字典學(xué)習(xí)方法能夠較為準(zhǔn)確地提取出海浪感應(yīng)磁噪聲，重構(gòu)后的大地電磁磁信號(hào)誤差也較小.

為了進(jìn)一步衡量去噪效果，引入三個(gè)評(píng)價(jià)指標(biāo)，即重構(gòu)信號(hào)信噪比(SNR)，均方根誤差(RMSE)和歸一化互相關(guān)系數(shù)(NCC)：

(14)

(15)

(16)

其中，c(i)和r(i)分別代表無噪信號(hào)c和去噪信號(hào)r的第i個(gè)元素，L代表信號(hào)長度.分別計(jì)算帶噪聲

圖5 模擬海浪感應(yīng)磁場X分量(a,b)和Y分量(c,d) (a),(c) 時(shí)間序列； (b),(d) 振幅譜.Fig.5 Simulated X component (a,b) and Y component (c,d) of the wave-induced magnetic field (a), (c) Time series; (b), (d) Amplitude spectrum.

圖6 模擬海洋大地電磁Hx分量(a,b,c)和Hy分量(d,e,f)加噪前后對(duì)比 (a),(d) 加噪前后時(shí)間序列； (b),(e) 加噪前時(shí)頻譜； (c),(f) 加噪后時(shí)頻譜.Fig.6 Simulated horizontal magnetic Hx component (a,b,c) and Hy component (d,e,f) (a), (d) Time series of clean signal and noisy signal; (b), (e) Spectrogram of clean signal; (c), (f) Spectrogram of noisy signal.

圖7 模擬海洋大地電磁Hx分量(a,b,c,d)和Hy分量(e,f,g,h)去噪前后時(shí)間序列細(xì)節(jié)對(duì)比Fig.7 Comparison of time series segment before and after denoising of simulated Hx (a,b,c,d) and Hy (e,f,g,h) components

圖8 去噪效果評(píng)價(jià)指標(biāo)對(duì)比 (a) 信號(hào)信噪比(SNR)； (b) 均方根誤差(RMSE)； (c) 歸一化互相關(guān)系數(shù)(NCC).Fig.8 Comparison of evaluation indicators (a) Signal noise ratio; (b) Root mean squared error; (c) Normalized cross correlation.

圖9 模擬海洋大地電磁數(shù)據(jù)視電阻率曲線和相位曲線 圖中灰線和黑線代表無噪數(shù)據(jù)Robust阻抗估計(jì)結(jié)果. (a) 含噪聲信號(hào)； (b) CEEMDAN方法； (c) 小波閾值去噪方法； (d) K-SVD字典學(xué)習(xí)方法.Fig.9 Apparent resistivity and phase curves of simulated MT data The gray and black lines in the Figure represent the Robust impedance estimation results of clean data. (a) Noisy signal; (b) CEEMDAN method; (c) Wavelet method; (d) K-SVD dictionary learning method.

圖10 南黃海測站1實(shí)測海洋大地電磁數(shù)據(jù) 灰圈部分為由海浪運(yùn)動(dòng)感應(yīng)磁噪聲引起的噪聲干擾. (a)—(d) 四分量海洋大地電磁數(shù)據(jù)時(shí)間序列； (e，f) 磁場水平分量振幅譜.Fig.10 Measured marine MT data from site 1 deployed in the South Yellow Sea The interference caused by the wave-induced magnetic noise can be seen clearly in the gray circle. (a)—(d) Time series of marine MT data; (e，f) Amplitude spectrum of horizontal component of magnetic field.

圖11 實(shí)測水平磁場Hx(a,b,c)和Hy(d,e,f)數(shù)據(jù)去噪效果對(duì)比 (a),(d) 去噪前后時(shí)間序列； (b),(e) 去噪前時(shí)頻譜； (c),(f) 去噪后時(shí)頻譜.Fig.11 Denoised results of measured horizontal magnetic Hx (a,b,c) and Hy (d,e,f) components (a),(d) Time series before and after denoising; (b),(e) Spectrogram of noisy signal; (c),(f) Spectrogram of denoised signal.

信號(hào)和經(jīng)過上述方法去噪后信號(hào)的信噪比、均方根誤差以及歸一化互相關(guān)系數(shù)，計(jì)算結(jié)果如圖8所示.在經(jīng)過CEEMDAN和小波閾值去噪方法處理后，MT數(shù)據(jù)的質(zhì)量均得到改善.但經(jīng)過K-SVD字典學(xué)習(xí)方法去噪后的信號(hào)具有更高的信噪比、更小的均方根誤差和更大的歸一化互相關(guān)系數(shù).其中，Hx分量的信噪比由加噪后的-1.3093 dB提高為12.1848 dB，均方根誤差由0.5116減小至0.1082，歸一化互相關(guān)系數(shù)由0.6573提高為0.9694.Hy分量的信噪比由加噪后的1.3777 dB提高為12.5642 dB，均方根誤差由0.3618減小至0.0998，歸一化互相關(guān)系數(shù)由0.7602提高為0.9720.

對(duì)加入強(qiáng)能量噪聲后的電磁數(shù)據(jù)和分別經(jīng)過CEEMDAN、小波閾值去噪和K-SVD字典學(xué)習(xí)方法去噪后的電磁數(shù)據(jù)進(jìn)行Robust阻抗估計(jì)，并計(jì)算視電阻率和相位，計(jì)算結(jié)果如圖9所示.在加入海浪感應(yīng)噪聲后，視電阻率和相位在受干擾的4～10 s周期范圍內(nèi)出現(xiàn)較大程度的畸變.CEEMDAN和小波閾值去噪方法可以在一定程度上抑制視電阻率曲線的畸變.但經(jīng)過K-SVD字典學(xué)習(xí)方法去噪處理后，視電阻率曲線和相位曲線變得更加光滑和連續(xù).可見在強(qiáng)能量海浪感應(yīng)噪聲干擾下，K-SVD字典學(xué)習(xí)方法能取得較好的去噪效果.

3 實(shí)測數(shù)據(jù)噪聲抑制

利用上述方法對(duì)南黃海三個(gè)測點(diǎn)(測站1、測站2和測站3)的實(shí)測海洋大地電磁數(shù)據(jù)進(jìn)行處理分析，進(jìn)一步驗(yàn)證方法的應(yīng)用效果.測區(qū)海域水深為21 m左右.以測站1為例，選擇平潮期采樣率為500 Hz、時(shí)長100 min的實(shí)測數(shù)據(jù)進(jìn)行去噪處理.選擇平潮期數(shù)據(jù)是為了盡量減少由于潮汐運(yùn)動(dòng)所造成的電磁噪聲干擾(徐震寰和李予國, 2019).圖10顯示的是實(shí)測電磁場水平分量的時(shí)間序列和磁場水平分量的振幅譜.由圖10e和圖10f可見，磁場水平分量的振幅譜在頻率0.1～0.3 Hz內(nèi)出現(xiàn)噪聲干擾，其中Hx分量尤為明顯.

利用上文方法對(duì)圖10所示的南黃海實(shí)測數(shù)據(jù)進(jìn)行處理分析，其中Hx分量的帶通濾波范圍為0.1～0.25 Hz，Hy分量的帶通濾波范圍為0.1～0.2 Hz，字典D的大小選為1024×2048.由于實(shí)測數(shù)據(jù)的采樣率為500 Hz，而所要提取的海浪感應(yīng)磁噪聲的主頻在10 s附近.此時(shí)字典中原子(列向量)的長度甚至不能覆蓋一個(gè)海浪感應(yīng)磁信號(hào)的周期.為了從訓(xùn)練集中尋找到與海浪感應(yīng)磁噪聲特征相符合的原子，有必要對(duì)實(shí)測數(shù)據(jù)進(jìn)行降采樣處理.

將經(jīng)過10 Hz降采樣處理后的數(shù)據(jù)作為訓(xùn)練集，在提取到海浪感應(yīng)磁噪聲信號(hào)后，對(duì)其再以500 Hz進(jìn)行升采樣處理.將升采樣數(shù)據(jù)從原始的磁場時(shí)間序列中減去，便可得到去噪后數(shù)據(jù)，去噪結(jié)果如圖11所示.由圖11c可見,去噪后的Hx分量中由海浪感應(yīng)磁場引起的強(qiáng)能量得到壓制，恢復(fù)到正常的數(shù)量級(jí).提取出的海浪感應(yīng)磁噪聲如圖12所示，其時(shí)間序列與模擬的海浪感應(yīng)磁噪聲相似，均具有窄帶譜的特征.而Hy分量由于受到海浪感應(yīng)磁場的影響較小，去噪前后的時(shí)頻譜差異不大.

此外，我們也選取了1號(hào)站位Hx分量部分?jǐn)?shù)據(jù)與南海實(shí)測大地電磁磁場數(shù)據(jù)進(jìn)行了對(duì)比，結(jié)果見圖13.南海大地電磁數(shù)據(jù)采集于深水區(qū)，水深約767 m.與淺水環(huán)境相比，深水區(qū)所獲取的大地電磁數(shù)據(jù)受海浪運(yùn)動(dòng)影響較小，在0.1～0.3 Hz頻率范圍內(nèi)沒有出現(xiàn)明顯的噪聲干擾，數(shù)據(jù)質(zhì)量較高.在經(jīng)過去噪后，南黃海實(shí)測數(shù)據(jù)時(shí)間序列中由海浪感應(yīng)磁噪聲所造成的規(guī)律性正弦型干擾得到壓制.

圖12 提取的海浪感應(yīng)磁噪聲 (a) 時(shí)間序列； (b) 振幅譜.Fig.12 The wave-induced magnetic noise extracted by the K-SVD method (a) Time series; (b) Amplitude spectrum.

同樣地，利用Robust阻抗估計(jì)方法對(duì)去噪前后的實(shí)測數(shù)據(jù)進(jìn)行阻抗估計(jì)并計(jì)算視電阻率和相位.由于相對(duì)較高頻部分(>0.1 Hz)存在隨機(jī)噪聲干擾，在處理中選擇小時(shí)窗進(jìn)行阻抗估計(jì)，從而可以獲得較多窗口疊加.而低頻部分則選擇大時(shí)窗進(jìn)行阻抗估計(jì).最后將這兩部分?jǐn)?shù)據(jù)進(jìn)行拼接處理，得到最終的視電阻率曲線和相位曲線，如圖14所示.圖14a展示的是去噪前的結(jié)果，可見在海浪感應(yīng)磁場的影響下，TE模式和TM模式的視電阻率曲線、相位曲線在周期3～10 s之間均出現(xiàn)了畸變.而其他時(shí)段因?yàn)槭艿降暮＠烁袘?yīng)噪聲干擾較小，視電阻率曲線和相位曲線相對(duì)比較光滑.經(jīng)過K-SVD算法去噪后的結(jié)果見圖14b，視電阻率曲線得到明顯改善，在周期3～10 s之間變得比較光滑.雖然相位曲線相較于原始數(shù)據(jù)的結(jié)果也有了改善，但是在周期4～5 s之間有一個(gè)頻點(diǎn)的結(jié)果不理想.利用式(10)對(duì)其校正，結(jié)果見圖14c.經(jīng)過相位校正后，相位曲線變得更加連續(xù).

利用相同的方法，對(duì)南黃海測站2和測站3的實(shí)測大地電磁數(shù)據(jù)進(jìn)行了去噪處理，經(jīng)Robust阻抗估計(jì)后得到的視電阻率和相位曲線如圖15所示.實(shí)測數(shù)據(jù)的處理結(jié)果表明，K-SVD字典學(xué)習(xí)方法可以有效地壓制海浪感應(yīng)磁噪聲干擾，提高M(jìn)T數(shù)據(jù)的質(zhì)量.同時(shí)能夠基于校正后的視電阻率信息實(shí)現(xiàn)相位校正，最終得到的視電阻率與相位曲線更加連續(xù)、光滑.

4 結(jié)語

淺水區(qū)的大地電磁勘探極易受到海浪感應(yīng)電磁場的影響，海浪對(duì)磁場的影響主要集中在1～30 s的周期范圍內(nèi)，呈現(xiàn)出貫穿全時(shí)段的強(qiáng)能量干擾.為壓制海浪感應(yīng)磁噪聲，本文將K-SVD字典學(xué)習(xí)方法與相位校正方法應(yīng)用于海洋MT數(shù)據(jù)處理中.通過仿真測試以及南黃海實(shí)測數(shù)據(jù)的處理分析，結(jié)果表明該方法能夠有效壓制海浪感應(yīng)磁噪聲，提高M(jìn)T觀測數(shù)據(jù)的質(zhì)量，并改善視電阻率及相位曲線.

圖13 南黃海(水深21 m)與南海(水深767 m)實(shí)測 大地電磁數(shù)據(jù)對(duì)比 (a) 時(shí)間序列； (b) 南黃海測站1部分?jǐn)?shù)據(jù)去噪前后的 振幅譜； (c) 南海數(shù)據(jù)的振幅譜.Fig.13 Comparison of marine MT data from the South Yellow Sea (21 m) and the South China Sea (767 m) (a) Time series; (b) Amplitude spectrum of data from South Yellow Sea; (c) Amplitude spectrum of data from South China Sea.

本文提出的方法適用于提取時(shí)間序列中規(guī)律性信號(hào).在實(shí)際的去噪工作中，通常需要經(jīng)過多次試驗(yàn)確定最優(yōu)噪聲增益系數(shù).如何高效、準(zhǔn)確地尋找最優(yōu)增益系數(shù)是需要進(jìn)一步研究的工作.

致謝感謝三位匿名審稿專家給出的寶貴的意見和建議.本文的研究工作受到國家自然科學(xué)基金項(xiàng)目的資助，在此表示誠摯的感謝.

圖14 測站1實(shí)測數(shù)據(jù)視電阻率和相位曲線 (a) Robust估計(jì)； (b) 經(jīng)K-SVD方法去噪后； (c) 經(jīng)相位校正后.Fig.14 Apparent resistivity and phase curves of measured data from site 1 (a) Robust estimation; (b) K-SVD method; (c) Phase correction.