劉 軍， 陳 磊， 李文燦， 孟憲武， 劉郁聰， 劉迪仁， 徐學(xué)恭， 方江雄， 楊 鳳

(1.東華理工大學(xué)核技術(shù)應(yīng)用教育部工程研究中心， 南昌 330013； 2.東華理工大學(xué)地球物理與測控技術(shù)學(xué)院， 南昌 330013；3.龍巖煙草工業(yè)有限責(zé)任公司， 龍巖 364000； 4.北京市京核鑫隆科技有限責(zé)任公司， 北京 100101； 5.天津市地震局， 天津 300201)

通過對地球磁場的測量，可用于探明礦產(chǎn)資源分布、分析地質(zhì)結(jié)構(gòu)和預(yù)測地震等。但地磁場是一個(gè)微弱的矢量場，由多個(gè)磁場疊加構(gòu)成，且隨時(shí)間緩慢變化，同時(shí)地磁臺站的觀測環(huán)境易受到附近電氣化基礎(chǔ)設(shè)施和實(shí)際運(yùn)行的影響(如鐵磁性干擾、車輛干擾等高頻干擾)，這些易造成地磁觀測數(shù)據(jù)的質(zhì)量下降。因此地磁勘探方法需要對所采集的地磁信號進(jìn)行降噪處理。

馮紅武等[1]使用FFT處理地磁信號去除信號中的部分噪聲，但當(dāng)噪聲頻帶與有效信號頻帶相同時(shí)，該方法將無法實(shí)現(xiàn)噪聲去除；Balster等[2]研究了基于圖像特征的小波壓縮去噪算法；康傳利等[3]對噪聲進(jìn)行多尺度分類組合濾波達(dá)到平滑去噪的目的；Shailesh等[4]采用了非局部均值經(jīng)驗(yàn)?zāi)B(tài)分解技術(shù)對心電圖信號進(jìn)行去噪，還存在模態(tài)混疊的問題；蘇小會(huì)等[5]研究了基于改進(jìn)小波閾值對遙測數(shù)據(jù)進(jìn)行噪聲壓制；Bhadauria等[6]利用曲線變換和總變異量的自適應(yīng)融合進(jìn)行醫(yī)學(xué)圖像去噪；汪偉明等[7]采用基于小波分析與數(shù)學(xué)形態(tài)學(xué)融合算法對地磁信號降噪處理；Hari等[8]基于小波變換和獨(dú)立分量分析的核磁共振圖像去噪混合自適應(yīng)算法；龍虹毓等[9]研究了蟻群優(yōu)化小波閾值法提取變電設(shè)備的狀態(tài)信號，但容易陷入局部最優(yōu)解。

以上方法對地磁信號噪聲壓制處理存在降噪不徹底、有效信息被濾除、適用性和泛化能力差等缺點(diǎn)?；谝陨媳尘埃疚难芯恳环N用于地磁信號降噪的混沌蟻群優(yōu)化小波軟閾值算法，將地磁信號進(jìn)行小波變換，利用GCV函數(shù)選取小波系數(shù)濾波閾值，再通過混沌蟻群優(yōu)化算法獲取最優(yōu)閾值，以達(dá)到較好的降噪效果。

1 小波軟閾值去噪方法

1.1 算法原理

假設(shè)有一含噪信號：

x(t)=s(t)+n(t)

(1)

式(1)中：x(t)為含噪信號;s(t)為真實(shí)信號;n(t)為噪聲信號。

利用小波分析方法具備的多尺度觀測能力，采用小波軟閾值方法去噪[10]，步驟如下：

(1)合理選取小波基和信號分解層數(shù)，對含噪信號進(jìn)行小波分解，得到含噪信號在各個(gè)尺度上的小波系數(shù)。

(2)選擇軟閾值函數(shù)[式(2)]，對各分解尺度的高頻小波系數(shù)進(jìn)行閾值化處理，保留有效信號的小波系數(shù)，濾除噪聲信號的小波系數(shù)。

(2)

(3)對閾值化處理后的小波系數(shù)進(jìn)行小波逆變換，重構(gòu)獲得去噪后的信號。

1.2 閾值函數(shù)

閾值λ的確定是閾值化處理的核心。若閾值選取過大，將導(dǎo)致部分噪聲無法濾除，去噪效果差；反之，若閾值選取過小，則將帶來有效信號被濾除的信號失真。

針對小波閾值去噪算法，常用的小波閾值選取法有Visu shrink閾值、Heursure閾值、Sure shrink閾值和Minmax閾值等。但這些方法皆有一重要前提，需要獲取噪聲統(tǒng)計(jì)特性，但實(shí)際待處理數(shù)據(jù)的噪聲是未知的，故噪聲統(tǒng)計(jì)特性亦無法獲得。

GCV閾值函數(shù)則避免了需要獲取噪聲信號統(tǒng)計(jì)特征估計(jì)的問題，在沒有先驗(yàn)信息的情況下，去除噪聲的同時(shí)仍能保留有效信號[11]。GCV函數(shù)表達(dá)式如式(3)所示。

(3)

式(3)中：N為高頻小波系數(shù)個(gè)數(shù)；N0為閾值化后被置0的小波系數(shù)個(gè)數(shù)；d為含噪信號的高頻小波系數(shù)；ωd為閾值化后的小波系數(shù)。

從式(3)可以看出GCV(λ)函數(shù)的求解只與輸入和輸出數(shù)據(jù)有關(guān)，無需獲取其他信息。

當(dāng)GCV(λ)函數(shù)值最小時(shí)，取得最優(yōu)小波閾值λ，因此將求閾值λ的最優(yōu)解的問題轉(zhuǎn)化為求GCV(λ)函數(shù)最小化問題。

2 混沌蟻群優(yōu)化小波閾值去噪算法

2.1 Kent混沌搜索算法

蟻群算法是模擬螞蟻根據(jù)信息素濃度搜索食物的行為，具有很強(qiáng)的正反饋特性，但易陷入局部最優(yōu)解。因此在蟻群優(yōu)化算法中加入混沌搜索，尋找全局最優(yōu)解。

針對蟻群算法局部最優(yōu)解問題，本文結(jié)合具有均勻遍歷性和隨機(jī)性的Kent混沌搜索，對蟻群優(yōu)化的最優(yōu)螞蟻個(gè)體進(jìn)行混沌搜索，提高蟻群算法的全局搜索能力[12]。Kent映射為

(4)

式(4)中：zk為混沌變量，是區(qū)間(0，1)隨機(jī)生成的數(shù)；δ∈(0,1)為控制參數(shù)，選擇δ=0.4，此時(shí)Kent映射概率密度函數(shù)在(0，1)內(nèi)服從均勻分布。

蟻群算法尋優(yōu)獲得最優(yōu)螞蟻個(gè)體位置后，按式(5)進(jìn)行混沌搜索，可以得到一個(gè)新的螞蟻位置。對式(5)迭代獲得一個(gè)混沌序列，將該混沌序列加入原最優(yōu)解集合作為螞蟻起點(diǎn)。

(5)

式(5)中：X0、Y0為蟻群優(yōu)化算法的最優(yōu)螞蟻個(gè)體位置；Vx、Vy為調(diào)節(jié)系數(shù)；z(j)為混沌變量。

Vx、Vy在迭代前期選取較大系數(shù)，加快收斂的速度；而隨著迭代次數(shù)P變大，搜索范圍變小，Vx、Vy也逐漸變小，可找到更精確的解，使算法快速跳出局部最優(yōu)解。Vx、Vy計(jì)算公式為

(6)

式(6)中：|d(j,k)|max為尺度j高頻小波系數(shù)絕對值的最大值。

2.2 混沌蟻群算法優(yōu)化算法

混沌蟻群優(yōu)化算法具有魯棒性較強(qiáng)，全局搜索性強(qiáng)，算法簡單等優(yōu)點(diǎn)，采用混沌蟻群優(yōu)化算法對GCV(λ)函數(shù)進(jìn)行尋優(yōu)[13]，算法中城市為閾值λ，螞蟻為閾值范圍[λmin,λmax]內(nèi)隨機(jī)生成的λk，具體步驟如下：

(1)設(shè)定初始參數(shù)：螞蟻數(shù)為m；最大迭代增加強(qiáng)度系數(shù)為Cmax；信息素重要程度參數(shù)為α；啟發(fā)式因子重要程度參數(shù)為β；信息素強(qiáng)度為Q；信息素蒸發(fā)系數(shù)為ρ。

(2)蟻群初始化，將λk賦予m只螞蟻?zhàn)鳛槌跏嘉恢?，初始信息素濃度為τij，混沌搜索迭代次數(shù)P。

(7)

(4)記錄本次迭代到訪城市和路徑，選擇最短距離的路徑作為最優(yōu)路徑。

(5)考慮信息素蒸發(fā)，此次循環(huán)在整個(gè)路徑上的信息素增量，更新信息素，禁忌表清零。

(8)

(9)

(6)對最優(yōu)螞蟻進(jìn)行Kent混沌搜索后，按照步驟(3)～(5)找到混沌搜索后的最優(yōu)螞蟻個(gè)體。

(7)判斷是否滿足終止條件，若滿足終止條件，輸出最優(yōu)閾值；若不滿足，則回到步驟(3)繼續(xù)迭代，直到最大迭代次數(shù)，輸出小波最優(yōu)閾值。該算法的流程如圖1所示。

圖1 混沌蟻群優(yōu)化小波閾值算法流程圖Fig.1 Flow chart of chaotic ant colony optimization wavelet threshold algorithm

3 實(shí)驗(yàn)與分析

3.1 合成正弦含噪信號處理

為驗(yàn)證該算法對地磁信號降噪效果，將不同強(qiáng)度噪聲加入一個(gè)已知信號中，采用本文算法對其進(jìn)行降噪處理，計(jì)算降噪后信號信噪比(SNR)和均方根誤差(RMSE)，分析降噪效果。

先構(gòu)造1個(gè)正弦信號，在正弦信號上加入信噪比為-1、2、5 dB的高斯白噪聲作為原始信號。正弦信號如式(10)所示。

y=sin(5x)

(10)

小波分析以雙正交小波基函數(shù)[如式(11)所示]作為小波基，分解層數(shù)為3層。蟻群算法中螞蟻數(shù)量為50，最大迭代次數(shù)為200。

(11)

圖2 合成正弦信號及各算法降噪效果Fig.2 Synthetic sinusoidal signal and noise reduction effect of each algorithm

選用兩種常用閾值法進(jìn)行去噪對比，閾值選取方法為固定形式閾值和Rigorous Sure(基于Stein無偏似然估計(jì)Sure閾值)，閾值選取為軟閾值方法。

去噪對比結(jié)果如圖2所示。其中Rigorous SURE閾值法去噪后信號與原信號相似度較好，但抖動(dòng)的幅度較大；固定形式閾值去噪信號更為平滑，抖動(dòng)較??；本文方法去噪信號最貼近真實(shí)信號且信號抖動(dòng)很小，去噪效果最好。

為了更為直觀地表達(dá)去噪效果，引入信噪比(SNR)和均方根誤差(RMSE)作為表征：

(12)

(13)

本文方法與對比方法對含噪的正弦信號處理后的信噪比和均方根誤差如表1所示。

由表1中的信噪比和均方根誤差可知，針對不同噪聲的處理，本文方法的SNR和RMSE都明顯優(yōu)于所對比的兩種常用方法，其中RMSE遠(yuǎn)遠(yuǎn)小于對比方法。Rigorous SURE方法的SNR和RMSE均優(yōu)于固定形式閾值；特別地，當(dāng)信噪比為5 dB時(shí)，經(jīng)固定形式閾值法處理后的信噪比反而降低，說明方法在去噪的同時(shí)，還在很大程度上濾除了有效信號。

表1 不同信噪比的正弦信號去噪后信噪比和均方根誤差

3.2 實(shí)測地磁信號降噪效果

將本文方法應(yīng)用于實(shí)測地磁數(shù)據(jù)去噪中，實(shí)驗(yàn)數(shù)據(jù)來源為北京時(shí)間2018年8月23日在靜海地震臺實(shí)際測量的地磁信號，該信號采樣率為1 min，該地磁信號包含大量的噪聲，嚴(yán)重影響了地磁資料后續(xù)的數(shù)據(jù)處理。

同樣將固定形式閾值和Rigorous SURE方法兩種方法作為對比，噪聲壓制處理結(jié)果如圖3所示，由圖3可以得出：對比算法濾除了部分噪聲，去噪波形與原始信號相近，但細(xì)節(jié)抖動(dòng)較大；本文算法去噪信號更平滑，抖動(dòng)較小，去噪效果優(yōu)于對比算法。

圖3 實(shí)測地磁信號降噪結(jié)果 Fig.3 Measured geomagnetic signal noise reduction results

將實(shí)測信號降噪結(jié)果60～80 min的局部波形放大，可以更清晰地看出：Rigorous SURE方法去噪信號跟隨原始信號上下抖動(dòng)較大；固定閾值法和本文方法去噪信號更為接近原始信號；但本文方法去噪信號更平滑，抖動(dòng)最小。

4 結(jié)論

針對地磁信號噪聲壓制處理存在降噪不徹底、有效信息被濾除、適用性和泛化能力差等缺點(diǎn)。研究一種用于地磁信號降噪的混沌蟻群優(yōu)化小波軟閾值算法，利用小波函數(shù)對地磁信號進(jìn)行分解，基于廣義交叉驗(yàn)證GCV閾值選取函數(shù)，結(jié)合混沌蟻群算法對閾值函數(shù)進(jìn)行尋優(yōu)獲取小波的最優(yōu)閾值。

將本文算法應(yīng)用于處理實(shí)際地磁信號的噪聲壓制中，通過對比本文方法與常用方法去噪效果，本文方法不僅較大程度地去除了噪聲，又較好地保留了信號的有效信息。