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

背景誤差尺度分離與多尺度混合濾波技術在CMA-MESO 3 km系統的應用*

2022-02-14 12:28:42徐枝芳王瑞春
氣象 2022年12期
關鍵詞:背景特征水平

徐枝芳 王瑞春

1 中國氣象局地球系統數值預報中心,北京 100081 2 中國氣象科學研究院災害天氣國家重點實驗室,北京 100081 3 國家氣象中心,北京 100081

提 要: 為提高CMA-MESO 3 km系統降水預報能力,采用二維離散余弦變換對2018年6月2日至8月31日3個月格點背景誤差樣本結合模式分辨率和天氣系統尺度范圍劃分進行3種尺度分離,并對這3種尺度背景誤差樣本分別進行水平協相關尺度擬合,通過采用3個不同水平特征相關尺度的遞歸濾波器在CMA-MESO 三維變分同化系統中實現3種擬合水平協相關尺度應用,替代業務測試系統單一水平特征相關尺度,開展個例和連續試驗分析。研究結果表明,采用二維離散余弦變換尺度分離背景誤差樣本的3種水平特征相關尺度垂直結構相似,水平特征尺度的水平尺度相隔幾十至幾百千米。擬合水平特征相關尺度在CMA-MESO 3 km系統應用結果顯示,3種水平特征相關尺度試驗對u風、v風、濕度分析有明顯正影響,分析更接近實況,對溫度分析影響較小;對降水預報有改善,冷啟動預報前6 h的TS評分提高明顯,偏差(Bias)減小向1靠近,暖啟動24 h逐 6 h降水預報TS評分都小幅提升,Bias差異不大。

引 言

隨著我國經濟高速發展,公眾及高敏感行業對暴雨等災害天氣的精準預報要求也越來越強烈。暴雨等災害天氣系統多是由多種尺度系統相互作用條件下產生,如導致暴雨產生的系統中,大尺度系統為暴雨提供了有利的環境條件,而中小尺度系統是造成暴雨的直接系統(陸漢城和楊國祥,2004)。采用數值模式預報暴雨等災害天氣系統,雖然動力與物理過程的描述存在著一定不足,但對暴雨等系統活動的精準預報具有一定的優勢。提高暴雨等災害天氣數值預報能力可從兩方面著手:一方面是提高模式分辨率,開發適合高分辨率模式的物理過程參數化方案和動力框架,增強模式預報精度;另一方面是改善模式初始條件,提高模式對中小尺度系統的刻畫(Vendrasco et al,2016)。資料同化給數值模式提供一個更為精準的初值場(陳東升等,2004),是改善數值模式預報技巧的一種有效手段,因而資料同化方案提供的初值需要在不破壞大尺度環流形勢場的同時,盡可能多地引入中小尺度信息。當前,各主要業務數值預報中心,主要采用變分資料同化方案為模式提供初值(Lorenc,2003;Bannister,2008)。在變分資料同化中,濾波特點以及觀測信息的傳播方式均由背景誤差協方差矩陣B決定,合理的B矩陣是做好變分同化的關鍵。背景誤差協方差矩陣在資料同化系統中控制信息從觀測位置向四周傳播的方式,并決定了模式變量之間在動力上是否協調一致。因此,同化分析增量的空間結構和多變量關系結構取決于背景誤差協方差的結構。

由于B是一個維數很高的矩陣,對其直接求逆不現實,在近似得到B-1的研究中,許多科研人員通過變量變換對目標函數進行預調節,B矩陣在預調節中被分解為幾個部分,其中水平背景場誤差協方差由各向同性遞歸濾波表示。遞歸濾波法不需要顯式構造B矩陣,可以顯著減小計算開銷和內存(李冬等,2011)。英國氣象局第一次將遞歸濾波方法引入到區域預報模式(Purser and McQuigg,1982)。目前,遞歸濾波被廣泛應用在區域三維變分系統 (Hayden and Purser,1995; Purser et al,2003a;2003b)。遞歸濾波中特征尺度是一個重要參數,特征尺度越大濾波影響范圍越大,對于觀測信息較少的地區可以擴大觀測信息的影響(Hayden and Purser,1995)。但是,特征尺度過大,同化結果將趨于平滑,無法提取觀測的短波信息;特征尺度過小,則不能準確地提取長波信息。何光鑫等(2011a;2011b)利用多元正態分布的可加性進行遞歸濾波擬合(Wilks,2006),將不同特征尺度遞歸濾波進行擬合,擬合后的多尺度遞歸濾波在保持原有大尺度信息的基礎上,清晰顯示出了更多中小尺度的信息。吳洋等(2018)在保證原遞歸濾波器水平特征尺度流函數、勢函數和非平衡質量場500 km和濕度場200 km 基礎上,依據統計結果進行多次線性擬合出3種水平特征尺度經驗性結果,然后采用3種特征尺度遞歸濾波器在分析和預報中獲得了更多的α中尺度信息,使得形勢場及降水預報技巧明顯改善。莊照榮和李興良(2021)通過分析高斯模型和尺度疊加高斯模型的空間特征,及拉普拉斯算子和譜響應函數的特征,以原方案流函數、勢函數和非平衡質量場水平特征尺度值為中心形成等差序列水平特征尺度組,研究了尺度疊加高斯相關模型的特征及其在三維變分同化系統中的應用效果。他們的工作主要研究了多尺度疊加對分析和預報影響,但如何客觀地從背景誤差協方差中提取和估計多種特征尺度長度的工作相對較少。莊照榮等(2018;2020)通過采用二維離散余弦轉換(2D-DCT)(Denis et al,2002;鄭永駿等,2008),將全球T639模式和區域GRAPES模式分析增量進行譜分解,獲得全球和區域模式產品不同尺度信息,然后將全球大尺度信息部分與區域中小尺度信息部分混合組成新分析增量,使得分析和降水預報明顯改善。本文參考莊照榮等(2018;2020)二維離散余弦轉換獲得全球和區域模式產品不同尺度信息方式,采用二維離散余弦變換對格點背景誤差進行波譜展開后,對背景誤差按照尺度大小分為3部分進行尺度分離和水平相關尺度擬合,然后將擬合的隨高度變化背景誤差水平相關特征尺度在3 km中尺度天氣數值預報系統CMA-MESO(原GRAPES區域中尺度數值預報系統GRAPES-MESO)進行個例和連續試驗,并與業務試驗系統進行對比分析。

2 方案設計

CMA-MESO 3 km 系統采用三維變分同化方法(馬旭林等,2009;薛紀善和陳德輝,2008),目標函數J為:

(1)

式中:x表示分析場,xb為背景場,B為背景場誤差協方差矩陣;yo為觀測向量;H是把大氣狀態投影到觀測空間的觀測算子;E為儀器觀測誤差協方差矩陣,F為代表性誤差協方差矩陣。背景場誤差協方差矩陣B、觀測算子H、觀測資料質量控制是三維變分同化三個核心問題,任何一個問題沒處理好,提供的分析場就會影響數值結果。上述目標函數最小化過程涉及背景誤差協方差矩陣B求逆。鑒于B是超大規模矩陣,直接求逆計算代價比較高,且有時無法求逆,導致極小化計算難以收斂,通過引入控制變量變換預處理則可以避免背景誤差協方差矩陣求逆,設x-xb=Uw,B=UUT,目標泛函J寫為如下形式:

(2)

式中:w為新的極小化控制變量,d=H(xb)-yo為新息向量,R=E+F為觀測誤差協方差。預處理變換中將誤差協方差矩陣U分裂為3個獨立的矩陣U=UpUhUv,Up為物理變換算子,為不同分析變量之間的誤差交叉協相關,控制著觀測信息在不同變量之間的傳播;Uv為誤差垂直變換算子,由經驗正交函數(EOF)方法實現;Uh為誤差水平變換算子。現有業務系統一般均假設背景水平誤差滿足各向同性的高斯分布,Uh采用遞歸濾波實現該誤差模型(薛紀善和陳德輝,2008;莊照榮等,2019)。CMA-MESO 3 km 三維變分同化系統中,誤差高斯分布函數定義為:

(3)

對于一個Ni×Nj的二維變量場f(i,j),2D-DCT的直接變換和反變換定義為:

(4)

(5)

其中

(6)

(7)

式中:f(i,j)為二維變量場在格點(i,j)的值,而F(m,n)為(m,n)維波數對應的譜系數。 2D-DCT對一個Ni×Nj維的物理量場f(i,j)的正變換可以獲得Ni×Nj維的譜數場F(m,n)。經過變換的每一個譜系數的位置,即二維波數(m,n)需要與單一波長L相聯系,因此對于一個Ni×Nj的二維譜系數場,則有:

(8)

式中:Δ為模式格點分辨率。波長L與對應的二維波數k為:

(9)

結合天氣系統尺度劃分系統和3 km模式系統關注的天氣系統,將背景誤差按照以下3種尺度區間(≥500 km,100~500 km,≤100 km)進行分離。

隨機選取一個時次某層溫度背景誤差,從圖1可見,3種尺度背景誤差分布特征及大值區與尺度分離前(圖1a)大體一致;≥500 km背景誤差(圖1d)分布相對光滑, 100~500 km區間背景誤差(圖1c)多了些小擾動,分布不如≥500 km的背景誤差光滑,小擾動少于≤100 km(圖1b)下尺度背景誤差。3種尺度分離后的結果基本反映了3種尺度信息,說明尺度分離合理。

圖1 隨機樣本量單層溫度誤差尺度分離前后的水平分布(a)原始,(b)≤100 km,(c)100~500 km,(d)≥500 kmFig.1 Horizontal distribution of a random layer temperature error sample before and after separation(a) original, (b) ≤100 km, (c) 100-500 km,(d) ≥500 km

將背景誤差尺度分離后3組樣本進行水平特征相關尺度擬合,結果見圖2。由圖2可見,擬合的3種水平特征相關尺度隨高度變化特征相似,水平特征尺度長度大小不同,低層水平特征尺度最大接近200 km,中間約為100 km,最小約為40 km。業務試驗系統的水平特征相關尺度隨高度變化特征與擬合的3種水平特征尺度不同,u、v風場高層、溫度中層以及濕度高層差異較大,這可能與業務試驗系統采用的背景誤差樣本量不同有關。在低層,業務試驗系統單一水平特征尺度長度和3種水平特征尺度長度中間的長度接近。

圖2 (a)u風場,(b)v風場,(c)溫度和(d)比濕的水平特征相關尺度長度隨氣壓變化(實線:業務測試系統/CTL,長虛線:≥500 km,中長虛線:100~500 km,短虛線:≤100 km)Fig.2 Horizontal correlation length changes with pressure for (a) u component,(b) v component,(c) temperature and (d) specific humidity(solid line: operational testing/CTL, long dash line: ≥500 km, medium dashed line: 100-500 km, short dashed line: ≤100 km)

3 CMA-MESO數值試驗結果

3.1 CMA-MESO簡介及試驗方案

本文試驗采用的是CMA-MESO 4.4版,模式系統水平分辨率為3 km,垂直不等間距51層,由于計算資源問題,試驗范圍為中國中東部區域(17°~36.8°N、102°~126.9°E)。 采用的三維同化系統控制變量為u、v風場,地面氣壓ps,溫度T和相對濕度RH。同化分析資料為探空報(u,v,T,RH),地面報(ps,RH),船舶報(ps,RH),飛機報(u,v,T),云導風,地基掩星反演可降水量(GPS/PW),雷達速度方位顯示反演風廓線(VAD風),雷達徑向風,風廓線雷達風,基于云分析系統(朱立娟等,2017)的雷達和衛星資料,松弛逼近方法(nudging)同化的地面自動站降水資料。試驗為每3 h同化分析24 h預報,00時和12時(世界時,下同)包含一次冷啟動和一次暖啟動分析預報。采用美國國家環境預報中心(NCEP) 0.5°×0.5° FNL 6 h預報場做背景和側邊界條件資料,同化分析預報試驗時段為2018年7月1—31日。

數值模擬試驗分為2組:(1)CTL試驗:采用業務試驗系統用水平特征相關尺度參與同化分析的控制試驗;(2)MF試驗:采用3種擬合水平特征相關尺度參與同化分析的敏感試驗。

3.2 數值試驗結果

3.2.1 個例試驗

2018年7月4—7日在低空急流和江淮氣旋的影響下,我國中東部大部地區出現較強降雨過程(張夕迪和孫軍,2018),河南東南部、安徽大部、湖北中東部、貴州中北部、重慶東部、湖南大部、江蘇大部、江西中北部、廣西北部、廣東中北部出現暴雨,其中河南東部、安徽大部、江蘇大部、江西北部、廣西北部出現大暴雨,江西北部個別站點特大暴雨。圖3為7月5日00時至6日00時逐6 h降水實況,5日00—06時30°N上下有一條雨帶,江蘇、安徽、湖北等地降水量較大,然后降雨帶逐步東移南壓,12 h后強降水自北向南逐步移至江西省北部地區。北京地區在5日06—18時有降水發生。

圖3 2018年7月5日00時至6日00時逐6 h降水實況(a)00—06時,(b)06—12時,(c)12—18時,(d)18—00時Fig.3 The 6 h accumulated rainfall from 00:00 UTC 5 to 00:00 UTC 6 July 2018(a) 00:00-06:00 UTC, (b) 06:00-12:00 UTC,(c) 12:00-18:00 UTC, (d) 18:00-00:00 UTC

沿120°E對MF試驗和CTL試驗7月5日00時的濕度、溫度和風場分析增量剖面進行對比分析發現(圖4),兩組試驗的濕度、溫度和風場分析增量差異較小,濕度分析增量差異不超過6%,風場分析增量差異不超過1.2 m·s-1,溫度分析增量差異不超過0.35℃。中低層(500 hPa以下),主要降水區(26°~34°N)的相對濕度分析增量MF試驗高于CTL試驗,而700 hPa以下溫度分析增量MF試驗小于CTL試驗, 2組試驗在29°~30°N及33°N附近u風場分析增量差為正偏差,在29°N附近和31°~32°N附近v風場分析增量為正偏差。

圖4 2018年7月5日00時兩組試驗分析增量差(MF試驗減CTL試驗)沿120°E剖面(a)相對濕度(單位:%),(b)溫度(單位:℃),(c)u風場(單位:m·s-1),(d)v風場(單位:m·s-1)Fig.4 The increment difference profile between MF experiment and CTL experiment along 120°E at 00:00 UTC 5 July 2018(a) relative humidity (unit: %), (b) temperature (unit: ℃),(c) u-component (unit: m·s-1), (d) v-component (unit: m·s-1)

從逐6 h降水圖來看(圖5,圖6),MF試驗和CTL試驗都預報出了此次降水過程,降水分布較實況強度偏強,范圍偏大。MF試驗和CTL試驗5日00—06時時段降水預報分布差異不明顯,江蘇和安徽境內降水較實況偏強,06—12時和12—18時在江蘇、安徽和湖北省區域降水預報逐漸出現細微差異,MF試驗更接近實況。5日18時至 6日00時預報雨帶較實況略偏南,江蘇北部地區預報偏大,MF試驗效果略差。不同預報時段的降水分布結果顯示,MF試驗在北京地區的降水預報更接近實況。由此可見,MF試驗和CTL試驗的分析差異雖小,但對降水預報有影響,MF試驗降水預報更接近實況。

圖5 同圖3,但為CTL試驗結果Fig.5 Same as Fig.3, but for CTL experiment

圖6 同圖3,但為MF試驗結果Fig.6 Same as Fig.3, but for MF experiment

3.2.2 連續試驗

上文通過個例進行了MF試驗和CTL試驗結果分析,下面進行連續試驗(2018年7月1—31日)分析評估。用探空新息向量(inno,背景減觀測)和分析殘差(ans,分析減觀測)分析MF試驗和CTL試驗分析場差異。由2018年7月1—31日探空u、v風場、溫度、相對濕度新息向量及分析殘差偏差(Bias)和標準差(Std)廓線圖(圖7)可見,00時冷啟動,MF試驗與CTL試驗新息向量偏差和標準差廓線完全重合,說明兩組試驗初始條件是一致,即初始的背景場(FNL資料)和同化分析資料(包括探空資料)一致。2組試驗分析殘差廓線存在差異,u、v風場偏差和標準差在400~200 hPa相差較大,溫度和濕度的偏差和標準差在低層相對較大。MF試驗風場分析更接近探空觀測,溫度和濕度則稍微偏離探空觀測。從相對濕度偏差和標準差來看,MF試驗比CTL試驗相對濕度更大(濕)。我國從2000年起應用L波段探空儀,采用的是碳濕敏元器件數字式電子探空儀,濕度觀測資料偏干(郝民等,2018;唐南軍等,2014)。因此,MF試驗的濕度分析更接近實況。2組試驗溫度分析差異小,探空溫度觀測在低層與ERA再分析比存在一些正偏差(田偉紅等,2019),因此MF試驗溫度分析更接近ERA再分析。總體而言,MF試驗分析風場和濕度更接近實況。

圖7 2018年7月1—31日冷啟動背景場、分析場與探空觀測的(a1~d1)Bias和(a2~d2)標準差(Std)(a)u風場,(b)v風場,(c)溫度,(d)相對濕度(實線:新息向量,虛線:分析殘差,黑線:CTL試驗,紅線:MF試驗)Fig.7 (a1-d1) Bias and (a2-d2) standard deviations of analyses and forecasts from radiosonde observations for the period from 1 to 31 July 2018 with cold start(a) u-component, (b) v-component, (c) temperature, (d) relative humidity(solid line: innovation, dash line: analysis residual, black line: CTL experiment, red line: MF experiment)

由2018年7月2—31日暖啟動探空u、v風場、溫度及相對濕度資料新息向量及分析殘差廓線(圖8)可見,00時暖啟動,MF試驗與CTL試驗的新息向量和分析殘差偏差和標準差廓線都不重合,說明兩組試驗背景場和分析場都存在差異。MF試驗和CTL試驗每3 h同化一次觀測資料并做3 h預報提供給下時次同化做背景場,暖啟動同化分析的背景場是CMA-MESO模式系統3 h預報,冷啟動時MF試驗與CTL試驗的分析已出現差異(圖7分析殘差),因此,多次同化循環暖啟動時2組試驗的背景場一定存在差異,雖然觀測資料一致,暖啟動時背景場不一致導致新息向量差廓線不重合。風場、溫度場和濕度場新息向量的標準差MF試驗較CTL試驗略有減小,說明暖啟動時MF試驗背景更接近觀測,即循環預報MF試驗更接近觀測。冷啟動和暖啟動2組試驗分析殘差廓線差異類似,400~200 hPa風場絕對誤差相差較大,溫度和濕度的均方差在低層更大。由此可見,暖啟動MF試驗分析的風場和濕度場更接近實況。

圖8 同圖7,但為2018年7月2—31日00時暖啟動背景場、分析場與探空觀測的Bias和標準差(Std)Fig.8 Same as Fig.7, but for the period from 2 to 31 July 2018 with warm start

降水常采用TS(threat score)和Bias(bias score)分析(韋青等,2020)。TS值越大,Bias值越接近1,表示對降水的預報越準確。6 h降水量劃分為小雨[0.1,4) mm,中雨[4,13) mm,大雨[13,25) mm,暴雨[25,60) mm和大暴雨(≥60 mm)5個量級。

從00時冷、暖啟動6 h降水檢驗TS評分(圖9、圖10)來看,00時冷啟動MF試驗對前6 h降水預報影響較大,各量級降水TS評分值明顯提升,隨著預報時效增加,對降水的影響逐漸減弱,TS評分值減小。MF試驗每3 h同化都能更新一次分析場,不斷改進模式的分析和預報效果(圖9顯示冷啟動結果前6 h基本為正貢獻),因而MF試驗暖啟動逐6 h降水TS評分檢驗基本表現為弱的正貢獻(圖10,大部分量級降水TS評分在各時段都有所提高)。6 h 降水檢驗偏差(圖9)則顯示00時冷啟動2組試驗偏差差異明顯,所有量級降水MF試驗的偏差向1減小,而00時暖啟動2組試驗主要是大雨和暴雨的偏差(圖10)有差異,前12 h MF試驗偏差更接近1,而后12 h則是CTL試驗更接近1。結合背景和分析與探空資料的偏差和標準差以及冷暖啟動的降水檢驗結果來看,MF試驗能改進同化分析場和降水預報,但改進的幅度有限,維持的時間不是很長。這也許是背景誤差做3種尺度分離時與模式的物理過程等結合不夠緊密有關,需要進一步深入研究。

圖9 2018年7月1—31日00時冷啟動逐6 h累計降水檢驗的(a~d)TS評分和(e~h)Bias(a,e)00—06時,(b,f)06—12時,(c,g)12—18時,(d,h)18—00時Fig.9 (a-d) TS and (e-h) Bias values of the 6 h accumulated rainfall simulated in the experimental region for the period from 00:00 UTC 1 to 00:00 UTC 31 July 2018 with cold start(a, e) 00:00-06:00 UTC,(b,f) 06:00-12:00 UTC,(c,g) 12:00-18:00 UTC,(d,h) 18:00-00:00 UTC

圖10 同圖9,但為2018年7月2—31日暖啟動結果Fig.10 Same as Fig.9, but for the period from 2 to 31 July 2018 with warm start

4 結論與討論

本文通過對2018年6—8月的背景誤差進行尺度分離,擬合出3種不同尺度的水平相關特征尺度在CMA-MESO 3Dvar系統中應用,與業務測試系統中單一水平特征相關尺度進行個例和連續試驗對比分析,結果顯示:

(1) 二維離散余弦變換尺度分離背景誤差樣本的3種水平特征相關尺度垂直結構相似,水平特征尺度的水平尺度相隔幾十至幾百千米。

(2) 3種水平特征相關尺度和單一水平特征相關尺度的濕度、溫度和風場分析增量差異較小,濕度分析增量差異不超過6%,風場分析增量差異不超過1.2 m·s-1,溫度分析增量差異不超過0.35℃。3種水平特征相關尺度試驗降水預報更接近實況。

(3) 3種水平特征相關尺度在CMA-MESO 3 km 系統應用連續試驗結果顯示,3種擬合的水平特征相關尺度試驗對u風、v風、濕度分析有明顯正影響,分析更接近實況,對溫度分析影響較小;對降水預報有改善,冷啟動預報前6 h TS評分提高明顯,偏差減小向1靠近,暖啟動24 h逐6 h降水預報TS評分值都有提升,Bias差異不大。

猜你喜歡
背景特征水平
張水平作品
“新四化”背景下汽車NVH的發展趨勢
《論持久戰》的寫作背景
當代陜西(2020年14期)2021-01-08 09:30:42
如何表達“特征”
不忠誠的四個特征
當代陜西(2019年10期)2019-06-03 10:12:04
加強上下聯動 提升人大履職水平
人大建設(2019年12期)2019-05-21 02:55:32
抓住特征巧觀察
晚清外語翻譯人才培養的背景
線性代數的應用特征
河南科技(2014年23期)2014-02-27 14:19:15
做到三到位 提升新水平
中國火炬(2010年8期)2010-07-25 11:34:30
主站蜘蛛池模板: 欧美日韩国产成人在线观看| 国产精品粉嫩| 亚洲精品免费网站| 亚洲国产理论片在线播放| 国产精品免费福利久久播放| 小说区 亚洲 自拍 另类| 在线精品欧美日韩| 日本妇乱子伦视频| 欧美丝袜高跟鞋一区二区| 日本高清有码人妻| 99热国产这里只有精品无卡顿"| 国产三区二区| 亚洲 日韩 激情 无码 中出| 无码国产伊人| 欧美成人午夜影院| 一本久道久久综合多人| 在线精品亚洲一区二区古装| 在线观看视频一区二区| 亚洲码一区二区三区| 国产成人福利在线视老湿机| 依依成人精品无v国产| 欧美中文字幕一区| 欧美午夜在线播放| 老熟妇喷水一区二区三区| 国产香蕉97碰碰视频VA碰碰看| 日韩性网站| 亚洲成人动漫在线| 国产91高跟丝袜| 久久人体视频| 亚洲男人的天堂在线观看| 天堂av高清一区二区三区| 亚洲精品波多野结衣| 欧美精品影院| 国产永久在线观看| 国产成人精品一区二区不卡| 成人午夜精品一级毛片| 成年女人a毛片免费视频| 一个色综合久久| 亚洲伊人电影| 亚洲国产无码有码| 黄色污网站在线观看| 性色在线视频精品| 国产精品毛片一区| 日本一区二区不卡视频| 国产精品第页| 巨熟乳波霸若妻中文观看免费| a毛片免费在线观看| 色综合天天综合中文网| 亚洲高清在线播放| 亚洲无线观看| 国产91视频免费观看| 好久久免费视频高清| 日本手机在线视频| 成人综合网址| 亚洲 日韩 激情 无码 中出| 国产成人欧美| 国产XXXX做受性欧美88| 国产在线欧美| 久久久四虎成人永久免费网站| 亚洲高清日韩heyzo| 亚洲美女一区二区三区| 国产综合欧美| 亚洲精品不卡午夜精品| 国产91线观看| 色综合成人| 欧洲亚洲欧美国产日本高清| 久久婷婷六月| 国产不卡一级毛片视频| 欧美色综合网站| 在线精品视频成人网| 国产尤物视频在线| 色综合久久88| 国产va欧美va在线观看| 熟妇丰满人妻| 国产69精品久久久久孕妇大杂乱| 亚洲综合第一区| 中文无码毛片又爽又刺激| 国产成人综合亚洲欧美在| 亚洲精品图区| 国产无套粉嫩白浆| 性激烈欧美三级在线播放| 天堂成人在线视频|