0 引言

地震数据采集需要考虑成本与质量之间的平衡。传统的地震勘探野外数据采集过程中，相邻炮的激发时间间隔较大，采集效率低，成本高。同时震源(simultaneous source)采集技术，也称多震源地震数据采集技术，改变了传统的采集方式。当震源数量一定时，同时震源采集能够缩短采集时间，增大覆盖范围；通过增加震源数量，提高覆盖次数，改善对地下构造的照明质量^[1-3]。

同时震源数据采集能够提高采集效率和改善照明质量，但是相邻震源之间会产生“串扰噪声”，因此同时震源数据不适用于传统的处理流程。目前此类数据的处理方法主要有两大类：①直接成像，即直接用同时震源数据进行偏移成像，这类方法受到串扰噪声的干扰，成像质量会降低^[4-10]；②分离法，即先对混合数据进行分离，得到类似单个震源激发的数据，再用常规地震方法进行处理。

分离法进一步划分为滤波类和反演类。从共炮域变换到其他域，在伪分离记录中主震源信号是连续的，混合噪声不再具有连续相干性。据此，Moore等^[11]和Akerberg等^[12]在共偏移距域进行滤波；Huo等^[13]利用多方向矢量中值滤波法对地层倾角进行扫描，求得最佳倾角后，再利用中值滤波去除共中心点域的混合噪声。Yang等^[14]提出先利用PWD(plane-wave destruction)求取地层的倾角信息，再结合矢量中值滤波对串扰噪声进行压制。反演类方法的思路是施加约束条件，将分离问题变为反演问题。Doulgeris等^[15]设计了迭代算子，在共检波点域压制噪声；Mahdad等^[16]在频率—波数域分选伪分离数据，设计迭代算法，分离混合波场；韩立国等^[17]结合多级中值滤波和Curvelet阈值迭代去噪算法，实现了混合数据的分离；Chen等^[18]、祖绍环等^[19]、Xue等^[20]和宋家文等^[21]分别将地震数据转换到Seislet域、Curvelet域、Radon域和频率—波数—波数域(FKK)，然后在稀疏域对混合地震数据进行正则化稀疏约束迭代，实现了混合数据的分离；Zu等^[22]引入两个卷积算子构造逆问题，先求得地震信号的倾角，然后通过共轭梯度算法求解该逆问题，实现混叠数据的分离；Zhou等^[23]在共检波点域利用秩缩减滤波器和相干滤波器，通过迭代压制非相干信号，也取得了不错的效果。

以上分离方法都要施加约束条件，比如地震信号的相干性和稀疏性。Wapenaar等^[24-26]在地震多维反褶积干涉的基础上，对同时震源数据的分离进行了研究，认为可以不施加这些额外的约束条件。在时间延迟算子的伪逆中，增加了基础点扩散矩阵，在满足一定的分离条件下，将混合数据按延迟时间进行归位，直接得到分离结果。

本文通过推广基础点扩散矩阵，将检波点和炮点规则排列采集的同时震源数据分离方法拓展应用于不规则排列采集的同时震源数据分离。

1 方法原理 1.1 数据混合的编码方式

同时震源数据采集中，原始未混合数据可表示为p(x_R^(k)，x_S⁽ⁱ⁾, t)，这里x_S⁽ⁱ⁾为第i个震源的位置，x_R^(k)为第k个检波器的位置，t代表时间。相邻单炮通过随机延迟时间编码的方式组合为一个超级炮，记为σ^(m)，m为超级炮编号。若该超级炮内的每一炮在相近的延迟时间t_i时刻激发，可以得到混合炮数据

$ {P_{{\rm{sim}}}}\left( {\mathit{\boldsymbol{x}}_{\rm{R}}^{\left( k \right)},{\mathit{\boldsymbol{\sigma }}^{\left( m \right)}},t} \right) = \sum\limits_{x_{\rm{S}}^{\left( i \right)} \in {\mathit{\boldsymbol{\sigma }}^{\left( m \right)}}} p \left( {\mathit{\boldsymbol{x}}_{\rm{R}}^{\left( k \right)},\mathit{\boldsymbol{x}}_{\rm{S}}^{\left( i \right)},t - {t_i}} \right) $

(1)

式中x_S⁽ⁱ⁾∈σ^(m)，表示将超级炮σ^(m)中的全部单炮x_S⁽ⁱ⁾加起来。通过傅里叶变换，式(1)变为

$ \begin{array}{l} {P_{{\rm{sim}}}}\left( {\mathit{\boldsymbol{x}}_{\rm{R}}^{\left( k \right)},{\mathit{\boldsymbol{\sigma }}^{\left( m \right)}},\omega } \right) = \\ \;\;\;\;\;\;\; = \sum\limits_{x_S^{\left( i \right)} \in {\mathit{\boldsymbol{\sigma }}^{\left( m \right)}}} p \left( {\mathit{\boldsymbol{x}}_{\rm{R}}^{\left( k \right)},\mathit{\boldsymbol{\sigma }}_{\rm{S}}^{\left( i \right)},\omega } \right)\exp \left( { - {\rm{i}}\omega {t_i}} \right) \end{array} $

(2)

式中ω为角频率。对每一个角频率分量，在频率域的混合炮数据可用矩阵的形式(图 1)表示^[1]为

$ {\mathit{\boldsymbol{P}}_{{\rm{sim}}}} = \mathit{\boldsymbol{P \boldsymbol{\varGamma} }} $

(3)

式中：矩阵P中的第k行、第i列表示原始数据p(x_R^(k), x_S⁽ⁱ⁾, ω)，如果有K个检波器和N个震源，P就是一个K×N的矩阵；Γ是时间延迟矩阵，如果是相邻n炮同时激发，Γ为N×(N/n)阶矩阵。n=2时的Γ如图 1所示，其中b_i=e^-iωt_i。此时混合后的矩阵P_sim的元素表示为p_sim(x_R^(k), σ^(m), ω)。可以看出，混合后共炮点道集数少于单炮下的共炮点道集数，因此式(3)是欠定的。

1.2 混合数据的伪分离

仅从方程的角度看，同时震源数据的分离就是通过式(3)求解原始数据矩阵P，所以形式上可以得到分离结果

$ \mathit{\boldsymbol{\hat P}} = {\mathit{\boldsymbol{P}}_{{\rm{sim}}}}{{\mathit{\boldsymbol{ \boldsymbol{\hat \varGamma} }}}_{{\rm{inv}}}} $

(4)

式中$\mathit{\boldsymbol{ \boldsymbol{\hat \varGamma} }}$_inv为分离算子。如果用欠定问题下的最小二乘法求解P，则$\mathit{\boldsymbol{ \boldsymbol{\hat \varGamma} }}$_inv表示为时间延迟矩阵Γ的伪逆

$ {{\mathit{\boldsymbol{ \boldsymbol{\hat \varGamma} }}}_{{\rm{inv}}}} = {\left( {{\mathit{\boldsymbol{ \boldsymbol{\varGamma} }}^\dagger }\mathit{\boldsymbol{ \boldsymbol{\varGamma} }}} \right)^{ - 1}}{\mathit{\boldsymbol{ \boldsymbol{\varGamma} }}^\dagger } $

(5)

式中：“†”代表复共轭转置；Γ^†Γ=nI，I为单位矩阵，n代表混合度；$\mathit{\boldsymbol{ \boldsymbol{\hat \varGamma} }}$_inv=1/nΓ^†。因此，通过最小二乘法求解可得到伪分离记录^[16-17]

$ \mathit{\boldsymbol{\hat P}} = \frac{1}{n}{\mathit{\boldsymbol{P}}_{{\rm{sim}}}}{\mathit{\boldsymbol{ \boldsymbol{\varGamma} }}^\dagger } $

(6)

式(6)表明伪分离法是先复制n份混合数据，再按主炮的延迟时间归位^[17]，但仍然保留其他炮的信号，形成串扰噪声。

1.3 构建基础点扩散矩阵

在伪分离中，干扰炮按照主炮的延迟时间处理，而非按照自己的延迟时间归位，形成所谓的“串扰噪声”。不同于常规加约束的迭代反演方法，这里不是简单地将这些干扰炮的“串扰”看成噪声，而是认为这些“串扰噪声”没有得到正确归位。通过讨论同时震源数据的分离与多维反褶积地震干涉技术在原理上的相似性，引入基础点扩散矩阵对混合数据进行处理^[25-26]，将“串扰噪声”按各自的延迟时间进行归位。

地震干涉技术可以将任意两个检波器接收到的数据合成为在若干检波器之间传播的波，相当于将其中的一个检波器看作虚拟震源(图 2)。

首先，对每次单个激发的震源构建如图 2a所示的地震卷积模型。炮点位于地表，检波点x位于界面T上，选取一个参考检波点记为x_B，F为目标层。对于单个震源x_S⁽ⁱ⁾在t_i时刻激发，地震波直接传播到检波器x，记录的波场表示为uⁱⁿ(x, x_S⁽ⁱ⁾, t)；地震波经过目标层反射后传播到检波器x_B，记录的波场表示为u^out(x_B, x_S⁽ⁱ⁾, t)；波场与相应的格林函数的关系可以表示为

$ {u^{{\rm{in}}}}\left( {\mathit{\boldsymbol{x}},\mathit{\boldsymbol{x}}_{\rm{S}}^{\left( i \right)},t} \right) = {G^{{\rm{in}}}}\left( {\mathit{\boldsymbol{x}},\mathit{\boldsymbol{x}}_{\rm{S}}^{\left( i \right)},t} \right) * {s^{\left( i \right)}}\left( {t - {t_i}} \right) $

(7)

$ {u^{{\rm{out}}}}\left( {{\mathit{\boldsymbol{x}}_{\rm{B}}},\mathit{\boldsymbol{x}}_{\rm{S}}^{\left( i \right)},t} \right) = {G^{{\rm{out}}}}\left( {{\mathit{\boldsymbol{x}}_{\rm{B}}},\mathit{\boldsymbol{x}}_{\rm{S}}^{\left( i \right)},t} \right) * {s^{\left( i \right)}}\left( {t - {t_i}} \right) $

(8)

式中：Gⁱⁿ(x, x_S⁽ⁱ⁾, t)和G^out(x_B, x_S⁽ⁱ⁾, t)分别为震源x_S⁽ⁱ⁾到x与x_B的格林函数；s⁽ⁱ⁾(t-t_i)表示延迟激发的震源子波。在此模型中，多维反褶积地震干涉法的直接求解表达式为^[24]

$ {G^{{\rm{out}}}}\left( {{\mathit{\boldsymbol{x}}_{\rm{B}}},\mathit{\boldsymbol{x}}_{\rm{S}}^{\left( i \right)},t} \right) = \int_{\rm{T}} {\bar G_{\rm{d}}^{{\rm{out}}}} \left( {{\mathit{\boldsymbol{x}}_{\rm{B}}},\mathit{\boldsymbol{x}},t} \right) * {G^{{\rm{in}}}}\left( {\mathit{\boldsymbol{x}},\mathit{\boldsymbol{x}}_{\rm{S}}^{\left( i \right)},t} \right){\rm{d}}\mathit{\boldsymbol{x}} $

(9)

式中$\bar G_{\rm{d}}^{{\rm{out }}}\left( {{{\boldsymbol{x}}_{\rm{B}}}, {\boldsymbol{x}}, t} \right)$为x_B与充当虚拟源的检波器x之间的待求格林函数。根据式(7)和式(8)，将式(9)两端同时与s⁽ⁱ⁾(t-t_i)进行褶积，得到

$ {u^{{\rm{out}}}}\left( {{\mathit{\boldsymbol{x}}_{\rm{B}}},\mathit{\boldsymbol{x}}_{\rm{S}}^{\left( i \right)},t} \right) = \int_{\rm{T}} {\bar G_{\rm{d}}^{{\rm{out}}}} \left( {{\mathit{\boldsymbol{x}}_{\rm{B}}},\mathit{\boldsymbol{x}},t} \right) * {u^{{\rm{in}}}}\left( {\mathit{\boldsymbol{x}},\mathit{\boldsymbol{x}}_{\rm{S}}^{\left( i \right)},t} \right){\rm{d}}\mathit{\boldsymbol{x}} $

(10)

在界面T上选取一个固定的参考检波器x_A，其记录的波场为uⁱⁿ(x_A, x_S⁽ⁱ⁾, t)，若在式(10)两端与uⁱⁿ(x_A, x_S⁽ⁱ⁾, t)做互相关计算，然后对全部炮求和，得到

$ {C_{{\rm{seq}}}}\left( {{\mathit{\boldsymbol{x}}_{\rm{B}}},{\mathit{\boldsymbol{x}}_{\rm{A}}},t} \right) = \int_{\rm{T}} {\bar G_{\rm{d}}^{{\rm{out}}}} \left( {{\mathit{\boldsymbol{x}}_{\rm{B}}},\mathit{\boldsymbol{x}},t} \right) * {R_{{\rm{seq}}}}\left( {\mathit{\boldsymbol{x}},{\mathit{\boldsymbol{x}}_{\rm{A}}},t} \right){\rm{d}}\mathit{\boldsymbol{x}} $

(11)

式中：C_seq(x_B, x_A, t)为相关函数；R_seq(x, x_A, t)为点扩散函数。其表达式分别为

$ {C_{{\rm{seq}}}}\left( {{\mathit{\boldsymbol{x}}_{\rm{B}}},{\mathit{\boldsymbol{x}}_{\rm{A}}},t} \right) = \sum\limits_i {{u^{{\rm{out}}}}\left( {{\mathit{\boldsymbol{x}}_{\rm{B}}},\mathit{\boldsymbol{x}}_{\rm{S}}^{\left( i \right)},t} \right) * {u^{{\rm{in}}}}\left( {{\mathit{\boldsymbol{x}}_{\rm{A}}},\mathit{\boldsymbol{x}}_{\rm{S}}^{\left( i \right)}, - t} \right)} $

(12)

$ {R_{{\rm{seq}}}}\left( {\mathit{\boldsymbol{x}},{\mathit{\boldsymbol{x}}_{\rm{A}}},t} \right) = \sum\limits_i {{u^{{\rm{in}}}}\left( {\mathit{\boldsymbol{x}},\mathit{\boldsymbol{x}}_{\rm{S}}^{\left( i \right)},t} \right) * {u^{{\rm{in}}}}\left( {{\mathit{\boldsymbol{x}}_{\rm{A}}},\mathit{\boldsymbol{x}}_{\rm{S}}^{\left( i \right)}, - t} \right)} $

(13)

若求得C_seq(x_B, x_A, t)与R_seq(x, x_A, t)并代入式(11)，通过傅里叶变换，式(11)变为乘积表达式，使用最小二乘法求解并变换到时域，可以求得$\bar G_{\rm{d}}^{{\rm{out }}}\left( {{{\boldsymbol{x}}_{\rm{B}}}, {\boldsymbol{x}}, t} \right)$。

上面只对单个点源x_S⁽ⁱ⁾进行了讨论，对于超级炮σ^(m)(图 2b)，若其中的全部单炮在相隔很短的时间内激发，则常规的地震干涉模型变为同时震源模型^[26]。与式(7)和式(8)相似，这些波场与相应的格林函数的关系可以表示为

$ {u^{{\rm{in}}}}\left( {\mathit{\boldsymbol{x}},{\mathit{\boldsymbol{\sigma }}^{\left( m \right)}},t} \right) = \sum\limits_{\mathit{\boldsymbol{x}}_{\rm{S}}^{\left( i \right)} \in {\mathit{\boldsymbol{\sigma }}^{\left( m \right)}}} {{G^{{\rm{in}}}}\left( {\mathit{\boldsymbol{x}},\mathit{\boldsymbol{x}}_{\rm{S}}^{(i)},t} \right) * {s^{\left( i \right)}}\left( {t - {t_i}} \right)} $

(14)

$ {u^{{\rm{out}}}}\left( {{\mathit{\boldsymbol{x}}_{\rm{B}}},{\mathit{\boldsymbol{\sigma }}^{\left( m \right)}},t} \right) = \sum\limits_{\mathit{\boldsymbol{x}}_{\rm{S}}^{\left( i \right)} \in {\mathit{\boldsymbol{\sigma }}^{\left( m \right)}}} {{G^{{\rm{out}}}}\left( {{\mathit{\boldsymbol{x}}_{\rm{B}}},\mathit{\boldsymbol{x}}_{\rm{S}}^{(i)},t} \right) * {s^{\left( i \right)}}\left( {t - {t_i}} \right)} $

(15)

在式(9)两边同时与s⁽ⁱ⁾(t-t_i)褶积，并对超级炮σ^(m)内的所有单炮求和，得到类似于式(10)的表达式

$ {u^{{\rm{out}}}}\left( {{\mathit{\boldsymbol{x}}_{\rm{B}}},{\mathit{\boldsymbol{\sigma }}^{\left( m \right)}},t} \right) = \int_{\rm{T}} {\bar G_{\rm{d}}^{{\rm{out}}}} \left( {{\mathit{\boldsymbol{x}}_{\rm{B}}},\mathit{\boldsymbol{x}},t} \right) * {u^{{\rm{in}}}}\left( {\mathit{\boldsymbol{x}},{\mathit{\boldsymbol{\sigma }}^{\left( m \right)}},t} \right){\rm{d}}\mathit{\boldsymbol{x}} $

(16)

对界面T上选取的一个固定参考检波器x_A，其记录的波场为uⁱⁿ(x_A, σ^(m), t)，在式(16)两端同时与uⁱⁿ(x_A, σ^(m), t)做互相关计算，并对所有超级炮σ^(m)求和，则式(11)可以改写为同时震源模型下的相关函数

$ {C_{{\rm{sim}}}}\left( {{\mathit{\boldsymbol{x}}_{\rm{B}}},{\mathit{\boldsymbol{x}}_{\rm{A}}},t} \right) = \int_{\rm{T}} {\bar G_{\rm{d}}^{{\rm{out}}}} \left( {{\mathit{\boldsymbol{x}}_{\rm{B}}},\mathit{\boldsymbol{x}},t} \right) * {R_{{\rm{sim}}}}\left( {\mathit{\boldsymbol{x}},{\mathit{\boldsymbol{x}}_{\rm{A}}},t} \right){\rm{d}}\mathit{\boldsymbol{x}} $

(17)

式中：C_sim(x_B, x_A, t)代表伪分离后的数据；R_sim(x, x_A, t)为点扩散函数。其表达式分别为

$ \begin{array}{l} {C_{{\rm{sim}}}}\left( {{\mathit{\boldsymbol{x}}_{\rm{B}}},{\mathit{\boldsymbol{x}}_{\rm{A}}},t} \right) = \sum\limits_m {\left[ {{u^{{\rm{out}}}}\left( {{\mathit{\boldsymbol{x}}_{\rm{B}}},{\mathit{\boldsymbol{\sigma }}^{\left( m \right)}},t} \right) * } \right.} \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\left. {{u^{{\rm{in}}}}\left( {{\mathit{\boldsymbol{x}}_{\rm{A}}},{\mathit{\boldsymbol{\sigma }}^{\left( m \right)}}, - t} \right)} \right] \end{array} $

(18)

$ \begin{array}{l} {R_{{\rm{sim}}}}\left( {\mathit{\boldsymbol{x}},{\mathit{\boldsymbol{x}}_{\rm{A}}},t} \right) = \sum\limits_m {\left[ {{u^{{\rm{in}}}}\left( {\mathit{\boldsymbol{x}},{\mathit{\boldsymbol{\sigma }}^{\left( m \right)}},t} \right) * } \right.} \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\left. {{u^{{\rm{in}}}}\left( {{\mathit{\boldsymbol{x}}_{\rm{A}}},{\mathit{\boldsymbol{\sigma }}^{\left( m \right)}}, - t} \right)} \right] \end{array} $

(19)

通过式(18)和式(19)分别求得C_sim(x_B, x_A, t)和R_sim(x, x_A, t)。在频率域式(17)是一个乘积表达式，通过最小二乘法求解并变换到时域，可以求出解混之后的响应$\bar G_{\rm{d}}^{{\rm{out }}}\left( {{{\boldsymbol{x}}_{\rm{B}}}, {\boldsymbol{x}}, t} \right)$。图 2b中，检波器排列在地下界面T上，若将其外推到地表，可以得到检波器位于地表条件下的同时震源解混之后的响应，此时求得的$\bar G_{\rm{d}}^{{\rm{out }}}\left( {{{\boldsymbol{x}}_{\rm{B}}}, {\boldsymbol{x}}, t} \right)$即为本文方法希望求取的解混后的数据$\hat P\left( {{{\boldsymbol{x}}_{\rm{B}}}, {\boldsymbol{x}}, t} \right)$。下面用矩阵的形式表示上述同时震源混合数据分离过程。

将式(9)变换到频率域，每个频率分量都有如下离散化后的矩阵表达式

$ {\mathit{\boldsymbol{G}}^{{\rm{out}}}} = \mathit{\boldsymbol{\bar G}}_{\rm{d}}^{{\rm{out}}}{\mathit{\boldsymbol{G}}^{{\rm{in}}}} $

(20)

式(14)和式(15)表示同时震源模型中的波场响应，其离散矩阵为

$ {\mathit{\boldsymbol{U}}^{{\rm{in}}}} = {\mathit{\boldsymbol{G}}^{{\rm{in}}}}\mathit{\boldsymbol{ \boldsymbol{\varGamma} }} $

(21)

$ {\mathit{\boldsymbol{U}}^{{\rm{out}}}} = {\mathit{\boldsymbol{G}}^{{\rm{out}}}}\mathit{\boldsymbol{ \boldsymbol{\varGamma} }} $

(22)

在式(20)两端同乘以Γ，得到式(16)的离散矩阵表达式

$ {\mathit{\boldsymbol{U}}^{{\rm{out}}}} = \mathit{\boldsymbol{\bar G}}_{\rm{d}}^{{\rm{out}}}{\mathit{\boldsymbol{U}}^{{\rm{in}}}} $

(23)

由于同时震源的分离问题是一个欠定的问题，在不加约束的情况下，利用欠定条件下的最小二乘法对式(23)进行求解，可以求得$\mathit{\boldsymbol{\bar G}}_{\rm{d}}^{{\rm{out }}}$

$ \mathit{\boldsymbol{\bar G}}_{\rm{d}}^{{\rm{out}}} = {\mathit{\boldsymbol{U}}^{{\rm{out}}}}{\left[ {{{\left( {{\mathit{\boldsymbol{U}}^{{\rm{in}}}}} \right)}^\dagger }{\mathit{\boldsymbol{U}}^{{\rm{in}}}}} \right]^{ - 1}}{\left( {{\mathit{\boldsymbol{U}}^{{\rm{in}}}}} \right)^\dagger } $

(24)

根据式(19)，R_sim在频率域的矩阵表达式为

$ {\mathit{\boldsymbol{R}}_{{\rm{sim}}}} = {\left( {{\mathit{\boldsymbol{U}}^{{\rm{in}}}}} \right)^\dagger }{\mathit{\boldsymbol{U}}^{{\rm{in}}}} $

(25)

将式(21)代入式(25)，得到

$ {\mathit{\boldsymbol{R}}_{{\rm{sim}}}} = {\mathit{\boldsymbol{ \boldsymbol{\varGamma} }}^\dagger }{\left( {{\mathit{\boldsymbol{G}}^{{\rm{in}}}}} \right)^\dagger }{\mathit{\boldsymbol{G}}^{{\rm{in}}}}\mathit{\boldsymbol{ \boldsymbol{\varGamma} }} = {\mathit{\boldsymbol{ \boldsymbol{\varGamma} }}^\dagger }\mathit{\boldsymbol{R \boldsymbol{\varGamma} }} $

(26)

式中(Gⁱⁿ)^†Gⁱⁿ被定义为基础点扩散矩阵(the basic point-spread function)，用符号R表示

$ \mathit{\boldsymbol{R}} = {\left( {{\mathit{\boldsymbol{G}}^{{\rm{in}}}}} \right)^\dagger }{\mathit{\boldsymbol{G}}^{{\rm{in}}}} $

(27)

在频率—波数域，远场响应的空间带宽是有限制的^[25]，范围通常在$-\frac{|\omega|}{v_{\mathrm{a}}} \sim \frac{|\omega|}{v_{\mathrm{a}}}$，其中$v_{\mathrm{a}}=\frac{v_{\mathrm{s}}}{\sin \alpha_{\max }}$，v_s表示介质中的最小速度，α_max表示最大入射角。频率—空间域的基础点扩散函数^[25]定义如下

$ \gamma \left( {{x_i},\omega } \right) = \frac{{\sin \frac{{\left| \omega \right|{x_i}}}{{{v_{\rm{a}}}}}}}{{{\rm{ \mathsf{ π} }}{x_i}}} $

(28)

式中：x_i为第i个震源到第1个震源的距离；γ(x_i, ω)与二维偏移中的分辨率函数类似^{[24, 27]}，与sinc插值函数成正比^{[25, 28]}。在信号重建中，重建函数^[29]

$ g\left( x \right) = \sum\limits_i {{g_{\rm{d}}}} \left( i \right){\rm sinc}\left( {x - i} \right) $

(29)

表示原始信号g(x)是待重建信号g_d(x)与sinc(x)的卷积，由此将sinc(x)离散化。同理将γ(x_i, ω)离散化并表示为矩阵R

$ \mathit{\boldsymbol{R}} = \left( {\begin{array}{*{20}{c}} {{\gamma _0}}&{{\gamma _1}}&{{\gamma _2}}& \cdots &{}&{}\\ {{\gamma _{ - 1}}}&{{\gamma _0}}&{{\gamma _1}}&{{\gamma _2}}& \cdots &{}\\ {{\gamma _{ - 2}}}&{{\gamma _{ - 1}}}&{{\gamma _0}}&{{\gamma _1}}&{{\gamma _2}}& \cdots \\ {}&{}& \vdots &{}&{}&{} \end{array}} \right) $

(30)

式中γ_i=γ(x_i, ω)=γ(iΔs, ω)，Δs为炮间距。

根据式(24)可以求得解混之后的响应$\mathit{\boldsymbol{\bar G}}_{\rm{d}}^{{\rm{out }}}$，如图 2b所示。在该模型中, 虚拟源x与检波点x_B都位于地下界面T上，引入外推算子W^-和W⁺，将$\mathit{\boldsymbol{\bar G}}_{\rm{d}}^{{\rm{out }}}$外推到地表，得到检波器位于地表条件下的同时震源分离数据$\mathit{\boldsymbol{\hat P}}$

$ \mathit{\boldsymbol{\hat P}} = {\mathit{\boldsymbol{W}}^ - }\mathit{\boldsymbol{\bar G}}_{\rm{d}}^{{\rm{out}}}{\mathit{\boldsymbol{W}}^ + } $

(31)

将式(24)代入上式，并用Gⁱⁿ替代W^+[26]，得到

$ \mathit{\boldsymbol{\hat P}} = {\mathit{\boldsymbol{W}}^ - }{\mathit{\boldsymbol{U}}^{{\rm{out}}}}{\left[ {{{\left( {{\mathit{\boldsymbol{U}}^{{\rm{in}}}}} \right)}^\dagger }{\mathit{\boldsymbol{U}}^{{\rm{in}}}}} \right]^{ - 1}}{\left( {{\mathit{\boldsymbol{U}}^{{\rm{in}}}}} \right)^\dagger }{\mathit{\boldsymbol{G}}^{{\rm{in}}}} $

(32)

其中W^-U^out为混合后的矩阵P_sim，分别将式(21)和式(27)代入上式，得到同时震源的直接反演分离公式

$ \begin{array}{l} \mathit{\boldsymbol{\hat P}} = {\mathit{\boldsymbol{P}}_{ sim}}{\left[ {{{\left( {{\mathit{\boldsymbol{U}}^{{\rm{in}}}}} \right)}^\dagger }{\mathit{\boldsymbol{U}}^{{\rm{in}}}}} \right]^{ - 1}}{\left( {{\mathit{\boldsymbol{U}}^{{\rm{in}}}}} \right)^\dagger }{\mathit{\boldsymbol{G}}^{{\rm{in}}}}\\ \;\;\; = {\mathit{\boldsymbol{P}}_{ sim}}{\left[ {{{\left( {{\mathit{\boldsymbol{G}}^{{\rm{in}}}}\mathit{\boldsymbol{ \boldsymbol{\varGamma} }}} \right)}^\dagger }{\mathit{\boldsymbol{G}}^{{\rm{in}}}}\mathit{\boldsymbol{ \boldsymbol{\varGamma} }}} \right]^{ - 1}}{\left( {{\mathit{\boldsymbol{G}}^{{\rm{in}}}}\mathit{\boldsymbol{ \boldsymbol{\varGamma} }}} \right)^\dagger }{\mathit{\boldsymbol{G}}^{{\rm{in}}}}\\ \;\;\; = {\mathit{\boldsymbol{P}}_{ sim}}{\left[ {{\mathit{\boldsymbol{ \boldsymbol{\varGamma} }}^\dagger }{{\left( {{\mathit{\boldsymbol{G}}^{{\rm{in}}}}} \right)}^\dagger }{\mathit{\boldsymbol{G}}^{{\rm{in}}}}\mathit{\boldsymbol{ \boldsymbol{\varGamma} }}} \right]^{ - 1}}{\mathit{\boldsymbol{ \boldsymbol{\varGamma} }}^\dagger }{\left( {{\mathit{\boldsymbol{G}}^{{\rm{in}}}}} \right)^\dagger }{\mathit{\boldsymbol{G}}^{{\rm{in}}}}\\ \;\;\; = {\mathit{\boldsymbol{P}}_{ sim}}{\left( {{\mathit{\boldsymbol{ \boldsymbol{\varGamma} }}^\dagger }\mathit{\boldsymbol{R \boldsymbol{\varGamma} }}} \right)^{ - 1}}{\mathit{\boldsymbol{ \boldsymbol{\varGamma} }}^\dagger }\mathit{\boldsymbol{R}} \end{array} $

(33)

式中点扩散矩阵为R_sim=Γ^†RΓ，则分离算子为$\mathit{\boldsymbol{ \boldsymbol{\hat \varGamma} }}$_inv=R_sim^-1 Γ^†R。现举列讨论矩阵R_sim在同时震源数据的直接反演分离中的作用，为了简单起见，假设混合度n=2，且每个超级炮的延迟时间相同，则时间延迟矩阵Γ为

$ \mathit{\boldsymbol{ \boldsymbol{\varGamma} }} = \left( {\begin{array}{*{20}{c}} {{b_1}}&0& \cdots &0\\ {{b_2}}&0& \cdots &0\\ 0&{{b_1}}& \cdots &0\\ 0&{{b_2}}& \cdots &0\\ {}&{}& \vdots &{}\\ 0&0& \cdots &{{b_1}}\\ 0&0& \cdots &{{b_2}} \end{array}} \right) $

(34)

式中b₁=exp(-iωt₁), b₂=exp(-iωt₂)。将式(30)和式(34)代入式R_sim=Γ^†RΓ中，得到

$ \begin{array}{l} {\mathit{\boldsymbol{R}}_{{\rm{sim}}}} = 2\left( {\begin{array}{*{20}{c}} {{\gamma _0}}&{{\gamma _2}}&{{\gamma _4}}& \cdots &{}&{}\\ {{\gamma _{ - 2}}}&{{\gamma _0}}&{{\gamma _2}}&{{\gamma _4}}& \cdots &{}\\ {{\gamma _{ - 4}}}&{{\gamma _{ - 2}}}&{{\gamma _0}}&{{\gamma _2}}&{{\gamma _4}}& \cdots \\ {}&{}& \vdots &{}&{}&{} \end{array}} \right) + \\ \;\;\;\;\;\exp \left( { - {\rm{i}}\omega \Delta t} \right)\left( {\begin{array}{*{20}{c}} {{\gamma _1}}&{{\gamma _3}}&{{\gamma _5}}& \cdots &{}&{}\\ {{\gamma _{ - 1}}}&{{\gamma _1}}&{{\gamma _3}}&{{\gamma _5}}& \cdots &{}\\ {{\gamma _{ - 3}}}&{{\gamma _{ - 1}}}&{{\gamma _1}}&{{\gamma _3}}&{{\gamma _5}}& \cdots \\ {}&{}& \vdots &{}&{}&{} \end{array}} \right) + \\ \;\;\;\;\;\exp \left( { + {\rm{i}}\omega \Delta t} \right)\left( {\begin{array}{*{20}{c}} {{\gamma _{ - 1}}}&{{\gamma _1}}&{{\gamma _3}}& \cdots &{}&{}\\ {{\gamma _{ - 3}}}&{{\gamma _{ - 1}}}&{{\gamma _1}}&{{\gamma _3}}& \cdots &{}\\ {{\gamma _{ - 5}}}&{{\gamma _{ - 3}}}&{{\gamma _{ - 1}}}&{{\gamma _1}}&{{\gamma _3}}& \cdots \\ {}&{}& \vdots &{}&{}&{} \end{array}} \right) \end{array} $

(35)

式中Δt=t₂-t₁。式(35)中右边第一项是式(30)中基础点扩散矩阵R的重采样变形式(采样因子为2)，第二和第三项也是基础点扩散矩阵的重采样变形式，但是对距离的改变量为±Δs，在空间上对串扰噪声进行了归位；因子exp $( \mp {\rm{i}}\omega \Delta t)$表示时间改变量±Δt，对串扰噪声在时间上进行了归位。式中的第二和第三项解释了对同时震源数据中的串扰噪声的处理方式。在对上式的分析中，假设一个超级炮是由相邻两个单炮混合组成的。如果一个超级炮由n个相邻震源组成，也能够得到类似于式(35)的表达式，只是基础点扩散矩阵R的重采样变形式在空间和时间的改变量会随之发生变化^[26]。

最后，对求得的频率域数据$\mathit{\boldsymbol{\hat P}}$进行傅里叶逆变换，并取实部，可得时间域的分离数据。

直接反演分离的具体流程如图 3所示。

1.4 分离适用条件

直接反演分离方法利用地震数据带宽有限的特点，使用与sinc插值函数类似的γ(x_i, ω)时，要求直接反演方法满足空间采样定理，不能产生空间假频，这意味着混合度n不能够无限制地增大。

根据远场理论，在共炮域，一个频率为f的平面简谐波入射到地面检波器上，空间假频出现的截止频率为

$ {f_{\max }} = \frac{v}{{2\Delta x\sin \theta }} $

(36)

式中：v为近地表地震波的传播速度；Δx为检波器

间距；θ为地震波的传播角度。根据波场互换原理，把共检波点道集中的检波点视为炮点，炮点视为检波点，则允许截止频率为

$ {f_{\max }} = \frac{v}{{2\Delta s\sin \theta }} $

(37)

在同时震源数据采集中，相邻n个单炮几乎同时激发，形成一个超级炮，超级炮的间距为

$ \Delta S = n\Delta s $

(38)

由全部超级炮形成的混合波场中，允许的截止频率为

$ {f_{\max }} = \frac{v}{{2\Delta S\sin \theta }} = \frac{v}{{2n\Delta s\sin \theta }} $

(39)

实际震源激发的频率f应不大于截止频率

$ f \le \frac{v}{{2n\Delta s\sin \theta }} $

(40)

可得

$ n\Delta s \le \frac{v}{{2f\sin \theta }} $

(41)

当满足式(41)时，空间假频就不会出现，才能使用γ(x_i, ω)，混合震源才会被正确解混为单个震源。

2 检波点与炮点的不规则排列

在前面原理部分，待分离的同时震源数据体是通过固定炮间距和固定检波点间距的观测系统混合采集得到。在实际生产中，遇到无法避开的障碍物时，不可避免地要改变检波点和炮点间距，形成不规则排列。为了使直接分离算法得到广泛应用，对检波点和炮点都不规则排列的采集条件进行研究。

2.1 检波点不规则排列

式(4)中每一个频率分量的待分离矩阵P_sim(如图 1)，每一个行向量表示一个待分离的共检波点道集。每一行的待分离共检波点道集与矩阵$\mathit{\boldsymbol{ \boldsymbol{\hat \varGamma} }}$_inv相乘后得到分离后的共检波点道集。以第k行共检波点数据为例，其计算公式为

$ {{\mathit{\boldsymbol{\hat P}}}_{k,.}} = {\mathit{\boldsymbol{P}}_{{\rm{sim}},k,.}}{{\mathit{\boldsymbol{ \boldsymbol{\hat \varGamma} }}}_{{\rm{inv}}}} $

(42)

其中$\mathit{\boldsymbol{\hat P}}$_{k, .}和P_{sim, k, .}分别是矩阵$\mathit{\boldsymbol{\hat P}}$和矩阵P_sim的第k行元素组成的一个行向量。式(42)说明基础点扩散矩阵R和点扩散矩阵R_sim都作用于共检波点域，且防止出现空间假频的分离适用条件(式(41))也是在共检波点域中推导得出，并与检波点间距无关，证明矩阵P_sim中各行之间没有关联性，因此检波点间距能够任意取值。即使在检波器个数减小到1，即当矩阵P_sim只含一个混合共检波点道集的行向量时，用该方法仍然能够得到分离后的波场数据。

2.2 炮点不规则排列

在信号的重建中，如果空间采样间隔是固定的，则利用sinc插值函数可对离散信号进行重建，插值函数sinc的离散化形式类似式(30)。同理，在共检波点域，将与sinc插值函数类似的基础点扩散函数$\gamma\left(x_{i}, \omega\right)=\frac{\sin \frac{|\omega| x_{i}}{v_{\mathrm{a}}}}{\pi x_{i}}$离散为式(30)的形式，可以通过点扩散矩阵对固定炮间距下采集的混合地震信号进行分离。如果炮点不规则排列，可对函数γ(x_i, ω)的离散方程(式(30))做适当的调整

$ \mathit{\boldsymbol{R}} = \left( {\begin{array}{*{20}{c}} {{\gamma _{{d_1} - {d_1}}}}&{{\gamma _{{d_2} - {d_1}}}}&{{\gamma _{{d_3} - {d_1}}}}& \cdots &{}&{}\\ {{\gamma _{{d_1} - {d_2}}}}&{{\gamma _{{d_2} - {d_2}}}}&{{\gamma _{{d_3} - {d_2}}}}&{{\gamma _{{d_4} - {d_2}}}}& \cdots &{}\\ {{\gamma _{{d_1} - {d_3}}}}&{{\gamma _{{d_2} - {d_3}}}}&{{\gamma _{{d_3} - {d_3}}}}&{{\gamma _{{d_4} - {d_3}}}}&{{\gamma _{{d_5} - {d_3}}}}& \cdots \\ {}&{}& \vdots &{}&{}&{} \end{array}} \right) $

(43)

式中d_i表示第i炮到原点的距离。上式将γ(x_i, ω)中的固定排列间距用实际采集的不规则排列间距替代。然后将式(43)代入式(33)中，可求得式(4)中分离后的矩阵$\mathit{\boldsymbol{\hat P}}$，最后对每一个频率分量的混合数据进行相同处理，通过傅里叶逆变换得到最终的分离波场。

在炮点不规则排列的观测系统中，防止出现空间假频的分离适用条件(式(41))变为

$ \Delta {s_1} + \Delta {s_2} + \Delta {s_3} + \cdots + \Delta {s_n} \le \frac{v}{{2f\sin \theta }} $

(44)

式中Δs_i为第i炮与第i+1炮之间的距离，Δs₁+Δs₂+Δs₃+…+Δs_n代表一个超级炮在地面的覆盖范围。

3 算列 3.1 规则排列下的模拟数据

图 4为传播速度具有横向和垂直变化的Marmousi模型。模型的最低速度为1500m/s，正演的震源截止频率为25Hz。在地表创建765个炮点，炮间距为10m，以3个相邻单炮组成一个超级炮，即混合度n=3，因此共有255个超级炮，每个单炮的延迟时间是随机的，取值为0~2s；用767个检波器接收信号，道间距为10m，采样间隔为4ms。模拟计算在Intel Xeon E7-4807、内存为128G的IBM服务器上运行，使用的编程语言为MATLAB。混合数据体(其中部分道集如图 5所示)经过一次处理可以得到全部解混后的数据。

如果仅采用最小二乘法进行分离，只对主炮进行了归位。图 5e为在共检波点道集中用最小二乘法求得的分离结果。对比图 5d，可以看出主炮已经归位，但是其他干扰炮形成的“串扰”特征依然存在，这些“串扰”数据携带其他干扰炮的信息，并不能简单地理解为噪声。所以在每个频率分量中，首先构建基础点扩散矩阵R(图 6)，然后计算点扩散矩阵R_sim。在共检波点域，矩阵R_sim把这些“串扰”在时间上按各炮的延迟时间进行归位，在空间上按炮间距进行归位，因此混叠的同时震源数据得到了分离。在实际计算中，为了提高矩阵R_sim求逆的稳定性，一般在矩阵R_sim中增加正则化参数β，即R_sim=Γ^†RΓ+βI，试验证明β=1.0×10^-6ε是一个合理的选择，其中ε是矩阵R_sim中的最大元素。

根据图 3完成后续步骤，得到全部765炮的分离数据。图 5中混合采集数据的分离结果见图 7。

定义如下信噪比公式定量描述数据分离的效果

$ {\rm{SNR}} = 10\lg \frac{{{{\left\| \mathit{\boldsymbol{P}} \right\|}^2}}}{{{{\left\| {\mathit{\boldsymbol{P}} - \mathit{\boldsymbol{\hat P}}} \right\|}^2}}} $

(45)

本次实验中，共炮点数据(图 5a与图 7a)和共检波点数据(图 5b与图 7b)的分离信噪比SNR分别为37.8dB和35.3dB。抽取第90道和第200道的原始数据(图 5b)与分离数据(图 7b)进行对比，结果见图 8。

在前述的计算机环境下，一次性分离全部混合数据体(由255个超级炮，767个检波点和每道1200个采样点组成)耗时381s，用时较短。对比分离后数据(图 7a、图 7b)与原始数据(图 5a、图 5b)，可以看出，不仅浅层数据得到了很好的恢复，深部细节也得到保留，残差较小，信噪比较高，且两条波形曲线的吻合度很高，这证明了本文方法的准确性。

3.2 规则排列下的实际数据

实际数据的炮间距为20m，采样间隔为1ms，观测时间为2.3s，混合度n=2。考虑到截止频率对本方法的影响较大，对实际数据滤波，去掉了高频分量。分离结果及残差见图 9c、图 9d。可以看出，虽然中间记录道的强“串扰噪声”没能完全被压制，但是分离后的深部反射层的同相轴信号得到了较好的恢复，最终的SNR为10.5dB。

3.3 不规则排列的模拟数据

不规则排列的正演速度模型仍为Marmousi模型(图 4)。炮点数为765，炮间距随机，但是需要满足分离适用条件(式(44))。检波点数为150，分为3段，每段包含50个检波点，分别布置在地面500~1500m、3000~4000m和5500~6500m，道间距也随机取值。炮点与检波点覆盖范围如图 10所示，混合度n=3，其他条件与前述规则排列正演相同。分别讨论当检波点或炮点不规则排列下，同时震源数据的直接反演分离效果。

首先，针对检波点不规则排列得到的同时震源数据进行分离测试(图 11)，分离后的信噪比SNR=39.5dB。然后，针对炮点不规则排列的情况进行了测试(图 12)，分离的信噪比SNR=38.1dB。对比分离结果(图 11c、图 12c)与原始采集数据(图 11a、图 12a)，可以看出，混合波场从浅层到深层都得到了很好的分离，在分离的残差(图 11d、图 12d)中几乎看不到连续的有效反射信号。证明本文方法不仅适用于规则排列采集的同时震源数据的分离，通过拓展后同样适用于不规则排列采集的同时震源数据的分离。

4 讨论

本文在求取矩阵R_sim=Γ^†RΓ+βI的逆的过程中采用了正则化方法，该求逆过程会引入较多的假频噪声。通过对正则化参数的研究，选择合适的参数β，可以将这种影响降到最低。图 13展示了取不同正则化参数β时，规则排列的Marmousi混合数据体中，第100个检波器记录的共检波点道集的分离信噪比变化曲线。可以看出当参数β取值较大或者较小时，信噪比都会降低，假频噪声增强，分离效果较差；当β取1.0×10^-4ε~1.0×10^-8ε时，可以得到较高的信噪比。

算例中，对于检波点与炮点规则排列采集的Marmousi正演混合数据体，分离后可以得到765个共炮点道集和767个共检波点道集，这些共炮点道集与共检波点道集各自的分离信噪比见图 14。可以看出，靠近观测系统中部采集的混合数据分离效果较好，处于两端的分离效果相对较差，具体原因还需进一步研究。

5 结论

一般而言，进行高密度震源采样时，若震源间距较小、激发的延迟时间很短，激发波场容易产生强相干性，这不利于同时震源数据的分离。本文采用的直接反演方法则不受这些条件限制，对分离激发延迟时间短、炮间距较小的高密度同时震源数据具有一定的优势。该方法只需在频率域进行一次性分离，不需要迭代，效率高、准确度高，对同时震源数据的分离有很好的应用前景。

通过适当改变基础点扩散矩阵，该方法被推广到检波点和炮点不规则排列的观测系统中，不仅满足陆地复杂地形的同时震源数据的分离，而且适用于海上拖缆的动态不规则采集的混合数据的分离。

需要注意的是，直接反演分离算法有一定的适用范围，要求混合度、炮间距和地震信号的频率满足本文中的分离适用条件。