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

    Prediction of Partial Ring Current Index Using LSTM Neural Network*

    2022-11-09 04:25:02LIHuiWANGRunzeWANGChi
    空間科學(xué)學(xué)報(bào) 2022年5期

    LI Hui WANG Runze WANG Chi

    Abstract The local time dependence of the geomagnetic disturbances during magnetic storms indicates the necessity of forecasting the localized magnetic storm indices. For the first time, we construct prediction models for the SuperMAG partial ring current indices (SMR-LT), with the advance time increasing from 1 h to 12 h by Long Short-Term Memory (LSTM) neural network. Generally, the prediction performance decreases with the advance time and is better for the SMR-06 index than for the SMR-00, SMR-12, and SMR-18 index. For the predictions with 12 h ahead, the correlation coefficient is 0.738, 0.608, 0.665, and 0.613, respectively. To avoid the over-represented effect of massive data during geomagnetic quiet periods, only the data during magnetic storms are used to train and test our models, and the improvement in prediction metrics increases with the advance time. For example,for predicting the storm-time SMR-06 index with 12 h ahead, the correlation coefficient and the prediction efficiency increases from 0.674 to 0.691, and from 0.349 to 0.455, respectively. The evaluation of the model performance for forecasting the storm intensity shows that the relative error for intense storms is usually less than the relative error for moderate storms.

    Key words Geomagnetic storm, Partial Ring Current Index (PRCI), Artificial Neural Network (ANN)

    0 Introduction

    Geomagnetic storms are violent global disturbances in the Earth’s magnetosphere, which are one of the consequences of enhanced solar wind-magnetosphere energy coupling through the magnetic reconnection mecha-[1–2]drops rapidly to its minimum value with intense fluctuations (main phase). Finally, along with the decrease of disturbances,Dstrises back to its quiet level slowly (recovery phase).

    Geomagnetic storms have been one of the most challenging and passionately explored topics in the geophysical community. This is not only because geomagnetic storms have a significant impact on the global geomagnetic field pattern, but also because they are the most important link in the solar-terrestrial energy coupling chain. Moreover, geomagnetic storms have serious effects on technological infrastructure such as communication systems, electric power systems, oil pipelines, and space vehicles. Therefore, the prediction of geomagnetic storms has important theoretical and practical value.

    Researches onDstforecasting have been conducted since the 1970s. Burtonet al.[3]developed an empirical algorithm for forecastingDstusing solar wind velocity and density and the north-south component of the IMF, which contains an approximate description of the effects of solar wind dynamic pressure, particle injection from the plasma sheet, and exponential decay, constrained by ordinary differential equations. Based on this model, several additional assumptions are made to characterize this nonlinear behavior[6–8].

    In addition to models constructed based on empirical formulas, there are also researches on forecasting from a data science perspective. One of the widely popular approaches is the construction of machine learning models through Artificial Neural Networks (ANN).Lundstedt and Wintoft[9] developed a feedforward neural network to predict 1 h ahead, using the density and the velocity of the solar wind and the z-component of IMF, which can predict the initial and main phases of magnetic storms. Later an ANN with a specific structure called Recurrent Neural Network (RNN) was introduced to further improve the forecasting effectiveness[10,11]. Different optimization techniques have also been employed for forecasting. Wei et al. [12] used multiscale radial basis networks by changing the activation function of the hidden layer neurons. Lazzús et al. [13]achieved forecasts up to 6 h in advance employing a swarm optimization algorithm. Gruet et al.[14] further achieved better results on forecasting 1~6 h in advance using the Long Short-Term Memory (LSTM) network. Besides, algorithms such ask-Nearest Neighbors(KNN), Nonlinear Auto Regressive Moving Average with eXogenous inputs (NARMAX), and Support Vector Machine (SVM) have been introduced toDstforecasting and other topics in space physics[15–17].

    Although theDstindex has a clear physical interpretation, its implicit simplified assumption of azimuthally invariantHperturbations from low or mid-latitude ground stations is not often consistent with observations. A significant dawn-dusk asymmetry ofHdepression tends to occur during storms (see Liet al.[18],and references therein). Given this observed strong Local Time (LT) variations in theH-component at mid-latitude stations, Newell and Gjerloev[19]introduced the SuperMAG Ring current index (SMR index), and Super-MAG partial Ring current indices (SMR-LT indices).SuperMAG is a worldwide collaboration of organizations and national agencies that currently operate more than 300 ground-based magnetometers. Using 98 low and midlatitude magnetometers, the SMR index is conceptually the same as theDstindex, and its local time version, SMR-LT, including SMR-00, SMR-06, SMR-12, and SMR-18, contains local time components named after their respective local time center point, representing midnight, dawn, noon and dusk. The SMR global index is the arithmetic mean of SMR-LT indices.

    In this work, we first analyze the asymmetry of SMR-LT indices during geomagnetic storms from several different perspectives. The study of the results shows the importance to predict the SMR region index. We then apply an artificial neural network approach, LSTM,to this topic to provide a forecast of the SMR-LT indices 1~ 12 h in advance. To improve the forecasting performance of the storm-time data, we train a new model using only the data during storms and analyze the model’s performance in predicting the intensity of magnetic storms.

    1 Importance of Partial Ring Current Index Prediction

    The hourly partial ring current indices, SMR-LT, and the corresponding averaged index, SMR, are obtained from SuperMAG**https://supermag.jhuapl.edu/indices/. During the concerned period from 1998 to 2019, 318 magnetic storms are analyzed.

    Liet al.[18]presented a significant dawn–dusk asymmetry ofDstindex during the storm main phase and early recovery phase due to the important contribution of the partial ring current during magnetic storm processes.Fig.1 shows the evolution of partial ring current indices during the magnetic storm on 26 August 2018. As expected, the four SMR-LT indices differ significantly in the main phase and early recovery phase. The depression ofH-components of geomagnetic field disturbances on the dusk side, such as SMR-18 and SMR-12,is generally larger than that on the dawn side, such as SMR-00 and SMR-06. The minimum values of SMR-00,SMR-06, SMR-12, and SMR-18 are –156.9 nT, –155.3 nT, –204.4 nT, and –209.3 nT, respectively. Besides, the SMR-18 index is the earliest to reach its minimum,while the SMR-00 index reaches its minimum one hour later. During the later recovery phase, all the indices are almost the same.

    Fig. 1 Partial ring current indices during the magnetic storm on 26 August 2018

    Fig. 2 Statistical characteristics of SMR-LT minimum and its time lag to SMR minimum for 318 magnetic storms

    Fig.2 gives the statistical characteristics of SMRLT minimum and its time lag to SMR minimum for 318 magnetic storms. The top panels present the relationship between the relative intensity of SMR-LT indices during storms and the storm intensity, which are defined by the ratio of SMR-LT minimum to SMR minimum and SMR minimum, respectively. Generally, the relative intensity of SMR-18 is greater than 1, indicating a larger geomagnetic depression, while the relative intensity of SMR-06 is less than 1, indicating a smaller geomagnetic depression. The results for SMR-00 and SMR-12 are similar, with the relative intensity greater than 1 in most cases. The mean minimum value of SMR-00, SMR-06,SMR-12, and SMR-18 for the 318 storms is –111.5 nT,–80.7 nT, –109.4 nT, and –133.8 nT, which agrees with the previous results. In addition, the relative intensity for all the four indices tends to be 1 as the storm intensity increases, indicating that the ring current tends to be more symmetric as suggested by Liet al.[18]. Instead of the storm intensity, the onset of storm minimum for the four partial ring current indices also differs from each other. The bottom panels give the distribution of the time lag of SMR-LT minimum to SMR minimum. For SMR-00, SMR-12, and SMR-18, the time lag is concentrated within 1 h. The median values of the time lag are the same, 0. And the mean value is 1.5, 1.6, and –0.7 h,respectively. For SMR-06, the time lag is greater than 0 in many cases, with the median value of 1 h and the mean value of 2.2 h.

    Thus, predicting the SMR-LT indices instead of the SMR index is supposed to be more appropriate for regional space weather warning or forecasting.

    2 Methodology and Data Set

    2.1 LSTM Networks

    LSTM is a special kind of RNN. Traditional neural networks cannot deal with time series forecasting. However, RNN can address this issue by designing loops in networks, allowing information to persist and be passed from one step of the network to the next step. An RNN can be regarded as multiple copies of the same network,with messages that can be passed from one copy to a successor copy. The chain-like nature reveals that RNN is intimately related to sequences and lists. However,when the gap between the relevant information and the point where it is needed grows too large, RNN becomes unable to learn to connect the information. The problem is called “l(fā)ong-term dependencies”, explored in depth by Bengioet al.[20], who found it arising from vanishing gradient problems during the training phase of RNN.

    LSTM is explicitly designed to avoid the long-term dependency problem. It was introduced by Hochreiter and Schmidhuber[21]and was later refined and popularized by many researchers, making it work tremendously well on a large variety of problems. LSTM also has a chain-like structure, but the repeating module is different. Instead of having a single neural network layer,there are four in LSTM, interacting in a very special approach.

    The key to LSTM is the cell state, which runs straight down the entire chain, with only some minor linear interactions. One end of the chain isCt-1, which represents the old cell state, and the other end isCt, which represents the new cell state. The LSTM is designed to remove or add information to the cell state, precisely regulated by structures called gates. Gates are a way to optionally let information through. They are composed of a sigmoid neural net layer and a pointwise multiplication operation. The sigmoid layer outputs numbers between 0 and 1, describing how much of each component should be let through. A value of 0 means completely rejected, while a value of 1 means perfectly accepted.

    An LSTM consists of three gates to protect and control the cell state, namely forget gate layer, input gate layer and output gate layer. The forget gate layer is a sigmoid layer to decide what information to be forgotten from the cell state. It receivesht-1andxt, which means the output of the precursor copy and the input of time stept, respectively, and outputs a number between 0 and 1, indicating the extent to which the cell stateCt-1is forgotten. For subsequent equations, the notation is kept thatw*andb*represent the corresponding weights and biases of the layer, respectively.

    The next step is to decide what new information to be stored in the cell state, which has two parts. First, a sigmoid layer called the “input gate layer” decides which values to be updated. Next, a tanh layer creates a vector of new candidate values,C~t, that could be added to the state.

    Finally, the output gate layer will produce a filtered version of cell stateCt. First, a sigmoid layer will be run, which decides what parts of the cell state to output. Then, we put the cell state through tanh (to push the values to be between –1 and 1) and multiply it by the output of the sigmoid gate, so that we only output the parts we decided to.

    2.2 Performance Metrics

    In this study, several performance metrics are used to evaluate the effectiveness of machine learning forecasting models. The root-mean-square error (Erms), is used to represent the error between the observed and predicted values.

    The Pearson correlation coefficient,C, is used to represent the linear relationship between the predictions and the observations. 1/–1 implies a perfect positive/negative correlation, while 0 indicates no linear correlation.

    In addition, the prediction efficiency,Pe, is also used and can be calculated as followed.

    Pe=1 indicates a perfect prediction andPe=0 indicates that the performance of the model is equivalent to the situation by the arithmetic mean of the observation.Pecan be negative values, which in this case means that the model prediction is no better than just taking the arithmetic mean value of the test data. In contrast toC,which only reflects the consistency of the trend of the time series,Pecan also indicate the amplitude of the modeled features, allowing to examine the accuracy of the forecast.

    2.3 Data Pre-processing

    The solar wind parameters are obtained from OMNI data**https://omniweb.gsfc.nasa.gov**https://supermag.jhuapl.edu/indices/maintained by the National Space Science Data Center (NSSDC) of the National Aeronautics and Space Administration (NASA). The SMR and SMR-LT indices are obtained from SuperMAG collaborators***https://omniweb.gsfc.nasa.gov**https://supermag.jhuapl.edu/indices/.

    A reasonable selection of feature parameters can make the learning efficiency of ANN greatly improved.In this study, the solar wind fieldEy, the solar wind dynamic pressurePd, the solar wind velocityVp, the solar wind densityNp, the solar wind temperatureTp, the IMF strengthBtand itsz-componentBzare used and denoted byDsw.

    The dataset we select from OMNI and SuperMAG spans the period from 1 January 1998, to 31 December 2019, with a temporal resolution of 1 h. The 22 years dataset is divided into three parts, namely the training set(1998- 2009), the validation set (2010-2014), and the test set (2015-2019). More than a complete solar activity cycle is considered by the training set and the validation set. A total of 257 magnetic storms are identified during these two intervals.

    We aim to forecast the SMR-LT indices several hours in advance, denoted asp, which in terms of data structure means developing an LSTM model for multivariate time series forecasting on multiple lag time steps.In this study,p=1,3,6,12 h are considered. The first step is to prepare the dataset for the LSTM, which involves framing the dataset as a supervised learning problem and normalizing the input variables. We frame the supervised learning problem as predicting the SMR-LT index at the current hour given the solar wind parameters on multiple prior time steps. The OMNI data have already been time-shifted from their location of observation to the Earth’s bow shock, which is desired to best support solar wind –magnetosphere coupling studies.Next, we normalize each feature separately to turn the physical quantities into a stream of data that can be processed by the ANN. After the ANN predicts the results,we then apply the scaling process to obtain the true SMR-LT indices.

    There are different examples of important design patterns for RNN. We adopt a neural network with recurrent connections between hidden units, that reads an entire sequence and then produces a single output, the SMR index or one of the SMR-LT indices. We construct two hidden layers, the first layer containing 50 neurons and the second layer containing 100 neurons.According to the structural characteristics of the data processed by LSTM, each parameter is a multidimensional vector on the time sequence. The length of the time sequence, denoted bys, can significantly affect the model performance. To evaluate this effect, the results obtained under differentsconditions are given in Fig.3.In general, the model performs best whens=3. Thus,we set the length of the time series to be 3 in the following study, and I nput is, therefore, a 3 ×8 matrix.

    Fig. 3 Performance of models in forecasting the SMR index using different sequence length ( s). The four panels represent the results when different p are considered

    3 Results

    3.1 Performances of the SMR-LT Indices

    Predictions

    Given 7 input features of the solar wind, a total of 127 combinations exist. To explore the optimal set of input parameters for eachp, all the possible parameter combinations are tried separately. After comparison, the parameter space,Ey,Pd,Tp,Bt,Bz, performs best and is applied to the prediction.

    Fig.4 gives the scatter plot of the model’s predictions of SMR-LT indices on the test set. The blue line represents the fitting result, and the black dotted line represents the exact predictions. The closer of the blueline to the black dotted line is, the better prediction results we obtain. It is clear that the predictions are very close to the observations whenp=1. Aspincreases, the deviation of the predictions, mostly undervaluation, to the observations increases accordingly.

    Fig. 4 Scatter plot of the model’s predictions of SMR-LT indices on the test set. The blue line represents the fitting result, and the black dotted line represents the exact prediction

    Table 2 gives the quantitative performances of our prediction models for SMR-00, SMR-06, SMR-12, and SMR-18 indices. Considering that there is no previous work on SMR-LT indices predictions, we compare the model predictions with the results of persistence operations to quantitatively represent the model ability. Generally, our model can predict the SMR-LT indices well from the perspectives ofErms,C, andPe. For the prediction of SMR-06 withp=1,3,6,12 h, theErmsis 5.752,6.592, 7.639, and 9.976 nT , respectively, which is 3.4%,17.7%, 22.1%, and 18.0% less than the results of persistence operations; theCis 0.925, 0.873, 0.825, and 0.738,respectively, which is 2.5%, 6.2%, 12.4%, and 25.1%better than the results of persistence operations; thePeis 0.817, 0.759, 0.677, 0.449, respectively, which is 1.6%,17.7%, 44.7%, and 149.4% better than the results of persistence operations. For the prediction of SMR-18 withp=1,3,6,12h, theErmsis 7.316, 10.692, 14.022, and 16.836 nT, respectively, which is 19.3%, 20.9%, 14.0%,and 15.2% less than the results of persistence operations;theCis 0.934, 0.844, 0.743, and 0.613, respectively,which is 4.4%, 10.0%, 12.4%, and 23.1% better than the results of persistence operations; thePeis 0.864, 0.708,0.499, 0.277, respectively, which is 9.4%, 32.6%,55.0%, and 5440.0% better than the results of persistence operations. For the predictions of SMR-00 and SMR-12 index, the results are similar. It is clear that the prediction performances decrease withp, however, the improvements to persistence significantly increase withp.

    Table 1 Comparison between Ref. [14], Ref. [22],Ref. [23], and our proposed model

    3.2 Performances of SMR-LT Index Forecast during Magnetic Storms

    The problem of the underestimated predictions shown in Fig.4 might come from the training data structure. In the previous section, all of the dataset is used to train our models. The whole dataset contains a total of 192839 data, of which only 65071 data, or about 33.74%, are during 318 magnetic storms. In the case of intense storms with an intensity of less than –100 nT, there are only 6452 data, accounting for about 3.35%. The overweighting of data during magnetic quiet periods, which is not the focus of the study, might affect the preference of the trained machine learning model, making it always biased to underestimate the intensity of magnetic storms and becomes more prominent for a longer forecast lead time. Thus, although governed by the same fundamental laws of physics, it is still necessary to separate the data during magnetic storms from all-time series from the perspective of data science. Note that, the parameter space with the best performance changes toEy,Pd,Vp,Np,Bt,Bz. Fig.5 shows an example of the predicted SMR-LT indices during the storm on 26 August 2018.Whenp=1,3 h, the predictions are similar to the observations, especially for the SMR-00 and SMR-06 index.Whenp=6,12 h, the deviations of the predictions from the observations are more significant.

    Fig. 5 Evolution of observed (in black) and predicted (in blue) SMR-LT indices during the storm on 26 August 2018

    Table 3 gives the quantitative performances of our prediction models for SMR-00, SMR-06, SMR-12, and SMR-18 indices during magnetic storms. “All-time”represents that all the dataset are used to train the model and the metrics are calculated based on the predictions during magnetic storms. However, “storm” represents that only the dataset during magnetic storms is used to train the model and the metrics are calculated based on the predictions during magnetic storms as well. Compared to the results shown in Table 2, theCandPeof both “persistence” and “all-time” have different degrees of reductions due to the removal of the data during magnetic quiet periods, and the reductions increase with the lead time of prediction. At the same time, the results of “storm” are better than the results of “alltime”, and the performance improvement increases with the lead time of prediction. For the prediction of SMR-06 withp=1,3,6,12 h, the improvement ofCis 0.2%,0.3%, 1.4%, and 2.5%, respectively; the improvement ofPeis 3.6%, 1.5%, 4.0%, and 30.4%, respectively. For the prediction of SMR-18 withp=1,3,6,12 h, the improvement ofCis 0.1%, 1.1%, 2.5%, and 8.3%, respectively; the improvement ofPeis –0.1%, –0.6%, 12.2%,and 101.9%, respectively. For SMR-00 and SMR-12 index, the improvements are similar. It is clear that the prediction performances decrease withp, however, the improvements significantly increase withp.

    Table 2 Performance of SMR-LT indices forecast

    Table 3 Performance of storm-time SMR-LT indices prediction

    3.3 Performance of Storm Intensity Prediction

    Table 4 gives the distribution of different magnetic storms intensity on the dataset, and Table 5 gives the mean values and standard deviations of the relative errors of storm intensity predictions. There are 61 magnetic storms during the test dataset. All these storms are divided into two groups, moderate storms with the minimum of SMR-LT index between –50 nT and –100 nT,and intense storms with the minimum of SMR-LT index less than –100 nT. Negative values represent underestimations of the storm intensity, while positive values represent overestimations. Note that, the storm intensities in this subsection are discussed in absolute values. It is clear that the relative errors tend to increase with the lead time of predictions. Taking the predictions for SMR-06 index during moderate storms, for example, the relative errors are –10.99%, –13.73%, and –16.88% forp=1,3,6h. During moderate storms, the smallest relative error forp=1 h is for SMR-12 index prediction,–9.81%. Forp=3,6 h, the smallest relative errors are–13.73% and –16.88%, both for SMR-06 index prediction. In contrast, the relative errors for SMR-00 and SMR-18 forecasts are relatively larger, as they tend to have a larger magnitude of variation and are therefore more difficult to predict accurately. In addition, the relative errors for intense storms are relatively less than those for moderate storms.

    Table 4 Magnetic storm intensity distribution

    Table 5 Relative errors (mean±standard deviation) of storm intensity predictions

    4 Summary

    In this paper, we construct models for predicting partial ring current index, SMR-LT, with an advance time from 1 h to 12 hviaLSTM neural networks. Although theDstindex is widely used to present the averaged magnetic disturbances on a global scale, the local time dependence of geomagnetic disturbances during magnetic storms indicates the importance of forecasting localized magnetic storm indices, especially for studying the effects of magnetic storms on a particular region.

    Using LSTM neural network, we construct prediction models of partial ring current index, SMR-LT, with an advance time from 1 h to 12 h for the first time. From 127 combinations of 7 solar wind parameters, the parameter space,Ey,Pd,Tp,Bt,Bz, performs best. The performances of the SMR index prediction are comparable with the results of the published models, which validates the reliability of our model when applying it to predict the SMR-LT indices. Generally, the prediction performance decreases with the advance time and is better for SMR-06 than for SMR-00, SMR-12, and SMR-18. For the prediction of SMR-00, SMR-06, SMR-12,and SMR-18 with 12 h ahead, the correlation coefficient is 0.608, 0.738, 0.665, and 0.613, respectively.

    The evaluation of the model results shows that the over-represented (over 66.26%) data during magnetic quiet periods might affect the learning preference of the ANN. To avoid the over-represented effect, we filter out the data during magnetic quiet periods from the dataset and train a new model, labeled as the “storm” model.The performance improvement increases with the advance time. Taking the prediction of storm-time SMR-06 index with a 12 h advance time, for example, the correlation coefficient and the prediction efficiency increa-ses from 0.674 to 0.691, and from 0.349 to 0.455, respectively.

    We also evaluate the model performance for forecasting storm intensity. The relative errors tend to increase with the lead time of predictions. The relative error of storm intensity prediction for intense storms is usually less than that for moderate storms.

    AcknowledgementThe authors thank the SuperMAG collaborators for the use of SMR and SMR-LT indices, which are accessible in https://supermag.jhuapl.edu/indices/ . The OMNI data were obtained from the GSFC/SPDF OMNIWeb interface at https://omniweb.gsfc.nasa.gov. Special thanks to Dr. Tang B B, Dr. Li W Y, and Prof. Guo X C for helpful discussions.

    久久精品国产亚洲av香蕉五月| .国产精品久久| 欧美激情在线99| 69人妻影院| 深夜a级毛片| av福利片在线观看| 久久精品国产鲁丝片午夜精品| 成人鲁丝片一二三区免费| 在线观看午夜福利视频| 日韩成人av中文字幕在线观看 | 久久人人爽人人爽人人片va| 日本五十路高清| 欧美三级亚洲精品| 夜夜看夜夜爽夜夜摸| 国产精品伦人一区二区| 夜夜夜夜夜久久久久| 97超级碰碰碰精品色视频在线观看| 日韩欧美免费精品| 日韩强制内射视频| 国产av不卡久久| 日本精品一区二区三区蜜桃| 免费电影在线观看免费观看| 久久中文看片网| 国产成年人精品一区二区| 免费av毛片视频| 国产高清激情床上av| 亚洲成人中文字幕在线播放| 亚洲精品日韩av片在线观看| 又粗又爽又猛毛片免费看| 麻豆久久精品国产亚洲av| 不卡视频在线观看欧美| 久久久久久久久久成人| 午夜福利在线在线| 国产亚洲精品综合一区在线观看| 成人精品一区二区免费| 午夜免费男女啪啪视频观看 | 亚洲成人中文字幕在线播放| 日韩中字成人| 男人舔女人下体高潮全视频| 久久久久久伊人网av| 国产成年人精品一区二区| 精品久久国产蜜桃| 男女之事视频高清在线观看| 亚洲aⅴ乱码一区二区在线播放| 久久午夜福利片| 国产亚洲精品久久久com| 久久九九热精品免费| 51国产日韩欧美| 久久精品国产99精品国产亚洲性色| 亚洲美女搞黄在线观看 | 精品人妻视频免费看| 久久精品夜夜夜夜夜久久蜜豆| av免费在线看不卡| 国产精品人妻久久久影院| 97在线视频观看| 最近视频中文字幕2019在线8| 校园人妻丝袜中文字幕| 日日摸夜夜添夜夜爱| 国产精品久久视频播放| 99在线视频只有这里精品首页| 99热全是精品| 亚洲一区高清亚洲精品| 亚洲不卡免费看| 国产乱人偷精品视频| 久久久久精品国产欧美久久久| 大又大粗又爽又黄少妇毛片口| 国产在线男女| 日韩国内少妇激情av| 国产精品1区2区在线观看.| 高清午夜精品一区二区三区 | 午夜福利高清视频| 亚洲成人av在线免费| 国产高潮美女av| 亚洲av中文av极速乱| www日本黄色视频网| 日本精品一区二区三区蜜桃| 久久久精品欧美日韩精品| 少妇猛男粗大的猛烈进出视频 | 久久久国产成人免费| 国产日本99.免费观看| 校园春色视频在线观看| 久久久精品大字幕| 看黄色毛片网站| 欧美日本视频| 欧美最新免费一区二区三区| 男人舔奶头视频| 国产精品久久久久久精品电影| 亚洲av成人av| 深爱激情五月婷婷| 国产成人a∨麻豆精品| 国产淫片久久久久久久久| 免费看av在线观看网站| 精品人妻一区二区三区麻豆 | 日本五十路高清| 99热网站在线观看| 久久久精品欧美日韩精品| 国产成人影院久久av| 秋霞在线观看毛片| 久久精品国产99精品国产亚洲性色| 精品久久久久久久久久免费视频| 中国美白少妇内射xxxbb| av福利片在线观看| 成人毛片a级毛片在线播放| 亚洲成人中文字幕在线播放| 精品午夜福利在线看| 狂野欧美激情性xxxx在线观看| 六月丁香七月| 少妇熟女欧美另类| 久久久久久久久久久丰满| 99精品在免费线老司机午夜| 精品少妇黑人巨大在线播放 | 99热6这里只有精品| 老司机福利观看| 国产精品野战在线观看| 看十八女毛片水多多多| 欧美一区二区精品小视频在线| 亚洲精品日韩av片在线观看| 国产男人的电影天堂91| 18禁在线播放成人免费| 久久午夜福利片| 成人鲁丝片一二三区免费| 亚洲三级黄色毛片| 欧美绝顶高潮抽搐喷水| 国产极品精品免费视频能看的| 自拍偷自拍亚洲精品老妇| 男女做爰动态图高潮gif福利片| 国产国拍精品亚洲av在线观看| 欧美zozozo另类| 欧美xxxx性猛交bbbb| 中文字幕人妻熟人妻熟丝袜美| 国产亚洲精品av在线| 国产午夜福利久久久久久| 狂野欧美白嫩少妇大欣赏| 一级黄片播放器| 日本三级黄在线观看| 亚洲性久久影院| 精品不卡国产一区二区三区| av在线观看视频网站免费| 亚洲美女视频黄频| 91久久精品国产一区二区三区| 狠狠狠狠99中文字幕| 久久久国产成人免费| 亚洲人成网站在线播| 赤兔流量卡办理| 女同久久另类99精品国产91| 久久国产乱子免费精品| 变态另类丝袜制服| 六月丁香七月| 免费看av在线观看网站| 黄色欧美视频在线观看| av天堂中文字幕网| 婷婷精品国产亚洲av在线| 黑人高潮一二区| 女人十人毛片免费观看3o分钟| 在线免费观看不下载黄p国产| 久久久久久久久久久丰满| 久久亚洲国产成人精品v| 色哟哟哟哟哟哟| 国产成人91sexporn| 中文字幕熟女人妻在线| 99久久精品国产国产毛片| 国产精品久久久久久久久免| 国内精品久久久久精免费| 欧美不卡视频在线免费观看| 日本五十路高清| 99在线人妻在线中文字幕| 18禁在线无遮挡免费观看视频 | 亚洲中文字幕一区二区三区有码在线看| 男人舔女人下体高潮全视频| 噜噜噜噜噜久久久久久91| 男人和女人高潮做爰伦理| 免费看日本二区| 国产伦精品一区二区三区四那| 国产精品1区2区在线观看.| av在线观看视频网站免费| 麻豆国产97在线/欧美| 18禁裸乳无遮挡免费网站照片| 麻豆av噜噜一区二区三区| 真人做人爱边吃奶动态| 身体一侧抽搐| 国产精品国产三级国产av玫瑰| 一个人观看的视频www高清免费观看| 偷拍熟女少妇极品色| 国产淫片久久久久久久久| 亚洲第一电影网av| 中文字幕免费在线视频6| 精品少妇黑人巨大在线播放 | 欧美一区二区国产精品久久精品| 18禁在线播放成人免费| 最近在线观看免费完整版| 内地一区二区视频在线| 欧美日韩国产亚洲二区| 精品人妻偷拍中文字幕| 在线免费观看的www视频| 亚洲高清免费不卡视频| 麻豆av噜噜一区二区三区| 看黄色毛片网站| 国内少妇人妻偷人精品xxx网站| 高清日韩中文字幕在线| 国产成人a∨麻豆精品| 极品教师在线视频| 国产精品人妻久久久久久| 国产免费一级a男人的天堂| 日韩精品有码人妻一区| 日韩欧美免费精品| 久久久久性生活片| 人妻丰满熟妇av一区二区三区| 女同久久另类99精品国产91| 日本 av在线| 超碰av人人做人人爽久久| 少妇人妻精品综合一区二区 | 在线免费观看不下载黄p国产| 久久久精品大字幕| 全区人妻精品视频| 久久99热这里只有精品18| 狠狠狠狠99中文字幕| 美女高潮的动态| 成人精品一区二区免费| 我的老师免费观看完整版| 国产精品,欧美在线| 午夜老司机福利剧场| 九九爱精品视频在线观看| 国产精品一区二区三区四区免费观看 | 国产亚洲精品久久久久久毛片| 搡女人真爽免费视频火全软件 | 精品久久久久久久久久久久久| 久久鲁丝午夜福利片| 美女大奶头视频| 激情 狠狠 欧美| 国产精品一区www在线观看| 青春草视频在线免费观看| av在线蜜桃| 亚洲av免费在线观看| 最后的刺客免费高清国语| 国产精品永久免费网站| 成年版毛片免费区| 18禁黄网站禁片免费观看直播| 久久综合国产亚洲精品| 亚洲成a人片在线一区二区| 美女被艹到高潮喷水动态| 亚洲一区二区三区色噜噜| 久久久久久久久久黄片| 午夜久久久久精精品| 在线观看66精品国产| 又粗又爽又猛毛片免费看| 变态另类丝袜制服| 人妻少妇偷人精品九色| 国产探花在线观看一区二区| 亚洲中文字幕一区二区三区有码在线看| 亚洲av美国av| 欧美色视频一区免费| 男女那种视频在线观看| 欧美日韩一区二区视频在线观看视频在线 | 国产成人a∨麻豆精品| 极品教师在线视频| 成人高潮视频无遮挡免费网站| 日本色播在线视频| 超碰av人人做人人爽久久| 亚洲一区高清亚洲精品| 久久精品国产亚洲av涩爱 | 男人舔女人下体高潮全视频| av视频在线观看入口| 免费看日本二区| 亚洲无线观看免费| 亚洲真实伦在线观看| 激情 狠狠 欧美| 成年女人毛片免费观看观看9| 久久鲁丝午夜福利片| 少妇丰满av| 性插视频无遮挡在线免费观看| 内地一区二区视频在线| 看十八女毛片水多多多| 日韩精品有码人妻一区| 69av精品久久久久久| 特大巨黑吊av在线直播| 国产毛片a区久久久久| 日韩精品有码人妻一区| 日日干狠狠操夜夜爽| 国产单亲对白刺激| 不卡视频在线观看欧美| 国产午夜精品久久久久久一区二区三区 | 国产真实伦视频高清在线观看| 欧美日韩乱码在线| 晚上一个人看的免费电影| 两性午夜刺激爽爽歪歪视频在线观看| 亚洲自拍偷在线| 偷拍熟女少妇极品色| 久久草成人影院| 白带黄色成豆腐渣| 一a级毛片在线观看| 久久韩国三级中文字幕| 亚洲精品乱码久久久v下载方式| 欧美激情久久久久久爽电影| 2021天堂中文幕一二区在线观| 日本免费a在线| 欧美极品一区二区三区四区| 伦理电影大哥的女人| 色综合亚洲欧美另类图片| 美女cb高潮喷水在线观看| 真实男女啪啪啪动态图| 寂寞人妻少妇视频99o| 亚洲熟妇中文字幕五十中出| 精品一区二区三区视频在线| 国产激情偷乱视频一区二区| 天天躁夜夜躁狠狠久久av| a级毛片a级免费在线| 蜜桃亚洲精品一区二区三区| 又黄又爽又免费观看的视频| 久久亚洲精品不卡| 亚洲五月天丁香| 精品久久久久久久久久久久久| 12—13女人毛片做爰片一| 97在线视频观看| 啦啦啦啦在线视频资源| 久久久久国产网址| 男女做爰动态图高潮gif福利片| 久久久久九九精品影院| 亚洲三级黄色毛片| 又粗又爽又猛毛片免费看| 日韩欧美免费精品| 蜜桃久久精品国产亚洲av| 熟女电影av网| 少妇的逼好多水| 中文字幕精品亚洲无线码一区| 中文字幕av在线有码专区| 此物有八面人人有两片| 综合色丁香网| 欧美国产日韩亚洲一区| 成人亚洲欧美一区二区av| 亚洲国产欧洲综合997久久,| 国产片特级美女逼逼视频| 赤兔流量卡办理| 蜜臀久久99精品久久宅男| 亚洲av第一区精品v没综合| 国产91av在线免费观看| 国产熟女欧美一区二区| 欧美精品国产亚洲| 精品人妻熟女av久视频| 欧美人与善性xxx| 精品国产三级普通话版| 男人和女人高潮做爰伦理| 久久亚洲国产成人精品v| avwww免费| 麻豆国产av国片精品| 在线国产一区二区在线| 3wmmmm亚洲av在线观看| 亚洲五月天丁香| 亚洲在线观看片| 国产午夜精品久久久久久一区二区三区 | 变态另类丝袜制服| 欧美xxxx性猛交bbbb| 99久久九九国产精品国产免费| 高清日韩中文字幕在线| 狂野欧美白嫩少妇大欣赏| 国内精品一区二区在线观看| 丰满人妻一区二区三区视频av| 99在线视频只有这里精品首页| 国产欧美日韩精品一区二区| 国产男靠女视频免费网站| 69人妻影院| 亚洲av成人精品一区久久| 人人妻人人澡欧美一区二区| 亚洲婷婷狠狠爱综合网| 国产在线精品亚洲第一网站| 日韩av不卡免费在线播放| 亚洲aⅴ乱码一区二区在线播放| 亚洲在线自拍视频| 国产av在哪里看| 欧美xxxx黑人xx丫x性爽| 99精品在免费线老司机午夜| 久久人人精品亚洲av| 亚洲高清免费不卡视频| 国产探花在线观看一区二区| 午夜影院日韩av| 性插视频无遮挡在线免费观看| 国产综合懂色| 久久热精品热| 精品免费久久久久久久清纯| 精品人妻一区二区三区麻豆 | 露出奶头的视频| 老司机午夜福利在线观看视频| 欧美三级亚洲精品| 在线观看美女被高潮喷水网站| 国产av一区在线观看免费| 日韩人妻高清精品专区| 日韩av在线大香蕉| 中文字幕久久专区| 日日啪夜夜撸| 18禁在线无遮挡免费观看视频 | 深夜精品福利| 秋霞在线观看毛片| 黄色日韩在线| 免费高清视频大片| 成人二区视频| 精品熟女少妇av免费看| 男女之事视频高清在线观看| 日本三级黄在线观看| 又粗又爽又猛毛片免费看| 99久久无色码亚洲精品果冻| 国产激情偷乱视频一区二区| 热99在线观看视频| 国产欧美日韩一区二区精品| 亚洲成av人片在线播放无| 99热网站在线观看| 99热这里只有是精品在线观看| 欧美最新免费一区二区三区| 精品久久国产蜜桃| 国产精品久久电影中文字幕| 久久热精品热| 一级毛片久久久久久久久女| 免费黄网站久久成人精品| 又粗又爽又猛毛片免费看| 一个人观看的视频www高清免费观看| 国产av麻豆久久久久久久| 99久久精品热视频| 秋霞在线观看毛片| 插逼视频在线观看| 深爱激情五月婷婷| 国产在线精品亚洲第一网站| 少妇人妻一区二区三区视频| 99久久精品国产国产毛片| 久久久久国产精品人妻aⅴ院| 国产亚洲91精品色在线| 蜜桃亚洲精品一区二区三区| 三级国产精品欧美在线观看| 丝袜美腿在线中文| 少妇的逼水好多| 国产私拍福利视频在线观看| 久久人人爽人人片av| 寂寞人妻少妇视频99o| 69人妻影院| 日日摸夜夜添夜夜添小说| 久久这里只有精品中国| 免费看av在线观看网站| 成人毛片a级毛片在线播放| 精品久久久久久久末码| 久久久久免费精品人妻一区二区| 色吧在线观看| 欧美成人精品欧美一级黄| av天堂在线播放| 好男人在线观看高清免费视频| 白带黄色成豆腐渣| 国产精品精品国产色婷婷| 国产免费一级a男人的天堂| 三级国产精品欧美在线观看| 日本a在线网址| 真人做人爱边吃奶动态| 久久久久久伊人网av| 又爽又黄无遮挡网站| 中文字幕av成人在线电影| 亚洲av免费高清在线观看| 一个人看视频在线观看www免费| 亚洲国产欧美人成| 国产大屁股一区二区在线视频| avwww免费| 欧美中文日本在线观看视频| 91久久精品电影网| 国产精品久久视频播放| 一个人看视频在线观看www免费| 女生性感内裤真人,穿戴方法视频| 尤物成人国产欧美一区二区三区| 亚洲乱码一区二区免费版| 人妻少妇偷人精品九色| 亚洲成人av在线免费| 99国产精品一区二区蜜桃av| 人人妻人人看人人澡| 久久久久久久久大av| 久久久国产成人精品二区| 国产乱人视频| 春色校园在线视频观看| 国语自产精品视频在线第100页| 美女xxoo啪啪120秒动态图| 一级毛片aaaaaa免费看小| 国产精品久久久久久精品电影| 特大巨黑吊av在线直播| 最好的美女福利视频网| 最近2019中文字幕mv第一页| 高清午夜精品一区二区三区 | 日韩高清综合在线| 伦理电影大哥的女人| 少妇裸体淫交视频免费看高清| 在线国产一区二区在线| 国产精品三级大全| 国产色爽女视频免费观看| 欧美bdsm另类| 桃色一区二区三区在线观看| 国内少妇人妻偷人精品xxx网站| 国产亚洲av嫩草精品影院| 女同久久另类99精品国产91| 99国产精品一区二区蜜桃av| 天天一区二区日本电影三级| 听说在线观看完整版免费高清| 黄色欧美视频在线观看| 日日撸夜夜添| 午夜福利成人在线免费观看| 无遮挡黄片免费观看| 亚洲人与动物交配视频| 悠悠久久av| 亚洲欧美中文字幕日韩二区| 国产视频一区二区在线看| 卡戴珊不雅视频在线播放| 91狼人影院| 欧美性感艳星| 久久久久免费精品人妻一区二区| 国产一区二区三区在线臀色熟女| 日韩一区二区视频免费看| 久久人人爽人人爽人人片va| 男女那种视频在线观看| 色综合亚洲欧美另类图片| 日本三级黄在线观看| 久久精品久久久久久噜噜老黄 | 亚洲熟妇熟女久久| 亚洲精品456在线播放app| 亚洲中文字幕日韩| 日韩欧美精品免费久久| 白带黄色成豆腐渣| 色播亚洲综合网| 久久鲁丝午夜福利片| 天美传媒精品一区二区| 国产免费男女视频| 黄色视频,在线免费观看| 国产亚洲精品av在线| 俄罗斯特黄特色一大片| 可以在线观看毛片的网站| 美女黄网站色视频| 亚洲第一电影网av| 国产精品久久久久久精品电影| 色尼玛亚洲综合影院| 日韩精品有码人妻一区| 欧洲精品卡2卡3卡4卡5卡区| 在线观看一区二区三区| 91精品国产九色| 深夜a级毛片| 欧美丝袜亚洲另类| 久久久久久九九精品二区国产| 精品欧美国产一区二区三| 高清午夜精品一区二区三区 | 插阴视频在线观看视频| 成人特级黄色片久久久久久久| 国产成人影院久久av| 成人特级黄色片久久久久久久| 99久国产av精品国产电影| 色哟哟哟哟哟哟| 欧美色视频一区免费| 12—13女人毛片做爰片一| 国产人妻一区二区三区在| 亚洲欧美精品自产自拍| 国产欧美日韩一区二区精品| 我的老师免费观看完整版| 午夜精品在线福利| 中文亚洲av片在线观看爽| 亚洲av免费高清在线观看| 免费观看在线日韩| 亚洲av免费高清在线观看| 国产免费一级a男人的天堂| 99在线视频只有这里精品首页| 欧美一区二区精品小视频在线| 人妻少妇偷人精品九色| 亚洲精品亚洲一区二区| 欧美色欧美亚洲另类二区| 最好的美女福利视频网| 99久国产av精品| 日韩精品有码人妻一区| 亚洲美女搞黄在线观看 | 亚洲图色成人| 久久婷婷人人爽人人干人人爱| av黄色大香蕉| 精品午夜福利在线看| 中文字幕久久专区| 成年女人永久免费观看视频| 国产伦精品一区二区三区四那| 久久久久久大精品| 亚洲无线观看免费| 最近中文字幕高清免费大全6| 国产精品一区二区三区四区免费观看 | 午夜激情欧美在线| 人人妻,人人澡人人爽秒播| 久久人人爽人人片av| 日产精品乱码卡一卡2卡三| 成人鲁丝片一二三区免费| 午夜影院日韩av| 午夜精品国产一区二区电影 | 精品久久久久久成人av| 男女下面进入的视频免费午夜| 久久精品国产鲁丝片午夜精品| 国产精品,欧美在线| 青春草视频在线免费观看| 国产一区二区亚洲精品在线观看| 波多野结衣巨乳人妻| 一级毛片久久久久久久久女| 婷婷亚洲欧美| 俺也久久电影网| 一区二区三区高清视频在线| 久久久国产成人免费| 男人舔奶头视频| 97超级碰碰碰精品色视频在线观看| 亚洲精品日韩av片在线观看| 欧美极品一区二区三区四区| 久久精品国产亚洲av香蕉五月| 最近最新中文字幕大全电影3| 国产高清视频在线观看网站| 看黄色毛片网站| 校园人妻丝袜中文字幕| 亚洲最大成人中文| 18禁在线无遮挡免费观看视频 | 日韩高清综合在线| av在线蜜桃| 麻豆av噜噜一区二区三区| 久久亚洲精品不卡| 啦啦啦观看免费观看视频高清| 淫秽高清视频在线观看| 黄色一级大片看看| 久久久久九九精品影院| 免费观看在线日韩| 麻豆乱淫一区二区|