Depth domain time-lapse seismic joint inversion method and device, equipment and storage medium

Abstract

The invention relates to a depth domain time-lapse seismic joint inversion method and device, equipment and a storage medium, and the method comprises the steps: extracting depth domain seismic wavelets according to obtained logging data, and building a depth domain seismic wavelet matrix; based on the depth domain seismic wavelet matrix, well-seismic calibration is carried out on the basic logging data and the monitoring logging data; solving a difference between the corrected basic logging data and the corrected monitoring logging data, determining differential logging data, performing full-band interpolation according to the differential logging data and the corrected basic logging data to construct an initial model, and obtaining a low-frequency model according to the initial model; establishing a time-lapse seismic simultaneous inversion objective function based on the forward operator of the time-lapse seismic simultaneous inversion; and substituting the low-frequency model into the target function for iterative solution, and outputting a time-lapse seismic inversion result meeting an iteration ending condition, thereby overcoming the limitation of the existing time-lapse seismic method while reducing the workload of time-domain time-lapse seismic inversion processing.

Description

Technical Field

[0001] The present disclosure relates to the technical field of seismic time-lapse reservoir monitoring, and particularly to a depth-domain seismic time-lapse joint inversion method, device, equipment and storage medium. Background Art

[0002] Seismic time-lapse monitoring technology is an effective reservoir dynamic monitoring technology, which has been maturely applied in multiple fields such as oilfield water injection development and CO2 sequestration monitoring. The results of time-lapse seismic interpretation can be used to locate the development of the reservoir, and at the same time, the remaining oil range can be predicted to deploy new well positions. Its principle is mainly that the changes in reservoir physical parameters (porosity, permeability, water saturation, etc.) caused by oil and gas production will produce obvious amplitude differences in time-lapse seismic data. Based on seismic inversion theory, this amplitude change can be used to obtain the corresponding reservoir parameter changes.

[0003] Time-lapse seismic inversion methods can be divided into separate inversion, differential inversion and simultaneous inversion. Time-lapse seismic separate inversion is to first invert two periods of seismic data separately, and then take the difference between the two inversion results to obtain reservoir changes. Since two data inversions are required, the differential parameters calculated by this method usually contain large calculation errors. In contrast, differential inversion is to directly calculate seismic differential parameters using the difference between two periods of seismic data, reducing the intermediate calculation process and thus reducing the calculation amount and calculation error. However, this method needs to assume that the two-period seismic wavelets are consistent and the time-lapse difference is small during calculation. Time-lapse seismic simultaneous inversion also uses two periods of seismic data at the same time, and can simultaneously invert differential parameters and basic data elastic parameters, and there is no assumption of wavelet and time-lapse difference, which has a wider applicability.

[0004] Currently, the processing and interpretation of time-lapse seismic, whether it is pre-stack or post-stack inversion, are basically realized in the time domain. However, due to the changes in underground reservoir parameters caused by development, which lead to changes in seismic elastic parameters (mainly seismic P-wave velocity), the time-lapse two-period seismic data in the time domain will produce a non-linear displacement amount with time at the reservoir position. The conventional approach is to first correct this displacement amount before performing inversion interpretation. Due to the development of depth migration technology, carrying out depth-domain time-lapse seismic technology can simplify the well-seismic calibration process, and since the change in P-wave velocity caused by reservoir changes in the depth domain will not cause reservoir position changes, the workload of matching time-lapse seismic data can be greatly reduced.

[0005] In summary, the following defects exist in the current research on time-lapse seismic inversion methods: First, the time-lapse seismic separate inversion method requires two seismic inversions, with a large calculation amount and may lead to large calculation errors; Second, time-domain time-lapse seismic processing and interpretation cannot avoid multiple well-seismic calibration problems and time-lapse seismic data matching problems, which may introduce human errors. Summary of the Invention

[0006] To solve the above technical problems or at least partially solve the above technical problems, embodiments of the present disclosure provide a deep-domain time-lapse seismic joint inversion method, apparatus, device, and storage medium.

[0007] In a first aspect, embodiments of the present disclosure provide a deep-domain time-lapse seismic joint inversion method, the method comprising:

[0008] Obtain logging data, based on the point spread function, extract a deep-domain seismic wavelet from the logging data, and establish a corresponding deep-domain seismic wavelet matrix;

[0009] Based on the deep-domain seismic wavelet matrix, perform well-seismic calibration on the basic logging data and the monitoring logging data to obtain the corrected basic logging data and the monitoring logging data;

[0010] Take the difference between the corrected basic logging data and the monitoring logging data to determine the differential logging data, construct an initial model by full-band interpolation according to the seismic horizon and logging horizon information, as well as the differential logging data and the corrected basic logging data, and perform low-pass filtering on the constructed initial model to obtain a low-frequency model, wherein the low-frequency model includes basic data elastic parameters and differential elastic parameters;

[0011] Determine the forward operator for time-lapse seismic simultaneous inversion, and establish a time-lapse seismic simultaneous inversion objective function based on this forward operator, wherein the first term of the objective function is the error term between the forward simulation record and the actual seismic data, and the second term is the low-frequency constraint term;

[0012] Substitute the low-frequency model into the time-lapse seismic simultaneous inversion objective function, and perform iterative solution of the objective function. Take the error of the first term being lower than a preset threshold as the iteration end condition, and output the deep-domain basic data elastic parameters and differential elastic parameters that meet the iteration end condition as the time-lapse seismic inversion result.

[0013] In a possible implementation manner, the method further comprises:

[0014] Perform an inversion test on the seismic trace near the well according to the time-lapse seismic inversion result, compare the logging results with the inversion results of the basic data elastic parameters and the differential elastic parameters respectively, and adjust at least one of the iteration times and the weight factor according to the comparison results;

[0015] Perform parallel time-lapse seismic inversion on the seismic data of the work area according to the adjusted iteration times and the weight factor to obtain the deep-domain basic data elastic parameters and differential elastic parameters of the work area.

[0016] In a possible implementation, based on the point spread function, extracting the depth-domain seismic wavelet from well logging data and establishing a corresponding depth-domain seismic wavelet matrix includes:

[0017] Extracting the time-domain seismic wavelet from well logging data by angle, and based on the point spread function, determining the corresponding depth-domain seismic wavelet according to the well logging data and the time-domain seismic wavelet, and constructing a depth-domain wavelet matrix.

[0018] In a possible implementation, based on the point spread function, determining the corresponding depth-domain seismic wavelet according to the well logging data and the time-domain seismic wavelet, and constructing a depth-domain wavelet matrix includes:

[0019] Under a one-dimensional velocity model, the following relationship is satisfied between the point propagation function and the time-domain seismic wavelet:

[0020]

[0021] Where T represents the period of the time-domain wavelet, which is used to represent the integral of the slowness of the spatial point spread function within one wavelength, λ is the wavelength, h is the depth, and v is the velocity corresponding to this depth;

[0022] Assuming that the wavelet is a zero-phase wavelet, when calculating the depth-domain wavelet at the velocity interface position, calculate the depth-domain wavelets corresponding to the upper and lower layer velocities respectively, intercept at the maximum value of the wavelet, and then splice and interpolate to calculate the depth-domain wavelet at this depth;

[0023] Based on the well logging data, construct a Toeplitz matrix for the depth-domain wavelets at each depth to obtain a depth-domain wavelet matrix.

[0024] In a possible implementation, based on the depth-domain seismic wavelet matrix, performing well-seismic calibration on the basic well logging data and the monitoring well logging data to obtain the corrected basic well logging data and monitoring well logging data includes:

[0025] Determine the reflection coefficient sequence according to the well logging data, and forward model the depth-domain pre-stack angle gathers based on the depth-domain seismic wavelet matrix and the reflection coefficient sequence, and perform well-seismic calibration on the depth-domain pre-stack angle gathers, the monitoring well logging data and the basic well logging data to obtain the corrected basic well logging data and monitoring well logging data.

[0026] In a possible implementation, determining the forward operator for time-lapse seismic simultaneous inversion and establishing a time-lapse seismic simultaneous inversion objective function based on this forward operator includes:

[0027] Linearize the non - linear Zoeppritz forward operator at the low - frequency model. Based on the linear equations, establish the inversion objective functions for the base data and the monitoring data respectively, and then take the derivatives of the two inversion objective functions and combine them to obtain the time - lapse seismic simultaneous inversion objective function.

[0028] In a possible implementation, the process of linearizing the non - linear Zoeppritz forward operator at the low - frequency model, establishing the inversion objective functions for the base data and the monitoring data respectively based on the linear equations, taking the derivatives of the two inversion objective functions and combining them to obtain the time - lapse seismic simultaneous inversion objective function includes:

[0029] The expressions for the base data and the monitoring data are as follows:

[0030] d1 = G(m1)=R(m1)*W1(m1)

[0031] d2 = G(m2)=R(m2)*W2(m2)

[0032] Where, G represents the non - linear operator of the Zoeppritz forward modeling; m represents the seismic elastic parameters; R represents the reflection coefficient sequence of the Zoeppritz forward modeling; W represents the wavelet matrix in the depth domain; d1 and d2 represent the base data and the monitoring data respectively.

[0033] The linearization of the non - linear operator adopts Taylor expansion. The expansion expression after expanding at the low - frequency model and retaining the linear part is:

[0034]

[0035]

[0036] Where, m0 and m'0 represent the low - frequency models of the elastic parameters of the base data and the monitoring data respectively.

[0037] Based on the expansion expression, establish the inversion equations for the base data and the monitoring data respectively:

[0038]

[0039]

[0040] Among them, the first term is the error term, and a, b, a', b' are the corresponding gradient terms and intercept terms in equations (4) - (5) respectively. The specific forms are:

[0041]

[0042]

[0043]

[0044]

[0045] The second term is the low-frequency constraint term, and the third term is the smoothing term. Here, ε and η respectively represent the weight factors of the constraint terms. Here, L represents the difference operator, and its expression is:

[0046]

[0047] Here, nt represents the number of vertical sampling points of the data.

[0048] Derive the inversion equations for the basic data and the monitoring data respectively, and set the derivatives equal to 0, the following forms can be obtained:

[0049]

[0050]

[0051] Rearrange and transpose the above expressions, and assume m2 = m1 + Δm, the simultaneous inversion equation for time-lapse seismic is obtained:

[0052]

[0053] Here, G b , G m , D b , D m respectively have the following forms:

[0054] G b = a T a + ε 2 + η 2 L T L, D b = a T (d1 - b) + ε 2 m0

[0055] G m = a' T a' + ε 2 + η 2 L T L, D m = a' T (d2 - b') + ε 2 m'0

[0056] The objective function for simultaneous inversion of time-lapse seismic in the depth domain is:

[0057]

[0058] Here, D newFor the sequence on the right side of the equal sign in the simultaneous inversion equation of time-lapse seismic, G new For the matrix on the left side of the equal sign in the simultaneous inversion equation of time-lapse seismic, m new For the matrix on the left side of the equal sign in the simultaneous inversion equation of time-lapse seismic, and γ represents the weight factor of the low-frequency constraint term.

[0059] In a second aspect, embodiments of the present disclosure provide a deep-domain time-lapse seismic joint inversion device, including:

[0060] An extraction module, configured to obtain logging data, extract a deep-domain seismic wavelet based on the point spread function according to the logging data, and establish a corresponding deep-domain seismic wavelet matrix;

[0061] A calibration module, configured to perform well-seismic calibration on the basic logging data and the monitoring logging data based on the deep-domain seismic wavelet matrix to obtain the corrected basic logging data and monitoring logging data;

[0062] A construction module, configured to find the difference between the corrected basic logging data and the monitoring logging data to determine the differential logging data, perform full-band interpolation construction on the initial model according to the seismic horizon and logging horizon information, the differential logging data, and the corrected basic logging data, and perform low-pass filtering on the constructed initial model to obtain a low-frequency model, where the low-frequency model includes basic data elastic parameters and differential elastic parameters;

[0063] An establishment module, configured to determine the forward operator of the simultaneous inversion of time-lapse seismic and establish a simultaneous inversion objective function of time-lapse seismic based on the forward operator, where the first term of the objective function is the error term between the forward simulation record and the actual seismic data, and the second term is the low-frequency constraint term;

[0064] A solution module, configured to substitute the low-frequency model into the simultaneous inversion objective function of time-lapse seismic, perform iterative solution of the objective function, use the error of the first term being lower than a preset threshold as the iteration end condition, and output the deep-domain basic data elastic parameters and differential elastic parameters that meet the iteration end condition as the time-lapse seismic inversion result.

[0065] In a third aspect, embodiments of the present disclosure provide an electronic device, including a processor, a communication interface, a memory, and a communication bus, where the processor, the communication interface, and the memory complete mutual communication through the communication bus;

[0066] The memory is used to store a computer program;

[0067] The processor is configured to implement the above-mentioned deep-domain time-lapse seismic joint inversion method when executing the program stored in the memory.

[0068] Fourthly, an embodiment of the present disclosure provides a computer-readable storage medium, on which a computer program is stored. The computer program, when executed by a processor, implements the above-mentioned depth-domain time-lapse seismic joint inversion method.

[0069] The above technical solutions provided by the embodiments of the present disclosure have at least some or all of the following advantages compared with the prior art:

[0070] For the depth-domain time-lapse seismic joint inversion method described in the embodiments of the present disclosure, logging data is acquired, and based on the point spread function, a depth-domain seismic wavelet is extracted from the logging data, and a corresponding depth-domain seismic wavelet matrix is established; based on the depth-domain seismic wavelet matrix, well-seismic calibration is performed on the basic logging data and the monitoring logging data to obtain the corrected basic logging data and the monitoring logging data; the difference between the corrected basic logging data and the monitoring logging data is calculated to determine the differential logging data, and an initial model is constructed by full-band interpolation according to the seismic horizon and logging horizon information, as well as the differential logging data and the corrected basic logging data, and the constructed initial model is subjected to low-pass filtering to obtain a low-frequency model, where the low-frequency model includes basic data elastic parameters and differential elastic parameters; a forward operator for time-lapse seismic simultaneous inversion is determined, and a time-lapse seismic simultaneous inversion objective function is established based on this forward operator, where the first term of the objective function is an error term between the forward simulation record and the actual seismic data, and the second term is a low-frequency constraint term; the low-frequency model is substituted into the time-lapse seismic simultaneous inversion objective function, and the objective function is iteratively solved, and the error of the first term being lower than a preset threshold is used as the iteration end condition, and the depth-domain basic data elastic parameters and differential elastic parameters that meet the iteration end condition are output as the time-lapse seismic inversion result. By performing depth-domain wavelet calculation, depth-domain non-steady-state convolution forward simulation, establishment of a low-frequency model, establishment of a time-lapse seismic simultaneous inversion objective function, and simultaneous solution of basic data elastic parameters and differential elastic parameters, while reducing the workload of time-domain time-lapse seismic inversion processing, the limitations of existing time-lapse seismic methods are overcome, and it can better indicate the oil reservoir development situation and the deployment arrangement of new well positions, and has great potential in the actual application field of oil reservoir reservoir characterization. BRIEF DESCRIPTION OF THE DRAWINGS

[0071] The accompanying drawings here are incorporated into the specification and form a part of this specification, showing embodiments consistent with the present disclosure, and are used together with the specification to explain the principles of the present disclosure.

[0072] In order to more clearly illustrate the technical solutions in the embodiments of the present disclosure or the prior art, the following will briefly introduce the accompanying drawings required for use in the description of the embodiments or related technologies. Obviously, for those of ordinary skill in the art, other drawings can also be obtained based on these drawings without creative efforts.

[0073] Figure 1 Schematically shows a schematic flow chart of a depth-domain time-lapse seismic joint inversion method according to an embodiment of the present disclosure;

[0074] Figure 2 Schematically shows a schematic diagram of a logging curve according to an embodiment of the present disclosure;

[0075] Figure 3 Schematically shows a schematic diagram of a differential elastic parameter curve obtained by calculating the difference between two periods of logging data according to an embodiment of the present disclosure;

[0076] Figure 4 Schematically shows a schematic diagram of a depth-domain wavelet matrix corresponding to basic logging data according to an embodiment of the present disclosure;

[0077] Figure 5 Schematically shows a schematic diagram of a depth-domain wavelet matrix corresponding to monitoring logging data according to an embodiment of the present disclosure;

[0078] Figure 6 Schematically shows a schematic diagram of the inversion result of basic data in noiseless depth-domain time-lapse seismic simultaneous inversion of logging data according to an embodiment of the present disclosure;

[0079] Figure 7 Schematically shows a schematic diagram of the inversion result of differential elastic parameters in noiseless condition according to an embodiment of the present disclosure;

[0080] Figure 8 Schematically shows a schematic diagram of the inversion result of basic data of logging data after adding random noise with a signal-to-noise ratio of 2 to the synthetic seismic record according to an embodiment of the present disclosure;

[0081] Figure 9 Schematically shows a schematic diagram of the inversion result of differential elastic parameters of logging data after adding random noise with a signal-to-noise ratio of 2 to the synthetic seismic record according to an embodiment of the present disclosure;

[0082] Figure 10 Schematically shows a block diagram of the structure of a depth-domain time-lapse seismic joint inversion device according to an embodiment of the present disclosure;

[0083] Figure 11 Schematically shows a block diagram of the structure of an electronic device according to an embodiment of the present disclosure. Detailed implementation manners

[0084] To make the objectives, technical solutions, and advantages of the embodiments of the present disclosure clearer, the technical solutions in the embodiments of the present disclosure will be clearly and completely described below with reference to the accompanying drawings in the embodiments of the present disclosure. Apparently, the described embodiments are some, but not all, of the embodiments of the present disclosure. All other embodiments obtained by those of ordinary skill in the art based on the embodiments in the present disclosure without creative efforts shall fall within the scope of protection of the present disclosure.

[0085] See Figure 1 , embodiments of the present disclosure provide a deep-time-lapse seismic joint inversion method, and the method includes:

[0086] S1. Obtain well logging data, extract a deep-domain seismic wavelet based on the point spread function according to the well logging data, and establish a corresponding deep-domain seismic wavelet matrix.

[0087] In this embodiment, the well logging data includes a longitudinal wave curve, a transverse wave curve, and a density curve.

[0088] S2. Based on the deep-domain seismic wavelet matrix, perform well-seismic calibration on the basic well logging data and the monitoring well logging data to obtain the corrected basic well logging data and the monitoring well logging data.

[0089] In this embodiment, the well-seismic calibration is to perform a small-range well-seismic calibration according to the seismic horizon information and the well logging horizon information.

[0090] S3. Take the difference between the corrected basic well logging data and the monitoring well logging data to determine the differential well logging data. Construct an initial model by full-band interpolation according to the seismic horizon and well logging horizon information, the differential well logging data, and the corrected basic well logging data, and perform low-pass filtering on the constructed initial model to obtain a low-frequency model, where the low-frequency model includes basic data elastic parameters and differential elastic parameters.

[0091] S4. Determine the forward operator for time-lapse seismic simultaneous inversion, and establish a time-lapse seismic simultaneous inversion objective function based on this forward operator, where the first term of the objective function is an error term between the forward simulation record and the actual seismic data, and the second term is a low-frequency constraint term.

[0092] S5. Substitute the low-frequency model into the time-lapse seismic simultaneous inversion objective function, and perform iterative solution of the objective function. Use the error of the first term being lower than a preset threshold as the iteration end condition, and output the deep-domain basic data elastic parameters and differential elastic parameters that meet the iteration end condition as the time-lapse seismic inversion result.

[0093] In this embodiment, the iterative solution of the objective function is implemented based on the Gauss-Newton algorithm.

[0094] In this embodiment, the method further includes:

[0095] Perform a wellside seismic trace inversion test based on the time-lapse seismic inversion results, compare the logging results with the elastic parameter inversion results of the basic data and the differential elastic parameter inversion results respectively, and adjust at least one of the iteration times and the weight factor according to the comparison results;

[0096] Perform parallel time-lapse seismic inversion on the seismic data in the work area according to the adjusted iteration times and weight factor to obtain the elastic parameters of the basic data and the differential elastic parameters in the depth domain of the work area.

[0097] In this embodiment, in step S1, based on the point spread function, extracting the seismic wavelet in the depth domain from the logging data and establishing the corresponding seismic wavelet matrix in the depth domain includes:

[0098] Extract the seismic wavelet in the time domain from the logging data by angle, and based on the point spread function, determine the corresponding seismic wavelet in the depth domain according to the logging data and the seismic wavelet in the time domain, and construct a depth domain wavelet matrix.

[0099] In this embodiment, based on the point spread function, determining the corresponding seismic wavelet in the depth domain according to the logging data and the seismic wavelet in the time domain, and constructing a depth domain wavelet matrix includes:

[0100] In a one-dimensional velocity model, the following relationship is satisfied between the point propagation function and the seismic wavelet in the time domain:

[0101]

[0102] where T represents the period of the seismic wavelet in the time domain, which is used to represent the integral of the slowness of the spatial point spread function within one wavelength, λ is the wavelength, h is the depth, and v is the velocity corresponding to this depth;

[0103] Since the underground velocity is non-uniform, assuming that the wavelet is a zero-phase wavelet, when calculating the seismic wavelet in the depth domain at the velocity interface position, calculate the seismic wavelets in the depth domain corresponding to the upper and lower layer velocities respectively, intercept at the maximum value of the wavelet, and then splice and interpolate to calculate the seismic wavelet in the depth domain at this depth;

[0104] Construct a Toeplitz matrix based on the seismic wavelets in the depth domain at each depth from the logging data to obtain a seismic wavelet matrix in the depth domain.

[0105] In this embodiment, in step S2, based on the seismic wavelet matrix in the depth domain, performing well-seismic calibration on the basic logging data and the monitoring logging data to obtain the corrected basic logging data and monitoring logging data includes:

[0106] Determine the reflection coefficient sequence according to the well logging data, forward model the pre-stack angle gathers in the depth domain based on the seismic wavelet matrix in the depth domain and the reflection coefficient sequence, and perform well-seismic calibration on the pre-stack angle gathers in the depth domain, the monitored well logging data, and the basic well logging data to obtain the corrected basic well logging data and monitored well logging data.

[0107] In this embodiment, the determining the reflection coefficient sequence according to the well logging data includes:

[0108] Substitute the well logging data into the Zoeppritz equation to calculate the reflection coefficient sequence.

[0109] In this embodiment, the forward modeling of the pre-stack angle gathers in the depth domain based on the seismic wavelet matrix in the depth domain and the reflection coefficient sequence is implemented based on the non-steady-state convolution algorithm.

[0110] In this embodiment, in step S4, the determining the forward operator of the time-lapse seismic simultaneous inversion and establishing the time-lapse seismic simultaneous inversion objective function based on the forward operator includes:

[0111] Linearize the non-linear Zoeppritz forward operator at the low-frequency model, establish the inversion objective functions for the basic data and the monitored data respectively based on the linear equations, and take the derivatives of the two inversion objective functions and combine them to obtain the time-lapse seismic simultaneous inversion objective function.

[0112] In this embodiment, the linearizing the non-linear Zoeppritz forward operator at the low-frequency model, establishing the inversion objective functions for the basic data and the monitored data respectively based on the linear equations, and taking the derivatives of the two inversion objective functions and combining them to obtain the time-lapse seismic simultaneous inversion objective function includes:

[0113] The expressions for the basic data and the monitored data are as follows:

[0114] d1 = G(m1) = R(m1) * W1(m1)

[0115] d2 = G(m2) = R(m2) * W2(m2)

[0116] Where, G represents the non-linear operator of the Zoeppritz forward modeling; m represents the seismic elastic parameter; R represents the reflection coefficient sequence of the Zoeppritz forward modeling; W represents the wavelet matrix in the depth domain; d1 and d2 represent the basic data and the monitored data respectively.

[0117] The linearization of the non-linear operator uses the Taylor expansion. The expansion expression after expanding at the low-frequency model and retaining the linear part is:

[0118]

[0119]

[0120] Among them, m0 and m'0 respectively represent the low-frequency models of the elastic parameters of the basic data and the monitoring data.

[0121] Based on the expansion expressions, the inversion equations for the basic data and the monitoring data are established respectively:

[0122]

[0123]

[0124] Among them, the first term is the error term, and a, b, a', and b' are the corresponding gradient terms and intercept terms in the expansion expression, and their expressions are:

[0125]

[0126]

[0127]

[0128]

[0129] The second term is the low-frequency constraint term, and the third term is the smoothing term. Among them, ε and η respectively represent the weight factors of the constraint term; among them, L represents the difference operator, and its expression is:

[0130]

[0131] Among them, nt represents the number of longitudinal sampling points of the data.

[0132] Derive the inversion equations for the basic data and the monitoring data respectively, and set the derivative equal to 0 to obtain the following expressions:

[0133]

[0134]

[0135] Rearrange and transpose the above expressions, and assume m2 = m1 + Δm to obtain the time-lapse seismic simultaneous inversion equation:

[0136]

[0137] Among them, G b , G m , D b , D m The expressions are as follows:

[0138] G b = a T a + ε 2 + η2 L T L, D b = a T (d1 - b)+ε 2 m0

[0139] G m = a' T a'+ε 2 +η 2 L T L, D m = a' T (d2 - b')+ε 2 m'0

[0140] The objective function for simultaneous inversion of time - lapse seismic in the depth domain is as follows:

[0141]

[0142] Where D new is the sequence on the right - hand side of the equal sign of the simultaneous inversion equation for time - lapse seismic, G new is the matrix on the left - hand side of the equal sign of the simultaneous inversion equation for time - lapse seismic, m new is the matrix on the left - hand side of the equal sign of the simultaneous inversion equation for time - lapse seismic, and γ represents the weight factor of the low - frequency constraint term.

[0143] In this embodiment, in step S5, the low - frequency model is substituted into the objective function for simultaneous inversion of time - lapse seismic through the following expression, and the objective function is iteratively solved for m new :

[0144]

[0145] Where k is the number of iterations, α is the iterative update step size, usually taking the value of 1, is the derivative of the objective function with respect to the model parameters, and its value is equal to G new , and the iteration stops when the error value is less than the given error threshold, and the current iteration result is output simultaneously

[0146] Through the joint inversion method of time - lapse seismic in the depth domain of the present disclosure, the time - lapse log curves shown in Figure 2 are inverted, and the specific process is as follows:

[0147] Substitute the Figure 2 shown time - lapse log curves into the point - spread function solution formula, and the depth wavelet matrices corresponding to the two - phase seismic data shown in Figure 4 and Figure 5 can be calculated. Figure 3 is the differential elastic parameter curve obtained by taking the difference between the two - phase log data. Where Figure 4 and Figure 5From left to right, they respectively represent the P-wave velocity, S-wave velocity, and density. The solid line represents the three basic logging parameter curves, and the dashed line represents the three monitored logging parameter curves. Figure 4 is the wavelet matrix in the depth domain corresponding to the basic logging data, Figure 5 is the wavelet matrix in the depth domain corresponding to the monitored logging data. As can be seen from the local magnification of the box in Figure 4 and Figure 5 , different from the seismic wavelet matrix in the time domain, due to the variation of the seismic wave velocity in the vertical direction, the seismic wavelet in the depth domain is distorted at the velocity interface, showing asymmetry. Figure 6 is the inversion result of the basic data when there is no noise in the depth-domain time-lapse seismic simultaneous inversion of logging data, Figure 7 is the inversion result of the differential elastic parameters when there is no noise. Among them, Figure 6 from left to right are the P-wave velocity (V p ) of the basic data, the S-wave velocity (V s ), and the density (ρ). Figure 6 from left to right are the differential P-wave velocity (ΔV p ) of the differential elastic parameters, the differential S-wave velocity (ΔV s ), and the differential density (Δρ). The curves shown in 601 are the true logging data and the differential logging data, and the curves shown in 602 are the inverted logging data and the differential logging data. As shown in Figure 8 and Figure 9 are the inversion results of the logging data after adding random noise with a signal-to-noise ratio of 2 to the synthetic seismic record. Among them, Figure 8 from left to right are the P-wave velocity (V p ) of the basic data, the S-wave velocity (V s ), and the density (ρ) when the signal-to-noise ratio is 2. Figure 9 from left to right are the differential P-wave velocity (ΔV p ) of the differential elastic parameters, the differential S-wave velocity (ΔV s ), and the differential density (Δρ) when the signal-to-noise ratio is 2. The curves shown in 701 are the true logging data and the differential logging data, and the curves shown in 702 are the inverted logging data and the differential logging data.

[0148] It can be seen from the above inversion test of the logging data that the depth-domain time-lapse seismic joint inversion method of the present disclosure can comprehensively and effectively utilize the two-phase data to calculate the elastic parameters and differential elastic parameters of the basic data. It can be seen from the anti-noise test that this method has a certain robustness and can be applied to seismic data inversion. And because the present disclosure overcomes the limitations of seismic wavelets and differential parameters, it can meet the more extensive time-lapse seismic reservoir prediction requirements while satisfying the inversion accuracy.

[0149] The joint inversion method of time-lapse seismic in the depth domain according to the present disclosure utilizes well logging data, calculates a seismic wavelet matrix in the depth domain based on the point spread function and the longitudinal wave travel-time curve of the well logging data, calculates the forward operator of the simultaneous inversion of time-lapse seismic based on the Zoeppritz equation, forward simulates the pre-stack angle gathers in the depth domain based on the non-steady convolution theory, performs well-seismic calibration in the depth domain, interpolates the well logging data and the differential well logging data in the depth domain in the full frequency band based on the horizon and well logging information, and obtains a low-frequency model through low-pass filtering processing, establishes the objective equation of the simultaneous inversion of time-lapse seismic and substitutes the low-frequency model into the forward operator of the simultaneous inversion of time-lapse seismic, and solves the objective function through an iterative method to obtain the elastic parameters of the basic data and the differential elastic parameters in the depth domain, which can be used for high-precision reservoir parameter estimation of time-lapse seismic in the depth domain.

[0150] See Figure 10 , an embodiment of the present disclosure provides a joint inversion device for time-lapse seismic in the depth domain, including:

[0151] An extraction module 11, configured to obtain well logging data, extract a seismic wavelet in the depth domain according to the well logging data based on the point spread function, and establish a corresponding seismic wavelet matrix in the depth domain;

[0152] A calibration module 12, configured to perform well-seismic calibration on the basic well logging data and the monitoring well logging data based on the seismic wavelet matrix in the depth domain to obtain the corrected basic well logging data and the monitoring well logging data;

[0153] A construction module 13, configured to find the difference between the corrected basic well logging data and the monitoring well logging data, determine the differential well logging data, perform full-frequency band interpolation construction on the initial model according to the seismic horizon and well logging horizon information, the differential well logging data and the corrected basic well logging data, and perform low-pass filtering processing on the constructed initial model to obtain a low-frequency model, where the low-frequency model includes the elastic parameters of the basic data and the differential elastic parameters;

[0154] An establishment module 14, configured to determine the forward operator of the simultaneous inversion of time-lapse seismic, and establish an objective function for the simultaneous inversion of time-lapse seismic based on the forward operator, where the first term of the objective function is an error term between the forward simulation record and the actual seismic data, and the second term is a low-frequency constraint term;

[0155] A solution module 15, configured to substitute the low-frequency model into the objective function of the simultaneous inversion of time-lapse seismic, and perform iterative solution of the objective function, use the error of the first term being lower than a preset threshold as the iteration end condition, and output the elastic parameters of the basic data and the differential elastic parameters in the depth domain that meet the iteration end condition as the time-lapse seismic inversion result.

[0156] The time-lapse seismic joint inversion device of the present disclosure simultaneously inverts through time-lapse seismic, while retaining the accuracy of time-lapse seismic differential inversion, overcomes the limitations of differential inversion on wavelets and differential data, and better meets the needs of practical applications; moreover, it is carried out in the depth domain, avoiding the time-depth conversion process that needs to be faced in the time domain, simplifying the well-seismic calibration process, and reducing the workload of time-lapse seismic data matching; compared with linear approximation, it overcomes the problem of reduced accuracy of large-angle incidence and can be applied to time-lapse seismic inversion in any angle range.

[0157] For the specific implementation process of the functions and roles of each unit in the above device, please refer to the implementation process of the corresponding steps in the above method, which will not be elaborated here.

[0158] For the device embodiment, since it basically corresponds to the method embodiment, the relevant parts can be referred to the partial description of the method embodiment. The device embodiments described above are only illustrative. The units described as separate components may or may not be physically separated, and the components shown as units may or may not be physical units, that is, they may be located in one place or distributed to multiple network units. Some or all of the modules can be selected according to actual needs to achieve the purpose of the solution of the present invention. Those of ordinary skill in the art can understand and implement it without creative efforts.

[0159] In the above embodiments, any combination of the extraction module 11, the calibration module 12, the construction module 13, the establishment module 14, and the solution module 15 can be combined and implemented in one module, or any one of the modules can be split into multiple modules. Or, at least part of the functions of one or more of these modules can be combined with at least part of the functions of other modules and implemented in one module. At least one of the extraction module 11, the calibration module 12, the construction module 13, the establishment module 14, and the solution module 15 can be at least partially implemented as a hardware circuit, such as a field programmable gate array (FPGA), a programmable logic array (PLA), a system on chip, a system on substrate, a system on package, an application specific integrated circuit (ASIC), or can be implemented by any other reasonable way of integrating or packaging circuits, etc., in hardware or firmware, or implemented in any one of the three implementation ways of software, hardware, and firmware, or in any appropriate combination of several of them. Or, at least one of the extraction module 11, the calibration module 12, the construction module 13, the establishment module 14, and the solution module 15 can be at least partially implemented as a computer program module, and when the computer program module is run, the corresponding functions can be executed.

[0160] Refer to Figure 11As shown in the figure, the electronic device provided by the embodiments of the present disclosure includes a processor 1110, a communication interface 1120, a memory 1130, and a communication bus 1140. Among them, the processor 1110, the communication interface 1120, and the memory 1130 complete mutual communication through the communication bus 1140;

[0161] The memory 1130 is used to store computer programs;

[0162] When the processor 1110 is used to execute the program stored on the memory 1130, the following deep domain time-lapse seismic joint inversion method is implemented:

[0163] Obtain well logging data, based on the point spread function, extract the deep domain seismic wavelet according to the well logging data, and establish a corresponding deep domain seismic wavelet matrix;

[0164] Based on the deep domain seismic wavelet matrix, perform well-seismic calibration on the basic well logging data and the monitoring well logging data to obtain the corrected basic well logging data and the monitoring well logging data;

[0165] Take the difference between the corrected basic well logging data and the monitoring well logging data to determine the differential well logging data. According to the seismic horizon and well logging horizon information, as well as the differential well logging data and the corrected basic well logging data, perform full-band interpolation to construct an initial model, and perform low-pass filtering on the constructed initial model to obtain a low-frequency model, where the low-frequency model includes the basic data elastic parameters and the differential elastic parameters;

[0166] Determine the forward operator for time-lapse seismic simultaneous inversion, and establish a time-lapse seismic simultaneous inversion objective function based on this forward operator. Among them, the first term of the objective function is the error term between the forward simulation record and the actual seismic data, and the second term is the low-frequency constraint term;

[0167] Substitute the low-frequency model into the time-lapse seismic simultaneous inversion objective function, and perform iterative solution of the objective function. Take the error of the first term being lower than the preset threshold as the iteration end condition, and output the deep domain basic data elastic parameters and the differential elastic parameters that meet the iteration end condition as the time-lapse seismic inversion result.

[0168] The above communication bus 1140 can be a Peripheral Component Interconnect (PCI) bus or an Extended Industry Standard Architecture (EISA) bus, etc. This communication bus 1140 can be divided into an address bus, a data bus, a control bus, etc. For the sake of simplicity, only a thick line is used to represent it in the figure, but it does not mean that there is only one bus or one type of bus.

[0169] The communication interface 1120 is used for communication between the above-mentioned electronic device and other devices.

[0170] The memory 1130 may include a random access memory (RAM), or may also include a non-volatile memory, such as at least one disk memory. Optionally, the memory 1130 may also be at least one storage device located far from the aforementioned processor 1110.

[0171] The above-mentioned processor 1110 may be a general-purpose processor, including a central processing unit (CPU), a network processor (NP), etc.; it may also be a digital signal processor (DSP), an application specific integrated circuit (ASIC), a field-programmable gate array (FPGA), or other programmable logic devices, discrete gate or transistor logic devices, discrete hardware components.

[0172] Embodiments of the present disclosure also provide a computer-readable storage medium. A computer program is stored on the above-mentioned computer-readable storage medium, and when the computer program is executed by a processor, the depth-domain time-lapse seismic joint inversion method described above is implemented.

[0173] The computer-readable storage medium may be included in the device / device described in the above embodiments; it may also exist alone without being assembled into the device / device. The above-mentioned computer-readable storage medium carries one or more programs, and when the above one or more programs are executed, the depth-domain time-lapse seismic joint inversion method according to the embodiments of the present disclosure is implemented.

[0174] According to the embodiments of the present disclosure, the computer-readable storage medium may be a non-volatile computer-readable storage medium, for example, it may include but is not limited to: portable computer disks, hard disks, random access memories (RAMs), read-only memories (ROMs), erasable programmable read-only memories (EPROMs or flash memories), portable compact disk read-only memories (CD-ROMs), optical storage devices, magnetic storage devices, or any suitable combination of the above. In the present disclosure, the computer-readable storage medium may be any tangible medium that contains or stores a program, and the program can be used by or combined with an instruction execution system, device, or device.

[0175] It should be noted that in this text, relational terms such as "first" and "second" are only used to distinguish one entity or operation from another entity or operation, and do not necessarily require or imply any actual relationship or order between these entities or operations. Moreover, the term "comprising", "including" or any other variant thereof is intended to cover non-exclusive inclusion, so that a process, method, article or device comprising a series of elements not only includes those elements, but also includes other elements not expressly listed, or further includes elements inherent to such process, method, article or device. Without further limitation, an element defined by the statement "comprising an..." does not exclude the existence of additional identical elements in the process, method, article or device comprising the said element.

[0176] The above are only specific embodiments of the present disclosure, enabling those skilled in the art to understand or implement the present disclosure. Various modifications to these embodiments will be obvious to those skilled in the art, and the general principles defined herein can be implemented in other embodiments without departing from the spirit or scope of the present disclosure. Therefore, the present disclosure will not be limited to these embodiments shown herein, but rather will be accorded the widest scope consistent with the principles and novel features claimed herein.

Claims

1. A depth-domain time-lapse seismic joint inversion method, characterized in that, The method includes: Obtain logging data, extract a depth-domain seismic wavelet based on the point spread function according to the logging data, and establish a corresponding depth-domain seismic wavelet matrix; Based on the depth-domain seismic wavelet matrix, perform well-seismic calibration on the basic logging data and the monitoring logging data to obtain the corrected basic logging data and monitoring logging data; Take the difference between the corrected basic logging data and the monitoring logging data to determine the differential logging data. Construct an initial model through full-band interpolation according to the seismic horizon and logging horizon information, as well as the differential logging data and the corrected basic logging data, and perform low-pass filtering on the constructed initial model to obtain a low-frequency model, where the low-frequency model includes basic data elastic parameters and differential elastic parameters; Determine the forward operator for time-lapse seismic simultaneous inversion, and establish a time-lapse seismic simultaneous inversion objective function based on this forward operator, where the first term of the objective function is the error term between the forward simulation record and the actual seismic data, and the second term is the low-frequency constraint term; Substitute the low-frequency model into the time-lapse seismic simultaneous inversion objective function, and perform iterative solution of the objective function. Take the error of the first term being lower than a preset threshold as the iteration end condition, and output the depth-domain basic data elastic parameters and differential elastic parameters that meet the iteration end condition as the time-lapse seismic inversion result.

2. The method according to claim 1, characterized in that, The method further includes: Perform well-side seismic trace inversion testing according to the time-lapse seismic inversion result, compare the logging results with the inversion results of the basic data elastic parameters and the differential elastic parameters respectively, and adjust at least one of the iteration times and the weight factor according to the comparison results; Perform parallel time-lapse seismic inversion on the seismic data in the work area according to the adjusted iteration times and weight factor to obtain the depth-domain basic data elastic parameters and differential elastic parameters of the work area.

3. The method according to claim 1, characterized in that, The step of extracting a depth-domain seismic wavelet based on the point spread function according to the logging data and establishing a corresponding depth-domain seismic wavelet matrix includes: Extract the time-domain seismic wavelet from the logging data by angle, and based on the point spread function, determine the corresponding depth-domain seismic wavelet according to the logging data and the time-domain seismic wavelet, and construct a depth-domain wavelet matrix.

4. The method according to claim 3, characterized in that, The step of determining the corresponding depth-domain seismic wavelet according to the logging data and the time-domain seismic wavelet based on the point spread function and constructing a depth-domain wavelet matrix includes: Under a one-dimensional velocity model, the following relationship is satisfied between the point propagation function and the time-domain seismic wavelet: where T represents the time-domain wavelet period, which is used to represent the integral of the slowness of the spatial point spread function within one wavelength, λ is the wavelength, h is the depth, and v is the velocity corresponding to this depth; Assume that the wavelet is a zero-phase wavelet. When calculating the depth-domain wavelet at the velocity interface position, calculate the depth-domain wavelets corresponding to the upper and lower layer velocities respectively, intercept at the maximum value of the wavelet, and then splice and interpolate to calculate the depth-domain wavelet at this depth; Based on the logging data, construct a Toeplitz matrix for the depth-domain wavelets at each depth to obtain a depth-domain wavelet matrix.

5. The method according to claim 1, characterized in that, The step of performing well-seismic calibration on the basic logging data and the monitoring logging data based on the depth-domain seismic wavelet matrix to obtain the corrected basic logging data and monitoring logging data includes: A reflection coefficient sequence is determined according to the logging data, and a depth domain prestack angle gather is forward modeled according to the depth domain seismic wavelet matrix and the reflection coefficient sequence. Well-seismic calibration is performed on the depth domain prestack angle gather, the monitoring logging data and the basic logging data to obtain corrected basic logging data and monitoring logging data.

6. The method according to claim 1, characterized in that, The step of determining a forward operator for time-lapse seismic simultaneous inversion and establishing a time-lapse seismic simultaneous inversion objective function based on the forward operator includes: The nonlinear Zoeppritz forward operator is linearized at the low-frequency model. The inversion objective functions of the basic data and monitoring data are established based on linear equations. The two inversion objective functions are derived and combined to obtain the time-lapse seismic simultaneous inversion objective function.

7. The method according to claim 6, characterized in that, The nonlinear Zoeppritz forward operator is linearized at the low-frequency model, and the inversion objective functions of the basic data and the monitoring data are established based on the linear equations, and the two inversion objective functions are derived and combined to obtain the time-lapse seismic simultaneous inversion objective function, including: The expressions of basic data and monitoring data are as follows: d1=G(m1)=R(m1)*W1(m1) d2=G(m2)=R(m2)*W2(m2) Where G represents the nonlinear operator of Zoeppritz forward modeling; m represents the seismic elastic parameter; R represents the reflection coefficient sequence of Zoeppritz forward modeling; W represents the wavelet matrix in the depth domain; d1 and d2 represent the basic data and monitoring data, respectively. The linearization of the nonlinear operator adopts Taylor expansion. The expansion expression after expansion at the low-frequency model and retaining the linear part is: Among them, m0 and m′0 represent the low-frequency models of elastic parameters of basic data and monitoring data, respectively. Based on the expanded expressions, the inversion equations for basic data and monitoring data are established respectively: Among them, the first term is the error term, a, b, a', b' are the corresponding gradient term and intercept term in equations (4) to (5), respectively. The specific form is: The second term is a low-frequency constraint term, and the third term is a smoothing term, where ε and η represent the weight factors of the constraint terms respectively; where L represents the difference operator, and its expression is: Among them, nt represents the number of vertical sampling points of the data, Deriving the inversion equations of basic data and monitoring data respectively and setting the derivatives equal to 0, the following forms can be obtained: Arrange and shift the terms in the above expressions, and assume that m2 = m1 + Δm, and we get the time-lapse seismic simultaneous inversion equation: Among them, G b , G m , D b , D m respectively have the following forms: G b = a T a + ε 2 + η 2 L T L, D b = a T (d1 - b) + ε 2 m0 G m = a' T a' + ε 2 + η 2 L T L, D m = a' T (d2 - b') + ε 2 m'0 The objective function of simultaneous time-lapse seismic inversion in the depth domain is: Among them, D new is the sequence on the right side of the time-lapse seismic simultaneous inversion equation, and G new is the matrix on the left side of the time-lapse seismic simultaneous inversion equation, m new is the matrix on the left side of the time-lapse seismic simultaneous inversion equation, and γ represents the weight factor of the low-frequency constraint term.

8. A device for joint inversion of time-lapse seismic data in the depth domain, characterized in that, include: An extraction module is used to obtain well logging data, extract depth domain seismic wavelets from the well logging data based on a point spread function, and establish a corresponding depth domain seismic wavelet matrix; A calibration module is used to calibrate the basic logging data and monitoring logging data based on the depth domain seismic wavelet matrix to obtain the corrected basic logging data and monitoring logging data; A building block is used to calculate the difference between the corrected basic logging data and the monitoring logging data to determine the differential logging data. Based on the seismic horizon and logging horizon information, as well as the differential logging data and the corrected basic logging data, full-band interpolation is performed to construct an initial model, and the constructed initial model is subjected to low-pass filtering to obtain a low-frequency model. The low-frequency model includes basic data elastic parameters and differential elastic parameters. An establishment module is used to determine the forward operator for time-lapse seismic simultaneous inversion and establish a time-lapse seismic simultaneous inversion objective function based on the forward operator. The first term of the objective function is the error term between the forward simulation record and the actual seismic data, and the second term is the low-frequency constraint term. A solution module is used to substitute the low-frequency model into the time-lapse seismic simultaneous inversion objective function and perform iterative solution of the objective function. Taking the error of the first term being lower than a preset threshold as the iteration end condition, and outputting the depth-domain basic data elastic parameters and differential elastic parameters that meet the iteration end condition as the time-lapse seismic inversion result.

9. An electronic device, characterized in that, It includes a processor, a communication interface, a memory, and a communication bus. Among them, the processor, the communication interface, and the memory complete communication with each other through the communication bus. The memory is used to store computer programs. The processor is used to implement the depth-domain time-lapse seismic joint inversion method according to any one of claims 1-7 when executing the program stored on the memory.

10. A computer-readable storage medium having a computer program stored thereon, characterized in that, The computer program, when executed by the processor, implements the depth-domain time-lapse seismic joint inversion method according to any one of claims 1-7.

Images