• <tr id="yyy80"></tr>
  • <sup id="yyy80"></sup>
  • <tfoot id="yyy80"><noscript id="yyy80"></noscript></tfoot>
  • 99热精品在线国产_美女午夜性视频免费_国产精品国产高清国产av_av欧美777_自拍偷自拍亚洲精品老妇_亚洲熟女精品中文字幕_www日本黄色视频网_国产精品野战在线观看 ?

    A two?step method to apply Xu–Payne multi?porosity model to estimate pore type from seismic data for carbonate reservoirs

    2020-06-22 06:04:36HongBingLiJiaJiaZhangShengJuanCaiHaoJiePan
    石油科學(xué)通報(bào) 2020年3期

    Hong?Bing Li · Jia?Jia Zhang · Sheng?Juan Cai · Hao?Jie Pan

    Abstract Carbonate reservoirs exhibit strong heterogeneity in the distribution of pore types that can be quantitatively characterized by applying Xu–Payne multi-porosity model. However, there are some prerequisites to this model the porosity and saturation need to be provided. In general, these application conditions are difficult to satisfy for seismic data. In order to overcome this problem, we present a two-step method to estimate the porosity and saturation and pore type of carbonate reservoirs from seismic data. In step one, the pore space of the carbonate reservoir is equivalent to a single-porosity system with an effective pore aspect ratio; then, a 3D rock-physics template (RPT) is established through the Gassmann’s equations and effective medium models; and then, the effective aspect ratio of pore, porosity and fluid saturation are simultaneously estimated from the seismic data based on 3D RPT. In step two, the pore space of the carbonate reservoir is equivalent to a triple-porosity system. Combined with the inverted porosity and saturation in the first step, the porosities of three pore types can be inverted from the seismic elastic properties. The application results indicate that our method can obtain accurate physical properties consistent with logging data and ensure the reliability of characterization of pore type.

    Keywords Carbonate reservoirs · Pore type · Physical properties · Rock-physics modeling · Rock-physics inversion

    1 Introduction

    A signi ficant portion of the world’s oil and gas reserves resides in carbonate rocks (Alotaibi et al. 2010; Bagrintseva 2015). Unlike conventional siliciclastic rocks, carbonate rocks undergo complicated depositional and diagenetic processes. They have strongly heterogeneous pore system,such as vugs, moldic, intergranular, intragranular, intercrystal pores and micro-cracks. Complex pore structure has a very signi ficant effect on elastic properties. Laboratory measurement results have demonstrated that the difference in pore type can lead to a great difference in P-wave velocity of up to 40% at certain same porosity (Anselmetti and Eberli 1993; Xu and Payne 2009). Such high heterogeneity in these porous rocks brings great uncertainties for carbonate reservoir prediction and makes oil and gas recovery difficult during development, which has been identi fied as the crucial parameter needed to resolve the scatter in seismic and petrophysical relationships (Anselmetti and Eberli 1999; Weger et al. 2009; Sun et al. 2015; Jin et al. 2017). Therefore, the pore-type estimation from seismic data is helpful to describe the heterogeneity of carbonate reservoirs, predict relatively high-quality reservoir, and further improve the accuracy of reservoir prediction.

    Quantitative estimation of pore type from geophysical measurements is frequently attained by using various rock-physics models that relate rock elastic properties to geophysical observations (Kumar and Han 2005; Xu and Payne 2009; Zhao et al. 2013; Yuan et al. 2018). The idea of an effective pore aspect ratio has been used in rockphysics modeling and shear velocity prediction for complex reservoirs (Sams and Andrea 2001; Li et al. 2013).Sams and Andrea (2001) found that the results of numerical modeling with a spectrum of pore aspect ratios can be matched using the same porosity and a single aspect ratio in a clastic reservoir by numerical simulation. Li et al.(2013) applied the differential effective medium (DEM)model to invert the effective pore aspect ratio from P-wave velocity and used it to predict shear wave velocity. Xu and Payne (2009) developed a new model to estimate the porosities of various pore types from well-logging data under the condition that the reservoir porosity and saturation are provided. These application conditions usually can be met for well-logging data because the porosity and saturation might be estimated from other logging data such as resistivity and neutron data. However, for seismic data there are not enough data to be applied. Some researchers attempt to apply Xu–Payne model (Xu and Payne 2009) to seismic data. Zhao et al. (2013) proposed to apply a neural network inversion method to only estimate the porosity and then estimate pore type in carbonate heavy oil reservoirs from seismic data. However, seismic elastic properties such as P-wave velocity are a function of porosity,saturation and pore type. That is, these reservoir parameters together have a comprehensive impact on the elastic parameters; therefore, one of them such as porosity cannot be recovered alone without concerning the others. Adelinet and Le Ravalec (2015) applied the effective medium theory to make an inference about the porosity and fracture from seismic wave impedance. Their method assumes that carbonate pore space contains two types: One is equant pore, and the other is crack, and these two porosity inclusions are filled with water, and then takes advantage of the available impedance to quantify the spatial variations in the equant porosity and fracture density. Therefore, there are some limitations in the current methods. For example,some methods assumed that pores are filled with water or heavy oil. Some methods simpli fied the pore system. At present, few scholars have proposed a solution to simultaneously invert the porosity, saturation and pore-structure parameters by means of the Xu–Payne model.

    In this paper, we present a two-step method to estimate porosity and saturation as well as pore-type distributions from seismic data by using only seismic rock-physics models. Firstly, we summarize our rock-physics modeling method for complex carbonate reservoirs. Then, we design a work flow on how to give the equivalent way dealing with the pore space of carbonate reservoirs and estimate the reservoir parameters. By treating the pore space of carbonate reservoirs as a single-porosity system, we simultaneously estimate the porosity, fluid saturation and effective pore aspect ratio; after that, by treating the pore space of carbonate reservoirs as a triple-porosity system, we apply Xu–Payne model to estimate the porosities of three pore types from seismic data; furthermore, by numerical experiments on theoretical models we prove that the elastic properties of porous rock can be equivalently simulated by a singleporosity system with an effective pore aspect ratio. Finally,the results of estimating pore type of a dolomite reservoir from real data are given to demonstrate the effectiveness of the proposed method.

    2 Rock?physics modeling in carbonate reservoirs

    2.1 Xu-Payne Model

    Xu–Payne model (Xu and Payne 2009) has been widely used in rock-physics modeling and inversion on carbonate rocks(Zhao et al. 2013; Neto et al. 2014). It partitions the pore space of carbonate rock into clay and non-clay pores. The non-clay pores are divided into three types of elliptical pores with different aspect ratios. These three representative pore types are de fined as follows: (a) Interparticle and intercrystal pores are de fined as reference pores; they have moderate pore aspect ratios which are between 0.12–0.15; (b) microfractures and micro-cracks are de fined as compliant cracks;they have lower aspect ratios which are usually below 0.02;(c) moldic pores and vugs are de fined as stiff pores; they have higher aspect ratios which are usually over 0.7. The porosity of clay pores that are regarded as micropores with bound water is the product of the total porosity and the shale volume.

    Xu and Payne (2009) proposed a rock-physics modeling work flow for their model. There are two main aspects. One is to estimate the elastic moduli of dry rock frame using the effective medium model. The other is to calculate the velocities of low-frequency, water-saturated rock using Gassmann’s formulae (Mavko et al. 2009). When clay micropores exist, they are added to the matrix using the differential effective medium as a part of “solid” (Xu and Payne 2009).

    2.2 Dry frame moduli for porous rock

    The elastic properties of carbonate rocks have been simulated using the effective medium models, such as Kuster–Toks?z theory and DEM theory, with different concentrations of pores and speci fied pore shapes or pore types (Kumar and Han 2005; Sayers 2008; Xu and Payne 2009; Zhao et al. 2013; Li and Zhang 2018). The classical DEM process and the Kuster–Toks?z theory applied in the original Xu–Payne model are used to calculate effective elastic moduli for dry rocks (Kuster and Toks?z 1974; Xu and Payne 2009). Compared with the Kuster–Toks?z theory,the DEM theory simulates the effective elastic properties of the composite more accurately and efficiently even if the volume percentages of the inclusion are large (Norris 1985;Berryman and Berge 1996; Li and Zhang 2014). At present,the DEM theory has been extended from the single-porosity medium to the multiple-porosity medium without the problem of the added order of the inclusions (Li and Zhang 2014). The ordinary differential equations for dry porous rocks can be written as (Li and Zhang, 2014):

    whereyis porosity,K?(y)andG?(y)are, respectively, the effective bulk and shear moduli of the composite,viis the volume percentage of theith inclusion in total porosityφ,P?iandQ?iare called as the polarization factors of the effective bulk and shear moduliK?andG?, respectively; they depend on the aspect ratios of the inclusions (Mavko et al. 2009).Nis the number of the inclusions. WhenN=1, they are the same as the classical DEM theory for dry rocks (Mavko et al.2009; Li and Zhang 2014). WhenN=3, they become the with valuesK?(0)=KmandG?(0)=Gmto the desired highest valuey=φ, but also from their analytical approximation formulae (Li and Zhang 2014).

    3 Rock?physics inversion in carbonate reservoirs

    When we use Xu–Payne model to estimate the pore type of carbonate reservoirs, we must firstly de fine the reference line(Xu and Payne 2009) which is calculated by applying rockphysics model and the known porosity and saturation under the assumption that the whole pore spaces of porous rock contain only reference pores. It means that we have to provide the porosity and saturation for Xu–Payne model when it is applied to seismic data. Since the complex pore system of the carbonate reservoir makes the relationship between the velocity and porosity widely scattered (Rafavich et al. 1984;Eberli et al. 2003; Sayers 2008; Baechle et al. 2009), it would lead to a large error in the estimation of porosity and saturation if the effect of pore structure was not taken into account in the adopted rock-physics model (Li and Zhang 2018). Li and Zhang (2018) proposed an inversion scheme for porosity and saturation to solve this problem based on three dimensions rock-physics template (3D RPT). Here, we draw on their research ideas and present a two-step work flow that triple-porosity model which is suitable to predict the effective elastic properties for Xu–Payne model.

    The effective elastic properties for dry rocks might be obtained not only from the numerical solutions of differential Eqs. 1 and 2 by integrating them from porosityy=0 only applies rock-physics model to estimate pore-type distribution of carbonate reservoirs from seismic data (Fig. 1).In step one, the pore space of carbonate reservoirs is treated as a single-porosity system with the effective aspect ratio of pore; then a 3D RPT on P- and S-wave impedances and density with the effective aspect ratio of pore is established through the Gassmann’s equations and effective medium models, porosity and fluid saturation; then, the effective aspect ratio of pore, porosity and fluid saturation are simultaneously estimated from the seismic data based on 3D rockphysics template (Li and Zhang 2018). In step two, the pore space of carbonate reservoirs is treated as a triple-porosity system just as the Xu–Payne model de fines; then combined with the inverted porosity and saturation obtained in the first step, three types of pores which include reference pores, stiffpores and cracks can be inverted by the means of Xu–Payne method. The seismic elastic data in Fig. 1 generally include P- and S-wave velocities and density, which are inverted from pre-stack gather data based on the amplitude variation with offset (AVO) theory.

    3.1 The effectiveness on porous rock taken as a single?porosity system

    Assuming that the weighted summation of the polarization factorsP?iandQ?iof inclusioniwith the pore aspect ratioαiin Eqs. 1 and 2 can be approximated by the average polarization factorsP?IandQ?Iwith an effective pore aspect ratiothat isthe differential Eqs. 1 and 2 become:

    Fig. 1 Diagram work flow of a two-step method to estimate porosity and saturation as well as pore-type distribution from seismic data for porous rock

    Fig. 2 Comparison of simulation results for elastic moduli between Eqs. 1–2 and Eqs. 3–4. a Bulk modulus versus porosity for dry rock; b shear modulus versus porosity for dry rock; c bulk modulus versus porosity for fully saturated oil rock; d bulk modulus versus porosity for fully saturated water rock

    The above equations indicate that the elastic properties for porous rocks can be simulated by DEM equations for a single-porosity system with an effective pore aspect ratioαI. We give a comparison of simulation results between Eqs. 1–4. Here, we imitate the Xu–Payne model and take dolomite as the host medium, letKm=89GPa andGm=37GPa. We assume that there are three pore types in the rock, which are the stiff pores with the pore aspect ratio of 0.9, the reference pores with the pore aspect ratio of 0.15 and the micro-cracks with the pore aspect ratio of 0.005, respectively. We also assume that three different materials exist within the pore spaces, which are oil,water and dry void, respectively. Then, the fourth-order Runge–Kutta scheme can be used to get numerical solutions of the differential Eqs. 1–4. The step size here used in numerical integration is set to Δy=0.1. In Fig. 2, the scattered points stand for the simulated elastic moduli using Eqs. 1 and 2 for different pore volumes. The red points represent the simulated values for the volume contents of the stiff pores, the reference pores and micro-cracks,which are 89.9%, 10% and 0.1%, respectively, indicating that the stiff pores in the pore space of porous rock are dominant, while the blue points represent the simulated values for the volume contents, which are 10%, 89.9% and 0.1%, respectively, indicating that the reference pores in porous rock are dominant; The pink points represent the simulated values for the volume contents, which are 45%,50% and 5%, respectively, indicating that micro-cracks in porous rock pores are dominant. The curves represent the simulated elastic moduli using Eqs. 3 and 4 for different effective pore aspect ratios. The red curve stands for the simulated values for the effective pore aspect ratioαeff=0.42; the blue curve stands for the simulated values forαeff=0.155; the pink curve stands for the simulated values forαeff=0.07. Note that these curves for both dry rocks and the water- and oil-saturated rocks match well with the scattered points. It is fully illustrated that the elastic properties of multiple-porosity medium can be simulated by an equivalent single-porosity medium with an effective aspect ratio.

    3.2 Porosity, saturation and effective pore aspect ratio inversion based on a single?porosity system

    The simultaneous inversion of porosity and saturation as well as effective pore aspect ratio can be carried out by using the 3D rock-physics templates (Li and Zhang 2018). Li and Zhang (2018) give concrete steps to realize the inversion:First determine the distribution ranges of the porosity, saturation and pore aspect ratio from logging data in the study area,and generate the 3D mesh in the reservoir parameters space,which is made up of the porosity, saturation and pore aspect ratio; ensure that all the grid nodes are assigned to a speci fic saturation, porosity and pore aspect ratio; then, employ the DEM Eqs. 3–4 and fluid distribution model for each grid node to compute its bulk and shear moduli of dry rock and fluid moduli under the given pore aspect ratio, porosity and saturation, and then calculate the P- and S-wave velocities and density of saturated rocks by means of the Gassmann’s equations, until all the nodes in the grid have been calculated;after that, build a 3D RPT in the elastic parameters space such as Ip-Is-ρ or λρ-μρ-ρ; superpose the elastic parameters measured by logging tools onto the built 3D RPT to verify the validity of the 3D RPT; finally, project these elastic parameters inverted from seismic gather onto the newly built 3D RPT, and use the least square method to search the closest grid point for the inverted elastic data; choose the pore aspect ratio and porosity as well as saturation corresponding to the closest mesh point as the inversion results of reservoir results.

    3.3 Pore?type inversion based on a multi?porosity system

    Once the porosity and saturation are determined, the volume percentages of three pore types (reference pores,stiff pores and cracks) can be estimated from the seismicinverted elastic properties by following the inversion scheme of Xu and Payne (Xu and Payne 2009; Zhao et al. 2013).The key inversion steps include: Firstly, de fine the aspect ratios of the three pore types that are reference pores, stiffpores and cracks in the carbonate reservoir. Then, assume that the whole pore spaces of carbonate rock contain reference pores alone and apply the DEM analytical models for single-porosity system and the Gassmann’s equations to compute P-wave velocitya t given porosity, saturationand aspect ratio.Next, determine whetherthe pore type of reservoir belongs to crack pore type or cavity pore type by comparing the measured P-wave velocity Vp with the estimated P-wave velocityFinally, apply nonlinear global optimization algorithm to find the best estimate for the pore volume percentages of both cracks and reference pores or both cavities and reference pores by minimizing the error between theoretical predictions and geophysical measurements based on the dual-porosity system.

    Fig. 3 Seismic pro file used for inversion tests in the study area

    4 Real application

    4.1 Survey of the study area

    Our work flow has been applied in porosity, saturation and pore-type prediction of Cambrian Longwangmiao formation dolomite reservoirs in Anyue gas field in southwestern China’s Sichuan Basin. Both matrix porosity and permeability of this dolomite reservoir are low. The porosity and permeability of the core are 2–15% and 0.001–10 mD, respectively. The reservoir porosity has strong heterogeneity. A large number of cavities and micro-cracks are developed in the reservoir. Pore types in the study area are classi fied into three types: intergranular and intragranular pores, dissolved cavities and fractures (Zhou et al. 2015; Li and Zhang 2018).Figure 3 shows a seismic pro file used for inversion tests in the study area. There are three wells in this cross-well section. Both wells X8 and X12 are industrial gas wells, while well X203 is a water-bearing well. They all locate at locally high points.

    4.2 RPT analysis of well log data

    Figure 4 shows the analysis of 3D rock-physics template on well X12. The parameters and display modes used in the rock-physics modeling are same as what Li and Zhang(2018) did. They took the bulk and shear moduli of rock matrix asKm= 89 GPa,Gm= 37 GPa, respectively, and the density of rock matrix is taken asρm= 2.8 g/cm3, and the bulk modulus of water is taken asKw= 2.2 GPa. The blue surface represents the trend of P- and S-wave impedances and density with the pore aspect ratio and porosity for water saturationSw= 100%, and the brown surface represents the trend of P- and S-wave impedances and density with the pore aspect ratio and porosity for water saturationSw= 0%. The yellow curve plotted on each surface denotes the one whose porosities are constant, and the light blue curve denotes the one whose pore aspect ratios are constant. Light green surface indicates the trend of P- and S-wave impedances and density with porosity and saturation for the effective pore aspect ratioα= 0.1. The blue surface for fully watersaturated rock is always situated above the brown for fully gas-saturated rock. This may due to the fact that the density of the water-bearing rock is always greater than that of the gas-bearing rock. The points on the cross-plot are color coded with water saturation in Fig. 4a and porosity in Fig. 4b, respectively (Li and Zhang 2018). Note that most of the samples whose water saturations that interpreted from logging data are below 20% situate near the brown surface,and are drawn in red colors. Only very few samples whose water saturations that inferred from well-logging data are more than 80% situate near the blue surface, and are drawn in blue colors. The reservoirs with high gas saturation have high porosities. More 3D RPT analysis results for other two wells which come from the same the study area are found in the article of Li and Zhang (2018). These results show that the predictions from the 3D RPTs are fully coincident with the log interpretation and real drilling results. Hence, they can be effectively used to characterize the P- and S-wave impedances and density with porosity and water saturation and pore aspect ratio (Li and Zhang 2018).

    Figure 5 shows the inversion results of the pore aspect ratio, porosity and saturation from logging data based on the 3D RPT in Fig. 4. In Fig. 5, the black curves indicate the measured P- and S-wave velocities, density, and the inferred porosity and saturation from logging data, respectively. The red curves represent, respectively, the inverted effective pore aspect ratio, porosity and water saturation. We also display the lithology derived from the field by the petrophysicist in the figure. It shows that most of the reservoir is dolomite,with only a small amount of clay whose volume content is almost constant and a very small amount of quartz minerals at its top and bottom. Observe that the porosity of well X12 is high, the maximum porosity is up to 15%, the gasbearing property is good, and the water saturation is generally less than 20%. Both the inverted porosities and water saturation based on the 3D RPT are in well coincidence with the inferred ones from logging data, and also coincide with those from the actual drilling. The inverted effective pore aspect ratios for high porosity reservoir are around 0.2.

    Fig. 4 Analysis of 3D rock-physics template on logging data for dolomite reservoir in gas well X12 with data points is coded with water saturation (a) and porosity (b), respectively

    Fig. 5 Depth pro files of Vp, Vs, density, gamma ray (GR), resistivity, neutron (CNL), clay content (Vsh), porosity, water saturation (Sw) from logging data and the predicted porosity, water saturation and pore aspect ratio (Pore AR) based on 3D rock-physics templates for dolomite reservoir in gas well X12 in Anyue gas field

    Fig. 6 Pore-type inversion results for well X12. Area of red, green and blue colors represent volume concentrations of crack-like, reference and stiff pores in the pore space, respectively

    Figure 6 shows the pore-type prediction results for gas well X12 following the solution of Xu and Payne (Xu and Payne 2009). The blue colors represent the inferred stiff pores, the green colors represent the inferred reference pores, and the red colors represent the inferred cracks. It can be seen that there are remarkable variations in both pore types and the corresponding effective pore aspect ratiosαeff, which are shown in Fig. 4 among samples. Reference pores account for a high proportion of the dolomite reservoir pore space in well X12.The rounded pores also dominate the pore space in the upper dolomite unit, and crack-like compliant pores are not well developed, only distribute in the lower dolomite unit. On the right side of Fig. 6, three thin slices at depths of 4622.1, 4644.5 and 4666.7 m are also shown. Note that the inverted results are consistent with thin slice observations. For example, there are several large pores in the top thin slice, and both the inverted reference porosity and stiff porosity indicate that there are high porosities at a depth of 4622.1 m. There are several cracks in the bottom thin slice, and the inverted crack porosity indicates that there are high porosities at a depth of 4666.7 m.

    4.3 Quantitative estimation of pore type from seismic data

    Figure 7a–c shows the prediction results of the porosity, saturation and effective pore aspect ratio from the inverted seismic elastic data through wells X12, X203 and X8 based on the 3D RPT for a single-porosity system. In addition to well X12 mentioned above, well X8 is also an industrial gas well, while well X203 is a water-bearing well whose reservoir has low porosity (Li and Zhang 2018). Figure 7a–c shows the inferred porosity, saturation and pore aspect ratio pro files, respectively.It can be found that both industrial gas wells X8 and X12 situate at locally high points and their porosities and gas saturations are relatively high, while water well X203 situates at the secondary high point and its porosities and gas saturations are relatively low. It is obvious that the pore aspect ratios of industrial gas well X8 are large, which indicates the pores are relatively harder, while the pore aspect ratios of well X203 are lower, which means the pores are relatively softer. The pore aspect ratios of well X12 are smaller than those of well X8.

    Figure 8a–c demonstrates the pore-type predictions from the inverted elastic properties combining with the inverted porosity and saturation displayed in Fig. 7a, b. One can observe that the dolomite reservoirs are extremely strong and heterogeneous. The inferred reference porosity distribution as shown in Fig. 8a indicates that the interparticle pores are the dominant pore type in this dolomite reservoir. The inferred stiff porosity distribution as shown in Fig. 8c indicates that the vuggy and moldic pores are mainly developed around wells X8 and X12. The porosities of reference pores and stiff pores of industrial gas well X12 are relatively larger than those of well X8,whereas the inferred compliant porosity distribution as shown in Fig. 8b indicates that cracks are not well developed in the whole cross-well pro file, only developed around well X203.

    Fig. 7 Sections of the predicted porosity (a), water saturation (b) and effective pore aspect ratio (c) from seismic data for dolomite reservoir

    5 Discussion

    It is important to note that the proposed method to infer the porosity and saturation as well as volume percentage or porosity of each pore type in pore system from seismic data is achieved by employing the Gassmann’s equations and the differential effective medium theory of multiple-porosity system. It is generally known that for a porous system,Gassmann’s equations have the following basic assumptions(Mavko et al. 2009; Adam et al. 2006): (a) The rock is macroscopically homogeneous and isotropic porous medium; (b) the pore space is connected; (c) a frictionless fluid is filled in the pores; (d) the rock- fluid system is closed; (e) when the rock is excited by a wave, the relative motion between the solid rock and the fluid is negligibly small compared to the motion of the whole saturated rock itself; (f) the pore fluid does not interact with the solid in a way that would soften or harden the frame. Therefore, when using Gassmann’s equations, the above assumptions must be taken into account. The assumption (a) on macroscopic homogeneity and isotropy of rock ensures that the wavelength is long compared to the grain and pore sizes. Most rocks can generally satisfy this assumption for seismic frequency band (20–200 Hz). On the other hand,the DEM theory assumes that the pores in porous medium disconnected each other and fluids cannot move through pores;this approach simulates very high-frequency saturated rock behavior suitable for ultrasonic laboratory condition (Mavko et al. 2009). However, some researchers have testi fied that there are almost no differences between elastic properties of rocks compared to their measurements at high and low frequencies in laboratory experiments (Mavko et al. 2009; Adam et al. 2006). For low-frequency seismic condition, a practical solution is to calculate first the dry rock bulk and shear moduli by applying effective medium model, then use the Gassmann’s equations to calculate the bulk and shear moduli of saturated rocks (Mavko et al. 2009; Kumar and Han 2005; Xu and Payne 2009; Zhao et al. 2013; Adelinet and Ravalec 2015).

    Carbonate rocks have strong elastic frames and very complex pore systems compared to siliciclastic rocks. Because a signi ficant number of pores in carbonate rocks is not well connected, carbonate rocks not always meet Gassmann’s basic assumption on pore connectivity. The applicability of Gassmann’s theory in carbonate rocks has been targetedly investigated, and an inspiring conclusion from those researches is that Gassmann’s equations are appropriate to simulate the saturated bulk and shear moduli of carbonate rocks (Sharma et al. 2006; Sharma et al. 2011; Grechka 2009; Gegenhuber and Pupos 2015; Yan and Han 2016).By using numerical tests, Grechka (2009) demonstrated that disconnected porosity is unlikely to cause signi ficant errors in in fill substitution when the pore aspect ratios are greater than about 0.2. If the pores are aligned or the bulk and shear moduli of the substituted in fills are similar, the limits of practical applicability of Gassmann’s predictions extend to smaller aspect ratios. Therefore, theoretically, Gassmann’s assumption (c) on pore connectivity is not a necessary condition (Grechka 2009; Yan and Han 2016).

    Fig. 9 Sections of the inverted P-wave impedance (a), S-wave impedance (b) and density (c) which are used to infer the reservoir parameters

    Our method to invert for pore types of reservoir from seismic data is based on P- and S-wave impedances and density inverted from pre-stack seismic data. Therefore, the quality of pore-type estimation depends strongly on the accuracy of pre-stack seismic inversion which is a relatively mature technique and has been widely applied in seismic reservoir prediction (Fatti et al 1994; Buland and More 2003; Hampson et al. 2005; Wang et al. 2019). It is generally believed that P- and S-wave data can be unambiguously resolved by pre-stack seismic inversion, but density, which is crucial for hydrocarbon detection, cannot (Tarantola 1986; Debski and Tarantola 1995). This is mainly because the amplitude is not sensitive enough to density, the offset limited data tend to result in unstable results of the three-parameter AVO inversion, and small noise may lead to a huge inversion error(Tarantola 1986; Debski and Tarantola 1995). In order to ensure the reliability of density inversion, a large number of studies have been carried out on the correlation between P-wave impedance and both S-wave impedance and density.Hampson et al. 2005 proposed a new method to invert P- and S-wave velocities and density from pre-stack seismic data simultaneously by assuming the linear relationship between the logarithm of P-wave impedance and both S-wave impedance and density. Their method is effective in both model data and actual seismic data Swan (2007) proposed a background trend correction method based on the correlation coefficient between P- and S-wave impedance. The P-wave re flectivity obtained by AVO inversion is used to correct the S-wave re flectivity obtained together, so that the correlation between them can match with the correlation of theoretical models or log data. Following their ideas, we correlate the inverted S-wave impedance and density by using the P-wave impedance obtained simultaneously to make sure that the correlation coefficients between the inverted P-wave impedance with both S-wave impedance and density are matched with the corresponding correlation coefficients in well log data. Figure 9 shows the result of pre-stack threeterm simultaneous inversion. The P- and S-wave impedances and densities of the two gas-bearing wells in the target interval (between the horizons E1g and E1l) have relatively low values, while the corresponding values of the water wells are relatively high. Figure 10 shows the cross-plots between P-wave impedance and both S-wave impedance and density for logging data of the three well in target strata and seismic inversion data adjacent to the borehole. The correlation coefficients among the inverted seismic data are close to those among the well-logging data. Our inversion method ensures the accuracy of the inverted reservoir parameters through correlation coefficient matching correction.

    Fig. 10 Cross-plots between P-wave impedance and S-wave impedance for log data (a), seismic inversion data (b) adjacent to the borehole, and between P-wave impedance and density for log data (c), seismic inversion data (d) adjacent to the borehole

    6 Conclusions

    In this paper, we proposed a new two-step method to estimate the porosity and saturation and quantitatively characterize spatial distribution of pore type from seismic data on carbonate reservoirs by using only seismic rock-physics models. In step one, the pore space of carbonate reservoirs was equivalent to a single-porosity system with an effective pore aspect ratio; then, a 3D RPT on P- and S-wave impedances and density with the effective aspect ratio of pore,porosity and fluid saturation was established through the Gassmann’s equations and effective medium models; and then, the effective aspect ratio of pore, porosity and fluid saturation were simultaneously estimated from the seismic data based on 3D RPT. In step two, the pore space of carbonate reservoirs was equivalent to a triple-porosity system just as what the Xu–Payne model de fines. Combined with the inverted porosity and saturation in the first step,the porosities of three pore types (reference pores, stiff pores and cracks) were inverted from the seismic elastic properties following the Xu–Payne’s scheme.

    The effectiveness on porous rock taken as a single-porosity system was testi fied by an effective equation derived from the differential effective medium theory and numerical experiments on theoretical models. The predicted elastic moduli for both dry rocks and the water- and oil-saturated rocks by the differential effective medium (DEM) model of multiple-porosity rock are well matched with the simulated results by the DEM model of single-porosity rock.

    The applications to well-logging data and seismic data both indicate that the proposed inversion scheme can obtain more accurate physical properties and ensure the reliability of characterization of pore-type distribution. Meanwhile, the estimated pore type shows that the dolomite reservoir has strongly heterogeneous pore space in the target study area.

    AcknowledgementsThis work was supported by the China National Key R D plan (2019YFC0605504) and Scienti fic Research & Technology Development Project of China National Petroleum Corporation(Grant Nos. 2017D-3504 and 2018D-4305). We thank the editors and reviewers for many helpful comments and suggestions that improved this paper.

    Open AccessThis article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source,provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creat iveco mmons.org/licen ses/by/4.0/.

    精品卡一卡二卡四卡免费| 大片免费播放器 马上看| 国产永久视频网站| 免费观看在线日韩| 制服丝袜香蕉在线| 最近中文字幕2019免费版| 国产永久视频网站| 婷婷色麻豆天堂久久| 中文天堂在线官网| 久久久久精品久久久久真实原创| 三级经典国产精品| 久久久久久久久久人人人人人人| 中国美白少妇内射xxxbb| 看免费成人av毛片| 麻豆成人av视频| 精品少妇黑人巨大在线播放| 热re99久久精品国产66热6| 啦啦啦视频在线资源免费观看| 人人澡人人妻人| 人人澡人人妻人| 久久久久久久精品精品| 亚洲色图综合在线观看| 国精品久久久久久国模美| 色视频在线一区二区三区| 成人特级av手机在线观看| 女人精品久久久久毛片| 在线看a的网站| 国产成人精品久久久久久| 纵有疾风起免费观看全集完整版| 黑人猛操日本美女一级片| 久久久久网色| 赤兔流量卡办理| 纵有疾风起免费观看全集完整版| 国产精品欧美亚洲77777| 免费观看的影片在线观看| 十八禁高潮呻吟视频 | 欧美最新免费一区二区三区| 成人影院久久| 欧美bdsm另类| 91aial.com中文字幕在线观看| 久久久午夜欧美精品| 亚洲欧美精品专区久久| 777米奇影视久久| 亚洲婷婷狠狠爱综合网| 观看美女的网站| 91久久精品国产一区二区成人| 国产在视频线精品| 热99国产精品久久久久久7| 国产精品一区www在线观看| 国产欧美日韩一区二区三区在线 | 午夜老司机福利剧场| 曰老女人黄片| 高清毛片免费看| 一级爰片在线观看| 国产毛片在线视频| 哪个播放器可以免费观看大片| 国产精品女同一区二区软件| 99热这里只有精品一区| 欧美区成人在线视频| 极品人妻少妇av视频| 国产精品久久久久久av不卡| 如何舔出高潮| 成年人免费黄色播放视频 | 欧美xxxx性猛交bbbb| 成年女人在线观看亚洲视频| 交换朋友夫妻互换小说| 青春草亚洲视频在线观看| 国产日韩欧美视频二区| 国产欧美亚洲国产| 99久久中文字幕三级久久日本| av卡一久久| av线在线观看网站| 亚洲第一区二区三区不卡| 日韩中文字幕视频在线看片| 简卡轻食公司| 久久av网站| 街头女战士在线观看网站| 色视频在线一区二区三区| 极品教师在线视频| 少妇高潮的动态图| 黑人高潮一二区| 国产精品熟女久久久久浪| 黄色欧美视频在线观看| 免费高清在线观看视频在线观看| 成人18禁高潮啪啪吃奶动态图 | 久久精品国产a三级三级三级| 亚洲av二区三区四区| 一区在线观看完整版| 成人国产麻豆网| 熟妇人妻不卡中文字幕| 我要看黄色一级片免费的| 中文乱码字字幕精品一区二区三区| 丝袜喷水一区| 亚洲人成网站在线观看播放| 伊人亚洲综合成人网| 国产极品粉嫩免费观看在线 | 国产成人精品福利久久| 美女大奶头黄色视频| 女人精品久久久久毛片| 国产av国产精品国产| 亚洲精品成人av观看孕妇| 国产日韩欧美亚洲二区| 十八禁网站网址无遮挡 | 国产黄色视频一区二区在线观看| 久久久国产欧美日韩av| 国产精品.久久久| 国产精品国产三级专区第一集| 大码成人一级视频| 在线观看免费高清a一片| 麻豆乱淫一区二区| 久久狼人影院| 女性生殖器流出的白浆| 国产亚洲5aaaaa淫片| 精品午夜福利在线看| 美女福利国产在线| 国产精品蜜桃在线观看| 欧美97在线视频| 亚洲欧美成人精品一区二区| 免费观看在线日韩| 国产一区二区在线观看日韩| 精品久久国产蜜桃| 久久精品久久精品一区二区三区| 91精品国产国语对白视频| 欧美 亚洲 国产 日韩一| 秋霞在线观看毛片| 中文字幕人妻熟人妻熟丝袜美| 在线观看人妻少妇| 国产av国产精品国产| 亚洲综合精品二区| 国产精品一区www在线观看| 如何舔出高潮| 大又大粗又爽又黄少妇毛片口| 蜜桃在线观看..| 色94色欧美一区二区| 欧美日韩视频精品一区| 人妻人人澡人人爽人人| 欧美高清成人免费视频www| 亚洲美女黄色视频免费看| 黄色毛片三级朝国网站 | 啦啦啦在线观看免费高清www| 最近的中文字幕免费完整| 亚洲性久久影院| 99热国产这里只有精品6| 午夜老司机福利剧场| 久久狼人影院| 久久久久精品久久久久真实原创| 欧美精品亚洲一区二区| 亚洲中文av在线| 18禁裸乳无遮挡动漫免费视频| 久久国内精品自在自线图片| 国产深夜福利视频在线观看| 久久午夜综合久久蜜桃| 日韩大片免费观看网站| 搡老乐熟女国产| 成人无遮挡网站| 国产中年淑女户外野战色| 午夜精品国产一区二区电影| 久久狼人影院| 免费看光身美女| 亚洲一区二区三区欧美精品| 日本vs欧美在线观看视频 | 亚洲av在线观看美女高潮| 80岁老熟妇乱子伦牲交| 人体艺术视频欧美日本| 国产成人精品无人区| 一区二区av电影网| 日本爱情动作片www.在线观看| 嘟嘟电影网在线观看| 日本黄色日本黄色录像| 亚洲熟女精品中文字幕| 国内精品宾馆在线| 人妻制服诱惑在线中文字幕| 亚洲自偷自拍三级| 久热久热在线精品观看| 在线看a的网站| 少妇裸体淫交视频免费看高清| 久久精品国产亚洲av天美| 国产免费一级a男人的天堂| 偷拍熟女少妇极品色| 伊人久久国产一区二区| 国产精品无大码| videossex国产| 80岁老熟妇乱子伦牲交| 最近最新中文字幕免费大全7| 亚洲欧美一区二区三区国产| 国产欧美亚洲国产| 狂野欧美激情性xxxx在线观看| 熟妇人妻不卡中文字幕| 男人舔奶头视频| 91精品国产国语对白视频| 久久精品熟女亚洲av麻豆精品| 亚洲熟女精品中文字幕| 男人爽女人下面视频在线观看| 亚洲欧美成人精品一区二区| 又粗又硬又长又爽又黄的视频| 久久精品久久久久久久性| 插逼视频在线观看| 中文字幕av电影在线播放| 免费少妇av软件| 中文欧美无线码| 一区二区三区乱码不卡18| 少妇的逼水好多| 日韩不卡一区二区三区视频在线| 欧美丝袜亚洲另类| 91久久精品电影网| 日韩免费高清中文字幕av| 熟女av电影| 成人国产av品久久久| 中文字幕人妻熟人妻熟丝袜美| 亚洲精品乱码久久久久久按摩| 日日爽夜夜爽网站| 蜜桃在线观看..| 国产一区二区在线观看av| 最近2019中文字幕mv第一页| 国产精品99久久久久久久久| 观看av在线不卡| 亚洲精品一二三| 三上悠亚av全集在线观看 | 日日啪夜夜撸| 亚洲av福利一区| 插阴视频在线观看视频| 韩国高清视频一区二区三区| 丝袜在线中文字幕| 欧美97在线视频| 大话2 男鬼变身卡| 亚洲欧洲日产国产| 国产乱来视频区| 亚洲无线观看免费| 3wmmmm亚洲av在线观看| 99九九线精品视频在线观看视频| 麻豆乱淫一区二区| 久久av网站| 日韩中字成人| 日本爱情动作片www.在线观看| 欧美xxxx性猛交bbbb| 亚洲精品国产av蜜桃| 午夜视频国产福利| 在线看a的网站| 18禁裸乳无遮挡动漫免费视频| 性色avwww在线观看| 中文字幕免费在线视频6| 国产在线一区二区三区精| 国产精品一区二区在线观看99| freevideosex欧美| 国产国拍精品亚洲av在线观看| 午夜老司机福利剧场| 内地一区二区视频在线| 亚洲人成网站在线观看播放| 边亲边吃奶的免费视频| 久久久久视频综合| av国产久精品久网站免费入址| 久久精品国产鲁丝片午夜精品| 午夜福利网站1000一区二区三区| 丰满少妇做爰视频| 女人久久www免费人成看片| 国产精品国产三级专区第一集| 国精品久久久久久国模美| 少妇被粗大的猛进出69影院 | 伦理电影大哥的女人| 91久久精品国产一区二区三区| 97在线人人人人妻| 久久精品久久久久久噜噜老黄| 亚洲欧美一区二区三区国产| 国产女主播在线喷水免费视频网站| 91精品伊人久久大香线蕉| 蜜桃在线观看..| 国产男人的电影天堂91| videos熟女内射| 欧美日韩精品成人综合77777| 最近2019中文字幕mv第一页| 久久99一区二区三区| 国产精品久久久久久av不卡| 黄色一级大片看看| 视频中文字幕在线观看| 亚洲精品一二三| 国产亚洲91精品色在线| 新久久久久国产一级毛片| 免费观看a级毛片全部| 欧美精品亚洲一区二区| 成人毛片a级毛片在线播放| 99九九在线精品视频 | 欧美 亚洲 国产 日韩一| 女人久久www免费人成看片| 狂野欧美激情性xxxx在线观看| 日本黄色片子视频| 亚洲av不卡在线观看| 久久久久久久精品精品| 精品人妻熟女av久视频| 成人午夜精彩视频在线观看| 欧美xxⅹ黑人| 国产亚洲欧美精品永久| 毛片一级片免费看久久久久| 99热这里只有是精品在线观看| 尾随美女入室| 一级,二级,三级黄色视频| 成人漫画全彩无遮挡| a级毛色黄片| 99国产精品免费福利视频| 汤姆久久久久久久影院中文字幕| 我要看日韩黄色一级片| 国产精品一二三区在线看| 99re6热这里在线精品视频| 久久99蜜桃精品久久| 最近中文字幕高清免费大全6| 亚洲国产精品国产精品| 国产欧美日韩综合在线一区二区 | 国产在线视频一区二区| 国产综合精华液| 午夜视频国产福利| 精品一区在线观看国产| 久久这里有精品视频免费| 日韩中文字幕视频在线看片| 亚洲av成人精品一区久久| 美女国产视频在线观看| 中文字幕人妻熟人妻熟丝袜美| 最新中文字幕久久久久| 18禁动态无遮挡网站| 国产真实伦视频高清在线观看| 国产中年淑女户外野战色| 欧美精品人与动牲交sv欧美| 亚洲国产精品成人久久小说| 欧美日韩一区二区视频在线观看视频在线| 夜夜骑夜夜射夜夜干| 老司机影院毛片| 在线观看免费高清a一片| 欧美日本中文国产一区发布| 一区二区三区四区激情视频| av有码第一页| 美女内射精品一级片tv| 大香蕉久久网| 99久久综合免费| 国产欧美亚洲国产| 国产黄片视频在线免费观看| 日本黄色日本黄色录像| 国产成人免费无遮挡视频| 在线精品无人区一区二区三| 青青草视频在线视频观看| 蜜桃在线观看..| 精品一区在线观看国产| 国产精品人妻久久久久久| 国产精品免费大片| 欧美变态另类bdsm刘玥| 免费久久久久久久精品成人欧美视频 | 日本-黄色视频高清免费观看| 99九九在线精品视频 | 乱系列少妇在线播放| 午夜福利在线观看免费完整高清在| 99九九在线精品视频 | 最近的中文字幕免费完整| 交换朋友夫妻互换小说| 97精品久久久久久久久久精品| 午夜福利网站1000一区二区三区| 精品人妻熟女av久视频| 日日啪夜夜撸| 国产无遮挡羞羞视频在线观看| 国产黄频视频在线观看| 亚洲国产欧美日韩在线播放 | 久久久久久久亚洲中文字幕| 国产精品人妻久久久影院| 国产高清三级在线| 2018国产大陆天天弄谢| 欧美变态另类bdsm刘玥| 偷拍熟女少妇极品色| 国产毛片在线视频| 黄色怎么调成土黄色| 在线精品无人区一区二区三| 久久6这里有精品| 看免费成人av毛片| 好男人视频免费观看在线| 夫妻性生交免费视频一级片| 一区二区三区乱码不卡18| 国产熟女欧美一区二区| 国产在视频线精品| 制服丝袜香蕉在线| 狂野欧美激情性bbbbbb| 日本与韩国留学比较| 青春草国产在线视频| 黄片无遮挡物在线观看| 亚洲av不卡在线观看| 欧美bdsm另类| 国产日韩欧美亚洲二区| 国产亚洲91精品色在线| 亚洲丝袜综合中文字幕| 亚州av有码| 丝袜在线中文字幕| av有码第一页| 午夜福利,免费看| 菩萨蛮人人尽说江南好唐韦庄| 久久久久网色| 亚洲性久久影院| 2018国产大陆天天弄谢| 久久久久国产精品人妻一区二区| 国产成人精品福利久久| 久久久久精品久久久久真实原创| 亚洲国产最新在线播放| 成人综合一区亚洲| 99热国产这里只有精品6| 黄色视频在线播放观看不卡| 久久精品国产亚洲av涩爱| 我的女老师完整版在线观看| 亚洲四区av| 免费大片黄手机在线观看| 97在线视频观看| 黑人巨大精品欧美一区二区蜜桃 | 午夜激情福利司机影院| 大陆偷拍与自拍| 国产综合精华液| 成人午夜精彩视频在线观看| 亚洲精品色激情综合| 亚洲国产毛片av蜜桃av| 成人特级av手机在线观看| 中文资源天堂在线| 国产黄片视频在线免费观看| 黑丝袜美女国产一区| 国产亚洲精品久久久com| 日韩欧美 国产精品| 久久久久国产精品人妻一区二区| 亚洲熟女精品中文字幕| 日日啪夜夜爽| 自拍欧美九色日韩亚洲蝌蚪91 | 色94色欧美一区二区| 精品卡一卡二卡四卡免费| 成人亚洲精品一区在线观看| 91在线精品国自产拍蜜月| 国产精品国产av在线观看| 国产精品女同一区二区软件| 三级经典国产精品| 国产精品嫩草影院av在线观看| 日韩av不卡免费在线播放| 国产精品99久久99久久久不卡 | 亚洲色图综合在线观看| 欧美区成人在线视频| 国产老妇伦熟女老妇高清| 精品人妻偷拍中文字幕| 成人二区视频| 亚洲国产av新网站| 精品视频人人做人人爽| 中国国产av一级| 日本黄色日本黄色录像| 日韩三级伦理在线观看| 蜜臀久久99精品久久宅男| 2021少妇久久久久久久久久久| 99re6热这里在线精品视频| 夫妻午夜视频| 久久国产精品男人的天堂亚洲 | 国产乱人偷精品视频| 亚洲av国产av综合av卡| 国产淫片久久久久久久久| 国产av码专区亚洲av| 啦啦啦视频在线资源免费观看| 亚洲欧美成人精品一区二区| 亚洲内射少妇av| 久久久久久久亚洲中文字幕| 国产一区二区在线观看av| 色视频www国产| 2022亚洲国产成人精品| 国产美女午夜福利| 亚洲av日韩在线播放| 免费黄色在线免费观看| 少妇丰满av| 91成人精品电影| 多毛熟女@视频| 国产探花极品一区二区| 日韩免费高清中文字幕av| 天天躁夜夜躁狠狠久久av| 丰满饥渴人妻一区二区三| 久久精品国产自在天天线| av在线app专区| 色哟哟·www| 亚洲无线观看免费| av国产久精品久网站免费入址| 国产女主播在线喷水免费视频网站| 国产视频内射| 精品国产一区二区久久| 日韩电影二区| 啦啦啦中文免费视频观看日本| 亚洲成人av在线免费| 精品一区二区三区视频在线| 街头女战士在线观看网站| 久久久久久久久久久丰满| 自拍欧美九色日韩亚洲蝌蚪91 | 能在线免费看毛片的网站| av网站免费在线观看视频| 欧美丝袜亚洲另类| 女人精品久久久久毛片| 亚洲精华国产精华液的使用体验| av卡一久久| 永久免费av网站大全| 国产无遮挡羞羞视频在线观看| www.色视频.com| 日韩欧美一区视频在线观看 | 亚洲第一区二区三区不卡| 午夜激情福利司机影院| 日本wwww免费看| 久久精品熟女亚洲av麻豆精品| 9色porny在线观看| 高清午夜精品一区二区三区| 99热6这里只有精品| 久久精品久久久久久噜噜老黄| 免费少妇av软件| 在线观看美女被高潮喷水网站| 人人妻人人澡人人爽人人夜夜| 简卡轻食公司| 欧美日韩精品成人综合77777| 国产高清三级在线| 在线观看美女被高潮喷水网站| 国产av码专区亚洲av| 少妇熟女欧美另类| 美女国产视频在线观看| 在线播放无遮挡| 欧美日韩国产mv在线观看视频| 美女大奶头黄色视频| 中文字幕亚洲精品专区| 国产男女内射视频| 成年人免费黄色播放视频 | 亚洲精品一区蜜桃| 街头女战士在线观看网站| 桃花免费在线播放| 欧美3d第一页| 亚洲国产精品一区二区三区在线| 夜夜骑夜夜射夜夜干| 国产精品久久久久久久电影| 国产精品熟女久久久久浪| 国内少妇人妻偷人精品xxx网站| 亚洲精品第二区| 国产一区二区三区av在线| 91精品国产九色| 久久99精品国语久久久| 色视频www国产| 欧美日韩精品成人综合77777| 欧美成人午夜免费资源| 日日撸夜夜添| 精品人妻一区二区三区麻豆| 精品久久久久久久久亚洲| 女的被弄到高潮叫床怎么办| 免费高清在线观看视频在线观看| 免费在线观看成人毛片| 美女大奶头黄色视频| 国产亚洲精品久久久com| 成人亚洲欧美一区二区av| 18禁在线播放成人免费| 亚洲,一卡二卡三卡| 伦精品一区二区三区| 亚洲精品国产成人久久av| 欧美高清成人免费视频www| 欧美另类一区| 久久精品熟女亚洲av麻豆精品| a级毛片免费高清观看在线播放| kizo精华| 久久青草综合色| 久久女婷五月综合色啪小说| 夜夜骑夜夜射夜夜干| 国产免费一级a男人的天堂| 97精品久久久久久久久久精品| 国产亚洲5aaaaa淫片| 91精品国产九色| 亚洲怡红院男人天堂| 亚洲av成人精品一区久久| 国产亚洲午夜精品一区二区久久| 人人妻人人看人人澡| 国产亚洲最大av| 欧美人与善性xxx| 午夜福利视频精品| 亚洲电影在线观看av| 精品人妻熟女av久视频| 在线观看免费日韩欧美大片 | 水蜜桃什么品种好| 国产精品福利在线免费观看| 亚洲精品456在线播放app| 99久久精品热视频| 精品久久久噜噜| 精品久久久久久久久av| 国产成人午夜福利电影在线观看| 亚洲va在线va天堂va国产| 午夜av观看不卡| 免费人成在线观看视频色| 91精品国产国语对白视频| 婷婷色麻豆天堂久久| 少妇精品久久久久久久| 日日爽夜夜爽网站| 久久久午夜欧美精品| 国产亚洲av片在线观看秒播厂| 黄色怎么调成土黄色| 在线观看av片永久免费下载| 高清在线视频一区二区三区| 综合色丁香网| 亚洲av欧美aⅴ国产| 精品人妻一区二区三区麻豆| 人人妻人人添人人爽欧美一区卜| 男的添女的下面高潮视频| 又爽又黄a免费视频| 校园人妻丝袜中文字幕| 日韩人妻高清精品专区| 国产69精品久久久久777片| 日韩精品有码人妻一区| 一级毛片我不卡| 亚洲av成人精品一区久久| 国产国拍精品亚洲av在线观看| 精品一品国产午夜福利视频| 国产日韩一区二区三区精品不卡 | 亚洲精品中文字幕在线视频 | 精品人妻熟女av久视频| 日韩伦理黄色片| 18禁动态无遮挡网站| 国产一区二区在线观看av| 男人添女人高潮全过程视频| 伦理电影免费视频| 国产在线免费精品| 精品国产国语对白av| 免费观看在线日韩| 最近2019中文字幕mv第一页| 国产精品无大码| 高清在线视频一区二区三区| 精品久久久久久久久av| 内射极品少妇av片p|