可控震源具有激发容易控制、观测系统分布完整合理、可以进行密集观测等优点, 已经成为地震资料采集时采用的主要激发方式之一。但这类人工震源普遍存在激发能量小, 接收信号弱等问题, 加之强噪声干扰和大地的吸收衰减作用, 因而地震资料信噪比很低。为了提高此类地震资料的信噪比, BODINE^[1]首先将地震波束形成方法应用于震源端, 设计出用于形成和控制地震波束传播方向的系统, 并通过实验验证了系统的有效性。1977年, ARNOLD^[2]通过液压式可控震源延时激发单频控制信号验证了定向地震波的存在。2003年, 廖成旺等^[3]通过实验验证了信号叠加技术抑制噪声的有效性。2004年, 姜弢等^[4]首次提出了相控震源激发定向地震波理论, 通过模拟仿真和野外试验, 证明了相控震源激发的定向地震波方法的可行性和优越性^[5-7]。由于震源端地震波束形成成本高, 难以保证震源间的一致性, 人们发展了更便捷、高效的基于接收阵列的时域地震波束形成(Time-domain Seismic Beam-forming Based on Receiver Array, TSBBRA)方法。2012年, 姜弢等^[8]依据震源端波束形成原理, 首次提出TSBBRA方法, 在接收端形成了等同于震源端的地震波束。2014年, 葛利华^[9]通过多种模型仿真实验和实测信号的处理, 证明了TSBBRA方法的有效性, 并分析了不同参数对地震主波束方向和宽度的影响。后来在TSBBRA方法的基础上, 宋健等^[10]通过在接收端对各单炮记录进行余弦振幅加权, 解决了因叠加数据量大而导致的波束变窄的问题。贾海青等^[11]提出了一种全方向波束定向方法, 通过不同方向有效波场的叠加, 实现全方向有效信号的增强。贾海青等^[12]依据主波束方向内外信号与原始信号的局部相关性差异又提出了一种基于局部相关加权的时域地震波束形成方法, 消除了畸变信号对有效信号的影响。由以上研究结果可知, TSBBRA方法通过控制波束方向和宽度可提高主波束方向上有效信号的信噪比。

SERGIO等^[13]首先将奇异值分解(Singular Value Decomposition, SVD)方法应用于地震勘探领域, 提出了基于KL变换的SVD方法, 并基于特征图像滤波的概念分离VSP地震波场的上、下行波。1994年, 陈遵德等^[14]提出了应用倾角扫描叠加方法, 解决了SVD滤波法仅能提高水平同相轴信噪比问题, 还提出了用神经网络法确定SVD滤波参数的方法。1998年, 赵德斌等^[15]提出了f-x域SVD预测滤波法, 在频率域实现了随机噪声的衰减。2004年, 李文杰等^[16]将SVD方法应用到直达波和折射波衰减处理中, 解决了采用切除方法会丢失浅层有效地质构造信息的问题。2007年, 苑益军等^[17]利用自动追踪SVD方法消除了叠前资料中的强线性干扰。近年, 沈鸿雁等^[18-19]利用SVD方法实现了地震记录中有效信号和随机噪声的波场分离。以上研究表明, SVD方法对相干噪声和随机噪声的压制有显著的效果。

在TSBBRA方法中, 决定信号信噪比提高程度的关键参数为叠加记录个数N, 当目标信号信噪比较低时, 采用较小的N即可增强信号的能量, 提高记录的信噪比。但当目标信号的信噪比很低时, TSBBRA方法需要较大的N进行叠加, 但已有研究结果^[10]表明, N越大波束越窄, 波束过窄则不能覆盖全部检波器排列, 导致有效波束宽度外的信号畸变, 无法满足地震解释的要求。另外, 对于相干噪声和随机噪声强的地震记录, TSBBRA方法虽能压制相干噪声和随机噪声, 但噪声压制不均匀, 当地震解释对记录的整体信噪比要求较高时, 应用效果不佳。因此, 为了改善相干噪声和随机噪声强、目标信号弱的地震记录的信噪比, 提出了基于奇异值分解的角度域地震波束形成方法(The Angle Domain Seismic Beam-forming method based on Singular Value Decomposition, ADSBSVD)。该方法先利用自动追踪SVD方法对记录进行预处理, 压制记录中相干噪声。然后利用特征值法估算有效信号信噪比, 当信号信噪比满足SVD适用条件时, 利用SVD直接压制随机噪声, 获得高信噪比地震记录; 当信号信噪比不满足SVD适用条件时, 利用角度域地震波束形成方法增强信号的能量和相关性, 再利用SVD方法压制随机噪声, 获得高信噪比地震记录。模拟数据和实际数据处理结果均验证了ADSBSVD方法的有效性和可行性。

1 方法原理

对于相干噪声和随机噪声强、目标信号弱的地震记录, 采用ADSBSVD方法时, 首先利用SVD方法压制记录中相干噪声, 再利用角度域地震波束形成方法增强目标信号的能量和相干性, 最后再由SVD方法重构相干信号压制记录中随机噪声。

1.1 角度域地震波束形成原理

角度域地震波束形成以传统地震单波束为基础, 首先将目标加强区域划分成若干个角度区间, 通过改变延时参数, 改变主波束方向, 在每个角度区间内均形成地震单波束, 然后提取每个角度区间内被加强的信号替换原来位置处信号, 实现整个目标区域信号能量的增强。角度域地震波束形成的工作原理如图 1所示。

为了清晰展示延时参数的求取过程, 以共接收点道集为例进行演示(图 2)。图 2中震源点等间距排列, 震源点与检波点共线。设震源激发的地震波与震源法向夹角为θ_max, 相邻震源激发的地震波到反射点的波程分别为r_i－1和r_i, 由于震源点距反射点的距离远大于震源间距, 所以近似认为相邻震源激发的地震波平行。由此相邻震源激发地震波的波程差为r_i－1－r_i=dsinθ_max, 则相邻震源间延时参数为:

$ \tau = \frac{{d\sin {\theta _{\max }}}}{v} $

(1)

式中:d为炮间距; θ_max为主波束方向; v为目标层上覆介质平均速度。

设采集的多炮地震记录为R=[R_－n, …, R_－1, R₀, R₁, …, R_n], n为奇数, R₀为中间炮地震记录, 对R中记录依次引入延时－nτ, …, －τ, 0, τ, …, nτ, 按照公式(2)进行叠加, 即形成中间震源处的定向地震记录:

$ {H_0} = \sum\limits_{j = - n}^n {{R_j}\left( {t + j\tau } \right)} $

(2)

式中:－n≤j≤n, j取整数。

在形成角度域地震波束之前, 需要将目标加强区域划分为若干个角度区间, 而决定划分角度区间个数的因素为主波束宽度, 主波束宽度通常用半功率宽度2α_0.5来表示。地震波场方向因子^[8]:

$ F\left( \alpha \right) = \frac{1}{N}\frac{{\sin \left[ {\frac{N}{2}kd\left( {\cos \theta - \cos {\theta _{\max }}} \right)} \right]}}{{\sin \left[ {\frac{1}{2}kd\left( {\cos \theta - \cos {\theta _{\max }}} \right)} \right]}} $

(3)

式中:k为自由空间波数; N为震源叠加个数。

令F(α)=0.707, 根据公式(4)可得半功率点α₁和α₂:

$ \begin{array}{l} {\alpha _1} = \arccos \left( {\frac{{2.78}}{{Nkd}} + \cos {\theta _{\max }}} \right)\\ {\alpha _2} = \arccos \left( { - \frac{{2.78}}{{Nkd}} + \cos {\theta _{\max }}} \right) \end{array} $

(4)

令2α_0.5=|α₁－α₂|, 定义角度区间Δα=2α_0.5。记目标加强区域的起始角度为α_s, 终止角度为α_f, 若α_s和α_f选取满足α_f－α_s大于反射同相轴的宽度, 那么, 欲划分的角度区间个数为:

$ n = \frac{{{\alpha _{\rm{f}}} - {\alpha _{\rm{s}}}}}{{2{\alpha _{0.5}}}} $

(5)

令每个角度区间的中心角为主波束方向, 则第i个角度区间的主波束方向为:

$ {\theta _{\max i}} = {\alpha _{\rm{s}}} + 2\left( {i - \frac{1}{2}} \right){\alpha _{0.5}} $

(6)

式中:1≤i≤n, i为整数。

根据公式(6)和公式(1)计算每个角度区间内的延时参数, 再按照公式(2)依次在每个角度区间内形成地震波束, 即完成角度域地震波束定向的工作。因为在每个主波束方向外均存在畸变波场, 因此, 为了防止不同地震波束所在波场相互干扰, 需要提取每个角度区间内的地震波场信号。

设每个角度区间内地震波场信号的提取函数为I[·], 则提取每个角度区间信号后的地震波场为:

$ H' = I\left[ {{H_{ - n}}} \right] + \cdots + I\left[ {{H_0}} \right] + \cdots + I\left[ {{H_n}} \right] $

(7)

1.2 SVD压制噪声技术

设地震记录为X, 道数为M, 每道的采样点数为N, 一般取M＜N^[13], 那么, 这M×N个数据可以构成二维矩阵, 即:

$ \mathit{\boldsymbol{X}} = \left[ {\begin{array}{*{20}{c}} {{x_{11}}}&{{x_{12}}}& \cdots &{{x_{1N}}}\\ {{x_{21}}}&{{x_{22}}}& \cdots &{{x_{2N}}}\\ \vdots&\vdots &{}& \vdots \\ {{x_{M1}}}&{{x_{M2}}}& \cdots &{{x_{MN}}} \end{array}} \right] $

(8)

X中存在两个正交矩阵U_M×M和V_N×N, 使得:

$ \mathit{\boldsymbol{X}} = {\mathit{\boldsymbol{U}}_{M \times M}}{\mathit{\boldsymbol{E}}_{M \times N}}\mathit{\boldsymbol{V}}_{N \times N}^{\rm{T}} = \sum\limits_{i = 1}^r {{\sigma _i}{\mathit{\boldsymbol{u}}_i}\mathit{\boldsymbol{v}}_i^{\rm{T}}} $

(9)

式中:1≤i≤r, r为X的秩; u_i是XX^T的特征值对应的特征向量; v_i是X^TX的特征值对应的特征向量; E_M×N=diag(σ₁, σ₂, …, σ_i, …, σ_r, 0, …, 0), 其中σ_i为XX^T(或X^TX)的特征值的非负平方根, 即σ_i≥0, 且按递减顺序排列。

依据所使用奇异值的分布范围, 可定义SVD重构的低通滤波器X_LP:

$ {\mathit{\boldsymbol{X}}_{{\rm{LP}}}} = \mathit{\boldsymbol{U}}{\mathit{\boldsymbol{E}}_{1,q}}{\mathit{\boldsymbol{V}}^{\rm{T}}} = \sum\limits_{i = 1}^q {{\sigma _i}{\mathit{\boldsymbol{u}}_i}\mathit{\boldsymbol{v}}_i^{\rm{T}}} $

(10)

式中:q一般选取奇异值分布曲线拐点处的值^[14], 1≤q≤r; E_{1, q}=diag(σ₁, σ₂, …, σ_q, 0, …, 0);X_LP代表相关性强的地震信号。按照公式(10)重构地震记录, 即可实现地震波场分离。

利用公式(10)重构相干噪声, 再从原记录中减去重构相干噪声, 即可实现相干噪声的压制。利用公式(10)重构具有相关性的目标信号, 即可实现随机噪声的压制。

1.3 方法与步骤

基于奇异值分解的角度域地震波束形成方法的实现过程如图 3所示, 其主要步骤如下。

1) 压制地震记录中的相干噪声。由于相干噪声能量强, 常将信号淹没, 所以首先利用SVD压制记录中相干噪声。因相干噪声和信号都具有线性相干性, 为了避免SVD压制相干噪声时破坏有效信号, 先采用自动追踪线性干扰斜率的方法^[17]提取包含相干噪声的小窗口数据, 再利用SVD压制相干噪声。

2) 利用特征值法对压制相干噪声后的记录进行信噪比分析。特征值法估算信噪比R_SN的计算公式^[20]:

$ {R_{{\rm{SN}}}} = \frac{{\lambda _1^2 - \frac{1}{{r - 1}}\sum\limits_{i = 2}^r {\lambda _i^2} }}{{\frac{r}{{r - 1}}\sum\limits_{i = 2}^r {\lambda _i^2} }} $

(11)

式中:r为奇异值个数; λ_i为奇异值。是否需要进行角度域地震波束定向, 取决于SVD方法是否仍然适合压制随机噪声。本文利用不同信噪比地震数据进行SVD实验, 得出决定SVD是否仍然适用于压制随机噪声的因素, 即信号信噪比的临界点R_SN0, 如图 4所示。图 4中, S₁为奇异值分布曲线拐点之前的奇异值之和, 代表具有相关性的有效信号; S₂为奇异值分布曲线拐点之后的奇异值之和, 代表不具有相关性的随机噪声。根据实验得出结论:当信号R_SN>R_SN0时, 可直接利用SVD压制随机噪声; 当R_SN≤R_SN0时, 需要进行角度域地震波束定向, 再利用SVD压制随机噪声。

3) 当信号R_SN>R_SN0时, 利用SVD直接压制随机噪声; 当信号的信噪比R_SN≤R_SN0, 对地震记录进行角度域地震波束定向。

4) 角度域地震波束定向后, 利用相关检测法计算有效信号间时移量, 将弯曲同相轴校平, 再利用SVD压制随机噪声, 最后将校平后的同相轴还原回弯曲同相轴。考虑到全局SVD会造成信号损失, 必要时采用小窗口进行校平和SVD处理。

2 仿真实验

为了验证ADSBSVD方法在提高信噪比方面的作用, 本文采用如图 5所示的水平层状介质速度模型进行数值模拟。第1层介质速度为1000m/s, 厚度为300m, 第2层介质速度为1200m/s, 厚度为700m。数值模拟的信号主频为45Hz, 炮间距为20m, 道间距为20m, 道数为150, 采用中间激发, 两边接收的观测系统。通过射线追踪模拟得到包含两组反射波信号的模拟地震记录如图 6a所示, 采样间隔为0.5ms, 记录时长为3.0s。在图 6a记录中添加相干噪声和高斯白噪声后的单炮记录如图 6b所示。

提取图 6b中相干噪声, 根据自动追踪线性干扰波斜率方法得到相干噪声的斜率为p=0.053。以p=0.053为基准, 上下截取Δt=0.1s长度的小窗口数据, 顺序存放, 得到水平相干噪声(图 7a)。对图 7a按照公式(9)进行SVD。图 7b为图 7a对应的奇异值谱, 令q=1, 按照公式(10)提取相干噪声, 再从图 6b中减去所提取的相干噪声, 即可得到压制相干噪声后的单炮记录, 如图 7c所示。

去除相干噪声后, 分析两层有效信号信噪比。根据特征值法计算出第1层有效信号的R_SN=0.3506dB, R_SN>R_SN0, 因此, 直接用SVD方法提取第1层有效信号。由于记录中存在两层反射信号, 为了避免利用传统的全局SVD提取第1层有效信号时损坏第2层有效信号, 采用小窗口SVD技术。针对第1层反射信号, 选取小窗口的尺寸为1800ms×150道, 利用相关检测法计算信号间时移量, 利用时移量校平第1层有效信号如图 8a所示。图 8b为校平后第1层有效信号对应的奇异值谱。令q=1, 按照公式(10)提取第1层有效信号如图 8c所示。图 8d为按照时移量还原后的第1层有效信号。

针对第2层反射信号, 选取小窗口的尺寸为1200ms×150道, 根据特征值法估算出第2层有效信号的R_SN=0.0136dB, 因R_SN＜R_SN0, 需要对第2层有效信号进行角度域波束定向处理。令N=9, 利用自动追踪SVD方法对9个单炮记录中相干噪声进行压制。利用公式(4)计算波束宽度2α_0.5=59.6°≈60°。令α_f－α_s=180°, 则角度区间个数n=3, 从而将反射同相轴所在空间划分为3个角度区间, 分别为(-90°, -30°), (-30°, 30°), (30°, 90°), 按照公式(6)计算出3个主波束方向θ_max=－60°, 0, 60°, 则3个主波束对应的延时参数分别为-14.5, 0, 14.5ms。然后按照公式(2)对9炮地震记录进行延时、叠加, 得到3个定向地震记录, 如图 9所示。

利用小窗口提取每个角度区间内有效信号替换原位置处有效信号, 得到图 10a记录, 计算其信噪比R_SN=0.3552dB, R_SN>R_SN0, 则可利用SVD提取第2层有效信号。首先利用相关检测法计算信号间时移量, 利用时移量对第2层有效信号校平, 图 10b为校平后第2层有效信号对应的奇异值谱, 令q=1重构第2层有效信号, 图 10c为按照时移量还原后的第2层有效信号。

经过ADSBSVD方法处理后的地震记录如图 11所示。由图 11可知, 处理后的地震记录不仅提高了有效信号的能量和相关性, 而且明显压制了相干噪声和随机噪声, 提高了整个地震记录的信噪比。根据图 8b和图 10b的奇异值谱, 按照公式(11)计算出两层有效信号信噪比分别改善14.20dB和14.18dB。

3 实测信号分析

为了进一步验证ADSBSVD方法的可行性, 利用ADSBSVD方法对野外采集的实际数据进行处理。图 12a是兴城地区一个相干噪声和随机噪声强、有效信号弱的地震记录, 道间距10m, 炮间距50m, 有效信号频率范围25~100Hz, 反射层上覆介质的平均速度为2500m/s。在1.0s左右可以看到一个模糊的反射同相轴。首先利用自动追踪SVD方法压制记录中相干噪声, 结果如图 12b所示。压制相干噪声后, 利用特征值法对记录中有效信号进行信噪比分析, 经估算, 有效信号信噪比R_SN=0.0638dB, R_SN＜R_SN0, 需对该地震记录进行角度域波束定向处理。令N=7, 利用公式(4)计算的波束宽度2α_0.5≈42.6°, 令α_f－α_s=90°, 计算角度区间个数n=2, 则反射同相轴所在空间划分为2个角度区间, 分别为(0, 45°), (45°, 90°)。两个主波束方向为θ_max=22.5°, 67.5°, 依据公式(1)计算的延时参数分别为7.65ms, 18.48ms, 然后按照公式(2)进行7元角度域波束定向, 提取不同主波束方向的有效信号后得到定向地震记录(图 12c), 再利用小窗口SVD方法提取定向地震记录中有效信号, 结果如图 12d所示。

为了更清楚地认识处理前、后记录的变化, 选取了3个小窗口数据H1, H2和H3分别定量计算相干噪声和随机噪声的压制程度以及有效信号信噪比的提高程度。H1和H3分别代表相干噪声和随机噪声, H2代表有效信号。利用处理后窗口内能量与处理前窗口内能量的百分比定义相干噪声和随机噪声压制程度, 经计算, 相干噪声压制了92%, 随机噪声压制了86%。图 13a和图 13b分别为图 12b和图 12d中H2矩形框内的数据放大显示(道数:148~189, 时间:0.9~1.1s)。从图 13中可以看出, 有效信号的连续性和信噪比都得到明显的提高。根据反射信号对应的奇异值谱, 按照公式(11)利用特征值法估算信噪比, 经计算, 处理后H2内数据的信噪比比处理前的提高了3.62dB。由此从实践角度证明了ADSBSVD方法的有效性。

4 结论

本文提出并实现了基于奇异值分解的角度域地震波束形成方法, 该方法通过角度域波束定向实现了目标区域整个反射同相轴能量的增强, 通过SVD方法实现了定向前相干噪声的预处理和定向后随机噪声的压制, 提高了地震记录的信噪比。根据模型试验和实测信号分析得出如下结论。

1) 角度域地震波束形成方法克服了传统单波束不能覆盖整个反射同相轴的缺陷, SVD方法的应用又解决了角度域地震波束形成方法噪声压制不均匀的问题。因此, 针对人工震源或药量较小的炸药震源激发的地震记录, ADSBSVD方法比传统地震波束形成方法和SVD方法应用效果更佳。

2) ADSBSVD方法适用于噪声强, 信号弱, 定向时所需叠加记录个数较大, 并且信号所在同相轴排列较长的情况。

3) 当地震记录中含有多个交叉反射波组时, 为了避免校正工作和奇异值分解过程中反射信号间相互干扰, 仍需进一步研究小窗口的截取方式。