郝玉潔， 郭子萱， 胡 兵， 唐昌兵， 王浩煜， 辛 勇

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

1 引 言

核反應(yīng)堆的燃耗深度[1]即單位質(zhì)量核燃料釋放的核裂變能量，其中乏燃料(即在反應(yīng)堆中經(jīng)受過(guò)輻射照射的核燃料)的燃耗深度是衡量核電站經(jīng)濟(jì)性和核燃料輻照性能的關(guān)鍵指標(biāo).測(cè)量乏燃料的燃耗深度可通過(guò)測(cè)量其釋放的γ射線(xiàn)計(jì)數(shù)[2]來(lái)實(shí)現(xiàn).在離散測(cè)量γ射線(xiàn)計(jì)數(shù)過(guò)程中，由于受到周?chē)c(diǎn)的干擾，探測(cè)器的讀數(shù)并不精準(zhǔn).這會(huì)對(duì)后續(xù)計(jì)算帶來(lái)誤差和困難.為了克服這個(gè)問(wèn)題，有必要通過(guò)建立數(shù)學(xué)模型、運(yùn)用數(shù)學(xué)方法來(lái)減小這種探測(cè)誤差，重構(gòu)實(shí)驗(yàn)數(shù)據(jù).

在實(shí)驗(yàn)數(shù)據(jù)的去噪處理中,均值濾波算法[3-4]和中值濾波算法[5-7]是常用的兩種方法，它們對(duì)不同的噪聲有不同的去噪特性. 多年來(lái),人們一直致力于改進(jìn)濾波器,并取得了明顯效果,參見(jiàn)文獻(xiàn)[5，7].近年來(lái)，隨著小波理論的發(fā)展，很多學(xué)者將小波變換應(yīng)用于信號(hào)去噪[8-9].這種方法雖然去噪效果較好，但容易丟失細(xì)節(jié)信息，導(dǎo)致邊緣模糊.最近，鄧和劉[10]在形態(tài)學(xué)噪聲濾除器基礎(chǔ)上提出了一種基于形態(tài)學(xué)成分分析(MCA)和K奇異值分析(K-SVD)的去噪方法，可以一定程度上解決邊緣模糊問(wèn)題.

在本文考慮的實(shí)驗(yàn)數(shù)據(jù)重構(gòu)問(wèn)題中，誤差源于離散探測(cè)點(diǎn)之間的相互干擾，不是服從某種隨機(jī)分布的噪聲，因而不能直接應(yīng)用上述去噪方法，需要建立新的重構(gòu)模型.本文首先建立了探頭掃描數(shù)據(jù)的正向計(jì)算模型，將工程問(wèn)題轉(zhuǎn)化為求解一種不定的非線(xiàn)性方程組[11].該方程組有無(wú)窮多個(gè)解，需要建立優(yōu)化問(wèn)題[12]形式的重構(gòu)模型，進(jìn)而求解其特定意義下的最優(yōu)解.本文建立了兩種重構(gòu)模型，第一種是根據(jù)觀測(cè)得到的相對(duì)誤差與參數(shù)之間的線(xiàn)性關(guān)系構(gòu)造優(yōu)化目標(biāo)函數(shù)；第二種根據(jù)Lagrange乘子法[13-15]建立重構(gòu)模型.數(shù)值實(shí)驗(yàn)結(jié)果表明，兩種重構(gòu)模型得到的重構(gòu)結(jié)果都比探測(cè)值更接近真實(shí)的精確值，相較而言第二種重構(gòu)模型的重構(gòu)結(jié)果比第一種更加精確，且具有可以直接把探測(cè)值作為迭代初值的優(yōu)點(diǎn).

本文結(jié)構(gòu)安排如下：我們?cè)诘诙?jié)中建立實(shí)驗(yàn)結(jié)果正向計(jì)算的矩陣模型；第三節(jié)分別基于優(yōu)化目標(biāo)函數(shù)和Lagrange乘子法建立兩種重構(gòu)模型；第四節(jié)通過(guò)數(shù)值實(shí)驗(yàn)對(duì)重構(gòu)效果進(jìn)行驗(yàn)證；第五節(jié)總結(jié)本文的工作.

本文考慮的離散測(cè)量γ射線(xiàn)計(jì)數(shù)實(shí)驗(yàn)場(chǎng)景如下：被探測(cè)物體為一塊矩形板，探測(cè)器從板的正上方進(jìn)行散點(diǎn)式的掃描探測(cè).假設(shè)有n1行n2列共n=n1×n2個(gè)探測(cè)點(diǎn)，將第i個(gè)探測(cè)點(diǎn)的真實(shí)值記為xi，此點(diǎn)的探測(cè)值記為yi.點(diǎn)的排列順序如圖1所示.

基于離散測(cè)量γ射線(xiàn)計(jì)數(shù)的特點(diǎn)[1-2]，我們對(duì)探測(cè)實(shí)驗(yàn)做以下基本假設(shè)：

(i) 探測(cè)點(diǎn)行與行之間的距離相同，列與列之間的距離相同，且列距不小于行距;

(ii) 目標(biāo)點(diǎn)xi周?chē)奶綔y(cè)點(diǎn)會(huì)影響目標(biāo)點(diǎn)的探測(cè)值yi，這種影響與兩點(diǎn)之間的距離負(fù)相關(guān)，并且距離相同的點(diǎn)的影響因子是相同的；

(iii) 目標(biāo)點(diǎn)的探測(cè)值yi只受到探測(cè)點(diǎn)本身xi和探測(cè)點(diǎn)左方xi-1,右方xi+1,上方xi-n2,下方xi+n2,左上方xi-n2-1,右上方xi-n2+1,左下方xi-1+n2,右下方xi+1+n2共9個(gè)點(diǎn)的影響，所有9個(gè)點(diǎn)的影響因子之和為1.

圖1 探測(cè)點(diǎn)的順序Fig.1 The order of detective points

圖2 目標(biāo)點(diǎn)的探測(cè)值可以簡(jiǎn)化為附近9個(gè)點(diǎn)(包括該點(diǎn))的真實(shí)值的線(xiàn)性組合

根據(jù)以上假設(shè)，在探測(cè)器掃描目標(biāo)點(diǎn)xi時(shí)，探測(cè)值yi受到目標(biāo)點(diǎn)上下兩個(gè)點(diǎn)的影響因子是相同的，記為a；探測(cè)值受到目標(biāo)點(diǎn)左右兩個(gè)點(diǎn)的影響因子是相同的，記為b；探測(cè)值受到掃描點(diǎn)左上，左下，右上，右下四個(gè)點(diǎn)的影響因子相同，記為c. 則目標(biāo)掃描點(diǎn)的影響因子為1-2a-2b-4c.探測(cè)值yi為附近9個(gè)點(diǎn)真實(shí)值的線(xiàn)性組合

yi=(1-2a-2b-4c)xi+

a(xi-n2+xi+n2)+b(xi-1+xi+1)+

c(xi-n2-1+xi-n2+1+xi-1+n2+xi+1+n2)

(1)

根據(jù)假設(shè)(ii),距離目標(biāo)點(diǎn)越近，影響因子越大，因而參數(shù)應(yīng)滿(mǎn)足a≥b>c>0.

記x=(xi)i=1,2,…,n∈Rn是每一個(gè)探測(cè)點(diǎn)的真實(shí)數(shù)據(jù)組成的列向量，y=(yi)i=1,2,…,n∈Rn是每個(gè)探測(cè)點(diǎn)的觀測(cè)數(shù)據(jù)組成的列向量.我們得到矩陣形式的數(shù)學(xué)模型

A(a,b,c)x=y

(2)

n階方陣A可以寫(xiě)成如下n1×n1塊的分塊三對(duì)角稀疏矩陣

其中C和D都是n2階的三對(duì)角稀疏矩陣,

C=a*diag(ones(n2))+

c*diag(ones(n2-1),-1)+

c*diag(ones(n2-1),1),

D=(1-2a-2b-4c)*diag(ones(n2))+

b*diag(ones(n2-1),-1)+

b*diag(ones(n2-1),1).

進(jìn)一步,對(duì)模型矩陣A化簡(jiǎn)，將其中的未知參數(shù)a、b、c分離出來(lái)，以使比較清晰地觀察A與a、b、c的關(guān)系.簡(jiǎn)單分析后有A=I-(aA1+bA2+cA3)，其中I為n階單位矩陣，A1、A2、A3為n階對(duì)稱(chēng)正定的矩陣.A1、A3可以寫(xiě)成n1×n1塊的分塊三對(duì)角稀疏矩陣，A2可以寫(xiě)成n1×n1塊的分塊對(duì)角稀疏矩陣，

E=2*diag(ones(n2))-

diag(ones(n2-1),-1)-

diag(ones(n2-1),1),

F=-diag(ones(n2-1),-1)-

diag(ones(n2-1),1).

通過(guò)前面建立的矩陣模型，我們已經(jīng)將這個(gè)工程問(wèn)題轉(zhuǎn)化成如下的數(shù)學(xué)問(wèn)題：

Find (x;a,b,c), s.t.A(a,b,c)x=y

(3)

這是一個(gè)不定的非線(xiàn)性方程組，有無(wú)窮多解. (y;0,0,0)是其中一個(gè)解，代表真實(shí)值就等于實(shí)驗(yàn)值，與實(shí)際實(shí)驗(yàn)中的相互干擾不符.因此我們需要構(gòu)建出優(yōu)化問(wèn)題形式的重構(gòu)模型，進(jìn)而尋找符合實(shí)際需求的非(y;0,0,0)的最優(yōu)解.

3.1 重構(gòu)模型Ⅰ

由于y=x-(aA1+bA2+cA3)x,則

‖x-y‖2=‖(aA1+bA2+cA3)x‖2，

(4)

圖與a2+b2+c2的比值Fig.3 The ratio of and a2+b2+c2

s.t.Ax=y,a≥b>c>0

(5)

在上述模型中，參數(shù)k需要事先確定.研究人員可以根據(jù)探測(cè)場(chǎng)景和自己的經(jīng)驗(yàn)選取合適的參數(shù)值.如果能夠利用工程手段得到若干組探測(cè)數(shù)據(jù)y和精確數(shù)據(jù)x，根據(jù)已有的若干組參考數(shù)據(jù)x和y，用最小二乘法擬合

可以訓(xùn)練得到最優(yōu)的模型參數(shù)k.

利用數(shù)據(jù)通過(guò)最小二乘訓(xùn)練得到參數(shù)k時(shí)，優(yōu)化函數(shù)在[x;a,b,c]處的取值并不能保證嚴(yán)格為0，但在[y;0,0,0]處是嚴(yán)格為0的，并且為目標(biāo)函數(shù)的一個(gè)全局最小值點(diǎn).所以，如果將[y;0,0,0]作為初值，那么程序并不會(huì)迭代.為了克服這個(gè)問(wèn)題，我們考慮給y加上一個(gè)小擾動(dòng)t，即將[y+t;0,0,0]作為初值.

3.2 重構(gòu)模型Ⅱ

第一種模型的構(gòu)造思路非常簡(jiǎn)單直接，但它不能將實(shí)驗(yàn)數(shù)據(jù)[y;0,0,0]直接作為重構(gòu)模型求解的迭代初值，而我們真正關(guān)心的是最終重構(gòu)結(jié)果，a，b，c只是引入的中間變量.為此我們基于Lagrange乘子法建立第二種重構(gòu)模型.

由于我們只有一個(gè)等式約束

x-(aA1+bA2+cA3)x=y

(6)

受Lagrange乘子法的啟發(fā)，我們可以將其視為KKT條件：

?xL(x;a,b,c)=

x-y-aA1x-bA2x-cA3x=0

(7)

于是有Lagrange乘子函數(shù)

(8)

其中K為常數(shù).對(duì)于乘子為a，b，c的拉格朗日乘子函數(shù)可以構(gòu)造很多種優(yōu)化問(wèn)題，我們選擇如下的最簡(jiǎn)單形式：

(9)

則根據(jù)KKT條件，優(yōu)化問(wèn)題(9)的解一定滿(mǎn)足式(3).

關(guān)于模型參數(shù)ki(i=1,2,3)的選取，我們有以下定理.

定理3.1若(x*;a*,b*,c*)是問(wèn)題(9)的解，且a*,b*,c*均為正數(shù)，則

證明 若(x*;a*,b*,c*)是問(wèn)題(9)的解，則(x*;a*,b*,c*)滿(mǎn)足以下KKT條件.

(10)

上述定理告訴我們k1，k2，k3的理論選取方式，但由于x*是重構(gòu)模型需要求解的未知向量，所以這三個(gè)參數(shù)需要在求解重構(gòu)模型之前用其他方式選定.類(lèi)似于重構(gòu)模型Ⅰ，研究人員可以根據(jù)探測(cè)場(chǎng)景和自己的經(jīng)驗(yàn)選取合適的參數(shù)ki.如果能夠借助工程手段得到若干組探測(cè)數(shù)據(jù)y和精確數(shù)據(jù)x，我們還可以根據(jù)這已有的若干組參考數(shù)據(jù)x和y，使用最小二乘法擬合x(chóng)TAix-kiyTAiy=0訓(xùn)練得到模型參數(shù)k1，k2，k3.

4 數(shù)值算例

本節(jié)中我們通過(guò)數(shù)值實(shí)驗(yàn)來(lái)檢驗(yàn)上述兩種重構(gòu)模型的效果.基于離散測(cè)量γ射線(xiàn)計(jì)數(shù)的一般規(guī)律，我們首先按以下方式隨機(jī)生成50組模擬數(shù)據(jù)x：矩形板上離散探測(cè)點(diǎn)最外一圈的取值位于區(qū)間[40，100]，再往內(nèi)一圈數(shù)值位于區(qū)間[100，200]，再往內(nèi)圈數(shù)值位于區(qū)間[500，600].

圖4 一組x數(shù)據(jù)的大致形狀

圖5 兩種重構(gòu)方法的相對(duì)誤差及其與探測(cè)值的相對(duì)誤差

經(jīng)過(guò)最小二乘訓(xùn)練，得到此次實(shí)驗(yàn)中重構(gòu)模型Ⅰ的參數(shù)為k=0.188 826 971,重構(gòu)模型Ⅱ的參數(shù)為k1=1.030 744 286,k2=1.024 777 026，k3=1.025 237 126.重構(gòu)效果見(jiàn)圖5.可以看出，兩種重構(gòu)方法都有效，而第二種重構(gòu)方法效果比第一種更好.

5 結(jié) 論

針對(duì)離散測(cè)量γ射線(xiàn)計(jì)數(shù)中的實(shí)驗(yàn)數(shù)據(jù)重構(gòu)問(wèn)題，本文建立了探頭掃描數(shù)據(jù)的正向計(jì)算模型.通過(guò)將工程問(wèn)題轉(zhuǎn)化為求解一種不定的非線(xiàn)性方程組的數(shù)學(xué)問(wèn)題，本文針對(duì)此種非線(xiàn)性方程組建立兩種優(yōu)化問(wèn)題形式的重構(gòu)模型.數(shù)值實(shí)驗(yàn)結(jié)果表明，兩種模型得到的重構(gòu)結(jié)果都比探測(cè)器的探測(cè)值更接近真實(shí)的精確值，相比較而言第二種模型的重構(gòu)結(jié)果比第一種更加精確，且允許直接把探測(cè)值作為迭代初值.

在一般的離散掃描實(shí)驗(yàn)中，離散探測(cè)點(diǎn)之間都存在局部的相互干擾，干擾方式可以被本文的正向計(jì)算模型近似.所以本文建立的重構(gòu)方法，也適用于大多數(shù)領(lǐng)域的離散掃描實(shí)驗(yàn)數(shù)據(jù)重構(gòu)，具有廣泛的適用性和工程實(shí)用價(jià)值.