[发明专利]一种起伏地表变偏移距VSP积分叠前深度偏移方法

权利要求书

1.一种起伏地表变偏移距VSP积分叠前深度偏移方法，其特征在于，包括以下步骤：

(1)输入变偏移距VSP共检波点数据、二维地质模型和偏移参数；

步骤(1)具体包括以下子步骤：

(1.1)输入变偏移距VSP共检波点波场数据{VSPData}，读取炮点坐标{SX,SY}、炮点深度SZ、检波点坐标{RX,RY}、检波点深度RZ、井口大地坐标{WMX,WMY}；

(1.2)输入二维网格层速度地质模型{ModVel}，模型横坐标MX、模型纵坐标MZ、井口模型坐标MWMX信息；

(1.3)输入偏移孔径aper、成像起始坐标xmig1、成像结束坐标xmig2、成像网格dxmig、成像终止深度zmig2、成像深度步长dzmig参数；

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

步骤(2)具体包括以下过程：

用步骤(1)中的炮点坐标、炮点深度、检波点坐标、检波点深度、井口大地坐标、井口模型坐标，建立炮点、检波点在二维地质模型中的坐标：

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

MSZ_i＝SZ_i

其中，MSX_i是第i个炮点在二维地质模型中的坐标，SX_i、SY_i是第i个炮点的大地坐标，SZ_i是第i个炮点的相对于二维地质模型的深度、表示起伏地表，WMX、WMY是井口的大地坐标，MWMX是井口模型坐标，sign是符号函数；

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

MRZ_j＝RZ_j

其中，MRX_j是第j个检波点在二维地质模型中坐标，RX_j、RY_j是第j个检波点的大地坐标，RZ_j是第j个检波点的深度，WMX、WMY是井口的大地坐标，MWMX是井口模型坐标，sign是符号函数；

(3)在步骤(2)的模型坐标系中，采用两点射线追踪计算散射点到炮点的走时，保存各个散射点的走时表；

步骤(3)具体包括以下子步骤：

(3.1)用步骤(1)输入的成像起始坐标、成像结束坐标、成像网格、成像终止深度、成像深度步长，计算散射点坐标：

其中，xmig1是成像起始坐标、xmig2是成像结束坐标、dxmig是成像网格，kxtr是散射点序号，xtr是第kxtr列散射点坐标，ztr是第kztr行散射点深度，zmig2是成像终止深度，dzmig成像深度步长；

(3.2)两点射线追踪第kxtr列第kztr行散射点到第i个炮点的走时：

其中，xtr是第kxtr列散射点坐标，ztr是第kztr行散射点深度，MSX_i、MSZ_i是第i个炮点在二维地质模型中的坐标，ModVel_xtr,ztr是xtr、ztr处的速度，max(ModVel_xtr,k·dzmig)是k∈[(ztr-MSZ_i)/dzmig,ktr]中的速度的最大值，dzmig成像深度步长；L是未知数，通过求解非线性方程得到；

其中，p是射线参数，L是上式计算得到，dzmig成像深度步长，xtr是第kxtr列散射点坐标；

其中，p是上式计算的射线参数，dzmig成像深度步长，xtr是第kxtr列散射点坐标，ztr是第kztr行散射点深度，ModVel_xtr,ztr是xtr、ztr处的速度，MSX_i、MSZ_i是第i个炮点在二维地质模型中的坐标，是第i个炮点的xtr、ztr处的走时；

(3.3)循环步骤(3.2)计算第kxtr列的散射点到第i个炮点的走时

(3.4)循环步骤(3.2)～(3.3)，计算第kxtr列的散射点到全部炮点的走时ts_kxtr，保存为第kxtr列散射点的走时表；

(3.5)循环步骤(3.1)～(3.4)，计算全部列散射点到全部炮点的走时表；

(4)在步骤(2)的模型坐标系中，采用两点射线追踪计算散射点到检波点的走时，保存各个检波点的走时表；

步骤(4)具体包括以下子步骤：

(4.1)用步骤(1)输入的成像起始坐标、成像结束坐标、成像网格、成像终止深度、成像深度步长，计算散射点坐标：

其中，xmig1是成像起始坐标、xmig2是成像结束坐标、dxmig是成像网格，kxtr是散射点序号，xtr是第kxtr列散射点坐标，ztr是第kztr行散射点深度，zmig2是成像终止深度，dzmig成像深度步长；

(4.2)两点射线追踪第kxtr列第kztr行散射点到第j个检波点的走时：

其中，xtr是第kxtr列散射点坐标，ztr是第kztr行散射点深度，MRX_j、MRZ_j是第j个炮点在二维地质模型中的坐标，ModVel_xtr,ztr是xtr、ztr处的速度，max(ModVel_xtr,k·dzmig)是中的速度的最大值，dzmig成像深度步长；L是未知数，通过求解非线性方程得到；

其中，p是射线参数，L是上式计算得到，dzmig成像深度步长，xtr是第kxtr列散射点坐标；

其中，p是上式计算的射线参数，dzmig成像深度步长，xtr是第kxtr列散射点坐标，ztr是第kztr行散射点深度，ModVel_xtr,ztr是xtr、ztr处的速度，MRX_j、MRZ_j是第j个检波点在二维地质模型中的坐标，是j个检波点的xtr、ztr处的走时；

(4.3)循环步骤(4.2)计算第kxtr列的散射点到第j个检波点的走时

(4.4)循环步骤(4.1)～(4.3)，计算全部列的散射点到第j个检波点的走时tr^j，保存为第j个检波点到散射点的走时表；

(5)在步骤(2)的模型坐标系中，用步骤(3)、步骤(4)计算共散射走时，实现起伏地表变偏移距VSP共检波点积分叠前深度偏移成像；

步骤(5)具体包括以下子步骤：

(5.1)用步骤(1)的偏移参数定义成像坐标，选择步骤(1)输入数据中第j个检波点的波场数据VSPData_j：

xmig_ktr＝xmig1+(ktr-1)·dxmig ktr＝[1,(xmig2-xmig1)/dxmig]

其中，xmig_ktr是第ktr道的成像坐标，xmig1是成像起始坐标，xmig2是成像结束坐标；

MigData_ktr,j＝{0}

其中，MigData_ktr,j第j个检波点第ktr道的成像数据，先置为零；

(5.2)用步骤(1)的偏移孔径、步骤(2)的炮点模型坐标、步骤(5.1)的成像坐标确定成像输入地震道范围：

xmig_ktr-aper≤MSX_kaper,j≤xmig_ktr+aper

其中，xmig_ktr是第ktr道的成像坐标，aper是偏移孔径，MSX_kaper,j是第j个检波点第ktr个成像道孔径内第kaper个炮点的模型坐标；

(5.3)选择步骤(3)第ktr个散射点的第kaper个炮点的走时加上步骤(4)第个j检波点第ktr个散射点的走时，即为xmig_ktr处kaper、j炮检对的共散射走时；

TSR＝ts_kaper,ktr+tr_j,ktr

其中，TSR是xmig_ktr处kaper、j炮检对的共散射走时，ts_kaper,ktr是第ktr个散射点的第kaper个炮点的走时，tr_j,ktr第个j检波点第ktr个散射点的走时；

(5.4)选择步骤(5.2)中第j个检波点第kaper个炮点的波场数据VSPData_kaper,j，把TSR对应的振幅映射为散射点深度，与xmig_ktr处的成像道MigData_ktr,j积分求和；

(5.5)循环步骤(5.2)～(5.4)，实现第j个检波点第ktr个成像道孔径内全部地震道积分叠前深度偏移成像，得到MigData_ktr,j；

(5.6)循环步骤(5.1)～(5.5)，实现第j个检波点全部成像道的积分叠前深度偏移成像得到MigData_j；

(5.7)并行循环步骤(5.1)～(5.6)，实现全部检波点全部成像道的积分叠前深度偏移成像得到MigData；

(6)步骤(5)的积分叠前深度偏移重排成共成像道集，共成像道集叠加，得到变偏移距VSP积分叠前深度偏移叠加成像。