不同地层渗透性聚类和分割方法研究

不同地层渗透性聚类和分割方法研究

1地层-地下水聚类分析聚类分析是一种使用多元统计理论对其进行分类的方法。假定有n个样品(对象、事件),每一个样品包含p个变量(指标),要求根据这n个样品的p个变量指标进行聚类。变量的选取依赖于聚类的目的,在对地层的水文地质分层中主要是依据土壤的渗透系数及其分布特性来聚类,即将各层的土壤渗透系数作为变量。土壤的渗透性不仅在平面上分布有差异,而且在垂向上分布也不均一,因此,水文地质分层既要对大量的地质剖面数据在平面上聚类(分区),又要对同一类剖面在垂向上进行分层,从而得到合理的概化结果。1.1p维欧氏空间中n个点的聚类设通过水文地质勘探和调查,可得到n个水文地质剖面资料,每一个剖面为一个样品,而每一个剖面又包含p层不同渗透系数的土壤,每一层的渗透系数即为一个变量。以xij表示第i个水文地质剖面中的第j层土壤的渗透系数,i=1,2,…,n;j=1,2,…,p。每个剖面(样品)包含有p个土层的渗透系数(变量)。将任一剖面的p层渗透系数(xi1,xi2,…,xip)视为p维欧氏空间中的一个点,n个剖面便组成p维欧氏空间中的n个点。在平面上对n个剖面进行分类就是对p维欧氏空间中的n个点进行聚类,自然是土壤渗透性较近的地层剖面应该聚为一类。各地层剖面的渗透系数是否相近,可以用欧氏空间中的距离来表示。其表达式为dij=√pΣk=1(xik-xjk)2。(1)对p维欧氏空间中的n个点进行聚类,常用的有系统聚类法、快速聚类法等。系统聚类法的基本思想是逐步将距离近的类合并在一起。首先将n个样品各自看成一类,即有n个类,此时类间距离即为样品间距离。把距离最小的两个类合并成一个新类,于是便得到n-1个类;重复上述步骤,每次至少合并一个类,如此一直进行到所聚的类别间的距离满足合适的精度要求且有最少的类别数为止。快速聚类法适用于对大样本进行快速聚类,尤其对类的特征(各变量值的范围)有了一定认识时,用此方法分类速度较快。这种方法需要事先确定要求聚成的类数,如要求聚为k类,则要首先选择k个样本作为聚类的种子,亦即k个聚类中心点,或称为初始类中心;按照分别距这k个类中心的距离(可用欧氏距离计算)最小的原则将样本分派到各类中心所在的类中,形成第一次迭代得出的k类;计算每一类中各样品的均值,作为第二次迭代的聚类中心点,再按照距类中心的距离最小的原则将样本分派到各类中心所在的类中,形成第二次迭代得出的k类;依次迭代下去,直到满足所要求的精度为止。1.2类内类直径l计算在地层剖面中,各层土壤的相对位置是固定不变的,因此,各层土壤的渗透系数是有序排列的样品,在分层时其次序不能打乱,属于按序聚类问题。记有序样品x1,x2,…,xm,其下标1,2,…,m反映了样品的次序。对有序样品聚类后,每一类的形状应是{xi,xi+1,…,xj},其中1≤i≤j≤m。由此可见,对有序样品的聚类分析实质上是把样品按原有次序分割成若干类(段),若将m个有序样品分成k类(段)(k=2,3,…,m),则有Ck-1m-1种不同的分法,目的是找到一种最优的分割方法,使该分割的方差最小。为此,需要给出类的均值和类的直径。若分割后某类(段)Gij由{xi,xi+1,…,xj}(i<j)组成,该类的均值为ˉxij=1j-i+1jΣl=ixl。(2)类Gij中j-i+1个样品的离差平方和,表示类Gij内各元素偏离其类平均值的程度,称之为类直径,可表示为D(i,j)=jΣl=i(xl-ˉxij)2。(3)设将m个有序样品分成k类(k=2,3,…,m)的某一种分法P(m,k)为{xi1,xi1+1,⋯,xi2-1},{xi2,xi2+1,⋯,xi3-1},⋯⋯⋯⋯⋯⋯{xik,xik+1,⋯,xik+1-1}。可简写成{i1,i1+1,⋯,i2-1},{i2,i2+1,⋯,i3-1},⋯⋯⋯⋯⋯⋯{ik,ik+1,⋯,ik+1-1}。其中分点为1=i1<i2<…<ik<ik+1-1=m。这种分类的目标函数为Ο[Ρ(m,k)]=kΣj=1D(ij,ij+1-1)。(4)如果确定了m与k的值,则O[P(m,k)]越小表示各类的离差平方和越小,亦即分类越合理。因此要找一种分法P(m,k),使目标函数达到最小。从目标函数的表达式出发,可导出如下递推公式:Ο[Ρ(m,2)]=min2≤j≤m{D(1,j-1)+D(j,m)};Ο[Ρ(m,k)]=mink≤j≤m{Ο[Ρ(j-1,k-1)]+D(j,m)}。(5)由公式(5)可以看出,要得到将m个样品分成k类的最优分割,必须先求出将m-1个样品分成k-1类的最优分割,而后者又必须先得出将m-2个样品分成k-2类的最优分割,依次类推。最优分割的具体计算步骤如下:(1)计算所有可能分类的直径D(i,j),其中i=1,2,…,m-1;j=2,3,…,m;(2)分别计算将i个样品(i=3,4,…,m)分成2类、3类、…时的最优分割目标函数,即Ο[Ρ(i,2)]=min2≤j≤i{D(1,j-1)+D(j,i)},i=3,4,⋯,m。Ο[Ρ(i,3)]=min3≤j≤i{Ο[Ρ(j-1,2)]+D(j,i)},i=4,5,⋯,m。⋯⋯⋯⋯⋯⋯Ο[Ρ(i,k)]=mink≤j≤i{Ο[Ρ(j-1,k-1)]+D(j,i)},i=k+1,k+2,⋯,m;k=3,4,⋯,m-1。(6)(3)合理确定分类数k。根据对样品分成k类、k-1类、…、2类等各种情况,计算最优分割目标函数值,满足上述最优分割目标函数值的分割方式,就是所要求的k类最优分割。2地层渗透系数计算实例对某堤防工程的堤基依据85个勘探剖面的资料进行地层渗透性聚类分析。所有剖面均垂直堤轴线布置,平均间距为0.5km,地质资料较详细,并具有较好的代表性。根据堤防加固工程设计的需要和当地的地质条件,勘探的重点是第四系和第三系松散沉积层,除少数钻孔外,堤身上的钻孔深度一般都在35m以内,堤身以外的钻孔深度一般都在30m以内。为了研究堤防加固以后对沿堤附近地下水运动的影响,需要对有防渗墙和无防渗墙条件下地下水渗流场进行计算分析。为此,要求对复杂多变的地层渗透性分布进行合理聚类分层,概化出几种有代表性的地层剖面类型用于渗流计算。实例分析中所采用的地层渗透系数,取自对堤防勘探所揭露的不同岩层的岩性分析,即在当地针对不同岩层而进行的渗透系数试验结果,见表1。由于不同剖面所揭露的深度不同,每层厚度各异,为方便聚类分析计算,将各剖面点的实测地层渗透系数分为若干个等间隔的离散值,每一间隔称为一层,离散层数的多少根据计算精度要求确定。设nd为离散后的总层数,每层厚度为dzi(i=1,2,…,nd)。对85个原始实测地层剖面的渗透系数,自地表向下分为nd=27层(按照实际勘探剖面揭露的天然地层的最深值确定),第i层厚度dzi=1m,总深度27m。3聚类分析结果所选样本(观测剖面)共85个,样本量较大,分类数也可预先估计设定,故采用快速样本聚类分析方法对地层渗透系数进行平面聚类。这里每个样本是一个观测剖面,它含有27层土壤渗透系数值,故每一个样本都是一列数组。确定了分类数以后,程序可自动生成各类的类中心数组,进行聚类计算。通过计算得出了将全部85个剖面聚为2类、3类、4类、5类、6类等情况下的平面聚类(分区)结果。从不同聚类数所对应的每一类内各变量(剖面)数组距类中心数组(层与层相对应)的均方差(表2)可以看出:分的类数越少,每一类的样本量越大,误差就越大,反之,误差就越小;聚为2类与聚为3类比较,前者的均方差明显大于后者,垂向各层的平均均方差相差达8.15;聚为3类与聚为4类比较,均方差的差别明显减小,为4.44;而将聚为5类与聚为6类比较,其均方差的差别就微乎其微了,只有0.49。可见,聚为2类时的精度太差,不可取;而聚为6类时的精度提高不明显,故没有必要分得太细;因此可认为3,4,5类的分法比较合理。对这3种分类的比较分析看出:按3类进行聚类,就可以得到渗透性的基本分类框架,当然,分的类数越多精度就越高,合理的分类数需要按照精度要求和实际情况确定。就该地区堤防沿线的地层渗透性特点而言,聚为3类或4类比较合理。在聚为同一类的所有观测剖面中,取各对应层的渗透系数的平均值作为聚类以后该层的渗透系数值,每一类仍为27层。经聚类以后各分类的渗透系数分布如图1～图7中的虚线带圆点的符号所示。4聚对层序地层渗透系数的影响经过平面聚类以后,对聚为3类、4类情况下的每一类剖面的岩层渗透性进行垂向分层。根据最优分割原理,分别对每一类进行按层最优分割计算,分割的层数根据实际情况分析确定。与平面聚类类似,分的层数越多,精度越高,但也不宜分层太多,否则给计算工作带来麻烦,同时层数过多,难于给出地层类型的宏观概念。经分析,认为分为5～8层比较合理。当在平面上聚为3类时,各类剖面岩层渗透性的最优分层结果见图1～图3中的折线。由图1～图3可以看出,垂向分层的结果较好地概化了原有岩层的分布:第一类剖面地层的渗透性沿深度上的变化呈现先小后大的趋势,15m深度以上,渗透系数在7.8m/d范围内变化,为细质砂层与粘壤土互层,具有相对隔水性;15m以下,渗透系数迅速增大,直到35m/d,为很好的含水砂层。第二类剖面属于明显的二元结构地层,从地面向下到19m深度范围内,渗透性较差(2～5.9m/d);19m深度以下,渗透系数呈阶梯状迅速增大。第三类剖面为单一结构地层,渗透系数从0.8m/d到7m/d逐渐变化,没有明显含水砂层出现。当在平面上聚为4类时,各类剖面最优分层结果绘于图4～图7。对比图1～图7,可以看出,聚为3类情况的第二类(图2)与第三类(图3)在这里基本保留(为