基于粒子濾波算法的源項反演及不確定性分析

鄭子豪 陳春花 朱婧嫻 陳黎偉 王世鵬

1（中國科學院合肥物質科學研究院 合肥 230031）

2（中國科學技術大學 合肥 230026）

3（合肥師范學院計算機學院 合肥 230601）

“十四五”規劃提出了以安全穩妥為前提，推動沿海核電項目建設，并提高核能綜合利用的目標［1］。核能的發展以核安全為前提，核應急又是核安全最后一道屏障，源項分析是否精確決定了核應急行動的有效性。核事故發生后，源項信息通常難以直接獲得，一般可利用廠外監測數據以及核素擴散模型推演得到源項信息，此過程為源項反演［2］。

近些年許多新的算法應用于源項反演，例如布谷鳥搜索算法、計算流體力學算法、粒子群優化算法和截斷總體最小二乘變分法。Wang等［3］利用改進的布谷鳥搜索算法預測未知釋放源的數量和坐標，結果表明，噪聲信噪比越大、監測范圍越廣時，對釋放源的追蹤結果越精確；Kovalets等［4］利用計算流體力學對未知泄漏源的位置進行了反演，相比結合了貝葉斯理論或卡爾曼濾波的大氣彌散模型，其結果的誤差減少了一個數量級；黎岢等［5］利用粒子群優化算法對氣態137Cs 進行了反演，并且設置了286個觀測點，在釋放率穩定的情況下，可以很好地追蹤到真值，非穩態情況估計結果較差；劉蘊等［6］通過建立截斷總體二乘變分核事故源項反演數值計算模型，在釋放率穩定的條件下，結果相對真值有20%的誤差。以上方法對背景條件設置較為理想，普遍側重研究穩態的源項釋放，對監測器數量和模擬數據設置過多，對噪音設定單一且不變。為了讓反演過程更加符合現實，還應針對變化的監測噪音、有限的監測數據和非穩態的源項釋放率進行研究。

考慮到監測噪音復雜多變，引入粒子濾波來降低噪音。由于粒子濾波算法需要將源項釋放狀態和源項擴散狀態作為輸入，本文在僅有兩組監測數據的條件下，結合大氣擴散模型，一組數據用于研究核事故初期源項瞬時釋放率的遞歸關系，將其作為粒子濾波的狀態方程；另一組數據作為測量數據，用作粒子濾波的測量方程。以此建立了在核事故發生初期，源項瞬時釋放率多變、監測器有限和多重噪音干擾條件下的基于粒子濾波的源項反演方法。最后對整個反演過程進行不確定性分析，并證實可以此方法可以得到較為精確的反演結果。

1 原理與方法

本文根據源項的不同釋放狀態和放射性核素在大氣中的擴散規律，建立了基于粒子濾波的實時反演方法。同時針對事故源項的特點，建立瞬時釋放率狀態方程，且對監測噪音變化規律進行了模擬假設。

1.1 源項狀態模型

核事故發生后，源項瞬時釋放量的狀態方程［7］為（1）式，測量方程為（2）式。

式中：Yk為監測向量；Hk(Q1，Q2，...，Qk)為觀測函數，代表煙團釋放至大氣后的擴散的規律；rk在本文中是監測儀器的測量相對誤差和其他干擾項的統一表達。

1.2 粒子濾波算法

粒子濾波基于蒙特卡羅模擬方法實現遞推貝葉斯濾波估計，對于系統的過程噪聲和量測噪聲沒有的限制［8］。粒子濾波實時追蹤源項的步驟如下四個步驟。

步驟一：設定初始值，分配初始權重。人為地設定一個初始值Q1代表釋放量的初始數據，采集n個粒子樣本，每個粒子疊加噪音q1，見（3）式。

用式（7）可以得到t= 2 時刻Q2的最小均方估計。

下一個時刻t=k+ 1，轉到步驟三中對粒子重采樣，再到步驟四中更新權值，并以此迭代。

1.3 大氣擴散模型

余琦等［9］模擬了穩定流場和非穩定流場中煙羽模型和煙團模型擴散變化，其結果表明，煙團模型的精度更高，同時對于連續性釋放，單位時間內釋放的煙團越多，精度越高。為了研究核事故初期源項瞬時釋放量，本文選擇高斯煙團模型，并設定每秒鐘釋放一個煙團，故每個煙團所含的放射總活度為源項的瞬時釋放率。假設泄漏源的坐標為原點，分別以正東和正北方向為X軸和Y軸，垂直于地面向上為Z軸，建立空間直角坐標系。開始泄漏時間T=0，當T=t秒時第i=t個煙團開始釋放，源項此時瞬時釋放率為Qt，此煙團經過τ秒后，用式（12）計算得到(x，y，z)處的活度濃度。

式中：Qi為需要反演的數據，Bq/s；σx(i)、σy(i)和σz(i)為大氣擴散參數，表示第i=t個煙團在水平和垂直方向的擴散參數，其數值可參考Briggs修正公式［10］；ux和uy為風速在X軸和Y軸上的分量，m/s；h0為煙團的有效釋放高度，m。為了使擴散更好地模擬各種氣象條件下的擴散過程，本文還考慮了煙團擴散時的干沉降、放射性核素的種類及其衰變修正［11］。

由于釋放率Qi與其他客觀條件相互獨立，可將其與時間τ進行變量分離，將點(x，y，z)處關于時間τ的函數設為Singlepuff(τ)。Singlepuff(τ)為高斯函數，反映了煙團擴散的規律，將式（12）簡化表示為式（13）。

當泄漏時間為T=t時，所有已釋放的煙團在(x，y，z)的活度濃度累計貢獻為式（14）。

1.4 瞬時釋放率狀態方程的構建

Singlepuff(τ)與源項的釋放率無關，在時間上呈現高斯分布，其特點是在某一時間內有著最大值，隨著此時間增大或減少，其數值下降快且趨近于0，且函數在時間上的積分快速收斂。根據公式（2）和（14），將觀測函數Hk(Q1，Q2，...，Qk)改寫成矩陣形式，見式（15）。

式中：rt為監測噪音。由式（15）可以得到，狀態矩陣[Q1，Q2，...，Qt]與Singlepuff(t)之間隨著時間呈現著有序的對應關系。因此可以通過監測數據Yt變化的趨勢判斷源項瞬時釋放率變化狀態，由此建立狀態方程。

（1）當源項釋放率連續恒定，可以認為每秒釋放的煙團總活度相等，由于設函數Singlepuff(t)在時間的積分等于HK，根據式（15）有式（16）。

式中：Q和HK均為固定值，當監測數據長時間在某個數值穩定波動時，可以認定此時釋放率穩定不變，此時狀態方程見式（17）。

（2）當泄漏源的瞬時釋放率為線性變化時，設每秒鐘釋放率變化k，由式（15）有方程組（18）。將式（18）兩式相減可得（19）式。

當t大于ux/x和uy/y時，第一個煙團開始往遠離監測點的方向擴散，第t+1時刻的煙團還沒有開始隨風擴散，公式（19）前兩項約等于0，得（20）式。

可以得到在tx時的平均釋放率為公式（23）。

此方法計算出的釋放率具有tx的偏移，具體的偏移范圍與氣象狀況和監測器相對地理位置有關。將平均釋放率--Qx代入高斯煙團模型正演其擴散規律，并與實際的監測數據比較后得到偏移量tx。

圖1（a）是核事故發生后600 s 內的監測數據；圖1（b）是計算出的平均釋放率，可見，其變化的趨勢與監測數據相同；圖1（c）是平均釋放率正演的結果與監測數據的對比，由于使用平均釋放率作為瞬時釋放率，正演數據相比監測數據平滑，根據圖像可以判斷偏移時間tx；圖1（d）表示時間修正后的釋放率。這樣求得的平均釋放率因監測誤差而在真值附近波動，對釋放率進行擬合，用擬合后的曲線作為狀態方程。

圖1 釋放率連續變化時的狀態方程建立流程：(a)核事故發生后600 s內的監測數據；(b)計算得到的平均釋放率；(c)平均釋放率正演的結果與監測數據的對比；(d)時間修正后的釋放率與擬合曲線Fig.1 The process of establishing the equation of state when the release rate changes continuously:(a)monitoring data within 600 s after the nuclear accident;(b)calculated average release rate;(c)comparison between the forward modeling results of average release rate and monitoring data;(d)time corrected release rate and fitting curve

1.5 監測噪音

一般認為，監測器的實時監測數據包含著此刻源項由大氣擴散后在周圍空氣中的活度濃度，并疊加監測器的觀測誤差。Morino等［12］模擬了福島核事故中131I 和137Cs 因干濕沉積累計到地面上的活度，其分別占釋放總活度的8.2%和11.7%。Hu［13］模擬計算了不同的源項釋放率的沉積過程，其結果表明，不同的釋放率、天氣、風向等因素，對放射性核素沉積到地面所污染的面積、污染的分布及其最大沉積活度均有影響。因此，具體某時刻的監測數據還應該包含核事故發生后、此時刻之前放射性核素的歷史沉積。由于本文源項瞬時釋放率是未知的，根據文獻［13］的結果可以對歷史沉積的總活度進行定性分析：在同一時刻，距離下風向和側風向越遠，其沉降的總活度越低。針對放射性物質沉積對監測數據的影響，在源項反演模擬過程中，設置為隨著時間逐漸增大的均勻噪聲，與監測器固有誤差進行疊加后，即為公式（2）中的rk。

2 仿真實驗和不確定性分析

2.1 仿真實驗設定

實驗模擬在核事故發生后t時刻，131I以3種不同的釋放狀態持續泄漏至環境中，瞬時釋放率分別為：穩定狀態，Q= 1 × 1010Bq/s；線性增長狀態，Q(t) = 108+t× 109Bq/s；正弦變化狀態，Q(t) = 109×(2 + sin(t/100)) Bq/s。131I 的 半 衰 期設置為8.07 d，僅布置兩個監測位點。

實驗的氣象條件設定：大氣穩定度為C類；風速為2 m/s；風向為東風；無降雨。地理信息設定：以泄漏源為坐標原點、正東方向為X軸、正北方向為Y軸、垂直于地面為Z軸建立空間直角坐標系，監測器A坐標為（50，20，10），監測器B的坐標為（100，20，10）；煙團有效釋放高度為10 m。反演時間設定：核事故發生后t=1 s 開始，監測器每秒鐘傳回監測數據，泄漏時間360 s 后開始反演。監測噪聲設定：兩個監測器相對測量誤差均為20%，監測器A 在t=60 s 后每隔60 s 疊加0～2%的均勻噪聲，監測器B 在t=120 s 后每隔60 s 疊加0～0.5%的均勻噪聲。粒子濾波的粒子數設置為100。

考慮到監測器B距離泄漏點更遠，其受到的噪音干擾相對更小，因此實驗中利用監測器B的監測數據設定狀態方程，過程誤差協方差矩陣用過程噪聲相對誤差10%計算，利用監測器A 的監測數據作為粒子濾波的觀測值。

2.2 計算狀態方程

根據§1.4瞬時釋放率狀態方程構建的方法。

（1）當釋放率連續穩定時，測量得到的釋放率如圖2（a）所示，設定釋放狀態方程為Q= 1.05 ×1010Bq/s。（2）當釋放率線性增長時，其測量釋放率如圖2（b）所示，t=1 s 時釋放率的計算測量值為Q(1) = 1.8 × 109Bq/s，將其設為釋放狀態方程的初始值，根據圖像狀態方程設定為Q(t) =Q(t?1) + 1.02 × 109Bq/s。（3）當釋放率正弦變化時，其測量釋放率如圖2（c）所示，經過時間修正后得到前300 s的平均釋放率如圖2（d）所示。對曲線進行擬合，選擇圖2（d）中最高處的點，并以此點為中心選擇另外兩個相對稱的點，計算經過這3個點的二次函數，并以此作擬合曲線，有Q(t)=?3.89×104×t2+1.26×107×t+2×109Bq/s，用此函數作為狀態方程。

圖2 由監測數據計算而得的釋放率：(a)釋放率穩定不變時的計算值；(b)釋放率線性變化時的計算值；(c)釋放率正弦變化時的計算值；(d)釋放率正弦變化時的修正數據和擬合曲線Fig.2 Calculated release rate from the monitoring data:(a)calculated value when the release rate is stable;(b)calculated value when the release rate changes linearly;(c)calculated value when the release rate changes sinusoidally;(d)correction data and fitting curve of sinusoidal change of release rate

2.3 實驗結果

三種釋放狀態的測量值和濾波值以及它們相對于真值的誤差如圖3 所示，根據圖3（a）、（c）和（e），濾波值相對測量值波動范圍減小，整體更加平滑；根據圖3（b）、（d）和（f）的誤差曲線，濾波值最終在真值附近波動，相對真值的誤差約5%，最大不超過10%。當釋放率為線性增長時，由于狀態方程設定初始值相比真值較大，前60 s濾波值誤差呈現逐漸減小的趨勢，在60 s后誤差在3%波動。由于監測器的相對誤差設定在20%，說明粒子濾波可以有效地減少誤差，反演結果優于計算測量值。

圖3 穩定(a)、(b)，線性變化(c)、(d)，和正弦變化(e)、(f)的源項釋放粒子濾波反演結果Fig.3 Particle filter inversion results of stable release(a),(b),linear release(c),(d),and sine release(e),(f)

2.4 粒子數和過程噪聲相對誤差對濾波結果的影響

粒子濾波算法的粒子數量代表著每次隨機采樣的次數，當采樣的范圍已被確定時，采樣次數越多，粒子落在真值處的概率越大。還需討論在特定的過程噪聲相對誤差下，粒子數設定對濾波結果的影響。

當釋放率線性增長、過程噪聲相對誤差為10%時，分別設置粒子數為50、100、1 000。粒子濾波反演如圖4所示，不同的粒子數有著幾乎相同的結果。由于較小的過程噪聲將粒子每次采樣限定在很小的范圍內，當過程噪聲相對誤差不大于監測器相對測量誤差時，粒子數不是減少反演誤差的主要因素。

圖4 過程噪聲相對誤差為10%時粒子數為50、100、1000的粒子濾波反演結果Fig.4 Particle filter inversion results with particle numbers of 50,100 and 1 000 when the process noise covariance matrix is 10%

設置過程噪聲相對誤差為50%，此時過程噪聲大于監測器噪聲，粒子濾波器相信監測數據，此時粒子濾波的作用不再是過濾噪聲，而是追蹤測量值。分別設置粒子數為50、100、1 000，粒子濾波反演結果的誤差如圖5 所示，當粒子數為1 000 時，反演值與測量值擬合結果優秀，幾乎沒有波動。以粒子數為1 000 時為基準，如圖5（a）所示，當粒子數為50 時，反演結果與測量值的擬合較差；當粒子數為100時，反演結果在測量值附近波動，波動范圍不小于20%。

圖5 以1 000為基準，粒子數分別為50(a)和100(b)時的結果相對誤差Fig.5 Based on 1 000,the relative error of the results when the number of particles is 50(a)and 100(b)

根據實驗結果可以看出，當狀態方程比監測值精確時，粒子數對反演結果影響較??；當狀態方程誤差比監測值更大時，濾波結果偏向于與監測值擬合，粒子數設置越多擬合程度越好。綜合兩種情況考慮，粒子數設置偏向于越多越好，但增加粒子數量會導致計算時間增長，因此粒子數設置需要同時兼顧狀態方程誤差和計算機性能。由于本文實驗設定狀態方程的數據來源于受到噪音干擾較小的監測器，故設置100個粒子時，即可獲得較佳的結果。

2.5 初始值對濾波結果的影響

粒子濾波算法中，狀態方程的初值是人為設定的。為了探討人為設定的初始值對濾波結果的影響，分別以源項初始釋放率真值的0.1%、1%、10%、100%、1 000%和10 000%作為粒子濾波的初始值，當源項釋放率線性增加時，設置Q(t) =108+t× 109Bq/s，反演結果相對真值的誤差，如圖6（a）所示，當初始值設定小于100%初值真值時，粒子濾波反演值可以快速追蹤到真值；當初始值設定大于100%真值時，粒子濾波反演值趨向于真值，隨著初始值增大，粒子濾波接近真值的速度越來越慢。

當釋放率線性減少時，設置Q(t) = 1011?t× 109Bq/s，分別設定源項初始釋放率真值的0.1%、1%、10%、100%、1 000%和10 000%作為粒子濾波的初始值。反演結果誤差的對數如圖6（b）所示，當設定的初始值遠小于或者遠大于真值時，結果的誤差極大；當初始值在真值附近時，相對誤差的對數在±1之間，即±10%。

圖6 線性增加(a)和線性減少(b)時不同初值對粒子濾波反演的影響（彩色見網絡版）Fig.6 Influence of different initial values with linear increase(a)and linear decrease(b)on particle filter inversion(color online)

根據以上分析可以得出結論，在粒子濾波反演源項中，初始值的設定不存在偏大或者偏小的偏好規律，為了得到最優的反演結果，應該始終以實際的監測數據作為參考。

2.6 監測器誤差對粒子濾波的影響

現實中不同的監測儀器、不同的環境、不同的使用方式等因素，均可導致不同的測量誤差［14］。為了探討監測器的相對測量誤差對濾波結果的影響，分別以5%、20%和50%為例，當釋放率恒定時，分別計算這三種不同的誤差所對應的濾波結果。

當相對監測誤差為5%時，測量釋放率與粒子濾波反演結果如圖7 所示。濾波值與測量值幾乎重合，由于預測初始值相對真值誤差比監測值大，前60 s的反演值相比測量值稍大，在180 s后由于均勻噪聲越來越大，反演值開始小于測量值，并且相比測量值更加接近于真值。由結果可知，當監測器相對測量誤差小于狀態方程的過程噪聲相對誤差時，濾波器相信監測數據，反演結果與測量釋放率大致相同。

圖7 監測相對誤差為5%時粒子濾波反演結果Fig.7 Particle filter inversion results when the observation noise covariance matrix is 5%

當相對監測誤差為20%時，根據圖3（a）和（b），粒子濾波反演值相比測量值更加平滑，精確度更高；當相對監測誤差為50%時，反演結果和測量值相對于真值的誤差如圖8所示，與圖3（b）的結果類似，粒子濾波的誤差波動相比測量值更小，更加靠近真值。由結果可知，當監測器相對測量誤差大于狀態方程的過程噪聲相對誤差時，粒子濾波器相信狀態方程，反演結果不被增大的相對監測誤差影響。

圖8 監測誤差為50%時的粒子濾波反演誤差與測量誤差Fig.8 Particle filter inversion error and measurement error when the noise covariance matrix is 50%

2.7 討論

考慮到現實中監測器布置有限，實驗只設定了兩組監測數據，分別用于建立狀態方程和作為粒子濾波的數據輸入，粒子濾波反演的結果相對真值誤差更小、數據波動幅度更低，說明粒子濾波在此反演模型中可以很好的過濾噪聲。

狀態方程的準確性是粒子濾波能有效地過濾噪聲的前提。當設定的狀態方程的誤差過大時，可以設置較大的過程噪聲，并設置大量的粒子數，濾波值才可以很好地與監測數據擬合；當設定狀態方程誤差過大，但設置較小的過程噪聲時，濾波值結果無法與監測數據擬合。實驗中狀態方程建立源自于監測數據，狀態方程近似于對監測結果的擬合，因此其相對誤差與監測數據的相對誤差近似。根據粒子濾波原理，當狀態方程和觀測方程的相對誤差較為接近時，濾波后的結果才可以很明顯地降低噪音。綜合以上實驗結果和不確定性分析可以得知，如果有其他的途徑獲得誤差更小的狀態方程，源項反演的濾波結果將越為準確，可見實施對源項變化趨勢的研究也應該是源項反演的重要工作。

3 結論

本文研究了在源項釋放率連續多變、監測器數量有限以及監測噪音多變的情況下，基于粒子濾波的源項反演方法。結果表明，僅使用兩組監測數據，選擇受噪音干擾較小的監測數據設定源項釋放率的狀態方程和初值；選擇受噪音干擾較大的監測數據設定為粒子濾波觀測數據；設置粒子濾波的粒子數100，當監測器均誤差為20%時，反演結果相對真值的誤差在±5%附近波動，最大不超過10%。此方法針對多種釋放場景有著類似的結果，原則上可以用于核事故發生早期對源項瞬時釋放率的反演計算。

作者貢獻說明 鄭子豪論文初稿撰寫，論文審閱與修訂，實驗仿真模擬和實驗結果分析；陳春花研究內容總體設計；朱婧嫻實際調查研究與項目管理；陳黎偉模型測試；王世鵬研究方法指導。全體作者均已閱讀并同意最終的文本。