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

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

(1.華北計(jì)算機(jī)系統(tǒng)工程研究所，北京 100083；2.清華大學(xué) 核能與新能源技術(shù)研究院，北京 100084)

0 引言

在核反應(yīng)堆中，控制棒(中子吸收體)吸收中子，通過(guò)控制棒插入或抽出核反應(yīng)堆的深度，控制核反應(yīng)的速度，當(dāng)遇到需要停堆的緊急情況時(shí)，控制棒以自由落體的方式快速落下，達(dá)到停止核反應(yīng)的安全狀態(tài)，因此對(duì)于控制棒棒位的測(cè)量是對(duì)反應(yīng)堆進(jìn)行控制和安全停堆的重要依據(jù)[1-2]。

電容傳感器的輸出/輸入關(guān)系按照曲線形狀可分為線性特性和非線性特性[3]。傳感器線性校正方法按照軟硬件構(gòu)成可分成兩大方向：一種是硬件電路補(bǔ)償方法，另一種軟件補(bǔ)償方法[4]。和軟件補(bǔ)償方法比較，硬件電路補(bǔ)償方法通常需要依據(jù)傳感器本身的電路設(shè)計(jì)來(lái)設(shè)計(jì)補(bǔ)償電路，由于設(shè)計(jì)電路復(fù)雜且難度大，實(shí)現(xiàn)較為困難，而軟件補(bǔ)償方法相對(duì)簡(jiǎn)單很多。

對(duì)于軟件補(bǔ)償法而言，使用頻率最高的方法有查表法、線性插值法、最小二乘法。其中，查表法[9]是根據(jù)精度要求對(duì)非線性曲線進(jìn)行分段，用多段折線逼近曲線，將折點(diǎn)坐標(biāo)值存入數(shù)據(jù)表中。但是查表法具有明顯缺點(diǎn)，當(dāng)測(cè)量范圍變大時(shí)，表的長(zhǎng)度隨之增加，將會(huì)占用大量的存儲(chǔ)空間，且造成查表速度變慢。線性插值法類似于模擬線性化方法中的非線性函數(shù)的折線近似逼近。但是線性插值大也需要占用大量的存儲(chǔ)空間，對(duì)傳感器性能具有較高的要求。相比較而言，最小二乘法簡(jiǎn)單快捷、易于實(shí)現(xiàn)、靈活性強(qiáng)且精度高。最小二乘法是在回歸分析中使用最多的方法之一，但是由于未考慮到誤差的影響，擬合效果并非特別理想。

由于測(cè)量數(shù)據(jù)會(huì)受到噪聲的干擾，測(cè)量值并非完全準(zhǔn)確，而卡爾曼濾波是一種根據(jù)測(cè)量值和經(jīng)驗(yàn)值去除噪音干擾，有效削弱誤差的影響[5]從而得到系統(tǒng)最優(yōu)估計(jì)真值的方法。因此提出了一種基于最小二乘的卡爾曼濾波方法，可以在控制棒棒位測(cè)量系統(tǒng)得到良好的應(yīng)用。

1 數(shù)據(jù)標(biāo)定系統(tǒng)

在進(jìn)行控制棒位移標(biāo)定時(shí)，控制棒套筒尾端端子(高低電平和地線)通過(guò)線纜與測(cè)量?jī)x接口連接，測(cè)量?jī)x通過(guò)Modbus RS-485接口轉(zhuǎn)換為Modbus RS-232接口線與個(gè)人計(jì)算機(jī)(控制棒棒位測(cè)量?jī)x上位機(jī)軟件)的USB口連接。建立控制棒棒位測(cè)量系統(tǒng)連接圖如圖1所示。

圖1 控制棒棒位測(cè)量系統(tǒng)圖

在進(jìn)行標(biāo)定實(shí)驗(yàn)時(shí)，總行程840 mm控制棒初始位置設(shè)置為0 mm行程，控制棒每前進(jìn)15 mm，控制棒測(cè)量當(dāng)前行程的電容值。由于電容值為非穩(wěn)定值，在每個(gè)行程，應(yīng)多測(cè)量幾次電容值取其平均值。本次實(shí)驗(yàn)測(cè)量次數(shù)設(shè)置為50次。將測(cè)量得到的電容值取平均值，得到行程與電容值對(duì)應(yīng)的數(shù)據(jù)如表1所示。

2 最小二乘法及其應(yīng)用

2.1 最小二乘法原理

最小二乘法曲線擬合是指：給定一組數(shù)據(jù)點(diǎn)(x1,y1),(x2,y2),…,(xn,yn)，由于無(wú)法求得唯一的函數(shù)多項(xiàng)式y(tǒng)=f(x)可以使所有數(shù)據(jù)點(diǎn)全部經(jīng)過(guò)，因此只能擬合一個(gè)近似曲線y=φ(x)使得對(duì)于所有的數(shù)據(jù)點(diǎn)來(lái)說(shuō)總的擬合誤差最小。即求：

(1)

表1 電容與行程關(guān)系

設(shè)擬合多項(xiàng)式為：

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

(2)

各個(gè)數(shù)據(jù)點(diǎn)到多項(xiàng)式曲線的距離之和的平方為偏差平方之和：

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

(3)

以上問(wèn)題等價(jià)于求式(4)中的li的值，對(duì)式(3)中l(wèi)j求偏導(dǎo)，其中，j=0,1,2,…,k，可以得到下面一組等式：

(4)

其中，j=0,1,2,…,k。將式(4)化簡(jiǎn)可得到：

(5)

根據(jù)式(5)可知X×L=Y，由于X和Y已知，可根據(jù)L=(X′X)-1X′Y求得L，從而得到擬合曲線[6]。

2.2 最小二乘法擬合分析

線性擬合又分為兩種，即單段線性擬合和分段線性擬合。單段線性擬合是整個(gè)設(shè)備采用一條線性公式，多段擬合是將行程平均分為多段，每段單獨(dú)進(jìn)行擬合，再合并在一起。線性擬合采用最小二乘法擬合曲線y=kx+b。當(dāng)進(jìn)行單段線性擬合時(shí)，分別求得擬合曲線的斜率k=0.013 2和截距b=71.857 6。

得到擬合曲線為：y=kx+b=0.013 2x+71.857 6。將行程x帶入擬合曲線中，可以求得擬合電容值y′，測(cè)量電容值y與擬合電容值y′的差得到殘差。利用MATLAB軟件繪制出擬合曲線圖和殘差圖，如圖2和圖3所示。

圖2 最小二乘法擬合曲線

圖3 最小二乘法線性擬合殘差圖

其中，圖2中data1表示實(shí)際測(cè)量電容值的散點(diǎn)值，data2表示最小二乘法擬合直線。從圖3可以看出，實(shí)際電容值與測(cè)量電容值的誤差大部分在0.25 pF 之內(nèi)，最大不超過(guò)0.317 6 pF。決定系數(shù)R2為：

(6)

最小二乘法誤差公式見(jiàn)式(7)：

(7)

在式(7)中，δL表示線性度，Δymax為最大電容值測(cè)量值與擬合直線對(duì)應(yīng)點(diǎn)的最大誤差絕對(duì)值，y為滿量程的電容測(cè)量值。根據(jù)以上測(cè)試數(shù)據(jù)可計(jì)算非線性誤差為：

(8)

3 卡爾曼濾波及其應(yīng)用

3.1 卡爾曼濾波原理

卡爾曼濾波是在已知系統(tǒng)狀態(tài)方程的情況下，使用系統(tǒng)的輸入和輸出測(cè)量數(shù)據(jù)來(lái)對(duì)系統(tǒng)狀態(tài)方程進(jìn)行最優(yōu)估計(jì)的方法[7-8]。

線性系統(tǒng)狀態(tài)預(yù)測(cè)方程如下所示：

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

(9)

線性系統(tǒng)測(cè)量方程：

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

(10)

其中，X(k)是k時(shí)刻系統(tǒng)狀態(tài)，U(k)是k時(shí)刻系統(tǒng)的控制量，Z(k)是k時(shí)刻傳感器測(cè)量值。A是狀態(tài)轉(zhuǎn)移系數(shù)，B是控制輸入的增益系數(shù)，H是指測(cè)量系數(shù)，在面對(duì)多模型系統(tǒng)時(shí)，A、B和H均為矩陣。W(k)是k時(shí)刻的過(guò)程激勵(lì)噪聲，Q是W(k)的協(xié)方差。

時(shí)間更新過(guò)程的公式如下所示[9]：

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

(11)

式(11)中，X(k-1|k-1)表示k-1時(shí)刻的最優(yōu)估計(jì)值，X(k|k-1)表示利用k-1時(shí)刻對(duì)k時(shí)刻的預(yù)測(cè)值，U(k)表示k時(shí)刻的控制量，當(dāng)位移測(cè)量系統(tǒng)沒(méi)有控制量時(shí)，它可以設(shè)置為0，大多數(shù)情況下沒(méi)有控制增益。

到現(xiàn)在為止，系統(tǒng)結(jié)果已更新，可是對(duì)應(yīng)于X(k|k-1)的協(xié)方差尚未更新。用P表示協(xié)方差：

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

(12)

式(12)中，A′表示A的轉(zhuǎn)置矩陣，Q是過(guò)程協(xié)方差，P(k|k-1)是X(k|k-1)的協(xié)方差，P(k-1|k-1)是與X(k-1|k-1)的協(xié)方差。已知預(yù)測(cè)值和測(cè)量值，便可以求出k時(shí)刻下的X(k|k)：

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

(13)

其中Kg是卡爾曼增益，為濾波器的中間結(jié)果?？柭鼮V波算法的關(guān)鍵就是求取這個(gè)Kg。此時(shí)，最小均方差就起到了作用。其中Kg計(jì)算方式如下：

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

(14)

到目前為止，已經(jīng)得到了k時(shí)刻下的最優(yōu)估計(jì)值X(k|k)。若想讓卡爾曼濾波器一直自回歸運(yùn)行，直到系統(tǒng)過(guò)程結(jié)束再停止，就還需要更新系統(tǒng)在k時(shí)刻下的最有估計(jì)值X(k|k)的協(xié)方差：

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

(15)

其中L為1的矩陣，對(duì)于單維數(shù)測(cè)量，L=1。當(dāng)系統(tǒng)進(jìn)入k+1狀態(tài)時(shí)，P(k|k)就是式(12)中的協(xié)方差P(k-1|k-1)。到現(xiàn)在為止，算法就可以一直運(yùn)行下去。

3.2 基于最小二乘法的卡爾曼濾波算法

當(dāng)對(duì)控制棒棒位測(cè)量系統(tǒng)進(jìn)行卡爾曼濾波處理時(shí)，需要先獲取系統(tǒng)的狀態(tài)方程和測(cè)量模型。將最小二乘法擬合得到的線性公式看作測(cè)量系統(tǒng)的先驗(yàn)估計(jì)，線性公式y(tǒng)=kx+b=0.013 2x+71.857 6。和大多數(shù)工業(yè)系統(tǒng)一樣，該系統(tǒng)沒(méi)有控制輸入，即U(k) = 0。首先建立控制棒系統(tǒng)的系統(tǒng)狀態(tài)方程：

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

(16)

X(k)=0.013 2k+71.857 6

(17)

則令

(18)

控制棒棒位測(cè)量?jī)x測(cè)量精度為0.000 1，選擇W(k)=1e-6。在實(shí)際現(xiàn)場(chǎng)應(yīng)用中，電容值時(shí)為步長(zhǎng)變化的量，根據(jù)上式得：

(19)

控制棒測(cè)量系統(tǒng)的測(cè)量方程為：

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

(20)

在數(shù)據(jù)標(biāo)定實(shí)驗(yàn)中，位移測(cè)量系統(tǒng)的測(cè)量值是狀態(tài)變量直接體現(xiàn)出來(lái)的，測(cè)量值中包括噪聲干擾，控制棒的測(cè)量值與理論值的誤差是由控制棒運(yùn)動(dòng)過(guò)程中的測(cè)量誤差引起的，也期望電容值測(cè)量值向理想值靠近，則：

(21)

由于在標(biāo)定測(cè)量時(shí)就是測(cè)量電容值，即Z(k)是已知的，且系統(tǒng)方程X(k)是已知的，則測(cè)量噪聲V(k)為：

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

(22)

在A、B、H、W(k)、V(k)已經(jīng)確定的情況下，把通過(guò)最小二乘法擬合得到的電容值視為經(jīng)驗(yàn)值，并且由測(cè)量?jī)x測(cè)量得到的電容值作為測(cè)量值，使用卡爾曼濾波算法計(jì)算電容最優(yōu)估計(jì)值，即得到一組新的電容值。濾波得到的最優(yōu)估計(jì)值與測(cè)量值之差最大為0.300 6。

由于得到最優(yōu)估計(jì)值依然是散點(diǎn)值，需要擬合出電容/位移公式，才能在控制棒落棒測(cè)量系統(tǒng)中應(yīng)用。則將經(jīng)過(guò)卡爾曼濾波的最優(yōu)真值使用最小二乘法進(jìn)行擬合。擬合圖像和散點(diǎn)圖如圖4和圖5所示。

圖4 帶卡爾曼濾波的擬合圖

圖5 帶卡爾曼濾波的最小二乘殘差圖

圖4中data1表示經(jīng)過(guò)卡爾曼濾波得到的散點(diǎn)值，data2表示實(shí)際測(cè)量電容值的散點(diǎn)值，data3表示最小二乘法擬合直線。從圖5中可以看出經(jīng)過(guò)卡爾曼濾波的線性最小二乘法殘差最大為0.300 6， 擬合結(jié)果優(yōu)于僅使用最小二乘法，新得到的最優(yōu)公式為y=kx+b=0.013 2x+71.830 8。根據(jù)以上測(cè)試數(shù)據(jù)可計(jì)算線性誤差為:

(23)

4 結(jié)論

本文使用MATLAB，分別對(duì)最小二乘法和經(jīng)過(guò)卡爾曼濾波之后的最小二乘法進(jìn)行數(shù)據(jù)擬合，根據(jù)擬合結(jié)果可知，經(jīng)過(guò)卡爾曼之后的數(shù)據(jù)擬合的線性誤差小于直接使用最小二乘法進(jìn)行數(shù)據(jù)擬合。綜合而言，在控制棒落棒時(shí)，對(duì)測(cè)量數(shù)據(jù)進(jìn)行卡爾曼濾波之后再次擬合，擬合效果會(huì)更好。