Direct well-tie method for depth-domain logging and seismic data

[0001] This application claims the benefit of priority from Chinese Patent Application No. 202210180144.4, filed on Feb. 25, 2022. The content of the aforementioned application, including any intervening amendments thereto, is incorporated herein by reference in its entirety.

[0002] This application relates to oil-gas seismic exploration, and more particularly to a direct well-tie method for depth-domain logging and seismic data.

[0003] During the geophysical logging process, various logging instruments manufactured on the basis of electricity, magnetism, sound, heat, and nuclear are deployed in the well to continuously collect a series of physical parameters varying with depth along the well, and variation curves of these physical parameters are plotted to help identify the subsurface lithology and mineral deposits such as oil, gas, water, coal, and metals. At the logging point, since the drill cores are obtained and the physical parameters are measured during logging, the underground geological conditions at this point are clearly understood. Currently, the seismic data obtained by artificial seismic exploration is mainly employed to infer the underground geological conditions of the site far from the logging point. The seismic data interpretation is to give the seismic data characteristics (such as reflection characteristics) a clear geological meaning and a conceptual model, so as to clarify the underground geological conditions. In the seismic data interpretation, the calibration of logging and seismic data not only plays a key role in correlating the subsurface geological conditions with observed seismic data, but also forms the foundations for the reliable stratigraphic interpretation and seismic inversion. The logging and seismic data calibration is performed through steps of: making a synthetic seismogram based on the velocity and density parameters obtained by logging, and aligning the obvious waveform features in the synthetic seismogram with the waveform features of the seismic data near the logging point (the alignment needs to be repeated many times), so as to correlate the features of the logging physical parameter curve with those of the seismic data to realize the extrapolation of geological conditions at the logging point, thereby inferring the subsurface geological conditions of a site far away from the logging point.

[0004] Since most of the seismic data employed for the previous interpretation are time-domain seismic data, it is required to convert the physical logging parameters varying with depth into physical logging parameters varying with time according to the logging velocity in the logging and seismic data calibration, and then synthesize the time-domain seismogram in the time domain based on the convolution model, so as to carry out calibration with the time-domain seismic data. Time-domain logging and seismic data calibration techniques have been well established. However, as the depth-domain seismic data with higher imaging accuracy of subsurface structures has been increasingly employed in the oil-gas seismic exploration, the seismic data interpretation is increasingly conducted relying on the depth-domain seismic data. The existing depth-domain logging and seismic data calibration technology corresponding to the depth-domain seismic data interpretation still remains to be improved. In the prior art, the depth-domain logging data and depth-domain seismic data are respectively converted into the time-domain logging data and time-domain seismic data to complete the calibration, and then the calibration results are converted back to the depth domain. Such calibration process is complex and requires multiple domain conversions, which will introduce conversion errors in the logging data and seismic data, resulting in low calibration accuracy, and even the matching of the logging data and seismic data at the wrong depth position.

[0005] An objective of this application is to provide a direct well-tie method for depth-domain logging and seismic data, where a direct depth-domain seismogram synthesis method and a seismic wavelet estimation method are integrated in the calibration process, namely, there is no domain conversion in the seismogram synthesis and the seismic wavelet estimation. Therefore, the method provided herein can not only make full use of the logging information, but also avoid the accumulation of errors generated in the conversion between different domains and the data information loss brought by resampling (existing in the indirect method). In the method provided herein, a relatively accurate logging depth-seismic depth corresponding relationship can be established between the logging data and the seismic data through the first logging-seismic calibration, thereby reducing the number of subsequent local adjustments of the depth correspondence between the logging data and the seismic data.

[0006] Technical solutions of this application are described as follows.

[0007] This application provides a direct well-tie method for depth-domain logging and seismic data, comprising:

• • [0008] (1) removing null values and outlier values from logging velocity data and logging density data to obtain valid logging velocity data and valid logging density data;
  • [0009] (2) based on the valid logging velocity data and the valid logging density data, calculating a logging reflection coefficient r at different depths, expressed as:
    r=[r(1), . . . ,r(l), . . . ,r(L)];
  • [0010] wherein L indicates the number of sampling points for the valid logging velocity data and the valid logging density data;
  • [0011] (3) estimating an initial zero-phase depth-domain seismic wavelet from a seismic section S, wherein the seismic section S is within a target depth range and passes through a well location, and the initial zero-phase depth-domain seismic wavelet is estimated through the steps of:
    • [0012] (3.1) calculating an average spectrum y of the seismic section S through the following equation:

[0013] $y=\frac{1}{N}{\sum }_{n=1}^N{\left(fft\left[S\left(n\right)\right]\right)}^2,n\in \left[1,2,\dots ,N\right];$

• • • [0014] wherein fft[⋅] represents vector's fast Fourier transform; the seismic section S contains N seismograms in total; and S(n) indicates an n-th seismogram in the seismic section S, and contains M sampling points;
    • [0015] (3.2) building an objective function

[0016] $\underset{\alpha ,\beta }{\min }{\exp \left[-{\left(\frac{K-\alpha }{\beta }\right)}^2\right]-y}_2,$

• • • [0017]  and solving the objective function by using least squares method to obtain a fitting parameter α and a fitting parameter β; wherein ∥∥₂ represents a vector L2 norm; exp[⋅] represents a vector exponent; and a seismic wavenumber vector K is expressed as K=[K(1), . . . , K(m), . . . , K(M)]; wherein

[0018] $K\left(m\right)=\frac{\left(m-1\right)}{M\Delta D},$

• • • [0019]  a unit of α is 1/km; ΔD indicates a depth sampling interval of the seismic section S, a unit of ΔD is m; and a unit of each of the fitting parameter α and the fitting parameter β is 1/km;
    • [0020] (3.3) constructing the initial zero-phase depth-domain seismic wavelet w₀ based on the fitting parameter β and a logging wavenumber vector k, expressed as follows:

[0021] $w_0=\frac{\left(45\pi -128\right)k^2}{9\pi {\beta }^2}\exp \left[-\left(\frac{\left(45\pi -128\right)k^2}{18{\mathrm{\pi \beta }}^2}\right)+1\right];$

• • • [0022] wherein k=[k(1), . . . , k(l), . . . , k(L)];

[0023] $k\left(l\right)=\frac{\left(l-1\right)}{L\Delta d},$

• • • [0024]  and Δd indicates a depth sampling interval of logging data, with a unit of m;
  • [0025] (4) making a depth-domain synthetic seismogram A based on the initial zero-phase depth-domain seismic wavelet w₀ and a logging reflection coefficient r, expressed as follows:
    A=real{P[rP)⊙w₀]};
  • [0026] wherein real{⋅} denotes an operation of taking a vector real part; ⊙ represents a vector Hadamard product; P is an L-order matrix; P is a conjugate matrix of the P; and an element P(l,l) in the P is expressed as follows:

[0027] $P\left(l,l\right)=\frac{P\left(l-1,l\right)}{L}\exp \left[-\frac{i2\pi \left(l-1\right)v_{\max }}{v\left(l\right)L}\right];$

• • [0028] wherein i=√{square root over (−1)} is an imaginary unit; and v_max is a maximum value in the valid logging velocity data v; a unit of v is m/s; at this time, a depth axis D^A corresponding to the depth-domain synthetic seismogram A is expressed as: D^A=[D^A(1), . . . , D^A(l), . . . , D^A(L)]; wherein D^A(1) is an initial depth of a depth-domain seismogram at the well location;
  • [0029] (5) according to the depth sampling interval Δd of the logging data, subjecting the depth-domain seismogram at the well location to interpolation by using a cubic spline interpolation method to obtain an interpolated depth-domain seismogram at the well location;
  • [0030] (6) performing cross-correlation operation between the depth-domain synthetic seismogram A and the interpolated depth-domain seismogram at the well location to obtain a cross-correlation result; updating the depth axis D^A corresponding to the depth-domain synthetic seismogram A into D^A′=D^A+[D₀−D^A(1)] according to a depth position D₀ corresponding to a maximum value in the cross-correlation result; at the same time, starting from D₀, extracting L sampling points from the interpolated depth-domain seismogram at the well location to obtain a depth-domain seismogram segment B; wherein a depth axis D^B corresponding to the depth-domain seismogram segment B is expressed as D^B=[D^B(1), . . . , D^B(l), . . . , D^B(L)];
  • [0031] (7) based on the depth-domain seismogram segment B and the logging reflection coefficient r, estimating a depth-domain seismic wavelet by using a direct depth-domain seismic wavelet estimation method; based on the depth-domain seismic wavelet and the logging reflection coefficient r, updating the depth-domain synthetic seismogram A by using a direct depth-domain seismogram synthesis method to obtain a depth-domain synthetic seismogram A′;
  • [0032] (8) selecting a wave peak or trough position from the depth-domain seismogram segment B, wherein the wave peak or trough position is similar to the depth-domain synthetic seismogram A′ in waveform characteristic; compressing or stretching the depth-domain seismogram segment B to align the wave peak or trough position with the depth-domain synthetic seismogram A′, and generating a corresponding depth deviation point set E, expressed as follows:
    E={E(1), . . . ,E(j), . . . ,E(J)},J<L;
  • [0033] wherein E(j)=D^B(j)−D^A′(j), representing a depth axis difference between the depth-domain seismogram segment B and the depth-domain synthetic seismogram A′ at a j-th point; and
  • [0034] subjecting elements in the depth deviation point set E to interpolation by using the cubic spline interpolation method to obtain a depth calibration curve C;
  • [0035] (9) subjecting the depth calibration curve C and the depth axis D^A′ to addition to obtain a depth axis D^A″=D^A′+C, and subjecting the depth calibration curve C and the depth axis D^B to addition to obtain a depth axis D^B′=D^B+C; according to the depth axis D^B′, intercepting a corresponding depth-domain seismogram segment B′ from the interpolated depth-domain seismogram at the well location; and repeating step (7); and
  • [0036] (10) determining whether the depth-domain seismogram segment B′ is aligned with the depth-domain synthetic seismogram A′ at the wave peak or trough position, and a correlation coefficient between the depth-domain seismogram segment B′ and the depth-domain synthetic seismogram A′ is higher than a correlation coefficient in a previous calibration process; if yes, completing a calibration operation, and generating a logging depth-seismic depth corresponding relationship and a depth-domain seismic wavelet; otherwise, repeating steps (8)-(9).

[0055] FIGS. 2a-2i schematically shows a process of calibrating logging data and depth-domain seismogram (shown in FIG. 1e) at the logging location, which is performed through the following steps.

• • [0056] (1) Null values and outlier values are removed from logging velocity data and logging density data to obtain valid logging velocity data v and valid logging density data ρ, where the logging velocity data is collected by means of a velocity sensor, and the logging density data is collected through a density sensor.
  • [0057] (2) Based on the valid logging velocity data v and the valid logging density data ρ, a logging reflection coefficient r at different depths, which is expressed as follows:
    r=[r(1), . . . ,r(l), . . . ,r(L)];
  • [0058] where

[0059] $r\left(l\right)=\frac{v\left(l+1\right)\rho \left(l+1\right)-v\left(l\right)\rho \left(l\right)}{v\left(l+1\right)\rho \left(l+1\right)+v\left(l\right)\rho \left(l\right)},$

• • [0060]  and L=1389 and L indicates the number of sampling points for the valid logging velocity data v and the valid logging density data ρ.
  • [0061] (3) An initial zero-phase depth-domain seismic wavelet is estimated from a seismic section S, where the seismic section S is within a target depth ranging from 2800 m to 3400 m and passes through an inline direction of a well location. The initial zero-phase depth-domain seismic wavelet is estimated through the following steps.
    • [0062] (3.1) An average spectrum y of the seismic section S is calculated through the following equation:

[0063] $\begin{array}{cc}y=\frac{1}{N}{\sum }_{n=1}^N{\left(\mathrm{fft}\left[S\left(n\right)\right]\right)}^2,& n\in \left[1,2,\dots ,N\right];\end{array}$

• • • [0064]  where fft[⋅] represents a vector's fast Fourier transform, the seismic section S contains a total of N=965 seismograms, and S(n) indicates an n-th seismogram in the seismic section S, and contains M=121 sampling points;
    • [0065] (3.2) An objective function

[0066] $\underset{\alpha ,\beta }{\min }{\exp \left[-{\left(\frac{K-\alpha }{\beta }\right)}^2\right]-y}_2$

• • • [0067]  is built and solved by using least squares method to obtain a fitting parameter α and a fitting parameter β. In the objective function, ∥∥₂ represents a vector L2 norm; exp[⋅] represents a vector exponent; and a seismic wavenumber vector K is expressed as K=[K(1), . . . , K(m), . . . , K(M)], where

[0068] $K\left(m\right)=\frac{\left(m-1\right)}{M\Delta D};$

• • • [0069]  ΔD indicates a depth sampling interval of the seismic section S, and ΔD=5 m; and a unit of each of the fitting parameter α and the fitting parameter β is 1/km.
    • [0070] (3.3) The initial zero-phase depth-domain seismic wavelet w₀ is constructed based on the fitting parameter β and a logging wavenumber vector k, expressed as follows:

[0071] $w_0=\frac{\left(45\pi -128\right)k^2}{9{\mathrm{\pi \beta }}^2}\exp \left[-\left(\frac{\left(45\pi -128\right)k^2}{18{\mathrm{\pi \beta }}^2}\right)+1\right];$

• • • [0072] where k=[k(1), . . . , k(l), . . . , k(L)],

[0073] $k\left(l\right)=\frac{\left(l-1\right)}{L\Delta d}$

• • • [0074]  and Δd indicates a depth sampling interval of logging data, and Δd=0.15 m.
  • [0075] (4) A depth-domain synthetic seismogram A is made based on the initial zero-phase depth-domain seismic wavelet w₀ and a logging reflection coefficient r, expressed as follows:
    A=real{P[(rP)⊙w₀]};
  • [0076] where real{⋅} denotes an operation of taking a vector real part; ⊙ represents a vector Hadamard product; P is an L-order matrix; P is a conjugate matrix of the P; and an element P(l,l) in the P is expressed as follows:

[0077] $P\left(l,l\right)=\frac{P\left(l-1,l\right)}{L}\exp \left[-\frac{i2\pi \left(l-1\right)v_{\max }}{v\left(l\right)L}\right];$

• • [0078] where i=√{square root over (−1)} is an imaginary unit; and v_max is a maximum value in the valid logging velocity data v, and v_max=3954.2 m/s; at this time, a depth axis D^A corresponding to the depth-domain synthetic seismogram A is expressed as: D^A=[D^A(1), . . . , D^A(l), . . . , D^A (L)], where D^A (1) is an initial depth of a depth-domain seismogram at the well location, and D^A (1)=2400 m.
  • [0079] (5) According to the depth sampling interval of the logging data, the depth-domain seismogram at the well location is subjected to interpolation by using a cubic spline interpolation method to obtain an interpolated depth-domain seismogram at the logging location.
  • [0080] (6) Cross-correlation operation between the depth-domain synthetic seismogram A and the interpolated depth-domain seismogram at the logging location is performed to obtain a cross-correlation result. The depth axis D^A corresponding to the depth-domain synthetic seismogram A is updated into D^A′=D^A+[D₀−D^A(1)], according to a depth position D₀ (i.e., 2941.2 m) corresponding to a maximum value in the cross-correlation result. At the same time, starting from D₀, L sampling points are extracted from the interpolated depth-domain seismogram at the logging location to obtain a depth-domain seismogram segment B, where a depth axis D^B corresponding to the depth-domain seismogram segment B is expressed as D^B=[D^B(1), . . . , D^B(l), . . . , D^B(L)].
  • [0081] (7) Based on the depth-domain seismogram segment B and the logging reflection coefficient r, a depth-domain seismic wavelet is estimated by using a model-based depth-domain seismic wavelet estimation method (referring to a method disclosed by Chinese Patent Application No. 202010504205.9). Based on the depth-domain seismic wavelet and the logging reflection coefficient r, the depth-domain synthetic seismogram A is updated by using a direct depth-domain seismogram synthesis method (referring to a method disclosed by Chinese Patent Application No. 202010504205.9) to obtain a depth-domain synthetic seismogram A′.
  • [0082] (8) A wave peak or trough position is selected from the depth-domain seismogram segment B, where the wave peak or trough position is similar to a corresponding wave peak or trough position in the depth-domain synthetic seismogram A′ in waveform characteristic, and has a large absolute amplitude (as shown in FIG. 3a). The depth-domain seismogram segment B is compressed and stretched to align the wave peak or trough position with the corresponding wave peak or trough position in the depth-domain synthetic seismogram A′, and a corresponding depth deviation point set E (shown as black mark points in FIG. 3b) is generated, which is expressed as follows:
    E={E(1), . . . ,E(j), . . . ,E(J)},J<L;
  • [0083] where E(j)=D^B(j)−D^A′(j), representing a depth axis difference between the depth-domain seismogram segment B and the depth-domain synthetic seismogram A′ at a j-th point.

[0084] Specifically, the “wave peak or trough position” is determined based on the combination of visual observation and correlation coefficient. Whether the depth-domain seismogram segment B is aligned with the depth-domain synthetic seismogram A′ is determined by human-computer interaction, and then the correlation coefficient is introduced to evaluate the similarity to select the desired wave peak or trough in the depth-domain seismogram segment B.

[0085] Elements in the depth deviation point set E are subjected to interpolation by using the cubic spline interpolation method to obtain a depth calibration curve C (shown as black solid line in FIG. 3b).

• • [0086] (9) The depth calibration curve C and the depth axis D^A′ are subjected to addition to obtain a depth axis D^A″=D^A′+C, and the depth calibration curve C and the depth axis D^B are subjected to addition to obtain a depth axis D^B′=D^B+C. According to the depth axis D^B′ (as shown in the abscissa in FIG. 3c), a corresponding depth-domain seismogram segment B′ is intercepted from the interpolated depth-domain seismogram at the logging location. Step (7) is repeated.
  • [0087] (10) Whether the depth-domain seismogram segment B′ is aligned with the depth-domain synthetic seismogram A′ at the wave peak or trough position, and a correlation coefficient between the depth-domain seismogram segment B′ and the depth-domain synthetic seismogram A′ is higher than a correlation coefficient in a previous calibration process (as shown in FIG. 2e) is determined. If yes, a calibration operation is completed, and a logging depth-seismic depth corresponding relationship and a depth-domain seismic wavelet are generated; otherwise, steps (8)-(9) are repeated.

[0088] It can be concluded that after the logging data and the depth-domain seismogram are subjected to calibration operation by using the method provided herein for two times (shown in FIGS. 2c and 2d), a reasonable logging depth-seismic depth corresponding relationship (shown in FIG. 4) and a depth-domain seismic wavelet are obtained (shown in FIG. 2i).

[0089] The Gamma Ray logging curve can be employed for lithology classification and location. For the sand and mudstone section, the pure sandstone has the lowest Gamma value, and the Gamma value rises with the increase of shale content. The Gamma Ray logging curve is projected onto the seismic section image based on the obtained logging depth-seismic depth corresponding relationship through the execution of a computer program on a processor. Referring to the through-well depth-domain seismic section shown in FIG. 5, low-value areas in the Gamma Ray logging curve are consistent with the two sets of sand bodies (predicted based on the available geological knowledge and available logging information), indicating that the method provided herein is effective and reliable, especially in the location of sand bodies (e.g., sandstone, which is applicable to the oil-gas exploration). By projecting the Gamma Ray logging curve onto the seismic section image based on the logging depth-seismic depth corresponding relationship obtained by the method provided herein, the logging data is introduced to predict the underground medium in a larger three-dimensional (3D) space, facilitating the generation of geological models (applied to the identification of geological characteristics and reservoir properties of oil and gas reservoirs, and calculation of oil and gas reserves) and location of sand bodies.

[0090] The depth-domain seismic wavelet obtained by the method provided herein is applicable to the subsequent seismic inversion to output an elastic impedance image of the underground medium.

[0091] Compared with the prior art, this application has the following beneficial effects.

[0092] In this application, a direct depth-domain seismogram synthesis method and a seismic wavelet estimation method are integrated in the calibration process, which can not only make full use of the logging information, but also avoid the accumulation of errors generated in the conversion between different domains (such as depth domain and time domain) and the data information loss brought by resampling (existing in the indirect method). In the method provided herein, a relatively accurate depth corresponding relationship can be established between the logging data and the seismic data through the first logging-seismic calibration, thereby reducing the number of subsequent local adjustments for the depth correspondence between the logging data and the seismic data. The logging depth-seismic depth corresponding relationship established based on the method provided herein can be employed in the building of a more accurate low-frequency (below 10 Hz) geological background model, which is further used for the identification of geological characteristics and reservoir properties of oil and gas reservoirs, and calculation of oil and gas reserves. Moreover, based on the logging depth-seismic depth corresponding relationship, the logging data and the seismic section image can be combined to locate the sand bodies (e.g., sandstone) in terms of depth and profile.

[0093] The above-mentioned embodiments are only illustrative of this application, and are not intended to limit this application. It should be understood that various variations, modifications and improvements made by those skilled in the art without departing from the spirit of this application shall fall within the scope of the present application defined by the appended claims.