0 引言

数据采集成本与数据处理难度一直在地震勘探领域中博弈。在传统地震数据采集过程中，震源与震源之间的激发时间间隔很大，足以避免炮记录之间产生交叠。这种采集方式虽然可以获得高质量的地震数据，但是时间和经济成本都很高。混采是指将不同位置的震源按照一定的时间编码同时或延时激发^[1]，从而获得互相干涉的混合炮集记录的地震数据采集技术。虽然混采数据的处理难度更大，但因为混采具有相对传统采集技术而言采集效率更高、经济成本更低的特点，近年来广受工业界的青睐。我们将混合炮集记录中由指定炮点激发产生的地震记录定义为有效信号，由其他炮点激发产生的地震记录定义为混叠噪声。混叠噪声的存在会降低地震数据的信噪比，为后续的数据处理造成困难，最终降低成像的质量。

混采数据的处理方法可以分为两类^[1]。第一类是直接对混采数据进行成像，同时将混叠噪声作为约束条件来压制成像过程中出现的串扰噪声。Tang等^[2]提出了针对混采数据的最小二乘偏移算法，可以在提高成像分辨率的同时压制串扰噪声。Verschuur等^[3]将含有自由表面多次波的混采数据直接应用到全波场偏移技术中，扩大了地下成像的照明范围。Anagaw等^[4]采用同时震源激发的策略，成功减少了频率域全波形反演的计算量。

第二类是先将混采数据分离成独立震源激发的地震记录，再对混采分离后的数据应用常规方法进行处理。分离混采数据主要有两种思路，一种基于滤波理论，另一种基于反演理论。基于滤波理论的混采分离方法利用了多震源激发时刻在共享时间间隔基础上存在随机抖动的特性。在实际混采过程中，不同位置的多个震源通常按照线性时间编码同时或者延时激发，且激发时间间隔在某一固定值的基础上存在随机抖动，因此混采数据中的有效信号和混叠噪声只在共炮域中存在运动学上的相关性，而在非共炮域(比如共检波点域、共中心点域、共偏移距域)中具有明显的运动学差异。比如在共检波点域中，有效信号是线性连续的，而混叠噪声表现为随机分布的非线性噪声。基于上述特征，我们可以利用去噪方法保留线性连续信号即有效信号，去除非线性噪声即混叠噪声，从而实现混采分离。Huo等^[5]采用多方向矢量中值滤波算法在共偏移距域对混叠噪声进行压制。张良等^[6]在共检波点域中使用基于时间窗边线和多级中值滤波的方法对混采数据进行分离。Chen^[7]使用动校正和变空间中值滤波法实现了有效信号和混叠噪声的分离。

基于反演理论的混采分离方法也在近年取得了极大的发展。我们将混采数据的采集过程看作正问题，那么其分离过程可以被看作反演问题。这种分离方法的基本理论是，在特定的混采方程中加入适当的约束条件，通过对正演模型算子的矩阵求逆，或者使用迭代框架来迭代估计有效信号并减去混叠噪声，从而解决反演问题，实现混采分离。Mahdad等^[8]在共检波点域采用结合阈值约束的最小二乘算法实现了混采分离。Lin等^[9]将近年来较为流行的压缩感知理论和混采分离进行结合，并在曲波域对混采数据的分离反演问题进行约束。Chen等^[10]提出了一种基于Seislet变换和整形正则化思想的混采数据分离算法。由于混叠噪声的出现会增大Hankel矩阵的秩，该特征也为使用矩阵降秩实现混采分离提供了新思路^[11-13]。值得注意的是，上述方法只能实现混采数据的炮集分离，分离后的地震数据依然需要进行诸如多次波压制、鬼波压制等后续处理。

多次波是地震勘探中的重要相干噪声，其压制效果直接关系到偏移成像质量与地质解释。由于海水表面和海底都是强反射界面，其在海上地震数据中尤为发育。如何在压制多次波的同时将混采数据分离成独立震源记录一直是地震勘探领域具有挑战性的热点问题。Verschuur等^[14]提出的表面相关多次波消除(SRME)方法是一种数据驱动的多次波压制方法，只需原始地震数据就能进行多次波预测，然后按照最小二乘原则将预测多次波与原始数据进行匹配相减，得到一次波结果。该方法在过去二十年中取得了长足发展，是商业软件中的常见模块。SRME方法只适用于传统地震数据，在混采流行后，工业界通常先对混采数据进行分离，然后再应用SRME方法以取得独立震源一次波估计结果。

2009年，van Groenestijn等^[15]提出了基于稀疏反演的一次波估计(EPSI)方法。该方法通过反演的形式对传统非混采数据中的地下一次波响应进行估计，克服了SRME在多次波匹配相减过程中对一次波能量产生伤害的缺点。2014年，van Groenestijn等^[16]改进了原有方法，将混采震源项线性编码引入了传统EPSI的一次波-多次波模型之中，提出了针对混采数据的EPSI方法，实现了混采分离与一次波估计的一体化处理。然而，EPSI理论要求使用者人为确定时窗来选取地下一次波响应的初值，这在地下地质情况复杂时很难实现。此外，在对混采数据应用改进EPSI方法时，由于混叠噪声的存在，人为设计时窗拾取地下一次波响应初值的方法并不可靠，与目标一次波交叠的混叠噪声极易被引入到估计结果中。

针对上述问题，本文在前人研究的基础上：首先，改进了混采数据EPSI方法的目标函数；其次，在估计地下一次波响应的过程中引入基于L1范数的双凸优化理论，并用基于L1范数的谱投影梯度(SPGL1)算法求解上述反演问题，避免陷入局部极值；最后，在地下一次波响应的估计过程中引入二维曲波变换和一维小波变换的联合作为三维稀疏约束的条件，在确保求解精度的前提下提高反演求解的速度。我们将本文所述的混采分离与一次波估计方法应用于模拟混采数据和海上实际混采数据，并与传统算法进行比较，以验证本文所述方法的有效性和优越性。

1 算法原理和实现 1.1 传统的混采数据EPSI方法

由传统地震勘探方法采集得到的地震波场包括一次波和表面相关多次波两部分。地震波上行波场、一次波和表面相关多次波在频率域中的关系可表示为^[14]

式中：P为采集得到的地震波上行波场；X₀为地下一次波响应；S⁺为震源子波函数；R为地下反射系数。X₀S⁺表示由X₀和S⁺在频率域相乘获得的一次波波场。值得注意的是，本文所指的一次波反射点不在地表，即一次波波场为没有在地表发生过反射的波场，包括真实一次波和层间多次波。将X₀、R、P进行多维褶积即可获得表面相关多次波X₀RP。

非混采情况下地震波场各元素之间的物理关系如图 1所示，地震波下行波场是S⁺和P被地表反射而向下传播所产生的波场的和。由此，P还可以表示为X₀和下行总波场在频率域的乘积：

混采数据是由两个及以上震源根据线性编码同时或延时激发所获得的地震记录^[1]。混采情况下的地震波总波场、一次波和表面相关多次波之间的关系如图 2所示。

由此，混采情况下的一次波-多次波模型关系可以表示为

式中：P_bl为混采波场；S_bl为所有混合震源的线性组合^[1]。

为了表示所有混合震源的线性编码过程，我们引入混采算子

由此可得

式中：ω为地震波场的频率；Δt₁和Δt₂表示不同震源之间的激发延迟时间。同理，混采数据也可以被表示为若干传统地震数据在频率域中的线性编码组合：

引入Γ_bl后，混采数据的一次波-多次波模型可以被表示为将非混采数据的一次波-多次波模型(式(2))等号两边同乘混采算子的形式：

对比式(7)和式(3)可知，混采波场也可以被表示为地下一次波响应和混采下行波场的乘积。注意，式(7)中的X₀和R依然是非混采情况下的地下震源子波和反射系数序列^[17]。

由于式(7)中同时存在X₀和S⁺两个未知数，我们利用交替反演的策略来对二者同时求解，进而对整个混采波场P_bl进行分解，继而得到目标函数J：

式中：J_i为第i次迭代的目标函数；和分别为通过反演获得的X₀和S⁺的估计值。此处假设R=- I(I为维度与P一致的单位矩阵)。

将和的初始值均设置为维度与P一致的0矩阵。在迭代的过程中，首先使用最速下降法的策略对进行更新：

式中：表示实际混采波场与估算波场之间未分解的残差；H表示复共轭转置。由于和的初始值设为0矩阵，因此在第一次更新的过程中，式(9)的计算结果等于混采波场与其自身的多维反褶积，即P _bl(RPΓ _bl)^H。

由式(9)可知，实际波场和估计波场之间的残差可通过两种途径更新到中，一是以混采数据一次波的形式，即，二是通过将残差转化为多次波的形式，即。

为稳定反演过程，获得准确的一次波响应估值，在每一次更新过程中，都在时间域选择时窗以尽可能多地选取强一次波响应，即在选定同相轴的时候尽可能只选取一次波同相轴、排除多次波同相轴，再将选取的一次波响应引入到的迭代更新过程中。然后，将更新好的加入到地下一次波响应的估计值中：

式中，α表示每一次更新的步长，确保目标函数残差随着迭代次数的增多而逐渐减小。

的更新过程与的更新过程原理相同，即

通过上述方式对X ₀和S ⁺交替更新，直到目标函数的残差满足设定要求，即可获得地下一次波响应X ₀和震源子波S ⁺，再将二者进行褶积便能获得独立震源一次波响应估计结果。

混采数据EPSI方法在输入含多次波的原始混采数据后，可以直接输出独立震源地下一次波响应估计值，能够实现混采分离和一次波估计的一体化处理，与先进行混采分离再进行多次波压制的非一体化处理方法相比，处理效率明显提高，为后续的数据处理流程节省了时间和经济成本。但是，正如前文所提，EPSI理论有其无法避免的局限性，在应用混采数据EPSI方法时，我们必须在的每一步更新过程中采用确定时窗的方法来选取同相轴，以稳定反演过程。因此，在地下地质情况复杂时，我们会因为难以选取合适的同相轴而无法保证反演结果的准确性。

1.2 基于三维稀疏反演的混采数据一次波估计方法

首先，我们将混采数据的一次波-多次波模型改写为线性算子相乘的形式，此时混采波场P _bl可表示为

式中：P _bl和X ₀分别为P _bl和X ₀的列向量；F_t和F_t^*分别为正、反傅里叶变换算子；BlockDiag_ω表示频率域分块对角矩阵；*表示矩阵的共轭转置；⊗为Kronecker乘积算符。式(13)表示地下一次波响应和上行波场在频率域进行多维褶积获得混采波场的物理过程。

为了提高反演的稳定性，将基于L1范数的双凸优化反演理论^[18]引入式(13)的求解中：

式中：为的列向量；σ为目标函数的跳出值，在反演过程中，一般将其设置为输入数据L2范数的1%~10%。当反演的目标函数满足设定要求后，即表示完成反演过程。式(14)中，目标函数和约束条件都是凸集合，这就确保了在求解过程中只会获得一个解，即全局最优解。因此，同传统混采数据EPSI方法采取的最速下降法反演策略相比，本文采取的基于L1范数的双凸优化反演策略具有明显的优势。

选用基于L1范数的谱投影梯度(SPGL1)^[19-21]算法求解上述反演问题。同时，在反演过程中引入三维稀疏变换作为约束条件，以提高最终反演结果的精度：

式中：t ₀为三维稀疏域的稀疏系数；T ^*为三维稀疏变换的逆算子，将数据由稀疏域变换回时间域。在三维稀疏域，一次波通常对应较大的稀疏系数，而多次波通常对应较小的稀疏系数；我们在反演过程中，通过对σ进行设置，保留较大的系数(一次波)，去掉较小的系数(多次波)，从而达到压制多次波的目的。

本文中采用的三维稀疏变换为二维曲波变换和一维小波变换的联合稀疏变换，即在炮检域使用二维曲波变换，在时间域使用一维小波变换。传统三维稀疏约束通常采用三维曲波变换，这种约束方法虽然可以获得高精度的求解结果，但是运算效率非常低，耗时很长。相较之下，本文采用的三维联合稀疏变换不仅能够保证较高的求解精度，而且可以大大提高求解的速度，达到二者兼得的目的^[22-23]。连续使用SPGL1算法对式(15)进行求解，即可获得地下一次波响应估计值，将它和子波估计值进行褶积，就能获得混采分离后的一次波估计值，达到一体化压制多次波和混采分离的目的。

具体算法流程如下：

Ⅰ.根据混合震源激发时差构建混采算子Γ _bl；

Ⅱ.输入混采数据P _bl并计算σ(输入数据L2范数的1%~10%)；

Ⅲ.用线性算子相乘的形式表示混采数据一次波-多次波模型；

Ⅳ.将和的初始值设置为0，通过式(13)求取地下一次波响应估计值；

Ⅴ.构建三维稀疏变换正算子T以及反算子T ^*(二维曲波变换-一维小波变换)；

Ⅵ.根据式(3)，基于最小二乘反演原则，获得S ⁺的估计值；

Ⅶ.将求得的震源子波项以及三维稀疏变换正反算子带入式(14)，使用SPGL1算法求解，获得地下一次波响应估计值；

Ⅷ.重复步骤Ⅴ—Ⅷ，直到实际混采数据与估算混采数据的残差小于σ；

Ⅸ.将最终获得的地下一次波响应估计值与震源子波项褶积，获得独立震源一次波估计结果。

2 应用实例 2.1 理论数据试算

为了验证本文方法的有效性，首先采用理论数据进行试算。理论数据试算的优点在于可以事先模拟出真实的地下一次波响应，这一点是实际数据试算无法做到的。通过将本文方法的处理结果与真实地下一次波响应进行对比，我们能够精确地评价本方法的处理效果。所有试算均会将本文方法与传统混采数据EPSI方法进行比较，以证明本文方法的先进性和优越性。

模拟放炮的观测系统中共有150个检波器，检波器间距为20 m，检波器采样时长为2 s，采样间隔为0.004 s，所用盐丘模型如图 3所示^[24-25]。放炮时，震源船1由第一个检波器处出发，驶向中间检波器，震源船2从中间检波器处出发，驶向最后一个检波器，且震源船2的放炮时间相对震源船1的放炮时间存在±0.5 s的随机抖动。将由上述观测系统采集得到的两例混采地震记录示于图 4中，图 4a、b中的目标炮记录分别为震源船1放炮顺序中的第20炮和震源船2放炮顺序中的第110炮，换言之，我们将由震源船2激发产生的炮记录视为混叠噪声。在传统地震数据采集系统中，两震源之间的激发时间间隔足够大，可以确保激发产生的两炮地震记录之间不存在同相轴的交叠。而本文中的两艘模拟震源船几乎同时放炮，其激发时刻之间只存在非常细微的差异，因此采集得到的地震记录表现为两炮记录相互交叠的形态，尤其是旅行时间1.2 s的部分，两炮记录完全重叠(图 4)。如果直接对该混采数据进行常规地震数据处理，混叠噪声(邻炮记录)会降低有效信号(目标炮记录)的信噪比，影响最终地震数据处理的精度。

为了准确评价本文方法与传统方法估计得到的地下一次波响应，我们模拟了上述两例混采数据的真实一次波，示于图 5中。并分别利用传统混采数据EPSI方法和本文方法进行混采分离与一次波估计，所得结果展示于图 6、图 7之中。可以看出，同地下真实一次波(图 5)相比，利用传统混采数据EPSI方法获得的独立震源一次波估计结果深部区域(1.2~2.0 s)存在许多由反演造成的假象(图 6)，而利用本文方法获得的独立震源一次波估计结果更接近真实一次波，尤其是深部区域的反演假象明显更少(图 7)；这证明本文方法能够更准确地反演出独立震源一次波响应。然而，不论是传统方法还是本文方法，都存在邻炮信息于反演过程中泄露的问题，导致部分混叠噪声残存在最终的独立震源一次波估计结果中。

2.2 实际数据试算

为了进一步验证本文方法的有效性和工业使用价值，我们选取某海域的海上实际数据进行试算，并与传统混采数据EPSI方法进行对比。该混采数据对应的观测系统中有两个几乎同时激发的震源，分别激发了170炮，系统中共有检波器170个，检波器记录时间为4 s，采样间隔为0.004 s。

海上实际混采数据如图 8a所示，与上一节中的模拟混采数据相比，实际混采数据的同相轴组成更加复杂，邻炮干扰更加强烈，两震源产生的炮记录自旅行时间0.8 s处开始产生交叠，在旅行时间大于1.8 s的区域中全部重叠在一起。左侧炮记录为有效炮记录，右侧炮记录为混叠噪声，处理目的是在分离消除混叠噪声的同时估计出有效炮记录中的一次波。图 8b、c分别为利用传统混采数据EPSI方法和本文方法获得的一次波估计结果。

通过对比图 8a、b、c可知，本文方法处理结果中的深层同相轴信息更加连续，远偏移距处一次波信息更加完整。将1.2 ~1.4 s、前20个检波器所示区域放大(图 9)后可以看出，本文方法局部假象获得了更好的压制(黑色箭头所示)，本文方法已经达到工业化的要求，具有一定的工业应用价值。

3 结论与展望

1) 本文将传统混采数据EPSI方法的目标函数改写为线性算子相乘的形式，将基于L1范数的双凸优化理论引入到混采分离与一次波估计方法的求解过程中，避免了传统方法反演过程陷入局部极值、需要人为调整时窗和选定同相轴的问题，稳定了求解过程，确保最终求取的一次波地下响应估计值为全局最优解。

2) 在反演过程中引入了三维联合稀疏约束，相比于传统混采数据EPSI方法，本文方法获得了更准确的一次波求解结果，证明了本文方法的先进性和优越性。

3) 本文方法对混采数据进行混采分离和一次波估计的处理结果与真实地下一次波响应十分接近，可以应对复杂的实际地质情况，具有向工业界进行推广的应用价值。文中提及的邻炮信息在一次波反演过程中泄露的问题，将会是笔者下一步的研究方向。