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

　　關鍵詞：卡爾曼濾波；數據處理；核脈沖全譜數據采集系統；解譜

0 引 言

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

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

　　 （1）

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

　　本文設計的卡爾曼濾波能夠消除譜圖的統計漲落，提高光滑度，遏制滯后，實現測井數據的有效濾波。

1 卡爾曼濾波的模型建立

　　1.1 建立系統的狀態方程

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

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

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

　　（2）

　　（3）

　　 （4）

　　 （5）

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

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

　　（6）

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

　　1.2 最優測量方程的建立

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

　　（7）

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

　　測量方程為：

　　 （8）

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

　　下面用狀態方程和系統測量方程結合它們的協方差來估算系統的最優化輸出。利用過程模型來預測下一狀態的系統。假設現在的系統狀態是k，根據系統的模型，可以基于系統的上一狀態而預測出當前狀態[4-5]：

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

　　 （10）

　　至此，系統結果已經更新了，對應于的協方差的更新如下（用P表示協方差）：

　　（11）

　　式(11)中Qn表示 W(k)動態噪聲的協方差矩陣；λ(h)T表示λ(h)的轉置矩陣，且：

　　（12）

　　其中，I為單位矩陣。

　　式（9）、式（11）即是對系統的預測。

　　結合當前狀態的預測值和當前狀態的測量值，得到現在狀態k的最優化估算值X(k,k)：

　　（13）

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

　　 （14）

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

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

　　（15）

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

2 測量方程的參數設計

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

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

　　（16）

　　（17）

　　（18）

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

3 實測濾波情況

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

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

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

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

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

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

4 結論

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

參考文獻

　　[1] 彭曉光，柏林，田彥民，等. 基于FPGA的全譜巖性密度數據采集系統[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.