An adaptive travel-time two-way fusion interpolation method and device and a storage medium
Patent Information
- Application Number
- CN202610951742.5
- Authority / Receiving Office
- CN · China
- Patent Type
- Applications(China)
- Current Assignee / Owner
- Filing Date
- 2026-06-29
- Publication Date
- 2026-11-03
AI Technical Summary
Existing technologies struggle to fill the missing travel time caused by high-velocity anomalies when processing seismic migration imaging, leading to blurred images and reduced resolution, especially below high-velocity anomaly layers where it is difficult to improve longitudinal resolution.
An adaptive bidirectional fusion interpolation method for travel time is adopted. Travel time is calculated by combining ray tracing with the Gaussian beam method. Multi-directional interpolation function fitting is used to select the optimal radial and tangential interpolation functions to fill in the missing values. The interpolation is then converted from Cartesian coordinates to polar coordinates. Weighted fusion is performed in combination with seismic wave propagation characteristics. Finally, the PSF is calculated analytically and deconvolution is performed.
It achieves efficient and high-fidelity least-squares migration imaging, which significantly improves imaging performance, especially the vertical resolution of the slice below high-velocity anomalies, reduces imaging blur, and enhances imaging clarity and resolution.
Smart Images
Figure CN123029612A_ABST
Abstract
Description
Technical Field
[0001] This invention relates to the field of seismic imaging technology, and more specifically, to an adaptive travel-time bidirectional fusion interpolation method, apparatus, and storage medium. Background Technology
[0002] Seismic migration imaging is a crucial tool in oil and gas exploration. It processes seismic data to generate seismic profiles, providing a foundation for subsequent seismic data interpretation and oil and gas reservoir prediction. However, in the theoretical solution of traditional migration imaging, migration algorithms typically use the transpose of the forward operator as the inverse operation, leading to blurred imaging results and a decrease in resolution. As seismic fine imaging technology gradually develops towards the least-squares migration method, seismic wave inversion imaging, as an extension of traditional migration imaging within the least-squares framework, constructs the Hessian matrix by analytically calculating the point spread function (PSF). This method can deblur the results of traditional migration imaging, thereby significantly improving the resolution of layer imaging and obtaining clearer imaging effects. Analytical characterization of PSF requires calculating the travel time of seismic waves from the excitation point to the subsurface and back to the receiver. Calculating travel time based on ray tracing combined with the Gaussian beam method is accurate and efficient. However, when high-velocity anomalies or high-velocity anomaly layers are present underground, the propagation direction based on ray tracing, following Snell's law, will exceed the critical angle, making it impossible to track the propagation of seismic waves across the entire area. This results in some exploration areas having missing values in the travel time calculated based on ray tracing. Furthermore, PSF based on analytical characterization with missing travel times suffers from poor illumination when applying deconvolution to the traditional migration imaging results for deblurring. Since it is difficult to obtain travel time below high-velocity anomalies or high-velocity anomaly layers, it is difficult to improve the vertical resolution of unclear layers.
[0003] How to efficiently and accurately fill in the gaps in a travel itinerary is an urgent problem to be solved. Summary of the Invention
[0004] The purpose of this invention is to provide an adaptive travel-time bidirectional fusion interpolation method, apparatus, and storage medium that can achieve efficient, high-fidelity, and non-iterative least-squares offset imaging.
[0005] This invention provides an adaptive bidirectional fusion interpolation method for travel time, comprising the following steps: S1: Based on the velocity model of the target work area, the travel time of different excitation points in the entire work area, including the missing values, is calculated using the method of ray tracing combined with Gaussian beam. S2: Based on the velocity model of the target work area and the travel time results, perform tangential and radial multi-directional interpolation function fitting based on different excitation points, and select the optimal radial and tangential interpolation functions; S3: The travel time results containing missing values are split according to different excitation points to obtain travel time data. The travel time of each excitation point is calculated based on the critical angle of each excitation point. The travel time is filled into the missing points to obtain the travel time data in the Cartesian coordinate system after tangential interpolation boundary filling. The travel time data after tangential interpolation boundary filling is converted from the Cartesian coordinate system to the polar coordinate system to obtain the travel time data in the polar coordinate system. S4: Based on the travel time data in the Cartesian and polar coordinate systems, calculate the travel distance and travel angle corresponding to each missing value, as well as the boundaries of radial and tangential interpolation; S5: Based on the boundaries of the radial and tangential interpolation, using the optimal radial and tangential interpolation functions and the travel distance and travel angle corresponding to each missing value, perform comprehensive interpolation to obtain the interpolation result. Fill the missing space with the interpolation result to obtain the global travel time. S6: Calculate the PSF using the global travel time analysis, and perform deconvolution processing on the offset imaging results based on the PSF to obtain the least squares offset results.
[0006] Furthermore, the calculation formula for the above multi-directional interpolation function fitting is as follows:
[0007] in, The sample value is the approximate travel time value; d is the ray path length from the excitation point to the target point; Indicates the root mean square velocity of the target point; This represents the multi-directional optimal interpolation function, i.e., different radial distances. and tangential angle The optimal interpolation function is composed of a basic function and a polynomial function; The computational notation for minimizing the function; It is the path distance r and the tangential angle. When traveling to a point; The interpolation error is represented by n; n is the number of global imaging points calculated during travel time. This represents the reference value for the true travel time of imaging point i; Based on different radial distances and tangential angle Calculated travel time.
[0008] Furthermore, the formulas for calculating the travel distance and travel angle mentioned above are as follows:
[0009] in, and These refer to travel distance and travel angle, respectively. and These are the coordinates corresponding to the missing values in the Cartesian coordinate system. The formulas for calculating the boundaries of the above radial and tangential interpolation are as follows:
[0010] in, and These are the boundaries for radial interpolation and tangential interpolation, respectively. When representing the nearest radial distance, This indicates the shortest travel time with the tangential value.
[0011] Furthermore, the calculation formula for the above comprehensive interpolation is as follows:
[0012] in, Indicates the travel time in the radial direction; Indicates the magnitude of travel time in the tangential direction; Represents the radial interpolation function that is optimal for the objective. Represents the optimal tangential interpolation function for the objective; Indicates the interpolation interval coefficients; These represent the nearest radial distance and the nearest tangential distance, respectively. The result is a comprehensive interpolation of the target.
[0013] The present invention also provides an adaptive travel-time bidirectional fusion interpolation device, comprising the following modules: The data interpolation preparation module is configured to: calculate the travel time of different excitation points in the entire work area, including missing values, based on the velocity model of the target work area using ray tracing combined with Gaussian beam method; The fitting function selection module is configured to: based on the velocity model of the target work area and the travel time results, perform tangential and radial multi-directional interpolation function fitting based on different excitation points, and select the optimal radial and tangential interpolation functions; The first data interpolation preprocessing module is configured to: split the travel time results containing missing values into different excitation points to obtain travel time data; calculate the travel time of each excitation point based on the critical angle of each excitation point; fill the missing travel time into the missing points to obtain travel time data in the Cartesian coordinate system after tangential interpolation boundary filling; and convert the travel time data after tangential interpolation boundary filling from the Cartesian coordinate system to the polar coordinate system to obtain travel time data in the polar coordinate system. The second data interpolation preprocessing module is configured to: calculate the travel distance and travel angle corresponding to each missing value, as well as the boundaries of radial interpolation and tangential interpolation, based on the travel time data in the Cartesian coordinate system and the polar coordinate system. The multi-directional interpolation fusion module is configured to: based on the boundaries of the radial and tangential interpolation, use the optimal radial and tangential interpolation functions and the travel distance and travel angle corresponding to each missing value to perform comprehensive interpolation to obtain the interpolation result, fill the missing space with the interpolation result, and obtain the global travel time; The travel time application module is configured to: use the global travel time analysis to calculate the PSF, and perform deconvolution processing on the migration imaging results based on the PSF to obtain the least squares migration result.
[0014] Furthermore, the calculation formula for the above multi-directional interpolation function fitting is as follows:
[0015] in, The sample value is the approximate travel time value; d is the ray path length from the excitation point to the target point; Indicates the root mean square velocity of the target point; This represents the multi-directional optimal interpolation function, i.e., different radial distances. and tangential angle The optimal interpolation function is composed of a basic function and a polynomial function; The computational notation for minimizing the function; It is the path distance r and the tangential angle. When traveling to a point; The interpolation error is represented by n; n is the number of global imaging points calculated during travel time. This represents the reference value for the true travel time of imaging point i; Based on different radial distances and tangential angle Calculated travel time.
[0016] Furthermore, the formulas for calculating the travel distance and travel angle mentioned above are as follows:
[0017] in, and These refer to travel distance and travel angle, respectively. and These are the coordinates corresponding to the missing values in the Cartesian coordinate system. The formulas for calculating the boundaries of the radial interpolation and tangential interpolation are as follows:
[0018] in, and These are the boundaries for radial interpolation and tangential interpolation, respectively. When representing the nearest radial distance, This indicates the shortest travel time with the tangential value.
[0019] Furthermore, the calculation formula for the above comprehensive interpolation is as follows:
[0020] in, Indicates the travel time in the radial direction; Indicates the magnitude of travel time in the tangential direction; Represents the radial interpolation function that is optimal for the objective. Represents the optimal tangential interpolation function for the objective; Indicates the interpolation interval coefficients; These represent the nearest radial distance and the nearest tangential distance, respectively. The result is a comprehensive interpolation of the target.
[0021] The present invention also provides a computer-readable storage medium having a computer program stored thereon, which, when executed by a processor, implements the steps of the above-described adaptive travel-time bidirectional fusion interpolation method.
[0022] The present invention also provides a computer program product, including a computer program that, when executed by a processor, implements the steps of the above-described adaptive travel-time bidirectional fusion interpolation method.
[0023] Implementing the adaptive travel-time bidirectional fusion interpolation method, apparatus, and storage medium provided by this invention has the following beneficial effects: This invention addresses the shortcomings of existing technologies by further suppressing the shielding effect of high-velocity anomalies and improving the longitudinal resolution of shielded layers. It aims to fill in missing travel times. For efficient and accurate filling of missing values, direct interpolation using known travel times is a better solution. Specifically, based on ray tracing and Gaussian beam methods, the seismic wave travel times at different excitation points within the work area are calculated (which may contain missing values due to the influence of high-velocity anomalies). A two-way fitting analysis is performed based on the seismic velocity model of the work area, and an appropriate interpolation function is selected for subsequent interpolation calculations. The calculated travel time results are converted from Cartesian coordinates to polar coordinates, and regions with zero travel time are marked. The travel times of different excitation points under critical angle excitation are directly calculated to fill the missing values in the polar coordinate system boundary region. The travel radius and travel angle of each excitation point within the blank region are calculated. Combining tangential and radial two-way adaptive interpolation functions, interpolation is performed in both angular and radial directions, and the closest non-null value is used as the distance-weighted interpolation result to fill the missing values.
[0024] Therefore, the interpolation method of this invention is adaptive, selecting the relatively optimal combination of interpolation functions based on the velocity model of the target work area and the interpolation direction. During the interpolation process, this invention achieves interpolation in polar coordinates through coordinate system transformation, which better reflects the propagation characteristics of seismic waves. By interpolating in multiple directions and weighted fusion of the results from different directions, a more reasonable interpolation result is obtained. When using the interpolation method of this invention, the entire process can be directly processed based on the target work area and the calculated true values, resulting in a smaller computational load. By applying the interpolation method of this invention, efficient, high-fidelity, and non-iterative least-squares migration imaging can be achieved throughout the entire process based on ray tracing and analytical calculation of the point spread function (PSF).
[0025] This invention employs a bidirectional weighted fusion method, combining travel time data with missing values in the Cartesian coordinate system, a seismic zone velocity model, seismic wave propagation characteristics, and travel time calculation methods to comprehensively calculate missing values, strictly adhering to the physical characteristics and kinematic theory of seismic waves. Based on the interpolated travel time data, it calculates the seismic imaging blur function and performs deblurring on the seismic imaging effect, thereby significantly improving the imaging quality. Furthermore, it calculates the global travel time using ray tracing and Gaussian beam methods, combining the seismic zone velocity model and seismic wave propagation characteristics, strictly adhering to physical characteristics and kinematic theory, and performs efficient and reasonable adaptive bidirectional interpolation to accurately fill travel time gaps, thus providing better support for subsequent PSF analytical characterization and migration imaging deblurring. This invention has promising application prospects in solving the interpolation processing of missing travel time values and also broadens the applicability of ray tracing combined with Gaussian beam methods for analytical PSF calculation and imaging correction. Attached Figure Description
[0026] The present invention will be further described below with reference to the accompanying drawings and embodiments. In the accompanying drawings: Figure 1 This is a flowchart of the adaptive travel-time bidirectional fusion interpolation method provided by the present invention; Figure 2 This is a schematic diagram illustrating the operation steps of the adaptive travel-time bidirectional fusion interpolation method provided by the present invention; Figure 3 This is a schematic diagram of the target detection zone velocity model provided by the present invention; Figure 4 This is a schematic diagram of the imaging results based on Gaussian beam shift provided by the present invention; Figure 5 This is a schematic diagram of partial null-valued travel time calculated based on ray tracing combined with Gaussian beam method in Cartesian coordinate system (left) and polar coordinate system (right) provided by the present invention; Figure 6This is a schematic diagram of the partial global travel time calculated based on the adaptive bidirectional fusion interpolation method provided by the present invention in the Cartesian coordinate system (left) and polar coordinate system (right); Figure 7 This is a schematic diagram of the PSF results (left) and wavenumber spectrum (right) based on pre-interpolation analytical calculation provided by the present invention; Figure 8 This is a schematic diagram of the PSF result (left) and wavenumber spectrum (right) based on analytical calculation after interpolation provided by the present invention; Figure 9 These are schematic diagrams of the correction of offset imaging results by PSF based on analytical calculation before interpolation (top) and the correction of offset imaging results by PSF based on analytical calculation after interpolation (bottom) provided by the present invention. Detailed Implementation
[0027] To provide a clearer understanding of the technical features, objectives, and effects of the present invention, specific embodiments of the present invention will now be described in detail with reference to the accompanying drawings.
[0028] Figure 1 A schematic diagram of the adaptive bidirectional travel-time fusion interpolation method of this embodiment is shown. In this embodiment, the adaptive bidirectional travel-time fusion interpolation method includes the following steps: S1: Based on the velocity model of the target work area, the travel time of different excitation points in the entire work area, including the missing values, is calculated using the method of ray tracing combined with Gaussian beam. As an exemplary embodiment, in step S1, based on the velocity model of the target work area, the travel time of different excitation points in the whole domain is calculated by using ray tracing combined with Gaussian beam method; the calculation results are traversed, and the positions of the missing values caused by high-speed anomaly shielding are marked to obtain a whole domain travel time dataset in Cartesian coordinate system with missing value position markings. S2: Based on the velocity model of the target work area and the travel time results, perform tangential and radial multi-directional interpolation function fitting based on different excitation points, and select the optimal radial and tangential interpolation functions; In one exemplary embodiment, the calculation formula for the multi-directional interpolation function fitting is:
[0029] in, The sample value is the approximate travel time value; d is the ray path length from the excitation point to the target point; Indicates the root mean square velocity of the target point; This represents the multi-directional optimal interpolation function, i.e., different radial distances. and tangential angle The optimal interpolation function is composed of a basic function and a polynomial function; The computational notation for minimizing the function; It is the path distance r and the tangential angle. When traveling to a point; The interpolation error is represented by n; n is the number of global imaging points calculated during travel time. This represents the reference value for the true travel time of imaging point i; Based on different radial distances and tangential angle Calculated travel time; As an exemplary embodiment, in step S2, tangential and radial multi-directional interpolation functions are fitted based on different excitation points of the work area velocity model. The travel sample value of the target point is directly calculated by using a fixed ray and root mean square velocity. Then, the interpolation function is fitted by changing the radial distance and tangential angle, and the interpolation function with the smallest error is selected as the corresponding interpolation function reserve. S3: The travel time results containing missing values are split according to different excitation points to obtain travel time data. The travel time of each excitation point is calculated based on the critical angle of each excitation point. The travel time is filled into the missing points to obtain the travel time data in the Cartesian coordinate system after tangential interpolation boundary filling. The travel time data after tangential interpolation boundary filling is converted from the Cartesian coordinate system to the polar coordinate system to obtain the travel time data in the polar coordinate system. As an exemplary embodiment, in step S3, the travel time results containing missing values are split according to different excitation points, and then the travel time of each excitation point is calculated according to its critical angle to fill the missing value. The split travel time data is then subjected to coordinate transformation, and converted from the Cartesian coordinate system to the polar coordinate system using the coordinate system transformation formula, and the position of the missing value is recorded. S4: Based on the travel time data in the Cartesian and polar coordinate systems, calculate the travel distance and travel angle corresponding to each missing value, as well as the boundaries of radial and tangential interpolation; In one exemplary embodiment, the formulas for calculating the travel distance and travel angle are as follows:
[0030] in, and These refer to travel distance and travel angle, respectively. and These are the coordinates corresponding to the missing values in the Cartesian coordinate system. In one exemplary embodiment, the formula for calculating the boundaries of the radial interpolation and the tangential interpolation is:
[0031] in, and These are the boundaries for radial interpolation and tangential interpolation, respectively. When representing the nearest radial distance, This indicates the shortest tangential distance; S5: Based on the boundaries of the radial and tangential interpolation, using the optimal radial and tangential interpolation functions and the travel distance and travel angle corresponding to each missing value, perform comprehensive interpolation to obtain the interpolation result. Fill the missing space with the interpolation result to obtain the global travel time. In one exemplary embodiment, the formula for calculating the comprehensive interpolation is:
[0032] in, Indicates the travel time in the radial direction; Indicates the magnitude of travel time in the tangential direction; Represents the radial interpolation function that is optimal for the objective. Represents the optimal tangential interpolation function for the objective; Indicates the interpolation interval coefficients; These represent the nearest radial distance and the nearest tangential distance, respectively. The result is a comprehensive interpolation of the target.
[0033] As an exemplary embodiment, in step S5, the travel time at the multi-directional gap is comprehensively calculated using boundary values and the selected interpolation function. Then, the weights of the multi-directional interpolation results are determined based on their interpolation distances. The comprehensive interpolation, which takes into account the velocity model of the work area, the seismic wave propagation characteristics, physical characteristics, and kinematic theory, is obtained through weighted calculation. The interpolation results are then filled into the gaps. S6: Calculate the PSF using the global travel time analysis, and perform deconvolution processing on the offset imaging results based on the PSF to obtain the least squares offset results.
[0034] This embodiment provides an adaptive travel-time bidirectional fusion interpolation device, comprising the following modules: The data interpolation preparation module is configured to: calculate the travel time of different excitation points in the entire work area, including missing values, based on the velocity model of the target work area using ray tracing combined with Gaussian beam method; The fitting function selection module is configured to: based on the velocity model of the target work area and the travel time results, perform tangential and radial multi-directional interpolation function fitting based on different excitation points, and select the optimal radial and tangential interpolation functions; The first data interpolation preprocessing module is configured to: split the travel time results containing missing values into different excitation points to obtain travel time data; calculate the travel time of each excitation point based on the critical angle of each excitation point; fill the missing travel time into the missing points to obtain travel time data in the Cartesian coordinate system after tangential interpolation boundary filling; and convert the travel time data after tangential interpolation boundary filling from the Cartesian coordinate system to the polar coordinate system to obtain travel time data in the polar coordinate system. The second data interpolation preprocessing module is configured to: calculate the travel distance and travel angle corresponding to each missing value, as well as the boundaries of radial interpolation and tangential interpolation, based on the travel time data in the Cartesian coordinate system and the polar coordinate system. The multi-directional interpolation fusion module is configured to: based on the boundaries of the radial and tangential interpolation, use the optimal radial and tangential interpolation functions and the travel distance and travel angle corresponding to each missing value to perform comprehensive interpolation to obtain the interpolation result, fill the missing space with the interpolation result, and obtain the global travel time; The travel time application module is configured to: use the global travel time analysis to calculate the PSF, and perform deconvolution processing on the migration imaging results based on the PSF to obtain the least squares migration result.
[0035] In one exemplary embodiment, the calculation formula for the above-mentioned multi-directional interpolation function fitting is as follows:
[0036] in, The sample value is the approximate travel time value; d is the ray path length from the excitation point to the target point; Indicates the root mean square velocity of the target point; This represents the multi-directional optimal interpolation function, i.e., different radial distances. and tangential angle The optimal interpolation function is composed of a basic function and a polynomial function; The computational notation for minimizing the function; It is the path distance r and the tangential angle. When traveling to a point; The interpolation error is represented by n; n is the number of global imaging points calculated during travel time. This represents the reference value for the true travel time of imaging point i; Based on different radial distances and tangential angle Calculated travel time.
[0037] In one exemplary embodiment, the formulas for calculating the travel distance and travel angle are as follows:
[0038] in, and These refer to travel distance and travel angle, respectively. and These are the coordinates corresponding to the missing values in the Cartesian coordinate system. The formulas for calculating the boundaries of the radial interpolation and tangential interpolation are as follows:
[0039] in, and These are the boundaries for radial interpolation and tangential interpolation, respectively. When representing the nearest radial distance, This indicates the shortest travel time with the tangential value.
[0040] In one exemplary embodiment, the calculation formula for the above-mentioned comprehensive interpolation is as follows:
[0041] in, Indicates the travel time in the radial direction; Indicates the magnitude of travel time in the tangential direction; Represents the radial interpolation function that is optimal for the objective. Represents the optimal tangential interpolation function for the objective; Indicates the interpolation interval coefficients; These represent the nearest radial distance and the nearest tangential distance, respectively. The result is a comprehensive interpolation of the target.
[0042] In some embodiments, the adaptive travel-time bidirectional fusion interpolation method described above can also be implemented in the following ways.
[0043] like Figure 2 As shown, in this embodiment, the adaptive travel-time bidirectional fusion interpolation method includes the following steps: Step 1: Obtain the velocity model of the target work area, and use the method of ray tracing combined with Gaussian beam to calculate the travel time of different excitation points in the entire work area. Due to the existence of an abnormally high-speed layer, there are gaps in the calculation results of the travel time of the entire area.
[0044] Step 2: Based on the velocity model of the work area, perform multi-directional interpolation function fitting for tangential and radial directions at different excitation points. Directly select fixed rays and root mean square velocities to calculate the travel sample value of the target point. Then, by changing the radial distance and tangential angle, fit the interpolation function and select the interpolation function with the smallest error as the reserve interpolation function for the corresponding direction.
[0045]
[0046] In the formula: Based on different radial distances and tangential angle Calculated travel time; This represents interpolation error; This represents the reference value for actual travel time; n is the number of travel times calculated for all imaging points in the global area. denoted as the multi-directional optimal interpolation function (where the interpolation function consists of a basic function and a polynomial function); d is the ray path length from the excitation point to the target point; This represents the root mean square velocity of the target point.
[0047] Step 3: Split the travel time results containing missing values according to different excitation points, and then calculate the travel time for each excitation point according to its critical angle to fill the missing value. Perform coordinate transformation on the split travel time matrix, convert it from Cartesian coordinate system to polar coordinate system using coordinate system transformation formula, and record the position of the missing value.
[0048] Step 4: After converting the travel time of the target work area into two coordinate systems, the travel time data in the Cartesian coordinate system is traversed to mark null values and calculate their travel distances. and travel perspective :
[0049] In the corresponding polar coordinate system, find the travel time with the nearest radial value as the boundary for radial interpolation and the travel time with the nearest tangential value as the boundary for tangential interpolation.
[0050]
[0051] Step 5: Calculate the travel time at the multi-directional gaps by combining boundary values and the selected interpolation function. Then, determine the weights of the multi-directional interpolation results based on their interpolation distances. Through weighted calculation, obtain the comprehensive interpolation that considers the velocity model of the work area, seismic wave propagation characteristics, physical characteristics, and kinematic theory. Fill the gaps with the interpolation results.
[0052]
[0053] In the formula: Indicates the travel time in the radial direction; Indicates the magnitude of travel time in the tangential direction; Represents the target optimal radial interpolation function. Represents the target optimal tangential interpolation function; Indicates the interpolation interval coefficients; These represent the nearest radial distance and the nearest tangential distance, respectively. Target integrated interpolation results.
[0054] Step 6: Obtain the global travel time by continuously filling in the travel time gaps in the Cartesian coordinate system, calculate the PSF using travel time analysis, and perform deconvolution processing on the migration imaging results based on the PSF to obtain the least squares migration result.
[0055] In some embodiments, the adaptive travel-time bidirectional fusion interpolation device described above can also be implemented in the following ways.
[0056] In this embodiment, the adaptive travel-time bidirectional fusion interpolation device includes the following modules: Data interpolation preparation module: used to calculate the travel time of different excitation points across the entire work area (there are gaps).
[0057] Fitting function selection module: Used to select a multi-directional interpolation function that fits the speed model of the target work area.
[0058] Data interpolation preprocessing module: used to fill in the critical angle ray travel and prepare for subsequent interpolation boundaries.
[0059] Multi-directional interpolation fusion module: used to calculate missing values by combining interpolation functions and interpolation distances in different coordinate systems.
[0060] Travel Time Application Module: Calculates PSF using global travel time analysis and corrects migration imaging results based on PSF.
[0061] In some embodiments, the adaptive travel-time bidirectional fusion interpolation method described above can also be implemented in the following ways.
[0062] In this embodiment, the adaptive travel-time bidirectional fusion interpolation method includes the following steps: The velocity model is imported; the simulated work area is a two-dimensional Hess model with a depth of 7000 meters and a length of 18085 meters. Figure 3 As shown, Gaussian beam migration imaging was first performed on the exploration area, and the results are as follows. Figure 4 As shown, a giant high-speed anomaly exists in the model, which causes the layers below the anomaly to be shielded. In this simulation, the shot point moves 20 meters each time for firing. Ray tracing calculations are performed for each shot point at a selected launch angle. During the firing process of a shot, the travel time is calculated using the Gaussian beam method. The travel times of all launch angles are superimposed to obtain the one-way global travel time of the corresponding shot point. Due to the frequency blocking effect of the giant high-speed anomaly, some areas lack travel time calculations, resulting in gaps.
[0063] Before performing interpolation calculations, the corresponding shot point location is selected, and the depth domain layer velocity model is transformed into the depth domain root mean square velocity model based on equation (1). Then, the approximate travel time sample value is calculated using equation (2). The calculated sample value is used to fill the entire exploration area of the corresponding shot point.
[0064] (1) (2) Based on this method, the travel time at all locations in the exploration area from 0 to 6000 meters, from 6000 meters to 12000 meters, and from 12000 meters to 18000 meters is calculated at 20-meter intervals of shot points, and a two-way interpolation function is fitted based on equation (3).
[0065] (3) Using the arithmetic mean of all values at corresponding locations, the interpolation function with the coefficient of determination closest to 1 was selected. Through experimental investigation, all results were obtained (Table 1). Based on the results in Table 1, it is clear that for the corresponding work area within the range of 0 to 18000 meters, a linear function was chosen for radial interpolation, and a [missing information - likely a specific interpolation function] was chosen for tangential interpolation. The exponential function provides the best fit. Therefore, all subsequent interpolation functions in this method adopt the above strategy.
[0066] Table 1: Average Fit Determination Coefficients for Multi-Function Multi-Source Excitation
[0067] After determining the interpolation function most suitable for the work area, this method begins to fill the missing travel time with interpolation. First, the travel time of the shot point firing to both sides at the critical angle is calculated in the whole domain to fill the missing critical angle. This operation can effectively fill the boundary of tangential interpolation. Then, the entire travel time of a certain shot point is traversed, and the travel radius and travel angle of each non-zero point relative to the shot point position are calculated using Equation (4).
[0068] (4) A new array is constructed to store non-zero travel times and positional relationships. In order to reduce errors, the selected travel radius interval is 100 meters and the travel angle interval is 1°, thus establishing a polar coordinate coefficient group of size (180, 180). After coordinate system transformation, the travel times of the Cartesian coordinate system with gaps are re-traversed to calculate the travel radius and travel time angle at the gap. Then, radial and tangential bidirectional interpolation is performed respectively. Before interpolation, the coordinates with the nearest value to the gap are obtained. Radial interpolation is performed by fixing the travel angle and finding the travel radius with the nearest value on both sides. The interpolation is performed accurately by combining the interpolation function (5) mentioned above. Tangential interpolation is performed by fixing the travel radius and finding the travel angle on both sides. The interpolation is performed by combining the selected interpolation function (6) mentioned above. After calculating the bidirectional interpolation results, the nearest angle and nearest radius of the two are combined as weights and the bidirectional fusion interpolation results at the gap are obtained by weighting by formula (7). Then, the filling value of the next gap is calculated, and the global travel time of the corresponding excitation point is filled in this way.
[0069] (5) (6) (7) The effects were compared by selecting different excitation points, before and after interpolation, as shown below. Figure 5 and Figure 6 As shown, the interpolation result is relatively smooth and complete. After filling in the global travel time of all excitation points, in order to investigate the interpolation fidelity, the PSF was analytically calculated based on the travel time with missing values and after interpolation, and its wavenumber spectrum was plotted as follows. Figure 7 and Figure 8 As shown in the figure, the PSF calculated analytically after interpolation is clearer and has a richer wavenumber spectrum. A comparison of the PSF of the two methods with the results of traditional migration imaging is shown below. Figure 9 It can be seen that the PSF correction effect of the analytical calculation of travel time after interpolation is better, and the layer of shielding under the anomalous high-velocity body is also revealed. Furthermore, it is more effective in suppressing imaging sidelobes, resulting in clearer layers. Therefore, the interpolation effect of this method is accurate and effective. This interpolation method, while ensuring efficiency, makes the technique of using ray tracing combined with Gaussian beam method to analytically calculate PSF for imaging correction in areas containing anomalous bodies more complete and has a wider range of applications.
[0070] This embodiment provides a computer-readable storage medium storing a computer program thereon. When executed by a processor, the computer program implements the steps of the aforementioned adaptive travel-time bidirectional fusion interpolation method. The storage medium may be a magnetic disk, optical disk, read-only memory (ROM), random access memory (RAM), flash memory, hard disk drive (HDD), or solid-state drive (SSD), etc.; the storage medium may also include combinations of the above-mentioned types of memory.
[0071] This embodiment provides a computer program product, including a computer program that, when executed by a processor, implements the steps of the above-described adaptive travel-time bidirectional fusion interpolation method.
[0072] The embodiments of the present invention have been described above with reference to the accompanying drawings. However, the present invention is not limited to the specific embodiments described above. The specific embodiments described above are merely illustrative and not restrictive. Those skilled in the art can make many other forms under the guidance of the present invention without departing from the spirit and scope of the claims. All of these forms are within the protection scope of the present invention.
Claims
1. An adaptive bidirectional fusion interpolation method for travel times, characterized in that, Includes the following steps: S1: Based on the velocity model of the target work area, the travel time of different excitation points in the entire work area, including the missing values, is calculated using the method of ray tracing combined with Gaussian beam. S2: Based on the velocity model of the target work area and the travel time results, perform tangential and radial multi-directional interpolation function fitting based on different excitation points, and select the optimal radial and tangential interpolation functions; S3: The travel time results containing missing values are split according to different excitation points to obtain travel time data. The travel time of each excitation point is calculated based on the critical angle of each excitation point. The travel time is filled into the missing points to obtain the travel time data in the Cartesian coordinate system after tangential interpolation boundary filling. The travel time data after tangential interpolation boundary filling is converted from the Cartesian coordinate system to the polar coordinate system to obtain the travel time data in the polar coordinate system. S4: Based on the travel time data in the Cartesian and polar coordinate systems, calculate the travel distance and travel angle corresponding to each missing value, as well as the boundaries of radial and tangential interpolation; S5: Based on the boundaries of the radial and tangential interpolation, using the optimal radial and tangential interpolation functions and the travel distance and travel angle corresponding to each missing value, perform comprehensive interpolation to obtain the interpolation result. Fill the missing space with the interpolation result to obtain the global travel time. S6: Calculate the PSF using the global travel time analysis, and perform deconvolution processing on the offset imaging results based on the PSF to obtain the least squares offset results.
2. The method according to claim 1, characterized in that, The calculation formula for the multi-directional interpolation function fitting is as follows: in, The sample value is the approximate travel time value; d is the ray path length from the excitation point to the target point; Indicates the root mean square velocity of the target point; This represents the multi-directional optimal interpolation function, i.e., different radial distances. and tangential angle The optimal interpolation function is composed of a basic function and a polynomial function; The computational notation for minimizing the function; It is the path distance r and the tangential angle. When traveling to a point; The interpolation error is represented by n; n is the number of global imaging points calculated during travel time. This represents the reference value for the true travel time of imaging point i; Based on different radial distances and tangential angle Calculated travel time.
3. The method according to claim 1, characterized in that, The formulas for calculating the travel distance and travel angle are as follows: in, and These refer to travel distance and travel angle, respectively. and These are the coordinates corresponding to the missing values in the Cartesian coordinate system. The formulas for calculating the boundaries of the radial interpolation and tangential interpolation are as follows: in, and These are the boundaries for radial interpolation and tangential interpolation, respectively. When representing the nearest radial distance, This indicates the shortest travel time with the tangential value.
4. The method according to claim 1, characterized in that, The formula for calculating the comprehensive interpolation is as follows: in, Indicates the travel time in the radial direction; Indicates the magnitude of travel time in the tangential direction; Represents the radial interpolation function that is optimal for the objective. Represents the optimal tangential interpolation function for the objective; Indicates the interpolation interval coefficients; These represent the nearest radial distance and the nearest tangential distance, respectively. The result is a comprehensive interpolation of the target.
5. An adaptive travel-time bidirectional fusion interpolation device, characterized in that, The device includes the following modules: The data interpolation preparation module is configured to: calculate the travel time of different excitation points in the entire work area, including missing values, based on the velocity model of the target work area using ray tracing combined with Gaussian beam method; The fitting function selection module is configured to: based on the velocity model of the target work area and the travel time results, perform tangential and radial multi-directional interpolation function fitting based on different excitation points, and select the optimal radial and tangential interpolation functions; The first data interpolation preprocessing module is configured to: split the travel time results containing missing values into different excitation points to obtain travel time data; calculate the travel time of each excitation point based on the critical angle of each excitation point; fill the missing travel time into the missing points to obtain travel time data in the Cartesian coordinate system after tangential interpolation boundary filling; and convert the travel time data after tangential interpolation boundary filling from the Cartesian coordinate system to the polar coordinate system to obtain travel time data in the polar coordinate system. The second data interpolation preprocessing module is configured to: calculate the travel distance and travel angle corresponding to each missing value, as well as the boundaries of radial interpolation and tangential interpolation, based on the travel time data in the Cartesian coordinate system and the polar coordinate system. The multi-directional interpolation fusion module is configured to: based on the boundaries of the radial and tangential interpolation, use the optimal radial and tangential interpolation functions and the travel distance and travel angle corresponding to each missing value to perform comprehensive interpolation to obtain the interpolation result, fill the missing space with the interpolation result, and obtain the global travel time; The travel time application module is configured to: use the global travel time analysis to calculate the PSF, and perform deconvolution processing on the migration imaging results based on the PSF to obtain the least squares migration result.
6. The apparatus according to claim 5, characterized in that, The calculation formula for the multi-directional interpolation function fitting is as follows: in, The sample value is the approximate travel time value; d is the ray path length from the excitation point to the target point; Indicates the root mean square velocity of the target point; This represents the multi-directional optimal interpolation function, i.e., different radial distances. and tangential angle The optimal interpolation function is composed of a basic function and a polynomial function; The computational notation for minimizing the function; It is the path distance r and the tangential angle. When traveling to a point; The interpolation error is represented by n; n is the number of global imaging points calculated during travel time. This represents the reference value for the true travel time of imaging point i; Based on different radial distances and tangential angle Calculated travel time.
7. The apparatus according to claim 5, characterized in that, The formulas for calculating the travel distance and travel angle are as follows: in, and These refer to travel distance and travel angle, respectively. and These are the coordinates corresponding to the missing values in the Cartesian coordinate system. The formulas for calculating the boundaries of the radial interpolation and tangential interpolation are as follows: in, and These are the boundaries for radial interpolation and tangential interpolation, respectively. When representing the nearest radial distance, This indicates the shortest travel time with the tangential value.
8. The apparatus according to claim 5, characterized in that, The formula for calculating the comprehensive interpolation is as follows: in, Indicates the travel time in the radial direction; Indicates the magnitude of travel time in the tangential direction; Represents the radial interpolation function that is optimal for the objective. Represents the optimal tangential interpolation function for the objective; Indicates the interpolation interval coefficients; These represent the nearest radial distance and the nearest tangential distance, respectively. The result is a comprehensive interpolation of the target.
9. A computer-readable storage medium having a computer program stored thereon, characterized in that, When the computer program is executed by a processor, it implements the steps of the adaptive travel-time bidirectional fusion interpolation method as described in any one of claims 1-4.
10. A computer program product, comprising a computer program, characterized in that, When the computer program is executed by a processor, it implements the steps of the adaptive travel-time bidirectional fusion interpolation method as described in any one of claims 1-4.