彭星杰，李 慶，王 侃

(1.清華大學(xué) 工程物理系，北京 100084；2.中國核動力研究設(shè)計院 核反應(yīng)堆系統(tǒng)設(shè)計技術(shù)重點(diǎn)實(shí)驗(yàn)室，四川 成都 610041)

為保證核電站運(yùn)行的經(jīng)濟(jì)性和安全性，必須對與燃料元件和包殼屏障完整性有關(guān)的安全參數(shù)進(jìn)行監(jiān)測，以保證其不超出設(shè)計限值，其中兩個最重要的參數(shù)是線功率密度(LPD)和偏離泡核沸騰比(DNBR)。第3代核電站能根據(jù)堆內(nèi)探測器(功率分布在線監(jiān)測系統(tǒng))信號精確監(jiān)測核電站相關(guān)安全參數(shù)。在線監(jiān)測系統(tǒng)的核心技術(shù)之一是堆芯功率重構(gòu)算法，即如何通過有限的堆內(nèi)探測器讀數(shù)重構(gòu)出全堆的功率分布。國內(nèi)外學(xué)者對這一問題進(jìn)行了大量研究。Combustion Engineering公司開發(fā)的耦合系數(shù)法(CECOR)[1]計算流程簡單且在工程上得到廣泛應(yīng)用，但其是建立在全堆分布變化不大的假設(shè)基礎(chǔ)上。Webb[2]對耦合系數(shù)法進(jìn)行了改進(jìn)，發(fā)展出Lagrange乘子法，其計算思想和計算精度與CECOR的差別不大。Lee等[3]發(fā)展了最小二乘法，通過求解中子擴(kuò)散方程與探測器響應(yīng)方程形成的超定方程組來獲得重構(gòu)通量分布，但這種算法必須嵌入到堆芯擴(kuò)散計算程序中，目前在商用堆芯核設(shè)計軟件中還未得到應(yīng)用。李富[4]開發(fā)了諧波綜合法，利用中子擴(kuò)散方程高階諧波的有限階線性組合來表示真實(shí)的堆芯通量分布。王常輝[5]通過建立諧波數(shù)據(jù)庫實(shí)現(xiàn)了諧波綜合法在線重構(gòu)堆芯實(shí)時功率分布，但這種方法需預(yù)先計算好若干參考工況的高階諧波，前期工作量巨大。Bryson等[6]將卡爾曼濾波方法應(yīng)用于堆芯功率分布的重構(gòu)，但未考慮觀測量的歸一化問題。在其他領(lǐng)域，如何通過有限的測量值得到物理場分布的估計也是一熱點(diǎn)問題[7-9]。本工作在不考慮軸向功率分布重構(gòu)的基礎(chǔ)上，研究普通克里金方法與卡爾曼濾波方法在堆芯徑向功率分布重構(gòu)中的應(yīng)用。

1 普通克里金方法

1.1 方法介紹

(1)

其中，權(quán)系數(shù)λi滿足無偏性條件且使估值方差最小，即：

(2)

通過拉格朗日法求解，可得權(quán)系數(shù)滿足如下方程組：

(3)

其中：μ為拉格朗日乘子；c(xi，xj)為觀測點(diǎn)xi與xj之間的空間協(xié)方差函數(shù)；c(x，xj)為待估值點(diǎn)x與xj之間的空間協(xié)方差函數(shù)。

1.2 空間協(xié)方差函數(shù)

空間協(xié)方差函數(shù)代表了一種距離相關(guān)量，空間統(tǒng)計性理論認(rèn)為估計值與觀測值的相似程度是通過點(diǎn)對距離相關(guān)量來度量的。點(diǎn)對距離相關(guān)量只與觀測場間的相互距離有關(guān)，而與它們的絕對位置無關(guān)。常用的距離相關(guān)量模型有指數(shù)模型、高斯模型等。根據(jù)具體情況可采用不同的相關(guān)函數(shù)，本工作采用二階自回歸模型[8](式(4))來描述距離相關(guān)量，其數(shù)學(xué)上的優(yōu)越性已得到了充分理論證明[10]。

(4)

其中：L為特征距離；rij為兩點(diǎn)間的距離。

1.3 普通克里金方法在功率重構(gòu)中的應(yīng)用

用探測器測量堆內(nèi)功率，其讀數(shù)可通過一定的轉(zhuǎn)化關(guān)系化為對應(yīng)組件的功率，中子擴(kuò)散程序計算在同等堆芯狀態(tài)下的功率分布。定義功率分布偏差函數(shù)為：

(5)

其中：n為探測器的數(shù)目；xi為放置了探測器的組件空間位置；Pmea(xi)為探測器測量的組件功率；Ppre(xi)為理論計算的組件功率。

將偏差函數(shù)按式(1)對未放置探測器的組件進(jìn)行克里金插值，可得全堆組件的偏差函數(shù)，與原理論計算功率值相乘，得到未放置探測器組件的重構(gòu)功率分布，再進(jìn)行一次歸一化操作后，便可得到全堆組件的重構(gòu)功率分布。數(shù)值實(shí)驗(yàn)表明，當(dāng)特征距離L取為2.5個組件長度后，使用普通克里金方法對秦山第二核電廠3號機(jī)組和大亞灣核電站1號機(jī)組功率進(jìn)行重構(gòu)時，重構(gòu)功率分布與參考功率分布符合較好，即重構(gòu)精度較高，故在程序中設(shè)定特征距離的默認(rèn)值為2.5個組件長度。

2 卡爾曼濾波方法

2.1 方法介紹

卡爾曼濾波方法[8]是一種典型的數(shù)據(jù)同化方法。數(shù)據(jù)同化方法是一種廣義上的最小二乘方法，求解原有理論計算值與實(shí)際測量值所構(gòu)成的廣義最小二乘問題。對于物理場Z，建立如下目標(biāo)函數(shù)：

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

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

(6)

其中：Zb為初始的物理場估計；Y0為物理場的相關(guān)觀測值；H為相應(yīng)的觀測算符；B為物理場的協(xié)方差矩陣；R為觀測量的協(xié)方差矩陣。

通過極小化上述目標(biāo)函數(shù)，可得到上述問題的最小二乘解。

當(dāng)觀測算符H為線性化算符時，可得到目標(biāo)函數(shù)對應(yīng)的最小二乘解為：

Za=Zb+K(Y0-HZb)

(7)

其中，K為增益矩陣。

K=BHT(HBHT+R)-1

(8)

當(dāng)觀測算符H為非線性時，式(7)變?yōu)椋?/p>

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

(9)

式(9)中的H變?yōu)槌跏脊烙嬑锢韴鯶b處的雅可比矩陣。

2.2 堆芯功率重構(gòu)中的應(yīng)用

將卡爾曼濾波方法應(yīng)用于堆芯功率重構(gòu)時，理論計算得到的功率分布即式(9)中的Zb，而探測器的歸一化讀數(shù)即為式(9)中Y0。應(yīng)用式(9)即可得到經(jīng)過重構(gòu)后的堆芯功率分布，之后還需進(jìn)行一次歸一化。

1) 物理場協(xié)方差矩陣B

關(guān)于B的選取，Bryson等將其表示為物理模型與計算誤差、輸入?yún)?shù)不確定性導(dǎo)致的誤差與中子通量的波動誤差之和[8]，這種表述較為繁瑣。B也可表示為：

B=σ2c

(10)

其中：σ為標(biāo)準(zhǔn)差系數(shù)，表示初始估計物理場與真實(shí)物理場之間的差異，一般取值為5%～7%，本文取σ2為0.005；矩陣c的元見式(4)。

數(shù)值實(shí)驗(yàn)表明，當(dāng)L取為1.5個組件長度，使用卡爾曼濾波方法對秦山第二核電廠3號機(jī)組和大亞灣核電站1號機(jī)組功率進(jìn)行重構(gòu)時，重構(gòu)功率分布與參考功率分布符合較好，即重構(gòu)精度較高，故在程序中設(shè)定特征距離的默認(rèn)值為1.5個組件長度。

2) 觀測量協(xié)方差矩陣R

可認(rèn)為各觀測值之間無明顯的關(guān)聯(lián)，從而R為一對角陣。通常觀測量的誤差用觀測值乘以一誤差因子得到，有：

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

(10)

其中，α為誤差因子，表示測量誤差絕對值與真實(shí)物理值的比值，本文假設(shè)該比值為固定值1%，即α取為0.01。

3) 觀測算符H

觀測算符H可由兩部分作用的乘積表述：第1部分為選取矩陣，也就是1個元素為0或1的矩陣，表示的意義為探測器組件在全部組件中的分布；第2部分為歸一化算符，也就是對探測器的讀數(shù)進(jìn)行歸一化操作，這種做法是符合實(shí)際工程應(yīng)用的，因?yàn)樘綔y器讀數(shù)的絕對值沒有意義，關(guān)注的僅是其相對值。由于歸一化算符為非線性算符，觀測算符H明顯也是一非線性算符，其雅可比矩陣可表示為歸一化算符的雅可比矩陣與選取矩陣的乘積。Bryson的工作未注意到這一問題而僅使用了選取矩陣[8]，導(dǎo)致觀測算符僅為線性算符，并不符合實(shí)際情況。

3 數(shù)值驗(yàn)證

以秦山第二核電廠3號機(jī)組第1循環(huán)和大亞灣核電站1號機(jī)組第13循環(huán)的實(shí)測堆芯功率分布作為參考功率分布，對上述兩種方法的重構(gòu)精度進(jìn)行評估。探測器的布置示于圖1。圖1中，灰色組件代表放置堆內(nèi)探測器的組件。實(shí)測堆芯功率分布由堆內(nèi)探測器實(shí)測數(shù)據(jù)經(jīng)一定的拓展方法得到，常用于功率分布重構(gòu)算法的驗(yàn)證[5]。

a——秦山第二核電廠3號機(jī)組；b——大亞灣核電站1號機(jī)組

利用普通克里金方法和卡爾曼濾波方法分別對秦山第二核電廠3號機(jī)組第1循環(huán)燃耗為437 MW·d/tU和4 219 MW·d/tU時的堆芯功率分布進(jìn)行重構(gòu)，并與耦合系數(shù)法的重構(gòu)結(jié)果進(jìn)行對比，1/4堆芯的功率分布重構(gòu)值與測量值的相對誤差示于圖2。

分別對大亞灣核電站1號機(jī)組第13循環(huán)燃耗為150 MW·d/tU和2 576 MW·d/tU時的堆芯功率分布進(jìn)行重構(gòu)，并與耦合系數(shù)法的重構(gòu)結(jié)果進(jìn)行對比，1/4堆芯的功率分布重構(gòu)值與測量值的相對誤差示于圖3。

由圖2、3可看出，普通克里金方法與卡爾曼濾波方法均能精確重構(gòu)堆芯功率分布。由圖2、3可得到如下結(jié)果。

1) 對放置探測器的組件，重構(gòu)功率的相對誤差基本上由重新歸一化引入，這里不再進(jìn)行分析；

2) 對未放置探測器的組件，普通克里金方法和卡爾曼濾波方法得到的功率重構(gòu)精度在絕大部分組件內(nèi)會高于耦合系數(shù)法，這是因?yàn)槠胀死锝鸱椒ê涂柭鼮V波方法均屬于空間統(tǒng)計性理論下的最優(yōu)插值算法類，均通過求解某一目標(biāo)函數(shù)的最優(yōu)值來確定插值結(jié)果，相比耦合系數(shù)法這種經(jīng)驗(yàn)性的插值算法更有數(shù)學(xué)上的合理性；

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

3) 對未放置探測器的組件，普通克里金方法和卡爾曼濾波方法得到的功率重構(gòu)精度在少部分組件內(nèi)并不高于耦合系數(shù)法，這是因?yàn)橹貥?gòu)相對誤差存在隨機(jī)性，有正有負(fù)，而為了保證整體重構(gòu)功率的歸一化，必然會導(dǎo)致某些組件的重構(gòu)精度不如耦合系數(shù)法。

從工程的角度考慮，重構(gòu)得到的功率分布與測量值之間的最大相對誤差(相對誤差絕對值的最大值)不應(yīng)大于某一限值，結(jié)果列于表1。

表1 重構(gòu)值與測量值相對誤差的限值

表2、3列出新提出的兩種重構(gòu)方法與耦合系數(shù)法的最大重構(gòu)相對誤差以及平均重構(gòu)相對誤差對比。由表2、3得到如下結(jié)果。

表2 秦山第二核電廠3號機(jī)組重構(gòu)方法的相對誤差對比

表3 大亞灣核電站1號機(jī)組重構(gòu)方法的相對誤差對比

1) 本工作提出的兩種新方法明顯滿足工程上對監(jiān)測堆芯功率的要求。

2) 對秦山第二核電廠3號機(jī)組的驗(yàn)證結(jié)果而言，兩種新方法得到的組件功率最大重構(gòu)相對誤差及組件功率平均重構(gòu)相對誤差均小于耦合系數(shù)法的相應(yīng)值。3種方法中重構(gòu)精度最高的是普通克里金方法。

3) 對大亞灣核電站1號機(jī)組的驗(yàn)證結(jié)果而言，普通克里金方法得到的組件功率最大重構(gòu)相對誤差及組件功率平均重構(gòu)相對誤差均小于耦合系數(shù)法；卡爾曼濾波方法得到的組件功率平均重構(gòu)相對誤差小于耦合系數(shù)法，而最大重構(gòu)相對誤差大于耦合系數(shù)法，但最大重構(gòu)相對誤差出現(xiàn)在外圍組件，影響不大。3種方法中重構(gòu)精度最高的是普通克里金方法。

普通克里金方法的最佳權(quán)重系數(shù)對于1座反應(yīng)堆只需計算1次，然后按式(1)進(jìn)行求和，明顯滿足在線監(jiān)測要求。對于卡爾曼濾波方法，在計算機(jī)硬件條件為CPU Intel(R) Core(TM) 3.00 GHz，內(nèi)存3.47 GB的條件下，1個軸向?qū)拥膹较蚬β释卣怪恍?.012 s，三維功率重構(gòu)預(yù)計小于1 s，也明顯滿足在線監(jiān)測要求。

4 結(jié)論與展望

本工作提出兩種基于空間統(tǒng)計性理論的功率重構(gòu)方法，并利用秦山第二核電廠3號機(jī)組和大亞灣核電站1號機(jī)組的實(shí)測功率分布進(jìn)行方法驗(yàn)證，由驗(yàn)證結(jié)果得到以下結(jié)論。

1) 兩種基于空間統(tǒng)計性理論的方法均能精確地進(jìn)行堆芯功率重構(gòu)，重構(gòu)精度滿足工程的要求；

2) 兩種方法均有很快的計算速度，滿足在線監(jiān)測的實(shí)時性要求；

3) 在實(shí)際的工程應(yīng)用中，兩種方法均可單獨(dú)使有，并可取代目前常用的耦合系數(shù)法進(jìn)行功率分布重構(gòu)。

在實(shí)際的反應(yīng)堆在線監(jiān)測應(yīng)用中，可直接獲取三維的堆內(nèi)探測器測量數(shù)據(jù)，且本工作提出的兩種方法均可拓展到三維，下一步工作將是兩種方法直接在三維功率分布重構(gòu)中的應(yīng)用，具體可采取下面兩種思路：

1) 通過三次樣條函數(shù)擬合[11]進(jìn)行軸向功率分布重構(gòu)，每一軸向?qū)拥膹较蚬β手貥?gòu)則采用本工作提出的兩種方法；

2) 將普通克里金方法和卡爾曼濾波方法直接用于三維功率分布重構(gòu)，其計算流程與二維徑向功率分布重構(gòu)一致。

參考文獻(xiàn)：

[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]李富. 重構(gòu)堆芯通量分布的諧波綜合法及其診斷應(yīng)用[D]. 北京：清華大學(xué)，1994.

[5]王常輝. 諧波展開法及其在反應(yīng)堆堆芯功率在線監(jiān)測中的應(yīng)用研究[D]. 西安：西安交通大學(xué)，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]王勇標(biāo). 克里金算法的改進(jìn)[D]. 荊州：長江大學(xué)，2012.

[8]翟進(jìn)乾. 克里金插值方法在煤層分布檢測中的應(yīng)用研究[D]. 太原：太原理工大學(xué)，2008.

[9]馬寨璞. 海洋流場數(shù)據(jù)同化方法與應(yīng)用的研究[D]. 杭州：浙江大學(xué)，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.