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

    Optimal Gridding Process for GMI Brightness Temperature Using the Backus?Gilbert Method

    2021-09-13 02:41:10GuangcanCHENandYunfeiFU
    Advances in Atmospheric Sciences 2021年11期

    Guangcan CHEN and Yunfei FU*,2,3

    1School of Earth and Space Sciences, University of Science and Technology of China, Hefei 230026, China

    2Hubei Key Laboratory for Heavy Rain Monitoring and Warning Research, Institute of Heavy Rain,China Meteorological Administration, Wuhan 430205, China

    3Anhui Meteorological Observatory, Hefei 230031, China

    ABSTRACT Satellite microwave instruments have different field of views (FOVs) in different channels.A direct average technique(“direct method”) is frequently used to generate gridded datasets in the earth science community.A large FOV will measure radiance from outside the area of a designated grid cell.Thus,the direct method will lead to errors in a measurement over a grid cell because some pixels covering areas outside of the cell are involved in the averaging process.The Backus?Gilbert method (BG method) is proposed and demonstrated to minimize those uncertainties.Three sampling resolutions (6.5 km ×6.0 km,11.5 km×6.0 km,13.0 km×6.0 km) are analyzed based on the scanning characteristics of the Global Precipitation Measurement (GPM) Microwave Imager (GMI) 18.9-GHz channel.Brightness temperatures (TBs) at 0.5 km×0.5 km resolution over eastern China are used to obtain synthetic 18.9-GHz TBs at the three sampling resolutions.The direct and BG methods are both applied to create a 25 km×25 km gridded dataset and their related uncertainties are analyzed.Results indicate the error variances with the direct method are 3.00,3.68 and 4.99 K2 at the three sampling resolutions,respectively.By contrast,the BG method leads to a much smaller error variance than the direct method,especially over areas with a large TB gradient.Two GMI orbital measurements are applied to verify the BG method for gridding process is reliable.The BG method could be utilized for general purpose of creating a gridded dataset.

    Key words:Backus?Gilbert method,optimization,satellite passive microwave imager

    1.Introduction

    Passive microwave remote sensing is a commonly used method to obtain meteorological parameters and land surface parameters.Microwaves are electromagnetic waves with long wavelength,which can penetrate cloud,and can be detected day and night.Signals obtained by microwave instruments are utilized in the retrieval of cloud liquid/ice water,precipitation,soil moisture,vegetation index,etc.(Liu et al.,2001;Liu and Fu,2007;Wang et al.,2009,2014,2019a,b),Over the last few decades,many microwave instruments have been launched—for instance,the Soil Moisture Active?Passive imager and the Micro-Wave Humidity Sounder,as well as other microwave imagers that use the conical scanning technique,such as the Special Sensor Microwave Imager/Sounder,the Global Precipitation Measurement (GPM) Microwave Imager (GMI),etc.,to monitor changes in the Earth’s atmosphere,oceans,and land.Although scan modes of passive microwave instruments are different,the retrieval and algorithm principles based on multi-band information are similar.

    Due to platform limitation,most instruments only have a paraboloidal antenna,and the antenna gain has an inverse association with wavelength when the size of the antenna is constant,meaning the fields of view (FOVs) of the high-frequency channels are smaller than those of the low-frequency channels (higher resolution).The different representativeness and beam-filling effect will cause error in the final products (Greenwald et al.,1997;Rapp et al.,2009).Up to now,the different spatial resolutions in different frequency channels,the beam-filling effect,and the overlapping effect of contiguous pixels are still the three main problems faced by passive microwave remote sensing.The overlapping effect can be used to obtain a high resolution for passive microwaves,especially at the lower frequencies (Fu et al.,2013),while the other problems need further study.

    The latest generation GMI is more advanced than the TMI (Tropical Rainfall Measuring Mission (TRMM)Microwave Imager).However,the low-frequency channels of GMI still have poor resolution.The FOV of GMI at 10.6-GHz is about 32.1 km×19.4 km,which is larger than the pixel sampling interval (6 km×12 km).Usually,the smaller the sampling interval of adjacent pixels,the larger the overlapping part.Furthermore,GMI uses the conical scanning technique that causes different sampling intervals between pixels in different parts of the scan line (see Fig.1).These different sampling intervals generate different overlapping effects.

    On the other hand,gridding meteorological data is very common and useful in the atmospheric sciences.For example,3A GPROF (GPM Profiling Algorithm) produces global 0.25°×0.25° (~25 km) gridded means using Level 2 GPROF data (Kummerow et al.,2001).Njoku et al.(2003)averaged the brightness temperatures onto 0.25°×0.25°grid cells before retrieving the soil moisture.Deeter and Vivekanandan (2006) also averaged the 37-and 89-GHz brightness temperatures observed by ASMR-E onto 0.25° ×0.25° grid cells,in order to retrieve the liquid water path.There are many benefits to using gridded data.For instance,the size of the dataset is reduced,which will benefit data storage.Furthermore,doing so can facilitate comparisons of retrieval product from different instrument platforms.

    Generally,researchers use the direct averaging technique (hereafter,“direct method”) to produce gridded datasets.The direct method averages all pixels in the grid box with equal weight,which is recognized as being both a simple and effective approach.Deeter and Vivekanandan(2006),for example,indicated that the direct method can reduce instrumental noise.At the same time,Heng and Fu(2014) studied TMI’s liquid water path errors at different grid resolutions.They found that the 0.25° gridded dataset retained local details in the instantaneous pixel data and that these data were suitable for analyzing weather activities from the meso to synoptic scale.However,there is a dense overlapping effect between pixels of microwave instruments.The direct method will count the overlapping part twice or more,but the non-overlapping parts only once.Therefore,the weights of different areas will vary.At the same time,the direct method does not consider the FOVs of pixels.The FOV of GMI at 10.6-GHz (32.1 km×19.4 km)is larger than the size of 0.25°×0.25° grid cells.When such pixels fall at the edge of the grid cells,half of the pixels’FOVs sample radiance lies outside the grid cells.With this in mind,in the present paper,a gridded dataset at higher resolution is obtained on the premise of ensuring these errors are as small as possible.

    Backus and Gilbert (1968) first proposed the Backus?Gilbert (BG) method in 1968.After many years of subsequent development,this method,which can reconstruct the brightness temperature of target pixels as observed by a virtual microwave imager,has since emerged as a complete theoretical system.It was first applied to data processing in the field of earth exploration.Green (1975) used the method to retrieve gravity profiles.Then,Stogryn (1978) attempted to use the BG method to match the FOVs of different channels of microwave imagers,in which a cost function was defined that balanced noise amplification and the matching effect.Since then,the BG method has been widely used in matching the FOVs of different channels of microwave imagers (Poe,1990;Farrar and Smith,1992;Long and Daum,1998;Rapp et al.,2009;Wang et al.,2011;Petty and Bennartz,2017;Chen et al.,2018;Zhou and Yang,2020).In these studies,the antenna patterns of pixels are considered.Obviously,the closer the antenna pattern and grid cell size,the smaller the error caused by the gridding method.This paper uses this as a basis to discuss the advantages of the BG method.

    This paper introduces the advantages of the BG method to the gridding process for low-frequency channels of microwave imagers and in so doing,promotes the BG method.The rest of the paper is organized as follows:section 2 introduces the GMI,antenna pattern,and BG method.Results from a simulation experiment are used in section 3 to quantitatively analyze the errors caused by the direct and BG methods.Then,some examples of applying the BG method to actual orbit data are presented in section 4.Finally,a discussion and conclusions are presented in section 5.

    2.Instrumentation and methods

    2.1.GMI

    GMI is a conical-scanning microwave imager that forms the main part of the GPM core platform.It has 13 channels,all of which—except 23-GHz and 183-GHz—are equipped with double polarization.The 13 channels are divided into two groups (S1 and S2,respectively).The lowfrequency group (S1) has an incidence angle of 52°,and the instrument samples 221 points in one scan cycle with a scan range of about 152°.The high-frequency group (S2) parameters are the same as those of the low-frequency group,except for the incident angle,which is about 49°.The FOVs of the 10-,18.9-and 89-GHz channels are 32.1 km (along-track direction)×19.8 km (cross-track direction),18.1 km× 11.7 km and 7.2 km× 6.4 km,respectively (Petty and Bennartz,2017).Since GPROF officially tries to match all channels’FOVs to that of the 18.9-GHz channel (Passive Microwave Algorithm Team Facility,2017),this paper is based on the characteristics of the 18.9-GHz channel.

    In order to intuitively understand the scanning characteristics of GMI,the spatial distribution of pixels observed by the low-frequency channel group (see Fig.1) is displayed.As indicated,distances between adjacent pixels vary.The distance at the edge is larger than that at the center.There are different overlapping parts at different scan positions because of different sampling intervals between adjacent pixels.In this paper,the shortest distance between adjacent pixels with different scan numbers on the scan line is calculated(see Fig.2).Figure 2a is a partly magnified plot of Fig.1.Dashed lines connect adjacent pixels (blue) with the same scan numbers,while solid lines connect pixels whose distances are shortest.Figure 2a also shows that the dashed and solid lines do not coincide in the along-track direction,which means that the distance between adjacent pixels with the same scan number is not the shortest.Curves of the shortest distance along the scan line are shown in Fig.2b.Black and red lines in Fig.2b represent the cross-and along-track direction,respectively.The black line is close to a straight line,indicating that the distances between adjacent pixels in the cross-track direction are similar,while the red line shows a large difference in distance with scan numberr.Five regions can be clearly distinguished by the different slopes of the red line,and the demarcation points are around the scan numbers of 52,93,130 and 172.The shortest distance also changes,from 3.2 km at the edge region to 13.0 km at the center region.In this paper,a scan line is divided into three regions—namely,“edge”,“subedge” and “center” regions.The sampling intervals of these three regions are 6.5 km (along-track direction)×6.0 km(cross-track direction),11.5 km×6.0 km and 13.0 km× 6.0 km,respectively.This paper is based on these three sampling intervals.

    Fig.2.(a) Partly magnified plot of Fig.1b,in which the pixels (blue) on the dashed line have the same scan number and the pixels on the solid line have the shortest distance.(b) Curves of shortest distance between pixels along the scan line,where the red line indicates the along-track direction and the black line the cross-track direction.

    2.2.Antenna pattern and projection

    The Sinek function or Gaussian gain function is generally used to approximate antenna pattern.This paper follows previous research (Wang et al.,2011) and uses a Gaussian gain function to describe the antenna pattern:

    In this application,w represents the beam width.In actual satellite observations,GMI uses the conical-scanning technique,wherein pixels at the surface approximate an ellipse,and the sizes of pixels are roughly the same.The parameters,such as beam width (w),observation location(H),azimuth (φ),and incident angle (θ),together determine the size and shape of pixels.Because the incident angle of GMI is about 52°,all locations need to be projected from the surface coordinate system to the instrument coordinate system to calculate the antenna response.In the projection process,it is assumed that the observation axis of the instrument is the negative direction of the z-axis,and the intersection of the axis and the ground is taken as the origin to establish a new Cartesian coordinate system.The orthogonal basis (A) of the new coordinate system is:

    Then,the coordinates of the point Vpon the ground in the new coordinate system can be expressed as:

    In Eq.(3),z equals zero,meaning the point is on the surface.The square of the distance from Vpto the new z-axis is:

    The projection is completed after replacing the x in Eq.(1) with d to obtain the antenna response.

    Figure 3 compares two antenna patterns with beam widths of 50 km and 11.7 km before and after projection (H=407 km,φ =0°,θ=52°),respectively.In Fig.3a,the curve is compressed in the direction of the satellite and stretched in the other direction.This situation is not evident in Fig.3b.It can be seen that the projected antenna pattern gives better agreement with the actual instrument measurement than the unprojected one,especially for channels with wide beam widths.Furthermore,for the 18.9-GHz channel of GMI,projection processing seems unnecessary.In this paper,however,the 18.9-GHz antenna pattern is projected.

    2.3.BG method

    A measured brightness temperature corresponds to a convolution of the original brightness temperature field using its antenna as a convolution kernel.The brightness temperature (TA) of an existing pixel can be written as

    where x0,y0represents the position of the center of the actual measured pixel.TB(x,y) is the original brightness temperature,and G is the sensor original antenna pattern.

    The brightness temperature at a target pixel can be represented by a linear combination of the adjacent pixels’ brightness temperatures,

    Fig.3.Normalized response curves of the antenna in the along-track direction before (red) and after (black)projection.The azimuth was set at 0°,the incident angle at 52°,the beam width at 50 km,and the satellite height at 407 km for (a),and the other parameters were kept the same but the beam width changed to 11.7 km for (b).

    where xd,ydrepresent the center position of a target pixel (virtual pixel); xi,yirepresent the center position of an actual measured pixel;qirepresents the weight coefficient,which is the solution to Eq.(15);and n is the number of adjacent pixels.In this paper,all pixels in the 50 km×50 km gridded cells are chosen.The gridded cells have the same center as the target gridded cells (25 km×25 km).

    Eq.(5) is substituted for Eq.(6) to give:

    It can be seen that the most important part of Eq.(7) is to solve all coefficients (qi) .Eq.(7) is simplified to:

    The function f can be viewed as an equivalent antenna pattern,and brightness temperature () would be observed by a designed imager with an equivalent antenna pattern.

    Since f is the equivalent antenna pattern,it has the constraints of an antenna pattern:

    where Eq.(12) describes the theoretical brightness temperature of the target pixel.F in Eq.(11) and Eq.(12) is a designed antenna pattern of the virtual instrument.For example,in applications of mapping high-frequency channels (high-resolution) to low-frequency channels (low-resolution),F represents the low-frequency channel antenna pattern and G in Eq.(5) as the high-frequency channel antenna pattern.In this paper,the defined principle of F is to ensure that the weight at any position within the grid cell is equal and zero outside the grid,and G is the 18.9-GHz antenna pattern.J in Eq.(11) is used to deal with the sidelobe,and it has been pointed out that the effect of this part is limited,so this parameter is always set to 1 (Wang et al.,2011).To minimize the cost function,only the equivalent antenna pattern[f in Eq.(9)] can be closest to the designed target one [F in Eq.(12)].

    Stogryn (1978) noted that there is noise in actual satel lite observations.The noise (e2) of BG method process can be expressed as

    where E is the error covariance matrix,andrepresent“native” instrument noise.The sum of squares of weight coefficients (qi) can be viewed as an error amplification factor.When the qiare all positive and equal,the error amplification factor is minimized,and the BG method corresponds to the direct method.If the cost function is considered as an evaluation criterion for improving the resolution,then Eq.(13)can be considered as an evaluation criterion for the effect of error amplification while improving the resolution.Finding the balance between error amplification and resolution improvement involves dealing with a trade-off,Petty and Bennartz (2017) proposed a new cost function (C):

    By adjusting γ,a trade-off between resolution improvement and error amplification is considered.In this paper,an error amplification factor no greater than 1 is required.The solution to Eq.(14) is:

    where δijis the Kronecker delta,G is the sensor original antenna pattern,and F is the designed antenna pattern.

    3.Simulation experiment

    In order to better understand the errors involved in the direct method,a simulation experiment was designed.Eastern China was selected as the study area in the experiment(Fig.4a),and the range was 1000 km×1000 km with a pixel resolution of 0.5 km×0.5 km.The presence of ocean and lakes in the study area was consistent with actual detection.The horizontal distribution of the synthesized brightness temperature was generated by one-to-one correspondence between surface vegetation attributes and brightness temperatures in the original grid cells (Fig.4a),and the vegetation-attribute data were obtained from the MCD12Q1 data product (available at https://lpdaac.usgs.gov/products/mcd 12q1v006/).This approach,as referenced in a previous study (Limaye et al.,2006),while somewhat arbitrary,but we believe is good enough for such a simulation experiment.The raw data were gridded by the direct method to obtain the 25 km×25 km spatial resolution (Fig.4b),and these results are considered as the “true” values.

    Fig.4.(a) Horizontal distribution of brightness temperature (0.5 km×0.5 km) in the simulated region.(b) Horizontal distribution of gridded brightness temperature (25 km×25 km),as processed by the direct method,from (a).The values in (b) are denoted as “true” values for the sake of clarity.

    To investigate influence of the three sampling intervals on the results using the direct method,the distribution of pixels and grid lines in the simulation area was generated.The horizontal distribution of pixels (dots) and grid lines(black lines) is shown in Fig.5 for the sampling interval of 6.5 km×6 km (edge).Ellipses in Fig.5 represent the FOV of 18.9-GHz pixels (18.1 km×11.7 km).Figure 5 shows that a pixel’s measurement is a large range convolution of brightness temperatures,and overlaps exist between adjacent pixels.The figure also shows that,when a pixel is at the grid edge,its measurement will include radiance from outside the grid cell.This implies that the direct method is not appropriate for processing data that are densely overlapped and have a large FOV,where the FOV of the processed data is clearly larger than the size of a grid cell.

    For the three sampling intervals defined in the previous section (i.e.,edge,sub-edge,and center),the horizontal distributions of pixels for 6.5 km (along-track direction) ×6 km (cross-track direction),11.5 km×6 km,and 13 km ×6 km were generated,respectively,using the generation approach in Fig.5.The antenna pattern of pixels is the same as that of the GMI 18.9-GHz channel and is projected;other observation parameters are set to w=11.7 km,H=407 km,φ=0°,and θ=52°,respectively.Each pixel’s brightness temperature is calculated by convolving the brightness temperature field (0.5 km×0.5 km) with the 18.9-GHz antenna pattern.The horizontal distributions of synthesized brightness temperatures for the three pixel sampling intervals are shown in Fig.6.The distribution of brightness temperatures in the edge region is finer than that of the center region.The brightness temperature distribution in all three regions loses detail compared to the raw data,which is due to resolution degradation.

    Fig.5.Horizontal distribution of pixels and grids when the pixel interval is about 6.5 km×6 km in the simulation experiment.Solid lines indicate the boundaries of the grid(25 km×25 km),ellipses represent the FOVs of 18.9-GHz,points (regardless of color) represent pixel locations,and blue points indicate pixels used by the BG method.

    Fig.6.Horizontal distribution of synthesized brightness temperature for three intervals:(a) edge,6.5 km×6.0 km;(b) subedge,11.5 km×6.0 km;(c) center,13.0 km×6.0 km.The method for generating pixel locations for the three sampling intervals is shown in Fig.5.GMI’s 18.9-GHz FOV is applied for each pixel of the three sampling intervals to obtain the synthesized brightness temperature.

    Fig.7.(a?c) Horizontal distribution of gridded brightness temperature (25 km×25 km),as processed by the direct method,for three intervals of pixels:(a) edge;(b) sub-edge;(c) center.(d?f) Horizontal distribution of the differences between the gridded brightness temperature and true values (see Fig.4):(d) edge;(e) sub-edge;(f) center.

    The gridded results at the 25-km resolution for the three synthesized distributions in Fig.6 are shown in Fig.7 using the direct method.To highlight the differences and facilitate comparison,a different colorbar scale is used in Fig.7.The three distributions in Figs.7a?c are different from Fig.4b,which shows the true values,but the entire distribution is the same.This shows that the direct method is a simple method when precision is not required.Results from comparing the gridded results with the true values are shown in Figs.7d?f.The errors in Fig.7d are smaller than those in Fig.7f,which indicates that the results from using the direct method have smaller error for the edge region than for the center region.In addition,the figure shows that there are large errors in regions of strong gradients in brightness temperatures,such as over the coastline,indicating that the direct method will produce larger error where the gradient of the grid cell is significant.This situation is frequently found in actual observations,because the sizes of mesoscale convective systems are smaller than the FOV of the GMI’s 18.9-GHz channel.In regions such as the edge of a convective system,a strong gradient in brightness temperature is generated by the non-uniform distribution of cloud phase and height.

    In essence,the errors of the direct method are mainly caused by the fact that the equivalent antenna pattern does not coincide with the designed antenna pattern.In other words,brightness temperature information from outside the grid cell is used in this method.The BG method is introduced to reduce this error.The advantage of the BG method is that it synthesizes an antenna pattern that is as consistent as possible with the target one by convolution and deconvolution.To compare the equivalent antenna patterns of the direct and BG methods,one grid cell in Fig.5 (50 km?75 km in the x-direction and 25 km?50 km in the y-direction) was selected as an example.16 pixels are used in the direct method.In this paper,all of the 72 blue pixels within the 50 km×50 km grid cell are involved in the BG method.The center of the 50 km×50 km grid cell coincides with the target grid cell center.So all of the 72 blue pixels have radiance contributions from this target grid cell.To strike a balance between resolution improvement and error amplification,γ is set to 3×10?5,so the error amplification factor is less than 1.0.After determining the weights of target grid cells from Eq.(15),an equivalent antenna pattern of the grid cells can be calculated using Eq.(9).

    The BG method’s equivalent antenna pattern is shown in Fig.8b.To discuss the advantages of the BG method,the antenna patterns of the direct and BG methods in the alongand cross-track directions are compared in Figs.8c and 8d.Black curves in Figs.8c and 8d are extracted from Fig.8a along the two solid lines (along-and cross-track),and have been normalized.Red and blue curves are the same as black curves but extracted from Fig.8b and the designed antenna pattern,respectively.The direct method’s antenna patterns in Figs.8c and 8d (black) are significantly different from the designed one (blue) in both the along-and cross-track direction.The distance between the two half-power (50% of max response) points of the BG method is closer to 25 km than that of the direct method.Therefore,the equivalent antenna pattern of the BG method (red) is reasonable.

    Fig.8.Equivalent antenna pattern [f in Eq.(9)] of the example grid in Fig.5 using the (a) direct and (b) BG method.Black curves in (c,d) are extracted from (a) along the two solid lines (along-and cross-track),and it has been normalized.Red and blue curves are the same as black curves but extracted from (b) and the designed antenna pattern,respectively.

    Figures 9a?c give the results of the three synthesized brightness temperatures in Fig.6 processed by the BG method to 25-km grid cells.Figures 9d?f,meanwhile,show the errors between the gridded results and the true values.It can be seen from Figs.9d?f that the grid cell error of the BG method is smaller.The error in Fig.9d is also smaller than that in Fig.9f,which indicates that the BG method’s effectiveness is limited when adjacent pixels are badly overlapped.When Petty and Bennartz (2017) matched the FOVs of GMI,they also found that the 89-GHz channels were badly overlapped in the cross-track direction,and the synthetic FOV fitted to the target effective FOV was poor in that direction.

    In order to quantify these errors,the mean,error variance,and correlation coefficients for two methods’ results were calculated,separately.As Table 1 shows,the three cases’ mean values from using the direct method are almost the same as the true values,but the error variances are 3.00,3.62 and 4.99 K2,respectively.This error is greater than the inherent instrumental noise of the GMI’s 18.9-GHz channel.Combined with the information from Fig.7,these errors are mainly associated with regions displaying strong gradients in brightness temperature.The correlation coefficients between the direct method and true values are greater than 0.99,which indicates that the direct method’s results do not change the spatial distribution of the data.The error variances of the BG method are 0.22,0.42 and 0.49 K2,respectively,which are only about 10% of those of the direct method.Because the BG method uses pixels whose center is outside the study region when computing the result of the grid edge of the study area,the average values of the BG method differ significantly from the true values.The correlation coefficients of the BG method are also higher than those of the direct method.These results show that the BG method produces smaller error and higher correlation coefficients than the direct method,which indicates that the data quality of the BG method is higher than that of the direct method.

    Fig.9.As in Fig.7 but using the BG method.

    Table 1.Mean,variance of error,and correlation coefficient of the direct method and BG method.

    In order to examine the consistency and reliability of the results,the simulation experiment was repeated at another location,and the results were consistent and reliable.For brevity,the related figures can be found in the electronic supplementary material (ESM,Figs.S1?S3,Table S1).

    4.Example of GMI gridded data

    The BG method,which can greatly reduce the error of gridded results,was applied to two actual orbits.The orbit numbers of these orbits are 14 216 (example 1) and 11 386(example 2),respectively.The 18.9-GHz vertically (V) polarized channel brightness temperatures are shown in Figs.10a(example 1) and 10d (example 2).The direct method’s gridded results are shown in Figs.10b and 10e,while Figs.10c and 10f show the gridded results of the BG method.Orbit 11 386 passes over the ocean,while orbit 14 216 passes over the Tibetan Plateau.Figures 10a and b show that there are some small-scale regions of high brightness temperature in each area.Because 18.9-GHz is near the water absorption line at 22.235-GHz,the channel is influenced by precipitating water content.The high brightness temperature regions indicate a high water content in example 2.The results from using the direct and BG methods have large differences in these regions.To investigate the relationship between the differences in the results of the two methods and the brightness temperature gradient,we calculated the horizontal gradient field for each orbit’s data (example 1:Fig.11c;example 2:Fig.11d),as well as the horizontal distribution of the differences between the direct and BG methods (example 1:Fig.11a;example 2:Fig.11b).The figures show that the errors of the two methods are larger in regions of large brightness temperature gradient,which is similar to the results from the simulation experiment.This finding demonstrates that the BG method can better preserve information in these regions,such as mesoscale convective systems and cloud boundaries,and reduce the errors of gridded data.The quality of subsequent retrieval products could then also be improved.

    5.Conclusion

    Fig.10.Horizontal distribution of (a) GMI 18.9-GHz V brightness temperature,and the gridded result (0.25°×0.25°) of the(b) direct method and (c) BG method.Orbit number:14 216.(d?f) As in (a?c) but for orbit number 11 386.

    Fig.11.Differences between the gridded results of the direct and BG methods for two examples:(a) example 1,orbit number 14 216;(b) example 2,orbit number 11 386.Horizontal distribution of the horizontal gradient of 18.9-GHz V brightness temperature for two examples:(c) example 1,orbit number 14 216;(d) example 2,orbit number 11 386.

    In this paper,the sampling intervals of different regions of the GMI scan line have been analyzed.Three regions on the scan line were computed—namely,the “edge”,“subedge” and “center” regions—according to the shortest distance between adjacent pixels.The sampling intervals of the three regions were 6.5 km×6.0 km,11.5 km×6.0 km and 13.0 km×6.0 km,respectively,which were all smaller than the FOV of the 18.9-GHz channel of GMI (18.1 km×11.7 km).Therefore,the overlapping parts of the pixels of these three regions were different.Because there were dense overlapping effects between adjacent pixels,overlapping parts had greater weights than others.Furthermore,the FOV of pixels was similar to the size of the 0.25° grid cells;half of the FOVs sampled radiance that lay outside the grid cell,if pixels fell at the grid cell edge.The direct method,which only averages the pixels in the grid cells,caused large errors.

    Next,three different sampling intervals of pixels were used to measure the synthesized brightness temperature of a 1000 km×1000 km area in eastern China.The direct method and BG method were used to process the synthesized brightness temperature of pixels to 25 km×25 km grid cells.By comparing the results with the true values,it was found that the direct method causes non-negligible errors,especially in regions of large brightness temperature gradient.The error variances between the direct method’s results and the true values were 3.00,3.62 and 4.99 K2for the three sampling intervals,respectively.Such errors will affect the quality of subsequent retrieval products.To reduce this error,the BG method has been introduced in this paper.The BG method reduces the error of the gridding process by convolving or deconvolving the antenna patterns of pixels onto grid cells’ equivalent antenna patterns.By applying the BG method,the error variance for the three sampling intervals was reduced by about 90%,to 0.22,0.42 and 0.49 K2,respectively.The errors in high brightness temperature gradient regions were also reduced.

    Subsequently,the method was applied to real GMI orbits.The results showed that the differences between the results of the BG method and direct method were mainly distributed in the region of large gradients in brightness temperature,which was the same as for the results of the simulation experiment.The results also showed that the BG method can better preserve the information of mesoscale convective systems,cloud boundaries and other regions with large brightness temperature gradients,reduce the error of gridded data,and improve the quality of subsequent retrieval products.However,the BG method also has the disadvantage of high computational expense.As computer hardware evolves,though,the impact of this disadvantage should diminish.

    It has been noted that the gridding process can reduce the noise of the instrument itself (Njoku et al.,2003;Deeter and Vivekanandan,2006).Also,the FOV of the GMI’s 18.9-GHz channel (18.1 km×11.7 km) is similar to the size of a 0.25° grid cell (25 km),indicating that inversion using the gridded dataset is essentially equivalent to inversion using the pixel.The data processed by the BG method not only reduce the noise of the instrument,but also the error of the gridding process,which is more suitable to inversion.This will be addressed in future research.

    Acknowledgements.This study was supported by the National Key R&D Program of China (Grant Nos.2018YFC1507200 and 2017YFC1501402),the Second Tibetan Plateau Scientific Expedition and Research (STEP) program (Grant No.2019QZKK0104),an NSFC Project (Grant Nos.91837310,41675041,and 41620104009),the Key Research and Development Projects in Anhui Province (Grant No.201904a07020099),and CLIMATE-TPE (ID 32070) under the framework of the ESAMOST Dragon 4 program.

    Electronic supplementary material:Supplementary material is available in the online version of this article athttps://doi.org/10.1007/s00376-021-0358-x.

    国产成人午夜福利电影在线观看| 亚洲精华国产精华液的使用体验| 一级a做视频免费观看| 免费看a级黄色片| 日本免费在线观看一区| 久久精品国产亚洲av涩爱| 精品视频人人做人人爽| 久久久久久久大尺度免费视频| 王馨瑶露胸无遮挡在线观看| 欧美一区二区亚洲| 久久久国产一区二区| 99热国产这里只有精品6| 国产真实伦视频高清在线观看| 十八禁网站网址无遮挡 | 亚洲人成网站在线观看播放| 99热这里只有是精品在线观看| 亚洲欧美日韩卡通动漫| 亚州av有码| 久久久成人免费电影| 男女边摸边吃奶| 国产免费福利视频在线观看| 国产精品一区二区三区四区免费观看| 国产色婷婷99| 免费看av在线观看网站| 人妻少妇偷人精品九色| 亚洲最大成人手机在线| 天天躁日日操中文字幕| 国产精品三级大全| 亚洲一级一片aⅴ在线观看| av在线观看视频网站免费| 国产有黄有色有爽视频| 尤物成人国产欧美一区二区三区| 麻豆成人午夜福利视频| 91精品伊人久久大香线蕉| 欧美性感艳星| 国产精品三级大全| 中国三级夫妇交换| 亚洲欧美成人综合另类久久久| 人妻夜夜爽99麻豆av| 内射极品少妇av片p| 国产伦理片在线播放av一区| 亚洲人成网站在线观看播放| av线在线观看网站| 99热这里只有是精品在线观看| 午夜老司机福利剧场| 日韩强制内射视频| 国产国拍精品亚洲av在线观看| eeuss影院久久| 欧美bdsm另类| 五月伊人婷婷丁香| 91久久精品电影网| 免费观看在线日韩| 三级男女做爰猛烈吃奶摸视频| 人妻 亚洲 视频| 国产乱人偷精品视频| 最后的刺客免费高清国语| 精品久久久久久久末码| 欧美日韩视频高清一区二区三区二| 久久久久久久久大av| 国产精品秋霞免费鲁丝片| 久久久亚洲精品成人影院| 午夜福利网站1000一区二区三区| 夫妻性生交免费视频一级片| 五月开心婷婷网| 看黄色毛片网站| 亚洲国产精品999| 亚洲天堂国产精品一区在线| 麻豆乱淫一区二区| 久久国内精品自在自线图片| 一二三四中文在线观看免费高清| 中文字幕亚洲精品专区| 欧美成人a在线观看| 欧美最新免费一区二区三区| 涩涩av久久男人的天堂| 亚洲高清免费不卡视频| 能在线免费看毛片的网站| 久久韩国三级中文字幕| 天天躁日日操中文字幕| 一级毛片久久久久久久久女| 精品一区在线观看国产| 精品一区在线观看国产| 国产高清国产精品国产三级 | 中文天堂在线官网| 午夜福利网站1000一区二区三区| 在线天堂最新版资源| 精品国产露脸久久av麻豆| 夜夜看夜夜爽夜夜摸| 亚洲精品一二三| 亚洲在久久综合| 99久久九九国产精品国产免费| 别揉我奶头 嗯啊视频| 在线亚洲精品国产二区图片欧美 | 狂野欧美激情性bbbbbb| 男人和女人高潮做爰伦理| 春色校园在线视频观看| 18禁裸乳无遮挡免费网站照片| 卡戴珊不雅视频在线播放| 亚洲欧美清纯卡通| 在现免费观看毛片| 下体分泌物呈黄色| 在线免费观看不下载黄p国产| 五月天丁香电影| 亚洲精品日本国产第一区| 久久精品夜色国产| 日韩一区二区三区影片| 国产精品久久久久久久久免| 成年av动漫网址| 真实男女啪啪啪动态图| 女的被弄到高潮叫床怎么办| 国产欧美日韩一区二区三区在线 | 成年人午夜在线观看视频| 在线免费观看不下载黄p国产| 国产乱来视频区| 亚洲精品456在线播放app| 熟妇人妻不卡中文字幕| 高清在线视频一区二区三区| 免费黄频网站在线观看国产| 亚洲无线观看免费| 久久久欧美国产精品| 成人无遮挡网站| 中国三级夫妇交换| 在线观看国产h片| av国产久精品久网站免费入址| 国产 精品1| 亚洲精品一区蜜桃| 麻豆成人av视频| 永久免费av网站大全| 日韩制服骚丝袜av| 久久久久久伊人网av| 有码 亚洲区| 99久国产av精品国产电影| 大香蕉97超碰在线| 国产亚洲精品久久久com| 男人舔奶头视频| 国产伦在线观看视频一区| 熟女av电影| 国产伦在线观看视频一区| 嘟嘟电影网在线观看| 性色avwww在线观看| 舔av片在线| 18禁在线播放成人免费| 久久精品久久精品一区二区三区| 免费观看性生交大片5| 亚洲va在线va天堂va国产| 狂野欧美激情性bbbbbb| 亚洲国产成人一精品久久久| 老女人水多毛片| 亚洲欧美成人精品一区二区| 免费大片黄手机在线观看| 国产午夜福利久久久久久| 寂寞人妻少妇视频99o| 日产精品乱码卡一卡2卡三| 高清日韩中文字幕在线| av又黄又爽大尺度在线免费看| 成人午夜精彩视频在线观看| 成人综合一区亚洲| 国产精品国产av在线观看| 亚洲色图av天堂| 两个人的视频大全免费| 亚洲精品乱码久久久久久按摩| 亚洲av中文av极速乱| 男女下面进入的视频免费午夜| 国产精品秋霞免费鲁丝片| 欧美少妇被猛烈插入视频| 免费大片18禁| 亚洲欧美日韩卡通动漫| 午夜老司机福利剧场| 国产69精品久久久久777片| 亚洲国产精品国产精品| 国产日韩欧美亚洲二区| 久热这里只有精品99| 看十八女毛片水多多多| av网站免费在线观看视频| 又大又黄又爽视频免费| 五月玫瑰六月丁香| 国产av不卡久久| 色视频在线一区二区三区| 国内精品宾馆在线| 又爽又黄无遮挡网站| 久久久久久久久久久丰满| av专区在线播放| 久久鲁丝午夜福利片| 欧美xxⅹ黑人| 97人妻精品一区二区三区麻豆| 欧美bdsm另类| 国产男人的电影天堂91| 国产高清国产精品国产三级 | 国产成人精品福利久久| 午夜福利网站1000一区二区三区| 只有这里有精品99| 日本wwww免费看| 97在线视频观看| 久久韩国三级中文字幕| 久久人人爽人人片av| 亚洲成人一二三区av| 亚洲av成人精品一区久久| 国产 一区 欧美 日韩| 久久久精品免费免费高清| 一个人看的www免费观看视频| 亚洲欧美一区二区三区黑人 | 最近中文字幕2019免费版| 亚洲精品视频女| 成人国产麻豆网| 久久久久久国产a免费观看| 在线观看国产h片| 亚洲国产最新在线播放| 精品一区二区三区视频在线| 蜜桃亚洲精品一区二区三区| 精品一区二区三卡| 国产毛片a区久久久久| 肉色欧美久久久久久久蜜桃 | 99久久精品国产国产毛片| 亚洲精品第二区| 中文在线观看免费www的网站| 精品视频人人做人人爽| 精品酒店卫生间| 亚洲真实伦在线观看| 亚洲成人中文字幕在线播放| 最近最新中文字幕免费大全7| 国产亚洲最大av| av国产久精品久网站免费入址| 国产美女午夜福利| 国产精品蜜桃在线观看| 日韩欧美一区视频在线观看 | 国产精品三级大全| 丝袜美腿在线中文| 国产精品一区www在线观看| 亚洲美女搞黄在线观看| 亚洲,欧美,日韩| 国产免费福利视频在线观看| 黄片无遮挡物在线观看| 国产亚洲91精品色在线| 丝袜美腿在线中文| 三级国产精品片| 成人综合一区亚洲| 三级经典国产精品| 亚洲真实伦在线观看| 国产真实伦视频高清在线观看| 菩萨蛮人人尽说江南好唐韦庄| 狂野欧美激情性xxxx在线观看| 黄色怎么调成土黄色| 一级黄片播放器| 男女无遮挡免费网站观看| 精品国产露脸久久av麻豆| 精品久久久久久久久亚洲| 狂野欧美激情性xxxx在线观看| 欧美成人一区二区免费高清观看| 久久人人爽人人片av| av在线老鸭窝| 成年版毛片免费区| 免费高清在线观看视频在线观看| 搡女人真爽免费视频火全软件| 欧美国产精品一级二级三级 | 亚洲在线观看片| 免费av观看视频| 国产成人精品久久久久久| 精品熟女少妇av免费看| av专区在线播放| 亚洲精华国产精华液的使用体验| 美女视频免费永久观看网站| 国产伦精品一区二区三区四那| 国产高清不卡午夜福利| 日本色播在线视频| 丰满乱子伦码专区| 男插女下体视频免费在线播放| 男人狂女人下面高潮的视频| 国产日韩欧美在线精品| 成人亚洲欧美一区二区av| 五月伊人婷婷丁香| 高清av免费在线| 男插女下体视频免费在线播放| 蜜桃亚洲精品一区二区三区| 另类亚洲欧美激情| 纵有疾风起免费观看全集完整版| 日本色播在线视频| 久久影院123| 人妻少妇偷人精品九色| 涩涩av久久男人的天堂| 人人妻人人澡人人爽人人夜夜| 婷婷色麻豆天堂久久| 国产成人精品婷婷| 一本久久精品| 久久久久久久国产电影| 精品少妇黑人巨大在线播放| 一区二区三区精品91| 国产成人精品婷婷| 国产在视频线精品| 99热这里只有精品一区| 久久精品国产亚洲av涩爱| av播播在线观看一区| 久久综合国产亚洲精品| 日韩一区二区三区影片| 精品视频人人做人人爽| 一级毛片黄色毛片免费观看视频| 亚洲精品乱码久久久久久按摩| 亚洲av.av天堂| 亚洲,一卡二卡三卡| 美女高潮的动态| 80岁老熟妇乱子伦牲交| 亚洲av日韩在线播放| 亚洲av.av天堂| 国产免费视频播放在线视频| 女人被狂操c到高潮| 99久久精品一区二区三区| av福利片在线观看| 婷婷色麻豆天堂久久| 精品久久久久久久久av| 黄色视频在线播放观看不卡| 美女被艹到高潮喷水动态| 内射极品少妇av片p| 午夜福利高清视频| 精品国产乱码久久久久久小说| 成人漫画全彩无遮挡| 国产 精品1| 国产淫语在线视频| 波野结衣二区三区在线| 欧美bdsm另类| 69av精品久久久久久| 成人二区视频| 久久亚洲国产成人精品v| 纵有疾风起免费观看全集完整版| 国产精品偷伦视频观看了| 深爱激情五月婷婷| 日韩一区二区三区影片| 亚洲av中文字字幕乱码综合| 精品国产露脸久久av麻豆| 一区二区av电影网| 丝袜美腿在线中文| 最近2019中文字幕mv第一页| 日韩视频在线欧美| 在线a可以看的网站| 亚洲成人中文字幕在线播放| 久久影院123| 成年免费大片在线观看| 欧美高清成人免费视频www| 汤姆久久久久久久影院中文字幕| 激情 狠狠 欧美| 亚洲久久久久久中文字幕| 国产精品国产av在线观看| 精品久久久久久久人妻蜜臀av| 国产女主播在线喷水免费视频网站| 少妇被粗大猛烈的视频| 婷婷色综合www| 国产伦理片在线播放av一区| 久久精品久久久久久噜噜老黄| 高清日韩中文字幕在线| 成人鲁丝片一二三区免费| 久久鲁丝午夜福利片| 国产黄a三级三级三级人| 男人爽女人下面视频在线观看| 看黄色毛片网站| 久久99蜜桃精品久久| 久久午夜福利片| 可以在线观看毛片的网站| 日本午夜av视频| 亚洲av免费高清在线观看| 国产精品人妻久久久久久| 视频区图区小说| 国产精品99久久99久久久不卡 | 国产有黄有色有爽视频| 久久精品综合一区二区三区| 免费看av在线观看网站| 国产精品久久久久久久电影| 日韩强制内射视频| 伊人久久国产一区二区| 神马国产精品三级电影在线观看| 啦啦啦中文免费视频观看日本| 久久精品国产鲁丝片午夜精品| 身体一侧抽搐| 欧美zozozo另类| 国产色婷婷99| 少妇人妻一区二区三区视频| 亚洲精品影视一区二区三区av| 尤物成人国产欧美一区二区三区| 高清欧美精品videossex| 亚洲婷婷狠狠爱综合网| 老师上课跳d突然被开到最大视频| 亚洲成色77777| 亚洲国产欧美在线一区| 熟女人妻精品中文字幕| 国产在线男女| 日韩伦理黄色片| 爱豆传媒免费全集在线观看| 欧美性猛交╳xxx乱大交人| 久久精品夜色国产| 伊人久久精品亚洲午夜| 一边亲一边摸免费视频| 中文字幕制服av| 亚洲av中文av极速乱| 国产亚洲av嫩草精品影院| 日韩视频在线欧美| 亚洲激情五月婷婷啪啪| 自拍欧美九色日韩亚洲蝌蚪91 | 青青草视频在线视频观看| 久久精品久久久久久久性| 日韩,欧美,国产一区二区三区| 在线观看一区二区三区| 久久国内精品自在自线图片| 看黄色毛片网站| 日韩中字成人| 日本欧美国产在线视频| 熟女av电影| 亚洲伊人久久精品综合| 黄色视频在线播放观看不卡| 免费看不卡的av| 九九久久精品国产亚洲av麻豆| 性色avwww在线观看| 国产成人a区在线观看| 免费高清在线观看视频在线观看| 亚洲av一区综合| 又爽又黄a免费视频| 在线观看av片永久免费下载| 在线观看美女被高潮喷水网站| 18禁动态无遮挡网站| 制服丝袜香蕉在线| 亚洲av电影在线观看一区二区三区 | 嫩草影院入口| 国产大屁股一区二区在线视频| 91久久精品电影网| 综合色丁香网| 国产日韩欧美在线精品| 久久国产乱子免费精品| 亚洲va在线va天堂va国产| 人人妻人人看人人澡| 精品人妻视频免费看| 亚洲va在线va天堂va国产| 国产在线男女| 久久久成人免费电影| 国产在视频线精品| 国产精品人妻久久久久久| 精品国产乱码久久久久久小说| 久久午夜福利片| 午夜免费男女啪啪视频观看| 国产 一区精品| 最近手机中文字幕大全| 亚洲无线观看免费| 视频区图区小说| 乱系列少妇在线播放| 久久久久国产网址| 国产亚洲午夜精品一区二区久久 | 白带黄色成豆腐渣| 国产精品秋霞免费鲁丝片| 国产中年淑女户外野战色| 国产白丝娇喘喷水9色精品| 亚洲自偷自拍三级| 韩国高清视频一区二区三区| 99视频精品全部免费 在线| 日韩强制内射视频| 精品亚洲乱码少妇综合久久| 中文天堂在线官网| 免费观看无遮挡的男女| 亚洲成人中文字幕在线播放| 麻豆成人av视频| 深夜a级毛片| 丝袜脚勾引网站| 99久久九九国产精品国产免费| 亚洲精品色激情综合| 少妇猛男粗大的猛烈进出视频 | 日韩欧美精品v在线| 99久久精品国产国产毛片| 国产伦在线观看视频一区| 亚洲美女视频黄频| 国产精品久久久久久精品电影小说 | 久久久色成人| 亚洲精华国产精华液的使用体验| 国产一区二区在线观看日韩| 久久久久网色| 国产成人freesex在线| 国产成人福利小说| 特大巨黑吊av在线直播| 1000部很黄的大片| 午夜亚洲福利在线播放| 大又大粗又爽又黄少妇毛片口| 少妇被粗大猛烈的视频| 建设人人有责人人尽责人人享有的 | 国产精品久久久久久精品古装| 亚洲经典国产精华液单| 天美传媒精品一区二区| 久久久欧美国产精品| xxx大片免费视频| 免费不卡的大黄色大毛片视频在线观看| 亚洲va在线va天堂va国产| 亚洲,欧美,日韩| 一级毛片电影观看| tube8黄色片| 亚洲欧美日韩无卡精品| 日本一二三区视频观看| 亚洲欧美成人精品一区二区| 一本—道久久a久久精品蜜桃钙片 精品乱码久久久久久99久播 | 毛片女人毛片| 天天一区二区日本电影三级| 在线观看三级黄色| 亚洲美女搞黄在线观看| 国产黄频视频在线观看| 国产精品嫩草影院av在线观看| 一个人看视频在线观看www免费| 亚洲av免费在线观看| 热re99久久精品国产66热6| 大话2 男鬼变身卡| 18禁在线无遮挡免费观看视频| 国产伦精品一区二区三区四那| 黄片wwwwww| 久久99热6这里只有精品| 永久网站在线| videossex国产| 国产伦在线观看视频一区| 菩萨蛮人人尽说江南好唐韦庄| 免费在线观看成人毛片| 欧美区成人在线视频| 69av精品久久久久久| 久久久久久久久久成人| 有码 亚洲区| 最近最新中文字幕大全电影3| 街头女战士在线观看网站| 水蜜桃什么品种好| 欧美最新免费一区二区三区| 国产高清有码在线观看视频| 一级毛片aaaaaa免费看小| 国产精品一及| 如何舔出高潮| 蜜臀久久99精品久久宅男| 久久久久网色| 九九在线视频观看精品| 色视频在线一区二区三区| 精品久久久久久久人妻蜜臀av| 一本—道久久a久久精品蜜桃钙片 精品乱码久久久久久99久播 | 久久精品国产a三级三级三级| 欧美亚洲 丝袜 人妻 在线| 人妻少妇偷人精品九色| av.在线天堂| 男男h啪啪无遮挡| 亚洲在久久综合| 久久影院123| 赤兔流量卡办理| 国产成人精品婷婷| 国产综合懂色| 国产真实伦视频高清在线观看| 久久久久精品性色| 天堂中文最新版在线下载 | 麻豆乱淫一区二区| 午夜福利网站1000一区二区三区| 亚洲婷婷狠狠爱综合网| 国产精品99久久久久久久久| 一级a做视频免费观看| 精品一区二区三卡| 一本—道久久a久久精品蜜桃钙片 精品乱码久久久久久99久播 | 国产欧美日韩一区二区三区在线 | 欧美精品人与动牲交sv欧美| 高清欧美精品videossex| 国产有黄有色有爽视频| av国产久精品久网站免费入址| 99re6热这里在线精品视频| 超碰av人人做人人爽久久| 亚洲精品一二三| 搡老乐熟女国产| 80岁老熟妇乱子伦牲交| 啦啦啦在线观看免费高清www| 中文字幕亚洲精品专区| 亚洲精品视频女| 丰满少妇做爰视频| 麻豆成人av视频| 欧美bdsm另类| 能在线免费看毛片的网站| 国产一区有黄有色的免费视频| 亚洲精品日韩av片在线观看| 少妇熟女欧美另类| 国产精品福利在线免费观看| 亚洲国产最新在线播放| 伦精品一区二区三区| 高清av免费在线| 亚洲人成网站在线播| 亚洲av日韩在线播放| 日韩成人av中文字幕在线观看| 久久精品国产a三级三级三级| 欧美日韩视频精品一区| 欧美丝袜亚洲另类| 边亲边吃奶的免费视频| 午夜福利视频1000在线观看| 国产成人精品福利久久| 亚洲欧美日韩卡通动漫| 欧美+日韩+精品| 免费观看性生交大片5| 九草在线视频观看| 一个人观看的视频www高清免费观看| 欧美日本视频| a级毛片免费高清观看在线播放| 黄色欧美视频在线观看| 国产成人精品婷婷| 久久这里有精品视频免费| 九色成人免费人妻av| 久久久久久久亚洲中文字幕| av又黄又爽大尺度在线免费看| av在线app专区| 一本色道久久久久久精品综合| 中文在线观看免费www的网站| 草草在线视频免费看| 国产精品蜜桃在线观看| 99热国产这里只有精品6| 精品亚洲乱码少妇综合久久| 国产成人免费观看mmmm| 又黄又爽又刺激的免费视频.| 在线观看一区二区三区| 九九爱精品视频在线观看| 亚洲怡红院男人天堂| 乱码一卡2卡4卡精品| 亚洲自偷自拍三级| 免费人成在线观看视频色| 国产精品三级大全| 国产成人一区二区在线| 超碰av人人做人人爽久久| 人妻制服诱惑在线中文字幕| 国产免费视频播放在线视频| 九九在线视频观看精品| 人人妻人人澡人人爽人人夜夜| 一级毛片久久久久久久久女|