一種用在控制棒棒位測量中的擬合方法*

張 菡，薄涵亮，楊文龍，王帥，傅一帆

(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是卡爾曼增益，為濾波器的中間結果?？柭鼮V波算法的關鍵就是求取這個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，分別對最小二乘法和經過卡爾曼濾波之后的最小二乘法進行數據擬合，根據擬合結果可知，經過卡爾曼之后的數據擬合的線性誤差小于直接使用最小二乘法進行數據擬合。綜合而言，在控制棒落棒時，對測量數據進行卡爾曼濾波之后再次擬合，擬合效果會更好。