CN103943158A - 一种消除自给能中子探测器延迟效应的方法 - Google Patents

Abstract

一种消除自给能中子探测器延迟效应的方法，写出探测器单位体积内各中间核素数量N_i(t)关于中子通量密度φ(t)的动态微分方程组以及探测电流I(t)与各中间核素数量及中子通量密度φ(t)的表达式，通过对方程组和表达式变换得到N_i(nT_s)和N_i((n-1)T_s)及I((n-1)T_s)的关系式，同时得到φ(nT_s)和N_i(nT_s)及I(nT_s)的关系式；给出N_i(0),I(0)的初始值；通过测量中子探测器的电流进行取样，得到电流采样的测量值I_m(nT_s)；对其处理得到I_s(nT_s)；得出下一步循环中探测器单位体积内各中间核素数量及当前时间段消除时间延迟效应的中子通量密度的瞬时值并输出；本发明克服了半衰期限制带来的电流信号延迟问题；引入噪声过滤使其更适合于反应堆堆芯的复杂的探测环境；无需拉氏变换、Z变换。

Description

一种消除自给能中子探测器延迟效应的方法

技术领域

本发明属于中子探测技术领域，具体涉及一种消除自给能中子探测器延迟效应的方法。

背景技术

核能是人类最具希望的未来能源。在核能反应堆中，中子通量密度是最能直观地体现反应堆功率以及反应堆状态的物理量，同时人们也是通过控制堆内中子通量密度的方式来控制反应堆。由于核反应堆的特殊性以及核反应堆安全运行的重要性，使得中子探测在反应堆内各种粒子以及射线的探测中处于至关重要的地位。

反应堆堆芯内部探测环境复杂，其对中子探测器的要求较高，要求耐高温、耐辐照，结构简单、小型化。目前常用的中子探测器按其工作机理可以分为气体探测器、半导体探测器、闪烁体探测器以及自给能探测器。其中气体探测器虽然耐高温、耐辐照，但是对于堆内高温高压的探测环境还是难以胜任。半导体探测器只适用于测量反应堆的快中子能谱，对现有的热中子反应堆应用价值不大。闪烁体探测器对高压电源的稳定性要求较高，在反应堆堆芯中难以实现。

而自给能探测器不需外加偏压、结构简单、体积小、全体固化、电子学设备简单等特性使之特别适宜于反应堆堆芯高中子通量的探测。然而当前的自给能探测器中，主要是¹⁰³Rh（铑）探测器、⁵¹V（钒）探测器和⁵⁹Co(钴）探测器，自给能探测器的探测原理如附图2所示，探测器放在堆芯中，其吸收中子后会经过几种不同途径放出电子，当电子被收集时将会在回路中产生电流；此电流强度与堆内中子通量密度有关系，即通过测量这一电流并将其经过某种处理将能够达到测量中子通量的目的。

在当前的自给能探测器中，¹⁰³Rh探测器应用较为广泛，但是由于¹⁰³Rh和⁵¹V元素在中子场中吸收中子后形成的同位素会以一定半衰期衰变，产生电子（或者伽玛射线，伽玛射线通过与物质相互作用转变为电子）形成探测器的电流信号。显然，由于半衰期的限制，电流信号不能及时反映中子通量的改变。比如，把自给能探测器突然放入一个恒定的中子场中，电流信号需要几分钟才能达到稳定值。这显然不符合反应堆堆芯中子通量的实时监控的要求。

下面以¹⁰³Rh为例具体说明在中子场中探测器电流信号的产生机制。如图2所示，在探测器电流的组成成分上，原理上应该包括三部分：1）第一部分来自¹⁰³Rh吸收中子瞬时发生(n,γ)反应，（(n,γ)反映释放出来的γ射线与物质发生光电效应或者康普顿效应产生电子。此为电流的瞬时成分；2）第二部分来自^104mRh退激成¹⁰⁴Rh放出的γ射线与物质发生光电效应或者康普顿效应所产生的电子；3）第三部分来自¹⁰⁴Rhβ衰变产生的电子。后两项因为Rh同位素半衰期的限制，属于延迟成分。

目前国内外很多学者对铑自给能探测器有过许多研究，得出了很多的成果，同时也有一些不足：1）相关文献上往往只包含了第一和第三部分，忽略了第二部分，这一部分对Rh虽不是主要成分，但是考虑到修正延迟的效果，其作用不容忽视；2）在修正延迟的时候没有考虑噪声的影响。实际上此类探测器的电流信号很小，噪声对其影响较大，所以噪声滤波是必须考虑的；3）文献中提到的方法牵扯拉氏变换、Z变换，具有一定的难度。

发明内容

为了解决上述现有技术存在的问题，本发明的目的在于提供一种消除自给能中子探测器延迟效应的方法，1）克服由于半衰期限制带来的电流信号延迟的问题；2）引入噪声过滤使得其更加适合于反应堆堆芯的复杂的探测环境；3）无需拉氏变换、Z变换；在一定程度上比较简单。

为达到以上目的，本发明采用如下技术方案：

一种消除自给能中子探测器延迟效应的方法，包括如下步骤：

步骤1：根据探测器材料在中子场中的反应物理过程画出其原理图；

步骤2：根据步骤1画出的反应机制原理图写出探测器单位体积内各中间核素数量N_i(t)关于中子通量密度φ(t)的动态微分方程组（1），写出探测电流I(t)与各中间核素数量及中子通量密度φ(t)的表达式（2）；并根据实际探测器刻度数据得到探测器单位体积内材料核素的个数、探测器对中子的瞬时成分的灵敏度、探测器单位体积内材料核素及各中间核素在中子场中的不同反应产生电流的效率参数；

dN_i(t)/dt＝F_i(N₁(t),....,N_i(t),....N_m(t),φ(t))     （1）

I(t)＝F(N₁(t),....,N_i(t),....N_m(t),φ(t))     （2）

式中：i表示第i个相关核素；m表示共有m个相关核素；t表示时间；

步骤3：将步骤2中得出的中子通量密度表达式（2）整理成表达式（3）的形式，即将中子通量密度表示为探测电流和各中间核素数量的关系，由于电流中缓发电流都是中间核素衰变产生的，则只要知道各中间核素数目即可将缓发电流消除得到能直接反应中子通量密度的瞬发电流的值，将表达式(3)带入动态微分方程组（1）中得到探测器单位体积内各中间核素数量N_i(t)关于探测电流I(t)的动态微分方程组（4）；

φ(t)＝G(N₁(t),....,N_i(t),....N_m(t),I(t))     （3）

dN_i(t)/dt＝G_i(N₁(t),....,N_i(t),....N_m(t),I(t))     （4）

步骤4：将步骤3中得到的动态微分方程组（4）中的dN_i(t)/dt写成相邻时间段内探测器单位体积内各中间核素的数量之差除以步长(N_i(nT_s)-N_i((n-1)T_s))/T_s的形式、N_i(t)写成N_i((n-1)T_s)的形式、I(t)写成I((n-1)T_s)的形式，T_s为探测器电流的取样时间间隔，即取样步长，n为步数；得到N_i(nT_s)和N_i((n-1)T_s)的关系式（5）；

N_i(nT_s)＝G_i(N₁((n-1)T_s),....,N_i((n-1)T_s),....N_m((n-1)T_s),I(nT_s))     （5）

步骤5：给出探测器单位体积内各中间核素数量N_i(t)和探测电流I(t)在时间t=0时的初始值N_i(0)和I(0)，当停堆时放入探测器并启用或者反应堆稳定功率运行时突然启用探测器时，N_i(0)=0、I(0)=0；当探测器已经启用一段时间时，初始值N_i(0)和I(0)依据不同状态取相应值；

步骤6：在实际测量中通过测量中子探测器的电流进行取样，得到电流采样的测量值I_m(nT_s)；

步骤7：对电流采样的测量值I_m(nT_s)使用改进型kalman滤波的自适应调整算法处理，降低噪声的影响并防止将突变平滑化；

所述的改进型kalman滤波的自适应调整算法由滤波和状态判断两部分组成；

（1）滤波：对电流采样的测量值I_m(nT_s)进行kalman滤波处理得到kalman滤波后电流I(nT_s)；

（2）状态判断：判断反应堆***是否处于突变状态来防止突变电流被kalman滤波***过滤掉，其判断规则为：

利用kalman滤波后电流I(nT_s)，对比I_m(nT_s)与I(nT_s)之差的绝对值，看其是否大于阈值；如果大于阈值，则表示电流信号出现比较大的波动，则判断反应堆内中子通量出现了突变，此时将电流采样的测量值I_m(nT_s)赋给输出电流I_s(nTs)并输出；如果小于阈值，则表示电流信号波动比较小，则判断此时波动是由噪声引起，则直接将kalman滤波后电流I(nTs)的值赋给输出电流I_s(nTs)并输出；

步骤8：将步骤7输出的电流值I_s(nTs)与当前步探测器单位体积内各中间核素数量N_i(nT_s)分别带入到N_i((n+1)T_s)＝H_i(N₁(nT_s),....,N_i(nT_s),....N_m(nT_s),I_s(nT_s))以及φ(nT_s)＝G(N₁(nT_s),....,N_i(nT_s),....N_m(nT_s),I_s(nT_s))中，得出下一步循环中探测器单位体积内各中间核素数量N_i((n+1)T_s)以及当前时间段消除时间延迟效应的中子通量密度φ(nT_s)的瞬时值；将φ(nT_s)输出，同时返回步骤6。

步骤7所述的阈值取噪声标准差的三倍。

与现有技术相比，本发明具有如下优点：

步骤2考虑了所有可能产生电流的反应，杜绝了由于忽略小电流导致修正延迟效果产生较大误差的后果，使修正延迟的效果更加真实。

步骤3-步骤5均为简单的数学处理及迭代过程，相对于进行拉氏变换以及Z变换处理，过程简单易操作。

步骤6考虑了噪声的影响能更真实地模拟堆内复杂的探测环境

步骤7进行了噪声滤波，采取了改进型kalman滤波算法，既能实现噪声滤波，又能在突变的时候防止突变状态被滤波导致信号延迟的发生。

附图说明

图1为本发明流程示意图。

图2为探测器探测中子相关反应过程图。

图3为归一化输出中子通量密度及归一化电流图。

图4为突变前后归一化中子通量密度及归一化电流图。

具体实施方式

以下结合附图及具体实施例，对本发明作进一步的详细描述。

下面以铑自给能中子探测器为例进行详细说明：具体方法过程如图1所示。

步骤1：根据探测器材料¹⁰³Rh在中子场中的反应物理过程画出其原理图，如图2所示。

步骤2：根据图2写出探测器单位体积内^104mRh的个数N₁(t)及¹⁰⁴Rh的个数N₂(t)关于中子通量密度φ(t)的动态微分方程（1-1）和（1-2），并写出探测电流I(t)与N₁(t)、N₂(t)及中子通量密度φ(t)的表达式(2-1)。

$\frac{{\mathrm{dN}}_1\left(t\right)}{\mathrm{dt}}={\mathrm{\σ}}_1{N}_{\mathrm{Rh}}^{103}\mathrm{\φ}\left(t\right)-{\mathrm{\λ}}_1{N}_1\left(t\right)...(1-1)$

$\frac{{\mathrm{dN}}_2\left(t\right)}{\mathrm{dt}}{\mathrm{\σ}}_2{N}_{\mathrm{Rh}}^{103}\mathrm{\φ}\left(t\right)+{\mathrm{\λ}}_1{N}_1\left(t\right)-{\mathrm{\λ}}_2{N}_2\left(t\right)...(1-2)$

$I\left(t\right)=(f_1{\mathrm{\σ}}_1{N}_{\mathrm{Rh}}^{103}+f_2{\mathrm{\σ}}_2{N}_{\mathrm{Rh}}^{103})\mathrm{\φ}\left(t\right)+j_1{\mathrm{\λ}}_1{N}_1\left(t\right)+j_2{\mathrm{\λ}}_2{N}_2\left(t\right)...(2-1)$

其中：

为探测器单位体积内¹⁰³Rh核素的个数；

N₁(t)为探测器单位体积内^104mRh的个数；

N₂(t)为探测器单位体积内¹⁰⁴Rh的个数；

σ₁为¹⁰³Rh与中子反应生成^104mRh的微观反应截面；

σ₂为¹⁰³Rh与中子反应生成¹⁰⁴Rh的微观反应截面；

λ₁为^104mRh的衰变常数；

λ₂为¹⁰⁴Rh的衰变常数；

f₁与f₂分别为探测器单位体积内产生(n,γ)反应生成^104mRh以及¹⁰⁴Rh的两个反应道及其光电效应产生瞬时电流的效率；

为探测器对中子的瞬时成分的灵敏度；

j₁为探测器单位体积内的^104mRh退激成¹⁰⁴Rh放出的γ射线与物质发生光电效应或者康普顿效应的电流产生效率；

j₂为探测器单位体积内的¹⁰⁴Rhβ衰变产生电流的效率；

φ(t)为中子通量密度；

I(t)为探测电流；

步骤3：将步骤2中得出的探测电流与中子通量密度表达式(2-1)整理成φ(t)＝G(N₁(t),....,N_i(t),....N_m(t),I(t))的形式，即将中子通量密度表示为探测电流和各中间核素数量的关系，由于电流中缓发电流都是中间核素衰变产生的，则只要知道各中间核素数目即可将缓发电流消除得到能直接反应中子通量密度的瞬发电流的值，如(3-1)；将表达式(3)带入步骤2所得到的动态微分方程组(1-1)(1-2)中得到探测器单位体积内各中间核素数量N_i(t)关于探测电流I(t)的动态微分方程组dN_i(t)/dt＝G_i(N₁(t),....,N_i(t),....N_m(t),I(t)),

如(4-1)、(4-2)。

$\mathrm{\φ}\left(t\right)=\frac{1}{f_1{\mathrm{\Σ}}_1+f_2{\mathrm{\Σ}}_2}[I\left(t\right)-j_1{\mathrm{\λ}}_1{N}_1\left(t\right)-j_2{\mathrm{\λ}}_2{N}_2\left(t\right)]...(3-1)$

$\frac{{\mathrm{dN}}_1\left(t\right)}{\mathrm{dt}}=\frac{{\mathrm{\Σ}}_1}{f_1{\mathrm{\Σ}}_1+f_2{\mathrm{\Σ}}_2}[I\left(t\right)-j_1{\mathrm{\λ}}_1{N}_1\left(t\right)-j_2{\mathrm{\λ}}_2{N}_2\left(t\right)]-{\mathrm{\λ}}_1{N}_1\left(t\right)...(4-1)$

$\frac{{\mathrm{dN}}_2\left(t\right)}{\mathrm{dt}}=\frac{{\mathrm{\Σ}}_2}{f_1{\mathrm{\Σ}}_1+f_2{\mathrm{\Σ}}_2}[I\left(t\right)-j_1{\mathrm{\λ}}_1{N}_1\left(t\right)-j_2{\mathrm{\λ}}_2{N}_2\left(t\right)]+{\mathrm{\λ}}_1{N}_1\left(t\right)...(4-2)$

其中：

为探测器单位体积内产生(n,γ)反应生成^104mRh的宏观截面；

为探测器单位体积内产生(n,γ)反应生成¹⁰⁴Rh的宏观截面；

N₁(t)为探测器单位体积内^104mRh的个数；

N₂(t)为探测器单位体积内¹⁰⁴Rh的个数；

λ₁为^104mRh的衰变常数；

λ₂为¹⁰⁴Rh的衰变常数；

f₁与f₂分别为探测器单位体积内产生(n,γ)反应生成^104mRh以及¹⁰⁴Rh的两个反应道及其光电效应产生瞬时电流的效率；

j₁为探测器单位体积内的^104mRh退激成¹⁰⁴Rh放出的γ射线与物质发生光电效应或者康普顿效应的电流产生效率；

j₂为探测器单位体积内的¹⁰⁴Rhβ衰变产生电流的效率；

φ(t)为中子通量密度；

I(t)为探测电流；

步骤4：将步骤3中得到的动态微分方程组（4-1）、（4-2）中的dN_i(t)/dt写成当前步所在时间段与上一步所在时间段探测器单位体积内各中间核素的数量之差除以步长(N_i(nT_s)-N_i((n-1)T_s))/T_s的形式、N_i(t)写成N_i((n-1)T_s)的形式、I(t)写成I((n-1)T_s)的形式，T_s为探测器电流的取样时间间隔，即取样步长，n为步数；得到N_i(nT_s)和N_i((n-1)T_s)的关系（5）；

N_i(nT_s)＝G_i(N₁((n-1)T_s),....,N_i((n-1)T_s),....N_m((n-1)T_s),I(nT_s))     （5）

具体过程如下：将步骤3中等式的微分形式表示成相邻的值除以步长的形式,如(3-1-1)、(4-1-1)、(4-2-1)；则步骤3中得到的关系式就可以表示相邻步长之间的相应参数之间的耦合关系,如(5-1)、(5-2)。

$\mathrm{\φ}\left((n-1)T_s\right)=\frac{1}{f_1{\mathrm{\Σ}}_1+f_2{\mathrm{\Σ}}_2}[I\left((n-1)T_s\right)-j_1{\mathrm{\λ}}_1{N}_1\left((n-1)T_s\right)-j_2{\mathrm{\λ}}_2{N}_2\left((n-1)T_s\right)]---(3-1-1)$

$\begin{array}{c}\frac{N_1\left(n{T}_s\right)-N_1\left((n-1)T_s\right)}{T_s}=\frac{{\mathrm{\Σ}}_1}{f_1{\mathrm{\Σ}}_1+f_2{\mathrm{\Σ}}_2}[I\left((n-1)T_s\right)-j_1{\mathrm{\λ}}_1{N}_1\left((n-1)T_s\right)\\ -j_2{\mathrm{\λ}}_2{N}_2\left((n-1)T_s\right)]-{\mathrm{\λ}}_1{N}_1\left((n-1)T_s\right)\end{array}...(4-1-1)$

$\begin{array}{c}\frac{N_1\left(n{T}_s\right)-N_1\left((n-1)T_s\right)}{T_s}=\frac{{\mathrm{\Σ}}_2}{f_1{\mathrm{\Σ}}_1+f_2{\mathrm{\Σ}}_2}[I\left((n-1)T_s\right)-j_1{\mathrm{\λ}}_1{N}_1\left((n-1)T_s\right)\\ -j_2{\mathrm{\λ}}_2{N}_2\left((n-1)T_s\right)]+{\mathrm{\λ}}_1{N}_1\left((n-1)T_s\right)-{\mathrm{\λ}}_2{N}_2\left((n-1)T_s\right)\end{array}...(4-2-1)$

得到

$\begin{array}{c}{N}_1\left(n{T}_s\right)=\frac{T_s{\mathrm{\Σ}}_1}{f_1{\mathrm{\Σ}}_1+f_2{\mathrm{\Σ}}_2}I\left((n-1)T_s\right)+[1-T_s{\mathrm{\λ}}_1(1+\frac{j_1{\mathrm{\Σ}}_1}{f_1{\mathrm{\Σ}}_1+f_2{\mathrm{\Σ}}_2})]N_1\left((n-1)T_s\right)\\ -\frac{T_s{\mathrm{\Σ}}_1{j}_2{\mathrm{\λ}}_2}{f_1{\mathrm{\Σ}}_1+f_2{\mathrm{\Σ}}_2}{N}_2\left((n-1)T_s\right)\end{array}...(5-1)$

$\begin{array}{c}{N}_2\left(n{T}_s\right)=\frac{T_s{\mathrm{\Σ}}_2}{f_1{\mathrm{\Σ}}_1+f_2{\mathrm{\Σ}}_2}I\left((n-1)T_s\right)+T_s{\mathrm{\λ}}_1(1-\frac{j_1{\mathrm{\Σ}}_2}{f_1{\mathrm{\Σ}}_1+f_2{\mathrm{\Σ}}_2})N_1\left((n-1)T_s\right)\\ -T_s{\mathrm{\λ}}_2(1+\frac{{\mathrm{\Σ}}_2{j}_2}{f_1{\mathrm{\Σ}}_1+f_2{\mathrm{\Σ}}_2})N_2\left((n-1)T_s\right)\end{array}...(5-2)$

其中：

为探测器单位体积内产生(n,γ)反应生成^104mRh的宏观截面；

为探测器单位体积内产生(n,γ)反应生成¹⁰⁴Rh的宏观截面；

T_s为取样步长；

N₁(nT_s)为当前步所在时间段探测器单位体积内^104mRh的个数；

N₂(nT_s)为当前步所在时间段探测器单位体积内¹⁰⁴Rh的个数；

I(nT_s)为当前步所在时间段的探测电流；

N₁((n-1)T_s)为上一步所在时间段探测器单位体积内^104mRh的个数；

N₂((n-1)T_s)为上一步所在时间段探测器单位体积内¹⁰⁴Rh的个数；

I((n-1)T_s)为上一步所在时间段的探测电流；

λ₁为^104mRh的衰变常数；

λ₂为¹⁰⁴Rh的衰变常数；

f₁与f₂分别为探测器单位体积内产生(n,γ)反应生成^104mRh以及¹⁰⁴Rh的两个反应道及其光电效应产生瞬时电流的效率；

j₁为探测器单位体积内的^104mRh退激成¹⁰⁴Rh放出的γ射线与物质发生光电效应或者康普顿效应的电流产生效率；

j₂为探测器单位体积内的¹⁰⁴Rhβ衰变产生电流的效率；

φ((n-1)T_s)为上一步所在时间段的中子通量密度；

则只要知道当前所在时间段^104mRh和¹⁰⁴Rh的数目以及当前所测得的电流即可得出下一个时间段的^104mRh和¹⁰⁴Rh的数目。

步骤5：给出探测器单位体积内各中间核素数量N_i(t)和探测电流I(t)在时间t=0时的初始值N_i(0)和I(0)，当停堆时放入探测器并启用或者反应堆稳定功率运行时突然启用探测器时，N_i(0)=0、I(0)=0；当探测器已经启用一段时间时，初始值N_i(0)和I(0)依据不同状态取相应值。

具体过程为：令t=0时刻参量N₁(t)、N₂(t)的值均为0，模拟停堆时放入探测器并启用或者反应堆稳定功率运行时突然启用探测器的情形。

步骤6：在实际测量中通过测量中子探测器的电流进行取样，得到电流采样的测量值I_m(nT_s)。

具体过程：取样步长间隔恒定；同时为了更好地模拟真实情况，则在取到的理论电流值的基础上加入了高斯白噪声作为实际采样的测量值，高斯白噪声为将所有频率的噪声叠加在一起形成的分布服从高斯分布的噪声。

步骤7：对电流采样的测量值I_m(nT_s)使用改进型kalman滤波的自适应调整算法处理，降低噪声的影响并防止将突变平滑化；

所述的改进型kalman滤波的自适应调整算法由滤波和状态判断两部分组成。

（1）滤波：对电流采样的测量值I_m(nT_s)进行kalman滤波处理得到kalman滤波后电流I(nT_s)。

kalman滤波的具体处理方式为：

(a)利用上一步滤波后的电流I((n-1)T_s)来预测当前步电流I_y(nT_s)；因为此处kalman***的作用是滤波，所以令当前步预测电流等于上一步滤波后的电流。

即I_y(nT_s)＝I((n-1)T_s)……………………………………………(6-1)

其中：I_y(nT_s)为当前步预测电流；I((n-1)T_s)为上一步滤波后的电流。

(b)利用上一步滤波后电流的协方差P((n-1)T_s)来预测当前步预测电流的协方差P_y(nT_s)。由于存在过程***噪声，所以当前步预测电流的协方差还需加上***过程噪声的协方差Q。

即P_y(nT_s)＝P((n-1)T_s)+Q…………………………………………(6-2)

其中：P_y(nT_s)为当前步预测电流的协方差；P((n-1)T_s)为上一步滤波后电流的协方差；Q为过程***噪声的协方差。

(c)已经得到了当前状态的预测结果I_y(nT_s)，也有当前电流采样的测量值I_m(nT_s)；可以得到kalman滤波之后的电流值I(nT_s)

I(nT_s)＝I_y(nT_s)+Kg(nT_s)×(I_m(nT_s)-I_y(nT_s))………………………(6-3)

其中I_m(nT_s)为电流采样的测量值；I_y(nT_s)为当前步预测电流；I(nT_s)为kalman

滤波之后的电流值；Kg(nT_s)为当前步的kalman增益；

且Kg(nT_s)＝P_y(nT_s)/(P_y(nT_s)+R)…………………………………(6-4)

其中P_y(nT_s)为当前步预测电流的协方差；R为由于测量***存在噪声而引入的测量***噪声的协方差。

(d)要使kalman滤波在下一步继续运行，还需要更新当前状态下滤波后电流的协方差P(nT_s)。

P(nT_s)＝(1-Kg(nT_s))×P_y(nT_s)………………………………………(6-5)

其中P_y(nT_s)为当前步预测电流的协方差；Kg(nT_s)为当前步的kalman增益；

P(nT_s)为当前状态下滤波后电流的协方差。

(2)状态判断：判断反应堆***是否处于突变状态来防止突变电流被kalman滤波***过滤掉，其判断规则为：

利用kalman滤波后得到的电流I(nT_s)，对比I_m(nT_s)与I(nT_s)之差的绝对值，看其是否大于阈值，其中阈值一般取噪声标准差的三倍，同时可以根据实际情况适当向前多比较若干步，以防止误判，并使得该方法具有更高的可靠性；如果大于阈值，则表示电流信号出现比较大的波动，则判断堆内中子通量出现了突变，此时将电流采样的测量值I_m(nT_s)赋给输出电流I_s(nTs)并输出；如果小于阈值，则表示电流信号波动比较小，则判断此时波动是由噪声引起，则直接将kalman滤波后电流I(nTs)的值赋给输出电流I_s(nTs)并输出。

步骤8：将步骤7输出的电流值I_s(nTs)以及当前步探测器单位体积内各中间核素数量N_i(nT_s)带入到N_i((n+1)T_s)＝H_i(N₁(nT_s),....,N_i(nT_s),....N_m(nT_s),I_s(nT_s))中得出下一步循环中的探测器单位体积内各中间核素数量N_i((n+1)T_s)；再带入到φ(nT_s)＝G(N₁(nT_s),....,N_i(nT_s),....N_m(nT_s),I_s(nT_s))中得到当前时间段消除时间延迟效应的中子通量密度φ(nT_s)的瞬时值。

具体代入公式如下：

$\begin{array}{c}{N}_1\left((n+1)T_s\right)=\frac{T_s{\mathrm{\Σ}}_1}{f_1{\mathrm{\Σ}}_1+f_2{\mathrm{\Σ}}_2}{I}_s\left({\mathrm{nT}}_s\right)+[1-T_s{\mathrm{\λ}}_1(1+\frac{j_1{\mathrm{\Σ}}_1}{f_1{\mathrm{\Σ}}_1+f_2{\mathrm{\Σ}}_2})]N_1\left({\mathrm{nT}}_s\right)\\ -\frac{T_s{\mathrm{\Σ}}_1{j}_2{\mathrm{\λ}}_2}{f_1{\mathrm{\Σ}}_1+f_2{\mathrm{\Σ}}_2}{N}_2\left(n{T}_s\right)\end{array}...(5-1-1)$

$\begin{array}{c}{N}_2\left((n+1)T_s\right)=\frac{T_s{\mathrm{\Σ}}_2}{f_1{\mathrm{\Σ}}_1+f_2{\mathrm{\Σ}}_2}{I}_s\left(n{T}_s\right)+T_s{\mathrm{\λ}}_1(1-\frac{j_1{\mathrm{\Σ}}_2}{f_1{\mathrm{\Σ}}_1+f_2{\mathrm{\Σ}}_2})N_1\left({\mathrm{nT}}_s\right)\\ -T_s{\mathrm{\λ}}_2(1+\frac{{\mathrm{\Σ}}_2{j}_2}{f_1{\mathrm{\Σ}}_1+f_2{\mathrm{\Σ}}_2})N_2\left({\mathrm{nT}}_s\right)\end{array}...(5-1-2)$

$\mathrm{\φ}\left(n{T}_s\right)=\frac{1}{f_1{\mathrm{\Σ}}_1+f_2{\mathrm{\Σ}}_2}[I_s\left(n{T}_s\right)-j_1{\mathrm{\λ}}_1{N}_1\left(n{T}_s\right)-j_2{\mathrm{\λ}}_2{N}_2\left(n{T}_s\right)]...(3-1-1-1)$

其中：

为探测器单位体积内产生(n,γ)反应生成^104mRh的宏观截面；

为探测器单位体积内产生(n,γ)反应生成¹⁰⁴Rh的宏观截面；

T_s为取样步长；

N₁(nT_s)为当前步所在时间段探测器单位体积内^104mRh的个数；

N₂(nT_s)为当前步所在时间段探测器单位体积内¹⁰⁴Rh的个数；

I_s(nT_s)为当前步所在时间段的输出电流；

N₁((n+1)T_s)为下一步所在时间段探测器单位体积内^104mRh的个数；

N₂((n+1)T_s)为下一步所在时间段探测器单位体积内¹⁰⁴Rh的个数；

λ₁为^104mRh的衰变常数；

λ₂为¹⁰⁴Rh的衰变常数；

f₁与f₂分别为探测器单位体积内产生(n,γ)反应生成^104mRh以及¹⁰⁴Rh的两个反应道及其光电效应产生瞬时电流的效率；

j₁为探测器单位体积内的^104mRh退激成¹⁰⁴Rh放出的γ射线与物质发生光电效应或者康普顿效应的电流产生效率；

j₂为探测器单位体积内的¹⁰⁴Rhβ衰变产生电流的效率；

φ(nT_s)为当前步所在时间段的中子通量密度；

取值说明：

^104mRh的半衰期为258.0秒，则^104mRh的衰变常数

¹⁰⁴Rh的半衰期为42.4秒；则¹⁰⁴Rh的衰变常数

¹⁰³Rh与中子反应生成^104mRh的微观反应截面σ₁=1.1×10^-23cm²；

¹⁰³Rh与中子反应生成¹⁰⁴Rh的微观反应截面σ₂=1.35×10^-22cm²。

由于所考察的时间长度远小于探测器的寿命长度，所以假设在考察时间段内探测器单位体积内¹⁰³Rh核素的个数保持不变，即为常数；由于铑的密度为12.7g/cm³,则单位体积内含有¹⁰³Rh的个数这里ρ=12.7g/cm³为铑的密度；M=103g/mol为铑的摩尔质量，N_A=6.02×10²³mol^-1为阿伏伽德罗常数。则探测器单位体积内产生(n,γ)反应生成^104mRh以及¹⁰⁴Rh的宏观截面分别为

铑自给能中子探测器的中子的灵敏度典型值为3.6×10^-20A/(cm*s)，其中瞬时成分占比约为5%～15%，此处取10%，即 j₁与j₂的值将在实际应用过程中根据探测器刻度值进行调谐确定;在这里暂且取j₁=2.9897×10^-22A/(cm³*s),j₂=2.9897×10^-21A/(cm³*s)；采样步长T_s可以根据具体情况进行优化，这里暂取为0.1秒。

在kalman滤波***中Q与R的取值因***不同而不同,对于P(0)来说如果取值为0的话就表示kalman***是最优的则会使算法不能收敛，且P(n)将随着n的增加而逐渐收敛于0，则可以取P(0)为一个随意的正整数，在这里取Q=0.0001，R=0.0004，I(0)=0,P(0)=10。

应用实例：

为了突出消除时间延迟的效果，考查将铑探测器突然放入一个稳定的中子场的过程，通常认为延迟时间是达到稳态值的95%所需的时间。以反应堆启堆过程为例，模拟反应堆中中子通量密度φ=0(cm²·s)^-1在t=100秒时到突然变为φ=5×10¹³(cm²·s)^-1的过程；整个过程前100秒中子通量密度φ=0(cm²·s)^-1，在第100秒时中子通量密度突变为φ=5×10¹³(cm²·s)^-1，并一直维持。在此模拟过程中，各相关参数参照具体实施方式以及取值说明部分中所述情形进行取值，电流取样时间步长T_s=0.1秒,同时为了更好地模拟真实情况，所以在取到的理论电流值的基础上加入了均值为0、标准差σ=1×10^-8A（约占稳态电流的1%）的高斯白噪声的影响作为实际所测电流；判断噪声状态还是突变状态的阈值为3σ。

消除延迟效应的结果如图3及图4所示，包括归一化中子通量密度和归一化电流随时间变化的曲线。归一化中子通量密度为用消除延迟算法得出的中子通量密度与突变后中子通量密度的比值（归一化分母取5×10¹³(cm²·s)^-1）；归一化电流为测得的存在噪声的电流的值与电流最终达到的稳定值的比值；在图3中可以看出电流需要三百秒以上的时间才能达到稳定，即如果仅仅依据电流数据来显示中子通量密度的话将存在很大的延迟，如果采用这种方法的话，对反应堆中子通量密度的监控也失去了意义；与之相对，采用消除延迟算法得出来的中子通量密度与实际中子通量密度吻合得非常好，同时将图3中突变前后的过程截取出来作为图4，在图4中可以看出延迟被控制在短短0.1秒以内(即使为了更高的可靠性，考虑多向前预测几步来判断噪声还是突变，也能够很好的控制在1秒之内)，可以实现了反应堆中子通量密度的实时监控，更利于反应堆安全控制。

虽然已参照典型实施例描述了本发明，但应该理解，所用的术语是说明和示例性、而非限制性的术语。由于本发明能够以多种形式具体实施而不脱离发明的精神和实质，所以应当理解，上述实施例不限于任何前述的细节，而应该在随附权利要求所限定的精神和范围内广泛地解释，因此落入权利要求或其等效范围内的全部变化和改型都应为随附权利要求所覆盖。

Claims (2)

1.一种消除自给能中子探测器延迟效应的方法，其特征在于：包括如下步骤：

步骤1：根据探测器材料在中子场中的反应物理过程画出其原理图；

步骤2：根据步骤1画出的反应机制原理图写出探测器单位体积内各中间核素数量N_i(t)关于中子通量密度φ(t)的动态微分方程组（1），写出探测电流I(t)与各中间核素数量及中子通量密度φ(t)的表达式（2）；并根据实际探测器刻度数据得到探测器单位体积内材料核素的个数、探测器对中子的瞬时成分的灵敏度、探测器单位体积内材料核素及各中间核素在中子场中的不同反应产生电流的效率参数；

dN_i(t)/dt＝F_i(N₁(t),....,N_i(t),....N_m(t),φ(t))     （1）

I(t)＝F(N₁(t),....,N_i(t),....N_m(t),φ(t))     （2）

式中：i表示第i个相关核素；m表示共有m个相关核素；t表示时间；

步骤3：将步骤2中得出的中子通量密度表达式（2）整理成表达式（3）的形式，即将中子通量密度表示为探测电流和各中间核素数量的关系，由于电流中缓发电流都是中间核素衰变产生的，则只要知道各中间核素数目即可将缓发电流消除得到能直接反应中子通量密度的瞬发电流的值，将表达式(3)带入动态微分方程组（1）中得到探测器单位体积内各中间核素数量N_i(t)关于探测电流I(t)的动态微分方程组（4）；

φ(t)＝G(N₁(t),....,N_i(t),....N_m(t),I(t))     （3）

dN_i(t)/dt＝G_i(N₁(t),....,N_i(t),....N_m(t),I(t))     （4）

步骤4：将步骤3中得到的动态微分方程组（4）中的dN_i(t)/dt写成相邻采样点时间段内探测器单位体积内各中间核素的数量之差除以步长(N_i(nT_s)-N_i((n-1)T_s))/T_s的形式、N_i(t)写成N_i((n-1)T_s)的形式、I(t)写成I((n-1)T_s)的形式，T_s为探测器电流的取样时间间隔，即取样步长，n为步数；得到N_i(nT_s)和N_i((n-1)T_s)的关系式（5）；

N_i(nT_s)＝G_i(N₁((n-1)T_s),....,N_i((n-1)T_s),....N_m((n-1)T_s),I(nT_s))     （5）

步骤5：给出探测器单位体积内各中间核素数量N_i(t)和探测电流I(t)在时间t=0时的初始值N_i(0)和I(0)，当停堆时放入探测器并启用或者反应堆稳定功率运行时突然启用探测器时，N_i(0)=0、I(0)=0；当探测器已经启用一段时间时，初始值N_i(0)和I(0)依据不同状态取相应值；

步骤6：在实际测量中通过测量中子探测器的电流进行取样，得到电流采样的测量值I_m(nT_s)；

步骤7：对电流采样的测量值I_m(nT_s)使用改进型kalman滤波的自适应调整算法处理，降低噪声的影响并防止将突变平滑化；

所述的改进型kalman滤波的自适应调整算法由滤波和状态判断两部分组成；

（1）滤波：对电流采样的测量值I_m(nT_s)进行kalman滤波处理得到kalman滤波后电流I(nT_s)；

（2）状态判断：判断反应堆***是否处于突变状态来防止突变电流被kalman滤波***过滤掉，其判断规则为：

利用kalman滤波后电流I(nT_s)，对比I_m(nT_s)与I(nT_s)之差的绝对值，看其是否大于阈值；如果大于阈值，则表示电流信号出现比较大的波动，则判断反应堆内中子通量出现了突变，此时将电流采样的测量值I_m(nT_s)赋给输出电流I_s(nTs)并输出；如果小于阈值，则表示电流信号波动比较小，则判断此时波动是由噪声引起，则直接将kalman滤波后电流I(nTs)的值赋给输出电流I_s(nTs)并输出；

步骤8：将步骤7输出的电流值I_s(nTs)与当前步探测器单位体积内各中间核素数量N_i(nT_s)分别带入到N_i((n+1)T_s)＝H_i(N₁(nT_s),....,N_i(nT_s),....N_m(nT_s),I_s(nT_s))以及φ(nT_s)＝G(N₁(nT_s),....,N_i(nT_s),....N_m(nT_s),I_s(nT_s))中，得出下一步循环中探测器单位体积内各中间核素数量N_i((n+1)T_s)以及当前时间段消除时间延迟效应的中子通量密度φ(nT_s)的瞬时值；将φ(nT_s)输出，同时返回步骤6。

2.根据权利要求1所述的方法，其特征在于：步骤7所述的阈值取噪声标准差的三倍。