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

(防化研究院，北京 102205)

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

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

(1)

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

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

(2)

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

(3)

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

(4)

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

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

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

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

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

1) 編碼

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

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

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

3) 選擇算子

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

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

(5)

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

4) 交叉算子

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

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

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

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

(6)

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

5) 變異算子

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

(7)

6) 終止條件

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

圖1 計(jì)算程序流程圖Fig.1 Computational flow chart

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

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

3.1 模擬實(shí)驗(yàn)驗(yàn)證

通過方程(4)產(chǎn)生1組實(shí)驗(yàn)數(shù)據(jù)，釋放點(diǎn)坐標(biāo)為(0,0)，監(jiān)測點(diǎ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í)驗(yàn)參數(shù)Table 1 Simulation experiment parameter

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

1) 反解參數(shù)

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

2) 反解結(jié)果

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

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

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

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

3.2 風(fēng)洞實(shí)驗(yàn)驗(yàn)證

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

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

1) 反解參數(shù)

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

σy=ayxbyσz=azxbz

(8)

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

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

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

φ∈0°,360°

ay>0,by∈0,20

az>0,bz∈0,20

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

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

H>0

u>0

2) 反解結(jié)果

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

首先用所有大于0的實(shí)驗(yàn)數(shù)據(jù)(全部數(shù)據(jù))來反解，結(jié)果列于表4。從適應(yīng)值看，擴(kuò)散模型與數(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)向影響不大，但擴(kuò)散規(guī)律變化較大。為比較其關(guān)系，對(duì)擴(kuò)散參數(shù)與下風(fēng)向的距離作圖，并與中性條件(D類大氣穩(wěn)定度)的Briggs擬合的大氣擴(kuò)散參數(shù)曲線進(jìn)行比較，結(jié)果示于圖4。由圖4可見，利用全部數(shù)據(jù)反解出的擴(kuò)散參數(shù)曲線與Briggs的曲線變化趨勢較接近，其中橫向擴(kuò)散在近距離(幾百米以內(nèi))時(shí)較接近城市地形規(guī)律，遠(yuǎn)距離時(shí)較接近標(biāo)準(zhǔn)地形，而垂直擴(kuò)散參數(shù)在相同的下風(fēng)向距離時(shí)較Briggs曲線大。

表4 基于所有大于0的實(shí)驗(yàn)數(shù)據(jù)的風(fēng)向和擴(kuò)散參數(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)向和擴(kuò)散參數(shù)反解結(jié)果Table 5 Inverse solution result of wind direction and diffusion parameters with prune data

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

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

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

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

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

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

4 總結(jié)

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