本發明屬于核與輻射突發事件后果大尺度模擬領域，具體涉及一種氣載放射性核素長距離遷移拉格朗日粒子擴散計算方法，用于對核設施在發生核事故時釋放或可能釋放到環境中的物質所造成的后果進行預測或評價。
背景技術：
一直以來，世界各國均非常重視和關注放射性核素大尺度大氣遷移問題，不僅是軍事、國防方面的需要，而且是核事故應急、環境危害問題等的需要。放射性物質在大氣中的長距離遷移和后果評價與分析技術涉及眾多的科學領域，屬于多學科交叉的技術，模擬技術和評價方法也在不斷地更新與發展，因此，包括美國、歐盟、日本在內的發達國家、地區也在不斷加強該技術的研發。另外，放射性核素大尺度遷移數值模擬技術已經作為一種國家的技術資源，隨著相關領域技術的不斷發展與更新，從戰略角度考慮，一個國家必須有技術儲備和可持續發展的能力。當前國際形勢變化多端，需要我們不斷建立與完善大氣輸運模擬技術，并提供技術支撐。
近些年來，在朝核試驗、日本福島核電泄漏等事件接連發生的情況下，建立并不斷完善一套可應對我國周邊數百至數千公里范圍內氣載放射性物質釋放情形應急響應系統的必要性已經成為共識。該系統的建立將為決策人員能夠面對類似事件，快速、科學、有效地提出決策建議與方案，最終達到保護我國公眾與環境安全的目標。
為此，中國輻射防護研究院自主開發、建立了核與輻射突發事件后果大尺度模擬系統，包括放射性核素大尺度大氣遷移數值模擬技術和latmes1.0計算程序，以及相關的模擬結果統計分析程序，可應對全球范圍的核與輻射突發事件下，氣載放射性物質的遷移、擴散及其環境安全后果的評價。針對朝鮮核試驗，根據總裝的要求，中輻院定期提供階段報告預測和分析了朝鮮核試驗場址以及寧邊核場址放射性核素釋放大氣遷移對我國的影響。目前已成功用于對日本福島核事故、朝鮮第三次核試驗等事件的有效響應。
境外核爆輻射后果評價系統是通過將區域氣象數值預報產品和中長距離放射性污染物遷移模擬和后果評價集成為一體。該系統設計為一個連接多種模型、數據庫和地理信息的綜合系統，用以模擬、預測和評估指定區域范圍內突發事件產生的影響。
技術實現要素：
針對現有技術中存在的問題，本發明提供一種氣載放射性核素長距離遷移軌跡拉格朗日粒子擴散計算方法，根據該方法能夠估算所有粒子在時間和空間上的總體分布，便于在發生突發事件時，決策人員能夠面對該類事件，快速、科學、有效的提出決策建議與對策，從而保護公眾與環境的安全。
為達到以上目的，本發明采用的技術方案是：提供一種氣載放射性核素長距離遷移拉格朗日粒子擴散計算方法，包括如下步驟：
1)獲取天氣預報數據、地形數據及模式參數；天氣預報數據兩個水平風分量u和v、一個垂直風分量w、溫度t、濕度的三維場；模式參數包括釋放源項參數和設定區域范圍；
2)空間坐標和時間坐標變換；
3)調用風場；
4)通過積分粒子運動的拉格朗日方程，計算粒子運行軌跡；xi(t+δt)＝xi(t)+vi(xi(t),t)δt+vi′(xi(t),t)δt+vi″(xi(t),t)δti＝1,2,3
其中，xi為粒子的三維坐標分量；vi為平均風速分量
5)計算放射性衰變和干、濕沉積；
6)統計粒子密度、計算網格濃度；
7)考慮粒子分裂條件，跟蹤所有粒子，檢查是否有新粒子產生，如果有新粒子產生，返回步驟3)，繼續計算；如果沒有新粒子產生，結束計算。
進一步，步驟2)中，所述空間坐標和時間坐標變換是將等壓層氣象數據轉換為隨地笛卡爾坐標。
進一步，步驟4)中，對運行軌跡中湍流參數和中尺度風脈動參數確定；湍流參數，采用hanna提出的參數化方法，即根據邊界層參數混合層高度h、莫寧-奧布霍夫長度l、對流速度尺度w*、粗糙長度z0和摩擦速度u*來計算；中尺度風脈動參數，根據maryon提出的方法，通過假定中尺度風速脈動與hanna參數化方法覆蓋的湍流脈動無關，解一個獨立的拉格朗日方程來計算，所用到的時間尺度和速度的方差是通過對一個測站的風觀測時間序列進行譜分析得到。
進一步，步驟5)中，所述放射性衰變通過下列公式計算：n(1)＝n(0)*0.5δt*λ，
其中，n(1)為衰減后的放射性物質量；n(0)為衰減前的放射性物質量；△t為衰變時間；λ為衰變常數。
進一步，步驟5)中，所述干、濕沉積通過下列公式計算：cd＝c*vd，
其中，cd為沉降后的濃度；c為空氣濃度；vd為沉降速度。
進一步，步驟6)中，所述網格濃度通過下列公式計算：
其中，q，釋放的放射性活度；n，釋放的粒子總數；t，所有粒子擴散時間的總和；nik，第i個粒子在k網格中穿行時間內的時間步數；δt，時間步長；vk，網格k的體積。
本發明的有益技術效果在于：本發明采用粒子隨機游走的方法模擬大氣擴散，把每個污染質點當成有標志的質點，通過釋放大量粒子，計算粒子的軌跡，而這些粒子描述了氣載污染物在大氣中的遷移擴散；本發明粒子在流場中按平均風輸送，同時又用一系列隨機位移來模擬湍流擴散，這樣就表達了平流和湍流擴散兩種作用，對于中長距離遷移問題還考慮了中尺度風脈動的影響，最后由這些粒子在空間和時間上的總體分布估算出污染物的分布，從而保證數據的準確性；由此，在發生突發事件時，決策人員能夠面對該類事件，快速、科學、有效的提出決策建議與對策，從而保護公眾與環境的安全。
附圖說明
圖1是本發明氣載放射性核素長距離遷移軌跡拉格朗日粒子擴散計算方法的流程圖。
具體實施方式
下面結合附圖，對本發明的具體實施方式作進一步詳細的描述。
如圖1所示，是本發明提供的氣載放射性核素長距離遷移軌跡拉格朗日粒子擴散計算方法，該方法包括如下步驟：
1)獲取天氣預報數據、地形數據、下墊面特征數據、模式參數；
天氣預報數據包括：(a)空間三維場數據：兩個水平風分量u和v，一個垂直風分量w，溫度t，特征濕度；(b)空間二維場數據：地表壓力，雪厚度，海面壓力，云狀，10m高u、v風速，2m高溫度，2m高露點溫度，大尺度降水，地面感應熱通量，太陽輻射，地面應力。
模式參數包括：(a)釋放源項參數，即爆炸當量、核素種類及放射性活度；(b)模擬區域范圍，并對該區域范圍進行網格劃分；(c)坐標系類型、時間長度、起止時間。
2)空間坐標和時間坐標變換，是通過已知的方法將等壓層氣象數據轉換為隨地笛卡爾坐標。
3)調用風場；按照文件給定順序讀取風場文件名稱、各風場時刻、以及各時刻各節點氣象數據信息。
4)通過積分粒子運動的拉格朗日方程，計算粒子運行軌跡；利用風場插值，獲得精確的風場數據。
拉格朗日粒子彌散模式是把每個污染質點當成有標志的質點，通過釋放大量粒子，計算粒子的軌跡，而這些粒子描述了氣載污染物在大氣中的遷移擴散。粒子在流場中按平均風輸送，同時又用一系列隨機位移來模擬湍流擴散，這樣就表達了平流和湍流擴散兩種作用，最后由這些粒子在空間和時間上的總體分布估算出污染物的分布。
積分粒子運動的拉格朗日方程，將粒子的運行軌跡寫成下列形式：
xi(t+δt)＝xi(t)+vi(xi(t),t)δt+vi′(xi(t),t)δti＝1,2,3(1)
其中，xi為粒子的三維坐標分量；vi為平均風速分量
u′(t+δt)＝u′(t)r+(1-r2)1/2σuξ
v′(t+δt)＝v′(t)r+(1-r2)1/2σvξ(2)
其中，式中右邊第二項代表速度漲落量中的隨機部分，ξ為符合高斯分布(平均值為0、標準差為σi)的隨機數；
方程(1)考慮了平均風輸送影響和大氣湍流風脈動的影響，湍流脈動反映了時間尺度小于1小時，對應于較短的長度尺度。中尺度運動可以使彌散的煙羽顯著增大(gupta等，1997)，對于大尺度模擬問題，需要考慮中尺度風脈動影響。因此，考慮中尺度風脈動影響的粒子運動方程為：
xi(t+δt)＝xi(t)+vi(xi(t),t)δt+vi′(xi(t),t)δt+vi″(xi(t),t)δti＝1,2,3(3)
其中，vi″為中尺度風脈動速度分量(u″,v″,w″)。
本發明運用拉格朗日粒子彌散模擬方法的關鍵在于確定決定粒子運動平均軌跡的平均風場、三個速度分量的拉格朗日時間尺度和漲落的標準差。
對平均風場，通過解二階差分形式的軌跡運動方程，
對于湍流參數，采用hanna(1982)提出的一種參數化方法，即根據邊界層參數混合層高度h,莫寧-奧布霍夫長度l,對流速度尺度w*,粗糙長度z0和摩擦速度u*來計算湍流參數。由于hanna的方法在整個行星邊界層并不總能得到光滑的σw廓線，導致粒子不能充分混合，因此采用ryall和maryon(1997)提出的一種修正方法確定σw。
對于邊界層參數(如l，u*)，利用模式第一層高度、地面10m和2m高度的風溫數據然后用廓線方法計算上述參數，采用迭代的方法解下列方程：
其中，κ，karman常數；zi，模式第一層高度；δu，模式第一層高度和10m高之間的風速差；δθ，模式第一層高度和2m高之間的位溫差；g，重力加速度；θ*，溫度尺度；l，莫寧-奧布霍夫長度；
其中，φ1和φ2分別為風速和溫度的廓線函數，其形式分別為
對于中尺度風脈動參數，采用maryon(1998)提出的方法：通過假定中尺度風速脈動與hanna參數化方法覆蓋的湍流脈動無關，解一個獨立的拉格朗日方程來解決，所用到的時間尺度和速度的方差是通過對一個測站的風觀測時間序列進行譜分析得到的。假定在網格尺度上觀測的風的方差也為亞網格風速標準差提供一些信息，這樣，用一種簡化方法確定風速標準差的取值方法是：計算粒子所在位置周圍16個網格點(在時間和空間上)的速度的標準差，然后取該標準差的一半作為解中尺度拉格朗日方程時所用的風速標準差。
5)采用源耗減的概念來計算干、濕沉積
給定物質的干沉積用沉積速度vd(m/s)來描述，可根據核素類型和下墊面特征來確定。雨水沖洗效應是造成地面高放射性污染水平的最重要的因素之一。濕沉積可以用類似干沉積的方法計算，不同之處僅在于用沖洗系數λ(s-1)代替干沉積速度，沖洗系數的大小取決于降雨強度。放射性衰變的估算：根據公式
n(1)＝n(0)*0.5δt*λ
其中，n(1)為衰減后的放射性物質量；n(0)為衰減前的放射性物質量；△t為衰變時間，s；λ為衰變常數，s-1。
干、濕沉積量的估算。只針對氣溶膠進行沉積計算。
cd＝c*vd
cd為沉降后的濃度，bq/m2；
c為空氣濃度，bq/m3；
vd為沉降速度，m-1。
當近地層有降水發生時，同時計算干沉積與濕沉積量；
當近地層無降水發生時，只計算干沉積。
6)統計粒子密度、計算網格濃度
為了提高擴散計算精度和效率，當采用粒子分裂技術和核函數方法。利用核函數方法對空間關注點濃度進行計算。核函數方法認為離散化后的煙團在遷移過程中，其自身也呈現高斯擴散的趨勢。對本模式而言，由于空間網格間距較大，而單純的粒子模式如果不使用網格嵌套技術，則網格內部呈現不合理的均勻濃度，并引起監測點濃度的系統誤差。所以，在模式中引入核函數方法，這樣不僅可以避免復雜的嵌套計算，而且可以快速合理地得到有限關注點的濃度。同時，對于一般意義的核函數方法，不強調單個質子團的空間擴散方向性，但由于大尺度模式中，質子團稀疏，對稱分布的核函數仍然會引起不真實的空間濃度分布。對此，模式根據氣象條件設定不同參數，使質子團的核函數擴散呈現更為真實的三維不對稱分布。
放射性核素空氣濃度的計算方法：每個網格中的濃度正比于質點通過該網格時所需時間的總和，因此每個網格的濃度ck(bq/m3)用下列公式計算：
其中，q，釋放的放射性活度，bq；n，釋放的粒子總數；t，所有粒子擴散時間的總和，s；nik，第i個粒子在k網格中穿行時間內的時間步數；δt，時間步長，s；vk，網格k的體積，m3。在濃度計算中，根據放射性核素的半衰期和粒子遷移的總時間對放射性衰變進行修正。
7)輸出時間間隔
如果達到設定的濃度輸出時間，計算該時刻網格的濃度，并輸出空氣濃度場；如果沒有達到設定的濃度輸出時間，考慮粒子分裂條件，跟蹤所有粒子、檢查是否有新粒子釋放；如果有新粒子釋放返回步驟3)，如果沒有新粒子釋放，結束程序。
其中：粒子分裂條件，是當某粒子所在空間的12個相鄰空間均無粒子分布時，對粒子進行分裂。即，分裂后的粒子質量減半，總模擬粒子數增加1。跟蹤所有的粒子，檢查是否有新粒子釋放，是判斷模擬時間是否處于釋放時間段內，如果處于釋放段內，則新增n個粒子的初始遷移位置，并與其它已有粒子一同進行遷移擴散模擬。
對于方程(3)求解，當給定初始條件，通過迭代計算就可以確定各隨機游走粒子的運動軌跡。
以下是本發明在計算過程中所涉及到的數據：
(1)坐標系
模式采用混合坐標系，即x、y、η坐標，其中η是一種對氣壓坐標進行變換的垂直坐標，η與氣壓坐標的變換方法為：
pk＝ak+bkps
ηk＝ak/p0+bk(8)
其中，ηk為模式第k層的η值；ps為地表壓力；p0為壓力常數(101325pa)。ak和bk為系數，由最接近地表(隨地坐標)的值和最大壓力高度層的值確定，中間高度層的系數值則根據隨地表層和壓力高度層之間的壓力梯度確定。
此外，如果緯度超過75°時考慮極地立體投影。
(2)時間、空間分辨率
模式的時間、空間分辨率取決于兩方面因素：一是輸入資料的時間、空間分辨率；二是模式計算的時間、空間分辨率，在輸入資料的時空分辨率一定的情況下，模式計算的時空分辨率主要取決于時間步長。
①輸入資料的時間、空間分辨率
在數值模式中，地形和氣象數據的空間分辨率是一致的，可接受的空間分辨率從0.25°到2.5°，時間分辨率從3h到12h。
數值計算分析結果表明，隨著風場網格分辨率的降低，數值模擬結果的準確程度也隨之降低。數值模擬實驗表明：風特征大約可以被0.5°(約45km)空間分辨、6小時時間分辨水平的預報風場解析，而利用2.5°(約225km)與12小時的時空分辨預報風場則無法滿足需要。因此，可以認為6小時的風場時間間隔是大尺度模式運行的基本條件。綜合考慮模擬精度和計算資源與計算速度，在本系統中采用0.5°的空間分辨率、3～6小時的時間分辨率水平完全滿足實際應用的需要。
②時間步長
時間步長的確定取決于網格的格距、風速和離散風場的時間間隔，從空間角度遵循：
其中，δxi為空間格距，i代表三個方向；vi分別代表三個方向的風分量。從時間角度遵循：
其中，δt風場為輸入風場的時間間隔，即風觀測或數值天氣預報風場的周期。粒子遷移計算模式的時間步長最大值要小于δt空間和δt時間的最小值，最小時間步長可為1秒。
本發明的氣載放射性核素長距離遷移軌跡拉格朗日粒子擴散計算方法并不限于上述具體實施方式，本領域技術人員根據本發明的技術方案得出其他的實施方式，同樣屬于本發明的技術創新范圍。