一种用于核辐射探测的数字滤波器优化设计方法

技术领域

本发明涉及核辐射探测领域，特别涉及一种用于核辐射探测中的数字滤波器优化 设计方法。

背景技术

盖革-米勒计数管(Geiger-Müllercounter，以下简称GM计数管)是核辐射探测领 域中常用的探测器，但其受到固有死时间与恢复时间的影响，线性范围与测量上限受到限 制。为了提高GM计数管的量程，近年来提出了一种通过控制GM计数管工作电压实现宽量程 测量的方法，文献《一种新的GM计数管Time-To-Count测量方法》(核电子学与探测技术， 2002年第03期)、《Time-to-Count技术测量辐射强度的原理》(核电子学与探测技术，2006年 第01期)、《Time-to-Count方法扩展GM计数管量程的实验研究》(原子能科学技术，2013年04 期)等均报道了此方法。通过GM计数管的工作电压在正常计数坪电压与低于起始电压下切 换，先计算GM计数管加上计数坪电压时间内的平均计数率(即GM计数管具有探测能力时间 内的平均计数率)，再由平均计数率来计算剂量率。该方法的关键在于计算平均等待时间 (即GM计数管加计数坪电压的平均时间)，但由于测量平均等待时间服从负指数分布，欲取 得较精确的平均值，需采用大量的测量结果进行平均或加权平均，这就要求测量系统必须 具有较大的存储量与较快的计算速度。文献《单管宽量程剂量仪关键技术分析》(第十六届 全国核电子学与核探测技术学术年会论文集(上册)，2012)从数字信号处理的角度，将平均 过程看作测量信号通过低通滤波器的结果，采用低通滤波器处理等待时间，使等待时间样 本通过滤波器后具有较小的方差，并建议通过改变设计指标的方式设计IIR数字低通滤波 器。虽然这种方法设计的低通滤波器可满足测量精度要求，但未考虑滤波器响应时间长短， 导致测量结果不能即时反映辐射强度或剂量率的变化，从而使得测量结果明显滞后于辐射 场的变化。

发明内容

为解决上述技术问题，本发明提供一种用于核辐射探测的数字滤波器优化设计方 法，该方法能在满足测量精度的前提下减少存储空间、提高计算效率与反应速度，达到减轻 测量装置的存储与计算压力、提升响应速度的目的。

本发明的基本构思是：将现有的平均或加权平均过程看作测量所得等待时间通过 FIR滤波器的结果，并采用IIR滤波器代替现有的平均或加权平均处理等待时间，通过粒子 群优化算法对IIR滤波器的系数[a₁,a₂,…,a_q]、[b₀,b₁,b₂,…,b_q]进行优化设计，使设计的 IIR滤波器在一定测量精度(由滤波器噪声增益表征)条件下，具有最快的响应速度。

本发明技术方案具体包括如下步骤：

步骤1：初始化参数

设置滤波器噪声增益g、滤波器阶数q、修正因子f(0<f<1)、粒子群规模M、最大迭代 次数K、权重c₀、c₁、c₂、粒子最大速度V_max，随机构造M个向量x_i＝{x_i1,x_i2,…,x_i(2q-1)}分别代 表M个粒子当前位置，设置迄今为止全局最优位置p_g＝{p_g1,p_g2,…,p_g(2q-1)}为任意值，对应 的滤波器系数a_g＝{a_g0,a_g1,a_g2,…,a_gq}、b_g＝{b_g0,b_g1,b_g2,…,b_gq}为任意值，响应时间T_g为 无穷大，各粒子的最优位置p_i＝{p_i1,p_i2,…,p_i(2q-1)}为任意值(其中i为粒子序号)，对应的 滤波器系数a_i＝{a_i0,a_i1,a_i2,…,a_iq}、b_i＝{b_i0,b_i1,b_i2,…,b_iq}为任意值，各最优位置对应 的响应时间T_i为无穷大；

步骤2：采用粒子群优化算法进行滤波器优化设计；

步骤2.1：令当前迭代次数k＝1；

步骤2.2：对各粒了进行设计；

步骤2.2.1：令当前粒子序号i＝1；

步骤2.2.2：利用式(1)更新粒子x_i：

v_id＝c₀v_id+c₁r₁(p_id-x_id)+c₂r₂(p_gd-x_id)

$\{\begin{array}{ccc}{v}_{id}=v_{\max },& if& v_{id}>v_{\max }\\ v_{id}=-v_{\max },& if& v_{id}<-v_{\max }\end{array}---\left(1\right)$

x_id＝x_id+v_id

式中，x_id代表粒子x_i的第d维变量，v_id代表粒子x_i的第d维速度；

步骤2.2.3：对当前粒子i，令a_j＝x_ij(j＝1,2,…,q)构造滤波器系数[1,a₁,a₂,…, a_q]，并计算滤波器的极点，若滤波器某极点位于单位圆外，则将该极点的倒数作为滤波器 的极点，并重新计算滤波器系数[1,a₁,a₂,…,a_q]，同时令x_ij＝a_j(j＝1,2,…,q)更新粒子的 前q维；

步骤2.2.4：对当前粒子i，令b_j＝x_i(q+j+1)(j＝0,1,2,…,q-2)，b_q-1＝0、b_q＝0，构造 滤波器系数[b₀,b₁,b₂,…,b_q]，根据系数[1,a₁,a₂,…,a_q]、[b₀,b₁,…,b_q]由式(2)计算滤波 器的单位响应h₁(n)：

$h_1\left(n\right)=-\underset{i=1}{\overset{q}{\Sigma }}{a}_i{h}_1(n-i)+\underset{i=0}{\overset{q}{\Sigma }}{b}_i\delta (n-i),(n=0,1,2,...,q)---\left(2\right)$

式中，δ(n)为单位函数，当n＝0时，其值为1，否则其值为0；

步骤2.2.5：按式计算滤波器最小噪声增益g_min：

$g_{\min }=\underset{n=0}{\overset{q}{\Sigma }}{h}_1{\left(n\right)}^2---\left(3\right)$

步骤2.2.6：如果g_min小于等于g，转步骤2.2.7，否则按式(4)更新粒子i的后q-1维， 并转步骤2.2.4；

$x_{i(q+j+1)}=x_{i(q+j+1)}f\sqrt{\frac{g}{g_{\min }}},(j=0,1,2,...,q-2)---\left(4\right);$

步骤2.2.7：对当前粒子i，由式(5)构建矩阵A₁、A₂，按式(6)计算向量c：

c＝[100…0](A₁+A₂)^-1(6)；

步骤2.2.8：由式(7)～(9)构建矩阵H₁、H₂、H₃，按式(10)将向量c拆为向量c₁、c₂，按 式(11)构建向量b₁；

$H_2=\left(\begin{array}{cc}0& h_1\left(0\right)\\ h_1\left(0\right)& h_1\left(1\right)\\ h_1\left(1\right)& h_1\left(2\right)\\ h_1\left(2\right)& h_1\left(3\right)\\ .& .\\ .& .\\ .& .\\ h_1(q-3)& h_1(q-2)\end{array}\right)---\left(8\right)$

$H_3=\left(\begin{array}{cccccc}0& h_1\left(0\right)& h_1\left(1\right)& h_1\left(2\right)& \mathrm{...}& h_1(q-3)\\ h_1\left(0\right)& h_1\left(1\right)& h_1\left(2\right)& h_1\left(3\right)& \mathrm{...}& h_1(q-2)\end{array}\right)---\left(9\right)$

$\left(\begin{array}{c}c=\left(\begin{array}{cccccc}{c}_0& c_1& c_2& c_3& \mathrm{...}& c_q\end{array}\right)\\ c_1=\left(\begin{array}{cccccc}{c}_0& c_1& c_2& c_3& \mathrm{...}& c_{q-2}\end{array}\right)\\ c_2=\left(\begin{array}{cc}{c}_{q-1}& c_q\end{array}\right)\end{array}\right)---\left(10\right)$

$b_1=\left(\begin{array}{c}{b}_0\\ b_1\\ .\\ .\\ .\\ b_{q-2}\end{array}\right)---\left(11\right);$

步骤2.2.9：由式(12)计算参数g₀、ch₁与ch₂，由式(13)计算参数u、s：

$\{\begin{array}{c}{g}_0=c_1{H}_1{b}_1+c_2{H}_3{b}_1\\ \left(\begin{array}{cc}{\mathrm{ch}}_1& {\mathrm{ch}}_2\end{array}\right)=c_1{H}_2\end{array}---\left(12\right)$

$\{\begin{array}{c}u=(\underset{i=0}{\overset{q}{\Sigma }}{a}_i-\underset{i=0}{\overset{q-2}{\Sigma }}{b}_i)\\ s=u+h_1\left(q\right)\end{array}---\left(13\right);$

步骤2.2.10：由式(14)计算一元二次方程系数A、B、C：

$\{\begin{array}{c}A=(a_1+2)c_q-c_{q-1}\\ B=ch-{\mathrm{ch}}_2+d[c_{q-1}-c_q(a_1+1)]+\\ (c_q-c_{q-1})h_1(q-1)+c_{q-1}{h}_1(q-2)-c_{q}s\\ C=u[{\mathrm{ch}}_2+c_{q-1}{h}_1(q-1)+c_{q}s]+g_0-g\end{array}---\left(14\right);$

步骤2.2.11：计算一元二次方程At²+Bt+C＝0的根，如果两根为实根，分别将两根 作为滤波器系数b_q-1，并按式(15)计算滤波器系数b_q，转步骤2.2.12；否则转步骤2.2.15；

$b_q=\underset{i=0}{\overset{q}{\Sigma }}{a}_i-\underset{i=0}{\overset{q-2}{\Sigma }}{b}_i-b_{q-1}---\left(15\right);$

步骤2.2.12：根据步骤2.2.2所得滤波器系数a₁,a₂,…,a_q、b₀,b₁,…,b_q-2，结合步骤 2.2.11计算所得b_q-1与b_q的两组值构建两个滤波器，按式(16)分别计算两滤波器的阶跃响 应，并以阶跃响应s(n)达到0.9时的n值，确定响应时间tn₁、tn₂，选响应时间短的滤波器作为 设计结果，取得滤波器系数[1,a₁,a₂,…,a_q]、[b₀,b₁,…,b_q]与响应时间tn；

$s\left(n\right)=-\underset{l=1}{\overset{q}{\Sigma }}{a}_{l}s(n-l)+\underset{l=1}{\overset{q}{\Sigma }}{b}_lϵ(n-l)---\left(16\right)$

其中，单位阶跃函数定义如下：

$ϵ\left(n\right)=\left(\begin{array}{cc}1,& n\ge 0\\ 0,& n<0\end{array}\right)---\left(17\right);$

步骤2.2.13：对于第i个粒子，将响应时间tn与其历史最优位置对应的响应时间T_i进行比较，如果tn小于T_i，令T_i＝tn、p_i＝x_i、a_i＝[1,a₁,a₂,…,a_q]、b_i＝[b₀,b₁,…,b_q]；

步骤2.2.14：将T_i与T_g进行比较，如果T_i小于T_g，令T_g＝T_i、p_g＝p_i、a_g＝a_i、b_g＝b_i；

步骤2.2.15：令i＝i+1，如果i<M，那么转步骤2.2.2，否则转步骤2.3；

步骤2.3：令k＝k+1，如果k<K，转步骤2.2.1，否则结束优化设计；

步骤3：将全局最优位置对应的滤波器输出，作为优化设计结果。

本发明的优势在于：

(1)可根据需要指定滤波器的噪声增益，控制测量误差；

(2)优化设计出的滤波器与传统平均或加权平均算法相比，在相同噪声增益情况 计算量很小；

(3)优化设计出的滤波器与现有滤波器设计工具设计的滤波器相比，在相同噪声 增益情况下具有更快的响应速度。

附图说明

图1数字滤波器优化设计方法程序流程图

图2传统方法设计滤波器与本方法设计滤波器阶跃响应比较图

具体实施方式

设核辐射测量装置要求测量所得辐射强度(剂量率)的误差小于10％，根据3σ法 则，要求平均等待时间测均方根误差小于1/30，那么平均时所需的测量数据应不少于900， 此处选择1000，即平均处理后，噪声功率变为原来的1/1000。按照文献《单管宽量程剂量仪 关键技术分析》(第十六届全国核电子学与核探测技术学术年会论文集(上册)，2012)修改 滤波器指标的方法，用传统滤波器设计方法，设计椭圆低通滤波器EF。设置滤波器EF的通带 归一化宽度为0.001，通带衰减为1dB，阻带归一化宽度为0.003，阻带衰减为40dB，设计出的 滤波器为3阶滤波器，该滤波器的噪声增益为0.00102，与平均方法噪声增益相当。采用1000 个测量数据进行平均、加权平均及采用滤波器EF处理三者的存储量、计算量对比见表1所 示。

表1平均、加权平均、滤波器EF处理三者的存储量、计算量对比

从表1可以看出，滤波器EF的计算量与存储量均明显少于平均或加权平均方法，也 就是说，在平滑性能相近情况下，核探测中选用IIR滤波器进行平滑处理可降低测量系统存 储与计算要求。滤波器EF的系数见表2所示。

本发明以上述设计低通滤波器EF为参照，设置噪声增益为0.00102(与滤波器EF相 同)，滤波器阶数同样为3阶。由于两滤波器阶数相同，其计算量与存储量相同，二者噪声增 益相同，其测量精度也相同，对比二者响应速度便能分析二者的优劣。具体步骤如下：

步骤1：初始化参数

设置滤波器噪声增益g＝0.00102、滤波器阶数q＝3、修正因子f＝0.9、粒子群规模 M＝100、最大迭代次数K＝10000、c₀＝1、c₁＝2、c₂＝2、粒子最大速度V_max＝100，随机构造M个 向量x_i＝{x_i1,x_i2,…,x_i(2q-1)}分别代表M个粒子当前位置，设置迄今为止全局最优位置p_g＝ {p_g1,p_g2,…,p_g(2q-1)}为任意值，对应的滤波器系数a_g＝{a_g0,a_g1,a_g2,…,a_gq}、b_g＝{b_g0,b_g1, b_g2,…,b_gq}为任意值，响应时间T_g为无穷大，各粒子的最优位置p_i＝{p_i1,p_i2,…,p_i(2q-1)}为 任意值(其中i为粒子序号)，对应的滤波器系数a_i＝{a_i0,a_i1,a_i2,…,a_iq}、b_i＝{b_i0,b_i1, b_i2,…,b_iq}为任意值，各最优位置对应的响应时间T_i为无穷大；

步骤2：采用粒子群优化算法进行滤波器优化设计；

步骤2.1：令当前迭代次数k＝1；

步骤2.2：对各粒了进行设计；

步骤2.2.1：令当前粒子序号i＝1；

步骤2.2.2：利用式(1)更新粒子x_i：

v_id＝c₀v_id+c₁r₁(p_id-x_id)+c₂r₂(p_gd-x_id)

$\{\begin{array}{ccc}{v}_{id}=v_{\max },& if& v_{id}>v_{\max }\\ v_{id}=-v_{\max },& if& v_{id}<-v_{\max }\end{array}---\left(1\right)$

x_id＝x_id+v_id

式中，x_id代表粒子x_i的第d维变量，v_id代表粒子x_i的第d维速度；

步骤2.2.3：对当前粒子i，令a_j＝x_ij(j＝1,2,…,q)构造滤波器系数[1,a₁,a₂,…, a_q]，并计算滤波器的极点，若滤波器某极点位于单位圆外，则将该极点的倒数作为滤波器 的极点，并重新计算滤波器系数[1,a₁,a₂,…,a_q]，同时令x_ij＝a_j(j＝1,2,…,q)更新粒子的 前q维；

步骤2.2.4：对当前粒子i，令b_j＝x_i(q+j+1)(j＝0,1,2,…,q-2)，b_q-1＝0、b_q＝0，构造 滤波器系数[b₀,b₁,b₂,…,b_q]，根据系数[1,a₁,a₂,…,a_q]、[b₀,b₁,…,b_q]由式(2)计算滤波 器的单位响应h₁(n)：

$h_1\left(n\right)=-\underset{i=1}{\overset{q}{\Sigma }}{a}_i{h}_1(n-i)+\underset{i=0}{\overset{q}{\Sigma }}{b}_i\delta (n-i),(n=0,1,2,...,q)---\left(2\right)$

式中，δ(n)为单位函数，当n＝0时，其值为1，否则其值为0；

步骤2.2.5：按式计算滤波器最小噪声增益g_min：

$g_{\min }=\underset{n=0}{\overset{q}{\Sigma }}{h}_1{\left(n\right)}^2---\left(3\right)$

步骤2.2.6：如果g_min小于等于g，转步骤2.2.7，否则按式(4)更新粒子i的后q-1维， 并转步骤2.2.4；

$x_{i(q+j+1)}=x_{i(q+j+1)}f\sqrt{\frac{g}{g_{\min }}},(j=0,1,2,...,q-2)---\left(4\right);$

步骤2.2.7：对当前粒子i，由式(5)构建矩阵A₁、A₂，按式(6)计算向量c：

c＝[100…0](A₁+A₂)^-1(6)；

步骤2.2.8：由式(7)～(9)构建矩阵H₁、H₂、H₃，按式(10)将向量c拆为向量c₁、c₂，按 式(11)构建向量b₁；

$H_2=\left(\begin{array}{cc}0& h_1\left(0\right)\\ h_1\left(0\right)& h_1\left(1\right)\\ h_1\left(1\right)& h_1\left(2\right)\\ h_1\left(2\right)& h_1\left(3\right)\\ .& .\\ .& .\\ .& .\\ h_1(q-3)& h_1(q-2)\end{array}\right)---\left(8\right)$

$H_3=\left(\begin{array}{cccccc}0& h_1\left(0\right)& h_1\left(1\right)& h_1\left(2\right)& \mathrm{...}& h_1(q-3)\\ h_1\left(0\right)& h_1\left(1\right)& h_1\left(2\right)& h_1\left(3\right)& \mathrm{...}& h_1(q-2)\end{array}\right)---\left(9\right)$

$\left(\begin{array}{c}c=\left(\begin{array}{cccccc}{c}_0& c_1& c_2& c_3& \mathrm{...}& c_q\end{array}\right)\\ c_1=\left(\begin{array}{cccccc}{c}_0& c_1& c_2& c_3& \mathrm{...}& c_{q-2}\end{array}\right)\\ c_2=\left(\begin{array}{cc}{c}_{q-1}& c_q\end{array}\right)\end{array}\right)---\left(10\right)$

$b_1=\left(\begin{array}{c}{b}_0\\ b_1\\ .\\ .\\ .\\ b_{q-2}\end{array}\right)---\left(11\right);$

步骤2.2.9：由式(12)计算参数g₀、ch₁与ch₂，由式(13)计算参数u、s：

$\{\begin{array}{c}{g}_0=c_1{H}_1{b}_1+c_2{H}_3{b}_1\\ \left(\begin{array}{cc}{\mathrm{ch}}_1& {\mathrm{ch}}_2\end{array}\right)=c_1{H}_2\end{array}---\left(12\right)$

$\{\begin{array}{c}u=(\underset{i=0}{\overset{q}{\Sigma }}{a}_i-\underset{i=0}{\overset{q-2}{\Sigma }}{b}_i)\\ s=u+h_1\left(q\right)\end{array}---\left(13\right);$

步骤2.2.10：由式(14)计算一元二次方程系数A、B、C：

$\{\begin{array}{c}A=(a_1+2)c_q-c_{q-1}\\ B=ch-{\mathrm{ch}}_2+d[c_{q-1}-c_q(a_1+1)]+\\ (c_q-c_{q-1})h_1(q-1)+c_{q-1}{h}_1(q-2)-c_{q}s\\ C=u[{\mathrm{ch}}_2+c_{q-1}{h}_1(q-1)+c_{q}s]+g_0-g\end{array}---\left(14\right);$

步骤2.2.11：计算一元二次方程At²+Bt+C＝0的根，如果两根为实根，分别将两根 作为滤波器系数b_q-1，并按式(15)计算滤波器系数b_q，转步骤2.2.12；否则转步骤2.2.15；

$b_q=\underset{i=0}{\overset{q}{\Sigma }}{a}_i-\underset{i=0}{\overset{q-2}{\Sigma }}{b}_i-b_{q-1}---\left(15\right);$

步骤2.2.12：根据步骤2.2.2所得滤波器系数a₁,a₂,…,a_q、b₀,b₁,…,b_q-2，结合步骤 2.2.11计算所得b_q-1与b_q的两组值构建两个滤波器，按式(16)分别计算两滤波器的阶跃响 应，并以阶跃响应s(n)达到0.9时的n值，确定响应时间tn₁、tn₂，选响应时间短的滤波器作为 设计结果，取得滤波器系数[1,a₁,a₂,…,a_q]、[b₀,b₁,…,b_q]与响应时间tn：

$s\left(n\right)=-\underset{l=1}{\overset{q}{\Sigma }}{a}_{l}s(n-l)+\underset{l=1}{\overset{q}{\Sigma }}{b}_lϵ(n-l)---\left(16\right)$

其中，单位阶跃函数定义如下：

$ϵ\left(n\right)=\left(\begin{array}{cc}1,& n\ge 0\\ 0,& n<0\end{array}\right)---\left(17\right);$

步骤2.2.13：对于第i个粒子，将响应时间tn与其历史最优位置对应的响应时间T_i进行比较，如果tn小于T_i，令T_i＝tn、p_i＝x_i、a_i＝[1,a₁,a₂,…,a_q]、b_i＝[b₀,b₁,…,b_q]；

步骤2.2.14：将T_i与T_g进行比较，如果T_i小于T_g，令T_g＝T_i、p_g＝p_i、a_g＝a_i、b_g＝b_i；

步骤2.2.15：令i＝i+1，如果i<M，那么转步骤2.2.2，否则转步骤2.3；

步骤2.3：令k＝k+1，如果k<K，转步骤2.2.1，否则结束优化设计；

步骤3：将全局最优位置对应的滤波器输出，作为优化设计结果。

由本发明设计滤波器与滤波器EF的阶跃响应见附图2所示，两滤波器的响应时间 见表2所示，两滤波器的系数见表3所示。

表2实施例设计滤波器与滤波器EF响应时间

表3实施例设计滤波器与滤波器EF系数

从上述分析可知，与传统的平均或加权平均相比，本发明设计的滤波器可有效减 少核辐射测量装置的计算量与存储量。从表2及附图2可以看出，在相同噪声增益情况下，传 统滤波器设计方法设计的低通滤波器响应速度低于本发明设计的低通滤波器，且传统方法 设计滤波器的阶跃响应还存在过冲现象，将其作为辐射探测装置中的滤波器使用效果较 差，而本发明设计的低通滤波器响应速度快，阶跃响应没有过冲现象，适合作为辐射探测装 置中的滤波器。