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

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

0 引 言

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

　　其信息處理技術的關鍵是實測譜的解析技術。解譜方法采用深度加權濾波法。公式如下：

　　 （1）

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

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

1 卡爾曼濾波的模型建立

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

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

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

　　首先設核脈沖個數(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)趨于零，則有當測井井段跨越一個地層界面時，λ1(h)、λ2(h)根據(jù)X1(h)、X2(h)的變化斜率取正值或負值。若把λ1(h)、λ2(h)設為增廣狀態(tài)變量X3(h)、X4(h)，則狀態(tài)方程可表示為：

　　（6）

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

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

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

　　（7）

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

　　測量方程為：

　　 （8）

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

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

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

　　 （10）

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

　　（11）

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

　　（12）

　　其中，I為單位矩陣。

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

　　結合當前狀態(tài)的預測值和當前狀態(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)過程結束，再對k狀態(tài)下X(k,k)的協(xié)方差進行如下更新：

　　（15）

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

2 測量方程的參數(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ù)，其中missing image file分別是W1(h)~W4(h)的方差，均值為0。missing image file分別是V1(k)~V4(k)的方差，均值為0。

3 實測濾波情況

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

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

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

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

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

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

4 結論

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

參考文獻

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

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

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

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

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

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