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

    Development and Testing of the GRAPES Regional Ensemble-3DVAR Hybrid Data Assimilation System

    2015-12-12 12:57:01CHENLianlii陳良呂CHENJing陳靜XUEJishan薛紀善andXIAYu夏宇
    Journal of Meteorological Research 2015年6期
    關鍵詞:陳靜

    CHEN Lianlii(陳良呂), CHEN Jing(陳靜),XUE Jishan(薛紀善),and XIA Yu(夏宇)

    1 Chongqing Institute of Meteorological Sciences,Chongqing 401147

    2 Chinese Academy of Meteorological Sciences,Beijing 100081

    3 Numerical Weather Prediction Center,China Meteorological Administration,Beijing 100081

    4 Chengdu University of Information Technology,Chengdu 610225

    Development and Testing of the GRAPES Regional Ensemble-3DVAR Hybrid Data Assimilation System

    CHEN Lianlii1,2(陳良呂), CHEN Jing3*(陳靜),XUE Jishan2(薛紀善),and XIA Yu4(夏宇)

    1 Chongqing Institute of Meteorological Sciences,Chongqing 401147

    2 Chinese Academy of Meteorological Sciences,Beijing 100081

    3 Numerical Weather Prediction Center,China Meteorological Administration,Beijing 100081

    4 Chengdu University of Information Technology,Chengdu 610225

    Based on the GRAPES(Global/Regional Assimilation and Prediction System)regional ensemble prediction system and 3DVAR(three-dimensional variational)data assimilation system,which are implemented operationally at the Numerical Weather Prediction Center of the China Meteorological Administration,an ensemble-based 3DVAR(En-3DVAR)hybrid data assimilation system for GRAPES?Meso(the regional mesoscale numerical prediction system of GRAPES)was developed by using the extended control variable technique to implement a hybrid background error covariance that combines the climatological covariance and ensemble-estimated covariance.Considering the problems of the ensemble-based data assimilation part of the system,including the reduction in the degree of geostrophic balance between variables,and the non-smooth analysis increment and its obviously smaller size compared with the 3DVAR data assimilation, corresponding measures were taken to optimize and ameliorate the system.Accordingly,a single pressure observation ensemble-based data assimilation experiment was conducted to ensure that the ensemble-based data assimilation part of the system is correct and reasonable.A number of localization-scale sensitivity tests of the ensemble-based data assimilation were also conducted to determine the most appropriate localization scale.Then,a number of hybrid data assimilation experiments were carried out.The results showed that it was most appropriate to set the weight factor of the ensemble-estimated covariance in the experiments to be 0.8.Compared with the 3DVAR data assimilation,the geopotential height forecast of the hybrid data assimilation experiments improved very little,but the wind forecast improved slightly at each forecast time,especially over 300 hPa.Overall,the hybrid data assimilation demonstrates some advantages over the 3DVAR data assimilation.

    GRAPES,GRAPES?Meso,hybrid data assimilation,regional ensemble prediction,extended control variable

    1.Introduction

    A three-dimensional variational(3DVAR)data assimilation system usually requires climate-statisticsbased,uniform,and isotropic background error covariance information.However,the real background error covariance changes with the weather.This conflict limits the quality of data assimilation to some extent, and further affects the ability to make improvements to the quality of weather forecasts.In contrast,in a four-dimensional variational(4DVAR)data assimilation system,the background error covariance is adjusted implicitly in the assimilation window.Thus, 4DVAR can,to a certain extent,capture the flowdependent characteristic.However,we still need to provide the climatological background error at the initial time in 4DVAR.The Ensemble Kalman filter (EnKF)data assimilation can produce a flow-depend-

    ent background error covariance by extracting the information from the samples of the short-range ensemble forecast(Xue et al.,2013). Unlike variational assimilation,however,the EnKF is incapable of dealing with large amounts of irregular observational data. Accordingly,researchers have developed hybrid schemes by applying the flow-dependent background error covariance derived from the ensemble forecast into the variational assimilation structure. This method is called ensemble-based variational hybrid data assimilation(En-VAR)(Lorenc, 2003;Buehner,2005;Wang et al.,2008a,b).En-VAR not only retains the flow-dependent characteristic of the background error covariance,but also takes advantage of variational assimilation in dealing with numerous irregular observational data.

    Over the past decade,hybrid data assimilation has emerged internationally as a hot topic in the data assimilation research community.Researchers abroad have carried out many studies in hybrid data assimilation on the basis of 3DVAR or 4DVAR.Lorenc(2003) first proposed to assimilate the ensemble-estimated background error covariance matrix into the structure of variational data assimilation by extending the control variable,and also demonstrated the theoretical possibility of this method.Wang et al.(2007) proved,theoretically,that the two hybrid schemes, i.e.,by extending the control variable and by weighting two types of background error covariance matrix directly,are equivalent.Many researchers(Buehner, 2005;Wang et al.,2008a,b;Wang,2010;Hamill et al.,2011)have performed hybrid data assimilation experiments based on different models.Their results illustrate that the application of ensemble forecast information in the variational structure can improve the shortcomings of climatological background error covariance,as well as the quality of the assimilation and forecast.This improvement is particularly obvious in regions where observations are lacking.Besides,hybrid data assimilation shows an obvious advantage in tropical storm forecasting(Wang,2011).Recently, basic and applied researches in hybrid data assimilation have become increasingly common.Internationally,a few advanced numerical weather prediction centers have already applied hybrid data assimilation into their operational systems(Clayton et al.,2013;Wang et al.,2013),and this application shows a good operational perspective.Research on hybrid data assimilation techniques in China is also quite mature,with many researchers having already carried out numerous studies on the basis of 4DVAR(Qiu and Chou,2006; Liu et al.,2008;Wang et al.,2010).Recently,a number of Chinese research centers have applied the WRF hybrid data assimilation system to analyses of tropical storm cases and in forecast experiments(Huang et al., 2014),with a few having already applied this system in an operational setting(personal communication).

    The NumericalWeather Prediction Center (NWPC)of the China Meteorological Administration (CMA)has been working on applying the EnKF to GRAPES(Global/Regional Assimilation Prediction System)since approximately 2000(Liu and Xue,2005; Zhuang et al.,2011a,b).However,thus far,the approach has been applied in only a few case studies, which were different from applying it to daily operational numerical weather prediction that involves vast quantities of observational assimilation.Meanwhile,research on hybrid data assimilation for operational numerical prediction has not yet started.Currently,the data assimilation used for operational numerical prediction in China still employs the single variational method(Liu and Xue,2014).Most researchers,however,tend to see the adoption of the hybrid data assimilation method as the next step for the GRAPES assimilation system.The NWPC has adopted the GRAPES regional 3DVAR system and GRAPES regional ensemble forecast system,laying a solid foundation for research on a regional ensemblebased 3DVAR(En-3DVAR)hybrid data assimilation system for GRAPES.Therefore,it is necessary to begin exploring hybrid data assimilation under the framework of GRAPES in preparation for the buildup of the operational GRAPES data assimilation system in future.

    In this study,the GRAPES regional En-3DVAR hybrid data assimilation system,based on the existing GRAPES regional 3DVAR system and GRAPES regional ensemble forecast system,was developed by using the widely adopted extended control variable method.A series of processing schemes for solving the

    problems caused by the ensemble prediction perturbation samples,which were obtained from the GRAPES regional ensemble prediction system,were developed. In addition,some preliminary experiments were carried out to further understand the remaining problems,which may serve as a reference for the future development of the GRAPES regional data assimilation system.

    2.The GRAPES regional 3DVAR data assimilation system and ensemble prediction system

    2.1The GRAPES regional 3DVAR data assimilation system

    This study was based on the existing GRAPES regional 3DVAR data assimilation system,which is implemented operationally at the NWPC of the CMA (Ma et al,2009b).The state variables of the system comprise the wind variable(zonal wind and meridional wind),mass variable(dimensionless pressure or potential temperature),and the humidity variable(relative humidity or specific humidity).The variational problem is to find the analysis fields xaby minimizing the cost function given below:

    Here,Bcand R are the background and observation error covariance matrices,respectively;H is the operator mapping from the model space to the observation space;y contains the observations;x is the real analysis field;and xbis the background field.The incremental analysis method is used in this system,where x′=xa?xbdenotes the analysis increment.Then, the cost function can be modified to

    where d=H(xb)?y denotes the innovation vector.

    The degrees of freedom of the climatological background error covariance matrix Bcis very large;generally,a value of 107can be achieved. In order to reduce the tremendous cost of minimization and make convergence easier,a preconditioning transform or control variable transform is used.First,the state variables are transformed to independent state variables through physical variable transform and balance transform;the independent state variables are the stream function(ψ),velocity potential(χ),unbalanced mass variable(πu/thu),and humidity variable (rh/q).Therefore,the background error covariance matrix becomes a block-diagonal matrix Bu,which is composed of the independent state variables,and can thus reduce the computation costs considerably.The square root of matrix Buis then used to precondition the cost function and define the control variable v for minimizing the cost function.In the actual calculation,the analysis increment x′is obtained through a series of transformation processes applied to the control variable v:

    where Uhis the horizontal precondition transform process accomplished via a recursive filter;Uvis the vertical precondition transform process accomplished via a projection from the eigenvector of the empirical orthogonal functions(EOFs)of the vertical component of the background error onto model levels;εbis a diagonal matrix denoting the root-mean-square errors of the independent state variables at the grid point; UKis the balance transform process used to transform πu/thuto π/th;and UPis the physical variable transform process used to transform ψ and χ to model state variables(U and V).

    Substituting Eq.(3)into Eq.(2),the cost function can be modified to

    The gradient of the cost function can be obtained after setting the first-order derivative of Eq.(4)with respect to v equal to zero:

    Equations(4)and(5)are the basic equations of the GRAPES regional 3DVAR data assimilation system.In this system,the Limited Broyden-Fletche-Goldfarb-Shanno method(LBFGS;initials of names of the scientists who invented this method)is used to solve the minimization problem of the control variables,and thus can obtain the optimal control variable v,then the analysis increment,through Eq.(4),and add the analysis increment to the background fields to obtain the analysis field.This is the basic calculation process of the GRAPES regional 3DVAR data assimilation system.

    2.2GRAPES regional ensemble prediction system

    The NWPC of the CMA began to implement the GRAPES regional ensemble prediction system operationally in May 2014(Zhang et al.,2014).This system is based on GRAPES?Meso(the regional mesoscale numerical prediction system of GRAPES),whose resolution is 0.15°×0.15°,and the model domain is set to(15°–65°N,70°–145°E).The initial conditions of the system are perturbed by using the ensemble transform Kalman filter(ETKF)scheme(Ma et al.,2009a). Each ensemble member uses different physical process parameterization schemes,and the background fields and lateral boundary conditions for each member are provided by the T639 global ensemble prediction system.This system consists of 15 members,including a control run and 14 perturbed ensemble members, and is initiated at 0600,1200,1800,and 0000 UTC each day.For each initiation time,the system provides 6-h forecast perturbations for the next ETKF cycle;specifically,for the 1200 and 0000 UTC initiation times,the model integrates to 72 h to provide ensemble prediction products.

    3.Design of the hybrid data assimilation system

    3.1 Theoretical design

    In the GRAPES regional En-3DVAR hybrid data assimilation system,flow-dependent ensemble covariance is incorporated in the variational minimization via the extended control variable method(Clayton et al.,2013).The normalized ensemble perturbation fieldsare taken from the 6-h forecast fields of the GRAPES regional ensemble prediction system.Here,K is the ensemble size,xkis the kth ensemble forecast,andˉx is the ensemble mean.The flow-dependent ensemble-estimated background error covarianceis then given by:

    Due to the limited ensemble size,the main problem of the ensemble-estimated background error covarianceis its low rank(at most,K)and the presence of sample error,which may lead to unreasonable correlation between variables in practical application. Based on this,a localization matrix C,whose size is the same asis designed.can be localized through the Schur product of C and itself,and then the ensemble-estimated background error covariance is replaced as Eq.(7):

    The hybrid data assimilation scheme can be achieved by using the method of directly combining the climatological and ensemble-estimated background error covariance:

    where x′is the analysis increment,andandare two factors that define the weights placed on the climatological and ensemble-estimated background error covariance. To conserve the total backgrounderror variance,andare constrained byWhenthe analysis returns to a 3DVAR analysis scheme,and whenthe analysis is mathematically equivalent to an ensemble-based analysis scheme.However, this scheme is difficult to implement directly in the existing 3DVAR scheme.Therefore,the extended control variable method is used to accomplish the scheme.

    In theory,it can be proven that the method of combining the climatological and ensemble-estimated background error covariance is equivalent to the method of combining the 3DVAR analysis increment

    and the ensemble-estimated analysis increment(Wang et al.,2007).

    In the hybrid data assimilation system,the analysis increment x′is the sum of two terms,defined as:

    The analysis increment x′can be obtained by minimizing the following cost function:

    Compared with the 3DVAR cost function,the sum of the J1and J2terms in Eq.(10)replaces the usual background term.J1is the traditional 3DVAR background term associated with the climatological background error covariance.The new term,J2,is associated with the ensemble-estimated background error covariance,and it is the same as J1.On this basis, J1is preconditioned by the square root of the climatological background error covariance Bc,and J2is preconditioned by the square root of the localization matrix C.

    For the first term,J1,as described in Section 3.1, the precondition process is composed of a series of transformations,as Eq.(3)shows.The variable v is the traditional 3DVAR control variable,which contains ψ,χ,πu/thu,and rh/q.The dimension of v is the number of mode variables multiplied by the number of three-dimensional(3D)grid points.

    For the second term,J2,since the localization matrix C is used for accomplishing the localization of the ensemble-estimated background error covariance, in this system,the same as the precondition process of the 3DVAR analysis scheme,the square root of matrix C is modeled by using a horizontal recursive filter and vertical EOF method,as Eq.(11)shows:

    Substituting the precondition process into Eq. (10),the analysis increment x′and the cost function of the hybrid data assimilation can be modified to

    and

    On this basis,the gradient of each control variable can be obtained:

    and

    Equations(12)–(15)are the basic equations of the GRAPES regional En-3DVAR hybrid data assimilation system.To accomplish this,the control variables are extended and the related calculation processes of the extended control variables are added based on the

    existing 3DVAR data assimilation system.The result is just a few changes to the procedural structure of the 3DVAR data assimilation system.

    3.2 Experimental design

    Considering the uncertainty in the humidity part of the ensemble perturbation fields provided by the GRAPES regional ensemble prediction system,the ensemble-estimated background error covariance information on humidity is not incorporated into the hybrid data assimilation system;the humidity part of the ensemble perturbation fields is all set to error.While calculating the eventual humidity analysis increment of the hybrid data assimilation system,the weight of the climatological background error covariance is set to one and the weight of the ensemble-estimated background error covariance is set to zero.Thus,the eventual humidity analysis increment is still the same as the 3DVAR data assimilation.

    In this study,following the design of the hybrid data assimilation system,a number of preliminary experiments at 0000 UTC 5 July 2014 were carried out. The forecast domain and resolution were the same as in the GRAPES?Meso mentioned above.The 6-h forecast of the T639 global prediction system was used as background fields,and the lateral boundary conditions were also provided by the T639 global prediction system.The default physical process parameterization schemes of GRAPES?Meso(the WSM6 microphysical process parameterization scheme,the RRTM longwave radiation scheme,the Dudhia shortwave radiation scheme,the Monin-Obukhov planetary boundary layer scheme,the Noah land surface process scheme, the MRF planetary boundary layer scheme,and the Betts-Miller-Janjic cumulus convection scheme)were used in the experiments.The FNL(resolution:1°×1°) analysis data of the NCEP were interpolated to the grid points of GRAPES?Meso,and then used as the real atmospheric state to verify the forecast of the experiments in this study.

    4.Optimization of ensemble-based data assimilation

    After setting up the above hybrid data assimilation structure,the correctness and rationality of the ensemble-based assimilation(referring to the analysis process while setting the weight of the ensembleestimated background error covariance to 1 in the hybrid data assimilation system)should be considered first.The weight coefficient of the system was defined as

    and a number of single-point idealized tests and real observation data assimilation tests were then conducted. The results showed that the degree of geostrophic balance between the model variables was reduced,and the analysis increment was non-smooth and obviously smaller than the 3DVAR data assimilation.To solve the above problems,their causes were carefully explored and corresponding measures were taken to optimize and ameliorate the system.

    4.1Optimization of the balance between analysis variables

    The results of the single-point idealized tests indicated that the degree of geostrophic balance between the variables is reduced in the ensemble-based data assimilation.There are probably two reasons for this. One is that some high frequency oscillations without meteorological significance were not removed effectively during the short-term integration procedures in the GRAPES regional model.This would lead to a worse degree of geostrophic balance in the ensemble prediction samples obtained from the 6-h ensemble forecast data.The other is that the localization transformation process of extended control variables may also reduce the degree of geostrophic balance between the variables.

    According to the method of Clayton et al.(2013), physical transformation and balance transformation were applied in the ensemble prediction perturbation samples(u,v,π),similar to the transformations applied in the 3DVAR data assimilation.The extended control variables were localized mainly for independent non-equilibrium variables(ψ,χ,πu).The atmospheric variables(u,v,π)were then obtained through the corresponding physical transformation and balance transformation algorithm,thus ensuring the bal-

    ance characteristics in the system,to a certain degree.

    Correspondingly,the equations for calculating the analysis increment and extended control variable gradient were then modified as follows:

    Here,the variable of the ensemble prediction samplesare ψ,χ,and non-equilibrium πu(unbalanced dimensionless pressure).andare the physical transformation and balance transformation algorithm, respectively.

    4.2 Smooth filtering of ensemble forecast perturbations

    Theensembleprediction samplesapplied in the hybrid data assimilation system are from the GRAPES regional ensemble prediction system(details in Section 2.2).The results of the single-point idealized tests and real observation data assimilation tests show that the analysis increment of the ensemble-based data assimilation is not smooth,especially for the wind variables,and there are much smaller scale increments.The reasons might be due to problems in the ensemble prediction perturbation samples.The ensemble prediction perturbations applied in the system are not smooth and there is much noise (Figs.1a1and 1b1),and the analysis increment of the ensemble-based data assimilation is actually the linear combination of the ensemble prediction perturbations. Thus,the obtained analysis increments of the ensemble prediction perturbations are generally disordered and non-smooth.To solve these problems,a typical five-grid filtering smooth algorithm is used in the system to smooth the ensemble prediction perturbation fields so that their noise is filtered and removed,and a more smooth and regular perturbation and analysis increment field is obtained.The detailed filtering formula is as follows: where s is the filtering coefficient(0–1).The larger the value of s,the more obvious the filtering effects. Meanwhile,the more filtering processes,the more obvious the filtering effects. After many contrasting tests,the coefficient of s was set to 0.7.The idealized effects after several filtering processes are shown in Figs.1a2and 1b2.

    4.3Amplification of ensemble prediction perturbations

    The results of the comparisons between the singlepoint idealized and real observation data assimilation tests indicate that the analysis increment of the ensemble-based data assimilation is obviously smaller than that of 3DVAR,especially in the upper layers of the atmosphere(at the 25th layer and above,the closer to the model top,the more obvious the difference). Plus,the analysis increment of dimensionless pressure is much more obvious than that of wind(U,V).This might be due to the smaller ensemble spread of the ensemble prediction samples.The ensemble spread is used as the estimation of the background error variance.If a low ensemble spread occurs,it means that the background error variance is estimated lower in the ensemble-estimated background error covariance matrix.This would lead to considerable background field fitting and less observation information to be counted during the data assimilation process.Finally,a smaller analysis increment of the ensemble-based data assimilation would occur.

    The vertical distribution characteristics for the 6-h ensemble prediction spread(dashed line)and RMSE of the 6-h ensemble prediction mean(solid line) at 1800 UTC 4 July 2014 are shown in Figs.2a–c. Similar vertical distribution characteristics are apparent for each variable spread and RMSE.Furthermore, each ensemble spread is smaller than the RMSE,especially in higher layers;meanwhile,for the 25th layer and above,the closer to the model top,the more obvious the difference.Figure 2d is the same as Figs.2a–c except that it shows the ratio of the ensemble spread and RMSE of the ensemble mean(U:solid line;V: dashed line;PI(dimensional pressure π):dotted line), which changes between 0 and 1.The closer to 1,the

    smaller the difference between the spread and RMSE.The differences between the spread and RMSE for dimensionless pressure are larger than those of wind(U,V).These characteristics correspond to those of the analysis increment of the ensemble-based data assimilation.

    Fig.1.The 20th layer ensemble perturbation fields for winds of the 7th ensemble member before and after smooth filtering:(a1)U wind before smooth filtering,(b1)V wind before smooth filtering,(a2)U wind after smooth filtering, and(b2)V wind after smooth filtering.

    In view of the above problems,the ensemble prediction perturbations were multiplied by an amplified coefficient in the ensemble-based data assimilation system to relieve the problem of a smaller ensemble spread.The principle of deciding the coefficient was that it should lead to a corresponding and appropriate analysis increment in both the ensemble-based data assimilation and the 3DVAR data assimilation.According to the prior vertical distribution characteristics for the ensemble spread and RMSE,the coefficient for the perturbations of dimensionless pressure should be larger than that of wind(U,V),but it should not be so large that it damages the balance constraints in the system.It should be larger for higher layers than lower ones,and it could be linearly increased from the 25th layer upwards.

    The results from a number of sensitivity experiments with the amplified coefficients showed high sensitivity of the ensemble-based data assimilation analysis to the coefficients.The coefficients should be selected appropriately,and not too big.Ultimately,the coefficients were selected as follows:below the 25th layer,it is 1.5 for perturbations of dimensionless pressure,1.3 for perturbations of wind(U,V);and for the 25th layer and above,it increases to 0.1 for perturbations of pressure and 0.05 for perturbations of wind (U,V)within each layer.

    5.Preliminary tests

    Fig.2.Vertical distribution characteristics for the ensemble prediction spread(dashed line)and RMSE of the ensemble prediction mean(solid line):(a)U wind;(b)V wind;(c)PI;and(d)the ratio of the ensemble spread to RMSE of the ensemble mean.

    Based on the above work,some preliminary tests (single-point idealized tests,sensitivity tests on the localization scale of the ensemble-based data assimilation,and real observation data assimilation tests) were implemented for the case introduced in Section 3.2,which occurred at 0000 UTC 5 July 2014.

    5.1 Idealized single-point tests

    After optimization of the ensemble-based data assimilation in the hybrid data assimilation system, single idealized tests for the 3DVAR,ensemble-based and hybrid assimilation(the weighted coefficients for the climatological and ensemble-estimated background error covariance were both 0.5)were implemented based on the updated constructed hybrid data assimilation system,to test the correctness and rationality of the hybrid data assimilation system.The horizontal and vertical localization scales were set the same as the horizontal and vertical correlation scales of the 3DVAR data assimilation.

    Figures 3 and 4 show the results of the horizontal and vertical analysis increment of the idealized tests for single-point pressure observation(33°N,100°E; 10th layer).Overall,the distribution characteristics for the ensemble-based data assimilation and 3DVAR data assimilation were similar.The wind and pressure variables satisfied the constraints of geostrophic balance.The analysis increment of the 3DVAR data assimilation was quite smooth and regular,and the analysis increment of the ensemble-based data assimilation was slightly irregular and presented a certain degree of anisotropy.After localization of the ensemble-estimated background error covariance matrix was conducted in the system,no obvious noise was apparent in the analysis increment of the ensemblebased data assimilation. The values of the analysis increment for the dimensionless pressure of the ensemble-based data assimilation were equal to those of 3DVAR.However,the analysis increments of wind (U,V)were obviously larger than those of 3DVAR. This might have been due to the differences between the background error variances,which were obtained from the climatological and ensemble-estimated error covariance matrix,respectively.In summary,under the condition of the same background error variance, the distribution characteristics of the analysis increment obtained from the ensemble-based data assimilation and 3DVAR data assimilation were basically the same.This therefore indicates the correctness and rationality of the ensemble-based data assimilation in the hybrid data assimilation system.

    The values of the analysis increment obtained from the updated hybrid data assimilation (the weighted coefficients for the climatologicaland ensemble-estimated background error covariance were both 0.5)were between those of the 3DVAR data assimilation and ensemble-based data assimilation.The distribution characteristics of the increments obtained from the above three sources presented generally the same pattern.Their variables were constrained by the geostrophic balance and behaved with a certain degree of anisotropy.This might have been because the hybrid data assimilation“absorbed”the background error variance information from the short-range ensemble prediction samples.This further indicates the correctness and rationality of the hybrid data assimilation system.

    Fig.3.Horizontal analysis increments of the 3DVAR,ensemble-based,and hybrid data assimilation:(a1)PI,3DVAR; (b1)U wind,3DVAR;(c1)V wind,3DVAR;(a2)PI,ensemble-based;(b2)U wind,ensemble-based;(c2)V wind, ensemble-based;(a3)PI,hybrid;(b3)U wind,hybrid;and(c3)V wind,hybrid.

    5.2 Sensitivity tests on the localization scales of the ensemble-based data assimilation

    The horizontal localization of the ensemble-based data assimilation is implemented through recursive filtering,and the scale of recursive filtering determines the scale of localization.The scale of horizontal recursive filtering for the 3DVAR data assimilation is 500 km.The vertical localization of the ensemble-based data assimilation is implemented through EOF decomposition,and the local scale is determined by a defined parameter.For ease of comparison,the horizontal and vertical localization scales of the ensemble-based data assimilation were set the same as the correlation scales of the 3DVAR data assimilation.Further investigation is still needed in terms of how sensitive the ensemblebased data assimilation is to the localization scales.

    A set of sensitivity tests on the localization scales were devised,in which the vertical localization scales of the ensemble-based data assimilation were set the same as the vertical correlation scales of the 3DVAR data assimilation.The horizontal localization scale of the ensemble-based data assimilation(control test) was set to 500 km,which is the same as the horizontal correlation scales of 3DVAR.In addition,four other sets of comparison tests were performed,in which the localization scales were 625,750,875,and 1000 km, respectively.All radiosonde data in the region were

    assimilated,and the ensemble-based assimilation data were numerically integrated.The total RMSE of the 3D grid field for the geopotential height and wind 12-h predictions were compared,thus enabling us to explore the sensitivity of the ensemble-based data assimilation to the horizontal localization scales.The test results are shown in Table 1.

    Fig.4.As in Fig.3,but for the vertical analysis increment.

    From Table 1,it is shown that the RMSE of the 12-h prediction of wind(U,V)and geopotential height exhibited a certain degree of sensitivity to the horizontal localization scale.A larger horizontal localization scale improved the test performance a little when compared with the control test;and,the larger the localization scale,the more obvious the improvement. The most improvement was achieved if the localization scale was set to 1000 km.Therefore,the localization scale was determined as 1000 km in the following tests.

    After the horizontal localization scale was determined,similar sensitivity tests for the vertical localization scale were implemented based on several sets

    of different vertical localization scale parameters.The test results showed that the ensemble-based data assimilation was not sensitive to the vertical localization scale.Therefore,the vertical localization scale parameter was set the same as the vertical correlation scale of the 3DVAR data assimilation.

    Table 1.Sensitivity test results(RMSE values)for the horizontal localization scale of the ensemble-based data assimilation

    5.3 Real observation data assimilation and prediction tests

    The above single-point idealized tests and localization scale sensitivity tests for the ensemble-based data assimilation indicated the correctness and rationality of the hybrid data assimilation system.With the optimal localization scales determined,assimilation and prediction tests by using real observational data were implemented in order to further analyze and verify the performance of the hybrid data assimilation system.Four sets of combinations with different weighting coefficients were applied in the hybrid data assimilation.All radiosonde data in the region were assimilated,and the numerical integration was carried out based on the hybrid data assimilation results.Based on a comparison of the characteristics of the total RMSE obtained from the 12-and 24-h predicted 3D grid fields for geopotential height and wind, the sensitivity of the hybrid data assimilation to the weighting coefficients was explored.The test results are shown in Table 2.

    From Table 2,it is shown that,based on the total RMSE of the 12-and 24-h prediction of wind (U,V)and geopotential height,there were obvious improvements in the hybrid assimilation compared with 3DVAR.The larger the weighting of the coefficient of the ensemble-estimated background error covariance,the more obvious the improvement.The most improvement was achieved when the coefficient of the ensemble-estimated background error covariance was defined as 0.8.The geopotential height of the ensemble-based data assimilation presented obvious negative effects,but there was a slight improvement in the wind(U,V).Overall,the optimal weighted coefficient for the ensemble-estimated background error covariance was determined as 0.8.

    Figure 5 shows the analysis increment fields of wind(U,V)at the 10th layer for the 3DVAR data assimilation and hybrid data assimilation(the weighting coefficient of the ensemble-estimated background error covariance was 0.8).It is seen that the analysis increment values from the hybrid data assimilation and 3DVAR data assimilation were similar.The distribution and locations of positive and negative centers were also similar.The maximum values from the hybrid assimilation were relatively smaller than those from the 3DVAR assimilation. The analysis increments from the 3DVAR assimilation were quite smooth and presented characteristics of uniformity and isotropy. However,the analysis increments from the hybrid assimilation presented characteristics of inhomogeneity and anisotropy.The hybrid assimilation basically retained the large-scale increment information of the 3DVAR data assimilation,and added some small-scale increment information.This might have been because it was integrated into some ensemble-estimated background error covariance information.The analysis increments for other model layers and other variables also presented the same characteristics(figure omitted).

    Figure 6 shows the 6–72-h RMSE for winds of the hybrid and 3DVAR data assimilation,in which the hybrid data assimilation performed slightly betterthan the 3DVAR data assimilation.The RMSE of the meridional wind(U)fields improved by 2.76%,generally,with the best one for the 6-h prediction(improvement of 7.39%).Whereas,the RMSE of the zonal wind(V)fields improved by 2.63%,generally,with the best one for the 6-h prediction(improvement of 5.77%).There was not much improvement forthe RMSE of geopotential height compared with that for the 3DVAR assimilation(figure omitted).

    Table 2.Sensitivity test results(RMSE values)for the weighting coefficient of the hybrid data assimilation

    Fig.5.The analysis increment for wind at the 10th level of the 3DVAR and hybrid data assimilation:(a1)U wind, 3DVAR;(b1)V wind,3DVAR;(a2)U wind,hybrid;and(b2)V wind,hybrid.

    Fig.6.The 6–72-h RMSE for winds of the hybrid and 3DVAR data assimilation:(a)U wind and(b)V wind.

    Figure 7 shows the vertical profiles of the 12-and

    24-h RMSEs for winds of the hybrid and 3DVAR data assimilation,in which the hybrid data assimilation performed slightly better than the 3DVAR assimilation for the 12-and 24-h wind prediction at all layers,especially at higher layers(300 hPa and above). For 300 hPa and above,the 12-h RMSE of the meridional wind(U)fields improved by 5.98%,generally, with the best one improving by 13.11%.Whereas, the 12-h RMSE of the zonal wind(V)fields improved by 6.04%,generally,with the best one improving by 12.45%.The 24-h RMSE of the meridional wind(U) fields improved by 2.43%,generally,with the best one improving by 5.11%.Whereas,the 24-h RMSE of the zonal wind(V)fields improved by 4.31%,generally, with the best one improving by 9.04%.There was a slight improvement for the 12-and 24-h RMSEs of geopotential height at very high layers(150 hPa and above),compared with that for the 3DVAR data assimilation.As for other layers,the 12-and 24-h RMSEs of geopotential height behaved with slight negative effects(figure omitted).

    Fig.7.Vertical profiles of 12-and 24-h RMSEs for winds of the hybrid(dashed line)and 3DVAR(solid line)data assimilation:(a)12-h forecast of U wind;(b)12-h forecast of V wind;(c)24-h forecast of U wind;and(d)24-h forecast of V wind.

    6.Summary and discussion

    This study performed the theoretical design of the GRAPES regional En-3DVAR hybrid data assimilation system by using the extended control variable method,and carried out preliminary experiments. Because the degree of geostrophic balance between variables was found reduced in the ensemble-based data assimilation,the system was optimized by introducing balance constraints between variables.It was also found that the analysis increment of the ensemblebased data assimilation was unsmooth,and this was solved through the five-point smoothing filtering of the ensemble perturbation fields.The ensemble perturbation fields were enlarged,thus relieving the problem caused by the low spread of ensemble perturbations. Then,a single pressure observation ensemble-based data assimilation experiment was conducted to ensure that the ensemble-based data assimilation is correct and reasonable.In addition,a number of localization scale sensitivity tests of the ensemble-based data assimilation were also conducted to determine the most appropriate localization scale.Finally,a number of hybrid data assimilation experiments were conducted. The results showed that the most appropriate setting for the weighting factor of the ensemble-estimated covariance was 0.8.Compared with the 3DVAR data assimilation,the geopotential height forecast in the hybrid data assimilation experiments improved very little,but the wind forecasts showed a slight improvement at each forecast time,especially over 300 hPa. Overall,the hybrid data assimilation has a number of advantages over the 3DVAR data assimilation.

    The following points discuss the problems and future development directions of the system:

    (1)In the experiments conducted in this study, the computational expense of the hybrid data assimilation was about five to eight times that of the 3DVAR data assimilation,and the computational expense increased the most for the ensemble-based data assimilation.Overall,the computational expense of the hybrid data assimilation increased substantially, but it was much less than the traditional EnKF data assimilation.

    (2)In the current system,the smoothing filter method for the ensemble perturbation fields is very simple,and an iterative smoothing filter may reduce the larger scale ensemble perturbation to some degree. In the future,a preferable smoothing filter method may be used in this system.

    (3)Because the spin-up time of GRAPES?Meso is long,some difficulties to solve this problem exist in the system.In the future,12-h ensemble perturbation fields will be used to replace the 6-h ensemble perturbation fields,which may increase the degree of geostrophic balance between the variables and reduce the noise of the following analysis increment.

    (4)Only one set of hybrid data assimilation experiments was carried out,based on which the appropriate inflation factors of the ensemble perturbations,the localization scale,and the weights of the climatological and ensemble-based background error covariance were selected.In the future,many more experiments will be carried out and analyzed to select the optimal sensitive factors mentioned above.There are some problems with the selection of inflation factors of the ensemble perturbations,and these problems interact with each other.Solving these problems and identifying the best inflation factors to reduce the error caused by the low spread of ensemble perturbations is a core issue of immediate priority in the future.In addition,selection of the localization scale is related to the scale of the weather situation(synoptic,meso, etc.)and the type of observation data.In this study, only radiosonde data were assimilated,and the resolution was just 15 km.The conclusions regarding the localization scale therefore have some limitations, and further in-depth studies related to the localization scale need to be carried out in the future.

    (5)In this study,the hybrid data assimilation experiment results showed that,in general,the hybrid data assimilation carries some advantages over the 3DVAR data assimilation.In the future,after further optimization of the system and identification of optimal sensitive parameters,many more experiments will be carried out to establish clearly whether or not the hybrid data assimilation is superior to the 3DVAR data assimilation.

    (6)There is a large gap between the current GRAPES regional ensemble prediction system and other ensemble prediction systems implemented operationally in advanced numerical weather prediction centers outside of China.In the future,studies on optimizing the GRAPES regional ensemble prediction system(e.g.,optimizing the ETKF initial condition perturbation,adding stochastic perturbations of the physical process tendency,etc.)will be carried out, to further improve the quality of the GRAPES hybrid data assimilation system.

    Buehner,M.,2005:Ensemble-derived stationary and flow-dependent background error covariances:Evaluation in a quasi-operational NWP setting.Quart. J.Roy.Meteor.Soc.,131,1013–1043.

    Clayton,A.M.,A.C.Lorenc,and D.M.Barker, 2013: Operational implementation of a hybrid ensemble/4D-Var global data assimilation system at the Met Office.Quart.J.Roy.Meteor.Soc., 139,1445–1461.

    Hamill,T.M.,J.S.Whitaker,D.T.Kleist,et al.,2011: Predictions of 2010's tropical cyclones using the GFS and ensemble-based data assimilation methods.Mon.Wea.Rev.,139,3243–3247.

    Huang Jiangping,Liu Jianwen,Dong Peiming,et al., 2014:Application of the hybrid data assimilation system for Typhoon Haiou.J.Trop.Meteor.,30, 1113–1118.(in Chinese)

    Lorenc,A.C.,2003:The potential of the ensemble Kalman filter for NWP—A comparison with 4DVar.Quart.J.Roy.Meteor.Soc.,129,3183–3204.

    Liu Chengsi and Xue Jishan,2005:The ensemble Kalman filter theory and method development.J.Trop.Meteor.,21,628–633.(in Chinese)

    Liu,C.S.,Q.N.Xiao,and B.Wang,2008:An ensemblebased four-dimensional variational data assimilation scheme.Part I:Technical formulation and preliminary test.Mon.Wea.Rev.,136,3363–3373.

    Liu Yan and Xue Jishan,2014:Assimilation of global navigation satellite radio occultation observations in GRAPES:Operational implementation.J.Meteor.Res.,28,1061–1074.

    Ma Xulin,Xue Jishan,and Lu Weisong,2009:Study on ETKF-based initial perturbation scheme for

    GRAPES global ensemble prediction.Acta.Meteor.Sinica,23,562–574.

    Ma Xulin,Zhuang Zhaorong,Xue Jishan,et al.,2009b: Developments of the three-dimensional variational data assimilation system for the nonhydrostatic GRAPES.Acta Meteor.Sinica,23,725–737.

    Qiu,C.,and J.Chou,2006:Four-dimensional data assimilation method based on SVD:Theoretical aspect. Theor.Appl.Climatol.,83,51–57.

    Wang Bin,Liu Juanjuan,Wang Shudong,et al.,2010: An economical approach to four-dimensional variational data.Adv.Atmos.Sci.,27,715–727.

    Wang,X.G.,2010:Incorporating ensemble covariance in the gridpoint statistical interpolation variational minimization:A mathematical framework.Mon. Wea.Rev.,138,2990–2995.

    Wang,X.G.,2011:Application of the WRF hybrid ETKF-3DVAR data assimilation system for hurricane track forecasts.Wea.Forecasting,26,868–884.

    Wang,X.G.,C.Snyder,and T.M.Hamill,2007:On the theoretical equivalence of differently proposed ensemble/3DVAR hybrid analysis schemes.Mon. Wea.Rev.,135,222–227.

    Wang,X.G.,D.M.Barker,C.Snyder,et al.,2008a:A hybrid ETKF-3DVAR data assimilation scheme for the WRF model.Part I:Observing system simulation experiment.Mon.Wea.Rev.,136,5116–5131.

    Wang,X.G.,D.M.Barker,C.Snyder,et al.,2008b:A hybrid ETKF-3DVAR data assimilation scheme for the WRF model.Part II:Real observation experiments.Mon.Wea.Rev.,136,5132–5147.

    Wang,X.G.,D.Parrish,D.Kleist,et al.,2013:GSI 3DVar-based ensemble-variational hybrid data assimilation for NCEP global forecast system:Singleresolution experiments.Mon.Wea.Rev.,141, 4098–4117.

    Xue Ming and Dong Jili,2013:Assimilating best track minimum sea level pressure data together with Doppler radar data using and ensemble Kalman filter for Hurricane Ike(2008)at a cloud-resolving resolution.Acta Meteor.Sinica,27,379–399.

    Zhang Hanbin,Chen Jing,Zhi Xiefei,et al.,2014:Study on the application of GRAPES regional ensemble prediction system.Meteor.Mon.,40,1076–1087. (in Chinese)

    Zhuang Zhaorong,Xue Jishan,and Li Xingliang,2011a: The GRAPES ensemble Kalman filter data assimilation system. Part I:Design and its tentative experiment.Acta Meteor.Sinica,69,620–630.(in Chinese)

    Zhuang Zhaorong,Xue Jishan,and Li Xingliang,2011b: The GRAPES ensemble Kalman filter data assimilation system.Part II:Regional analysis and ensemble prediction.Acta Meteor.Sinica,69,860–871.(in Chinese)

    :Chen Liangl¨u,Chen Jing,Xue Jishan,et al.,2015:Development and testing of the GRAPES regional ensemble-3DVAR hybrid data assimilation system.J.Meteor.Res.,29(6),981–996,

    10.1007/s13351-015-5021-y.

    Supported by the National Natural Science Foundation of China(91437113 and 41275111)and China Meteorological Administration Special Public Welfare Research Fund(GYHY201506005).

    ?Corresponding author:chenj@cma.gov.cn.

    ?The Chinese Meteorological Society and Springer-Verlag Berlin Heidelberg 2015

    (Received May 2,2015;in final form October 26,2015)

    猜你喜歡
    陳靜
    醉人的春夜
    氫氧燃料電池演示實驗的改進
    化學教學(2022年4期)2022-05-07 22:54:09
    指向深度學習的高中生物實驗教學探究
    高中學生數(shù)學解題能力培養(yǎng)探究
    醉人的春夜
    醉人的春夜
    陳靜:寧靜致遠,香飄南美
    華人時刊(2020年13期)2020-09-25 08:21:36
    呂乘乘、陳靜、黃峰、尹臻作品
    醉人的春夜
    太陽
    99热这里只有精品一区| 国产精品人妻久久久久久| 大陆偷拍与自拍| 国产极品天堂在线| 国产国拍精品亚洲av在线观看| 热re99久久国产66热| 欧美亚洲日本最大视频资源| av专区在线播放| 天天影视国产精品| 欧美最新免费一区二区三区| 免费播放大片免费观看视频在线观看| 美女福利国产在线| 超碰97精品在线观看| 午夜视频国产福利| 超色免费av| 夫妻午夜视频| 午夜激情av网站| 日本欧美视频一区| 亚洲第一区二区三区不卡| 九九久久精品国产亚洲av麻豆| 久久97久久精品| 妹子高潮喷水视频| 亚洲无线观看免费| 一区在线观看完整版| 丰满迷人的少妇在线观看| 久久精品夜色国产| 成年女人在线观看亚洲视频| 一级毛片 在线播放| 在线观看免费高清a一片| 国产精品久久久久久久久免| 国产精品一区二区三区四区免费观看| 精品国产一区二区久久| 婷婷色综合www| 伊人久久国产一区二区| 久久女婷五月综合色啪小说| 我的女老师完整版在线观看| 亚洲av免费高清在线观看| a 毛片基地| 国产精品蜜桃在线观看| 热re99久久精品国产66热6| 国产精品成人在线| 寂寞人妻少妇视频99o| 日韩熟女老妇一区二区性免费视频| 免费av不卡在线播放| 日日摸夜夜添夜夜添av毛片| 不卡视频在线观看欧美| 国产高清不卡午夜福利| 狠狠婷婷综合久久久久久88av| 亚洲av日韩在线播放| av线在线观看网站| 精品人妻一区二区三区麻豆| 成年女人在线观看亚洲视频| 亚洲久久久国产精品| 人人妻人人爽人人添夜夜欢视频| 日韩欧美一区视频在线观看| 日韩成人伦理影院| 91精品一卡2卡3卡4卡| 久久久久国产网址| 久久精品久久久久久噜噜老黄| 久久精品久久精品一区二区三区| 国产深夜福利视频在线观看| 91精品一卡2卡3卡4卡| 三级国产精品欧美在线观看| 黄片无遮挡物在线观看| 亚洲欧美精品自产自拍| 大话2 男鬼变身卡| 亚洲av二区三区四区| 视频在线观看一区二区三区| 免费观看a级毛片全部| 亚洲欧洲国产日韩| 久久ye,这里只有精品| 国产精品嫩草影院av在线观看| 亚洲图色成人| 国产成人一区二区在线| 大陆偷拍与自拍| 日韩伦理黄色片| 特大巨黑吊av在线直播| 丰满少妇做爰视频| 国产成人免费无遮挡视频| 全区人妻精品视频| 伦理电影大哥的女人| 高清午夜精品一区二区三区| 少妇被粗大猛烈的视频| 国产精品国产三级国产av玫瑰| 9色porny在线观看| av不卡在线播放| 九色亚洲精品在线播放| 久久久久久久国产电影| 亚洲第一av免费看| 亚洲久久久国产精品| 欧美日韩综合久久久久久| 亚洲精品自拍成人| 国产精品一国产av| 色吧在线观看| 制服人妻中文乱码| 国语对白做爰xxxⅹ性视频网站| 99久久精品国产国产毛片| 人妻少妇偷人精品九色| 欧美成人午夜免费资源| 在线播放无遮挡| 永久网站在线| 黄片无遮挡物在线观看| 国产成人免费观看mmmm| 亚洲欧洲日产国产| 丝袜美足系列| 国产无遮挡羞羞视频在线观看| 午夜福利视频在线观看免费| h视频一区二区三区| 久久女婷五月综合色啪小说| 我的老师免费观看完整版| 亚洲欧洲国产日韩| 国产高清国产精品国产三级| 99九九在线精品视频| 性色av一级| 免费黄频网站在线观看国产| 国产精品99久久久久久久久| 国产免费福利视频在线观看| 欧美日韩视频精品一区| 一级a做视频免费观看| 成人无遮挡网站| 欧美3d第一页| 国产在线视频一区二区| 桃花免费在线播放| 曰老女人黄片| 看免费成人av毛片| 国产成人免费观看mmmm| 人妻系列 视频| 国产一区二区三区综合在线观看 | 熟女电影av网| 成人漫画全彩无遮挡| 飞空精品影院首页| 美女国产高潮福利片在线看| 亚洲精华国产精华液的使用体验| 成人国产麻豆网| 国产亚洲精品久久久com| 秋霞伦理黄片| 日韩伦理黄色片| 亚洲欧美清纯卡通| 人妻一区二区av| 久久精品国产自在天天线| 久久午夜福利片| 欧美一级a爱片免费观看看| 久久免费观看电影| 18禁在线无遮挡免费观看视频| 免费少妇av软件| 久久免费观看电影| 久久 成人 亚洲| 久久久久久久亚洲中文字幕| 18在线观看网站| 国产亚洲精品第一综合不卡 | 搡老乐熟女国产| 韩国高清视频一区二区三区| 成人亚洲精品一区在线观看| 久久久精品94久久精品| 乱人伦中国视频| 国产成人精品无人区| av天堂久久9| 国产精品蜜桃在线观看| 国产免费一级a男人的天堂| 99九九线精品视频在线观看视频| 国产色爽女视频免费观看| 人体艺术视频欧美日本| 最黄视频免费看| 又粗又硬又长又爽又黄的视频| 99九九在线精品视频| 欧美成人午夜免费资源| 在线精品无人区一区二区三| 亚洲天堂av无毛| 免费观看性生交大片5| 亚洲不卡免费看| 欧美最新免费一区二区三区| 日韩不卡一区二区三区视频在线| 国精品久久久久久国模美| 免费观看的影片在线观看| 亚洲精品中文字幕在线视频| 大又大粗又爽又黄少妇毛片口| 老司机影院毛片| 在线天堂最新版资源| 在线观看免费日韩欧美大片 | 少妇人妻精品综合一区二区| 国产一区有黄有色的免费视频| 一个人看视频在线观看www免费| 日本黄大片高清| 国产伦精品一区二区三区视频9| 久久久久久久精品精品| 老司机亚洲免费影院| av线在线观看网站| 亚洲国产日韩一区二区| 亚洲精品亚洲一区二区| 国产成人91sexporn| 国产精品女同一区二区软件| 99热6这里只有精品| 考比视频在线观看| 亚洲精品亚洲一区二区| 视频区图区小说| 岛国毛片在线播放| 欧美三级亚洲精品| 男女啪啪激烈高潮av片| 一本色道久久久久久精品综合| 一二三四中文在线观看免费高清| 人妻少妇偷人精品九色| 中文字幕久久专区| 国产熟女欧美一区二区| av免费观看日本| 久久久久久久精品精品| 国产日韩欧美视频二区| 亚洲三级黄色毛片| 国产高清不卡午夜福利| 亚洲欧美成人综合另类久久久| 91精品国产国语对白视频| 欧美 日韩 精品 国产| 日韩视频在线欧美| 在线看a的网站| 亚洲国产欧美在线一区| 久久久久久人妻| 母亲3免费完整高清在线观看 | 国产成人精品无人区| 亚洲色图综合在线观看| 亚洲人成网站在线观看播放| 国产一区亚洲一区在线观看| 美女视频免费永久观看网站| 91午夜精品亚洲一区二区三区| 9色porny在线观看| 午夜影院在线不卡| 日韩中字成人| 美女国产高潮福利片在线看| 国产精品久久久久久久电影| 91精品国产国语对白视频| 久久久久久伊人网av| 少妇人妻精品综合一区二区| 日本-黄色视频高清免费观看| 视频中文字幕在线观看| 欧美精品高潮呻吟av久久| 在线观看免费日韩欧美大片 | 欧美日韩成人在线一区二区| 18禁在线无遮挡免费观看视频| 日韩欧美一区视频在线观看| 国产一区亚洲一区在线观看| av黄色大香蕉| av福利片在线| 丰满饥渴人妻一区二区三| 大话2 男鬼变身卡| 性色avwww在线观看| 伦精品一区二区三区| 大片电影免费在线观看免费| 国产色爽女视频免费观看| 汤姆久久久久久久影院中文字幕| 久久久久网色| 亚洲成人av在线免费| 欧美性感艳星| 少妇被粗大的猛进出69影院 | 激情五月婷婷亚洲| 国产免费又黄又爽又色| 国产成人91sexporn| 婷婷成人精品国产| 免费人妻精品一区二区三区视频| 97精品久久久久久久久久精品| 国产高清不卡午夜福利| 美女福利国产在线| 岛国毛片在线播放| 国产精品女同一区二区软件| 久久久欧美国产精品| 成人国语在线视频| 亚洲精品一区蜜桃| 日韩伦理黄色片| 国产日韩欧美视频二区| 亚洲精品亚洲一区二区| 日韩亚洲欧美综合| 日本免费在线观看一区| 亚洲婷婷狠狠爱综合网| 久久精品国产亚洲av涩爱| 97在线视频观看| 涩涩av久久男人的天堂| 亚洲国产精品专区欧美| 91精品国产九色| 亚洲精品日韩在线中文字幕| 99久国产av精品国产电影| av有码第一页| 中文字幕av电影在线播放| a级片在线免费高清观看视频| 亚州av有码| 中文字幕免费在线视频6| 国产精品国产三级专区第一集| 色婷婷久久久亚洲欧美| 国产精品久久久久久久电影| 最后的刺客免费高清国语| 黄色配什么色好看| 久久99热这里只频精品6学生| 日韩一区二区视频免费看| 一本色道久久久久久精品综合| 久久久久久久国产电影| 精品国产乱码久久久久久小说| 秋霞在线观看毛片| 国产免费福利视频在线观看| 国内精品宾馆在线| 精品熟女少妇av免费看| 少妇猛男粗大的猛烈进出视频| 人人妻人人爽人人添夜夜欢视频| 9色porny在线观看| 人人妻人人爽人人添夜夜欢视频| 欧美性感艳星| 人人妻人人爽人人添夜夜欢视频| 久久综合国产亚洲精品| 人人妻人人爽人人添夜夜欢视频| 欧美性感艳星| 伦精品一区二区三区| av一本久久久久| 男女无遮挡免费网站观看| 99热网站在线观看| 十八禁高潮呻吟视频| 免费观看a级毛片全部| 麻豆精品久久久久久蜜桃| 国产精品 国内视频| 日韩亚洲欧美综合| 永久免费av网站大全| 成人无遮挡网站| 国产精品.久久久| 午夜免费鲁丝| 国产精品99久久久久久久久| 午夜免费鲁丝| 看免费成人av毛片| 日日爽夜夜爽网站| 女人久久www免费人成看片| 视频中文字幕在线观看| 免费看光身美女| 18禁在线播放成人免费| 国产伦理片在线播放av一区| 国产亚洲最大av| 国产淫语在线视频| 九色成人免费人妻av| 免费少妇av软件| 亚洲精品乱久久久久久| 精品久久久噜噜| 一级毛片 在线播放| 亚洲国产精品国产精品| 一级毛片 在线播放| 伊人久久国产一区二区| 日韩大片免费观看网站| 久久久久久久亚洲中文字幕| 又粗又硬又长又爽又黄的视频| 亚洲av成人精品一区久久| 99热网站在线观看| av在线老鸭窝| 久久人人爽人人爽人人片va| 青春草视频在线免费观看| 一本色道久久久久久精品综合| 涩涩av久久男人的天堂| 国产一区亚洲一区在线观看| 免费人成在线观看视频色| 久久久久视频综合| 久久狼人影院| 久久久国产一区二区| 精品久久久噜噜| 国产精品人妻久久久影院| 少妇精品久久久久久久| 高清av免费在线| 日韩视频在线欧美| 欧美国产精品一级二级三级| 99久久精品一区二区三区| 精品久久久久久久久av| 国产精品人妻久久久久久| 男女国产视频网站| 一本一本久久a久久精品综合妖精 国产伦在线观看视频一区 | 欧美日韩一区二区视频在线观看视频在线| 我要看黄色一级片免费的| 美女脱内裤让男人舔精品视频| 最后的刺客免费高清国语| 人体艺术视频欧美日本| 免费日韩欧美在线观看| av网站免费在线观看视频| a级毛片黄视频| 亚洲精品aⅴ在线观看| 欧美日韩av久久| 久久久久久久久久久久大奶| 久久人人爽人人片av| freevideosex欧美| 国产一级毛片在线| 国产欧美日韩综合在线一区二区| 丝瓜视频免费看黄片| 极品人妻少妇av视频| av在线app专区| 国产成人91sexporn| 国产亚洲av片在线观看秒播厂| 我的女老师完整版在线观看| 亚洲三级黄色毛片| 亚洲欧洲精品一区二区精品久久久 | 国产色婷婷99| 国产精品三级大全| 午夜视频国产福利| 制服人妻中文乱码| 九草在线视频观看| 国模一区二区三区四区视频| 大话2 男鬼变身卡| 免费黄色在线免费观看| 欧美日韩在线观看h| 欧美成人精品欧美一级黄| 免费播放大片免费观看视频在线观看| 精品亚洲成a人片在线观看| 各种免费的搞黄视频| 亚洲av成人精品一区久久| 国产国拍精品亚洲av在线观看| 一区二区日韩欧美中文字幕 | 热re99久久国产66热| 国产av精品麻豆| 99久久精品国产国产毛片| 亚洲精品aⅴ在线观看| 狂野欧美白嫩少妇大欣赏| 超色免费av| 91精品国产九色| 日韩精品有码人妻一区| 啦啦啦啦在线视频资源| 中文欧美无线码| 国产精品国产三级国产专区5o| 欧美日韩在线观看h| 亚洲国产成人一精品久久久| av免费在线看不卡| 青春草视频在线免费观看| 一级毛片电影观看| 欧美+日韩+精品| 久久 成人 亚洲| 精品国产一区二区三区久久久樱花| 99久久精品国产国产毛片| 欧美精品一区二区大全| 建设人人有责人人尽责人人享有的| 中文欧美无线码| 三上悠亚av全集在线观看| 亚洲综合色惰| 日韩一区二区视频免费看| 成人18禁高潮啪啪吃奶动态图 | 激情五月婷婷亚洲| 国产一区二区三区综合在线观看 | 亚洲国产欧美日韩在线播放| 少妇人妻久久综合中文| 国产淫语在线视频| 天堂8中文在线网| 黄色一级大片看看| 插阴视频在线观看视频| 交换朋友夫妻互换小说| 国产亚洲一区二区精品| 26uuu在线亚洲综合色| 亚洲精品一区蜜桃| 亚洲国产欧美在线一区| 久久亚洲国产成人精品v| 免费观看性生交大片5| 亚洲成人av在线免费| 五月天丁香电影| 最新的欧美精品一区二区| 亚洲成人一二三区av| 国产伦精品一区二区三区视频9| 亚洲国产毛片av蜜桃av| 精品一区在线观看国产| 亚洲国产日韩一区二区| 免费观看的影片在线观看| 国产永久视频网站| 欧美成人精品欧美一级黄| 色婷婷av一区二区三区视频| 一级毛片黄色毛片免费观看视频| 精品人妻熟女av久视频| 五月伊人婷婷丁香| 乱码一卡2卡4卡精品| 色哟哟·www| 中国美白少妇内射xxxbb| 免费观看无遮挡的男女| 在线看a的网站| 国产视频内射| 少妇人妻久久综合中文| 国产av一区二区精品久久| 亚洲av.av天堂| 亚洲国产毛片av蜜桃av| 久久午夜综合久久蜜桃| 波野结衣二区三区在线| 亚洲国产色片| 欧美日韩在线观看h| 久久精品国产鲁丝片午夜精品| av又黄又爽大尺度在线免费看| 成人午夜精彩视频在线观看| 一级二级三级毛片免费看| 成人亚洲精品一区在线观看| 下体分泌物呈黄色| 亚洲天堂av无毛| 日本91视频免费播放| 国产片内射在线| 亚洲国产日韩一区二区| 一本一本综合久久| 91国产中文字幕| 亚洲国产精品国产精品| 激情五月婷婷亚洲| 亚洲经典国产精华液单| a级毛片黄视频| 制服诱惑二区| 欧美精品高潮呻吟av久久| 青春草亚洲视频在线观看| 亚洲国产日韩一区二区| 五月天丁香电影| 国产精品无大码| 新久久久久国产一级毛片| 丰满乱子伦码专区| 亚洲成色77777| 黄色视频在线播放观看不卡| 欧美bdsm另类| 一级毛片aaaaaa免费看小| 欧美日韩视频高清一区二区三区二| 亚洲成人av在线免费| 免费黄色在线免费观看| 日本wwww免费看| 日韩一区二区视频免费看| 国精品久久久久久国模美| 亚洲欧美精品自产自拍| 999精品在线视频| 黑人高潮一二区| 高清视频免费观看一区二区| 美女cb高潮喷水在线观看| 亚洲一区二区三区欧美精品| 久久精品国产亚洲av天美| 一边亲一边摸免费视频| 人妻制服诱惑在线中文字幕| 国产片特级美女逼逼视频| 免费观看a级毛片全部| 国产日韩欧美在线精品| 精品少妇久久久久久888优播| 国产av码专区亚洲av| 秋霞在线观看毛片| 亚洲av男天堂| 亚洲性久久影院| 我要看黄色一级片免费的| 黄片无遮挡物在线观看| 亚洲av成人精品一区久久| 综合色丁香网| 满18在线观看网站| 中文字幕制服av| 国产精品免费大片| 久久久久久久久久久免费av| 最近的中文字幕免费完整| 亚洲怡红院男人天堂| 欧美变态另类bdsm刘玥| 一区二区三区免费毛片| 欧美亚洲 丝袜 人妻 在线| 极品少妇高潮喷水抽搐| 搡老乐熟女国产| 九九在线视频观看精品| 丰满饥渴人妻一区二区三| 18禁裸乳无遮挡动漫免费视频| 人人妻人人爽人人添夜夜欢视频| 欧美 亚洲 国产 日韩一| 国产 一区精品| 黄色欧美视频在线观看| 亚洲精品久久午夜乱码| 两个人的视频大全免费| 免费黄网站久久成人精品| 久久久精品94久久精品| 99精国产麻豆久久婷婷| 国产69精品久久久久777片| 日韩中文字幕视频在线看片| 日韩,欧美,国产一区二区三区| 久久久久网色| 亚洲欧洲国产日韩| 亚洲精华国产精华液的使用体验| 日韩免费高清中文字幕av| 在线观看免费日韩欧美大片 | 18禁在线播放成人免费| 国产成人精品久久久久久| 欧美日韩亚洲高清精品| 亚洲欧洲日产国产| 色网站视频免费| 欧美3d第一页| 午夜激情久久久久久久| 看十八女毛片水多多多| 日韩一区二区三区影片| 国产精品人妻久久久影院| 免费大片黄手机在线观看| 最近2019中文字幕mv第一页| 精品国产一区二区久久| 亚洲精品一二三| 日韩免费高清中文字幕av| 亚洲天堂av无毛| 久热这里只有精品99| 日韩欧美精品免费久久| 日本av免费视频播放| 日韩精品免费视频一区二区三区 | 午夜激情av网站| 日本vs欧美在线观看视频| 伦精品一区二区三区| 亚洲国产色片| 少妇人妻久久综合中文| 国产精品99久久久久久久久| 丝袜脚勾引网站| 亚洲伊人久久精品综合| 欧美人与善性xxx| 久久久久精品性色| 国产伦精品一区二区三区视频9| 中国美白少妇内射xxxbb| 色婷婷久久久亚洲欧美| 精品国产一区二区久久| 亚洲精华国产精华液的使用体验| 午夜福利影视在线免费观看| 欧美日韩精品成人综合77777| 久久这里有精品视频免费| 王馨瑶露胸无遮挡在线观看| 91久久精品电影网| 啦啦啦啦在线视频资源| 日本猛色少妇xxxxx猛交久久| 夫妻性生交免费视频一级片| 超碰97精品在线观看| 多毛熟女@视频| 成人亚洲欧美一区二区av| 亚洲精品中文字幕在线视频| 多毛熟女@视频| 99精国产麻豆久久婷婷| 晚上一个人看的免费电影| 亚洲欧洲精品一区二区精品久久久 | 18禁裸乳无遮挡动漫免费视频| 欧美激情极品国产一区二区三区 | 国产成人午夜福利电影在线观看| 亚洲四区av| 亚洲国产av新网站|