王 琦，杜海龍，張 蒙，韋 健

（1.吉林大學(xué)通信工程學(xué)院，長(zhǎng)春 130026；2.中國(guó)能源建設(shè)集團(tuán)遼寧電力勘測(cè)設(shè)計(jì)院有限公司，沈陽(yáng) 110000）

0 引言

磁共振探測(cè)（Magnetic Resonance Sounding，MRS）技術(shù)以其直接、定量、準(zhǔn)確和高效等優(yōu)點(diǎn)，已應(yīng)用于地下水探測(cè)、區(qū)域水資源調(diào)查和地質(zhì)災(zāi)害預(yù)警等領(lǐng)域［1-3］。然而，由于城市和村莊等環(huán)境中存在大量的隨機(jī)噪聲，導(dǎo)致數(shù)據(jù)的信噪比較低，難以為后續(xù)的水文解釋提供可靠結(jié)果［4-5］。因此本文針對(duì)隨機(jī)噪聲干擾下磁共振信號(hào)提取方法展開研究。

時(shí)頻分析方法將一維時(shí)間域信號(hào)拓展到二維時(shí)頻域，實(shí)現(xiàn)信號(hào)的時(shí)頻分布信息的定位，反映信號(hào)的能量分布狀態(tài)，利用時(shí)間和頻率的聯(lián)合函數(shù)來表征信號(hào)［6-7］。其中，短時(shí)傅里葉變換（Short-time Fourier Transform，STFT）作為一種基本的時(shí)頻分析方法，是在傅里葉變換的基礎(chǔ)上，利用窗函數(shù)截?cái)喾瞧椒€(wěn)信號(hào)，再對(duì)該信號(hào)做傅里葉變換的方法［8］。但是，該方法無法同時(shí)滿足時(shí)間和頻率的高分辨率要求。因此，于剛等［9-10］提出了同步提取變換（Synchro-extracting Transform，SET）方法，實(shí)現(xiàn)了基于SET 的挖掘機(jī)振聲信號(hào)時(shí)頻分析。同步提取變換是在STFT 的基礎(chǔ)上，通過構(gòu)造同步提取算子，僅利用瞬時(shí)頻率位置的時(shí)頻系數(shù)生成新的時(shí)頻譜，達(dá)到提高能量集中度的目的。而且，SET信號(hào)恢復(fù)需要較少的參數(shù)，因此可以更方便地實(shí)現(xiàn)信號(hào)重建［11-12］。因此，本文提出將SET 應(yīng)用于地面磁共振信號(hào)的提取。

本文介紹了MRS的基本特征，基于SET 原理，使MRS數(shù)據(jù)的時(shí)頻系數(shù)集中在拉莫爾頻率位置，并利用同步提取算子（Synchro-extracting Operator，SEO）得到MRS數(shù)據(jù)的SET 時(shí)頻譜。進(jìn)一步，基于脊重建原理，實(shí)現(xiàn)MRS信號(hào)的重構(gòu)。最后，本文研究了不同噪聲水平及不同窗寬對(duì)MRS信號(hào)提取結(jié)果的影響。

1 磁共振信號(hào)提取方法

1.1 MRS信號(hào)特征

地面核磁共振通過向地面上的線圈回路中通入Larmor頻率的交變電流，激發(fā)地下水中的氫質(zhì)子，使其旋轉(zhuǎn)軸發(fā)生偏轉(zhuǎn)，即發(fā)生核磁共振現(xiàn)象。當(dāng)交變電流停止后，氫質(zhì)子自旋逐漸恢復(fù)到地磁場(chǎng)方向，稱為弛豫過程。此時(shí)，地面上的線圈回路接收到MRS信號(hào)為

式中：e0表示FID信號(hào)的初始振幅；表示橫向弛豫時(shí)間；ω ＝2πfL，fL表示拉莫爾頻率；φ0表示初始相位。

MRS信號(hào)的傅里葉變換為

通常，線圈接收到的FID 信號(hào)的幅度約為500 nV，而環(huán)境噪聲可達(dá)100 nV 以上。因此，采集到的MRS數(shù)據(jù)的信噪比較低，導(dǎo)致信號(hào)提取準(zhǔn)確率低。

1.2 磁共振信號(hào)的同步提取變換

同步提取變換是作為短時(shí)傅里葉變換的后處理過程，因此，首先對(duì)MRS數(shù)據(jù)進(jìn)行STFT變換［13-14］：

式中：g（u-t）為窗函數(shù)，由于MRS信號(hào)隨時(shí)間呈指數(shù)衰減趨勢(shì)，本文利用高斯窗g（t）＝。設(shè)MRS 信號(hào)的頻率fL＝1.81 kHz，e0＝50 nV，＝0.5 s，φ0＝π／4 rad，窗函數(shù)中a＝0.25。得到時(shí)頻譜如圖1（a）所示，可以看出MRS 能量主要分布在1.79～1.83 kHz之間，并在1.81 kHz處達(dá)到最大值，STFT不能精確地表示信號(hào)的時(shí)頻特征。

圖1 MRS信號(hào)的時(shí)頻譜結(jié)果

對(duì)式（3）增加附加相移eiωt，則

根據(jù)帕斯瓦爾定理，式（3）可以寫為［15］

式中，V（ξ）表示V（u）的傅里葉變換。

將式（2）代入式（5），得到FID信號(hào)的STFT為

可知，F(xiàn)ID的能量集中在ω0附近，并在ω0處取得最大值。為了提高FID 的時(shí)頻分辨率，使能量集中在ω0處，引入同步提取算子SEO，得到FID信號(hào)的同步提取變換為［9］：

1.3 磁共振信號(hào)的重構(gòu)

SET變換可近似為時(shí)頻譜的脊提?。?6］，因此可利用脊重建方法實(shí)現(xiàn)MRS信號(hào)的重構(gòu)。由于MRS 時(shí)頻系數(shù)的瞬時(shí)頻率為

則SET變換可以寫為

即當(dāng)ω等于瞬時(shí)頻率時(shí)，

因此，重構(gòu)得到的MRS信號(hào)為

2 實(shí)驗(yàn)結(jié)果與分析

2.1 實(shí)驗(yàn)步驟

為了驗(yàn)證基于SET方法提取MRS信號(hào)的有效性，本文仿真構(gòu)造一組MRS 信號(hào)：fL＝1.81 kHz，e0＝100 nV，＝0.4 s，φ0＝π／3 rad，并加入10 nV的高斯噪聲，采樣頻率為10 kHz，采樣時(shí)間為1 s，如圖2 中黑色曲線所示，對(duì)其進(jìn)行SET 變換，時(shí)頻譜如圖3 所示?？梢钥闯?，時(shí)頻譜能量主要集中在1.81 kHz，受隨機(jī)噪聲影響，時(shí)頻譜在0.7～1 s間存在微小波動(dòng)。利用式（11）得到MRS信號(hào)重構(gòu)結(jié)果如圖2 中藍(lán)色曲線所示，可以看出重構(gòu)結(jié)果包絡(luò)光滑，不存在明顯的隨機(jī)噪聲。進(jìn)一步，將重構(gòu)后的MRS信號(hào)與仿真的無噪信號(hào)相減，得到殘差如圖2 中灰色曲線所示，計(jì)算得到平均誤差為0.40 nV。因此可以得出，SET 方法適用于MRS信號(hào)的提取。

圖2 基于SET的MRS信號(hào)提取結(jié)果

圖3 MRS數(shù)據(jù)SET的時(shí)頻譜

2.2 窗函數(shù)對(duì)MRS信號(hào)提取結(jié)果的影響

研究了不同窗函數(shù)對(duì)MRS信號(hào)提取結(jié)果的影響，給出了高斯窗函數(shù)中a以步長(zhǎng)0.1 從0.1 變化至0.5時(shí)對(duì)應(yīng)的MRS 信號(hào)提取的統(tǒng)計(jì)結(jié)果（重復(fù)100 次實(shí)驗(yàn)），如圖4 所示。圖4（a）為5 組不同的高斯窗函數(shù)，圖4（b）為不同高斯窗函數(shù)在3 組不同噪聲水平數(shù)據(jù)下對(duì)應(yīng)的MRS信號(hào)的統(tǒng)計(jì)平均誤差。由圖4（b）可以看出，隨著噪聲水平的增加，基于SET 的MRS 信號(hào)提取精度降低。對(duì)于同一噪聲水平情況，當(dāng)噪聲水平為5 nV時(shí)，高斯窗函數(shù)a＝0.3 時(shí)，平均誤差取得最小值，0.24 nV；當(dāng)噪聲水平為10 nV時(shí)，高斯窗函數(shù)a＝0.4 時(shí)，平均誤差取得最小值，0.42 nV；當(dāng)噪聲水平為10 nV時(shí)，高斯窗函數(shù)a＝0.4 時(shí)，平均誤差取得最小值，0.57 nV?？梢缘贸?，隨著噪聲水平的增加，適當(dāng)增大窗寬，可以得到更加理想的MRS信號(hào)提取結(jié)果。

圖4 MRS信號(hào)提取結(jié)果隨高斯窗函數(shù)的變化規(guī)律

3 結(jié)語

本文針對(duì)磁共振信號(hào)的隨機(jī)噪聲干擾，提出了基于同步提取變換（SET）的MRS 信號(hào)提取方法。該方法在STFT的基礎(chǔ)上，通過構(gòu)造同步提取算子，達(dá)到提高時(shí)頻譜系數(shù)能量集中度的目的。進(jìn)一步，通過實(shí)驗(yàn)結(jié)果分析，得出隨著環(huán)境噪聲水平的增加，適當(dāng)增大SET窗函數(shù)寬度，可以提高M(jìn)RS 信號(hào)提取結(jié)果精度。本文內(nèi)容可以作為“信號(hào)與系統(tǒng)”“現(xiàn)代信號(hào)處理”等課程的提升性實(shí)驗(yàn)內(nèi)容。