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

施加松，常蕓芬，吳涢暉，李發明

(防化研究院，北京 102205)

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

1 放射性煙羽擴散反問題數學建模

(1)

其中，θ為反解參數集。

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

(2)

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

(3)

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

(4)

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

至此，放射性煙羽擴散反問題已完全轉換為求解目標函數(2)的最小值的數學問題。

2 遺傳算法實現

遺傳算法是由美國Holland于1975年首先提出的，該算法是通過模擬達爾文生物進化論的自然選擇過程搜索最優解。遺傳算法已廣泛應用于組合優化、機器學習、信號處理、自適應控制和人工生命等領域，是現代有關智能計算中的關鍵技術。

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

1) 編碼

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

2) 適應度選取

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

3) 選擇算子

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

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

(5)

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

4) 交叉算子

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

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

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

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

(6)

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

5) 變異算子

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

(7)

6) 終止條件

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

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

3 實驗驗證

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

3.1 模擬實驗驗證

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

表1 模擬實驗參數Table 1 Simulation experiment parameter

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

1) 反解參數

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

2) 反解結果

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

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

圖2 模擬實驗數據渲染圖Fig.2 Contour map of simulation experiment data

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

3.2 風洞實驗驗證

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

表3 風洞實驗等效設計參數Table 3 Design parameter of wind tunnel experiment

1) 反解參數

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

σy=ayxbyσz=azxbz

(8)

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

圖3 風洞實驗數據渲染圖Fig.3 Contour map of wind tunnel experiment data

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

φ∈0°,360°

ay>0,by∈0,20

az>0,bz∈0,20

其中，by和bz取上限的目的是避免計算搜索時浮點溢出，可根據計算機軟硬件環境調整。

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

H>0

u>0

2) 反解結果

(1) 風向和擴散參數的反解結果

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

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

表4 基于所有大于0的實驗數據的風向和擴散參數反解結果Table 4 Inverse solution result of wind direction and diffusion parameters with all experimental data greater than 0

表5 基于部分數據的風向和擴散參數反解結果Table 5 Inverse solution result of wind direction and diffusion parameters with prune data

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

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

(2) 釋放參數的反解結果

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

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

表7 基于部分數據的擴散參數反解結果Table 7 Inverse solution result of release parameters with prune data

4 總結

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