随着“两宽一高”地震勘探的推进, 单点高密度地震采集技术越来越受到重视。刘欣欣等^[1]详细论述了国内外单点高密度地震勘探技术的发展。钱荣钧^[2]指出, 高密度、均匀空间采样技术既是解决当前地震勘探所面对的复杂问题的技术对策, 也是今后地震勘探的发展趋势。李燕燕等^[3]指出, 全方位高密度、超多道数接收仍是不变的地震采集技术发展方向, 预期2020年后形成全方位高密度技术、超多道数固定排列采集技术。曲寿利^[4]结合国内外应用实例分析了高密度采集技术的关键问题, 指出高密度三维地震技术是老油区二次勘探的关键技术之一。正是由于单点高密度采集技术的重要作用, 多年来国内学者对单点高密度数据特点进行了详细分析。曹务祥^[5]、杨贵祥等^[6]、杨照海等^[7]、CAI等^[8]、冯刚等^[9]研究了数字检波器单点采集地震资料的特点, 并结合实际数据分析了不同组合间距对地震信号的影响, 同时指出接收道组合对资料前期的预处理十分必要。于世焕等^[10]对比分析了数字检波器单点采集与组合接收野外试验资料, 发现单点采集有利于保护有效波频率成分, 但会降低地震数据的信噪比。张永刚等^[11]、胡莲莲等^[12]结合数值模拟和实测数据, 从分辨率、信噪比和数据处理的角度剖析了单点高密度地震数据与常规地震数据的区别。近年来单点高密度采集技术在川东碳酸盐岩^[13]、塔里木盆地缝洞型碳酸盐岩^[14]、辽河青龙台^[15]等地区发挥了重要作用。

目前国内针对单点高密度地震数据的室内组合, 主要利用数据一致性处理技术消除组合数据间的差异, 然后采用算术平均法进行组合叠加。曹务祥^[5]提出对地震资料进行严格的去相干噪声、道间振幅和相位的差异调整以及道间时差的校正等处理后, 再进行接收道组合的方法。孙成禹等^[16]利用自适应匹配滤波一致性校正方法对组合数据进行一致性校正, 改善了单点高密度地震数据室内组合的效果。刘志鹏等^[17]利用海上高密度地震数据对比分析了相位谱校正组合与自适应匹配组合方法。在对地震数据进行去噪、静校正、一致性校正等预处理后, 虽然尽可能地消除了组合内波形间的差异, 但仍然无法达到完全一致, 此时采用算术平均法进行组合, 会导致高频信号的损失。

本文通过引入波束形成理论, 分析了基于窗函数的波束形成算法和基于PCA的鲁棒自适应MVDR波束形成算法, 并将其应用于点高密度地震数据的室内组合中, 拓宽了室内组合方法的研究思路, 并且利用Marmousi2模型正演模拟数据对本文组合方法进行了应用测试。

1 方法原理

对单点高密度采集的地震数据, 组合地震道经过一致性处理后, 采用算术平均法合成数据, 如图 1所示。地震道的组合可以视为地震信号经过波束形成器的过程, 其输出信号为^[18-20]:

$ \mathit{\boldsymbol{y}}\left( t \right) = {\mathit{\boldsymbol{W}}^{\rm{T}}}\mathit{\boldsymbol{X}}\left( t \right) $

(1)

其中, X(t)=[x₁(t)  …  x_N(t)]^T是N×1的矩阵, 表征在t时刻N道需要组合的实际数据, 满足:X(t)=S(t)+I(t), S(t)=[s₁(t)  …  s_N(t)]^T是N×1的有效信号矩阵, I(t)=[i₁(t)  …  i_N(t)]^T是N×1的噪声矩阵。W^T=[w₁  w₂  …  w_N]称为权重矢量, 上标T表示矩阵的转置。

由公式(1)可以发现, 地震道的组合受权重矢量的影响。在数据处理中, 各种窗函数的应用十分广泛, 将窗函数引入权重矢量可以拓宽理论研究组合的图 1常规波束形成器示意思路。目前常用的窗函数有均匀窗、余弦窗、汉明窗和布莱克曼窗^[21]等, 分别如下:

1) 均匀加权。

$ {w_n} = \frac{1}{N} $

(2)

其中, n=1, 2, …, N。

2) 余弦加权。

$ {w_{\tilde n}} = \sin \left( {\frac{{\rm{ \mathit{ π} }}}{{2N}}} \right)\cos \left( {{\rm{ \mathit{ π} }}\frac{{\tilde n}}{N}} \right) $

(3)

其中, $\tilde n = n - N/2$ ; n=1, 2, …, N。

3) 汉明加权。

$ {w_{\tilde n}} = 0.54 + 0.46\cos \left( {\frac{{2{\rm{ \mathit{ π} }}\tilde n}}{N}} \right) $

(4)

4) 布莱克曼加权。

$ {w_{\tilde n}} = 0.42 + 0.5\cos \left( {\frac{{2{\rm{ \mathit{ π} }}\tilde n}}{N}} \right) + 0.08\cos \left( {\frac{{4{\rm{ \mathit{ π} }}\tilde n}}{N}} \right) $

(5)

分别将公式(2)~公式(5)代入公式(1)中, 可得到不同窗函数加权下的输出信号。

自适应MVDR波束形成算法能在有效保障信号保真度的前提下提高地震数据的信噪比, 建立二者之间的平衡关系, 但需已知较为准确的标准道。考虑到单点高密度地震数据一致性校正处理后, 已尽可能地消除了组合内波形间的差异, 而PCA可以得到组合数据的主成分, 其第一主成分反映了有效信号的信息, 因此, 可以利用第一主成分重构地震数据, 以此作为组合数据的标准道。针对单点高密度地震数据的特点, 本文引入自适应MVDR波束形成算法, 进一步利用PCA和基于线性组合的对角载入算法对其进行改进, 构建鲁棒自适应MVDR加权函数。

1.1 MVDR波束形成器

利用公式(1)构建自适应MVDR滤波器, 即:在保证信号不失真的情况下, 使得信噪比最大^[22-23]。

$ \mathop {\min }\limits_\mathit{\boldsymbol{W}} {\mathit{\boldsymbol{W}}^{\rm{T}}}{\mathit{\boldsymbol{R}}_i}\mathit{\boldsymbol{W}}\;且\;{\mathit{\boldsymbol{W}}^{\rm{T}}}{\mathit{\boldsymbol{R}}_s}\mathit{\boldsymbol{W}} = 1 $

(6)

其中, R_s=E{S(t)S(t)^T}和R_i=E{I(t)I(t)^T}分别是有效信号的协方差矩阵和噪声的协方差矩阵。

利用拉格朗日乘数法构建目标函数:

$ H\left( {\mathit{\boldsymbol{W}},\lambda } \right) = {\mathit{\boldsymbol{W}}^{\rm{T}}}{\mathit{\boldsymbol{R}}_i}\mathit{\boldsymbol{W}} + \lambda \left( {1 - {\mathit{\boldsymbol{W}}^{\rm{T}}}{\mathit{\boldsymbol{R}}_s}\mathit{\boldsymbol{W}}} \right) $

(7)

其中, λ为拉格朗日乘数。利用(7)式对W^T求偏导并令其为零, 可得到:

$ {\mathit{\boldsymbol{R}}_i}\mathit{\boldsymbol{W}} = \lambda {\mathit{\boldsymbol{R}}_s}\mathit{\boldsymbol{W}} $

(8)

由于矩阵R_i和R_s均是半正定矩阵, 因此λ是一个非负的实数, 由公式(8)可得出W^TR_iW=λW^TR_sW。如果 $\mathop {\min }\limits_\mathit{\boldsymbol{W}} {\mathit{\boldsymbol{W}}^{\rm{T}}}{\mathit{\boldsymbol{R}}_i}\mathit{\boldsymbol{W}}$ 成立, 那么所求出的λ值也应是最小的。

对公式(8)进行简单变换, 即:

$ \mathit{\boldsymbol{R}}_i^{ - 1}{\mathit{\boldsymbol{R}}_s}\mathit{\boldsymbol{W}} = \frac{1}{\lambda }\mathit{\boldsymbol{W}} $

(9)

由公式(9)可以发现, 需要求取的最优权重矢量W_opt即为矩阵R_i^－1R_s的最大特征值1/λ_min对应的特征向量W_max, 即:

$ {\mathit{\boldsymbol{W}}_{\max }} = P\left\{ {\mathit{\boldsymbol{R}}_i^{ - 1}{\mathit{\boldsymbol{R}}_s}} \right\} $

(10)

在实际应用中, 噪声的协方差矩阵R_i难以获取, 当利用实际数据的协方差矩阵R_x=E{X(t)X(t)^T}=R_i+R_s代替R_i时, (10)式的最优解不会发生改变, 此时:

$ {\mathit{\boldsymbol{W}}_{\max }} = P\left\{ {\mathit{\boldsymbol{R}}_x^{ - 1}{\mathit{\boldsymbol{R}}_s}} \right\} $

(11)

由公式(11)可以发现, 要求得权重矢量W_max需已知有效信号的协方差矩阵R_s, 因此, 需要利用PCA提取组合信号的有效信息^[24-25]。具体实现步骤如下:

1) 设X_N×M表示N道地震数据, 每道地震数据的点数为M, 对X_N×M进行归一化(零均值或者去均值处理), 得到δX_N×M;

2) 计算δX_N×M的协方差矩阵D_N×M;

3) 对协方差矩阵进行特征值分解, 即D·V=Λ·V。其中, Λ为特征值矩阵, Λ=diag(λ₁, λ₂, …, λ_N), 且λ₁≥λ₂≥…≥λ_N≥0;V是N×N的特征向量矩阵, 且VV^T=1;

4) 计算δX_N×M在V中的投影, 即Y=V^TδX, 得到δX_N×M矩阵的主成分y₁, y₂, …, y_N;

5) 选择代表有效信号的主成分并使用公式 $\mathit{\boldsymbol{X'}} = \sum\limits_{i = 1}^r {{\mathit{\boldsymbol{v}}_i}{\mathit{\boldsymbol{y}}_i}} $ , r≤N进行重构, 利用重构数据作为组合数据的有效信号。

1.2 基于对角载入的鲁棒MVDR自适应波束形成算法

在实际应用中, 准确的协方差矩阵R_x和R_s无法获得, 但可以估计得到, 因此有效信号协方差矩阵和信号协方差矩阵可以近似为^[22]:

$ {{\mathit{\boldsymbol{\hat R}}}_s} = \frac{1}{M}\sum\limits_{k = 1}^M {\mathit{\boldsymbol{S}}\left( k \right)\mathit{\boldsymbol{S}}{{\left( k \right)}^{\rm{T}}}} $

(12)

$ {{\mathit{\boldsymbol{\hat R}}}_x} = \frac{1}{M}\sum\limits_{k = 1}^M {\mathit{\boldsymbol{X}}\left( k \right)\mathit{\boldsymbol{X}}{{\left( k \right)}^{\rm{T}}}} $

(13)

式中:M表示每道数据的采样点数。

此时, 公式(11)变为:

$ {\mathit{\boldsymbol{W}}_{\max }} = P\left\{ {\mathit{\boldsymbol{\hat R}}_x^{ - 1}{{\mathit{\boldsymbol{\hat R}}}_s}} \right\} $

(14)

但是, 实际数据中M≫N, 使得近似 ${{\mathit{\boldsymbol{\hat R}}}_x}$ 和 ${{\mathit{\boldsymbol{\hat R}}}_s}$ 常为病态矩阵且rank( ${{\mathit{\boldsymbol{\hat R}}}_s}$ )≠1。为了更加准确地估计协方差矩阵R_x和R_s且增强算法的鲁棒性, 我们利用PCA算法重构数据获取 ${{\mathit{\boldsymbol{\hat R}}}_s}$ , 并引入基于线性组合的对角载入算法对 ${{\mathit{\boldsymbol{\hat R}}}_x}$ 和 ${{\mathit{\boldsymbol{\hat R}}}_s}$ 进行修正, 使得修正后矩阵 ${{\mathit{\boldsymbol{\tilde R}}}_x}$ 和 ${{\mathit{\boldsymbol{\tilde R}}}_s}$ 在满秩的情况下逼近R_x和R_s, 再利用公式(14)求取最优权重。

对 ${{\mathit{\boldsymbol{\hat R}}}_x}$ 的一个线性收缩估计为^[26]:

$ {{\mathit{\boldsymbol{\tilde R}}}_x} = \alpha \mathit{\boldsymbol{J}} + \beta {{\mathit{\boldsymbol{\hat R}}}_x} $

(15)

式中: ${{\mathit{\boldsymbol{\tilde R}}}_x}$ 是 ${{\mathit{\boldsymbol{\hat R}}}_x}$ 的递进估计; J是单位矩阵; α∈(0, 1), β∈(0, 1)是收敛系数。

为使估计矩阵 ${{\mathit{\boldsymbol{\tilde R}}}_x}$ 尽量接近矩阵R_x, 利用 ${{\mathit{\boldsymbol{\tilde R}}}_x}$ 与R_x的最小均方误差求取α, β, 即:

$ {\rm{MSE}}\left( {{{\mathit{\boldsymbol{\tilde R}}}_x}} \right) = {\rm{E}}\left\{ {{{\left\| {{{\mathit{\boldsymbol{\hat R}}}_x} - {\mathit{\boldsymbol{R}}_x}} \right\|}^2}} \right\} $

(16)

式中:‖·‖表征矢量范数。

将公式(15)代入公式(16)中并分别对α, β求偏导, 经过推导得到最优解α₀, β₀^[26]:

$ \left\{ \begin{array}{l} {\alpha _0} = \frac{{\left( {1 - {\beta _0}} \right){\rm{tr}}\left( {{\mathit{\boldsymbol{J}}^{\rm{T}}}{\mathit{\boldsymbol{R}}_x}} \right)}}{{{{\left\| \mathit{\boldsymbol{J}} \right\|}^2}}}\\ {\beta _0} = \frac{\gamma }{{\rho + \gamma }} \end{array} \right. $

(17)

其中, $\rho = {\rm{E}} \left\{ {{{\left\| {{{\mathit{\boldsymbol{\hat R}}}_x} - {\mathit{\boldsymbol{R}}_x}} \right\|}^2}} \right\}$ , γ=[‖R_x‖²‖J‖²－tr²(JR_x)]/‖J‖²。

由于准确的协方差矩阵R_x无法获得, 但可以利用估计值 ${{\mathit{\boldsymbol{\hat R}}}_x}$ 进行替换, 因此得到估计值 ${{\hat \alpha }_0}$ , ${{\hat \beta }_0}$ :

$ \left\{ \begin{array}{l} {{\hat \alpha }_0} = \hat v\left( {1 - {{\hat \beta }_0}} \right)\\ {{\hat \beta }_0} = \frac{{\hat \gamma }}{{\hat \rho + \hat \gamma }} \end{array} \right. $

(18)

其中, $\hat \rho = 1/{N^2}\sum\limits_{n = 1}^N {{{\left\| {\mathit{\boldsymbol{X}}\left( n \right)\mathit{\boldsymbol{X}}{{\left( n \right)}^{\rm{T}}} - {{\mathit{\boldsymbol{\hat R}}}_x}} \right\|}^2}} $ , $\hat v = \left[ {{\rm{tr}}\left( {{\mathit{\boldsymbol{J}}^{\rm{T}}}{{\mathit{\boldsymbol{\hat R}}}_x}} \right)} \right]/{\left\| \mathit{\boldsymbol{J}} \right\|^2}$ , $\hat \gamma = {\left\| {\hat v\mathit{\boldsymbol{J}} - {{\mathit{\boldsymbol{\hat R}}}_x}} \right\|^2}$ 。

此时, 协方差矩阵R_x由 ${{\mathit{\boldsymbol{\tilde R}}}_x}$ 代替^[26]:

$ {{\mathit{\boldsymbol{\tilde R}}}_x} = {{\hat \alpha }_0}\mathit{\boldsymbol{J}} + {{\hat \beta }_0}{{\mathit{\boldsymbol{\hat R}}}_x} $

(19)

同理可以得到有效信号的最优估计 ${{\mathit{\boldsymbol{\tilde R}}}_s}$ , 这里不再赘述。

此时, 公式(14)变为:

$ {\mathit{\boldsymbol{W}}_{\max }} = P\left\{ {\mathit{\boldsymbol{\tilde R}}_x^{ - 1}{{\mathit{\boldsymbol{\tilde R}}}_s}} \right\} $

(20)

将公式(20)代入公式(1)中, 可以得到经过鲁棒自适应MVDR波束形成器后的输出信号。

2 数值分析

为了进一步对比分析不同加权组合情况下地震记录的特点, 分别从组合后输出数据的信噪比、幅频特性、同相轴连续性等几个方面进行分析比较。

根据地震波传播的特征, 选取雷克子波作为有效信号。雷克子波的表达式如下^[27]:

$ f\left( t \right) = \left( {1 - 2{{\rm{ \mathit{ π} }}^2}f_p^2{t^2}} \right)\exp \left( { - {{\rm{ \mathit{ π} }}^2}f_p^2{t^2}} \right) $

(21)

式中:t表示时间; f_p为子波主频。

为了分析不同加权组合对地震信号信噪比的影响, 利用公式(21)设计10道主频f_p=50Hz的地震信号S(t), 每道地震信号加入随机噪声I(t), 则第n道输入信号的信噪比为^[21]:

$ {R_{{\rm{SNin}}}} = \frac{{{\mathit{\boldsymbol{R}}_s}}}{{{\mathit{\boldsymbol{R}}_i}}} $

(22)

有效信号协方差和噪声协方差近似为 ${\mathit{\boldsymbol{R}}_s} = \left( {1/M} \right)\sum\limits_{k = 1}^M {{s_n}\left( k \right){s_n}{{\left( k \right)}^{\rm{T}}}} $ 和 ${\mathit{\boldsymbol{R}}_i} = \left( {1/M} \right)\sum\limits_{k = 1}^M {{i_n}\left( k \right){i_n}{{\left( k \right)}^{\rm{T}}}} $ , n取值范围为0~10, 表示第n道检波器。其中, S(t)=[s₁(t)  …  s₁₀(t)]^T和I(t)=[i₁(t)  …  i₁₀(t)]^T分别表示有效信号矢量和噪声矢量。

当采用加权组合合成地震信号时, 组合后输出信号的信噪比满足^[22]:

$ {R_{{\rm{SNout}}}} = \frac{{{\mathit{\boldsymbol{W}}^{\rm{T}}}{\mathit{\boldsymbol{R}}_s}\mathit{\boldsymbol{W}}}}{{{\mathit{\boldsymbol{W}}^{\rm{T}}}{\mathit{\boldsymbol{R}}_i}\mathit{\boldsymbol{W}}}} $

(23)

其中, ${\mathit{\boldsymbol{R}}_s} = {\rm{E}}\{ \mathit{\boldsymbol{S}}\left( t \right)\mathit{\boldsymbol{S}}{\left( t \right)^{\rm{T}}}\} \approx \left( {1/M} \right)\sum\limits_{k = 1}^M {\mathit{\boldsymbol{S}}\left( k \right)\mathit{\boldsymbol{S}}{{\left( k \right)}^{\rm{T}}}} $ , ${\mathit{\boldsymbol{R}}_i} = {\rm{E}}\{ \mathit{\boldsymbol{I}}\left( t \right)\mathit{\boldsymbol{I}}{\left( t \right)^{\rm{T}}}\} \approx \left( {1/M} \right)\sum\limits_{k = 1}^M {\mathit{\boldsymbol{I}}\left( k \right)\mathit{\boldsymbol{I}}{{\left( k \right)}^{\rm{T}}}} $ 分别为N×N的有效信号协方差矩阵和噪声协方差矩阵。

图 2为10道含随机噪声且不含时差的组合数据, 表 1为不同加权接收对该组合数据信噪比的改善情况。图 3为10道含随机噪声和时差的组合数据, 其中第2道和第8道有较小的时差, 表 2为不同加权接收对该组合数据信噪比的改善情况。可以发现加权组合可以提高地震数据的信噪比:①在信号不存在时差的情况下, 各加权算法对信噪比的改善差异较小; ②当信号中存在时差时, 相比于其它几种加权组合, 基于PCA的鲁棒自适应MVDR加权算法能进一步提高信号的信噪比。

为了分析不同加权对地震数据幅频特性及地震记录同相轴连续性的影响, 在不考虑噪声的情况下, 构建含气的Marmousi2模型, 如图 4所示。模型大小为1001×451, 网格大小为3m×3m, 时间采样间隔为1ms, 记录时间为2s, 激发主频为30Hz, 炮点在第470道距地面4个采样点(12m)处。对模型进行正演模拟, 得到炮记录如图 5所示。分别对炮记录进行均匀加权组合、余弦加权组合、汉明加权组合、布莱克曼加权组合和基于PCA的鲁棒自适应MVDR加权组合并显示在同一坐标系下, 其正演记录如图 6所示。

由图 6可知, 相比于常规的均匀加权组合, 基于其它窗函数的加权组合和基于PCA的鲁棒自适应MVDR加权组合均能使炮记录上同相轴的能量更加均衡, 更易识别追踪, 其中基于PCA的鲁棒自适应MVDR加权组合效果更为明显。同时, 抽取第90道地震记录进行频谱分析(图 7)发现, 相比于均匀加权, 其它几种加权算法均能更好地保护中高频成分, 取得较好的保真效果。

为了进一步对比分析各加权组合对信号的影响和它们之间的差异, 对图 6加权组合正演的地震记录进行F-K变换并计算各加权地震记录F-K变换结果与均匀加权变换结果的差值(图 8~图 12)。由图 8~图 12可以发现:与常规均匀加权相比, 基于窗函数的加权组合能保护部分高频信号; 基于PCA的鲁棒自适应MVDR加权组合对地震记录各频率均有一定的保护作用, 图 12a中20~40Hz频段范围内信号的保真效果更为明显, 而进行正演模拟时采用的子波主频为30Hz, 体现了MVDR波束形成器的优势, 公式(6)中W^TR_sW=1保证了有效信号在不失真的情况下进行组合, 同时, 对比差值剖面可以发现, 基于PCA的鲁棒自适应MVDR加权组合对低波数信号的保护更加明显。由于地层速度随着深度的增加呈增大趋势, 因此在频率一定的情况下, 深层信号的波数越小。可见, 相比于其它组合方式, 基于PCA的鲁棒自适应MVDR加权组合对有效信号乃至深层信号能起到更为显著的保护作用。

3 结论与认识

本文重点研究了基于窗函数的波束形成算法和基于PCA的鲁棒自适应MVDR波束形成算法, 将其与目前常用的均匀加权(算术平均)算法对比, 得到以下认识:①组合能提高地震数据的信噪比, 余弦加权、汉明加权、布莱克曼加权能进一步保护地震数据的高频成分; ②基于PCA的鲁棒自适应MVDR加权组合实现了信号保真与信噪比的有效统一, 对深层信号的保护作用较为明显。

针对单点高密度地震数据的特点, 对地震数据的去噪、静校正和一致性校正等预处理很有必要, 基于波束形成理论的数字组合方法能够对当前单点高密度地震数据的室内组合处理技术进行适当的补充, 下一步将研究该方法在单点高密度实际资料处理中的应用。