第五章 地下水污染預(yù)測(cè)

1、13/7/20221第五章第五章 地下水污染預(yù)測(cè)地下水污染預(yù)測(cè) 主要內(nèi)容主要內(nèi)容l 概述（了解）概述（了解）l 近似解析法近似解析法（了解）（了解）l 數(shù)值法數(shù)值法（熟悉）（熟悉）l 數(shù)理統(tǒng)計(jì)法數(shù)理統(tǒng)計(jì)法（熟悉）（熟悉）l 灰色預(yù)測(cè)法（了解）灰色預(yù)測(cè)法（了解）l 國(guó)外地下水模型軟件國(guó)外地下水模型軟件23/7/20222第一節(jié)第一節(jié) 概述概述 地下水污染預(yù)測(cè)預(yù)報(bào)包括：地下水污染預(yù)測(cè)預(yù)報(bào)包括：l 預(yù)報(bào)污染地下水分布預(yù)報(bào)污染地下水分布邊界的推進(jìn)情況邊界的推進(jìn)情況（即經(jīng)過(guò)（即經(jīng)過(guò)一定時(shí)間一定時(shí)間t后邊界推進(jìn)的距離后邊界推進(jìn)的距離x，或預(yù)測(cè)推進(jìn)一定，或預(yù)測(cè)推進(jìn)一定距離需要多長(zhǎng)時(shí)間）；距離需要多長(zhǎng)時(shí)間）；

2、l 確定污染地下水中污染物質(zhì)的確定污染地下水中污染物質(zhì)的濃度濃度在空間的分在空間的分布及隨時(shí)間的變化情況；布及隨時(shí)間的變化情況；l 預(yù)測(cè)污染地下水能否侵入附近的預(yù)測(cè)污染地下水能否侵入附近的水源地水源地，污染，污染物質(zhì)達(dá)到水源地所需的時(shí)間；物質(zhì)達(dá)到水源地所需的時(shí)間；l 預(yù)測(cè)防治地下水污染的預(yù)測(cè)防治地下水污染的措施措施的效果等。的效果等。33/7/20223第一節(jié)第一節(jié) 概述概述 預(yù)測(cè)評(píng)價(jià)的基本步驟如下：預(yù)測(cè)評(píng)價(jià)的基本步驟如下：l 建立預(yù)測(cè)模型；建立預(yù)測(cè)模型；l 模型識(shí)別（平衡性、穩(wěn)定性、敏感性）；模型識(shí)別（平衡性、穩(wěn)定性、敏感性）；l 模型參數(shù)估計(jì)；模型參數(shù)估計(jì)；l 模型方程的驗(yàn)證；模型方程的驗(yàn)

3、證；l 地下水污染預(yù)測(cè)。地下水污染預(yù)測(cè)。43/7/20224第二節(jié)第二節(jié) 近似解析解近似解析解前面介紹了一些解析解公式，但應(yīng)用起來(lái)仍十前面介紹了一些解析解公式，但應(yīng)用起來(lái)仍十分復(fù)雜，不便于實(shí)際應(yīng)用。分復(fù)雜，不便于實(shí)際應(yīng)用。實(shí)際中可以進(jìn)一步簡(jiǎn)化：實(shí)際中可以進(jìn)一步簡(jiǎn)化：l完全忽略彌散作用，將污染地下水的運(yùn)動(dòng)視為完全忽略彌散作用，將污染地下水的運(yùn)動(dòng)視為“活活塞式塞式”的推擠淡水的運(yùn)動(dòng)，假設(shè)兩者始終保持明顯的推擠淡水的運(yùn)動(dòng)，假設(shè)兩者始終保持明顯的鉛直分界面；的鉛直分界面；l將地下水動(dòng)力運(yùn)移和彌散分開(kāi)考慮，先以上述方法將地下水動(dòng)力運(yùn)移和彌散分開(kāi)考慮，先以上述方法計(jì)算，然后計(jì)算由于彌散所形成的過(guò)渡混合帶，

4、修計(jì)算，然后計(jì)算由于彌散所形成的過(guò)渡混合帶，修正按平均鋒面運(yùn)移的距離。正按平均鋒面運(yùn)移的距離。53/7/20225第二節(jié)第二節(jié) 近似解析解近似解析解 計(jì)算平均滲透鋒面的運(yùn)移距離計(jì)算平均滲透鋒面的運(yùn)移距離 ； 求出平均滲透鋒面處污染物的濃度求出平均滲透鋒面處污染物的濃度 ； 確定分界面彌散帶確定分界面彌散帶 和變形帶和變形帶 ，在平，在平均滲透鋒面距離上加上這兩個(gè)帶就可確定污染均滲透鋒面距離上加上這兩個(gè)帶就可確定污染水實(shí)際鋒面的運(yùn)移位置：水實(shí)際鋒面的運(yùn)移位置： nKmtLnDtLLLLLpdpd/cos/)(5 . 0LCdLpL相對(duì)密度差相對(duì)密度差巖層傾角巖層傾角63/7/20226第二節(jié)第二

5、節(jié) 近似解析解近似解析解l 承壓（或潛水）含水層承壓（或潛水）含水層 污染地下水在承壓含水層的運(yùn)動(dòng)，如果補(bǔ)給來(lái)源污染地下水在承壓含水層的運(yùn)動(dòng)，如果補(bǔ)給來(lái)源穩(wěn)定，可視為似穩(wěn)定運(yùn)動(dòng)，屬于一維平面平行流，則穩(wěn)定，可視為似穩(wěn)定運(yùn)動(dòng)，屬于一維平面平行流，則滲透速度為：滲透速度為：而實(shí)際流速為：而實(shí)際流速為：常數(shù)mqLHHKVex)(12nVuxx地下水的天然單寬流量地下水的天然單寬流量L73/7/20227第二節(jié)第二節(jié) 近似解析解近似解析解l承壓（或潛水）含水層承壓（或潛水）含水層eeeeCCCtxCCqmnLttmnqxxL00),(初始時(shí)刻污染鋒面的位置初始時(shí)刻污染鋒面的位置某時(shí)刻污染鋒面的位置某時(shí)

6、刻污染鋒面的位置天然地下水中污染物的濃度天然地下水中污染物的濃度初始分界面以?xún)?nèi)，地下水中污染物初始分界面以?xún)?nèi)，地下水中污染物的濃度的濃度若是潛水，則取平均厚度代替上述各式中的若是潛水，則取平均厚度代替上述各式中的m。L83/7/20228第三節(jié)第三節(jié) 數(shù)值解數(shù)值解由于含水層以及污染物質(zhì)遷移的復(fù)雜性，由于含水層以及污染物質(zhì)遷移的復(fù)雜性，解析解不易求得，通常采用數(shù)值解法。數(shù)解析解不易求得，通常采用數(shù)值解法。數(shù)值法在剖分后認(rèn)為每一單元內(nèi)均質(zhì)同性，值法在剖分后認(rèn)為每一單元內(nèi)均質(zhì)同性，將方程線性化。將方程線性化。 有限差分法有限差分法 有限單元法有限單元法 邊界元法邊界元法93/7/20229一、一、

7、有限差分法有限差分法有限差分法的有限差分法的基本思路基本思路是按照時(shí)間步長(zhǎng)和空間步是按照時(shí)間步長(zhǎng)和空間步長(zhǎng)將時(shí)間和空間區(qū)域剖分成若干網(wǎng)格，用未知函長(zhǎng)將時(shí)間和空間區(qū)域剖分成若干網(wǎng)格，用未知函數(shù)在網(wǎng)格節(jié)點(diǎn)上的值所構(gòu)成的差分商數(shù)在網(wǎng)格節(jié)點(diǎn)上的值所構(gòu)成的差分商近似代替近似代替所所用偏微分方程中出現(xiàn)的各階導(dǎo)數(shù)，從而把表示變用偏微分方程中出現(xiàn)的各階導(dǎo)數(shù)，從而把表示變量連續(xù)變化關(guān)系的偏微分方程離散為有限個(gè)代數(shù)量連續(xù)變化關(guān)系的偏微分方程離散為有限個(gè)代數(shù)方程，然后解此方程，然后解此線性代數(shù)方程線性代數(shù)方程，以求出污染物濃，以求出污染物濃度在各度在各網(wǎng)格節(jié)點(diǎn)網(wǎng)格節(jié)點(diǎn)上不同時(shí)刻的解。上不同時(shí)刻的解。 103/7/2

8、02210導(dǎo)數(shù)的差分近似導(dǎo)數(shù)的差分近似設(shè)設(shè)f(x)為任一足夠光滑的函數(shù)，把為任一足夠光滑的函數(shù)，把f(x)沿沿x的正向的正向展開(kāi)為泰勒級(jí)數(shù)：展開(kāi)為泰勒級(jí)數(shù)：上式稱(chēng)為上式稱(chēng)為f(x)的的一階向前差分一階向前差分 333222! 3! 2)()(dxfdxdxfdxdxdfxxfxxf)()()(xRxxfxxfdxdfxxfxxfdxdf)()(113/7/202211導(dǎo)數(shù)的差分近似導(dǎo)數(shù)的差分近似如果把如果把f(x)對(duì)對(duì)x的負(fù)向展開(kāi)為泰勒級(jí)數(shù)：的負(fù)向展開(kāi)為泰勒級(jí)數(shù)：上式稱(chēng)為上式稱(chēng)為f(x)的的一階向后差分一階向后差分 333222! 3! 2)()(dxfdxdxfdxdxdfxxfxxf)()

9、()(xOxxxfxfdxdfxxxfxfdxdf)()(123/7/202212導(dǎo)數(shù)的差分近似導(dǎo)數(shù)的差分近似把上面兩個(gè)差分方程相減，可得：把上面兩個(gè)差分方程相減，可得：上式稱(chēng)為上式稱(chēng)為f(x)的的一階中心差分一階中心差分同理可得同理可得f(x)的二階導(dǎo)數(shù)的差分公式：的二階導(dǎo)數(shù)的差分公式：xxxfxxfdxdf2)()(222)()(2)(xxxfxfxxfdxfd133/7/202213導(dǎo)數(shù)的差分近似導(dǎo)數(shù)的差分近似若把上述一元函數(shù)若把上述一元函數(shù)f(x)換為濃度函數(shù)，并按網(wǎng)格換為濃度函數(shù)，并按網(wǎng)格節(jié)點(diǎn)編號(hào)，則對(duì)節(jié)點(diǎn)編號(hào)，則對(duì)x的一階向前、向后和中心差分的一階向前、向后和中心差分分別為：分別為

10、：xCCxcxCCxcxCCxcnkjinkjinkjinkjinkjinkjinkjinkjinkji2|,1,1),(,1,),(,1),(同樣可以寫(xiě)出同樣可以寫(xiě)出對(duì)時(shí)間對(duì)時(shí)間t的向的向前、向后和中前、向后和中心差分公式心差分公式143/7/202214導(dǎo)數(shù)的差分近似導(dǎo)數(shù)的差分近似C(x,y,z,t)對(duì)對(duì)x的二階差分公式為：的二階差分公式為：2,1,1),(222|xCCCxcnkjinkjinkjinkji因此，在任何一個(gè)節(jié)點(diǎn)處，濃度因此，在任何一個(gè)節(jié)點(diǎn)處，濃度C(x,y,z,t)的一階和二階偏導(dǎo)數(shù)均可用該節(jié)點(diǎn)及其的一階和二階偏導(dǎo)數(shù)均可用該節(jié)點(diǎn)及其相鄰節(jié)點(diǎn)上濃度值的線性組合表示。相鄰節(jié)點(diǎn)

11、上濃度值的線性組合表示。153/7/202215一維彌散方程的差分格式一維彌散方程的差分格式設(shè)在無(wú)限含水層中，存在一維均勻流場(chǎng)，其滲透設(shè)在無(wú)限含水層中，存在一維均勻流場(chǎng)，其滲透速度速度V=nu，流動(dòng)方向?yàn)椋鲃?dòng)方向?yàn)閤軸正向。則一維彌散方軸正向。則一維彌散方程為：程為： Tt,0 Lx0 22xcuxcDtcL163/7/202216一維彌散方程的差分格式一維彌散方程的差分格式對(duì)時(shí)間區(qū)域?qū)r(shí)間區(qū)域0,T和空間區(qū)域和空間區(qū)域0,L都作等距剖分，都作等距剖分，設(shè)時(shí)間步長(zhǎng)為設(shè)時(shí)間步長(zhǎng)為t，空間步長(zhǎng)為，空間步長(zhǎng)為x，把第，把第i個(gè)節(jié)點(diǎn)個(gè)節(jié)點(diǎn)xi處在處在tn時(shí)刻的濃度記為時(shí)刻的濃度記為Cin。173/7

12、/202217一維彌散方程的差分格式一維彌散方程的差分格式顯示差分格式顯示差分格式 將將 及及 差分公式中的濃度取為差分公式中的濃度取為tn時(shí)刻的值，時(shí)刻的值，便可得到其顯示差分格式：便可得到其顯示差分格式：xC22xCxCCuxCCCDtCCnininininiLnini22112111niLniLniLniCxtuxtDCxtDCxtuxtDC122121)2()21()2(183/7/202218一維彌散方程的差分格式一維彌散方程的差分格式顯示差分格式的穩(wěn)定條件顯示差分格式的穩(wěn)定條件 在應(yīng)用顯示差分格式時(shí)，剖分的步長(zhǎng)必須滿(mǎn)足：在應(yīng)用顯示差分格式時(shí)，剖分的步長(zhǎng)必須滿(mǎn)足：1212xtuxtD

13、L及193/7/202219一維彌散方程的差分格式一維彌散方程的差分格式隱示差分格式隱示差分格式 將將 及及 差分公式中的濃度取為差分公式中的濃度取為tn1時(shí)刻的時(shí)刻的值，便可得到其隱示差分格式：值，便可得到其隱示差分格式：xC22xCxCCuxCCCDtCCnininininiLnini2211112111111niniLniLniLCCxtuxtDCxtDCxtuxtD11212112)2()21()2(203/7/202220一維彌散方程的差分格式一維彌散方程的差分格式隱示差分格式的穩(wěn)定條件隱示差分格式的穩(wěn)定條件 在應(yīng)用隱示差分格式時(shí)，剖分的步長(zhǎng)必須滿(mǎn)足：在應(yīng)用隱示差分格式時(shí)，剖分的步長(zhǎng)

14、必須滿(mǎn)足：xDuL2|213/7/202221一維彌散方程的差分格式一維彌散方程的差分格式隱示差分格式的矩陣形式隱示差分格式的矩陣形式1nCA223/7/202222一維彌散方程的差分格式一維彌散方程的差分格式Crank-Nicolson差分格式差分格式 取顯示差分格式和隱示差分格式的平均，便可得取顯示差分格式和隱示差分格式的平均，便可得到到C-N差分格式：差分格式：)22(21)22(211111211111112111xCCuxCCCDxCCuxCCCDtCCnininininiLnininininiLnini233/7/202223一維彌散方程的差分格式一維彌散方程的差分格式Crank-

15、Nicolson差分格式差分格式 整理后可得：整理后可得：niLniLniLniLniLniLCxtuxtDCxtDCxtuxtDCxtuxtDCxtDCxtuxtD1221211212112)2()1 ( 2)2()2()1 ( 2)2(1nAC243/7/202224追趕法：追趕法：iniiniiniiCCC11111方程的統(tǒng)一形式：方程的統(tǒng)一形式：253/7/202225追趕法：追趕法：102111121111nnnCxDCC引入邊界條件：引入邊界條件：當(dāng)當(dāng)i=1時(shí)：時(shí)：當(dāng)當(dāng)i=M時(shí)：時(shí)：MMMMMMMnMMnMMnMnMnMCCCCC2211111111其中則令263/7/202226

16、追趕法：追趕法：令令A(yù)C=MMnMMMMMMMCCCCCA121112111133322211273/7/202227追趕法：追趕法：追趕法解方程組的方法：追趕法解方程組的方法：先把對(duì)角矩陣先把對(duì)角矩陣A分解為：分解為：A=BW其中：其中：MMMMbabababB112211111121MwwwW283/7/202228追趕法：追趕法：由于由于A=BW1112111332333221222111MMMMMMMMMMwabawbwabawbwabawbwabawbbBWA293/7/202229追趕法：追趕法：由此可得：由此可得： b1=1 b1w1=1 w1=1/b1 b2+a2w1=2 bi

17、wi=i wi=i/bi (i=1M-1) bi+aiwi-1=i bi=i-aiwi-1 bM=M-aMwM-1303/7/202230追趕法：追趕法：由于由于A=BW B(WC)=又令：又令：WC=Y 其中：其中：Y=(y1, y2, , yM)T BY=MMMMMMMMyyyybabababB12112111221313/7/202231追趕法：追趕法：由以上等式可知：由以上等式可知： b1y1=1 y1=1 /b1 a2y1+b2y2=2 y2=(2-a2y1)/ b2 aiyi-1+biyi=i yi=(i-aiyi-1)/ bi aMyM-1+bMyM=M yM=(M-aMyM-1

18、)/ bM至此可以求出所有的至此可以求出所有的Y(y1, y2, , yM)T值。值。323/7/202232追趕法：追趕法：又因?yàn)橛忠驗(yàn)閃C=Y，即：，即：MMnMMMyyyyCCCCwwwCW12111211211111333/7/202233追趕法：追趕法：根據(jù)上述等式，可得：根據(jù)上述等式，可得：MnMMnMMnMniiininnnnyCyCwCCwyCyCwCyCwC111111111213212112111343/7/202234追趕法：追趕法：追：追：根據(jù)根據(jù)BY=，求出，求出bi和和wi，并求出，并求出yi；趕：趕：根據(jù)根據(jù)WC=Y，求出，求出Cin+1。353/7/202235

19、二維彌散方程的差分格式二維彌散方程的差分格式若多孔介質(zhì)是各向同性的，流場(chǎng)為均勻的若多孔介質(zhì)是各向同性的，流場(chǎng)為均勻的一維流場(chǎng)，滲透速度一維流場(chǎng)，滲透速度V=nu，流動(dòng)方向取，流動(dòng)方向取x軸軸的正向，與流速垂直為的正向，與流速垂直為y方向，對(duì)流擴(kuò)散方方向，對(duì)流擴(kuò)散方程可簡(jiǎn)化為：程可簡(jiǎn)化為： xcuycDxcDtcTL2222363/7/202236二維彌散方程的差分格式二維彌散方程的差分格式將研究區(qū)域網(wǎng)格化將研究區(qū)域網(wǎng)格化tntyjcyxiaxNttMcdyMabxnji,/,/ )(,/ )(373/7/202237二維彌散方程的差分格式二維彌散方程的差分格式建立矩形網(wǎng)差分格式建立矩形網(wǎng)差分格

20、式l隱式差分格式隱式差分格式tCCtcxCCxcnjinjinjinji,1,1,11,12211,1,11,2221,11,1,12222yCCCycxCCCxcnjinjinjinjinjinji383/7/202238二維彌散方程的差分格式二維彌散方程的差分格式于是，內(nèi)節(jié)點(diǎn)的差分方程為于是，內(nèi)節(jié)點(diǎn)的差分方程為xCCVyCCCDxCCCDtCCnjinjinjinjinjiTnjinjinjiLnjinji2221,11,1211,1,11,21,11,1,1,1,393/7/202239二維彌散方程的差分格式二維彌散方程的差分格式整理后可得整理后可得njinjiTnjiTnjiLnjiT

21、LnjiLCCyTDCyTDCxtVxtDCytDxtDCxtVxtD,11,211,21,121,221,12)2()221()2(403/7/202240二維彌散方程的差分格式二維彌散方程的差分格式邊界條件邊界條件初始條件初始條件),(),(),(),(1,10,1,1,0ninMinininjnjMnjnjtdxCCtcxCCtybCCtyaCC),(00,jijiyxCC413/7/202241二維彌散方程的差分格式二維彌散方程的差分格式將網(wǎng)格取為正方形，即將網(wǎng)格取為正方形，即x=y=h，方程就，方程就可簡(jiǎn)化為：可簡(jiǎn)化為：njinjiTnjiTnjiLnjiTLnjiLCTCDhTCD

22、hCtVhtDhCtDtDhhCtVhtDh,11,211,21,121,221,1211)(21)22(1)(21423/7/202242二維彌散方程的差分格式二維彌散方程的差分格式交替方向差分格式（交替方向差分格式（ADI） ADI法的實(shí)質(zhì)是法的實(shí)質(zhì)是將二維隱式差分格式化為一將二維隱式差分格式化為一維的隱式差分格式維的隱式差分格式。 在第在第n層與第層與第n+1層兩層正中間引入一個(gè)過(guò)層兩層正中間引入一個(gè)過(guò)渡層渡層n+1/2。433/7/202243二維彌散方程的差分格式二維彌散方程的差分格式交替方向差分格式（交替方向差分格式（ADI） 分兩步走：分兩步走：在前半時(shí)間步長(zhǎng)內(nèi)，對(duì)在前半時(shí)間步長(zhǎng)

23、內(nèi)，對(duì)x偏導(dǎo)用隱式偏導(dǎo)用隱式表示，對(duì)表示，對(duì)y的偏導(dǎo)用顯示表示；在后半時(shí)間步長(zhǎng)的偏導(dǎo)用顯示表示；在后半時(shí)間步長(zhǎng)內(nèi)，對(duì)內(nèi)，對(duì)x偏導(dǎo)用顯示表示，對(duì)偏導(dǎo)用顯示表示，對(duì)y的偏導(dǎo)用隱式表的偏導(dǎo)用隱式表示。示。443/7/202244二維彌散方程的差分格式二維彌散方程的差分格式交替方向差分格式（交替方向差分格式（ADI） 前半時(shí)間步長(zhǎng)：前半時(shí)間步長(zhǎng)： 21,1,1,12,1,1,2222/212121212121yCCCDxCCVxCCCDtCCnjinjinjiTnjinjinjinjinjiLnjinji453/7/202245二維彌散方程的差分格式二維彌散方程的差分格式交替方向差分格式（交替方向差分

24、格式（ADI） 后半時(shí)間步長(zhǎng)：后半時(shí)間步長(zhǎng)： 211,1,11,1,12,1,1,1,2222/212121212121yCCCDxCCVxCCCDtCCnjinjinjiTnjinjinjinjinjiLnjinji463/7/202246二維彌散方程的差分格式二維彌散方程的差分格式計(jì)算步驟：計(jì)算步驟：（1）根據(jù)初始條件或前一時(shí)段計(jì)算結(jié)果，可知）根據(jù)初始條件或前一時(shí)段計(jì)算結(jié)果，可知tn時(shí)刻的濃度值，那么第一個(gè)方程中只有三個(gè)未時(shí)刻的濃度值，那么第一個(gè)方程中只有三個(gè)未知數(shù)（知數(shù)（n+1/2時(shí)刻）。如果固定時(shí)刻）。如果固定j值，對(duì)值，對(duì)i的聯(lián)立的聯(lián)立方程組就是三對(duì)角線的（追趕法），即方程組就是三對(duì)

25、角線的（追趕法），即沿沿x方向方向追趕追趕，可求出過(guò)渡層（，可求出過(guò)渡層（n+1/2）層上的全部濃度）層上的全部濃度值。值。 473/7/202247二維彌散方程的差分格式二維彌散方程的差分格式計(jì)算步驟：計(jì)算步驟：（2）根據(jù)前一步求出的過(guò)渡層（）根據(jù)前一步求出的過(guò)渡層（n+1/2）層上的）層上的濃度值，再按同樣的方法求出第二個(gè)方程中的濃度值，再按同樣的方法求出第二個(gè)方程中的（n+1）時(shí)刻的濃度值。即對(duì)于一個(gè)固定的）時(shí)刻的濃度值。即對(duì)于一個(gè)固定的i值，值，方程也是一個(gè)三對(duì)角方程組（追趕法），方程也是一個(gè)三對(duì)角方程組（追趕法），沿沿y方方向追趕向追趕。 483/7/202248有限單元法基本原理有

26、限單元法基本原理在在有限差分法有限差分法中，研究區(qū)域被一系列的網(wǎng)格節(jié)點(diǎn)中，研究區(qū)域被一系列的網(wǎng)格節(jié)點(diǎn)所代替，在各個(gè)節(jié)點(diǎn)上以差分代替導(dǎo)數(shù)而得到求所代替，在各個(gè)節(jié)點(diǎn)上以差分代替導(dǎo)數(shù)而得到求解相應(yīng)偏微分方程的線性代數(shù)方程組。解相應(yīng)偏微分方程的線性代數(shù)方程組。有限單元法有限單元法是把研究區(qū)域剖分為有限個(gè)子區(qū)域，是把研究區(qū)域剖分為有限個(gè)子區(qū)域，在每個(gè)子區(qū)域上用某種插值函數(shù)來(lái)近似待求解的在每個(gè)子區(qū)域上用某種插值函數(shù)來(lái)近似待求解的未知函數(shù)，從而得到求解相應(yīng)偏微分方程的線性未知函數(shù)，從而得到求解相應(yīng)偏微分方程的線性代數(shù)方程組。代數(shù)方程組。 493/7/202249邊界元法基本原理邊界元法基本原理通過(guò)把求解區(qū)域

27、的邊界剖分為若干單元，將求函數(shù)通過(guò)把求解區(qū)域的邊界剖分為若干單元，將求函數(shù)的解簡(jiǎn)化為求節(jié)點(diǎn)的函數(shù)值，求解積分方程就化為的解簡(jiǎn)化為求節(jié)點(diǎn)的函數(shù)值，求解積分方程就化為求解一組線性代數(shù)方程，在求出邊界上的物理量以求解一組線性代數(shù)方程，在求出邊界上的物理量以后，計(jì)算區(qū)域內(nèi)部的任一點(diǎn)的未知量可通過(guò)邊界上后，計(jì)算區(qū)域內(nèi)部的任一點(diǎn)的未知量可通過(guò)邊界上的已知量求出。的已知量求出。與有限元法相比，具有降低維數(shù)、輸入數(shù)據(jù)準(zhǔn)備簡(jiǎn)與有限元法相比，具有降低維數(shù)、輸入數(shù)據(jù)準(zhǔn)備簡(jiǎn)單、計(jì)算工作量小、精度比有限元法高等優(yōu)點(diǎn)。單、計(jì)算工作量小、精度比有限元法高等優(yōu)點(diǎn)。在處理非均質(zhì)、非線性、非穩(wěn)定等復(fù)雜問(wèn)題時(shí)計(jì)算在處理非均質(zhì)、非線

28、性、非穩(wěn)定等復(fù)雜問(wèn)題時(shí)計(jì)算效率低，不能將對(duì)流項(xiàng)化成沿邊界的積分。效率低，不能將對(duì)流項(xiàng)化成沿邊界的積分。 503/7/202250第四節(jié)第四節(jié) 數(shù)理統(tǒng)計(jì)法數(shù)理統(tǒng)計(jì)法地下水水質(zhì)受多種因素綜合影響，具有不地下水水質(zhì)受多種因素綜合影響，具有不確定性，可用數(shù)理統(tǒng)計(jì)法預(yù)測(cè)水質(zhì)變化；確定性，可用數(shù)理統(tǒng)計(jì)法預(yù)測(cè)水質(zhì)變化；主要包括：時(shí)間序列分析法、相關(guān)分析法、主要包括：時(shí)間序列分析法、相關(guān)分析法、概率統(tǒng)計(jì)法。概率統(tǒng)計(jì)法。513/7/202251（1）時(shí)間序列分析法）時(shí)間序列分析法基本思路基本思路：分析時(shí)間序列的變化特征，選：分析時(shí)間序列的變化特征，選擇適當(dāng)?shù)哪Ｊ浇㈩A(yù)測(cè)模型，利用模型進(jìn)擇適當(dāng)?shù)哪Ｊ浇㈩A(yù)測(cè)模型，

29、利用模型進(jìn)行趨勢(shì)外推預(yù)測(cè)，對(duì)模型預(yù)測(cè)值進(jìn)行評(píng)價(jià)行趨勢(shì)外推預(yù)測(cè)，對(duì)模型預(yù)測(cè)值進(jìn)行評(píng)價(jià)和修正，得到預(yù)測(cè)結(jié)果。和修正，得到預(yù)測(cè)結(jié)果。523/7/202252（1）時(shí)間序列分析法）時(shí)間序列分析法-平均數(shù)平均數(shù)平均數(shù)預(yù)測(cè)平均數(shù)預(yù)測(cè)：對(duì)動(dòng)態(tài)數(shù)列中若干個(gè)時(shí)期的數(shù)值，：對(duì)動(dòng)態(tài)數(shù)列中若干個(gè)時(shí)期的數(shù)值，求其算術(shù)平均值或加權(quán)平均值，作為下一個(gè)時(shí)求其算術(shù)平均值或加權(quán)平均值，作為下一個(gè)時(shí)期的預(yù)測(cè)值。期的預(yù)測(cè)值。niinCnC111niiinCWnC11)(1533/7/202253（1）時(shí)間序列分析法）時(shí)間序列分析法-指數(shù)平滑指數(shù)平滑一次指數(shù)平滑預(yù)測(cè)一次指數(shù)平滑預(yù)測(cè)：對(duì)動(dòng)態(tài)數(shù)列中若干個(gè)時(shí)期：對(duì)動(dòng)態(tài)數(shù)列中若干個(gè)時(shí)期的數(shù)值

30、，求其算術(shù)平均值或加權(quán)平均值，作為的數(shù)值，求其算術(shù)平均值或加權(quán)平均值，作為下一個(gè)時(shí)期的預(yù)測(cè)值。下一個(gè)時(shí)期的預(yù)測(cè)值。 221)1()1 ()1 (ttttCaaCaaaCS)1(121)1()1 ()1 ()1 ( ttttttSaaCCaaaCaaCS543/7/202254（1）時(shí)間序列分析法）時(shí)間序列分析法-指數(shù)平滑指數(shù)平滑一次指數(shù)平滑法是以最近周期的一次指數(shù)平滑一次指數(shù)平滑法是以最近周期的一次指數(shù)平滑值作為下一周期的預(yù)測(cè)值，即：值作為下一周期的預(yù)測(cè)值，即：其中平滑系數(shù)其中平滑系數(shù)a的選擇直接影響預(yù)測(cè)效果。的選擇直接影響預(yù)測(cè)效果。如如果時(shí)間序列長(zhǎng)期趨勢(shì)穩(wěn)定，應(yīng)取較小的果時(shí)間序列長(zhǎng)期趨勢(shì)穩(wěn)定

31、，應(yīng)取較小的a值；值；若時(shí)間序列不穩(wěn)定，應(yīng)取較大的若時(shí)間序列不穩(wěn)定，應(yīng)取較大的a值。值。)1(1)1(1)1 (ttttSaaCSC553/7/202255（1）時(shí)間序列分析法）時(shí)間序列分析法-指數(shù)平滑指數(shù)平滑a越小，對(duì)數(shù)據(jù)的平滑能力越強(qiáng)，但對(duì)數(shù)據(jù)變化的敏感越小，對(duì)數(shù)據(jù)的平滑能力越強(qiáng)，但對(duì)數(shù)據(jù)變化的敏感性越差；性越差；a越大，對(duì)數(shù)據(jù)的平滑能力越差，但對(duì)數(shù)據(jù)變?cè)酱?，?duì)數(shù)據(jù)的平滑能力越差，但對(duì)數(shù)據(jù)變化的敏感性越強(qiáng)?；拿舾行栽綇?qiáng)。563/7/202256（1）時(shí)間序列分析法）時(shí)間序列分析法-指數(shù)平滑指數(shù)平滑初始值初始值S0(1)對(duì)預(yù)測(cè)結(jié)果的影響很大，一般可取第一個(gè)數(shù)對(duì)預(yù)測(cè)結(jié)果的影響很大，一般可取第

32、一個(gè)數(shù)據(jù)（據(jù)（x1）作為初始值；若總體數(shù)據(jù)量較少，則可取最）作為初始值；若總體數(shù)據(jù)量較少，則可取最初幾個(gè)數(shù)據(jù)的平均值作為初始值。初幾個(gè)數(shù)據(jù)的平均值作為初始值。573/7/202257（1）時(shí)間序列分析法）時(shí)間序列分析法-指數(shù)平滑指數(shù)平滑二次指數(shù)平滑預(yù)測(cè)二次指數(shù)平滑預(yù)測(cè)：對(duì)一次指數(shù)平滑值再做一：對(duì)一次指數(shù)平滑值再做一次指數(shù)平滑：次指數(shù)平滑：)2(1)1()2()1 (tttSaaSS583/7/202258（1）時(shí)間序列分析法）時(shí)間序列分析法-指數(shù)平滑指數(shù)平滑二次指數(shù)平滑值不直接用于預(yù)測(cè)，而是根據(jù)滯二次指數(shù)平滑值不直接用于預(yù)測(cè)，而是根據(jù)滯后偏差建立線性預(yù)測(cè)模型：后偏差建立線性預(yù)測(cè)模型：其中：其中

33、：TbaCttTt)2()1(2tttSSa1)2()1(tttSSaab593/7/202259（1）時(shí)間序列分析法）時(shí)間序列分析法-指數(shù)平滑指數(shù)平滑若一、二次指數(shù)平滑值不滿(mǎn)足預(yù)測(cè)要求，應(yīng)采若一、二次指數(shù)平滑值不滿(mǎn)足預(yù)測(cè)要求，應(yīng)采用三次指數(shù)平滑法建立非線性預(yù)測(cè)模型：用三次指數(shù)平滑法建立非線性預(yù)測(cè)模型：(3)1)2()3()1 (tttSaaSS603/7/202260（1）時(shí)間序列分析法）時(shí)間序列分析法-指數(shù)平滑指數(shù)平滑三次指數(shù)平滑法建立的非線性預(yù)測(cè)模型：三次指數(shù)平滑法建立的非線性預(yù)測(cè)模型：)3()2()1(333ttttSSSa)34()810(2)56()12()3()2()1(2ttt

34、tSaSaSaaab2TtttTtcTbaC2)12()3()2()1(2ttttSSSaac613/7/202261（2）相關(guān)分析法相關(guān)分析法線性相關(guān)預(yù)測(cè)線性相關(guān)預(yù)測(cè)：建立直角坐標(biāo)下數(shù)據(jù)散點(diǎn)圖，：建立直角坐標(biāo)下數(shù)據(jù)散點(diǎn)圖，用一元線性回歸方程近似描述數(shù)據(jù)點(diǎn)的分布規(guī)用一元線性回歸方程近似描述數(shù)據(jù)點(diǎn)的分布規(guī)律：律：bxay623/7/202262（2）相關(guān)分析法相關(guān)分析法根據(jù)最小二乘法根據(jù)最小二乘法(要求殘差要求殘差RSS最小最小)可求得：可求得：用相關(guān)系數(shù)用相關(guān)系數(shù)r來(lái)判斷回歸方程實(shí)用性：來(lái)判斷回歸方程實(shí)用性：22222)()(iiiiiiiiiiiiixnxyxnyxbxnxxyxyxa22)

35、()()(yyxxyyxxriiii2)(iibxayRSS633/7/202263（2）相關(guān)分析法相關(guān)分析法曲性相關(guān)預(yù)測(cè)曲性相關(guān)預(yù)測(cè)：將直角坐標(biāo)下散點(diǎn)圖擬合形成：將直角坐標(biāo)下散點(diǎn)圖擬合形成的有規(guī)律的曲線進(jìn)行函數(shù)轉(zhuǎn)換，變?yōu)榍笆龅囊坏挠幸?guī)律的曲線進(jìn)行函數(shù)轉(zhuǎn)換，變?yōu)榍笆龅囊辉€性回歸方程進(jìn)行計(jì)算。元線性回歸方程進(jìn)行計(jì)算。如：如：baxy bXaYyYxXlglg,lg則令643/7/202264國(guó)外地下水模擬軟件的發(fā)展現(xiàn)狀與趨勢(shì)國(guó)外地下水模擬軟件的發(fā)展現(xiàn)狀與趨勢(shì) 利用數(shù)值模型對(duì)地下水流和溶質(zhì)運(yùn)移問(wèn)題進(jìn)行利用數(shù)值模型對(duì)地下水流和溶質(zhì)運(yùn)移問(wèn)題進(jìn)行模擬的方法以其有效性、靈活性和相對(duì)廉價(jià)性模擬的方法以其有

36、效性、靈活性和相對(duì)廉價(jià)性逐漸成為地下水研究領(lǐng)域的一種不可或缺的重逐漸成為地下水研究領(lǐng)域的一種不可或缺的重要方法，并受到越來(lái)越大的重視和廣泛的應(yīng)用。要方法，并受到越來(lái)越大的重視和廣泛的應(yīng)用。一個(gè)完整的地下水模擬過(guò)程包含一個(gè)完整的地下水模擬過(guò)程包含3個(gè)部分：個(gè)部分：前處前處理、模型計(jì)算和后處理理、模型計(jì)算和后處理。653/7/202265國(guó)外地下水模擬軟件的發(fā)展現(xiàn)狀與趨勢(shì)國(guó)外地下水模擬軟件的發(fā)展現(xiàn)狀與趨勢(shì)前處理前處理是指在進(jìn)行模擬計(jì)算之前對(duì)計(jì)算過(guò)程中是指在進(jìn)行模擬計(jì)算之前對(duì)計(jì)算過(guò)程中所需數(shù)據(jù)的整理、組織、輸入及計(jì)算網(wǎng)格的編所需數(shù)據(jù)的整理、組織、輸入及計(jì)算網(wǎng)格的編號(hào)與生成。號(hào)與生成。模型計(jì)算模型計(jì)算

37、是進(jìn)行地下水流動(dòng)或水質(zhì)運(yùn)移正反演是進(jìn)行地下水流動(dòng)或水質(zhì)運(yùn)移正反演計(jì)算，常用的方法主要有：有限差分法、有限計(jì)算，常用的方法主要有：有限差分法、有限元法、邊界元法等。元法、邊界元法等。后處理后處理是將計(jì)算所產(chǎn)生的結(jié)果數(shù)據(jù)，用圖形或是將計(jì)算所產(chǎn)生的結(jié)果數(shù)據(jù)，用圖形或表格顯示或存放起來(lái)，以供研究人員方便地進(jìn)表格顯示或存放起來(lái)，以供研究人員方便地進(jìn)行分析和使用。行分析和使用。663/7/202266國(guó)外地下水模擬軟件的發(fā)展現(xiàn)狀與趨勢(shì)國(guó)外地下水模擬軟件的發(fā)展現(xiàn)狀與趨勢(shì)通過(guò)近二十年的研究與發(fā)展，國(guó)際上已經(jīng)形成通過(guò)近二十年的研究與發(fā)展，國(guó)際上已經(jīng)形成了一批非常有影響的地下水模擬的軟件，它們了一批非常有影響的地

38、下水模擬的軟件，它們今天在國(guó)際地下水模擬研究領(lǐng)域依舊非?；钴S，今天在國(guó)際地下水模擬研究領(lǐng)域依舊非?；钴S，如如MODFLOW、FEFLOW、MT3DMS、MT3D99、PEST、MODPATH、UCODE等。等。 673/7/202267MODFLOW MODFLOW是由美國(guó)地質(zhì)調(diào)查局的是由美國(guó)地質(zhì)調(diào)查局的McDonald和和Harbaugh于于80年代開(kāi)發(fā)出來(lái)的一套專(zhuān)門(mén)用年代開(kāi)發(fā)出來(lái)的一套專(zhuān)門(mén)用于孔隙介質(zhì)中于孔隙介質(zhì)中三維有限差分三維有限差分地下水流數(shù)值模擬地下水流數(shù)值模擬的軟件。的軟件。l程序結(jié)構(gòu)的程序結(jié)構(gòu)的模塊化模塊化（模擬抽水引起地面沉降的子（模擬抽水引起地面沉降的子程序包、模擬水平流動(dòng)

39、障礙的子程序包等）。程序包、模擬水平流動(dòng)障礙的子程序包等）。l離散方法的離散方法的簡(jiǎn)單化簡(jiǎn)單化（有限差分法（有限差分法36036018個(gè)個(gè)網(wǎng)格單元）。網(wǎng)格單元）。lMODFLOW引進(jìn)了引進(jìn)了應(yīng)力期應(yīng)力期(Stress Period)概念。概念。 l求解方法的多樣化（強(qiáng)隱式法、逐次超松弛迭代求解方法的多樣化（強(qiáng)隱式法、逐次超松弛迭代法、預(yù)調(diào)共軛梯度法）。法、預(yù)調(diào)共軛梯度法）。 683/7/202268MT3D99 MT3D99是美國(guó)是美國(guó)S.S.Papodopulos & Associates公司的公司的Zheng博士設(shè)計(jì)的模擬博士設(shè)計(jì)的模擬三維地下水溶質(zhì)三維地下水溶質(zhì)運(yùn)移運(yùn)移程序程序M

40、T3D(1990)的增強(qiáng)版，成為目前世的增強(qiáng)版，成為目前世界上首屈一指的溶質(zhì)運(yùn)移模擬軟件。界上首屈一指的溶質(zhì)運(yùn)移模擬軟件。lMT3D99能夠模擬地下水系統(tǒng)中的平流、擴(kuò)散、衰能夠模擬地下水系統(tǒng)中的平流、擴(kuò)散、衰減、溶質(zhì)化學(xué)反應(yīng)、線性與非線性吸附作用等現(xiàn)象，減、溶質(zhì)化學(xué)反應(yīng)、線性與非線性吸附作用等現(xiàn)象，能夠?qū)Τ袎汉畬?，不承壓含水層，承壓與不承壓能夠?qū)Τ袎汉畬?，不承壓含水層，承壓與不承壓交替的含水層以及傾斜的和單元厚度變化的含水層交替的含水層以及傾斜的和單元厚度變化的含水層進(jìn)行空間離散。進(jìn)行空間離散。 lMT3D99提供了豐富的求解方法（標(biāo)準(zhǔn)有限差分法、提供了豐富的求解方法（標(biāo)準(zhǔn)有限差分法、基于

41、基于Eulerian-Lagrangian的粒子跟蹤方法和高階有的粒子跟蹤方法和高階有限體積限體積TVD方法）。方法）。 693/7/202269PEST2000 PEST2000是由澳大利亞是由澳大利亞Watermark Computing公公司開(kāi)發(fā)的司開(kāi)發(fā)的PEST新版本，是功能強(qiáng)大的模型獨(dú)立的新版本，是功能強(qiáng)大的模型獨(dú)立的參數(shù)估計(jì)參數(shù)估計(jì)程序程序。lPEST2000利用一個(gè)強(qiáng)有力的數(shù)值反演算法來(lái)利用一個(gè)強(qiáng)有力的數(shù)值反演算法來(lái)“控制控制”運(yùn)行中的模型，程序在每次模擬之后自運(yùn)行中的模型，程序在每次模擬之后自動(dòng)調(diào)整所選擇的模型參數(shù)，直到將校正的目標(biāo)動(dòng)調(diào)整所選擇的模型參數(shù)，直到將校正的目標(biāo)最小化為止。最小化為止。lPEST2000引進(jìn)了參數(shù)估計(jì)的最新技術(shù)引進(jìn)了參數(shù)估計(jì)的最新技術(shù)預(yù)測(cè)預(yù)測(cè)分析。分析。 703/7/202270FEFLOW FEFLOW是加拿大是加拿大Wate