一种大斜井炮检移动VSP自由表面多次波成像方法和装置与流程

本发明涉及地球物理勘探中地震数据偏移成像，特别是涉及一种大斜井炮检移动vsp自由表面多次波成像方法和装置。

背景技术：

炮检移动vsp是针对大斜井设计的观测系统，即炮线布设于地表、位于井中检波器的正上方，炮线随着检波器移动。炮检移动vsp的优势是覆盖次更加均匀。这种观测系统主要用于斜井钻进过程中的导向入靶。通常，自由表面多次波作为干扰波被衰减了。与反射相比，自由表面多次波具有易识别、照明宽、反射角小、传播路径长的特点。

吴世萍等研究了《walkawayvsp多次波成像技术研究》，文献将地震相干成像用于walkawayvsp多次波成像。李建国等研究了《walkawayvsp自由表面多次波叠前深度偏移成像》，文献采用单程波延拓实现了变偏移距vsp的自由表面多次波成像；但是由于针对炮检移动vsp的自由表面多次波成像的研究还不成熟，在地址研究、钻进导向等领域还存在着诸多不便。

技术实现要素：

本发明的目的在于克服现有技术的不足，提供一种大斜井炮检移动vsp自由表面多次波成像方法和装置，有效利用自由表面多次波所携带的丰富信息，实现了炮检移动vsp自由表面多次波叠前深度偏移成像，为地质研究、钻井导向提供了准确的数据基础。

本发明的目的是通过以下技术方案来实现的：一种大斜井炮检移动vsp自由表面多次波成像方法，包括以下步骤：

s1.输入大斜井炮检移动vsp全波场共检波点道集、2维地质模型、偏移参数；

s2.根据输入数据，建立炮点、检波点在2维地质模型中的坐标；

s3.对输入的共检波点道集做快速傅里叶变换到频率波数域，设置频率域震源；

s4.将频率域震源从检波点深度向上延拓至自由表面，频率波数域炮检移动vsp共检波点全波场、震源波场从自由表面向下延拓，实现单检波点的自由表面多次波单程波叠前深度偏移成像；

s5.循环执行步骤s4，完成所有检波点的炮检移动vsp自由表面多次波叠前深度偏移成像；

s6.将步骤s5的成像重排成共成像道集，共成像道集叠加，得到炮检移动vsp自由表面多次波叠前深度偏移叠加成像。

进一步地，所述步骤s1包括：

s101.输入大斜井炮检移动vsp全波场共检波点道集{vspdata}，读取炮点坐标{sx,sy}、炮点深度sz、检波点坐标{rx,ry}、检波点深度rz、井口大地坐标{wmx,wmy}；

s102.输入纵波2维网格层速度地质模型{v}，模型横坐标mx、模型横纵坐标mz、井口模型坐标mwmx信息；

s103.输入成像算子最大倾角dip、成像终止深度zmig2、成像深度步长dzmig、波场{vspdata}最小频率fl、波场{vspdata}最大频率参数fh。

所述步骤s2包括：

s201.利用步骤s1中的炮点坐标、炮点深度、检波点坐标、检波点深度、井口大地坐标、井口模型坐标，建立炮点、检波点在2维地质模型中的坐标：

炮点在2维地质模型中的坐标为：

mszi＝szi

其中，msxi是第i个炮点在2维地质模型中的坐标，sxi、syi是第i个炮点的大地坐标，szi是第i个炮点的相对于2维地质模型的深度，wmx、wmy是井口的大地坐标，mwmx是井口模型坐标，sign是符号函数；

检波点在2维地质模型中的坐标为：

mrzj＝rzj

其中，mrxj是第j个检波点在2维地质模型中坐标，rxj、ryj是第j个检波点的大地坐标，rzj是第j个检波点的深度，wmx、wmy是井口的大地坐标，mwmx是井口模型坐标，sign是符号函数。

所述步骤s3包括：

s301.将步骤s1输入的第j个检波点的时间空间域全波场vspdataj做快速傅里叶变换转换到频率波数域fvspj；

s302.在步骤s2的坐标系下，将频率域震源fsouj设置在第j个检波点处。

所述步骤s4包括：

s401.计算频率域震源向上延拓网格：

kzmig＝(k-1)·dzmig

其中，mrzj是第j个检波点的模型深度，dzmig是成像深度步长，kzmig是第k个向上延拓的深度；

s402.频率域震源用向上延拓一个深度步长：

f＝[fl:df:fh]

kk＝f/v

gs0＝1

gs1＝-i·π·dzmig/kk

gs2＝(-i·π·dzmig·kk/4-π²·dzmig²·kk²/2)/kk⁴

gs3＝(-i·π·dzmig·kk/8+π²·dzmig²·kk²/4-i·π³·dzmig³·kk³/6)/kk⁶

gs4＝(-i·π·dzmig·kk·5/64-π²·dzmig²·kk²·5/32+i·π³·dzmig³·kk³/8+π⁴·dzmig⁴·kk⁴/24)/kk⁸

其中，dip是最大倾角，fl是最小频率，fh是最大频率，f＝[fl:df:fh]是频率向量，表示频率从fl到fh间隔为df，v是速度，kx是波数，dzmig是成像深度步长，i是虚数单位，fft是傅里叶变换函数，ifft是傅里叶逆变换函数，km是最大倾角dip对应的最大波数，mutes是波数切除因子，kk、gs0-gs4、ts0-ts4是中间变量，是向上延拓一个深度步长的波场；

s403.循环执行步骤s401～s402实现频率域震源向上延拓至自由表面；

s404.计算向下延拓成像网格：

kzmig＝(k-1)·dzmig

其中，zmig2是成像终止深度，dzmig是成像深度步长，kzmig是第k个成像深度；

s405.将频率波数域炮检移动vsp共检波点全波场用向下延拓一个深度步长：

f＝[fl:df:fh]

kk＝f/v

g0＝1

g1＝-i·π·dzmig/kk

g2＝(-i·π·dzmig·kk/4-π²·dzmig²·kk²/2)/kk⁴

g3＝(-i·π·dzmig·kk/8+π²·dzmig²·kk²/4-i·π³·dzmig³·kk³/6)/kk⁶

g4＝(-i·π·dzmig·kk·5/64-π²·dzmig²·kk²·5/32+i·π³·dzmig³·kk³/8+π⁴·dzmig⁴·kk⁴/24)/kk⁸

其中，dip是最大倾角，fl是最小频率，fh是最大频率，f＝[fl:df:fh]是频率向量，v是速度，kx是波数，dzmig是成像深度步长，i是虚数单位，fft是傅里叶变换函数，ifft是傅里叶逆变换函数，km是最大倾角dip对应的最大波数，mute是波数切除因子，kk、g0-g4、t0-t4是中间变量，是向下延拓一个深度步长的波场；

s406.将频率域震源向下延拓一个深度步长：

f＝[fl:df:fh]

kk＝f/v

gs0＝1

gs1＝i·π·dzmig/kk

gs2＝(i·π·dzmig·kk/4-π²·dzmig²·kk²/2)/kk⁴

gs3＝(i·π·dzmig·kk/8+π²·dzmig²·kk²/4+i·π³·dzmig³·kk³/6)/kk⁶

gs4＝(i·π·dzmig·kk·5/64-π²·dzmig²·kk²·5/32-i·π³·dzmig³·kk³/8+π⁴·dzmig⁴·kk⁴/24)/kk⁸

其中，dip是最大倾角，fl是最小频率，fh是最大频率，f＝[fl:df:fh]是频率向量，v是速度，kx是波数，dzmig是成像深度步长，i是虚数单位，fft是傅里叶变换函数，ifft是傅里叶逆变换函数，是向下延拓一个深度步长的波场；

s407.相关成像条件提取成像值：

其中，是向下延拓一个深度步长的波场，是向下延拓一个深度步长的波场，是提取的kzmig+dzmig深度的成像值，conj是复共轭函数；

s408.循环执行步骤s404～s407，至zmig2是成像终止深度，完成第j个检波点的炮检移动vsp自由表面多次波叠前深度偏移成像。

一种大斜井炮检移动vsp自由表面多次波成像装置，包括：

数据输入模块，用于输入大斜井炮检移动vsp全波场共检波点道集、2维地质模型、偏移参数；

坐标建立模块，用于根据输入数据，建立炮点、检波点在2维地质模型中的坐标；

多次波叠前成像模块，用于对输入的共检波点道集做快速傅里叶变换到频率波数域，设置频率域震源；将频率域震源从检波点深度向上延拓至自由表面，频率波数域炮检移动vsp共检波点全波场、震源波场从自由表面向下延拓，实现单检波点的自由表面多次波单程波叠前深度偏移成像；并按照同样的方式完成所有检波点的炮检移动vsp自由表面多次波叠前深度偏移成像；

叠加成像模块，用于将得到的多次波叠前深度偏移成像重排成共成像道集，共成像道集叠加，得到炮检移动vsp自由表面多次波叠前深度偏移叠加成像。

本发明的有益效果是：本发明输入炮检移动vsp共检波点全波场、成像速度模型，给定频率域震源波场从检波点处向上延拓至自由表面，再从自由表面向下延拓共检波点全波场、上述的震源波场，采用相关成像条件提取成像值，实现炮检移动vsp自由表面多次波叠前深度偏移成像。

附图说明

图1为本发明的方法流程图；

图2为实施例中大斜井炮检移动vsp全波场共检波点道集示意图；

图3为实施例中大斜井炮检移动vsp的检波点分布示意图；

图4为实施例中大斜井炮检移动vsp的炮点分布示意图；

图5为实施例中成像速度模型示意图；

图6为实施例中炮检移动vsp自由表面多次波叠前深度偏移成像示意图。

图7为实施例中炮检移动vsp自由表面多次波叠前深度偏移叠加成像示意图；

图8为本发明的装置原理框图。

具体实施方式

下面结合附图进一步详细描述本发明的技术方案，但本发明的保护范围不局限于以下所述。

如图1所示，一种大斜井炮检移动vsp自由表面多次波成像方法，包括以下步骤：

s1.输入大斜井炮检移动vsp全波场共检波点道集、2维地质模型、偏移参数：

s101.输入大斜井炮检移动vsp全波场共检波点道集{vspdata}，读取炮点坐标{sx,sy}、炮点深度sz、检波点坐标{rx,ry}、检波点深度rz、井口大地坐标{wmx,wmy}；

s102.输入纵波2维网格层速度地质模型{v}，模型横坐标mx、模型横纵坐标mz、井口模型坐标mwmx信息；

s103.输入成像算子最大倾角dip、成像终止深度zmig2、成像深度步长dzmig、波场{vspdata}最小频率fl、波场{vspdata}最大频率参数fh。

在本申请的实施例中，输入的大斜井炮检移动vsp全波场共检波点道集如图2所示，图2中横坐标为道号，纵坐标为时间(单位：毫秒)；输入的大斜井炮检移动vsp的检波点分布如图3所示，图3中三角形为检波点所在位置，横坐标为长度(单位：米)，纵坐标为深度(单位：米)；输入的大斜井炮检移动vsp的炮点分布如图4所示，图4中，点线为炮线位置；横坐标为长度(单位：米)；纵坐标表示第几次激发；输入的成像速度模型如图5所示，图5中，横坐标为长度(单位：米)；纵坐标为深度(单位：米)。

s2.根据输入数据，建立炮点、检波点在2维地质模型中的坐标；

s201.利用步骤s1中的炮点坐标、炮点深度、检波点坐标、检波点深度、井口大地坐标、井口模型坐标，建立炮点、检波点在2维地质模型中的坐标：

炮点在2维地质模型中的坐标为：

mszi＝szi

其中，msxi是第i个炮点在2维地质模型中的坐标，sxi、syi是第i个炮点的大地坐标，szi是第i个炮点的相对于2维地质模型的深度，wmx、wmy是井口的大地坐标，mwmx是井口模型坐标，sign是符号函数；

检波点在2维地质模型中的坐标为：

mrzj＝rzj

其中，mrxj是第j个检波点在2维地质模型中坐标，rxj、ryj是第j个检波点的大地坐标，rzj是第j个检波点的深度，wmx、wmy是井口的大地坐标，mwmx是井口模型坐标，sign是符号函数。

s3.对输入的共检波点道集做快速傅里叶变换到频率波数域，设置频率域震源：

所述步骤s3包括：

s301.将步骤s1输入的第j个检波点的时间空间域全波场vspdataj做快速傅里叶变换转换到频率波数域fvspj；

s302.在步骤s2的坐标系下，将频率域震源fsouj设置在第j个检波点处。

s4.将频率域震源从检波点深度向上延拓至自由表面，频率波数域炮检移动vsp共检波点全波场、震源波场从自由表面向下延拓，延拓算子为ω-x单程波算子，相关成像条件提取成像值，实现单检波点的自由表面多次波单程波叠前深度偏移成像：

s401.计算频率域震源向上延拓网格：

kzmig＝(k-1)·dzmig

其中，mrzj是第j个检波点的模型深度，dzmig是成像深度步长，kzmig是第k个向上延拓的深度；

s402.频率域震源用向上延拓一个深度步长：

f＝[fl:df:fh]

kk＝f/v

gs0＝1

gs1＝-i·π·dzmig/kk

gs2＝(-i·π·dzmig·kk/4-π²·dzmig²·kk²/2)/kk⁴

gs3＝(-i·π·dzmig·kk/8+π²·dzmig²·kk²/4-i·π³·dzmig³·kk³/6)/kk⁶

gs4＝(-i·π·dzmig·kk·5/64-π²·dzmig²·kk²·5/32+i·π³·dzmig³·kk³/8+π⁴·dzmig⁴·kk⁴/24)/kk⁸

其中，dip是最大倾角，fl是最小频率，fh是最大频率，f＝[fl:df:fh]是频率向量，表示频率从fl到fh间隔为df，v是速度，kx是波数，dzmig是成像深度步长，i是虚数单位，fft是傅里叶变换函数，ifft是傅里叶逆变换函数，km是最大倾角dip对应的最大波数，mutes是波数切除因子，kk、gs0-gs4、ts0-ts4是中间变量，是向上延拓一个深度步长的波场；

s403.循环执行步骤s401～s402实现频率域震源向上延拓至自由表面；

s404.计算向下延拓成像网格：

kzmig＝(k-1)·dzmig

其中，zmig2是成像终止深度，dzmig是成像深度步长，kzmig是第k个成像深度；

s405.将频率波数域炮检移动vsp共检波点全波场用向下延拓一个深度步长：

f＝[fl:df:fh]

kk＝f/v

g0＝1

g1＝-i·π·dzmig/kk

g2＝(-i·π·dzmig·kk/4-π²·dzmig²·kk²/2)/kk⁴

g3＝(-i·π·dzmig·kk/8+π²·dzmig²·kk²/4-i·π³·dzmig³·kk³/6)/kk⁶

g4＝(-i·π·dzmig·kk·5/64-π²·dzmig²·kk²·5/32+i·π³·dzmig³·kk³/8+π⁴·dzmig⁴·kk⁴/24)/kk⁸

其中，dip是最大倾角，fl是最小频率，fh是最大频率，f＝[fl:df:fh]是频率向量，v是速度，kx是波数，dzmig是成像深度步长，i是虚数单位，fft是傅里叶变换函数，ifft是傅里叶逆变换函数，km是最大倾角dip对应的最大波数，mute是波数切除因子，kk、g0-g4、t0-t4是中间变量，是向下延拓一个深度步长的波场；

s406.将频率域震源向下延拓一个深度步长：

f＝[fl:df:fh]

kk＝f/v

gs0＝1

gs1＝i·π·dzmig/kk

gs2＝(i·π·dzmig·kk/4-π²·dzmig²·kk²/2)/kk⁴

gs3＝(i·π·dzmig·kk/8+π²·dzmig²·kk²/4+i·π³·dzmig³·kk³/6)/kk⁶

gs4＝(i·π·dzmig·kk·5/64-π²·dzmig²·kk²·5/32-i·π³·dzmig³·kk³/8+π⁴·dzmig⁴·kk⁴/24)/kk⁸

其中，dip是最大倾角，fl是最小频率，fh是最大频率，f＝[fl:df:fh]是频率向量，v是速度，kx是波数，dzmig是成像深度步长，i是虚数单位，fft是傅里叶变换函数，ifft是傅里叶逆变换函数，是向下延拓一个深度步长的波场；

s407.相关成像条件提取成像值：

其中，是向下延拓一个深度步长的波场，是向下延拓一个深度步长的波场，是提取的kzmig+dzmig深度的成像值，conj是复共轭函数；

s408.循环执行步骤s404～s407，至zmig2是成像终止深度，完成第j个检波点的炮检移动vsp自由表面多次波叠前深度偏移成像。

s5.循环执行步骤s4，完成所有检波点的炮检移动vsp自由表面多次波叠前深度偏移成像；

在本申请的实施例中，炮检移动vsp自由表面多次波叠前深度偏移成像如图6所示，与图2对应，横坐标为道号；纵坐标为深度(单位：米)；

s6.将步骤s5的成像重排成共成像道集，共成像道集叠加，得到炮检移动vsp自由表面多次波叠前深度偏移叠加成像。

在本申请的实施例中，炮检移动vsp自由表面多次波叠前深度偏移叠加成像，如图7所示，横坐标为长度(单位：米)；纵坐标为深度(单位：米)。

如图8所示，一种大斜井炮检移动vsp自由表面多次波成像装置，包括：

数据输入模块，用于输入大斜井炮检移动vsp全波场共检波点道集、2维地质模型、偏移参数；

坐标建立模块，用于根据输入数据，建立炮点、检波点在2维地质模型中的坐标；

多次波叠前成像模块，用于对输入的共检波点道集做快速傅里叶变换到频率波数域，设置频率域震源；将频率域震源从检波点深度向上延拓至自由表面，频率波数域炮检移动vsp共检波点全波场、震源波场从自由表面向下延拓，实现单检波点的自由表面多次波单程波叠前深度偏移成像；并按照同样的方式完成所有检波点的炮检移动vsp自由表面多次波叠前深度偏移成像；

叠加成像模块，用于将得到的多次波叠前深度偏移成像重排成共成像道集，共成像道集叠加，得到炮检移动vsp自由表面多次波叠前深度偏移叠加成像。

综上，本发明输入炮检移动vsp共检波点全波场、成像速度模型，给定频率域震源波场从检波点处向上延拓至自由表面，再从自由表面向下延拓共检波点全波场、上述的震源波场，采用相关成像条件提取成像值，实现炮检移动vsp自由表面多次波叠前深度偏移成像。

以上所述是本发明的优选实施方式，应当理解本发明并非局限于本文所披露的形式，不应该看作是对其他实施例的排除，而可用于其他组合、修改和环境，并能够在本文所述构想范围内，通过上述教导或相关领域的技术或知识进行改动。而本领域人员所进行的改动和变化不脱离本发明的精神和范围，则都应在本发明所附权利要求的保护范围内。