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

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

0 引 言

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

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

　　 （1）

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

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

1 卡爾曼濾波的模型建立

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

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

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

　　首先設(shè)核脈沖個數(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)是動態(tài)噪聲，方差為，均值為零。

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

　　（6）

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

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

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

　?。?）

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

　　測量方程為：

　　 （8）

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

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

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

　　 （10）

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

　　（11）

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

　　（12）

　　其中，I為單位矩陣。

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

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

　　（13）

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

　　 （14）

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

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

　　（15）

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

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

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

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

　　（16）

　?。?7）

　　（18）

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

3 實測濾波情況

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

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

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

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

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

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

4 結(jié)論

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

參考文獻

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

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

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

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

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

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