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

    基于伴隨敏感性的風(fēng)廓線雷達(dá)和地基微波輻射計觀測對模式預(yù)報的影響評估研究

    2022-08-01 23:29:28唐兆康鮑艷松顧英杰范水勇齊亞杰崔偉陳強
    大氣科學(xué) 2022年4期
    關(guān)鍵詞:貢獻(xiàn)站點觀測

    唐兆康 鮑艷松 顧英杰 范水勇 齊亞杰 崔偉 陳強

    1 南京信息工程大學(xué)氣象災(zāi)害預(yù)報預(yù)警與評估協(xié)同創(chuàng)新中心, 南京 210044

    2 南京信息工程大學(xué)中國氣象局氣溶膠—云—降水重點開放實驗室, 南京 210044

    3 北京城市氣象研究院, 北京 100089

    4 上海衛(wèi)星工程研究所, 上海 201109

    5 南京南瑞水利水電科技有限公司,南京 211000

    1 引言

    數(shù)值天氣預(yù)報是已知一個現(xiàn)時大氣狀態(tài)的估計,用數(shù)值計算來模擬其未來一定時段的大氣運動狀態(tài)和天氣現(xiàn)象。數(shù)值模式預(yù)報的誤差主要來自初始條件的誤差和模式本身的誤差,Rabier et al. (1996)的敏感性試驗表明,許多較大的預(yù)報誤差可能是由分析中的缺陷造成的。因此,獲取更優(yōu)的初始條件對于數(shù)值模式精度的提高極其重要。

    基于數(shù)據(jù)同化理論,綜合考慮觀測值和背景值各自的誤差,最終可以獲得一個更優(yōu)的初始條件(分析場)(成巍等, 2012; 陸續(xù)等, 2015; 李佳等, 2017;王丹等, 2019)。同化大量觀測資料雖然可以改進(jìn)分析和預(yù)報結(jié)果,但是不同觀測對預(yù)報的影響必然存在著差異,并非所有被同化到分析中的觀測值都會提高預(yù)報性能。因此,有必要監(jiān)測和定量評估觀測在數(shù)據(jù)同化和預(yù)報系統(tǒng)中的應(yīng)用情況。

    觀測系統(tǒng)實驗(Observing System Experiments,簡稱OSEs)是評估觀測對模式預(yù)報影響的傳統(tǒng)方法,該方法是通過將控制分析(同化所有觀測)的預(yù)報結(jié)果和其他分析(增加或減去特定的觀測)的預(yù)報結(jié)果進(jìn)行比較,進(jìn)而評估增加或減去的觀測對預(yù)報的影響(Jung et al., 2010, 2012; Yang et al., 2014; 張斌等, 2014; 丁錦鋒, 2015; 王業(yè)桂等, 2018)。因此,如果要全面評估所有觀測對預(yù)報的影響,需要進(jìn)行大量試驗,將耗費大量的計算資源,代價昂貴。

    除OSEs 方法外,基于伴隨的預(yù)報對觀測的敏感性(Forecast Sensitivity to Observation,簡稱FSO)方法和集合預(yù)報敏感性方法也是評估觀測對預(yù)報影響的有效方法。前者最早由Baker and Daley(2000)提出,隨后Langland and Baker(2004)證明了利用伴隨來估計所有同化的觀測對短期預(yù)報誤差的影響是一種靈活且高效率的方法。Errico(2007)、Gelaro et al.(2007)和Trémolet(2007)等人討論了衡量預(yù)報誤差的高階近似及其在基于伴隨計算觀測影響背景下的特性,擴展了該方法。Auligné(2010)在NCAR(National Center for Atmospheric Research)開發(fā)了WRFDA(Weather Research and Forecasting model’s Data Assimilation)框架下基于伴隨的工具WRFDA-FSO,并用于檢驗東亞地區(qū)熱帶氣旋季節(jié)的觀測對預(yù)報的影響(Jung et al.,2013)。后者由Liu and Kalnay(2008)基于卡爾曼濾波提出,并證明集合敏感性法估計的觀測對預(yù)報的影響與基于伴隨的FSO 方法的結(jié)果相似。隨后Kunii et al.(2012)和Wu et al.(2015)利用該方法分別評估了探空觀測和AMVs(Atmospheric Motion Vectors)觀測對臺風(fēng)預(yù)報的影響,其結(jié)果表明,探空觀測對12 小時臺風(fēng)預(yù)報的正影響最大,而中層和上層的AMVs觀測對于改善臺風(fēng)的初始位置、強度和三維風(fēng)結(jié)構(gòu)及其預(yù)報尤為重要。而本文主要利用基于伴隨的FSO 方法開展觀測對預(yù)報的影響評估研究。

    與OSEs 方法相比,基于伴隨的方法在評估觀測對短期預(yù)報的影響時,結(jié)果與之相似(Gelaro and Zhu, 2009; Jung et al., 2013)。但其優(yōu)勢在于,利用該方法可以直接評估被同化到分析中的任何觀測對選定的短期預(yù)報誤差的影響,提供有關(guān)觀測對預(yù)報影響的更多細(xì)節(jié)。雖然該方法存在伴隨切線性近似的限制,但是由于其計算時間少,且可用于診斷具體目標(biāo)觀測的有效性和全面評估短期天氣預(yù)報誤差對觀測的敏感性,從而幫助設(shè)計數(shù)據(jù)同化中的觀測選擇,以及為未來開發(fā)更優(yōu)化的觀測系統(tǒng)提供參考,潛力巨大,目前已在國外眾多數(shù)值天氣預(yù)報中心里作為OSEs 的替代或補充(Zhu and Gelaro, 2008;Gelaro et al., 2010; Kim and Kim, 2014)。

    目前,國內(nèi)對于基于伴隨評估觀測對模式預(yù)報的影響的研究還處于探索和實踐階段,研究較少。其中,王曼等(2015),韓峰等(2018)利用NCAR開發(fā)的WRFDA-FSO 系統(tǒng),分析了地面和高空觀測對WRF模式預(yù)報誤差的貢獻(xiàn)。同時,已有的研究主要針對觀測對水平分辨率為幾十至上百公里模式預(yù)報的影響,缺少觀測對高分辨率模式預(yù)報的影響的研究,以及缺少定量評估中國地基微波輻射計和風(fēng)廓線雷達(dá)觀測對模式預(yù)報的影響的研究。本文通過構(gòu)建WRFDA-FSO 系統(tǒng),同化北京地區(qū)地基微波輻射計和風(fēng)廓線雷達(dá)觀測及常規(guī)觀測資料,評估觀測對水平分辨率為5 km 的WRF模式12 h 預(yù)報的影響,并分析觀測對預(yù)報貢獻(xiàn)存在差異的可能性原因,以期為北京地區(qū)觀測資料的使用及未來的觀測工作提供參考。

    2 方法與原理

    2.1 模式基本概念

    在數(shù)值天氣預(yù)報模式中,預(yù)報的大氣狀態(tài)由一個非線性模型導(dǎo)出,可以表示為

    其中,M是模型在時間段0 ≤t≤f的非線性傳播算子,xf是t=f時刻模式預(yù)報的大氣狀態(tài)向量,x0是t= 0時刻初始大氣狀態(tài)向量。初始時刻給定擾動 δx0,當(dāng)t=f時線性演變?yōu)?δxf,

    其中,M是從初始x0沿預(yù)測軌跡的切線性傳播算子。

    在三維變分資料同化中,分析增量xa?xb用 最佳線性無偏估計方程表示:

    2.2 預(yù)報誤差定義

    為了研究觀測對預(yù)報的影響,以t=f時刻的大氣參考狀態(tài)(假設(shè)為大氣真實狀態(tài))為衡量標(biāo)準(zhǔn),預(yù)報誤差定義為

    其中, 〈, 〉表示兩個向量的歐式內(nèi)積結(jié)果,C是具有預(yù)報誤差分量加權(quán)系數(shù)的對角矩陣。通常,C使用Rabier et al.(1996)的干總能量標(biāo)準(zhǔn)。g、、、和cs分別為重力加速度(單位:m s?2)、Brunt–V?is?l?頻率(單位:s?1)、位溫(單位:K)、大氣密度(單位:kg m?3)和聲速(單位:m s?1),u′和v′、θ′、p′分別表示風(fēng)變量(單位:m s?1)、溫度變量(單位:K)、氣壓變量(單位:Pa)的擾動,則預(yù)報誤差的單位為J kg?1,從而統(tǒng)一以干能量范數(shù)的標(biāo)準(zhǔn)衡量預(yù)報值與參考值的誤差。

    以背景場(xb) 和分析場(xa)為初始條件的預(yù)報場分別為xg和xf,給定大氣參考狀態(tài)xt,則預(yù)報誤差分別為

    2.3 觀測對預(yù)報誤差的影響

    為了衡量觀測對預(yù)報誤差的影響,以ea和eb之間的差定義標(biāo)量函數(shù) ?e:

    為了表示由于初始場擾動 δx0而導(dǎo)致預(yù)報誤差變化 ?e,Errico(2007)利用泰勒級數(shù)獲得了公式(7)近似:

    結(jié)合公式(2)、(3)和(4),可得 δx0=xa?xb,?e/?x0=(?xf/?x0)(?e/?xf)=2MTC(xf?),則公式(8)的一級、二級和三級近似分別為

    Gelaro et al.(2007)對三級近似進(jìn)行變形,提出了準(zhǔn)確性更高的三階近似增強形式:

    公式(12)形式上可以表示為 δe=〈δy,?e/?y〉,其中δy是由WRFDA 三維變分同化系統(tǒng)進(jìn)行計算的,預(yù)報誤差對觀測的敏感性?e/?y則是通過WRF模式及其伴隨模式和WRFDA 三維變分同化的伴隨模式進(jìn)行計算的。本文采用公式(12)衡量觀測對預(yù)報誤差的影響,計算觀測對預(yù)報的貢獻(xiàn)。而當(dāng) δe為負(fù)值時,表明同化的觀測減小了預(yù)報誤差,觀測對預(yù)報起正貢獻(xiàn);反之, δe為正值則表示同化的觀測增大了預(yù)報誤差,觀測對預(yù)報起負(fù)貢獻(xiàn)。

    3 數(shù)據(jù)與模式

    3.1 觀測資料

    同化的觀測資料來自地面氣象站(Synop)、無線電探空儀(Sound)、靜止衛(wèi)星大氣運動風(fēng)矢量(Geoamv)、風(fēng)廓線雷達(dá)(Wind Profile Radar,簡稱WPR)和地基微波輻射計(Microwave Radiometer,簡稱MWR)觀測。其中,Synop、Sound 和Geoamv觀測來自NCEP(National Centers for Environmental Prediction)集成的全球觀測資料,數(shù)據(jù)格式為PREPBUFR。WPR 和MWR觀測來自國家重點研發(fā)計劃項目課題“超大城市觀測試驗數(shù)據(jù)融合、評估和應(yīng)用示范”在北京地區(qū)的觀測,其中WPR 為L 波段的邊界層風(fēng)廓線雷達(dá),探測高度0~10 km。其產(chǎn)品數(shù)據(jù)包括:實時的采樣高度上的產(chǎn)品數(shù)據(jù),半小時平均的采樣高度上的產(chǎn)品數(shù)據(jù)和1 小時平均的采樣高度上的產(chǎn)品數(shù)據(jù);MWR 為Airda-HTG3型地基多通道微波輻射計,由輻射計14 通道微波接收機測得大氣輻射亮溫,分別為7 個溫度接收通道(51.26~58.00 GHz)和7 個濕度接收通道(22.24~31.40 GHz),探測高度0~10 km。輸出的探測數(shù)據(jù)有三級:LV0、LV1 和LV2。

    本文選取2019年9月00時、12時(協(xié)調(diào)世界時,下同)的北京地區(qū)7 個站點的WPR 的1 小時平均數(shù)據(jù)和MWR 的LV2 數(shù)據(jù),同化前依次對其進(jìn)行以下處理:(1)使用“三倍標(biāo)準(zhǔn)差準(zhǔn)則”進(jìn)行簡單質(zhì)量控制;(2)將每個觀測站點同一高度的觀測作為一個數(shù)據(jù)集,依次統(tǒng)計背景值同相應(yīng)的觀測的線性擬合系數(shù),然后利用該系數(shù)對相應(yīng)的觀測進(jìn)行偏差訂正,使得觀測與背景無偏,且觀測與背景的差值分布滿足高斯分布,進(jìn)而計算每個觀測站點的觀測誤差廓線;(3)根據(jù)背景場的模式層高度,選取與其最接近的觀測,進(jìn)行廓線觀測的垂直稀疏化;(4)將處理后的WPR 的U、V風(fēng)場觀測和MWR 的溫度T、比濕Q觀測數(shù)據(jù)轉(zhuǎn)換成同化系統(tǒng)可識別的PREPBUFR 格式。

    3.2 模式設(shè)置

    本文依托南京信息工程大學(xué)高性能計算機群構(gòu)建WRFDA-FSO 系統(tǒng)(Auligné, 2010),該系統(tǒng)主要由WRFDA_V3.8.1、WRFV3.8.1 及其伴隨模式WRFPLUS_V3.8.1 構(gòu)成。模式設(shè)置的試驗區(qū)域中心經(jīng)緯度為40°N、116°E,網(wǎng)格點數(shù)為121×151,水平分辨率為5 km,基本覆蓋整個京津冀地區(qū)。模式垂直分層從地面到模式頂高50 hPa 共51 層。模式采用的物理參數(shù)化方案如下:微物理方案為Thompson(Thompson et al., 2008);邊界層方案為ACM2(Pleim, 2007);陸面過程方案為Noah(Chen and Dudhia, 2001);長、短 波 輻 射 方 案 為RRTMG(Iacono et al., 2008);積云參數(shù)化方案為Kain-Fritsch(Kain,2004)。

    觀測對預(yù)報的影響試驗的時間為2019年9月1~30日,采用WRFDA 3D-Var 作為同化系統(tǒng),同化時刻為每天的00時和12時。模式的背景場由0.25°×0.25°的ERA5 再分析資料積分6 h 得到,分析場由背景場同化PREPBUFR 格式的觀測資料得到。背景場和分析場分別積分12 h 后,得到兩個預(yù)報場。使用ERA5 再分析資料作為預(yù)報場時刻的大氣參考狀態(tài),通過WRFDA-FSO 系統(tǒng)計算觀測對WRF模式12 h 預(yù)報的貢獻(xiàn)。同時,考慮到模式頂層附近伴隨結(jié)果可能存在較大誤差(Lorenc and Marriott, 2014; Kim et al., 2017),以及同化的觀測基本上位于10 km 以下,對模式高層的預(yù)報結(jié)果影響較小,所以選取地面至模式30 層的數(shù)據(jù)計算觀測對預(yù)報的貢獻(xiàn)。

    4 觀測對預(yù)報的貢獻(xiàn)評估

    4.1 FSO 線性近似準(zhǔn)確性分析

    由于基于伴隨的預(yù)報對觀測的敏感性方法的計算涉及模式伴隨以及泰勒級數(shù)近似,所以在評估觀測對預(yù)報的貢獻(xiàn)之前,需要驗證近似結(jié)果的準(zhǔn)確性,即公式(7)的結(jié)果與公式(12)的結(jié)果相比,符號是否一致及數(shù)值是否差距懸殊,從而保證結(jié)果和分析的可靠性。圖1 分別為2019年9月00時和12時觀測對預(yù)報誤差的影響及其近似估計結(jié)果。

    圖1 2019年9月(a)00時(協(xié)調(diào)世界時,下同)和(b)12時觀測對預(yù)報的影響(? e, 單位:103 J kg?1)及其近似估計(δ e,單位:103 J kg?1)的時間序列Fig. 1 Time series of the impact of observation on the forecast (? e, units: 103 J kg?1) and its approximate estimation (δ e, units: 103 J kg?1) at (a) 0000 UTC and (b) 1200 UTC in September 2019

    從上圖可以看出,除個別時次外,多數(shù)時次的?e均是負(fù)值,說明同化的觀測在多數(shù)時次均是改進(jìn)了WRF模式12 h 預(yù)報,減小了預(yù)報誤差,對預(yù)報為正貢獻(xiàn)。同時, δe和 ?e相比,雖然多數(shù)時次的近似結(jié)果符號一致,但仍有少數(shù)時次的近似結(jié)果并不理想,具體原因還在探索。其中,00時有6 個時次的結(jié)果近似不準(zhǔn)確,包括符號不一致的9月6、7、11、25 和26日以及線性近似明顯偏大的9月16日,而12時有10 個時次的結(jié)果近似不準(zhǔn)確,包括符號不一致的9月4、13、15、20、22~26日以及線性近似明顯偏大的9月9日。而各觀測資料對預(yù)報誤差貢獻(xiàn)的全部時次統(tǒng)計結(jié)果和剔除近似結(jié)果較差(符號不一致和數(shù)值差距懸殊)的時次統(tǒng)計結(jié)果的對比如表1 所示。

    由表1 可知,00時,剔除前的結(jié)果相比于剔除后的結(jié)果,觀測對預(yù)報的正貢獻(xiàn)作用有所減弱,這是由于00時剔除的時次的近似結(jié)果與實際結(jié)果符號相反,且多為正值。所以,若全部時次結(jié)果參與統(tǒng)計,那么必然削弱觀測對預(yù)報的正貢獻(xiàn)。而對于12時,剔除前的結(jié)果相比于剔除后的結(jié)果,觀測對預(yù)報的正貢獻(xiàn)是增強的,且WPR觀測尤其突出,這是由于9月9日12時的近似結(jié)果顯著偏高,虛假地增強了觀測對預(yù)報的正貢獻(xiàn)。若只剔除9月9日12時一個時次的結(jié)果,則WPR觀測對預(yù)報的貢獻(xiàn)將為?17.78×104J kg?1。而MWR觀測之所以對預(yù)報為負(fù)貢獻(xiàn)(11.89×104J kg?1),是由于9月9日12時并沒有MWR觀測,且其他剔除的錯誤時次的結(jié)果多為正值,從而削弱了MWR觀測對預(yù)報的正貢獻(xiàn),最終使得其對預(yù)報為虛假的負(fù)貢獻(xiàn)。

    表1 剔除近似結(jié)果不準(zhǔn)確的時次前后各觀測資料對預(yù)報誤差的貢獻(xiàn)統(tǒng)計(單位:104 J kg?1)Table 1 Statistics of contributing various observations to the forecast error before and after eliminating inaccurate approximate results (units: 104 J kg?1)

    因此,為避免近似不準(zhǔn)確時次的結(jié)果影響觀測對預(yù)報真實的貢獻(xiàn),本文對近似不準(zhǔn)確的時次結(jié)果不進(jìn)行統(tǒng)計分析。這樣雖然使得參與統(tǒng)計的試驗結(jié)果減少了,但是保證了試驗結(jié)果的可信度,且保留了至少20 個時次的試驗結(jié)果進(jìn)行累計統(tǒng)計,足以合理地評估觀測對預(yù)報的影響,統(tǒng)計結(jié)果如下所示。

    4.2 不同觀測對預(yù)報的貢獻(xiàn)分析

    圖2 為2019年9月00時和12時同化的所有觀測資料的站點觀測對預(yù)報的貢獻(xiàn)情況及站點位置分布。由圖2 可知,總體上,同化的5 種觀測資料的多數(shù)站點觀測減小了預(yù)報誤差,對預(yù)報的改進(jìn)效果為正貢獻(xiàn)。其中,12時所有WPR 的站點觀測對預(yù)報均為正貢獻(xiàn),優(yōu)于00時的WPR站點觀測。對于Geoamv觀測,00時的Geoamv觀測正貢獻(xiàn)站點數(shù)明顯多于12時的Geoamv觀測正貢獻(xiàn)站點,這可能是由于00時的Geoamv觀測質(zhì)量優(yōu)于12時的Geoamv觀測所致(劉志明等, 2002)。而MWR、Sound 和Synop觀測的正貢獻(xiàn)站點數(shù)在00時和12時差別不大。

    圖2 2019年9月(a)00時和(b)12時不同觀測的站點觀測對預(yù)報的貢獻(xiàn)及其位置分布(紅色表示該站點觀測對預(yù)報負(fù)貢獻(xiàn),藍(lán)色表示該站點觀測對預(yù)報正貢獻(xiàn);圖例中“/”前后的數(shù)字分別表示該觀測的正貢獻(xiàn)站點數(shù)和總站點數(shù))Fig. 2 Contribution of the station observation with various observations to the forecast and its location distribution at (a) 0000 UTC and (b) 1200 UTC in September 2019 (red indicates a negative contribution of the station observation to the forecast, and the blue indicates a positive contribution of the station observation to the forecast; the numbers before and after “/” in the legend represent the number of positive contribution stations and total stations of the observation respectively)

    為了評估不同觀測資料對預(yù)報的貢獻(xiàn),統(tǒng)計分析2019年9月00時和12時各觀測資料對預(yù)報誤差的貢獻(xiàn)及其有益觀測百分比,結(jié)果分別如圖3a,b所示。其中,有益觀測百分比為觀測資料中對預(yù)報為正貢獻(xiàn)的觀測數(shù)在其總觀測數(shù)中所占百分比。

    由圖3a 可知,00時和12時的5 種觀測資料對預(yù)報誤差的貢獻(xiàn)值均為負(fù)值,說明觀測減小了預(yù)報誤差,對預(yù)報的改進(jìn)效果為正貢獻(xiàn)。其中,00時的MWR觀測對預(yù)報的影響最大,其次分別為Sound、Synop、WPR 和Geoamv觀測。而12時依然是MWR觀測對預(yù)報的影響最大,其次分別為Synop、WPR 、Sound 和Geoamv觀測。而由圖3b可知,除12時的Geoamv觀測外,其余觀測的有益觀測百分比均超過了50%,與Lorenc and Marriott(2014)和Kim et al.(2017)的研究結(jié)果基本一致。其中,Synop 的有益觀測百分比較高,超過了65%,WPR、MWR 和Sound 的有益觀測百分比處于55%~60%。

    圖3 2019年9月00時和12時的(a)不同觀測對預(yù)報誤差的貢獻(xiàn)(單位:104 J kg?1)和(b)不同觀測的有益觀測百分比Fig. 3 (a) Contribution (units: 104 J kg?1) of various observations to the forecast error and (b) percentage of beneficial observations of various observations at 0000 UTC and 1200 UTC in September 2019

    4.3 WPR 和MWR觀測對預(yù)報的貢獻(xiàn)分析

    本文統(tǒng)計2019年9月00時和12時觀測對預(yù)報的影響試驗結(jié)果中WPR 和MWR觀測對預(yù)報誤差的貢獻(xiàn),分析其不同觀測要素、不同站點和不同高度層的觀測對WRF模式12 h 預(yù)報的影響。

    4.3.1 不同觀測要素觀測對預(yù)報的貢獻(xiàn)分析

    為了評估WPR 和MWR 各觀測要素觀測對預(yù)報的影響,選取Sound觀測作參考。同時,考慮到觀測資料的觀測數(shù)存在較大差異,本文統(tǒng)計分析各觀測要素觀測對預(yù)報誤差的平均貢獻(xiàn)(總貢獻(xiàn)除以總觀測數(shù)),結(jié)果如表2 所示。

    由表2 可知,除00時的WPR 的U觀測,WPR 和MWR 各觀測要素觀測對預(yù)報誤差的平均貢獻(xiàn)均為負(fù)值,表明各觀測要素觀測減小了預(yù)報誤差,對預(yù)報為正貢獻(xiàn)。同時,對比各觀測要素的平均貢獻(xiàn)值的大小,可知,WPR 的V觀測對預(yù)報的改進(jìn)效果優(yōu)于U觀測,而MWR觀測對預(yù)報的正貢獻(xiàn)主要來自T觀測。同時,相比于Sound 的觀測,WPR 的風(fēng)場觀測對預(yù)報的改進(jìn)效果總體上優(yōu)于Sound 的風(fēng)場觀測,MWR 的T觀測對預(yù)報的改進(jìn)效果優(yōu)于Sound的T觀測,而MWR 和Sound 的Q觀測對預(yù)報的平均貢獻(xiàn)值相當(dāng)。由此可見,WPR 和MWR觀測對預(yù)報的改進(jìn)效果顯著。在目前Sound觀測數(shù)據(jù)較少的情況下,WPR 和MWR觀測數(shù)據(jù)對于預(yù)報的改進(jìn)有著較大的意義。

    表2 2019年9月00時和12時廓線觀測對預(yù)報誤差的平均貢獻(xiàn)(單位:J kg?1)統(tǒng)計Table 2 Statistics of the average contribution of profile observations to the forecast error (units: J kg?1) at 0000 UTC and 1200 UTC in September 2019

    4.3.2 不同觀測站點觀測對預(yù)報的貢獻(xiàn)分析

    圖4 為00時和12時的WPR 和MWR 各站點觀測對預(yù)報誤差的平均貢獻(xiàn)情況,7 個站點分別是海淀站(HD),延慶站(YQ),懷柔站(HR),密云站(MY),平谷站(PG),大興站(DX)和霞云嶺站(XYL)。由圖4 可知,WPR 的站點觀測中,在00時,海淀站、密云站、平谷站和大興站4 個站點的WPR觀測對預(yù)報均為正貢獻(xiàn),而剩余3 個站點的觀測則增大了預(yù)報誤差,其中延慶站負(fù)貢獻(xiàn)較大。而在12時,7 個站點的WPR觀測均減小了預(yù)報誤差,對預(yù)報為正貢獻(xiàn)??傮w上,7 個站點中,大興站的WPR觀測對預(yù)報的改進(jìn)效果最優(yōu)。

    圖4 2019年9月00時和12時北京地區(qū)7站點(HD、YQ、HR、MY、PG、DX 和XYL)的(a)WPR觀測和(b)MWR觀測對預(yù)報誤差的平均貢獻(xiàn)(單位:102 J kg?1)。HD:海淀,YQ:延慶,HR:懷柔,MY:密云,PG:平谷,DX:大興,XYL:霞云嶺Fig. 4 Average contribution (units: 102 J kg?1) of (a) WPR observations and (b) MWR observations to the forecast error at seven stations, which are Haidian (HD), Yanqing (YQ), Huairou (HR), Miyun (MY), Pinggu (PG), Daxing (DX), and Xiayunling (XYL) in Beijing at 0000 UTC and 1200 UTC in September 2019

    對于MWR觀測,除12時的海淀站觀測對預(yù)報產(chǎn)生了較低的負(fù)貢獻(xiàn),其余站點的MWR觀測對預(yù)報的貢獻(xiàn)總體上均為正貢獻(xiàn),其中大興站和霞云嶺站的MWR觀測對預(yù)報的改進(jìn)效果較優(yōu)。

    4.3.3 不同高度層觀測對預(yù)報的貢獻(xiàn)分析

    為了評估不同高度層WPR 和MWR觀測對預(yù)報的貢獻(xiàn),分別統(tǒng)計2019年9月00時和12時的WPR 和MWR 各個高度層觀測對預(yù)報誤差的平均貢獻(xiàn),結(jié)果如表3、4 所示。由表3 可知,00時,WPR觀測對預(yù)報的負(fù)貢獻(xiàn)主要位于1000~2000 m之間,且主要來自U觀測,而其余高度層觀測對預(yù)報總體上均為正貢獻(xiàn)。而在12時,總體上WPR各個高度層觀測均改進(jìn)了預(yù)報,減小了預(yù)報誤差。

    表3 2019年9月00時和12時不同高度層WPR觀測對預(yù)報誤差的平均貢獻(xiàn)(單位:J kg?1)統(tǒng)計Table 3 Statistics of the average contribution of WPR observations at various altitudes to the forecast error (units:J kg?1) at 0000 UTC and 1200 UTC in September 2019

    由表4 可知,MWR 的多數(shù)高度層觀測均減小了預(yù)報誤差,對預(yù)報為正貢獻(xiàn)。00時,MWR觀測對預(yù)報的負(fù)貢獻(xiàn)主要位于500~800 hPa 之間,且主要來自T觀測,而其余高度層觀測對預(yù)報均為正貢獻(xiàn)。12時,MWR觀測對預(yù)報的負(fù)貢獻(xiàn)主要位于600~700 hPa,其余高度層的觀測總體上均減小了預(yù)報誤差,對預(yù)報為正貢獻(xiàn)。同時,T觀測對預(yù)報的貢獻(xiàn)主要位于近地面800 hPa 以下,且對預(yù)報為正貢獻(xiàn),600 hPa 以上的T觀測對預(yù)報的影響較小。

    表4 2019年9月00時和12時不同高度層MWR觀測對預(yù)報誤差的平均貢獻(xiàn)(單位:J kg?1)統(tǒng)計Table 4 Statistics of the average contribution of MWR observations at various altitudes to the forecast error (units:J kg?1) at 0000 UTC and 1200 UTC in September 2019

    4.4 觀測對預(yù)報貢獻(xiàn)差異的原因分析

    雖然4.2 和4.3 節(jié)分析了不同觀測資料對WRF模式12 h 預(yù)報的貢獻(xiàn)情況,但是不同觀測類型或觀測要素或觀測高度的觀測對預(yù)報的貢獻(xiàn)存在明顯的差異,尤其是2019年9月00時的試驗結(jié)果中,WPR 的U觀測對預(yù)報為負(fù)貢獻(xiàn),WPR觀測的正貢獻(xiàn)站點數(shù)也相對較少,并且MWR觀測對預(yù)報的正貢獻(xiàn)非常突出(圖3a),對預(yù)報的改進(jìn)效果遠(yuǎn)遠(yuǎn)超出了Sound觀測,且其對預(yù)報為正貢獻(xiàn)的觀測主要是集中在近地面800 hPa 以下的T觀測。所以選取2019年9月00時觀測對預(yù)報的影響試驗結(jié)果中WPR 和MWR觀測相應(yīng)的數(shù)據(jù),進(jìn)行分類統(tǒng)計分析,探究可能導(dǎo)致觀測對預(yù)報的貢獻(xiàn)存在差異的原因。

    4.4.1 WPR觀測對預(yù)報貢獻(xiàn)差異性分析

    針對00時的WPR 的U、V觀測數(shù)據(jù),選取新息增量(觀測與背景的差異)的絕對值|δy|和觀測誤差σ兩個參數(shù),分別以4 m s?1和3 m s?1為閾值對結(jié)果進(jìn)行分類,統(tǒng)計各類區(qū)間的正負(fù)貢獻(xiàn)值和正負(fù)貢獻(xiàn)觀測數(shù),結(jié)果如下表5 所示。其中,觀測對預(yù)報誤差的貢獻(xiàn)值為正值時,表示觀測增大了預(yù)報誤差,對預(yù)報為負(fù)貢獻(xiàn),而貢獻(xiàn)值為負(fù)值則表示觀測減小了預(yù)報誤差,對預(yù)報為正貢獻(xiàn)。

    由表5 可知,當(dāng)觀測誤差σ>3時,WPR觀測數(shù)較少,僅約占總觀測數(shù)的7.04%,其對預(yù)報的貢獻(xiàn)值(正、負(fù)貢獻(xiàn))也較低;而當(dāng)觀測誤差σ≤3時,對應(yīng)的觀測數(shù)和觀測對預(yù)報的貢獻(xiàn)值(正、負(fù)貢獻(xiàn))均較高,表明2019年9月00時的WPR觀測中,多數(shù)觀測數(shù)據(jù)精度較高,且對預(yù)報產(chǎn)生了較大影響,只有少數(shù)觀測數(shù)據(jù)精度較差,其對預(yù)報的影響也較低。

    表5 2019年9月00時WPR觀測對預(yù)報誤差的貢獻(xiàn)(單位:J kg?1)分類統(tǒng)計Table 5 Classified statistics of the contribution of WPR observations to the forecast error (units: J kg?1) at 0000 UTC in September 2019

    當(dāng)觀測誤差σ≤3 且|δy|>4時,雖然觀測數(shù)較少,但從平均貢獻(xiàn)上來看,其對預(yù)報的正、負(fù)貢獻(xiàn)值偏高,說明σ≤3 且|δy|>4 的觀測對預(yù)報產(chǎn)生了較大的影響,這與其計算方法有直接的聯(lián)系。綜合公式(3)和公式(12)可以看出,觀測對預(yù)報的貢獻(xiàn)一定程度上與δy成正比,而與觀測誤差σ成反比。同時,針對負(fù)貢獻(xiàn)的觀測數(shù)較少卻產(chǎn)生了更高的負(fù)貢獻(xiàn)值的現(xiàn)象,本文詳細(xì)地查看了觀測誤差σ≤3 且|δy|>4 的觀測點,結(jié)果如下圖5 所示,圖中“YQ_10”標(biāo)記表示延慶站9月10日觀測。

    圖5 2019年9月00時觀測誤差σ≤3 且新息增量|δy|>4 的WPR(a)U觀測和(b)V觀測對應(yīng)的對預(yù)報誤差的貢獻(xiàn)(單位:102 J kg?1)與新息增量(單位:m s?1)的散點圖(黑色表示該觀測對預(yù)報負(fù)貢獻(xiàn),灰色表示該觀測對預(yù)報正貢獻(xiàn))Fig. 5 Scatter plot of forecast error contribution (units: 102 J kg?1) and innovation increment (units: m s?1) corresponding to (a) U observations and(b) V observations of WPR with observation error less than or equal to 3 and innovation increment greater than 4 (σ≤ 3 and |δy|>4) at 0000 UTC in September 2019 (black indicates a negative contribution of the observation to the forecast, and gray indicates a positive contribution of the observation to the forecast)

    由圖易知,負(fù)貢獻(xiàn)的觀測數(shù)較少卻產(chǎn)生了更高的負(fù)貢獻(xiàn)值,主要是由于延慶站9月10日個別觀測產(chǎn)生了極高的負(fù)貢獻(xiàn)所致。9月10日延慶站U、V觀測對預(yù)報的負(fù)貢獻(xiàn)值,分別約占該類觀測對預(yù)報總的負(fù)貢獻(xiàn)的49.01%和55.36%,說明少數(shù)個別結(jié)果對總體的結(jié)果產(chǎn)生了極大的影響。而且該部分觀測通過了雙權(quán)重法(Lanzante, 1996)的離群值檢驗,說明可能并不是觀測資料導(dǎo)致了觀測對預(yù)報為負(fù)貢獻(xiàn)。而在觀測質(zhì)量較高(σ≤3),背景場質(zhì)量較低(|δy|>4)的情況下,該少部分觀測對預(yù)報的貢獻(xiàn)卻為負(fù)貢獻(xiàn),也說明觀測對預(yù)報為負(fù)貢獻(xiàn)的原因并不是觀測和背景場的問題,而是可能與同化系統(tǒng)和模式伴隨的局限性有關(guān),比如本文背景、觀測誤差協(xié)方差的計算和同化方案中一些參數(shù)的設(shè)計并不完美,無法適用于所有時次(Lorenc and Marriott, 2014),從而使得一些觀測質(zhì)量較好的觀測對預(yù)報為負(fù)貢獻(xiàn)。這也說明利用基于伴隨的預(yù)報對觀測的敏感性方法評估觀測對預(yù)報的影響時,必須使用長時間的觀測數(shù)據(jù)集,本文試驗所用數(shù)據(jù)還是較少,未來需要進(jìn)一步改進(jìn)。

    而當(dāng)σ≤3 且|δy|≤4時,該類觀測數(shù)較多,且超過50%的觀測對預(yù)報的貢獻(xiàn)為正貢獻(xiàn),但也存在較多的觀測對預(yù)報產(chǎn)生了負(fù)貢獻(xiàn),U、V觀測中對預(yù)報為負(fù)貢獻(xiàn)的觀測數(shù)分別為730 和697,分別約占觀測總數(shù)的37.49%,35.80%。這可能仍是本文背景、觀測誤差協(xié)方差的計算和同化方案中一些參數(shù)的設(shè)計并不完美,導(dǎo)致其無法適用于所有時次。但是具體的原因,需要未來通過更多的試驗進(jìn)行更深入地研究與分析。

    4.4.2 MWR觀測對預(yù)報貢獻(xiàn)差異性分析

    針對00時的MWR 的T、Q觀測數(shù)據(jù),選取新息增量的絕對值|δy|和觀測誤差σ,分別以觀測誤差的1 K 和1 g kg?1以及|δy|的2 K 和2 g kg?1為閾值對結(jié)果進(jìn)行分類,統(tǒng)計各類區(qū)間的正負(fù)貢獻(xiàn)值和正負(fù)貢獻(xiàn)觀測數(shù),結(jié)果如下表6 所示。其中,貢獻(xiàn)值為正值表示觀測增大了預(yù)報誤差,對預(yù)報為負(fù)貢獻(xiàn),而貢獻(xiàn)值為負(fù)值則表示觀測減小了預(yù)報誤差,對預(yù)報為正貢獻(xiàn)。

    由表6 可知,當(dāng)觀測誤差σ≤1時,T觀測數(shù)僅為562,約占T觀測總數(shù)的23.57%,但是該部分觀測對模式12 h 預(yù)報誤差的貢獻(xiàn)值(正、負(fù)貢獻(xiàn)絕對值)卻極高,達(dá)到了T觀測總貢獻(xiàn)值(正、負(fù)貢獻(xiàn)絕對值)的73.34%。其中,相比于新息增量絕對值|δy|≤2 的觀測,|δy|>2 的觀測對預(yù)報誤差的貢獻(xiàn)值總體上更高。這表明T觀測中觀測誤差較小且與背景差異較大的觀測,對預(yù)報產(chǎn)生了極大影響,這與WPR觀測的分析結(jié)果一致。

    表6 2019年9月00時 MWR觀測對預(yù)報誤差的貢獻(xiàn)(單位:J kg?1)分類統(tǒng)計Table 6 Classified statistics of the contribution of MWR observations to the forecast error (units: J kg?1) at 0000 UTC in September 2019

    而Q觀測的統(tǒng)計結(jié)果依然體現(xiàn)了觀測誤差較小且與背景差異較大的觀測,對預(yù)報的影響更大。雖然對于觀測誤差σ≤1 的觀測數(shù)為1274,約占Q觀測總數(shù)的55.29%,而對12 h 預(yù)報誤差的貢獻(xiàn)值(正、負(fù)貢獻(xiàn)絕對值)僅達(dá)到了總貢獻(xiàn)值(正、負(fù)貢獻(xiàn)絕對值)的54.23%。并沒有像U、V和T觀測那樣,突出質(zhì)量較高的觀測對預(yù)報有著更大的影響。這主要是由于Q觀測中觀測誤差σ≤1 的觀測約89.64%的觀測來自700 hPa 以上的觀測,而700 hPa 以上的Q觀測數(shù)值非常小,導(dǎo)致統(tǒng)計得到的觀測誤差也非常小,所以該部分觀測的觀測誤差σ≤1 并不能代表其觀測質(zhì)量較高。

    與表5 中WPR觀測的結(jié)果一樣的是,MWR觀測中同樣存在少量的觀測誤差σ≤1 但對預(yù)報為較高的負(fù)貢獻(xiàn)的觀測,而且該部分T、Q觀測同樣通過了雙權(quán)重法的離群值檢驗。所以,這可能仍是本文背景、觀測誤差協(xié)方差的計算和同化方案中一些參數(shù)的設(shè)計并不完美,導(dǎo)致其無法適用于所有時次。

    同時,考慮到2019年9月00時的T觀測對預(yù)報的正貢獻(xiàn)主要位于800 hPa 以下(表4),所以統(tǒng)計T觀測對預(yù)報誤差的貢獻(xiàn)值以及分析增量隨高度的分布,以探究造成T觀測對預(yù)報的正貢獻(xiàn)較高的可能因素,結(jié)果如圖6a,b 所示。而由公式(3)可知,分析增量又受新息增量和觀測誤差的影響,分析增量一定程度上可能與新息增量成正比,而與觀測誤差成反比,所以同樣統(tǒng)計其隨高度的分布,探究其對分析增量的影響,結(jié)果如圖6c,d 所示。

    由圖6a,b 可知,T觀測對預(yù)報誤差的貢獻(xiàn)值和分析增量的高值區(qū)均主要位于800 hPa 以下,而800 hPa 以上的分析增量和T觀測對預(yù)報誤差的貢獻(xiàn)值均較小,這說明分析增量可能影響了T觀測對預(yù)報誤差貢獻(xiàn)程度。而對比圖6b、c 和d 可知,MWR 的T觀測的分析增量同新息增量并不是正比關(guān)系,而同觀測誤差一定程度上成反比關(guān)系,說明MWR 的T觀測的分析增量大小主要受觀測誤差的影響。低層新息增量雖然較小,但是對應(yīng)的觀測誤差也較小,反而使得新息增量的權(quán)重較大,而高層新息增量雖然較大,但是觀測誤差也較大,于是使得新息增量的權(quán)重較小,從而導(dǎo)致了分析增量低層高值,高層低值的現(xiàn)象。

    圖6 2019年9月00時MWR 的T觀測對應(yīng)的(a)對預(yù)報誤差的貢獻(xiàn)(單位:103 J kg?1)、(b)分析增量(單位:K)、(c)新息增量(單位:K)以及(d)觀測誤差(單位:K)隨高度的分布Fig. 6 The distribution of (a) a forecast error contribution (units: 103 J kg?1), (b) an analysis increment (units: K), (c) an innovation increment (units:K), and (d) an observation error (units: K) with a height corresponding to temperature observations of MWR at 0000 UTC in September 2019

    綜合以上分析,對于MWR 的T觀測,觀測對預(yù)報的貢獻(xiàn)差異更多地受觀測的觀測誤差的影響。800 hPa 以下的T觀測的觀測誤差較小,使得分析增量較高,對應(yīng)的MWR 的T觀測對預(yù)報的正貢獻(xiàn)在800 hPa 以下最為顯著;而600 hPa 以上的T觀測的觀測誤差較大,使得分析增量也較低,對應(yīng)的600 hPa 以上的T觀測對預(yù)報主要為較低的負(fù)貢獻(xiàn)。

    5 總結(jié)與討論

    隨著越來越多的觀測資料被應(yīng)用于數(shù)據(jù)同化,改進(jìn)數(shù)值模式預(yù)報,定量評估觀測在數(shù)據(jù)同化和模式預(yù)報系統(tǒng)中的應(yīng)用情況將變得尤為重要。本文通過2019年9月1日 至30日00時 和12時 的FSO試驗,驗證了本文構(gòu)建的WRFDA-FSO 系統(tǒng)可以較好地定量評估觀測對預(yù)報的貢獻(xiàn),并對試驗結(jié)果進(jìn)行了分析,得到以下結(jié)論:

    (1)同化的觀測資料(MWR、WPR、Sound、Synop 和Geoamv)能夠較好地改進(jìn)WRF模式水平分辨率5 km 的12 h 預(yù)報,5 種觀測資料的有益觀測百分比基本上均超過了50%。其中,MWR觀測對預(yù)報的改進(jìn)效果最優(yōu),WPR 的風(fēng)場觀測對預(yù)報的改進(jìn)效果也優(yōu)于Sound 的風(fēng)場。

    (2)WPR觀測和MWR觀測的7 個觀測站點中,大興站的WPR觀測和大興站以及霞云嶺站的MWR觀測對預(yù)報的改進(jìn)效果較優(yōu)。WPR 和MWR 的觀測要素中,WPR 的V觀測對預(yù)報的改進(jìn)效果優(yōu)于其U觀測,而MWR 的T觀測對預(yù)報的改進(jìn)效果優(yōu)于其Q觀測,且T觀測中對預(yù)報為正貢獻(xiàn)的觀測主要位于近地面800 hPa 以下。

    (3)2019年9月00時的WPR觀測和MWR觀測的分類統(tǒng)計分析表明,觀測數(shù)據(jù)的觀測誤差及其與背景值的偏差的不同,會使得觀測對預(yù)報的貢獻(xiàn)存在差異。其中觀測誤差較?。ㄓ^測質(zhì)量較高),且背景場質(zhì)量相對較差(|δy|較大)時,對應(yīng)的觀測對預(yù)報有著明顯的影響,易對預(yù)報產(chǎn)生較高的貢獻(xiàn)值(正、負(fù)貢獻(xiàn))。

    本文對利用基于伴隨的方法評估觀測對高分辨率模式預(yù)報的影響進(jìn)行了初步研究,對該方法以及各觀測資料對預(yù)報的具體影響有了一定的了解。關(guān)于本文出現(xiàn)的近似結(jié)果準(zhǔn)確性的問題,未來需要對其進(jìn)行更深入地研究和改進(jìn)。而觀測對預(yù)報貢獻(xiàn)存在差異的原因也是未來值得研究的問題。

    猜你喜歡
    貢獻(xiàn)站點觀測
    觀測到恒星死亡瞬間
    軍事文摘(2023年18期)2023-11-03 09:45:42
    中國共產(chǎn)黨百年偉大貢獻(xiàn)
    為加快“三個努力建成”作出人大新貢獻(xiàn)
    基于Web站點的SQL注入分析與防范
    電子制作(2019年14期)2019-08-20 05:43:42
    2017~2018年冬季西北地區(qū)某站點流感流行特征分析
    貢獻(xiàn)榜
    海洋貢獻(xiàn)2500億
    商周刊(2017年6期)2017-08-22 03:42:37
    天測與測地VLBI 測地站周圍地形觀測遮掩的討論
    首屆歐洲自行車共享站點協(xié)商會召開
    中國自行車(2017年1期)2017-04-16 02:53:52
    怕被人認(rèn)出
    故事會(2016年21期)2016-11-10 21:15:15
    日韩欧美在线二视频| 亚洲一区高清亚洲精品| 狂野欧美激情性xxxx| 欧美日韩一级在线毛片| 女警被强在线播放| 色av中文字幕| 最新美女视频免费是黄的| 亚洲精品久久国产高清桃花| 亚洲内射少妇av| 日韩欧美三级三区| 精华霜和精华液先用哪个| 一区二区三区高清视频在线| 精品一区二区三区视频在线观看免费| 一区二区三区激情视频| 久久精品亚洲精品国产色婷小说| av天堂中文字幕网| 日韩成人在线观看一区二区三区| 久久久久久久午夜电影| 午夜免费激情av| 老司机福利观看| 亚洲国产欧美网| 人妻夜夜爽99麻豆av| 亚洲美女视频黄频| 一夜夜www| 午夜激情欧美在线| 一级毛片女人18水好多| 在线a可以看的网站| 国产真实伦视频高清在线观看 | 久久久精品大字幕| 波野结衣二区三区在线 | 99国产精品一区二区蜜桃av| 免费在线观看亚洲国产| 日本三级黄在线观看| 亚洲av中文字字幕乱码综合| 国产亚洲精品一区二区www| 老司机在亚洲福利影院| 在线观看66精品国产| 久久精品人妻少妇| 九九久久精品国产亚洲av麻豆| 亚洲av日韩精品久久久久久密| 久久亚洲真实| aaaaa片日本免费| 高清在线国产一区| 精品国内亚洲2022精品成人| 变态另类丝袜制服| 国产一区二区在线av高清观看| 99国产综合亚洲精品| 亚洲av熟女| 免费看美女性在线毛片视频| 极品教师在线免费播放| 亚洲精品色激情综合| 黄色日韩在线| 国产亚洲精品综合一区在线观看| 人妻丰满熟妇av一区二区三区| 精品久久久久久久久久久久久| 尤物成人国产欧美一区二区三区| 国产伦一二天堂av在线观看| 男女午夜视频在线观看| 搡老妇女老女人老熟妇| 婷婷丁香在线五月| 欧美日韩精品网址| 国产精品久久久人人做人人爽| 亚洲精品影视一区二区三区av| 亚洲男人的天堂狠狠| 99久久九九国产精品国产免费| www日本黄色视频网| 一夜夜www| 老司机福利观看| 激情在线观看视频在线高清| 亚洲美女视频黄频| 嫁个100分男人电影在线观看| 最近最新中文字幕大全电影3| 欧美区成人在线视频| 丰满人妻一区二区三区视频av | 国产伦精品一区二区三区四那| 亚洲精品亚洲一区二区| 亚洲性夜色夜夜综合| 欧美丝袜亚洲另类 | 一级毛片女人18水好多| ponron亚洲| 欧美日韩中文字幕国产精品一区二区三区| 欧美日韩福利视频一区二区| 在线观看舔阴道视频| av中文乱码字幕在线| 免费观看人在逋| 亚洲精品456在线播放app | 9191精品国产免费久久| 日韩欧美精品v在线| 亚洲av第一区精品v没综合| 99久国产av精品| 99久久九九国产精品国产免费| 国产精品久久久人人做人人爽| 俺也久久电影网| 精品一区二区三区av网在线观看| 老熟妇仑乱视频hdxx| 久久精品综合一区二区三区| 桃色一区二区三区在线观看| 欧美精品啪啪一区二区三区| 久久精品国产亚洲av香蕉五月| 国产精品久久视频播放| 欧美一级a爱片免费观看看| 三级国产精品欧美在线观看| 亚洲无线在线观看| 亚洲欧美精品综合久久99| 国产黄片美女视频| 欧美黄色片欧美黄色片| 欧美日韩福利视频一区二区| 在线播放无遮挡| 国产午夜精品久久久久久一区二区三区 | 天堂√8在线中文| 欧美黑人欧美精品刺激| 小蜜桃在线观看免费完整版高清| 真人做人爱边吃奶动态| 国产精品一区二区三区四区免费观看 | 丁香六月欧美| 婷婷精品国产亚洲av| 免费看a级黄色片| 国产精品一区二区三区四区久久| 91久久精品国产一区二区成人 | 国产精品亚洲美女久久久| 91久久精品电影网| 一进一出好大好爽视频| 男女之事视频高清在线观看| 国产成人aa在线观看| 黑人欧美特级aaaaaa片| 99久久九九国产精品国产免费| 色综合亚洲欧美另类图片| 亚洲av免费高清在线观看| 中文字幕高清在线视频| eeuss影院久久| 欧美xxxx黑人xx丫x性爽| 国内精品久久久久精免费| 好男人电影高清在线观看| 亚洲av一区综合| 内地一区二区视频在线| 色哟哟哟哟哟哟| 舔av片在线| 亚洲av不卡在线观看| 欧美最新免费一区二区三区 | 嫩草影院精品99| aaaaa片日本免费| 久久性视频一级片| 在线观看66精品国产| 白带黄色成豆腐渣| 精品一区二区三区视频在线观看免费| 日韩欧美 国产精品| 久久这里只有精品中国| 亚洲欧美日韩高清专用| 亚洲无线观看免费| 午夜a级毛片| 成年人黄色毛片网站| 波多野结衣巨乳人妻| 男女午夜视频在线观看| 久久这里只有精品中国| 中文资源天堂在线| 成人高潮视频无遮挡免费网站| 欧美日本视频| 免费在线观看日本一区| 亚洲国产日韩欧美精品在线观看 | 叶爱在线成人免费视频播放| 精品久久久久久久久久久久久| 国产高清视频在线播放一区| 日韩欧美免费精品| 久久国产精品人妻蜜桃| 女人高潮潮喷娇喘18禁视频| 日韩高清综合在线| a级一级毛片免费在线观看| 村上凉子中文字幕在线| 中亚洲国语对白在线视频| 一级毛片女人18水好多| 久久精品国产自在天天线| 午夜免费激情av| 国产亚洲精品综合一区在线观看| 男女做爰动态图高潮gif福利片| 毛片女人毛片| 国产精品亚洲美女久久久| 他把我摸到了高潮在线观看| 国产成人aa在线观看| 国产毛片a区久久久久| 黄片小视频在线播放| 亚洲人与动物交配视频| 亚洲专区国产一区二区| 午夜福利免费观看在线| 欧美成人免费av一区二区三区| 久久精品91蜜桃| 深爱激情五月婷婷| 制服人妻中文乱码| av中文乱码字幕在线| 日本 欧美在线| 国产激情欧美一区二区| 日本 av在线| 欧美日韩福利视频一区二区| 成年版毛片免费区| 亚洲av电影在线进入| 久久伊人香网站| 久久久久久九九精品二区国产| 亚洲欧美日韩卡通动漫| 国产探花在线观看一区二区| 中文字幕人成人乱码亚洲影| 一级作爱视频免费观看| 搞女人的毛片| 天堂√8在线中文| 成人亚洲精品av一区二区| 色老头精品视频在线观看| 精品久久久久久,| 高清毛片免费观看视频网站| 制服人妻中文乱码| 亚洲 国产 在线| 亚洲五月天丁香| 亚洲中文日韩欧美视频| 久久中文看片网| 性欧美人与动物交配| 亚洲中文字幕日韩| 亚洲欧美一区二区三区黑人| 日日夜夜操网爽| 亚洲精品日韩av片在线观看 | 亚洲色图av天堂| АⅤ资源中文在线天堂| 亚洲国产精品久久男人天堂| 午夜免费成人在线视频| 桃红色精品国产亚洲av| 久久精品国产亚洲av涩爱 | 久久九九热精品免费| 两性午夜刺激爽爽歪歪视频在线观看| 国产成人福利小说| 高清在线国产一区| 丁香欧美五月| 三级毛片av免费| 又紧又爽又黄一区二区| 香蕉丝袜av| 夜夜躁狠狠躁天天躁| 神马国产精品三级电影在线观看| 18禁裸乳无遮挡免费网站照片| 在线播放无遮挡| 亚洲美女视频黄频| 久久精品91无色码中文字幕| 亚洲精品久久国产高清桃花| 亚洲无线在线观看| 俺也久久电影网| 色综合欧美亚洲国产小说| 夜夜躁狠狠躁天天躁| 黄色成人免费大全| 丁香欧美五月| 90打野战视频偷拍视频| aaaaa片日本免费| 又黄又粗又硬又大视频| 亚洲精品成人久久久久久| 亚洲aⅴ乱码一区二区在线播放| 亚洲avbb在线观看| 一级黄色大片毛片| 精品乱码久久久久久99久播| 男插女下体视频免费在线播放| 久久国产精品人妻蜜桃| 欧美一区二区精品小视频在线| 久久久久久大精品| 亚洲一区高清亚洲精品| 国产精品一区二区三区四区免费观看 | 亚洲专区中文字幕在线| 偷拍熟女少妇极品色| 99热这里只有是精品50| 亚洲国产精品sss在线观看| 在线观看免费午夜福利视频| 国产精品1区2区在线观看.| 成人特级av手机在线观看| 精品熟女少妇八av免费久了| 亚洲国产色片| 久久精品综合一区二区三区| 欧美一区二区亚洲| 在线免费观看不下载黄p国产 | 国产精品久久久久久精品电影| 男插女下体视频免费在线播放| 国产精品电影一区二区三区| 亚洲精品一区av在线观看| av视频在线观看入口| 国产精品 欧美亚洲| 久久精品91蜜桃| 国产精品日韩av在线免费观看| 亚洲成人久久爱视频| 亚洲一区二区三区不卡视频| 九九久久精品国产亚洲av麻豆| 18禁裸乳无遮挡免费网站照片| 国产在线精品亚洲第一网站| 一区二区三区高清视频在线| 麻豆国产97在线/欧美| 操出白浆在线播放| 一级毛片高清免费大全| 中文字幕av成人在线电影| 日本免费a在线| 欧美丝袜亚洲另类 | 18禁美女被吸乳视频| 国内揄拍国产精品人妻在线| 欧美zozozo另类| 国内久久婷婷六月综合欲色啪| 婷婷精品国产亚洲av| АⅤ资源中文在线天堂| 欧美乱妇无乱码| 亚洲av中文字字幕乱码综合| 老熟妇仑乱视频hdxx| 综合色av麻豆| 色综合欧美亚洲国产小说| 熟女电影av网| 九九热线精品视视频播放| 99久久综合精品五月天人人| 欧美日韩黄片免| 性色av乱码一区二区三区2| 国产在视频线在精品| 亚洲中文字幕日韩| 精品久久久久久久毛片微露脸| 成人鲁丝片一二三区免费| 婷婷亚洲欧美| 亚洲精品成人久久久久久| 欧美成人a在线观看| 老师上课跳d突然被开到最大视频 久久午夜综合久久蜜桃 | 亚洲精品乱码久久久v下载方式 | 69人妻影院| www国产在线视频色| 国产精品99久久99久久久不卡| 亚洲熟妇中文字幕五十中出| 18禁裸乳无遮挡免费网站照片| 日本一本二区三区精品| 一进一出抽搐gif免费好疼| 一个人免费在线观看的高清视频| 午夜激情欧美在线| 好看av亚洲va欧美ⅴa在| 国产三级在线视频| 99久久无色码亚洲精品果冻| 九色国产91popny在线| 午夜免费成人在线视频| 成人特级av手机在线观看| 蜜桃亚洲精品一区二区三区| 波野结衣二区三区在线 | 亚洲欧美一区二区三区黑人| www.熟女人妻精品国产| 国产亚洲精品久久久久久毛片| 欧美黄色片欧美黄色片| 久久久国产成人免费| 日韩大尺度精品在线看网址| 两个人的视频大全免费| 国产成人福利小说| 一级黄片播放器| 小蜜桃在线观看免费完整版高清| 嫩草影视91久久| 日韩av在线大香蕉| 亚洲男人的天堂狠狠| 亚洲中文字幕日韩| 亚洲精品色激情综合| 两人在一起打扑克的视频| 久久久久国产精品人妻aⅴ院| 91麻豆精品激情在线观看国产| tocl精华| 国产激情偷乱视频一区二区| 午夜福利在线在线| а√天堂www在线а√下载| 精品日产1卡2卡| 手机成人av网站| 欧美午夜高清在线| 最近在线观看免费完整版| 99久久精品一区二区三区| 波多野结衣高清无吗| 深爱激情五月婷婷| 99久久综合精品五月天人人| 婷婷精品国产亚洲av| 19禁男女啪啪无遮挡网站| 亚洲国产色片| 小蜜桃在线观看免费完整版高清| 国产乱人视频| 亚洲狠狠婷婷综合久久图片| 女人被狂操c到高潮| www日本在线高清视频| 久久久久久久久大av| 国产成人aa在线观看| 韩国av一区二区三区四区| 午夜影院日韩av| 69av精品久久久久久| 老师上课跳d突然被开到最大视频 久久午夜综合久久蜜桃 | 色播亚洲综合网| 尤物成人国产欧美一区二区三区| 国产一级毛片七仙女欲春2| 天堂动漫精品| 久9热在线精品视频| 国产在线精品亚洲第一网站| 熟女电影av网| 欧美日韩中文字幕国产精品一区二区三区| 97超级碰碰碰精品色视频在线观看| 757午夜福利合集在线观看| 3wmmmm亚洲av在线观看| 美女cb高潮喷水在线观看| 久久久久免费精品人妻一区二区| 久久精品国产自在天天线| 成年人黄色毛片网站| 中文资源天堂在线| avwww免费| 嫁个100分男人电影在线观看| 在线国产一区二区在线| 久久亚洲精品不卡| 亚洲欧美激情综合另类| 欧美最新免费一区二区三区 | 熟妇人妻久久中文字幕3abv| 男女下面进入的视频免费午夜| 一本久久中文字幕| 国产伦在线观看视频一区| 美女大奶头视频| 黄片大片在线免费观看| 免费人成在线观看视频色| 男人舔女人下体高潮全视频| 亚洲国产色片| www日本在线高清视频| 免费高清视频大片| 亚洲av不卡在线观看| 国产精品久久久久久亚洲av鲁大| 亚洲精品影视一区二区三区av| 国产精品亚洲美女久久久| 国产91精品成人一区二区三区| 亚洲成人久久爱视频| 又紧又爽又黄一区二区| 美女高潮喷水抽搐中文字幕| 国产精品久久久久久久久免 | 亚洲中文日韩欧美视频| 亚洲国产精品999在线| 男人的好看免费观看在线视频| 国产激情偷乱视频一区二区| 亚洲成人中文字幕在线播放| 在线观看午夜福利视频| 在线观看舔阴道视频| 草草在线视频免费看| 亚洲av一区综合| 亚洲国产中文字幕在线视频| 久久久久久人人人人人| 欧美成人a在线观看| 精品久久久久久久人妻蜜臀av| 国产野战对白在线观看| 青草久久国产| 99久久99久久久精品蜜桃| 在线视频色国产色| 久久性视频一级片| 内射极品少妇av片p| 国产欧美日韩精品亚洲av| 欧美成人免费av一区二区三区| 亚洲一区二区三区色噜噜| 男女之事视频高清在线观看| 日韩大尺度精品在线看网址| 国产精品香港三级国产av潘金莲| 18禁黄网站禁片午夜丰满| 两个人看的免费小视频| h日本视频在线播放| 乱人视频在线观看| 一级毛片女人18水好多| 久久99热这里只有精品18| 欧美丝袜亚洲另类 | 免费看美女性在线毛片视频| 久久国产精品人妻蜜桃| 日韩成人在线观看一区二区三区| 日本一二三区视频观看| 搞女人的毛片| 成年女人永久免费观看视频| 日本在线视频免费播放| 精品国产三级普通话版| 亚洲av不卡在线观看| 亚洲成人精品中文字幕电影| 51国产日韩欧美| 激情在线观看视频在线高清| 露出奶头的视频| 中文字幕高清在线视频| 一二三四社区在线视频社区8| 国产又黄又爽又无遮挡在线| 在线观看一区二区三区| 婷婷丁香在线五月| 村上凉子中文字幕在线| 97超级碰碰碰精品色视频在线观看| 大型黄色视频在线免费观看| 在线观看免费午夜福利视频| 欧美一区二区精品小视频在线| 熟女人妻精品中文字幕| 国产精品综合久久久久久久免费| 午夜老司机福利剧场| av天堂中文字幕网| 国产精华一区二区三区| 色哟哟哟哟哟哟| 男人和女人高潮做爰伦理| 欧美不卡视频在线免费观看| 制服丝袜大香蕉在线| 香蕉久久夜色| 香蕉丝袜av| 国产亚洲精品av在线| 亚洲熟妇熟女久久| 九色国产91popny在线| 男女做爰动态图高潮gif福利片| 亚洲成a人片在线一区二区| 悠悠久久av| 精品一区二区三区视频在线 | 日本a在线网址| 日日夜夜操网爽| 长腿黑丝高跟| 欧美最黄视频在线播放免费| 一进一出好大好爽视频| 黄色丝袜av网址大全| 成人国产一区最新在线观看| 91麻豆精品激情在线观看国产| 久9热在线精品视频| 精品国产超薄肉色丝袜足j| 99久久无色码亚洲精品果冻| 又黄又爽又免费观看的视频| 亚洲av第一区精品v没综合| 少妇熟女aⅴ在线视频| 观看美女的网站| 国产午夜精品论理片| 国产三级黄色录像| 午夜视频国产福利| 老汉色av国产亚洲站长工具| 制服人妻中文乱码| 亚洲av免费在线观看| 亚洲久久久久久中文字幕| 午夜福利欧美成人| 欧美黄色淫秽网站| 亚洲精品一区av在线观看| 在线免费观看不下载黄p国产 | 麻豆国产av国片精品| 99久久久亚洲精品蜜臀av| 一进一出抽搐动态| 男人舔奶头视频| 99久久无色码亚洲精品果冻| 精品一区二区三区视频在线 | 色综合站精品国产| 成人18禁在线播放| 日日夜夜操网爽| 高清日韩中文字幕在线| 精品人妻一区二区三区麻豆 | 国产精品免费一区二区三区在线| 亚洲人成电影免费在线| 淫妇啪啪啪对白视频| 国内揄拍国产精品人妻在线| 国产单亲对白刺激| 国产三级在线视频| 亚洲成a人片在线一区二区| 亚洲欧美日韩无卡精品| 蜜桃久久精品国产亚洲av| 日韩欧美一区二区三区在线观看| 精品99又大又爽又粗少妇毛片 | 狂野欧美白嫩少妇大欣赏| 亚洲av成人不卡在线观看播放网| aaaaa片日本免费| 日韩欧美精品v在线| 久久草成人影院| 热99在线观看视频| 国产aⅴ精品一区二区三区波| 美女黄网站色视频| 中文字幕久久专区| 久久久国产精品麻豆| 搡老熟女国产l中国老女人| 制服丝袜大香蕉在线| 国产毛片a区久久久久| 日本 欧美在线| 成人特级黄色片久久久久久久| 欧美色视频一区免费| 中文字幕av成人在线电影| xxxwww97欧美| 亚洲成人中文字幕在线播放| 蜜桃久久精品国产亚洲av| 宅男免费午夜| 精品一区二区三区av网在线观看| 哪里可以看免费的av片| 欧美高清成人免费视频www| 少妇裸体淫交视频免费看高清| 成人国产综合亚洲| 久久天躁狠狠躁夜夜2o2o| 国产极品精品免费视频能看的| 操出白浆在线播放| 欧美日本亚洲视频在线播放| 婷婷丁香在线五月| 亚洲激情在线av| 脱女人内裤的视频| 嫩草影院精品99| 激情在线观看视频在线高清| 嫁个100分男人电影在线观看| 99久久久亚洲精品蜜臀av| 99久久精品国产亚洲精品| 亚洲精品成人久久久久久| 黄色丝袜av网址大全| 亚洲熟妇中文字幕五十中出| 人人妻人人看人人澡| 性色av乱码一区二区三区2| 亚洲专区国产一区二区| 真人一进一出gif抽搐免费| 欧美又色又爽又黄视频| 国产精品久久视频播放| 99国产精品一区二区三区| 中国美女看黄片| 精品一区二区三区人妻视频| 在线看三级毛片| 亚洲国产精品sss在线观看| 精品人妻一区二区三区麻豆 | 国产一区二区三区在线臀色熟女| 国产精品99久久久久久久久| 亚洲aⅴ乱码一区二区在线播放| 日本黄色视频三级网站网址| 国产成人aa在线观看| 18禁黄网站禁片免费观看直播| 亚洲不卡免费看| www.999成人在线观看| 色综合婷婷激情| 亚洲欧美日韩卡通动漫| 久久精品91蜜桃| 在线观看日韩欧美| 欧美av亚洲av综合av国产av| 亚洲男人的天堂狠狠| АⅤ资源中文在线天堂| 国产亚洲精品一区二区www| 久久精品91无色码中文字幕| 精品人妻偷拍中文字幕| 国产精品久久久久久久久免 | 国产精品久久久久久亚洲av鲁大| 成年女人毛片免费观看观看9| 国产精品香港三级国产av潘金莲| 久久国产精品人妻蜜桃| a级毛片a级免费在线|