基于Unscented卡爾曼濾波器的反應(yīng)堆周期計(jì)算算法研究

熊華勝，李 鐸

（清華大學(xué) 核能與新能源技術(shù)研究院 先進(jìn)反應(yīng)堆工程與安全教育部重點(diǎn)實(shí)驗(yàn)室，北京 100084）

核反應(yīng)堆周期T 直接關(guān)系到核反應(yīng)堆的運(yùn)行和安全，需專門儀器來監(jiān)測(cè)［1］。周期太短表示反應(yīng)堆內(nèi)中子通量或功率水平增長(zhǎng)太快，如果周期超過整定值就需觸發(fā)保護(hù)系統(tǒng)緊急停堆動(dòng)作，關(guān)閉反應(yīng)堆。若不能及時(shí)準(zhǔn)確地計(jì)算出反應(yīng)堆周期，有潛在地引起核事故的可能性。

本文基于Unscented卡爾曼濾波器（UKF）對(duì)反應(yīng)堆周期的計(jì)算算法進(jìn)行研究。

1 反應(yīng)堆周期計(jì)算的離散隨機(jī)狀態(tài)空間模型

1.1 反應(yīng)堆周期的定義

反應(yīng)堆周期是指核反應(yīng)堆內(nèi)中子密度增加到初始時(shí)的e倍所用的時(shí)間。

式中：n0為t=0 時(shí)刻反應(yīng)堆內(nèi)的中子密度；n為t時(shí)刻的中子密度。

反應(yīng)堆周期反映了核反應(yīng)堆所處狀態(tài)。

1.2 離散隨機(jī)狀態(tài)空間模型推導(dǎo)

在數(shù)字化核測(cè)量系統(tǒng)中，通常在系統(tǒng)的前置放大模塊中將探測(cè)器輸出的弱電流信號(hào)I進(jìn)行放大，然后轉(zhuǎn)變?yōu)殡妷盒盘?hào)，再使用A／D 模塊實(shí)現(xiàn)對(duì)電壓信號(hào)的采集。根據(jù)核反應(yīng)堆工程原理可知：反應(yīng)堆功率與中子密度呈正比；探測(cè)器輸出的電流信號(hào)與反應(yīng)堆功率呈正比［2］。因此，電壓信號(hào)V 直接反應(yīng)了核反應(yīng)堆的功率水平，即有下面的關(guān)系式：

其中，V0為0時(shí)刻的電壓信號(hào)。

將式（3）離散化得到：

其中：Tc為V（k+1）與V（k）之間的時(shí)間間隔，s；T（k）為k時(shí)刻的反應(yīng)堆周期。

在實(shí)際的工程中，Tc一般較小，為ms級(jí)，而T（k）一般至少為s級(jí)，通常滿足Tc／T（k）?1。對(duì)式（3）進(jìn)行泰勒展開，得到：

其中：θ（k）=Tc／T（k）。

令w1（k）=V（k）·ο（θ（k）），可得：

并假設(shè)反應(yīng)堆周期T（k+1）與T（k）滿足T（k+1）=T（k）+ο（T（k）），o（T（k））為遠(yuǎn)小于T（k）的增量，令：

可得：

在數(shù)字化核測(cè)量系統(tǒng)中，將探測(cè)器輸出的電流信號(hào)轉(zhuǎn)變?yōu)殡妷盒盘?hào)，可表示為：

其中：y（k）為觀測(cè)向量；v（k）為測(cè)量噪聲。

由式（6）、（7）、（8）得到離散的非線性模型：

其中：

w（k）是非線性模型的系統(tǒng)噪聲，是對(duì)系統(tǒng)建模產(chǎn)生的誤差進(jìn)行補(bǔ)償。式（9）即為用于計(jì)算核反應(yīng)堆周期的離散隨機(jī)狀態(tài)空間模型。

2 基于UKF的反應(yīng)堆周期的計(jì)算方法

2.1 反應(yīng)堆周期的計(jì)算公式

對(duì)于式（9），如果根據(jù)測(cè)量電壓y（k）估計(jì)出離散非線性模型的狀態(tài)向量x（k），那么就可計(jì)算出核反應(yīng)堆周期T（k），即：

因此，反應(yīng)堆周期的計(jì)算問題就轉(zhuǎn)化為離散隨機(jī)狀態(tài)空間模型的狀態(tài)估計(jì)問題。

2.2 UKF算法

由于擴(kuò)展卡爾曼濾波器（EKF）方法是將非線性函數(shù)利用泰勒公式展開，并保留一階項(xiàng)［3］，從而實(shí)現(xiàn)非線性函數(shù)的線性化并保留一階精度。針對(duì)EKF 的缺陷，Julier等［4］提出了基于Unscented變換的UKF 方法。該方法直接利用非線性模型進(jìn)行遞推估計(jì)，避免了線性化誤差的引入，且遞推估計(jì)過程不需計(jì)算Jacobian矩陣，比EKF 更為簡(jiǎn)單，并在執(zhí)行濾波時(shí)可有二階以上的精度。

假設(shè)w（k）、v（k）的統(tǒng)計(jì)特征為：

選取系統(tǒng)增廣狀態(tài)變量為：

記xak為L(zhǎng) 維列向量，則L=2+2+1=5。增廣后的系統(tǒng)狀態(tài)方差為：

UKF狀態(tài)估計(jì)算法如下。

1）初始化

2）計(jì)算Sigma點(diǎn)

其中，χk為L(zhǎng)×（2L+1）維矩陣。

3）時(shí)間更新方程

4）測(cè)量更新方程

算法中對(duì)應(yīng)的權(quán)系數(shù)為：

其中，λ=α2（L+κ）-L。α和κ 用于控制Sigma點(diǎn)的傳播，β 是一與x 的分布相關(guān)的系數(shù)。這里取α=0.001，κ=0和β=2。

利用式（11）～（13）給出的UKF算法可估計(jì)出式（9）離散隨機(jī)狀態(tài)空間模型的狀態(tài)，再依據(jù)式（10）可以計(jì)算出反應(yīng)堆周期。

3 反應(yīng)堆周期仿真計(jì)算與實(shí)驗(yàn)驗(yàn)證

利用建立起來的離散隨機(jī)狀態(tài)空間模型和設(shè)計(jì)好的UKF 算法，分別利用仿真數(shù)據(jù)和實(shí)驗(yàn)數(shù)據(jù)來驗(yàn)證本文提出的反應(yīng)堆周期的計(jì)算算法的正確性和有效性。

3.1 仿真計(jì)算

在仿真過程中利用式（4）產(chǎn)生反應(yīng)堆周期為28.85s 的仿真測(cè)量信號(hào)，并對(duì)測(cè)量信號(hào)V（k）施加具有均值為0、方差Rk為0.000 01統(tǒng)計(jì)特征的高斯噪聲信號(hào)。

利用UKF 方法，對(duì)生成的仿真數(shù)據(jù)進(jìn)行處理得到相應(yīng)的結(jié)果如圖1所示。從仿真計(jì)算結(jié)果來看，噪聲雖然幅值很小，但對(duì)計(jì)算的影響卻非常大，造成計(jì)算結(jié)果的不確定性大。應(yīng)用UKF可較好地復(fù)現(xiàn)信號(hào)的本來特征。

圖1 反應(yīng)堆周期的仿真計(jì)算結(jié)果Fig.1 Calculation result of reactor period using simulation data

3.2 實(shí)驗(yàn)驗(yàn)證

在數(shù)字化核測(cè)量系統(tǒng)的原理樣機(jī)中，利用信號(hào)發(fā)生器生成的信號(hào)來模擬實(shí)際的電離室探測(cè)器輸出信號(hào)，并通過這種半物理實(shí)驗(yàn)的方式來獲取實(shí)驗(yàn)數(shù)據(jù)（圖2）。

圖2 數(shù)字化核測(cè)量系統(tǒng)原理圖Fig.2 Schematic of digital neutron monitoring system

信號(hào)發(fā)生器用于模擬探測(cè)器的輸出信號(hào)，產(chǎn)生反應(yīng)堆周期為T 的弱電流信號(hào)。前置處理模塊用于信號(hào)的調(diào)理和隔離，對(duì)微電流信號(hào)進(jìn)行放大，并將電流信號(hào)調(diào)理為電壓信號(hào)。信號(hào)采集模塊在計(jì)算機(jī)系統(tǒng)中用于對(duì)電壓信號(hào)進(jìn)行A／D 轉(zhuǎn)換；信號(hào)處理及反應(yīng)堆周期計(jì)算模塊用于對(duì)采集的數(shù)字信號(hào)進(jìn)行變換、濾波處理，然后用本文提出的方法計(jì)算得出反應(yīng)堆周期。

用信號(hào)發(fā)生器生成反應(yīng)堆周期為28.85s的電流信號(hào)。實(shí)驗(yàn)過程中，信號(hào)發(fā)生器產(chǎn)生的電流信號(hào)經(jīng)前置處理模塊后電壓范圍為0.003～3V。數(shù)字化核測(cè)量系統(tǒng)的信號(hào)處理模塊通過采集、處理和計(jì)算后得到的計(jì)算結(jié)果示于圖3。從圖3可看出，輸入到核測(cè)量系統(tǒng)計(jì)算機(jī)的信號(hào)經(jīng)過采集模塊后，采集的電壓信號(hào)中總是伴隨著測(cè)量噪聲，測(cè)量噪聲對(duì)計(jì)算結(jié)果有較大影響。反應(yīng)堆周期的計(jì)算結(jié)果中包含著測(cè)量噪聲對(duì)周期計(jì)算的貢獻(xiàn)，導(dǎo)致計(jì)算后的反應(yīng)堆周期存在較大的不確定性和隨機(jī)性，使得周期的計(jì)算結(jié)果嚴(yán)重失真，與實(shí)際的期望值存在很大的差別。利用本文基于UKF提出的反應(yīng)堆周期計(jì)算方法可較好地抑制噪聲影響，得出較為理想的計(jì)算結(jié)果，從而復(fù)現(xiàn)了信號(hào)的本來特征。

圖3 反應(yīng)堆周期的半物理實(shí)驗(yàn)計(jì)算結(jié)果Fig.3 Calculation result of reactor period using experimental data

4 結(jié)論

數(shù)字化核測(cè)量系統(tǒng)將電離室探測(cè)器輸出的弱電流信號(hào)經(jīng)過前置放大、調(diào)理和隔離后將信號(hào)轉(zhuǎn)變?yōu)榕c反應(yīng)堆功率呈正比的電壓信號(hào)，核測(cè)量系統(tǒng)計(jì)算機(jī)采集的電壓信號(hào)總是伴隨著測(cè)量噪聲，從而對(duì)反應(yīng)堆周期的計(jì)算帶來很大影響。本文基于UKF 對(duì)反應(yīng)堆周期的計(jì)算進(jìn)行了研究，推導(dǎo)出適于使用UKF的離散隨機(jī)狀態(tài)空間模型，提出了反應(yīng)堆周期估計(jì)算法，并利用仿真數(shù)據(jù)和半物理實(shí)驗(yàn)數(shù)據(jù)對(duì)提出的算法進(jìn)行了計(jì)算驗(yàn)證。計(jì)算結(jié)果表明：本文提出的反應(yīng)堆周期計(jì)算算法可較好地抑制噪聲的影響，得到較為理想的計(jì)算結(jié)果，從而復(fù)現(xiàn)信號(hào)的本來特征。

［1］ 陳小軍，李保祥，李凱.數(shù)字化核測(cè)儀表的設(shè)計(jì)與算法分析［J］.核電子學(xué)與探測(cè)技術(shù)，2003，23（2）：163-165，175.CHEN Xiaojun，LI Baoxiang，LI Kai.The design＆arithmetic analysis of a digital nuclear instrument［J］.Nuclear Electronics ＆Detection Technology，2003，23（2）：163-165，175（in Chinese）.

［2］ 胡守印.反應(yīng)堆周期監(jiān)測(cè)裝置的研制［J］.工業(yè)儀表與自動(dòng)化裝置，2000（4）：47-48，6.HU Shouyin.The development of a period monitor for nuclear reactors［J］.Industrial Instrumentation＆Automation，2000（4）：47-48，6（in Chinese）.

［3］ NEMRA A，AOUF N.Robust INS／GPS sensor fusion for UVA localization using SDRE nonlinear filtering［J］.IEEE Sensors Journal，2010，11（4）：789-798.

［4］ JULIER S，UHLMANN J，WHYTE H F D.A new method for the nonlinear transformation of means and covariances in filters and estimators［C］∥IEEE Transactions on Automatic Control.USA：IEEE Press，2000：477-482.