基于自適應和小波模極大值重構的地面核磁共振信號噪聲壓制方法

孫淑琴，孟慶云，方秀成，林 君，王應吉，田寶鳳，楊 楠

（吉林大學 儀器科學與電氣工程學院，長春130061）

0 引 言

核磁共振測水方法是利用地層中豐度最高的氫原子核的順磁性產生的核磁共振（Magnetic resonance sounding，MRS）信號，通過改變激發電流的大小完成對不同深度地下水的測量［1－2］。對采集的信號經過濾波處理之后，即可反演解釋得到地下含水層深度、厚度、含水量大小、含水層平均孔隙度等信息［3－4］。

國內外專家學者提出了多種方法來抑制噪聲，主要有：儀器采集數據的疊加方法［5－6］，缺點在于疊加次數多，系統的工作效率低；探測線圈采用八字形線圈［7－8］方式，與圓形探測線圈相比，探測深度降低一半；采用同步檢波和FIR 低通濾波技術［9－10］的方法會引起有用信號的失真；此外還有滑動平均濾波法、區塊對消法，正弦對消法、陷波濾波器法［11－18］。這幾種方法主要針對電磁噪聲的濾除，由于噪聲的無規律性，以上濾波方法均存在一定的局限性。

上述濾波方法主要針對MRS信號中振幅提取的效果進行分析和處理，然而未對MRS 信號中頻率、平均衰減時間、相位等其他幾個特征參數提取情況進行分析和總結。本文研究了自適應和小波模極大值重構方法，自適應濾波方式主要是針對固定頻率噪聲的濾波處理，小波模極大值重構方法主要是針對白噪聲的濾波處理［19－22］。將兩者結合起來能夠突破傳統傅里葉變換在時域沒有任何分辨率的限制，具有良好的時域分析特性，能夠從強干擾的信號中提取有用成分，彌補了其他方法在非平穩信號處理上的不足。將本文方法與自適應濾波方法和小波模極大值重構方法進行了對比，可知用本文方法所得到的信噪比更高，信號曲線與噪聲曲線能夠得到明顯的分開，且信號和噪聲的曲線都變得更光滑，對水文地質參數的解釋更準確。

1 核磁共振測試信號及數據處理

1.1 數據處理流程

利用JLMRS儀器所采集的信號中通常混有復雜的噪聲成分，表達式如下：

式中：第一項為MRS信號表達式；第二項為工頻諧波噪聲表達式，由多個單頻的諧波分量疊加而成，其中每個單頻率諧波的參數都是獨立的；第三項為隨機噪聲，εreal（t）中包含平穩的隨機噪聲和非平穩的隨機噪聲。

儀器測試信號的數據處理過程如圖1所示。測試信號經過正交分解、噪聲濾除、提取信號的初始振幅E0，平均衰減時間T＊2，接收頻率f 和初始相位φ0 后，經過反演計算即可得到地下含水量及含水層孔隙度等信息。MRS 數據參數為E0、T＊2、φ0，相應的水文地質參數為單位體積含水量（有效孔隙度）、空隙大小（滲透性）和含水層導電性（電阻率）。MRS信號的幾個特征參數對于含水地層的解釋均有直接關系。

圖1 MRS信號數據處理流程圖Fig.1 Flowchart of MRS signal data processing

1.2 自適應濾波算法原理

當原始信號中包含頻率為ω的工頻諧波干擾時，利用自適應濾波器可以消除這種特定頻率成分的干擾，自適應濾波器結構圖［23］如圖2 所示。圖中核磁共振信號與干擾信號疊加，經采樣后給dj端，參考輸入是與單頻干擾相關的正弦信號，采樣送到d1j、d2j端，d2j是由參考信號經過90°相移得來，由兩個參考信號得到ω1j、ω2j兩個權值，即自由度，使得組合后正弦波振幅、相角與原始輸入干擾分量的振幅、相角相同，從而使輸出的ej中的ω 被抵消掉，達到濾波的目的。

圖2 自適應濾波器結構圖Fig.2 Structure chart of adaptive trap

1.3 二進小波變換

設信 號f（t）是 平 方 可 積 的，即f（t）∈L2（R），則這一信號的二進小波變換定義為［24］：

式中：j∈z；τ為平移參數；ψ（t）為母小波。

假設χ（t）為尺度函數，應用式（2）可知χ（t）、ψ（t）滿足以下尺度方程：

式中：α為低通濾波器：β為高通濾波器。

如果對式（2）的τ 進行離散化，并假設χ（t）∈L2（R）、ψ（t）∈L2（R），二進小波變換的小波系數和尺度系數為：

經過濾波器級聯結構求出各尺度上的小波系數。令αj（k）和βj（k）分別是在濾波器α（k）和β（k）每兩個系數之間添加2j－1個零后所得到的濾波器，記。二進小波變換的快速分解算法為：

快速重構算法為：

1.4 小波模極大值重構消噪的原理

將含有噪聲的信號進行多尺度小波分解，求取相應的模極大值，然后根據信號和噪聲在不同尺度因子下的不同傳播特性來區分信號和噪聲，這就是小波模極大值重構消噪的理論依據。

通過Lipschitz指數α 可以度量信號在某個時間區間或某一時刻的正規性，即刻畫信號的奇異性。奇異性也可以通過跟蹤小波變換在細尺度下的模極大值來檢測。而且Lipschitz指數和二進小波變換模極大值之間存在以下關系：如果信號的Lipschitz指數α＞0，該信號的二進小波變換模極大值將隨著尺度的增大而增大，這樣模極大值點是由需要的信號所產生；反之，α ＜0，該信號的二進小波變換模極大值將隨著尺度的增大而減小，這樣的模極大值由噪聲產生。

根據以上原理可知，隨著尺度因子的不斷變大，噪聲小波系數的模極大值點的幅度值將明顯變小。經過幾次二進小波變換之后，由噪聲產生的模極大值點的幅度值已經變得非常小，甚至很多由噪聲產生的模極大值點的值都變為零，這樣剩下的模極大值點主要是由信號產生的。一般情況下，尺度越大，信號產生的小波系數模極大值點的比例就越大，大分解尺度上的模極大值點主要是由信號產生的。所以就可以利用大尺度上的小波系數模極大值點逐步確定由信號產生的小尺度上的小波系數模極大值點，最后根據剩余的模極大值點重構信號，這樣就能夠將噪聲削弱。

1.5 小波模極大值重構消噪過程

根據小波模極大值重構消噪原理，消噪算法過程如下：

（1）對含有噪聲的信號進行二進小波變換，小波分解尺度的總個數選擇應該適中，一般小波分解尺度數為4或5，大于此值有用信號可能會丟失，小于此值不能保證產生的小波系數模極大值點的優勢程度。最后，求出小波系數的模極大值點及相對應的模極大值。

（2）設定小波系數的閾值為T，最大尺度上的小波系數模極大值如果大于閾值T，則保留相對應的小波模極大值點，否則剔除相對應的小波模極大值點，從而使最大尺度上小波系數的小波模極大值點更新。根據核磁共振實際處理，假設閾值T ＝0.1max，其中max 表示小波模極大值中的最大值。

（3）假設最大尺度為j。利用第（2）步中的小波模極大值點的位置構造U（xjn），xjn為尺度j 上的第n個小波模極大值點，之后在U（xjn）內找到尺度j－1上的小波系數模極大值點，使尺度j－1上位于U（xjn）內的模極大值點保留，其余的小波模極大值點進行清除，尺度j－1上的小波模極大值點就進行了更新。假設j＝j－1，當尺度數為2時，停止循環第（3）步。

（4）尺度j＝1或j＝2位置上相對應的極大值點進行保留，其他位置的極大值點置為零。

（5）應用各個尺度上余下的小波模極大值點進行小波系數的重構，最后利用重構的小波系數重構信號。

2 仿真模型的建立

在地下水探測技術中，核磁共振儀器采用4倍頻采樣的數字正交技術來獲得包絡信號。模擬這樣的一個理想信號并加入干擾噪聲來驗證本文研究方法的效果，將處理后的信號與原始信號作對比，以信噪比（信號強度與噪聲強度的比值）高低為評價指標，驗證本文方法的可行性效果。

定義信號的信噪比如下：

式中：Sx、Sy分別為含噪信號的x 分量和y 分量；Nx、Ny分別為噪聲的x 分量和y 分量。

2.1 理想信號仿真

核磁共振信號中的噪聲主要是白噪聲和工頻噪聲，工頻噪聲頻率一般為50Hz的整數倍，為了模擬野外的實際信號，在理想信號中疊加白噪聲和工頻噪聲。建立仿真模型，信號初始振幅E0設為400nV，初始相位φ0 設為165°，弛豫時間T＊2設為150ms，頻率f 設為2325Hz的理想衰減正弦波信號，其包絡信號和頻譜圖如圖3所示。

圖3 理想包絡信號和頻譜圖Fig.3 Ideal signal envelope and spectrum diagram

2.2 信號加入白噪聲

首先在包絡信號上加入高斯白噪聲，其幅值為150nV。此時信號的信噪比為1.43，見圖4。

圖4 白噪聲包絡信號和頻譜圖Fig.4 Signal envelope and spectrum diagram of white noise

圖5 去噪前后模極大值點對比Fig.5 Contrast before and after of modulus maxima

對加入白噪聲的MRS信號與經小波模極大值重構處理后的MRS信號進行模極大值點的對比如圖5所示。由圖可知：原始數據中噪聲大部分集中于第1層和第2層，經小波模極大值重構處理后保存下來的系數較少。隨著分解尺度的增加噪聲部分越來越少，有用信號成為主要部分，所以保留下來的系數也越來越多。小波模極大值重構算法的主要優點是根據小波變換系數模極大值的跨尺度傳播規律，可以明確區分每一分解層次上代表噪聲和信號的小波變換系數，既不會去掉微弱的有用信號，也不會保留強度較大的噪聲［25］。

自適應濾波主要針對固定頻率噪聲的濾除，小波模極大值重構主要針對白噪聲的濾除，為了反映基于自適應和小波模極大值重構濾波方法的優越性，下面給出了幾種濾波方式的對比。經過不同濾波方式所得到的包絡信號和頻譜圖如圖6所示。圖6（b）和（c）所得到的曲線變得比較光滑，信號與噪聲明顯地被分開。在3種濾波方法中，經自適應和小波模極大值重構所得到的信噪比最高，所以經此種濾波方式所得到的信號是最可信的。

圖6 處理之后的包絡信號和頻譜圖Fig.6 Envelope signal and the spectrum diagram after filtering

對加入白噪聲的信號進行濾波，比較濾波前后情況，白噪聲在整個頻域范圍內均勻分布，無固定頻率值。比較結果如表1所示。

由表1可以看出：經過自適應和小波模極大值重構處理的數據與另外兩種方法對比，所得到的特征參數最接近于理想信號且信噪比最高，信噪比可達到2.61。經過自適應和小波模極大值重構處理后的MRS特征參數的擬合誤差E0%達到0.82%，達到6.3%，φ0%達到1.75%，頻譜幅度衰減率達到0.12%，滿足實際應用的要求［14－15］。MRS數據處理之后，所得到的特征參數越接近于理想信號說明處理的結果越準確，反演解釋結果的可信度越高。同時由圖6可得，經自適應濾波和小波模極大值重構方法所得到的信號曲線和噪聲曲線變得比較光滑且信號與噪聲明顯分開，獲得了可信的包絡信號。

2.3 信號加入工頻噪聲

理想信號中加入工頻諧波干擾，其中頻率f1＝2350Hz，幅值為80nV。信噪比為2.88，如圖7所示。

對加入工頻噪聲的MRS 信號進行處理前、后模極大值點的對比如圖8所示。與加入白噪聲的對比結果類似，隨著分解尺度的增加噪聲部分越來越少，有用信號成為主要部分，所以保留下來的系數也越來越多。

為了反映基于自適應和小波模極大值重構濾波方法的優越性，下面給出了幾種經過不同濾波方式得到的包絡信號和頻譜圖，如圖9所示。由圖可以看出：經自適應和小波模極大值重構濾波后所得到的信號曲線最光滑，且噪聲壓制得最低，信號與噪聲分開的效果最好。

對加入噪聲的信號進行濾波，比較濾波前、后情況。核磁共振含噪信號經過不同濾波方法得到特征參數的擬合誤差與信噪比的大小以及頻譜幅度衰減率的比較結果見表2。

表1 數據經不同濾波方式后的對比（加入白噪聲）Table 1 Contrast data with different filtering way（white noise）

圖7 工頻噪聲包絡信號和頻譜圖Fig.7 Signal envelope and spectrum diagram of industrial frequency noise

圖8 去噪前后模極大值點對比Fig.8 Contrast before and after of modulus maxima

由表2可以看出，經過基于自適應的小波模極大值重構處理的數據與另外兩種方法相對比，所得到的特征參數最接近于理想信號，且信噪比最高。MRS特征參數（振幅E0，弛豫時間，相位φ0）的擬合誤差控制在了3%以內，且頻譜幅度衰減率達到0.68%。由加入白噪聲所處理的結果分析可知，MRS數據處理之后，所得到的特征參數接近于理想信號所對應的特征參數，說明處理的結果準確，反演解釋結果的可信度高。

圖9 處理之后的包絡信號和頻譜圖Fig.9 Envelope signal and the spectrum diagram after filtering

表2 數據經不同濾波方式后的對比（加入工頻噪聲）Table 2 Contrast data with different filtering way（power noise）

結合處理后的信號曲線圖與表2可知：經自適應和小波模極大值重構處理之后所得到的信號曲線的平滑度最好，可以很好地將噪聲與信號分離開來，信噪比達到最高且最可信。

3 實測數據處理

2005年4月，核磁組在內蒙古烏拉特地區進行了野外實驗，對野外采集數據應用不同的濾波方法進行數據處理。首先對每個測點的單個脈沖矩進行濾波處理，再利用處理后的數據得出反演結果。

（1）測點WLTQ100。當地的拉莫爾頻率為2353Hz，對原始數據進行處理，然后提取MRS信號的特征參數E0、T＊2、φ0、f。由表1可知MRS信號的特征參數所對應的水文地質參數。圖10為MRS信號的不同特征參數所對應的曲線。由圖可得出，經自適應和小波模極大值重構濾波后所得到的曲線的平滑度和穩定性最好。

圖10 不同濾波方式下的MRS信號特征參數曲線Fig.10 MRS signal characteristic parameter curves after different filter methods

圖11為經過不同濾波方式反演所得到的含水量柱狀圖。測點 WLTQ100 距離打井點XF277很近，為2.03km，認為兩者的水文地質情況是相似的。XF277的水文地質情況如圖12所示。由圖可知：XF277打井深度為1.575～50m，含水層厚度為48.425m。飽和含水層從5m 開始，5～6 m 為粗砂，6～17.9 m 為中砂，17.9～34.9m 為粗砂，34.9～50m 為細砂。由圖11、圖12及與含水層的巖性關系分析可得，經自適應和小波模極大值重構反演所得到的含水量柱狀圖與鉆探結果基本相符，更接近于實際情況。

圖11 不同濾波后反演所得到的含水量柱狀圖Fig.11 MRS inversion results with the different methods

圖12 XF277水文地質情況Fig.12 XF277hydrological geology

（2）測點WLTQ94。當地的拉莫爾頻率為2354Hz，與WLTQ100處理過程相同。圖13為MRS信號相應的特征參數對應的曲線。從圖13中可知：經自適應和小波模極大值重構濾波后所得到的特征曲線的平滑度提高的程度最大，穩定性最好。

圖13 不同濾波方式下的MRS信號特征參數曲線Fig.13 MRS signal characteristic parameter curves after different filter methods

圖14 為經過不同濾波方式反演所得到的含水量柱狀圖。測點WLTQ94 距離打井點HL11很近，為3.17km，同樣兩者的水文地質情況是相似的。HL11的水文地質情況如圖15 所示。由圖可知：HL11打井深度為1.5～42m，含水層厚度為40.5m，出水層從1.5m 開始，1.5～5m 為粗砂，5～21.7 m 為中砂，21.7～42 m 為粗砂。經圖14、圖15及T＊2與含水層的巖性關系的分析得出，經自適應和小波模極大值重構反演得到的含水量柱狀圖更接近于鉆井的實際情況。

圖14 不同濾波后反演所得到的的含水量柱狀圖Fig.14 MRS inversion results with the different methods

圖15 HL11水文地質情況Fig.15 HL11hydrological geology

4 結 論

（1）MRS特征參數（振幅E0，弛豫時間相位φ0）的擬合誤差控制在了3%以內。由處理的結果分析可知，MRS數據處理之后，所得到的特征參數接近于理想信號所對應的特征參數，說明處理的結果準確，反演解釋結果的可信度高。

（2）實測數據處理表明，在地下水含量未知的條件下對于不同強度的噪聲干擾基于自適應和小波模極大值重構算法獲得了較好的去噪效果。

（3）通過理想信號加入兩種噪聲，經過不同濾波方法的比較，得出了不同噪聲情況下的消噪策略，實測信號也驗證了上述策略的適用性，為野外試驗的數據處理提供了指導，縮短了尋找最佳方法所用時間，提高了效率，對以后消噪方法的研究給出了參考。

［1］Legchenko A，Baltassat J M，Bobachev A，et al.Magnetic resonance sounding applied to aquifer characterization［J］.Ground Water，2004，42（3）：363－373.

［2］Legchenko A，Baltassat，J M，Beauce A，et al.Nuclear magnetic resonance as a geophysical tool for hydrogeologists［J］.Journal of Applied Geophysics，2002，50（1）：21－46.

［3］Müller－Petke M，Hiller T，Herrmann R，et al.Reliability and limitations of surface NMR assessed by comparison to borehole NMR［J］.Near Surface Geophysics，2011，9（2）：123－134.

［4］Vouillamoz J M，Legchenko A，Nandagiri L.Characterizing aquifers when using magnetic resonance sounding in a heterogeneous geomagnetic field［J］.Near Surface Geophysics，2011，9（2）：135－144.

［5］Juan Plata，Felix Rubio.MRS experiments in a noisy area of a detrital aquifer in the south of Spain［J］.Journal of Applied Geophysics，2002，50（1）：83－94.

［6］Jiang C D，Lin J，Duan Q M，et al.Statistical stacking and adaptive notch filter to remove high－level electromagnetic noise from MRS measurements［J］.Near Surface Geophysics，2011，9（5）：459－468.

［7］Konstantinos Chalikakis，Mette Ryom Nielsen，Anatoly Legchenko.MRS applicability for a study of glacial sedimentary aquifers in Central Jutland，Denmark［J］.Journal of Applied Geophysics，2008，66（3）：176－187.

［8］Trushkin D V，Shushakov O A，Legchenko A V.The potential of a noise－reducing antenna for surface NMR groundwater surveys in the earth＇s magnetic field［J］.Geophysical Prospecting，1994，42（8）：855－862.

［9］Anatoly Legchenko，Pierre Valla.A review of the basic principles for proton magnetic resonance sounding measurements［J］.Journal of Applied Geophysics，2002，50（1）：3－19.

［10］Legchenko A，Valla P.Processing of surface proton magnetic resonance signals using non－linear fitting［J］.Journal of Applied Geophysics，1998，39（2）：77－83.

［11］Legchenko A，Valla P.Removal of power－line harmonics from proton magnetic resonance measurements［J］.Journal of Applied Geophysics，2003，53（2）：103－120.

［12］Legchenko A.MRS measurements and inversion in presence of EM noise［J］.Boletín Geológico y Minero，2007，118（3）：489－508.

［13］Walsh D O.Multi－channel surface NMR instrumentation and software for 1D／2Dgroundwater investigations［J］.Journal of Applied Geophysics，2008，66（3）：140－150.

［14］田寶鳳，段清明.核磁共振信號工頻諧波的自適應濾除方法［J］.吉林大學學報：信息科學版，2009，27（3）：223－228.Tian Bao－feng，Duan Qing－ming.Removal method of industrial frequency harmonics in nuclear magnetic resonance signal based on adaptive filter［J］.Journal of Jilin University（Information Science Edition），2009，27（3）：223－228.

［15］田寶鳳，林君，段清明，等.基于參考線圈和變步長自適應的磁共振信號噪聲壓制方法［J］.地球物理學報，2012，55（7）：2462－2472.Tian Bao－feng，Lin Jun，Duan Qing－ming，et al.Variable step adaptive noise cancellation algorithm for magnetic resonance sounding signal with a reference coil［J］.Chinese Journal of Geophysics，2012，55（7）：2462－2472.

［16］曾亮，李振宇，王鵬.小波分析在提高核磁共振找水信號信噪比中的應用探討［J］.CT 理論與應用研究，2007，15（2）：1－5.Zeng Liang，Li Zhen－yu，Wang Peng.Improveing S／N ratio of MRS signal of detecting ground water with wavelet analysis［J］.CT Theory and Applications，2007，15（2）：1－5.

［17］Strehl S，Rommel I，Hertrich M，et al.New strategies for filtering and fitting of MRS signals［C］∥Proc Proceedings 3rd International MRS Workshop，Madrid，Spain，2006：65－68.

［18］Strehl S.Development of strategies for improved filtering and fitting of SNMR－Signals［D］.Berlin：Department of Applied Geophysics Diplomarbeit，Institute of Applied Geosciences，Technical University of Berlin，2006.

［19］Abdullah Al Jumah.Denoising of an image using discrete stationary wavelet transform and various thresholding techniques［J］.Journal of Signal and Information Processing，2013，4（1）：33－41.

［20］Nair S S，Joseph K P.Wavelet based electroretinographic signal analysis for diagnosis［J］.Biomedical Signal Processing and Control，2014（9）：37－44.

［21］張金榜，孫藝笑，王潤典，等.改進的閾值函數去噪算法［J］.電子科技，2014，27（2）：17－20.Zhang Jin－bang，Sun Yi－xiao，Wang Run－dian，et al.An improved threshold denoising algorithm［J］.Electronic Science and Technology，2014，27（2）：17－20.

［22］李文剛，張振勇，王艷波，等.小波分析在礦井高密度電法數據處理中的應用［J］.中國煤炭地質，2011，23（6）：52－55.Li Wen－gang，Zhang Zhen－yong，Wang Yan－bo，et al.Application of wavelet analysis in coalmine highdensity electric prospecting data processing［J］.Coal Ceology of China，2011，23（6）：52－55.

［23］Sun S，Yang N，Lin J，et al.Adaptive analysis of filter methods for magnetic resonance sounding［C］∥In Digital Manufacturing and Automation（ICDMA），IEEE，2013：106－111.

［24］秦毅，王家序，毛永芳.基于軟閾值和小波模極大值重構的信號降噪［J］.振動測試與診斷，2011，31（5）：543－547.Qin Yi，Wang Jia－xu，Mao Yong－fang.Based on soft threshold and wavelet modulus maxima reconstruction signal noise reduction［J］.Journal of Vibration，Measurement ＆Diagnosi，2011，31（5）：543－547.

［25］何永紅，文鴻雁，靳鵬偉.基于小波模極大值改進算法的變形模型研究［J］.測繪科學，2007，32（4）：18－20.He Yong－hong，Wen Hong－yan，Jin Peng－wei.Deformation of the improved algorithm based on wavelet modulus maxima model research［J］.Science of Surveying and Mapping，2007，32（4）：18－20.