王振宇，黃偉奇,*，孫 健，彭廣宇

(1.陸軍防化學(xué)院，北京 102205；2.湖南交通職業(yè)技術(shù)學(xué)院，湖南 長(zhǎng)沙 410132)

掌握輻射劑量場(chǎng)的分布對(duì)實(shí)施輻射防護(hù)行動(dòng)具有重要意義，掌握輻射劑量場(chǎng)分布的方法主要可分為計(jì)算和監(jiān)測(cè)兩大類(lèi)：計(jì)算類(lèi)有點(diǎn)核積分、蒙特卡羅、離散縱標(biāo)等方法；監(jiān)測(cè)一般指部署固定或移動(dòng)輻射測(cè)量設(shè)備來(lái)獲取輻射場(chǎng)數(shù)據(jù)。當(dāng)放射源的活度、位置等信息未知時(shí)，輻射劑量場(chǎng)分布的調(diào)查主要依靠輻射探測(cè)器監(jiān)測(cè)，監(jiān)測(cè)方法的缺點(diǎn)是監(jiān)測(cè)數(shù)據(jù)分布離散、成本高。而插值算法能在少量監(jiān)測(cè)數(shù)據(jù)的基礎(chǔ)上，通過(guò)構(gòu)建數(shù)學(xué)模型對(duì)未知輻射劑量場(chǎng)分布進(jìn)行重構(gòu)，研究插值算法有助于快速、低成本地解決輻射劑量場(chǎng)分布問(wèn)題。

插值算法是一種常見(jiàn)的數(shù)據(jù)處理方法，在各領(lǐng)域均有不同程度的應(yīng)用。在進(jìn)行輻射監(jiān)測(cè)時(shí)，合理運(yùn)用插值算法可提高計(jì)算效率。Ródenas等[1]和Wang等[2-3]分別采用線性插值和網(wǎng)函數(shù)插值完成了監(jiān)測(cè)數(shù)據(jù)呈網(wǎng)格均勻分布時(shí)的輻射劑量場(chǎng)的插值重建。但實(shí)際上監(jiān)測(cè)數(shù)據(jù)呈離散分布的情況更普遍，此時(shí)相對(duì)復(fù)雜的反距離權(quán)重法[4-6]、克里金(Kriging)法[7]和徑向基函數(shù)法[8-9]能進(jìn)行輻射劑量場(chǎng)的插值重建。此外還有一些考慮到輻射場(chǎng)特性的改進(jìn)插值方法，如張敏[10]在反距離平方加權(quán)法的基礎(chǔ)上提出的輻射劑量值占優(yōu)構(gòu)造還原法。

泊松克里金法屬于改進(jìn)的插值算法，最早由Monestiez等[11]在調(diào)查長(zhǎng)須鯨數(shù)量分布時(shí)提出，Zhao等[12]將其應(yīng)用于輻射場(chǎng)的插值重構(gòu)。實(shí)際上，泊松克里金法并沒(méi)有改變克里金法的求解方法，全局插值優(yōu)化時(shí)仍存在最優(yōu)解無(wú)法收斂的問(wèn)題，此外泊松克里金法在輻射監(jiān)測(cè)領(lǐng)域的應(yīng)用較少，參數(shù)影響效果、使用范圍與效果有待進(jìn)一步研究討論。本文嘗試在泊松克里金優(yōu)化計(jì)算中采用代理模型的方法，通過(guò)實(shí)驗(yàn)分析參數(shù)影響的效果，并用不同條件下的輻射場(chǎng)進(jìn)行驗(yàn)證，為相關(guān)研究提供參考。

1 泊松克里金算法

1.1 泊松克里金算法原理

(1)

(2)

泊松克里金算法是克里金插值算法的一種特殊應(yīng)用，主要區(qū)別是泊松克里金算法用觀測(cè)數(shù)據(jù)取代真實(shí)數(shù)據(jù)進(jìn)行無(wú)偏估計(jì)。泊松克里金算法主要應(yīng)用于真實(shí)樣本數(shù)據(jù)難以獲取的情況，其假設(shè)是觀測(cè)值與真實(shí)值之間符合泊松分布的關(guān)系。在輻射場(chǎng)測(cè)量中，由于放射性衰變規(guī)律和輻射探測(cè)器原理的限制，真實(shí)的劑量率無(wú)法測(cè)得，只能獲得其觀測(cè)值x。根據(jù)輻射場(chǎng)計(jì)數(shù)泊松分布理論，輻射探測(cè)器獲得的劑量率觀測(cè)值也服從泊松分布關(guān)系[12]，即：

x|d～Poisson(d)

(3)

故本文采用泊松克里金算法，即用樣本的劑量率觀測(cè)值xi對(duì)輻射場(chǎng)真實(shí)劑量率進(jìn)行無(wú)偏估計(jì)，估計(jì)值的表達(dá)式示于式(4)。

(4)

求解出滿足條件的一組最優(yōu)權(quán)重系數(shù)ci，即可求得輻射場(chǎng)真實(shí)劑量率的估計(jì)值。

1.2 基于代理模型的泊松克里金算法求解

克里金插值算法的權(quán)重系數(shù)方程組的求解通常采用拉格朗日乘數(shù)法，該方法只能得到局部最優(yōu)解。前期計(jì)算中發(fā)現(xiàn)，當(dāng)輻射場(chǎng)劑量率數(shù)據(jù)多時(shí)，會(huì)出現(xiàn)計(jì)算時(shí)間長(zhǎng)、最優(yōu)解無(wú)法收斂的情況。本文將克里金代理模型用于泊松克里金算法的求解。代理模型是指可替代復(fù)雜真實(shí)模型的一類(lèi)簡(jiǎn)化近似數(shù)學(xué)模型，克里金代理模型是代理模型的一種。本文在泊松克里金理論的基礎(chǔ)上對(duì)克里金代理模型進(jìn)行適當(dāng)修改，同時(shí)將模式搜索算法用于最優(yōu)解的收斂?jī)?yōu)化。

在構(gòu)建的泊松克里金代理模型中，某點(diǎn)的輻射場(chǎng)真實(shí)劑量率d可表示為回歸模型mTβ與隨機(jī)誤差e的和[13-14]：

d=mTβ+e

(5)

式中：mT為回歸模型基函數(shù)，為行向量；β為回歸系數(shù)，為列向量；e為誤差函數(shù)。同時(shí)在泊松克里金算法中，輻射場(chǎng)劑量率觀測(cè)值x服從期望為真實(shí)值d的泊松分布，根據(jù)全期望公式E(x)=E[E(x│d)]與總方差公式var(x)=E[var(x│d)]+var[E(x│d)]，可求得劑量率觀測(cè)值x的期望與方差：E(x)=mTβ、var(x)=mTβ+e。因此劑量率觀測(cè)值x也可表示為回歸模型與整體誤差之和(式(6))，其中的mTβ和e的取值與式(5)不同。

x=mTβ+e

(6)

因此，在泊松克里金代理模型中，真實(shí)值的估計(jì)值可用矩陣表示為式(7)。

(7)

式中：回歸模型的基函數(shù)mT和誤差函數(shù)e為可選擇的已知參數(shù)，基函數(shù)mT為可變化階次的多項(xiàng)式函數(shù)，后續(xù)計(jì)算中用矩陣M表示一組基函數(shù)mT；誤差函數(shù)e之間的相關(guān)性用相關(guān)系數(shù)表示，其中矩陣R表示樣本數(shù)據(jù)的相關(guān)系數(shù)；向量cT為權(quán)重系數(shù)?；貧w系數(shù)向量β與權(quán)重系數(shù)向量cT為未知參數(shù)，需單獨(dú)求解。

回歸系數(shù)β：為確定最優(yōu)參數(shù)β，可根據(jù)樣本點(diǎn)輻射劑量率觀測(cè)值X進(jìn)行最大似然估計(jì)(X為x的矩陣形式)。在模型中，輻射場(chǎng)劑量率是一個(gè)n維二元高斯分布問(wèn)題，樣本觀測(cè)劑量的似然函數(shù)可構(gòu)建為式(8)。

(8)

式中：L為似然函數(shù)；σ2為高斯分布的方差；n為坐標(biāo)向量的維度(本文取n=2)，求解式(8)后可得：

(9)

(10)

將式(9)、(10)代入似然函數(shù)(式(8))中，忽略常數(shù)項(xiàng)后可將似然函數(shù)化簡(jiǎn)為：

(11)

在最大似然估計(jì)中，似然函數(shù)(式(11))為目標(biāo)函數(shù)，相關(guān)系數(shù)矩陣的行列式|R|為自變量。通過(guò)模式搜尋算法，可逐步逼近最優(yōu)相關(guān)系數(shù)矩陣R，進(jìn)而也可求出回歸系數(shù)β。通過(guò)相關(guān)系數(shù)矩陣R也能確定相關(guān)系數(shù)模型的具體參數(shù)，利用該相關(guān)系數(shù)模型可求解待估計(jì)點(diǎn)的相關(guān)系數(shù)rT。

(12)

同樣通過(guò)拉格朗日乘數(shù)法求解后，輻射場(chǎng)真實(shí)劑量率的估計(jì)值表達(dá)式為：

(MTR-1M)-1MTR-1X

(13)

(14)

2 泊松克里金法參數(shù)影響分析

為研究泊松克里金算法中參數(shù)的影響，基于SuperMC軟件構(gòu)建一個(gè)室內(nèi)多放射源虛擬輻射場(chǎng)。該虛擬輻射場(chǎng)尺寸為40 m(長(zhǎng))×30 m(寬)，含兩個(gè)活度為3.7×108Bq 的137Cs和活度為1.85×108Bq的60Co放射源。通過(guò)計(jì)算得到距地面1 m高度處的當(dāng)量劑量率，如圖1所示。采用datasample函數(shù)從圖1中隨機(jī)均勻抽取100個(gè)數(shù)據(jù)作為樣本數(shù)據(jù)。該樣本數(shù)據(jù)用于分析泊松克里金代理模型中各參數(shù)變量對(duì)輻射場(chǎng)插值還原結(jié)果的影響。

2.1 多項(xiàng)式基函數(shù)的影響

多項(xiàng)式基函數(shù)的階次影響回歸模型擬合的整體趨勢(shì)，擬合結(jié)果通過(guò)判定系數(shù)(R2)來(lái)評(píng)價(jià)，R2是線性回歸中對(duì)估計(jì)的回歸方程擬合優(yōu)度的度量，其計(jì)算公式如式(15)所示。不同階次多項(xiàng)式的擬合結(jié)果列于表1。

(15)

式中：SSE為殘差平方和；SST為總離差平方和。

表1 回歸模型擬合結(jié)果Table 1 Regression model fitting result

由表1可知，多項(xiàng)式階次越高，判定系數(shù)R2越大，即擬合效果越好，但系數(shù)個(gè)數(shù)增加，計(jì)算更加復(fù)雜。其中二次多項(xiàng)式擬合與三次多項(xiàng)式擬合結(jié)果中R2相差較小，系數(shù)個(gè)數(shù)卻相差較大。綜合上述因素，選擇二次多項(xiàng)式為回歸模型的基函數(shù)。

2.2 誤差相關(guān)模型的影響

不同點(diǎn)位處的誤差具有相關(guān)性，這也是不同位置輻射劑量率的相互關(guān)聯(lián)在數(shù)學(xué)模型上的具體體現(xiàn)。常用的相關(guān)模型列于表2。

圖1 距地面1 m高度處的當(dāng)量劑量率Fig.1 Equivalent dose rate at height of 1 m above ground

表2 相關(guān)模型函數(shù)Table 2 Correlation model function

圖2 各相關(guān)系數(shù)模型函數(shù)曲線Fig.2 Correlation coefficient model function curve

由圖2可看出，總體上相關(guān)系數(shù)R隨距離d的增大而減小，參數(shù)θ主要影響距離d與相關(guān)系數(shù)R之間的函數(shù)關(guān)系。當(dāng)θ<0.1時(shí)，無(wú)論距離d大小，相關(guān)系數(shù)始終接近于1；當(dāng)θ增加時(shí)，呈現(xiàn)距離近處的相關(guān)系數(shù)大，距離遠(yuǎn)處的相關(guān)系數(shù)小，當(dāng)θ增加至一定值時(shí)，相關(guān)系數(shù)隨距離的增大迅速減小，且距離稍遠(yuǎn)處的相關(guān)系數(shù)為0，此時(shí)繼續(xù)增加θ對(duì)模型函數(shù)的影響較小。如在EXP、LIN、SPHERICAL、CUBIC、SPLINE模型中，θ=20的函數(shù)曲線與θ=50的函數(shù)曲線十分接近，GAUSS模型中θ=50的函數(shù)曲線與θ=100的函數(shù)曲線接近。因此在后續(xù)泊松克里金代理模型的最優(yōu)相關(guān)系數(shù)求解中，可將模式搜索算法中θ的取值限定在一定范圍內(nèi)，在EXP、LIN、SPHERICAL、CUBIC、SPLINE模型中θ∈(0.1,20)，GAUSS模型中θ∈(0.1,50)。在此基礎(chǔ)上，采用不同的相關(guān)系數(shù)模型對(duì)原輻射場(chǎng)進(jìn)行插值還原。主要參數(shù)如下：采用二階多項(xiàng)式回歸，GAUSS模型θ∈(0.1,50)，初始θ=25，其余模型θ∈(0.1,20)，初始θ=10。各模型插值還原的效果示于圖3，圖中黑點(diǎn)為樣本數(shù)據(jù)。

圖3 不同相關(guān)模型對(duì)輻射場(chǎng)插值還原效果Fig.3 Radiation field reconstruction effect of different related models

在CUBIC和SPLINE兩種模型結(jié)果中，相關(guān)系數(shù)隨距離的增加迅速減小，與其他模型相比，在同一距離處有更小的相關(guān)系數(shù)，因此同一峰上的不同數(shù)據(jù)被識(shí)別成多個(gè)單峰。EXP、LIN和SPHERICAL三種模型形成的數(shù)據(jù)不夠平滑，不符合實(shí)際輻射場(chǎng)連續(xù)場(chǎng)的特征，其中EXP模型形成的輻射峰值底部較大，不符合輻射峰較尖銳的特征，會(huì)錯(cuò)誤地增加輻射場(chǎng)區(qū)域面積。而GAUSS模型形成的峰與原始數(shù)據(jù)較一致，且平滑效果較好。此外，對(duì)比樣本數(shù)據(jù)點(diǎn)可發(fā)現(xiàn)，GAUSS模型還原的輻射場(chǎng)最高劑量率為2.5 mSv/h，明顯高于樣本數(shù)據(jù)2 mSv/h，且與原始數(shù)據(jù)趨勢(shì)一致。在實(shí)際應(yīng)用中，這種尋找樣本數(shù)據(jù)范圍外的熱點(diǎn)區(qū)域的功能十分重要。綜上所述，GAUSS相關(guān)系數(shù)模型較適用于γ射線在空氣中形成的輻射場(chǎng)插值還原。

3 實(shí)驗(yàn)驗(yàn)證

3.1 真實(shí)輻射場(chǎng)的構(gòu)建與測(cè)量

實(shí)驗(yàn)場(chǎng)地選在室內(nèi)，場(chǎng)地尺寸為15 m×20 m×4 m。以1 m為間隔在實(shí)驗(yàn)場(chǎng)地面上建立x-y坐標(biāo)系，用標(biāo)記貼在網(wǎng)格點(diǎn)處順次編號(hào)標(biāo)記。本實(shí)驗(yàn)共設(shè)2組，放射源均為137Cs。單源組放射源活度為1.746 MBq；雙源組中一號(hào)放射源活度為0.786 MBq，二號(hào)放射源活度為0.96 MBq。

采用6150AD-b閃爍體探測(cè)器(卡迪諾科技貿(mào)易有限公司，國(guó)防科技工業(yè)電離輻射一級(jí)計(jì)量站于2013年12月26日校準(zhǔn))進(jìn)行測(cè)量，探測(cè)器能量響應(yīng)范圍為23 keV～7 MeV，最大誤差為±30%。實(shí)驗(yàn)測(cè)量高度為0.5 m，與放射源之間的測(cè)量角度保持在±80°以內(nèi)，測(cè)量過(guò)程中保證儀器與放射源間無(wú)障礙物。

經(jīng)過(guò)對(duì)原始數(shù)據(jù)進(jìn)行判棄篩選，取某點(diǎn)處連續(xù)3次所讀取數(shù)據(jù)的平均值作為該點(diǎn)的最終測(cè)量劑量率。經(jīng)過(guò)數(shù)據(jù)處理，單源輻射場(chǎng)、雙源輻射場(chǎng)的三維分布和等高線圖分別示于圖4、5，其中單源輻射場(chǎng)最高劑量率為5.7 μSv/h，雙源輻射場(chǎng)2號(hào)峰劑量率為2.8 μSv/h、3號(hào)峰劑量率為3.3 μSv/h。

圖4 單源輻射場(chǎng)分布Fig.4 Radiation field of single source

圖5 雙源輻射場(chǎng)分布Fig.5 Radiation field of dual sources

3.2 基于小范圍測(cè)量數(shù)據(jù)的算法驗(yàn)證

小范圍輻射場(chǎng)數(shù)據(jù)取自上述實(shí)驗(yàn)測(cè)量數(shù)據(jù)。為研究樣本點(diǎn)數(shù)量對(duì)輻射場(chǎng)的還原效果，在樣本數(shù)據(jù)量N=20、40、60、80、100、150、200、300時(shí)，進(jìn)行插值計(jì)算并繪制等劑量率線圖。單源輻射場(chǎng)等劑量率線圖示于圖6，雙源輻射場(chǎng)等劑量率線圖示于圖7。

平均誤差E用于定量對(duì)比效果的影響，所有網(wǎng)格點(diǎn)位處誤差e的平均值為平均誤差E，誤差e計(jì)算公式如式(16)所示。為消除樣本數(shù)據(jù)隨機(jī)選取過(guò)程對(duì)誤差帶來(lái)的影響，對(duì)每種情況選取100個(gè)樣本數(shù)據(jù)，并進(jìn)行計(jì)算分析，結(jié)果示于圖8。

(16)

圖6 單源輻射場(chǎng)等劑量率線圖Fig.6 Contour map of single source radiation field

圖7 雙源輻射場(chǎng)等劑量率線圖Fig.7 Contour map of dual sources radiation field

由圖6、7可知，隨著樣本點(diǎn)數(shù)量的增加，輻射場(chǎng)劑量率熱點(diǎn)區(qū)域逐漸收斂，同時(shí)輻射劑量率峰也逐漸尖銳。圖8中，當(dāng)樣本點(diǎn)數(shù)量為10時(shí)，單源和雙源輻射場(chǎng)的平均誤差分別為27.63%和26.30%；當(dāng)樣本點(diǎn)數(shù)量增加到150時(shí)，平均誤差下降較快，分別降至4.53%和4.24%；樣本點(diǎn)數(shù)量為150～300時(shí)，平均誤差下降較慢，最終仍分別存在2.56%和2.46%的平均相對(duì)誤差。綜上，在10%的相對(duì)誤差允許范圍內(nèi)，泊松克里金算法還原輻射場(chǎng)的最小樣本量為網(wǎng)格點(diǎn)數(shù)的1/10。

3.3 基于大范圍(福島周邊)數(shù)據(jù)的算法驗(yàn)證

為驗(yàn)證該算法在大范圍輻射場(chǎng)的插值重構(gòu)效果，從日本原子能研究開(kāi)發(fā)機(jī)構(gòu)(JAEA)網(wǎng)站獲取了福島第一核電站80 km范圍內(nèi)空間劑量率的監(jiān)測(cè)結(jié)果[15]，監(jiān)測(cè)時(shí)間為2016年8—10月，如圖9a所示。在插值重構(gòu)中分別以監(jiān)測(cè)數(shù)據(jù)坐標(biāo)經(jīng)緯度范圍(東經(jīng)140.125 3°～141.037 9°、北緯36.710 6°～38.165 3°)為x、y坐標(biāo)范圍，對(duì)該區(qū)域重新劃分成200×300的網(wǎng)格。然后通過(guò)隨機(jī)函數(shù)選取600個(gè)劑量率數(shù)據(jù)點(diǎn)作為樣本數(shù)據(jù)。采用算法插值重構(gòu)的輻射場(chǎng)分布如圖9b所示。

由圖9可知，采用算法插值重構(gòu)的劑量率分布(圖9b)中峰值的分布走向與原始劑量率分布(圖9a)一致，說(shuō)明在大范圍輻射場(chǎng)應(yīng)用中，輻射場(chǎng)熱點(diǎn)區(qū)域能基本識(shí)別出來(lái)。計(jì)算顯示，其平均相對(duì)誤差為141.69%，表明熱點(diǎn)定位效果不如小范圍輻射場(chǎng)準(zhǔn)確。此外，在大范圍輻射場(chǎng)計(jì)算中，形成可靠的輻射場(chǎng)劑量分布所需的樣本點(diǎn)數(shù)量較多，分析其主要原因是福島周邊大范圍輻射場(chǎng)中放射源的位置分布不規(guī)則，放射性核素濃度分布不均勻，形成的真實(shí)輻射場(chǎng)較復(fù)雜。

圖8 平均誤差與樣本點(diǎn)數(shù)量的關(guān)系Fig.8 Influence of sample data size on error

圖9 原始劑量率分布(a)[14]和采用算法插值重構(gòu)的劑量率分布(b)Fig.9 Original dose rate distribution (a)[14] and dose rate distribution reconstructed by algorithm (b)

4 結(jié)論

本文對(duì)代理模型求解的泊松克里金算法進(jìn)行了參數(shù)分析，并通過(guò)實(shí)測(cè)輻射場(chǎng)數(shù)據(jù)、福島核電站周邊輻射場(chǎng)數(shù)據(jù)進(jìn)行了算法應(yīng)用驗(yàn)證。結(jié)果表明，該算法在小范圍簡(jiǎn)單輻射場(chǎng)重構(gòu)中效果較好，由30個(gè)監(jiān)測(cè)點(diǎn)數(shù)據(jù)構(gòu)建的輻射場(chǎng)的相對(duì)誤差在10%左右；在大范圍復(fù)雜輻射場(chǎng)重構(gòu)中效果略差，由600個(gè)監(jiān)測(cè)點(diǎn)數(shù)據(jù)構(gòu)建的輻射場(chǎng)的相對(duì)誤差為141.69%。本文計(jì)算結(jié)果初步驗(yàn)證了泊松克里金算法在不同條件下輻射劑量場(chǎng)重構(gòu)應(yīng)用中的可行性，證明該方法在快速、低成本解決未知放射源輻射場(chǎng)的重構(gòu)方面中具有一定的潛力和研究?jī)r(jià)值。