堆芯功率分布重構(gòu)的兩種空間統(tǒng)計(jì)性算法

彭星杰，李 慶，王 侃

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

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

1 普通克里金方法

1.1 方法介紹

(1)

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

(2)

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

(3)

其中：μ為拉格朗日乘子；c(xi，xj)為觀測(cè)點(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)計(jì)性理論認(rèn)為估計(jì)值與觀測(cè)值的相似程度是通過(guò)點(diǎn)對(duì)距離相關(guān)量來(lái)度量的。點(diǎn)對(duì)距離相關(guān)量只與觀測(cè)場(chǎng)間的相互距離有關(guān)，而與它們的絕對(duì)位置無(wú)關(guān)。常用的距離相關(guān)量模型有指數(shù)模型、高斯模型等。根據(jù)具體情況可采用不同的相關(guān)函數(shù)，本工作采用二階自回歸模型[8](式(4))來(lái)描述距離相關(guān)量，其數(shù)學(xué)上的優(yōu)越性已得到了充分理論證明[10]。

(4)

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

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

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

(5)

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

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

2 卡爾曼濾波方法

2.1 方法介紹

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

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

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

(6)

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

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

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

Za=Zb+K(Y0-HZb)

(7)

其中，K為增益矩陣。

K=BHT(HBHT+R)-1

(8)

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

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

(9)

式(9)中的H變?yōu)槌跏脊烙?jì)物理場(chǎng)Zb處的雅可比矩陣。

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

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

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

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

B=σ2c

(10)

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

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

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

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

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

(10)

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

3) 觀測(cè)算符H

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

4 結(jié)論與展望

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

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

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

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

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

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

2) 將普通克里金方法和卡爾曼濾波方法直接用于三維功率分布重構(gòu)，其計(jì)算流程與二維徑向功率分布重構(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]王常輝. 諧波展開(kāi)法及其在反應(yīng)堆堆芯功率在線監(jiān)測(cè)中的應(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]. 荊州：長(zhǎng)江大學(xué)，2012.

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

[9]馬寨璞. 海洋流場(chǎng)數(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.