基于压缩感知的地震观测系统设计及数据重建

2022-11-04崔庆辉芮拥军秦宁刘鹏翔隆文韬

科学技术与工程 2022年27期

崔庆辉，芮拥军，秦宁，刘鹏翔，隆文韬

(中国石油化工股份有限公司胜利油田分公司物探研究院，东营 257022)

随着社会经济的快速发展，当前在东部老油区城镇、工厂、村庄、道路等障碍物密集分布，地震勘探的地表条件发生了巨大变化，造成基于规则采样理论即Shannon采样定理[1]设计的观测系统受密集障碍物影响，在实际布设时点位发生严重偏离，特别是炮点位置变化更大，呈现出随机分布特点，导致成像结果偏离预期。另一方面，一些老油田勘探开发程度较高，需要解决的地质目标表现出“深、小、薄、碎”等特征，在常规规则采样理论下只能不断提高采样密度即炮道密度，这无疑会造成采集成本的急剧上升。

压缩感知理论是一种新的采样理论，基于信号的稀疏特性，在远小于Nyquist采样率的条件下，通过随机采样获取信号的离散样点值，然后进行信号的重建[2-3]。压缩感知理论一经提出，就引起学术界和工业界的广泛关注，在各个领域得到迅速发展。

国外较早将压缩感知理论应用于海上地震勘探，取得了较好的应用效果[4-5]。马坚伟[6-7]较早研究了压缩感知理论在地震勘探领域应用的可行性，周松等[8]和吕公河等[9]在西部沙漠地区开展了基于压缩感知的地震采集试验，为中国陆上基于压缩感知的地震勘探实际应用奠定了基础。黄小刚[10-11]提出了一套海上压缩感知采集设计及数据重构方法，并用模型数据展示了海上压缩感知采集的效果。总之，在当前低油价背景下，压缩感知理论在地震勘探领域的研究得到了广泛的重视[12]。

利用压缩感知理论实现随机地震采样并重构规则数据，首先要分析地震数据的稀疏性，地震信号(或波场)的时空复杂度(与地下及地表速度结构、噪音特征等相关)会显著影响地震信号的稀疏性，进而导致很多情形下随机采集的实际数据的重建效果无法达到规则采样数据的质量[13]。中国陆上油田地表噪声发育，地下波场复杂，基于压缩感知理论进行随机地震数据采集及规则化重建具有更大难度。

首先研究基于压缩感知理论的随机地震观测系统设计方法及基于地震数据高维空间低秩约束的随机地震数据规则化重建方法，并通过三维模型数值模拟证明了以上方法在同等采样密度下比常规规则采样具有更好的成像效果。

1 随机地震观测系统设计基本理论

随机地震采集需要满足三个条件：地震数据具有稀疏性、合理的随机采样方法和有效的数据规则化重建算法。关于地震数据的稀疏性，很多学者的研究表明实际地震数据在时间域或其他域具有稀疏特征，即地震数据可稀疏表达[13-15]。现重点研究随机地震观测系统设计及数据规则化方法。

1.1 随机地震采样矩阵

野外地震数据观测实际上是对信号的一种采样，地震数据采集过程可以表示为

d=Df

(1)

式(1)中：f∈Rm为原始地震数据；d∈Rn为经过采样得到的地震数据；D∈Rn×m为采样矩阵；通常有n

若f本身不稀疏，则可以使用f在变换域中的稀疏性来得到一种等价的表示。这里用C表示f的稀疏变换，则式(1)就可以转换为

d=Ak

(2)

A=DCH

(3)

式中：H代表反变换矩阵；k∈Rn为f在稀疏域中的系数；n为变换域的维度(若C为正交变换，则n=m；若C为冗余变换，n>m)，A∈Rn×n。

要满足A为随机矩阵的要求，在变换矩阵C确定的情况下，采样矩阵R也必须是随机分布矩阵，那么就可以利用采样数据高概率地恢复完整地震数据。

1.2 随机地震观测系统设计方法

通过高斯随机函数进行完全的随机采样，被基于压缩感知的技术广泛采用。但是，这种方法容易造成采样点高度聚集或者高度分散的情况，导致对信号的某些部分采样冗余或者采样不足，不仅不符合高精度地震数据采集的要求，而且也会给后续数据恢复处理带来困难，难以达到理想的重建效果。

计算机图像学领域应用的泊松碟随机采样方法[16]、抖动采样方法(Jittered)[17]以及分段采样方法不仅能保证测点的随机性，还能保证测点的空间均匀分布，并具有明显的“蓝色噪声”频谱特征[14]。“蓝色噪声”频谱特征是描述随机类型特征的一个统计模型，其特点是在主要频率周围分布的其他频率值都非常小，所以真值或者主频率就很容易被检测到。

结合实际地震数据特点，基于Jittered采样原理提出了以下随机地震观测系统设计技术流程。

步骤1按照常规的三维地震观测系统设计方法确定理想的观测参数，在正交观测系统下可以用炮点和检波点的点距、线距、点数和线数表示，即{x0,y0,Δx,Δy,nx,ny}，其中(x0,y0)为采样点的起始坐标，Δx和Δy分别为采样点之间的横向和纵向间隔，nx和ny分别为横向和纵向采样点数。基于以上参数可得到规则的采样点坐标{x,y}。

步骤2根据工区地表及地下复杂程度设置纵横向压缩采样比例{rx,ry}，基于该压缩采样比例计算随机采样点坐标的期望值{xe,ye}。

步骤3设置随机采样点坐标相对于期望值的扰动量，即高斯随机分布的标准差(σx,σy)，一般取(σx,σy)=(3Δx,3Δy)。

最后，根据步骤2中期望值，以及步骤3中标准差利用随机函数生成随机观测网格{xrand,yrand}。

2 高维空间低秩约束的随机地震数据重建

地震数据的偏移对数据规则性要求很高，实际采集的地震数据受观测系统、地表变观等影响需要进行规则化处理，因此产生了很多的地震数据插值算法[18-19]，主要包括基于函数变换的重建、基于波场延拓的重建、低秩约束的地震数据重建。低秩约束的地震数据重建算法应用最为广泛，其基本原理为：完整地震数据在频率域可以被一个低秩矩阵或者数组表示，不规则缺道数据被看作随机噪声，而随机噪声会增加矩阵或者高维数组的秩。地震数据重建就可以归结为一个高秩矩阵或者数组的降秩问题，通过降秩实现重建。

地震数据中的反射波同相轴在局部表现为线性特征，利用数据的局部相似性，抽取局部五维数据，在高维空间进行“稀疏-低秩联合反演重建”，从而在压缩感知理论下保证了重建数据的精度，算法主要流程如图1所示。

图1 低秩约束的高维空间地震数据重建算法流程Fig.1 Low-rank constrained high-dimensional spatial seismic data reconstruction algorithm flow

步骤1定义一个理想的规则工区网格G，按照规则观测系统设计方法[20]在网格G上定义规则观测系统，记为Ωreg。

步骤2获取原始输入数据集Dirr的非规则地震数据的观测系统。以炮集为单位获取炮检点坐标，取整到临近网格，得到Dirr对应的非规则观测系统Ωirr。

步骤3以窗口滑动方式从原始数据集取得一个局部高维数据子集。通过在炮端检端分别取局部窗口，抽取原始道集中落在窗口中的所有地震道，从而获得一个不完备的五维局部数据体。通过滑动中心炮检点，可以循环获取不同位置的局部高维数据。

步骤4在步骤3获取的不完备五维局部数据体上，在频率域进行规则数据重建，输出一个完备且规则的局部五维局部数据体。计算过程即求解优化问题，即

(4)

该优化问题实现了对重建结果的三种约束。

(1)要求在观测点位上的重建结果与原始观测相吻合(L2项)。

(2)要求重建数据在四维波数具有稀疏性(L1项)。

(3)要求重建数据在进行Hankel排列后具有低秩性(核范数项)。

式(4)所示优化问题的求解，可通过以下迭代算法实现。

(1)关于L2项的梯度更新为

(5)

(6)

(7)

(2)关于L1项的软阈值收缩为

(8)

(3)关于核范数项的奇异值收缩：

(9)

式(9)中：H(·)和H-1(·)分别是Hankel矩阵化正、反变换算子；Sμ(x)表示对x进行奇异值软阈值收缩，阈值取μ。

步骤5对步骤4重建的局部五维规则数据体做余弦窗函数加权，然后以累加方式更新到规则数据集Dreg中，并记录输出数据集中各个地震道的累加总权重Csum(xs,xr)。

更新规则数据集之前需要对重建数据的余弦加权操作，即每一道都需要乘以一个标量系数，对于震源位于xs，检波器位于xr的地震道d(xs,ys,xr,yr,t) ，它对应的加权系数C计算公式为

(10)

(11)

同时，对式(10)进行更新，每一次更新之后都需要将各个地震道的加权系数累加到总权重数组Csum中，以便后续消除加权幅值不均匀导致的畸变。更新公式为

(12)

步骤6重复步骤3～步骤5，直到完成遍历整个数据集。通过滑动局部窗遍历数据的过程中，保持数据窗尺寸不变，各方向的滑动步长取相应方向的窗口尺寸的一半。

步骤7通过除以各个地震道的累加总权重Csum(xs,xr)，将计算得到的规则数据集Dreg权重归一化。将重建的规则数据的各个地震道Dreg(xs,xr,t)除以其相应的累加总幅值Csum(xs,xr)，以消除加权幅值不均匀导致的畸变。

3 模型数值模拟

选取东部某油田一个三维模型(图2)，通过炮集数据的数值模拟及叠前深度偏移成像，验证文中随机地震观测系统设计及数据重建算法的有效性。

图2 东部油田某工区三维地质模型Fig.2 Three-dimensional geological model of a work area in the eastern oilfield

首先设计如图3(a)所示的三维规则观测系统，炮点横向间隔50 m。纵向间隔400 m，检波点横向间隔400 m，纵向间隔50 m。考虑到目前实际生产中单方向随机采集可行性最高，因此本次随机观测系统设计对炮点只在横向上按照压缩率50%进行随机采集设计[图3(b)]，对检波点只在纵向上按照压缩率50%进行随机采集设计[图3(c)]。

图3 用于数值模拟的观测系统Fig.3 Geometry used in numerical simulation

利用三维波动方程正演得到模拟的规则炮集和随机炮集数据[图4(a)]，分别加入高斯噪声与椒盐噪声[图5(a)和图6(a)]，利用本文提出的高维空间低秩约束的随机地震数据重建算法得到重建后的规则地震数据，对比重建地震数据的波数谱。在波数谱中可以看出，模型正演地震记录的波数谱展布区域较大[图5(a)]，表明了地震信号的复杂度较高。此外，无噪声及加入不同噪声情况下，随机地震数据重建结果较好且波数谱也得到了较好的恢复，并未在波数谱中引入明显的空间假频(图4～图6)，充分展示了高维空间低秩约束的随机地震数据重建的性能。

将规则的地震数据按照地震道均匀抽取25%，得到稀疏的规则地震数据，基于完全相同的处理流程和速度模型进行叠前深度偏移，对比在同样采集密度下稀疏规则采样和随机采样得到的数据成像效果(图7、图8)。