Surface velocity acquisition method and device based on dual-well microlog

Abstract

The embodiment of the application provides a method and device for obtaining surface velocity based on double-well micro logging. The method for obtaining surface velocity based on double-well micro logging comprises the following steps: obtaining first arrival time records of double-well micro logging in a target work area, wherein the double-well micro logging comprises a shooting well and a receiving well; dividing a surface layer of the target work area into a plurality of stratum units according to positions of each shooting point in the shooting well; each stratum unit corresponds to a shooting point; based on initial layer velocities of each stratum unit, iteratively calculating corrected layer velocities of each stratum unit by using a total travel time constraint equation and a horizontal projection distance constraint equation, and obtaining surface velocities of the target work area; wherein the layer velocity is a velocity of a seismic wave propagating in the stratum unit.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of exploration geophysics, and in particular to a method and device for obtaining surface velocity based on double-well micro-logging. BACKGROUND

[0002] Surface geological structure investigation is an important part of geophysical exploration, which directly affects the data acquisition, data processing and geological comprehensive interpretation of seismic exploration. The important results of surface geological structure investigation are the thickness and velocity of low-velocity layer, the thickness and velocity of velocity drop layer and the velocity of high-velocity layer, and a high-precision near-surface structure velocity model plays a crucial role in the calculation of datum surface static correction and the establishment of surface Q field.

[0003] In related technologies, there are various methods for obtaining shallow surface velocity, such as shallow refraction method, micro-logging method, VSP layer velocity inversion method, etc. However, for complex surface geological structure areas, there is a lack of shallow surface velocity calculation scheme with high precision and wide application range.

[0004] Therefore, how to accurately obtain the shallow surface velocity of the work area is a technical problem to be solved. SUMMARY

[0005] Embodiments of the present application provide a method and device for obtaining surface velocity based on double-well micro-logging, which can be used to accurately obtain the shallow surface velocity of the work area.

[0006] One of the embodiments of the present application provides a method for obtaining surface velocity based on double-well micro-logging, which comprises: obtaining the first arrival time record of double-well micro-logging in a target work area, wherein the double-well micro-logging comprises a shooting well and a receiving well; dividing the surface layer of the target work area into a plurality of stratum units according to the positions of each shooting point in the shooting well; each stratum unit corresponds to a shooting point; obtaining a total travel time constraint equation of a seismic wave according to the time taken by the seismic wave excited by any shooting point in the shooting well to pass through each stratum unit and the first arrival time corresponding to the seismic wave in the first arrival time record; obtaining a horizontal projection distance constraint equation according to the horizontal projection distance between the transmission point of the seismic wave excited by any shooting point in the shooting well and the shooting point and the horizontal projection distance between the shooting well and the receiving well; based on the initial layer velocity of each stratum unit, iteratively calculating the corrected layer velocity of each stratum unit by using the total travel time constraint equation and the horizontal projection distance constraint equation to obtain the surface velocity of the target work area; wherein the layer velocity is the velocity of the seismic wave propagating in the stratum unit.

[0007] In some embodiments, the total travel time constraint equation is as follows:

[0008]

[0009] wherein t d,j is the total travel time of the seismic wave excited by the target shot point to reach the geophone in the receiving well; j is the number of the target shot point; n is the number of the stratum units through which the seismic wave excited by the target shot point reaches the geophone in the receiving well; t i is the time taken by the seismic wave to pass through the i-th stratum unit, which is calculated by using the stratum unit travel time formula.

[0010] In some embodiments, the stratum unit travel time formula is as follows:

[0011]

[0012] wherein t i is the time taken by the seismic wave to pass through the i-th stratum unit, v i is the interval velocity of the i-th stratum unit, h i is the height of the i-th stratum unit; v1 is the interval velocity of the first stratum unit adjacent to the i-th stratum unit through which the seismic wave passes to enter the i-th stratum unit; θ1 is the incidence angle of the seismic wave when it passes out of the first stratum unit.

[0013] In some embodiments, the seismic wave horizontal projection distance constraint equation is as follows:

[0014]

[0015] wherein L is the horizontal projection distance between the shot well and the receiving well; n is the number of stratum units through which the seismic wave excited by the target shot point reaches the geophone in the receiving well; l i is the horizontal projection distance between the transmission point of the seismic wave excited by the target shot point when it passes out of the i-th stratum unit and the target shot point, which is calculated by using the stratum unit horizontal projection distance calculation formula.

[0016] In some embodiments, the stratum unit horizontal projection distance calculation formula is as follows:

[0017]

[0018] wherein l i is the horizontal projection distance between the transmission point of the seismic wave excited by the target shot point when it passes out of the i-th stratum unit and the target shot point, h i is the height of the i-th stratum unit; v1 is the interval velocity of the first stratum unit; v iis the initial layer velocity of the i-th stratum unit; θi is the incident angle when the seismic wave passes through the i-th stratum unit.

[0019] In some embodiments, the initial layer velocity of each stratum unit is used to iteratively calculate the corrected layer velocity of each stratum unit by using the total travel time constraint equation and the horizontal projection distance constraint equation, so as to obtain the surface velocity of the target work area; wherein the layer velocity is the velocity of the seismic wave propagating in the stratum unit, and the calculation formula of the curved ray horizontal projection is obtained according to the stratum unit horizontal projection distance calculation formula and Snell's law as follows:

[0020]

[0021] wherein s m is the target horizontal projection distance between the target curved ray passing through the target shot point and the target shot point; m is the number of stratum units; M is the total number of stratum units through which the seismic wave passes from the target shot point to the receiver in the receiving well; h i is the height of the i-th stratum unit; θ i is the incident angle when the seismic wave passes through the i-th stratum unit; according to the stratum unit travel time formula, the curved ray travel time calculation formula is obtained as follows:

[0022]

[0023] wherein t 0,m is the target total travel time of the seismic wave from the target shot point to the receiver; m is the number of stratum units; M is the total number of stratum units through which the seismic wave passes from the target shot point to the receiver in the receiving well; h i is the height of the i-th stratum unit; is the initial layer velocity of the i-th stratum unit; θ i is the incident angle when the seismic wave passes through the i-th stratum unit; v i-1 is the corrected layer velocity of the i-1-th stratum unit;

[0024] The initial layer velocity of each stratum unit is used to iteratively calculate the corrected layer velocity of each stratum unit according to the curved ray horizontal projection calculation formula and the curved ray travel time calculation formula.

[0025] In some embodiments, the iterative calculation of the corrected interval velocity of each stratum unit based on the initial interval velocity of each stratum unit according to the curved ray horizontal projection calculation formula and the curved ray travel time calculation formula comprises: taking each shot point as a target shot point from the bottom to the top of the shot well, for each target shot point, performing the following operations: adjusting the value of the incident angle θ i in the curved ray horizontal projection calculation formula, calculating the value of the target horizontal projection distance according to the value of the incident angle θ i after each adjustment by using the curved ray horizontal projection calculation formula until the value of the target horizontal projection distance satisfies the total travel time constraint equation of the seismic wave; adjusting the value of the initial interval velocity in the curved ray travel time calculation formula, calculating the value of the target total travel time according to the value of the adjusted initial interval velocity by using the curved ray travel time calculation formula; in the case that the value of the target total travel time does not satisfy the total travel time constraint equation of the seismic wave, re-performing the step of adjusting the value of the incident angle θ i in the curved ray horizontal projection calculation formula until the value of the target total travel time satisfies the total travel time constraint equation of the seismic wave; taking the initial interval velocity in the case that the value of the target total travel time satisfies the total travel time constraint equation of the seismic wave as the corrected interval velocity of the i-th stratum unit.

[0026] One of the embodiments of the present application provides a device for obtaining surface velocity based on double-well micrologging, which comprises: a first obtaining module, configured to obtain first arrival time records of double-well micrologging in a target work area, wherein the double-well micrologging comprises a shot well and a receiving well; a division module, configured to divide a surface of the target work area into a plurality of stratum units according to positions of each shot point in the shot well; wherein each stratum unit corresponds to a shot point; a second obtaining module, configured to obtain a total travel time constraint equation of a seismic wave according to a time taken by the seismic wave excited by any shot point in the shot well to pass through each stratum unit and a first arrival time corresponding to the seismic wave in the first arrival time records; a third obtaining module, configured to obtain a horizontal projection distance constraint equation according to a horizontal projection distance between a transmission point of the seismic wave excited by any shot point in the shot well and the shot point and a horizontal projection distance between the shot well and the receiving well; and a fourth obtaining module, configured to obtain a corrected interval velocity of each stratum unit based on an initial interval velocity of each stratum unit by using the total travel time constraint equation of the seismic wave and the horizontal projection distance constraint equation, and obtain a surface velocity of the target work area; wherein the interval velocity is a velocity of the seismic wave propagating in the stratum unit.

[0027] The embodiment of the present application provides an electronic device, the electronic device comprises a memory and a processor, the memory stores a computer program, and the processor executes the program to perform the method.

[0028] The embodiment of the present application provides a storage medium for storing a computer readable program, and the computer readable program is executed to perform the method.

[0029] The above technical solution provided by the embodiment of the present application has at least the following advantages compared with the prior art.

[0030] In the embodiment provided by the present application, the shallow layer of the target work area is divided into a plurality of stratum units according to the positions of the excitation points in the excitation well; each stratum unit corresponds to an excitation point; based on the initial layer velocity of the stratum unit, the corrected layer velocity of each stratum unit is iteratively calculated by using a total travel time constraint equation and a horizontal projection distance constraint equation, and the shallow layer velocity of the target work area is obtained. Therefore, the shallow layer velocity of the work area can be accurately obtained. BRIEF DESCRIPTION OF DRAWINGS

[0031] The present application will be further illustrated in the form of exemplary embodiments, which will be described in detail with reference to the accompanying drawings. These embodiments are not limiting, and in these embodiments, the same reference numbers represent the same structures, wherein:

[0032] Figure 1 is an exemplary flowchart of a method for obtaining a shallow layer velocity based on double-well micrologging according to some embodiments of the present application;

[0033] Figure 2 is an exemplary flowchart of a method for iteratively calculating the corrected layer velocity of each stratum unit according to the curved ray horizontal projection calculation formula and the curved ray travel time calculation formula according to some embodiments of the present application;

[0034] Figure 3a is an exemplary schematic diagram of double-well micrologging according to some embodiments of the present application;

[0035] Figure 3b is an exemplary schematic diagram of the propagation of seismic waves in double-well micrologging according to some embodiments of the present application;

[0036] Figure 4 is an exemplary schematic diagram of the seismic record picked up by the bottom and mouth hole detectors according to some embodiments of the present application;

[0037] Figure 5 is an exemplary schematic diagram of the first arrival picked up by the mouth hole detector of double-well micrologging according to some embodiments of the present application;

[0038] Figure 6ais an exemplary schematic diagram of multiple formation units according to some embodiments of the present application;

[0039] Figure 6b is an exemplary schematic diagram of rays formed by seismic waves passing through multiple formation units according to some embodiments of the present application;

[0040] Figure 7 is an exemplary schematic diagram of average velocity and iterative interval velocity calculation according to some embodiments of the present application;

[0041] Figure 8 is an exemplary schematic diagram of 52915229 double well microlog first arrival, average velocity, initial interval velocity in Huanxian area according to some embodiments of the present application;

[0042] Figure 9 is an exemplary flow chart of interval velocity update iteration according to some embodiments of the present application;

[0043] Figure 10a is an exemplary schematic diagram of seismic wave ray path acquired according to embodiments of the present application and seismic wave ray path acquired according to prior art;

[0044] Figure 10b is an exemplary schematic diagram of surface interval velocity curve of a work area acquired according to embodiments of the present application and surface interval velocity curve of the work area acquired according to prior art;

[0045] Figure 11 is an exemplary schematic diagram of surface interval velocity acquisition device based on double well microlog according to some embodiments of the present application;

[0046] Figure 12 is an exemplary structural schematic diagram of an electronic device according to some embodiments of the present application. DETAILED DESCRIPTION

[0047] In order to more clearly illustrate the technical solutions of the embodiments of the present application, the following will briefly introduce the drawings needed to be used in the embodiment description. Obviously, the drawings in the following description are only some examples or embodiments of the present application, and for those skilled in the art, the present application can also be applied to other similar scenarios without creative labor on the basis of these drawings. Unless it is obvious from the language environment or otherwise stated, the same reference numbers in the drawings represent the same structures or operations.

[0048] It should be understood that the "system", "device", "unit" and / or "module" used herein is a method for distinguishing different components, elements, parts, portions or assemblies at different levels. However, if other words can achieve the same purpose, the words can be replaced by other expressions.

[0049] As shown in the present application and claims, unless the context clearly indicates otherwise, the words "a," "an," "the," and / or "this" are not limited in scope to the singular, but rather include the plural. Generally, the terms "including," "includes," and "comprising," "comprises" are used in the sense of "including at least" and not the sense of "including only," such that these terms permit, but do not require, the steps or elements that follow the respective term to be present.

[0050] Flowcharts are used in the present application to illustrate the operations performed by the system according to embodiments of the present application. It should be understood that the preceding or following operations are not necessarily performed in sequence. Instead, the steps can be processed in reverse order or simultaneously. Other operations can also be added to these processes, or one or more steps can be removed from these processes.

[0051] Figure 1 is an exemplary flowchart of a method for obtaining surface velocity based on double-well micrologging according to some embodiments of the present application.

[0052] In step S110, the first arrival time record of double-well micrologging including the shooting well and the receiving well in the target work area is obtained.

[0053] Double-well micrologging technology is a relatively mature near-surface measurement technology. It can accurately understand the depth of the water table at the well location, the near-surface velocity structure, and provide a basis for the best shooting lithology and well depth design. At the same time, it can qualitatively analyze the absorption and attenuation of seismic waves in the near-surface, and provide help for accurately obtaining the datum surface static correction and establishing the surface Q field. Therefore, it is often used in the investigation of the surface geological structure of seismic exploration.

[0054] In the specific implementation process, double-well micrologging includes two wells penetrating the low velocity zone: the shooting well and the receiving well, and the interval between the two wells is about several meters.

[0055] As shown in the Figure 3a , multiple shooting points (shot points) are arranged in the shooting well, and the seismic wave is shot at intervals from the well bottom to the well mouth. The geophone is buried at the bottom of the receiving well for receiving downgoing waves, and the geophone is buried at the well mouth for receiving upgoing waves.

[0056] As shown in the Figure 3b , for the bottom geophone of the receiving well, the shooting point A shoots the seismic wave which propagates to the bottom geophone B through the underlying strata. For the well mouth geophone, the seismic wave propagates to the well mouth geophone C through the overlying strata.

[0057] In the specific implementation process, the first arrival time record of double-well micrologging including the shooting well and the receiving well in the target work area can be obtained by the following method.

[0058] The firing wells are fired in turn from deep to shallow, and each time a shot is fired in a firing well, the wellhead and the bottomhole geophone of the receiving well simultaneously record one trace, and a double-well microlog record is formed according to the firing sequence, as shown in FIG. 2. Figure 4 The first arrival time of the double-well microlog record represents the travel time of the seismic wave from the firing point position to the wellhead or the bottomhole geophone. In actual implementation, the first arrival travel time is picked up according to the negative first break point, and the first arrival time record of the double-well microlog is obtained, as shown in FIG. 3. Figure 5

[0059] In step S120, the surface layer of the target work area is divided into a plurality of stratigraphic units according to the positions of the firing points in the firing well.

[0060] In actual implementation, the double-well microlog data can be obtained and analyzed to obtain the well depth, well spacing, and firing point interval of the double-well microlog.

[0061] In order to maximize the simulation of the actual propagation path of the seismic wave and obtain a relatively accurate surface layer velocity of the target work area, in some embodiments, the surface layer of the target work area can be divided into a plurality of stratigraphic units according to the positions of the firing points in the firing well, as shown in FIG. 4, and each stratigraphic unit can be regarded as being composed of a horizontally uniform stratified viscoelastic medium. Figure 6a

[0062] As shown in FIG. 5, each stratigraphic unit corresponds to a firing point according to the firing sequence of the firing points in the firing well. Figure 6a

[0063] According to the ray theory, the seismic wave excited by each firing point in the firing well is received by the bottomhole geophone and the wellhead geophone in the receiving well, and needs to satisfy the following two conditions:

[0064] (1) The total travel time of the seismic wave excited by any firing point from the firing point to the geophone is equal to the first arrival time corresponding to the seismic wave in the first arrival time record.

[0065] (2) The sum of the horizontal projection distances of each ray path through which the seismic wave excited by any firing point travels from the firing point to the geophone is equal to the well spacing between the firing well and the receiving well.

[0066] Step S130 of the present specification establishes a seismic wave total travel time constraint equation according to condition (1).

[0067] Step S140 of the present specification establishes a horizontal projection distance constraint equation according to condition (2).

[0068] ​​​In step S130, the total travel time constraint equation for seismic waves is obtained based on the time taken for the seismic wave excited at any excitation point in the excitation well to travel through each stratigraphic unit and the first arrival time corresponding to the seismic wave in the first arrival time record.

[0069] In some embodiments, the total travel time constraint equation for seismic waves can be expressed by the following formula:

[0070]

[0071] Among them, t d,j The total travel time of the seismic wave excited at the target excitation point to the geophone located in the receiving well; j is the number of the target excitation point; n is the number of the formation units through which the seismic wave travels from the target excitation point to the geophone in the receiving well; t i The time taken for the seismic wave to pass through the i-th stratigraphic unit is calculated using the stratigraphic unit travel time formula.

[0072] The travel time of a seismic wave through the i-th stratigraphic unit can be calculated using the following formula:

[0073]

[0074] According to Snell's Law:

[0075] Substituting formula (3) into formula (2), we obtain the following formula for the travel time of a stratigraphic unit:

[0076]

[0077] In formulas (2) to (4), such as Figure 6b As shown: t i v is the time taken for the seismic wave to pass through the i-th stratigraphic unit. i h is the layer velocity of the i-th formation unit. i θ represents the height of the i-th stratigraphic unit; v1 represents the layer velocity of the first stratigraphic unit, which is adjacent to the i-th stratigraphic unit. Seismic waves pass through the first stratigraphic unit and then enter the i-th stratigraphic unit; i θ1 is the incident angle of the seismic wave when it exits the i-th stratigraphic unit, and θ2 is the incident angle of the seismic wave when it exits the first stratigraphic unit.

[0078] In summary, the total travel time constraint equation for seismic waves can be expressed as follows:

[0079]

[0080] In step S140, the horizontal projection distance constraint equation is obtained based on the horizontal projection distance between the transmission point and the excitation point when the seismic wave excited by any excitation point in the excitation well passes through each formation unit, and the horizontal projection distance between the excitation well and the receiving well.

[0081] In some embodiments, the horizontal projection distance constraint equation for seismic waves is as follows:

[0082]

[0083] Where L is the horizontal projected distance between the excitation well and the receiving well; n is the number of stratigraphic units traversed by the seismic wave excited at the target excitation point from the target excitation point to the detector in the receiving well; for example Figure 6b As shown, l i The horizontal projection distance between the transmission point of the seismic wave excited by the target excitation point and the target excitation point when it passes through the i-th stratigraphic unit is calculated using the formula for calculating the horizontal projection distance of the stratigraphic unit.

[0084] like Figure 6b As shown, according to ray theory, the seismic wave generated at point A travels through n horizontal strata to reach point B, with an offset distance of L at point A. n The depth of point B is H. n The path of a seismic wave through the i-th stratigraphic unit can be expressed by the following formula:

[0085]

[0086] The horizontal projection distance l between the transmission point of the seismic wave excited by the target excitation point and the target excitation point when it exits the i-th stratigraphic unit. i This can be expressed by the following formula:

[0087]

[0088] In summary, in some embodiments, the formula for calculating the horizontal projection distance of a stratigraphic unit can be specifically expressed as follows:

[0089]

[0090] Among them, h i v1 is the height of the i-th stratigraphic unit; v1 is the layer velocity of the first stratigraphic unit; v i θ1 is the layer velocity of the i-th stratigraphic unit; θ1 is the incident angle of the seismic wave when it exits the first stratigraphic unit.

[0091] In step S150, based on the initial layer velocity of the stratum unit, the corrected layer velocity of each stratum unit is iteratively calculated by using the total travel time constraint equation and the horizontal projection distance constraint equation, and the surface velocity of the target work area is obtained; wherein the layer velocity is the velocity of the seismic wave propagating in the stratum unit.

[0092] According to the principle of seismic exploration, based on the straight-line path, the first arrival time record of the dual-well microlog is the total travel time of the seismic wave from the shooting point to the wellhead or bottom geophone after propagating through different stratum units, then the average velocity of the upgoing wave above the upper boundary of the i-th stratum unit and the downgoing wave below the lower boundary of the i-th stratum unit is and which can be expressed as follows:

[0093]

[0094]

[0095] In the specific implementation process, for the upgoing wave, the average velocity of the n-th stratum unit can be taken as the initial layer velocity of the stratum unit; for the downgoing wave, the average velocity of the first stratum unit can be taken as the initial layer velocity of the stratum unit; the initial layer velocities corresponding to the upgoing wave and the downgoing wave of other stratum units can be obtained according to formula (12) and formula (13), that is, the initial layer velocities of the upgoing wave and the downgoing wave are and

[0096]

[0097]

[0098] Due to the difference in observation mode and ray path between the upgoing wave and the downgoing wave, the initial layer velocities of the stratum units of the upgoing wave and the downgoing wave are different, and the difference can be reduced by weighted average to obtain the initial layer velocity as shown below:

[0099]

[0100] In formula (10)-(14), L is the distance between the shooting well and the receiving well, H i is the depth of the i-th shooting point, t u,i is the first arrival time recorded by the wellhead geophone of the i-th shot, and t d,i is the first arrival time recorded by the well bottom geophone of the i-th shot.

[0101] In some embodiments, the following curved ray horizontal projection calculation formula can be obtained according to the stratum unit horizontal projection distance calculation formula and Snell's law:

[0102]

[0103] wherein s m is a target horizontal projection distance between a target bending ray of a seismic wave excited by the target excitation point and the target excitation point; m is a number of a stratum unit; M is a total number of stratum units through which the seismic wave passes from the target excitation point to a detector in a receiving well; h i is a height of the i-th stratum unit; θ i is an incident angle of the seismic wave when passing through the i-th stratum unit.

[0104] According to the stratum unit travel time formula, the following bending ray travel time calculation formula is obtained:

[0105]

[0106] wherein t 0,m is a target total travel time of the seismic wave excited by the target excitation point to reach the detector, m is a number of a stratum unit; M is a total number of stratum units through which the seismic wave passes from the target excitation point to a detector in a receiving well; h i is a height of the i-th stratum unit; is an initial layer velocity of the i-th stratum unit; θ i is an incident angle of the seismic wave when passing through the i-th stratum unit; v i-1 is a corrected layer velocity of the i-1-th stratum unit.

[0107] According to the bending ray horizontal projection calculation formula and the bending ray travel time calculation formula, the corrected layer velocities of the stratum units are iteratively calculated.

[0108] For detailed description of the corrected layer velocities of the stratum units being iteratively calculated according to the bending ray horizontal projection calculation formula and the bending ray travel time calculation formula, refer to the related content in Figure 2 , which will not be repeated here.

[0109] In the specific implementation process, according to the principle of seismic exploration, the corrected layer velocities of the upgoing wave and the downgoing wave when passing through each stratum unit can be calculated respectively based on the straight ray path according to the above method; for each stratum unit, the average corrected layer velocity of the upgoing wave and the downgoing wave when passing through the stratum unit is calculated, and the average corrected layer velocity is taken as the corrected layer velocity of the stratum unit.

[0110] In the embodiments provided in the present application, according to the positions of each excitation point in the excitation well of the dual-well micro-logging in the target work area, the surface layer of the target work area is divided into multiple stratum units; by using the first arrival time record of the dual-well micro-logging, the total travel time constraint equation and the horizontal projection distance constraint equation are established according to the ray theory; and according to the total travel time constraint equation and the horizontal projection distance constraint equation, the corrected layer velocity of each stratum unit is iteratively calculated. Since the scheme of using the bent ray to simulate the seismic wave propagation path is realized by dividing the stratum units, the corrected layer velocity of each stratum unit obtained is closer to the true value, thereby providing a precise basis model for surface velocity modeling, static correction, near-surface Q value estimation and the like.

[0111] Figure 2 is an example flowchart of a method for iteratively calculating the corrected layer velocity of each stratum unit according to the bent ray horizontal projection calculation formula and the bent ray travel time calculation formula according to some embodiments of the present application.

[0112] In the specific implementation process, each excitation point can be sequentially taken as a target excitation point from the bottom to the top of the excitation well, and for each target excitation point, the following operations are performed (see the detailed flow in Figure 9 ).

[0113] In step S210, the value of the incident angle θ i in the bent ray horizontal projection calculation formula is adjusted sequentially, and according to the value of the incident angle θ i after each adjustment, the value of the target horizontal projection distance is calculated by using the bent ray horizontal projection calculation formula, until the value of the target horizontal projection distance satisfies the total travel time constraint equation.

[0114] The initial layer velocity of the stratum unit can be obtained by formula (10)-(14).

[0115] For example only, the calculation of the layer velocity of the downgoing wave propagating in the second stratum unit adjacent to the first stratum unit shown in Figure 7 is described below.

[0116] According to Snell's law, in the case where the layer velocities of two adjacent layers are known, a bent ray AB2 simulating the path of the seismic wave passing through the initial stratum unit and the second stratum unit can be obtained according to formula (15), and the length of the horizontal projection s2 of the bent ray can be calculated by the following formula:

[0117] s2=h1tanθ1+h2tanθ2 (17)

[0118] wherein the initial value of θ1 is a preset initial incident angle θ 1,0 (through transmission to obtain a transmission angle θ 2,0 , i.e. θ2)

[0119] According to the horizontal projection distance constraint equation, the size relationship between s2 and the horizontal well distance L between the shot well and the receiver well is compared:

[0120] If the difference between s2 and L is greater than the horizontal distance error, a small angle increment Δθ is added to the incident angle θ1, s2 is recalculated, and the size relationship between s2 and L is compared until the difference between s2 and L is less than the horizontal distance error.

[0121] If the difference between s2 and L is less than the horizontal distance error, the horizontal projection distance constraint equation is satisfied at this time, that is, the horizontal projection distance of the curved ray formed by the seismic wave from the shot point to the geophone is equal to the well distance between the shot well and the receiver well.

[0122] In step S220, the value of the initial interval velocity in the curved ray travel time calculation formula is adjusted, and the value of the target total travel time is calculated by using the curved ray travel time calculation formula according to the adjusted value of the initial interval velocity.

[0123] According to formula (16), for the curved ray AB2 in the example of step S210, the corresponding target total travel time t 0,2 may be expressed as the following formula:

[0124]

[0125] where θ i is the adjusted incident angle in step S210.

[0126] The size relationship between t 0,2 and the first arrival time t d,2 in the double-well microseismic first arrival time record is compared, and if the difference between the two is less than the travel time error, it jumps to step S240 for execution.

[0127] If the difference between t 0,2 and t d,2 is greater than the travel time error, a small velocity update amount Δv is added to the initial interval velocity , that is,

[0128] In step S230, in the case where the value of the target total travel time does not satisfy the seismic wave total travel time constraint equation, the step of adjusting the value of the incident angle θ i of the curved ray horizontal projection calculation formula is re-executed until the value of the target total travel time satisfies the seismic wave total travel time constraint equation.

[0129] In step S240, the initial interval velocity in the case where the value of the target total travel time satisfies the seismic wave total travel time constraint equation is taken as the corrected interval velocity of the i-th stratum unit.

[0130] Continuing to describe the above example, in the case that the value of the target total travel time satisfies the seismic wave total travel time constraint equation, the initial layer velocity is corrected to obtain a corrected layer velocity of the second stratum unit.

[0131] After obtaining the corrected layer velocity of the second stratum unit, the operations in steps S210-S240 can be continued to be performed with the corrected layer velocity of the second stratum unit as a known quantity, until a corrected layer velocity of a third stratum unit is obtained. Figure 7 As shown in FIG. 6, the third stratum unit corresponds to a target shot point, and the corrected layer velocity of the second stratum unit is taken as a known quantity, the operations in steps S210-S240 are re-executed until a corrected layer velocity of the third stratum unit is obtained.

[0132] In a specific implementation process, in order to obtain more accurate corrected layer velocities of the stratum units, the corrected layer velocities of the stratum units corresponding to the downgoing wave and the upgoing wave can be calculated respectively according to steps S210-S240, and are denoted as and respectively for the convenience of distinction.Then, a more accurate layer velocity is obtained by weighted average:

[0133]

[0134] After the layer velocities of the stratum units are obtained by formula (19), the transmission time of the seismic wave in each stratum unit can be calculated according to the layer velocities of the stratum units, and a set Figure 9 is obtained as shown in FIG. 7. i i Further, the surface model can be established, the static correction calculation and the near-surface Q value estimation can be performed according to the set.

[0135] In the related art, the travel time of the seismic wave passing through the stratum is calculated by using the straight ray propagation path, which causes a large error in the establishment of the near-surface model.

[0136] In the embodiments provided in the present application, the stratum units are divided, the real curved ray path of the seismic wave passing through the surface is simulated to estimate the surface layer velocity of the target work area, so that a more reasonable and accurate surface layer velocity can be obtained. The scheme provided in the embodiments of the present application can greatly improve the accuracy of the near-surface model establishment, the static correction calculation and the near-surface Q value estimation, and provides an important basis and guarantee for further improving the quality of the seismic data processing.

[0137] The technical effects of the technical scheme of the present application are described below with specific application examples.

[0138] The 3D survey of Huanxian is located at the junction of Tianhuan depression in Ordos basin and the slope in northern Shaanxi, with an exploration area of 1007km 2 , and the overall stratum is relatively flat. The surface of the work area is loess mountainous area with crisscrossing gullies and large changes in topography, and the shallow velocity structure is complex. In order to further find out the variation law of the shallow velocity, the method for obtaining surface velocity based on double-well micrologging provided by the application is applied in this example. The following will take well 52915229 in the 3D work area of Huanxian as an example for specific description.

[0139] Step one: Obtain the relevant data of double-well micrologging in the work area and the first arrival time record.

[0140] The double-well micrologging well 52915229 is located in the middle of the work area, the maximum excitation well and the receiving well therein have a depth of 200 meters, the well spacing between them is 5 meters, the excitation point distance is 1-4 meters (the excitation point distance is equally spaced by 1 meter at the depth of 1-20 meters, the excitation point distance is equally spaced by 2 meters at the depth of 20-40 meters, and the excitation point distance is equally spaced by 4 meters at the depth of 40-200 meters) (as shown in Figure 3a ); the first arrival time of the wellhead detector gradually decreases by trace, and the first arrival time of the well bottom detector gradually increases by trace (as shown in Figure 3b ). Fine picking of the first arrival of the double-well micrologging lays the foundation for subsequent interval velocity calculation.

[0141] Step two: divide the surface layer of the work area into multiple stratum units according to the position of the excitation point.

[0142] The double-well micrologging excitation well and the receiving well have a spacing of 5 meters and a depth of 200 meters, and there are 70 excitation points. A horizontal layered medium model composed of multiple stratum units is established according to the point distance of the excitation points (as shown in Figure 6a );

[0143] Step three: establish the ray propagation constraint equation

[0144] The horizontal projection distance between the excitation well and the receiving well is 5m for calculation, and the total travel time constraint equation and the horizontal projection distance constraint equation are established by using the first arrival time record of the double-well micrologging.

[0145] Step four: calculate the initial interval velocity of the double-well micrologging stratum unit based on the straight ray theory.

[0146] Based on the straight ray theory, the average velocity from the excitation point to the wellhead detector and the well bottom detector is calculated according to the average velocity formula and the excitation depth and well spacing, which is used as the initial interval velocity of the stratum unit.

[0147] Step five: calculate the double-well micrologging interval velocity based on the bending ray and using the layer stripping method.

[0148] Based on the bending ray theory, a preset initial incident angle is set as 0°, an angle increment Δθ of the incident angle is 0.1°, a velocity updating amount Δv is 1 m / s, a travel time error is 2 ms, a horizontal distance error is 0.01 m, and according to a layer velocity acquisition procedure shown in Figure 9 , a corrected layer velocity of each stratum unit is obtained through iterative calculation.

[0149] A ray path of a seismic wave calculated by an embodiment of the present application from a shot point to a geophone and a corrected layer velocity of each stratum unit are compared and analyzed. A shot point A is used to excite a receiving point B, and a corrected ray path is a scatter point, and a path before correction is a straight line (as shown in Figure 10a ). A layer velocity curve calculated by the embodiment of the present application and a layer velocity curve calculated by a method used in the related art (as shown in Figure 10b ) are compared. A local layer velocity is corrected from 2097.86 m / s to 1122.52 m / s, a velocity error is 975 m / s, and the accuracy is higher.

[0150] Figure 11 is an exemplary schematic diagram of a surface layer velocity acquisition device based on double-well micrologging according to some embodiments of the present application.

[0151] As shown in Figure 11 , the surface layer velocity acquisition device based on double-well micrologging includes a first acquisition module 1110, a division module 1120, a second acquisition module 1130, a third acquisition module 1140, and a fourth acquisition module 1150.

[0152] The first acquisition module 1110 is configured to acquire a first arrival time record of double-well micrologging of a target work area including a shot well and a receiving well.

[0153] The division module 1120 is configured to divide a surface layer of the target work area into a plurality of stratum units according to positions of each shot point in the shot well; wherein each stratum unit corresponds to a shot point.

[0154] The second acquisition module 1130 is configured to obtain a total travel time constraint equation of a seismic wave according to a time used by the seismic wave excited by any shot point in the shot well to pass through each stratum unit and a first arrival time corresponding to the seismic wave in the first arrival time record.

[0155] The third acquisition module 1140 is configured to obtain a horizontal projection distance constraint equation according to a horizontal projection distance between a transmission point of the seismic wave excited by any shot point in the shot well and the shot point when the seismic wave passes out of each stratum unit, and a horizontal projection distance between the shot well and the receiving well.

[0156] The fourth obtaining module 1150 is configured to obtain the corrected interval velocity of each stratum unit by iteratively calculating the corrected interval velocity of each stratum unit based on the initial interval velocity of each stratum unit, the total travel time constraint equation and the horizontal projection distance constraint equation, and obtain the surface velocity of the target work area; the interval velocity is the velocity of the seismic wave propagating in the stratum unit.

[0157] In the embodiments of the surface velocity obtaining device based on the dual-well microlog, the specific processing of each module and the technical effects brought by the processing can be referred to the related description in the corresponding method embodiments, which will not be repeated here.

[0158] Figure 12 FIG. 1 is an exemplary structural schematic diagram of an electronic device according to some embodiments of the present application.

[0159] As shown in FIG. 1, the electronic device includes at least one processor 1201, at least one communication interface 1202, at least one memory 1203 and at least one communication bus 1204. Figure 12 The processor 1201 can be a processor CPU, or a specific integrated circuit ASIC (Application Specific Integrated Circuit), or one or more integrated circuits configured to implement embodiments of the present application. The memory 1203 can include a high-speed RAM memory, and can also include a non-volatile memory such as at least one disk memory. The memory 1203 stores a program, and the processor 1201 invokes the program stored in the memory 1203 to execute part or all of the method embodiments described above.

[0160] The present application relates to a storage medium for storing a computer readable program, which, when executed, performs part or all of the method embodiments described above.

[0161] Optionally, the storage medium can be a non-transitory computer readable storage medium, for example, the non-transitory computer readable storage medium can be a ROM, a random access memory (RAM), a CD-ROM, a magnetic tape, a floppy disk and an optical data storage device, etc.

[0162] Based on the same inventive concept, the embodiments of the present application also provide a computer program product including a computer program, which, when executed by a processor, implements part or all of the method embodiments described above.

[0163] Having described the basic concepts, it is obvious to those skilled in the art that the foregoing detailed disclosure is intended to be illustrative only and not limiting of the scope of the application. Although specific modifications, improvements, and alterations to the application have been discussed herein, other modifications, improvements, and alterations are also possible. Such modifications, improvements, and alterations are suggested to be within the spirit and scope of the application as described in the foregoing disclosure.

[0164] Also, the use of "an" or "one" to refer to an item of an embodiment also is understood as meaning "at least one" or "one or more" unless otherwise indicated. Similarly, the use of "another" to refer to an item of "at least one" or "one or more" is understood as meaning "at least one other" or "one or more other". Also, the use of "a" or "an" when used to modify the term "comprising" is understood to mean "comprising at least one of" or "comprising one or more of" unless otherwise indicated.

[0165] In addition, the order of processing elements or sequence of steps in which elements of the application are carried out, the use of numbers to refer to embodiments of the application, or the use of other designations by themselves, are not intended to limit the disclosure in any way. Although the foregoing disclosure has discussed various exemplary embodiments of the application, it should be apparent that modifications, improvements and alterations to the above-disclosed embodiments can occur to those skilled in the art. Such modifications, improvements and alterations are suggested to be within the spirit and scope of the application.

[0166] Similarly, it is to be noticed that the term "comprising", used in the description, does not exclude other elements or steps than the ones described, but it does not mean these optional elements or steps are essential, used in the claims, the use of the verb "comprise" means that other elements can also be present in the process, and they do not mean that they are essential to the technical solution described in each claim. Also, the term "one embodiment" is not used to refer to the same embodiment in different locations as used in the description and the claims unless otherwise indicated. The terms "the embodiment" or "an embodiment" are used to describe certain feature, structure, or characteristic in one or more embodiments. Thus, the embodiments are contemplated to encompass a variety of features, structures, or characteristics that can be implemented in any of the embodiments.

[0167] In some embodiments, numbers that describe amounts, dimensions, and so forth, are used in the description of the embodiments. It should be understood that such numbers are used only to illustrate certain embodiments and that the application is not limited to the numbers. In some examples, such numbers are modified by the modifier "about" or "approximately." Unless otherwise specified, "about" or "approximately" means ±20% of the value of the measured quantity that the term describes. Accordingly, in some embodiments, the numerical parameters in the description and claims are approximations that can vary depending upon the desired properties sought to be obtained by the individual embodiments. In some embodiments, numerical parameters are determined by the use of standard techniques. Although the numerical ranges and parameters setting forth the broad scope of the application are approximations, the numerical values set forth in the specific examples are reported as precisely as practicable. The numerical values set forth in the specific examples are provided to be as precise as reasonably possible. However, some variations may occur depending on the choice of input used to develop or derive the numerical values in the examples.

[0168] Each patent, patent application, publication, and other material cited in this application is hereby incorporated by reference in its entirety. In the event of inconsistencies between the disclosure of this application and the materials incorporated by reference, the disclosure of this application is intended to prevail. Nothing in this application is intended to be forgo the broadest scope of a claim in this application. It is noted that, as used in this application, the singular form "a", "an", and "the" include plural references unless the context clearly dictates otherwise. It is further noted that the use of "or" in the context of describing items in a list is intended to encompass all possible combinations of the items in the list.

[0169] Finally, it should be understood that the embodiments described herein are merely illustrative of the principles of this application. Numerous modifications and adaptations will be readily apparent to those skilled in the art. Accordingly, the scope of the application is not limited to the embodiments described herein, but rather is to be accorded the full scope of the claims.

Claims

1. A method for obtaining surface velocity based on dual-well micrologging, characterized in that, The method comprises: acquiring a first arrival time record of double-well micrologging in a target work area including a shooting well and a receiving well; dividing a surface layer of the target work area into a plurality of formation units according to positions of each shooting point in the shooting well; wherein each formation unit corresponds to a shooting point; obtaining a total travel time constraint equation of a seismic wave according to a time taken by the seismic wave excited by any shooting point in the shooting well to pass through each formation unit and a first arrival time corresponding to the seismic wave in the first arrival time record; obtaining a horizontal projection distance constraint equation according to a horizontal projection distance between a transmission point of the seismic wave excited by any shooting point in the shooting well and passing through each formation unit and a horizontal projection distance between the shooting well and the receiving well; iteratively calculating a corrected interval velocity of each formation unit based on an initial interval velocity of each formation unit, the total travel time constraint equation and the horizontal projection distance constraint equation to obtain a surface velocity of the target work area; wherein the interval velocity is a velocity of the seismic wave propagating in the formation unit.

2. The method of claim 1, wherein, The total travel time constraint equation of the seismic wave is as follows: wherein t d,j is the total travel time of the seismic wave excited by the target excitation point to reach the geophone located in the receiving well; j is the number of the target excitation point; n is the number of the stratum units through which the seismic wave travels from the target excitation point to reach the geophone in the receiving well; t i is the time taken by the seismic wave to pass through the i-th stratum unit, which is calculated by using a stratum unit travel time formula.

3. The method of claim 2, wherein, The formation unit travel time formula is as follows: where t i is the time taken by the seismic wave to pass through the i-th said stratum unit, v i is the interval velocity of the i-th said stratum unit, h i is the height of the i-th stratum unit; v1 is the interval velocity of a first stratum unit which is adjacent to the i-th stratum unit and through which the seismic wave passes before entering the i-th stratum unit; θ1 is the angle of incidence of the seismic wave as it exits the first stratum unit.

4. The method of claim 3, wherein, The horizontal projection distance constraint equation of the seismic wave is as follows: Wherein, L is the horizontal projection distance between the excitation well and the receiving well; n is the number of stratum units through which the seismic wave excited by the target excitation point reaches the geophone in the receiving well from the target excitation point; l i is the horizontal projection distance between the transmission point when the seismic wave excited by the target excitation point penetrates out of the i th stratum unit and the target excitation point, which is calculated by using the stratum unit horizontal projection distance calculation formula.

5. The method of claim 4, wherein, The formation unit horizontal projection distance calculation formula is as follows: wherein, l i is the horizontal projection distance between the transmission point when the seismic wave excited by the target excitation point penetrates through the i-th stratum unit and the target excitation point, h i is the height of the i-th stratum unit; v1 is the layer velocity of the first stratum unit; v i is the layer velocity of the i-th stratum unit; θ1 is the incidence angle when the seismic wave penetrates through the first stratum unit.

6. The method of claim 5, wherein, The iteratively calculating a corrected interval velocity of each formation unit based on an initial interval velocity of each formation unit, the total travel time constraint equation and the horizontal projection distance constraint equation to obtain a surface velocity of the target work area; wherein the interval velocity is a velocity of the seismic wave propagating in the formation unit, comprises: The curved ray horizontal projection calculation formula is obtained according to the formation unit horizontal projection distance calculation formula and Snell's law as follows: wherein s m is a target horizontal projection distance between a curved ray of a seismic wave excited from a target excitation point and the target excitation point, which the seismic wave excited from the target excitation point reaches a geophone; m is a number of a stratum unit; M is a total number of stratum units through which a seismic wave excited from a target excitation point reaches a geophone in a receiving well; h i is a height of an i-th stratum unit; θ i is an incident angle of a seismic wave when the seismic wave penetrates through the i-th stratum unit; The curved ray travel time calculation formula is obtained according to the formation unit travel time formula as follows: wherein t 0,m is the target total travel time for the seismic wave excited by the target excitation point to reach the geophone, m is the number of the stratum unit; M is the total number of stratum units through which the seismic wave travels from the target excitation point to reach the geophone in the receiving well; h i is the height of the i-th stratum unit; is the initial layer velocity of the i-th stratum unit; θ i is the incidence angle of the seismic wave when it passes out of the i-th stratum unit; v i-1 is the corrected layer velocity of the i-1-th stratum unit; The iteratively calculating a corrected interval velocity of each formation unit based on an initial interval velocity of each formation unit, the curved ray horizontal projection calculation formula and the curved ray travel time calculation formula.

7. The method of claim 6, wherein, The iteratively calculating a corrected interval velocity of each formation unit based on an initial interval velocity of each formation unit, the curved ray horizontal projection calculation formula and the curved ray travel time calculation formula, comprises: from the bottom to the top of the shooting well, each shooting point is sequentially taken as a target shooting point, and for each target shooting point, the following operations are performed: successively adjusting a value of the incident angle θ i in the curved ray horizontal projection calculation formula, calculating a value of the target horizontal projection distance by using the curved ray horizontal projection calculation formula according to a value of the incident angle θ i after each adjustment, until the value of the target horizontal projection distance satisfies the seismic wave total travel time constraint equation; adjusting a value of the initial interval velocity in the curved ray travel time calculation formula, calculating a value of the target total travel time by using the curved ray travel time calculation formula according to the adjusted value of the initial interval velocity; If the total travel time of the target does not satisfy the total travel time constraint equation for the seismic wave, the incident angle θ in the formula for calculating the horizontal projection of the curved ray is re-executed. i The steps of calculating the value of the target total travel time continue until the value of the target total travel time satisfies the seismic wave total travel time constraint equation; taking the initial interval velocity when the value of the target total travel time satisfies the total travel time constraint equation of the seismic wave as the corrected interval velocity of the i-th formation unit.

8. A device for obtaining surface velocity based on double-well micrologging, characterized in that, The device comprises: a first acquisition module configured to acquire a first arrival time record of double-well micrologging in a target work area including a shooting well and a receiving well; The division module is configured to divide the surface layer of the target work area into a plurality of stratum units according to the positions of the respective excitation points in the excitation well; each stratum unit corresponds to an excitation point; The second acquisition module is configured to obtain a total travel time constraint equation of a seismic wave according to the time taken by the seismic wave excited by any excitation point in the excitation well to pass through each stratum unit and the first arrival time corresponding to the seismic wave in the first arrival time record; The third acquisition module is configured to obtain a horizontal projection distance constraint equation according to the horizontal projection distance between the transmission point of the seismic wave excited by any excitation point in the excitation well when the seismic wave passes out of each stratum unit and the excitation point and the horizontal projection distance between the excitation well and the receiving well; The fourth acquisition module is configured to iteratively calculate a corrected stratum velocity of each stratum unit by using the total travel time constraint equation and the horizontal projection distance constraint equation based on an initial stratum velocity of each stratum unit, so as to obtain a surface velocity of the target work area; the stratum velocity is the velocity of the seismic wave propagating in the stratum unit. 9.An electronic device, comprising a memory and a processor, the memory storing a computer program, and the processor executing the program to perform the method of any one of claims 1 to 7. 10.A storage medium for storing a computer readable program, the computer readable program being executed to perform the method of any one of claims 1 to 7.

Images