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

    集合均方根濾波同化地面自動(dòng)站資料的技術(shù)研究

    2015-12-05 07:49:12邵長(zhǎng)亮閔錦忠
    大氣科學(xué) 2015年1期

    邵長(zhǎng)亮 閔錦忠

    ?

    集合均方根濾波同化地面自動(dòng)站資料的技術(shù)研究

    邵長(zhǎng)亮1, 2, 3閔錦忠1, 3

    1南京信息工程大學(xué)氣象災(zāi)害預(yù)報(bào)預(yù)警與評(píng)估協(xié)同創(chuàng)新中心,南京210044;2中國(guó)氣象局氣象探測(cè)中心,北京100081;3氣象災(zāi)害教育部重點(diǎn)實(shí)驗(yàn)室,南京210044

    模式地形與觀測(cè)站地形高度差異一直是地面資料同化面臨的棘手問(wèn)題,合理的同化方案能夠?qū)⒌孛孀詣?dòng)站資料有效的同化到中尺度數(shù)值模式中。本文首先采用Guo et al.(2002)的方案實(shí)現(xiàn)了在WRF模式中應(yīng)用集合Kalman濾波方法同化地面自動(dòng)站資料;然后對(duì)方案進(jìn)行調(diào)整,對(duì)10 m高度風(fēng)場(chǎng)、2 m高度位溫、2 m高度露點(diǎn)和地表氣壓進(jìn)行同化。通過(guò)均方根誤差分析,模擬結(jié)果和同化增量分析來(lái)確定集合平方根濾波(EnSRF)同化地面自動(dòng)站資料的有效性,并進(jìn)行敏感性試驗(yàn)分析檢驗(yàn)?zāi)J綄?duì)各要素物理量的響應(yīng)狀況。結(jié)果表明:在EnSRF同化系統(tǒng)中應(yīng)用Guo et al.(2002)的方案將地面自動(dòng)站資料進(jìn)行同化到數(shù)值模式中,能夠部分改善模擬結(jié)果;地面觀測(cè)資料(溫度、濕度、風(fēng)場(chǎng)、地表氣壓)中各物理量分別同化到數(shù)值模式都能影響18小時(shí)降水預(yù)報(bào),但各物理量所起作用大小不同,其中對(duì)結(jié)果影響最大的是露點(diǎn);使用位溫、露點(diǎn)分別代替溫度、比濕進(jìn)行同化模擬效果更好,對(duì)自動(dòng)站資料的同化也更加有效。

    資料同化 集合Kalman濾波 自動(dòng)站資料 敏感性試驗(yàn)

    1 引言

    中尺度對(duì)流天氣如暴雨、冰雹、龍卷等往往會(huì)造成比較嚴(yán)重的自然災(zāi)害, 給人們的生命財(cái)產(chǎn)帶來(lái)重大損失。而我國(guó)又是一個(gè)暴雨頻發(fā)的國(guó)家, 因此對(duì)中尺度天氣系統(tǒng)的預(yù)報(bào)就顯得尤為重要。當(dāng)前主要通過(guò)數(shù)值預(yù)報(bào)模式來(lái)實(shí)現(xiàn)對(duì)中尺度天氣的預(yù)報(bào),而要得到準(zhǔn)確的數(shù)值天氣預(yù)報(bào),需要具備兩個(gè)條件:一個(gè)是準(zhǔn)確的初值, 另一個(gè)是準(zhǔn)確的大氣運(yùn)動(dòng)規(guī)律。隨著模式的不斷發(fā)展完善,對(duì)初始條件的精確性要求也日趨提高,因而資料同化的作用也愈發(fā)凸顯。

    資料同化發(fā)展至今經(jīng)歷了從屬于三維同化的多項(xiàng)式擬合、逐步訂正、最優(yōu)插值(OI)、三維變分(3DVAR)同化到屬于四維同化的四維變分(4DVAR)和集合卡爾曼濾波(EnKF)幾個(gè)不同的階段。3DVAR因其能夠直接同化非常規(guī)資料的優(yōu)勢(shì),已經(jīng)取代OI成為當(dāng)前主流的同化方法,在國(guó)內(nèi)外主要的中尺度數(shù)值預(yù)報(bào)模式平臺(tái)上都已構(gòu)建了相應(yīng)的3DVAR同化系統(tǒng)(Gao et al., 2001;Barker et al., 2004;薛紀(jì)善等,2008)。盡管取得了較好的效果,但因?yàn)?DVAR同化方法本身所具有的局限性,所以為了給數(shù)值模式提供更好的初始場(chǎng),人們已將當(dāng)前研究的重心轉(zhuǎn)向了四維同化,其中4DVAR和EnKF是當(dāng)前主流的四維同化方法。在主要的中尺度模式上, 包括雷達(dá)和衛(wèi)星等非常規(guī)資料在內(nèi)的4DVAR同化得到了一定的研究(De Pondeca and Zou,2001;Martinelli et al., 2003;Huang et al., 2009),此外,為了回避4DVAR需要復(fù)雜伴隨模式的問(wèn)題,一種顯式的4DVAR(Qiu and Chou,2006)方法在近幾年也得到了一定的研究。EnKF作為與4DVAR并列的未來(lái)可能替換已有業(yè)務(wù)3DVAR同化系統(tǒng)的一種選擇, 從1994年Evensen(1994)提出至今, 已經(jīng)得到了廣泛的研究,對(duì)包括EnKF在風(fēng)暴尺度領(lǐng)域應(yīng)用的可行性, 同化參數(shù)的設(shè)置等方面在內(nèi)的許多問(wèn)題有了一定的研究(J?rvinen et al., 1999;Thompson et al., 2003;Hacker and Snyder,2005)。相比于4DVAR同化方法,EnKF并不需要復(fù)雜的伴隨模式,而且可以和模式相互獨(dú)立,易于構(gòu)建和維護(hù)。同時(shí),EnKF所使用的背景誤差協(xié)方差隨模式預(yù)報(bào)而變化,具有流依賴的性質(zhì),這使得EnKF相對(duì)于3DVAR具有明顯的優(yōu)勢(shì)。但正如所有的方法都會(huì)有自身的缺點(diǎn)一樣,EnKF也存在一些尚待解決的問(wèn)題,其中有限集合數(shù)是集合資料同化存在的本質(zhì)性問(wèn)題。目前的解決方法多為采用誤差協(xié)方差的局地化和誤差協(xié)方差膨脹。

    同時(shí),隨著計(jì)算機(jī)計(jì)算能力及數(shù)值計(jì)算水平不斷提高,數(shù)值模式的分辨率越來(lái)越高,數(shù)值預(yù)報(bào)和模擬研究有較大程度的提高。中尺度數(shù)值模式空間分辨率的提高使得常規(guī)探空(高空溫、壓、濕、風(fēng)場(chǎng))資料已經(jīng)越來(lái)越不能滿足其需要。因?yàn)?00 km間距或以上探空網(wǎng)站的探測(cè)點(diǎn)稀疏,用于客觀分析的常規(guī)觀測(cè)資料密度不夠,使物理量場(chǎng)分析往往過(guò)于平滑,不能提供有足夠精度的水平物理量梯度,從而導(dǎo)致了輻散場(chǎng)、非絕熱加熱和濕度場(chǎng)之間的初始場(chǎng)缺少一致性,因而在資料客觀分析中往往喪失掉一些很重要的中尺度特征。因此若僅用常規(guī)探空資料,則可能漏掉一些中小尺度系統(tǒng)。Zhang and Fritcsh(1986)認(rèn)為若在中尺度數(shù)值模式的初始場(chǎng)中同化進(jìn)更多的中尺度信息,則能夠在一定程度上克服上述缺陷。為了適應(yīng)中尺度、短時(shí)效的精細(xì)天氣預(yù)報(bào)的需求,現(xiàn)代化的觀測(cè)手段是獲取能描述中尺度、短時(shí)天氣現(xiàn)象的觀測(cè)資料所必不可少的。除衛(wèi)星遙感、飛機(jī)觀測(cè)外,地面自動(dòng)站(AWS)觀測(cè)資料也是重要資料源之一。目前,世界上許多地區(qū)或國(guó)家都建立了分布密集的AWS觀測(cè)網(wǎng)(Shafer et al., 2000;Vejen et al., 2002)。2006 年初,全國(guó)共有5000多個(gè)地面自動(dòng)站。一些地區(qū)自動(dòng)站間的間距已小于10 km,地面觀測(cè)的空間覆蓋率遠(yuǎn)大于探空,而且自動(dòng)氣象站還提供1小時(shí)一次基本要素觀測(cè),觀測(cè)頻率非常之高,蘊(yùn)含著豐富的中尺度信息。我國(guó)的AWS觀測(cè)網(wǎng)正在受到高度重視和日趨完善,許多省市都有AWS觀測(cè)網(wǎng)。至今全國(guó)有約3萬(wàn)個(gè)AWS。AWS觀測(cè)資料具有站點(diǎn)分布密集、地形差異大、測(cè)站環(huán)境惡劣、數(shù)據(jù)采集和傳輸自動(dòng)化程度高、資料實(shí)時(shí)性強(qiáng)、中小尺度天氣現(xiàn)象明顯等特點(diǎn)(陶士偉等,2009)。

    在EnKF同化地面自動(dòng)站資料領(lǐng)域,J?rvinen et al.(1999)指出在四維同化方案中同化地面觀測(cè)可以提高預(yù)報(bào)的準(zhǔn)確性。但是,模式地形與實(shí)際地形高度不一致所帶來(lái)的問(wèn)題是不可否認(rèn)的。Hacker and Snyder(2005)闡明地面觀測(cè)富含的信息應(yīng)該被更有效地應(yīng)用在完美模式的資料同化中。Zhang and Fritcsh(1986)指出單獨(dú)同化地面觀測(cè)可以減小冬季對(duì)美國(guó)東海岸氣旋爆發(fā)的預(yù)報(bào)誤差。Stensrud et al.(2009)針對(duì)一次冷池過(guò)程,在WRF-DART系統(tǒng)中使用EnKF方法直接同化地面觀測(cè),并指出,EnKF不僅可以反映實(shí)際的地面狀況,而且可以反映邊界層以及邊界層以上,由這些觀測(cè)到的地表特征所反映的垂直運(yùn)動(dòng)和垂直結(jié)構(gòu)。這些研究說(shuō)明應(yīng)用EnKF同化地面資料對(duì)改進(jìn)邊界層預(yù)報(bào)是個(gè)富有成效的方法。

    由于地面觀測(cè)資料受地形、地貌的影響較大, 且一般模式地形與實(shí)際觀測(cè)站地形存在一定的高度差異,因此將地面觀測(cè)資料應(yīng)用到數(shù)值模式中的研究工作相對(duì)雷達(dá)衛(wèi)星等其他非常規(guī)資料的同化工作少許多(Miller and Benjamin,1992;Ruggiero et al., 1996;Urban,1996),進(jìn)展也不大。Guo et al.(2002)設(shè)計(jì)的方案沒(méi)有考慮實(shí)際觀測(cè)站地形與模式地形高度的差異,而是假定所有測(cè)站的資料(除地面氣壓)都是位于模式面,然后利用相似理論建立10 m高度風(fēng)場(chǎng)(10,10)和2 m高度溫度(2)、濕度(2)的觀測(cè)算子及相應(yīng)的切線和伴隨模式,同時(shí)在進(jìn)行極小化運(yùn)算前將地面氣壓(sfc) 折算到模式最低層。實(shí)現(xiàn)該方案采用的同化分析模式是MM5-3DVAR。此方案充分利用了觀測(cè)資料,但卻沒(méi)有考慮模式地形與實(shí)際測(cè)站地形的高度差異,對(duì)于地形分布復(fù)雜,模式地形與觀測(cè)站地形高度差異較大的區(qū)域,由于溫度、氣壓、風(fēng)場(chǎng)等觀測(cè)資料是隨著海拔高度和地形分布而變化的,差異較大就會(huì)造成模式中各物理量梯度不協(xié)調(diào)。Fujita et al.(2006)在MM5模式中應(yīng)用EnKF方法,同化風(fēng)場(chǎng)(,)、位溫()和露點(diǎn)(d),地表氣壓()只用來(lái)計(jì)算位溫(),而不進(jìn)行同化,從而降低高度差異帶來(lái)的影響,但這樣做對(duì)氣壓場(chǎng)信息的利用并不充分,同時(shí)在一定程度上造成氣壓場(chǎng)與風(fēng)場(chǎng)、溫度場(chǎng)和濕度場(chǎng)不匹配。徐枝芳等(2007a,2007b)認(rèn)為我國(guó)的地形比較復(fù)雜,模式地形與實(shí)際觀測(cè)站地形在許多地區(qū)存在較大差異,地面觀測(cè)資料同化方案設(shè)計(jì)中有必要考慮模式與實(shí)際觀測(cè)站地形高度差異,她基于 MM5_3DVAR系統(tǒng)對(duì)Guo et al. (2002)采用的方法提出了改進(jìn):在地面觀測(cè)誤差中增加由于模式地形與觀測(cè)站地形高度差異引起的地形代表性誤差。這個(gè)同化方案有效合理地將地面資料同化到了數(shù)值模式中,改進(jìn)了暴雨模擬結(jié)果,但是增加地形代表性誤差會(huì)導(dǎo)致地面觀測(cè)誤差不滿足無(wú)偏假定,分析場(chǎng)也不能達(dá)到最優(yōu)。Stensrud et al.(2009)在WRF-DART系統(tǒng)中使用EnKF方法直接同化10 m高度風(fēng)場(chǎng)(10,10)和2 m高度位溫(2)、露點(diǎn)(d2),沒(méi)有同化地表氣壓(sfc),但是同樣沒(méi)有考慮高度差異。

    本文在WRF模式中應(yīng)用EnKF方法同化自動(dòng)站資料,針對(duì)我國(guó)的復(fù)雜地形和地面自動(dòng)站的特點(diǎn),以一次暴雨個(gè)例作為研究對(duì)象,對(duì)同化方案進(jìn)行研究。首先采用Guo et al.(2002)的方案實(shí)現(xiàn)對(duì)自動(dòng)站資料的同化并檢驗(yàn)其在EnKF中表現(xiàn);然后對(duì)其進(jìn)行調(diào)整,使用位溫和露點(diǎn)代替溫度和比濕進(jìn)行同化形成新的同化方案,并檢驗(yàn)新方案的同化效果。

    2 試驗(yàn)方案設(shè)計(jì)

    2.1 集合均方根濾波(EnSRF)方法介紹

    本文應(yīng)用比較成熟的WRF模式,應(yīng)用集合平方根濾波方法。EnSRF較之基于Monte Carlo思想的傳統(tǒng)EnKF能夠避免由于擾動(dòng)觀測(cè)帶來(lái)的采樣誤差而導(dǎo)致的低估分析誤差協(xié)方差的問(wèn)題(Whitaker and Hamill,2002),由于對(duì)觀測(cè)進(jìn)行擾動(dòng)會(huì)引入額外的采樣誤差,導(dǎo)致分析誤差協(xié)方 差被低估,而如果不對(duì)觀測(cè)進(jìn)行擾動(dòng),原分析誤 差協(xié)方差就會(huì)變?yōu)?,同樣?huì)低估分析誤差協(xié)方差,其中,卡爾曼增益,為背景誤差協(xié)方差,經(jīng)典EnKF中觀測(cè)誤差協(xié)方差矩陣,為觀測(cè)算子,其作用是將模式變量轉(zhuǎn)變?yōu)橛^測(cè)變量。為了解決這個(gè)矛盾,在EnSRF方法中引入了一個(gè)小參數(shù),令,使得能夠滿足原EnKF中分析誤差協(xié)方差的公式。在單一觀測(cè)條件下的解:(Whitaker and Hamill,2002),因此,EnSRF的更新方程變?yōu)椋海?,其中代表集合的擾動(dòng),代表集合平均,代表觀測(cè)。該方法并不需要對(duì)觀測(cè)進(jìn)行擾動(dòng),更新過(guò)程中,對(duì)于集合平均采用經(jīng)典集合Kalman濾波的更新方程,而集合成員則采用所謂的“減”增益的更新方法。觀測(cè)采用順序同化方法,對(duì)觀測(cè)逐個(gè)分析,因而不涉及到矩陣轉(zhuǎn)置計(jì)算的問(wèn)題,同時(shí)EnSRF方法避免了觀測(cè)采樣誤差的引入,相對(duì)于對(duì)觀測(cè)加擾動(dòng)的“隨機(jī)”方法,這種方法也稱為“確定性”的集合Kalman濾波方法。較之同屬于確定性方法的EAKF,EnSRF方法在計(jì)算量上有著較大的優(yōu)勢(shì)(Tippett et al., 2003)。

    2.2 自動(dòng)站資料同化方案

    本文所使用的自動(dòng)站要素為10 m風(fēng)場(chǎng)(10,10)、2 m溫度(2)、2 m相對(duì)濕度(2)和地表氣壓(sfc)。設(shè)計(jì)兩種同化方案,方案一即為Guo et al.(2002)方案,方案二為對(duì)Guo et al.(2002)調(diào)整后的新方案。方案一同化的觀測(cè)類型為:10 m高度風(fēng)場(chǎng)(10,10)、2 m高度溫度(2)、2 m高度比濕(2)和地表氣壓(sfc)。方案二同化的觀測(cè)類型為10 m高度風(fēng)場(chǎng)(10,10)、2 m高度位溫(2)、2 m高度露點(diǎn)(d2)和地表氣壓(sfc)。新方案的調(diào)整主要有兩部分內(nèi)容:(1)用位溫()代替溫度()進(jìn)行同化。在白天,位溫在充分混合的邊界層,垂直分布較簡(jiǎn)單,有利于減少內(nèi)插造成的誤差。(2)用露點(diǎn)溫度(d)代替混合比()。試驗(yàn)表明,當(dāng)較小時(shí),分析階段對(duì)的更新會(huì)造成濕度變量產(chǎn)生較大的虛假增量。這很可能是由于當(dāng)發(fā)生微小變化時(shí),發(fā)生迅速變化,導(dǎo)致較強(qiáng)的非高斯誤差分布。相反露點(diǎn)溫度變化更為平滑。

    剔除-大于5倍觀測(cè)誤差的觀測(cè)資料,作為對(duì)自動(dòng)站資料的質(zhì)量控制。其中代表觀測(cè)值,代表背景場(chǎng)值。

    方案一,對(duì)10 m高度風(fēng)場(chǎng)(10,10)、2 m高度溫度(2)、2 m高度比濕(2)的同化,觀測(cè)算子包含兩部分,由垂直外推和水平差值組成。(1)垂直外推以Monin-Obukhov相似理論為基礎(chǔ),由最低半σ層預(yù)報(bào)變量得到2、2、10、10。(2)水平差值應(yīng)用線性內(nèi)插法得到相應(yīng)變量在觀測(cè)位置上的值。其中,垂直外推在WRF內(nèi)部完成并與地表物理過(guò)程參數(shù)化方案一致。

    方案二,對(duì)10 m高度風(fēng)場(chǎng)(10,10)、2 m高度位溫(2)、2 m高度露點(diǎn)(d2)的同化,觀測(cè)算子包含三部分,由垂直外推、變量變換和水平差值組成。(1)垂直外推以Monin-Obukhov相似理論為基礎(chǔ),由最低半σ層預(yù)報(bào)變量得到2、2、10、10。(2)變量變換將2、2變換為2、d2。(3)水平差值應(yīng)用線性內(nèi)插法得到相應(yīng)變量在觀測(cè)位置上的值。其中,垂直外推在WRF內(nèi)部完成并與地表物理過(guò)程參數(shù)化方案一致,變量變換中由2到2的變換在WRF內(nèi)部完成。d2由2、2、sfc計(jì)算得到。

    方案一和方案二對(duì)地表氣壓(sfc)的同化一致,模式中的sfc由WRF中直接得到,應(yīng)用靜力平衡方程將sfc由模式地形高度訂正到實(shí)際測(cè)站地形高度。同化方案對(duì)比見(jiàn)表1。

    表1 同化方案對(duì)比

    2.3 暴雨過(guò)程

    2005年9月18~21日, 副熱帶高壓邊緣暖濕氣流和北方冷空氣在山東交匯,山東出現(xiàn)了歷史上罕見(jiàn)的秋季連續(xù)暴雨過(guò)程。特別是18日夜間出現(xiàn)了以濟(jì)南為中心的東西向暴雨帶,18小時(shí)[9月18日12時(shí)至19日06時(shí)(協(xié)調(diào)世界時(shí),下同)] 降水實(shí)況如圖1所示。由于這次暴雨過(guò)程沒(méi)有出現(xiàn)明顯的強(qiáng)輻合系統(tǒng)(低渦、鋒面、氣旋、切變線等)等典型的暴雨形勢(shì)特征,屬突發(fā)區(qū)域性暴雨,預(yù)報(bào)難度較大,于是本文選這次暴雨過(guò)程作為研究個(gè)例。

    圖1 2005年9月18日12時(shí)至19日06時(shí)18小時(shí)降水觀測(cè)(單位:mm)

    2.4 試驗(yàn)設(shè)計(jì)

    本文采用WRF模式單重區(qū)域,水平格局10 km,試驗(yàn)區(qū)中心點(diǎn)取為(37°N,117.5°E),格點(diǎn)數(shù)取為200×180,垂直方向共有31層。全部試驗(yàn)的物理過(guò)程均選取了Ferrier(new Eta)對(duì)流參數(shù)化過(guò)程,MRF邊界層參數(shù)化方案,Kain-Fritsch(new Eta)微物理過(guò)程參數(shù)化方案。模式起始時(shí)間為2005年9月18日06時(shí),積分24小時(shí)?;驹囼?yàn)資料為2005年9月18日06時(shí)至19日06時(shí)每6小時(shí)一次的NCEP1°×1° FNL背景場(chǎng)資料和18日09時(shí)至18日12時(shí)逐時(shí)的山東省地面自動(dòng)站資料。

    圖2 2005年9月18日09時(shí)自動(dòng)站觀測(cè)分布情況(a)和自動(dòng)站實(shí)際地形高度與模式地形高度差(單位:m)(b)

    同化方法為EnSRF。集合數(shù)40,生成初始集合時(shí),全場(chǎng)加均值為零的隨機(jī)擾動(dòng),的擾動(dòng)標(biāo)準(zhǔn)差分別為2 m/s、2 m/s、1 m/s,的擾動(dòng)標(biāo)準(zhǔn)差為2 K,Q的擾動(dòng)標(biāo)準(zhǔn)差為0.0005 kg/kg(此處Q為背景場(chǎng)要素值)。協(xié)方差膨脹使用松弛膨脹法,將預(yù)報(bào)集合擾動(dòng)與分析集合擾動(dòng)按一定比例相加Zhang et al.(2004),背景場(chǎng)系數(shù)為0.7,分析場(chǎng)系數(shù)為0.3;局地化使用的相關(guān)函數(shù)為Gaspari and Cohn(1999)方案,使用Schur算子。水平和垂直局地化距離分別為40 km和10 km;方案一中,(10,10)、2、2、sfc觀測(cè)誤差分別為22、d2、sfc觀測(cè)誤差分別為2 m/s、2 K、2 K、1 hPa。模式起始時(shí)間為18日06時(shí),08時(shí)加入擾動(dòng),09時(shí)開(kāi)始同化第一次自動(dòng)站資料,共同化四個(gè)時(shí)次到m/s、2 K、1 g/kg、1 hPa;方案二中,(10,10)、18日12時(shí),然后開(kāi)始預(yù)報(bào)至19日06時(shí)結(jié)束。圖2為9月18日09時(shí)自動(dòng)站觀測(cè)分布情況和自動(dòng)站實(shí)際地形高度與模式地形高度差。

    為檢驗(yàn)數(shù)值預(yù)報(bào)模式對(duì)單要素物理量的響應(yīng)情況,比較不同要素物理量對(duì)模式初始場(chǎng)的改進(jìn)作用,以及比較同化單要素物理量和同化全要素物理量對(duì)數(shù)值預(yù)報(bào)模式結(jié)果的改善狀況。本文設(shè)計(jì)了一組試驗(yàn)對(duì)僅同化單要素物理量做一探討研究,為進(jìn)一步深入地開(kāi)展地面資料同化研究工作奠定基礎(chǔ)。具體設(shè)計(jì)(見(jiàn)表2)為:試驗(yàn)1為控制試驗(yàn),不做同化;試驗(yàn)2~7分別同化風(fēng)場(chǎng)、溫度、比濕、位溫、露點(diǎn)和氣壓;試驗(yàn)8為使用方案一同化所有物理量(10,10,2,2,sfc);試驗(yàn)9為使用方案二同化所有物理量(10,10,2,d2,sfc)。

    表2 試驗(yàn)方案設(shè)計(jì)

    3 EnSRF資料同化結(jié)果分析

    3.1 敏感性試驗(yàn)結(jié)果分析

    圖3為同化各要素改進(jìn)模式初始場(chǎng)后積分18小時(shí)的降水模擬結(jié)果。試驗(yàn)1(控制試驗(yàn))存在三個(gè)強(qiáng)降水中心,與實(shí)況不符,并且降水強(qiáng)度均偏小。與試驗(yàn)1對(duì)比,試驗(yàn)2~9中,除試驗(yàn)3(同化溫度)出現(xiàn)兩個(gè)強(qiáng)降水中心外,其他各試驗(yàn)均只有一個(gè)強(qiáng)降水中心,其中試驗(yàn)6(同化露點(diǎn))模擬出了120 mm的暴雨中心,其強(qiáng)度與實(shí)況較為接近。試驗(yàn)9(方案二同化所有物理量)與試驗(yàn)6相比,暴雨中心范圍有所擴(kuò)大,更加接近實(shí)況;與試驗(yàn)8(方案一同化所有物理量)相比,暴雨中心范圍與降水強(qiáng)度均更加接近實(shí)況。從以上分析可以看出,模式對(duì)同化各要素物理量均有響應(yīng),但是敏感程度各不相同。模式對(duì)露點(diǎn)最為敏感,對(duì)風(fēng)場(chǎng)、位溫和氣壓的敏感程度相似,與露點(diǎn)比相對(duì)較弱。綜合同化所有物理量對(duì)改善數(shù)值預(yù)報(bào)效果略好一些,且方案二好于方案一。

    圖3 2005年9月18日12時(shí)至19日06時(shí)18小時(shí)降水:(a)實(shí)況;(b)試驗(yàn)1;(c)試驗(yàn)2;(d)試驗(yàn)3;(e)試驗(yàn)4;(f)試驗(yàn)5;(g)試驗(yàn)6;(h)試驗(yàn)7;(i)試驗(yàn)8;(j)試驗(yàn)9

    3.2 均方根誤差分析

    為了定量分析EnSRF的同化效果,計(jì)算了各時(shí)次EnSRF分析前后背景場(chǎng)與所同化的地面自動(dòng)站要素之間的均方根誤差(RMSE)。圖4給出了使用方案一同化地面自動(dòng)站資料(10,10,2,2,sfc)和使用方案二同化地面自動(dòng)站資料(10,10,2,d2,sfc)過(guò)程中,10 m風(fēng)場(chǎng)(10)、10 m風(fēng)場(chǎng)(10)、2 m溫度(2)、2 m比濕(2)、2 m位溫(2)、2 m露點(diǎn)(d2)和地表氣壓(sfc)的RMSE隨時(shí)刻的變化。

    由圖4可以看出,方案一對(duì)2的同化,分析后的RMSE比預(yù)報(bào)結(jié)果的RMSE要大,這種負(fù)效果可能是由于2的劇烈變化產(chǎn)生的較強(qiáng)非高斯誤差分布(Fujita et al., 2006)導(dǎo)致的;而方案二,在同化過(guò)程中使用d2代替2,能較好的改進(jìn)對(duì)濕度觀測(cè)的同化。除方案一中的2外,每次分析結(jié)果的RMSE都比預(yù)報(bào)結(jié)果的RMSE小,這表明每次同化地面自動(dòng)站資料后的結(jié)果都比同化前更接近實(shí)際觀測(cè),體現(xiàn)了EnSRF同化的有效性。同時(shí),RMSE_2與RMSE_2相比,值相對(duì)較小,下降趨勢(shì)更明顯,說(shuō)明使用2代替2能夠使結(jié)果更接近實(shí)際觀測(cè)。方案二中的RMSE_10和RMSE_10比方案一中的要小,說(shuō)明使用方案二比使用方案一對(duì)風(fēng)場(chǎng)資料的同化更有效。方案一與方案二的RMSE_sfc差別較小,說(shuō)明對(duì)氣壓場(chǎng)的同化效果差別不大。在方案二中,RMSE_10和RMSE_10隨著時(shí)間的增加變化不明顯,09時(shí)預(yù)報(bào)的RMSE分別為1.17 m/s和1.32 m/s,12時(shí)分析后分別為0.95 m/s和0.99 m/s;RMSE_2、RMSE_d2和RMSE_sfc隨著時(shí)間的增加迅速減小,09時(shí)預(yù)報(bào)的RMSE分別為1.846434 K、3.04816 K和1.142898 hPa, 12時(shí)分析后分別為0.970432 K、1.32523 K和0.317531 hPa。

    3.3 觀測(cè)有效利用情況

    在EnSRF系統(tǒng)同化地面自動(dòng)站資料時(shí),做了簡(jiǎn)單的質(zhì)量控制,當(dāng)觀測(cè)值與模式值之間的差值大于5倍觀測(cè)誤差的時(shí)候,剔除該觀測(cè)資料,不參與同化分析。圖5分別為為方案一和方案二觀測(cè)參與同化分析的情況。

    圖5 觀測(cè)參與同化分析的情況:(a)方案一;(b)方案二

    在4個(gè)時(shí)次的自動(dòng)站資料中,共包含455組(group)資料,方案一同化過(guò)程中,比濕(2)、風(fēng)場(chǎng)(10,10)觀測(cè)全部參與了同化分析,參與同化分析的溫度(2)觀測(cè)有454個(gè),氣壓(sfc)觀測(cè)有441個(gè)。方案二同化過(guò)程中位溫(2)、露點(diǎn)(d2)、風(fēng)場(chǎng)(10,10)觀測(cè)全部參與了同化分析,氣壓(sfc)觀測(cè)有441個(gè)參與了同化分析。

    3.4 同化增量分析

    同化增量為同化分析結(jié)果減去背景場(chǎng)之差。圖6、圖7、圖8分別為模式第二層()關(guān)于風(fēng)場(chǎng)的散度場(chǎng)增量、溫度場(chǎng)增量和濕度場(chǎng)增量。從散度場(chǎng)增量圖上可以看出,試驗(yàn)8山東北部地區(qū)及山東中部及偏西南地區(qū)為輻合區(qū),試驗(yàn)9較試驗(yàn)8在(35°N,116.5°E)附近輻合中心較強(qiáng)。在溫度場(chǎng)增量圖上,試驗(yàn)8在山東中北部地區(qū)為增溫區(qū),試驗(yàn)9較試驗(yàn)8增溫中心位置偏西,且增溫中心強(qiáng)度偏大。在濕度場(chǎng)增量圖上,試驗(yàn)8山東中北部地區(qū)為增濕區(qū),在山東西北部、南部及西南部為減濕 區(qū),試驗(yàn)9較試驗(yàn)8增濕區(qū)位置偏南,且增濕幅度偏大,西北部減濕區(qū)減濕幅度偏小。由同化單要素物理量的敏感性試驗(yàn)結(jié)論可知濕度場(chǎng)對(duì)初始場(chǎng)的影響作用最大,因此試驗(yàn)8和試驗(yàn)9降水區(qū)偏東北的主要原因可能是初始場(chǎng)中增濕區(qū)偏東北。

    圖6 模式第二層水平散度場(chǎng)增量(單位:10?5 s?1):(a)試驗(yàn)8;(b)試驗(yàn)9

    圖7 模式第二層溫度場(chǎng)增量(單位:K):(a)試驗(yàn)8;(b)試驗(yàn)9

    圖8 模式第二層濕度場(chǎng)增量(單位:g/kg):(a)試驗(yàn)8;(b)試驗(yàn)9

    4 總結(jié)與討論

    本文將Guo et al.(2002)的方案應(yīng)用在EnSRF系統(tǒng)中,并對(duì)其進(jìn)行了調(diào)整形成了新方案,針對(duì)一次暴雨過(guò)程進(jìn)行了自動(dòng)站資料同化分析研究。結(jié)果表明:地面觀測(cè)資料(溫度、濕度、風(fēng)場(chǎng)、地面氣壓)中各物理量分別同化到數(shù)值模式都能影響18小時(shí)降水預(yù)報(bào),但各物理量所起作用大小不同,使用位溫代替溫度參與同化以及使用露點(diǎn)代替比濕進(jìn)行同化都能改數(shù)值結(jié)果,其中對(duì)結(jié)果影響最大的是露點(diǎn);在EnSRF同化系統(tǒng)中應(yīng)用Guo et al.(2002)的方案將地面自動(dòng)站資料進(jìn)行同化到數(shù)值模式中,能夠部分改善模擬結(jié)果;新方案較Guo et al.(2002)的方案10 m風(fēng)場(chǎng)的均方根誤差更小,說(shuō)明對(duì)地面溫度、濕度資料同化分析過(guò)程中除影響自身的分析場(chǎng),同時(shí)還影響風(fēng)場(chǎng)的分析;新方案對(duì)溫度資料的利用也更加充分??偟膩?lái)說(shuō),新方案較Guo et al.(2002)的方案的模擬結(jié)果更加接近實(shí)況,對(duì)自動(dòng)站資料的同化也更加有效。

    本文使用了風(fēng)、溫度、濕度與氣壓數(shù)據(jù)進(jìn)行同化研究,并分析了各類數(shù)據(jù)的貢獻(xiàn),新方案對(duì)促進(jìn)地面自動(dòng)站觀測(cè)在數(shù)值預(yù)報(bào)中的應(yīng)用有一定意義。EnSRF方法及新方案對(duì)于同化地面資料有較好的應(yīng)用前景,也為充分合理地利用中尺度信息、解決中尺度問(wèn)題提供了參考和建議。

    本文雖然取得了一定的成果,但是仍存在不足。在我國(guó)地形復(fù)雜,模式地形與實(shí)際測(cè)站地形在一些地區(qū)存在較大差異的背景下,除對(duì)氣壓的同化之外,Guo et al.(2002)的方案和本文中的新方案均沒(méi)有考慮這種差異,需更深入的研究。

    致謝 感謝兩位匿名審稿專家及編委對(duì)本文提出的寶貴意見(jiàn)。

    (References:)

    Barker D M, Huang W, Guo Y R, et al. 2004. A Three Dimensional Variational (3DVAR) data assimilation system for use with MM5: Implementation and initial results [J]. Mon. Wea. Rev., 132 (4): 897–914.

    De Pondeca M S F V, Zou X L. 2001. A case study of the variational assimilation of GPS zenith delay observations into a mesoscale model [J]. J. Appl. Meteor., 40 (9): 1559–1576.

    Evensen G. 1994. Sequential data assimilation with a non-linear quasi-geostrophic model using Monte Carlo methods to forecast error statistics [J]. J. Geophys. Res., 99 (C5): 10143–10162.

    Fujita T, Stensrud D J, Dowell D C. 2006. Surface data assimilation using an ensemble Kalman filter approach with initial condition and model physics uncertainties [J]. Mon. Wea. Rev., 135 (5): 1846–1868.

    Gao J, Xue M, Brewster K, et al. 2001. A 3DVAR data assimilation scheme for storm-scale model [C]. 14th Conference on NWP, AMS, J: 72–74.

    Guo Y R, Shin D H, Lee J H, et al. 2002. Application of the MM5 3DVAR system for a heavy rain case over the Korean Peninsula [C]// Papers Presented at the Twelfth PSU/NCAR Mesoscale Model Users’ Workshop NCAR, June 24–25, 2002.

    Hacker J P, Snyder C. 2005. Ensemble Kalman filter assimilation of fixed screen-height observations in a parameterized PBL [J]. Mon. Wea. Rev., 133 (11): 3260–3275.

    Huang X Y, Xiao Q, Barker D M, et al. 2009. Four-dimensional variational data assimilation for WRF: Formulation and preliminary results [J]. Mon. Wea. Rev., 137 (1): 299–314.

    J?rvinen H, Andersson E, Bouttier F. 1999. Variational assimilation of time sequences of surface observations with serially correlated errors [J]. Tellus, 51 (4): 469–488.

    Martinelli J T, Pasken R W, Lin Y J, et al. 2003. A high resolution numerical simulation of a linear mesoscale convective system utilizing the MM5 4DVar system and single WSR-88D data [C]. 31st International Conference on Radar Meteorology American Meteorology Society, P1C. 3, Abstract and paper.

    Miller P A, Benjamin S G. 1992. A system for the hourly assimilation of surface observations in mountainous and flat terrain [J]. Mon. Wea. Rev., 120 (10): 2342–2359.

    Qiu C J, Chou J F. 2006. Four-dimensional data assimilation method based on SVD: Theoretical aspect [J]. Theor. Appl. Climatol., 83 (1–4): 51–57.

    Ruggiero F H, Sashegyi K D, Madala R V, et al. 1996. The use of surface observations in four dimensional data assimilation using a mesoscale model [J]. Mon. Wea. Rev., 124 (5): 1018–1033.

    Shafer M A, Fiebrich C A, Arent D S, et al. 2000. Quality assurance procedures in the Oklahoma Mseonet-work [J]. J. Atmos. Oceanic Techenol., 17 (4): 474–494.

    Stensrud D J, Yussouf N, Dowell D C, et al. 2009. Assimilating surface data into a mesoscale model ensemble: Cold pool analyses from spring 2007 [J]. Atmospheric Research, 93 (1–3): 207–220.

    陶士偉, 仲躋芹, 徐枝芳, 等. 2009. 地面自動(dòng)站資料控制方案及應(yīng)用 [J]. 高原氣象, 28 (5): 1202–1210. Tao Shiwei, Zhong Qiqin, Xu Zhifang, et al. 2009. Quality control schemes and its application to automatic surface weather observation system [J]. Plateau Meteorology (in Chinese), 28 (5): 1202-1210.

    Thompson R L, Edwards R, Hart J A, et al. 2003. Close proximity soundings within supercell environments obtained from the rapid update cycle [J]. Wea. Forecasting,18: 1243–1261.

    Tippett M K, Anderson J L, Bishop C H, et al. 2003. Ensemble square-root filters [J]. Mon. Wea. Rev., 131: 1485–1490.

    Urban B. 1996. Coherent observation operators for surface data assimilation with application to snow depth [J]. J. Appl. Meteor., 35 (2): 258–270.

    Vejen F, Jacobsson C, Fredriksson U, et al. 2002. Quality control of meteorological observations automatic methods used in the Nordic countries [R]. Climate Report, No. 8, KLIMA.

    Whitaker J S, Hamill T M. 2002. Ensemble data assimilation without perturbed observations [J]. Mon. Wea. Rev., 130 (7): 1913–1924.

    徐枝芳, 龔建東, 王建捷, 等. 2007a. 復(fù)雜地形下地面觀測(cè)資料同化Ⅰ. 模式地形與觀測(cè)站地形高度差異對(duì)地面資料同化的影響研究 [J]. 大氣科學(xué), 31 (2): 222–232. Xu Zhifang, Gong Jiandong, Wang Jianjie, et al. 2007a. A study of assimilation of surface observational data in complex terrain. Part Ⅰ: Influence of the elevation difference between model surface and observation site [J]. Chinese J. Atmos. Sci. (in Chinese), 31 (2): 222-232.

    徐枝芳, 龔建東, 王建捷, 等. 2007b. 復(fù)雜地形下地面觀測(cè)資料同化Ⅱ. 模式地形與觀測(cè)站地形高度差異代表性誤差 [J]. 大氣科學(xué), 31 (3): 449–458. Xu Zhifang, Gong Jiandong, Wang Jianjie, et al. 2007b. A study of assimilation of surface observational data in complex terrain. Part Ⅱ: Representative error of the elevation difference between model surface and observation site [J]. Chinese J. Atmos. Sci. (in Chinese), 31 (3): 449-458.

    薛紀(jì)善, 莊世宇, 朱國(guó)富, 等. 2008. GRAPES新一代全球/區(qū)域變分同化系統(tǒng)研究 [J]. 科學(xué)通報(bào), 53 (20): 2408–2417. Xue Jshan, Zhuang Shiyu, Zhu Guofu, et al. 2008. Research of new generation global/region variational assimilation system GRAPES [J]. Chin. Sci. Bull. (in Chinese), 53 (20): 2408-2417.

    Zhang D L, Fritcsh J M. 1986. A case study of the sensitivity of numerical model simulation of mesoscale convective systems to varying initial condition [J]. Mon. Wea. Rev., 114 (12): 2481–2431.

    邵長(zhǎng)亮, 閔錦忠. 2015. 集合均方根濾波同化地面自動(dòng)站資料的技術(shù)研究[J]. 大氣科學(xué), 39 (1): 1?11, doi:10.3878/j.issn.1006-9895.1406.13263. Shao Changliang, Min Jinzhong. 2015. A study of the assimilation of surface automatic weather station data using the ensemble square root filter [J]. Chinese Journal of Atmospheric Sciences (in Chinese), 39 (1): 1?11

    A Study of the Assimilation of Surface Automatic Weather Station Data Using the Ensemble Square Root Filter

    SHAO Changliang1, 2, 3and MIN Jinzhong1, 3

    1,,210044;2,100081;3,210044

    Handling the difference in elevation between a model surface and an observation site is always a challenge in surface data assimilation. However, a reasonable assimilation scheme can efficiently assimilate surface automatic weather station (AWS) data into a mesoscale model. In this paper, surface AWS data are first assimilated into a weather research and forecasting (WRF) model through an ensemble Kalman filter using the Guo et al. (2002) scheme. Then an adjusted scheme is proposed that assimilates 10-m wind observations, 2-m potential temperature, 2-m dew point temperature, and surface pressure. This scheme is then validated by mean square root error analysis, simulated result and assimilation increment analysis, and sensitive experiments to check the assimilation response of each AWS meteorological parameter. Results show that the assimilation of surface AWS data through the ensemble square root filter (EnSRF) using the Guo et al. (2002) scheme can improve the simulation results. The separate assimilation of any element of the surface observation data (temperature, humidity, wind, surface pressure) can affect the forecast of 18 h accumulated rainfall. However, different elements have different impacts, and the one having most influence is the dew point temperature. The use of 2-m potential temperature and 2-m dew point temperature, instead of 2-m temperature and 2-m specific humidity, leads to better simulation results.

    Data Assimilation, Ensemble Kalman filter, AWS (automatic weather station) data, Sensitive experiments

    1006?9895(2015)01?0001?11

    P456.7

    A

    10.3878/j.issn.1006-9895.1406.13263

    2013?09?16;網(wǎng)絡(luò)預(yù)出版日期 2014?07?08

    國(guó)家重點(diǎn)基礎(chǔ)研究發(fā)展計(jì)劃(973計(jì)劃)2013CB430102,江蘇省普通高校研究生科研創(chuàng)新計(jì)劃項(xiàng)目KYLX_0824

    邵長(zhǎng)亮,男,1986年出生,博士研究生,工程師,主要從事自動(dòng)站資料同化的研究。E-mail:shchl1@163.com

    閔錦忠,E-mail: minjz@nuist.edu.cn

    丰满迷人的少妇在线观看| 涩涩av久久男人的天堂| 露出奶头的视频| 亚洲五月婷婷丁香| 91老司机精品| 超碰97精品在线观看| 欧美变态另类bdsm刘玥| 国产成人影院久久av| 老司机午夜福利在线观看视频 | 热re99久久精品国产66热6| 天天影视国产精品| 91av网站免费观看| 精品欧美一区二区三区在线| 色婷婷久久久亚洲欧美| 久久国产精品男人的天堂亚洲| 国产亚洲欧美在线一区二区| 精品熟女少妇八av免费久了| 99精品在免费线老司机午夜| 少妇粗大呻吟视频| 国产午夜精品久久久久久| 大香蕉久久成人网| 黄色怎么调成土黄色| 中文亚洲av片在线观看爽 | 久久久久久久久久久久大奶| 免费av中文字幕在线| 欧美激情久久久久久爽电影 | 国产高清videossex| 女性被躁到高潮视频| 91成人精品电影| 日韩熟女老妇一区二区性免费视频| 久久九九热精品免费| 搡老乐熟女国产| 黄色成人免费大全| a级片在线免费高清观看视频| 久久亚洲真实| 久久国产精品人妻蜜桃| 丰满饥渴人妻一区二区三| 在线永久观看黄色视频| 久久精品91无色码中文字幕| 两个人看的免费小视频| 热re99久久精品国产66热6| 久久中文看片网| 中国美女看黄片| 多毛熟女@视频| 一边摸一边抽搐一进一出视频| 欧美av亚洲av综合av国产av| 高清黄色对白视频在线免费看| 成人国语在线视频| 欧美 亚洲 国产 日韩一| 亚洲av成人不卡在线观看播放网| 99精国产麻豆久久婷婷| av又黄又爽大尺度在线免费看| 成人永久免费在线观看视频 | 黄色丝袜av网址大全| 国产精品一区二区精品视频观看| 国产欧美日韩一区二区三| a级毛片在线看网站| 国产精品免费大片| 亚洲欧美激情在线| 久久久水蜜桃国产精品网| 中文字幕人妻熟女乱码| 在线观看舔阴道视频| 99国产极品粉嫩在线观看| 老司机影院毛片| 老司机午夜十八禁免费视频| 亚洲av成人一区二区三| 国产日韩欧美在线精品| 国产伦人伦偷精品视频| 一级毛片电影观看| 777久久人妻少妇嫩草av网站| 好男人电影高清在线观看| 老司机靠b影院| 亚洲avbb在线观看| 黑人猛操日本美女一级片| 亚洲人成伊人成综合网2020| 日本撒尿小便嘘嘘汇集6| 在线av久久热| 777米奇影视久久| 黄色a级毛片大全视频| 美国免费a级毛片| 亚洲欧美激情在线| 亚洲美女黄片视频| 18禁国产床啪视频网站| 老司机影院毛片| 麻豆av在线久日| 久久这里只有精品19| 多毛熟女@视频| 一本大道久久a久久精品| 欧美+亚洲+日韩+国产| 成人免费观看视频高清| 777米奇影视久久| 下体分泌物呈黄色| 国产在线视频一区二区| 亚洲国产av新网站| 手机成人av网站| 国产欧美日韩一区二区三| 91国产中文字幕| 日本av手机在线免费观看| 国产高清视频在线播放一区| 成年动漫av网址| 香蕉丝袜av| 欧美人与性动交α欧美精品济南到| 成人18禁高潮啪啪吃奶动态图| 一边摸一边抽搐一进一出视频| av免费在线观看网站| 新久久久久国产一级毛片| 国产91精品成人一区二区三区 | 亚洲avbb在线观看| 中文字幕精品免费在线观看视频| 一区二区三区激情视频| 午夜福利欧美成人| 国内毛片毛片毛片毛片毛片| 在线观看免费午夜福利视频| 丰满迷人的少妇在线观看| netflix在线观看网站| 婷婷成人精品国产| 男女无遮挡免费网站观看| 亚洲精品美女久久av网站| 人成视频在线观看免费观看| 亚洲人成电影免费在线| 亚洲精品国产色婷婷电影| 一本久久精品| 国产一区二区三区视频了| 国产精品香港三级国产av潘金莲| 19禁男女啪啪无遮挡网站| 亚洲黑人精品在线| 精品国产一区二区久久| 最近最新中文字幕大全电影3 | 99久久国产精品久久久| 香蕉久久夜色| 考比视频在线观看| 久久久精品国产亚洲av高清涩受| 美女国产高潮福利片在线看| 日本黄色日本黄色录像| 亚洲国产欧美在线一区| 黑人巨大精品欧美一区二区mp4| 欧美黑人欧美精品刺激| 亚洲成国产人片在线观看| www日本在线高清视频| 日韩欧美一区二区三区在线观看 | 欧美日韩福利视频一区二区| 欧美激情久久久久久爽电影 | 脱女人内裤的视频| 日本欧美视频一区| 别揉我奶头~嗯~啊~动态视频| 王馨瑶露胸无遮挡在线观看| 成年动漫av网址| 欧美精品一区二区免费开放| 色播在线永久视频| 国产免费av片在线观看野外av| 国产av又大| 午夜91福利影院| 丁香欧美五月| 最新的欧美精品一区二区| 国产99久久九九免费精品| 亚洲精品在线观看二区| 18在线观看网站| 91精品国产国语对白视频| 午夜福利视频在线观看免费| 亚洲色图 男人天堂 中文字幕| 亚洲国产欧美网| 亚洲第一av免费看| 中文欧美无线码| 一区二区三区国产精品乱码| 啦啦啦免费观看视频1| 亚洲精华国产精华精| 欧美久久黑人一区二区| 久久精品国产a三级三级三级| 肉色欧美久久久久久久蜜桃| √禁漫天堂资源中文www| 777久久人妻少妇嫩草av网站| 捣出白浆h1v1| 久久久国产一区二区| 国产精品99久久99久久久不卡| 日韩欧美一区视频在线观看| 香蕉国产在线看| 国产亚洲午夜精品一区二区久久| 欧美亚洲日本最大视频资源| 国产免费现黄频在线看| 免费观看a级毛片全部| 成人永久免费在线观看视频 | 国产免费视频播放在线视频| 天天躁狠狠躁夜夜躁狠狠躁| 久9热在线精品视频| 一区二区三区乱码不卡18| 考比视频在线观看| 操出白浆在线播放| 国产亚洲av高清不卡| 视频区欧美日本亚洲| 欧美另类亚洲清纯唯美| 亚洲少妇的诱惑av| 国产精品av久久久久免费| 日本wwww免费看| 色综合婷婷激情| 国产在线一区二区三区精| 黄频高清免费视频| 成年女人毛片免费观看观看9 | 脱女人内裤的视频| 精品少妇一区二区三区视频日本电影| 99re在线观看精品视频| 久久国产精品大桥未久av| 精品熟女少妇八av免费久了| 国产精品一区二区免费欧美| 久热爱精品视频在线9| 国产又色又爽无遮挡免费看| 一个人免费在线观看的高清视频| 青草久久国产| 男女高潮啪啪啪动态图| 一本综合久久免费| 嫩草影视91久久| 视频区欧美日本亚洲| 日本一区二区免费在线视频| 午夜免费鲁丝| 九色亚洲精品在线播放| 巨乳人妻的诱惑在线观看| 亚洲av日韩在线播放| av网站在线播放免费| 一区二区三区乱码不卡18| av天堂在线播放| 国产av精品麻豆| 女人精品久久久久毛片| 免费在线观看视频国产中文字幕亚洲| 亚洲少妇的诱惑av| 久久精品91无色码中文字幕| 热99国产精品久久久久久7| 99re6热这里在线精品视频| 亚洲成国产人片在线观看| 超碰97精品在线观看| 亚洲欧美一区二区三区黑人| 亚洲av成人一区二区三| 80岁老熟妇乱子伦牲交| 亚洲成a人片在线一区二区| 成年人午夜在线观看视频| 美女高潮到喷水免费观看| 国产亚洲午夜精品一区二区久久| 国产深夜福利视频在线观看| 香蕉国产在线看| 久久久久网色| 国产精品一区二区在线观看99| 国产精品欧美亚洲77777| 亚洲 国产 在线| av网站免费在线观看视频| 精品高清国产在线一区| 亚洲精品美女久久av网站| 国产成人啪精品午夜网站| 久久人人97超碰香蕉20202| 亚洲精华国产精华精| 欧美日韩国产mv在线观看视频| 脱女人内裤的视频| avwww免费| 久久久精品国产亚洲av高清涩受| 欧美 日韩 精品 国产| 欧美精品人与动牲交sv欧美| 久久精品国产综合久久久| 欧美在线黄色| av片东京热男人的天堂| 国产日韩欧美视频二区| 中文字幕人妻丝袜制服| 男男h啪啪无遮挡| 久热这里只有精品99| 一区二区av电影网| 国产成人欧美| 电影成人av| 美国免费a级毛片| 成年人免费黄色播放视频| 丰满人妻熟妇乱又伦精品不卡| 精品亚洲成国产av| 亚洲国产av影院在线观看| 精品国产超薄肉色丝袜足j| 精品免费久久久久久久清纯 | 亚洲成人免费电影在线观看| 人成视频在线观看免费观看| 色婷婷久久久亚洲欧美| 黄频高清免费视频| 日韩免费高清中文字幕av| 日韩三级视频一区二区三区| 黑人欧美特级aaaaaa片| 欧美老熟妇乱子伦牲交| 熟女少妇亚洲综合色aaa.| 涩涩av久久男人的天堂| 性少妇av在线| 久久精品国产a三级三级三级| 成年人黄色毛片网站| 欧美日韩精品网址| 免费观看av网站的网址| 成人永久免费在线观看视频 | 曰老女人黄片| 一夜夜www| 大片免费播放器 马上看| 深夜精品福利| 精品一区二区三区av网在线观看 | 女人精品久久久久毛片| 午夜精品国产一区二区电影| 最黄视频免费看| 建设人人有责人人尽责人人享有的| 69av精品久久久久久 | 一级,二级,三级黄色视频| 香蕉久久夜色| 丝袜美腿诱惑在线| 91精品国产国语对白视频| 欧美黑人欧美精品刺激| 人人澡人人妻人| 一区二区日韩欧美中文字幕| 欧美另类亚洲清纯唯美| 国产精品香港三级国产av潘金莲| av网站在线播放免费| 丁香六月天网| av片东京热男人的天堂| 涩涩av久久男人的天堂| 女人精品久久久久毛片| 在线观看人妻少妇| 亚洲成a人片在线一区二区| 天天躁夜夜躁狠狠躁躁| 亚洲人成电影免费在线| 国产精品久久久久久精品古装| 麻豆av在线久日| 欧美日韩av久久| 亚洲欧美一区二区三区黑人| 一本大道久久a久久精品| 精品久久久久久久毛片微露脸| 91老司机精品| 一二三四社区在线视频社区8| 久久久久精品国产欧美久久久| 男女无遮挡免费网站观看| 一本—道久久a久久精品蜜桃钙片| 老熟女久久久| 日本vs欧美在线观看视频| 自线自在国产av| 免费观看a级毛片全部| 动漫黄色视频在线观看| 国产av又大| 一区二区三区激情视频| 国产精品 国内视频| 大片电影免费在线观看免费| 最新美女视频免费是黄的| 热99re8久久精品国产| svipshipincom国产片| 亚洲熟女毛片儿| 性色av乱码一区二区三区2| 亚洲国产欧美网| 国产精品.久久久| 日日夜夜操网爽| 一二三四社区在线视频社区8| 亚洲情色 制服丝袜| 欧美老熟妇乱子伦牲交| 淫妇啪啪啪对白视频| 久久av网站| 久久热在线av| 精品少妇黑人巨大在线播放| 久久99热这里只频精品6学生| 美女高潮喷水抽搐中文字幕| 脱女人内裤的视频| av有码第一页| 中文字幕av电影在线播放| 精品少妇黑人巨大在线播放| 日日摸夜夜添夜夜添小说| 精品少妇黑人巨大在线播放| 黄色毛片三级朝国网站| 久久精品人人爽人人爽视色| 精品国产国语对白av| 国产男女超爽视频在线观看| 女人爽到高潮嗷嗷叫在线视频| 国产av精品麻豆| 在线观看舔阴道视频| 久久精品人人爽人人爽视色| 香蕉久久夜色| 高清黄色对白视频在线免费看| 成年人午夜在线观看视频| 国产福利在线免费观看视频| 手机成人av网站| 久久久久久久久久久久大奶| 麻豆乱淫一区二区| 成人黄色视频免费在线看| 欧美日韩黄片免| 日本av免费视频播放| 老汉色∧v一级毛片| 男人操女人黄网站| 日韩免费av在线播放| 精品熟女少妇八av免费久了| 国产精品影院久久| 精品久久蜜臀av无| 国产精品香港三级国产av潘金莲| 国产高清videossex| 91麻豆精品激情在线观看国产 | 久久九九热精品免费| 婷婷成人精品国产| 无人区码免费观看不卡 | 国产精品免费大片| 国产不卡av网站在线观看| 美女扒开内裤让男人捅视频| 亚洲伊人色综图| 欧美成狂野欧美在线观看| 国产成人精品久久二区二区免费| 我要看黄色一级片免费的| 国产欧美日韩精品亚洲av| 男女高潮啪啪啪动态图| av在线播放免费不卡| 亚洲成a人片在线一区二区| 欧美亚洲日本最大视频资源| 亚洲精品自拍成人| 欧美性长视频在线观看| 久久国产精品大桥未久av| 脱女人内裤的视频| 热99久久久久精品小说推荐| 精品久久蜜臀av无| 欧美+亚洲+日韩+国产| 亚洲精品自拍成人| av电影中文网址| 亚洲成人国产一区在线观看| 精品欧美一区二区三区在线| 日韩制服丝袜自拍偷拍| 欧美乱妇无乱码| 国产精品熟女久久久久浪| 岛国毛片在线播放| 中文字幕最新亚洲高清| 免费黄频网站在线观看国产| 精品熟女少妇八av免费久了| 亚洲欧美激情在线| 69av精品久久久久久 | 国产欧美日韩一区二区三| 狂野欧美激情性xxxx| www.自偷自拍.com| 欧美成狂野欧美在线观看| 亚洲av欧美aⅴ国产| 9191精品国产免费久久| 天堂动漫精品| 男女高潮啪啪啪动态图| 亚洲国产毛片av蜜桃av| 亚洲伊人久久精品综合| 久久人人97超碰香蕉20202| 日韩一卡2卡3卡4卡2021年| 国产在线观看jvid| 一级毛片电影观看| 91精品三级在线观看| 一级毛片精品| 国内毛片毛片毛片毛片毛片| 亚洲精品一卡2卡三卡4卡5卡| 国产精品偷伦视频观看了| 免费在线观看影片大全网站| 精品一区二区三区四区五区乱码| 美女午夜性视频免费| 国产精品欧美亚洲77777| 美女高潮喷水抽搐中文字幕| 亚洲成人手机| 99国产精品一区二区蜜桃av | kizo精华| 亚洲色图综合在线观看| 国产亚洲精品久久久久5区| 777久久人妻少妇嫩草av网站| 久久精品亚洲av国产电影网| 日韩欧美国产一区二区入口| 欧美一级毛片孕妇| 国产精品 欧美亚洲| 一区福利在线观看| 亚洲精品国产色婷婷电影| 精品少妇黑人巨大在线播放| 欧美乱码精品一区二区三区| 天天操日日干夜夜撸| 亚洲五月色婷婷综合| 午夜精品久久久久久毛片777| 日韩制服丝袜自拍偷拍| 2018国产大陆天天弄谢| 亚洲成av片中文字幕在线观看| 亚洲中文日韩欧美视频| cao死你这个sao货| 激情视频va一区二区三区| 热99re8久久精品国产| 日本欧美视频一区| 国产在线精品亚洲第一网站| 大陆偷拍与自拍| 精品福利观看| 中文字幕人妻熟女乱码| 日本欧美视频一区| 青青草视频在线视频观看| 91成年电影在线观看| 亚洲va日本ⅴa欧美va伊人久久| 这个男人来自地球电影免费观看| 久久人人97超碰香蕉20202| 久久久久国产一级毛片高清牌| 国产欧美日韩一区二区三| 中文字幕av电影在线播放| 久久人妻熟女aⅴ| 变态另类成人亚洲欧美熟女 | 日韩 欧美 亚洲 中文字幕| 精品第一国产精品| 大陆偷拍与自拍| 九色亚洲精品在线播放| 1024视频免费在线观看| 黄色成人免费大全| 国产av精品麻豆| 国产91精品成人一区二区三区 | 精品国产超薄肉色丝袜足j| 亚洲国产欧美一区二区综合| a级毛片黄视频| 母亲3免费完整高清在线观看| 欧美人与性动交α欧美软件| 久久精品人人爽人人爽视色| 91精品国产国语对白视频| 久久国产精品大桥未久av| 精品久久久精品久久久| 日本精品一区二区三区蜜桃| 欧美中文综合在线视频| 免费黄频网站在线观看国产| 国产极品粉嫩免费观看在线| 999久久久国产精品视频| 精品熟女少妇八av免费久了| 午夜视频精品福利| 中文字幕人妻熟女乱码| 少妇粗大呻吟视频| 99精品久久久久人妻精品| 亚洲精品国产区一区二| 久久精品亚洲av国产电影网| 精品国产亚洲在线| 嫁个100分男人电影在线观看| 日韩一区二区三区影片| 免费人妻精品一区二区三区视频| 亚洲五月色婷婷综合| 嫁个100分男人电影在线观看| 丰满饥渴人妻一区二区三| 欧美乱妇无乱码| 飞空精品影院首页| 欧美日韩福利视频一区二区| 自线自在国产av| 日本av免费视频播放| 侵犯人妻中文字幕一二三四区| 午夜视频精品福利| 久久精品aⅴ一区二区三区四区| 女人精品久久久久毛片| 亚洲国产成人一精品久久久| 高清毛片免费观看视频网站 | 亚洲成av片中文字幕在线观看| 最新美女视频免费是黄的| 高清视频免费观看一区二区| 精品一区二区三区av网在线观看 | 日本一区二区免费在线视频| 欧美+亚洲+日韩+国产| 国产区一区二久久| 啦啦啦免费观看视频1| 国产成人啪精品午夜网站| 91麻豆精品激情在线观看国产 | 中文字幕人妻丝袜制服| 亚洲成人免费电影在线观看| 一本色道久久久久久精品综合| 久久av网站| 高清欧美精品videossex| 如日韩欧美国产精品一区二区三区| 日韩成人在线观看一区二区三区| 91老司机精品| 满18在线观看网站| 999久久久国产精品视频| √禁漫天堂资源中文www| 久久久国产精品麻豆| 侵犯人妻中文字幕一二三四区| 极品人妻少妇av视频| 国产精品免费一区二区三区在线 | 精品国产亚洲在线| 亚洲第一欧美日韩一区二区三区 | 正在播放国产对白刺激| 国产精品av久久久久免费| 亚洲精品在线美女| 久久精品国产99精品国产亚洲性色 | 99国产精品免费福利视频| 色婷婷久久久亚洲欧美| 国产欧美日韩综合在线一区二区| 老熟女久久久| 亚洲国产欧美在线一区| 欧美变态另类bdsm刘玥| 女人精品久久久久毛片| 日日夜夜操网爽| 国产欧美亚洲国产| 热99国产精品久久久久久7| 电影成人av| 成人18禁高潮啪啪吃奶动态图| 亚洲成av片中文字幕在线观看| 99九九在线精品视频| 国产福利在线免费观看视频| 精品国产亚洲在线| 99香蕉大伊视频| 亚洲 国产 在线| 操美女的视频在线观看| 涩涩av久久男人的天堂| 少妇的丰满在线观看| 久久亚洲精品不卡| 久久人妻av系列| 高清视频免费观看一区二区| 国精品久久久久久国模美| 日本精品一区二区三区蜜桃| 亚洲五月婷婷丁香| 高清欧美精品videossex| 嫁个100分男人电影在线观看| 97在线人人人人妻| av一本久久久久| 99re6热这里在线精品视频| 国产av国产精品国产| 色94色欧美一区二区| 丝袜美足系列| 久久99一区二区三区| tocl精华| 久久婷婷成人综合色麻豆| 无人区码免费观看不卡 | 纵有疾风起免费观看全集完整版| 日本黄色视频三级网站网址 | av网站在线播放免费| 欧美激情 高清一区二区三区| 看免费av毛片| 欧美日韩亚洲综合一区二区三区_| 777久久人妻少妇嫩草av网站| 男人操女人黄网站| videos熟女内射| 国产成人啪精品午夜网站| 亚洲熟妇熟女久久| 久久精品亚洲精品国产色婷小说| 电影成人av| 纯流量卡能插随身wifi吗| 狂野欧美激情性xxxx|