基于變分同化的化學(xué)危害擴散預(yù)測研究?

蔣金利 諸雪征 顧 進 韓朝帥， 李 廷

（1.陸軍防化學(xué)院 北京 102205）（2.軍事科學(xué)院防化研究院 北京 102205）

1 引言

當(dāng)戰(zhàn)場發(fā)生化學(xué)襲擊或者化工廠發(fā)生泄漏，爆炸事故后，毒劑被分散成為蒸氣、氣溶膠進入大氣中，形成毒劑云團。而毒劑云團將隨氣流一起運動，同時向四周擴散，形成化學(xué)危害區(qū)域［1］。指導(dǎo)部隊進行有效化學(xué)防護的前提是對戰(zhàn)場化學(xué)危害區(qū)域進行準確的界定、檢測及預(yù)測，只有掌握毒劑濃度的時空分布規(guī)律，準確描繪毒劑云團的危害地域，才能預(yù)測未來一段時間窗口內(nèi)的毒劑濃度分布，進而分析出該段時間內(nèi)的毒害作用效應(yīng)。當(dāng)前對于化學(xué)危害的預(yù)測主要通過大氣擴散理論和數(shù)學(xué)擴散模型［2］來推導(dǎo)未來一段時間的毒劑濃度分布，從而對化學(xué)危害的發(fā)展態(tài)勢做出評估，但是基于大氣擴散模式的化學(xué)云團濃度預(yù)測具有一定的偏差，且戰(zhàn)場氣象、地形條件復(fù)雜多變，計算結(jié)果往往會有較大的誤差。

數(shù)據(jù)同化是一種最初來源于數(shù)值天氣預(yù)報的數(shù)據(jù)處理技術(shù)，現(xiàn)在已成功用于對大氣污染物濃度的預(yù)測中來。在化學(xué)危害預(yù)測領(lǐng)域引入數(shù)據(jù)同化算法，以更快的速度和更高的精準度得出危害場未來一段時間的濃度走勢及分布圖，在保證速度的同時盡最大可能地為防護和救援提供準確的決策依據(jù)。當(dāng)前關(guān)于污染物擴散數(shù)據(jù)同化主要有逐步迭代法［3］，最優(yōu)插值法［4］、卡爾曼濾波法［5］、集合卡爾曼濾波法［6］、遺傳算法［7］、神經(jīng)網(wǎng)絡(luò)算法［8］和機器學(xué)習(xí)［9］等，但逐步迭代法和最優(yōu)插值法精度偏低，屬于早期的同化算法；卡爾曼濾波系列算法對集合數(shù)目的要求偏高，難以快速應(yīng)用于危害環(huán)境；遺傳算法，神經(jīng)網(wǎng)絡(luò)和機器學(xué)習(xí)是新興智能優(yōu)化算法，用于數(shù)據(jù)同化領(lǐng)域還處于初期探索階段。所以，本文選用變分同化［10～11］方法，既避免了精度不夠的問題，也沒有集合卡爾曼濾波對數(shù)據(jù)樣本這么高的要求。

因此，本文結(jié)合變分同化理論，結(jié)合經(jīng)典的高斯擴散模式，對化學(xué)危害擴散中數(shù)據(jù)同化算法進行構(gòu)建和案例分析。

2 化學(xué)危害擴散模型

高斯擴散模式［12］是經(jīng)典的污染物擴散模型。在大氣中，幾乎處處時時都存在著各種不同尺度的湍流運動，常見的湍流理論主要有梯度輸送理論和統(tǒng)計理論兩種。但是無論從哪一種理論出發(fā)，高斯擴散模型是可推導(dǎo)的。由湍流擴散梯度輸送理論建立起來的擴散方程，當(dāng)擴散系數(shù)為常數(shù)時，可以得出正態(tài)分布形式的解；當(dāng)氣流平穩(wěn)均勻時，擴散物質(zhì)的粒子分布也遵循正態(tài)分布。盡管實際大氣中，并不能嚴格滿足相關(guān)條件，但是在不知道實際物質(zhì)濃度分布的情況下，高斯分布起碼是一種較好的近似。在實際操作中，高斯擴散模式所需要的氣象條件少，計算速度快，也可以較好地反映濃度的時空分布。

2.1 瞬時擴散模型

瞬時釋放又稱為高斯煙團模式［13～14］，可以看成是將高斯煙羽模式分成數(shù)個煙團。瞬時點源排放的遷移模式的基本方程:

同時，假定在原點（0，0，0），即x=y=z=0處有一個瞬時點源排放；當(dāng)時間趨近零時，濃度也趨近零，當(dāng)時間趨向無窮大時，濃度也趨近于零；整個擴散過程，滿足質(zhì)量守恒定律，且假定時間為零時，排放物的總質(zhì)量為M。

在滿足上述條件之后，微分方程式（1）的解為

式中:C（x，y，z，t）為 t時刻，空間中任一點（x，y，z）的濃度，單位為（g/m3）；Qi，t為單個煙團的濃度；σx，i，t，σy，i，t，σz，i，t為t時刻各方向上的擴散參數(shù)；xi，t，yi，t，zi，t為某個毒劑煙團的位置參數(shù)。

2.2 連續(xù)擴散模型

點源擴散是高斯擴散模式的主要研究對象之一，它是指在擴散進行了一段持續(xù)時間且氣象狀況平穩(wěn)（風(fēng)向恒定）的條件下所進行的擴散。在這個過程中，我們把釋放源的位置作為坐標(biāo)原點，且源強連續(xù)、均衡地釋放毒劑氣體或氣溶膠，釋放出來的氣體延風(fēng)向不斷進行擴散運動。如果以風(fēng)向方向為x軸，垂直方向為y軸建立坐標(biāo)系的話，不難預(yù)見，毒劑氣體將在下風(fēng)方向形成圓錐形的云團，并且在x軸的任一截面處都符合二維正態(tài)分布。連續(xù)擴散下，高斯擴散模式又叫高斯煙羽模型［14］。

地面點源連續(xù)排放的粒子遷移工程為

上式中各參數(shù)含義與瞬時擴散時相同。在滿足質(zhì)量守恒定律的前提下，該微分方程的通解為

按照所排放毒劑的性質(zhì)的不同，標(biāo)準解一般可分為兩種情況，即地面反射與地面不反射。

在擴散過程中，當(dāng)毒劑沉降到地面時，存在反射回大氣層或者不沉降到地面的情況，數(shù)學(xué)上的表達為滿足積分:

將式（5）代入式（6），可得:

其中，h為源項與地面的高度，如果是地面源，那么其值為零。其他項的含義和瞬時擴散方程相同。

2.3 擴散參數(shù)

在上述擴散過程中，確定擴散參數(shù)σx，σy，σz其實是很困難的，一般來說，需要特定條件下的氣象觀測，并進行大量的計算工作。但是在實際操作中，往往只有常規(guī)的觀測資料。為了能方便高效地根據(jù)已有資料推算出擴散參數(shù)，研究人員進行了深入的研究，得出了一些行之有效的方法，包括帕斯奎爾擴散曲線法、GB/T 公式以及 Briggs［15］公式等方法。這些方法將大氣劃分A-F為六個穩(wěn)定度級別，穩(wěn)定度劃分如表1所示。

表1 大氣穩(wěn)定度分級

Pasquill和Gifford也給出的不同穩(wěn)定度下的σy與σz隨下風(fēng)距離變化的經(jīng)驗曲線。在P-G曲線的基礎(chǔ)上，我國根據(jù)實際情況制訂了國家標(biāo)準《制定地方大氣污染物排放標(biāo)準的技術(shù)方法》［16］，確定大氣穩(wěn)定度后，通過下面兩個公式來計算擴散參數(shù):

式中，x為下風(fēng)方向距離；a1，a2及γ1，γ2為擴散參數(shù)，都是無量綱常數(shù)。

3 化學(xué)危害擴散數(shù)據(jù)同化模型

根據(jù)變分同化理論，結(jié)合高斯擴散模式，對戰(zhàn)場化學(xué)危害問題進行數(shù)學(xué)建模。變分同化方案的基本功能:在化學(xué)危害事件發(fā)生之后，結(jié)合當(dāng)時的氣象條件，觀測數(shù)據(jù)等，輸入變分同化的工作流程。當(dāng)誤差統(tǒng)計減小到可接受的范圍后，輸出危害場的濃度分析值，反映當(dāng)時狀態(tài)下的危害程度，為后續(xù)的救援行動以及事后評估提供決策依據(jù)。下面依據(jù)變分理論，進行方案設(shè)計，構(gòu)建基于高斯擴散模型的戰(zhàn)場化學(xué)危害同化模型。

3.1 變分算法原理

四維變分方法（Four—Dimensional Variational Algorithm:4D-Var）最早由我國學(xué)者顧震潮［17］提出，是以三維變分法方法（3D-Var）為基礎(chǔ)發(fā)展起來的，是3D-Var在時間維上的擴展，實現(xiàn)了整個同化時間窗口觀測數(shù)據(jù)的同化，預(yù)測精度更高，適用范圍更廣。

四維變分同化的公式:

式中Hi表示第i時刻的觀測算子，表示第i時刻的觀測值，Mt0，ti表示將預(yù)報模式從時刻積分到ti時刻的模式預(yù)報算子，故:

指由模式算子積分得到的第i時刻的模式變量。同理可得:

指將第i時刻的模式變量經(jīng)第i時刻的觀測算子Hi映射到觀測場。

1994年，Courtier等為了避免大量的計算，降低運算成本，提出增量方法［18］，它不直接求解分析變量，而是求解分析增量并修正預(yù)測場。

定義分析增量:

同時，對預(yù)報模式算子進行泰勒展開:

代表觀測算子的切線模式或其他線性近似模式。舍去式中的高階無窮小項，并令向量:

從而得到增量變分公式:

3.2 化學(xué)危害擴散變分同化方案設(shè)計和模型建立

3.2.1 化學(xué)危害變分同化方案設(shè)計

依據(jù)變分同化原理，結(jié)合高斯擴散模式以及化學(xué)危害初生云形式，設(shè)計了變分數(shù)據(jù)同化方案，對工作流程進行了梳理。

該方案主要包括六個部分:外部條件獲取模塊，源項信息估算模塊，觀測數(shù)據(jù)處理模塊，高斯擴散預(yù)測模塊，變分同化模塊以及結(jié)果分析模塊。

外部條件獲取模塊:獲取相關(guān)的風(fēng)場信息，包括風(fēng)向、風(fēng)速、溫度等；同時確定危害地域的地形地貌信息，有無大型山體、農(nóng)田等。

源項信息估算模塊:主要是為高斯擴散模型提供源項信息，生成同化背景場。

觀測數(shù)據(jù)處理模塊:處理各類裝備收集到的各種觀測資料，主要對濃度以及風(fēng)場的數(shù)據(jù)進行處理，在同化時間窗口內(nèi)生成觀測向量以及數(shù)據(jù)文件。

高斯擴散預(yù)測模塊:將所獲得的信息條件輸入，進行預(yù)測，生成背景場信息，傳輸給同化模塊。并不斷更新時間域，不斷向前預(yù)測。

變分同化模塊:綜合處理背景資料，觀測資料，通過所建立的變分同化模型得出分析場。

結(jié)果分析模塊:輸出分析結(jié)果，并分析過程誤差。

3.2.2 化學(xué)危害變分數(shù)據(jù)同化模型構(gòu)建

1）變分同化梯度函數(shù)

變分同化算法的方案設(shè)計部分，是在空間場內(nèi)求解最優(yōu)分析解的數(shù)據(jù)同化算法。下式是上一章節(jié)提到的變分法目標(biāo)函數(shù)的基礎(chǔ)公式:

對標(biāo)量函數(shù)J關(guān)于自變量x求導(dǎo)，可以得到其梯度:

求解梯度公式，可得分析值:

式中，x為需要求解的狀態(tài)參量，xa稱為分析值（analysis value）；xb是模型的背景場，在本文中，它是由高斯擴散模型輸出的參量；B是背景誤差協(xié)方差矩陣；y為觀測場數(shù)據(jù)，H是觀測算子，它的作用是將觀測場數(shù)據(jù)的類型映射到與狀態(tài)量相一致；R為觀測誤差協(xié)方差矩陣。

誤差是整個同化過程中不可避免的話題，從數(shù)學(xué)的角度來說，可以把同化過程理解為最小化誤差的過程。誤差在同化開始前主要體現(xiàn)在背景誤差協(xié)方差矩陣B上，它是整個同化過程的重要組成部分，和觀測誤差協(xié)方差共同決定了分析場的誤差。因此，在同化過程開始之前，必須對誤差進行修正，盡可能使分析值誤差降低，提高準確度。

一般情況下，確定誤差協(xié)方差矩陣前要有一個假設(shè)，即觀測誤差與背景誤差不相關(guān)，即:

式中:εb為背景誤差；εo為觀測誤差，上標(biāo)T表示矩陣的轉(zhuǎn)置。在實際分析中，這也是完全合理的假設(shè)，背景場參量與觀測場參量相互獨立。

2）背景誤差協(xié)方差矩陣B

背景誤差協(xié)方差矩陣表示了預(yù)報誤差的概率分布函數(shù)（Probability Distribution Function，PDF），是一個高斯型的函數(shù)分布。它的準確與否直接決定，其同化系統(tǒng)所產(chǎn)生的分析場參量大概只有15%來自觀測資料，而剩余的85%均來自背景場［19］。在物理意義上，背景場是同化前的該地域的狀態(tài)場，它將同化循環(huán)中的觀測數(shù)據(jù)不斷融合進最新的分析場。所以說，如何對分析場進行準確有效的更新，不僅取決于模型的質(zhì)量，也取決于背景場質(zhì)量。其定義如下:

式中xtrue代表的是真值，εb為背景誤差，由于很多時刻真值未知，所以，很多時候B矩陣也表示為

式中，上橫線代表均值，其余含義與上式相同。由此可以看出，B矩陣是正定對稱矩陣，而且特征值均為非負值。

假設(shè)背景誤差為向量(e1,e2,e3)，那么其背景誤差協(xié)方差矩陣可以寫成:

從上面這種形式下的B矩陣可以看出，其中的每個元素反映的都是兩個位置上背景誤差的相關(guān)性，和兩個變量之間的相關(guān)性。對角線元素或者對角矩陣代表的是背景誤差方差，而非對角元素和矩陣代表的則是協(xié)方差。

定義x，x0兩點間的背景場相關(guān)系數(shù)b(x,x0)為兩點間距離dist(x-x0)的函數(shù):

3）觀測誤差協(xié)方差矩陣R

R是觀測誤差協(xié)方差矩陣，它包含了關(guān)于觀測誤差的統(tǒng)計信息。對于很多觀測系統(tǒng)來說，不同位置上的觀測誤差在統(tǒng)計上是相互獨立的。同樣地，在確定觀測誤差協(xié)方差矩陣之前，假定觀測誤差是無偏的，εo的數(shù)學(xué)期望為零。

觀測誤差一般包括觀測算子誤差、儀器誤差和代表性誤差。觀測算子誤差主要是因為對兩種不同種類的數(shù)據(jù)之間聯(lián)系的理解不足或者理解過于片面導(dǎo)致的，主要出現(xiàn)在非直接觀測情況中。代表性誤差一般來源于兩個方面:一是模型不完美，得不到完全準確的測量；二是來自有限的模型分辨率，分辨率3km的觀測系統(tǒng)肯定無法反映10km范圍的區(qū)域狀態(tài)。這類誤差的減小可以利用高密度觀測，提高模式分辨濾等手段。儀器誤差是出廠時帶有的，一般會帶有數(shù)據(jù)參數(shù)，考慮設(shè)備老化，使用時間長等實際情況可以對參數(shù)進行適當(dāng)調(diào)整。除此之外，就目前的研究水平而言，國內(nèi)外也鮮有好的解決辦法。一般而言，對R的處理相當(dāng)簡化，在觀測點之間相互獨立的情況下，可以將其簡化為對角矩陣處理，方便下一步的計算。

3.2.3 最優(yōu)化算法

共軛梯度法是介于最速下降法與牛頓法之間的一個方法［20］，它僅需利用一階導(dǎo)數(shù)信息，但克服了最速下降法收斂慢的缺點，又避免了牛頓法需要存儲和計算Hessen矩陣并求逆的缺點。本文結(jié)合化學(xué)危害問題的研究背景以及變分同化的計算要求，選定共軛梯度法為最優(yōu)算法。

對于一般二次型目標(biāo)函數(shù)來說:

其中G是正定矩陣。由擴展子空間原理可知，若果矩陣G有n個共軛方向，那么從任意的初始點出發(fā)，至多做n次精確一維搜索，就可以得到目標(biāo)函數(shù)唯一的極小點。根據(jù)這一理論，可以導(dǎo)出對于正定二次矩陣型函數(shù)的共軛梯度法。用g表示梯度，d表示下降方向，算法步驟如下:

4 案例分析

為了驗證所建模型的同化能力以及各參數(shù)對同化質(zhì)量的影響，本節(jié)利用數(shù)值模擬實驗進行分析。

首先運行危害物質(zhì)的擴散模型，得到危害區(qū)域內(nèi)各個網(wǎng)格點的濃度，將其作為背景濃度場，并將此時即為記為同化窗口的0時刻；繼續(xù)向前運行擴散模型，簡單起見，將所得到的模型預(yù)測值作為“真值”；其次，在危害區(qū)域內(nèi)設(shè)定觀測點位置，將觀測點位置處預(yù)測濃度輸出并記為觀測數(shù)據(jù)；同時，保持參數(shù)不變，運行另一組擴散模型，并給初始預(yù)測結(jié)果加上10%隨機誤差，用上述得到觀測值同化帶有誤差的預(yù)測值，檢驗同化效果以及收斂性。

條件設(shè)定:假設(shè)某平原地區(qū)化工廠發(fā)生煙囪內(nèi)單一化學(xué)物質(zhì)泄漏，并且已獲取當(dāng)時的源項信息以及氣象信息，廠內(nèi)為預(yù)防事故已經(jīng)配備了具有輸出功能的觀測設(shè)備，其他條件如表2所示。

表2 實驗條件設(shè)定

為了使實驗盡量簡單，設(shè)定釋放方式為均勻釋放，并且理想化外部條件，不考慮降雨、雷電等惡劣天氣影響。除同化效果分析之外，本節(jié)還將分析收斂性以及過程誤差。

首先，給出背景場，真實場以及模型預(yù)測場的擴散圖。

圖3 濃度場比較圖

從圖中可以分析出來，剛剛開始擴散時所預(yù)測出來的濃度圖，本試驗的背景場，在整個擴散區(qū)域內(nèi)都有著或多或少的濃度分布，但是30min擴散后，除危害物質(zhì)初生云之外，其他地域的物質(zhì)濃度似乎為零；同時，加上隨機誤差后，可以明顯的看出預(yù)測值與真實值之間有了差異。

取迭代次數(shù)為1000次，下圖是試驗的同化效果圖，收斂性以及誤差變化。

圖4 同化實驗結(jié)果圖

從圖中可以看出，通過同化，分析值比觀測值跟更接近真值，而背景值在初始狀態(tài)下與其他三條曲線有著明顯偏差，在一定迭代次數(shù)后，逐漸與其他幾條曲線擬合；由于觀測誤差協(xié)方差矩陣R為對角矩陣，且對角線元素為均方差，為固定值，所以不管迭代次數(shù)是多少次，同化時間窗口有多長，觀測場方差都不變；同化過程的單次誤差雖然各次不同，但是在迭代200次左右開始，便趨于收斂，變的穩(wěn)定；同樣地，同化的背景場方差也是隨著迭代次數(shù)的增多，開始逐漸收斂，次數(shù)也大概在200次左右。

圖5是迭代次數(shù)20，50，100，200次時的同化圖。

圖5 不同迭代次數(shù)下的同化效果

從圖5中對比分析可以看出，在200次迭代之后，同化基本達到效果，各曲線之間擬合程度加大，開始收斂。這也與圖5中分析200次迭代后同化過程開始收斂的結(jié)果相一致。危害區(qū)域的濃度信息經(jīng)過同化后，危害物質(zhì)濃度逐步逼近真值，圍繞其上下波動。并且隨著時間的推進和觀測數(shù)據(jù)的增加，同化結(jié)果越來越穩(wěn)定，逐漸消除了初始時刻所加的10%的隨機誤差，達到收斂狀態(tài)。

5 結(jié)語

本文主要對四維變分數(shù)據(jù)同化在化學(xué)危害擴散預(yù)測中的應(yīng)用進行了系統(tǒng)研究和驗證。首先，分別對高斯擴散模式的瞬時釋放模式和連續(xù)釋放模式進行了研究，提出適用于化學(xué)危害預(yù)測的擴散參數(shù)；其次，構(gòu)建了適用于化學(xué)危害預(yù)測的變分同化算法，對同化方案、背景誤差協(xié)方差矩陣、觀測誤差協(xié)方差矩陣以及最優(yōu)化算法進行了建模，明確了尋優(yōu)步驟；最后，通過Matlab仿真，對算法的同化質(zhì)量和收斂性進行了試驗和對比分析。但是，由于為證明參數(shù)的影響，示例分析部分應(yīng)用了理想化數(shù)值實驗，簡化了實驗條件，接下來的研究還要更加細化實驗條件，做更細致的分析。