一种用于核辐射探测的数字滤波器优化设计方法
【技术领域】
[0001] 本发明涉及核辐射探测领域，特别涉及一种用于核辐射探测中的数字滤波器优化 设计方法。
【背景技术】
[0002] 盖革-米勒计数管(Geiger-Μ??Ier counter，以下简称GM计数管)是核福射探测领 域中常用的探测器，但其受到固有死时间与恢复时间的影响，线性范围与测量上限受到限 制。为了提高GM计数管的量程，近年来提出了一种通过控制GM计数管工作电压实现宽量程 测量的方法，文献《一种新的GM计数管Time-To-Count测量方法》（核电子学与探测技术， 2002年第03期）、《Time-t〇-Count技术测量辐射强度的原理》(核电子学与探测技术，2006年 第01期）、《Time-t 〇-C〇Unt方法扩展GM计数管量程的实验研究》（原子能科学技术，2013年04 期）等均报道了此方法。通过GM计数管的工作电压在正常计数坪电压与低于起始电压下切 换，先计算GM计数管加上计数坪电压时间内的平均计数率（即GM计数管具有探测能力时间 内的平均计数率），再由平均计数率来计算剂量率。该方法的关键在于计算平均等待时间 (即GM计数管加计数坪电压的平均时间），但由于测量平均等待时间服从负指数分布，欲取 得较精确的平均值，需采用大量的测量结果进行平均或加权平均，这就要求测量系统必须 具有较大的存储量与较快的计算速度。文献《单管宽量程剂量仪关键技术分析》(第十六届 全国核电子学与核探测技术学术年会论文集(上册），2012)从数字信号处理的角度，将平均 过程看作测量信号通过低通滤波器的结果，采用低通滤波器处理等待时间，使等待时间样 本通过滤波器后具有较小的方差，并建议通过改变设计指标的方式设计IIR数字低通滤波 器。虽然这种方法设计的低通滤波器可满足测量精度要求，但未考虑滤波器响应时间长短， 导致测量结果不能即时反映辐射强度或剂量率的变化，从而使得测量结果明显滞后于辐射 场的变化。
【发明内容】
[0003] 为解决上述技术问题，本发明提供一种用于核辐射探测的数字滤波器优化设计方 法，该方法能在满足测量精度的前提下减少存储空间、提高计算效率与反应速度，达到减轻 测量装置的存储与计算压力、提升响应速度的目的。
[0004] 本发明的基本构思是:将现有的平均或加权平均过程看作测量所得等待时间通过 FIR滤波器的结果，并采用IIR滤波器代替现有的平均或加权平均处理等待时间，通过粒子 群优化算法对IIR滤波器的系数[ &1，&2，···，&(1] JboAib2,…，!^]进行优化设计，使设计的 IIR滤波器在一定测量精度（由滤波器噪声增益表征)条件下，具有最快的响应速度。
[0005] 本发明技术方案具体包括如下步骤：
[0006]步骤1:初始化参数
[0007]设置滤波器噪声增益g、滤波器阶数q、修正因子f (0〈f〈l)、粒子群规模M、最大迭代 次数K、权重CQ、C1、C2、粒子最大速度Vmax，随机构造 M个向量Xi= {xil，Xi2,…，Xi(2q-1)}分别代 表M个粒子当前位置，设置迄今为止全局最优位置Pg= {pgl，Pg2,…，pg(2q-1)}为任意值，对应 的滤波器系数ag= {ag〇，agi，ag2，…，agq}、bg= {bg〇，bgi，bg2，…，bgq}为任思值，响应时间Tg为 无穷大，各粒子的最优位置Pi= {pil，Pi2,…，Pi(2q-1)}为任意值(其中i为粒子序号），对应的 滤波器系数ai = {ai〇，aii，ai2，…，aiq}、bi = {bi〇，bii，bi2，…，biq}为任意值，各最优位置对应 的响应时间Ti为无穷大；
[0008] 步骤2:采用粒子群优化算法进行滤波器优化设计；
[0009] 步骤2.1:令当前迭代次数k = l;
[0010] 步骤2.2:对各粒了进行设计；
[0011] 步骤2.2.1:令当前粒子序号i = l;
[0012] 步骤2.2.2:利用式（1)更新粒子Xi:
[0013] ViH = CdViH+Ciri (DiH-XiHHcm(Pgd-Xid)
[0014](1)
[0015] Xid = Xid+Vid
[0016] 式中，Xid代表粒子Xi的第d维变量，Vid代表粒子Xi的第d维速度；
[0017] 步骤2.2.3:对当前粒子i，令aj = xij( j = l，2，···，q)构造滤波器系数[I，ai，a2，···， aq]，并计算滤波器的极点，若滤波器某极点位于单位圆外，则将该极点的倒数作为滤波器 的极点，并重新计算滤波器系数[I，ai，a2，…，aq]，同时令Xij = aj(j = 1，2，…，q)更新粒子的 前q维；
[0018] 步骤2·2·4:对当前粒子i，令bj = Xi(q+j+l)(j = 0，l，2，···，q-2)，bq-l = 0、bq = 0，构造 滤波器系数[bo，bi，b2，…，bq]，根据系数[I，ai，a2，…，aq]、[bo，bi，…，bq]由式(2)计算滤波 器的单位响应hi (η):
[0019]
(2)
[0020] 式中，δ(η)为单位函数，当η = 〇时，其值为1，否则其值为〇;
[0021 ]步骤2.2.5:按式计算滤波器最小噪声增益gmin:
[0022]
(3)
[0023] 步骤2.2.6:如果gmin小于等于g，转步骤2.2.7,否则按式(4)更新粒子i的后q-1维， 并转步骤2.2.4;
[0024]
(4).;：
[0025] 步骤2.2.7:对当前粒子i，由式(5)构建矩阵心、A2，按式(6)计算向量c:
[0026]
(5)
[0027] c = [l 0 0 ··· OKA1+A2)-1 (6)；
[0028] 步骤2.2.8:由式(7)~(9)构建矩阵HH3,按式（10)将向量c拆为向量 C1、c2,按 式(11)构建向量b1;
[0034] 步骤2.2.9:由式（12)计算参数8〇、(：111与(3112，由式（13)计算参数11、8 :
[0029]
[0030]
[0031]
[0032]
[0033]
[0035] (12)
[0036] (13)?
[0037] 步骤2.2.10:由式（14)计算一元二次方程系数A、B、C:
[0038]
(14),
[0039]步骤2.2.11:计算一元二次方程4#+8丨+0 = 0的根，如果两根为实根，分别将两根 作为滤波器系数bq-i，并按式(15)计算滤波器系数bq，转步骤2.2.12;否则转步骤2.2.15;
[0040]
(15)；
[00411 步骤2.2.12:根据步骤2.2.2所得滤波器系数ai，a2，…，a q、bo，bi，…，bq-2，结合步骤 2.2.11计算所得^^与^的两组值构建两个滤波器，按式（16)分别计算两滤波器的阶跃响 应，并以阶跃响应s (η)达到0.9时的η值，确定响应时间tm、tn2，选响应时间短的滤波器作为 设计结果，取得滤波器系数[I，ai，a2，…，aq]、[bo，bi，…，bq]与响应时间tn;
[0042]
(16)
[0043] 其中，单位阶跃函数定义如下：
[0044]
(17),
[0045] 步骤2.2.13:对于第i个粒子，将响应时间tn与其历史最优位置对应的响应时间Ti 进行比较，如果tn小于Ti，令Ti = tn、pi = Xi、ai = [I，ai，a2,…，aq]、bi = [bo，bi，…，bq];
[0046] 步骤2·2· 14:将Ti 与 Tg 进行比较，如果Ti 小于Tg，令 Tg = Ti、pg = pi、ag = ai、bg = bi;
[0047] 步骤2.2.15:令i = i+1，如果i〈M，那么转步骤2.2.2，否则转步骤2.3;
[0048] 步骤2.3:令k = k+l，如果k〈K，转步骤2.2.1，否则结束优化设计；
[0049] 步骤3:将全局最优位置对应的滤波器输出，作为优化设计结果。
[0050] 本发明的优势在于：
[0051 ] (1)可根据需要指定滤波器的噪声增益，控制测量误差；
[0052] (2)优化设计出的滤波器与传统平均或加权平均算法相比，在相同噪声增益情况 计算量很小；
[0053] (3)优化设计出的滤波器与现有滤波器设计工具设计的滤波器相比，在相同噪声 增益情况下具有更快的响应速度。
【附图说明】
[0054] 图1数字滤波器优化设计方法程序流程图
[0055] 图2传统方法设计滤波器与本方法设计滤波器阶跃响应比较图
【具体实施方式】
[0056] 设核辐射测量装置要求测量所得辐射强度（剂量率）的误差小于10%，根据30-法 贝1J，要求平均等待时间测均方根误差小于1/30,那么平均时所需的测量数据应不少于900， 此处选择1000,即平均处理后，噪声功率变为原来的1/1000。按照文献《单管宽量程剂量仪 关键技术分析》(第十六届全国核电子学与核探测技术学术年会论文集(上册），2012)修改 滤波器指标的方法，用传统滤波器设计方法，设计椭圆低通滤波器EF。设置滤波器EF的通带 归一化宽度为〇. 001，通带衰减为IdB，阻带归一化宽度为0.003，阻带衰减为40dB，设计出的 滤波器为3阶滤波器，该滤波器的噪声增益为0.00102,与平均方法噪声增益相当。采用1000 个测量数据进行平均、加权平均及采用滤波器EF处理三者的存储量、计算量对比见表1所 不。
[0057] 表1平均、加权平均、滤波器EF处理三者的存储量、计算量对比
[0059]从表1可以看出，滤波器EF的计算量与存储量均明显少于平均或加权平均方法，也 就是说，在平滑性能相近情况下，核探测中选用IIR滤波器进行平滑处理可降低测量系统存 储与计算要求。滤波器EF的系数见表2所示。
[0060] 本发明以上述设计低通滤波器EF为参照，设置噪声增益为0.00102(与滤波器EF相 同），滤波器阶数同样为3阶。由于两滤波器阶数相同，其计算量与存储量相同，二者噪声增 益相同，其测量精度也相同，对比二者响应速度便能分析二者的优劣。具体步骤如下：
[0061 ]步骤1:初始化参数
[0062] 设置滤波器噪声增益g = 0.00102、滤波器阶数q = 3、修正因子f = 0.9、粒子群规模 M= 100、最大迭代次数K= 10000、Co = I、C1 = 2、C2 = 2、粒子最大速度Vmax= 100，随机构造 M个 向量11={如412，"_，叫2(1-1)}分别代表1个粒子当前位置，设置迄今为止全局最优位置 1^ = (Pgi，Pg2，…，Pg(2q-1)}为任思值，对应的滤波器系数ag= {ago，agi，ag2，…，agq}、bg = {bg〇，bgi， bg2，…，bgq}为任意值，响应时间Tg为无穷大，各粒子的最优位置Pi = {pil，Pi2，···，P