胡 煜,張文俊,任華堂,夏建新

(中央民族大學環(huán)境科學系,北京 100081)

隨著工業(yè)化進程的不斷加快,近年來我國進入了環(huán)境高風險時期,環(huán)境污染事件以其突發(fā)性和不確定性成為環(huán)境保護的重大威脅。在水環(huán)境事故應急處理過程中,截斷污染源、評估污染事故的環(huán)境影響以及污染事故的責任認定是核心內容[1],而這三者都需要根據(jù)監(jiān)測數(shù)據(jù)確定污染源泄漏的位置和所泄漏的污染負荷量[2]。因此,基于污染事故發(fā)生中的監(jiān)測數(shù)據(jù)進行反演確定初始時刻污染源的分布特性具有重大意義。

利用不同時刻、不同位置的監(jiān)測數(shù)據(jù)確定污染源相當于求解對流-擴散方程的逆過程,該問題在一般意義上是不適定的，且極易導致數(shù)值計算發(fā)散[3]。雖然相關研究人員對正問題進行了大量研究并取得一定進展,但目前針對逆過程的研究相對較少。通常需要對源項做出必要的先驗假設才能使反問題的求解具有一定的可行性。如,Ling等[4]從數(shù)學理論的角度證明了二維熱傳導方程在不考慮測量誤差條件下未知源項位置具有唯一性,并且在基于未知源函數(shù)是若干已知源函數(shù)之和的先驗假設下證明了3個測量點數(shù)據(jù)對于確定未知源函數(shù)的必要性并進行了數(shù)值驗證；王澤文等[5]在對污染源及測量點的先驗假設下,采用截斷奇異值分解正則化方法求解離散的病態(tài)方程組,確定非零邊界條件下解的唯一性,并進一步獲得源項識別的局部穩(wěn)定結果；而Bruckner等[6]針對波動方程的初邊值問題,在監(jiān)測指標向量處于有限范數(shù)的必要條件下,證明了點源識別的唯一性和穩(wěn)定性。

對于反問題的求解主要有3類方法:①概率統(tǒng)計方法,如朱嵩等[7]基于貝葉斯推理建立了水動力-水質耦合數(shù)學模型的點源強度與位置的聯(lián)合識別方法,通過馬爾科夫鏈蒙特卡羅后驗抽樣獲得了污染源位置和強度的后驗概率分布和估計量。②遺傳算法或人工神經網(wǎng)絡方法,如閔濤等[8]應用遺傳算法將源項反問題和參數(shù)反問題轉化為優(yōu)化問題,研究一維對流-擴散方程源項識別反問題,辛小康等[9]將遺傳算法與數(shù)學分析相結合建立了水污染事故污染源識別模型。這兩類方法均存在物理意義不夠明確,且所需樣本量比較大的問題。③伴隨方法,該方法源于氣象和海洋預報領域的數(shù)據(jù)同化研究[10],Bennet等[11]最早將其應用于物理海洋學領域,Tziperman等[12]相繼利用該方法反演水位和海洋初始狀態(tài),國內學者馬繼瑞等[13]利用伴隨方法反演了水位的初始場。

近年來,伴隨方法逐步被引入到環(huán)境領域解決污染源反演問題,但是多數(shù)相關研究集中于邊界初值或相關參數(shù)的反演,對于濃度場反演甚少涉及。如谷藝等[14]利用伴隨方法根據(jù)水位觀測資料反演渤海二維非線性潮汐模型的底摩擦系數(shù);Akcelik等[15]采用Tikhonov正則化和總變差正則化[16]的方法求解了對流-擴散輸運系統(tǒng)的邊界初值反演問題;吳自庫等[17]在利用伴隨同化方法求解一維對流-擴散方程逆過程反問題中,從物理上對代價函數(shù)增加約束以克服不適定性與計算不穩(wěn)定性,實現(xiàn)了污染源初始條件的反演。在實際環(huán)境污染事故中,污染物隨時間在空間場中迅速遷移擴散,空間位置的不確定性極大地增加了污染源識別的難度[18]。采用伴隨方法實現(xiàn)空間場中任意位置的污染物濃度分布時間序列反演,從而獲取污染源的污染特性信息對于對流-擴散方程反問題研究具有重要意義。

本文基于伴隨方法,將對流-擴散方程逆過程污染源反演轉化為模型初始場最優(yōu)化問題,建立目標泛函,差分求解原方程與伴隨方程,確定目標泛函的下降梯度,迭代求解使得泛函達極小值,實現(xiàn)了污染物遷移轉化逆過程反問題的求解。

1 對流-擴散方程逆過程源項反演的伴隨方法

1.1 對流-擴散方程反向求解過程

污染源反演的實質是確定一個最可能的污染物初始空間分布,從而得到污染源的時空特性信息?；诖?首先建立目標泛函J:

(1)

式中:C為污染物質量濃度的計算值;Cobs為污染物質量濃度的觀測值;W為權重函數(shù),在觀測點有監(jiān)測資料時,其值取1,反之則取0。

污染物的遷移變化受對流-擴散方程的約束,滿足該約束條件的拉格朗日函數(shù)為

(2)

式中:λ為拉格朗日乘子;H為哈密頓算符;V為流速;K為縱向離散系數(shù);f為污染物的源(匯)項。

當泛函L在約束條件下取得極值時,即被積函數(shù)F取得極值時,目標函數(shù)J最小。目標泛函極值條件如下:

(3)

2W(C-Cobs)+ψ*(λ)=0

(4)

(5)

式中：ψ*(λ)為ψ(C)的伴隨函數(shù)。

將式(5)代入式(4)中,得到關于污染物質量濃度C的伴隨方程:

(6)

關于時間反向積分上式,得到初始時刻λi,0的全場分布,用于推算污染源相關特性。這即為反向時間積分求解對流-擴散方程的途徑。

根據(jù)式(2),可以得到目標泛函相對污染物初始濃度分布C0的梯度:

?Ω-ψ*(λ)0dΩ

(7)

其中,下標“0”表示t=T初始時刻。將式(7)寫成離散形式后,得到J的增量為

(8)

式中：T*表示時間反向,即T*=T-t。

將根據(jù)式(6)得到的初始時刻λi,0的全場分布帶入式(8),即可優(yōu)化污染物濃度分布初始場:

(9)

1.2 數(shù)值求解方法和計算過程

在一維條件下,對以上方程進行差分求解。對流-擴散方程、伴隨方程以及式(8)的離散差分格式為

(10)

(11)

(12)

式中：u為污染物遷移速度；K為擴散系數(shù)；Ciobs,n+1為n+1時刻固定空間網(wǎng)格點i所監(jiān)測到的污染物質量濃度值；Δt及Δx分別為時間步長和空間步長。

求解中,首先依據(jù)監(jiān)測資料給定初始場的猜測值。一般而言,猜測值僅影響迭代次數(shù),不影響計算結果,猜測值愈接近真值,所需迭代次數(shù)越少?；诓聹y的污染物初始場正向積分式(10),得到污染物濃度分布空間場,并由該濃度空間分布場計算目標泛函,若泛函值小于設定誤差ε,則該初始場猜測值可視為初始場實際值,否則,通過時間反向積分式(11),計算得到初始值λi,0;并進一步利用式(12)計算目標泛函的下降梯度,將其代入式(9)對污染物濃度初始場進行優(yōu)化,得到濃度初始場新的猜測值,利用優(yōu)化后初始場再次計算目標泛函,直至目標泛函滿足誤差要求,停止計算確定污染物濃度初始場。具體流程如圖1所示。

圖1 污染源反演計算流程

2 數(shù)值試驗與結果分析

2.1 計算工況與參數(shù)設置

為驗證該方法的適用性,設計兩種不同類型工況進行數(shù)值試驗。工況1與工況2針對污染事故發(fā)生后根據(jù)污染物影響范圍內的濃度空間分布資料,反演之前某時刻污染源的濃度空間分布;工況3根據(jù)固定監(jiān)測站的污染物濃度時間序列值,反演上游某位置排放的污染物時間序列。3種工況計算的污染源均為瞬時點源形成的濃度分布,污染物為保守物質且假定在河道中瞬間充分混合,污染物在水體中表現(xiàn)為對流和擴散兩種物理過程。

假設河道為一維,3種工況參數(shù)設置如下:污染物遷移速度均為u=1.0 m/s,擴散系數(shù)均為300 m2/s,空間步長均為500 m。為保證數(shù)值計算穩(wěn)定性,工況1與工況2條件下計算所采用的時間步長為1 s,工況3條件下所采用的時間步長為2 s。

工況1已知污染源為單源且污染物質量濃度分布滿足:

(13)

據(jù)此反演6.0×104s之前污染物的空間分布情況,該問題的解析解為

(14)

工況2已知雙源形成的質量濃度空間分布為

(15)

據(jù)此反演9.0×104s之前污染物的空間分布情況,該問題的解析解為

(16)

工況3給定下游固定監(jiān)測站點的質量濃度分布時間序列監(jiān)測值:

(17)

反演污染事故對于該站點上游相對距離為7.5 km的位置上游固定監(jiān)測站點的污染物質量濃度分布時間序列。該時間序列的解析解如下式:

(18)

2.2 計算結果與誤差分析

根據(jù)3種工況的數(shù)值實驗結果進行分析,探究伴隨方法對于污染源特性反演的有效性。為檢驗本文源項反演方法的精度,采用平均相對誤差M進行衡量,定義如下:

%

(19)

式中:XP為計算值;XR為監(jiān)測值。平均相對誤差反映了計算值與監(jiān)測值之間的偏離情況,其值越接近0,說明計算值與監(jiān)測值越接近,計算準確性越高。表1列出了3種工況下目標函數(shù)的變化與滿足誤差要求下的迭代步數(shù)以及上游污染物濃度分布計算值與解析解、下游污染物質量濃度分布計算值與監(jiān)測值之間的平均相對誤差。

表1 不同工況條件下的平均相對誤差

2.2.1 污染源質量濃度空間分布初始場反演

工況1條件下,上游污染源質量濃度空間分布初始場反演結果(圖2)與解析解十分接近,污染負荷計算值為20 010.375 kg/m3,解析解為19 976 kg/m3,相對誤差為0.2%。并且反演得到的初始場形成污染物濃度空間分布與監(jiān)測值相差甚小,平均相對誤差在5%之內。該結果說明本文的伴隨方法在單污染源釋放條件下質量濃度空間分布初始場反演中具有較好的可靠性與精確性。

圖2 工況1污染源質量濃度空間分布反演

從圖3可知,工況2條件下,污染源質量濃度空間分布初始場與解析解具有較好的一致性,初始場污染負荷計算值為39 910.44 kg/m3,解析解為39 928 kg/m3,相對誤差為-0.04%。同時反演得到的初始場形成污染物質量濃度空間分布與監(jiān)測值相差甚小,平均相對誤差僅為1.12%,顯示了伴隨方法在雙污染源釋放條件下質量濃度空間分布初始場反演的精確性。需要指出的是,當存在多個污染源時且反演時段過長可能會導致污染源混淆,此時需要結合其他資料確定污染源的數(shù)量。

在計算量方面,工況1和工況2均較小。計算過程中,目標泛函下降十分迅速,僅需迭代一百余步,即可滿足精度要求,顯示出伴隨方法的高效性和在反演污染源質量濃度空間分布初始場時效性上的獨特優(yōu)勢。需要指出的是,當所需反演的污染源空間質量濃度梯度過大時,容易出現(xiàn)解的不適定現(xiàn)象,計算值和真值之間會出現(xiàn)差異,這是由于現(xiàn)有的對流-擴散微分方程對濃度函數(shù)所做的光滑性假定導致,而非伴隨方法本身的錯誤。

圖3 工況2污染源質量濃度空間分布反演

2.2.2 污染物質量濃度分布時間序列反演

工況3條件下,上游站點污染物質量濃度分布時間序列反演值與解析解具有較好一致性(圖4),平均相對誤差最大不超過15%。根據(jù)反演值得到的下游染物濃度分布時間序列與監(jiān)測值基本一致,進一步說明該方法的可靠性。

圖4 工況3上游固定位置污染物質量濃度分布時間序列反演

3 結 論

a. 以污染物質量濃度監(jiān)測值與計算值之間的誤差的平方在時間、空間上的積分作為目標泛函,通過構建拉格朗日函數(shù),進一步將條件約束下的泛函極值問題轉化為無約束極值問題,通過求解泛函極值將污染源反演問題轉化為優(yōu)化問題,建立了污染源初始場的反演優(yōu)化方法。

b. 應用伴隨函數(shù)構建伴隨方程,通過反向時間積分伴隨方程并對其進行差分離散求解,避免了對污染物濃度方程直接反向積分導致的計算發(fā)散,實現(xiàn)了對流-擴散方程逆過程的數(shù)值求解,為反問題的求解提供了新的思路與方法。

c. 基于伴隨方法,利用下游監(jiān)測站點污染物濃度信息對污染源的空間初始場以及上游監(jiān)測站點的濃度時間序列進行反演,能夠反映污染物對流-擴散逆過程的物理機制,且在迭代計算過程中誤差逐漸減小,使得反演所得污染源特性與真解十分接近,平均相對誤差最大不超過15%,滿足工程應用的精度要求。