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

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

1（中國科學院合肥物質(zhì)科學研究院 合肥 230031）

2（中國科學技術(shù)大學 合肥 230026）

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

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

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

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

1 原理與方法

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

1.1 源項狀態(tài)模型

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

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

1.2 粒子濾波算法

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

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

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

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

1.3 大氣擴散模型

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

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

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

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

1.4 瞬時釋放率狀態(tài)方程的構(gòu)建

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

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

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

式中：Q和HK均為固定值，當監(jiān)測數(shù)據(jù)長時間在某個數(shù)值穩(wěn)定波動時，可以認定此時釋放率穩(wěn)定不變，此時狀態(tài)方程見式（17）。

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

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

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

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

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

圖1 釋放率連續(xù)變化時的狀態(tài)方程建立流程：(a)核事故發(fā)生后600 s內(nèi)的監(jiān)測數(shù)據(jù)；(b)計算得到的平均釋放率；(c)平均釋放率正演的結(jié)果與監(jiān)測數(shù)據(jù)的對比；(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 監(jiān)測噪音

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

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

2.1 仿真實驗設定

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

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

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

2.2 計算狀態(tài)方程

根據(jù)§1.4瞬時釋放率狀態(tài)方程構(gòu)建的方法。

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

圖2 由監(jiān)測數(shù)據(jù)計算而得的釋放率：(a)釋放率穩(wěn)定不變時的計算值；(b)釋放率線性變化時的計算值；(c)釋放率正弦變化時的計算值；(d)釋放率正弦變化時的修正數(shù)據(jù)和擬合曲線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 實驗結(jié)果

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

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

2.4 粒子數(shù)和過程噪聲相對誤差對濾波結(jié)果的影響

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

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

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

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

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

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

2.5 初始值對濾波結(jié)果的影響

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

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

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

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

2.6 監(jiān)測器誤差對粒子濾波的影響

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

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

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

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

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

2.7 討論

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

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

3 結(jié)論

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

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