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

    EVALUATION OF THE SHADING METHOD FOR REDUCING IMAGE BLOOMING IN THE PIV MEASUREMENT OF OPEN CHANNEL FLOWS*

    2012-06-27 05:54:10ZHONGQiangLIDanxunCHENQigangWANGXingkui

    ZHONG Qiang, LI Dan-xun, CHEN Qi-gang, WANG Xing-kui

    State Key Laboratory of Hydroscience and Engineering, Tsinghua University, Beijing 100084, China,

    E-mail:zhongq08@mails.tsinghua.edu.cn

    EVALUATION OF THE SHADING METHOD FOR REDUCING IMAGE BLOOMING IN THE PIV MEASUREMENT OF OPEN CHANNEL FLOWS*

    ZHONG Qiang, LI Dan-xun, CHEN Qi-gang, WANG Xing-kui

    State Key Laboratory of Hydroscience and Engineering, Tsinghua University, Beijing 100084, China,

    E-mail:zhongq08@mails.tsinghua.edu.cn

    The shading method is a simple but effective way of reducing image blooming in the measurement of open channel flows with the Particle Image Velocimetry (PIV). The current paper proposes a simplified analytical model for light attenuation using this method. The model is verified against experimental data, and the influence of several parameters is illustrated numerically. The possible adverse effect due to the light attenuation is shown to be limited when the parameters in the shading method are in an adequate range, as shown by processing standard images of Case B in PIV Challenge 03. A simple criterion for setting the shade in experiment is given for controlling the errors caused by the shading technique within an acceptable range.

    Particle Image Velocimetry (PIV), blooming, shade, open channel flows

    Introduction

    The Particle Image Velocimetry (PIV) is a measurement technique widely used in fluid mechanics study. Because it is a 2-D (or 3-D) and undisturbed optical method, it is especially used in water flow researches in recent years, for example, the gas-liquid or solid-liquid two-phase flow[1,2], and the open channel flow with submerged rigid vegetation[3]. Because the PIV is a technique based on the information from images of the flow, the image blooming, which frequently occurs in PIV experiments, becomes an issue.

    The image blooming refers to a situation in which the neighboring pixels are saturated with excess charges, which would lead to streaking, flare, or even complete “whiteout” of the recorded images, ruining the cross-correlation analysis in the PIV measurement[4]. Basically, the output of a pixel in a Charge-Coupled Device (CCD) is related with the amount of electrons and is therefore proportional to the amount of light falling onto it. Each pixel cell can store a maximum number of electrons, as its full well capacity[5], which depends on the CCD type. When this number is exceeded during exposure, the additional electrons migrate to the neighboring pixels, leading to image blooming. In many PIV measurement cases, the objects in the test section always have highly reflective surfaces, so the blooming would frequently occur in places close to these surfaces. Examples are aluminum or steel surfaces in wind tunnels and fluid machineries, glass surfaces in the channel, and solid particles or bubbles in multiphase flows, among others. Information of the blooming area is missing, so the vectors calculated from these areas must be interpreted carefully because with the blooming, their values can significantlydistorted[5], and sometimes the blooming may cause measurement failures near the surface[6]. In addition, the blooming may also damage the sensors of cameras[7]. So the light reflection and the blooming control are the key factors of the PIV measurement[8].

    Several methods are used commonly in experiments to prevent blooming. First is that the PIV recording techniques might be improved. In modern CCD sensors, the blooming can be reduced significantly through incorporating specialized anti-blooming architectures[4], an example is the prevention of the electrons of saturated pixels from polluting adjacent cells by the anti-blooming gate between pixels[9]. However, the use of such devices may limit the pixel size and possibly the resulting pixel sensitivity. Another option is the use of Complementary Metal Oxide Semiconductor (CMOS) cameras[10,11]. One of the main advantages of most CMOS sensors is their capability to record images with high contrast without blooming. However, the image quality and signal-to-noise ratios of CMOS cameras are generally lower than those of CCD cameras[10]. The second method is the limited use of reflective surfaces. Arnott et al.[12]used a matt black, self-adhesive plastic foil to cover the reflective surface for eliminating the area for blooming. In the work of Yu et al.[8], Di et al.[13], Wu et al.[14]and Mula et al.[15], the surfaces in the model were painted with coating to reduce laser blooming effects. The last and a commonly used method is the adjustment of various parameters such as intensity of laser, aperture value, shutter time, ISO, and gain, among others. This method reduces the blooming probability by finding an appropriate combination of parameters. It is more like an art, and the image quality is generally affected because fewer photons hit the pixel and the sensitivity of CCD is also affected.

    In open channel flow experiments, all aforementioned “anti-blooming” methods have serious limitations. The first method requires an increased investment in hardware. The second one fails completely as the channel bed, usually made of glass, must remain clear and transparent in the test section to allow passage of laser light for the bottom-to-top illuminetion. The third method is difficult to apply in practice, as the optimization of illumination parameters requires tremendous experience and patience. Other approaches are not commonly used. For example, highly reflective objects could be placed outside the measurement windows to avoid camera blooming[16]. However, this method is unsuitable because the PIV test section is usually determined by the problems to be studied.

    Another commonly used method, termed as the“shading method” in the current study, proves to be a simple but effective way of reducing the image blooming. The idea of the shading method is not really new. In fact, it is also a widely used approach, mostly as one of the last resorts when all other options have failed. For example, Voges et al.[17]used this method in studying the rotor blade tip interaction with a casing treatment in a transonic compression stage by PIV. The reflected light from the moving blades led to image blooming, and they shaded the blade with a piece of sheet metal immediately in front of the observation window, resulting in the drastic reduction of blooming. Therefore, the shading method is more convenient than many other methods.

    However, the influence of the shading method in the PIV experiment and its evaluation are lack in literature. Therefore, the current study aims to examine the effectiveness of the shading method and to discuss its possible adverse effects on image quality and PIV calculation, and to provide some guidance in setting the shade in experiment.

    Fig.1 Sketch of the shading method

    1. Light attenuation by the shading method

    1.1 Sketch of the shading method

    In an open channel flow measurement with PIV, the image blooming often occurs in the vicinity of the channel bed because of the intensive light reflections caused by the glass bed or by the clumping of particles on the bed, as illustrated in Fig.1, where the plane AC is the target and is illuminated by a laser sheet. In this article, the streamwise direction is defined as x, and the wall-normal direction is defined as y, and the coordinate origin point is set on the bed and at the middle of the channel. The lens is reduced to a single convex lens; its diameter is designated by Φ, and the distance from the lower edge to the plane of the channel bed is designated byH. By attaching a small shade to the flume wall, the high-intensity light coming from the near-bed region is blocked from entering a part of the CCD lens, i.e., the DE part of the lens is prevented from “seeing” the channel bed (point A). Before the shade is set, the reflective light from point A can go into the camera through the dotted line AD and the solid line AF. All parts of the lens can receive light from A. When the shade is set, the light can only go through the solid line AE and AF, so the shaded part of the lens cannot receive light from A anymore. Therefore, the gray value of point A in the images will be reduced, so will the risk of blooming. In addition, the effect of shade is limited. Figure 1 clearly shows that B is the demarcation point. The height of DE is reduced to zero as A goes up to point B, beyond which the influence of the shade vanishes. Therefore, the shading method is a local influence method, it only reduces the amount of light near the bed but not other parts falling into the camera, and the near-bed band is the only area where the blooming takes place frequently. The risk of blooming will be reduced near the bed, and at the same time, the image quality will be kept as good as before.

    The distance of point B above the channel bed iswhere h is the height of the shade, d1is the distance of the shade to the target plane AC, d2is the distance between the shade and the lens, andν= 1/1.3333 is the refractive index. Note that for simplicity, the influence ofν is not shown in Fig.1.

    When y<LAB, the height of the “blocked” part of the lens is

    When y>LAB, the entire lens is free from the blocking effect, i.e., m=0 and whenν2d2+1(h-y)2(ν2-1)<0, m=Φ. Formula (2) shows that mincreases with h and decreases with H. A combination of h and H should be determined carefully to ensure thatm <Φ, otherwise, the entire CCD lens will be blocked.

    1.2 Analysis of the light-attenuation factor

    Let OV(y) represent the original image gray value at point y, which changes to CV(y) when the shading is in effect, i.e.,

    where α(y) is the light intensity attenuation factor of the shading method. α(y) is assumed to be proportional to the area of the CCD lens that receives light, i.e.,

    where θ=arccos(1-2m/Φ). Expanding the second term of the numerator in the right part of formula (4),

    As 1/2Φmsinθ>0,θ>sinθ when 0≤θ≤π, therefore, 1/4Φ2sinθ-Φ/2msinθ>0. According to formula (4),α≤1. Formula (4) can be simplified to the following form

    Based on Formula (4), the effect of the light attenuation can be estimated approximately. Therefore, the value of α at every point is

    Fig.2 Influence of hon the light attenuation factor

    1.3 Numerical illustration of the light attenuation factor

    The light attenuation factor, as quantified in Formula(5), involves some complexity. The following shows the numerical variations of α(y) with d1/d2and h. Under the condition of Φ=0.02m, H= 0.003m, d1+d2=0.4m , d1/d2=1 and with the water level is 0.03 m, the characteristics ofα are shown inFig .2(a). α=1 ata very small h, which indicates that the shade has no effect in blocking thelight when it is small. Ash increases,α decreases, whereas L, the size of the influence area, increases. The variation in a gray value gradient factor (1-α)/LABwith h is shown in Fig.2(b). The gray value gradient at h=0.01m is about five times as great asthatat h=0.0015m.

    Fig.3 Influence of d1/d2on the light attenuation factor

    The variation in α with d1is shown in Fig.3 for the caseof Φ=0.02m, H=0.003m,h= 0.006m, d1+d2=0.4m, and with the water level is 0.03 m. The value of αincreases ac cordingly as d1increases as shown inFig.3(a). As expected, the influence area also increases. Figure 3(b) shows a dramatic decrease in (1-α)/LABwith d1.

    In practical PIV experiments, parameters h and d1/d2can be determined ba s ed on Formula (5) with the above illustrations in mind.

    Fig.4 Schematic diagram of the experimental setup

    1.4 Experimental illustration of the light attenuation factor

    Formula (5) is checked against experimental data fromPIV measurements conducted in a 10.8 m long, 0.25m wide hydraulic flume with a glass bed and glass side walls. A 2W Nd-YAG laser provides the illumination. Two identical CCD cameras (AVT Pike F032B ASG16 with lens of EF 0.05 m f/1.2L USM) were placed symmetrically on both sides of the flume (Fig.4). Their heights and other parameters were kept the same for capturing the same series of images from the same illuminated sheet. The difference is that the shading method was used only on the left side.

    Images were captured under different sizes of the shade, e.g., h=0, 0.003 m, 0.005 m, 0.008 m, and 0.01m, and the gray level of each pixel in an image was obtained.Based on the linear relationship between the photon hits and the pixel saturation, the effect of light attenuation is quantified by comparing the gray values to the benchmark values at h=0. To ensure the statistical significance, a total of 8 000 images were averaged for each case (frame rate is 50 Hz and sampling time is 160 s).

    Table 1 Parameters of the experiment

    The parameters of the CCD, as well as those of the flow, are listed in Table 1. Note that the gain and thegamma value of the CCD were set to zero to ensure the linear relationship between the photon hits and the pixel saturation. The shutter time was 10 000 μs. To prevent the seeding particles from becoming streamwise bands, or the existence of dagger-like stripes in the images at high velocity flows, the experiments were conducted under a very small flow velocity.

    Figure 5(a) shows the average gray value of the imagescaptured without the shading effect (h=0m). The gray value reaches its maximum on the be d an d decre ases dramatica lly as the dista nce abovethechannelbedincreases,indicatingthattheimage blooming occurs most likely in the near-bed region. Figure 5(b) shows the experimental results of the light attenuation factor for each case. The experimental results agree well with Formula (5).

    Fig.5Experimental results of light attenuation

    Table 2 Experiment parameters

    2. Effect of the shading method

    The effect of the shading method is investigated in the experiments. The flow parameters, the aperture value, the shutter time, the ISO value, and the gain in the experiment are all set within the normal range and are listed in Table 2. Figure 6 shows the comparison of the sampled images and the corresponding vector fields calculated with the same PIV algorithm. The image pair consists of two images taken simultaneously using two cameras (see Fig.4). Due to the lasersheet thickness, the shape, the size, and the brightness of the corresponding particles shown in Figs.6(a) and 6(b) are not exactly the same. Figure 6(a) shows that the blooming image is contaminated by several bright vertical streaks and a flaring bright layer in the vicinity of the channel bed, which affects greatly the velocity calculation. In contrast, the image quality shown in Fig.6(b) (with shading, h=0.003m) is enhanced substantially, and the corresponding vector field appears faultless. Figure 6 shows that the shading method can reduce the risk of blooming and improve the quality of PIV data near the bed.

    Fig.6 Comparison of images and calculated velocity vectors. The field of view is 0.0225 m×0.0181 m

    To quantify the effect of the shadingmethod, the risk of blooming is defined as a parameter in the presentstudy. When the image blooming occurs, the gray level of a pixel reaches its maximum value of 255 in an eight-bit BMP image. Therefore, counting the number of points whose gray value is equal to 255 is an accurate way of showing the blooming influence. The risk of blooming at a certain pixel can be defined as the ratio of N/NT, where N is the number of images in which the gray level of the pixel is equal to 255, and NTis the total number of images in a series, and NT=8000 in this article.

    The size of the shade h is an important parameter in the shading method. An increase of h leads to adecrease of the risk of blooming. The effects of the shading method under four different sizes of h are compared. The other parameters of the experiments are listed in Table 2. The image capturing rate is 50 Hz, and a total of 8 000 images are gathered for each experimental run.

    Figure 7(a) shows the percentage R of blooming along the y-axis.For clarity, the axes are set inlogarithmic scales. When there is no shade (solid line, h=0m), about 10% of the points at the bed are suffered with blooming, and this ratio decreases rapidly to ab out 0.04% at y=0.0005m. The risk of blooming drops rapidly when the shading method is used, e.g., from 10% to 0.5%, 0.002%, and nearly zero at h= 0.003 m, 0.005 m, and 0.008 m, respectively. At a higher flow area (y>0.01m), however, the shading method has no effect at all. This finding is further confirmed by the risk ratio at various h values over that without shading h=0m (Fig.7(b)). Figure 7 shows that the shading method can reduce the risk of blooming near the bed and does not affect the upper part as expected.

    Fig.7 Risk of image blooming at various shading sizes

    3. Influence of shading method on PIV calculation and a rule for setting the shade

    3.1 Influence of shading method on PIV calculation

    The shading method is effective inreducing the risk of image blooming. However, there remains one major question to be anwered: Does the light attenuation have an adverse effect on the PIV correlation and calculation? This question is answered by applying thefrom a Direct Numerical Simulation (DNS) simulation of a turbulent open channel flow with the following variables: the flow depth is 0.3 m, the maximum velocity Umax=1.19m/s, the shear velocity uτ= 0.0341m/s, the viscosity ν= 0.01616m2/s, and the Reynolds number Reτ= huτ/ν=640.Theimage diameter is 1.3 pixels for a unit standard deviation and 2.6 pixels for two standard deviations (the standard deviation is 0.65 pixels).

    Table 3 Various parameters (units: m)

    Fig.8Variation in the light attenuation factor imposed on the original images

    The shading effect is simulated by transformingthe o riginal images through mu ltiplying them with the light attenuation factor, as in Formula (5). The parameters are listed in Table 3, and the attenuation factors are shown in Fig.8. The original and transformed images are analyzed by the window displacement iterative multigrid algorithm[19-21], which uses an iterative multi-grid approach for cross-correlation and a three-point Gaussian peak fit for the determination of displacement with a sub-pixel accuracy. The grid size is 16×4. The iteration number is 2, and the IW sizes are 64×16, 32×8 and 16×4.

    Figure 9(a) shows the comparison of the mean velocity profile with the original data for the five cases. The agreement between the different terms is very good. Figure 9(b) shows the difference ε between the shading values and the original value for the mean streamwise velocity. For y>100 pixels, the agreement among the results of all cases is very good. In fact, good results are also obtainednear the bed, with deviations all within ±0.015 pixels except for the first case (No. 1). The result of the first case

    Fig.9 Mean velocity profile obtained from the five cases compared with the original data

    Fig.10 Streamwise and wall-normal turbulence intensities obtained from the five cases compared withthe originnal data

    has a more significant deviation near the bed. Figures 10(a) and 10(b) show the profiles of the results for streamwise and wall-normal turbulence intensities, Urmsand Vrms. The agreement is particularly good except for the points near the bed, where one sees some deviations due to the gray level gradient and the reduction of the particles, and due to the low value of α.

    Fig11 Cumulative histogram of the error on the streamwise and wall-normal velocities near the bed

    Another way of assessing the accuracy is to plothe cumulative histogram of the velocity differences near the bed. From Fig.8, the shortest influence area is about 60 pixels for No. 1, therefore, the velocity differences of 14 points near the bed for all five cases are taken into account. The results are shown in Fig.11, which shows significant differences among the different cases. The results of No. 5, 4, 3 and 2 are particularly good, where more than 80% errors are less than 0.01 pixels for both streamwise and wall-normal velocities. The result of No.1 is relatively poor, and only about 60% errors are less than 0.01 pixels for both streamwise and wall-normal velocities. A few large errors reaching dozen pixels appear in No.1 because the gray value is very low near the bed due to the low value of α (see Fig.8), with a limited amount of information. Figure 12 shows the profiles of the root mean square (rms) errors ε for the five groups. The conclusion is similar to that of the other figures. The rms errors are all within about 0.08 pixelfor streamwise velocities and 0.035 pixel for wallnormal velocities except No. 1.

    Fig.12 RMS errors on the streamwise and wall-normal velocity

    3.2 A ru le for setting the shade

    The above discussion shows that the shading will introduce some small errors in the PIV calculation.of the pixel size (see Fig.11), which is smaller than that in the experimental system. For example, the errors of many PIV processing methods are in the order of 1 of the pixel size[16]. Therefore, this error is acceptable in actual experiments. When the parameters are not adequate, errors will increase, for example, as in the case of No.1.

    The error in the shading method is due to the gray level gradient in the interrogation window and the reduction in the number of particles, and more accurately, is due to the valueofαat the bed, as shown in Fig.8. The value of α at the bed is a most important parameter that affects the characteristics of errors. From foregoing discussions and Fig.8, we may give a simple criterian for determining α at the bed, to control the error characteristics in the PIV calculation

    Based on Formula (5), this criterion can be written as

    In this criterion, all parameters can be obtained easily except forΦ. Φis the diameter of the single convex lens, which isdetermined by the real lens used inthe experiment. It may be estimated roughly through the diameter of the lens aperture, or can be obtained by fitting the experimental data with Formula (5). After all the parameters in the Formula (7) are known, one may determine whether the experiment conditions meet this criterion. If this criterion is satisfied, the errors caused by the shading technique will be controlled within an acceptable range; if not, it would be best to adjust the test parameters to meet this simple criterion.

    4. Conclusion

    The present article has evaluated the shading method for reduc ing the image blooming in the PIV measurement. The effect of the method on the light attenuation isquantified by a simplified analytical model and is tested by experimental data. The results show that the risk of blooming near the channel bed is reduced substantially when the shading is implemented. To facilitate its practical use, the influence of parameters on the image gray level is discussed. The standard images of Case B in PIV Challenge 03 are used for evaluating the adverse effects of the conversion of the light attenuation on the PIV calculation. The result shows that when the parameters in the shading method are adequate, the shading method causes negligible errors on the PIV calculations. Based on the discussion, a simple criterion for setting the shade in the experiment is given, and when the criterion is satisfied, the errors coming from the shading technique will be controlled within an acceptable range. Generally speaking, the shading method is a convenient and efficient method for the reduction of image blooming in the PIV measurement of open channel flows. The proposed model of the light attenuation can be used when applying the shading method in practice.

    [1]LI Shi-rong, CHENG Wen and WANG Meng et al. The flow patterns of bubble plume in an MBBR[J]. Journal of Hydrodynamics, 2011, 23(4): 510-515.

    [2]ZHANG Jin-feng, ZHANG Qing-he and LU zhao. Lattice Boltzmann simulation of particle settling and experimental validation by PIV[J]. Advances in Water Science, 2009, 20(4): 480-484(in Chinese).

    [3]WU Fu-sheng, JIANG Shu-hai and YANG Xue-lin. Characateristics of 2D-vortex field in open channel flow with submerged rigid vegetation[J]. Chinese Journal of Hydrodynamics, 2010, 38(2): 256-259(in Chinese).

    [4]SENTER J., SOLLIEC C. Flow field analysis of a turbulent slot air jet impinging on a moving flat surface[J]. International Journal of Heat and Fluid Flow, 2007, 28(4): 708-719.

    [5]RAFFER M., WILLERT C. E. and WERELEY S. T.et al. Particle image velocimetry a practical guide[M]. Germany: Springer-Verlag, 2007, 69-71.

    [6]MA Chong, LI Hua-xing and MENG Xuan-shi et al. PIV study of induced flowfield on steady and periodicpulsed dielectric-barrier-discharges[J]. Chinese Journal of Applied Mechanics, 2011, 28(3): 308-312(in Chinese).

    [7]DI FELICE F., PEREIRA F. Developments and applications of PIV in naval hydrodynamics particle image velocimetry: New developments and recent applications[M]. Germany, Berlin: Springer-Verlag, 2008, 475-503.

    [8]YU Xiao, LUO Xiang and XU Guo-qing et al. Particle Image Velocimetry(PIV) measurements of the flow in a rotating cavity with a radial inflow[J]. Journal of Aero- space Power, 2009, 24(11): 2483-2488(in Chinese).

    [9]FISHER C. W., MIDDLETON N. T. and LINNE M. A. et al. Phase sensitive imaging system development for laser dlagnostic systems with inherently large background constraints[C]. Lasers and Electro-Optics, Conference on Lasers and Electrooptics. Washington DC, USA,1999, 15.

    [10]HAIN R., KAHLER C. J. and TROPEA C. Comparison of CCD, CMOS and intensified cameras[J]. Experiments in Fluids, 2007, 42(3): 403-411.

    [11]EKOTO I. W., BOWERSOX R. D. W. and SRINIVASAN R. et al. Examination of cameras for near-wall PIV measurements in high speed flows[C]. 44th AIAA Aerospace Sciences Meeting and Exhibit. Reno, USA, 2006, 1-10.

    [12]ARNOTT A., SCHNEIDER G. and NEITZKE K. P. et al. Multi-window PIV for high-lift measurements[C]. 20th International Congress on Instrumentation in Aerospace Simulation Facilities. G?ttingen, Germany, 2003, 44-53.

    [13]DI Yi-bo, DING Jue and WENG Pei-fen. Experiments on viscous flows around low-aspect-ratio wing at low Reynolds numbers[J]. Journal of Shanghai University (Natural Science Edition), 2009, 15(4): 399-403(in Chinese).

    [14]WU P., S?LLSTR?M E. and UKEILEY L. et al. Structural dynamics, Volume 3[M]. New York: Springer, 2011, 1441-1451.

    [15]MULA S. M., STEPHENSON J. and TINNEY C. et al. Vortex jitter in hover[C]. American Helicopter Society Southwest Region Technical Specialists Meeting. Fort Worth, Texas, USA, 2011, 182-206.

    [16]DI FELICE F., DI FLORIO D. and FELLI M et al. Experimental investigation of the propeller wake at different loading conditions by particle image velocimetry[J]. Journal of Ship Research, 2004, 48(2): 168- 190.

    [17]VOGES M., SCHNELL R. and WILLERT C. et al. Investigation of blade tip interaction with casing treatment in a transonic compressor–Part I: Particle image velocimetry[J]. Journal of Turbomachinery–Transa- ctionsof the ASME, 2011(1), 133: 1-11.

    [18]STANISLAS M., OKAMOTO K. and KAHLER C. J. et al. Main results of the second international PIV challe- nge[J]. Experiments in Fluids, 2005, 39(2): 170-191.

    [19]SCARANO F., RIETHMULLER M. L. Iterative multigrid approach in PIV image processingwith discrete window offset[J]. Experiments in Fluids, 1999, 26(6): 513-523.

    [20]SCARANO F., RIETHMULLER M. L. Advances in iterative multigrid PIV image processing[J]. Experiments in Fluids, 2000, 29(Suppl. 1): S51-S60.

    [21]SCARANO F. Iterative image deformation methods in PIV[J]. Measurement Science and Technology, 2002, 13(1): 1-19.

    September 24, 2011, Revised January 5, 2012)

    * Project supported by the National Natural Science Foundation of China (Grant No. 50779023).

    Biography: ZHONG Qiang (1985-), Male. Ph. D. Candidate

    LI Dan-xun,

    E-mail: lidx@mail.tsinghua.edu.cn

    亚洲国产欧美人成| 久久精品国产亚洲网站| 欧美一级a爱片免费观看看| 人妻制服诱惑在线中文字幕| 欧美日韩在线观看h| 天堂√8在线中文| 午夜福利在线观看吧| 欧美3d第一页| 国产精品电影一区二区三区| 不卡一级毛片| 爱豆传媒免费全集在线观看| 久久国内精品自在自线图片| 久久精品国产亚洲av天美| 人人妻人人看人人澡| 国产一级毛片七仙女欲春2| 亚洲欧美成人综合另类久久久 | 久久久久久久午夜电影| 国产精品一区www在线观看| 中文在线观看免费www的网站| 日本熟妇午夜| 亚洲美女视频黄频| 18禁在线无遮挡免费观看视频| 国产单亲对白刺激| 日本熟妇午夜| .国产精品久久| 日本五十路高清| 久久精品影院6| 97热精品久久久久久| 99久国产av精品| 国产亚洲精品av在线| 亚洲激情五月婷婷啪啪| 日韩强制内射视频| 夜夜夜夜夜久久久久| 日韩一区二区视频免费看| 免费观看在线日韩| 国产成人精品婷婷| 国产精品一区www在线观看| 青春草视频在线免费观看| 九九在线视频观看精品| 国产精品av视频在线免费观看| 观看美女的网站| 日韩一区二区三区影片| kizo精华| 国产一区二区三区av在线 | 成人高潮视频无遮挡免费网站| 久久99精品国语久久久| 日韩成人av中文字幕在线观看| 久久精品国产亚洲网站| 不卡视频在线观看欧美| 精品国产三级普通话版| 能在线免费看毛片的网站| 精品久久久久久成人av| 变态另类成人亚洲欧美熟女| 99热只有精品国产| 国产片特级美女逼逼视频| 真实男女啪啪啪动态图| 精品人妻偷拍中文字幕| 最近视频中文字幕2019在线8| 热99re8久久精品国产| 国产激情偷乱视频一区二区| 日本五十路高清| 免费看美女性在线毛片视频| 国产91av在线免费观看| 夫妻性生交免费视频一级片| 亚州av有码| 春色校园在线视频观看| 日韩一区二区视频免费看| 日韩一本色道免费dvd| 亚洲内射少妇av| 久久久国产成人免费| 亚洲人成网站高清观看| 久久精品91蜜桃| 丰满人妻一区二区三区视频av| 99riav亚洲国产免费| 亚洲av成人精品一区久久| 日本成人三级电影网站| 国产精品国产三级国产av玫瑰| а√天堂www在线а√下载| 欧美极品一区二区三区四区| 免费av不卡在线播放| 久久99热这里只有精品18| 成人无遮挡网站| 中文字幕人妻熟人妻熟丝袜美| 看十八女毛片水多多多| 亚洲av免费在线观看| 男人的好看免费观看在线视频| 久久人妻av系列| 久久精品91蜜桃| 少妇的逼好多水| 可以在线观看的亚洲视频| 精品久久久久久久人妻蜜臀av| 国内精品美女久久久久久| 在线观看一区二区三区| 色5月婷婷丁香| 国产又黄又爽又无遮挡在线| 午夜视频国产福利| 精品国内亚洲2022精品成人| 黄色配什么色好看| 国产精品久久久久久久电影| 三级经典国产精品| 国产女主播在线喷水免费视频网站 | 欧美精品一区二区大全| 成人国产麻豆网| 亚洲国产精品sss在线观看| 久久精品夜夜夜夜夜久久蜜豆| 深夜精品福利| 日韩三级伦理在线观看| 亚洲内射少妇av| 男女视频在线观看网站免费| 又爽又黄无遮挡网站| 欧美变态另类bdsm刘玥| 一边亲一边摸免费视频| 综合色av麻豆| 午夜福利在线观看吧| 亚洲国产色片| 久久草成人影院| 如何舔出高潮| 老师上课跳d突然被开到最大视频| 亚洲第一电影网av| 亚洲不卡免费看| 亚洲成a人片在线一区二区| 日韩一区二区三区影片| 成年免费大片在线观看| 国产人妻一区二区三区在| 国产真实乱freesex| 最近的中文字幕免费完整| 嘟嘟电影网在线观看| 成人永久免费在线观看视频| 99久久精品热视频| 色播亚洲综合网| 日本三级黄在线观看| 中文亚洲av片在线观看爽| 亚洲av成人av| 蜜桃亚洲精品一区二区三区| 夫妻性生交免费视频一级片| 淫秽高清视频在线观看| 亚洲美女搞黄在线观看| 成人毛片60女人毛片免费| 精品久久久久久久久av| 蜜桃久久精品国产亚洲av| 看免费成人av毛片| 精品午夜福利在线看| 熟女人妻精品中文字幕| 欧美人与善性xxx| 少妇人妻一区二区三区视频| 亚洲高清免费不卡视频| 亚洲三级黄色毛片| 99国产极品粉嫩在线观看| 青青草视频在线视频观看| 午夜精品国产一区二区电影 | 久久中文看片网| 欧美激情在线99| 又黄又爽又刺激的免费视频.| 中文亚洲av片在线观看爽| 亚洲美女视频黄频| 1024手机看黄色片| 黄色配什么色好看| 久久久精品94久久精品| 国产精品久久久久久久久免| 欧美日韩一区二区视频在线观看视频在线 | 一级毛片久久久久久久久女| 国产高清三级在线| 日本与韩国留学比较| 一个人观看的视频www高清免费观看| 国产成人精品一,二区 | 国产色爽女视频免费观看| 日韩一本色道免费dvd| 岛国在线免费视频观看| 国产欧美日韩精品一区二区| 亚洲欧洲国产日韩| 精品一区二区三区视频在线| 两性午夜刺激爽爽歪歪视频在线观看| 国产男人的电影天堂91| 欧美成人免费av一区二区三区| 欧美日韩综合久久久久久| 精品久久久久久久久av| 又爽又黄a免费视频| 成人毛片a级毛片在线播放| av在线天堂中文字幕| 亚洲不卡免费看| 日韩人妻高清精品专区| 欧美xxxx黑人xx丫x性爽| 国产精品电影一区二区三区| 一区福利在线观看| 少妇人妻精品综合一区二区 | 青春草亚洲视频在线观看| 欧美高清性xxxxhd video| 日本黄大片高清| 99九九线精品视频在线观看视频| 黄色欧美视频在线观看| 亚洲精品成人久久久久久| 久久久色成人| 99久久九九国产精品国产免费| 五月玫瑰六月丁香| 亚洲图色成人| 亚洲精品乱码久久久久久按摩| 欧美激情在线99| 免费人成在线观看视频色| 91午夜精品亚洲一区二区三区| 69人妻影院| 国产在线精品亚洲第一网站| 少妇被粗大猛烈的视频| 亚洲欧美中文字幕日韩二区| 人人妻人人澡人人爽人人夜夜 | 男女啪啪激烈高潮av片| 国产成人影院久久av| 国产伦理片在线播放av一区 | 国产精品蜜桃在线观看 | 欧美成人a在线观看| 成人午夜精彩视频在线观看| 99久久中文字幕三级久久日本| 少妇人妻精品综合一区二区 | 国产高清有码在线观看视频| 看片在线看免费视频| 国产不卡一卡二| 高清毛片免费观看视频网站| 熟妇人妻久久中文字幕3abv| 亚洲内射少妇av| 精品久久久久久成人av| 日韩精品青青久久久久久| 天天躁日日操中文字幕| 高清在线视频一区二区三区 | 免费av毛片视频| 国产精品久久视频播放| 亚洲人成网站高清观看| 麻豆一二三区av精品| 色播亚洲综合网| 91精品国产九色| 亚洲在线自拍视频| 国产黄片美女视频| 精品熟女少妇av免费看| 国产高清激情床上av| 久久午夜福利片| 99久久成人亚洲精品观看| 色综合色国产| 久久国产乱子免费精品| av免费在线看不卡| 日韩强制内射视频| 午夜免费激情av| 国产在视频线在精品| 色5月婷婷丁香| 免费看a级黄色片| 丰满乱子伦码专区| 性色avwww在线观看| 国产91av在线免费观看| 久久精品国产亚洲av涩爱 | 免费人成视频x8x8入口观看| 亚洲av不卡在线观看| 国产成人freesex在线| 精品一区二区三区人妻视频| 免费电影在线观看免费观看| 六月丁香七月| 男人狂女人下面高潮的视频| 国产亚洲精品av在线| 91午夜精品亚洲一区二区三区| 国产黄a三级三级三级人| 男人的好看免费观看在线视频| 蜜桃久久精品国产亚洲av| 99久久九九国产精品国产免费| 最好的美女福利视频网| 长腿黑丝高跟| 国模一区二区三区四区视频| 91精品国产九色| 久久久久久久久久久免费av| 深夜a级毛片| 免费看美女性在线毛片视频| 中出人妻视频一区二区| 国产av在哪里看| 国产精华一区二区三区| 小说图片视频综合网站| 女人十人毛片免费观看3o分钟| 亚洲精华国产精华液的使用体验 | kizo精华| 欧美一区二区国产精品久久精品| 精品久久久噜噜| 国内揄拍国产精品人妻在线| 久久久久性生活片| 国产亚洲精品久久久久久毛片| 中国国产av一级| 精品久久久久久成人av| 久久人人精品亚洲av| 可以在线观看的亚洲视频| 伦理电影大哥的女人| 国产精品,欧美在线| 最后的刺客免费高清国语| 国产精品久久电影中文字幕| 中文亚洲av片在线观看爽| 欧美一区二区亚洲| 亚洲不卡免费看| 波野结衣二区三区在线| 日韩制服骚丝袜av| 欧美一区二区亚洲| 成人二区视频| 婷婷亚洲欧美| 菩萨蛮人人尽说江南好唐韦庄 | 精品午夜福利在线看| 毛片女人毛片| 日韩欧美在线乱码| 悠悠久久av| 精品久久国产蜜桃| 少妇人妻精品综合一区二区 | 日韩欧美一区二区三区在线观看| 如何舔出高潮| 91aial.com中文字幕在线观看| 国产片特级美女逼逼视频| 亚洲人成网站在线播放欧美日韩| 性插视频无遮挡在线免费观看| 欧美日韩国产亚洲二区| 麻豆国产97在线/欧美| 午夜久久久久精精品| 国产一区二区在线观看日韩| 中文字幕制服av| 免费看a级黄色片| 51国产日韩欧美| 不卡一级毛片| 1024手机看黄色片| 久久精品91蜜桃| 能在线免费看毛片的网站| 老熟妇乱子伦视频在线观看| 亚洲欧洲日产国产| 精品无人区乱码1区二区| 久久久欧美国产精品| 精品人妻视频免费看| 超碰av人人做人人爽久久| 干丝袜人妻中文字幕| 国产免费一级a男人的天堂| 最后的刺客免费高清国语| 亚洲av二区三区四区| 搡老妇女老女人老熟妇| 亚洲在久久综合| 欧美性感艳星| 亚洲无线观看免费| 我要看日韩黄色一级片| 91精品一卡2卡3卡4卡| 亚洲无线在线观看| 看黄色毛片网站| 久久人人爽人人爽人人片va| 亚洲欧美日韩高清在线视频| 一边亲一边摸免费视频| 精品久久久久久久人妻蜜臀av| 一边亲一边摸免费视频| 成人综合一区亚洲| 一边亲一边摸免费视频| 啦啦啦韩国在线观看视频| 国产亚洲91精品色在线| 别揉我奶头 嗯啊视频| 边亲边吃奶的免费视频| 久久精品人妻少妇| 久久亚洲国产成人精品v| av在线观看视频网站免费| 嫩草影院新地址| 国产高清三级在线| 最好的美女福利视频网| 午夜福利在线观看吧| 免费av观看视频| 日韩欧美精品v在线| 18禁在线无遮挡免费观看视频| 午夜久久久久精精品| 日韩人妻高清精品专区| 26uuu在线亚洲综合色| 熟女人妻精品中文字幕| 99久久精品一区二区三区| 天堂中文最新版在线下载 | 久久人人爽人人爽人人片va| 性欧美人与动物交配| 人人妻人人澡人人爽人人夜夜 | 国产午夜精品论理片| 亚洲精品粉嫩美女一区| 啦啦啦啦在线视频资源| 男人舔奶头视频| 日韩 亚洲 欧美在线| 成人毛片a级毛片在线播放| 午夜爱爱视频在线播放| 亚洲真实伦在线观看| 色视频www国产| 99国产极品粉嫩在线观看| 日本撒尿小便嘘嘘汇集6| 日本与韩国留学比较| 成人二区视频| 国产高清激情床上av| 久久精品国产鲁丝片午夜精品| 日韩欧美一区二区三区在线观看| 国产精品麻豆人妻色哟哟久久 | 日韩欧美一区二区三区在线观看| 又黄又爽又刺激的免费视频.| 深夜精品福利| 午夜视频国产福利| 我要搜黄色片| 国产白丝娇喘喷水9色精品| 国产精品久久视频播放| 天堂中文最新版在线下载 | 亚洲人成网站在线播放欧美日韩| 久久久色成人| 亚洲精品日韩在线中文字幕 | 男女做爰动态图高潮gif福利片| 日本在线视频免费播放| 丰满的人妻完整版| 亚洲人成网站高清观看| 久久久久久久午夜电影| 搞女人的毛片| 久久午夜福利片| 国产精品久久久久久久电影| 精品一区二区三区视频在线| 日本五十路高清| 免费看日本二区| 亚洲精品成人久久久久久| 内地一区二区视频在线| 久久久久久久久久黄片| 国产午夜精品一二区理论片| 嫩草影院精品99| 熟妇人妻久久中文字幕3abv| 亚洲高清免费不卡视频| 国产精品一区二区性色av| 日韩 亚洲 欧美在线| 久久99精品国语久久久| 亚洲人成网站在线播| 久久久久久伊人网av| 成人美女网站在线观看视频| 久久精品久久久久久噜噜老黄 | 内地一区二区视频在线| 日本与韩国留学比较| 午夜激情福利司机影院| 国内少妇人妻偷人精品xxx网站| 性插视频无遮挡在线免费观看| 成人一区二区视频在线观看| 99久久精品热视频| 亚洲一区高清亚洲精品| 亚洲中文字幕一区二区三区有码在线看| 欧美zozozo另类| 在线国产一区二区在线| 亚洲美女视频黄频| 久久99热6这里只有精品| 两个人的视频大全免费| av视频在线观看入口| 在线国产一区二区在线| 国产精品免费一区二区三区在线| 人妻久久中文字幕网| 精品熟女少妇av免费看| 免费搜索国产男女视频| 成人无遮挡网站| 一个人看的www免费观看视频| 午夜免费激情av| 国产高潮美女av| 日韩一区二区三区影片| 久久精品夜夜夜夜夜久久蜜豆| 观看免费一级毛片| 色哟哟·www| 国产伦在线观看视频一区| 狠狠狠狠99中文字幕| 日本成人三级电影网站| 亚洲欧美日韩高清专用| 日日啪夜夜撸| 中文字幕制服av| 国产精品免费一区二区三区在线| 看黄色毛片网站| 亚洲av第一区精品v没综合| 国产精品.久久久| 99久久精品国产国产毛片| 91久久精品电影网| 成年免费大片在线观看| 日韩欧美一区二区三区在线观看| 成人欧美大片| 观看美女的网站| 在线a可以看的网站| 国产精品久久久久久精品电影小说 | 亚洲欧美精品专区久久| 99久久人妻综合| 国产成年人精品一区二区| 国产成人午夜福利电影在线观看| 国产高清有码在线观看视频| 免费观看在线日韩| 国产视频内射| 免费av毛片视频| 国产精品野战在线观看| 国产一区二区亚洲精品在线观看| 国产日韩欧美在线精品| 99在线人妻在线中文字幕| 国产亚洲av片在线观看秒播厂 | 国产精品爽爽va在线观看网站| 99热全是精品| 久久久久久九九精品二区国产| 国产精品久久视频播放| 亚洲av二区三区四区| 亚洲欧美日韩卡通动漫| 伊人久久精品亚洲午夜| 在线观看美女被高潮喷水网站| 国产精品不卡视频一区二区| 男女边吃奶边做爰视频| 亚洲丝袜综合中文字幕| 久久人人爽人人爽人人片va| 黄色欧美视频在线观看| 亚洲熟妇中文字幕五十中出| 久久人妻av系列| 一本—道久久a久久精品蜜桃钙片 精品乱码久久久久久99久播 | 午夜激情福利司机影院| 久久精品人妻少妇| 国产中年淑女户外野战色| 此物有八面人人有两片| 久久精品夜色国产| 国产乱人视频| 成人性生交大片免费视频hd| 看免费成人av毛片| 身体一侧抽搐| 一本久久中文字幕| 欧美区成人在线视频| 亚洲婷婷狠狠爱综合网| 两性午夜刺激爽爽歪歪视频在线观看| 免费观看人在逋| 麻豆成人av视频| 国产精品女同一区二区软件| 91精品一卡2卡3卡4卡| 乱人视频在线观看| av在线亚洲专区| 观看美女的网站| 亚洲最大成人中文| 日本av手机在线免费观看| 欧美日韩国产亚洲二区| 午夜亚洲福利在线播放| 天天一区二区日本电影三级| av视频在线观看入口| 亚洲色图av天堂| 国产精品久久久久久av不卡| 白带黄色成豆腐渣| 久久韩国三级中文字幕| 小说图片视频综合网站| 亚洲精品国产av成人精品| 亚洲欧美精品自产自拍| 国产精华一区二区三区| 1000部很黄的大片| 亚洲人与动物交配视频| 蜜臀久久99精品久久宅男| 日韩一区二区三区影片| 97人妻精品一区二区三区麻豆| 又爽又黄无遮挡网站| 日韩欧美一区二区三区在线观看| 日韩在线高清观看一区二区三区| 性色avwww在线观看| 自拍偷自拍亚洲精品老妇| 久久人妻av系列| 少妇裸体淫交视频免费看高清| 亚洲第一区二区三区不卡| 国产亚洲欧美98| 久久久久久久久久黄片| 免费看光身美女| 麻豆成人午夜福利视频| 3wmmmm亚洲av在线观看| 国产伦精品一区二区三区四那| 夜夜爽天天搞| 2022亚洲国产成人精品| 久久久a久久爽久久v久久| 青青草视频在线视频观看| 亚洲中文字幕一区二区三区有码在线看| 国产精品伦人一区二区| 亚洲内射少妇av| 一个人看的www免费观看视频| 非洲黑人性xxxx精品又粗又长| 久久久久久伊人网av| 欧美日韩国产亚洲二区| 在线观看午夜福利视频| 色哟哟哟哟哟哟| 精品一区二区三区视频在线| 国产精品人妻久久久久久| 亚洲乱码一区二区免费版| 国产精品伦人一区二区| 少妇被粗大猛烈的视频| 成年av动漫网址| 亚洲av.av天堂| 国产精品电影一区二区三区| 国产男人的电影天堂91| 日韩欧美在线乱码| 国产黄色小视频在线观看| 黄色配什么色好看| 亚洲七黄色美女视频| 夫妻性生交免费视频一级片| 久久九九热精品免费| 我的老师免费观看完整版| 人妻少妇偷人精品九色| 精品人妻偷拍中文字幕| 久久精品国产亚洲网站| 精华霜和精华液先用哪个| 精品久久久久久成人av| 3wmmmm亚洲av在线观看| av黄色大香蕉| 国国产精品蜜臀av免费| 看片在线看免费视频| 亚洲欧洲日产国产| 18+在线观看网站| 69人妻影院| 插逼视频在线观看| 成人性生交大片免费视频hd| 日韩精品青青久久久久久| 成人鲁丝片一二三区免费| 亚洲欧洲国产日韩| 又黄又爽又刺激的免费视频.| 一级毛片电影观看 | 久久午夜亚洲精品久久| 大型黄色视频在线免费观看| 免费看日本二区| 国产私拍福利视频在线观看| 99久久精品一区二区三区| 精品免费久久久久久久清纯| 一个人看的www免费观看视频| 久久午夜亚洲精品久久| 亚洲精品久久久久久婷婷小说 | 国产精品美女特级片免费视频播放器| 亚洲真实伦在线观看| av黄色大香蕉| а√天堂www在线а√下载| 国内少妇人妻偷人精品xxx网站| 午夜久久久久精精品| 成人午夜高清在线视频| 亚洲精品亚洲一区二区| 真实男女啪啪啪动态图| 久久韩国三级中文字幕| 亚洲成人精品中文字幕电影|