卡尔曼滤波反演核设施核事故中核素释放率的研究

唐秀欢，包利红，李 华，万俊生

(西北核技术研究所，陕西 西安 710024)

估算核设施核事故源项一般有3种方法：监测手段、安全分析中的事故分析及事故中核设施工况的详细分析[1]。后两种方法往往过于保守，或在事故时难以实施，例如福岛核事故早期全厂断电，详细分析难以完成。因此，根据核设施周围的监测数据来反推源项的方法即源项反演得到重视。目前核事故源项反演主要分为以下4类[2]：最优插值法、遗传算法、人工神经网络、卡尔曼滤波及由其发展的扩展卡尔曼滤波和集合卡尔曼滤波。Jeong等[3]使用最优插值法估计了Yeoung-Kwang核设施示踪气体实验的点源释放率。陈晓秋等[4]采用线性化处理反推了福岛核事故的核素释放量。宁沙沙等[5]通过遗传编码搜索了最佳的核素释放率、释放高度及释放位置参数；陈媛华[6]建立了环境污染事件源项反演的遗传算法程序；陈军明等[7]使用遗传算法和多点源扩散模式，从控制点浓度成功地反演了点源排放清单。Astrup等[8-9]采用扩展卡尔曼方法跟踪了单烟团的扩散过程，反演了烟团位置和烟团活度参数。朱江等[10]、Zheng等[11]研究了集合卡尔曼方法对污染源的反演应用，讨论了污染源反演问题的统计方法。

核事故应急监测数据往往是核设施周边环境固定、移动的剂量或地面空气活度浓度，以监测数据反演源项，必然涉及测量和预测的误差；同时，由于监测数据具有实时和连续性，反演过程应具有连续数据递归处理的能力。卡尔曼滤波是一种最优递归资料处理算法，它综合利用一切可能的观测信息以及模式和观测信息的误差统计特征，使估计的统计误差达到最小。本文研究卡尔曼滤波反演核事故核素释放率算法，实现其在核设施核事故应急中的有效应用。

1 原理和方法

1.1 状态和测量方程

对于非移动的核设施，核应急重点关注的核素释放率状态方程如下：

Xk+1=Xk+qk

(1)

式中：X为状态向量；k为采样系列；qk为过程激励噪声。

测量方程为线性方程：

Zk=HkXk+rk

(2)

式中：Zk为源项经大气扩散后在核设施周围环境地面空气的瞬时活度浓度测量值；rk为观测噪声；Hk为观测矩阵。采用文献[12]描述的高斯多烟团模型计算放射性核素的大气扩散，在t时刻，地面空气的瞬时活度浓度与事故源项的释放率、释放时间等有关，假设在t=0时，第i烟团释放出来(时段释放率为Xi)，第i烟团在tw时刻即第w时段在某网格(x，y，0)上的地面空气活度浓度为：

(3)

(4)

(5)

扩散参数计算采用Briggs修正公式[12]，考虑干沉积过程造成的烟羽耗损和大气混合层对烟羽扩散的影响。

1.2 算法流程实现

将核事故核素释放率状态噪声和地面活度浓度的观测噪声均视为高斯白噪声，状态方程和观测方程的误差协方差矩阵分别为Q和R，则卡尔曼滤波算法步骤如下。

1) 计算核事故核素释放率预估计值X：

X(k|k-1)=A·X(k-1|k-1)

(6)

其中：A为状态转移矩阵，假定释放时释放率恒定，A=1；X(k|k-1)为利用上一状态预测的核素释放率结果；X(k-1|k-1)为上一状态最优的核素释放率结果。

2) 计算预估计协方差矩阵P：

P(k|k-1)=A·P(k-1|k-1)·A′+Q

(7)

其中，Q=E(qk，qk)为状态方程误差协方差矩阵，不同时刻的观测结果其平行性和准确性均不同，在迭代中其值亦不同。

3) 计算卡尔曼增益矩阵Kg：

(8)

其中，R=E(rk，rk)为观测方程误差协方差矩阵，当测量噪声极大时，残余中可用的信息量极低，增益很小，卡尔曼滤波器无法发挥作用，此时应使用先验估计值。在大气扩散计算中Hk是变量，每次测量数值的激励系数均不同。

4) 根据预测值和测量值计算核事故核素释放率更新估计值X：

X(k|k)=X(k|k-1)+Kg(k)·

(Z(k)-Hk·X(k|k-1))

(9)

5) 计算更新后估计值的协方差矩阵P：

P(k|k)=(I-Kg(k)·Hk)·P(k|k-1)

(10)

综合以上各式，可得卡尔曼滤波的计算流程如图1所示。

先验估计值采用源项经验判断结果，或用最优插值法计算第1批次数据作为先验估计值。

图1 卡尔曼滤波算法流程

2 数值实验与结果分析

2.1 数值实验

为验证卡尔曼滤波方法实时反演核事故核素释放率的可行性及其精度，本文设计了以下两种方案的反演数值实验。第1种，假定核事故核素释放率是稳定均匀的，即不考虑核事故核素释放率随时间的变化；第2种，考虑核事故核素释放率随采样点数而变化。实验过程设计如下：1) 假设释放高度为35 m，释放时大气稳定度为D级，10 m高度风速为4 m/s，风向为SE，释放时间为10 min；2) 核素监测点放置在核设施外半径200、300、400、500、600 m处，每个半径按π/4角度平均布置8个，共计40个，计算时假设只有下风向及其邻近的监测点响应核事故的核素释放，以核设施为原点，取包络下风向的π扇面监测点数据为输入，共计25个；3) 以真值加入白噪声后所得数据为监测数据，在事故后2 min时第1批25个监测数据传回，每隔2 min下一批数据传回，每个数据点作一次反演，由卡尔曼滤波器进行数据处理并反演事故释放率；4) 状态方程和观测方程的误差协方差矩阵Q和R用过程噪声相对误差10%和观测相对误差50%来计算。

2.2 释放率不随时间变化条件

在核设施发生核事故，核素释放率不随时间变化，其为一常数时，设核素释放率为105Bq/s。采用环境监测数据及不同初始X0值对核素释放率参数进行卡尔曼滤波反演，2～10 min内真值、测量值及卡尔曼滤波反演结果示于图2。

由图2可见，核设施核事故核素释放率不随时间变化时，卡尔曼滤波法可跟踪核素释放率，其反演结果较测量数据平滑。在卡尔曼滤波迭代初始值较真值低时，卡尔曼滤波法反演约10次后追踪到真值，随反演点数的加大，反演标准差逐渐比测量标准差小(图2a、b)；在卡尔曼滤波迭代初始值较真值高(真值的105倍)时，卡尔曼滤波法经60次反演后追踪到真值，由于初始误差较大，总体反演标准差比测量标准差高(图2c、d)；由反演效果看，较小的初始X0值具有较快的追踪速度和较小的标准差(图2a、b)。可见，当核设施核事故核素释放率不随时间变化时，根据环境监测数据进行的卡尔曼滤波法很快收敛，可有效跟踪核素释放率。

2.3 释放率随采样点数变化条件

核设施发生核事故时，核素释放率往往随时间发生变化，先假设核素释放率真值随反演点数线性变化，即X=105t(t=1，2，…，125)；再假设核素释放率真值随点数非线性变化，即X=|105-102(t-10)2+103t|(t=1，2，…，125)。采用环境监测数据及不同初始P0值对核素释放率参数进行卡尔曼滤波反演，真值、测量值及卡尔曼滤波反演结果示于图3。

由图3可见：在核素释放率真值线性变化条件下，根据环境监测数据进行卡尔曼滤波反演可跟踪真值的线性变化；同样在核素释放率真值非线性变化条件下，根据环境监测数据进行卡尔曼滤波反演也能追踪核素释放率的变化。在线性和非线性两种变化条件下，卡尔曼滤波均经过约10次反演后即可追踪到核素释放率真值。在误差方面，不同初始P0值对反演结果影响很小，经过约20次反演后卡尔曼滤波的标准差比测量值的标准差小，总体标准差随真值的增大而增大。可见，当核设施核事故核素释放率随时间变化时，根据环境监测数据进行的卡尔曼滤波法约10次追踪到真值，可有效跟踪核素释放率的线性或非线性变化。

2.4 讨论

核设施核事故放射性核素排放可能出现很复杂的情况：时间上不连续性，如短时连续释放、长时间间断释放；空间上不均匀性，如有效释放高度变化、气象变化造成的紊流扩散等。核事故源项的状态向量可能需要包含释放位置、释放高度、释放率等。本文设计的卡尔曼滤波法主要对线性系统进行分析，适用于固定高度和位置、短时连续排放的核事故核素释放率参数反演，且前提条件是假设环境监测传回的数据经过预先处理后形成核素释放率数据的测量值。常用源项反演最佳插值法是卡尔曼滤波法的简化版，主要用于线性系统的单变量分析，但对连续批量数据的处理能力不足。实际情况是释放高度、释放位置等多参数向量的观测方程和状态转移方程均具有很强的非线性。对于非线性特性，可采用扩展卡尔曼滤波进行基于泰勒级数展开的线性化处理，得到包括一阶近似项的状态方程和测量方程的近似表述形式。本文在不同释放高度位置对释放高度参数展开泰勒级数，分析一阶近似项截断误差，结果示于图4。由图4可见，在不同展开高度，当释放高度真值不同时，截断相对误差不同，真值距展开点越远，截断相对误差越大，且数量级可达1015，远超出了反演精度要求。

图2 释放率不随时间变化时卡尔曼滤波反演结果

图3 释放率随时间变化时卡尔曼滤波反演结果

图4 释放高度扩展卡尔曼滤波截断相对误差

大多核事故的位置地点是已知的，本文在状态向量选取中只考虑释放率，文献[8-9]的状态向量则更加关注移动的核事故释放源，考虑的源项参数有位置、释放率、高度等，例如单烟团反演系统状态向量(x，y，z，q)[8]。对于源项的其他参数如泄漏时间、泄漏总量等特征量尚未见反演的研究报道。另外，本文测量方程中测量值为核素地面空气的瞬时活度浓度，而非容易获得的地面空气的γ吸收剂量率。原因在于γ空气吸收剂量率随核素组成、浓度大小、烟云空间分布等变化而变化，如何通过γ空气吸收剂量率测量数据来确定每个核素的含量仍是复杂的技术问题。对于卡尔曼滤波法，未来的研究重点是协方差数据和γ空气吸收剂量率与测量方程的关系，特别是核事故应急测量的活度以及大气扩散模型的误差，关键核素的γ空气吸收剂量率与观测方程的物理关系。通过减少模型误差，简化输入参数，从而提高卡尔曼滤波法核素释放率的反演能力。

3 结论

根据核设施核事故环境监测数据的实时性和连续性特点，采用与高斯多烟团扩散模型结合的卡尔曼滤波法可用于固定高度和位置、短时连续排放的核事故核素释放率参数的反演，这为核设施核事故应急提供一种新的源项反演手段。由于释放高度参数的源项反演观测方程和状态转移方程具有较强的非线性，近似处理后截断误差不能满足反演要求，基于近似处理的扩展卡尔曼滤波法在释放高度参数源项反演的应用中有着一定的局限性。本文重点关注核事故核素释放率参数，在将来的工作中，卡尔曼滤波法观测方程和状态转移方程实际协方差数据的获取和分析，以及γ空气吸收剂量率与测量方程物理关系的研究将是决定该方法能否得到实际应用的关键。

参考文献：

[1] 王醒宇，康凌. 核事故后果评价方法及其新发展[M]. 北京：原子能出版社，2003：11.

[2] 徐志新，奚树人，曲静原. 核事故源项反演技术及其研究现状[J]. 科技导报，2007，25(5)：16-20.

XU Zhixin, XI Shuren, QU Jingyuan. Review on source inversion technology in analyzing nuclear accidents[J]. Science & Technology Review, 2007, 25(5): 16-20(in Chinese).

[3] JEONG H J, KIM E H, SUH K S, et al. Determination of the source rate released into the environment from a nuclear power plant[J]. Radiation Protection Dosimetry, 2005, 113(3): 308-313.

[4] 陈晓秋，杨端节，李冰. 利用福岛第一核电厂事故期间环境监测资料反推事故释放源项[J]. 核化学与放射化学，2012，34(2)：83-87.

CHEN Xiaoqiu, YANG Duanjie, LI Bing. Reverse estimation of accidental release amounts from Fukushima Daiichi Nuclear Power Plant by environment monitoring date[J]. Journal of Nuclear and Radiochemistry, 2012, 34(2): 83-87(in Chinese).

[5] 宁莎莎，蒯琳萍. 混合遗传算法在核事故源项反演中的应用[J]. 原子能科学技术，2012，46(增刊)：469-472.

NING Shasha, KUAI Linping. Back calculation of source terms by hybrid genetic algorithm in nuclear power plant accident[J]. Atomic Energy Science and Technology, 2012, 46(Suppl.): 469-472(in Chinese).

[6] 陈媛华. 河流突发环境污染事件源项反演及程序设计[D]. 哈尔滨：哈尔滨工业大学，2011.

[7] 陈军明，徐大海，朱蓉. 遗传算法在点源扩散浓度反演排放源强中的应用[J]. 气象，2002，28(9)：12-16.

CHEN Junming, XU Dahai, ZHU Rong. Application of genetic algorithms to point-source inversion[J]. Meteorological Monthly, 2002, 28(9): 12-16(in Chinese).

[8] ASTRUP P, TURCANU C, PUCH R O, et al. Data assimilation in the early phase: Kalman filtering RIMPUFF, RODOS(RA5)-TN(04)-01[R]. Denmark: Risø National Laboratory Information Service Department, 2004.

[9] PUCH R O, ASTRUP P. An extended Kalman filter methodology for the plume phase of a nuclear accident, RODOS(RA5)-TN(01)-07[R]. Denmark: Risø National Laboratory Information Service Department, 2002.

[10] 朱江，汪萍. 集合卡尔曼平滑和集合卡尔曼滤波在污染源反演中的应用[J]. 大气科学，2006，30(5)：871-882.

ZHU Jiang, WANG Ping. Ensemble Kalman smoother and ensemble Kalman filter approaches to the joint air quality state and emission estimation problem[J]. Chinese Journal of Atmospheric Sciences, 2006, 30(5): 871-882(in Chinese).

[11] ZHENG D Q, LEUNG J K C, LEE B Y, et al. Data assimilation in the atmospheric dispersion model for nuclear accident assessments[J]. Atmospheric Environment, 2007, 41: 2 438-2 446.

[12] 胡二邦，陈家宜. 核电厂大气扩散及其环境影响评价[M]. 北京：原子能出版社，1999.