基于變分模態分解的地面磁共振諧波消噪方法

摘要： 針對地面磁共振信號非常弱， 極易受電磁噪聲干擾的問題， 提出一種基于變分模態分解的地面磁共振諧波消噪方法. 該方法采用基于改進變分模態分解的工頻諧波消除方式， 并根據頻譜分析設定模態分量數與初始中心頻率， 解決了常規諧波建模消噪方法僅能處理單次采集數據而導致的運算效率慢等問題. 實驗結果表明， 該方法在多基頻或基頻隨時間變化等復雜噪聲場景中， 得到了良好的諧波分量估計效果， 并可快速、 有效地消除工頻諧波干擾， 大幅度提升了地面磁共振探測數據信噪比.

關鍵詞： 變分模態分解； 地面磁共振； 諧波干擾； 基頻變化

中圖分類號： TP631.2文獻標志碼： A文章編號： 1671-5489（2025）02-0559-08

Harmonic Noise Reduction Method for Surface MagneticResonance Based on Variational Mode Decomposition

WANG Qi1， LIU Zhaowen2， DU Hailong1， XUAN Yubo1， DIAO Shu3

（1. College of Communication Engineering， Jilin University， Changchun 130012， China;2. College of Instrumentation and Electrical Engineering， Jilin

University， Changchun 130061， China;3. School of Control Engineering， Wuxi Institute of Technology， Wuxi 214121， Jiangsu Province， China）

Abstract： Aiming at the problem of very weak surface magnetic resonancesignalsandsusceptibility to electromagnetic noise interference， we proposed a harmonic noise reduction method for surface magnetic resonance based on variational mode decomposition.This method adopted an improved variational mode decomposition-based method for power frequency harmonic elimination， and set the mode number and initial center frequency according to spectral analysis， solving the problem of slow computational efficiency caused by "conventional harmonic modeling denoising methods， which could only handle single-acquisition data. The experimental results show that the method achieves good effect ofharmonic component estimation in complex noise scenarios such asmultiple fundamentalfrequencies or fundamental frequency variations over time，and canquickly and effectively eliminate power frequencyharmonic interference， significantly improving the signal-to-noise ratio of surfacemagnetic resonance detection data.

Keywords： variational mode decomposition; surface magnetic resonance; harmonic interference; fundamental frequency variation

收稿日期： 2024-03-22. 網絡首發日期： 2024-11-25.

第一作者簡介： 王 琦（1988—）， 女， 蒙古族， 博士， 高級工程師， 從事信號與信息處理的研究， E-mail： wangqi0115@jlu.edu.cn.

通信作者簡介： 刁 庶（1986—）， 女， 漢族， 博士， 講師， 從事磁共振探測技術應用的研究， E-mail： diaoshu@jlu.edu.cn.

基金項目： 吉林省教育廳科學技術研究項目（批準號： JJKH20221006KJ）.

網絡首發地址： https：//link.cnki.net/urlid/22.1340.O.20241122.1534.001.

地面磁共振是一種直接探測地下水的地球物理探測方法， 具有直接高效、 定量準確和解釋唯一等優點， 目前廣泛應用于水資源調查、 海水入侵監測和地質災害預警等領域^［1-2］.但磁共振信號非常弱， 通常僅為納伏級， 極易受噪聲干擾， 導致儀器采集的數據信噪比低， 從而嚴重影響了磁共振數據解釋的準確度.

環境噪聲根據其特征和來源的不同， 可分為隨機噪聲、 尖峰噪聲和諧波噪聲3種^［3］.隨機噪聲是一種高斯白噪聲， 通常可用疊加方法或時頻峰值濾波的方法有效抑制. 尖峰噪聲是由云層放電、 汽車打火等現象引起的大幅度瞬時干擾， 可用壓縮小波變換和非線性能量算子等方法去除. 諧波噪聲由50 Hz頻率為基頻的多個諧波組成， 目前最常用的方法是利用諧波建模法消除穩定的單基頻諧波噪聲^［4］.針對基頻隨時間變化以及多基頻的諧波噪聲， 目前主要使用分段諧波建模和多維諧波建模等消噪方法^［5-6］.這些方法可有效消除諧波噪聲， 但缺陷是計算效率低， 對野外實驗測得的大量地面磁共振數據很難實現實時處理.

經驗模態分解（empirical mode decomposition， EMD）是一種自適應信號處理方法， 可將原始信號迭代分解為多個中心頻率的模態， 并提取出有用信號的特征^［7-8］.該方法已被使用在地面磁共振數據的噪聲消除領域中， 但其模態混疊的現象會導致工頻諧波的頻域分離欠佳^［9-11］.文獻［12］提出了一種變分模態分解（variational mode decomposition， VMD）方法， 該方法通過搜尋變分模型最優解確定每個分量的頻率中心及帶寬， 有效分離非平穩信號的各模態分量. 本文將該方法首次應用到磁共振數據中工頻諧波強噪聲的去除， 以大幅度提升數據處理的計算效率. 不同于EMD直接對時域信號的濾波處理， 變分模態分解將信號分解成具有一定帶寬的頻域分量^［13］.因此， 當磁共振數據中含有多基頻諧波或諧波基頻隨時間變化時， 該方法仍有良好的諧波分量估計效果.

進一步， VMD可自適應地將信號分解為一系列模態信號， 但需人為調整先驗參數^［14］.其中， 模態數對信號的分解效果影響較大. 模態數偏小會導致欠分解， 即信號中部分模態不能被完全提取； 而模態數偏大又會導致過分解， 即分解產生信號中并不包含的虛假模態. 同時這些模態分量的初始中心頻率也會影響VMD方法的收斂速度和分解準確性^［15］， 先驗參數的合理設定對VMD分解效果至關重要. 但對地面磁共振數據中的諧波噪聲， 其頻率均為基頻的整數倍， 在將磁共振數據經過低通濾波后， 其諧波分量的個數與中心頻率是完全已知的. 基于此， 本文提出一種基于改進變分模態分解的地面磁共振工頻諧波噪聲消除方法. 通過頻譜分析設立先驗模態參數， 優化諧波噪聲分量的分解效果并快速實現磁共振的信噪比提升.

1 原理與方法

1.1 地面磁共振信號和數據噪聲類型

地下水中的氫質子受到激發后會由低能態躍遷至高能態， 在停止激發后會向外輻射能量， 產生自由感應衰減（free induction decay， FID）信號^［16］， 用公式表示為

VFID=V0e^-t/T2cos（2πfLt+φ），（1）

其中： V0表示初始振幅， 反映含水量大小； T2表示弛豫時間， 反映含水體孔隙度大小； fL表示Larmor頻率， 與地磁場強度相關; φ表示初始相位， 與地下電阻率、 發射接收

線圈電動勢相位差相關. 地面磁共振FID信號通常會淹沒在復雜環境噪聲中. 采集的磁共振數據VMRS用公式表示為

VMRS=VFID+VHarmonic+VSpike+VRandom，（2）

其中VHarmonic為工頻諧波噪聲， VSpike為尖峰噪聲， VRandom為隨機環境噪聲. 諧波噪聲由（50±0.1）Hz為基頻f0的多個諧波累加組成， 表示為

VHarmonic=∑Nn=1Ancos（2πf0t+φn），（3）

其中An和φn分別是第n個諧波的幅值與相位.

1.2 工頻諧波消除

利用地面磁共振地下水探測儀器采集實測含噪數據VMRS， 以激發頻率fT作為頻譜搬移頻率， 將采集的全波數據轉換為復包絡信號V*=V*real+V

*imag. 將復包絡信號通過低通濾波器限制其帶寬， 濾除通帶外的諧波分量. 為避免模態分解的端點效應導致信號失真， 本文對復包絡信號的實部V*real和虛部V*imag分別進行時域延拓， 并對延拓后的數據進行變分模態分解.

對于磁共振數據工頻諧波消噪過程， 變分模態分解涉及的模態數和初始頻率等先驗參數可根據低通濾波器的頻域特征確定. 諧波噪聲分為單基頻、 多基頻和時變基頻， 假設諧波基頻在（f0±0.1）Hz內， 則全波數據的各諧波頻率可粗略表示為fi=nf0， 其中n為正整數， 復包絡信號諧波頻率為f*i=f0-fT. 當低通濾波截止頻率為fLP時， 復包絡信號中只保留了f*ilt;fLP的諧波分量. 因此， 可精準設定殘留諧波的中心頻率和模態分量數， 優化后續VMD分解的去噪效果與計算效率.

將復包絡信號的實部V*real和虛部V*imag分別變分模態分解成多個諧波模態與磁共振信號模態. 分解目標是各模態分量uk的和等于輸入信號， 且各模態分量中心頻率ωk的帶寬和最小， 用公式表示為

min{uk，ωk}∑ktδ（t）+jπt

×uk（t）e^-jωkt22，（4）

s.t. ∑kuk=VVMD，（5）

其中VVMD為用于VMD分解的輸入信號， 即復包絡信號的實部V*real或虛部V*imag.

引入Lagrange乘子λ（t）和二階懲罰因子α， 將上述約束變分問題轉化為無約束變分問題：

L（{uk}，{ωk}，λ）=α∑ktδ（t）+jπt×uk（t）e^-jωkt22+VVMD（t）-∑kuk（t）22+〈λ（t）

，VVMD（t）-∑kuk（t）〉.（6）

利用交替方向乘子法迭代更新得到各分量uk及其中心頻率ωk， 最終求解到無約束模型的解， 迭代過程如下：

uⁿ⁺¹k（ω）←VVMD（ω）-∑i≠kui（ω）+λ（ω）/21+2α（ω-ωk）2，（7）

ωⁿ⁺¹k←∫∞0ωuⁿ⁺¹k（ω）2dω∫∞0uⁿ⁺¹k（ω）2dω，（8）

其中中心頻率不在0附近的模態函數均為諧波噪聲， 將原始含噪數據減去上述模態函數即可得到去除工頻噪聲后的磁共振信號.

2 仿真結果

2.1 VMD消噪過程

生成一組仿真磁共振數據用于展示本文方法的消噪過程. FID信號參數如下： V0= 100 nV， T*2= 0.3 s， fL=2 330 Hz， φ=2.62 rad. 為模擬環境強噪聲， 本文在該數據中添加了工頻諧波和隨機噪聲. 諧波基頻為50 Hz， 起始階數為40階， 共15次諧波. 幅值隨機分布在500～1 000 nV內， 相位分布在-π～π內. 隨機高斯白噪聲水平為500 nV. 假設儀器采樣率為50 kHz， 采樣時間為1 s， 發射頻率為2 331 Hz， 重復采集64次并疊加數據. 通過Hilbert變換、 頻譜搬移和低通濾波， 將該疊加數據轉換為復包絡信號. 含噪數據與仿真信號的實部和虛部如圖1所示. 由圖1可見， 仿真FID信號被淹沒在含噪數據中.

使用變分模態分解方法處理復包絡數據的實部和虛部. 根據低通濾波截止頻率特性， 確定分解模態數和初始中心頻率. 其中， 實部和虛部分別可被分解成16個模態分量， 包括15次諧波噪聲和磁共振信號. 磁共振信號分量的初始中心頻率設為0， 諧波分量初始頻率為諧波頻率減去2 331 Hz（頻譜搬移頻率）. 圖2（A）為實部分解后的磁共振信號實部分量， 圖2（B）～（H）為多個諧波模態分量.

由于每個模態均為具有一定帶寬的頻域分量， 這些諧波分量包括了64次疊加數據中相近頻率的諧波. 因此本文方法可直接處理多次采集的疊加數據， 而不需耗時地分別單獨處理每次采集的數據. 得到的磁共振信號與擬合曲線的實部和虛部如圖3所示， 擬合結果為V0=98.5 nV， T*2=0.3 s， fL= 2 330 Hz， φ= 2.65 rad.

去噪前后的數據頻譜對比如圖4所示. 由圖4可見， 本文方法可有效濾除諧波噪聲和大部分隨機噪聲， 并且不會導致磁共振信號失真.

2.2 與諧波建模消噪結果對比

為展示本文方法在工頻諧波消噪效果上的優勢， 本文在兩種噪聲情形下對比該方法與諧波建模消噪方法的處理結果. 第一種情形： 諧波基頻隨時間緩慢變化， f0（t）=50+N（0，0.02）+0.05t2N（0，0.02）；

第二種情形： 每次采集數據的諧波基頻不同，f0（t）=50+N（0，0.02）.

本文將疊加64次的數據用于去噪.

圖5（A）和（B）分別為當諧波基頻隨時間變化時， 實部和虛部的含噪數據與消噪結果. 由圖5（A）和（B）可見，諧波建模不能準確估計諧波干擾， 導致去噪結果仍含較多噪聲. 而VMD方法可去除所有諧波分量， 并有效壓制隨機噪聲， 分解得到的磁共振信號的實部和虛部均與仿真FID基本重合. 這是因為諧波建模消噪方法只針對數據中只含有單一基頻的噪聲情形. 圖5（C）和（D）分別為當數據含有多個基頻時， 實部和虛部的含噪數據與消噪結果. 由圖5（C）和（D）可見， 諧波建模對多次諧波的估計是錯誤的. 而變分模態分解是將信號分解成具有一定帶寬的頻域分量. 當磁共振數據中諧波基頻隨時間變化時， 該方法可將相近頻率分量分解成一個模態， 因此能準確估計諧波噪聲. 進一步， 當每次采集數據的諧波基頻不同時， VMD去噪效果仍明顯優于諧波建模方法， 且不會使磁共振信號失真. 假設單次采集數據中諧波基頻固定， 使用諧波建模對每次采集數據分別單獨處理也能有良好的去噪效果， 但會需要更多的計算時間.

2.3 不同噪聲水平的測試

為進一步驗證VMD方法的去噪能力， 本文對比了多組不同強度噪聲的處理結果. 仿真FID信號參數設為V0=100 nV， T*2=0.3 s. 圖6（A）和（B）給出了當單次采集諧波噪聲分別為500，1 000，1 500，2 000 nV， 隨機噪聲固定為500 nV時， 100組數據去噪后磁共振信號參數的非線性擬合統計結果. 初始振幅和弛豫時間誤差均集中在2%以下， 少量分布在5%以下. 即便諧波噪聲水平達到2 000 nV， 本文方法仍有較理想的去噪效果， 磁共振參數擬合誤差基本不變. 注意到， 隨著噪聲水平提升， 誤差離群點反而逐漸減少， 這是因為諧波模態的能量集中會使VMD分解方法更清晰地識別諸多模態. 圖6（C）和（D）給出了隨機噪聲為500，1 000，1 500，2 000 nV， 諧波噪聲固定為500 nV時的擬合統計結果.

由圖6（C）和（D）可見， 當隨機噪聲超過1 000 nV時， 部分初始振幅擬合誤差超過10%， 弛豫時間也有較大偏差. 表明VMD方法主要對信號的頻率差異敏感， 而隨機噪聲水平會影響多個模態分量的頻域帶寬尋優， 導致諧波估計存在誤差.

3 實測結果

為驗證本文提出的VMD去噪方法對實際采集磁共振信號進行隨機噪聲壓制的有效性， 在吉林省長春市燒鍋鎮進行探測實驗， 探測點附近噪聲較小， 且有豐富的地下水資源.

野外實驗數據使用JLMRS型磁共振地下水探測儀獲得， 如圖7所示. 測量地點的地磁場強度為54 908 nT， Larmor頻率為2 330 Hz， 地磁傾角為60°， 收發一體線圈尺寸為50 m×50 m， 使用20組不同的脈沖矩激發^[17-18].

圖8為其中一組測量信號的噪聲壓制結果. 圖8（A）為去噪前后的頻譜對比. 由圖8（A）可見， 其他頻點諧波分量被有效消除， Larmon頻率處幅值突出. 由式（1）可知， 理想包絡信號應為平滑曲線. 圖8（B）和（C）分別為實測數據與去噪結果的實部與虛部對比. 由圖8（B）和（C）可見， 消噪后的曲線平滑， 諧波噪聲大幅度衰減. 該去噪結果表明， 本文方法在反演獲得地下水分布信息的過程中有重要優勢.

綜上所述， 針對地面磁共振數據中強噪聲干擾的問題， 本文提出了一種基于改進變分模態分解的工頻諧波消除方法， 大幅度提升了去噪效果與計算效率. 當磁共振數據中含有多基頻諧波或諧波基頻隨時間變化時， 該方法仍有良好的諧波分量估計效果， 顯著優于常規諧波建模. 針對先驗分解參數影響去噪效果的問題， 本文提出了根據頻譜分析設定模態分量數與初始中心頻率的改進方式. 大量仿真測試與實測數據驗證了本文去噪算法的有效性和實用性.

參考文獻

［1］林君. 核磁共振找水技術的研究現狀與發展趨勢［J］.地球物理學進展， 2010， 25（2）： 681-691. （LIN J. Situation and Progress of Nuclear Magnetic Resonance Technique for Groundwater Investigations［J］.Progress in Geophysics， 2010， 25（2）： 681-691.）

［2］張榮， 胡祥云， 楊迪琨， 等. 地面核磁共振技術發展述評［J］.地球物理學進展， 2006， 21（1）： 284-289. （ZHANG R， HU X Y， YANG D K， et al. Review of Development of Surface Nuclear Magnetic Resonance［J］.Progress in Geophysics， 2006， 21（1）： 284-289.）

［3］BEHROOZMAND A， KEATING K， AUKEN E. A Review of the Principles and Applications of the NMR Technique for Near-Surface Characterization［J］.Surveys in Geophysics， 2015， 36（1）： 27-85.

［4］田寶鳳， 朱慧， 易曉峰， 等. 基于諧波建模和自相關的磁共振信號消噪與提取方法研究 ［J］.地球物理學報， 2018， 61（2）： 767-780. （TIAN B F， ZHU H， YI X F， et al. Denoising and Extraction Method of Magnetic Resonance Sounding Signal Based on Adaptive Harmonic Modeling and Autocorrelation［J］.Chinese Journal of Geophysics， 2018， 61（2）： 767-780.）

［5］YI X F， ZHANG J， TIAN B F， et al. Design of Magnetic Resonance Sounding Antenna and Matching Circuit for the Risk Detection of Tunnel Water-Induced Disasters［J］.IEEE Transactions on Instrumentation and Measurement， 2019， 68（12）： 4945-4953.

［6］LARSEN J J， DALGAARD E， AUKEN E， et al. Noise Cancelling of Magnetic Resonance Sounding Signals Combining Model-Based Removal of Powerline Harmonics and Multichannel Wiener Filtering［J］.Geophysical Journal International， 2014， 196（2）： 828-836.

［7］黃建， 胡曉光， 鞏玉楠. 基于經驗模態分解的高壓斷路器機械故障診斷方法［J］.中國電機工程學報， 2011， 31（12）： 108-113. （HUANG J， HU X G， GONG Y N. Mechanical Fault Diagnosis Method for High-Voltage Circuit Breakers Based on Empirical Mode Decomposition［J］.Proceedings of the CSEE， 2011， 31（12）： 108-113.）

［8］BOUDRAA A O， CEXUS J C. EMD-Based Signal Filtering［J］.IEEE Transactions on Instrumentation and Measurement， 2007， 56（6）： 2196-2202.

［9］YU D J， CHENG J S， YANG Y， et al. Application of EMD Method and Hilbert Spectrum to the Fault Diagnosis of Roller Bearings［J］.Mechanical Systems and Signal Processing， 2005， 19（2）： 259-270.

［10］LIN L， JI H B. Signal Feature Extraction Based on an Improved EMD Method［J］.Measurement， 2009， 42（5）： 796-803.

［11］CHENG J S， YU D J， YANG Y， et al. Research on the Intrinsic Mode Function （IMF） Criterion in EMD Method［J］.Mechanical Systems and Signal Processing， 2006， 20（4）： 817-824.

［12］DRAGOMIRETSKIY K， ZOSSO D. Variational Mode Decomposition［J］.IEEE Transactions on Signal Processing， 2013， 62（3）： 531-544.

［13］NAZARI M， SAKHAEI S M. Successive Variational Mode Decomposition［J］.Signal Processing， 2020， 174： 107610-1-107610-10.

［14］GUO Y F， ZHANG Z S. Generalized Variational Mode Decomposition： A Multiscale and Fixed-Frequency Decomposition Algorithm［J］.IEEE Transactions on Instrumentation and Measurement， 2021， 70： 1-13.

［15］CHEN X J， YANG Y M， CUI Z X， et al. Vibration Fault Diagnosis of Wind Turbines Based on Variational Mode Decomposition and Energy Entropy［J］.Energy， 2019， 174： 1100-1109.

［16］JIANG C D， MIKE M P， WANG Q， et al. Two-Dimensional QT Inversion of Complex Magnetic Resonance Tomography Data［J］.Geophysics， 2018， 83（6）： 65-75.

［17］林婷婷， 史文龍， 齊鑫， 等. 核磁共振地下水探測全波收錄系統 [J]. 吉林大學學報（工學版）， 2014， 44（4）： 1024-1030. （LIN T T， SHI W L， QI X， et al. Full-Wave Recording of MRS Ground Water Exploration System [J]. Journal of Jilin University （Engineering and Technology Edition）， 2014， 44（4）： 1024-1030.）

［18］林君， 張洋. 地面磁共振探水技術的研究現狀與展望 [J]. 儀器儀表學報，2016， 37（12）： 2657-2670. （LIN J， ZHANG Y. Research Status and Expectation of Surface Nuclear Magnetic Resonance Technique for Groundwater Detection [J]. Chinese Journal of Scientific Instrument， 2016， 37（12）： 2657-2670.）

（責任編輯： 韓 嘯）