CN112034519A - Iterative method for calculating average velocity of seismic wave rays in horizontal layered medium - Google Patents

Iterative method for calculating average velocity of seismic wave rays in horizontal layered medium Download PDF

Info

Publication number
CN112034519A
CN112034519A CN202011064021.1A CN202011064021A CN112034519A CN 112034519 A CN112034519 A CN 112034519A CN 202011064021 A CN202011064021 A CN 202011064021A CN 112034519 A CN112034519 A CN 112034519A
Authority
CN
China
Prior art keywords
ray
velocity
seismic
horizontal
offset
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Granted
Application number
CN202011064021.1A
Other languages
Chinese (zh)
Other versions
CN112034519B (en
Inventor
李启成
苑树鹏
郑新娟
闵也
贺翔
吴奎
肖艳妃
王明成
徐伊豪
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Liaoning Technical University
Original Assignee
Liaoning Technical University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Liaoning Technical University filed Critical Liaoning Technical University
Priority to CN202011064021.1A priority Critical patent/CN112034519B/en
Publication of CN112034519A publication Critical patent/CN112034519A/en
Application granted granted Critical
Publication of CN112034519B publication Critical patent/CN112034519B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
    • G01V1/28—Processing seismic data, e.g. for interpretation or for event detection
    • G01V1/30—Analysis
    • G01V1/303—Analysis for determining velocity profiles or travel times
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
    • G01V1/28—Processing seismic data, e.g. for interpretation or for event detection
    • G01V1/30—Analysis
    • G01V1/306—Analysis for determining physical properties of the subsurface, e.g. impedance, porosity or attenuation profiles
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V2210/00—Details of seismic processing or analysis
    • G01V2210/60—Analysis
    • G01V2210/62—Physical property of subsurface
    • G01V2210/622—Velocity, density or impedance
    • G01V2210/6222—Velocity; travel time

Landscapes

  • Engineering & Computer Science (AREA)
  • Remote Sensing (AREA)
  • Physics & Mathematics (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Acoustics & Sound (AREA)
  • Environmental & Geological Engineering (AREA)
  • Geology (AREA)
  • General Life Sciences & Earth Sciences (AREA)
  • General Physics & Mathematics (AREA)
  • Geophysics (AREA)
  • Geophysics And Detection Of Objects (AREA)

Abstract

The invention provides an iterative method for calculating the average speed of seismic wave rays in a horizontal layered medium, which comprises the steps of drilling a hole in an exploration area of a horizontal layered isotropic continuous medium and drilling a target layer; determining the layering thickness and the layering speed of the horizontal layered medium by using a seismic logging method; arranging a longitudinal measuring line, and determining the offset of a recording point; determining a ray constant through iterative operation; the ray average velocity is determined by the ray constant. At present, the average velocity in seismic exploration only reflects the change along with the depth, and the change of the average velocity in the horizontal direction is not considered, so that the exploration result generates errors due to the change of offset in the seismic exploration, and when the offset is large, the errors are often generated. According to the method, under the premise of determining the offset, iterative operation is carried out by taking a ray constant as a variable, and after the ray constant is determined, the average speed of the ray determined under the determined offset can be calculated. The average speed of the rays is closer to the real propagation speed of seismic waves, so that the exploration result is more accurate.

Description

Iterative method for calculating average velocity of seismic wave rays in horizontal layered medium
Technical Field
The invention belongs to the technical field of seismic exploration, and particularly relates to an iteration method for calculating the average velocity of seismic wave rays in a horizontal laminar medium.
Background
Seismic wave propagation velocity is one of the most important and fundamental parameters in seismic exploration. During time-depth conversion, the average speed is used, and the average speed is defined as the ratio of the total thickness of the layered medium to the travel time of seismic waves in the direction vertical to the stratum under the condition of a plurality of layers of media; when the processing such as migration is carried out, the average speed of the ray is used, and the definition of the average speed of the ray is the ratio of the total distance traveled by seismic waves along a certain ray to the used time; when dynamic correction processing is carried out or hyperbolic curve fitting is carried out on an actually measured time distance curve, the superposition speed is used, and the superposition speed is defined as a certain average value of different speeds of a certain wave recorded by various wave detection points in a row; in lithology interpretation, the interval velocity is used, and the interval velocity is defined as the propagation velocity of seismic waves in a uniform stratum. For the homogeneous medium with an inclined interface, the stacking speed is the equivalent speed; for a horizontal laminar medium, the stacking velocity is the root mean square velocity. The physical nature of the propagation velocity of seismic waves is a very complex parameter, and the propagation velocity is closely related to factors such as lithology, density, porosity, filling materials in pores, burial depth and geological age of rocks.
The seismic wave velocity can change in the horizontal direction and the vertical direction of propagation, in sedimentary rock areas, the stratum can be approximately regarded as a layered medium, and the seismic wave propagation velocity change along with the depth can be obtained through various methods such as seismic logging and the like. In seismic exploration, even in a sedimentary rock region of a laminated medium, the average speed of rays can be changed due to the change of the offset, if the change is ignored, the exploration result error can be caused, and when the offset is larger, the larger error can be caused.
With the exploitation of a large amount of oil and gas resources and the increase of the demand of the economic society on the oil and gas resources, the seismic exploration target is changed from a simple lithologic-structure oil and gas reservoir to a complex lithologic-structure oil and gas reservoir, a hidden oil and gas reservoir, a composite oil and gas reservoir and an unconventional oil and gas reservoir, and the main difficulty of the work is that the size of a target geologic body is reduced, the complexity and the heterogeneity are increased, and the horizontal velocity change is serious. The problem of inaccurate average speed caused by offset change is solved. The invention researches the problem that the average speed of seismic wave rays changes along with the offset in a horizontal layered medium.
Generally, the propagation velocities of seismic waves along an isotropic continuous medium are the same. In a horizontal laminar medium, the average velocity of the medium is generally calculated using equation (1). The formula is obtained by assuming that the seismic waves propagate perpendicular to the stratum, but in practical application, the average velocity in the horizontal laminar medium is calculated by using the formula (1) even if the seismic waves propagate along the direction which is not perpendicular to the stratum, the horizontal direction change of the average velocity of the seismic waves is not considered, and the average velocity is essentially taken as the ray average velocity, which is the main reason for the error of exploration in the horizontal laminar medium at present.
Figure RE-893414DEST_PATH_IMAGE001
（1）
In the formula (1), v is the average velocity of seismic waves in vertical stratum propagation, n is the number of strata, and h isiIs the thickness of the ith formation, ViIs the layer velocity of the ith layer. In seismic exploration, offset variations in the horizontal direction can cause seismic rays to vary in propagation path through the layers of the medium, resulting in changes in the mean ray velocity of the seismic waves throughout the earth. This ray mean velocity is related to the path traversed by the seismic ray. Establishing a medium model consisting of three horizontal layered stratums as shown in figure 1, wherein the parameters are layered thickness h1、h2And h3Layer velocity v1、v2And v3There are specific values in the figures. The model has three reflecting surfaces R1、R2And R3. At R1At the interface, the incident angles are respectively alpha1、α2And alpha3Etc.; at R2At the interface, the incident angles are respectively beta1、β2And beta3Etc.; at R3At the interface, the incident angles are respectively gamma1、γ2And gamma3And the like.
The average velocity in the three-layer medium was calculated to be 4286m/s according to equation (1). Accurately calculating the mean ray velocity requires calculating the refraction angles β and γ based on the known incidence angle α according to Snell's law. Snell's law content is shown in equation (2).
Figure RE-647743DEST_PATH_IMAGE002
（2）
P in equation (2) is a ray constant that depends on the initial propagation direction of the ray, regardless of the particular medium of the formation.
Figure RE-322438DEST_PATH_IMAGE003
Is the incident angle of seismic waves in the ith layer of medium,
Figure RE-845823DEST_PATH_IMAGE004
and n is the maximum stratum serial number. From equation (2), assume α1=10 °, β can be obtained1=16°32′，γ1=20 ° 10', which rays are incident obliquely downward and have respective propagation path lengths L in the three-layer medium1=1016m，L2=1044m，L3=1065m, the seismic ray from point O arrives at point S1The total distance traveled by the point was 6250m and the total time taken for propagation was 1.45s, and the mean velocity of the ray was calculated to be 4310 m/s. The offset x can be obtained according to the geometric relation in FIG. 11= 1684; similarly, when α is2When the angle is not less than 20 degrees, the average speed of the ray is 4420m/s and the offset x2=3977 m; when alpha is3When the radiation intensity is not less than 25 ℃, the average ray speed is 4560m/s, x3=6080 m; when alpha is4When the angle is not less than 27 degrees, the average speed of the ray is 4670m/s and the offset x4=7570 m; when alpha is5When the angle is not less than 29 degrees and is 30', the average speed of the ray is 5160m/s and the offsetx5=15458 m; when alpha is6When the angle is not less than 30 degrees, the average speed of the ray is 5450m/s, and the offset x6=27025m。
According to the above calculation, when the overall medium is not uniform, the average velocities of the rays propagating along different paths are different, the larger the offset distance is, the larger the average velocity of the seismic rays is, because the larger the offset distance is, the longer the path of the seismic rays propagating in the high-speed layer is, the larger the average velocity of the rays is, and the better the average velocity of the seismic waves can be represented only by considering the seismic wave propagation path, that is, the average velocity of the rays, which represents the velocity change with the depth in the layered medium and the velocity change in the horizontal direction, while the currently used average velocity only represents the velocity change with the depth.
Reference documents: land-based Mone, Wangyonggang, seismic exploration principle [ M ]. third edition, eastern Shandong, China Petroleum university Press, 2011, 144-.
Disclosure of Invention
Based on the fact that the problem of horizontal direction change is not considered in the process of calculating the average velocity in the horizontal layered medium at present, an iteration method for calculating the average velocity of seismic wave rays in the horizontal layered medium is provided. An iteration method for calculating the average velocity of seismic wave rays in a horizontal layered medium comprises the following specific processes:
step 1: drilling a hole in a horizontal layered medium exploration area to reach a target layer;
step 2: determining the layer thickness of a horizontal layered medium using seismic logging
Figure RE-728197DEST_PATH_IMAGE005
And speed of stratification
Figure RE-704244DEST_PATH_IMAGE006
I =1,2,3., n, i is the stratum serial number, and n is the maximum stratum serial number;
and step 3: arranging longitudinal measuring lines and determining the offset of recording points
Figure RE-917050DEST_PATH_IMAGE007
J =1,2,3.. m, j is the number of the recording points, m is the maximum order of the recording pointsNumber;
and 4, step 4: the calculation is iterated with the formula (3),
Figure RE-747472DEST_PATH_IMAGE008
the p value corresponding to the minimum value is the ray constant corresponding to the offset, the step length of the p value is 0.001, and the range is from 0 to sin (89 degrees) per v1，V1Is the minimum seismic interval velocity in the formation, typically the uppermost formation velocity, p is the ray constant,
Figure RE-82638DEST_PATH_IMAGE008
representing the iterative function of the jth recording point;
Figure RE-14822DEST_PATH_IMAGE009
（3）
and 5: substituting the P value into formula (4) to obtain a geophone offset of
Figure RE-93637DEST_PATH_IMAGE007
The ray mean velocity of the horizontal laminar medium;
Figure 100002_2
（4）
in the formula (4), the first and second groups,
Figure DEST_PATH_IMAGE002
the mean velocity variation in the horizontal and vertical directions in the horizontal laminar medium is taken into account at the same time for the ray mean velocity.
The theoretical derivation of the method for calculating the mean velocity of rays in a stratified homogeneous medium is described below with particular reference to FIG. 2. The ray travel time t in FIG. 2 is:
Figure RE-438533DEST_PATH_IMAGE012
（5）
Figure RE-321039DEST_PATH_IMAGE013
substituting the formula (2) into the formula (5) to obtain the formula (6) for the incident angle of the seismic wave in the ith layer of medium:
Figure RE-594894DEST_PATH_IMAGE014
（6）
when the formula (5) is substituted into the formula (2), the ray constant is taken as known, and the sine of the incident angle is calculated first, and then the cosine is calculated.
The total path s that the ray passes through is:
Figure RE-537442DEST_PATH_IMAGE015
（7）
the general expression of offset x is:
Figure RE-913060DEST_PATH_IMAGE016
（8）
the mean ray velocity can be found to be:
Figure RE-786207DEST_PATH_IMAGE017
（4）
in the formula (4), as long as p is obtained, the ray average speed can be calculated
Figure RE-258777DEST_PATH_IMAGE018
。
The ray constant P can be obtained by iterative computation using equation (3).
Figure RE-778751DEST_PATH_IMAGE019
（3）
Equation (3) is essentially derived from equation (8) because the ray constants are the only variables that can be solved for, but because implicit, are calculated iteratively. Formula (3) is the essence of the invention, firstly, the ray constant is determined as the only variable, then the ray constant is determined by an iteration method according to the offset, and further the ray average speed is calculated, which is a technical scheme not used at present. The current technical scheme is to determine the offset by using the incident angle, and the technical scheme is to determine the ray constant by using the offset, determine the incident angle and the like on the basis, and finally calculate the average speed of the ray.
In the formula (4), the first and second groups,
Figure RE-703981DEST_PATH_IMAGE007
j =1,2,. for the offset of the jth receiving point, m, j is the number of the recording point, and m is the maximum serial number of the recording point. Let VminIs the minimum seismic wave velocity in the formation, typically the uppermost formation wave velocity.
Figure RE-115240DEST_PATH_IMAGE020
Is an error function, and the p value takes 0.001 as a step size ranging from 0 to sin (89 °)/vminCalculating
Figure RE-442316DEST_PATH_IMAGE020
Taking the numerical value of
Figure RE-929929DEST_PATH_IMAGE008
And the ray constant P corresponding to the minimum value is the ray constant corresponding to the offset. Substituting the ray constant into equation (3) can obtain the average velocity of the ray.
The true ray mean velocity and the iteratively calculated ray mean velocity are listed in table 1.
TABLE 1 true ray average velocity and ray average velocity calculated by iterative method
Offset (m) True average velocity (m/s) P value Average velocity (m/s) obtained by iteration method
1684 4310 .5745E-04 4311
3977 4420 .1138E-03 4417
6080 4560 .1409E-03 4559
7570 4670 .1509E-03 4668
15458 5160 .1642E-03 5119
27025 5450 .1660E-03 5431
The first column in table 1 is the offset, the second column is the true ray average velocity at the corresponding offset, the third column is the ray constant P value obtained by iteration at the corresponding offset, and the fourth column is the ray average velocity calculated by using the P value. As can be seen from Table 1, the true ray mean velocity is very close to the ray mean velocity calculated by the iterative method, which preliminarily shows that the method for calculating the ray mean velocity in the horizontal lamellar medium proposed by us is effective.
The beneficial technical effects are as follows:
the invention provides an iterative method for calculating the average velocity of seismic wave rays in a horizontal layered medium. At present, the average velocity in seismic exploration only reflects the change along with the depth, and the change of the average velocity in the horizontal direction is not considered, so that the exploration result has errors due to the change of offset in the seismic exploration, and when the offset is larger, the larger error is often caused. The method firstly determines the ray constant as the only variable, then determines the ray constant by an iteration method according to the offset, and further calculates the ray average speed which is closer to the real propagation speed of seismic waves, so that the exploration result is more accurate.
Drawings
FIG. 1 is a schematic diagram of the calculation of horizontal lamellar ray average velocity according to the present invention;
FIG. 2 is a schematic diagram of the formula for deriving the mean ray velocity according to the present invention;
FIG. 3 is a theoretical model of a formation in a certain plain area according to an embodiment of the present invention;
FIG. 4 is a flow chart of an embodiment of the present invention;
FIG. 5 is a comparison of actual reflective surfaces and survey results according to the present invention;
FIG. 6 is a schematic diagram of the method for determining reflection points by bi-ellipse.
Detailed Description
The average speed of the seismic wave rays in the horizontal lamellar medium is calculated by an iterative method for the new generation stratum of a certain plain area, the depth of a target layer is explored, the embodiment of determining the average speed of the rays is explained, the exploration effect of the average speed of the rays is also researched, and the stratum profile is shown in figure 3.
An iterative method for calculating the average velocity of seismic rays in a horizontal laminar medium is shown in the flowchart of an embodiment of fig. 4. The specific process is as follows:
step 1: drilling a hole in a horizontal layered isotropic continuous medium exploration area, and drilling a target layer;
step 2: determining the layer thickness of a horizontal layered medium using seismic logging
Figure RE-342456DEST_PATH_IMAGE021
And speed of stratification
Figure RE-26247DEST_PATH_IMAGE006
Obtaining layered thickness and layered speed data, as shown in table 2, i =1,2,3, n, i is a stratum serial number, and n is a stratum maximum serial number;
and step 3: arranging longitudinal measuring lines and determining the offset of recording points
Figure RE-207830DEST_PATH_IMAGE007
J =1,2,3., m, j is the number of the recording point, and m is the maximum serial number of the recording point;
and 4, step 4: the calculation is iterated with the formula (3),
Figure RE-866344DEST_PATH_IMAGE022
the p value corresponding to the minimum value is the ray constant corresponding to the offset, and the p value takes 0.001 as the step length and ranges from 0 to sin (89 degrees) per vmin，VminIs the minimum seismic interval velocity in the formation, typically the uppermost formation velocity, p is the ray constant,
Figure RE-710976DEST_PATH_IMAGE022
representing the iterative computation function of the jth recording point;
Figure RE-277087DEST_PATH_IMAGE023
（3）
and 5: substituting the P value into formula (4) to obtain a geophone offset of
Figure RE-250859DEST_PATH_IMAGE007
The average velocity of the horizontal laminar medium;
Figure RE-877012DEST_PATH_IMAGE024
（4）
in the formula (4), the first and second groups,
Figure 964565DEST_PATH_IMAGE002
the calculated ray average velocities are shown in table 3 as ray average velocities.
The correctness of the method can be proved by calculating the average speed of the rays by using the law of refraction.
Further verification of the correctness of the ray average speed:
the average speed of seismic wave rays under different offset distances in a horizontal layered medium is calculated by an iterative method, but the influence of the average speed of the rays calculated by the iterative method on the exploration result is not proved. To further verify the effectiveness of the proposed method for calculating the mean ray velocity, the most ideal method is to separately explore the known subsurface reflection surface by the mean velocity calculation method adopted in the formula (1) at present and the proposed method for calculating the mean ray velocity by the iterative method, and compare the calculated mean ray velocity with the known reflection surface. However, the incorporability of the crust of earth determines that the underground reflecting surface cannot be mastered accurately, so that a theoretical model is established, the theoretical model is respectively explored by two methods for calculating the average speed, and exploration results are compared.
The new-generation stratum in a certain plain area is typical and can be regarded as a horizontal laminar medium, and the parameters of the medium with 9 layers on the surface are listed in the table 2. We build theoretical models with the parameters in Table 2, the first 8 layers are horizontal, the 9 th layer has dip, two methods are used to survey the 9 th layer for earthquake, and the survey results are compared with the actual results. Fig. 4 is a theoretical model of 9 strata in a certain plain area, the first 8 strata are horizontal, for convenience of drawing, all 8 strata are not drawn, and the 9 th stratum is an inclined stratum and is a target stratum for exploration. The key point of the verification is to calculate the theoretical seismic ray travel time as seismic exploration record.
TABLE 2 typical formation velocity parameters for new generation in certain plain area
Sequence of layers Depth of field Thickness (m) Lithology of stratum Layer velocity (m/s)
1 0 ̴2 2 The fourth series is surface soil, loose sand and loam 360
2 2 ̴5 3 Sand and clay of the fourth system 600
3 5 ̴15 10 Fourth series of water-containing sand and soil 1050
4 15 ̴50 35 The fourth is the upper part 1800
5 50 ̴200 150 The fourth part of the series 2000
6 200 ̴ 1000 800 Recent line N 2300
7 1000 ̴ 2000 1000 Superior in the near system E3 2800
8 2000 ̴ 4000 2000 Lower part of the ancient system E2 3500
9 4000 ̴ 6000 2000 E1 ̴ Ek transition layer 4500
Assuming that point S in fig. 3 is a seismic wave excitation point on the horizontal ground, an incident seismic wave SO is refracted by the multiple layers of horizontal layered media and then enters point O, the horizontal displacement of a seismic ray in the first 8 layers of horizontal layered media is KO, and the magnitude of KO can be calculated by equation (9):
Figure RE-820884DEST_PATH_IMAGE025
（9）
formula (9) has a similar meaning to formula (8)), with n =8 in formula (9). The difference between the formula (9) and the formula (8) is that the formula (8) includes the horizontal displacement of the incident wave and the reflected wave, and the formula (9) is only the horizontal displacement of the incident wave in the first 8 layers of horizontal layered media. P is calculated iteratively from equation (3), and n =9 in equation (3) during the iteration.
The horizontal stratum EC is 8 horizontal strata above, and is calculated by the formula (6) when rays travel in the strata, and n = 8; the layer 9 below the horizontal stratum EC is an inclined stratum with the inclination angle theta, as shown in FIG. 4. When the rays OB and BC travel need to be calculated, the total ray travel time is equal to the sum of the first 8 layer travel time and the 9 th layer travel time obtained by the formula (6).
Let the refraction angle λ of the point O in the 9 th inclined layer, the sine of the refraction angle is calculated by formula (10), and the normal depth OM of the point O is known, the length of OA can be calculated by formula (11).
Figure RE-649162DEST_PATH_IMAGE026
i=9 （10）
Figure RE-446217DEST_PATH_IMAGE027
（11）
Both OA and OM represent length. The refracted ray OB enters point B and reaches point C by reflection. Applying the sine theorem in Δ OAB yields:
Figure RE-242004DEST_PATH_IMAGE028
（12）
in equation (12), OB can be calculated from OA.
In Δ OBC:
Figure RE-415496DEST_PATH_IMAGE029
（13）
from equation (13), OC can be calculated. The wave velocity of the medium above the inclined reflection inclined plane below the EC plane is V9The travel time of the reflected wave in FIG. 3 below the EC plane and above the inclined reflection plane can be obtained by the equation (14)
Figure RE-629440DEST_PATH_IMAGE030
I.e. by
Figure RE-597396DEST_PATH_IMAGE030
Is the travel time of the seismic ray in the layer 9 medium.
Figure RE-614899DEST_PATH_IMAGE031
（14）
Thus, the formula (6) can be used for calculating the travel time of the seismic ray in the first 8 layers of horizontal strata
Figure RE-998607DEST_PATH_IMAGE032
And calculating the travel time T of the seismic ray in the 9 th layer medium by using the formula (14) to obtain the travel time T of the seismic ray in all the 9 th layer media.
Figure RE-129374DEST_PATH_IMAGE033
（15）
In the formula (15), the compound represented by the formula (I),
Figure RE-455182DEST_PATH_IMAGE034
is the travel of the seismic ray incident and reflected in the 1 st to 8 th layers of media, the formula is a specific application of formula (6), where n =8, and T is the travel of the seismic ray in all 9 layers of media.
Handle maleFormulas (6) and (14) into
Figure RE-38610DEST_PATH_IMAGE033
Obtaining:
Figure RE-226009DEST_PATH_IMAGE035
（16）
again, it is emphasized that in equations (15) and (16)
Figure RE-211282DEST_PATH_IMAGE036
n=8。
Given that the normal depth OM in fig. 4 is known, the EK length can be calculated by equation (17):
Figure RE-973571DEST_PATH_IMAGE037
（17）
KO in formula (17) has the meaning shown in FIG. 3. If the seismic ray from source S strikes point O ', the normal depth of the point O ' M ' is calculated using equation (18):
Figure RE-44295DEST_PATH_IMAGE038
（18）
the meaning of O 'M' in formula (18) is the same as OM in formula (17), and all represent the normal depth of the incident point of the layer 9 medium, but for different incident points, the determination method of the first incident point position is different from that of other incident points, and different letters are used for descriptive convenience. Similarly, KO' has the same meaning as KO in formula (17), as shown in FIG. 3. By knowing the normal depth of the incident point, the travel time of the reflected wave below the EC plane and above the inclined interface can be calculated by formula (14)
Figure RE-769806DEST_PATH_IMAGE030
And then calculates the travel time of the entire seismic ray.
Summarizing the steps of calculating the theoretical travel time of the 9 layers of media;
step one, establishing a model shown in figure 3, and obtaining medium layering thickness and layer velocity data through seismic logging, wherein the number of media in the model is 9, and the specific data are shown in table 2;
secondly, determining the offset between the excitation point and the receiving point
Figure RE-609586DEST_PATH_IMAGE039
J =1,2,3.. m, m is the maximum serial number of the recording point;
thirdly, iterating by using a formula (3) to obtain a ray constant P under the offset, wherein the number of strata in the formula (3) is n =9, and n is the maximum stratum serial number;
the fourth step, assuming that 9 layers of media are all horizontally layered, calculating the average speed in the 9 layers of media by using formula (4), wherein n = 9;
fifthly, calculating the horizontal displacement of the front 8 layers of horizontal stratums by using a formula (9), wherein n = 8;
and sixthly, calculating the refraction angle lambda in the 9-layer medium calculated by the formula (10), setting the inclination angle theta =10 degrees of the reflecting surface in the 9-layer medium, setting the normal depth OM of the first incidence point considered in the 9-layer medium to be known, and setting OM =1990m, and calculating the theoretical travel time of the seismic ray in the 9-layer medium as the seismic wave travel time of the seismic record by the formula (16).
Thus, the calculation of the average velocity of the first incidence point and the theoretical travel time of the seismic ray is completed, and the correlation calculation of other incidence points needs to recalculate the normal depth of the incidence point by using the formula (18). After the above-mentioned calculation of all recording points is completed, adopting existent method for calculating average speed in horizontal laminar isotropic continuous medium and method for calculating average speed of ray in horizontal laminar medium by using iteration method to respectively utilize double-ellipse method to make exploration, comparing the exploration result with actual one and analyzing error.
The average velocity obtained by calculation using the formula (1) was 3269 m/s. Table 3 shows the mean ray velocity calculated by equation (4) and the travel time calculated by equation (16) for different offsets.
TABLE 3 different shot-geophone distances shot-downLine average speed and travel time
Offset (m) Ray mean velocity (m/s) Traveling time(s)
100 .3269E+04 .3712E+01
200 .3269E+04 .3717E+01
300 .3269E+04 .3721E+01
400 .3269E+04 .3725E+01
500 .3269E+04 .3730E+01
600 .3270E+04 .3735E+01
700 .3270E+04 .3741E+01
800 .3270E+04 .3746E+01
900 .3270E+04 .3752E+01
1000 .3271E+04 .3759E+01
1200 .3271E+04 .3772E+01
1500 .3272E+04 .3794E+01
2000 .3275E+04 .3838E+01
3000 .3282E+04 .3947E+01
4000 .3291E+04 .4090E+01
5000 .3303E+04 .4269E+01
6000 .3318E+04 .4488E+01
7000 .3335E+04 .4754E+01
8000 .3354E+04 .5079E+01
9000 .3375E+04 .5478E+01
10000 .3397E+04 .5979E+01
12000 .3446E+04 .7498E+01
FIG. 5 is a survey result with lines and survey positions represented on the horizontal scale and depth represented on the vertical scale. Black bold line indicates the actual reflecting surface position; the chain line represents the position of the reflecting surface calculated by the average speed of the ray, and the average speed of the ray changes along with the offset; the black thin lines indicate the position of the reflecting surface obtained by the current average velocity, which does not vary with the offset. It can be seen that the position of the reflecting surface obtained by the average velocity of the ray is close to the actual reflecting surface no matter under a small offset or a large offset. And when the offset is small, the obtained position of the reflecting surface is in accordance with the actual position, but when the offset is large, the obtained reflecting surface has a great difference from the actual position. The maximum offset adopted in the calculation process is 12000m, and the position coordinates of the reflection point obtained by calculation do not exceed 4200m and are smaller than half of the maximum offset. It is explained that the reflection points are all concentrated in the tilt-up direction of the reflection surface.
The reason that errors exist in the positions of the reflecting surfaces obtained by exploration through an iterative method for calculating the average velocity of seismic wave rays in the horizontal layered medium is that the layered thicknesses of the layered medium in the horizontal direction are assumed to be consistent. We must make this assumption because we cannot increase the borehole density indefinitely. Nevertheless, the accuracy is significantly increased compared to the current methods of calculating average velocity.
The exploration process described above uses the double ellipse method we propose. For convenience of reading, the bi-elliptic method is briefly described as follows:
in reflected wave exploration, the coordinates of figure 6 are established in the ray plane, from the shot point
Figure RE-277196DEST_PATH_IMAGE040
Excited seismic waves, after reflection at the reflection point O
Figure RE-772900DEST_PATH_IMAGE041
The point is received and the travel time of the seismic wave is
Figure RE-629997DEST_PATH_IMAGE042
The wave velocity in medium 1 is v. The path length of the seismic ray from emission, reflection and final reception is vt1Possible reflecting surfaces are
Figure RE-245655DEST_PATH_IMAGE040
And
Figure RE-897216DEST_PATH_IMAGE043
as a focus, in vt1The ellipse is a constant elliptical surface, and the ellipse is set to ellipse 1. Of course, the ellipsoid should be located below the ground, i.e. the ordinate representing the position of the reflecting surface is negative. Then setting a slave shot point
Figure RE-880216DEST_PATH_IMAGE040
Excited seismic waves, after reflection at another reflection point
Figure RE-275425DEST_PATH_IMAGE044
Received, possible reflecting surfaces are on the ellipse 2 in fig. 1. The assumption that the reflecting surface is a plane with a fixed inclination is consistent with the understanding of the reflecting surface in current seismic exploration. The reflection points on the ellipses 1 and 2 are always tangent to the reflection surface at the same time, so that the reflection surface is in the ray plane and is always positioned on the common tangent line of the two ellipses. Two common tangent points below the ground are determined, namely the actual reflection points.
As with fig. 6, we re-assume for ease of discussion later. Order to
Figure RE-11169DEST_PATH_IMAGE040
As a shot point, from
Figure RE-99211DEST_PATH_IMAGE040
The seismic waves emitted by the points are reflected by the points O and arrive
Figure RE-303927DEST_PATH_IMAGE045
Is received with seismic wave travel time of
Figure RE-689778DEST_PATH_IMAGE046
，
Figure RE-358657DEST_PATH_IMAGE047
N is the recording point serial number, and N is the maximum recording point serial number. A reflecting surface exists, wherein a medium 1 is arranged above the reflecting surface, the wave velocity is V, and a medium 2 is arranged below the reflecting surface. The point of intersection of the seismic ray with the reflecting surface is not uniquely defined, the reflecting surfaces being respectively
Figure RE-289703DEST_PATH_IMAGE040
And
Figure RE-778454DEST_PATH_IMAGE048
is a focus point of
Figure RE-967995DEST_PATH_IMAGE049
Figure RE-491381DEST_PATH_IMAGE050
For a fixed-length ellipsoid, the possible reflection points are on the ellipse. A coordinate system is established as shown in FIG. 6, and horizontal coordinates are established on the vertical survey line of the common shot exploration
Figure RE-858908DEST_PATH_IMAGE040
And a receiving point
Figure RE-100534DEST_PATH_IMAGE048
The midpoint P of the connecting line is the origin of coordinates and the ordinate
Figure RE-562608DEST_PATH_IMAGE051
Indicating the ordinate of the reflecting surface in the ray plane. The major axis of the ellipse of the possible reflecting surface is anMinor axis of bnFocal length of cn。
The general form of the ellipse equation is:
Figure RE-940499DEST_PATH_IMAGE052
handle
Figure RE-478928DEST_PATH_IMAGE053
The shaft translates to
Figure RE-129221DEST_PATH_IMAGE054
Is to make as
Figure RE-473615DEST_PATH_IMAGE055
Axis, then the ellipse equation becomes:
Figure RE-909275DEST_PATH_IMAGE056
（19）
Figure RE-415343DEST_PATH_IMAGE057
Figure RE-84091DEST_PATH_IMAGE058
Figure RE-966596DEST_PATH_IMAGE059
equation (19) alone does not allow accurate determination of the reflection point. If the reflecting surface is assumed to be a continuous inclined surface and has a determined inclination angle, an equation is established by using two receiving points, and the accurate position of the reflecting point can be determined. As shown in FIG. 6, respectively
Figure RE-991184DEST_PATH_IMAGE060
And
Figure RE-668153DEST_PATH_IMAGE061
an equation is established for the focus point,
Figure RE-293038DEST_PATH_IMAGE061
is defined differently than for descriptive convenience
Figure RE-979234DEST_PATH_IMAGE048
The other receiving point of (2). Obtaining:
Figure RE-389487DEST_PATH_IMAGE062
（20）
Figure RE-424308DEST_PATH_IMAGE063
Figure RE-349539DEST_PATH_IMAGE064
Figure RE-511530DEST_PATH_IMAGE065
the letter meaning in equation (20) is similar to equation (19).
The reflecting surface is a continuous inclined plane with a fixed inclination angle between two reflecting points, and the linear equation of the reflecting surface in the seismic wave discovery plane is as follows:
Figure RE-838606DEST_PATH_IMAGE066
（21）
where K is the slope of the line and e is the intercept.
The straight line and the two ellipses respectively have an intersection point, namely the straight line is on the common tangent line of the two ellipses. Substituting equation (21) into equation (19) yields:
Figure RE-309908DEST_PATH_IMAGE067
（22）
Figure RE-456855DEST_PATH_IMAGE068
（23）
the straight line represented by equation (21) is tangent to the ellipse represented by equation (19) with only one intersection, and the discriminant at the root of the quadratic equation of unity (23) should be equal to 0, so:
Figure RE-422537DEST_PATH_IMAGE069
（24）
solving equation (24) yields:
Figure RE-604120DEST_PATH_IMAGE070
（25）
the reason for taking a negative value before the sign in equation (24) is that the reflection section is underground and its intercept is negative.
Similarly, the following formula (20) gives:
Figure RE-246323DEST_PATH_IMAGE071
Figure RE-146145DEST_PATH_IMAGE072
setting:
Figure RE-649939DEST_PATH_IMAGE073
（26）
and (3) setting the K value to be 0 to tg89 degrees and-tg 89 degrees to 0 to change, wherein the step size of K is 0.0001, and performing iterative operation on the equation (26) to obtain the slope of the reflecting surface, wherein K corresponding to the minimum value of m in the equation is obtained. When the value of K is substituted into the formula (25), the intercept e of the reflecting surface is obtained.
Solving equation (23) to obtain the horizontal coordinate of the first reflection point
Figure RE-686028DEST_PATH_IMAGE074
Substituting the value into equation (21) to obtain the corresponding ordinate
Figure RE-499132DEST_PATH_IMAGE075
。
Figure RE-823934DEST_PATH_IMAGE076
（27）
Figure RE-193736DEST_PATH_IMAGE077
（28）
Similarly, the horizontal coordinate of another reflection point can be obtained
Figure RE-271282DEST_PATH_IMAGE078
Ordinate and ordinate of the
Figure RE-68337DEST_PATH_IMAGE079
。
Figure RE-614856DEST_PATH_IMAGE080
（29）
Figure RE-522769DEST_PATH_IMAGE081
（30）
Accordingly, a method for determining the accurate position of a reflection point by using a double ellipse is provided, which comprises the following specific steps:
step 1: establishing a longitudinal measuring line, determining position coordinates of a shot point and a receiving point, and establishing a coordinate system by taking the measuring line as a horizontal coordinate axis, the depth as a longitudinal coordinate and the shot point as an origin of coordinates in a ray plane;
step 2: exciting seismic waves, and recording arrival times of the seismic waves to obtain a time distance curve of reflected waves;
and step 3: starting from the first receiving point, two adjacent receiving points are taken
Figure RE-720401DEST_PATH_IMAGE043
And
Figure RE-688357DEST_PATH_IMAGE044
and (3) substituting the recorded data into equation (8), setting the slope K value of the reflecting surface to be changed between 0 and tg89 degrees and-tg 89 degrees and 0, wherein the step size of K is 0.0001, and performing iterative operation on the equation (8) to obtain the slope of the reflecting surface corresponding to the minimum value of m in the equation. Substituting the K value into a formula (7) to obtain an intercept e of the reflecting surface; calculating the position coordinates of the two reflection points by equations (9) to (12) ((
Figure RE-722172DEST_PATH_IMAGE082
，
Figure RE-433776DEST_PATH_IMAGE083
），（
Figure RE-751494DEST_PATH_IMAGE084
，
Figure RE-890351DEST_PATH_IMAGE085
) (ii) a Taking two adjacent receiving points by the same principle
Figure RE-145883DEST_PATH_IMAGE086
And
Figure RE-848129DEST_PATH_IMAGE087
....,
Figure RE-833402DEST_PATH_IMAGE088
and
Figure RE-346423DEST_PATH_IMAGE089
substituting into related formula to calculate to obtain coordinate sequence of reflection point
Figure RE-151568DEST_PATH_IMAGE090
；
And 4, step 4: and drawing a reflecting surface position coordinate curve by using a coordinate sequence q.

Claims (2)

1. An iterative method for calculating the average velocity of seismic wave rays in a horizontal laminar medium is characterized by comprising the following specific processes:
step 1: drilling a hole in a horizontal layered medium exploration area to reach a target layer;
step 2: determining the layer thickness of a horizontal layered medium using seismic logging
Figure RE-DEST_PATH_IMAGE002
And speed of stratification
Figure RE-DEST_PATH_IMAGE004
I =1,2,3., n, i is the stratum serial number, and n is the maximum stratum serial number;
and step 3: arranging longitudinal measuring lines and determining the offset of recording points
Figure RE-DEST_PATH_IMAGE006
J =1,2,3., m, j is the serial number of the recording point, and m is the maximum serial number of the recording point;
and 4, step 4: iterative computation using equation (1), where p is the ray constant,
Figure RE-DEST_PATH_IMAGE008
representing the iterative function of the j-th recording point
Figure RE-666058DEST_PATH_IMAGE008
The p value corresponding to the minimum value is the ray constant corresponding to the offset, the step length of the p value is 0.001, and the range is from 0 to sin (89 degrees) and/vmin，VminThe minimum seismic wave layer velocity in the stratum is generally the uppermost stratum wave velocity;
Figure RE-DEST_PATH_IMAGE009
（1）
and 5: substituting the P value into the formula (2) to obtain the offset of
Figure RE-753837DEST_PATH_IMAGE006
The ray mean velocity of the horizontal laminar medium;
Figure 1
（2）
in the formula (2), the first and second groups,
Figure 2
vertical and horizontal velocity variations in the horizontal laminar medium are taken into account for ray mean velocity.
2. The iterative method for calculating the average velocity of seismic wave rays in a horizontal lamellar medium according to claim 1, characterized in that the ray constant to be solved in formula (1) is implicit and is used as an iterative operation for variables.
CN202011064021.1A 2020-09-30 2020-09-30 Iterative method for calculating average velocity of seismic wave rays in horizontal layered medium Active CN112034519B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202011064021.1A CN112034519B (en) 2020-09-30 2020-09-30 Iterative method for calculating average velocity of seismic wave rays in horizontal layered medium

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202011064021.1A CN112034519B (en) 2020-09-30 2020-09-30 Iterative method for calculating average velocity of seismic wave rays in horizontal layered medium

Publications (2)

Publication Number Publication Date
CN112034519A true CN112034519A (en) 2020-12-04
CN112034519B CN112034519B (en) 2022-09-30

Family

ID=73573046

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202011064021.1A Active CN112034519B (en) 2020-09-30 2020-09-30 Iterative method for calculating average velocity of seismic wave rays in horizontal layered medium

Country Status (1)

Country Link
CN (1) CN112034519B (en)

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20030208321A1 (en) * 2001-07-31 2003-11-06 Ruben Martinez Relative true amplitude migration
CN1611964A (en) * 2003-10-31 2005-05-04 中国石油化工股份有限公司 Method for determining underground speed structure for oil exploration
CN111239807A (en) * 2020-04-30 2020-06-05 辽宁工程技术大学 Method for determining accurate position of reflection point by using double ellipses

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20030208321A1 (en) * 2001-07-31 2003-11-06 Ruben Martinez Relative true amplitude migration
CN1611964A (en) * 2003-10-31 2005-05-04 中国石油化工股份有限公司 Method for determining underground speed structure for oil exploration
CN111239807A (en) * 2020-04-30 2020-06-05 辽宁工程技术大学 Method for determining accurate position of reflection point by using double ellipses

Non-Patent Citations (4)

* Cited by examiner, † Cited by third party
Title
刘骁: "高精度广角反射动校正时距关系方法研究", 《中国优秀博硕士学位论文全文数据库(硕士)基础科学辑》 *
李可恩: "含高速屏蔽层的地震数据采集及分析", 《中国优秀博硕士学位论文全文数据库(硕士)基础科学辑》 *
杜世通等: "迭偏速度的研究——基本定义和计算方法", 《华东石油学院学报》 *
田仁飞等: "射线平均速度对反射波时距方程简化的效果分析", 《内蒙古石油化工》 *

Also Published As

Publication number Publication date
CN112034519B (en) 2022-09-30

Similar Documents

Publication Publication Date Title
Cordier Velocities in reflection seismology
CN107783187B (en) Method for establishing three-dimensional velocity field by combining logging velocity and seismic velocity
CN106597533A (en) Depth domain velocity modeling method for piedmont zone seismic data processing
CN103513277B (en) Seismic stratum fracture crack density inversion method and system
CN106353793A (en) Cross-well seismic tomography inversion method on basis of travel time incremental bilinear interpolation ray tracing
CN104570102B (en) Method for combining near-surface velocity model with middle-deep stratum velocity model
CN103176211B (en) Based on gas-bearing reservoir prediction method and the device of many sensibility elasticities parameter
CN101634717A (en) Fine shear-wave (S-wave) impedance access technology based on logging and prestack channel set seismic data
CN105549082B (en) Method and system for establishing three-dimensional geomechanical field of ultra-deep carbonate reservoir
CN101487898A (en) Method for identifying oil, gas and water by longitudinal wave seismic exploration post-stack data
CN109738945A (en) A method of structural map is directly generated using pre-stack depth migration achievement
CN104834003B (en) Phased compression coefficient earthquake prediction method for unconventional tight reservoir
CN101630018A (en) Seismic exploration data processing method for controlling full acoustic wave equation inversion process
CN109470187A (en) Reservoir thickness prediction method based on three attribute of earthquake
CN104122582B (en) The method that accurately seismic velocity is asked for using stack velocity
CN103149588B (en) Method and system for calculating VTI anisotropic parameters by using well seismic calibration
CN107544093A (en) The structure interpretation layer depth system compensation method of borehole restraint
WO2015053659A1 (en) Method of producing an a priori hodograph for carrying out lithostratigraphic correlation
CN102901984A (en) Method for constructing true earth surface dip angle trace gathers of seismic data
CN110646840B (en) Angle gather extraction method and system
CN112034519B (en) Iterative method for calculating average velocity of seismic wave rays in horizontal layered medium
CN105510958A (en) Three-dimensional VSP observation system designing method suitable for complex medium
CN106842291A (en) A kind of unconformity trap reservoir lithology Forecasting Methodology based on pre-stack seismic ray Impedance Inversion
CN106990433B (en) Identification method for micro erosion groove in bump area
CN111596355B (en) Zero offset VSP time frequency analysis stratum division and layer velocity determination method

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant