文 秀 秀

(四川大學(xué)水力學(xué)與山區(qū)河流開發(fā)保護(hù)國家重點實驗室，成都 610065)

0 引 言

自20世紀(jì)80年代起，許多學(xué)者對河道突發(fā)事故污染源反問題進(jìn)行了研究。金忠青等利用脈沖譜-優(yōu)化方法解決了一維以及二維情況的對流-擴散方程源項識別以及邊界條件控制反問題[3,4]；WANGER B J[5]將水流和污染物輸運模擬與非線性極大似然估計相結(jié)合，識別處置污染源及過程參數(shù)；閔濤[6,7]等給出了利用遺傳算法求解對流-擴散方程源項識別反問題的一種新方法；SINGH R M[8]等將遺傳算法與仿真模型相結(jié)合，并用于污染源識別;朱嵩[9]等利用貝葉斯推理建立了二維水源對流擴散方程參數(shù)估計的數(shù)學(xué)模型,通過貝葉斯定理,獲得了模型參數(shù)的后驗分布,從而獲得反問題的解。

本文通過寬淺河道瞬時點源二維水質(zhì)模擬分析的數(shù)學(xué)模型的解析解的變形，在降解系數(shù)已知的情況下，根據(jù)最小二乘原理求出待求參數(shù)的綜合值，然后根據(jù)參數(shù)間的內(nèi)在關(guān)系，聯(lián)立方程得出待求參數(shù)。這種方法易于理解，所需數(shù)據(jù)簡單，能夠同時反演得到河道離散系數(shù)、河道平均流速、污染源初始位置、污染源排放時間、污染源強度等參數(shù)。

1 方 法

由相關(guān)資料[10]可以查閱二唯瞬時點源排放的解析解為：

(1)

式中：C(x,y,t)代表在t=t0時刻，在平面二維坐標(biāo)(x0,y0)處瞬時有質(zhì)量為M的污染物進(jìn)入水體，則因此產(chǎn)生的不同時間不同位置的濃度，kg/m3；ux和uy分別代表水流在縱向、橫向上的流速分量，m/min；M為瞬時點源派發(fā)的污染物數(shù)量，kg；K為污染物的降解系數(shù)，min；Dx、Dy分別代表污染物的縱向、橫向離散系數(shù)，m2/min。

我們假設(shè)事故是發(fā)生在水域為河流寬深比足夠大的所謂寬淺河道，并且假定河道之前未被污染，即污染物的背景濃度值為0。若初始監(jiān)測時刻不為0，設(shè)為t0，則解析解變?yōu)椋?/p>

C(x,y,t)=

(2)

對式(2)兩邊取對數(shù)，因為實際中uy一般取0，所以ux為水流平均速度，對上式進(jìn)行變形：

如今,小學(xué)數(shù)學(xué)的主要教學(xué)目的仍是提高升學(xué)率。任課教師也把教學(xué)的主要精力放在培養(yǎng)學(xué)生的解題能力，從而在考試中取得較好的成績。很多教師在課上忽視了數(shù)學(xué)文化的滲透，注重在有限的時間內(nèi)給學(xué)生灌輸更多的知識，做更多的練習(xí)題。因此，加強小學(xué)課堂數(shù)學(xué)文化的滲透應(yīng)該成為小學(xué)教學(xué)課堂的重要內(nèi)容之一。我們對數(shù)學(xué)的教育應(yīng)該為學(xué)生數(shù)學(xué)文化素養(yǎng)的養(yǎng)成提供有效途徑,培養(yǎng)學(xué)生解決問題的思維模式,為學(xué)生解決問題提供方法。

ax2+bx+cy2+dy+e

(3)

其中：

(4)

如果我們間隔一定時間對濃度進(jìn)行再次采樣，可以得到另一組a、b、c、d、e這5個表征著未知參數(shù)綜合信息的擬合參數(shù)。設(shè)第1組擬合參數(shù)的下標(biāo)為1，第2組擬合參數(shù)的下標(biāo)為2。

我們根據(jù)2個不同時刻的a、b、c、d、e可以得到：

(5)

經(jīng)由以上分析可以發(fā)現(xiàn)，對二維寬淺河道瞬時點源的參數(shù)估計問題就轉(zhuǎn)化為了求解2個不同時刻擬合參數(shù)a、b、c、d、e的問題。

根據(jù)最小二乘法進(jìn)行參數(shù)擬合的步驟，首先要選出相應(yīng)于式(3)的基函數(shù)，本文中為線性多項式系：{1,x,y,x2,y2}。

若假設(shè)有n個不同位置的監(jiān)測點,那么根據(jù)這些監(jiān)測點的坐標(biāo)，我們可以得到一個n行5列的矩陣，即矩陣A：

(6)

式中：x，y分別為這n個監(jiān)測點各自對應(yīng)的橫縱坐標(biāo)。

同時，為與式(3)中所作的對數(shù)變形相對應(yīng)，對某一采樣時刻下這n個監(jiān)測點的濃度數(shù)據(jù)進(jìn)行對數(shù)處理，得到矩陣B：

(7)

式中：C1，C2，C3，…，Cn為各個監(jiān)測點的濃度值。

根據(jù)最小二乘法進(jìn)行曲線擬合的步驟，此時問題轉(zhuǎn)化為求解線性方程組：

(ATA)S=ATB

(8)

式中：S=(S1,S2,S3,S4,S5)T，為待求的擬合參數(shù)。

通過Doolitle分解將ATA分解為一個上三角矩陣L,一個下三角矩陣U,進(jìn)行求解，結(jié)果{S1、S2、S3、S4、S5}分別為式(4)中的e、b、d、a、c的值。

2 算例描述

假定有一寬淺河道研究區(qū)域，平均水深h=1.0 m,縱向流速ux=17 m/min,橫向流速uy=0,離散系數(shù)Dx=20 m2/min,Dy=30 m2/min。在x0=-2 100 m,y0=25 m處發(fā)生瞬時點源污染事故，在t0=-120 min時刻該點源釋放污染物總量M=5 000 kg,降解系數(shù)K=0.001 min-1。

由公式(2)計算得到第i個監(jiān)測點在j時刻的污染物濃度的理論值Cobs(xi,yi,tj)(i=1,2,…,n;j=1,2,…,k)，為了使結(jié)果更貼近實際情況，考慮監(jiān)測誤差，對理論值附加一個高斯數(shù)據(jù)噪聲[11]，即:

Cobs(xi,yi,tj)=Cobs(xi,yi,tj)+Ewi,jCobs(xi,yi,tj)

(9)

式中：wi,j為在(-1,1)區(qū)間內(nèi)的標(biāo)準(zhǔn)正態(tài)分布隨機數(shù)；E為監(jiān)測誤差水平。

在此算例中，我們選取2個不同時刻：t1=10 min,t2=20 min，每個時刻選取16個監(jiān)測點，即n=16,k=2，詳細(xì)監(jiān)測點坐標(biāo)數(shù)據(jù)見表1。

表1 監(jiān)測點坐標(biāo) mTab.1 Monitoring point coordinates

3 算例結(jié)果與影響因素分析

3.1 濃度監(jiān)測誤差的影響

選擇不同的監(jiān)測誤差水平，即不同的E值， 本次算例中取E為5%、10%、20%，得到的不同監(jiān)測誤差水平下參數(shù)反演結(jié)果見表2。

表2 不同監(jiān)測誤差下污染源的參數(shù)反演結(jié)果Tab.2 Parameter inversion results for pollution sources under different monitoring errors

通過不同監(jiān)測誤差下污染物參數(shù)反演的絕對相對誤差(Absolute Relative Error)ARE來衡量反演的結(jié)果,詳見圖1。相關(guān)計算公式如下:

(10)

式中：rE為不同監(jiān)測誤差水平下的參數(shù)反演值；r0為參數(shù)的給定值。

圖1 不同E值下各參數(shù)反演結(jié)果的ARE值Fig.1 ARE for the inversion result of each parameter under different monitoring error levels

由表2可以看出：監(jiān)測誤差與參數(shù)的反演結(jié)果總的趨勢表現(xiàn)為監(jiān)測誤差越大，參數(shù)反演值A(chǔ)RE就越大，也就是反演結(jié)果越不精確。但是即使監(jiān)測誤差達(dá)到20%，相對誤差也在10%以內(nèi)。并且由圖1可以看出，在同樣的誤差水平下，Dy、M的反演結(jié)果最差，平均流速ux的反演結(jié)果最好。這是由于參數(shù)Dy、M是在其余參數(shù)的基礎(chǔ)上計算的，相對于其他的參數(shù)并不獨立，所以Dy、M的反演結(jié)果受到誤差累積的影響，相對誤差較大。

為了進(jìn)一步對監(jiān)測誤差帶來的影響進(jìn)行分析，利用公式(2)和給定的各參數(shù)取值得到解析解濃度值，并將其在平面坐標(biāo)系中畫出，得到不同位置在t1=10 min時的濃度值。分別將E=5%、10%、20%時得到的反演參數(shù)帶入公式(2)，得到相應(yīng)的濃度值，并與參數(shù)給定值得到的濃度值做差，將結(jié)果在平面直角坐標(biāo)系中表現(xiàn)出來，見圖2、圖3。

圖2 E=0時濃度分布Fig.2 Concentration distribution under the condition of E=0

圖3 不同E值下濃度值與E=0時濃度值差值分布Fig.3 Concentration difference distribution between different value of E and E=0

可以看出，不同E值下，根據(jù)反演參數(shù)求得的濃度值均與給定參數(shù)下求得的濃度值存在差距，且表現(xiàn)為濃度值大的地方偏差大，濃度值小的地方偏差?。徊⑶矣捎趨?shù)反演存在誤差，濃度值偏差存在偏態(tài)分布的特性，而不是嚴(yán)格符合正態(tài)分布特性。

3.2 監(jiān)測點數(shù)量對反演結(jié)果的影響

從理論上講，隨著監(jiān)測點數(shù)量的增加，反演結(jié)果也會更加精確。為了探究這一點，設(shè)計了一系列算例進(jìn)行研究。算例中各參數(shù)給定值與上述條件相同，只是將監(jiān)測點數(shù)量分別減少到12組以及增加到24組，監(jiān)測誤差設(shè)定為E=5%。各參數(shù)的反演結(jié)果詳見表3。

采用上文提及的ARE(絕對相對誤差)對以上的反演結(jié)果進(jìn)行分析，將分析比較結(jié)果直觀地顯示在圖4中。

由圖4反演結(jié)果對比分析可得：減小測點數(shù)量，反演參數(shù)的相對誤差較之前發(fā)生了較為明顯的增大，而增加測點數(shù)量，反演參數(shù)的精度略有提高，但是提高不明顯，這主要是由于本身使用測點數(shù)量已經(jīng)足夠反演濃度分布信息，因此再增加監(jiān)測點數(shù)目對反演結(jié)果的精確性影響不大，但是如果減少測點數(shù)目，就造成了測點數(shù)量較少，在具有監(jiān)測誤差時不能真實反映濃度分布信息。同時，我們通過橫向?qū)Ρ?，不難發(fā)現(xiàn)，依然是Dy、M的反演結(jié)果最差，平均流速 的反演結(jié)果最好。可見，該方法在流速的反演上具有較高精度，而在Dy的反演上精度較低，但是也不超過10%的相對誤差。

表3 E=5%時不同測點數(shù)量下參數(shù)反演結(jié)果Tab.3 Parameter inversion results of different number of monitoring points under the condition of E=5%

圖5 不同測點數(shù)目下各參數(shù)反演結(jié)果的ARE比較Fig.5 ARE comparison of parameters inversion results under different number of monitoring points

4 結(jié)論與展望

通過以上的算例，可以得到：對于寬淺河道中的污染源參數(shù)反演問題，本文提出的基于最小二乘法的采用同時刻不同空間點的濃度數(shù)據(jù)的反演方法是可行的，該方法能夠精確地反演出各個源項參數(shù)(釋放時間、釋放位置、釋放量)以及水流參數(shù)(流速、縱向離散系數(shù))；該方法對監(jiān)測誤差具有一定的敏感性，誤差越大，反演結(jié)果越差，但是即使監(jiān)測誤差為20%時，反演結(jié)果相對誤差也不超過10%，具有實際應(yīng)用的價值；監(jiān)測點數(shù)量對參數(shù)反演的結(jié)果具有一定影響，在合理范圍內(nèi)，隨著監(jiān)測點數(shù)量的增加，參數(shù)反演的精度也在提高，但是監(jiān)測點數(shù)量增加所帶來的參數(shù)反演結(jié)果的精確性提升是逐漸減弱的。本研究由于缺乏實測資料，故采用假想算例進(jìn)行驗證，在下一步的工作中將搜集實測資料對方法進(jìn)行驗證并作出改進(jìn)。