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

基于卷積神經網絡的雷電臨近預警模型*

2021-04-20 00:38:40張燁方馮真禎
氣象 2021年3期
關鍵詞:模型

張燁方 馮真禎 劉 冰

1 福建省災害天氣重點實驗室,福州 350001 2 福建省氣象災害防御技術中心,福州 350001

提 要: 從研究人工智能雷電臨近預警模型的目的出發,以卷積神經網絡模型為基礎,結合多個時間序列的雷達產品(組合反射率、液態水含量、回波頂高)與閃電數據,對雷電臨近預報方法進行基于卷積神經網絡結構的應用,以福建省2017—2018年雷達、閃電數據為樣本完成了模型的訓練與預測研究。訓練結果顯示,15~30 min模型訓練樣本測試集準確率為0.798 5;選取福建省2019年20個雷電過程驗證分析表明,15~30 min模型對動力抬升型雷電過程預警TS評分為0.716,夏季局地熱雷暴預警TS評分為0.694,與常規采用雷達、閃電閾值控制的雷電預警算法相比,準確率有一定的提高,具有一定的實踐意義。

引 言

人工智能是近幾年來發展特別迅速的一門科學技術,普遍被認為是下一場科技革命的生產力代表,2017年8月國務院印發《新一代人工智能發展規劃》,提出了面向2030年我國新一代人工智能發展的戰略目標。逐步推進氣象業務工作與人工智能技術的對接,拓展未發掘的預報模式,滿足多角度氣象專業服務需求,是實現氣象工作自動化、智能化、現代化的重要途徑。

我國在20世紀90年代初已開始氣象領域人工智能方面的預報研究,如中國氣象局在“八五”期間就成立了《人工智能、模式識別在天氣預報業務中的應用》課題項目,運用神經網絡模式對暴雨預報進行試驗(王寶榮,1993)。進入21世紀后,國內出現更多基于機器學習的氣象人工智能研究,如賀佳佳等(2017)采用支持向量機(SVM)開展局部短時臨近降雨預測研究;孫全德等(2019)對數值天氣預報模式ECMWF預測的華北地區近地面10 m 風速進行訂正,指出三種機器學習算法的訂正效果均好于傳統訂正方法模式輸出統計〈MOS〉的結果;熊亞軍等(2015)利用KNN(KNearest Neighbor)數據挖掘算法構建等級預報分類器,開展霾等級客觀識別實驗;陳勇偉等(2013)、趙旭寰等(2009)使用BP(back propagation)神經網絡模型,選用7個對流參數對雷暴活動進行潛勢預報。

在雷電臨近預警領域,目前主要是使用雷達組合反射率、液態水含量、回波頂高、閃電數據進行預警,如:張其林等(2010)對閃電發生進行高密度區域識別并外推,得出雷電預警產品;呂偉濤等(2009)借鑒美國NCAR的TITAN(Thunderstorm Identification,Tracking,Analysis,and Nowcasting)算法對雷達數據進行重點區域識別與外推,研究格點雷電發生概率;秦微等(2016)在深圳地區采用TITAN算法對某次強風暴過程中的強回波區進行分析,得出該算法識別效果理想的結果;劉維成等(2015)對甘肅中部雷達回波單體與雷電活動之間的對應關系進行分析,得出回波強度、回波頂高、垂直累積液態含水量與雷電發生的閾值關系,并根據該結論建立了雷電預警方案;常越等(2010)對湖南省閃電的發生和雷達關系進行研究,認為閃電發生與回波強度、回波頂高、速度場特征、垂直液態含水量等有著密切的關系;此外還有研究使用衛星、探空、大氣電場等數據對雷電預警關系進行分析。國外在雷電臨近預警研究方面,Hondl and Eilts(1994)根據雷暴發展到不同高度層的回波特征與閾值進行短時臨近預報,Schmeits et al(2008)運用數理統計方法處理和開發雷電臨近預報系統,Jacobson et al(2011)統計分析衛星、閃電定位系統等提供的雷電數據,計算雷電發生概率。近幾年國外在雷電臨近預報領域也有了一些新的成果,如Ivanova(2019)應用地面微波輻射計進行雷電預警;Baldini et al(2018)利用雙偏振雷達識別graupel粒子,研究與雷電活動的相關性。在人工智能的氣象應用上,還沒有明確應用于雷電臨近預警的成果,但在氣象預報的其他領域已經有較多的人工智能應用產品,如:Andreev et al(2019)將卷積神經網絡應用于云識別與檢測,Saha et al(2016)采用自編碼器對印度季風進行預測應用;Geng et al(2019)引入ConvLSTM建立包含不同編碼器、預測解碼器的LightNet模型對雷電進行預測等。探索一種新的雷電預警思路與當前人工智能技術的興起不謀而合,利用人工智能技術來進行雷電臨近預警具有很好的探索意義。

現行的雷電臨近預警技術方案主要以線性或簡單的非線性函數為主,非線性化程度低,以非變形外推為主,較少看到變形外推的應用模型,或者以中尺度分析技術中的“配料法”為基礎,對雷達、閃電數據進行“配料”閾值控制預警,而實際雷電的發生與各影響要素之間不僅只存在線性的關系,可能還有我們未發掘的更為復雜的非線性規則,因此采用具有強非線性表征能力的神經網絡技術來處理雷電臨近預警問題是提升預警效果的有效途徑。

本文選用卷積神經網絡(convolutional neural network,CNN)為基礎,構建基于CNN的雷電臨近預警模型,將雷達組合反射率、液態水含量、回波頂高、云地閃電數據等指標按CNN所需求的“圖片”模式進行網格“切片”化處理,對數據進行時間、空間維度上的延伸,建立相應的神經網絡模型,采用谷歌公司人工智能線性代數編譯器(TENSORFLOW)為計算工具,對該模型進行神經網絡訓練,將得到的結果進行實例的檢驗和應用。

1 模型設計

1.1 CNN概述

CNN是一種專門用來處理類似網格結構數據的神經網絡,它采用具有一定大小、包含權重矩陣的格點網(卷積核,也稱過濾器)對輸入圖像逐個位置進行掃描和卷積,利用卷積實變函數運算實現對圖像局部特征的信息提取并保存在權重矩陣中,加上中間添加的非線性激活函數,使整個神經網絡具有很強的學習和非線性處理能力,實踐表明該計算方式可以很好地獲取圖片的總體與局部特征,其一般流程是對輸入的多通道圖像(可以把氣象格點數據看成是圖像)經過多次的卷積、池化操作,提取到圖像深層次的高維特征,再將卷積、池化操作后的矩陣展開成神經網絡層,經過一定層數的神經網絡層計算后,根據需求輸出對輸入圖像某一種類別屬性的預測結論,例如判斷輸入的圖像是哪一種動物、哪一個數字等。

盡管CNN已經在計算機視覺、圖像分析等領域取得了積極的成果,但也存在著一定的不足和缺點,主要表現在對訓練數據的需求量大、自身模型計算量龐大及訓練時間長、硬件要求高、生物學基礎支持不足、無法掌握其內部的計算機理、過多的參數設計使得模型設計的最優解難以精確得到,等等。尤其在氣象預報應用領域,由于常規的CNN處理的是圖片輸入、有限類別輸出的情形,而在格點化氣象預報領域,則要求輸入圖片,針對圖片的所有格點都要求有結果輸出的情況,即圖片輸入、圖片輸出,在人工智能領域稱這樣的問題為像素級的預測模型,因此常規的CNN模型無法直接應用到氣象格點預報中。為了解決氣象預報這樣像素級的人工智能CNN模型,本文提出了一種采用切片式的數據處理方法,將輸入的氣象格點圖片根據像素預報點轉換成多量的“小切片”,作為新的輸入圖片進行CNN預測;為了規避這種方法可能帶來的大計算量,本文設計了適宜的數據壓縮模型,有效降低了整個模型的計算量,使模型的預測計算可以在1~2 min內實現福建省區域內的格點化預警預報,保證預警產品的時效性,具體模型內容討論如下。

1.2 參數選取

深度神經網絡是一個封閉的黑盒子計算單元,目前在人工智能領域對于神經網絡的工作原理尚不清楚,只知道“神經元+激活函數”可以實現復雜的非線性數學變換,本文根據目前雷電臨近預警預報工作常用的思路與方法,結合雷電發生、發展的規律以及預報產品的時效、精度需求,選用下面幾個參量作為神經網絡的輸入與輸出單元。

(1)雷達產品數據,包括組合反射率、液態水含量、回波頂高三個參數。 這是目前雷電臨近預警預報工作最常用的三個指標。國內多項研究也表明這三個參數可以很好地適配雷電臨近預警工作,結合雷達數據質量控制與SWAN雷達拼圖的處理等相關研究,如對有源干擾回波識別(黃小玉等,2019)、雷達CAPPI格點數據均一性處理(葉飛等,2020)、SWAN雷達拼圖產品處理(張勇等,2019)。對所有雷達產品數據進行質量控制和處理,剔除可能存在誤差或錯誤的雷達產品數據。

(2)云地閃數據。 閃電定位監測數據是目前表征雷電發生位置、時間最直接、精度最高的觀測產品。利用不同時次地閃發生位置的變化來預測下一時次的雷電發生位置,是當前地閃數據在雷電臨近預警工作中的主要應用思路。本文通過對云地閃數據按雷達時次(本文規定每個小時的第00、06、12、…、54 min,即雷達數據產生的時間為雷達時次,以下同)和設計的柵格大小進行每個柵格云地閃頻數的統計(時間分辨率:6 min,空間分辨率:0.01°),實現將一條一條的云地閃數據轉變成和雷達產品數據一樣格式的格點化數據,應用到神經網絡模型中,設待計算時刻為t1,計算t1往前6 min的時刻t0,獲取[t0,t1]時間、柵格區域內所有發生的云地閃數據,計算其總頻數作為該t1時刻、該柵格區域的標示值。

(3)預警輸出標簽設計。按二分類對預警輸出結果進行標簽編制,考慮到雷達、閃電數據的實際到達時間存在一定的延遲(SWAN拼圖要等到所有單站雷達數據全到位才能拼圖,因此延遲在10 min左右),加上模型計算、結果傳輸所消耗的時間,雖然也可以做到0~15 min的預警,但該預警結果基本失去其時效意義,因此預報結果設定為當前時刻之后15~30、30~45、45~60 min內雷電發生位置的網格化產品,限于篇幅關系,本文僅對15~30 min 內的預報模型進行討論、介紹。具體規則:對于某一個預警格點,當該格點在15~30 min的時間區間里監測到有一次以上的云地閃數據時,即標記該預警格點標簽值為1,否則則記為0。

1.3 數據切片處理

針對氣象預報對預報結果的需求特點(大輸入、大輸出),本文把每個預報格點看成是一個實例,把一張300×300像素點的預報圖片,看成是預報300×300個實例的樣式;對每一個預報格點,分別向四個方向拓展指定的網格,形成以預報格點為中心,具有一定長、寬的網格切片,對每個時次、每個參數的網格圖層進行一次切片,得到指定大小的網格矩陣,類比為圖片的一個“通道”,這樣每個時次的每個參數就可以得到一個通道,n個參數、m個時次就可以得到n×m個通道的輸入層,再經過CNN的訓練和計算后,就可以得到該網格(像素點)的預報結果了,具體思路如圖1所示。

圖1 單個格點的切片處理方法及雷電臨近預警CNN模型示意圖Fig.1 Slice processing method for single grid points and lightning nowcasting warning CNN model diagram

1.4 數據壓縮

在上述采用切片提取與該網格預報結果有關參數的矩陣時,通過向四個方向分別拓展一定距離的網格來形成該參數一個“通道”的“圖片”。經驗表明,福建省內中小尺度天氣系統移動速度最大在0.1°·6 min-1以內,也就是說,如果要做15~30 min預報的話,需要拓展網格的距離在0.25°~0.5°,保證對中心點有影響的數據包含在這個“圖片”內,換算成0.01°的個數就是25~50個柵格,取40格來計算的話,那么按上述拓展方法得到的一個時次、一個指標的切片通道“圖片”的分辨率是(40×2+1,40×2+1),即(81,81)。由于這個圖片的大小比較大,為了降低計算量、提高預警速率,需要對原始網格化圖片進行壓縮、池化處理。原來網格單元為0.01°×0.01°,選取k×k個相鄰的網格,將這k×k個網格的所有數值提取出來作為一個樣本序列,取這個樣本序列的最大值、中位數表征這個網格壓縮后的信息。以本文所討論15~30 min為例,對于15~30 min 預報模型,選取k=3,即把某參數、某通道切片的每9個網格壓縮成1個網格,每個網格儲存原9個網格數值的最大值、中位數兩個數值信息,如圖2所示,對雷達組合反射率、液態水含量、回波頂高提取這兩個值,對于地閃參數則僅用這3×3個網格的總地閃個數表示,即實現了2/9的壓縮、提取,大大減少了計算量,滿足了臨近預報對預報產品時效性的需求。

圖2 某參數某通道數據3×3壓縮方法示意圖Fig.2 A diagram of a channel data 3×3 compression method for a parameter

本文實際模型中先對各指標按3×3網格進行格點壓縮,對于每個預報格點按東10、西13、南8、北6個格點距離拓展進行切片處理,這樣處理后每個參數通道的切片分辨率為(10+13+1,8+6+1),即24×15=360,對每個參數選取當前預報時刻往前2個雷達時次的值,最終本文設計模型的每個預報格點輸入情況總結如下:

(1)預報格點基準大小為0.03°×0.03°;

(2)共包含雷達組合反射率、液態水含量、回波頂高以及云地閃4個參量;

(3)每個參量按0.03°×0.03°網格進行格點壓縮,壓縮后雷達組合反射率、液態水含量、回波頂高包含最大值、中位數值兩個圖片通道,云地閃包含總頻次一個圖片通道,4個參量合計有3×2+1×1=7個圖片通道;

(4)每個參量包含當前時刻以及當前時刻往前6 min時刻的2個時間序列,即按上述處理后,共計有7×2=14個0.03°×0.03°大小的圖片通道;

(5)每個待預報格點在每個通道上按東10、西13、南8、北6個格點距離拓展進行切片處理,每個切片大小為24×15,14個通道得到14個通道的24×15分辨率的小圖片,實際輸入大小為24×15×14=5040。

1.5 CNN模型設計

按上述方法對雷達、閃電數據進行切片及切片壓縮處理后,對于每一個預報格點,可以得到14個通道、24×15分辨率的小圖片,將這14個通道的圖片作為神經網絡模型的輸入,按“兩次卷積-池化-兩次卷積-池化-展開-全連接層-全連接層-激活函數(sigmoid)激活輸出”的流程建立深度學習模型,具體如圖3所示,模型中各個處理單元設計如下:

圖3 CNN設計模型示意圖Fig.3 Diagram of the model design of convolutional neural network

(1)卷積層:選用2×2濾波器,移動步長為1,采用邊界不填0的卷積處理方式。這主要是考慮使用的4個參數的數值存在較多0值的情況,如果采用0值填邊的話會使神經網絡在訓練時以為邊界沒有天氣活動,進而引起較大的誤差。通道數按輸入通道數的1.25倍數遞增(取整),采用層標準化技術(batch normalization)對每層輸入值進行規范化處理,最后按激活函數(relu)對卷積層做非線性激活,隨意節點丟棄率(dropout)取值為0.15,防止模型的過擬合,將結果傳遞到下一層模型中。

(2)池化層:按最大池化(max pooling)模式進行池化操作,移動步長為2,采用邊界不填0的卷積處理方式,理由同卷積層的設計。

(3)全連接層:在兩次卷積-池化操作后,將輸出的多通道圖片按單列格點展開,進行全連接操作。全連接共兩層,全連接節點設置為64、32,每個全連接層都采用層標準化技術對輸入值進行規范化處理,選用激活函數(relu),隨意節點丟棄率(dropout)取值為0.15,防止模型的過擬合,最后使用激活函數(sigmoid)對前面神經網絡運算結果進行最后預測,定義最后的預測值>0.5為Y=1類標簽,預測值≤0.5為Y=0類標簽,得到0、1的二分類結果。

對于Y標簽按如下規則確定:以本文討論的15~30 min預報為例,若預報網格在預報時刻往后的15~30 min監測到云地閃,則不管該時段內發生多少次閃電,都標記為1,沒有監測到云地閃即標記為0。

1.6 數據與模型訓練

以2017年3—10月、2018年3—7月福建省雷達、閃電監測數據為樣本,對每個月選取當月閃電定位監測數據文本文件大小最大的兩個日期,對每個選取日期按每個小時、每個雷達時次、每個經緯度網格進行切片選取(X輸入)與對應15~30 min后該經緯度網格雷電發生情況(Y二分類標簽)進行提取。 由于Y=0且X大部分數據為0值(即各個網格的雷達、閃電數據基本為0,對應預報時效內沒有閃電)的情況占所遴選樣本的比例非常大,因此在數據提取時判定當X數據中雷達組合反射率的值大于20 dBz的個數占比小于10%時,本次樣本提取放棄; 此外由于閃電定位監測數據存在一定的偏差,經常出現在同一時間段,某個閃電數據的發生位置與當時的閃電集中區域偏離較大的情況,為了使訓練樣本Y=1的標簽降低噪聲、干擾的情況,樣本提取時剔除閃電定位方法為2站、雷電流幅值絕對值≤2 kA、某個閃電前后10 min內距離本次閃電0.5°范圍內沒有其他閃電監測記錄的閃電數據。為了避免訓練樣本出現0、1標簽個數差別太大的情況,對當前得到的所有訓練樣本進行隨機遴選,根據計算機的計算負荷,設置總量為80 000個樣本,遴選0、1樣本個數各40 000個,進入后續的網絡訓練。 具體網絡訓練采用Adam優化器,即在梯度下降進行迭代求最優解的基礎上,根據一階梯度矩估計、二階梯度矩估計計算更新步長的優化模型,設置初始學習率為0.001,模型學習率自動衰減為0.9,設置單次訓練樣本數(batch size)為120,選取10%的樣本容量作為測試數據,選擇二分類損失函數(binary crossentropy)為損失計算方式,設置最大訓練次數為500次,將訓練過程中準確率最好的參數矩陣保存為最終訓練結果。

按上述方法與規則對模型進行訓練,繪制整個訓練過程的損失(losses)與準確率(accuracy)曲線(圖4)。由圖4可知,損失曲線下降的速度適均,說明模型設計的初始學習率可行;最終測試集損失控制在0.430 1,測試集準確率為0.798 5,測試集與訓練集的損失、準確率曲線基本貼合,訓練后期也沒有出現明顯的分開,表明模型沒有出現明顯的過擬合現象,也說明按本文設計的訓練數據提取、遴選方法可行,保證了數據集數據質量的一致性與噪聲的隨機分布。整個訓練模型在大概訓練350個訓練次數后開始趨于平穩,模型學習率經自動衰減100個訓練次數仍沒有明顯下降和變化,可見模型在這時基本已訓練到飽和,整個模型的訓練過程合理、可行。

圖4 模型訓練損失曲線(a)和準確率曲線(b)Fig.4 Curves of model training losses (a) and accuracy (b)

2 應用與檢驗

本文所建立與訓練的模型在福建省2019年的汛期工作(3—10月)中得到完整的應用。從總體上看本文討論的15~30 min雷電預警效果總體表現穩定和良好,取2019年4月22日14時的預報、實況圖(圖5)為例,由圖可見,模型預報的幾個雷電發生區域在15~30 min后的都有閃電發生。需要說明的是,預報圖中對預報格點按福建省邊界外圍一定距離做了裁剪,因此實況圖中左下角的閃電并非沒有預報出來,而是被裁剪。

圖5 2019年4月22日14時的15~30 min雷電預報(a)、實況(b)對比Fig.5 The 15-30 min lightning nowcasting (a) and observation (b) at 14:00 BT 22 April 2019

為了更精確了解本模型的實際準確率情況,本文按系統動力抬升型雷電過程、副高邊緣局地熱雷暴過程這兩種福建省主要雷電天氣類別,選取了2019年主要雷電天氣(兩種類型各選取10 d)作為樣本進行TS評分計算,結果表明:

動力抬升型雷電過程, TS平均評分為0.716,平均漏報率、空報率分別為0.095、0.190;

局地熱雷暴型雷電過程,TS平均評分為0.694,平均漏報率、空報率分別為0.112、0.194。

檢驗結果表明,模型對前汛期動力抬升型的雷電過程預警效果略微好于夏季局地熱雷暴的雷電過程,考慮到閃電定位系統本身存在的誤差,且福建省多山地地形,這樣的預報結果基本可以接受。有研究采用雷達與閃電閾值控制算法,對福建省進行格點化雷電臨近預警應用,按雷達組合反射率37 dBz、垂直液態水含量1.5 kg·m-2、回波頂高11 km作為雷達預警控制的閾值,研究中直接閃電預警判定距離為10 km,間接閃電預警判定距離為15 km,閃電預警解除判定距離為15 km,最后分析得出其TS評分為0.671(張燁方等,2019)。相比之下,本文設計與訓練的CNN模型的準確率有了一定提升。

3 結論與展望

本文以雷達產品(組合反射率、液態水含量、回波頂高)和閃電數據為指標,對CNN模型進行了適應雷電臨近預警需求的改進和調整,對改進后的模型進行了訓練、檢驗,得出以下結論:

(1)本文設計的切片式處理方法可以解決CNN在氣象格點預報應用中遇到的像素級分類問題,設計的格點壓縮算法可以在短時間、一般計算硬件的條件下,得到一個省范圍內的格點預報產品,滿足臨近預警服務對時效性的要求。檢驗表明設計的模型在預警效果上比常規采用“配料法-閾值控制”的雷電臨近預警模型有了一定的提高,模型對動力抬升型雷電過程的預警效果稍微好于副熱帶高壓邊緣局地熱雷暴過程的預警效果。

(2)原始訓練數據的質量對模型的訓練效果有很大影響,原始數據噪聲的大小、不同類別標簽的數量比例遴選不合理甚至可能會使訓練的模型矩陣出現預報結果全是0或1的“垃圾輸出”;此外,在雷電臨近預警領域,對閃電定位數據進行適宜的修訂可以有效提高模型的訓練準確率,本文按“剔除2站定位、雷電流幅值絕對值≤2 kA、某個閃電前后10 min 內距離本次閃電0.5°范圍沒有其他閃電監測記錄的閃電數據”的條件對訓練的閃電數據進行了剔除和修訂后,模型的訓練準確率提高了近0.10,是模型訓練過程中各個環節、參數調整中準確率變化最大的地方。

(3)1~2 h預報思路展望。相較傳統的雷達回波相關性跟蹤(TREC)的外推算法,本文所建立基于切片式處理、CNN模型的雷電臨近預警模型在預報區域的外推上,初步具有變形外推的效果,但在不少方面仍有改進和不足的地方。如:模型沒有根據不同的天氣系統進行分類訓練、預警;模型在參數的選擇上沒有把高程、土壤、云閃及潛勢預報的因子列入其中;限于計算硬件設備的影響,本文輸入數據的時間序列時段、切片大小有限等。需要看到的是,單純依靠上述模型只能滿足15~60 min的雷電預警需求,要做1~2 h的預警產品,需要再引入其他強天氣潛勢預報涉及的物理量,但這些物理量在格點精度、時間序列點上與雷達、閃電的數據無法同步,這也是影響1~2 h臨近預報格點產品效果的重要因素,本文研究項目也正嘗試對這些物理量進行基于深度學習網絡模型的處理,努力將其融入到本文所研究模型的1~2 h預報產品中,以更好地實現人工智能技術在雷電臨近預警工作中的有效應用。

猜你喜歡
模型
一半模型
一種去中心化的域名服務本地化模型
適用于BDS-3 PPP的隨機模型
提煉模型 突破難點
函數模型及應用
p150Glued在帕金森病模型中的表達及分布
函數模型及應用
重要模型『一線三等角』
重尾非線性自回歸模型自加權M-估計的漸近分布
3D打印中的模型分割與打包
主站蜘蛛池模板: 国产麻豆91网在线看| 色老头综合网| 美女免费黄网站| 伊人91视频| 日韩国产综合精选| 亚洲精品无码人妻无码| 亚洲IV视频免费在线光看| 国产成人高清精品免费软件 | 亚洲婷婷丁香| 国产精品2| 国产亚洲精品无码专| 国产一级精品毛片基地| 国产成人艳妇AA视频在线| 精品少妇三级亚洲| 亚洲成年人片| 欧美性色综合网| 99视频在线观看免费| 精品国产成人国产在线| 狂欢视频在线观看不卡| 亚洲一区二区三区国产精华液| 天天综合网亚洲网站| 亚洲精品老司机| 国产欧美日韩视频一区二区三区| 亚洲熟妇AV日韩熟妇在线| 亚洲系列无码专区偷窥无码| 好紧太爽了视频免费无码| 国产永久免费视频m3u8| 国产亚洲精品97AA片在线播放| 无码精品福利一区二区三区| 18禁色诱爆乳网站| 国产va视频| 国产精品第一区| 午夜福利视频一区| 亚洲自偷自拍另类小说| 久久久久夜色精品波多野结衣| 视频二区国产精品职场同事| 国产精品任我爽爆在线播放6080 | 亚洲不卡无码av中文字幕| 亚洲第一av网站| 亚洲不卡无码av中文字幕| 亚洲欧美精品一中文字幕| 成年女人a毛片免费视频| 中文字幕2区| 99人体免费视频| 人人91人人澡人人妻人人爽| 久久久久久久久18禁秘| 国产精品所毛片视频| 欧美一级特黄aaaaaa在线看片| 久久综合丝袜长腿丝袜| 欧美精品不卡| 国产精品爽爽va在线无码观看| 尤物特级无码毛片免费| 亚洲人成在线精品| 91精品专区| av天堂最新版在线| 亚洲男人天堂网址| 国产日韩欧美黄色片免费观看| 日韩免费毛片视频| 国产精品人莉莉成在线播放| 日本黄色a视频| 久久一级电影| 日韩麻豆小视频| 成人免费午夜视频| 自拍中文字幕| 久久夜色精品| 国产福利2021最新在线观看| 亚欧乱色视频网站大全| a毛片基地免费大全| 欧美a在线视频| 操美女免费网站| 亚洲日韩在线满18点击进入| 999精品视频在线| 美女一级免费毛片| 亚洲精品成人片在线播放| 中文字幕在线欧美| 国产欧美日韩18| 欧美精品H在线播放| 精品成人一区二区| 亚洲一区二区三区国产精品| 久久99精品久久久久纯品| 亚洲免费毛片| 亚洲欧美国产高清va在线播放|