金文龍，彭 輝，高俊杰，潘曼曼，韓龍喜

(1.河海大學環(huán)境學院，江蘇南京 210098;2.水利部水文局，北京 100053)

隨著我國社會經(jīng)濟的迅速發(fā)展，突發(fā)性水污染事故頻頻發(fā)生，由于突發(fā)水污染事故的不確定性及發(fā)現(xiàn)時間的滯后性，一般情況下很難及時監(jiān)測到完整的污染源排放歷史，給事故發(fā)生后的預警預報與處置工作帶來困難。因此，如何在突發(fā)性水污染事故發(fā)生后，根據(jù)下游污染物濃度分布等信息反演完整的污染源排放歷史具有重要意義。

在環(huán)境水力學反問題領域，地表水污染源的反演問題通常被提為源項識別反問題，其主要任務是根據(jù)環(huán)境水力學測量數(shù)據(jù)，來反推污染源的信息，如污染源的強度和位置等。目前相關方面的研究尚不成熟，Akcelika 等［1］采用變分有限元方法求解了對流擴散輸運的源項反演問題;Wang 等［2］提出了多孔介質流動中的源項識別的馬爾科夫隨機場模型;Ling 等［3］從理論上研究了二維對流擴散方程的源項位置識別問題;閔濤等［4］將對流-擴散方程源項識別反問題轉化為優(yōu)化問題，并用遺傳算法求解;韓龍喜［5］以污水處理費用最小為目標函數(shù)，采用遺傳算法求解了一維水質模擬帶約束條件的污染物排放量控制反問題;并構造了河網(wǎng)地區(qū)帶約束條件的污染源控制反問題［6］，并采用約束最優(yōu)化算法中的簡約梯度法成功求解。王澤文等［7］研究了流域中單個點污染源識別問題的唯一性、穩(wěn)定性和反演算法。

地表水污染源的反演問題也可被提為邊界條件反問題，金忠青等［8］分別應用脈沖譜- 優(yōu)化法和Green 函數(shù)法求解了一維、二維對流-擴散方程邊界條件控制反問題，基本解決了恒定水體中上游污染控制問題;韓龍喜等［9］首次給出水網(wǎng)地區(qū)無約束的水質邊界濃度反問題的提法，并采用局部基本解展開算法進行穩(wěn)態(tài)條件下邊界濃度反演;李子［10］將識別突發(fā)性污染源以及模擬污染物濃度時空變化過程提為邊界條件控制反問題，并對以GST-MQ(global space time-multi quadric)配點法為基礎的反演方法進行了研究。

總體而言，關于污染源反演問題的研究主要集中在污染源排放位置或排放區(qū)域已知情況下的排污(污染源)的控制問題以及對邊界上濃度變化過程的反演，很少涉及直接對污染源源強變化過程的反演。因此，筆者在污染源位置已知的情況下，根據(jù)事故下游某斷面一組污染物濃度的時間序列觀測資料，提出了中小型河道一維突發(fā)性水污染事故源強變化過程的反問題。在反問題的求解方法上，筆者運用了徑向基配點法，該方法作為一種新興的無網(wǎng)格方法具有不依賴于控制方程的形式、multi-quadric(MQ)函數(shù)的插值精度高、求解反問題與正問題計算過程相同、不需要迭代優(yōu)化過程以及對測量數(shù)據(jù)有較強的適應性等獨特優(yōu)勢。

1 問題的數(shù)學描述

1.1 中小型河道一維水質模型

假定研究水域是長度為L 的中小型河道，污染源位于研究水域上邊界(x =0)處，河道寬度和深度均遠小于河段的長度，當污染事故發(fā)生時，污染物在對流、擴散過程的作用下在水域［0，L］內(nèi)隨水體遷移擴散，對于中小型河道，河段一般較為順直，因此污染物的斷面平均濃度C 沿河段的分布可用源項為零的一維對流-擴散方程來描述:

式中:D 為擴散系數(shù);t 為時間;x 為空間坐標;u 為流速;K 為污染物降解系數(shù)。

其初始條件為:

上邊界條件用源強的變化過程S(t)來表示:

式中:Q 為流量;C0為上游邊界(x =0)處的背景濃度。

下邊界為Neumann 邊界條件:

當D、u、K、C0(x)、C0、S(t)已知時，式(1)～(4)構成河流一維水質模擬分析的適定問題即正問題，可采用常規(guī)的計算方法預測水域中任何時刻的污染物濃度分布C(x，t)。為研究簡便計，不妨構造一個污染源，根據(jù)突發(fā)性水污染事故的特點，假定污染源為分布于某一時段的非連續(xù)源，其源強變化過程的數(shù)學描述為:

其中H 為階躍函數(shù):

式中:As為源強排放峰值;Ts為排放的持續(xù)時間。

假定初始條件C0(x)=0，則定解問題式(1)～(4)的解析解為:

1.2 反問題的提法

由于突發(fā)水污染事故的不確定性，一般情況下很難及時監(jiān)測到完整的源強變化過程S(t)，從而導致邊界條件C(0，t)的缺失，對求解正問題造成困難。為初步判斷污染源的真實情況，一類有效的方法就是在污染源下游某斷面x =l 處設置觀測點，通過采樣觀測，獲得一組污染物濃度的時間序列觀測資料:

式中，W 為觀測數(shù)據(jù)的個數(shù)。

通過這組觀測數(shù)據(jù)來反演源強的變化過程，這就構成了一類源強變化過程的反問題。這類反問題的實質是給出方程的部分解來確定控制方程的邊界條件，根據(jù)環(huán)境水力學反問題的觀點，即為邊界條件控制反問題。應當指出的是，從物理角度出發(fā)，邊界條件和污染源項是兩個不同的概念，但在水質模型中，當污染源位于邊界時邊界條件反問題與源項反問題常是可以相互轉換的，兩者從數(shù)學角度看處理方法是一致的［11］。

2 反問題的求解方法

2.1 最小二乘GST-MQ 配點法簡介

徑向基函數(shù)配點法是求解偏微分方程的一種配點型無網(wǎng)格法，它的基本思想是把未知函數(shù)用一組徑向基函數(shù)的線性組合來表示，然后再令域內(nèi)的各節(jié)點滿足控制方程，邊界上各節(jié)點滿足一致邊界條件［12］。最小二乘-徑向基函數(shù)配點法的基本思想是用徑向基函數(shù)逼近未知函數(shù)后，令控制方程在包括區(qū)域內(nèi)和邊界上的所有節(jié)點上滿足，而初始條件、邊界條件等則在其各自相應的節(jié)點上滿足，最后利用最小二乘法求解超定線性方程組［13-14］。相對于傳統(tǒng)徑向基配點法，最小二乘-徑向基配點法可以松弛邊界上各節(jié)點，滿足邊界條件這一要求，從而改善線性方程組中系數(shù)矩陣的病態(tài)條件，提高解的精度。

MQ 函數(shù)是由Hardy 于1968 年提出來的一種徑向基函數(shù)［15］。對于一維空間，MQ 函數(shù)的形式為

式中:c 為形狀參數(shù);xk為空間基點坐標。

逆MQ 函數(shù)的形式為

由于MQ 函數(shù)和逆MQ 函數(shù)都只是空間距離的函數(shù)，只能用于求解穩(wěn)態(tài)問題，而在水質模擬分析中，污染物的濃度變化是個非穩(wěn)態(tài)的過程，因此引入全時空的概念，將時間維看作偽空間，進而將時間變量直接插入基函數(shù)形式中構造新的徑向基函數(shù)。以逆MQ 函數(shù)為例進行構造，新構造的GST-MQ 徑向基函數(shù)形式為:

其中

式中:ω 為引入的時空比例因子，以平衡時間和空間的屬性差異;tk為基點時間;Δx 為空間上的配點距離;Δt 為時間上的配點距離;c0為形狀參數(shù)相對值。

2.2 計算步驟

a. 將空間區(qū)域［0，L］和計算時間范圍［0，T］分別用N 和M 個節(jié)點離散，則污染物濃度C(x，t)可由一組GST-MQ 徑向基函數(shù)近似為:

其中

式中:αij為待定系數(shù);N、M 分別為空間區(qū)域［0，L］和計算時間范圍［0，L］上的配點總數(shù)。

b. 將式(11)代入式(1)，并在時空區(qū)域及邊界上的所有節(jié)點上滿足，可得一個關于αij、有N ×M個方程、N×M 個未知數(shù)的方程組:

c. 將式(11)代入式(2)，并在空間區(qū)域［0，L］的所有配點及t =0 上滿足，可得一個關于αij、有N個方程、N×M 個未知數(shù)的方程組:

d. 將式(11)代入式(4)，并在x =L 及時間范圍［0，T］的所有配點上滿足，可得一個關于αij、有M個方程、N×M 個未知數(shù)的方程組:

e. 將式(11)代入式(7)，并在x=l 及觀測時間{tq;q=1，2，…，W}上滿足，可得一個關于αij、有W個方程、N×M 個未知數(shù)的方程組:

則式(10)～(13)構成了一個關于αij，含有N×M 個未知數(shù)，N×M +N +M +W 個方程的超定線性方程組，即最小二乘問題:min{‖Aα -b‖2}。用奇異值分解法解這個最小二乘問題可得系數(shù)αij，再利用式(3)可反演得源強變化過程S(t)。

3 數(shù)值算例

假定事故水域范圍L =5 000 m，計算時間T =10000 s，水體流動速度u=0.5 m/s，流量Q=20 m3/s，擴散系數(shù) D = 1.0 m2/s，污染物降解系數(shù)K=0.000 1s-1。假定源強峰值As=200 g/s，污染物排放時間Ts=5000s，觀測斷面位于l=1000m 處，觀測起訖時間為事發(fā)后400 ～10000s 時間段，觀測頻次為每隔400 s 觀測一次，對應的觀測數(shù)據(jù)個數(shù)W=25 個。

首先由正問題的解析解式(6)給出一組觀測數(shù)據(jù)的理論值，再由式(16)生成觀測誤差并附加到理論值，作為實際的觀測數(shù)據(jù)值:

式中:h(tq)為構造的觀測數(shù)據(jù)的實際值;h'(tq)為觀測數(shù)據(jù)的理論值;e 為誤差水平;ωq為一組符合標準正態(tài)分布的隨機數(shù)。

用式(7)和式(16)分別計算觀測數(shù)據(jù)的理論值h'(tq)和實際值h(tq)，誤差水平e 取0.02，部分計算結果見表1。

表1 計算的觀測數(shù)據(jù)理論值和實際值的部分結果

最小二乘GST-MQ 配點法的計算參數(shù)有:空間上的配點總數(shù)N、時間上的配點總數(shù)M、形狀參數(shù)c0、和時空比例因子ω。對于空間和時間上的配點總數(shù)N、M，一般情況下，N 和M 越大，計算結果越精確，但計算機所需內(nèi)存越多，計算時間越長，本文在計算時取N =51、M =26。對于形狀參數(shù)c0和時空比例因子ω，反演誤差對這兩個參數(shù)的取值較為敏感，本文經(jīng)試算分析，在計算時取c0=1.5，ω=1.4。

根據(jù)表1 中的觀測數(shù)據(jù)的實際值，使用最小二乘GST-MQ 配點法程序對源強變化過程S(t)進行反演。另外，為減小計算結果的波動，提高反演精度，對源強增加約束條件，使源強不小于0，即S(t)≥0。經(jīng)計算，反演計算結果與真實值的比較見圖1。

由圖1 可知:在t =0 s 和t =Ts=5 000 s 附近，S(t)的反演解與真實值偏差較大，這是由于t =0 s時刻污染源排放突然開始，t=5 000 s 時刻污染源排放突然停止造成的。在數(shù)學上，GST-MQ 徑向基函數(shù)連續(xù)可導，因此用連續(xù)光滑的反演解描述真實排放過程在t=0 s 和t =5 000 s 處的間斷性。在其他時段，S(t)的反演解在真實值附近波動，大致反映了污染源的排放過程?？傮w而言，S(t)的反演解在大多數(shù)時段比較接近污染源的真實情況，反演結果較為準確。

4 結 語

根據(jù)中小型河道特點以及突發(fā)性水污染事故的特征，用污染源下游某斷面污染物濃度的時間序列觀測資料作為附加條件，提出了突發(fā)性水污染事故污染源強變化過程反問題，并用最小二乘GST-MQ配點法對反問題進行了求解。數(shù)值算例表明，反問題的計算結果準確合理。

圖1 S(t)的反演解和真實值比較

［1］AKCELIKA V，BIROS G，GHATTAS O，et al.A variational finite element method for source inversion for convectivediffusive transport［J］. Finite Elements in Analysis and Design，2003，39:683-705.

［2］WANG J，ZABARAS N. Using bayesian statistics in the estimation of heat source in radiation［J］. International Journal of Heat and Mass Transfer，2005，48:15-29.

［3］LING L，YAMAMOTO M，HON Y C，et al.Identification of source locations in two-dimensional heat equations ［J］.Inverse Problems，2006，22:1289-1305.

［4］閔濤，周孝德，張世梅，等.對流-擴散方程源項識別反問題的遺傳算法［J］. 水動力學研究與進展:A 輯，2004，19 (4):520-524. (MIN Tao，ZHOU Xiaode，ZHANG Shimei，et al. Genetic algorithm to an inverse problem of source term identification for convectiondiffusion equation ［J ］.Chinese Journal of Hydrodynamics，2004，19(4):520-524. (in Chinese))

［5］韓龍喜.河道一維污染源控制反問題［J］.水科學進展，2001，12(1):39-44. (HAN Longxi. Inverse problem on amount of pollutant into natural channel［J］. Advances in Water Science，2001，12(1):39-44. (in Chinese))

［6］韓龍喜，朱黨生.河網(wǎng)地區(qū)水環(huán)境規(guī)劃中的污染源控制方法［J］. 水利學報，2001(10):28-31. (HAN Longxi，ZHU Dangsheng. The pollutant source control method for network water environment management ［J］. Journal of Hydraulic Engineering， 2001 (10 ): 28-31. (in Chinese))

［7］王澤文，邱淑芳. 一類流域點污染源識別的穩(wěn)定性與數(shù)值模擬［J］. 水動力學研究與進展:A 輯，2008，23(4):364-371. (WANG Zewen，QIU Shufang. Stability and numerical simulation of pollution point source identification in a watershed ［J］. Chinese Journal of Hydrodynamics，2008，23(4):364-371. (in Chinese))

［8］金忠青，陳夕慶.用脈沖譜-優(yōu)化法求解對流-擴散方程邊界條件控制反問題［J］.河海大學學報:自然科學版，1991，19 (1):1-8. (JIN Zhongqing，CHEN Xiqing.Numerical solution to an inverse problem of boundary condition control for convection-diffusion equations by PST-Optimization ［J］. Journal of Hohai University:Natural Sciences，1991，19(1):1-8. (in Chinese))

［9］韓龍喜，蔣莉華，朱黨生.組合單元水質模型中的邊界條件及污染源項反問題［J］.河海大學學報:自然科學版，2001，29(5):23-26. (HAN Longxi，JIANG Lihua，ZHU Dangsheng. Inverse problem on boundary condition and pollutant source in water quality model of river network ［J］. Journal of Hohai University:Natural Sciences，2001，29(5):23-26. (in Chinese))

［10］李子.基于GST-MQ 配點法的突發(fā)水污染事故反演模型研究［D］.北京:清華大學，2010.

［11］劉曉東，姚琪，薛紅琴，等.環(huán)境水力學反問題研究進展［J］. 水科學進展，2009，20 (6):885-890. (LIU Xiaodong，YAO Qi，XUE Hongqin，et al. Advance in inverse problems of environmental hydraulics ［J］.Advances in Water Science，2009，20(6):885-890. (in Chinese))

［12］劉芳.徑向基函數(shù)配點法的應用及其誤差估計［D］.蘇州:蘇州大學，2008.

［13］ROCCA A L，POWER H. A double boundary collocation hermitian approach for the solution of steady state convection-diffusion problems ［J］. Computers and Mathematics with Applications，2008，55:1950-1960.

［14］CHEN W，TANAKA M. New insights into boundary-only and domain-type RBF methods［J］. International Journal of Nonlinear Sciences and Numerical Simulation，2000，1(3):145-152.

［15］魏義坤.徑向基函數(shù)插值法解偏微分方程及計算滲流問題［D］.成都:成都理工大學，2009.