Method for characterizing gas-bearing reservoir based on logging constraint

Abstract

The invention relates to petroleum seismic exploration, and more specifically to a method for gas-bearing reservoir characterization using logging information. In the method, logging information is used as a constraint to indirectly characterize the distribution range of the gas-bearing reservoir by determining the upper and lower boundaries. In addition, this method enables the automatic determination of the optimal calculation parameters according to the characteristics of input data, allowing for more accurate results.

Claims

1. A method for gas-bearing reservoir characterization based on logging constraint, comprising: (1) inputting a seismic record of a 3D seismic data set, and smoothing two ends of the seismic record to obtain a smoothed seismic record x(t); (2) substituting the smoothed seismic record x(t) to $\mathrm{FB}={\left[\frac{{\mathrm{SP}}_{\mathrm{opt}}}{{}_{\mathrm{opt}}}\right]}^2;$ to obtain an optimal fluid mobility interface curve FB; wherein .sub.opt is optimal frequency parameter, and SP.sub.opt is energy of optimal instantaneous amplitude spectrum; calculating SP.sub.opt according to
SP.sub.opt=[real(S.sub.opt)].sup.2+[imag(S.sub.opt)].sup.2; wherein S.sub.opt is optimal instantaneous amplitude spectrum; real(g) represents the real part of S.sub.opt, imag(g) represents the imaginary part of S.sub.opt; calculating S.sub.opt according to $S_{\mathrm{opt}}=\mathrm{ift}\left[X\exp \left(-\frac{2{}^2{f}^2}{{}_{\mathrm{opt}}^2{}_{\mathrm{opt}}^{2{}_{\mathrm{opt}}}}\right)\right];$ wherein X is derived from a Fourier transform of x(t); f is a vector constructed by index numbers of sampling points; .sub.opt and .sub.opt are optimal adjustment parameters, ift(g) represents an inverse Fourier transform; setting a one-dimension vector =[4, 5,L,.sub.u], wherein .sub.u is an upper limit frequency; setting a one-dimension vector =[0.5,0.6,L,1.5]; and setting subscripts of elements in the vector to be i=1, 2,L,11; and determining optimal parameters .sub.opt, .sub.opt and .sub.opt according to the following steps: a) for (i), setting subscripts of elements in the vector to be j=1,2,L,n, and establishing an objective function T((i), (j))=(f.sub.1+f.sub.2) for each element (j) of ; wherein f.sub.1=sum(fb.sub.n), f.sub.2=sum(g), sum(g) represents a summation operation; finding out , which maximizes the value of objective function T((i), (j)), by using a one-dimension grid search method; storing into (j) of a one-dimensional vector , and storing the value of the objective function into T(j) of a one-dimension vector T; b) finding out a maximum value T(maxid) of T, and finding out (maxid) and (maxid) at a location of subscript (maxid); and storing T(maxid), (maxid), (i) and (maxid) into an i-th row of a two-dimension vector with eleven rows and four columns; c) setting i=i+1, if i>12, performing step d), otherwise repeating step a) and step b) again; d) finding out a maximum value in a first column of to obtain last three values in a row corresponding to a subscript of the maximum value in the first column of , wherein the obtained last three values are the optimal frequency parameter .sub.opt and the optimal adjustment parameters .sub.opt and .sub.opt, respectively; in the objective function, setting fb.sub.n as a normalization of the fluid mobility interface curve fb, and calculating fb according to $\mathrm{fb}={\left[\frac{\mathrm{SP}}{\left(i\right)}\right]}^2;$ wherein SP is energy of instantaneous amplitude spectrum; calculating SP according to
SP=[real(S)].sup.2+[imag(S)].sup.2; wherein S is instantaneous amplitude spectrum, and is calculated according to $S=\mathrm{ift}\left[X\exp \left(-\frac{2{}^2{f}^2}{\left(i\right){\left(j\right)}^2}\right)\right];$ in the objective function, calculating g according to g=fb.sub.ngI.sub.G; wherein the symbol g represents a Hadamard product operation; determining a threshold a according to characteristics of input seismic data and logging information, and calculating I.sub.G according to $\{\begin{array}{c}{I}_G\left(I_G<a\right)=0\\ I_G\left(I_{G}a\right)=1\end{array}\mathrm{or}\{\begin{array}{c}{I}_G\left(I_{G}a\right)=1\\ I_G\left(I_G>a\right)=0\end{array};$ (3) normalizing the optimal fluid mobility interface curve FB to obtain FB.sub.n, and standardizing x(t) to obtain x.sub.z(t); (4) determining an upper boundary point set and a lower boundary point set of a gas-bearing reservoir according to extrema on the curve x.sub.z(t) and FB.sub.n, and then determining indication results of the gas-bearing reservoir by using the upper boundary point set and the lower boundary point set according to the following steps: a) determining other threshold b according to information about gas saturation in the logging information, picking up extrema bigger than the threshold b on the curve FB.sub.n; and forming a contrast point set CP through index numbers of these extrema on the curve FB.sub.n; b) standardizing the x(t) to obtain x.sub.z(t), and picking up extrema on the x.sub.z(t); and forming an information point set I through index numbers of these extrema on the curve x.sub.z(t); c) matching the information point set I with the contrast point set CP to form a candidate point set CA; d) selecting a maximum point on the curve FB.sub.n, finding out a reference point corresponding to an index number of the maximum point from the candidate point set CA and taking the reference point as a boundary to divide CA into an upper candidate point set CA.sub.L and a lower candidate point set CA.sub.B, sorting points in the upper candidate point set CA.sub.L in descending order, and sorting points in the lower candidate point set CA.sub.B in ascending order; e) for a point in the upper candidate point set CA.sub.L or in the lower candidate point set CA.sub.B, if its corresponding value in x.sub.z(t) is smaller than an average of x.sub.z(t), putting the point in the CA.sub.L or in the CA.sub.B into the upper boundary point set UB; otherwise if its corresponding value in x.sub.z(t) is equal or bigger than the average of x.sub.z(t), putting the point in the CA.sub.L or in the CA.sub.B into the lower boundary point set LB; wherein there is an one-to-one match between the upper boundary point set UB and the lower boundary point set LB; f) calculating the normalized gas-bearing reservoir indication GI of the input seismic record according to $\mathrm{GI}\left(j\right)=\{\begin{array}{cc}1,& \mathrm{UB}\left(i\right)j\mathrm{LB}\left(i\right)\\ 0,& \mathrm{others}\end{array};$ wherein j is a subscript of an element in GI, and i is a subscript of an element in UB and LB; g) performing inverse proportional transformation on GI to obtain the gas-bearing reservoir indication result of the input seismic record; (5) repeating steps (1) to (4) until all the seismic records in the input 3D seismic data set are completely processed to obtain the gas-bearing reservoir indication results of the whole 3D seismic data set.

2. The method of claim 1, wherein in the use of the logging information as constraints, a distribution range of the gas-bearing reservoir is indirectly characterized by determining upper and lower boundaries of the gas-bearing reservoir, and optimal calculation parameters are automatically determined according to the characteristics of the input data.

3. The method of claim 1, wherein a plurality of wells are used as constraints to obtain a plurality of sets of the optimal frequency .sub.opt, and the optimal adjustment parameters .sub.opt and .sub.opt, and each set of the optimal parameters is used to calculate the gas-bearing reservoir indication result of the input seismic data set, and the respective calculation results are superimposed to obtain a more accurate result.

Description

BRIEF DESCRIPTION OF THE DRAWINGS

(1) FIG. 1 is a section passing through well A of high gas saturation reservoir indication result in the target area according to an embodiment of the invention; where the abscissa indicates trace number; and the ordinate indicates time in second (i.e. s);

(2) FIGS. 2A-2D show the comparison between the normalized high gas saturation reservoir indication result and the logging information of well A;

(3) FIG. 2A shows five seismic traces nearby well A, where a borehole-side trace of well A is marked by arrow, the abscissa indicates trace number, the ordinate indicates time(s);

(4) FIG. 2B shows a logging gas saturation curve, where the abscissa indicates logging gas saturation (%), the ordinate indicates time(s);

(5) FIG. 2C shows a normalized high gas saturation reservoir indication result (solid line) and the optimal fluid mobility interface curve (dash line), where the ordinate indicates time(s);

(6) FIG. 2D shows a logging interpretation curve, where the abscissa indicates the numerical codes, the ordinate indicates time(s);

(7) FIG. 3 is another section passing through well B of high gas saturation reservoir indication result in the same target area according to an embodiment of the invention; where the abscissa indicates trace number; the ordinate indicates time(s);

(8) FIGS. 4A-4D show the comparison between a normalized high gas saturation reservoir indication result and the logging information of well B;

(9) FIG. 4A shows five seismic traces near well B, where a borehole-side trace of well B is marked by an arrow, the abscissa indicates trace number; the ordinate indicates time(s);

(10) FIG. 4B shows a logging gas saturation curve, where the abscissa indicates logging gas saturation (%), the ordinate indicates time(s);

(11) FIG. 4C shows a normalized high gas saturation reservoir indication result (solid line) and the optimal fluid mobility interface curve (broken line), where the ordinate indicates time(s);

(12) FIG. 4D shows a logging interpretation curve, where the abscissa indicates numerical code, the ordinate indicates time(s);

(13) FIG. 5 is a horizontal slice extracted at the target interval, which shows a 3D high-gas saturation reservoir indication result in the target area.

DETAILED DESCRIPTION OF EMBODIMENTS

(14) The method for gas-bearing reservoir characterization based on logging constraint according to the present application is now described below.

(15) (1) inputting a seismic record of a 3D seismic data set, and smoothing two ends of the seismic record to obtain the smoothed seismic record x(t);

(16) (2) substituting the smoothed seismic record x(t) to the following formula to obtain the optimal fluid mobility interface curve FB:

(17) $\mathrm{FB}={\left[\frac{{\mathrm{SP}}_{\mathrm{opt}}}{{}_{\mathrm{opt}}}\right]}^2;$

(18) where .sub.opt is optimal frequency parameter, and SP.sub.opt is energy of optimal instantaneous amplitude spectrum; SP.sub.opt is calculated according to the following formula:
SP.sub.opt=[real(S.sub.opt)].sup.2+[imag(S.sub.opt)].sup.2;

(19) where S.sub.opt is optimal instantaneous amplitude spectrum; real(g) represents the real part of S.sub.opt, imag(g) represents the imaginary part of S.sub.opt;

(20) S.sub.opt is calculated according to the following formula:

(21) $S_{\mathrm{opt}}=\mathrm{ift}\left[X\exp \left(-\frac{2{}^2{f}^2}{{}_{\mathrm{opt}}^2{}_{\mathrm{opt}}^{2{}_{\mathrm{opt}}}}\right)\right];$

(22) where X is derived from Fourier transform of x(t); f is a vector constructed by the index numbers of sampling points; .sub.opt and .sub.opt are optimal adjustment parameters, ift(g) represents the inverse Fourier transform;

(23) based on a one-dimension vector =[4, 5,L,.sub.u], where .sub.u, is an upper limit frequency; and a one-dimension vector =[0.5, 0.6,L,1.5]; and setting the subscripts of the elements in the vector to be i=1,2,L,11, optimal parameters .sub.opt, .sub.opt and .sub.opt are determined according to the following steps:

(24) a) for (i), setting the subscripts of the elements in the vector to be j=1, 2,L,n, and establishing the following objective function for each element (j) of :
T((i),(j))=(f.sub.1+f.sub.2)

(25) where f.sub.1=sum(fb.sub.n), f.sub.2=sum(g), sum(g) represents a summation operation;

(26) finding out , which maximizes the value of objective function T((i), (j)), by using a one-dimension grid search method; storing into (j) of a one-dimensional vector , and storing the corresponding value of the objective function into T(j) of a one-dimension vector T;

(27) b) finding out the maximum value T(maxid) of T, and finding out the corresponding (maxid) and (maxid) at the location of subscript (maxid); and storing T(maxid), (maxid), (i) and (maxid) into the i-th row of a two-dimension vector with eleven rows and four columns;

(28) c) setting i=i+1, if i>12, performing step d), otherwise repeating step a) and step b) again;

(29) d) finding out the maximum value in the first column of , and the last three values in the row corresponding to the subscript of the maximum value are the optimal frequency parameter .sub.opt and the optimal adjustment parameters .sub.opt and .sub.opt, respectively;

(30) in the objective function, fb.sub.n is the normalization of the fluid mobility interface curve fb, and fb is calculated according to the following formula:

(31) $\mathrm{fb}={\left[\frac{\mathrm{SP}}{\left(i\right)}\right]}^2;$

(32) where SP is energy of instantaneous amplitude spectrum, which is calculated according to the following formula:
SP=[real(S)].sup.2+[imag(S)].sup.2;

(33) where S is instantaneous amplitude spectrum, which is calculated according to the following formula:

(34) 0$S_{\mathrm{opt}}=\mathrm{ift}\left[X\exp \left(-\frac{2{}^2{f}^2}{\left(i\right){\left(j\right)}^2}\right)\right];$

(35) in the objective function, g is calculated as follow:
g=fb.sub.ngI.sub.G;

(36) where the symbol g represents the Hadamard product operation;

(37) a threshold a is determined according to characteristics of input seismic data and logging information, and I.sub.G is calculated as the following formula:

(38) $\{\begin{array}{c}{I}_G\left(I_G<a\right)=0\\ I_G\left(I_{G}a\right)=1\end{array}\mathrm{or}\{\begin{array}{c}{I}_G\left(I_{G}a\right)=1\\ I_G\left(I_G>a\right)=0\end{array}$

(39) (3) normalizing the optimal fluid mobility interface curve FB to obtain FB.sub.n, and standardizing x(t) to obtain x.sub.z(t);

(40) (4) determining an upper boundary point set and a lower boundary point set of the gas-bearing reservoir according to extrema on the curve x.sub.z(t) and FB.sub.n, and then determining indication results of the gas-bearing reservoir by using the upper boundary point set and the lower boundary point set according to the following steps:

(41) a) determining other threshold b according to the information about gas saturation in the logging information, picking up extrema bigger than the threshold b on the curve FB.sub.n, and the corresponding index numbers of these extrema on the curve FB.sub.n form a contrast point set CP;

(42) b) standardizing the x(t) to obtain x.sub.z(t), and picking up extrema on the x.sub.z(t), and the corresponding index numbers of these extrema on the curve x.sub.z(t) form an information point set I;

(43) c) matching the information point set I with the contrast point set CP to form the candidate point set CA;

(44) d) selecting the maximum point on the curve FB.sub.n, finding out a point corresponding to the index number of the maximum point from the candidate point set CA and taking the point as a reference point, taking the reference point as a boundary to divide CA into an upper candidate point set CA.sub.L and a lower candidate point set CA.sub.B, sorting points in the upper candidate point set CA.sub.L in descending order, and sorting points in the lower candidate point set CA.sub.B in ascending order;

(45) e) for the point in the upper candidate point set CA.sub.L and the lower candidate point set CA.sub.B, if its corresponding value in x.sub.z(t) is smaller than the average of x.sub.z(t), putting the point into the upper boundary point set UB; otherwise if its corresponding value in x.sub.z(t) equal or bigger than the average of x.sub.z(t), putting the point into the lower boundary point set LB; finally, there is an one-to-one match between the upper boundary point set UB and the lower boundary point set LB;

(46) f) calculating the normalized gas-bearing reservoir indication GI of the input seismic record according to the following formula:

(47) $\mathrm{GI}\left(j\right)=\{\begin{array}{cc}1,& \mathrm{UB}\left(i\right)j\mathrm{LB}\left(i\right)\\ 0,& \mathrm{others}\end{array};$

(48) where j is the subscript of an element in GI, and i is the subscript of an element in UB and LB;

(49) g) performing inverse proportional transformation on GI to obtain the gas-bearing reservoir indication result of the input seismic record;

(50) (5) repeating steps (1) to (4) until all the seismic records in the input 3D seismic data set are completely processed to obtain the gas-bearing reservoir indication results of the whole 3D seismic data set.

(51) a plurality of wells are used as constraints to obtain a plurality of sets of the optimal frequency parameters .sub.opt, and the optimal adjustment parameters .sub.opt and .sub.opt each set of the optimal parameters is used to calculate the gas-bearing reservoir indication results of the input seismic data set, and the respective calculation results are superimposed to obtain a more accurate result.

(52) FIG. 1 is a section passing through well A of high gas saturation reservoir indication result in the target area according to an embodiment of the invention, while there are the optimal frequency .sub.opt=16 Hz, and the optimal adjusting parameters .sub.opt=0.7 and .sub.opt=21.11. In this case, the logging gas saturation and the interpretation of well A are used as constraints to obtain the optimal parameters .sub.opt, .sub.opt and .sub.opt. In the target area, the main object is to find out the distribution range of subsurface high gas saturation reservoirs. Although low gas saturation reservoirs contain gas, they are not profitable for the hydrocarbon exploration in deep water.

(53) FIGS. 2A-D shows the comparison between the normalized indication result of the high gas saturation reservoir and the logging information at well A, while there are the optimal frequency .sub.opt=16 Hz, the optimal adjusting parameters .sub.opt=0.7 and .sub.opt=21.11. In the logging interpretation curve of this area, the numerical code of the high gas saturation reservoir is defined as 11 and the numerical code of the low gas saturation reservoir is defined as 12. The upper limit frequency .sub.u in step (2) of the embodiment is set to 100 Hz. According to the logging interpretation curve, the threshold a is set to 12, and I.sub.G is calculated according to

(54) $\{\begin{array}{c}{I}_G\left(I_{G}a\right)=1\\ I_G\left(I_G>a\right)=0\end{array}.$
The threshold b in step (6) of this embodiment is set to 0.4. It can be seen in the FIG. 2 that the normalized indication result (the solid line in FIG. 2C) of high gas saturation reservoir is highly consistent with the logging gas saturation curve (FIG. 2B) and the logging interpretation curve (FIG. 2D).

(55) FIG. 3 is another section passing through well B of the indication result of high gas saturation reservoir according to an embodiment of the invention, while there are optimal frequency .sub.opt=16 Hz and the optimal adjusting parameters .sub.opt=0.7 and .sub.opt=21.11. In this case, the logging gas saturation and the interpretation of well A are used as constraints to obtain the optimal parameters .sub.opt, .sub.opt and .sub.opt. Well B is used as a reference.

(56) FIGS. 4A-D shows the comparison between the normalized indication result of high gas saturation reservoir and the logging information at well B, while there are optimal frequency .sub.opt=16 Hz, optimal adjusting parameters .sub.opt=0.7 and .sub.opt=21.11. It also can be seen that the normalized indication result (the solid line in FIG. 4C) of high gas saturation reservoir is highly consistent with both the logging gas saturation curve (FIG. 4B) and the logging interpretation curve (FIG. 4D) Thus, the result verifies the accuracy of the invention for the prediction of the distribution range of high gas saturation reservoirs.

(57) FIG. 5 is a horizontal slice extracted at the target interval, which shows a 3D high-gas saturation reservoir indication result in the target area. It can be seen that well A and well B are both located in high gas saturation zones. More importantly, in undrilling zones, the slice may provide reference for determining new drilling locations for profitable wells.

(58) The advantages of the invention are described as follows. 1) The logging information is used as constraints, so that the optimal frequency .sub.opt and the optimal adjustment parameters .sub.opt and .sub.pt can be automatically determined according to the characteristics of the input data. Thus it does not need to manually adjust the parameters. 2) The distribution range of the gas-bearing reservoirs is indirectly obtained by determining the upper and lower boundaries of the gas-bearing reservoirs, which contributes to a more accurate result. 3) The invention has a high computational efficiency, and thus it can be applied to the large-scale 3D seismic data set.

(59) The above embodiments are merely illustrative of the invention, and various variations can be made to the steps in the above embodiments. Any equivalent changes and improvements made without departing from the spirit of the invention should fall within the scope of the invention.