在被动源地震勘探中，激发利用的是地下的微震动，例如，天然地震，构造破裂等，震源类型比较复杂，激发位置也比较随机.地震波干涉技术(Claerbout，1968)，作为一种非常有效的方法，能够将被动源数据进行转化，形成地表激发、地表接收的虚拟炮集记录，从而获得地下的构造信息.利用被动源地震勘探，不但对地表条件没有较多的要求，而且节约了大量的生产成本，很快成为了人们关注的热点研究问题.互相关算法(朱恒等，2012)是地震波干涉技术最早使用的算法，也是最常用的一种算法，但它依赖于地表被均匀照明的假设，在激发角度和强度方面需要具有一致性.在实际采集中，这种假设很难实现，互相关算法就不能给出振幅准确的近地表响应.随后，地震波干涉技术又发展了反褶积算法(Vasconcelos和Snieder， 2008a，2008b)，相比于互相关算法，在成像效果上有所提高.Wapenaar等(2008)为了克服这种影响，在被动源直达波已知的前提下，利用多维反褶积算法进行被动源地震波干涉技术，成功地对近地表响应的振幅进行了校正.2012年，Wapenaar(2011)又从理论推导，对应的物理意义，以及理论与实际数据的试算等方面，对互相关算法和多维反褶积算法进行了一个较为系统的对比.通过互相关算法得到的虚拟炮集记录，不仅包含一次波信息，同时也包含了表面相关多次波和层间多次波，这就对后续虚拟炮道集的处理解释，造成了一定的影响.

Van Groenestijn和Verschuur(2009a，2009b)首次提出常规主动源稀疏反演一次波估计(Estimating Primaries by Sparse Inversion，EPSI)方法，和表面相关多次波消除方法(Surface-Related Multiple Elimination，SRME)一样，都是基于表面多次波模型，完全数据驱动，不需要知道地下的任何先验信息，在一次波反射系数稀疏的假设条件下，对数据驱动的波场信息进行反演，直接对一次波进行估计，这就避免了从原始数据中匹配减去多次波的过程，取而代之的是一个大规模的反演过程.通过迭代更新一次波响应，利用稀疏约束，一次波和它们相关的多次波分别与原始数据进行匹配，同时进行一次波和多次波的分离，完全避免了自适应减去的过程.2010年，Van Groenestijn和Verschuur(2011)在传统主动源稀疏反演EPSI方法的基础上，将其算法进行了一定程度修改，给出了对被动源数据进行一次波估计的方法.这种方法在被动源数据稀疏假设条件成立的前提下，直接获得不含表面相关多次波的一次波响应数据，避免了对虚拟炮集记录进行表面相关多次波匹配相减的过程.

稀疏反演EPSI方法得到了快速发展和广范应用.EPSI算法不断地被进行修改，以便应用于各种方法的地震数据. Baardman等(2010)将EPSI应用于双检波器数据处理；Van Groenestijn等(2011)又将EPSI应用于OBC数据处理.

Lin和Herrmann(2013)提出对常规主动源数据利用L1范数约束稀疏反演一次波估计，把求取稀疏解作为明确的目标，使传统的EPSI方法得到了进一步的增强.这种方法仍然属于梯度求解范畴，但它是从一个双凸优化结构推导过来的，基于一个广义的基本追踪(BP)去噪算法(Van，et al，2008).这样就避免了在传统反演中需要过多地对参数进行调整，而且还提高了成像质量.同年，冯飞等(2013)将三维曲波变换(Ying，et al，2005)结合到主动源L1范数约束稀疏反演一次波估计当中，在三维曲波域中，数据变得更加稀疏，所求得的解变得更加准确，但由于在计算过程中，需要进行三维曲波正、反变换，所以这种方法的计算速率往往较低.

本文首先回顾了传统主动源数据稀疏反演EPSI方法及被动源数据稀疏反演EPSI方法的基本理论；然后将双凸L1范数约束的最优化求解方法，引入到了被动源数据稀疏反演一次波估计问题中，替代了原始L0范数约束的最速下降求解方法；其次，在L1范数约束的最优化求解过程中，结合了二维曲波变换和小波变换(Lin，et al，2011)，将数据变到了更加稀疏的二维曲波与小波结合的域中进行求解，相比于三维曲波变换，二维曲波变换和小波变换结合的正、反变换，在速度上有所提高；最后，运用二维模型进行验证，分别求取了L1范数约束和结合了二维曲波变换及小波变换的L1范数约束的一次波响应.同时，运用了传统的最小二乘进行求解，与上述求解的一次波结果进行了对比.

在常规地表主动源地震勘探中，采集到的地震记录通常含有一次波和表面相关多次波，上行波场P^-在频率域可以表示为

在这里，令S⁺=S(ω)I，就意味着所用震源都采用一个常数震源子波，并且忽略了震源的方向特性，因此方程(1)被改写为

首先，更新一次波反射系数X₀.利用最速下降方法更新ΔX₀：

为了约束这个反演过程，在逐次迭代过程中，使₀的更新加上稀疏性约束.

在时间域中，假设X₀可以用有限数目具有大振幅的尖脉冲响应(来自于主要反射层)，和一些小振幅脉冲响应(来自于其他)来表示.在时间域更新₀的时候，需要设置一个时间窗，选取每一道的强反射响应.在每次迭代中，通过调整窗的大小，使得目标函数得到收敛.窗口里面尽可能地不包含与一次波无关的响应.出现在一次反射之前的人工影响也不能在窗口内出现.然后，将稀疏更新的ΔX₀加到一次脉冲响应中

通过反复地交替优化，进行₀和的更新，
直到目标函数值达到一定的要求为止，就得到了一次反射系数和震源子波的最终估计值.

被动源数据与主动源数据最大的区别，就是在接收到的被动源数据中，没有类似于常规数据的一次波响应，取而代之的是一个来自于被动源的直达波响应，如图 1所示.这个直达波响应引起了一系列与之相关的表面相关多次波响应.可以用一个与方程(2)相似的表达式描述:

图 1

Fig. 1

图 1 主动源一次反射波响应(a)及由被动源到地表检波器的直接波响应(b)Fig. 1 The primary of control source(a) and the direct wave of passive seismic data(b)

第2节所述的主动源EPSI，是一个L0范数约束下求解一次波反射系数的最优化问题：

为了解决L0范数在迭代过程中收敛和稳定性问题，在足够稀疏的条件下，可以利用L1范数约束最优化问题求解方法代替传统L0范数约束最速下降方法，进行上述问题的求解.如:

将三维曲波变换引入到L1范数约束的EPSI:

本文将L1范数稀疏约束双凸优化问题求解方法引入到被动源数据稀疏反演EPSI求解问题中，解决了原始被动源数据稀疏反演EPSI过程中加时窗问题.根据被动源数据一次波估计的目标函数方程(8)，利用L1范数稀疏约束最优化问题进行求解，得

首先利用一个二维的简单模型进行验证，模型如图 2所示，其中包含了一个倾斜界面和一个水平界面.模型长为3000 m，深为1000 m，检波点设置在地表，道间距为15 m，脉冲震源位于水平层下面，位置随机，进行激发. 图 3(a—c)分别为随机震源横向位于350 m，1500 m，2650 m，深度位置不同的被动源地震记录.图 3(d—f)为表面设置为吸收表面的标准采集的地震记录，不含有表面相关多次波，与图 3(a—c)的横向位置相对应，已对直达波进行切 除处理.用常规的互相关算法地震波干涉技术对接 收到的被动源记录处理，得到与之对应的虚拟炮记录(图 3(g—i))，可以看出，除了一次波信息以外，虚拟炮记录中还包含了表面相关多次波信息和其他一些人为计算产生的伪波等.

图 2

Fig. 2

图 2 简单模型示意图Fig. 2 Simple model

图 3

Fig. 3

图 3 随机震源记录(a—c)；原始数据记录(d—f)；虚拟炮记录(g—i)Fig. 3Rand om-source records(a—c); original-source records(d—f); virtual-source records(g—i)

利用本文提出的L1范数稀疏约束最优化问题求解方法对被动源数据一次波进行求解，得到的结果如图 4(d—f)所示.再将2D Curvelet变换和小波变换结合到L1范数稀疏约束最优化问题求解过程中，得到的结果如图 4(g—i)所示.为了验证L1范数稀疏约束的准确性及优越性，本文又用L2范数约束的最小二乘方法对被动源稀疏反演一次波估计进行求解，得到的结果如图 4(a—c)所示.

与标准采集的炮记录进行对比，由于被动震源照明不均匀，远偏移距照明不足，震源孔径有限，远偏移距处没有震源覆盖的相位点，导致估计出的一次波记录在远偏移距处的走时误差较大.分析以上算法的结果可以发现，三种方法都比较准确地估计出了虚拟炮记录的一次波响应，但也都含有计算过程中产生的人为影响；从计算速度方面看，L2范数稀疏约束速度最快，依次为L1范数稀疏约束和结合了2D Curvelet变换和小波变换的L1范数约束方法；从成像质量方面看，图 3中对应位置的箭头所示，结合了2D Curvelet变换与小波变换的L1范数约束方法对人为影响的压制最好，成像质量最高，然后依次为L1范数稀疏约束和L2范数稀疏约束.可见，在数据稀疏的条件下，L1范数稀疏约束要优越于L2范数稀疏约束，并且数据越稀疏，L1范数稀疏约束求得结果越精确，但计算的时间也就越长.

图 4

Fig. 4

图 4 最小二乘求解结果(a—c)； L1范数约束求解结果(d—f)；结合2D Curvelet和小波变换的L1范数约束求解结果(g—i)Fig. 4 Results of least-square(a—c); results pf L1-norm constraint(d—f); results of L1-norm constraint with 2D Curvelet and Wavelet Transform(g—i)

在经过简单模型验证的前提下，本文又利用一个较为复杂的模型进行验证，如图 5.模型的长为 5000 m，深度为1800 m，其中包含了高速盐丘和断层构造.检波点设置在地表，道间距为20 m，脉冲震源位于模型底层，深度随机分布，进行激发.如图 6(a—c)所示，分别为随机震源横向位于1000 m，2500 m，4000m，深度位置不同的被动源地震记录，从图中可以发现，由于高速盐丘的影响，对应位置处的透射波和反射波能量较弱.图 6(d—f)为表面设置为吸收表面的标准采集的地震记录，不含有表面相关多次波，与图 6(a—c)的横向位置相对应.利用基于互相关算法的地震波干涉技术对接收到的被动源记录进行处理，得到对应的虚拟炮记录，如图 6(g—i)，可以看出，除了一次波信息之外，虚拟炮记录中包含了许多多次波信息.同时，由于高速盐丘和断层构造的影响，虚拟炮记录中反射波反应的信息比较复杂，而且不同位置处能量存在着一定的差异，靠近盐丘凸起的地方能量较弱，其他地方相对较强.

图 5

Fig. 5

图 5 复杂模型示意图Fig. 5 Complex model

图 6

Fig. 6

图 6 随机震源记录(a—c)；原始数据记录(d—f)；虚拟炮记录(g—i)Fig. 6 Random-source records (a—c); Original-source records (d—f); Virtual-source records (g—i)

同样，利用最小二乘，L1范数约束和结合了2D Curvelet变换和小波变换的L1范数约束优化方法求解，分别得到了如图 7(a—c)、图 7(d—f)、图 7(g—i)的结果.从成像结果来看，与标准采集的单炮记录进行对比，三种方法得到的虚拟炮记录中一次波结果都比较准确.同时，从图 7(a—c)中箭头位置可以看出，由于模型结构的复杂性，L2范数稀疏约束在计算过程中引人的人为影响是较为严重的.通过L1范数约束方法及结合了2D Curvelet变换和小波变换的L1范数约束方法，使人为影响得到减 弱，成像质量得到提高，如图 7(d—f)和7(g—i)所示.

图 7

Fig. 7

图 7 小二乘求解结果(a—c)； L1范数约束求解结果(d—f)；结合2Dcurvelet和小波变换的L1范数约束求解结果(g—i)Fig. 7 Results of least-square (a—c); Results pf L1-norm constraint (d—f); Results of L1-norm constraint with 2D Curvelet and Wavelet Transform(g—i))

通过上述的分析可知，模型比较复杂的时候，接收到的被动源数据也比较复杂，其中可能含有透射波，反射波，绕射波等，经过被动源稀疏反演一次波估计得到的结果中，引入的人为影响就会比较严重，但进行一定的变换，使数据变得更加稀疏，求得的一次波结果准确性就会有所提高，同时减小了引入的人为影响.

被动源稀疏反演一次波估计方法，直接利用被动源数据估计一次波，给出了只含有一次波，不含表面相关多次波的虚拟炮记录.在数据为稀疏的假设条件下，本文利用L1范数约束最优化问题求解方法替代了原始被动源稀疏反演一次波的求解方法，解决了原始算法中加时窗函数的问题，并且使求得的解变得更加准确，在一定程度上压制了人为影响，使成像质量得到提高.将其与L2范数约束的最小二乘求解结果进行了对比，体现了L1范数约束最优化求解方法的准确性和优越性.在L1稀疏约束求解过程中，又结合了2D Curvelet变换和小波变换，使数据变得更加稀疏，解的精确程度得到了提高，同时对人为影响的压制也得到了进一步改善.与主动源L1范数约束稀疏反演一次波估计方法相比，被动源L1范数稀疏反演一次波估计方法只进行一次波反射系数的估计，无需对震源子波进行估计，更不需要了解震源的属性，简化了求解过程.但是，由于被动震源照明不均匀，远偏移距照明不足，震源孔径有限，远偏移距处没有震源覆盖的相位点，造成了估计的一次波记录远偏移距走时误差较大.

巨大的勘探成本和复杂的地表条件已经成为了目前常规地震勘探的难题.然而被动源地震勘探，无需人工激发，利用地下天然源进行激发，合成虚拟炮记录，对地下的地质构造进行成像，既保护环境，又节约大量生产成本.直接对被动源数据进行一次波估计，避免了对复杂的被动源数据进行表面相关多次波的预测减去过程，同样节省生产成本，存在着一定的实际意义.