堆芯功率分布重构的两种空间统计性算法

彭星杰，李 庆，王 侃

(1.清华大学 工程物理系，北京 100084；2.中国核动力研究设计院 核反应堆系统设计技术重点实验室，四川 成都 610041)

为保证核电站运行的经济性和安全性，必须对与燃料元件和包壳屏障完整性有关的安全参数进行监测，以保证其不超出设计限值，其中两个最重要的参数是线功率密度(LPD)和偏离泡核沸腾比(DNBR)。第3代核电站能根据堆内探测器(功率分布在线监测系统)信号精确监测核电站相关安全参数。在线监测系统的核心技术之一是堆芯功率重构算法，即如何通过有限的堆内探测器读数重构出全堆的功率分布。国内外学者对这一问题进行了大量研究。Combustion Engineering公司开发的耦合系数法(CECOR)[1]计算流程简单且在工程上得到广泛应用，但其是建立在全堆分布变化不大的假设基础上。Webb[2]对耦合系数法进行了改进，发展出Lagrange乘子法，其计算思想和计算精度与CECOR的差别不大。Lee等[3]发展了最小二乘法，通过求解中子扩散方程与探测器响应方程形成的超定方程组来获得重构通量分布，但这种算法必须嵌入到堆芯扩散计算程序中，目前在商用堆芯核设计软件中还未得到应用。李富[4]开发了谐波综合法，利用中子扩散方程高阶谐波的有限阶线性组合来表示真实的堆芯通量分布。王常辉[5]通过建立谐波数据库实现了谐波综合法在线重构堆芯实时功率分布，但这种方法需预先计算好若干参考工况的高阶谐波，前期工作量巨大。Bryson等[6]将卡尔曼滤波方法应用于堆芯功率分布的重构，但未考虑观测量的归一化问题。在其他领域，如何通过有限的测量值得到物理场分布的估计也是一热点问题[7-9]。本工作在不考虑轴向功率分布重构的基础上，研究普通克里金方法与卡尔曼滤波方法在堆芯径向功率分布重构中的应用。

1 普通克里金方法

1.1 方法介绍

(1)

其中，权系数λi满足无偏性条件且使估值方差最小，即：

(2)

通过拉格朗日法求解，可得权系数满足如下方程组：

(3)

其中：μ为拉格朗日乘子；c(xi，xj)为观测点xi与xj之间的空间协方差函数；c(x，xj)为待估值点x与xj之间的空间协方差函数。

1.2 空间协方差函数

空间协方差函数代表了一种距离相关量，空间统计性理论认为估计值与观测值的相似程度是通过点对距离相关量来度量的。点对距离相关量只与观测场间的相互距离有关，而与它们的绝对位置无关。常用的距离相关量模型有指数模型、高斯模型等。根据具体情况可采用不同的相关函数，本工作采用二阶自回归模型[8](式(4))来描述距离相关量，其数学上的优越性已得到了充分理论证明[10]。

(4)

其中：L为特征距离；rij为两点间的距离。

1.3 普通克里金方法在功率重构中的应用

用探测器测量堆内功率，其读数可通过一定的转化关系化为对应组件的功率，中子扩散程序计算在同等堆芯状态下的功率分布。定义功率分布偏差函数为：

(5)

其中：n为探测器的数目；xi为放置了探测器的组件空间位置；Pmea(xi)为探测器测量的组件功率；Ppre(xi)为理论计算的组件功率。

将偏差函数按式(1)对未放置探测器的组件进行克里金插值，可得全堆组件的偏差函数，与原理论计算功率值相乘，得到未放置探测器组件的重构功率分布，再进行一次归一化操作后，便可得到全堆组件的重构功率分布。数值实验表明，当特征距离L取为2.5个组件长度后，使用普通克里金方法对秦山第二核电厂3号机组和大亚湾核电站1号机组功率进行重构时，重构功率分布与参考功率分布符合较好，即重构精度较高，故在程序中设定特征距离的默认值为2.5个组件长度。

2 卡尔曼滤波方法

2.1 方法介绍

卡尔曼滤波方法[8]是一种典型的数据同化方法。数据同化方法是一种广义上的最小二乘方法，求解原有理论计算值与实际测量值所构成的广义最小二乘问题。对于物理场Z，建立如下目标函数：

J(Z)=(Z-Zb)TB-1(Z-Zb)+

(H(Z)-Y0)TR-1(H(Z)-Y0)

(6)

其中：Zb为初始的物理场估计；Y0为物理场的相关观测值；H为相应的观测算符；B为物理场的协方差矩阵；R为观测量的协方差矩阵。

通过极小化上述目标函数，可得到上述问题的最小二乘解。

当观测算符H为线性化算符时，可得到目标函数对应的最小二乘解为：

Za=Zb+K(Y0-HZb)

(7)

其中，K为增益矩阵。

K=BHT(HBHT+R)-1

(8)

当观测算符H为非线性时，式(7)变为：

Za=Zb+K(Y0-H(Zb))

(9)

式(9)中的H变为初始估计物理场Zb处的雅可比矩阵。

2.2 堆芯功率重构中的应用

将卡尔曼滤波方法应用于堆芯功率重构时，理论计算得到的功率分布即式(9)中的Zb，而探测器的归一化读数即为式(9)中Y0。应用式(9)即可得到经过重构后的堆芯功率分布，之后还需进行一次归一化。

1) 物理场协方差矩阵B

关于B的选取，Bryson等将其表示为物理模型与计算误差、输入参数不确定性导致的误差与中子通量的波动误差之和[8]，这种表述较为繁琐。B也可表示为：

B=σ2c

(10)

其中：σ为标准差系数，表示初始估计物理场与真实物理场之间的差异，一般取值为5%～7%，本文取σ2为0.005；矩阵c的元见式(4)。

数值实验表明，当L取为1.5个组件长度，使用卡尔曼滤波方法对秦山第二核电厂3号机组和大亚湾核电站1号机组功率进行重构时，重构功率分布与参考功率分布符合较好，即重构精度较高，故在程序中设定特征距离的默认值为1.5个组件长度。

2) 观测量协方差矩阵R

可认为各观测值之间无明显的关联，从而R为一对角阵。通常观测量的误差用观测值乘以一误差因子得到，有：

Rii=(α(Y0)i)2i=1,…,n

(10)

其中，α为误差因子，表示测量误差绝对值与真实物理值的比值，本文假设该比值为固定值1%，即α取为0.01。

3) 观测算符H

观测算符H可由两部分作用的乘积表述：第1部分为选取矩阵，也就是1个元素为0或1的矩阵，表示的意义为探测器组件在全部组件中的分布；第2部分为归一化算符，也就是对探测器的读数进行归一化操作，这种做法是符合实际工程应用的，因为探测器读数的绝对值没有意义，关注的仅是其相对值。由于归一化算符为非线性算符，观测算符H明显也是一非线性算符，其雅可比矩阵可表示为归一化算符的雅可比矩阵与选取矩阵的乘积。Bryson的工作未注意到这一问题而仅使用了选取矩阵[8]，导致观测算符仅为线性算符，并不符合实际情况。

3 数值验证

以秦山第二核电厂3号机组第1循环和大亚湾核电站1号机组第13循环的实测堆芯功率分布作为参考功率分布，对上述两种方法的重构精度进行评估。探测器的布置示于图1。图1中，灰色组件代表放置堆内探测器的组件。实测堆芯功率分布由堆内探测器实测数据经一定的拓展方法得到，常用于功率分布重构算法的验证[5]。

a——秦山第二核电厂3号机组；b——大亚湾核电站1号机组

利用普通克里金方法和卡尔曼滤波方法分别对秦山第二核电厂3号机组第1循环燃耗为437 MW·d/tU和4 219 MW·d/tU时的堆芯功率分布进行重构，并与耦合系数法的重构结果进行对比，1/4堆芯的功率分布重构值与测量值的相对误差示于图2。

分别对大亚湾核电站1号机组第13循环燃耗为150 MW·d/tU和2 576 MW·d/tU时的堆芯功率分布进行重构，并与耦合系数法的重构结果进行对比，1/4堆芯的功率分布重构值与测量值的相对误差示于图3。

由图2、3可看出，普通克里金方法与卡尔曼滤波方法均能精确重构堆芯功率分布。由图2、3可得到如下结果。

1) 对放置探测器的组件，重构功率的相对误差基本上由重新归一化引入，这里不再进行分析；

2) 对未放置探测器的组件，普通克里金方法和卡尔曼滤波方法得到的功率重构精度在绝大部分组件内会高于耦合系数法，这是因为普通克里金方法和卡尔曼滤波方法均属于空间统计性理论下的最优插值算法类，均通过求解某一目标函数的最优值来确定插值结果，相比耦合系数法这种经验性的插值算法更有数学上的合理性；

a——437 MW·d/tU；b——4 219 MW·d/tU

3) 对未放置探测器的组件，普通克里金方法和卡尔曼滤波方法得到的功率重构精度在少部分组件内并不高于耦合系数法，这是因为重构相对误差存在随机性，有正有负，而为了保证整体重构功率的归一化，必然会导致某些组件的重构精度不如耦合系数法。

从工程的角度考虑，重构得到的功率分布与测量值之间的最大相对误差(相对误差绝对值的最大值)不应大于某一限值，结果列于表1。

表1 重构值与测量值相对误差的限值

表2、3列出新提出的两种重构方法与耦合系数法的最大重构相对误差以及平均重构相对误差对比。由表2、3得到如下结果。

表2 秦山第二核电厂3号机组重构方法的相对误差对比

表3 大亚湾核电站1号机组重构方法的相对误差对比

1) 本工作提出的两种新方法明显满足工程上对监测堆芯功率的要求。

2) 对秦山第二核电厂3号机组的验证结果而言，两种新方法得到的组件功率最大重构相对误差及组件功率平均重构相对误差均小于耦合系数法的相应值。3种方法中重构精度最高的是普通克里金方法。

3) 对大亚湾核电站1号机组的验证结果而言，普通克里金方法得到的组件功率最大重构相对误差及组件功率平均重构相对误差均小于耦合系数法；卡尔曼滤波方法得到的组件功率平均重构相对误差小于耦合系数法，而最大重构相对误差大于耦合系数法，但最大重构相对误差出现在外围组件，影响不大。3种方法中重构精度最高的是普通克里金方法。

普通克里金方法的最佳权重系数对于1座反应堆只需计算1次，然后按式(1)进行求和，明显满足在线监测要求。对于卡尔曼滤波方法，在计算机硬件条件为CPU Intel(R) Core(TM) 3.00 GHz，内存3.47 GB的条件下，1个轴向层的径向功率拓展只需0.012 s，三维功率重构预计小于1 s，也明显满足在线监测要求。

4 结论与展望

本工作提出两种基于空间统计性理论的功率重构方法，并利用秦山第二核电厂3号机组和大亚湾核电站1号机组的实测功率分布进行方法验证，由验证结果得到以下结论。

1) 两种基于空间统计性理论的方法均能精确地进行堆芯功率重构，重构精度满足工程的要求；

2) 两种方法均有很快的计算速度，满足在线监测的实时性要求；

3) 在实际的工程应用中，两种方法均可单独使有，并可取代目前常用的耦合系数法进行功率分布重构。

在实际的反应堆在线监测应用中，可直接获取三维的堆内探测器测量数据，且本工作提出的两种方法均可拓展到三维，下一步工作将是两种方法直接在三维功率分布重构中的应用，具体可采取下面两种思路：

1) 通过三次样条函数拟合[11]进行轴向功率分布重构，每一轴向层的径向功率重构则采用本工作提出的两种方法；

2) 将普通克里金方法和卡尔曼滤波方法直接用于三维功率分布重构，其计算流程与二维径向功率分布重构一致。

参考文献：

[1]KARLSON C F. Continuing advancements in in-core power distribution measurement methods using SIMULATE-3 and CECOR3.4[J]. Nuclear Science and Engineering, 1995, 121(1): 57-66.

[2]WEBB R J. Comparison of CECOR algorithm to Lagrange multiplier method to estimate reactor power distributions[J]. Nuclear Technology, 2000, 132(2): 206-213.

[3]LEE K, KIM C H. The least-squares method for three-dimensional core power distribution monitoring in pressurized water reactors[J]. Nuclear Science and Engineering, 2003, 143(3): 268-280.

[4]李富. 重构堆芯通量分布的谐波综合法及其诊断应用[D]. 北京：清华大学，1994.

[5]王常辉. 谐波展开法及其在反应堆堆芯功率在线监测中的应用研究[D]. 西安：西安交通大学，2011.

[6]BRYSON J W, LEE J C, HASSBERGER J A. Optimal flux map generation through parameter estimation techniques[J]. Nuclear Science and Engineering, 1993, 114(3): 238-251.

[7]王勇标. 克里金算法的改进[D]. 荆州：长江大学，2012.

[8]翟进乾. 克里金插值方法在煤层分布检测中的应用研究[D]. 太原：太原理工大学，2008.

[9]马寨璞. 海洋流场数据同化方法与应用的研究[D]. 杭州：浙江大学，2005.

[10] GASPARI G, COHN S E. Construction of correlation functions in two and three dimensions[J]. Quarterly Journal of the Royal Meteorological Society, 1999, 125(5): 723-757.

[11] HAN S, KIM U, SEONG P. A methodology for benefit assessment of using in-core neutron detector signals in core protection calculator system for Korea standard nuclear power plants[J]. Annals of Nuclear Energy, 1999, 26(6): 471-488.