基于序贯高斯模拟法的水利工程坝址区速度建模

李伟, 周怀来*, 刘兴业, 王元君, 马慧莲

(1.成都理工大学地球物理学院, 成都 610059; 2.成都理工大学地球勘探与信息技术教育部重点实验室, 成都 610059)

全球对水利工程设施的建设需求的不断增加,勘探力度逐步加大。中国拥有众多河流,发展前景良好。在追求发展速度的同时,越来越注重水利等工程项目的质量。水库作为水利工程中的基础设施之一,为区域居民的生活和生产提供了重要的水资源保障[1]。在修建水库之前,坝址的选择至关重要,必须对影响水库后期投入使用的各种因素进行必要的分析,也包括地质灾害分析工作等[2],在全面符合条件的情况下才能进行规划设计和施工,以确保水库的正常投入使用。在坝址的选择中,分析坝址区的速度和岩性是重要的两个参考因素[3]。通常,坝址区岩层速度的横向变化范围大,仅依靠测井资料的垂向速度信息很难作为依据,需要利用其他条件加以约束[4]。通过建立速度模型或岩性模型可用于分析坝址区地下结构的稳定性,然而这对建模成像精度提出了更高的要求[5]。

水库堤坝进行选址勘察时,会对速度及地下岩石的性质进行重点勘察,速度分析是必不可少的一步。岩石越坚硬致密,测井上的速度响应相对越高,而在岩石松散或孔隙较大部位,速度呈现异常低值,这对水利工程的修建以及后期的运行形成垮塌的风险隐患。分析不同类型岩石的速度对堤坝类型和选址的影响也是重要的一步[6]。对坝址区建立地下岩石的速度模型,是分析坝址稳定性的一种方法。然而这对建模精度有很高的要求,常规的建模方法很难实现,此问题待解决。为得到更高精度的近地表速度模型,学者们开展了许多的研究。赵玲芝等[7]运用多信息融合的方法建立近地表速度模型,将近地表调查资料与大炮初值层析反演结果进行协克里格插值技术融合,得到更高精度的速度模型并在静校正及叠前深度偏移中得到了很好的效果。刘灵等[8]提出一种相控速度建模方法,综合地震速度、测井速度、地震属性及沉积相等多种信息,通过多属性融合及神经网络聚类分析来建立低频相控速度模型,并用于约束反演结果,极大提高了速度建模的精度及分辨率。

水利工程设施的建设对地下速度建模精度要求严格[9],而常规的速度建模思路是利用测井资料得到数据进行插值,多为线性插值方法[10],对于水利坝址区此类复杂地层结构,常规插值方法难以解决问题,此外井数据只包含纵向速度信息,对于井之间的横向大范围空白区域缺乏约束条件。综上所述常规插值方法比较简单[11],建模效果不佳,需寻找高精度的速度建模方法。为此,有学者提出了利用地震层位约束来进行高精度建模的方法。在地震勘探领域,孟会杰等[12]为解决地震偏移成像精度不高的问题,应用层位约束的网格层析建模方法,通过不断迭代、反演最终建立地下深度域高精度速度模型,并消除了实际工区中的构造成像假象。后来,有学者在对地质统计学建模的研究中发现,相比于常规线性插值建模,基于地质统计学理论的建模方法在效果上更好,精度更高,能够充分利用变量的空间变化特征[13]。

为提高建模精度,准确分析水利坝址区内部构造特征,服务于坝址的选取及水库的建设,鉴于前人专家学者的研究,现提出一种基于地质统计学理论的地震与测井相结合的速度建模方法,对常规的测井插值建模的方法以及地质统计学建模方法进行完善和补充,应用地质统计学序贯高斯模拟随机建模,以地震层位信息作为横向约束条件联合测井资料建立速度模型,为坝址区岩石稳定性提供数据支撑,对选坝提供重要参考。同时,在前人对速度建模方法研究的基础之上,创新性地将本文提出的基于地质统计学理论的建模方法应用于水利工程领域的坝址选取工作上。

1 方法原理

岩石速度是坝址选取的重要影响因素,因此以速度作为关键切入点对水利坝址区建立速度模型,进而评估待选坝址区的稳定性。常规速度建模方法较为简单,模型精度无法满足水利工程建设的要求,而使用基于地质统计学理论的地震与测井相结合的方法建立高精度速度模型,旨在解决坝址区速度建模精度低的问题,在高精度建模结果上圈定有利于修建坝址的区域,更加准确且具有说服力,为选坝修坝提供重要依据。

地质统计学作为一门发展于20世纪60年代的数学地质边缘学科,包含空间相关性分析、克里金插值以及随机模拟三部分[14]。将研究区域在未知位置的函数值作为一组相关随机变量,通过随机变量理论分析来预测未知点的信息。相较于线性的插值方法有独特的优势特征,且超越了普通的插值问题[15],这也是利用地质统计学进行速度建模的理论依据。

地质统计学提出区域化变量[16]这一概念,能用空间分布来表达具有结构性的随机变量。同时,为度量区域化变量的空间变异性,提出了变差函数的概念,它兼顾区域化变量的随机变化和空间结构性规律[17],表示空间中数据方差随着两点之间相对距离的变化而在特定方向上的变异性[18]。假设随机变量ξ满足本征假设,则在任一方向α,相对距离为|h|的两个随机变量ξ(u)和ξ(u+h)的增量的方差称为变差函数γ(h),表示为

(1)

式(1)中:u为空间中的任意点;h为两点之间的距离,在地质统计学中称为滞后距;E为数学期望。在平稳假设条件下变差函数γ(h)的变化只与滞后距h有关,而γ(h)随h的变化称为变差函数曲线,如图1所示。图1中c0表示滞后距h很小时,两点之间具有一定的变异性。a表示变程,指区域化变量在空间中的相关范围,变程越大说明观测数据在更大范围内相关。超出变程范围后,两点之间无关;在变程内,变量越近,变异性越小。c为基台值,代表变量在空间中总变异性的大小。在速度建模过程中,对变差函数参数进行分析是至关重要的一步,其结果直接影响建模效果,所以需要对参数进行细致的分析。

图1 变差函数曲线Fig.1 Variation function curve

通常实际观测获取的样品数目有限,例如在速度建模过程中能提供的钻井曲线信息较少,而通过有限观测值经式(1)计算获得的变差函数作为实验变差函数γ*(h),表示为

(2)

式(2)中:N(h)为滞后距为h的点对个数;ξ为随机变量;ui为空间中的任意点。

以h为横坐标,γ*(h)为纵坐标即可得到实验变差函数图。实际应用中,可供利用的点对数量越多,获取的实验变差函数就具有越强的代表性。

实际的地质统计学研究用到的都是理论变差函数,因此需要在获得实验变差函数之后进行变差函数拟合以得到理论变差函数的参数[19]。其中球状模型、指数模型、幂模型及高斯模型等是一些基本理论变差函数。根据实际原始数据特点选取相应的变差函数模型,不同模型的插值结果也会不同。本文研究使用球型变差函数模型,其灵活性大、数据适应性强,模型在原点处呈线性变化,在变程处达到基台值。图2为球型、高斯型、指数型这三类常用理论模型的变差函数曲线。

图2 常用变差函数模型Fig.2 Common variation function model

使用球形模型进行变量空间估计时,需定义球体变差函数的长轴、短轴和中长轴的变程以及块金值、拱高和基台值等。同时定义球体的空间走向,作为变差函数的方位参数。本文的速度建模方法将测井得到的速度数据视作地质统计学中的区域化变量,对其进行变差函数分析,研究其数据特征选择球型变差函数模型,求取最佳的变差函数参数。

确定变差函数重要参数后,即可进行地质统计学序贯高斯模拟。它是随机模拟的一种方法,以已知信息为约束,使用随机函数建立多个不同地质模型的方法。在模拟过程中考虑区域化变量随机性的同时,还考虑了其空间相关性。且要求模拟结果忠实于观测点的值[20]。假设随机变量符合高斯分布,结合高斯概率理论与序贯模拟算法对连续分布的区域化变量进行模拟,以像元为单位按照选定的随机路径进行,并且将已经模拟过的节点的状态值也直接加入条件数据中。模拟步骤如图3所示[21]。序贯高斯模拟方法建立速度模型,相比于常规的建模方法在原理上有明显优势,它能充分利用实际观测数据的空间信息,并对其进行变差函数分析以此建立更为精准的速度模型,为坝区提供高精度的地下速度信息,作为实际选坝的主要参考。

序贯高斯模拟可将地震层位作为约束条件加入模拟结果进行随机建模。而无层位约束的模拟结果无法反映地下岩层的横向变化特征,会得到一个相对光滑、连续的模型结果。水利坝址区地下构造因其速度在横向变化范围大,若无横向方向的约束将缺失建模的精度,最终模型结果无说服性。因此文章将地震层位作为约束进行速度建模,提高建模的精度与准确性。

序贯高斯模拟法在多次模拟时会生成多个等概率的模型,其属于条件模拟方法,因此模拟结果忠于观测点处的数据,而对于非观测点区域,模型之间差异较大。可对多次模拟得到的速度模型之间有明显差异的区域做不确定性分析,选坝时则应避开此类位置。

2 建模步骤

(1)对测井速度曲线进行粗化处理,提高垂向采样使数据更加连续,同时进行井曲线的标准化,使其更加符合地质规律,便于后续的速度插值与建模。

(2)分析变差函数参数。根据已有数据特征和精度,选取球形模型变差模型进行建模,确定球形模型主、次、垂向变程,以及最优的块金值、拱高和基台值等重要参数。

(3)对选取的变差函数参数进行序贯高斯模拟,首先进行无层位约束模拟,再以地震层位作为横向约束进行模拟,完成不同层位的变差函数分析和分层速度建模。

(4)使用协克里金法将上述有、无层位约束情况下建立的速度模型进行加权[22],增强模型的纵向连续性以及准确性,形成最终的速度模型。并与常规的建模方法进行对比、评价。

(5)将最终得到的模型使用Garden公式计算出密度,得出反射系数,与雷克子波进行褶积得到合成地震记录,与实际地震记录相对比,评价速度建模的效果。

3 应用实例

在对以上方法研究的基础上,对位于四川省蒙顶山某一处的近地表钻井资料进行实际数据测试,该区域地层速度横向变化范围大,难以获得精确的地下速度分布,因此需在高精度速度建模的结果上进行选址分析。工区范围内分布11口可用井(W1～W11),井位信息如图4所示,井位分布比较均匀。首先对工区内井曲线做批量标准化处理,得到的速度直方图如图5所示。图6为井W6～W10 5口井标准化后的结果。此外对井速度数据做“粗化”处理,取速度值的算术平均对测井曲线进行重采样,增强垂向上的连续性。

图4 坝区井位信息Fig.4 Well location information of dam area

图5 井曲线标准化前后直方图Fig.5 Histogram of well curve before and after standardization

图6 井曲线(W6～W10)标准化前后结果Fig.6 Well curve (W6～W10 ) results before and after standardization

完成上述两步骤后进行序贯高斯模拟,若以地震层位作为横向约束,则需要对每一层变差函数进行分析。表1为有层位无层位约束的变差函数参数的分析结果。本文变差函数模型选择球形模型,由于水平方向变差函数参数对模拟结果影响很小,因此将其视为各向同性,只对垂向变差函数参数进行分析。

表1 有层位、无层位约束变差函数参数Table 1 Variation function parameters with horizon constraints and without horizon constraint

为分析不同变差函数对模拟结果的影响,下文使用无层位约束的模拟方式并通过设置不同的块金值n(0、0.15、0.30、0.45)及垂向变程v(5、25、45、65)进行测试,如图7和图8所示。

图7 不同块金值的模拟结果Fig.7 The simulation results of different nugget values

图8 不同垂直变程的模拟结果Fig.8 The simulation section results of different vertical ranges

通过分析,块金值越大,变量的随机性表现得越强烈,当块金值为0时,模拟结果平滑,能较好体现出研究区域的速度变化趋势,随机性较小;对于垂向变差函数分析结果,不同垂直变程对模拟结果的差别较大。取较小的垂直变程能够提高垂向分辨率,但对区域整体趋势刻画不足;若取较大的垂直变程,能体现出区域整体的变化趋势,却缺失垂向分辨率。因此选取最佳垂向变差参数才能达到较好的模拟结果,此处选取25作为最佳垂直变程。

完成变差函数分析后,使用序贯高斯模拟方法得到速度模型。首先图9(a)展示无地震层位约束的速度模型过井W6～W10的剖面,图中5条黑色横线为地震层位,由上至下依次命名为S1～S5(S1、S5分别对应t=0 ms和t=90 ms时刻)。从图中得出由于没有层位的约束,在分层处速度横向变化特征未能突出。图9(b)为有地震层位约束后的模拟结果,在对每个层位的井数据都做了独立的变差函数分析后,得到了相比无层位约束更高精度的结果,同时能反映出横向速度变化特征,但其速度模拟结果在局部层位分界处之间变化剧烈,与常规地质规律认识不符。

图9 无层位约束与有层位约束模拟结果过井W6～W10剖面Fig.9 Simulation section results without and with horizon constraint cross well W6～W10

将无层位约束的速度模拟结果与有层位约束的模拟结果进行协同克里金加权插值,加权步骤中将有层位约束的结果设置为主变量,无层位约束的模拟结果为次变量,通过分析确定权值后进行插值最终得到图10所示的速度模型。对模拟结果进行分析,模型既能体现出速度的横向变化的复杂性,也符合地质认识,且相较于上文结果精度进一步提升,为水库的坝址选取分析以及安全性评价提供重要的数据支撑。

图10 协同克里金序贯高斯结果过井W6～W10剖面Fig.10 Co-kriging sequential gaussian section results cross well W6～W10

本文研究中使用的序贯高斯模拟属于随机模拟方法,多次模拟能得到多个等概率的模拟结果,但在钻井位置均能忠实于观测数据。多次模拟结果的差异体现了建模结果的不确定性。每次模拟后可对速度值差异明显的区域进行分析和取舍,以此缩小范围确定最佳的选坝位置。图11展示使用协克里金序贯高斯再次模拟两次的结果。

图11 两次协克里金序贯高斯模拟结果过井W6～W10剖面Fig.11 Twice co-kriging sequential gaussian simulation section results cross well W6～W10

结合图10,在剖面上红色区域圈定出了一处三次模拟结果未发生明显变化且速度响应较高的区域,说明该位置的模拟结果响应更趋近于真实地下构造。在水库修建中,基于安全考虑选坝位置在速度高值响应特征的岩层中最为合适,而此分析结果为坝址选取提供了一项重要的参考依据。

下文分析常规速度建模方法效果。图12(a)为克里金插值法得到的速度模型,其结果得到的是一种相对较平滑的结果,局部的速度异常点未得到突出,只能用作反映全局特征。其次克里金插值是一种确定性方法,每次插值得到的模型相同,无法反映模型之间的差异性以及不确定性。

图12 克里金插值与井速度散点插值结果过井W6～W10剖面Fig.12 Kriging interpolation and well velocity scatter interpolation section results cross well W6～W10

图12(b)是井速度散点直接插值为速度体后的结果,图中的黑色竖线为钻井提供的速度观测值。分析后得出插值效果并不理想,速度呈明显条状分布,且插值结果不完整,模型出现大片空白区域以及速度值为0的异常位置,不宜作为最终速度模型的参考。

通过以上分析,使用地震层位约束的地质统计学协克里金序贯高斯模拟得到的速度模型精度最佳,能较好地展示坝址区地下岩石的速度信息。图13为t=63.60 ms时刻的速度模型平面图,通过分析速度在平面上的分布特征,发现速度高值响应大多分布于坝区靠南区域,图13中红色实线为上文展示的过井剖面图在此平面的投影,蓝色框线部分圈定了一处有利修建水库的坝址区,区域内速度值较高,同时在过井剖面上为高值特征响应,其结果为水利工程坝址区的选址提供了重要的参考依据。

图13 平面上t=63.60 ms时刻协克里金序贯高斯法速度模型Fig.13 Co-kriging sequential gaussian velocity model on the plane at time t=63.60 ms

以图10速度模型作为最终结果,利用Garden公式计算密度及反射系数,并与主频为50 Hz的雷克子波进行褶积得到合成地震记录,检验建模效果,如图14(a)所示。

图14 合成地震记录与实际地震记录过井W6～W10剖面Fig.14 Synthetic seismograms and the actual seismic record section cross well W6～W10

将计算得到的合成记录与实际地震记录图14(b)进行对比,从中分析得出:合成记录与实际资料的吻合程度较好,地震层位对应地震剖面上的同相轴部分连续性较好,通过本文的方法建立的速度模型对坝址区的稳定性分析及选址都能够提供重要的依据,达到了速度建模的目的及实际效果。

本文提出的基于地质统计学理论的高精度速度建模方法圈定预测出的实际区域有利坝址区,已被当地相关水利施工单位认证,在实际现场勘测中,区域内地下岩性普遍致密,速度呈现高值响应。坝址位置选取正确,建设方案已被采取,现正进入施工阶段。

4 结论

(1)地层的速度信息是影响坝址区选择以及后期水利工程安全投入使用的重要因素。因此,对速度进行分析是必要的一步,有必要运用近地表速度模型的精细构建技术来确立精确近地表速度模型。

(2)针对水利工程坝址的选取,本文提出基于地质统计学随机模拟的速度建模方法,并引入地震层位横向约束建模过程。根据坝址区近地表速度分布的特点选取合理的变差函数模型是确保模型精度的重要前提,若参数选取不当则得到的速度模型无法用来评估坝址区的稳定性以及后期投入使用。利用测井数据作为硬约束,地震层位作为横向约束的地质统计学序贯高斯模拟方法建立的速度模型符合地质规律。

(3)多次模拟可以得到多个等概率的速度模型,不同模型之间的差异可以反映建模结果的不确定性,为坝址选择的风险决策提供依据。将提出的方法应用于实际数据资料上并取得了较好的成果,依据建立的速度模型分析并圈定出了有利于修建水利工程的坝址区,圈定方案被采用,实际正在施工当中,为水库的选坝的合理性及后期的投入使用提供了重要的数据支撑与保障。