一种基于辅助坐标系的起伏地表波形反演方法技术

本发明专利技术公开了一种基于辅助坐标系的起伏地表波形反演方法，属于石油地球物理勘探技术领域，本发明专利技术首先采用辅助坐标系下的早至波波形反演方法对曲网格区域进行反演，以得到更新的近地表速度，再采用辅助坐标系下的全波形反演全局速度场进行速度更新，以克服近地表速度不准确对深部速度场反演的影响；采用时间域多尺度反演方法从低频到高频进行速度反演，以克服波形反演方法对初始速度模型的依赖性。本发明专利技术能够很好地反演剧烈起伏地表的速度场，为高精度成像方法提供准确的偏移速度场。

【技术实现步骤摘要】
一种基于辅助坐标系的起伏地表波形反演方法
本专利技术属于石油地球物理勘探领域，具体涉及一种基于辅助坐标系的起伏地表波形反演方法。
技术介绍
剧烈起伏地表对地震资料采集、处理带来了严重的影响，地球物理研究人员发展了一系列方法来克服这个问题。目前，针对起伏地表的处理主要有两种策略：一是对表层波场进行校正，二是基于起伏地表进行深度偏移成像。但复杂条件下静校正量难以准确计算，且静校正无法彻底消除地表起伏对地震波场造成的扭曲，因此直接基于起伏地表的深度偏移成像方法逐渐成为研究的热点。由于全波形反演高精度以及高分辨率的特点，使其成为速度建模的一种有力工具，逐渐成为研究的热点。全波形反演是一个非线性数据拟合的过程，通过减少观测数据与预测数据的之间的差值来更新参数模型，这个过程以迭代的方式重复下去，直到数据差值足够小为止。
技术实现思路
针对现有技术中存在的上述技术问题，本专利技术提出了一种基于辅助坐标系的起伏地表波形反演方法，设计合理，具有良好的推广价值。为了实现上述目的，本专利技术采用如下技术方案：一种基于辅助坐标系的起伏地表波形反演方法，按照如下步骤进行：步骤1：输入初始全局速度场、常规叠前炮记录、起伏高程及震源子波，并建立观测系统；步骤2：根据初始全局速度场及起伏高程进行网格剖分，将近地表附近的网格剖分成曲网格，深层的网格剖分成矩形网格；步骤3：将初始全局速度场变换到辅助坐标系下的矩形网格，采用下式所示的变换格式：其中，x和z表示笛卡尔坐标系下的横纵坐标；ξ和η表示辅助坐标系下的横纵坐标；zi-1(ξ)和zi(ξ)是笛卡尔坐标系下第i层顶界面、底界面的高程，定义最深层高程为零；ηi-1(ξ)和ηi(ξ)是辅助坐标系第i层顶界面、底界面的纵向采样点数，定义最深层纵向采样点数为零；步骤4：对常规叠前炮记录划分时窗，在辅助坐标系下对近地表附近的曲网格区域速度场应用早至波波形反演，更新近地表的速度场，早至波波形反演的梯度方向如下：其中，g为梯度；xs为炮点坐标；v为介质速度；p为声压；p*是残差波场的反向传播；t为时间；Te表示早至波时窗；步骤5：判断使用近地表速度场正演模拟的炮记录与步骤4中划分时窗后的常规叠前炮记录之差是否满足误差条件；若：判断结果是使用近地表速度场正演模拟的炮记录与步骤4中划分时窗后的常规叠前炮记录之差满足误差条件，则近地表速度场更新完成，然后执行步骤6；或判断结果是使用近地表速度场正演模拟的炮记录与步骤4中划分时窗后的常规叠前炮记录之差不满足误差条件，则执行步骤4；步骤6：将步骤1中输入的常规叠前炮记录分解成不同主频的多尺度炮记录，采用如下式所示的分解公式：其中，Fw是维纳滤波器；Wo表示初始炮记录子波；Wt是生成的炮记录子波；ω是角频率；ε为一个小数；*表示共轭转置；步骤7：将更新过近地表速度的全局速度场作为初始速度场，在辅助坐标系下应用全波形反演从低频到高频更新全局速度场，全波形反演的梯度方向如下：其中，f表示主频；f1和fmax为多尺度分解的最低频率和最高频率；Tmax表示常规叠前炮记录的最大记录时间；表示主频是f的残差波场反传；步骤8：判断使用全局速度场正演模拟的炮记录与常规叠前炮记录之差是否满足误差条件；若：判断结果是使用全局速度场正演模拟的炮记录与常规叠前炮记录之差满足误差条件，则全局速度场更新完成，然后执行步骤9；或判断结果是使用全局速度场正演模拟的炮记录与常规叠前炮记录之差不满足误差条件，则执行步骤7；步骤9：将更新完成的全局速度场反变换到笛卡尔坐标系下，反变换公式如下：步骤10：输出反演的速度场。优选地，在步骤4中，具体包括步骤4.1：定义目标函数：其中，u(t,xr,xs)代表模拟波场u＝(u,w,p)T，其中T表示转置；R为限定算子；dobs(t,xr,xs)是常规叠前炮记录；xs和xr表示震源点和检波点的位置坐标；t表示时间；E为目标函数值；步骤4.2：将目标函数变分，得到变分表达式：步骤4.3：定义变换格式：通过变换格式(3)与锁链法则得到下面的映射公式：通过映射公式(4)，得到辅助坐标系下的一阶方程：其中，p是声压；u和w分别是水平方向和垂直方向的质点速度；v是介质速度；s表示震源项；ρ为密度；速度的扰动δv可以引起地震波场的扰动δu，δu＝(δu,δw,δp)T，将v+δv→u+δu代入一阶方程(5)后与一阶方程(5)相减得到如下方程：进一步可得：其中，L代表正演算子；步骤4.4：将方程(7)代入变分表达式(2)，可得：其中，L*R*(Ru-dobs)表示波场的逆时传播；利用伴随状态法，L*R*(Ru-dobs)可由下式求得：则L*R*(Ru-dobs)＝p*；步骤4.5：求取目标函数对速度模型的梯度，得到早至波波形反演的梯度方向：其中，g为梯度；xs为炮点坐标；v为介质速度；p为声压；p*是残差波场的反向传播；t为时间；Te表示早至波时窗。本专利技术所带来的有益技术效果：本专利技术提出了一种基于辅助坐标系的起伏地表波形反演方法，与现有技术相比，一种基于辅助坐标系的起伏地表波形反演方法，采用一阶方程变密度方程克服传统常密度二阶方程在处理密度变化比较大地区的速度反演不准确的缺点；同时在起伏地表处采用曲网格剖分，在深部区域采用矩形网格剖分，并统一变换到辅助坐标系下的矩形网格进行计算，克服起伏地表对地震波传播的影响；首先采用辅助坐标系下的早至波波形反演方法对曲网格区域进行反演，以得到更新的近地表速度，再采用辅助坐标系下的全波形反演全局速度场进行速度更新，以克服近地表速度不准确对深部速度场反演的影响；采用时间域多尺度反演方法从低频到高频进行速度反演，以克服波形反演方法对初始速度模型的依赖性。本专利技术能够很好地反演剧烈起伏地表的速度场，为高精度成像方法提供准确的偏移速度场。附图说明图1为本专利技术的一种基于辅助坐标系的起伏地表波形反演方法的流程图。图2为本专利技术使用的真实起伏地表速度模型。图3为变换到辅助坐标系下的速度模型。图4为采用统一曲网格变换到辅助坐标系下的速度模型。图5为初始全局速度场。图6为输入的常规叠前炮记录。图7为采用统一曲网格正演模拟得到的炮记录。图8为多尺度分解得到的5Hz叠前炮记录。图9为多尺度分解得到的15Hz叠前炮记录。图10为采用本专利技术方法得到的第一尺度反演速度场(主频为5Hz)。图11为采用本专利技术方法得到的第二尺度反演速度场(主频为15Hz)。图12为采用本专利技术方法得到的第三尺度反演速度场(主频为25Hz)。图13为采用单尺度波形反演方法得到的速度场。图14为采用统一曲网格得到的波形反演速度场。具体实施方式下面结合附图以及具体实施方式对本专利技术作进一步详细说明：一种基于辅助坐标系的起伏地表波形反演方法的流程图(如图1所示)，包括如下步骤：输入初始全局速度场、常规叠前炮记录、起伏高程及震源子波，并建立观测系统；根据初始全局速度场及起伏高程进行网格剖分，近地表附近的网格剖分成曲网格，深层的网格剖分成矩形网格；将初始全局速度场变换到辅助坐标系下的矩形网格；对常规叠前炮记录划分时窗，在辅助坐标系下对近地表附近的曲网格区域速度场应用早至波波形反演，更新近地表的速度场，判断使用近地表速度场正演模拟的炮记录与划分时窗的常规叠前炮记录的差是否满足误差条件，如果满足则近地表速度场更新本文档来自技高网...

【技术保护点】
一种基于辅助坐标系的起伏地表波形反演方法，其特征在于：按照如下步骤进行：步骤1：输入初始全局速度场、常规叠前炮记录、起伏高程及震源子波，并建立观测系统；步骤2：根据初始全局速度场及起伏高程进行网格剖分，将近地表附近的网格剖分成曲网格，深层的网格剖分成矩形网格；步骤3：将初始全局速度场变换到辅助坐标系下的矩形网格，采用下式所示的变换格式：{x(ξ,η)=ξz(ξ,η)=zi-1(ξ)-zi(ξ)ηi-1(ξ)-ηi(ξ)(η-ηi)+zi(ξ);]]>其中，x和z表示笛卡尔坐标系下的横纵坐标；ξ和η表示辅助坐标系下的横纵坐标；zi‑1(ξ)和zi(ξ)是笛卡尔坐标系下第i层顶界面、底界面的高程，定义最深层为零；ηi‑1(ξ)和ηi(ξ)是辅助坐标系第i层顶界面、底界面的纵向采样点数，定义最深层为零；步骤4：对常规叠前炮记录划分时窗，在辅助坐标系下对近地表附近的曲网格区域速度场应用早至波波形反演，更新近地表的速度场，早至波波形反演的梯度方向如下：g(v)=2v3ΣxsΣt=0Te∂p∂t·p*;]]>其中，g为梯度；xs为炮点坐标；v为速度场的速度；p为声压；p*是残差波场的反向传播；t为时间；Te表示早至波时窗；步骤5：判断使用近地表速度场正演模拟的炮记录与步骤4中划分时窗后的常规叠前炮记录之差是否满足误差条件；若：判断结果是使用近地表速度场正演模拟的炮记录与步骤4中划分时窗后的常规叠前炮记录之差满足误差条件，则近地表速度场更新完成，然后执行步骤6；或判断结果是使用近地表速度场正演模拟的炮记录与步骤4中划分时窗后的常规叠前炮记录之差不满足误差条件，则执行步骤4；步骤6：将步骤1中输入的常规叠前炮记录分解成不同主频的多尺度炮记录，采用如下式所示的分解公式：f(ω)=Wt(ω)Wo*(ω)|Wo(ω)|2+ϵ2]]>其中，f是维纳滤波器；Wo表示初始炮记录子波；Wt是生成的炮记录子波；ω是角频率；ε为一个小数；*表示共轭转置；步骤7：将更新过近地表速度的全局速度场作为初始速度场，在辅助坐标系下应用全波形反演从低频到高频更新全局速度场，全波形反演的梯度方向如下：g(v)=2v3ΣxsΣt=0TmaxΣf=f1fmax∂p∂t·pf*]]>其中，f表示主频；f1和fmax为多尺度分解的最低频率和最高频率；Tmax表示常规叠前炮记录的最大记录时间；表示主频是f的残差波场反传；步骤8：判断使用全局速度场正演模拟的炮记录与常规叠前炮记录之差是否满足误差条件；若：判断结果是使用全局速度场正演模拟的炮记录与常规叠前炮记录之差满足误差条件，则全局速度场更新完成，然后执行步骤9；或判断结果是使用全局速度场正演模拟的炮记录与常规叠前炮记录之差不满足误差条件，则执行步骤7；步骤9：将更新完成的全局速度场反变换到笛卡尔坐标系下，反变换公式如下：ξ(x,z)=xη(x,z)=ηi-1(ξ)-ηi(ξ)zi-1(ξ)-zi(ξ)(z-zi)+ηi(ξ);]]>步骤10：输出反演的速度场。...

【技术特征摘要】
1.一种基于辅助坐标系的起伏地表波形反演方法，其特征在于：按照如下步骤进行：步骤1：输入初始全局速度场、常规叠前炮记录、起伏高程及震源子波，并建立观测系统；步骤2：根据初始全局速度场及起伏高程进行网格剖分，将近地表附近的网格剖分成曲网格，深层的网格剖分成矩形网格；步骤3：将初始全局速度场变换到辅助坐标系下的矩形网格，采用下式所示的变换格式：其中，x和z表示笛卡尔坐标系下的横纵坐标；ξ和η表示辅助坐标系下的横纵坐标；zi-1(ξ)和zi(ξ)是笛卡尔坐标系下第i层顶界面、底界面的高程，定义最深层高程为零；ηi-1(ξ)和ηi(ξ)是辅助坐标系第i层顶界面、底界面的纵向采样点数，定义最深层纵向采样点数为零；步骤4：对常规叠前炮记录划分时窗，在辅助坐标系下对近地表附近的曲网格区域速度场应用早至波波形反演，更新近地表的速度场，早至波波形反演的梯度方向如下：其中，g为梯度；xs为炮点坐标；v为介质速度；p为声压；p*是残差波场的反向传播；t为时间；Te表示早至波时窗；步骤5：判断使用近地表速度场正演模拟的炮记录与步骤4中划分时窗后的常规叠前炮记录之差是否满足误差条件；若：判断结果是使用近地表速度场正演模拟的炮记录与步骤4中划分时窗后的常规叠前炮记录之差满足误差条件，则近地表速度场更新完成，然后执行步骤6；或判断结果是使用近地表速度场正演模拟的炮记录与步骤4中划分时窗后的常规叠前炮记录之差不满足误差条件，则执行步骤4；步骤6：将步骤1中输入的常规叠前炮记录分解成不同主频的多尺度炮记录，采用如下式所示的分解公式：其中，Fw是维纳滤波器；Wo表示初始炮记录子波；Wt是生成的炮记录子波；ω是角频率；ε为一个小数；*表示共轭转置；步骤7：将更新过近地表速度的全局速度场作为初始速度场，在辅助坐标系下应用全波形反演从低频到高频更新全局速度场，全波形反演的梯度方向如下：其中，f表示主频；f1和fmax为多尺度分解的最低频率和最高频率；Tmax表示常规叠前炮记录的最大记录时间；表示主频是f的残差波场反传；步骤8：判断使用全局速度场正演模拟的炮记录与常规叠前炮记录之差是否满足误差条件；若：判断结果是使用全局速度场正演模拟的炮记录与常规叠前炮记录之差满足误差条件，则全局速度场更新完成，然后执行步骤9；或判断结果是使用全局速度场正演模拟的炮记录与常规叠前炮记录之差不满足误差条件，则执行步骤7；步骤9：将更新完成的全局速度场反变换到笛卡尔坐标系下，反变换公式如下：步骤10：输出反演的速度场。2.根据权利要求1所述的基于辅助坐标系的起伏地表波形反演方法，其特征在于：在步骤4中，具体包括步骤4.1：定义目标函数：其中，u(t,xr,xs)代表模拟波场u＝(u,w,p)T，其中T表示转置；R为限定算子；dobs(t,xr,xs)是常规叠前炮记录；xs和xr表示震源点和检波点的位置坐标；t表示时间；E为目标函数值；步骤4.2：将目标函数变分，得到变分表达式：步骤4.3：定义变换格式：通过变换格式(3)与锁链法则得到下面的映射公式：通过映射公式(4)，得到辅助坐标系下的一阶方程：

【专利技术属性】
技术研发人员：曲英铭，李振春，李金丽，黄建平，
申请(专利权)人：中国石油大学华东，曲英铭，
类型：发明
国别省市：山东;37