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

基于非線性多尺度模型的黃河三角洲降水量預測

2021-09-02 08:23:36穆玉珠王燕鵬張鳳華喬冬梅
灌溉排水學報 2021年8期
關鍵詞:模型

穆玉珠,王燕鵬,張鳳華,喬冬梅

基于非線性多尺度模型的黃河三角洲降水量預測

穆玉珠1,王燕鵬1,張鳳華1,喬冬梅2*

(1.河南省新鄉(xiāng)水文水資源勘測局,河南 新鄉(xiāng) 453000;2.中國農業(yè)科學院 農田灌溉研究所,河南 新鄉(xiāng) 453002)

【】提高降水量的預測精度,反映降水量的實際特征。基于經驗模態(tài)分解對非線性時間序列的分析和處理的優(yōu)勢,對黃河三角洲氣象站點1954—2018年連續(xù)65a月均降水量數據進行經驗模態(tài)分解(Empirical Mode Decomposition, EMD),得到了系列本征模態(tài)函數(Intrinsic Mode Function,IMF),對IMF進行Hilbert變換,在此基礎上建立了2種黃河三角洲降水量多尺度預報模型。黃河三角洲降水量存在著9、13、23、76、135月左右的周期,并以9個月的波動為主;65 a月均降水量數據預測結果顯示:模型一的相對誤差在0.9%~9.8%之間,模型二的相對誤差在1.6%~11.8%之間,在建模時不考慮初相位的模型一平均預測誤差為2.70%,整體預測精度要優(yōu)于考慮初相位的模型二。2種模型的擬合精度及顯著性均符合要求。

降水量;時間序列;多尺度;EMD;預測

0 引言

【研究意義】在氣候變化和人類活動雙重因素的影響下,各地區(qū)的降水量也發(fā)生了變化[1-3],降水與人類生產、生活及生態(tài)息息相關,降水量的變化關系到區(qū)域水資源的可持續(xù)利用、生態(tài)環(huán)境的保護與經濟社會的發(fā)展,降水量的變化特征與演變趨勢研究已成為水資源領域研究的熱點,相關學者和研究人員都非常關注降水量的精確預測[3-4]?!狙芯窟M展】國內外針對降水量的預測主要為概率統(tǒng)計法和時間序列分析法[5-7]。付明明[8]運用ARIMA模型對新疆喀什地區(qū)的降水量進行預測,利用喀什地區(qū)2個子流域1950—2015年實測年降水量數據分析其模型的適用性和預測精度,模擬降水量和實測降水量之間的誤差相對值低于20%。降雨時間序列往往隨著時間呈現一定的周期性趨勢,這種現象受多種隨機因素的影響和制約,導致周期性趨勢的偏離、跳躍和擺動,使得降水量存在一定的隨機性和非平穩(wěn)性。楊沛羽等[9]基于超閾值(POT)的極端降水事件抽樣、變異分析、趨勢分析等方法對黃河流域77個氣象站點1957—2014年日降水數據進行分析,結果表明黃河流域極端降水量存在顯著的非平穩(wěn)性。概率統(tǒng)計無法適應數據的隨機性和不平穩(wěn)性,傳統(tǒng)的時間序列法常常對數據自身的平穩(wěn)性要求較高。也有部分學者采用多種方法對降水量非線性時間序列進行建模和預測[10-11],何慧等[12]應用PCA-BP神經網絡模型、遺傳算法優(yōu)化神經網絡、RBF神經網絡預測模型、小波神經網絡模型、粒子群神經網絡模型等對廣西月降水量預測,得到較好的應用效果。但降水量時間序列是非線性、非穩(wěn)定的和多尺度的,這些變化有月、年和代等多種,其演變特征是由很多錯綜復雜因素共同作用的結果,上述方法并未在本質上降低時間序列的非平穩(wěn)性。20世紀90年代末,Huang等[13]、任博等[14]提出了一種時間序列分析方法經驗模態(tài)分解(Empirical Mode Decomposition,EMD),該方法是一種基于時間序列局部特征的分解方法,EMD適用于對非線性、非穩(wěn)定過程的復雜數據進行線性化和平穩(wěn)化處理,其將復雜的時間序列分解為有限的、少量的固有模態(tài)函數(Intrinsic Mode Function,IMF)。通過對這些IMF進行Hilbert變換,得到各IMF的時變振幅和頻率,能夠準確地用波內調制機制反映出系統(tǒng)的非線性特性,這也是經驗模態(tài)分解的另一優(yōu)點,因此,可以有效提高降水量的預測精度。【切入點】目前主流的研究方法多為單一尺度的,無法反映降水量實際特征。小波分析雖然作為常用的多尺度分析方法,但是它的不足之處是存在較多噪點干擾,數據分辨率一般,缺少從降低時間序列非平穩(wěn)性和多尺度層面的降水量預測方法研究?!緮M解決的關鍵問題】研究通過對黃河三角洲1954-2018年月降水量的EMD分解,得到IMF固有模態(tài)函數,并對IMF進行Hilbert變換分析,得到相應的瞬時頻率和振幅,分別建立針對研究區(qū)域的是否考慮初相位的2種降水量的多尺度數學模型,為黃河流域降水量的定量研究、水資源的利用及河口治理提供科學依據。

1 材料與方法

1.1 研究區(qū)概況與資料來源

黃河三角洲是指黃河入??跀y帶泥沙在渤海凹陷處沉積形成的沖積平原,該區(qū)域以山東省東營市利津縣為頂點,北鄰吐海河口,南鄰小青河口,面積5 450 km2,其中東營市境內面積5 200 km2。黃河三角洲屬暖溫帶半濕潤大陸性季風氣候區(qū),冬寒夏熱,四季分明,多年平均氣溫為11.7~12.6 ℃,年平均日照時間為2 590~2 830 h,無霜期211 d,年均降水量530~630 mm。三角洲頂部和中部海拔10 m以下的土壤脫鹽較好,為農業(yè)主產區(qū),海拔4 m以下沿海低地地下水位高,土壤鹽堿化強,大部分仍是荒地。

研究所需資料來源于黃河三角洲氣象站點1954-2018年連續(xù)65 a實測資料,在降水量數據處理時,部分年份缺測資料利用相鄰站點采用相關分析進行插補。

1.2 研究方法

1.2.1 經驗模態(tài)分解

EMD算法適合于分析非線性平穩(wěn)信號,可以根據信號本身的特點提取信號的本征模函數,該算法不同于傳統(tǒng)的基于線性穩(wěn)定假設的傅立葉變換和小波變換,是一種新型的自適應信號時頻分析方法。

EMD算法分解的IMF分量應滿足2個條件:第一個條件是這些分量的最大值、最小值和零交叉數必須相等或最多有一個差;第二個條件是由最大值和最小值確定的上下包絡的平均值在任何時間都為0。分解步驟如下:

1)尋找時間序列()的局部極大值和極小值,極大值形成的上包絡為()、極小值形成的下包絡為(),上下包絡的均值()=[()+()]/2。

2)設1()=()-(),檢驗1()是否滿足IMF的2個條件,如果滿足1()就是IMF,如果不滿足,對1()重復1),直至得到k()滿足IMF的2個條件,k()就是第一個IMF分量,1()。記序列的剩余部分為1()=()-1()。

3)對1()重復1)和2),直至剩余部分為一單調序列,分解結束。通常情況下,IMF停止的標準是前后2個()的標準差,即:

當為0.2~0.3時IMF有較好的穩(wěn)定性,同時可以使IMF的物理意義更加清晰。

1.2.2 Hilbert變換

對所有分解出的IMF進行Hilbert變換,C()的Hilbert定義如下:

式中:是柯西主值,根據(2)式,()與()形成復共軛,得到一個解析信號():

()為瞬時振幅,()為相位,即:

相應的瞬時頻率為:

1.2.3 多尺度預測模型

1)模型一

根據EMD和Hilbert變換的結果,不考慮各IMF初相位的情況下建立如下模型:

2)模型二

考慮IMF的初相位,建立如下模型:

式中:表示各IMF分量的初相位,其余參數與模型一的意義相同。

2 結果與分析

2.1 EMD分解結果

對黃河三角洲氣象站點1954-2018年連續(xù)65 a月平均降水量數據,進行EMD分解,得出6個IMF分量和1個殘余分量(圖1),每個IMF分量代表1個尺度上的波動變化特性。由圖1可以看出,IMF1表達了降水量序列最高頻波動的信息,從IMF1到IMF6,振幅逐漸減小,頻率變低;其中殘差項是單調的,體現序列的趨勢,即降水量有上升的趨勢。

2.2 Hilbert變換的統(tǒng)計值

對每個IMF組件執(zhí)行Hilbert變換,以獲取每個組件的平均幅度和平均頻率,如表1所示。IMF1具有最高振幅,并且包含最高頻率的振動信息。這表明研究所用時間序列的年降水量的時間序列振蕩期約為9月。IMF1-IMF6的振幅逐漸減小,頻率減小。除了9個月左右的周期外,還有13、15、23、76和135個月左右的周期。

表1 年降水量各分量統(tǒng)計值

2.3 模型參數的確定

模型一不考慮初相位,各項參數為:IMF1平均周期8.8,平均振幅1.71;IMF2平均周期12.9,平均振幅1.92;IMF3平均周期15.1,平均振幅1.55;IMF4平均周期22.7,平均振幅0.58;IMF5平均周期76.3,平均振幅0.72;IMF6平均周期135.2,平均振幅0.64(表1),代入式(7),建立的降水量多尺度模型為:

模型二較模型一考慮到各IMF分項的初始相位,根據表1數據可得各項參數為:IMF1平均周期8.8,平均振幅1.71,初相-1.21;IMF2平均周期12.9,平均振幅1.92,初相-0.62;IMF3平均周期15.1,平均振幅1.55,初相-0.38;IMF4平均周期22.7,平均振幅0.58,初相0.20;IMF5平均周期76.3,平均振幅0.72,初相1.19;IMF6平均周期135.2,平均振幅0.64,初相1.53,建立的降水量多尺度模型為:

2.4 降水量預測

根據所建立的模型,對黃河三角洲2019年月平均降水量進行預測,并與實際降水量進行對比,結果見表2、圖2和圖3。

表2 2019年降水量預測值與實測值對比表

圖2 模型一預測值與實測值對比圖

圖3 模型二預測值與實測值對比圖

由表2可以看出,模型一對黃河三角洲降水量的預測誤差均在10%以下,相對誤差平均值為3.3%,最大相對誤差為9.8%,模型二對黃河三角洲降水量的預測誤差均在12%以下,相對誤差平均值為5.6%,最大相對誤差為11.8%。2種模型預測精度均表現出較高水平,無論是否考慮初相位的影響,2個模型在1月預測誤差均為年內最大,這可能與該時期為干旱期降水量量級較小有一定關系。

本次選用非特征年型2019年分別對模型一、模型二的預測數據進行誤差檢驗分析,并取得了不錯的效果,這在一定程度上說明2種模型在黃河下游降水量預測上具有一定的適用性。

3 討論

降水的形成易受地形、氣候、下墊面及人類活動等眾多因素的相互影響,降水量時間序列一般都具有非線性、非平穩(wěn)性的特點。降水量預測為降水量的定量研究、水資源的利用及河口治理提供科學依據,因此降水量非線性時間序列的預測十分重要。許多專家和學者采用各種方法對非線性時間序列進行建模預測,但是這些方法大部分是單尺度的[15]。近些年已有一些學者將EMD分解用于時間序列預測分析,取得了不錯的成果[16],但多為各分量累加預測,未能完全考慮高頻IMF分量的預測誤差。本研究對時間序列EMD分解后的不同尺度下的各分項進行Hilbert變換,得到其瞬時頻率、周期以及初相建立非線性多尺度區(qū)域降水量預測模型,考慮了降水量數據在不同尺度上的演變特性對預測結果的影響,同時也降低了高頻分量的預測難度,模型有較高精度。

模型一、模型二的具體參數是根據對黃河三角洲1954-2018年連續(xù)65 a月降水數據分解分析得到,因此2種模型對黃河三角洲降水預測具有較強的針對性,有助于提高預測的可信度。模型一相對誤差平均值為3.3%,模型二相對誤差平均值為5.6%,均小于國際水文組織認可的20%誤差范圍。與傳統(tǒng)的降水量預測方法相比,時間序列模型精度為7.5%[15],加權馬爾科夫鏈法精度為6.9%[7],ARIMA模型精度13.1%[8],人工神經網絡模型精度多在10%~20%之間[12],本模型預測精度均較高,這表明黃河三角洲降水量時間序列具有顯著的非線性特征,對其進行EMD分解,并進行Hilbert變換,新生成的序列在一定程度上能反映降水量內在的變化特征與演變機制,從而使得預測精度較高。

從研究結果看,在建模時不考慮初相位的模型一平均預測誤差為2.70%,整體預測精度要優(yōu)于考慮初相位的模型二,平均誤差5.63%,說明對于初相的引入雖在一定程度上完善了模型,但由于初相位確定的過程中存在的誤差可能影響整個模型的預測精度。在對12個月降水量預測精度的橫向對比中發(fā)現,模型對4-7月的預測精度要優(yōu)于其他8個月,與降水量年內分布具有負相關性。

本文僅僅針對已有實際降水量數據進行分析,從數據序列本身特征來揭示其演變規(guī)律,未來從下墊面變化等影響降水量的物理機制方面進行預測是需要進一步研究的重點。

4 結論

1)研究運用非線性多尺度模型對黃河三角洲1954—2018年連續(xù)65 a 780個月的降水量時間序列進行了EMD分解,得到6個IMF分量和1個殘差分量,在一定程度上反映了研究區(qū)降水量變化的影響機制和驅動因素。

2)本文所建模型一預測誤差在10%以下,相對誤差平均值為3.3%,最大相對誤差為9.8%;模型二預測誤差均在12%以下,相對誤差平均值為5.6%,最大相對誤差為11.8%。

3)2種模型預測精度均較高,在建模時不考慮初相位的模型一平均預測誤差為2.70%,整體預測精度要優(yōu)于考慮初相位的模型二,平均誤差5.63%,2個模型在1月預測誤差均為年內最大,這可能與該時期為干旱期,降水量絕對值相對較小有一定關系。

本研究方法為降水量的快速較高精度的預測提供了一次理論實踐探討,有利于促進黃河流域生態(tài)系統(tǒng)長久健康發(fā)展。

[1] LIEBMANN B, MARENGO J A. Interannual variability of the rainy season and rainfall in the Brazilian Amazon basin[J]. Journal of Climate, 2001, 14(22): 4 308-4 318.

[2] ZENG W, YU Z, WU S H, QIN J B. Changes in annual, seasonal and monthly precipitation events and their link with elevation in Sichuan Province, China[J]. International Journal of Climatology, 2016, 36(5): 2 303-2 322.

[3] 王英, 曹明奎, 陶波, 等. 全球氣候變化背景下中國降水量空間格局的變化特征[J]. 地理研究, 2006, 25(6): 1 031- 1 040, 1 148.

WANG Ying, CAO Mingkui, TAO Bo, et al. Variation characteristics of spatial pattern of precipitation in China under the background of global climate change[J]. Geographic Research, 2006, 25 (6): 1 031-1 040, 1 148.

[4] 劉田, 陽坤, 秦軍, 等. 青藏高原中、東部氣象站降水資料時間序列的構建與應用[J]. 高原氣象, 2018, 37(6): 1 449-1 457.

LIU Tian, YANG Kun, QIN Jun, et al. Construction and applications of time series of monthly precipitation at weather stations in the central and eastern Qinghai-Tibetan Plateau[J]. Plateau Meteorology, 2018, 37(6): 1 449-1 457.

[5] 劉向培, 王漢杰, 何明元. 應用統(tǒng)計降尺度方法預估江淮流域未來降水[J]. 水科學進展, 2012, 23(1): 29-37.

LIU Xiangpei, WANG Hanjie, HE Mingyuan. Estimation of precipitation under future climate scenarios in the Yangtze-Huaihe region using statistical downscaling[J]. Advances in Water Science, 2012, 23(1): 29-37.

[6] 劉健文, 周小剛. SSA方法在氣候時間序列分析和預測中的應用[J].氣象科技, 1996, 24(3): 18-22.

LIU Jianwen, ZHOU Xiaogang. Application of SSA method in climate time series analysis and prediction [J]. Meteorological Science and Technology, 1996, 24(3): 18-22.

[7] 梁顯麗, 寶秋利, 代海燕. 加權馬爾科夫鏈在鄂爾多斯市年降水量預測中的應用[J]. 數學的實踐與認識, 2021, 51(4): 161-171.

LIANG Xianli, BAO Qiuli, DAI Haiyan. Application of weighted Markov chain in prediction of annual precipitation in Ordos region[J]. Mathematics in Practice and Theory, 2021, 51(4): 161-171.

[8] 付明明. ARIMA模型在新疆喀什地區(qū)中長期降水量預測中的應用研究[J]. 地下水, 2019, 41(3): 142-144.

FU Mingming. Application of ARIMA model in prediction of medium and long-term precipitation in Kashi area, Xinjiang[J]. Ground Water, 2019, 41(3): 142-144.

[9] 楊沛羽, 張強, 史培軍, 等. 黃河流域極端降水時空分布特征及其影響因素[J]. 武漢大學學報(理學版), 2017, 63(4): 368-376.

YANG Peiyu, ZHANG Qiang, SHI Peijun, et al. Spatiotemporal distribution of precipitation extremes and related implications across the Yellow River Basin, China[J]. Journal of Wuhan University (Natural Science Edition), 2017, 63(4): 368-376.

[10] 杜懿, 龍鎧豪, 王大洋, 等. 基于機器學習方法的安徽省年降水量預測[J]. 水電能源科學, 2020, 38(7): 5-7, 41.

DU Yi, LONG Kaihao, WANG Dayang, et al. Annual precipitation prediction in Anhui Province based on machine learning[J]. Water Resources and Power, 2020, 38(7): 5-7, 41.

[11] 孔德萌, 李維德, 吳金冉. 基于協(xié)整理論的極限學習機模型在降水預測中的應用[J]. 水電能源科學, 2017, 35(9): 1-3, 12.

KONG Demeng, LI Weide, WU Jinran. Application of extreme learning machine model in precipitation prediction based on cointegration theory[J]. Water Resources and Power, 2017, 35(9): 1-3, 12.

[12] 何慧, 陸虹, 覃衛(wèi)堅, 等. 人工神經網絡在月降水量預測業(yè)務中的研究和應用綜述[J]. 氣象研究與應用, 2021, 42(1): 1-6.

HE Hui, LU Hong, QIN Weijian, et al. Research and application of artificial neural network in monthly precipitation forecast[J]. Journal of Meteorological Research and Application, 2021, 42(1): 1-6.

[13] HUANG N E, SHEN Z, LONG S R, et al. The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis[J]. Proceedings of the Royal Society of London Series A: Mathematical, Physical and Engineering Sciences, 1998, 454(1 971): 903-995.

[14] 任博, 薛澤宇, 任全志, 等. 基于EMD的凌河流域降水徑流預測模型研究[J]. 人民黃河, 2016, 38(6): 63-65.

REN Bo, XUE Zeyu, REN Quanzhi, et al. Rainfall runoff forecast application of BP prediction models based on EMD in linghe river basin[J]. Yellow River, 2016, 38(6): 63-65.

[15] 于?;? 時間序列模型在遼西降水量動態(tài)預測的應用[J]. 東北水利水電, 2019, (3): 42-44, 54.

YU Baohui. Application of time sequence model in precipitation dynamic forecast in western Liaoning province[J]. Water Resources & Hydropower of Northeastr, 2019, (3): 42-44, 54.

[16] 馬軍, 陸甲, 趙金彪. EMD方法在廣西夏季降水量預測中的應用[J]. 氣象研究與應用, 2014, 35(3): 31-35.

MA Jun, LU Jia, ZHAO Jinbiao. Application of EMD Method in Prediction of Summer Precipitation in Guangxi[J]. Journal of Meteorological Research and Application, 2014, 35(3): 31-35.

Precipitation Forecast in Yellow River Delta Based on Nonlinear Multi-scale Mode

MU Yuzhu1,WANG Yanpeng1,ZHANG Fenghua1,QIAO Dongmei2*

(1. Henan Xinxiang Hydrology and Water Resources Survey Bureau, Xinxiang 453000, China;2. Farmland Irrigation Research Institute, Chinese Academy of Agricultural Sciences, Xinxiang 453002, China)

【】Precipitation is closely related to human production, life and ecology. The change of precipitation is related to the sustainable utilization of regional water resources, the protection of ecological environment and the development of economy and society, The research on the variation characteristics and evolution trend of precipitation has become a hot topic in the field of climate and water resources. Scholars and researchers has paid much attention on the accurate prediction of precipitation.【】The purpose of this paper is to improve the prediction accuracy of precipitation, and reflect the actual characteristics of precipitation. 【】Based on the advantages of empirical mode decomposition in the analysis and processing of nonlinear time series and other fields, Empirical Mode Decomposition (EMD) was carried out for the monthly average precipitation data of the Yellow River Delta Meteorological Station from 1954 to 2018, and a series of eigenmode functions were obtained. Hilbert transform was performed on IMF, and on this basis, two multi-scale forecast models of precipitation in the Yellow River Delta were established.【】The results showed that there were periods of 9, 13, 23, 76 and 135 months in precipitation in the Yellow River Delta, and 9-month fluctuations were the main ones;The 65-year monthly average precipitation data was predicted. The relative error of the model 1 was between 0.9% and 9.8%, and the relative error of the model 2 was between 1.6% and 11.8%. When modeling, the average prediction error of model 1 without considering the initial phase was 2.70%, and the overall prediction accuracy was better than that of model 2 considering the initial phase.【】The fitting accuracy and significance of the two models meet the requirements.

precipitation; time series; multiscale; EMD; predict

S161.4

A

10.13522/j.cnki.ggps.2021149

1672 – 3317(2021)08 - 0123 - 06

穆玉珠, 王燕鵬, 張鳳華, 等. 基于非線性多尺度模型的黃河三角洲降水量預測[J]. 灌溉排水學報, 2021, 40(8): 123-128.

MU Yuzhu, WANG Yanpeng, ZHANG Fenghua, et al. Precipitation Forecast in Yellow River delta Based on Nonlinear Multi-scale Mode[J]. Journal of Irrigation and Drainage, 2021, 40(8): 123-128.

2021-04-14

穆玉珠(1974-),女,河南新鄉(xiāng)人。工程師,研究方向為水文水資源。E-mail: muyuzhu1020@126.com

喬冬梅(1978-),女。副研究員,碩士生導師,博士,研究方向為水資源與水環(huán)境。E-mail: qiaodongmei78@163.com

責任編輯:趙宇龍

猜你喜歡
模型
一半模型
一種去中心化的域名服務本地化模型
適用于BDS-3 PPP的隨機模型
提煉模型 突破難點
函數模型及應用
p150Glued在帕金森病模型中的表達及分布
函數模型及應用
重要模型『一線三等角』
重尾非線性自回歸模型自加權M-估計的漸近分布
3D打印中的模型分割與打包
主站蜘蛛池模板: 91无码国产视频| 99草精品视频| 日日碰狠狠添天天爽| 成人亚洲天堂| 国产欧美日韩综合在线第一| 亚洲第一色视频| 亚洲首页在线观看| 人妻免费无码不卡视频| 亚洲精品视频免费| 99re在线观看视频| 美美女高清毛片视频免费观看| 亚洲无码日韩一区| 国产在线精品99一区不卡| 丁香综合在线| 少妇精品网站| 国产福利免费观看| 亚洲 欧美 偷自乱 图片| 亚洲天堂网站在线| 亚洲国语自产一区第二页| 国产在线观看精品| 欧美人人干| 麻豆精品在线| 怡春院欧美一区二区三区免费| 99人体免费视频| 国模沟沟一区二区三区| 永久成人无码激情视频免费| 99久久精品免费看国产免费软件 | 视频二区国产精品职场同事| 亚洲九九视频| 第九色区aⅴ天堂久久香| 亚洲精品自在线拍| 美女国产在线| 99在线视频免费观看| 日韩欧美网址| 99资源在线| 国产一区二区三区日韩精品| 亚洲视频免| 亚洲欧美在线综合一区二区三区| 亚洲中久无码永久在线观看软件 | 伊人色在线视频| 国内精品久久人妻无码大片高| 国产丰满大乳无码免费播放| 国产成在线观看免费视频| 色综合久久无码网| 一区二区自拍| 久久人与动人物A级毛片| 伊人久久久久久久| 亚洲Aⅴ无码专区在线观看q| 国产一级特黄aa级特黄裸毛片 | 少妇精品在线| 免费a在线观看播放| 欧美日韩中文国产va另类| 一级黄色网站在线免费看| 成人欧美在线观看| 91免费国产在线观看尤物| 亚洲第一区在线| 国产美女自慰在线观看| 国产亚洲精品无码专| 婷婷开心中文字幕| 欧美色综合网站| 久久久久88色偷偷| 69av免费视频| 午夜丁香婷婷| 成年A级毛片| 青青草一区二区免费精品| 熟女成人国产精品视频| 在线日本国产成人免费的| 成人在线观看不卡| 精品综合久久久久久97超人| 婷婷综合在线观看丁香| 亚洲第一av网站| 亚洲第一精品福利| 亚洲精品久综合蜜| 亚洲综合久久成人AV| 亚洲国产精品美女| 午夜在线不卡| 美女无遮挡被啪啪到高潮免费| 69国产精品视频免费| 免费视频在线2021入口| 久草网视频在线| 综合天天色| 欧美不卡视频在线观看|