地震正演与高斯束偏移

地震成像与速度反演李振春中国石油大学（华东）地球科学与技术学院提纲一、双变网格正演模拟二、高斯束偏移三、波动方程偏移四、绕射波成像五、多次波成像六、层析速度分析与全波形反演双变网格正演模拟算法主要包括：基于声波、弹性波、各向异性介质的2D/3D有限差分算法，基于粘声、粘弹介质的2D正演模拟算法。适用于起伏地表模拟的变网格算法、2D谱元素法及伪谱法；适用于裂缝型介质的双变网格正演模拟算法；以及较为完善的边界处理方法包括：吸收衰减边界条件、PML以及CPML边界条件等。2D3D声波粘弹性/各向异性弹性介质模型传统有限差分地震波模拟方法基于规则的矩形网格，是处理非均匀介质，特别是连续非均匀介质波动问题的有效方法！双变网格正演模拟原理

难题2难题1为保证计算精度和计算稳定性，需要精细网格剖分，导致计算量增加和局部过采样需要发展等效模型，将非均质体的物性参数与波动方程联系起来速度场离散方式1速度场离散方式2双变网格算法即根据地震地质模型将速度场进行不同尺度的网格剖分，即变网格，同时又结合局部变时间算法，提高精度和效率。模型网格划分的优劣与波场正演的效率密切相关。

速度场多尺度离散策略Virieux（1986）和Graves（1996b）提出的一阶速度应力方程，各向同性非均匀介质方程可写为：变网格弹性波方程离散化原理在传统交错网格技术中，变量关于的2n阶空间偏导数差分近似式为交错网格高阶差分格式

差分系数是常数差分算子是固定步长弹性波方程离散化在空间变网格技术中，变量关于的2n阶空间偏导数差分近似式为变交错网格高阶差分格式

差分系数是网格步长的函数差分算子是空间变化的算法优点：

在全局时间计算时刻粗细网格过渡区域不需要重复设置计算区域空间的差分精度不需要有任何改变只需要前期求解出差分系数，便可以在每一计算时刻引用弹性波方程离散化

空间网格的变化需要时间采样满足精细网格采样对应的稳定性条件！全局时间采样不需要满足最短波长原则，同传统交错网格算法的稳定性条件是一致的。但是在全局时间层内，局部精细网格的时间计算步长需要满足最短波长原则。双变网格算法局部可变的时间采样步长首先给出一些假设：已知时刻的波场值已知时刻的波场值

已知时刻的波场值已知时刻的波场值

具体实现步骤如下：1：粗网格平面的求解时刻的波场值时刻的波场值2：利用过渡区域的波场值判断是否进入精细时间采样局部时间采样变化（LVTS）首先给出一些假设：已知时刻的波场值已知时刻的波场值

已知时刻的波场值已知时刻的波场值

具体实现步骤如下：3：精细时间采样平面内的精细网格点波场值求解时刻的波场值利用双线性公式：获得精细时间采样平面边界处的波场值4：波场传递局部时间采样变化（LVTS）速度场网格分为粗糙网格、精细网格，中间存在过渡区域时间层波场分为全局时间层波场和精细时间层波场总波场满足弹性波波动方程时间层局部波场求解空间精细网格剖分变交错网格求解时间层全局波场求解空间粗糙网格剖分传统交错网格求解波场传递总波场算法流程模型参数：水平距离：401*6m纵向深度：400*6m精细网格采样：2m全局时间采样：0.3ms局部时间采样：0.1ms炮点位置：（201*6m，5*6m）模拟时间长度：2.5s传统算法x分量双变算法x分量传统算法z分量双变算法z分量典型模型一：低降速带模型2、双变网格模型试算

对比可知：基于相同的常网格大小，而对目标低速区域采用精细网格剖分的双变网格算法，能够较好的压制频散，得到振幅连续性更好的同相轴分布。速度场模型典型模型二：局部低速体绿色框内表示双变区域模型参数：水平距离：930*6m纵向深度：550*6m精细网格采样：2m全局时间采样：0.3ms局部时间采样：0.06ms炮点位置：（101*6m，3*6m）模拟时间长度：2.4s左图双变算法波场快照中图常网格波场快照右图二者波场快照之差

波场快照对比图

由于采用了局部精细化的双变网格技术，在局部目标区域网格可以从米级变化到毫米级（而在全局采用毫米级的常网格算法将导致巨大的计算量，因而是不可行的），从波场快照的对比图可以看出：双变网格算法能够较好模拟局部小尺度低速构造。Z分量X分量左图双变算法波场快照中图常网格波场快照右图二者波场快照之差

炮记录对比图

对比可知：两种网格剖分方式的炮记录差异主要来自于局部小尺度的低速异常构造，为实现高精度的正演模拟及偏移成像，基于变网格剖分的双变网格算法，具有较好的应用前景。Z分量X分量双变网格比精细网格、精细采样节省80.1%的时间比局部精细网格，全局采样节省60.3%的时间从计算耗时对比可知：在保证计算精度的前提下，模型网格非均匀剖分的双变网格算法，能够极大的提高正演模拟的效率。双变区域模型参数：水平距离：601*8m纵向深度：700*6m精细网格采样：2m全局时间采样：0.5ms局部时间采样：0.16ms炮点位置：（301*8m，10*6m）模拟时间长度：3.5s典型模型三：陡坡带砂砾岩体薄层图示速度场最薄层厚30-60m左图双变算法波场快照中图常网格波场快照右图二者波场快照之差

波场快照对比图Z分量X分量双变网格算法能够较好的模拟出来自陡坡带砂砾岩体的地震波响应，振幅更加均匀、同相轴一致性更好。左图双变炮记录中图常网格炮记录右图二者炮记录之差Z分量X分量双变网格算法能够展示更多的来自于地下细微构造的波场信息，局部分辨率较常网格大大提高。左图全局速度场右图含裂缝局部速度场图示含裂缝速度场裂缝参数：裂缝位于3000m深度，密度10m/条，共200个裂缝，缝宽2mm典型模型五：裂缝发育层左图无缝炮记录中图裂缝延伸150m炮记录右图裂缝对应某时刻波场快照裂缝顶部裂缝底部上界面反射P波上界面反射S波裂缝的波场响应从波场快照及单炮记录可以看出：双变网格高阶有限差分算法，能够较好的模拟毫米级（2mm）的裂缝型储层，在裂缝密度达到一定程度时，裂缝地震响应特征较为清晰。基于双变网格算法的小结1、双变网格有限差分方法，依据模型特征采用非均匀网格剖分，可以在不显著增加计算成本的基础上，极大提高目标研究区域的正演模拟精度并抑制计算频散；2、双变网格有限差分方法，能够较好的对毫米级裂缝进行正演模拟，裂缝波场响应特征清晰；3、双变网格有限差分方法，较容易与物理模拟方法结合，并为物理模拟提供技术支持；4、双变网格有限差分方法，可以进一步发展为双变网格的RTM，极大提高RTM的成像效率，为RTM的生产实用化提供帮助。提纲一、双变网格正演模拟二、高斯束偏移三、波动方程偏移四、绕射波成像五、多次波成像六、层析速度分析与全波形反演地震偏移成像方法高斯束偏移成像波动方程的数值解法单程波动方程偏移双程波动方程偏移几何射线类方法Kirchhoff偏移（复杂构造难以成像）高斯束偏移（高效灵活精确）缺陷：计算量巨大主要内容射线类偏移高斯束偏移★高斯束的基本性质★高斯束偏移成像技术★高斯束偏移共成像点道集

高斯束逆时偏移★声波高斯束逆时偏移★弹性波高斯束逆时偏移高斯束偏移成像垂向传播的高斯束为虚数代表空间局部化、方向局部化的地震波场1、高斯束的基本性质高斯束偏移成像高斯束偏移成像1、高斯束的基本性质高斯束表示球面波高斯束表示平面波高斯束偏移成像高斯束偏移发展历程GaussianbeammigrationRossHill,1990,GeophysicsPrestackGaussian-beamdepthmigrationRossHill,2001,GeophysicsGaussianbeammigrationofcommonshotrecordingSamuelGray,2005,GeophysicsTrueamplitudeGaussianbeammigrationSamuelGray,NormanBleistein,2009,Geophysics高斯束偏移成像2、高斯束偏移基本过程1）不同道集成像公式2）保真振幅高斯束偏移3）复杂地表高斯束偏移4）弹性波多分量高斯束偏移5）3D高斯束偏移高斯束偏移成像高斯束叠后偏移加窗局部倾斜叠加叠后偏移高斯束偏移成像高斯束叠后偏移洼陷模型叠后记录

30个波束最终成像结果高斯束偏移成像常规高斯束偏移计算流程震源成像点高斯束叠前偏移高斯束偏移成像不同道集成像过程共炮道集成像公式共接收点道集成像公式高斯束偏移成像不同道集成像过程共偏移距道集成像公式共偏移距道集倾斜叠加高斯束偏移成像不同道集成像过程共中心点道集成像公式共中心点道集倾斜叠加高斯束偏移成像共偏移距道集高斯束偏移Marmousimodel共炮集高斯束偏移Kirchhoff偏移不同道集成像过程高斯束偏移成像对于每个束中心位置{

加窗局部倾斜叠加并计算由震源和接收点出射的高斯束对于每个震源射线参数{

对于每个接收点射线参数{

对于局部孔径内的每个成像点{

根据高斯束的走时以及振幅利用所对应的倾斜叠加道进行成像

}成像点循环结束

}接收点高斯束循环结束

}震源高斯束循环结束}束中心循环结束共炮集偏移局部孔径的成像运算偏移计算量为：束中心为nb，震源高斯束为nps,接收点高斯束为npr全波至算法计算流程全波至算法高效算法高斯束偏移成像对于每个束中心位置{

加窗局部倾斜叠加并计算由震源和接收点出射的高斯束对于每个中心点射线参数{

对于每个粗网格点{

扫描并确定使虚值走时最小的偏移距射线参数

}粗网格循环结束对于孔径内的每个成像点{

根据高斯束的走时振幅信息利用所对应的倾斜叠加道进行成像

}成像点循环结束

}中心点射线参数循环结束}束中心循环结束局部孔径的成像运算偏移计算量为：高频复值振荡积分最速下降法高效算法计算流程(Hill,Geophysics,2001)算法比较高斯束偏移成像Sigsbee2amodel全波至算法（182.73s）高效算法（51.03s）算法比较高斯束偏移成像全波至算法高效算法常规Kirchhoff偏移单程波偏移高斯束偏移成像2、高斯束偏移基本过程1）不同道集成像公式2）保真振幅高斯束偏移3）复杂地表高斯束偏移4）弹性波多分量高斯束偏移5）3D高斯束偏移保幅高斯束偏移高斯束偏移成像保幅偏移定义

偏移过程中消除地震波几何扩散造成的能量损失，使得成像振幅正比于界面的反射系数。保幅偏移应用条件

保幅成像条件+保幅的延拓算子。反褶积成像条件(Claerbout,Geophysics1971)互相关成像条件(Claerbout,Geophysics1971)高斯束偏移成像炮域保幅角度域保幅（SamuelGray,NormanBleistein,Geophysics,2009)保幅高斯束偏移高斯束偏移成像单炮记录提取的反射系数四层模型（常速变密度）保幅高斯束偏移保幅高斯束偏移高斯束偏移成像保幅高斯束偏移单炮成像结果（互相关）角度域共成像道集（互相关）归一化反射系数归一化反射系数高斯束偏移成像保幅高斯束偏移2D盐丘模型常规高斯束偏移保幅高斯束偏移（反褶积）高斯束偏移成像保幅高斯束偏移常规高斯束偏移保幅高斯束偏移（互相关）CDP=350CDP=450CDP=530CDP=640高斯束偏移成像30~40°20~30°0~10°10~20°不同角度范围的叠加成像结果保幅高斯束偏移高斯束偏移成像2、高斯束偏移基本过程1）不同道集成像公式2）保真振幅高斯束偏移3）复杂地表高斯束偏移4）弹性波多分量高斯束偏移5）3D高斯束偏移复杂地表条件下地震成像方法波动方程基准面校正静校正法波场沿拓法射线类偏移方法准确性差计算量大，灵活性差，对数据要求高灵活，高效复杂地表高斯束偏移高斯束偏移成像高斯束偏移成像局部静校正法基本原理（Gray，Geophysics，2005）复杂地表高斯束偏移RayleighII积分（岳玉波等，地球物理学报，2011第三期）复杂地表保幅成像公式起伏地表平面波分解

起伏地表模型单炮记录复杂地表高斯束偏移高斯束偏移成像沿反射层提取的归一化振幅曲线保幅高斯束偏移常规高斯束偏移复杂地表高斯束偏移高斯束偏移成像加拿大逆掩断层模型单炮记录地表高程高斯束偏移成像复杂地表高斯束偏移Kirchhoff偏移局部静校正法保幅延拓法“直接下延”法波动方程偏移高斯束偏移成像复杂地表高斯束偏移速度模型地表高程保幅延拓法“直接下延”法波动方程偏移复杂地表高斯束偏移高斯束偏移成像高斯束偏移成像2、高斯束偏移基本过程1）不同道集成像公式2）保真振幅高斯束偏移3）复杂地表高斯束偏移4）弹性波多分量高斯束偏移5）3D高斯束偏移（1）弹性波Kirchhoff偏移具有极高的计算效率；成像噪声压制；容易解决转换波极性问题。精度不高，不能对多次波至进行成像。（2）弹性波逆时偏移具有极高的成像精度；计算量巨大；难以压制不同波型的噪声串扰；难以解决转换波极性问题。（3）弹性波高斯束偏移兼具较高成像精度与计算效率。弹性波高斯束偏移高斯束偏移成像公式推导弹性波Kirchhoff-Helmohltz积分弹性动力学高斯束表示格林张量P波位移S波位移权系数高斯束偏移成像弹性波高斯束偏移P-P成像公式P-S成像公式沿不同波型极性方向投影互相关成像条件求取标量成像值高斯束偏移成像弹性波高斯束偏移

速度模型（纵横波速度比:1.5/1）垂直分量水平分量高斯束偏移成像弹性波高斯束偏移单分量高斯束偏移弹性波高斯束偏移P-SP-SP-PP-PZ分量X分量高斯束偏移成像弹性波高斯束偏移(纵横波速度比为：1.7/1.0)斜层模型高斯束偏移成像弹性波高斯束偏移弹性波高斯束偏移P垂直分量水平分量S高斯束偏移成像弹性波高斯束偏移弹性波高斯束偏移P-SP-P高斯束偏移成像弹性波高斯束偏移高斯束偏移成像2、高斯束偏移基本过程1）不同道集成像公式2）保真振幅高斯束偏移3）复杂地表高斯束偏移4）弹性波多分量高斯束偏移5）3D高斯束偏移SEG/EAGEC3-NA模型速度模型观测系统第1300炮记录模型大小：676*676*210截取：451*481*210数据大小：6.25G偏移孔径：Inline2000m,Crossline1000m计算时间：100个计算节点，炮域并行算法，30h高斯束偏移成像3D高斯束偏移高斯束偏移成像3D高斯束偏移高斯束偏移成像3D高斯束偏移Inline337常规Kirchhoff偏移高斯束偏移高斯束偏移成像3D高斯束偏移Depth=2080m常规Kirchhoff偏移高斯束偏移高斯束偏移成像3D高斯束偏移高斯束偏移成像3、高斯束偏移共成像道集1）偏移距域共成像点道集2）角度域共成像点道集二维反射角域共成像点道集三维方位角/反射角域共成像点道集多分量反射角域共成像点道集3）平面波域共成像点道集

映射转换法空移成像条件时移成像条件L-ODCIGsTSCIGs高斯束偏移成像角道集的提取直接提取法走时抛物展开水平及垂直慢度传播角度高斯束偏移成像角道集的提取传播方向矢量3D反射张角3D反射方位角高斯束偏移成像角道集的提取P、S波传播角度Snell定律P波反射张角S波反射张角高斯束偏移成像角道集的提取洼陷模型L-ODCIGsTSCIGs成像结果高斯束偏移成像角道集的提取空移成像条件时移成像条件直接提取高斯束偏移成像角道集的提取速度模型成像结果CDP=500CDP=850高斯束偏移成像角道集的提取28°21°49°15°22°38°Azimuth=0°Azimuth=135°ADCIGs高斯束偏移成像角道集的提取Inline337成像结果ADCIGs高斯束偏移成像角道集的提取胜利模型P-P成像结果P-S成像结果PSPSPSPSCDP=300CDP=500CDP=700高斯束偏移成像角道集的提取高斯束偏移成像3、高斯束偏移共成像道集1）偏移距域共成像点道集2）角度域共成像点道集二维反射角域共成像点道集三维方位角/反射角域共成像点道集多分量反射角域共成像点道集3）平面波域共成像点道集高斯束偏移平面波偏移高斯束偏移成像平面波共成像道集CDP=300CDP=400CDP=500CDP=600Marmousi成像结果高斯束偏移成像平面波共成像道集不同平面波成像结果-15~0°0~15°15~30°30~45°高斯束偏移成像平面波共成像道集主要内容高斯束偏移成像1、声波高斯束逆时偏移2、弹性波高斯束逆时偏移声波高斯束逆时偏移高斯束偏移成像对于边界，取地面水平面，边界的侧面忽略(如下图)，Kirchhoff方程中的格林函数由高斯束渐近形式代替，其满足Kirchhoff近似中地震边界条件。结果有以下近似方程，即偏移域内一点时间间隔为的反向延拓波场，地下一点，时间时的波场为：上行波场声波高斯束逆时偏移高斯束偏移成像下行波场由地面一点源产生的直达波场