Prediction of Shale Gas Reservoirs Using Fluid Mobility Attribute Driven by Post-Stack Seismic Data: A Case Study from Southern China

Paleozoic marine shale gas resources in Southern China present broad prospects for exploration and development. However, previous research has mostly focused on the shale in the Sichuan Basin. The research target of this study is expanded to the Lower Silurian Longmaxi shale outside the Sichuan Basin. A prediction scheme of shale gas reservoirs through the frequency-dependent seismic attribute technology is developed to reduce drilling risks of shale gas related to complex geological structure and low exploration level. Extracting frequency-dependent seismic attribute is inseparable from spectral decomposition technology, whereby the matching pursuit algorithm is commonly used. However, frequency interference in MP results in an erroneous time-frequency (TF) spectrum and affects the accuracy of seismic attribute. Firstly, a novel spectral decomposition technology is proposed to minimize the effect of frequency interference by integrating the MP and the ensemble empirical mode decomposition (EEMD). Synthetic and real data tests indicate that the proposed spectral decomposition technology provides a TF spectrum with higher accuracy and resolution than traditional MP. Then, a seismic fluid mobility attribute, extracted from the post-stack seismic data through the proposed spectral decomposition technology, is applied to characterize the shale reservoirs. The application result indicates that the seismic fluid mobility attribute can describe the spatial distribution of shale gas reservoirs well without well control. Based on the seismic fluid mobility attribute section, we have learned that the shale gas enrich areas are located near the bottom of the Longmaxi Formation. The inverted velocity data are also introduced to further verify the reliability of seismic fluid mobility. Finally, the thickness map of gas-bearing shale reservoirs in the Longmaxi Formation is obtained by combining the seismic fluid mobility attribute with the inverted velocity data, and two favorable exploration areas are suggested by analyzing the thickness, structure, and burial depth. The present work can not only be used to evaluate shale gas resources in the early stage of exploration, but also help to design the landing point and trajectory of directional drilling in the development stage.


Introduction
The recoverable resources of shale gas in China are about 25.1 × 10 12 m 3 [1] with the technological innovations of directional drilling and hydraulic fracturing. The Paleozoic marine shale in southern China provides the largest shale gas potential in China [2,3]. At present, almost all commercial shale gas production in China comes from the Paleozoic marine shale in Sichuan Basin, southern China [4]. Previous research has also predominantly focused on the Sichuan Basin [5][6][7][8]. To meet the growing demand for shale gas, it is necessary to expand the exploration scope beyond the Sichuan Basin and find new shale gas producing areas. The western Hunan-Hubei (WHH) fold belt, adjacent to the Sichuan Basin, contains rich black marine shale. However, the strong heterogeneous reservoir, complex geological structure, and low exploration degree in the WHH fold belt greatly increase the investment risks of shale gas drilling. Seismic attribute technology is a reliable method for describing the spatial distribution of hydrocarbon reservoirs, especially in the case with few drilling data [9,10]. It can help design the direction and trajectory of directional drilling, reduce the exploration risk, and optimize the production.
Seismic attributes commonly used for hydrocarbon detection include two categories. One is the seismic attribute based on rock physics theory, such as amplitude versus offset (AVO) attribute estimated from the pre-stack data through AVO inversion [11]. Nevertheless, in areas lacking enough well controls, AVO attribute may not be adequate to assist hydrocarbon identification [12]. The other is the seismic attribute based on seismic attenuation theory, such as low-frequency shadow [13] and high-frequency attenuation [14] calculated from post-stack data. Compared with the pre-stack attribute, the post-stack one has the advantages of no need for strict well control, low computational cost, and high signal-to-noise ratio [15]. However, the post-stack attributes mentioned above can directly reflect the magnitude of seismic attenuation, but not always directly reflect the property of hydrocarbon reservoir. It is necessary to find a universal seismic attribute indicating hydrocarbon directly. In recent years, some post-stack fluid mobility attributes have been applied to hydrocarbon identification in deep-water sandstone reservoirs [7,10,15,16]. Fluid mobility is the direct reflection of the rock permeability and fluid viscosity of hydrocarbon reservoirs. Silin et al. [17] and Goloshubin and Silin [18] approximated the low-frequency seismic reflection coefficient as a function of fluid mobility in fractured reservoir. Therefore, we extract fluid mobility from low-frequency seismic reflection data to directly indicate gas in (fractured) shale reservoir.
Extracting information from low-frequency seismic data is inseparable from spectral decomposition technology [19], which is a powerful tool for obtaining time-frequency (TF) spectrum of seismic data. The bilinear Wigner-Ville distribution (WVD) can achieve a higher-resolution TF localization than the linear spectral decomposition techniques [20], but suffers from the notorious cross-term interferences. The integration of WVD and MP algorithm can produce a stable high-solution TF spectrum without cross-term interferences [10,21]. MP is a greedy algorithm for compressive sensing [22]. However, the traditional MP algorithm is difficult to determine the precise TF locations of seismic signal due to the interference effect between different seismic frequency contents [23]. More specifically, the seismic low-frequency contents of interest may be suppressed by the seismic high-frequency contents during the process of MP algorithm, which may result in an average effect in the final TF spectrum and thereby probably reduces the accuracy of hydrocarbon detection. EEMD is a nonlinear technique for adaptively decomposing a seismic signal into a finite number of narrow-band signals [24]. Performing MP over the narrow-band signals can reduce the effect of frequency interference. Therefore, a spectral decomposition technology without frequency interference is proposed by modifying MP based on EEMD.
In this paper, the spatial distribution of shale gas reservoirs in complex structural areas lacking well control is predicted by using the fluid mobility attribute based on post-stack seismic data and spectral decomposition technique without frequency interference. We first propose a new spectral decomposition technology by modifying MP based on EEMD, and the synthetic data test shows that the new technology reduces the interference of different frequency content and obtains more accurate spectral information, especially in the low frequency domain. We then define a seismic fluid mobility attribute to detect shale gas by extracting fluid mobility measurements from the low-frequency post-stack seismic data, and real data application shows that the seismic fluid mobility attribute can describe the spatial distribution of shale gas reservoirs well. We also introduce inverted velocity to further verify the reliability of the seismic fluid mobility attribute. Finally, we predict the thickness map of shale gas reservoirs in study area by integrating seismic fluid mobility attribute and velocity data. This work not only can help design the direction and trajectory of the directional drilling, but can also be used for resource evaluation of shale gas reservoirs in the study area.

Target Overview
Commercial gas production has been achieved from the Lower Silurian Longmaxi shale in the eastern Sichuan Basin (ESB) [4,25]. In this paper, the research target is the Lower Silurian Longmaxi shale in the WHH fold belt, adjacent to the ESB ( Figure 1a). Structurally, the WHH fold belt is composed of a series of anticlinoria and synclinoria formed by multistage tectonic deformation. Compared with eastern Sichuan Basin, the WHH fold belt belongs to a tectonically active area, which results in the strong heterogeneity of shale gas distribution [26]. The Longmaxi Formation in both the WHH and ESB areas is integrated on top of the nodular limestone of the Ordovician Baota Formation and under the yellowgreen shale of the Luoreping Formation, which indicates their sealing capacities to be comparable [27]. The sedimentary environment in both WHH and ESB areas in the early Silurian was shallow water semi-closed retention basin, which made that the WHH and ESB Longmaxi shales have similar mineral composition, such as clay minerals accounting for 34.7% and 34.6% on average, feldspar accounting for 7.3% and 8.3% on average, and quartz accounting for 56.5% and 44.4% on average, respectively [25,27]. Field investigations and laboratory measurements have confirmed that the hydrocarbon generation ability and organic matter maturity of the Longmaxi shale in the WHH area are comparable to those in the ESB area [27]. The thicknesses of black Longmaxi shales with total organic carbon (TOC) greater than 2% range from 15 to 30 m in the WHH area, and these in the ESB area range from 20 to 50 m. The average kerogen reflectance of Longmaxi shales in the WHH and ESB areas are 2.68% and 2.69%, respectively, which indicates that they are both in the high maturity stage. In addition, the burial depths of the Longmaxi Formation in the WHH area range from 0 to 3300 m, which is basically in line with the favored burial depths (100-3500 m) for shale gas production [28]. Therefore, Longmaxi shale gas in the WHH fold belt is of great potential to exploration.
Specifically, the study area (blue rectangle in Figure 1a), covering an area of about 1189.7 km 2 , is located between the Yidu-Hefeng anticlinorium and the Sangzhi-Shimen synclinorium. Figure 1b is the elevation plan of the bottom interface of Longmaxi Formation in study area. Only the target formation in the northwest and middle of study area is retained due to erosion. The retained formation is characterized by well-developed faults and high-steep structures (elevation range −2500-0 m). Owing to erosion and fluctuating terrain, the surface elevation of study area varies dramatically between 0 and 1300 m, which is not conducive to three-dimension (3D) seismic acquisition. Now the study area is covered by twelve 2D seismic lines (dashed lines L1-L8, C1-C4 in Figure 1b). Six wells are drilled through the Longmaxi Formation in the study area, of which S-3 is a gas-free well, S-1 is a poor gas well, and the other four are gas wells (circles with different fill colors in Figure 1b). Deviated well SX-1 and vertical wells S-2, S-3, and S-4 are located near the seismic lines. Figure 2 displays the logging data of the vertical wells S-2, S-3, and S-4. The thickness of target formation is about 100 m, and its lithology is mainly carbonaceous shale with intercalations of argillaceous siltstone in the Longmaxi Formation. Here, shale is characterized by high acoustic time difference (AC) and high natural gamma (GR). Conversely, siltstone is characterized by low AC and low GR. As marked by red rectangles, the gas zones are located at the positions with high TOC value (about >2%). Therefore, TOC data can used as an indicator of a potential gas enrichment zone. TOC data from the vertical wells implies that the vertical gas enrichment zone is located near the bottom of Longmaxi Formation. Comparing the TOC data of three wells, we can also see that the lateral distribution of shale gas enrichment zones is different. The lateral gas enrichment zone in the Longmaxi Formation is the target of interest.  Logging data of wells S-2 (left), S-3 (middle) and S-4 (right). GR is the natural gamma curve, AC is the acoustic time difference curve, DEN is the density curve, CNL is the compensated neutron logging curve, PEF is the photoelectric absorption cross-section index curve, and TOC is the total organic carbon curve. The black dotted lines represent the interfaces of Longmaxi Formation, and the red rectangles denote the gas zones.

Spectral Decomposition Technology without Frequency Interference
The key of using low-frequency seismic attribute to characterize hydrocarbon reservoir is the TF resolution of spectral decomposition algorithm. MP integrated with WVD is a high TF resolution spectral decomposition technique without the limit of Heisenberg's uncertainty principle, where MP is the crucial signal decomposition algorithm and WVD is the spectral calculation algorithm. MP decomposes a seismic signal into a linear superposition of a series of pre-defined wavelet atoms.
where s is the (real-valued) seismic signal, n is the iteration number, c n is the complex coefficient at the nth iteration, m n is the optimal wavelet atom at the nth iteration, r N is the residual error at the Nth iteration, which is only noise and all significant features of seismic signal are recorded in the first term on the right-hand side of Equation (1), and Re{•} is the real part operator. The wavelet atom m n , containing the time and frequency features of seismic signal, is extracted from the atomic dictionary. Here, we choose the Morlet wavelet as the atom because it is suitable for frequency quantification and attenuation studies of seismic signal [21].
where t is the time (s), u n is the center time (s), k n is the (dimensionless) scale, f n is the center frequency (Hz), and i is the imaginary unit. The key parameters (u n , k n , f n ) of m n are estimated by the following equations [10] u n = arg max (k n , f n ) = arg max where r n−1 is the residual error at the (n-1)th iteration, r 0 = s + is H , superscript H is the Hilbert transform operator, superscript env is the envelope operator, • is the dot product operator, D k and D f is the searching interval of k n and f n , respectively. In general, D k and D f are given empirically, such as D k = [0. 3,6] and where f u is instantaneous frequency of r n-1 at the time u n , and R f is the search radius of f n , which is always set to 20 Hz [10,21]. The complex coefficient c n , containing the amplitude and phase features of seismic signal, is expressed as follows where superscript T is the transpose operator and ε is the damping constant. According to the WVD, the TF spectrum of seismic signal is expressed as [10] A where Im{•} is the imaginary part operator. As follows from Equation (4), the accuracy of f n affects the frequency resolution of MP algorithm. However, the interference of multi-frequency components in nonstationary seismic signals causes inaccuracy in searching f n , and then introduce error in the TF spectrum. To reduce the effect of frequency interference, a modified MP scheme based on EEMD (MP-EEMD) is proposed in this work. MP-EEMD mainly includes four stages. In stage 1, performing the EEMD algorithm over seismic signal yields where j is the serial number of intrinsic mode function (IMF), ϕ j is the jth IMF, and η is the residual. In stage 2, decomposed components with similar frequency-band characteristics are adaptively stacked according to the following principle where ϕ j+1 = η, l = 1, 2, . . . , L (L ≤ J + 1) is the serial number of the stacked components, is the effective frequency band of ϕ j , estimated through fast Fourier transform (FFT), and p is a constant threshold. Here, we set the threshold p to be 80%. In stage 3, the searching interval D l f for the stacked components ψ l is set to be where f l is the instantaneous frequency of ψ l at the time is the effective frequency band of ψ l , and r = f max l − f min l is the bandwidth of ψ l . In stage 4, we perform MP algorithm on the ψ l to obtain the final TF spectrum. After EEMD and adaptive stack, the narrow frequency-band component ψ l can reduce the effect of frequency interference of original seismic signal to some extent. It can be seen from Equation (9)  ), which give more constraints to search for optimal frequency parameter f n than MP.
To clearly illustrate the advantages of MP-EEMD, we design a synthetic signal with multi-frequency interference (blue in Figure 3a) by superimposing four wavelets W1-W4 with different frequency (black in Figure 3a). Figure 3b,c show the TF spectrum of the synthetic signal calculated by MP and MP-EEMD algorithms, respectively. As shown in Figure 3b, the spectral amplitude corresponding to the low-amplitude wavelet W1 (near dotted line) is vague and divergent due to frequency interference, and there are some significant deviations between the TF centers calculated by MP (white squares) and the true ones (black squares) of spectral amplitudes corresponding to the wavelets W2-W4. The larger the difference between interfering frequencies in synthetic signal, the larger the frequency deviation in TF spectrum. As shown in Figure 3c, the spectral amplitude corresponding to the weak wavelet W1 (near white dotted line) is clear and convergent, and the TF centers calculated by MP-EEMD (black points) are close to the true ones (white points) of spectral amplitudes corresponding to the wavelets W2-W4. Comparing Figure 3b,c, we can observe that the accurate and TF resolution of MP-EEMD is higher than that of MP. Table 1. Parameters of the four wavelets (Figure 3a).
In panels (b,c), the white dotted line is the true TF location of W1, the three white points are the true TF centers of W2-W4, the three black points are the TF centers of W2-W4 estimated by spectral decomposition algorithm, the cool colors represent small spectral amplitudes, and the warm colors represent large spectral amplitudes.
To further verify the feasibility of MP-EEMD, the seismic trace data near well S-2 ( Figure 4a) is selected for real data test. Figure 4b,c displays the components decomposed by EEMD and their corresponding frequency spectra obtained by FFT. It can be seen that the five components have different dominant frequency and frequency band. By comparing with the frequency spectra of the real seismic signal s and that of the five components, we regard the component ϕ 1 as a high-frequency (>200 Hz) low-amplitude noise and remove it in the subsequent steps. In practice, the frequency band of seismic signal is limited and narrow (about 10-60 Hz), so we set the threshold of dominant frequency as 80 Hz in this study. The components with dominant frequency greater than this threshold are considered to be noises, and the others are the valid signals. Therefore, MP-EEMD can eliminate certain noise from the real seismic signal by analysis the spectral features of the components decomposed by EEMD. Then, after adaptive stacking (Equation (8)), the four valid components are stacked to be two sub-signals (Figure 4d), which represent the high-frequency and low-frequency information of the real seismic signal, respectively. Figure 5 displays the TF spectra of the real seismic signal using MP and MP-EEMD. As indicated in the white rectangles, MP-EEMD can avoid the obvious average effect, resulting from frequency interference of ψ 1 .and ψ 2 , and obtain a more focused TF spectrum than MP. In addition, the dominant frequency of MP-EEMD-based TF spectrum decreases gradually from shallow to deep, which is conforms to the propagation theory of seismic waves. Therefore, the TF spectrum obtained by MP-EEMD performs better.

Seismic Fluid Mobility Attribute
A simplified asymptotic representation of low-frequency reflection coefficient R for the normal-incidence plane P-wave in a fluid-saturated media is expressed as [17,18] where m is the fluid mobility, R 0 and R 1 are dimensionless coefficients, which depend on neither frequency nor fluid mobility. Taking the derivative of Equation (10) with respect to frequency yields Equation (11) implies the fluid mobility is proportional to the product of the square of the first derivative over frequency of the reflection coefficient spectrum and the frequency. According to seismic convolution theory, the seismic amplitude spectrum A(f ) is proportional to the reflection coefficient spectrum R(f ) under the assumption of fixed seismic wavelet. Studies on the seismic amplitude spectrum change in fluid included [10,29] also show that fluid mobility reflects the first derivative of the seismic amplitude spectrum over frequency. Therefore, we define the seismic fluid mobility M as where ∆f is the frequency increment, and Ω is range of the calculational (low) frequency, which is estimated from the seismic trace data near wells and the logging interpretation data. Choosing a frequency range instead of a single frequency can improve the reliability of fluid mobility. Figure 6 displays the seismic amplitude spectra and the corresponding frequency gradients for different lithologies with different gas-bearing conditions near S-2 versus frequency. Noted that the spectrum amplitude of gas shale is higher than that of gas-free shale and siltstone. We also can note that the frequency gradients of the gas zone 13-23 Hz (blue rectangle) are also higher than that of the gas-free zone. This means that the seismic fluid mobility attribute can effectively distinguish the gas zone and the gas-free zone.   (12)) near S-2, S-3, and S-4. The seismic fluid mobility attribute extracted by using the MP-EEMD algorithm corresponds to the blue curve and extracted by using the MP algorithm corresponds to the blue dotted curve. From Figure 7, we can see that the high value of seismic fluid mobility attribute indicates the location of gas zone (red rectangle) and the low value corresponds to the gas-free zone. As shown with black arrows in Figures 8 and 9, the seismic fluid mobility based on the MP-EEMD algorithm (blue curves) can exhibit less interference for the gas-free zone at well S-3 and better indicate the gas zone at well S-4 than that based on the MP algorithm (blue dotted curves). These factors not only prove the effectiveness of seismic fluid mobility attribute, but also prove the advanced nature of the MP-EEMD algorithm.

Results
The 1D seismic fluid mobility attributes have indicated that the MP-EEMD-based seismic fluid mobility technology is effective in distinguishing the gas wells from the gas-free wells. Here, the MP-EEMD-based seismic fluid mobility attribute technology is applied to predict the spatial distribution of shale gas accumulation zone in Longmaxi Formation of the whole study area (indicated with blue rectangle in Figure 1a). As shown in Figure 1b, wells S-2, S-3, and SX-1 are located near the seismic line L4 (blue dashed line), and well S-4 is located near the seismic line C2 (green dashed line). Therefore, we next focus on showing and analyzing the seismic fluid mobility attribute extracted from seismic lines L4 and C2. Figure 10a displays the seismic section of seismic line L4 (blue dashed line in Figure 1b) near the gas wells SX-1 and S-2, Figure 10b displays the normalized seismic fluid mobility section extracted from Figure 10a, and Figure 10c displays the logging data of the deviated well SX-1. The actual well data in Figure 3 indicates a good shale gas display at the bottom section of well S-2, and the TOC curve in Figure 10c indicates a good and discontinuous shale gas display in the horizontal section (from A to B) of well SX-1. However, little information about shale gas can be seen from the post-stack seismic section (Figure 10a). It is clear from the seismic fluid mobility attribute section (Figure 10b) that the shale gas near S-2 and SX-1 is vertically concentrated at the bottom of the Longmaxi Formation, but laterally distributed in a discontinuous manner. Specifically, for the vertical well S-2, the gas-accumulation zone predicted from seismic fluid mobility attribute section is located at the bottom section with high GR and high TOC values (Figure 2). For the deviated well SX-1, the predicted gas-bearing condition from A to B has a trend of getting better, which is basically consistent with the horizontal variation of TOC curve in Figure 10c. Based on the predicted result, we have learned that the shale gas enrich areas in the study area are located near the bottom of the Longmaxi Formation. In addition, the seismic fluid mobility attribute is driven only by the post-stack seismic data, eliminating the need for excessive well control. Therefore, the spatial distribution of shale gas reservoirs predicted from seismic fluid mobility attribute can help researchers design the direction and trajectory of directional drilling, reduce the exploration risk, and optimize the production, even in the early stages of exploration and development. Same as the gas wells SX-1 and S-2, the gas-free well S-3 is also located near the seismic line L4 (blue dashed line in Figure 1b). Figure 11 displays the comparison between the seismic section and the corresponding seismic fluid mobility section connecting well S-3. Compared with seismic section, seismic fluid mobility attribute section clearly describes the distribution characteristics of shale gas. Predicted result based on seismic fluid mobility attribute indicates that the area near well S-3 displays no gas, while the deeper areas far away from well S-3 (indicated with red arrows in Figure 11b) display good gas-bearing properties. To further verify the accuracy of the above shale gas prediction result, we introduce another petrophysical parameter, i.e., velocity, through post-stack seismic inversion technology (see Appendix A for detailed algorithm). The inverted velocity section connecting well S-3 is shown in Figure 12. Here the blue color represents the low-velocity pure shale, the yellow represents the middle-velocity pure siltstone, the green represents the transition between shale and siltstone in the target formation, and the red color represent the high-velocity limestone outside the target formation. The inverted velocity section reveals that the target formation mainly contains two sets of low-velocity shale reservoirs (marked by black arrows in Figure 12), which is consistent with the lithology data in Figure 2. In addition, it can be seen from the TOC data of Figure 2 that the gas content of the bottom shale reservoir is greater than that of the top shale reservoir. We can also see that the predicted two shale gas zones (marked by red arrows) in Figure 11b are located at the bottom shale reservoir in Figure 12 by comparing the inverted velocity and the seismic fluid mobility.  S-4 is a gas well located near the seismic line C2 (green dashed line in Figure 1b). Figure 13 displays the comparison between seismic fluid mobility section and inverted velocity section connecting well S-4. Figure 13a indicates that the gas zone in well S-4 is located near the bottom interface of the Longmaxi Formation, which is consistent with the actual well data. Figure 13b describes the distribution characteristics of shale reservoirs with low-velocity, where the blue color inside the black ellipse represents a set of lowvelocity, thick, and pure shale, i.e., favorable shale reservoir. Comparing Figure 13a,b, we can see that the gas accumulation area predicted from the seismic fluid mobility attribute is located in the favorable shale reservoir interpreted from the inverted velocity data. These results confirm the validity and stability of the seismic fluid mobility attribute technique in identifying the spatial distribution of shale gas enrichment areas.  The thickness map of shale gas reservoirs is obtained in three steps: first, the time location of shale gas reservoir is estimated using a threshold value of seismic fluid mobility (i.e., m > 0.3), then the thickness data at the twelve seismic lines is calculated by using another threshold value of velocity (i.e., 3900 < v < 5300 m/s), and finally the thickness in the entire study area can be obtained by interpolation technology. As shown in Figure 14, well S-3 has no gas (zero thickness), well S-1 has poor gas content (thickness less than 6 m), and the remaining four wells have good gas content (thickness greater than 10 m), which is matched with the actual drilling result. In addition, the sedimentary facies of the target formation is the transition facies from shallow water shelf in the southeast to deep water shelf in the northwest of the study area. According to the sedimentology theory, with the increase of water depth, the thickness of shale deposition and the maturity of organic matter increase as well as the thickness of shale gas reservoir. Thus, predicted thickness map is consistent with sedimentary background of marine shale. We also give two favorable exploration areas (inside the cyan dotted curves) with large thickness, few faults, and large burial depth (more than 1000 m), which can be used to guide subsequent exploration of shale gas in the Longmaxi Formation.

Conclusions
We have predicted the spatial distribution of shale gas reservoirs in Longmaxi Formation of the study area in Southern China using a frequency-dependent seismic fluid mobility attribute, extracted from the post-stack seismic data through a novel spectral decomposition technique (i.e., MP-EEMD scheme), and we also have demonstrated its applicability for the shale reservoirs in complex structural area lacking well control. There are some significant deviations between the TF locations estimated by using the traditional MP algorithm and the true ones due to the interference of multi-frequency components in seismic signal. The synthetic data test indicates that the larger the difference between interfering frequencies, the larger the deviation between the estimated and true TF locations, and the low-amplitude component cannot be recognized on the TF spectrum. These deviations in the TF spectrum will affect the computational accuracy of frequency-dependent seismic attribute. Based on the EEMD algorithm, we propose a novel MP-EEMD scheme to minimizes the effect of frequency interference. The synthetic data test indicates that the MP-EEMD scheme provides a TF spectrum with high accuracy and resolution, which also recognize the weak component. The real data test indicates that the MP-EEMD scheme not only achieves a more focused time-frequency spectrum, but also adaptively removes noise from real seismic data. By introducing the asymptotic low-frequency reflection coefficient, the seismic fluid mobility attribute can be extracted from the low-frequency post-stack seismic data through the spectral decomposition technique. The 1D attributes near wells show that the seismic fluid mobility based on the MP-EEMD algorithm exhibits less interference at the gas-free zone and better indication at the gas zone than that based on the MP algorithm. Thus, the MP-EEMD-based seismic fluid mobility attribute technology is applied to characterize the shale reservoirs with complex subsurface structure and low exploration level in Southern China. The spatial distribution of shale gas reservoirs predicted from the seismic fluid mobility attribute section is consistent with actual well data. Based on the predicted result, we have learned that the shale gas enrich areas in the study area are located near the bottom of the Longmaxi Formation. Referring to the predicted spatial distribution pattern of shale gas reservoirs in the area can help engineers design the directions and trajectories of directional drilling during the development stage. By combining the seismic fluid mobility attribute with the inverted velocity data, the thickness map of shale gas reservoirs in target formation is obtained, which is consistent with sedimentary background of marine shale. Referring to the thickness map of shale gas reservoirs in the area can help researchers to evaluate shale gas resources in the early stage of exploration. We have also given two potential area of future research by analyzing the thickness, structure, and burial depth. This work can provide some references for the subsequent exploration and research of shale gas reservoirs. the Longmaxi Formation, the number of time sampling points during the inversion process is about 50 (Figures 10a and 11a). Due to these properties, we approximate density of Longmaxi Formation to be a constant and then the Equation (A3) can be can be converted to be r(v) = r 1 , r 2 , . . . , r g−1 where v = v 1 , v 2 , . . . , v g T is the velocity vector.
Due to limitation in the frequency bandwidth of seismic data, directly solving for the objective function (13) is strongly unstable. Based on the Tikhonov regularization theory and the Bayes theory, a stable and reliable impedance inversion result can be obtained [30]. The impedance update during iterative process of seismic inversion is expressed as v (y+1) = v (y) + G (y) T G (y) + λ (y)