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

基于多周期趨勢分解和兩級融合策略的浪高預測方法

2023-07-29 11:47:24鄭小羅李其超鄧小東
海洋科學進展 2023年3期
關鍵詞:融合模型

鄭小羅,李其超,姜 浩,宋 巍*,鄧小東

(1. 上海海洋大學 信息學院,上海 201306;2. 國家海洋局 東海預報中心,上海 200081)

浪高是全球大氣系統中的一個重要參數,由于海洋環境極其復雜,預測浪高是一項非常具有挑戰性的任務[1]。一般地,浪高的預測方法大致可以分為3 類:數值波浪模型、經典線性時間序列模型和機器學習模型。

數值波浪模型基于能量平衡方程,例如第三代波浪模型WAM(Wave Model)、基于WAM 的近岸波浪模擬SWAN(Simulating Waves Nearshore)和WW3(WAVE WATCH Ⅲ)[2],其利用預報的風場資料,將第三代生成波模式應用于大尺度、長時間的浪高預測。然而,數值波浪模型是物理規律驅動的數值逼近模型,具有輸入復雜度高、對風場資料的要求高及計算過程耗時長等特點,比較適用于大區域海浪的反演和預報。

由數據驅動的方法中,經典線性時間序列模型包括自回歸模型(Autoregressive,AR)、自回歸移動平均模型(Autoregressive Moving Average,ARMA)和差分整合移動平均自回歸模型(Autoregressive Integrated Moving Average,ARIMA)等[3]。Soares 和Ferreira[4]將AR 用于浪高預測;Ge 和Kerrigan[5]對比了AR 和ARMA 用于海浪預報的性能,在10 種不同的海浪數據測試下,ARMA 的誤差和效率均優于AR。Agrawal 和Deo[6]采用ARMA 和ARIMA 在不同預測區間進行在線浪高預測,結果表明ARMA 和ARIMA 的性能非常相似。但是ARMA 模型對序列的平穩性有要求,而ARIMA 能夠通過差分運算將非平穩數據轉換成平穩數據,解決平穩性依賴問題,所以ARIMA 模型是應用更為廣泛的經典線性時間序列模型。

近年來,以大數據為基礎的機器學習算法表現出強大的預測能力。Makarynskyy[7]和Cornegobueno 等[8]利用人工神經網絡(Artificial Neural Networks,ANN)和極限學習機(Extreme Learning Machine,ELM)預測浪高,在準確率和一致性方面都有較好的表現。由于支持向量機(Support Vector Machine,SVM)和支持向量回歸(Support Vector Regression,SVR)具有強大的泛化能力和完整的數學理論,在海浪預測方面也得到廣泛應用[9-11]。Malekmohamadi 等[12]對SVM、ANN、貝葉斯網絡和自適應神經模糊推理系統(Adaptive neural fuzzy Inference System,ANFIS)的浪高預測效果進行了深入研究,結果表明,ANN、ANFIS 和SVM 的預測結果都在可接受范圍,而貝葉斯網絡的預測結果相對不可靠。陸小敏等[13]通過集成棧式自編碼器和XGBoost(eXtreme Gradient Boosting)[14]算法,進行有效波高預測,但是并沒有考慮到數據分布的不平衡和數據集的規模大小。黃心裕等[15]基于Prophet 算法[16]對海南近海波高和周期進行了分析和預測。Prophet 算法是基于時間序列分解(Seasonal、Trend 及Residual)和機器學習的擬合實現的,適用于具有強烈季節性影響的時間序列分析,對數據異常和趨勢變化具有穩健性。深度學習作為機器學習的子領域,極大地促進了計算機視覺和自然語言處理的發展。長短時記憶網絡(Long Short Term Memory,LSTM)是最流行的時間序列預測深度學習網絡之一,被廣泛地用于預測風速、交通、太陽能和股票價格等[17]領域。莫旭濤和李自立[18]將卷積神經網絡(Convolutional Neural Network,CNN)和LSTM 相結合,在中國南海北部灣的浪高預測取得較高的預測精度。顧興健等[19]基于LSTM 建立的海浪環境短期預報模型,在我國海域上的預測結果取得了較好的精度,但是隨著時間長度的推移,誤差會逐漸增大。此外,許多新的深度學習網絡,例如DeepAR[20]也被廣泛應用到時間序列預測相關領域,但是在浪高預測中應用較少。

上述方法嘗試從多個角度進行浪高預測,但未能充分考慮浪高序列的周期性、非平穩性和非線性,一定程度地限制了浪高預測的精度。基于局部加權回歸(Locally Weighted Regression,Loess)的周期趨勢分解(Seasonal and Trend decomposition using Loess,STL)[21]能夠探索復雜的時序數據的規律,適用于時序預測任務,且具有較強的魯棒性。為了實現對多站點未來12 h 的浪高預測,本文針對多源浪高時序數據的多周期性、非平穩性和非線性等特點,提出了一種基于多周期STL 分解和兩級融合策略的浪高預測方法MSTL-WH。首先通過周期圖法估計多源浪高數據集中的主要周期,然后利用STL 將復雜的浪高序列分解為簡單的趨勢項、周期項和余項,最后使用LSTM 結合兩級融合策略以融合長中短時序下的浪高特征并進行浪高預測。應用實測浪高數據開展實驗,并與經典時間序列模型ARIMA、加性時間序列預測模型Prophet、集成學習模型XGBoost、長短時記憶網絡LSTM 和基于概率預測深度學習模型DeepAR。實驗結果表明,MSTL-WH 能夠對不同站點未來12 h 的浪高進行準確預測,綜合表現優于ARIMA、Prophet、XGBoost、LSTM 和DeepAR 方法,各級浪高下MAE 均低于20%,符合業務化運行標準。

1 方 法

MSTL-WH 的流程如圖1 所示。首先對浪高數據集進行預處理,包括缺失值填充、異常值處理等。使用周期圖法提取浪高數據集中的4 個主要周期,以代表不同監測站點的浪高周期,并根據這4 個主要周期對浪高序列進行4 次STL 分解,將復雜的浪高序列分解為12 個簡單分量;然后采用結合兩級融合策略的LSTM 網絡進行特征提取與融合,有效學習不同周期、不同時間海浪信息對未來浪高的影響;最后使用全連接層并結合注意力機制輸出浪高預測值,實現對多站點浪高的精確預測。

圖1 方法流程圖Fig. 1 Flow chart of the proposed method

1.1 基于STL 的多周期-趨勢分解

近岸浪高受多重因素影響,使得浪高難以預測,近岸浪高具有一定的周期性,可以使用STL 將浪高時間序列分解為3 個簡單分量:趨勢分量、周期分量和剩余分量。不同站點浪高時間序列一般具有不同的周期,因此本文利用從數據集中提取的4 個主要周期,對每一條浪高時間序列進行4 次STL 分解,提取其不同周期下的分量,從而提高預測方法的泛化性能,實現對多源浪高時間序列的精準預測。高序列X(n) ,對其進行離散傅里葉變換為X(ejω),該序列的周期圖可以定義為:

浪高數據主要周期提取的方法采用經典的周期圖法,即信號功率譜密度估計方法。對于一個浪

根據式(1)得到的功率譜估值,找出最大功率對應的頻率,并取其倒數得到該序列的主要周期。

STL 是由Cleveland 等[21]提出的一種經典的時間序列分解方法,其優點在于,對帶有異常值的時間序列分解出的分量有更強的魯棒性,本文中,STL 將浪高序列分解為3 個分量的分解表達如下:

式中: Xi為 浪高序列; Ci為趨勢分量; Si為周期分量; Ri為剩余分量;N 為序列長度。STL 是一種迭代方法,其分解過程主要分為內循環和外循環,內循環嵌套在外循環中。內循環用于更新趨勢分量和周期分量,外循環用于計算穩健的權值,通過在下一次內循環中使用這些權值,以減少異常值對更新后續內循環中趨勢分量和周期分量的影響。內循環是STL 的主要部分,具體循環過程如下。

第1 步:去趨勢性。在內循環迭代 k+1次 ,使用在第 k 次 迭代時得到的趨勢分量對原始浪高序列進行去趨勢:。

第2 步:平滑周期性。使用Loess 對去趨勢浪高序列 Xi′的 每個周期子序列進行平滑,得到浪高序列的臨時周期分量。

第4 步:臨時周期分量去趨勢性。將第2 步中得到的臨時周期分量減去防止低頻信號進入周期分量,得到第 k+1次 循環的周期分量,表示為。

第5 步:去周期性。原始浪高序列 Xi減 去周期分量得到去周期序列,表示為。

第2、3、4 步是內循環的周期平滑部分,第6 步是趨勢平滑部分。在外循環中,利用內循環所得的浪高趨勢分量和浪高周期分量來計算浪高剩余分量。浪高序列 4 個主要周期經過STL 分解,被分解為12 個簡單分量,接下來以這12 個分量作為神經網絡預測模型的輸入。

1.2 兩級融合策略

對于不同周期STL 分解的12 個分量,直接利用LSTM 進行特征學習可能會破壞各周期的獨立性,造成多周期之間的干擾。為此,結合兩級融合的思想設計特征學習網絡,網絡第一層使用4 個LSTM 網絡分別提取4 個不同周期輸入子序列的特征。

LSTM 網絡是RNN 的一種改進模型,相較于RNN 模型,它具有特殊的記憶和遺忘模式,每一個LSTM 網絡擁有一個記憶單元(cell),在 t 時刻的狀態記為 ct,其值通過輸入門、遺忘門和輸出門更新,它們一般使用sigmoid 或者tanh 函數進行激活。記憶單元的工作流程為:在時刻t,記憶單元通過3 個門接收當前狀態 xt與 上一時刻記憶單元的隱藏狀態 ht-1,除此之外,每一個門還接收記憶單元的狀態 ct-1。接收到輸入信息之后,每一個門對不同來源的輸入進行計算,并且由其激活函數決定其是否激活。輸入門的輸入經非線性變換后,與經過遺忘門處理過的記憶單元狀態進行組合,產生新的記憶單元狀態 ct。最后,記憶單元狀態 ct通過非線性運算和輸出門的控制產生記憶單元的輸出 ht。

從浪高時間序列的長、中、短周期考慮,結合兩級融合的思想設計預測網絡,網絡第一層使用4 個LSTM 網絡分別提取輸入序列中4 個子序列的特征,然后將輸出融合到一個LSTM 網絡中,之后再連接2 層LSTM 網絡進一步提取特征(LSTM 網絡中的隱藏單元數量均為128)。

1.3 自注意力機制

對浪高特征進行充分提取之后,每個特征已經獲得了一定的權重,為進一步強化特征,提升模型預測精度,使用自注意力模塊對權重進一步分配。

注意力機制本質為一個查詢對應多個鍵值對的映射。注意力機制的原理可分為3 個階段:第1階段計算每一個查詢Q 和各個鍵K 的相似度,獲得每個鍵對應值V 的權重;第2 階段使用Softmax激活函數歸一化權重,從而獲得權重系數;第3 階段對權重系數和對應的鍵值V 進行加權求和,獲得最終注意力權值,其計算式為:

式中: Q ∈Rn×dk; K ∈Rn×dk; V ∈Rn×dk, n為 輸入序列長度, dk為輸入序列特征維度。

自注意力機制也被稱作內部注意力機制,是一種特殊的注意力機制。一般情況下,注意力機制中的輸入和輸出是不一樣的,比如在翻譯任務下,輸入序列是需要翻譯的英文語句,輸出序列則是翻譯后的中文語句,注意力機制發生在輸出的Q 和輸入的所有元素之間。對于自注意力機制則是發生在輸入序列或輸出序列的內部元素之間,也就是式(3)中Q、K 和V 三個矩陣都來源于相同的輸入序列。自注意力機制可以減少對外部信息的依賴,將輸入序列不同位置的信息關聯,擅長提取數據的內部相關性。

2 實驗驗證

實驗驗證了MSTL-WH 在測試集和新數據集上的預測誤差,并通過與ARIMA、XGBoost、Prophet 和LSTM 這4 個具有代表性的預測方法對比,證明本文所提方法的優越性。本文實驗硬件環境如下:處理器為英特爾酷睿i9-10885H,CPU 頻率為2.4 GHz,內存為16 GB;操作系統為Windows10(64 位);程序設計語言為Python3.6.9(64 位),集成開發環境為PyCharm Professional 2019.3.3;程序設計過程中的LSTM 網絡由Keras2.1.5 程序包實現。整個浪高預測實驗主要分為4 部分:①數據集構建;②消融實驗;③對比實驗;④泛化性能驗證。

2.1 評價指標

為驗證MSTL-WH 的有效性,選擇平均絕對百分誤差(Mean Absolute Percentage Error,MAPE)和平均絕對誤差(Mean Absolute Error,MAE)作為評價指標,對比不同預測模型預測浪高的性能。MAPE 和MAE 的公式如下所示:

式中: yi為 真實浪高;為預測浪高; Nt為預測長度。模型在預測浪高時,MAPE 和MAE 的數值越小,代表預測結果越準確。

2.2 數據集構建

本文的浪高實測數據來自國家海洋科學數據中心[22]。由于浪高監測設備異常或網絡異常等原因,實測數據會包含異常值,因此需要對實測數據進行異常值處理,具體處理方式為:對于異常值連續長度不超過3 的浪高序列,使用異常點前后的3 個有效值的均值進行填充;對于異常值連續長度超過3 的浪高序列進行截斷保存。對所有站點的實測數據進行上述處理后,得到94 條長度不一的浪高序列,這些序列用于后續數據集的制作。

數據集中的觀測站點分布在24°~39°N,包含小長山(XCS)、小麥島(XMD)、連云港(LYG)、大陳(DCH)、老虎灘(LHT)、芝罘島(ZFD)、南麂(NJD)七個海洋站的實測數據,每條浪高數據時間跨度為2017 年1 月至2020 年4 月,采樣頻率為1 次/h,共29 160 h。

由于浪高實測數據跨越了不同海區,其浪高數據的周期性變化也大不相同。為了確定對浪高數據進行STL 分解的周期值,使用周期圖法對每一條浪高序列進行功率譜密度估計。提取功率最大的頻率,取其倒數作為候選周期,共獲得94 個候選周期,各觀測站點浪高時序數據的主要周期如圖2 所示。結合候選周期的直方圖,最終確定140、221、314 和460 h 作為周期值。根據上述4 個周期,STL 將每一條復雜的浪高序列分解成12 個簡單分量,將每個分量歸一化至區間(0,1)后用于后續的數據集制作。為了捕捉不同站點的海浪周期,從而更加準確地預測浪高,在制作數據集時,滑動窗口的長度應大于海浪序列數據的周期。因此,在使用滑動窗口切分浪高序列時,窗口大小設置為500 h,步長為2 h,前488 個時間點的浪高值及其分量作為輸入,后12 個時間點的浪高值作為真值,對所有序列進行滑動窗口切分后,共得到69 836 個樣本。按照8∶2 劃分訓練集和測試集,得到55 868 個樣本的訓練集、13 967 個樣本的測試集,數據形狀分別為(55 868,488,12)和(13 967,12),訓練集中各等級浪高樣本的數量如圖3 所示。

圖3 訓練集中各等級浪高的樣本數量Fig. 3 Number of samples at different levels of wave height in the training data set

2.3 消融實驗

首先通過對比實驗,確定一個具有3 層LSTM 的預測網絡,基于此網絡對比采用單周期STL 分解前后和多周期STL 分解的浪高預測誤差,驗證多周期STL 分解的有效性,然后分別在是否采用多周期STL 分解的基礎上,對比采取兩級融合前后的浪高預測誤差,證明兩級融合的有效性,最后通過實驗確定融合層后LSTM 的層數。

為了驗證多周期STL 分解的有效性,向一個3 層LSTM 預測網絡分別輸入原始浪高時間序列、單周期STL 分解后的浪高時間序列和多周期STL 分解后的浪高時間序列,對比浪高預測誤差,其中每層LSTM 的隱藏神經元個數為128。未采用STL 分解時,模型輸入數據為原始浪高序列,維度為(488,1);采用單周期STL 分解時,使用周期圖法得到數據集中功率最大周期,并據此對輸入的浪高序列進行STL 分解,模型的輸入數據維度為(488,3);采用多周期STL 分解時,模型的輸入數據維度為(488,12)。模型的輸出均為未來12 h 的浪高值,結果如表1 所示。

表1 多周期STL 分解對預測性能的影響Table 1 Impact of multi-period STL decomposition on prediction performance

從表1 可以看到,使用單周期進行STL 分解在一定程度上可以減小預測誤差,相比之下,多周期進行STL 分解可以更進一步減小預測誤差,使MAE 降至0.15 m,MAPE 降至30.97%,這表明多周期STL 分解對于減小預測誤差發揮了重要作用。

接下來為了驗證兩級融合的有效性并進一步說明多周期STL 分解的有效性,我們分別對單周期STL 分解及多周期STL 分解下是否采用兩級融合對預測誤差的影響進行了對比。在融合層之后連接1 層LSTM 網絡和1 層全連接層來輸出浪高預測值,其中每層LSTM 的隱藏神經元個數為128。實驗結果如表2 所示。

表2 單/多周期STL 分解下兩級LSTM 特征融合對預測性能的影響Table 2 Impact of two-stage LSTM feature fusion on prediction performance under single/multi-period STL

從表2 可以看出,無論在單周期STL 分解下還是在多周期STL 分解下,采用兩級融合均可明顯減小預測誤差。相對單周期STL 分解,采用多周期STL 分解具有更小的預測損失,結合兩級融合可將MAE 降低至0.11 m,MAPE 降低至20.56%。因此,多周期STL 分解及兩級融合對于多源浪高時間序列預測任務具有明顯提升效果。

最后,在采用多周期STL 分解及兩級融合的基礎上,我們在特征融合層后連接不同層數的LSTM,然后使用自注意力層對權重再分配,最后連接1 個全連接層輸出浪高預測值,其中每層LSTM 的隱藏神經元個數為128,結果如表3所示。

表3 LSTM 層數對預測性能的影響Table 3 Impact of the number of layers of LSTM on prediction performance

由表3 可以看出,隨著網絡中LSTM 層數的增加,網絡的預測性能也隨之提升,但當LSTM層數增加至4 層時,由于網絡參數量過多導致過擬合,網絡性能開始下降,當LSTM 的添加層數為3 時,模型的預測誤差最小,MAE 低至0.07 m,MAPE 低至11.87%,因此,在MSTL-WH 中,使用3 層LSTM 進一步提取融合后的特征。

2.4 對比實驗

基于時間序列的預測模型主要包括經典線性時間序列模型、機器學習模型和深度學習模型,本文從中選取了ARIMA、Prophet、XGBoost、LSTM 和DeepAR 這5 個具有代表性的預測模型與MSTLWH 進行對比。

根據最新發布的海浪等級國家標準[23],將模型預測結果計算平均值后分為4 個等級區間,即2級對應浪高0.10~0.50 m、3 級對應浪高0.50~1.25 m、4 級對應浪高1.25~2.50 m、5 級對應浪高2.50~4.00 m,評價指標結果如表4 所示。同時,為直觀觀察各模型在不同等級浪高下的預測性能,在不同浪高等級范圍隨機選擇了4 個時間序列,使用各種對比模型進行了24 h 預測,預測結果如圖4 所示。

表4 各模型在不同浪高等級下的評價指標結果Table 4 Evaluation metrics of different models at different levels of wave height

圖4 各模型不同等級浪高下的預測效果Fig. 4 Results of different models at different levels of wave height

根據24 h 的浪高值預測結果對比(圖4),MSTL-WH(紅色)與實際浪高值(黑色)的擬合度較高。集成學習方法XGBoost(黃色)和經典的ARIMA 模型(藍色)整體上比其他3 種方法性能更好。DeepAR 模型(粉色)整體性能與ARIMA 模型類似。Prophet 模型(綠色)對低等級的浪高預測值偏高,但是高等級的浪高預測值偏低,說明其對數據變化適應性較差。同樣,LSTM 網絡(灰色)預測網絡效果也不佳,總體預測值偏低,這對海浪預警非常不利。最后,大部分模型在海況4 級的浪高范圍內預測趨勢較好,這可能與這一等級的訓練數據較充分有關,但也存在波動峰值滯后的情況,本文模型滯后性較小。

表4 從誤差指標上進一步給出客觀評價。MSTL-WH 在各等級浪高下的預測指標均優于其他模型,尤其對于3、4 級浪高,即在0.50~1.25 m、1.25~2.50 m 區間的浪高序列,MSTL-WH 的平均絕對誤差僅為0.10 m 和0.17 m,在平均絕對誤差±0.20 m 的誤差范圍內,浪高等級不會發生變化。對于2 級浪高,MSTL-WH 的平均絕對誤差僅為0.03 m,幾乎可以忽略不計。因此,本文提出的MSTLWH 在各等級浪高下有著良好的預測性能,滿足浪高預報業務化的要求(MAPE≤20%)。

2.5 泛化性能驗證

盡管MSTL-WH 在驗證集上的準確性和穩定性表現良好,為了進一步驗證其泛化能力,本文在不同站點以及新數據集上分別進行了實驗。各模型針對不同站點(多源)的MAE 評價指標結果如表5 所示。

表5 各模型在不同站點(多源)的MAE 評價指標結果(m)Table 5 MAE evaluation metrics of different models at different stations (m)

由表5 可以看出,MSTL-WH 在7 個不同站點中有5 站MAE 評價指標上優于其他模型。在DCH站點,MSTL-WH 和XGBoost 的預測指標相同,均優于其他模型。在XMD 站點中,DeepAR 的MAE 評價指標更小,預測性能相對更好,但MSTL-WH 仍然是所有預測模型中排名第二的模型。總體來看,MSTL-WH 在不同站點(多源)的預測誤差變化小,模型性能穩定。

對國家海洋科學數據中心與訓練集一致的7個站點2020 年4 月至5 月的浪高實測數據進行預處理(詳見2.2 節描述)后,共得到3 620 條浪高序列,覆蓋了2~5 級的浪高。MSTL-WH 在這一新數據集上的各項指標如表6 所示。

表6 MSTL-WH 在新數據集上的各項指標Table 6 Results of MSTL-WH in new dataset

由表6 可見,MSTL-WH 在新數據集上各等級浪高的平均相對誤差均小于20%,對比表4 中的性能指標,MSTL-WH 除了在2 級浪高預測有較大誤差外,其余等級浪高情況下最大增加了1.36%的平均絕對誤差。雖然2 級浪高預測的MAPE 為19.16%,但其平均絕對誤差僅為0.05 m,對實際預報業務幾乎沒有影響。由此可見,MSTL-WH 在不同等級浪高預測上具有良好的泛化性。

3 結 語

本文基于浪高實測數據,針對多源浪高時序序列難以預測問題,提出了基于多周期STL 分解和兩級融合策略預測方法。首先使用周期圖法從多源浪高時序數據集中提取主要周期,并據此對浪高序列進行STL 分解,然后使用LSTM 從短、中和長期不同的視野提取浪高特征,采用兩級融合策略進行特征融合,最后使用全連接層輸出未來12 h 的浪高值。

本文所提出的浪高預測方法MSTL-WH 對2、3、4 和5 級浪高序列預測的MAPE 均低于20%,滿足浪高預報業務化運行的要求(MAPE≤20%)。相較于LSTM、ARIMA、Prophet、XGBoost 和DeepAR,MSTL-WH 能更好地預測具有不同周期的浪高序列,同時解決了對不同地區(多源)預測適應性差的問題,降低了預測模型的滯后性,且在新的浪高數據集下,具有良好的泛化性。

在充分分析時間序列數據特征的前提下,本文提出的方法不僅適用于浪高預測,還可以用于其他領域的時間序列預測任務,具有廣闊的應用前景。

猜你喜歡
融合模型
一半模型
一次函數“四融合”
村企黨建聯建融合共贏
今日農業(2021年19期)2022-01-12 06:16:36
融合菜
從創新出發,與高考數列相遇、融合
重要模型『一線三等角』
寬窄融合便攜箱IPFS500
《融合》
現代出版(2020年3期)2020-06-20 07:10:34
重尾非線性自回歸模型自加權M-估計的漸近分布
3D打印中的模型分割與打包
主站蜘蛛池模板: 国产成人亚洲无吗淙合青草| 秋霞国产在线| 亚洲a级毛片| 福利一区在线| 风韵丰满熟妇啪啪区老熟熟女| 尤物精品视频一区二区三区| 尤物特级无码毛片免费| 最近最新中文字幕在线第一页| 毛片免费视频| 伊人AV天堂| 91亚洲视频下载| 伊人91在线| 青青青国产视频| 国产农村1级毛片| 日本爱爱精品一区二区| 孕妇高潮太爽了在线观看免费| 中文字幕有乳无码| 丁香五月婷婷激情基地| 麻豆国产在线观看一区二区| 亚洲Av激情网五月天| 亚洲天堂免费观看| 九九九久久国产精品| 97无码免费人妻超级碰碰碰| 国产麻豆91网在线看| 蝴蝶伊人久久中文娱乐网| 亚洲最大情网站在线观看| 亚洲无线视频| 亚洲欧州色色免费AV| 成年人福利视频| 精品国产91爱| 99视频只有精品| 波多野结衣一区二区三区88| 性欧美在线| 亚洲天堂首页| 无码电影在线观看| 亚洲综合二区| 国产小视频免费| 欧美精品色视频| 国产精品自在在线午夜区app| 亚洲a级毛片| 国产成人AV男人的天堂| 蜜芽国产尤物av尤物在线看| 女人毛片a级大学毛片免费| 成人福利在线观看| 日本妇乱子伦视频| 在线精品欧美日韩| 日本成人在线不卡视频| 亚洲国产清纯| 激情综合五月网| 欧美日韩一区二区三| 日韩国产精品无码一区二区三区| 青青国产视频| 国产区人妖精品人妖精品视频| 欧美无专区| 日本精品影院| 国产在线一区视频| 熟女成人国产精品视频| 91在线国内在线播放老师| 99久久精彩视频| 2021精品国产自在现线看| 91国内视频在线观看| 亚洲热线99精品视频| 国产欧美精品午夜在线播放| 欧美激情一区二区三区成人| 亚洲天堂网视频| 欧美精品一区在线看| 永久免费精品视频| 久久精品国产电影| 国产一级特黄aa级特黄裸毛片| 国产人成午夜免费看| 伊人久久大香线蕉成人综合网| 成人av专区精品无码国产| 欧美午夜在线观看| 色国产视频| 国产精品专区第一页在线观看| 亚洲第一综合天堂另类专| 免费a级毛片18以上观看精品| 999国产精品永久免费视频精品久久| 国产欧美精品一区aⅴ影院| 99re在线免费视频| 成年人视频一区二区| 国产视频大全|