999精品在线视频,手机成人午夜在线视频,久久不卡国产精品无码,中日无码在线观看,成人av手机在线观看,日韩精品亚洲一区中文字幕,亚洲av无码人妻,四虎国产在线观看 ?

安徽省一次強對流降水過程閃電資料同化數值模擬試驗

2021-12-21 00:54:02字俁丞王京景
氣象與環境科學 2021年6期

劉 輝,章 昊,字俁丞,麥 立,王京景 ,李 志

(1.國網安徽省電力有限公司,合肥 230000;2.中國科學院大氣物理研究所,北京 100029;3.安徽繼遠軟件有限公司,合肥 230000)

引 言

強對流天氣具有突發性強、空間尺度小、生命史短和破壞力大等特點,經常伴有雷暴、大風等劇烈天氣現象[1-6]。由于電網輸變電設備普遍分布于野外、100 m高度以下,氣象環境復雜多變,因此強對流天氣常給電網的輸變電設備造成嚴重的影響[7-11]。到目前為止,強對流天氣的預報準確率仍然較低,致使電力部門應對強對流等氣象災害的應急方案不夠完善。為了提高數值模式對強對流天氣的預報能力,減少其對電網設備帶來的損失,近年來許多研究將多種非常規觀測資料(如雷達和閃電等)同化到數值模式中,用以改善模式初始場的質量,進而提升模式的預報技巧[12-17]。

閃電活動對深對流系統有很好的指示作用[18-26]。根據閃電活動情況,可以大致判斷雷暴云的發展階段、影響區域及強度等。強對流天氣過程經常伴有頻繁的閃電活動,并且閃電是雷暴云對流活動的產物。因此,閃電資料同化是提高模式初值中對流活動信息準確性的有力手段。當有閃電發生時,模式在對應空間位置上未能模擬出相當強度的對流,則可根據閃電活動的狀況,利用公式來調整該格點上的微物理量,如降水[27-31]、冰相粒子含量[13]、雷達反射率[32]等,利用初始化方法將中尺度信息引入模式,改變初始條件,進一步改進強對流系統的位置與強度,使初始場更加接近于真實場,降低初始誤差。

目前,有一些閃電資料同化的研究。Jones等[33]利用閃電資料通過張弛逼近法對中尺度模式的潛熱進行調整,提高了對亞熱帶氣旋的模擬水平。Chang等[26]分析了1998年Groundhog暴雨中對流降水和閃電的關系,并同化到模式中,其結果使短時預報(9-18 h)效果得到提高。Papadopoulos等[27]把歐洲長期的閃電資料進行同化,發展了一種把閃電實時觀測資料同化到氣象數值模式的技術。Fierro等[28]利用閃電定位資料指示對流活動的空間位置,在觀測到閃電但是模式未能模擬出深對流的格點上,根據閃電頻數和軟雹含量增加混合相態層(-20-0 ℃)的水汽,激發對流產生。郄秀書等[13]從非感應起電機制出發,認為冰相粒子濃度(冰晶、雪晶和軟雹)對云內閃電過程有重要影響,通過華北地區雷暴過程的觀測與模擬,建立了閃電頻數同冰相粒子濃度的經驗公式,根據閃電資料調整混合相態層內冰相粒子的濃度,使其更接近實際情況。張榮[14]把Fierro等[28]的方法進一步改進,首先利用閃電觀測來調整相對濕度,然后再將調整后的相對濕度同化到WRFDA系統中,進而提升了對降水的預報效果。

就目前的研究狀況來看,閃電資料同化的研究仍然較少。本文在前人研究的基礎上,以2020年6月15日安徽地區一次強對流降水天氣過程為例,利用安徽省78站逐小時站點降水資料、江淮閃電定位數據、安徽電網地面微氣象資料及中國氣象局MICAPS系統的探空資料、美國全球預測系統(GFS)逐3 h預報場資料和歐洲ERA5再分析資料,對模式啟動時刻前后10 min的閃電頻次資料進行統計,利用其與相對濕度的經驗關系式[14]和WRF-3DVAR資料同化系統,間接地把閃電資料同化到WRF模式的初值場中,并設計了地面和探空資料同化與否的對比試驗,分析閃電資料同化在模式強對流降水天氣過程模擬預報中的作用,為閃電資料同化應用提供科學參考。

1 強對流過程概況

2020年6月15日00:00(UTC,下同)環流形勢場中,中高緯500 hPa高度場為兩槽一脊的分布(圖略),低壓槽分別位于烏拉爾山以東地區和我國東北地區,高壓脊位于貝加爾湖以西附近,這種環流形勢有利于冷空氣的南下。至15日12:00,環流形勢維持兩槽一脊分布(圖1a),但低壓槽和高壓脊有所東移,西北太平洋副熱帶高壓(以下簡稱副高)315 dagpm特征線異常偏西,西脊點西伸至我國西南地區,安徽處于副高西北側外圍水汽輸送的強輻合區中和冷暖空氣的交匯地帶,安徽南部對流有效位能(CAPE)超過1600 J·kg-1(圖1b),存在明顯的對流不穩定能量。

圖1 2020年6月15日12時500 hPa高度場(a)和700 hPa高度場、水平風場及對流有效位能(b)、安徽省16日12時24 h累積降水(c)和15日12時1 h累積降水及前后各10 min閃電定位分布(d)(b)中紅色粗線為副高315 dagpm特征線;(d)中“+”為閃電定位標記

在上述環流背景下,6月15日12:00至16日12:00,安徽發生了一次區域性強對流天氣過程,在安徽西部形成一條自北向南的狹長強雨帶,強度超過25 mm,中心則超過50 mm;在安徽北部和東北部也分別出現一個范圍較小的強雨帶,而安徽東南部降水則明顯偏弱(圖1c)。此次降水過程伴有明顯的閃電活動。如15日12:00前后各10 min共捕捉到760次閃電,主要出現在安徽西南部,并且與1 h累積降水空間分布具有很好的一致性(圖1d)。鑒于此,本文選取12:00前后各10 min的閃電頻次進行同化模擬,具有較好的典型性和代表性。

2 數值模式簡介和閃電資料同化方法

2.1 數值模式簡介

本文采用WRF模式4.2版本非靜力方案進行兩重雙向嵌套網格模擬,試驗中心位于32.2°N、117.3°E。考慮到計算資源有限,水平分辨率分別為18 km、6 km,水平方向格點數分別為401×401、121×121;模式所有區域垂直方向共51個σ層,模式層頂50 hPa;時間步長設為90 s。模式物理過程采用Morrison雙參數化云微物理方案、Kain-Fritsch積云參數化方案、MYJ邊界層方案、Noah LSM 陸面模式、RRTMG長波和短波輻射方案,輻射方案每5 min計算一次。側邊界條件由GFS 0.25°×0.25° 每3 h一次的預報場提供。采用WRFDA(WRF Data Assimilation)三維變分同化系統對觀測資料進行同化,同化在兩個區域內進行。文中所有數值模擬結果的分析均在D02區域進行。

2.2 WRF-3DVAR同化系統

三維變分同化方案把同化歸結為如下目標函數的極小化問題[34]:

(1)

式中,X是分析變量,Xb是背景場,yo是觀測場,H是觀測算子(代表模式空間向觀測空間的一種映射);B-1是背景誤差協方差矩陣的逆,O-1是觀測誤差協方差矩陣的逆,上標T表示矩陣轉置。

背景誤差協方差在同化過程中起著關鍵性作用,它決定著背景場和觀測場之間的相對重要性。雖然WRF-3DVAR同化系統中包含了由NCEP統計得到的全球適用的背景誤差協方差(CV3),但不宜直接用于中小尺度天氣過程。研究表明,基于實際模式框架和物理過程統計得到的背景誤差協方差(CV5)對同化性能改進明顯[35-36]。為此,本文利用2020年7月1日兩次的預報結果(共獲得62個統計樣本),采用NMC方法[37],用同一時刻的24 h預報減去12 h預報作為背景誤差,統計獲得了背景誤差協方差(CV5)。

2.3 閃電資料同化方法

由于閃電觀測量不是模式常規變量,因此閃電資料不能直接用于模式初始化。同化閃電資料前,首先采用經驗或半經驗的關系式,將閃電觀測反演成模式變量或相關的診斷量。利用Fierro等[28]建立的閃電頻數與水汽及霰混合比的關系式,并將其作為閃電同化的觀測算子,將閃電資料反演成水汽混合比。該關系式為

(2)

式中,Qv為調整后得到的水汽混合比,Qsat為飽和水汽混合比,Qg為霰混合比;tanh為雙曲正切函數;X為閃電頻數;A、B、C、D、a為常參數,A=85%,B=0.20,C=0.02,D=0.25,a=2.20。

然后利用RH=Qv/Qsat對(1)式進行數學變換[14],得到關系式:

(3)

式中,RH為相對濕度。將關系式(3)作為閃電同化的觀測算子引入WRF-3DVAR同化系統。首先利用閃電觀測來調整相對濕度,然后再將調整后的相對濕度同化到WRF-3DVAR同化系統中。常數A用來設置關系式(3)適用的相對濕度閾值條件,只有當模式網格柱上的相對濕度小于A,且霰混合比Qg小于3 g·kg-1時,才對該網格柱上0-20 ℃溫度層之間的水汽進行調整。

2.4 不同資料同化試驗方案設計

利用WRF-3DVAR同化系統對2020年6月15日12:00的地面、探空和閃電資料進行同化,根據同化的不同觀測資料類別,設計了4種同化試驗方案。

(1)控制試驗(CTL):不同化任何資料,僅以6月15日12:00背景場作為初值驅動WRF模式進行24 h預報。其中背景場的產生過程如下:以GFS模式6月15日00:00起報的6 h預報場作為初值,驅動WRF模式進行6 h預報至6月15日12:00,以該時刻預報場作為控制試驗的背景場。

(2)閃電資料同化試驗(LTN):其背景場同CTL的背景場一致,采用WRF-3DVAR同化系統同化6月15日12:00閃電頻數資料(圖1a),以此作為初值驅動WRF模式進行24 h預報。其中閃電頻數統計的是同化時刻前后各10 min內的閃電發生次數。

(3)地面探空資料同化試驗(SH):基于6月15日12:00的背景場,采用WRF-3DVAR同化系統同化該時刻的地面和探空資料,以此作為初值驅動WRF模式進行24 h預報。該試驗背景場與前2個試驗略有不同,其產生過程如下:以GFS模式6月15日00:00起報的6 h預報場作為背景場,同化15日06:00地面探空資料后作為初值驅動WRF模式進行6 h預報至6月15日12:00,以該時刻預報場作為SH試驗的背景場。

(4)地面探空和閃電資料同化試驗(SHL):在SH基礎上增加6月15日12:00閃電頻數資料同化。

2.5 評估檢驗方法

首先采用線性插值法,將模式模擬的結果插值到安徽地區的自動觀測站上,形成一個與觀測場對應的站點預報場。根據24 h的降水量,將降水量分為小雨(0.1 mm≤R<10.0 mm)、中雨(10.0 mm≤R<25.0 mm)、大雨(25.0 mm≤R<50.0 mm)、暴雨以上(R>50.0 mm)等4個等級,計算TS評分,來檢驗不同模式在降水分級預報中的表現。此外,采用空間相關系數和均方根誤差 RMSE 對各預報時次累積降水進行評估。TS、RMSE計算公式分別為

(4)

(5)

式中,i代表不同降水量級,NAi為預報正確次數,NBi為空報次數,NCi為漏報次數;Fi為預報值,Oi為實況值,N為總站點數。

3 結果分析

3.1 同化對模式初值的改進

本文采用的閃電資料同化方法,是利用閃電觀測來調整模式相對濕度,然后再將調整后的相對濕度進行同化的方法,以期盡可能準確地為模擬強對流提供有利的環境條件。為清晰反映閃電資料同化對模式初值的改進,對初始時刻(2020年6月15日12:00)4個試驗的物理量場分布進行對比分析。

對比LTN和CTL試驗發現,兩者主要區別在于同化閃電資料之后,在閃電活動區相對濕度明顯增加,而在其他區域則變化不明顯(圖2 a、b);SH試驗在CTL試驗基礎上同化了地面微氣象資料和探空資料,其相對濕度水平分布較CTL試驗的也有所增加,但增加的幅度較LTN同化試驗的小(圖2c);SHL試驗在SH試驗基礎上同化閃電資料之后,在閃電活動區相對濕度則明顯增加(圖2d)。4個試驗與ERA(把ERA再分析資料視為觀測資料)對比發現(圖2e),觀測場的相對濕度皆大于4個試驗的數據,但同化閃電資料之后(圖2b、d),模式相對濕度場更加接近觀測。這說明閃電資料同化能在同化區域顯著改進該區域的相對濕度。

圖2 安徽省2020年6月15日12時模式初始場和再分析資料的500 hPa相對濕度分布(a)CTL試驗,(b)LTN試驗,(c)SH試驗,(d)SHL試驗,(e)ERA;矩形框為閃電活動區

從閃電活動區相對濕度和溫度的緯度-高度垂直剖面也進一步證明了這一現象(圖3)。由LTN和CTL試驗發現,同化閃電資料之后,在閃電活動區域(29.5°-31.5°N)700-200 hPa高度層間,相對濕度明顯增強,其中相對濕度在500-400 hPa調整最明顯,最大調整量達24%(圖3a、b),而溫度調整量很小(圖3f、g)。在無閃電區域(31°N以北)相對濕度和溫度的調整沒有閃電活動區的明顯。

圖3 安徽省2020年6月15日12時初始場沿閃電活動區(117.3°E)的相對濕度(a-e)和溫度(f-j)經向-垂直剖面分布(a)(f)CTL試驗,(b)(g)LTN試驗(等值線為LTN與CTL的差值),(c)(h)SH試驗(等值線為SH與CTL的差值),(d)(i)SHL試驗(等值線為SHL與SH的差值),(e)(j)ERA(等值線為ERA與SHL的差值)

對比SH和CTL試驗發現,同化地面探空資料后,31°N以北的500 hPa以下高度層相對濕度明顯減小,說明地面探空資料同化之后使大部分區域高度層的濕度差減小(圖3 b、c),這可能是造成SH試驗降水普遍減少的原因。SH試驗中溫度在中高層400-300 hPa和低層900-850 hPa調整最明顯,較LTN試驗的高1.0 ℃(圖3g、h)。

對比SHL和SH試驗發現,同化閃電資料之后,在閃電活動區(30.5°N)以南700-200 hPa高度之間相對濕度均明顯增加,并且500-400 hPa增加最明顯,而在無閃電區域(30.5°N以北)相對濕度有所減小(圖3c、d),溫度則在30.5°N中上層500-300 hPa調整最明顯,超過0.6 ℃(圖3h、i)。

從觀測場和SHL試驗結果來看,SHL試驗多數高度層中相對濕度都明顯低于觀測值。同化閃電資料后,在閃電活動區(29°-31°N)的800-400 hPa高度層中,相對濕度和觀測值的差異相對較小(圖3d、e),溫度則在400-300 hPa中和觀測值的差異最小(圖3i、j)。說明在同化閃電資料以后,相對濕度和溫度得到顯著改善。其他一些物理量(如u、v風場)的初始場水平和垂直分布也有類似的狀況(圖略)。

3.2 模式的物理量變化特征

考察模式積分輸出結果可以發現,4個試驗均能模擬出西南氣流,但對造成安徽強對流降水的短波槽系統的模擬存在差異。

對比LTN和CTL試驗發現,模式積分6 h,兩種模式模擬的水平風速差異較小,降水強度均小于5 mm(圖4a1、b1)。模式積分12 h,安徽西北部有短波槽東移,北部均出現強降水帶,強度超過30 mm(圖4a2、b2)。LTN試驗在西北部降水大值區內,模擬的風速明顯偏大(圖4b2),但位溫明顯減小(圖略),暖濕氣流較弱,降水強度較CTL的偏弱。模式積分18 h,短波槽加深,但位溫仍偏低(圖略),因此降水強度弱于CTL試驗模擬的降水強度(圖4a3、b3)。此外,LTN試驗在31°N附近出現一強降水中心。模式積分24 h,短波槽持續東移,降水均開始減少,LTN試驗在北部模擬的風速大于CTL試驗模擬的風速(圖4a4、b4),模擬的降水較CTL試驗的偏少;CTL試驗分別在北部和中東部地區模擬出多個強降水中心。

對比SH和CTL試驗發現,模式積分6 h,SH試驗在安徽西部和東北部模擬的低空風速較CTL試驗的明顯增加(圖4 c1),模擬的33°N以北的位溫顯著增高(圖略),因此出現一降水中心(圖4c1)。模式積分12 h,SH試驗在西北部也模擬出一強降水中心,但模擬的短波槽明顯偏弱(圖4a2、c2),西北部風速明顯偏大,導致該區域降水偏弱。模式積分18 h,SH試驗模擬的短波槽較弱,水平風速大于CTL試驗模擬的風速(圖4a3、c3),模擬的位溫卻低于CTL試驗的(圖略),因此造成大部分區域降水偏少。模式積分24 h,SH試驗模擬的水平風速整體偏強(圖4a4、c4),但短波槽仍然較弱,導致降水偏少。

SHL和SH試驗與ERA對比發現,模式積分6 h,SHL和SH試驗模擬的水平風場(圖4c1、d1)和位溫(圖略)差異較小,降水強度也偏弱。ERA的風速均大于SHL試驗模擬的風速(圖4d1、e1),在安徽南部SHL試驗沒有模擬出強降水。模式積分12 h,ERA和SHL試驗在北部出現了短波槽,有風場的輻合,SHL試驗模擬的西北部風速明顯大于SH試驗的模擬值,但又小于ERA的值(圖4c2、d2、e2),且模擬的位溫偏低(圖略),導致模擬的降水強度弱于SH試驗和實況的強度,但SHL試驗能夠在西部模擬出一條南-北向的狹長雨帶。模式積分18 h,短波槽東移并且明顯減弱,SHL和SH試驗模擬的降水都逐漸減少,模擬的雨帶位置與實況基本重合。其中,SHL試驗模擬的降水量明顯小于ERA值(圖4d3、e3),也較SH試驗的少,大部分地區降水量小于15 mm。模式積分24 h,隨著短波槽的東移,北部的降水減少,SHL試驗模擬的全區域風速小于ERA的值(圖4d4、e4)。在中南部31°N附近,SHL試驗模擬結果和實況均出現了降水(圖4c4、d4、e4)。

圖4 安徽省2020年6月15日12時-16日12時模式模擬和實況逐6 h累積降水分布(等值線)、700 hPa風場(矢量)及風速差值場(陰影)(a1)-(a4)CTL試驗,(b1)-(b4)LTN試驗(陰影為LTN與CTL的差值),(c1)-(c4)SH試驗(陰影為SH與CTL的差值),(d1)-(d4)SHL試驗(陰影為SHL與SH的差值),(e1)-(e4)ERA(陰影為ERA與SHL的差值);降水單位:mm

3.3 降水模擬的改進

3.3.1 24 h雨帶的模擬

對比實況(圖1c)和4次試驗模擬的24 h累積降水(圖5),CTL試驗能夠模擬出安徽西北部的雨帶(圖5a),但模擬的降水強度明顯偏大,強度超過35 mm,且模擬的降水范圍也較實況的大,擴展至整個安徽北部;模擬的安徽西南部降水較少,不能很好地模擬出安徽西南部的雨帶,在中南部和東南部模擬出一較強的雨帶,與實際降水的空間相關系數僅為0.15。說明CTL對24 h累積降水的模擬能力較弱。LTN試驗(圖5b)模擬出的雨帶分布特征與CTL試驗的雨帶較為類似,模擬的雨帶強度偏弱,空間相關系數提升至0.3,但仍沒有模擬出西南部的降水,表明單獨同化閃電資料對改善模式降水模擬有一定正效果,但改進程度相對有限。

從SH試驗雨帶的分布情況來看(圖5c),模式能夠較好地模擬出安徽西北部的雨帶中心,中心降水強度超過50 mm,優于前兩個試驗模擬效果,與實況降水的空間相關系數也增至0.37,超過了0.05的顯著性檢驗,但仍然不能模擬出安徽西南部的雨帶。SHL試驗能模擬出西北部的雨帶中心位置(圖5d),但模擬的降水強度和面積略偏小,對于最北部和東北的雨帶及降水強度也有一定的模擬能力,與前三個試驗相比較,模式還較好地模擬出了西南部的雨帶中心位置,與實況降水的空間相關系數為0.38,超過了0.05的顯著性檢驗。這說明在地面和探空資料同化的基礎上增加閃電資料的同化,能夠顯著地改善降水的預報效果。

圖5 安徽省2020年6月16日12時模式模擬的24 h累積降水量分布(a)CTL試驗,(b)LTN試驗,(c)SH試驗,(d)SHL試驗

3.3.2 24 h的降水評估檢驗

由各等級降水的 TS評分結果可知(圖6a),在所有降水等級中,CTL試驗的評分值均為最低,說明無同化的模擬效果最差;而SHL試驗的均為最高,說明閃電和地面探空資料的同化能夠提高不同量級的降水預報能力,其中對中雨及其以上等級降水的預報效果提升尤為明顯;LTN的小雨和大雨評分值要稍高于SH試驗的評分,說明僅同化閃電資料對小雨和大雨的預報能力也有提升。

采用均方根誤差RMSE對模式模擬的24 h累積降水進行檢驗(圖6b)。CTL試驗在模式開始積分時,RMSE為17.24 mm,隨著時間推移,均方根誤差呈準線性逐漸增加。LTN試驗開始時,RMSE較小,9 h后顯著增加,但均較CTL試驗的低,說明僅進行閃電資料同化,即可提升模式對降水的模擬能力;SH試驗開始時RMSE較小,6-24 h RMSE保持在15~20 mm以內,這表明地面和探空的資料同化能夠顯著改進模式的降水模擬性能。SHL試驗開始時RMSE幾乎與SH試驗的一致,但是在6 h后均小于SH試驗值,說明在地面和探空資料同化上增加閃電資料同化可以進一步提升降水模擬的性能。

圖6 安徽省2020年6月15日12時-16日12時模擬的24 h累積降水TS評分(a)和模擬的逐小時降水量的RMSE演變(b)

4 結論與討論

本文以2020年6月15日12:00-16日12:00安徽一次伴有雷暴的強對流天氣為研究對象,設計閃電資料、地面和探空資料對比同化試驗,對比分析了閃電資料同化對模式初值改進和24 h累積降水的模擬能力,主要結論如下:

(1)單獨同化閃電資料的LTN試驗能夠在一定程度上對初值進行改進,使初始場大氣相對濕度增加,溫度降低,但改進能力有限;而同化地面和探空資料的SH試驗則使初始場中相對濕度垂直差異減小,溫度增加,改善效果好于LTN試驗的;同時同化地面、探空和閃電資料的SHL試驗則能夠使同化區域的相對濕度和溫度等更加趨近于觀測值。

(2)CTL試驗未能較好地模擬出降水的空間分布結構,LTN試驗對此有所改善,但效果有限;SH試驗模擬的降水空間分布改進明顯,但仍未模擬出安徽西南部的一個強雨帶;而SHL試驗則進一步改進了降水空間結構的模擬,并較好地模擬出這一雨帶。

(3)對安徽地區不同等級的24 h累積降水模擬結果進行TS評分,中雨和暴雨的評分值CTL的LTN的>SH的>SHL的。因此,在地面探空同化基礎上增加閃電資料的同化,能夠顯著減小模式降水模擬的誤差,使模擬結果更加接近實況。

本文從一次強對流天氣過程的閃電同化模擬出發,對比分析了閃電資料、地面和探空資料同化模擬試驗,發現閃電資料同化對強對流降水的模擬有一定改進,并對結果進行了檢驗。但本文僅同化了一次閃電過程,結論是否具有普遍性,還需要今后進行批量對比試驗來進一步檢驗。

主站蜘蛛池模板: 国产成人一区在线播放| 成年人福利视频| 亚洲电影天堂在线国语对白| 美女内射视频WWW网站午夜| 国产全黄a一级毛片| 在线观看免费黄色网址| 国产拍揄自揄精品视频网站| 高潮毛片免费观看| 91九色视频网| 人妻一区二区三区无码精品一区 | 亚洲综合久久一本伊一区| 日本不卡在线| 亚洲AⅤ波多系列中文字幕| 曰AV在线无码| 麻豆AV网站免费进入| 国产h视频在线观看视频| 99精品国产电影| 国产理论最新国产精品视频| 福利在线不卡一区| 色综合天天操| 精品中文字幕一区在线| 久久香蕉国产线看观看式| 激情综合激情| 国产精品网址在线观看你懂的| 国产视频久久久久| 国产成人亚洲日韩欧美电影| 精品久久久久无码| 欧美人与牲动交a欧美精品| 国产精品流白浆在线观看| 国产SUV精品一区二区6| 大学生久久香蕉国产线观看| 国产杨幂丝袜av在线播放| 91精品专区| 精品一区二区三区自慰喷水| 国产黄网站在线观看| 亚洲第一网站男人都懂| 巨熟乳波霸若妻中文观看免费 | 国产系列在线| 成人午夜亚洲影视在线观看| 国产va欧美va在线观看| 欧美中文字幕在线视频| 日韩毛片基地| 精品自窥自偷在线看| 亚洲视频a| 亚洲 欧美 偷自乱 图片| 国产精品久久久久无码网站| 亚洲首页国产精品丝袜| 在线看片免费人成视久网下载| 精品国产网| 波多野结衣国产精品| 亚洲AV无码精品无码久久蜜桃| 婷婷亚洲天堂| 精品五夜婷香蕉国产线看观看| 永久免费AⅤ无码网站在线观看| 国产av一码二码三码无码| 久久精品嫩草研究院| a级免费视频| 国产视频一区二区在线观看| 亚洲专区一区二区在线观看| 亚洲av中文无码乱人伦在线r| 亚洲欧美精品在线| 五月激激激综合网色播免费| 欧洲亚洲欧美国产日本高清| 欧美一级专区免费大片| 久久精品无码中文字幕| 激情综合婷婷丁香五月尤物| 91精品人妻互换| 国产精品视频3p| 国产精品人人做人人爽人人添| 日韩av电影一区二区三区四区| 欧美精品亚洲日韩a| 伊人久久婷婷| 成人精品午夜福利在线播放| 精品视频在线一区| 国产九九精品视频| 一级爆乳无码av| 亚洲bt欧美bt精品| 亚洲国产欧美中日韩成人综合视频| 久久综合亚洲鲁鲁九月天| 啪啪啪亚洲无码| 狠狠色香婷婷久久亚洲精品| a毛片在线|