馬爾科夫鏈蒙特卡洛地面磁共振信號(hào)參數(shù)提取

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

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

0 引 言

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

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

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

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

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

(1)

(2)

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

(3)

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

Vnoi=VMRS+Vp+Vm

(4)

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

2.1 馬爾可夫鏈的蒙特卡洛

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

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

(5)

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

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

(6)

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

(7)

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

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

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

2) 令s∶=s+1;

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

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

5) 從均勻分布U～(0,1)生成一個(gè)隨機(jī)數(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ù)提取方法

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

(8)

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

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

(9)

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

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

(10)

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

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

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

(11)

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

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

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

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

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

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

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

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

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

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

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

3 數(shù)值模擬

3.1 仿真實(shí)驗(yàn)

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

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

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

3.2 穩(wěn)定性分析

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

圖2 100組仿真數(shù)據(jù)的MRS包絡(luò)信號(hào)參數(shù)提取結(jié)果統(tǒng)計(jì)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)行多組測(cè)量記錄結(jié)果后, 繪制參數(shù)的箱圖如圖3所示。由圖3可見(jiàn), 在加性噪聲水平一定的情況下, MCMC參數(shù)提取方法可保證結(jié)果準(zhǔn)確度。隨著乘性噪聲系數(shù)增大,兩種方法的參數(shù)準(zhǔn)確度逐漸降低, 但MCMC方法準(zhǔn)確度仍然可以保持穩(wěn)定。

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

4 結(jié) 語(yǔ)

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

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