基于浅地层剖面的海底浅表层沉积物物理性质参数反演技术研究——以渤海海底管线路由区为例

黄必桂，李家钢，周庆杰，李西双，刘乐军，高珊，周航，张承艺

(1. 中海油研究总院有限责任公司，北京 100028；2. 自然资源部第一海洋研究所 自然资源部海洋地质与成矿作用重点实验室，山东 青岛266061)

1 引言

海底浅表层沉积物位于海底海水与沉积地层交界处，主要由砂、粉砂、黏土和孔隙流体（海水）等物质组成，处于生物的、化学的和物理的不断作用中，承载了大量且丰富的沉积环境信息[1]。对于海底浅表层沉积物物理性质的了解，目前常用的方法是采集足够多的沉积物样品，通过测试分析得到沉积物的物理参数或力学指标，或者使用原位测量的方式获得海底浅层土体的力学性质等参数[2]。上述方法可以直接准确地得到海底表层沉积物的物理性质，但无疑具有较高的成本，经济和时间上的限制导致只能进行有限的离散站位式采样。对于范围较大或者各向异性较强的区域来说，有限的采样可能无法提供足够的信息，进而造成对沉积物物理性质了解的偏差。因此，探索一种经济可行的、快速了解较大范围内海底浅层土物理性质的方法显得尤为重要。

浅地层剖面是利用声学方法了解沉积层的结构，声学信号包含了地下介质的弹性参数信息。反射强度、反射系数、波阻抗等声学属性与沉积物的物理性质密切相关[3]。在海床处（水和沉积物的分界面），浅地层剖面通常对应于第一个正极性的强反射，其强度取决于声阻抗和海床处的声反射系数[4-6]。而反射系数和声阻抗又与沉积物密度、含水量、孔隙比、粒度等物理性质密切相关[7]。因此，这些声学参数和沉积物性质之间的关系为利用浅地层剖面数据得到海底表层沉积物的物理性质架起了桥梁。

近年来，国外一些学者开始利用浅地层剖面数据进行反演沉积物物理属性的尝试性工作。Kim等[8]基于Chirp浅地层剖面和沉积物测试数据，建立了Ulleung盆地海底沉积物声学特征和岩土结构特征之间的关系，用于分析沉积物的分布格局和沉积过程。Schock[9]基于Biot模型利用Chirp浅地层剖面和实验室测量的沉积物物理性质提出了一种评估海底沉积物物理性质的方法。此外，国内也有学者应用 Biot-Stoll模型和原位测量数据探讨了声波在海底沉积物介质中的传播规律[10-13]，并将该模型用于反演沉积物粒径、孔隙度、密度等物性参数[14-15]。然而，Biot-Stoll模型在不同的海域采用不同的参数所预测的声学性质和物理性质之间的关系存在差异，因此，针对特定的海域利用实测物理参数构建Biot-Stoll模型，才能获得较为准确的声学和物理参数预测关系。

渤海是中国重要的陆架海区，受黄河等入海河流和渤海环流等多种因素的影响，其沉积物类型组成比较复杂，是我国近海沉积物特征研究的热点海域[16-18]。我们在渤海LD16-3CEPA至LD10-1PAPD路由段积累了丰富的浅地层剖面数据和表层沉积物取样测试数据（图1）。为此，本文以该区的海底表层沉积物的密度、孔隙度等取样测试数据为基础，建立Biot-Stoll模型计算海底反射系数与沉积物物理性质之间的相关关系。基于该区浅地层剖面计算得到的海底反射系数和Biot-Stoll模型建立的海底反射系数与沉积物物理性质之间的相关关系反演浅表层沉积物的孔隙度、密度、平均粒径等物理性质，并与取样点处实测物理性质结果进行对比，探讨该方法的可行性与准确度，为快速获取连续的海底沉积物物理性质提供新的方法参考。

图1 研究区位置及浅剖与取样站位分布Fig. 1 Location of study area and distribution of sub-bottom profiles and sampling stations

2 数据与方法

2.1 浅地层剖面数据

本文开展海底浅表层沉积物物理性质参数反演方法研究所使用浅地层剖面数据，是在渤海LD16-3CEPA至LD10-1PAPD路由段开展工程勘察过程中获取的。该浅剖数据是由Innomar SES2000-70参量阵浅地层剖面仪采集，主频率为2～8 kHz（图2），浅剖数据格式为SUN工作站的数据格式，4字节IBM浮点数记录每道采样点数为4500，采样间隔为0.014 ms。此外，为有效检验反演方法的准确性和适用性，选取浅剖测线区域内表层沉积物取样的工程地质测试数据进行对比，取样站位处的海底浅表层沉积物物理性质参数见表1。

表1 海底浅表层沉积物物理性质Table 1 Physical properties of seabed shallow surface sediments

图2 浅地层剖面频谱分析特征Fig. 2 Spectral analysis feature of sub-bottom profileb和c中黑色线为剖面整体频谱特征，蓝色线为其中1道的频谱特征The black lines are the overall spectral characteristics of the profile, and the blue lines are the spectral characteristics of one channel in b and c

2.2 Biot-Stoll模型法

Biot[19]于1956年提出了一种可以同时考虑介质的孔隙性与弹性的经典沉积物介质特性描述理论模型，在各向异性和具有黏弹性的双相饱和多孔介质流体计算及声学特性分析中应用广泛。1974年，Stoll[20]对Biot模型进行了改进，提出了Biot-Stoll模型，该模型可应用于海底沉积物介质的声速和声衰减计算，并推导出了简谐平面波在多孔介质中传播的方程为

式中，ρ为体密度；ρf为孔隙流体密度；ω=2πf为角频率；参数为宏观压力梯度下流体流动的相位（其中，c为曲折度，n为孔隙度）；j为虚数；Fη是用于解释泥沙孔隙中振荡流的频率依赖性黏性损失的一种黏性修正因子；k为复波数；H、C、M分别为Biot理论弹性模量、附加弹性模量和复弹性模量。

Schock[9]在Biot-Stoll模型基础上，假定孔隙流体在海底表层沉积物与海水的分界面处是沿与界面垂直的方向连续进出骨架的，并且界面总应力及界面流体压力平衡，在这些边界条件的限定下，可以计算得到海床的反射系数为

式 中， ρw和cw为 海水的密 度 和 声速；k1和k2为 快 慢纵波的波数；R为海床反射系数；Dr和Di分别为入射波和反射波的复位移振幅值；A1和A2分别为沉积物骨架在快纵波和慢纵波作用下的复位移振幅值；B1和B2分别为快纵波和慢纵波作用下孔隙流体相对于骨架运动的复相对位移；G1和G2分别为快纵波和慢纵波的复相位与位移比。

Biot-Stoll模型涉及的参数较多，参数的选取对计算结果有不同程度的影响。通过物理性质测试获取了颗粒密度、孔隙度、湿密度、粒径等参数，其他参数，如渗透率、沉积物孔隙因子、颗粒体积模量和剪切模量、孔隙水体积模量等，均需要利用经验公式计算或从文献资料中获取，详细参数取值如表2所示。

表2 Biot-Stoll模型输入的沉积物物理参数Table 2 The input sediment physical parameters of the Biot-Stoll model

2.3 反射系数求取方法

在浅地层剖面上，海底地震记录可以简单的认为是震源子波与海底反射系数相褶积的结果，可以利用地震反褶积技术进行子波处理消除子波的影响，得到实际反射系数序列[21]。由于对数分解法提取子波，避免了对子波进行假设，不需要知道子波的类型，也不需要假设反射系数是白噪，因此，本文中，我们利用对数分解法（同态求取法）计算浅地层剖面的海底反射系数，子波提取过程如图3所示。

图3 对数分解法求取子波的过程图解Fig. 3 Process diagram of wavelet extraction by logarithmic decomposition method

在求得地震子波b(t)后，进而计算其反子波b′(t)，然后将反子波b′(t)与地震记录x(t)褶积，即可求出反射系数r(t)为

海底反射系数计算结果如图4所示。

图4 典型浅地层剖面图（a）及该剖面计算得到的海底反射系数（b）Fig. 4 Typical sub-bottom profile (a) and the calculated seabed reflection coefficients (b)

3 结果讨论

3.1 基于Biot-Stoll模型建立的海底反射系数与物理性质关系

基于Biot-Stoll模型研究了海底反射系数随频率的变化关系，计算了频率为5 kHz时（Chirp浅剖主频为5 kHz）海底反射系数与沉积物孔隙度、密度、平均粒径之间的相关关系，建立拟合方程（图5）。图5a可以看出在低频（f＜103Hz）时，纵波波速受频率变化的影响较小，在中高频（103Hz＜f＜105Hz）时，纵波波速受频率的变化影响较大，几乎呈线性关系，而在频率大于一定值以后，纵波速度又会趋于稳定；图5b为海底反射系数与沉积物孔隙度的变化关系，随着孔隙度的减小，反射系数呈增大的趋势；图5c为海底反射系数与沉积物密度的变化关系，可以看出其近似呈线性关系，为提高拟合的准确度，本文以二次方程对其进行拟合；图5d为海底反射系数与沉积物平均粒径的变化关系，可以看出，反射系数与平均粒径呈负相关，随平均粒径的增大而减小。

图5 海底反射系数与沉积物物理性质的相关关系（主频为5 kHz）Fig. 5 Correlation between seafloor reflection coefficient and sediment physical properties, (the dominant frequency is 5 kHz)

3.2 物理性质反演结果与实测结果对比

对研究区所有的浅地层剖面进行计算，得到了该区海床浅表层沉积物物理性质如图6所示。对反演结果与取样站位处平均粒径、孔隙度和密度等测试结果的对比统计见表3，偏差对比见图7。

图6 海底表层沉积物物性反演结果与取样测试结果对比Fig. 6 Comparison between inversion results and sampling test results of seafloor surface sediment physical properties图中坐标为TM投影，中央经线为120°EThe coordinates in the figure are TM projection and the central meridian is 120°E

图7 反演结果与样品测试结果对比Fig. 7 Comparison between inversion results and sample test results

表3 反演结果与取样测试数据对比信息Table 3 Comparison information between inversion results and sampling test data

（1）平均粒径Ф反演结果总体较实测值偏大，最大偏差约为1.87（偏差率为-30.83%），总体偏差率在-30.83%～17.46%范围内。

（2）孔隙度反演结果偏差值最大约为0.087（偏差率为-13.26%），总体偏差在-13.26%～9.64%范围内。

（3）表层沉积物密度反演结果与取样测试结果的最大偏差值约为309.7 g/m3（偏差率为17.48%），总体偏差率在-2.64%～17.48%范围内。

4 结论

基于Biot-Stoll模型和渤海LD16-3CEPA至LD10-1PAPD路由段的浅地层剖面数据反演了该区海底浅表层沉积物的平均粒径、孔隙度、密度等工程地质参数，得到以下几点认识：

（1）结合研究区海底表层沉积物的工程地质特性，建立Biot-Stoll模型，得出主频为5 kHz时海底反射系数与沉积物孔隙度、密度、平均粒径之间的关系方程，且方程拟合度较高，可决系数R2均大于0.99，为利用海底反射系数反演沉积物物理性质提供了可靠的转换依据。

（2）对研究区的浅地层剖面数据利用子波反褶积的方法计算海底反射系数，该方法的关键是子波的提取。本文是利用对数分解法（同态求取法）提取浅地层剖面的子波，对数分解法提取子波的优点在于避免了对子波进行假设，不需要知道子波的类型，也不需要假设反射系数是白噪。

（3）对渤海LD16-3CEPA至LD10-1PAPD路由段的浅地层剖面数据计算其海底反射系数，并反演了海底表层沉积物孔隙度、密度、平均粒径等工程地质参数。反演结果与取样测试结果进行对比，可以较好地反映海底表层沉积物类型的变化。平均粒径反演结果的总体偏差率在-30.83%～17.46%范围内，孔隙度反演结果的总体偏差在-13.26%～9.64%范围内，表层沉积物密度反演结果的总体偏差率在-2.64%～17.48%范围内。

综合本文研究内容，基于浅地层剖面的海底表层沉积物物性反演结果可以有效地弥补离散取样点之间的空白，从平面分布上可直接看出海底表面沉积物物性参数的变化规律，可以较好地反映海底表层沉积物类型的变化。