2016, Vol. 59
2. Broadgeo, Inc, Houston, TX 770994;
3. 中海石油研究中心, 北京 100027
2. Broadgeo, Inc, Houston, TX 770994;
3. CNOOC Research Center, Beijing 100027, China
立体层析反演是层析反演类速度建模方法中极具特色的一种方法.其将一个局部相干同相轴在炮道集与检波点道集内的射线参数水平分量(Psx，Prx)纳入到反演的数据空间之中，不仅使得数据空间的准备相比传统反射层析更为简便，也使得立体层析反演成为运动学层析反演方法中唯一一种可以同时反演速度结构与反射点位置的方法(Billette and Lambaré，1998；Lambare，2008；倪瑶等，2013；任丽娟等，2014；李振伟等，2014).近年来亦有尝试将立体层析反演与有限频理论相结合，以期获得更好的反演结果(Ni et al.，2012).
在立体层析反演中，数据空间可以写为d=[Sx，Rx，T，Psx，Prx].其中Sx，Rx为炮检点坐标；T为总走时；Psx，Prx为炮检点处的射线参数水平分量.注意对有效数据点的确定完全是根据局部波包的Psx，Prx是否足够精确、可靠作为标准.换言之，数据点的所有信息中，Psx，Prx首先得到确定，之后再确定该局部波包对应的炮检点坐标Sx，Rx和走时T.因此，射线参数水平分量的精度问题是立体层析反演数据准备过程中首要关心的问题.
由于射线参数的确定实质就是同相轴的局部倾角估计问题，前人大多使用局部倾斜叠加或平面波摧毁滤波器来获得这类信息.本文引入图像处理界中发展成熟的边缘检测算法——结构张量算法来估计同相轴的局部倾角.石油工业界中已有人将该类算法用于地震解释中的自动层位追踪或者提高信噪比的构造预测滤波(Wu and Hale，2013; Hale，2009).不同于上述两个文献，本文使用的是基于梯度平方的结构张量算法.实际测试表明该算法具有很好的抗噪性，将其用于实际地震数据的局部倾角估计是非常合适的.引入该算法的另一个重要理由是它的计算效率远高于刚才提到的两种常规算法.该算法的引入将使得立体层析数据点的分布密度大幅提高而无需担心计算成本的大幅上升.这是实施高密度立体层析反演的一个基本前提保证.
本文首先简要回顾立体层析反演的方法原理并阐明射线参数水平分量对立体层析数据空间中的重要地位.然后详细介绍梯度平方结构张量算法，并将其与局部倾斜叠加和平面波摧毁滤波器的计算效率和性能进行对比，讨论如何将该算法应用于叠前地震数据、从而实现智能化的自动倾角拾取.在此基础上，作者提出了一种基于物理规律的严格质量监控(QC)方法来过滤非物理的、冗余的数据点.需要指出，正是由于梯度平方结构张量算法的高效以及合乎物理规律的严格QC的引入，使得适用于实际数据的高密度立体层析反演成为可能.基于二维理论数据与南海某二维深水地震数据的测试表明，上述策略的使用不但使得数据点的数量大幅度增加，同时也保证了数据点的有效性和可靠性.这对于降低层析矩阵的病态性、提高反演的精度是至关重要的.
2 立体层析反演方法原理介绍首先介绍立体层析反演的模型空间、数据空间以及层析矩阵的建立过程.为简明起见，以二维为例.图 1a显示了一根从炮点S出发、到地下反射点X反射、回到地表R的射线.从透射的角度，不妨将 其理解为从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代表地下介质速度信息.
 | 图 1 二维立体层析数据空间各分量与模型空间各分量 (a) 立体层析的模型与数据分量; (b) 同相轴斜率与射线理论的慢度矢量.Fig. 1 The model components and the data components in 2D stereotomography data space (a) The model components and the data components in 2D stereotomography data space; (b) The relation between the slope of the seismic events and the slowness vector in ray theory. |
注意图 1b中所示的Psx可以在共检波点道集内搜索同相轴的局部倾角得到，Prx可以在共炮点道集内搜索同相轴的局部倾角得到. 原理如下：在图 1b所显示的炮记录中，如果考察其中某一局部相干同相 轴，由几何关系以及射线理论中慢度矢量的定义得
Billette和Lambare(1998)指出，正是由于Prx，Psx(共炮道集和共检波点道集数据内同相轴的斜率)被引入层析数据空间，使得在层析反演中反 射/散射波在炮检点处的局部传播方向得到强有力 的约束.具体的优势体现在立体层析反演中无需拾取连续的反射同相轴，只要选择连续性好、信噪比高的局部波包即可.同时，将反射问题化为透射问题，也回避了传统反射层析中速度与反射层深度耦合的问题.由于在数据空间中引入了Prx，Psx，模型空间必然也随之扩大.因此在反演速度的同时必须还要反演反射点的深度位置以及反射层的倾角.这使得立体层析反演成为运动学层析反演方法中唯一一种可以同时反演速度、反射点位置以及局部倾角的层析反演方法.而常规反射层析反演需要在速度更新之后再单独更新反射层深度，这通常是一个需要分两步实施的过程.
公式(2)展示了二维立体层析反演的层析矩阵(或FRECHET偏导数矩阵)，其中σ为针对层析矩阵每一行物理量的数量级不同设置的均衡系数.可以看到常规层析反演由于其模型空间一般仅由速度V构成，其数据空间一般仅考虑走时T，可认为其层析矩阵仅相当于公式(2)中的右上角那一部分.可以看出立体层析矩阵无论从规模还是其稀疏性方面都超过了常规的走时层析矩阵.如果能够实现高密度、高质量的数据空间提取，对于改善层析矩阵的条件数，降低求解时对规则化的要求，进而提高求解精度显然都是有意义的.
Lambare(2008)指出，选择立体层析数据空间的标准就是优先考虑连续性好、信噪比高的一次反射局部波包.从数据空间的特点不难看出，一旦选择好了Psx和Prx参数，其他参数如炮检点横坐标和走时都可以轻易随之得到.如前所述，Psx和Prx参数估计等价于同相轴局部倾角估计.换言之，如何准确、 快速地估计共炮点道集与共检波点道集内同相轴的 局部倾角是立体层析反演中最为重要的数据准备工作.
3 梯度平方结构张量算法原理及其特点工业界两种最常用的局部倾角提取手段当属局部倾斜叠加(Ottolini，1983)与平面波摧毁滤波器(Fomel，2002).前者在时空域或者频率空间域针对目标同相轴实施局部叠加，通过叠加能量最大找到最合适的局部倾角；后者主要通过对局部平面波方程进行有限差分分解，在频率域内构建预测算子，然后再转换到时间域时对频率域相移算子时进行近似得到.这两种算法都被工业界大量使用.但是如果从计算成本的角度来考量，图像处理领域的梯度平方结构张量算法则具有巨大的优势.
图 2显示了南海某二维深水测线中的一个单炮记录，其信噪比不是很理想.从图像处理的角度，这个信噪比不高的地震记录可以视作一张反差不太强烈的二维图像.而针对一张信噪比不高的地震记录如何去估算其每一个样点处的局部倾角和图像处理中对于反差不够强烈的图像如何实施边缘检测，本质上是同一个问题.梯度平方结构张量算法正是针 对这个问题而提出的(Van Vliet and Verbeek，1995).
 | 图 2 某海洋二维测线原始炮记录Fig. 2 A raw shot record in a marine 2D seismic line |
对于任何一张图像，为检测其边缘，必须首先计算出每一个样点处X方向与Y方向的梯度gx，gy.而高效计算每一样点的梯度在图像处理界早已有成熟算法，这就是结构张量算法中使用的迭代高斯滤波.需要指出，这正是结构张量算法具有高计算效率的原因所在.关于迭代高斯滤波的技术细节可以参考相关文献(Van Vliet et al.，1998).
针对低信噪比图像，为了提高方向预测的稳定性，避免梯度矢量在图像边缘两侧互相抵消，Van Vliet和Verbeek(1995)引入了所谓的梯度平方张量矩阵，对图像中某一样点而言，其梯度平方张量矩阵的定义如下：
获得光滑的梯度平方张量矩阵 G′ 之后，针对半正定矩阵 G′，求解其特征值与特征向量，具体可由求解特征方程｜ G′ - λI ｜=0得到：
特征向量则描述了图像局部线性结构的方向性.如图 2上的白色箭头就是作者使用该算法计算中的两个特征方向，特征向量 V 1正交于图像的主结构方向，特征向量 V 2平行于图像的主结构方向.图 3则是Van Vliet和Verbeek(1995)对一张信噪比欠佳的岩芯照片实施该算法后的结果.可以看到各类参数表达的图像信息具有良好的一致性，表现出该算法具有不错的稳定性和抗噪性.
 | 图 3 梯度平方张量算法应用效果(Van Vliet和Verbeek(1995)) (a) 输入图像; (b) 特征值λ1剖面; (c) 特征值λ2剖面; (d) 局部走向剖面; (e) 线性度(同相性)剖面.Fig. 3 The illustration of the gradient square tensors (Van Vliet和Verbeek(1995)) (a) Input image; (b) Eigen value λ1; (c) Eigen value λ2; (d) Local orientation; (e) Local linearity. |
图 4展示了对于图 2所示的实际单炮记录，应用该算法提取的四种属性参数剖面.可以看出由于梯度平方结构张量算法的稳定和高效，它能够胜任针对实际地震数据的局部倾角估计工作.
 | 图 4 梯度平方张量算法在某单炮记录上的应用效果 (a) 输入图像; (b) 特征值λ1剖面; (c) 特征值λ2剖面; (d) 局部走向剖面; (e) 线性度(同相性)剖面.Fig. 4 The illustration of the gradient square tensors (a) Input image; (b) Eigen value λ1; (c) Eigen value λ2; (d) Local orientation; (e) Local linearity. |
对于立体层析而言，似乎只要知道特征方向 V 2就足够了.实际上，我们同样需要知道特征方向 V 1，以及相应的最大与最小特征值λ1，λ2.因为线性度是非常重要的一个参数，其计算依赖于最大与最小特征值λ1，λ2.在后面的章节将会提及，线性度的大小正是评价局部波包同相性的一个重要指标，也是我们评价数据点质量的第一个标准.
图 5a显示了针对一个纵横向样点数为200×200的测试数据，其中含有少量随机噪声.应用局部倾斜叠加、平面波摧毁滤波器与梯度平方结构张量方法估算出的图像局部倾角信息分别如图 5b，5c，5d所示.三种算法得到的局部倾角信息大致相似，但是计算效率则有明显不同.表 1中是三种算法的计算时间对比.
 | 图 5 三种倾角估计算法结果比较 (a) 理论沉积数据模型; (b) 使用局部倾斜叠加获得的局部倾角剖面; (c) 使用平面波摧毁获得的局部倾角剖面; (d) 使用梯度平方结构张量算法获得局部倾角剖面.Fig. 5 The comparison among three dip-angle estimation algorithms (a) Synthetic sedimentary model; (b) Local slopes estimated with the local slant stack; (c) Local slopes estimated with the plane-wave destructor; (d) Local slopes estimated with the gradient based structure tensor algorithm. |
 | 表 1 三种算法的计算时间对比(单位：s) Table 1 Comparison among the computational time of three algorithms (unit: s) |
从表 1可见梯度平方结构张量算法的计算效率高于平面波摧毁近一个数量级，高于局部倾斜叠加达两个数量级.注意梯度平方结构张量算法获得的局部倾角信息是一个邻域内的平均效应.因此其变化比较光滑连续，不像局部倾斜叠加获得的局部倾角信息似有突变.但是这并不表示其分辨率低于局部倾斜叠加，相反其倾角估计的稳定性恰恰来自这种邻域平均效应.需要指出，即便考虑了这个邻域平均效应，它在计算成本方面相对于另外两种常规算法依然存在巨大优势.这使得它成为立体层析反演实际应用中，提取射线参数水平分量的首选算法.
4 基于严格质量监控的自动拾取流程梯度平方结构张量算法的高效、稳定使得它成为自动化拾取流程中的首选.其设计思路如下：首先，对于共偏移距剖面中线性度比较大的样点位置，利用某些准则判断其是否位于波峰位置.如果确实位于波峰位置，则读取其局部倾角，然后换算为射线参数水平分量信息Pmx；其次，根据走时信息与炮检点坐标位置，计 算该局部同相轴在炮道集与检波点道集中的射线参数水平分量 Prx，Psx. 然后基于三者间应满足的定量关系 进行质量监控(QC)(图 6).该定量关系见4.2节详细介绍.
 | 图 6 基于梯度平方结构张量算法的自动拾取流程Fig. 6 The gradient-square-structure-tensors algorithm based workflow |
如前所述，对共偏移距剖面进行梯度平方结构张量计算后得到了每个点的线性度与结构向量.注意线性度实质上可以理解为局部波包的同相性指标，可以通过该指标判断数据点与相邻样点的相关度.从方便拾取的角度出发，我们同时希望立体层析所需的数据点都采集在波峰上.为了判断波峰位置，进一步通过已经获得的方向梯度(gx，gy)与梯度平方结构张量向量 V 1来计算沿同相轴垂直方向的梯度d：
 | 图 7 梯度平方结构张量算法获得的波峰位置 (a) 原始共偏移距数据(理论数据); (b)捕捉到的波峰位置 (线性度大于0.7得到保留).Fig. 7 The wave-crest position picked by the gradient-square-structure-tensors algorithm (a) The raw synthetic common offset section； (b) The wave-crest position picked by some coherency criterion (Linearity is bigger than 0.7). |
4.1节中数据点的选择策略可以认为是基于图 像信息本身的特点实施了一次初选.但是叠前地震数据具有自己固有的物理波动规律.注意立体层析所需要的是反射或绕射信息，仅仅根据数据本身的线性度(或同相性)作为入选标准，在实际应用中几乎可以肯定会有不少非反射、非绕射的信息入选.这些数据点对立体层析而言缺乏物理意义，属于应该剔除的无效信息.
我们制定的物理QC准则是：对通过初选的每个局部反射同相轴，首先统计其Pmx信息，再根据走时信息、炮检点坐标信息计算其在相关的炮道集、检波点道集中的射线参数信息Prx 与Psx.然后与共偏移距剖面中的斜率Pmx进行比较，如果满足Psx+Prx=Pmx(推导见附录)的精度准则，则为有效拾取，保存拾取结果，不满足则剔除.这个物理QC准则对实际数据的反演尤为重要.在实际应用中，如果没有上述物理准则的严格QC，必然会有相当数量的不满足物理规律的、冗余的数据点入选.这些数据点对反演算法而言完全属于噪声，将会对反演结果造成负面影响.
5 二维理论数据算例二维理论数据基于图 8a所示的一个六层盐丘模型得到.我们采用了三角剖分方式建模，基于反射射线正演的方法获得炮数据.共计模拟401炮，每炮401道，道间距10 m，炮间距20 m.图 8a显示了原始速度模型以及某一炮对应的射线路径.图 8b显示了反演所用的初始模型.这是一个常规的垂直梯度模型.
 | 图 8 基于射线正演的理论数据及用于反演的初始模型 (a) 真实速度模型及某一炮的射线路径; (b) 初始速度模型.Fig. 8 The synthetic data based on ray modeling and the initial model for stereotomography (a) The true velocity model and the ray paths related to one shot; (b) The true velocity model and the ray paths related to one shot. |
图 9显示了正演得到的共偏移距道集、共炮道集、共检波点道集以及基于这些记录利用梯度平方结构张量算法获得的数据点信息(与局部倾角信息联合显示).搜索时，首先利用结构张量算法提供的各种信息，结合波峰判断准则和线性度准则(这里依然选择线性度大于0.7的同相轴信息)获得地震记录上的有效数据点，图 9(a，b，c)展示了基于上述原则获得的数据点和局部倾角信息.可以看到这种方式获得的数据点信息是相当准确的.
 | 图 9 理论数据梯度平方结构张量数据空间显示 (a) 0 m共偏移距数据点与倾角信息联合显示; (b) 共炮道集数据点与倾角信息联合显示; (c) 共检波点道集数据点与倾角信息联合显示; (d) 0 m共偏移距数据点与倾角信息联合显示(红点代表目标数据点); (e) 目标数据点在对应的共炮道集(第1炮的 201道)上与倾角信息(蓝色线条)联合显示; (f) 目标数据点在对应的检波点道集(第1个共检波点道集的第1道)上与倾角信息 (蓝色线条)联合显示. Fig. 9 The data space obtained by the presented workflow (a) A joint display of all data points in zero-offset section and their dip-bar (in blue color)； (b) A joint display of all data points in shot No.1 and their dip-bar (in blue color)； (c) A joint display of all data points in receiver No.1 and their dip-bar (in blue color)； (d) A joint display of one data point (in read color) in zero-offset section and its dip-bar； (e) A joint display of one data point at channel 201 of shot No.1 and its dip-bar (in blue color); (f) A joint display of one data point at channel 1 of receiver No.1 and its dip-bar (in blue color). |
图 9(d，e，f)显示了在获得数据点信息之后，如何实施严格物理QC的过程.注意在0 m共偏移距剖面上用红色标注了一个目标数据点(如图 9d上红点所示)，其搜索到的Pmx=-0.054300 ms·m-1.根据道头信息，可以找到在对应共炮道集的相应位置(第1炮的第201道，蓝色线条代表倾角条)和对 应的共检波点道集(第1个共检波点道集的第1道，蓝色线条代表倾角条)的位置，利用结构张量 算法可以分别计算得到 Psx=-0.031207 ms·m-1和Prx=-0.023089 ms·m-1.计算Psx+Prx-Pmx=0.000004 ms·m-1，发现差值小于事先指定的精度门槛0.00001 ms·m-1，符合精度要求，因此这个数据点信息得到了保留.
为了测试梯度平方结构张量算法与物理QC的作用，在这个正演记录中还加入了一定量的随机噪声.可以看到无论是共偏移距记录还是共炮点记录或共检波点记录，通过线性度指标和物理QC准则的联合作用，梯度平方结构张量算法展示了高精度搜索同相轴局部倾角的能力.更重要的是，搜索这些倾角信息的计算成本相比局部倾斜叠加下降近两个数量级.本测线共401炮，采用200 m偏移距的间隔，在48 个核的并行计算，仅用了40 min就获得了所有的数据点信息.如果采用局部倾斜叠加，需接近20 h才能完成同等规模的计算.即便应用同样的线性度标准，局部倾斜叠加获得的数据点密度也明显少于结构张量算法.
对这条理论二维数据，以200 m偏移距为间隔，在不同的偏移距剖面上初选了共3万个数据点用于反演.图 10a显示了没有实施物理QC的反演结果.由于数据点密度很大，也获得了大致合理的速度结构.但是由于线性噪声的加入，使得部分数据点缺乏物理意义，导致靠近盐丘顶部的第三、四、五、六层反射层倾角反演结果出现异常，盐丘结构发生了变异.这时加入了物理QC的过滤准则.给出了非常苛刻的精度门槛(0.00001 ms·m-1).这种严格物理QC的加入使得合乎要求的数据点骤降到1万个左右，但是其反演结果是相当稳定的.如图 10b所示，虽然在盐丘两侧的倾角信息出现了一些不连续，这显然与过于严格的过滤准则有关.但剩下的数据点质量很高，反演结果中没有出现异常的、违反地质规律的现象.
 | 图 10 理论数据反演结果对比 (a) 未实施物理QC； (b) 实施物理QC后.Fig. 10 The comparison between inverted velocity before and after physical QC (a) Before the implementation of physical QC; (b) After the implementation of physical QC. |
选择南海某二维深水测线进行立体层析反演处理.整条测线长150 km，共2935炮，最大偏移距8275 m，炮间距50 m，检波点间距25 m.整条测线包括浅水、深水两个部分，海水深度变化从200 m到2000 m.反演之前的数据做了精细的前处理：如压制与海面有关的多次波、实施预测反褶积提高分辨率等等.在数据拾取之前我们进行了10～45 Hz的窄带滤波器并实施短时窗的自动增益以保证同相轴可以识别，并且切除了直达波以上的数据来减少拾取工作量，使得预处理后的数据有较高的信噪比来保证拾取的可靠性.
由于梯度平方结构张量算法的高效，使得本次反演的数据空间提取效率大为提高.反射波拾取范围划定在100~6200 m偏移距，时间范围在8 s以内.每隔200 m偏移距进行一次拾取.图 11展示了在南海某二维深水数据的300 m共偏移距剖面上，应用物理QC之后的数据点前后的变化.可以看到，在应用物理QC过滤准则后，一部分拾取点被删除了.由于采取了比较严格的过滤条件，图 11c中展示的数据点信息都是非常可靠的.
 | 图 11 共偏移距剖面上基于结构张量算法实现数据空间提取 (a) 原始数据的300 m偏移距的共偏移距道集； (b) 仅利用数据信息自动拾取出的位置，绿色小点处为波峰； (c) 物理QC准则应用于图b之后留下的数据点，绿色小点为拾取点，蓝色斜线为该点处局部倾角.Fig. 11 The extraction of data space with structure tensors algorithm in common offset section (a) Raw 300 m CO section; (b) Data points picked by the linearity threshold (0.7) and their dip-bar information (in blue color); (c) Data points picked after physical QC and their dip-bar information (in blue color). |
立体层析的目标泛函是基于拟合数据域信息建立的，因此拾取数据的精度十分重要，这是反演能否成功的最基本保证.结合梯度平方结构张量算法，实现了针对二维实际数据的高密度立体层析反演.在MPI多进程并行技术支持下，采用120进程并行计算2.5 h就全自动完成数据空间准备工作.得到了关于整条测线的14万个数据点的拾取数据量.如果使用局部倾斜叠加的方法得到同等的数据空间准备需要6~8天，而且还达不到目前结构张量算法所能提供的数据密度.梯度平方结构张量算法的应用意味着在更短的生产周期内得到了更密集均匀的数据拾取覆盖率.这对于立体层析反演的实际应用能力无疑具有积极意义.
6.2 反演策略与结果分析反演模型的规模为横向150 km，深度为10 km，为了压缩模型空间，采用了B样条表达.采用的B样条系数的横向间距和纵向间距均为200 m.即便采用了B样条表达，由于模型空间很大，射线覆盖往往不均匀，层析矩阵依然是稀疏、欠定、病态的.为了改善层析矩阵的稀疏性，降低规则化对反演结果的平滑效应，高密度的数据空间准备是非常重要的.
初始模型使用最普通的梯度模型.在初始模型上通过叠前深度偏移可以得到初始深度成像剖面，由于初始模型中对应于海水的速度是准确的，因此在初始深度成像剖面上可以拾取出海底高程信息.因此在反演过程中无需计算数据空间对海水速度的Frechet导数.同时针对海水部分给出很强的阻尼规则化参数，使得海水层内速度保持稳定.此外，在规则化参数设置方面，在没有射线或低射线覆盖区域设置了较高的阻尼系数，在高射线覆盖区域则设置较低的阻尼系数，以期在反演的稳定与精度之间获得一种平衡.
为求解这个大规模稀疏立体层析矩阵方程，我们采用了多波前多线程的稀疏矩阵QR分解算法(Davis，2011)来求解方程.在MPI环境下进行多进程并行计算射线正演以及Frechet导数求取等.在48个进程并行计算的情况下，在4 h内即完成24次迭代，得到了满意的反演结果.
图 12a显示了立体层析在24次迭代更新之后的反演结果.注意其中的蓝色线条为模型空间中的局部反射倾角信息，将其贴合到反演速度上联合显示的是立体层析中展示反演成果最常用的一种方式.由于实际数据的信噪比在不同的构造区域内有比较大的差别，因此拾取数据在信噪比高的区域比较密集，在信噪比低的区域则相对稀疏.图 12b将图 12a中20～60 km段的剖面放大显示.可以看到高密度的数据点拾取所对应的反射点位置在反演过程中自动聚集形成了有效的层位信息，和实际的成像层位位置符合的非常好.需要强调，类似的高密度反演结果是利用局部倾斜叠加提取数据空间时无法做到的.图 12c展示了基于该模型的全测线叠前深度成像结果，可以看到整体构造形态得到了精细的刻画，有利于后续的地质解释.
 | 图 12 通过结构张量算法获得的数据空间反演得到的速度模型 (a) 完整测线速度模型与模型空间倾角; (b) 局部展示20000～60000 m范围的速度反演结果与模型空间倾角; (c) 基于完整测线速度模型偏移得到的叠前深度偏移结果.Fig. 12 The stereotomography inverted model based the data space extracted by the structure tensors algorithm (a) The stereotomography inverted model with structure-tensor； (b) Zoom in from Fig.12a; (c) The PSDM image with model shown in Fig.12a. |
由于局部倾斜叠加准备全测线的数据空间过于耗时，因此这里无法提供全测线的、两种算法的反演结果对比.图 13展示了在测线末尾抽取了900炮记录后在其中分选得到的800 m偏移距剖面.其中图 13a和图 13b分别显示了局部倾斜叠加与梯度平方结构张量算法提取得到的数据空间分步对比.其中黄色线条代表射线覆盖密度，可以看出后者的数据点密度与均匀度均明显超过前者，形成了更为均匀的射线覆盖强度，这对于提高反演的稳定性显然是有利的.图 13c，图 13d展示了基于这两种数据空间，最终反演得到的速度模型.显然结构张量算法提供了更多的有效数据点，使得反演得到的结构更为精细.这种差异也体现在图 13e，13f所示的两张叠前深度成像点道集以及图 13g，13h所示叠前深度偏移剖面上.可以看出基于结构张量算法数据空间得到的反演结果在叠前深度成像之后，无论对于能量的聚焦(CDP 3800-CDP 4000)、或对于基底(CDP 5200-CDP 5500)的归位、或对于构造形态的正确刻画(CDP 4100-CDP 4400)，其精度都超过了倾斜叠加数据空间反演结果提供的深度成像剖面.
 | 图 13 测线局部(110~150 km)的倾斜叠加数据空间与梯度平方结构张量数据空间的反演与成像结果对比 (a) 局部倾斜叠加在800 m单偏移距剖面内获得的数据空间; (b) 梯度平方结构张量在800 m单偏移距剖面获得的数据空间; (c) 局部倾斜叠加数据空间反演的速度模型(110～150 km); (d) 梯度平方结构张量数据空间反演的速度模型(110～150 km); (e) 用图c所示模型实施PSDM后，在模型130 km处的若干CIG显示; (f) 用图d所示模型实施PSDM后，在模型130 km处的若干CIG显示; (g) 局部 倾斜叠加数据空间的叠前深度偏移结果(110～150 km); (h) 梯度平方结构张量数据空间的叠前深度偏移结果(110～150 km).Fig. 13 The comparison between the PSDM image with the model shown in Fig.13c and the PSDM image with the model shown in Fig.13d (a) The model space information generated in 800 m CO section with local slant stack; (b) The model space information generated in 800 m CO section with gradient square tensors algorithm; (c) The inverted model based on the data space extracted with local slant stack (110～150 km); (d) The inverted model based on the data space extracted with structure tensors (110～150 km); (e) The PSDM cigs with inverted model shown in Fig.13c; (f) The PSDM cigs with inverted model shown in Fig.13d; (g) The PSDM image with local slant stack inverted model (110～150 km); (h) The PSDM image with the gradient square tensors inverted model (110～150 km). |
本文将图像处理中用于边缘检测的梯度平方结构张量算法引入立体层析数据空间的准备.实践表明，该算法能够高效率、高密度的获得局部同相轴的倾角信息，这是实施高密度立体层析反演的基本保证.同时，作者提出一种基于物理规律的质量监控方法来过滤非物理的冗余数据点.由于梯度平方结构张量算法以及合乎物理规律的严格QC的引入，使得自动化的高密度立体层析反演成为可能.基于二维理论数据与南海某二维深水地震数据的测试表明，与常规局部倾角估计算法准备数据空间的反演结果相比，高密度的有效数据点的使用降低了层析矩阵的病态性、提高了反演的精度，改善了叠前深度成像的品质.上述特点使得立体层析反演具备了应对实际数据的能力.同时，由于梯度平方结构张量算法很容易推广到三维，因此基于三维梯度平方结构张量算法的三维立体层析反演的全自动实现也具有很好的应用潜力.
致谢 作者感谢科罗拉多矿业学院伍新民博士研究生在梯度平方结构张量算法方面的有益讨论；作者同时感谢国家自然科学基金面上项目(41274117、 41574098)、国家科技重大专项子课题(2011ZX05025-001-03)、 以及海洋地质国家重点实验室自主课题(MG20130304)对于这项研究的支持. 附录A 质量监控公式推导对原始数据集，可看做炮域数据和共偏移距域数据，关系如下：
根据半偏移距h的定义h=(s-r)/2和中点m的定义m=(s+r)/2，我们有
| [1] | Billette F, Lambaré G. 1998. Velocity macro-model estimation from seismic reflection data by stereotomography. Geophysical Journal International, 135(2):671-690. |
| [2] | Davis T A. 2011. Algorithm 915, Suite Sparse QR:Multifrontal multithreaded rank-revealing sparse QR factorization. ACM Transactions on Mathematical Software(TOMS), 38:8. |
| [3] | Fomel S. 2002. Applications of plane-wave destruction filters. Geophysics, 67(6):1946-1960. |
| [4] | Hale D. 2009. Structure-oriented smoothing and semblance. CWP Report 635, Center Wave Phenomena. |
| [5] | Ni Yao, Yang Kai. 2012. Slope tomography assisted by finite-frequency sensitivity kernel.//82th SEG Expanded Abstract. |
| [6] | Lambare. 2008. Stereotomgraphy. Geophysics, 73(5):VE25-VE34. |
| [7] | Wu X M, Hale D. 2013. Extracting horizons and sequence boundaries from 3D seismic images.//83th SEG Expanded Abstract. |
| [8] | Ottolini R. 1983. Signal/noise separation in dip space. Stanford Exploration Project, 37:143-149. |
| [9] | Van Vliet L J. Verbeek P W. 1995. Estimators for orientation and anisotropy in digitized images.//Proc. of the First Annual Conf. of the Advanced School for Computing and Imaging(ASCI'95), Heijen, NL, May 16-18, 442-450. |
| [10] | Van Vliet L J, Young I T, Verbeek P W. 1998. Recursive Gaussian derivative filters.//Proceedings of the Fourteenth International Conference on Pattern Recognition. Brisbane, Qld:IEEE, 509-514. |
| [11] | Li Z W, Yang K, Ni Y, et al. 2014. Migration velocity analysis with stereo-tomography inversion. Geophysical Prospecting for Petroleum(in Chinese), 53(4):444-452. |
| [12] | Ni Y, Yang K, Chen B S. 2013. Stereotomography inversion method theory and application testing. Geophysical Prospecting for Petroleum(in Chinese), 52(2):430-436. |
| [13] | Ren L J, Sun X D, Li Z C. 2014. The stereotomography based on eigen-wave attributes.//The expanded abstract of the CPS/SEG Beijing 2014 International Geophysical Conference(in Chinese). |
| [14] | 李振伟, 杨锴, 倪瑶等. 2014. 基于立体层析反演的偏移速度建模应用研究. 石油物探, 53(4):444-452. |
| [15] | 倪瑶, 杨锴, 陈宝书. 2013. 立体层析反演方法理论分析与应用测试. 石油物探, 52(2):121-130. |
| [16] | 任丽娟, 孙小东, 李振春. 2014. 基于特征波属性参数的立体层析速度反演方法研究.//CPS/SEG北京2014国际地球物理会议摘要. |