以木结构为主的古建筑是我国的瑰宝，具有很高的历史、文化和艺术价值。然而由于木材本身的性质以及外界因素的影响，长期处在自然环境下的古建筑木结构会受到不同类型和程度的损伤，严重影响结构的安全性和稳定性。因此，为了实现对古建筑木结构的预防性保护^[1]，首先需要掌握构件的力学性能时变规律，然后对一定环境下构件的残余强度和剩余寿命等功能函数进行分析，最后跨越到以安全性为主的结构服役状态的评估。自19世纪四五十年代以来，研究人员就对木材强度随时间变化的规律展开研究。在得到木材强度受长期荷载作用随时间衰减的经验曲线的基础上，Gerhards等^[2-3]基于累积损伤概念建立了长期荷载下木构件强度模型。考虑影响木构件承载力的其他因素，又分别建立了不同残损影响下的木构件损伤时变模型：李铁英^[4]对中国古建筑木结构中的腐朽数据进行分析和整理，研究得到木构件腐朽层厚度随时间变化的发展规律，国内关于腐朽对木构件力学性能影响的众多研究均以此为分析基础；谢启芳等^[5]通过去除腐朽部分来模拟木柱的局部残损，基于轴压试验和有限元模拟来分析局部损伤木柱的受力性能退化规律。考虑木材虫蛀，澳大利亚的CSIRO团队^[6]根据现场调研获得的大量数据，建立了白蚁攻击木构件的风险演化模型，瞿伟廉等^[7]基于该模型提出了理想环境下木材虫蛀与时间的数学模型，并据此推导出木构件的抗力衰减模型。而在木材干裂方面，基于对木材试验数据的拟合分析，彭磊等^[8]提出了木构件开裂深度和含水率之间的变化关系；陈孔阳^[9]通过力学试验、有限元数值模拟和理论分析，提出考虑纵向裂缝影响的梁、柱承载力计算公式；谢启芳等^[10]对不同干裂程度的木柱进行轴压试验，考虑稳定系数、截面面积损失和初始偏心距3个影响参数，建立了自然干裂木柱的受压承载力和刚度退化模型。综合考虑多因素，王雪亮^[11]和秦术杰^[12]建立了受长期荷载、虫蛀、腐朽和干缩裂缝共同影响的木构件强度退化随机时变模型。由上所述，国内外对残损木构件强度退化的研究多针对单个影响因素，综合考虑多因素的较少。

同时，影响构件残损的因素，如环境温湿度、构件所用材种材性、荷载大小、承载面积等均具有随机不确定性，仅用“安全”或“失效”来评定其服役状态是不合理的，基于可靠度理论对木构件进行安全性评估是较为科学的分析方法。可靠度分析的传统算法有一次二阶矩法、二次二阶矩法、蒙特卡罗法和响应面法等^[13]。对传统算法进行改进和结合，又衍化出适应不同情况的新算法。蒙特卡罗法（Monte Carlo method，MC）是可靠度分析较为常用的方法，但在参数较多的情况下，计算效率低，且不能保证结果的收敛性。近年来，针对随机结构动力系统分析，李杰等^[14]基于有限单元法基本原理提出了概率密度演化方法（probability density evolution method，PDEM），弥补了传统可靠度算法存在的不足，在多个领域得到广泛应用，为可靠性分析开辟了新的道路。

如前所述，古建筑木构件的损伤是多因素综合作用的结果，需要考虑影响因素的随机性进行可靠度分析。因此，本研究考虑多种损伤，参考已有的单因子损伤时变模型建立木构件多因子损伤时变模型，并推导木构件强度退化模型；同时以一古建筑木结构的典型柱为例，基于概率密度演化方法对柱的残余强度进行可靠度分析，与蒙特卡罗法的分析结果进行对比，验证概率密度演化方法用于残损木构件可靠度分析的可行性。具体技术路线如图1所示。

图  1  技术路线图

Figure  1.  Technology roadmap

下载: 全尺寸图片 幻灯片

1.   古建筑残损木构件强度退化时变模型

古建筑木构件在多因素的共同影响下产生损伤，而长期荷载、腐朽、虫蛀和干缩裂缝是对木构件力学性能影响较大的损伤^[12]，本章将基于前人研究建立的单因子损伤时变模型，构建综合考虑这4种残损的多因子损伤时变模型及强度退化模型。

1.1   单因子损伤时变模型

本节分别列出了仅考虑单一因素的木构件损伤时变模型，得到各因素影响下构件剩余有效截面积的时变模型，可作为构建多因子损伤时变模型的基础。

1.1.1   长期荷载

以Gerhards提出的累积损伤模型^[2]为基础，推得木材强度随时间变化的预测模型^[15]。

式中：${\varphi }_{\mathrm{M}}$为构件残余强度比，t为服役时间（a），f为长期荷载下木材的残余强度（MPa），${f}_{0}$为短期静强度（MPa），$\alpha \left(t\right)$为随时间变化残余寿命的损伤率，其取值范围为[0,1]，$\alpha \left(t\right)=0$时为无损状态，$\alpha \left(t\right)=1$时完全失效；$\sigma$为外力作用下截面的最大应力（MPa），${{{B}}}_{1}$、${{{B}}}_{2}$为参数，需要通过试验进行标定。由式可知木材残余强度比是短期强度${f}_{0}$、应力比$\dfrac{\sigma }{{f}_{0}}$和时间t的函数。

1.1.2   腐　朽

根据李铁英总结出的古建筑木构件腐朽层厚度估计式^[4]，得到腐朽影响下构件截面积随时间变化率的估计式。

式中：$\xi$为木构件的截面积变化率，${\delta }_{0}$为构件当前腐朽层厚度（mm），d为构件截面直径或边长（mm）， ${T}_{0}$为预期已服役时间（a），t为时间（a），$\alpha$为考虑腐朽层厚度发展速度的指数参数：当$t\leqslant 400$时，$\alpha =1$；当$400 < t\leqslant 800$时，$\alpha =1.5$；当$t > 800$，暂无数据。构件截面积变化率是木构件当前腐朽层厚度${\delta }_{0}$与截面直径（或边长）d之比、预期服役时间${T}_{0}$和时间t的函数。

1.1.3   虫　蛀

基于CSIRO团队建立的随机虫蛀模型^[6]，将虫蛀过程分为初始期${t}_{\mathrm{l}\mathrm{a}\mathrm{g}}$和发展期${t}_{\mathrm{d}}$两个阶段，得到虫洞引起构件的截面积损伤变量^[16]。

式中，${D}_{\mathrm{i}\mathrm{n}\mathrm{s}\mathrm{e}\mathrm{c}\mathrm{t}}\left(t\right)$为虫蛀引起的截面积损伤变量；a为建筑开始被虫蛀的时间（a）；b为系数；${t}_{\mathrm{l}\mathrm{a}\mathrm{g}}$为虫蛀初始期总时间（a），包括在建筑50 m范围内建立一个成熟蚁巢的时间${t}_{1}$ （a）；白蚁形成虫道并到达建筑的时间${t}_{2}$ （a）；白蚁突破或绕开所设置防蚁屏障的时间${t}_{3}$ （a），分别通过专家对各阶段内影响虫蛀事件的评估以及环境的优劣综合计算得到。

1.1.4   干缩裂缝

假设裂缝均为图2所示的张开TR型^[17]规则三角形裂缝，基于木构件相对开裂深度的时变模型^[11,18]以及损伤变量的定义，得到受干缩裂缝影响有效截面积随时间变化的预测模型。

图  2  干缩裂缝类型

TR指裂纹面的方向为弦向，且裂纹沿径向扩展。图片引自参考文献[17]。TR means the direction of the crack surface is tangential and the crack extends along the radial direction. The picture is cited from reference [17].

Figure  2.  Shrinkage crack type

下载: 全尺寸图片 幻灯片

当木柱截面为圆形，则

当木柱截面为方形，则

式中：${D}_{\mathrm{c}\mathrm{r}\mathrm{a}\mathrm{c}\mathrm{k}}\left(t\right)$为干缩开裂引起的截面积损伤变量，是随含水率和时间变化的函数。$\lambda$为裂缝相对宽度和相对深度的比值，其取值范围一般为$\lambda \in [0.15, 0.25]$^[12]；E 、F 为简化系数， E = C₁ + C₂C_MC，F = C₁C₂C_MC；C₁、C₂、k为参数，需通过试验标定；${C}_{\mathrm{M}\mathrm{C}}$为相对平衡状态下的含水率。

1.2   多因子强度退化模型

由上述单因子模型的分析可知，在各种因素的共同作用下，损伤累积使得构件有效截面积减小，应力随之发生变化。因此，基于累积损伤理论^[2]，通过逐年累加得到多因素作用下古建筑木结构残损构件的残余寿命损伤率随时间变化的关系式。

式中：$\tau \left(i\right)$为木柱服役第i年时的应力比$\dfrac{\sigma \left(i\right)}{{f}_{0}}$。此外，本节所有公式各参数含义均与上节相同。

本文以古建筑木结构中柱的抗压强度退化为例，假设各影响因素相互独立，且受损伤部份完全失效，综合以上单因子模型，对考虑长期荷载、腐朽、虫蛀、干缩裂缝的残损木柱损伤时变模型及强度退化模型进行推导。

假设木柱（图3）轴心受压，在顺纹方向承受轴向压力N （N），服役初期初始截面积为${A}_{0}$ （mm²），在服役时间为t年时有效截面积的值为$A\left(t\right)$ （mm²）。

图  3  木柱轴心受压示意图

Figure  3.  Diagram of a wooden column under axial compression

下载: 全尺寸图片 幻灯片

根据上文所构建腐朽影响木构件剩余有效横截面积的时变模型（式(2)）和虫蛀、干缩裂缝影响木构件有效截面积损伤变量的时变模型（式(3) ~ (6)），叠加3种类型的损伤，得到木柱有效截面积A随时间变化的模型。

则木柱有效截面积随时间的变化率$\varphi \left(t\right)$为

由以上推导得到木柱在竖向荷载N作用下截面正应力$\sigma \left(t\right)$。

式中：${\sigma }_{0}$为构件的初始应力（MPa）。

木柱服役t年时应力比$\tau \left(t\right)$为

式中：${\tau }_{0}$为构件的初始应力比。

将式（11）代入式（7）得木柱在长期荷载、腐朽、虫蛀和干缩裂缝多因素影响下累积损伤变量的时变模型。

代入式（2）得到木柱相应的强度退化模型。

式（13）所示模型分别以初始应力比、腐朽层厚度和截面边长、虫蛀3个阶段的时间、裂缝相对深度作为量化长期荷载、腐朽、虫蛀和干缩裂缝损伤的参数，通过建立各参数与有效截面积的关系，实现不同损伤相同形式的累积，进一步构建有效截面积与抗压强度的关系式，得到多因素共同作用下木柱抗压强度随时间衰减的预测模型。该式可用于计算受到长期荷载、腐朽、虫蛀和干缩裂缝损伤的木柱在服役时间为t年时抗压强度的残余强度比，模型后续将作为构件可靠度分析的基础。

2.   基于概率密度演化方法的残损木构件可靠度分析

由前文所述，需要对古建筑残损木构件进行可靠度分析。在可靠度分析过程中，随着影响参数的增加，计算量大、结果精度不足、概率分布拟合不准确等问题就会随之出现，因此本文采用取值量要求小且计算精度高的概率密度演化方法^[19-20]，对古建筑残损木构件的残余强度进行可靠度评估。

以木柱为例，其抗压强度的残余强度预测方程如式（13）。式中不确定性参数较多，会增加可靠度分析时的计算量，为了提高计算效率，需要对参数进行初步降维。灵敏度分析^[21]是量化各参数对模型影响程度的重要方法，通过非关键参数常数化可以在不影响计算精度的前提下有效减少可靠度分析的计算量。因此进行可靠度分析之前，先进行影响参数的灵敏度计算，对灵敏度排序后确定影响构件残损的关键参数，关键参数在后续可靠度分析中仍为随机变量，则$\theta =\left({\theta }_{1},{\theta }_{2},\cdot \cdot \cdot ,{\theta }_{S}\right)$多维随机向量，前文已假设各影响因素之间相互独立，设其联合概率密度函数为${p}_{\theta }\left(\theta \right)$。

基于随机变量的分布类型和初始条件，采样得到随机变量$\theta$的代表点。带入目标函数式（13），利用差分法计算残余强度比变化率。由式（13）知，上述方程包括了影响木构件残损的主要随机变量，可看作一个随机系统，满足概率密度守恒原理，则可构建如下所示概率密度演化方程。

式中：x为残余强度比，$\dot{X}\left(\theta ,t\right)$为残余强度比变化率，${p}_{x\theta }\left(x,\theta ,t\right)$为残余强度比与随机参数的联合概率密度函数。

对应的初始条件为

式中：$\delta \left(\cdot \right)$为狄拉克函数，${t}_{0}$和${x}_{0}$分别为初始时间和初始强度比。

基于差分法对概率密度演化方程式（14）进行求解，得到残余强度比与随机参数的联合概率密度函数${p}_{x\theta }\left(x,\theta ,t\right)$。

在求解时，引入吸收边界条件

将解得的${p}_{x\theta }\left(x,\theta ,t\right)$在随机参数域内积分得残余强度比随时间变化的概率密度函数。

概率密度函数在规定的安全域内积分，得到可靠度

式中：${x}_{\mathrm{B}}$为构件失效界限值。

一般情况下，当构件的强度衰减至极限时，

或变形达到极限时，

判定构件失效。而当木构件因过度变形而失效时，其残余强度通常不低于构件在运行荷载作用下的应力比^[18]，因此将式（20），即$\alpha \left(t\right)\geqslant 1$作为失效界限，超过界限值1则失效。在后文中将应用本节所述方法对所构建多因子强度退化模型进行可靠度分析，并对比蒙特卡罗法的分析结果，以验证该方法的可行性。

3.   古建筑木结构典型柱可靠度评价

为验证前文所提出模型和方法的可行性，以文献[12]中从一座藏式古建筑木结构中维修替换下来的一组木柱为例，进行柱抗压强度退化的可靠度分析。共有32根材质相同的木柱，服役历史约为350 a，经测量后可知：方柱截面边长均值约为400 mm，腐朽层厚度均值约为4 mm，未见明显虫蛀损伤。

为实现本例中木柱强度退化可靠度的计算以及分析方法的验证，具体的步骤如下：

（1）根据所构建强度退化模型确定各随机变量，明确随机变量统计特征及其他参数的取值；

（2）对随机变量进行灵敏度计算并排序后设定阈值，确定关键参数，仅保留关键参数的随机性；

（3）基于概率密度演化方法对木柱的残余强度进行可靠度分析，并通过对比蒙特卡罗法的计算结果，验证概率密度演化方法的有效性。

3.1   随机变量的确定

基于前文构建残损木构件强度退化模型，定义构件初始应力比、腐朽层厚度、截面直径（或截面边长）、干缩裂缝宽度与深度的比值、白蚁在建筑50 m范围内建立一个成熟蚁巢所用时间（阶段1）、白蚁形成虫道并到达建筑的时间（阶段2）和白蚁突破或绕开所设置防蚁屏障的时间（阶段3）为量化损伤的随机变量。

针对本算例中随机变量及其他影响参数的情况做出以下说明与假设。

（1）本建筑采取物理防虫措施，外界环境可评定为中等；

（2）参考文献[12]和[18]，将本模型中的相关参数设置为${{{B}}}_{1}=13.952\;6$，${{{B}}}_{2}=25.049\;8$，${{k}}=-2.74\times 1{0}^{-5}$，${{b}}=0.000\;2$；此外，根据算例实际情况对C_MC取值，设定简化系数$E=1.82$，$F=-4.134$；

（3）根据以上基本情况可设定对木柱进行灵敏度分析时基准状态下各参数取值：初始应力比${\tau }_{0}=0.3$，腐朽层厚度${\delta }_{0}=4$ mm，截面边长$d=400$ mm，裂缝宽度和深度的比值$\lambda =0.2$，虫蛀初始期阶段1时间${t}_{1}=16.7$ a，虫蛀初始期阶段2时间${t}_{2}=4.9$ a，虫蛀初始期阶段3时间${t}_{3}=37.9$ a；

（4）根据相关研究^[22-23]可知，古建筑木构件尺寸、材料强度、构件应力比、腐朽层厚度和干缩开裂深度均为服从正态分布的随机变量，同时假设虫蛀模型中${t}_{1}$、${t}_{2}$和${t}_{3}$也服从正态分布。

基于以上假设，结合木柱材性试验^[12]及对虫蛀过程进行的调查分析和数据统计^[24]，得到各参数统计特征（表1）。

表  1  参数统计特征

Table  1.  Parameter statistical characteristics

参数 Parameter	均值 Mean value	变异系数 Variation coefficient
初始应力比 Initial stress ratio $({\tau }_{0})$	0.3	0.1
腐朽层厚度 Thickness of decay layer $({\delta }_{0})$/mm	4	0.1
截面边长 Section side length (d)/mm	400	0.05
虫蛀初始期阶段1/a Initial stage of insect infestation 1 (${t}_{1}$)/year	16.7	0.646
虫蛀初始期阶段2/a Initial stage of insect infestation 2 (${t}_{2}$)/year	4.9	0.813
虫蛀初始期阶段3/a Initial stage of insect infestation 3 (${t}_{3}$)/year	37.9	0.675
裂缝宽度和深度的比值 Ratio of crack width to depth ($\lambda$)	0.2	0.08

下载: 导出CSV 

| 显示表格

3.2   参数灵敏度分析

依前文述，通过灵敏度分析来判定各参数对残余强度影响的重要程度。本算例中残余强度由7个影响参数$x=\left\{{x}_{1},{x}_{2},\cdot \cdot \cdot ,{x}_{7}\right\}$决定，设${\varphi }_{\mathrm{M}}=f\left(x\right)$。对于某一设定基准状态${x}^{*}=\left\{{x}_{1}^{*},{x}_{2}^{*},\cdot \cdot \cdot ,{x}_{n}^{*}\right\}$，残余强度比为${\varphi }_{{\rm{M}}}^{*}$。在分析某一参数${x}_{i}$对残余强度比的影响时，其余参数的值保持不变，仅使${x}_{i}$在其参数域内发生变化。

根据文献[21]，将灵敏度函数${S}_{i}\left({x}_{i}\right)$定义为参数变化率、参数值与残余强度比的函数，即

取${x}_{i}={x}_{i}^{*}$，代入式（21）可得参数${x}_{i}$在基准状态${x}^{*}=\left\{{x}_{1}^{*},{x}_{2}^{*},\cdot \cdot \cdot ,{x}_{n}^{*}\right\}$下的灵敏度因子${S}_{i}^{*}$，

由此可计算得各参数灵敏度因子对应比重${\gamma }_{i}^{*}$，

以表1中各参数均值作为基准状态，分别代入式（22）和式（23）对各参数进行灵敏度计算，得到灵敏度因子${S}^{*}$和对应比重${\gamma }^{*}$如表2所示。

表  2  参数灵敏度因子（${S}^{*}$）及其对应比重（${\gamma }^{*}$）

Table  2.  Parameter sensitivity factor (${S}^{*}$) and specific gravity (${\gamma }^{*}$)

参数 Parameter	${\tau }_{0}$	${\delta }_{0}$	d	${t}_{1}$	${t}_{2}$	${t}_{3}$	$\lambda$
${S}^{\mathrm{*}}$	4.339 2 × 10−4	0.505 4 × 10−4	0.505 4 × 10−4	0.043 1 × 10−4	0.012 6 × 10−4	0.097 7 × 10−4	0.017 2 × 10−4
${\gamma }^{\mathrm{*}}$	0.786 0	0.091 5	0.091 5	0.007 8	0.002 3	0.017 7	0.003 1
注：${\tau }_{0}$为木柱初始应力比，${\delta }_{0}$为腐朽层厚度（mm），d为截面边长（mm），${t}_{1}$为虫蛀初始期第1阶段的时间（a），${t}_{2}$为虫蛀初始期第2阶段的时间（a），${t}_{3}$为虫蛀初始期第3阶段的时间（a），$\lambda$为干缩裂缝宽度和深度的比值，下同。Notes: ${\tau }_{0}$ is the initial stress ratio of the wooden columns, ${\delta }_{0}$ is the thickness of decay layer (mm), d is the section side length (mm), ${t}_{1}$ is the time of the first stage of the initial period of insect infestation (year), ${t}_{2}$ is the time of the second stage of the initial period of insect infestation (year), ${t}_{3}$ is the third stage of the initial period of insect infestation (year), $\lambda$ is the ratio of crack width to depth. The same below.

下载: 导出CSV 

| 显示表格

将各参数灵敏度因子对应比重进行排序并绘制直方图，如图4所示。

图  4  参数灵敏度排序

Figure  4.  Parameter sensitivity ranking

下载: 全尺寸图片 幻灯片

灵敏度分析应在筛选出关键参数的前提下保证后续计算的准确性，并尽可能考虑多类损伤。综合以上计算结果，灵敏度排序前4的参数（${\tau }_{0}$、${\delta }_{0}$、d、${t}_{3}$）灵敏度比重相加大于95%，并包含3种主要损伤，满足灵敏度分析的要求，因此本文以${t}_{3}$的灵敏度为界限，设定灵敏度因子${S}^{*}=0.099\;7\times 1{0}^{-4}$为阈值，则大于等于阈值的参数（${\tau }_{0}$、${\delta }_{0}$、d、${t}_{3}$）为关键参数，在后续计算中仍视为随机变量，而低于阈值的其余参数将不再考虑其随机性，视为常量处理，以此实现参数的降维。

3.3   构件可靠度分析

以上述对参数进行的灵敏度计算为基础，基于概率密度演化方法对本例木柱残余强度的可靠度进行如下分析。

（1）已知各变量的分布类型和统计信息如表1，利用合适的抽样方法在随机变量的参数空间${\varOmega }_{\theta }$内选择${N}_{\mathrm{s}\mathrm{e}\mathrm{l}}$个代表点集${\theta }_{q}\left(q=\mathrm{1,2},\cdots ,{N}_{\mathrm{s}\mathrm{e}\mathrm{l}}\right)$，本文取${N}_{\mathrm{s}\mathrm{e}\mathrm{l}}= 1 \;000$，并求得代表点相应的赋得概率${P}_{q}$。在众多选点方法中，拉丁超立方抽样法^[25]是可靠度分析常用的抽样方法之一，在抽样时能避免重复抽样，并以较小的样本量反映总体的变异规律，有较高的抽样效率，同时与针对概率密度演化方法提出的映射降维法^[26]、数论选点法^[27]和切球选点法^[28]相比更简单易懂，因此选择该方法对参数进行选点。

（2）将每一个代表点${\theta }_{q}$带入构件强度退化模型即式（13）进行计算，得到${N}_{\mathrm{s}\mathrm{e}\mathrm{l}}$组不同数据下的木柱在1 ~ 1 000 a的损伤变量及残余强度比随时间变化的趋势，如图5和图6，差分法计算得到对应残余强度比变化率$\dot{X}\left(\theta ,t\right)$。

图  5  损伤变量随时间变化规律

损伤变量 ≥ 1时，构件失效。The component failures when the damage variable is greater than or equal to 1.

Figure  5.  Variation pattern of damage variable with time

下载: 全尺寸图片 幻灯片

图  6  残余强度比随时间变化规律

Figure  6.  Variation law of residual strength ratio with time

下载: 全尺寸图片 幻灯片

图5和图6中曲线分别代表1 000组不同参数取值的木柱在服役时间1 ~ 1 000 a时，损伤变量和残余强度比随时间的变化。两图表明，随着服役时间增加，损伤量逐渐增大，残余强度比减小，损伤变量增长速率和残余强度衰减速率均逐渐增大。根据失效界限的定义，当构件的变形或强度衰减达到极限时判定为失效，反映在图中即损伤变量大于1，残余强度比小于0时构件失效。对比两图，在服役时间为1 000 a时，模拟得到的1 000组构件基于损伤变量来判断几乎都已达到失效界限值，而基于残余强度比来衡量则均未达到失效界限值，进一步说明了当木构件因变形过大发生失效的时间要早于残余强度达到极限时的时间。

此外，由文献[11]已知算例中木柱根据材性试验得到抗压强度的残余强度比约为55.65%，将各参数数据代入所构建强度退化模型进行计算，木柱服役500 a时的预测残余强度比为73.40%，两者相比较模型预测结果明显高于材性试验结果。说明本文构建模型存在一定局限性，还需要基于更多的新旧材试验数据对模型进行进一步的修正。

（3）根据概率守恒原理，将（2）求得的变化率代入概率密度演化方程式（14），并基于式（15）和式（16），通过差分法解得联合概率密度函数${p}_{x\theta }\left(x,\theta ,t\right)$。

（4）根据式（17），在随机参数域${\varOmega }_{\theta }$内积分，得到残损木柱残余强度比随时间变化的概率密度函数${p}_{x}\left(x,t\right)$，以概率密度为z轴的概率密度演化曲面如图7所示。

图  7  残余强度比的概率密度演化曲面

Figure  7.  Probability density evolution curved surface of residual strength ratio

下载: 全尺寸图片 幻灯片

通过上述概率密度演化的计算，建立起了如图7所示时间、残余强度比和概率密度的三维图像。图7表明，不同服役时间，发生概率最高的残余强度比不同，随着服役时间增加，发生概率最高的残余强度比逐渐减小。同时，由图7所示的概率密度演化曲面可以预测任意时间时任意残余强度比的发生概率，例如曲面上的点（0.11，548，0.61），表明当木柱服役时间为548 a时，其抗压强度的残余强度比为0.11的概率密度是0.61。

（5）由式（18）将概率密度函数${p}_{x}\left(x,t\right)$在安全域[0,1]内积分，得到相应的残余强度可靠度$R\left(t\right)$并求出失效概率$F\left(t\right)$，失效概率随时间变化的趋势如图8所示。

图  8  失效概率随时间变化规律

Figure  8.  Changing law of failure probability with time

下载: 全尺寸图片 幻灯片

图8表明：随着服役时间的增加，失效概率逐渐增高，则其可靠度就越低。此外，由图8所示模型可知1 ~ 1 000 a内任意时间木柱的失效概率，如图中点（700，0.885），表明当服役时间为700 a时，木柱的失效概率为0.885。

上文基于所构建多因子强度退化模型，利用概率密度演化理论对算例中的残损木柱进行可靠度分析，分别得到了服役时间为1 ~ 1 000 a时损伤变量、残余强度比、失效概率的变化趋势，以及残余强度比的概率密度函数。为了验证概率密度演化方法的分析结果，使用蒙特卡洛法在相同参数条件下对本例残损木柱进行可靠度分析。基于蒙特卡罗法随机抽样产生$n=1{0}^{4}$组参数数据，以式（13）为基础计算得到1 ~ 1 000 a的失效概率，两种方法在$t=100$、$t=300$、$t=500$、$t=700$、$t=900$时的失效概率计算结果对比如表3所示。

表  3  概率密度演化方法（PDEM）与蒙特卡罗法（MC）失效概率计算结果对比

Table  3.  Comparison of failure probability calculation results between probability density evolution method (PDEM) and Monte Carlo method (MC)

时间/a Time/year	失效概率 Failure probability/%		差值 Difference value/%
	PDEM	MC
100	2.90	12.38	9.48
300	38.50	34.58	3.92
500	67.20	61.10	6.10
700	91.60	100.00	8.40
900	95.60	100.00	4.40

下载: 导出CSV 

| 显示表格

表3所示，随着时间变化，PDEM和MC两种方法计算所得失效概率逐渐增加，在几个服役时间节点的失效概率差值均小于10%，证明概率密度演化理论适用于古建筑木结构残损木构件残余强度的可靠度分析。同时，表4表明：在参数条件相同的前提下，PDEM参数取样量较MC参数取样量更少，计算时间更短，但得到了差值较小的计算结果，进一步证明了PDEM的高效性。

表  4  PDEM与MC计算效率对比

Table  4.  Comparison of PDEM and MC calculation efficiency

研究方法 Research method	抽样次数 Number of sampling	运算时间 Operation time/s
PDEM	1 000	1.28
MC	10 000	25.83

下载: 导出CSV 

| 显示表格

4.   结　　论

本研究考虑长期荷载、腐朽、虫蛀和干缩裂缝4个重要影响因素，建立了古建筑残损木构件的强度退化时变模型。以一古建筑木结构中的柱为例，利用灵敏度方法实现了参数的降维，基于概率密度演化方法进行了构件残余强度的可靠度分析，并与蒙特卡罗法的计算结果进行了对比，主要结论如下。

（1）建立了古建筑木柱的多因子强度退化模型，可用于其残余强度的评估和预测；

（2）所建立的基于PDEM的古建筑残损木构件可靠度分析方法具有较好的适用性，且理论基础明确，计算效率高。