基于自适应优化有限差分的全波场vsp逆时偏移成像

基于自适应优化有限差分的全波场vsp逆时偏移成像

1有限差分正演算法近年来，随着vsp技术的发展和地震勘探对象的日益复杂，vsp图像可以适应波浪表面、陡峭结构和复杂速度区域的成像。由于地球地震数据的处理通常无法提供全面可靠的目标层岩丘下方的结构数据，vsp技术在澄清井侧结构、探索岩丘下方结构、提取纵向和横向波信息方面发挥着重要作用。VSP逆时偏移中四个关键问题为波场延拓、吸收边界条件、成像条件以及低频噪声的压制．波场延拓的实质是求解波动方程，由于有限差分方法具有计算效率高、精度较高、稳定性较好的优点，因而被广泛应用于求解波动方程．目前主要有两种基于频散关系计算空间导数有限差分系数的方法，一种是基于泰勒级数展开的方法，另一种是基于优化的方法．这两种方法都有一定的优越性和局限性，泰勒级数展开法求解有限差分系数的计算效率高，但是只能在较小的波数或者频率范围内达到较高的精度．对于优化的方法，局部优化方法得到的差分系数难以达到全局优化，而全局优化方法计算效率一般较低，且不一定能得到全局最优解．目前，VSP逆时偏移主要以上行反射波作为有效波进行逆时延拓，但仅用上行反射波做逆时延拓会导致偏移成像结果照明度低、成像范围窄（2方法原理2.1优化有限差分方法、时空域多元振动方程的频散曲线对比其中，p为标量波场，v为速度，t为时间，x、z为空间坐标．VSP逆时偏移的震源正向延拓和地震记录逆时延拓的表达式分别为其中，p其中，时空域泰勒级数展开有限差分方法是以平面波理论和泰勒级数展开为基础，基于时空域频散关系推导的，其差分系数可以通过求解下面的方程组得到（二维声波波动方程数值模拟的相速度频散δ表示为其中，v当δ图1是优化有限差分方法（optimal-basedFDmethod）、空间域泰勒级数展开有限差分方法（spaceTE-basedFDmethod）和时空域泰勒级数展开有限差分方法（time-spaceTE-basedFDmethod）的频散曲线对比图．由图可见：优化有限差分方法在区间[0，π]内，总体比空间域泰勒级数展开有限差分方法和时空域泰勒级数展开有限差分方法的数值频散都小，并且空间域泰勒级数展开有限差分方法和时空域泰勒级数展开有限差分方法的频散相近．2.2自适应变空间差分算子长度方案在求解波动方程时，通常采用固定的算子长度计算空间导数，而自适应变空间差分算子长度方案采用变化的算子长度计算空间导数．自适应变差分算子长度方案是基于算法的精度控制不同速度网格点的M值，其表达式为（其中，ε为单位网格内传播时间的相对误差，f为验证自适应变空间差分算子长度方案的有效性，设计一个速度模型，其速度的变化范围为1500～5000m·s(1）对于相同的η，速度越小，M值越大；(2）对于优化有限差分方法，η越小，则精度越高，M越大；(3）当η=102.3平滑过渡区的基本原理在VSP逆时偏移中，一个不可避免的问题是如何有效压制由计算区域边界引起的边界反射能量，本文采用混合吸收边界条件进行压制．混合吸收边界条件的基本原理是把计算区域分为内部区、过渡区、边界区．第一步，采用双程波波动方程计算内部区和过渡区波场值；第二步，采用单程波波动方程计算边界区与过渡区波场值；第三步，采用双程波波场值与单程波波场值的线性加权得到过渡区最终波场值．过渡区起到了对波场平滑过渡的作用，避免由内部区的双程波波动方程到边界区的单程波波动方程的急剧变化而导致的边界干扰反射无法得到有效压制的问题．以模型的上边界和左上角为例，其单程波波动方程表达式分别为（ClaytonandEngquist,1997;LiuandSen,2010)当过渡区的网格厚度为1时，混合吸收边界条件等效于Clayton-Engquist吸收边界条件；过渡区网格厚度为10时能达到较好的吸收效果．以重新采样后的Marmousi速度模型（图3）为例，速度模型大小为5000m×3500m，震源是位于（2500m,0m）的主频为30Hz的Ricker子波，h=10m，τ=1ms，记录时间为4s．图4为Marmousi速度模型的无边界条件与边界网格点数为10的混合吸收边界条件模拟地震记录，由图可见混合吸收边界条件能有效地压制边界反射．2.4震源归一化零延迟互相关成像条件t、pb模型震源归一化零延迟互相关成像条件的表达式为(ChattopadhyayandMcmechan,2008)其中，S(x,z,t）表示震源波场，R(x,z,t）表示检波点波场，T反褶积成像条件的表达式为（其中，S(x,z，ω）表示震源波场，珚S(x,z，ω）表示S(x,z，ω）的共轭，R(x,z，ω）表示检波点波场，R珚（x,z，ω）表示R(x,z，ω）的共轭，ω为频率．为对比震源归一化零延迟互相关成像条件与反褶积成像条件的成像效果，设计一个双层层状模型，模型大小为1000m×1000m,h=10m，τ=1ms，震源为主频为30Hz的Ricker子波，混合吸收边界过渡区网格点数为10，自适应变空间算子长度方案参数为f3模型计算3.1优化有限差分方法验证为了验证自适应优化有限差分方法的可行性，设计一个倾斜界面模型如图6a所示，模型大小为4000m×4000m,h=20m，τ=1ms，数值模拟的震源采用主频为20Hz的Ricker子波，且位于（2000m,0m），混合吸收边界过渡区网格点数为10．分别采用泰勒级数展开有限差分法与优化有限差分方法实现波场模拟．图6b—d为1.2s时刻数值模拟的波场快照，可见图6b的频散比较大，图6c与d频散相近，其结果表明：(1）对于相同算子长度（M=8），优化有限差分方法（图6d）模拟精度高于泰勒级数展开有限差分方法（图6b);(2)M=8时的优化有限差分方法（图6d）与M=32时的泰勒级数展开有限差分方法（图6c）模拟精度相近．倾斜界面模型试算结果表明：在相同算子长度情况下，优化有限差分方法模拟精度高于泰勒级数展开有限差分方法模拟精度，表明了优化有限差分方法适用于求解波动方程空间差分系数，同时试算结果也表明混合吸收边界条件能有效地压制边界反射．3.2反偏移vsp3.2.1逆时偏移试验设计一个层状模型，其速度模型如图7所示，模型大小为1000m×1200m,h=10m，τ=1ms，记录时间为2s，地表为自由边界，震源是主频为30Hz的Ricker子波．观测系统为：炮点位于（500m,0m）处，检波点位于（0m,400m）处．图8为检波点处地震记录，地震记录中可见直达波、一次反射波、一阶多次波、二阶多次波以及高阶多次波．分别对这5种波进行逆时偏移成像，偏移结果如图9所示，其中图9中深度z=650m处的黑线为界面位置．图9表明：直达波只会造成炮点到检波点路径上的噪声（图9a）；一次反射波在R1点处成像（图9b）；一阶多次波在R2点处成像，与R1点相比，R2点距离井的位置更远（图9c），所以多次波可以拓宽成像区域；二阶多次波可以贡献两个成像点（图9d），一个比R1点距井源更近，一个比R2点距井源更加远，说明了多次波可以增加照明度及拓宽成像范围；高阶多次波的成像范围比R1点到R2点相距范围更广（图9e中箭头所示）．以上5种波中，直达波和一阶多次波为下行波；一次反射波和二阶多次波为上行波；高阶多次波中既有上行也有下行波，如图8所示，其能量比一次反射波能量弱．若只用上行反射波作为逆时偏移的有效波场则会损失部分有效信息．理论分析表明：除直达波以外的所有波场对界面的成像均有贡献，所以多次波可以作为有效波进行逆时延拓，以增加照明度、拓宽偏移成像范围．3.2.2井旁层模型比较首先，我们针对层状介质模型分别实现自适应优化有限差分方法和自适应泰勒级数展开有限差分方法的VSP逆时偏移，以验证自适应优化有限差分方法的有效性和效率．其次，我们采用自适应优化有限差分方法完成三个模型试算，分别用分离后的上行波、分离后的下行波以及分离前的全波实现VSP逆时偏移，下文中逆时偏移采用的波场均不包括直达波．虽然上行波包括一次反射波以及多次波，但是直达波除外的下行波都是多次波．我们采用下行波偏移即多次波偏移，偏移效果可以证明多次波偏移的有效性，而上行波偏移结果与全波偏移结果的对比也可以验证全波偏移的优越性．(I）层状模型模型如图10a所示，模型大小为4000m×4000m,h=20m，τ=1ms，震源为主频为20Hz的Ricker子波，记录长度为2.5s．混合吸收边界过渡区网格点数为10．自适应变空间算子长度方案参数为f我们设计另一个层状模型（图11）分析多次波在VSP逆时偏移中的影响．模型大小为2000m×3000m,h=10m，τ=1ms，记录总时长为2.5s，震源是主频为30Hz的Ricker子波，混合吸收边界过渡区网格数为10．自适应变空间算子方案参数为f(1）上行波作为有效波时（图13a），成像位置准确，井旁界面成像清晰；(2）下行波作为有效波时（图13b），远井源距的界面成像比上行波偏移成像结果更清晰，而井旁界面成像不如上行波偏移好，成像结果有一定噪声，但其噪声能量较弱；(3）全波作为有效波（图13c）时，井旁界面和远井源距界面均准确成像．由于下行波成像的噪声能量与上行波的成像能量差异大，在全波成像结果中几乎看不到成像噪声．此层状模型试算结果表明，全波作为有效波时成像效果比仅用上行波或者仅用下行波作为有效波的成像效果更好．(II）断层模型断层速度模型（图14）大小为2000m×2000m,h=10m，τ=1ms，记录总时间为2s，震源是主频为30Hz的Ricker子波，混合吸收边界过渡区网格数为10，自适应变空间算子长度方案参数为f(1）上行波、下行波、全波波场分别作为有效波时，VSP逆时偏移成像位置都准确；(2）上行波作为有效波时（图15a），井旁界面成像较清晰；(3）下行波作为有效波时（图15b），噪声比上行波偏移成像噪声严重；(4）全波VSP逆时偏移成像结果（图15c）比仅用上行波或者下行波逆时偏移成像界面形态更清晰．对于此断层模型，全波VSP逆时偏移成像的结果比仅用上行波或下行波成像结果更理想．(III)Marmousi模型层状模型与断层模型试算表明了全波在简单构造中VSP逆时偏移中的有效性，下面针对复杂构造模型进行试算．以Marmousi速度模型（图3）为例，τ=1ms，记录长度为4s，震源为主频为30Hz的Ricker子波，混合吸收边界条件的过渡带点数为10，自适应变空间算子长度方案参数为f(1）上行波、下行波及全波VSP逆时偏移成像位置准确，构造清晰；(2）上行波作为有效波时（图17a），井旁构造成像较清晰，但远井源距处成像结果不清晰，部分构造无法显现；(3）下行波作为有效波时（图17b），远井源距构造成像效果较好，但近井源距成像结果不如上行波成像结果连续；(4）全波作为有效波时（图17c），井旁构造更为清晰、连续，并且远井源距构造也能较好成像．波场分离后的上行波和下行波，都有部分波场信息的丢失，同时在波场分离中难免会损害部分有效信号，而全波波场信息较全，所以采用全波进行VSP逆时偏移效果更好．为验证自适应优化有限差分算法的效率性，我们做了一个计算效率测试，CPU的型号是Intel(R)Xeon(R)CPUE5-26400@2.5GHz，分别采用固定算子长度的泰勒级数展开有限差分方法（fixedlengthTE-basedFDmethod）、自适应泰勒级数展开有限差分方法（adaptiveTE-basedFDmethod）、自适应优化有限差分方法（adaptiveoptimal-basedFDmethod）三种有限差分方法实现Marmousi模型全波VSP逆时偏移．以上三种有限差分方法的M值如图18所示．图19a—c分别为以上三种有限差分方法在时刻t=1200ms时的波场快照．由图19可见，这三种有限差分方法的波场快照模拟精度相近，图19d为固定算子长度的泰勒级数展开有限差分方法M=7（图19a）与自适应优化有限差分方法M=2～4（图19c）波场快照的差值，由图19d可知，两者的误差很小，即两者的精度相近．本文对比了以上三种方法的效率（表1），可得如下结论：(1）自适应泰勒级数展开有限差分方法与固定算子长度的泰勒级数展开有限差分方法对比，效率提高了15%左右；(2）自适应优化有限差分方法与固定算子长度的泰勒级数展开有限差分方法对比，计算量约节省了28%.图19表明了自适应优化有限差分方法的精度，图18和表1表明了该方法的效率性．综上所述，相同精度前提下，在固定算子长度的泰勒级数展开有限差分方法中引入自适应变空间差分算子长度方案可以提高计算效率；再用优化有限差分方法替代泰勒级数展开有限差分方法可以进一步提高计算效率．3.2.3实际vsp资料为进一步验证自适应优化有限差分方法的VSP逆时偏移的有效性与正确性，本文将该有限差分方法应用于某地区的实际VSP资料逆时偏移中．原始地震记录为395炮，炮间距不固定，单个共检波点道集如图20a所示，其炮点没有位于整网格点上．本文采用抗泄露Fourier变换规则化重建方法（4模型实现方案相对于地面地震资料而言，VSP资料能更好地反映井旁信息，所以理论上对井旁的成像处理有着天然的优势．如果能结合VSP资料与地面资料进行井地联合逆时偏移，理论上可以更加清晰地描述地下构造．以图3所示的Marmousi速度模型为例，分别实现地面逆时偏移和VSP逆时偏移，地面勘探的观测系统为地面放炮，炮点x坐标范围为0～4900m，炮间距为100m，共50炮，地面接收，接收点x坐标范围为