遺傳算法動態(tài)修正核素大氣擴(kuò)散模型的適應(yīng)度函數(shù)研究

吉志龍 馬元巍 王德忠

（上海交通大學(xué) 核科學(xué)與工程學(xué)院 上海 200240）

核事故應(yīng)急支持系統(tǒng)中，核素擴(kuò)散模型為應(yīng)急處理提供技術(shù)支持，其準(zhǔn)確性影響到應(yīng)急決策是否及時有效。擴(kuò)散模型中，傳統(tǒng)的擴(kuò)散系數(shù)通常采用基于大量擴(kuò)散實(shí)驗(yàn)得到的經(jīng)驗(yàn)參數(shù)。由于事故狀態(tài)下的風(fēng)場、大氣湍流、地表等與實(shí)驗(yàn)條件存在差異，導(dǎo)致經(jīng)驗(yàn)擴(kuò)散參數(shù)難以準(zhǔn)確反映實(shí)際的擴(kuò)散過程。以IAEA和USNRC推薦的Pasquill-Gifford (PG)擴(kuò)散系數(shù)為例，其基于Prairie-Grass-Field (PGF)實(shí)驗(yàn)得到[1?2]，適用于多植被地形，摩擦層高度為3cm，近地釋放源的情形，通常用于保守預(yù)測，如應(yīng)用于核事故后果評價，其準(zhǔn)確性有待提高[3]。

例如在2002年福建惠安核電廠的示蹤實(shí)驗(yàn)中，根據(jù)PG系數(shù)得到的擴(kuò)散參數(shù)要明顯高于根據(jù)實(shí)際擴(kuò)散結(jié)果擬合得到的擴(kuò)散參數(shù)[4]。如果核電廠進(jìn)行了大氣擴(kuò)散實(shí)驗(yàn)則以在地的擴(kuò)散參數(shù)為準(zhǔn)，但是該擴(kuò)散參數(shù)忽略了大氣湍流強(qiáng)度以及溫度梯度變化等時序上的變化。

為了使經(jīng)驗(yàn)擴(kuò)散參數(shù)能更好地反映實(shí)際擴(kuò)散過程，一個可行的辦法是以經(jīng)驗(yàn)參數(shù)作為先驗(yàn)值，使用觀測數(shù)據(jù)對其進(jìn)行動態(tài)修正，使得擴(kuò)散模型計(jì)算值與觀測值盡可能相近，將經(jīng)過修正的擴(kuò)散參數(shù)代入原擴(kuò)散模型后，可以更準(zhǔn)確地預(yù)測下一階段核素的擴(kuò)散濃度分布，這類技術(shù)稱為源項(xiàng)反演或參數(shù)修正。Haupt等[5]首次提出使用遺傳算法耦合擴(kuò)散模型和受體模型，進(jìn)行源項(xiàng)反演，以達(dá)到增強(qiáng)污染物擴(kuò)散模擬準(zhǔn)確度的目的；Jeong等[6]使用類似的技術(shù)對核電廠泄露物的擴(kuò)散進(jìn)行了分析；Allen等[7]提出了動態(tài)修正風(fēng)場信息，以提高源項(xiàng)反演精確度的技術(shù)；Qin等[8]使用源項(xiàng)反演技術(shù)分析了英國Dundee地區(qū)的空氣污染源，對污染物擴(kuò)散過程進(jìn)行了模擬。研究表明，源項(xiàng)反演和參數(shù)動態(tài)修正是增強(qiáng)擴(kuò)散模型準(zhǔn)確性的有效方法，但對于適應(yīng)度函數(shù)和觀測誤差等影響算法性能的重要方面還缺少進(jìn)一步研究。本文將通過數(shù)值模擬的方法研究適應(yīng)度函數(shù)對參數(shù)修正的影響，找出更適合PG參數(shù)修正的適應(yīng)度函數(shù)。然后，使用最優(yōu)適應(yīng)度函數(shù)結(jié)合 Kincaid實(shí)驗(yàn)數(shù)據(jù)對擴(kuò)散模型進(jìn)行修正，研究修正后的擴(kuò)散模型預(yù)測能力的變化。

1 研究方法

1.1 拉格朗日煙團(tuán)模型

研究大氣擴(kuò)散模型的參數(shù)估計(jì)問題，需要擴(kuò)散模型具有能夠模擬變動流場，計(jì)算效率高等特點(diǎn)。相對于高斯煙羽模型和蒙特卡洛模型[9?10]，拉格朗日煙團(tuán)模型更符合此要求[11]，本文采用其作為核素擴(kuò)散模型。其基本思想為將連續(xù)釋放的氣態(tài)流出物簡化為連續(xù)釋放的煙團(tuán)，每個煙團(tuán)內(nèi)的污染物濃度分布符合高斯分布，環(huán)境內(nèi)某點(diǎn)的污染物濃度則為各個煙團(tuán)在該點(diǎn)濃度分布值的疊加。單個煙團(tuán)內(nèi)的濃度分布公式為：

式中，Q為該煙團(tuán)內(nèi)放射性物質(zhì)總活度，Bq·m?3；xc、yc、zc為煙團(tuán)中心的坐標(biāo)，m；zinv為逆溫層頂?shù)母?，m；x、y、z為空間任一點(diǎn)的坐標(biāo)，m；σx、σy、σz分別為煙團(tuán)的水平和垂直擴(kuò)散參數(shù)，m。存在多種根據(jù)大氣穩(wěn)定度和下風(fēng)向距離計(jì)算這些擴(kuò)散參數(shù)的經(jīng)驗(yàn)公式，其中最常用的是PG擴(kuò)散曲線[7]：

式中，x為煙團(tuán)中心距離釋放源的下風(fēng)向距離，m；py、qy、pz、qz為PG擴(kuò)散系數(shù)，由當(dāng)時的風(fēng)速、風(fēng)向、大氣穩(wěn)定度等信息決定。模型中假設(shè)兩個水平方向符合相同的濃度分布規(guī)律。

1.2 PG參數(shù)修正

擴(kuò)散參數(shù)σ是擴(kuò)散模型中的重要參數(shù)，而常用的PG擴(kuò)散曲線是基于Parairie Grass場地實(shí)驗(yàn)數(shù)據(jù)推導(dǎo)而來的，其使用條件存在一定局限性，當(dāng)實(shí)際情況與PG場地實(shí)驗(yàn)條件相差較大時，模型誤差較大[8]。為了使擴(kuò)散模型更好地反映實(shí)際擴(kuò)散情況，可以將原始的PG系數(shù)作為先驗(yàn)值，結(jié)合觀測數(shù)據(jù)，對原PG系數(shù)做修正：

1.3 使用遺傳算法修正PG系數(shù)

大氣擴(kuò)散模型具有非線性化程度高，參數(shù)多等特點(diǎn)，近年來，在基于大氣擴(kuò)散模型的源項(xiàng)反演研究中，使用遺傳算法作為最優(yōu)化過程求解的方法已經(jīng)取得了一些進(jìn)展[5?8,12]。

遺傳算法是模仿自然界生物進(jìn)化機(jī)制發(fā)展起來的隨機(jī)全局搜索和優(yōu)化方法，它借鑒了達(dá)爾文的進(jìn)化論和孟德爾的遺傳學(xué)說，由 Holland教授[13]于1975年首先提出。其本質(zhì)是一種高效、并行、全局搜索的方法。典型的遺傳算法執(zhí)行流程如圖1所示。

圖1 遺傳算法流程圖Fig.1 Flow chart of genetic algorithm.

圖1中遺傳算法參數(shù)是指種群規(guī)模、交叉率、突變率等，其影響了算法的資源消耗率和穩(wěn)定性。“選擇”、“交叉”、“變異”對應(yīng)于進(jìn)化論與遺傳學(xué)說中的自然選擇、交配、突變等概念。

2 適應(yīng)度函數(shù)選擇

適應(yīng)度函數(shù)對模型計(jì)算值和實(shí)際觀測值之間的差距給出度量標(biāo)準(zhǔn)，是遺傳算法修正擴(kuò)散系數(shù)的核心問題。

2.1 可選適應(yīng)度函數(shù)分析

作為模型計(jì)算結(jié)果與實(shí)際測量結(jié)果差距的度量，適應(yīng)度函數(shù)值越小則說明兩者越接近；另外，不同的觀測點(diǎn)，由于所處位置和探測條件等不同，所提供的觀測信息的可信度不同，所以適應(yīng)度函數(shù)應(yīng)符合以下兩個原則：

(1) 對于任一監(jiān)測站i，當(dāng)│Ci?Qi│增大，即觀測值與計(jì)算值差距增大時，適應(yīng)度函數(shù)值應(yīng)增大。

(2) 適應(yīng)度函數(shù)中，不同監(jiān)測站點(diǎn)的貢獻(xiàn)值所占權(quán)重可能不同，觀測數(shù)據(jù)越可信，所占權(quán)重越大。

對于原則(2)，如果引起監(jiān)測站探測信息置信度變化的僅為儀器誤差，且各站所用儀器的精度一致，則觀測值越小，相對誤差越大，探測信息的可信度越小。所以原則(2)可簡化為對探測值越小的站點(diǎn)，應(yīng)在適應(yīng)度函數(shù)中給予更小的權(quán)重，探測值越大，說明站點(diǎn)越接近擴(kuò)散流中心，更能反映擴(kuò)散的趨勢，應(yīng)在適應(yīng)度函數(shù)中給予更大的權(quán)重。極端情況下，對于探測數(shù)據(jù)為零的點(diǎn)，在適應(yīng)度函數(shù)中應(yīng)不予考慮，因?yàn)檫@樣的讀數(shù)跟離釋放源無窮遠(yuǎn)的點(diǎn)的讀數(shù)是一致的，該點(diǎn)讀數(shù)不能為修正PG系數(shù)提供有效信息。

基于以上兩點(diǎn)，本文在最小二乘法基礎(chǔ)上構(gòu)造并嘗試了使用以下4種遺傳算法的適應(yīng)度函數(shù)來對PG系數(shù)進(jìn)行修正：

式中，N表示觀測點(diǎn)總數(shù)；i表示第i個觀測點(diǎn)；Ci為擴(kuò)散模型在第i點(diǎn)的濃度計(jì)算值；Oi為第i點(diǎn)的實(shí)際測量值；di為觀測點(diǎn)i到釋放源的距離。ft2中α為常數(shù)，用于保證αCi和αOi均大于1。

在這4種適應(yīng)度函數(shù)中，ft1方法簡單，對各站點(diǎn)信息進(jìn)行平權(quán)處理，沒有考慮觀測點(diǎn)誤差的影響。ft2由于對濃度進(jìn)行了對數(shù)運(yùn)算，所以ft2對濃度的指數(shù)級變化更敏感。ft3和ft4是對ft1的修正，通常探測濃度越小，或者離釋放源越遠(yuǎn)，則意味著該觀測點(diǎn)離釋放流的中心距離更遠(yuǎn)，對于校正模型所提供的信息也越小，故在適應(yīng)度函數(shù)中給予更小的比重。

2.2 數(shù)值模擬

本文將采用數(shù)值模擬的方法比較4種使用適應(yīng)度函數(shù)對PG系數(shù)進(jìn)行修正的效果，通過設(shè)置擴(kuò)散模型參數(shù)，使用擴(kuò)散模型生成濃度場，并在該濃度場上加入一定噪聲以模擬真實(shí)測量值，使用“測量值”反演設(shè)定的模型參數(shù)。

2.2.1 擴(kuò)散模型參數(shù)設(shè)置

為簡化計(jì)算，本文將釋放源簡化為位于坐標(biāo)(0,0)的單點(diǎn)源，其距離地表高度為187 m，以9.4 g·s?1的速度釋放示蹤氣體，風(fēng)場穩(wěn)定，風(fēng)向?yàn)?00°，風(fēng)速為 3 km·h?1，PG 系數(shù)設(shè)置為 PGorig：

使用拉格朗日擴(kuò)散模型計(jì)算示蹤氣體持續(xù)釋放1 h后的地表濃度分布，結(jié)果如圖2所示。

圖2 擴(kuò)散1 h后的地表濃度分布Fig.2 Concentrationdistributionsimulated after one hour’s dispersion.

2.2.2 觀測點(diǎn)設(shè)定

由圖2，1 h后由于風(fēng)向原因，擴(kuò)散區(qū)域主要集中在釋放源的東南方向。從該區(qū)域均勻抽取24個點(diǎn)作為“觀測點(diǎn)”，對這些點(diǎn)上的濃度加入噪聲，得到信噪比(signal-to-noise ratio, SNR)分別為無限大（無噪聲）、50、5、2.5、0.25 的“觀測數(shù)據(jù)”，分別對應(yīng)無噪聲、小噪聲、噪聲與信號數(shù)量級一致、噪聲大于信號情形。信噪比越小，說明監(jiān)測站的觀測誤差越大。SNR無限大表示沒有加入噪聲。

2.2.3 PG系數(shù)修正

針對以上得到的“觀測值”，使用遺傳算法配合不同的適應(yīng)度函數(shù)對PG參數(shù)進(jìn)行修正，為模擬真實(shí)情況，對PGorig加入擾動，作為PG參數(shù)的先驗(yàn)值PGpre：

針對PGpre，使用遺傳算法計(jì)算得到使得各適應(yīng)度函數(shù)值最小的最優(yōu)PG修正系數(shù)，將這些修正系數(shù)乘以原PG系數(shù)后，得到修正后的PG系數(shù)PGcorr。比較修正結(jié)果 PGcorr與本文設(shè)定的 PGorig之間的差異，差異越小，則說明修正效果越好。

2.3 結(jié)果分析

通過數(shù)值模擬，得到對應(yīng)于4種適應(yīng)度函數(shù)的4組PQ擴(kuò)散參數(shù)，各自的修正誤差由式(10)得到：

由于q在PG參數(shù)公式中處于下風(fēng)向距離x的冪的位置，對擴(kuò)散參數(shù)生成的影響更大，故在式(10)中對p進(jìn)行取對數(shù)處理。理論上，系數(shù)修正效果越好，則式(10)的值越小。

4組數(shù)值模擬的結(jié)果圖3所示。由圖3，當(dāng)觀測誤差很小時，4種適應(yīng)度函數(shù)修正 PG參數(shù)的效果相近，隨著觀測誤差增大，ft3和ft4要明顯好于前兩者，使用ft2的修正效果最差。

圖3 不同適應(yīng)度函數(shù)修正效果隨觀測誤差的變化Fig.3 Correction result of different fitness functions.

結(jié)果表明：

(1) 當(dāng)某個觀測點(diǎn)的計(jì)算濃度值與觀測值差距增大時，4個適應(yīng)度函數(shù)的值都增大，所以當(dāng)觀測條件為理想狀態(tài)，即觀測數(shù)據(jù)都絕對準(zhǔn)確時，它們的修正效果基本一致。

(2)ft1對各監(jiān)測站數(shù)據(jù)作平權(quán)處理，即假設(shè)各站監(jiān)測數(shù)據(jù)的可信度是相等的。這只有在測量儀器和測量條件在理想狀況下時才能滿足。

(3)ft2的關(guān)鍵部分是 log10(aCi+1.0)?log10(aOi+1.0)=log10[(aCi+1.0)/(aOi+1.0)，模型計(jì)算誤差一定時，Oi越小，適應(yīng)度函數(shù)越大，這違反了原則二，所以當(dāng)觀測誤差較大時，ft2的修正效果最差。

(4)ft3和ft4是對ft1的修正，ft3中使用探測值大小做加權(quán)，是原則(2)的直觀體現(xiàn)。ft4中以觀測站到釋放源的距離倒數(shù)做加權(quán)，也體現(xiàn)了原則(2)的思想，因?yàn)橐话汶x釋放源越遠(yuǎn)，則擴(kuò)散濃度值越小。所以當(dāng)觀測誤差增大時，這兩個適應(yīng)度函數(shù)的表現(xiàn)要好于ft1。

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

從數(shù)值模擬的結(jié)果中可知，使用遺傳算法和考慮了觀測誤差的適應(yīng)度函數(shù)，可以有效地對PG系數(shù)進(jìn)行修正，使其更符合真實(shí)的擴(kuò)散情況。這有助于使用擴(kuò)散模型預(yù)測氣載核素的下一步擴(kuò)散過程，為核應(yīng)急預(yù)警提供信息。本文將使用 Kincaid實(shí)驗(yàn)數(shù)據(jù)來驗(yàn)證PG系數(shù)修正對拉格朗日煙團(tuán)模型預(yù)測能力的影響。

3.1 Kincaid實(shí)驗(yàn)數(shù)據(jù)

Kincaid實(shí)驗(yàn)是于 1980?1981年在美國伊利諾斯州 Kincaid電廠開展的一項(xiàng)大氣擴(kuò)散實(shí)驗(yàn)。實(shí)驗(yàn)以六氟化硫?yàn)槭聚櫄怏w，共進(jìn)行了23組釋放實(shí)驗(yàn)，每次釋放實(shí)驗(yàn)持續(xù) 3?9 h，分布在距離釋放點(diǎn)不同距離段上的監(jiān)測站在實(shí)驗(yàn)中實(shí)時記錄探測到的六氟化硫濃度[14]。

3.2 驗(yàn)證方法

以1981年5月28日的Kincaid實(shí)驗(yàn)結(jié)果為例，這組實(shí)驗(yàn)共持續(xù)了9 h，一共有101個觀測站，每個觀測站都記錄了每個小時末探測到的六氟化硫濃度，同時EPRI 空氣質(zhì)量數(shù)據(jù)中心提供了實(shí)驗(yàn)時詳盡的氣象數(shù)據(jù)，本文通過如下方法來驗(yàn)證PG參數(shù)校正對模型預(yù)測能力的影響：

(1) 對1?8 h內(nèi)每階段得到的監(jiān)測數(shù)據(jù)，使用標(biāo)準(zhǔn)PG參數(shù)作為先驗(yàn)值，ft3或ft4作為適應(yīng)度函數(shù)，經(jīng)遺傳算法對各階段的 PG參數(shù)進(jìn)行修正，得到 8組修正后的PG參數(shù)值。

(2) 將前面得到的8組PG參數(shù)修正值做平均，代入拉格朗日擴(kuò)散模型，計(jì)算第9個小時各觀測站的濃度值。

(3) 將標(biāo)準(zhǔn)PG參數(shù)代入拉格朗日擴(kuò)散模型，計(jì)算第9個小時各觀測站的濃度值。

(4) 將第(2)、(3)步中計(jì)算得到的第9個小時后各觀測站的濃度值與 Kincaid實(shí)驗(yàn)記錄值做比較。觀察使用修正后的PG參數(shù)和未經(jīng)修正的PG參數(shù)的擴(kuò)散模型計(jì)算值與真實(shí)測量值之間的差異。這里使用統(tǒng)計(jì)學(xué)量 FA2、BIAS、FB等表征計(jì)算值與觀測值的相近程度。

3.3 結(jié)果分析

通過對Kincaid實(shí)驗(yàn)中的20組實(shí)驗(yàn)數(shù)據(jù)作上述驗(yàn)證，將每組實(shí)驗(yàn)的統(tǒng)計(jì)結(jié)果作平均后得到的結(jié)果如表1所示。

表1 PG參數(shù)修正對模型預(yù)測值的影響Table 1 Influence of PG correction on dispersion model’s forecast ability.

FA2的統(tǒng)計(jì)意義為 0.5≤C/O≤2，即擴(kuò)散模型預(yù)測值與觀測值相差在兩倍以內(nèi)的站點(diǎn)在總站點(diǎn)數(shù)中的比例，F(xiàn)A5、FA10類似，這些統(tǒng)計(jì)數(shù)字越大，則說明模型預(yù)測值與觀測值更接近；BIAS為平均偏差，其定義為即模型計(jì)算平均值與真實(shí)測量平均值之差，該值絕對值越小，則說明預(yù)測更準(zhǔn)確；是平化偏差(Flatted Bias)，其絕對值越小，則預(yù)測更準(zhǔn)。從表1中可以看出，與直接使用PG擴(kuò)散參數(shù)相比，在拉格朗日擴(kuò)散模型中使用經(jīng)遺傳算法修正后的PG參數(shù)，模型的計(jì)算結(jié)果與實(shí)測結(jié)果更相符，模型預(yù)測的準(zhǔn)確性得到了提高。

4 結(jié)語

本文提出了使用遺傳算法，利用實(shí)時觀測數(shù)據(jù)對放射性核素大氣擴(kuò)散模型中的經(jīng)驗(yàn)擴(kuò)散參數(shù)進(jìn)行實(shí)時修正，以增強(qiáng)擴(kuò)散模型準(zhǔn)確性的方法。通過數(shù)值模擬，對比分析了使用遺傳算法對拉格朗日大氣擴(kuò)散模型中的 Pasquill-Gifford擴(kuò)散系數(shù)進(jìn)行修正時，4種適應(yīng)度函數(shù)的修正效果。結(jié)果表明，使用考慮了觀測誤差和觀測站距離等影響的適應(yīng)度函數(shù)

ft3和ft4對PG系數(shù)修正的結(jié)果更接近于真實(shí)值。將在數(shù)值模擬中表現(xiàn)占優(yōu)的ft3和ft4應(yīng)用到Kincaid實(shí)驗(yàn)數(shù)據(jù)集發(fā)現(xiàn)，使用遺傳算法對PG系數(shù)進(jìn)行修正后，擴(kuò)散模型的預(yù)測能力得到了提高。

1 IAEA. Atmospheric dispersion model for application in relation to radionuclide release[M]. Austria: International Atomic Energy Agency, 1986

2 USNRC. Methods for estimating atmospheric transport and dispersion of gaseous effluents in routine release from light water cooled reactor[M]. Washington: Nuclear Regulatory Commision, 1977

3 Venkatram Akula. An examination of the Pasquill-Giffrd-Turner dispersion scheme[J]. Atmospheric Environment, 1996, 30(8): 1283?1290

4 胡二邦, 幸存田, 宣仁義, 等. 福建惠安核電廠址 SF6示蹤實(shí)驗(yàn)研究[J]. 環(huán)境科學(xué)學(xué)報, 2004, 24(2): 320?325 HU Erbang, XING Cuntian, XUAN Renyi,et al. The SF6 tracer study in the Fujian Huian nuclear plant[J]. Acta Scientiae Circumstantiae, 2004, 24(2): 320?325

5 Haupt S E, Young G S, Allen G T. Validation of a receptor-dispersion model coupled with a genetic algorithm using synthetic data[J]. Journal of Applied Meteorology and Climatology, 2006, 45(3): 476?490

6 Jeong H J, Kim E H, Suh K S,et al. Determination of the source rate released into the environment from a nuclear power plant[J]. Radiation Protection Dosimetry, 2005,113(3): 308?313

7 Allen C T, Young G S, Haupt S E. Improving pollutant source characterization by better estimating wind direction with a genetic algorithm[J]. Atmospheric Environment,2007, 41(11): 2283?2289

8 Qin Y, Oduyemi K. Atmospheric aerosol source identification and estimates of source contributions to air pollutions in Dundee, UK[J]. Atmospheric Environment,2003, 37(13): 1799?1809

9 閆政, 吳信民, 鄧?yán)? 等. 使用蒙特卡羅方法模擬核事故氣載放射性污染物大氣擴(kuò)散[J]. 核技術(shù), 2011, 34(3):193?198 YAN Zheng, WU Xinmin, DENG Lei,et al. An M-C model to simulate airborne radioactive material dispersion in nuclear accident[J]. Nuclear Techniques, 2011, 34(3):193?198

10 Turner D B. Workbook of atmospheric dispersion estimates: an introduction to dispersion modelling[M].Florida: Lewis Publishers, Boca Raton, 1994

11 Raza S, Avila R, Cervantes J. A 3D Lagrangian particle model for the atmospheric dispersion of toxic pollutants[J].International Journal of Energy Research, 2002, 26(2):94?95

12 寧莎莎, 蒯琳萍. 混合遺傳算法在核事故源項(xiàng)反演中的應(yīng)用[J]. 原子能科學(xué)技術(shù), 2012, 46(z1): 469?472 NING Shasha, KUAI Linping. Back-calculation of source terms by hybrid genetic algorithm in nuclear power plant accident[J]. Atomic Energy Science and Technology, 2012,46(zl): 469?472

13 Haupt R L, Sue Ellen Haupt. Practical genetic algorithms[M]. New Jersey: John Wiley & Sons, Inc, 2004:23?24

14 Olesen H R. User’s guide to the model validation kit[M].Denmark: National Research Institute, 2005: 18?26