戚庭野，衛(wèi)會汝，馮國瑞，張新軍，余傳濤，趙德康，杜孫穩(wěn)

(1. 太原理工大學礦業(yè)工程學院，山西太原，030024；2. 山西省綠色采礦工程技術研究中心，山西太原，030024)

隨著我國礦產資源的開發(fā)利用，地下資源開采后遺留的采空區(qū)面積逐年增加。地下采空區(qū)坍塌往往會造成巨大經濟損失和人員傷亡[1?2]，因此，準確探測并掌握采空區(qū)的空間分布對礦山的發(fā)展和礦區(qū)城市化進程具有重要意義[3]。目前，采空區(qū)探測手段眾多，其中瞬變電磁法(transient electromagnetic method，TEM)因其勘探深度大、體積效應小、對低阻層分辨率高和對高阻層探測能力強被廣泛應用[4?5]。瞬變電磁法通過對接收的二次感應電動勢計算反演得到相關的地層信息。因此，接收到的感應電動勢信號質量直接關系到解釋工作的可靠性。然而，在野外勘探環(huán)境中，瞬變電磁信號不可避免地被周圍的高壓電線、儀器設備、作業(yè)車輛以及自身傳感器產生的電磁噪聲影響，導致很難從接收的信號中提取出真實的感應電動勢的衰減特征，目標區(qū)域的地質信息得不到有效反映[5]。因此，有必要提出一種濾除信號中的噪聲成分的方法，這成為采空區(qū)準確預測的前提條件。

當前研究中，研究者根據(jù)瞬變電磁衰減曲線的特點，提出了諸多降噪方法，其中，時頻分析方法是有效識別和分離瞬變電磁信號噪聲的主要手段之一。DAI 等[6?8]根據(jù)理論電磁響應的瞬時特性和實測噪聲特點，研究了小波變換方法在抑制電磁干擾方面的優(yōu)勢。雖然單獨使用小波變換可以顯著降低噪聲，但細節(jié)系數(shù)中保留噪聲仍需要聯(lián)合其他方法進行解決。解海軍等[9]采用混合小波去噪方法克服了傳統(tǒng)閾值和單一閾值去噪方法的缺陷，在提高瞬變電磁信號信噪比的同時保證了數(shù)據(jù)質量，但沒有解決小波基選取的問題。此外，經驗模態(tài)分解(EMD)方法能自適應對非平穩(wěn)、非線性信號進行分解，在瞬變電磁信號去噪領域得到廣泛應用。

WANG 等[10]利用EMD 方法將瞬變電磁信號自適應分解為一系列的IMF 分量，通過功率檢測自適應識別噪聲獲得有效信號，但沒有考慮EMD 存在的端點效應和模態(tài)混疊現(xiàn)象。朱通等[11]在EMD算法的基礎上，利用多項式擬合延拓來壓制信號的端點失真問題，可有效地識別電磁信號中的諧波噪聲。LIU 等[12]將改進的EMD 算法即集成經驗模態(tài)分解算法(EEMD)應用到瞬變電磁信號降噪領域，成功去除了信號中的運動噪聲，但忽略了添加的輔助高斯白噪聲的分解殘留問題。

鑒于 EMD 及其改進算法的不足，DRAGOMIRETSKIY 等[13]提出一種非遞歸式自適應分解時間序列信號的變分模態(tài)分解(variational mode decomposition，VMD)算法。該算法具有堅實的數(shù)學理論，魯棒性強，避免了端點效應和模態(tài)混疊問題，但分解結果受分量個數(shù)K和懲罰因子α的限制。為避免人為經驗選取參數(shù)造成的時間損失和效果偏差，本文將MIRJALILI 等[14]提出的鯨魚優(yōu)化算法(whale optimization algorithm，WOA)用于VMD 參數(shù)的優(yōu)化。相比傳統(tǒng)的優(yōu)化算法[15]，鯨魚優(yōu)化算法不僅需要調節(jié)的參數(shù)少，尋優(yōu)精度高，收斂速度快，而且全局尋優(yōu)能力強，在許多領域得到廣泛應用[16]?；诖耍疚奶岢鲆环N基于WOA-VMD 的瞬變電磁信號噪聲處理方法，通過仿真模擬和現(xiàn)場實測數(shù)據(jù)處理，利用該方法對瞬變電磁仿真信號和實測信號進行自適應分解，并且與基于EEMD 方法的分解結果對比，以檢驗本文所提出的方法的有效性和可行性。

1 降噪算法

1.1 VMD算法

VMD(variational mode decomposition)是一種將信號非遞歸分解成數(shù)個具有準正交性的帶限本征模態(tài)函數(shù)的算法，其基本原理可以表示為一個約束變分問題的求解，構造的約束變分問題表示為

采用乘子交替方向算法交替更新得到K個IMF及其對應的中心頻率。

從式(2)可以看出：K和α影響著VMD 的分解性能。若設置的K較小，則信號的多個分量可能同時包含在1 個模態(tài)中；若K較大，則會導致1 個分量包含在多個模態(tài)中，迭代得到的中心頻率也會重疊。對α而言，若α很大，則帶寬限制就會很窄，從而導致有用的頻率成分被消除；反之，冗余頻率成分將會被保留下來。因此，本文提出鯨魚優(yōu)化算法來優(yōu)化確定最佳參數(shù)組合(K，α)。

1.2 WOA算法

鯨魚優(yōu)化算法是一種基于座頭鯨種群迭代進化搜索的元啟發(fā)式優(yōu)化算法。通過模擬座頭鯨獨特的氣泡網攻擊狩獵行為，實現(xiàn)復雜優(yōu)化問題尋優(yōu)，具體優(yōu)化機理包括局部開發(fā)階段和全局探索階段2個部分[14]。

1)在局部開發(fā)階段，鯨魚假設當前的獵物位置為目標位置，在收縮包圍機制和螺旋更新機制之間以50%的概率不斷靠近獵物獲得最優(yōu)解，其數(shù)學模型如下：

2)在全局探索階段，鯨魚以其他鯨魚位置為參考，不斷更新自身位置進行隨機搜索，其數(shù)學模型可表示為

式中：A的絕對值大于1；Xrand為隨機選擇的鯨魚位置。

1.3 基于WOA的VMD參數(shù)優(yōu)化

本文提出采用鯨魚優(yōu)化算法實現(xiàn)變分模態(tài)分解算法的參數(shù)優(yōu)化。該算法中，參數(shù)優(yōu)化的關鍵是選擇適應度函數(shù)。排列熵作為一種度量混沌時間序列復雜性的平均熵參數(shù)，時間序列復雜度參數(shù)計算更簡單，抗干擾能力更強，具有更好地魯棒性的特點，尤其適用于存在動態(tài)噪聲和觀測噪聲的信號[17]?；诖耍疚倪x擇排列熵作為WOAVMD算法的適應度函數(shù)，具體構建步驟如下。

1)給定1個離散時間序列{x(i)，i=1～N}，對該序列進行相空間重構，得到如下重構矩陣：

式中：r=N-(m- 1)τ，r為重構行向量的個數(shù)；j為重構矩陣的第j行分量，j=1～r；m為嵌入維數(shù)；τ為延遲時間。

2)將重構矩陣中的r個行向量分別按照元素數(shù)值進行升序排列，得到各自的符號序列S(q)=(j1,j2,…,jm)。其中，q=1～r，r≤m！；m！為m維相空間映射的符號序列的總數(shù)；j1，j1，…，jm表示各元素在原重構分量中的索引號。

3) 計算r種不同符號序列S(q)出現(xiàn)的概率P1，P2，…，Pr，則時間序列{x(i)，i=1～N}的排列熵定義為

基于上述提出的適應度函數(shù)，鯨魚優(yōu)化算法和變分模態(tài)分解算法聯(lián)合后的信號優(yōu)化分解的流程圖如圖1所示，具體步驟如下。

圖1 WOA算法優(yōu)化VMD參數(shù)的流程圖Fig.1 Flow chart of WOA algorithm to optimize VMD parameters

1)輸入信號，設置VMD 算法中K和α的參數(shù)范圍，初始化WOA模型中的各項參數(shù)，包括種群規(guī)模、最大迭代次數(shù)、空間維度以及初始種群個體。

2)對信號進行VMD分解，利用式(6)計算初始種群中每個個體的適應度。作為該方法的適應度函數(shù)，排列熵被用來衡量參數(shù)組合的分解效果。熵越小，則時間序列分布越有規(guī)律，表示VMD處理得到的IMF 包含更多的有效信息；反之，時間序列越接近隨機分布，IMF 中噪聲成分更多。因此，當排列熵最小時，對應的參數(shù)最優(yōu)。

3)當|A|≤1時，選擇最小排列熵對應的鯨魚位置作為局部開發(fā)的目標值，然后根據(jù)p，選擇式(3)更新鯨魚個體的位置；當|A|>1時，隨機選擇一個鯨魚的位置，根據(jù)式(4)更新鯨魚個體的位置，保留最優(yōu)適應度及對應的參數(shù)組合。

4)保留更新后的鯨魚種群位置作為新一輪的初始種群，循環(huán)迭代，直到達到所設定的最大迭代次數(shù)為止。

5)輸出最優(yōu)鯨魚個體及對應的適應度。

1.4 相關模態(tài)的選擇

為了有效識別和分離經VMD分解后的IMF分量中的噪聲分量，在對IMF 分量進行時頻分析的基礎上引入相關系數(shù)，進一步判斷IMF 與原信號之間的相關程度。相關系數(shù)C的定義如下[18]：

式中：ui(t)和x(t)為分別表示VMD 分解后的第i個IMF分量和原始瞬變電磁響應電壓信號；E和D表示數(shù)學上的期望和方差。

C越大，對應的IMF 分量與原信號相關性越好。本文首先根據(jù)式(7)計算各個分量與原始信號的相關系數(shù)，然后將相關系數(shù)按照分量順序進行排列，取第1個極值點位置作為有效信號分量和噪聲分量的分界點。由于VMD 分解得到1 組從低頻到高頻的IMF 分量，而EEMD 的頻率隨著IMF 階數(shù)增加而逐漸減小，因此，對于VMD方法選取極值點前的分量將作為有效分量，EEMD 方法極值點前對應的則為噪聲分量，最后將有效分量重構得到去噪后的信號。

2 仿真分析

2.1 仿真信號的建立

為了測試基于WOA-VMD 算法對瞬變電磁含噪信號的降噪性能，進行仿真試驗。瞬變電磁現(xiàn)場探測經常受到背景噪聲、工頻噪聲及其諧波干擾的影響。

首先，為了仿真出接近真實工況的含噪瞬變電磁信號，基于Gaver-Stehfest 算法構建電阻率為20 Ω·m 的均勻半空間模型，模擬線圈邊長為30 m中心回線接收的瞬變電磁信號理想衰減曲線。

其次，根據(jù)背景噪聲隨機分布的特點和工頻噪聲及其諧波干擾為固定頻率特征，信號采用高斯函數(shù)和正弦函數(shù)產生3組強度分別為?130，?147和?158 dB的噪聲。

最后，將3 組噪聲添加到瞬變電磁理想信號中產生信噪比(SNR)分別為?5.858，1.642 1 和5.852 6 dB 的含噪信號中，具體的信號計算模型如下：

式中：x(t)為含噪聲干擾的瞬變電磁信號的電壓；f(t)表示理想的瞬變電磁信號的電壓；g(t)為高斯白噪聲；f1為工頻干擾固定頻率50 Hz；f2和f3分別為二次諧波和三次諧波的固定頻率，分別為100 Hz和150 Hz；A1，A2和A3分別為信號振幅的系數(shù)。

2.2 WOA-VMD去噪處理

以3 種不同信噪比的仿真信號驗證WOAVMD 算法的有效性。當含噪信號的信噪比為?5.858 7 dB 時，仿真信號的時域波形及頻譜圖如圖2 所示。從圖2(a)可見：晚期波形出現(xiàn)劇烈震蕩，信號的衰減趨勢被噪聲嚴重干擾，難以準確解譯地層真實地電信息。從圖2(b)可見，相應的頻域變化曲線在50，100 和150 Hz 處發(fā)生異常突起，且高頻區(qū)被毛刺覆蓋，頻譜的分布特征被破壞。

圖2 瞬變電磁含噪信號時域波形及頻譜Fig.2 Time-domain waveform and frequency spectrum of TEM noisy signal

對建立的含噪信號利用WOA-VMD 算法進行降噪處理，使用WOA 對VMD 的參數(shù)(K，α)進行尋優(yōu)，設置種群規(guī)模為50，最大迭代次數(shù)為30，K的迭代范圍為(3，15)，α迭代范圍為(1 000，15 000)，最小排列熵迭代曲線如圖3 所示。由圖3可見：最小排列熵在迭代第9次達到最小值6.426，得到最優(yōu)參數(shù)組合(K，α)=(3，14 886)。

圖3 VMD參數(shù)的WOA優(yōu)化結果(仿真信號)Fig.3 WOA optimization results of VMD parameters(simulation signal)

設置優(yōu)化后VMD的參數(shù)，分解仿真信號得到3 個IMF 分量，如圖4 所示。由圖4 可見：TEM 信號被成功提取到IMF1分量中，工頻及其諧波干擾和高斯白噪聲則混合在其余分量中。計算IMF 分量與TEM 含噪信號之間的相關系數(shù)，結果見表1。

表1 仿真信號的IMF分量相關系數(shù)Table 1 Correlation coefficient of IMFs of simulated signal

圖4 WOA-VMD分解后IMF分量的時域波形Fig.4 Time domain waveform of IMFs after WOA-VMD decomposition

根據(jù)IMF分量的相關系數(shù)在IMF2處首次取得極小值0.072 7，進一步確定IMF1為有效分量，重構后的信號如圖5所示。由圖5可見：降噪后的信號時域波形與原信號波形相符，隱藏于仿真信號中的噪聲信息基本被剔除，TEM 信號的特征信息被很好地保留下來。頻譜圖中固定頻率處的異常突起和高頻處的毛刺也被消除，信號頻域曲線變得更加平滑。但由于TEM 信號為以低頻為主的寬頻信號，頻帶的重合導致VMD分解時仿真信號發(fā)生部分能量損失。值得注意的是，從圖5(a)中可以看出，能量損失對信號分解影響很小，信號重構效果比較理想。

圖5 瞬變電磁含噪信號經WOA-VMD降噪后的時域波形及頻譜Fig.5 Time-domain waveform and frequency spectrum of TEM noisy signal after WOA-VMD denoising

2.3 不同降噪方法對比

為了驗證本文提出的WOA-VMD 算法方法在處理復雜信號的優(yōu)勢，將其與EEMD 方法進行對比。采用EEMD 方法對相同仿真信號進行分解，得到8 個IMF 分量如圖6 所示。VMD 只有3 個分量，說明VMD方法的計算效率更高，信息聚集性能更好[19]。從圖6(a)中僅能分辨出IMF1 為噪聲分量，EEMD 的IMF2～IMF8 分量為無趨勢信號，無法區(qū)分有效信號和噪聲。結合圖6(b)中IMF 分量晚期時域波形及相關系數(shù)的計算結果(見表1)中首個極值點的位置，可以確定IMF1～IMF3 為噪聲分量。從仿真信號中將這些噪聲信號剔除獲得的降噪后信號如圖7所示。從圖7(b)可見：去噪后的頻域曲線較光滑，噪聲成分在一定程度上被去除。然而，觀察圖7(a)發(fā)現(xiàn)仍有部分噪聲殘留在信號中，信號晚期發(fā)生上下波動，主要原因是EEMD在基于極值點定義的局部包絡函數(shù)進行反復遞歸分解過程中發(fā)生包絡估計誤差[20?21]，另外，受端點效應影響，信號早期出現(xiàn)解調誤差，波形發(fā)生輕微凹陷。

圖6 EEMD分解后IMF分量的時域波形Fig.6 Time domain waveform of IMFs after EEMD decomposition

圖7 瞬變電磁含噪信號經EEMD降噪后的時域波形及頻譜Fig.7 Time-domain waveform and frequency spectrum of TEM noisy signal after EEMD denoising

同樣地，采用WOA-VMD 和EEMD 方法的仿真信號進行降噪，結果如圖8 所示。由圖8 可見：雖然2種方法都保留了無干擾信號的衰減趨勢，但經EEMD 去噪后的信號不僅在早期出現(xiàn)了畸變，而且在晚期依舊存在波動，顯然，WOA-VMD 去噪后的平滑效果優(yōu)于EEMD方法。

圖8 瞬變電磁含噪信號降噪結果Fig.8 Denoising results of TEM noisy signals

2.4 仿真試驗結論

為了定量評價WOA-VMD對TEM信號的去噪效果，利用信噪比、均方誤差和平均相對誤差對晚期數(shù)據(jù)進行衡量，2種方法的降噪評價指標見表2。其中，信噪比越大，均方差越小，平均相對誤差越小，去噪效果越好。對比不同噪聲強度情況下VMD 和EEMD 的客觀降噪指標發(fā)現(xiàn)，VMD 方法的總體效果均優(yōu)于EEMD 的分解效果，WOAVMD 方法得到的信噪比最大，得到的均方差和平均相對誤差最小。結合上述仿真結果，證實本文提出的WOA-VMD 算法在處理瞬變電磁信號方面具有優(yōu)越性。

表2 WOA-VMD和EEMD的降噪評價指標Table 2 Denoising evaluation index of WOA-VMD and EEMD

3 現(xiàn)場應用

3.1 采空區(qū)地質條件

為驗證WOA-VMD 降噪算法的實用效果，將該方法用于袁家村鐵礦實測瞬變電磁信號處理。袁家村鐵礦位于山西省呂梁市北部，地層由淺至深依次為新生代第四系、古近系和新近系，下部富集有殘坡積鐵礦礫石、奧陶系下統(tǒng)白云巖、寒武系礫巖鐵礦、灰?guī)r、頁巖交互層以及主要的含鐵巖層上太古界呂梁山群。由于礦山前期被無序開采，導致地下存在大規(guī)模的采空區(qū)，為此采用ATEM-II型瞬變電磁系統(tǒng)探查采空區(qū)分布。

3.2 不同方降噪分析

瞬變電磁探測中，接收到的二次衰減曲線信號受作業(yè)設備以及高壓電力線電磁噪聲的影響發(fā)生劇烈振蕩，晚期信號特征被噪聲完全掩蓋。本文選取某一條測線作為降噪對象，采用WOAVMD方法和EEMD方法分別對該測線上的每個測點進行降噪處理。其中，測點3 實測信號的VMD參數(shù)的WOA 優(yōu)化迭代過程如圖9 所示，在迭代至第5 次時排列熵取得最小值4.958，輸出最優(yōu)參數(shù)組合(K，α)=(2，2 044)。設置相應參數(shù)對信號進行EEMD分解和VMD分解，結果分別如圖10和圖11所示。由圖10 可知：IMF4～IMF8 為虛假分量，IMF1～IMF3 之間發(fā)生模態(tài)混疊，通過對分量進行局部放大判定IMF1 和IMF2 為噪聲分量。由圖11可見：噪聲分量的中心頻率在50 Hz左右，波形變化與工頻噪聲相似。分別剔除這些噪聲分量，重構后的響應數(shù)據(jù)衰減曲線如圖12 所示，信號晚期的噪聲明顯降低。EEMD 降噪曲線在早期由于模態(tài)混疊造成信號畸變，WOA-VMD 的降噪效果更有優(yōu)勢，基本還原了信號的衰減趨勢。

圖9 VMD參數(shù)的WOA優(yōu)化結果(實測信號)Fig.9 WOA optimization results of VMD parameters(measured signal)

圖10 EEMD分解后IMF分量的時?頻圖Fig.10 Time?frequency diagram of IMFs after EEMD decomposition

圖11 WOA-VMD分解后IMF分量的時?頻圖Fig.11 Time?frequency diagram of IMFs after WOA-VMD decomposition

圖12 瞬變電磁實測信號的降噪效果Fig.12 Denoising effect of measured TEM signal

為更直觀獲知本文提出方法的有效性和優(yōu)越性，將測線的11 個測點實測數(shù)據(jù)分別用EEMD 算法和WOA-VMD 算法進行降噪處理，抽取實測信號和降噪信號的26 個測道的感應電壓繪制成感應電動勢多測道曲線，并進行對比，結果如圖13 所示。由圖13(a)可見：信號第2測道數(shù)據(jù)畸變嚴重，晚期測道數(shù)據(jù)由于受噪聲干擾相互交錯嚴重，從圖中無法正確識別采空區(qū)的位置。由圖13(b)可見：經EEMD 處理后的測道曲線的信號早期測道數(shù)據(jù)由于算法無法有效分離不同振蕩特性的信息出現(xiàn)失真，引入了更多的干擾信息，嚴重模糊了信號原有的變化趨勢。由圖13(c)可見：經WOA-VMD處理后，早期突變與晚期的交叉現(xiàn)象都被去除，信號中的噪聲信息得到有效過濾，整體剖面曲線變得更平滑。以上分析表明，WOA-VMD 算法能夠有效抑制信號中的噪聲，較好地恢復信號中的有效信息。

圖13 感應電動勢剖面曲線圖Fig.13 Profile curve of induced electromotive force

3.3 實測降噪數(shù)據(jù)反演

視電阻率等值線剖面直觀反映探測數(shù)據(jù)的解譯結果，因此，對上述實測數(shù)據(jù)和降噪后的數(shù)據(jù)進行反演計算，生成視電阻率剖面圖，如圖14 所示。結合地質資料分析可知，所選取的測線賦存10號礦體，礦體最高標高為1 780 m，最低標高為1 450 m，延伸140～830 m，傾角70°～80°，鉆進中測不到水位，地表亦未見泉水露頭，為透水不含水層。根據(jù)瞬變電磁接收信號可知，在早期，實測數(shù)據(jù)的反演結果(見圖14(a))在距離20 m 和40 m處由于噪聲干擾出現(xiàn)高阻假異?，F(xiàn)象，經EEMD處理后得到的視電阻率剖面圖(見圖14(b))中，在1 480～1 520 m 處出現(xiàn)多處低阻假異常區(qū)，40 m 和80 m 處的圖像質量因模態(tài)混疊而無參考意義，而經WOA-VMD 處理后的剖面圖(見圖14(c))很好地還原了早期的層狀結構。隨著勘察深度增加，晚期數(shù)據(jù)受噪聲的影響增大，從實測數(shù)據(jù)的剖面圖可以看出高阻異常區(qū)散亂分布于1 420～1 500 m，經WOA-VMD 處理后的剖面圖像清晰地反映了高阻異常區(qū)呈層狀分布于1 460～1 500 m。根據(jù)處理結果，在距離60 m 處布設了鉆孔位置，結果顯示在1 470～1 495 m 范圍內出現(xiàn)了掉鉆現(xiàn)象，證明了采空區(qū)的存在，這驗證了本文提出方法的有效性和準確性。

圖14 視電阻率等值線剖面Fig.14 Contour profile of apparent resistivity

4 結論

1)提出一種基于鯨魚優(yōu)化算法改進變分模態(tài)分解的信號降噪方法，該方法很好地解決了瞬變電磁法探測采空區(qū)時存在噪聲干擾的問題。

2)在不同信噪比條件下，改進的VMD方法能夠準確地濾除仿真信號中的工頻干擾和隨機噪聲，更好地恢復整體信號的衰減特征。特別是在定量評價分析中，經WOA-VMD方法去噪后的TEM信號的信噪比分別從?5.858 7，1.642 1和5.852 6 dB提高到?1.251 2，13.546 1和25.101 5 dB，信號質量得到較大改善。

3)WOA-VMD 在瞬變電磁法實測采空區(qū)數(shù)據(jù)降噪中的應用表明，該方法能夠有效去除強干擾信號中的環(huán)境噪聲，提高高阻異常區(qū)的分辨率，得到更準確的采空區(qū)位置及范圍，具有良好現(xiàn)場適用性。