随机场离散的Karhunen-Loéve展开法误差分析

刘雪莹, 谭晓慧, 费锁柱, 侯晓亮

(合肥工业大学 资源与环境工程学院,安徽 合肥 230009)

在进行岩土工程的随机有限元分析时,通常利用随机场来表示岩土体参数的空间变异性。一个随机场H(x)就是一个随机过程,其自变量是连续空间位置参数x。给定一个空间位置x0,H(x0)为一个随机变量。随机场的实现就是一个确定的连续函数。随机场模型一般不能在随机有限元分析中直接使用,应先将随机场用有限个随机变量表示,这个过程称为随机场的离散。常用的随机场离散方法有中心点法、局部平均法、Karhunen-Loéve级数展开法(简称“KL展开法”)等[1-3]。其中,采用KL展开法进行随机场的离散能得到连续的随机场(即随机场的离散值在网格边界不会发生突变),而且具有离散效率高、得到的随机变量数较少等优点。因此,KL展开法是一种较为理想的随机场离散方法。

在采用KL展开法进行随机场离散时,必须确定级数的展开项数(又称“截断项数”),KL展开法的展开项数即离散后的随机变量数。理论上,展开项数越多,离散结果越精确,但是,由于随机场离散时计算量也随着展开项数的增加而增加,以及后续的可靠度分析计算量也随之增大,因此,必须合理确定级数的展开项数。文献[4-6]采用随机场期望能的比率因子来确定展开项数;文献[7]采用逐点均方误差的平均值来确定离散误差及展开项数。本文通过理论推导及实例计算,分析了离散误差的影响因素,探究了这2类误差计算公式之间的关系,提出了KL展开法误差计算的新公式。

1 随机场离散的KL展开法

1.1 KL展开法的基本原理

记平稳正态随机场的均值与均方差分别为μ、σ,自相关函数为ρ(x1,x2),则随机场H(x)可展开为:

(1)

其中,ξi为M个互为独立的标准正态变量;M为随机场的展开项数(截断项数);λi为自相关函数ρ(x1,x2)前M个由大到小排列的特征值;fi(x)为与λi对应的特征函数。

将fi(x)用正交基函数hk(x)展开[1,7],即

(2)

其中,dik为特征向量(即特征向量矩阵D的第i列向量)的第k个元素。经变换,λi、dik计算公式为:

AD=DΛ

(3)

其中,Λ为M个特征值组成的对角线矩阵;A为M×M阶方阵,其元素计算公式为:

(4)

其中,j为第n个离散点;hk(x)为fi(x)的正交基函数展开项,可用Legendre多项式函数表示。

k=1,2,…,∞

(5)

其中,Lk-1(·)为标准Legendre多项式函数;T、a分别为平移参数与缩放参数,T=(xmax+xmin)/2,a=(xmax-xmin)/2。此时,(4)式可简写为:

(6)

对于空间二维问题,特征值、特征函数可分别由x、y2个坐标轴方向的特征值乘积、特征函数乘积表示[7-9],即

fi(x)=fi(x,y)=fi1(x)fi2(y)

(7)

常用的自相关函数型式有指数型(exp)与指数平方型(exp2)[1,7],其表达式分别为:

(8)

(9)

其中,Lh、Lv分别为水平、竖直方向上随机场的相关长度(自相关距离);x1=(x1,y1),x2=(x2,y2),表示随机场的空间位置。

1.2 KL展开法的误差

理论上,(1)式中的M应该取无穷大,但在保证计算精度的前提下,通常只截取前M项来提高计算效率。因此,如何定义随机场的离散误差,对于M正确选取及随机场离散精度具有重要影响。

1.2.1 比率法

文献[4-6]建议采用随机场期望能的比率因子确定截断项数,认为比率因子越接近1时离散误差越小,对应的离散误差(本文记为ε1)为:

(10)

其中,S为离散区域的面积。

1.2.2 逐点均方误差法

当采用KL展开法进行随机场的离散时,逐点均方误差[7]可表示为:

(11)

展开(11)式,可得:

(12)

由(12)式可得整个区域上逐点误差的平均值(本文记为ε2)为:

(13)

其中,Ne为离散点的个数。

理论上,由(11)式可得逐点均方误差ε(x)应为0或正值。但计算表明,部分离散点处ε(x)可能为负,这使得采用(13)式求解ε(x)的平均值时,正、负误差中和可能使求得的平均误差ε2偏小,不能反映整个离散区域内误差的真实情况。因此,本文先对每点的离散误差取绝对值,再对所有离散点误差的绝对值取平均(本文记为ε3)来表示整个区域离散误差的平均值,即

(14)

1.2.3 3种误差的对比分析

上述计算误差的3种方法中,ε(x)是逐点误差,与离散点的具体位置有关;ε1、ε2、ε3是离散区域上的平均误差,它们能更好地代表整个离散区域上随机场的离散精度,是决定M的重要依据。

ε1与ε2的关系为:当随机场离散的网格密度较大时,(13)式中的求和符号将变为积分符号;将(12)式代入(13)式,由特征函数的正交性可知,ε1与ε2的计算公式相同,因此,当随机场离散的网格密度足够大时,ε2趋近于ε1,即ε1是ε2的极限形式。

ε2与ε3的关系为:若ε(x)全为非负数值,则ε2=ε3;否则,ε2<ε3。

2 算例分析

2.1 算例简介

本文用2个算例来进行计算分析。

例1 对一个杆件进行随机场离散,属于一维空间问题。该杆件长l1=10 m。

例2 对一个矩形域进行随机场离散,属于二维空间问题。该矩形区域的长度l1=14 m,高度l2=6 m。该离散区域对应于平面应变问题的地基承载力分析,其长度及高度方向分别表示土体中的水平及竖直方向。

2个算例对应的各参数基准值见表1所列。其中,对于二维空间问题,根据文献[10-13]的研究结果,令长度、高度方向(即水平、竖直方向)的自相关距离各不相同,分别为30、3 m。

表1 随机场的统计特征及离散条件基准值

2.2 随机场的逐点离散误差

由(10)～(14)式可知,ε1只与特征值有关;ε2、ε3只与特征值、特征函数有关。因此,当离散区域大小、网格密度、展开项数、自相关函数类型及自相关距离等条件确定时,随机场的逐点离散误差与平均误差是定值,不随随机场离散次数的变化而变化。在表1基准值条件下,例1与例2随机场的逐点离散误差如图1所示。

由图1可知,无论是一维问题还是二维问题,随机场的逐点离散误差都是在区域边界上相对较大。由(11)式可知,ε(x)应为非负数值。从图1可以看出,由(12)式求得的逐点均方误差ε(x)有部分值为负。产生这种现象的原因是:

(1) 在KL展开法中,λi是自相关函数ρ(x1,x2)的前M个由大到小排列的特征值;但fi(x)只是与λi对应的特征函数,它并非由大到小排列,这使得在只截取无穷级数的前M项作为随机场离散值的近似值时,若前M项中含有的特征函数值fi(x)较大,而相应的λi又不是很小,则逐点均方误差ε(x)就有可能为负。

(2) 由(2)式可知,特征函数值与正交基展开函数h(x)有关,而后者又与标准Legendre多项式函数及离散点的具体位置x有关。例如,M=4,k为1～M时,标准Legendre多项式Lk-1(u)=Lk-1((x-T)/a)的函数图形如图2所示。

从图2可以看出,对于离散域边界上的点,Lk-1(u)恒为1;对于离散域边界附近的点,Lk-1(u)的绝对值也较大。因此,离散域边界上及边界附近,hk(x)的绝对值随着k的增加而增加(见(5)式),导致特征函数fi(x)也较大(见(2)式),这使得离散区域边界及其附近离散点上的逐点均方误差ε(x)极易为负值。

图1 随机场的逐点离散误差

图2 k取不同值时标准Legendre多项式曲线

2.3 随机场离散误差的影响因素

以表1的参数值为基准值,分别改变离散区域大小、离散网格密度、展开项数M、自相关函数类型及自相关距离,研究其对随机场平均离散误差ε1、ε2、ε3的影响。分析某个影响因素时,只改变此影响因素的相应值,而其他计算条件仍取表1的基准值。

2.3.1 离散区域大小及离散网格密度的影响

分别取离散区域的边长为基准条件的1～3倍(对于二维算例,矩形区域的2条边长等比例增加),再对2种离散区域,分别取随机场网格边长为0.25、0.50、1.00 m。离散区域大小、离散网格密度对随机场离散误差的影响分别如图3、图4所示。

图3 随机场离散误差与离散区域大小的关系

图4 随机场离散误差与网格密度的关系

图3中,横坐标l1/Lh表示离散区域的相对大小。从图3可以看出,无论是一维算例还是二维算例,3种误差均随着离散区域的增大而快速增加;ε1与ε2较为接近(ε1略大于ε2),ε3较大。图4中,横坐标所示的网格边长越小,则网格密度越大。从图4可以看出,ε1与网格密度无关;ε2随网格密度的增加而有所增加,网格越密,ε2越接近于ε1;ε3随网格密度的增加无明显变化规律。图4结果与1.2.3节的理论分析结果一致。

对比图3、图4的纵坐标可以发现,图3纵坐标的数值变化很大,图4纵坐标的数值变化很小。因此,离散区域大小对随机场离散误差的影响很大,网格密度对随机场离散误差的影响则很小,相比之下可以忽略不计。实际上,随机场离散的KL展开法是一种级数展开法,其离散结果是连续的,与随机场的网格密度无关。上述计算结果中,ε2、ε3与网格密度有关,这是由于ε2、ε3只计算了特定离散点处的平均离散误差(见(13)式、(14)式)。

2.3.2 展开项数的影响

在采用KL展开法进行随机场的离散时,展开项数M是关系到离散精度的重要因素。M越大,离散误差越小;但是,随着M增加,随机场离散的计算工作量也大幅增加,计算速度变慢。因此,需要确定合理的M。基准值条件下,2个算例的随机场离散误差与M的关系如图5所示。

图5 随机场离散误差与M的关系

(1) 当M较小时,ε1、ε2明显小于ε3,因此,根据ε1、ε2、ε3是否小于某个误差限值来确定M时,3种方法得到的M值并不相同。例如,设误差限值为ε0=5%,对于例1,由ε1<ε0及ε2<ε0可得M=3,由ε3<ε0可得M=4;对于例2,由ε1<ε0可得M=4,由ε2<ε0可得M=3,由ε3<ε0可得M=4。

(2)ε1、ε2、ε3总体上随M增加而快速减小,但二维算例的ε2在M=3时明显偏小,这是由于此条件下逐点误差的过度正负中和而导致的,因此,ε2不适用于判断M。由于ε1是ε2的极限形式,因此,ε1也不适用于判断M。ε3能综合反映离散域上各点误差的算术平均值,是判断M的较好判别标准。

2.3.3 自相关函数类型、自相关距离的影响

自相关函数类型、自相关距离大小对随机场离散误差的影响如图6所示。

图6 随机场离散误差与自相关函数及自相关距离的关系

因为ε2不适用于判断M,而ε1是ε2的极限形式,所以只绘了ε3与自相关函数类型、自相关距离大小的关系。由图6可知,exp2对应的离散误差明显小于exp对应的误差,因此对于一个具体的研究对象,如果对随机场的自相关函数类型没有明确的试验数据,可以采用exp2,其对应的离散误差偏小;自相关函数类型相同时,自相关距离越大,随机场离散误差越小,这是由于自相关距离越大时,同样条件下对应的自相关函数值越大(见(8)式、(9)式),由(4)式得到的矩阵A的元素值越大,其对应的特征值也越大,因此离散面积不变时,随机场的离散误差就越小。

3 结 论

(1) 对于随机场离散的KL展开法,可以采用3种方法求解离散域上的平均误差。误差ε1是误差ε2的极限形式;误差ε3可以合理反映逐点均方误差为负的情况,是判断随机场展开项数M的一种较好的误差判别标准。

(2) 离散区域大小增加时,离散误差ε1、ε2、ε3增加;离散网格密度增加时,ε1不变,ε2趋近于ε1。

(3) 指数平方型自相关函数对应的离散误差明显小于指数型自相关函数对应的误差。

(4) 自相关函数类型相同时,自相关距离越大,随机场离散误差越小。