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

    SAR Tomography with Improved Non-Local Means Filtering Based on Adaptive Window

    2024-01-12 14:07:12ShengleiWangZhiyangChenYuanhaoLiChengHu

    Shenglei Wang, Zhiyang Chen, Yuanhao Li, Cheng Hu

    Abstract: In order to mitigate speckle noise in synthetic aperture radar (SAR) images and enhance the accuracy of SAR tomography, non-local means (NL-means) filtering has been proven to be an effective method for improving the quality of SAR interferograms.Apart from considerations like noise type and the definition of similarity, the size and shape of filtering windows are critical factors influencing the efficacy of NL-means filtering, yet there has been limited research on this aspect.This paper introduces an enhanced NL-means filtering method based on adaptive windows,allowing for the automatic adjustment of filtering window size according to the amplitude information of the SAR interferogram.Simultaneously, a directional window is incorporated to align SAR interferograms, achieving the dual objective of preserving filtering standards and retaining detailed information.Experimental results on interferogram filtering and tomography, based on TerraSAR-X data, demonstrate that the proposed method effectively reduces phase noise while maintaining texture accuracy, thereby improving tomography quality.

    Keywords: NL-means filter; adaptive window; SAR interferogram filtering; SAR tomography

    1 Introduction

    Synthetic aperture radar tomography (TomoSAR)extends the principle of synthetic aperture to elevation direction, utilizing multiple SAR images from different observation angle to achieve threedimensional high-resolution imaging [1].It can be applied to fields such as urban three-dimensional imaging and forest parameter inversion [2, 3].The quality of 3D imaging results is significantly contingent upon the quality of the SAR images employed in the imaging process.

    SAR images are often contaminated by speckle noise due to the imaging principles, which can reduce the quality of SAR images and further affect the results of tomography.Therefore,before conducting tomography, the SAR images need to be denoised first.Common filtering methods include mean filtering, median filtering,wavelet transform filtering, Goldstein filtering and so on [4, 5].

    The common speckle processing methods mainly include four categories: spatial speckle suppression algorithm, transform domain speckle suppression algorithm, anisotropic diffusion algorithm and non-local means algorithm.The spatial speckle suppression algorithm is the earlier algorithm used in speckle suppression, including the classical Frost filter [6], Kuan filter [7], etc.Transform domain speckle suppression algorithms are also widely used in speckle processing,including wavelet transform based speckle suppression methods [8], Shearlet transform based speckle suppression methods [9] and so on.anisotropic diffusion (AD) filter [10] is an anisotropic diffusion filter that considers both edge-preserving and speckle suppression.After improvement, it has also been applied to SAR image speckle suppression and has a good speckle suppression effect [11, 12].Compared with the above methods, the non-local means method is a new method for SAR speckle suppression, which has better edge and detail preservation effect,and is suitable for SAR images with complex structure.

    Non-local means filtering is an emerging image denoising method proposed by Buades et al.in 2005 [13].The main idea of this algorithm is to utilize the high redundancy of natural images and the decorrelation of noise.The principle of NL-means filtering is to find patches obtain the filtered patch by weighted averaging of these similar patches.The NL-means algorithm excels in preserving image edges, textures,and other features.Over recent years, it has emerged as a prominent research focus in the realm of image denoising.Through numerous refinements, this algorithm has demonstrated exceptional efficacy in enhancing image quality by reducing noise.

    Due to its effective denoising capabilities,the NL-means filter is increasingly applied in SAR image denoising.In the context of SAR image denoising, the enhancements to the traditional NL-means filter primarily focus on the following directions:

    1) Improved the filtering algorithm based on the noise characteristics of SAR images.For example, Zhong et al.[14] adopted a pixel classification strategy to effectively reduce the influence of the multiplicative speckle model and improve the effectiveness of similar patch search,and designed a method to search rotation invariant similar patches by using orientation information.

    2) Based on the statistical characteristics of SAR images, redefine the similarity in the algorithm according to different probability density functions.For example, based on the assumption that the speckle of SAR multi-view images follows the generalized gamma distribution,Deledalle et al.[15] used the probability similarity of image patches instead of Euclidean distance, and proposed the PPB speckle reduction algorithm based on the probability similarity of image patches.Zhong et al.[16] derived a SAR image speckle reduction algorithm based on Bayesian non-local means based on Bayesian framework, and proposed a statistical distance to measure the similarity of image patches, which has certain speckle reduction effect.

    3) Based on the structural characteristics of SAR images, redefine the similarity in the algorithm according to characteristic parameters.For example, Yi et al.[17] considered the structural similarity of images and introduced the structural similarity index (SSIM) into the non-local means to suppress the speckle of SAR images and realize the optimization of edge-preserving characteristics.Chen et al.[18] proposed a SAR image speckle reduction algorithm based on nonlocal means and coefficient of variation, which also has good filtering and edge preservation characteristics.

    In addition to the above several improvement ideas, the size and shape of the filtering window are also important factors affecting the denoising effect, but there has not been a systematic discussion and analysis.In terms of windows,the current problems of NL-means filtering mainly include:

    1) The window size was fixed, and could not be adjusted according to the image characteristics.

    2) The shape of the window is a fixed square window, which cannot accurately match similar patches according to the texture information of the image.

    To solve these problems, this paper proposes an improved NL-means filtering method,which can adaptively adjust the size of the filtering window according to the amplitude information of the SAR interferogram.At the same time,the direction window is introduced to match the phase fringes, so as to achieve the purpose of preserving the filtering standard while retaining the detail information.

    This paper is organized as follows.Section 2 mainly introduces the basic principle of NLmeans filtering.Section 3 introduces the specific process of the improved NL-means filtering algorithm proposed in this article, as well as the confirmation method of key parameters.Section 4 shows the filtering effect of the proposed NLmeans filtering method for real SAR interferograms, which verifies the good effect of the proposed method in detail preservation of SAR interferograms.Section 5 mainly shows the application of the improved algorithm proposed in this paper to SAR tomography, which verifies that using the SAR interferogram sequence filtered by the proposed method for tomography can indeed improve the quality of 3D imaging.

    2 The Principle of NL-means Filtering

    The main idea of non-local means filter is to utilize the intrinsic redundancy of natural images and the decorrelation of noise.In a typical natural image, various elements such as edges, points,and other features exhibit similarity in multiple regions within the same image [19].The process of NL-means filtering is to traverse all pixels of the whole image once for each pixel to be filtered.During this traversal, the algorithm assesses and quantifies the similarity, assigns weights based on this similarity, and subsequently computes a weighted average to derive the filtered pixel.The basic expression of the algorithm is as follows:

    wherev?(x) denotes the denoised image,v(y)denotes the noisy image,w(x,y) denotes the weight, andΩxdenotes the neighborhood of the pixelx.The determination of weight is the key of the algorithm.In the traditional NL-means algorithm, the weight is related to the Euclidean distance between two image patches, which can be expressed as [13]

    In practical applications, the computation involved in traversing the entire image for similarity determination is notably extensive, and there are numerous patches exhibiting low similarity throughout the image, the computational efficiency is adversely affected.Hence, it is necessary to set the search window and similarity window to reduce the range of pixel traversal.

    Fig.1 is the windows diagram of NL-means filtering.The blue block represents the patch to be filtered, the green block represents the patch used to compute similarity, and the large white window is the search window.During the filtering process, the green block systematically slides,pixel by pixel, within the search window.It calculates similarity with the blue block at each position and executes a weighted average to accomplish the filtering.In practical terms, the window size significantly influences the filtering outcome.According to the non-local means filtering principle, the size of the search window dictates the distance between the similar patch and the patch to be filtered.A larger search window encompasses patches from more distant regions,resulting in a more pronounced filtering effect.Generally, the recommended sizes for the search window and similarity window are denoted by|W|=21×21 and |Δ|=7×7 respectively.This window size not only exhibits robust noise-filtering capabilities but also excels in handling intricate details and fine structures [13].

    Fig.1 The windows diagram of NL-means filtering

    3 Improved NL-means Filtering Method

    While the original NL-means filtering method demonstrates effectiveness in optical image processing, several enhancements are imperative for SAR image processing.The key areas for improvement include: 1) The original NL-means filter employs a fixed window size, devoid of adaptability to the unique characteristics of the image.Consequently, the filtering effect may not be consistently pronounced, and in certain instances, essential details may be inadvertently filtered, especially for images with varying precision; 2) During the computation of weights based on similarity, the substantial difference between the pixel slated for filtering and those within the search window can result in weight instability,thereby compromising the overall filtering effect;3) The fixed shape of the filtering window induces transitional filtering for structures in the image with a distinct directional stripe pattern.This can lead to the loss of detailed information due to blurring.Furthermore, the noise characteristics inherent in SAR interferograms differ from the additive noise found in optical images.Consequently, directly applying the original NLmeans filtering method to SAR interferograms may not yield optimal filtering results.

    Addressing the aforementioned limitations,this paper proposes an improved NL-means filtering method based on the distinctive characteristics of the interferogram itself.The objective is to implement more profound filtering in regions with low signal-to-noise ratios, preserve original detail information in high signal-to-noise ratio areas, and mitigate the interaction between pixels exhibiting significant differences.The definition of the weight in the original NL-means filter is adjusted to enhance the weight stability.Additionally, the size and shape of the search window are adaptively adjusted, ensuring the retention of intricate structural details.This adaptive approach contributes to the overall improvement in the quality of the filtered SAR interferogram.The specific process of the improved NL-means filtering method is as follows:

    Step 1 Input the SAR interferogram, and calculate the SNR for the entire image using the amplitude values derived from the complex image, referencing a designated noise region.

    Step 2 Identify pixels for filtering based on the calculated SNR, the pixels whichSNR ≥15 dBdo not need to be filtered.

    Step 3 Compute the correlation coefficient for the identified pixels according to the SNR,and determine the size of the search window according to the correlation coefficient, that is,the numberLof pixels participating in the final weighted average process.

    Step 4 Select the appropriate direction window to align with the fringe direction and texture characteristics of the interferogram.

    Step 5 Calculate the similarity between each pixel within the direction window and the target pixel slated for filtering.the target pixel slated for filtering, and filter out the topLmost similar pixels for weighted average to complete the filtering of a single pixel.

    Step 6 Traverse all identified pixels across the entire image, executing the complete image filtering process.

    3.1 Definition of Similarity

    In actual SAR interferograms, the distribution of strong and weak scatterers is frequently non-uniform.During NL-means filtering, the pixel intensities within the search window can markedly differ from the intensity of the pixel slated for filtering.This discrepancy is particularly pronounced when filtering weak scatterers, where introducing excessive information from strong scatterers is undesirable and may lead to interference.Consequently, this paper introduces a modification to the weight calculation formula in the original NL-means filter.The aim is to address the issue of elevated weight assigned to strong scatterers due to the lower intensity of the central pixel slated for filtering.To mitigate the issue of excessively high weights assigned to strong scatterers, an enhancement is introduced to the original NL-means filter weight definition.This improvement incorporates a normalization process using the intensity of the central pixel‖v(x)‖2.The refined weight calculation formula, designed to alleviate the influence of the center pixel, is expressed as follows:

    3.2 Introduction of Directional Windows

    In traditional NL-means filtering, the prevalent window shape is predominantly rectangular.However, due to the pronounced directional characteristics of interference fringes in interferometric phase, the conventional rectangular filtering window may disrupt the continuity of fringes,particularly in densely packed fringe regions.The interferometric phase filtering algorithm based on directional window should adopt a directional window that fits the direction of the interference fringe.The pixel values within a directional window exhibit proximity, suggesting approximate homogeneity.Consequently, employing a directional window for filtering serves the dual purpose of noise reduction without compromising the integrity of interference fringes.

    In the classical Lee filter proposed by Lee et al.[20], 16 kinds of directional window templates are used for filtering, and the quality of filtered images is successfully improved.The templates of these 16 directional windows are also introduced in the improved NL-means filtering algorithm proposed in this paper, as shown in Fig.2.During the filtering process, calculations exclusively involve the white pixels within these templates.

    Fig.2 Directional window templates in classical Lee filtering

    The selection of the directional window can be based on the complex interferometric phase,and the specific procedure is to calculate the mean value based on 16 directional windows in the window.Subsequently, the window with the highest amplitude is selected for further processing.

    3.3 Adaptive Windows Size

    3.3.1 Similarity Window Size Selection

    For determining the similarity window, a simulation experiment approach is employed.Firstly,construct a simulation plane scene with phase 0,and add noise with different signal-to-noise ratios.Then, filter the noise by the combination of search window and similarity window with different sizes.Next calculate the RMSE value of the filtered phase to compare the filtering effect.

    Fig.3 The trend of interferogram phase RMSE with filtering window for different SNR: (a) 0 dB; (b) 5 dB; (c) 10 dB

    Fig.3 shows the changes of phase RMSE after filtering in the simulation plane with three kinds of noises added.According to the above experimental results, it can be seen that in NLmeans filtering, the influence of similarity window on filtering effect is much less than that of search window, and when the similarity window is increased to 3×3, the filtering effect is no longer significantly improved.Therefore, the method proposed in this paper does not focus on the change of the similarity window, which is fixed as a 3×3 rectangular window, and the subsequent adjustment of the size of the search window can meet the requirements of adaptive filtering.

    3.3.2 Search Window Size Selection

    The search window is determined by the probability density distribution of the interferogram phase.The size of the search window can be deduced by the probability density function of the interferogram phase when the threshold is set.In the SAR interferogram, the probability density function of the interferometric phase can be expressed as [21]

    The simulation results of the variation trend of the phase standard deviation with the correlation coefficient under different looks are shown in Fig.4.

    Fig.4 Trend of phase standard deviation with correlation coefficient

    It can be seen from Eq.(4) that the probability density function of the interferometric phase is related to the looksLand the correlation coefficientγ, where the looksLis a key parameter used in multi-look processing.The multi-look process is the process of averaging the range or azimuth directions of the SAR image, so the looksLcan be considered as the number of pixels involved in the average, which is similar to the search window size in NL-means filtering.In NL-means filtering, the size of the search window is the number of pixels that eventually participate in the weighted average; therefore, the looksLcan be considered equivalent to the size of the search window in NL-means filtering.

    In the algorithm proposed in this paper, the function of adaptively improving the window size need to be realized by the amplitude information of the SAR interferogram.Therefore, the correlation coefficient in the above derivation needs to be associated with the SNR of the SAR image.The correlation coefficientγhas the relationship with SNR as following functional:

    On the basis of the above derivation, the phase standard deviation at SNR=15 dB is used as the threshold to simulate the variation trend of the phase standard deviation with the looks in the case of different SNR, as shown in Fig.5.It can be seen that in the case of different SNR, as long as the looks is increased to a certain value,the phase standard deviation can reach the standard of SNR=15 dB.According to this feature,the looksLcan be equivalent to the size of the search window, that is, the number of pixels participating in the weighted average, and finally the purpose of adaptively determining the search window can be achieved.

    Fig.5 Trend of phase standard deviation with looks at different SNR

    In order to verify the equivalence between the size of the search window and the looksL,the simplest weighted average filtering method,mean filtering, is used to filter the simulation noise plane.The value of the mean filtering window is set to the value of the looksLwhen the SNR reaches the standard of SNR=15 dB, and test whether the filtering effect of noisy images with different SNR is up to the standard.

    The results of the mean filtering are shown in Tab.1.The simulation shows that when SNR =15 dB, RMSE = 0.404 1.Through the experiment, it can be seen that the phase RMSE can basically reach the standard after the filtering window size is selected according to the simulation looks, and the RMSE after filtering is slightly higher than the standard value because the direction of the filtering window is fixed, butit can still be judged that the looksLis equivalent to the size of the filtering window.

    Tab.1 Comparison of mean filtering results for noisy images with different SNR

    3.4 Determination of Filter Coefficients h

    The main parameter in filtering also includes the filtering coefficienth, which needs to be further analyzed on the basis of the improved weight calculation formula.

    As we know from the previous derivation, if the intensity of the center pixel is much larger than that of the surrounding points, it will not be filtered.Therefore, this paper focuses on the case that the center pixel is much smaller than the surrounding pixels, that is, the weight calculation method whenb ?1, and the improved weight formula can be simply written asw2=exp(-b2/h2).

    Whenb ?1, the intensity of the central pixel is much smaller than that of the surrounding pixels.In this case, the product of the corresponding weight of the surrounding pixels and the pixel should be much smaller than that of the central pixel, so as to reduce the influence of the surrounding pixels on the central pixel.The expression is as follows:

    After simplification, the parametershandbactually need to meet the following conditions:

    The effect of changing the parameterhon the weight calculation is shown in Fig.6.Through the simulation results of the influence ofthe change of parametershandbon the weight,it can be seen that when the parameterh≤2,the condition of Eq.(9) can be satisfied.In NLmeans filtering, the parameterhrepresents the filtering degree, therefore, the larger the value ofh, the more obvious the filtering effect.In the case of satisfying the condition, we still want to obtain more obvious filtering effect.Therefore, in this paper, the value of parameterhis determined as 2.

    4 Interferogram Filtering Experiments

    In order to verify the effectiveness of the proposed method for filtering SAR interferograms, in this section, a measured interferogram is filtered and the filtering results are compared.The methods adopted include the original NL-means filtering, the NL-means filtering with directional window only and the improved NL-means filtering proposed in this paper.

    The original data of the experiment is SAR interferogram generated by two pre-processed(registration, deramping, phase compensation)TerraSAR-X images.The imaging area is the Regent Hotel and some surrounding buildings.The slant range resolution is 0.45 m, azimuth resolution is 0.85 m, range bandwidth is 300 MHz,and perpendicular baseline length is 21.24 m.The specific acquisition parameters are given in Tab.2.

    Tab.2 TerraSAR-X acquisition parameters

    4.1 Interferometric Phase Filtering Effect

    The before-after comparison of the interferometric phase obtained by different filtering methods is shown in Fig.7.From the comparison of interferometric phase before and after filtering, it can be seen that the filtering effect of the original NL-means filtering is very obvious, but it will lead to the destruction of structural information.At the location of the building, the NL filter completely destroyed the details, and the results were improved after adding Lee window, but the details could not be completely preserved.On this basis, adding window filter could achieve the purpose of preserving details as much as possible under the condition of meeting the filtering effect requirements (phase RMSE of 15 dB).

    Fig.7 Comparison of phase filtering effect of SAR interferogram in the vicinity of Regent Hotel: (a) before filtering; (b) original NLmeans filter (Search window 9×9); (c) NL-means filter with directional windows; (d) NL-means filter proposed in this paper

    In order to more clearly show the effect of the proposed method on the interferometric phase filtering, the zoom of line feature of the interferometric phase before and after filtering are intercepted for comparison.As we can see in Fig.8, a line target is circled in the black box.The phase standard deviation of the line target is before filtering is 1.78, and which after filtering is 1.75.It can be seen that the phase continuity on the filtered line target is better, and the phase noise is effectively filtered without destroying the structural information of other surrounding targets.

    4.2 Suppression of Interference from Strong Scatterers

    In addition to has a better filtering effect on the interferometric phase, the proposed method can also suppress the strong scatterers point spread phenomenon.In order to verify the suppression effect of the proposed method on the interference of strong scatterers, the amplitude maps under different filtering cases were obtained for comparison.The before-after comparison of the interferogram amplitude is shown in Fig.9.As can be seen from the unfiltered amplitude, there are many strong scatterers points within the scope of the Regent Hotel building.The reason is that the building has many glass Windows on the facade, and the metal structure of the glass window frame will form the corner reflector effect, so the scattering coefficient is very high.The strong scatterer interference suppression performance of the algorithm can be evaluated by comparing the weak scatterer filtering case within the scope of the Regent Hotel building.

    Fig.8 Zoom of line feature before and after filtering: (a) before filtering; (b) filter proposed in paper

    By observing the SAR interferogram amplitude of the Regent Hotel before and after filtering, it can be seen that the original NL-means filtering will lead to a relatively serious diffusion effect.The addition of Lee window has a certain suppression effect on the diffusion effect, but it is also quite serious.However, the diffusion of strong scattering points is effectively suppressed in the results of the improved filtering method,which also indicates that it is effective to select filter points according to the search window.

    Fig.10 shows the comparison of coherence map before and after filtering.The mean of the correlation coefficient before filtering is 0.48, and which after filtering is 0.53.It can be seen that after filtering, the correlation of all imaging areas is improved, especially the building area, which also proves the effectiveness of the filtering method proposed in this paper.

    Fig.10 Coherence map before and after filtering: (a) before filtering; (b) filter proposed in paper

    In order to quantitatively analyze the filtering effect of the algorithm, Tab.3 shows the comparison of the edge-preserving index (EPI)and the equivalent number of looks (ENL) in the flat region of the phase before and after filtering by different filtering methods.

    Tab.3 Comparison of evaluation indexes of different filtering methods

    It can be seen from the results that although the ENL in the flat region after filtering by the proposed method is lower than that of the original NL-means filtering method, the EPI of the phase is much higher than that of the original NL-means filtering method, which also reflects the optimization effect of the proposed method on phase structure preservation.

    5 Comparison of Tomographic Results

    In order to verify the improvement of the proposed filtering method in SAR tomography, in this section, interferometric filtering-based tomography is performed on a set of TerraSAR-X data to compare the imaging results.This set of data contains 19 interferograms processed from 20 SAR images, and the specific imaging parameters are the same as in the above section.The imaging area is the Regent Hotel Beijing and some nearby buildings.When generating the interferogram, the first SAR image is used as the master image, and the remaining images are used as slave images in turn.The acquisition time and baseline settings of the 20 SAR images are shown in Tab.4, and the spatio-temporal baseline is shown in Fig.11.

    Tab.4 SAR image acquisition time and baseline

    Fig.11 Spatio-temporal baseline for tomography

    The proportion of pixels with different SNR in the scene and the corresponding window sizes are given in Tab.5.

    Tab.5 Window sizes corresponding to different SNRs

    The tomographic results before and after interferogram filtering are shown in Fig.12.From the following three sets of tomographic results, it can be seen that the use of filtered interferograms for tomographic imaging can effectively reduce the number of spurious points in the scene, can make the point cloud more concentrated in the building area, and improve the quality of tomographic imaging.However, the original NL-means filtering will cause building structures besides the main building to be filtered out, and the overall number of targets is reduced, which is effectively improved by the improved method proposed in this paper.

    Fig.12 Comparison of tomographic results before and after interferogram filtering: (a) before filtering; (b) aftering filtering (Original NL-means filter); (c) after filtering (filter proposed in this paper)

    In order to prove the effectiveness of the proposed method in SAR tomography, the imaging effect of the Regent Hotel area is analyzed.Tab.6 shows the points in the complete point cloud before and after filtering and the points in the part of the point cloud of the Regent Hotel building.It can be seen that although the number of points in the filtered complete point cloud and the point cloud of the Regent hotel have decreased, the proportion of the points in the hotel area in the complete point cloud has increased, which indicates that the spurious points have been successfully filtered after filtering.Although the original NL-means filter can increase the proportion of the main buildings in the whole scene, it will lead to a large reduction of the number of points in the complete point cloud, which indicates that the original NLmeans filter will filter out most of the detail information.Compared with the original NLmeans filter, the filtering method proposed in this paper retains more structural information of other buildings on the premise of increasing the proportion of the main buildings.This reflects that the filtering method proposed in this paper has improved the quality of SAR tomography.

    Tab.6 Comparison of points in the point cloud before and after filtering

    6 Conclusion

    To enhance the quality of SAR tomography, it is essential to filter the SAR interferogram beforehand.While the original NL-means filter demonstrates commendable filtering efficacy, it grapples with challenges such as a fixed window size and suboptimal adaptability to speckle noise in SAR interferograms.In response to these limitations, this paper introduces an improved NLmeans filter, leveraging the distinctive features of SAR interferograms.Through adjustments to the similarity definition in the original NL-means filter and adaptive resizing of the search window,this approach facilitates precise and refined filtering of SAR interferograms.

    Filtering experiments conducted on measured interferograms and subsequent tomography experiments with the acquired data affirm the effectiveness of the proposed method.It demonstrates not only a robust denoising capability but also adept detail retention.Nevertheless, the method presented in this paper still harbors areas for improvement.For instance, the fixed shape of the similarity window lacks adaptability to varying image features, and the directional window directly borrows the template from the Lee filter without more innovative modifications.In future research, these aspects will be subject to further exploration and enhancement.

    嫩草影视91久久| 一区二区日韩欧美中文字幕| 亚洲九九香蕉| bbb黄色大片| 国产aⅴ精品一区二区三区波| 高清欧美精品videossex| 自线自在国产av| 日本黄色日本黄色录像| 久久九九热精品免费| 欧美国产精品一级二级三级| av不卡在线播放| 老熟女久久久| 露出奶头的视频| av福利片在线| 中文字幕av电影在线播放| 久久久久久免费高清国产稀缺| 国产精品久久久人人做人人爽| 精品免费久久久久久久清纯 | a级毛片黄视频| 午夜影院日韩av| 国产成人精品在线电影| 欧美日韩一级在线毛片| 热re99久久国产66热| 国产精品久久久久久精品古装| 1024香蕉在线观看| 成年人黄色毛片网站| av一本久久久久| 亚洲色图av天堂| 国产一区在线观看成人免费| 天堂中文最新版在线下载| 操美女的视频在线观看| 天天躁日日躁夜夜躁夜夜| 999久久久精品免费观看国产| 久久精品aⅴ一区二区三区四区| 欧美日韩乱码在线| 亚洲avbb在线观看| 少妇 在线观看| 国产高清激情床上av| 男女下面插进去视频免费观看| 777米奇影视久久| 91在线观看av| 亚洲美女黄片视频| 国产精品久久久久久人妻精品电影| 精品久久蜜臀av无| 亚洲自偷自拍图片 自拍| videos熟女内射| 日韩视频一区二区在线观看| 成人黄色视频免费在线看| 黄片播放在线免费| 成人av一区二区三区在线看| 一个人免费在线观看的高清视频| 国产单亲对白刺激| 亚洲欧美精品综合一区二区三区| 国产精品二区激情视频| 中文字幕人妻熟女乱码| 啦啦啦在线免费观看视频4| 免费在线观看影片大全网站| 捣出白浆h1v1| 久久香蕉激情| 咕卡用的链子| 丁香欧美五月| 最近最新中文字幕大全免费视频| 一区在线观看完整版| 在线观看免费午夜福利视频| 国产蜜桃级精品一区二区三区 | 黑人欧美特级aaaaaa片| 日日夜夜操网爽| 国产乱人伦免费视频| 久久国产乱子伦精品免费另类| 日韩三级视频一区二区三区| tocl精华| 9色porny在线观看| 一级a爱视频在线免费观看| 日韩成人在线观看一区二区三区| 成年动漫av网址| 香蕉丝袜av| 午夜精品久久久久久毛片777| 丰满人妻熟妇乱又伦精品不卡| 91在线观看av| 日韩免费av在线播放| 午夜免费观看网址| 欧美日韩乱码在线| 一区在线观看完整版| 正在播放国产对白刺激| 91在线观看av| 看片在线看免费视频| 精品卡一卡二卡四卡免费| 国产精品免费视频内射| 欧美精品高潮呻吟av久久| 捣出白浆h1v1| 狠狠狠狠99中文字幕| 中文字幕人妻熟女乱码| xxx96com| 水蜜桃什么品种好| 男女午夜视频在线观看| 久久久水蜜桃国产精品网| 国产精品 国内视频| 精品一区二区三区四区五区乱码| 中文字幕最新亚洲高清| 99国产精品一区二区三区| 欧美成人午夜精品| 欧美精品高潮呻吟av久久| 亚洲精品美女久久久久99蜜臀| 免费高清在线观看日韩| 麻豆av在线久日| 99热网站在线观看| 精品午夜福利视频在线观看一区| svipshipincom国产片| 午夜免费鲁丝| 热99国产精品久久久久久7| 久久婷婷成人综合色麻豆| 日韩 欧美 亚洲 中文字幕| 国产精品.久久久| 欧美黑人欧美精品刺激| 久久天堂一区二区三区四区| 免费高清在线观看日韩| 男人操女人黄网站| 午夜久久久在线观看| 精品国产乱子伦一区二区三区| 99精品久久久久人妻精品| 女警被强在线播放| 丰满人妻熟妇乱又伦精品不卡| 美女福利国产在线| 黄网站色视频无遮挡免费观看| 搡老岳熟女国产| 人人妻人人澡人人爽人人夜夜| 久久久久国内视频| 日韩一卡2卡3卡4卡2021年| 国产精品免费大片| 啦啦啦在线免费观看视频4| 激情在线观看视频在线高清 | 极品教师在线免费播放| 国产精品成人在线| x7x7x7水蜜桃| 久久人妻av系列| 欧美国产精品va在线观看不卡| 国产无遮挡羞羞视频在线观看| 热re99久久精品国产66热6| 每晚都被弄得嗷嗷叫到高潮| 亚洲精品在线观看二区| 国产精品久久久人人做人人爽| 满18在线观看网站| 中出人妻视频一区二区| 国产欧美日韩综合在线一区二区| 人人妻人人澡人人爽人人夜夜| 大香蕉久久成人网| 欧美日韩福利视频一区二区| 久久青草综合色| 97人妻天天添夜夜摸| av在线播放免费不卡| 91av网站免费观看| 可以免费在线观看a视频的电影网站| 91精品国产国语对白视频| 国产精品1区2区在线观看. | 久久久精品免费免费高清| 欧美日韩亚洲综合一区二区三区_| 黄片大片在线免费观看| 欧美另类亚洲清纯唯美| 在线国产一区二区在线| 午夜两性在线视频| 免费少妇av软件| 成人亚洲精品一区在线观看| 国产精华一区二区三区| √禁漫天堂资源中文www| 香蕉久久夜色| 婷婷精品国产亚洲av在线 | 王馨瑶露胸无遮挡在线观看| 日韩欧美一区二区三区在线观看 | 亚洲精品一卡2卡三卡4卡5卡| 亚洲av成人av| 久久久国产一区二区| 久久久精品区二区三区| 久久久久久久午夜电影 | 手机成人av网站| 国产精品二区激情视频| 一级片'在线观看视频| 不卡一级毛片| 日韩欧美免费精品| 国产成人av激情在线播放| 婷婷精品国产亚洲av在线 | 999久久久国产精品视频| 少妇裸体淫交视频免费看高清 | 一边摸一边抽搐一进一小说 | 18禁黄网站禁片午夜丰满| 自拍欧美九色日韩亚洲蝌蚪91| 午夜成年电影在线免费观看| 国产aⅴ精品一区二区三区波| 日日夜夜操网爽| 丝袜美足系列| 黄网站色视频无遮挡免费观看| 国产又爽黄色视频| 亚洲精品一二三| 女人久久www免费人成看片| a级毛片在线看网站| 精品久久久久久电影网| 青草久久国产| 久久这里只有精品19| www.精华液| 久久人人爽av亚洲精品天堂| 一本综合久久免费| 精品电影一区二区在线| 亚洲,欧美精品.| 国产精品av久久久久免费| 欧美人与性动交α欧美精品济南到| 欧美国产精品一级二级三级| 成年版毛片免费区| 搡老岳熟女国产| 亚洲七黄色美女视频| 无人区码免费观看不卡| 国产激情欧美一区二区| 悠悠久久av| 捣出白浆h1v1| 一区在线观看完整版| 夜夜躁狠狠躁天天躁| 一级黄色大片毛片| 久久精品aⅴ一区二区三区四区| 99国产精品一区二区蜜桃av | 国产精品二区激情视频| 国产成人啪精品午夜网站| 纯流量卡能插随身wifi吗| 亚洲精品中文字幕一二三四区| 精品国产乱子伦一区二区三区| 啦啦啦免费观看视频1| 一区在线观看完整版| 丰满的人妻完整版| 日韩一卡2卡3卡4卡2021年| 一级a爱视频在线免费观看| 91成年电影在线观看| 高清在线国产一区| 亚洲成人免费av在线播放| 免费看十八禁软件| 女性被躁到高潮视频| 精品人妻熟女毛片av久久网站| 大片电影免费在线观看免费| 亚洲成人免费av在线播放| 两性夫妻黄色片| 亚洲一卡2卡3卡4卡5卡精品中文| 欧美日韩瑟瑟在线播放| 色94色欧美一区二区| 国内久久婷婷六月综合欲色啪| 午夜精品国产一区二区电影| 久久精品91无色码中文字幕| 亚洲一区二区三区不卡视频| 91精品三级在线观看| 成熟少妇高潮喷水视频| 午夜免费成人在线视频| 中文欧美无线码| av不卡在线播放| 久久香蕉国产精品| 久久人人爽av亚洲精品天堂| 成年动漫av网址| 亚洲人成伊人成综合网2020| 精品一区二区三区视频在线观看免费 | 好看av亚洲va欧美ⅴa在| 欧美日韩一级在线毛片| 最新美女视频免费是黄的| 变态另类成人亚洲欧美熟女 | svipshipincom国产片| 侵犯人妻中文字幕一二三四区| 午夜福利欧美成人| 人成视频在线观看免费观看| 精品午夜福利视频在线观看一区| 亚洲专区字幕在线| 新久久久久国产一级毛片| 欧美激情高清一区二区三区| 少妇的丰满在线观看| 又黄又爽又免费观看的视频| 两个人免费观看高清视频| 久久亚洲真实| 制服人妻中文乱码| 久久久国产欧美日韩av| 国产欧美日韩综合在线一区二区| 免费在线观看影片大全网站| 香蕉国产在线看| 美女国产高潮福利片在线看| 男女之事视频高清在线观看| 搡老熟女国产l中国老女人| 国产精品二区激情视频| 最新美女视频免费是黄的| 久久国产精品男人的天堂亚洲| 国产成人一区二区三区免费视频网站| 国产成人欧美| 嫁个100分男人电影在线观看| 免费日韩欧美在线观看| 久久国产精品大桥未久av| 人妻丰满熟妇av一区二区三区 | 少妇的丰满在线观看| 男人舔女人的私密视频| 欧美色视频一区免费| 又紧又爽又黄一区二区| 欧美精品人与动牲交sv欧美| 欧美不卡视频在线免费观看 | 十分钟在线观看高清视频www| 少妇粗大呻吟视频| 国产精品偷伦视频观看了| 曰老女人黄片| 高清黄色对白视频在线免费看| 亚洲av片天天在线观看| 人人妻人人澡人人爽人人夜夜| 午夜免费鲁丝| 宅男免费午夜| av天堂久久9| 亚洲片人在线观看| 久久草成人影院| 日韩欧美国产一区二区入口| 精品无人区乱码1区二区| 国产一区二区三区视频了| 这个男人来自地球电影免费观看| 亚洲精品一卡2卡三卡4卡5卡| 男女之事视频高清在线观看| 又黄又粗又硬又大视频| 精品国产一区二区三区久久久樱花| 免费观看精品视频网站| 99在线人妻在线中文字幕 | 国产精品av久久久久免费| 少妇猛男粗大的猛烈进出视频| 亚洲成av片中文字幕在线观看| svipshipincom国产片| 国产一区在线观看成人免费| 人妻久久中文字幕网| 日本五十路高清| 国产精品美女特级片免费视频播放器 | 大片电影免费在线观看免费| 国产不卡一卡二| 久久草成人影院| 国产日韩欧美亚洲二区| 夜夜夜夜夜久久久久| 亚洲全国av大片| 免费观看a级毛片全部| 黄网站色视频无遮挡免费观看| 色精品久久人妻99蜜桃| 少妇 在线观看| 欧美精品人与动牲交sv欧美| 国产精品一区二区精品视频观看| av超薄肉色丝袜交足视频| 国产精品国产高清国产av | 亚洲熟女毛片儿| 一级,二级,三级黄色视频| 曰老女人黄片| 久久精品91无色码中文字幕| 亚洲精品成人av观看孕妇| 免费人成视频x8x8入口观看| 欧美性长视频在线观看| 亚洲七黄色美女视频| 亚洲自偷自拍图片 自拍| 欧美日韩亚洲综合一区二区三区_| 亚洲国产欧美网| 国产亚洲欧美在线一区二区| 免费在线观看日本一区| 国产精品九九99| 国产激情欧美一区二区| 国产免费现黄频在线看| 91在线观看av| 男女之事视频高清在线观看| 久久久久国产精品人妻aⅴ院 | 丰满饥渴人妻一区二区三| a在线观看视频网站| 亚洲欧美日韩另类电影网站| 亚洲午夜精品一区,二区,三区| 亚洲熟妇熟女久久| 日日爽夜夜爽网站| 久久久国产精品麻豆| 真人做人爱边吃奶动态| 久久久国产精品麻豆| 91精品国产国语对白视频| 国产三级黄色录像| 亚洲午夜理论影院| 国产亚洲精品久久久久5区| 成年人午夜在线观看视频| 国产伦人伦偷精品视频| 国产精品国产高清国产av | 欧美日韩黄片免| 如日韩欧美国产精品一区二区三区| 黄网站色视频无遮挡免费观看| 国产成人影院久久av| 亚洲欧美激情在线| 亚洲色图 男人天堂 中文字幕| 精品国产亚洲在线| av福利片在线| 午夜免费成人在线视频| 国产单亲对白刺激| 精品国产超薄肉色丝袜足j| 他把我摸到了高潮在线观看| 亚洲 欧美一区二区三区| 欧美av亚洲av综合av国产av| 色精品久久人妻99蜜桃| 国产在视频线精品| 欧美日本中文国产一区发布| 国产精品久久电影中文字幕 | 日韩 欧美 亚洲 中文字幕| 免费在线观看亚洲国产| 亚洲中文字幕日韩| 国产有黄有色有爽视频| 最新美女视频免费是黄的| 日韩 欧美 亚洲 中文字幕| 亚洲免费av在线视频| 老熟妇仑乱视频hdxx| 两性夫妻黄色片| 日本黄色视频三级网站网址 | 国产精品1区2区在线观看. | 99香蕉大伊视频| 九色亚洲精品在线播放| 久久精品人人爽人人爽视色| 欧美大码av| 搡老岳熟女国产| 国产欧美日韩综合在线一区二区| 一二三四在线观看免费中文在| 欧美日韩亚洲综合一区二区三区_| 亚洲欧美日韩高清在线视频| 精品一品国产午夜福利视频| 精品一区二区三区四区五区乱码| 精品国产亚洲在线| 男女下面插进去视频免费观看| 王馨瑶露胸无遮挡在线观看| 91精品三级在线观看| 多毛熟女@视频| 免费黄频网站在线观看国产| 国产精品美女特级片免费视频播放器 | 村上凉子中文字幕在线| aaaaa片日本免费| 国产高清视频在线播放一区| 亚洲精华国产精华精| 欧美日韩国产mv在线观看视频| xxx96com| ponron亚洲| 在线观看午夜福利视频| 欧美激情久久久久久爽电影 | 黄片播放在线免费| 新久久久久国产一级毛片| 久久久国产成人免费| 国产精品亚洲一级av第二区| 叶爱在线成人免费视频播放| 看免费av毛片| 他把我摸到了高潮在线观看| 老司机午夜福利在线观看视频| av线在线观看网站| 后天国语完整版免费观看| 精品福利观看| а√天堂www在线а√下载 | 国产精品.久久久| 精品一区二区三区av网在线观看| 色婷婷久久久亚洲欧美| 在线永久观看黄色视频| 久久青草综合色| 天天添夜夜摸| 欧美日韩精品网址| 亚洲色图 男人天堂 中文字幕| 久久久久久久国产电影| 亚洲第一欧美日韩一区二区三区| 国产欧美日韩一区二区三| 精品少妇久久久久久888优播| 人人妻人人澡人人爽人人夜夜| 免费在线观看视频国产中文字幕亚洲| 日韩人妻精品一区2区三区| 激情视频va一区二区三区| 亚洲人成电影免费在线| 国产精品久久久av美女十八| www.熟女人妻精品国产| 岛国毛片在线播放| 亚洲va日本ⅴa欧美va伊人久久| 狠狠狠狠99中文字幕| 在线观看66精品国产| 后天国语完整版免费观看| 亚洲伊人色综图| 自拍欧美九色日韩亚洲蝌蚪91| 久久天堂一区二区三区四区| 99国产精品99久久久久| 天堂中文最新版在线下载| 老司机在亚洲福利影院| 制服诱惑二区| 精品久久久久久,| www.熟女人妻精品国产| 亚洲中文av在线| 久久精品亚洲熟妇少妇任你| 国产精华一区二区三区| 女人高潮潮喷娇喘18禁视频| 成人18禁在线播放| 黄色视频不卡| 岛国在线观看网站| 如日韩欧美国产精品一区二区三区| 国产精品综合久久久久久久免费 | 在线天堂中文资源库| 亚洲一区二区三区不卡视频| 麻豆av在线久日| 午夜精品国产一区二区电影| 久久精品国产99精品国产亚洲性色 | 欧美成人免费av一区二区三区 | 一级片'在线观看视频| 日韩欧美三级三区| 王馨瑶露胸无遮挡在线观看| 亚洲aⅴ乱码一区二区在线播放 | 一二三四社区在线视频社区8| 咕卡用的链子| 国产黄色免费在线视频| 日韩欧美国产一区二区入口| 亚洲伊人色综图| 欧洲精品卡2卡3卡4卡5卡区| 久久久水蜜桃国产精品网| 女人爽到高潮嗷嗷叫在线视频| 国产精品 国内视频| 999精品在线视频| 国产区一区二久久| 亚洲精品一二三| 成人精品一区二区免费| 久久亚洲真实| 亚洲 国产 在线| 精品午夜福利视频在线观看一区| 欧美黄色片欧美黄色片| 精品久久久精品久久久| 嫁个100分男人电影在线观看| 国产精品影院久久| 国产精品永久免费网站| 91成年电影在线观看| 97人妻天天添夜夜摸| 久99久视频精品免费| 精品一区二区三卡| 99精品在免费线老司机午夜| 亚洲国产精品一区二区三区在线| 亚洲欧美日韩另类电影网站| 精品一区二区三区四区五区乱码| 国产精品久久久人人做人人爽| 久久国产亚洲av麻豆专区| 国产男女超爽视频在线观看| 好看av亚洲va欧美ⅴa在| 天堂中文最新版在线下载| 人妻丰满熟妇av一区二区三区 | 欧美日韩精品网址| 大香蕉久久成人网| 婷婷成人精品国产| 久9热在线精品视频| 丰满饥渴人妻一区二区三| 欧美人与性动交α欧美精品济南到| 亚洲专区中文字幕在线| 日韩一卡2卡3卡4卡2021年| 国产野战对白在线观看| 精品午夜福利视频在线观看一区| 国产精品国产av在线观看| 久久久久久久久免费视频了| 精品乱码久久久久久99久播| 制服人妻中文乱码| 免费在线观看完整版高清| 亚洲色图综合在线观看| 久久久久精品人妻al黑| 国产激情欧美一区二区| 91成年电影在线观看| 午夜福利影视在线免费观看| 精品国产乱码久久久久久男人| 国产亚洲精品一区二区www | 王馨瑶露胸无遮挡在线观看| 妹子高潮喷水视频| 久久久久久人人人人人| 免费观看a级毛片全部| 美国免费a级毛片| 人妻 亚洲 视频| 国产精品欧美亚洲77777| 国产不卡av网站在线观看| 国产aⅴ精品一区二区三区波| 一本一本久久a久久精品综合妖精| 欧美最黄视频在线播放免费 | 欧美日韩亚洲高清精品| 欧美日韩福利视频一区二区| 婷婷精品国产亚洲av在线 | 在线观看一区二区三区激情| 成熟少妇高潮喷水视频| 国产av一区二区精品久久| 国产精品一区二区免费欧美| 大香蕉久久网| 成人永久免费在线观看视频| 免费人成视频x8x8入口观看| 久久香蕉精品热| 麻豆国产av国片精品| 亚洲伊人色综图| 成人三级做爰电影| 少妇被粗大的猛进出69影院| 国产99久久九九免费精品| 人成视频在线观看免费观看| 99久久人妻综合| 久久人人爽av亚洲精品天堂| 99国产精品一区二区蜜桃av | 国产精品免费大片| 亚洲人成77777在线视频| 久久久精品区二区三区| av网站在线播放免费| 五月开心婷婷网| 国产亚洲av高清不卡| 亚洲免费av在线视频| 老汉色av国产亚洲站长工具| 国产乱人伦免费视频| 性少妇av在线| www.999成人在线观看| 777久久人妻少妇嫩草av网站| 色婷婷av一区二区三区视频| 99国产精品免费福利视频| 欧美不卡视频在线免费观看 | 国产成+人综合+亚洲专区| 久久久国产成人免费| 国产精品久久久人人做人人爽| av电影中文网址| 欧美人与性动交α欧美精品济南到| 久久久精品免费免费高清| 午夜福利影视在线免费观看| 99香蕉大伊视频| 一进一出好大好爽视频| 两性午夜刺激爽爽歪歪视频在线观看 | 黄色片一级片一级黄色片| a级片在线免费高清观看视频| 热99国产精品久久久久久7| 操出白浆在线播放| 久久精品国产清高在天天线| av不卡在线播放| 国产精品美女特级片免费视频播放器 | 在线观看66精品国产| 国产成人欧美在线观看 | 又黄又爽又免费观看的视频| 很黄的视频免费|