劉 軍, 陳 磊, 李文燦, 孟憲武, 劉郁聰, 劉迪仁, 徐學恭, 方江雄, 楊 鳳
(1.東華理工大學核技術應用教育部工程研究中心, 南昌 330013; 2.東華理工大學地球物理與測控技術學院, 南昌 330013;3.龍巖煙草工業有限責任公司, 龍巖 364000; 4.北京市京核鑫隆科技有限責任公司, 北京 100101; 5.天津市地震局, 天津 300201)
通過對地球磁場的測量,可用于探明礦產資源分布、分析地質結構和預測地震等。但地磁場是一個微弱的矢量場,由多個磁場疊加構成,且隨時間緩慢變化,同時地磁臺站的觀測環境易受到附近電氣化基礎設施和實際運行的影響(如鐵磁性干擾、車輛干擾等高頻干擾),這些易造成地磁觀測數據的質量下降。因此地磁勘探方法需要對所采集的地磁信號進行降噪處理。
馮紅武等[1]使用FFT處理地磁信號去除信號中的部分噪聲,但當噪聲頻帶與有效信號頻帶相同時,該方法將無法實現噪聲去除;Balster等[2]研究了基于圖像特征的小波壓縮去噪算法;康傳利等[3]對噪聲進行多尺度分類組合濾波達到平滑去噪的目的;Shailesh等[4]采用了非局部均值經驗模態分解技術對心電圖信號進行去噪,還存在模態混疊的問題;蘇小會等[5]研究了基于改進小波閾值對遙測數據進行噪聲壓制;Bhadauria等[6]利用曲線變換和總變異量的自適應融合進行醫學圖像去噪;汪偉明等[7]采用基于小波分析與數學形態學融合算法對地磁信號降噪處理;Hari等[8]基于小波變換和獨立分量分析的核磁共振圖像去噪混合自適應算法;龍虹毓等[9]研究了蟻群優化小波閾值法提取變電設備的狀態信號,但容易陷入局部最優解。
以上方法對地磁信號噪聲壓制處理存在降噪不徹底、有效信息被濾除、適用性和泛化能力差等缺點。基于以上背景,本文研究一種用于地磁信號降噪的混沌蟻群優化小波軟閾值算法,將地磁信號進行小波變換,利用GCV函數選取小波系數濾波閾值,再通過混沌蟻群優化算法獲取最優閾值,以達到較好的降噪效果。
假設有一含噪信號:
x(t)=s(t)+n(t)
(1)
式(1)中:x(t)為含噪信號;s(t)為真實信號;n(t)為噪聲信號。
利用小波分析方法具備的多尺度觀測能力,采用小波軟閾值方法去噪[10],步驟如下:
(1)合理選取小波基和信號分解層數,對含噪信號進行小波分解,得到含噪信號在各個尺度上的小波系數。
(2)選擇軟閾值函數[式(2)],對各分解尺度的高頻小波系數進行閾值化處理,保留有效信號的小波系數,濾除噪聲信號的小波系數。

(2)

(3)對閾值化處理后的小波系數進行小波逆變換,重構獲得去噪后的信號。
閾值λ的確定是閾值化處理的核心。若閾值選取過大,將導致部分噪聲無法濾除,去噪效果差;反之,若閾值選取過小,則將帶來有效信號被濾除的信號失真。
針對小波閾值去噪算法,常用的小波閾值選取法有Visu shrink閾值、Heursure閾值、Sure shrink閾值和Minmax閾值等。但這些方法皆有一重要前提,需要獲取噪聲統計特性,但實際待處理數據的噪聲是未知的,故噪聲統計特性亦無法獲得。
GCV閾值函數則避免了需要獲取噪聲信號統計特征估計的問題,在沒有先驗信息的情況下,去除噪聲的同時仍能保留有效信號[11]。GCV函數表達式如式(3)所示。

(3)
式(3)中:N為高頻小波系數個數;N0為閾值化后被置0的小波系數個數;d為含噪信號的高頻小波系數;ωd為閾值化后的小波系數。
從式(3)可以看出GCV(λ)函數的求解只與輸入和輸出數據有關,無需獲取其他信息。
當GCV(λ)函數值最小時,取得最優小波閾值λ,因此將求閾值λ的最優解的問題轉化為求GCV(λ)函數最小化問題。
蟻群算法是模擬螞蟻根據信息素濃度搜索食物的行為,具有很強的正反饋特性,但易陷入局部最優解。因此在蟻群優化算法中加入混沌搜索,尋找全局最優解。
針對蟻群算法局部最優解問題,本文結合具有均勻遍歷性和隨機性的Kent混沌搜索,對蟻群優化的最優螞蟻個體進行混沌搜索,提高蟻群算法的全局搜索能力[12]。Kent映射為

(4)
式(4)中:zk為混沌變量,是區間(0,1)隨機生成的數;δ∈(0,1)為控制參數,選擇δ=0.4,此時Kent映射概率密度函數在(0,1)內服從均勻分布。
蟻群算法尋優獲得最優螞蟻個體位置后,按式(5)進行混沌搜索,可以得到一個新的螞蟻位置。對式(5)迭代獲得一個混沌序列,將該混沌序列加入原最優解集合作為螞蟻起點。

(5)
式(5)中:X0、Y0為蟻群優化算法的最優螞蟻個體位置;Vx、Vy為調節系數;z(j)為混沌變量。
Vx、Vy在迭代前期選取較大系數,加快收斂的速度;而隨著迭代次數P變大,搜索范圍變小,Vx、Vy也逐漸變小,可找到更精確的解,使算法快速跳出局部最優解。Vx、Vy計算公式為

(6)
式(6)中:|d(j,k)|max為尺度j高頻小波系數絕對值的最大值。
混沌蟻群優化算法具有魯棒性較強,全局搜索性強,算法簡單等優點,采用混沌蟻群優化算法對GCV(λ)函數進行尋優[13],算法中城市為閾值λ,螞蟻為閾值范圍[λmin,λmax]內隨機生成的λk,具體步驟如下:
(1)設定初始參數:螞蟻數為m;最大迭代增加強度系數為Cmax;信息素重要程度參數為α;啟發式因子重要程度參數為β;信息素強度為Q;信息素蒸發系數為ρ。
(2)蟻群初始化,將λk賦予m只螞蟻作為初始位置,初始信息素濃度為τij,混沌搜索迭代次數P。

(7)

(4)記錄本次迭代到訪城市和路徑,選擇最短距離的路徑作為最優路徑。
(5)考慮信息素蒸發,此次循環在整個路徑上的信息素增量,更新信息素,禁忌表清零。

(8)


(9)

(6)對最優螞蟻進行Kent混沌搜索后,按照步驟(3)~(5)找到混沌搜索后的最優螞蟻個體。
(7)判斷是否滿足終止條件,若滿足終止條件,輸出最優閾值;若不滿足,則回到步驟(3)繼續迭代,直到最大迭代次數,輸出小波最優閾值。該算法的流程如圖1所示。

圖1 混沌蟻群優化小波閾值算法流程圖Fig.1 Flow chart of chaotic ant colony optimization wavelet threshold algorithm
為驗證該算法對地磁信號降噪效果,將不同強度噪聲加入一個已知信號中,采用本文算法對其進行降噪處理,計算降噪后信號信噪比(SNR)和均方根誤差(RMSE),分析降噪效果。
先構造1個正弦信號,在正弦信號上加入信噪比為-1、2、5 dB的高斯白噪聲作為原始信號。正弦信號如式(10)所示。
y=sin(5x)
(10)
小波分析以雙正交小波基函數[如式(11)所示]作為小波基,分解層數為3層。蟻群算法中螞蟻數量為50,最大迭代次數為200。


(11)

圖2 合成正弦信號及各算法降噪效果Fig.2 Synthetic sinusoidal signal and noise reduction effect of each algorithm
選用兩種常用閾值法進行去噪對比,閾值選取方法為固定形式閾值和Rigorous Sure(基于Stein無偏似然估計Sure閾值),閾值選取為軟閾值方法。
去噪對比結果如圖2所示。其中Rigorous SURE閾值法去噪后信號與原信號相似度較好,但抖動的幅度較大;固定形式閾值去噪信號更為平滑,抖動較小;本文方法去噪信號最貼近真實信號且信號抖動很小,去噪效果最好。
為了更為直觀地表達去噪效果,引入信噪比(SNR)和均方根誤差(RMSE)作為表征:

(12)

(13)

本文方法與對比方法對含噪的正弦信號處理后的信噪比和均方根誤差如表1所示。
由表1中的信噪比和均方根誤差可知,針對不同噪聲的處理,本文方法的SNR和RMSE都明顯優于所對比的兩種常用方法,其中RMSE遠遠小于對比方法。Rigorous SURE方法的SNR和RMSE均優于固定形式閾值;特別地,當信噪比為5 dB時,經固定形式閾值法處理后的信噪比反而降低,說明方法在去噪的同時,還在很大程度上濾除了有效信號。

表1 不同信噪比的正弦信號去噪后信噪比和均方根誤差
將本文方法應用于實測地磁數據去噪中,實驗數據來源為北京時間2018年8月23日在靜海地震臺實際測量的地磁信號,該信號采樣率為1 min,該地磁信號包含大量的噪聲,嚴重影響了地磁資料后續的數據處理。
同樣將固定形式閾值和Rigorous SURE方法兩種方法作為對比,噪聲壓制處理結果如圖3所示,由圖3可以得出:對比算法濾除了部分噪聲,去噪波形與原始信號相近,但細節抖動較大;本文算法去噪信號更平滑,抖動較小,去噪效果優于對比算法。

圖3 實測地磁信號降噪結果 Fig.3 Measured geomagnetic signal noise reduction results
將實測信號降噪結果60~80 min的局部波形放大,可以更清晰地看出:Rigorous SURE方法去噪信號跟隨原始信號上下抖動較大;固定閾值法和本文方法去噪信號更為接近原始信號;但本文方法去噪信號更平滑,抖動最小。
針對地磁信號噪聲壓制處理存在降噪不徹底、有效信息被濾除、適用性和泛化能力差等缺點。研究一種用于地磁信號降噪的混沌蟻群優化小波軟閾值算法,利用小波函數對地磁信號進行分解,基于廣義交叉驗證GCV閾值選取函數,結合混沌蟻群算法對閾值函數進行尋優獲取小波的最優閾值。
將本文算法應用于處理實際地磁信號的噪聲壓制中,通過對比本文方法與常用方法去噪效果,本文方法不僅較大程度地去除了噪聲,又較好地保留了信號的有效信息。