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

    Development of a 2D spatial displacement estimation method for turbulence velocimetry of the gas puff imaging system on EAST

    2024-04-06 07:15:52LetianLI李樂天ShaochengLIU劉少承NingYAN顏寧XiaojuLIU劉曉菊andXiangGAO高翔
    Plasma Science and Technology 2024年3期
    關(guān)鍵詞:顏寧高翔樂天

    Letian LI (李樂天) ,Shaocheng LIU (劉少承) ,Ning YAN (顏寧) ,Xiaoju LIU (劉曉菊) and Xiang GAO (高翔)

    1 Anhui University,Hefei 230039,People’s Republic of China

    2 Institute of Plasma Physics,Chinese Academy of Sciences,Hefei 230031,People’s Republic of China

    Abstract A gas puff imaging (GPI) diagnostic has been developed and operated on EAST since 2012,and the time-delay estimation (TDE) method is used to derive the propagation velocity of fluctuations from the two-dimensional GPI data.However,with the TDE method it is difficult to analyze the data with fast transient events,such as edge-localized mode (ELM).Consequently,a method called the spatial displacement estimation (SDE) algorithm is developed to estimate the turbulence velocity with high temporal resolution.Based on the SDE algorithm,we make some improvements,including an adaptive median filter and super-resolution technology.After the development of the algorithm,a straight-line movement and a curved-line movement are used to test the accuracy of the algorithm,and the calculated speed agrees well with preset speed.This SDE algorithm is applied to the EAST GPI data analysis,and the derived propagation velocity of turbulence is consistent with that from the TDE method,but with much higher temporal resolution.

    Keywords: gas puff imaging,spatial displacement estimation,SDE,edge turbulence velocity,TDE,EAST tokamak

    1.Introduction

    The edge turbulence can contribute to remarkable cross-field transport in magnetically confined plasmas,such as blobs or filaments [1-3].During the propagation of blobs from the plasma edge to the scrape-off layer (SOL),a large amount of heat and numerous particles across the last closed flux surface (LCFS) can reach the divertor plates in a short time and lead to heat deposition on the first plasma-facing materials.Some fast transient events,such as edge-localized mode(ELM) can cause significant energy loss in the core plasma and high heat flux at the divertor target,which is a crucial issue for large fusion devices such as ITER.Consequently,the fine dynamic evolution of the edge turbulence is important for understanding edge cross-field transport and core plasma confinement.

    Two-dimensional (2D) imaging technology is applied to the measurement of edge turbulence structure and its propagation velocity.Several 2D diagnostics for turbulence structure have been developed,including beam emission spectroscopy (BES) [4-6 ],Langmuir probe [7-9],electron cyclotron emission imaging (ECEI) [10-14] and gas puff imaging (GPI) [15-17].The Langmuir probe array inserted into the edge plasma can measure floating potentials and ion saturation current directly,providing the turbulence information with high temporal evolution [18].The BES and ECEI can measure the edge 2D turbulence structure via the light emission and microwave,respectively.The GPI measures the 2D edge turbulence by injecting neutral gas (helium or deuterium) and capturing the line emission with a high-speed camera.Many tokamaks are equipped with GPI diagnostics,such as NSTX [19],Alcator C-Mod [20],TEXTOR [21],HL-2A [17] and J-TEXT [22].

    Since the turbulence flow velocity is an essential parameter for the understanding of turbulence instability and radial transport,various methods are used in velocity estimation,such as the time-delay method (TDE) [23],optical flow [24],orthogonal dynamics programming (ODP) [25] and spatial displacement estimation (SDE) [26].For example,velocity analysis for the BES diagnostic on DIII-D is performed using the ODP method [27].Velocity calculations for the Langmuir probe on NSTX and ECEI diagnostic on KSTAR are done using TDE estimation [28,29].In the TDE process,cross-correlation is performed between a reference point and the surrounding points in a series of time sequences,and the displacement is given by the distance between the two points with the maximum cross-correlation coefficient.Finally,the speed is estimated by the ratio of displacement and delay time.The velocity calculated by the TDE method has high precision,and its spatial resolution can reach optical resolution.However,time series data are required for TDE analysis,and it is difficult for this analysis to provide precise velocity in a situation with large-scale motion during a short time period.The optical flow method estimates the speed by detecting the light change on a target object.This method requires constant brightness on the target object,which is a challenge for the GPI image.In addition,if the movement is too fast and the specific structure exceeds the observation window,the optical flow method will fail to attain the correct speed,which is called the aperture problem.The SDE method calculates the spatial displacement from the spatial lag of the 2D spatial cross-correlation coefficient function.Since the SDE can be operated for two adjacent images,it has much higher temporal resolution than the TDE method,which is an advantage for the velocity estimation of fast structures.Due to the SDE method needing enough large image pixels,its spatial resolution is limited.

    Some improvements are added to the SDE algorithm that is applied to the GPI data analysis on EAST.First,we introduce an adaptive median filter [30],which can effectively solve the filtering problem when the noise is relatively high.The window size of the median filter can be dynamically modified according to the preset conditions,aiming to denoise and simultaneously preserve the detailed information of the signal.Second,the super-resolution reconstruction by bicubic upsampling [31] is used to acquire high spatial resolution image data and easily extract the detailed structure.Third,the location of a specific structure is captured automatically using an appropriate window size according to the structure's peak.When we finish the development of the SDE algorithm,it is tested by two given motion trajectories,exhibiting accurate moving speed.Finally,the SDE method is applied to the GPI data analysis,and the calculated velocity agrees well with that from the TDE method.The imaged frame can be split into 20×30 pixel2sub-frames.Therefore,the velocity distribution in the poloidal-radial plane is obtained.The rest of this paper is organized as follows.Section 2 is the development of the SDE algorithm,and the corresponding test and application on GPI data are presented in section 3 .Section 4 is a summary.

    2.Development of SDE algorithm

    In this section,we introduce the SDE algorithm in detail.The algorithm includes the following steps: data pre-processing,image enhancement,speed estimation.The purpose of data pre-processing is to remove background noise and the motion trend.Image enhancement can improve the quality and visibility of images that are easy to analyze and understand.Adaptive median filtering can help us to reasonably process the signal noise under different noise levels,preserve image details and make the subsequent calculations more accurate.Super-resolution processing can increase the image’s resolution,thereby improving the details and clarity of the picture.Using super-resolution processing technology,we can convert a low-resolution image into a high-resolution image,consequently enhance the quality of the image.The third step is to convert the time-domain signal into a frequency-domain signal through fast Fourier transform(FFT) [32],perform cross-correlation operations,and obtain the displacement from the polynomial fitting.Finally,the velocity is calculated by the relative displacement and the delay time.

    2.1.Pre-processing

    In this part,we preprocess the GPI data.First,since the background radiation will disturb the turbulence motion in GPI images,the background emission should be removed in the pre-processing stage.The background emission is given by the average of a series of frames before the helium gas reaches the edge plasma,and then every frame subtracts the background emission to remove the background [15].The original GPI image is shown in figure 1(a),and the background light radiation intensity before gas puffing is presented in figure 1(b).When we remove the background emission intensity from the original 2D signal,we obtain figure 1(c).Finally the fluctuations of the GPI emission intensity are derived by subtracting the averaged emission intensity of the interested GPI frames (t=3.466925-3.467121s) from figure 1(c),as shown in figure 1(d).

    Detrending refers to removing some information in the image sequence that affects the speed estimation [33],such as signal offset.Detrending is an essential part of speed estimation because it can further improve the velocity estimation accuracy.After the background removal step,we detrend the data usingnth-order polynomial fitting,as presented in reference [26].The images before and after the polynomial detrending are shown in figures 2(a) and (b) for discharge #102993,respectively.The trend is much lower than the signal level,as illustrated in figure 2(c).

    Figure 1.(a) Original GPI image measured in discharge #102993,(b) background emission intensity before gas puffing,(c) the same image after removing the background emission,(d) the image with the averaged emission intensity subtracted.

    Figure 2.(a) The image before detrending,(b) the image after detrending,(c) the trend of image calculated by the nth-order polynomial fitting.

    2.2.Image enhancement

    2.2.1.Adaptive median filteringThe median filter is widely used in image processing [34].The median filter works by calculating the median value of the pixels in the neighborhood of each pixel and assigning this value to the central pixel.Compared to other filtering methods,the median filter is an effective way to remove noise from images and simultaneously preserve edge information without remarkable blurring.This makes it a widely-used tool in the field of image processing.

    If a larger filter window is used for the median filter,it can effectively remove noise,but it also leads to a fuzzy image that is unfavorable in the subsequent calculations.Conversely,using a smaller filter window can help preserve more details and enhance the clarity of the image,but it is not an effective way to remove noise.Therefore,an adaptive median filter is adopted to dynamically modify the window size based on the noise level [30].

    Figure 3.Workflow diagram of the adaptive median filter.

    The adaptive median filter can dynamically change the window size of the filter according to the preset conditions.The filtering is performed on a point (x,y) of the image and inside a rectangular windowWx,y,with an allowed maximum window size ofWmax.Gx,yrepresents the gray value of point (x,y) .Gmed,GminandGmaxare the median,minimum gray value and maximum gray value in the window (excluding point (x,y)),respectively.The working principle of the adaptive median filter is illustrated in figure 3.In the first diamond,we use the intensity ofG1andG2to determine whether the gray median value in the window is too large or too small.IfGmedis betweenGminandGmax,the window size is appropriate.Next,we determine whether the valueGx,yis noise;ifGx,yis beyond the range fromGmintoGmax,the window size should be expanded.As shown in the second diamond,if the gray valueGx,yis also within a reasonable range,we consider that this point is not noise,and outputGx,y;if the point is not within a reasonable range,we consider it to be noise,and outputGmed.If the window size needs to be enlarged after the first diamond,we will increase the window sizeWx,y,then determine whether the new window size exceeds the maximum valueWmax.If not,we return to “Start”,otherwise we outputGmed.

    2.2.2.Super-resolutionThe upgrade of the GPI system on EAST was completed in 2021 [35].In the upgraded optical system,the high-speed camera has a speed of up to 531645 frames/s and the captured image has a resolution of 128×64 pixel2.When estimating the velocity of the turbulence structure,we extract the image matrix of the structure and its surroundings (approximately twice the structure length in total) and minimize the influence of movement from the outside area.To achieve both high accuracy and smoothness of the image detail,bicubic upsampling super-resolution processing is adopted to extract image data [31].In the bicubic interpolation algorithm,we assume that the interpolation at a certain pointp(x,y) is obtained by weighting the four points around it,which can be written as follows:

    The interpolation kernel function is:

    The interpolation result ofp(x,y) is:

    Assume that pointp(x,y) is located in the original image at the position shown in figure A1,where Q is the matrix of pixel values of the 16 nearest points top. Wxand Wyare calculated by the interpolation kernel function in thexdirection andydirection ofp:

    Rpis the interpolation value of pointpand calculated by multiplying the three matrices of Wx,Q and.In equation (2),ais set as a constant and the default value is -0.5 [31].The specific derivation process is shown in the Appendix.

    For example,we extract the turbulence structure of one GPI frame for super-resolution testing,as shown in figure 4.The reconstructed image with bicubic upsampling exhibits smoother shape edge and a clearer structure compared to the original image.

    Figure 4.(a) The original image,(b) the reconstructed image enlarged ten times with bicubic upsampling.

    2.3.Speed estimation

    2.3.1.Fast Fourier transformFourier transform has been used in various fields of engineering and science.Fourier transform converts images from the spatial domain to the frequency domain,enabling efficient and accurate image analysis.Fast Fourier transform (FFT) [32] is an efficient algorithm of discrete Fourier transform (DFT) [36] and is more suitable for handling large-scale signals or images than DFT.Compared to DFT,FFT incorporates a recursive algorithm based on the divide-and-conquer strategy,which enhances the speed of Fourier transform calculations.The FFT algorithm reduces the computational complexity of DFT fromO(n2) toO(nlogn),and this reduction in complexity allows for faster computation of Fourier transforms.In image processing,FFT plays a significant role in transforming an image from the spatial to the frequency domain.For the discrete 2D Fourier transform of an image with sizeM×N,the equation is as follows:

    whereF(u,v) is the value in the frequency domain,Mis the width of the image,Nis the height of the image,f(x,y) is the value in the spatial domain,iis the imaginary unit,xandyare pixel coordinates,uandvrepresent the horizontal and vertical coordinates in the frequency domain,the value range ofuis 0 ≤u≤M-1,and the value range ofvis 0 ≤v≤N-1.

    2.3.2.2D cross-correlationCross-correlation has many applications in signal and image processing.2D cross-correlation,which computes the correlation coefficient between two images in the spatial domain,is commonly used in searching for similar structures among images,image alignment and image registration.In the SDE method,2D crosscorrelation is used to estimate the relative shift [26],as defined by:

    whereA,B*are 2D signals in frequency-domain obtained using the FFT process,B*is the complex conjugate ofBandF-1denotes the inverse Fourier transform operation.φ(x,y)is the circular 2D cross-correlation function.Since convolution in the time or spatial domain can be expressed by frequency domain multiplication,and correlation is a form of convolution,we express correlation by frequency-domain multiplication [37].The location of the largest correlation can be found at the maximums of φ(x,y) [38].

    Taking two consecutive frames as an example,first,the cross-correlation operation is carried out between the original frame and itself,in order to obtain the relative initial position of the target structure,namedPori.Next,the crosscorrelation between the two frames is calculated to find the point with the maximum coherence value,namedPc.Finally,the relative displacement is given by the Euclidean distance betweenPoriandPc[39].

    The EAST GPI images of discharge #102993 after the pre-processing and adaptive median filtering are shown in figure 5.A bright eddy structure is in the middle of the images,and clear propagation is observed in the time sequence of the images.In order to illustrate the 2D crosscorrelation method,two images from figures 5(i) and (p) are selected to perform the FFT and 2D cross-correlation analysis.The 2D cross-correlation spectrum is converted back to the spatial domain by inverse FFT,and the results are shown in figure 6.Comparing the self-correlation of the first frame with the cross-correlation between frames (i) and (p) in figure 5,the coordinates of the pointsPoriandPccan be fixed,and a clear relative displacement between them is observed,which is consistent with the direction of movement in figure 5.

    2.3.3.Sub-pixel precision processingThe accuracy of the relative displacement derived fromPoriandPcin the previous section is at pixel level.The image recorded by the GPI high-speed camera is 64×128 pixel2on EAST,and the pixel area occupied by the turbulence structure is small.In order to improve the calculation accuracy,a second-order polynomial fitting is used to fit the location with the maximum correlation coefficient and its surrounding area.Consequently,the maximum correlation position can be obtained with sub-pixel accuracy [26].The fitting equation is:

    Figure 5.Time sequences of GPI image after the pre-processing and adaptive median filtering.The corresponding time slices are annotated in panels (a) to (p).

    Figure 6.(a) The cross-correlation operation between the first frame (panel (i) of figure 5) and itself,(b) the cross-correlation calculation between the second frame (panel (p) of figure 5) and the first frame.

    whereai(i=0,1,2,3,4,5) is the fitting coefficient andx,yare coordinates in the cross-correlation spectrum.The position with the maximum correlation coefficient is calculated by the partial differentiation ofxandywith polynomial fitting in equation (9),and the calculated position is given by:

    2.3.4.Velocity calculationFrom the obtainedxmaxandymax,a more accurate relative displacement betweenPoriandPccan be calculated.The coordinates of pointsPoriandPcare(xmax1,ymax1) and (xmax2,ymax2) respectively,and the Euclidean distance between these two points can be expressed as:

    Finally,the speed of movement can be estimated from the relative displacementdand the delay time τ.

    3.Application of SDE algorithm

    3.1.Algorithm test

    In this section,the validity and accuracy of the SDE algorithm are tested.Two types of motion trajectories are used in the test,with a linear motion and a curvilinear motion,as illustrated in figure 7.The frames are set at 137×159 pixel2,which is similar to the GPI rotated image of EAST.The time period of these two movements is set at 5 s,with an interval of 0.1 s between two adjacent frames.For the linear motion,the elliptical structure moves along the diagonal straight line,with a constant shape and direction.For the curvilinear motion,the trajectory is a curve,and the shape and tilted angle of the red structure change gradually frame by frame.The moving velocity calculated by the SDE algorithm is presented in figure 8.For both movements,the calculated velocities agree well with the preset velocities.The velocity of the linear motion has symmetrical distribution,i.e.the red structure accelerates first and then decelerates with the same accelerated speed.The results of both cases demonstrate that the SDE algorithm is valid and its accuracy is very high.

    Figure 7.Motion trajectories.(a) A linear motion,(b) a curvilinear motion.The speed in both cases is first accelerated and then decelerated.

    Figure 8.The structure velocity for the straight and curvilinear motions was calculated by the SDE method.(a) Two-dimensional velocity,(b) velocity in the x-direction,(c) velocity in the y-direction.

    3.2.GPI data analysis

    The SDE algorithm is applied to the GPI data analysis in EAST for discharge #103882.The GPI diagnostic is operated with a frame rate of 531645 Hz and an interval of 1.88μs between two adjacent frames.As a reference,the turbulence propagation velocities in the radial and poloidal directions are calculated using the TDE method,with a velocity point derived from 50 continuous frames,as shown in figure 9.Around the LCFS,the turbulence radial velocityVris mainly directed outwards with a maximum value of about 2 km/s.The turbulence poloidal velocityVθis directed in the ion-diamagnetic drift direction in the SOL and in the electro-diamagnetic drift direction inside the LCFS,with a maximum value of over 3 km/s.During the L-H transition on NSTX,turbulence poloidal velocity was also observed to propagate towards the electron diamagnetic drift direction inside the LCFS [40,41].Note that the event at 8.6256 s is the eruption of an ELM burst.The turbulence velocities are also calculated using the SDE algorithm in the same time period as the TDE method,as illustrated in figure 10.The GPI image is divided into 6×6 sub-frames,with 20×30 pixel2for each sub-frame.In contrast with the TDE method,the radial velocityVrin figure 10(a) exhibits clear outward propagation around the LCFS,which is similar to that in figure 9(a).Note that theVrfrom the SDE method reveals more details of the variation due to its high temporal resolution.The poloidal velocityVθin figure 10(b) reveals the same radial distribution as the TDE method in figure 9(b).Due to the high temporal resolution of the SDE algorithm,the dynamic evolution ofVθis observed within the 550μs period,and the significant changes inVθare consistent with the peaks of the divertor Dαsignal,indicating the correlation between upstream turbulence and downstream particle deposition.During the ELM burst,a significant increase is detected in the Dαsignal,with numerous particles moving outwards across the magnetic field lines.At the same time,high turbulence propagation velocity is observed onVrandVθin figures 9 and 10.In summary,the SDE algorithm reveals similar velocity estimation results to the TDE method,but has much higher temporal resolution.

    4.Summary

    The SDE algorithm is developed and applied to the GPI diagnostic on EAST tokamak.The SDE algorithm includes pre-processing,adaptive median filtering,super-resolution,2D cross-correlation via FFT,and displacement estimation.In these procedures,we remove the background and noise from the image data and preserve the detailed information of the structure.The computation efficiency is high because FFT is introduced into the 2D spatial cross-correlation technology.The temporal resolution of the SDE algorithm is as high as the temporal resolution of the GPI diagnostic,which is a significant advantage for the velocity analysis of fast events.The developed SDE algorithm is tested with a linear motion and a curvilinear motion,and the derived velocities agree with the preset velocities with high accuracy.The new SDE algorithm is applied to the GPI data analysis on EAST,and the calculated propagation velocity of turbulence is consistent with that from the TDE method,but has much higher temporal resolution.Due to its precision and high temporal resolution,the SDE method can estimate the propagation velocity of turbulence in fast transient cases,such as the eruption of ELM and other instabilities.Consequently,the SDE method provides an effective way to analyze the 2D image data with fast variations,which is important for the understanding of edge turbulence in fusion devices.

    Figure 9.The turbulence propagation velocities calculated using the TDE method for discharge #103882.(a) Radial velocity,with positive value directed outside,(b) poloidal velocity,with positive value directed upwards (electron-diamagnetic drift direction).The purple dashed line represents LCFS.

    Figure 10.Turbulence propagation velocities calculated using the SDE algorithm.(a) Radial velocity,(b) poloidal velocity.The temporal variations of a Dα signal are shown in panel (b).The time traces of GPI signal and Dα signal are aligned by some fast events.The purple dashed line represents LCFS.

    Acknowledgments

    This work was supported by the National Magnetic Confinement Fusion Energy R&D Program of China (Nos.2022YFE03030001,2022YFE03020004 and 2022YFE 03050003),National Natural Science Foundation of China(Nos.12275310,11975275,12175277 and 11975271),the Science Foundation of Institute of Plasma Physics,Chinese Academy of Sciences (No.DSJJ-2021-01),the Collaborative Innovation Program of Hefei Science Center,Chinese Academy of Sciences (No.2021HSC-CIP019) and the Users with Excellence Program of Hefei Science Center,Chinese Academy of Sciences (Nos.2021HSC-UE014 and 2021HSCUE012).

    Appendix.Bicubic upsampling

    We developed a coordinate system to demonstrate the calculation process by the kernel function in equation (2).

    Figure A1.Interpolation process in super-resolution processing.

    The coordinates of pointpare (xp,yp),and the interpolation result isRpand the distance from the point to nearby pixels is Δxand Δy.If 0<Δx<1,then 1<Δx+1<2,-1<Δx-1<0,-2<Δx-2<-1.Therefore,from equation (2) we obtain:

    Equation (5) can be derived from the equations listed above.

    猜你喜歡
    顏寧高翔樂天
    清 高翔 錄書七言詩軸
    中國書法(2023年4期)2023-08-28 06:02:08
    High-fidelity quantum sensing of magnon excitations with a single electron spin in quantum dots
    Scaled Preconditioned Splitting Iterative Methods for Solving a Class of Complex Symmetric Linear Systems
    《樂天》
    “清華學(xué)術(shù)女神”在線打假
    東西南北(2019年19期)2019-12-12 06:10:24
    陽關(guān)故人
    飛魔幻A(2019年11期)2019-02-06 03:58:09
    一段苦澀又奇特的成長經(jīng)歷
    花山我的故鄉(xiāng)
    歌海(2016年5期)2016-11-15 09:29:30
    Research survey and review of the effect of Compound Danshen Dripping Pills on the uric acid metabolism of patients with coronary heart disease
    最美科學(xué)家的最美事業(yè)
    做人與處世(2013年6期)2013-06-24 09:38:20
    99久久精品热视频| 大香蕉97超碰在线| 亚洲国产精品国产精品| 国产成人freesex在线| 久久精品人妻少妇| 国产男人的电影天堂91| 在现免费观看毛片| 观看免费一级毛片| 免费观看av网站的网址| 国产精品国产三级国产av玫瑰| 人人妻人人添人人爽欧美一区卜 | 少妇高潮的动态图| 在线观看免费日韩欧美大片 | 国产一区二区在线观看日韩| 国产精品女同一区二区软件| 干丝袜人妻中文字幕| 久久这里有精品视频免费| 全区人妻精品视频| 九九久久精品国产亚洲av麻豆| 少妇裸体淫交视频免费看高清| 国产黄色视频一区二区在线观看| 又黄又爽又刺激的免费视频.| 看十八女毛片水多多多| 国产精品秋霞免费鲁丝片| 国产黄片视频在线免费观看| 青春草国产在线视频| 亚洲自偷自拍三级| 免费观看的影片在线观看| 中文字幕亚洲精品专区| 亚洲av二区三区四区| 久久精品人妻少妇| 深爱激情五月婷婷| 国产精品一二三区在线看| 成人综合一区亚洲| 亚洲av欧美aⅴ国产| 亚洲自偷自拍三级| 水蜜桃什么品种好| 国产精品国产av在线观看| 亚洲成人av在线免费| 日韩伦理黄色片| 人人妻人人添人人爽欧美一区卜 | 欧美精品一区二区大全| 有码 亚洲区| 欧美成人一区二区免费高清观看| 国产大屁股一区二区在线视频| 成人高潮视频无遮挡免费网站| 亚洲美女黄色视频免费看| 狂野欧美激情性xxxx在线观看| 在线免费十八禁| 国产成人a区在线观看| 在线天堂最新版资源| www.色视频.com| 六月丁香七月| 丝瓜视频免费看黄片| 亚洲国产日韩一区二区| 精品久久久久久久久亚洲| 国产精品欧美亚洲77777| 性色avwww在线观看| 国产在线一区二区三区精| 国产成人精品福利久久| 国产黄色免费在线视频| 欧美日韩视频高清一区二区三区二| 精品一区二区三区视频在线| 久久婷婷青草| 亚洲激情五月婷婷啪啪| a级毛片免费高清观看在线播放| 在现免费观看毛片| 赤兔流量卡办理| 日韩欧美 国产精品| 精品久久久久久电影网| 高清欧美精品videossex| 日日啪夜夜爽| 日韩成人伦理影院| 纵有疾风起免费观看全集完整版| 一级a做视频免费观看| 一本一本综合久久| 精品一品国产午夜福利视频| 国产乱人视频| 天美传媒精品一区二区| 纯流量卡能插随身wifi吗| 亚洲,欧美,日韩| 欧美日韩在线观看h| 噜噜噜噜噜久久久久久91| 97精品久久久久久久久久精品| 国产av一区二区精品久久 | 热99国产精品久久久久久7| 涩涩av久久男人的天堂| 边亲边吃奶的免费视频| 中国国产av一级| 免费看av在线观看网站| 大陆偷拍与自拍| 久久久成人免费电影| 欧美xxxx黑人xx丫x性爽| 亚洲精品一区蜜桃| 三级国产精品欧美在线观看| 国产精品免费大片| 国产无遮挡羞羞视频在线观看| 成人综合一区亚洲| 国产亚洲5aaaaa淫片| 一级二级三级毛片免费看| 亚洲精品456在线播放app| 欧美激情极品国产一区二区三区 | 精品一区二区三卡| 91久久精品电影网| 18禁在线播放成人免费| 麻豆精品久久久久久蜜桃| 日日摸夜夜添夜夜添av毛片| 欧美bdsm另类| 插阴视频在线观看视频| 久热这里只有精品99| 中文天堂在线官网| 欧美另类一区| 晚上一个人看的免费电影| 亚洲国产日韩一区二区| 国产欧美另类精品又又久久亚洲欧美| 午夜激情福利司机影院| 亚洲,欧美,日韩| 又大又黄又爽视频免费| 一本一本综合久久| 欧美+日韩+精品| 美女主播在线视频| 女性被躁到高潮视频| 久久精品熟女亚洲av麻豆精品| 午夜福利高清视频| 丰满迷人的少妇在线观看| 欧美3d第一页| 黄色日韩在线| 日韩av免费高清视频| 亚洲精品第二区| 日韩国内少妇激情av| 精华霜和精华液先用哪个| 最近中文字幕高清免费大全6| 国产欧美日韩一区二区三区在线 | 久久久久国产网址| 亚洲精品色激情综合| 丝袜脚勾引网站| 另类亚洲欧美激情| 26uuu在线亚洲综合色| 国产免费视频播放在线视频| 日韩一区二区视频免费看| 午夜日本视频在线| 亚洲国产精品一区三区| 亚洲精品456在线播放app| av天堂中文字幕网| 一级片'在线观看视频| 国产精品一区www在线观看| 久热这里只有精品99| 少妇人妻一区二区三区视频| 王馨瑶露胸无遮挡在线观看| 亚洲av成人精品一区久久| 身体一侧抽搐| 日韩成人av中文字幕在线观看| 国产在线一区二区三区精| 亚洲欧美一区二区三区黑人 | 老司机影院成人| 日韩在线高清观看一区二区三区| 日韩制服骚丝袜av| 99久久综合免费| 又黄又爽又刺激的免费视频.| 国产乱人视频| 欧美日韩国产mv在线观看视频 | 舔av片在线| 极品教师在线视频| 在线看a的网站| 91久久精品国产一区二区三区| 成人二区视频| 99精国产麻豆久久婷婷| 免费大片18禁| 精品一区二区三卡| 免费观看a级毛片全部| 国产免费视频播放在线视频| 另类亚洲欧美激情| 丰满人妻一区二区三区视频av| 久久久久精品性色| 亚洲无线观看免费| 久久人人爽人人爽人人片va| 亚洲av不卡在线观看| 国产高清不卡午夜福利| 舔av片在线| 内地一区二区视频在线| 国产爱豆传媒在线观看| 亚洲国产最新在线播放| 欧美日本视频| 搡女人真爽免费视频火全软件| 免费观看的影片在线观看| 亚洲欧美成人精品一区二区| 我的女老师完整版在线观看| 欧美 日韩 精品 国产| 国产有黄有色有爽视频| 午夜免费鲁丝| 国产一区有黄有色的免费视频| 男女啪啪激烈高潮av片| 欧美成人一区二区免费高清观看| 99国产精品免费福利视频| 亚洲一级一片aⅴ在线观看| 国产精品一区www在线观看| 人妻制服诱惑在线中文字幕| 男女边吃奶边做爰视频| 涩涩av久久男人的天堂| 亚洲成人中文字幕在线播放| 在现免费观看毛片| 人妻制服诱惑在线中文字幕| 精品一区在线观看国产| 最近最新中文字幕大全电影3| 这个男人来自地球电影免费观看 | 青青草视频在线视频观看| 高清在线视频一区二区三区| 日本色播在线视频| 国产精品99久久99久久久不卡 | 国产成人一区二区在线| 最黄视频免费看| 2021少妇久久久久久久久久久| kizo精华| 日韩av免费高清视频| 中文字幕制服av| 亚洲欧美一区二区三区黑人 | 中文字幕制服av| 成人毛片a级毛片在线播放| 亚洲精品久久午夜乱码| 日韩欧美一区视频在线观看 | 91精品伊人久久大香线蕉| 99久久精品一区二区三区| 一区二区av电影网| 国产精品久久久久久精品古装| 九九在线视频观看精品| 我要看日韩黄色一级片| 精品久久久噜噜| 精品人妻视频免费看| 国产在视频线精品| 国产白丝娇喘喷水9色精品| 亚洲内射少妇av| 欧美成人午夜免费资源| 日韩一本色道免费dvd| 卡戴珊不雅视频在线播放| 日韩电影二区| 五月伊人婷婷丁香| 亚洲精品成人av观看孕妇| 黄色日韩在线| 一区二区三区四区激情视频| 午夜免费男女啪啪视频观看| 亚洲伊人久久精品综合| 99热全是精品| 久久热精品热| 日本一二三区视频观看| www.av在线官网国产| 2018国产大陆天天弄谢| av黄色大香蕉| 精品国产乱码久久久久久小说| 久久久久精品性色| 日本-黄色视频高清免费观看| 欧美一区二区亚洲| a级一级毛片免费在线观看| 色婷婷av一区二区三区视频| 狠狠精品人妻久久久久久综合| 亚洲av日韩在线播放| 嫩草影院新地址| 大又大粗又爽又黄少妇毛片口| 亚洲欧美成人综合另类久久久| 久久久久性生活片| 交换朋友夫妻互换小说| 一级毛片我不卡| 国产毛片在线视频| 免费av中文字幕在线| 直男gayav资源| 汤姆久久久久久久影院中文字幕| 精品人妻熟女av久视频| 黄色日韩在线| 91精品国产国语对白视频| 亚洲欧美成人精品一区二区| 国产精品欧美亚洲77777| 新久久久久国产一级毛片| 国产精品偷伦视频观看了| 大香蕉97超碰在线| 高清午夜精品一区二区三区| 国产av国产精品国产| 久热这里只有精品99| av免费在线看不卡| 蜜桃久久精品国产亚洲av| 亚洲精品乱码久久久久久按摩| 性色avwww在线观看| 久久久久国产网址| 搡老乐熟女国产| tube8黄色片| 国产一区二区在线观看日韩| 国产v大片淫在线免费观看| 成人毛片a级毛片在线播放| 97热精品久久久久久| 亚洲国产精品999| 五月玫瑰六月丁香| 99久久人妻综合| 日韩人妻高清精品专区| 简卡轻食公司| 成人毛片a级毛片在线播放| 亚洲综合精品二区| 女人久久www免费人成看片| 国产视频首页在线观看| 国产在线男女| 有码 亚洲区| 亚洲av成人精品一区久久| 99久久精品热视频| 亚洲av在线观看美女高潮| 91午夜精品亚洲一区二区三区| 亚洲国产最新在线播放| 亚洲成人av在线免费| 亚洲欧美日韩卡通动漫| 国产探花极品一区二区| 校园人妻丝袜中文字幕| 联通29元200g的流量卡| 久久精品国产自在天天线| 啦啦啦啦在线视频资源| av在线观看视频网站免费| 麻豆成人午夜福利视频| 国产成人精品久久久久久| 一区二区三区免费毛片| 久久久久久久久久人人人人人人| 99久久精品热视频| 丝瓜视频免费看黄片| 哪个播放器可以免费观看大片| 国产成人aa在线观看| 亚洲欧美一区二区三区黑人 | 亚洲人成网站高清观看| 最近最新中文字幕免费大全7| 欧美 日韩 精品 国产| 国产精品福利在线免费观看| 欧美丝袜亚洲另类| 国产黄色免费在线视频| 亚洲精品成人av观看孕妇| 边亲边吃奶的免费视频| 乱系列少妇在线播放| av福利片在线观看| 能在线免费看毛片的网站| 亚洲无线观看免费| 观看美女的网站| 在线观看免费视频网站a站| 国产精品人妻久久久影院| 伊人久久国产一区二区| 精品久久久久久电影网| 涩涩av久久男人的天堂| 亚洲一级一片aⅴ在线观看| 亚洲国产精品国产精品| 久久久久久久久大av| 狠狠精品人妻久久久久久综合| 一区在线观看完整版| 一级黄片播放器| 国产免费福利视频在线观看| 2021少妇久久久久久久久久久| 熟女电影av网| 国产精品熟女久久久久浪| av天堂中文字幕网| 一级毛片我不卡| 老司机影院毛片| 啦啦啦视频在线资源免费观看| 国产一区二区三区综合在线观看 | 色网站视频免费| 亚洲欧美日韩卡通动漫| 三级经典国产精品| 国产成人a∨麻豆精品| 大香蕉久久网| 91精品国产国语对白视频| 黄色一级大片看看| 一二三四中文在线观看免费高清| 波野结衣二区三区在线| 国模一区二区三区四区视频| 亚洲色图av天堂| 日本爱情动作片www.在线观看| 亚洲国产精品999| 成人漫画全彩无遮挡| 小蜜桃在线观看免费完整版高清| 成人特级av手机在线观看| 国产av一区二区精品久久 | 国产一区有黄有色的免费视频| 我的女老师完整版在线观看| 日本欧美视频一区| 国产精品免费大片| 日韩在线高清观看一区二区三区| 少妇裸体淫交视频免费看高清| 久久久久久久久久成人| 最近中文字幕高清免费大全6| 不卡视频在线观看欧美| 国精品久久久久久国模美| 黄色一级大片看看| 97热精品久久久久久| 最近中文字幕2019免费版| 一本久久精品| 国产av码专区亚洲av| 精品少妇久久久久久888优播| 国产成人一区二区在线| 亚洲欧美日韩卡通动漫| 色网站视频免费| 亚洲va在线va天堂va国产| 免费看光身美女| 高清不卡的av网站| 亚洲成色77777| 精品久久久久久久末码| 亚洲av二区三区四区| 欧美成人一区二区免费高清观看| 午夜福利在线观看免费完整高清在| 日本欧美国产在线视频| 日韩欧美 国产精品| 狂野欧美白嫩少妇大欣赏| 男女啪啪激烈高潮av片| 妹子高潮喷水视频| 国内揄拍国产精品人妻在线| 观看免费一级毛片| 国产精品三级大全| kizo精华| 久久6这里有精品| 九草在线视频观看| 国精品久久久久久国模美| 久久久久久久久久久免费av| 国产在视频线精品| 男的添女的下面高潮视频| 少妇 在线观看| 日本与韩国留学比较| 妹子高潮喷水视频| 一区二区av电影网| 欧美一级a爱片免费观看看| 一级a做视频免费观看| 国产又色又爽无遮挡免| 亚洲av福利一区| 女的被弄到高潮叫床怎么办| 欧美bdsm另类| 亚洲av在线观看美女高潮| 中文欧美无线码| 熟女人妻精品中文字幕| 亚洲欧美精品专区久久| 亚洲精品国产av蜜桃| 亚洲欧美成人综合另类久久久| 精品人妻一区二区三区麻豆| 精品久久久久久久久av| 中文资源天堂在线| 国产成人a∨麻豆精品| 在线天堂最新版资源| 国产片特级美女逼逼视频| 亚洲av日韩在线播放| 久久久久国产精品人妻一区二区| av播播在线观看一区| 亚洲精品国产av蜜桃| 亚洲av男天堂| 黄色欧美视频在线观看| 在线观看国产h片| 亚洲国产日韩一区二区| 日韩免费高清中文字幕av| 寂寞人妻少妇视频99o| 久久久久精品性色| 91aial.com中文字幕在线观看| 久久久久久人妻| 最近最新中文字幕免费大全7| 国产淫语在线视频| 在线观看国产h片| 黄色一级大片看看| av不卡在线播放| 男女啪啪激烈高潮av片| 欧美成人精品欧美一级黄| 亚洲国产精品专区欧美| 国产精品久久久久久久电影| 乱系列少妇在线播放| 街头女战士在线观看网站| 极品教师在线视频| 亚洲一级一片aⅴ在线观看| 午夜免费观看性视频| 青青草视频在线视频观看| 国产黄片美女视频| 麻豆成人午夜福利视频| 日本av免费视频播放| 亚洲精品乱码久久久v下载方式| 成人国产av品久久久| 干丝袜人妻中文字幕| 亚洲av成人精品一二三区| 免费少妇av软件| 久久久精品94久久精品| 国产精品一区二区三区四区免费观看| 午夜福利影视在线免费观看| 22中文网久久字幕| 成人二区视频| 国产精品国产三级国产av玫瑰| 亚洲人成网站在线观看播放| 亚洲aⅴ乱码一区二区在线播放| av.在线天堂| 亚洲欧美一区二区三区国产| 国产一区二区在线观看日韩| 国产免费又黄又爽又色| 欧美日韩视频高清一区二区三区二| av福利片在线观看| 最近2019中文字幕mv第一页| 观看av在线不卡| 国产色爽女视频免费观看| av一本久久久久| 一个人看视频在线观看www免费| 国产精品久久久久久av不卡| 国产爱豆传媒在线观看| 亚洲av不卡在线观看| 国产一区亚洲一区在线观看| 国产伦在线观看视频一区| 午夜福利视频精品| 男的添女的下面高潮视频| 美女视频免费永久观看网站| 99久久中文字幕三级久久日本| 国产淫语在线视频| 国产高清有码在线观看视频| 男女免费视频国产| 成年免费大片在线观看| 日本wwww免费看| 精品久久久噜噜| 国语对白做爰xxxⅹ性视频网站| 纵有疾风起免费观看全集完整版| 国产精品久久久久久久久免| 一本色道久久久久久精品综合| 亚洲精品一二三| 亚洲高清免费不卡视频| 成人毛片60女人毛片免费| 亚洲精品456在线播放app| 一本色道久久久久久精品综合| 欧美日本视频| 欧美日韩国产mv在线观看视频 | 免费看av在线观看网站| 国产精品国产三级国产av玫瑰| 日韩在线高清观看一区二区三区| av国产久精品久网站免费入址| 在线看a的网站| 一级av片app| 色吧在线观看| 欧美高清成人免费视频www| 18禁在线无遮挡免费观看视频| 日韩一本色道免费dvd| 免费黄网站久久成人精品| 一级a做视频免费观看| 国产成人精品福利久久| 国产在视频线精品| 亚洲欧美精品专区久久| 亚洲人成网站高清观看| 久久久久久久精品精品| 国产免费一级a男人的天堂| 蜜桃亚洲精品一区二区三区| 久久久久视频综合| 精品少妇黑人巨大在线播放| 一本一本综合久久| 日韩中字成人| 大又大粗又爽又黄少妇毛片口| av国产精品久久久久影院| 国产免费一级a男人的天堂| 特大巨黑吊av在线直播| 亚洲欧洲日产国产| 久久女婷五月综合色啪小说| 99热6这里只有精品| 精品国产三级普通话版| 亚洲av中文字字幕乱码综合| 一区二区三区免费毛片| 免费人妻精品一区二区三区视频| videos熟女内射| 亚洲精品国产成人久久av| 国产高清不卡午夜福利| 天堂俺去俺来也www色官网| 亚洲国产精品专区欧美| 亚洲无线观看免费| 日本欧美国产在线视频| 少妇被粗大猛烈的视频| 一级毛片aaaaaa免费看小| 91精品伊人久久大香线蕉| 91午夜精品亚洲一区二区三区| 久久精品国产鲁丝片午夜精品| 我要看日韩黄色一级片| 精品熟女少妇av免费看| 女的被弄到高潮叫床怎么办| 国国产精品蜜臀av免费| 亚洲经典国产精华液单| 麻豆精品久久久久久蜜桃| 爱豆传媒免费全集在线观看| 免费看光身美女| 丝袜脚勾引网站| 九九久久精品国产亚洲av麻豆| 国产精品不卡视频一区二区| 2021少妇久久久久久久久久久| av国产免费在线观看| 国产伦精品一区二区三区四那| 久久午夜福利片| 精品熟女少妇av免费看| 少妇熟女欧美另类| 成人国产麻豆网| 国产高清三级在线| 国国产精品蜜臀av免费| 看十八女毛片水多多多| 水蜜桃什么品种好| 视频区图区小说| 久久这里有精品视频免费| 日日摸夜夜添夜夜爱| 午夜福利在线观看免费完整高清在| 99热这里只有是精品在线观看| 国产成人91sexporn| 午夜免费男女啪啪视频观看| 六月丁香七月| 亚洲av欧美aⅴ国产| 欧美极品一区二区三区四区| 亚洲内射少妇av| 国产在线免费精品| 免费观看的影片在线观看| 蜜臀久久99精品久久宅男| 久久精品国产亚洲av天美| 有码 亚洲区| 日韩,欧美,国产一区二区三区| 女性被躁到高潮视频| 国产午夜精品久久久久久一区二区三区| 丝袜喷水一区| 在线观看国产h片| 波野结衣二区三区在线| 亚洲欧美中文字幕日韩二区| 日韩伦理黄色片| 人妻 亚洲 视频| 久久97久久精品| 国产爱豆传媒在线观看| 亚洲精品乱码久久久v下载方式| a级毛片免费高清观看在线播放| 国产精品一区二区在线观看99| 秋霞在线观看毛片| 免费不卡的大黄色大毛片视频在线观看|