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

Abstract

The present disclosure relates to a deep domain time-lapse seismic joint inversion method and device, equipment and storage medium, the method comprises: extracting a deep domain seismic wavelet according to the obtained logging data, establishing a deep domain seismic wavelet matrix; based on the deep domain seismic wavelet matrix, calibrating the basic logging data and the monitoring logging data well-seismic; difference is obtained by subtracting the corrected basic logging data and monitoring logging data, and the difference logging data is determined, the initial model is constructed by full-band interpolation according to the difference logging data and the corrected basic logging data, and the low-frequency model is obtained according to the initial model; the forward operator of time-lapse seismic simultaneous inversion is used to establish the objective function of time-lapse seismic simultaneous inversion; the low-frequency model is substituted into the objective function to solve iteratively, and the time-lapse seismic inversion result meeting the iteration end condition is output, which reduces the workload of time-lapse seismic inversion processing, and overcomes the limitations of the existing time-lapse seismic method.

Description

TECHNICAL FIELD

[0001] The present disclosure relates to the technical field of seismic time-lapse seismic reservoir monitoring, in particular to a deep domain time-lapse seismic joint inversion method and device, equipment and storage medium. BACKGROUND

[0002] Seismic time-lapse seismic monitoring technology is an effective reservoir dynamic monitoring technology, which has been maturely applied in many fields such as oilfield water injection development and CO2 storage monitoring. The time-lapse seismic interpretation result can be used to locate the development situation of the reservoir, and can also be used to predict the range of remaining oil to deploy new well sites. The principle is mainly that the change of reservoir physical parameters (porosity, permeability, water saturation, etc.) caused by oil and gas exploitation will produce obvious amplitude difference on time-lapse seismic data, and based on the seismic inversion theory, the corresponding reservoir parameter change can be obtained by using the amplitude change.

[0003] The time-lapse seismic inversion method can be divided into separate inversion, difference inversion and simultaneous inversion. The time-lapse seismic separate inversion is to first invert two periods of seismic data respectively, and then to calculate the reservoir change by subtracting the two periods of inversion results. Since two data inversions are needed, the difference parameters calculated by using this method usually contain large calculation errors. In contrast, the difference inversion directly calculates the seismic difference parameters by using the difference of two periods of seismic data, which reduces the intermediate calculation process and thus reduces the calculation amount and calculation error. However, this method needs to assume that the two periods of seismic wavelets are consistent and the time-lapse difference is small. The time-lapse seismic simultaneous inversion also uses two periods of seismic data simultaneously, and can simultaneously invert the difference parameters and the basic data elastic parameters, and does not have the wavelet and time-lapse difference assumption, and has more extensive usability.

[0004] At present, the time-lapse seismic processing and interpretation, whether pre-stack or post-stack inversion, are basically realized in the time domain. However, the change of underground reservoir parameters caused by development will cause the change of seismic elastic parameters (mainly seismic P-wave velocity), and the time-lapse two periods of seismic data will produce a time-dependent nonlinear displacement amount at the reservoir position in the time domain. The conventional method needs to correct this displacement amount before inversion and interpretation. Due to the development of depth migration technology, the development of depth domain time-lapse seismic technology can simplify the well-seismic calibration process, and since the P-wave velocity change caused by reservoir change in the depth domain will not cause the change of reservoir position, the matching workload of time-lapse seismic data can be greatly reduced.

[0005] In summary, the current time-lapse seismic inversion method research has the following defects: first, the time-lapse seismic separate inversion method needs to perform two seismic inversions, which has large calculation amount and may cause large calculation error; second, the time domain time-lapse seismic processing and interpretation cannot avoid the problem of multiple well-seismic calibration and the problem of time-lapse seismic data matching, which may introduce artificial error. SUMMARY

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

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

[0008] obtaining logging data, extracting a depth domain seismic wavelet from the logging data based on a point spread function, and establishing a depth domain seismic wavelet matrix corresponding thereto;

[0009] based on the depth domain seismic wavelet matrix, calibrating the base logging data and the monitoring logging data to obtain corrected base logging data and monitoring logging data;

[0010] determining the difference logging data by differencing the corrected base logging data and the monitoring logging data, constructing an initial model by full-band interpolation according to the seismic horizon and logging horizon information and the difference logging data and the corrected base logging data, and performing low-pass filtering processing on the constructed initial model to obtain a low-frequency model, wherein the low-frequency model includes base data elastic parameters and difference elastic parameters;

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

[0012] substituting the low-frequency model into the time-lapse seismic simultaneous inversion objective function, and iteratively solving the objective function, taking the error of the first term being lower than a preset threshold as an iteration end condition, and outputting the depth domain base data elastic parameters and the difference elastic parameters satisfying the iteration end condition as the time-lapse seismic inversion result.

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

[0014] performing wellside seismic trace inversion testing according to the time-lapse seismic inversion result, comparing the logging result with the base data elastic parameter inversion result and the difference elastic parameter inversion result respectively, and adjusting at least one of the iteration number and the weight factor according to the comparison result;

[0015] performing parallel time-lapse seismic inversion on the work area seismic data according to the adjusted iteration number and weight factor to obtain the depth domain base data elastic parameters and the difference elastic parameters of the work area.

[0016] In a possible implementation, the depth domain seismic wavelet is extracted from the logging data based on the point spread function, and a depth domain seismic wavelet matrix corresponding to the depth domain seismic wavelet is established, including:

[0017] The time domain seismic wavelet is extracted from the logging data by angle division, and the corresponding depth domain seismic wavelet is determined based on the point spread function and the logging data and the time domain seismic wavelet, and a depth domain wavelet matrix is constructed.

[0018] In a possible implementation, the depth domain seismic wavelet is determined based on the point spread function and the logging data and the time domain seismic wavelet, and a depth domain wavelet matrix is constructed, including:

[0019] Under a one-dimensional velocity model, the point spread function and the time domain seismic wavelet satisfy the following relationship:

[0020]

[0021] Wherein, T represents a 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 the depth;

[0022] Supposing that the wavelet is a zero-phase wavelet, the depth domain wavelet at the position of the velocity boundary is calculated by respectively calculating the depth domain wavelets corresponding to the upper and lower layer velocities, cutting at the maximum value of the wavelet, and splicing and interpolating to calculate the depth domain wavelet at the depth;

[0023] The depth domain wavelets at each depth are constructed into a Toeplitz matrix based on the logging data, to obtain a depth domain wavelet matrix.

[0024] In a possible implementation, the depth domain seismic wavelet matrix is used to calibrate the basic logging data and the monitoring logging data to obtain corrected basic logging data and monitoring logging data, including:

[0025] The reflection coefficient sequence is determined according to the logging data, and the depth domain prestack angle gather is forward calculated based on the depth domain seismic wavelet matrix and the reflection coefficient sequence, and the depth domain prestack angle gather is calibrated with the monitoring logging data and the basic logging data to obtain corrected basic logging data and monitoring logging data.

[0026] In a possible implementation, the forward operator of the time-lapse seismic simultaneous inversion is determined, and a time-lapse seismic simultaneous inversion objective function is established based on the forward operator, including:

[0027] The nonlinear Zoeppritz forward operator is linearized at a low-frequency model, and inversion objective functions of the basic data and the monitoring data are respectively established based on linear equations, and the inversion objective functions are derived and combined to obtain a simultaneous inversion objective function of time-lapse seismic.

[0028] In a possible implementation, the nonlinear Zoeppritz forward operator is linearized at a low-frequency model, and inversion objective functions of the basic data and the monitoring data are respectively established based on linear equations, and the inversion objective functions are derived and combined to obtain a simultaneous inversion objective function of time-lapse seismic, including:

[0029] Expressions of the basic 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] Wherein, G represents a nonlinear operator of Zoeppritz forward; m represents a seismic elastic parameter; R represents a reflection coefficient sequence of Zoeppritz forward; W represents a depth domain wavelet matrix; d1 and d2 represent the basic data and the monitoring data respectively,

[0033] Linearization of the nonlinear operator adopts Taylor expansion, and an expansion expression after expansion at a low-frequency model and linear part reservation is as follows:

[0034]

[0035]

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

[0037] Based on the expansion expression, inversion equations of the basic data and the monitoring data are respectively established:

[0038]

[0039]

[0040] Wherein, the first term is an error term, a, b, a', b' are gradient terms and intercept terms corresponding to equations (4)-(5) respectively, and specific forms are as follows:

[0041]

[0042]

[0043]

[0044]

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

[0046]

[0047] Wherein nt represents the number of longitudinal sampling points of data,

[0048] The inversion equations of the basic data and the monitoring data are differentiated respectively, and the derivatives are equal to 0, so that the following forms are obtained:

[0049]

[0050]

[0051] The above expressions are rearranged and moved, and it is assumed that m2=m1+Δm, so that the simultaneous inversion equation of time-lapse seismic is obtained:

[0052]

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

[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 of the time-lapse seismic simultaneous inversion in the depth domain is:

[0057]

[0058] Wherein D newG is the sequence on the right side of the simultaneous inversion equation of time-lapse seismic new m is the matrix on the left side of the simultaneous inversion equation of time-lapse seismic new m is the matrix on the left side of the simultaneous inversion equation of time-lapse seismic, and γ represents a weight factor of a low-frequency constraint term.

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

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

[0061] A calibration module is configured to calibrate, based on the depth domain seismic wavelet matrix, base well logging data and monitoring well logging data to obtain corrected base well logging data and monitoring well logging data;

[0062] A construction module is configured to determine difference well logging data by subtracting the corrected base well logging data from the monitoring well logging data, construct an initial model by full-band interpolation according to seismic horizons and well logging horizon information and the difference well logging data and the corrected base well logging data, and perform low-pass filtering on the constructed initial model to obtain a low-frequency model, wherein the low-frequency model includes base data elastic parameters and difference elastic parameters;

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

[0064] A solution module is configured to substitute the low-frequency model into the time-lapse seismic simultaneous inversion objective function, and iteratively solve the objective function, and take an error of the first term being lower than a preset threshold as an iteration end condition, and output depth domain base data elastic parameters and difference elastic parameters satisfying the iteration end condition as time-lapse seismic inversion results.

[0065] In a third aspect, embodiments of the present disclosure provide an electronic device, comprising a processor, a communication interface, a memory and a communication bus, wherein the processor, the communication interface and the memory complete communication with each other through the communication bus;

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

[0067] The processor is configured to execute the program stored on the memory to implement the depth domain time-lapse seismic joint inversion method described above.

[0068] In a fourth aspect, the embodiments of the present disclosure provide a computer readable storage medium, having stored thereon a computer program, wherein the computer program, when executed by a processor, implements the deep domain time-lapse seismic joint inversion method described above.

[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] The deep domain time-lapse seismic joint inversion method provided by the embodiments of the present disclosure acquires logging data, extracts a deep domain seismic wavelet from the logging data based on a point spread function, and establishes a deep domain seismic wavelet matrix corresponding thereto; calibrates the base logging data and the monitoring logging data based on the deep domain seismic wavelet matrix to obtain corrected base logging data and monitoring logging data; determines difference logging data by subtracting the corrected base logging data from the monitoring logging data, constructs an initial model by full-band interpolation based on the seismic horizon and logging horizon information and the difference logging data and the corrected base logging data, and performs low-pass filtering on the constructed initial model to obtain a low-frequency model, wherein the low-frequency model includes base data elastic parameters and difference elastic parameters; determines a forward operator of time-lapse seismic simultaneous inversion, and establishes a time-lapse seismic simultaneous inversion objective function based on the forward operator, wherein a first term of the objective function is an error term between a forward simulation record and actual seismic data, and a second term is a low-frequency constraint term; substitutes the low-frequency model into the time-lapse seismic simultaneous inversion objective function, and iteratively solves the objective function, takes the error of the first term being lower than a preset threshold as an iteration end condition, and outputs the deep domain base data elastic parameters and the difference elastic parameters satisfying the iteration end condition as a time-lapse seismic inversion result. Through deep domain wavelet calculation, deep domain non-stationary convolution forward simulation, low-frequency model establishment, time-lapse seismic simultaneous inversion objective function establishment, and simultaneous solving of the base data elastic parameters and the difference elastic parameters, the workload of time domain time-lapse seismic inversion processing is reduced, the limitations of the existing time-lapse seismic method are overcome, the development of the oil reservoir can be better indicated, and new well deployment arrangements can be better arranged, which has great potential in the field of actual oil reservoir characterization applications. BRIEF DESCRIPTION OF DRAWINGS

[0071] The accompanying drawings, which are incorporated into and form a part of the specification, illustrate preferred embodiments consistent with the present disclosure and, together with the description, serve to explain the principles of the disclosure.

[0072] In order to more clearly illustrate the technical solutions in the embodiments of the present disclosure or the prior art, the accompanying drawings needed to be used in the embodiments or related description will be briefly introduced. Obviously, for those skilled in the art, other drawings can also be obtained without creative labor based on these drawings.

[0073] Figure 1 A schematic diagram of a deep domain time-lapse seismic joint inversion method is shown according to an embodiment of the present disclosure;

[0074] Figure 2 A schematic diagram of a well logging curve is shown according to an embodiment of the present disclosure;

[0075] Figure 3 A schematic diagram of a difference elastic parameter curve calculated by differencing two periods of well logging data is shown according to an embodiment of the present disclosure;

[0076] Figure 4 A schematic diagram of a deep domain wavelet matrix corresponding to base well logging data is shown according to an embodiment of the present disclosure;

[0077] Figure 5 A schematic diagram of a deep domain wavelet matrix corresponding to monitoring well logging data is shown according to an embodiment of the present disclosure;

[0078] Figure 6 A schematic diagram of a base data inversion result of well logging data in a deep domain time-lapse seismic simultaneous inversion noise-free time is shown according to an embodiment of the present disclosure;

[0079] Figure 7 A schematic diagram of an inversion result of a difference elastic parameter in a noise-free time is shown according to an embodiment of the present disclosure;

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

[0081] Figure 9 A schematic diagram of a difference elastic parameter inversion result of well logging data after adding random noise with a signal-to-noise ratio of 2 to a synthetic seismic record is shown according to an embodiment of the present disclosure;

[0082] Figure 10 A structural block diagram of a deep domain time-lapse seismic joint inversion device is shown according to an embodiment of the present disclosure;

[0083] Figure 11 A structural block diagram of an electronic device is shown according to an embodiment of the present disclosure. DETAILED DESCRIPTION

[0084] To make the purposes, technical solutions, and advantages of the embodiments of the present disclosure clearer, the technical solutions in the embodiments of the present disclosure will be described clearly and completely below with reference to the drawings in the embodiments of the present disclosure. Obviously, the described embodiments are only some but not all of the embodiments of the present disclosure. Based on the embodiments in the present disclosure, all other embodiments obtained by those of ordinary skill in the art without creative effort are within the scope of the present disclosure.

[0085] Referring to Figure 1 The embodiments of the present disclosure provide a deep domain time-lapse seismic joint inversion method, which comprises the following steps:

[0086] S1, obtaining logging data, extracting a deep domain seismic wavelet based on a point spread function according to the logging data, and establishing a deep domain seismic wavelet matrix corresponding to the deep domain seismic wavelet.

[0087] In the present embodiment, the logging data comprises a P-wave curve, a S-wave curve, and a density curve.

[0088] S2, calibrating the base logging data and the monitoring logging data based on the deep domain seismic wavelet matrix to obtain corrected base logging data and monitoring logging data.

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

[0090] S3, calculating the difference between the corrected base logging data and the monitoring logging data to determine difference logging data, constructing an initial model by full-band interpolation according to the seismic horizon and logging horizon information and the difference logging data and the corrected base logging data, and performing low-pass filtering processing on the constructed initial model to obtain a low-frequency model, wherein the low-frequency model comprises base data elastic parameters and difference elastic parameters.

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

[0092] S5, substituting the low-frequency model into the time-lapse seismic simultaneous inversion objective function, and iteratively solving the objective function, taking the error of the first term being lower than a preset threshold as an iteration end condition, and outputting the deep domain base data elastic parameters and the difference elastic parameters satisfying the iteration end condition as a time-lapse seismic inversion result.

[0093] In the present embodiment, the iterative solving of the objective function is realized based on a Gauss-Newton algorithm.

[0094] In the present embodiment, the method further comprises:

[0095] According to the time-lapse seismic inversion result, well seismic trace inversion testing is performed, logging results are compared with the basic data elastic parameter inversion result and the difference elastic parameter inversion result respectively, and at least one of the iteration number and the weight factor is adjusted according to the comparison result;

[0096] According to the adjusted iteration number and weight factor, parallel time-lapse seismic inversion is performed on the work area seismic data, and the depth domain basic data elastic parameter and the difference elastic parameter of the work area are obtained.

[0097] In this embodiment, in step S1, the depth domain seismic wavelet is extracted from the logging data based on the point spread function, and a depth domain seismic wavelet matrix corresponding thereto is established, including:

[0098] The time domain seismic wavelet is extracted from the logging data by angle, and the corresponding depth domain seismic wavelet is determined based on the point spread function according to the logging data and the time domain seismic wavelet, and a depth domain wavelet matrix is constructed.

[0099] In this embodiment, the corresponding depth domain seismic wavelet is determined based on the point spread function according to the logging data and the time domain seismic wavelet, and a depth domain wavelet matrix is constructed, including:

[0100] Under a one-dimensional velocity model, the point spread function and the time domain seismic wavelet satisfy the following relationship:

[0101]

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

[0103] Since the underground velocity is not uniform, assuming that the wavelet is a zero-phase wavelet, the depth domain wavelet at the velocity interface position is calculated by respectively calculating the depth domain wavelets corresponding to the upper and lower layer velocities, and cutting at the maximum value of the wavelet, and then splicing and interpolating to calculate the depth domain wavelet at the depth;

[0104] Based on the logging data, the depth domain wavelets at each depth are constructed into a Toeplitz matrix to obtain a depth domain wavelet matrix.

[0105] In this embodiment, in step S2, the depth domain seismic wavelet matrix is used to calibrate the basic logging data and the monitoring logging data to obtain corrected basic logging data and monitoring logging data, including:

[0106] The reflection coefficient sequence is determined according to the logging data, and a depth domain pre-stack angle gather is forward calculated according to a depth domain seismic wavelet matrix and the reflection coefficient sequence, and the depth domain pre-stack angle gather is calibrated with the monitoring logging data and the basic logging data to obtain corrected basic logging data and monitoring logging data.

[0107] In the embodiment, the reflection coefficient sequence is determined according to the logging data, including:

[0108] The logging data is substituted into the Zoeppritz equation to calculate the reflection coefficient sequence.

[0109] In the embodiment, the depth domain pre-stack angle gather is forward calculated according to the depth domain seismic wavelet matrix and the reflection coefficient sequence based on a non-stationary convolution algorithm.

[0110] In the embodiment, in step S4, the forward operator of the time-lapse seismic simultaneous inversion is determined, and a target function of the time-lapse seismic simultaneous inversion is established based on the forward operator, including:

[0111] The non-linear Zoeppritz forward operator is linearized at a low frequency model, and the inversion target functions of the basic data and the monitoring data are respectively established based on linear equations, and the derivatives of the two inversion target functions are respectively solved and combined to obtain the target function of the time-lapse seismic simultaneous inversion.

[0112] In the embodiment, the non-linear Zoeppritz forward operator is linearized at a low frequency model, and the inversion target functions of the basic data and the monitoring data are respectively established based on linear equations, and the derivatives of the two inversion target functions are respectively solved and combined to obtain the target function of the time-lapse seismic simultaneous inversion, including:

[0113] The expressions of the basic data and the monitoring data are respectively as follows:

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

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

[0116] Wherein, G represents a non-linear operator of Zoeppritz forward calculation; m represents a seismic elastic parameter; R represents a reflection coefficient sequence of Zoeppritz forward calculation; W represents a depth domain wavelet matrix; d1 and d2 respectively represent the basic data and the monitoring data,

[0117] The linearization of the non-linear operator adopts Taylor expansion, and the expansion expression after expansion at a low frequency model and retaining a linear part is as follows:

[0118]

[0119]

[0120] where m0 and m'0 represent the low-frequency model of the elastic parameters of the base data and the monitoring data, respectively,

[0121] Based on the expansion expression, the inversion equations of the base data and the monitoring data are respectively established:

[0122]

[0123]

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

[0125]

[0126]

[0127]

[0128]

[0129] 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 term, respectively; where L represents a difference operator, and the expression is:

[0130]

[0131] where nt represents the number of longitudinal sampling points of the data,

[0132] The inversion equations of the base data and the monitoring data are differentiated, respectively, and the derivatives are equal to 0, to obtain the following expressions:

[0133]

[0134]

[0135] The above expressions are rearranged and moved, and it is assumed that m2 = m1 + Δm, to obtain the simultaneous inversion equation of time-lapse seismic:

[0136]

[0137] where the expressions of G b , G m , D b , D m 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-shifted seismic data in the depth domain is:

[0141]

[0142] Among them, D new For the right-hand side sequence of the time-shifted earthquake simultaneous inversion equation, G new Let m be the matrix on the left-hand side of the time-shifted earthquake simultaneous inversion equation. new γ represents the left-hand side matrix of the time-shifted earthquake simultaneous inversion equation, and γ represents the weighting 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-shifted seismic data using the following expression, and the objective function is iteratively solved m. new :

[0144]

[0145] Where k is the number of iterations, and α is the iteration update step size, which is usually taken as 1. Let G be the derivative of the objective function with respect to the model parameters, which is numerically equal to G. new The iteration stops when the error value is less than the given error threshold, and the current iteration result is output.

[0146] Using the depth-domain time-shifted seismic joint inversion method disclosed herein, the inversion is performed. Figure 2 The time-lapse logging curves shown are illustrated in the following flowchart:

[0147] Will Figure 2 Substituting the time-lapse logging curve shown into the point diffusion function solution formula, the result can be calculated as follows: Figure 4 and Figure 5 The depth wavelet matrices corresponding to the two earthquake data periods shown are: Figure 3 To calculate the differential elastic parameter curve from the difference between two well logging data, where, Figure 4 and Figure 5From left to right respectively represent the longitudinal wave velocity, transverse wave velocity and density, the solid line represents the logging basic three parameter curve, the dotted line represents the logging monitoring three parameter curve, Figure 4 is the depth domain wavelet matrix corresponding to the basic logging data, Figure 5 is the depth domain wavelet matrix corresponding to the monitoring logging data, through Figure 4 and Figure 5 It can be known from the local enlargement in the box that, unlike the time domain seismic wavelet matrix, due to the change of seismic wave velocity in the longitudinal direction, the depth domain seismic wavelet produces distortion at the velocity boundary surface, and shows asymmetry. Figure 6 is the basic data inversion result of the depth domain time-lapse seismic simultaneous inversion without noise, Figure 7 is the inversion result of the difference elastic parameter at the noise time, wherein, Figure 6 From left to right respectively are the longitudinal wave velocity (V p ), transverse wave velocity (V s ), density (ρ) of the basic data, Figure 6 From left to right respectively are the difference longitudinal wave velocity (ΔV p ), difference transverse wave velocity (ΔV s ), difference density (Δρ) of the difference elastic parameter, the curve shown in 601 is the real logging data and the difference logging data, and the curve shown in 602 is the inverted logging data and the difference logging data. As Figure 8 and Figure 9 shown are the logging data inversion results after adding random noise with a signal-to-noise ratio of 2 into the synthetic seismic record, wherein, Figure 8 From left to right respectively are the longitudinal wave velocity (V p ), transverse wave velocity (V s ), density (ρ) of the basic data at the signal-to-noise ratio of 2, Figure 9 From left to right respectively are the difference longitudinal wave velocity (ΔV p ), difference transverse wave velocity (ΔV s ), difference density (Δρ) of the difference elastic parameter at the signal-to-noise ratio of 2, the curve shown in 701 is the real logging data and the difference logging data, and the curve shown in 702 is the inverted logging data and the difference 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 two-period data to calculate the elastic parameters and the difference elastic parameters of the basic data, and through the anti-noise test, it can be seen that the method has a certain robustness, can be applied to seismic data inversion, and since the present disclosure overcomes the limitation of the seismic wavelet and the difference parameter, it can meet the more extensive time-lapse seismic reservoir prediction demand while meeting the inversion accuracy.

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

[0150] Referring to Figure 10 The embodiment of the present disclosure provides a depth domain time-lapse seismic joint inversion device, comprising:

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

[0152] The calibration module 12 is configured to calibrate the basic logging data and the monitoring logging data based on the depth domain seismic wavelet matrix, and obtain corrected basic logging data and monitoring logging data;

[0153] The construction module 13 is configured to determine difference logging data by subtracting the corrected basic logging data from the monitoring logging data, construct an initial model by performing full-band interpolation on the difference logging data and the corrected basic logging data according to seismic horizons and logging horizon information, and obtain a low-frequency model by performing low-pass filtering processing on the constructed initial model, wherein the low-frequency model comprises basic data elastic parameters and difference elastic parameters;

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

[0155] The solving module 15 is configured to substitute the low-frequency model into the objective function of simultaneous inversion of time-lapse seismic, and solve the objective function iteratively, take an error of the first term being lower than a preset threshold as an iteration end condition, and output the depth domain basic data elastic parameters and the difference elastic parameters satisfying the iteration end condition as a time-lapse seismic inversion result.

[0156] The device for depth domain time-lapse seismic joint inversion of the present disclosure overcomes the limitation of wavelet and difference data in difference inversion while retaining the precision of time-lapse seismic difference inversion, and meets the needs of practical application more; and is developed in the depth domain, avoiding the time-depth conversion process in the time domain, simplifying the well-seismic calibration process, and reducing the matching workload of time-lapse seismic data; compared with linear approximation, it overcomes the problem of reduced precision for large-angle incidence, and can be applied to time-lapse seismic inversion in any angle range.

[0157] The implementation process of the functions and roles of each unit in the above device is specifically described in the implementation process of the corresponding steps in the above method, which will not be repeated here.

[0158] For the device embodiment, since it basically corresponds to the method embodiment, the relevant part can be referred to the part of the method embodiment. The device embodiments described above are only illustrative, and the units described as separate components can or can not be physically separated, and the components displayed as units can or can not be physical units, i.e. they can be located in one place or distributed on multiple network units. Part or all of the modules can be selected to achieve the purpose of the present application scheme according to actual needs. Those skilled in the art can understand and implement without creative labor.

[0159] In the above embodiments, any multiple of the extraction module 11, the calibration module 12, the construction module 13, the establishment module 14 and the solving module 15 can be combined in one module, or any one of them can be split into multiple modules. Alternatively, at least part of the function of one or more of these modules can be combined with at least part of the function 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 solving 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 any other reasonable way of integrating or packaging a circuit, etc. hardware or firmware, or any one of software, hardware and firmware or any appropriate combination of several of them. Alternatively, at least one of the extraction module 11, the calibration module 12, the construction module 13, the establishment module 14 and the solving module 15 can be at least partially implemented as a computer program module which can perform corresponding functions when running.

[0160] Reference Figure 11As shown, the electronic device provided by the embodiment of the present disclosure includes a processor 1110, a communication interface 1120, a memory 1130 and a communication bus 1140, wherein 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 for storing a computer program.

[0162] The processor 1110 is used for executing the program stored in the memory 1130 to realize the depth domain time-lapse seismic joint inversion method as shown below:

[0163] Logging data is acquired, a depth domain seismic wavelet is extracted from the logging data based on a point spread function, and a depth domain seismic wavelet matrix corresponding to the depth domain seismic wavelet is established;

[0164] 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 corrected basic logging data and monitoring logging data;

[0165] Difference logging data is determined by subtracting the corrected basic logging data from the monitoring logging data, an initial model is constructed by full-band interpolation based on seismic horizons and logging horizon information and the difference logging data and the corrected basic logging data, and low-frequency model is obtained by low-pass filtering processing on the constructed initial model, wherein the low-frequency model includes basic data elastic parameters and difference elastic parameters;

[0166] A forward operator of time-lapse seismic simultaneous inversion is determined, and a time-lapse seismic simultaneous inversion objective function is established based on the forward operator, wherein a first term of the objective function is an error term between forward simulation records and actual seismic data, and a second term is a low-frequency constraint term;

[0167] The low-frequency model is substituted into the time-lapse seismic simultaneous inversion objective function, and iterative solution of the objective function is performed, the error of the first term is lower than a preset threshold as an iteration end condition, and the depth domain basic data elastic parameters and the difference elastic parameters meeting the iteration end condition are output as time-lapse seismic inversion results.

[0168] The communication bus 1140 described above can be a peripheral component interconnect (PCI) bus or an extended industry standard architecture (EISA) bus, etc. The communication bus 1140 can be divided into an address bus, a data bus, a control bus, etc. For the convenience of representation, only one thick line is shown 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 configured to enable communication between the electronic device described above and other devices.

[0170] The memory 1130 can include a random access memory (RAM) and can further include a non-volatile memory, e.g., at least one disk storage. Optionally, the memory 1130 can also be at least one memory storage located remotely from the aforementioned processor 1110.

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

[0172] Embodiments of the present disclosure further provide a computer-readable storage medium. The computer-readable storage medium described above stores a computer program, and the computer program is executed by a processor to implement the depth domain time-lapse seismic joint inversion method described above.

[0173] The computer-readable storage medium can be included in the device / apparatus described in the above embodiments; or can exist separately and not be assembled into the device / apparatus. The computer-readable storage medium described above carries one or more programs, and when the 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 embodiments of the present disclosure, the computer-readable storage medium can be a non-volatile computer-readable storage medium, which can include, but is not limited to, a portable computer diskette, a hard disk, a random access memory (RAM), a read-only memory (ROM), an erasable programmable read-only memory (EPROM or flash memory), a portable compact disc read-only memory (CD-ROM), an optical storage device, a magnetic storage device, or any appropriate combination thereof. In the present disclosure, the computer-readable storage medium can be any tangible medium that contains or stores a program, which can be used by or in connection with an instruction execution system, apparatus, or device.

[0175] It should be noted that, in this document, relational terms such as "first" and "second" are used merely to distinguish one entity or operation from another, and do not necessarily require or imply any such actual relationship or order between these entities or operations. Furthermore, the terms "comprising," "including," or any other variations thereof are intended to cover non-exclusive inclusion, such that a process, method, article, or apparatus that comprises a list of elements includes not only those elements but also other elements not expressly listed, or elements inherent to such a process, method, article, or apparatus. Without further limitations, an element defined by the phrase "comprising one..." does not exclude the presence of other identical elements in the process, method, article, or apparatus that includes said element.

[0176] The above description is merely a specific embodiment of this disclosure, enabling those skilled in the art to understand or implement it. Various modifications to these embodiments will be readily apparent to those skilled in the art, and the general principles defined herein may be implemented in other embodiments without departing from the spirit or scope of this disclosure. Therefore, this disclosure is not to be limited to the embodiments shown herein, but is to be accorded the widest scope consistent with the principles and novel features claimed herein.

Claims

1. A method for joint inversion of time-shifted seismic data in the depth domain, characterized in that, The method includes: Acquire well logging data, extract depth domain seismic wavelets based on point spread function, and establish the corresponding depth domain seismic wavelet matrix; Based on the depth domain seismic wavelet matrix, well-seismic calibration is performed on the basic logging data and monitoring logging data to obtain the corrected basic logging data and monitoring logging data. The difference between the corrected basic logging data and the monitoring logging data is calculated to determine the differential logging data. Based on the seismic and logging horizon information, as well as the differential logging data and the corrected basic logging data, a full-band interpolation is performed to construct an initial model. The constructed initial model is then subjected to low-pass filtering to obtain a low-frequency model, wherein the low-frequency model includes the elastic parameters of the basic data and the differential elastic parameters. A forward modeling operator for simultaneous time-shifted earthquake inversion is determined, and an objective function for simultaneous time-shifted earthquake inversion is established based on this forward modeling operator. The first term of the objective function is the error term between the forward modeling simulation record and the actual earthquake data, and the second term is the low-frequency constraint term. The low-frequency model is substituted into the objective function of the time-shifted seismic inversion, and the objective function is solved iteratively. The error of the first term is lower than a preset threshold as the iteration termination condition. The elastic parameters and differential elastic parameters of the depth domain basic data that meet the iteration termination condition are output as the time-shifted seismic inversion result.

2. The method according to claim 1, characterized in that, The method further includes: Based on the time-shifted seismic inversion results, a well-side seismic trace inversion test is performed. The well logging results are compared with the elastic parameter inversion results of the basic data and the differential elastic parameter inversion results, respectively. Based on the comparison results, at least one of the iteration number and weighting factors is adjusted. Parallel time-shifted seismic inversion of the seismic data in the work area is performed based on the adjusted iteration number and weighting factor to obtain the elastic parameters and differential elastic parameters of the basic depth domain data of the work area.

3. The method according to claim 1, characterized in that, The step of extracting depth-domain seismic wavelets from well logging data based on the point spread function and establishing the corresponding depth-domain seismic wavelet matrix includes: Time-domain seismic wavelets are extracted from well logging data by angle, and based on the point spread function, the corresponding depth-domain seismic wavelets are determined according to the well logging data and the time-domain seismic wavelets, and the depth-domain wavelet matrix is ​​constructed.

4. The method according to claim 3, characterized in that, The method based on the point spread function determines the corresponding depth-domain seismic wavelet according to well logging data and time-domain seismic wavelet, and constructs a depth-domain wavelet matrix, including: Under the one-dimensional velocity model, the point propagation function satisfies the following relationship with the time-domain seismic wavelet: in, The time-domain wavelet period is used to represent the integral of the velocity of the spatial point spread function over one wavelength, where λ is the wavelength. h For depth, v The velocity corresponding to this depth; Assuming the wavelet is a zero-phase wavelet, the depth domain wavelet at the velocity interface is calculated by separately calculating the depth domain wavelets corresponding to the velocities of the upper and lower layers, truncating at the maximum value of the wavelet, and then splicing and interpolating to calculate the depth domain wavelet at that depth. Based on well logging data, Toeplitz matrices are constructed from depth domain wavelets at various depths to obtain depth domain wavelet matrices.

5. The method according to claim 1, characterized in that, The method involves performing well-seismic calibration on basic logging data and monitoring logging data based on the depth-domain seismic wavelet matrix, resulting in corrected basic logging data and monitoring logging data, including: Based on the well logging data, the reflection coefficient sequence is determined, and based on the depth domain seismic wavelet matrix and the reflection coefficient sequence, the depth domain pre-stack angle gather is forward modeled. Well-seismic calibration is then performed on the depth domain pre-stack angle gather, monitoring well logging data, and basic well logging data to obtain the corrected basic well logging data and monitoring well logging data.

6. The method according to claim 1, characterized in that, The determination of the forward modeling operator for simultaneous time-shifted earthquake inversion, and the establishment of the objective function for simultaneous time-shifted earthquake inversion based on the forward modeling operator, includes: The nonlinear Zoeppritz forward modeling operator is linearized in the low-frequency model. Inversion objective functions for basic data and monitoring data are established based on the linear equations. The derivatives and simultaneous equations of the two inversion objective functions are obtained to obtain the time-shifted earthquake simultaneous inversion objective function.

7. The method according to claim 6, characterized in that, The nonlinear Zoeppritz forward modeling operator is linearized in the low-frequency model. Inversion objective functions for basic data and monitoring data are established based on linear equations. The derivatives and simultaneous equations of the two inversion objective functions are then obtained to yield the simultaneous inversion objective function for time-shifted earthquakes, including: The expressions for basic data and monitoring data are as follows: in, This represents a nonlinear operator for the Zoeppritz forward model; Represents seismic elastic parameters; This represents the Zoeppritz forward modeling sequence of reflection coefficients; The depth domain wavelet matrix is ​​represented; d1 and d2 represent the basic data and monitoring data, respectively. Linearization of the nonlinear operator is achieved using Taylor expansion. The expansion expression after expanding at the low-frequency model and preserving the linear part is as follows: in, and Low-frequency models representing the elasticity parameters of basic data and monitoring data, respectively. Based on the expanded expression, inversion equations are established for both basic data and monitoring data: The first term is the error term. , To expand the corresponding gradient and intercept terms in expression d1, , To expand the gradient and intercept terms in expression d2, the specific form is as follows: The second term is a low-frequency constraint term, and the third term is a smoothing term. and Let represent the weighting factors of the constraint terms; where L represents the difference operator, whose expression is: Where nt represents the number of vertical sampling points of the data. Differentiating the inversion equations for the basic data and monitoring data respectively, and setting the derivatives to 0, we obtain the following form: Rearrange and rearrange the terms in the above expression, and assume... The simultaneous inversion equation for time-shifted earthquakes is obtained: in, , , , They each have the following forms: ， ， The objective function for simultaneous inversion of time-shifted seismic data in the depth domain is: in, This is the right-hand side sequence of the simultaneous inversion equation for time-shifted earthquakes. The left-hand side matrix of the time-shifted earthquake simultaneous inversion equation is... The left-hand side matrix of the time-shifted earthquake simultaneous inversion equation is... This represents the weighting factor for the low-frequency constraint term.

8. A depth-domain time-shifted seismic joint inversion device, characterized in that, include: The extraction module is used to acquire well logging data, extract depth domain seismic wavelets based on the point spread function, and establish the corresponding depth domain seismic wavelet matrix. The calibration module is used to perform well-seismic calibration on basic logging data and monitoring logging data based on the depth domain seismic wavelet matrix, so as to obtain the corrected basic logging data and monitoring logging data. The construction module 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 and logging horizon information, as well as the differential logging data and the corrected basic logging data, a full-band interpolation is performed to construct an initial model. The constructed initial model is then subjected to low-pass filtering to obtain a low-frequency model. The low-frequency model includes basic data elastic parameters and differential elastic parameters. A module is established to determine the forward modeling operator for simultaneous time-shifted earthquake inversion, and an objective function for simultaneous time-shifted earthquake inversion is established based on the forward modeling operator. The first term of the objective function is the error term between the forward modeling simulation record and the actual earthquake data, and the second term is the low-frequency constraint term. The solution module is used to substitute the low-frequency model into the objective function of time-shifted seismic inversion and perform iterative solution of the objective function. The error of the first term is lower than a preset threshold as the iteration termination condition, and the depth domain basic data elastic parameters and differential elastic parameters that meet the iteration termination condition are output as the time-shifted seismic inversion results.

9. An electronic device, characterized in that, It includes a processor, a communication interface, a memory, and a communication bus, wherein the processor, the communication interface, and the memory communicate with each other through the communication bus; Memory, used to store computer programs; The processor, when executing the program stored in the memory, implements the depth-domain time-shifted seismic joint inversion method as described in any one of claims 1-7.

10. A computer-readable storage medium having a computer program stored thereon, characterized in that, When the computer program is executed by the processor, it implements the depth-domain time-shifted seismic joint inversion method as described in any one of claims 1-7.

Images