基于贝叶斯理论面波频散曲线随机反演

doi:10.11720/wtyht.2021.1100

基于贝叶斯理论面波频散曲线随机反演

刘辉^,¹, 李静^,¹^,², 曾昭发¹, 王天琪¹

1.吉林大学地球探测科学与技术学院,吉林长春 130021

2.地球信息探测仪器教育部重点实验室 (吉林大学),吉林长春 130026

Stochastic inversion of surface wave dispersion curves based on Bayesian theory

LIU Hui^,¹, LI Jing^,¹^,², ZENG Zhao-Fa¹, WANG Tian-Qi¹

1. College of Geo-exploration Science and Technology,Jilin University,Changchun 130021,China

2. Key Laboratory of Geophysical Exploration Equipment,Ministry of Education (Jilin University),Changchun 130026,China

通讯作者: 李静(1988-),男,博士,教授,博士生导师,主要从事应用地球物理新方法技术研究工作。Email:ljwy1209@163.com

责任编辑: 叶佩

收稿日期: 2020-03-9 修回日期: 2021-01-13

基金资助:

国家自然科学基金项目.  41874134
吉林省优秀青年基金.  20190103142JH
吉林省自然科学基金.  20200201216JC
中国科协第五届青年托举人才项目.  2019QNRC001

Received: 2020-03-9 Revised: 2021-01-13

作者简介 About authors

刘辉(1996-),男,硕士研究生,主要从事面波频散曲线反演方法研究工作。Email: 177344303@qq.com

摘要

面波频散曲线反演是获得地下横波速度结构的重要地球物理方法。常规基于迭代最小二乘等线性反演方法依赖于初始模型,且存在多极值、容易陷入局部最小、反演精度低等问题。基于贝叶斯理论的随机反演方法是一种可以融合先验信息的非线性反演方法,该方法无需人为给定初始模型,仅利用先验信息对模型进行随机采样,根据概率分布筛选接受合适的后验概率密度估计结果,可达到对细节信息的准确估计。本文针对瑞利面波频散曲线,提出了基于GPR数据先验资料约束的贝叶斯马尔科夫蒙特卡洛(MCMC)随机反演方法,通过随机改变模型参数并计算其频散曲线与实际频散曲线的似然函数来选择是否接受新的模型参数,不断重复此过程,最终得到与实际频散曲线拟合效果最佳的最优解以及横波速度解的后验概率密度分布。通过理论模型以及实际数据反演测试,验证了该方法与常规无约束的随机反演相比,可以有效地提高反演速度和反演精度。

关键词： 贝叶斯理论 ; 马尔科夫链随机反演 ; 瑞利波频散曲线 ; 约束反演

Abstract

Surface wave dispersion curve inversion is an important geophysical method for obtaining the velocity and thickness distribution of underground shear wave.Conventional linear inversion methods,such as iterative least squares,relying on the initial model and multiple solution,are easy to fall into local minimum and low inversion accuracy.The stochastic inversion method based on Bayesian theory is a nonlinear inversion method which can integrate prior information.This method does not need initial model,only uses prior information to sample the model randomly,and selects and accepts the appropriate inversion model according to the probability distribution.It achieves the accurate estimation of the detail information.In this paper,the authors present a Bayesian Markov Monte Carlo (MCMC) stochastic inversion method based on GPR data constraints to invert the Rayleigh-waves dispersion curve.In the inversion process,by randomly changing the model parameters and calculating the likelihood function of the dispersion curve and the actual dispersion curve,researchers can choose whether to accept the new model parameters,repeat this process continuously,and finally get the best fitting result with the actual dispersion curve and the posterior probability density distribution of the VS solution.The typical numerical model test and field seismic data demonstrate that,compared with the conventional unconstrained stochastic inversion,the proposed method can effectively reduce the multiple solution and improve the efficiency and accuracy.

Keywords： Bayesian theory ; Monte Carlo Markov Chain (MCMC) stochastic inversion ; Rayleigh wave dispersion curve ; constraint inversion

PDF (4903KB) 元数据多维度评价相关文章导出 EndNote| Ris| Bibtex 收藏本文

本文引用格式

刘辉, 李静, 曾昭发, 王天琪. 基于贝叶斯理论面波频散曲线随机反演. 物探与化探[J], 2021, 45(4): 951-960 doi:10.11720/wtyht.2021.1100

LIU Hui, LI Jing, ZENG Zhao-Fa, WANG Tian-Qi. Stochastic inversion of surface wave dispersion curves based on Bayesian theory. Geophysical and Geochemical Exploration[J], 2021, 45(4): 951-960 doi:10.11720/wtyht.2021.1100

0 引言

瑞利面波是一种由纵波与横波干涉形成沿自由表面传播的地震波。1887年,英国学者Rayleigh首先发现并证明均匀半空间中瑞利面波的存在^[1,2]。面波频散曲线反演是获取浅地表地下空间、深部岩石圈和地球构造横波速度的重要方法技术^[3,4,5]。Park等^[6]提出的多道面波分析方法(MASW)广泛应用于浅地表面波数据采集处理;Xia等^[7]采取奇异值分解和L-M方法,实现了瑞利波基阶频散曲线反演;Song等^[8]使用Occam算法进行一维频散曲线反演后获得了一个含低速层的拟二维横波速度剖面;鲁来玉等^[9]和罗银河等^[10]分别采用遗传算法和最小二乘方法实现了不同尺度多阶面波频散联合反演,可有效提高单纯只依靠基阶面波的反演精度。

常规基于迭代的最小二乘线性反演方法依赖于初始模型且存在多极值、容易陷入局部最小、反演精度低等问题。基于全局最优的非线性反演方法由于不依赖于初始模型,已广泛应用于地球物理反演。Beaty等^[11]利用模拟退火法在一定概率下可以接受劣解的搜索方式实现了多阶瑞利波联合反演;宋先海等将模式识别算法^[12]和粒子群算法^[13]应用到瑞利波频散曲线反演中;蔡伟等^[14]将萤火虫和蝙蝠群智能算法应用于瑞利波频散曲线反演;于东凯等分别利用改进蜂群算法^[15]和蚱蜢算法^[16]实现了瑞利波频散曲线反演。在天然地震被动源面波成像研究中,Aki^[17]提出利用背景噪声研究地下结构,从台阵信号中提取瑞利波频散信息的空间自相关法(SPAC);徐佩芬等^[18]利用SPAC法从微动信号中提取频散曲线并利用遗传算法进行反演获得S波速度结构;张宝龙等^[19]利用扩展空间自相关法(ESPAC)提取频散曲线结合面波层析成像技术获得了地壳浅层三维S波速度结构。虽然传统算法如遗传算法具有对初始模型要求宽松、无需求梯度、易于实现并行化等优点,较适合于解决非线性、多参数的瑞利波频散曲线反演问题,但也存在计算量大、局部搜索能力差和收敛稳定性差等缺陷,使其在频散曲线反演中受到一定的限制。

针对上述算法中存在的问题,结合前人研究成果,本文提出一种基于贝叶斯马尔科夫蒙特卡洛(MCMC)理论的随机反演方法,通过先验信息建立解空间,采用马尔科夫链扰动模拟算法,计算模型频散曲线与实际曲线的似然函数,再利用转移概率规则引导搜索过程,可得到与实际频散曲线匹配最佳的最优解,获得大概率的横波速度后验概率分布,从而减少反演的多解性,提高分辨率。本文通过典型理论模型和实测面波数据验证基于贝叶斯理论的随机反演方法的可行性。实测数据应用中,利用GPR解释剖面进行深度约束作为先验模型信息,在此基础上开展面波随机反演,从而获得近地表横波速度结构。

1 贝叶斯随机反演原理

根据贝叶斯理论可以用d表示实际数据,用m表示横波(V_s)速度模型:

(1)

\begin{matrix} p (m │ d) = p (d │ m) p (m) / p (d), \propto p (d│m) p (m), \end{matrix}

其中,p(m│d)是模型在给定数据的情况下的后验概率,p(m)是引入数据之前有关模型的先验信息,p(d│m)是似然函数(在给定特定模型(m)的情况下观察到实测数据的概率),p(d)是实际数据的准确性,与m无关可以视为常数。

实际数据由N对频率(f)和瑞利波相速度(PV)组成的频散曲线以及偏差(σ)组成:

(2)

d = [f_{1} f_{2} \dots f_{N}, P V_{1} P V_{2} \dots P V_{N}],

(3)

σ = [σ_{1} σ_{1} \dots σ_{N}] 。

相速度偏差(σ)是表示所拾取的频散曲线的不确定性的量,其大小可取频散能量图的各频率的相速度能量集中的区间长度。

使用Voronoi核来表示模型参数^[20]。如图1所示。灰色为v_s模型空间,每一层都包含一个约束核和多个浮动核,在每个深度的v_s值由离其最近的核决定,Voronoi核在模型空间中移动,其中浮动核可以在模型空间的深度范围内移动,而约束核只能在其固定的层移动,于是模型参数为:

(4)

\begin{array}{l} m = [k, d p_{1}, d p_{2}, \dots, d p_{k}, v_{s_{1}}, v_{s_{2}}, \dots, v_{s_{k}}, I, \\ d p_{c_{1}}, d p_{c_{2}}, \dots, d p_{c_{I}}, v_{s_{c_{2}}}, \dots, v_{s_{c_{I}}}], \end{array}

图1

新窗口打开| 下载原图ZIP| 生成PPT

图1 使用Voronoi核进行模型参数化

a—无先验信息的v_s均匀半空间模型;b—有先验信息的三层v_s空间模型

Fig.1 Model parameterization using Voronoi nuclei

a—1 layer model without constraints;b—3 layer model with GPR constraints

式中,k是浮动核的数量,I为约束核的数量,dp_k为浮动核的深度, $v_{s_{k}}$ 是浮动核的横波速度,d $p_{c_{i}}$ 是约束核的深度, $v_{s_{c_{i}}}$ 是约束核的横波速度。

先验信息p(m)可以分为两部分:

(5)

p (m) = p (m | n) p (n) 。

其中,p(n)是Voronoi核数的先验,设区间I={n∈N|n_min<n≤n_max},因此:

(6)

p (n) = \{\begin{array}{l} 1 / (Δ n), n \in I; \\ 0, n \notin I 。 \end{array}

其中Δn=(n_max-n_min),由于Voronoi核的深度c和速度v相互独立,所以p(m|n)可拆分为:

(7)

p (m | n) = p (c | n) p (v | n) 。

速度的先验区间为J={v_i∈R|v_min<v_i≤v_max},因此:

(8)

p (v_{i} | n) = \{\begin{array}{l} 1 / (Δ v), v_{i} \in J; \\ 0, n \notin J \end{array} 。

其中,Δv=(v_max-v_min)。由于每个Voronoi核的速度相互独立,则:

(9)

p (v | n) = \overset{n}{\prod_{i = 1}} p (v_{i} | n) 。

设深度的网格数为N,则n个Voronoi核在N个网格上有 $[\frac{N!}{n! (N - n)!}]$ 种可能性,每种可能性的概率相等,因此:

(10)

p (c | n) = {[\frac{N!}{n! (N - n)!}]}^{- 1} 。

结合式(5)~(10)先验概率密度函数:

(11)

p (m) = \{\begin{array}{l} \frac{n! (N - n)!}{N! {(Δ v)}^{n} Δ n}, n \in I 且 \forall i \in [1, n], v_{i} \in J; \\ 0, n \notin I 或 \exists \forall i \in [1, n], v_{i} \notin J 。 \end{array}

似然函数p(d│m)表示在给定模型m的情况下观测到实际数据的可能性,实际上是通过对模型m进行正演模拟得到频散曲线G(m)并计算其与实际频散曲线的相似度来实现的。假设数据d_i (i=1,2,…,N_data)的每对频率f_i与相速度PV_i相互独立,则似然函数:

(12)

p (d │ m) = \overset{N_{data}}{\prod_{i = 1}} \exp \{- \frac{[P V_{i} - G_{f_{i}} {(m)]}^{2}}{2 σ (f_{i})^{2}}\},

使用马尔科夫链蒙特卡洛法(MCMC)对后验概率分布进行采样^[21],该方法会产生一个长的模型链,每次迭代都会产生一个新的模型m',该模型与当前模型的差异包括以下4个扰动方式(图2):①随机改变某个核的v_s值,v_s扰动概率密度q_v₁( $v_{i}^{'}$ |v_i)= $\frac{1}{δ_{v_{s}} \sqrt[]{2 π}}$ exp $[- \frac{(v_{i}^{'} - v_{i})^{2}}{2 δ_{v_{s}}^{2}}]$ ; ②随机改变某个核的深度,深度扰动概率密度q_c( $c_{i}^{'}$ |c_i)= $\frac{1}{δ_{dp} \sqrt[]{2 π}}$ exp $[- \frac{(c_{i}^{'} - c_{i})^{2}}{2 δ_{dp}^{2}}]$ ,浮动核可以在整个深度范围内移动,但约束核只能在其固定层移动;③产生一个新的浮动核,在 $c_{n + 1}^{'}$ 深度产生一个新核的概率q( $c_{i}^{'}$ │m)=1/(N-n),新核的速度为 $v_{n + 1}^{'}$ 的概率q(v'│m)= $\frac{1}{δ_{birth} \sqrt[]{2 π}}$ exp $[- \frac{(v_{n + 1}^{'} - v_{i})^{2}}{2 δ_{birth}^{2}}]$ ; ④随机移除一个浮动核,移除深度为c_i的浮动核的概率q( $c_{i}^{'}$ │m)=1/k,移除的深度为c_i的浮动核的速度为v_i的概率q(v'│m)=1。式中: $δ_{v_{s}}$ 是v_s的最大扰动值,δ_dp是深度的最大扰动值,v_i是基于当前Voronoi核分布的v_s值,c_i为当前Voronoi核深度,N为深度的网格数,n为当前Voronoi核数量,k为当前浮动核数量,σ_birth为新浮动核的v_s最大扰动值。

图2

新窗口打开| 下载原图ZIP| 生成PPT

图2 模型的4个扰动方式

a—改变核的v_s;b—改变核的深度;c—产出一个浮动核;d—移除一个浮动核

Fig.2 Illustration of four possible perturbations to a current model

a—change v_s of a nuclcus;b—move a nuclcus to a different depth;c—give birth to a new floating nucleus;d—remove a floating nucleus

对每个新模型进行正演计算以确定其是否达到可接受的标准:

(13)

a = \min [1, \frac{(p (m') q (m | m') p (d | m')}{p (m) q (m'| m) p (d | m)} \cdot | J |],

式中:q(m'|m)是从模型m变化到模型m'的概率,q(m|m')是从模型m'变化到模型m的概率,J为雅可比矩阵。对于不改变模型尺度,改变核的v_s(图2a)和改变核的深度(图2b),其扰动是对称的,即从m变化到m'的概率等于从m'变化到m的概率:

(14)

q_{v 1} (v_{i}^{'} │ v_{i}) = q_{v 1} (v_{i} │ v_{i}^{'}),

(15)

q_{c} (c_{i}^{'} │ c_{i}) = q_{c} (c_{i} │ c_{i}^{'}),

所以这两种扰动:

(16)

q (m'| m) = q (m | m') 。

对于产生一个新的浮动核(图2c)和移除一个浮动核(图2d),模型尺寸发生了变化,且扰动不对称。由于Voronoi核的深度和速度相互独立,对于产生新核:

(17)

\frac{q (m | m')}{q (m'| m)} = \frac{q (c | m')}{q (c'| m)} \frac{q (v | m')}{q (v'| m)},

(18)

{[\frac{q (m | m')}{q (m'| m)}]}_{birth} = \frac{\sqrt[]{2 π} (N - n) δ_{birth}}{(k + 1)} \exp [\frac{(v_{n + 1}^{'} - v_{i})^{2}}{2 δ_{birth}^{2}}] 。

从式(18)可以看出,产生一个新的浮动核时,随着新核的速度 $v_{n + 1}^{'}$ 远离原模型同深度的速度v_i, ${[\frac{q (m | m')}{q (m'| m)}]}_{birth}$ 呈指数增加。因此,产生新核的扰动促进速度变化,与式(13)中似然函数比项结合时,这种指数增长的概率密度最终会受到压制,以惩罚没有使数据拟合改善的较大速度扰动。

对于移除浮动核:

(19)

\begin{matrix} {[\frac{q (m | m')}{q (m'| m)}]}_{death} = \frac{k}{\sqrt[]{2 π} (N - n + 1) δ_{birth}} \exp [- \frac{(v_{j}^{'} - v_{i})^{2}}{2 δ_{birth}^{2}}], \end{matrix}

式中: $v_{j}^{'}$ 为移除一个浮动核后相同深度的速度值。随着 $v_{j}^{'}$ 远离原模型同深度的速度v_i, ${[\frac{q (m | m')}{q (m'| m)}]}_{death}$ 呈指数减小,因此,对于移除后模型速度变化小的浮动核将更有可能被移除。

式(13)中,J是从模型m变化到模型m'的雅可比矩阵,表示m与m'的尺度变化,即只有在两个不同尺度的模型之间扰动(产生或移除Voronoi核)时才需计算雅克比矩阵,如果当前模型与扰动模型具有相同尺度,则|J|=1,可以忽略。对于产生Voronoi核,从m到m'的双射映射h:

(20)

(c, v, u_{c}, u_{v}) \leftrightarrow (c, v, c_{n + 1}^{'}, v_{n + 1}^{'}) = m',

随机变量u_c∈[1,2,…,N-n]用于产生新核的位置 $c_{n + 1}^{'}$ ,随机变量u_v~N(0, $δ_{birth}^{2}$ )用于分配给新核的速度:

(21)

v_{n + 1}^{'} = v_{i} + u_{v},

v_i是新核位置扰动前的速度值。模型空间分为离散空间(核深度)和连续空间(速度),u_c是用于离散空间转换的离散变量,而u_v是用于连续空间转换的连续变量。Denison等^[22]研究表明对于离散变换,雅克比项|J|=1,可以忽略,所以雅克比项仅考虑以下变量:

(22)

(v, u_{v}) \leftrightarrow (v, v_{n + 1}^{'}) = v',

因此:

(23)

\begin{matrix} | J |_{birth} = |\frac{δ (v')}{δ (v, u_{v})}| = |\frac{δ (v, v_{n + 1}^{'})}{δ (v, u_{v})}| = |\frac{δ (v_{i}, v_{n + 1}^{'})}{δ (v_{i}, u_{v})}| = |\begin{array}{l} 1 & 0 \\ 1 & 1 \end{array}| = 1 。 \end{matrix}

对于移除Voronoi核,|J|_death=|J^-1|_birth=1,所以对于上述4种扰动,雅克比项都等于1,可以忽略。

在求得新模型的a值后,将其与随机数u~U[0,1]进行比较:如果a>u,则接受新模型m'并将其添加到链中,如果a<u,则保留现有模型m并寻找新的模型。似然函数p(d|m')大的扰动通常会被接受,反之则极有可能会被拒绝,但偶尔也会被接受。

面波随机反演的基本流程如图3所示,其基本实现流程如下:

图3

新窗口打开| 下载原图ZIP| 生成PPT

图3 面波随机反演基本流程

Fig.3 Basic flowchart of surface wave stochastic inversion

1)根据地震面波资料提取频散曲线;

2)根据先验信息建立解空间;

3)对解空间内的随机模型进行预反演,其结果作为反演的初始模型m;

4)计算初始模型m的理论频散曲线与似然函数;

5)随机扰动产生新模型m',并计算理论频散曲线及似然函数;

6)计算a,若a>u(0,1),则将新模型m'加入到马尔科夫链中,否则保留原模型m;

7)重复步骤5)、6),直至满足迭代停止条件。

2 面波频散随机反演模型测试

为验证本文所述随机反演方法,分别采用4层v_s递增模型和4层含低速夹层模型开展反演测试。表1给出了4层v_s递增模型的基本参数。图4为4层v_s递增模型基阶频散曲线的反演结果。图4a为无约束反演结果,反演结果(最佳模型)与真实模型在深度7 m的第2、3层分界面处有所差异,v_s的后验概率分布与真实模型基本拟合;图4b为带约束反演结果,v_s的先验约束根据频散曲线设为300~1 300m/s(v_s的下边界为0.9倍的相速度最小值;v_s的上边界为1.1倍的相速度最大值)。深度约束为3层,分别为0~2 m、3~15 m、16~30 m,相比于无约束反演结果,带约束反演最佳模型和v_s的后验概率密度分布与真实模型更好地吻合。图4c为无约束反演的反演结果(最佳模型)与实际模型频散曲线,两条曲线在所有频率范围内基本重合。图4d为有约束反演的最佳模型与实际模型频散曲线。图4e为两种模式下的迭代误差曲线,有约束反演(红线)迭代1 000次左右就能将误差降低到无约束反演(蓝线)迭代8 000次的误差大小,且目标函数最终的迭代误差也更小。该模型测试验证了对于速度递增模型,本文提出的随机反演策略即使在无约束条件下也能得到较好的反演结果。

表1 4层v_s递增模型参数

Table 1 4-layer v_s incremental model parameters

层序号	层厚度/m	v_s/(m·s^-1)	v_p/v_s	密度/(g·cm^-3)	v_s约束/(m·s^-1)
1	2	400	2	1.8	300~1300
2	5	600	2	1.8	300~1300
3	8	800	2	1.8	300~1300
4	15	1200	2	1.8	300~1300

新窗口打开| 下载CSV

图4

新窗口打开| 下载原图ZIP| 生成PPT

图4 4层v_s递增模型基阶频散曲线反演

a—无约束反演结果;b—有先验约束反演结果;c—无约束反演的最佳模型与实际模型频散曲线;d—有先验约束反演的最佳模型与实际模型频散曲线;e—迭代误差曲线

Fig.4 Fundamental dispersion curve inversion of four-layer v_s increasing model

a—inversion results without constraints;b—inversion results with prior constraints;c—dispersion curves of the best fitting model and true model for unconstrained inversion;d—dispersion curves of the best fitting model and the true model with prior constraint inversion;e—misfit iteration curve

表2给出了4层含低速夹层测试模型的基本参数。图5为4层含低速夹层模型基阶频散曲线的反演结果。图5a为无约束反演结果,最佳模型(反演结果)与真实模型误差较大;图5b为有约束反演结果,此处v_s的先验约束范围同样设为300~1 300 m/s,深度约束为3层,分别为0~2 m、3~17 m、18~30 m,最佳模型和v_s的后验概率密度分布与真实模型基本吻合。图5c为无约束反演的最佳模型与实际模型频散曲线,两曲线基本重合,最佳模型与真实模型的反演误差主要来自于基阶波对低速层的不敏感以及无约束随机反演局部搜索能力差。图5d为无约束反演频散曲线误差,可以看出在不同频率范围内存在一定的频散曲线拟合误差。图5e和图5f为带约束反演的最佳模型与实际模型频散曲线和频散曲线误差,相比无约束反演,有约束反演误差较小且整体几乎等于零,表明了约束随机反演具有较高的反演精度。图5g为两种模式下的迭代误差曲线,无约束反演(蓝线)在迭代8 000次左右时误差达到最小且趋于稳定,而有约束反演(红线)迭代3 000次左右就能达到同样的误差大小且最终的迭代误差也更小。上述模型测试表明,对于含低速夹层的速度结构,由于面波频散信息复杂,存在模态混叠现象,需要在一定约束条件下才能得到满意的反演结果。

表2 4层含低速夹层模型参数

Table 2 4-layer model with low velocity interlayer

层序号	层厚度/m	v_s/(m·s^-1)	v_p/v_s	密度/(g·cm^-3)	v_s约束/(m·s^-1)
1	2	600	2	1.8	300~1300
2	15	800	2	1.8	300~1300
3	8	400	2	1.8	300~1300
4	5	1200	2	1.8	300~1300

新窗口打开| 下载CSV

图5

新窗口打开| 下载原图ZIP| 生成PPT

图5 4层含低速夹层模型基阶频散曲线反演

a—无约束反演结果;b—有先验约束反演结果;c—无约束反演的最佳模型与实际模型频散曲线;d—无约束反演频散曲线误差;e—有先验约束反演的最佳模型与实际模型频散曲线;f—带约束反演频散曲线误差;g—迭代误差曲线

Fig.5 Fundamental dispersion curve inversion of low-velocity layer model

a—stochastic inversion results without constraints;b—stochastic inversion results with prior constraints;c—dispersion curve of the best fitting model and true model for unconstrained inversion;d—error of dispersion curve for unconstrained inversion;e—dispersion curves of the best fitting model and the true model with prior constraint inversion;f—error of dispersion curve with prior constraint inversion;g—misfit iteration curve

3 实测地震数据测试

为进一步验证随机反演方法的可行性与可靠性,采用浅层主动源面波实测数据开展横波速度反演测试。实测数据来自于中国安徽省淮南市某校园操场^[23],开展的地球物理方法包括浅层主动源面波与探地雷达(GPR)方法。面波数据沿测线共24炮,采用24个4 Hz检波器,道间距为1 m,最小偏移距10 m,震源采用锤击震源;GPR数据包含以第13炮集记录测量点为中心的共中心点道集(CMP)记录以及覆盖整个地震测线的共偏移记录。图6a为处理探地雷达数据得到的GPR剖面。图6b为基于卷积神经网络(CNN)^[24]进行GPR层界面自动识别结果,根据GPR剖面与CNN识别结果将地下划分为4层,黑色实线为层界面。图7a为经过处理仅保留面波成分的第13炮地震记录。图7b为由炮集记录提取的频散能量,根据能量的最大值提取出基阶频散曲线。图7c无约束反演结果,根据v_s的后验概率分布可以看出,在10 m深度范围内地下共4层,其中第3层为低速夹层,深度为4 ~8 m,反演的最佳模型在浅层与v_s后验概率分布基本匹配,但深层误差较大。图7d为GPR结果深度约束反演结果,v_s约束根据频散曲线设为100~500 m/s,与无约束反演相比,第3层的v_s后验概率分布更加集中,且与反演的最佳模型更加匹配。图 7e为迭代误差曲线,可见在有约束条件下,随着迭代次数的增加,GPR约束反演的收敛速度更快。

图6

新窗口打开| 下载原图ZIP| 生成PPT

图6 探地雷达处理剖面(a)与基于卷积神经网络的GPR层界面识别结果(b)

Fig.6 Processed GPR profile(a)and GPR layer recognition results based on Convolutional Neural Networks(b)

图7

新窗口打开| 下载原图ZIP| 生成PPT

图7 实测数据反演

a—第13炮地震记录;b—由炮集记录提取的频散能量;c—无约束反演结果;d—GPR深度约束反演结果;e—迭代误差曲线

Fig.7 Inversion results of real surface wave data

a—one shot gather of 13^th shot;b—extracted dispersion curve;c—unconstrained stochastic inversion result;d—stochastic inversion with GPR constrained;e—iterations

将上述一维反演结果按照空间位置分布组合成拟二维剖面,如图8所示。图8a为无约束反演结果组合成的拟二维剖面,黑色实线为GPR分层界面,从图中可以看出:面波频散反演的横波速度结果前两层分界面与GPR数据基本匹配,但后3层分界面误差较大。图8b为GPR深度约束反演结果组合成的拟二维剖面,前3层分界面与GPR数据基本匹配,虽然底层界面仍然存在误差,但相比无约束反演,误差已经有所降低。

图8

新窗口打开| 下载原图ZIP| 生成PPT

图8 拟二维剖面反演结果对比

a—无约束反演结果组合成的拟2D剖面;b—GPR深度约束反演结果组合成的拟2D剖面

Fig.8 Comparison of quasi-two-dimensional profile inversion results

a—2D velocity profile of unconstrained stochastic inversion;b—2D velocity profile of GPR constrained stochastic inversion

结合前人钻孔研究成果^[25],地下结构分层解释结论为:第1层由松软的沙子、砾石、土壤组成,深度为1~2.2 m,横波速度为100~150 m/s;第2层由紧凑且致密的沙子、砾石、土壤组成,深度为2.2~4.2 m,横波速度为350~500 m/s;第3层为黏土,深度为4.2~8 m,横波速度为120~180 m/s;基岩层由风化程度不同的砂岩组成,深度大于8 m,横波速度为400~500 m/s。前人钻孔深度解释与本文约束随机反演横波速度结构有较好的吻合性。

4 结论

本文提出了基于贝叶斯反演理论的随机反演方法并应用于面波频散曲线反演。利用先验信息深度约束分别对理论和实际面波数据进行无约束、有约束随机频散反演对比测试,结果表明:

1)随机反演方法是一种高效、稳健的反演方法,具有较高的全局寻优能力,且不依赖于初始模型。

2)带约束的随机反演可以有效且可靠地应用于面波频散曲线反演中,与传统无约束反演相比,随机反演可以很好地融合先验信息,减少多解性,提高反演效率与反演精度。

本文所使用的先验约束为GPR深度约束,而v_s速度范围约束是由频散曲线确定的,如果实际测区有地质资料可以确定地下各层的v_s速度范围,将其加入到先验约束中可以进一步减少多解性。

参考文献

原文顺序

文献年度倒序

文中引用次数倒序

被引期刊影响因子

[1]

李庆春, 邵广周, 刘金兰, 等.

瑞雷面波勘探的过去、现在和未来

[J]. 地球科学与环境学报, 2006,28(3):74-77.