基于貝葉斯推斷的時變流場下污染源反演

朱建杰，周晅毅，顧 明

（同濟大學土木工程防災國家重點實驗室，上海 200092）

當公共場所發(fā)生毒氣襲擊或者化工廠發(fā)生?；沸孤r，需要快速對泄漏源進行定位和強度評估，從而給疏散民眾和制定應急措施提供指導。由于事故發(fā)生時，流場和濃度場往往具有隨時間變化的特性，因此在時變流場下進行污染源參數(shù)反演具有重要的應用價值。

污染源參數(shù)反演的主要方法有直接求解法［1-2］、優(yōu)化方法［3-4］和概率方法［5-6］。直接求解法通過構造控制方程的反問題，直接求解源參數(shù)。Zhang等［1］利用Tikhonov 正則化法［7］反演了污染源強度，得到了污染物釋放速率的變化曲線。Wei 等［2］將Tikhonov正則化法和貝葉斯方法相結合，得到了多個污染源的污染物釋放速率和污染源位置的概率分布。優(yōu)化方法通過構造目標函數(shù)，求解使目標函數(shù)最小的源參數(shù)。江思珉等［3］利用模擬退火算法反演了地下水污染源的源強度。曾令杰和高軍［4］通過遺傳算法對空調系統(tǒng)引發(fā)的污染進行了快速溯源。貝葉斯概率方法考慮了測量濃度和模擬濃度之間的誤差，將誤差量化為似然概率，通過貝葉斯公式計算相應源參數(shù)的概率分布。Keats 等［5］利用伴隨方程求解傳感器的伴隨濃度場，使模擬濃度的計算量大大降低，得到了非時變濃度場下污染源位置和強度的概率分布。Guo 等［6］考慮了污染物濃度場從零濃度狀態(tài)到穩(wěn)定狀態(tài)的變化過程，利用MCMC（Markov chain Monte Carlo）抽樣算法直接對符合后驗概率的樣本參數(shù)進行抽樣，獲得了污染源的位置和強度信息，然而該研究的污染源所在的流場仍然是不隨時間變化的。

本文在Keats［5］、Guo［6］的研究基礎上，將基于貝葉斯推斷的污染源參數(shù)反演方法擴展至時變流場的場景?？紤]了流場變化對濃度場的影響，求解伴隨方程得到了不同污染源參數(shù)下傳感器位置處的模擬濃度，再通過貝葉斯公式計算相應污染源參數(shù)的后驗概率。比較了不同位置污染源的參數(shù)反演效果，還討論了利用污染物擴散不同階段的測量濃度數(shù)據(jù)時，不同測量數(shù)據(jù)對反演效果的影響。

1 貝葉斯推斷

1.1 貝葉斯公式

貝葉斯公式可表示為

式中：P（m|d）為后驗概率，是已知測量濃度為d條件下，污染源參數(shù)的概率分布；P（m）為先驗概率，是未知測量信息情況下的污染源參數(shù)概率分布；P（d|m）為似然概率，表示在污染源參數(shù)為m的條件下，模擬濃度c與測量濃度d之間的相似程度。測量濃度d為傳感器測量得到的濃度，模擬濃度c通過污染物擴散數(shù)值模擬獲得。

對于本文的污染源參數(shù)反演問題，污染源參數(shù)為m=｛xs，ys，zs，qs｝，其中xs，ys，zs為污染源的空間坐標，qs為污染物的釋放速率；測量濃度d=｛d1，1，d1，2，…，di，j，…，dN，M｝，其中di，j表示第i個傳感器第j個測量濃度值，M為每個傳感器記錄數(shù)據(jù)的總數(shù)，N為傳感器的數(shù)量。

1.2 先驗概率

先驗概率是在未得到任何測量信息情況下污染源參數(shù)的概率。這里認為先驗概率服從均勻分布，即：

式中：C為常數(shù)。

1.3 似然概率

似然概率表示污染源參數(shù)為m的情況下，模擬濃度c與測量濃度d之間的相似程度，通常用它們之間誤差的概率來表達。許多學者［5-6，8-13］認為模擬濃度c與測量濃度d之間的誤差符合高斯分布，也有一部分學者［14-16］考慮到不同傳感器之間的濃度值存在量級差，認為利用對數(shù)高斯分布描述誤差可以降低量級差對反演準確性的影響?？紤]到高斯分布在污染源反演研究領域應用更加廣泛，并且具有良好的反演效果，故本文采用高斯分布描述模擬濃度與測量濃度之間的誤差，似然概率為

式中：ci，j（m）表示污染源參數(shù)為m的情況下第i個傳感器第j個模擬濃度值；σi，j為模擬濃度ci，j（m）和測量濃度di，j之間誤差的標準差。假設誤差與測量濃度處于同一量級［11-13］，可以將誤差標準差σi，j設為與測量濃度di，j相同的值，即：

由于反演需要利用不同時刻多個傳感器的濃度信息，假設不同測量值之間相互獨立，可以將所有測量濃度對應的似然概率相乘。于是，式（1）中的似然概率P（d|m）表示為

2 模擬濃度計算方法

利用貝葉斯推斷解決污染源參數(shù)反演問題，往往要計算特定源參數(shù)條件下的模擬濃度ci，j，傳統(tǒng)的模擬濃度計算方法需要解對流擴散方程。為了獲得完整的污染源參數(shù)后驗概率分布，須求解所有可能污染源參數(shù)條件下的對流擴散方程。這種方法的計算量十分龐大，反演的效率非常低。Keats 等［5］提出利用伴隨方程計算模擬濃度，使得模擬濃度計算問題變得高效且結果精確。下面將對時變流場中的伴隨方程進行推導。

2.1 對流擴散方程

對流擴散方程為

式中：Cm為源參數(shù)為m條件下的模擬濃度場；U為污染源所在的流場；K為擴散系數(shù)；S（m）為源項。當污染物釋放速率恒定且污染源為點源時，S（m）=qsh（x-xs，y-ys，z-zs），其中h（·）為delta函數(shù)，表示僅在坐標為xs，ys，zs的位置存在釋放速率為qs的污染源。式（6b）和式（6c）分別為對流擴散方程的邊界條件和初始條件。模擬濃度場Cm與傳感器模擬濃度ci，j（m）的關系為

式中：xi=｛xi，yi，zi｝表示第i個傳感器的空間坐標；tj表示第j個測量數(shù)據(jù)的測量時刻。

由于污染物擴散對流場的影響較小，故對流擴散方程（6a）的流場U不會受到源項變化的影響。在這種情況下，濃度場Cm與污染物釋放速率qs呈線性關系；而濃度場Cm與污染源位置xs，ys，zs呈非線性關系。當污染源位置發(fā)生變化時，需要求解新的對流擴散方程，再通過式（7）計算傳感器位置處的模擬濃度。對流擴散方程的求解是非常耗時的，由于在一次反演中污染源往往存在非常多的潛在位置，若對每一個潛在位置都求解對流擴散方程，反演計算的效率將變得非常低。

2.2 伴隨濃度

為了快速計算任意參數(shù)條件下傳感器位置處的模擬濃度，引入了伴隨濃度C*這一變量。將伴隨濃度C*與源項S（m）相乘，將其乘積在時間和空間上積分得到：

利用分部積分法和高斯公式，可以將式（8）轉化為

伴隨方程（10a）可以理解成以傳感器為源項，在tj時刻存在脈沖的污染物釋放，同時將時間推進方向和流場速度方向進行反向處理后得到的對流擴散方程。當對流擴散方程（6a）滿足邊界條件（6b）和初始條件（6c）并且伴隨方程（10a）滿足邊界條件（10b）和初始條件（10c）時，式（9）中ndSdt=0，可得到伴隨濃度場C*與傳感器模擬濃度ci，j（m）的關系式：

以上為利用伴隨方程計算模擬濃度的推導過程。由伴隨方程（10a）可知，伴隨濃度場C*i，j跟源參數(shù)m無關，它僅跟第i個傳感器的坐標和測量時刻tj相關。當傳感器位置xi和測量時刻tj不變時，計算不同源參數(shù)條件下的模擬濃度只需將相應源項S（m）代入式（11）即可，式（11）只涉及積分運算，而利用對流擴散方程時需要求解偏微分方程，前者的計算量遠遠小于后者，這樣就可以極大地加快反演效率。設每個傳感器測量的次數(shù)為M，傳感器的數(shù)量為N，則所需求解伴隨方程的總數(shù)為M×N。

3 算例分析

3.1 污染物擴散數(shù)值模擬

3.1.1 模擬工況

如圖1 所示，建筑模型尺寸為0.1 m×0.1 m×0.2 m（L·B·H），其周圍存在A、B、C、D、E 5 個不同位置的污染源，均位于地面，用空心方塊表示，對5個污染源分別做污染物擴散數(shù)值模擬。來流剖面遵循指數(shù)律，建筑頂部高度處風速為uH=4.2 m·s-1，指數(shù)α=0.25。10個傳感器位于建筑背風側，用實心黑點表示，編號為1～10，位于距地面0.062 5H高度處。

圖1 污染源與傳感器布置方式Fig.1 Layout of sources and sensors

3.1.2 計算域與模擬設置

根據(jù)日本建筑學會的計算流體動力學（computational fluid dynamics，CFD）指南［17］，計算域（圖2）尺寸設為3.4 m×2.0 m×0.8 m（17H·10H·4H）；計算域入口與建筑模型的距離為4.5H；計算域的側面和出口與建筑模型的距離分別為4.75H和12H。

圖2 計算域Fig.2 Computational domain

流場模擬以及污染物擴散模擬參數(shù)設置如表1所示，湍流模型采用剪切應力輸運（shear stress transfer，SST）k-ω；流場的計算時間為50 s（1 050H/uH）。當流場已經達到穩(wěn)定（t=630H/uH）時，污染物開始釋放。由于該模擬工況與日本東京工藝大學于2006年做的污染物擴散風洞試驗工況類似，故來流平均速度u和湍動能k根據(jù)試驗數(shù)據(jù)確定，釋放速率為qs=5.83×10-6m3·s-1。

表1 邊界條件與數(shù)值模擬設置Tab.1 Boundary conditions and solver settings for simulation case

3.1.3 網格無關性檢驗

為了確保模擬結果不受網格劃分方案影響，需要進行網格無關性檢驗（表2）。計算域采用結構化網格劃分，建筑模型周圍以及傳感器附近區(qū)域進行了網格加密，粗、中、細三種網格的最小單元尺寸分別為H/20、H/30和H/40，網格的伸縮比設置為1.05，總網格數(shù)分別為571 400、1 254 123和2 105 408。

表2 網格無關性檢驗參數(shù)設置Tab.2 Parameter settings of grid independent analysis test

圖3 為污染源C 處釋放的污染氣體在建筑模型下風向0.5H和H位置且y=0 處的模擬濃度剖面，可見僅在建筑下風向0.5H且高度大于1.2H的位置處，粗網格條件下模擬得到的氣體濃度小于中網格和細網格條件下的模擬結果，其他位置處三種網格的模擬濃度曲線重合度較高，故本文選取中網格對污染物擴散進行數(shù)值模擬。

圖3 不同網格劃分方案的平均濃度剖面Fig.3 Profile of mean concentration for different meshing schemes

3.1.4 自保持性驗證

為了保證計算域入口至建筑模型處的來流剖面不會發(fā)生顯著變化，需要進行自保持性驗證。圖4a為空風場條件下，數(shù)值模擬和風洞試驗中的x向平均風速剖面從計算域入口至建筑位置處的變化情況。在數(shù)值模擬結果中，x向平均風速從計算域入口至建筑位置處的變化較小，與風洞試驗的測量結果也較為吻合。圖4b為空風場條件下，數(shù)值模擬和風洞試驗中的湍動能剖面從計算域入口至建筑位置處的變化情況。數(shù)值模擬和風洞試驗得到的湍動能均出現(xiàn)了一定的衰減，相比風洞試驗結果，數(shù)值模擬的湍動能衰減程度更大。Kim和Baik［18］研究了數(shù)值模擬中湍動能對污染物濃度的影響，根據(jù)研究結果可知本文數(shù)值模擬中的湍動能衰減對污染物濃度場的影響不大。

圖4 空風場條件下的x向平均風速與湍動能剖面（y/H=0）Fig.4 Profile of mean velocity in x direction and turbulent kinetic energy in empty wind field (y/H=0)

3.1.5 流場

圖5 為3 號傳感器位置處的x向的風速（u）時程，流場穩(wěn)定后的風速存在有規(guī)律的周期性波動。

圖5 3號傳感器處x向瞬時風速時程Fig.5 Temporal velocity in x direction for Sensor 3

3.1.6 對流擴散方程計算測量濃度

當流場已達穩(wěn)定狀態(tài)（t=630H/uH）時，釋放污染物氣體。圖6 給出了分別由污染源A～E 釋放的污染氣體，擴散至3 號傳感器位置處隨時間變化的濃度，圖中的濃度值經過了量綱一處理，其中c0=qs/（uH·H2）=34.7×10-6。從污染物釋放初期（t=630H/uH）開始，濃度由零逐漸發(fā)展至穩(wěn)定，受到周期脈動的流場影響，污染物的濃度在發(fā)展過程中和穩(wěn)定后均表現(xiàn)出有規(guī)律的周期波動特性。這里通過對流擴散方程（6a）對濃度場進行求解，將式（7）的計算結果作為傳感器的測量濃度d。

圖6 3號傳感器濃度時程Fig.6 Concentration history for Sensor 3

為了在后續(xù)研究中探究利用污染物擴散不同階段的測量數(shù)據(jù)反演源參數(shù)時不同測量數(shù)據(jù)對反演結果的影響，以t=660H/uH為分界將污染物擴散過程劃分成兩個階段。在t=660H/uH之前為發(fā)展階段，在t=660H/uH之后為穩(wěn)定階段。由于3號傳感器的濃度時程最具有代表性，其他傳感器的濃度變化與之類似，故不再畫出其他傳感器的濃度時程圖。

3.1.7 伴隨方程計算模擬濃度

2.2 節(jié)介紹了在計算任意源參數(shù)條件下某一傳感器位置處某一時刻的模擬濃度時，基于伴隨方程的計算方法在節(jié)省計算量方面具有優(yōu)勢，故本文通過求解伴隨方程（10a），再利用式（11）計算模擬濃度。圖7 對比了污染源分別為A～E 時，3 號傳感器位置處的測量濃度和模擬濃度的計算結果。當污染源為A時，模擬濃度與測量濃度相比偏?。划斘廴驹礊锽時，模擬濃度與測量濃度相比偏大；對于其他位置的污染源，兩種濃度的計算結果基本保持一致。雖然污染源分別為A 和B 時，測量濃度和模擬濃度在3 號傳感器位置的結算結果存在差異，但是兩者均表現(xiàn)出了一致的周期性波動趨勢，數(shù)值上的誤差可能是由網格劃分、微分方程離散和四舍五入等因素導致的，這表明基于伴隨方程的模擬濃度方法可以在忽略數(shù)值誤差的情況下達到與基于對流擴散方程時相同的效果。在反演的過程中，這里的誤差可以作為數(shù)值模擬與實際測量的差異。值得注意的是，在基于對流擴散方程計算5 個不同位置污染源產生的污染物擴散至3號傳感器位置處在某一時刻的模擬濃度時，須求解5 種源參數(shù)條件下的對流擴散方程；而基于伴隨方程的方法只需針對3 號傳感器在同一測量時刻求解一次，即可達到相同效果，大大節(jié)省了計算量。

圖7 3號傳感器位置處測量濃度與模擬濃度比較Fig.7 Comparison of measurements and simulated concentrations at location of Sensor 3

3.2 污染源參數(shù)的后驗概率

假設污染源位于地面，即污染源的豎向坐標zs=0，因此僅對污染源的xs、ys坐標和污染源的污染物釋放速率qs進行反演。

3.2.1 污染源參數(shù)的傳感器與污染源相對位置的影響

圖8、圖9 給出了當污染源分別位于A～E 時通過貝葉斯推斷得到的污染源位置和強度的后驗概率。在污染源位置的后驗概率云圖中，星號處為污染源的真實位置，空心圓為傳感器；在污染源強度的后驗概率柱狀圖中，虛線表示污染源的真實強度。污染源A位置的后驗概率分布在x方向呈較寬的分布，后驗概率均值為-1.29H，與真實值之間的誤差為0.54H，后驗概率極值點為-0.75H，與真實值一致，但其后驗概率標準差高達0.40H，表明污染源位置x坐標的反演結果具有較大的不確定性；污染源A的強度后驗概率均值和極值點分別為0.72qs和0.62qs，反演結果與真實值相比偏小。雖然污染源B位置的反演誤差和不確定性與污染源A相比有所減少，但是在其強度反演結果中，后驗概率分布較為平坦，后驗概率均值為1.87qs，與真實值之間的誤差高達87%，后驗概率標準差為0.80qs，表明污染源B強度的預測可靠性較低。

圖8 污染源A～E的位置概率云圖Fig.8 Probability contour of location for pollution source A-E

圖9 污染源A～E的強度概率密度Fig.9 Probability density of strength for pollution source A-E

污染源A 和污染源B 的反演誤差較大，是因為以A 和B 為污染源做污染物擴散數(shù)值模擬時，傳感器位置處的測量濃度和模擬濃度之間存在較大的誤差，所以當數(shù)值模擬不能準確反映真實的測量結果時，會導致反演出現(xiàn)較大的偏差。同時，在反演結果中可以發(fā)現(xiàn)，污染源A 和污染源B 的源參數(shù)后驗概率標準差也大于其他位置的污染源。這是因為這兩個污染源與傳感器之間的距離較遠。由似然概率式（3）可知，當污染源參數(shù)xs，、ys和qs能正確描述真實污染源時，似然概率達到最大值，即概率極值。在計算距傳感器較遠的污染源所產生的污染物在傳感器位置的模擬濃度c（xs，ys，qs）時，模擬濃度值對于源參數(shù)的變化不敏感，即和較小，后驗概率在概率極值點附近的變化不明顯，導致后驗概率標準差較大，因此污染源A 和污染源B 的反演結果具有較大的不確定性。為了減少反演結果的不確定性，可以根據(jù)歷史信息（曾經發(fā)生過泄露的位置）和經驗設計傳感器布置方式。

在計算污染源C～E 產生的污染物擴散至傳感器位置的濃度值時，測量濃度和模擬濃度的計算結果吻合較好，因此源參數(shù)反演的誤差較小。同時，傳感器的濃度值對于污染源C、污染源D 和污染源E的參數(shù)變化靈敏度較高，導致了源參數(shù)的后驗概率標準差較低，反演的不確定性較小。

3.2.2 利用不同擴散階段測量數(shù)據(jù)對反演結果的影響

如圖6 所示，以t=660H/uH為分界線將污染物擴散過程分為兩個階段。第一個階段為發(fā)展階段，包含了污染物濃度從零濃度達到穩(wěn)定狀態(tài)的發(fā)展過程；第二個階段為穩(wěn)定階段，完全由周期性波動的測量數(shù)據(jù)構成。下面分別利用污染源在C處時傳感器位置處發(fā)展階段、穩(wěn)定階段和所有時刻的測量數(shù)據(jù)對污染源參數(shù)進行反演，反演結果分別如圖10、圖11所示。相比利用穩(wěn)定階段的測量數(shù)據(jù)的污染源位置反演結果，利用發(fā)展階段的測量數(shù)據(jù)可以降低位置誤差和減小后驗概率標準差，但穩(wěn)定階段數(shù)據(jù)提升的效果并不顯著。圖表中的結果也反映出利用所有時刻的測量數(shù)據(jù)時，其反演結果的精確程度提高非常有限。

圖10 利用不同階段測量數(shù)據(jù)的污染源A～E的位置概率云圖Fig.10 Probability contour of location for pollution source A-E by using measurement from different stages

圖11 利用穩(wěn)定階段測量數(shù)據(jù)的污染源A～E的強度概率密度Fig.11 Probability density of strength for pollution source A-E by using measurement from different stages

4 結論

本文基于貝葉斯推斷對時變流場下的污染源參數(shù)進行了反演，比較了不同位置污染源的參數(shù)反演結果，利用不同擴散階段的測量數(shù)據(jù)探討了其對反演結果的影響，得到了以下結論：

（1）模擬濃度與測量濃度之間的誤差是影響反演精確性的主要因素，當數(shù)值模擬能準確地反映污染物擴散時，污染源反演的準確性能夠得到極大的提高。

（2）當污染源與傳感器距離較近時，污染源參數(shù)的后驗概率分布比較集中，反演結果的不確定性較小，真實的污染源參數(shù)可能出現(xiàn)在較窄的區(qū)間，污染源所需搜索的范圍較??；當污染源與傳感器距離較遠時，污染源參數(shù)的后驗概率分布較寬，反演結果的不確定性較大，真實的污染源參數(shù)可能出現(xiàn)在較寬的區(qū)間，需要較大范圍搜索污染源。

（3）不同階段的污染物擴散測量數(shù)據(jù)對反演結果不會造成顯著的影響。

本研究的不足之處在于污染源參數(shù)反演是在已知流場的前提下進行的，流場的未知性會影響模擬濃度的計算，從而可能造成反演不準確。如何減少流場未知性對反演的影響，需要進一步深入研究。