US 12,287,441 B1
System for updating clastic rock lithology model by well-to-seismic integration with deep earth oil and gas precision navigation
Fei Tian, Beijing (CN); Qingyun Di, Beijing (CN); Wenhao Zheng, Beijing (CN); Yongyou Yang, Beijing (CN); and Wenjing Cao, Beijing (CN)
Assigned to Institute of Geology and Geophysics, Chinese Academy of Sciences, Beijing (CN)
Filed by Institute of Geology and Geophysics, Chinese Academy of Sciences, Beijing (CN)
Filed on Jan. 14, 2024, as Appl. No. 18/412,597.
Claims priority of application No. 202311353155.9 (CN), filed on Oct. 19, 2023.
Int. Cl. G01V 1/28 (2006.01); G01V 1/40 (2006.01)
CPC G01V 1/282 (2013.01) [G01V 1/40 (2013.01)] 2 Claims
OG exemplary drawing
 
1. A method for updating geological prediction model while drilling, the method is configured to be performed on a system comprising a data acquisition and standardization module, a variance attribute acquisition module, a formation thickness prediction module, the frequency-division seismic data acquisition module, an optimal seismic frequency band acquisition module, a geological prediction model acquisition module and a model updating module, the method comprising:
acquiring historical seismic data and historical logging data by the data acquisition and standardization module; and collecting original seismic data and original logging data while drilling, and acquiring logging data from which outliers are eliminated through an isolated forest algorithm; eliminating the outliers from the logging data to obtain standardized logging data, the standardized logging data comprising standardized gamma ray logging parameters, obtaining the standardized logging data comprising the following steps:
A10: using the original logging data as a data set X to be processed;
A20: randomly selecting φ data points from the data set X to be processed to form a data subset X′ to be processed, and storing in a root node the data subset x′ to be processed;
A30: randomly selecting a dimension q from the data set X to be processed, and randomly generating a break point p in the dimension q, where the break point p satisfies min(xij,j=q,xij∈X′)<p<max(xij,j=q,xij∈X′), j represents a serial number of the break point, and i represents a data category;
A40: generating a hyperplane, which divides data in the dimension q into two subspaces, according to the break point p; and specifying data points of the dimension q with values less than p to be placed into a first leaf node, and data points with values greater than or equal to p to be placed into a second leaf node;
A50: performing recursion on the methods in A30 to A40 until the leaf nodes have only one data point or an isolated tree has reached a preset height;
A60: repeating the methods in A30 to A50 until T isolated trees are generated, wherein the T isolated trees represent: isolated trees having no external nodes of leaf nodes, or having two leaf nodes (Ti, Tr) and one internal node test; the internal node test of the T isolated trees consists of a dimension q and a break point p, and points of q<p belong to Ti, and points of q>p belong to Tr;
A70: forming the T isolated trees into an isolated tree forest, setting each data point xi to traverse each isolated tree, and calculating a height h(xi) of a data point xi in each isolated tree, that is, a number of edges passed by the data point xi from a root node to a leaf node of the isolated tree; calculating an average height of the data point xi in the isolated tree forest, and normalizing the average height of the data point Xi to obtain a normalized average height h(xi) of the data points;
A80: calculating an outlier score s(x, φ) based on the normalized average height h(xi) of the data points:

OG Complex Work Unit Math
where c(φ) represents an average value of path lengths of a binary tree constructed by φ data points, and E(*) represents an expectation;

OG Complex Work Unit Math
where H(i) represents a harmonic number, and is estimated by ln(i)+0.5772156649, 0.5772156649 being a Euler's constant;
eliminating corresponding data points in response to the outlier score s(x, φ) being less than a preset outlier threshold s to obtain logging data C={c1, . . . , cα, . . . , cm} from which the outliers are eliminated, 1≤α≤m, cα∈C, and m represents a number of data points in the logging data from which the outliers are eliminated;
A90: continuously standardizing the logging data from which the outliers are eliminated to obtain the standardized logging data:

OG Complex Work Unit Math
where gz represents a data value of a zth sampling point of the logging data from which the outliers are eliminated, and Average represents a calculated average value; g represents all data sample points of the logging data from which the outliers are eliminated; v represents a variance of the logging data from which the outliers are eliminated; cz represents a data value of a zth sampling point of the standardized logging data from which the outliers are eliminated;
extracting variance attribute characteristics of standardized gamma ray logging parameters by the variance attribute acquisition module to obtain a variance attribute curve while drilling;
the variance attribute acquisition module being specifically configured to:
set a time window size w of a curve sequence with the variance attribute characteristics of the standardized gamma ray logging parameters according to a target formation thickness, and calculate a variance value var of the standardized gamma ray logging parameters within the time window range of the curve sequence:

OG Complex Work Unit Math
where n represents a GR curve sampling point, GRn represents a value of an nth GR sampling point, and w represents the time window size; and
calculate each sampling point to obtain the variance attribute curve while drilling;
acquiring by the formation thickness prediction module, a predicted formation thickness value of each set based on the variance attribute curve while drilling;
the formation thickness prediction module being specifically configured to:
take the derivative of the variance attribute curve while drilling to acquire a transient attribute of the variance while drilling;
extract maximum points in the transient attribute of the variance while drilling to acquire a mutation point detection result;
input new sampling points, repeatedly acquire transient attributes of the variance while drilling to obtain a plurality of mutation point detection results, and determine that different mutation point detection results are of separate sets to obtain real-time updated set interface depth judgment results; and
count a minimum variance value varmin and a corresponding maximum formation thickness hmax, as well as a maximum variance value varmax and a corresponding minimum formation thickness hmin according to the variance attribute curve while drilling and the real-time updated set interface depth judgment results, establish a linear formation thickness prediction model, and further calculate a predicted formation thickness value hr of each set:

OG Complex Work Unit Math
acquiring by the frequency-division seismic data acquisition module, frequency-division seismic data based on the original seismic data by means of a frequency-division processing method based on wavelet transform;
the frequency-division seismic data acquisition module being specifically configured to:
design a processed wavelet basis function, which matches the original seismic data better, based on a Morlet wavelet basis function:

OG Complex Work Unit Math
where ψg represents the wavelet basis function; t represents time; k represents a variable that satisfies a normalization condition of the processed wavelet basis function; σ represents a time-frequency two-dimensional resolution interval; μ represents a constant that controls a central position of a frequency after the wavelet transform; ω represents an imaginary unit;
allow the processed wavelet basis function to satisfy the normalization condition ∥ψg2=1, and set:

OG Complex Work Unit Math
input all the seismic trace data, and perform frequency-division transform on a seismic trace through the processed wavelet basis function to obtain frequency-division seismic data:

OG Complex Work Unit Math
where a represents a scale factor; b represents a time-shift factor; Wg(a, b) represents an amplitude value at a specific frequency a and a specific time b; fp represents a frequency of a wavelet basis; α and β are distribution parameters in a normal distribution N(β/2α, 1/2α), respectively;
acquiring by the optimal seismic frequency band acquisition module, optimal seismic frequency band information based on the predicted formation thickness value of each set;
the optimal seismic frequency band acquisition module being specifically configured to:
calculate a ratio of the predicted formation thickness value of each set to a wavelength of seismic waves, superpose reflected waves generated by a formation top and a formation bottom to form a strong amplitude signal in response to the ratio of the predicted formation thickness value of each set to the wavelength of seismic waves being ½, and calculate a ratio of a seismic wave velocity to a predicted formation thickness value of the corresponding set to obtain the optimal frequency band information;
establishing by the geological prediction model acquisition module of a current formation, a geological prediction model of the current formation based on the frequency-division seismic data, the optimal seismic frequency band information, historical gamma ray logging parameters in the historical logging data, and the standardized gamma ray logging parameters;
the geological prediction model acquisition module of the current formation being specifically configured to:
determine a formation time window top interface, a main frequency f(Hz) and a complete seismic wave wavelength Δt=1/f(s) according to three-dimensional curved surface data of a target horizon t1(x,y) determined according to the historical seismic data, extend the three-dimensional curved surface data of the target horizon downward by Δt to obtain a formation time window bottom interface, and obtain a formation time window and a formation time window size;
intercept seismic waveform data of an interest interval based on the formation time window, obtain seismic waveform sampling data of the interest interval, resample the seismic waveform sampling data of the interest interval, and obtain seismic waveform resampling data of the interest interval;
cluster the seismic waveform resampling data of the interest interval to obtain a waveform clustering result of each frequency band;
establish a geological prediction model of the current formation by means of waveform indication inversion based on the waveform clustering result of each frequency band, the frequency-division seismic data, the optimal seismic frequency band information, the historical gamma ray logging parameters and the standardized gamma ray logging parameters; and
repeating functions of the data acquisition and standardization module, the variance attribute acquisition module, the formation thickness prediction module, the frequency-division seismic data acquisition module, the optimal seismic frequency band acquisition module, and the geological prediction model acquisition module of the current formation, in response to each drill to a preset depth interval, to generate an updated geological prediction model of the current formation; and
guiding a drilling trajectory through the updated geological prediction model of the current formation.