地层格架混合总则化及其应用

地层格架混合总则化及其应用

0基于层析的地质结构模型求解应用于偏移成像的宏观速度模型构建是反射地震图像处理的重要步骤。基于传统的成像环境分析反演速度建模，已进入行业标准偏移速度建模流程。与传统的成像环境分析方法不同，三维anova方法是数据域分析的一种方法。与以往基于时间的反射分析算法相比，该算法充分利用地震波场中所有运动信息。将本地相位轴在炮道收集和检测点收集的波束参数（以下简称p参数）的水平分量、炮点坐标和检测点坐标纳入分析反演的数据空间，并重新排列模型空间和数据空间。三维anova反演是唯一一种抗弯折射方法，可以同时反演速度、反射点位置和局部地层倾角。对于工业界用户来说,最经典的获得地质一致性块状模型的方式当属层剥离方法.该方法从20世纪90年代初即进入商业软件,至今仍为工业界广泛应用.但是层剥离方法的问题是浅层的速度和界面位置一旦确定就无法再改动,误差会从浅层慢慢积累至深层,导致深层的速度建模产生大的偏差.因此有学者开始探索是否可以在射线类层析方法中同时实现界面和速度的更新.近年来在射线类层析的研究进展方面,有两类获得地质一致性块状模型的方法值得关注.第一种方法是使用混合正则化项,比如在地质界面处使用L1范数或其他稀疏类范数使得反演结果保留模型的突变特征,而在其他地方则使用L2范数或其他光滑类范数保持反演结果的平滑(基于射线类层析方法获得地质一致性块状结构速度的目标非常明确:反演得到的速度模型应在保证地质块内足够光滑的同时尽可能维持边界的突变.对于立体层析而言,这个目标可以通过更简洁的方法来实现.本文提出一种实用的模型正则化方法,该方法的核心思想是将封闭后的地层格架信息融入到立体层析Fréchet导数矩阵中,在通过融入地层格架信息简化了立体层析Fréchet导数矩阵之后,使得更新后的速度模型自动呈现出与地质构造较为一致的的层状或者块状结果.其实现过程如下:(1)通过非均匀B样条方法对更新后的反射点位置进行自动拟合得到更新后地质层位信息;(2)利用地质块体之间的拓扑关系对地质层位信息实施自动封闭得到封闭后的地层格架信息;(3)利用封闭后的地层格架信息对Fréchet导数矩阵进行分块,每块内的速度更新量以及速度梯度更新量均为常数,由于一个模型内的地质块数目一般不会很多,这样处理的结果使得模型空间中的速度分量数目被大幅度压缩,并且在更新过程中模型自动呈现出层状或者块状地质特征.设计上述实现策略的原因是:反射点位置本身即属于立体层析模型空间的一部分,反射点位置的更新在立体层析反演的迭代过程中是可以自动得到的.相比在每一轮反演后实施偏移成像再进行人工拾取,这是一种不仅效率更高、而且更为合理的方法.同时相比于本文第1节首先介绍传统立体层析方法的数据空间和模型空间以及本文提出的地层格架约束下立体层析的数据空间和模型空间;第2节给出了地层格架正则化下的立体层析Fréchet导数矩阵的具体算法;第3节给出了基于散乱反射点的非均匀B样条拟合算法与地层格架自动封闭的实现过程;第4节给出一个实用的工作流程;第5节给出了理论与实际数据应用算例.作者期望这种特殊的正则化方式获得的立体层析反演结果将可以用于逆时偏移处理或作为全波形反演的高质量初始模型.1有限层框架下的二维三维层析反演1.1维立体层析freget偏导数矩阵公式首先介绍常规二维立体层析的模型空间、数据空间及Fréchet导数矩阵矩阵的建立过程.图1a显示了一根从炮点S出发、到地下反射点X反射、回到地表R的射线.图1b显示了射线中心坐标系中射线起点和端点扰动的几何关系.从透射的角度,不妨将其理解为从X出发、分别以入射角θ任何一种层析反演方法都必须建立数据空间扰动与模型空间扰动的线性关系,即:这样利用观测到的数据残差就可以计算出模型空间的扰动量,达到更新初始模型的目的.其中关键是建立方程(1)中所示的F矩阵,即Fréchet偏导数矩阵公式.Fréchet偏导数矩阵的物理意义是数据空间任一分量关于模型空间任一分量的敏感度.在常规走时层析中,走时t关于速度V的偏导数其实就是射线经过每一个速度单元的弧长.(2)式展示了二维常规立体层析Fréchet偏导数矩阵,除了第一行走时关于模型空间的偏导数根据走时积分方程即可得到外.其余元素是都通过射线扰动理论(Farraand1.2层位优化设计在传统立体层析中速度模型一般由B样条基函数的权系数来描述,这种表达可以达到压缩立体层析Fréchet导数矩阵规模的目的.然而,B样条描述方式的一个副作用是会使得反演结果过于平滑.与此同时,传统的正则化技术一般都建议加强解的平滑性通过一个概念实验即可清楚阐述地层格架约束立体层析的核心思想.图2显示了一个6层背斜模型.如果按照均匀B样条基函数表达这个模型,在节点纵横向间距都为400m(一般取200m或400m)时也需要480个B样条系数,也就是说模型空间中关于速度的未知数有480个.如果说用常规的规则网格表达图2所示的速度模型,在网格纵横向间隔都设为100m时,会有120×48=5760个网格,也就是说模型空间中关于速度的未知数将有5760个.然而,如果我们假定每一块内的速度都应该是常数时,其模型空间内的速度未知数其实应该是6个.基于一个由6个地质界面封闭得到的地层格架实施上述规则化就相当于将速度未知量的数目从一个很大的规模压缩为6个.同时,由于每一块内的速度更新量各有差异,因此更新后的模型在地质层位两侧必然会呈现出强反差.这正是地层格架约束立体层析所想要得到的地质一致性块状特征.一个关键问题是:地层格架如何在反演中的到正确的更新?注意地层格架是由许多地质层位封闭成的,而这些地质层位是由许多反射点构建的.而这些反射点的位置本身就是立体层析模型空间的一部分,在每一次迭代中它们都和速度、散射角等模型信息一起被更新.因此完全可以利用更新后的反射点位置来拟合出更新后的地质层位.这里我们使用非均匀B样条(NUBS)(2维中心坐标系下地表格架内射线扰动的一阶关系那么二维地层格架约束立体层析的Fréchet导数如何计算?注意二维地层格架约束立体层析的数据空间与传统二维立体层析并无差别,依然可以描述为d=(S炮点一侧的数据空间对于检波点一侧的散射角显然没有敏感度,反之亦然.同时根据费马原理,走时T在速度、炮点与接收点坐标都固定时已经是最小走时,因此散射角对走时亦不存在一阶扰动,有:根据Farra和Madariaga(1987),二维中心坐标系下立体层析所用的线性扰动关系如(4)、(5)式所示.(4)式等号左边为地表观测的坐标扰动(Δq(4)式中的Δq将(7)式代入(4)式容易推出:方括号中后两项为显然,在射线到达地表后,需要将扰动量从中心坐标系换算到地表水平观测面.因此需要建立射线中心坐标系下的位移扰动Δq将(6)、(7)、(8)、(9)、(10)式代入到(4)式中,可得出二维中心坐标系下地表观测数据信息与初始射线扰动扰动信息的一阶关系为如前所述,当速度模型被地层格架信息分块之后,速度更新量在一个地质块内将是一个常数,因此(12)式中括号里的第一项,即速度更新量关于空间坐标的导数就应该是零;其次一个地质块内的速度定义为根据链式法则,二维中心坐标系下坐标扰动Δq、慢度扰动Δp关于v最终我们得到:类似地,将链式法则应用于(5)式,就可以得到走时关于v至此,在引入地层格架信息约束后,走时t,地表观测坐标x,地表观测水平慢度矢量p3基于非均匀b样条的地质建模对地层格架约束立体层析而言,在反演实施之前对散乱反射点的拟合和自动封闭至关重要.这里我们采用了非均匀B样条拟合算法将散乱反射点拟合为连续的地质层位.相比均匀B样条拟合,该算法能够高效拟合多值函数,其稳健性在其他行业中已经得到了充分体现,证明完全能够适应二维、三维复杂地质建模的需求.这里对其实现原理做一简介,更多细节请读者参阅(Piegland3.1未知控制点求解B样条基函数的定义如下:令U={u而对于p次B样条曲线的定义则为对于数据拟合,需要预先计算好数据点的参数值和节点矢量,然后,建立并求解线性最小二乘问题来求解未知控制点.假定p≥1,n≥p,并且给定数据点Q其满足条件:Q求取拟合曲线主要分为以下几个步骤:(1)预先通过弦长参数化得到数据点Q由于节点的分布应该反映节点矢量的分布,假定c是一个正实数,用i=int(c)表示小于或等与d的最大整数,则总共需要n+p+2个节点,因此有n-p个内节点和n-p-1个内部节点区间,令c=(m+1)/(n-p+1),然后按公式(20)定义内节点,即:(2)建立目标函数,公式为3.2层位组合格架封闭第3.1节中的非均匀B样条拟合仅仅是针对某一个地质层位的离散反射点完成了连续地质层位的拟合.但是要实施地层格架的自动封闭必须考虑地质层位之间的相交或者令其自然延伸到剖面的边界处,这些在实施自动封闭时都需要非常仔细的考虑.此外,对于含有多个地质块体的地震剖面,地层的封闭不但需要考虑块体顶底面之间的关系,还要考虑块体与地层以及地层与地层之间的关系.这要求离散的反射点数据必须含有地质层位信息,即在反演之前,在自动拾取主要地质层位的立体层析数据空间时就需要加入层位编号.实施地层格架封闭分为以下5个步骤:(1)获取层位编号信息,获取块体与地层、地层与地层之间的几何拓扑关系,即确定块体与层位之间、层位与层位之间的拓扑关系(比如一个交点由哪些层位组成,一个块体由哪些层位包围).(2)通过反射点信息拟合各个层位:对每个层位的数据点根据横坐标进行排序后,然后对同一个层位的数据进行NUBS拟合,利用结构张量算法提取拟合后界面的切向量和法向量信息(Wuand(3)处理界面与边界之间的关系:在确定层位之间的关系时,统计每个层位的交点个数(只有0、1、2三种情况),当只有一个交点时,判断哪个端点需要延拓到边界,当没有交点时,地层两端根据当前的界面端点处的斜率延拓到边界,斜率信息是通过结构张量算法获得的.(4)处理界面之间的拓扑关系:对于每个交点,确定哪些界面相交于该点,然后通过一些处理使得这些层位相交于一点,处理过程中同样需要结构张量获取的界面的切向量和法向量信息.(5)根据块体与层位之间的拓扑关系求取块体顶底面的深度坐标,每一个块体内部一个横坐标的顶底面深度是实现分块求取Fréchet导数的重要参数.图4显示了基于离散的反射点实施地层格架封闭的全过程.图4a显示了块体与界面之间、界面与界面之间的拓扑关系;图4b显示了如何通过反射点信息拟合各个反射界面;图4c显示了将界面延伸到边界后的情况;图4d显示了最为重要也是相对复杂的一步:处理好界面与界面之间的交点关系.4基于离散反射点的连续界面拟合到目前为止我们已经实现了:(1)推导了二维地层格架约束立体层析的Fréchet导数;(2)利用NUBS插值实现基于离散反射点的连续界面拟合;(3)利用初始成像界面提供的地质块与界面的拓扑结构实现地层格架的自动封闭.根据上述这些算法可以得到一个实用的、适用于二维实际数据的工作流程.该流程实现过程如下:(1)对输入的叠前地震数据利用结构张量算法提取高密度立体层析数据空间[S5理论数据和真实数据的计算5.1初始速度模型的建立和反演的验证图6显示这个工作流程在图2所示的六层背斜模型理论数据上是如何工作.基于图2显示的模型,共正演了6000个射线对并得到对应的立体层析数据空间.首先运行一个传统立体层析得到初始反射点位置和初始速度模型(如图6a所示).图6b显示基于图6a中获得的初始反射点位置实施NUBS拟合可以获得6个连续的反射界面.图6c显示根据图6b所示反射界面实施地层格架封闭后的结果.地层格架约束立体层析的初始模型是将图6c所示的地层格架内充填一个初始速度后得到的.初始充填速度的获得方式和图5的建议流程完全一致,即采取在图6a所示的初始模型内每一个地质块内做一个简单的算术平均后得到的.图6d显示了地层格架约束二维立体层析在这个模型上的最终反演结果.我们可以看出图6d所示的反演结果和图2所示的真实模型已经非常接近,证明这个流程在图2所示的层状模型上是有效的.5.2基于深度成像的地层格架和控制律设计的优化模型图7显示这个工作流程在另一个更为复杂的理论数据上如何工作.基于图7a显示的模型,正演了一条2D理论测线,共激发201炮,炮间距为40m,每炮401道接收,道间距为10m,最大偏移距2km.在自动搜索到5700个数据点后,首先运行一个传统立体层析得到初始反射点位置和初始速度模型(如图7b).基于图7b显示的反射点位置通过NUBS拟合算法可以实现地层格架自动封闭(如图7d所示).图7e显示了基于图7b显示的初始模型和图7d所示的封闭后地层格架在每一块内做一个简单的平均后得到了用于二维地层格架约束立体层析的初始模型.在图7e基础上,我们运行了地层格架约束立体层析.图7f、g分别显示地层格架约束立体层析第1次和第10次迭代更新的速度模型.我们可以看出图7g所示的最终反演结果和图7a所示的真实模型已经非常接近.图7h为目标函数的下降曲线,可以看出最终的二范数残差已经非常小.图8a、b显示了用图7b所示传统立体层析模型实施偏移成像得到的共成像点道集(CIG)以及用图7g所示地层格架约束立体层析反演模型实施偏移成像获得的CIG之间的对比.图8c、d显示了用图7b所示传统立体层析模型实施偏移成像得到的成像剖面以及用图7g所示地层格架约束立体层析反演模型实施偏移成像获得的成像剖面的对比.可以看出,本文方法获得的深度成像剖面的聚焦程度和成像品质相比传统立体层析有明显提高.5.3基于数据信噪比的地层格架约束立体层析二维实际数据是2011年在南海深水区采集获得,共选取了该测线中950炮,炮间距为50m,道间距为25m,最大偏移距为8275m.考虑到实际模型的复杂性,每一层内的速度设置为随纵、横向线性变化.图9a显示了传统立体层析得到的最终反演模型,其中反演得到的倾角条用蓝色线画在速度模型之上.应用同样的工作流程,在NUBS拟合算法和结构张量算法的帮助下,我们得到了一个封闭的地质格架.由于该数据信噪比的原因,仅仅基于了几个关键层位建立了一个比较粗放的地层格架,没有针对更精细的反射点信息实施地层格架封闭.尽管如此,地层格架约束立体层析对于深部的速度结构反演依然有明显改善.图9b、c显示地层格架约束立体层析的第1次和第7次迭代结果.图8d显示迭代过程中目标函数的下降曲线.由于实际数据存在信噪比以及数据点提取的密度和精度问题,因此其目标函数的下降幅度不可能像理论数据的目标函数那样,可以下降到接近于零的水平.但是在引入地层格架之后,对于测线中部5000～8000m深度的内幕的成像依然有明显改进.图9e、f显示了用图9a所示传统立体层析模型实施偏移成像得到的共成像点道集(CIG)以及用图9c所示地层格架约束立体层析反演模型实施偏移成像获得的CIG之间的对比.图9g