土性参数各向异性随机场的表征与建模方法

王占盛，陈 健，戎虎仁，杨定强

(1.中铁二院工程集团有限责任公司，成都 610083； 2.中国科学院武汉岩土力学研究所，武汉 430071；3.山西大学 土木工程系，太原 030006； 4.中铁第六勘察设计院集团有限公司，天津 300308)

1 引 言

土性参数具有非均值特性，即在空间上具有一定的变异性[1]。常见土性参数(如强度参数、变形参数和水力参数等)的空间变异性根源于土体的沉积过程或后沉积过程，后沉积过程包括局部成岩作用以及剪切带或不连续变形带的形成过程等[2]。因此土性参数的空间变异性是一种固有特性，这种特性对岩土工程问题的影响非常值得探讨。

现阶段岩土工程问题的分析多是基于均值参数展开，取得了丰硕的成果。但考虑参数空间变异性之后，会对问题分析结果带来怎样的影响，这方面的研究还有很多的工作要做。尤其是当土性参数在高度非均值的情况下，自然土体或工程岩土结构的稳定性可能会受材料的不连续或材料性质的弱区域主导，此时如果再采用均值参数进行分析，可能结果就会与实际情况有很大的偏差。

自认识到土性参数具有空间变异性以来，学者们开展了大量的研究工作。Vanmarcke[3,4]提出土体剖面的随机场模型，采用波动范围(scale of fluctuation)来描述土体参数的空间变异性，之后又完善了随机场理论。Fenton等[5]认为土体参数相关距离取决于土体沉积过程中的地质作用。文献[6-10]对岩土参数的波动范围进行了大量的实验和统计研究，对工程实践极具参考作用。冷伍明[11]改进了相关距离计算的递推空间法。傅旭东等[12-14]探讨了相关距离不同计算方法的差异。程强等[15,16]对相关函数法计算相关距离的方法和相关函数形式的选择、拟合范围等问题进行了探讨，并通过大量CPT资料发现了一般土层相关距离的分布范围。李小勇等[17,18]根据工程勘察资料和室内强度试验，分析了太原典型粘土强度指标的空间变异性，并对几种不同方法求得的参数相关距离进行了对比分析。韩宪军[19]在土性相关距离的计算中考虑了统计不确定性，并探讨了它对计算结果的影响。

另一方面，为分析参数空间变异性对工程的影响，学者们引入了可靠度分析方法，获得了大量的成果,常见的可靠度分析有估计法(一次二阶矩法[20]、概率矩点估计法[21,22]和响应面法[23]等)和模拟法。随着计算机技术的快速发展，有限元蒙特卡洛分析方法开始越来越广泛地应用于可靠度分析[23-25]。

在土性参数相关距离研究方面，由于垂直向地层参数的获取较水平向地层参数的获取更容易，因此相关距离的研究多集中在一个方向的讨论上，而事实上在任意方向上都应有对应的相关距离值，即在空间上变异性具有各向异性。这在很多工程现场都能有一定的直观感受。图1所示为文献[2]给出的几种常见相关模式的实拍图，其中图1(a)为各向同性，即参数在每个方向上的变异特征是一致的；图1(b)为横观各向同性，即参数空间相关性具有两正交主轴，其中沿长主轴方向的参数相关程度较高，变异程度较低，而沿短主轴方向的参数相关程度较低，变异程度较高，常见于水平沉积、河床演变后的自然土体，以及水平固结后的工程土体；图1(c)为旋转横观各向同性，是横观各向同性两相关主轴旋转一定角度后形成的新模式；图1(d)为一般各向异性，该模式下两相关主轴不正交，且其中一相关主轴平行于坐标轴，常见于受断层运动影响的沉积土、土体崩积、颠覆、端倾斜坡处土体以及具有两组或多组弱面的土体等；图1(e)为旋转一般各向异性，是旋转各向异性的相关主轴旋转一定角度后形成的新模式，无相关主轴平行于坐标轴；图1(f)为复合型相关模式，即土体不同区域相关模式不相同。

为反映参数各向异性对岩土工程的影响，必须建立合理的参数反演模型。而随机场理论为此目的实现提供了理论基础。目前，适用于各向异性随机场建模的方法主要有局部平均划分法LAS和矩阵分解法MD,本文将在反演精度与反演效率两方面分析的基础上对比两种方法的优劣，并基于此对矩阵分解法进行有益的改进。

图1 土性参数常见相关模式

2 土性参数各向异性随机场的表征

对于图1的几种土体相关模式，如何从数学上进行精细的描述,需从相关等值线(面)谈起。

2.1 相关等值线(面)

空间某一点参数值与一定点参数值间的相关程度可用相关系数来衡量，参数x(t)的自相关系数函数定义为

(1)

式中τ为空间相对位置，E为数学期望，mx和σx分别为参数x(t)的均值和标准差。根据相关系数函数的定义，有

(2)

式(2)表明，自相关系数函数是偶函数。在自相关系数函数基础上，定义相关距离θ为

(3)

相关距离是衡量空间变异性的一个特征值，对于各向异性相关模式，其值在不同方位上会有所差异。

2.2 相关距离函数

(4)

(5)

对于横观各向同性情况，有

(6)

(7)

对于其他的相关模式，也可采用类似的方法，求得相关距离函数。因此采用相关距离函数可以对各向异性问题进行很好地描述。

表1 五种自相关结构

图2 空间相对位置

3 各向异性随机场常用反演方法分析

3.1 局部平均划分法LAS

文献[27]给出了横观各向同性随机场LAS法反演的流程，以图3为例，具体如下。

(1) 选取水平向相关距离作为任意方向上的相关距离，生成各向同性随机场，此时异性系数ξ=1，如图3(a)所示。

(2) 将步骤(1)生成的随机场垂直方向压缩到原来的1/3，此时各向异性系数ξ=3，如图3(b)所示。

(3) 将步骤(2)生成的随机场水平方向拉升到原来的3倍，此时各向异性系数ξ=9，如图3(c)所示。

在步骤(1～3)中，也可先基于垂直向相关距离生成各向同性随机场，然后再在水平向和垂直向进行或拉升或压缩操作得到相应的各向异性随机场。

图3 LAS方法横观各向同性随机场生成过程

3.2 矩阵分解法MD

该方法将空间域内所有网格的参数值集合成一随机向量X，则其反演模型可采用方程

X=m+CU

(8)

式中m为X的均值向量，U为脉冲场，C为待确定的系数矩阵。则X的自相关矩阵为

B=Cov(X,X)=CCΤ

(9)

根据自相关系数函数性质，B为对称阵，可采用Cholesky分解方法求得C，故此法称为矩阵分解法。

在矩阵分解法中，相关矩阵包含了所有相关信息，因此极为重要，以图4所示的二维随机场网格为例，其对应的相关矩阵为

(10)

图4 随机场网格

3.3 LAS与MD的比较分析

(1) LAS与MD各向异性随机场反演方法对比

如前所述，LAS方法反演各向异性随机场是在同性随机场基础上进行拉伸或压缩方法得到的，由于拉伸与压缩操作只能在水平或竖直方向进行，因此该方法适用于横观各向同性随机场，即适合近似水平地层的模拟，而对于模拟其他相关模式的各向异性随机场并不适用。

MD方法中，各向异性特征集中体现在相关矩阵上，因此在构建此矩阵时，只要采用相应的各向异性相关函数对矩阵对应位置元素赋值即可。因此，MD方法能处理任意相关模式的各向异性问题。

(2) LAS与MD反演效率与精度分析

从反演方程的形式上看，LAS只考虑相邻网格对目标网格的相关性影响，MD则需考虑所有网格对目标网格的相关性影响。以图4所示的网格为例，当对目标网格(3,3)进行反演时，LAS考虑的相关性影响范围内的网格有(2,2)，(2,3)，(2,4)，(3,2)，(3,3)，(3,4)，(4,2)，(4,3)和(4,4)等9个，而MD考虑的相关性影响范围内的网格包括了从(1,1)到(m，n)的所有网格。也正因此差异性，两者在反演效率与精度上有显著差异。

在反演效率上，相关矩阵Cholesky分解起到最关键的影响。一个网格数为n的区域构建的相关矩阵，其Cholesky分解可能进行的浮点运算次数为[29]

Ο(Cholesky)=n3/3

可见，相关矩阵Cholesky分解运算次数随着网格规模的扩大呈几何级数上升，而LAS方法由于只考虑相邻网格相关性的影响，因此其系数求解过程中的运算次数并不随网格规模的扩大而上升。因此，MD较LAS反演效率更低，且这种差异会随着反演规模的扩大而急剧增大。

在反演精度上，由于MD方法中相关矩阵的构建考虑了所有网格对目标网格的相关性影响，故经此方法反演得到的随机场相关性统计值与理论值是吻合的。LAS由于只考虑相邻网格对目标网格的相关性影响，故经此方法反演得到的随机场相关性统计值与理论值有一定偏差，图5为Fenton[28]给出的LAS方法相关性统计结果。该结果是尺寸5 m×5 m、网格数256×256、相关距离θ=0.5 m的各向同性随机场反演100次统计得到的，可以看出统计值与理论值之间有一定的偏离且不收敛。

图5 LAS方法相关性统计结果

4 改进的矩阵分解法

4.1 反演效率与反演精度平衡的关键

由前文可知，无论反演效率还是反演精度都与相关性考虑范围密切相关，因此选择合适的相关性考虑范围是平衡两者的关键。为找此平衡点，需从相关距离的物理意义与相关函数的特性着手。

式(3)为相关距离的定义式，可以看出，相关距离是相关系数函数在全空间域上的积分，故其是空间相关性整体的量化值。常用的相关函数有高斯型、指数型函数、二阶自回归函数和指数余弦函数等。以最常用到的高斯型函数为例，其具体表达式为

(11)

图6所示为高斯型相关函数曲线，可以看出，相关系数随τ/θ的增大而快速衰减，当τ/θ=3时，相关系数值已趋近于0。其他类型相关函数也具有类似的快速衰减特性。

相关函数的快速衰减特性说明，空间相关性集中体现在相关特征尺度(相关距离)一定倍数范围内，如果随机场反演时相关性考虑范围选择合适，或者说相关性考虑程度取舍合适，则不仅能保证随机场的反演精度，还能保证随机场的反演效率。

图6 高斯型相关函数曲线

4.2 矩阵分解法的改进

从上面的分析可知，相关性考虑范围是平衡反演效率与反演精度的关键。以图7所示的网格为例，当采用高斯型相关函数时，如果2倍网格尺寸大于3θ时，则只考虑图7阴影部分的相关性影响已足够，否则需进一步扩大考虑的网格范围。

如果2倍网格尺寸大于3θ，参数随机场X1可采用如下反演方程

(12)

对于{C11}Τ的求解，一方面可基于图7阴影部分构建部分相关矩阵，并对此矩阵进行乔列斯基分解获得转化系数向量，当相关考虑范围较大时，采用此种方法误差很小。

图7 相关考虑范围

(13)

以此类推，可建立如下相关性条件方程组

(14)

此方程组为多元二次方程组，可采用牛顿迭代法获取系数向量{C11}Τ。

4.3 各向异性随机场反演流程

根据本文所述相关结构数学表达与随机场反演模型，确定各向异性随机场反演流程如下。

(1) 根据反演精度要求，确定相关考虑范围。

(2) 求解系数向量{C11}Τ。

(3) 生成脉冲场U1。

(4) 利用式(14)生成随机场X1。

5 改进的矩阵分解法反演结果分析

5.1 改进的矩阵分解法与LAS方法反演精度

采用与图5所示统计结果相同的模型和参数，采用改进的矩阵分解法，得到相应的随机场相关性统计结果如图8所示。可以看出，改进的矩阵分解法相关性统计值与理论值的吻合度较LAS法更高，这也说明本文的改进思路是可行且有效的。

5.2 各向异性随机场反演结果

反演精度统计结果表明采用改进的矩阵分解法生成随机场是适合的，本节采用此方法生成了横观各向同性随机场与旋转横观各向同性随机场，结果分别如图9和图10所示。

(1) 图9中，随机场均值为0，方差为1，模型尺寸50 m×50 m，网格数100×100，垂直向相关距离2 m。图9对应的各向异性系数分别为1.0,3.0 和6.0，水平向相关距离分别为2 m，6 m和12 m。图中颜色越浅的区域，其值越小，颜色越深的区域，其值越大。可以看出，在垂直向相关距离确定的情况下，随着各向异性系数的增大，水平向带状分布越明显，即参数在大范围空间内取相近值的可能性提高。

图9 横观各向同性随机场反演结果

(2) 图10中，各向异性系数为3.0，其余随机场基本参数与横观各向同性随机场基本参数相同。

图10 旋转横观各向同性随机场反演结果

图10对应的旋转角度分别为0°，30°，60°，90°，120°和150°。可以看出，在垂直向相关距离、各向异性系数确定的情况下，旋转角度决定了带状分布的角度，这适合非水平地层随机场的反演。

6 结 语

综上所述，本文结论如下。

(1) 从空间相关等值线(面)出发，可采用相关距离函数描述参数各向异性相关的问题。

(2) LAS方法只适合地层水平或地层垂直情况的随机场反演；而常规的矩阵分解方法虽然能处理非水平地层的情况，但其会受到相关矩阵规模的限制，并进而影响到反演的效率。

(3) 基于相关函数的快速衰减特性，本文改进了矩阵分解法，该法首先要根据反演精度的要求确定相关考虑范围，这也是该方法质量的保证。相关性统计结果也表明，改进的矩阵分解法较LAS方法的反演精度要高。

(4) 采用改进的矩阵分解法，生成了横观各向同性与旋转横观各向同性两种模式的随机场，表明此方法能处理随机场较复杂的各向异性问题，有较好的适用性。