概率密度演化方程差分格式的计算精度及初值条件改进

朱志辉，刘禹兵，高雪萌，周高扬，余志武

(1. 中南大学土木工程学院，湖南，长沙 410075；2. 中南大学高速铁路建造技术国家工程研究中心，湖南，长沙 410075)

工程结构服役期内会受到各类随机动力作用[1]，结构响应也表现出明显的随机性。开展随机动力作用下的结构可靠度分析对结构性能和安全性评估至关重要。近年来，由李杰和陈建兵[2-3]提出的广义概率密度演化理论(Generalized probability density evolution theory, GPDET)为结构的随机响应预测及其可靠度分析提供了新的途径。

由于实际工程的复杂性，GPDET 的偏微分方程很难获得精确理论解析解，因此常采用数值方法求解[4]。这些数值算法包括有限元法[5]、路径积分法[6]、等效线性化法[7]、有限差分法[8]等。当采用适当的差分格式时，有限差分法在概率密度演化分析中表现出良好性能[9]。目前常用的概率密度演化方程(PDEE)有限差分法按照差分近似微分的处理方式分为单边差分、L-W 双边差分和总变差减小(TVD)性质差分[10]三种数值差分方法。单边差分和L-W 双边差分法分别利用单向差分和中心差分来近似微分结果，两者分别具有一阶和二阶计算精度。而TVD 差分法由于引入了通量限制器，避免了L-W 双边差分法可能存在的非物理负数概率解。因此，李杰和陈建兵[10]构造了一类适用于PDEE 求解的TVD 差分格式并将其沿用于其后开展的各项概率密度演化分析。

GPDET 分为概率分配和PDEE 求解两大关键步骤。概率分配对计算精度影响的研究[11-12]已较为成熟。在概率分配满足计算精度的前提下，YU等[13]基于GPDET 开展铁路桥梁在随机荷载条件下的结构响应分析时，通过对比数论选点法与蒙特卡洛法的激励样本计算结果，验证了车轨桥随机振动模型的正确性；XU 等[14]基于GPDET 将车辆—轨道相互作用模型与轨道不平顺概率模型有效耦合，同样利用数论选点法与蒙特卡洛法计算结果对比，验证了耦合模型的正确性。以往研究中，GPDET 与蒙特卡洛法的对比结果仅体现了概率分配步骤中代表性样本对样本空间的表征性，而两种方法对比的均值和标准差均采用TVD 差分法计算。由于该差分方法计算精度有限，且在局部极值点附近精度降低[15]，这也导致概率密度演化分析中[16-17]，在荷载作用的后半程(卸载阶段)，结构响应的概率密度演化标准差结果依旧持续增大，这样的计算结果显然不满足精度要求，甚至不符合实际情况。因此，PDEE 求解步骤中产生的计算误差不容忽视，开展基于TVD 差分法的PDEE数值求解精度分析非常重要。

通常情况下，PDEE 的初值条件为在空间方向可导性较差的脉冲型狄拉克(Dirac)函数，这将影响到概率密度数值解和可靠度分析的精确性。为提升TVD 差分方法的求解精度，石晟等[18]提出了余弦型初值条件，余弦型初值虽比脉冲型初值平滑，但空间方向的初始概率赋值范围受限于余弦函数周期，且未给出时间差分步长选取原则。

本文针对PDEE，首先对单边差分、L-W 双边差分和TVD 差分三种差分方法开展单样本和多样本的差分精度分析；进而提出了TVD 差分的时间步长合理选取方法和正态分布型初值条件，通过与蒙特卡洛法的数理统计均值和方差作对比，验证了本文方法的合理性。

1 概率密度演化方程及数值求解方法

1.1 概率密度演化方程

结构响应量Z和随机变量 Θ组成的增广系统(Z,Θ)是一个保守的随机系统，随机性完全来自于Θ，整个过程中，没有随机因素的增减。将(Z,Θ)的联合概率密度函数记为pZ,θ(Z,θ,t)，根据概率守恒原理[19]，对于给定的离散参数代表点集θl(l=1,2,···,nsel)，PDEE 可写为一维偏微分方程形式：

第l组离散代表点集对应的边界条件和初始条件为：

式中： dz为 空间离散步长；nz为空间离散点数量。

1.2 数值求解方法

式(1)的偏微分方程很难给出理论解析解，故通常采用数值差分方法进行求解计算。通过对概率密度值p进 行时间t方向的一阶和二阶泰勒展开[9]，利用差分方法近似p的一阶和二阶偏微分，可得：

1.3 差分格式的数值耗散

由于式(5)和式(6)利用差分代替微分得到线性算子Aγ，由式(17)可知，差分格式、结构响应方向上各时刻概率密度函数p的可导性、结构响应一阶时间导数Z˙ 、差分求解次数k(或时间长短)和概率初值条件p0均对数值差分结果有直接影响。且任一时刻的概率密度值都可由初值条件p0j迭代获得。

整理得到单边差分与双边差分概率密度演化算子Bo和Bd的表达式分别为：

当 |λZ˙|＜1时，双边差分的概率密度演化算子Bd比单边差分的Bo更接近1。由式(18)可知，只考虑一阶泰勒展开的单边差分法，其数值耗散大于考虑二阶泰勒展开项的双边差分法。结合式(14)可知，TVD 差分法的数值耗散情况介于单边差分与双边差分之间。

2 数值差分时间步长及正态分布型初值条件

2.1 时间步长选取

数值差分的稳定性条件为CFL 条件(Courant-Friedrichs-Lewy condition)：

由式(9)和式(19)可得空间离散步长的要求为：

式(20)表明：只要响应离散步长 dZ足够大，便可保证概率密度演化计算收敛。而当工程结构受地震[20]、轨道不平顺[21]和风载荷[22]等剧烈变化的随机荷载时，则需要更小的时间步长与空间步长来提升概率密度演化结果计算精度。为呈现出概率密度演化结果的峰值变化，并保证结果误差符合工程要求，必须使各样本的概率演化曲线在“时间—空间”网格上有明显的重叠和分离部分，这就对网格离散精度提出了要求。

有限差分法求解过程示意图如图1 所示，竖轴代表响应量，样本A 和样本B 分别代表样本空间中某两条以通量速度为Z˙A和Z˙B在“时间t—空间Z”网格上进行概率密度演化的曲线。

图1 中，区域1 的两样本差异小于 dZ，两样本的概率密度值将在同一空间离散步左右相邻网格节点上表征，致使两样本结果差异性较低，随着类似情况的增加，误差不断累积，最终的概率密度分布范围也将比实际情况更宽(标准差结果失真)。同时，由于概率守恒(各时刻概率之和恒为1)，概率密度峰值也将比实际情况更小。因此，只有当空间离散步长 dZ较小时，各样本的变化差异性才能很好的体现出来(如区域2)。

使用TVD 差分法求解PDEE 时，由于通量限制器φ取值范围为[0, 1]，故TVD 差分法的单样本概率密度演化曲线的离散性只能随差分迭代次数的增加逐渐增大或维持离散性不变，这也将导致最终的概率密度演化结果(或标准差结果)误差较大，甚至是错误的。

式中： |Zˆ|为关注时间段的各样本结构响应差值范围；T为结构响应分析总时长。

为同时满足计算稳定性与精度要求，改进时间步长 dt*应为：

式中：T为概率密度演化分析总时长； Ψ为保守系数，各样本响应结果差异越小，保守系数取值越大。

2.2 正态分布型初值条件

在评价可靠度时，通常将获得关注量的响应与规定限值比较，令“W=规定限值 -关注量响应”。通过假设等价极值事件Wsin(t)，将该虚拟三角函数的一阶时程导数Wcos(t)代入式(1)，即可得到虚拟随机过程的概率密度演化结果。当t=π/2 时 ，等价极值事件与实际过程W相等[23]。

由式(23)可知，概率密度pZ的求解精度会直接影响后续的可靠度分析结果，通常各样本对应的PDEE 初始条件为式(3)所示的Z0处的脉冲型Dirac 函数，而Dirac 函数初值条件在空间方向振荡明显，这导致初值条件关于空间方向的可导性对数值计算结果产生影响[18]。

为进一步提高TVD 差分的计算精度，本文基于正态分布函数提出了一种如图2 所示的初值条件：

3 差分格式计算精度的影响因素

3.1 样本精度影响分析

3.1.1 单样本数值差分精度

为对比不同差分格式对结构可靠度计算精度的影响，构造目标三角函数y(t)：

式中：A为三角函数幅值，取值为2；t为时间，选取范围为0 s～3.15 s。

选取差分时间步长0.0001 s，按式(4)计算单样本的均值和标准差(该统计特征值无实际物理意义，仅用于判断演化结果的离散程度)，结果见图3。算要求(0.01 s)的概率结果离散性明显大于0.0001 s的标准差结果。当差分步长小于0.004 s 时，TVD差分在第一个周期内的标准差均小于5%。按式(22)计算的时间步长为 dt*≤0.05/(5×2×5)=0.001 s，该步长能很好地保证TVD 差分的单样本计算精度。图5(b)中，双边差分标准差均符合误差限值要求，这表明相同误差限值要求下，双边差分可采用更大的差分步长，提升计算效率。

3.1.2 多样本数值差分精度

图6 为不同差分时间步长的多样本统计特征值。图6(a)中，差分时间步长对多样本均值不产生影响。图6(b)中，TVD 差分的标准差结果受差分时间步长影响显著，0.0008 s 和0.01 s 时间步长的 π/2时刻处标准差值分别为0.098 56 和0.1147，与精确解0.1 的误差分别为1.44%和14.7%，随着概率的演化，各差分时间步长的标准差相对于精确解均出现了增大现象，其中 π时刻精确解为0，而0.01 s 步长计算的标准差为0.071 22，计算结果严重偏离精确解。图6(c)中，由于双边差分的单样本的数值耗散远小于TVD 差分，故其多样本标准差结果也更接近于精确解。双边差分法以0.0008 s 和0.01 s 时间步长计算的 π/2时刻处标准差值分别为0.099 69 和0.0987，与精确解0.1 的误差分别为0.31%和1.3%。因此，当计算结果仅关注均值和标准差等统计指标时，使用数值耗散较小的双边差分法，可在较大步长条件下获得高精度的统计指标，并大幅缩减计算时间。

3.2 改进初值条件

型初值条件零时刻的概率分布较为尖锐，概率密度函数(Probability density function, PDF)峰值为12.85，而图7 的精确解初始时刻PDF 峰值为7.97。

选取原始的时间步长为 dt=0.01 s，根据式(27)和式(28)即可得到100 条目标样本曲线。这100 条随机样本的均值和标准差分别按传统数理统计方法和数值差分法计算，结果见图9。

图9(a)中， dt=0.01 s的时间步长能够满足构造样本精度要求，而对于 dZ≥1.0075的CFL 条件，过大的空间步长将导致无法进行差分计算。这进一步表明当结构响应一阶时间导数过大时，仅满足结构响应精度要求的时间步长并不能满足概率密度演化方程的差分求解精度。图9(b)中，初值条件在分布范围方向上的平滑性很差，这也导致图9(b)的标准差结果在初始时刻大幅震荡，进而导致之后0 s～2 s 内的标准差都严重偏离精确解；数值差分标准差曲线后半段的离散性逐渐增大是因为通量限制器的存在导致概率密度分布范围只能增大不能减小引起的。

图10 正态分布型初值条件标准差 σ=dZ时的概率密度分布结果。由于采用了正态分布型初值条件，其初始时刻PDF 峰值为11.78，与图8 的Dirac 函数型初值条件相比，PDF 峰值降低，分布范围方向上的可导性提高，使概率密度分布结果更接近图7 的精确解。

以差分解计算误差的时间平均值Erro作为计算精度指标：

式中，Ei为第i时间点均值或标准差的误差。

使用不同 σ值( σ 以0.2 dZ为间隔)的正态分布型初值条件开展数值计算，结果见图11。图11 结果表明：改进的正态分布型初值条件可通过调节标准差 σ达到减小均值和标准差计算误差的目的。随着 σ的增大，均值累计误差一直减小；而标准差累计误差先减小后增大， σ=dZ时达到最小值。与式(27)精确解的初始概率相比，按式(3)Dirac 函数型初值条件累加求和得到的初始概率分布更集中(即峰值更大，分布范围更窄)，因此，正态分布型初值条件在 σ=dZ之前， σ的增大导致初始概率分布的峰值逐渐降低，分布范围逐渐变宽，初始概率分布向精确解逼近；随着 σ的进一步增大，虽然初始概率分布的峰值依旧向精确解靠近，但分布范围却比精确解更大，导致标准差累计误差不再减小。

4 结论

本文分析了单边差分、L-W 双边差分和TVD差分三种差分方法在不同差分时间步长条件下的样本离散性，指出了差分精度的影响因素，提出了TVD 差分法的时间步长合理选取方法和正态分布型初值条件，通过数值算例验证了时间步长选取方法和正态分布型初值条件的准确性及普适性，并得到如下结论：

(1) L-W 双边差分法的单样本概率演化离散性最小，误差允许范围内可采用的时间步长最大，当计算结果仅关注均值和标准差等统计指标，而不考虑概率密度分布情况时，建议使用L-W 双边差分法开展高效的数值求解。

(2) 通过对比三种差分方法的数值耗散情况，得到数值差分计算精度的影响因素包括：结构响应方向上各时刻概率密度函数的可导性、结构响应一阶时间导数、差分求解次数 (或时间长度)和概率初值条件。

(3) 基于TVD 差分开展PDEE 数值求解时，空间离散步长不仅要满足CFL 条件要求，还应满足各样本差异性的最大值限制，本文提出的差分时间步长选取方法有效降低了TVD 差分的单样本离散性和多样本概率密度求和导致的累加误差，并且差分时间步长可根据样本差异程度灵活调整。

(4) 在传统的脉冲函数型初值条件基础上提出了一种更平滑且可以通过标准差 σ控制的正态分布型初值条件形式。当 σ取值较小时，该分布条件近似于脉冲型Dirac 函数初值分布。合适的标准差σ可以提高正态分布型初值条件的数值差分精度。