999精品在线视频,手机成人午夜在线视频,久久不卡国产精品无码,中日无码在线观看,成人av手机在线观看,日韩精品亚洲一区中文字幕,亚洲av无码人妻,四虎国产在线观看 ?

Application of machine learning in predicting the rate-dependent compressive strength of rocks

2022-10-09 12:51:42MingdongWeiWenzhoMengFengDiWeiWu

Mingdong Wei,Wenzho Meng,Feng Di,Wei Wu,*

a School of Civil and Environmental Engineering,Nanyang Technological University,Singapore

b State Key Laboratory of Hydraulics and Mountain River Engineering,College of Water Resource and Hydropower,Sichuan University,Chengdu,610065,China

Keywords:Machine learning Rock dynamics Compressive strength Strain rate

A B S T R A C T Accurate prediction of compressive strength of rocks relies on the rate-dependent behaviors of rocks,and correlation among the geometrical,physical,and mechanical properties of rocks.However,these properties may not be easy to control in laboratory experiments,particularly in dynamic compression experiments.By training three machine learning models based on the support vector machine(SVM),backpropagation neural network(BPNN),and random forest(RF)algorithms,we isolated different input parameters,such as static compressive strength,P-wave velocity,specimen dimension,grain size,bulk density,and strain rate,to identify their importance in the strength prediction.Our results demonstrated that the RF algorithm shows a better performance than the other two algorithms.The strain rate is a key input parameter influencing the performance of these models,while the others(e.g.static compressive strength and P-wave velocity)are less important as their roles can be compensated by alternative parameters.The results also revealed that the effect of specimen dimension on the rock strength can be overshadowed at high strain rates,while the effect on the dynamic increase factor(i.e.the ratio of dynamic to static compressive strength)becomes significant.The dynamic increase factors for different specimen dimensions bifurcate when the strain rate reaches a relatively high value,a clue to improve our understanding of the transitional behaviors of rocks from low to high strain rates.

1.Introduction

Rate dependence of rock strength is a critical issue in rock dynamics,as rock behaviors vary drastically at different strain rates.Understanding the rate-dependent strengths of rocks is important for rock engineering design and construction(Zhao et al.,1999;Hudson,2012).Our current approaches in determining the compressive strength of rocks rely basically on laboratory experiments,such as the static compression test(ISRM,1979)and the split Hopkinson pressure bar(SHPB)test(Zhou et al.,2012).These test methods investigate the rate dependence of compressive strength of rocks at low and high strain rates,respectively.A general understanding is that the compressive strength slightly increases with a higher strain rate at a low strain rate(e.g.<5 s-1)and grows rapidly at a high strain rate(e.g.>500 s-1).The test methods are difficult to be used directly to investigate the compressive strength of rocks at intermediate strain rates(5-500 s-1),which are of interest to understand the mechanisms of excavation-induced geohazards such as rockburst(Duan et al.,2019a;Zhong et al.,2021).Efforts have been made to modify the test setups to control the intermediate strain rates(e.g.Li et al.,2005;Whittington et al.,2015).However,a few issues remain challenging,such as inertia effects from different test methods,specimen dimensions used in different test methods,inconstant strain rates in different individual tests.Numerical simulations can control these effects and have been used extensively(e.g.Li et al.,2009;Du et al.,2018;Huang et al.,2018;Duan et al.,2019b),but an accurate prediction of rock strength depends on a reasonable constitutive model and a reliable model calibration.Therefore,a simple,fast,and data-driven method is needed to predict the compressive strength of rocks at various strain rates and to promote our understanding of the rate-dependent strength of rocks.

In recent years,with the development of computer and artificial intelligence technology,machine learning models have been widely employed to address the complexity of engineering issues,including the nonlinear behaviors of rocks.For example,the support vector machine(SVM)algorithm was applied to forecast the static compressive strength and shear strength of rocks by inputting rock density,point load strength,P-wave velocity,and slake durability index(Li and Tan,2017).The back-propagation neural network(BPNN)algorithm was used to predict the static compressive strength and the elasticity modulus of rocks based on slake durability index,Schmidt hammer rebound number,effective porosity,and point load index(Y?lmaz and Yuksek,2008).Based on the random forest(RF)algorithm,the static compressive strength of rocks was estimated based on the available drilling data,such as weight on bit,drill string rotating speed,drilling torque,stand-pipe pressure,mud pumping rate,and penetration rate(Gamal et al.,2021).Given these successful applications to rock mechanics problems,the application of machine learning models can further be extended to rock dynamics to improve our ability to forecast the rate-dependent strength of rocks.

This study aimed to use machine learning models to predict the compressive strength of rocks and to provide new insights into the rate-dependent strength of rocks.We first collected the laboratory data of the static compression tests and the SHPB tests on various types of hard rocks as well as the corresponding geometrical,physical,and mechanical properties.We then introduced the machine learning algorithms(e.g.SVM,BPNN,and RF)and established three models based on these algorithms.We then trained these models using three quarters of the data and validated the models using the remaining data.Finally the importance of input parameters and the effects of strain rate and specimen dimension on the compressive strength of rocks were discussed.

2.Data and methods

2.1.Data collection

To train and validate the machine learning models in predicting the compressive strength of rocks,we first collected a few data sets of the static compression tests and the SHPB tests on eleven types of hard rocks from published papers(e.g.Frew et al.,2001;Li et al.,2005;Xia et al.,2008;Wang and Tonon,2011;Zhang and Zhao,2013;Zhou et al.,2017,2020;Weng et al.,2020;Yu et al.,2021).Typical rock types in rock engineering were included,such as granite,marble,sandstone,and siltstone.The geometrical,physical,and mechanical properties of these rocks are listed in Table A1,including specimen length and diameter,grain size,bulk density,Pwave velocity,as well as compressive strength with the corresponding strain rate.The static compressive strength was considered as rate-independent and approximated as a constant value,not specific for each specimen.Selection of the input parameters mainly relies on data availability,and other parameters,such as porosity,mineralogy,and microscale features(e.g.Fan et al.,2018,2020),can also be used if available in the literature.

To understand the simple correlation between any two parameters in Table A1,we calculated the Pearson correlation coefficient(Edwards,1976),a ratio of the covariance of two parameters to the product of their standard deviations.A correlation coefficient of 1 indicates a strong positive correlation,while-1 means a strong negative correlation.As illustrated in Fig.1,for example,the P-wave velocity shows the strongest correlation with the static compressive strength,as they are highly related to the integrity of rocks.The P-wave velocity is also positively correlated with the grain size and the bulk density.The strain rate is poorly correlated with the other parameters,as it is controlled by the loading condition.The specimen length and diameter are strongly correlated with each other but poorly correlated with the other parameters due to a length-todiameter ratio of 0.5-1 suggested by the SHPB test standard(Zhou et al.,2012)and widely adopted in many studies.

Fig.1.Correlations between input parameters based on Pearson correlation coefficients.Values 1 and-1 mean strong positive and negative correlations,respectively.

2.2.Machine learning algorithms

Three well-proved machine learning algorithms,i.e.SVM,BPNN,and RF,in rock mechanics were used.However,their applicability to rock dynamics has rarely been reported yet.Fig.2 illustrates the specific steps of applying the three algorithms to predict the ratedependent compressive strength of rocks.The basic features of these algorithms are summarized in Figs.3-5.

2.2.1.SVM algorithm

The SVM algorithm was originally proposed by Boser et al.(1992)for data classification and subsequently developed by Drucker et al.(1996)for regression and prediction analysis.The SVM based on the principle of structural risk minimization can minimize empirical risk and confidence range,as well as obtain excellent statistical laws even for small sample sizes(Samui,2008).Fig.3 shows a schematic diagram of the basic principle of the SVM to solve the regression and prediction problems.The goal of the SVM training is to find the best hyperplane closest to all data points.The training process is to continuously adjust the hyperplane parameters to minimize the sum of the distances of all data points from the hyperplane plane.In this study,the SVM is used to construct the mapping relationship between the compressive strength and the other parameters in Table A1.Thus,the feature vector of the compressive strength can be written as{xi,yi}(i=1,2,…,m;xi∈Rd;yi∈R),where x represents the input,ydenotes the output,mis the number of data points,RandRdare the onedimensional andd-dimensional vector space,respectively,anddis the number of input variables.The value ofy(i.e.compressive strength)can be solved as

where w is a vector of weight coefficient,φ(x)is a nonlinear mapping function,andbis a bias constant.

The parameters w andbcan be determined by minimizingunder the following premise(Samui 2008):

Fig.2.A flowchart summarizing how to apply SVM,BPNN,and RF algorithms to establish machine learning models.

Fig.3.A schematic diagram of the principle of SVM algorithm.

Fig.4.A schematic diagram of the principle of BPNN algorithm.

where ε is a free parameter considered as a threshold.By solving Eq.(4),the function model for support vector regression can be determined as

whereK(xi,x)is the kernel function,adopting the radial basis function to address nonlinear problems(Samui,2008):

As summarized in Fig.2,the two parameters(i.e.cand γ)in the radial basis function should be optimized.The SVM-based model aims to find the optimumcand γ values through the grid search method(Cheung et al.,1997)and to achieve good generalization capacity and high calculation speed.

2.2.2.BPNN algorithm

The BPNN algorithm developed by Rumelhart et al.(1986)is a multi-layer feed-forward algorithm.The process of constructing the BPNN includes signal forward propagation and error backpropagation.The structure of the BPNN consists of an input layer,one or more hidden layers,and an output layer(Fig.4).In this study,the compressive strength is set as the output neuron,and the others(e.g.strain rate,P-wave velocity)are the inputs.To accelerate the convergence of the neural network,the raw data(X)are normalized using Eq.(7)and mapped into a range of 0-1:

whereXiis the normalized data;andXmaxandXminare the maximum and minimum values,respectively.

The number of hidden-layer nodes has a significant influence on the BPNN training and is determined as

wherehandkare the number of nodes in the hidden and output layers,respectively,andais a constant with a value of 0-10.In our model,d=7 andk=1.By takinga=8,hcan be determined as 10.Information transfer between the layers relies on the matrix operation:

In the error back-propagation algorithm,the Sigmoid function is used to introduce nonlinearity to the BPNN-based model to better solve complex nonlinear problems:

In the initial BPNN-based model,the weights and biases are randomly initialized.After being fed with the input data,the model compares the output data using the loss function:

whereflossis the prediction error,and^yandyare the measured and the predicted strength values,respectively.Because of the random parameter setting,a learning mechanism should be implemented to optimize the model.

The BPNN-based model has the two parameters(i.e.to be adjusted during the training process(Fig.2).Based on the prediction error,the error term of each hidden layer is calculated in reverse order through the error back-propagation algorithm,and the gradient descent method is applied to adjust the two parameters of each hidden layer(Hochreiter et al.,2001).The adjustment of the weight and bias values is determined by the gradient of each parameter calculated using the loss function.The dropout regularization method is used to prevent the model from overfitting(Srivastava et al.,2014).An unacceptable error is propagated backwards again,and the weight and bias values are adjusted accordingly.This process is repeated until finding an acceptable error or reaching the maximum cycle number.

2.2.3.RF algorithm

The RF algorithm proposed by Breiman(2001)can integrate multiple decision trees into an ensemble algorithm and average the results of these decision trees as the output(Fig.5).The RF has a fast convergence speed and can prevent data from overfitting to improve the prediction accuracy.The main features of the RF include(i)bootstrap sampling,(ii)random feature selection,(iii)out-of-bag(OOB)error estimation,and(iv)full-depth decision tree growth(Li et al.,2018).The training data set(Sm)containingmobservations can be written as

Fig.5.A schematic diagram of the principle of RF algorithm.

wheremis the number of training data points.

Also shown in Fig.2,the bootstrap sampling is performed to generatetsets of training data from the original data set.Each new data set corresponds to a decision tree that is fully grown and not pruned by the classification and regression trees(CART)methodology(Duda et al.,2000).The RF can find the best split among the randomly selected features at each node.It is common to artificially control the number of branches in the decision tree or let it grow freely for a limited number of times and subsequently to seek the best single decision tree.After the training process,a global prediction functiony=h(x,Sm)is constructed over Sm.The RF produces the prediction result by averaging the output of these decision trees(i.eTherefore,the final output obtained from a new input vector x is:

An optimal RF-based model is realized by tuning the number of decision trees and the depth of each decision tree.The model can be more robust and accurate with a larger number of decision trees.However,the number is also associated with the calculation time,and the prediction accuracy is prone to stabilizing after hitting a critical point.The depth of the decision tree affects its own strength and the correlation between the trees.The purpose of the RF is to find a tradeoff among the calculation time,decision tree strength,and correlation between the trees.When the terminal conditions are reached,the model is considered acceptable,and the training is completed.

3.Model training and validation

We trained three machine learning models based on the SVM,the BPNN,and the RF algorithms by randomly selecting three quarters of the data in Table A1 and checked the performance of these models based on the remaining data.The compressive strength was the output,and the other parameters in Table A1 were used as the inputs(i.e.static compressive strength,P-wave velocity,specimen dimension,grain size,bulk density,and strain rate).

We first qualitatively evaluated the performance of the three models based on the comparison between the predicted and the measured compressive strengths(Fig.6).A data point on or close to the diagonal line indicates that the predicted compressive strength matches with the measured compressive strength.The data points with a small deviation from the diagonal line are mainly distributed in the region less than 120 MPa compressive strength.This is because the data points are densely distributed in this region,as shown in the histograms next to the upper and right axes.The results indicate that the amount of data points is critical for the accuracy of machine learning prediction.

We then quantitatively compared the performance of these models based on the following metrics:

Fig.6.Comparison between the measured and predicted strengths from(a)SVMbased,(b)BPNN-based,and(c)RF-based models.The histograms on the upper and right axes show the density distributions of strength values.

whereMAE,MXAE,MAPE,andRMSEdenote the mean absolute error,the maximum absolute error,the mean absolute percentage error,and the root mean square error,respectively;nrepresents the number of data points;andyiand^yiare the predicted and the measured values for theith strength value,respectively.

The performance metrics plotted in Fig.7 demonstrate that the compressive strengths predicted by the three machine learning models are mostly similar.The RF-based performance metrics are lower than the SVM-and BPNN-based performance metrics,manifesting the superiority of the RF over the SVM and the BPNN in predicting the compressive strength of rocks.The BPNN-based performance metrics are lower than the SVM-based performance metrics,except the MAPE.The high MAPE value is attributed to the significant absolute error obtained from the SVM-based model when predicting a high strength value.Given the better performance of the RF-based model,we used this model to investigate the rate dependence of rock strength in Section 4.

4.Discussion

4.1.Sensitivity analysis

Our study demonstrates that the machine learning models perform well to predict the compressive strength of rocks using the six input parameters(i.e.static compressive strength,P-wave velocity,specimen dimension,grain size,bulk density,and strain rate).In rock engineering,one may not be able to collect sufficient data points to train the three models and to ensure the accuracy of machine learning prediction.Insufficient data points of single input parameter may influence significantly or negligibly the performance of these models.Here,we performed a sensitivity analysis to evaluate the importance of each input parameter in the strength prediction.In each model,we excluded one of the six input parameters in turn.In other words,we inputted the other five parameters in the model and compared the performance difference between the model with five input parameters and the model with six input parameters.The procedure for the machine learning prediction of compressive strength followed that in Fig.2.We quantitatively assessed the sensitivity of the strength prediction to each excluded parameter based on the percentage increments of four performance metrics(P1j,P2j,P3j,andP4j):

Fig.7.Percentage metrics,including(a)MAE,(b)MXAE,(c)MAPE,and(d)RMSE,to evaluate the performance of the SVM-based,BPNN-based,and RF-based models,respectively.

where the subscriptsjand 0 refer to the performance metrics obtained from the models with five and six input parameters,respectively;andj(1,2,…and6)corresponds to staticcompressive strength,P-wave velocity,specimen dimension,grain size,bulk density,and strain rate,respectively.For example,P26denotes the percentage incrementof the MXAEvaluewhen the strain rate isexcludedinthese models.IfP26>P25,the strain rate plays a more critical role in predicting the compressive strength than the bulk density.

The percentage increments of the performance metrics plotted in Fig.8 exhibit that the strength prediction result is more sensitive to the strain rate than the other input parameters.When the strain rate is excluded from the input parameters of the machine learning models,the performance metrics except the MXAE show larger percentage increments,compared with those excluding the other five parameters.The exception of the MXAE is most likely attributed to one or a few data points with unusual values.The strain rate controlled by the loading condition is less dependent on the geometrical,physical,and mechanical properties of rocks and should be included in the strength prediction.This observation is similar to that observed from the Pearson correlation coefficient.The static compressive strength is correlated with P-wave velocity,grain size,and bulk density(Fig.1).If the static compressive strength is excluded,the impact of the static compressive strength on the strength prediction can be partially reflected by P-wave velocity and bulk density.That is probably why the performance metricsP11,P21,P31,andP41approach zero.The results indicate that the influence of specimen dimension is relatively minor,which may be addressed insufficiently,particularly the specimen diameter,as the bar diameter in the SHPB setup is mostly limited to a range of 25-100 mm.Also,the aspect ratio as a related factor is not well studied due to the requirement of dynamic stress equilibrium in rock specimens with appropriate aspect ratios in the SHPB tests.

4.2.Rate dependence of rock strength

Our results demonstrated that the RF-based model shows a better performance than the SVM-and the BPNN-based models(Fig.7),and the exclusion of the static compressive strength has a negligible impact on the prediction accuracy(Fig.8).Thus the RFbased model with five input parameters was adopted,except the static compressive strength,to further discuss the rate dependence of rock strength following the procedure in Fig.2.This model considers a series of fine-grained sandstone specimens(Yu et al.,2021)with the same aspect ratio and three different specimen diameters(12.5 mm,25 mm,and 50 mm).The prediction result shows that the compressive strength is relatively constant under the strain rate lower than 20 s-1and drastically increases under the strain rate higher than 20 s-1(Fig.9a).The increase in compressive strength is mainly attributed to the transition from inter-crystalline fracture to intra-crystalline fracture,the amplified energy absorption,and the inertial effect(Zhang et al.,2000).The trend of strength increase is similar for different specimen diameters,and the compressive strength decreases with increasing specimen diameter.Then,we used the dynamic increase factor(DIF),defined as the ratio of the dynamic to static strength,to highlight the effect of specimen diameter on the compressive strength.The three curves start to fluctuate when the strain rate exceeds 20 s-1(Fig.9b).From 20 s-1to 50 s-1,the DIF values of these specimens with the three diameters increase in a similar trend.However,when the strain rate exceeds 50 s-1,the DIF value of the specimen with 50 mm diameter increases faster than those of the other specimens.The effect of specimen diameter can be ascribed to an amplified inertia effect(Chen and Song,2011).The results additionally reveal that the compressive strength and the DIF value exhibit different responses to an increase in strain rate.The difference among the three curves for compressive strength decreases with a higher strain rate,while the three curves for the DIF values bifurcate when the strain rate is greater than 50 s-1.This is because the static compressive strength decreases with a larger specimen diameter(Masoumi et al.,2015).

Fig.8.Percentage increments of performance metrics:(a)MAE,(b)MXAE,(c)MAPE,and(d)RMSE,obtained from excluding each input parameters in the SVM-(red solid line),BPNN-(black dotted line),and RF-based(blue dashed line)models,respectively.

The RF-based model also considers similar fine-grained sandstone specimens with a fixed diameter(50 mm)and three different aspect ratios(0.5,0.75,and 1).The results exhibit that the compressive strength slightly decreases with increasing aspect ratio at a given strain rate,probably due to more defects in a longer specimen and weaker end frictioneffect(Fig.10a).The negative correlation between the compressive strength and specimen length is also observed in other rocks and solid materials(like concrete)(Krauthammer et al.,2003;Munoz and Taheri,2017).Fig.10b shows that the aspect ratio has a minor effect on the DIF value for a strain rate lower than 50 s-1,but obviously influences the DIF value beyond the strain rate.The fluctuation of thecurvesislikely due tothediscrepancies of input data points,which can be related to the system errors from the test setup and the data acquisition unit.Combining Figs.9b and 10b,it suggests that the influence of specimen dimension on the strength enhancement can beoverlookedindetermination of thecompressive strength at low strain rates(e.g.<50 s-1),but should be examined at intermediate and high strain rates.Both the figures also illustrate that the DIF values for different specimen dimensions increase at the same increment rate at low strain rates.However,the increment rate becomes different at higher strain rates.This observation provides new insights into the transitional behaviors of rocks from the low to high strain rates.The strength enhancement at the beginning of the transitional process is less dependent of the specimen dimension,and the dependence of the specimen dimension becomes obvious with increasing strain rate.

Fig.9.Effect of strain rate on(a)Compressive strength,and(b)Dynamic increase factor for different specimen diameters(12.5mm,25mm,and 50 mm),predicted by the RF-based model.

Our study highlights that the correlations among various geometrical,physical,and mechanical properties of rocks are nonlinear and very different from the simple correlations based on the Pearson correlation coefficient(Fig.1).Additionally,the nonlinear correlations can be even more complex if the three models can be trained by the data points obtained from laboratory and field experiments on rocks with larger sizes under higher strain rates.In other words,this study is limited with the data points in small ranges of specimen dimensions and strain rates available in the literature.The established models are applicable to forecast the compressive strength of hard rocks at a strain rate ranging from 5×10-6s-1to 223 s-1(Table A1).For other strain rates,additional training with the corresponding data is needed to extend the applicability of these models.The study demonstrates that machine learning models could be an effective method to explore the complexity of rock dynamics,particularly those difficult to control in laboratory experiments.

5.Conclusions

Fig.10.Effect of strain rate on(a)Compressive strength,and(b)Dynamic increase factor for different aspect ratios(0.5,0.75,and 1),predicted by the RF-based model.

This study addressed the feasibility of applying machine learning models to predict the rate-dependent strength of rocks.We first collected the strength data of various types of hard rocks as well as the corresponding geometrical,physical,and mechanical properties in the literature.We then established three machine learning models based on SVM,the BPNN,and the RF algorithms.In the machine learning models,the compressive strength was the output,and the other properties(i.e.static compressive strength,Pwave velocity,specimen dimension,grain size,bulk density,and strain rate)were the inputs.After training and validating the models,we discussed the importance of each input parameter in the strength prediction and the transitional behaviors of rocks with different specimen dimensions.

Our results indicated that the machine learning models can achieve reliable prediction for the rate-dependent strength of rocks.The RF algorithm manifests slight superiority to SVM and BPNN in predicting the rock strength.The strain rate as input is more important than the other parameters to ensure the excellent performance of the models.The dynamic increase factor is independent of specimen dimension at low strain rates but bifurcate for the strain rate over a relatively high level.This observation can be used to study the transitional behaviors of rocks from low to high strain rates,which may yield more insights into the specimen dimension effect on rock strength and further highlight the strain rate effect.This study introduced a new avenue to study the complex mechanical behaviors of rocks at various strain rates.

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgments

This research is supported by National Research Foundation,Singapore under its Virtual Singapore R&D Programme(Award No.NRF2019VSG-GMS-001).

Appendix A.Supplementary data

Supplementary data to this article can be found online at https://doi.org/10.1016/j.jrmge.2022.01.008.

主站蜘蛛池模板: 国产精品2| 精品视频福利| 国产精品永久免费嫩草研究院| 成人免费一区二区三区| 污污网站在线观看| 草逼视频国产| 中文国产成人久久精品小说| 欧美一级片在线| 2021精品国产自在现线看| 欧美翘臀一区二区三区| 色婷婷啪啪| 日韩欧美国产成人| 大陆国产精品视频| 91精品最新国内在线播放| 2018日日摸夜夜添狠狠躁| 五月天在线网站| 国产成人1024精品| 久久网欧美| 欧美日韩国产综合视频在线观看 | 香蕉蕉亚亚洲aav综合| 亚洲九九视频| 亚洲综合极品香蕉久久网| 国产精品一区在线麻豆| 狠狠综合久久久久综| 亚洲国产精品人久久电影| 国产成人精品亚洲日本对白优播| 在线中文字幕网| 久久人体视频| 国产资源站| 日韩毛片在线播放| 人妻丰满熟妇αv无码| 欧美日韩一区二区在线播放| a级毛片免费网站| 欧美黑人欧美精品刺激| 欧美精品啪啪一区二区三区| 亚洲日韩国产精品综合在线观看| 中文字幕伦视频| 中文字幕人妻av一区二区| 国产精品美人久久久久久AV| 亚洲久悠悠色悠在线播放| 久久天天躁狠狠躁夜夜2020一| www.av男人.com| 国产区免费| 国产综合欧美| 国产第二十一页| 精品人妻无码区在线视频| 亚洲精品制服丝袜二区| 久久黄色毛片| 国产最爽的乱婬视频国语对白| 91九色视频网| 激情综合五月网| 国产在线精品美女观看| 黄色网在线免费观看| 亚洲成人网在线观看| 狠狠色香婷婷久久亚洲精品| 日韩在线成年视频人网站观看| 精品久久久久久中文字幕女 | 亚洲午夜福利在线| 国产成年女人特黄特色大片免费| 日韩午夜福利在线观看| 国内精品91| 欧美一级在线| 久久www视频| 日韩欧美中文| 在线观看热码亚洲av每日更新| 欧美精品在线视频观看| 国产精品露脸视频| 日本不卡免费高清视频| 亚洲精品天堂自在久久77| 亚洲IV视频免费在线光看| 亚洲国产亚综合在线区| 久久黄色毛片| 有专无码视频| 久草国产在线观看| 亚洲国产系列| 国产精品30p| 日本成人福利视频| 亚洲欧美在线精品一区二区| 狠狠做深爱婷婷久久一区| 中文字幕有乳无码| 激情成人综合网| 亚洲精品黄|