条件模拟原理和技术进展及在水资源系统中的应用

1、条件模拟原理和技术进展及在水资源系统中的应用陈亚新,徐 英,史海滨,魏占民(内蒙古农业大学水资源与水土工程研究所,内蒙古呼和浩特 010018摘要:在简要介绍条件模拟基本原理基础上,重点评述了国外20世纪90年代新的数学方法和计算技术进展及在国内外水资源系统应用的主要成果,并对今后学科和技术发展作了展望。关 键 词:地质统计学;条件模拟;水资源系统;应用;进展中图分类号:P 628+ 2 文献标识码:A 文章编号:1001 6791(200202 0246 07条件模拟(Conditional Simulation 简写CS是新兴边缘学科地质统计学(Geostatistics 简写GS的一个具

2、有近代科学内容的重要组成部分。自G Matheron(1970倡导的GS 问世以来,对地球科学中有关资源信息空间变异性用区域化变量理论进行表征,不仅对地质预测和评估理论有很大发展,也对许多非地质科学如水土资源、水文气象、生物环境和工程技术1产生广泛影响,目前已经形成研究空间信息分布特征或模拟其离散性和波动性,均可用GS 及相应CS 理论进行研究,成为评估各种区域性自然现象、自然资源及再现其波动过程新的工程科学。20世纪70年代中期数值计算方法和新的概率统计方法,在计算机的支持下获得广泛应用,改变了工程技术领域依靠物理模型和数学解析的传统方法。许多具有随机性的自然现象及其过程,常用统计数值方法模

3、拟试验,统计模拟试验方法虽可使其结果逐渐接近真值,但总要先用计算机产生具有一定统计特性的伪随机数,服从一定概率分布具有给定数学期望和方差,再用数学运算方法解决各种数学统计模型,即一般所称Monte Carlo 试验模拟(简写M-C2。在GS 中,因研究对象是一种具有结构特征的随机变量,对它进行数学模拟时除要用到M-C 方法外,模拟时要求更高,还要求保持一定的空间相关性,即保持实验信息具有相同的协方差函数C(h或变差函数r (h ,在GS 中称非条件模拟;如果再增加一个更严格的约束条件,将模拟条件化,即令在各观测点处的模拟值等于该点的实测值则称CS 。CS 比传统M-C 模拟有一些新特点:它不但

4、能保持变量的空间相关性如C (h 或r (h 不变,能使观测点处的模拟值等于实测值,当控制的观测点越多时模拟就越接近客观实际,并易于实现从一维模拟扩展到三维空间模拟;在分析品位空间变异过程,若用CS 重现真实试验值的离散方差 2或变差系数C V ,能再现真实资源特性空间变异曲线的波动。自A G Journel (1974、1978发表GS 的著名论著3,4,完善并系统总结了CS 的理论方法和应用,常见用于新矿区的开采生产设计,可在开发和大规模钻孔开钻前对资源作出开采预测;第二、三届国际GS收稿日期:2001 01 02;修订日期:2001 05 30基金项目:国家自然科学基金项目(501790

5、13作者简介:陈亚新(1936-,男,湖南长沙人,内蒙古农业大学教授,博士生导师。主要从事溶质运动及其空间变异性、节水灌溉等研究。第13卷第2期2002年3月 水科学进展ADVANCES I N WATER SCIENCE Vol 13,No 2 Mar ,2002大会期间(1983-19885,6,为了适应矿山与选矿厂设计阶段的配矿要求Journel 和Isaacs 开发了适合混合型矿层协同区域化变量联合估计CS;Guarassim 和Musso 等并将C S 作为矿床配套设计中的重要组成部分完善了选矿的生产管理;H M Parker 还用CS 指导特定环境中矿带的各种不同勘探、设计的储量估

6、计方法;20世纪80年代在国外水资源系统获得广泛应用。第28届国际地质大会(1989上还提出了非线性GS 的CS 研究成果,使CS 理论和方法开始从线性向非线性领域推进,CS 在不同条件获得推广。经多年实践以后M W Davis(1987评价:C S 是一种十分有效的方法,是对研究GS 的最新理论和实践应用的发展,还可以进行多元条件模拟。我国从20世纪90年代在地质与非地质(水资源、水-土环境710领域开始应用研究,受到广泛关注。1 条件模拟的基本原理和应用方法1 1 基本原理设Z ( 是满足二阶平稳假设的区域化变量,E Z( =m 并存在C (h或r(h。要求Z( 的条件模拟为Z cs (

7、,即要求得与Z( 具有 同构 的区域化变量Z cs ( 的一个现实,且在实测点X 上的模拟值等于实测值,即Z cs ( =Z( , , =1,2,!,n (1为了求得Z cs ( 的理论方程,在此引入OK 估计值Z *k ( 和非条件模拟值Z s ( (即其他条件均满足,只不满足条件(1式的模拟值,由于Z ( 在任一点 处的真实值Z ( 可表示为OK 估计值Z *k ( 与其误差R ( 之和,即Z( =Z *k ( +Z( -Z *k ( =Z *k ( +R( (2式中 R( 为Z( OK 估计的未知误差。由于OK 估计误差Z ( -Z *k ( 只取决于信息的数据构形并正交(独立于OK 估

8、值,在任意两点x 和y 上则有E Z *k (y Z ( -Z *k ( =0; x ,y (见文献4。利用这一特性,用一个与此误差 同构 且独立的非条件模拟OK 误差Z s ( -Z *sk ( 来代替(3式中的未知OK 误差Z (x -Z *k ( 即可得条件模拟理论方程:Z cs ( =Z *k ( +Z s ( -Z *sk ( (3式中 Z *sk ( 是非条件模拟的OK 估计值。因为:从OK 估计量的无偏性,可知随机函数Z sc( 和Z ( 具有相同的期望值;被模拟的OK 估计差值Z s ( -Z *sk ( 同构于真实误差Z( -Z *k ( ,且独立于Z *k ( ,随机函数Z

9、 sc ( 同构于Z ( ,有相同的变差函数;被模拟的现实Z sc ( 是以已知点的实验值为条件,即Z sc (x =Z ( ,故(4式可满足条件模拟的要求。因此可用Z s ( -Z *sk ( 代替Z( -Z *k ( 。由于(4式中Z s ( -Z *sk ( 是可求出的,若用它代替(2式中未知的真实OK 误差Z ( -Z *k ( ,(2式右端值即可求得,但求得的不是Z( ,而是Z ( 的条件模拟值Z cs ( ,可证明它满足条件模拟 同构 的要求4。由于Z s ( 与Z( 有相同的变差函数,且求OK 估计值Z *sk 和Z *k ( 时的数据构形又相同,故普通Kriging (简写OK

10、 方程组也一样,方程组的解也一样,即有相同的权系数 , =1,2,!,n 。于是将Z *k =!n =1 Z ( ,Z *Sk ( =!n=1 Z S ( 代入(4式得条件模拟实用方程:Z cs ( =Z s ( +!n =1 Z( -Z s ( (5(5式表明要计算CS 值Z cs ( ,只要先求出一个非条件模拟值Z s ( ,再对实测点 上247 第2期陈亚新等:条件模拟原理和技术进展及在水资源系统中的应用248水科学进展第13卷的差值Z( -Z s( , =1,2,!,n进行OK估计,最后相加即得Z cs( ,此方程可减少一次Kriging方程组的运算,简单实用称CS基本方程,若对(5式

11、反复进行模拟就得到任意多个CS的现实,成为一种新的M C方法。1 2 原理应用的基本方法为应用CS实用方程,要先进行一维空间的非条件模拟求Z s( ,如采用滑动平均法实现离散逼近;再进行三维空间的非条件模拟求Z cs( ,将一维空间非条件模拟扩展到三维空间,往往要耗费较多机时,常采用G Matheron(1981首创的转向带(Turning band简写TB法4,将所有三维模拟约简为沿若干(N直线的几个独立的一维模拟在三维空间(R3内进行旋转,进行带权叠加得到多维随机场的模拟现实,取来自N条直线N个贡献之和最后得空间任一点 的一个最后值Z s( :最终完成(5式的三维条件模拟Z cs( 。条件

12、模拟计算流程框图见文献10。2 条件模拟技术进展20世纪80年代后期至90年代,C S基本理论在应用数学方法和计算技术等方面获得了很快的发展,使其适用性更加广泛,可针对各种复杂的生产需求,计算技术也更加简捷1527。(1为了克服TB法对垂向各向异性差别大的问题和要求有较多网格采样点及快速Fourier (简写Fo变换对模拟尺度大小的局限,M W Davis(1987提出LU三角分解协方差矩阵快速、优质、简便的CS方法,完成Dallas的一个污染区Pb浓度模拟三维网状图和等值线图,与传统TB法(1981和Fo变换法(1984比较有明显优点,如模拟和条件化可同时进行并通过运用协方差矩阵可进行大量运

13、算,数据点没有和模拟点重合。对于中等尺度(750个网格节点条件下CS提供一种简单算法,尤其适用于大数量小网格分析研究。(2GS的应用一般要求符合二阶平稳假设其CS也多限于平稳问题。而对于非平稳领域含有漂移问题的求解由于具有循环证的矛盾,虽产生了K阶本征函数方法(IRE#K企图绕过漂移从非平稳的Z i( 直接过渡到平稳的Z( ,但少见用于C S。Dimierakopoulos R(1990提出规则K阶本征函数的C S,形成模拟非平稳现象的一种随机方法,它可以产生广义协方差和将取样位置得到的有效数据,直接求得贴近模拟结果并深入考虑了各向异性问题。(3传统的CS一般均是在已知信息(控制点间进行内插模

14、拟,Hoshiya(1995提出在Gauss 离散域的Kriging估计和逐步向信息域以外的有效外推C S方法,避免协方差矩阵的C holesky和Modal分解,是基于最优估值和方差正交基础上具潜在用途。表明当有多组观测值时,其观测值更新的逐步外推扩展过程与Karlman滤波算法相当。(4过去人们很少注意到原始数据含有观测误差的C S方法,D Marcotte(1995提供了观测数据具有噪音的过滤处理与CS方法,用模拟前过滤和模拟后过滤两种方法进行检验比较,结果表明先模拟后过滤效果较好,过滤方法是用多元地质统计学中的因子Kriging(Fk法和修正因子Kriging(MFK法,可确保基本变异

15、函数很好再现。(5为了解决简单Fo条件模拟的Sinc卷函数在小尺度网格上不能精确与期望的统计模型相匹配,J A Goff和J W Jemings(1999提出极限带部分统计模拟对CS的Fo变换的改进方法。用协方差法创造人工高波数含量改进Fo变换,可弥补Sinc平滑效应,评价和改进了非条件和条件模拟Fo变换的统计精度及其性能,这种快速Fo法尽管其精度不如OK法,克服了条件化数据点数量大、需花费机时过长缺限且不产生明显的数据失真,对Fo 条件模拟的后验统计也不产生坏影响。可适用于有大量数据条件具备各向异性模型的有1001个点以上大区域。这些新方法和技术开拓了CS 基本方法在更加广泛和复杂领域中的应

16、用,也提供了高效、精确、适用的新技术。为推动CS 在各个资源领域的应用创造了更严格和适用的数学基础。3 条件模拟在水资源系统中的应用(1国外水资源系统中的应用 多见于地下水和水文地质两方面2527,已经进入用多种方法综合比较研究阶段。D M Cliffon 和S P Neuman (1982最早在Arizona 南部的Arra 含水层用Kriging 估计效果和逆问题的CS,研究和比较了三种不确定水平条件下的随机模型,其最高水平含水层水头是通过多元正态随机函数的非条件模拟所测定,协方差用CS 方法确定;中等水平条件下的不确定性通过非条件模拟测定,认为这种情况下的CS 是一个全新的概念,虽然CS

17、 过程的效果比Kriging 估计好,但它的水头预测值是Kriging 方差( 2k 1 0倍,而逆问题模型方差比 2k 可缩小14 3倍,1992年UNESCO 召开了GS 在地面水#地下水水文方面最新进展和应用国际讨论会,D Weber 对前人:(1Cliffon 和Newman (1982,(2Hoeksena 和Kitanidis (1985,(3Strecke 和Lettemairer (1987,(4Wagner 和Gorelick (1989,(5Weber 和Easley (1991等五项即近10年来利用CS 和参数估数、逆问题模拟方法和Kriging,Co Kriging 等

18、方法获得的各种方案在各向同性假设下作了系统比较见表1。表1 五个代表性研究成果的比较Table 1Comparisons of results o f five representative studies方 案(作者代号测量资料估计方法估计参数生成技术方案或结果改良系数(MSEH(1(1(2(2(3(3(4(5T H T ,T ,H T ,C H ,C H K ,H Kriging 逆问题模拟Kriging CO Kriging 表面逆问题模拟逆问题模拟Caus sion 条件中值T T log T log T T T K log K 条件模拟条件模拟数值模拟数值模拟数值模拟数值模拟数值模拟

19、条件模拟或数值模拟1 014 31 015 01 06 2综合认为:导水系数(T 的随机组成能经CS 引进再现水头(H 域的空间变异性,但要用CS 精确确定估计网格位置处H 值的协方差和均值仍有一定困难;从方案(1中H 的逆问题模拟与T 的OK 法相比较其改良系数为14 3倍(协方差;方案(2中的T 和H 的C O -K 估计比T 的O K 法相比改良系数为15 0倍,方案(3H 域逆模拟与测量T 域估计相比改良系数为62;由H 水头域值可得到最优拟合值(代表精确值;用C S 技术生成的许多反映H 和导水率(K ,可用于地下水质矿化度(C 管理模型。虽然在5个方案间还不能进行横向直接比较确定那

20、种方案最有效,但对利用C S 和参数估计方法对地下水观测井设置方案具有提高监测效果和治理质量并产生一定开发潜力的值价给于了肯定。Journel(1994已开展了水文地质平均地段渗透性的CS,使CS 研究从常见的观测网格点的预测发展到区域性块(段平均的CS,进入了更大范围的应用。(2我国水资源系统中的应用 20世纪90年代中期以来,主要集中在包气带的土壤水、249 第2期陈亚新等:条件模拟原理和技术进展及在水资源系统中的应用土壤资源与环境科学的监测与评估研究方面810。张泽浦、王学军(19988将CS用于北京惠通河污水灌溉小尺度(46%10m%(75%10m=34 5hm2土壤重金属Cr,Cu,

21、Pb!等8种微量元素的空间分布,研究了区域化变量存在明显漂移的CS,对残差进行模拟时由于模拟函数呈近似对数动态分布、经Gauss变换计算协方差,再对条件模拟结果进行逆Gauss变换,结果表明CS可以很好地再现土壤的微量元素的空间变异结构特征为水-土环境质量评估提供了有力工具。徐英、陈亚新等9(1999在对黄河灌区冻土大量定位观测基础上,发现季节性冻土层的形成可能是土壤水盐信息空间变异的地学结构背景,特别是最大冻深期,信息变异性比较稳定,可用CS 和OK估值,认为如果水盐预测或监测目标不需求得垂向某点的未知信息,只需预测其空间分布的总趋势时,可用较少的采样点时的已知信息,用CS揭示全剖面的未知信

22、息波动规律,再现不同时期的水-盐时空变异的波动,可以减少冻土外业采样,有实用价值。此外还对C S过程的M-C试验模拟技术进行了较深入的探讨,由于CS用随机数系列求解,是一种新的M-C方法,为了实现一维协方差函数离散和C(1(S向积分式基台值C(1(0的渐近趋势,需要高密度的离散逼近,虽Journel(1978曾建议离散逼近时随机数系列长度R(自然数的选取一般R&20,但并未指明其上限值,也未涉及随机种子数R d的取值及信息特征对模拟效果的影响,应用有一定困难。研究认为,土壤水分在不同定位点和不同冻融时期有相似规律,即20R100时均值和C v都模拟比较稳定并接近实验值;当模拟信息目标不

23、同,R值的变化随水分和盐分取值而应有不同,当盐分R为定值而R d变化时,随R d变化均值和C v跳动都比较大,但当R&100后R d也渐趋平稳,可作为盐分R d取值的首选范围。上述成果可供离散逼近的M-C过程应用的补充参考。陈亚新、史海滨等10(2000在黄河内蒙古灌区(1993-1999实验基础上,探讨了大面积农田防治土壤盐渍化的水盐监测,在中小尺度采样810hm2(50m%50m,24hm2(10m %10m两个实验区的控制点间(内插进行了CS与OK估值的全面比较,如表2。表2 条件模拟和Kriging估计与实验值比较Table2Comparison of the conditio

24、nal simulations with estimative and observed values项 目 秋(1993 春(1994 tes t CS O K test CS O K水分#mean/%20 9420 8320 9028 6228 4428 7C V0 130 120 100 070 060 05C V/C V tes t10 910 7810 910 67盐溶质TDSmean/%0 220 210 220 210 230 20C V0 320 290 250 460 450 42C V/C V tes t10 930 7910 980 91成果表明CS的模拟方差接近实验值及

25、其空间变异的波动过程,在掌握一定距离若干控制点的采样值以后,可以较少钻孔模拟孔间并再现空间变异规律及其离散性和波动性有突出优点;但模拟值不是最优估计值,不如OK估计值能更接近平均实验值,因而模拟和估计是两种不同概念和从不同角度揭示数据的空间特性,各有不同优缺点,不能相互代替,应根据不同目标和需求选择不同方法或结合使用;由于C S可再现真实值的方差 2(或C V还可通过模拟的 2(或C V,能较好地再现水盐特性和待估的田间监测系统的宏观模拟合理采样数N CS,可由疏到密求得不同设计尺度的估计采样精度和评估经济效益。250水科学进展第13卷第 2期 陈亚新等: 条件模拟原理和技术进展及在水资源系统

26、中的应用 251 4 展 望 ( 1 在水资源系统中引进 GS 时要考虑其时空变异特点, 既要符合 GS 的基本假定又要具 有水资源变量的时空变异特点, 不同于地质变量多是在宏观时间大尺度下形成背景, 一般只研 究空间分布特性; 而在水资源中的水文变量还具有中观和微观尺度的时间变化特征, 研究目标 是时- 空域的随机变量, 要将 GS 的 CS 引进 到水资源领域 应用, 需要 开发多元地质统 计学 ( MVG , 利用时- 空协同区域化理论进行 Co Kriging 估计和因子 Kriging 法( FKA 的 CS。 ( 2 目前 CS 研究与应用还是以中小尺度为主, 大尺度的预测虽是热点

27、, 但又是难点, 单 纯依靠采样钻孔进行预测或检验有一定难度, 将 CS 的数据采集技术与 3S 技术( GIS、 GPS、 RS 结合, 有可使大尺度时- 空的 CS 获得突破, 有待进一步实行多学科的合作攻关。 ( 3 GS 及其 CS 理论正在带动地球科学与相关学科的相互促进, 向现代化多学科的综合方 向发展, 正在从单变量向多变量, 连续型变量向离散型变量、平稳型向非平稳型、线性向非线 性方向发展。从内插模拟方法补充外推模拟, 国内外大量成果表明, CS 有能力直接为国民经 济建设和管理服务, 可逐步解决 21 世纪的资源和环境问题。 ( 4 尽管我国在地质与水资源领域应用 CS 尚处

28、于起步阶段, 但从学科的特色优势和国内 外大量成功事例可预计综合运用 GS 的 CS 理论、模型、方法与近代计算机技术结合, 围绕与 社会长远发展目标密切相关的学科前沿, 协同解决水资源系统的智能化, 定量化和模拟、仿真 前景十分广阔。 参考文献: 1 陈亚新, 史海滨 地质统计学在水资源系统的 应用和发展 J 内蒙古水利, 1997, ( 1 : 12- 16 2 Shi Haibin, Chen Yaxin A rational sampling number and Monte Carlo simulation test in information management of soil

29、 water resources in irrigated area A Press 1995, 2( 3 : 171- 178 3 Journel A G Geostatics for conditional simulations of ore bodies J Ecom Geol, 1974, 69( 5 : 673- 687 4 Journel A G, Huijbregts C J Mining Geostatistics M New York: Academic Press, 1978, 600 5 侯景儒, 黄竞先 地质统计学的理论与方法 M . 北京 : 地质出版社, 1990

30、, 6- 7 6 Journel A G, Isaaks E H Conditional indicator simulation: application to a Sakatchewanuranium deposit J M ath Geol, 1984, 16( 7 : 685- 718 7 韩金焱, 孙小彦, 等 163- 171. 8 张泽浦, 王学军 土壤微量元素含量空间分布的条件模拟 J 土壤学报, 1998, 35( 3 : 423- 429 Beijing, China, 9 Xu Ying, Chen Yaxin, et al Study on application of

31、 conditional simulation to time space variability of water salt during soil freezing and thawing A Proceedings of International Conference on Agriculurdl Engineering C December, 1999, II- 129- 133 10 陈亚新, 史海滨, 田圃德, 等 水盐空间变异性监 测的条件模拟 J 水利学报, 2000, ( 6 : 67- 73 11 Dimitrakopoulos R Conditional simulat

32、ion of intrinsic random functions of order K J Math Geol, 1990,12( 3 : 361- 380 12 Oliver D S On conditional simulation to inaccurate data J Math Geol, 1990, 28( 6 : 811- 812 13 Easley D H, Borgman L E, Weber D Monitoring well placement using conditional simulation of hydraulic head J Mathematical G

33、eology, 1991, 23( 8 : 1 059- 1 080 14 陈亚新, 史海滨 渠床土壤入渗率局部 估计点的 Kriging 插值问题 J 水利学报, 1991, ( 2 : 11- 18 矿床条件模拟在矿山生产中的应用 A 中国数学地质 C . 北 京: 地质出版社, 1995, Advances in Hydro Science and Engineering C Beijing : Tsinghua University 252 水 科 学 进 展 第 13 卷 15 Weber A An evaluation of probabilistic procedures for

34、 monitoring well placement that use conditional simulation and pa rameter estimation geostatistical methods A Paris, 1992, 21- 30 16 Dietich C R Computationally efficient generation of ion conditional simulations over regular sample grids J Math Geol, 1993, 25( 4 : 434- 451 17 Hoshiya M, Member Krig

35、ing and conditional simulation of gaussian field J 12( 2 : 181- 186 18 Marcotte M Comditional simulation with data subject to measurement error: Post simulation filtering with modified facto rial Kringing J Math Geol, 1995, 27( 6 : 749- 762 19 Oliver D S Moring averages for gaussian simulation in tw

36、o and three dimensions J M ath Geol, 1995, 27( 8 : 763774. 20 Journel A G Conditional simulation of geologically averaged block permeability J 23- 35 21 陈亚新, 史 海滨, 陈慧新, 等 盐溶 质及其 离子 品位 空间变 异性 的 Co Kriging 估 计 J ( 6 : 35- 40 22 Painter S Numerical method for conditional simulation of levy random field

37、s J M ath Geol, 1998, 30( 2 : 163- 179 23 Cliftun P M , Neuman S P Effects of Kringing and inverse modeling on conditional simulation of the Avra volley aquifer in Southern Arizona J Water Resources Research, 1982, 18( 4 : 1 215- 1 234 24 Dagan G Stochastic modeling of groundwater flow by unconditio

38、nal and conditional probabilities 1 Conditional Simula tion and the direct problem J Water Resour Res, 1982, 18( 4 : 813- 833 25 Davis M W Production of conditional simulations via the LU decomposition of the covariance matrix J Math Geol, 1987, 19: 91- 98 26 Tompson A F, Ababou R, Gelhar L W Implementation of three dimensional turning band random field generator J Water Resour Res 1989, 25( 10 : 2 227- 2 243 水 利学 报,