彭星杰，李 慶，王 侃

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

為保證核電站運行的經(jīng)濟性和安全性，必須對與燃料元件和包殼屏障完整性有關(guān)的安全參數(shù)進行監(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é)者對這一問題進行了大量研究。Combustion Engineering公司開發(fā)的耦合系數(shù)法(CECOR)[1]計算流程簡單且在工程上得到廣泛應(yīng)用，但其是建立在全堆分布變化不大的假設(shè)基礎(chǔ)上。Webb[2]對耦合系數(shù)法進行了改進，發(fā)展出Lagrange乘子法，其計算思想和計算精度與CECOR的差別不大。Lee等[3]發(fā)展了最小二乘法，通過求解中子擴散方程與探測器響應(yīng)方程形成的超定方程組來獲得重構(gòu)通量分布，但這種算法必須嵌入到堆芯擴散計算程序中，目前在商用堆芯核設(shè)計軟件中還未得到應(yīng)用。李富[4]開發(fā)了諧波綜合法，利用中子擴散方程高階諧波的有限階線性組合來表示真實的堆芯通量分布。王常輝[5]通過建立諧波數(shù)據(jù)庫實現(xiàn)了諧波綜合法在線重構(gòu)堆芯實時功率分布，但這種方法需預(yù)先計算好若干參考工況的高階諧波，前期工作量巨大。Bryson等[6]將卡爾曼濾波方法應(yīng)用于堆芯功率分布的重構(gòu)，但未考慮觀測量的歸一化問題。在其他領(lǐng)域，如何通過有限的測量值得到物理場分布的估計也是一熱點問題[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)為觀測點xi與xj之間的空間協(xié)方差函數(shù)；c(x，xj)為待估值點x與xj之間的空間協(xié)方差函數(shù)。

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

空間協(xié)方差函數(shù)代表了一種距離相關(guān)量，空間統(tǒng)計性理論認(rèn)為估計值與觀測值的相似程度是通過點對距離相關(guā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為兩點間的距離。

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

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

(5)

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

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

2 卡爾曼濾波方法

2.1 方法介紹

卡爾曼濾波方法[8]是一種典型的數(shù)據(jù)同化方法。數(shù)據(jù)同化方法是一種廣義上的最小二乘方法，求解原有理論計算值與實際測量值所構(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)后的堆芯功率分布，之后還需進行一次歸一化。

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

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

B=σ2c

(10)

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

數(shù)值實驗表明，當(dāng)L取為1.5個組件長度，使用卡爾曼濾波方法對秦山第二核電廠3號機組和大亞灣核電站1號機組功率進行重構(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è)該比值為固定值1%，即α取為0.01。

3) 觀測算符H

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

3 數(shù)值驗證

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

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

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

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

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

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

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

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

3) 對未放置探測器的組件，普通克里金方法和卡爾曼濾波方法得到的功率重構(gòu)精度在少部分組件內(nèi)并不高于耦合系數(shù)法，這是因為重構(gòu)相對誤差存在隨機性，有正有負，而為了保證整體重構(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號機組重構(gòu)方法的相對誤差對比

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

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

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

3) 對大亞灣核電站1號機組的驗證結(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)進行求和，明顯滿足在線監(jiān)測要求。對于卡爾曼濾波方法，在計算機硬件條件為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號機組和大亞灣核電站1號機組的實測功率分布進行方法驗證，由驗證結(jié)果得到以下結(jié)論。

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

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

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

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

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

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

參考文獻：

[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). 克里金算法的改進[D]. 荊州：長江大學(xué)，2012.

[8]翟進乾. 克里金插值方法在煤層分布檢測中的應(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.