An isochronous stratigraphic framework is the premise of seismic sedimentology analysis, and high-quality seismic horizon interpretation is the key to establish an isochronous stratigraphic framework. From the perspective of seismic sedimentology, seismic events are not strictly isochronous, and their geological significance is related to frequency. Generally speaking, seismic events with low frequencies reflect lithological information, and seismic events with high frequencies reflect isochronous sedimentary information ^[21-22]. The events which do not change significantly and maintain stable as the seismic frequency increases, are tend to be more isochronous. Therefore, it is a good way to find isochronous seismic events by analyzing the variation of seismic events with the change of frequency. In practice, we extracted the dip attributes of the original data and the seismic data after frequency division processing, respectively, and calculated the difference of two dip attributes to judge whether the seismic events change with frequency. More specifically, the larger of the dip difference, the more likely the seismic events are isochronous. In addition, appropriate frequency division within the effective frequency range can improve the resolution of seismic data, especially for thin layers. Choosing a suitable frequency division method for a specific geological object is important. In the research of seismic structural interpretation, the configuration of seismic events is the key point. Sedimentary reservoir prediction emphasizes amplitude maintenance, and hydrocarbon prediction focuses on the maintenance of the relative relationship of frequency bands. Generalized spectral decomposition ^[23-24] is a common post-stack frequency division technique that can improve the accuracy of seismic structural interpretation.

Lithology changes rapidly in the fluvial-deltaic deposits, and the continuity of the seismic events is poor. The dip difference was used to assist seismic horizon interpretation. At the same time, the frequency-division seismic data was referred to determine the interpretation scheme and avoid the diachronism in seismic events. Taking the H6 layer as an example (Fig. 3a), the resolution of the original seismic data in the study area is low, and the main frequency is 28 Hz (Fig. 3b). On the basis of accurate well-seismic calibration, the wave peak was traced from Well B7-1W-1 to Well A-3. While wave peak bifurcates up and down in the middle of the seismic profile to the east, it is easy to trace downward on the original seismic profile, which is the H7 layer in Well A-3, resulting in diachronism. By overlaying the conventional seismic profile with the frequency-division dip difference, it can be clearly seen that the dip difference in the middle of the seismic profile is relatively large (Fig. 3a), indicating that the event is poorly diachronic. Therefore, it is necessary to refer to the frequency-division seismic data to improve the accuracy of seismic interpretation. The seismic resolution is higher on the seismic profile at 35 Hz dominant frequency, and the peak event in the middle of the seismic profile becomes narrower (Fig. 3c, 3d). Combined with the geological stratification in Well A-3, it is more reasonable to interpret the H6 layer upward from the bifurcation point of the wave peak.

Rock physical analysis is the premise of seismic lithology prediction. −90° phase shifted seismic data can represent the lithology when the impedance contrast is obvious between sandstone and mudstone, and the reservoir thickness is moderate (less than 1/2 seismic wave-length). Taking low-impedance sandstone, negative polarity and zero-phase seismic data as an example, the top of the thin sandstone is generally calibrated at the maximum of the trough through accurate well-seismic calibration. After −90° phase shifting, the top of the thin sandstone is at zero-phase, and the sandstone is completely corresponding to the seismic trough events. In general, it has been proved to be the most convenient way to characterize thin layers by shifting the seismic phase by −90° ^[11-12]. However, it is not effective for predicting thick sandstone layers or intervals with small impedance contrast between sandstone and mudstone. A thick sandstone interval has multiple sets of trough reflections on the −90° phase shifted seismic profile when the reservoir interval is thicker than half the wavelength. In addition, when the impedance contrast between sandstone and mudstone is not obvious, conventional seismic data cannot identify sandstone and mudstone. It is necessary to find sensitive parameters through petrophysical analysis, and conduct seismic inversion to represent the lithology. However, inversion also has uncertainties in the prediction of thin layers. Therefore, it is necessary to select appropriate seismic data to establish regional seismic lithologic bodies, especially when the petrophysical characteristics change rapidly in different areas. Multi-attribute splicing is a must to get accurate and reliable seismic facies in large areas. In theory, the splicing boundary of attributes is consistent with the point where the petrophysical characteristics change, but it is difficult to determine the specific point on plane. Buried depth has a great influence on rock physics. With the increase of depth, the lithology would transition from sandstone with low impedance to sandstone and mudstone with overlapped impedance, and finally to sandstone with high impedance. Therefore, the splicing boundaries can be determined by the assistance of structural depth. In practice, splicing attributes can be carried out by referring to the continuity of the geological morphology of plane attributes. Conventional seismic interpretation describes “morphology”, while seismic inversion delineates “lithology”Take H3 layer as an example. The shallow sandstone in the west of the study area (Well D11) is of low impedance (Fig. 4a) while the relatively deep sandstone in the east (Well B1-1-2) is of overlapped impedance with small difference between sandstone and mudstone (Fig. 4b). Conventional −90° phase shifted seismic attributes (Fig. 5a) can better characterize the low-impedance sandstone in the west of the H3 layer, and the geological morphology of the channels is clear. But conventional seismic amplitude slices are not effective for distinguishing the sandstone from the mudstone in the east, and the geological morphology is not typical. Compared with post-stack conventional seismic data, pre-stack gradient seismic data is better for characterizing overlapped impedance sandstone. The eastern meandering rivers are typical in the minimum amplitude slices of gradient seismic data (Fig. 5b). The corresponding plane attributes with clear channel morphology and high good fitting degree were spliced according to the continuity of the shape of the meandering river channel (Fig 5c). The spliced seismic attributes of H3a sublayer (H3 is divided into H3a, H3b, and H3c) have the advantages of both conventional seismic slices and gradient slices, laying a solid foundation for regional sedimentary interpretation. In addition, for the formation with special lithology, appropriate inversion attributes should be selected. Taking H7b sublayer (H7 is divided into H7a and H7b) as an example, the reservoir in Well A-3 in the south of the study area has thin coal seams reflected as strong trough amplitude on the seismic profile, and it is difficult to distinguish from sandstone (Fig. 5d). Pre-stack synchronous inversion (Fig. 5e) can better characterize the reservoir in the south of H7b sublayer. The whole zone of the H7b attribute can be obtained by splicing the conventional seismic attribute with the sand content from pre-stack synchronous inversion (Fig. 5f).

Slicing is the core technology for seismic geomorphology characterization in seismic sedimentology. Compared with the traditional seismic interpretation of sand bodies, slicing has certain advantages. Affected by limited seismic resolution, not all sand bodies can be traced or interpreted. Slicing adopts the interpolation of the top and bottom of the layer, and extracts the attributes within a small time window. Combined with the sand body development in well sites, we can analyze the sedimentary evolution by slicing attributes. Currently, there are three slicing methods: time slicing, horizon slicing and stratal slicing. Stratal slicing is the method most used for the formation that is neither sheetlike nor flat-lying. But stratal slicing doesn't consider the change of depositional rate with time and location ^[25], because it is linear. So it is unsuitable for characterizing some special sedimentary bodies that the depositional rate varies with time and location. For example, the depositional rate of the main body is higher than the edge of the channel. In addition, it is difficult to get isochronous stratal slices in some situations, such as overlapping and downlap at top and bottom.

Non-linear strata slicing is an extension of conventional stratal slicing, and it can solve the problem of diachronism. Non-linear strata slicing not only consider the change of the depositional rate with the isochronous reference plane, but also consider the impact of the change of the termination point of the isochronous reference plane. At present, there are two types of nonlinear slicing, namely nonlinear chronological stratal slicing based on Wheeler transform ^[26] and locally optimized nonlinear stratal slicing ^[25]. The former can form relative geologic chronology intuitively, and without stretching or compressing. It has higher requirements on isochronous seismic unit, and is difficult to implement in cases of complex stratigraphic contact relations and poor data quality. The latter extracts the amplitude of adjacent slices for the delineated area that does not match the depositional model by adopting a bell-shaped function or elliptic function, and form a series of nonlinear stratal slices, which can quickly restore the discontinuous channels and missing fan bodies. The latter method was used in this study. Taking H3b sublayer as a case, the conventional stratal slices showed NS-oriented channel deposits. Due to the incision of the channel and the high depositional rate, the attributes of the conventional stratal slices only represented the mudstone deposits at the top of the channel in some area, which made local channel deposits lost and discontinuous (Fig. 6a). By using local nonlinear slices, the channels were re-imaged, and continuous NS-oriented distributary channels were obtained (Fig. 6b).

Attribute fusion generally includes RGB frequency-division fusion and neural network attributes fusion ^[10]. RGB frequency fusion is mostly used to characterize the plane distribution of reservoirs with different thickness, and it can improve the geological morphology of attributes. With the development of machine learning, various neural network algorithms are applied to the optimization of attributes, such as error backpropagation neural network (BP), probabilistic neural network (PNN), convolutional neural network (CNN) and support vector machine (SVM). Neural network attributes fusion has two aspects of contribution: (1) It can greatly improve the geological morphology of the attribute, highlight the boundary, and assist the interpretation of sedimentary microfacies. For example, on the basis of attribute optimization, the unsupervised clustering of BP neural network can be used to quickly and automatically divide sedimentary microfacies, and reduce the uncertainties of single attribute. (2) It can quantitatively predict the plane thickness of sand bodies. For example, convolutional neural network can improve the prediction accuracy of the thickness of the sand body with many faults and a few samples taken from wells. The original attributes input for fusion should well correlate with the thickness of the sand body. Multi-attribute fusion has advantages over single attribute in sedimentary interpretation, but only act as auxiliary tool.

As the background for the development of sedimentary strata, paleogeomorphology not only reproduces the original tectonic framework, but also constructs the ancient source system ^[27]. At present, many scholars adopt the impression method to approximately reflect paleogeomorphology ^[28]. Paleogeomorphology reconstruction does not simply restore the paleo-thickness of strata, but also consider geological processes, such as the tectonic evolution, missing strata, sediment gravity subsidence, and compaction. The restoration process should include paleo-thickness restoration of sequence strata, missing strata restoration, compaction correction, crust isostatic restoration and fault restoration. This paper takes measures such as compaction correction, denudation recovery and differential subsidence correction to restore the paleogeomorphology of the target layer on the basis of the residual thickness (Fig. 7).

The procedures of paleogeomorphology restoration are as follows: (1) Subtract the structural depth between the reference planes to obtain the residual thickness of the target layer (Fig. 7a). (2) Considering the formation dip (Fig. 7b), calculate the true formation thickness after dip correction (Fig. 7c). (3) Conduct compaction correction to restore the sedimentary thickness before buried and remove the impact of compaction. By analyzing the relationship between the porosity and the depth of different lithofacies, calculate the compaction coefficient and initial porosity. Analyze the burial history and compaction ratio by using basin simulation software. The planar distribution of compaction coefficient (Fig. 7d) is obtained by interpolation, which is used for the restoration of paleogeomorphology after compaction correction (Fig. 7e). (4) There exists local denudation in the west of H5 layer, which should be recovered. The denuded thickness in a single well can be quantitatively restored according to the acoustic transit time of mudstone. According to the structural trend, extrapolate the denudation volume in a single well to the entire H5 layer (Fig. 7f). (5) Calculate the structural subsidence at initial deposition period (Fig. 7g) and restore the paleowater depth (Fig. 7h) through gravity balance and correction. After five steps above, the final paleogeomorphology map is obtained (Fig. 7i). The final paleogeomorphology (Fig. 7i) has the same trend as the residual thickness (Fig. 7c), but there are some local differences in details. The uplifted scope and absolute elevation of the slope are more reasonable.

Paleogeomorphology can assist the sedimentary interpretation based on plane stratal slice attributes. Using multiple methods can enhance the reliability of the results. Taking the paleogeomorphology of the H5 layer (Fig. 7i) as a case, the west is higher than the east. From west to east, H5 layer can be divided into a denudation zone (pass-by zone), a slope zone and a low-lying zone. The paleogeomorphology well matches the distribution of sand bodies. Sand bodies are enriched in the low-lying zone while less developed in the slope zone with limited accommodation space.

H4b (H4 is divided into H4a, H4b, H4c, and H4d) is a thin interbed. Frequency-division seismic data was used for horizon tracking and interpretation to improve the isochronism of the stratigraphic framework. Attribute analysis found the minimum amplitude, average energy and arc length had the highest contribution to BP neural network fusion. The result shows a clear delta shape of H4b attribute (Fig. 8a). On the plane, it is a narrow belt and potato-like shape developed in a shallow water delta environment. There are NW-SE distributary channels (Fig. 8b) whose main bodies are located at the west of the study area. The sandstone-to-formation ratio is 38%, and the thickness of single sand body is small as 2.8 m. The logging curves are in the shapes of bell and superimposed boxes. On the whole, the sand body is thicker in the south and thinner in the north. Wells A-A7 and A-A8 were drilled into the main bodies of the channel in H4b sublayer, where the sand bodies are 20 m thick (10 m per sand body). Seismic prograding reflections of delta can be seen along the source direction. The river channel on the delta plain shows continuous trough reflections (Fig. 8c). The seismic section perpendicular to the delta front channels shows point-like lenticular reflections which are seismic responses of distributary channels (Fig. 8d).

For H5 with overlapped impedance of sandstone and mudstone, the minimum amplitude of pre-stack gradient seismic data (Fig. 9a) is better for lithology characterization while the sedimentary microfacies can be analyzed by clustering attributes (Fig. 9b). The channels with the clear minimum amplitude of H5 layer are north-south braided rivers (Fig. 9c), and the sand bodies are distributed in a large scale in the study area. H5 layer was developed in the early-middle stage of the mid-term base-level ascending cycle, with sand bodies superimposed vertically. The sand bodies have been drilled in 30-50 m thick, with high sand-to-formation ratio more than 50%. The braided river deposits are mainly braided channels and channel bars. The logging curve of the braided river channels is box- or bell-shaped, with a sudden change at the bottom, indicating sufficient source supply and strong hydrodynamic energy. The down-cutting effect is obvious, and the composite sand bodies with large scale and good lateral continuity are formed. The width of the braided river is about 2-8 km, showing strong trough reflection axis on the seismic section. Lenticular sand bodies can be seen in the direction of cross-cutting the river channel (Fig. 9d). The channel bar is the product of vertical accretion of multiple sand bodies. The sand bodies are distributed in a large scale, and the GR curve is box-shaped (Well B1-1-2). The amplitude of the channel bar is stronger than the braided channel. The channel bar can be distinguished by attributes clustering (Fig. 9b).

Due to the development of coal seams, the lithology in H7b sublayer (H7 is divided into H7a and H7b) can no longer be effectively distinguished using a single conventional seismic attribute (Fig. 5d). But pre-stack synchronous inversion can better characterize the distribution of the sand bodies in the coal-bearing background (Fig. 5e). By splicing the minimum conventional seismic amplitude with the minimum amplitude of the southern pre-stack synchronous inversion, we got the final planar attribute of H7b. The geological morphology of H7b sublayer is clear (Fig. 10a), and shows a large north-south meandering river (Fig. 10b). Wells B7-3-1 and B2-1-1 drilled in the point bars where the sandstone is 40 m thick. Well A-B7 drilled in the edge of a point bar. The rest wells drilled in the mudstone in the flooded plain. The point bar complex has a span of about 10 km and shows multiple lenticular troughs on the seismic profile (Fig. 10c). Referring to the weak amplitude between the point bar complexes and the thickness (20 m) of a single point bar in Well B7-3-1, the point bar complex in the middle can be divided into three stages, spanning about 6 km per stage. In the south (Well A-3), thin coal seams show high frequency and strong trough reflection on the seismic profile (Fig. 10d).

The study area is a fluvial-deltaic sedimentary system where the sedimentary evolution of the Huagang Formation (Fig. 11a) is as follows from the bottom up: H12-H11 are braided rivers, H10-H9 are shallow water deltas, H8 is a braided river, H7 is a meandering river, H6 is a shallow water delta, H5 is a braided river, H4 is a shallow water delta, H3 is a braided river, and H2-H1 are meandering rivers.

The north-south axial source is dominant in the study area. The development of sand bodies is controlled by the mid-term base-level cycle and paleogeomorphology. With the rise and fall of the lake level, fast lake transgression and regression events occurred many times during the deposition of the Huagang Formation in the study area, resulting in a series of vertically successive fluvial-shallow water delta systems. (1) From the early to middle ascending stages of the mid-term base-level cycle (Fig. 11b), the lake level was low, the accumulation space was large and source supply was sufficient, and braided channel deposits developed mainly (H5 layer). The channel was straight, the sand bodies were continuous, and superimposed vertically. The fluctuating and steep paleogeomorphology (Fig. 11c) and restricted valleys imposed control on the distribution of the sand bodies. Most sand bodies are concentrated on the low-lying valleys, but not developed on the slope. (2) From the late ascending stage to the early descending stage of the mid-term base-level cycle (Fig. 11d), the lake level rose rapidly to a high level, the accommodating space was moderate, and the source supply became less, resulting in meandering rivers (i.e., H7b sublayer). The sand bodies are distributed in narrow strips, and isolated in the vertical direction. The paleogeomorphology (Fig. 11e) is relatively flat, with unrestricted valleys, and the meandering rivers migrate quickly laterally. (3) From the middle to late descending stages of the mid-term base-level cycle (Fig. 11f), the lake level fell, entering the lake regression. The accommodating space became smaller, and the source supply was moderate, resulting in shallow water deltas, such as H4b. The sand bodies in the delta plain are distributed in locally narrow, continuous and interwoven strips. Vertically, they are migratory and in medium-large thickness. They are lobe-like underwater distributary channels and mouth bars at the mouth of the delta front. The paleogeomorphology is relatively flat (Fig. 11g), and the delta front slope is less than 0.1°. The deposit is mainly shallow water delta. In the south of the study area, fore-delta-shallow lacustrine deposits can be observed. The lake shoreline changes frequently.

On the basis of reservoir characterization, structural features and fault development, favorable targets were selected in the study area. On a favorable structural background, favorable structural-lithological traps may develop in the case of rational spatial arrangement of trap-forming elements such as sand bodies and faults. Sealing conditions in high parts decide the effectiveness of lithologic traps. Through the analysis of the sealing condition in the study area, there are mainly two types of structural-lithological traps, namely lateral reservoir pinch-out and sealing fault structural-lithological traps (Fig. 12a) and reservoir updipped pinch-out structural-lithologic traps (Fig. 12b).

Under the guidance of the types of trap models, two potential targets were selected. The delta plain channel in H4b sublayer pinches out laterally, and the high part is blocked by faults, forming a good lateral pinch-out and sealing fault structural-lithologic trap (Fig. 12c). On the typical seismic section across the H4b target layer, the local amplitude is strong (Fig. 12e) and the structural low of the sand body connects the oil-source fault, demonstrating a high possibility of hydrocarbon accumulation. The meandering river point bar in H7b sublayer pinches out toward the structural high, forming an updipped pinch-out structural-lithologic trap (Fig. 12d). On the seismic section across the H7b target layer, the phase of the seismic events shifted to the structural high, the pinch-out point is clear (Fig. 12f). The amplitude of the structural high is stronger than that of the structural low, which may be caused by gas. There are many similar updipped pinch-out structural-lithologic traps in the structural wings of the Xihu sag, which have a broad prospect for hydrocarbon exploration and development in the future.