郜澤霖, 魏 晉, 蔣川東, 刁 庶

(吉林大學(xué) 儀器科學(xué)與電氣工程學(xué)院, 長春 130061)

0 引 言

地面磁共振技術(shù)(MRS: Magnetic Resonance Sounding)以其無損、 定量、 直接等優(yōu)勢, 廣泛應(yīng)用于地下水探測和水文地質(zhì)調(diào)查等領(lǐng)域[1]。在測量過程中, 獲得的MRS信號十分微弱, 易受到各種環(huán)境噪聲的干擾, 低信噪比是MRS信號提取面臨的最大挑戰(zhàn)。

環(huán)境噪聲主要分為尖峰噪聲、 工頻諧波噪聲和隨機噪聲。尖峰噪聲是由天空中打雷放電、 開關(guān)突然的開閉等引起的瞬時大幅度干擾, 會對MRS信號的后續(xù)濾波疊加處理造成非常嚴(yán)重的影響, 一般是首先去除的一類噪聲。目前, 比較成熟的去噪方法有能量運算方法[2]、 尖峰建模方法[3]、 統(tǒng)計疊加方法和小波模極大值方法[4]。工頻諧波噪聲是由高壓輸電線、 變壓器和其他電器設(shè)備產(chǎn)生的, 在測量數(shù)據(jù)中最常出現(xiàn), 影響最大。由于工頻諧波噪聲在多通道采集時具有較強的相關(guān)性, 因此可以采用自適應(yīng)噪聲抵消[5]等方法進(jìn)行消除。隨機噪聲是一種沒有規(guī)律, 但統(tǒng)計規(guī)律服從高斯分布的噪聲, 分為加性隨機噪聲與乘性隨機噪聲兩類。對加性噪聲有基于L1范數(shù)的低場核磁共振T2譜稀疏反演方法[6]、 稀疏表示法[7]、 時頻峰值濾波法[8]等方法降低對MRS信號的干擾, 但對乘性隨機噪聲還沒有成熟的處理算法。

筆者針對儀器采集數(shù)據(jù)中的加性和乘性隨機噪聲, 使用馬爾科夫鏈蒙特卡洛(MCMC: Markov Chain Monte Carlo)算法[9-14]提取MRS復(fù)包絡(luò)信號的參數(shù)。常規(guī)的非線性擬合方法不適用于乘性隨機噪聲, 而傳統(tǒng)的蒙特卡洛算法的計算量大, 筆者將MCMC方法運用于MRS信號提取, 不僅可以獲得最優(yōu)的算法性能, 而且其計算復(fù)雜度也不會隨著樣本次數(shù)的增加而明顯增加, 僅與迭代次數(shù)成線性關(guān)系。

1 磁共振數(shù)據(jù)信噪特征

根據(jù)文獻(xiàn)[1], 通過對采集到的MRS數(shù)據(jù)進(jìn)行去噪處理, 消除數(shù)據(jù)中的尖峰噪聲和工頻諧波噪聲, 再經(jīng)過希爾伯特變換和低通濾波處理后, 可獲得MRS信號的復(fù)包絡(luò)表達(dá)式

(1)

(2)

在采集MRS信號時, 環(huán)境噪聲的影響很大, 其中隨機噪聲會淹沒MRS信號。一般隨機噪聲服從高斯分布, 其噪聲水平根據(jù)環(huán)境不同達(dá)到200 nV以上。對于在MRS數(shù)據(jù)Vi中的隨機噪聲, 可以分為加性噪聲Vp∈N(μp,σp)和乘性噪聲Vm∈N(μm,σm)兩類, 分布模型為

(3)

其中乘性噪聲的噪聲水平與MRS信號幅度VMRS相關(guān),σmi=|eln(α)Vi|, ln(α)是乘性噪聲的系數(shù)。因此在受到隨機噪聲影響后的MRS復(fù)包絡(luò)數(shù)據(jù)為

Vnoi=VMRS+Vp+Vm

(4)

2 馬爾可夫鏈蒙特卡洛參數(shù)提取

2.1 馬爾可夫鏈的蒙特卡洛

MH(Metropolis Hasting)算法是常用的MCMC采樣方法之一。MH算法的基本思路是在每步狀態(tài)由x移動到新的狀態(tài)x′時, 其概率為q(x′|x),q稱為MH算法的建議概率密度分布(也稱為轉(zhuǎn)換核), 同時滿足以下條件

π(x)q(x′|x)=π(x′)q(x|x′)

(5)

其中π(x)是目標(biāo)分布(即待采樣的分布)。滿足式(5)的條件也稱為由一條馬爾科夫鏈固定的概率密度函數(shù), 但容易造成該馬爾科夫鏈性質(zhì)不相等, 因此需要添加額外的參數(shù)條件, 即

π(x)q(x′|x)α(x′|x)=π(x′)q(x|x′)α(x|x′)

(6)

其中α(x|x′)定義為狀態(tài)從x′移動到x的概率;q(x|x′)α(x|x′)稱為更新后的建議概率密度分布。為了保證這個式子成立, 可以設(shè)α(x|x′)≡1, 得到一次移動的概率為

(7)

其中α(x′|x)是狀態(tài)從x移動到x′的概率, 也稱為接受概率。如果建議被接受, 則新狀態(tài)被賦值為x′, 否則新狀態(tài)仍為x。

MH算法的實現(xiàn)步驟為：

1) 設(shè)s=0, 生成一個初始值x, 即初始狀態(tài)x(0)=x;

2) 令s∶=s+1;

3) 從建議分布中生成一個建議值x′;

4) 根據(jù)式(7)計算接受概率α(x′|x);

5) 從均勻分布U～(0,1)生成一個隨機數(shù)u;

6) 如果u<α(x′|x), 則令x(s+1)=x′; 反之, 則令x(s+1)=x(s)=x;

7) 循環(huán)執(zhí)行步驟2)～步驟6), 直到s=N(N代表迭代的次數(shù), 保證在馬爾科夫鏈下的目標(biāo)分布達(dá)到平穩(wěn)分布)。

2.2 參數(shù)提取方法

在貝葉斯條件概率的框架下, 通過MRS復(fù)包絡(luò)數(shù)據(jù)推算參數(shù)的問題可以表示為： 在給定觀測數(shù)據(jù)V和參數(shù)先驗信息M條件下, 估計關(guān)于參數(shù)的后驗分布概率密度函數(shù)

(8)

其中p(V|M)代表似然函數(shù), 表示MRS信號復(fù)包絡(luò)數(shù)據(jù)與未知參數(shù)之間的相互關(guān)系;p(M)是參數(shù)的先驗概率分布;p(V)可認(rèn)為是一個僅取決于觀測數(shù)據(jù)的常數(shù), 僅對后驗分布起到歸一化的作用, 因此可以忽略不計。進(jìn)而在p(V)為定值的條件下有

p(M|V)∝p(V|M)p(M)

(9)

其中后驗分布概率函數(shù)p(M|V)就是各參數(shù)的解。

根據(jù)文獻(xiàn)[4], 由于噪聲服從高斯分布, 而且各參數(shù)之間相互獨立, 因此似然函數(shù)p(V|M)可定義為一種多維的正態(tài)分布

(10)

其中V是MRS復(fù)包絡(luò)觀測數(shù)據(jù);f(M)是給定參數(shù)的MRS復(fù)包絡(luò)表達(dá)式的預(yù)測數(shù)據(jù);N是觀測數(shù)據(jù)的個數(shù);Cd是觀測數(shù)據(jù)的協(xié)方差矩陣, 由每個時間點對應(yīng)的預(yù)測值與觀測數(shù)據(jù)誤差方差組成的一個對角矩陣; |Cd|表示協(xié)方差矩陣Cd的行列式。

通過貝葉斯框架的建立, 可求解出后驗概率分布p(M|V), 但由于p(M|V)表達(dá)式較為復(fù)雜, 抽樣維度高, 需要抽樣的樣本量大, 采用蒙特卡洛的方法完全列舉后驗分布或?qū)篁灧植歼M(jìn)行均勻采樣都是不符合實際的, 因此這里使用MCMC方法實現(xiàn)后驗分布采樣。

在計算后驗分布時, 其后驗概率分布p(M|V)就是上述MH算法中的目標(biāo)分布π(x), 先驗概率分布p(M)即為上述建議概率密度分布。由于建立的先驗概率分布模型p(M)服從高斯分布, 具有對稱性, 滿足p(M′|M)=p(M|M′)。因此化簡得到接受概率為

(11)

通過MCMC模型實現(xiàn)MRS復(fù)包絡(luò)參數(shù)提取的步驟如下:

1) 對采集得到的MRS數(shù)據(jù)通過已有方法去除尖峰噪聲和工頻噪聲;

2) 通過希爾伯特變換以及低通濾波處理得到MRS復(fù)包絡(luò)數(shù)據(jù)V;

3) 設(shè)非線性擬合提取MRS復(fù)包絡(luò)各參數(shù), 確定參數(shù)的先驗信息概率密度分布p(M);

5) 設(shè)置迭代次數(shù)為N,迭代步長為Δs, 令s∶=s+Δs;

6) 從先驗概率分布模型p(M)中生成一組建議參數(shù)向量H;

7) 根據(jù)式(11)計算接受概率α(H|M);

8) 從均勻分布U～(0,1)生成一個隨機數(shù)u;

9) 如果u<α(H|M), 則令M(s+Δ s)=H; 反之, 則令M(s+Δ s)=M(s)=M;

10) 循環(huán)執(zhí)行步驟2)～步驟6), 直到s=N, 得到采樣后的后驗分布函數(shù)p(M|V)。

使用MCMC方法采樣后驗分布的收斂條件,是指采樣得到的樣本達(dá)到穩(wěn)定狀態(tài)。在馬爾可夫鏈達(dá)到穩(wěn)定后, 其之前狀態(tài)得到的數(shù)據(jù)都會被拋棄。當(dāng)?shù)螖?shù)N一定, 如果迭代步長Δs過小, 則會出現(xiàn)無法收斂的情況。在文獻(xiàn)[2]中提到了一種初始設(shè)置方案, 可選擇取鏈條數(shù)為3, 采樣迭代次數(shù)為1 000 000次, 默認(rèn)丟棄前50%樣本后收斂, 達(dá)到穩(wěn)定。

3 數(shù)值模擬

3.1 仿真實驗

a 包絡(luò)實部提取結(jié)果對比 b 包絡(luò)虛部提取結(jié)果對比圖1 非線性擬合、 MCMC擬合結(jié)果與無噪聲MRS復(fù)包絡(luò)信號對比Fig.1 Comparison of nonlinear fitting, MCMC fitting and noiseless complex envelope of MRS

通過圖1可以發(fā)現(xiàn), 當(dāng)引入乘性噪聲后, 非線性擬合效果下降。而通過MCMC方法擬合參數(shù)受到的影響較小, 依然維持在較好的擬合效果(見表1)。

表1 參數(shù)提取結(jié)果對比

3.2 穩(wěn)定性分析

選取加性噪聲水平分別為100 nV、 150 nV、 200 nV和500 nV, 乘性噪聲系數(shù)為lg(1.5)。在進(jìn)行多組測量記錄結(jié)果后, 繪制參數(shù)的箱圖如圖2所示。由圖2可見, 在乘性噪聲系數(shù)一定的情況下, MCMC參數(shù)提取方法可以保證結(jié)果準(zhǔn)確度。隨著加性噪聲水平增大, 兩種方法的參數(shù)準(zhǔn)確度逐漸降低, 但MCMC方法準(zhǔn)確度仍然可以保持穩(wěn)定。

圖2 100組仿真數(shù)據(jù)的MRS包絡(luò)信號參數(shù)提取結(jié)果統(tǒng)計Fig.2 100 sets simulation data of parameters extraction from MRS envelope signals

選取乘性噪聲水平系數(shù)為lg(1.5)、lg(2.0)、lg(3.0)、lg(5.0), 加性噪聲水平為100 nV,在進(jìn)行多組測量記錄結(jié)果后, 繪制參數(shù)的箱圖如圖3所示。由圖3可見, 在加性噪聲水平一定的情況下, MCMC參數(shù)提取方法可保證結(jié)果準(zhǔn)確度。隨著乘性噪聲系數(shù)增大,兩種方法的參數(shù)準(zhǔn)確度逐漸降低, 但MCMC方法準(zhǔn)確度仍然可以保持穩(wěn)定。

圖3 100組仿真數(shù)據(jù)的MRS包絡(luò)信號參數(shù)提取結(jié)果統(tǒng)計Fig.3 100 sets of simulation data of parameters extraction from MRS envelope signals

4 結(jié) 語

筆者針對已有利用非線性擬合提取MRS信號復(fù)包絡(luò)表達(dá)式參數(shù)方法的基礎(chǔ)上, 利用非線性擬合的方法, 實現(xiàn)了建立各參數(shù)的先驗信息分布和參數(shù)與數(shù)據(jù)的似然函數(shù), 通過MCMC方法中的MH算法采樣并求解參數(shù)的后驗分布, 找到后驗分布中出現(xiàn)次數(shù)最多、 權(quán)值最大的數(shù)值作為參數(shù)最優(yōu)估計值。在不同加性噪聲水平和不同乘性噪聲系數(shù)的條件下仿真參數(shù)提取實驗, 將MCMC參數(shù)提取方法與非線性擬合結(jié)果對比, 數(shù)值模擬結(jié)果證明, 在乘性與加性隨機噪聲的影響下, MCMC參數(shù)提取方法優(yōu)于非線性擬合參數(shù)提取結(jié)果。

綜上, 基于馬爾科夫鏈蒙特卡洛的地面磁共振信號參數(shù)提取方法可以有效減少了由于非線性擬合無法處理乘行噪聲所帶來的誤差。雖然需要的運算時間變長, 但保證了參數(shù)的估計精度與穩(wěn)定性, 即解決了乘性噪聲和加性噪聲共同影響參數(shù)提取的問題。