基于福島核事故監測數據的輻射劑量場重構方法研究

徐 瑤，陳曉雷，黃廣偉，吳涢暉，陳 林

(軍事科學院防化研究院，北京 102205)

準確評定核事故后危險區等級是維護國家核安全和國防安全的重要課題，有助于緩解和遏制污染源，同時預測污染物時空傳輸和擴散情況，為防護、洗消和行動提供決策基礎。當前對于核擴散的預測主要依靠放射性煙羽在大氣中的擴散原理和數學擴散模型[1]，在掌握源項信息的基礎上建立物理模型，計算輻射劑量場的分布，從而對核危害的發展態勢做出評估，常用的模擬計算方法有蒙特卡羅方法[2]、點核積分方法[3]等。實際上，基于大氣擴散模式的模型具有一定偏差，事故點地形、氣象等條件復雜多變，且在核事故發生早期源項信息往往很難準確獲得，導致計算結果往往會有較大誤差。隨著信息技術快速發展，核輻射傳感器信息化能力得到大幅提升，將核輻射傳感器廣泛布設于監測區域內，可以實現劑量數據的實時采集和動態監測，為基于實測數據進行輻射劑量場重構奠定了技術基礎。然而由于核事故后污染范圍極大，很難獲取整個污染區的精細數據，只能利用有限采樣數據對輻射劑量場進行重構，在這種情況下，通過少量實測數據重構核輻射劑量分布的研究變得尤為重要?？臻g插值方法可以有效地基于少量樣本實現大范圍輻射劑量場重構：Thevenaz等人[4]研究表明插值計算在輻射劑量場重建中具有重要作用；李廷等人[5]利用反距離權重方法估算三維空間輻射劑量率場，實現真實場景下的光子劑量場與核素活度場耦合；張敏等人[6]基于MCNP程序模擬特定空間內存在中子源時，獲取了輻射場內網格節點的濃度值，通過立方樣條插值恢復網格線的濃度值，建立了一個基于網格節點的反演模型；賽雪等人[7]將徑向基函數散亂數據插值方法應用于γ輻射場，對具有不同分布對稱性的γ輻射場進行了數據重構；陳雷等人[8]采用克里金插值模型進行沉降量內插，獲得可靠的估算數據，但選取的樣本點的數量較少，影響克里金插值模型整體精度?？梢钥闯?，沒有一種算法適用于所有的場景，不同的插值方法對計算結果影響不同。為比較不同插值方法重構輻射劑量場的精度，本文利用反距離權重插值法、克里格插值法、線性三角網插值法以及徑向基函數插值法研究散亂節點的二維輻射劑量場的重構問題，探討不同插值方法對輻射劑量值的插值效果，進一步通過交叉驗證對比分析各種方法插值結果的精度。

1 問題和算法描述

在大量的應用領域中，通常面臨利用某一個函數來描述測量值的任務。對于這個問題一般有插值和回歸兩種解決辦法，插值法更適用于假定數據正確的環境。對于需要尋找多元解析函數來描述高維空間中多維數據的任務，如天氣預報風云高壓圖、山川河流等高線測量等，二維插值是最常用的方法，輻射劑量場的重構也屬于這類問題，其基本理論假設是空間位置上越靠近的點越可能具有相似的特征值，距離越遠的點其特征值相似的可能性越小。實際問題中節點數據網格往往是散亂無規律的，為插值增加了難度。散亂節點網格插值的數學過程描述如下:

已知n個節點(xi,yi,zi),i=1,2,…,n，其中(xi,yi)互不相同。需要構造一個二元函數z=f(x,y)，通過全部已知節點，即f(xi,yi)=zi，再用f(x,y)計算插值，即z*=f(x*,y*)。

通過上述方法可以對不足或缺失數據進行補充插值，實現數據網格化，并且可以進一步獲得等值線以顯示數據的空間分布特征。

1.1 反距離權重插值法

反距離權重插值法(inverse distance weighted，IDW)是由 Shepard[9]提出的一種全局插值算法，其利用插值點周圍樣本點值的信息來計算插值點的值，以插值點與樣本點間的距離為權重進行加權平均，離插值點越近的樣本點賦予的權重越大，即[10]：

(1)

式中，n為樣本點個數；Z*(x*,y*)為點(x*,y*)處的插值計算值；λi為插值點對于樣本點的權重，與插值點與樣本點之間的距離相關，Zi(xi,yi)為在樣本點(xi,yi)處的測量值。

權重λi的表達式為：

(2)

式中，di為插值點(x*,y*)與各樣本點(xi,yi)之間的歐氏距離；α為指數項，反映樣本點之間的權重，通常取值1～3，α值影響插值圖像的光滑度，α值越大，插值的光滑性越差。通過選取不同α值進行網格化分析，當α取3時得到一個較滿意的結果，誤差最小，因此本文中α取3。

1.2 克里金插值法

克里金插值(kriging)法又稱空間自協方差最佳插值法[11]，綜合考慮了監測點和被估計點的之間的相對位置關系和各監測點之間的相對位置關系，是一種基于變異函數理論對變量值進行無偏、最優估計的空間數據插值方法[12]。克里金插值法包括普通克里金法、簡單克里金法、指示克里金法、協同克里金法、泛克里金法等。

設Z*(x*,y*)是區域化變量，滿足二階平穩或本征的，其數學期望為常數，Zi(xi,yi)是定義在樣本點(xi,yi)上的一組離散的信息樣品數據，在任一非樣本點(x*,y*)處的克里金估計值為[13]：

(3)

克里金插值法的原則需要保證估計值無偏且估計方差最小，根據拉格朗日原理，求解得出n個權值系數應滿足以下條件[13]：

(4)

式中，μ是拉格朗日乘子；c是協方差函數。通過求解方程組(4)得到權重系數λi，再根據式(3)可得出最終的插值結果。求解上述過程的關鍵一步是選擇合適的變異函數模型，克里金提供了不同的理論變異函數模型，常用的有球體模型、高斯模型、指數模型、線性模型等[14]。本文通過實驗對比驗證，選擇普通克里金模型以及線性變異函數模型使函數曲線與監測點進行最優匹配。

1.3 線性三角網插值

線性三角網插值(linear triangulation interpolation)又叫線性插值，線性三角網法采用最佳Delaunay三角化方法，將離散數據點直接相連構成三角形，區域內所有三角形的邊互不相交構成一張三角形網[15]。

線性三角網插值的原理是根據區域內任意待插值點P(x,y)與其所在三角形的三個頂點之間的線性關系形成權系數，利用此權系數計算出待插值點的函數值，其插值函數為[16]：

(5)

式中，λi為點Zi(xi,yi)的權系數；Zi(xi,yi)，i=1,2,3為三角形頂點的函數值；Z*(x*,y*)為插值點P(x,y)的函數值，權系數λi須在待插值點位于三角形內的情況下進行插值。線性三角網插值方法充分利用了已知點的信息，反映數據點與其鄰近點間的拓撲連接關系。

1.4 徑向基函數插值

徑向基函數(radial basis function，RBF)是一種基于距離的函數，具有結構簡單、各向同性、尺寸無關和無網格等特點，適用于對大量離散數據的表面插值和擬合[17]。RBF插值的基本理論基礎是未知函數在高維空間中，任一點的函數值Z*(x*,y*)可以用徑向基函數的線性組合來逼近，計算公式為[18]：

(6)

式中，n為插值點個數；αi為第i個插值節點的權重；φ為徑向基函數，函數值取決于n維空間的歐式距離；di為插值點(x*,y*)與各采樣點(xi,yi)之間的歐氏距離。

徑向基函數包括五種基本核函數：多重二次曲面函數(MQF)、多重對數函數(MLF)、反多重二次曲面函數(IMQF)、薄板樣條曲面函數(TPSF)和自然三次樣條曲面函數(NCSF)[19]。核函數是影響插值精度的重要因素，選用五種核函數對樣本進行插值分析，采用標準偏差來分析不同核函數的插值效果。五種核函數插值精度從高到低排序依次為：多重二次曲面函數>自然三次樣條函數>薄板樣條函數>多重對數函數>反多重二次曲面函數，因此本文選擇的核函數為多重二次曲面函數。

2 算例計算

2.1 數據來源

本文選用的數據來自日本原子能研究開發機構(Japan Atomic Energy Agency，JAEA)[20]在日本福島第一核電站事故后發布的監測數據，事故發生地點地理位置為北緯37°25′17″，東經141° 01′57″，選用距福島第一核電站約80 km區域內的空氣劑量率數據，該數據根據空中監測調查的測量結果創建，通過轉換在每個監測點上方空氣中測得的計數率，計算得到地面上1 m處空氣劑量率數據。數據庫包括監測點的經緯度、劑量率數據、監測起止時間等信息。

2.2 輻射劑量插值實現

本次研究共選取2 000個樣本點作為已知點進行計算，這些監測點的地理位置及對應位置的劑量率值如圖1所示，隨機分布于監測范圍內。

圖1 監測點地理位置及對應劑量率值Fig.1 The geographic location of the monitoring pointand the corresponding dose rate value

該算例的計算過程通過美國Golden Software公司Surfer軟件實現，該軟件可用于制作基面圖、數據點位圖、等值線圖等，并提供多種數據網格化和數據計算方法[21]。計算過程首先需要對數據進行網格化處理，網格間距往往對結果的準確性有影響：網格間距過大導致結果不準確，網格間距過小會導致計算時間長。本文分別選取不同大小的網格間距對輻射劑量數據進行網格化，計算插值數值與測量數值之間的標準偏差，同時綜合考慮準確性以及計算機的容量和運行時間，最終選擇行列節點數為200×200進行網格化，該間距下準確度較高且耗時較小。

將非均勻分布的原始劑量數據，通過反距離權重插值法、克里金插值法、線性三角網插值法以及徑向基函數插值法計算得出200×200的網格節點數據。

依據得到的細網格數據，使用surfer軟件進行輻射劑量場的顯示。圖2是利用2 000個樣本點分別采用上述4種方法繪制生成的輻射劑量趨勢面，從圖2中可以看出，4種插值方法均能很好地呈現出輻射劑量場的空間分布狀態，呈現效果相似。在監測區域內部，相比較于克里金插值法和徑向基函數插值法，反距離權重插值法與線性三角網插值法得到的曲面峰值區域圖形更加尖銳，更易受到臨近點劑量值的影響。

圖2 基于四種方法生成的輻射劑量率趨勢面Fig.2 Radiation dose rate trend surface generated based on four methods

為了更直觀的比較網格化的結果，將四種插值方法生成的細網格數據在二維平面擬合生成輻射劑量率等值線分布，如圖3所示。從圖3中可以看出，4種方法都能較好地插值擬合出輻射劑量率的大致分布情況，且等值線趨勢相近。但是，克里金插值法和徑向基函數插值法得到的等值線分布圖較另外兩種方法在局部更加精細，比較右上角區域(經度140.7°以東，緯度37.6°以北區域)，可以看出克里金插值法和徑向基函數插值法在劑量率0.13～0.25 μSv/h區間，層次仍然比較分明。線性三角網插值法得到的等值線分布圖在整個區域內線條清晰，但其嚴格控制了實測數據的邊界，插值范圍有限，空白區域較大。而反距離權重插值法等值線比較粗糙，在局部監測區域更易出現“牛眼狀”分布，插值結果會與真實值有較大偏差，無法真實地反映出這些點位的劑量變化。

圖3 四種方法生成的輻射劑量率等值線分布圖Fig.3 Distribution map of radiation dose rate contour generated by four methods

2.3 誤差分析

為了進一步對比4種插值方法的準確性，采用交叉驗證對插值方法進行比較。將采用計算的2 000個樣本點劃分為兩部分，一部分為參與插值計算的訓練集，一部分為未參與插值計算的測試集。對于測試集的選取，劃分不同數量的樣本作為測試集對結果進行評價，測試集的樣本數量分別隨機選取1, 10, 20, 50, 100, 200, 300, 400, 500, 600, 700, 800, 900, 1 000個，即假設這些測試集點位的數據未知，用訓練集的數據使用不同的插值方法計算出測試集所在位置處的模擬值，然后再與測試集中實際的測量值進行對比，計算模擬值與測量值之間的誤差。此步驟重復進行20次，最終計算誤差的平均值。誤差評價指標主要選擇常用的標準偏差(σ)、平均絕對誤差(MAE)和均方根誤差(RMSE)[22]。標準偏差用來衡量數據自身的離散程度，平均絕對誤差用來衡量模擬值的誤差范圍，均方根誤差用來衡量模擬值與真值之間的偏差程度。這3個指標的計算方法如式(7)～(9)[22]，誤差計算結果如圖4所示。

圖4 誤差分析結果Fig.4 Error analysis result

(7)

(8)

(9)

式中，μ表示樣本總體xi的均值；fi表示模擬值；yi表示測量值；n為樣本數目。

從圖4計算結果可以看出，不同插值方法之間的標準偏差、平均絕對誤差以及均方根誤差并不相同。測試集樣本數逐漸增多，意味著訓練集樣本數減少，各插值方法精度均有所下降，σ、MAE和RMSE均隨著訓練樣本數減少而波動提升。在相同的測試集樣本數下，徑向基函數插值法的預測誤差最小，反距離權重插值法預測誤差范圍最大，線性三角網插值法和克里金插值法介于兩者之間。徑向基函數插值法的RMSE和σ均保持在1.7以下，MAE保持在0.41以下，在插值精度上優于其他3種方法，說明徑向基函數插值法的插值精度更高，預測效果更好。對插值點與測量點對比分析，4種方法產生偏差主要分布在極值附近和劑量變化明顯的地方，四周邊緣區域由于缺乏測量數據作參考，插值填充后該區域失真較嚴重，也是偏差主要分布位置。但徑向基函數插值法和克里金插值法能很好地表現出細節處的劑量變化情況，偏差相對較小，說明該兩種插值算法的誤差比較穩定。造成以上結果的原因是因為反距離權重插值法僅考慮了距離的影響，根據周圍測量點的距離權重進行賦值，忽略了各樣本點之間的空間關聯性，導致計算值容易受數據點集群的影響，易出現孤立數據點現象，反映在輻射劑量等值線分布圖上即“牛眼狀”分布。線性三角網插值法構造三角形的速率會導致計算速度慢，且無法保證曲面的光滑性要求，采樣數據點個數較少且分布不均勻的情況下重構精度不高，易出現明顯的邊界?？死锝鸩逯捣ú粌H考慮了距離權重，還考慮了各監測點之間的相對位置關系，使得求取權系數更合理，利用變異函數所揭示的測量值的內在聯系來估計插值點數值，插值精度高。徑向基函數插值法具有很強的擬合數據點、產生光滑曲面的能力，能夠突出輻射劑量的變化，在劑量數據變化較大時插值精度較高，模擬曲面能較好地保持劑量局部變化細節與趨勢。這也解釋了為什么圖3中徑向基函數插值方法得到的等值線分布圖更精細且誤差最小。

3 結論

本文基于有限數量的輻射監測數據，利用4種插值算法實現了輻射劑量場的重構，探索不同插值算法在輻射劑量場重構的可行性，并比較不同方法的精度。通過對插值結果與預測誤差的分析，得到的結論如下：本文所涉及到的插值方法均能表達出輻射劑量場的空間分布特征，根據不同插值方法的誤差統計和效果顯示比較可以發現，徑向基函數插值法無論是插值效果上還是插值精度上都優于其余插值方法，能更好地呈現出輻射劑量明顯的變化趨勢。同時，為提高插值精度，可以適當增加采樣點數量和劃分網格數量。

本研究為核事故早期源項信息不明的情況下基于實測數據構建輻射劑量分布提供一種有效的預測方法，對輻射劑量場的空間插值應用具有一定參考意義。下一步，可以考慮改進完善徑向基函數插值法，嘗試將該方法應用于其他輻射劑量場重構與分析中，同時將地形、氣候等不同因素加入到插值模型，結合以上提出輻射劑量場分布特點的基函數，以此構造出更精確且符合真實場景的插值計算結果。