UCODE反演程序的原理及应用

第17卷第6期2010年11月

地学前缘(中国地质大学(北京);北京大学)

Earth S cien ce Frontiers (Ch ina University of Geosciences(Beijing);Peking University)Vol.17No.6Nov.2010

收稿日期:2010-07-01;修回日期:2010-08-01基金项目:国家自然科学基金重点项目(50639090)

作者简介:夏 强(1982 ),男,博士研究生,地下水科学与工程专业,主要从事地下水数值模拟方面的研究。E -mail:qian gwa@

U CODE 反演程序的原理及应用

夏 强

1,2

, 万 力1, 王旭升1, E Poeter

2

1 中国地质大学(北京)水资源与环境学院,北京100083

2 科罗拉多矿业学院国际地下水模型中心,美国戈尔登80401

Xia Qiang 1,2, Wan Li 1, Wang Xusheng 1, E Po eter 2

1 S ch ool of W ater Resource s and E nv ir onme nt,Ch ina Unive rsity of Ge oscience s(B eij ing ),Beij ing 100083,China

2 International Gr ound W ate r M od eling Center ,Color ad o S chool of M ine s,Gold en,Colorad o 80401,US A

Xia Qiang,Wan Li,Wang Xusheng,et al.Principles and applications of the inverse problem program:UC ODE.Earth Science Frontiers ,2010,17(6):147-151

Abstract:T his paper illustrates the sig nificance o f calibration fo r g r oundwater mo deling,and demo nstr ates that aut omated calibration techniques using inver se problem prog r am ar e super ior to manual tria-l and -er ro r meth -ods.T he w idely used U CODE is o ne of such pr og rams w hich o pt imizes the parameter v alues by Gauss -New ton methods.T he initial par ameter values play an impor tant ro le to ca librat ion.A synthetic transient model is co n -str ucted,and six numerical ex periments ar e perfo rmed to v erify t he practicability of U CODE prog ram.T he r e -sult s show that althoug h the initial v alues of parameter s w ould influence the pr ocedur e of calibr atio n,g iv en ap -pr opriat e v alues,U COD E co uld achieve the objective for optimizatio n.Key words:U CO DE;inv erse problem;g roundw ater;modeling

摘 要:对地下水模型进行反演是模拟过程中的一个必要步骤,使用反演程序自动校正模型可快速确定最佳拟合的参数值,分析参数对模拟结果的敏感性,比人工试算-调整法更为优越。U CO DE 是一款被广泛应用的地下水模型反演程序,它使用高斯牛顿法进行参数优化,反演结果对参数初值有一定的依赖性。通过建立假想的非稳定流模型,进行6组数值试验,验证了U CO DE 程序的实用性。尽管参数的初始取值会影响反演的进程,但只要取值适当,U CODE 就能实现优化参数的目的。关键词:U CO DE;反演;地下水;模拟

中图分类号:P 641 2 文献标志码:A 文章编号:1005-2321(2010)060147-05

在地下水系统的研究分析中,数值模拟正在被越来越广泛地使用。大多数的地下水模型都是分布式参数模型。从本质上讲,这些具有一定物理意义的数值模型其实是使用有限差或有限元的方法近似求解的[1]。建立数值模型之后,通常需要校正模型。校正模型的过程就是要调整输入模型的参数,直到模型输出的结果与野外观测数据达到一

定程度的拟合。

在实际应用中,从来不能很完善地定义模型输入参数。无论对参数进行了多少次测量,也无论对野外条件刻画得多么详尽,它们总是具有一定的不确定性。因为,模型中某些参数虽然是通过实际测量获取,但实测和模型之间尺度上的差异,还是可能导致模拟结果与观测数据较大的偏离。所以,我们

很难利用初始参数就使模拟结果和野外观测高度一致[2-3]。校正是优化参数值的重要方法,从这个意义上讲,模型校正(mo del calibration)也相当于参数估计(param eter estimation)和参数识别(parameter identification)。通常认为参数估计与识别主要关注参数的取值,而地下水模型的属性大多具有空间变异性,对于此类模型,校正其参数的取值和结构(如渗透系数如何分区)都是很重要的。

相对于通过模型得到系统状态的正向问题而言,校正模型属于使用观测数据而求解系统参数的反问题(inverse problem),所以模型校正也称之为 反演 。对模型进行反演可以采用人工试算-调整的方法(tria-l and-er ror),也可以应用反演程序自动校正(autom ated calibration)。后者具有明显的优势,如快速确定最佳拟合的参数值,量化校正结果,以及量化数据缺陷等[4]。

基于非线性最小二乘理论的地下水反演程序有MODINV(Doherty,1990)、M ODFLOWP(H ill, 1992)、PEST(Doherty,1994)和U CODE(Po eter, 1998)等,此类程序的开发、应用使模型的自动校正分析更加方便、快捷[4]。U CODE是一款常用的地下水模型参数反演程序,在商业化的地下水数值模拟软件Pr ocessing M ODFLOW、V isual M ODF-LOW和GM S等中都有镶套。该程序由美国科罗拉多矿业学院国际地下水模型中心的Poeter教授于1998年开发,到目前已发展到UCODE_2005[5]。

UCODE正被越来越多地应用于流域模拟的实践之中,如Geza等(2009)使用U CODE程序对美国科罗拉多州T urkey Creek流域的WARMF(流域风险管理分析)模型进行敏感度分析,优化了20个参数中敏感度较大的7个参数,观测与模拟流量拟合情况较好。Fog lia等(2009)使用U CODE校正了瑞士Magg ia流域的T OPKAPI模型,利用基于误差的观测数据与先验信息的权重分配、局部敏感度分析以及单个目标函数非线性回归等方法,得到模型35个参数敏感度的量化评价,识别出对模型校正最重要的数据类型,以及影响数值解的非唯一性的参数的相关性。

1 U CO DE的优化原理

模型参数反演方法分为两类,直接解法和间接解法[6]。直接法需要充足的数据以确定模型的水位分布或溶质浓度分布等,把它们作为已知变量,而把模型参数作为因变量,直接求解。间接法则不需要直接法那么多的数据,它基于最优化目标函数的原理,不断执行正向问题,使用最优化方法改进参数取值,直到相邻迭代之间参数取值无明显改变[2]。

Yeh(1986)综述了十余种直接参数反演方法和二十余种间接参数反演方法。用于识别水文地质参数的最优化方法主要有最速下降法、逐个修正法、高斯牛顿法、Pow ell方法、单纯形法、线性规划法、二次规划法、拟线性化方法和罚函数法等[7]。

国内学者在这方面也进行了深入的研究。姚磊华等(2003)提出了一种改进的遗传算法进行参数识别,并进一步充分利用遗传算法善于进行全局搜索和高斯牛顿法善于进行局部搜索的优点,于2005年用改进的遗传算法和高斯牛顿法联合反演地下水数值模型参数,克服了两种方法各自的不足。

1 1 高斯牛顿法

UCODE中使用高斯牛顿法进行参数优化。相对于其他的方法,高斯牛顿法最大的优点在于求解局部最优解时,只要参数初值选择适当,收敛速度非常快,收敛精度非常高[7]。

地下水模型可简记为下面的非线性表达式[8]:

Y=f( ,b)+e(1)其中:f为n维列向量,这里可认为是模拟得到的水头值; 为自变量矩阵,表示模型中的已知量;b是参数向量,为p维列向量(必须满足n p),b也是反演问题中的求解对象;e是误差向量;Y为水头观测值,它包含了模拟结果f和误差e两部分。

而高斯牛顿法就是要使误差平方和最小,由此确定最优的参数向量^b,令:

S(b)=e T e=(Y-f( ,b))T(Y-f( ,b))(2) f是非线性的,首先使用泰勒公式在初始参数b0处,将f改写为线性形式:

f( ,b) f( ,b0)+X0(b-b0)(3)其中:X0={x0ij}=

f i

b j b=b0(n p)

为敏感度矩阵,f i是第i个观测点的模拟值表达式, X0的各个元素被称为敏感度系数。

将式(3)代入式(2)得到:

S(b) [Y-f( ,b0)-X0(b-b0)]T

[Y-f( ,b0)-X0(b-b0)]

此时,求S(b)对单个参数b j的导数,并令导数为0,经过简化得到: