蔣金利 諸雪征 顧 進(jìn) 韓朝帥， 李 廷

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

1 引言

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

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

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

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

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

2.1 瞬時擴(kuò)散模型

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

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

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

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

2.2 連續(xù)擴(kuò)散模型

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

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

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

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

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

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

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

2.3 擴(kuò)散參數(shù)

在上述擴(kuò)散過程中，確定擴(kuò)散參數(shù)σx，σy，σz其實是很困難的，一般來說，需要特定條件下的氣象觀測，并進(jìn)行大量的計算工作。但是在實際操作中，往往只有常規(guī)的觀測資料。為了能方便高效地根據(jù)已有資料推算出擴(kuò)散參數(shù)，研究人員進(jìn)行了深入的研究，得出了一些行之有效的方法，包括帕斯奎爾擴(kuò)散曲線法、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)準(zhǔn)《制定地方大氣污染物排放標(biāo)準(zhǔn)的技術(shù)方法》［16］，確定大氣穩(wěn)定度后，通過下面兩個公式來計算擴(kuò)散參數(shù):

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

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

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

3.1 變分算法原理

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

四維變分同化的公式:

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

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

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

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

定義分析增量:

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

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

從而得到增量變分公式:

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

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

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

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

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

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

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

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

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

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

3.2.2 化學(xué)危害變分?jǐn)?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是模型的背景場，在本文中，它是由高斯擴(kuò)散模型輸出的參量；B是背景誤差協(xié)方差矩陣；y為觀測場數(shù)據(jù)，H是觀測算子，它的作用是將觀測場數(shù)據(jù)的類型映射到與狀態(tài)量相一致；R為觀測誤差協(xié)方差矩陣。

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

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

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

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

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

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

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

假設(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)在非直接觀測情況中。代表性誤差一般來源于兩個方面:一是模型不完美，得不到完全準(zhǔn)確的測量；二是來自有限的模型分辨率，分辨率3km的觀測系統(tǒng)肯定無法反映10km范圍的區(qū)域狀態(tài)。這類誤差的減小可以利用高密度觀測，提高模式分辨濾等手段。儀器誤差是出廠時帶有的，一般會帶有數(shù)據(jù)參數(shù)，考慮設(shè)備老化，使用時間長等實際情況可以對參數(shù)進(jìn)行適當(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是正定矩陣。由擴(kuò)展子空間原理可知，若果矩陣G有n個共軛方向，那么從任意的初始點出發(fā)，至多做n次精確一維搜索，就可以得到目標(biāo)函數(shù)唯一的極小點。根據(jù)這一理論，可以導(dǎo)出對于正定二次矩陣型函數(shù)的共軛梯度法。用g表示梯度，d表示下降方向，算法步驟如下:

4 案例分析

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

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

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

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

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

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

圖3 濃度場比較圖

從圖中可以分析出來，剛剛開始擴(kuò)散時所預(yù)測出來的濃度圖，本試驗的背景場，在整個擴(kuò)散區(qū)域內(nèi)都有著或多或少的濃度分布，但是30min擴(kuò)散后，除危害物質(zhì)初生云之外，其他地域的物質(zhì)濃度似乎為零；同時，加上隨機(jī)誤差后，可以明顯的看出預(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次迭代之后，同化基本達(dá)到效果，各曲線之間擬合程度加大，開始收斂。這也與圖5中分析200次迭代后同化過程開始收斂的結(jié)果相一致。危害區(qū)域的濃度信息經(jīng)過同化后，危害物質(zhì)濃度逐步逼近真值，圍繞其上下波動。并且隨著時間的推進(jìn)和觀測數(shù)據(jù)的增加，同化結(jié)果越來越穩(wěn)定，逐漸消除了初始時刻所加的10%的隨機(jī)誤差，達(dá)到收斂狀態(tài)。

5 結(jié)語

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