传统陆上和海上地震勘探中, 为了避免地震记录之间相互干扰, 通常采用逐炮激发的方式, 且在相邻炮之间设置足够的激发时间间隔, 这易造成采集周期过长, 成本过高的问题。近年来发展的同步震源采集技术可以以较小的时间间隔激发处于不同空间位置的多个震源, 接收器连续接收相互干扰的混合炮记录^[1-3]。该方法可以大幅降低采集时间, 显著增加震源覆盖密度, 从而提高地下照明度和地震数据的质量, 因此在地震勘探中, 尤其在海上地震勘探中受到广泛关注。

目前常规的地震资料处理及解释方法均基于单炮记录, 因此将同步震源混合记录分离为传统的单炮记录是进行后续资料处理的关键。基于滤波的方法是将共接收点道集或者其它道集内的混合干扰近似为随机噪声, 利用频率-波数域倾角滤波、空间-时间域中值滤波等去除噪声的方法恢复道集记录^{[1, 4-6]}。这类方法具有计算速度快、易实现的优点。然而, 当混合采集系统中有多个震源混合时, 则会产生较强的干扰噪声, 有效信号易受到损伤, 进而影响分离效果。基于稀疏反演的方法是将多震源混合地震记录作为观测数据, 待恢复的有效地震记录作为模型变量, 将混合记录分离作为反问题, 然后利用稀疏促进的迭代反演方法进行求解^[7-14]。与滤波类方法相比, 基于反演方法分离的有效地震记录具有更高的信噪比和保真度。这类方法还可以与滤波类方法结合进一步用于提高混合记录的分离精度^[13]。

基于稀疏反演的同步震源混合地震记录分离方法的核心是字典基的选取, 目前采用的主要变换有Radon变换^[7]、Curvelet变换^[9-10]和Seislet变换^[11]等。这类基具有较为明确的表达式和几何特征, 但是形式较为固定, 因此很难用于对复杂地震信号进行最优稀疏表示。近些年发展起来的字典学习方法基于机器学习得到自适应于分析信号的字典, 使其能够更好地匹配待分析信号, 分解系数更为稀疏, 因此有望实现地震信号的更优稀疏表示^[14-17]。本文基于机器学习构造出能够对有效地震记录进行最优稀疏表示的学习型字典; 然后在同步震源地震采集和稀疏反演的理论框架下, 将待恢复的有效地震记录基于学习型字典的稀疏表示作为正则化约束, 形成一种同步震源地震数据分离方法; 最后用复杂模型和实际资料验证了算法的有效性。

1 算法原理

将同步震源混合地震数据的分离当作一个反问题, 该反问题通常为欠定的, 因此需要额外的先验信息进行约束, 以期得到较为满意的分离结果。本算法主要基于学习型字典构造稀疏约束项, 通过迭代算法进行求解。

1.1 字典学习

稀疏表示理论通常假定信号可以用少数几个字典原子的线性组合来近似表示。基于块结构字典学习理论, 利用机器学习能够训练一个用于对信号进行稀疏表示的字典。相比目前常用的固定变换(字典), 该学习型字典是数据驱动的, 能够自适应于分析信号, 从而提供信号的更优稀疏表示。字典学习的基本步骤是给定一组固定大小的训练信号样本集, 将该信号集的稀疏表示作为反问题, 采用交替迭代更新的策略更新字典和稀疏表示系数。经过数次迭代后, 我们可以得到学习型字典, 进一步用于非训练样本的稀疏表示。基于L0范数稀疏约束的字典学习问题可以表示为^[16]:

$ \begin{array}{l} \;\;\;\;\;\;\;\;\;\mathop {{\rm{min}}}\limits_{\mathit{\boldsymbol{D, A}}} \left\| {\mathit{\boldsymbol{Y}}{\rm{ - }}\mathit{\boldsymbol{DA}}} \right\|_2^2\\ {\left\| {{\mathit{\boldsymbol{\alpha }}_i}} \right\|_0} \le {\theta _0}\;\;\;\;\;\;\;i = 1, 2, \cdots, L \end{array} $

(1)

式中:Y=[y₁, y₂, …, y_L]∈R^m×L包含L个信号; m表示向量化训练信号的维数; 字典D=[d₁, d₂, …, d_p]∈R^m×p包含p个字典原子, 通常需要对每个原子进行模值归一化以保证算法的收敛性; A=[α₁, α₂, …, α_L]∈R^p×L为分解系数矩阵; ‖α_i‖₀为分解系数向量的L0范数, 表示该向量中不为零的元素个数; θ₀为一正整数。在本文研究的地震信号分离问题中, 我们基于共炮点和共接收点记录的相似性, 从未混合的共炮点数据中抽取训练样本集进行字典学习, 进一步将该字典用于分离共接收点记录。公式(1)为非凸问题, 通常利用交替迭代优化方法进行求解, 本文采用K奇异值分解(K-singular value decomposition, K-SVD)的方法更新得到学习型字典^[17]。

1.2 基于稀疏约束的同步震源数据分离

若采用随机放炮或者随机时间抖动的同步震源地震采集技术, 则对一个固定位置的接收器, 接收到的混合地震记录可表示为^[9]:

$ \mathit{\boldsymbol{b}} = {\mathit{\boldsymbol{T}}_{\rm{r}}}\mathit{\boldsymbol{f}} $

(2)

式中:向量f表示向量化的常规共接收点道集, 若连续放炮N_s次, 常规记录中每道包含有N_t个时间采样点, 则f的维数为N(N=N_sN_t); 矩阵T_r表示混合算子, 它的元素由震源的激发方式决定^[9]; 向量b为接收的混合地震记录, 其维数n满足$n \ll N$。

基于稀疏约束的同步震源数据分离反问题可以表示为:

$ \begin{array}{l} \mathop {{\rm{min}}}\limits_{f,\mathit{\boldsymbol{A}}} v\;\left\| {{\mathit{\boldsymbol{T}}_{\rm{r}}}\mathit{\boldsymbol{f - b}}} \right\|\;_2^2 + \sum\limits_{ij} {\left\| {{\mathit{\boldsymbol{R}}_{ij}}\mathit{\boldsymbol{f - D}}{\mathit{\boldsymbol{\alpha }}_{ij}}} \right\|\;_2^2} \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\left\| {{\mathit{\boldsymbol{\alpha }}_{ij}}} \right\|{\;_0} < {\theta _0},\forall i,j \end{array} $

(3)

式中:v为正整数, 其大小与接收数据的信噪比有关, 若接收数据中的随机噪声较小, 则适当增大v的取值, 以减小数据拟合项$\left\| {{\mathit{\boldsymbol{T}}_{\rm{r}}}\mathit{\boldsymbol{f - b}}} \right\|\;_2^2$的值, 反之, 若噪声较强, 则应适当减小v值, 使得$\left\| {{\mathit{\boldsymbol{T}}_{\rm{r}}}\mathit{\boldsymbol{f - b}}} \right\|\;_2^2$接近噪声方差水平; ${\mathit{\boldsymbol{R}}_{ij}} \in {{\bf{R}}^{m \times N}}(m = {n_{\rm{t}}} \times {n_{\rm{s}}})$表示从待恢复的常规地震记录中抽取向量化大小为n_t×n_s的块数据, n_t和n_s分别表示抽取数据块包括的时间和空间点数, i, j分别为块数据在整个共接收点记录剖面上的时间和空间位置; A为分解系数矩阵, α_ij为其列向量。

本文采用交替迭代的方法求解公式(3), 其求解步骤为:

1) 固定f, 此时公式(3)改写为

$ \begin{array}{l} \mathop {{\rm{min}}}\limits_\mathit{\boldsymbol{A}} \sum\limits_{ij} {\left\| {{\mathit{\boldsymbol{R}}_{ij}}\mathit{\boldsymbol{f}} - \mathit{\boldsymbol{D}}{\mathit{\boldsymbol{\alpha }}_{ij}}} \right\|\;_2^2} \\ \left\| {{\mathit{\boldsymbol{\alpha }}_{ij}}} \right\|{\;_0} < {\theta _0}, \forall i, j \end{array} $

(4)

利用正交匹配追踪方法^[17]求解公式(4)可以得到稀疏系数矩阵A;

2) 固定A, 公式(3)变为

$ \mathop {{\rm{min}}}\limits_f v\;\left\| {{\mathit{\boldsymbol{T}}_{\rm{r}}}\mathit{\boldsymbol{f - b}}} \right\|\;_2^2 + \sum\limits_{ij} {\left\| {{\mathit{\boldsymbol{R}}_{ij}}\mathit{\boldsymbol{f - D}}{\mathit{\boldsymbol{\alpha }}_{ij}}} \right\|\;_2^2} $

(5)

该问题为最小二乘问题, 本文采用高效的共轭梯度法进行求解, 进而更新f。上述两个步骤交替迭代进行, 当公式(4)中的数据拟合项小于设定的数值或者迭代达到设定的次数时, 迭代终止, 此时将得到的A和f作为最终结果。

2 实际算例 2.1 参数选取

算例包括一个复杂合成记录和一个实际炮集数据数值合成的混合记录, 均假定为双边采集, 并且接收器固定。本文采用的同步震源混合采集策略为每7炮数据采用随机时间抖动的震源激发方式进行混合采集后, 记录一炮未混合的常规单炮数据, 两者交替进行。

f的初始值选取为T_r^Tb。对于合成地震记录, 选取的地震数据块包括32个时间采样点和16个空间采样点, 字典原子个数为512。对于实际地震数据, n_t和n_s分别为16和8, 字典包括128个原子。参与训练的数据块为20000个, 均从未参与混合的单炮记录中选取。另外, 为了提高计算效率, 可将字典的初始值选取为局部余弦基或者Cuverlet基。θ₀选取为32, v选取为10。在本文的算例中, 由于(未混合的)常规地震记录已知, 因此采用如下的信噪比公式来评价最终的分离效果:

$ {S_{{\rm{NR}}}}\left( {\mathit{\boldsymbol{f}}, \mathit{\boldsymbol{\tilde f}}} \right) = 10{\rm{lg}}\frac{{\left\| \mathit{\boldsymbol{f}} \right\|\;_2^2}}{{\left\| {\mathit{\boldsymbol{\tilde f - f}}} \right\|\;_2^2}} $

(6)

其中, ${\mathit{\boldsymbol{\tilde f}}}$为分离的共接收点地震记录。

2.2 模型数据

基于Marmousi速度模型, 利用有限差分方法生成了256炮地震数据, 每炮256道接收, 每道有1024个采样点, 炮点间隔和接收器间隔均为8m, 时间采样间隔为2ms。图 1a为接收到的部分同步震源混合地震记录, 其中有一炮数据未参与混合, 其余7炮数据的记录相互干扰。抽取若干未混合的单炮地震记录用于字典学习, 得到的部分学习型字典原子如图 1b所示。可以看出, 这些原子反映了地震数据的局部波形特征。图 2给出了基于学习型字典稀疏反演的模拟记录分离结果。其中, 图 2a为混合的共接收点道集, 存在明显的炮间干扰; 图 2b为采用本文方法恢复的共接收点道集; 图 2c为图 2b对应的误差剖面。可以看出, 采用本文方法可以很好地去除干扰噪声, 恢复常规的共接收点道集数据, 且误差剖面上基本没有有效信号残留。基于公式(6)计算的信噪比值为16.7。进一步采用固定的二维局部离散余弦变换(local discrete cosine transform, LDCT)作为字典进行迭代求解, 恢复结果的误差如图 2d所示。对比图 2c和图 2d可以看出, 利用固定的二维LDCT恢复的共接收点道集误差较大, 对应的信噪比值为12.8。

2.3 实际数据合成的混合记录

实际数据包括128炮, 每炮有256个接收点, 炮间隔和检波器间隔分别为25.0m和12.5m, 每道记录的时间样点为512个, 采样间隔为4ms。采用与模型算例同样的同步震源采集方式合成混合地震记录(图 3a), 然后将若干未混合的单炮记录采样为25m间隔的数据, 进而进行字典学习, 得到的部分学习型字典原子如图 3b所示。图 4给出了实际数据分离结果。其中, 图 4a为待分离的共接收点记录; 图 4b为采用本文方法恢复的共接收点道集; 图 4c为图 4b对应的误差剖面; 图 4d为采用二维LDCT方法得到的误差剖面。对比图 4c和图 4d可以看出, 因为实际资料波形的复杂性, 误差剖面上均有部分有效信号残留, 但是图 4d上有较多的有效信号没有被恢复, 因此不能保证恢复结果的精度。

3 讨论与结论

在基于时间随机抖动的同步震源混合采集框架下, 给出了一种基于块状学习型字典稀疏反演的同步震源混合地震记录分离方法。基于K-SVD算法, 利用未混合的共炮点记录训练的学习型字典反映了地震记录的局部波形特征, 能够提供有效地震信号的高稀疏度表示, 进一步基于共炮点和共接收点记录的相似性, 可将其用于稀疏表示待恢复的共接收点道集数据。模型数据和实际数据的恢复结果表明, 相比固定二维LDCT字典基下的稀疏表示, 将待分离记录基于学习型字典的稀疏表示作为正则化约束项, 可以有效抑制干扰噪声, 较好地恢复共接收点记录, 从而提高分离精度。

另外, 为了得到训练样本集, 本文所采用的同步震源采集方式较为特殊, 其采集了若干炮未混合的地震数据。下一步需要针对一般的同步震源采集系统研究混合地震记录的分离方法。