2017, Vol. 60
2. 西安交通大学电子与信息工程学院, 西安 710049;
3. 海洋石油勘探国家工程实验室, 西安 710049;
4. 中海石油研究中心, 北京 100029
2. Xi'an Jiaotong University, School of Electronic and Information Engineering, Xi'an 710049, China;
3. National Engineering Laboratory for Offshore Oil Exploration, Xi'an 710049, China;
4. Research Center of CNOOC, Beijing 100029, China
反射地震勘探得到的反射地震数据反映的是地下介质中波阻抗分界面对地表人工激发的地震波的反射效应.为了进一步刻画地下介质的岩性变化通常需要利用反射地震数据进行地层参数反演，如，波阻抗反演等.本文将地下介质视为声学介质，下文提到的波阻抗均指声波阻抗.通过声波阻抗反演把反射地震数据转换成声波阻抗数据以后，可以直接与已知地质、测井信息对比标定，便于解释人员有效地对地层岩性参数的变化进行研究，得到岩性参数在空间的分布规律，指导油气的勘探开发.可见，发展高精度的声波阻抗反演方法对于精细刻画地层岩性具有重要意义.
声波阻抗反演技术在20世纪70年代出现，80年代开始得到蓬勃发展.经历了从直接反演到模型反演，从叠后反演到叠前反演，从线性反演到非线性反演的发展过程(张永刚，2002).按照介质模型的不同，声波阻抗反演可以分为基于波动方程的反演和基于褶积模型的反演两大类.前者由于算法结构复杂、计算量大、抗干扰能力差，还没有在实际中大规模应用.而基于褶积模型的声波阻抗反演方法由于算法简单易行、对地震噪声敏感性小，因此在生产中得到了广泛应用(崔岩等，2009；邸海滨等，2011).本文只关注基于褶积模型的声波阻抗反演方法.基于褶积模型的声波阻抗反演方法可分为三类，即：道积分，基于模型的反演和递推反演.道积分反演是最简单的声波阻抗反演方法，它通过积分方法把地震道记录直接转换为声波阻抗数据，反演过程具有积分误差小、计算简单、对计算机要求较低的优点，但是道积分反演过程仅仅依赖于地震资料，地震资料的品质和带宽对反演结果影响很大.基于模型的反演方法需要先建立一个初始地层声波阻抗模型，然后由此模型进行地震正演，求得合成地震记录，并与实际地震记录相比较，根据比较结果修改阻抗模型，再做模型正演，如此反复多次迭代直至合成地震记录与实际记录接近，从而得到地下的声波阻抗模型.模型反演方法是在模型基础上进行的正演计算，反演过程对地质模型有很高的依赖性.递推反演方法首先利用实际地震记录估计反射系数序列，然后用递推方法求取地层声波阻抗，方法计算速度快、对计算机要求较低.递推反演结果能够反映地下声波阻抗的分布规律.利用测井数据校正地层反射系数振幅，然后对校正后的反射系数递推反演得到的波阻抗做低频校正可得到绝对波阻抗，能够用于储层预测与砂体雕刻.
递推反演中最重要的部分是要由地震资料正确地估算地层反射系数.要从反射地震道中得到反射系数序列，就是要从反射地震道中消除子波效应.由于子波是带限的，这使得上述反射系数反演问题是一个病态问题，需要加约束来消除问题的多解性.由于较大的反射系数在声波阻抗反演中占主要地位，而这些反射系数往往是稀疏的，因此通常会对反射系数加稀疏性约束，从而就有了稀疏脉冲反褶积方法(Latimer et al., 2000; Nguyen, 2008).稀疏脉冲反褶积是实际生产中常用的求取递推反演反射系数的方法.目前已有许多稀疏脉冲反褶积方法：其中一些方法采用不同的搜索策略来定位稀疏脉冲，然后依据反射系数的统计模型构造目标函数求解反射系数(Kormylo and Mendel, 1978)；一些方法采用随机反演，通过非线性优化算法随机搜索脉冲反射系数的位置并计算反射系数的振幅，使得最少数量的脉冲构成目标函数的解(Velis, 2008)；其他一些方法通过加范数约束得到稀疏反射系数(Oldenburg et al., 1983; Debeye and Riel, 1990; Sacchi et al., 1994; Wang, 2011).然而，当上述逐道处理方法应用于地层结构不规则的多维地震数据时，往往会出现横向连续性差，反褶积效果因为噪声干扰和子波估计有误差而大打折扣的现象.虽然学者们提出了一些方法来缓解这一问题(Wang et al., 2006; Lu, 2009; Kazemi and Sacchi, 2014; Guitton and Claerbout, 2015; Nose-Filho, 2016; Yuan et al., 2016)，但是在面对具有复杂结构和横向变化较大的地震数据时，上述方法仍然难以取得令人满意的效果.这也就导致基于这些数据递推反演得到的波阻抗效果不佳.
为此，本文提出了一种矩阵Toeplitz稀疏分解算法，通过对地震数据剖面进行矩阵Toeplitz稀疏分解可同时反演子波和反射系数剖面.这种基于矩阵分解的稀疏脉冲反褶积方法在计算的过程中同时使用了多个地震道估计反演参数.这些参数适应于多道地震数据的特性，从而使得反演得到的反射系数剖面能够保持地震数据的横向连续性并且能够减少随机噪声的影响.然后，由反演得到的反射系数剖面进一步递推反演可以得到高分辨率的相对波阻抗剖面.
2 方法原理 2.1 矩阵Toeplitz稀疏分解令Y表示一个n×m的矩阵，它可被因式分解为一个n×n的矩阵A和一个n×m的矩阵R的乘积，相应的矩阵分解的基本模型可表示为
|
(1) |
这里F表示F范数，即Frobenius范数.为了实现矩阵的Toeplitz稀疏分解，需要对A加Toeplitz结构约束，对R加稀疏约束.
若矩阵A具有Toeplitz结构，则A可表示成下列形式
|
(2) |
令Ik表示除了第|k|条对角线上的元素为1(当k=0时，主对角线元素为1；k>0时，主对角线上方的第k条对角线元素为1；k < 0时，主对角线下方的第|k|条对角线元素为1)，其他元素都为0的n×n矩阵.此时，(2)式可以表示为
|
(3) |
可见，Toeplitz矩阵A可用一组向量α=(a-(n-1), …, a0…, an-1)T表示.
稀疏问题通常可以转化为L0正则化问题来求解.然而，L0正则化问题往往是一个NP难题，很多NP难题没办法在多项式时间内求解，这意味着用任何现有的算法求解该问题所用的时间都会随着问题规模的增大而迅速增加(Ge et al., 2011).为此，许多学者建议松弛L0正则化，转而用L1正则化代替.L1正则化问题能够被转化为二次凸优化问题，求解计算效率高(Xu et al., 2012).因此，本文采用L1范数(Tibshirani, 1996; Tibshirani et al., 2005)来约束矩阵R的列，从而得到稀疏矩阵R.相应的矩阵Toeplitz稀疏分解的基本模型公式如下：
|
(4) |
这里，rj∈ℝn是反射系数矩阵R的第j列，‖ · ‖ 1是向量的L1范数.在实际问题中，A中元素除了具有Toeplitz结构外，也可能具有其他性质，例如，光滑性或者稀疏性.因此，我们有必要给α另加一个约束项，从而得到矩阵Toeplitz稀疏分解的更一般形式：
|
(5) |
为了便于计算，正则化项Ψ(·)选为变量α的凸函数.
公式(5)并不是对所有的变量都是凸的，因此较难求解.但是，当A(或α)或者R确定后，原来的问题就变成两个简单的凸子问题，通过交替迭代求解这两个子问题就可以逼近原来问题的解，这实际上是块坐标下降方法(Tseng, 2001; Xu and Yin, 2013)，Xu和Yin (2013)对该方法的收敛性做了详细的讨论和证明.
当A确定后，公式(5)的优化问题与A无关，此时(5)式可变为
|
(6) |
由于Y和R是可分的，(6)式可继续分为一系列的子问题：
|
(7) |
这里j=1, …, m，yj∈ℝn表示矩阵Y的第j列.(7)式对应的L1范数正则化问题已经被广泛的研究过，有很多算法可以求解(Efron et al., 2004; Hale et al., 2007; Kim et al., 2007; Beck and Teboulle, 2009).实际中，我们可以根据具体问题选择合适的算法.
同样，当R确定后，(5)式可变为
|
(8) |
经过一些简单的代数运算，(8)式可变为
|
(9) |
这里
算法1总结了矩阵Toeplitz稀疏分解的整个求解过程:
算法1 矩阵Toeplitz稀疏分解
输入:观测矩阵Y∈ℝn×m，稀疏参数λ，Ψ(α)中的参数，初始子波矩阵A(0)∈ℝn×n.
输出: Toeplitz矩阵A和稀疏矩阵R.
S1: k=1;
S2: i=1;
S3:求解
S4: i=i+1;
S5:如果i≤m，返回S3，否则执行下一步S6;
S6: R(k)=(r1(k), …, rm(k))，
S7:求解
S8:根据公式(3)由α(k)得到A(k);
S9: k=k+1;
S10:如果收敛(或者计算次数达到给定的值)，执行下一步S11，否则返回S2;
S11: A=A(k), R=R(k).
2.2 基于矩阵Toeplitz稀疏分解的相对波阻抗反演令w(t)表示震源子波，r(t)表示反射系数序列，n(t)表示噪声，则反射地震道y(t)可用褶积模型表示为(Yilmaz, 1987)
|
(10) |
这里*表示褶积。(10)式相应的矩阵表示式为
|
(11) |
这里A表示具有Toeplitz结构的子波褶积矩阵，矩阵中的元素为Aij=wi-j+1，其中wi-j+1表示子波w(t)的第i-j+1个采样点。y, r和n分别表示y(t), r(t)和n(t)对应的离散采样列向量.对于一个二维的地震剖面，公式(11)变为
|
(12) |
这里Y表示地震剖面，R和N分别为相应的反射系数剖面和噪声剖面.实际中，我们很难观测到准确的子波，然而我们可以利用矩阵A的Toeplitz结构特性和矩阵R的稀疏特性来得到子波和反射系数.因此，上一节提出的矩阵Toeplitz稀疏分解算法可用来从地震剖面估计反射系数.
我们知道，地震子波通常是光滑的并且是紧支集的，其长度往往比地震数据短很多.也就是说，如果我们假设地震子波的长度为l个采样点长，则通常有l«n.这意味着在某种程度上α是稀疏的，并且非零值集中在矩阵A的对角线附近.因此，本文选择Fused Lasso (Tibshirani et al., 2005)惩罚项作为公式(5)的正则化项，即
|
(13) |
其中第一项使得α是稀疏的(即，使得α具有紧支集)，第二项和第三项使得α是光滑的.如果我们仅仅用第二项来作为α的光滑性约束项，估计出的子波可能会出现连续相等的值的情况，因此需要加入第三项来避免这种情况.为了使得参数的选取尽量简单，在算法迭代的过程中检测到子波出现连值时可用平均平滑来代替第三项.本文在后面的算例中就用平均平滑代替了第三项.将(13)式去掉第三项后代入(9)式可得
|
(14) |
(14)式有很多标准求解算法(Hoefling, 2010; Liu et al., 2010; Ye and Xie, 2011)，本文选用SPAMS (a SPArse Modeling Software)(Mairal, 2012)软件中的mexFistaFlat函数来求解该式.在(14)式中，参数β与地震道的长度和子波主频有关.实际中，我们先估计地震数据子波的主频，然后用具有相同地震道长度和子波主频的合成记录做实验来选取β.平均平滑可以减小(14)式中第三项的权重，当信噪比很高时，β1可取为0，当地震记录含有噪声时，β1可以根据情况选择0~1之间的值.
对于(7)式，本文用基于回溯的快速萎缩阈值迭代算法来求解(Beck and Teboulle, 2009)，该迭代算法用稀疏度K，也就是非零反射系数个数的上界代替参数λ.反射系数的稀疏度K可以基于地震道包络局部峰值的个数来估计.地震道包络局部峰值包含有较大的地震反射系数的信息，在每个峰值附近会有几个较大的反射系数.例如，对一道地震数据，我们用K0来表示地震道包络局部峰值的个数，然后以某一峰值附近最窄的波形作为单元去估计平均每个峰值周围可能有的反射系数的个数，并以稍大于这一平均个数的值作为c (c为常数)，则可取稀疏度K=cK0.可见，基于回溯的快速萎缩阈值迭代算法可将公式(7)中L1约束项参数λ的选择问题，转化为反射系数稀疏度K的估计问题，使得参数选择简便易行.基于回溯的快速萎缩阈值迭代算法如算法2所示.
算法2 基于回溯的快速萎缩阈值迭代算法
输入:观测地震道y∈ℝn×1, 稀疏度K, 初始子波矩阵A(0) ∈ℝ n×n, 初始反射系数序列r(0) ∈ℝ n×1, t(0)=1，z(0)=r(0).
输出:稀疏反射系数序列r ∈ℝ n×1.
S1: k=0，任选小常数ε∈(0, 1)设
S2: k=k+1，μ(k)=μ(0)；
S3: e=y-A(k-1)*r(k-1), u(k)=(A(k-1))Te,
b(k)=r(k-1)+μ(k)*u(k)；
S4:将b(k)按降序排列，令
S5: λμ=
S6:令JK=J(1:K), z(k)=zeros (n, 1),
|
S7: 
S8: 
S9:如果收敛(或者计算次数达到给定的值)，执行下一步S10，否则返回S2；
S10: r=r(k).
这一算法的基本原理在Beck和Teboulle (2009)和Xu等(2012)的文章中有论述.
本文用地震数据所有道振幅谱平均值的逆傅氏变换得到的零相位子波生成Toeplitz子波矩阵A的初值.对(14)式涉及的各参数进行了敏感性分析(Wang et al., 2016)，详见附录A，由此给出了表 1所示本文出现的所有参数的影响及其选择方法.
|
表 1 参数的影响及其选择方法 Table 1 Influence and selection method for the parameters |
此外，附录A中图A4给出的分析算例也表明，本文方法对于上述子波初值具有很好的容错度，在无噪情况下，即使初始子波相位误差达到-90°或接近90°，本文方法也都能很好地恢复子波.
当地震波垂直入射到两种介质的平面分界面上时，反射系数的表达式为(马劲风等，1999)
|
(15) |
其中Zi=ρivi为第i层的波阻抗，Zi+1=ρi+1vi+1为第i+1层的波阻抗.因此，如果知道了反射系数序列，就可以借助递推公式
|
(16) |
把反射系数转换为波阻抗.当已知第m(0≤m < n)层波阻抗Zm时，第n层波阻抗为
|
(17) |
同理，当已知第n层波阻抗Zn时，第m层波阻抗为
|
(18) |
如果我们可以得到某一标志层的波阻抗，根据(17)和(18)式就可推出其他层的波阻抗.
本文在反演得到反射系数矩阵后，设阻抗的初值Z1=1，根据(16)式递推反演可得到这种初值下的波阻抗.由于实际中各道地震数据的初始阻抗并不为1，因此本文设计了一个高通滤波器，滤掉递推的累计误差所造成的低频分量(8 Hz以下)，就可以得到具有横向连续性的相对波阻抗剖面.然后利用测井资料校正相对波阻抗振幅并加入低频分量后，得到高分辨率的绝对波阻抗剖面.
3 数值算例我们将用合成地震记录模型和实际地震资料验证方法的有效性.
3.1 模型算例大多数的稀疏脉冲反褶积方法需要预知子波(Nguyen, 2008).尽管目前有许多子波估计方法，其中有些方法可以估计非最小相位子波，但是子波估计的效果往往依赖于数据质量，不准确的子波将会影响稀疏脉冲反褶积的效果(Velis，2008).本文方法不需要预先知道准确的子波，而是用常用方法粗略预估一个子波作为初值，在算法迭代的过程中自动校正子波，对初始估计的子波有很好的容错度.
本算例所用的合成地震记录剖面由图 1所示的Marmousi2模型中矩形框圈出的区域生成，区域数据大小为650×200(每道采样点数×道数).为了简单起见，本文假设该区域的密度为常数1，每道有650个采样点，采样间隔为1 ms.即这块数据是一块有200道，道间隔为12.5 m，每道有650个采样点，采样间隔为1 ms的数据.由这块数据生成反射系数后与带有30°相位旋转的主频为30 Hz的Ricker子波褶积后，生成图 2a所示的无噪合成地震记录剖面.在图 2a的无噪合成地震记录剖面上添加标准差为地震道绝对值最大值的10%的高斯白噪声后，得到图 2b所示含噪合成地震记录剖面，计算可得该剖面的信噪比为2.2 dB.图 3a和图 3b分别给出了这块Marmousi2模型数据对应的相对阻抗和绝对阻抗.由于该块数据的密度设为常数1，图 3b中的绝对阻抗即为图 1b中的速度.图 3a中的相对阻抗是用过渡带为4~12 Hz的高通滤波器对图 3b中的绝对波阻抗做滤波后得到的.
|
图 1 本文合成地震记录算例所用的Marmousi2模型 (a) Marmousi2 P-波速度模型；(b)矩形框中区域放大图，用来生成合成地震记录算例中的反射系数、相对波阻抗和绝对波阻抗剖面. Fig. 1 Marmousi2 model for synthetic examples (a) Marmousi2 P-wave velocity model and (b) an enlarged portion in the rectangle area that will be used to generate the reflectivity profile, relative acoustic impedance profile, and acoustic impedance profile for the synthetic examples. |
|
图 2 Marmousi2模型生成的合成地震数据 (a)无噪合成地震数据；(b)含有噪声的合成地震数据，含有标准差为地震道绝对值最大值的10%的加性高斯白噪声，信噪比为2.2 dB. Fig. 2 The synthetic seismic profile from the Marmousi2 model (a) The noise-free synthetic data; (b) The noisy synthetic data with Gaussian noise approximately 10% of the maximum trace absolute value, the SNR is 2.2 dB. |
|
图 3 Marmousi2模型理论波阻抗 (a)相对波阻抗；(b)绝对波阻抗. Fig. 3 Acoustic impedance from Marmousi2 (a) Relative acoustic impedance; (b) Acoustic impedance. |
图 4a给出了从图 2a所示的无噪合成地震数据中估计出的相对波阻抗.在这个算例中，我们对图 2a中合成地震数据的所有地震道的振幅谱的平均值做逆傅氏变换，得到零相位子波，并用这个零相位子波生成Toeplitz子波矩阵A的初值.根据表 1给出的参数选择方法，本列选取参数K=40，β=6，β1=0.由算法1、公式(7)和(14)可反演得到反射系数剖面.然后取阻抗初值为1，根据反演递推公式(16)可以从反射系数剖面得到相对波阻抗剖面，再用过渡带为4~12 Hz的高通滤波器滤波并做振幅校正后可以得到图 4a所示的相对波阻抗剖面.对比图 4a和图 3a中的真实相对波阻抗剖面可见，本文方法估计得到的相对波阻抗剖面与真实相对波阻抗剖面差别不大，对地层纵向和横向波阻抗变化都刻画的比较细致.用过渡带为4~12 Hz的低通滤波器对图 3b中的绝对阻抗做滤波得到低频剖面，添加到图 4a所示的相对波阻抗剖面，可以得到图 4b所示的绝对波阻抗剖面.对比图 4b和图 3b中的真实绝对波阻抗剖面可见，本文方法得到的绝对波阻抗剖面与真实绝对波阻抗剖面差别不大，进一步证明了本方法的有效性.
|
图 4 从图 2a所示的无噪合成地震数据中估计的波阻抗 (a)相对波阻抗；(b)绝对波阻抗. Fig. 4 Acoustic impedance from the noise-free synthetic data in Fig 2a (a) Relative acoustic impedance; (b) Acoustic impedance. |
图 5给出了从图 2b所示的含噪合成地震数据中估计出的波阻抗.本算例也用图 2b中数据所有道振幅谱平均值的逆傅氏变换得到的零相位子波生成Toeplitz子波矩阵A的初值.根据表 1给出的参数选择方法，选取参数K=30，β=6，β1=0.01.然后同上述图 4a和4b的生成方法，可得到图 5a所示的相对波阻抗和图 5b所示的绝对波阻抗.对比图 5a和图 3a，以及图 5b和图 3b可见，即使地震数据含有较严重的随机噪声，本文方法估计得到的相对波阻抗，以及对其加低频后生成的绝对波阻抗仍然能够较好地刻画主要岩层展布.
|
图 5 从图 2b所示的含噪合成地震数据中估计的波阻抗 (a)相对波阻抗；(b)绝对波阻抗. Fig. 5 Acoustic impedance from the noisy synthetic data in Fig. 2b (a) Relative acoustic impedance; (b) Acoustic impedance. |
图 6给出了图 4无噪和图 5含噪情况下子波的对比.其中，点线表示初始估计得到的子波，是地震数据所有道振幅谱平均值做逆傅氏变换后得到的零相位子波；实线表示最终估计得到的子波，是将地震数据做Toeplitz稀疏分解得到的子波矩阵对应的子波；虚线对应的是真实子波，即带有30°相位旋转的主频为30 Hz的Ricker子波.由图 6a可见，虽然初始估计的子波并不准确，但是用本文Toeplitz稀疏分解方法反演估计得到的子波可以很好地逼近真实子波.由图 6b可见，噪声的影响使得初始估计的子波振幅和相位都存在误差，用本文方法反演得到的子波的振幅和相位都得到了很好地恢复.这充分说明本文方法对初始子波有很好的容错度，并且能够减少随机噪声对结果的影响.
|
图 6 子波对比 (a) 图 4无噪和(b)图 5含噪情况下，初始估计得到的子波(点线)，最终估计得到的子波(实线)和真实子波(虚线)的对比. Fig. 6 Wavelet comparison Comparison of the initial wavelet (dot line), the finally estimated wavelet (solid line), and the true wavelet (dashed line) corresponding to (a) the noise-free case in Fig. 1 and (b) the noisy case in Fig. 5. |
图 7抽取了无噪情况下从合成地震数据反演得到的相对波阻抗(图 4a)和真实相对波阻抗(图 3a)的第10道，第100道和第150道分别进行对比.图中红色实线对应反演得到的相对波阻抗，蓝色实线对应真实的相对波阻抗，三道对应的均方根误差分别为0.0161，0.0155，0.0163.可见，本文方法反演得到的相对波阻抗能够很好地逼近真实相对波阻抗.图 8抽取了无噪情况下得到的绝对波阻抗(图 4b，红线)和真实绝对波阻抗(图 3b，蓝线)的第10道，第100道和第150道分别进行对比.对比可见，得到波阻抗的低频信息以后，将它与反演得到的相对波阻抗结合可以很好地恢复绝对波阻抗.
|
图 7 从图 3a和图 4a所示的相对波阻抗剖面中抽3道对比(无噪.蓝线：真实；红线：估计) (a)第10道；(b)第100道；(c)第150道. Fig. 7 Three traces extracted from the relative acoustic impedance in Fig. 3a and Fig. 4a.(Corresponding to the noise-free data. Blue solid is the true impedance from Fig. 3a and red solid is the estimated impedance from Fig. 4a) (a) The 10th trace; (b) The 100th trace; (c) The 150th trace. |
|
图 8 从图 3b和图 4b所示的绝对波阻抗剖面中抽3道对比(无噪.蓝线：真实；红线：估计). (a)第10道；(b)第100道；(c)第150道. Fig. 8 Three traces extracted from the acoustic impedance in Fig. 3b and Fig. 4b.(noise-free, Blue: the true impedance from Fig. 3b, Red: the estimated impedance from Fig. 4b). (a) The 10th trace; (b) The 100th trace; (c) The 150th trace. |
图 9抽取了含噪情况下从合成地震数据中反演得到的相对波阻抗(图 5a)和真实相对波阻抗(图 3a)的第10道，第100道和第150道分别进行对比.图中红色实线对应反演得到的相对波阻抗，蓝色实线对应真实的相对波阻抗，三道对应的均方根误差分别为0.0175，0.0216，0.0228.对比可见，虽然地震数据含有较严重的随机噪声，本文方法反演得到的相对波阻抗仍然能够较好地描述主要岩层的阻抗特征.图 10抽取了含噪情况下得到的绝对波阻抗(图 5b，红线)和真实绝对波阻抗(图 3b，蓝线)的第10道，第100道和第150道分别进行对比.对比可见，在地震数据含有较严重的随机噪声的情况下，本文方法得到的绝对波阻抗仍然可以较好地刻画主要岩层展布.
|
图 9 从图 3a和图 5a所示的相对波阻抗剖面中抽3道对比(含噪.蓝线：真实；红线：估计) (a)第10道；(b)第100道；(c)第150道. Fig. 9 Three traces extracted from the relative acoustic impedance in Fig. 3a and Fig. 5a.(noisy. Blue: the true impedance from Fig. 3a, Red: the estimated impedance from Fig. 5a). (a) The 10th trace; (b) The 100th trace; (c) The 150th trace. |
|
图 10 从图 3b和图 5b所示的绝对波阻抗剖面中抽3道对比(含噪.蓝线：真实；红线：估计) (a)第10道；(b)第100道；(c)第150道. Fig. 10 Three traces extracted from the acoustic impedance in Fig. 3b and Fig. 5b.(noisy. Blue: the true impedance from Fig. 3b, Red: the estimated impedance from Fig. 5b) (a) The 10th trace; (b) The 100th trace; (c) The 150th trace. |
图 11给出了一个叠后实际地震资料算例.图 11a是一段中海石油研究中心提供的叠后地震数据，数据的大小为300×801(每道采样点数×道数)，时间采样间隔为2 ms.井1是该段地震数据中的一口验证井.图 11b给出了用Jason软件反演得到的相对波阻抗.与合成地震记录模型算例同理，我们仍采用由所有地震道振幅谱平均值的逆傅氏变换得到的零相位子波生成Toeplitz子波矩阵A的初值.该段地震数据的地震道包络局部峰值的平均个数为27.实际中地下介质层数通常比观测到的多很多，根据表 1给出的参数选择方法，本文选择稀疏度K=80，这一数值约为27的3倍.该地震剖面的主频为30 Hz，以合成地震数据算例(附录图A3)作参考，本例选择β=2.3.由于该段剖面含有噪声，因此选择β1=0.2来减少噪声的影响.将上述参数代入算法1、公式(7)和(14)可反演得到反射系数剖面.取阻抗初值为1，根据反演递推公式(16)可以从反射系数剖面得到相对波阻抗剖面，再用过渡带为4~12 Hz的高通滤波器滤波可以得到图 11c所示的相对波阻抗剖面.由图 11b和图 11c可见，本文方法反演得到的相对波阻抗剖面与Jason软件反演得到的相对波阻抗剖面相比，分辨率显著提高，且与测井资料得到的相对波阻抗有更好的对应关系.图 11d给出了用Jason软件反演得到的绝对波阻抗.图 11e给出了用Jason软件从测井和地震数据提取低频，然后与图 11c中的相对波阻抗合并得到的绝对波阻抗.对比图 11d和图 11e可见，本文方法得到的绝对波阻抗剖面对地层的刻画更加精细，并且与测井资料得到的绝对波阻抗有更好的对应关系.此外，本文方法反演得到的波阻抗剖面能够很好地保持横向连续性.
|
图 11 实际地震数据算例 (a)一段具有复杂地层结构的叠后实际地震数据；(b)用Jason软件反演得到的相对波阻抗；(c)用本文矩阵分解方法反演得到的相对波阻抗；(d)用Jason软件反演得到的绝对波阻抗；(e)从测井和地震数据提取低频，然后与(c)中的相对波阻抗合并得到的绝对波阻抗. Fig. 11 Field data example (a) A portion of poststack field data with complex structures; (b) The relative acoustic impedance obtained by Jason; (c) The relative acoustic impedance obtained by the proposed method; (d) The acoustic impedance obtained by Jason; (e) The high resolution acoustic impedance achieved by adding the low-frequency components derived from the well to the relative acoustic impedance in Fig. 11c. |
为了观察更加清晰，放大显示图 11b和11c矩形框中的区域，如图 12a和12b所示.由图可见，与Jason软件反演得到的相对波阻抗剖面相比，本文方法得到的相对波阻抗剖面分辨率更高，且与测井资料得到的相对波阻抗吻合更好.另外，图 12c和12d也给出了测井对应的岩性柱图，由图可见，Jason软件反演得到的相对波阻抗不能区分厚度分别为9.3 m和11.7 m的两个油层，而本文方法得到的波阻抗能够很好地区分这两个油层，如图中箭头指示区域所示.与测井对应的岩性柱图对比可见，本文方法将储层的可区分厚度从原来的20 m以上提升到6~10 m.
|
图 12 (a) 图 11b矩形框中区域放大图；(b) 图 11c矩形框中区域放大图；(c)在图 12a所示剖面中插入测井对应的岩性柱图；(d)在图 12b所示剖面中插入测井对应的岩性柱图 Fig. 12 (a) An enlarged portion in the rectangle area of Fig. 11b; (b) An enlarged portion in the rectangle area of Fig. 11c; (c) The profile in Fig. 12a with lithologic column; (d) The profile in Fig. 12b with lithologic column |
本文提出了一种基于矩阵Toeplitz稀疏分解的相对波阻抗反演方法.当地震记录满足褶积模型，并且反射系数满足稀疏性假设时，矩阵Toeplitz稀疏分解算法通过交替迭代求解两个子反演问题，其一以Toeplitz子波矩阵元素为待反演的参数，其二以稀疏反射系数矩阵元素为待反演参数，将地震剖面同时分解为一个子波矩阵和一个反射系数矩阵.Fused Lasso方法被用来求解Toeplitz子波矩阵中的元素，该方法既可保证模型系数的稀疏性(紧支集)，又使得邻近系数间更加平滑.基于回溯的快速萎缩阈值迭代算法被用来反演L1范数约束下的稀疏反射系数矩阵，从而将L1约束项的参数选择问题，转化为反射系数稀疏度的估计问题，使得参数选择简便易行.由反射系数矩阵递推反演可以得到高分辨率的相对波阻抗剖面，然后利用测井资料加入低频分量后，也可以得到高分辨率的绝对波阻抗剖面.矩阵分解方法的引入，使得反演过程有多道数据参与运算，能够有效减少随机噪声干扰影响.Marmousi2模型生成的合成记录算例和实际地震资料算例均表明：本文方法可以从带限地震数据中有效地反演出高分辨率的相对波阻抗，即使地震资料被噪声污染也能取得较好的效果，并且该方法对初始子波有很好的容错度，在处理具有复杂结构的地震数据时，得到的阻抗反演结果仍然能够很好地保持横向连续性.与Jason软件反演得到的相对波阻抗和绝对波阻抗剖面对比表明，本文方法得到的波阻抗剖面对地层的刻画更加精细，与测井资料得到的波阻抗有更好的对应关系，可将储层的可区分厚度从原来的20 m以上提升到6~10 m之间.
需要指出的是本文方法中的矩阵运算比较耗内存.对于大规模的三维地震数据，可以先在数据中选取一些地震数据块，将各个数据块展开成二维剖面用本方法处理，然后将所有估计得到的子波平均，再用算法2计算整个三维数据体的反射系数，并递推反演得到三维数据体的波阻抗.
附录A为了方便阅读，将参数敏感性分析算例(Wang et al., 2016)复述如下：
图A1给出了无噪情况下稀疏度K的敏感性分析算例.图A1a为真实反射系数序列，反射系数的个数为22.图A1b依次给出了稀疏度K等于22，25，30，40以及50时估计得到的反射系数序列，这里取β=6，β1=0.从这一算例，我们可以看到，当K增大到一定程度时，本文方法会产生一些小的假反射系数.尽管如此，当K的值达到真实反射系数个数的2倍以上时，本文方法仍旧能够取得好的效果.反射系数的稀疏度K可以基于地震道包络局部峰值的个数来估计.地震道包络局部峰值包含有较大的地震反射系数的信息，在每个峰值附近会有几个较大的反射系数.假设地震道包络局部峰值的个数为K0，估计平均每个峰值周围可能有的反射系数的个数c(通常c可取1.5~3之间的数)，则可设K=cK0.
|
图 A1 稀疏度K的敏感性分析 (a)真实的反射系数序列; (b)从上到下依次为K等于22, 25, 30, 40和50时估计得到的反射系数，这里取β=6，β1=0. Fig. A1 Sensitivity analysis for the choice of sparsity K (a) True reflectivity; (b) Estimated reflectivity when K is 22, 25, 30, 40 and 50 from top to bottom, β=6 and β1=0 in all cases. |
图A2给出了参数β的敏感性分析算例.图A2b给出了β依次等于0.01，2，6，10，20，50以及100时估计得到的子波和反射系数，这里K=30，β1=0.从图中可以看到，随着β的增加，子波的支集变的更加紧凑.当β取为2~20之间的值时，本文方法都能取得好的效果，这意味着β的取值范围较宽.
|
图 A2 参数β的敏感性分析 (a)真实子波和真实反射系数序列; (b)从上到下依次为β等于0.01, 2, 6, 10, 20, 50和100时估计得到的子波和反射系数，这里取K=30，β1=0. Fig. A2 Sensitivity analysis for the choice of β (a) True wavelet and true reflectivity; (b) Estimated wavelet and estimated reflectivity when β is 0.01, 2, 6, 10, 20, 50 and 100 from top to bottom, K=30 and β1=0 in all cases. |
图A3给出了参数β与地震道长度和子波主频关系分析算例.其中地震道的长度从图A2中的800个采样点缩减为前300个采样点.图A3b给出了β依次等于1，2.3，5和10时估计得到的子波和反射系数，这里子波仍为主频30 Hz的Ricker子波，K=10，β1=0.从图中可以看到，当β等于2.3和5时，本文方法能够取得好的效果，但当β增加到10时，方法的效果开始变差.原因在于，当子波主频不变，而地震数据长度变短时，子波矩阵A的稀疏度变弱，因此β的值相应变小.与图A2中的算例相比，本例中β可取为图A2中β的3/8(即数据长度之比).图A3c给出了β依次等于6，15，20以及60，子波主频为60 Hz的Ricker子波，K=10，β1=0时估计得到的子波和反射系数.图中可以看到，这里β的取值比图A3b大.这是因为当子波的主频增加，而地震数据长度不变时，子波矩阵A的稀疏度增强，因此β的值相应变大.对比真实反射系数和估计得到的反射系数可见，β可由地震数据的长度，子波的主频来确定.实际中，我们可从与实际数据具有相同长度和相同子波主频的合成数据实验中给出β的参考值.
|
图 A3 参数β与地震数据长度和子波主频关系分析 与上一算例相比，本算例中地震数据的长度减为300个采样点.(a)真实的反射系数; (b)从上到下依次为β等于1, 2.3, 5和10时估计得到的子波和反射系数，这里子波依旧为主频为30 Hz的Ricker子波; (c)从上到下依次为β等于6, 15, 20和60时估计得到的子波和反射系数，这里子波变为主频为60 Hz的Ricker子波.本算例中取K=10，β1=0. Fig. A3 Analysis of the relation of β to the length of the seismic trace and the dominant frequency of the wavelet Comparing with the above examples, the length of the seismic trace is reduced to 300 sampling points in this example. (a) True reflectivity; (b) Estimated wavelet and estimated reflectivity when β is 1, 2.3, 5 and 10 from top to bottom, where the wavelet is still a 30 Hz Ricker wavelet; (c) Estimated wavelet and estimated reflectivity when β is 6, 15, 20 and 60 from top to bottom, where the wavelet is changed to a 60 Hz Ricker wavelet. Here, K=10 and β1=0 in all cases. |
图A4给出了本文方法对地震数据初始子波相位误差容错度的分析算例.本文用地震数据所有道振幅谱平均值的逆傅氏变换得到的子波生成Toeplitz子波矩阵A的初值，也就是说初始子波是零相位.图A3b给出了真实地震子波的相位依次从-90°变化到90°时估计得到的子波和反射系数，这里K=30，β=6，β1=0.从图中可以看到，除了真实子波相位为90°时，估计得到的子波和反射系数的极性与真实情况相反外，其他相位情况下，本文方法都能很好的估计出子波和反射系数.因此，在实际中，我们首先需要知道地震数据的极性，然后赋予初始零相位子波正确的极性才能取得好的处理效果.
|
图 A4 初始子波相位容错度分析 (a)真实反射系数序列; (b)估计子波(实线)，真实子波(虚线)和估计得到的反射系数序列，数据相位从-90°变到90°，这里取K=30, β=6, β1=0. Fig. A4 Analysis of the tolerance to phase error of the initial wavelet (a) True reflectivity; (b) estimated wavelet (solid line), true wavelet (dash line), and estimated reflectivity, data phase changes in the range from -90° to 90°, K=30, β=6, and β1=0 in all cases. |
| Beck A, Teboulle M. 2009. A fast iterative shrinkage-thresholding algorithm for linear inverse problems. SIAM Journal on Imaging Sciences, 2(1): 183-202. DOI:10.1137/080716542 | |
| Cui Y, Wang Y F, Yang C C. 2009. Regularizing method with a priori knowledge for seismic impedance inversion. Chinese J. Geophys.(in Chinese), 52(8): 2135-2141. DOI:10.3969/j.issn.0001-5733.2009.08.023 | |
| Debeye H W J, van Riel P. 1990. Lp-norm deconvolution. Geophysical Prospecting, 38(4): 381-403. DOI:10.1111/j.1365-2478.1990.tb01852.x | |
| Di H B, Guo Y Q, Liu X W. 2011. Relative acoustic impedance inversion based on Cauchy sparseness constraint Bayesian estimation. Geophysical Prospecting for Petroleum (in Chinese), 50(2): 124-128. DOI:10.3969/j.issn.1000-1441.2011.02.003 | |
| Efron B, Hastie T, Johnstone I, et al. 2004. Least angle regression. The Annals of Statistics, 32(2): 407-499. DOI:10.1214/009053604000000067 | |
| Ge D D, Jiang X Y, Ye Y Y. 2011. A note on the complexity of Lp minimization. Math. Program., Ser. B, 129(2): 285-299. DOI:10.1007/s10107-011-0470-2 | |
| Guitton A, Claerbout J. 2015. Nonminimum phase deconvolution in the log domain:A sparse inversion approach. Geophysics, 80(6): WD11-WD18. DOI:10.1190/geo2015-0016.1 | |
| Hale E T, Yin W, Zhang Y. 2007. A fixed-point continuation method for l1-regularized minimization with applications to compressed sensing. CAAM Technical Report 07-07, Rice University. | |
| Hoefling H. 2010. A path algorithm for the fused lasso signal approximator. Journal of Computational and Graphical Statistics, 19(4): 984-1006. DOI:10.1198/jcgs.2010.09208 | |
| Kaaresen K F, Taxt T. 1998. Multichannel blind deconvolution of seismic signals. Geophysics, 63(6): 2093-2107. DOI:10.1190/1.1444503 | |
| Kazemi N, Sacchi M D. 2014. Sparse multichannel blind deconvolution. Geophysics, 79(5): V143-V152. DOI:10.1190/geo2013-0465.1 | |
| Kim S J, Koh K, Lustig M, et al. 2007. An interior-point method for large-scale l1-regularized least squares. IEEE Journal of Selected Topics in Signal Processing, 1(4): 606-617. DOI:10.1109/JSTSP.2007.910971 | |
| Kormylo J, Mendel J M. 1978. On maximum-likelihood detection and estimation of reflection coefficients.//48th Ann. Internat. Mtg., Soc. Expi. Geophys.. Expanded Abstracts, 45-46. | |
| Latimer R B, Davidson R, van Riel P. 2000. An interpreter's guide to understanding and working with seismic-derived acoustic impedance data. The Leading Edge, 19(3): 242-256. DOI:10.1190/1.1438580 | |
| Liu J, Yuan L, Ye J P. 2010. An efficient algorithm for a class of fused lasso problems.//Proceedings of the 16th ACM SIGKDD International Conference on Knowledge Discovery and Data Mining. Washington, DC, USA:ACM, 323-332. | |
| Lu W K. 2009. Frequency recovery of band-limited seismic data based on sparse spike train deconvolution and lateral coherence constraint.//Beijing 2009 International Geophysical Conference and Exposition. Beijing, China, 117, doi:10.1190/1.3603642. | |
| Ma J F, Wang X F, Zhao S L, et al. 1999. Choice of initial wave impedance value and the multisolution of seismic trace inversion. Oil Geophysical Prospecting (in Chinese), 34(3): 332-340. DOI:10.13810/j.cnki.issn.1000-7210.1999.03.012 | |
| Ma J, Wang X, Zhao S, et al. 1999. Choice of initial wave impedance value and the multisolution of seismic trace inversion. OGP, 34(3): 332-340. DOI:10.13810/j.cnki.issn.1000-7210.1999.03.012 | |
| Mairal J. 2012. SPAMS:a SPArse Modeling Software, v2.3, http://spams-devel.gforge.inria.fr/index.html, released May 23, 2012. | |
| Nguyen T H. 2008. High resolution seismic reflectivity inversion[Ph. D. thesis]. Houston:University of Houston. https://www.researchgate.net/publication/288777891_High-resolution_reflectivity_inversion | |
| Nose-Filho K, Takahata A K, Lopes R, et al. 2016. A fast algorithm for sparse multichannel blind deconvolution. Geophysics, 81(1): V7-V16. DOI:10.1190/geo2015-0069.1 | |
| Oldenburg D W, Scheuer T, Levy S. 1983. Recovery of the acoustic impedance from reflection seismograms. Geophysics, 48(10): 1318-1337. DOI:10.1190/1.1441413 | |
| Sacchi M D, Velis D R, Cominguez A H. 1994. Minimum entropy deconvolution with frequency-domain constraints. Geophysics, 59(6): 938-945. DOI:10.1190/1.1443653 | |
| Tibshirani R, Saunders M, Rosset S, et al. 2005. Sparsity and smoothness via the fused lasso. Journal of the Royal Statistical Society, Series B (Statistical Methodology), 67(1): 91-108. DOI:10.1111/rssb.2005.67.issue-1 | |
| Tibshirani R. 1996. Regression shrinkage and selection via the lasso. Journal of the Royal Statistical Society, Series B (Methodological), 58(1): 267-288. | |
| Tseng P. 2001. Convergence of a block coordinate descent method for nondifferentiable minimization. Journal of Optimization Theory and Applications, 109(3): 475-494. DOI:10.1023/A:1017501703105 | |
| Velis D R. 2008. Stochastic sparse-spike deconvolution. Geophysics, 73(1): R1-R9. DOI:10.1190/1.2790584 | |
| Wang J F, Wang X S, Perz M. 2006. Structure preserving regularization for sparse deconvolution.//76th Ann. Internat. Mtg., Soc. Expi. Geophys., Expanded Abstracts, 2072-2076. | |
| Wang L L, Zhao Q, Gao J H, et al. 2016. Seismic sparse-spike deconvolution via Toeplitz-sparse matrix factorization. Geophysics, 81(3): V169-V182. DOI:10.1190/GEO2015-0151.1 | |
| Wang Y F. 2011. Seismic impedance inversion using l1-norm regularization and gradient descent methods. Journal of Inverse and Ill-Posed Problems, 18(7): 823-838. DOI:10.1515/jiip.2011.005 | |
| Xu Y Y, Yin W. 2013. A block coordinate descent method for regularized multiconvex optimization with applications to nonnegative tensor factorization and completion. SIAM Journal of Imaging Sciences, 6(3): 1758-1789. DOI:10.1137/120887795 | |
| Xu Z B, Chang X Y, Xu F M, et al. 2012. L1/2 regularization:A thresholding representation theory and a fast solver. IEEE Transactions on Neural Networks and Learning Systems, 23(7): 1013-1027. DOI:10.1109/TNNLS.2012.2197412 | |
| Yilmaz Ö. 1987. Seismic data processing. Tulsa Oklahoma:Society of Exploration Geophysicists. | |
| Ye G B, Xie X H. 2011. Split bregman method for large scale fused lasso. Computational Statistics & Data Analysis, 55(4): 1552-1569. DOI:10.1016/j.csda.2010.10.021 | |
| Yuan S Y, Wang S X, Luo C M, et al. 2015. Simultaneous multitrace impedance inversion with transform-domain sparsity promotion. Geophysics, 80(2): R71-R80. DOI:10.1190/GEO2014-0065.1 | |
| Zhang Y G. 2002. The present and future of wave impedance inversion technique. Geophysical Prospecting for Petroleum (in Chinese), 41(4): 385-390. | |
| 崔岩, 王彦飞, 杨长春. 2009. 带先验知识的波阻抗反演正则化方法研究. 地球物理学报, 52(8): 2135–2141. DOI:10.3969/j.issn.0001-5733.2009.08.023 | |
| 邸海滨, 郭玉清, 刘喜武. 2011. 基于稀疏约束贝叶斯估计的相对波阻抗反演. 石油物探, 50(2): 124–128. DOI:10.3969/j.issn.1000-1441.2011.02.003 | |
| 马劲风, 王学军, 赵圣亮, 等. 1999. 波阻抗初值选择与地震道反演的多解性. 石油地球物理勘探, 34(3): 332–340. DOI:10.13810/j.cnki.issn.1000-7210.1999.03.012 | |
| 张永刚. 2002. 地震波阻抗反演技术的现状和发展. 石油物探, 41(4): 385–390. | |