拉格朗日煙團模型的大氣擴散系數自適應修正

談文姬，王德忠，馬元巍，吉志龍

（上海交通大學 機械與動力工程學院，上海 200240）

自1986年發生切爾諾貝利事故以來，國際上對核事故后果評價和決策支持系統的研究有了很大進展，特別是在2011 年福島核事故之后，放射性核素釋放到環境中的影響引起了國際社會的廣泛關注。核事故后果決策支持系統能對事故發生后放射性核素泄漏到空氣中的濃度分布進行實時和預測性計算，為保護周圍公眾的安全提供數據依據。

目前，用于核事故應急決策支持系統中的大氣擴散模型主要有3類：高斯煙羽模型、拉格朗日煙團模型和蒙特卡羅粒子模型［1］。高斯煙羽模型適用于穩定均勻流動及平坦地形的條件，雖模型簡單、運算快捷，但難以適用于復雜流場下的實際核事故。蒙特卡羅粒子模型能較為真實地模擬粒子在大氣中的輸運過程，但計算時間較長，并存在邊界稀疏的問題，往往不能滿足實時核事故應急響應［2］。煙團模型能適用于非穩定流動及復雜地形的條件，在局地尺度上有較好的精度，應用較為廣泛，如歐洲RODOS等核事故應急決策支持系統均采用該模型［3］。

描述煙團擴散能力的大氣擴散系數，是影響拉格朗日煙團模型的重要參數。傳統的拉格朗日煙團模型采用的擴散參數是由Pasquill-Gifford（P-G）曲線確定的。P-G 曲線是根據理想條件下的實驗結果擬合得到的，與實際情況的不符將給煙團模型帶來偏差。本文針對拉格朗日煙團模型不能實時反映局地擴散能力的缺點，提出大氣擴散系數自適應的改進方法，并使用實際大氣擴散示蹤實驗的數據，對大氣擴散系數自適應修正方法進行檢驗。

1 煙團模型的大氣擴散系數自適應修正

1.1 煙團模型簡介

煙團模型是將釋放的污染物假想成離散的煙團，對每個煙團中心的輸運過程進行模擬。每個煙團在某一特定時間間隔內固定不動，根據此刻固定住的煙團計算濃度。然后煙團在下一個時間步長內繼續移動，大小和強度繼續變化，直到下次采樣時間再次固定。在基本時間步長內，接受點的濃度為周圍所有煙團采樣時間內的平均濃度總和。有風的條件下，第i個煙團在某一點的濃度［4］為：

式中：Ci（x，y，z）為第i個煙團（x，y，z）點的濃度，g／m3；Q（i）為第i個煙團的強度，g／s；σx、σy、σz分別為沿風向、水平方向、垂直方向的擴散參數，m；（xc（i），yc（i），zc（i））為第i個煙團的中心位置，m。

計算污染物濃度及分布時，必須知道源強、風場、有效源高和大氣擴散參數。其中，源強往往是通過測量獲得或由設計基準事故給出，氣象信息由當地氣象測量站獲取。大氣擴散參數對擴散后果有較大的影響，但目前仍為理想條件下的示蹤實驗數據計算得到，這些數據在核事故后果評價中與實際情況有較大偏差，應提出新的方法對其進行修正。

1.2 大氣擴散系數自適應修正方法

大氣擴散參數（σx、σy、σz）表示污染物在隨機湍流場中的擴散能力和散布范圍，它與下風向距離x、大氣穩定度和下墊面粗糙度等相關。大氣擴散參數隨x 及下墊面粗糙度的增大而增大，隨大氣穩定度的增大而減小。

大多數模型均通過Pasquill氣象穩定度與P-G 曲線來確定大氣擴散參數，該曲線是由大量擴散實驗和理論分析得到的，具體確定擴散參數的方法為：首先根據云況和日射及地面風速選擇大氣穩定度，然后根據擴散曲線讀出不同下風向距離處的擴散參數。RIMPUFF模型中使用的Karlsruhe-Jülich（K-J）大氣擴散系數即采用此方法，K-J大氣擴散系數是根據大氣穩定度及源項的釋放高度將大氣擴散系數進行分類，共分為18套大氣擴散系數［5］。但由于P-G曲線是依據平坦理想條件下的實驗得到的，對于不同源高、下墊面等條件下的擴散曲線還有待進一步修正。

本文提出一種根據當地的實測數據（源項釋放數據、觀測點位置和測得的濃度數據、氣象數據）及煙團模型的濃度計算公式，對大氣擴散系數進行實時估計。這種將上一步大氣擴散系數的估計結果作為修正的大氣擴散系數進行下一個時間步長運算的方法，稱為大氣擴散系數的自適應修正。將得到的適用于當地的實時大氣擴散系數作為修正后的系數繼續進行擴散計算，較使用理想條件下的實驗得到的擴散參數更符合實際情況，能提高拉格朗日煙團模型的準確性。

假設擴散參數為下風向距離x 的冪函數，即：

式中：py、qy、pz、qz為大氣擴散系數。

若確定py、qy、pz、qz，就 可 得 到 適 合 當 時當地地形氣象條件的σy和σz。這4個大氣擴散系數可利用最小二乘法確定。為了權衡觀測數據的精度，引入一標志測量精度的權重g，使地面濃度的計算值Cmi（Cmi=C（xi，yi，0））與實際觀測值COi之差的平方和S 最小，則S 可表示為：

式中：n為實際觀測的數據個數；gi為采樣點的權重，gi=COi／CO，max，CO，max為觀測點采集數 據中的最大濃度。

采用梯度迭代法計算求解使S 為最小的非線性方程組，可得到實時的大氣擴散系數。本文根據上述方法，針對拉格朗日煙團模型進行修正，開發了相應的大氣擴散系數自適應模型。

2 模型驗證

2.1 Kincaid擴散實驗

MVK 是荷蘭國家公共健康與環境研究院開發的大氣擴散模型驗證對比工具［6］，包括Kincaid、Indianapolis、Copenhagen、Lillestrom 4個實驗數據集，定量的統計分析工具（BOOT程序）和定性的圖形分析工具（SIGPLOT 程序），以及可用于評價模型的準確性、適用性和揭示模型的錯誤和數據的壞點。目前，MVK 已廣泛用于驗證和比較自主開發的大氣擴散模型［7］。

Kincaid實驗是EPRI煙羽模型驗證和改進項目的一部分。EPRI空氣質量數據中心提供了實驗時的氣象數據，包括氣壓、溫度、濕度、風速和風向。Kincaid電站坐落在美國伊利諾斯州（北緯39.59°，西經89.49°），平均海拔約為180m，其周圍是平坦的農田及湖泊。圖1為以Kincaid電站為中心的80km×80km 范圍的等高地形圖。圖1中為UTM 坐標，x正方向為東，y正方向為北，網格劃分為100m×100m。

圖1 以Kincaid電站為中心的80km×80km 范圍的等高地形圖Fig.1 80km×80km wide topographic contour centered at Kincaid Power Plant

實驗時，六氟化硫將從電站中的一座高187m、直徑9 m 的煙囪頂部釋放。與釋放點距離為1、2、3、5、7和10km 的圓弧上設有觀測點，采樣并記錄數據。實驗數據集包含約171h的數據。

2.2 數據驗證方法

BOOT程序是MVK中的統計工具，能將模型預測值與實驗觀測值進行比較，且可計算出各種統計參數。除了模型預測與實驗觀測的平均值MEAN、標準差SIGMA 外，其他數據還包括平均偏差BIAS、標準化均方差NMSE、相關系數CORR、比例偏差FB、FA2等，這些統計參數能為模型評價提供依據［8］，其計算公式為：

式中：CM為模型預測濃度，μg／m3；CO為觀測濃度，μg／m3；σM為模型預測濃度的標準差，μg／m3；σO為觀測濃度的標準差，μg／m3。FA2為0.5≤CM／CO≤2的幾率，即模型計算值誤差在50%以內的幾率。

3 計算結果與討論

3.1 大氣擴散系數最優擬合方法

擬合大氣擴散系數所需數據包括：觀測站數據（觀測站相對于源的位置坐標和測得的濃度數據）、氣象數據（當時的風速、風向，由氣象站提供）、源項釋放率（實驗設計提供）。1980年4月20日14：00的測量數據列于表1。其中源項釋放處的高度為源高處，源高處風向以正北方向為0°。

表1 Kincaid實驗1980年4月20日測量數據Table 1 Measurement data of 1980-04-20in Kincaid experiment

大氣擴散系數初值的設置會對計算得到的大氣擴散系數帶來較大影響。不同的初值求解得到的大氣擴散系數不同，需對解得的大氣擴散系數進行篩選。本文提出一種大氣擴散系數最優擬合方法，根據表1中的數據，當大氣擴散系數py、qy、pz、qz初值為（1，1，1，1）時，求解得到大氣擴散系數A=（0.739，0.881，0.838，0.959）；當初值為（0.8，0.8，0.8，0.8）時，解得大氣擴散系數B=（3.84，0.681，0.553，1.007）。將A 和B代入拉格朗日煙團模型時，兩者計算得到的觀測點坐標處的濃度相同，無法立即判斷這兩組大氣擴散系數的適用性。在這種情況下，將A、B 兩組大氣擴散系數作為初值代入下一個時間步長的迭代計算。下一步計算結果表明，使用A 計算的濃度較使用B 計算的濃度更接近觀測值，則選擇A 代入下一步的計算。重復該過程即可得到大氣擴散系數的最優解。

3.2 示蹤實驗驗證結果

用MVK 中的Kincaid數據集對大氣擴散系數自適應修正方法進行檢驗，并對大氣擴散系數自適應修正方法的煙團模型與K-J大氣擴散系數煙團模型的計算結果進行比較，以判斷該方法的可行性。根據Kincaid實驗記錄，從K-J大氣擴散系數中選擇與源項釋放高度為180m 和D 類大氣穩定度相對應的一組大氣擴散系數（0.208，0.903，0.307，0.734）［9］。

對實驗值（OBS）、K-J大氣擴散系數煙團模型的計算值和大氣擴散系數自適應修正的煙團模型（APM）的計算值進行對比的結果列于表2～6。

表2 1980年4月20日3h18個有效觀測站數據檢驗結果Table 2 Result of comparison of 18valid experiment data in 3hon 1980-04-20

表3 1980年4月25日6h45個有效觀測站數據檢驗結果Table 3 Result of comparison of 45valid experiment data in 6hon 1980-04-25

表4 1980年5月1日6h20個有效觀測站數據檢驗結果Table 4 Result of comparison of 20valid experiment data in 6hon 1980-05-01

表5 1980年5月4日6h19個有效觀測站數據檢驗結果Table 5 Result of comparison of 19valid experiment data in 6hon 1980-05-04

由表2、3可看出，與大氣擴散系數煙團模型相比，大氣擴散系數自適應修正方法計算得到的平均值、標準差更接近于觀測值，比例偏差FB更接近于0、FA2更接近于1。采取大氣擴散系數自適應方法的煙團模型的各項統計參數均優于KJ大氣擴散系數煙團模型的。就模型計算值與觀測值誤差在50%以內的幾率而言，準確性提高約20%。表2、3、6的數據修正效果優于表4、5的，這是由于風向變化較大對修正效果產生了負影響。1980年4月20日、4月25日及5月7日3天的風向變化較?。?°～40°），而5月1日、4日兩天的風向變化較大，有時達180°。

表6 1980年5月7日9h36個有效觀測站數據檢驗結果Table 6 Result of comparison of 36valid experiment data in 9hon 1980-05-07

圖2、3分別為大氣擴散系數自適應修正方法模型及K-J大氣擴散系數模型的預測值與觀測值的分位數圖（Q-Q 圖）。

圖2 APM 模型的Q-Q 圖Fig.2 Q-Q plot of APM model

圖3 K-J模型的Q-Q 圖Fig.3 Q-Q plot of K-J model

從圖2和3可看出，兩個模型的預測值較實驗值均偏大，這種保守的估計，在核事故濃度預測上是可被接受的。APM 模型點的個數多于K-J模型的，這是由于APM 模型預測的觀測點處濃度的非零值較多。修正模型的預測值的Q-Q 圖更接近y=x 直線，表明APM 模型的預測值較K-J模型的更符合觀測值。

3.3 討論

由于傳統煙團模型中的擴散參數是根據大氣穩定度及釋放高度的條件選取的，可供選取的值是以理想條件下的實驗結果為依據，因此，擴散參數在復雜流動和地形條件下，適用性受到了很大局限。為改善擴散參數的適用性，利用當時當地的氣象數據、源項數據及觀測數據，以最小二乘法求解大氣擴散系數，將符合實際情況的大氣擴散系數代入模型中，以達到系數自適應修正的效果。通過計算，將傳統擴散參數的煙團模型與系數自適應修正的煙團模型的計算結果進行比較，后者與觀測數據一致性的統計結果得到明顯的提高，證明大氣擴散系數自適應修正方法是有效的。

在擬合求解大氣擴散系數時，一個時間步長中風速風向的影響會累積到下一個計算時間步長內?？s短計算時間步長可減少風速風向變化帶來的影響，更準確地求解大氣擴散系數。但在計算過程中要求有與計算時間步長相匹配的氣象數據。

由計算過程可發現，大氣擴散系數自適應修正方法需要有較為準確的觀測值和釋放數據。由于對觀測值的數量和精度均有一定的要求，因此，在廠區內如何合理有效地布置監測站是需要考慮的問題。事故發生時，破口出現的位置及泄漏的數據應盡可能地詳細，因此安全殼內各部件的監測系統需靈敏、準確，且需發展源項反演技術，以便獲得準確的釋放數據。

4 結論

本文針對拉格朗日煙團模型提出了一種大氣擴散系數自適應的方法，建立了動態的大氣擴散系數自適應修正的拉格朗日煙團模型。通過Kincaid實驗數據驗證表明，大氣擴散系數自適應修正方法的煙團模型優于傳統P-G 曲線擴散參數的煙團模型?；谛路椒ń⒌膭討B大氣擴散系數自適應修正煙團模型可提高大氣擴散模型計算結果的準確性。

［1］ 胡二邦，姚仁太，宣儀仁，等.核事故后果評價系統的進展與比較［J］.輻射防護通訊，2003，23（1）：6-12.HU Erbang，YAO Rentai，XUAN Yiren，et al.Comparison and advance of consequences assessment system for nuclear accident［J］.Radiation Protection Bulletin，2003，23（1）：6-12（in Chinese）.

［2］ 曲靜原，曹建主，劉磊，等.我國核應急決策支持系統研究開發的現狀與展望［J］.原子能科學技術，2001，35（3）：283-288.QU Jingyuan，CAO Jianzhu，LIU Lei，et al.Status and perspective on the research and development of the Chinese decision support system for nuclear emergency management［J］.Atomic Energy Science and Technology，2001，35（3）：283-288（in Chinese）.

［3］ EHRHARDT J，WEIS A.Development of a comprehensive decision support system for nuclear emergency in europe following an accidental release to atmosphere，FZKA-5772［R］.Karlsruhe：Forshungszentrum Karlsruhe，1996.

［4］ 蔣維楣.空氣污染氣象學［M］.南京：南京大學出版社，2003：102-111.

［5］ 曲靜原，曹建主.RODOS4.0與RODOS的未來發展［J］.輻射防護，2002，22（4）：198-206.QU Jingyuan，CAO Jianzhu.The future development of RODOS4.0and RODOS［J］.Radiation Protection，2002，22（4）：198-206（in Chinese）.

［6］ OLESEN H R.Model validation kit-status and outlook［J］.International Journal of Environment and Pollution，2000，14（12）：65-76.

［7］ OLESEN H R.User’s guide to the model validation kit［R］.Denmark：National Environment Research Institute，2005.

［8］ ELEVELD H.Application of a methodology to validate atmospheric dispersion models［C］∥Proceedings of the 7th International Conference on Harmonization Within Atmospheric Dispersion Modelling for Regulatory Purposes.Italy：［s.n.］，2001.

［9］ THYKIER-NIELSEN S，DEME S，MIKKELSEN T.Description of the atmospheric dispersion module RIMPUFF［R］.Denmark：Risoe National Laboratory，1999.