基于马尔可夫链理论考虑土层变异性的边坡随机场模型

2021-07-21柳彤晖豆红强纪歆雅

水利与建筑工程学报 2021年3期

柳彤晖，豆红强，纪歆雅

(福州大学环境与资源学院, 福建福州 350116)

在岩土工程中，前期勘察对于后续设计以及施工有着举足轻重的影响，通过勘察的钻孔数据获取地下土层的分布情况，才能有效指导并解决后续的岩土工程问题。而获得可靠的土层分布必须依靠大量的钻孔数据，随之而来的便是高昂的勘探成本[1-3]。与此同时，即使岩土工程勘察范围再广，钻孔测井布置再密集，都很难获取连续的地下土层分布，最终只能依据有限的钻孔数据借助当地经验或其它手段进行土层交界面的推断，这会带来很大程度的不确定性。为了利用少量的勘探数据解析地下研究区域的土层分布，刘发全等[4]运用Delaunay三角网和网格生成计算技术建立了三维工程地质模型，Calcagno等[5]基于插值法利用地层交界面的位置和方位数据进行地质建模。Guillen等[6]基于地质数据集合构建了先验模型，提出利用反演方法以概率方式选择出最合适的三维地质模型。祁小辉等[7-8]根据已有钻孔资料设计了不同的钻孔布置方案，建立了表征地层变异的马尔可夫链模型。但这些方法只能得到整个地下区域土层分布的一个最大后验概率估计，而未考虑地层的变异性。换而言之，虽然通过这些方法可以得到具有最大可能性出现的地层分布，但模拟出的地层剖面具有多大的置信度却不得而知。为了经济高效地得到可靠的地层分布情况，需要定量估计地层的不确定性。其中，Elfeki等[9]提出利用单耦合马尔可夫链模型来表征地层结构的非均质性。Strebelle[10]利用多点统计法从训练图像中导出信息用于建立复杂的地质模型。张伟等[11]改良了传统的地质统计学，通过与相控建模的结合，开展了形态和结构复杂的地质体的建模研究。

综上，已有的研究大多以一些特定的地质假设为前提，不能很好地还原真实地层的特性，因此本文将基于马尔可夫链理论建立描述地下土层结构的随机模型，该模型能够反映地层固有的非均质性、各向异性和非平稳性，并在此基础上进一步量化地层的不确定性，从而通过有限的勘探资料定量分析整个模型的不确定性，并进一步探讨了钻孔位置对边坡土层结构随机模型的影响，从而为优化模型提供科学依据。

1 马尔科夫随机场理论

1.1 邻域系统

在马尔可夫随机场中，所研究的地质区域被离散成若干个小的四边形单元，假设S={i1i=1,2,…,n}代表单元集合，其中i为单元索引。S中的单元通过邻域系统相互关联影响，用N={Ni1i∈S}来表示，Ni即称为单元i的邻域，是所有i的相邻单元的集合，如图1所示。单元i的邻域系统Ni包含除了它以外的八个相邻单元{j1,…,j8}，即i∉Ni。邻域系统代表单元之间的相互关系，因此可以利用邻域系统来表示真实地层的空间相关性。

图1 区域邻域系统

1.2 邻域系统

假设R={Ri,i∈S}是定义在S上的随机变量集合，其中每一个Ri在状态集合L={1,2,…,m,…,l}中都有对应的取值ri，这里ri代表岩性单元标记，如碎石土，砂土，黏土等等。记(Ri=ri,i∈S)为随机场的联合分布，简化表示为联合分布事件R=r，其中r={r1,…,rn}为随机场R的一个实现，即随机场的一个地层结构。用R={r=(r1,…,rn)1ri∈L,i∈S}表示所有可能的地层结构。R被称为集合S上关于邻域系统N的马尔可夫随机场，满足：

P(R=r)>0，r∈R

(1)

P(Ri=ri1rs-{i})=P(Ri=ri1rNi)

(2)

式中：S-{i}表示除去单元i的所有单元；rNi表示与单元i相邻的所有单元的岩性标记集合。公式(1)、公式(2)分别体现了马尔可夫随机场的正定性和马尔可夫性，只有满足这两个约束条件的随机场才能称为马尔可夫随机场。其中马尔可夫性是随机场R的局部特性，简单来说就是一个单元的土体类型标记只与其邻域系统内的单元标记有关，其余单元的土体状态均不会对该单元产生影响。

1.3 势函数

为了计算吉布斯分布的概率测度，必须先得到配分函数Z，但由于涉及到大量单元，即使只是解决中等规模的问题，直接计算配分函数Z也是难以实现的。因此利用马尔可夫链蒙特卡罗模拟方法从配置空间R中提取地层实现来得到最大后验概率估计。如前所述，区域邻域系统的条件概率可以基于吉布斯分布通过局部特征的显式形式进行计算：

P(Ri=ri1rs-{i})=P(Ri=ri1rNi)=

(3)

(4)

式中:Λ⊂L表示在区域领域系统中包含的岩性单元标记集合。

势函数Vc(ri,rj)反映了相邻单元的区域空间相关性，且代表了一个吉布斯分布的能量，具有多种不同的形式，本节只考虑双点基团，其结构简单，便于计算。另外考虑到相邻单元有处于相同土体状态的倾向，假定ri=rj有助于降低势。综上，将势函数定义如下：

(5)

式中：β(i,j)为衡量两个相邻单元i和j间区域空间相关性的系数，其大小受单元间相关程度的影响。

考虑到随机场中的单元数量十分庞大，需要一种能在基于区域领域系统的马尔可夫链蒙特卡罗模拟中重复使用的采样方法，使得计算更加简单易行。在此，采用区域转换采样器生成一个初始的地层结构，之后依据土层的非均质性和各向异性不断进行更新，使其更符合真实地质体中固有的空间相关性。

区域邻域系统中的空间相关性根据地质构造的平面走向分为法向相关和切向相关两部分。单元i与其相邻单元j的标准几何关系如图2所示，其中θ为单元i,j的形心连线与笛卡尔坐标系中X轴的交会角。笛卡尔坐标系中的X轴旋转φ可以得到一个以单元i的形心为中心的椭圆，如图3所示，其半径长度ρ(θ)代表单元i和j区域空间相关性，即β(i,j)=ρ(θ)。ρ(θ)越大，相关强度越大，相邻单元j对单元i的影响越大，即两者更易具有相同的岩性单元标记。ρ(θ)可以通过下式进行计算：

(6)

式中：旋转角φ与地层的方位信息密切相关，能够反映区域的非平稳特性。参数a为切向相关强度与法向相关强度的比值，代表了区域各向异性的程度。因此通过选择ρ(θ)函数中不同的参数φ和a就可以代表不同的空间相关性，且该空间相关性即为马尔科夫随机场相邻单元具有相同土体类型的相关强度。

图2 单元i,j的标准几何关系

图3 区域空间相关性的椭圆模型

2 边坡场地随机场的实现

通过随机建模方法模拟边坡土层结构，且要体现真实地层的变异性，就需要一个客观的度量方法来定量分析土层的不确定性，从而确定合理有效的地层模拟次数。

2.1 土层不确定性量化

如前文所述，由于不能对整个地下研究区域进行连续的观测，通过有限的钻孔数据得到的地层分布不可避免地会具有一定的不确定性，对不确定性进行量化才能检验基于研究区域建立的随机地质模型的质量，从而用最低的钻探成本得到最可靠的地层分布估计。目前已有的评价地质模型的方法有很多，包括简单的不确定性测量、统计评价和建模模拟等[12-14]，但大部分都具有一定的局限性，比如评价主要依靠主观估计，或只适用于较为简单的模型。

将地层不确定性可视化实际上需要一种客观量测的方法，这种方法可以度量随机地质模型中任一单元的不确定性，且不受模型条件的限制。因此Wellmann等[15]提出一种利用信息熵量化不确定性，从而进行地质模型质量评估的方法。在此，本文也采用信息熵作为一种客观度量方法来量化地层的不确定性，评估不同模型参数的变化对模型质量的影响。为了定量分析地层的不确定性，引入下式对信息熵进行计算：

(7)

显然，单元的信息熵取决于其可能出现的岩性单元标记数以及其对应的概率。对于已知土体类型的单元，由于该岩土体类型出现的概率为1，显然可以得到其信息熵为0。而不确定单元可能出现多种不同的岩性标记，各个标记出现的概率也不尽相同。基于式(3)—式(7)可以计算出地质模型内每个单元i可能出现的岩性单元标记l以及对应的概率Pl(i)(l∈L)，从而得到其信息熵的值。计算出每个单元的信息熵后可以绘制出信息熵图，从而直观地观察到研究区域内各个单元的不确定性。

2.2 地层模拟的实现次数

由前文可知，单次地层实现的偶然性很大，需要通过多次模拟才能体现地层的变异性，更好地还原真实的地层分布，因此确定合适的地层模拟实现次数极为关键。本文拟考虑通过观察不同模拟次数下总信息熵的波动程度来判断是否能够较好地反映真实的地层分布，从而确定合理的地层模拟实现次数。整个地质体的总信息熵可以通过下式得到：

(8)

式中：1S1表示集合S中包含单元的数量。与单元的信息熵类似，当模型内所有单元都处于相同的土体状态时总信息熵的值为0，若每个单元出现各种标记的概率相同时总信息熵达到最大值。通过计算总信息熵可以对整个研究区域内的地层结构不确定性进行定量分析，从而确定地层模拟的次数。

2.3 马尔可夫随机场的生成

在采用区域转换采样器得到初始地层结构后，其生成的具有最小后验能量的地层分布被称为地层的最大后验概率估计。此后，在所建立的模型框架下所生成的一系列初始结构及其相对应的最大后验概率估计即为可能的真实地层结构

其中，地层最大后验概率估计可以利用吉布斯采样法通过马尔可夫链蒙特卡罗模拟不断更新初始地层结构而得到。通过已有的观测数据生成地层实现的主要步骤如下：

(1) 将所研究的地质区域离散成合适大小的四边形单元，构造邻域系统，建立马尔可夫随机场。

(2) 根据勘探资料赋予钻孔单元和地表单元相应的岩性单元标记，同时利用Kriging法根据已知点的方位信息计算各个单元的方向参数θj→i和φj→i。

(3) 根据单元形心到钻孔位置的距离对相邻单元进行分类，建立吉布斯采样的扫描顺序。

(4) 依照前述所建立的扫描顺序，利用吉布斯采样生成初始的地层结构。

(5) 通过吉布斯采样进一步利用马尔可夫链蒙特卡罗模拟得到初始结构的最大后验概率估计。

3 算例分析

假设有一高10 m，坡度为1∶1的边坡场地，边坡区域内有一虚拟钻孔，根据钻孔资料可知该边坡共分为三层，由上至下分别为残积粉质黏土、全风化花岗岩和强风化花岗岩。将研究区域离散成600个四边形单元，除通过勘探资料可知钻孔内单元的土体类型外，还能直接观察到天然地表单元的土体类型，将以上已知数据赋予对应的单元，可以得到建模区域如图4所示。需要说明的是，由于地层模拟是通过邻域系统来表示相邻单元间的关系，因此统一将所研究的边坡区域离散成较为均匀的四边形单元，避免在坡面附近出现不规则三角形单元，影响计算结果的精度。这种单元离散化使得边坡上部的部分钻孔单元只通过节点相连接，整体钻孔形状没有呈现出规则的矩形长条，但在实际尺寸边坡中的影响可以忽略不计。

图4 建模区域示意图

3.1 迭代次数的确定

地层模拟过程中将最大后验概率估计视为一次实现，其对应的地层结构处于相对稳定的状态。前文已经提到，由于马尔可夫随机场的概率分布难以直接求取，考虑利用马尔可夫随机场和吉布斯分布的一致性，在模拟过程中引入吉布斯分布理论，借助吉布斯分布的能量函数得到马尔可夫随机场的条件概率。因此总能量的波动情况可以作为能量收敛的指标，来判断生成的地层实现是否足够稳定。图5即为边坡土层结构的总能量随迭代次数增长的变化情况。

图5 总能量与迭代次数关系图

从图5可以看出，随着迭代次数的增长，边坡地层的总能量不断降低，且速率较大。当迭代次数达到30次左右之后，总能量依旧保持降低的趋势，但速率明显减小，总能量的波动趋于平稳。这是由于边坡内位于不同土层交界面的单元的岩性标记改变受到了一定程度的限制。值得注意的是，虽然总能量越小，其代表对应的地层分布越稳定，但若总能量继续降低，模拟生成的地下土层结构将趋于具有高度线性的地层交界面，这与真实的边坡地层结构不相符，无法体现出真实地层的变异性。因此本文利用总能量的变异系数作为能量收敛的判定标准，若后续10次采样的总能量变异系数小于0.1%，则认为总能量的变化不明显，趋于平稳，当前迭代次数得到的边坡地层分布可以作为最大后验概率估计，即一次实现。基于上述能量收敛判定准则，可以确定1 000次地层实现对应需要的迭代次数。

3.2 地层模拟次数的确定

单次实现得到的是地层剖面的最大后验概率估计，其没有考虑到边坡地层的变异性，因此需要进行多次地层模拟才能生成所有可能出现的土层分布情况，还原真实的边坡地层结构。通过计算总信息熵可以量化整个边坡地质模型的不确定性，作为判定收敛的指标。与总能量与迭代次数的关系类似，当总信息熵的波动趋于稳定，则认为对应的地层模拟次数是足够的。具体地说，当新增的500次地层实现总信息熵的变异系数小于0.5%，则说明在当前的地层模拟次数下各单元可能出现各种土体状态的概率趋于平稳，即生成的地层实现可以代表边坡场地所有可能出现的地层结构。图6给出了模型的总信息熵随地层模拟次数的变化图，从图中可以看出，随着地层模拟次数的增加，整个边坡地层结构的总信息熵先是急剧增长，之后逐渐趋于平稳。同时为判定模型收敛所需的地层模拟次数，图6也给出了总信息熵变异系数随地层实现次数变化的曲线图。由图可知，总信息熵的变异系数随着地层实现次数的增加不断减小，最后几乎保持不变。当地层模拟次数达到900次左右时，总信息熵的变异系数满足前述的收敛标准，因此取1 000次地层实现来模拟该人工边坡场地的土层分布。

图6 总信息熵与其变异系数

图7则为根据单元信息熵绘制的信息熵图，利用信息熵图可以很好地达到地层不确定性可视化的目的。单元信息熵的大小在模型网格图中用连续变化的颜色来表示。单元信息熵的值越大，说明单元的不确定性越高，对应的网格颜色越接近深色。相反地，信息熵的值越小，单元的岩性单元标记就越固定，对应的颜色也越接近淡灰色。如钻孔单元和地表单元对应信息熵图中的网格颜色均为淡灰色，这就是由于它们都是已知土体类型的固定单元，信息熵为0。而深色的单元都集中在地层交界面的位置，随着到地层交界面的距离增大，网格颜色逐渐从深色过渡到淡灰色。这说明处于土层界面处的单元不确定性最大，距离地层界面的距离越大，单元的不确定性越低。可见，信息熵图可直观地反映单元的不确定性，利用总信息熵作为模拟收敛的判定标准具有可行性。

图7 信息熵图

3.3 钻孔位置的影响

为研究钻孔位置对土层结构的影响，拟在该人工边坡场地坐标2.25 m、7.25 m以及12.75 m处分布布设虚拟钻孔，如图8所示。钻孔布设方案则如表1所示。

图8 钻孔位置示意图

表1 不同的钻孔布置方案

图9即为不同钻孔布置方案下该人工边坡土层结构的信息熵图。由图9可知，三种钻孔方案对应的信息熵图中已知单元的数量和区域由于钻孔位置不同存在明显差异。可见，地层剖面中的发散区域各不相同，不确定性集中的区域也不同。通过对比图9可以发现，钻孔2的位置正好位于方案A对应的信息熵图中不确定性集中的区域。经过计算可以得到方案A、B与C模型总信息熵分别为0.376 3，0.304 6，0.426 1，通过客观度量可以得到三种方案的不确定性大小关系为：方案C>方案A>方案B，即钻孔方案C的不确定性最大，钻孔方案B相对具有最小的不确定性。综合三种方案的信息熵图和总信息熵的情况，推测不同的钻孔布置方案对边坡模型质量的影响可能来源于两个方面。其一是将钻孔设置在不确定性较大的区域可以降低模型的模糊性，其二是钻孔内部单元的数量。因此在建立随机场模型时，钻孔的选择不同也会影响模拟结果的准确性。