SEISMIC TRAVEL TIME TOMOGRAPHIC INVERSION METHOD BASED ON TWO POINT RAY TRACING

Abstract

The present application provides a seismic travel time tomographic inversion method based on two-point ray tracing comprising: collecting seismic data including direct wave travel time and reflected wave travel time; establishing an initial one-dimensional continuously layered model having continuously a varying intraformational velocity; representing a ray parameter p by a variable q, representing a source-receiver distance X by a function X=f(q) of the variable q, solving the function X=f(q) using a Newton iteration method; calculating a theoretical direct wave travel time and reflected wave travel time according to the ray parameter p; comparing the calculated theoretical arrival time with actual arrival time, using an optimal algorithm to adjust velocity parameters of the initial one-dimensional continuously layered model, until a difference between the theoretical direct wave travel time and reflected wave travel time and the actual direct wave travel time and reflected wave travel time complies with a predetermined error standard.

Claims

1. A seismic travel time tomographic inversion method based on two-point ray tracing, comprising: obtaining direct wave travel time data and reflected wave travel time data by collecting seismic data in a research area; performing model parameterizing for the research area, establishing an initial one-dimensional continuously layered model having a continuously varying intraformational velocity; representing a ray parameter p by a variable q according to the one-dimensional continuously layered model having the continuously varying intraformational velocity, representing a source-receiver distance X by a function X=f(q) of the variable q, solving the function X=f(q) by using a Newton iteration method, thereby obtaining the ray parameter p, a ray path is determined solely by the ray parameter p; calculating a theoretical direct wave travel time and reflected wave travel time after the parameter p is obtained; comparing the calculated theoretical direct wave travel time and reflected wave travel time with actual direct wave travel time and reflected wave travel time obtained by collecting the seismic data; determining whether a difference between the theoretical direct wave travel time and reflected wave travel time and the actual direct wave travel time and reflected wave travel time obtained by collecting the seismic data complies with a predetermined error standard; if yes, outputting the one-dimensional continuously layered model having the continuously varying intraformational velocity, if no, performing a next step; adjusting the one-dimensional continuously layered model having the continuously varying intraformational velocity by an optimal algorithm, until the calculated difference between the theoretical direct wave travel time and reflected wave travel time and the actual direct wave travel time and reflected wave travel time obtained by collecting the seismic data complies with the predetermined error standard, and outputting the one-dimensional continuously layered model having the continuously varying intraformational velocity.

2. The seismic travel time tomographic inversion method based on two-point ray tracing according to claim 1, wherein, in the one-dimensional continuously layered model having the continuously varying intraformational velocity, V.sub.k is a function of the depth z, the velocity function of the k-th layer is represented as:
V.sub.k=a.sub.kz+b.sub.k, wherein, a subscript k represents the k-th layer, a.sub.k and b.sub.k are model parameters that need to be inverted, and represent the gradient and the intercept of the velocity function of the k-th layer respectively, when an underground compressional wave velocity model is inverted, the intraformational velocity V.sub.k is a compressional wave velocity; when an underground shear wave velocity model is inverted, the intraformational velocity V.sub.k is the shear wave velocity; when an underground converted wave velocity model is inverted, the intraformational velocity V.sub.k is determined to be a compressional wave velocity or a shear wave velocity by the attribute of the converted wave, which can be a shear-to-compressional converted wave or a compressional-to-shear converted wave.

3. The seismic travel time tomographic inversion method based on two-point ray tracing according to claim 2, wherein, when a.sub.k=0, the k-th layer of the one-dimensional continuously layered model having the continuously varying intraformational velocity is a homogeneous layer having a constant intraformational velocity.

4. The seismic travel time tomographic inversion method based on two-point ray tracing according to claim 1, wherein, the source-receiver distance X is represented as a function of the ray parameter p according to Snell's law, and is formulated as: $X = {.Math.}_{k = 0}^{n} .Math. .Math. [_{k}^{(a)} .Math. \frac{_{k}}{a_{k}} .Math. (\sqrt{p^{2} - {.Math.}_{k}} - \sqrt{p^{- 2} -_{k}}) + (1 -_{k}^{(a)}) .Math. \frac{_{k} .Math. h_{k}}{\sqrt{p^{- 2} - {.Math.}_{k}}}]$ wherein, each of coefficients is defined as follow: ${.Math.}_{k} = {\begin{matrix} {(a_{k} .Math. z^{(k - 1)} + b_{k})}^{2} & k = 1, 2, .Math. .Math., n \\ {(a_{s} .Math. z^{(s - 1)} + b_{s})}^{2} & k = 0 \end{matrix}; .Math._{k} = {\begin{matrix} {(a_{k} .Math. z^{(k)} + b_{k})}^{2} & k = 1, 2, .Math. .Math., n \\ {(a_{s} .Math. z_{s} + b_{s})}^{2} & k = 0 \end{matrix}; .Math. h_{k} = {\begin{matrix} (a_{k} .Math. z^{(k - 1)} + b_{k}) .Math. (z^{(k)} - z^{(k - 1)}) & k = 1, 2, .Math. .Math., n \\ (a_{s} .Math. z^{(s - 1)} + b_{s}) .Math. (z_{s} - z^{(s - 1)}) & k = 0 \end{matrix}; .Math._{k} = {\begin{matrix} 1 & k = 1, .Math. .Math., s - 1 & (all .Math. .Math. waves) \\ 0 & k = s, .Math. .Math., n & (direct .Math. .Math. wave) \\ 2 & k = s, .Math. .Math., n & (reflected .Math. .Math. wave) \\ 1 -_{s} & k = 0 \end{matrix}; .Math._{k}^{(a)} = {\begin{matrix} 1 & a_{k} 0 \\ 0 & a_{k} = 0 \end{matrix}$ wherein, .sub.k, .sub.k, h.sub.k, .sub.k and .sub.k are intermediate parameters, a subscripts represents the layer where the seismic source is located and has a value range of (1, n), n represents a label of the layer where a reflection of the reflected wave occurs, Z.sub.s is the depth of the seismic source, z.sup.(k) represents the depth of the k-th layer, a subscript k represents the k-th layer, the term k=0 is a correction term regarding the location of the seismic source.

5. The seismic travel time tomographic inversion method based on two-point ray tracing according to claim 4, wherein, the ray parameter p is represented by a variable q as $p = \sqrt{\frac{q^{2}}{V_{M}^{2} (1 + q^{2})}},$ wherein, V.sub.M represents the fastest velocity of the ray path passing through the layers.

6. The seismic travel time tomographic inversion method based on two-point ray tracing according to claim 5, wherein, representing the source-receiver distance X as a function of the variable q, the source-receiver distance X is represented as: $X = {.Math.}_{k = 0}^{n} .Math. .Math. [_{k}^{(a)} .Math. {\tilde{}}_{k} (\sqrt{q^{- 2} + {\tilde{.Math.}}_{k}} - \sqrt{q^{- 2} + {\tilde{}}_{k}}) + (1 -_{k}^{(a)}) .Math. \frac{{\tilde{h}}_{k}}{\sqrt{q^{- 2} + {\tilde{.Math.}}_{k}}}]$ wherein, each of coefficients is defined as follow: ${\tilde{}}_{k} = \frac{_{k} .Math. V_{M}}{a_{k}};$ ${\tilde{h}}_{k} = \frac{_{k} .Math. h_{k}}{V_{M}};$ ${\tilde{.Math.}}_{k} = 1 - \frac{{.Math.}_{k}}{V_{M}^{2}};$ ${\tilde{}}_{k} = 1 - \frac{_{k}}{V_{M}^{2}},$ presetting the source-receiver distance X, using a Newton iteration method to solve X=f(q) so as to obtain a value of the parameter q, and substituting the parameter q back to a relational equation $p = \sqrt{\frac{q^{2}}{V_{M}^{2} (1 + q^{2})}}$ of the ray parameter p, thereby obtaining the value of the ray parameter p.

7. The seismic travel time tomographic inversion method based on two-point ray tracing according to claim 6, wherein, when the Newton iteration method is used to solve the equation X=f(q), an initial value of q can be obtained by a method as follows: approximating the source-receiver distance X to the variable q using a rational function formula, the rational function formula is expressed as: $X = \frac{_{1} .Math. q +_{2} .Math. q^{2}}{1 +_{1} .Math. q +_{2} .Math. q^{2}}$ wherein, coefficients .sub.1, .sub.2 .sub.1 and .sub.2 are obtained by equating the coefficients of the Taylor series of the equation of source-receiver distance X and the rational function formula at q.fwdarw.0 and q.fwdarw.. An initial estimate of q is obtained as: $q = \frac{_{1} .Math. X -_{1} + \sqrt{(_{1}^{2} - 4 .Math._{2}) .Math. X^{2} + 2 .Math. (2 .Math._{2} -_{1} .Math._{1}) .Math. X +_{1}^{2}}}{2 .Math. (_{2} -_{2} .Math. X)},$ wherein, expressions of coefficients .sub.1, .sub.2 .sub.1 and .sub.2 are respectively as follows: $_{1} = {.Math.}_{k = 0}^{n} .Math. .Math. [_{k}^{(a)} .Math. \frac{1}{2} .Math. {\tilde{}}_{k} ({\tilde{.Math.}}_{k} - {\tilde{}}_{k}) + (1 -_{k}^{(a)}) .Math. {\tilde{h}}_{k}];$ $_{2} = \frac{c_{0} (c_{0}^{2} + {dc}_{- 1})}{c_{- 1}^{2} - c_{0} .Math. c_{- 2}};$ $_{1} = \frac{c_{0} .Math. c_{- 1} + {dc}_{- 2}}{c_{0} .Math. c_{- 2} - c_{- 1}^{2}};$ $_{2} = \frac{c_{0}^{2} + {dc}_{- 1}}{c_{- 1}^{2} - c_{0} .Math. c_{- 2}};$ $c_{0} = {.Math.}_{k = 0}^{n} .Math. .Math. [_{k}^{(a)} .Math. {\tilde{}}_{k} (_{k}^{(.Math.)} .Math. \sqrt{{\tilde{.Math.}}_{k}} -_{k}^{()} .Math. \sqrt{{\tilde{}}_{k}}) + (1 -_{k}^{(a)}) .Math._{k}^{(.Math.)} .Math. \frac{{\tilde{h}}_{k}}{\sqrt{{\tilde{.Math.}}_{k}}}];$ $c_{1} = {.Math.}_{k = 0}^{n} .Math. .Math. [(1 -_{k}^{(a)}) .Math. (1 -_{k}^{(.Math.)}] .Math. {\tilde{h}}_{k}];$ $c_{- 1} = {.Math.}_{k = 0}^{n} .Math. .Math. [_{k}^{(a)} (_{k}^{()} -_{k}^{(.Math.)}) .Math. {\tilde{}}_{k}];$ $c_{- 2} = {.Math.}_{k = 0}^{n} .Math. .Math. [_{k}^{(a)} .Math. \frac{1}{2} .Math. {\tilde{}}_{k} (\frac{_{k}^{(.Math.)}}{\sqrt{{\tilde{.Math.}}_{k}}} - \frac{_{k}^{()}}{\sqrt{{\tilde{}}_{k}}}) - (1 -_{k}^{(a)}) .Math._{k}^{(.Math.)} .Math. \frac{{\tilde{h}}_{k}}{2 .Math. {\tilde{.Math.}}_{k}^{3 .Math. / .Math. 2}}];$ $_{k}^{(.Math.)} = {\begin{matrix} 1 & {\tilde{.Math.}}_{k} 0 \\ 0 & {\tilde{.Math.}}_{k} = 0 \end{matrix}; .Math._{k}^{()} = {\begin{matrix} 1 & {\tilde{}}_{k} 0 \\ 0 & {\tilde{}}_{k} = 0 \end{matrix};$ using the initial value estimation formula of q to obtain initial value, performing iterative computation, and obtaining an accurate solution of q after the iteration.

8. The seismic travel time tomographic inversion method based on two-point ray tracing according to claim 1, wherein, a calculation formula of the direct wave travel time and the reflected wave travel time is expressed as: $T = {.Math.}_{k = 0}^{n} .Math. .Math. [_{k}^{(a)} .Math. \frac{_{k}}{a_{k}} .Math. \ln (\sqrt{\frac{_{k}}{{.Math.}_{k}}} \frac{1 + \sqrt{1 - {.Math.}_{k} .Math. p^{2}}}{1 + \sqrt{1 -_{k} .Math. p^{2}}}) + (1 -_{k}^{(a)}) .Math. \frac{_{k} .Math. h_{k}}{{.Math.}_{k} .Math. \sqrt{1 - {.Math.}_{k} .Math. p^{2}}}]$

9. The seismic travel time tomographic inversion method based on two-point ray tracing according to claim 1, wherein, the collecting seismic data in the research area is particularly: prospecting seismic signals on the earth's surface, wherein both a seismic source and an array of receivers are arranged at the earth's surface; or vertical seismic profiling, wherein a seismic source is located at the earth's surface, an array of receivers is located in a borehole; or vertical seismic profiling, wherein a seismic source is located in a borehole, a receiver array is located at the earth's surface; or cross-well tomographic imaging, wherein a seismic source and a receiver array are located in different wells.

10. The seismic travel time tomographic inversion method based on two-point ray tracing according to claim 1, wherein a seismic signal used for seismic travel time tomographic inversion can be a compressional wave, or a shear wave, or a compressional-to-shear converted wave, or a shear-to-compressional converted wave.

11. A system for seismic travel time tomographic inversion based on two-point ray tracing, comprising: at least one sensor that collects direct wave travel time seismic data and reflected wave travel time seismic data in a research area; and at least one processor operatively communicating with the sensor, the at least one processor performing the following: parameterize a model for the research area and establishing an initial one-dimensional continuously layered model having a continuously varying intraformational velocity; represent a ray parameter p by a variable q according to the one-dimensional continuously layered model, including representing a source-receiver distance X by a function of the variable q, and solving the function by using a Newton iteration method, thereby obtaining the ray parameter p, determine a ray path solely from the ray parameter p; calculate a theoretical direct wave travel time and reflected wave travel time once the parameter p is obtained; compare the calculated theoretical direct wave travel time and reflected wave travel time with direct wave travel time and reflected wave travel time collected by the sensor; determine whether a difference between the theoretical direct wave travel time and reflected wave travel time and the actual direct wave travel time and reflected wave travel time complies with a predetermined error standard; if the determining determines a difference between the theoretical direct wave travel time and reflected wave travel time and the actual direct wave travel time and reflected wave travel time complies with a predetermined error standard, output the one-dimensional continuously layered model having the continuously varying intraformational velocity, if the determining determines a difference between the theoretical direct wave travel time and reflected wave travel time and the actual direct wave travel time and reflected wave travel time does not comply with a predetermined error standard, adjust the one-dimensional continuously layered model; repeat the comparison, the determining and the adjusting until the calculated difference between the theoretical direct wave travel time and reflected wave travel time and the actual direct wave travel time and reflected wave travel time obtained by collecting the seismic data complies with the predetermined error standard, and output the one-dimensional continuously layered model having the continuously varying intraformational velocity.

12. A system comprising: a seismic data sensor that collects direct wave travel time and reflected wave travel time seismic data; and at least one processor connected to receive the seismic data collected by the sensor, the at least one processor performing the following: establish a one-dimensional continuously layered model having a continuously varying intraformational velocity; represent a ray parameter p by a variable q according to the one-dimensional continuously layered model, including representing a source-receiver distance X by a function of the variable q, and solving the function by using a Newton iteration method, thereby obtaining the ray parameter p, determine a ray path solely from the ray parameter p; calculate a theoretical direct wave travel time and reflected wave travel time based on the ray parameter p and the determined ray path; compare the calculated theoretical direct wave travel time and reflected wave travel time with direct wave travel time and reflected wave travel time seismic data collected by the sensor; and iteratively adjust the one-dimensional continuously layered model based on an iterated comparison until a difference between the theoretical direct wave travel time and reflected wave travel time and the actual direct wave travel time and reflected wave travel time is within a predetermined error tolerance.

Description

BRIEF DESCRIPTION OF THE DRAWINGS

[0035] These and other features and advantages will be better and more completely understood by referring to the following detailed description of non-limiting example embodiments in conjunction with the drawings, of which;

[0036] FIG. 1 illustrates a flow chart of a seismic travel time tomographic inversion method based on two-point ray tracing provided by one example embodiment of the present application;

[0037] FIG. 2 illustrates a parameter definition of a continuously layered model of the seismic travel time tomographic inversion method based on two-point ray tracing provided by the example embodiment of the present application;

[0038] FIG. 3 illustrates a schematic view of a stratum model of the seismic travel time tomographic inversion method based on two-point ray tracing provided by the example embodiment of the present application;

[0039] FIG. 4 illustrates a schematic view of collecting seismic exploration surface data in the seismic travel time tomographic inversion method based on two-point ray tracing provided by the example embodiment of the present application;

[0040] FIG. 5 illustrates a schematic view of earth surface shooting vertical seismic profiling in the seismic travel time tomographic inversion method based on two-point ray tracing provided by the example embodiment of the present application;

[0041] FIG. 6 illustrates a schematic view of downhole shooting vertical seismic profiling in the seismic travel time tomographic inversion method based on two-point ray tracing provided by the example embodiment of the present application;

[0042] FIG. 7 illustrates a schematic view of cross-well imaging in the seismic travel time tomographic inversion method based on two-point ray tracing provided by the example embodiment of the present application.

DETAILED DESCRIPTION OF THE PREFERRED EMBODIMENTS

[0043] In order to make the purpose, the technical solution and the advantages of the present application be clearer and more understandable, the present application will be further described in detail below with reference to accompanying figures and example embodiments. It should be understood that the specific example embodiments described herein are merely intended to illustrate but not to limit the present application.

[0044] Referring to FIGS. 1-3, an example embodiment of the present application provides a seismic travel time tomographic inversion method based on two-point ray tracing comprising following steps:

[0045] step S1, collecting seismic data in a research area. Seismic data is collected in the research area to obtain direct wave travel time data and reflected wave travel time data;

[0046] step S2, model parameterizing. Model parameterizing is performed for the research area, and an initial one-dimensional continuously layered model having a continuously varying intraformational velocity is established;

[0047] particularly, the intraformational velocity can be a compressional wave velocity or a shear wave velocity, and can be increasing or decreasing with depth;

[0048] step S3, computing a ray parameter. A ray parameter p is represented by a variable q according to the one-dimensional continuously layered model having the continuously varying intraformational velocity, a ray parameter p is represented by a variable q, a source-receiver distance Xis represented as a function X=f(q) of the variable q, the function X=f(q) is solved by using a Newton iteration method, thereby obtaining the ray parameter p, a ray path is determined solely by the ray parameter p;

[0049] step S4, obtaining theoretical travel time. Direct wave travel time and reflected wave travel time are calculated according to the ray parameter p;

[0050] step S5, comparing the theoretical travel time with actual travel time. The calculated direct wave travel time and reflected wave travel time are compared with the actual direct wave travel time and reflected wave travel time obtained by collecting the seismic data; it is determined whether a difference between the theoretical direct wave travel time and reflected wave travel time and the actual direct wave travel time and reflected wave travel time obtained by collecting the seismic data complies with a predetermined error standard; if yes, step S7 is performed, otherwise, a step S6 is performed;

[0051] S6, optimizing the model. The one-dimensional continuously layered model having the continuously varying intraformational velocity is adjusted by an optimal algorithm, until the calculated difference between the theoretical direct wave travel time and reflected wave travel time and the actual direct wave travel time and reflected wave travel time complies with the predetermined error standard.

[0052] Particularly, the optimal algorithm used can be a steepest descent method, conjugate gradient inversion, a Newton iteration method, a random search method, etc.

[0053] S7, outputting the model. When the calculated difference between the theoretical travel time and the actual travel time complies with the predetermined error standard, the model is outputted.

[0054] The seismic travel time tomographic inversion method based on two-point ray tracing establishes the one-dimensional continuously layered model having the continuously varying stratum velocity, and represents the intraformational velocity as a function of depth, this model allows intraformational velocities to be continuously varying, reduces the number of required model layers greatly, thereby reducing the amount of inversion variables, such that an actual stratum velocity structure can be described more accurately, the computational efficiency of inversion can be improved significantly. Moreover, the ray parameter p is represented by the variable q, in this way, the ray parameter p can be solved by the variable q indirectly, such that the iteration solving process converges rapidly and stably, and a problem of iteration non-convergence occurring under the condition of a large incident angle can be avoided effectively.

[0055] Furthermore, turning now to FIG. 2, in the one-dimensional continuously layered model having the continuously varying intraformational velocity, the intraformational velocity V.sub.k is a function of depth z, the velocity function of the k-th layer is formulated as:

V.sub.k=a.sub.kz+b.sub.k,

[0056] wherein, a subscript k represents the k-th layer, a.sub.k and b.sub.k are model parameters that need to be inverted, and represent the gradient and the intercept of the velocity function of the k-th layer respectively, when an underground compressional wave velocity model is inverted, the intraformational velocity V.sub.k is the compressional wave velocity; when an underground shear wave velocity model is inverted, the intraformational velocity V.sub.k is the shear wave velocity; when an underground converted wave velocity model is inverted, the intraformational velocity V.sub.k is determined to be the compressional wave or the shear wave velocity by the attribute of the converted wave, which can be a shear-to-compressional converted wave or a compressional-to-shear converted wave.

[0057] When a.sub.k=0, it means that the k-th layer is a homogeneous layer having a constant intraformational velocity.

[0058] Furthermore, the source-receiver distance X is represented as a function of the ray parameter p according to Snell's law, and is formulated as:

[00011] $X = {.Math.}_{k = 0}^{n} .Math. .Math. [_{k}^{(a)} .Math. \frac{_{k}}{a_{k}} .Math. (\sqrt{p^{2} - {.Math.}_{k}} - \sqrt{p^{- 2} -_{k}}) + (1 -_{k}^{(a)}) .Math. \frac{_{k} .Math. h_{k}}{\sqrt{p^{- 2} - {.Math.}_{k}}}]$

[0059] wherein, each of coefficients is defined as follow:

[00012] ${.Math.}_{k} = {\begin{matrix} {(a_{k} .Math. z^{(k - 1)} + b_{k})}^{2} & k = 1, 2, .Math. .Math., n \\ {(a_{s} .Math. z^{(s - 1)} + b_{s})}^{2} & k = 0 \end{matrix}; .Math._{k} = {\begin{matrix} {(a_{k} .Math. z^{(k)} + b_{k})}^{2} & k = 1, 2, .Math. .Math., n \\ {(a_{s} .Math. z_{s} + b_{s})}^{2} & k = 0 \end{matrix}; .Math. h_{k} = {\begin{matrix} (a_{k} .Math. z^{(k - 1)} + b_{k}) .Math. (z^{(k)} - z^{(k - 1)}) & k = 1, 2, .Math. .Math., n \\ (a_{s} .Math. z^{(s - 1)} + b_{s}) .Math. (z_{s} - z^{(s - 1)}) & k = 0 \end{matrix}; .Math._{k} = {\begin{matrix} 1 & k = 1, .Math. .Math., s - 1 & (all .Math. .Math. waves) \\ 0 & k = s, .Math. .Math., n & (direct .Math. .Math. wave) \\ 2 & k = s, .Math. .Math., n & (reflected .Math. .Math. wave) \\ 1 -_{s} & k = 0 \end{matrix}; .Math._{k}^{(a)} = {\begin{matrix} 1 & a_{k} 0 \\ 0 & a_{k} = 0 \end{matrix}$

[0060] wherein, .sub.k, .sub.k, h.sub.k, .sub.k and .sub.k are intermediate parameters, a subscripts represents a layer where the seismic source is located and has a value range of (1, n), n represents a label of the layer where a reflection of the reflected wave occurs, Z.sub.s is the depth of the seismic source, z.sup.(k) represents the depth of the k-th layer, a subscript k represents the k-th layer, the term k=0 is a correction term regarding the location of the seismic source.

[0061] Furthermore, in order to avoid a problem that the iteration non-convergence occurs under the condition of a large incident angle, the ray parameter p is represented by a variable q as

[00013] $p = \sqrt{\frac{q^{2}}{V_{M}^{2} (1 + q^{2})}},$

wherein, V.sub.M represents the fastest velocity of the ray path passing through the layers.

[0062] Furthermore, presenting the source-receiver distance X as a function of the variable q, the source-receiver distance X is represented as:

[00014] $X = {.Math.}_{k = 0}^{n} .Math. .Math. [_{k}^{(a)} .Math. {\tilde{}}_{k} (\sqrt{q^{- 2} + {\tilde{.Math.}}_{k}} - \sqrt{q^{- 2} + {\tilde{}}_{k}}) + (1 -_{k}^{(a)}) .Math. \frac{{\tilde{h}}_{k}}{\sqrt{q^{- 2} + {\tilde{.Math.}}_{k}}}]$

, wherein, each of coefficients is defined as follow:

[00015] ${\tilde{}}_{k} = \frac{_{k} .Math. V_{M}}{a_{k}};$ ${\tilde{h}}_{k} = \frac{_{k} .Math. h_{k}}{V_{M}};$ ${\tilde{.Math.}}_{k} = 1 - \frac{{.Math.}_{k}}{V_{M}^{2}};$ ${\tilde{}}_{k} = 1 - \frac{_{k}}{V_{M}^{2}},$

[0063] presetting the source-receiver distance X, using a Newton iteration method to solve X=f(q) so as to obtain a value of the parameter q, and substituting the parameter q back to a relational equation

[00016] $p = \sqrt{\frac{q^{2}}{V_{M}^{2} (1 + q^{2})}}$

of the ray parameter p, thereby obtaining a value of the ray parameter p.

[0064] Furthermore, when the Newton iteration method is used to solve the equation X=f(q), an initial value of q can be obtained by the following method:

[0065] approximating the source-receiver distance X to the variable q using a rational function formula, the rational function formula is expressed as:

[00017] $X = \frac{_{1} .Math. q +_{2} .Math. q^{2}}{1 +_{1} .Math. q +_{2} .Math. q^{2}}$

[0066] wherein, coefficients .sub.1, .sub.2 .sub.1 and .sub.2 are obtained by equating the coefficients of the Taylor series of the equation of source-receiver distance X and the rational function formula at q.fwdarw.0 and q.fwdarw.. An initial estimate of q is obtained as:

[00018] $q = \frac{_{1} .Math. X -_{1} + \sqrt{(_{1}^{2} - 4 .Math._{2}) .Math. X^{2} + 2 .Math. (2 .Math._{2} -_{1} .Math._{1}) .Math. X +_{1}^{2}}}{2 .Math. (_{2} -_{2} .Math. X)},$

[0067] wherein, expressions of coefficients .sub.1, .sub.2 .sub.1 and .sub.2 are respectively as follows:

[00019] $_{1} = {.Math.}_{k = 0}^{n} .Math. .Math. [_{k}^{(a)} .Math. \frac{1}{2} .Math. {\tilde{}}_{k} ({\tilde{.Math.}}_{k} - {\tilde{}}_{k}) + (1 -_{k}^{(a)}) .Math. {\tilde{h}}_{k}];$ $_{2} = \frac{c_{0} (c_{0}^{2} + {dc}_{- 1})}{c_{- 1}^{2} - c_{0} .Math. c_{- 2}};$ $_{1} = \frac{c_{0} .Math. c_{- 1} + {dc}_{- 2}}{c_{0} .Math. c_{- 2} - c_{- 1}^{2}};$ $_{2} = \frac{c_{0}^{2} + {dc}_{- 1}}{c_{- 1}^{2} - c_{0} .Math. c_{- 2}};$ $c_{0} = {.Math.}_{k = 0}^{n} .Math. .Math. [_{k}^{(a)} .Math. {\tilde{}}_{k} (_{k}^{(.Math.)} .Math. \sqrt{{\tilde{.Math.}}_{k}} -_{k}^{()} .Math. \sqrt{{\tilde{}}_{k}}) + (1 -_{k}^{(a)}) .Math._{k}^{(.Math.)} .Math. \frac{{\tilde{h}}_{k}}{\sqrt{{\tilde{.Math.}}_{k}}}];$ $c_{1} = {.Math.}_{k = 0}^{n} .Math. .Math. [(1 -_{k}^{(a)}) .Math. (1 -_{k}^{(.Math.)}] .Math. {\tilde{h}}_{k}];$ $c_{- 1} = {.Math.}_{k = 0}^{n} .Math. .Math. [_{k}^{(a)} (_{k}^{()} -_{k}^{(.Math.)}) .Math. {\tilde{}}_{k}];$ $c_{- 2} = {.Math.}_{k = 0}^{n} .Math. .Math. [_{k}^{(a)} .Math. \frac{1}{2} .Math. {\tilde{}}_{k} (\frac{_{k}^{(.Math.)}}{\sqrt{{\tilde{.Math.}}_{k}}} - \frac{_{k}^{()}}{\sqrt{{\tilde{}}_{k}}}) - (1 -_{k}^{(a)}) .Math._{k}^{(.Math.)} .Math. \frac{{\tilde{h}}_{k}}{2 .Math. {\tilde{.Math.}}_{k}^{3 .Math. / .Math. 2}}];$ $_{k}^{(.Math.)} = {\begin{matrix} 1 & {\tilde{.Math.}}_{k} 0 \\ 0 & {\tilde{.Math.}}_{k} = 0 \end{matrix}; .Math._{k}^{()} = {\begin{matrix} 1 & {\tilde{}}_{k} 0 \\ 0 & {\tilde{}}_{k} = 0 \end{matrix};$

[0068] using the initial value estimation formula of q to obtain an initial value, performing iterative computation, and obtaining an accurate solution of q after the iteration. The accuracy of the initial value obtained by this method can be very high, for obtaining the accurate solution of the variable q, the iterative computation only needs to be performed for 1 to 3 times.

[0069] Furthermore, after the value of the ray parameter p is obtained, a calculation formula of the direct wave travel time and the reflected wave travel time is expressed as:

[00020] $T = {.Math.}_{k = 0}^{n} .Math. .Math. [_{k}^{(a)} .Math. \frac{_{k}}{a_{k}} .Math. \ln (\sqrt{\frac{_{k}}{{.Math.}_{k}}} \frac{1 + \sqrt{1 - {.Math.}_{k} .Math. p^{2}}}{1 + \sqrt{1 -_{k} .Math. p^{2}}}) + (1 -_{k}^{(a)}) .Math. \frac{_{k} .Math. h_{k}}{{.Math.}_{k} .Math. \sqrt{1 - {.Math.}_{k} .Math. p^{2}}}]$

[0070] Furthermore, referring to FIGS. 4-7, the collecting seismic data in the research area comprises:

[0071] prospecting seismic signals on the earth's surface as shown in FIG. 4, wherein both a seismic source and an array of receivers are arranged at the earth's surface; or

[0072] vertical seismic profiling as shown in FIG. 5, wherein a seismic source is located at the earth's surface, an array of receivers is located in a borehole; or

[0073] vertical seismic profiling as shown in FIG. 6, wherein a seismic source is located in a borehole, a receiver array is located at the earth's surface; or

[0074] cross-well tomographic imaging as shown in FIG. 7, wherein a seismic source and a receiver array are located in different wells.

[0075] Furthermore, the seismic signal used for seismic travel time tomographic inversion can be a compressional wave, or a shear wave, or a compressional-to-shear converted wave, or a shear-to-compressional converted wave.

[0076] The aforementioned example embodiments are only preferred embodiments of the present invention, and should not be regarded as limiting the present invention. Rather, the invention(s) is/are defined by the claims. Any modification, equivalent replacement, improvement, and so on, which are made within the spirit and the principle of the present invention, should be included in the protection scope of the present invention.

SEISMIC TRAVEL TIME TOMOGRAPHIC INVERSION METHOD BASED ON TWO POINT RAY TRACING

Inventors

Cpc classification

Classification Explorer

G01V2210/671

PHYSICS

Classification Explorer

G01V1/303

PHYSICS

Classification Explorer

G01V1/305

PHYSICS

Classification Explorer

G01V2210/6222

PHYSICS

International classification

Classification Explorer

G01V1/30

PHYSICS

Abstract

Claims

Description