CN113435022B - 基于多道面波勘探的岩土体参数二维空间变异性表征方法 - Google Patents

Abstract

本发明涉及一种基于多道面波勘探的岩土体参数二维空间变异性表征方法，属于岩土工程勘察领域。建立表征岩土体参数空间变异性的二维空间线性动态系统；在给定多道面波勘探的二维剪切波速数据和文献中勘探数据与岩土体参数的经验关系条件下，根据贝叶斯理论构建识别二维空间线性动态系统参数的贝叶斯方程；采用基于子集模拟的贝叶斯更新方法求解最优模型参数；采用向前递归算法和向后递归算法获得最可能的岩土体参数二维空间分布。合理地考虑剪切波速数据的二维空间相关性对岩土体参数空间变异性表征识别影响，充分利用二维空间勘探数据表征岩土体参数非稳态空间变异性，克服传统稳态随机场模型局限性，为提高岩土工程场地勘探准确性提供有效途径。

Description

基于多道面波勘探的岩土体参数二维空间变异性表征方法

技术领域

本发明涉及岩土工程勘察领域，具体涉及一种基于多道面波勘探的岩土体参数二维空间变异性表征方法。

背景技术

岩土体在其形成过程中，经受了各种复杂的地质作用，因而有着复杂的结构。因此，工程师在岩土工程设计中不仅要面对复杂的地层构造，并且还要兼顾到岩土类型、物理力学特性在空间上呈现出的变异性。岩土体参数(如土体容重、有效摩擦角)的准确性和可靠性将直接影响岩土工程设计结果。然而，工程师根据有限的勘察数据描述的岩土体参数和地层剖面存在不确定性。为了能够在岩土工程分析中合理地考虑岩土体参数空间变异性等不确定性，则需要定量表征勘察场地的岩土体参数空间变异性。

目前，广泛应用的岩土体参数空间变异性表征模型是随机场模型，随机场模型表征方法依赖原位测试手段(如静力触探试验、标准贯入试验)收集的一维勘察数据，难以根据勘察的二维剖面数据直接识别二维平稳随机场模型表征岩土体参数空间变异性。此外，以上所述的原位测试手段都属于侵入型试验，会对场地的岩土体产生扰动。工程师无法对场地的每个位置进行测试，对于岩土工程场地中没有任何勘察信息的空间位置，仅能根据有限的勘察信息进行推断和预测。尽管存在很多插值方法和拟合方法，但是仍然难以预测可能出现的未知土层，比如两个相邻钻孔或测点之间存在与勘探点处不同的岩土体类型或者局部软弱夹层透镜体等。因此，仅基于原位测试勘察信息无法得到全面的场地岩土体空间变异性。而非侵入式的多道面波勘探方法可提供二维空间数据，但是传统的稳态随机场模型又难以捕捉到岩土体参数的二维空间变化趋势。本发明基于工程地球物理探测收集的二维空间勘探数据，采用二维空间线性动态系统合理表征岩土体参数空间变异性。

基于上述情况，本发明提出了一种基于多道面波勘探的岩土体参数二维空间变异性表征方法，可有效解决以上问题。

发明内容

本发明的目的在于克服传统稳态随机场模型的局限性，提高岩土工程场地勘察的准确性。

为了达到上述目的，本发明提供一种基于多道面波勘探的岩土体参数二维空间变异性表征方法，包括如下步骤：

S1、利用多道面波勘探方法收集岩土体的空间二维剪切波速数据，用矩阵E＝(e_i,j)∈R^m×n表示；

S2、根据文献中的剪切波速e与岩土体参数f的经验关系公式推导二维空间线性动态系统中观测变量β和隐变量α之间的线性转换关系：

β＝Cα+b+ν (1)

式中：观测变量β是剪切波速数据e的函数，隐变量α是岩土体参数f的函数，C和b是已知的常数；v表示转换模型的误差，假设服从均值为零方差为Σ的高斯分布，即ν～N(0,Σ)，Σ即为根据隐变量α估计观测变量β时估计误差的方差；ψ_α,β＝{C,b,Σ}作为已知的模型参数；

S3、将勘探数据E转换成模型中的观测变量β用矩阵X_β＝(x_i,j)∈R^m×n表示，其中x_i,j表示矩阵中第i行第j列的元素且{i＝1,2,…,m；j＝1,2,…,n}，对应的隐变量α用矩阵Z_α＝(z_i,j)∈R^m×n表示，其中z_i,j表示矩阵中第i行第j列的元素且{i＝1,2,…,m；j＝1,2,…,n}，由公式(1)可得观测变量x_i,j与隐变量z_i,j之间的关系：

x_i,j＝Cz_i.j+b+v v～N(v|0,Σ)；1≤i≤m,1≤j≤n (2)

S4、建立二维空间线性动态系统，确定未知模型参数Φ_α,β＝{μ₀,V₀,A^V,Γ^V,A^H,Γ^H,A₁,A₂,Γ}；

S5、在给定二维观测变量数据X_β和模型参数ψ_α,β＝{C,b,Σ}条件下，基于贝叶斯理论确定未知模型参数的后验分布的贝叶斯公式：

式中：p(Φ_α,β|X_β,ψ_α,β)表示给定观测变量数据X_β和模型参数ψ_α,β条件下，未知模型参数的后验分布；p(X_β|Φ_α,β,ψ_α,β)表示给定所有模型参数条件下X_β的概率密度函数，称作似然函数；p(Φ_α,β|ψ_α,β)表示未知模型参数的先验分布；p(X_β|ψ_α,β)为归一化常数；

S6、根据先验信息确定未知模型参数Φ_α,β＝{μ₀,V₀,A^V,Γ^V,A^H,Γ^H,A₁,A₂,Γ}的先验分布p(Φ_α,β|ψ_α,β)；

S7、基于向前递归算法推导似然函数p(X_β|Φ_α,β,ψ_α,β)的计算公式；

S8、采用基于子集模拟的贝叶斯更新方法求解贝叶斯方程即公式(3)，得到未知模型参数Φ_α,β的后验样本，根据参数的后验样本确定最可能的模型参数Φ_α,β ^*＝{μ₀ ^*,V₀ ^*,A^V*,Γ^V*,A^H*,Γ^H*,A₁ ^*,A₂ ^*,Γ^*}；

S9、基于向后递归算法计算给定ψ_α,β、Φ_α,β ^*和X_β条件下的隐变量Z_α的边缘后验高斯分布γ(z_i,j)＝p(z_i,j|X_β,ψ_α,β,Φ_α,β ^*)的分布参数

和

为隐变量边缘分布的方差，

为隐变量边缘分布的均值，

组成最可能的隐变量矩阵Z_α ^*；

S10、根据步骤S2中隐变量α与岩土体参数f的函数关系，得到最可能的岩土体参数的二维空间分布矩阵F^*＝(f_i,j ^*)∈R^m×n。

优选的，步骤S4中建立的二维空间线性动态系统，具体为：假设二维空间中隐变量Z_α的每一个位置的状态z_i,j与水平和垂直方向上相邻位置的状态具有如下线性关系：

z_1，1＝μ₀+u u～N(u|0,V₀) (4)

式中：z_1,1表示第一行第一列隐变量，μ₀表示隐变量的初始估计值，u表示估计误差，假设误差u服从均值为0方差为V₀的高斯分布；

z_1,j＝A^Hz_1,j-1+w^H w^H～N(w^H|0,Γ^H)；1＜j≤n (5)

式中：z_1,j表示第一行第j列隐变量，A^H表示第一行相邻隐变量线性关系系数，w^H表示根据z_1,j-1估计z_1,j的估计误差，假设误差w^H服从均值为0方差为Γ^H的高斯分布；

z_i,1＝A^Vz_i-1,1+w^V w^V～N(w^V|0,Γ^V)；1＜i≤m (6)

式中：z_i,1表示第i行第一列隐变量，A^V表示第一列相邻隐变量线性关系系数，w^V表示根据z_i-1,1估计z_i,1的估计误差，假设误差w^V服从均值为0方差为Γ^V的高斯分布；

z_i,j＝A₁z_i-1,j+A₂z_i,j-1+w w～N(w|0,Γ)；1＜i≤m,1＜j≤n (7)

式中：z_i,j表示第i行第j列隐变量，A₁和A₂分别表示竖直方向和水平方向的线性关系系数，w表示根据z_i-1,j和z_i,j-1估计z_i,j的估计误差，假设误差w服从均值为0方差为Γ的高斯分布；

由公式(4)～(7)可以推导出隐变量的条件高斯分布：

p(z_1，1|Φ_α,β)＝N(z_1，1|μ₀,V₀) (8)

p(z_1,j|z_1,j-1,Φ_α,β)＝N(z_1,j|A^Hz_1,j-1,Γ^H) (9)

p(z_i,1|z_i-1,1,Φ_α,β)＝N(z_i,1|A^Vz_i-1,1,Γ^V) (10)

p(z_i,j|z_i-1,j,z_i,j-1,Φ_α,β)＝N(z_i,j|A₁z_i-1,j+A₂z_i,j-1,Γ) (11)

由公式(2)可以推导出观测变量的条件高斯分布：

p(x_i,j|z_i,j,ψ_α,β)＝N(x_i,j|Cz_i,j+b,Σ) (12)

由公式(4)～(7)建立的二维空间线性动态系统可以确定未知模型参数为Φ_α,β＝{μ₀,V₀,A^V,Γ^V,A^H,Γ^H,A₁,A₂,Γ}。

优选的，步骤S7中基于向前递归算法推导似然函数p(X_β|Φ_α,β,ψ_α,β)的计算公式，具体步骤为：

S71、由条件概率可知p(X_β|Φ_α,β,ψ_α,β)可以分解成多个条件边缘分布乘积的形式：

S72、根据条件概率理论，由公式(8)和(12)的条件高斯分布可以推导出x_1,1条件高斯分布p(x_1,1|Φ_α,β,ψ_α,β)为：

p(x_1,1|Φ_α,β,ψ_α,β)＝N(x_1,1|Cμ₀+b,Σ+CV₀C^T) (14)

S73、根据条件概率理论，由公式(9)和(12)可以推导条件高斯分布p(x_1,j|x_1,1,x_1,2:j-1,Φ_α,β,ψ_α,β)为：

式中：

S74、根据条件概率理论，由公式(10)和(12)可以推导条件高斯分布p(x_i,1|x_1,1,x_1,2:n,x_2:i-1,1,Φ_α,β,ψ_α,β)为：

式中：

S75、根据条件概率理论，由公式(11)和(12)可以推导条件高斯分布p(x_i,j|x_1,1,x_1,2:n,x_2:m,1,x_2:i-1,2:n,x_i,2:j-1,Φ_α,β,ψ_α,β)为：

式中：

S76、由公式(13)～(17)可以推导出似然函数计算公式为：

优选的，步骤S9中基于向后递归算法计算给定ψ_α,β、Φ_α,β ^*和X_β条件下最可能的隐变量矩阵Z_α ^*＝(z_i,j ^*)∈R^m×n，由于隐变量Z_α的边缘后验分布p(z_i,j|X_β,ψ_α,β,Φ_α,β ^*)服从高斯分布，因此，隐变量边缘后验分布的均值

即为最可能隐变量z_i,j ^*。

优选的，计算

的具体步骤为：

S91、由步骤S7中向前递归算法可以得隐变量矩阵Z_α中第i行第j列隐变量z_i,j的条件边缘分布步骤如下所示：

S911、当i＝1且j＝1时：

式中：

S912、当i＝1且1<j≤n时：

式中：

S913、当1<i≤m且j＝1时：

式中：

S914、当1<i≤m且1<j≤n时：

式中：

S92、由公式(22)可知第m行第n列隐变量z_m,n的后验高斯分布为：

式中：

S93、结合公式(22)由向后递归算法可以得出z_m,n-1和z_m-1,n的联合后验高斯分布为：

式中：

S94、当1≤i≤m-2且j＝n时，z_i,n的后验高斯分布为：

式中：

S95、当i＝m且1≤j≤n-2时，z_m,j的后验高斯分布为：

式中：

S96、当2≤i≤m-1且2≤j≤n-1时，z_i,j的后验高斯分布为：

式中：

S97、当i＝1且2≤j≤n-1时，z_1,j的后验高斯分布为：

式中：

S98、当2≤i≤m-1且j＝1时，z_i,1的后验高斯分布为：

式中：

S99、当i＝1且j＝1时，z_1,1的后验高斯分布为：

式中：

本发明与现有技术相比，具有以下优点及有益效果：

与现有方法相比，基于多道面波勘探的岩土体参数二维空间变异性表征方法，为二维空间数据分析和反演提供了新思路，充分利用多道面波勘探的二维空间剪切波速数据表征岩土体参数空间异性，克服了传统稳态随机场模型的局限性，为提高岩土工程场地勘探的准确性提供了有效途径。

附图说明

图1是本发明的流程框图。

图2是二维空间线性动态系统示意图。

图3是二维空间线性动态系统与模型参数示意图。

图4是向前递归算法计算似然函数流程图。

图5是向前递归算法计算隐变量条件边缘分布流程图。

图6是向后递归算法计算最可能的隐变量矩阵Z_α ^*流程图。

图7是多道面波勘探和数据采集示意图。

图8是一套天然边坡的剪切波速二维空间数据。

图9是文献中剪切波速与天然密度的经验关系。

图10是观测变量数据。

图11是最可能的隐变量数据。

具体实施方式

下面结合附图和具体实例对本发明作详细说明。

如图1所示，本发明所述的基于多道面波勘探的岩土体参数二维空间变异性表征方法，具体实施过程为：

(1)、利用多道面波勘探方法收集岩土体的空间二维剪切波速数据，用矩阵E＝(e_i,j)∈R^m×n表示。

(2)、根据历史文献中的剪切波速e与岩土体参数f的经验关系公式推导二维空间线性动态系统中观测变量β和隐变量α之间的线性转换关系：

β＝Cα+b+ν (1)

式中：观测变量β是剪切波速e的函数，隐变量α是岩土体参数f的函数，C和b是已知的常数；v表示转换模型的误差，假设服从均值为零方差为Σ的高斯分布，即ν～N(0,Σ)，Σ即为根据隐变量α估计观测变量β时估计误差的方差。ψ_α,β＝{C,b,Σ}作为已知的模型参数。

(3)、将勘探数据E转换成模型中的观测变量β用矩阵X_β＝(x_i,j)∈R^m×n表示，其中x_i,j表示矩阵中第i行第j列的元素且{i＝1,2,…,m；j＝1,2,…,n}，对应的隐变量α用矩阵Z_α＝(z_i,j)∈R^m×n表示，其中z_i,j表示矩阵中第i行第j列的元素且{i＝1,2,…,m；j＝1,2,…,n}，由公式(1)可得观测变量x_i,j与隐变量z_i,j之间的关系：

x_i,j＝Cz_i.j+b+v v～N(v|0,Σ)；1≤i≤m,1≤j≤n (2)

(4)、建立二维空间线性动态系统，确定未知模型参数Φ_α,β＝{μ₀,V₀,A^V,Γ^V,A^H,Γ^H,A₁,A₂,Γ}。

(5)、在给定二维观测变量数据X_β和模型参数ψ_α,β＝{C,b,Σ}条件下，基于贝叶斯理论确

定未知模型参数的后验分布的贝叶斯公式：

式中：p(Φ_α,β|X_β,ψ_α,β)表示给定观测变量数据X_β和模型参数ψ_α,β条件下，未知模型参数的后验分布；p(X_β|Φ_α,β,ψ_α,β)表示给定所有模型参数条件下X_β的概率密度函数，称作似然函数；p(Φ_α,β|ψ_α,β)表示未知模型参数的先验分布；p(X_β|ψ_α,β)为归一化常数。

(6)、根据先验信息确定未知模型参数Φ_α,β＝{μ₀,V₀,A^V,Γ^V,A^H,Γ^H,A₁,A₂,Γ}的先验分布p(Φ_α,β|ψ_α,β)。

(7)、基于向前递归算法推导似然函数p(X_β|Φ_α,β,ψ_α,β)的计算公式。

(8)、采用基于子集模拟的贝叶斯更新方法求解贝叶斯方程即公式(3)，得到未知模型参数Φ_α,β的后验样本，根据参数的后验样本确定最可能的模型参数Φ_α,β ^*＝{μ₀ ^*,V₀ ^*,A^V*,Γ^V*,A^H*,Γ^H*,A₁ ^*,A₂ ^*,Γ^*}。

(9)、基于向后递归算法计算给定ψ_α,β、Φ_α,β ^*和X_β条件下的隐变量Z_α的边缘后验高斯分布γ(z_i,j)＝p(z_i,j|X_β,ψ_α,β,Φ_α,β ^*)的分布参数

和

为隐变量边缘分布的方差，

为隐变量边缘分布的均值，

组成最可能的隐变量矩阵Z_α ^*。

(10)、根据步骤(2)中隐变量α与岩土体参数f的函数关系，得到最可能的岩土体参数的二维空间分布矩阵F^*＝(f_i,j ^*)∈R^m×n。

如图2、3所示，上述步骤(4)中建立的二维空间线性动态系统，具体为：假设二维空间中隐变量Z_α的每一个位置的状态z_i,j与水平和垂直方向上相邻位置的状态具有如下线性关系：

z_1，1＝μ₀+u u～N(u|0,V₀) (4)

式中：z_1,1表示第一行第一列隐变量，μ₀表示隐变量的初始估计值，u表示估计误差，假设误差u服从均值为0方差为V₀的高斯分布；

z_1,j＝A^Hz_1,j-1+w^H w^H～N(w^H|0,Γ^H)；1＜j≤n (5)

式中：z_1,j表示第一行第j列隐变量，A^H表示第一行相邻隐变量线性关系系数，w^H表示根据z_1,j-1估计z_1,j的估计误差，假设误差w^H服从均值为0方差为Γ^H的高斯分布；

z_i,1＝A^Vz_i-1,1+w^V w^V～N(w^V|0,Γ^V)；1＜i≤m (6)

式中：z_i,1表示第i行第一列隐变量，A^V表示第一列相邻隐变量线性关系系数，w^V表示根据z_i-1,1估计z_i,1的估计误差，假设误差w^V服从均值为0方差为Γ^V的高斯分布；

z_i,j＝A₁z_i-1,j+A₂z_i,j-1+w w～N(w|0,Γ)；1＜i≤m,1＜j≤n (7)

式中：z_i,j表示第i行第j列隐变量，A₁和A₂分别表示竖直方向和水平方向的线性关系系数，w表示根据z_i-1,j和z_i,j-1估计z_i,j的估计误差，假设误差w服从均值为0方差为Γ的高斯分布；

由公式(4)～(7)可以推导出隐变量的条件高斯分布：

p(z_1，1|Φ_α,β)＝N(z_1，1|μ₀,V₀) (8)

p(z_1,j|z_1,j-1,Φ_α,β)＝N(z_1,j|A^Hz_1,j-1,Γ^H) (9)

p(z_i,1|z_i-1,1,Φ_α,β)＝N(z_i,1|A^Vz_i-1,1,Γ^V) (10)

p(z_i,j|z_i-1,j,z_i,j-1,Φ_α,β)＝N(z_i,j|A₁z_i-1,j+A₂z_i,j-1,Γ) (11)

由公式(2)可以推导出观测变量的条件高斯分布：

p(x_i,j|z_i,j,ψ_α,β)＝N(x_i,j|Cz_i,j+b,Σ) (12)

由公式(4)～(7)建立的二维空间线性动态系统可以确定未知模型参数为Φ_α,β＝{μ₀,V₀,A^V,Γ^V,A^H,Γ^H,A₁,A₂,Γ}。

如图4所示，上述步骤(7)中基于向前递归算法推导似然函数p(X_β|Φ_α,β,ψ_α,β)的计算公式，具体步骤为：

(7.1)由条件概率可知p(X_β|Φ_α,β,ψ_α,β)可以分解成多个条件边缘分布乘积的形式：

(7.2)根据条件概率理论，由公式(8)和(12)的条件高斯分布可以推导出x_1,1条件高斯分布p(x_1,1|Φ_α,β,ψ_α,β)为：

p(x_1,1|Φ_α,β,ψ_α,β)＝N(x_1,1|Cμ₀+b,Σ+CV₀C^T) (14)

(7.3)根据条件概率理论，由公式(9)和(12)可以推导条件高斯分布p(x_1,j|x_1,1,x_1,2:j-1,Φ_α,_β,ψ_α,β)为：

式中：

(7.4)根据条件概率理论，由公式(10)和(12)可以推导条件高斯分布p(x_i,1|x_1,1,x_1,2:n,x₂:_i-1,1,Φ_α,β,ψ_α,β)为：

式中：

(7.5)根据条件概率理论，由公式(11)和(12)可以推导条件高斯分布p(x_i,j|x_1,1,x_1,2:n,x_2:m,1,x_2:i-1,2:n,x_i,2:j-1,Φ_α,β,ψ_α,β)为：

式中：

(7.6)由公式(13)～(17)可以推导出似然函数计算公式为：

如图6所示，上述步骤(9)中基于向后递归算法计算给定ψ_α,β、Φ_α,β ^*和X_β条件下最可能的隐变量矩阵Z_α ^*＝(z_i,j ^*)∈R^m×n，由于隐变量Z_α的边缘后验分布p(z_i,j|X_β,ψ_α,β,Φ_α,β ^*)服从高斯分布，因此，隐变量边缘后验分布的均值

即为最可能隐变量z_i,j ^*。计算

的具体步骤为：

(9.1)如图5所示，由步骤(7)中向前递归算法可以得隐变量矩阵Z_α中第i行第j列隐变量z_i,j的条件边缘分布步骤如下所示：

(9.1.1)当i＝1且j＝1时：

式中：

(9.1.2)当i＝1且1<j≤n时：

式中：

(9.1.3)当1<i≤m且j＝1时：

式中：

(9.1.4)当1<i≤m且1<j≤n时：

式中：

(9.2)由公式(22)可知第m行第n列隐变量z_m,n的后验高斯分布为：

式中：

(9.3)结合公式(22)由向后递归算法可以得出z_m,n-1和z_m-1,n的联合后验高斯分布为：

式中：

(9.4)当1≤i≤m-2且j＝n时，z_i,n的后验高斯分布为：

式中：

(9.5)当i＝m且1≤j≤n-2时，z_m,j的后验高斯分布为：

式中：

(9.6)当2≤i≤m-1且2≤j≤n-1时，z_i,j的后验高斯分布为：

式中：

(9.7)当i＝1且2≤j≤n-1时，z_1,j的后验高斯分布为：

式中：

(9.8)当2≤i≤m-1且j＝1时，z_i,1的后验高斯分布为：

式中：

(9.9)当i＝1且j＝1时，z_1,1的后验高斯分布为：

式中：

Q_1,1＝((Q_i′_,j)^-1+(A^H*)^T(Γ^H*)^-1A^H*)^-1；

下面结合多道面波方法勘探的剪切波速数据反演岩土天然密度进行实例说明。

本发明的具体实施过程如下：

(1)利用多道面波勘探方法收集岩土体的空间二维剪切波速数据：

图7为多道面波勘探和数据采集示意图，如图所示一共设置了12个检波器，每个检波器之间等间隔，即道间距为2m，偏移距为10m，表示震源和第一个检波器之间的距离，移动步距为4m，表示下一次试验震源位置和上一次试验震源位置的距离，下一次试验时将震源和检波器一起平移不改变偏移距和道间距。每一次试验测得第6个和第7个检波器之间中间位置地下一定埋深的剪切波速，从坡顶到坡底顺次测量得到等间隔的一维剪切波速数据，然后合并得到沿边坡A-B方向断面的二维剪切波速数据。

图8为多道面波方法勘探的一个天然边坡场地的剪切波速V_s(m/s)二维空间数据E∈R^40×350，探测水平方向总长度为350m，数据间隔为1m，竖直方向总探测深度为20m，数据间隔为0.5m。

(2)如图9所示，根据文献中的勘探数据V_s(m/s)与岩土体天然密度ρ_d(g/cm³)的经验关系：ρ_d＝0.277+0.648log₁₀(V_s)，根据V_s估计ρ_d的估计误差的标准差为0.126，由此推导出观测变量β＝log₁₀(V_s)和隐变量α＝ρ_d之间的线性转换关系：

根据公式(31)可知二维空间线性动态系统中参数ψ_α,β＝{C,b,Σ}的取值，即C＝1/0.648＝1.54；b＝-0.277/0.648＝-0.43；Σ＝(0.126/0.648)²＝0.04。图10为由勘探数据E计算得到的观测变量数据X_β。

(3)确定二维空间线性动态系统未知参数Φ_α,β＝{μ₀,V₀,A^V,Γ^V,A^H,Γ^H,A₁,A₂,Γ}的先验分布：在模型参数信息不足时，先验分布通常在合理的范围内取为均匀分布。μ₀和V₀分别表示初始位置隐变量密度ρ_d的均值和方差，可以根据图9中纵坐标密度的实测数据的统计范围来确定，因此，确定均值μ₀的最大范围为0.6～2.6(g/cm³)，方差的最大值假设为ρ_d服从均匀分布U[0.6～2.6]时的方差，即V_0,max＝0.33。由于相邻土体参数相关性较强性质接近，假设关系系数A^V和A^H不会超过2，假设根据相邻状态估计的隐变量的误差方差Γ^V、Γ^H和Γ不会超过1。对于除第一行和第一列以外的位置的隐变量受到相邻两个方向的同时影响，假设两个方向所占的权重一样，z_i,j最大不会超过(z_i,j-1+z_i-1,j)，即A₁和A₂均不会超过1。并假设每一个参数的先验分布为独立的均匀分布。综上所述，可得未知模型参数的先验分布范围如表1所示。

表1未知模型参数Φ_α,β的先验范围

(4)建立未知模型参数后验分布的贝叶斯方程：

(5)采用基于子集模拟的贝叶斯更新方法求解公式(32)贝叶斯方程，识别最可能的未知模型参数，结果如表2所示。

表2最可能的模型参数Φ_α,β ^*

(6)采用先前递归算法和向后递归算法计算给定ψ_α,β、Φ_α,β ^*和X_β条件下的最可能的隐变量数据Z_α ^*，结果如图11所示。本算例中隐变量即为天然密度，因此，图11同时表示最可能天然密度二维空间分布。

通过一套二维剪切波速勘探数据反演土体天然密度二维空间分布，实例表明本发明所提出的方法合理地考虑了观测数据和岩土体参数之间的转换模型误差以及岩土体参数空间相关性，反演的土体密度二维剖面表征了天然密度的空间变异性以及地下土体密度的空间变化趋势，更加符合工程实际。

Claims (5)

1.基于多道面波勘探的岩土体参数二维空间变异性表征方法，其特征在于，包括如下步骤：

S1、利用多道面波勘探方法收集岩土体的空间二维剪切波速数据，用矩阵E＝(e_i,j)∈R^{m ×n}表示；

S2、根据文献中的剪切波速e与岩土体参数f的经验关系公式推导二维空间线性动态系统中观测变量β和隐变量α之间的线性转换关系：

β＝Cα+b+ν (1)

式中：观测变量β是剪切波速数据e的函数，隐变量α是岩土体参数f的函数，C和b是已知的常数；v表示转换模型的误差，假设服从均值为零方差为Σ的高斯分布，即ν～N(0,Σ)，Σ即为根据隐变量α估计观测变量β时估计误差的方差；ψ_α,β＝{C,b,Σ}作为已知的模型参数；

S3、将勘探数据E转换成模型中的观测变量β用矩阵X_β＝(x_i,j)∈R^m×n表示，其中x_i,j表示X_β矩阵中第i行第j列的元素，对应的隐变量α用矩阵Z_α＝(z_i,j)∈R^m×n表示，其中z_i,j表示Z_α矩阵中第i行第j列的元素，由公式(1)可得观测变量x_i,j与隐变量z_i,j之间的关系：

x_i,j＝Cz_i.j+b+v v～N(v|0,Σ)；1≤i≤m,1≤j≤n (2)

S4、建立二维空间线性动态系统，确定未知模型参数Φ_α,β＝{μ₀,V₀,A^V,Γ^V,A^H,Γ^H,A₁,A₂,Γ}，其中，μ₀和V₀分别表示初始隐变量z_1,1的均值和方差；A^V和Γ^V分别表示第一列相邻隐变量z_i-1,1和z_i,1的线性关系系数和线性关系误差的方差；A^H和Γ^H分别表示第一行相邻隐变量z_1,j-1和z_1,j的线性关系系数和线性关系误差的方差；A₁和A₂分别表示隐变量z_i,j、z_i-1,j和z_i,j-1的线性关系系数；Γ表示z_i,j、z_i-1,j和z_i,j-1之间线性关系误差的方差；

S5、在给定二维观测变量数据X_β和模型参数ψ_α,β＝{C,b,Σ}条件下，基于贝叶斯理论确定未知模型参数的后验分布的贝叶斯公式：

式中：p(Φ_α,β|X_β,ψ_α,β)表示给定观测变量数据X_β和模型参数ψ_α,β条件下，未知模型参数的后验分布；p(X_β|Φ_α,β,ψ_α,β)表示给定所有模型参数条件下X_β的概率密度函数，称作似然函数；p(Φ_α,β|ψ_α,β)表示未知模型参数的先验分布；p(X_β|ψ_α,β)为归一化常数；

S6、根据先验信息确定未知模型参数Φ_α,β＝{μ₀,V₀,A^V,Γ^V,A^H,Γ^H,A₁,A₂,Γ}的先验分布p(Φ_α,β|ψ_α,β)；

S7、基于向前递归算法推导似然函数p(X_β|Φ_α,β,ψ_α,β)的计算公式；

S8、采用基于子集模拟的贝叶斯更新方法求解贝叶斯方程即公式(3)，得到未知模型参数Φ_α,β的后验样本，根据参数的后验样本确定最可能的模型参数Φ_α,β ^*＝{μ₀ ^*,V₀ ^*,A^V*,Γ^V*,A^H*,Γ^H*,A₁ ^*,A₂ ^*,Γ^*}；

S9、基于向后递归算法计算给定ψ_α,β、Φ_α,β ^*和X_β条件下的隐变量Z_α的边缘后验高斯分布γ(z_i,j)＝p(z_i,j|X_β,ψ_α,β,Φ_α,β ^*)的分布参数

和

为隐变量边缘分布的方差，

为隐变量边缘分布的均值，

组成最可能的隐变量矩阵Z_α ^*；

S10、根据步骤S2中隐变量α与岩土体参数f的函数关系，得到最可能的岩土体参数的二维空间分布矩阵F^*＝(f_i,j ^*)∈R^m×n。

2.根据权利要求1所述的基于多道面波勘探的岩土体参数二维空间变异性表征方法，其特征在于：步骤S4中建立的二维空间线性动态系统，具体为：假设二维空间中隐变量Z_α的每一个位置的状态z_i,j与水平和垂直方向上相邻位置的状态具有如下线性关系：

z_1，1＝μ₀+u u～N(u|0,V₀) (4)

式中：z_1,1表示第一行第一列隐变量，μ₀表示隐变量的初始估计值，u表示估计误差，假设误差u服从均值为0方差为V₀的高斯分布；

z_1,j＝A^Hz_1,j-1+w^H w^H～N(w^H|0,Γ^H)；1＜j≤n (5)

式中：z_1,j表示第一行第j列隐变量，A^H表示第一行相邻隐变量线性关系系数，w^H表示根据z_1,j-1估计z_1,j的估计误差，假设误差w^H服从均值为0方差为Γ^H的高斯分布；

z_i,1＝A^Vz_i-1,1+w^V w^V～N(w^V|0,Γ^V)；1＜i≤m (6)

式中：z_i,1表示第i行第一列隐变量，A^V表示第一列相邻隐变量线性关系系数，w^V表示根据z_i-1,1估计z_i,1的估计误差，假设误差w^V服从均值为0方差为Γ^V的高斯分布；

z_i,j＝A₁z_i-1,j+A₂z_i,j-1+w w～N(w|0,Γ)；1＜i≤m,1＜j≤n (7)

式中：z_i,j表示第i行第j列隐变量，A₁和A₂分别表示竖直方向和水平方向的线性关系系数，w表示根据z_i-1,j和z_i,j-1估计z_i,j的估计误差，假设误差w服从均值为0方差为Γ的高斯分布；

由公式(4)～(7)可以推导出隐变量的条件高斯分布：

p(z_1，1|Φ_α,β)＝N(z_1，1|μ₀,V₀) (8)

p(z_1,j|z_1,j-1,Φ_α,β)＝N(z_1,j|A^Hz_1,j-1,Γ^H) (9)

p(z_i,1|z_i-1,1,Φ_α,β)＝N(z_i,1|A^Vz_i-1,1,Γ^V) (10)

p(z_i,j|z_i-1,j,z_i,j-1,Φ_α,β)＝N(z_i,j|A₁z_i-1,j+A₂z_i,j-1,Γ) (11)

由公式(2)可以推导出观测变量的条件高斯分布：

p(x_i,j|z_i,j,ψ_α,β)＝N(x_i,j|Cz_i,j+b,Σ) (12)

由公式(4)～(7)建立的二维空间线性动态系统可以确定未知模型参数为Φ_α,β＝{μ₀,V₀,A^V,Γ^V,A^H,Γ^H,A₁,A₂,Γ}，其中，μ₀和V₀分别表示初始隐变量z_1,1的均值和方差；A^V和Γ^V分别表示第一列相邻隐变量z_i-1,1和z_i,1的线性关系系数和线性关系误差的方差；A^H和Γ^H分别表示第一行相邻隐变量z_1,j-1和z_1,j的线性关系系数和线性关系误差的方差；A₁和A₂分别表示隐变量z_i,j、z_i-1,j和z_i,j-1的线性关系系数；Γ表示z_i,j、z_i-1,j和z_i,j-1之间线性关系误差的方差。

3.根据权利要求2所述的基于多道面波勘探的岩土体参数二维空间变异性表征方法，其特征在于：步骤S7中基于向前递归算法推导似然函数p(X_β|Φ_α,β,ψ_α,β)的计算公式，具体步骤为：

S71、由条件概率可知p(X_β|Φ_α,β,ψ_α,β)可以分解成多个条件边缘分布乘积的形式：

S72、根据条件概率理论，由公式(8)和(12)的条件高斯分布可以推导出x_1,1条件高斯分布p(x_1,1|Φ_α,β,ψ_α,β)为：

p(x_1,1|Φ_α,β,ψ_α,β)＝N(x_1,1|Cμ₀+b,Σ+CV₀C^T) (14)

S73、根据条件概率理论，由公式(9)和(12)可以推导条件高斯分布p(x_1,j|x_1,1,x_1,2:j-1,Φ_α,β,ψ_α,β)为：

式中：

S74、根据条件概率理论，由公式(10)和(12)可以推导条件高斯分布p(x_i,1|x_1,1,x_1,2:n,x_2:i-1,1,Φ_α,β,ψ_α,β)为：

式中：

S75、根据条件概率理论，由公式(11)和(12)可以推导条件高斯分布p(x_i,j|x_1,1,x_1,2:n,x_2:m,1,x_2:i-1,2:n,x_i,2:j-1,Φ_α,β,ψ_α,β)为：

式中：

S76、由公式(13)～(17)可以推导出似然函数计算公式为：

4.根据权利要求1所述的基于多道面波勘探的岩土体参数二维空间变异性表征方法，其特征在于：步骤S9中基于向后递归算法计算给定ψ_α,β、Φ_α,β ^*和X_β条件下最可能的隐变量矩阵Z_α ^*＝(z_i,j ^*)∈R^m×n，由于隐变量Z_α的边缘后验分布p(z_i,j|X_β,ψ_α,β,Φ_α,β ^*)服从高斯分布，因此，隐变量边缘后验分布的均值

即为最可能隐变量z_i,j ^*。

5.根据权利要求4所述的基于多道面波勘探的岩土体参数二维空间变异性表征方法，其特征在于，计算

的具体步骤为：

S91、由步骤S7中向前递归算法可以得隐变量矩阵Z_α中第i行第j列隐变量z_i,j的条件边缘分布步骤如下所示：

S911、当i＝1且j＝1时：

式中：

S912、当i＝1且1<j≤n时：

式中：

S913、当1<i≤m且j＝1时：

式中：

S914、当1<i≤m且1<j≤n时：

式中：

S92、由公式(22)可知第m行第n列隐变量z_m,n的后验高斯分布为：

式中：

S93、结合公式(22)由向后递归算法可以得出z_m,n-1和z_m-1,n的联合后验高斯分布为：

式中：

S94、当1≤i≤m-2且j＝n时，z_i,n的后验高斯分布为：

式中：

S95、当i＝m且1≤j≤n-2时，z_m,j的后验高斯分布为：

式中：

S96、当2≤i≤m-1且2≤j≤n-1时，z_i,j的后验高斯分布为：

式中：

S97、当i＝1且2≤j≤n-1时，z_1,j的后验高斯分布为：

式中：

S98、当2≤i≤m-1且j＝1时，z_i,1的后验高斯分布为：

式中：

S99、当i＝1且j＝1时，z_1,1的后验高斯分布为：

式中：

Q_1,1＝((Q′_i,j)^-1+(A^H*)^T(Γ^H*)^-1A^H*)^-1；

Images