CN105676283A

CN105676283A - Method for calculating seismic wave speed of stratum by using travel time of through seismic waves of inclined shaft

Info

Publication number: CN105676283A
Application number: CN201610040822.1A
Authority: CN
Inventors: 刘池阳; 黄翼坚
Original assignee: NORTHWEST UNIVERSITY
Current assignee: NORTHWEST UNIVERSITY
Priority date: 2015-03-03
Filing date: 2016-01-21
Publication date: 2016-06-15
Anticipated expiration: 2036-01-21
Also published as: CN105676283B

Abstract

The invention discloses a method for calculating the formation seismic wave velocity by using the travel time of the direct seismic wave from an inclined well. The method assumes that the position where the seismic wave is excited is the shot point S, and its two-dimensional coordinates are (x _s , z _s ), and there are N along the drilling trajectory. seismic wave receiving points, these receiving points are recorded as R ₁ , R ₂ , R ₃ ,...,R _N in ascending order from top to bottom, and the corresponding receiving point coordinates are marked as (x ₁ , z ₁ ),(x ₂ ,z ₂ ),(x ₃ ,z ₃ ),...,(x _N ,z _N ), the actual travel time of direct seismic waves received by each receiving point is recorded as t ₁ ,t ₂ , t ₃ ,...,t _N . The invention is suitable for deviated wells and can estimate smooth formation seismic wave velocity.

Description

A Method of Calculating Formation Seismic Wave Velocity Using Direct Seismic Travel Time of Inclined Wells

技术领域technical field

本发明属于地震勘探技术领域，涉及一种计算地层地震波速度的方法，尤其是一种利用斜井直达地震波走时计算地层地震波速度的方法。The invention belongs to the technical field of seismic exploration, and relates to a method for calculating the seismic wave velocity of a formation, in particular to a method for calculating the seismic wave velocity of a formation by using the travel time of the direct seismic wave from an inclined shaft.

背景技术Background technique

地震勘探领域，井中地震方法在钻井井壁放置检波器接收人工激发的地震波，根据所接收的地震波传播时间(也称为走时)可估计地层的地震波传播速度或用于反射波数据处理分析，进而可研究地层的孔隙度、含油气、泊松比等物理性质，应用于油气勘探和开发。目前利用斜井直达地震波走时估计地层速度的方法一般基于垂直井假设。In the field of seismic exploration, in the borehole seismic method, geophones are placed on the drilling wall to receive artificially excited seismic waves. According to the received seismic wave propagation time (also called travel time), the seismic wave propagation velocity of the formation can be estimated or used for reflected wave data processing and analysis, and then It can study the physical properties of the formation such as porosity, oil and gas, Poisson's ratio, and is used in oil and gas exploration and development. At present, the method of estimating formation velocity by using the direct seismic travel time of deviated wells is generally based on the assumption of vertical wells.

发明内容Contents of the invention

本发明的目的在于克服上述现有技术的缺点，提供一种利用斜井直达地震波走时计算地层地震波速度的方法。The purpose of the present invention is to overcome the above-mentioned shortcomings of the prior art, and provide a method for calculating formation seismic wave velocity by using the direct seismic wave travel time of a deviated well.

本发明的目的是通过以下技术方案来实现的：The purpose of the present invention is achieved through the following technical solutions:

这种利用斜井直达地震波走时计算地层地震波速度的方法，包括以下步骤：The method for calculating the formation seismic wave velocity by using the travel time of the direct seismic wave from the inclined well comprises the following steps:

1)设激发地震波的位置为炮点S，其二维坐标为(x_s,z_s)，沿钻井井轨迹有N个地震波接收点，这些接收点从上到下按从小到大的顺序依次记为R₁,R₂,R₃,...,R_N，相应的接收点坐标记为(x₁,z₁),(x₂,z₂),(x₃,z₃),...,(x_N,z_N)，各接收点所接收到直达地震波的实际走时依次记为t₁,t₂,t₃,...,t_N；1) Assume that the position where the seismic wave is excited is the shot point S, and its two-dimensional coordinates are (x _s , z _s ), and there are N seismic wave receiving points along the drilling trajectory, and these receiving points are in ascending order from top to bottom Denoted as R ₁ , R ₂ , R ₃ ,...,R _N , the corresponding receiving point coordinates are marked as (x ₁ ,z ₁ ),(x ₂ ,z ₂ ),(x ₃ ,z ₃ ),. ..,(x _N , z _N ), the actual travel time of direct seismic waves received by each receiving point is recorded as t ₁ , t ₂ , t ₃ ,...,t _N in sequence;

设地下有N个水平地层，从上到下按从小到大的顺序从1开始依次编号，各地层底界面纵坐标从上到下依次对应接收点的纵坐标z₁,z₂,z₃,...,z_N，各地层速度从上到下依次为v₁,v₂,v₃,...,v_N；Assume that there are N horizontal strata underground, numbered from 1 from top to bottom in ascending order, and the vertical coordinates of the bottom interface of each stratum correspond to the vertical coordinates z ₁ , z ₂ , z ₃ of the receiving points from top to bottom. ...,z _N , the velocity of each layer from top to bottom is v ₁ ,v ₂ ,v ₃ ,...,v _N ;

2)按地震波以直线传播的假设从上到下逐层计算地层速度，计算公式为：2) According to the assumption that the seismic wave propagates in a straight line, the formation velocity is calculated layer by layer from top to bottom, and the calculation formula is:

${v v}_{i i} = = (({t t}_{i i} - - {Σ Σ}_{j j = = 11}^{i i - - 11} {L L}_{j j} \cdot \cdot {v v}_{j j})) / / {L L}_{i i} - - - - - - ((11))$

其中，L_j表示从炮点S到接收点R_j的地震波传播路径(也称为地震射线)在第j个地层里的长度，其计算公式为：Among them, L _j represents the length of the seismic wave propagation path (also called seismic ray) from the shot point S to the receiving point R _j in the jth stratum, and its calculation formula is:

3)对以上计算得到的层速度v_i(i＝1,2,3,...,N)采用(2M+1)点滑动窗口平均进行光滑处理，计算公式如下：3) The layer velocity v _i (i=1,2,3,...,N) calculated above is smoothed by (2M+1) point sliding window average, and the calculation formula is as follows:

${v v}_{i i} = = \{\begin{matrix} {Σ Σ}_{j j = = 11}^{11 + + 22 M m} {L L}_{j j} \cdot &Center Dot; {v v}_{j j} / / ((22 M m + + 11)),, & 11 \leq \leq i i \leq \leq M m \\ {Σ Σ}_{j j = = i i - - M m}^{i i + + M m} {L L}_{j j} \cdot \cdot {v v}_{j j} / / ((22 M m + + 11)),, & M m < < i i < < N N - - M m \\ {Σ Σ}_{j j = = N N - - 22 M m}^{N N} {L L}_{j j} \cdot &Center Dot; {v v}_{j j} / / ((22 M m + + 11)),, & N N - - M m \leq \leq i i \leq \leq N N \end{matrix} - - - - - - ((33))$

其中，M的大小反映计算结果的光滑程度，M的值在M<(N-1)/2条件下任意选取，计算时根据实际情况确定；令迭代次数iter的值为1；Among them, the size of M reflects the smoothness of the calculation result, and the value of M is arbitrarily selected under the condition of M<(N-1)/2, and is determined according to the actual situation during calculation; the value of the iteration number iter is 1;

4)根据地震波Snell定律计算在速度v_i(i＝1,2,3,...,N)下的理论直达地震波走时T_i(i＝1,2,3,...,N)，并计算理论走时T_i(i＝1,2,3,...,N)和实际走时t_i(i＝1,2,3,...,N)的均方根误差rms_error；理论直达地震波走时的计算方法采用射线追踪法，追踪过程采用密集射线打靶法；首先从炮点S开始，在沿着以炮点为圆心、以炮点S到接收点R_i方向为中心的90°扇形范围内，以小角度(如0.01°)间隔发射密集的地震波射线，射线穿过各地层界面后的出射角θ_i按以下Snell定律计算4) Calculate the theoretical direct seismic wave travel time T _i (i=1,2,3,...,N) at the velocity v _i (i=1,2,3,...,N) according to the seismic wave Snell's law, And calculate the root mean square error rms_error of the theoretical travel time T _i (i=1,2,3,...,N) and the actual travel time t _i (i=1,2,3,...,N); the theoretical direct The calculation method of seismic wave travel time adopts the ray tracing method, and the tracing process adopts the dense ray shooting method; firstly, starting from the shot point S, along the 90° sector with the shot point as the center and the direction from the shot point S to the receiving point R _i as the center Within the range, dense seismic wave rays are launched at intervals of small angles (such as 0.01°), and the exit angle _θi after the rays pass through the interface of each layer is calculated according to the following Snell law

${θ θ}_{i i} = = a a r r c c s the s i i n no ((\frac{{v v}_{i i} {sinθ sinθ}_{i i - - 11}}{{v v}_{i i - - 11}})),, i i = = 11,, 22,, 33,, ... ...,, N N - - - - - - ((44))$

其中，θ_i-1表示射线从第i-1层进入第i层的入射角；射线在各地层界面的交点坐标(X_i,Z_i)采用以下两式计算：Among them, θ _i-1 represents the incident angle of the ray entering the i-th layer from the i-1th layer; the intersection coordinates (X _i , Z _i ) of the ray at the interface of each layer are calculated by the following two formulas:

${X x}_{i i} = = \{\begin{matrix} {x x}_{s the s} - - (({z z}_{i i} - - {z z}_{S S})) {tanθ tanθ}_{i i - - 11},, & ((i i = = 11)) \\ {x x}_{s the s} - - (({z z}_{i i} - - {z z}_{i i - - 11})) {tanθ tanθ}_{i i - - 11},, & ((i i = = 22,, 33,, 44,, ... ...,, N N)) \end{matrix} - - - - - - ((55))$

Z_i＝z_i,i＝1,2,3,...,N(6)Z _i = z _i , i = 1,2,3,...,N(6)

然后选取这些射线与垂线x＝x_i的交点到R_i最近的射线作为成功射线，并根据下式计算地震波走时Then select the nearest ray from the intersection point of these rays and the vertical line x= _xi to R _i as the successful ray, and calculate the travel time of the seismic wave according to the following formula

${T T}_{i i} = = {Σ Σ}_{j j = = 11}^{i i} \frac{{L L}_{j j}}{{v v}_{j j}},, i i = = 11,, 22,, 33,, ... ...,, N N - - - - - - ((77))$

其中，L_j由下式计算where L _j is calculated by

${L L}_{j j} = = \{\begin{matrix} \sqrt{{(({X x}_{11} - - {x x}_{s the s}))}^{22} + + {(({Z Z}_{11} - - {z z}_{s the s}))}^{22}},, & ((j j = = 11)) \\ {Σ Σ}_{k k = = 11}^{j j} \sqrt{{(({X x}_{k k} - - {X x}_{k k - - 11}))}^{22} + + {(({Z Z}_{k k} - - {Z Z}_{k k - - 11}))}^{22}},, & ((j j = = 22,, 33,, 44,, ... ...,, N N)) \end{matrix} - - - - - - ((88))$

由上可计算得到各接收点直达地震波走时以及相应的射线路径在各地层界面上的交点坐标。理论走时T_i(i＝1,2,3,...,N)和实际走时t_i(i＝1,2,3,...,N)的均方根误差rms_error用下式计算From the above, the direct arrival seismic wave travel time of each receiving point and the intersection coordinates of the corresponding ray paths on the interface of each layer can be obtained. The root mean square error rms_error of theoretical travel time T _i (i=1,2,3,...,N) and actual travel time t _i (i=1,2,3,...,N) is calculated by the following formula

$r r m m s the s__e e r r r r o o r r = = \sqrt{\frac{{Σ Σ}_{i i = = 11}^{N N} {(({T T}_{i i} - - {T T}_{i i}))}^{22}}{N N}} - - - - - - ((99))$

5)设定一个正数eps以及一个最大迭代次数Nmax，例如eps＝0.01及Nmax＝100。若计算的理论走时T_i(i＝1,2,3,...,N)和实际走时t_i(i＝1,2,3,...,N)的均方差rms_error>eps，并且迭代次数iter<Nmax，则按下式更新速度：5) Set a positive number eps and a maximum number of iterations Nmax, for example, eps=0.01 and Nmax=100. If the mean square error rms_error>eps between the calculated theoretical travel time T _i (i=1,2,3,...,N) and the actual travel time t _i (i=1,2,3,...,N), and The number of iterations iter<Nmax, then update the speed according to the formula:

${dv dv}_{i i} = = \frac{((i i - - 11)) (({T T}_{i i} - - {t t}_{i i}))}{N N {Σ Σ}_{j j = = 11}^{i i} {L L}_{j j}},, i i = = 11,, 22,, 33,, ... ...,, N N - - - - - - ((1010))$

然后令迭代次数iter的值增1并返回步骤4)；否则，输出v_i(i＝1,2,3,...,N)作为最终计算结果，计算结束。Then increase the value of the iteration number iter by 1 and return to step 4); otherwise, output v _i (i=1, 2, 3, . . . , N) as the final calculation result, and the calculation ends.

本发明具有以下有益效果：The present invention has the following beneficial effects:

本发明公开一种利用斜井直达地震波走时、适用于斜井、并可估算出光滑的地层地震波速度的方法。The invention discloses a method for directly reaching the travel time of seismic waves by using inclined wells, which is suitable for inclined wells and can estimate the velocity of smooth stratum seismic waves.

附图说明Description of drawings

图1为井中地震波观测示意图；Figure 1 is a schematic diagram of borehole seismic wave observation;

图2为井中实测地震波记录图；Fig. 2 is the record diagram of seismic waves measured in the well;

图3(a)为拾取的实际直达地震波走时，图3(b)为计算的地层速度曲线；Fig. 3(a) is the actual direct seismic travel time picked up, and Fig. 3(b) is the calculated formation velocity curve;

图4为利用所计算地层速度处理得到的反射地震波成像剖面图。Fig. 4 is an imaging profile of reflected seismic waves obtained by processing the calculated formation velocities.

具体实施方式detailed description

下面结合附图对本发明做进一步详细描述：The present invention is described in further detail below in conjunction with accompanying drawing:

如图1所示，假设激发地震波的位置为炮点S，其二维坐标为(x_s,z_s)，沿钻井井轨迹有N个地震波接收点，这些接收点从上到下按从小到大的顺序依次记为R₁,R₂,R₃,...,R_N，相应的接收点坐标记为(x₁,z₁),(x₂,z₂),(x₃,z₃),...,(x_N,z_N)，各接收点所接收到直达地震波的实际走时依次记为t₁,t₂,t₃,...,t_N。As shown in Fig. 1, assuming that the position where the seismic wave is excited is the shot point S, its two-dimensional coordinates are (x _s , z _s ), and there are N seismic wave receiving points along the drilling trajectory, and these receiving points are arranged from small to small from top to bottom The larger sequence is marked as R ₁ , R ₂ , R ₃ ,...,R _N , and the corresponding receiving point coordinates are marked as (x ₁ ,z ₁ ),(x ₂ ,z ₂ ),(x ₃ ,z ₃ ),...,(x _N , z _N ), the actual travel time of the direct seismic waves received by each receiving point is recorded as t ₁ , t ₂ , t ₃ ,...,t _N in turn.

1.假设地下有N个水平地层，从上到下按从小到大的顺序从1开始依次编号，各地层底界面纵坐标从上到下依次对应接收点的纵坐标z₁,z₂,z₃,...,z_N，各地层速度从上到下依次为v₁,v₂,v₃,...,v_N。1. Assuming that there are N horizontal strata underground, they are numbered from small to large from top to bottom, and the vertical coordinates of the bottom interface of each stratum correspond to the vertical coordinates z ₁ , z ₂ , z of the receiving points from top to bottom. ₃ ,...,z _N , the velocities of each layer from top to bottom are v ₁ ,v ₂ ,v ₃ ,...,v _N .

2.按地震波以直线传播的假设从上到下逐层计算地层速度，计算公式为2. According to the assumption that the seismic wave propagates in a straight line, the formation velocity is calculated layer by layer from top to bottom, and the calculation formula is

${v v}_{i i} = = (({t t}_{i i} - - {Σ Σ}_{j j = = 11}^{i i - - 11} {L L}_{j j} \cdot &Center Dot; {v v}_{j j})) / / {L L}_{i i} - - - - - - ((11))$

其中，L_j表示从炮点S到接收点R_j的地震波传播路径(也称为地震射线)在第j个地层里的长度，其计算公式为Among them, L _j represents the length of the seismic wave propagation path (also called seismic ray) from the shot point S to the receiving point R _j in the jth stratum, and its calculation formula is

3.对以上计算得到的层速度v_i(i＝1,2,3,...,N)采用(2M+1)点滑动窗口平均进行光滑处理，计算公式如下3. The layer velocity v _i (i=1,2,3,...,N) calculated above is smoothed by (2M+1) point sliding window average, and the calculation formula is as follows

${v v}_{i i} = = \{\begin{matrix} {Σ Σ}_{j j = = 11}^{11 + + 22 M m} {L L}_{j j} \cdot &Center Dot; {v v}_{j j} / / ((22 M m + + 11)),, & 11 \leq \leq i i \leq \leq M m \\ {Σ Σ}_{j j = = i i - - M m}^{i i + + M m} {L L}_{j j} \cdot &Center Dot; {v v}_{j j} / / ((22 M m + + 11)),, & M m < < i i < < N N - - M m \\ {Σ Σ}_{j j = = N N - - 22 M m}^{N N} {L L}_{j j} \cdot &Center Dot; {v v}_{j j} / / ((22 M m + + 11)),, & N N - - M m \leq \leq i i \leq \leq N N \end{matrix} - - - - - - ((33))$

其中，M的大小反映计算结果的光滑程度，其值可在M<(N-1)/2条件下任意选取，计算时根据实际情况确定。令迭代次数iter的值为1。Among them, the size of M reflects the smoothness of the calculation results, and its value can be selected arbitrarily under the condition of M<(N-1)/2, and it is determined according to the actual situation during calculation. Let the value of iteration count iter be 1.

4.根据地震波Snell定律计算在速度v_i(i＝1,2,3,...,N)下的理论直达地震波走时T_i(i＝1,2,3,...,N)，并计算理论走时T_i(i＝1,2,3,...,N)和实际走时t_i(i＝1,2,3,...,N)的均方根误差rms_error。理论直达地震波走时的计算方法采用射线追踪法，追踪过程采用密集射线打靶法。首先从炮点S开始，在沿着以炮点为圆心、以炮点S到接收点R_i方向为中心的90°扇形范围内，以小角度(如0.01°)间隔发射密集的地震波射线，射线穿过各地层界面后的出射角θ_i按以下Snell定律计算4. Calculate the theoretical direct seismic wave travel time T _i (i=1,2,3,...,N) at the velocity v _i (i=1,2,3,...,N) according to the seismic wave Snell's law, And calculate the root mean square error rms_error of theoretical travel time T _i (i=1,2,3,...,N) and actual travel time t _i (i=1,2,3,...,N). The calculation method of the theoretical direct seismic wave travel time adopts the ray tracing method, and the dense ray shooting method is used in the tracing process. Firstly, starting from the shot point S, within a 90° fan-shaped range with the shot point as the center and the direction from the shot point S to the receiving point R _i as the center, emit dense seismic wave rays at intervals of small angles (such as 0.01°), The exit angle _θi of the ray after passing through the interface of each layer is calculated according to the following Snell's law

其中，θ_i-1表示射线从第i-1层进入第i层的入射角，入射角与出射角θ_i的图示见图1。射线在各地层界面的交点坐标(X_i,Z_i)采用以下两式计算Among them, θi _-1 represents the incident angle of the ray entering the i-th layer from the i-1th layer, and the diagram of the incident angle and the outgoing angle _θi is shown in Figure 1. The intersection coordinates (X _i , Z _i ) of the ray on the interface of each layer are calculated using the following two formulas

Z_i＝z_i,i＝1,2,3,...,N(6)Z _i = z _i , i = 1,2,3,...,N(6)

其中，L_j由下式计算where L _j is calculated by

5.设定一个正数eps以及一个最大迭代次数Nmax，例如eps＝0.01及Nmax＝100。若计算的理论走时T_i(i＝1,2,3,...,N)和实际走时t_i(i＝1,2,3,...,N)的均方差rms_error>eps，并且迭代次数iter<Nmax，则按下式更新速度5. Set a positive number eps and a maximum number of iterations Nmax, for example, eps=0.01 and Nmax=100. If the mean square error rms_error>eps between the calculated theoretical travel time T _i (i=1,2,3,...,N) and the actual travel time t _i (i=1,2,3,...,N), and The number of iterations iter<Nmax, then update the speed according to the formula

然后令迭代次数iter的值增1并返回第4步骤。否则，输出v_i(i＝1,2,3,...,N)作为最终计算结果，计算结束。Then increase the value of the iteration number iter by 1 and return to step 4. Otherwise, output v _i (i=1, 2, 3, . . . , N) as the final calculation result, and the calculation ends.

如图2是在斜井Hb001中接收到的地震记录。该记录相应的炮点位于地表，距离Hb001井井口136m，共有206个接收点，接收点沿斜井井壁以20m等间距分布于斜深100～4200m之间。据该地震记录可拾取到各接收点的实际直达地震波走时(见图3a)。由拾取的实际走时按本发明计算得到井下各地层地震波速度(见图3b)。该速度曲线较为光滑，可直接用于反射地震数据成像处理，所得成像结果见图4。据实钻情况，该成像结果成功预测了石炭系云岩含气层的横向变化，显示了本发明的应用效果。Figure 2 is the seismic record received in the inclined well Hb001. The corresponding shot points of this record are located on the surface, 136m away from the wellhead of Hb001, and there are 206 receiving points in total. The receiving points are distributed at 20m equal intervals along the inclined shaft wall between the inclined depths of 100 and 4200m. According to the seismic records, the actual travel time of direct seismic waves at each receiving point can be picked up (see Fig. 3a). From the actual travel time picked up, the velocity of seismic waves in each formation downhole is calculated according to the present invention (see Fig. 3b). The velocity curve is relatively smooth and can be directly used for reflection seismic data imaging processing. The obtained imaging results are shown in Fig. 4. According to actual drilling conditions, the imaging result successfully predicts the lateral variation of the Carboniferous dolomite gas-bearing layer, showing the application effect of the present invention.

Claims

1. A method for calculating formation seismic wave velocity by direct seismic wave travel time utilizing inclined shaft, is characterized in that, comprises the following steps:

1) Assume that the position where the seismic wave is excited is the shot point S, and its two-dimensional coordinates are (x _s , z _s ), and there are N seismic wave receiving points along the drilling trajectory, and these receiving points are in ascending order from top to bottom Denoted as R ₁ , R ₂ , R ₃ ,...,R _N , the corresponding receiving point coordinates are marked as (x ₁ ,z ₁ ),(x ₂ ,z ₂ ),(x ₃ ,z ₃ ),. ..,(x _N , z _N ), the actual travel time of direct seismic waves received by each receiving point is recorded as t ₁ , t ₂ , t ₃ ,...,t _N in turn;

Assume that there are N horizontal strata underground, numbered from 1 from top to bottom in ascending order, and the vertical coordinates of the bottom interface of each stratum correspond to the vertical coordinates z ₁ , z ₂ , z ₃ of the receiving points from top to bottom. ...,z _N , the velocity of each layer from top to bottom is v ₁ ,v ₂ ,v ₃ ,...,v _N ;

2) According to the assumption that the seismic wave propagates in a straight line, the formation velocity is calculated layer by layer from top to bottom, and the calculation formula is:

{v v}_{i i} = = (({t t}_{i i} - - {Σ Σ}_{j j = = 11}^{i i - - 11} {L L}_{j j} \cdot &Center Dot; {v v}_{j j})) / / {L L}_{i i} - - - - - - ((11))

Among them, L _j represents the length of the seismic wave propagation path from the shot point S to the receiving point R _j in the jth stratum, and its calculation formula is:

3) For the above-calculated layer velocity v _i , use (2M+1) point sliding window average for smoothing, where i=1,2,3,...,N, the calculation formula is as follows:

{v v}_{i i} = = \{\begin{matrix} {Σ Σ}_{j j = = 11}^{11 + + 22 M m} {L L}_{j j} \cdot &Center Dot; {v v}_{j j} / / ((22 M m + + 11)),, & 11 \leq \leq i i \leq \leq M m \\ {Σ Σ}_{j j = = 11 i i - - M m}^{i i + + M m} {L L}_{j j} \cdot &Center Dot; {v v}_{j j} / / ((22 M m + + 11)),, & M m < < i i < < N N - - M m \\ {Σ Σ}_{j j = = N N - - 22 M m}^{N N} {L L}_{j j} \cdot \cdot {v v}_{j j} / / ((22 M m + + 11)),, & N N - - M m \leq \leq i i \leq \leq N N \end{matrix} - - - - - - ((33))

Among them, the size of M reflects the smoothness of the calculation results, and the value of M is arbitrarily selected under the condition of M<(N-1)/2; the value of the iteration number iter is 1;

4) According to the seismic wave Snell's law, calculate the theoretical direct seismic wave travel time T _i at the velocity v _i , and calculate the root mean square error rms_error between the theoretical travel time T _i and the actual travel time t _i , where i=1,2,3,... N; The calculation method of the theoretical direct seismic wave travel time T _i adopts the ray tracing method;

5) Set a positive number eps and a maximum number of iterations Nmax, such as eps=0.01 and Nmax=100, if the calculated theoretical travel time T _i (i=1,2,3,...,N) and the actual travel time t The mean square error of _i (i=1,2,3,...,N) rms_error>eps, and the number of iterations iter<Nmax, update the speed according to the following formula:

{dv dv}_{i i} = = \frac{((i i - - 11)) (({T T}_{i i} - - {t t}_{i i}))}{N N {Σ Σ}_{j j = = 11}^{i i} {L L}_{j j}},, i i = = 11,, 22,, 33,, ... ...,, N N - - - - - - ((1010))

Then increase the value of the iteration number iter by 1 and return to step 4); otherwise, output v _i (i=1, 2, 3, . . . , N) as the final calculation result, and the calculation ends.

2. the method for calculating the formation seismic wave velocity utilizing the direct seismic wave travel time of an inclined shaft according to claim 1, is characterized in that, in step 4), the tracking process adopts the dense ray shooting method: at first starting from the shot point S, along the following The shot point is the center of the circle, and within the 90° fan-shaped range centered on the direction from the shot point S to the receiving point R _i , dense seismic wave rays are emitted at small angular intervals, and the outgoing angle _θi of the rays after passing through the interface of each layer follows the following Snell law calculate

{θ θ}_{i i} = = arcsin arcsin ((\frac{{v v}_{i i} {sinθ sinθ}_{i i - - 11}}{{v v}_{i i - - 11}})),, i i = = 11,, 22,, 33,, ... ...,, N N - - - - - - ((44))

Among them, θ _i-1 represents the incident angle of the ray entering the i-th layer from the i-1th layer; the intersection coordinates (X _i , Z _i ) of the ray at the interface of each layer are calculated by the following two formulas:

{X x}_{i i} = = \{\begin{matrix} {x x}_{s the s} - - (({z z}_{i i} - - {z z}_{S S})) {tanθ tanθ}_{i i - - 11},, & ((i i = = 11)) \\ {x x}_{s the s} - - (({z z}_{i i} - - {z z}_{i i - - 11})) {tanθ tanθ}_{i i - - 11},, & ((i i = = 22,, 33,, 44,, ... ...,, N N)) \end{matrix} - - - - - - ((55))

Z _i = z _i , i = 1,2,3,...,N(6)

Then select the nearest ray from the intersection point of these rays and the vertical line x= _xi to R _i as the successful ray, and calculate the travel time of the seismic wave according to the following formula

{T T}_{i i} = = {Σ Σ}_{j j = = 11}^{i i} \frac{{L L}_{j j}}{{v v}_{j j}},, i i = = 11,, 22,, 33,, ... ...,, N N - - - - - - ((77))

where L _j is calculated by

{L L}_{j j} = = \{\begin{matrix} \sqrt{{(({X x}_{11} - - {x x}_{s the s}))}^{22} + + {(({Z Z}_{11} - - {z z}_{s the s}))}^{22}},, & ((j j = = 11)) \\ {Σ Σ}_{k k = = 11}^{j j} \sqrt{{(({X x}_{k k} - - {X x}_{k k - - 11}))}^{22} + + {(({Z Z}_{k k} - - {Z Z}_{k k - - 11}))}^{22}},, & ((j j = = 22,, 33,, 44,, ... ...,, N N)) \end{matrix} - - - - - - ((88))

From the above, the travel time of the direct seismic wave at each receiving point and the intersection coordinates of the corresponding ray paths on the interface of each layer can be calculated; the root mean square error rms_error of the theoretical travel time T _i and the actual travel time t _i is calculated by the following formula:

r r m m s the s__e e r r r r o o r r = = \sqrt{\frac{{Σ Σ}_{i i = = 11}^{N N} {(({T T}_{i i} - - {t t}_{i i}))}^{22}}{N N}} - - - - - - ((99)) . .