0.   引言

土水特征曲线（SWCC）代表了非饱和土体积含水量随基质吸力的变化规律。针对非饱和土，确定SWCC对于估算非饱和土特性（如非饱和抗剪强度和渗透系数）至关重要^[1-2]。通过指定在SWCC试验期间的测点数量及其相应的控制变量值（例如采用压盘法时的基质吸力^[3]）设定相应的试验方案，进一步按照预先设定的试验方案可通过各种现场或室内试验确定土样SWCC。确定SWCC通常十分耗时，例如，Gatabin等^[4]研究后认为，确定膨润土SWCC的单个数据点最多需要50 d。因此，从单个SWCC试验获得的试验数据的数量通常是有限的。美国农业部Salinity试验室开发的非饱和土数据库（UNSODA）包含了730个土体样品的脱湿SWCC的试验数据，其中95%的土体样品（即691个土体样品）的数据点数少于20个。尽管可以通过有限数量的试验数据和参数模型^[5-6]确定SWCC，但估计SWCC的不确定性不可避免，其不确定性直接影响非饱和土参数（例如抗剪强度和渗透系数）和岩土工程可靠性^[7]。基于特定地点的SWCC数据，已有研究采用贝叶斯方法探讨了SWCC的不确定性^[8-9]。

如图 1所示，估计SWCC的不确定性大小取决于从试验方案获得的试验数据。利用不同基质吸力处的SWCC数据（参见图 1中的空心圆，图 1（a）~1（c）），估计SWCC的95%置信区间（credibility intervals，CIs）可能明显不同，反映了估计SWCC的不同置信度（或不确定性）。随着试验设备和技术的进步，在SWCC试验过程中更容易获得更多的数据点，这提升了SWCC表征的准确性和可靠性。然而，SWCC试验仍存在一些实际问题，即多少个测量点以及如何指定控制变量的值能精确和可靠地估计SWCC，尽可能地缩小估计SWCC的CIs。因此，相比于盲目地增加数据点（需要花费大量时间和精力），通过试验设计确定最优试验方案更加合理。由于在试验之前没有试验数据，并且候选试验方案的数目非常多，甚至有无限个，确定最优的试验方案并不简单，试验设计优化（experimental design optimization，EDO）提供了一个确定最优试验方案的合理工具^[10]。文献中已经发展了多种EDO方法，包括基于经典统计学的传统试验设计优化（classical experimental design optimization，CEDO）方法^[11]和基于贝叶斯推理或信息理论的贝叶斯试验设计优化（Bayesian experimental design optimization，BEDO）方法^[12]。与CEDO相比，BEDO能够量化各种不确定性（如固有变异性、测量误差和统计不确定性），这些不确定性会影响估计SWCC的CIs。同时，BEDO可以对比不同试验方案所提供信息的预期价值，进而评估SWCC不确定性的折减程度。BEDO目前已应用于岩土工程和地质工程中现场取样^[13]和勘察方案^[14]设计。然而BEDO在岩土室内试验设计（例如SWCC试验）中的应用较少。目前仍缺乏一种过程简单、应用方便的SWCC试验设计方法。

图  1  试验方案对SWCC置信区间的影响

Figure  1.  Effects of experimental schemes on CIs of SWCC

下载: 全尺寸图片 幻灯片

本文提出了一种SWCC试验单阶段贝叶斯试验设计优化（OBEDO）方法，该方法应用十分方便，采用单阶段优化的方式直接确定SWCC试验的最优试验方案，以减少估计SWCC的不确定性。所提方法首先对SWCC分区，并对各个区间离散，获取所有候选试验方案的集合（即设计空间）。候选试验方案分为控制点及附加点两种，控制测点确定SWCC的主要趋势，附加点精确描述SWCC。由于在试验设计阶段没有可用的测量数据，所提出的OBEDO方法基于给定参数模型（例如VGB）的先验信息，模拟每个候选方案的试验数据，并采用模拟数据计算候选方案的期望效用。进一步利用SSO^[15]搜索不同测点数目条件下具有最大化期望效用的候选方案，并将其确定为最优设计方案。值得指出的是，在SSO优化过程中，直接生成包含控制点及附加点的样本（即候选试验方案），采用单阶段优化的方式确定最优试验方案，简化试验设计过程。

应注意到，SWCC试验的控制变量可以为含水率或吸力值，本文所提OBEDO方法以基质吸力为控制变量，最优试验方案由基质吸力测点构成。为了达到最佳的试验精度，最大程度减少所估计SWCC的不确定性，本文采用期望效用量化SWCC试验结果的不确定性，利用SSO搜索具有最大期望效用的试验方案，并将该方案确定为最优试验方案，包括最优测点数目及对应的基质吸力值。该最优试验方案在保障最佳试验精度（估计SWCC的不确定性最小）前提下确定的最优测点数目。本文首先描述了OBEDO的框架，然后介绍了如何生成候选试验方案的设计空间及量化候选试验方案的期望效用，并进一步说明了如何采用SSO搜索具有最大期望效用的试验方案。最后以SWCC试验设计为例，说明了所提方法的有效性。

1.   SWCC试验单阶段贝叶斯优化设计框架

本文提出的SWCC试验单阶段贝叶斯优化设计方法（OBEDO）的框架如图 2所示：首先，针对研究的土体类型，在试验前收集SWCC的先验信息（即主要的SWCC模型及其模型参数的典型范围）以及采用试验设备的信息。后续设计过程包括确定候选试验方案设计空间、计算候选试验方案的期望效用、利用SSO算法确定最优试验设计方案3个步骤。需要指出的是，利用收集到的信息确定设计空间时，设计空间应与SWCC试验所采用的设备和技术提供的控制变量的精度和范围相适应。下面将详细介绍OBEDO框架的3个步骤。

图  2  SWCC试验单阶段贝叶斯优化设计（OBEDO）框架

Figure  2.  One-stage Bayesian experimental design optimization (OBEDO) framework for SWCC tests

下载: 全尺寸图片 幻灯片

2.   基于VGB模型的候选试验方案设计空间

对于给定的SWCC参数模型，进气值${\psi _{\text{b}}}$，残余基质吸力${\psi _{\text{r}}}$和基质吸力反弯点${\psi _{\text{i}}}$及其相应的饱和度将SWCC（即AB）划分为AP_b，P_bP_i，P_iP_r，P_rB4个部分，如图 3所示。为了确定SWCC，需要4个控制测量点C₁，C₂，C₃，C₄，且4个点对应的基质吸力应位于4个区间$(0,{\psi _{\text{b}}}]$，$({\psi _{\text{b}}},{\psi _{\text{i}}}]$，$({\psi _{\text{i}}},{\psi _{\text{r}}}]$，$({\psi _{\text{r}}},{\text{ + }}\infty )$之内。但是，4个控制测点C₁~C₄只能确定SWCC的主要趋势，不能精确地刻画SWCC局部细节。为了进一步提高SWCC的精度，需要附加测点D₁~D_n-4，其分布见图 3所示。令总测点数目为n（n为4~N的整数），其中N为试验设计中最大的总测点数目。则控制测点数目为4，附加测点数目则为n-4。候选试验方案E_n由两部分组成，一部分为控制测点C₁~C₄对应的基质吸力值（即${\psi _{{C_1}}}$，${\psi _{{C_2}}}$，${\psi _{{C_3}}}$，${\psi _{{C_4}}}$），另一部分为附加测点D₁~D_n-4对应的基质吸力值，即${\psi _{{D_j}}}$(j=1, 2, …, n-4)。

图  3  附加点D₁，D₂，…，D_n-4分布

Figure  3.  Positions of additional points D₁, D₂, …, D_n-4

下载: 全尺寸图片 幻灯片

在试验设计阶段，无试验数据，因此与给定的SWCC模型相对应的${\psi _{\text{b}}}$，${\psi _{\text{i}}}$，${\psi _{\text{r}}}$的值均是未知的。在本文的研究中，基于给定的SWCC模型，利用变量${\psi _{\text{b}}}$，${\psi _{\text{i}}}$，${\psi _{\text{r}}}$的先验信息（例如SWCC模型参数的典型范围）确定其期望值（或平均值）。例如对于给定的VGB模型，其对应的${\psi _{\text{b}}}$，${\psi _{\text{i}}}$，${\psi _{\text{r}}}$取值满足公式^[16-18]：

式中：${\alpha _{vb}}$，${n_{vb}}$为VGB模型参数；${S_{{\text{e,i}}}}$为基质吸力${\psi _{\text{i}}}$对应的有效饱和度；k₁为基质吸力${\psi _{\text{i}}}$处的斜率；ψ'为基质吸力大于${\psi _{\text{r}}}$时SWCC开始线性下降时的基质吸力；${S'_{\text{e}}}$为基质吸力$\psi '$对应的有效饱和度；k₂为点$\left( {\psi ',{{S'}_{\text{e}}}} \right)$的斜率。上述符号如图 4所示。

图  4  典型土水特征曲线

Figure  4.  Typical SWCC

下载: 全尺寸图片 幻灯片

利用式（1）~（3）确定${\psi _{\text{b}}}$，${\psi _{\text{i}}}$，${\psi _{\text{r}}}$时，需要提供参数${\alpha _{vb}}$，${n_{vb}}$的相关信息。在试验设计阶段，${\alpha _{vb}}$，${n_{vb}}$为未知参数，无法获取${\psi _{\text{b}}}$，${\psi _{\text{i}}}$，${\psi _{\text{r}}}$的具体值。因而，采用它们各自的平均值（即${\bar \psi _{\text{i}}}$，${\bar \psi _{\text{b}}}$，${\bar \psi _{\text{r}}}$）替代${\psi _{\text{i}}}$，${\psi _{\text{b}}}$，${\psi _{\text{r}}}$。${\bar \psi _{\text{i}}}$，${\bar \psi _{\text{b}}}$，${\bar \psi _{\text{r}}}$利用参数${\alpha _{vb}}$，${n_{vb}}$的先验信息估计。基于参数α_vb，n_vb已知的典型范围（即[α_{vb, min}, α_{vb, max}]和[n_{vb, min}, n_{vb, max}]），根据最大熵原理定义它们的均匀分布，并将该均匀分布作为参数的先验分布。在此基础上，然后通过Monte Carlo模拟，从参数${\alpha _{vb}}$，${n_{vb}}$的均匀先验分布中模拟N_p个随机样本，利用这些随机样本和式（1）~（3）计算${\psi _{\text{i}}}$，${\psi _{\text{b}}}$和${\psi _{\text{r}}}$的N_p个估计值。根据${\psi _{\text{i}}}$，${\psi _{\text{b}}}$和${\psi _{\text{r}}}$的N_p个估计值，计算它们各自的平均值（即${\bar \psi _{\text{i}}}$，${\bar \psi _{\text{b}}}$和${\bar \psi _{\text{r}}}$）。控制测点C₁，C₂，C₃，C₄的基质吸力${\psi _{{C_1}}}$，${\psi _{{C_2}}}$，${\psi _{{C_3}}}$，${\psi _{{C_4}}}$的可能区间分别取为$(0,{\bar \psi _{\text{b}}}]$，$({\bar \psi _{\text{b}}},{\bar \psi _{\text{i}}}]$，$({\bar \psi _{\text{i}}},{\bar \psi _{\text{r}}}]$，$({\bar \psi _{\text{r}}},{\psi _m})$，如图 5所示。其中，${\psi _m}$为试验设备能够施加的最大基质吸力，其值大小取决于试验条件。如前所述，4个控制测点只能确定SWCC的总体趋势，为了进一步提高预测的SWCC的精度，需要增加n-4个附加测点D₁~D_n-4。采用${\bar \psi _{\text{i}}}$，${\bar \psi _{\text{b}}}$，${\bar \psi _{\text{r}}}$将SWCC分区后，附加测点D₁~D_n-4分布见图 6所示。附加测点的基质吸力${\psi _{{D_j}}}$(j=1, 2, …, n-4)的取值范围为排出已选控制测点基质吸力${\psi _{{C_1}}}$，${\psi _{{C_2}}}$，${\psi _{{C_3}}}$，${\psi _{{C_4}}}$后，区间$\left( {\left. {0,{{\bar \psi }_{\text{b}}}} \right]} \right.$，$\left( {\left. {{{\bar \psi }_{\text{b}}},{{\bar \psi }_{\text{i}}}} \right]} \right.$，$\left( {\left. {{{\bar \psi }_{\text{i}}},{{\bar \psi }_{\text{r}}}} \right]} \right.$，$\left( {{{\bar \psi }_{\text{r}}},{\psi _{\text{m}}}} \right)$内的剩余值。

图  5  基于${\bar \psi _{\text{b}}}$，${\bar \psi _{\text{i}}}$和${\bar \psi _{\text{r}}}$的SWCC分区

Figure  5.  Partition of SWCC

下载: 全尺寸图片 幻灯片

图  6  控制点和附加点分布

Figure  6.  Positions of control points and additional point

下载: 全尺寸图片 幻灯片

当给定试验仪器时，基质吸力可能值${\psi _{\text{h}}}$的范围可表示为${\mathit\Omega _O} = \left( {{\text{0:}}\Delta {\psi _1}:{{\bar \psi }_{\text{b}}}} \right]$$\cup \left( {{{\bar \psi }_b}{\text{:}}\Delta {\psi _2}:{{\bar \psi }_i}} \right]$ $\cup \left( {{{\bar \psi }_{\text{i}}}{\text{:}}\Delta {\psi _3}:{{\bar \psi }_{\text{r}}}} \right]$ $\cup \left( {{{\bar \psi }_{\text{r}}}{\text{:}}\Delta {\psi _4}:{\psi _m}} \right]$，其中$\Delta {\psi _1}$，$\Delta {\psi _2}$，$\Delta {\psi _3}$，$\Delta {\psi _4}$分别为区间$\left( {\left. {0,{{\bar \psi }_{\text{b}}}} \right]} \right.$，$\left( {\left. {{{\bar \psi }_{\text{b}}},{{\bar \psi }_{\text{i}}}} \right]} \right.$，$\left( {\left. {{{\bar \psi }_{\text{i}}},{{\bar \psi }_{\text{r}}}} \right]} \right.$，$\left( {{{\bar \psi }_{\text{r}}},{\psi _m}} \right)$的离散精度（例如，可取试验仪器最小基质吸力增量）。在实际SWCC试验过程中，采用仪器基质吸力测量范围不同时，也应该相应地采用不同的$\Delta {\psi _{\text{i}}}$(i=1, 2, 3, 4)。在本文中，$\Delta {\psi _{\text{i}}}$ (i=1, 2, 3, 4)的取值随区间长度变化。上述离散方法生成${N_O}$个基质吸力可能值${\psi _{\text{h}}}$。假设${N_O}$个${\psi _{\text{h}}}$值中有${N_{{C_{\text{1}}}}}$，${N_{{C_{\text{2}}}}}$，${N_{{C_{\text{3}}}}}$，${N_{{C_{\text{4}}}}}$个${\psi _{\text{h}}}$的值分别位于区间$\left( {\left. {0,{{\bar \psi }_{\text{b}}}} \right]} \right.$，$\left( {\left. {{{\bar \psi }_{\text{b}}},{{\bar \psi }_{\text{i}}}} \right]} \right.$，$\left( {\left. {{{\bar \psi }_{\text{i}}},{{\bar \psi }_{\text{r}}}} \right]} \right.$，$({\bar \psi _{\text{r}}},{\psi _m})$之中，位于这4个区间之中的${\psi _h}$的值构成的集合分别表示为${\mathit\Omega _{{C_{\text{1}}}}}$，${\mathit\Omega _{{C_2}}}$，${\mathit\Omega _{{C_3}}}$，${\mathit\Omega _{{C_4}}}$。4个控制测点C₁，C₂，C₃，C₄的基质吸力取值范围分别满足${\psi _{{C_{\text{1}}}}} \in {\mathit\Omega _{{C_{\text{1}}}}}$，${\psi _{{C_2}}} \in {\mathit\Omega _{{C_2}}}$，${\psi _{{C_3}}} \in {\mathit\Omega _{{C_3}}}$，${\psi _{{C_4}}} \in {\mathit\Omega _{{C_4}}}$。同时，令${\mathit\Omega _{D{\text{|}}C}}$表示每个附加测点D_j(j=1, 2, …, n-4)对应的基质吸力${\psi _{{D_j}}}$(j=1, 2, …, n-4)的取值集合，其表示为

任意4个控制点及n-4个附加点基质吸力构成一个完整的候选试验方案E_n，且所有的E_n组成了设计空间${\mathit\Omega _n}$。给定测点数目n时，任意候选试验方案E_n表示为

设计空间${\mathit\Omega _n}$中候选试验方案E_n的总数为${N_{{C_{\text{1}}}}} \cdot {N_{{C_2}}} \cdot {N_{{C_3}}} \cdot {N_{{C_4}}} \cdot C_{{N_o} - 4}^{n - 4}$，本研究主要解决的问题之一为如何在试验空间${\mathit\Omega _n}$中确定具有最高期望效用的候选方案，期望效用的具体定义及优化算法如下文所论述。

3.   候选试验方案期望效用

3.1   期望效用的定义

根据候选试验方案E_n中包含的基质吸力值获取相应的SWCC试验数据S_e={S_{e, k}, k=1, 2, …, n}。利用获取的S_e更新给定的SWCC模型参数分布。相应地，根据SWCC模型参数分布的更新效果（相比于模型参数先验分布）评估数据S_e的信息价值。理论上，采用相对熵$R({E_n})$量化模型参数分布的更新效果。相对熵$R({E_n})$揭示更新后的参数分布$p(\mathit\Theta \left| {{S_\text{e}},{E_n}} \right.)$与参数先验分布$p(\mathit\Theta \left| {{E_n}} \right.)$之间的区别，表示为^[10]

式中：$R({E_n})$为相对熵；$p(\mathit\Theta \left| {{S_\text{e}},{E_n}} \right.)$为SWCC参数后验分布；$p(\mathit\Theta \left| {{E_n}} \right.)$为先验分布。注意到，利用式（6）计算$R({E_n})$需要利用试验数据S_e，但是仍处于试验设计阶段，尚未获取数据S_e。采用候选方案${E_n}$对应试验数据的期望效用$U({E_n})$评估方案的有效性，期望效用$U({E_n})$表示为^[19]

式中：$U({E_n})$为期望效用；$p({S_\text{e}}\left| {{E_n}} \right.)$为方案${E_n}$对应试验数据S_e的概率密度函数(probability density function，PDF)。对于给定的SWCC模型及相应的模型参数$\mathit\Theta$，$p({S_\text{e}}\left| {{E_n}} \right.)$表示为

式中：$p({S_\text{e}}\left| {\mathit\Theta ,{E_n}} \right.)$为给定SWCC模型参数$\mathit\Theta$的条件下，方案${E_n}$对应试验数据S_e的概率密度函数（PDF）；S_{e, k}为方案${E_n}$中第k个测点${\psi _k}$对应的饱和度，其表示为

式中：$M(\mathit\Theta ,{\psi _{\text{k}}})$为基质吸力为${\psi _{\text{k}}}$时的模型预测值；高斯随机变量$\varepsilon$的期望为0，标准差为${\sigma _\varepsilon }$，代表了由于测量误差及模型误差诱发的不确定性，通常假设不同测量值处的$\varepsilon$相互独立^[20]。根据式（9），$p({S_{\text{e}}}\left| {\mathit\Theta ,{E_n}} \right.)$可表示为

根据式（7）~（9），（10），可以利用模拟数据计算$U(E_n^{})$，详见3.2节。

3.2   基于随机模拟计算期望效用

基于式（7）~（9），（10）计算候选方案${E_n}$的期望效用$U(E_n^{})$需要模拟数据。模拟数据依据候选试验方案${E_n}$生成，代表了现场SWCC试验可能获取的试验数据值。在本研究中，使用祖先取样^[21]模拟候选方案${E_n}$对应的SWCC试验数据。具体地，从先验分布$p(\mathit\Theta \left| {{E_n}} \right.)$中抽取N_r个$\mathit\Theta$的样本，在无信息条件下$p(\mathit\Theta \left| {{E_n}} \right.)$为参数$\mathit\Theta$的联合均匀分布^[8]。本研究中$\mathit\Theta$不仅包括SWCC参数，还包含$\varepsilon$的标准差。利用式（10）及$\mathit\Theta$的每一个样本生成候选方案${E_n}$中n个基质吸力对应的S_e，最终获取N_r个S_e（即${S_{{\bf{e}},r}}$(r=1, 2, …, N_r)）的样本值。利用N_r组$\mathit\Theta$和S_e的样本值计算U(${E_n}$)的值。将式（6）代入式（7）后得到

利用$\mathit\Theta$和S_e的样本值，式（11）写为

式中：$p({\mathit\Theta _r}\left| {{S_{{\text{e}}\,,r}},{E_n}} \right.){\kern 1pt}$为给定${S_{{\text{e}},r}}$条件下，${\mathit\Theta _r}$的后验分布。在贝叶斯框架下，$p({\mathit\Theta _r}\left| {{S_{{\text{e}}\,,r}},{E_n}} \right.){\kern 1pt}$可表示为

将式（13）代入式（12）后得

分别利用式（8），（10）以及样本值$\mathit\Theta$=${\mathit\Theta _r}$和${S_{\text{e}}} = {S_{{\text{e}},r}}$ (r=1, 2, …, N_r)计算${\kern 1pt} p({S_{{\text{e}},r}}\left| {{E_n}} \right.)$和${\kern 1pt} p({S_{{\text{e}},r}}\left| {{\mathit\Theta _r},{E_n}} \right.)$。确定${\kern 1pt} p({S_{{\text{e}},r}}\left| {{E_n}} \right.)$时，使用模拟的方法计算所涉及的积分。在本文中，利用直接蒙特卡洛方法^[19]计算${\kern 1pt} p({S_{{\text{e}},r}}\left| {{E_n}} \right.)$

式中：${\mathit\Theta _t}$为利用先验分布$p(\mathit\Theta \left| {{E_n}} \right.)$抽样获取的样本，t=1, 2, …, N_t。

在给定测点数目n的条件下，最优试验设计方案$E_n^*$具有最大的期望效用U(${E_n}$)。

第4节将论述如何采用SSO在设计空间识别设计方案$E_n^*$，使其具有最大期望效用。

4.   基于SSO的SWCC试验方案优化

在候选试验空间中确定方案$E_n^*$，使其具有最高的期望效用$U(E_n^*)$，这一问题可表达为

式中：${\mathit\Omega _{{C_i}}}$，${\mathit\Omega _{D{\text{|}}C}}$分别为${\psi _{{{\text{C}}_i}}}$(i=1, 2, 3, 4)及${\psi _{{D_j}}}$(j=1, 2, …, n-4)的设计空间（见第2节）。如前所述，确定设计方案$E_n^*$需要求解式（17），在本文中，采用SSO在设计空间搜索最优方案$E_n^*$。为了搜索U(${E_n}$)的最大值，${\psi _{{C_i}}}$和${\psi _{{D_j}}}$均视为可行域中的均匀随机变量。每一组基质吸力样本构成的集合确定任意一个试验方案$E_n^{}$，同时计算其相应的$U(E_n^{})$。令$U(E_n^*)$表示$U(E_n^{})$的最大值，$U(E_n^*)$即为最优设计方案$E_n^*$的期望效用$U(E_n^{})$。由于$E_n^*$是全局优化的试验方案，理论上其被超越的概率为0。因此，通过求解如下可靠性分析问题，在候选试验方案中搜索最优方案$E_n^*$^[15]。

式中：$F = \{ U({E_n}) > U(E_n^*)\}$为辅助目标事件。为了搜索候选试验方案的可行域（即${\mathit\Omega _{{C_i}}}{\kern 1pt}$(i=1, 2, 3, 4)和${\mathit\Omega _{D{\text{|}}C}}$）内的事件F，子集模拟逐步生成一系列满足中间失效事件的条件样本，且这些中间事件满足${F_1} \supset {F_2} \supset \cdots$ $\supset {F_{{N_s}}}{\text{ = }}F$。P(F)表示为

式中：${F_m} = {\text{\{ }}U({E_n}) > {U_m}{\text{\} }}$(m=1, 2, …, N_s)；$P({F_1})$等于$P\left\{ {U({E_n}) > {U_1}} \right\}$；$P(F_m^{}\left| {F_{m - 1}^{}} \right.)$等于$P(U({E_n}) > {U_m}$ $\left| {U({E_n}) > {U_{m - 1}}} \right.)$ (m=2, 3, …, N_s)；${U_1} < {U_2} < \cdots < {U_{{N_s}}} =$ $U(E_n^*)$为N_s个逐渐递增的中间阈值。样本估计值$P\left( {{F_1}} \right)$和$P(F_m^{}\left| {F_{m - 1}^{}} \right.)$始终对应条件概率的指定值p₀（例如，0.1）。SSO算法和实现的更多细节可参考文献[15]。

5.   算例

5.1   候选试验方案

在本例中，模型参数的典型范围^[8]${\alpha }_{vb}\in (0\text{ kPa}，\text{ }$ $20{\text{ kPa}}]$，${n_{vb}} \in (2,10]$，${\sigma _\varepsilon } \in (0,1]$作为各自参数的先验信息。如第2节所述，根据模型参数典型范围定义的均匀先验分布（由其典型范围定义）模拟了N_p个模型参数的样本。同时，使用式（1）~（3）计算N_p个${\psi _{\text{b}}}$，${\psi _{\text{i}}}$和${\psi _{\text{r}}}$的估计值，并在此基础上计算${\bar \psi _{\text{b}}}$，${\bar \psi _{\text{i}}}$，和${\bar \psi _{\text{r}}}$的值。图 7所示为${\bar \psi _{\text{b}}}$，${\bar \psi _{\text{i}}}$和${\bar \psi _{\text{r}}}$随N_p的变化。当N_p大于2000时，${\bar \psi _{\text{b}}}$，${\bar \psi _{\text{i}}}$和${\bar \psi _{\text{r}}}$的估计值趋于稳定，因而在本例中取${\bar \psi _{\text{b}}}$=16 kPa，${\bar \psi _{\text{i}}}$=24 kPa，${\bar \psi _{\text{r}}}$=94 kPa。

图  7  利用先验信息确定${\bar \psi _{\text{b}}}$，${\bar \psi _{\text{i}}}$和${\bar \psi _{\text{r}}}$

Figure  7.  Determination of ${\bar \psi _{\text{b}}}$, ${\bar \psi _{\text{i}}}$ and ${\bar \psi _{\text{r}}}$ based on prior knowledge

下载: 全尺寸图片 幻灯片

例如，假设SWCC试验设备允许的基质吸力测量范围为0 kPa到2000 kPa，即${\psi _m}$=2000 kPa。依据本文第1节所述，${\bar \psi _{\text{b}}}$=16 kPa，${\bar \psi _{\text{i}}}$=24 kPa，${\bar \psi _{\text{r}}}$=94 kPa将SWCC划分为$(0,16]$，$(16,24]$，$(24,94]$和$(94,2000)$4个区间，这4个区间的离散尺寸分别取为$\Delta {\psi _1}$=2 kPa，$\Delta {\psi _2}$=2 kPa，$\Delta {\psi _3}$=2 kPa，$\Delta {\psi _4}$=50 kPa。如前所述，控制点的基质吸力为${\psi _{{C_i}}}$，其可能值分别为${\mathit\Omega _{{C_1}}} =$ $\left\{2，4，\cdots ，14，16\right\}$，${\mathit\Omega }_{{C}_{2}}=\left\{18，20，22，24\right\}$，${\mathit\Omega _{{C_3}}} =$ $\{26，$ $28$，$\cdots ，$$92，94\}$和${\mathit\Omega _{{C_4}}} =$$\left\{144，194，\cdots ，1944，1994\right\}$（单位：kPa）。附加点的基质吸力为${\psi _{{D_j}}}$可能值为${\mathit\Omega _{D{\text{|}}C}} =$ $\{ {\psi _{{D_j}}}\left| {{\psi _{{D_j}}} \in } \right.{\mathit\Omega _O}{\kern 1pt} {\text{and}}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\psi _{{D_j}}} \ne {\psi _{{C_i}}}\}$。

5.2   SWCC试验最优试验方案

确定${\mathit\Omega _{{C_i}}}$(i=1, 2, 3, 4)及${\mathit\Omega _{D|{\text{C}}}}$之后，基于SSO优化最优试验方案，其中p₀，N_s分别取0.1，20，每层模拟采用2000个样本。图 8所示为当测点数目为n=8时，随着模拟层数m的增加，期望效用的中间阈值（即${U_m}$）的变化。在m=11时，数值增加并稳定于1.72。为了保证${U_m}$收敛，将SSO执行到第20层（即m=20），之后将最后一层的${U_m}$视为$U(E_8^*)$的估计值，并将相应的$E_8^*$视为n=8时具有最大期望效用值的设计方案。

图  8  SSO过程中期望效用中间阈值的变化

Figure  8.  Evolution of intermediate threshold value of expected utility during SSO

下载: 全尺寸图片 幻灯片

为了考虑测点数目n对试验方案所包含数据价值的影响，考虑了一系列n值，包括4，5，6，7，8，9，10，11，12，13，14，15，20，25，对于每一个n值，以${p_0}$=0.1，N_s=20执行子集模拟优化，每层选取2000个样本。在不同模拟层下期望效用的中间阈值（即${U_m}$）的变化如图 9所示。图 10展示了在不同测点数目n的条件下最大化期望效用$U(E_n^*)$随着n的变化规律。结果表明，当n小于13时，$U(E_n^*)$的值迅速增加。当n > 13时，继续增加更多的测量点，$U(E_n^*)$几乎不变。最优测点数目n确定为13。当n=13时，最大期望效用值为2.25，对应的基质吸力值为{6，8，14，18，22，30，58，86，294，394，744，944，1444}（单位为kPa）。因此，测点数为13的方案为最优设计方案。若测点数目小于13，则不能达到最佳的试验精度；若将测点数目确定为13，但依据工程经验确定试验方案，根据试验数据的效用值该方案很难优于所提方法确定的最优设计方案，最终也难以达到最佳试验精度；若测点数目大于13，由于试验精度已经趋于稳定，浪费了试验的时间及人力成本。因而测点数目13是保障最佳试验精度前提下所需的最优测点数目。

图  9  测点数目n对应的期望效用中间阈值U_m的变化

Figure  9.  Evolution of intermediate threshold value of expected utility U_m for different measuring points n

下载: 全尺寸图片 幻灯片

图  10  不同测点数目n时的最大期望效用

Figure  10.  Expected utilities of optimal experimental schemes with different numbers of measuring points n

下载: 全尺寸图片 幻灯片

5.3   实测数据校核

本小节采用黄土SWCC试验数据验证所提方法的有效性。分别在文献[22，23]中选择一组黄土数据，如图 11所示，文献[22，23]中土体样本分别有15，16个数据点。注意到，文献中的原始试验数据包括体积含水量$\theta$及其相应的基质吸力$\psi$，将体积含水量$\theta$转换为有效饱和度${S_{\text{e}}} = (\theta - {\theta _{\text{r}}})/({\theta _{\text{s}}} - {\theta _{\text{r}}})$，其中${\theta _{\text{s}}}$是饱和体积含水量，${\theta _{\text{r}}}$是残余体积含水量。

图  11  SWCC实测数据

Figure  11.  Measured data of SWCC

下载: 全尺寸图片 幻灯片

如5.1节和5.2节所述，在给定VGB模型和模型参数（即${\alpha _{vb}}$，${n_{vb}}$）的先验信息的条件下，13个测量点足以表征SWCC。本文将利用所提方法获得的最优试验方案（即{6，8，14，18，22，30，58，86，294，394，744，944，1444}（单位kPa））称之为贝叶斯最优设计方案（bayesian optimal scheme，BOS），并将BOS与实测的测量点比较。为了对比，从文献[22]中的15个数据点及文献[23]中的16个数据点中分别随机选择13个点，用13个测量点模拟试验方案，分别得到$C_{15}^{13}$=105及$C_{16}^{13}$=560个随机试验方案（random experimental schemes，RES）。每个RES包含13组${S_\text{e}}$，$\psi$的值，在此基础上使用式（7）及均匀先验分布计算RES的效用（即R（RES））。图 12所示为文献[22，23]的RESs的效用值。在文献[22，23]的效用值中，RES最大的效用分别为1.77，1.83，该方案为随机最优方案（random optimal scheme，ROS），如图 12中的点划线所示。ROS的效用小于本文所提方法获得的BOS的期望效用（即2.25），表明依据BOS获取的试验数据含有更高的信息价值。为了对比基于BOS及ROS获得的SWCC曲线，基于给定的模型参数的典型范围（${\alpha _{vb}}$∈(0 kPa, 20k Pa]，n_vb∈(2, 10]），随机抽取500组样本（${\alpha _{vb}}$和n_vb），计算500组样本在BOS对应基质吸力点（{6，8，14，18，22，30，58，86，294，394，744，944，1444}（单位kPa））的500个有效饱和度模拟值，并统计相应基质吸力点处有效饱和度的平均值，将其作为BOS对应的模拟数据。限于文章篇幅，本文采用贝叶斯方法确定了基于BOS模拟数据及ROS实测数据（文献[22]）估计的SWCC，见图 13所示。基于BOS模拟数据及ROS实测数据估计的SWCC模型参数后验最可能值分别为${\alpha _{vb}}$=17.76，n_vb=2.29及${\alpha _{vb}}$=5.614，n_vb=1.45。图 13对比了基于BOS模拟数据和ROS实测数据获得的SWCC置信区间。由图 13可知，在高基质吸力区间（图中阴影部分），基于ROS实测数据获得的SWCC置信区间较宽，表明此区间内估计的SWCC不确定性较大。相反地，此区间内基于BOS模拟数据获得的SWCC不确定性相对较小。ROS在高基质吸力区间测点分布少，估计的SWCC不确定性较大；而BOS在该区间合理地布置了测点，有效降低了估计的SWCC在该基质吸力范围内的不确定性，验证了所提方法的有效性。因此，与任意布置测点相比，本文提出的OBEDO方法利用先验信息和仪器信息为SWCC试验设计提供了一种合理的工具。

图  12  随机试验方案的效用

Figure  12.  Utilities of random experimental schemes

下载: 全尺寸图片 幻灯片

图  13  基于BOS模拟数据和ROS实测数据获得的SWCC置信区间对比

Figure  13.  CIs of SWCC with simulated data of BOS and measuring data of ROS

下载: 全尺寸图片 幻灯片

5.4   先验信息对试验设计结果的影响

采用本文所提方法确定最优试验方案需要指定模型参数的先验范围，本文第5.1节中算例的先验范围（α_vb∈(0 kPa, 20 kPa]，n_vb∈(2, 10]）依据文献指定（先验范围Ι），为了进一步探究先验范围对试验结果的影响，本文增加了另外两组具有不同先验范围的试验设计算例，其中一组扩大文献指定的先验范围（先验范围Ⅱ），参数先验范围满足α_vb∈(0 kPa, 30 kPa]，n_vb∈(2, 13]；另外一组缩小文献指定的先验范围（先验范围Ⅲ），参数先验范围满足α_vb∈(5 kPa, 15 kPa]，n_vb∈(3.5, 8.5]。当给定先验范围Ⅱ时，进气值、基质吸力反弯点、残余基质吸力分别为18，27，98 kPa。若给定先验范围Ⅲ，进气值、基质吸力反弯点、残余基质吸力分别为13，22，88 kPa。针对新增的两组试验设计算例，采用期望效用（式（14））量化候选试验方案所包含测点的信息价值，并且基于SSO优化最优试验方案，其中p₀和N_s分别取0.1，20，每层模拟采用2000个样本。图 14展示了3种不同先验范围情况下测点数目n对应的最大期望效用值。由图 14可知，当给定测点数目n时，先验范围Ⅱ对应的最大期望效用值最大，这是因为与先验范围Ⅰ及先验范围Ⅲ相比，先验范围Ⅱ所提供的信息量相对较少、参数不确定性较高，因此相同数量的测点信息价值更高，不确定性缩减更显著。先验范围Ⅰ、先验范围Ⅱ及先验范围Ⅲ确定的最优测点数目分别为13，15和10，其对应的最优设计方案见图 15所示。先验范围Ⅲ在试验之前提供了更确定的SWCC有效信息，因而能够使用更少的测点最大程度降低SWCC的不确定性。

图  14  不同先验信息条件下测点数目n对应的最大期望效用值

Figure  14.  Maximum expected utilities using different prior knowledge

下载: 全尺寸图片 幻灯片

图  15  不同先验范围对应的最优设计方案BOS

Figure  15.  Optimal design scheme BOS using different prior knowledge

下载: 全尺寸图片 幻灯片

6.   结论

本文提出了一种SWCC试验单阶段贝叶斯优化设计（OBEDO）方法，该方法利用先验信息和试验装置信息确定SWCC试验的最优试验方案。所提方法通过离散试验控制变量生成试验方案设计空间，空间中的候选试验方案中的测点由控制点及附加点构成，其作用分别为控制SWCC的主要趋势和降低SWCC的不确定性。然后利用期望效用量化候选试验方案对应数据的价值，并基于子集模拟优化（SSO）方法搜索具有最大期望效用的候选试验方案（即最优试验方案），在SSO优化过程中，直接生成包含控制点及附加点的样本（即候选试验方案），采用单阶段优化的方式直接确定最优试验方案（包括最优测点数目及基质吸力测量值）。文中对所提方法进行了理论推导，并用一个SWCC设计实例说明了所提方法得到3点结论。

（1）对于给定的测点数目n，在SSO过程中，搜索到的试验方案的期望效用随着子集模拟层数的增加而逐渐增加，最终收敛到一个最大值，此时可确定给定n条件下的最大化期望效用值及相应试验方案。

（2）逐渐增大测点数目n，试验方案的最大期望效用值也逐步增加，当测量点的数目n足够大（例如在本文实例中为13）时，最大期望效用基本保持稳定。此时的测点数目n及最大期望效用对应的试验方案即为最优设计方案。

（3）文中采用两组实测黄土SWCC试验数据验证所提方法的有效性，利用试验数据生成随机试验方案，模拟SWCC试验中测点的任意组合。结果表明：SWCC试验中测点的任意组合的效用值很难超越本文确定的最优试验方案（以期望效用为评估指标）。

基于先验知识和试验仪器信息，所提出的OBEDO方法为SWCC试验方案设计提供了一个合理的工具，有助于获得具有较高信息价值的SWCC试验数据，降低估计SWCC的不确定性。