基于vi介质的最小二乘ffd叠前深度偏移

基于vi介质的最小二乘ffd叠前深度偏移

1基于vti介质模型的权能:基于分散剂的波场延拓算子随着高精度地震勘探技术的发展，分辨率高、振幅高、保幅效果好的最小二乘偏移成像方法引起了人们的关注。与一次偏移相比,虽然最小二乘偏移在高精度成像方面具有独特优势,但也会引入一系列问题。首先反复迭代会显著增加计算量针对计算量过大目前提出的解决思路大致分为两种:(1)采用预条件或正则化进行约束,加快收敛速度,通过减少迭代次数,达到降低计算量的目的,例如照明预条件对于沉积地层广泛发育的中国中东部地区来说,地下介质存在明显的各向异性特征。忽略介质的各向异性会导致模拟波场和实际传播波场之间存在较大偏差,引入成像噪声,降低成像精度。本文首先通过推导VTI介质高阶FFD波场延拓算子,将常规声波各向同性LSFFD偏移方法拓展到VTI介质;针对最小二乘偏移计算量大、成像噪声明显这两个问题,通过引入平面波编码策略和PWC构造滤波算子进行优化,实现了一种优化的VTI介质平面波最小二乘FFD叠前深度偏移方法;最后,通过修改的各向异性Marmousi2模型进行数值测试,验证方法的有效性和适应性。2vti-ffd高阶波场延拓算子和平面波合成正演算法本文将从三个方面进行论述:首先,介绍时间域和频率域平面波合成的基本原理;其次,将VTI-FFD高阶波场延拓算子和平面波合成原理纳入最小二乘反演成像中,给出叠加和非叠加平面波最小二乘偏移的详细计算流程;最后,在原有的高效算法基础上引入PWC算子,作为构造约束预条件,并在CIG道集上滤除相干成像噪声。2.1差异共炮点道集的成像结果近年来,随着大炮检距、高密度和宽方位角等采集技术的发展,地震数据的数据量显著增加,传统逐炮偏移算法面临计算量过大的挑战。对迭代优化的反演类成像方法,这一问题更加严重。通过将成百上千炮的共炮点道集(CSG)记录压缩成一个平面波道集(CPG)记录,可大量减少地震数据,提高计算效率。图1是平面波合成原理示意图,从中可知平面波合成可分为以下两步完成。一、对不同炮点位置炮记录做线性时移,时移量Δt与平面波射线参数p和震源位置x式中:i表示炮号,i=1,2,…,N二、对时移后的炮记录进行线性叠加式中:x表示检波点坐标;t表示记录时间;d(x,t;p)是合成的平面波道集;d(x,t;s式中:ω为圆频率;i为虚数单位。由式(3)可知,平面波合成(编码)是一种特殊相位编码方式,其相位改变量Δφ=ωΔt利用平面波编码可以将多炮地震记录压缩成一炮平面波记录,通过大量压缩地震数据,实现减小计算量、提高计算效率的目的。由于合成平面波的过程是地震信号相干叠加的过程,因此对平面波道集进行偏移并不会在偏移成像结果中引入串扰噪声。但当使用的平面波道集数目较少时会产生偏移假象,此问题可通过叠加不同入射角度的平面波道集的成像结果进行消除式中:W是检波点的排列长度;f在进行平面波合成时,为了满足地震信号相干叠加这一条件,要求炮间隔尽量小一些,炮间隔的约束条件如下式中:Δs表示炮间隔;p到目前为止,本文推导了地震数据由共炮点道集向平面波道集的转换公式。对于全接收固定观测系统,若采用同一偏移或反演算法,则共炮点道集(N最小二乘偏移通过最小化观测数据和拟合数据之间的残差,对传统一次偏移中存在的振幅不均衡,同相轴分辨率低等问题进行校正以获得高分辨率高保幅的偏移结果。根据Born线性正演理论,频率域一次反射波可表示成式中:f(ω)是频率域震源子波;m(x)是x点的模型参数扰动,即最小二乘偏移成像的目标;U式中:Δz是波场由当前层z式中:v是波场沿VTI介质对称轴方向的传播相速度;ε和δ是Thomsen式(6)可写成更简洁的矢量矩阵相乘的形式式中:L是线性正演算子(反偏移算子),因此线性正演过程也是由模型空间m向数据空间U单程波偏移算子是反偏移算子的伴随式中上标“*”表示伴随,即对复数取共轭,对矩阵取共轭转置。所以偏移结果为式中H被称为汉森矩阵。由式(12)可知,常规偏移成像是由真实模型参数扰动m经过汉森矩阵H滤波得到的结果,由于H的主对角线元素互不相等且不等于单位元素1,非对角线元素不全为0。因此模型空间中某一点的成像结果不仅受到该点对应的真实模型参数扰动的影响,还受到整个模型其他点的影响,传统一次偏移成像结果H基于L式中:Lm表示模拟数据;d式中g表示求取的梯度。其含义为:对模拟数据和观测数据的残差进行偏移成像,将残差由数据空间投影到模型空间,反映数据扰动对模型扰动的影响。当g=0时表示观测数据与模型数据最大程度吻合,偏移结果达到对真实模型参数扰动的最佳逼近,并且随着迭代的进行不再发生变化。利用式(14),给出预条件共轭梯度法的迭代公式为式中:k为迭代次数;β为共轭梯度方向修正参数;z表示共轭梯度方向;α为最优化迭代步长;C为预条件算子。根据目的不同,在具体计算过程中可加入不同的预条件。如为了加速收敛,可加入照明补偿预条件,此时C表示对角汉森矩阵的逆对于式(15),当输入不同道集的观测数据时,对应不同的最小二乘算法,如共炮点道集或者平面波道集。对于平面波道集,根据在迭代计算过程中,是否对不同入射角度平面波道集的偏移结果进行叠加,又可分为叠加平面波最小二乘偏移和非叠加平面波最小二乘平面波偏移由表1可知,叠加和非叠加平面波最小二乘的计算量大致相同,均小于传统最小二乘偏移。与叠加算法相比,非叠加算法对偏移结果的存储量增加N预条件和正则化理论作为反演问题中相对独立的一部分,近年来发展迅速。最小二乘的主要目的是为了拟合得到的观测数据。预条件和正则化则是为了使反演得到的模型符合已知的先验信息(例如得到的更加精确的测井信息,地质信息或者先验假设)而分别对模型引入的直接约束和间接约束。如果先验信息与最小二乘目的一致则通常会加速收敛,如为了使成像振幅更加均衡而引入的照明补偿预条件;当先验信息与最小二乘目的发生冲突时,会使得收敛速度变慢,如为了使模型变化更加平缓而引入的Tikhonov正则化。引入约束条件更重要的目的是改善反演问题的适定性,减少多解性,使得反演结果更加精确合理,而不仅仅是为了改变迭代收敛速度。在一定条件下预条件算子和正则化算子是可互逆转换的。对式(13)引入正则化算子R约束,得目标泛函式中λ是为了均衡正则化项引入比例的一个标量因子。该因子的具体数值通常通过实验测试得到,并且随着迭代次数的增加而逐渐减小。式(16)的最小二乘解为式(16)还可通过引入预条件约束条件进行求解式中P是模型重参数化算子,又称预条件算子,Pp=m。重参数化算子P将反演结果由原始模型空间m变换到另一个更加符合先验信息的模型空间p。该式的解为式中I是单位矩阵。求解式(19)时,标量因子λ则式(19)与式(17)是等价的。本文在求解过程中引入了PWC预条件算子进行约束。原因在于,在偏移速度相对准确的假设条件下,在CIG道集上,任一成像点对应的一次反射波成像同相轴都应该是水平的,即同相轴的局部倾角等于零。而多次波或其他相干噪声对应的同相轴通常具有一定的倾角分布范围,且倾角非零。因此可利用相干噪声和有效构造在CIG道集上的局部倾角差异借助PWC算子进行构造滤波,压制成像噪声,提高信噪比。PWC算子是PWD算子的逆,假设CIG道集共由N道组成,s=[s式中:r是输出的数据残差;D是平面波解构滤波器PWD,其定义如下式中d对预测得到的PWC算子D取逆,可得到平面波构造滤波器C。在最小二乘迭代过程中,PWC算子的应用,对更新梯度在CIG道集上引入了局部平面波平滑约束。通过限定CIG道集上的有效构造倾角趋近于零,可有效地对大倾角相干成像噪声进行滤除。3合成平面波合成为了验证方法的准确性和有效性,本文对修改后Marmousi2模型分别用传统共炮点道集最小二乘偏移算法LSFFD、叠加平面波最小二乘偏移算法PLSFFD方法和非叠加平面波最小二乘偏移算法PPLSFFD进行成像测试,并对不同算法的成像结果进行对比分析。测试所用的Marmousi2模型速度场和各向异性参数场如图2所示。模型尺寸为651×281网格点,纵横向网格间距相同,均为10m。地震数据由时间二阶空间十阶有限差分正演得到,震源子波为主频20Hz的雷克子波。采用全接收的观测系统,共得到326炮地震数据,每炮651个检波点接收,炮间距20m,满足方程5,道间隔10m。采样时间为4s,时间采样间隔为4ms。p(θ)=sinθ/v的取值范围为-0.3~0.3s/km,近地表速度为1500m/s,平面波最大入射角约为28°,由于随着深度增加,速度增大,波场入射角也会随之增大。因此,深层的照明角度范围将大于浅层。对于2000m/s的地层,照明倾角可以达到近37°,对于Marmousi2模型可实现对地下构造的充足照明。共合成了31个平面波记录,虽然利用式(4)计算所需平面波数目更多,但由于最小二乘使用迭代方法反复成像,因此所需平面波数目通常小于理论公式图3a是由Marmousi2模型得到的正演炮记录,由左到右炮点位置分别位于0、1620、3260、4880、6500m处。利用得到的炮记录进行平面波合成,合成结果如图3b所示,由左到右平面波道集对应的射线参数分别为-0.3、-0.16、0、0.16、0.3s/km。由于在平面波记录合成过程中炮点数不是无限大,炮点间隔不是无限小,因此可能会引入端点效应和空间假频。观察平面波记录可以发现,对于不同射线参数,记录中端点效应较弱,无空间假频现象,说明本文炮点间隔选择合理。数值实验测试发现在进行单程波偏移成像时,具有强振幅的直达波会对浅层构造的成像结果产生较大影响。为了避免浅层成像结果受到污染,同时进一步减弱端点效应,本文在平面波合成之前对直达波进行了切除。切除之后合成的平面波如图1所示,经过切除处理后,浅层强直达波噪声被消除,端点效应被完全消除。记录中的低频噪声主要是在直达波和反射波耦合区域直达波切除不完全造成的。由于低频噪声随机分布且能量较弱,对最终偏移成像剖面影响较小。对共炮点道集分别进行ISO-FFD、ISO-LSFFD、VTI-FFD、VTI-LSFFD偏移得到对应成像结果如图4所示。对比各向同性偏移成像结果(图4a、图4b)与各向异性偏移成像结果(图4c、图4d)可发现,在中深层强各向异性区域处,各向同性偏移成像结果中出现偏移画弧,同相轴连续性差、能量弱,偏移噪声明显等问题。这主要是该区域模拟波场与真实波场传播速度不一致造成的。模拟波场的非垂直方向传播速度小于真实波场速度,在地下真实成像点处震源波场与检波波场不满足零延时互相关成像条件,因而反射波无法准确归位,绕射波不能完全收敛,偏移质量变差。更严重的是浅层波场旅行时不准,这一效应会逐层累积,进一步影响深层成像质量,因此即使深层为各向同性介质,成像质量也会变差。经过各向异性矫正之后的成像结果成像精度显著提高,尤其是强各向异性盖层下的构造,同相轴振幅和连续性均得到明显恢复。对比传统一次偏移成像结果(图4a、图4c)和最小二乘10次迭代偏移成像结果(图4b、图4d)可发现,一次偏移结果剖面整体分辨率较低,高陡构造以及弱照明区域同相轴能量太弱以致难以有效识别。而最小二乘偏移可以显著压缩地震子波,拓宽频带宽度,提高成像分辨率,并对弱照明区域进行有效补偿,均衡成像振幅。但是,与一次偏移相比,最小二乘偏移剖面上偏移噪声也更加明显。这是由于在迭代过程中,传统最小二乘偏移算法无法将有效信号和偏移噪声分辨开,一次偏移剖面上原本并不明显的成像噪声被统一看成微弱有效信号进行补偿。随着迭代的进行,一次有效反射波在初始的几次迭代过程中已经基本被完全拟合,而由多次波等引入的相干偏移噪声在迭代过程中无法准确拟合,经过反复补偿之后变得越来越明显。为了消除这一部分噪声必须在成像过程中有效识别偏移噪声并进行迭代压制,这也是本文引入PWC构造滤波算子的初衷。对平面波道集分别进行ISO-PFFD、ISO-PLS-FFD、VTI-PFFD、VTI-PLSFFD处理得到对应偏移成像结果如图5所示。对比图5与图4可见,除了可得到与图4相同的上述结论之外,平面波道集所得成像结果与共炮点道集所得成像结果的精度几乎一致。在VTI-PFFD和VTI-PLSFFD成像剖面上并无空间假频现象出现,这也证明了31个平面波道集对于Marmousi2模型进行成像是足够的。在同等成像精度条件下,平面波偏移计算效率提高了326/31=10.5倍。如果共炮点道集持续增加,计算效率还会进一步提高,这对于工业化大规模应用和三维条件下的最小二乘偏移来说具有重要意义。对比图4与图5偏移结果的浅层可发现,共炮点道集的单程波偏移成像方法浅层同相轴容易受到低频噪声干扰,而平面波道集对应的偏移成像结果对浅层构造具有更高的成像精度。原因在于,对于单程波类的偏移成像方法来说,其波场延拓算子精确的拟合角度难以超过90°,在浅层大炮检距位置即波场近似水平方向传播位置,由于入射角度过大,模拟的震源波场和检波波场是不准确的,如图6a所示,因而应用互相关成像条件之后会呈现出低频/低波数同相轴。而对于平面波道集来说震源是平面波震源,波场传播方向为垂直或近似垂直向下传播,波场入射角度较小(图6b)。震源波场和检波波场能被精确模拟出来,应用互相关成像条件后,成像点清晰准确。由于平面波道集数目较少,在最小二乘偏移迭代过程中针对不同射线参数对应的偏移成像结果可独立更新,而不做叠加处理,在迭代之后叠加。本文称这种成像方法为非叠加平面波最小二乘傅里叶有限差分偏移(PPLSFFD)。这样做有两个好处:(1)扩大了模型空间的维度,使模型具有更大的空间和灵活度,能针对不同射线参数对应的平面波记录进行针对性拟合,从而进一步减小数据残差,提高成像精度;(2)将同一位置不同射线参数对应的成像结果按射线参数由负到正依次排列起来可得到共成像点道集CIG。该道集可用于做成像质控、噪声分析与压制和速度分析等。图7是将PWC构造滤波算子应用于第一次迭代梯度的CIG道集得到的结果。观察图7a(原始道集)可见,在速度和各向异性参数准确的条件下,CIG道集上除了被拉平的一次反射波对应成像同相轴之外,还存在具有较大倾角分布范围的多次波等相干成像噪声;从图7b(滤除的成像噪声)可知应用PWC算子滤除的主要是高陡相干噪声及部分随机高频噪声;图7c(去噪后道集)表明,CIG道集经过构造滤波之后高陡相干噪声已基本被完全压制。图8是平面波道集VTI-PPLSFFD迭代10次偏移成像结果。从图8a(迭代过程中未应用PWC算子)可见盖层以下仍然有部分多次波相干成像噪声残留,但与图6d(应用了PWC算子)相比能量弱了很多。原因在于图6d对应的PLSFFD算法在迭代过程中仅对叠加后的成像剖面进行更新,每次迭代之后剖面上的相干噪声都会增强。而对于PPLSFFD算法来说,不同射线参数对应的成像剖面独立更新,虽然剖面上的相干噪声随着迭代的进行也会增强,但不同剖面上的噪声是不相关的,而有效信号是相关的。因此最终的叠加处理会极大地压制相干成像噪声。然而仅靠叠加并不能彻底将相干噪声去除,应用构造滤波或者其他有效的滤波器将噪声分离出成像剖面才是解决该问题最直接有效的方法。如图8b所示,经过PWC构造滤波该剖面上相干成像噪声已被完全滤除。图9给出了三条典型的归一化数据残差收敛曲线,分别对应平面波道集ISO-PLSFFD算法、VTI-PLSFFD算法、VTI-PPLSFFD算法。观察可以发现,随着迭代的进行,数据残差均逐渐变小,且收敛速度不断变慢,迭代10次之后残差趋于稳定,这也是本文取10次迭代成像结果进行对比的原因。其次,对比前两条收敛曲线可以发现,忽略介质的各向异性会显著增加最终稳定的数据残差。由于共炮点道集ISO-LSFFD和VTI-LSFFD算法对应的残差收敛曲线与这两条曲线基本一致,所以本文不再列出。对比后两条曲线可发现,增加模型的维度,可针对性地拟合不同射线参数对应的平面波道集数据,从而进一步减小数据残差。将图5b、图5d、图8b与三条曲线联系起来可发现,偏移剖面上的成像噪声越弱,对应收敛曲线的稳定残差越小。4新方法的认识借助VTI介质高阶FFD波场延拓算子,本文实现了VTI介质条件下的单程波最小二乘傅里叶有限差分偏移成像方法。针对共炮点道集最小二乘偏移计算量大及相干成像噪声明显这两个问题,通过引入平面波编码技术和PWC构造滤波算子对原有算法进行优化。最后,将三种算法应用于Marmousi2模型进行成像测试,得出以下认识和结论。(