可控震源具有安全、环保、高效、低能耗的特点, 这使得可控震源采集技术得到了广泛应用。国内可控震源采集技术的利用率逐年提高, 尤其是在西部大沙漠地区。随着可控震源采集技术的不断发展和应用, 可控震源已成为油气勘探领域的研究热点^[1-6]。地球物理学家们在确保施工效率和采集数据质量的前提下, 发展了一系列高效率的采集方法, 如交替扫描、滑动扫描、独立同步扫描(independent simultaneous sweeping, ISS)、远距离同步扫描(distance separated simultaneous sweeping, DSSS)、V1扫描、远距离同步滑动扫描(distance separated simultaneous and slip sweeping, 简称DS4)等。高保真可控震源地震(high fidelity vibratory seismic, HFVS)采集因采用多组经过相位编码后的震源在不同炮点处同时激发, 故优于常规采集, 具有“高效”的特点。相对于高效采集, 高保真采集得到的地震数据具有良好的保真度。广义的混叠地震数据也称超道集地震数据, 是采用多震源组合激发获得的混合地震波场, 能否对混叠地震数据进行有效的分离是后续地震数据处理的关键^[7]。

SILVERMAN^[8]最早提出震源同时激发的概念。随后, BEASLEY等^[9]将组合激发的概念推广至炸药震源, 并引起了勘探地球物理学界的关注。针对可控震源高效采集的弊端, SALLAS等^[10]首次提出了对地面力信号进行反褶积处理获得可控震源高保真地震记录的方法。随后, ALLEN等^[11]介绍了HFVS采集方法的基本原理和现场试验结果。KROHN等^[12]对此进行深入研究并提出了基于多台可控震源的高保真数据采集方法。BERKHOUT^[13]从理论上详细讨论了组合震源的数学表征, 并基于同步激发理论提出了随机时延激发的混叠地震采集方法。

目前主要有两种方法实现混叠地震数据成像:第1种是先进行混叠地震数据分离^[14-17]再成像; 第2种是不进行混叠地震数据分离而直接成像的混叠地震数据直接成像方法, 在此过程中采用编码或者约束方法消除串扰噪声的影响^[18-21]。目前已有的混叠地震数据分离方法主要包括两类:①将混叠地震数据中的邻炮干扰当作噪声, 采用滤波方法对其进行压制^[22-23], 此类方法通常针对人工分选道集如共接收点道集、共中心点道集及共偏移距道集^[13]; ②将混叠地震数据的分离转化为反演问题, 由于混叠方程通常为病态方程, 所以需要对其施加约束条件进行求解, 此类反演问题可以转化为直接求解正演算子的矩阵^[24]或者采用迭代方法进行迭代相减^[25-26]求解。编码可以减少混叠地震数据之间的相干性, 提高混叠地震数据的信噪比。因可控震源激发子波具有可控性, 故编码尤其适用于可控震源采集得到的地震数据。

近年来, 混叠地震数据直接成像方法因处理流程简单、高效得到了广泛的关注、研究与应用, 在混叠地震数据直接成像方法中, 我们引入了3种压制可控震源混叠地震数据串扰成像噪声的方法, 并将这3种方法联合应用于混叠地震数据直接成像。本文首先介绍了针对多台震源多次激发记录的可控震源混叠地震数据分离方法; 然后介绍了可控震源混叠噪声压制方法; 最后提出了可控震源混叠地震数据直接成像方法, 并将其应用于模型试算和实际地震资料测试。

1 可控震源高保真混叠地震数据分离 1.1 高保真混叠地震数据采集原理

可控震源高保真地震数据采集方法是可控震源野外采集方法之一, 它弥补了高效采集中数据不保真(受扫描信号影响)、谐波干扰严重、初至难以拾取等不足, 是今后高精度可控震源勘探的重要发展方向。利用可控震源高保真地震数据采集方法所获得的野外地震采集资料为多震源混叠波场, 需要对其进行数据分离以获得单炮记录。

高保真可控震源混叠地震数据分离方法只适用于利用高保真可控震源地震数据采集方法得到的地震数据。利用高保真可控震源地震数据采集方法进行地震数据采集时需满足以下条件:①每组震源车间隔距离相同; ②每组震源车激发次数一定不少于震源车的组数, 如震源车共4组, 那么每组震源车激发次数不能少于4次; ③一般情况下, 每组震源车初始激发相位不完全相同; ④每组震源车每次激发的初始相位不同。合理的相位编码矩阵可以提升高保真地震数据的分离效果。

1.2 高保真混叠地震数据的分离原理

以4台震源车多次扫描为例阐述高保真可控震源采集混叠地震数据的分离原理。在野外4台震源车以不同初始相位在不同炮点同时激发(组合震源激发的一种方式), 在检波点处记录混叠波场得到可控震源高保真地震数据, 图 1是HFVS采集示意。

地震波场在观测道集中满足线性叠加原理, 因此多震源波场可视作单个震源激发(或模拟)的地震波场的线性组合, 在波动方程表达式中表现为震源的线性组合。相较于单个震源激发, 多个震源激发可得到如下的地震波场波动方程:

$ {\nabla ^2}u - \frac{1}{{{v^2}}}\frac{{{\partial ^2}u}}{{\partial {t^2}}} = {s_{\rm{s}}}(t)\varGamma ({X_{\rm{s}}})\delta (X - {X_{\rm{s}}}) $

(1)

式中:X=(x, y, z), 为波场某一点的坐标值; X_s为炮点坐标; s_s(t)为各炮点处震源函数; v为速度函数; u为波场; Γ(X_s)为混合算子; t为采样时间, δ(X－X_s)代表脉冲信号。

由组合震源概念^[7]可知, 混合算子在频率域表达式为:

$ \varGamma ({X_{\rm{s}}}) = \alpha ({X_{\rm{s}}}){{\rm{e}}^{{\rm{j}}\omega t({X_{\rm{s}}})}} $

(2)

式中:α(X_s)为高保真采集的权系数, 本文令α(X_s)=1, 混合算子存在初始相位的差异。波动方程虽然精确, 但很难推算出准确的多震源混合算子, 因此本文根据褶积模型理论推导混合算子, 并采用有限差分法实现合成HFVS采集数据的分离, 具体如下:首先将振动记录作为地面力信号与地层脉冲响应的褶积, 然后考虑吸收衰减和激发接收响应等因素的影响, 再利用地面力信号对振动记录进行反褶积, 最终求出包含衰减的地层脉冲响应。

以N台震源车M次扫描为例, 可得到混合波场的数学表达式为:

$ \begin{array}{l} {x_i}(t) = \sum\limits_{j = 1}^N {{g_{ij}}} (t) * {{\rm{e}}_j}(t)\\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} i = 1,2, \cdots ,M \end{array} $

(3)

式中:x_i(t)为第i次扫描的混合波场, 该波场形成的记录由来自第i次扫描中每个震源-检波点路径反射和折射后所有振动信号组成; g_ij(t)为第j组震源车第i次激发扫描时记录的近源力信号(也被称为地面力信号); e_j(t)为第j组震源车在同一点激发的地层脉冲响应, 可以认为震源每次激发的e_j(t)不变。由信号分析理论可知, 时域褶积在频率域表现为乘积形式, (3)式在频域可改写为:

$ \begin{array}{*{20}{l}} {{x_i}(f) = \sum\limits_{j = 1}^N {{g_{ij}}} (f){e_j}(f)}\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} i = 1,2, \cdots ,M} \end{array} $

(4)

(4) 式可用矩阵形式表示为如下方程:

$ \mathit{\boldsymbol{GE}} = \mathit{\boldsymbol{X}} $

(5)

式中: G表示力信号系数矩阵; E表示频域脉冲响应矩阵; X表示频域振动记录矩阵。G、E、X具体可表示为:

$ \begin{array}{l} \mathit{\boldsymbol{G}} = \left[ {\begin{array}{*{20}{c}} {{g_{11}}}&{{g_{12}}}& \cdots &{{g_{1N}}}\\ {{g_{21}}}&{{g_{22}}}& \cdots &{{g_{2N}}}\\ \vdots & \vdots &{}& \vdots \\ {{g_{M1}}}&{{g_{M2}}}& \cdots &{{g_{MN}}} \end{array}} \right]\quad \mathit{\boldsymbol{E}} = \left[ {\begin{array}{*{20}{c}} {{e_1}}\\ {{e_2}}\\ \vdots \\ {{e_N}} \end{array}} \right]\\ \mathit{\boldsymbol{X}} = \left[ {\begin{array}{*{20}{c}} {{x_1}}\\ {{x_2}}\\ \vdots \\ {{x_M}} \end{array}} \right] \end{array} $

(6)

从数学角度分析可知, 根据系数矩阵G的秩与N的大小不同, 方程(5)存在如下3类情况:①当rank(G | X)=rank(G)=N时, 方程存在唯一解; ②当rank(G | X)=rank(G) < N时, 方程有无穷多个解, 求解此类欠定反问题需要事先施加特定的先验约束条件, 如解的L₂范数为最小范数解(也称为解的能量最小约束); ③当rank(G | X)≠rank(G)时, 方程无解, 需要求取某种条件下的近似解。以上3种情况下, E可以统一表示为:

$ \mathit{\boldsymbol{E}} = {\mathit{\boldsymbol{G}}^ + }\mathit{\boldsymbol{X}} $

(7)

式中: G⁺代表G的广义逆。考虑到噪声的影响, 还可以采用基于奇异值分解的方法求解方程(5), 结果分别为:

$ {\mathit{\boldsymbol{G}} = \mathit{\boldsymbol{UD}}{\mathit{\boldsymbol{V}}^{\rm{T}}}} $

(8)

$ {\mathit{\boldsymbol{E}} = \mathit{\boldsymbol{V}}{\mathit{\boldsymbol{D}}^{ - 1}}{\mathit{\boldsymbol{U}}^{\rm{T}}}\mathit{\boldsymbol{X}}} $

(9)

式中: U和V是正交矩阵; D是由降序排列的非负对角线元素组成的对角矩阵。将根据公式(9)求得的每台震源车对应的传输矩阵E与扫描信号褶积, 可以得到对应震源的单炮振动记录。分离混叠地震数据的步骤如下:①对近源力信号和每次扫描获得的混叠地震数据进行傅里叶变换; ②在目标频率域内的每个频率点处利用频率域力信号构造力信号系数矩阵, 并求取广义逆算子; ③将广义逆算子作用于对应频率点处的地震记录, 得到此频率点处的E; ④对所有的频率点重复步骤②和③, 即可求出所有频率点处的E, 进而可分离出每台震源对应的包含衰减的地层脉冲响应; ⑤对步骤④求出的E进行时间域反傅里叶变换, 再将其与扫描信号褶积即可得到分离的单炮记录。

在实际地震勘探中, 由于X存在扰动, G也存在扰动, 这必然引起E的扰动, 由数值分析可知: E的扰动大小与矩阵G的条件数有关。G的条件数Cond(G)为:

$ {\rm{Cond}} (\mathit{\boldsymbol{G}}) = {\left\| {{\mathit{\boldsymbol{G}}^{ - 1}}} \right\|_p}{\left\| \mathit{\boldsymbol{G}} \right\|_p} $

(10)

式中:‖ ‖_p为范数表达式。在野外HFVS采集过程中, 对每个炮点处的扫描信号设置不同的相位, 令扫描次数M不小于震源台数N, 可以保证地面力信号系数矩阵G的条件数较小。上述情况下, 即使炮记录存在扰动也能得到良好的分离结果。但考虑到噪声的影响, 即使M>N也可能方程欠定, 因此实际采集过程中应该进行工区分析试验, 在此基础上, 选择一组较好的相位编码矩阵用于混叠噪声的压制。

1.3 高保真地震数据的试验

对实际采集的高保真地震数据进行分离, 采集参数如下:扫描频率范围5~90 Hz, 扫描时间为16 s, 记录时间为8 s, 4台可控震源车5次扫描, 总接收时间为24 s, 采用如表 1所示的相位编码。图 2为第1次扫描得到的高保真混叠地震记录, 可以看出4炮数据混叠在一起, 互相干扰, 对每个炮记录而言, 其它3炮的结果都可看作混叠噪声。采用混叠地震数据分离方法分离后的高保真单炮记录如图 3所示, 4炮记录均被很好地分离, 而且无谐波噪声的影响。矩阵的条件数越小, 分离结果越接近真实情况, 分离效果越好。本数据在生产中采用4台震源车5次扫描, 因此编码矩阵的条件数小, 分离效果良好。

2 基于地表一致性原理的可控震源混叠噪声压制 2.1 基本原理

混叠噪声是一种炮域相干性强、压制难度极大的可控震源特征噪声, 它属于外来信号, 不符合观测系统对地质目标的激发接收规律, 也不符合地表一致性原理。因混叠噪声在共中心点域无任何相干性, 表现为随机出现的振幅能量较强的噪声, 故可以在共中心点域对其进行压制。如图 4所示, 混叠噪声在炮域表现为强能量的振幅以及强相干性的同步交涉干扰, 转换到共中心点域后相干性变弱。

利用混叠噪声在共中心点域随机分布和强能量的特征, 提出了基于地表一致性原理的混叠噪声压制方法, 并采用随机异常振幅衰减技术对混叠噪声进行压制, 步骤如下。

1) 在特定的数据域(非炮域)从小到大划分时窗, 时窗需有一定的重叠以避免边界效应。

2) 以时窗内的数据为计算单元, 拾取时窗内地震道的振幅属性, 分析其均方根振幅、平均绝对振幅、最大振幅等属性。

3) 将拾取的振幅属性按照地表一致性原理, 分解至炮域、共检波点域、共偏移距域和共中心点域, 然后对各个域的分量进行动态的时空变加权处理, 加权因子的计算公式为:P(i)=S_i·R_i·G_i·M_i, 其中, P(i)表示时空加权因子; S_i, R_i, G_i和M_i分别为该时窗下的炮域、共检波点域、共偏移距域和共中心点域振幅属性。采用矢量中值滤波^[26]或者其它随机噪声压制方法压制随机分布的强能量振幅。

4) 应用加权因子, 将各域地震数据合并成整体地震数据, 完成一个时窗内的数据处理。

5) 重复上述所有步骤, 直到所有时窗内的地震数据全部处理完。

可控震源混叠地震数据压制方法的适用条件为不同时激发多个炮, 激发时间尽可能随机。

2.2 ISS地震数据混叠噪声压制

对ISS地震数据进行混叠噪声压制测试, 可控震源的扫描参数如下:线性升频扫描信号, 频率范围为5~80 Hz, 扫描时间为16 s, 记录时间为8 s。图 5a为混叠噪声压制前, 与扫描信号互相关后的ISS单炮记录, 该记录存在明显的混叠噪声(图 5a中的蓝色箭头), 这些混叠噪声分别来自另外两次扫描。采用基于地表一致性原理的可控震源混叠噪声压制方法对混叠噪声进行压制, 压制后的单炮记录如图 5b所示, 可以看出, 来自另外两次扫描的混叠噪声得到了很好的压制。图 6为从图 5a和图 5b中抽取的第50道的地震波波形, 图中红色实线代表抽取自混叠噪声压制前的ISS单炮记录, 蓝色虚线代表抽取自混叠噪声压制后的ISS单炮记录, 可以看出, 中深层的混叠地震数据能量得到了准确的压制, 从浅层与中、深层混叠地震数据振幅对比可知, 混叠噪声压制后, 有效能量几乎没有任何损失, 说明了基于地表一致性原理的可控震源混叠噪声压制方法可用于压制ISS地震数据的混叠噪声。

3 可控震源混叠地震数据直接成像 3.1 最小二乘逆时偏移

由于可控震源扫描信号可控, 因此相较于炸药震源, 对可控震源混叠地震数据进行最小二乘逆时偏移处理时易出现的子波问题可以得到很好的解决, 在反演过程中, 输入的子波信号可近似为扫描信号自相关后的结果, 通过最小二乘逆时偏移处理在一定程度上可以削弱混叠噪声的影响。最小二乘逆时偏移处理可被视作一种求解最优反演问题的方法, 其目的是在最小二乘意义下寻找扰动模型, 将(11)式视作最小化目标函数来拟合观测数据:

$ \begin{array}{*{20}{l}} {\mathit{\boldsymbol{J}}(\mathit{\boldsymbol{m}}) = \frac{1}{2}\left\| {{\mathit{\boldsymbol{d}}_{{\rm{cal}}}} - {\mathit{\boldsymbol{d}}_{{\rm{obs}}}}} \right\|_2^2}\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} = \frac{1}{2}\left\| {{\mathit{\boldsymbol{d}}_{{\rm{cal}}}} - {\mathit{\boldsymbol{d}}_{{\rm{obs}}}}} \right\|_2^2 \to {\rm{min}}} \end{array} $

(11)

式中: d_cal和d_obs分别表示合成地震数据和观测地震数据; L表示线性正演算子; m表示反射系数模型, d_cal= Lm。一阶声波速度-应力方程为:

$ \left\{ {\begin{array}{*{20}{l}} {\rho \frac{{\partial u}}{{\partial t}} = \frac{{\partial p}}{{\partial x}}}\\ {\rho \frac{{\partial w}}{{\partial t}} = \frac{{\partial p}}{{\partial z}}}\\ {\frac{1}{{{v^2}}}\frac{{\partial p}}{{\partial t}} = \rho \left( {\frac{{\partial u}}{{\partial x}} + \frac{{\partial w}}{{\partial z}}} \right) + \sum\limits_{i = 1}^N {{f_i}} } \end{array}} \right. $

(12)

式中:p是声压波场; u和w分别是水平分量和垂直分量的质点速度波场; ρ为密度; v为速度; x, z分别为水平和垂直方向的空间坐标; t为时间; f_i为震源项。

基于伯恩近似理论, 反偏移算子方程为:

$ \left\{ {\begin{array}{*{20}{l}} {\rho \frac{{\partial \delta u}}{{\partial t}} = \frac{{\partial \delta p}}{{\partial x}}}\\ {\rho \frac{{\partial \delta w}}{{\partial t}} = \frac{{\partial \delta p}}{{\partial z}}}\\ {\frac{1}{{{v^2}}}\frac{{\partial \delta p}}{{\partial t}} = \rho \left( {\frac{{\partial \delta u}}{{\partial x}} + \frac{{\partial \delta w}}{{\partial z}}} \right) + \frac{{m(v)}}{{v_0^2}}\frac{{\partial {p_0}}}{{\partial t}}} \end{array}} \right. $

(13)

式中:p₀为背景声压波场; δu, δw, δp分别为速度扰动产生的扰动质点水平分量和垂直分量速度波场以及扰动声压波场; v₀为背景速度场; m(v)为速度反射系数场。m(v)可表示为:

$ m(v) = \frac{{2\delta v}}{v} $

(14)

根据伴随状态理论得到伴随状态方程为:

$ \left\{ {\begin{array}{*{20}{l}} {\rho \frac{{\partial {p^*}}}{{\partial x}} - \rho \frac{{\partial {u^*}}}{{\partial t}} = 0}\\ {\rho \frac{{\partial {p^*}}}{{\partial z}} - \rho \frac{{\partial {w^*}}}{{\partial t}} = 0}\\ {\frac{{\partial {u^*}}}{{\partial x}} + \frac{{\partial {u^*}}}{{\partial z}} - \frac{1}{{{v^2}}}\frac{{\partial {p^*}}}{{\partial t}} = ({\mathit{\boldsymbol{d}}_{{\rm{cal}}}} - \sum\limits_{i = 1}^N {{E_i}} {\mathit{\boldsymbol{d}}_i})} \end{array}} \right. $

(15)

式中:u^*, w^*, p^*分别为伴随质点的水平分量和垂直分量速度波场以及伴随声压波场。最小二乘逆时偏移最速下降法采用的梯度方程为:

$ \mathit{\boldsymbol{g}} = \int {{p_0}} {p^*}{\rm{d}}t $

(16)

式中: g为梯度。

3.2 静态编码压制串扰噪声

可控震源产生的可控扫描信号为混叠地震数据编码提供了保障, 炸药震源因子波信号不可控, 在激发过程中无法实现编码, 因此在混叠地震数据直接成像时, 可控震源占有绝对的优势。但可控震源特有的采集方式常常会导致混叠噪声的产生, 因此可控震源采用编码方式进行激发, 可以使相邻的地震道数据相干性降到最低, 从而在最小二乘逆时偏移成像过程中减少混叠噪声。常用的可控震源编码方式如下。

1) 伪随机扫描编码^[27]:震源激发时采用了一种由恒频载波信号与伪随机序列信号^[28]相乘得到的伪随机非线性扫描信号, 该信号作为一种较为新型的扫描信号因振动幅度变化大, 故长期使用会对震源车造成损害。因此本文对其不做详细描述。

2) 可控震源混合编码:其原理是将极性编码、随机时间延迟编码、频率编码和相位编码结合起来进行混合编码, 应用于不同震源车激发的扫描信号中, 使得多震源炮记录之间的相干性尽可能地低。该编码方式可提高混叠地震数据分离与直接成像的效率。将混合编码应用于可控震源混叠地震数据采集, 对混叠地震数据分别采用各自震源的扫描信号进行处理, 将其它震源车的炮记录视为噪声, 这方便了混叠地震数据噪声的压制。如果对可控震源混合编码的炮记录进行最小二乘逆时偏移成像处理, 即可很好地消除混叠噪声干扰, 提高成像剖面的信噪比。可控震源混合编码的扫描信号s_i为:

$ \begin{array}{*{20}{l}} {{s_i}(t) = }\\ {\left\{ {\begin{array}{*{20}{c}} 0&{t < E_i^{\rm{d}}}\\ {E_i^{\rm{p}}{a_1}(t - E_i^{\rm{d}}){\rm{sin}}[2\pi E_i^\varphi (t) + E_i^\psi ]}&{E_i^{\rm{d}} \le t \le T} \end{array}} \right.} \end{array} $

(17)

式中:i表示震源车的编号; a₁为振幅函数; E_i^p、E_i^d、E_i^φ和E_i^ψ分别为第i台震源车的极性编码函数、随机时间延迟编码函数、频率编码函数和相位编码函数。相位编码函数E_i^ψ的设计原则已经在1.2节中给出, E_i^p、E_i^d和E_i^φ的表达式为:

$ {E_i^{\rm{p}} = {\rm{ rand }}( \pm 1)} $

(18)

$ {E_i^{\rm{d}} = U( - {t_{{\rm{max}}}},{t_{{\rm{max}}}})} $

(19)

$ {E_i^\varphi (t) = \left( {E_{i,1}^f + \frac{{E_{i,2}^f - E_{i,1}^f}}{{2T}}t} \right)t} $

(20)

式中:rand(±1)为随机生成的为1或-1的函数; U为均等概率分布函数; t_max为最大计算时间; E_{i, 1}^f和E_{i, 2}^f分别为起始频率和终止频率的编码函数。

3.3 利用整形正则化编码压制串扰噪声

FOMEL^[29]提出将整形正则化编码引入共轭梯度算法并进行了试验, 试验结果表明相较于常规共轭梯度算法, 整形正则化编码得到的结果更好, 且迭代收敛速度更快^[30]。本文将整形正则化编码应用于可控震源混叠地震数据直接成像, 并且引入静态编码技术, 得到了成像效果良好的剖面。

整形正则化采用的目标函数J (m)为:

$ \mathit{\boldsymbol{J}}(\mathit{\boldsymbol{m}}) = \frac{1}{2}\left\| {\mathit{\boldsymbol{Lm}} - {\mathit{\boldsymbol{d}}_{{\rm{obs}}}}} \right\|_2^2 + {\varepsilon ^2}\left\| {\mathit{\boldsymbol{Dm}}} \right\|_2^2 $

(21)

式中: D是正则化算子; ε是正则化参数。上述目标函数的最小二乘解为:

$ \mathit{\boldsymbol{\hat m}} = {({\mathit{\boldsymbol{L}}^{\rm{T}}}\mathit{\boldsymbol{L}} + {\varepsilon ^2}{\mathit{\boldsymbol{D}}^{\rm{T}}}\mathit{\boldsymbol{D}})^{ - 1}}{\mathit{\boldsymbol{L}}^{\rm{T}}}\mathit{\boldsymbol{d}} $

(22)

FOMEL^[29]定义整形正则化编码后, 上述目标函数的最小二乘解为:

$ \mathit{\boldsymbol{\hat m}} = {\lambda ^2}\mathit{\boldsymbol{I}} + \mathit{\boldsymbol{S}}{({\mathit{\boldsymbol{L}}^{\rm{T}}}\mathit{\boldsymbol{L}} - {\lambda ^2}\mathit{\boldsymbol{I}})^{ - 1}}\mathit{\boldsymbol{S}}{\mathit{\boldsymbol{L}}^{\rm{T}}}\mathit{\boldsymbol{d}} $

(23)

式中: S为模型的整形算子; λ为尺度参数; I是单位矩阵。正则化编码的目的是对估计的模型施加约束:传统正则化编码是对模型沿着水平方向或者垂直方向进行平滑; 而整形正则化编码则是由平面波解构滤波器求得沿着构造倾角方向进行滤波构造倾角方法^[30-31]。将整形正则化编码引入混叠地震数据最小二乘逆时偏移, 可以很好地压制混叠地震数据间串扰成像噪声。(22)式和(23)式之间的关系可以表示为:

$ {\varepsilon ^2}{\mathit{\boldsymbol{D}}^{\rm{T}}}\mathit{\boldsymbol{D}} = {\mathit{\boldsymbol{S}}^{ - 1}} - \mathit{\boldsymbol{I}} $

(24)

可控震源混叠地震数据直接成像方法适用于所有可控震源采集得到的混叠地震数据, 当相邻道混叠地震数据之间的相干性越低时, 混叠地震数据直接成像导致的串扰噪声越弱。因此在可控震源采集时, 需尽可能采用多种编码方法进行多炮激发。

3.4 模型试算

利用Marmousi模型(图 7)测试最小二乘逆时偏移(least-squares reverse time migration, LSRTM)方法的效果, 采集参数如下:8组震源车同时激发, 每组震源车以滑动扫描方式激发10次, 每次激发移动32 m, 8组震源车的起始激发位置分别为328, 656, 984, 1 312, 1 640, 1 968, 2 296, 2 624 m。8组震源车采用相同的线性升频扫描信号, 频率范围为5~80 Hz, 扫描时间为10 s, 每次激发的记录时间为2 s, 每组震源车两次激发的滑动时间为2 s。图 8为可控震源混叠地震数据第1次激发得到的正演模拟记录, 存在明显的混叠噪声; 对其进行逆时偏移(RTM)得到的成像剖面如图 9a所示, 可以看出图中存在严重的串扰成像噪声, 且地震资料信噪比低; 对其进行LSRTM 50次迭代后的成像结果如图 9b所示, 可以看出多次迭代后, 大部分串扰成像噪声得到了很好的压制。

对8台震源车的扫描信号生成的可控震源混叠地震数据, 采用基于静态极性编码的可控震源LSRTM方法成像, 50次迭代后的成像剖面如图 10a所示, 可以看出串扰成像噪声被进一步压制, 信噪比得到了提高。对采用混合编码函数生成的可控震源混叠地震数据应用基于整形正则化编码的可控震源LSRTM方法(本文方法)进行成像, 而后利用平面波解构滤波器得到构造倾角方向, 再将该构造倾角方向作为约束条件得到第2次迭代的剖面结果, 以此类推得到的50次迭代后的成像剖面如图 10b所示, 可以看出整形正则化编码使得串扰成像噪声被更好地压制。采用本文方法可以很好地压制串扰成像噪声, 提高成像质量, 同时简化了处理流程, 避免了混叠地震数据分离与压制过程中的信号失真及有效能量损失。

3.5 实际资料试算

可控震源DS4的采集方式决定了其混叠地震数据满足随机时间延迟编码的要求, 因此可以采用本文方法对可控震源DS4激发的实际生产数据进行成像。图 11a为常规RTM的成像剖面, 可以看出存在明显的串扰成像噪声, 同相轴不连续, 特别在中深层区域, 成像质量差。采用基于整形正则化编码的可控震源LSRTM方法50次迭代后成像剖面如图 11b所示, 对比可知, 后者的成像质量明显更佳, 同相轴清晰且连续, 能量均衡, 且混叠地震数据采集造成的串扰成像噪声被很好地压制。综上可知, 本文提出的可控震源混叠地震数据直接成像方法具有良好的适用性和有效性。

4 讨论

随着可控震源采集效率的提高, 混叠噪声成为影响成像剖面质量不可忽略的因素, 通常采用以下方法思路对其进行压制:

1) 利用反演的思想, 对混叠地震数据进行分离。对于可控震源HFVS和ISS扫描等采集方法得到的地震数据可采用该思路进行分离。首先对可控震源HFVS采集得到的混叠地震数据, 利用每次扫描不改变地层响应函数的原理进行计算, 从试算结果来看, 近道和远道混叠地震数据的分离效果均良好; 然后通过近源力信号来设计系数矩阵, 在此过程中, 因考虑了谐波畸变的影响, 故分离后的炮记录中不存在谐波干扰, 混叠地震数据分离效果取决于相位编码的优劣; 最后对可控震源ISS扫描得到的混叠地震数据, 基于随机时间延迟编码, 若对不同震源车的扫描信号进行混合编码, 如极性编码、相位编码、频率编码等, 均可以减少不同震源车炮记录之间的相关性, 改善分离效果。利用上述两种方法得到的混叠地震数据分离均需引入广义逆算子、奇异值分解等手段来改善分离效果。

2) 利用去噪的思想, 对混叠噪声进行压制。利用DS4和ISS扫描等采集方法得到的地震数据可采用这种思路进行混叠噪声压制。混叠噪声在炮域中与有效信号的特征相同, 变换到人工分选道集如共中心点域则呈现随机噪声形态, 采用常规随机噪声去除方法如矢量中值滤波可以较好地压制混叠噪声。若在采集时引入编码函数, 可提高混叠噪声压制效率, 并且改善混叠噪声压制效果。

3) 目前比较热门的方法是利用基于整形正则化编码的可控震源LSRTM方法直接成像。该方法通过简化处理流程, 提高处理速度, 但在成像过程中易引入串扰成像噪声, 对此需要进行特殊的处理, 如对成像剖面进行滤波等。滤波对串扰成像噪声压制效果不佳, 而且保幅效果差。因此本文首先通过最小二乘逆时偏移方法压制串扰成像噪声, 然后引入静态编码技术进一步压制串扰成像噪声, 最后利用整形规则化滤波技术消除串扰成像噪声, 同时压制低频噪声、改善振幅均衡性、提高深部能量、提高成像剖面的分辨率。

在实际资料处理中, 如果相邻混叠地震数据的相干性较强, 单独使用混叠地震数据分离方法和混叠地震数据成像方法均无法很好地对混叠噪声进行压制, 因此需要同时使用两种方法, 以便更好地压制串扰成像噪声、提高成像质量。

5 结论

可控震源高效采集技术提高了采集效率, 但却引入了大量的混叠噪声, 造成地震成像信噪比低。

针对可控震源混叠采集引入的混叠噪声, 本文采用3种方法进行处理。对于HFVS数据, 采用基于广义逆算子和奇异值分解的可控震源混叠地震数据分离方法, 基于反演思想进行混叠地震数据分离, 分离效果取决于相位编码矩阵的选取, 通过设计条件数低的相位编码矩阵提高高保真数据的分离效果; 对于ISS数据, 采用基于地表一致性原理的可控震源混叠噪声压制方法, 基于去噪思想进行混叠噪声压制, 该方法将混叠地震数据变换到人工分选道集, 采用矢量中值滤波随机噪声压制技术压制混叠噪声。以上两种方法对数据采集的要求较高, 且在混叠地震数据分离与混叠噪声压制过程中会造成有效能量的损失。因此, 本文采用基于整形正则化编码的可控震源LSRTM方法对可控震源混叠地震数据进行成像, 既简化了处理流程, 压制混叠噪声后又得到了高信噪比、高分辨率、高保真度的成像结果。

模型数据和实际数据测试表明, 本文中3种针对混叠地震数据的处理方法和思路, 特别是基于整形正则化编码的可控震源LSRTM具有广阔的应用前景。但基于整形正则化编码的可控震源LSRTM方法具有明显的局限性:该方法利用最小二乘反演思想通过不断迭代改善成像质量, 计算量大, 在实际应用时需要考虑计算成本。