Well-to-seismic precise calibration method based on step-by-step matching optimization

Abstract

The application discloses a well-seismic precise calibration method based on step-by-step matching optimization, inputs seismic data and well logging data, extracts the nearest seismic data of a well site as a well-side seismic trace participating in calibration; for a well to be calibrated, pre-processes well logging data and calculates reflection coefficients; extracts a seismic wavelet and the reflection coefficients to perform convolution, and produces a synthetic seismic record of the well in a time domain; calculates the correlation of the well-side seismic trace and the synthetic seismic record, finds a global time shift amount of the two, and preliminarily aligns the synthetic seismic record by using the time shift amount; calculates a local time shift amount of the well-side seismic trace and the synthetic seismic record, and precisely calibrates the synthetic seismic record by using the time shift amount. The application innovatively proposes a step-by-step matching and point-by-point optimization idea, uses the global time shift amount to constrain a preliminary alignment range of signals, and realizes high-precision alignment of two signals by using the local time shift amount on the basis, so that the method has high precision, strong noise resistance and high efficiency.

Description

Technical Field

[0001] This invention provides a method for accurate well-seismic calibration based on step-by-step matching optimization, which belongs to the field of oil and gas exploration. Background Technology

[0002] In the field of oil and gas exploration, time-domain seismic data is currently more commonly used. However, since well logging data is depth-domain, it is necessary to establish the relationship between time-domain seismic data and depth-domain well logging data through well-seismic calibration, so as to impart geological and well logging information to the seismic data. Therefore, well-seismic calibration is also one of the key steps in subsequent seismic data interpretation and reservoir prediction.

[0003] Well-seismic calibration first utilizes velocity and density information from well logging data to create synthetic seismic records, which are then shifted, aligned, stretched, and compressed. Currently, well-seismic calibration is divided into manual and automatic calibration. Manual calibration introduces the subjective factors of the interpreter, with two main drawbacks: firstly, the accuracy of the calibration heavily relies on the interpreter's experience; secondly, many oilfields have a large number of wells logged, making manual calibration inefficient. These factors severely restrict the efficiency and accuracy of well-seismic calibration. Automatic calibration calculates the time shift between the synthetic seismic record and the well-side seismic trace using algorithms, and then stretches and compresses the synthetic seismic record based on this time shift. Commonly used automatic calibration algorithms fall into two main categories: one is correlation-based, which finds the time shift by calculating the maximum correlation value between the two sequences. The drawback of this method is that it requires manual selection of the correlation window size, leading to inaccurate time shift calculations. The second method is based on dynamic time warping, which calculates the cumulative Euclidean distance matrix of well seismic records through dynamic programming to find the optimal time-shift path. The disadvantages of this method are that it requires manual alignment of the well seismic records, a given time window size for calculation, and is susceptible to noise interference. Summary of the Invention

[0004] To address the shortcomings of the aforementioned well-seismic calibration methods, this invention proposes a step-by-step matching optimization-based precise well-seismic calibration method. This method aims to first find an optimal global time-shift to perform a preliminary matching between the synthetic seismic record and the well-side seismic trace, and then calculate the local fine-grained time-shift to achieve high-precision automatic calibration. This method not only improves the accuracy and stability of the calibration but also increases the efficiency while maintaining accuracy.

[0005] The specific technical solution of this invention is as follows:

[0006] A well-seismic accurate calibration method based on step-by-step matching optimization includes the following steps:

[0007] S1: Input seismic data and well logging data, extract the seismic data closest to the well location coordinates as the wellside seismic trace for calibration, and perform standardization processing on the seismic data trace.

[0008] S2: For the well to be calibrated, preprocess its logging data and calculate the reflection coefficient sequence.

[0009] S3: Extract the seismic wavelet and reflection coefficient sequence from the seismic trace near the well, convolve them to create a synthetic seismic record of the well in the time domain, and then standardize the synthetic seismic record.

[0010] S4: Calculate the correlation between well-side seismic traces and synthetic seismic records using a high-precision cross-correlation scanning algorithm, find the global time shift between synthetic seismic records and well-side seismic traces, and use this global time shift to perform preliminary alignment of the synthetic seismic records.

[0011] S5: Calculate the local time shift of the well-side seismic trace and the synthetic seismic record using a path-constrained dynamic time warping algorithm, and use this local time shift to perform precise calibration of the synthetic seismic record.

[0012] Preferably, S1 includes the following sub-steps:

[0013] S11: Input data, including time-domain seismic data and seismic interpretation horizons, of which well logging data include sonic transit time logging data, density logging data, resistivity logging data, well location coordinates, and geological strata.

[0014] S12: Calculate the Euclidean distance based on the well location coordinates and the coordinates of nearby seismic traces, extract the seismic trace with the smallest Euclidean distance as the well-side seismic trace to participate in the calibration, denoted as seis(t), and standardize the trace, denoted as Seis(t).

[0015] Preferably, S2 includes the following sub-steps:

[0016] S21: Preprocess the sonic logging data. Sonic logging data may contain outliers due to environmental influences. Replace the outlier segments in the sonic logging data using the following formula:

[0017]

[0018] Wherein, DT(d) represents the sonic logging data before preprocessing; R(d) represents the resistivity logging curve; This represents the reconstructed sonic logging data; d is the depth value; m is a constant.

[0019] S22: Due to changes in well diameter, density logging data showed anomalies. The Gardner formula was used to replace the outliers in the density logging curves. The reconstructed calculation formula is as follows:

[0020]

[0021] in This represents the reconstructed density logging curve;

[0022] S23: The moving average digital filtering method is used to eliminate high-frequency noise in the acoustic time-of-flight logging curve. The specific algorithm formula is as follows:

[0023]

[0024] in, This represents the sonic logging data after smoothing by moving average filtering; n is a constant representing the number of points involved in the filtering.

[0025] S24: Calculate the P-wave velocity sequence using the processed sonic logging data:

[0026]

[0027] Vp(d) represents the longitudinal wave velocity calculated using sonic logging data.

[0028] S25: The depth domain reflection coefficient sequence Ref(d) is calculated using the P-wave velocity sequence and density logging curve.

[0029] Preferably, S3 includes the following sub-steps:

[0030] S31: Calculate the correspondence between the depth domain and the time domain using sonic logging data, and convert the reflection coefficient sequence in the depth domain to the time domain, denoted as Ref(t).

[0031] S32: Using the reflection coefficient sequence combined with well-side seismic trace records, the seismic wavelet is obtained from the convolution model. Based on the least squares principle, the error energy is established, and the seismic wavelet ω(t) that minimizes the error energy is obtained.

[0032]

[0033] Where ε represents the error energy; ω(t) represents the obtained seismic wavelet; Ref(t) represents the time-domain reflection coefficient sequence; M represents the total number of sampling points; * represents the convolution operation; seis(t) represents the unnormalized well-side seismic trace;

[0034] S33: The synthesized seismic record is obtained by performing a convolution operation on the obtained seismic wavelet and the time-domain reflection coefficient sequence, and then the synthesized seismic record is standardized.

[0035] Syn(t) = Ref(t) * ω(t)

[0036] Syn(t) represents the synthetic seismic record.

[0037] Preferably, S4 includes the following sub-steps:

[0038] S41: The total number of sampling points for both the well-side seismic trace Seis(t) and the synthetic seismic record Syn(t) is M. For each sampling point, a Gaussian window with a time window size of N and a standard deviation of σ is used to truncate the synthetic seismic record Syn(t) and the well-side seismic trace Seis(t) to generate a truncated signal. and

[0039]

[0040]

[0041] in, This indicates the truncated signal after cutting off the well-side seismic trace Seis(t) using a Gaussian time window; σ represents the truncated signal of the synthesized seismic trace Syn(t) using a Gaussian time window; σ represents the standard deviation; t represents the time value, ranging from [1, M]; k represents the sampling point of the time window, ranging from...

[0042] S42: Convert the geological stratification in the depth domain to the time domain, calculate the time deviation between the geological stratification and the seismic horizon, and take the absolute value, denoted as U.

[0043] S43: Calculation of truncated synthetic seismic records at different time shifts Seismic channel next to the cut-off well The similarity is calculated as follows:

[0044]

[0045] Among them, cor t，p This represents the truncated signal formed by cutting off the data at the t-th sampling point when the drift is p. With truncated signal The correlation; p represents the time shift, with a value range of [-2U, 2U];

[0046] S44: From cor t，p Find the relevant maximum value in the range and record its corresponding p as the global time shift p1:

[0047] p1 = argmax p (cor t，p )

[0048] Where argmax p This indicates that the function cor t，p The value of variable p when it reaches its maximum value.

[0049] S45: Initial alignment of the synthetic seismic record is performed using the global time shift p1 to form the synthetic seismic record. The formula is:

[0050]

[0051] in, This represents the preliminary aligned synthetic seismic record.

[0052] Preferably, S5 includes the following sub-steps:

[0053] S51: Construct an energy error matrix using preliminarily aligned synthetic seismic records and well-side seismic traces;

[0054] E l，t =|Seis(t+l)-Syn(t)| 2 ;

[0055] Where E represents the energy error matrix constructed from the well-side seismic traces and synthetic seismic records. This error matrix is ​​a two-dimensional array with 2L+1 rows and M columns, i.e., l = -L, ..., -1, 0, 1, ..., L, t = 1, ..., T; l，t This represents the value in row l, column t.

[0056] S52: Utilize path-constrained dynamic programming to achieve optimal iterative accumulation of the energy error matrix E, establishing the accumulated energy error matrix D, as follows: Figure 5 As shown, the algorithm is as follows:

[0057]

[0058] Where D represents the cumulative energy error matrix established by the well-side seismic traces and synthetic seismic data. This error table is a two-dimensional array with 2L+1 rows and T columns, i.e., l = -L, ..., -1, 0, 1, ..., L, t = 1, ..., T; D l，t This represents the value in row l, column t;

[0059] S53: Perform a path scan on the accumulated energy error matrix to find an optimal path q(t), where q(t) is a local time shift that stores the time shift of each sampling point in the synthetic seismic record. The path scan algorithm is as follows:

[0060] q(t)=argmin l (D l，t )

[0061] Where, q t q represents the sequence storing time shifts. t ={q1, ...,q M};argmin l This represents the value of variable l that makes the function reach its minimum value.

[0062] S54: Calculate new synthetic records using local time shifts The formula is:

[0063]

[0064] This invention solves the problems of low accuracy and poor stability in current well seismic calibration, and improves calibration efficiency while ensuring the accuracy of automatic calibration. Attached Figure Description

[0065] Figure 1 This is a flowchart of the present invention;

[0066] Figure 2 The synthetic seismic records and well-side seismic traces to be calibrated in this embodiment;

[0067] Figure 3 In this embodiment, a high-precision cross-correlation scanning algorithm is used to calculate the correlation between the synthetic seismic record and the seismic trace near the well.

[0068] Figure 4 This is a comparison diagram with the seismic traces near the well after the synthetic seismic record is initially aligned using global time shift in this embodiment;

[0069] Figure 5 The cumulative energy error matrix established between the initially aligned synthetic seismic record and the well-side seismic trace in this embodiment, and the local time shift q sought;

[0070] Figure 6 This is a comparison diagram with the well-side seismic trace, showing the calibration of the preliminarily aligned synthetic seismic record using local time shifts in this embodiment. Detailed Implementation

[0071] The specific technical solutions of the present invention will be described with reference to the embodiments.

[0072] The present invention will be further described below with reference to the accompanying drawings and specific embodiments:

[0073] Taking a well in an actual work area as an example, this paper provides a method for accurate well seismic calibration based on step-by-step matching optimization, such as... Figure 1 As shown, it consists of five steps:

[0074] The first step is to input seismic data and well logging data, extract the seismic traces near the wells, and complete the standardization, denoted as Seis(t). Figure 2 As shown by the dashed line.

[0075] The second step involves inputting acoustic, density, and resistivity logging data, completing the preprocessing of the acoustic and density logging data, and calculating the depth-domain reflection coefficient sequence of the well.

[0076] The third step involves converting the depth-domain reflection sequence to the time domain, extracting seismic wavelets from the well-side seismic traces, and convolving them with the time-domain reflection coefficients to complete the synthesis of the seismic record and its standardization, denoted as Syn(t). Figure 2 As shown by the solid line in the middle.

[0077] The fourth step involves automatically calculating the absolute deviation U between geological strata and seismic horizons, and then using a high-precision cross-correlation scanning algorithm to calculate the correlation between well-side seismic traces and synthetic seismic records. Figure 3 As shown, the global time shift p1 is extracted and used for initial alignment to form a preliminary aligned synthetic seismic record. like Figure 4 As shown.

[0078] Step 5: Input the preliminarily aligned synthetic seismic records. The local time shift q is calculated using a path-constrained dynamic time warping algorithm based on the seismic trace Seis(t) near the well. Figure 5 As shown, and aligned with q, as Figure 6 As shown.

Claims

1. A well-to-seismic precise calibration method based on step-by-step matching optimization, characterized in that, It comprises the following steps: S1: input seismic data and logging data, extract the seismic data closest to the well coordinate as the well seismic trace for calibration, and standardize the seismic data; S2: preprocess the logging data of the well to be calibrated, and calculate the reflection coefficient sequence; S3: extract the seismic wavelet of the well seismic trace and the reflection coefficient sequence, make a synthetic seismogram of the well in the time domain, and standardize the synthetic seismogram; S4: calculate the correlation between the well seismic trace and the synthetic seismogram using a high-precision cross-correlation scanning algorithm, find the global time shift between the synthetic seismogram and the well seismic trace, and use the global time shift to perform preliminary alignment on the synthetic seismogram; S5: calculate the local time shift between the well seismic trace and the synthetic seismogram using a path-restricted dynamic time warping algorithm, and use the local time shift to perform accurate calibration on the synthetic seismogram.

2. The well-to-seismic precise calibration method based on step-by-step matching optimization according to claim 1, characterized in that, S1 comprises the following sub-steps: S11: input data, including time-domain seismic data, seismic interpretation horizon, and well coordinate, and well logging data including acoustic travel time logging data, density logging data, resistivity logging data, well coordinate, and geological layering; S12: calculate the Euclidean distance between the well coordinate and the nearby seismic trace coordinate, extract the seismic trace with the minimum Euclidean distance as the well seismic trace for calibration, and standardize the seismic trace, denoted as S(t).

3. The well-to-seismic precise calibration method based on step-by-step matching optimization according to claim 2, characterized in that, S2 comprises the following sub-steps: S21: preprocess the acoustic logging data, which may be affected by the environment and have data outliers, replace the abnormal section of the acoustic logging data, and the replacement formula is: Wherein, DT(d) represents the acoustic logging data before pretreatment; R(d) represents the resistivity logging curve; represents the reconstructed acoustic logging data; d is the depth value; m is a constant; S22: due to the change of well diameter, the density logging data has abnormal values, Gardner formula is used to replace the abnormal values in the density logging curve, and the reconstruction calculation formula is: wherein denotes the reconstructed density log curve; S23: use sliding mean digital filtering method to eliminate high-frequency noise in the acoustic travel time logging curve, and the specific algorithm formula is as follows: wherein, S represents the acoustic logging data after smoothing by sliding average filtering; n is a constant, representing the number of points participating in filtering; S24: calculate the P-wave velocity sequence using the processed acoustic logging data: Vp(d) represents the P-wave velocity calculated using the acoustic logging data; S25: calculate the depth domain reflection coefficient sequence Ref(d) using the P-wave velocity sequence and the density logging curve.

4. The well-to-seismic precise calibration method based on step-by-step matching optimization according to claim 3, characterized in that, S3 comprises the following sub-steps: S31: calculate the correspondence between the depth domain and the time domain using the acoustic logging data, convert the reflection coefficient sequence in the depth domain to the time domain, denoted as Ref(t); S32: use the reflection coefficient sequence and the well seismic trace record to obtain the seismic wavelet from the convolution model, establish the error energy according to the least squares principle, and obtain the seismic wavelet ω(t) that minimizes the error energy: Where ε represents the error energy; ω(t) represents the obtained seismic wavelet; Ref(t) represents the reflection coefficient sequence in the time domain; M represents the total number of sampling points; * represents convolution operation; seis(t) represents the well seismic trace without standardization; S33: perform convolution operation using the obtained seismic wavelet and the reflection coefficient sequence in the time domain to obtain the synthetic seismogram, and standardize the synthetic seismogram: Syn(t) = Ref(t) * ω(t) wherein Syn(t) represents a synthetic seismic record.

5. The well-to-seismic precise calibration method based on step-by-step matching optimization according to claim 4, characterized in that, S4 comprises the following sub-steps: S41: the total number of sampling points of the normalized wellside seismic trace S(t) and the synthetic seismic record Syn(t) is M, for each sampling point, a Gaussian window with a window size of N and a standard deviation of σ is selected to truncate the synthetic seismic record Syn(t) and the wellside seismic trace S(t) respectively, to generate truncated signals with wherein, represents the truncated signal after truncating the seismic trace beside the well by using the Gaussian time window; represents the truncated signal after truncating the synthetic seismic trace Syn(t) by using the Gaussian time window; σ represents the standard deviation; t represents the time value, and the value range is [1, M]; k represents the sampling point of the time window, and the value range is S42: converting the geological layering in depth domain to time domain, and calculating the time deviation amount of the geological layering from the seismic horizon and taking absolute value, denoted as U; S43: Calculate truncated synthetic seismic records at different time shifts Similarity to the truncated near-borehole seismic traces is calculated as follows: wherein cor t，p represents a correlation between the truncated signal and the truncated signal ; p represents a shift amount, and the value range of p is [-2U, 2U]. S44: Find the correlation maximum from cor t,p and record its corresponding p as the global time shift p1: p1 = argmax p (cor t，p ) where argmax p denotes the value of the variable p for which the function cor t,p takes its maximum value; S45: preliminary alignment of the synthetic seismic record using the global time shift p1, forming a synthetic seismic record The formula is: wherein, represents a preliminary aligned synthetic seismogram.

6. The well-to-seismic precise calibration method based on step-by-step matching optimization according to claim 5, characterized in that, S5 comprises the following sub-steps: S51: constructing an energy error matrix by using the preliminary aligned synthetic seismic record and the well seismic trace; where E represents an energy error matrix constructed from the well- by- seismic trace and the synthetic seismic record, the error matrix being a two- dimensional array of 2L+1 rows and M columns, i.e., l = -L,..., -1, 0, 1,..., L, t = 1, 2,..., M; E l，t represents the value of the lth row, tth column; S52: implementing optimal iterative accumulation on the energy error matrix E by using path-limited dynamic programming, establishing an accumulated energy error matrix D, and the algorithm is as follows: where D represents the cumulative energy error matrix established by the wellside seismic trace and the synthetic seismic, which is a two-dimensional array with 2L+1 rows and M columns, i.e. l=-L,...,-1,0,1,...,L,t=1,2,...,M; D l，t represents the value of the lth row and the tth column. S53: performing path scanning on the accumulated energy error matrix to find an optimal path q(t), q(t) is a local time shift amount, which stores the time shift amount of each sampling point in the synthetic seismic record, and the path scanning algorithm is as follows: q(t) = argmin l (D l，t ) where q t represents a sequence of storage time shifts, q t = {q1,..., q M}; argmin l represents the value of the corresponding variable l at which the function takes its minimum value; S54: Calculate new synthetic record using local time shift The formula is:

Images