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

    基于CMA模式體系的京津冀地區(qū)復(fù)雜地形下冬季的精細(xì)化地面要素多模式集成預(yù)報(bào)研究*

    2022-02-14 12:28:44張玉濤齊倩倩王遠(yuǎn)哲王大鵬
    氣象 2022年12期
    關(guān)鍵詞:格點(diǎn)方根時(shí)效

    佟 華 張玉濤 齊倩倩 王遠(yuǎn)哲 王大鵬

    1 中國(guó)氣象局地球系統(tǒng)數(shù)值預(yù)報(bào)中心,北京 100081 2 災(zāi)害天氣國(guó)家重點(diǎn)實(shí)驗(yàn)室,北京 100081

    提 要: 基于CMA模式體系的四個(gè)模式(CMA-GFS、CMA-REPS、CMA-MESO 3 km、CMA-MESO 1 km)和2020年12月1日至2021年3月15日的近地面要素2 m溫度、10 m風(fēng)速、2 m相對(duì)濕度預(yù)報(bào),對(duì)京津冀地區(qū)復(fù)雜地形下冬季誤差訂正后的各要素進(jìn)行基于貝葉斯模型平均(BMA)方法的多模式集成試驗(yàn)。結(jié)果表明,每個(gè)模式各要素誤差訂正后的均方根誤差都有明顯的減小。BMA方法多模式集成后預(yù)報(bào)效果優(yōu)于每一個(gè)參加模式,2 m溫度BMA預(yù)報(bào)較幾個(gè)模式原始誤差的改進(jìn)在0.5~1.4℃,均方根誤差減少了20%~40%,10 m風(fēng)速和2 m相對(duì)濕度的均方根誤差分別減少了12%~45%和25%~35%。各要素均方根誤差水平分布表明不同要素在不同地形高度的地區(qū)誤差分布明顯不同,此方法使得京津冀地區(qū)的誤差顯著減小。此外,BMA預(yù)報(bào)的概率分布情況可定量地預(yù)測(cè)各要素的不確定性。

    引 言

    隨著數(shù)值預(yù)報(bào)技術(shù)不斷進(jìn)步,準(zhǔn)確率不斷提高,數(shù)值預(yù)報(bào)結(jié)果已經(jīng)成為預(yù)報(bào)員的重要參考。而來(lái)自不同研發(fā)單位、不同分辨率、不同預(yù)報(bào)特點(diǎn)的海量數(shù)值預(yù)報(bào)產(chǎn)品對(duì)預(yù)報(bào)員的高效使用產(chǎn)生了較大挑戰(zhàn),迫切需要對(duì)多種預(yù)報(bào)特色的預(yù)報(bào)產(chǎn)品基于檢驗(yàn)進(jìn)行多模式統(tǒng)計(jì)后處理集成,從而獲得效果更優(yōu)的預(yù)報(bào)結(jié)果。在多模式集成技術(shù)中,貝葉斯模式平均(Bayesian model averaging,BMA)方法是一種非常有效的提高模式預(yù)報(bào)準(zhǔn)確率的方法(Raftery et al,2005;代刊等,2018),此方法將觀測(cè)與不同模式得出的預(yù)報(bào)結(jié)果作為先驗(yàn)信息,通過(guò)求解參數(shù),計(jì)算各模式相對(duì)最優(yōu)的權(quán)重等參考值,此權(quán)重就是預(yù)報(bào)變量后驗(yàn)概率分布,代表著每個(gè)模型在訓(xùn)練階段相對(duì)的預(yù)報(bào)技巧。再經(jīng)過(guò)對(duì)偏差校正后的單個(gè)模型概率密度函數(shù)(probability density function,PDF)加權(quán)平均,得到多模式成員預(yù)報(bào)的連續(xù)的PDF,它不偏好也不摒棄各個(gè)模型,而是對(duì)各個(gè)模型結(jié)果進(jìn)行綜合,融合更多信息,以發(fā)揮各模型優(yōu)勢(shì),因此其預(yù)測(cè)均方根誤差通常小于單個(gè)預(yù)測(cè)的誤差。自從Raftery et al(2005)、Gneiting and Raftery(2005)指出BMA方法對(duì)天氣預(yù)測(cè)非常有效后, BMA方法在氣溫、風(fēng)速、降水、能見(jiàn)度等天氣和氣候預(yù)測(cè)中得到廣泛應(yīng)用。劉建國(guó)等(2013)基于TIGGE多模式集合預(yù)報(bào)資料進(jìn)行24 h氣溫的概率預(yù)報(bào),極大地提高了地面氣溫的預(yù)報(bào)技巧。智協(xié)飛等(2018)基于BMA方法對(duì)地面氣溫延伸期概率預(yù)報(bào),提高了延伸期溫度預(yù)報(bào)技巧,進(jìn)一步展示出超級(jí)集合產(chǎn)品進(jìn)行BMA概率預(yù)報(bào)的應(yīng)用潛力。陳朝平等(2010)、韓焱紅等(2013)利用BMA方法對(duì)區(qū)域和全球集合預(yù)報(bào)模式的降水概率預(yù)報(bào)產(chǎn)品進(jìn)行修正,更新的概率分布有效增強(qiáng)了暴雨可能出現(xiàn)的信號(hào)。祁海霞等(2020)基于降水集合預(yù)報(bào)資料,對(duì)清江流域的降水預(yù)報(bào)進(jìn)行試驗(yàn)分析,BMA模型預(yù)報(bào)比原始集合預(yù)報(bào)有更高預(yù)報(bào)技巧。Sloughter et al(2010;2013)的研究表明BMA訂正方法能夠明顯地改善北美洲西北部地區(qū)的最大風(fēng)速和風(fēng)向的預(yù)報(bào)。石嵐等(2017)利用BMA模型對(duì)ECMWF集合風(fēng)速預(yù)報(bào)產(chǎn)品進(jìn)行預(yù)報(bào)校準(zhǔn),減小了風(fēng)速預(yù)報(bào)誤差,MAE降低了11.7%。Chmielecki and Raftery(2011)對(duì)能見(jiàn)度要素進(jìn)行BMA概率預(yù)報(bào)研究,改善了原始集合模式能見(jiàn)度預(yù)報(bào)的概率密度函數(shù)。

    中國(guó)新一代數(shù)值天氣預(yù)報(bào)系統(tǒng)GRAPES(Global and Regional Assimilation and PrEdiction System-Mesoscales Ensemble Prediction System)是中國(guó)氣象局自主研究并初步建立的多尺度通用資料同化與數(shù)值天氣預(yù)報(bào)系統(tǒng),建成了從區(qū)域3~10 km 到全球25~50 km分辨率的確定性與集合預(yù)報(bào)的完整數(shù)值天氣預(yù)報(bào)業(yè)務(wù)技術(shù)體系(沈?qū)W順等,2020;薛紀(jì)善和陳德輝,2008),現(xiàn)統(tǒng)一更名為CMA模式。2020年CMA數(shù)值預(yù)報(bào)業(yè)務(wù)體系主要包括CMA-GFS全球中期預(yù)報(bào)系統(tǒng)(Zhang et al,2019;Chen et al,2020)、CMA-MESO區(qū)域高分辨率預(yù)報(bào)系統(tǒng)(黃麗萍等,2017)、CMA-GEPS全球集合預(yù)報(bào)系統(tǒng)(李曉莉等,2019)、CMA-REPS區(qū)域集合預(yù)報(bào)系統(tǒng)(陳靜和李曉莉,2020;張涵斌等,2014)、CMA-TYM區(qū)域臺(tái)風(fēng)預(yù)報(bào)系統(tǒng)(麻素紅和陳德輝,2018)等,為氣象預(yù)報(bào)服務(wù)和防災(zāi)減災(zāi)提供重要的支撐作用。為更好地提供氣象保障服務(wù),發(fā)揮CMA各模式的優(yōu)點(diǎn)及整個(gè)模式體系的綜合效益,需要對(duì)CMA體系產(chǎn)品進(jìn)行誤差訂正及多模式集成技術(shù)進(jìn)行模式融合,得到集成預(yù)報(bào)產(chǎn)品,為預(yù)報(bào)員提供更準(zhǔn)確的天氣預(yù)報(bào)服務(wù)。

    2022年北京冬奧會(huì)首次在冬季季風(fēng)性大陸氣候地區(qū)舉行,大部分戶(hù)外雪上賽事在地形復(fù)雜的延慶賽區(qū)和張家口賽區(qū)舉行,氣象條件在不同地點(diǎn)的差異非常大,而且賽事對(duì)氣象要素的要求也非常嚴(yán)格,因此對(duì)京津冀高精細(xì)化次公里分辨率的近地面氣象要素預(yù)報(bào),而以往對(duì)北方冬季復(fù)雜地形下近地面氣象要素的精細(xì)化預(yù)報(bào)研究較少。另一方面,之前的集成研究主要針對(duì)ECMWF集合預(yù)報(bào)產(chǎn)品以及TIGGE超級(jí)集合產(chǎn)品,對(duì)同一個(gè)模式體系不同分辨率模式進(jìn)行集成以提高模式體系的預(yù)報(bào)技巧的研究還不多見(jiàn)。本文旨在針對(duì)CMA模式體系不同分辨率模式對(duì)不同要素預(yù)報(bào)技巧存在差異并各有優(yōu)勢(shì)的特點(diǎn),對(duì)CMA不同分辨率模式進(jìn)行集成,以獲得復(fù)雜地形下更準(zhǔn)確的近地面要素預(yù)報(bào)及其概率特點(diǎn),為北京冬奧會(huì)氣象預(yù)報(bào)服務(wù)提供更精細(xì)化更有準(zhǔn)確的要素預(yù)報(bào)產(chǎn)品。

    1 多模式集成方法介紹

    本文中采用首先對(duì)CMA體系中的每一個(gè)模式分別進(jìn)行誤差訂正去除系統(tǒng)誤差,然后再進(jìn)行多模式集成的方法,得到各近地面要素的集成預(yù)報(bào)結(jié)果。

    1.1 誤差訂正方法介紹

    采用一階自適應(yīng)的卡爾曼濾波方法,這是一種基于卡爾曼濾波思想,通過(guò)不斷對(duì)模式誤差進(jìn)行更新,獲得當(dāng)前時(shí)刻的誤差估計(jì)值,來(lái)降低偏差尺度的方法(Cui et al, 2012;佟華等,2014)。這種方法既考慮了氣候平均預(yù)報(bào)誤差特征,保證了估計(jì)誤差整體的穩(wěn)定性,又加入了臨近時(shí)刻誤差信息,融入了天氣系統(tǒng)連續(xù)性特點(diǎn),將二者以權(quán)重系數(shù)相結(jié)合,共同估計(jì)遞減平均誤差。公式如下:

    Bi(t)=(1-ω)Bi(t-1)+

    ω[fi(t-1)-Oi(t-1)]

    (1)

    Fi(t)=fi(t)-Bi(t)

    (2)

    式中:t代表當(dāng)前預(yù)報(bào)時(shí)間,i代表各站點(diǎn),fi和Oi分別代表各站點(diǎn)的預(yù)報(bào)值和對(duì)應(yīng)時(shí)間的觀測(cè)值,Bi(t)定義為當(dāng)前預(yù)報(bào)時(shí)間經(jīng)過(guò)權(quán)重平均后的模式預(yù)報(bào)偏差估算值,通過(guò)權(quán)重系數(shù)ω的選擇,對(duì)兩個(gè)不同時(shí)段(前一天和歷史累積)的模式預(yù)報(bào)偏差加權(quán)平均而獲得。對(duì)于權(quán)重系數(shù),本文選用0.05,該系數(shù)是根據(jù)中國(guó)氣象局2017年參加ICE-POP 2018冬季試驗(yàn)時(shí)做的敏感性試驗(yàn)確定的(張玉濤等,2020)。利用系統(tǒng)平均加權(quán)誤差Bi(t)對(duì)預(yù)報(bào)值fi(t)進(jìn)行誤差訂正,得到各站點(diǎn)當(dāng)前預(yù)報(bào)時(shí)間訂正后的預(yù)報(bào)值Fi(t)。

    1.2 BMA多模式集成方法介紹

    BMA是一種結(jié)合多個(gè)模型進(jìn)行聯(lián)合推斷和預(yù)測(cè)的統(tǒng)計(jì)后處理方法,是多個(gè)模型后驗(yàn)分布的加權(quán)平均(Raftery et al,2005)。令f=f1,…,fK,分別假定f是K個(gè)不同數(shù)值模式的預(yù)報(bào)結(jié)果,y代表需要預(yù)報(bào)的變量,YT代表訓(xùn)練數(shù)據(jù)。BMA預(yù)報(bào)的概率密度函數(shù)PDF可表示為如下的多模式預(yù)報(bào)加權(quán)平均的形式:

    (3)

    y|(fk,YT)~gk(ak+bkfk,σ2)

    (4)

    對(duì)于風(fēng)的預(yù)報(bào),由于風(fēng)速的概率密度分布大多呈現(xiàn)偏態(tài)的特點(diǎn),很多研究使用Gamma分布作為風(fēng)速的分布函數(shù)(Sloughter et al,2010;Eide et al,2017),也有研究采用截?cái)嗾龖B(tài)分布(Frogneret al,2016;石嵐等,2017)。本研究對(duì)于研究區(qū)域大量的觀測(cè)分布顯示,風(fēng)場(chǎng)和濕度場(chǎng)可近似滿(mǎn)足截?cái)嗾龖B(tài)分布,為了方便各個(gè)要素統(tǒng)一處理,假定地面風(fēng)和相對(duì)濕度服從截?cái)嗾龖B(tài)分布,其中風(fēng)速:g(x|μ,σ)=0,x<0;相對(duì)濕度:g(x|μ,σ)=0,x<0 或x>100。

    由式(4)和式(3)可知BMA預(yù)報(bào)均值,即BMA模型多模式集成的確定性預(yù)報(bào)結(jié)果表示為:

    (5)

    然后使用訓(xùn)練數(shù)據(jù),利用EM算法(expectation maximization algorithm)進(jìn)行極大似然估計(jì),通過(guò)極大似然原則求解對(duì)數(shù)似然函數(shù)來(lái)求解BMA預(yù)報(bào)模型參數(shù)wk(k=1,2,…,K):

    L(w1,…,wk,σ2)=

    (6)

    假設(shè)預(yù)報(bào)誤差在空間s與時(shí)間t上相互獨(dú)立,上述方程由于不存在解析的極大值解,參數(shù)wk(k=1,2,…,K)與σ2需要進(jìn)行迭代求解獲得。

    BMA方法進(jìn)行多模式集成預(yù)報(bào)主要包括以下幾個(gè)步驟:(1)確定訓(xùn)練期長(zhǎng)度;(2)對(duì)訓(xùn)練期數(shù)據(jù)采用EM算法進(jìn)行模型參數(shù)率定,獲得相應(yīng)的BMA模型;(3)設(shè)計(jì)訓(xùn)練期為一個(gè)滑動(dòng)窗口進(jìn)行滾動(dòng)預(yù)報(bào),即BMA采用先前的N天作為訓(xùn)練期進(jìn)行訓(xùn)練,訓(xùn)練出的BMA系數(shù)應(yīng)用到下一天的BMA模型預(yù)報(bào)中,從而每一天動(dòng)態(tài)建立研究區(qū)域內(nèi)各點(diǎn)的BMA模型。

    2 研究區(qū)域和資料

    將試驗(yàn)區(qū)域定為東西方向約700 km、南北方向約680 km的京津冀地區(qū)范圍,分辨率為0.01°×0.01°。圖1為試驗(yàn)區(qū)域的地形高度,整個(gè)區(qū)域呈現(xiàn)西北地形高、東部和南部地形低的特點(diǎn)。西北部山區(qū)地形高度都在1000 m以上,東南部平原則基本在100 m以下。采用CMA模式體系的CMA-GFS、CMA-REPS、CMA-MESO 3 km、CMA-MESO 1 km模式系統(tǒng)進(jìn)行多模式集成試驗(yàn),試驗(yàn)要素包含2 m溫度和相對(duì)濕度、10 m風(fēng)。其中CMA-MESO 1 km 模式是面向精細(xì)化預(yù)報(bào)需求,基于CMA-MESO 3 km 高分辨率快速循環(huán)同化預(yù)報(bào)業(yè)務(wù)系統(tǒng),通過(guò)提高模式動(dòng)力框架計(jì)算精度和穩(wěn)定性,引進(jìn)更為精細(xì)的下墊面資料,選擇合適的物理過(guò)程參數(shù)化方案組合等,所建立的覆蓋京津冀的重點(diǎn)區(qū)域1 km分辨率的CMA-MESO高分辨率數(shù)值預(yù)報(bào)系統(tǒng)。首先對(duì)四個(gè)模式分別進(jìn)行空間和時(shí)間降尺度,空間降尺度采用雙線性插值,時(shí)間降尺度采用線性插值,形成京津冀地區(qū)1 km分辨率和每小時(shí)間隔的統(tǒng)一的格點(diǎn)預(yù)報(bào)場(chǎng)。再基于國(guó)家氣象信息中心1 km分辨率多源氣象數(shù)據(jù)融合格點(diǎn)實(shí)況產(chǎn)品進(jìn)行誤差訂正。這個(gè)格點(diǎn)實(shí)況產(chǎn)品是利用數(shù)據(jù)融合和數(shù)據(jù)同化技術(shù),綜合地面自動(dòng)氣象站、雷達(dá)、衛(wèi)星等多種來(lái)源觀測(cè)資料及多模式模擬數(shù)據(jù),獲得的高精度、高質(zhì)量、時(shí)空連續(xù)的多源數(shù)據(jù)融合氣象格點(diǎn)產(chǎn)品(師春香等,2018;2019)。采用一階自適應(yīng)的卡爾曼濾波方法,分別得到四個(gè)模式誤差訂正后的預(yù)報(bào)結(jié)果。其中CMA-REPS區(qū)域集合預(yù)報(bào)模式的誤差訂正采用對(duì)控制預(yù)報(bào)進(jìn)行卡爾曼濾波得到誤差文件,再應(yīng)用到每個(gè)成員進(jìn)行訂正的方法。然后使用BMA方法進(jìn)行四個(gè)模式的集成,根據(jù)選取的訓(xùn)練期,計(jì)算每個(gè)模式的權(quán)重系數(shù)w和方差,每天動(dòng)態(tài)建立BMA模型。使用各模式00時(shí)(世界時(shí),下同)起報(bào)的0~24 h預(yù)報(bào),從2020年12月1日開(kāi)始進(jìn)行訂正和集成試驗(yàn),對(duì)2021年2月1日至3月15日進(jìn)行預(yù)報(bào)并檢驗(yàn)其預(yù)報(bào)效果。

    圖1 模擬區(qū)域地形海拔高度Fig.1 Terrain height in simulated domain

    3 預(yù)報(bào)試驗(yàn)結(jié)果分析

    3.1 訓(xùn)練期長(zhǎng)度確定

    使用BMA方法進(jìn)行多模式集成,不同的訓(xùn)練期長(zhǎng)度會(huì)對(duì)集成結(jié)果產(chǎn)生影響,因此需要對(duì)BMA方法的訓(xùn)練期長(zhǎng)度確定進(jìn)行試驗(yàn)。表1分別是訓(xùn)練期為20、25、30 、35、40 d的2 m溫度、10 m風(fēng)速和2 m相對(duì)濕度24 h平均的均方根誤差比較,綜合對(duì)比各要素各時(shí)效不同訓(xùn)練期的誤差情況可見(jiàn),訓(xùn)練期為30 d時(shí)綜合誤差最小,因此使用此訓(xùn)練期長(zhǎng)度進(jìn)行后續(xù)的分析和評(píng)估。

    表1 不同訓(xùn)練期均方根誤差比較Table 1 Comparison of RMSE in different training periods

    3.2 多模式集成預(yù)報(bào)試驗(yàn)及檢驗(yàn)結(jié)果分析

    對(duì)2021年2月1日至3月15日00時(shí)起報(bào)的24 h內(nèi)的逐小時(shí)多模式集成預(yù)報(bào)進(jìn)行評(píng)估檢驗(yàn)。圖2分別為CMA-GFS、CMA-REPS、CMA-MESO 3 km、CMA-MESO 1 km四個(gè)模式2 m溫度(圖2a)、10 m風(fēng)速(圖2b)和2 m相對(duì)濕度(圖2c)在原始、訂正后和BMA集成的0~24 h預(yù)報(bào)均方根誤差格點(diǎn)平均。由圖可見(jiàn),三個(gè)要素各個(gè)預(yù)報(bào)時(shí)效誤差訂正后都較模式原始預(yù)報(bào)有明顯改進(jìn),而模式集成后的誤差在各個(gè)時(shí)效都明顯優(yōu)于每個(gè)模式的訂正結(jié)果。

    圖2 各模式不同預(yù)報(bào)時(shí)效的(a)2 m溫度,(b)10 m風(fēng)速,(c)2 m相對(duì)濕度預(yù)報(bào)的原始、訂正后和集成后均方根誤差對(duì)比Fig.2 The root-mean-square error (RMSE) comparison of the original, corrected and blended forecast about (a) 2 m temperature, (b) 10 m wind speed, (c) 2 m relative humidity of each model in 0-24 h lead time

    2 m溫度方面,四個(gè)模式的原始誤差、訂正和集成后的誤差都呈隨著預(yù)報(bào)時(shí)效的增加逐漸增大的趨勢(shì)。對(duì)于不同時(shí)效,訂正后的溫度預(yù)報(bào)誤差較訂正前有明顯減小,均方根誤差減小在0.2~1.0℃。四個(gè)模式總體效果比較來(lái)看,CMA-MESO 1 km溫度的原始預(yù)報(bào)與訂正后的預(yù)報(bào)效果優(yōu)于其他三個(gè)模式,其中原始誤差較其他模式的優(yōu)勢(shì)在0.5℃左右,訂正后優(yōu)勢(shì)為0.2℃,可見(jiàn)原始誤差越大的訂正效果越明顯。多模式集成后雖然仍像原始結(jié)果一樣存在一定日變化,總體均方根誤差在1.8~2.2℃,較訂正后預(yù)報(bào)誤差最小的模式又有0.3℃的改善,較幾個(gè)模式原始誤差的改善在0.5~1.4℃,誤差改善率為20%~40%。

    10 m風(fēng)速方面,總體看四個(gè)模式的預(yù)報(bào)效果相差明顯,尤其在9 h時(shí)效以后。區(qū)域格點(diǎn)平均的均方根誤差從小到大分別是CMA-GFS、CMA-MESO 1 km、CMA-MESO 3 km、CMA-REPS,而且誤差小的模式原始誤差甚至優(yōu)于誤差稍大的模式的訂正后誤差。各個(gè)模式中,CMA-GFS的風(fēng)速誤差最小且優(yōu)勢(shì)明顯,原始誤差在2 m·s-1左右,訂正后誤差約為1.8 m·s-1,并且不同時(shí)效的預(yù)報(bào)誤差很穩(wěn)定,并沒(méi)有其他三個(gè)模式誤差隨時(shí)效延長(zhǎng)而增長(zhǎng)的趨勢(shì)。這些都和在日常的風(fēng)速預(yù)報(bào)中全球模式如ECMWF、CMA-GFS和NCEP_GFS模式誤差穩(wěn)定小于區(qū)域模式的情況相符,分析原因認(rèn)為全球模式對(duì)于大的形勢(shì)場(chǎng)預(yù)報(bào)較區(qū)域模式更穩(wěn)定。四個(gè)模式集成后,結(jié)果仍略?xún)?yōu)于CMA-GFS訂正結(jié)果,均方根誤差在1.75 m·s-1左右,并且預(yù)報(bào)誤差穩(wěn)定,隨預(yù)報(bào)時(shí)效的增加誤差增大不明顯。BMA 集成后較幾個(gè)模式原始誤差的改善在0.3~1.4 m·s-1,誤差改善率約為12%~45%。

    2 m相對(duì)濕度的訂正和集成效果非常明顯,每個(gè)模式經(jīng)過(guò)訂正后均方根誤差都減小了4%左右,多模式集成之后誤差又減小了1%~4%。四個(gè)模式隨著預(yù)報(bào)時(shí)效的增加均方根誤差表現(xiàn)有所不同,其中CMA-REPS、CMA-MESO 3 km、CMA-MESO 1 km的原始誤差非常接近,訂正后的誤差結(jié)果也非常相似,都是在1~6 h時(shí)效誤差有所下降,之后開(kāi)始逐漸上升,訂正后的誤差從10%左右增加到24 h預(yù)報(bào)的15%。而CMA-GFS的誤差隨時(shí)效的變化明顯與其他三個(gè)模式不同,其在9 h時(shí)效之前的誤差是穩(wěn)定上升的趨勢(shì),均方根誤差比其他三個(gè)模式大3%左右,而9 h時(shí)效后誤差開(kāi)始逐漸下降,總體誤差在12%~14%,到15 h時(shí)效后均方根誤差開(kāi)始比其他三個(gè)模式更低。多模式集成預(yù)報(bào)同樣隨著時(shí)效增加誤差改善得越明顯,最終誤差在9%~12%,改善率為25%~35%。

    3.3 各模式權(quán)重對(duì)比

    圖3分別為相應(yīng)的四個(gè)模式各時(shí)效2 m溫度(圖3a)、10 m風(fēng)速(圖3b)、2 m相對(duì)濕度(圖3c)的集成權(quán)重系數(shù)對(duì)比,由圖可見(jiàn),各要素各時(shí)效的集成權(quán)重是根據(jù)模式的預(yù)報(bào)效果動(dòng)態(tài)改變的。如圖3b的10 m風(fēng)速預(yù)報(bào)中,CMA-GFS模式的風(fēng)速誤差明顯小于其他三個(gè)模式,它的權(quán)重系數(shù)達(dá)到0.4~0.55。CMA-MESO 1 km預(yù)報(bào)誤差位列次席,且隨著時(shí)效增加誤差穩(wěn)定并有減小的趨勢(shì),與CMA-GFS的誤差逐漸接近,因此在9 h時(shí)效以后權(quán)重占比明顯增加,達(dá)到0.3~0.4。CMA-MESO 3 km和CMA-REPS因誤差相對(duì)最大,權(quán)重系數(shù)很小,為 0.1 左右。對(duì)于2 m相對(duì)濕度(圖3c),各模式的集成權(quán)重系數(shù)與模式預(yù)報(bào)誤差的變化趨勢(shì)非常一致,在預(yù)報(bào)前期CMA-GFS誤差相對(duì)較大的時(shí)效,集成權(quán)重系數(shù)相對(duì)最小,隨著后期預(yù)報(bào)誤差減小,權(quán)重系數(shù)也迅速增加,到24 h預(yù)報(bào)時(shí)達(dá)到了0.4,其他三個(gè)模式則逐漸下降到0.25左右。2 m溫度(圖3a)的權(quán)重系數(shù)也基本符合上述規(guī)律,可見(jiàn)BMA多模式集成方法能夠根據(jù)訓(xùn)練期各模式的預(yù)報(bào)效果動(dòng)態(tài)調(diào)整各模式所占的權(quán)重,從而得到比所有模式預(yù)報(bào)效果更優(yōu)的確定性預(yù)報(bào)結(jié)果。

    圖3 各模式不同預(yù)報(bào)時(shí)效的(a)2 m溫度,(b)10 m風(fēng)速,(c)2 m相對(duì)濕度預(yù)報(bào)權(quán)重系數(shù)對(duì)比Fig.3 The weight coefficient comparison of each model about (a) 2 m temperature, (b) 10 m wind speed, (c) 2 m relative humidity in 0-24 h lead time

    圖4為2021年2月1日至3月15日平均CMA-MESO 1 km原始模式(圖4a,4c,4e)及多模式集成預(yù)報(bào)(圖4b,4d,4f)的京津冀地區(qū)2 m溫度(圖4a,4b)、10 m風(fēng)速(圖4c,4d)和2 m相對(duì)濕度(圖4e,4f)24 h預(yù)報(bào)均方根誤差水平分布??傮w來(lái)說(shuō)整個(gè)地區(qū)的誤差都有明顯減小,但不同要素在不同地形高度的地區(qū)誤差分布明顯不同。溫度誤差較大的地區(qū)分布在北京城區(qū)西北向地形較高的山區(qū),原始誤差可達(dá)到4℃,其他平原地區(qū)誤差在2℃左右,多模式集成后北京西部和北部山區(qū)的誤差減小到2~3℃,城區(qū)及以南地區(qū)則減小到1~2℃。10 m風(fēng)速和2 m相對(duì)濕度的誤差分布與溫度分布正相反,誤差相對(duì)較大的地區(qū)分布在城區(qū)及以南地區(qū),其原始CMA-MESO 1 km和集成后山區(qū)風(fēng)速誤差從非常不均勻的0~5 m·s-1下降到0~1 m·s-1;而在北京城區(qū)及以南的平原地區(qū),風(fēng)速誤差從2~5 m·s-1下降到1.5~2.5 m·s-1。相對(duì)濕度誤差從20%~30%降到10%左右,且在進(jìn)行集成后山區(qū)和平原地區(qū)的誤差分布非常均勻,基本不存在地形方面的差異了。綜上分析,多模式集成后三個(gè)要素誤差較大的地區(qū)預(yù)報(bào)效果都有明顯減小,訂正和集成的技術(shù)方案和效果都是較為理想的。

    圖4 2021年2月1日至3月15日平均(a,c,e)CMA-MESO 1 km原始模式與(b,d,f)多模式集成預(yù)報(bào)的京津冀地區(qū)(a,b)2 m溫度,(c,d)10 m風(fēng)速,(e,f)2 m相對(duì)濕度24 h預(yù)報(bào)均方根誤差分布Fig.4 Comparison of the 24 h forecast RMSE distribution of (a, b) 2 m temperature, (c, d) 10 m wind speed, (e, f) 2 m relative humidity between (a, c, e) the original CMA-MESO 1 km and (b, d, f) multi-model blending forcast averaged from 1 Febuary to 15 March 2021

    3.4 個(gè)例分析

    3.4.1 確定性預(yù)報(bào)分析

    對(duì)2021年3月9日00時(shí)起報(bào)24 h預(yù)報(bào)進(jìn)行了個(gè)例分析。考慮到統(tǒng)計(jì)后處理主要的優(yōu)勢(shì)在于對(duì)較穩(wěn)定存在的系統(tǒng)誤差進(jìn)行修正,因此沒(méi)有選擇較為極端的天氣個(gè)例,而選擇了較為普通的一天:3月9日為個(gè)例。圖5為3月9日00時(shí)起報(bào)24 h預(yù)報(bào)的2 m溫度(圖5a~5d)、10 m風(fēng)速(圖5e~5h)、2 m 相對(duì)濕度(圖5i~5l)的水平分布對(duì)比,其中每行從左至右,分別為觀測(cè)、CMA-MESO 1 km模式原始預(yù)報(bào)、CMA-MESO 1 km訂正結(jié)果、多模式集成結(jié)果。由圖可見(jiàn),經(jīng)過(guò)誤差訂正的模式預(yù)報(bào),都較原始預(yù)報(bào)更接近觀測(cè)實(shí)況。而再經(jīng)過(guò)多模式集成后,則較誤差訂正后的預(yù)報(bào)更進(jìn)一步接近觀測(cè)實(shí)況。如2 m 溫度預(yù)報(bào),原始預(yù)報(bào)在京津冀地區(qū)北部地形較高的山區(qū)和在北京城區(qū)及以南的平原地區(qū)較實(shí)況分析都偏高4℃左右,經(jīng)過(guò)誤差訂正后都降低了2~3℃,集成后北部山區(qū)和平原北部溫度進(jìn)一步降低,與觀測(cè)非常接近。對(duì)于10 m風(fēng)速,模式原始預(yù)報(bào)中平原地區(qū)的風(fēng)速比觀測(cè)明顯偏強(qiáng),經(jīng)過(guò)訂正和集成,這一地區(qū)風(fēng)速偏強(qiáng)的程度逐步減弱,集成之后平原地區(qū)的風(fēng)速誤差由3~4 m·s-1降低到1~2 m·s-1。對(duì)于2 m相對(duì)濕度,在大部分預(yù)報(bào)區(qū)域較觀測(cè)明顯偏低,總體偏低20%~30%,誤差訂正后誤差明顯減小,但在預(yù)報(bào)中南部地區(qū)有過(guò)訂正的問(wèn)題,相對(duì)濕度比觀測(cè)略微偏大,而集成后又重新減小,基本分布和數(shù)值都與觀測(cè)非常相似,說(shuō)明訂正集成后的結(jié)果都明顯地減小了整個(gè)區(qū)域的預(yù)報(bào)誤差。

    圖5 2021年3月9日00時(shí)起報(bào)24 h預(yù)報(bào)的(a~d)2 m溫度,(e~h)10 m風(fēng)速,(i~l)2 m相對(duì)濕度的水平分布對(duì)比(a,e,i)觀測(cè),(b,f,j)CMA-MESO 1 km模式原始預(yù)報(bào),(c,g,k)CMA-MESO 1 km訂正結(jié)果,(d,h,j)多模式集成結(jié)果Fig.5 Comparison of horizontal distribution of (a-d) 2 m temperature, (e-h) 10 m wind speed, (i-l) 2 m relative humidity in the 24 h forecast started from 00:00 UTC 9 March 2021(a, e, i) observation, (b, f, j) CMA-MESO 1 km model original forecast, (c, g, k) CMA-MESO 1 km correction, (d, h, j) multi-model blending

    3.4.2 BMA概率預(yù)報(bào)特征分析

    BMA方法不但可以得到確定性預(yù)報(bào)值,它的優(yōu)勢(shì)還在于可以得到預(yù)報(bào)場(chǎng)的概率預(yù)報(bào)特征。以2021年3月9日00時(shí)起報(bào)24 h溫度預(yù)報(bào)為例進(jìn)行BMA概率預(yù)報(bào)分析。圖6為逐小時(shí)區(qū)域格點(diǎn)平均的2 m溫度BMA百分位預(yù)報(bào)和觀測(cè)對(duì)比。BMA方法通過(guò)分析百分位預(yù)報(bào)數(shù)據(jù),可對(duì)其不確定性給出定量的估計(jì)。由圖可見(jiàn),逐小時(shí)的觀測(cè)基本落在40~60百分位上,說(shuō)明BMA預(yù)報(bào)的40~60百分位的2 m溫度能將實(shí)況比較好地包含進(jìn)去,而B(niǎo)MA確定性預(yù)報(bào)能較好的預(yù)報(bào)當(dāng)天各預(yù)報(bào)時(shí)效的溫度狀況。

    圖6 2021年3月9日00時(shí)起報(bào)逐小時(shí)區(qū)域格點(diǎn)平均的2 m溫度BMA百分位預(yù)報(bào)和觀測(cè)對(duì)比Fig.6 The area-averaging 2 m temperature comparison of BMA percentile forecasts and observations started from 00:00 UTC 9 March 2021

    對(duì)研究時(shí)段內(nèi)單個(gè)格點(diǎn)的各要素進(jìn)行BMA概率預(yù)報(bào)研究。BMA概率預(yù)報(bào)可以給出全概率的PDF結(jié)果,表征對(duì)集合預(yù)報(bào)不確定性的概率估計(jì)。如果PDF曲線越尖,區(qū)間范圍越窄,說(shuō)明BMA預(yù)報(bào)的不確定性更小。圖7分別給出地形高度較低的平原格點(diǎn)(38°N、117°E)和北部地形高度較高的山地格點(diǎn)(41.5°N、117°E)2 m溫度和10 m風(fēng)速的PDF曲線。其中圖7a和7b分別為平原格點(diǎn)和山地格點(diǎn)的2 m溫度預(yù)報(bào)PDF,圖7c和7d為兩個(gè)格點(diǎn)的10 m風(fēng)速預(yù)報(bào)PDF。圖7中黑色曲線為BMA預(yù)報(bào)的PDF,黑色垂直實(shí)線為BMA確定性預(yù)報(bào),藍(lán)色垂直實(shí)線為同一格點(diǎn)的觀測(cè),兩條紅色垂直虛線分別為BMA預(yù)報(bào)PDF的10、90百分位預(yù)報(bào)值。從2 m溫度預(yù)報(bào)的PDF曲線可見(jiàn),平原格點(diǎn)2 m 溫度的PDF曲線峰值為0.34,BMA的確定性預(yù)報(bào)位于曲線峰值的位置,與觀測(cè)之間的誤差為0.73℃,而PDF的10~90百分位溫度范圍在5.6~9.3℃,表明溫度預(yù)報(bào)最大可能落在此范圍內(nèi)。山地格點(diǎn)的2 m溫度PDF曲線峰值為0.23,確定性預(yù)報(bào)與格點(diǎn)觀測(cè)誤差為0.74℃,10~90百分位間的溫度范圍為-7.5~0.86℃,與平原格點(diǎn)比較,平原格點(diǎn)的曲線明顯更尖更窄,說(shuō)明山地格點(diǎn)溫度預(yù)報(bào)的不確定性較平原站點(diǎn)要大,與實(shí)際的預(yù)測(cè)相符。綜合上述分析,BMA方法能給出合理的不確定性的預(yù)測(cè)。

    圖7 2021年3月9日00時(shí)起報(bào)(a,b)2 m溫度和(c,d)10 m風(fēng)速的24 h預(yù)報(bào)PDF(a,c)平原格點(diǎn)(38°N、117°E),(b,d)山地格點(diǎn)(41.5°N、117°E)(黑色曲線:BMA預(yù)報(bào)PDF,黑色垂直實(shí)線:BMA確定性預(yù)報(bào),藍(lán)色垂直實(shí)線:觀測(cè),紅色垂直虛線BMA預(yù)報(bào)PDF的10和90百分位預(yù)報(bào))Fig.7 BMA predictive PDF of (a, b) 2 m temperature and (c, d) 10 m wind speed in 24 h forecast started from 00:00 UTC 9 March 2021 at (a, c) Plain grid (38°N, 117°E) and (b, d) Hilly Area grid (41.5°N, 117°E) (black solid curve: BMA predictive PDF, black solid vertical line: BMA deterministic forecast, blue solid vertical line: observation, red dashed vertical lines: the 10th percentile and 90th percentile forecasts for BMA predictive PDF)

    從10 m風(fēng)速預(yù)報(bào)的PDF看,兩個(gè)格點(diǎn)的PDF曲線也有明顯差異。山區(qū)格點(diǎn)較平原格點(diǎn)相比, PDF曲線峰值較高,曲線又尖又窄,表明此平原站點(diǎn)預(yù)報(bào)的不確定性更大。這個(gè)特點(diǎn)與CMA模式在平原的風(fēng)速預(yù)報(bào)偏大明顯,預(yù)報(bào)誤差較大有關(guān),在應(yīng)用此預(yù)報(bào)結(jié)果時(shí)需要采用一定的概率值。

    從單個(gè)格點(diǎn)的概率預(yù)報(bào)可了解溫度和風(fēng)速在某一格點(diǎn)上預(yù)報(bào)值的不確定性,但是單個(gè)格點(diǎn)不能完全反映區(qū)域整體特點(diǎn)。接著從區(qū)域的角度分析CMA模式BMA概率預(yù)報(bào)的特點(diǎn)。圖8給出了CMA多模式集成的2 m溫度以BMA 確定性結(jié)果為中心,區(qū)間為2℃的概率分布(圖 8a),以及10 m風(fēng)速以BMA確定性結(jié)果為中心,區(qū)間為2 m·s-1的概率分布(圖8b)。由此判斷不同地區(qū)選取哪一個(gè)概率下的預(yù)報(bào)結(jié)果更合理。從2 m溫度預(yù)報(bào)的概率分布可見(jiàn),概率分布主要受地形高度和海陸分布的影響,在右側(cè)渤海灣地區(qū)的概率值最大,陸地較海洋的概率值是減小的。陸地上在北京城區(qū)及以南的河北省的平原地區(qū),概率值較大,為50%~80%,而在區(qū)域西北部山區(qū)概率值是減小的,概率在20%~50%。說(shuō)明山區(qū)的溫度預(yù)報(bào)方差和預(yù)報(bào)不確定性是較大的。這與BMA的2 m溫度預(yù)報(bào)的誤差分布很相似,預(yù)報(bào)不確定性較大的地區(qū)其BMA確定性預(yù)報(bào)的誤差較大。

    圖8 CMA模式(a)2 m溫度以BMA確定性結(jié)果為中心,區(qū)間為2℃的概率分布;(b)10 m風(fēng)速以BMA確定性結(jié)果為中心,區(qū)間為2 m·s-1的概率分布Fig.8 The probability distribution of (a) 2 m temperature centered as BMA deterministic hindcast with the interval length of 2℃, and (b) 10 m wind speed centered as BMA deterministic hindcast with the interval length of 2 m·s-1

    對(duì)于10 m風(fēng)速的概率分布,概率值較大的地區(qū)位于試驗(yàn)區(qū)域的西北部山區(qū)以及與平原地區(qū)的交界處,概率可達(dá)70%~90%,說(shuō)明預(yù)報(bào)值的不確定性很小,風(fēng)速的預(yù)報(bào)是較為準(zhǔn)確的。而在其他地區(qū)概率有所減小,在30%~50%。這一分布與原始CMA模式風(fēng)場(chǎng)預(yù)報(bào)誤差相似。概率分布的情況可以很好地說(shuō)明預(yù)報(bào)的不確定性,因此可根據(jù)區(qū)域差異合理選取某一概率下的預(yù)報(bào)結(jié)果。

    4 結(jié)論與討論

    基于CMA模式體系的CMA-GFS、CMA-REPS、CMA-MESO 3 km、CMA-MESO 1 km四個(gè)模式的2 m溫度、10 m風(fēng)速、2 m相對(duì)濕度等近地面要素預(yù)報(bào),對(duì)京津冀地區(qū)進(jìn)行誤差訂正和貝葉斯模式平均(BMA)方法多模式集成。對(duì)比了每個(gè)模式原始誤差、訂正后的誤差以及多模式集成誤差;比較了不同訓(xùn)練期結(jié)果確定最優(yōu)訓(xùn)練期;分析了模式集成過(guò)程中各模式的權(quán)重特點(diǎn),以及誤差在預(yù)報(bào)區(qū)域內(nèi)的水平分布特點(diǎn)等,得到以下主要結(jié)論:

    (1)通過(guò)最優(yōu)訓(xùn)練期長(zhǎng)度選取試驗(yàn),發(fā)現(xiàn)在20~40 d,不同訓(xùn)練期對(duì)BMA方法的要素預(yù)報(bào)的誤差先減小后增大,訓(xùn)練期長(zhǎng)度在30 d可得到綜合最優(yōu)的預(yù)報(bào)效果。

    (2)通過(guò)對(duì)比每個(gè)模式原始誤差、訂正后的誤差以及多模式集成誤差,結(jié)果顯示每個(gè)模式訂正后的均方根誤差都較原始模式預(yù)報(bào)有明顯的減小,而之后多模式集成的預(yù)報(bào)效果比其中任一模式訂正后的誤差更優(yōu)。2 m溫度集成預(yù)報(bào)誤差較原始模式的改善在0.5~1.4℃,改善率為20%~40%;10 m風(fēng)速和2 m相對(duì)濕度的均方根誤差改善率分別為12%~45%和25%~35%。

    (3)通過(guò)BMA方法四個(gè)模式各時(shí)效的集成權(quán)重系數(shù)對(duì)比,顯示各要素各時(shí)效的集成權(quán)重根據(jù)各模式的預(yù)報(bào)效果動(dòng)態(tài)改變,當(dāng)某一模式預(yù)報(bào)誤差相對(duì)較小時(shí),集成權(quán)重系數(shù)相對(duì)較大,反之亦然。如CMA-GFS模式的風(fēng)速誤差明顯小于其他三個(gè)模式,其權(quán)重系數(shù)達(dá)到0.4~0.55。

    (4)通過(guò)對(duì)比CMA-MESO 1 km原始模式及訂正、集成后的各要素預(yù)報(bào),從均方根誤差的多日平均水平分布和個(gè)例分析,不同要素在不同地形高度處誤差分布明顯不同,溫度誤差較大的地方分布在北京城區(qū)西北地形較高的山區(qū),10 m風(fēng)速和2 m相對(duì)濕度的誤差分布與溫度分布正相反,相對(duì)較大的地區(qū)分布在城區(qū)及以南地區(qū)。經(jīng)過(guò)訂正集成后整個(gè)地區(qū)的誤差都有明顯減小,要素水平分布通過(guò)訂正和集成逐漸向觀測(cè)靠近,溫度誤差減小1~2℃,風(fēng)速誤差減小0.5~4 m·s-1,相對(duì)濕度降低10%~20%。

    (5)BMA方法通過(guò)分析其百分位預(yù)報(bào)數(shù)據(jù),對(duì)不確定性給出定量的預(yù)計(jì)。BMA預(yù)報(bào)的概率分布情況能較好地說(shuō)明預(yù)報(bào)的不確定性,并將實(shí)際大氣可能發(fā)生狀態(tài)縮小到一個(gè)更小的區(qū)間范圍,預(yù)報(bào)的不確定性減小。

    通過(guò)將每個(gè)模式分別進(jìn)行誤差訂正后再進(jìn)行多模式集成的方法,既能夠明顯減小單個(gè)模式的預(yù)報(bào)誤差,又能對(duì)各個(gè)模型結(jié)果進(jìn)行綜合集成,發(fā)揮各模型優(yōu)勢(shì),使模式集成預(yù)報(bào)效果優(yōu)于任一單個(gè)模式。因此這一方法切實(shí)有效,既匯集了CMA模式體系總體優(yōu)勢(shì),獲得可靠的確定性預(yù)報(bào)結(jié)果,也提供完整的概率密度函數(shù)(PDF),對(duì)極端天氣事件的概率預(yù)報(bào)技巧進(jìn)行研究和評(píng)估。

    猜你喜歡
    格點(diǎn)方根時(shí)效
    方根拓展探究
    帶有超二次位勢(shì)無(wú)限格點(diǎn)上的基態(tài)行波解
    一種電離層TEC格點(diǎn)預(yù)測(cè)模型
    帶可加噪聲的非自治隨機(jī)Boussinesq格點(diǎn)方程的隨機(jī)吸引子
    均方根嵌入式容積粒子PHD 多目標(biāo)跟蹤方法
    J75鋼的時(shí)效處理工藝
    一種新型耐熱合金GY200的長(zhǎng)期時(shí)效組織與性能
    上海金屬(2016年3期)2016-11-23 05:19:47
    環(huán)保執(zhí)法如何把握對(duì)違法建設(shè)項(xiàng)目的追責(zé)時(shí)效?
    揭開(kāi)心算方根之謎
    格點(diǎn)和面積
    热99国产精品久久久久久7| 一边摸一边抽搐一进一出视频| 他把我摸到了高潮在线观看 | 国产精品九九99| 日韩欧美国产一区二区入口| 在线 av 中文字幕| 美国免费a级毛片| 成年女人毛片免费观看观看9 | a 毛片基地| 亚洲国产欧美网| 美女福利国产在线| 一级毛片精品| 久久久精品国产亚洲av高清涩受| 91av网站免费观看| 久久久久视频综合| 欧美精品啪啪一区二区三区 | 久久国产精品影院| 国产真人三级小视频在线观看| 久久国产精品大桥未久av| 久9热在线精品视频| 中国国产av一级| 亚洲一区中文字幕在线| 午夜福利乱码中文字幕| 亚洲欧美一区二区三区黑人| 高清黄色对白视频在线免费看| 热99re8久久精品国产| 黄色片一级片一级黄色片| 国产日韩一区二区三区精品不卡| 久久久精品区二区三区| 亚洲久久久国产精品| 日韩大码丰满熟妇| 久久影院123| 精品福利观看| 亚洲视频免费观看视频| 丝袜人妻中文字幕| 人人澡人人妻人| 亚洲熟女毛片儿| av视频免费观看在线观看| 99国产精品一区二区蜜桃av | 在线观看免费午夜福利视频| 精品国内亚洲2022精品成人 | 看免费av毛片| 在线亚洲精品国产二区图片欧美| 午夜福利免费观看在线| 久久这里只有精品19| 亚洲 欧美一区二区三区| 亚洲中文字幕日韩| 欧美97在线视频| 性色av乱码一区二区三区2| 少妇人妻久久综合中文| 日韩免费高清中文字幕av| 精品人妻1区二区| 国产亚洲欧美在线一区二区| 99久久国产精品久久久| 亚洲精品国产av成人精品| 亚洲伊人色综图| 成人免费观看视频高清| 亚洲国产av影院在线观看| 一级毛片精品| 乱人伦中国视频| 爱豆传媒免费全集在线观看| 国产老妇伦熟女老妇高清| 韩国精品一区二区三区| 丁香六月天网| 国产xxxxx性猛交| 午夜福利一区二区在线看| 欧美黄色片欧美黄色片| av视频免费观看在线观看| 精品人妻1区二区| 精品福利观看| 亚洲欧美激情在线| 中文字幕高清在线视频| 不卡av一区二区三区| 岛国在线观看网站| 一区二区三区乱码不卡18| 在线十欧美十亚洲十日本专区| 国产精品久久久久成人av| 国产精品熟女久久久久浪| 精品少妇一区二区三区视频日本电影| 热re99久久精品国产66热6| 午夜影院在线不卡| 中国美女看黄片| 中国国产av一级| 啦啦啦啦在线视频资源| 人人妻人人澡人人爽人人夜夜| 国产精品成人在线| 国产精品免费大片| 精品一区二区三区av网在线观看 | 亚洲av美国av| 成在线人永久免费视频| 日本a在线网址| 久久久久国内视频| 亚洲天堂av无毛| 成年人黄色毛片网站| 操出白浆在线播放| 国产成人一区二区三区免费视频网站| 欧美日韩亚洲国产一区二区在线观看 | 精品国内亚洲2022精品成人 | 亚洲男人天堂网一区| 91精品国产国语对白视频| 亚洲中文日韩欧美视频| 在线十欧美十亚洲十日本专区| 国产精品免费视频内射| a 毛片基地| 婷婷色av中文字幕| 成人国产一区最新在线观看| 男人爽女人下面视频在线观看| 99热全是精品| 男人舔女人的私密视频| 丝瓜视频免费看黄片| 在线十欧美十亚洲十日本专区| 精品少妇内射三级| 91精品伊人久久大香线蕉| 欧美久久黑人一区二区| 久久影院123| 高清av免费在线| 成人三级做爰电影| 人人妻人人澡人人看| 亚洲人成电影免费在线| 菩萨蛮人人尽说江南好唐韦庄| 亚洲一卡2卡3卡4卡5卡精品中文| 成人手机av| 免费人妻精品一区二区三区视频| av网站免费在线观看视频| 中亚洲国语对白在线视频| 永久免费av网站大全| 亚洲黑人精品在线| 99精品久久久久人妻精品| 丰满迷人的少妇在线观看| 亚洲精品第二区| 久久毛片免费看一区二区三区| 男女下面插进去视频免费观看| 夜夜夜夜夜久久久久| 99国产精品一区二区三区| 亚洲国产毛片av蜜桃av| 最近最新中文字幕大全免费视频| www.精华液| 亚洲午夜精品一区,二区,三区| 亚洲第一av免费看| 满18在线观看网站| 久久久国产一区二区| 国产成人精品久久二区二区免费| av在线老鸭窝| 国产日韩欧美在线精品| 亚洲欧美一区二区三区黑人| 一级毛片电影观看| 色老头精品视频在线观看| 美女中出高潮动态图| 成年美女黄网站色视频大全免费| 欧美国产精品一级二级三级| 黑人巨大精品欧美一区二区蜜桃| 色94色欧美一区二区| 女警被强在线播放| 日本av手机在线免费观看| 婷婷色av中文字幕| 免费高清在线观看日韩| 自拍欧美九色日韩亚洲蝌蚪91| 亚洲专区字幕在线| 国产真人三级小视频在线观看| 亚洲精品久久久久久婷婷小说| 两个人看的免费小视频| 飞空精品影院首页| 精品久久久久久久毛片微露脸 | 精品国产乱码久久久久久小说| 国产精品久久久久久人妻精品电影 | 久久99热这里只频精品6学生| av有码第一页| 久久亚洲国产成人精品v| 国产精品久久久久久精品古装| 精品福利观看| 男人爽女人下面视频在线观看| 国产视频一区二区在线看| 精品人妻1区二区| 蜜桃在线观看..| 久久精品人人爽人人爽视色| 纵有疾风起免费观看全集完整版| 视频在线观看一区二区三区| 别揉我奶头~嗯~啊~动态视频 | 狂野欧美激情性bbbbbb| 久久久精品94久久精品| 日韩大片免费观看网站| 欧美日韩视频精品一区| 久久精品成人免费网站| 久久性视频一级片| 99热网站在线观看| 三级毛片av免费| 日本vs欧美在线观看视频| 国产一区二区在线观看av| 99精品久久久久人妻精品| 9191精品国产免费久久| 国产在线观看jvid| 69av精品久久久久久 | 色播在线永久视频| 亚洲国产看品久久| 一进一出抽搐动态| 日韩免费高清中文字幕av| 精品亚洲成国产av| 亚洲人成电影观看| 又黄又粗又硬又大视频| 精品国产一区二区三区四区第35| 一进一出抽搐动态| 水蜜桃什么品种好| 欧美激情极品国产一区二区三区| 久久久久精品国产欧美久久久 | 香蕉丝袜av| 国产极品粉嫩免费观看在线| 美女扒开内裤让男人捅视频| 大片电影免费在线观看免费| 精品久久久久久久毛片微露脸 | 成人国语在线视频| 99久久99久久久精品蜜桃| 12—13女人毛片做爰片一| 亚洲中文av在线| 欧美日韩成人在线一区二区| 久久久久视频综合| 亚洲欧美精品综合一区二区三区| 精品亚洲成国产av| 亚洲视频免费观看视频| 精品免费久久久久久久清纯 | 亚洲国产精品成人久久小说| 日本欧美视频一区| 久久人妻福利社区极品人妻图片| 日韩有码中文字幕| 精品国产乱子伦一区二区三区 | a在线观看视频网站| 久久久精品94久久精品| 亚洲精品一卡2卡三卡4卡5卡 | 黄频高清免费视频| 国产精品秋霞免费鲁丝片| 国精品久久久久久国模美| 美女视频免费永久观看网站| 久久av网站| 久久精品国产亚洲av高清一级| 在线观看免费日韩欧美大片| 正在播放国产对白刺激| 夜夜骑夜夜射夜夜干| 亚洲专区中文字幕在线| 亚洲欧美色中文字幕在线| 精品一区二区三卡| 久久天堂一区二区三区四区| 欧美在线一区亚洲| 大型av网站在线播放| 热99re8久久精品国产| 国产男人的电影天堂91| 成人18禁高潮啪啪吃奶动态图| 三级毛片av免费| 欧美黄色淫秽网站| e午夜精品久久久久久久| 欧美一级毛片孕妇| 国产精品久久久久久人妻精品电影 | 正在播放国产对白刺激| 国产精品 欧美亚洲| 五月天丁香电影| 日本欧美视频一区| 免费看十八禁软件| 成年人午夜在线观看视频| 人妻一区二区av| 亚洲伊人久久精品综合| 97人妻天天添夜夜摸| 亚洲国产日韩一区二区| 国产区一区二久久| 蜜桃国产av成人99| 肉色欧美久久久久久久蜜桃| 久久人人爽av亚洲精品天堂| 欧美精品高潮呻吟av久久| 久久久久国内视频| 一边摸一边抽搐一进一出视频| 99国产精品免费福利视频| 首页视频小说图片口味搜索| 亚洲熟女毛片儿| 岛国毛片在线播放| 欧美国产精品一级二级三级| 日韩电影二区| 欧美人与性动交α欧美精品济南到| 国产在视频线精品| 亚洲久久久国产精品| 一本色道久久久久久精品综合| 五月天丁香电影| 亚洲第一青青草原| 女人精品久久久久毛片| 超碰成人久久| 美女中出高潮动态图| 久久人人爽人人片av| 午夜福利一区二区在线看| 日韩欧美一区二区三区在线观看 | 秋霞在线观看毛片| 91字幕亚洲| 亚洲第一青青草原| 国产一区二区三区在线臀色熟女 | 久久久水蜜桃国产精品网| 免费观看av网站的网址| 中国美女看黄片| 最新在线观看一区二区三区| 自拍欧美九色日韩亚洲蝌蚪91| 亚洲国产精品一区三区| 超碰成人久久| 国产亚洲av片在线观看秒播厂| 国产欧美亚洲国产| 欧美中文综合在线视频| 韩国精品一区二区三区| 久久久久国内视频| 国产精品熟女久久久久浪| 亚洲精品美女久久久久99蜜臀| 精品国产一区二区久久| 在线观看www视频免费| av网站免费在线观看视频| 黄色毛片三级朝国网站| 一区二区日韩欧美中文字幕| 国产区一区二久久| 亚洲精品国产区一区二| 69av精品久久久久久 | av在线播放精品| 国产免费av片在线观看野外av| 国产国语露脸激情在线看| 欧美精品亚洲一区二区| 90打野战视频偷拍视频| 正在播放国产对白刺激| 日本欧美视频一区| 成人黄色视频免费在线看| 一级黄色大片毛片| 久久女婷五月综合色啪小说| 国产成人av激情在线播放| 国产一区有黄有色的免费视频| 亚洲美女黄色视频免费看| 欧美性长视频在线观看| 精品久久久久久久毛片微露脸 | 自线自在国产av| 老鸭窝网址在线观看| 亚洲国产欧美网| 黄片小视频在线播放| 亚洲精品乱久久久久久| 欧美+亚洲+日韩+国产| 波多野结衣一区麻豆| av免费在线观看网站| 午夜福利免费观看在线| 亚洲情色 制服丝袜| 巨乳人妻的诱惑在线观看| 别揉我奶头~嗯~啊~动态视频 | 国产色视频综合| 午夜福利影视在线免费观看| 欧美黑人欧美精品刺激| 美女扒开内裤让男人捅视频| 久久久国产欧美日韩av| av福利片在线| 亚洲精品成人av观看孕妇| 亚洲七黄色美女视频| 91精品国产国语对白视频| 国产欧美日韩一区二区三 | 丰满少妇做爰视频| 成年人免费黄色播放视频| √禁漫天堂资源中文www| 91字幕亚洲| 99热国产这里只有精品6| 日韩 亚洲 欧美在线| 久久人妻福利社区极品人妻图片| 啦啦啦在线免费观看视频4| 久久毛片免费看一区二区三区| 欧美黄色淫秽网站| 成人国产一区最新在线观看| 日本撒尿小便嘘嘘汇集6| 妹子高潮喷水视频| 亚洲国产欧美一区二区综合| 免费日韩欧美在线观看| 美女高潮喷水抽搐中文字幕| 亚洲精品国产色婷婷电影| 无遮挡黄片免费观看| 69精品国产乱码久久久| 动漫黄色视频在线观看| 视频在线观看一区二区三区| 我要看黄色一级片免费的| 免费高清在线观看日韩| 精品少妇久久久久久888优播| 美女中出高潮动态图| 少妇人妻久久综合中文| 嫩草影视91久久| 亚洲欧美成人综合另类久久久| 激情视频va一区二区三区| av网站免费在线观看视频| 手机成人av网站| 99国产综合亚洲精品| 丰满少妇做爰视频| 国产男女超爽视频在线观看| 啦啦啦中文免费视频观看日本| 日韩欧美一区二区三区在线观看 | 成年人午夜在线观看视频| 99精国产麻豆久久婷婷| 大片免费播放器 马上看| 成人亚洲精品一区在线观看| 免费少妇av软件| 午夜激情久久久久久久| 精品一区二区三卡| av一本久久久久| 久久久国产欧美日韩av| 国产成人一区二区三区免费视频网站| 每晚都被弄得嗷嗷叫到高潮| 成年人黄色毛片网站| 亚洲国产欧美日韩在线播放| 波多野结衣一区麻豆| 亚洲三区欧美一区| 欧美成狂野欧美在线观看| 亚洲欧美日韩另类电影网站| 亚洲av日韩精品久久久久久密| 日日摸夜夜添夜夜添小说| 美女高潮喷水抽搐中文字幕| 国产成人av激情在线播放| 少妇被粗大的猛进出69影院| 美女扒开内裤让男人捅视频| 美女午夜性视频免费| 18禁裸乳无遮挡动漫免费视频| 久久人人97超碰香蕉20202| 欧美日韩视频精品一区| 精品国产超薄肉色丝袜足j| 1024视频免费在线观看| 欧美精品亚洲一区二区| 亚洲精品久久久久久婷婷小说| 午夜激情av网站| 俄罗斯特黄特色一大片| 欧美激情久久久久久爽电影 | cao死你这个sao货| 国产精品熟女久久久久浪| 啪啪无遮挡十八禁网站| 午夜激情av网站| 日韩大码丰满熟妇| 亚洲熟女毛片儿| 男人操女人黄网站| 黄色 视频免费看| 欧美成狂野欧美在线观看| 中文字幕制服av| 91国产中文字幕| 亚洲欧美激情在线| 丝袜人妻中文字幕| 国产日韩欧美视频二区| 国产熟女午夜一区二区三区| 国产高清视频在线播放一区 | 欧美另类一区| 99re6热这里在线精品视频| 日日爽夜夜爽网站| 黄色视频不卡| 人人澡人人妻人| 女性生殖器流出的白浆| 黄色毛片三级朝国网站| 国产精品九九99| 成年动漫av网址| 久久久久视频综合| 99国产极品粉嫩在线观看| 18禁黄网站禁片午夜丰满| 91成人精品电影| 亚洲综合色网址| 午夜福利一区二区在线看| 色视频在线一区二区三区| 欧美变态另类bdsm刘玥| 色视频在线一区二区三区| 免费在线观看完整版高清| 久久 成人 亚洲| 女性生殖器流出的白浆| 美女福利国产在线| 最近最新中文字幕大全免费视频| 老熟妇乱子伦视频在线观看 | 国产av国产精品国产| 日韩,欧美,国产一区二区三区| 久久影院123| av有码第一页| 91精品伊人久久大香线蕉| 国产三级黄色录像| 两性午夜刺激爽爽歪歪视频在线观看 | a级毛片在线看网站| 建设人人有责人人尽责人人享有的| 汤姆久久久久久久影院中文字幕| 秋霞在线观看毛片| 国产成人av教育| 99久久人妻综合| 成人三级做爰电影| 久久精品aⅴ一区二区三区四区| 国产成人欧美在线观看 | 91精品三级在线观看| 我的亚洲天堂| 国产在线免费精品| 女性被躁到高潮视频| 黄片大片在线免费观看| 激情视频va一区二区三区| 黄片大片在线免费观看| 乱人伦中国视频| 日本vs欧美在线观看视频| 一区福利在线观看| 国产精品 欧美亚洲| 成年人免费黄色播放视频| 又黄又粗又硬又大视频| 人成视频在线观看免费观看| 精品国产超薄肉色丝袜足j| 久久九九热精品免费| 97在线人人人人妻| 国产成人免费观看mmmm| 亚洲全国av大片| 美女主播在线视频| 丝袜在线中文字幕| 国产国语露脸激情在线看| 天天躁夜夜躁狠狠躁躁| 如日韩欧美国产精品一区二区三区| 午夜精品国产一区二区电影| 亚洲人成电影免费在线| 亚洲国产看品久久| 黄色怎么调成土黄色| 国产三级黄色录像| 丝袜脚勾引网站| 色精品久久人妻99蜜桃| 91字幕亚洲| 欧美另类亚洲清纯唯美| 久久中文看片网| 国产免费一区二区三区四区乱码| 最近中文字幕2019免费版| av免费在线观看网站| av福利片在线| 欧美av亚洲av综合av国产av| 久久影院123| 欧美日本中文国产一区发布| 王馨瑶露胸无遮挡在线观看| 五月开心婷婷网| 99久久综合免费| 美女大奶头黄色视频| 飞空精品影院首页| 国产亚洲欧美精品永久| 一二三四社区在线视频社区8| 亚洲精品国产av成人精品| 超色免费av| 新久久久久国产一级毛片| 国产人伦9x9x在线观看| 老司机影院毛片| 侵犯人妻中文字幕一二三四区| 一区二区av电影网| 亚洲国产看品久久| 黄色毛片三级朝国网站| av超薄肉色丝袜交足视频| 亚洲欧美精品综合一区二区三区| 一区二区av电影网| 久久天堂一区二区三区四区| 一级毛片女人18水好多| cao死你这个sao货| 亚洲全国av大片| 一本色道久久久久久精品综合| 亚洲国产精品999| 精品一区二区三区av网在线观看 | 久久精品aⅴ一区二区三区四区| 亚洲精品av麻豆狂野| 一个人免费在线观看的高清视频 | 亚洲中文av在线| 9热在线视频观看99| 久久久久久久久久久久大奶| 一级a爱视频在线免费观看| 一进一出抽搐动态| 美女午夜性视频免费| 日日夜夜操网爽| 人人妻人人澡人人爽人人夜夜| 一二三四在线观看免费中文在| 国产麻豆69| 在线观看一区二区三区激情| 精品国产乱子伦一区二区三区 | 又紧又爽又黄一区二区| 国产欧美日韩综合在线一区二区| avwww免费| 日日爽夜夜爽网站| 午夜福利在线观看吧| 亚洲少妇的诱惑av| 人人澡人人妻人| 亚洲av日韩在线播放| 久久久久久人人人人人| 国产精品欧美亚洲77777| 国产男女超爽视频在线观看| 国产亚洲一区二区精品| 爱豆传媒免费全集在线观看| 久久久水蜜桃国产精品网| 黄色视频在线播放观看不卡| 成人手机av| 亚洲国产看品久久| 精品少妇内射三级| 免费观看av网站的网址| 久久久久久人人人人人| 高清在线国产一区| 亚洲精品美女久久久久99蜜臀| 男人爽女人下面视频在线观看| 视频区图区小说| 天天躁狠狠躁夜夜躁狠狠躁| av在线播放精品| 久久天躁狠狠躁夜夜2o2o| 欧美日韩亚洲高清精品| av一本久久久久| 少妇裸体淫交视频免费看高清 | 99久久综合免费| cao死你这个sao货| 成人av一区二区三区在线看 | 少妇精品久久久久久久| 我的亚洲天堂| 久久久久久久久久久久大奶| 新久久久久国产一级毛片| 午夜免费鲁丝| av超薄肉色丝袜交足视频| 国产野战对白在线观看| 亚洲精品久久午夜乱码| a级毛片在线看网站| 免费黄频网站在线观看国产| 窝窝影院91人妻| 最新的欧美精品一区二区| 日本欧美视频一区| 热re99久久国产66热| 51午夜福利影视在线观看| 久久人人97超碰香蕉20202| 久久午夜综合久久蜜桃| 国产av精品麻豆| 悠悠久久av| 欧美少妇被猛烈插入视频| 久久精品人人爽人人爽视色| 性色av一级| 人妻久久中文字幕网| 国产欧美亚洲国产| 一边摸一边做爽爽视频免费| 91麻豆av在线|