随着我国油气资源勘探向中西部发展, 地球物理的研究也逐渐转向地表构造复杂的区域。在这些区域进行勘探, 不论是资料采集, 还是后期的资料处理和解释, 都对地球物理人员提出了更高的要求。

对于地震勘探来说, 在复杂地区施工, 往往会遇到峭壁、河流、沟壑、村庄和工业区等不易甚至无法布置检波器的情况, 从而给数据的采集带来了一定的困难, 甚至导致无法采集的情况发生。如果因为这些原因导致原始数据的大量缺失, 就会严重影响整个工区的勘探质量和成像效果, 从而给资料解释带来困难。

经典信号采集理论一般都会遵循Shannon-Nyquist采样定理, 即:若要使采集到的数据能不失真地保持原信号中的信息, 采样频率必须是原信号频带宽度的两倍以上。地球物理的野外资料采集也都遵循该定理。为了描述目标体的精细特征, 往往需要采集更精确的信号, 因此, 需要通过在空间上加密采样点数和在时间上减小采样间隔来实现。这样的采集方式大大增加了采集的数据量, 从而导致采集成本的急剧增高。

基于信号稀疏性的采样理论^[1-2]的压缩感知(Compressed Sensing, CS)在图像压缩、无线通信、模式识别和医疗成像等研究领域受到了格外的关注。2006年, CANDES等^[3]提出基于压缩感知的信号恢复技术, 大大推动了压缩感知在地球物理勘探方面的应用。

压缩感知技术的发展和应用基于大部分信号都可以进行稀疏表示。在稀疏域对信号进行重建, 可以获得原始的真实信号。压缩感知技术突破了Shannon-Nyquist采样定理的局限, 大大节省了采集数据所占的空间, 成为提高地震勘探效率的重要方法。

目前, 压缩感知理论在地球物理勘探领域得到广泛应用^[4-15]。基于压缩感知理论, HERRMANN等^[4]提出了一种多震源的采集框架; MOLDOVEANU^[5]探索了海上数据采集的随机观测方法; MOSHER等^[6-7]提出并完善了一种在观测系统设计中基于约束条件的更优选择激发点和接收点位置的非均匀采样方法; 陈生昌等^[8]提出了一种针对地球物理数据的高效采集方法, 并通过归一化各列数据的感知矩阵的互相关系数最小作为约束条件来设计最优的采集方案。

基于压缩感知的地震采集技术, 获得的数据需要进行重建^[9-10]。HENNENFENT等^[11]提出了利用抖动欠采样进行波场重建的方法, 并基于曲波变换进行了地震数据的重建; WU等^[12]提出了Dreamlet域地震成像和数据重建的方法; LI等^[13]提出了基于插值压缩感知的数据重建方法(interpolated Compressive Sensing, interpolated CS), 设计了非均匀优化采样观测系统(Non-Uniform Optimal Sampling, NUOS), 并通过地震数据重建, 分别在海上和陆上取得了良好的试验效果; 王华忠等^[14]详细讨论了压缩感知技术在地震勘探应用中的问题和难点, 从不同方面展示了压缩感知的有效应用。本文通过降低感知矩阵的最大互相关值, 从而进行非规则观测系统的优化设计; 并利用L₀正则化和L₁正则化混合迭代的方法进行欠定采样矩阵的求解, 完成了地震数据的重建。

1 基于压缩感知的非均匀观测系统设计方法

压缩感知理论的基本假设是目标信号具有稀疏性或者可压缩性, 也就是说, 目标信号或者其在某个变换域中只有有限个成分不等于0(对应稀疏性), 或者只有有限个成分远大于0(对应可压缩性)。若目标信号有K个成分不等于0, 则称该信号为K稀疏的。

假设具有稀疏性或者可压缩性的目标信号为x, 稀疏变换采用傅里叶变换F, 则正交稀疏基为傅里叶基函数或分窗傅里叶基函数, 于是有:

$ \mathit{\boldsymbol{x}} = {\mathit{\boldsymbol{F}}^{\rm{H}}}\mathit{\boldsymbol{s}} $

(1)

式中:s为信号x在傅里叶域的稀疏表示; (·)^H表示共轭转置。

这样采集得到的数据可以看作是采样函数或采样矩阵与目标信号相乘的结果。记采样矩阵为Φ, 当进行满采样时, 有Φ=I, 其中, I表示单位矩阵。当不能满足满采样条件(如采样数不足)时, Φ为I抽出的若干列组成的矩阵, 只有采样位置对应的列向量被保留, 则有采样数据:

$ \mathit{\boldsymbol{y}} = \mathit{\boldsymbol{ \boldsymbol{\varPhi} x}} = \mathit{\boldsymbol{ \boldsymbol{\varPhi} }}{\mathit{\boldsymbol{F}}^{\rm{H}}}\mathit{\boldsymbol{s}} = \mathit{\boldsymbol{ \boldsymbol{\varPsi} s}} $

(2)

式中:Ψ=ΦF^H记为感知矩阵。

根据压缩感知理论, 成功重建K稀疏的目标信号的充分必要条件是感知矩阵Ψ满足有限等距准则(Restricted Isometry Property, RIP)条件, 即对于任意2K稀疏向量v均有(1-δ)‖v‖₂²≤‖Ψ·v‖₂²≤(1+δ)‖v‖₂²成立。然而, 求解满足RIP条件的矩阵是个NP-hard问题, 取而代之, 可采用一种可计算的准则来进行采样设计。设矩阵Ψ的列向量为ψ_i, 列向量间的最大互相关值为:

$ \mu \mathit{\boldsymbol{ = }}\mathop {\max }\limits_{i \ne j} \frac{{\left| {\mathit{\boldsymbol{\psi }}_i^{\rm{H}} \cdot {\mathit{\boldsymbol{\psi }}_j}} \right|}}{{\left( {{{\left\| {{\mathit{\boldsymbol{\psi }}_i}} \right\|}_2} \cdot {{\left\| {{\mathit{\boldsymbol{\psi }}_j}} \right\|}_2}} \right)}} $

(3)

JAMAIL-RAD等^[16]指出, 感知矩阵列向量间的最大互相关值即为非规则采样归一化频谱的最大非零频率振幅。在频率域, 采样数据y的频谱是感知矩阵Ψ的频谱和目标信号的稀疏表达s的频谱之间卷积作用的结果。μ就是因为非规则采样造成傅里叶基(稀疏基)的正交性被破坏所引起的最大频谱泄漏, 压制μ使感知矩阵Ψ的频谱近似Delta函数, 否则, 感知矩阵Ψ的频谱出现多个峰, 与稀疏表达s卷积作用后, 将在频率域(稀疏域)产生“假频”噪声。CANDES等^[17]指出, μ越小, 则信号在非规则采样后能够重建的概率就越高。由于Ψ=ΦF^H, 其中F是固定的, 因此, 在采样数目不足时, 可以通过改变采样矩阵Φ, 即改变非规则采样点的分布来降低最大互相关值, 从而优化采样矩阵, 提高目标信号重建的概率, 即:

$ \mathit{\boldsymbol{ \boldsymbol{\hat \varPhi} }} = \arg \min \mu \left( \mathit{\boldsymbol{ \boldsymbol{\varPhi} }} \right) $

(4)

JAMAIL-RAD等^[16]给出了μ达到下界时的一种确定性采样方式, 然而这种采样方式只在一些特定的条件下才能成立。注意到该优化问题是非凸的, 具体优化时只要能够找到局部最优解就足够满足需求^[18-19]。这里我们采用贪心序贯策略来优化采样点的位置, 进而构建采样矩阵。具体策略如下:

1) 首先确定目标信号的网格点数目和间距, 将所有目标信号的网格点作为候选采样点;

2) 遍历所有候选采样点, 计算将其加入采样矩阵后的μ值, 取使μ值最小的候选点作为新的采样点, 更新采样矩阵;

3) 重复步骤2)的操作, 直到采样点数达到规定数目;

4) 在初步选定规定数目的采样点位置后, 利用时间抖动(Jitter)方法^[18-19], 以(4)式为目标, 进行采样微调, 最终确定采样矩阵。

以一维简单情况为例, 假设在751个网格点上, 选取375个点(近50 %)作为采样点。当规则采样时, 采样点分布及其归一化后的傅里叶谱如图 1所示。可以看出, 欠采样的规则采样函数会产生很强的频谱泄漏, 造成假频混叠。同样的网格和采样点数, 当随机采样时, 采样点分布及其归一化后的傅里叶谱如图 2所示。可以看到, 随机产生的非规则采样降低了采样矩阵与傅里叶矩阵的相关性, 没有假频产生, 但是在整个频率范围都出现了很低的频谱泄漏, 最大幅值为0.7949×10^-1。保持网格和采样点数目不变, 采用贪心序贯策略获得优化的非规则采样, 采样点分布及其归一化后的傅里叶谱如图 3所示。图 3中采样点非均匀分布, 压制了假频, 降低了整个频率范围内的频谱泄漏, 频谱泄漏的最大幅值0.5382×10^-1, 比随机采样的频谱泄漏小。

由压缩感知原理可知, 更小的频谱泄漏表明由非规则采集而造成的漫布在整个频带内的杂乱噪声更弱, 这样采样频谱就更加逼近规则满采样时的频谱, 更有利于目标信号的重构, 优化后的非规则采样比随机采样所得到的数据能够完美重建的概率更高。当采样点数较少或者采样率较低时, 这种差异更加明显。而当采样点较多时, 采样条件往往会受其它因素如地形、实际施工条件限制等影响, 此时可在采样点优化设计过程中加入这些约束, 从而得到满足施工要求的采集方案, 这是随机采样很难实现的。另外, 对于某些二次施工工区, 采样设计时还需要考虑其它指导目标如地下照明、数值模拟恢复信号的信噪比等, 此时可在优化目标(4)式中加入其它约束项, 然后优化完成满足要求的非规则采样点设计, 这也是随机采样很难实现的。

地震信号的稀疏性已经被广泛应用。对于三维地震勘探的数据采集情况来说, 假设真实信号为D_real, 其时间维度为N_t, 与非规则采集的时间维度相同, 而空间维度为N=N_x×N_y, 与非规则采集数据重建后的规则数据所具有的空间维度对应。选择傅里叶变换作为稀疏变换, 先将均匀满采样的时间方向做傅里叶变换, 真实信号D_real由t-xy域变换至f-xy域。由于采样网格只与空间有关, 故对所有频率切片使用相同的非规则采样设计。为了降低公式的维度, 不失一般性, 取其第q个频率切片, 记为S^q, 再对S^q做空间方向上的傅里叶变换, 将其变换至f-k域, 记为 ${\mathit{\boldsymbol{\tilde S}}^q}$ , 则有:

$ {\mathit{\boldsymbol{S}}^q} = \mathit{\boldsymbol{F}}_x^{\rm{H}} \cdot {{\mathit{\boldsymbol{\tilde S}}}^q} \cdot \mathit{\boldsymbol{F}}_y^{\rm{H}} $

(5)

式中:F_x和F_y分别为x方向和y方向上的傅里叶变换矩阵。

引入向量化操作^[16], 将二维目标信号转化为一维信号。向量化算子:

$ {\rm{vec}}\left( \mathit{\boldsymbol{A}} \right) = \mathit{\boldsymbol{a}} $

(6)

表示将任意矩阵A各列叠加排列成一个很长的列向量a。向量化算子具有如下性质:

$ {\rm{vec}}\left( {\mathit{\boldsymbol{A}} \cdot \mathit{\boldsymbol{B}} \cdot \mathit{\boldsymbol{C}}} \right) = \left( {{\mathit{\boldsymbol{C}}^{\rm{T}}} \otimes \mathit{\boldsymbol{A}}} \right) \cdot {\rm{vec}}\left( \mathit{\boldsymbol{B}} \right) $

(7)

式中:A, B, C是任意3个连乘有意义的矩阵; ⊗表示克罗内克积。通过以上计算, 可以将f-k域中第q个频率切片的数据向量化, 得到真实信号为D_real向量化的稀疏表示 ${\mathit{\boldsymbol{\tilde s}}^q}$ 。

在实际非规则采样时, 时间上均匀满采样, 与真实信号的时间维度保持一致, 均为N_t。在空间上, 使用基于压缩感知的非规则优化采集方案, 实际的采样点数为M, M＜N。也就是说, 要应用一个采样算子G, 对稀疏表示 ${\mathit{\boldsymbol{\tilde s}}^q}$ 进行采样, 采样算子G实际上就是从N维单位矩阵I_N中取M行而得到, 其决定了实际非规则采集的采样点(检波点和炮点)位置, 且G的选择与这些非规则位置一一对应。于是有:

$ \mathit{\boldsymbol{G}} \cdot {\mathit{\boldsymbol{s}}^q} = \mathit{\boldsymbol{G}} \cdot \left( {\mathit{\boldsymbol{F}}_y^{\rm{H}} \otimes \mathit{\boldsymbol{F}}_x^{\rm{H}}} \right) \cdot {{\mathit{\boldsymbol{\tilde s}}}^q} $

(8)

目标就是根据估计 ${\mathit{\boldsymbol{\tilde s}}^q}$ 重建D_real。通过向量化, 将二维数据转化为一维数据求解, 则此时感知矩阵Ψ可表示为:

$ \mathit{\boldsymbol{ \boldsymbol{\varPsi} }} = \mathit{\boldsymbol{G}} \cdot \left( {\mathit{\boldsymbol{F}}_y^{\rm{H}} \otimes \mathit{\boldsymbol{F}}_x^{\rm{H}}} \right) $

(9)

于是, 可以通过最小化最大互相关值来确定采样算子G, 即:

$ \mathit{\boldsymbol{\hat G}} = \arg \min \mu \left( \mathit{\boldsymbol{G}} \right) $

(10)

采用贪心序贯策略解决该问题时, 需要增加候选网格的维度。考虑到非规则采样在滚动排列上的复杂性, 实际施工的试验区采用了非滚动全检波点接收的方案, 降低了炮点和检波点的联系。这里, 将炮点和检波点的分布进行独立设计, 减少了计算量。同时, 假设没有地下构造的先验信息, 仅凭目标信号的采样网格精度和采样数目来约束设计方案。

在实际观测系统设计优化时, 需要根据实际情况添加一定的约束条件。特别是在设计检波点位置时, 考虑到设备的实际情况, 需要添加约束条件。第一是检波线距非规则, 检波点非规则地分布在这些检波线上, 而不是所有检波点都完全分散分布; 第二是相邻检波器之间的最大距离受实际施工时电缆线长度的限制。考虑到这些因素, 在设计检波点位置时, 先优化检波线的位置, 以(4)式为优化目标; 然后固定检波线的位置, 将M个采样点位置平均分布在N_x×N_y的网格上, 然后利用Jitter采样^[17]的方式, 让采样点在初始位置附近抖动取值, 以(10)式为目标进行迭代优化。每次迭代只抖动其中一个采样点, 而其它点保持不变, 取使μ最小的位置为新的采样点, 更新该点位置, 然后选择下一个点进行抖动, 进行下一次迭代, 直到抖动任何一个采样点均不能使μ变小为止^[18]。为加速迭代收敛, 可以先让各采样点在平均分布位置附近随机抖动, 然后以此为初始条件进行迭代优化, 同样可以达到满足要求的优化设计结果。

对炮点进行设计时, 考虑最大炮间距和最小炮间距, 在炮间距限制范围之内进行优化设计, 使得非规则炮点的位置在目标网格上完全打散分布。此时可以根据序贯思想, 利用贪心算法, 逐个添加炮点位置以使得当前的μ最小。每步迭代先固定之前已有的所有采样点, 然后在剩余未被选中的网格点中, 选出使μ最小的网格点作为采样点。要加速运算, 可以控制变量, 每次只在一个方向取得最小值, 然后换方向接着搜索, 直到收敛。在所有M个采样点全部确定之后, 再以(10)式为目标, 进行全局微调, 最终得到炮点的优化非规则采样设计。

2 TFT工区的实际应用

我们利用TFT工区的一块区域进行实际工区的非规则观测系统设计, 所选工区范围为14km×8km, 检波点布设到整个工区范围, 炮点区域为10km×2km。规则检波点设计为32条接收线, 线距240m, 每条线456道, 道间距30m。非规则设计是在原始规则检波点设计的基础上, 保持大致的覆盖范围, 但是只取总采样点数的75%, 即线数32条不变, 平均每条线取342道, 平均道间距40m, 共10944道。重建的目标网格为42条接收线, 线距180m, 每条线911道, 道间距15m。优化后的检波点分布如图 4所示, 每个蓝色的十字型代表一个检波点的位置, 最小点距为15m, 最大点距为75m。同时, 线距也进行了非规则优化设计, 最小线距为150m, 最大线距为705m。这里由于检波线数较少, 为保证重建质量, 降低频谱泄漏, 采样点没有落在实际的整数网格上。检波点位置归一化后的频谱如图 5所示, 在垂直于检波线方向上由于检波线条数较少, 频谱泄漏较强, 最大值为0.0930, 但仍然要比规则采样导致的假频小很多。

实际工区的非规则炮点设计为将1800个炮点分散在10km×2km的范围内, 目标重建网格为炮点距15m, 炮排距90m。优化后的炮点分布如图 6所示, 每个红色的方框代表一个炮点位置, 其归一化后的频谱如图 7所示, 可以看到, 频谱泄漏很小, 最大值为0.0444。

最终得到的该工区优化的非规则观测系统如图 8所示。

3 地震记录的模拟与重建

对于地球物理问题来说, 基于压缩感知的数学重建模型可表示为:

$ \mathit{\boldsymbol{y}} = \mathit{\boldsymbol{ \boldsymbol{\varPhi} x}} $

(11)

式中:x∈R^N表示原始地震数据; Φ表示采样矩阵; y∈R^M表示观测地震数据。这里, M＜N, 即采样矩阵是一个欠定矩阵, 意味着需要从欠采样的地震数据y中重建出原始地震数据。但是, 如果当原始数据x在某个域C中的系数是稀疏的, 即s=Cx, 那么方程(11)可写为:

$ \mathit{\boldsymbol{y}} = \mathit{\boldsymbol{ \boldsymbol{\varPsi} s}} $

(12)

式中:Ψ表示感知矩阵, Ψ=ΦC^H。

对此, 可以通过求解反问题:

$ \min \left\| {\mathit{\boldsymbol{y}} - \mathit{\boldsymbol{ \boldsymbol{\varPsi} s}}} \right\|_2^2 $

(13)

来获取稀疏系数s, 进而通过x=C^Hs求取原始数据x。但是, 由于采样矩阵是一个欠定矩阵, 可知上述问题是一个典型的NP-hard问题。对此, 需要对上述问题增加约束条件, 本文利用L₀正则化和L₁正则化混合迭代的方法求解上述问题。首先通过软阈值算法求解L₁范数约束问题:

$ \min \left\| {\mathit{\boldsymbol{y}} - \mathit{\boldsymbol{ \boldsymbol{\varPsi} s}}} \right\|_2^2 + \lambda {\left\| \mathit{\boldsymbol{s}} \right\|_1} $

(14)

通过对(14)式的求解, 可在进行下一步算法前获得一个好的迭代初值, 并估算信号的稀疏水平。然后, 通过硬阈值算法求解L₀范数约束问题:

$ \min \left\| {\mathit{\boldsymbol{y}} - \mathit{\boldsymbol{ \boldsymbol{\varPsi} s}}} \right\|_2^2 + \lambda {\left\| \mathit{\boldsymbol{s}} \right\|_0} $

(15)

反复迭代, 直至获得最终结果。

为了验证该算法重建地震记录的有效性, 采用主频为30Hz的Ricker子波合成500道的单炮地震记录, 道间距10m, 每道3001个采样点, 采样间隔为2ms, 如图 9a所示。图 9b为采样50%地震道后的不完整地震记录, 图 9c为使用本文算法恢复的地震记录, 图 9d为本文算法恢复的地震记录与原始模型记录的误差。

从图 9c的重建地震记录与原始记录(图 9a)的对比以及两者的误差(图 9d)可以看出, 两者同相轴清晰, 连续性相同, 振幅一致性准确, 重建误差小。因此, 该算法具有良好的数据恢复效果。

图 10是该地震记录的频率-波数谱。从原始地震记录和重建地震记录的频率-波数谱对比可以看出, 重建数据的频率-波数谱与原始记录的频率-波数谱非常接近, 两者的误差很小(图 10d), 进一步表明了该算法良好的重建效果。以上实验均表明, 利用基于压缩感知理论的重建方法, 可以较好地恢复缺失的地震数据。

4 结论与建议

基于压缩感知理论, 本文进行了地震勘探非规则观测系统的优化设计和地震数据的模拟重建。通过降低感知矩阵的最大互相关值的方法优化设计采样点的位置, 并应用于实际工区, 该过程通过添加部分约束条件来使得最终的设计在实际可行的情况下达到最优解。利用L₀正则化和L₁正则化混合迭代的方法进行地震数据的重建, 取得了良好的效果。

随着勘探精度的提高, 野外地震资料的采集都要求能够达到“两宽一高”的标准, 这就给野外采集带来了相当大的工作量。基于压缩感知的非规则观测系统的地震数据采集和信号重建, 可以达到“相同野外采集工作量的条件下获得更高分辨率的地震资料”或者“获得相同分辨率的地震资料需要较少的野外工作量”的目的。

因此, 利用压缩感知技术进行地震数据的采集和处理, 可以大大降低野外资料采集的工作量, 对目前“两宽一高”要求的地震勘探具有重要意义。

致谢: 感谢中石化西北石油分公司在正演模型建立方面给予的大力支持, 感谢中石化休斯敦研发中心的罗明秋博士在目标模型的建立和观测系统的设计方面给出的一些好的建议！