Chromatographic velocity modeling method, device and equipment

Abstract

The invention relates to the field of petroleum geophysical exploration, in particular to a chromatography velocity modeling method, device and equipment. The tomography velocity modeling method comprises the following steps: acquiring depth migration imaging gathers of a plurality of azimuth angles of a research work area according to an initial tomography velocity model; determining a dominant azimuth angle based on the depth migration imaging gather of the plurality of azimuth angles of the research work area; constructing an azimuth angle weight matrix according to the dominant azimuth angle; and updating the initial chromatography velocity model according to the azimuth angle weight matrix. According to the embodiment of the invention, the dominant azimuth angles can be determined according to the depth migration imaging gathers of the plurality of azimuth angles of the research work area, the weight matrix is constructed according to the dominant azimuth angles so as to determine the tomography speed updating amount, more stable and higher-precision tomography speed inversion is realized, the tomography speed modeling efficiency is improved, and the modeling cost is reduced. And the actual production application of the complex exploration orientation distinguishing grid chromatography velocity inversion technology is promoted.

Description

Technical Field

[0001] The embodiments of the present invention relate to the field of petroleum geophysical exploration, and specifically to a tomographic velocity modeling method, device and equipment. Background Art

[0002] With the rapid development of the seismic exploration industry, in recent years, the seismic data acquisition method has developed from narrow azimuth data acquisition to multi-azimuth data acquisition, wide azimuth data acquisition, and even full-angle data acquisition. Some more complex geological structures have been discovered, and the accuracy requirements for exploration imaging in complex areas have also been continuously improved. However, in terms of modeling the tomographic velocity model of complex exploration areas, due to the complex surface and underground geological structures in complex exploration areas, the variable lithology and lithofacies of the formations, and other factors, the velocity distribution patterns are variable both vertically and horizontally. Conventional velocity modeling methods are difficult to achieve tomographic velocity modeling of non-uniform media.

[0003] Traditional azimuth grid tomography distinguishes the azimuth of input data related to complex exploration areas, can perform targeted processing on different azimuths during correlation analysis and ray tracing calculations, and can use the velocity azimuth anisotropy characteristics of the medium to update the tomographic velocity.

[0004] In the above-mentioned traditional azimuth grid tomography method, although the ray tracing calculation and the correlation analysis of the residual depth of the phase axis based on the offset gather are both performed based on different azimuths, in the subsequent solution of the tomography velocity equation, the solution areas of different azimuths will overlap. This reduces the inversion accuracy and modeling efficiency of the azimuth grid tomography, and seriously restricts the advantages of multi-azimuth and full-angle grid tomography technology in complex exploration areas and its production promotion and application. Therefore, how to improve the inversion accuracy and modeling efficiency of tomography velocity modeling and promote the advantages of multi-azimuth and full-angle grid tomography technology in complex exploration areas and its production promotion and application is a key issue that needs to be solved urgently. Summary of the invention

[0005] The purpose of the embodiments of this specification is to provide a tomographic velocity modeling method, device and equipment to achieve more stable and higher-precision tomographic velocity inversion, improve the efficiency of tomographic velocity modeling, and promote the practical production application of tomographic velocity inversion technology for complex exploration differentiation azimuthal grids.

[0006] On the one hand, an embodiment of the present specification proposes a tomographic velocity modeling method, which includes: obtaining depth migration imaging gathers at multiple azimuths of a study area based on an initial tomographic velocity model; determining a dominant azimuth based on the depth migration imaging gathers at multiple azimuths of the study area; constructing an azimuth weight matrix based on the dominant azimuth; and updating the initial tomographic velocity model based on the azimuth weight matrix.

[0007] On the other hand, a tomographic velocity modeling device comprises: an acquisition module for acquiring depth migration imaging gathers of multiple azimuths of a study area according to an initial tomographic velocity model; a determination module for determining a dominant azimuth based on the depth migration imaging gathers of multiple azimuths of the study area; a construction module for constructing an azimuth weight matrix according to the dominant azimuth; and an update module for updating the initial tomographic velocity model according to the azimuth weight matrix.

[0008] In another aspect, a computer device is provided, comprising a memory, a processor, and a computer program stored in the memory and executable on the processor, wherein the processor executes the above-mentioned tomographic velocity modeling method.

[0009] It can be seen from the technical solutions provided in the above embodiments of this specification that the embodiments of this specification obtain depth migration imaging gathers of multiple azimuths of the research area according to the initial tomographic velocity model; determine the dominant azimuth based on the depth migration imaging gathers of multiple azimuths of the research area; construct an azimuth weight matrix according to the dominant azimuth; and update the initial tomographic velocity model according to the azimuth weight matrix. Compared with the existing methods, the dominant azimuth is determined according to the depth migration imaging gathers of multiple azimuths of the research area, and the weight matrix is ​​constructed according to the dominant azimuth to determine the tomographic velocity update amount, which realizes a more stable and higher-precision tomographic velocity inversion, improves the efficiency of tomographic velocity modeling, and promotes the actual production application of azimuth grid tomographic velocity inversion technology for complex exploration. BRIEF DESCRIPTION OF THE DRAWINGS

[0010] In order to more clearly illustrate the embodiments of this specification or the technical solutions in the prior art, the drawings required for use in the embodiments or the description of the prior art are briefly introduced below.

[0011] Figure 1 is a flow chart of the tomographic velocity modeling method provided in the embodiments of this specification;

[0012] Figure 2 It is a schematic diagram of different ray paths of multi-azimuth and omni-directional seismic wave data at different azimuths and different perceptions of the inhomogeneous body when a local inhomogeneous body exists provided in an embodiment of this specification;

[0013] Figure 3 is a schematic diagram of a superimposed cross section divided into discrete grid spaces provided in an embodiment of this specification;

[0014] Figure 4 It is the theoretical tomographic velocity model of the three-dimensional salt dome model of region A provided in the embodiments of this specification;

[0015] Figure 5It is the initial tomographic velocity model of the region A containing a local high-speed body provided in the embodiments of this specification;

[0016] Figure 6 It is the depth migration imaging gather of 6 azimuth angles in area A provided in the embodiment of this specification;

[0017] Figure 7 is the residual depth curve corresponding to the first and third azimuths of region A provided in the embodiment of this specification;

[0018] Figure 8 It is a tomographic velocity field profile of the initial tomographic velocity model of the region A provided in the embodiments of this specification after being updated based on the tomographic velocity modeling method of the present invention and a tomographic velocity field profile of the initial tomographic velocity model after being updated based on the traditional multi-directional grid tomography;

[0019] Fig. 9 It is the tomographic velocity field of the initial tomographic velocity model after updating by the tomographic velocity modeling method of the present invention in the A region provided in the embodiments of this specification and the tomographic velocity field of the initial tomographic velocity model after updating based on the traditional multi-directional grid tomography;

[0020] Fig.10 It is a comparison diagram of the single-channel velocity curves of the heterogeneous body extracted from the theoretical chromatography velocity model of region A, the initial chromatography velocity model, the chromatography velocity model after one round of updating of the present invention, and the chromatography velocity model after one round of updating of the traditional multi-directional chromatography method provided in the embodiments of this specification;

[0021] Fig.11 The depth migration imaging gathers of the initial tomographic velocity model after updating by the tomographic velocity modeling method of the present invention in the A region provided in the embodiments of this specification and the depth migration imaging gathers of the initial tomographic velocity model after updating based on the traditional multi-azimuth grid tomography;

[0022] Fig.12 It is the vertical spectrum of the initial tomographic velocity model after updating by the tomographic velocity modeling method of the present invention in the A region provided in the embodiments of this specification, the vertical spectrum of the initial tomographic velocity model, and the vertical spectrum of the initial tomographic velocity model after updating based on the traditional multi-directional grid tomography;

[0023] Fig.13 It is a depth migration stacking section before and after the initial tomographic velocity model is updated in area B using the tomographic velocity modeling method of the present invention, as provided in the embodiments of this specification;

[0024] Fig.14 It is a schematic diagram of the structural composition of the tomographic velocity modeling device provided in the embodiments of this specification;

[0025] Fig.15 It is a schematic diagram of the structural composition of a computer device provided in an embodiment of this specification. DETAILED DESCRIPTION

[0026] The following will be combined with the drawings in the embodiments of this specification to clearly and completely describe the technical solutions in the embodiments of this specification. Obviously, the described embodiments are only part of the embodiments of this specification, not all of the embodiments. Based on the embodiments in this specification, all other embodiments obtained by ordinary technicians in this field without creative work should fall within the scope of protection of this specification.

[0027] Figure 1 The flow chart of the tomographic velocity modeling method is shown in FIG. When implemented, the method includes the following steps:

[0028] S101: Based on the initial tomographic velocity model, depth migration imaging gathers of multiple azimuths of the study area are obtained.

[0029] In some examples, based on the multi-azimuth and omni-directional seismic wave excitation and acquisition observation system, multi-azimuth and omni-directional seismic wave data related to the research area can be recorded, and the recorded multi-azimuth and omni-directional seismic wave data can be preprocessed to obtain time domain multi-azimuth and omni-directional seismic data related to the research area. Based on the time domain multi-azimuth and omni-directional seismic data related to the research area, an underground velocity model corresponding to the research area can be determined.

[0030] By collecting multi-azimuth and full-angle seismic wave data related to the study area and performing corresponding data preprocessing, multi-azimuth and full-angle seismic data in the time domain can be obtained, providing an important basis for subsequent geological interpretation and oil and gas exploration. Multi-azimuth and full-angle seismic data in the time domain not only contain detailed information on the underground geological structure of the study area, but also reflect the propagation characteristics of seismic waves in the study area in different directions and angles, which helps to more comprehensively understand the underground geological conditions of the study area.

[0031] Azimuthal anisotropy is a common problem in underground media. Figure 2 , Figure 2It is shown that when local heterogeneous bodies exist, the ray paths of multi-azimuth and omni-directional seismic wave data at different azimuths are different, and the perception of heterogeneous bodies is different. Multi-azimuth and omni-directional seismic wave excitation points and acquisition receiving points can be designed and arranged. Multi-azimuth seismic wave excitation points and acquisition receiving points mean that the seismic wave excitation points and acquisition receiving points are set in multiple different directions and angles, so that seismic wave data can be recorded from multiple directions and angles. Omni-directional seismic wave excitation points and acquisition receiving points mean that the seismic wave excitation points and acquisition receiving points are set in all possible directions and angles, so that seismic wave data can be recorded from all possible directions and angles. Explosives, air guns or other source equipment can be used to generate seismic waves at the excitation points. These seismic waves will propagate to the underground research area and interact with the underground medium in the research area. Seismic detectors or other sensors can be used at the acquisition receiving points to record the propagation of seismic waves. These sensors convert seismic waves into electrical signals and record them for subsequent processing.

[0032] The collected multi-azimuth and all-directional seismic wave data can be sorted out, including noise removal, time delay correction and superposition processing. At the same time, the integrity and accuracy of the data are verified to ensure data quality. Filtering technology can be applied to remove high-frequency noise and low-frequency interference in seismic wave data to improve the signal-to-noise ratio of seismic wave data. The seismic wave data can be statically corrected to eliminate the influence of surface topography and underground medium velocity changes on the propagation time of seismic waves. Seismic wave data from different excitation points and receiving points can be superimposed to enhance effective signals and suppress random noise. Each seismic wave data in the time domain multi-azimuth and all-directional seismic data is based on the time axis, showing the propagation and reflection characteristics of seismic waves in different strata.

[0033] In some embodiments, an initial tomographic velocity model can be established based on geological data, logging data and time domain migration velocity data; based on the seismic data of the study area and the initial tomographic velocity model, depth migration imaging gathers of multiple azimuths of the study area are generated.

[0034] An initial tomographic velocity model can be established through geological data, logging data and time domain migration velocity data, which lays the foundation for subsequent refined tomographic velocity modeling and tomographic velocity optimization.

[0035] The geological data comes from the preliminary stage of geological exploration, mainly including geological survey, geological mapping, geological profile measurement and other data. These data provide basic information such as geological structure, stratigraphic distribution, lithological characteristics, etc. related to the study area, and provide geological background for subsequent logging and time domain migration imaging. Well logging data are obtained through logging operations, mainly including resistivity, natural gamma, acoustic wave and other logging data. These data reflect the geophysical properties of the rock formations in the study area, such as conductivity, radioactivity, acoustic properties, etc. Time domain migration velocity data is obtained through seismic migration imaging technology, which reflects the velocity distribution of underground media in the study area. Time domain migration imaging technology includes post-stack time domain migration and pre-stack time domain migration, etc. These technologies can be applied to geological structures of different complexity. The geological data, logging data and time domain migration velocity data related to the study area can be obtained through literature research and public related databases.

[0036] After obtaining the geological data, logging data and time-domain migration velocity data related to the study area, the initial tomographic velocity model of the study area can be established based on these data. The geological data, logging data and time-domain migration velocity data can be integrated to analyze the information such as strata, lithology and velocity. This helps to understand and obtain the relationship between the geological structure and velocity distribution of the study area, and provides a basis for subsequent modeling. Specifically, according to the geological data, a geological framework is established, including information such as stratum distribution and fault location, so as to determine the basic form and scope of the initial tomographic velocity model. On the basis of the time-domain migration velocity data, combined with the velocity information in the logging data, the initial tomographic velocity model is preliminarily constructed. This can be achieved through interpolation, fitting and other methods to ensure the continuity and accuracy of the initial tomographic velocity model within the geological framework. After the initial tomographic velocity model is preliminarily constructed, it can be compared with the actual logging data, geological interpretation can be performed and corresponding preliminary corrections can be made. The initial tomographic velocity model after preliminary correction can better reflect the velocity distribution and geological structural characteristics of the underground medium in the study area.

[0037] In some embodiments, the initial tomographic velocity model of the study area can be subjected to multi-azimuth and omnidirectional prestack depth migration calculations based on the obtained time-domain multi-azimuth and omnidirectional seismic data of the study area, thereby generating depth migration imaging gathers of multiple azimuths of the study area.

[0038] By performing multi-azimuth and all-round pre-stack depth migration calculations on the initial tomographic velocity model, we can obtain depth migration imaging gathers at multiple azimuths of the study area, achieving all-round imaging of the underground structure of the study area, which helps to better understand the geological characteristics of the depth domain of the study area.

[0039] Based on the initial tomographic velocity model, the velocity change in the study area can be obtained, and the propagation path of the seismic wave and the anisotropy of the medium can be considered. The prestack depth migration algorithm can be used to calculate the multi-azimuth and all-directional seismic data in the time domain. The prestack depth migration algorithm can be the Kirchhoff integral method or the wave equation method. For example, the Kirchhoff integral method is used to calculate the seismic data in each azimuth time domain to generate the corresponding depth migration imaging gather. The Kirchhoff integral method is based on the Huygens principle, that is, each point on the wave front can be regarded as a secondary wave source, and the waves emitted by these secondary wave sources will form a new wave front at subsequent times. The underground structure is reconstructed by calculating the integral of the propagation path of the seismic wave in the underground medium. Specifically, the Kirchhoff integral method generates the corresponding depth migration imaging gather by calculating the amplitude and phase of the seismic wave at each receiving point in each azimuth time domain. For another example, the wave equation method is used to calculate the seismic data in each azimuth time domain to generate the corresponding depth migration imaging gather. The wave equation method is based on the principle that seismic waves propagating in underground media satisfy the wave equation, and the underground structure is reconstructed by solving the wave equation. Specifically, numerical methods such as the finite difference method, the finite element method, or the spectral element method are used to solve the wave equation corresponding to each azimuth time domain seismic data, and then generate the corresponding depth migration imaging gathers.

[0040] S102: Determine the dominant azimuth based on the depth migration imaging gathers of multiple azimuths in the research area.

[0041] In some embodiments, the depth-migrated imaging gathers of multiple azimuths of the study work area can be superimposed to obtain a stacked profile of the study work area; based on the stacked profile of the study work area, the velocity error sensitivity of the multiple azimuths of the study work area can be determined; based on the velocity error sensitivity of the multiple azimuths of the study work area, the dominant azimuth can be determined.

[0042] By determining the dominant azimuth among multiple azimuths, the foundation is laid for constructing an azimuth weight matrix based on the dominant azimuth, which in turn helps to improve the optimization accuracy and modeling efficiency of azimuth tomography speed.

[0043] The stacking section is to stack multiple geological sections together to form an overall geological model. By analyzing the changes in different layers on the stacking section, the changes in the underground geological structure can be more comprehensively inferred. The depth migration imaging gathers of multiple azimuths can be superimposed to obtain the stacking section of the study area. The stacking method can be arithmetic stacking, weighted stacking and phase stacking. Arithmetic stacking is to directly add the imaging gathers of each azimuth, which is suitable for situations with high data quality and relatively simple underground structures. Weighted stacking gives different weights according to the reliability and importance of each azimuth data to improve the accuracy of the stacking results. Phase stacking utilizes the coherence principle of seismic waves, that is, when seismic waves propagate in the underground medium, their phase information can reflect the characteristics of the underground structure. Phase stacking pays more attention to the phase information of seismic waves and is suitable for situations where the morphology and position of underground structures need to be accurately described. For example, the depth migration imaging gathers of multiple azimuths can be phase stacked. The depth migration imaging gathers of each azimuth can be phase-matched, that is, the phase information of the seismic wave can be corrected to ensure the phase consistency of the depth migration imaging gathers of each azimuth. The depth migration imaging gathers that have undergone phase matching are stacked. During the stacking process, weighted stacking can be performed according to the phase information to highlight the imaging effect of the underground structure. By stacking the depth migration imaging gathers of multiple azimuths, seismic data from different azimuths can be integrated, the lateral resolution and signal-to-noise ratio of seismic imaging can be improved, and the geological characteristics of the underground structure can be more accurately revealed.

[0044] After obtaining the stacked profile of the study area after stacking multiple azimuth depth migration imaging gathers, the sensitivity of each azimuth to the velocity disturbance of the abnormal body in the study area can be analyzed based on the stacked profile of the study area, that is, the velocity error sensitivity of each azimuth can be determined based on the stacked profile. A velocity error sensitivity threshold is set in advance. For each azimuth, if its velocity error sensitivity is greater than or equal to the preset velocity error sensitivity threshold, then this azimuth can be used as the dominant azimuth.

[0045] In some embodiments, the inclination field of multiple azimuths of the study work area can be calculated based on the stacked profile; the quality factor field of multiple azimuths of the study work area can be determined based on the inclination field of multiple azimuths of the study work area; the residual curvature of the depth migration imaging track set of multiple azimuths of the study work area can be calculated using similarity analysis based on the quality factor field of multiple azimuths of the study work area; the velocity error sensitivity of multiple azimuths of the study work area can be determined based on the residual curvature of the depth migration imaging track set of multiple azimuths of the study work area.

[0046] By determining the velocity error sensitivity in multiple azimuths, the foundation is laid for determining the dominant azimuth based on the velocity error sensitivity of the azimuth and constructing the azimuth weight matrix, which in turn helps to improve the optimization accuracy and modeling efficiency of the azimuth tomography velocity.

[0047] The superimposed cross-section of the study area can be divided into discrete grid spaces. Figure 3 , Figure 3 The schematic diagram of the superimposed profile of the study area is shown as a discrete grid space. It can be agreed that the grid index parallel to the main survey line direction is denoted as i=1,...,I, the grid index orthogonal to the main survey line direction is denoted as j=1,...,J, and the grid index in the depth direction is denoted as k=1,...,K. The inclination field of each azimuth of the study area describes the degree of inclination of the underground geological structure of the study area at each azimuth relative to the horizontal plane. The inclination field of the study area includes the inclination of each discrete grid in the superimposed profile of the study area. By calculating the inclination of each discrete grid in the superimposed profile of the study area at each azimuth, the inclination field corresponding to each azimuth can be obtained.

[0048] When the azimuth is determined, the local scanning algorithm can be used to determine the inclination corresponding to each discrete grid in the study area. Specifically, for each discrete grid in the study area, a scanning window centered on the discrete grid can be set. The scanning window includes multiple discrete grids. Scan lines of different angles are emitted at the center of the scanning window. These scan lines can be straight lines, curves, or other shapes, depending on the geological structure characteristics of the superimposed profile, which will not be described here. For each scan line, the characteristic values ​​of the pixels or voxels it passes through are extracted. These characteristic values ​​can be amplitude, frequency, phase, etc. Compare the similarity between the extracted characteristic values ​​on each scan line and the reference characteristic values ​​(which can be adjacent scan lines, known geological structures or models, etc.). The similarity can be measured by calculating the Euclidean distance, correlation coefficient, mutual information, and other metrics between the characteristic values.

[0049] For example, when the azimuth angle is determined, for each scan line, the amplitude value of each discrete grid through which the scan line passes is extracted. The amplitude value of the discrete grid can be the amplitude value of the center point of the discrete grid. Calculate the similarity between the amplitude values ​​of all discrete grids involved in the scan line and the amplitude values ​​of all discrete grids in the scan window, and take the angle corresponding to the maximum similarity as the inclination angle of the discrete grid. n is the nth azimuth; (i,j,k) represents the discrete grid with index (i,j,k) in the study area; B is the number of discrete grids in the scanning window; A tis the amplitude corresponding to the discrete grid. Repeat the above steps for each discrete grid in the stacked section, calculate the inclination of each discrete grid, and then determine the inclination field corresponding to the current azimuth of the study area. Repeat the above steps for each azimuth to determine the inclination field corresponding to multiple azimuths of the study area.

[0050] After determining the inclination fields of multiple azimuths in the study area, the quality factor fields of multiple azimuths in the study area can be determined. The quality factor field of the study area includes the quality factor values ​​of each discrete grid in the superimposed profile of the study area. By calculating the quality factor values ​​of each discrete grid in the superimposed profile of the study area at each azimuth, the quality factor field corresponding to each azimuth can be obtained. When the azimuth is determined, the quality factor value of each discrete grid in the superimposed profile of the study area can be the maximum similarity value corresponding to the inclination of each discrete grid. Repeat the above steps for each azimuth to determine the quality factor field of multiple azimuths in the study area.

[0051] Based on the quality factor field of multiple azimuths in the study area, a quality threshold can be set in advance, and discrete grids with quality factor values ​​greater than or equal to the quality threshold can be selected as seed points. The residual curvature of the depth migration imaging gather at each azimuth in the study area includes the residual curvature of various sub-points in the stacked section of the study area. By calculating the residual curvature of all seed points in the stacked section of the study area at each azimuth, the residual curvature of the depth migration imaging gather at each azimuth in the study area can be obtained. According to the stacked section of the study area, a hyperbola model can be established using the corresponding geological data. The hyperbola model is used to describe the velocity change and reflection interface morphology of seismic waves propagating in underground media. The establishment of the hyperbola model requires the determination of some key parameters, including the propagation velocity of seismic waves in the stacked section of the study area, the morphology and position of the reflection interface, etc. The geological data provides information on the physical properties and geometric morphology of the underground medium, which is the basis for establishing the hyperbola model. Specifically, the relevant information in the geological data can be input into the hyperbola model, and the reflection wave event axis in the stacked section can be fitted by adjusting the model parameters. After multiple iterations and adjustments, the best hyperbola model can be obtained. This hyperbola model can be used to match data from depth migration imaging gathers at multiple azimuths to find the most similar hyperbola. Specifically, when the azimuth is determined, for each seed point in the study area, a local scanning window is defined with the seed point as the center. Features are extracted from the depth migration imaging gathers within the scanning window, and a hyperbola is fitted based on the extracted features. The extracted features can be amplitude, frequency, phase, etc. The hyperbola fitted based on the extracted features is matched with the hyperbola model to calculate the similarity. The similarity can be measured by calculating the Euclidean distance, correlation coefficient, mutual information, and other metrics between the feature values. Traverse all possible hyperbolas and find the hyperbola with the highest similarity as the remaining depth curve of the seed point.

[0052] For example, when the azimuth angle is determined, for each scan line, the amplitude value of each seed point that the scan line passes through is extracted. Calculate the similarity between the amplitude values ​​of all seed points involved in the scan line and the amplitude values ​​of all seed points in the scan window, and take the hyperbola corresponding to the maximum similarity value as the residual depth curve of the seed point. n is the nth azimuth; (i', j', k') represents the seed point with index (i', j', k'); B' is the number of seeds in the scanning window; A t is the amplitude corresponding to the seed point. Repeat the above steps for each seed point in the stacked profile to obtain the residual depth curve of each seed point. Take the second derivative of the residual depth curve of each seed point to obtain the residual curvature corresponding to each seed point, and then determine the residual curvature corresponding to the current azimuth. Repeat the above steps for each azimuth to determine the residual curvature of multiple azimuth depth migration imaging gathers in the study area.

[0053] Based on the residual curvature of the depth migration imaging gathers at multiple azimuths in the study area, the velocity error sensitivity at multiple azimuths in the study area can be determined. The residual curvature can characterize the curvature deviation of the depth migration imaging gathers obtained based on the initial tomographic velocity model relative to the ideal plane or reference plane. The event axis curvature deviation of the depth migration imaging gathers of the seed points in the stacked profile of the study area at each azimuth characterizes the sensitivity of each azimuth to the velocity disturbance of the abnormal body, that is, the velocity error sensitivity. The residual curvature at the seed point corresponding to each azimuth can be used as the velocity error sensitivity at the seed point, and the velocity error sensitivity of each azimuth can be obtained by fusion based on the velocity error sensitivities of all seed points corresponding to each azimuth. The fusion method can be to calculate the mean, weighted mean, variance, standard deviation, etc. For example, the standard deviation of the velocity error sensitivity of all seed points corresponding to each azimuth can be calculated as the velocity error sensitivity of each azimuth.

[0054] S103: constructing an azimuth weight matrix according to the dominant azimuth.

[0055] In some embodiments, the fracture breakpoint responsiveness of the study area at multiple azimuths can be calculated based on geological data; the dominant azimuth is determined based on the velocity error sensitivity and fracture breakpoint responsiveness of the study area at multiple azimuths; and the azimuth weight matrix is ​​constructed based on the velocity error sensitivity and fracture breakpoint responsiveness of the dominant azimuth.

[0056] By constructing an azimuth weight matrix based on the velocity error sensitivity of the dominant azimuth and the fracture breakpoint responsiveness, the sensitivity of different azimuths to the initial tomographic velocity model error is more reasonably and comprehensively characterized, which helps to achieve higher-precision tomographic velocity modeling and optimization.

[0057] According to the geological data related to the research area, the fracture breakpoint responsivity information of each azimuth can be obtained. The fracture breakpoint responsivity of each azimuth includes the relationship information between the azimuth and the fracture direction and the width information of each azimuth. When the azimuth is consistent or nearly consistent with the fracture direction, the fracture is displayed more clearly on the imaging section, and the fracture breakpoint responsivity of the azimuth obtained according to the geological data will be higher. This is because the propagation characteristics of seismic waves in the fracture direction may make the fracture easier to identify at a specific azimuth. On the contrary, when the azimuth is perpendicular or nearly perpendicular to the fracture direction, the fracture is displayed more blurred on the imaging section, and the fracture breakpoint responsivity of the azimuth obtained according to the geological data will be lower. This is because in this case, the propagation characteristics of seismic waves may make it difficult to identify fractures. In the case of wide azimuth (for example, the aspect ratio of the research area is greater than 0.5), seismic exploration can obtain higher spatial imaging resolution, which is helpful to identify fracture breakpoints. At this time, the fracture breakpoint responsivity of wide azimuth obtained according to the geological data will be higher. This is because wide-azimuth seismic data can provide more information about underground structures and help reduce ambiguity on imaging sections. In contrast, narrow-azimuth seismic exploration may be limited by the complexity of underground structures and the propagation characteristics of seismic waves, resulting in a more blurred display of fracture breakpoints on imaging sections. At this time, the narrow-azimuth fracture breakpoint response obtained based on geological data will be lower. The relationship information between each azimuth and the fracture direction in the geological data and the width of each azimuth can be obtained to determine the fracture breakpoint response at each azimuth.

[0058] The velocity error sensitivity and crack breakpoint responsiveness corresponding to each azimuth can be obtained. The corresponding velocity error sensitivity threshold and crack breakpoint responsiveness threshold are pre-set. For each azimuth, if its velocity error sensitivity is greater than or equal to the pre-set velocity error sensitivity threshold or its crack breakpoint responsiveness is greater than or equal to the pre-set crack breakpoint responsiveness threshold, the azimuth is considered to be the dominant azimuth. After determining the dominant azimuth among multiple azimuths, an azimuth weight matrix can be constructed according to the velocity error sensitivity and crack breakpoint responsiveness of the dominant azimuth. The weight corresponding to the dominant azimuth can be determined according to the velocity error sensitivity and crack breakpoint responsiveness of the dominant azimuth. The mapping relationship between the dominant azimuth and the weight corresponding to the dominant azimuth can be constructed by weighted regression or machine learning algorithm. For example, the weighted least squares method can be used to construct the mapping relationship between the velocity error sensitivity and crack breakpoint responsiveness of the dominant azimuth and the weight corresponding to the dominant azimuth. The velocity error sensitivity and fracture breakpoint responsiveness of the dominant azimuth can be used as input parameters, and the weighted least square method can assign different weights to the velocity error sensitivity and fracture breakpoint responsiveness of the dominant azimuth, and output the corresponding weight of the dominant azimuth under the weight addition. For the non-dominant azimuths in the azimuth weight matrix, the weights of these azimuths can be uniformly set to a constant, such as 1.

[0059] S104: updating the initial tomographic velocity model according to the azimuth weight matrix.

[0060] In some embodiments, the following tomographic velocity inversion matrix equations may be established based on the residual moveout at multiple azimuths of the study area, the ray paths at multiple azimuths of the study area, and the azimuth weight matrix:

[0061]

[0062] Among them, W n (n=1,...N) represents the weight of the nth azimuth; L n (n=1,...N) represents the ray path of the nth azimuth angle in the study area; ΔS n (n=1,...N) represents the residual slowness of the nth azimuth in the study area; ΔT n (n=1,...,N) represents the residual time difference of the nth azimuth of the study area; according to the tomographic velocity inversion matrix equation group, the residual slowness of multiple azimuths of the study area is solved; based on the residual slowness of multiple azimuths of the study area, the tomographic velocity update amount of multiple azimuths of the study area is determined; according to the tomographic velocity update amount of multiple azimuths of the study area, the initial tomographic velocity model is updated.

[0063] By constructing a tomographic velocity inversion matrix and updating the tomographic velocity field of the study area according to the solved tomographic velocity update value, a more stable and higher-precision azimuthal grid tomographic velocity inversion is achieved, thereby improving the efficiency of tomographic velocity modeling.

[0064] In mathematical expression, tomographic velocity inversion is to solve a set of linear equations: L·ΔS=ΔT. Where L is a matrix whose elements represent the ray paths in the discrete grids of the study area; ΔS represents the slowness change in each grid; ΔT is the difference between the theoretical travel time and the actual travel time calculated using the initial tomographic velocity model. The ray paths of multiple azimuths, the residual time differences of multiple azimuths, and the azimuth weight matrix can be obtained, and the tomographic velocity inversion matrix equations constrained by the dominant azimuth are established by the following formula:

[0065]

[0066] Where W n (n=1,...N) represents the weight of the nth azimuth angle; Indicates the azimuth angle as α n , the offset is h m The ray path length of the shot detection pair in the discrete grid cell indexed i, j, k, and denote the residual slowness and residual time difference of the discrete grid cells indexed as i, j, and k respectively.

[0067] The L2 norm iteration can be used to solve the tomographic velocity inversion equations, and regularization constraints can be performed to obtain the residual slowness. The above tomographic inversion equations can be simplified to:

[0068] W α L α ΔS=ΔT α

[0069] The system of equations after regularization can be expressed as:

[0070]

[0071] Where λ is the regularization parameter used to balance the data fitting term (i.e. ) and the regularization term (i.e. ).

[0072] Based on the regularized system of equations, construct the objective function:

[0073]

[0074] Taking its derivative we can get:

[0075] ΔS=((W α Lα ) T (W α L α )+λI) -1 (W α L α ) T ΔT α

[0076] The above residual slowness can be solved iteratively using an iterative method, which may be a conjugate gradient method, a least square method, etc., which will not be described in detail here.

[0077] Slowness represents the reciprocal of the velocity at which seismic waves propagate in a medium. Residual slowness represents the reciprocal of the velocity error in the initial tomographic velocity model. The residual slowness obtained by the solution can be taken as the reciprocal to obtain the tomographic velocity update values ​​for multiple azimuths of the study area. Based on the tomographic velocity update values ​​of the tomographic velocity field of the study area, the tomographic velocity field and the initial tomographic velocity model of the study area can be updated.

[0078] In some embodiments, multi-azimuth and omnidirectional pre-stack depth migration calculations can be completed based on the initial tomographic velocity model updated in the above steps and multi-azimuth and omnidirectional seismic data. Evaluate the straightness of each azimuth migration gather, the imaging quality of the depth migration stacked profile, and the rationality of the updated tomographic velocity model. If the tomographic velocity model meets the requirements, the modeling process is terminated and the tomographic velocity model is output. If the requirements are not met, the research area is reselected, and new multi-azimuth and omnidirectional seismic data and a new initial tomographic velocity model are obtained based on the reselected research area, and then the azimuth weight matrix of the dominant azimuth is adjusted, and the process in the above steps is repeated until a tomographic velocity model that meets the requirements is output.

[0079] In some embodiments, the residual moveout at multiple azimuths of the study area may be calculated based on the residual curvature of the depth-migrated imaging gathers at multiple azimuths of the study area.

[0080] The residual time differences at multiple azimuths are calculated by using the residual curvature of the depth-migrated imaging gathers at multiple azimuths in the study area, laying a data foundation for the construction of the tomographic velocity inversion matrix.

[0081] According to the residual curvature at the seed point of the study area at each azimuth, the residual depth curve of the seed point corresponding to the residual curvature can be determined. The seed point is taken as a reflection point, and the residual depth of the seed point can be calculated by the residual depth curve of the seed point. The residual depth refers to the vertical distance from the surface (or a reference surface) to the reflection point, minus the theoretical depth calculated according to the initial tomographic velocity model. This difference reflects the error degree of the initial tomographic velocity model. The residual time difference at each azimuth includes the residual time difference of all seed points in the study area at each azimuth. The residual time difference refers to the difference between the time required for seismic waves to propagate in the underground medium and the time required to propagate in the ideal uniform medium, that is, the difference between the theoretical travel time and the actual travel time calculated using the initial tomographic velocity model. The theoretical travel time of seismic waves in the discrete grid where each seed point is located can be calculated based on prior seismic data. According to the initial velocity model and the residual depth at each seed point, the actual travel time of seismic waves in the discrete grid where each seed point is located can be calculated. The residual time difference of each seed point can be obtained by subtracting the two, and then the residual time difference of the current azimuth can be determined. Repeat the above steps for each azimuth to determine the residual time difference for multiple azimuths in the study area.

[0082] In some embodiments, the tomographic velocity fields at multiple azimuths of the study area may be determined based on the initial tomographic velocity model; and the ray paths at multiple azimuths of the study area may be calculated according to the tomographic velocity fields and inclination fields at multiple azimuths of the study area.

[0083] Based on the initial tomographic velocity model, the tomographic velocity field of multiple azimuths in the study area can be determined. The tomographic velocity field of each azimuth in the study area includes the tomographic velocity of each discrete grid in the study area. The ray path of each azimuth in the study area includes the ray paths of all seed points in the study area at each azimuth. The seed point can be used as a reflection point, and the normal direction of the reflection interface at the seed point can be obtained from the corresponding inclination of each seed point. Starting from the seed point, based on the tomographic velocity corresponding to various sub-points, ray tracing is performed on both sides of the normal at the same angle as the normal angle. The exit points of these two rays on the ground are recorded as the excitation point and the receiving point, respectively, so as to obtain a pair of rays of the reflection point and the corresponding offset distance. Select the ray pairs corresponding to the offset distance in the remaining depth at the same seed point to obtain the ray path in the discrete grid where the seed point is located. When the azimuth is determined, repeat the above steps for each seed point at the current azimuth to obtain the ray path of each seed point, and then determine the ray path corresponding to the current azimuth. Repeat the above steps for each azimuth to determine the ray paths for multiple azimuths of the study area.

[0084] The following introduces a specific scenario example of the embodiment of this specification.

[0085] 1. Arrange a multi-azimuth seismic wave excitation and acquisition observation system in area A, record multi-azimuth seismic wave data, remove noise, correct time delay and perform stacking processing to obtain multi-azimuth seismic data in the time domain. Figure 4 Part (a) shows the theoretical tomographic velocity model of the three-dimensional salt dome model. Figure 4 Part (b) shows the multi-azimuth seismic data based on the theoretical tomographic velocity model.

[0086] 2. Based on geological data, logging data and time domain migration velocity data, an initial tomographic velocity model is established in the depth domain. Figure 5 Part (a) shows the initial tomographic velocity model constructed with local high-velocity bodies. Figure 5 Part (b) shows the multi-azimuth seismic data based on the initial tomographic velocity model.

[0087] 3. Use multi-azimuth seismic data to perform multi-azimuth prestack depth migration calculations on the initial tomographic velocity model containing local high-velocity bodies, and output depth migration imaging gathers at multiple azimuths. Figure 6 , Figure 6 The depth migration imaging gathers of 6 azimuth angles are shown. Figure 6 Parts (a), (b), (c), (d), (e), and (f) represent the depth migration imaging gathers at azimuths of 0°, 30°, 60°, 90°, 120°, and 150°, respectively. It can be seen that the prestack depth migration gathers at different azimuths have different sensitivities to the velocity error of the anomaly.

[0088] 4. Stack the depth migration imaging gathers of 6 azimuths to generate stacked sections. Calculate the inclination field, quality factor field, residual curvature, residual moveout and ray path corresponding to the 6 azimuths. Figure 7 Part (a) shows the residual depth curve corresponding to an azimuth of 0°. Figure 7 Part (b) shows the corresponding residual depth curve when the azimuth angle is 60°. It can be seen that different azimuth angles have different recognition of velocity errors and different residual depth curves.

[0089] 5. Calculate the velocity error sensitivity and crack breakpoint responsiveness corresponding to each of the six azimuths. According to the pre-set velocity error sensitivity threshold and crack breakpoint responsiveness threshold, azimuths 2, 3, 4, and 5 are determined to be dominant azimuths. According to the velocity error sensitivity and crack breakpoint responsiveness of the dominant azimuths, the weights of azimuths 3 and 4 are 10, the weights of azimuths 2 and 5 are 8, and the weights of non-dominant azimuths are uniformly set to 1.

[0090] 6. Establish the tomographic velocity inversion matrix equation group and iteratively solve the tomographic velocity update. Figure 8 Part (a) shows the tomographic velocity field profile of the initial tomographic velocity model after updating based on the tomographic velocity modeling method of the present invention. Figure 8 Part (b) shows the tomographic velocity field profile of the updated initial tomographic velocity model obtained based on traditional multi-azimuth grid tomography.

[0091] Table 1

[0092]

[0093] 7. Complete multi-azimuth prestack depth migration calculation based on the updated tomographic velocity model and multi-azimuth seismic data. Fig. 9 and Table 1 shows Figure 8 The two tomographic velocity fields are shown, compared with the theoretical model ( Figure 3 ) are compared with the RMS amplitude error. Fig. 9 Part (a) shows the theoretical tomographic velocity model, Fig. 9 Part (b) shows the initial tomographic velocity model. Fig. 9 Part (c) shows the tomographic velocity model after an update of the present invention. Fig. 9 Part (d) shows the chromatographic velocity model after one round of updating of the traditional multi-directional chromatography method. It can be seen that the chromatographic velocity modeling method of the present invention can converge the error faster for the updated chromatographic velocity model of the heterogeneous body than the traditional multi-directional chromatography method. Table 1 shows that after the chromatographic velocity model is updated based on the present invention, the velocity matching rate is improved by 24% (the relative error is reduced from 3.65% to 2.77%, the reduction percentage is 0.88%, and the matching rate is improved by: 0.88 / 3.65% = 24%).

[0094] 8. Fig.10 A comparison chart of single-channel velocity curves of heterogeneous bodies extracted by the theoretical chromatography velocity model (line 1), the initial chromatography velocity model (line 2), the chromatography velocity model after an updated round of the present invention (line 3), and the chromatography velocity model after an updated round of the traditional multi-directional chromatography method (line 4) is shown. Fig.11 Part (a) shows the depth migration imaging gathers corresponding to the updated tomographic velocity model based on the tomographic velocity modeling method of the present invention. Fig.11 Part (b) shows the depth migration imaging gathers corresponding to the updated tomographic velocity model based on the traditional multi-azimuth tomography method.

[0095] 9. Fig.12 Part (a) shows the depth migration imaging gathers and vertical spectra corresponding to the tomographic velocity model after the tomographic velocity modeling method of the present invention is updated. Fig.12 Part (b) shows the depth migration imaging gathers and vertical spectra corresponding to the initial velocity model. Fig.12Part (c) shows the depth migration imaging gathers and vertical spectra corresponding to the updated tomographic velocity model based on the traditional multi-azimuth tomography method. It can be seen that the velocity model updated by the tomographic velocity modeling method of the present invention can effectively improve the quality of the depth migration gathers.

[0096] 10. Fig.13 Part (a) shows the depth migration stack section corresponding to the initial tomographic velocity model in area B. Fig.13 Part (b) shows the depth migration stacking section after the initial tomographic velocity model is updated using the tomographic velocity modeling method of the present invention in area B. From the comparison diagram, it can be seen that after the initial tomographic velocity model is updated using the tomographic velocity modeling method of the present invention, the quality of the depth migration data is improved, indicating the feasibility and effectiveness of the actual production of the present invention.

[0097] Based on the above-mentioned tomographic velocity modeling method, this specification also proposes an embodiment of a tomographic velocity modeling device. Fig.14 As shown, the tomographic velocity modeling device may specifically include the following modules:

[0098] The acquisition module 1401 can be used to acquire depth migration imaging gathers of multiple azimuth angles in the research area according to the initial tomographic velocity model;

[0099] A determination module 1402 may be used to determine a dominant azimuth based on depth migration imaging gathers of multiple azimuths in the study area;

[0100] A construction module 1403 may be used to construct an azimuth weight matrix according to the dominant azimuth;

[0101] The updating module 1404 may be configured to update the initial tomographic velocity model according to the azimuth weight matrix.

[0102] In some embodiments, the acquisition module 1401 can be specifically used to establish an initial tomographic velocity model based on geological data, logging data and time domain migration velocity data; based on the seismic data and initial tomographic velocity model of the study area, generate depth migration imaging gathers of multiple azimuths of the study area.

[0103] In some embodiments, the above-mentioned determination module 1402 can be specifically used to superimpose the depth-migrated imaging gathers of multiple azimuths of the study work area to obtain a stacked profile of the study work area; based on the stacked profile of the study work area, determine the velocity error sensitivity of multiple azimuths of the study work area; based on the velocity error sensitivity of multiple azimuths of the study work area, determine the dominant azimuth.

[0104] In some embodiments, the above-mentioned determination module 1402 can also be specifically used to calculate the inclination field of multiple azimuths of the study work area according to the stacked profile; determine the quality factor field of multiple azimuths of the study work area based on the inclination field of multiple azimuths of the study work area; calculate the residual curvature of the depth migration imaging track set of multiple azimuths of the study work area using similarity analysis based on the quality factor field of multiple azimuths of the study work area; determine the velocity error sensitivity of multiple azimuths of the study work area based on the residual curvature of the depth migration imaging track set of multiple azimuths of the study work area.

[0105] In some embodiments, the above-mentioned construction module 1403 can be specifically used to calculate the fracture breakpoint responsiveness of multiple azimuths in the study area based on geological data; determine the dominant azimuth based on the velocity error sensitivity and fracture breakpoint responsiveness of multiple azimuths in the study area; and construct an azimuth weight matrix based on the velocity error sensitivity and fracture breakpoint responsiveness of the dominant azimuth.

[0106] In some embodiments, the updating module 1404 can be specifically used to establish the following tomographic velocity inversion matrix equations according to the residual time difference of multiple azimuths of the study area, the ray paths of multiple azimuths of the study area, and the azimuth weight matrix:

[0107]

[0108] Among them, W n (n=1,...N) represents the weight of the nth azimuth; L n (n=1,...N) represents the ray path of the nth azimuth angle in the study area; ΔS n (n=1,...N) represents the residual slowness of the nth azimuth in the study area; ΔT n (n=1,...,N) represents the residual time difference of the nth azimuth of the study area; according to the tomographic velocity inversion matrix equation group, the residual slowness of multiple azimuths of the study area is solved; based on the residual slowness of multiple azimuths of the study area, the tomographic velocity update amount of multiple azimuths of the study area is determined; according to the tomographic velocity update amount of multiple azimuths of the study area, the initial tomographic velocity model is updated.

[0109] In some embodiments, the updating module 1404 may be specifically used to calculate the residual time difference of multiple azimuths of the study work area according to the residual curvature of the depth-migrated imaging gathers of multiple azimuths of the study work area.

[0110] In some embodiments, the above-mentioned updating module 1404 can also be specifically used to determine the tomographic velocity field of multiple azimuths of the study area based on the initial tomographic velocity model; and calculate the ray paths of multiple azimuths of the study area according to the tomographic velocity field and inclination field of multiple azimuths of the study area.

[0111] It should be noted that the units, devices or modules described in the above embodiments can be implemented by computer chips or entities, or by products with certain functions. For the convenience of description, the above devices are described separately by functions divided into various modules. Of course, when implementing this specification, the functions of each module can be implemented in the same or more software and / or hardware, or the modules that implement the same function can be implemented by a combination of multiple sub-modules or sub-units. The device embodiments described above are only schematic. For example, the division of the units is only a logical function division. There may be other division methods in actual implementation, such as multiple units or components can be combined or integrated into another system, or some features can be ignored or not executed. Another point is that the mutual coupling or direct coupling or communication connection shown or discussed can be through some interfaces, indirect coupling or communication connection of devices or units, which can be electrical, mechanical or other forms.

[0112] As can be seen from the above, based on the tomographic velocity modeling device provided in the embodiments of this specification, the dominant azimuth can be determined based on the depth-migrated imaging gathers of multiple azimuths in the research area, and a weight matrix can be constructed based on the dominant azimuth to further determine the tomographic velocity update amount, thereby achieving more stable and higher-precision tomographic velocity inversion, improving the efficiency of tomographic velocity modeling, and promoting the actual production application of tomographic velocity inversion technology for azimuth grids in complex exploration.

[0113] The embodiments of the present specification also provide a computer device for a tomographic velocity modeling method, comprising a processor and a memory for storing instructions executable by the processor, wherein the processor can perform the following steps according to the instructions when implemented: obtaining depth migration imaging gathers of multiple azimuths of a study area according to an initial tomographic velocity model; determining a dominant azimuth based on the depth migration imaging gathers of multiple azimuths of the study area; constructing an azimuth weight matrix according to the dominant azimuth; and updating the initial tomographic velocity model according to the azimuth weight matrix.

[0114] In order to complete the above instructions more accurately, refer to Fig.15 As shown, the embodiment of this specification also provides a computer device based on the tomographic velocity modeling method, wherein the computer device includes a network communication port 1501, a processor 1502 and a memory 1503, and the above structures are connected by internal cables so that each structure can perform specific data interaction.

[0115] In the embodiments of this specification, the network communication port 1501 can be a virtual port that is bound to different communication protocols to send or receive different data. For example, the network communication port 1501 can be a port responsible for web data communication, or it can also be a port responsible for FTP data communication, or it can also be a port responsible for email data communication. In addition, the network communication port 1501 can be a physical communication interface or communication chip. For example, the network communication port 1501 can be a wireless mobile network communication chip. For example, the network communication port 1501 can be a wireless mobile network communication chip that supports protocols such as GSM and CDMA. The network communication port 1501 can also be a Wifi chip or a Bluetooth chip, etc.

[0116] The processor 1502 can be specifically used to obtain depth migration imaging gathers of multiple azimuths of the study area based on the initial tomographic velocity model; determine the dominant azimuth based on the depth migration imaging gathers of multiple azimuths of the study area; construct an azimuth weight matrix based on the dominant azimuth; and update the initial tomographic velocity model based on the azimuth weight matrix.

[0117] In the embodiments of this specification, the processor 1502 may be implemented in any appropriate manner. For example, the processor 1502 may take the form of a microprocessor or processor and a computer-readable medium storing a computer-readable program code (such as software or firmware) executable by the (micro)processor, a logic gate, a switch, an application-specific integrated circuit (ASIC), a programmable logic controller, and an embedded microcontroller, etc. This specification does not limit this.

[0118] The memory 1503 is specifically used to store corresponding instruction programs.

[0119] In the embodiments of this specification, the memory 1503 includes volatile memory and non-volatile memory. In digital systems, anything that can store binary data can be a memory; in integrated circuits, a circuit with storage function without a physical form is also called a memory, such as RAM, FIFO, etc.; in computer systems, a storage device with a physical form is also called a memory, such as a memory stick, TF card, etc.

[0120] Those skilled in the art will appreciate that embodiments of the present invention may be provided as methods, systems, or computer program products. Therefore, the present invention may take the form of a complete hardware embodiment, a complete software embodiment, or an embodiment combining software and hardware. Moreover, the present invention may take the form of a computer program product implemented on one or more computer-usable storage media (including but not limited to disk storage, CD-ROM, optical storage, etc.) containing computer-usable program code.

[0121] The present invention is described with reference to flowcharts and / or block diagrams of methods, devices (systems), and computer program products according to embodiments of the present invention. It should be understood that each process and / or block in the flowchart and / or block diagram, as well as the combination of processes and / or blocks in the flowchart and / or block diagram, can be implemented by computer program instructions. These computer program instructions can be provided to a processor of a general-purpose computer, a special-purpose computer, an embedded processor, or other programmable data processing device to produce a machine, so that the instructions executed by the processor of the computer or other programmable data processing device generate instructions for implementing the processes in the flowchart and / or block diagram. Figure 1 A process or multiple processes and / or boxes Figure 1 A device that provides the functions specified in a block or multiple blocks.

[0122] These computer program instructions may also be stored in a computer-readable memory capable of directing a computer or other programmable data processing device to operate in a specific manner, so that the instructions stored in the computer-readable memory produce an article of manufacture comprising an instruction device, which implements the process Figure 1 A process or multiple processes and / or boxes Figure 1 A function specified in one or more boxes.

[0123] These computer program instructions can also be loaded onto a computer or other programmable data processing device so that a series of operating steps are executed on the computer or other programmable device to produce a computer-implemented process, thereby providing instructions for implementing the process. Figure 1 A process or multiple processes and / or boxes Figure 1 The steps for the functions specified in one or more boxes.

[0124] The specific embodiments described above further illustrate the objectives, technical solutions and beneficial effects of the present invention in detail. It should be understood that the above description is only a specific embodiment of the present invention and is not intended to limit the scope of protection of the present invention. Any modifications, equivalent substitutions, improvements, etc. made within the spirit and principles of the present invention should be included in the scope of protection of the present invention.

Claims

1. A tomographic velocity modeling method, characterized in that: include: According to the initial tomographic velocity model, depth migration imaging gathers of multiple azimuths in the study area are obtained; Determining a dominant azimuth based on depth migration imaging gathers at multiple azimuths in the research area; According to the dominant azimuth, construct an azimuth weight matrix; The initial tomographic velocity model is updated according to the azimuth weight matrix.

2. The method according to claim 1, characterized in that: The method of obtaining depth migration imaging gathers of multiple azimuth angles in the research area according to the initial tomographic velocity model includes: An initial tomographic velocity model is established based on geological data, well logging data and time domain migration velocity data; Based on the seismic data and initial tomographic velocity model of the study area, depth migration imaging gathers of multiple azimuths in the study area are generated.

3. The method according to claim 1, characterized in that The determining of the dominant azimuth angle based on the depth migration imaging gathers of multiple azimuth angles in the research area comprises: Superimposing depth migration imaging gathers of multiple azimuths of the study area to obtain a stacked profile of the study area; Based on the stacked profiles of the study area, determining the velocity error sensitivity of multiple azimuths of the study area; Based on the velocity error sensitivity of multiple azimuths in the research work area, the dominant azimuth is determined.

4. The method according to claim 3, characterized in that The determining of the velocity error sensitivity of multiple azimuth angles of the study area based on the stacked profile of the study area includes: According to the stacked profile, the dip fields at multiple azimuths of the study area are calculated; Determining the quality factor fields of the multiple azimuth angles of the study work area based on the inclination fields of the multiple azimuth angles of the study work area; Based on the quality factor fields of multiple azimuth angles in the study area, similarity analysis is used to calculate the residual curvature of the depth migration imaging gathers at multiple azimuth angles in the study area; Based on the residual curvature of the depth-migrated imaging gathers at multiple azimuth angles in the study area, the velocity error sensitivity at multiple azimuth angles in the study area is determined.

5. The method according to claim 1, characterized in that: The step of constructing an azimuth weight matrix according to the dominant azimuth includes: Based on geological data, the fracture breakpoint responsiveness at multiple azimuths in the study area is calculated; Based on the velocity error sensitivity and fracture breakpoint responsiveness at multiple azimuths in the study area, the dominant azimuth is determined; According to the velocity error sensitivity and fracture breakpoint responsiveness of the dominant azimuth, an azimuth weight matrix is ​​constructed.

6. The method according to claim 5, characterized in that The updating of the initial tomographic velocity model according to the azimuth weight matrix comprises: According to the residual time difference of multiple azimuths in the study area, the ray paths of multiple azimuths in the study area and the azimuth weight matrix, the following tomographic velocity inversion matrix equations are established: Among them, W n (n=1,...N) represents the weight of the nth azimuth; L n (n=1,...N) represents the ray path of the nth azimuth angle in the study area; ΔS n (n=1,...N) represents the residual slowness of the nth azimuth in the study area; ΔT n (n=1,...,N) represents the residual time difference of the nth azimuth in the study area; According to the tomographic velocity inversion matrix equation group, the residual slowness of multiple azimuths in the study area is solved; Determining tomographic velocity update amounts for multiple azimuths of the study area based on the residual slownesses for multiple azimuths of the study area; The initial tomographic velocity model is updated according to the tomographic velocity update amounts at multiple azimuth angles in the research work area.

7. The method according to claim 6, characterized in that The method further comprises: According to the residual curvature of the depth-migrated imaging gathers at multiple azimuths in the study area, the residual moveout at multiple azimuths in the study area is calculated.

8. The method according to claim 6, characterized in that The method further comprises: Based on the initial tomographic velocity model, determining the tomographic velocity fields at multiple azimuth angles in the research area; According to the tomographic velocity field and the inclination field at multiple azimuths of the research area, the ray paths at multiple azimuths of the research area are calculated.

9. A tomographic velocity modeling device, characterized in that: The device comprises: An acquisition module is used to acquire depth migration imaging gathers of multiple azimuths in the study area according to the initial tomographic velocity model; A determination module, used for determining a dominant azimuth based on depth migration imaging gathers of multiple azimuths in the research area; A construction module, used for constructing an azimuth weight matrix according to the dominant azimuth; An updating module is used to update the initial tomographic velocity model according to the azimuth weight matrix.

10. A computer device comprising a memory, a processor, and a computer program stored in the memory and executable on the processor, characterized in that: When the processor executes the computer program, the method according to any one of claims 1 to 8 is implemented.

Images