第五章 地下水污染预测

1、11 共81页1 第五章第五章 地下水污染预测地下水污染预测 主要内容主要内容 l 概述（了解）概述（了解） l 近似解析法近似解析法（了解）（了解） l 数值法数值法（熟悉）（熟悉） l 数理统计法数理统计法（熟悉）（熟悉） l 灰色预测法（了解）灰色预测法（了解） l 国外地下水模型软件国外地下水模型软件 22 共81页2 第一节第一节 概述概述 地下水污染预测预报包括：地下水污染预测预报包括： l 预报污染地下水分布预报污染地下水分布边界的推进情况边界的推进情况（即经过（即经过 一定时间一定时间t后边界推进的距离后边界推进的距离x，或预测推进一定，或预测推进一定 距离需要多长时间）；距离

2、需要多长时间）； l 确定污染地下水中污染物质的确定污染地下水中污染物质的浓度浓度在空间的分在空间的分 布及随时间的变化情况；布及随时间的变化情况； l 预测污染地下水能否侵入附近的预测污染地下水能否侵入附近的水源地水源地，污染，污染 物质达到水源地所需的时间；物质达到水源地所需的时间； l 预测防治地下水污染的预测防治地下水污染的措施措施的效果等。的效果等。 33 共81页3 第一节第一节 概述概述 预测评价的基本步骤如下：预测评价的基本步骤如下： l 建立预测模型；建立预测模型； l 模型识别（平衡性、稳定性、敏感性）；模型识别（平衡性、稳定性、敏感性）； l 模型参数估计；模型参数估计；

3、 l 模型方程的验证；模型方程的验证； l 地下水污染预测。地下水污染预测。 44 共81页4 第二节第二节 近似解析解近似解析解 前面介绍了一些解析解公式，但应用起来仍十前面介绍了一些解析解公式，但应用起来仍十 分复杂，不便于实际应用。分复杂，不便于实际应用。 实际中可以进一步简化：实际中可以进一步简化： l完全忽略弥散作用，将污染地下水的运动视为完全忽略弥散作用，将污染地下水的运动视为“活活 塞式塞式”的推挤淡水的运动，假设两者始终保持明显的推挤淡水的运动，假设两者始终保持明显 的铅直分界面；的铅直分界面； l将地下水动力运移和弥散分开考虑，先以上述方法将地下水动力运移和弥散分开考虑，先以

4、上述方法 计算，然后计算由于弥散所形成的过渡混合带，修计算，然后计算由于弥散所形成的过渡混合带，修 正按平均锋面运移的距离。正按平均锋面运移的距离。 55 共81页5 第二节第二节 近似解析解近似解析解 计算平均渗透锋面的运移距离计算平均渗透锋面的运移距离 ； 求出平均渗透锋面处污染物的浓度求出平均渗透锋面处污染物的浓度 ； 确定分界面弥散带确定分界面弥散带 和变形带和变形带 ，在平，在平 均渗透锋面距离上加上这两个带就可确定污染均渗透锋面距离上加上这两个带就可确定污染 水实际锋面的运移位置：水实际锋面的运移位置： nKmtLnDtL LLLL pd pd /cos/ )(5 . 0 L C

5、d L p L 相对密度差相对密度差 岩层倾角岩层倾角 66 共81页6 第二节第二节 近似解析解近似解析解 l 承压（或潜水）含水层承压（或潜水）含水层 污染地下水在承压含水层的运动，如果补给来源污染地下水在承压含水层的运动，如果补给来源 稳定，可视为似稳定运动，属于一维平面平行流，则稳定，可视为似稳定运动，属于一维平面平行流，则 渗透速度为：渗透速度为： 而实际流速为：而实际流速为： 常数 m q L HHK V e x )( 12 n V u x x 地下水的天然单宽流量地下水的天然单宽流量 L 77 共81页7 第二节第二节 近似解析解近似解析解 l承压（或潜水）含水层承压（或潜水）含

6、水层 e e e e CC CtxC C q mnL t t mn q xxL 0 0 ),( 初始时刻污染锋面的位置初始时刻污染锋面的位置 某时刻污染锋面的位置某时刻污染锋面的位置 天然地下水中污染物的浓度天然地下水中污染物的浓度 初始分界面以内，地下水中污染物初始分界面以内，地下水中污染物 的浓度的浓度 若是潜水，则取平均厚度代替上述各式中的若是潜水，则取平均厚度代替上述各式中的m。 L 88 共81页8 第三节第三节 数值解数值解 由于含水层以及污染物质迁移的复杂性，由于含水层以及污染物质迁移的复杂性， 解析解不易求得，通常采用数值解法。数解析解不易求得，通常采用数值解法。数 值法在剖分

7、后认为每一单元内均质同性，值法在剖分后认为每一单元内均质同性， 将方程线性化。将方程线性化。 有限差分法有限差分法 有限单元法有限单元法 边界元法边界元法 99 共81页9 一、一、 有限差分法有限差分法 有限差分法的有限差分法的基本思路基本思路是按照时间步长和空间步是按照时间步长和空间步 长将时间和空间区域剖分成若干网格，用未知函长将时间和空间区域剖分成若干网格，用未知函 数在网格节点上的值所构成的差分商数在网格节点上的值所构成的差分商近似代替近似代替所所 用偏微分方程中出现的各阶导数，从而把表示变用偏微分方程中出现的各阶导数，从而把表示变 量连续变化关系的偏微分方程离散为有限个代数量连续变

8、化关系的偏微分方程离散为有限个代数 方程，然后解此方程，然后解此线性代数方程线性代数方程，以求出污染物浓，以求出污染物浓 度在各度在各网格节点网格节点上不同时刻的解。上不同时刻的解。 10 10 共81页10 导数的差分近似导数的差分近似 设设f(x)为任一足够光滑的函数，把为任一足够光滑的函数，把f(x)沿沿x的正向的正向 展开为泰勒级数：展开为泰勒级数： 上式称为上式称为f(x)的的一阶向前差分一阶向前差分 3 33 2 22 ! 3! 2 )()( dx fdx dx fdx dx df xxfxxf )( )()( xR x xfxxf dx df x xfxxf dx df )()(

9、 11 11 共81页11 导数的差分近似导数的差分近似 如果把如果把f(x)对对x的负向展开为泰勒级数：的负向展开为泰勒级数： 上式称为上式称为f(x)的的一阶向后差分一阶向后差分 3 33 2 22 ! 3! 2 )()( dx fdx dx fdx dx df xxfxxf )( )()( xO x xxfxf dx df x xxfxf dx df )()( 12 12 共81页12 导数的差分近似导数的差分近似 把上面两个差分方程相减，可得：把上面两个差分方程相减，可得： 上式称为上式称为f(x)的的一阶中心差分一阶中心差分 同理可得同理可得f(x)的二阶导数的差分公式：的二阶导数的

10、差分公式： x xxfxxf dx df 2 )()( 22 2 )()(2)( x xxfxfxxf dx fd 13 13 共81页13 导数的差分近似导数的差分近似 若把上述一元函数若把上述一元函数f(x)换为浓度函数，并按网格换为浓度函数，并按网格 节点编号，则对节点编号，则对x的一阶向前、向后和中心差分的一阶向前、向后和中心差分 分别为：分别为： x CC x c x CC x c x CC x c n kji n kji nkji n kji n kji nkji n kji n kji nkji 2 | | | ,1,1 ),( ,1, ),( ,1 ),( 同样可以写出同样可以

11、写出 对时间对时间t的向的向 前、向后和中前、向后和中 心差分公式心差分公式 14 14 共81页14 导数的差分近似导数的差分近似 C(x,y,z,t)对对x的二阶差分公式为：的二阶差分公式为： 2 ,1,1 ),( 2 2 2 | x CCC x c n kji n kji n kji nkji 因此，在任何一个节点处，浓度因此，在任何一个节点处，浓度C(x,y,z,t) 的一阶和二阶偏导数均可用该节点及其的一阶和二阶偏导数均可用该节点及其 相邻节点上浓度值的线性组合表示。相邻节点上浓度值的线性组合表示。 15 15 共81页15 一维弥散方程的差分格式一维弥散方程的差分格式 设在无限含水

12、层中，存在一维均匀流场，其渗透设在无限含水层中，存在一维均匀流场，其渗透 速度速度V=nu，流动方向为，流动方向为x轴正向。则一维弥散方轴正向。则一维弥散方 程为：程为： Tt,0 Lx0 2 2 x c u x c D t c L 16 16 共81页16 一维弥散方程的差分格式一维弥散方程的差分格式 对时间区域对时间区域0,T和空间区域和空间区域0,L都作等距剖分，都作等距剖分， 设时间步长为设时间步长为t，空间步长为，空间步长为x，把第，把第i个节点个节点 xi处在处在tn时刻的浓度记为时刻的浓度记为Cin。 17 17 共81页17 一维弥散方程的差分格式一维弥散方程的差分格式 显示差

13、分格式显示差分格式 将将 及及 差分公式中的浓度取为差分公式中的浓度取为tn时刻的值，时刻的值， 便可得到其显示差分格式：便可得到其显示差分格式： x C 2 2 x C x CC u x CCC D t CC n i n i n i n i n i L n i n i 2 2 11 2 11 1 n i L n i L n i L n i C x tu x tD C x tD C x tu x tD C 1 2 2 1 2 1 ) 2 ( ) 2 1() 2 ( 18 18 共81页18 一维弥散方程的差分格式一维弥散方程的差分格式 显示差分格式的稳定条件显示差分格式的稳定条件 在应用显示差

14、分格式时，剖分的步长必须满足：在应用显示差分格式时，剖分的步长必须满足： 1 2 1 2 x tu x tDL 及 19 19 共81页19 一维弥散方程的差分格式一维弥散方程的差分格式 隐示差分格式隐示差分格式 将将 及及 差分公式中的浓度取为差分公式中的浓度取为tn 1时刻的 时刻的 值，便可得到其隐示差分格式：值，便可得到其隐示差分格式： x C 2 2 x C x CC u x CCC D t CC n i n i n i n i n i L n i n i 2 2 1 1 1 1 2 1 1 11 1 1 n i n i L n i L n i L CC x tu x tD C x

15、tD C x tu x tD 1 1 2 1 2 1 1 2 ) 2 ( ) 2 1() 2 ( 20 20 共81页20 一维弥散方程的差分格式一维弥散方程的差分格式 隐示差分格式的稳定条件隐示差分格式的稳定条件 在应用隐示差分格式时，剖分的步长必须满足：在应用隐示差分格式时，剖分的步长必须满足： x D u L 2 | 21 21 共81页21 一维弥散方程的差分格式一维弥散方程的差分格式 隐示差分格式的矩阵形式隐示差分格式的矩阵形式 1n CA 22 22 共81页22 一维弥散方程的差分格式一维弥散方程的差分格式 Crank-Nicolson差分格式差分格式 取显示差分格式和隐示差分格

16、式的平均，便可得取显示差分格式和隐示差分格式的平均，便可得 到到C-N差分格式：差分格式： ) 2 2 ( 2 1 ) 2 2 ( 2 1 1 1 1 1 2 1 1 11 1 11 2 11 1 x CC u x CCC D x CC u x CCC D t CC n i n i n i n i n i L n i n i n i n i n i L n i n i 23 23 共81页23 一维弥散方程的差分格式一维弥散方程的差分格式 Crank-Nicolson差分格式差分格式 整理后可得：整理后可得： n i L n i L n i L n i L n i L n i L C x tu

17、 x tD C x tD C x tu x tD C x tu x tD C x tD C x tu x tD 1 22 1 2 1 1 2 1 2 1 1 2 ) 2 ()1 ( 2) 2 ( ) 2 ()1 ( 2) 2 ( 1n AC 24 24 共81页24 追赶法：追赶法： i n ii n ii n ii CCC 1 1 11 1 方程的统一形式：方程的统一形式： 25 25 共81页25 追赶法：追赶法： 1 0 2 1 1 1 1 21 1 11 n nn C x D CC 引入边界条件：引入边界条件： 当当i=1时：时： 当当i=M时：时： MMMMMM M n MM n M

18、M n M n M n M CC CCC 2 2 11 1 1 1 11 1 其中 则 令 26 26 共81页26 追赶法：追赶法： 令令AC= M M n M M MM MMM C C C C CA 1 2 1 1 1 2 1 111 333 222 11 27 27 共81页27 追赶法：追赶法： 追赶法解方程组的方法：追赶法解方程组的方法： 先把对角矩阵先把对角矩阵A分解为：分解为：A=BW 其中：其中： MM MM ba ba ba b B 11 22 1 1 1 1 1 1 2 1 M w w w W 28 28 共81页28 追赶法：追赶法： 由于由于A=BW 1 112111

19、332333 221222 111 MMMM MMMMMM waba wbwaba wbwaba wbwaba wbb BWA 29 29 共81页29 追赶法：追赶法： 由此可得：由此可得： b1=1 b1w1=1 w1=1/b1 b2+a2w1=2 biwi=i wi=i/bi (i=1M-1) bi+aiwi-1=i bi=i-aiwi-1 bM=M-aMwM-1 30 30 共81页30 追赶法：追赶法： 由于由于A=BW B(WC)= 又令：又令：WC=Y 其中：其中：Y=(y1, y2, , yM)T BY= M M M M MM MM y y y y ba ba ba b B 1

20、 2 1 1 2 1 11 22 1 31 31 共81页31 追赶法：追赶法： 由以上等式可知：由以上等式可知： 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)/ bM 至此可以求出所有的至此可以求出所有的Y(y1, y2, , yM)T值。值。 32 32 共81页32 追赶法：追赶法： 又因为又因为WC=Y，即：，即： M M n M MM y y y y C C C C w w w CW 1 2 1 1 1 2 1 1 2 1

21、1 1 1 1 33 33 共81页33 追赶法：追赶法： 根据上述等式，可得：根据上述等式，可得： M n M M n MM n M n iii n i nn nn yC yCwC CwyCyCwC yCwC 1 1 1 1 1 1 1 1 1 2 1 32 1 2 1 1 21 1 1 34 34 共81页34 追赶法：追赶法： 追：追：根据根据BY=，求出，求出bi和和wi，并求出，并求出yi； 赶：赶：根据根据WC=Y，求出，求出Cin+1。 35 35 共81页35 二维弥散方程的差分格式二维弥散方程的差分格式 若多孔介质是各向同性的，流场为均匀的若多孔介质是各向同性的，流场为均匀的

22、 一维流场，渗透速度一维流场，渗透速度V=nu，流动方向取，流动方向取x轴轴 的正向，与流速垂直为的正向，与流速垂直为y方向，对流扩散方方向，对流扩散方 程可简化为：程可简化为： x c u y c D x c D t c TL 2 2 2 2 36 36 共81页36 二维弥散方程的差分格式二维弥散方程的差分格式 将研究区域网格化将研究区域网格化 tntyjcyxiax NttMcdyMabx nji , /,/ )(,/ )( 37 37 共81页37 二维弥散方程的差分格式二维弥散方程的差分格式 建立矩形网差分格式建立矩形网差分格式 l隐式差分格式隐式差分格式 t CC t c x CC

23、 x c n ji n ji n ji n ji , 1 , 1 ,1 1 ,1 2 2 1 1, 1 , 1 1, 2 2 2 1 ,1 1 , 1 ,1 2 2 2 2 y CCC y c x CCC x c n ji n ji n ji n ji n ji n ji 38 38 共81页38 二维弥散方程的差分格式二维弥散方程的差分格式 于是，内节点的差分方程为于是，内节点的差分方程为 x CC V y CCC D x CCC D t CC n ji n ji n ji n ji n ji T n ji n ji n ji L n ji n ji 2 2 2 1 ,1 1 ,1 2 1

24、1, 1 , 1 1, 2 1 ,1 1 , 1 ,1, 1 , 39 39 共81页39 二维弥散方程的差分格式二维弥散方程的差分格式 整理后可得整理后可得 n ji n ji T n ji T n ji L n ji TL n ji L CC y TD C y TD C x tV x tD C y tD x tD C x tV x tD , 1 1, 2 1 1, 2 1 ,1 2 1 , 22 1 ,1 2 ) 2 ( ) 22 1() 2 ( 40 40 共81页40 二维弥散方程的差分格式二维弥散方程的差分格式 边界条件边界条件 初始条件初始条件 ),(),( ),(),( 1,10

25、, 1,1,0 ni n Mini n i nj n jMnj n j tdxCCtcxCC tybCCtyaCC ),( 0 0 ,jiji yxCC 41 41 共81页41 二维弥散方程的差分格式二维弥散方程的差分格式 将网格取为正方形，即将网格取为正方形，即x=y=h，方程就，方程就 可简化为：可简化为： n ji n jiT n jiT n jiL n jiTL n jiL CTCD h TCD h CtVhtD h CtDtDh h CtVhtD h , 1 1, 2 1 1, 2 1 ,1 2 1 , 2 2 1 ,1 2 11 )( 2 1 )22( 1 )( 2 1 42 4

26、2 共81页42 二维弥散方程的差分格式二维弥散方程的差分格式 交替方向差分格式（交替方向差分格式（ADI） ADI法的实质是法的实质是将二维隐式差分格式化为一将二维隐式差分格式化为一 维的隐式差分格式维的隐式差分格式。 在第在第n层与第层与第n+1层两层正中间引入一个过层两层正中间引入一个过 渡层渡层n+1/2。 43 43 共81页43 二维弥散方程的差分格式二维弥散方程的差分格式 交替方向差分格式（交替方向差分格式（ADI） 分两步走：分两步走：在前半时间步长内，对在前半时间步长内，对x偏导用隐式偏导用隐式 表示，对表示，对y的偏导用显示表示；在后半时间步长的偏导用显示表示；在后半时间步

27、长 内，对内，对x偏导用显示表示，对偏导用显示表示，对y的偏导用隐式表的偏导用隐式表 示。示。 44 44 共81页44 二维弥散方程的差分格式二维弥散方程的差分格式 交替方向差分格式（交替方向差分格式（ADI） 前半时间步长：前半时间步长： 2 1,1,1,1 2 ,1,1, 2 2 2 2/ 2 1 2 1 2 1 2 1 2 1 2 1 y CCC D x CC V x CCC D t CC n ji n ji n ji T n ji n ji n ji n ji n ji L n ji n ji 45 45 共81页45 二维弥散方程的差分格式二维弥散方程的差分格式 交替方向差分格式（

28、交替方向差分格式（ADI） 后半时间步长：后半时间步长： 2 1 1, 1 , 1 1,1,1 2 ,1,1, 1 , 2 2 2 2/ 2 1 2 1 2 1 2 1 2 1 2 1 y CCC D x CC V x CCC D t CC n ji n ji n ji T n ji n ji n ji n ji n ji L n ji n ji 46 46 共81页46 二维弥散方程的差分格式二维弥散方程的差分格式 计算步骤：计算步骤： （1）根据初始条件或前一时段计算结果，可知）根据初始条件或前一时段计算结果，可知tn 时刻的浓度值，那么第一个方程中只有三个未时刻的浓度值，那么第一个方程中

29、只有三个未 知数（知数（n+1/2时刻）。如果固定时刻）。如果固定j值，对值，对i的联立的联立 方程组就是三对角线的（追赶法），即方程组就是三对角线的（追赶法），即沿沿x方向方向 追赶追赶，可求出过渡层（，可求出过渡层（n+1/2）层上的全部浓度）层上的全部浓度 值。值。 47 47 共81页47 二维弥散方程的差分格式二维弥散方程的差分格式 计算步骤：计算步骤： （2）根据前一步求出的过渡层（）根据前一步求出的过渡层（n+1/2）层上的）层上的 浓度值，再按同样的方法求出第二个方程中的浓度值，再按同样的方法求出第二个方程中的 （n+1）时刻的浓度值。即对于一个固定的）时刻的浓度值。即对于一个

30、固定的i值，值， 方程也是一个三对角方程组（追赶法），方程也是一个三对角方程组（追赶法），沿沿y方方 向追赶向追赶。 48 48 共81页48 有限单元法基本原理有限单元法基本原理 在在有限差分法有限差分法中，研究区域被一系列的网格节点中，研究区域被一系列的网格节点 所代替，在各个节点上以差分代替导数而得到求所代替，在各个节点上以差分代替导数而得到求 解相应偏微分方程的线性代数方程组。解相应偏微分方程的线性代数方程组。 有限单元法有限单元法是把研究区域剖分为有限个子区域，是把研究区域剖分为有限个子区域， 在每个子区域上用某种插值函数来近似待求解的在每个子区域上用某种插值函数来近似待求解的 未知

31、函数，从而得到求解相应偏微分方程的线性未知函数，从而得到求解相应偏微分方程的线性 代数方程组。代数方程组。 49 49 共81页49 边界元法基本原理边界元法基本原理 通过把求解区域的边界剖分为若干单元，将求函数通过把求解区域的边界剖分为若干单元，将求函数 的解简化为求节点的函数值，求解积分方程就化为的解简化为求节点的函数值，求解积分方程就化为 求解一组线性代数方程，在求出边界上的物理量以求解一组线性代数方程，在求出边界上的物理量以 后，计算区域内部的任一点的未知量可通过边界上后，计算区域内部的任一点的未知量可通过边界上 的已知量求出。的已知量求出。 与有限元法相比，具有降低维数、输入数据准备

32、简与有限元法相比，具有降低维数、输入数据准备简 单、计算工作量小、精度比有限元法高等优点。单、计算工作量小、精度比有限元法高等优点。 在处理非均质、非线性、非稳定等复杂问题时计算在处理非均质、非线性、非稳定等复杂问题时计算 效率低，不能将对流项化成沿边界的积分。效率低，不能将对流项化成沿边界的积分。 50 50 共81页50 第四节第四节 数理统计法数理统计法 地下水水质受多种因素综合影响，具有不地下水水质受多种因素综合影响，具有不 确定性，可用数理统计法预测水质变化；确定性，可用数理统计法预测水质变化； 主要包括：时间序列分析法、相关分析法、主要包括：时间序列分析法、相关分析法、 概率统计法

33、。概率统计法。 51 51 共81页51 （1）时间序列分析法）时间序列分析法 基本思路基本思路：分析时间序列的变化特征，选：分析时间序列的变化特征，选 择适当的模式建立预测模型，利用模型进择适当的模式建立预测模型，利用模型进 行趋势外推预测，对模型预测值进行评价行趋势外推预测，对模型预测值进行评价 和修正，得到预测结果。和修正，得到预测结果。 52 52 共81页52 （1）时间序列分析法）时间序列分析法-平均数平均数 平均数预测平均数预测：对动态数列中若干个时期的数值，：对动态数列中若干个时期的数值， 求其算术平均值或加权平均值，作为下一个时求其算术平均值或加权平均值，作为下一个时 期的预

34、测值。期的预测值。 n i in C n C 1 1 1 n i iin CW n C 1 1 )( 1 53 53 共81页53 （1）时间序列分析法）时间序列分析法-指数平滑指数平滑 一次指数平滑预测一次指数平滑预测：对动态数列中若干个时期：对动态数列中若干个时期 的数值，求其算术平均值或加权平均值，作为的数值，求其算术平均值或加权平均值，作为 下一个时期的预测值。下一个时期的预测值。 2 2 1 )1( )1 ()1 ( tttt CaaCaaaCS )1( 1 21 )1( )1 ( )1 ()1 ( tt tttt SaaC CaaaCaaCS 54 54 共81页54 （1）时间序

35、列分析法）时间序列分析法-指数平滑指数平滑 一次指数平滑法是以最近周期的一次指数平滑一次指数平滑法是以最近周期的一次指数平滑 值作为下一周期的预测值，即：值作为下一周期的预测值，即： 其中平滑系数其中平滑系数a的选择直接影响预测效果。的选择直接影响预测效果。如如 果时间序列长期趋势稳定，应取较小的果时间序列长期趋势稳定，应取较小的a值；值； 若时间序列不稳定，应取较大的若时间序列不稳定，应取较大的a值。值。 )1( 1 )1( 1 )1 ( tttt SaaCSC 55 55 共81页55 （1）时间序列分析法）时间序列分析法-指数平滑指数平滑 a越小，对数据的平滑能力越强，但对数据变化的敏感

36、越小，对数据的平滑能力越强，但对数据变化的敏感 性越差；性越差；a越大，对数据的平滑能力越差，但对数据变越大，对数据的平滑能力越差，但对数据变 化的敏感性越强。化的敏感性越强。 56 56 共81页56 （1）时间序列分析法）时间序列分析法-指数平滑指数平滑 初始值初始值S0(1)对预测结果的影响很大，一般可取第一个数对预测结果的影响很大，一般可取第一个数 据（据（x1）作为初始值；若总体数据量较少，则可取最）作为初始值；若总体数据量较少，则可取最 初几个数据的平均值作为初始值。初几个数据的平均值作为初始值。 57 57 共81页57 （1）时间序列分析法）时间序列分析法-指数平滑指数平滑 二

37、次指数平滑预测二次指数平滑预测：对一次指数平滑值再做一：对一次指数平滑值再做一 次指数平滑：次指数平滑： )2( 1 )1()2( )1 ( ttt SaaSS 58 58 共81页58 （1）时间序列分析法）时间序列分析法-指数平滑指数平滑 二次指数平滑值不直接用于预测，而是根据滞二次指数平滑值不直接用于预测，而是根据滞 后偏差建立线性预测模型：后偏差建立线性预测模型： 其中：其中： TbaC ttTt )2()1( 2 ttt SSa 1 )2()1( ttt SS a a b 59 59 共81页59 （1）时间序列分析法）时间序列分析法-指数平滑指数平滑 若一、二次指数平滑值不满足预测

38、要求，应采若一、二次指数平滑值不满足预测要求，应采 用三次指数平滑法建立非线性预测模型：用三次指数平滑法建立非线性预测模型： (3) 1 )2()3( )1 ( ttt SaaSS 60 60 共81页60 （1）时间序列分析法）时间序列分析法-指数平滑指数平滑 三次指数平滑法建立的非线性预测模型：三次指数平滑法建立的非线性预测模型： )3()2()1( 333 tttt SSSa )34()810(2)56( )12( )3()2()1( 2 tttt SaSaSa a a b 2 T tttTt cTbaC 2 )12( )3()2()1( 2 tttt SSS a a c 61 61 共

39、81页61 （2）相关分析法相关分析法 线性相关预测线性相关预测：建立直角坐标下数据散点图，：建立直角坐标下数据散点图， 用一元线性回归方程近似描述数据点的分布规用一元线性回归方程近似描述数据点的分布规 律：律： bxay 62 62 共81页62 （2）相关分析法相关分析法 根据最小二乘法根据最小二乘法(要求残差要求残差RSS最小最小)可求得：可求得： 用相关系数用相关系数r来判断回归方程实用性：来判断回归方程实用性： 22 22 2 )( )( ii iiii ii iiiii xnx yxnyx b xnx xyxyx a 22 )()( )( yyxx yyxx r ii ii 2 )

40、( ii bxayRSS 63 63 共81页63 （2）相关分析法相关分析法 曲性相关预测曲性相关预测：将直角坐标下散点图拟合形成：将直角坐标下散点图拟合形成 的有规律的曲线进行函数转换，变为前述的一的有规律的曲线进行函数转换，变为前述的一 元线性回归方程进行计算。元线性回归方程进行计算。 如：如： b axy bXaY yYxX lg lg,lg 则 令 64 64 共81页64 国外地下水模拟软件的发展现状与趋势国外地下水模拟软件的发展现状与趋势 利用数值模型对地下水流和溶质运移问题进行利用数值模型对地下水流和溶质运移问题进行 模拟的方法以其有效性、灵活性和相对廉价性模拟的方法以其有效性

41、、灵活性和相对廉价性 逐渐成为地下水研究领域的一种不可或缺的重逐渐成为地下水研究领域的一种不可或缺的重 要方法，并受到越来越大的重视和广泛的应用。要方法，并受到越来越大的重视和广泛的应用。 一个完整的地下水模拟过程包含一个完整的地下水模拟过程包含3个部分：个部分：前处前处 理、模型计算和后处理理、模型计算和后处理。 65 65 共81页65 国外地下水模拟软件的发展现状与趋势国外地下水模拟软件的发展现状与趋势 前处理前处理是指在进行模拟计算之前对计算过程中是指在进行模拟计算之前对计算过程中 所需数据的整理、组织、输入及计算网格的编所需数据的整理、组织、输入及计算网格的编 号与生成。号与生成。

42、模型计算模型计算是进行地下水流动或水质运移正反演是进行地下水流动或水质运移正反演 计算，常用的方法主要有：有限差分法、有限计算，常用的方法主要有：有限差分法、有限 元法、边界元法等。元法、边界元法等。 后处理后处理是将计算所产生的结果数据，用图形或是将计算所产生的结果数据，用图形或 表格显示或存放起来，以供研究人员方便地进表格显示或存放起来，以供研究人员方便地进 行分析和使用。行分析和使用。 66 66 共81页66 国外地下水模拟软件的发展现状与趋势国外地下水模拟软件的发展现状与趋势 通过近二十年的研究与发展，国际上已经形成通过近二十年的研究与发展，国际上已经形成 了一批非常有影响的地下水模

43、拟的软件，它们了一批非常有影响的地下水模拟的软件，它们 今天在国际地下水模拟研究领域依旧非常活跃，今天在国际地下水模拟研究领域依旧非常活跃， 如如MODFLOW、FEFLOW、MT3DMS、 MT3D99、PEST、MODPATH、UCODE等。等。 67 67 共81页67 MODFLOW MODFLOW是由美国地质调查局的是由美国地质调查局的McDonald 和和Harbaugh于于80年代开发出来的一套专门用年代开发出来的一套专门用 于孔隙介质中于孔隙介质中三维有限差分三维有限差分地下水流数值模拟地下水流数值模拟 的软件。的软件。 l程序结构的程序结构的模块化模块化（模拟抽水引起地面沉降

44、的子（模拟抽水引起地面沉降的子 程序包、模拟水平流动障碍的子程序包等）。程序包、模拟水平流动障碍的子程序包等）。 l离散方法的离散方法的简单化简单化（有限差分法（有限差分法36036018个个 网格单元）。网格单元）。 lMODFLOW引进了引进了应力期应力期(Stress Period)概念。概念。 l求解方法的多样化（强隐式法、逐次超松弛迭代求解方法的多样化（强隐式法、逐次超松弛迭代 法、预调共轭梯度法）。法、预调共轭梯度法）。 68 68 共81页68 MT3D99 MT3D99是美国是美国S.S.Papodopulos & Associates 公司的公司的Zheng博士设计的模拟博士

45、设计的模拟三维地下水溶质三维地下水溶质 运移运移程序程序MT3D(1990)的增强版，成为目前世的增强版，成为目前世 界上首屈一指的溶质运移模拟软件。界上首屈一指的溶质运移模拟软件。 lMT3D99能够模拟地下水系统中的平流、扩散、衰能够模拟地下水系统中的平流、扩散、衰 减、溶质化学反应、线性与非线性吸附作用等现象，减、溶质化学反应、线性与非线性吸附作用等现象， 能够对承压含水层，不承压含水层，承压与不承压能够对承压含水层，不承压含水层，承压与不承压 交替的含水层以及倾斜的和单元厚度变化的含水层交替的含水层以及倾斜的和单元厚度变化的含水层 进行空间离散。进行空间离散。 lMT3D99提供了丰富

46、的求解方法（标准有限差分法、提供了丰富的求解方法（标准有限差分法、 基于基于Eulerian-Lagrangian的粒子跟踪方法和高阶有的粒子跟踪方法和高阶有 限体积限体积TVD方法）。方法）。 69 69 共81页69 PEST2000 PEST2000是由澳大利亚是由澳大利亚Watermark Computing公公 司开发的司开发的PEST新版本，是功能强大的模型独立的新版本，是功能强大的模型独立的 参数估计参数估计程序程序。 lPEST2000利用一个强有力的数值反演算法来利用一个强有力的数值反演算法来 “控制控制”运行中的模型，程序在每次模拟之后运行中的模型，程序在每次模拟之后 自动调整所选择的模型参数，直到将校