Various types of fluvial facies such as straight, meandering and braided fluvial facies have become the research hotspot in the field of oil and gas resources for many years. Many large- and medium-sized oil and gas fields of fluvial facies have been found in West Siberia of Russia, Wyoming of the United States, Ordos Basin and Sichuan Basin of China. As oil and gas exploration and development is deepened in China, the target for fluvial resources has gradually shifted from large- and medium-sized to micro-sized oil and gas reservoirs. Microscopic ancient river channels are widely distributed in various types of formations, such as shallow, medium and deep transition facies between land and sea. They have rich oil and gas reserves and huge exploration potential, but they have disadvantages of narrow channel width, thin sand bodies and complicated oil and gas accumulation^[1,2]. In addition, with the increase of burial depth, seismic reflection exhibits unfavorable characteristics such as weak energy, dominant low frequency, narrow bandwidth and low vertical and lateral resolutions, which restrict the accuracy of narrow ancient river channels, thin reservoirs and other oil and gas information extract. The identification of microscopic ancient river channels has become a prominent problem in increasing well sites, horizontal well drilling and other exploration and development risks.

Aiming at the identification problem of microscopic ancient river channels, generally we start from the perspectives of seismic acquisition, processing and interpretation, and use high-resolution seismic data to obtain key parameters and spatial distribution characteristics such as the lateral width and vertical thickness of ancient river channels. As far as high-resolution seismic acquisition is concerned, seismic data collected by advanced recording systems such as high density, "2W and 1H" (wide azimuth, wide frequency band, and high density) can improve the accuracy of prediction and description of shallow and middle reservoirs. However, for deep geological targets, it has so far failed to break through bottlenecks in excitation and reception, such as the rapid energy attenuation of the seismic source and the difficulty of receiving high-frequency effective signals by the geophone. In terms of high-resolution seismic data processing and interpretation, methods like inverse Q compensation^[3], spectrum shaping^[4], bandwidth expansion^[5], thin-layer reflection coefficient inversion^[6], low-frequency extension and deep learning^[7], high-resolution impedance inversion^[8,9] can make more reflection events, wave impedance interfaces, etc. in seismic data. However, they may reduce the SNR (signal-to-noise ratio) and make it difficult to evaluate the wave group effectiveness and with false reflectors. As a result, the reliability of geological prediction is reduced. In addition, the deeper the burial and the smaller the scale, the more difficult it is to identify microscopic ancient river channels, especially those extremely narrow and with complex flow directions. It is necessary to make breakthroughs to seismic resolution limits from the perspectives of seismic energy propagation and reflection frequency, and improve the identification of microscopic ancient river channels.

The key to identifying microscopic ancient river channels is to extract the high-resolution time-frequency responses of seismic signals using the high-resolution time-frequency analysis method, and then use the energy, frequency, phase and other information of seismic waves to characterize the width and flow direction of microscopic ancient river channels effectively. At present, many high-resolution processing and interpretation methods, such as inverse Q compensation and bandwidth expansion, are based on time-frequency analysis, but they have not achieved high-resolution time-frequency response feature extraction, making it difficult to accurately identify small geological targets like microscopic ancient river channels. In fact, since Gabor first proposed the spectrogram of acoustic signals in 1946^[10], three types of time-frequency analysis methods (linear, bilinear and nonlinear) with different kernel functions have been quickly developed^[11]. Among them, Gabor transform, short-time Fourier transform (STFT), continuous wavelet transform (CWT), S transform (ST) and generalized S transform (GST) are linear time-frequency analysis methods. Although they have the advantage of higher computational efficiency, most of them are window-dependent and controlled by the uncertainty principle of Heisenberg W^[12], the resolution ability for time and frequency cannot be improved at the same time. Adaptive STFT^[13], modal decomposition^[14,15], synchronous extrusion wavelet transform^[16], matching pursuit (MP)^[17] are non-linear methods. Although they have the advantage of higher resolution, the amount of calculation is large, the robustness is insufficient, and it is difficult to adapt to the time-frequency analysis of low SNR signals. Wigner-Ville distribution (WVD)^[18], pseudo Wigner-Ville distribution (PWVD), smooth pseudo Wigner-Ville distribution (SPWVD) and Choi-Williams distribution (CWD)^[19] are bilinear methods (also called quadratic). They are not restricted by Heisenberg's uncertainty principle^[12], and their excellent mathematical properties lay a good theoretical foundation for high-resolution time-frequency analysis.

Based on the maximum entropy principle of seismic signals, we studied the methods of enhancing the spectral focusing characteristics of WVD and improving the time-frequency resolution of seismic signals on the premise of overcoming the interference of WVD cross-terms to highlight the response characteristics of narrow and thin ancient river channels. We also used emulational signals, forward numerical simulation and other theoretical data to carry out experiments and demonstration. We developed a new method for microscopic ancient river channels identification, which is helpful to locating new wells and drilling horizontal wells.

Entropy refers to the probability of occurrence of discrete random events or specific information. It was first defined by Shannon C E in 1948, which was mainly used to characterize the uncertainty of random variables or characteristic information^[20]. In 1957, Jaynes E T proposed the maximum entropy criterion, that is, under the premise that the prior information is accurate, the probability distribution that best represents the current information state is the probability distribution with the largest entropy. In other words, when the probability distribution of prior information is uncertain, the information distribution with the largest entropy is reasonable^[21].

In 1967, Burg introduced the maximum entropy criterion into signal spectrum analysis, and Burg algorithm for maximum entropy time-frequency analysis was developed^[22]. Burg believes that entropy is the distribution function of signal $X$, which is only related to the probability of $X$ and has nothing to do with its actual value. When the value of the signal $X$ is x, the distribution density function is $p(x)$, and its entropy is:

Then Burg defines the power spectrum entropy of signal $X$ as:

When using Eq. (2) to obtain the power spectrum, it is necessary to recursively solve the autocorrelation function of signal $X$ at the maximum of $H[P(f)]$ to minimize the power spectrum error and improve the resolution. In fact, Burg algorithm is equivalent to a linear autoregressive AR (Auto Regressive) prediction model, and the calculation formula of the maximum entropy power spectrum^[22] is:

where $0\le f\le \frac{1}{2\Delta t}$

Boashash defined^[23] that when the analytical signal of a continuous time series signal $x(t)$ is $z(t)$, the Wigner-Ville distribution (WVD) is:

However, the measured signal is not continuous. In practical applications, it is necessary to discretize $W(t,f)$. Supposing the discrete expressions of $x(t)$ and $z(t)$ are $x(n)$ and $z(n)$ ($0\le n\le N-1$), respectively, and performing Hilbert transformation on $x(n)$, after calculating $z(n)$, the discrete expression of $W(t,f)$ can be obtained:

where $0\le |l|\le N$

According to the nature of WVD, discrete WVD is the discrete Fourier transform of instantaneous autocorrelation function. The instantaneous autocorrelation function of the discrete signal is:

For each sample point n, ${{k}_{n}}(l)$ will form a convolution kernel sequence ${{k}_{z}}(n)$, as in Eq. (7). When $n=0$ and $n=N-1$, the convolution kernel sequence ${{k}_{z}}(n)={{k}_{z}}(0)$ and ${{k}_{z}}(n)={{k}_{z}}(N-1)$ is the shortest; when $n=\frac{N}{2}$, the convolution kernel sequence ${{k}_{z}}(n)={{k}_{z}}\left( \frac{N}{2} \right)$is the longest. Taking ${{k}_{n}}(0)$ as the center, adding value 0 to the left and right sides of ${{k}_{z}}(n)$, ${{k}_{z}}(n)$ can be expanded to a convolution kernel sequence with length N.

Discrete Fourier transform is performed on the convolution kernel sequence ${{k}_{\text{z}}}(n)$, which is from the n^th sampling point, to obtain the WVD instantaneous power spectrum, as follows:

where the sub-item of ${{W}_{\text{z}}}(n,f)$ is expressed as:

By calculating the instantaneous power spectrum corresponding to each sampling point n, finally the Wigner-Ville power spectrum of discrete signal $x(n)$ can be obtained, namely:

Of course, affected by the bi-linearity of WVD, when using the instantaneous autocorrelation function ${{k}_{n}}(l)$ to multiply the data before and after, cross-term interference will be generated, which is not conducive for high-precision power spectrum feature analysis.

When the first derivative of instantaneous autocorrelation function ${{k}_{n}}(l)$ is 0, the Wigner-Ville distribution calculated is the maximum entropy spectrum. Since the one-dimensional maximum entropy spectrum is equivalent to the AR model spectrum, only the AR model parameters and the autoregression (AR) model spectrum are required to obtain the maximum entropy spectrum^[24]. In this way, it is possible not to use Eq. (10) to calculate the power spectrum, and eliminate the cross-term interference generated when the data before and after the instantaneous autocorrelation function is multiplied.

When calculating the maximum entropy power spectrum $P(f)$ based on Burg algorithm, Eq. (3) and Eq. (10) are equivalent. Based on the Levinson-Durbin recursion rule^[25], after obtaining AR model parameters such as prediction error $E$ and autoregressive coefficients $\alpha $, and substituting them into Eq. (3), the WVD power spectrum constrained by the maximum entropy criterion can be obtained. Burg algorithm recurs the autoregressive coefficient of the AR model through forward prediction error and backward prediction error^[26,27]. Set ${{E}_{\text{F},0}}(n)=$ ${{E}_{\text{B},m-1}}(n)=z(n)z'(n)$ as the initial condition, and the recurrence formula is:

Using discrete signal $x(n)$ and its analytical signal $z(n)$, the initial power of prediction error is calculated:

Using the m^th order forward prediction error and backward prediction error, the average power of prediction errors can be calculated, namely:

According to the Levinson-Durbin recursive constraint rule, let ${{\alpha }_{m}}(j)$ be the j^th autoregressive coefficient of the m^th order AR model when the order is M, and ${{E}_{m}}$ is the minimum prediction error of the m^th order. Then the recurrence relationship is:

In this way, by substituting ${{\alpha }_{m}}(1)$, ${{\alpha }_{m}}(2)$ … ${{\alpha }_{m}}(m-1)$ and ${{E}_{m}}$ into Eq. (3), the maximum entropy power spectrum of WVD can be calculated. The distribution of WVD maximum entropy power spectrum in time domain and frequency domain is the Wigner-Ville maximum entropy time spectrum (referred to as MEWVD).

To sum up, the calculating steps of WVD maximum entropy power spectrum of Burg algorithm^[28] include: (1) Setting $m=1$, calculate ${{P}_{0}}(f)$, ${{E}_{\text{F},0}}(n)$ and ${{E}_{\text{B},0}}(n)$; (2) Calculate ${{\alpha }_{m}}(m)$; (3) Calculate the j^th autoregressive coefficient ${{\alpha }_{m}}(j)$ when the order is m, $j=1,2,\cdots \cdots,$ $m-1$; (4) Calculate ${{P}_{m}}(f)$; (5) Calculate ${{E}_{\text{F},m}}(n)$ and ${{E}_{\text{B},m}}(n)$ when the order is m; (6) Let $m=m+1$, repeat (2)-(5) until ${{P}_{m}}(f)$ is no longer significantly reduced. Using the Levinson Durbin recurrence relationship, AR model parameters can be obtained, and finally calculate the maximum entropy power spectrum of WVD.

In time domain, MEWVD does not multiply the data before and after the instantaneous autocorrelation function. Instead, it only uses the data near the current point n to avoid cross-term interference. In frequency domain, subject to the maximum entropy criterion, the power spectrum corresponding to point n is more accurate and the energy is more focused. At the same time, under the premise of maximum entropy, the Levinson-Durbin constraint is used to recur unknown instantaneous autocorrelation function, so that the calculated power spectrum error is reduced, and the signal time-frequency resolution capability is improved.

It should be pointed out that the maximum entropy criterion is widely used in mathematics, economics, information and other fields. Its major advantage is that it not only retains the uncertainty of information distribution, but also minimizes prediction risk. In oil and gas field, it is difficult to identify micro-geological bodies such as ancient river channels, sand dams, and point reefs that are narrow horizontally and thin vertically because of seismic wave attenuation, harmonic waves and noises. Accurate and highly focused seismic power spectrum is an important basis for identifying micro-geological bodies. MEWVD can best represent the seismic spectrum response characteristics of micro-geological bodies when the entropy is the largest. It can not only achieve the optimization of time-frequency characteristics classification, but also improve the robustness of signal analysis in a noisy environment^[29], which is conducive to improving the reliability of recognition of micro-geological bodies.

When seismic waves propagate underground, they are affected by factors such as formation depth, porosity, filling fluid, etc., which cause energy attenuation, and changes in dynamics, kinematics and geometric characteristics. Using parameters such as amplitude, absorption coefficient, frequency, phase and travel time, the propagating process of seismic waves can be emulated:

In the above formula, the emulational seismic signal ${{x}_{\text{F}}}(t)$ is formed by superimposing the attenuation component ${{x}_{1}}(t)$ and the stationary component ${{x}_{2}}(t)$. Among them, ${{x}_{1}}(t)$ includes parameters such as absorption coefficient, amplitude, frequency, phase and travel time; ${{x}_{2}}(t)$ does not include the absorption coefficient, and the transmitted energy will not be attenuated.

In order to compare and analyze the time-frequency characteristics of seismic waves during propagation process, using parameter ${{x}_{\text{F}}}(t)$, a emulational seismic signal which is strongly correlated with different high- and low-frequency energy is constructed. As shown in Table 1, the initial energy ${{A}_{1}}$ and the propagation time t of signals ${{x}_{\text{F}1}}$ and ${{x}_{\text{F}2}}$ are the same. In ${{x}_{\text{F}1}}$, ${{A}_{2}}$, ${{f}_{2}}$, and ${{\varphi }_{2}}$ are 0, indicating that ${{x}_{F1}}$ is only composed of attenuated components, while ${{A}_{1}}$, ${{f}_{1}}$ and Q are larger, similar to seismic signals with strong energy, high frequency and strong attenuation.${{x}_{\text{F}2}}$ contains ${{x}_{\text{F}1}}$. Only by changing ${{f}_{2}}$, the steady component of low frequency and weak amplitude is added to ${{x}_{\text{F}2}}$, which means that the propagating energy and amplitude of ${{x}_{\text{F}2}}$ within t will change with frequency. The first column in Fig. 1a and Fig. 1b shows the waveform characteristics of ${{x}_{\text{F}1}}$ and ${{x}_{\text{F}2}}$, respectively. It can be seen that the amplitude of ${{x}_{\text{F}1}}$ decreases with t, and the energy decays quickly; as t changes, the amplitude of ${{x}_{\text{F}2}}$ decreases, and the energy is attenuated too. However, the control effect of the attenuation component is getting weaker and weaker. After the signal propagates about 0.25 s, the signal is mainly controlled by the stationary component. And due to the difference in ${{f}_{2}}$, if the oscillation period of ${{x}_{\text{F}1}}$ is short, the waveform is "skinny"; if the oscillation period of ${{x}_{\text{F}2}}$ is long, the waveform becomes "fat".

At the same time, after comparing different time-frequency analysis methods, it is found that MEWVD has high-quality time-frequency focusing characteristics and resolution. As shown in Fig. 1a and Fig. 1b, the second column shows the fast Fourier transform (FFT) spectrum of the emulational seismic signal, revealing that the energy peak of signal ${{x}_{\text{F}1}}$ is concentrated at 80 Hz, and the energy peak of ${{x}_{\text{F}1}}$ is concentrated at 10 Hz and 80 Hz. Columns 3 to 10 respectively show the time-frequency effects of STFT, GST, CWT, WVD, PWVD, SPWVD, CWD and MEWVD. Comparative analysis shows that although the strongest energy of ${{x}_{\text{F}1}}$ and ${{x}_{\text{F}2}}$ are both concentrated on the 80 Hz high-frequency point, different time-frequency analysis methods show differences in the degree of energy focusing and the capability of frequency resolution. Moreover, the focusing effect and resolution of MEWVD are the best, followed by WVD, and STFT is the weakest. GST and CWT are similar in time-frequency resolution. Both of them can effectively identify signals with 80 Hz high frequency and strong amplitude and 10 Hz low frequency and weak amplitude, but the energy focus and time resolution are worse than MEWVD. In addition, WVD time-frequency characteristics also have cross-term interference. Although PWVD suppresses cross-term interference to a certain extent, its energy focusing is weakened, and the resolution ability to identify the low-frequency and weak-amplitude information of ${{x}_{\text{F}2}}$ is reduced. SPWVD and CWD basically suppress cross-term interference, but they can hardly identify the middle- and low-frequency and weak-amplitude response of ${{x}_{\text{F}2}}$, and the resolution loss is more serious than PWVD.

In order to further confirm the spectral focusing and high time-frequency resolution of MEWVD, a two-dimensional numerical model was designed using the parameters shown in Table 2, including the depth, width and thickness of a narrow and thin geological body. Its "lens-like" geological feature is similar to the cross-section of a microscopic ancient river channel. Based on the forward numerical simulation method of the two-dimensional seismic wave equation, the 35 Hz Rick wavelet was used as a point source, and a forward simulation record was obtained at sampling interval of 10 m and sampling rate of 1 ms. Models ④-⑨ gradually become thinner and narrower in thickness of sand bodies. Models ④-⑥ are wider and thicker than 1/4 wavelength (about 30 m), and the reflection events at the top and bottom are obvious. Models ⑦ and ⑨ are narrower and thinner than 1/4 wavelength, and they show harmonic reflections. Models ⑦ and ⑧ show strong harmonic reflections, while the extremely thin and narrow Model ⑨ shows weakly harmonic reflections (Fig. 2).

Separately extracting single-channel seismic signals representing models ④-⑨, and after calculating their MEWVD, the time-frequency difference of the models can be analyzed. Fig. 2 shows the single-channel seismic signals and MEWVD spectrum characteristics of models ④-⑨. It can be seen that MEWVD shows strong spectral focusing and high-resolution characteristics for the reflected signals of the numerical models with different width and thickness. The single-channel seismic signals of the wider and thicker models ④-⑥ show clear reflections from tops and bottoms, and MEWVD shows narrower spectrum distribution in the frequency axis. For the narrower and thinner model ⑦-⑨, their single-channel seismic signals only show the reflection from the top interface, and the bottom interface is difficult to identify due to the effect of wave field tuning. MEWVD shows wider spectrum distribution in the direction of the frequency axis, and the thinner the thickness, the wider the spectrum distribution.

On the single-frequency section, the spectrum focusing of MEWVD is very significant, and the power spectrum at different frequencies exhibit different time-frequency resolutions. Fig. 3 shows the instantaneous amplitude of the forward seismic signal and the MEWVD spectrum characteristics at 25, 35, and 45 Hz. Fig. 3a shows that the instantaneous amplitude is limited by the time-frequency resolution capability. Only the thicker Model ④ shows clear "biaxial" strong instantaneous amplitude at the top and bottom. The gradually thinner models ⑤-⑨ show strong instantaneous amplitude distribution on "single axis". Even the tops and bottoms of models ⑤ and ⑥ greater than 1/4 wavelength are difficult to identify. Fig. 3b to Fig. 3c show the MEWVD spectrum distribution at three frequency points. They all have higher resolution than the instantaneous amplitude. But the MEWVD spectrum responses at different frequency points are different. The 45 Hz MEWVD has the highest spectral resolution, which can accurately identify the top and bottom interfaces and lateral widths of the narrow and thin models ⑤ and ⑥. The 35 Hz MEWVD spectrum resolution is worse, but it can still effectively identify the top and bottom interfaces and lateral widths of ⑤ and ⑥, and has excellent spectral focusing responses to extremely narrow and thin models ⑧ and ⑨. And affected by different widths, the spectral thickness of the two is similar, and the width difference is significant. Further analysis of the extremely narrow and extremely thin models ⑧ and ⑨ shows that the thickness of the two is the same, the width is 5 times different, the lateral difference of their instantaneous amplitude is not obvious, and the MEWVD spectrum anomaly is marked. Especially in Fig. 3d, spectral focusing is excellent at 45 Hz, showing the characteristics of the MEWVD spectrum anomaly of the models with the same thickness and 5 times the width difference. At the same time, although the 25 Hz MEWVD spectrum resolution is generally not as good as the 35 Hz and 45 Hz spectra, there is a significant difference in the spectral focusing response of the extremely narrow and extremely thin models ⑧ and ⑨. The frequency spectra at 35 Hz and 45 Hz are in "clusters", while the spectrum at 25 Hz is "layered". It shows that the maximum entropy spectrum of Wigner-Ville has nothing to do with the strength of the seismic signal. When the tuning effect makes the power spectrum distribution at the top and bottom interfaces of the model uncertain, the power spectrum with the largest entropy is the most reasonable distribution.

In a word, the emulational seismic signals and the narrow and thin forward numerical models show that MEWVD has a good ability of spectral focusing and time-frequency resolution. MEWVD can accurately characterize the time-frequency distribution of signals, and MEWVD at different frequencies can accurately obtain the spectral response on different scales of narrow and thin models. In addition, since the maximum entropy spectrum of Wigner-Ville is only related to the probability of the spectrum distribution, when a geological body is narrow or thinner than 1/4 wavelength with tuning effect, the seismic reflection event cannot be directly identified. The power spectrum with the largest entropy is the most reasonable distribution of the geological body. Therefore, combining the MEWVD spectrum characteristics at different frequencies can identify the spatial distribution of the microscopic ancient river channel similar to the narrow and thin model.

The Sichuan Basin is a multi-cycle sedimentary basin with rich gas resource. As of 2020, in 27 formations including Sinian, Cambrian, Silurian, Carboniferous, Triassic, Jurassic, etc., four types of gas reservoirs have been discovered (structural, structural-lithologic, lithologic- structural and lithologic reservoirs). They are distributed in four major industrial gas regions in the northeast, west, south and north of the Sichuan Basin^[30,31]. Among them, the Jurassic Shaximiao Formation gas reservoir in the western region is mainly distributed in the western Sichuan Depression, which is adjacent to the Longmenshan nappe belt in the west and the central Sichuan uplift in the east. In the east slope of the western Sichuan Basin, a large lithologic-structural gas field buried above 3100 m, Zhongjiang Gas Field, was discovered.

Zhongjiang gas field covers an area of 2350 km². The Jurassic Shaximiao Formation in the field is dominated by shallow delta deposits, with sedimentary microfacies such as distributary channels, distributary bays, estuary dams, and rupture fans in the delta plain. Among them, distributary channels mainly develop sandstone reservoirs such as fine-medium lithic feldspar sandstone, feldspar lithic sandstone, lithic sandstone and lithic quartz sandstone, with average porosity of 8.94% and average permeability of 0.5×10^-3 μm². At present, in the Shaximiao Formation gas reservoir with rich resources, about 150 industrial gas wells have been put into production, with daily gas production exceeding 220×10⁴ m³. However, with deeper development, micro-distributary channel gas reservoirs became more concealed. The reservoirs are dominated by strip-like superimposed thin-layer channel sand bodies, narrow horizontally and thin vertically, making it extremely difficult to finely describe. To support gas reservoir description, horizontal drilling and hydraulic fracturing, it is necessary to develop a precise method for identifying the narrow and thin microscopic ancient river channels.

The Middle Jurassic Shaximiao Formation in Zhongjiang Gas Field is divided into three sections: upper (J₂s₁), middle (J₂s₂), and lower (J₂s₃). Among them, J₂s₃ is composed of J₂s₃¹, J₂s₃² and J₂s₃³. According to the sedimentary characteristics of the ancient channels and sandstones, J₂s₃³ is subdivided into small J₂s₃^3-1 and J₂s₃^3-2. At present, J₂s₃^3-2 is the primary gas-producing interval. The intra-layer networked and banded ancient river channels mainly flow from northeast to southwest, extending 18-58 km long, 0.2-3.5 km wide, and 6.5-45.2 m thick. The reservoirs are distributed in sheet, box and bell shapes, and show seismic reflections at low frequency, with "bright spots" and "blanks". Along the ancient river channels, the seismic reflections are generally continuous and stable, perpendicular to the river channels, and show “lens-like” or weak reflection anomalies.

In recent years, with in-depth development of the J₂s₃^3-2 gas reservoir in Zhongjiang Gas Field, narrow and thin ancient river channel sandstone with width less than 500 m and thickness less than 35 m has become a key target. And accurate identification of microscopic ancient river channels has also become a key link in the development of gas reservoirs. However, affected by many factors such as burial depth, width, thickness, lithology, physical properties, etc., different microscopic ancient river channels show obvious differences in seismic reflection. Fig. 4a shows that the dominant seismic frequency of J₂s₃^3-2 is about 30 Hz, and the frequency band is 8 to 70 Hz; Fig. 4b shows that the ancient channel in J₂s₃^3-2 in Well E has low GR anomalies, with the seismic reflection showing local double-peak "lens-like" anomaly. The channel is predicted to be about 354 m wide and 35 m thick. Fig. 4c shows that there is a low GR value anomaly in ancient channel in J₂s₃^3-2 in Well D. It is predicted that the channel sandstone is about 310 m wide and 30 m thick, which is narrower and thinner than the channel in Well E, and the "lens-like" reflection is weaker. Fig. 4d shows that there is a low GR anomaly in the J₂s₃^3-2 ancient channel in Well C. The channel sandstone is predicted to be about 248 m wide and 30 m thick. Compared with wells D and E, the "lens-like" anomaly disappears, showing a weak reflection response. Obviously, by comparing the seismic reflection characteristics of the microscopic ancient river channels in wells C, D, and E in Fig. 4, as they become narrower and thinner, they are more concealed; and the weaker the seismic reflection, the more difficult to identify.

The MEWVD method was used to extract the spectral characteristics of the Shaximiao Formation in Zhongjiang Gas Field. It’s found that the microscopic ancient river channels with different widths and thicknesses showed different MEWVD anomalies, revealing the spatial distribution of narrow and thin ancient channels. On the instantaneous amplitude profile shown in Fig. 5a, the J₂s₃^3-2 ancient river channel in Well E has a strong instantaneous amplitude anomaly with a large anomaly range, poor focusing feature and low resolution. On the MEWVD spectral profile shown in Fig. 5b-5d, there is a strong anomaly in the ancient channel in J₂s₃^3-2 in Well E. The anomaly range decreases with the increase of frequency, the focusing ability and the resolution have been improved. Comparing the instantaneous amplitude with the 25 Hz and 35 Hz MEWVD spectral characteristics, on the 45 Hz MEWVD profile, the resolution is higher, the longitudinal spectra have been separated, the transverse spectra have been narrowed, and the spectrum anomaly changes from "cluster" to "double-layer shape". This is more consistent with the low GR anomaly of the sand body, effectively revealing the top and bottom interfaces of the 35 m thin ancient channel sandstone in J₂s₃^3-2 in Well E. It can be seen that MEWVD spectrum anomalies can effectively describe the longitudinal distribution of narrow and thin ancient river channels.

In addition, MEWVD spectrum anomalies can characterize the horizontal distribution of microscopic ancient river channels. As shown in Fig. 6a, the anomaly with high instantaneous amplitude depicts relatively wide main river channels Ⅰ-Ⅳ, Ⅵ, Ⅷ. But in the black box, the channel is very narrow with thin channel sandstone, the seismic reflection energy is strong, and the resolution is low. Therefore, the instantaneous amplitude cannot be used to identify the spatial position, width, and flow direction of the narrow and thin ancient channel. In Fig. 6b, the 25 Hz spectrum obtained by the MEWVD method has high value anomalies, depicting the wide main river channels I-IV, VI, and VIII. Moreover, the low and medium anomalies portray the spatial distribution of narrower V and VII river channels. Especially in the black box on the left, the width and flow direction of V and VII river channels are very clear. Compared with microscopic river channels Ⅰ-Ⅷ, river channel Ⅸ is narrower with great change in flow direction, and is more concealed. It is difficult to identify with the instantaneous amplitude and 25 Hz frequency spectrum. From the 35 Hz and 45 Hz MEWVD spectrum, the layer characteristics of the meandering river channel IX are very clear. It can be seen from Fig. 6 that MEWVD has higher resolution than instantaneous amplitude, and MEWVD spectrum characteristics at different frequencies can identify narrow and thin river channels of different widths. Of course, since any object has a corresponding sensitive frequency band, river channels V and VII are more sensitive to low frequency, so the MEWVD spectrum characteristics at 35 Hz and 45 Hz are not as clear as that at 25 Hz (the black box in Fig. 6). But overall, high-frequency spectra have better characterization effect.

Using RGB fusion (red, green, and blue mixed to generate other colors), the MEWVD spectra at different frequencies merged can more clearly show the spatial distribution of river channels with different widths and thicknesses. As shown in Fig. 7a, the RGB fusion method was used to fuse the MEWVD spectra at 25, 35 and 45 Hz, and effectively depicted the distribution of ancient river channels Ⅰ to IX on the J₂s₃^3-2 sub-layer. Especially the extremely narrow river channels V, VII and IX, and more microscopic ancient river channels (indicated by red arrows) in the east of river channel Ⅷ and west of river channel Ⅵ, can also accurately reveal the width and flow direction and other spatial information. Of course, combined with the sedimentary structure in the area, using core analysis and logging interpretation, the reliability of microscopic ancient river channels identification by MEWVD spectrum anomalies can be confirmed. As shown in Fig. 7b, in the vertical direction, delta fronts micro-facies are developed in Well D, such as distributary channels, interdistributary bays, and estuary bars. Among them, the J₂s₃^3-2 sublayer is dominated by distributary channel sedimentary microfacies, with the AC and GR logging curves showing low values, and the deep lateral and shallow lateral resistivity curves showing high values. Cores taken at 2869.57 to 2871.49 m in the middle of the distributary channel show that the ancient channel has developed horizontal beddings. And under strong hydrodynamic action, a tight sandstone reservoir about 30 m thick was formed. Affected by lithology, physical properties, thickness and other factors, the impedance difference between the reservoir and the surrounding rock is not large. The "lens-like" river channel shown in Fig. 4c has weaker reflections, but there are prominent MEWVD spectrum anomalies in Figs. 6 and 7. It can be seen that considering seismic reflections, cores, logging response data, using MEWVD spectrum anomalies to describe microscopic ancient river channels is of high reliability.

In addition, MEWVD spectrum anomalies can also characterize the distribution of faults, but faults and ancient channels have significant MEWVD differences. As shown in Fig. 6d and Fig. 7a, the overall environment of sedimentary structure in Zhongjiang Gas Field is relatively stable. There are only a few near-east-west faults in the southern part of the study area, and no fractures are developed. The MEWVD spectral anomaly of the faults is stronger at the boundary between the hanging wall and the bottom wall, and weaker in the middle part. Although the faults and the ancient river channels are in slender strips, the frequency spectrum at the boundary of the ancient channel is weaker and the stronger in the middle. It can be seen that although the MEWVD spectrum anomalies of the faults and the ancient river channels are similar, there are still significant differences. In essence, the similarities and differences between the two are determined by the geological environment, and these seismic responses are accurately described by MEWVD anomalies.

In short, MEWVD is able to focus on frequency spectra and improve time-frequency resolution, and MEWVD at different frequencies has different sensitivity degrees to different channels. This effectively reveals spatial information such as the width, thickness and flow direction of the microscopic ancient river channels in the Shaximiao Formation in Zhongjiang Gas Field.

The spectral characteristics at the maximum entropy can best represent the spectral distribution of seismic signals. MEWVD not only avoids the cross-term interference generated by WVD, but also enhances the spectral focusing to the greatest extent, as such it effectively improves the time-frequency resolution of seismic signals.

Emulational seismic signals and forward simulation signals on narrow and thin models revealed that MEWVD can extract spectrum characteristics more accurately than STFT, GST, CWT, WVD, PWVD, SPWVD, CWD and other methods. Especially on numerical models with different width and thicknesses, they provide special spectral responses due to the difference in the maximum entropy. MEWVD at different frequencies can accurately identify narrow and thin geological bodies at different scales.

The microscopic ancient river channels identification method based on the maximum entropy criterion and Wigner-Ville distribution plays the full role of MEWVD such as strong spectral focusing and high time-frequency resolution. Microscopic ancient river channels at different scales have been effectively identified in Zhongjiang Gas Field. In addition, MEWVD at different frequencies has different degrees of sensitivity to different channels, so the spatial information such as the width, thickness and flow direction of the microscopic ancient river channels can be accurately described.

$\alpha $—autoregressive coefficient, also known as prediction error factor, dimensionless;

${{\alpha }_{m}}$—the autoregressive coefficient with serial number m, dimensionless;

A₁—the amplitude of attenuating component, dB;

A₂—the amplitude of stationary component, dB;

E—prediction error, W/Hz;

E_B,m—the backward prediction error with sequence number m, W/Hz;

E_F,m—the forward prediction error of serial number m, W/Hz;

f—the frequency of signal, Hz;

f₁—the frequency of attenuation component, Hz;

f₂—the frequency of stationary component, Hz;

GR—natural gamma, API;

H—entropy function, dimensionless;

i—unit of the imaginary number;

j—the serial number of autoregressive coefficient α_m, dimensionless;

k_n—autocorrelation function, dimensionless;

k_z—convolution kernel of k_n, dimensionless;

l—the serial number of the extended sampling point of discrete signal, dimensionless;

m—the serial number of α_m, E_F,m, E_B,m, w_n and other functions, dimensionless;

M—the order of AR model, dimensionless;

n—the serial number of discrete signal sampling point, dimensionless;

N—the length of discrete signal, dimensionless;

p—function of distribution density, dimensionless;

P—power spectrum, W/Hz;

P₀—initial power spectrum, W/Hz;

P_m—average power spectrum, W/Hz;

Q—absorption coefficient, dimensionless;

R_lld—deep lateral resistivity, Ω•m;

R_lls—shallow lateral resistivity, Ω•m;

t—time, s;

W—continuous WVD function, W/Hz;

W_z—discrete WVD function, W/Hz;

w_n—a sub-item of W_z, W/Hz;

X—signal, dimensionless;

x—signal in a certain period, dB;

x₁—attenuation component, dB;

x₂—stationary component, dB;

Δt—time sampling rate, s;

x_F—emulational signal, dB;

φ₁—the phase of attenuation component, rad;

φ₂—the phase of stationary component, rad;

τ—time delay, s;

z—discrete signal, dB.