技术领域
[0001] 本
发明属于叠前深度偏移的技术领域，具体涉及一种基于矩阵乘法的双程波叠前深度偏移的方法。
背景技术
[0002]
地震偏移成像是指在一定的数学模型
基础上(一般以弹性波或
声波方程作为偏移成像的数学模型)建立合适的波场外推算子，用该算子在已知的宏观速度场中对地表观测的数据进
行波场外推(波的反向传播)，并把地表观测到的波场反向外推到产生该反射(或绕射)的反射界面或绕射点上，利用适当的成像条件，计算成像剖面。地震偏移成像的本质是利用数学的手段使地表观测到的地震数据反向传播，消除
地震波的传播效应得到地下结构的图像，其中涉及两个关键问题：波场外推(或称为波场延拓)和成像条件。地震偏移方法大致可以分为三类，即射线积分类偏移方法，单程波方程类偏移方法和双程波类偏移方法。
[0003] 单程波类偏移
算法是以单程波方程为基础展开研究的。单程波方程是将双程声波方程进行上下行波分解得到的，单程波偏移方法就只是对下行波或者是下行波进行深度延拓之后，采用成像条件实现成像计算。常规的单程波方程在上下行波分解过程中，舍弃了振幅矫正项，因此不具备保幅特性(Zhang et al.,2003,2005)。单程波偏移算法作为一种
波动方程类偏移算法，与Kirchhoff积分类偏移方法相比，单程波偏移无高频近似假设，能够处理多值走时和由速度变化所引起的焦散问题，在地震偏移方法发展过程中形成了丰富的发展成果，例如在常规单程波偏移算法中应用广泛的单程波傅里叶有限差分偏移方法和广义屏偏移方法(Gazdag and Sguazzero,1984；Ristow and Ruhl,1994；Wu,1994；Le Rousseau and de Hoop,2001)。但是常规的单程波偏移方法常常都是基于Taylor 级数或者是Pade级数等理论对垂直
波数进行近似处理的，导致了常规单程波偏移算法对强横向速度变化介质和复杂介质的成像能
力有限。
[0004] 要实现对复杂介质的准确成像，特别是实现保幅深度偏移成像，基于准确双程波声波方程的偏移是达到该目的的理论基础，例如逆时偏移算法在复杂构造中的准确成像正说明了该问题。显然，基于单程波方法的偏移方法在成像
精度并不满足该要求。因此，在理论上和实践上都完全有必要开展基于双程波声波方程的波场深度偏移成像的研究。但是，当前开展关于双程波声波方程的波场深度延拓及成像的理论和实践报道还比较少，具有理论研究价值。众所周知，声波方程是一个关于时间变量和空间变量的二阶偏导数微分方程，在理论上，需要两个边界条件方能在深度域上进行准确求解。常规的地震
数据采集系统只是在地表处记录波场值，一般只能提供一个边界条件，即地表初记录的波场值，因此，在理论上缺乏求解双程波声波方程充足的边界条件。Kosloff and Baysal (1983)在地表是均匀介质的假设条件下，利用单程波方法估算出波场值对深度方向的偏导数，结合地表处的波场值实现双程声波方程的波场深度延拓。在Kosloff and Baysal提出的双程声波方程深度偏移算法中使用最大速度建立的低通
滤波器消除在深度延拓过程中的倏逝波，同时也滤掉了一部分陡倾
角构造的反射波
能量，这导致了该算法对陡倾角构造的成像能力不足。Sandberg and Beylkin (2009)在分析了倏逝波和有效波的特征值差异之后，提出了利用谱投影算法消除在双程声波方程深度延拓过程中的倏逝波，实现了利用双程声波方程深度偏移算法对盐丘模型的准确成像。Wu et al.(2012)建立了基于垂直波数及其三角函数表达式的双程波深度延拓算法，并利用单程波传播算子中的Beamlet算子对双程波深度偏移，但是该方法对复杂介质的成像表现出不
稳定性。Pan(2015)建立的双程波深度延拓算法与Wu et al.的方程类似，但是Pan选择利用与常规
相位加插值类似的方法实现双程波深度偏移成像。
发明内容
[0005] 本发明的目的在于针对
现有技术中的上述不足，提供一种基于矩阵乘法的双程波叠前深度偏移的方法，以解决现有技术对双程波深度偏移方案求解的问题。
[0006] 为达到上述目的，本发明采取的技术方案是：
[0007] 一种基于矩阵乘法的双程波叠前深度偏移的方法，其包括：
[0008] S1、将二维常
密度介质下的
频率域声波方程改写成关于深度变量的一阶方程系统：
[0009]
[0010]
[0011] 其中，为在z深度处波场值，H为二阶矩阵， 为在z+Δz深度处波场的偏导数值，x和z分别表示
水平坐标和垂直坐标，垂直坐标即表示深度，v是介质的速度，ω为角频率；
[0012] S2、根据步骤S1中关于深度变量的一阶方程系统，计算得到矩阵H的特征值分解或谱分解：
[0013]
[0014] 其中，I为单位矩阵，λ表示矩阵的特征值 矩阵H 有两个特征值，都是复值；若垂直波数kz为虚数，则矩阵H的特征值为实数，与振幅指数增长的倏逝波有关，引起失稳，故在波场深度外推时应滤掉；
[0015] S3、基于S1中关于深度变量的一阶方程系统计算得到波场深度外推的数学解，并根据矩阵指数函数eHΔz和频域深度变量的初始条件计算得到z到z+Δz 的递归双程波波场深度延拓方程：
[0016]
[0017] 其中， 为在z+Δz深度处波场的偏导数值， 为在 z+Δz深度处波场值，sin(kzΔz)为垂直波数的正弦函数，cos(kzΔz)为垂直波数的余弦函数，为在z深度处的波场值， 为在z深度处的波场偏导数值；
[0018] S4、将含有介质任意速度变化的Helmholtz算子L以矩阵形式离散到某一频率点，并根据矩阵乘法计算垂直波数及垂直波数的三角函数；
[0019] S5、将所述垂直波数及其三角函数带回至步骤S3中，求解得到一般形式下的递归双程波波场深度延拓方程，并利用所得递归双程波波场深度延拓方程对z 深度处的检波点波场 和
震源波场 计算 z+Δz深度处检波点波场和震源波场
[0020] S6、在求得下一个深度的延拓波场之后，利用互相关成像原理计算偏移结果其中*表示共轭运算；
[0021] S7、判断偏移程序是否计算到最大模型深度，如果计算到最大深度处，算法结束；否则重复执行步骤S3-S6；
[0022] S8、S8、基于盐丘模型验证步骤S1-S7在盐丘结构成像的稳定性和准确性。
[0023] 优选地，步骤S3的具体步骤包括：
[0024] 基于S1中关于深度变量的一阶方程系统计算得到深度波场外推的数学解：
[0025]
[0026] 其中，矩阵指数函数eHΔz为：
[0027]
[0028] 其中，I为单位矩阵，λ表示矩阵H的特征值；
[0029] 根据频域的初始条件得到z到z+Δz的递归双程波波场深度延拓方程：
[0030]
[0031] 优选地，步骤S4中求解垂直波数的具体步骤包括：
[0032] 将具有任意速度变化介质的Helmholtz算子L以矩阵形式离散到某一频率点：
[0033]
[0034]
[0035] 其中，ωi是第i个角频率点；
[0036] 将单程波波动方程的垂直波数或平方根算子写成频域-空间域：
[0037]
[0038] 即
[0039] 采用谱投影算法消除矩阵中的负特征值：
[0040] S0＝L/||L||2
[0041]
[0042] 其中，S0为初始矩阵，Sk+1为第k+1次
迭代计算的矩阵，Sk为第k次迭代计算的矩阵；
[0043] 经过有限次迭代计算，Sk→sign(L)，谱投影算子定义为：P＝(I-sign(L))/2，其中sign(L)表示符号矩阵，将谱投影算子作用于原矩阵：
[0044]
[0045] 其中，λk表示矩阵PL的第k个特征值；
[0046] 滤除矩阵的负特征值之后，采用Schulz迭代算法计算矩阵的平方根：
[0047]
[0048]
[0049] 其中，I表示单位矩阵，Y0为初始矩阵，Z0为初始矩阵，Yk+1为第k+1次迭代计算的矩阵，Yk为第k次迭代计算的矩阵，Zk+1为第k+1次迭代计算的矩阵， Zk为第k次迭代计算的矩阵，经过有限次数的迭代计算，Yk→L1/2，Zk→L-1/2，最后Yk的计算结果即为要求解的矩阵L的平方根或者称为垂直波数。
[0050] 优选地，步骤S4中求解垂直波数的三角函数的具体步骤包括：
[0051] 根据计算得到垂直波数，利用泰勒级数展开近似计算垂直波数的三角函数：
[0052]
[0053]
[0054]
[0055] Where x＝kzΔz.。
[0056] 本发明提供的基于矩阵乘法的双程波叠前深度偏移的方法，具有以下有益效果：
[0057] 本发明基于深度变量的一阶方程系统计算得到深度波场外推的数学解，并根据矩阵指数函数和频域深度变量的初始条件计算得到递归双程波波场深度延拓方程；进而采用矩阵乘法计算并构建了关于垂直波数及其三角函数的一般形式下的递归双程波波场深度延拓方程，并进一步利用矩阵乘法理论对双程波深度偏移方程进行准确求解，解决了现有文献中对双程波深度偏移方案求解的问题。
[0058] 除此，本发明通过对盐丘模型的成像，验证了本算法对强速度变化介质准确、稳定的成像能力，体现了本算法在成像上的优越性。
附图说明
[0059] 图1为利用基于Bleamlet传播算子的双程波深度偏移算法计算成像剖面。
[0060] 图2(a)原始速度模型；(b)利用带3×3窗的高斯滤波器生成的平滑速度模型；(c)利用带5×5窗的高斯滤波器生成的平滑速度模型。
[0061] 图3为利用图2b中的速度模型通过不同方法获得的成像剖面，其中，(a) 传统的单程波广义屏偏移算法(GSP)方法；(b)传统的逆时偏移(RTM)方法；(c)本方案提出的双程波深度偏移方法。
[0062] 图4为利用图2c中的速度模型通过不同方法获得的成像剖面，其中，(a) 传统的单程波广义屏偏移算法(GSP)方法；(b)传统的逆时偏移(RTM)方法；(c)本方案提出的双程波深度偏移方法。
具体实施方式
[0063] 下面对本发明的具体实施方式进行描述，以便于
本技术领域的技术人员理解本发明，但应该清楚，本发明不限于具体实施方式的范围，对本技术领域的普通技术人员来讲，只要各种变化在所附的
权利要求限定和确定的本发明的精神和范围内，这些变化是显而易见的，一切利用本发明构思的发明创造均在保护之列。
[0064] 根据本
申请的一个
实施例，本方案的基于矩阵乘法的双程波叠前深度偏移的方法，包括：
[0065] 步骤S1、构建关于深度变量的一阶方程系统；
[0066] 在二维常密度介质中，密度恒定的频域波方程可以表示为：
[0067]
[0068] 其中，x和z分别表示水平坐标和垂直坐标，垂直坐标即为深度，v(x,z)是介质的速度， 是波场压力的频率域表示形式。
[0069] 为了进行深度偏移，将公式(1)写成关于深度变量的一阶方程系统：
[0070]
[0071]
[0072] 其中，为在z深度处波场值，H为二阶矩阵， 为在z+Δz深度处波场的偏导数值，x和z分别表示水平坐标和垂直坐标，垂直坐标即表示深度，v是介质的速度，ω为角频率。
[0073] 与单程波方程不同，可以用一个初始条件来求解，如记录在表面的波场，但方程(2)需要两个边界条件才能求解。因此，除了传统的压力波场外，还要提供波场对深度的偏导数值。根据目前的研究，有两种方法来解决这个问题。一种解决方案是使用单程波波方程估计波场对深度的偏导数值(Kosloff和Baysal, 1983)。该方法适用于海洋勘探，因其具有恒定的
海水速度；另一种方法是使用双
检波器地震采集系统，设置两个不同深度的检波器来记录波场，然后使用有限差分法计算这些波场的偏导数(You et al.，2016)。
[0074] 本实施例中，利用第一种解决方案给出了波场偏导数的初始条件。
[0075] 步骤S2、矩阵H的特征值分解或谱分解；
[0076] 根据步骤S1中关于深度变量的一阶方程系统，计算得到矩阵H的特征值分解或谱分解：
[0077]
[0078] 其中，I为单位矩阵，λ表示矩阵的特征值， 从式 (3)可以看出，矩阵H有两个特征值，都是复值。若垂直波数kz为虚数，则矩阵H的特征值为实数，与振幅指数增长的倏逝波有关，引起不稳定性问题，故在波场深度外推时应滤掉。
[0079] 步骤S3、得到z到z+Δz的递归双程波波场深度延拓方程；
[0080] 基于式(2)深度波场外推的数学解可定义为：
[0081]
[0082] 矩阵指数函数eHΔz为：
[0083]
[0084] 最后，利用频域深度变量的初始条件，实现z到z+Δz的递归双程波波场深度延拓方程。
[0085]
[0086] 在利用式(6)进行深度波场延拓时，传统的地震采集系统往往只提供压力数据，而不考虑压力对深度的偏导数。双
传感器地震采集系统采集两层不同深度的压力数据，可以估算压力对深度的偏导数。对于一般应用，由于在海洋勘探中速度相对稳定，用单程波传播子来估计压力对深度的偏导数。
[0087] 步骤S4、采用矩阵乘法计算垂直波数及垂直波数的三角函数；
[0088] 经过上面介绍的理论推导，需要解决一个核心问题：如何准确地解出(6)式，当速度模型为各向同性或与深度相关的速度模型v＝v(z)时，可采用
相移(PS)方法，在频率波数域中实现式(6)如下式所示：
[0089]
[0090]
[0091] 然而，在实际情况下，很难满足上述假设。实际速度模型在水平和垂直方向上的变化，在某些情况下速度变化非常强烈，比如盐丘模型。直接使用方程(7) 进行波场深度延拓无法满足实际应用需求。
[0092] 解决这一问题的关键仍然在于如何计算垂直波数kz及其三角函数，为了解决这一问题，需要将含有介质任意速度变化的Helmholtz算子L以矩阵形式离散到某一频率点：
[0093]
[0094]
[0095] 其中，ωi是第i个角频率点。
[0096] 在L算子的离散化中，采用二阶有限差分算子进行计算；显而易见，利用高阶有限差分算子可以获得更高的计算精度。在此基础上，提出了计算垂直波数及其三角函数的纯矩阵乘法方法。
[0097] 单程波动方程的垂直波数或平方根算子可以写成频域-空间域写为：
[0098]
[0099] 垂直波数可以看作算子L的平方根，即 而对于一个对称矩阵，如矩阵L，一般对称矩阵包含有负特征值和正特征值，矩阵L中的负特征值对应的是波场传播中的倏逝波，倏逝波在波场深度延拓过程中会呈现指数增长，不利于波场深度延拓算法的稳定性，必须在波场延拓过程中予于消除。在本方案中，采用谱投影算法消除矩阵中的负特征值：
[0100] S0＝L/||L||2
[0101]
[0102] 其中，S0为初始矩阵，Sk+1为第k+1次迭代计算的矩阵，Sk为第k次迭代计算的矩阵；
[0103] 经过有限次迭代计算，Sk→sign(L)，谱投影算子定义为：P＝(I-sign(L))/2，其中，sign(L)表示符号矩阵，将谱投影算子作用于原矩阵：
[0104]
[0105] 其中，λk表示矩阵PL的第k个特征值。
[0106] 在滤除了矩阵的负特征值之后，采用Schulz迭代算法计算矩阵的特质值，[0107]
[0108]
[0109] 其中，I表示单位矩阵，Y0为初始矩阵，Z0为初始矩阵，Yk+1为第k+1次迭代计算的矩阵，Yk为第k次迭代计算的矩阵，Zk+1为第k+1次迭代计算的矩阵， Zk为第k次迭代计算的矩阵，经过有限次数的迭代计算，Yk→L1/2，Zk→L-1/2，最后Yk的计算结果即为要求解的矩阵L的平方根或者称为垂直波数。
[0110] 在得到垂直波数后，利用泰勒级数展开近似计算垂直波数的三角函数[0111]
[0112]
[0113]
[0114] Where x＝kzΔz.
[0115] 步骤S5、将垂直波数及其三角函数带回至步骤S3中，求解得到一般形式下的递归双程波波场深度延拓方程，并利用所得递归双程波波场深度延拓方程对z 深度处的检波点波场 和震源波场 计算 z+Δz深度处检波点波场 和震源波场
[0116] 步骤S6、在求得下一个深度的延拓波场之后，利用互相关成像原理计算偏移结果其中*表示共轭运算。
[0117] 步骤S7、判断偏移结果是否计算到最大模型深度，如果计算到最大深度处，算法结束；否则重复执行步骤S3-S6。
[0118] 步骤S8、验证步骤S1-S7的稳定性和准确性。
[0119] 盐丘模型是复杂模型成像评估的经典模型。由于高速盐体的存在，对盐丘体下方构造的成像成为偏移算法研究的重点。
[0120] 在之前的研究中，Wu等人(2012)采用基于Bleamlet传播算子的双程波动方程的波场深度外推方法来成像盐模型。然而，Wu等(2012)提出的双程波bleamlet 传播器未能获得满意的结果，成像结果如图1所示。在他们的工作中，他们只提供了一个叠后成像剖面，在盐丘下有很强的成像噪音。
[0121] 在本方案中，使用矩阵乘法的方法来进行双程波的波场深度延拓，保证了稳定性和准确性。为了研究在速度反演不准确的情况下偏移算法的成像性能，本发明采用了两个平滑的速度模型。一种平滑的速度模型是通过使用带有3×3 个窗口的高斯平滑来产生的，另一种是通过使用带有5×5个窗口的高斯平滑来产生的，原始的和平滑的速度模型如图2所示。
[0122] 在数值实验中，分别采用常规的单程波GSP方法、常规的RTM方法和本方案提出的偏移方法，利用图2b和图2c中的速度模型对盐丘模型进行成像，成像结果如图3和图4所示。
[0123] 对比图1、图3c和图4c的成像剖面，发现本方案提出的利用矩阵乘法的双程波偏移方法成功地解决了之前算法的不稳定性问题，特别是对盐丘下方构造的成像。传统的单程波GSP方法很难对盐丘下的结构提供准确的成像结果，特别是在大角度传播过程中或具有较强横向速度变化的介质中。传统的单程波偏移方法容易带来成像噪声，不利于地震解释，如图3a、3a所示。
[0124] 对于RTM方法，高速盐体也给其带来了不可回避的问题，如盐体内部的散射问题，导致了盐丘体下方的多次被虚假成像问题。仔细对比图3b和图3c的偏移剖面，本方案提出的双波偏移方法与传统的RTM方法相比，似乎能提供更清晰的结构
接触关系，成像噪声更小，即使在图2b中，盐丘下的反射轴几乎不可见。然而，当速度模型变得更加模糊和不准确时，如图2c所示，利用矩阵乘法的双程波偏移方法所提供的偏移剖面变得越来越不准确，这说明在速度模型不太精确的情况下，很难描述细尺度结构。
[0125] 虽然结合附图对发明的具体实施方式进行了详细地描述，但不应理解为对本
专利的保护范围的限定。在权利要求书所描述的范围内，本领域技术人员不经创造性劳动即可做出的各种
修改和
变形仍属本专利的保护范围。
