A multi-information fusion speed modeling method and device

Abstract

The application discloses a kind of multi-information fusion's speed modeling method and device, the method includes, after seismic migration, the imaging data A of acquisition, data phase control type mask B, seismic structure interpretation horizon and fault data C, well logging velocity curve data D and seismic velocity field data E;Edge feature of imaging data A is extracted by edge detection technique, and the structural skeleton data F of seismic image is obtained;Phase control type mask B is fused with structural skeleton data F, and phase control structure guide data G is obtained;Seismic structure interpretation data C is added in phase control structure guide data G, and space interpolation guide data H is formed;Through two-step interpolation extrapolation, the space interpolation distance weight field Q corresponding to space interpolation guide data H and well logging velocity interpolation filling model P are obtained;Classification mask M is established for well logging velocity interpolation filling model P;Two-step filling is carried out to class mask M, and multi-information fusion model R is obtained.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of seismic data processing, and particularly relates to a multi-information fusion velocity modeling method and device. BACKGROUND

[0002] The core method to realize complex structure and reservoir imaging is prestack depth migration technology. The accuracy of prestack depth migration imaging greatly depends on the reliability of depth domain velocity modeling. There are complex seismic and geological characteristics in China's complex exploration area, such as dramatic terrain undulation, fast change of surface exposed lithology, large lateral change of near-surface thickness and velocity, broken underground structure, developed faults, large stratigraphic dip, various occurrence, dramatic lateral change of velocity, and local vertical velocity inversion, and so on. It is urgent to play the advantages of different exploration data (geology, seismic, and non-seismic) and effectively cope with complexity through the technical means of multi-information fusion modeling.

[0003] Geological understanding (structure, lithology), seismic velocity field, and logging velocity fusion modeling is a technical challenge. Geological understanding is composed of horizon, fault, and lithological body boundary characteristics, and provides large-scale geological framework structure information. Seismic velocity field has medium-low frequency characteristics, good lateral continuity, obvious large suite velocity structure and form characteristics, but low vertical resolution, and mainly provides velocity information greater than the Fresnel zone scale. Logging velocity has high frequency characteristics and high vertical resolution, but small lateral detection range, and mainly provides small-scale velocity high frequency detail information in a local range. It can be seen that in the process of geological understanding (structure, lithology), seismic velocity field, and logging velocity fusion modeling, the fusion of information of different scales is a key technical problem to be solved.

[0004] There are two kinds of conventional ways to realize the fusion of geological structure, seismic velocity field, and logging velocity: one is to smooth the logging velocity to the scale comparable to the seismic velocity, and then to fill the model by extrapolating the smoothed logging velocity along the layer, and in the filling area, a velocity model completely consistent with the smoothed logging velocity details is obtained; the other is to learn the mapping relationship between the seismic velocity and the logging velocity by using the neural network, and then to apply the local mapping relationship to the whole model to obtain a velocity model with detailed features. These two technical ideas have three common problems: first, the role of seismic velocity field is not fully played; second, the complex situation at different positions under the ground is expressed by the "one-hole view" of the logging velocity at the well site, resulting in a large difference between the details of the model formed by direct interpolation and filling and the details of the model formed by global application of local mapping relationship; and third, neither of the two ideas considers the control of lithofacies belt. SUMMARY

[0005] The application provides a multi-information fusion velocity modeling method and device, which fully considers geological understanding (structure, lithology), seismic velocity field and logging velocity data characteristics, comprehensively utilizes geological structure and lithofacies belt characteristics to control logging velocity extrapolation range and form of well site, fully considers the role of seismic velocity field, and realizes adaptive weighted fusion of geological understanding (structure and lithofacies belt), logging velocity and seismic velocity field.

[0006] A multi-information fusion velocity modeling method, comprising the following steps,

[0007] Obtaining imaging data A after seismic migration processing;

[0008] Obtaining lithofacies classification interpretation data B, referred to as facies-controlled type mask plate;

[0009] Obtaining seismic structure interpretation horizon and fault data C;

[0010] Obtaining logging velocity curve data D and performing smoothing processing;

[0011] Obtaining seismic velocity field data E;

[0012] Extracting edge features of the imaging data A through edge detection technology to obtain structure skeleton data F of the seismic image;

[0013] Fusing the facies-controlled type mask plate B and the structure skeleton data F to obtain facies-controlled structure guide data G;

[0014] Adding the seismic structure interpretation data C to the facies-controlled structure guide data G to form spatial interpolation guide data H;

[0015] Obtaining spatial interpolation distance weight field Q and logging velocity interpolation filling model P corresponding to the spatial interpolation guide data H through two-step interpolation extrapolation;

[0016] Establishing a classification mask plate M for the logging velocity interpolation filling model P;

[0017] Obtaining a multi-information fusion model R through two-step filling of the classification mask plate M.

[0018] Preferably, the imaging data A is obtained through user-provided prestack seismic trace data volume and user-provided migration velocity field.

[0019] Preferably, the lithofacies classification interpretation data B is obtained through logging lithofacies interpolation extrapolation method, logging lithofacies neural network intelligent prediction method or seismic data seismic facies intelligent prediction method according to the imaging data A and user-provided work area core logging columnar chart / logging lithofacies chart data.

[0020] Preferably, the seismic structure interpretation horizon and fault data C is obtained by manually interpreting or neural network intelligent interpreting the imaging data A.

[0021] Preferably, the well logging velocity curve data D is obtained by smoothing the user-provided well logging velocity data.

[0022] Preferably, the smoothing process comprises median filtering the well logging velocity data; and Gaussian recursive filtering and smoothing the data obtained by the median filtering.

[0023] Preferably, the median filtering comprises sorting the values in a sliding window containing an odd number of points, and assigning the median value to the center point; and the shape of the sliding window of the median filtering comprises a line, a square, a circle, a cross or a special geometric shape provided by the user.

[0024] Preferably, the Gaussian recursive filtering step comprises forward filtering the median filtered data; and backward filtering the forward filtered data; wherein,

[0025] The forward filtering formula is:

[0026] GSForward(n) = a1*Input1(n) + a2*Input1(n-1) - a3*GSForward(n-1) - a4*GSForward(n-2)

[0027] wherein GSForward is the data generated by the forward filtering, Input1 is the input data for the forward filtering, which is the median filtered data herein, n is the data sample position, and a1, a2, a3 and a4 are filter coefficients, which are set according to the actual needs of different work areas;

[0028] The backward filtering formula is:

[0029] GSBackward(n) = b1*Input2(n) + b2*Input2(n+1) - b3*GSBackward(n+1) - b4*GSBackward(n+2)

[0030] wherein GSBackward is the data generated by the backward filtering, Input2 is the input data for the backward filtering, which is the forward filtered data herein, n is the data sample position, and b1, b2, b3 and b4 are filter coefficients, which are set according to the actual needs of different work areas.

[0031] Preferably, the seismic velocity field data E is migration velocity field data corresponding to the imaging data A after seismic migration processing, and is obtained by processing the user-provided pre-stack seismic trace data volume Data.

[0032] Preferably, the spatial interpolation distance field Q corresponding to the spatial interpolation guide data H and the logging velocity interpolation filling model P obtained by two-step interpolation extrapolation include,

[0033] Adding the logging velocity curve data D in the spatial interpolation guide data H, judging the position relationship of the logging data space and time range and the nearest neighbor facies-controlled area of the facies-controlled type mask B.

[0034] Preferably, the specific method of the two-step interpolation extrapolation is as follows: according to the local structure skeleton features provided by the spatial interpolation guide data H,

[0035] The logging data falling into the facies-controlled area is processed by the radial basis function method;

[0036] The extrapolation relationship of all logging data is established by the Kriging interpolation method, and the logging velocity along the structure and the related facies belt is extrapolated and filled.

[0037] Preferably, the formula of the radial basis function method is as follows:

[0038]

[0039] Wherein, ε represents the scale factor of the radial basis function, which is 1 here; r represents the Euclidean distance between the interpolation point and the well position.

[0040] Preferably, the region division of the classification mask M is divided into a 0-value indicating area and a 1-value indicating area;

[0041] The 0-value indicating area represents an unfilled area of the logging velocity interpolation filling model P;

[0042] The 1-value indicating area represents a filling area of the velocity value after the weighted fusion of the seismic velocity model and the logging velocity interpolation filling model by the weighted least square method.

[0043] A device of a multi-information fusion velocity modeling method, the device comprising,

[0044] The information collection device is used to obtain the imaging data A after seismic migration processing, the lithofacies classification interpretation data B, the seismic structure interpretation horizon and fault data C, the logging velocity curve data D, and the seismic velocity field data E;

[0045] The information processing device is used to extract the edge features of the imaging data A by an edge detection technology, obtain the structure skeleton data F of the seismic image, fuse the facies-controlled type mask B with the structure skeleton data F, obtain the facies-controlled structure guide data G, add the seismic structure interpretation data C in the facies-controlled structure guide data G, and form the spatial interpolation guide data H;

[0046] The information modeling device is used for obtaining a spatial interpolation distance field Q corresponding to the spatial interpolation guide data H and a logging velocity interpolation filling model P by two-step interpolation extrapolation, establishing a classification mask M for the logging velocity interpolation filling model P, and obtaining a multi-information fusion model R by two-step filling on the classification mask M.

[0047] Technical effects and advantages of the present application:

[0048] The present application fully considers the geological understanding (structure, lithology), seismic velocity field and logging velocity data characteristics, on the one hand, comprehensively utilizes the geological structure and lithofacies belt characteristics to control the logging velocity extrapolation range and shape at the well site, and on the other hand, fully considers the role of the seismic velocity field, realizes the adaptive weighted fusion of the geological understanding (including structure and lithofacies belt), logging velocity and seismic velocity field.

[0049] Other features and advantages of the present application will be set forth in the following description, and in part will become apparent to those skilled in the art from the description, or can be learned by practice of the present application. The objects and other advantages of the present application will be realized and attained by the structure particularly pointed out in the description, claims and drawings. BRIEF DESCRIPTION OF DRAWINGS

[0050] Figure 1 It is a multi-information fusion velocity modeling method flowchart provided by the present application;

[0051] Figure 2 It is a seismic imaging data diagram used in the present embodiment;

[0052] Figure 3a It is a schematic diagram of a seismic attribute representing lithofacies provided for a user in the present embodiment;

[0053] Figure 3b It is a schematic diagram of a closed area representing lithofacies formed by artificial intelligence geologic body sculpture in the present embodiment;

[0054] Figure 4 It is a seismic structure interpretation data diagram used in the present embodiment;

[0055] Figure 5 It is a position relationship diagram of four logging velocity curve data and seismic imaging data used in the present embodiment;

[0056] Figure 6 It is a seismic velocity field data diagram used in the present embodiment;

[0057] Figure 7a It is a superimposed display diagram of seismic velocity, original logging velocity and smoothed logging velocity at the 1# logging position in the present embodiment;

[0058] Figure 7bFigure 6 is a superimposed display diagram of the seismic velocity, the original logging velocity and the smoothed logging velocity at the 2# logging position in the present embodiment;

[0059] Figure 7c Figure 7 is a superimposed display diagram of the seismic velocity, the original logging velocity and the smoothed logging velocity at the 3# logging position in the present embodiment;

[0060] Figure 7d Figure 8 is a superimposed display diagram of the seismic velocity, the original logging velocity and the smoothed logging velocity at the 4# logging position in the present embodiment;

[0061] Figure 8 Figure 9 is a spatial interpolation distance weight field diagram obtained in the present application;

[0062] Figure 9 Figure 10 is a logging velocity interpolation filling model diagram obtained in the present application;

[0063] Figure 10 Figure 11 is a relative medium frequency fusion velocity model diagram obtained in the present application;

[0064] Figure 11 Figure 12 is a relative high frequency fusion velocity model diagram obtained in the present application;

[0065] Figure 12 Figure 13 is a comparison diagram of the seismic velocity, the logging interpolation velocity and the relative medium frequency fusion velocity at the logging position in the present application;

[0066] Figure 13 Figure 14 is a comparison diagram of the seismic velocity, the logging interpolation velocity and the relative high frequency fusion velocity at the logging position in the present application. DETAILED DESCRIPTION

[0067] The technical solutions in the embodiments of the present application will be clearly and completely described below with reference to the drawings in the embodiments of the present application. Obviously, the described embodiments are only part of the embodiments of the present application, rather than all the embodiments of the present application. Based on the embodiments in the present application, all the other embodiments obtained by those skilled in the art without creative labor fall within the protection scope of the present application.

[0068] To solve the problems of the prior art, the application discloses a multi-information speed fusion method and device, the device comprises: an information collection device for obtaining imaging data A after seismic migration processing, lithofacies classification interpretation data B, seismic structure interpretation horizon and fault data C, well logging velocity curve data D, and seismic velocity field data E; an information processing device for extracting edge features of the imaging data A through edge detection technology to obtain structure skeleton data F of a seismic image, fusing the phase control type mask B with the structure skeleton data F to obtain phase control structure guiding data G, adding the seismic structure interpretation data C in the phase control structure guiding data G to form spatial interpolation guiding data H; and an information modeling device for obtaining spatial interpolation distance weight field Q corresponding to the spatial interpolation guiding data H and well logging velocity interpolation filling model P through two-step interpolation extrapolation, establishing a classification mask M for the well logging velocity interpolation filling model P, and obtaining a multi-information fusion model R through two-step filling of the classification mask M.

[0069] As shown in Figure 1 the method of the application comprises the following steps: obtaining imaging data A after seismic migration processing; obtaining lithofacies classification interpretation data B, referred to as a phase control type mask; obtaining seismic structure interpretation horizon and fault data C; obtaining well logging velocity curve data D and performing smoothing processing; obtaining seismic velocity field data E; extracting edge features of the imaging data A through edge detection technology to obtain structure skeleton data F of a seismic image; fusing the phase control type mask B with the structure skeleton data F to obtain phase control structure guiding data G; adding the seismic structure interpretation data C in the phase control structure guiding data G to form spatial interpolation guiding data H; obtaining spatial interpolation distance weight field Q corresponding to the spatial interpolation guiding data H and well logging velocity interpolation filling model P through two-step interpolation extrapolation; establishing a classification mask M for the well logging velocity interpolation filling model P; and obtaining a multi-information fusion model R through two-step filling of the classification mask M.

[0070] Further, the imaging data A after seismic migration processing is obtained, as shown in Figure 2As shown in the figure, Xline (trace) represents seismic traces, time represents time, Amplitude represents amplitude, and the position of the light and dark change in the image indicates the position of the underground reflection layer and the response amplitude. The response amplitude of the reflection layer at different positions in the figure is related to the corresponding geological lithology of the reflection layer, and different geological lithologies have different absorption and attenuation of seismic waves, so the imaging response amplitudes are different. Commonly used seismic prestack depth migration methods include ray-based migration and wave equation migration. The ray-based method mainly includes Kirchhoff migration, Gaussian beam migration, and control beam migration. The wave equation-based method mainly includes one-way wave migration and reverse time migration. Both ray-based migration and wave equation-based migration are based on wave equation theory, but the difference lies in that the ray-based migration uses geometric ray theory to calculate the amplitude and phase information of the wave field, thereby realizing the extension imaging of the wave field, while the wave equation-based migration is based on the numerical solution of the wave equation. Generally speaking, the wave equation-based migration has higher imaging accuracy, while the ray-based migration has higher computational efficiency and flexibility. Mainstream seismic data processing commercial software, such as CGG, Omega, AGT, Paradigm, etc. provides processing modules for two types of migration technology. At the same time, open source software Seismic Unix and Madagascar also provide processing modules for two types of migration technology. Under the condition that the user provides prestack seismic trace data body Data and migration velocity field Velocity, the user can choose the migration algorithm module of the processing software at will to obtain the imaging data A after the seismic migration processing, as follows:

[0071] A = ImageModule (Data, Velocity) (1)

[0072] Where ImageModule is the imaging processing module of commercial software or open source software, Data is the prestack seismic trace data body provided by the user, Velocity is the migration velocity field provided by the user, and A is the imaging data body obtained after the seismic migration processing.

[0073] Further, the lithofacies classification interpretation data B is obtained, which is called facies-controlled type mask. The methods for obtaining lithofacies classification interpretation data include well logging lithofacies interpolation and extrapolation method, well logging lithofacies neural network intelligent prediction method, and seismic data seismic facies intelligent prediction method. According to the actual data of different work areas, the appropriate method module can be selected to generate lithofacies classification interpretation data B. The work flow is described as follows: first, the core logging columnar chart / well logging lithofacies chart of the drilling / logging position in the work area Well_Lithology_Model is calibrated and matched with the imaging data A after the seismic migration processing, then the lithofacies classification interpretation data B is obtained by the well logging lithofacies interpolation and extrapolation method or the well logging lithofacies neural network intelligent prediction method or the seismic data seismic facies intelligent prediction method, as follows:

[0074] B=LithologyModule(A,Well_Lithology_Model) (2)

[0075] In the formula, LithologyModule is the processing module for well logging lithofacies interpolation and extrapolation, well logging lithofacies neural network intelligent prediction, or seismic data seismic facies intelligent prediction; A is the imaging data volume obtained after seismic migration processing; Well_Lithology_Model is the core logging columnar section / well logging lithofacies map data provided by the user; and B is the lithofacies classification interpretation data volume obtained by the lithofacies classification interpretation.

[0076] Phase-controlled masks can be such as Figure 3a The image shows the seismic properties of the lithofacies. The facies-controlled mask can also be as follows: Figure 3b The engraved geological body shown represents a closed region representing lithofacies. It can also be various types of lithofacies zone features or attribute features provided by non-seismic methods (gravity, magnetic, electrical, geochemical methods). In the figure, Xline (trace) represents a seismic trace, and time represents time. The facies-controlled mask mentioned in this example is only used to provide key geological lithofacies boundary features, and the implementation scheme is not limited.

[0077] Obtain seismic tectonic interpretation horizon and fault data C. There are two main categories of methods for obtaining seismic tectonic interpretation horizon and fault data: manual interpretation and neural network intelligent interpretation. The appropriate method module can be selected to generate the seismic tectonic interpretation horizon and fault data C based on the actual data conditions of different work areas. The workflow is described as follows: The imaging data A obtained after seismic migration processing is interpreted using either manual interpretation or neural network intelligent interpretation to obtain the seismic tectonic interpretation horizon and fault data C. The method is as follows:

[0078] C = StructureModule(A) (3)

[0079] In the formula, StructureModule is a processing module for manual interpretation or neural network intelligent interpretation; A is the imaging data volume obtained after seismic migration processing; and C is the seismic tectonic interpretation horizon and fault data obtained from the seismic tectonic interpretation. For example... Figure 4 As shown in this embodiment, there are a total of 2 faults and 10 layers. In the figure, Xline (trace) represents the seismic trace, time represents time, and Amplitude represents amplitude.

[0080] Furthermore, the well logging velocity curve data D is acquired and appropriately smoothed; for example... Figure 5As shown in the figure, Xline (trace) represents seismic trace, time represents time, and Amplitude represents amplitude. The figure shows the position relationship of the 4 well logging (from left to right, 1# well logging, 2# well logging, 3# well logging and 4# well logging) velocity curve data and the calibrated seismic imaging data in the area. Horizontally, the 4 well logs are located at different positions. Vertically, through well-to-seismic calibration interpretation, it is determined that the starting positions of the 4 well logging velocity curves are different. At the same time, the well logging velocity curve is represented by a vertical bar graph, and the light and dark relationship thereon reflects the change of the velocity value with depth. After the user obtains the well logging velocity data WellVelocity, the user can obtain the smoothed well logging velocity data D through the filter processing module of the mainstream commercial seismic data processing software CGG, Omega, AGT, Paradigm, etc. or open source software Seismic Unix and Madagascar, etc. The smoothing degree should be selected according to the actual requirements of different work areas, which is not limited here. The implementation method is as follows:

[0081] D = FilterModule (WellVelocity) (4)

[0082] In the formula, FilterModule is a filter processing module of commercial software or open source software, WellVelocity is well logging velocity data provided by the user, and D is smoothed well logging velocity data obtained after filtering.

[0083] Generally, well logging velocity data is discretely distributed in space without any geometric structure, and may not even be located on a uniformly sampled seismic grid. In order to obtain a smooth well logging velocity data, the well logging velocity data is smoothed. Figure 5 The well logging (from left to right, 1#, 2#, 3# and 4#) velocity curve data is smoothed, and the relationship between the well logging velocity before and after smoothing and the seismic velocity is shown in Figure 7a-7d The figure shows that the horizontal label velocity represents velocity, and Time represents time. According to Figure 7a-7d It can be seen that the curve with greater fluctuation represents the original well logging velocity data, the smooth line represents the seismic velocity data at the corresponding position of the well point, and the curve with smaller fluctuation represents the smoothed well logging velocity data. The curve with smaller fluctuation is in the middle of the curve with greater fluctuation. The comparison of the three curves shows that the original well logging velocity data and the seismic velocity data (compare the curve with greater fluctuation and the smooth line) have a large difference in frequency component. The seismic velocity data presents a low frequency feature, and the original well logging velocity data presents an extremely high frequency feature. The smoothed well logging velocity data (curve with smaller fluctuation) contains some velocity detail information, and the frequency component is not too high to deviate too much from the seismic velocity.

[0084] The smoothing process described herein includes two steps, first, median filtering is performed on the logging velocity data, and then Gaussian recursive filtering is performed on the median filtered data to smooth the data. The median filtering generally uses a sliding window with an odd number of points, and the median value of the values in the window is used to replace the value of the center point, that is, the values in the window are sorted, and then the median value is assigned to the center point. Commonly used median filtering window shapes include linear, square, circular, cross-shaped or special geometric shapes provided by the user, and the working method is described as follows:

[0085] Median filtered data = MedianFilter (window shape, window size) (5)

[0086] where MedianFilter is a median filter provided by commercial software or open source software.

[0087] The specific steps of the Gaussian recursive filter are divided into two steps: first, forward filtering is performed on the median filtered data obtained by MedianFilter, and second, backward filtering is performed on the forward filtered data obtained in the first step. The forward filtering formula is:

[0088] GSForward(n) = a1*Input1(n) + a2*Input1(n-1) - a3*GSForward(n-1) - a4*GSForward(n-2) (6)

[0089] where GSForward is the data generated by forward filtering, Input1 is the forward filtering input data, which is the median filtered data obtained by formula 5, n is the data sample position, and a1, a2, a3, a4 are filter coefficients, which are set according to the actual needs of different work areas;

[0090] The backward filtering formula is:

[0091] GSBackward(n) = b1*Input2(n) + b2*Input2(n+1) - b3*GSBackward(n+1) - b4*GSBackward(n+2) (7)

[0092] where GSBackward is the data generated by backward filtering, Input2 is the backward filtering input data, which is the forward filtered data obtained by formula 6, n is the data sample position, and b1, b2, b3, b4 are filter coefficients, which are set according to the actual needs of different work areas.

[0093] Further, the seismic velocity field data E is obtained, as shown in Figure 6 Figure 6 ​A seismic velocity model is shown, which is a layer velocity model, presenting low velocity in shallow part, high velocity in deep part, medium-high velocity in the left side of the fault, and medium-low velocity in the right side of the fault. The seismic velocity field data is migration velocity field data corresponding to the imaging data A after seismic migration processing. Under the condition of user-provided pre-stack seismic trace data volume Data, the migration velocity modeling processing module of mainstream commercial seismic data processing software, such as CGG, Omega, AGT, Paradigm, etc., or open source software Seismic Unix and Madagascar, etc., can be directly obtained, and the method is as follows:

[0094] E = VelocityModule(Data) (8)

[0095] In the formula, VelocityModule is the migration velocity modeling processing module of commercial software or open source software, Data is the pre-stack seismic trace data volume provided by the user, and E is the seismic velocity field data volume obtained by migration velocity modeling processing.

[0096] Further, the edge features of the imaging data A are extracted by edge detection technology to obtain the structural skeleton data F of the seismic image; the phase-controlled type mask B is fused with the structural skeleton data F to obtain the phase-controlled structural guide data G; and the seismic structural interpretation data C is added in the phase-controlled structural guide data G to form the spatial interpolation guide data H.

[0097] Further, in order to obtain the spatial interpolation distance weight field as shown in Figure 8 and the well logging velocity interpolation filling model as shown in Figure 9 , in which Xline (trace) represents seismic trace, Distance field represents distance field, interpolated velocity represents well logging velocity interpolation, and time represents time, Figure 8 The light and dark changes in the spatial interpolation distance weight field diagram shown in Figure 9 represent the interpolation weight size at different positions from the well point. Generally, the maximum weight is set near the well site, and the farther from the well point, the smaller the interpolation weight; Figure 8The weight relationship extrapolation interpolation is obtained by the space interpolation distance weight field. After judging the space and time range of the logging data and the position relationship with the main phase control area of the phase control type mask, firstly, the logging data falling into the main phase control area is processed by extrapolation, that is, in the phase control area, according to the local structure skeleton characteristics provided by the space interpolation guide data H, the logging velocity at the well site is locally extrapolated and filled by the radial basis function method, and the extrapolation process is extended to the boundary of the phase control area and terminated; secondly, according to the overall structure guide characteristics provided by the space interpolation guide data H, the extrapolation relationship of all logging data is established by the Kriging interpolation method, and the logging velocity is extrapolated and filled along the structure and the related phase belt. Through the above two steps of interpolation and extrapolation, the space interpolation distance weight field Q corresponding to the space interpolation guide data H and the logging velocity interpolation filling model P are obtained.

[0098] The radial basis function method formula is as follows:

[0099]

[0100] Wherein, ε represents the scale factor of the radial basis function, which is 1 here; r represents the Euclidean distance between the interpolation point and the well position. The velocity interpolated by the method has structure information.

[0101] Further, the radial basis function here can have many choices, such as choosing a conventional radial basis function, choosing a radial basis function neural network, or designing a new radial basis function according to specific needs. In the embodiment, a distance-related radial basis function is designed, and the expression is as follows:

[0102]

[0103] In formula (10), it is assumed that there are N wells in the study area, the i-th well falls into the main phase control area, and there are K i sampling points on the velocity curve segment of the i-th well, (x k , y k , z k ) represents the coordinates of each sampling point on the velocity curve segment, (x, y, z) represents the coordinates of the i-th well falling into the main phase control area, w i (x, y, z) represents the local space interpolation distance weight field generated by the radial basis function (10) at the i-th well (x, y, z). ∈ is a radial parameter used to control the calculation accuracy and balance the relationship between different weighting schemes.

[0104] Secondly, according to the overall structure guiding features provided by the spatial interpolation guide data, an extrapolation relationship between all logging data is established by using the Kriging interpolation method, and a spatial interpolation distance weight field of a point (x, y, z) to be interpolated outside the main facies control area and related to the structure and facies belt is generated, and the spatial interpolation distance weight field generated by the radial basis function and the Kriging interpolation function is shown as w(x, y, z) in Figure 8 The spatial interpolation distance weight field is related to the position of the well and is mainly used to consider the influence of the distance from the well. This is because the velocity far away from the well is unreliable when extrapolating the well. Using the spatial interpolation distance weight field can better control the influence of the velocity far away from the well on the model fusion. Figure 8 It can be seen from the figure that the distribution of the distance weight field is consistent with the distribution of the well.

[0105] Suppose that the i-th well logging velocity data is represented by v i (x, y, z), and the well logging velocity interpolation filling model of the whole area N wells is calculated and generated according to the following formula:

[0106]

[0107] Figure 9 The well logging velocity interpolation filling model is calculated by formula (11), and it can be seen from the figure that the well logging velocity interpolation filling model can well reflect the comprehensive features of the horizon, fault, seismic facies and lithofacies, has high vertical resolution and rich model details, and provides a good data basis for subsequent multi-information velocity fusion.

[0108] Further, a classification mask M is established for the well logging velocity interpolation filling model P, in which the area filled with the well logging velocity is indicated by 1, and the area not filled with the well logging velocity is indicated by 0; according to the area division of the well logging velocity interpolation filling model classification mask M, the 0 value indicated area of the M mask is filled by using the seismic velocity field based on the well logging velocity interpolation filling model P, and in the 1 value indicated area of the M mask, the weighted fusion of the seismic velocity model and the well logging velocity interpolation filling model is performed by using the weighted least square method, and the area is filled with the fused velocity value. Through the above two steps of filling, the multi-information fusion model R proposed in the present application is finally obtained.

[0109] Further, the seismic velocity field is obtained, the well logging velocity interpolation filling model and the seismic velocity model are fused, a relatively medium frequency fused velocity model as shown in Figure 10 and a relatively high frequency fused velocity model as shown in Figure 11 are obtained, in which Xline (trace) represents a seismic trace, time represents time, and Merged velocity represents fused velocity. Figure 10The fusion velocity model shown presents medium frequency, which reflects both the information of the background velocity field and the details of the logging velocity. Figure 11 The fusion velocity model shown presents high frequency, which mainly reflects the characteristics of the logging velocity interpolation and filling model and can well reflect the details. Figure 10 And Figure 11 The two velocity models shown are consistent with the spatial variation trend and the structural features.

[0110] We express the fusion of multiple velocity models as the process of estimating a unified velocity model that fits the known multiple velocity models. The fusion method provided by the present application does not limit the number and dimensions of the models, can be any multiple velocity models, and can be one-dimensional, two-dimensional, three-dimensional or multi-dimensional models. The fusion method can be similarly extended according to the number and dimensions of the velocity models participating in the fusion. Meanwhile, the present application proposes to implement velocity fusion by using a weighted least square method, which can be a conventional weighted least square method, an improved weighted least square method or a weighted least square method implemented by using a neural network, and the technical effects of the three are the same, which are not limited here. In the embodiment, we design an improved weighted least square method, and the process of generating the fusion model v(x) from two one-dimensional velocity models v1(x) and v2(x) is expressed by the following formula:

[0111]

[0112] In formula (12), it is assumed that v1(x) represents a seismic velocity model ( Figure 6 ) and v2(x) represents a logging velocity interpolation and filling model ( Figure 9 ). w1(x) and w2(x) are space-varying quality maps, and the value ranges are both [0, 1], which respectively represent the confidence of the velocity models v1(x) and v2(x) in space. For example, in the v2(x) velocity model, the confidence of the velocity in the area far from the well or the area where the well data is missing is relatively low, and w2(x) should be set to a value close to 0 in the area, giving a smaller weight in the fusion of multiple velocity models. Therefore, w2(x) can be defined as a function of the distance from a point x in space to the nearest well, and the greater the distance, the smaller w2(x). In addition, in the fusion process of multiple velocity models, in order to eliminate artificial artifacts or discontinuities generated in the fusion process, we introduce a regularization term D. In the area where the two known velocity models v1 and v2 conflict (or have large differences), the regularization term helps to obtain a final fusion velocity model v(x) with a smooth transition. Finally, λ1, λ2 and ∈ are all constants in the range of 0 to 1, which are used to balance the equations.

[0113] Let Equation (12) can be simplified as:

[0114]

[0115] To solve the least squares solution of the multiple equations in the above equation (13), we solve its corresponding normal equations as follows:

[0116] (L T L+∈ 2 DTD)v=L T d (14)

[0117] The solution of this equation can be written as follows:

[0118] v=(L T L+∈ 2 D T D) -1 L T d (15)

[0119] To better introduce the existing geological understanding in the velocity fusion process and obtain the fusion velocity model we expect, we introduce a regularization term S in the above solving process:

[0120] S=(I+∈ 2 D T D) -1 (16)

[0121] This regularization term can be implemented as a general smoothing operator that is used to effectively introduce our existing understanding. For example, the velocity model we want to estimate should be discontinuous on both sides of the fault and should have certain spatial continuity along the stratigraphic structure direction. At this time, we can implement S as a spatially variant smoothing operator that is used to enhance and maintain the discontinuity of the velocity model at the fault location and ensure the spatial continuity of the velocity model along the stratigraphic structure direction. The above regularization term can be rewritten as

[0122] ∈ 2 D T D=S -1 -I (17)

[0123] Bringing it into equation (15), the following solving form based on the regularization operator can be obtained:

[0124] v=(L T L+S -1 -I) -1 L T d=[I+S(L T L-I)] -1 SL T d (18)

[0125] We find that when S = I, v = (L T L)-1 L T d, where the fusion velocity is the unregularized least squares solution. When (L T L) = I, v = SL T d, where no inversion is needed. When S = λI, in the limit λ→0, v ≈ λL T d. Since in practical problems L can have physical units, a scaling factor 1 / λ is needed, and equation (18) becomes:

[0126] v = [λ 2 I + S(L T L - λ 2 I)] -1 SL T d (19)

[0127] Since in iterative inversion with the conjugate gradient method, one usually needs a positive definite and symmetric operator, let S = HH T , where

[0128] v = H[λ 2 I + H T (L T L - λ 2 1)] -1 H T L T d (20)

[0129] We can invert equation (20) with the conjugate gradient method to obtain the velocity fusion result.

[0130] Figure 6 A seismic velocity model is shown in the figure, where Xline (trace) represents a seismic trace, time represents time, and vint represents the layer velocity in the time domain. The velocity model is a layer velocity model, which presents low velocity characteristics in the shallow part, high velocity characteristics in the deep part, medium-high velocity characteristics in the left range of the fault, and medium-low velocity characteristics in the right range of the fault. The velocity difference near the fault is obvious. We take two kinds of weights to fuse the seismic velocity model shown in Figure 6 and the well logging velocity interpolation filling model shown in Figure 9 , and obtain a relatively medium frequency fusion velocity model ( Figure 10 ) and a relatively high frequency fusion velocity model ( Figure 11 ). Figure 10 The fusion velocity model shown in Figure 11 presents medium frequency, which not only reflects the background velocity field information, but also embodies the details of the well logging velocity. Figure 10 The fusion velocity model shown in Figure 11 presents high frequency, mainly embodies the characteristics of the well logging velocity interpolation filling model, and can well reflect the details.The two velocity models shown are consistent with the spatial variation trend and the structural features.

[0131] We also compared the two fusion velocity models with the seismic velocity model and the well-logging velocity interpolation filling model at the well location. Figure 12 is the comparison of the seismic velocity, the well-logging interpolation velocity and the relative medium frequency fusion velocity at the well location in the present application, in which the horizontal label Velocity represents velocity, the vertical label Time represents time, Figure 13 is the comparison of the seismic velocity, the well-logging interpolation velocity and the relative high frequency fusion velocity at the well location in the present application, in which the horizontal label Velocity represents velocity, the vertical label Time represents time, and Figure 12 and Figure 13 We can find in and that: in the range of the upper and lower sides of the well, the fusion velocity (the longest line) is consistent with the seismic velocity (the smooth line), which can effectively ensure the stability of the fusion velocity. In the well-logging curve area, the fusion velocity (the longest line) is between the well-logging interpolation velocity (the line with the largest fluctuation) and the seismic velocity (the smooth line). And the overall trend change of the fusion velocity (the longest line) is consistent with the seismic velocity (the smooth line). In the detail part, the waveform change of the fusion velocity (the longest line) is consistent with the well-logging interpolation velocity (the line with the largest fluctuation). Figure 12 and Figure 13 The velocity waveform comparison provided at the well location better confirms the reliability of the fusion method proposed by us.

[0132] Finally, it should be noted that: the above only for the preferred embodiments of the present application, and not for the limitation of the present application, although the foregoing detailed description of the present application is made with reference to the foregoing embodiments, for the person skilled in the art, it still can be modified to the technical solutions recorded in the foregoing each embodiment, or to the equivalent replacement of part of the technical features, any modification, equivalent replacement, improvement, etc. made within the spirit and principles of the present application, should be included in the protection scope of the present application.

Claims

1. A speed modeling method based on multi-information fusion, characterized in that, The method Includes the following steps, Acquire the seismic migration-processed imaging data A; The lithofacies classification interpretation data B is obtained and is called the facies type mask. Obtain seismic tectonic interpretation horizon and fault data C; Acquire the logging velocity curve data D and perform smoothing processing; Obtain earthquake velocity field data E; Edge features of the imaging data A are extracted using edge detection technology to obtain the structural skeleton data F of the seismic image; The phase control type mask B is fused with the structural skeleton data F to obtain phase control structural guidance data G; Seismic tectonic interpretation data C is added to the phased structure guidance data G to form spatial interpolation guidance data H; Through two-step interpolation extrapolation, the spatial interpolation distance weight field Q and the logging velocity interpolation filling model P corresponding to the spatial interpolation guidance data H are obtained; where... The process of obtaining the spatial interpolation distance weight field Q and logging velocity interpolation filling model P corresponding to the spatial interpolation guidance data H through two-step interpolation extrapolation includes adding logging velocity curve data D to the spatial interpolation guidance data H, and determining the spatial and temporal range of the logging data and its positional relationship with the nearest neighbor phase control area of ​​the phase control type mask B. The specific method for the two-step interpolation extrapolation is as follows: The local skeleton features are constructed based on the spatial interpolation-guided data H. Extrapolation of logging data falling into the phase control region is performed using the radial basis function method; Extrapolation relationships for all logging data are established using the Kriging interpolation method, enabling the extrapolation and filling of logging velocities along structural and related facies zones. A classification mask M is established for the well logging velocity interpolation filling model P; A multi-information fusion model R is obtained by performing a two-step filling process on a mask-like surface M; where... The classification mask M is divided into a 0-value indication area and a 1-value indication area; The 0-value indication area represents the unfilled area of ​​the logging velocity interpolation filling model P; The 1-value indication area represents the area filled with velocity values ​​after weighted fusion of the seismic velocity model and the well logging velocity interpolation filling model using the weighted least squares method.

2. The speed modeling method for multi-information fusion according to claim 1, characterized in that, The imaging data A was acquired using pre-stack seismic gather data volume and migration velocity field provided by the user.

3. The speed modeling method for multi-information fusion according to claim 1, characterized in that, The lithofacies classification interpretation data B is obtained based on imaging data A and core logging columnar section / well logging lithofacies data provided by the user, through well logging lithofacies interpolation and extrapolation, well logging lithofacies neural network intelligent prediction, or seismic data seismic facies intelligent prediction.

4. The speed modeling method for multi-information fusion according to claim 1, characterized in that, The seismic tectonic interpretation horizon and fault data C are obtained by interpreting the imaging data A through manual interpretation or neural network intelligent interpretation.

5. The speed modeling method for multi-information fusion according to claim 1, characterized in that, The logging velocity curve data D is obtained by smoothing the logging velocity data provided by the user.

6. The speed modeling method for multi-information fusion according to claim 5, characterized in that, The smoothing process includes: performing median filtering on the logging velocity data; and performing Gaussian recursive filtering smoothing on the data obtained from the median filtering.

7. The speed modeling method for multi-information fusion according to claim 6, characterized in that, The median filtering involves sorting the values ​​of a sliding window containing an odd number of points and assigning the median value to the center point; the shape of the sliding window for median filtering includes linear, square, circular, cross-shaped, or user-provided special geometric shapes.

8. The speed modeling method for multi-information fusion according to claim 6, characterized in that, The Gaussian recursive filtering step includes: performing forward filtering on the median-filtered data; and performing backward filtering on the forward-filtered data; wherein... The forward filtering formula is: GSForward(n)=a1 Input1(n)+a2 Input1(n-1)-a3 GSForward(n-1) -a4 GSForward(n-2) In the formula, GSForward is the data generated by forward filtering, Input1 is the input data of forward filtering, which is median filtering data here, n is the data sample point position, and a1, a2, a3, and a4 are filter coefficients, which are set according to the actual needs of different work areas. The backward filtering formula is as follows: GSBackward (n)=b1 Input2(n)+b2 Input2(n+1)-b3 GSBackward (n+1)-b4 GSBackward(n+2) In the formula, GSBackward is the data generated by backward filtering, Input2 is the input data for backward filtering (here, it is the data for forward filtering), n is the data sample point position, and b1, b2, b3, and b4 are filter coefficients, which are set according to the actual needs of different work areas.

9. The speed modeling method for multi-information fusion according to claim 1, characterized in that, The seismic velocity field data E is the migration velocity field data corresponding to the seismic migration processed imaging data A, obtained by processing the pre-stack seismic gather data volume Data provided by the user.

10. The speed modeling method for multi-information fusion according to claim 1, characterized in that, The formula for the radial basis function method is as follows: Where ε represents the scaling factor of the radial basis function, which is set to 1 here; r represents the Euclidean distance between the interpolation point and the well location.

11. An apparatus for a speed modeling method based on multi-information fusion, characterized in that, The device includes, The information collection device is used to acquire the following data: A) imaging data after seismic migration processing; B) lithofacies classification and interpretation data; C) seismic tectonic interpretation of stratigraphic and fault data; D) well logging velocity curve data; E) seismic velocity field data. The information processing device is used to extract the edge features of the imaging data A through edge detection technology to obtain the structural skeleton data F of the seismic image, fuse the lithofacies classification interpretation data B with the structural skeleton data F to obtain the facies-controlled structural guidance data G, and add the seismic structural interpretation data C to the facies-controlled structural guidance data G to form spatial interpolation guidance data H. The information modeling device is used to obtain the spatial interpolation distance weight field Q and the logging velocity interpolation filling model P corresponding to the spatial interpolation guidance data H through two-step interpolation extrapolation. A classification mask M is established for the logging velocity interpolation filling model P, and a multi-information fusion model R is obtained by performing two-step filling on the classification mask M. The process of obtaining the spatial interpolation distance weight field Q and logging velocity interpolation filling model P corresponding to the spatial interpolation guidance data H through two-step interpolation extrapolation includes adding logging velocity curve data D to the spatial interpolation guidance data H, and determining the spatial and temporal range of the logging data and its positional relationship with the nearest neighbor phase control area of ​​the phase control type mask B. The specific method for the two-step interpolation extrapolation is as follows: The local skeleton features are constructed based on the spatial interpolation-guided data H. Extrapolation of logging data falling into the phase control region is performed using the radial basis function method; Extrapolation relationships for all logging data are established using the Kriging interpolation method, enabling the extrapolation and filling of logging velocities along structural and related facies zones. The classification mask M is divided into a 0-value indication area and a 1-value indication area; The 0-value indication area represents the unfilled area of ​​the logging velocity interpolation filling model P; The 1-value indication area represents the area filled with velocity values ​​after weighted fusion of the seismic velocity model and the well logging velocity interpolation filling model using the weighted least squares method.

Images