地球物理学报  2015, Vol. 58 Issue (2): 674-684   PDF    
L1范数约束被动源数据稀疏反演一次波估计
程浩1, 王德利1, 冯飞1,2, 王通1    
1. 吉林大学地球探测科学与技术学院, 长春 130026;
2. 中海油服物探事业部数据处理中心, 天津 300451
摘要:对于被动源地震数据,运用常规的互相关算法得到的虚拟炮记录中,不仅含有一次波反射信息,还包括了表面相关多次波.然而,通过传统的被动源数据稀疏反演一次波估计(EPSI)方法,可以求得只含有一次波,不含表面相关多次波的虚拟炮记录.本文改进了传统的被动源数据稀疏反演一次波估计问题的求解方法,将被动源稀疏反演一次波估计求解问题转化为双凸L1范数约束的最优化求解问题,避免了在传统的稀疏反演一次波估计过程中用时窗防止反演陷入局部最优化的情况.在L1范数约束最优化的求解过程中,又结合了2D Curvelet变换和小波变换,在2D Curvelet-wavelet域中,数据变得更加稀疏,从而使求得的结果更加准确,成像质量得到了改善.通过简单模型和复杂模型,验证了本文提出方法的有效性.
关键词被动源     稀疏反演     L1正则化     凸优化    
Estimating primaries by sparse inversion of passive-source seismic data with L1-norm constraint
CHENG Hao1, WANG De-Li1, FENG Fei1,2, WANG Tong1    
1. College of Geo-Exploration Science and Technology, Jilin University, Changchun 130026, China;
2. Data Processing Center of Geophysics, China Oilfield Services Limited, Tianjin 300451, China
Abstract: This work has improved the original algorithm to estimate primaries by sparse inversion of passive-source seismic data through replacing the original algorithm to solve the convex optimization problem with L1-norm constraint. It avoids using a time-window to prevent the inversion from into local optimization situations when estimating primaries by sparse inversion. Moreover, during the solving of the optimization problem with L1-norm constraint, 2D Curvelet transform and wavelet transform are used at the same time. In 2D Curvelet and wavelet domains, the data become more sparse, then the results obtained are more accurate and the quality of imaging is improved.
First, the method of convex optimization problem with L1-norm constraint is introduced to solve the problem of estimating primaries by sparse inversion of passive-source seismic data, instead of the steepest descent method under L0-norm constraint. Second, 2D Curvelet transform and wavelet transform are combined during the sparse inversion. In the 2D Curvelet-wavelet domain, the data become more sparse. Comparing with 3D Curvelet transform, the velocity of 2D Curvelet-wavelet transform is improved. Third, a simple model and a complex model are used to simulate the passive seismic data. The method of convex optimization problem with L1-norm constraint and that combined with 2D Curvelet transform and wavelet transform are used to estimate primaries from the passive seismic data, respectively. At last, comparison with the results obtained by the traditional LSQR algorithm illustrates that the method proposed is feasible and effective.
The method of estimating primaries by sparse inversion can directly estimate primaries from the passive seismic data, and obtain the virtual-shot gathers which are free of the surface-related multiples. Under the assumption that the data is sparse, this work uses the method of convex optimization problem with L1-norm constraint to replace the traditional one to estimate primaries, which avoids using a time-window to prevent the inversion from into local optimization situations during sparse inversion, and improves the precision of the primaries estimated. It also suppresses artificial influence and improves the imaging quality.
Comparing with the result obtained by lsqr algorithm shows the accuracy and superiority of the convex optimization problem with L1-norm constraint. During sparse inversion by L1-norm inversion, 2D Curvelet transform and wavelet transform are combined to make the data more sparse, and improve the precision of primaries estimated. At the same time, the artificial influence suppression is improved further.Comparing with the traditional method to estimate primaries, the convex optimization problem with L1-norm constraint can avoid using a time-window to prevent the inversion from into local optimization situations, and the primaries estimated become more accurate. 2D Curvelet transform and wavelet transform introduced make the data sparse, and improve the precision of primaries estimated.
Key words: Passive sources     Sparse inversion     L1 regularization     Convex optimization    
1 引言

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

Van Groenestijn和Verschuur(2009a2009b)首次提出常规主动源稀疏反演一次波估计(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范数约束的一次波响应.同时,运用了传统的最小二乘进行求解,与上述求解的一次波结果进行了对比.

2 一次波估计(EPSI)

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

式中,X0表示一次波反射系数序列,S+表示震源子波特性函数,两者在频率域相乘给出了一次波记录P0=X0S+.注意,在本文中,所谓的一次波指的是那些没有在地表发生过反射的响应(包含层间多次波).一次波反射系数X0与自由表面反射算子R和总的上行波数据P-进行多维相乘,就得到了地表相关多次波M=X0RP-.

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

分析方程(2),同时存在两个未知数X0和S.利用一个反演处理的过程,同时对两个未知数进行估计0,进而分解总的上行波场P-.由此,可以得到目标函数J的表达式

式中,i表示迭代次数,表示整个矩阵元素的相加,表示所有频率的相加.用EPSI方法通过迭代优化处理来求解这样的问题.假设R=- I .在这种算法的第一次迭代中,均设置0的初始值为0.

首先,更新一次波反射系数X0.利用最速下降方法更新ΔX0

式中,(i I +RP-)H是(i I +RP-)的复杂共轭.(P --0,ii0,iRP-)可以被看做是未分解的数据或者残差.由于在第一次迭代中X0的值为0,所以,第一步等于数据和它本身的一个多维褶积,即P-(R P-)H.

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

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

式中,α是一个正的频率独立因子,来衡量迭代步长.选择步长因子α是为了让目标函数值逐渐
减小. 以更新0同样的方式,
更新

从全矩阵ΔS之中,选择对角元素,求平均获得步长ΔS.此处,将ΔS变回时间域,它的长度是有限的.在此之后,为了确保下一次目标函数方程(3)的值比上一次的低,每次都需计算更新以后的目标函数.

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

3 被动源数据EPSI

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

式中,Ppas是在被动源数据采集中观测到的上行波 场数据,Pdir是来自于近地表被动源到检波器的直达波响应.方程(7)右侧的第一项,与方程(2)右侧第一项不同,不能用一次波脉冲响应表示,而两式的第二项一样,都利用一次波反射系数序列与表面反射算子和总体上行波场进行相乘.将方程(7)进行移项处理,可得Pdir=Ppas-X0RPpas.在时间域,由于对X0的稀疏约束,不能同时用来求解PdirX0两个未知数.因此,就需要一个额外的假设条件对其进行约束.在这里,假设Pdir具有最小的能量.由此,可以得到一个趋于最小的目标函数:

式中,令R=-I.开始第一次迭代的时候,同样设置X0为0,就得到X0的更新,如下式:

第一次迭代的时候,更新就等于数据本身的多维相关Ppas(RPpas)H,得到的结果和地震波干涉的结果相同.为了避免反演过程中,陷入局部最优解的问题,需要在迭代过程中利用一个时窗进行约束.
图 1 主动源一次反射波响应(a)及由被动源到地表检波器的直接波响应(b)Fig. 1 The primary of control source(a) and the direct wave of passive seismic data(b)
4 L1范数稀疏约束

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

式中,A是一个线性算子,BlockDiagω表示频率域分块对角矩阵,x0为一次波反射系数,“~”表示估计值,p为总的上行波场,x0和p都为列向量,S为震源子波,P为总的上行波场,I表示自由表面算子,ft为时间域傅里叶变换,ft*为傅里叶反变换,为Kronecker积,k为迭代次数,τ为0范数值.

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

方程(12)中,τ为1范数值.求解方程(12)的最优化问题,可以转化为

式中,σ为实测数据与估计出的数据的差值.利用谱梯度投影算法(SPGL1)求解对方程(12)和(13)的最优化反演问题.

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

式中,c0为一次波反射系数的Curvelet系数,C*为3D Curvelet反变换.在三维曲波域中,一次波反射系数具有较大的Curvelet系数,使数据变得更加稀疏和具体.在反演过程中,通过求取一次波反射系数的Curvelet系数,进而求得一次波反射系数,从而可以获得更加精确的解.但由于需要进行三维曲波正变换和反变换,耗时较长,在精度提高的同时,计算速率大大降低.

5 L1范数稀疏约束被动源数据EPSI

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

式中各参数的物理意义同前.为了提高求解的精度,本文将二维曲波变换和小波变换结合,代替三维曲波变换,引入到求解过程当中.这样,在使数据变得更加稀疏的同时,运算速度有所提高,得

式中,s0为经过2D Curvelet和小波变换后的系数,S为结合的2D Curvelet变换和小波变换,S*S的反变换.通过逐次的迭代,得出最终的估计值 0,将其与任意震源子波褶积,就能得到一次波记录.

6 模型试算 6.1 简单模型

首先利用一个二维的简单模型进行验证,模型如图 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 Simple model

图 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 最小二乘求解结果(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)
6.2 复杂模型

在经过简单模型验证的前提下,本文又利用一个较为复杂的模型进行验证,如图 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 Complex model

图 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 小二乘求解结果(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))

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

7 结论与展望

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

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

参考文献
[1] Baardman R H, Verschuur D J, Van Borselen R G, et al. 2010. Estimation of primaries by sparse inversion using dual-sensor data. 80th Annual International Meeting, SEG, Expanded Abstracts,3468-3472.
[2] Claerbout J F. 1968. Synthesis of a layered medium from its acoustic transmission response. Geophysics, 33(2): 264-269.
[3] Feng F, Wang D L, Zhu H, et al. 2013. Estimating primaries by sparse inversion of the 3D curvelet transform and the L1-norm constraint. Applied Geophysics (in Chinese), 10(2): 201-209.
[4] Lin T T Y, Herrmann F J. 2011. Estimating primaries by sparse inversion in a curvelet-like representation domain. 73rd EAGE Conference and Exhibition, EAGE, Extended Abstracts, H043.
[5] Lin T T Y, Herrmann F J. 2013. Robust estimation of primaries by sparse inversion via one-norm minimization. Geophysics, 78(3): R133-R150.
[6] Van Den Berg E, Friedlander M P. 2008. Probing the pareto frontier for basis pursuit solution. SIAM Journal on Scientific Computing, 31(2): 890-912.
[7] Van Groenestijn G J A, Verschuur D J. 2009a. Estimating primaries by sparse inversion and application to near-offset data reconstruction. Geophysics, 74(3): A23-A28.
[8] Van Groenestijn G J A, Verschuur D J. 2009b. Estimation of primaries and near-offset reconstruction by sparse inversion: marine data application. Geophysics, 74(6): R119-R128.
[9] Van Groenestijn G J A, Verschuur D J. 2010. Estimation of primaries by sparse inversion from passive seismic data. Geophysics, 75(4): SA61-SA69.
[10] Van Groenestijn G J A, Ross W. 2011. Primary estimation on OBC data by sparse inversion. 81th Annual International Meeting, SEG, Expanded Abstracts, 3531-3535.
[11] Vasconcelos I, Snieder R. 2008a. Interferometry by deconvolution, Part 1—Theory for acoustic waves and numerical examples. Geophysics, 73(3): S115-S128.
[12] Vasconcelos I, Snieder R. 2008b. Interferometry by deconvolution: Part 2—Theory for elastic waves and application to drill-bit seismic imaging. Geophysics, 73(3): S129-S141.
[13] Wapenaar K, van der Neut J, Ruigrok E. 2008. Passive seismic interferometry by multidimensional deconvolution. Geophysics, 73(6): A51-A56.
[14] Wapenaar K, van der Neut J, Ruigrok E, et al. 2011. Seismic interferometry by cross-correlation and by multi-dimensional deconvolution: a systematic comparison. Geophysical Journal International, 185(3): 1335-1364.
[15] Ying L X, Demanet L, Candes E J. 2005. 3-D discrete curvelet transform. Proceedings SPIE 5914, Wavelets XI. San Diego: SPIE, 344-354.
[16] Zhu H, Wang D L, Shi Z A, et al. 2012. Passive seismic imaging of seismic interferometry. Progress in Geophysics (in Chinese), 27(2): 496-502, doi: 10.6038/j.issn.1004-2903.2012.02.012.
[17] 冯飞, 王德利, 朱恒等. 2013. 三维曲波变换L1范数约束稀疏反演一次波估计方法研究.应用地球物理, 10(2): 201-209.
[18] 朱恒, 王德利, 时志安等. 2012. 地震干涉技术被动源地震成像. 地球物理学进展, 27(2): 496-502, doi: 10.6038/j.issn.1004-2903.2012.02.012.