孫公毅， 裴建新,2,3??， 陳家林， 張學慶

(1. 中國海洋大學海洋地球科學學院, 山東 青島 266100； 2. 中國海洋大學海底科學與探測技術教育部重點實驗室, 山東 青島 266100；3. 青島海洋科學與技術試點國家實驗室 海洋礦產(chǎn)資源評價與探測技術功能實驗室， 山東 青島 266237；4. 中國海洋大學工程學院, 山東 青島 266100； 5. 中國海洋大學環(huán)境科學與工程學院, 山東 青島 266100)

中國海域遼闊，海洋資源豐富，海底地質構造類型齊全，隨著對海洋的深入探索，各種水下探測技術不斷進步，對于多物理場探測的要求不斷提高，目前成熟的水下探測方式主要有聲學、光學、和電磁探測等[1-6]。水下電磁探測具有精度高、穩(wěn)定性好，受水文天氣限制相對較小等優(yōu)點，在特定場合具有不可替代的作用[1]。近年來，隨著深遠海探測的需要，我國極低頻電磁信號探測技術作為一個新興的方向，在水下電磁探測中得到快速發(fā)展，未來的極低頻電磁探測在深海探測、深海開發(fā)方面將具有非常廣泛的應用前景。然而，海水運動產(chǎn)生的感應電磁噪聲與水下極低頻電磁信號存在部分頻率范圍的重疊，這成為水下電磁探測信噪比降低的主要因素。為了提高探測質量，在進一步提高發(fā)射天線輻射效率和等效帶寬的同時，壓制海水運動感應電磁噪聲，并提高水下電磁信號信噪比，則是一種可行的方案，且具重要的研究意義[7]。

魏文博等[8]、于彩霞等[9]、Chen Kai等[10]進行了海洋大地電磁法(Magnetotelluric，MT)中的海水運動感應電磁噪聲壓制研究，通過壓制海浪感應電磁噪聲提高了海洋MT數(shù)據(jù)質量；在有源信號探測方面，羅忠濤、常浩等人開展了水下電磁通信中的噪聲壓制相關研究，這些工作一般將噪聲源簡化為低頻噪聲或高斯白噪聲[11-12]，沒有考慮海洋動力環(huán)境產(chǎn)生的影響。另外，因海流運動其振幅及周期特性隨時間變化而變化，產(chǎn)生的海流感應電磁場為非平穩(wěn)信號，常用的頻域濾波或時域濾波方法不利于壓制復雜的海洋環(huán)境噪聲。小波分析方法因其具有多分辨率分析的特性和時頻局部化的優(yōu)勢，適于處理非平穩(wěn)信號，閾值去噪更適于處理有源信號中的非平穩(wěn)噪聲。本文在分析海流感應電磁噪聲對電磁信號影響的基礎上，合成包含MT信號、海流感應電磁噪聲的極低頻探測信號數(shù)據(jù)，通過小波閾值去噪方法壓制實測數(shù)據(jù)中海流感應電磁噪聲，提高資料信噪比，改善水下探測質量。

1 海流感應電磁場噪聲壓制方法

按照海水運動的類型及頻率特征，海水運動感應電磁場信號可簡單分為海浪感應電磁場和海流感應電磁場。Fraser[13]、Andrew[14]、張自力等[15]、張寶強[16]研究表明風生海浪的感應電磁場其影響頻帶范圍較窄，大多分布在0.03～0.15 Hz，海浪運動越強烈其感應電磁場頻率范圍越低[16]，根據(jù)《中華人民共和國無線電頻率劃分規(guī)定》，極低頻頻率范圍是3～30 Hz(Extremely Low Frequency，ELF)，所以一般海況下的海浪運動對水下極低頻信號影響較小。引起海洋表面波動的要素除了風浪還有地震、潮汐等多種要素，這些波動產(chǎn)生多種形式的海流運動，如潮流、風生環(huán)流、密度流、大洋環(huán)流分支等，其周期從0.01 s到百年尺度[17]，其頻率特征使其成為影響水下極低頻電磁信號的重要因素，如潮流產(chǎn)生的感應電磁場其能量強弱及影響頻率范圍與漲落潮密切相關，漲潮與落潮期，潮流速度較大，引起的感應電磁場能量較強[10, 18]，平潮期流速較小，感應電磁場能量較弱。

海水運動的復雜性導致了海水運動感應電磁場的復雜性，通過正演模擬可以有針對性地進行噪聲和信號分析，基于正演模擬和小波系數(shù)分析確定閾值，考慮實際極低頻電磁信號的頻率范圍，結合帶通濾波和小波閾值去噪進行極低頻電磁信號中的海流感應電磁噪聲壓制，以信噪比衡量噪聲壓制的效果，設計的處理流程如圖1所示。

圖1 去噪流程圖Fig.1 Denoising flow chart

結合正演模擬結果，根據(jù)小波系數(shù)密度分布確定閾值及小波基函數(shù)[19]，對小波變換后的小波系數(shù)進行處理，設置閾值界定有效信號和噪聲，將小于閾值的小波系數(shù)充零或收縮，大于閾值的小波系數(shù)保留，最后反變換就得到了去噪后的數(shù)據(jù)[20]。

采用的離散小波變換公式為：

(1)

在小波基函數(shù)不變的情況下，某一分解層小波系數(shù)的大小與信號片段fi在頻率域的振幅ai呈正相關Wi∝ai，即頻率域振幅ai越大，對應小波系數(shù)Wi越大。

某一分解層小波系數(shù)的閾值量化處理過程可以表示為：

Wi(j)=sign(Wi(j))·FT(Wi(j))。

(2)

式中:Wi(j)是第j分解層的小波系數(shù);Wi(j)是閾值量化后的小波系數(shù);sign是符號函數(shù);FT(Wi(j))是閾值函數(shù);T是閾值。本文采用圖2所示的軟硬折衷閾值函數(shù)。

圖2 常用閾值函數(shù)曲線Fig.2 Common threshold function curves

根據(jù)Mallat算法和小波多尺度分解思想[21-22]，分解層數(shù)與頻率的對應關系為：

(3)

式中:fs為信號采樣頻率;Range(j)為第j層細節(jié)系數(shù)對應頻率范圍，所以分解層數(shù)可由目標信號的頻率范圍確定。

在此，為了分析去噪效果，設去噪前后信噪比SNR為：

(4)

其中:Esignal為信號能量;Enoise為一定帶寬內(nèi)背景噪聲的平均能量。

2 合成數(shù)據(jù)實驗

2.1 合成包含海流運動感應電磁場的數(shù)據(jù)

為了進行合成數(shù)據(jù)實驗，需要獲得包含大地電磁信號的海流感應電磁場數(shù)據(jù)。首先，構建大地電磁模擬所需的地電模型，在此以膠州灣海域地質情況為例：膠州灣地區(qū)第四系地層厚度在0.5～30 m不等，主要為海相沉積淤泥及淤泥質地層、河流沖洪積形成粘性土地層及砂層。下部為白堊系沉積巖及青山群火山巖地層，埋深在20～65 m，灣口處為花崗巖，埋深在30～65 m[23-24]，參考孫永福給出的海底土的電阻率[25]，設置一維地電模型如圖3所示?；谠摰仉娔Ｐ?，利用構造系統(tǒng)函數(shù)法[26]，獲得大地電磁正演模擬數(shù)據(jù)水平分量(見圖4(a)、圖5(a))，結合實測數(shù)據(jù)情況，考慮海水運動對海洋磁場信號的影響較大[6]，本文以磁場分量為例進行闡述。

(σ為電導率，h為層厚度。σ indicates conductivity, h indicates layer thickness.)圖3 正演模擬一維層狀電導率模型Fig.3 1D layered conductivity model for forward modeling

為了進一步引入海流感應電磁場的影響，本文根據(jù)該區(qū)域長期觀測電磁數(shù)據(jù)分析結果，選取漲潮期連續(xù)觀測1 h、受海水運動影響明顯、頻率范圍為0.1～10 Hz左右的海洋電磁場數(shù)據(jù)磁場水平分量作為合成數(shù)據(jù)中的海流感應電磁場數(shù)據(jù)，進而可得到極低頻電磁探測的背景場數(shù)據(jù)，如圖4(b)、圖5(b)所示。為了突出極低頻段信號處理效果，在此設置時頻圖的頻率顯示范圍為0.5～10 Hz。

在以上背景電磁場數(shù)據(jù)基礎上加入模擬的極低頻電磁信號，其中包括：發(fā)射頻率為1 Hz的模擬源信號(時長10 min)和發(fā)射頻率為3Hz的模擬源信號(時長30 min)，信號時頻圖如圖4(c)、圖5(c)所示。由圖5(c)可見：由于作為模擬源信號的方波信號是奇諧信號，會產(chǎn)生奇數(shù)倍的諧波信號，即1和3 Hz的極低頻信號分別產(chǎn)生了各自的諧頻信號(3、5、7、9 Hz)，這些信號均受到了海洋環(huán)境電磁場影響，其中1Hz基頻信號所受影響較為嚴重。

((a)正演模擬的MT數(shù)據(jù)磁場水平分量 Magnetic horizontal component of MT data from forward modeling；(b)引入海流感應磁場噪聲的背景磁場數(shù)據(jù) Background magnetic field data with current induced magnetic noise；(c)加入極低頻信號的合成數(shù)據(jù) Synthetic data with ELF signal.)圖4 合成數(shù)據(jù)時間序列Fig.4 Time series of the synthetic data

((a)正演模擬的MT數(shù)據(jù)磁場水平分量 Magnetic horizontal component of MT data from forward modeling；(b)引入海流感應磁場噪聲的背景磁場數(shù)據(jù) Background magnetic field data with current induced magnetic noise；(c)加入極低頻信號的合成數(shù)據(jù)時頻圖 Synthetic data with ELF signal. 振幅取對數(shù)表示。The Amplitude is expressed in logarithm.)圖5 合成數(shù)據(jù)時頻圖Fig.5 Time-frequency spectrum of the synthetic data

2.2 小波系數(shù)分析

根據(jù)公式(3)確定分解層數(shù)為第3、4、5層；按照基于正演模擬確定小波閾值的思路[19]，通過分別對包含有模擬MT信號與模擬極低頻信號的數(shù)據(jù)、海流感應電磁場數(shù)據(jù)及由這些數(shù)據(jù)合成的數(shù)據(jù)進行小波系數(shù)分析，如圖6所示，可見：模擬源信號在MT資料背景中很容易被識別；而海流感應電磁噪聲所對應的小波系數(shù)遠大于MT數(shù)據(jù)和源信號所對應的小波系數(shù)，這表明MT信號和極低頻信號均會被較強的海流噪聲影響所淹沒，那么數(shù)據(jù)處理的關鍵環(huán)節(jié)就集中于海流感應電磁噪聲的壓制?；诤Ａ鞲袘姶旁肼暸cMT數(shù)據(jù)和源信號之和在小波系數(shù)密度分布上的差別，進行噪聲壓制處理。

((a)MT+極低頻信號MT+ELF； (b)海流感應電磁場Current induced electromagnetic field；(c)包含前述數(shù)據(jù)的合成數(shù)據(jù)Synthetic data containing the above data.)圖6 模擬數(shù)據(jù)在目標分解層的小波系數(shù)Fig.6 Wavelet coefficients of simulated data at target decomposition level

在小波系數(shù)分析結果(見圖6)基礎上結合小波系數(shù)密度分布情況(見圖7)，確定小波基函數(shù)及閾值：首先，以不同小波基對信號做小波變換，可得到源信號所代表的小波系數(shù)和背景信號所代表的小波系數(shù)，求取二者之比，在待選用的55個小波基中選取比值最大的小波基作為數(shù)據(jù)處理的小波基函數(shù)[16]，在此選用‘dmey’小波基；然后，閾值的確定以合成數(shù)據(jù)實驗中第3分解層為例(見圖8)，閾值確定方法為：由于MT和極低頻信號相對較弱，小波系數(shù)絕對值最大值為0.676 6，而同一分解層海流數(shù)據(jù)小波系數(shù)絕對值最大值為12.404 2；MT和極低頻信號對應數(shù)據(jù)在第3層的小波系數(shù)密度分布中，大多數(shù)的小波系數(shù)絕對值都是在0.6以內(nèi)，且閾值設置不應大于數(shù)據(jù)在相應分解層數(shù)小波系數(shù)最大值，綜合圖6和圖8信息，確定第3層的閾值為0.6；同理確定第4、第5層閾值分別為1和2。

圖7 小波系數(shù)密度在(a)第3層；(b)第4層；(c)第5層的分布Fig.7 The density distribution of wavelet coefficients in (a) level 3，(b) level 4，(c) level 5

((a)代表MT和極低頻信號 MT and ELF signal;(b)代表海流感應電磁場 Current induced electromagnetic field.)圖8 第3層小波系數(shù)密度分布圖Fig.8 The density distribution of wavelet coefficients in level 3

2.3 效果分析

按照本文小波閾值去噪方法壓制噪聲，對包含1和3 Hz以及其倍頻信號的合成數(shù)據(jù)進行處理，得到處理前后效果對比，如圖9和圖10所示，可見：3 Hz以下的海流感應電磁噪聲得到明顯壓制，淹沒于噪聲的1 Hz頻率的信號也被較好地識別，兩種信號在時頻域的識別能力均有所提高。表1列出了去噪前后3和1 Hz頻率的信號信噪比變化情況，3 Hz信號附近的噪聲被壓制，信噪比從21.81 dB提升到29.33 dB，淹沒在海流感應電磁噪聲中的1 Hz頻率的信號得以顯現(xiàn)，其信噪比也有所提升，從-5.00 dB提升到1.47 dB。

((a)噪聲壓制前數(shù)據(jù)Data before noise suppression；(b)噪聲壓制后結果The result after noise suppression. 振幅取對數(shù)表示。The Amplitude is expressed in logarithm.)圖9 合成數(shù)據(jù)時頻圖Fig.9 Time-frequency spectra of the synthetic data

圖10 噪聲壓制前后振幅譜Fig.10 Amplitude spectrum before and after noise suppression

表1 噪聲壓制前后信噪比(不同頻率信號)Table 1 SNR before and after noise suppression (Different frequency signals) /dB

為進一步分析該方法的有效性，針對1 Hz的極低頻信號，計算了兩種更低信噪比條件下的噪聲壓制結果，處理前后振幅譜以及信噪比如圖11和表2所示，可見：在較低信噪比條件下，該方法仍具備一定的噪聲壓制能力。

((a)信噪比-6.46 dB SNR -6.46 dB; (b)信噪比-8.07 dB SNR -8.07 dB.)圖11 低信噪比合成數(shù)據(jù)噪聲壓制效果(極低頻信號頻率為1 Hz)Fig.11 Noise suppression effect of low SNR synthetic data (ELF signal of 1 Hz)

表2 噪聲壓制前后信噪比(處理前信號的信噪比較低)Table 2 SNR before and after noise suppression (Low SNR before processing) /dB

3 極低頻電磁信號海上試驗

為驗證該方法的有效性，在南黃海海域開展極低頻電磁信號探測試驗。首先，需要了解試驗的背景環(huán)境，利用海底電磁場采集站(OBEM)開展5 d背景場觀測試驗，從試驗數(shù)據(jù)中截取25 h數(shù)據(jù)，獲得如圖12所示背景場時頻圖，可見：周期約為6 h的強能量干擾重復出現(xiàn)，這與試驗區(qū)正規(guī)半日潮的潮流特征相吻合，判斷其為海流運動引起的感應電磁場周期性的體現(xiàn)，其影響頻率范圍從0.01至幾Hz，漲潮期接收到的電磁場信號比平潮期大1個數(shù)量級以上，足以影響水下極低頻電磁信號的頻率。

(振幅取對數(shù)表示。The Amplitude is expressed in logarithm.)圖12 OBEM采集的25 h背景場(磁場)信號時頻圖Fig.12 Time-frequency spectrogram of 25 h background field (magnetic field) signal collected by OBEM

在此基礎上開展水下極低頻電磁信號探測試驗，在海流流速開始增強的時段發(fā)射3、5和7 Hz三個頻率的極低頻電磁信號，各頻率信號均穩(wěn)定發(fā)射5 min，如圖13(a)所示。對于此實測數(shù)據(jù)，以不同小波基對信號做小波變換，可得到源信號所代表的小波系數(shù)和背景信號所代表的小波系數(shù)，求取二者之比，在待選用的55個小波基中選取比值最大的‘bior3.7’小波基作為此次海上試驗數(shù)據(jù)處理的小波基函數(shù)。然后，根據(jù)公式(3)及小波系數(shù)分析方法，確定分解層數(shù)為第4、5層，閾值分別為0.8和2，得到噪聲壓制后結果，如圖13(b)所示，可見，影響極低頻頻段內(nèi)的噪聲已被明顯壓制，信號未受明顯影響。利用振幅譜(見圖14)和公式(4)計算噪聲壓制前后信噪比，得到處理前信噪比約為16.05 dB,處理后信噪比為20.11 dB，信噪比得到了較明顯提升，同時極低頻信號在時頻域內(nèi)的識別能力得到進一步提高。

((a)噪聲壓制前數(shù)據(jù)Data before noise suppression；(b)噪聲壓制后結果The result after noise suppression. 振幅取對數(shù)表示。The Amplitude is expressed in logarithm.)圖13 實測數(shù)據(jù)時頻圖Fig.13 Time-frequency spectra of measured data

圖14 噪聲壓制前后的振幅譜Fig.14 Amplitude spectrum before and after noise suppression

4 結語

在地磁場作用下海流運動產(chǎn)生的感應電磁場是極低頻電磁探測中的重要噪聲來源，其嚴重影響著極低頻水下電磁探測信號的信噪比。本文將海流感應電磁場特征引入小波分析的閾值優(yōu)化中，基于構造系統(tǒng)函數(shù)法獲得MT信號，合成包含海流感應電磁噪聲的極低頻信號，并用于模型測試；通過合成數(shù)據(jù)試驗和南黃海海域電磁信號探測試驗先后對該方法進行了測試，結果表明該方法可有效壓制海流感應電磁場噪聲的影響，并提高了極低頻探測信號的信噪比。