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

Abstract

The invention discloses a method for calculating the seismic wave speed of a stratum by using travel time of through seismic waves of an inclined shaft. According to the method, a position which excites seismic waves is supposed to be shot point S; 2D coordinates of the shot point is (xs, zs); N seismic wave reception points are arranged along a track of a drilling well and marked as R1, R2, R3 to RN successively from top to bottom and from small to large; the reception points are marked in coordinates of (x1, z1), (x2, z2), (x3, z3) to (xN, zN); and practical travel time of the through seismic waves received by the reception points are marked as t1, t2, t3 to tN. The method is suitable for the inclined shaft and can be used to estimate the seismic wave speed of the smooth stratum.

Description

A kind of method utilizing the through seimic travel time of inclined shaft to calculate formation seismic wave velocity

Technical field

The invention belongs to seismic exploration technique field, relate to a kind of method calculating formation seismic wave velocity, especially a kind of method utilizing the through seimic travel time of inclined shaft to calculate formation seismic wave velocity.

Background technology

Field of seismic exploration, in well, seismic method is placed cymoscope at drilling shaft lining and is received the seismic wave of artificial excitation, seismic travel time (also referred to as when walking) according to receiving can estimate the seismic wave propagation speed on stratum or for reflected waveform data Treatment Analysis, and then the physical propertys such as the porosity on stratum, oily, Poisson's ratio can be studied, it is applied to oil-gas exploration and exploitation. Estimate that the method for formation velocity is generally basede on peupendicular hole and assumes currently with the through seimic travel time of inclined shaft.

Summary of the invention

It is an object of the invention to the shortcoming overcoming above-mentioned prior art, it is provided that a kind of method utilizing the through seimic travel time of inclined shaft to calculate formation seismic wave velocity.

It is an object of the invention to be achieved through the following technical solutions:

This method utilizing the through seimic travel time of inclined shaft to calculate formation seismic wave velocity, comprises the following steps:

1) setting the position of earthquake-wave-exciting as shot point S, its two-dimensional coordinate is (x_s,z_s), there is N number of seismic receiving point along drilling well well track, these receive point and are designated as R successively by order from small to large from top to bottom₁,R₂,R₃,...,R_N, the corresponding point coordinates that receives is designated as (x₁,z₁),(x₂,z₂),(x₃,z₃),...,(x_N,z_N), each receive the received through seismic wave of point actual walk time be designated as t successively₁,t₂,t₃,...,t_N;

If there is N number of acline underground, starting number consecutively by order from small to large from 1 from top to bottom, each Bottom surfaces of strata vertical coordinate is corresponding in turn to the vertical coordinate z receiving point from top to bottom₁,z₂,z₃,...,z_N, each formation velocity is followed successively by v from top to bottom₁,v₂,v₃,...,v_N;

2) successively calculating formation velocity by seismic wave from top to bottom with the hypothesis of straightline propagation, computing formula is:

$v_i=(t_i-\underset{j=1}{\overset{i-1}{\Σ}}{L}_j\·v_j)/L_i---\left(1\right)$

Wherein, L_jRepresent from shot point S to receiving some R_jThe seismic propagation path (also referred to as seismic ray) length in jth stratum, its computing formula is:

3) to above calculated interval velocity v_i(i=1,2,3 ..., N) adopt (2M+1) to put sliding window and on average carry out smooth treatment, computing formula is as follows:

$v_i=\left\{\begin{array}{cc}\underset{j=1}{\overset{1+2M}{\&Sigma;}}{L}_j\&CenterDot;v_j/(2M+1),& 1\&le;i\&le;M\\ \underset{j=i-M}{\overset{i+M}{\&Sigma;}}{L}_j\&CenterDot;v_j/(2M+1),& M<i<N-M\\ \underset{j=N-2M}{\overset{N}{\&Sigma;}}{L}_j\&CenterDot;v_j/(2M+1),& N-M\&le;i\&le;N\end{array}\right.---\left(3\right)$

Wherein, the smooth degree of the size reflection result of calculation of M, the value of M is arbitrarily chosen under M < (N-1)/2 condition, determines according to practical situation during calculating; The value making iterations iter is 1;

4) base area seismic wave Snell law calculates in speed v_i(i=1,2,3 ..., N) under the through seimic travel time T of theory_i(i=1,2,3 ..., N), and T when the theory of computation is walked_i(i=1,2,3 ..., N) and actual walk time t_i(i=1,2,3 ..., N) root-mean-square error rms_error;The computational methods of theoretical through seimic travel time adopt ray casting, and tracing process adopts intensive ray shooting method; First, from shot point S, it is being the center of circle along with shot point, is putting a R with shot point S to reception_iIn 90 ° of fan-shaped ranges centered by direction, launch intensive seismic ray with low-angle (such as 0.01 °) interval, the output angle θ after the ray each bed boundary of traverse_iCalculate by following Snell law

${\mathrm{\&theta;}}_i=arcsin\left(\frac{v_i{\mathrm{sin\&theta;}}_{i-1}}{v_{i-1}}\right),i=1,2,3,...,N---\left(4\right)$

Wherein, θ_i-1Represent the ray angle of incidence from the i-th-1 layer entrance i-th layer; Ray is at the intersecting point coordinate (X of each bed boundary_i,Z_i) adopt following two formulas to calculate:

$X_i=\left\{\begin{array}{cc}{x}_s-(z_i-z_S){\mathrm{tan\&theta;}}_{i-1},& (i=1)\\ x_s-\left(z_i-z_{i-1}\right){\mathrm{tan\&theta;}}_{i-1},& \left(i=2,3,4,\mathrm{...},N\right)\end{array}\right.---\left(5\right)$

Z_i=z_i, i=1,2,3 ..., N (6)

Then these rays and vertical line x=x are chosen_iIntersection point to R_iNearest ray is as successful ray, and calculates seimic travel time according to following formula

$T_i=\underset{j=1}{\overset{i}{\&Sigma;}}\frac{L_j}{v_j},i=1,2,3,...,N---\left(7\right)$

Wherein, L_jCalculated by following formula

$L_j=\left\{\begin{array}{cc}\sqrt{{\left(X_1-x_s\right)}^2+{\left(Z_1-z_s\right)}^2},& \left(j=1\right)\\ \underset{k=1}{\overset{j}{\mathrm{\&Sigma;}}}\sqrt{{\left(X_k-X_{k-1}\right)}^2+{\left(Z_k-Z_{k-1}\right)}^2},& \left(j=2,3,4,\mathrm{...},N\right)\end{array}\right.---\left(8\right)$

Can be calculated from above and obtain each reception through seimic travel time of point and corresponding ray path intersecting point coordinate on each bed boundary. T when theory is walked_i(i=1,2,3 ..., N) and actual walk time t_i(i=1,2,3 ..., N) root-mean-square error rms_error following formula calculate

$rms\_error=\sqrt{\frac{\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{\left(T_i-T_i\right)}^2}{N}}---\left(9\right)$

5) an a positive number eps and maximum iteration time Nmax is set, for instance eps=0.01 and Nmax=100. T when if the theory calculated is walked_i(i=1,2,3 ..., N) and actual walk time t_i(i=1,2,3 ..., N) mean square deviation rms_error>eps, and iterations iter<Nmax, then press following formula renewal speed:

${\mathrm{dv}}_i=\frac{(i-1)(T_i-t_i)}{N\underset{j=1}{\overset{i}{\&Sigma;}}{L}_j},i=1,2,3,...,N---\left(10\right)$

Then the value making iterations iter increases 1 and returns step 4); Otherwise, output v_i(i=1,2,3 ..., N) as final calculation result, calculating terminates.

The method have the advantages that

Disclosure is a kind of utilize inclined shaft go directly seimic travel time, suitable in inclined shaft method that smooth formation seismic wave velocity can be estimated.

Accompanying drawing explanation

Fig. 1 is seismic wave observation schematic diagram in well;

Fig. 2 is actual measurement seismic record figure in well;

The through seimic travel time of reality that Fig. 3 (a) is pickup, Fig. 3 (b) is the formation velocity curve calculated;

Fig. 4 utilizes computed formation velocity to process the reflection seismic waves imaging section figure obtained.

Detailed description of the invention

Below in conjunction with accompanying drawing, the present invention is described in further detail:

As shown in Figure 1, it is assumed that the position of earthquake-wave-exciting is shot point S, its two-dimensional coordinate is (x_s,z_s), there is N number of seismic receiving point along drilling well well track, these receive point and are designated as R successively by order from small to large from top to bottom₁,R₂,R₃,...,R_N, the corresponding point coordinates that receives is designated as (x₁,z₁),(x₂,z₂),(x₃,z₃),...,(x_N,z_N), each receive the received through seismic wave of point actual walk time be designated as t successively₁,t₂,t₃,...,t_N。

1. having down N number of acline hypothetically, start number consecutively by order from small to large from 1 from top to bottom, each Bottom surfaces of strata vertical coordinate is corresponding in turn to the vertical coordinate z receiving point from top to bottom₁,z₂,z₃,...,z_N, each formation velocity is followed successively by v from top to bottom₁,v₂,v₃,...,v_N。

2. successively calculating formation velocity by seismic wave from top to bottom with the hypothesis of straightline propagation, computing formula is

$v_i=(t_i-\underset{j=1}{\overset{i-1}{\&Sigma;}}{L}_j\&CenterDot;v_j)/L_i---\left(1\right)$

Wherein, L_jRepresent from shot point S to receiving some R_jThe seismic propagation path (also referred to as seismic ray) length in jth stratum, its computing formula is

3. calculated interval velocity v more than couple_i(i=1,2,3 ..., N) adopt (2M+1) to put sliding window and on average carry out smooth treatment, computing formula is as follows

$v_i=\left\{\begin{array}{cc}\underset{j=1}{\overset{1+2M}{\&Sigma;}}{L}_j\&CenterDot;v_j/(2M+1),& 1\&le;i\&le;M\\ \underset{j=i-M}{\overset{i+M}{\&Sigma;}}{L}_j\&CenterDot;v_j/(2M+1),& M<i<N-M\\ \underset{j=N-2M}{\overset{N}{\&Sigma;}}{L}_j\&CenterDot;v_j/(2M+1),& N-M\&le;i\&le;N\end{array}\right.---\left(3\right)$

Wherein, the smooth degree of the size reflection result of calculation of M, its value arbitrarily can be chosen under M < (N-1)/2 condition, determines according to practical situation during calculating. The value making iterations iter is 1.

4. base area seismic wave Snell law calculates in speed v_i(i=1,2,3 ..., N) under the through seimic travel time T of theory_i(i=1,2,3 ..., N), and T when the theory of computation is walked_i(i=1,2,3 ..., N) and actual walk time t_i(i=1,2,3 ..., N) root-mean-square error rms_error.The computational methods of theoretical through seimic travel time adopt ray casting, and tracing process adopts intensive ray shooting method. First, from shot point S, it is being the center of circle along with shot point, is putting a R with shot point S to reception_iIn 90 ° of fan-shaped ranges centered by direction, launch intensive seismic ray with low-angle (such as 0.01 °) interval, the output angle θ after the ray each bed boundary of traverse_iCalculate by following Snell law

${\mathrm{\&theta;}}_i=arcsin\left(\frac{v_i{\mathrm{sin\&theta;}}_{i-1}}{v_{i-1}}\right),i=1,2,3,...,N---\left(4\right)$

Wherein, θ_i-1Represent the ray angle of incidence from the i-th-1 layer entrance i-th layer, angle of incidence and output angle θ_iDiagram see Fig. 1. Ray is at the intersecting point coordinate (X of each bed boundary_i,Z_i) adopt following two formulas to calculate

$X_i=\left\{\begin{array}{cc}{x}_s-(z_i-z_S){\mathrm{tan\&theta;}}_{i-1},& (i=1)\\ x_s-\left(z_i-z_{i-1}\right){\mathrm{tan\&theta;}}_{i-1},& \left(i=2,3,4,\mathrm{...},N\right)\end{array}\right.---\left(5\right)$

Z_i=z_i, i=1,2,3 ..., N (6)

Then these rays and vertical line x=x are chosen_iIntersection point to R_iNearest ray is as successful ray, and calculates seimic travel time according to following formula

$T_i=\underset{j=1}{\overset{i}{\&Sigma;}}\frac{L_j}{v_j},i=1,2,3,...,N---\left(7\right)$

Wherein, L_jCalculated by following formula

$L_j=\left\{\begin{array}{cc}\sqrt{{\left(X_1-x_s\right)}^2+{\left(Z_1-z_s\right)}^2},& \left(j=1\right)\\ \underset{k=1}{\overset{j}{\mathrm{\&Sigma;}}}\sqrt{{\left(X_k-X_{k-1}\right)}^2+{\left(Z_k-Z_{k-1}\right)}^2},& \left(j=2,3,4,\mathrm{...},N\right)\end{array}\right.---\left(8\right)$

Can be calculated from above and obtain each reception through seimic travel time of point and corresponding ray path intersecting point coordinate on each bed boundary. T when theory is walked_i(i=1,2,3 ..., N) and actual walk time t_i(i=1,2,3 ..., N) root-mean-square error rms_error following formula calculate

$rms\_error=\sqrt{\frac{\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{\left(T_i-T_i\right)}^2}{N}}---\left(9\right)$

5. set an a positive number eps and maximum iteration time Nmax, for instance eps=0.01 and Nmax=100. T when if the theory calculated is walked_i(i=1,2,3 ..., N) and actual walk time t_i(i=1,2,3 ..., N) mean square deviation rms_error>eps, and iterations iter<Nmax, then press following formula renewal speed

${\mathrm{dv}}_i=\frac{(i-1)(T_i-t_i)}{N\underset{j=1}{\overset{i}{\&Sigma;}}{L}_j},i=1,2,3,...,N---\left(10\right)$

Then the value making iterations iter increases 1 and returns the 4th step. Otherwise, output v_i(i=1,2,3 ..., N) as final calculation result, calculating terminates.

If Fig. 2 is the earthquake record received in inclined shaft Hb001. The corresponding shot point of this record is positioned at earth's surface, distance Hb001 well well head 136m, has 206 and receives point, receives point and is equally distributed in tiltedly between deep 100～4200m with 20m along slopes wall. Each through seimic travel time (see Fig. 3 a) of reality receiving point can be picked up according to this earthquake record. By pick up actual walk time by the present invention calculate obtain each stratum, down-hole seimic wave velocity (see Fig. 3 b). This rate curve is relatively smooth, can be directly used for seismic reflection data imaging processing, and gained imaging results is shown in Fig. 4. Bore situation factually, this imaging results success prediction cross directional variations of Carboniferous System cloud rock gas-bearing bed, it is shown that the application effect of the present invention.

Claims (2)

1. one kind utilizes the method that the through seimic travel time of inclined shaft calculates formation seismic wave velocity, it is characterised in that comprise the following steps:

1) setting the position of earthquake-wave-exciting as shot point S, its two-dimensional coordinate is (x_s,z_s), there is N number of seismic receiving point along drilling well well track, these receive point and are designated as R successively by order from small to large from top to bottom₁,R₂,R₃,...,R_N, the corresponding point coordinates that receives is designated as (x₁,z₁),(x₂,z₂),(x₃,z₃),...,(x_N,z_N), each receive the received through seismic wave of point actual walk time be designated as t successively₁,t₂,t₃,...,t_N;

If there is N number of acline underground, starting number consecutively by order from small to large from 1 from top to bottom, each Bottom surfaces of strata vertical coordinate is corresponding in turn to the vertical coordinate z receiving point from top to bottom₁,z₂,z₃,...,z_N, each formation velocity is followed successively by v from top to bottom₁,v₂,v₃,...,v_N;

2) successively calculating formation velocity by seismic wave from top to bottom with the hypothesis of straightline propagation, computing formula is:

$v_i=(t_i-\underset{j=1}{\overset{i-1}{\&Sigma;}}{L}_j\&CenterDot;v_j)/L_i---\left(1\right)$

Wherein, L_jRepresent from shot point S to receiving some R_jSeismic propagation path length in jth stratum, its computing formula is:

3) to above calculated interval velocity v_i, adopt (2M+1) to put sliding window and on average carry out smooth treatment, wherein i=1,2,3 ..., N, computing formula is as follows:

$v_i=\left\{\begin{array}{cc}\underset{j=1}{\overset{1+2M}{\&Sigma;}}{L}_j\&CenterDot;v_j/(2M+1),& 1\&le;i\&le;M\\ \underset{j=1i-M}{\overset{i+M}{\&Sigma;}}{L}_j\&CenterDot;v_j/(2M+1),& M<i<N-M\\ \underset{j=N-2M}{\overset{N}{\&Sigma;}}{L}_j\&CenterDot;v_j/(2M+1),& N-M\&le;i\&le;N\end{array}\right.---\left(3\right)$

Wherein, the smooth degree of the size reflection result of calculation of M, the value of M is arbitrarily chosen under M < (N-1)/2 condition;The value making iterations iter is 1;

4) base area seismic wave Snell law calculates in speed v_iUnder the through seimic travel time T of theory_i, and T when the theory of computation is walked_iWith actual walk time t_iRoot-mean-square error rms_error, wherein i=1,2,3 ... N; Theoretical through seimic travel time T_iComputational methods adopt ray casting;

5) an a positive number eps and maximum iteration time Nmax is set, for instance eps=0.01 and Nmax=100, if T when the theory calculated is walked_i(i=1,2,3 ..., N) and actual walk time t_i(i=1,2,3 ..., N) mean square deviation rms_error>eps, and iterations iter<Nmax, then press following formula renewal speed:

${\mathrm{dv}}_i=\frac{(i-1)(T_i-t_i)}{N\underset{j=1}{\overset{i}{\&Sigma;}}{L}_j},i=1,2,3,\mathrm{...},N---\left(10\right)$

Then the value making iterations iter increases 1 and returns step 4); Otherwise, output v_i(i=1,2,3 ..., N) as final calculation result, calculating terminates.

2. the method utilizing the through seimic travel time of inclined shaft to calculate formation seismic wave velocity according to claim 1, it is characterized in that, step 4) in, tracing process adopts intensive ray shooting method: first from shot point S, is being the center of circle along with shot point, is putting a R with shot point S to reception_iIn 90 ° of fan-shaped ranges centered by direction, launch intensive seismic ray with low-angle interval, the output angle θ after the ray each bed boundary of traverse_iCalculate by following Snell law

${\mathrm{\&theta;}}_i=\arcsin \left(\frac{v_i{\mathrm{sin\&theta;}}_{i-1}}{v_{i-1}}\right),i=1,2,3,...,N---\left(4\right)$

Wherein, θ_i-1Represent the ray angle of incidence from the i-th-1 layer entrance i-th layer; Ray is at the intersecting point coordinate (X of each bed boundary_i,Z_i) adopt following two formulas to calculate:

$X_i=\left\{\begin{array}{cc}{x}_s-(z_i-z_S){\mathrm{tan\&theta;}}_{i-1},& (i=1)\\ x_s-(z_i-z_{i-1}){\mathrm{tan\&theta;}}_{i-1},& (i=2,3,4,...,N)\end{array}\right.---\left(5\right)$

Z_i=z_i, i=1,2,3 ..., N (6)

Then these rays and vertical line x=x are chosen_iIntersection point to R_iNearest ray is as successful ray, and calculates seimic travel time according to following formula

$T_i=\underset{j=1}{\overset{i}{\&Sigma;}}\frac{L_j}{v_j},i=1,2,3,...,N---\left(7\right)$

Wherein, L_jCalculated by following formula

$L_j=\left\{\begin{array}{cc}\sqrt{{(X_1-x_s)}^2+{(Z_1-z_s)}^2},& (j=1)\\ \underset{k=1}{\overset{j}{\&Sigma;}}\sqrt{{(X_k-X_{k-1})}^2+{(Z_k-Z_{k-1})}^2},& (j=2,3,4,...,N)\end{array}\right.---\left(8\right)$

Can be calculated from above and obtain each reception through seimic travel time of point and corresponding ray path intersecting point coordinate on each bed boundary; T when theory is walked_iWith actual walk time t_iRoot-mean-square error rms_error following formula calculate:

$rms\_error=\sqrt{\frac{\underset{i=1}{\overset{N}{\&Sigma;}}{(T_i-t_i)}^2}{N}}---\left(9\right).$