2019, Vol. 62
2. 中石化上海海洋油气分公司, 上海 200120;
3. 中石化石油物探技术研究院, 南京 210014
2. Sinopec Shanghai Offshore Oil and Gas Company, Shanghai 200120, China;
3. Sinopec Geophysical Research Institute, Nanjing 210014, China
适用于偏移成像的宏观速度模型建立是反射地震成像处理中的一个关键步骤，基于层析反演的偏移速度建模已经进入了工业界的标准偏移速度建模流程.不同于传统的成像域层析成像方法，立体层析是一种数据域层析成像方法；也不同于传统的基于走时的反射层析算法，该方法充分利用了地震波场中所有的运动学信息.它将局部相干同相轴在炮道集与检波点道集内的射线参数(以下简称P参数)水平分量、炮点坐标以及检波点坐标都纳入到层析反演的数据空间之中，重新定义了层析反演的模型空间和数据空间，使得立体层析反演成为射线类层析方法中唯一一种可以同时反演速度、反射点位置与局部地层倾角的方法(Billette and Lambaré, 1998).国内外研究者在关于立体层析的反演理论与应用实践方面已经有大量工作(Billette et al., 2003；Alerin et al., 2007; 倪瑶等，2013；Prieux et al., 2013；李振伟等，2014；任丽娟等，2014；王宇翔等，2016；杨锴等，2016；邵炜栋等，2017).近年来，基于射线类层析方法建立具有地质一致性特征的块状速度模型成为一个研究热点，虽然这个目标听起来似乎超出了射线类层析方法的能力，但是学术界和工业界一直在专注于这方面的工作，也获得了相当多的成果(Farra and Madariaga, 1988; Zelt and Smith, 1992; Kosloff et al., 1996; Clapp et al., 2004; Zhou, 2006).
对于工业界用户来说，最经典的获得地质一致性块状模型的方式当属层剥离方法.该方法从20世纪90年代初即进入商业软件，至今仍为工业界广泛应用.但是层剥离方法的问题是浅层的速度和界面位置一旦确定就无法再改动，误差会从浅层慢慢积累至深层，导致深层的速度建模产生大的偏差.因此有学者开始探索是否可以在射线类层析方法中同时实现界面和速度的更新.近年来在射线类层析的研究进展方面，有两类获得地质一致性块状模型的方法值得关注.第一种方法是使用混合正则化项，比如在地质界面处使用L1范数或其他稀疏类范数使得反演结果保留模型的突变特征，而在其他地方则使用L2范数或其他光滑类范数保持反演结果的平滑(Zhang and Zhang, 2012).第二种方法是Lambaré等(2013)提出的多层斜率层析，作者声称引入了一种混合速度格式，这种混合速度格式可以同时定义块内的光滑速度信息和地质层位信息，并在迭代中结合运动学偏移来更新地质层位.该方法同时提供高分辨率层内反演结果和强速度反差的地质界面，实际应用中表现非常稳健.然而这种混合速度格式的细节从未正式发表.邵炜栋等(2017)基于三角网格剖分发展了一种新的立体层析方法，理论数据测试表明该方法能够获得地质一致性块状模型，对实际数据的应用研究仍在进行之中.
基于射线类层析方法获得地质一致性块状结构速度的目标非常明确：反演得到的速度模型应在保证地质块内足够光滑的同时尽可能维持边界的突变.对于立体层析而言，这个目标可以通过更简洁的方法来实现.本文提出一种实用的模型正则化方法，该方法的核心思想是将封闭后的地层格架信息融入到立体层析Fréchet导数矩阵中，在通过融入地层格架信息简化了立体层析Fréchet导数矩阵之后，使得更新后的速度模型自动呈现出与地质构造较为一致的的层状或者块状结果.其实现过程如下：(1)通过非均匀B样条方法对更新后的反射点位置进行自动拟合得到更新后地质层位信息；(2)利用地质块体之间的拓扑关系对地质层位信息实施自动封闭得到封闭后的地层格架信息；(3)利用封闭后的地层格架信息对Fréchet导数矩阵进行分块，每块内的速度更新量以及速度梯度更新量均为常数，由于一个模型内的地质块数目一般不会很多，这样处理的结果使得模型空间中的速度分量数目被大幅度压缩，并且在更新过程中模型自动呈现出层状或者块状地质特征.
设计上述实现策略的原因是：反射点位置本身即属于立体层析模型空间的一部分，反射点位置的更新在立体层析反演的迭代过程中是可以自动得到的.相比在每一轮反演后实施偏移成像再进行人工拾取，这是一种不仅效率更高、而且更为合理的方法.同时相比于Lambaré等(2013)建议的利用运动学偏移来更新地质层位的方式，本文提出的方法在实践中表现也更为稳定.由于每个地质块体内的速度将设为线性变化甚至是常速，因此模型空间中关于速度信息的变量个数显然大大减少了，这对于层析线性方程组的求解无论是精度和效率都极有益处.
本文第1节首先介绍传统立体层析方法的数据空间和模型空间以及本文提出的地层格架约束下立体层析的数据空间和模型空间；第2节给出了地层格架正则化下的立体层析Fréchet导数矩阵的具体算法；第3节给出了基于散乱反射点的非均匀B样条拟合算法与地层格架自动封闭的实现过程；第4节给出一个实用的工作流程；第5节给出了理论与实际数据应用算例.作者期望这种特殊的正则化方式获得的立体层析反演结果将可以用于逆时偏移处理或作为全波形反演的高质量初始模型.
1 地层格架约束下的二维立体层析反演 1.1 二维立体层析的数据空间与模型空间及立体层析Fréchet导数矩阵的建立过程
首先介绍常规二维立体层析的模型空间、数据空间及Fréchet导数矩阵矩阵的建立过程.图 1a显示了一根从炮点S出发、到地下反射点X反射、回到地表R的射线.图 1b显示了射线中心坐标系中射线起点和端点扰动的几何关系.从透射的角度，不妨将其理解为从X出发、分别以入射角θs与反射角θr发射至炮点和检波点的两根透射射线.当速度正确时，两根透射射线分别以正确的出射方向、正确的走时达到正确的炮点和检波点位置.这样构成的立体层析数据空间各分量为：d=[Sx, Rx, T, Psx, Prx].其中Sx、Rx为炮检点坐标；T为走时；Psx、Prx为炮检点处的射线参数水平分量.地下的模型空间m =(X0, Z0, θs, θr, V)，其中X0、Z0为地下反射点X的坐标；θs、θr分别代表从炮点一侧与检波点一侧出射的地下张角, V代表地下介质速度信息.注意Psx可以在共检波点道集内搜索同相轴的局部倾角得到，Prx可以在共炮点道集内搜索同相轴的局部倾角得到.假设Pslope为局部相干同相轴的斜率，根据几何关系以及射线理论中慢度矢量的定义得：
|
图 1 二维立体层析数据空间各分量与模型空间各分量 (a)立体层析的模型与数据分量；(b)射线中心坐标系下射线起点和端点扰动的几何关系 Fig. 1 Model and data components in 2D stereo tomography data space (a) The model components and the data components in 2D stereotomography data space; (b) Geometrical relation between perturbations of ray starting and end points in ray-centered coordinates |
任何一种层析反演方法都必须建立数据空间扰动与模型空间扰动的线性关系，即：
|
(1) |
这样利用观测到的数据残差就可以计算出模型空间的扰动量，达到更新初始模型的目的.其中关键是建立方程(1)中所示的F矩阵，即Fréchet偏导数矩阵公式.Fréchet偏导数矩阵的物理意义是数据空间任一分量关于模型空间任一分量的敏感度.在常规走时层析中，走时t关于速度V的偏导数其实就是射线经过每一个速度单元的弧长.(2)式展示了二维常规立体层析Fréchet偏导数矩阵，除了第一行走时关于模型空间的偏导数根据走时积分方程即可得到外.其余元素是都通过射线扰动理论(Farra and Madariaga, 1987)计算得到的，其推导细节请查阅Billette和Lambaré(1998)或杨锴等(2016)，这里不再赘述.其中σ为针对矩阵中每一行的数量级不同设置的均衡系数.可以看出立体层析矩阵无论从规模还是其稀疏性方面都远远超过了常规的走时层析矩阵.如果能够通过某种有效的正则化手段压缩模型空间、改善层析矩阵的条件数，降低求解时对光滑规则化的要求，对于获得地质一致性的模型无疑是非常有益的.公式(2)为
|
(2) |
在传统立体层析中速度模型一般由B样条基函数的权系数来描述，这种表达可以达到压缩立体层析Fréchet导数矩阵规模的目的.然而，B样条描述方式的一个副作用是会使得反演结果过于平滑.与此同时，传统的正则化技术一般都建议加强解的平滑性(Costa et al., 2008)，这两个因素叠加在一起会使得反演结果愈发平滑.与上述模型描述方式和正则化方式不同，本文提出的地层格架正则化约束下的立体层析试图强化反演结果在地质界面处的突变特征.这个目的是通过将一个闭合的二维地质格架信息作用于Fréchet导数矩阵，使得更新结果体现出地质一致性块状特征而实现的.在建立地层格架信息并对其实施封闭之后，每个地质块中的速度被设置为是关于横坐标x和纵坐标z的线性变量，在某些简单情况下(比如海底之上的水层速度)甚至可以设为常速.将地层格架信息作用于Fréchet导数矩阵相当于施加了一个很强的正则化手段，其优点是在反演时每一地质块内的速度更新量和速度变化率的更新量都将是一个常数，更新之后的模型将自动呈现出具有地质一致性的块状特征.从某种程度上，它也可以理解为对模型实施了重新参数化，使得模型以用更为稀疏的方式得到表达(Harlan, 1995).
通过一个概念实验即可清楚阐述地层格架约束立体层析的核心思想.图 2显示了一个6层背斜模型.如果按照均匀B样条基函数表达这个模型，在节点纵横向间距都为400 m(一般取200 m或400 m)时也需要480个B样条系数，也就是说模型空间中关于速度的未知数有480个.如果说用常规的规则网格表达图 2所示的速度模型，在网格纵横向间隔都设为100 m时，会有120×48=5760个网格，也就是说模型空间中关于速度的未知数将有5760个.然而，如果我们假定每一块内的速度都应该是常数时，其模型空间内的速度未知数其实应该是6个.基于一个由6个地质界面封闭得到的地层格架实施上述规则化就相当于将速度未知量的数目从一个很大的规模压缩为6个.同时，由于每一块内的速度更新量各有差异，因此更新后的模型在地质层位两侧必然会呈现出强反差.这正是地层格架约束立体层析所想要得到的地质一致性块状特征.
|
图 2 地层格架约束立体层析概念实验 Fig. 2 The concept experiment of stereo tomography constrained by stratum framework |
一个关键问题是：地层格架如何在反演中的到正确的更新？注意地层格架是由许多地质层位封闭成的，而这些地质层位是由许多反射点构建的.而这些反射点的位置本身就是立体层析模型空间的一部分，在每一次迭代中它们都和速度、散射角等模型信息一起被更新.因此完全可以利用更新后的反射点位置来拟合出更新后的地质层位.这里我们使用非均匀B样条(NUBS)(Piegl and Tiller, 1997)来实施拟合.NUBS是一种发展完备的数值拟合算法，在计算机辅助设计(CAD)、计算机辅助制造(CAM)领域中得到广泛应用.它所特有的高效灵活以及能够应付多值函数的特点非常适于将反射点拟合为连续界面.在反射点拟合为连续层位之后，对地质块进行实施封闭则需要利用图像处理中的结构张量算法(Van Vliet and Verbeek, 1995)，这是一种发展完备的图像边缘检测算法.在利用NUBS拟合出连续的地质层位后，再结合结构张量算法提供的界面梯度信息，就可以快速实现地层格架的自动封闭.基于反射点位的地层更新以及地层格架的自动封闭是实现本文地层格架约束立体层析反演算法的一个重要前提条件.
2 二维地层格架约束立体层析Fréchet导数矩阵求取那么二维地层格架约束立体层析的Fréchet导数如何计算？注意二维地层格架约束立体层析的数据空间与传统二维立体层析并无差别，依然可以描述为d=(Sx, Rx, T, Psx, Prx).但是由于模型空间被划分为若干个地质块，而且在每块中其速度可以设为v=v0+kxx+kzz，这里v0为参考速度，kx和kz分别用来描述速度在水平和垂直方向上的线性变化率，其模型空间变为m =(x0, z0, θs, θR, v0, Kx, Kz).因此二维地层格架约束立体层析的Fréchet导数矩阵可以写为
|
(3) |
炮点一侧的数据空间对于检波点一侧的散射角显然没有敏感度，反之亦然.同时根据费马原理，走时T在速度、炮点与接收点坐标都固定时已经是最小走时，因此散射角对走时亦不存在一阶扰动，有：
|
根据Farra和Madariaga(1987)，二维中心坐标系下立体层析所用的线性扰动关系如(4)、(5)式所示.(4)式等号左边为地表观测的坐标扰动(Δq1)和射线参数扰动(Δp1)，其实就是地表观测的立体层析数据空间(注意不包括走时t)；(4)式等号右边第一项是地下初始相空间扰动与传播矩阵П相乘，代表了当射线出射点坐标发生扰动(Δq0)和射线出射方向发生扰动(Δp0)时对立体层析数据空间的一阶效应；(4)式等号右边第二项代表当地下速度参数发生扰动(ΔW)时对立体层析数据空间造成的一阶效应.(5)式代表地表观测的走时扰动和地下模型空间参数扰动的线性关系，这里以炮点一侧为例，S0代表地下射线出射点坐标，S1代表地表观测点坐标，θ代表炮点一侧散射角，各参数具体含义如图 1b所示.公式(4)、(5)为
|
(4) |
|
(5) |
(4) 式中的Δq0、Δp0代表二维中心坐标系下的初始坐标扰动和初始慢度扰动，Δq1、Δp1代表二维中心坐标系下射线每一点上的坐标扰动和射线参数扰动(几何关系如图 1b所示)，П代表传播矩阵，传播矩阵中的四个元素Q1、Q2、P1、P2可通过二维动力学射线追踪得到.这里的传播矩阵 П、П-1以及Δw算子定义为
|
(6) |
|
(7) |
将(7)式代入(4)式容易推出：
|
(8) |
方括号中后两项为
|
(9) |
显然，在射线到达地表后，需要将扰动量从中心坐标系换算到地表水平观测面.因此需要建立射线中心坐标系下的位移扰动Δq0与水平观测面上的位移扰动Δξ之间、射线初始起飞角扰动Δθ和初始慢度扰动Δp0之间、射线中心坐标系下初始位移扰动Δq0与直角坐标系下两个位移扰动Δx0、Δz0之间、射线中心坐标系下慢度扰动Δp1与水平观测面慢度扰动Δpξ之间的几何关系(如图 1b所示).这些关系表达式为
|
(10) |
将(6)、(7)、(8)、(9)、(10)式代入到(4)式中，可得出二维中心坐标系下地表观测数据信息与初始射线扰动扰动信息的一阶关系为
|
(11) |
|
(12) |
如前所述，当速度模型被地层格架信息分块之后，速度更新量在一个地质块内将是一个常数，因此(12)式中括号里的第一项，即速度更新量关于空间坐标的导数就应该是零；其次一个地质块内的速度定义为
|
(13) |
根据链式法则，二维中心坐标系下坐标扰动Δq、慢度扰动Δp关于v0, kx, kz的一阶导数可以推导为
|
(14) |
最终我们得到：
|
(15) |
类似地，将链式法则应用于(5)式，就可以得到走时关于v0、kx、kz的一阶偏导数为
|
(16) |
至此，在引入地层格架信息约束后，走时t, 地表观测坐标x, 地表观测水平慢度矢量px关于v0、kx、kz的偏导数已经全部得到，公式(3)所示Fréchet导数矩阵中的每一个元素都可以计算得到.
3 基于散乱反射点的非均匀B样条拟合与地层格架自动封闭对地层格架约束立体层析而言，在反演实施之前对散乱反射点的拟合和自动封闭至关重要.这里我们采用了非均匀B样条拟合算法将散乱反射点拟合为连续的地质层位.相比均匀B样条拟合，该算法能够高效拟合多值函数，其稳健性在其他行业中已经得到了充分体现，证明完全能够适应二维、三维复杂地质建模的需求.这里对其实现原理做一简介，更多细节请读者参阅(Piegl and Tiller, 1997).
3.1 非均匀B样条拟合B样条基函数的定义如下：令U={u0, u1, …, um}是一个单调不减的实数序列，即ui≤ui+1, i=0, 1, …, m－1.其中ui称为节点，U称为节点矢量，用Ni, p(u)表示第i个p次B样条基函数，其定义为
|
(17) |
而对于p次B样条曲线的定义则为
对于数据拟合，需要预先计算好数据点的参数值和节点矢量，然后，建立并求解线性最小二乘问题来求解未知控制点.假定p≥1, n≥p，并且给定数据点Q0，…, Qr(r>n).我们试图需求一条p次非均匀B样条曲线，公式为
|
(18) |
其满足条件：Q0=C(0), Qr=C(1).
求取拟合曲线主要分为以下几个步骤：
(1) 预先通过弦长参数化得到数据点Qk对应的参数值uk.令u0=0, ur=1，则：
|
(19) |
由于节点的分布应该反映节点矢量的分布，假定c是一个正实数，用i=int(c)表示小于或等与d的最大整数，则总共需要n+p+2个节点，因此有n-p个内节点和n-p-1个内部节点区间，令c=(m+1)/(n-p+1)，然后按公式(20)定义内节点，即：
|
(20) |
(2) 建立目标函数，公式为
|
(21) |
(3) 根据目标函数求法方程，进行求解.令f关于n-1个未知控制点的偏导都等于零，则得到含n-1个未知量和n-1个方程的线性方程组，公式为
|
(22) |
这里，N是由标量组成的(r-1)×(n-1)的矩阵，R′是由(n-1)个点组成的列向量.通过对公式(22)通过进行求解可以得到所求控制点的值.
(4) 将所求的控制点值代入(18)式中，便可得到根据数据点拟合的曲线.图 3b显示了使用3阶非均匀B样条基函数的插值效果.
|
图 3 B样条基函数以及拟合结果 (a) 3阶B样条基函数；(b) 10个数据点(红点)的非均匀B样条拟合曲线(蓝线)；(c)中为大量数据点(红点)的B样条拟合曲线(蓝线)；(d)一个背斜顶面的拟合结果(蓝线) Fig. 3 B-Spline basis functions and curve fitting (a) Nonzero third-degree basis functions; (b) Non-uniform B-Spline curve fitting (blue line) over 10 data points (red dots); (c) Curve fitting (blue line) over a lot of data points (red dots); (d) Curve fitting of an anticline |
第3.1节中的非均匀B样条拟合仅仅是针对某一个地质层位的离散反射点完成了连续地质层位的拟合.但是要实施地层格架的自动封闭必须考虑地质层位之间的相交或者令其自然延伸到剖面的边界处，这些在实施自动封闭时都需要非常仔细的考虑.此外，对于含有多个地质块体的地震剖面，地层的封闭不但需要考虑块体顶底面之间的关系，还要考虑块体与地层以及地层与地层之间的关系.这要求离散的反射点数据必须含有地质层位信息，即在反演之前，在自动拾取主要地质层位的立体层析数据空间时就需要加入层位编号.实施地层格架封闭分为以下5个步骤：
(1) 获取层位编号信息，获取块体与地层、地层与地层之间的几何拓扑关系，即确定块体与层位之间、层位与层位之间的拓扑关系(比如一个交点由哪些层位组成，一个块体由哪些层位包围).
(2) 通过反射点信息拟合各个层位：对每个层位的数据点根据横坐标进行排序后，然后对同一个层位的数据进行NUBS拟合，利用结构张量算法提取拟合后界面的切向量和法向量信息(Wu and Hale, 2015).
(3) 处理界面与边界之间的关系：在确定层位之间的关系时，统计每个层位的交点个数(只有0、1、2三种情况)，当只有一个交点时，判断哪个端点需要延拓到边界，当没有交点时，地层两端根据当前的界面端点处的斜率延拓到边界，斜率信息是通过结构张量算法获得的.
(4) 处理界面之间的拓扑关系：对于每个交点，确定哪些界面相交于该点，然后通过一些处理使得这些层位相交于一点，处理过程中同样需要结构张量获取的界面的切向量和法向量信息.
(5) 根据块体与层位之间的拓扑关系求取块体顶底面的深度坐标，每一个块体内部一个横坐标的顶底面深度是实现分块求取Fréchet导数的重要参数.
图 4显示了基于离散的反射点实施地层格架封闭的全过程.图 4a显示了块体与界面之间、界面与界面之间的拓扑关系；图 4b显示了如何通过反射点信息拟合各个反射界面；图 4c显示了将界面延伸到边界后的情况；图 4d显示了最为重要也是相对复杂的一步：处理好界面与界面之间的交点关系.
|
图 4 块状地层封闭的过程 (a)统计块体与界面之间、界面与界面之间的拓扑关系；(b)通过反射点信息拟合各个反射界面；(c)处理界面与边界关系的情况；(d)处理界面与界面之间关系(交点)情况 Fig. 4 Process of the closure of block-like formation (a) Topological relationship between the block and the interface and the interface and the interface; (b) Curve fitting by information of reflection points; (c) Processing of relationship between interface and boundary; (d) Processing of relationship between interface and interface |
到目前为止我们已经实现了：(1)推导了二维地层格架约束立体层析的Fréchet导数；(2)利用NUBS插值实现基于离散反射点的连续界面拟合；(3)利用初始成像界面提供的地质块与界面的拓扑结构实现地层格架的自动封闭.根据上述这些算法可以得到一个实用的、适用于二维实际数据的工作流程.该流程实现过程如下：(1)对输入的叠前地震数据利用结构张量算法提取高密度立体层析数据空间[Sx, Rx, T, Psx, Prx]；(2)然后实施传统立体层析得到初始速度模型，这个过程是全自动的，无需人工干预(王宇翔等，2016)；(3)基于初始立体层析提供的反射点位置拟合关键地质层位并对其实施封闭得到初始地层格架信息；(4)利用初始地层格架信息对每一个地质块内的速度利用简单的算术平均得到每一个块内初始速度；(5)运行本文提出的地层格架约束立体层析得到最终的地质一致性块状模型.以上实现思路可总结为如图 5所示的流程图.
|
图 5 保持边界二维立体层析反演实现流程 Fig. 5 Flow chart of 2D stereo-tomography inversion with kept boundary |
图 6显示这个工作流程在图 2所示的六层背斜模型理论数据上是如何工作.基于图 2显示的模型，共正演了6000个射线对并得到对应的立体层析数据空间.首先运行一个传统立体层析得到初始反射点位置和初始速度模型(如图 6a所示).图 6b显示基于图 6a中获得的初始反射点位置实施NUBS拟合可以获得6个连续的反射界面.图 6c显示根据图 6b所示反射界面实施地层格架封闭后的结果.地层格架约束立体层析的初始模型是将图 6c所示的地层格架内充填一个初始速度后得到的.初始充填速度的获得方式和图 5的建议流程完全一致，即采取在图 6a所示的初始模型内每一个地质块内做一个简单的算术平均后得到的.图 6d显示了地层格架约束二维立体层析在这个模型上的最终反演结果.我们可以看出图 6d所示的反演结果和图 2所示的真实模型已经非常接近，证明这个流程在图 2所示的层状模型上是有效的.
|
图 6 二维地层格架约束立体层析应用于六层背斜模型 (a)常规初始立体层析反演结果；(b)基于常规初始立体层析反演提取初始地层格架信息；(c)对初始格架信息进行封闭后的结果；(d)最终地层格架约束立体层析反演结果. Fig. 6 Application of 2D stereo-tomography constrained by stratum framework to a 6-layer anticline model (a) Inversion result of conventional stereo-tomography; (b) Stratum framework information from conventional stereo-tomography; (c) The result after the closure of initial framework data; (d) Final inversion result. |
图 7显示这个工作流程在另一个更为复杂的理论数据上如何工作.基于图 7a显示的模型，正演了一条2D理论测线，共激发201炮，炮间距为40 m，每炮401道接收，道间距为10 m，最大偏移距2 km.在自动搜索到5700个数据点后，首先运行一个传统立体层析得到初始反射点位置和初始速度模型(如图 7b).基于图 7b显示的反射点位置通过NUBS拟合算法可以实现地层格架自动封闭(如图 7d所示).图 7e显示了基于图 7b显示的初始模型和图 7d所示的封闭后地层格架在每一块内做一个简单的平均后得到了用于二维地层格架约束立体层析的初始模型.在图 7e基础上，我们运行了地层格架约束立体层析.图 7f、g分别显示地层格架约束立体层析第1次和第10次迭代更新的速度模型.我们可以看出图 7g所示的最终反演结果和图 7a所示的真实模型已经非常接近.图 7h为目标函数的下降曲线，可以看出最终的二范数残差已经非常小.
|
图 7 带边界的立体层析应用于合成数据 (a)真实模型；(b)传统的立体层析得到的初始模型；(c)从传统的立体层析中反演得到的反射点位置；(d)用NURBS得到的封闭的地质格架；(e)在每一块中通过简单的平均重新分配速度；(f)第一次迭代；(g)第十次迭代；(h)成本函数的收敛曲线. Fig. 7 Application of stereo-tomography with boundaries to synthetic data (a) True model; (b) Initial model inverted from conventional stereo tomography; (c) Reflector positions inverted from conventional stereo tomography; (d) A closed geological framework by NUBS; (e) Re-assigned velocity via a simple average within each block; (f) Updated velocity in 1st iteration; (g) Updated velocity in 10th iteration; (h) Convergence curve of the cost function. |
图 8a、b显示了用图 7b所示传统立体层析模型实施偏移成像得到的共成像点道集(CIG)以及用图 7g所示地层格架约束立体层析反演模型实施偏移成像获得的CIG之间的对比.图 8c、d显示了用图 7b所示传统立体层析模型实施偏移成像得到的成像剖面以及用图 7g所示地层格架约束立体层析反演模型实施偏移成像获得的成像剖面的对比.可以看出，本文方法获得的深度成像剖面的聚焦程度和成像品质相比传统立体层析有明显提高.
|
图 8 应用传统立体层析模型和地层格架约束立体层析模型的叠前深度偏移成像对比 (a)应用传统方法获得的CIG；(b)应用本文方法获得的CIG；(c)应用传统方法得到的最后深度成像；(d)应用本文方法得到的最终深度成像. Fig. 8 Comparison of PSDM images between conventional stereo tomography and stratum-frame-constrained stereo tomography (a) CIGs with conventional method; (b) CIGs with presented method; (c) Final stacked image with conventional method; (d) Final stacked image with presented method. |
二维实际数据是2011年在南海深水区采集获得，共选取了该测线中950炮，炮间距为50 m, 道间距为25 m，最大偏移距为8275 m.考虑到实际模型的复杂性，每一层内的速度设置为随纵、横向线性变化.图 9a显示了传统立体层析得到的最终反演模型，其中反演得到的倾角条用蓝色线画在速度模型之上.
|
图 9 地层格架约束立体层析应用于南海的实际数据 (a)传统立体层析反演模型；(b)带地层格架的初始模型；(c) 7次迭代后的模型；(d)反演过程中的成本函数下降曲线；(e)由图(a)的模型实施叠前深度偏移得到的CIG; (f)由图(c)的模型实施叠前深度偏移得到的CIG；(g)由图(a)模型得到的最终深度偏移图像；(h)由图(c)模型得到的最终深度偏移图像. Fig. 9 Application of stratum-frame-constrained stereo tomography to real data from the South China Sea (a) Conventional stereo tomographic inversion model; (b) Initial model with stratum frame; (c) The model after 7th iteration; (d) The cost function curve during the inversion; (e) Some CIGs generated with model in (a); (f) Some CIGs generated with model in (c); (g) PSDM stacked image with model in (a); (h) PSDM stacked image with model in (c). |
应用同样的工作流程，在NUBS拟合算法和结构张量算法的帮助下，我们得到了一个封闭的地质格架.由于该数据信噪比的原因，仅仅基于了几个关键层位建立了一个比较粗放的地层格架，没有针对更精细的反射点信息实施地层格架封闭.尽管如此，地层格架约束立体层析对于深部的速度结构反演依然有明显改善.图 9b、c显示地层格架约束立体层析的第1次和第7次迭代结果.图 8d显示迭代过程中目标函数的下降曲线.由于实际数据存在信噪比以及数据点提取的密度和精度问题，因此其目标函数的下降幅度不可能像理论数据的目标函数那样，可以下降到接近于零的水平.但是在引入地层格架之后，对于测线中部5000~8000 m深度的内幕的成像依然有明显改进.图 9e、f显示了用图 9a所示传统立体层析模型实施偏移成像得到的共成像点道集(CIG)以及用图 9c所示地层格架约束立体层析反演模型实施偏移成像获得的CIG之间的对比.图 9g、h显示了用图 9a所示传统立体层析模型实施偏移成像得到的成像剖面以及用图 9c所示地层格架约束立体层析反演模型实施偏移成像获得的成像剖面的对比.从CIG对比可以看到本文方法反演模型得到的深度成像有更好的聚焦；从成像剖面对比可以看出成像品质确有所提高.再一次验证了算法的正确性和工作流程的合理性.
6 结论本文提出一种在地层格架正则化约束下的二维立体层析反演方法，该方法的关键之处在于把地层格架信息作用于到立体层析Fréchet导数矩阵中，使得更新后的速度模型自动呈现出符合地质规律的块状特征.地质层位信息基于反演迭代中得到的反射点位置进行B样条插值拟合得到，因此地层格架信息会随着反射点位置的更新而更新.与前人算法相比，本文所提出的正则化方法更为直接也更容易实现.理论和实际数据算例证明该正则化方法的稳健性和可靠性，能够得到与地质认识一致的反演结果.三维地层格架约束立体层析的实用性，从理论上没有任何问题.实践中只要解决了三维地质层位的自动封闭，本文提出的工作流程就可以顺利推广到三维.该方法能够自然的扩展到三维，进一步提升三维立体层析速度建模的品质.
Alerini M, Lambaré G, Baina R, et al. 2007. Two-dimensional PP/PS-stereotomography:P-and S-waves velocities estimation from OBC data. Geophysical Journal International, 170(2): 725-736. DOI:10.1111/gji.2007.170.issue-2 |
Billette F, Lambaré G. 1998. Velocity macro-model estimation from seismic reflection data by stereotomography. Geophysical Journal International, 135(2): 671-690. DOI:10.1046/j.1365-246X.1998.00632.x |
Billette F, Le Bégat S, Podvin P, et al. 2003. Practical aspects and applications of 2D stereotomography. Geophysics, 68(3): 1008-1021. DOI:10.1190/1.1581072 |
Clapp R G, Biondi B L, Claerbout J F. 2004. Incorporating geologic information into reflection tomography. Geophysics, 69(2): 533-546. DOI:10.1190/1.1707073 |
Costa J C, da Silva F J C, Gomes E N S, et al. 2008. Regularization in slope tomography. Geophysics, 73(5): VE39-VE47. DOI:10.1190/1.2967499 |
Farra V, Madariaga R. 1987. Seismic waveform modeling in heterogeneous media by ray perturbation theory. Journal of Geophysical Research, 92(B3): 2697-2712. DOI:10.1029/JB092iB03p02697 |
Farra V, Madariaga R. 1988. Non-linear reflection tomography. Geophysical Journal International, 95(1): 135-147. DOI:10.1111/j.1365-246X.1988.tb00456.x |
Harlan W S. 1995. Regularization by model reparameterization. http://www.billharlan.com/pub/papers/regularization.pdf.
|
Kosloff D, Sherwood J, Koren Z, et al. 1996. Velocity and interface depth determination by tomography of depth migrated gathers. Geophysics, 61(5): 1511-1523. DOI:10.1190/1.1444076 |
Lambaré G, Guillaume P, Zhang X, et al. 2013. Multi-layernon-linear slope tomography.//75th EAGE Conference & Exhibition. EAGE.
|
Li Z W, Yang K, Ni Y, et al. 2014. Migration velocity analysis with stereo-tomography. Geophysical Prospecting for Petroleum (in Chinese), 53(4): 444-452. |
Ni Y, Yang K, Chen B S. 2013. Stereo-tomography inversion method:theory and application testing. Geophysical Prospecting for Petroleum (in Chinese), 52(2): 121-130, 111. |
Piegl L, Tiller W. 1997. B-spline programming concepts.//The NURBS Book. Berlin, Heidelberg: Springer, 593-628.
|
Prieux V, Lambaré G, Operto S, et al. 2013. Building starting models for full waveform inversion from wide-aperture data by stereotomography. Geophysical Prospecting, 61(S1): 109-137. |
Ren L J, Sun X D, Li Z C, 2014. The stereo-tomography based on eigen-wave attributes.//Proceedings of the CPS/SEG Beijing 2014 International Geophysical Conference (in Chinese). Beijing: CPS/SEG.
|
Shao W D, Yang K, Xing F Y, et al. 2017. Stereo-tomography in squared-slowness triangulated model. Chinese Journal of Geophysics (in Chinese), 60(9): 3574-3586. DOI:10.6038/cjg20170923 |
Van Vliet L J, Verbeek P W. 1995. Estimators for orientation and anisotropy in digitized images.//ASCI Imaging Workshop. Venray, NL, 442-450.
|
Wang Y X, Yang K, Yang X C, et al. 2016. A high-density stereo-tomography method based on the gradient square structure tensors algorithm. Chinese Journal of Geophysics (in Chinese), 59(1): 263-276. DOI:10.6038/cjg20160122 |
Wu X M, Hale D. 2015. Horizon volumes with interpreted constraints. Geophysics, 80(2): IM21-IM33. DOI:10.1190/geo2014-0212.1 |
Yang K, Xing F Y, Li Z W, et al. 2016. 3D stereo-tomography based on the non-reduced Hamiltonian. Chinese Journal of Geophysics (in Chinese), 59(9): 3366-3378. DOI:10.6038/cjg20160920 |
Zelt C A, Smith R B. 1992. Seismic traveltime inversion for 2-D crustal velocity structure. Geophysical Journal International, 108(1): 16-34. DOI:10.1111/gji.1992.108.issue-1 |
Zhang X, Zhang J. 2012. Edge preserving regularization for seismic traveltime tomography.//2012 SEG Annual Meeting. Las Vegas, Nevada: SEG, 2830-3834.
|
Zhou H W. 2006. Multiscale deformable-layer tomography. Geophysics, 71(3): R11-R19. DOI:10.1190/1.2194519 |
李振伟, 杨锴, 倪瑶, 等. 2014. 基于立体层析反演的偏移速度建模应用研究. 石油物探, 53(4): 444-452. DOI:10.3969/j.issn.1000-1441.2014.04.010 |
倪瑶, 杨锴, 陈宝书. 2013. 立体层析反演方法理论分析与应用测试. 石油物探, 52(2): 121-130. |
任丽娟, 孙小东, 李振春. 2014.基于特征波属性参数的立体层析速度反演方法研究.//CPS/SEG北京2014国际地球物理会议摘要.北京.
|
邵炜栋, 杨锴, 邢逢源, 等. 2017. 慢度平方三角网格立体层析反演方法. 地球物理学报, 60(9): 3574-3586. DOI:10.6038/cjg20170923 |
王宇翔, 杨锴, 杨小椿, 等. 2016. 基于梯度平方结构张量算法的高密度二维立体层析反演. 地球物理学报, 59(1): 263-276. DOI:10.6038/cjg20160122 |
杨锴, 邢逢源, 李振伟, 等. 2016. 基于非降阶汉密尔顿算子的三维立体层析反演. 地球物理学报, 59(9): 3366-3378. DOI:10.6038/cjg20160920 |