CN110489838B - 一种基于贝叶斯推理的湿木材热参数反演方法 - Google Patents

Abstract

本发明公开了一种基于贝叶斯推理的湿木材热参数反演方法，首先，基于阿伦尼乌斯定律建立湿木材热分解本构模型,详细考虑了木材热解过程中的多个阶段；其次，基于贝叶斯方法建立木材热参数反演模型；再次，基于马尔科夫链蒙特卡罗方法(MCMC)给出所建立模型的求解算法，提出混合反演算法；最后为了更好的实际应用，还提出一种基于MATLAB、ABAQUS和Python的求解算法通用计算方案，该计算方案灵活调用多种计算软件，实现了复杂运算和数据的自动化处理。本发明与传统方法相比，所提方法与实际热解过程拟合较好，同时简单高效，具有较高的热参数反演精度，具备巨大的工程应用价值。

Description

一种基于贝叶斯推理的湿木材热参数反演方法

技术领域

本发明属于土木工程、能源与环境技术领域，特别涉及一种基于贝叶斯推理的湿木材热参数反演方法。

背景技术

木材料在国内外现代建筑和古建筑中均得到了广泛的应用，木结构的安全性和耐久性需要着重考虑，关乎到人民的生命财产安全和历史文化价值。由于木材是一种可燃材料，火灾是影响木建筑安全性的重要的因素，木结构的防火问题和正确的防火设计是需要重点研究的关键问题。预测暴露在火灾下的木结构的结构行为和力学性能的退化需要得到更多的关注，这涉及到对影响木结构在火灾中性能的物理现象的正确认识。现如今木结构的材性建模和防火设计主要参考欧洲规范(Eurocode5)，但是该建模方法过于简单，有较多简化，无法准确表征木材在加热分解过程中的复杂物理化学变化。同时，目前ABAQUS和ANSYS等主流通用有限元软件中缺少现有的本构模型，可以对木材等可燃材料进行较好的热分析。为了更好的研究火灾下木结构的结构行为，需要更加准确的木材热分解建模方法。

为了更好的进行木结构的防火设计，不仅需要准确的木材热分解建模，更为关键的是在正确建模的基础上，探究木结构在火灾中的热参数的变化特征。当前木结构的防火的主要措施是在木结构表面涂防火涂料，这种方式可以延缓火灾中木材的升温。探究木结构在有防火涂料保护下热参数的变化特征，并以此来指导防火设计，也是值得关注的重要问题。这类问题往往非常复杂，是一种典型的热传导反问题。人们也提出了诸如采用粒子群算法、遗传算法、伴随方程法、共轭梯度法和一些非迭代方法进行反问题的研究，但是在参数反演的精度和效率方面还存在诸多问题，许多方法存在效率过低，反演精度不高等问题。所以迫切需要一种兼具精度和效率的火灾中，复杂环境下木材热参数反问题求解方法。

发明内容

发明目的：为了克服现有技术中存在的不足，提供一种基于贝叶斯推理的湿木材热参数反演方法，可以有效描述木材热解过程的多阶段反应模型，并且可以有效反演木材热分解反应中的相关热参数，有较高的反演精度和计算效率。

技术方案：为实现上述目的，本发明提供一种基于贝叶斯推理的湿木材热参数反演方法，包括如下步骤：

S1：基于阿伦尼乌斯定律建立湿木材热分解本构模型；

S2：基于贝叶斯方法建立木材热参数反演模型；

S3：基于马尔科夫链蒙特卡罗方法给出步骤S2所建立反演模型的求解算法，提出混合反演算法。

进一步的，所述步骤S1中湿木材热分解本构模型的建立具体为：

由于湿木材在有氧环境下通过导热和对流换热的热分解过程分为三个阶段，第一阶段是湿木材的干燥和水分的挥发，第二阶段是成炭过程，第三阶段是木炭的缓慢氧化放热最终变成灰烬，具体的反应如下，每个反应都由一阶阿伦尼乌斯定律模拟，假定T1为湿木材，T2为干木材，W为湿木材内部水分，V为水蒸气，C为木炭，G为挥发物，A为灰烬，则多阶段热解反应可表示为：

W→V

T2→C+G1

C→A+G2

湿木材热解的守恒方程由下式表示：

其中，T表示温度，λ_x、λ_y和λ_z表示x、y和z三个方向随温度变化的导热系数，ρ_w、ρ_c和ρ_l表示木材、木炭和液体的密度，c_w、c_c和c_l表示木材、木炭和液体的比热容，t是时间变量，Q″为在温度T下三种热解反应的反应热之和，其与每个热解反应的反应速率、活化能和热反应有关，下式给出了反应热Q″的求解表达式：

木材、木炭、挥发性气体、湿木材内部水分和水蒸汽的质量守恒方程表示为下式：

进一步的，所述步骤S1中木材、木炭、挥发性气体、湿木材内部水分和水蒸汽的质量守恒方程式中k用阿伦尼乌斯公式计算，其具体表达式如下

其中，k为速率常数，R为摩尔气体常量，T为热力学温度，E_a为表观活化能，A为指前因子，也成为频率因子。

进一步的，所述步骤S1中k的求解需要确定初始温度分布和边界条件，其具体的确定方式如下：

初始条件：t＝0时，T(x,y,z,0)＝T₀(x,y,z)，T₀表示初始温度场；

边界条件：热流通过对流和辐射与木材的外露表面交换热量，由以下边界条件表示：

其中，n表示木材表面的外法线方向，σ_S-B是Stefan-Boltzmann常数，ε是辐射系数，h_c是对流系数，T_f表示炉内测量温度，表征外加热流。

进一步的，所述步骤S2中木材热参数反演模型的建立具体为：

S2-1：设随机事件组A₁,A₂,...,A_n,...是样本空间Ω的一组可列无穷划分，p(A_i)＞0,i＝1,2,...，

则

称为贝叶斯公式，上式是离散分布的情形，对于连续分布的情形，贝叶斯公式表达为下式：

S2-2：上两式是贝叶斯公式的一般形式，将贝叶斯公式改写作下式：

X表示未知随机变量Y表示与该随机变量相对应的相关“观测”值数据，p_X|Y(x|y)为后验概率密度函数，p_Y|X(y|x)为似然函数，边缘分布密度p_X(x)为先验概率密度函数；

在对后验状态空间进行数值求解时，不需要计算式中的归一化常数c＝p_Y(y)，其未知变量X的概率分布反演无影响，因此可以将后验概率密度函数表达式简写为：

p_X|Y(x|y)∝p_Y|X(y|x)p_X(x)

S2-3：建立马尔可夫链蒙特卡罗方法进行后验概率密度函数反演算法，在实际反演过程中，不同的先验分布最终均能够收敛到目标概率分布，区别在于算法效率有高低，因此通常情况下，先验分布有多种表达形式。均匀分布、高斯分布、马尔可夫随机场都被用于先验分布建模中。例如先验分布为高斯分布时，有

S2-4：似然函数反映了未知模型参数与实际观测数据之间的拟合程度，包括了所有观测数据下未知的随机变量参数的似然度信息，设定测量误差与观测数据之间相互独立，用ω_m表示测量噪声，则可以将统计模型表示为下式：

Y＝F(X)+ω_m

S2-5：将随机噪声ω_m建模为具有高斯分布的平稳白噪声，即平稳高斯白噪声，假定高斯白噪声的均值和方差已知，设均值为0，方差为v_T，因此可以将似然函数表示为

得到后验概率密度函数求解公式：

将常数项省略，得到简化后的后验概率密度求解公式：

S2-6：将上式模型修改为下式：

p(k|Y)＝0,k≤0

其中k表示木材的热参数，F(k)表示之前建立的正问题模型的求解算子，Y表示观测温度；

如果热参数的均值k和方差ν_k未知，则可以得到分级后验概率密度函数公式：

如果似然函数中，高斯白噪声的方差ν_T也未知，则可以得到增强分级后验概率密度函数公式：

进一步的，所述步骤S3中求解算法具体为：

S3-1：初始化马尔可夫链：

和

S3-2：循环采样，For t＝0:Nmcmc-1，Nmcmc表示马尔科夫链长度；

S3-3：高维Gibbs循环采样—反演高维未知参数θ^(t)(t≥0)：

采样：

采样：

采样：

采样：

S3-4：高维Metropolis-Hastings采样—反演高维未知参数X^(t)(t≥0)：

循环X的维度，For j＝1:n；

采样

从均匀分布中采样u_x～Uniform(0,1)；

如果

接受转移

否则拒绝转移：

其中

S3-5：一维Metropolis-Hastings采样—高斯白噪声均方差ν_T：

采样

从均匀分布中采样

如果

接受转移

否则拒绝转移：

其中

S3-6：

即为满足需要的样本集。

进一步的，所述步骤S3中基于ABAQUS、MATLAB和Python对模型求解算法进行通用计算，其方法具体为：

MATLAB编写MCMC算法，需求信息写入TXT文件，Python读取TXT文件，修改INP文件，并调用ABAQUS，计算结果反馈回DAT文件，Python读取DAT文件，处理输出数据写入TXT文件中，MATLAB读取计算结果并进行下一步计算。

有益效果：本发明与现有技术相比，与实际热解过程拟合较好，同时简单高效，并且能够兼具非常高的热参数反演精度和计算效率，具备巨大的工程应用价值。

附图说明

图1是本发明提出的反演方法思路概况图；

图2是湿木材加热分解过程及产物示意图；

图3是本发明提出的计算方案流程图；

图4是具体的计算思路示意图；

图5是快速计算方法流程图；

图6是增量步自适应算法流程图；

图7是基于ABAQUS和UMATHT用户子程序的计算流程图；

图8是实施例1中的一维热传导模型示意图；

图9是实施例1中的马尔可夫链采样示意图；

图10是实施例1中的反演温度和真实温度比较示意图；

图11是实施例2中的马尔可夫链采样示意图；

图12是实施例2中的目标温度与反演温度对比图。

具体实施方式

下面结合附图和具体实施例，进一步阐明本发明。

如图1所示，本发明提供一种基于贝叶斯推理的湿木材热参数反演方法，该方法主要有四部分构成，依次分别是：一种木材热分解本构模型、基于本构模型的热参数反演模型、一种用于求解模型的模型求解算法、一种适合求解算法的通用计算方案。本实施例中分别就这四个部分进行详细说明。

一、木材热分解模型(正问题)

首先建立木材热分解模型，研究热量、种类和水分的变化的影响。目前，ABAQUS和ANSYS等通用有限元软件中缺少现有的本构模型，可以对木材等可燃材料进行热分析。通常人们使用欧洲规范(Eurocode5)进行木材材性建模，但是该建模方法过于简单，有较多简化，无法准确表征木材在加热分解过程中的物理化学变化。

这一部分本发明采用一种基于阿伦尼乌斯公式的木材热分解模型建模方法，与传统方案相比，所提方案可以自定义材料、可以包含更多的物理化学过程，因此更为贴近实际，有非常高的应用价值。

如图2所示，给出了湿木材在有氧环境下加热分解的过程，这是一个复杂的物理化学反应过程。

由图2可知，湿木材在有氧环境下的通过导热和对流换热的热分解过程大致分为三个阶段，第一阶段是湿木材的干燥，水分的挥发(包括自由水和结合水)，第二阶段是成炭过程，第三阶段是木炭的缓慢氧化放热最终变成灰烬。具体的反应如下，每个反应都由一阶阿伦尼乌斯定律模拟。假定T1为湿木材，T2为干木材，W为湿木材内部水分，V为水蒸气，C为木炭，G为挥发物，A为灰烬，则多阶段热解反应可表示为：

W→V

T2→C+G1

C→A+G2

湿木材热解的守恒方程可以由下式表示：

其中，T表示温度，λ_x、λ_y和λ_z表示x、y和z三个方向随温度变化的导热系数，ρ_w、ρ_c和ρ_l表示木材、木炭和液体的密度，c_w、c_c和c_l表示木材、木炭和液体的比热容，t是时间变量。Q″为在温度T下三种热解反应的反应热之和，与每个热解反应的反应速率、活化能和热反应有关，下式给出了反应热Q″的求解表达式。

木材(timber)、木炭(char)、挥发性气体(gas)、湿木材内部水分(water)和水蒸汽(vapor)的质量守恒方程可以表示为下式：

因为湿木材热解的化学反应可以用一阶阿伦尼乌斯公式进行描述。因此上式中，k可以用阿伦尼乌斯公式计算。具体表达式如下

其中，k为速率常数，R为摩尔气体常量，T为热力学温度，E_a为表观活化能，A为指前因子(也成为频率因子)。

上述所提的瞬态热传导控制方程的求解需要确定初始温度分布和一定的边界条件，其具体方案如下：

(1)初始条件：t＝0时，T(x,y,z,0)＝T₀(x,y,z)，T₀表示初始温度场。

(2)热流通过对流和辐射的与木材的外露表面交换热量，可以由以下边界条件表示：

其中，n表示木材表面的外法线方向；σ_S-B是Stefan-Boltzmann常数；ε是辐射系数；h_c是对流系数；T_f表示炉内测量温度(f:furnace)，表征外加热流。

二、基于贝叶斯方法的反演模型

导热系数求解本质为反问题，此部分提出了一种基于贝叶斯方法的反问题求解模型。贝叶斯估计的主要目的是根据当前获得的信息或先前的知识中当中推导出未知量的概率分布，贝叶斯方法的基础是贝叶斯公式。设随机事件组A₁,A₂,...,A_n,...是样本空间Ω的一组可列无穷划分，p(A_i)＞0,i＝1,2,...。

则

称为贝叶斯公式。上式是离散分布的情形，对于连续分布的情形，贝叶斯公式可表达为下式：

上两式是贝叶斯公式的一般形式，本专利中，将贝叶斯公式改写作下式：

其中，X表示未知随机变量Y表示与该随机变量相对应的相关“观测”值数据。p_X|Y(x|y)被称为后验概率密度函数，贝叶斯方法的主要目标就是根据事情的果(即已有的观测数据)去推导事情的因，即计算出未知随机变量的概率分布也就是后验概率密度函数(posterior probability density function,PPDF)，p_Y|X(y|x)为似然函数，边缘分布密度p_X(x)为先验概率密度函数。基于上式可看出，给定一定观测值(证据)的后验概率与其似然概率和先验(无条件)概率的乘积成正比。贝叶斯推理为估计、假设检验、模型选择和决策问题提供了强有力的技术支撑。在对后验状态空间进行数值求解时，不需要计算式中的归一化常数c＝p_Y(y)，其未知变量X的概率分布反演无影响，因此可以将后验概率密度函数表达式简写为：

p_X|Y(x|y)∝p_Y|X(y|x)p_X(x)

本发明中建立了马尔可夫链蒙特卡罗方法进行后验概率密度函数反演算法，在实际反演过程中，不同的先验分布最终均能够收敛到目标概率分布，区别在于算法效率有高低，因此通常情况下，先验分布有多种表达形式。均匀分布、高斯分布、马尔可夫随机场都被用于先验分布建模中。例如先验分布为高斯分布时，有

似然函数反映了未知模型参数与实际观测数据之间的拟合程度，包括了所有观测数据下未知的随机变量参数的似然度信息。通常观测的数据在测量仪器的不确定度(传感器的噪声)、材料自身因素和测量方法的局限性会产生测量误差即观测噪声。可以假定测量误差与观测数据之间相互独立，用ω_m表示测量噪声，则可以将统计模型表示为下式：

Y＝F(X)+ω_m

本发明中将随机噪声ω_m建模为具有高斯分布的平稳白噪声，即平稳高斯白噪声。假定高斯白噪声的均值和方差已知，设均值为0，方差为v_T，因此可以将似然函数表示为

因此，可以得到后验概率密度函数求解公式：

将常数项省略，可以得到简化后的后验概率密度求解公式：

因为本发明提出的是一种木材热参数的反演方法，因此上式模型可以修改为下式：

p(k|Y)＝0,k≤0

其中k表示木材的热参数，F(k)表示之前建立的正问题模型的求解算子，Y表示观测温度。

如果热参数的均值

和方差ν_k未知，则可以得到分级后验概率密度函数(hierarchical PPDF)公式：

如果似然函数中，高斯白噪声的方差ν_T也未知，则可以得到增强分级后验概率密度函数(augmented and hierarchical PPDF)公式：

不同的先验分布和似然函数构成的后验概率密度函数在一定意义上都是正确的，但是不同的后验概率密度函数在收敛速度和精度有所不同，需要经过多次尝试选取最优的组合搭配形式，达到高效率的贝叶斯反演。

三、基于MCMC方法的模型求解算法

此部分提出一种基于MCMC方法的反问题模型求解算法，MCMC全称是马尔可夫链蒙特卡罗方法(Markov Chain Monte Carlo)，是一种随机模拟方法，其利用马尔可夫链可以收敛到平稳分布的性质，构造一个平稳分布是目标分布的马尔可夫链，收敛到平稳分布之后的马尔可夫链就是目标概率分布的样本。如果分布形式非常复杂，如高维分布、隐式的高维联合分布，MCMC方法依然适用。

MCMC方法中，广泛采用的是Metropolis-Hastings算法和Gibbs算法，Metropolis-Hastings算法往往针对低维问题，在高维问题中，接受率较低，往往采用Gibbs算法。Gibbs算法总是接受的。此部分给出一种Metropolis-Hastings算法+Gibbs算法的混合反演算法，可以解决多反演目标且反演目标维度不同的情况。对于单目标、或者单维度等更简单的问题，可以使用所提算法的一部分或者几部分。

假定待反演高维参数1有m维，设为

可以是m维的导热系数也可以是m维的未知边界热流密度，假定待反演高维参数2有n维，设为

可以是n维的导热系数也可以是n维的未知边界热流密度，显然与θ不能为同一未知参数。还可以带有未知参数ν_T，一般用来表征高斯白噪声的均方差，如果反演了马尔可夫随机场建模的热流密度也可以用来表征马尔可夫随机场的归一化参数。为了算法表述方便，将

简写为

同理。算法中，反演高维未知参数θ^(t)(t≥0)采用的是高维Gibbs算法，反演高维未知参数X^(t)(t≥0)采用的是高维Metropolis-Hastings算法，反演一维未知参数ν^T采用的是一维Metropolis-Hastings算法。在高维反演问题中，如果完整的条件分布为常规分布(如

)则建议使用Gibbs算法，反之，则只能使用高维Metropolis-Hastings算法。在一维反演问题中，则使用一维Metropolis-Hastings算法。

1)初始化马尔可夫链：

和

2)循环采样，For t＝0:Nmcmc-1，Nmcmc表示马尔科夫链长度；

3)高维Gibbs循环采样—反演高维未知参数θ^(t)(t≥0)：

采样：

采样：

采样：

采样：

4)高维Metropolis-Hastings采样—反演高维未知参数X^(t)(t≥0)：

循环X的维度，For j＝1:n；

采样

从均匀分布中采样u_x～Uniform(0,1)；

如果

接受转移

否则拒绝转移：

其中

5)一维Metropolis-Hastings采样—高斯白噪声均方差ν_T：

采样

从均匀分布中采样

如果

接受转移

否则拒绝转移：

其中

6)

即为满足需要的样本集。

上述算法中，

表示建议分布。通常建议分布可以取正态分布或均匀分布。本发明首次引入多种建议分布结合思想。先大范围采样，可以有效避免收敛于局部最优的情况，后进行小范围精细化采样，提高收敛精度。多种建议分布结合思想与单一建议分布相比，可以加速收敛至全局最优解，提高采样效果。一种可行的建议分布函数如下式所示：

木材热流密度或导热系数等相关参数反演问题中，后验概率密度函数往往包含了与高斯白噪声相关的似然函数。高斯白噪声的均值通常为0，方差为v_T。通常的方法是按照经验给定方差v_T为一个较小的常数，或者按照上述算法所提的，将其建模为未知变量进行反演。但是该方法涉及多次正问题计算，效率较低，而给定为固定常数的方法收敛速度较低。此部分针对这个问题，给出一种多阶段似然函数方法。在采样初期，σ_T较小

则后验概率密度较大，提高接受率，可以加速收敛；在采样进行一段时间后，更换似然函数，使得σ_T较大，则后验概率密度减小，降低接受率，可以提高采样精度。下式为一种三阶段似然函数改变方法，通过改变标准差σ_T改变似然函数数值，也可以通过扩大分子乘数的方法进行调整。

四、基于ABAQUS、MATLAB和Python的模型求解算法通用计算方案

为了更好的实际应用，此部分还提出一种基于ABAQUS、MATLAB和Python的模型求解算法的通用计算方案。

该计算方案基于ABAQUS、MATLAB和Python，灵活调用多种计算软件。其中，MATLAB软件进行贝叶斯方法架构，ABAQUS有限元软件进行木材热解正问题计算，Python用来控制MATLAB与ABAQUS之间的交互，从而实现复杂运算和数据的自动化处理。

如图3所示，其计算方法为：

MATLAB编写MCMC算法，需求信息写入TXT文件，Python读取TXT文件，修改INP文件，并调用ABAQUS，计算结果反馈回DAT文件，Python读取DAT文件，处理输出数据写入TXT文件中。MATLAB读取计算计算结果并进行下一步计算。

MATLAB编写的MCMC求解算法进行导热系数反演的流程，具体如图4所示。反演热流密度等其他参数与下类似。

由于导热系数范围通常较小，本发明提出一种快速计算方法，设Hash数组储存每次计算的结果，Hash数组下标表示导热系数，数值表示后验概率密度。在计算时，如果当前导热系数k已经被计算过(Hash[k]≠0)，则直接使用Hash[k]的数值，如果没有计算过(Hash[k]＝0)，才进行有限元计算。显然，此时决定计算用时的不再是马尔可夫链长度，而是导热系数k状态空间的数目。由于之前已经假设0＜k≤1，兼顾到精度和计算需求，将k状态划分为1000种，在区间[0,1]内，均匀划分1000段，每一段用一个概率密度来表示。因为一次有限元计算耗时大约9秒，这样，遍历整个状态空间的耗时为1000×9＝9000s＝2.5h，考虑到一次MCMC采样不一定可以遍历所有状态空间，忽略非ABAQUS计算耗时，采样结束总时间T≤2.5h，对于长度为10000的马尔可夫链，计算用时降低为原来的1/10。对于长度为100000的马尔可夫链，则计算用时降低为原来的1/100。精度可以达到0.001，效果非常显著。考虑到数组下表必须为正整数，实际编程中，设k＝fix(k×1000)，fix为MATLAB内置库函数，数学上向下取整操作。同时考虑到当后验概率密度很小时，Hash[k]＝0，但是k已经被计算过，设立Flag数组，其数值由下式确定。其具体的算法流程如图5所示。

由于有限元计算中存在同样的增量步设置下，不同的导热系数下收敛情况不同的现象，本发明还提出一种增量步自适应算法，用于提高计算收敛概率。在自适应算法下，对于不同的导热系数，基本上可以得到100％的收敛效果。基本思想是在无法计算收敛时，自动调整初始增量步和最大增量步，直至计算收敛。调整策略为，首先将初始增量步和最大增量步设为足够小，计算收敛则跳出循环，如果不收敛则每次成比例扩大最大增量步，直至到一个比较大的值，如果超过这个值会带来较大插值误差，如果依然不收敛则按比例初始增量步(乘2)，最大增量步和前类似，从较小的值按比例扩大，直至计算收敛。这种策略基本保证了对所有的k均能计算收敛。总的数据插值+增量步子适应处理流程具体如图6所示。

针对模型求解算法中正问题求解算子F(x)，本发明还提出一种基于ABAQUS通用有限元软件和UMATHT用户子程序的计算方案。首先在ABAQUS/CAE界面中拟定所研究问题的几何模型、分析步、边界条件、单元与网格；基于Fortran语言编写湿木材热分解本构模型的UMATHT用户子程序，并将其引入ABAQUS有限元程序中，将描述模型的偏微分方程离散到网格层上，形成一个强耦合的代数方程组。基于Newton-Raphson算法，在每次增量结束时，使用代数方程组迭代计算未知变量；有限元模型建立完成，生成计算用的INP文件；调用ABAQUS运行计算INP文件；保存生成的DAT文件的计算结果，以备后续使用。因为在INP文件中修改模型十分方便，同时基于INP文件调用ABAQUS进行有限元计算效率高，不需要打开CAE界面即可完成，满足程序自动化调用、自动数据处理和多次有限元计算的需求。ABAQUS热分析正问题求解算子的计算流程具体如图7所示。

本发明中分别利用上述技术方案分别进行导热系数反演和热流密度反演，其具体如下：

实施例1：

导热系数反演

木材一维热传导问题中，导热系数反演。热电偶位置为距离自由端40、50、60、65和70mm，图8给出了一维热传导模型。

假定导热系数为一维参数，真值为0.5，采样收敛之后的后验均值估计

与真实值误差为0.1564％。图9为马尔可夫链采样图示，可以看到采样效果较好，在采样初期就第一次收敛，在1000次左右，收敛到真值。

如图10所示，为了直观看出反演出的导热系数反演效果，将k＝0.499218带入ABAQUS计算，绘制出反演温度和真实温度的97％置信区间，可以直观看出良好的反演效果。

实施例2：

热流密度反演

反演3维的热流密度，热流密度真值为50、100和200W·m²，图11给出了不同热流密度的采样图示，马尔可夫链长度为1500，可以发现热流密度均有较好的反演效果，结合采样结果，最终收敛以后(选取为1450-1500次)热流密度参数的后验均值估计分别为，

和

3个反演结果与目标热流密度之间的相对误差分别为1.93％、1.22％和1.23％。误差相当小，并且采样次数只有1500次，算法效率非常高。

图12给出了反演温度与目标温度的对比图，通过图12可以看到两者几乎没有差别，侧面反映了所提木材热参数反演算法具备极高的精度。

Claims (2)

1.一种基于贝叶斯推理的湿木材热参数反演方法，其特征在于：包括如下步骤：

S1：基于阿伦尼乌斯定律建立湿木材热分解本构模型；

S2：基于贝叶斯方法建立木材热参数反演模型；

S3：基于马尔科夫链蒙特卡罗方法给出步骤S2所建立反演模型的求解算法，提出混合反演算法；

所述步骤S1中湿木材热分解本构模型的建立具体为：

由于湿木材在有氧环境下通过导热和对流换热的热分解过程分为三个阶段，第一阶段是湿木材的干燥和水分的挥发，第二阶段是成炭过程，第三阶段是木炭的缓慢氧化放热最终变成灰烬，具体的反应如下，每个反应都由一阶阿伦尼乌斯定律模拟，假定T1为湿木材，T2为干木材，W为湿木材内部水分，V为水蒸气，C为木炭，G为挥发物，A为灰烬，则多阶段热解反应可表示为：

W→V

T2→C+G1

C→A+G2

湿木材热解的守恒方程由下式表示：

其中，T表示温度，λ_x、λ_y和λ_z表示x、y和z三个方向随温度变化的导热系数，ρ_w、ρ_c和ρ_l表示木材、木炭和液体的密度，c_w、c_c和c_l表示木材、木炭和液体的比热容，t是时间变量，Q″为在温度T下三种热解反应的反应热之和，其与每个热解反应的反应速率、活化能和热反应有关，下式给出了反应热Q″的求解表达式：

木材、木炭、挥发性气体、湿木材内部水分和水蒸汽的质量守恒方程表示为下式：

所述步骤S1中木材、木炭、挥发性气体、湿木材内部水分和水蒸汽的质量守恒方程式中k用阿伦尼乌斯公式计算，其具体表达式如下

其中，k为速率常数，R为摩尔气体常量，T为热力学温度，E_a为表观活化能，A为指前因子，也称为频率因子；

所述步骤S1中k的求解需要确定初始温度分布和边界条件，其具体的确定方式如下：

初始条件：t＝0时，T(x,y,z,0)＝T₀(x,y,z)，T₀表示初始温度场；

边界条件：热流通过对流和辐射与木材的外露表面交换热量，由以下边界条件表示：

其中，n表示木材表面的外法线方向，σ_S-B是Stefan-Boltzmann常数，ε是辐射系数，h_c是对流系数，T_f表示炉内测量温度，表征外加热流；

所述步骤S2中木材热参数反演模型的建立具体为：

S2-1：设随机事件组A₁,A₂,...,A_n,...是样本空间Ω的一组可列无穷划分，p(A_i)＞0,i＝1,2,...，

则

称为贝叶斯公式，上式是离散分布的情形，对于连续分布的情形，贝叶斯公式表达为下式：

S2-2：上两式是贝叶斯公式的一般形式，将贝叶斯公式改写作下式：

X表示未知随机变量Y表示与该随机变量相对应的相关“观测”值数据，p_X|Y(x|y)为后验概率密度函数，p_Y|X(y|x)为似然函数，边缘分布密度p_X(x)为先验概率密度函数；

在对后验状态空间进行数值求解时，不需要计算式中的归一化常数c＝p_Y(y)，其未知变量X的概率分布反演无影响，因此将后验概率密度函数表达式简写为：

p_XY(x|y)∝p_Y|X(y|x)p_X(x)

S2-3：建立马尔可夫链蒙特卡罗方法进行后验概率密度函数反演算法；

S2-4：似然函数反映了未知模型参数与实际观测数据之间的拟合程度，包括了所有观测数据下未知的随机变量参数的似然度信息，设定测量误差与观测数据之间相互独立，用ω_m表示测量噪声，则将统计模型表示为下式：

Y＝F(X)+ω_m

S2-5：将随机噪声ω_m建模为具有高斯分布的平稳白噪声，即平稳高斯白噪声，假定高斯白噪声的均值和方差已知，设均值为0，方差为v_T，因此可以将似然函数表示为

得到后验概率密度函数求解公式：

将常数项省略，得到简化后的后验概率密度求解公式：

S2-6：将上式模型修改为下式：

p(k|Y)＝0,k≤0

其中k表示木材的热参数，F(k)表示之前建立的正问题模型的求解算子，Y表示观测温度；

如果热参数的均值

和方差ν_k未知，则得到分级后验概率密度函数公式：

如果似然函数中，高斯白噪声的方差ν_T也未知，则得到增强分级后验概率密度函数公式：

所述步骤S3中求解算法具体为：

S3-1：初始化马尔可夫链：

i＝1,2,...,n和

S3-2：循环采样，For t＝0:Nmcmc-1，Nmcmc表示马尔科夫链长度；

S3-3：高维Gibbs循环采样—反演高维未知参数θ^(t)(t≥0)：

S3-4：高维Metropolis-Hastings采样—反演高维未知参数X^(t)(t≥0)：

循环X的维度，For j＝1:n；

采样

从均匀分布中采样u_x～Uniform(0,1)；

如果

接受转移

否则拒绝转移：

其中

S3-5：一维Metropolis-Hastings采样—高斯白噪声均方差ν_T：

采样

从均匀分布中采样

如果

接受转移

否则拒绝转移：

其中

S3-6：

即为满足需要的样本集。

2.根据权利要求1所述的一种基于贝叶斯推理的湿木材热参数反演方法，其特征在于：所述步骤S3中基于ABAQUS、MATLAB和Python对模型求解算法进行通用计算，其方法具体为：

MATLAB编写MCMC算法，需求信息写入TXT文件，Python读取TXT文件，修改INP文件，并调用ABAQUS，计算结果反馈回DAT文件，Python读取DAT文件，处理输出数据写入TXT文件中，MATLAB读取计算结果并进行下一步计算。

Images