基于局部均值分解的地面核磁共振信號包絡提取方法

田寶鳳，孫士聰，劉隆昌，蔣川東

1.地球信息探測儀器教育部重點實驗室(吉林大學)，長春 130026 2.吉林大學儀器科學與電氣工程學院，長春 130026

0 引言

地面核磁共振(surface nuclear magnetic resonance, SNMR)又稱為磁共振測深(magnetic resonance sounding, MRS)，是一種先進的能夠定量進行地下水測量的地球物理勘探方法。與其他的地球物理方法相比，SNMR是一種定量、直接的找水方法，通過對地下水中氫質子的響應來實現地下水探測，無需鉆孔就能夠反演解釋地下水的含量、位置以及儲水介質的孔隙大小等水文地質信息[1]。因此，近年來，該方法在水資源探測、評價以及地質災害隱患監測預警等方面均獲得了廣泛應用[2-3]。但在實際的信號采集與包絡提取過程中，噪聲環境較為復雜，往往會淹沒納伏量級的SNMR信號，導致目標信號特征無法辨識，嚴重制約了SNMR信號包絡的提取[4]。

針對SNMR有效信號檢測的問題，通常的技術手段是首先進行噪聲壓制，然后再通過數據擬合方法提取包絡。普遍采用尖峰建模方法去除尖峰噪聲[5]，采用基于最小二乘法理論的工頻諧波建模等方法消除工頻諧波噪聲[6-7]，采用雙奇異值分解等方法同時抑制工頻諧波和隨機噪聲干擾[8]，實現參數提取。此外，田寶鳳[9]等提出了以獨立成分分析為基礎的快速固定點算法實現信噪分離，王琦等[10]對隨機噪聲背景下多弛豫SNMR信號提出了基于稀疏表示的參數提取方法。盡管如此，由于噪聲環境的高復雜性、不確定性與隨機性，SNMR信號的包絡提取問題依然是地面磁共振探測領域的前沿核心技術難題[11]。

因此，根據信號特征直接獲得SNMR信號的包絡提取方法吸引了諸多學者的目光。Ghanati等[12]基于經驗模態分解(empirical mode decomposition，EMD)以及完全集合經驗模態分解方法對有效SNMR信號衰減趨勢進行提取，但只有信噪比較高時才能獲得較好的包絡提取效果。Liu Lichao等[13]通過對SNMR信號使用譜分析和滑動窗口的方式提取信號包絡。上述方法雖然在一定前提條件下獲得了一定效果，但依然存在一定的局限性。

為實現SNMR信號與噪聲的有效分離，進行多種工況下信號包絡的高精度提取，本文提出了基于局部均值分解(local mean decomposition, LMD)的SNMR信號包絡提取方法[14]。LMD算法能夠將信號分解成一系列乘積(product function，PF)分量，PF分量由包絡信號與純調頻信號相乘獲得，因而瞬時頻率是正的、連續的、具有物理意義的[15]。相較于常用的EMD算法，LMD算法采用滑動平均法來處理目標信號，避免了EMD算法采用三次樣條函數形成上下包絡時產生的過包絡和欠包絡現象[16]。此外，LMD算法比EMD算法的端點效應小且擴散較慢[17]。本文首先采用LMD算法將原始含噪SNMR信號逐級分離，得到一系列的PF分量；其次通過PF分量逐一判定去除其他高頻成分；然后通過余量提取SNMR信號的包絡信息，進一步采用非線性擬合方法提取關鍵參數，進行誤差分析[18]；再通過仿真實驗，論證方法的提取效果與可行性；最后，通過野外實測數據進一步證明所提LMD算法的有效性和實用性。

1 方法原理

1.1 局部均值分解

LMD算法能夠自適應地將非平穩、非線性的信號分解為一系列包含信號特征的PF分量，凸顯信號不同時間尺度的波動特征，保留原始信號的幅頻特性。

首先，對目標信號f(t)進行處理，得到第一個PF分量Z1(t)；接著，從原始信號f(t)中分離出Z1(t)，得到余量q1(t)，并將其作為新的待分解信號，循環分解；然后直至經過k次循環后，分離出全部的PF分量，得到殘余分量qk(t)為一個單調函數，LMD算法分解結束。

(1)

最終，原始信號f(t)分解得到k個PF分量和殘余分量：

(2)

式中，Zi(t)為信號的第i個PF分量。

1.2 PF分量計算

PF分量是包絡信號與純調頻函數的乘積，這些乘積函數構成的有限集合中元素的總和就是復雜信號的時頻分布，每一個元素稱為一個PF分量。PF分量的本質是一個幅度、頻率調制信號，它的瞬時頻率可以用純頻率調制信號直接計算得出，反映了每個波動特征對應的物理意義。而包絡信號就是PF分量的幅度調制函數。

對于頻率調制函數和幅度調制函數的計算，則通過以下過程實現。

首先，需要依次比較信號中某一采樣點與跟它相鄰的兩個采樣點的大小關系，得到信號包絡局部極值點ni，計算相鄰兩個極值點(包含極大值點和極小值點)的算術平均值mi為

(3)

計算相鄰兩個極值點的局部包絡函數值ai為

(4)

再連接mi序列得到一條折線，對ai采用相同方式得到另一條折線。

隨后，采用滑動平均法分別對2條折線進行平滑處理，計算出信號Z1(t)第一次循環時的局部均值函數和包絡估計函數，記為m11(t)與a11(t)。從信號f(t)中分離出m11(t)得到余量h11(t)為

h11(t)=f(t)-m11(t)。

(5)

利用a11(t)對h11(t)進行解調，計算得頻率調制信號為

(6)

式中：s11(t)為信號Z1(t)第一次循環下的頻率調制信號。

接著，以s11(t)作為新的待分解信號，進行多次循環分解。當分解了n次后，得到n個頻率調制函數與n個包絡估計函數。Z1(t)的分解計算過程如下：

(7)

(8)

當Z1(t)的第n個頻率調制信號s1n(t)的包絡估計函數a1n+1(t)滿足式(9)條件時停止迭代，此時，s1n(t)為Z1(t)的純調頻信號。

1-Δ≤a1n+1(t)≤1+Δ。

(9)

(10)

Z1(t)=s1n(t)a1(t)。

(11)

最后，從信號中分離出得到的PF分量，對剩余信號重復上述計算過程，循環迭代，最終得到k個PF分量和單調的殘余分量qk(t)。

1.3 SNMR信號包絡提取流程

在分解SNMR信號的過程中，信號頻率與環境中摻雜噪聲的頻率往往十分接近。為此，首先對采集到SNMR信號采用Hilbert變換，轉換成復包絡信號，提取信號的實部與虛部；再對信號的實部與虛部分別進行局部均值分解處理，得到各自的PF分量與余量。分解后信號f(t)可以表示為

(12)

式中：Re[·]表示實部數據；Im[·]表示虛部數據。

根據SNMR信號頻率的特征，分離出原始信號中的高頻成分，可以得到SNMR信號的包絡曲線：

x(t)=Re[qk(t)]+Im[qk(t)]·i。

(13)

對得到的包絡曲線進行非線性擬合提取特征參數，以此還原出理想的SNMR信號。

我國國庫集中支付制度起步較晚，從2001年開始確立，緊接著在幾所本科高校實行改革試點，并在幾年內在全國高校進行推廣。發展到今天，國庫集中支付制度已經得到了很多高校的認可，各項工作的開展正有條不紊的進行，對高校財務管理系統產生了較大影響。筆者根據實際工作和學習情況，對國庫集中支付對高校財務管理體系的構建進行研究，希望促進我國高校財務管理體系的健康發展。

2 仿真實驗與結果分析

針對采集SNMR信號時噪聲環境較為復雜、隨機性高的問題，本文基于LMD算法對仿真SNMR信號進行包絡提取，并量化提取結果的信噪比、初始振幅擬合誤差和平均橫向弛豫時間擬合誤差，來評估方法的有效性與優越性。信噪比RSN計算方式為

(14)

式中：fE為理想SNMR信號；fx為含噪信號；PE為理想SNMR信號的功率；PN為噪聲的功率。

對提取到的理想SNMR信號趨勢項進行非線性擬合，得特征參數的初始振幅和平均橫向弛豫時間的擬合誤差分別為：

(15)

(16)

式中：ΔE0為初始振幅擬合誤差；E0′為所提取包絡結果的初始振幅；E0為理想SNMR信號的初始振幅；ΔT2*為平均橫向弛豫時間擬合誤差；T2*′為所提取包絡結果的平均橫向弛豫時間；T2*為平均橫向弛豫時間。

2.1 迭代精度設置

在對信號進行局部均值分解的過程中，內層循環通過解調源數據與局部均值函數之差得到調頻信號，分析調頻信號的包絡是否為1來判斷該信號是否為純調頻信號，進而決定是否跳出此次循環進行下一個PF分量的求取。但在實際計算過程中，為了避免迭代次數過多，造成運算量增大、程序運行時間過長的問題，采取式(9)進行判斷，即純調頻信號包絡信號的值域區間為[1-Δ，1+Δ]。Δ取值為10-7。

同時，在外層循環中，還需要判斷原始含噪信號分離出所有PF分量后的余量qk(t)是否為一單調函數，決定是否繼續迭代分解下一個PF分量。當qk(t)的極值點數目小于2時，可認為分解過程結束，原始含噪信號表示為所有PF分量與余量之和。

2.2 基于LMD的SNMR信號包絡提取結果

理想的全波核磁共振信號表達式為

(17)

式中：fL為采集地點的拉莫爾頻率；φ0為初始相位。

利用本文方法對目標信號進行處理。首先，對仿真含噪信號進行Hilbert變換，將得到的復包絡信號分解為實部與虛部信號(圖2)。然后，采用LMD算法，分別對實部信號和虛部信號進行分解，分別得到其各階PF分量(圖3)。以實部和虛部信號的最后2個PF分量作為特征分量合成目標信號包絡(圖4)。在前0.05 s內由于端點效應導致包絡信號失真，但是在后續時間里包絡信號與理想SNMR信號保持一致的趨勢，包絡信號在理想包絡線處上下波動，經過算法處理后提取出的信號損失與失真程度很小。

a. 幅值-時間圖；b. 幅值-頻率圖。

圖2 含噪SNMR包絡信號實部與虛部

對LMD算法處理結果采用最小二乘法進行非線性擬合，目標函數公式為

(18)

式中，fobs為 LMD算法處理后提取的包絡信號。圖1中也給出了LMD算法提取信號的擬合結果。由圖1可見，擬合信號SNMR呈現明顯的e指數衰減趨勢，與理想SNMR信號重合程度較高，能量損失較小。局部均值分解后提取的SNMR信號初始振幅擬合誤差為1.52%，平均橫向弛豫時間擬合誤差為1.93%，信噪比提升了29.0 dB。

為了進一步討論LMD算法的適用范圍，分析算法在不同環境條件下的處理效果，開展了不同工頻諧波幅值、不同噪聲幅度、不同平均橫向弛豫時間的算法實驗，結果如表1、表2、表3所示。用LMD算法處理SNMR信號，當信號初始信噪比大于-18.0 dB時，初始振幅和平均橫向弛豫時間的擬合誤差分別控制在±4.17%和±5.63%范圍內，信噪比提升30.3～37.2 dB;當初始信噪比低于-18.0 dB時，參數擬合誤差偏大，如當初始信噪比為-21.4 dB時，初始振幅的擬合誤差可達15.02%，平均橫向弛豫時間的擬合誤差達到17.41%。

2.3 方法對比

為了更加全面、具象地說明算法的處理效果，通過仿真實驗將LMD算法與經典的EMD算法和工頻諧波建模方法進行對比，結果如圖5所示。從圖5中可以看出，對于3種方法的處理結果，LMD算法提取的SNMR信號包絡可以更好地逼近理想SNMR信號，而且受到端點效應和模態混疊的影響更小。

經過3種方法處理前所得后的參數提取誤差如

a. 實部；b. 虛部。

圖4 LMD算法提取SNMR信號包絡

表4所示，相比之下，LMD算法的初始振幅擬合誤差、平均橫向弛豫時間擬合誤差均為最小，信噪比提升程度最大，證明了LMD算法的處理效果更加準確，提取的包絡與理想SNMR信號包絡趨勢最為接近。

3 實測數據處理結果

為驗證LMD算法對于實際復雜噪聲環境中SNMR信號包絡的提取效果，于長春市文化廣場實地采集信號源產生的SNMR信號和環境噪聲的合成數據開展算法實測數據實驗驗證。采用吉林大學自主研發全波磁共振探測儀器系統采集了4組不同參數的實測數據，采集線圈尺寸為25 m×25 m，采樣率為25 kHz，采樣時長為0.5 s，數據1、數據2與數據4的拉莫爾頻率為2 355 Hz，數據3的拉莫爾頻率為2 360 Hz。從圖6中可以看出，含噪信號中都存在較大工頻和隨機噪聲，且在拉莫爾頻率點附近有噪聲的密集分布。經過LMD算法處理后，由擬合后SNMR信號與含噪SNMR信號的對比可以發現，本文方法可以有效去除噪聲，提取目標信號包

表1 不同工頻諧波幅值下LMD算法處理結果的信噪比與參數擬合誤差

表2 不同隨機白噪聲幅值下LMD算法處理結果的信噪比與參數擬合誤差

表3 不同平均橫向弛豫時間下LMD算法處理結果的信噪比與參數擬合誤差

圖5 LMD算法、EMD算法和工頻諧波建模方法處理結果對比

表4 LMD、EMD與工頻諧波建模參數擬合誤差與信噪比提升程度

絡。4組數據處理后平均橫向弛豫時間擬合誤差分別為1.53%、4.80%、3.65%、2.30%，可見，基于LMD的SNMR信號包絡提取方法有著很好應用性與精確性，可以有效地從復雜噪聲中提取目標的SNMR信號，反映水體特征。

a. 數據1；b. 數據2；c. 數據3；d. 數據4。

4 結論

1) 基于LMD算法實現了對不同信噪比下含噪信號復包絡的實部與虛部的分解，基于余量組合實現了對SNMR信號包絡的有效提取。算法仿真結果表明，該方法能夠有效去除隨機噪聲、工頻諧波等環境噪聲，初始振幅的擬合誤差在±4.17%之內，平均橫向弛豫時間的擬合誤差在±5.63%之內，信噪比提高了30.3～37.2 dB。

2) 對比分析了LMD算法與經典EMD包絡提取算法和諧波建模方法的性能，表明LMD算法有著較好的優越性與準確程度，提取的包絡與理想信號包絡更加吻合。實測數據處理結果進一步驗證了方法的有效性。該方法可以顯著提升參數提取精度，降低運算時間，減少環境噪聲對勘查結果的影響。 對于拓展磁共振探測儀器的應用范圍具有借鑒意義。