（固体力学专业论文）土柱试验模型及参数反演.pdf

山东理t 人学硕i ：学位论文 摘要 摘要 土柱试验是研究土壤和地下水中污染物迁移转化规律的重要手段土柱试验的主要 目的是研究包气带岩土特性、厚度及其结构，通过数值方法等手段测定土层的含水性及 溶质迁移转化参数，建立符合实际情况的溶质迁移模型，进而用以研究相同水文地质条 件下的土壤及地下水污染问题 本文基于土壤及地下水中溶质运移理论，以土柱试验为背景，主要对一维土柱试验 模型及反问题进行研究和讨论，包括试验机理的分析、试验模型的构建、数值求解及参 数反演算法及应用问题其主要工作是： 1 ) 通过引进不动点同伦方法的思想，建立了一种改进的最佳摄动量算法，并应用于 以硫酸盐为主的多组分溶质运移模型的参数反演问题，验证了改进的最佳摄动量算法在 扩大收敛范围方面的优越性 2 ) 完成了一个实验室扰动土柱试验，通过对试验数据的分析，建立了一个带有时间 相关源( 汇) 项的溶质运移模型，来描述硫酸根离子、钾离子等八大离子在土层中的迁 移行为过程，并应用改进的最佳摄动量迭代算法对模型中的未知参数进行了有效的数值 反演，反演计算结果与实际观测数据基本吻合 3 ) 考察了孝妇河一张冉原状土柱试验，将溶质的水动力弥散理论与水化学分析结合 起来，建立了一个基本符合实际情况的数学模型，并应用最佳摄动量算法探讨了多参数的 联合反演问题，数值模拟结果较好 关键词：土柱试验；同伦方法；溶质运移模型；最佳摄动量算法；反演算法；数值模拟 山东理t 大学硕 ：学何论文 a b s t r a c t a b s t r a c t a ni m p o r t a n tm e t h o df o rd e s c r i b i n ga n dr e s e a r c h i n gs o l u t et r a n s p o r ti nt h es o i la n d g r o u n d w a t e ri sb yd o i n gs o i l c o l u m ne x p e r i m e n t t h et a r g e to ft h es o i l c o l u m ne x p e r i m e n ti s t os t u d yt h eg e o t e c h n i c a lp r o p e r t y , t h i c k n e s sa n ds t r u c t u r eo ft h es o i l ，d e t e r m i n et h ea q u o s i t y a n dt h es o l u t et r a n s p o r t a t i o np a r a m e t e r s ，a n dc o n s t r u c tat r a n s p o r tm o d e lw h i c hs a t i s f i e st h e p r a c t i c a ls i t u a t i o n b a s e do nt h es o l u t et r a n s p o r tt h e o r y , t h i sp a p e ri st os t u d ys o l u t et r a n s p o r ta c t i o n si nt h e s o i l sa n dg r o u n d w a t e ri no n ed i m e n s i o nb ya n a l y s z i n gt w ok i n g so fs o i l - c o l u m ne x p e r i m e n t s t h em a i ns t u d ye n d e a v o r st od ot h ef o l l o w i n gt h r e et h i n g s f i r s t l y , a ni m p r o v e do p t i m a lp e r t u r b a t i o n i t e r a t i o na l g o r i t h mw a sc o n s t r u c t e db y r e c o m m e n d i n gah o m o t o p ym e t h o d t h ei m p r o v e di n v e r s i o na l g o r i t h mw a st e s t i f i e dt oh a v e w i d ec o n v e r g e n tr e g i o n sb yp e r f o r m i n go na ni n v e r s ep r o b l e mo fd e t e r m i n i n gp a r a m e t e r si n m u l t i - c o m p o n e n ts o l u t e st r a n s p o r to fs u l p h a t e s e c o n d , ad i s t u r b e ds o i l - c o l u m ne x p e r i m e n tw a sc o m p l e t e d ，a n ds o l u t e st r a n s p o r tm o d e l d e s c r i b i n gt r a n s f e rb e h a v i o r so fe i g h ti o n s 谢mn o n l i n e a rs o u r c e ( s i n k ) t e r m sr e l a t i n gw i t ht i m e w a se s t a b l i s h e db a s e do nt h ea n a l y s i so ft h ee x p e r i m e n t a lp r o c e s s f u r t h e r m o r e ，t h ei m p r o v e d o p t i m a lp e r t u r b a t i o ni t e r a t i o na l g o r i t h mw a sa p p l i e dt od e t e r m i n eu n k n o w np a r a m e t e r si nt h e m o d e l ，a n dt h en u m e r i c a li n v e r s i o nr e s u l t sb a s i c a l l yc o i n c i d ew i t ht h eb r e a k t h r o u g hd a t a f i n a l l y , a nu n d i s t u r b e ds o i l - e o l u n me x p e r i m e n tw a si n v e s t i g a t e d , a n dam a t h e m a t i cm o d e l w h i c hs a t i s f i e st h ep r a c t i c a ls i t u a t i o nw a sc o n s t r u c t e db yc o m b i n i n gh y d r o d y n a m i cd i s p e r s i o n t h e o r yw i t hh y d r o c h e m i s t r ya n a l y s i s n u m e r i c a ls i m u l a t i o nr e s u l t ss h o wt h a tt h eo p t i m a l p e r t u r b a t i o ni t e r a t i v ea l g o r i t h mi se f f e c t i v ea tl e a s tt ot h ei n v e r s ep r o b l e mo fm u l t i p a r a m e t e r i d e n t i f i c a t i o na r i s i n gi nt h eu n d i s t u r b e ds o i l - c o l u n me x p e r i m e n t k e yw o r d s ：s o i l - c o l u m ne x p e r i m e n t ；h o m o t o p ym e t h o d ；s o l u t et r a n s p o r tm o d e l ；o p t i m a l p e r t u r b a t i o ni t e r a t i o na l g o r i t h m ；i n v e r s i o na l g o r i t h m ；n u m e r i c a ls i m u l a t i o n u 独创性声明 本人声盟所呈交的论文是我个人在导师指导下进行的研究工作及取得的研究成果。 尽我所知，除了文中特别加以标注和致谢的地方外，论文中不包含其他人已经发表或撰 写过的研究成果，也不包含为获得山东理工大学或其它教育机构的学位或证书而使用过 的材料。与我一同工作的同志对本研究所做的任何贡献均已在论文中作了明确的说明并 表示了谢意。 研究生妣瓿锪时间：二吖年月石日 ” 关于论文使用授权的说明 本人完全了解山东理工大学有关保留、使用学位论文的规定，即：学校有权保留送 交论文的复印件和磁盘，允许论文被查阅和借阅；学校可以用不同方式在不同媒体上发 表、传播学位论文的全部或部分内容，可以采用影印、缩印或扫描等复制手段保存、汇 编学位论文。 ( 保密的学位论文在解密后应遵守此协议) 研究生签名： 甄镶 时间：2 哆年月乡日 导师签名：弓遂 时间： 每 跫拥 学位论文出版授权书 本人完全同意中国优秀博硕士学位论文全文数据库出版章程( 以下简称 “章程”) ，愿意将本人的学位论文提交“中国学术期刊( 光盘版) 电子杂志社 在中国优秀博硕士学位论文全文数据库中全文发表。中国优秀博硕士学位 论文全文数据库可以以电子、网络及其它数字媒体形式公开出版，并同意编入 中国知识资源总库，在中国博硕士学位论文评价数据库中使用和在互联 网上传播，同意按“章程”规定享受相关权益。 导师签名：氛咐汤如咐 丝生竺 ： 一 名 年 隧 堕 作 文一 山东理t 人学硕 ：学位论文第一章绪论 1 1 研究背景 第一章绪论 我国是一个水资源并不丰富的国家，水资源人均占有量只占世界平均数的1 4 左右 随着社会经济的迅速发展和人口的不断增长，水的供需矛盾日益突出地下水作为水资 源的重要组成部分，是举足轻重的供水源和生态系统的重要支撑，尤其在我国北方，饮 用水基本由地下水供给但是，人类活动及工业化进程的负面效应造成了相当严重的水环 境污染，对我国国民经济的持续发展和社会稳定都产生了较为严重的影响已有调查显 示，土壤及地下水环境污染已呈现由点向面、由城市向农村扩展的趋势，水环境恶化和 水质污染已到了极为严重的程度这使得对污染控制与修复的理论、方法及技术创新的 研究变得极为迫切 作为土壤一水环境系统中的研究热点，污染物质在土壤介质中的运移模型问题一直没 有得到很好的解决由于土壤是一个极其复杂的体系，它不仅包括固相、液相和气相， 而且各相又具有有机和无机的组分，这些特点影响着土壤的物理化学性质，对土壤中污 染物溶质的运移具有直接的影响污染物进入土壤，将会发生水文的、物理的、化学的 以及生物的反应物理作用过程有对流、弥散、扩散、压缩、衰变( 放射性污染物) 等化 学反应包括水相络合作用、表面络合作用、氧化一还原作用、酸碱反应、离子交换、吸附 作用、以及沉淀一溶解作用，有机物质还会发生降解反应这些作用直接影响着污染物 在土壤中的运移行为 显然，搞清楚污染物溶质在土壤中的吸附、水解、溶解、离子交换以及对流弥散等 物理化学过程，对于制定土壤及地下水污染的防治对策具有重要的指导意义对于溶质 运移行为的研究，常用的数学模型是对流弥散方程，该方程考虑了对流、分子扩散、机 械弥散、吸附以及源( 汇) 项等的作用，是研究溶质在土壤及地下水中迁移行为的基本 模型不过，由于土壤及地下水介质具有不可直接测量的的特点，模型参数一般很难直 接得到因而，依据某些附加观测数据，借助于数值方法来反演确定难以直接获得的模型 参数，就构成了土壤及地下水中溶质运移的参数识别反问题随着计算机的发展和计算 方法的进步，反演方法在水文地质等问题的研究中得到了广泛的应用 污染物在土壤及地下水中的迁移转化是一个复杂的过程，要真正掌握其迁移转化机 理必须借助于渗流试验研究土柱渗流试验是研究土壤和地下水中污染物迁移转化规 律的重要手段i l j 通过试验，利用数值反演方法等手段可以确定土层的含水性及溶质迁 移转化参数，建立符合实际情况的溶质运移模型，进而研究相同水文地质条件下的土壤 及地下水的污染问题 t l l 东理t 大学硕l ：学位论文 第。章绪论 1 2 研究现状及其进展 国内外学者对污染物在土壤及地下水中迁移转化规律的研究高度重视，且不断取得 新的成果近五十年来，土壤溶质运移理论及应用研究取得了很大进展从上世纪5 0 年代 开始，t a y l o r ，b e a r ，n i e l s e n ，b i g g e r 和b r e a s l e r 等从试验和理论上进一步说明了土 壤溶质运移过程中对流、扩散和化学反应的祸合性质，开创了应用数学模型来说明和解 释溶质运移过程的新局面 t a y l o r ( 1 9 5 3 ) 首先定量研究了毛管中的纵向弥散作用；a x i s ( 1 9 5 6 ) 将t a y l o r 的方法 应用于不规则形状的毛管，认为局部的速度分布不是抛物线的上世纪六十年代初期，土 壤学家n i e l s e n 和b i g g e r 将石油科学中的“混合置换理论”原理应用到土壤科学领域， 推动了混合置换理论在农业生产实践中的应用研究溶质运移中的混合置换现象实际上 是由对流、分子扩散和机械弥散三个物理过程以及溶质在运移过程中所发生的众多化学、 物理化学过程和其他过程综合的结果随后n i e l s e n 和b i g g e r 从理论上推导并建立了对 流一弥散方程，较好地描述了保守性溶质在多孔介质中的迁移规律n i e l s o n 首次系统论 述了c d e ( c o n v e c f i o nd i s p e r s i o ne q u a t i o n ) 的科学性和合理性从n i e l s e n 建立c d e 方程以 来，该方程一直应用到现在，被称为土壤溶质运移理论研究的经典方程，在溶质运移研 究史上，建立了第一个丰碑 以色列学者j b e a r 分别于1 9 7 2 与1 9 7 9 年发表的多孔介质流体运动领域的经典著 作d y n a m i c so ff l u i d si np o r o u s m e d i a 和h y d r a u l i c so fg r o u n d w a t e r ，极大 推动了土壤及地下水溶质迁移理论的研究以多孔介质中的物质输运理论为基础，以计算 机为工具建立溶质的数值模型来预报污染物质的扩散趋势，已经成为环境保护、城市规 划、污染治理等领域的必要手段 描述土壤及地下水中溶质迁移行为特征的模型不外乎两类，一类是溶质的水动力迁 移( 对流一弥散) 模式，这种模式只考虑水流对溶质迁移的影响，把溶质的迁移行为看作是 随水的流动，基本不考虑水化学因素的作用；另一类是溶质的多组分水化学平衡模式， 这种模式研究的是在各种复杂的水化学作用下溶质的形态分布，几乎没有考虑水动力作 用对溶质迁移的影响事实上，溶质在地下介质中的运动除了受地下水中对流一弥散作用 影响外，还受到地下水多组分系统中的络合作用、吸附一解吸作用、溶解一沉淀作用、氧 化一还原作用和酸一碱作用的制约仅仅考虑水动力因素或水化学因素的溶质迁移模式都 是片面的，都无法真实地反映客观世界的真实情况，只有将两者合理地耦合起来，才能 进一步揭示多孔介质中溶质的运动规律，从而为土壤及地下水污染的预报、防治提供依 据【2 1 土柱试验是研究土壤和地下水中污染物迁移转化规律的一个重要手段为了揭示土 壤溶质运移的客观机理和运移规律，进行了大量的室内外土柱试验。如澳大利亚土地研 究中心的r o s e f 3 1 ( 1 9 7 1 ) 在室内一维土柱上进行了溶质运移的试验研究，测出了溶质穿透 2 山东理t 人学硕 ：学位论文 第。章绪论 曲线；美国加州大学j u r y l 4 1 等( 1 9 7 4 ) 在砂土中拌盐，用灌水入渗研究溶质运移问题；美 国国家盐改实验室v a ng e n u c h t e n 5 j 进行了一系列室内土柱试验后，系统的论述了室内土 柱试验的边界条件问题；1 9 8 6 年，n i e l s e n 等 6 1 从试验和理论两方面阐述了土壤溶质运 移过程中对流、扩散和化学反应的耦合性质，并给出了说明和解释溶质运移过程的数学 模型 随着计算技术的进步，数值方法与计算软件在土柱试验及溶质运移研究中得到广泛 应用k e n n e t h 7 以钼元素为示踪剂在砂砾石潜水层中进行迁移试验，建立了反应性溶质 运移模型并进行了数值模拟；p a n g 掣8 l 通过b a t c h 实验和土柱试验，利用c x t f i t 软件研 究计算了镉离子在沙砾层中的非平衡运移行为国内对于土柱试验的研究相对较晚，主要 侧重于应用方面如史海滨等【9 j 通过吸附实验与水动力弥散实验建立了吸附模型与不流 动水模型，并应用土柱试验讨论了可动区域的溶质运移模型；王东海等【l o l 应用土柱试验 讨论了石油类污染物在砂砾石层中的迁移行为；崔凯等【l l 】通过微机控制的土柱试验装 置，研究了镉离子在非饱和土中的迁移转化规律，并利用梯度正则化算法对模型参数进 行了估算；李功胜等【1 2 塘过对一个污水淋滤原状土柱试验数据的分析，获得了描述硫酸 根离子在酸性土壤中的迁移模型，并对模型中的源项系数进行了有效的数值反演；马昱 等 2 1 考察了一个矿坑水淋滤原状土柱试验，根据水化学分析与水动力弥散方程，建立了 一个描述多组分溶质运移的数学模型，并提出了同时反演多个参数的最佳摄动量迭代算 法 可以看出，研究土柱试验的关键之一在于数学模型的构建与未知参数的反演目前， 虽然有很多土柱试验的数学模型，但一个具体的土柱试验模型的构建仍然比较困难这对 应着对实际问题所发生一系列复杂的化学、物理变化机理的认识和把握，或者说是对试 验土层中发生的吸附、溶解、离子交换、对流弥散、氧化一还原等物理化学作用的理解和 描述在解决了试验机理及其模型的构建之后，应用已有软件直接对问题进行分析计算和 讨论在对问题进行讨论时存在一个问题，就是未知参数的确定问题，属于一类对流弥散 一反应扩散方程( 组) 的系数反演问题其核心是使得模型计算值与土柱实验观测值尽量 吻合，或者使计算值与穿透数据在某种意义上最接近这种解决问题的途径在数学上就是 数理方程反问题的数值解法 通过以上分析，我们应该更加关注土柱试验反问题理论的探讨，加大对模型机理的 研究，将水动力弥散模式与水化学作用结合起来，建立耦合模型；不仅要对模型进行数 值模拟研究，而且对于模型中难以直接测量得到的模型参数，探讨应用反演算法进行优 化确定这将是今后对于土壤及地下水溶质运移行为研究的主要方面 1 3 本文所做工作 本文依据两类土柱试验数据，探讨了土壤及地下水中溶质运移模型及参数反演问题 3 l i i 东理t 大学硕i j 学位论文 第一币绪论 主要丌展了试验机理分析、数学模型构建、模型求解以及参数反演算法等方面的研究工 作所做的主要工作是： ( 1 ) 基于不动点同伦方法对一般最佳摄动量算法进行改进，提出了一种对初值依赖性 较小的最佳摄动量算法，并应用于一个土柱试验模型参数的反演问题，获得了较好的反演 结果 ( 2 ) 完成了一个扰动土柱试验，通过对试验过程及试验机理的分析，建立了一个描述 土层中主要组分对流弥散过程的数学模型，并应用最佳摄动量迭代算法对模型中的参数 进行了数值反演求解，反演重建结果与实际观测数据基本吻合 ( 3 ) 考察了一个用矿坑水淋滤酸性土壤的原状土柱试验，将水动力弥散模式与水化学 分析结合起来，建立了一个描述土层中各组分运移行为的耦合模型，并应用最佳摄动量 算法对模型中的多个未知参数进行了数值反演模拟，模拟结果较好 4 i ij 东珲t 大学硕t 学位论文 第二章模型数值解法及j 反演算法的改进 第二章模型数值解法及其反演算法的改进 2 1 模型数值解法 如前所述，对流一弥散方程是揭示溶质在多孔介质中运移行为i b j 9 1 的基本数学模型 不过，在建立具体区域中的溶质运移模型时，还应该综合考虑对流、弥散、离子吸附与交 换以及不同化学组分之间可能发生的化学反应等各种因素的影响通常，一维对流一弥散 方程表示为： r t a u ( x , t ) - d 挈一y t o u ( x , t ) 州咖) ( 2 1 ) 西 融2融 、7 ( z 1 j ( o x ，0 t t t 0 1 ) 其中：u ( x ，f ) 表示为t 时刻、x 处的溶质浓度，为溶质运移区域的长度，为所考虑的 运移时间； r 为阻滞因子，其值与土样密度成，吸附系数k d 以及体积含水率秒等相关，可表示 为 欠：l + 华； ( 2 2 ) d 为弥散系数，v 为污染物穿过土层的达西流速； r ( x ，f ；甜) 是源( 汇) 项，它反映了溶质运移中可能发生的物理、化学作用，其值的变 化与运移空间、时间及浓度等相关 由泛定方程式( 2 1 ) ，再结合适当的初边值条件，即构成一个数理方程定解问 题，可用于研究多孔介质中溶质的运移行为过程对于方程( 2 1 ) 结合初边值条件构 成的定解问题可以用m a t l a b 计算软件自带的偏微分方程工具箱中的p d e p e 命令 函数求解下面给出应用p d e p e 函数求解方程( 2 1 ) 的具体方法【2 吻1 p d 印e 用于求解如下形式的一维椭圆形或抛物型方程的初边值问题 c ( 蹦双詈= x 1 昙( 妒厂( 列双罢) ) + s ( 列双面o u )( 2 3 ) t o ，t ，口x b 初值条件，= t ou ( x ，t o ) = u o ( x ) ( 2 4 ) 边界条件石：b p ( x ，z f ) + g ( x ，) 厂( x ，f ，甜，争：0 ( 2 5 ) 其中 a ，b 为有限区域，a 0 ( 五f ，z f ，罢) 是通量项，j ( 墨，甜，i o u ) 为源项 i i l 东理t 大学硕i j 学位论文 第二章模型数值解法及其反演算法的改进 p d e p e 命令调用格式为： s o l2p d e p e ( m ，p d e f u n ，i c f u n ，b c f u n ，x m e s h ，t s p a n ，o p t i o n s ) 各参数含义分别为： m m = o ，对应平面对称性；m = l ，为柱对称性；m - - 2 ，是球对称性； l x i e f u n 定义微分方程形式的调用函数，其调用格式为 【c ，f ，s = p d e f u n ( x ，t ，u ，d u d x ) 参数c 、f 、s 都是列矢量，为微分方程函数p d e f u n 的各项系数； i c f u n 定义初始条件函数文件，其形式为 u o = i c f u n 通过返回参数( u 0 ) ，来确定初始时刻的取值； b c f u n 定义边界条件函数文件，其格式为 嘲，q l ，p r ，q r = b c f u n ( x l ，u l ，x r ，u r ，0 对于边界条件，只有p 与u 有关p 、q 在a 、b 两点处的取值将在边界条件( 2 5 ) 中给与 赋值确定 x m e s h 矢量【秭，五，】是空间变量x 的网格划分取值，而 五 且元素数 目3 ： t s p a n 矢量，f l ，厶】是时间变量f 的网格划分取值，要求气 f l 0 e u ( x ，0 ) = g ( x ) ，x q ( 2 1 2 ) i b u ( x ，) = f ( x ，f ) ，x 孢，f 0 其中z f 为一向量函数，三为矩阵微分算子，b 为矩阵边值条件算子，e 为矩阵初值条件 算子，c ( x ) 为未知的待定向量函数，三线性依赖于c ( x ) ，q 为空间区域，a q 为q 的边 界 该类问题的反演问题为：从( 2 1 2 ) 及附加条件 材( 堋i 耐2 烈x ，) ( 2 1 3 ) ( x a q ca q ，0 t 丁) 1 3 山东理t 大学顾i j 学位论文 第章模型数位解泫及其反演算法的改进 来确定算子三中的未知向量函数f ( x ) 该类反演问题很容易建立一个求下面非线性算子方程的解： 彳p ( x ) 】= 缈( z ，) ( 2 1 4 ) 利用正则化方法可将反演问题转化为下列非线性泛函 以( c ) = l l 么【c 】一伊i 巳( 削【o ，r 】) + e t d ( c ) ( 2 1 5 ) 的极值问题，其中口为正则化参数，d 为l 2 ( t a ) 上的稳定泛函 利用数值方法求解非线性泛函( 2 1 5 ) 的极小问题，便可求得反演问题的数值解最佳 摄动量算法是求解上述极小问题的一种方法最佳摄动量算法根据算子识别的摄动法、线 性化技术和函数逼近论提出的一种数值迭代方法，通过给定未知参数经验值，求解正问 题得到与观测值对应的计算值，根据计算值与观测值误差最小的思想来求解最佳扰动量 配( x ) ，进而确定出勺“； c n + 1 ( x ) = c ；( x ) + 万厶( x ) ( 2 1 6 ) 其中，摄动量。( x ) 由下列非线性最优化问题来确定： 以( 万) 爿l 么【厶+ 眈卜缈| | ( 触，x l o ，d + 口d ( 万) ( 2 1 7 ) 下面根据上节同伦方法的基本思想，将( 2 1 7 ) 做不动点同伦修正，构造新的最佳 摄动量算法定义泛函 以( 万乌) = ( 1 一允) l l 么【c 。+ 万巳卜缈l | ( m o ，r d + 允d ( 万巳) ( 2 1 8 ) 同伦参数允( 0 ，1 ) 的选取可以参考人工神经网络的传递函数，即 1 允= a 伽) = 石, 枭, - j - - p t v 0 - - n j ， ( 2 l9 ) 1 其中，为下降速率参数，0 为初值选取参数，玎为迭代步数 一般地，传递函数满足三条原则：该函数介于( o ，1 ) 之间，且随着迭代地进行不断减 小；该函数是连续光滑的，这样可使得同伦参数稳定的下降；随着迭代的深入，五会 缓慢的向0 靠近，但不等于零，这样可以在迭代结束时使名伽) 停留于一个较合适的位置 对于参数卢，若采用较大的值，可以增大同伦参数的修正步长，加速求解的收敛速度， 但过大的值可能会引起迭代的发散对于另外一个参数o ，取值越大，则同伦参数的 初始估计值越接近1 对于本文的问题，取值为0 5 ，0 选取3 当名为o 时，泛函( 2 1 8 ) 转化为拟解法求解极小化问题从问题( 2 1 8 ) 出发，在求解 过程中，按照某种方式不断减小旯的数值，最后可以通过一个合理而又稳定的求解路径 获得反问题的近似解，通过这种方式，可以有效的扩大原问题的计算收敛区域 假设方程( 2 1 2 ) 反演问题的精确解为c 。( x ) ，相对应于原方程定解问题的解为 “。( x ，f ) c ( x ) 为线性完备实函数空间k 上的一个元素，足上的一个基函数族为 1 4 山东理e 大学形! f j 学位论义筇币筏掣数值解法及j e 及演耳法的改进 办( x ) ，晚( x ) ，则 c + ( x ) = f 仍( 力 ( 2 2 0 ) 取有限项进行逼近，得 c ( x ) f 仍( x ) ( 2 2 1 ) 疗的大小取决于逼近精度的要求因此，这类问题就是确定一个n 维实向量 k 7 = ( 毛，k ) r ” ( 2 2 2 ) 使得函数 c ( x ) = i = lt 仍( x ) = k r 矽( x ) ( 2 2 3 ) 矽( x ) = ( 仍( x ) ，纯( x ) ) 1 满足式( 2 1 2 ) 、( 2 1 3 ) 令c ( x ) 对应于式( 2 1 2 ) 边值问题的解为“( x ，；c ( x ) ) ，( x ) = 砰仍( x ) = 耐痧( x ) 为 真解c + o ) 附近的一个函数，是优化搜索的初始点，对c o ( x ) 增加一个微小的扰动量 ( x ) = 万矽仍( x ) = 万瑶( x ) ( 2 2 4 ) 对应与c 0 ( 功+ 5 c o ( x ) ，式( 2 1 2 ) 的初边值问题的解记为u ( x ，t ；c o ( x ) + 。( x ) ) 所以，c ( x ) 的确定问题可以转化为万民的确定问题并_ 且s k o 可由下列目标函数的局部极小值来确 定： f 6 k o = ( 1 一五) i l “( x ，t ；c o ( x ) + s c o ( x ) ) 一缈( x ，) i i ( 耐。l o 7 t j ) + 名s ( s k o ) ( 2 2 5 ) 其中，s ( 6 k o ) 是万的稳定化泛函 5 c o ( x ) 是一个微小的扰动量，把u ( x ，t ；c o ( x ) + 8 c o ( x ) ) 用多元泰勒公式展开得到： 于是有： u ( x ，f ；( x ) + 万( x ) ) = 甜( x ，f ；o ) ) + v 乙“( x ，f ；( x ) ) 8 i c o + o ( 1 18 c o ( x ) 1 1 ) ，陋c o = ( 1 一名) 1 1 甜( x ，f ；c o ( x ) ) 一q g ( x , t ) + v t 甜( 而，；( 石) ) 万甄峨删咿。( 2 2 6 ) + 兄s ( 万k ) 若在区间硷7 【o ，丁】上有m 个离散点( ，t m ) ( m = l ，2 ，m ) ，而取s ( 万) = 万碍万r o 则 研讽】= ( 1 一句【z f ( ，乙；( 瑚一贴。，乙) + 吒z ，( ，乞；( 砌6 k o 2 + 2 5 k r g k o ( 2 2 7 ) 令 1 5 山东理t 人学硕l j 学位论文第二审模型数值解法及j e 反演算法的改进 u = 甜( ( x ) ；五，f 1 ) u ( c o ( x ) ；x 2 ，t 2 ) u ( c o ( x ) ；x m ，0 ) v = 缈( 而，) 矽( 叠，2 ) 缈( ，0 ) a 2 ( f ) 掰。 = 熹氓，t m ；c o ( 砌u 二 则上式简化为： f f s k o 】= ( 1 一名) 万磁彳r 彳万k + 2 ( 1 一彳) 万碍彳r ( u 一矿) + ( 1 一允) ( u y ) 7 ( u 一矿) + 硒k ：8 k o 不难求得f s k o 】的局部极小值万k 满足下列线性方程组： o - 允 ) 么r 彳+ 五 i 6 k o = ( 1 一名n ) 么r ( 矿一u ) 其中，名( ”) 为第刀步迭代的同伦参数，么为导数矩阵 2 4 算法迭代步骤及流程图 ( 2 2 8 ) ( 2 2 9 ) 改进最佳摄动量算法的迭代步骤： ( 1 ) 给定初始向量k o ，参数夕，n o 以及迭代终止误差e p s ，迭代步歹= o ； ( 2 ) 求解正问题( 2 1 2 ) 的数值解，得到甜( 靠，乙；巳( x ”； ( 3 ) 利用数值微分计算 导甜( 毛，乞；q ( x ”= 三陋( ，t i n , o n ( x ) + 磁( x ) ) - - u ( 靠，o ；巳( x ) ) 】，求解迭代时的导数矩阵 o k i i a ； ( 4 ) 求解( 2 2 9 ) 得到参数万k ，从而得到阮( x ) ； ( 5 ) + 。( x ) = 厶( x ) + 氓( x ) ，并刀= n + l ，返回( 2 ) 重复上述过程，直到满足精度要 求为止 用m a t l a b 编写计算程序，程序框图如下： 1 6 山东理t 大学硕卜学化论文 第_ 章模型数值解法及j e 反演算法的改进 2 5 算例 图2 5 改进摄动量算法流程图 算例3 考察某溶质在给定土层中的运移问题 ( 1 + p b k a ) 0 优c = 吼v 嘉一v 塞，o 刚，0 m c ( x ，0 ) = 0 ( 2 3 0 ) c ( o ，f ) = c o c a t ，) = 0 1 7 山东理t 大学硕i ：学位论丈 第_ 币模，弘数伉堵i 法及e 反演算泫的改进 其中c ：c f x , t ) 为所考察溶质的浓度，已知的参数是土层的十容重岛= 1 4 5 ，体积含水率 秒= 0 2 ，运移距离，= 6 2 ，平均流速1 ，= 2 4 x l o - 3 ，而吸附系数钇与纵向弥散度日：是未 知的参数这时，需要根据一定的观测数据，来反演求解髟，口r 这里用到的附加数据是 x = ，处的出流数据 先进行模型的无量纲化另c = 三，z = ，t = ；，可得 (1+pb6，kd)ac=一al万ac一瓦ac，oz1，o丁 m g ( o h ) 2 其中，乞表示土层中硫酸钙的分解系数；后g 、k 分别表示石膏( 含2 个水分子的结 晶硫酸钙) 与氢氧化镁的沉淀系数；后。表示土层中含若干水分子的结晶硫酸镁的裂解系 数 记液相中钙离子、硫酸根、镁离子以及氯根的浓度分别为q 、乞、c 3 与c 4 ，并用墨 与矗分别表示土层中固有的硫酸钙与结晶硫酸镁的浓度 联合水动力弥散方程与水化学反应式( 2 3 0 ) ，可得： 1 9 山东理t 大学硕l ：学位论文 第二章模型数值解法及其反演算法的改进 鲁= 。鲁一v 鲁蝴电，乞 鲁= 噜一v 鲁瑚如m 文 亿3 3 ， 鲁= 。等一v 芸+ 哎是一吒乞 0 c 4 ：d 1 a z c 4 一，堕 a t苏2叙 同时，注意到岛= s t o e 一和s 2 = s 2 0 e 一婶，其中而o = 毛( o ) ，龟o = 岛( 0 ) ，再令= 包， r z = 包而。，r 3 = 镌，r 4 = k s ，r s = 屯，r 6 = k ，则上述方程组( 2 3 3 ) 改写为： 鲁= 。鲁一v 鲁岬叫唧。巳 鲁= 。鲁a 以c 2 + r 2 e - ；_ r 3 q c 2 + r 5 矿 鲁= 。等一v 鲁坶叫啪 q 门4 a c 4 ：d _ a z c 4 一v 亟 a t耐瓠 f o x l ：0 t z ，、 相应地有初始条件： e ( x ，0 ) = 0 ( 2 3 5 ) 以及边界条件： q ( o ，于) = q o a( 2 3 6 ) e ( j ，) = 0 t t x 其中，f = 1 ，2 ，3 ，4 方程组( 2 3 4 ) 联合初边值条件( 2 3 5 ) 一( 2 3 6 ) 即为描述矿坑水中四种主要离子 c a 2 + ( q ) ，孵一( 乞) ，m 9 2 + ( g ) 以及c 厂( 乞) 在稳定流速下穿越匀质土层时浓度时空 变化的数学模型该模型实际上是一个带有源( 汇) 项的对流弥散方程组，其中的源( 汇) 项恰恰刻化了多组分溶质运移过程发生的物理化学作用模型中的六个参数，；、眨、吩、 ，= i 、吩与r 6 是难以直接观测得到的，这些参数通过数值反演得到 为便于数值计算，对上述方程( 2 3 4 ) 及初边值条件( 2 3 5 ) 一( 2 3 6 ) 进行无量纲化处 理令e = q q o ( i = l ，2 ，3 ，4 ) 以及z = x l ，丁= v t l l ，p = v l ld ，并记弓= r s l v ( j = l ，2 ，6 ) ，如此有： 山东理t 人学硕l j 学位论丈第一二章模犁数值角i f 法及j e 反演算法的改进 万a c , = 万1 万0 2 c i 一瓦a c , + 1 q 。r 2 e _ r i r _ c 2 0 b c l c 2 璺三：丢要鲁一坠+上恐p一焉rci。恐clc2+上bp一即at尸a z 2a z c 2 。2 一2 c 2 0 ( 2 3 7 ) 坠：三堡宴一坠+上马p一即一rc30tpa za z z g n o， a c 4 1a 2 c 4 a c 4 - - - - - = 一- - - - - - i 一- - - - - - a tpa z 2a z 相应的，初边值条件变为 e ( z ，o ) = 0 ( 2 3 8 ) 及 g ( o ，丁) = l 昙c ：f ( 1 = 。 q 。 如前所述，上述模型中所包含的六个参数r ，( ，= 1 ，2 ，6 ) 是反映矿坑水淋滤渗入 土层时可能发生的物理化学作用的模型参数，它们是难以通过直接测量获得的根据土 柱淋滤试验的实际情况，一般是在土柱出流端定时接收样本溶液，通过对其分析化验而 获得各组分溶质的浓度穿透曲线数据换言之，为了确定模型中的未知参数，我们给定 如下一组附加数据( 无量纲化) ： e ( 1 ，瓦) = c ：( 瓦) ，k = 1 ，2 ，k ，f = 1 ，2 ，3 ，4 ( 2 4 0 ) 其中k 表示试验所得的样本个数 这样，由方程组( 2 3 7 ) 与初边值条件( 2 3 8 ) 一( 2 3 9 ) ，联合附加数据( 2 4 0 ) 就构成 一个确定模型参数尺，( = 1 ，2 ，6 ) 的参数识别反问题 对模型进行数值模拟本次土柱淋滤试验所用矿坑水含硫酸根以及钙离子、镁离 子等盐类溶质，矿坑水中钙离子( q ) 、硫酸根( 乞) 、镁离子( 巳) 以及氯根( 乞) 的含 量分别为q o = 3 3 8 2 8 ( r a g l ) ，c 2 0 = 1 0 6 2 9 ( m g l ) ，c 3 0 = 1 0 4 4 2 ( m g l ) 以及 c 4 。= 2 1 9 1 ( m g l ) ；所用的淋滤土柱长，= 0 4 5 ( m ) ，土柱中稳定达西流速为 v = 5 1 x l o _ 6 ( 朋i s ) ，弥散系数d = 1 0 _ 9 ( 朋2i s ) ，矿坑水淋滤总时间为= 1 2 1 ( h ) 反 演时所用的附加数据由假定的参数真值代入模型中计算而产生 根据上述参数值、模型( 2 3 7 ) 一( 2 3 9 ) ，将无量纲的时间t ( o ，v t o , 1 ) 与空间 z ( o ，1 ) 作2 0 等分，取参数真值为屉= ( 1 , 0 5 ，0 0 1 ，1 ，2 ，1 5 ) ，分别应用一般最佳摄动量算 法与改进的最佳摄动量算法进行计算数值实验结果比较如表2 3 其中，方法( 1 ) 为原最 佳摄动量法，方法( 2 ) 为修正的最佳摄动量法初始值为( o ，0 ，0 ，0 ，0 ，o ) ，应用改进最佳摄 动量算法计算，计算结果见图2 7 2 1 山东理t 大学硕i j 学位论文第_ 幸模型数值解法及j e 反演算法的改进 p o o o 0 c o 3 2 1 o 2460 t ( a ) 寸 c ) 246 t 彻 02460 t 【c ) 246 t 曲 图2 7 ：初始值为( 0 ，0 ，0 ，0 。0 。o ) ，方法( 2 ) 的计算结果 表2 3 初始值相同时两种算法的比较 山东理t 大学硕i 学位论文第二章模型数值解法及j c 反演算法的改进 2 6 小结 最佳摄动量算法属于迭代正则化方法，是一种局部收敛、依赖初值