• <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热全是精品| 啦啦啦在线观看免费高清www| 在线观看国产h片| 免费日韩欧美在线观看| 男人舔女人的私密视频| 国产探花极品一区二区| 18在线观看网站| 色婷婷av一区二区三区视频| 好男人视频免费观看在线| 国产亚洲一区二区精品| 亚洲,欧美精品.| 欧美日韩一区二区视频在线观看视频在线| 色吧在线观看| 亚洲欧美日韩卡通动漫| 下体分泌物呈黄色| 成年动漫av网址| 免费高清在线观看日韩| 制服诱惑二区| 成人18禁高潮啪啪吃奶动态图| 在现免费观看毛片| 男人添女人高潮全过程视频| 熟妇人妻不卡中文字幕| 波野结衣二区三区在线| 欧美人与性动交α欧美精品济南到 | 国产不卡av网站在线观看| 精品熟女少妇av免费看| 色吧在线观看| 精品卡一卡二卡四卡免费| 久久久久久人人人人人| 久久久久精品人妻al黑| 又大又黄又爽视频免费| 卡戴珊不雅视频在线播放| 亚洲精品456在线播放app| 国产精品人妻久久久久久| 精品人妻偷拍中文字幕| 伦理电影免费视频| 侵犯人妻中文字幕一二三四区| 久久狼人影院| 最近最新中文字幕大全免费视频 | 久久久久久久亚洲中文字幕| 熟女电影av网| 99久久中文字幕三级久久日本| 亚洲高清免费不卡视频| 亚洲欧美色中文字幕在线| 精品国产露脸久久av麻豆| 菩萨蛮人人尽说江南好唐韦庄| 91精品国产国语对白视频| 国产精品蜜桃在线观看| av有码第一页| 久久精品久久久久久噜噜老黄| 国精品久久久久久国模美| 国产精品久久久久久久久免| 国产免费现黄频在线看| 日韩制服骚丝袜av| 精品国产国语对白av| 黑人巨大精品欧美一区二区蜜桃 | 久久久久网色| 男女无遮挡免费网站观看| 国产免费一级a男人的天堂| 久久久久视频综合| 亚洲欧美一区二区三区国产| 91精品三级在线观看| 中文欧美无线码| 两个人看的免费小视频| 丰满迷人的少妇在线观看| 亚洲av福利一区| 免费av中文字幕在线| 欧美精品国产亚洲| 在线观看三级黄色| 成人影院久久| 国产免费视频播放在线视频| 亚洲国产最新在线播放| 亚洲欧美日韩另类电影网站| 国产精品久久久久久精品古装| 久久99热6这里只有精品| 国产精品99久久99久久久不卡 | 97人妻天天添夜夜摸| 国产精品久久久久成人av| 在现免费观看毛片| 欧美老熟妇乱子伦牲交| 波野结衣二区三区在线| 蜜臀久久99精品久久宅男| 国产精品一区二区在线不卡| 久久久久精品性色| 丝袜在线中文字幕| 国产成人精品在线电影| 国产极品粉嫩免费观看在线| 99久久综合免费| 麻豆精品久久久久久蜜桃| 国产免费又黄又爽又色| 深夜精品福利| 日韩av不卡免费在线播放| 精品久久久久久电影网| 夜夜骑夜夜射夜夜干| 999精品在线视频| 日本黄色日本黄色录像| 只有这里有精品99| 欧美少妇被猛烈插入视频| 久久精品夜色国产| 精品一品国产午夜福利视频| 人体艺术视频欧美日本| 五月开心婷婷网| 亚洲第一av免费看| 97在线人人人人妻| av在线app专区| 三上悠亚av全集在线观看| 国产日韩欧美在线精品| 亚洲人成77777在线视频| 久久女婷五月综合色啪小说| 赤兔流量卡办理| 自线自在国产av| 国产精品久久久久久久久免| 三级国产精品片| a级毛片黄视频| 女的被弄到高潮叫床怎么办| 婷婷色av中文字幕| 韩国高清视频一区二区三区| 亚洲欧美日韩另类电影网站| 亚洲国产色片| 亚洲四区av| 三上悠亚av全集在线观看| 91成人精品电影| 中文字幕亚洲精品专区| 久久久国产精品麻豆| 最新的欧美精品一区二区| 亚洲成国产人片在线观看| 久久99热这里只频精品6学生| 欧美日韩av久久| 国产深夜福利视频在线观看| 久久久久久久久久久久大奶| 国产一区二区在线观看av| 热99国产精品久久久久久7| 亚洲图色成人| 人人妻人人爽人人添夜夜欢视频| 少妇的逼水好多| 国产精品 国内视频| 又粗又硬又长又爽又黄的视频| 九色亚洲精品在线播放| 美女中出高潮动态图| 欧美人与善性xxx| 久久这里只有精品19| 色视频在线一区二区三区| 亚洲精品日本国产第一区| 极品少妇高潮喷水抽搐| 99国产综合亚洲精品| 一级a做视频免费观看| a级毛片黄视频| 少妇精品久久久久久久| 中文欧美无线码| 黑人猛操日本美女一级片| 中文字幕av电影在线播放| 日日啪夜夜爽| 伦理电影大哥的女人| av又黄又爽大尺度在线免费看| 亚洲国产最新在线播放| 男人舔女人的私密视频| 久久99一区二区三区| 久久久久久伊人网av| 午夜免费鲁丝| 伊人久久国产一区二区| 色94色欧美一区二区| 欧美精品国产亚洲| 国产精品一区www在线观看| 97在线人人人人妻| 91aial.com中文字幕在线观看| 97在线视频观看| 亚洲成人一二三区av| 欧美+日韩+精品| 日本-黄色视频高清免费观看| 夫妻午夜视频| 国产老妇伦熟女老妇高清| 香蕉精品网在线| 久久鲁丝午夜福利片| 91国产中文字幕| 美女内射精品一级片tv| 黑人巨大精品欧美一区二区蜜桃 | 激情五月婷婷亚洲| a级毛色黄片| 久久精品人人爽人人爽视色| 色吧在线观看| 久久人人爽av亚洲精品天堂| 国产成人精品一,二区| 人体艺术视频欧美日本| 在线观看美女被高潮喷水网站| 少妇人妻精品综合一区二区| 熟女av电影| 一区二区av电影网| 精品久久久精品久久久| 久久精品国产亚洲av涩爱| 亚洲国产色片| 久久国产精品男人的天堂亚洲 | 熟女人妻精品中文字幕| 免费看光身美女| 色吧在线观看| 国产老妇伦熟女老妇高清| 国产午夜精品一二区理论片| 国产有黄有色有爽视频| 国产极品粉嫩免费观看在线| 丝袜人妻中文字幕| 国产成人一区二区在线| 乱人伦中国视频| 久久99热6这里只有精品| 国语对白做爰xxxⅹ性视频网站| 亚洲,一卡二卡三卡| 丝袜在线中文字幕| 日韩欧美一区视频在线观看| 久久久久久久久久久免费av| 男女啪啪激烈高潮av片| 亚洲精品视频女| 亚洲av男天堂| 一级片'在线观看视频| 亚洲欧洲国产日韩| 国产综合精华液| 中文乱码字字幕精品一区二区三区| 99九九在线精品视频| 人妻人人澡人人爽人人| 精品久久久精品久久久| 秋霞在线观看毛片| 国产av国产精品国产| 日本av手机在线免费观看| 精品一区二区免费观看| 亚洲av.av天堂| 国产精品成人在线| 午夜av观看不卡| 精品福利永久在线观看| 国产熟女午夜一区二区三区| 日韩伦理黄色片| 免费观看a级毛片全部| 一级,二级,三级黄色视频| 国产极品粉嫩免费观看在线| 黑人高潮一二区| √禁漫天堂资源中文www| 久久 成人 亚洲| 999精品在线视频| 亚洲在久久综合| 国产精品久久久久久久电影| 日韩人妻精品一区2区三区| 交换朋友夫妻互换小说| 一本一本久久a久久精品综合妖精 国产伦在线观看视频一区 | 寂寞人妻少妇视频99o| 在线观看免费日韩欧美大片| 国产一区二区在线观看av| 黄色毛片三级朝国网站| 侵犯人妻中文字幕一二三四区| 国产免费又黄又爽又色| 只有这里有精品99| 午夜福利在线免费观看网站| 国内毛片毛片毛片毛片毛片| 日本一区二区免费在线视频| 久久久久精品人妻al黑| 久久中文字幕人妻熟女| 国产精品成人在线| 91麻豆av在线| 又黄又粗又硬又大视频| 无限看片的www在线观看| 老司机靠b影院| 久久精品熟女亚洲av麻豆精品| 老汉色av国产亚洲站长工具| 电影成人av| 精品欧美一区二区三区在线| 欧美乱码精品一区二区三区| 国产99白浆流出| 一区二区日韩欧美中文字幕| 一二三四社区在线视频社区8| 久9热在线精品视频| 国产精品免费大片| 99国产精品99久久久久| 欧美 亚洲 国产 日韩一| 久久久精品区二区三区| 国产真人三级小视频在线观看| 国产男女超爽视频在线观看| 1024视频免费在线观看| 免费高清在线观看日韩| 国产亚洲精品久久久久5区| 午夜福利在线免费观看网站| 成人亚洲精品一区在线观看| 久久久国产精品麻豆| а√天堂www在线а√下载 | 亚洲av成人av| 国产成人精品久久二区二区91| 大片电影免费在线观看免费| 久久青草综合色| 久久久久久免费高清国产稀缺| 久久精品国产清高在天天线| 国产精品美女特级片免费视频播放器 | 日韩一卡2卡3卡4卡2021年| 黑人猛操日本美女一级片| 亚洲av成人一区二区三| 午夜影院日韩av| 欧美日本中文国产一区发布| 黑丝袜美女国产一区| 一级作爱视频免费观看| 99热国产这里只有精品6| 制服人妻中文乱码| 欧美不卡视频在线免费观看 | 精品无人区乱码1区二区| 亚洲色图综合在线观看| 久久久久国产精品人妻aⅴ院 | 高清毛片免费观看视频网站 | 中文字幕制服av| 亚洲成人国产一区在线观看| 天天操日日干夜夜撸| 国产成人免费观看mmmm| 色婷婷久久久亚洲欧美| 老司机靠b影院| 精品乱码久久久久久99久播| 久久99一区二区三区| 国产高清国产精品国产三级| 亚洲国产精品合色在线| 99久久人妻综合| 大陆偷拍与自拍| 亚洲色图综合在线观看| 新久久久久国产一级毛片| 欧美精品一区二区免费开放| 在线观看舔阴道视频| 久久国产乱子伦精品免费另类| 一级作爱视频免费观看| 欧美不卡视频在线免费观看 | 国产97色在线日韩免费| 手机成人av网站| 精品免费久久久久久久清纯 | a级片在线免费高清观看视频| 亚洲人成电影观看| 国产精品永久免费网站| 亚洲片人在线观看| 色精品久久人妻99蜜桃| 色婷婷久久久亚洲欧美| 免费在线观看完整版高清| 欧美日韩成人在线一区二区| 亚洲中文日韩欧美视频| 精品少妇一区二区三区视频日本电影| 人成视频在线观看免费观看| 精品午夜福利视频在线观看一区| 亚洲精品中文字幕在线视频| 男女之事视频高清在线观看| 亚洲成a人片在线一区二区| 天天躁夜夜躁狠狠躁躁| 久久天躁狠狠躁夜夜2o2o| 手机成人av网站| 免费高清在线观看日韩| 一级毛片精品| 午夜福利乱码中文字幕| 男女高潮啪啪啪动态图| 十八禁高潮呻吟视频| 18禁国产床啪视频网站| 无人区码免费观看不卡| 少妇 在线观看| 免费黄频网站在线观看国产| 欧美日韩视频精品一区| av电影中文网址| 搡老岳熟女国产| 久久精品国产99精品国产亚洲性色 | 亚洲熟女毛片儿| 中亚洲国语对白在线视频| 免费一级毛片在线播放高清视频 | 欧美中文综合在线视频| 午夜福利乱码中文字幕| 欧美成狂野欧美在线观看| 麻豆成人av在线观看| 免费高清在线观看日韩| 精品午夜福利视频在线观看一区| 日韩欧美国产一区二区入口| 大型av网站在线播放| 亚洲国产精品合色在线| 日韩熟女老妇一区二区性免费视频| 99精品在免费线老司机午夜| 在线观看免费视频网站a站| 99久久精品国产亚洲精品| 国产激情欧美一区二区| 飞空精品影院首页| 一级a爱片免费观看的视频| 欧美亚洲 丝袜 人妻 在线| 成人国语在线视频| 一区二区三区国产精品乱码| av在线播放免费不卡| 亚洲av日韩精品久久久久久密| 亚洲第一欧美日韩一区二区三区| 亚洲自偷自拍图片 自拍| 午夜免费成人在线视频| 91精品三级在线观看| 狂野欧美激情性xxxx| 欧美乱码精品一区二区三区| 国产又爽黄色视频| 午夜精品久久久久久毛片777| 夜夜躁狠狠躁天天躁| 国产成人欧美| 午夜福利一区二区在线看| 精品无人区乱码1区二区| 香蕉久久夜色| 精品一区二区三卡| 亚洲第一欧美日韩一区二区三区| www.自偷自拍.com| 18在线观看网站| 久久香蕉精品热| 黑丝袜美女国产一区| 亚洲专区字幕在线| 亚洲精品国产精品久久久不卡| 亚洲 欧美一区二区三区| 午夜福利乱码中文字幕| 精品第一国产精品| av在线播放免费不卡| 777米奇影视久久| 超碰成人久久| 亚洲五月天丁香| 91精品国产国语对白视频| 一进一出好大好爽视频| 首页视频小说图片口味搜索| 日韩 欧美 亚洲 中文字幕| 免费观看a级毛片全部| 免费观看精品视频网站| 天天躁狠狠躁夜夜躁狠狠躁| 国产成人欧美在线观看 | 亚洲,欧美精品.| 久久国产精品人妻蜜桃| 久久精品熟女亚洲av麻豆精品| svipshipincom国产片| 新久久久久国产一级毛片| 久久草成人影院| 国产精品综合久久久久久久免费 | av不卡在线播放| 国产亚洲精品久久久久5区| 亚洲色图综合在线观看| 精品国产美女av久久久久小说| xxx96com| 精品一区二区三卡| 久久人妻av系列| 亚洲欧美色中文字幕在线| 午夜精品久久久久久毛片777| av免费在线观看网站| 一级,二级,三级黄色视频| 国产在视频线精品| 国产成人系列免费观看| 久久久久久久久久久久大奶| 久久久久国内视频| 正在播放国产对白刺激| 国产精品综合久久久久久久免费 | 免费黄频网站在线观看国产| 亚洲av日韩精品久久久久久密| av天堂久久9| 999久久久精品免费观看国产| 大香蕉久久网| 成人影院久久| 好男人电影高清在线观看| 一区二区三区精品91| 欧美精品啪啪一区二区三区| 欧美激情 高清一区二区三区| www.精华液| 丝瓜视频免费看黄片| 精品国产一区二区三区四区第35| 狠狠婷婷综合久久久久久88av| 成人影院久久| 欧美在线黄色| 欧美精品人与动牲交sv欧美| av国产精品久久久久影院| 国产精品国产av在线观看| 久久久久视频综合| 十八禁高潮呻吟视频| 国产熟女午夜一区二区三区| 久热爱精品视频在线9| 久久国产乱子伦精品免费另类| 久久九九热精品免费| 老熟妇乱子伦视频在线观看| 久久精品人人爽人人爽视色| 最新美女视频免费是黄的| 午夜福利欧美成人| 欧美丝袜亚洲另类 | 精品少妇一区二区三区视频日本电影| 最近最新中文字幕大全免费视频| 黑人欧美特级aaaaaa片| 午夜福利欧美成人| 精品欧美一区二区三区在线| 每晚都被弄得嗷嗷叫到高潮| 国产精品久久视频播放| 精品久久久精品久久久| 国产三级黄色录像| 国产精品一区二区免费欧美| 国产xxxxx性猛交| av片东京热男人的天堂| 美女视频免费永久观看网站| 天天操日日干夜夜撸| 天天躁日日躁夜夜躁夜夜| 午夜福利欧美成人| 在线观看日韩欧美| 欧美精品亚洲一区二区| 色94色欧美一区二区| 99精品在免费线老司机午夜| 午夜福利视频在线观看免费| 亚洲人成电影观看| www日本在线高清视频| 黑丝袜美女国产一区| 91精品国产国语对白视频| 80岁老熟妇乱子伦牲交| 美女扒开内裤让男人捅视频| 欧美日韩亚洲高清精品| 18在线观看网站| 国产在视频线精品| 99国产精品免费福利视频| 法律面前人人平等表现在哪些方面| 男人舔女人的私密视频| 青草久久国产| 搡老熟女国产l中国老女人| 宅男免费午夜| 中出人妻视频一区二区| 欧美日韩亚洲国产一区二区在线观看 | 国内毛片毛片毛片毛片毛片| 看黄色毛片网站| 国产精品久久久av美女十八| 日韩成人在线观看一区二区三区| 亚洲av熟女| 午夜91福利影院| 国产午夜精品久久久久久| 欧美精品一区二区免费开放| 免费在线观看视频国产中文字幕亚洲| 女人高潮潮喷娇喘18禁视频| av网站免费在线观看视频| 人成视频在线观看免费观看| 91老司机精品| 一a级毛片在线观看| 麻豆国产av国片精品| 亚洲一区中文字幕在线| 国产不卡av网站在线观看| 免费看十八禁软件| 亚洲国产欧美网| 久久久国产精品麻豆| svipshipincom国产片| 涩涩av久久男人的天堂| 这个男人来自地球电影免费观看| 亚洲 欧美一区二区三区| 80岁老熟妇乱子伦牲交| 精品国产一区二区久久| 18禁观看日本| 少妇裸体淫交视频免费看高清 | 亚洲欧洲精品一区二区精品久久久| 极品少妇高潮喷水抽搐| 一区二区三区精品91| 美女国产高潮福利片在线看| 精品人妻熟女毛片av久久网站| www.熟女人妻精品国产| 亚洲国产欧美日韩在线播放| 精品国产一区二区三区久久久樱花| 免费观看精品视频网站| 99re在线观看精品视频| 三级毛片av免费| 动漫黄色视频在线观看| 精品国产一区二区三区久久久樱花| 50天的宝宝边吃奶边哭怎么回事| 国产精品久久久人人做人人爽| tocl精华| 一本大道久久a久久精品| 免费少妇av软件| 国产麻豆69| 午夜两性在线视频| 免费女性裸体啪啪无遮挡网站| 国产精品一区二区在线不卡| 99久久人妻综合| 久9热在线精品视频| 国产精品国产高清国产av | 国产亚洲欧美在线一区二区| 人人妻人人澡人人看| 国产精品国产高清国产av | 欧美成狂野欧美在线观看| 一本综合久久免费| 久久精品国产清高在天天线| 操美女的视频在线观看| 丰满饥渴人妻一区二区三| 久久久久精品人妻al黑| 成熟少妇高潮喷水视频| 欧美日韩黄片免| 狠狠狠狠99中文字幕| 伦理电影免费视频| 久久精品国产a三级三级三级| 最近最新免费中文字幕在线| 精品国产国语对白av| 日韩欧美一区视频在线观看| 国产精品久久久久久人妻精品电影| 欧美中文综合在线视频| 人人妻人人添人人爽欧美一区卜| 中文字幕最新亚洲高清| 飞空精品影院首页| 久久亚洲真实| 日韩免费高清中文字幕av| 99在线人妻在线中文字幕 | 免费日韩欧美在线观看| 老司机影院毛片| 丰满人妻熟妇乱又伦精品不卡| 久久狼人影院| 欧美日韩视频精品一区| a级毛片黄视频| 久久久久精品人妻al黑| 国产成人av激情在线播放| 一级a爱片免费观看的视频| 精品国产一区二区久久| 成人国产一区最新在线观看| 成人手机av| 欧美精品亚洲一区二区| 久久天躁狠狠躁夜夜2o2o| 99久久99久久久精品蜜桃| 十八禁高潮呻吟视频| 别揉我奶头~嗯~啊~动态视频| 亚洲专区国产一区二区| 中国美女看黄片| 国产精品1区2区在线观看. | 无人区码免费观看不卡| 黄色 视频免费看| 丝袜人妻中文字幕| 国产亚洲一区二区精品| av福利片在线| 丰满饥渴人妻一区二区三| 一边摸一边抽搐一进一小说 | 人人妻人人澡人人爽人人夜夜| 欧美成人免费av一区二区三区 | 十八禁网站免费在线| 9色porny在线观看| 人人妻人人爽人人添夜夜欢视频| 操出白浆在线播放|