1 引言

随着油气地震勘探目标的复杂程度日益提高，地震数据采集工作正朝着高覆盖次数、高密度、宽方位角的方向发展，致使地震数据采集系统的道容量不断的得到提升(Heath, 2012; Keho and Kelamis, 2012).据有关统计，地震数据采集系统的道容量自20世纪70年代初起，就随着勘探仪器的更新换代成指数增长：平均每3.5年，地震数据采集系统的带道能力就会提高一倍，这也被称为地震勘探仪器中的“摩尔定律”(Monk, 2006；Keho and Kelamis, 2012).目前，已有地球物理服务公司推出了百万道级的地震数据采集系统，并在实际工作中进行了应用(Lansley, 2013；Ellis，2014).在此趋势下，本文提出了一种基于压缩感知的地震数据采集技术与方法，通过软件上的改进，以引入一定范围内的噪声为代价，实现提高地震数据采集系统的道容量的目的，从而为突破地震数据采集系统大规模实时数据传输的瓶颈提供了一种新的途径.

压缩感知理论自2006年被提出来以后，在图像处理、核磁共振成像、遥感成像、传感网络、地震数据处理等领域得到了广泛的关注(Donoho, 2006; Baraniuk, 2011; Li et al., 2013; Herrmann et al., 2012；Zhang et al., 2015).压缩感知打破了传统的Nyquist-Shannon采样理论的限制，为更高效率地信号采集提供了理论依据(Candes et al., 2006; Donoho, 2006).在地震数据采集方面，Herrmann(2010)将压缩感知理论运用在地震勘探工作中，提出了根据地震数据稀疏性的随机采样方法，以减少数据维度，避免数据的维度灾难，从而提高地震采集的效率；陈生昌等(2015)和王汉闯(2015)已经使用压缩感知理论，以随机采样的方式，对地震信号的高效采集进行了初步研究.在一般化的信号采集方面，刘小洋等(Liu et al., 2015)根据压缩感知理论，提出了应用于无线传感网络中压缩形式的数据回收方法；Ravelomanantsoa等(2015)将压缩感知的思想应用于无线体域网络中，以实现低功耗、高效率的多种信号采集任务.

在地震勘探领域，压缩感知理论最初只是应用于地震数据去噪、重建等方面.自2014年以后，马坚伟、陈生昌等逐渐尝试通过随机采样的方式进行高效率地震数据采集工作的相关研究.近几年，虽然已有学者提出利用压缩感知理论，以数据压缩传感的方式完成无线传感网络、无线体域网络中的信号采集工作，但是地震勘探中，结合数据传输的网络拓扑结构使用压缩感知理论进行地震数据压缩传感的相关技术与方法还没有被提出过.因此，本文借鉴无线传感网络中应用压缩感知理论的采集框架，提出了基于压缩感知的多跳地震数据采集技术与方法，以压缩传感的方式，降低系统中所需要传输的数据量，从而提高整个系统的道容量.此外，通过在吉林大学研制的无缆自定位地震勘探系统中进行软件上的修改，完成了该方法的应用.

2 压缩感知原理

在一个有N个节点的传感网络中，若将编号为j的节点在一段时间内所采集到的数据表示为一维向量x_j(j=1，2，…，N)，则采集到的数据则为一个二维矩阵X=[x₁，x₂，…，x_N]^T.传统的信号采集方式是将原始信号X进行记录，而后对X进行非线性数据压缩，以便于传输或存储.压缩感知理论则提出通过随机化次采样技术将信号采样和压缩编码合并为一个操作，从而避免了高分辨率的信号采集过程和大量的非线性数据压缩计算，同时随机化的信号采样过程有利于压制数据采集中的随机噪声(Herrmann, 2010).压缩感知的前提条件是被采样信号具有较好的稀疏性.通过利用这种稀疏性，压缩感知能够实现低采样率下的高效率信号采集.假设传感网络中每个节点在第i时刻采集到的信号组成x_i(x_i∈R^N)，且在Ψ变换域上的具有K-稀疏(KN)的性质(‖θ_i‖=K)，即信号x_i在Ψ变换域上的系数θ_i(x_i=Ψθ_i)最多只含有K个非零值或主要的元素，其余的(N-K)个元素均可以在信号没有明显失真的前提下被舍去，那么该信号可以通过压缩感知的方式进行高效采集.在基于压缩感知的数据采集过程中，原始采集到的信号x_i乘以观测矩阵Φ得到测量值y_i(y_i∈R^M)，其中Φ∈R^M×N(N>M>2K)，从而使原被测信号的维度从N降为M，进而实现压缩的效果(Candes and Wakin, 2008; Leinonen et al., 2015; Liu et al., 2015).该信号采集过程可表示为：

(1)

当完成信号采集后，通过对测量值进行恢复操作即可获得原始信号.该过程具体形式为：

(2)

(3)

其中，为恢复后的被测信号在Ψ变换域上的系数，为恢复后的被测信号.

为了最终在解决压缩感知中的逆问题时能够获得唯一解，从而恢复出原始信号，观测矩阵和被测信号的稀疏变换矩阵要满足有限等距限制(Restricted Isometry Property)条件(Candes and Tao, 2005).观测矩阵可以在稀疏变换矩阵确定后，由RIP条件求得.常用的观测矩阵可以分为：随机性观测矩阵和确定性观测矩阵两类(Candes and Wakin, 2008; Chartrand and Staneva, 2008; Duarte-Carvajalino and Sapiro, 2009).典型的随机性观测矩阵有高斯随机观测矩阵和伯努利随机观测矩阵两种；典型确定性观测矩阵有基于混沌理论构造的循环矩阵、基于Legendre序列的双极性观测矩阵等(Zhang et al., 2016；Zhao et al., 2016).稀疏变换矩阵可以分为固定变换(fixed transform)和学习字典(dictionary learning)两类(Candes and Wakin, 2008; Rubinstein et al., 2010；Ravishankar and Bresler, 2013).应用于压缩感知的固定变换主要有curvelet变换、ripplet变换、dreamlet变换等；学习字典一般是通过机器学习等算法根据先验数据计算得到的自适应冗余字典(Geng et al., 2009; Rubinstein et al., 2010；Ghahremani and Ghassemian, 2015).在原始信号恢复过程中，公式(2)可以通过贪婪算法(直接利用L₀范数)、凸松弛优化(使用L₁范数代替L₀范数)等方法进行求解(Figueiredo et al., 2007；Donoho et al., 2012).

3 基于压缩感知的地震数据采集框架

目前，地震数据采集系统的每条测线中，无论是通过有线通信技术还是无线通信技术完成地震数据传输，主要使用的都是如图 1所示的多跳接力式架构(张林行，2007；Crice, 2014；张晓普，2016).图 1中d_j(j=1，2，…，N)代表着第j个采集节点在地震信号采集过程中采集到的地震数据.采用多跳接力式架构组成的地震数据传输网络存在着严重的通信带宽瓶颈问题，即越靠近根节点处需要传输的数据量越大，从而在实时数据回收的系统中需要的通信速率越高.因此，常规数据采集模式下地震数据采集系统中的道容量与系统中通信带宽的大小成正比，表达式为：

(4)

其中，N为系统的道容量；w为系统的通信带宽.

针对常规数据采集模式下的通信带宽瓶颈问题，本文提出了如图 2所示的基于压缩感知的地震数据采集框架，其中的ϕ_ij(i=1，2，…，M; j=1，2，…，N)为观测矩阵Φ中的元素.该数据采集方法通过在不同观测系数(权重值)下的传递采集系统整体的测量值代替每个采集站独自的测量值，从而使系统中传输的数据量减小，增加系统的道容量，同时消除了数据回收过程中的瓶颈效应，使整个系统的通信负载达到均衡.该方法的具体实现思路如下：第一个节点先将其采集到的数据乘以ϕ_i1，然后把得到的乘积发送给第二个节点；之后，第二个节点将其采集到的数据乘以ϕ_i2，再将ϕ_i1d₁+ϕ_i2d₂的值发送给第三个节点.以此类推，每个节点将自己采集到的数据乘以观测矩阵Φ中对应的列，然后将这个“加权”了的测量值与上一个节点发来的测量值作和，并将得到的新的测量值发送给下一个节点.最终，根节点接收到测量值是所有节点采集到的数据的加权和，可以表示为.在该过程中每个节点都使用了M个不同的观测系数(权重)，所以根节点最终能够收到M个测量值.根据压缩感知理论，当采集到的信号是可压缩(即在某个域上是稀疏信号)，并且信号的稀疏域与观测矩阵满足RIP条件时，根节点能够在M＜N的条件下恢复出原始采集到的数据d_j.因此，采用基于压缩感知的数据采集模式，系统中各节点间的数据传输量为一个恒定的常数，达到了数据传输的负载平衡，消除了系统中通信带宽瓶颈的问题.此外，更重要的是与常规采集模式下相比，系统中节点间的最大数据传输量由N降为了M，即系统的道容量最多能够提高到N/M倍.

在本文提出的基于压缩感知的地震数据采集框架中，观测矩阵、恢复算法及稀疏变换矩阵的设计如下：

随机观测矩阵已经被证明具有许多优越性.例如，随机矩阵作为观测矩阵时，可以在较高概率的条件下满足RIP条件；当观测矩阵中元素的分布具有零均值和有限方差时，随着M的增大，相关性会得到收敛；随机结构使观测值是一般化的，所以随机化的观测矩阵对于测量值丢失等情况还具有鲁棒性(Baraniuk et al., 2008；Candes and Wakin, 2008).但是由于完全的随机矩阵在硬件实现时是困难的，所以在目前常见的无线传感网络中，没有采用完全随机矩阵作为压缩感知中的观测矩阵的应用.然而，在地震数据采集网络中，每个采集节点往往还具有着较强的计算功能.例如：吉林大学自主研制的无缆自定位地震仪使用ARM9处理器，而且应用了嵌入式Linux操作系统，拥有较为丰富的软硬件资源，为实现生成随机数提供了便利的平台，可以通过软件编程对完全随机矩阵进行实现(杨泓渊，2009).而且，在实际地震数据采集工作中，各节点独自工作，易于生成独立同分布的随机矩阵.因此，在本文提出的基于压缩感知的数据采集框架中，选择高斯随机矩阵作为观测矩阵，观测矩阵中的所有元素都是独立同分布的，且服从N(0, 1/M).

用于压缩感知的恢复算法有很多种，主要可以分为两大类：贪婪迭代算法和凸优化算法(李珅等，2013).在计算精度要求不高的前提下，贪婪迭代算法的复杂度较低，计算速度快；凸优化算法的精度比较高，但是计算速度比较慢.由于地震数据采集系统需要对地震数据进行实时回收以便进行数据质量监控，因此应该选择具有复杂度较低的恢复算法，以便确保系统的实时性，从而提高采集工作的效率.所以，本文在复杂度较低的贪婪迭代算法中进行选择，将其中最具代表性的正交匹配追踪算法(Orthogonal Matching Pursuit)作为基于压缩感知的地震数据采集中的恢复算法(Tropp and Gilbert, 2007).

如何获得地震数据的稀疏变换矩阵是基于压缩感知的地震数据采集中最关键的问题.因为稀疏变换矩阵直接影响着地震数据能够到达的稀疏水平K，而根据压缩感知理论，M的取值至少应大于2K(Candes et al., 2006).所以，稀疏变换矩阵的选择直接影响着该方法能够提高道容量的倍数.此外，也有相关文献提出稀疏变换后的系数应该满足幅值衰减条件才能获得较好的压缩效果(Cohen et al., 2009).该幅值衰减条件可表示为：

(5)

其中，|θ|_(i)为地震数据稀疏变换后系数中绝对值第i大的值，C为常数，s为地震数据的压缩性(s越大，数据可压缩性越大)，一般认为s应大于等于0.5.

因此，为了实现较好的压缩效果，必须得到地震数据的稀疏域.目前，常用的稀疏变换方法主要可以分为两类：固定变换和学习字典.固定变换的方法具有易于实现、一般情况下计算速度快的特点，但是由于没有充分利用地震数据间的特点以及一些先验知识，所以无法针对具体数据实现其最佳的稀疏表示结果；字典学习的方法具有自适应能力强的特点，能够根据数据自身的特点，通过机器学习算法发掘出适合数据本身的稀疏表示方法，从而达到较好的稀疏结果，但是该类方法的计算量通常较大，不适于实时地震数据采集的工作情景.针对以上情况，本文提出了基于有序并行原子更新的字典学习算法，通过将一小部分地震数据作为先验知识(初始字典)，对其进行复杂度较低的处理，从而得到适合该数据的稀疏表示方法.

4 基于有序并行原子更新的字典学习算法

常用的经典字典学习算法主要有：广义主成分分析(Generalized Principal Component Analysis)，最优方向法(Method of Optimal Directions)和K奇异值分解算法(K-SVD Algorithm)三种(Engan et al., 2000; Vidal et al., 2005; Aharon et al., 2006；Rubinstein et al., 2010).GPCA算法主要是针对解决数据聚类和降维问题而提出的，其独特具有的一个性质是它能在给定训练集合中探测到需要的字典维度.但是，该算法具有指数计算复杂度，因此在子空间的数量和维度增加时，算法的开销变得十分大，不适于处理地震数据.MOD算法是用于解决稀疏表示问题的最早算法之一.该算法仅仅要求几次迭代过程就能够收敛，而且总体上看是一个非常有效的算法.然而，该方法因需要求矩阵的逆而有着较高的复杂度.K-SVD算法主要目的与MOD算法一致.该算法最显著的特点是其通过每个原子逐个进行简单高效的更新从而实现整个字典的更新，该过程无需求矩阵的逆.此外，在更新字典中各原子的过程中，训练集合的稀疏系统也同时进行着更新，因此该算法与MOD算法相比具有较快的计算速度.因此本文主要基于K-SVD算法的框架对地震数据的稀疏算法进行设计.

K-SVD算法对输入的初始字典中每个原子进行一个接一个的优化，从而得到适合于训练集合的稀疏变换矩阵.在每个原子更新过程中，该原子对应的所有非零系数也得到了更新.换句话说，只有与该原子相关的信号参与了该原子的更新过程.设{t_i}_{i = 1}^L为一个由L个训练信号组成的P维数据集合，将所有训练信号作为列向量放入矩阵T中，T∈R^P×L.设A为稀疏变换矩阵(A∈R^P×K)，X为T在A下对应的稀疏系数(Y≈AX).那么需要更新的就是A中第i个原子a_i，和其对应X中第i行稀疏系数x_i^T的非零系数(i=1，2，…，L).定义ω_i={j:x_i^T(j) ≠0}，Ω_i为一个大小为L×|ω_i|的矩阵，只有第(ω_i (j), j)项元素为1，其余元素均为0.那么，更新的过程就可以描述为解决如下最小化问题，公式为：

(6)

其中，E′_i=E_i(:, ω_i)，而E_i=Y-=x_i^TΩ_i.该最小化问题实际上可通过对E′_i进行奇异值分解从而得到近似解.

在K-SVD算法中，每一次原子更新的过程中都需要进行一次奇异值分解.然而，在对高维矩阵进行奇异值分解时，通常都是基于Jacobi旋转或Hestenes-Jacbobi等算法进行实现，二者算法复杂度均为O(n³)(Serra et al., 2017).因此，K-SVD算法整体复杂度约为O(n⁴)，算法所需的计算时间较长.为了进一步加快K-SVD算法的计算速度，Sadeghi等提出了PAU-DL(Parallel Atom-Updating dictionary learning)算法，并实现了加速字典收敛速率和提高变换后系数稀疏水平的效果(Sadeghi et al., 2013, 2014).该算法采用和K-SVD算法一致的思想和计算框架，但是没有通过奇异值分解的方法对E′_i进行求解，而是利用迭代的方法，交替求得a_i和的近似解，即， =a_i^TE′_i.该近似求解方法的理论依据和收敛性在相关文献中得到了证明(Rubinstein et al., 2008).然而，在PAU-DL算法中，迭代求近似解时并没有通过多次迭代逐一求出各原子最终的优化值，而是采用并行的思想，先将所有的原子统一更新后再进入下一次的迭代过程，从而提高算法的收敛速率.在基于3.8 GHz CPU、8G RAM的硬件条件，运行在Windows 7中的MATLAB R2010b环境下，PAU-DL算法能够在达到与K-SVD算法相同效果的前提下，将运行时间缩短80%(Sadeghi et al., 2014).

本文在PAU-DL算法的基础上，对原子更新的次序进行优化，提出了基于有序并行原子更新的字典学习算法.该算法不再按照任意次序，而是根据每个原子a_i所对应稀疏系数x_i^T的零范数进行排序，先对零范数较大的x_i^T所对应的原子a_i进行更新.这样设计的原因主要有两点：其一，可以将初始字典中与训练集合间的残差尽可能的优先分配到被频繁使用到的原子中，有助于将样本空间中的能量集中于少数几个代表主要成分的分量，从而加快字典的收敛速率；其二，更新过程中积累下来的残差应分解到某些比较少用到的原子中，这样有利于字典的表示系数达到较高的稀疏水平.基于有序并行原子更新的字典学习算法的实现步骤如下：

步骤1：输入初始字典A，训练集合T；

步骤2：根据OMP算法求得训练集合T的对应于字典A的稀疏系数X；

步骤3：计算字典A的残差E=T-AX；

步骤4：将字典A中各原子a_i按照其对应稀疏系数x_i^T的零范数进行排序，x_i^T的零范数越大，a_i的位置越靠前(即i值越小)；

步骤5：由E_i=E+a_ix_i^T，a_i =E_iΩ_i， =a_i^TE_iΩ_i，E=E_i-a_ix_i^T四个式子依次按照顺序对A中每个原子a_i及其对应的、E进行更新；

步骤6：重复步骤4至步骤5直至字典A达到收敛条件.

5 仿真数据及实测数据分析

为了验证基于压缩感知的多跳地震数据采集方法的可行性与有效性，本节分别对一组仿真数据和一组辽宁兴城地区的实际地震数据进行实验.在这两组实验中，选择高斯随机矩阵作为观测矩阵，OMP算法作为恢复算法.此外，将每组实验的前一部分地震数据做为训练集合，并把训练集合等分为两大小相同的两个子集合，然后对两个子集合进行奇异值分解，将分解后得到的两组特征向量合并成一个矩阵，作为基于有序并行原子更新的字典学习算法的原始字典.本文用于衡量基于压缩感知的地震数据采集质量的指标为信噪比(SNR)，即：

(7)

式中，f₀为原始地震数据；f为经过采集后恢复出的地震数据.

5.1 仿真地震数据实验

图 3为一个有限差分法模拟的炮记录，含有200道地震数据(即N=200)，采样率为1000 Hz，数据采集时间长度为8 s.该地震数据的前20%的数据被作为训练集合，即前1600个采样时刻采集到的地震数据.观测矩阵的行数取为列数的16%，即M=0.16N，换而言之，地震数据的压缩率设为16%.

为了充分验证基于有序并行原子更新的字典学习算法能够有效地得到地震数据的稀疏变换矩阵，本小节在含有较为丰富地下结构信息的第6400至第6700采样时刻的地震数据中，随机的选择某一个采样时刻的全部数据进行验证.将第6530采样时刻的地震数据作为测试数据进行验证.首先，将训练集合中前800个采样时刻的数据和后800个采样时刻的数据分别进行奇异值分解，得到两组秩为200的正交矩阵；然后将这两组正交矩阵合并，作为原子数量为400的初始字典.图 4a为初始字典的系数及其幅值衰减条件对比图，图中红线为第6530采样时刻的地震数据在初始字典下的变换系数(已按照其绝对值的大小进行排序)；蓝线为初始字典系数所对应的幅值衰减条件(见公式(5)，其中s=0.5)，用来验证其是否满足应用压缩感知理论所要求的稀疏性.图 4b为字典学习稀疏域系数及其幅值衰减条件对比图，其中红线为第6530采样时刻的地震数据在通过基于有序并行原子更新的字典学习算法后得到的稀疏变换域下的稀疏系数；蓝线为该稀疏变换域系数所对应的幅值衰减条件(其中s=0.5).由图 4a可以看出，使用该初始字典对地震数据进行稀疏时，变换后的能量没有足够集中，从而使得其变换后的系数在数值较大的部分衰减速度小于幅值衰减条件(在图 4a中绝对值大于5的部分，红线高于蓝线)，没有严格满足稀疏性的要求.由图 4b可以看出，即使蕴含较丰富信息的第6530采样时刻的地震数据，在通过字典学习后得到的稀疏域下也具有良好的稀疏性，且其系数的衰减速度大于幅值衰减条件.由图 4a和图 4b对比也可以看出，基于有序并行原子更新的字典学习算法能够有效的发掘出地震数据的稀疏域，可以作为基于压缩感知的地震数据采集框架中自适应地获取稀疏变换矩阵的方法.此外，基于有序并行原子更新的字典学习算法具有较低的计算复杂度，其计算复杂度不仅低于其他常用的字典学习算法，而且在实际的大数据量计算中也明显比curvelet变换、seislet变换快，有助于实现实时恢复出采集到的地震数据，进而完成实时的大道容量的质量监控.

在得到合适的稀疏变换矩阵后，通过OMP算法能够恢复出原始的地震数据.图 5a为经压缩感知采集后恢复出的地震数据，图 5b为恢复地震数据与原始地震数据的误差图.由图 5a和图 5b对比可以发现，地质构造基本上都得到了保留，剩余的误差不明显.用该方法进行数据采集的信噪比为18.8 dB，原始数据的均方根‖f₀‖₂的值为146.6326，误差的均方根‖f-f₀‖₂的值为16.9232.

在基于压缩感知的数据采集模式下的数据传输量占传统模式下数据传输量的百分比R_cs可由下式计算得到：

(8)

其中T_t为训练集合原始地震数据量的百分比.在本次实验中，经计算R_cs =32.8%，即在此情况下系统的道容量可以提高到原来的3倍多.

5.2 实测地震数据实验

为进一步验证该方法，本小节又随机选择了一组辽宁兴城地区的实际地震数据进行实验.该实际地震数据含有150道的地震数据(即N=150)，采样率为1000 Hz，数据采集时间长度为3 s.实验过程中，将该组地震数据的前20%的数据作为训练集合，即前600个采样时刻采集到的地震数据.观测矩阵的行数取为列数的16%，即M=0.16N，换而言之，地震数据的压缩率设为16%.

图 6a为该组在兴城地区实测得到的原始地震数据，图 6b为数据经压缩感知采集后恢复出的地震数据.由图 6a和图 6b对比可以发现，主要地震信息都能够得到采集.该原始地震数据的均方根‖f₀‖₂的值为210180，误差的均方根‖f-f₀‖₂的值为38244，数据经压缩感知采集后恢复出的地震数据与原始数据的信噪比为14.8 dB.由公式(8)可计算得整体数据的压缩比为R_cs =32.8%，道容量可提高到原来的3倍多.

6 结论

本文根据压缩感知理论，结合地震数据传输网络的拓扑结构，首先提出了基于压缩感知的多跳地震数据采集框架，以压缩传感的方式，减少系统中的数据传输量，从而在不做硬件修改的前提下提高系统的带道能力.为了能够充分的减少传输的数据量，同时保证数据采集的信噪比，本文还提出了通过自适应方式获得地震数据的稀疏变换矩阵的方法——基于有序并行原子更新的字典学习算法.经过模拟实验，证明该算法能够在快速的前提下，得到满足幅值衰减条件的稀疏变换矩阵.基于压缩感知的多跳地震数据采集方法也通过仿真数据和实测数据进行了验证，信噪比控制在14 dB以上时，即引入噪声的幅值在原始数据幅值的18%左右时，可以最多将系统的带道能力提高3倍以上.

采用本文提出的基于压缩感知的多跳地震数据采集框架还可在观测矩阵、恢复算法以及字典学习算法三方面进行改进.其中，字典学习算法为最主要的影响因素.本文提出的基于有序并行原子更新的字典学习算法虽然能够获得满足幅值衰减条件的稀疏变换矩阵，但是没有充分利用地震数据的物理意义(例如，同相轴的斜率、振幅的衰减等).因此，下一步研究可以设计结合地震波场具体物理意义的字典学习算法，获得更稀疏的表示方法，从而更大幅度地提高系统的道容量.

致谢

感谢国土资源部地球探测技术及仪器重点实验室对本文提供的数据支持，感谢国土资源部地质环境监测技术重点实验室张青教授对本文的指导和建议.