文 秀 秀

(四川大學(xué)水力學(xué)與山區(qū)河流開(kāi)發(fā)保護(hù)國(guó)家重點(diǎn)實(shí)驗(yàn)室，成都 610065)

0 引 言

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

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

1 方 法

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

(1)

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

我們假設(shè)事故是發(fā)生在水域?yàn)楹恿鲗捝畋茸銐虼蟮乃^寬淺河道，并且假定河道之前未被污染，即污染物的背景濃度值為0。若初始監(jiān)測(cè)時(shí)刻不為0，設(shè)為t0，則解析解變?yōu)椋?/p>

C(x,y,t)=

(2)

對(duì)式(2)兩邊取對(duì)數(shù)，因?yàn)閷?shí)際中uy一般取0，所以u(píng)x為水流平均速度，對(duì)上式進(jìn)行變形：

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

ax2+bx+cy2+dy+e

(3)

其中：

(4)

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

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

(5)

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

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

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

(6)

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

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

(7)

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

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

(ATA)S=ATB

(8)

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

通過(guò)Doolitle分解將ATA分解為一個(gè)上三角矩陣L,一個(gè)下三角矩陣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ā)生瞬時(shí)點(diǎn)源污染事故，在t0=-120 min時(shí)刻該點(diǎn)源釋放污染物總量M=5 000 kg,降解系數(shù)K=0.001 min-1。

由公式(2)計(jì)算得到第i個(gè)監(jiān)測(cè)點(diǎn)在j時(shí)刻的污染物濃度的理論值Cobs(xi,yi,tj)(i=1,2,…,n;j=1,2,…,k)，為了使結(jié)果更貼近實(shí)際情況，考慮監(jiān)測(cè)誤差，對(duì)理論值附加一個(gè)高斯數(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)分布隨機(jī)數(shù)；E為監(jiān)測(cè)誤差水平。

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

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

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

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

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

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

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

(10)

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

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

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

圖3 不同E值下濃度值與E=0時(shí)濃度值差值分布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)測(cè)點(diǎn)數(shù)量對(duì)反演結(jié)果的影響

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

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

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

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

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

4 結(jié)論與展望

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