摘 要： 針對(duì)核脈沖全譜數(shù)據(jù)采集技術(shù)中采用的深度加權(quán)濾波法在數(shù)據(jù)量少時(shí)平滑效果差、滯后大的缺點(diǎn)，設(shè)計(jì)了卡爾曼濾波算法替代傳統(tǒng)濾波算法進(jìn)行解譜的方案。濾波算法首先根據(jù)信號(hào)與噪聲的狀態(tài)空間模型，建立狀態(tài)方程和測(cè)量方程，然后根據(jù)廣義卡爾曼濾波對(duì)測(cè)量方程進(jìn)行更新，最后根據(jù)現(xiàn)場(chǎng)的標(biāo)準(zhǔn)刻度井測(cè)量數(shù)據(jù)確定測(cè)量矩陣，對(duì)測(cè)量數(shù)據(jù)進(jìn)行濾波。實(shí)測(cè)結(jié)果表明，該濾波算法相對(duì)于傳統(tǒng)濾波算法，能夠消除譜圖的統(tǒng)計(jì)漲落，提高光滑度，遏制滯后，并保證濾波后曲線不失真。

　　關(guān)鍵詞：卡爾曼濾波；數(shù)據(jù)處理；核脈沖全譜數(shù)據(jù)采集系統(tǒng)；解譜

0 引 言

　　核脈沖全譜數(shù)據(jù)采集技術(shù)是新一代放射性測(cè)井儀器的核心技術(shù)。由系統(tǒng)的放射源發(fā)出662 kev γ射線，在地層中經(jīng)過一次或多次的康普頓散射和光電效應(yīng)，經(jīng)物理探測(cè)窗口進(jìn)入NaI晶體，產(chǎn)生與入射能量成正比例的光電信號(hào)。 光電信號(hào)經(jīng)光電倍增管形成核電子脈沖信號(hào)。核脈沖信號(hào)A/D轉(zhuǎn)換，分道測(cè)量計(jì)數(shù)[1]，上傳地面計(jì)算機(jī)，形成全譜數(shù)據(jù)譜線。由于計(jì)數(shù)器存在嚴(yán)重的統(tǒng)計(jì)漲落，必須去偽存真，以便資料分析。

　　其信息處理技術(shù)的關(guān)鍵是實(shí)測(cè)譜的解析技術(shù)。解譜方法采用深度加權(quán)濾波法。公式如下：

　　 （1）

　　其中，y是n個(gè)采樣值的加權(quán)平均值，yi是加權(quán)因子，ai是加權(quán)系數(shù)，各個(gè)加權(quán)系數(shù)均為小于1的小數(shù)，且滿足總和等于1的約束條件。計(jì)算雖然簡(jiǎn)單，但是需確定歸一化因子的值，即需要確定ai的值，且在數(shù)據(jù)量少時(shí)平滑效果差、滯后大。

　　本文設(shè)計(jì)的卡爾曼濾波能夠消除譜圖的統(tǒng)計(jì)漲落，提高光滑度，遏制滯后，實(shí)現(xiàn)測(cè)井?dāng)?shù)據(jù)的有效濾波。

1 卡爾曼濾波的模型建立

　　1.1 建立系統(tǒng)的狀態(tài)方程

　　卡爾曼濾波以最小均方誤差為估計(jì)的最佳準(zhǔn)則，以此來尋求一套遞推估計(jì)的算法，其基本思想是：采用信號(hào)與噪聲的狀態(tài)空間模型，利用前一時(shí)刻的估計(jì)值和現(xiàn)時(shí)刻的觀測(cè)值來更新對(duì)狀態(tài)變量的估計(jì)，求出現(xiàn)時(shí)刻的估計(jì)值。它適合實(shí)時(shí)處理和計(jì)算機(jī)運(yùn)算[2]。

　　卡爾曼濾波注重對(duì)象的物理過程描述，要求先建立系統(tǒng)的狀態(tài)方程，由于地層的變化是未知的，很難用一組確定的微分方程加以描述，因此采用一階高斯馬爾科夫過程描述地層變化。

　　首先設(shè)核脈沖個(gè)數(shù)α和每道脈沖能量β為狀態(tài)變量(X1，X2)[3]，有：

　　（2）

　　（3）

　　 （4）

　　 （5）

　　式（2）~式（5）中，h表示井的深度，λ1(h)、λ2(h)是非奇異狀態(tài)轉(zhuǎn)移矩陣，是隨井深變化的斜率參數(shù)。W1(h)~W4(h)是動(dòng)態(tài)噪聲，方差為，均值為零。

　　參數(shù)λ1(h),λ2(h)是局部隨機(jī)變化的。若測(cè)井井段類似于均勻的地層區(qū)， f取得足夠小時(shí)，λ1(h)、λ2(h)趨于零，則有當(dāng)測(cè)井井段跨越一個(gè)地層界面時(shí)，λ1(h)、λ2(h)根據(jù)X1(h)、X2(h)的變化斜率取正值或負(fù)值。若把λ1(h)、λ2(h)設(shè)為增廣狀態(tài)變量X3(h)、X4(h)，則狀態(tài)方程可表示為：

　　（6）

　　f表示式（2）~式（5）的非線性關(guān)系。式（6）描述了地層隨井深的變化，同時(shí)還需要建立反映測(cè)量信息與狀態(tài)變量間的對(duì)應(yīng)關(guān)系的測(cè)量方程。

　　1.2 最優(yōu)測(cè)量方程的建立

　　設(shè)長(zhǎng)、短道源采集的核脈沖的個(gè)數(shù)X1(h)和每個(gè)脈沖的能量X2(h)的測(cè)量計(jì)數(shù)率為Z(h)，設(shè)，則其關(guān)系可表示為:

　　（7）

　　由于增廣的狀態(tài)為4個(gè)，最后測(cè)量矩陣為4×4階矩陣。

　　測(cè)量方程為：

　　 （8）

　　式（7）反映了狀態(tài)變量式（6）與測(cè)量方程式（8）之間的對(duì)應(yīng)關(guān)系。由于狀態(tài)方程式（6）是非線性的，而測(cè)量方程式（8）是線性的，所以需要采用廣義卡爾曼濾波對(duì)測(cè)量方程（8）進(jìn)行濾波更新。

　　下面用狀態(tài)方程和系統(tǒng)測(cè)量方程結(jié)合它們的協(xié)方差來估算系統(tǒng)的最優(yōu)化輸出。利用過程模型來預(yù)測(cè)下一狀態(tài)的系統(tǒng)。假設(shè)現(xiàn)在的系統(tǒng)狀態(tài)是k，根據(jù)系統(tǒng)的模型，可以基于系統(tǒng)的上一狀態(tài)而預(yù)測(cè)出當(dāng)前狀態(tài)[4-5]：

　　式(9)中，狀態(tài)微分方程組右端項(xiàng)；T表示系統(tǒng)的采樣間隔；其中：

　　 （10）

　　至此，系統(tǒng)結(jié)果已經(jīng)更新了，對(duì)應(yīng)于的協(xié)方差的更新如下（用P表示協(xié)方差）：

　　（11）

　　式(11)中Qn表示 W(k)動(dòng)態(tài)噪聲的協(xié)方差矩陣；λ(h)T表示λ(h)的轉(zhuǎn)置矩陣，且：

　　（12）

　　其中，I為單位矩陣。

　　式（9）、式（11）即是對(duì)系統(tǒng)的預(yù)測(cè)。

　　結(jié)合當(dāng)前狀態(tài)的預(yù)測(cè)值和當(dāng)前狀態(tài)的測(cè)量值，得到現(xiàn)在狀態(tài)k的最優(yōu)化估算值X(k,k)：

　　（13）

　　式（13）中，kg表示卡爾曼增益(Kalman Gain)：

　　 （14）

　　式（13）、（14）中，H表示測(cè)量矩陣；Rm表示測(cè)量噪聲的協(xié)方差矩陣。

　　為使卡爾曼濾波器不斷地運(yùn)行下去直到系統(tǒng)過程結(jié)束，再對(duì)k狀態(tài)下X(k,k)的協(xié)方差進(jìn)行如下更新：

　　（15）

　　式（9）~式（15）是對(duì)測(cè)量方程（8）的濾波更新，即得到最優(yōu)測(cè)量方程。結(jié)合式（6）~式（8）以及設(shè)定的參數(shù)即可對(duì)譜圖進(jìn)行濾波。

2 測(cè)量方程的參數(shù)設(shè)計(jì)

　　以上遞推公式中，需要設(shè)計(jì)式（7）以及式（9）~（15）中的參數(shù)。式（7）中，V1~V4為測(cè)量噪聲，H矩陣為4×2階測(cè)量矩陣，可由現(xiàn)場(chǎng)的標(biāo)準(zhǔn)刻度井測(cè)量數(shù)據(jù)來確定。需要說明的是，每個(gè)測(cè)井儀器都有一套自己的H矩陣，它可認(rèn)為是不隨其他因素變化而變化的常數(shù)矩陣。本文中，H1取為：

　　根據(jù)所建立的系統(tǒng)模型式（2）~式（6）以及式（8），對(duì)方程進(jìn)行離散線性化處理，則有：

　　（16）

　　（17）

　　（18）

　　式（16）~式（18）是式（9）~式（15）的設(shè)計(jì)參數(shù)，其中missing image file分別是W1(h)~W4(h)的方差，均值為0。missing image file分別是V1(k)~V4(k)的方差，均值為0。

3 實(shí)測(cè)濾波情況

　　結(jié)合以上遞推公式，在給定濾波初值（即X(0)）后，估計(jì)誤差協(xié)方差矩陣初值P(0,0)、采樣間隔T以及協(xié)方差陣(Qn,Rm)，即可以濾波。

　　在核脈沖全譜數(shù)據(jù)采集系統(tǒng)中，其控制量為儀器中直接貼于NaI晶體端面的一種低強(qiáng)度137Cs源，137Cs源發(fā)射的662 kev射線直接進(jìn)入晶體，γ衰變是核現(xiàn)象，與溫度無關(guān)，故與之相對(duì)應(yīng)的電子脈沖信號(hào)應(yīng)維持恒定幅度，并以其為參考值調(diào)整高壓，在譜圖處理中，它可以認(rèn)為是常數(shù)[6]，譜圖對(duì)比不受影響。

　　由于在測(cè)井過程中，儀器運(yùn)動(dòng)速度快，所以核脈沖采集周期一般為20 ms或40 ms。在本文設(shè)計(jì)中，核脈沖的采集周期為40 ms ，即ΔT=40 ms。測(cè)試取一幀的數(shù)據(jù)，即256道譜圖。

　　本文對(duì)實(shí)驗(yàn)數(shù)據(jù)進(jìn)行了處理，在相同的實(shí)驗(yàn)數(shù)據(jù)下，分別對(duì)硫、鋁進(jìn)行測(cè)試對(duì)比，通過長(zhǎng)、短源譜圖濾波，對(duì)卡爾曼濾波譜圖與深度加權(quán)濾波譜圖進(jìn)行比較，具體譜圖波形見圖1~圖4。

　　圖1~圖4中，齊的黑線如15.90、27.30等是參考道址，為了方便觀察譜圖是否滯后，便于對(duì)比。道址192.30~265之間是穩(wěn)譜峰圖，主要是為了增加分布在0~192.30道址之間核脈沖的能量強(qiáng)度，便于突出異同，在計(jì)算解釋譜圖時(shí)會(huì)去掉，計(jì)算不受影響。

　　圖1~圖4中，上圖是深度加權(quán)濾波得到的譜圖，下圖是卡爾曼濾波得到的譜圖，通過上下圖的對(duì)比可以看出，下圖的包絡(luò)明顯趨向平滑，能夠有效消除統(tǒng)計(jì)漲落的影響。通過兩組對(duì)比可以看出，在相同的實(shí)驗(yàn)數(shù)據(jù)下，此卡爾曼濾波得到的譜圖比深度加權(quán)濾波得到的譜圖包絡(luò)線更接近于平滑，有力消除了譜圖的統(tǒng)計(jì)漲落，更利于測(cè)井?dāng)?shù)據(jù)的分析。

4 結(jié)論

　　文中通過對(duì)核脈沖全譜數(shù)據(jù)采集系統(tǒng)的解譜需求，設(shè)計(jì)了一種卡爾曼函數(shù)濾波器，通過實(shí)驗(yàn)譜圖對(duì)比可以看出，其包絡(luò)線幾乎接近于平滑，消除了譜圖的統(tǒng)計(jì)漲落，提高了光滑度。說明了本文所設(shè)計(jì)的卡爾曼濾波對(duì)譜圖處理的有效性和可行性。

參考文獻(xiàn)

　　[1] 彭曉光，柏林，田彥民，等. 基于FPGA的全譜巖性密度數(shù)據(jù)采集系統(tǒng)[J]. 核電子學(xué)與探測(cè)技術(shù), 2012,32(7):758-760.

　　[2] 邵玉華. 卡爾曼波形估計(jì)在雷達(dá)信號(hào)處理中的應(yīng)用[J]. 黑龍江科技信息, 2012(2):77.

　　[3] 彭丁聰. 卡爾曼濾波的基本原理及應(yīng)用[J]. 軟件導(dǎo)刊, 2009(11):33-34.

　　[4] 張永偉，楊鎖昌，張敏，等. 卡爾曼濾波在落點(diǎn)偏差預(yù)測(cè)算法中的應(yīng)用[J]. 中國(guó)測(cè)試, 2012,38(s1):112-113.

　　[5] 黃亞萍. 基于卡爾曼濾波的抑制NLOS定位算法研究[J]. 電信快報(bào), 2012(4):48-51.

　　[6] 彭曉光. 新型全譜巖性密度測(cè)井儀的研制[D]. 邯鄲:中船重工718研究所, 2012,44-47.