阅读说明：本技术 一种地震波旅行时间计算方法 (Seismic wave travel time calculation method ) 是由 魏建 孙祥娥 王文松 梁晏翔 于 2020-05-21 设计创作，主要内容包括：本发明涉及地震勘探数据处理技术领域,尤其涉及一种地震波旅行时间计算方法；包括以下步骤：S1,确定基本旅行时间计算公式基本形式；S2,根据所述泰勒级数原理和连续分式近似旅行时间平方公式对所述基本旅行时间计算公式进行更高阶的近似；S3,对更高阶的所述基本旅行时间计算公式中偏移距的4次方和6次方进行同时降幂处理,输出旅行时间；本发明实施例通过对基本旅行时间计算公式进行更高阶的近似,然后对更高阶的所述基本旅行时间计算公式降幂处理输出旅行时间；操作简单,运算高效。(The invention relates to the technical field of seismic exploration data processing, in particular to a seismic wave travel time calculation method; the method comprises the following steps: s1, determining the basic form of the basic travel time calculation formula; s2, according to the Taylor series principle and continuous fraction approximate travel time square formula, carrying out higher-order approximation on the basic travel time calculation formula; s3, simultaneously performing power reduction processing on the 4 th power and the 6 th power of the offset in the higher-order basic travel time calculation formula, and outputting travel time; the embodiment of the invention carries out higher-order approximation on a basic travel time calculation formula, and then lowers the power of the higher-order basic travel time calculation formula to output the travel time; the operation is simple and the operation is efficient.)

1. A seismic travel time calculation method is characterized by comprising the following steps:

s1, determining the basic form of the basic travel time calculation formula; the basic travel time calculation formula is an expression of the square of the basic travel time obtained by expanding according to the Taylor series principle and utilizing an anisotropic medium model; the basic travel time calculation formula is as follows:

T²(x)＝c₀+c₁x²+c₂x⁴+c₃x⁶+c₄x⁸+...

c₀＝t²

wherein x is an offset distance; c. C₀、c₁、c₂、c₃、c₄Coefficients that are Taylor series; t ═ T (0) is the time of spontaneous excitation and spontaneous reception, v ═ v_NMOIndicates NMO speed, S_k(k ═ 2,3,4) is a heterogeneity parameter;

s2, according to the Taylor series principle and the continuous fraction approximate travel time square formula, carrying out higher-order approximation on the basic travel time calculation formula to obtain the following formula:

wherein, the parameter B₁Is a kind of coefficient approximate expression based on 6 th order Taylor expansion and 4 th order Taylor expansion, parameter B₂Is based on a coefficient approximation of the 8 th order Taylor expansion and the 6 th order Taylor expansion;

s3, in x²Simultaneously performing power reduction processing on the 4 th power and the 6 th power of the offset distance in the higher-order basic travel time calculation formula for the cardinal number, and outputting travel time; the basic travel time calculation formula after the power reduction processing is as follows:

wherein, the common approximation can be takenTo simplify the parameter c in the formula₃、B₁、B₂The simplified parameters are as follows:

B₂＝0

wherein S₂1+8 η, anisotropy parameterShowing the strength of anisotropy, showing the change in NMO velocity and reflection amplitude, the sum of the parameters being the Thomsen anisotropy parameter.

2. The method of claim 1, wherein the step S3 is further followed by forming an exponentiation compensation formula for the exponentiation-processed basic travel time calculation formula in combination with an approximate compensation method:

wherein fb is the first arrival time, T_DE(x) And the travel time calculated by the calculation formula of the basic travel time after the power reduction is calculated, wherein n is the number of model layers.

3. A seismic travel time computing system, comprising:

the calculation optimization module determines the basic form of the basic travel time calculation formula; the basic travel time calculation formula is an expression of the square of the basic travel time obtained by expanding according to the Taylor series principle and utilizing an anisotropic medium model; the basic travel time calculation formula is as follows:

T²(x)＝c₀+c₁x²+c₂x⁴+c₃x⁶+c₄x⁸+...

c₀＝t²

wherein x is an offset distance; c. C₀、c₁、c₂、c₃、c₄Coefficients that are Taylor series; t ═ T (0) is the time of spontaneous excitation and spontaneous reception, v ═ v_NMOIndicates NMO speedDegree, S_k(k ═ 2,3,4) is a heterogeneity parameter;

and according to the Taylor series principle and the continuous fraction approximate travel time square formula, carrying out higher-order approximation on the basic travel time calculation formula to obtain the following formula:

wherein, the parameter B₁Is a kind of coefficient approximate expression based on 6 th order Taylor expansion and 4 th order Taylor expansion, parameter B₂Is based on a coefficient approximation of the 8 th order Taylor expansion and the 6 th order Taylor expansion;

calculation output module, with x²Simultaneously performing power reduction processing on the 4 th power and the 6 th power of the offset distance in the higher-order basic travel time calculation formula for the cardinal number, and outputting travel time; the basic travel time calculation formula after the power reduction processing is as follows:

wherein, the common approximation can be takenTo simplify the parameter c in the formula₃、B₁、B₂The simplified parameters are as follows:

wherein S₂1+8 η, anisotropy parameter

4. The seismic travel time calculation system of claim 3, wherein the calculation output module further comprises a power reduction compensation formula formed by combining the base travel time calculation formula after power reduction processing with an approximate compensation method:

wherein fb is the first arrival time, T_DE(x) And the travel time calculated by the calculation formula of the basic travel time after the power reduction is calculated, wherein n is the number of model layers.

Technical Field

The invention relates to the technical field of seismic exploration data processing, in particular to a seismic wave travel time calculation method.

Background

The seismic wave travel time plays an important role in the seismic data processing process; and the limitation of the number of rays participating in calculation, reflectors in different shapes and the like is added while the innovation of part of formulas is sought, so that the subsequent application of the formula to actual seismic data is influenced, and therefore, how to acquire more accurate and effective seismic wave travel time is one of key contents needing to be researched currently.

The currently applied methods for calculating the travel time are a shift hyperbolic method and a generalized time difference approximation method, which can effectively calculate the travel time.

Disclosure of Invention

In order to overcome the defects of the prior art, the embodiment of the invention provides a method for calculating the travel time of the seismic wave, which can reduce the difference between the calculated travel time and the first arrival time to a certain extent, and has the advantages of simple operation and high calculation efficiency.

In one aspect, an embodiment of the present invention provides a seismic travel time calculation method, including the following steps:

s1, determining the basic form of the basic travel time calculation formula; the basic travel time calculation formula is an expression of the square of the basic travel time obtained by expanding according to the Taylor series principle and utilizing an anisotropic medium model; the basic travel time calculation formula is as follows:

T²(x)＝c₀+c₁x²+c₂x⁴+c₃x⁶+c₄x⁸+...

c₀＝t²

wherein x is an offset distance; c. C₀、c₁、c₂、c₃、c₄Coefficients that are Taylor series; t ═ T (0) is the time of spontaneous excitation and spontaneous reception, v ═ v_NMOIndicates NMO speed, S_k(k ═ 2,3,4) is a heterogeneity parameter;

s2, according to the Taylor series principle and the continuous fraction approximate travel time square formula, carrying out higher-order approximation on the basic travel time calculation formula to obtain the following formula:

wherein, the parameter B₁Is a kind of coefficient approximate expression based on 6 th order Taylor expansion and 4 th order Taylor expansion, parameter B₂Is based on a coefficient approximation of the 8 th order Taylor expansion and the 6 th order Taylor expansion;

s3, in x²Simultaneously performing power reduction processing on the 4 th power and the 6 th power of the offset distance in the higher-order basic travel time calculation formula for the cardinal number, and outputting travel time; the basic travel time calculation formula after the power reduction processing is as follows:

wherein, the common approximation can be takenTo simplify the parameter c in the formula₃、B₁、B₂The simplified parameters are as follows:

wherein S₂1+8 η, anisotropy parameter

Showing the strength of anisotropy, showing the change in NMO velocity and reflection amplitude, the sum of the parameters being the Thomsen anisotropy parameter.

In another aspect, an embodiment of the present invention provides a seismic wave travel time calculation system, including:

the calculation optimization module determines the basic form of the basic travel time calculation formula; the basic travel time calculation formula is an expression of the square of the basic travel time obtained by expanding according to the Taylor series principle and utilizing an anisotropic medium model; the basic travel time calculation formula is as follows:

T²(x)＝c₀+c₁x²+c₂x⁴+c₃x⁶+c₄x⁸+...

c₀＝t²

wherein x is an offset distance; c. C₀、c₁、c₂、c₃、c₄Coefficients that are Taylor series; t ═ T (0) is the time of spontaneous excitation and spontaneous reception, v ═ v_NMOIndicates NMO speed, S_k(k ═ 2,3,4) is a heterogeneity parameter;

and according to the Taylor series principle and the continuous fraction approximate travel time square formula, carrying out higher-order approximation on the basic travel time calculation formula to obtain the following formula:

wherein, the parameter B₁Is a kind of coefficient approximate expression based on 6 th order Taylor expansion and 4 th order Taylor expansion, parameter B₂Is based on a coefficient approximation of the 8 th order Taylor expansion and the 6 th order Taylor expansion;

calculation output module, with x²Simultaneously performing power reduction processing on the 4 th power and the 6 th power of the offset distance in the higher-order basic travel time calculation formula for the cardinal number, and outputting travel time; the basic travel time calculation formula after the power reduction processing is as follows:

wherein, the common approximation can be taken

To simplify the parameter c in the formula₃、B₁、B₂The simplified parameters are as follows:

wherein S₂1+8 η, anisotropy parameter

Showing the strength of anisotropy, showing the change in NMO velocity and reflection amplitude, the sum of the parameters being Thompsen anisotropy parameter.

The embodiment of the invention provides a seismic wave travel time calculation method and a system, wherein higher-order approximation is carried out on a basic travel time calculation formula, and then the travel time is output by lowering power of the higher-order basic travel time calculation formula; the operation is simple and the operation is efficient.

Drawings

In order to more clearly illustrate the technical solution of the present invention, the drawings needed to be used in the technical description of the present invention will be briefly introduced below, and it is apparent that the drawings in the following description are only some embodiments of the present invention, and it is obvious for those skilled in the art that other drawings can be obtained according to the drawings without inventive labor.

FIG. 1 is a schematic flow chart of a seismic travel time calculation method according to an embodiment of the present invention;

FIG. 2 is a schematic view of a planar horizontal layered medium model;

FIG. 3 is a diagram of a model used in the experiment;

FIG. 4 is a travel time comparison graph of different time algorithms of the DEC method and the shift hyperbolic method, the generalized time difference approximation method;

FIG. 5 is a time difference analysis comparison diagram of different travel time algorithms of the DEC method, the shift hyperbolic method and the generalized time difference approximation method.

Detailed Description

In order to make the objects, technical solutions and advantages of the embodiments of the present invention clearer, the technical solutions in the embodiments of the present invention will be clearly and completely described below with reference to the drawings in the embodiments of the present invention, and it is obvious that the described embodiments are some, but not all, embodiments of the present invention. All other embodiments, which can be derived by a person skilled in the art from the embodiments given herein without making any creative effort, shall fall within the protection scope of the present invention.

FIG. 1 is a schematic flow chart of a seismic travel time calculation method according to an embodiment of the present invention; as shown in fig. 1, S1, determining a basic travel time calculation formula basic form; the basic travel time calculation formula is an expression of the square of the basic travel time obtained by expanding according to the Taylor series principle and utilizing an anisotropic medium model; the basic travel time calculation formula is as follows:

T²(x)＝c₀+c₁x²+c₂x⁴+c₃x⁶+c₄x⁸+...

c₀＝t²

wherein x is an offset distance; c. C₀、c₁、c₂、c₃、c₄Coefficients that are Taylor series; t ═ T (0) is the time of spontaneous excitation and spontaneous reception, v ═ v_NMOIndicates NMO speed, S_k(k ═ 2,3,4) is a heterogeneity parameter;

s2, according to the Taylor series principle and the continuous fraction approximate travel time square formula, carrying out higher-order approximation on the basic travel time calculation formula to obtain the following formula:

wherein, the parameter B₁Is a kind of coefficient approximate expression based on 6 th order Taylor expansion and 4 th order Taylor expansion, parameter B₂Is based on a coefficient approximation of the 8 th order Taylor expansion and the 6 th order Taylor expansion;

s3, in x²Simultaneously performing power reduction processing on the 4 th power and the 6 th power of the offset distance in the higher-order basic travel time calculation formula for the cardinal number, and outputting travel time; the basic travel time calculation formula after the power reduction processing is as follows:

wherein, the common approximation can be taken

To simplify the parameter c in the formula₃、B₁、B₂The simplified parameters are as follows:

wherein S₂1+8 η, anisotropy parameterShowing the strength of anisotropy, showing the change in NMO velocity and reflection amplitude, the sum of the parameters being the Thomsen anisotropy parameter.

Specifically, FIG. 2 is a diagram of a planar horizontal layered media model; the basic travel time calculation formula is an expression of the square of the basic travel time obtained by expanding the formula according to the Taylor series principle and utilizing an anisotropic medium model; as shown in connection with fig. 2, shot S and receiver G are separated by an offset x. The travel time of a wave reflected from a shot point S to a receiving point G through N layers of interfaces is denoted as t (x), and can be represented by the following parameterized formula:

and

wherein the ray parameter p is p sin β_k/v_kTo give, β_kIs the angle of incidence of the ray on the bottom surface of the k-th layer, v_kLayer velocity of the k-th layer, Δ t_kIs the single pass travel time of layer k;

as can be seen from equation (2): for the first layer reflection, we can eliminate p and obtain the familiar hyperbolic travel-time formula:

wherein T (0) is the self-exciting and self-receiving time, x represents the offset distance, v₁Is the layer velocity of the first layer.

The travel time curve t (x) after multi-layer reflection can be represented in the form:

T²(x)＝c₀+c₁x²+c₂x⁴+c₃x⁶+...

wherein c is₀、c₁、c₂、c₃Is a coefficient of the Taylor series, c₀＝T²(0) T (0) is the self-excited self-receiving time, and x represents the offset.

Then, as can be seen from equation (1): the formula (1) can be expanded into the formula (4) by using Taylor series expansion,

or can be represented as

Wherein the content of the first and second substances,

from equation (5), we can derive x²，x⁴，x⁶Absolute convergence progression to the power of the equal power; for x²We get:

wherein the content of the first and second substances,

and is

..................

For x²ⁿComprises the following steps:

wherein the content of the first and second substances,

by the above method, we can deduce T²Power series expansion of (1):

wherein the content of the first and second substances,

J_j+1and V_(2j)The definition of (A) is the same as that of the formula (7) or (8); in addition, what is involved in the derivation procedure described above

J_j、J_j+1、V_(2j)Are all algebraic representations.

Finally, we take equation (9) and equation (13) into equation (3), and compare the power series of p on both sides of the equation to determine the coefficient c_i(i＝1,2,...)；

To write out T from the above theory²(x) As shown in formula (15):

T²(x)＝c₀+c₁x²+c₂x⁴+c₃x⁶+c₄x⁸+... (15)

wherein the content of the first and second substances,

c₀＝t²

in the above, T ═ T (0) is the self-excited time, and v ═ v is_NMOIndicates NMO speed, S_k(k ═ 2,3,4) is a heterogeneity parameter defined as:

secondly, carrying out higher-order approximation on the basic travel time calculation formula according to the Taylor series principle and a continuous fraction approximation travel time square formula; the continuous fraction approximation travel time squared formula that has been proposed in 2006 is:

we perform a higher order approximation of equation (18) according to the taylor series principle, as follows:

in the formula, the parameter B₁Is a kind of coefficient approximate expression based on 6 th order Taylor expansion and 4 th order Taylor expansion, parameter B₂Is based on a coefficient approximation of the 8 th order Taylor expansion and the 6 th order Taylor expansion;

finally, with x²4 and 6 times of offset in the formula for radix vs. higher order base travel timeThe method carries out simultaneous power reduction processing and outputs travel time; since Hubral stated in 1980: in the actual seismic data test result, some redundant information with little practical value can be obtained by the high-order terms of T (x). In theory, higher order terms have been bound by mathematically complex forms, and further approximation processing using higher order terms is generally not recommended.

Based on the above analysis, it is now necessary to handle the higher order terms in the formula. Now, the power-down thinking is used for processing, and the formula (21) is obtained:

wherein a general approximation can be taken according to the formula (17)

To simplify the parameter c in the formula (21)₃、B₁、B₂The simplified parameters are as follows:

in the formula, S₂Can continue to be represented by equations (23) - (24):

S₂＝1+8η (23)

where η is an anisotropy parameter, which can be expressed as a sum, where the strength of the anisotropy is expressed, the change in NMO velocity and reflection amplitude is expressed, and all are Thomsen anisotropy parameters.

The embodiment of the invention provides a seismic wave travel time calculation method, which comprises the steps of performing higher-order approximation on a basic travel time calculation formula, and then performing power reduction processing on the higher-order basic travel time calculation formula to output travel time; the operation is simple and the operation is efficient.

Further, the step S3 is followed by forming an exponentiation compensation formula for the base travel time calculation formula after the exponentiation processing in combination with an approximation compensation method:

wherein fb is the first arrival time, T_DE(x) And the travel time calculated by the calculation formula of the basic travel time after the power reduction is calculated, wherein n is the number of model layers.

In particular, from the point of view of the taylor series high order approximation, this "hard" processing method for the formula would change the meaning of "high" in the high order terms, leaving the calculated travel time less than fine accurate. Therefore, compensating for the calculated travel time is a way to compensate for the lack of detail, and an approximate compensation method is used herein as shown in equation (25):

finally, combining equation (21) with the compensation method can form a power down compensation (DEC) equation, denoted as DEC method, as shown in equation (26):

wherein fb is the first arrival time, T_DE(x) The travel time calculated for equation (21), n being the number of model layers.

For example, fig. 3 is a diagram of a model used in the experiment. As shown in FIG. 3, the experiment used a 13-layer horizontal laminar anisotropy model with track spacing of 40m, 226 receivers, maximum depth of 6000m, and maximum offset of 9000 m. The first arrival time is recorded by the detector above the first layer of the model.

FIG. 4 is a travel time comparison diagram of different time algorithms of the DEC method and the shift hyperbolic method and the generalized time difference approximation method. As shown in fig. 4, in which (a) a travel time calculation value at a depth of 800m and (b) a travel time calculation value at a depth of 2600 m; (c) travel time calculation value when the depth is 3500m (d) travel time calculation value when the depth is 4600 m;

FIG. 5 is a time difference analysis comparison diagram of different travel time algorithms of the DEC method, the shift hyperbolic method and the generalized time difference approximation method. As shown in fig. 5, in which (a) the time difference contrast is 800m in depth; (b) time difference contrast at a depth of 2600 m; (c) time difference comparison with a depth of 3500 m; (d) time difference contrast with a depth of 4600 m; the equation (27) is used in the time difference comparison experiment to verify the difference.

Wherein, the first arrival time fb is used as a reference value, and tc represents the time calculated by the travel time formula.

As can be seen from fig. 4: the DEC calculation method can effectively calculate the travel time as the other 2 methods. Although the calculated values of the DEC method show little difference at the near offset distance 2000m of 4 different depths, the advantages can be obviously seen at the middle offset distance 4500m and the far offset distance 8000 m; and the DEC method is better in travel time calculation comparison experiments.

As can be seen from fig. 5: the DEC method still shows small difference from other 2 methods in the vicinity of the near offset distance of 2000 m; considering the absolute value of the time difference calculated by the 3 methods, the time difference obtained by the DEC method is the minimum near the middle offset distance 4500m and the far offset distance 8000m, which also shows that the difference between the time difference and the first arrival time is smaller than that of other methods; the DEC method provides high calculation accuracy. The embodiment of the invention provides a seismic wave travel time calculation method, which comprises the steps of performing higher-order approximation on a basic travel time calculation formula, then combining the higher-order basic travel time calculation formula with a compensation method after power reduction processing to form power reduction compensation, and outputting travel time; the operation is simple, and the operation is efficient; the DEC method not only shortens the gap with the first arrival time in the calculation of the stratum travel time problem, but also improves the calculation accuracy.

Based on the above embodiments, an embodiment of the present invention provides a seismic travel time calculation system, including:

the calculation optimization module determines the basic form of the basic travel time calculation formula; the basic travel time calculation formula is an expression of the square of the basic travel time obtained by expanding according to the Taylor series principle and utilizing an anisotropic medium model; the basic travel time calculation formula is as follows:

T²(x)＝c₀+c₁x²+c₂x⁴+c₃x⁶+c₄x⁸+...

c₀＝t²

wherein x is an offset distance; c. C₀、c₁、c₂、c₃、c₄Coefficients that are Taylor series; t ═ T (0) is the time of spontaneous excitation and spontaneous reception, v ═ v_NMOIndicates NMO speed, S_k(k ═ 2,3,4) is a heterogeneity parameter;

and according to the Taylor series principle and the continuous fraction approximate travel time square formula, carrying out higher-order approximation on the basic travel time calculation formula to obtain the following formula:

wherein, the parameter B₁Is a kind of coefficient approximate expression based on 6 th order Taylor expansion and 4 th order Taylor expansion, parameter B₂Is a kind of coefficient approximate expression based on 8 th order Taylor expansion and 6 th order Taylor expansion

Calculation output module, with x²Simultaneously performing power reduction processing on the 4 th power and the 6 th power of the offset distance in the higher-order basic travel time calculation formula for the cardinal number, and outputting travel time; the basic travel time calculation formula after the power reduction processing is as follows:

wherein, the common approximation can be takenTo simplify the parameter c in the formula₃、B₁、B₂The simplified parameters are as follows:

wherein S₂1+8 η, anisotropy parameter

The anisotropy intensity is expressed, the NMO velocity and the change in reflection amplitude are expressed, and the Thomsen anisotropy parameters are all.

The embodiment of the invention provides a seismic wave travel time calculation system for executing the method, which is characterized in that higher-order approximation is carried out on a basic travel time calculation formula, then the higher-order basic travel time calculation formula is subjected to power reduction processing and then is combined with a compensation method to form power reduction compensation, and travel time is output; the operation is simple, and the operation is efficient; the DEC method not only shortens the gap with the first arrival time in the calculation of the stratum travel time problem, but also improves the calculation accuracy.

The above-described embodiments of the apparatus are merely illustrative, and the units described as separate parts may or may not be physically separate, and parts displayed as units may or may not be physical units, may be located in one place, or may be distributed on a plurality of network units. Some or all of the modules may be selected according to actual needs to achieve the purpose of the solution of the present embodiment. One of ordinary skill in the art can understand and implement it without inventive effort.

Through the above description of the embodiments, those skilled in the art will clearly understand that each embodiment can be implemented by software plus a necessary general hardware platform, and certainly can also be implemented by hardware. With this understanding in mind, the above-described technical solutions may be embodied in the form of a software product, which can be stored in a computer-readable storage medium such as ROM/RAM, magnetic disk, optical disk, etc., and includes instructions for causing a computer device (which may be a personal computer, a server, or a network device, etc.) to execute the methods described in the embodiments or some parts of the embodiments.

Finally, it should be noted that: the above examples are only intended to illustrate the technical solution of the present invention, but not to limit it; although the present invention has been described in detail with reference to the foregoing embodiments, it will be understood by those of ordinary skill in the art that: the technical solutions described in the foregoing embodiments may still be modified, or some technical features may be equivalently replaced; and such modifications or substitutions do not depart from the spirit and scope of the corresponding technical solutions of the embodiments of the present invention.