放射性煙羽擴散反問題解模型的初步研究

施加松，常蕓芬，吳涢暉，李發(fā)明

(防化研究院，北京 102205)

放射性煙羽擴散反問題是利用輻射環(huán)境監(jiān)測數(shù)據(jù)反解放射性煙羽特征問題的總稱(本文簡稱反問題)。由于解決反問題的技術(shù)方法可用于信息深度融合、擴散參數(shù)擬合、源項反演、核危害趨勢預(yù)測等方面，因此在輻射防護、輻射環(huán)境保護、核與輻射應(yīng)急等領(lǐng)域，關(guān)于這些方法和模型的研究一直是一個熱點，尤其是2011年福島核事故以后，反解核事故源項的技術(shù)研究達到了新的高度，大量研究方法涌現(xiàn)[1-7]，如最優(yōu)插值法(OI)[8-10]、遺傳算法(GA)[2-3]、模擬退火法(SA)[5]、人工神經(jīng)網(wǎng)絡(luò)法(ANN)[7]、卡爾曼濾波法(KF)及由其發(fā)展得到的擴展卡爾曼濾波法(EKF)和集合卡爾曼濾波法(EnKF)[4,6]、逆擴散模擬技術(shù)(IDT)[11]、貝葉斯推導(dǎo)算法[12-13]。這些方法各有其特點，在不同的應(yīng)用領(lǐng)域各有其優(yōu)勢。美國國家大氣釋放咨詢中心(NARAC)在充分考察多種優(yōu)化算法后認(rèn)為，啟發(fā)式算法(包括GA、SA、ANN等)、逆向擴散算法和貝葉斯推導(dǎo)算法在解決源項反演問題方面都具有良好的應(yīng)用條件。相對于其他算法，遺傳算法具有很好的收斂性，在計算精度要求不是很高的情況下，計算時間少、魯棒性高，其缺點是編程實現(xiàn)較復(fù)雜、依賴參數(shù)較多，易產(chǎn)生早熟收斂現(xiàn)象。本文擬通過優(yōu)化的遺傳算法來構(gòu)建放射性煙羽反問題求解模型，并通過模擬數(shù)據(jù)和風(fēng)洞數(shù)據(jù)的反解問題對模型進行驗證，從而為擴散參數(shù)擬合和源項反演等反問題性質(zhì)的工程應(yīng)用提供一種高效而精確的技術(shù)方法。

1 放射性煙羽擴散反問題數(shù)學(xué)建模

(1)

其中，θ為反解參數(shù)集。

當(dāng)然也可對方程(1)進行一些細(xì)致的處理，如加權(quán)重、取對數(shù)等[14]。這種基于最小二乘法思想的構(gòu)建方法對數(shù)據(jù)要求較高，對于誤差分布滿足高斯分布的情況，計算結(jié)果一般較好，但對于其他分布或由于監(jiān)測數(shù)據(jù)中有不可靠數(shù)據(jù)而造成存在非常大偏差的情況，計算結(jié)果會很不理想。實際上，由于監(jiān)測數(shù)據(jù)的誤差和質(zhì)量問題與煙羽模型的計算誤差之間的耦合關(guān)系，使得整體誤差分布很難滿足高斯分布，此時采用方程(1)作為目標(biāo)函數(shù)，其求解結(jié)果很難保證。為解決該問題，本文通過Steiner[15]提出的最頻值理論來構(gòu)建目標(biāo)函數(shù)：

(2)

其中：n為監(jiān)測值的個數(shù)；ε為分散系數(shù)，表征誤差分布分散程度；k為無量綱系數(shù)，當(dāng)誤差分布服從柯西分布時k=1，當(dāng)分布未知時k=2，當(dāng)分布中存在尖峰形態(tài)(如杰佛里斯分布)時k=3，本文取k=2。ε可通過方程(3)的迭代計算獲得：

(3)

目前，放射性煙羽擴散模型主要有直線高斯煙羽模型、煙團模型、隨機擴散模型等，本文采用直線高斯煙羽模型作為擴散模型，其數(shù)學(xué)表達式如下：

(4)

其中：C為下風(fēng)向放射性煙羽的濃度；(x0,y0)為釋放點位置；H為釋放高度；Q為釋放率；λ為放射性核素的衰變常量；u為釋放高度為H時的風(fēng)速；(σy,σz)為下風(fēng)向距離x-x0處橫向和垂直方向的擴散參數(shù)。

至此，放射性煙羽擴散反問題已完全轉(zhuǎn)換為求解目標(biāo)函數(shù)(2)的最小值的數(shù)學(xué)問題。

2 遺傳算法實現(xiàn)

遺傳算法是由美國Holland于1975年首先提出的，該算法是通過模擬達爾文生物進化論的自然選擇過程搜索最優(yōu)解。遺傳算法已廣泛應(yīng)用于組合優(yōu)化、機器學(xué)習(xí)、信號處理、自適應(yīng)控制和人工生命等領(lǐng)域，是現(xiàn)代有關(guān)智能計算中的關(guān)鍵技術(shù)。

遺傳算法從代表問題潛在解集的一個種群開始，按照適者生存和優(yōu)勝劣汰的原理，逐代演化產(chǎn)生出越來越好的近似解。在每一代演化過程中，需根據(jù)個體的適應(yīng)度大小選擇個體(選擇算子)，并借助于基因進行組合交叉(交叉算子)和變異(變異算子)，產(chǎn)生出代表新的解集的種群。這個過程將導(dǎo)致種群像自然進化一樣，后代種群較前代更適應(yīng)環(huán)境，末代種群中的最優(yōu)個體可作為問題近似最優(yōu)解。遺傳算法其實是一種搜索最優(yōu)解的思想，它沒有給出如何選擇、如何交叉、如何遺傳，甚至沒有給出采用什么形式來編碼(形成種群中的個體)，具有非常寬泛的適用性。為使其更好地適用于放射性煙羽擴散的反問題，采用Fortran語言編寫適用于方程(4)反解問題的遺傳算法。具體如下。

1) 編碼

以反解參數(shù)作為基因，反解參數(shù)集形成染色體(種群中個體)。

2) 適應(yīng)度選取

將式(2)的倒數(shù)作為適應(yīng)度，則min(F2(θ))的求解問題即變成最大值求解問題。

3) 選擇算子

為在充分體現(xiàn)適者生存的進化法則的同時，也兼顧個體廣泛性，本文采用了一種優(yōu)化處理方法，具體如下。

(1) 通過式(5)計算個體在下一代群體中生存的期望數(shù)目：

(5)

其中：N為個體數(shù)；Fi為第i個個體的適應(yīng)度；int函數(shù)表示取整；Ni為第i個個體進入下一代遺傳選擇操作的次數(shù)，通過此方法可確定∑Ni個個體進入下一代。

4) 交叉算子

為加快收斂性，本文采用精英交叉策略對選擇出的下一代個體進行交叉操作。精英交叉策略的基本思想如下。

(1) 首先選擇1個精英e，在此處即為適應(yīng)度最大的個體。

(2) 對于每個個體mi，均隨機產(chǎn)生1個大于0小于1的數(shù)p，如果滿足p

本文對精英交叉概率pc進行改進，提出了一種自適應(yīng)的精英交叉概率，即每個個體可根據(jù)其適應(yīng)度來調(diào)整其概率選取值，具體公式如下：

(6)

此處，F(xiàn)max取1。通過方程(6)可看出，對于適應(yīng)度較大的個體，其交叉的概率變小了，這樣做的好處是保護優(yōu)秀基因的個體，進一步加快收斂的速度。

5) 變異算子

變異算子是維持種群多樣性、避免過早收斂到局部最優(yōu)值(即“早熟”)的核心操作。本文采用當(dāng)前基因值為中心、求解域為方差的高斯隨機抽樣的方法來獲取基因突變后的個體。其中變異概率同樣采用自適應(yīng)的方法動態(tài)調(diào)整，其公式如下：

(7)

6) 終止條件

遺傳算法終止一般有以下幾種方法：(1) 采用設(shè)定最大代數(shù)終止算法運行；(2) 判定群體的收斂程度終止算法運行；(3) 通過算法離線性能和在線性能判斷是否終止算法。本算法采用第1種和第2種的混合方法，設(shè)定最大代數(shù)為20 000代，或適應(yīng)度值為1時結(jié)束搜索。計算程序流程圖示于圖1。

圖1 計算程序流程圖Fig.1 Computational flow chart

3 實驗驗證

實驗驗證的目的是驗證該算法在放射性煙羽擴散反問題領(lǐng)域應(yīng)用的可行性和可靠性，為進一步研究提供參考經(jīng)驗。驗證采用2種方法：模擬實驗驗證和風(fēng)洞實驗驗證。其中模擬實驗驗證是通過模擬實驗數(shù)據(jù)來驗證該算法解決放射性煙羽擴散反問題的可行性；風(fēng)洞實驗驗證是通過風(fēng)洞實驗數(shù)據(jù)來驗證算法的可靠性。

3.1 模擬實驗驗證

通過方程(4)產(chǎn)生1組實驗數(shù)據(jù)，釋放點坐標(biāo)為(0,0)，監(jiān)測點在x方向(西東方向)坐標(biāo)50～4 000 m和y方向(南北方向)坐標(biāo)-1 000～1 000 m矩形區(qū)域內(nèi)以x方向80 m、y方向40 m的間隔均勻分布。釋放核素選取137Cs，其他模擬參數(shù)列于表1。

表1 模擬實驗參數(shù)Table 1 Simulation experiment parameter

當(dāng)擴散參數(shù)的選取滿足Briggs擬合曲線[16](標(biāo)準(zhǔn)曲線)時，可模擬計算出監(jiān)測點的放射性濃度，其渲染圖如圖2所示，這是一個標(biāo)準(zhǔn)的放射性煙羽擴散的模擬結(jié)果(或稱理論結(jié)果)，本文用該數(shù)據(jù)反解釋放和氣象參數(shù)來驗證算法的可行性。

1) 反解參數(shù)

從方程(4)的直線高斯煙羽模型可知，釋放參數(shù)與氣象參數(shù)一起反解，會出現(xiàn)重解問題。因此，本文固定風(fēng)速和風(fēng)向來反解釋放點坐標(biāo)(x0,y0)、釋放高度H、釋放率Q和大氣穩(wěn)定度s(s為1～6的整數(shù)，s=1時大氣強不穩(wěn)定，即穩(wěn)定度等級為A，s=2時大氣不穩(wěn)定，即穩(wěn)定度等級為B，…，s=6時大氣穩(wěn)定，即穩(wěn)定度等級為F)。因此，反解數(shù)為5個，其矢量表達形式為(x0,y0,H,Q,s)。

2) 反解結(jié)果

模擬實驗中不同模擬監(jiān)測點的濃度差別很大，最小濃度為2.39×10-29Bq/m3，最大濃度達2.79×106Bq/m3，為快速計算，忽略影響較小的監(jiān)測點數(shù)據(jù)，使用大于1 000 Bq/m3的數(shù)據(jù)進行反解，結(jié)果列于表2。

表2中的適應(yīng)值反映了擴散模型與模擬數(shù)據(jù)的適應(yīng)程度，若等于1則為完全適應(yīng)。從表2可看出，適應(yīng)值非常接近1，這與本次模擬實驗監(jiān)測數(shù)據(jù)是由擴散模型計算得到的事實相符，說明該算法在解決放射性煙羽擴散反問題方面是可行的。

圖2 模擬實驗數(shù)據(jù)渲染圖Fig.2 Contour map of simulation experiment data

表2 模擬實驗反解結(jié)果Table 2 Reverse solution result of simulation experiment

3.2 風(fēng)洞實驗驗證

通過1組風(fēng)洞實驗數(shù)據(jù)來驗證反解算法的可靠性。風(fēng)洞實驗的等效參數(shù)列于表3，其中的風(fēng)速為釋放高度處的風(fēng)速。風(fēng)洞實驗數(shù)據(jù)的渲染圖示于圖3，其中xg為下風(fēng)向距離(本文為由西向東方向)，y為垂直下風(fēng)向距離。本文對監(jiān)測數(shù)據(jù)進行歸一化處理，即將濃度除以釋放率，并定義為擴散因子，單位為s/m3。

表3 風(fēng)洞實驗等效設(shè)計參數(shù)Table 3 Design parameter of wind tunnel experiment

1) 反解參數(shù)

實驗設(shè)計的風(fēng)向是西風(fēng)(風(fēng)向角270°)，但由圖3可看出，實驗煙羽中心線偏離了實驗設(shè)計方案，因此，在反解參數(shù)中要考慮風(fēng)向(φ)的變化。但對于煙羽模型，不同的風(fēng)向會出現(xiàn)重解的問題，因此必須單獨反解。另外，由于風(fēng)洞實驗的擴散特性未知，所以加入擴散參數(shù)σy和σz的擬合公式系數(shù)的反解，擬合公式如下：

σy=ayxbyσz=azxbz

(8)

其中：(ay,by)和(az,bz)分別為橫向(y方向)和垂直方向(z方向)擴散參數(shù)的擬合系數(shù)。

圖3 風(fēng)洞實驗數(shù)據(jù)渲染圖Fig.3 Contour map of wind tunnel experiment data

所以，本文考慮2層反解：第1層，反解風(fēng)向和擴散參數(shù)；第2層，用反解出的風(fēng)向和擴散參數(shù)來反解釋放點坐標(biāo)(x0,y0)、釋放高度H和風(fēng)速u。因此，第1層的反解參數(shù)為5個，其矢量表達形式為(φ,ay,az,by,bz)，其求解條件為：

φ∈0°,360°

ay>0,by∈0,20

az>0,bz∈0,20

其中，by和bz取上限的目的是避免計算搜索時浮點溢出，可根據(jù)計算機軟硬件環(huán)境調(diào)整。

第2層的反解參數(shù)為4個，其矢量表達式為(x0,y0,H,u)，其求解條件為：

H>0

u>0

2) 反解結(jié)果

(1) 風(fēng)向和擴散參數(shù)的反解結(jié)果

首先用所有大于0的實驗數(shù)據(jù)(全部數(shù)據(jù))來反解，結(jié)果列于表4。從適應(yīng)值看，擴散模型與數(shù)據(jù)的適應(yīng)性并不理想，從數(shù)據(jù)分布圖也可看出有些數(shù)據(jù)不滿足高斯煙羽分布，為此去除大于1.0×10-5和小于1.0×10-8的數(shù)據(jù)(裁剪數(shù)據(jù))，重新反解，結(jié)果列于表5。

從表5可看出，適應(yīng)性有了較大提高，風(fēng)向影響不大，但擴散規(guī)律變化較大。為比較其關(guān)系，對擴散參數(shù)與下風(fēng)向的距離作圖，并與中性條件(D類大氣穩(wěn)定度)的Briggs擬合的大氣擴散參數(shù)曲線進行比較，結(jié)果示于圖4。由圖4可見，利用全部數(shù)據(jù)反解出的擴散參數(shù)曲線與Briggs的曲線變化趨勢較接近，其中橫向擴散在近距離(幾百米以內(nèi))時較接近城市地形規(guī)律，遠距離時較接近標(biāo)準(zhǔn)地形，而垂直擴散參數(shù)在相同的下風(fēng)向距離時較Briggs曲線大。

表4 基于所有大于0的實驗數(shù)據(jù)的風(fēng)向和擴散參數(shù)反解結(jié)果Table 4 Inverse solution result of wind direction and diffusion parameters with all experimental data greater than 0

表5 基于部分?jǐn)?shù)據(jù)的風(fēng)向和擴散參數(shù)反解結(jié)果Table 5 Inverse solution result of wind direction and diffusion parameters with prune data

圖4 橫向(a)和垂直方向(b)擴散參數(shù)擬合曲線 Fig.4 Transverse (a) and vertical (b) diffusion parameter fitting curve

從以上反解結(jié)果看，風(fēng)向結(jié)果受數(shù)據(jù)的裁剪與否影響不大，但對于擴散參數(shù)，數(shù)據(jù)的多少和數(shù)據(jù)的質(zhì)量對結(jié)果影響很大。綜合起來，全部數(shù)據(jù)反解的結(jié)果較好，本文采用其反解出的風(fēng)向和大氣擴散參數(shù)來反解釋放點、釋放高度和風(fēng)速。

(2) 釋放參數(shù)的反解結(jié)果

同理，首先用所有大于0的實驗數(shù)據(jù)來反解，再用裁剪后的數(shù)據(jù)進行反解，結(jié)果分別列于表6、7。從反解結(jié)果看，二者有一定的差別，其中釋放高度的差別較大，通過數(shù)據(jù)裁剪后，這種差別有一定的改善，這也可從適應(yīng)值變大看出，由于裁剪后的監(jiān)測數(shù)據(jù)與煙羽模型符合程度變高，所以反解的結(jié)果也隨著改善，這說明反解結(jié)果的好壞對煙羽模型和監(jiān)測數(shù)據(jù)的符合程度有較高的要求。

表6 基于所有大于0的實驗數(shù)據(jù)的釋放參數(shù)反解結(jié)果Table 6 Inverse solution result of release parameters with all experimental data greater than 0

表7 基于部分?jǐn)?shù)據(jù)的擴散參數(shù)反解結(jié)果Table 7 Inverse solution result of release parameters with prune data

4 總結(jié)

本文采用最頻值理論來建立目標(biāo)函數(shù)，并與遺傳算法相結(jié)合，構(gòu)建了一種放射性煙羽擴散反問題求解模型；采用Fortran語言編寫相應(yīng)的程序，其中為提高遺傳算法的搜索效率，對遺傳算法中的選擇、交叉和變異算子進行了優(yōu)化；最后利用模擬實驗數(shù)據(jù)和風(fēng)洞實驗數(shù)據(jù)進行反問題求解，以驗證該模型的可行性和可靠性。驗證結(jié)果表明：1) 該模型應(yīng)用范圍很廣，可用于不同數(shù)據(jù)類型的多參數(shù)的反問題求解；2) 在擴散模型與監(jiān)測數(shù)據(jù)適應(yīng)性較高的情況下，該模型精度很高，但當(dāng)適應(yīng)性降低時，解的精度也降低，這是由于高斯煙羽模式本身在很多具體情況下會帶來很大誤差。