1 研究背景

地震波在非完全弹性介质中传播时，除了受到与地震波传播特性有关的衰减外，还会受到地层的吸收作用^[1]，导致地震波能量发生变化，影响地震资料的信噪比和分辨率。地震勘探中，通常使用品质因子Q表征介质的吸收衰减性质^[2-3]，并将其倒数定义为地震波传播一个波长(或周期)距离时，因介质的非弹性所消耗的能量与在该周期内介质所积累的最大弹性应变能之比，用公式可表示为

$ \frac{1}{Q}=\frac{1}{2 \pi} \frac{\Delta E}{E} $

(1)

式中：E是处在最大应力和应变状态下的弹性能；ΔE是谐振激励下，每一个周期内的振动能耗。Q值越小，地层对地震波的吸收衰减量越大，Q值越大，地层对地震波的吸收衰减量越小。

由于压实作用弱、结构松散，近地表地层的Q值一般很小^[4]，因此尽管厚度相对于深部地层要小得多，但地震波穿过表层介质时遭受的吸收衰减却十分严重。研究表明，表层介质对地震波高频成分的吸收约占地层吸收总量的80%，是降低地震资料分辨率的重要因素之一^[5]。另外，复杂地表区表层的岩性、速度、厚度在横向上往往变化很大，不同位置接收的地震信号所遭受的吸收衰减程度不同，会导致地震记录的振幅、波形存在横向不一致。这种表层吸收带来的不一致可能远大于储层因素及油气因素引起的空间变化，解决不好会影响叠加成像，甚至引起地质空间的假象，进而导致错误的解释。

为了消除近地表吸收衰减对地震资料的影响，研究人员相继开展了大量近地表吸收补偿方法的研究，并在应用中见到了良好的效果^[6-9]。这些近地表吸收补偿技术无一例外，都以准确调查近地表Q值并建立高精度近地表吸收模型为前提。目前，近地表Q调查大多采用微测井法采集数据，利用频谱比率法、质心频率频移法、峰值频率频移法等估算Q值^[10]。Q值估算精度会受到震源子波与检波器耦合响应不一致、近场效应、(面波、浅层折射波和虚反射等)多种噪声以及旅行时间拾取误差或速度估算误差等因素影响，很难准确求取^[11-14]。

针对提高近地表Q值估算精度问题，李伟娜等^[14]借鉴微测井分层速度回归分析思想，提出基于微测井非相邻道谱比对数双线性回归Q估计方法，可降低时间拾取误差或速度估计误差等影响，使Q值估计结果更稳定。Wang等^[15]基于谱比法原理，通过积分运算，提出对数谱面积差(LSDA)Q值估算方法，与谱比法和质心频移法相比，在Q值估算的抗噪性和稳定性方面均具有明显优势。李君君等^[16]基于质心频移法的变式(线性衰减式)，提出一种通过拟合质心频率—传播时间的斜率估算Q值的方法，一定程度上克服了激发震源的非一致性，具有一定的理论可行性和鲁棒性。罗勇等^[17]设计了一种新的单井微测井观测系统，通过模型试算确认质心频率及其方差与炮检距满足三次多项式关系，提出一种基于多项式拟合的质心频移法估算近地表等效Q值，得到与模型等效Q值接近的结果。Li等^[18]提出一种井地联合近地表Q因子层析反演方法，利用所提出的几何结构降低接收端耦合对Q估计的影响，然后通过在炮道集中选择参考道消除源特征对Q反演的影响。翟桐立等^[19]针对常规微测井近地表Q值估计时面临的激发与接收耦合的非一致性、虚反射及近场问题等，提出一种深井激发、浅井和地面短排列联合接收的井地一体化测量的近地表Q值估算方法，避免了低速层底部界面虚反射对子波的干涉，最大程度地保证激发与接收的耦合一致性。于承业等^[20]提出一种利用双井微测井资料估算近地表Q值的方法，通过建立、求解地面和井下检波器的峰值频率变化方程组，获得近地表Q值的解析解。

目前应用的近地表Q值调查方法，一般采用雷管或炸药激发，但在一些民爆物品使用敏感地区，这种激发方式受到限制。基于重锤激发、多级井中检波器接收的地面激发—井中接收微测井，能够规避民爆物品使用的安全风险，已在近地表速度调查中得到应用。但由于存在激发子波不一致和激发噪声干扰强等问题，直接应用于近地表Q值调查，其估算精度难以保证。针对上述问题，本文在常见地面激发—井中接收观测系统的基础上，尝试在井口位置增设延迟校正检波器、在相邻炮次接收范围的交接处设置一个重复接收点的改进观测系统，并提出一套基于该改进观测系统的近地表Q值估算方法。该方法利用地面校正检波器消除记录延迟时差，利用重复接收道进行整形滤波消除不同炮次激发子波非一致性，通过多项式拟合减小初至时间拾取误差和干扰波的影响，从而提高Q值估算结果的稳定性，获得更可靠的近地表Q值。

2 方法原理 2.1 近地表Q值调查野外施工方法

微测井观测方法能够准确和有针对性地对近地表不同层位进行观测，是一种精细、可靠的近地表波场观测方法，广泛应用于近地表速度结构调查和近地表吸收衰减调查。地面激发—井中接收是常见的微测井观测方法之一。根据施工方式，可分为常规井中微测井和基于多级井中检波器的井中微测井。常规地面激发—井中接收微测井如图 1a所示，激发点与接收井之间的距离为x，接收井中布设n个接收点。施工时，在全井段布设检波器后，地面进行一次激发，完成数据采集。基于多级井中检波器的井中微测井，采用地面重复位置多次激发、井中分段接收的方式进行数据采集。图 1b为采用多级井中检波器接收的井中微测井观测系统示意图，井中n个接收点自下而上从1开始顺序编号。第N炮激发时，井中检波器置于编号为(N－1)×i+1~N×i检波点对应的井段接收。第N炮采集完成后，将检波器提升i倍道距的距离，再进行下一炮采集，直至完成全井段数据采集。

采用多级井中检波器施工不需要封井，检波器可以重复利用，而且完成表层调查后的井孔可用作常规地震的激发井，具有明显的经济优势。但采用多级井中检波器接收时，不同激发炮次之间会存在一定的随机记录延迟时差和激发子波非一致性问题。为了能在后续处理中消除不同炮次之间的记录延迟时差和激发子波非一致性问题，提高该方法Q值估算的精度，本文对采用多级井中检波器的井中微测井观测系统进行了如下改进。

(1) 在地面增加一个激发延迟时差校正检波器(图 1c)，每次激发时井中检波器与校正检波器同时接收。地面延迟时校正检波器可以位于井口，也可置于附近的其他位置，但需要在所有炮次激发完成前保持位置固定不变。

(2) 在相邻炮次接收范围的交接处，设置一个重复接收点，用于消除相邻炮次激发子波非一致性。第N炮激发时，井中检波器置于编号为(i－1)(N－1)+1~(i－1)(N－1)+i的检波点位置接收(图 1c)。第N炮采集完成后，将检波器提升i－1倍道距后进行下一炮采集，直至遍历全井段。

2.2 锤击激发记录延迟时差校正

利用延迟时差校正检波器接收到的地震波初至记录，计算出标准的系统延迟时间，进而完成对各炮次初至时间的校正，消除不同炮次之间的记录延迟时差，具体步骤如下。

(1) 求取相对延迟时

以第1炮延迟校正检波器接收的初至记录为标准道，对其他炮延迟验证检波器接收的初至记录按照下式进行互相关，得出相关后各炮的峰值时间，即为各炮的相对延迟时

$ R(n)=\frac{1}{N} \sum\limits_{m=0}^{N-1} X_1(m) X_i(m+n) $

(2)

式中：N为采样点数；n和m为时间采样点序号，其中n在－(N－1)~N－1之间变化；X₁、X_i分别为第1炮和任意第i炮激发相应延迟校正检波器接收的初至记录；R(n)为相关后的地震记录。

(2) 标准系统延迟时求取

根据炮次数量，分别采用不同计算方法剔除延迟时的相对奇异点，在剔除奇异点后计算标准系统延迟时。首先，分别利用下面两式

$ \mu=\sum\limits_{i=1}^N \tau_i $

(3)

$ \sigma=\sqrt{\frac{\sum\limits_{i=1}^N\left(\tau_i-\mu\right)^2}{N}} $

(4)

计算均值和标准差。式中：τ_i为第i炮的相对延迟时；μ为各炮相对延迟时的均值；σ为各炮相对延迟时标准差。当炮数大于15时，利用拉依达(PauTa)准则迭代剔除各炮次相对延迟时的奇异点，当炮数未超过15时，利用格拉布斯(Grubbs)准则法迭代剔除各炮次相对延迟时的奇异点。然后，计算剔除奇异点后的延迟时均值，即为标准系统延迟时。

(3) 延迟时差求取及校正

利用式(5)求取各炮的记录延迟时差，再根据式(6)对各炮不同接收道初至时间进行校正，消除不同炮次之间的记录延迟时差

$ \tau_i^{\prime}=\tau_i-\mu $

(5)

$ T_{0 i j}^{\prime}=T_{0 i j}-\tau_i^{\prime} $

(6)

式中：τ′_i为第i炮的系统延迟时；τ_i为第i炮的相对延迟时；μ为标准系统延迟时；T′_0ij、T_0ij分别为校正前、后的第i炮第j道的初至时间。

2.3 锤击激发子波一致性校正

通过整形滤波可以消除不同源记录的子波一致性差异^[21]。基本原理为：输入待校正数据和指定的期望输出数据，在一定的时间和空间范围内，利用维纳滤波算法求出整形滤波算子，然后将该算子应用于待校正数据，以消除该数据与指定的期望输出数据之间的差异。

本文采用的微测井调查方法，不同炮次之间存在子波非一致性问题。为了利于消除子波非一致性，在观测系统设计时，为相邻炮次设计一个重复接收点，即相邻两炮中前炮最后一个接收道与后炮第一个接收道位于相同的接收位置。根据整形滤波处理方法原理，将后炮第一道数据作为待校正数据、前炮最后一道数据作为期望输出数据，求出整形滤波算子，然后将该算子应用于后炮其他道数据，从而消除前、后两炮数据之间的激发子波非一致性。

2.4 Q值估算方法

根据地震波吸收衰减过程中高频成分的能量衰减速度远大于低频成分的特点，Quan等^[22-23]提出了质心频移法，定义α₀为吸收因子，即α₀=π/Qv，其中v为地震波在介质中传播的速度。

吸收因子与吸收系数的关系为α=α₀f。

定义震源信号振幅谱S(f)的质心频率f_S和方差σ_S²分别为

$ f_S=\frac{\int_0^{\infty} f S(f) \mathrm{d} f}{\int_0^{\infty} S(f) \mathrm{d} f} $

(7)

$ \sigma_S^2=\frac{\int_0^{\infty}\left(f-f_s\right) S(f) \mathrm{d} f}{\int_0^{\infty} S(f) \mathrm{d} f} $

(8)

则对应的接收信号振幅谱R(f)的质心频率f_R及其方差σ_R²分别为

$ f_R=\frac{\int_0^{\infty} f R(f) \mathrm{d} f}{\int_0^{\infty} R(f) \mathrm{d} f} $

(9)

$ \sigma_R^2=\frac{\int_0^{\infty}\left(f-f_R\right) R(f) \mathrm{d} f}{\int_0^{\infty} R(f) \mathrm{d} f} $

(10)

假设入射波信号振幅谱为高斯频谱，经介质衰减后f_R与f_S及σ_S²之间有以下关系

$ f_R=f_S-\sigma_S^2 \int_{\text {ray }} \alpha_0 \mathrm{~d} l $

(11)

经推导，可得

$ Q=\frac{\pi t \sigma_S^2}{f_S-f_R} $

(12)

式中t为地震波从激发点到接收点的传播时间。

如果不考虑初至拾取误差及干扰波等其他因素对初至波频谱的影响，采用本文方法获得的资料，经过记录延迟时校正和激发子波一致性校正后，将近炮点道初至波视作震源信号，远炮点道初至波视作接收信号，将两道初至时差视作地震波在两个地震道接收点对应地层中的传播时间，可以利用式(12)计算出任意两道之间地层的Q值。

然而，在实际应用中，地面激发—井中接收的近炮点道初至波波形会受到面波干扰的严重影响，初至时间也会因传播路径复杂化出现不正常的时—深关系，无法用于Q值估算。图 2为采用地面激发—井中接收观测方式得到的炮集记录，激发点与接收井相距7.5 m。图中6 m以上深度接收的信息受到了面波干扰，其初至波形与6 m以下深度接收的信号有明显差异。用6 m和18 m深度接收的初至波估算Q值为3.026，用2 m和18 m深度接收的初至波估算Q值为－0.399，该结果显然不正确。

针对上述问题，基于黄土塬、沙漠区和戈壁砾石区近地表为连续介质的表层结构特点，提出采用多项式拟合求取浅表层初至时间、初至波频谱质心频率和质心频率方差，进而估算浅表层Q值的方法。具体如下。

(1) 拾取初至。去掉存在干扰影响地震道的初至，进行多项式拟合求得时—深关系式。用拟合公式，代入深度，计算出各道初至时间，用于Q值估算。

(2) 截取初至波。计算各道初至波质心频率和质心频率方差，去掉存在干扰影响问题地震道的相关信息，进行多项式拟合得出初至波质心频率与深度、初至波质心频率方差与深度的关系公式。用拟合公式，代入深度，计算得出各道的初至波质心频率和质心频率方差，用于Q值估算。

3 基于正演的Q值估算方法验证 3.1 正演模拟

根据沙漠和黄土山区近地表结构特点，设计了如图 3a所示的低降速层具有近连续介质结构特征的正演模型，其浅层60 m范围内划分为12个5 m厚的小层，层速度从550 m/s等差递增到1100 m/s，60 m以下深度地层的速度设定为1600 m/s，地层Q按照经验公式^[4]给定，即Q=14v^2.2，其中v为地震波纵波速度(km/s)。正演采用黏滞弹性波动方程，激发子波参考实际资料选择为100 Hz主频的雷克子波，观测系统为地面激发—井中接收方式。接收井位于模型中部，井深为80 m，井中自下而上按2.5 m点距布设33个接收点，激发点布设在距离接收井5 m的位置。

图 3b为正演得到的模拟地震记录，从图中可看出，地面附近接收时初至波受到面波干扰，45~55 m深度接收时初至波受到来自1600 m/s速度层顶界面(深度60 m)反射波的干扰。

3.2 Q值估算方法分析及验证

利用傅里叶变换求得图 3b中各道初至波振幅谱，再利用式(9)和式(10)求取各道质心频率和方差。图 4展示了各道初至波质心频率和方差与接收点深度的对应关系以及剔除异常值前后进行多项式拟合得到的曲线。其中，剔除异常值前后的拟合曲线分别对应式(13)和式(14)。

从图 4可见，初至波质心频率和方差与接收深度存在非线性关系，与多项式整体有较高的拟合度，受面波和反射波干扰影响的地震道初至波质心频率和方差随接收点深度变化的趋势与其他道有明显差异，剔除异常值后拟合的曲线形态能够更客观的反映初至波质心频率和方差随接收点深度变化的趋势。

$ \left\{\begin{aligned} f_S= & 0.00000002 h^5-0.00000187 h^4-0.00014389 h^3+0.03067443 h^2- \\ & 1.70346890 h+94.03400180 \quad R^2=0.99480993 \\ \sigma^2= & 0.00000244 h^5-0.00044119 h^4+0.02232670 h^3+0.11263277 h^2- \\ & 35.62080949 h+1419.62915352 \quad R^2=0.99446251 \end{aligned}\right. $

(13)

$ \left\{\begin{array}{c} f_S=-0.00000003 h^5+0.00000912 h^4-0.00106632 h^3+0.06515671 h^2- \\ 2.24686674 h+96.64488506 \quad R^2=0.99922103 \\ \sigma^2=0.00000025 h^5+0.00010938 h^4-0.01743190 h^3+1.31412110 h^2- \\ 50.00586177 h+1469.07739113 \quad R^2=0.99935532 \end{array}\right. $

(14)

拾取图 3b中地震记录各道初至(初至波主瓣峰值时间)，并进行多项式拟合，校正初至拾取误差后计算出前后道传播时差。图 5展示了各道初至时间及时差与接收点深度的对应关系以及剔除异常值前后进行多项式拟合得到的曲线。其中，剔除异常值前后的拟合曲线分别对应式(15)和式(16)。

从图 5中不拟合初至时间及时差随深度变化的关系可看出，在接收点深度小于10 m和大于60 m时，初至时间及时差随深度的变化趋势与其他深度有一定的差异，在深度50~57.5 m时，初至时差存在异常抖动变化。分析认为，深度小于10 m时的异常是受非零偏移距观测带来的影响，深度大于60 m以后的异常是由于该范围地层不再具备连续介质特征带来的影响，50.0~57.5 m深度的异常是反射波干涉初至波带来的影响。

$ t=0.00000163 h^4-0.00031405 h^3+0.01410365 h^2+1.02309818 h+9.01303657, R^2=0.99980057 $

(15)

$ t=0.00000067 h^4-0.00008406 h^3-0.00217205 h^2+1.43895419 h+5.89562472, R^2=0.99999212 $

(16)

将存在异常的初至剔除后进行拟合得到的结果与不剔除异常值拟合的结果存在较大差异。前者更能反映出深度小于60 m近连续介质初至时间随深度变化特点，但在深度大于60 m范围拟合初至时间与拾取的初至时间存在较大偏差，不能正确地反映不足60 m的浅层初至时间随深度变化的趋势。

利用式(13)和式(15)求得各道质心频率方差和初至时间，命名为不剔除异常拟合方案。利用式(14)和式(16)求得各道质心频率方差和初至时间，命名为剔除异常拟合方案。以前道接收信号为震源信号，后道接收信号为衰减后的接收信号，分别将不拟合、不剔异常拟合及剔除异常拟合方案得到的质心频率、方差及初至时差代入式(12)计算得到3套不同深度地层的Q值(图 6、表 1)。图表中理论Q值为正演模型的地层Q值，用于对比不同方法计算Q值的精度。

从图 6中可以看出，在60 m深度范围内，剔除异常值拟合方案得到的Q值与理论Q值接近，不剔除异常值拟合方案得到的Q值与理论Q值在45~60 m深度段存在较大偏差，不拟合方案得到的Q值与理论Q值的偏差整体较大，且在45~60 m深度段偏差达到最大值。深度超过60 m，三套方案得到的Q值与理论Q值均存在较大偏差。

结合图 3、图 4和图 5分析认为：地表附近接收的信号受到干扰波影响、50~60m深度段接收的信号受反射波影响，导致它们的质心频率、方差和初至时间随深度变化的趋势异于整体趋势，是造成不拟合方案Q值与理论Q值的偏差大的主要原因；不剔除异常值拟合不能完全消除干扰影响，虽能减小Q值估算误差，但还存在不足；剔除异常值后拟合，能够消除近连续介质范围内(0~60 m深度)干扰波对质心频率、方差和初至时间的影响，达到提高Q值估算精度的目的；在60~80 m深度范围内，地层速度不再呈连续变化特征，且速度与上覆地层之间存在较大差异，质心频率、方差和初至时间随深度的变化关系与0~60 m深度范围内的变化关系不一致，用0~60 m的拟合公式求得的60~80 m范围内的质心频率、方差和初至与真实值存在较大偏差，导致拟合法不能正确估算60~80 m深度范围内的Q值。

基于上述认识，笔者认为，通过多项式拟合质心频率、方差和初至时间可以消除干扰影响，提高质心频率偏移法Q值估算精度，可适用于沙漠和黄土塬等近地表地层具有连续介质特征区域的近地表Q值调查，但当地层结构由连续介质转换为高速层状介质后该方法的适用性会下降。

4 实际应用

实际应用选择在塔里木盆地沙漠区，该区地表被起伏风成沙丘覆盖，近地表存在较稳定的潜水面。以潜水面为界可将近地表划分为上部疏松沙层和下部含水沙层。疏松沙层具有近连续介质特征，速度随深度增加逐渐增大，但总体较低；含水沙层在表层调查的探测范围内速度稳定在1750 m/s左右。

4.1 资料采集方案

观测系统如图 1c所示，接收井深50 m，接收点距1 m，激发点与接收井相距2 m。采用6级井中检波器接收，12 lb重锤地面激发10炮次。每次激发完成后井中检波器向上提升5 m，前炮最浅接收点与后炮最深接收点位置重复。

4.2 资料预处理

按照上述方法采集，在理想状态下延迟校正检波器接收到的初至波起跳时间相同、波形应一致。但在实际采集的资料上，不同炮次延迟校正检波器接收记录的起跳时间和波形存在一定的差异(图 7a)，反映了不同炮次之间存在着记录延迟时差和激发子波非一致性问题。

采用2.2节和2.3节所述方法对资料进行记录延迟时校正和激发子波一致性校正后，延迟校正检波器接收的记录起跳时间和波形一致性得到明显改善(图 7b)。

4.3 Q值估算

截取预处理后数据的初至波，进行傅里叶变换并计算得到各道初至波质心频率及方差。选择低降速层中接收道的信息，剔除异常值后拟合得到质心频率、方差和初至时间与深度的关系曲线(图 8)。利用拟合值，并对照近地表结构层位解释结果，计算得到各层层Q(表 2)。从表中可以看出，除地表下2 m深度内地层外，本文方法得到的Q值与经验公式^[4]得到的Q值接近。

4.4 应用效果

将求取的近地表等效Q值应用于实际地震资料，进行近地表吸收补偿处理。图 9为补偿前、后的叠前时间偏移剖面。从图中可以看出，经过表层吸收补偿的剖面分辨率得到了提高，能够更清楚地刻画薄层的地震响应。

5 结论及认识

(1) 基于锤击和井中多级检波器的井中接收微测井用于近地表Q值调查存在重复激发带来的记录延迟时差、激发子波非一致性以及干扰波场的影响等问题，采用该方法得到的资料直接用于Q值估算，得到的Q值误差较大；

(2) 通过采用增设了地面延迟校正检波器和在相邻炮次接收范围的交界处设置一个重复接收点的改进观测系统，以及基于该观测方法的激发延迟时差和激发子波一致性校正方法，能够消除重复激发带来的记录延迟时差、激发子波非一致性问题；

(3) 在具有连续介质特征的低降速层中，地震波质心频率和质心频率方差随深度的变化关系存在非线性关系，能够通过多项式拟合逼近，剔除异常值拟合的效果更好；

(4) 通过多项式拟合获得不同深度的质心频率、质心频率方差和初至时间，可解决井中接收微测井Q质估算面临的干扰影响问题，能够提高Q值估算精度，该方法可适用于沙漠和黄土塬等近地表地层具有连续介质特征区域的近地表Q值调查。