随着计算机技术的发展，数值模拟已经成为复杂系统研究及设计的重要手段。以往，研究人员常常将材料参数、几何尺寸等作为确定性的输入参数，代入程序进行数值模拟，从而得到确定性的结果。然而实际上，数值模拟过程中存在各种不确定性来源。近十几年，不确定度量化(uncertainty quantification, UQ)方法发展迅速，逐渐成为数值模拟研究的热点^[1-2]。不确定度量化主要关注如何表征、量化、减小数值模拟过程中的各种不确定性来源，是实现高置信度数值模拟的有效手段。在进行不确定度量化分析时，为了研究输入参数引入的不确定性，每个输入参数常常需要人为给定概率分布。输入参数概率分布的选择很可能会对分析结果有直接影响，为了减少分布选择的主观性，输入参数的概率分布最好能够由试验数据提供。

目前，圆筒试验是评估炸药的做功能力、确定爆轰产物状态方程参数的基准试验，已经得到了广泛应用^[3-8]。标准的圆筒试验方法主要采用电离探针测量炸药爆速，用高速扫描相机记录定常滑移爆轰驱动下圆筒壁的径向膨胀距离随时间的变化关系，并利用经验公式计算圆筒壁的膨胀速度和比动能等表征炸药做功能力的特征量，最后利用这些试验结果标定JWL(Jones-Wilkins-Lee)状态方程的参数^[9-10]。本文中，利用贝叶斯方法对圆筒试验JWL状态方程参数进行标定，给出标定参数的估计值，得到标定参数的后验分布，然后将该后验分布作为不确定度量化分析时输入参数的初始分布，以减小输入参数分布选择引入的不确定性。

1.   圆筒试验

圆筒试验的试验装置示意图如图 1所示，高速相机VISAR从狭缝扫描记录的是圆筒壁外表面的径向膨胀过程，通过对照片底片的处理获得圆筒壁的径向膨胀位移和时间的关系。得到的观测数据如图 2所示。炸药爆轰产物JWL状态方程形式和等熵条件如下：

图  1  圆筒试验的试验装置示意图

Figure  1.  Device of cylinder test

下载: 全尺寸图片 幻灯片

图  2  圆筒试验径向半径和时间的关系

Figure  2.  Radius vs. time in cylinder test

下载: 全尺寸图片 幻灯片

式中：p为压力，下标“S”表示熵，V为相对比容，e为初始体积能量，A、B、C、R₁、R₂和ω为待标定参数。利用JWL状态方程在C-J点的特性，由圆筒试验和相关试验的观测数据得到爆压、爆速、密度、内能和参数ω，再构造代数方程组，给定R₁和R₂，联立求解A、B和C^[11]。

2.   贝叶斯标定方法

M.Kennedy和A.O‘Hagan^[12]提出用贝叶斯方法进行参数标定和预测。由此发展出来一系列方法，统称为KOH方法^[13-15]。KOH方法用高斯过程(Gaussian process, GP)对计算结果和观测结果进行建模，同时对模型偏差进行建模。KOH方法的主要思想是用高斯过程建立代理模型，然后用代理模型计算后验分布和标定参数。下面简单介绍KOH方法的框架，更多技术细节参见文献[12]。

2.1   定义和统计模型

令X=(X⁽¹⁾, …, X^(q_X))为q_X维控制变量，θ=(θ⁽¹⁾, …, θ^(q_θ))为q_θ维标定参数，Y为数值模拟计算的输出结果，Z为试验观测结果。KOH模型如下：

式中：f_C(·)为计算模型，b(·)为模型偏差，ε_j为观测误差，ρ为回归系数，M为数值模拟样本量，N为试验样本量，x_i^*为控制变量在第i次数值模拟中的取值，x_j为控制变量在第j次实验中的取值(按照数理统计中约定，采用大写字母表示随机变量，相应的小写字母表示该随机变量的取值，本文中，“X”、“Y”、“Z”、“T”、“R₁”、“R₂”均采用该约定)。为了得到数值模拟的结果，需要给标定参数赋值，令s表示数值模拟时输入的标定参数，θ表示标定参数的真实值，以便于区分。

假设ε_j服从均值为零、方差为λ的正态分布；f_C(·)用高斯过程建模，均值函数为m₁(x, s)，协方差函数为c₁(·, ·)；模型偏差b(x)也用高斯过程建模，均值函数为m₂(x)，协方差函数为c₂(·, ·)。假设均值函数和协方差函数的形式如下：

式中：h_i(·)为统计建模时任意给定的连接函数，σ_i²、ω_Xi、ω_θi为协方差函数建模中需要估计的未知参数。令ψ=(ψ₁, ψ₂)，ψ_i是协方差函数c_i(·, ·)中的超参数(σ_i², ω_Xi, ω_θi)，i=1, 2。KOH模型的超参数包括β=(β₁, β₂)和φ=(ρ, λ, ψ), 则全部要估计的参数为(θ, β, φ)。

2.2   数据结构

数据包括观测数据z和数值模拟结果y, 即d^T=(z^T, y^T)。数值模拟过程中输入变量的取值记为D₁={(x₁^*, s₁), …, (x_M^*, s_M)}，试验数据控制变量的取值记为D₂={x₁, …, x_N}，与D₁相对应，可以将D₂扩充为D₂(θ)={(x₁, θ), …, (x_N, θ)}。

2.3   后验分布

令先验分布P(β)∝1。根据贝叶斯理论，在给定先验分布的条件下，后验分布为：

式中：f(d; m_d(θ), V_d(θ))是多元正态分布N(m_d(θ), V_d(θ))的密度函数，且

式中：H₁、H₂为由连接函数h_i构成的向量，V₁(D₁)为由集合D₁中数据构成的方差矩阵，V₂(D₂)为由集合D₂中数据构成的方差矩阵，C₁为D₁和D₂中数据构成的协方差矩阵，I_N为N维单位矩阵。

2.4   计算步骤

直接计算式(6)中的后验分布较困难，因此将计算过程分为4个步骤：步骤1和步骤2是建立GP模型，步骤3是用建立的GP模型计算后验分布，步骤4是利用步骤3得到的后验分布标定参数。

步骤1：建立GP代理模型，估计β₁和c₁(·, ·)的超参数ψ₁。用数值模拟所得数据y估计β₁和ψ₁，得到估计值$\hat{\beta}_{1}$和$\hat{\psi}_{1}$。步骤2：建立观测数据的GP代理模型，估计ρ、λ、β₂和c₂(·, ·)的超参数ψ₂。在$\hat{\beta}_{1}$和$\hat{\psi}_{1}$给定的条件下，用观测数据z估计ρ、λ、β₂和c₂(·)的超参数ψ₂，从而得到估计值$\mathit{\boldsymbol{\widehat \beta }} = \left( {{{\hat \beta }_1}, {{\hat \beta }_2}} \right)$和$\mathit{\boldsymbol{\widehat \varphi }} = (\hat \rho, \hat \lambda, \hat \psi )$。步骤3：计算标定参数θ的后验分布。由于$P(\mathit{\boldsymbol{\theta }}|\mathit{\boldsymbol{\varphi }} = \mathit{\boldsymbol{\widehat \varphi }}, \mathit{\boldsymbol{d}}) \propto P(\mathit{\boldsymbol{\theta }}, \mathit{\boldsymbol{\widehat \varphi }}|\mathit{\boldsymbol{d}})$，因此通过计算$P(\mathit{\boldsymbol{\theta }}, \mathit{\boldsymbol{\widehat \varphi }}|\mathit{\boldsymbol{d}})$即可得到θ的后验分布。步骤4：利用θ的后验分布标定参数。标定参数的估计可以定义为后验分布的均值${\mathit{\boldsymbol{\widehat \theta }}^{({\rm{mean}})}}$、中位数${\mathit{\boldsymbol{\widehat \theta }}^{({\rm{mid}})}}$或最大概率值${\mathit{\boldsymbol{\widehat \theta }}^{({\rm{max}})}}$，其中${\mathit{\boldsymbol{\widehat \theta }}^{({\rm{max}})}}$也称为最大后验估计(maximum posterior estimates, MPE)。

3.   JWL状态方程参数的标定

采用KOH贝叶斯方法标定圆筒试验JWL状态方程参数。将标定参数记为(R₁, R₂)；将时间视为设计变量，记为T；同时令M=30, N=10。数值计算过程中，R₁和R₂的抽样区间分别为[4.3, 5.5]和[0.8, 1.5]，用拉丁超立方抽样(latin hypercube sample, LHS)方法得到M组参数并代入模拟程序进行计算，得到M条径向半径随时间变化的曲线。对设计变量T进行LHS抽样得到M个样本。将T和(R₁, R₂)的样本随机配对得到模拟计算的输入数据D₁={(t₁^*, R₁₁, R₂₁), …, (t_M^*, R_1M, R_2M)}。D₁对应的计算输出结果为y={y₁, …, y_M}。需要注意的是，每次数值计算的结果都是径向半径随时间变化的曲线，这里只是从曲线上随机取一个抽样点作为输出结果。对时间变量T随机抽样得到N个样本，即D₂={t₁, …, t_N}。将D₂扩充为D₂(r₁, r₂)={(t₁, r₁, r₂), …, (t_N, r₁, r₂)}，其中(r₁, r₂)是标定参数的真实值。D₂(r₁, r₂)对应的径向半径观测数据为z={z₁, …, z_N}。假设(R₁, R₂)的先验分布为均匀分布，即对(R₁, R₂)没有任何知识；超参数ψ的先验分布为标准正态分布。数值模拟使用的计算程序为由北京应用物理与计算数学研究所自主研制的相容拉氏流体动力学工具包CHAP (compatible hydrodynamics analysis program)^[16]。

下面分别采用两个模型进行建模，其主要区别在于模型中是否包含偏差。模型1是KOH模型的通用形式，如式(2)~(3)所示；模型2是KOH模型的特殊形式，相当于在模型1的基础上令ρ=1，b(x_j)≡0，此时式(3)可以简化为：

按照第2.4节中的4个步骤依次进行计算。第1步是建立数值模拟的GP代理模型。GP代理模型中参数的估计值分别为${\hat \beta _1} = 3.090, \hat \sigma _1^2 = 2.296, {\hat \omega _{X1}} = 0.603, {\hat \omega _{{R_1}1}} = 0.064, {\hat \omega _{{R_2}1}} = 0.059$。为检验GP代理模型预测的能力，重新随机抽样4组(R₁, R₂)样本代入数值模拟程序进行计算，将计算得到的结果和采用GP模型预测的结果进行比较，如图 3所示。从图 3可以看出：代理模型预测曲线和数值计算的曲线基本重合，说明代理模型的预测足够精确。第2步是建立观测数据的GP代理模型。模型1中模型参数的估计值分别为$\hat \rho = 1.057, \hat \lambda = 1.810, {\hat \beta _2} = - 0.607, \hat \sigma _2^2 = 1.585, {\hat \omega _{X2}} = 2.220$。由于模型2忽略了模型偏差并且ρ=1，所以该步骤中只需要估计λ，估计值为$\hat{\lambda}=2.608 \times 10^{-5}$。第3步是计算标定参数的后验分布。图 4给出了标定参数(R₁, R₂)的后验概率分布等高线，其中后验概率共分为9个概率区间，分别由9种颜色进行区分，从小到大依次用P_k (k=1, …, 9)表示，P₉表示最大概率区间的值域。图 4(a)为采用模型1得到的后验概率分布，可以看出，概率分布从右下角向左上角逐级递增，概率最大区域P₉位于左上角；图 4(b)为采用模型2得到的后验概率分布，可以看出，概率分布从右下角和左边向中间区域逐级递增，P₉近似位于中间位置。第4步是标定参数。由图 4可知：模型1的MPE为$\hat R_1^{(\max )} = 4.68258, \hat R_2^{(\max )} = 1.31066$；模型2的MPE为$\hat R_1^{(\max )} = 4.74985, \hat R_2^{(\max )} = 1.10314$。两个模型的MPE均近似位于P₉区域的中心位置。

图  3  GP代理模型预测的结果和数值计算结果的比较

Figure  3.  Comparison of predictions of the GP surrogate model with numerical simulations

下载: 全尺寸图片 幻灯片

图  4  两模型的后验概率等高线和最大后验估计

Figure  4.  Posterior distribution and MPE with different models

下载: 全尺寸图片 幻灯片

图 5为模型1和模型2的预测效果，其中标定结果(红线)是将MPE代入数值模拟程序计算得到，预测结果(绿线)采用观测数据的GP代理模型计算得到，试验结果用黑线表示。图 5显示：两个模型的预测结果均与试验结果符合很好，说明建立的代理模型都能够比较精确地描述试验数据。模型2中试验结果和标定结果符合得比较好，表明标定的参数值是最佳拟合的估计；但是，当t>39μs时，模型1的标定结果和试验结果差异逐渐增大。由于数值模拟和基于观测数据代理模型的预测(即${\hat f_{\rm{C}}}\left( {{\mathit{\boldsymbol{x}}_j}, \mathit{\boldsymbol{\widehat \theta }}} \right)$和${\hat f_{\rm{C}}}\left( {{\mathit{\boldsymbol{x}}_j}, \mathit{\boldsymbol{\widehat \theta }}} \right) + \hat b\left( {{\mathit{\boldsymbol{x}}_j}} \right)$)都足够精确，所以两者的差异即为模型偏差的预测$\hat b(\mathit{\boldsymbol{x}})$，说明模型1标定的参数值是“调和”了模型偏差$\hat b(\mathit{\boldsymbol{x}})$和标定参数$\widehat{\mathit{\boldsymbol{ }}\!\!\theta\!\!\text{ }}$的最佳估计。

图  5  径向半径的预测、试验和模拟结果的比较

Figure  5.  Comparison of radius histories obtained by prediction, experiment and simulation

下载: 全尺寸图片 幻灯片

4.   结论

(1) 用KOH贝叶斯方法标定圆筒试验JWL状态方程的参数R₁和R₂，得到了标定参数的估计值和后验分布，并将后验分布作为子系统级(或系统级)复杂系统数值模拟不确定度量化分析的输入参数的初始分布，以尽可能地减少分布选择引入的认知不确定性。(2)分析了两种KOH统计模型假设对标定参数估计值和后验分布的影响。从参数标定的角度分析，模型2得到的估计值更好；从不确定度量化的角度分析，由于模型1还分析了模型偏差引入的不确定性，因此模型1更合适。