地震数据中含有丰富的有效信号，然而时空域的信号表达方式无法直接描述地震数据的特征.为此，可以用某种数学变换抓取地震信号的主要特征，实现地震数据的稀疏表达.这有利于地震信号处理、波场分析和反演成像.事实上，地震数据可以在变换域中稀疏表达，常用的变换域有Fourier域、小波域、Curvelet域、Radon域、Gabor域等.广义上讲，可以通过选取一组正交基、或框架函数、甚至基函数字典，使得变换域中的地震数据具有稀疏性质.利用数据在变换域中的稀疏性，可以实现地震数据的压缩与重构插值等.Xu等(2005)在Fourier域反演非规则地震数据频谱，实现数据重构.Herrmann和Hennenfent(2008)、白兰淑等(2014)，利用地震数据在Curvelet域中的稀疏性，实现地震数据重建.

然而，Curvelet变换无法直观描述地震波的传播特征，因而没有明确的地球物理含义.本文中，我们把叠前地震数据抽象为不同出射角度(即出射慢度矢量)以及不同到达时的局部平面波的线性叠加，而接收到的地震记录可以看作局部平面波波前在空间离散的检波器处的采样.在立体层析理论(Billette和Lambaré，1998)中，地震记录被简化为炮点、检波点坐标，炮点、检波点处的射线参数以及射线的双程走时.上述参数构成了立体层析中的数据空间，但忽略了地震子波的波形信息.

若用Gabor框架作为基函数，则可以地震数据精确重构.Žáček(2003)用相干态变换(coherent-state transform)方法构造一组完备的Gabor框架函数，实现了地震数据的高斯波包分解.由于Gabor框架函数的冗余度高(Sondergaard，2007)，该方法无法实现数据的稀疏表达.Geng(2011)忽略了高斯波包中一些不重要的参数，压缩了重构观测波场的高斯波包泛函空间，从而减少框架函数的数目.李辉等(2014)提出了地震数据的特征高斯波包(CGP)分解方法，在Radon域选取出地震数据的时空和方向特征，大幅减少了高斯波包框架的数量，实现了Gabor域的稀疏化表达.

基于线性Radon变换(LRT)的局部平面波分解的方法是主要的地震信号的线性表达方法，以其明确的物理含义(可以描述局部平面波的传播方向)而被广泛应用.非线性的Radon变换，如抛物Radon变换(Sacchi和Ulrych，1995)，双曲Radon变换(Liu和Sacchi，2002)，多项式Radon变换(牛滨华等，2001)等，则主要用于多次波压制、去噪、数据重构等方面.不同于Fourier变换和正交小波变换，LRT的基函数是非完备正交的，因而导致LRT结果难以实现稀疏表达.王雄文和王华忠(2014)将常规的LRT转化为压缩感知问题，在压缩感知的理论框架下实现了高分辨率的LRT.在Radon域中，可以测量局部平面波在炮点处的入射射线参数 P_s 及检波点处的出射射线参数 P_r .若引入了局部平面波的方向信息，叠前地震记录 d = d(x_s，x_r，t)表现为立体地震数据，可描述为 d = d(x_s，p_s，x_r，p_r，t).

在天然地震学中，通常用检波器排列的某种组合来分离相干信号和噪声.最基本的方法就是波束合成方法(Beam Forming，Rost和Thomas，2002；Krüger等，1993)，它利用局部平面波波前到达不同检波器的时差，对检波器接收到的地震记录进行时移叠加，使局部平面波相干叠加而其他信号相互抵消，从而提取某种特定的波型(震相).

波束合成的思想也可应用于勘探地震学中，如高斯束偏移方法(Hill，1990；Hale，1992)，在炮集中对不同检波器中心进行加窗局部倾斜叠加(即波束合成)，利用炮集的 τ-p 谱进行偏移.在射线束偏移(Sun等，2001)和控制束偏移(Gao，2006)中，通过选择 τ-p 谱的能量，实现有选择性的波场传播与成像.

波束合成算法需要估计局部平面波的方向信息.本文中，我们考虑Radon谱约束下的局部平面波反演方法(胡江涛等，2013)：

式中，X(ω)为局部空间窗内的地震信号，X(ω)=[ x₁(ω)x₂(ω)… x_nx(ω)]^T，S(ω)为待求解的局部平面波，S(ω)= [s₁(ω)s₂(ω)… s_np(ω)]^T，A(ω)为局部平面波传播的方向矩阵，A_i，j=exp(iωp_j(x_i－x₀)). R 的元素为狄里克雷函数： r_i，j=，其中 E(i)是第i个平面波在Radon谱上的能量.

上述Radon谱约束下的局部平面波反演方法可以提高Radon谱的精度并压制泄漏的噪声.更进一步地，为了得到更稀疏的局部平面波解(王雄文和王华忠，2014)，可以把式(1)改写为压缩感知问题，在L0范数意义下求解如下反问题：

利用匹配追踪方法求解(2)式，可以得到局部空间窗内的稀疏平面波解.

通过传统的线性Radon变换或者高分辨率的平面波反演方法，可以得到局部平面波及其慢度矢量(或射线参数 p).此时，我们可以考虑对炮集的波束合成.

在共炮点集中，对中心检波点 x_r0 的波束合成可以写为

其中，d(x_r，t，x_s)为地震数据； Φ(x_s，x_r0)为中心检波点处的倾斜叠加孔径，它代表炮点坐标位于 x_s 处，且检波点坐标位于中心检波点坐标 x_r0 的一个局部空间窗 W_r 内地震数据集合(|x_r－ x _r0| <W_r)； p_r 为中心检波点处的出射射线参数(代表局部平面波在地表的出射慢度的水平分量)，可以通过传统的Radon变换或本文中的Radon谱约束下的局部平面波反演方法得到；上标BF代表波束合成(Beam Forming).

我们把在中心检波点处进行一次波束合成之后的地震数据记为 d^BF(x_r0，p_r，t，x_s). 显然，经过波束合成后，实现了出射方向为 p_r 的局部平面波的同相叠加，而其他信号相互干涉.

若将地震数据分选为共检波点道集，此时可以再次在中心炮点 x_s0 处进行波束合成：

其中，Φ(x_s0)为中心炮点处的倾斜叠加孔径，代表炮点坐标 x_s 位于中心炮点坐标 x_s0 的一个局部空间窗 W_s 内地震数据集合(|x_s－ x_s0| <W_s)； p_s 为中心炮点处的入射射线参数(代表局部平面波在地表的入射慢度的水平分量)；上标DBF代表两次波束合成(Double Beam Forming).

同公式(3)，经过中心炮点处的再次波束合成后，实现了入射方向为 p_s 的局部平面波的同相叠加.相应的，我们把两次波束合成之后的地震数据记为 d^DBF(x_r0，p_r，x_s0，p_s，t). 其中，(x_s0，x_r0)为波束合成的中点坐标，而(p_r，p_s)为局部平面波的慢度矢量.

将(3)式中的波束合成代入(4)式中，可得同时在中心炮点和检波点处的波束合成算法：

在三维时，上述公式中的坐标与射线参数均为矢量，表示为 x_s=(s_x，s_y)，x_r=(g_x，g_y)，p_s=(p_sx，p_sy)，p_r=(p_rx，p_ry).经过特征波场合成，五维地震数据 d = d(x_s，x_r，t)被投影到九维的立体数据空间 d = d(x_s，p_s，x_r，p_r，t)中.由于(5)式波束合成的地震数据可以描述局部平面波的传播方向特征，因此我们称(5)式合成的地震数据为特征波场.

从公式(5)可以看出，特征波场合成与局部平面波的入射与出射射线参数有关，其本质是利用局部平面波的入射和出射慢度计算道间时差(Δt= p ·Δ x)，并将波束合成孔径内的地震道按照道间时差进行时移叠加.在时移叠加的过程中，只有给定慢度矢量(局部平面波入射与出射射线参数)的局部平面波波前实现了同相叠加，而其他的信号相互抵消，从而实现不同震相(波型)的分离，即特征波场合成.

事实上，上述方法也可以推广到任意线性排列的观测系统中(如VSP或井间观测系统等)，只要满足局部的炮点/检波点排列具有线性特征.相应的，局部平面波的慢度矢量(射线参数)定义为地震波走时对地震排列(炮点或检波点坐标)的方向导数.

此外，在频率域特征波场合成公式(5)对应的相移矩阵为 A _i，j=exp(iω p _j·(x _i－ x ₀))，其中，射线参数向量 p =(p_s，p_r)，炮检坐标向量 x =(x_s，x_r). 此时，基于反演理论的特征波场合成方法可以在L2或L0范数意义下求解(形式与公式(1)和(2)相同)，无需预先估计射线参数再进行特征波场合成.

地震数据偏移成像的过程可以描述为地表波场的反传播加上成像条件.以运动学射线追踪为传播算子的Kirchhoff叠前深度偏移(PSDM)技术在横向变速剧烈情况下存在多路径和振幅焦散等问题，基于射线束的成像方法(刘少勇等，2012)在一定程度上解决了这些问题.高斯束偏移技术(Hill，1990；Hale，1992)已经成为工业界一个重要的叠前深度偏移成像工具.Sun和Schuster(2001)提出菲涅尔带控制的波路径偏移方法，将积分法偏移在成像空间中沿等时面的投影转化为沿胖射线路径的投影，在计算效率上有一定优势.Liu和Palacharla(2011)对高斯束偏移进行简化，只保留其运动学特征，推导出简单的解析表达式来控制射线束的宽度和振幅，理论上说该实现方式效率要优于高斯束偏移，成像精度基本有保障.快速射线束(Fast-beam)偏移(Gao等，2006)便是该思想的一种实现.特征波偏移成像利用射线束波传播算子来反传局部特征波场，既保留了射线类成像的灵活性，又保证了成像精度.

特征波场是地震波场在高维空间中的稀疏表达，且局部平面波的入射与出射慢度矢量已知.利用这个性质，射线束类的偏移方法可以实现从“画弧”向“搬家”的转变，进一步提高成像效率.与上节的特征波场合成方法结合，匹配合适的局部射线束传播 算子，可以实现特征波叠前深度偏移成像(characteristic wave imaging，简记为CWI). 更进一步的，“搬家”或一对一映射的偏移方法与基于角度道集的层析速度反演可以自然的融合成一体，共享同样的数据空间和成像空间的数据体.成像与层析速度反演一体化也是特征波成像的一大优势.

传统基于炮道集的射线束(包括高斯束)成像条件可以描述为

其中： x 表示成像点坐标，W^BF 表示由动力学射线追踪得到的振幅加权函数，W 表示成像孔径，即对地下成像点有贡献的地表炮检点范围.类似的，特征波成像条件可以表示为

其中，W^DBF 表示一个更为简单的振幅加权函数，它是成像空间坐标的函数，在不同的成像方法中可有不同的表达.不论在传统射线束成像还是特征波成像过程中，由于可以方便的根据走时梯度计算射线参数，进而通过射线参数求出张角(三维时是张角和方位角).张角 θ 可以由下式表示：

传统成像条件和特征波成像条件的几何解释可由图 1表示.图中红色反射层为地层真实位置，传统射线束成像方法是将局部平面波投影到整个椭圆(三维情况为椭球)上.由于特征波场在地表的入射与出射慢度已知，因而特征波成像只对特征波场进行“能量搬家”.红色小椭圆的位置代表一对特征波射线束在成像域的菲涅尔带范围.从图 1可以看出，基于特征波场表达，实现了成像方式从“画弧”到“搬家”的转变.传统积分类成像的画弧成像过程必然会带来画弧噪声，而搬家的方式可最大限度减少此类噪声同时能兼顾菲涅尔带的成像分辨率，其成像效率也会得到很大的提高.

图1

Fig. 1

图 1 特征波场偏移成像示意图 Fig. 1 sketch of characteristic wave imaging

传统的波动方程走时反演方法(Luo和Schuster，1991)利用互相关方法估计观测数据与模拟数据的时差，需要对两者进行加窗处理，从而提取初至波或者折射波.该方法需要拾取特定的波型并在时间域加窗，繁重的拾取工作限制了该方法的应用.

利用特征波场合成，单道地震记录被分离为一系列单独的震相.此时，可直接利用独立的震相进行反演，无需对地震记录加窗处理.更重要的是，特征数据使得多尺度、逐层反演成为可能，如先利用折射波反演浅层速度，再利用由浅到深的反射波反演中深层速度.

由于特征波场只含有单独的波型(震相)，可用互相关函数(Luo和Schuster，1991)测量观测数据与模拟数据的时差，无需对地震数据加窗.本文中，我们将特征波场合成算法应用初至波场分解，并应用于透射波走时反演中.至于特征波场分解在反射波走时反演中的应用，可以参考冯波(2015)，本文不再赘述.

若观测数据 d 与模拟数据 u 经过特征波场合成后，得到的初至震相分别可以表示为 d ^DBF 与 u ^DBF .此时，观测数据 d ^DBF 与模拟数据 u ^DBF 的互相关函数可以表示为

观测数据与模拟数据的初至走时残差 Δt 定义为

基于走时残差L2范数的目标函数可以表示为

经过推导，目标泛函梯度(冯波，2015)为

其中，g₀(x，－t| x_r，0)代表逆时传播的格林函数，u ₀ 为入射场，m 为慢度平方，y(t)为伴随源：

对于初至走时反演，d ^DBF与 u ^DBF 代表对初至波的特征波场合成，因此无需对数据进行加窗处理.

梯度类的迭代算法为：

其中，α^k 为步长，P 为梯度预条件算子.

根据上文的分析，我们可以利用特征波场合成实现不同波型(震相)的分离，而分离后的单独震相可以用于后续的成像及反演.

在数值实验中，我们首先分析不同信噪比的数据对射线参数估计以及特征波场合成的影响，然后，用一个简单模型(凹陷模型)展示特征波场成像的效果及其优势.最后，用特征波场合成方法提取初至波场，并用于波动方程初至走时反演中.

我们用Marmousi模型正演观测数据.速度模 型为：X方向网格数目nx=737，Z方向网格数目nz=250； 网格间隔均为12.5 m.观测系统设置为：第一炮的坐标为625 m，炮间隔25 m，总共301炮.检波器为固定排列，第一道的坐标为625 m，检波器间隔12.5 m，总共601道.用10 Hz主频的Ricker子波作为震源.

首先，选取中心地震道的炮点和检波点坐标为(1850，3000).图 2a中，从左到右分别是x_s0=1850 m 处的共炮点道集，局部空间窗内(x_r0=3000 m)的地震道及其倾斜叠加得到的 τ-p 谱.从 τ-p 谱中可以拾取局部平面波的出射射线参数p_r，它代表了局部平面波在地表的出射方向.我们拾取了能量最大的4个包络对应的出射射线参数.

图2

Fig. 2

图 2 (A)和(B)分别为共炮点/检波点道集，局部窗内的地震道以及局部平面波的入射/出射线参数谱 Fig. 2 (A) and (B) is the common shot/receiver gather, local seismic traces and the ray parameter spectrum, respectively

为了测量局部平面波的入射射线参数，需要抽取共检波点道集.图 2b中，从左到右分别是x_r0=3000 m处的共检波点道集，局部空间窗内(x_s0=1850 m)的地震道及其倾斜叠加得到的 τ-p 谱.从 τ-p 谱中可以拾取局部平面波的出射射线参数p_s，它代表了局部平面波在地表的入射方向.与图 2a中对应，我们拾取能量最大的四组入射射线参数.

利用图 2中拾取的入射和出射射线参数，根据公式(5)，可进行特征波场合成.图 3a中，第1—4道为合成的特征波场.其中，第1道代表对直达波的特征波场合成，它将直达波从地震记录中提取.第2—4道代表对反射波的特征波场合成，它将反射波从地震记录中提取.第5道为重构的中心道(第1—4道叠加)，与原始中心地震道(第6道)相比，只保留了我们提取的震相，提高了特征波的信噪比.同时，考虑到特征波场(第1—4道)的时间局部特征，可以(在变换域)进一步压缩存储.

图3

Fig. 3

图 3 Marmousi模型正演记录的特征波场合成对比 (a) 无噪声; (b) 加入高斯白噪声，S/N=100; (c) 加入高斯白噪声，S/N=20 Fig. 3 Comparison of characteristic wavefield decomposition

接着，为了测试噪声对特征波场合成的影响，我们对原始数据加入了高斯白噪声(S/N=100).图 4中展示了相应的道集及其 τ-p 谱(所有参数含义同图 2).可以看出，高斯随机噪声(S/N=100)降低了 τ-p 谱的分辨率，因而可能会对 τ-p 谱的自动拾取有一定的影响，但对人工拾取的影响较小.图 3b中，展示了对图 4中的数据进行特征波场合成的结果(所有参数同图 2a).显然，特征波场合成之后，仍然能够得到高信噪比的单独震相.这是由于特征波场合成的过程中，同相的信号相互增强，而其他信号相互干涉.若进一步增强噪声强度(S/N=20)，此时，随机噪声完全压制了有效信号，即使在 τ-p 谱也无法拾取射线参数.若我们仍然用真实的射线参数(图 2中拾取的结果)，特征波场合成结果如图 4c所示.其中，直达波的波形仍然具有较为保真的形态，而后续的几个反射波则受到了一定程度的噪声干扰，略有震荡.

图4

Fig. 4

图 4 (A)和(B)分别为原始数据加入高斯随机噪声后的共炮点/检波点道集，局部窗内的地震道以及局部平面波的入射/出射线参数谱 Fig. 4 (A) and (B) is the common shot/receiver gather with Gaussian random noise, local seismic traces and the ray parameter spectrum, respectively

上述试验说明，本文中的特征波场合成算法十分稳健，随机噪声对其影响较小.当噪声能量逐渐增强时，τ-p 谱的分辨率会逐渐降低，直至无法拾取准确的射线参数，进而导致特征波场合成算法失效.

本文用洼陷模型测试特征波场成像方法的有效性，其速度模型如图 5a所示.模拟数据共201炮，炮间距20 m，每一炮301道接收，道间距10 m.本节分 别对该模型数据进行传统的Kirchhoff PSDM和特征波成像来进行效率效果的对比.图 5b是Kirchhoff PSDM的成像剖面，由于其成像过程是通过等时面叠加收敛绕射波，成像剖面上留下了画弧噪声.图 5c是特征波场成像得到的成像剖面，由于其通过对特征域数据进行“搬家”，其画弧范围控 制在局部范围内，其画弧噪声明显弱于传统Kirchhoff PSDM. 由于特征波的成像过程中可以得到射线方向，张角道集可以方便的输出.传统Kirchhoff PSDM一般只能输出地表偏移距道集，在对后续的速度分析或层析成像的适应性方面上来说，张角道集质量要优于地表偏移距道集.图 5(d，e)分别为洼陷模型3500 m处的Kirchhoff PSDM地表偏移距道集和特征波成像的张角道集，对比两图可见，张角道集的成像质量要高于地表偏移距道集.

图5

Fig. 5

图 5 速度模型(a)，Kirchhoff PSDM的成像剖面(b)，特征波场成像剖面(c), Kirchhoff PSDM的成像点道集(d)，特征波场成像点道集(e)及面向目标成像结果(f) Fig. 5 Velocity model (a), image produced by Kirchhoff PSDM (b),characteristic wavefield PSDM (c), Kirchhoff PSDM common imaging point gathers in offset domain, characteristic wavefield PSDM angle-domain CIGs, target- oriented migration result (f)

特征波成像由于只对特征数据进行反投影(特征数据是对原地震数据的一个压缩处理)，其成像效率要远远高于传统的Kirchhoff PSDM或者传统的射线束成像.本文对比了Kirchhoff PSDM和特征波成像在洼陷模型上的成像效率，具体数据见表 1，数据得到了有效的压缩，同时成像效率也有近一个数量级的提高.由于特征波场包含了数据方向信息，特征波场成像方法可以方便的推广到目标体成像，图 5f展示了仅仅对最后一个反射层的面向目标成像结果.对特定目标成像将会进一步提高效率，缩短处理周期.

表1(Table 1)

表 1 K-PSDM 和 CWI 效率对比 Table 1 Comparison of the computational time of K-PSDM and CWI KPSDM CWI 运行时间 13.94(min) 1.46(min) 加速倍数: 9.55 输入数据规模 591(M) 11(M) 数据压缩倍数: 53.7

表 1 K-PSDM 和 CWI 效率对比 Table 1 Comparison of the computational time of K-PSDM and CWI

在本节中，我们将特征波场合成算法分离初至波型，并用于波动方程初至走时层析中.速度模型(BP2004速度模型的右端部分，横向约26 km，深度约12 km)如图 6a所示.观测系统为：炮点深度12.5 m，炮点间隔为200 m，共107炮.检波器深度12.5 m，间隔400 m，每炮54道.观测系统总共5778道.初始速度模型采用一维模型，其中，水体速度以及水底深度已知.总共进行18次反演迭代，反演得到的速度模型见图 6b.反演得到了整体背景速度的变化趋势，其中模型右侧的高速隆起反演得较好，但没有反演出浅层小尺度的低速异常.图 7中为初至波走时对比曲线.其中，红色、绿色和蓝色曲线分别是真实模型、初始模型和第18次反演的速度模型中的初至走时.可以看出，真实的初至走时与反演得到的速度模型中的初至走时基本吻合，实现了初至走时反演的目标.

图6

Fig. 6

图 6 速度模型(a)及反演得到的速度模型(b) Fig. 6 The true velocity model (a) and the inverted result (b)

图7

Fig. 7

图 7 初至波走时对比，背景为观测的炮集((a)和(b)分别为两个不同位置的炮集).红色、绿色和蓝色曲线分别是真实模型、初始模型和第18次反演的速度模型中的初至走时 Fig. 7 The comparison of first-arrival travel-times. The red, green and blue curves are the true, initial and inverted first-arrival travel-time. The background of (a) and (b) is the observed shot gather from different location

本文提出了一种地震数据的稀疏表达方法.不同于传统的变换域中的数据稀疏表达理论，该方法利用局部平面波的传播方向(慢度矢量)，在中心炮检点处同时进行波束合成，从而将地震数据投影到局部平面波域(高维空间，二维观测对应5D数据；三维观测对应9D数据)中.由于波束合成后的地震数据描述了局部平面波的方向特征，因此我们称之为特征波场.相应的，我们称这种波场合成(分离)算法为特征波场合成方法.

特征波场合成的思想是利用局部平面波到达不同检波器的时差，对地震记录进行时移叠加，使得给定的入射/出射慢度的局部平面波相干加强而其他信号相互抵消，从而分离地震记录中的不同的波型(震相).因此，该方法对含有噪声的地震数据有较好的适应能力.文中的特征波场合成正问题需要估计射线参数，当地震数据受噪声干扰严重时，难以在常规 τ-p 谱中自动估计局部平面波的射线参数(慢度矢量).此时，可以考虑基于反演理论的特征波场合成方法.在压缩感知理论框架下同时反演局部平面波及其传播方向，从而提高特征波合成的自动化程度并保持方法的稳健性.

然而，特征波场合成正问题要求局部炮/检排列在同一水平面上.因此，对于地表高程剧烈变化的起伏地表，该方法不再成立.此外，若炮点(检波器)采样过于稀疏，也会影响射线参数估计的精度.此时，可以仅考虑单次波束合成.

通过特征波场合成，可以将地震数据分解为单独的震相.此时，传统的Kirchhoff及高斯束偏移方法可以实现从“沿等时面的画弧”到向“向反射点归位搬家投影”的转变.相应的，常速介质中的偏移脉冲响应由“椭圆/球”退化为点响应.同时，由于在偏移过程中不需要对炮/检射线参数扫描，因而叠前偏移效率将会大幅提高.若炮点的慢度矢量(入射射线参数)难以估计，特征波场的偏移成像可以自动退化为高斯束类的成像方法.

更进一步地，特征波偏移成像这种一对一映射的偏移方法非常适合与基于角度道集的速度层析反演方法融合成一体.它们具有同样的数据空间和成像空间.这种特征波场驱动的高效偏移成像和速度层析反演融为一体的技术体系会改变目前的地震成像流程，与大规模计算机集群和现代图像处理方法结合，在很大程度上可以实现半自动或自动的地震波成像.