吉志龍 馬元巍 王德忠

（上海交通大學 核科學與工程學院 上海 200240）

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

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

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

1 研究方法

1.1 拉格朗日煙團模型

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

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

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

1.2 PG參數(shù)修正

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

2.2 數(shù)值模擬

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

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

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

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

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

2.2.2 觀測點設(shè)定

由圖2，1 h后由于風向原因，擴散區(qū)域主要集中在釋放源的東南方向。從該區(qū)域均勻抽取24個點作為“觀測點”，對這些點上的濃度加入噪聲，得到信噪比(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ù)進行修正，為模擬真實情況，對PGorig加入擾動，作為PG參數(shù)的先驗值PGpre：

針對PGpre，使用遺傳算法計算得到使得各適應(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擴散參數(shù)，各自的修正誤差由式(10)得到：

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

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

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

結(jié)果表明：

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

(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)，模型計算誤差一定時，Oi越小，適應(yīng)度函數(shù)越大，這違反了原則二，所以當觀測誤差較大時，ft2的修正效果最差。

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

3 示蹤實驗驗證

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

3.1 Kincaid實驗數(shù)據(jù)

Kincaid實驗是于 1980?1981年在美國伊利諾斯州 Kincaid電廠開展的一項大氣擴散實驗。實驗以六氟化硫為示蹤氣體，共進行了23組釋放實驗，每次釋放實驗持續(xù) 3?9 h，分布在距離釋放點不同距離段上的監(jiān)測站在實驗中實時記錄探測到的六氟化硫濃度[14]。

3.2 驗證方法

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

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

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

(3) 將標準PG參數(shù)代入拉格朗日擴散模型，計算第9個小時各觀測站的濃度值。

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

3.3 結(jié)果分析

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

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

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

4 結(jié)語

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

ft3和ft4對PG系數(shù)修正的結(jié)果更接近于真實值。將在數(shù)值模擬中表現(xiàn)占優(yōu)的ft3和ft4應(yīng)用到Kincaid實驗數(shù)據(jù)集發(fā)現(xiàn)，使用遺傳算法對PG系數(shù)進行修正后，擴散模型的預(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示蹤實驗研究[J]. 環(huán)境科學學報, 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 閆政, 吳信民, 鄧磊, 等. 使用蒙特卡羅方法模擬核事故氣載放射性污染物大氣擴散[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 寧莎莎, 蒯琳萍. 混合遺傳算法在核事故源項反演中的應(yīng)用[J]. 原子能科學技術(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