一种用在控制棒棒位测量中的拟合方法*

张 菡，薄涵亮，杨文龙，王帅，傅一帆

(1.华北计算机系统工程研究所，北京 100083；2.清华大学 核能与新能源技术研究院，北京 100084)

0 引言

在核反应堆中，控制棒(中子吸收体)吸收中子，通过控制棒插入或抽出核反应堆的深度，控制核反应的速度，当遇到需要停堆的紧急情况时，控制棒以自由落体的方式快速落下，达到停止核反应的安全状态，因此对于控制棒棒位的测量是对反应堆进行控制和安全停堆的重要依据[1-2]。

电容传感器的输出/输入关系按照曲线形状可分为线性特性和非线性特性[3]。传感器线性校正方法按照软硬件构成可分成两大方向：一种是硬件电路补偿方法，另一种软件补偿方法[4]。和软件补偿方法比较，硬件电路补偿方法通常需要依据传感器本身的电路设计来设计补偿电路，由于设计电路复杂且难度大，实现较为困难，而软件补偿方法相对简单很多。

对于软件补偿法而言，使用频率最高的方法有查表法、线性插值法、最小二乘法。其中，查表法[9]是根据精度要求对非线性曲线进行分段，用多段折线逼近曲线，将折点坐标值存入数据表中。但是查表法具有明显缺点，当测量范围变大时，表的长度随之增加，将会占用大量的存储空间，且造成查表速度变慢。线性插值法类似于模拟线性化方法中的非线性函数的折线近似逼近。但是线性插值大也需要占用大量的存储空间，对传感器性能具有较高的要求。相比较而言，最小二乘法简单快捷、易于实现、灵活性强且精度高。最小二乘法是在回归分析中使用最多的方法之一，但是由于未考虑到误差的影响，拟合效果并非特别理想。

由于测量数据会受到噪声的干扰，测量值并非完全准确，而卡尔曼滤波是一种根据测量值和经验值去除噪音干扰，有效削弱误差的影响[5]从而得到系统最优估计真值的方法。因此提出了一种基于最小二乘的卡尔曼滤波方法，可以在控制棒棒位测量系统得到良好的应用。

1 数据标定系统

在进行控制棒位移标定时，控制棒套筒尾端端子(高低电平和地线)通过线缆与测量仪接口连接，测量仪通过Modbus RS-485接口转换为Modbus RS-232接口线与个人计算机(控制棒棒位测量仪上位机软件)的USB口连接。建立控制棒棒位测量系统连接图如图1所示。

图1 控制棒棒位测量系统图

在进行标定实验时，总行程840 mm控制棒初始位置设置为0 mm行程，控制棒每前进15 mm，控制棒测量当前行程的电容值。由于电容值为非稳定值，在每个行程，应多测量几次电容值取其平均值。本次实验测量次数设置为50次。将测量得到的电容值取平均值，得到行程与电容值对应的数据如表1所示。

2 最小二乘法及其应用

2.1 最小二乘法原理

最小二乘法曲线拟合是指：给定一组数据点(x1,y1),(x2,y2),…,(xn,yn)，由于无法求得唯一的函数多项式y=f(x)可以使所有数据点全部经过，因此只能拟合一个近似曲线y=φ(x)使得对于所有的数据点来说总的拟合误差最小。即求：

(1)

表1 电容与行程关系

设拟合多项式为：

y=φ(x)=l0+l1x+…+lkxk

(2)

各个数据点到多项式曲线的距离之和的平方为偏差平方之和：

σ2(l0,l1,…,lk)=

(3)

以上问题等价于求式(4)中的li的值，对式(3)中lj求偏导，其中，j=0,1,2,…,k，可以得到下面一组等式：

(4)

其中，j=0,1,2,…,k。将式(4)化简可得到：

(5)

根据式(5)可知X×L=Y，由于X和Y已知，可根据L=(X′X)-1X′Y求得L，从而得到拟合曲线[6]。

2.2 最小二乘法拟合分析

线性拟合又分为两种，即单段线性拟合和分段线性拟合。单段线性拟合是整个设备采用一条线性公式，多段拟合是将行程平均分为多段，每段单独进行拟合，再合并在一起。线性拟合采用最小二乘法拟合曲线y=kx+b。当进行单段线性拟合时，分别求得拟合曲线的斜率k=0.013 2和截距b=71.857 6。

得到拟合曲线为：y=kx+b=0.013 2x+71.857 6。将行程x带入拟合曲线中，可以求得拟合电容值y′，测量电容值y与拟合电容值y′的差得到残差。利用MATLAB软件绘制出拟合曲线图和残差图，如图2和图3所示。

图2 最小二乘法拟合曲线

图3 最小二乘法线性拟合残差图

其中，图2中data1表示实际测量电容值的散点值，data2表示最小二乘法拟合直线。从图3可以看出，实际电容值与测量电容值的误差大部分在0.25 pF 之内，最大不超过0.317 6 pF。决定系数R2为：

(6)

最小二乘法误差公式见式(7)：

(7)

在式(7)中，δL表示线性度，Δymax为最大电容值测量值与拟合直线对应点的最大误差绝对值，y为满量程的电容测量值。根据以上测试数据可计算非线性误差为：

(8)

3 卡尔曼滤波及其应用

3.1 卡尔曼滤波原理

卡尔曼滤波是在已知系统状态方程的情况下，使用系统的输入和输出测量数据来对系统状态方程进行最优估计的方法[7-8]。

线性系统状态预测方程如下所示：

X(k)=AX(k-1)+BU(k)+W(k)

(9)

线性系统测量方程：

Z(k)=HX(k)+V(k)

(10)

其中，X(k)是k时刻系统状态，U(k)是k时刻系统的控制量，Z(k)是k时刻传感器测量值。A是状态转移系数，B是控制输入的增益系数，H是指测量系数，在面对多模型系统时，A、B和H均为矩阵。W(k)是k时刻的过程激励噪声，Q是W(k)的协方差。

时间更新过程的公式如下所示[9]：

X(k|k-1)=AX(k-1|k-1)+BU(k)

(11)

式(11)中，X(k-1|k-1)表示k-1时刻的最优估计值，X(k|k-1)表示利用k-1时刻对k时刻的预测值，U(k)表示k时刻的控制量，当位移测量系统没有控制量时，它可以设置为0，大多数情况下没有控制增益。

到现在为止，系统结果已更新，可是对应于X(k|k-1)的协方差尚未更新。用P表示协方差：

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

(12)

式(12)中，A′表示A的转置矩阵，Q是过程协方差，P(k|k-1)是X(k|k-1)的协方差，P(k-1|k-1)是与X(k-1|k-1)的协方差。已知预测值和测量值，便可以求出k时刻下的X(k|k)：

X(k|k)=X(k|k-1)+
Kg(k)(Z(k)-HX(k|k-1))

(13)

其中Kg是卡尔曼增益，为滤波器的中间结果。卡尔曼滤波算法的关键就是求取这个Kg。此时，最小均方差就起到了作用。其中Kg计算方式如下：

Kg(k)=P(k|k-1)H′/(HP(k|k-1)H′+R)

(14)

到目前为止，已经得到了k时刻下的最优估计值X(k|k)。若想让卡尔曼滤波器一直自回归运行，直到系统过程结束再停止，就还需要更新系统在k时刻下的最有估计值X(k|k)的协方差：

P(k|k)=(L-Kg(k)H)P(k|k-1)

(15)

其中L为1的矩阵，对于单维数测量，L=1。当系统进入k+1状态时，P(k|k)就是式(12)中的协方差P(k-1|k-1)。到现在为止，算法就可以一直运行下去。

3.2 基于最小二乘法的卡尔曼滤波算法

当对控制棒棒位测量系统进行卡尔曼滤波处理时，需要先获取系统的状态方程和测量模型。将最小二乘法拟合得到的线性公式看作测量系统的先验估计，线性公式y=kx+b=0.013 2x+71.857 6。和大多数工业系统一样，该系统没有控制输入，即U(k) = 0。首先建立控制棒系统的系统状态方程：

X(k-1)=0.013 2(k-1)+71.857 6

(16)

X(k)=0.013 2k+71.857 6

(17)

则令

(18)

控制棒棒位测量仪测量精度为0.000 1，选择W(k)=1e-6。在实际现场应用中，电容值时为步长变化的量，根据上式得：

(19)

控制棒测量系统的测量方程为：

Z(k)=HX(k)+V(k)

(20)

在数据标定实验中，位移测量系统的测量值是状态变量直接体现出来的，测量值中包括噪声干扰，控制棒的测量值与理论值的误差是由控制棒运动过程中的测量误差引起的，也期望电容值测量值向理想值靠近，则：

(21)

由于在标定测量时就是测量电容值，即Z(k)是已知的，且系统方程X(k)是已知的，则测量噪声V(k)为：

V(k)=Z(k)-X(k)

(22)

在A、B、H、W(k)、V(k)已经确定的情况下，把通过最小二乘法拟合得到的电容值视为经验值，并且由测量仪测量得到的电容值作为测量值，使用卡尔曼滤波算法计算电容最优估计值，即得到一组新的电容值。滤波得到的最优估计值与测量值之差最大为0.300 6。

由于得到最优估计值依然是散点值，需要拟合出电容/位移公式，才能在控制棒落棒测量系统中应用。则将经过卡尔曼滤波的最优真值使用最小二乘法进行拟合。拟合图像和散点图如图4和图5所示。

图4 带卡尔曼滤波的拟合图

图5 带卡尔曼滤波的最小二乘残差图

图4中data1表示经过卡尔曼滤波得到的散点值，data2表示实际测量电容值的散点值，data3表示最小二乘法拟合直线。从图5中可以看出经过卡尔曼滤波的线性最小二乘法残差最大为0.300 6， 拟合结果优于仅使用最小二乘法，新得到的最优公式为y=kx+b=0.013 2x+71.830 8。根据以上测试数据可计算线性误差为:

(23)

4 结论

本文使用MATLAB，分别对最小二乘法和经过卡尔曼滤波之后的最小二乘法进行数据拟合，根据拟合结果可知，经过卡尔曼之后的数据拟合的线性误差小于直接使用最小二乘法进行数据拟合。综合而言，在控制棒落棒时，对测量数据进行卡尔曼滤波之后再次拟合，拟合效果会更好。