The Chengdao Oilfield is geographically in Shengli extra-shallow sea area with a water depth of 2-18 m. It is located structurally in the intersection of Bozhong Depression and Jiyang Depression, in the southeast of Chengbei Low bulge of Chengning Uplift. It is adjacent to Chengbei Sag by Chengbei Fault in the west and dips in Bozhong Sag eastward (Fig. 1). The area is an inherited overlapping structure formed on the background of pre-Tertiary buried hill after receiving Tertiary sediments.

The upper member of Neogene Guantao Formation (Upper Guantao Member for short) is one of the main oil-bearing strata in the area, in which the reservoir is typical meandering river deposit featuring shallow burial depth, small thickness of single layers, multi-types of combinations, and fast changes in vertical and horizontal directions^[1]. Practice shows that the key to exploration and development of this kind of oil reservoir is the fine identification and description of thin meandering river reservoirs.

Because the extra-shallow sea is affected by surface tide, silt and sea water reverberation, the seismic data acquired by 3D dual sea and land geophones are characterized by low signal-to-noise (S/N) ratio and low resolution after processed by conventional fusion technique, by which it is difficult to predict thin reservoirs correctly as the seismic event generally reflects the seismic response of an interval of certain thickness made up of thin interbeds or thick sand body. The existing technologies such as horizontal slice, poststack attribute, poststack frequency bandwidth expanding and poststack inversion^[2,3,4,5] can only delineate channel sandbodies more than 10 m thick in this area, but there is no breakthrough in seismic delineation and fluid identification technology for thin reservoirs less than 10 m thick, which restricts the efficient exploration and development of this area.

In view of the characteristics of seismic data and difficulties in fine delineation of thin reservoirs in the Shengli Chengdao extra-shallow sea oilfield, we have worked out a multi-disciplinary technology series to realize the intelligent recognition and delineation of seismic sedimentary microfacies and thin reservoirs, further the semi-quantitative prediction of oil-bearing sandbodies in this area. In this technology series, first, sea and land dual-sensor prestack two-step high resolution frequency bandwidth expanding processing methods, including sea and land dual-sensor consistency fusion processing technology and prestack large angle gather amplitude preservation and frequency bandwidth expanding technology, are used to improve S/N ratio and resolution of seismic data and the ability to identify thin reservoirs of seismic data effectively, and lay foundation for later thin reservoir prediction. Then, the intelligent horizon-facies controlled recognition technology developed by ourselves on the basis of seismic fine horizon constraint is used to establish seismic facies models of three types of meandering river sedimentary microfacies, i.e. floodplain, natural levee and point bar, to realize the intelligent recognition and delineation of seismic sedimentary microfacies and thin reservoirs, so as to improve the delineation accuracy and efficiency of thin reservoirs of meandering river facies. Finally, the prestack multi-parameter fluid probability technique is employed to realize the semi-quantitative prediction of oil-bearing sandbodies in this area. This set of technology provides sound geophysical technical support for exploration and development of thin reservoirs in the oilfield.

At present, seismic data acquisition by dual-sensor Ocean Bottom Cable (OBC) has started in Shengli extra-shallow sea area, Chengdao Oilfield. The seismic data, processed by conventional sea and land dual-sensor fusion method, is characterized by low S/N ratio and low resolution, and cannot meet the requirements of thin reservoir identification and delineation, hindering the exploration and development of this kind of oil reservoir. By strengthening tests of processing parameters and quality control, and following the principle of fidelity and amplitude preservation, a two-step frequency bandwidth expanding processing method coring on sea and land dual- sensor consistency fusion processing and prestack large angle gather frequency bandwidth expanding while preserving amplitude has been worked out, which can improve the S/N ratio and resolution of seismic data greatly.

The conventional sea and land dual-sensor fusion processing technology is based on the assumption that there are the same wave field for water sensor and land sensor. The early conventional fusion just fused directly two types of seismic data or adjust their events for fusion after matching their energy, which commonly resulted in the damage of effective waves and poor suppression of sea water reverberation ^[6,7,8,9,10]. Hence, the seismic data obtained was characterized by low S/N ratio and low resolution, and couldn’t meet the needs of thin reservoir delineation. The consistency fusion processing technology emerging later solves this problem well. In this technology, firstly, the land sensor data is differentiated, and the scale transformation factor is determined according to the energy of the water sensor data, and then the land sensor data after differential is merged with the water sensor data. It can be seen from the seismic sections that the seismic data after sea and land dual-sensor consistency fusion processing has improved significantly in S/N ratio and resolution (Fig. 2).

This method solves the problem of poor consistency of two types of sea and land dual-sensor seismic data. With notch waves compensated, the seismic data is improved in dominant frequency and effective frequency band (Fig. 3), and thus has much higher resolution of sand bodies. The original seismic data before processing often has a dominant frequency of 30 Hz and superior frequency band of 10-50 Hz, and is able to show sand bodies of 10 m thick. The seismic data after processing has a dominant frequency of 40 Hz and superior frequency band of 10-70 Hz, and is able to show sand bodies of 7.5 m thick.

Although the seismic data after sea and land dual- sensor consistency fusion processing has improved in S/N ratio and resolution, it still does not meet the requirements of seismic delineation of thin reservoirs in extra-shallow sea. The poststack frequency bandwidth expanding methods such as conventional anti-Q filtering, frequency expanding and frequency division cannot ensure data fidelity.

Due to the interference tuning effect of prestack large angle common imaging point (CIP) gather, seismic data after stacked in full gather would become narrower in frequency bandwidth and lower in resolution, not conducive to fine delineation of thin fluvial reservoirs. To preserve amplitude and expand frequency bandwidth of seismic data, the technique of prestack large angle gather frequency bandwidth expanding and amplitude preserving has been developed. The main principle is based on data homology, the reflected waves of large angle gathers are stretched and corrected by working out firstly the matching factor between large angle data and small angle data, and then the gather after matched is stacked and imaged to make the final imaging data keeping the same wavelet and frequency bandwidth characteristics of small angle gather^[11]. This technology makes the spectra of data of different angle gathers consistent, and enables in-phase superposition imaging. The imaging result obtained has amplitude well preserved and high resolution, avoids the accumulation effect caused by large angle data interference tuning and the problems caused by traditional methods, such as wavelet deformation, narrow frequency bandwidth, aliasing and low resolution. Finally, the seismic data after processing has frequency band effectively broadened and dominant frequency increased. The seismic data processed by this technology has dominant frequency increased from 40 Hz to 50 Hz, and superior frequency bandwidth of 10-90 Hz, and allow identification of 6 m thick sand bodies, which basically meets the requirement on seismic data resolution in the exploration of thin reservoirs in this oilfield. Taking the seismic section through Well CB 208 as an example, the Layer 2 of the upper of Guantao Formation (with top at 1400 m depth) and the Layer 3 of the upper of Guantao Formation (with the top at 1520 m depth) (Fig. 4a) appear as blank and weak reflections on the original seismic section (Fig. 4b). After amplitude preserving and band expanding processing, the seismic section show medium- strong amplitude reflections at these two layers (Fig. 4c), and has much higher resolution of sand bodies.

Horizon control refers to the establishment of isochronous constraint surface for sedimentary mi-crofacies prediction. Facies control refers to three seismic microfacies of meandering (point bar, flood-plain and natural levee) as constraints for seismic characteristics analysis. Horizon-facies control refers to the joint constraints for intelligent prediction of meandering sedimentary microfacies and sand body delineation.

The intelligent horizon-facies controlled recognition technique was used to identify and delineate thin meandering river sand bodies in Chengdao extra-shallow sea oilfield on the seismic section after sea and land dual-sensor prestack two-step high resolution frequency bandwidth expanding processing. Based on the establishment of isochronal constraint surfaces, the seismic facies characteristics of different sedimentary microfacies were examined and the seismic facies of thin meandering river reservoir were recognized intelligently, which mainly included sample point picking, seismic attribute specification, construction of intelligent sedimentary microfacies recognition sample set, and intelligent prediction of sedimentary microfacies sand bodies.

The Guantao Formation in Chengdao extra-shallow sea oilfield is a long-term base level cycle from bottom to top, bounded by the unconformity between the Guantao Formation and the Paleogene at the bottom and the interface between Guantao Formation and Minghuazhen Formation at the top. Within the long-term base level cycle, the Upper Guantao Member is subdivided into three medium-term base level cycles^[12,13,14], Layer N₁g₁1, N₁g₁2 and N₁g₁3 (Fig. 5a, 5b). The target layer N₁g₁2⁵ is located in the middle cycle of Layer N₁g₁2, which is characterized by alternate intermittent and strong amplitude reflections under the background of weak reflection, representing the deposit of mudstone and isolated channel sand body in floodplain (Fig. 5c)

On the basis of the division of medium-term base level cycles, the sand layers were correlated bed by bed based on marker beds and auxiliary marker beds by considering factors such as vertical characteristics of sedimentary sequence, spatial variation of sedimentary facies and superposition patterns of sandstone. On the basis of calibration of synthetic seismic record, the horizons calibrated in single wells were correlated with horizons on seismic section to ensure consistency between the seismic horizons and geological horizons, and then the isochronal constraint surfaces for intelligent prediction of sedimentary microfacies were established (Fig. 5c) to ensure the interpretation accuracy of the top of N₁g₁2⁵ sand layer.

In view of the meandering river deposits in the Upper Guantao Member of Chengdao extra-shallow sea oilfield, 147 sedimentary microfacies intervals (including 57 point bars, 32 natural levees and 58 floodplains) in Well CB25, CB27, CB252, CBG10 and SHH2 with complete logging data, conventional logging curves, time-depth relationship and sedimentary microfacies interpretation, were selected as the samples for analysis to sort out the seismic facies patterns of floodplain, natural levee, point bar sedimentary microfacies (Fig. 6). The floodplain deposit is dominated by argillaceous sediment with single lithology, has large reflection coefficient at the interfaces with point bar and levee sandstone, and appear as trough reflection on the positive seismic section (Fig. 6a). The natural levee deposit is mainly composed of interbedded thin sand beds. The destructive interference is not obvious because of its thin thickness, the natural levee still appears as one event in seismic section, has seismic amplitude and frequency slightly stronger than the floodplain, and often shows weak amplitude on the positive seismic section (Fig. 6b). The point bar is resulted from lateral erosion of riverbed and lateral accretion of sediment, which is characterized by bedding scale decreasing and positive rhythm of grain size from coarse to fine from bottom to top. Its lithologic interface is the abrupt interface from sandstone to mudstone which has large wave impedance difference. The point bar is characterized by high frequency medium strong-strong amplitude reflection on the positive seismic section, with reflection intensity much higher than natural levee (Fig. 6c).

Based on single well sedimentary microfacies division, same types of sedimentary microfacies around the wells were picked according to different seismic facies models of different sedimentary microfacies to increase sample points. The seismic attribute time window selection and eigenvalues in the vertical range of the sample points were specified to get the intelligent recognition sample set of meandering river sedimentary microfacies. The Random Forest Method^[15] was used to train and predict the sample set to realize the intelligent recognition of sedimentary microfacies of sand body.

Random Forest Algorithm is to build multiple decision trees through random repeated sampling technology and node random splitting technology based on one decision tree, and finally combine the prediction results of a large number of decision trees and output them as a whole. Integrated learning by multiple decision trees can overcome effectively the problems of overfitting and low classification accuracy resulted from one decision tree, and reduce effectively the generalization error of the learning system^[16].

2.3.1. Picking sample points

Based on the interpretation of sedimentary microfacies at well sites, sedimentary microfacies interpreted from all intervals at the wells were extrapolated to the surrounding strata near the wells according to the plane development law of sedimentary microfacies of fluvial facies and the vertical seismic response characteristics of different sedimentary microfacies. New sedimentary microfacies samples were picked up again from volumes of the seismic amplitude and waveforms by slicing analysis to increase samples for machine leaning. 43426 samples of sedimentary microfacies were extracted in time domain, including 16635 point bar samples, 7491 natural levee samples, 19338 floodplain samples.

2.3.2. Seismic attribute specification

The original seismic data were divided into five frequency volumes, i.e. 5-100 Hz, 5-50 Hz, 50-100 Hz, 5-25 Hz and 25-50 Hz. The five frequency bands represent the related waveform information of a typical frequency band. Some high-frequency information would be covered up in low frequency section, while low-frequency information would be covered up in high frequency section. The frequency bands have little overlapping with each other, and no obvious Gibbs phenomenon in time domain, meeting both resolution requirements of time domain and frequency domain, so they can accurately characterize the response characteristics of seismic wave in the corresponding frequency bands. Waveform, instantaneous amplitude and wave impedance attribute were extracted respectively from the original seismic record and five volumes of frequency division bands to obtain 18 volumes including different seismic attributes as subsequent feature input data.

The different properties of different sedimentary microfacies strata not only reflect in the characteristics of the seismic attributes of their current positions, but also the seismic responses of the strata above and below their current positions. It is necessary to specify the time windows and eigenvalues of seismic attributes vertically, so as to fully characterize the seismic response characteristics of different sedimentary microfacies. In this study, one seismic attribute was specified in the time window with fixed length (w) in three intervals with five eigenvalues. For example the eigenvalues of seismic attribute (X) at time t include attribute mean (X_C) at the current interval (range of time window from t-0.5w to t + 0.5w); attribute mean (X_U) above the current interval (range of time window from t-1.5w to t-0.5w); attribute mean (X_L) below the current interval (range of time window from t + 0.5w to t + 1.5w); the mean gradient of the attribute in the three intervals (X_D); the gradient ratio of the upper and lower intervals (X_R). The seismic attribute X in time window t-1.5w to t + 1.5w has five eigenvalues X_C, X_U, X_L, X_D and X_R, which can effectively characterize the variation trend of the seismic attribute in this range, and the possible morphological characteristics such as increasing, decreasing, convex and concave, and have a good corresponding relationship with the seismic facies response of different sedimentary microfacies. For example, the typical wave group characteristics of point bar microfacies combination features negative polarity at the current interval and positive polarity at the top and bottom, and negative wave magnitude of the bottom is significantly higher than that of the top. For the original waveform W and the instantaneous amplitude A, 10 eigenvalues were obtained according to the aforementioned specification method, and their relationships are W_L>W_U>0>W_C, A_U<A_L, W_R<1, A_R<1, and W_D>0 etc.

For the seismic facies characterization of a certain type of sedimentary facies, it is very important to determine the appropriate time window specification. In this study, 2, 4, 6, 8 and 10 ms were selected respectively as the specified window lengths for the original seismic record sampling at 2 ms intervals, and the classification prediction performance of each sedimentary microfacies within different specification lengths were estimated by random forest algorithm^{[15,16,17,18,19]}. It is found the optimal specification time windows for point bar, natural levee and floodplain are 2, 4 and 6 ms, respectively.

2.3.3. Construction of sample set for sedimentary microfacies intelligent recognition

Based on the 43426 sedimentary microfacies sample intervals mentioned above, the 18 seismic attributes extracted from the three sedimentary microfacies of meandering river, i.e. point bar, natural levee and floodplain, were specified at time window of 2, 4 and 6 ms, respectively to obtain their respective training sample sets. Five features are extracted from each seismic attribute, together with three coordinates of inline, crossline and time representing the space-time position, the whole sample set contains 93 dimensions.

2.3.4. Intelligent prediction of sedimentary microfacies sand body

Because only three types of sedimentary microfacies of meandering river, i.e. point bar, floodplain and natural levee, were involved in the machine learning training. But in fact the microfacies of meandering river sediment are far more than these. Deterministic method would result in inconsistency with the actual situation. Hence, in this study, the fuzzy classification intelligent prediction method was used for the prediction to calculate the probability values of the three types of sedimentary microfacies, point bar, floodplain and natural levee, in the target interval respectively. If the probability value of one sedimentary microfacies is particularly large while the probability values of other sedimentary microfacies are small, for example the probability value of point bar is 0.88, that of natural levee is 0.10 and that of floodplain is 0.02, it can be considered that the target is indeed the sedimentary microfacies with the superior probability. If the probability values of two types of sedimentary microfacies are both larger, for example, the probability values of point bar, natural levee and floodplain are respectively 0.50, 0.45 and 0.05, the target is probable at the boundary of the two superior sedimentary microfacies. Due to the unavoidable errors resulted from the original seismic survey, processing and calculation, the recognition results at this point have some uncertainty. If the probability values of the three types of sedimentary microfacies are similarly low, for example, the probability values of point bar, natural levee and floodplain are respectively 0.33, 0.33 and 0.34, the target may be sedimentary microfacies other than the three. In the specific prediction process, probability prediction models of sedimentary microfacies for point bar, natural levee and floodplain, were established by different specification window lengths. For a given data point to be predicted, the probability of each sedimentary microfacies was predicted through the three models, and then the intelligent discrimination was done by above logic.

Through intelligent prediction, the 3D probability prediction volumes of the point bar, natural levee and floodplain in Layer N₁g₁2⁵ meandering river facies in Chengdao extra-shallow sea oilfield, were obtained. Then, the 3D prediction volumes were sliced for the three types of sedimentary microfacies prediction constrained with Layer N₁g₁2⁵. The sedimentary microfacies with smaller probability were set as transparent in the color scale. Then, the prediction probability slices of three types of sedimentary microfacies were displayed by superposition to obtain finally the intelligent prediction results of Layer N₁g₁2⁵ sandbodies (Fig. 7). The predicted point bar has clear edge, shape in accord with the natural shape of channel, and the predicted area of point bar is basically at the edge of river channel, and tally in the plane contact relationship with the actual fluvial facies sedimentary law. The top and bottom of sandbody in Layer N₁g₁2⁵ were traced and interpreted under the constraint of Layer N₁g₁2⁵ top by using point bar probability prediction volume, realizing dual horizon-facies constraint intelligent recognition and delineation of meandering river sandbodies in Chengdao extra-shallow sea oilfield (Fig. 8). The statistical results of sand body thickness from actual drilling and prediction (Table 1) show that the prediction error is less than 1.5 m, marking significant improvement in the accuracy of sand body delineation of meandering river thin reservoir.

In view of the low accuracy of fluid prediction by poststack seismic data in Chengdao extra-shallow sea oilfield in the early stage, fluid in channel sand reservoir in the Upper Guantao Member was identified by prestack seismic data through petrophysical analysis under the condition of ensuring the good quality of prestack gather data.

Fluid elastic impedance is a function of P-wave velocity, S-wave velocity, density and incident angle. In order to connect fluid elastic impedance with seismic data, seismic data volume must be in the form of partial angle stacking. Only common middle point (CMP) gathers reflecting the relationship between amplitude and offset can be obtained after conventional seismic acquisition and processing, so angle gather stacking is a key link in prestack inversion^[20,21]. The purpose of partial stack processing of gather is to provide seismic data for prestack amplitude variation with offset (AVO) or elastic impedance inversion, so CMP gathers need to be specially processed: (1) Fine wave front diffusion processing; (2) correction of source array and geophone array effect; (3) inverse Q filtering; (4) surface consistent processing, including surface consistent deconvolution, surface consistent amplitude correction and surface consistent static correction; (5) prestack denoising; (6) prestack residual amplitude compensation; (7) fine first break cutting. These processes affect directly the AVO (or amplitude variation with incident angle (AVA)) attribute of seismic gathers.

CMP gathers processed were transformed further into angle gathers. According to the principle that the maximum angle cannot exceed the maximum offset and ensure the highest illumination of the target interval, partial stack volumes in three angles of the study area were obtained, which provided base data for prestack fluid detection.

Selecting appropriate parameters can reduce the uncertainty and multiple solutions of fluid detection. Rock elastic parameters interpreted from logging data (P-wave velocity, density and S-wave velocity) were used to calculate attribute parameters such as Poisson's ratio, Lame parameter×density, shear modulus×density and fluid factor. Combined with the estimated logging curves with shear wave curve, the properties of fluids were determined according to cross plotting of multiple parameters^[22,23,24]. Fig. 9 is the cross plot of elastic parameters of Layer N₁g₁2⁵ in Chengdao extra-shallow sea oilfield, it can be seen that oil layers have the characteristics of low Poisson's ratio, low fluid factor, and low product of Lame parameter and density, so fluid can be identified to some extent. But the oil, water and dry layers have some overlap, it is impossible to use a simple cut-off value or polygonal boundary to interpret oil-bearing property, so it is necessary to carry out fluid probability analysis^[25,26].

Fluids at wells were counted and simulated stochastically to establish the probability density distribution functions of elastic properties of different fluids. And then the probability mapping transformation of the product of Lame parameter by density and Poisson's ratio was carried out in oil and water by using Bayesian discriminant criterion to obtain the corresponding oil-bearing probability volume and water bearing probability volume. The method makes full use of the elastic parameters obtained by inversion, the fluid data resulted has more clear geological meaning compared with data resulted from direct inversion, and the fluid prediction risk can be evaluated semi-quantitatively according to the probability distribution trend of the fluid. The technical process of fluid probability analysis mainly includes the following steps: (1) Relationships between seismic inversion elastic parameters under conditions of various fluid-bearing conditions are established through the analysis of logging data. (2) The elastic parameters, which are sensitive to fluid facies and weak linear correlation, are selected based on logging analysis results. The probability density functions of various fluids are established and adjusted according to the logging analysis results. (3) Prestack seismic data are used for prestack synchronous inversion to obtain fluid elastic attribute volumes. (4) The defined probability density function is used to combine logging data with seismic inversion results to analyze the probability distribution and semi-quantitatively predict the type of fluid. Through the above steps, the multiple attribute data volumes from prestack inversion are converted into fluid probability volumes to further improve the accuracy of fluid prediction.

The final fluid probability results integrate responses of multiple prestack elastic parameters to the oil layers, and the prediction results are highly consistent with the actual drilling results. The profile of fluid probability crossing Wells CB252, SHHG2 and CB206 predicted by seismic data shows very clear characteristics of oil-water relationship (Fig. 10). Fig. 11 is the predicted fluid probability of N₁g₁2⁵ sand on the plane, which shows clear distribution of oil and water. Comparison with drilling data shows the coincidence rate of reservoir fluid identification is more than 90% (Table 2).

Three techniques, sea and land dual-sensor prestack two-step high resolution frequency bandwidth expanding processing, intelligent horizon-facies controlled recognition and prestack multi-parameter fluid probability prediction, have been used together to provide geophysical support for exploration and development of thin reservoirs in the Shengli Chengdao extra-shallow sea oilfield.

The sea and land dual-sensor prestack two-step high resolution frequency bandwidth expanding processing consists of sea and land dual-sensor consistency fusion and prestack large angle gather frequency bandwidth expanding while preserving amplitude. It can improve S/N ratio and resolution of seismic data. The dominant frequency of target layer can often be increased from 30 Hz to 50 Hz, as a result, 6-m-thick sandbodies can be identified. The seismic data after processing provides data base for thin reservoir prediction.

The intelligent horizon-facies controlled recognition technique includes establishment of seismic facies models of three types of sedimentary microfacies of meandering river, floodplain, natural levee and point bar on the basis of seismic fine horizon constraint, setup of seismic attribute specification, and intelligent identification and delineation of seismic sedimentary microfacies and thin reservoirs through massive data analysis and intelligent recognition based on the orderliness and invariability of sedimentary facies distribution and the regularity and correlation of seismic attributes. This technology can greatly improve the accuracy and efficiency of thin reservoir delineation, with the prediction error of reservoir thickness of less than 1.5 m.

On the basis of prestack gather processing and petrophysical analysis, the fluid probability semi-quantitative identification technique by multiple prestack seismic parameters has a coincidence rate of identified reservoir fluid with drilled fluid of more than 90%.