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

淮河流域極端氣候事件非平穩(wěn)特征研究

2021-04-29 07:53:58王懷軍肖明賢
中國農村水利水電 2021年4期
關鍵詞:趨勢模型

王懷軍,曹 蕾,肖明賢,馮 如

(1.淮陰師范學院城市與環(huán)境學院,江蘇淮安223300;2.南京水利科學研究院水文水資源與水利工程科學國家重點實驗室,南京210029;3.南京水利科學研究院水利部應對氣候變化研究中心,南京210029)

氣候變化導致極端氣候事件的發(fā)生頻率、強度、空間范圍及持續(xù)時間發(fā)生改變,并可能引發(fā)前所未有的極端事件[1-3]。中國是世界上受自然災害影響最為嚴重的國家之一,也是極端氣候事件發(fā)生頻率與強度較高的國家之一[4]。隨著全球氣候變化以及中國社會經(jīng)濟快速發(fā)展對資源環(huán)境和生態(tài)產(chǎn)生的壓力加劇,使得極端氣候事件防范應對形勢更加嚴峻復雜[5]。據(jù)民政部和國家減災委員會等部門統(tǒng)計,1990-2010年氣象災害造成的直接經(jīng)濟損失高達2 000 億元/a,約相當于當年GDP 的1%~3%[6]。并且近年來氣象災害損失有日益加重趨勢,氣象災害的頻繁發(fā)生已成為制約我國國民經(jīng)濟持續(xù)穩(wěn)定發(fā)展的主要因素之一。

目前,在極端氣候事件時空變化研究方面,主要從以下方面進行研究:①極端氣候事件的表征。基于單一臺站監(jiān)測技術的定義,采用“氣候極值”量化異常天氣氣候現(xiàn)象,當氣候要素達到定義的閾值時,便可認為極端事件發(fā)生[1,7,8]。顯然,閾值的確定受人為因素的影響較大;②極端氣候事件時空變化特征。主要利用趨勢檢驗、突變分析、小波變換以及地統(tǒng)計學技術分析某流域或者地區(qū)的極端氣候事件時空分布規(guī)律[9,10];例如,尹義星等[11]對1951-2013年江蘇省極端最高和最低氣溫變化趨勢進行研究,結果表明冰凍日數(shù)和冷晝日數(shù)呈下降趨勢,夏季日數(shù)和暖晝日數(shù)呈上升趨勢;③基于平穩(wěn)極值函數(shù)的時空分布特征。極端氣候事件是一個隨機變量,利用廣義極值分布(GEV)和廣義帕累托分布(GPD)對其進行擬合以及重現(xiàn)期水平提取[12,13]。如趙玲玲等[14]利用GEV 對華南強降雨特性進行分析,余錦華等[15]采用GPD 對中國東部夏季極端降水頻率特征進行分析,結果均表明GEV和GPD模型可以很好的擬合極端氣候事件。然而,氣候變化和人類活動導致極端氣候事件呈非平穩(wěn)性,傳統(tǒng)極值函數(shù)分析中,極端氣候事件的模擬基于平穩(wěn)的邊際分布函數(shù),這往往忽視氣候變化背景下水文變量非平穩(wěn)特性的影響[16]。溫慶志等[17]指出淮河流域年最高氣溫表現(xiàn)出非平穩(wěn)性,增幅達1.5 ℃。Kalai 等[18]研究結果表明由于土地使用方式、人類干預和氣候變化的迅速變化,非平穩(wěn)區(qū)域洪水頻率分析(RFFA)能夠更有效地捕獲洪水的時變行為。Gu等[19]指出城市化導致局部地區(qū)(例如中國北方)的極端降水呈現(xiàn)非平穩(wěn),城市化地區(qū)氣候發(fā)生嚴重的非平穩(wěn)性的可能性更高。Yang 等[20]則指出在平穩(wěn)和非平穩(wěn)條件下檢測到的城市洪水存在差異,需要將更多的精力進一步研究非平穩(wěn)條件下的城市洪水變化。因此,有必要將數(shù)據(jù)的非平穩(wěn)性引入到GEV 和GPD 函數(shù)中,以便更好地揭示極端氣候事件的時空分布特征。

淮河流域地處我國南北氣候過渡帶,是中國重要的地理生態(tài)分界線和生態(tài)環(huán)境脆弱區(qū),對全球氣候變化十分敏感[21]。由于淮河流域本身特殊的地理環(huán)境條件,氣候系統(tǒng)、水循環(huán)系統(tǒng)的不穩(wěn)定性增大,導致極端氣候事件發(fā)生的頻率和強度都有所增大,這對人民生產(chǎn)生活產(chǎn)生了不可忽視的影響[22]。本研究選擇位于南北氣候過渡帶這一特殊地理環(huán)境下的淮河流域為研究區(qū),借助非平穩(wěn)GPD 和GEV 模型擬合極端氣候事件,揭示非平穩(wěn)條件下極端氣候事件時空演變規(guī)律,以期為淮河流域減災防災提供決策參考。

1 研究區(qū)概況、數(shù)據(jù)與方法

1.1 研究區(qū)概況

淮河流域地處我國東部,介于中國兩大流域-長江與黃河之間,其地理位置位于30°55′~36°36′N,111°55′~121°25′E 之間(見圖1),流域面積27 萬km2。淮河流域包括鄂、豫、皖、魯、蘇五省40 個地(市)181 個縣(市),流域內總人口1.65 億人,平均人口密度居各大江大河流域人口密度之首,為611 人/km2。流域地處我國南北氣候過渡地帶,年平均氣溫在11~16 ℃之間。最高月平均氣溫出現(xiàn)在7月份,約為25 ℃,最低月平均氣溫出現(xiàn)在1月份,約為0 ℃。極端氣溫最高達44.5 ℃,最低達-24.1 ℃,是中國極端氣溫事件的高發(fā)區(qū)之一。多年平均降水量約為920 mm,其分布情況大致是由南向北遞減,平原少,山區(qū)多,內陸少于沿海。

圖1 研究區(qū)和氣象站點分布圖Fig.1 Study area and distributed weather stations

1.2 數(shù) 據(jù)

1960-2018年淮河流域內70 個站點的每日降水和氣溫數(shù)據(jù)來自于中國氣象科學數(shù)據(jù)共享服務網(wǎng)(http://www.cdc.cma.gov.cn),數(shù)據(jù)經(jīng)過重重檢查,質量良好。本文首先借助R 語言提取年最大值序列(AM)和超門限序列(POT)。AM 序列包括1日最低氣溫(日最低氣溫的年最低值,Tnn)、1日最高氣溫(日最高氣溫的年最大值,Txx)和1日最大降水量(日降水量的年最大值,Rx1day)。POT 序列閾值取第95%百分位的日最高氣溫、日最低氣溫和日降水量,分別命名為超門限最低氣溫(Tmingpd)、超門限最高氣溫(Tmaxgpd)和超門限降水(Pregpd)。

1.3 方 法

1.3.1 平穩(wěn)性檢驗

時間序列的平穩(wěn)性可以表示如下:①均值E(X)=μ是與時間t無關的常數(shù);②方差Var(Xt)=σ2是與時間t無關的常數(shù);③協(xié)方差Cov(Xt,Xt+k)=γ k是只與時期間隔k有關,與時間t無關的常數(shù)。常用的判斷方法包括圖示法、單位根檢驗法、自相關函數(shù)(ACF)判斷法。圖示法通過時間序列的路徑圖來粗略地判斷它是否是平穩(wěn)的,平穩(wěn)序列表現(xiàn)出一種圍繞其均值不斷波動的過程,非平穩(wěn)序列則往往表現(xiàn)出一定的上升或者下降趨勢。為方便進行處理,我們采用趨勢和突變檢驗來進行檢驗,如果時間序列存在顯著趨勢和突變點,則說明該序列非平穩(wěn)。趨勢檢驗采用Man-Kendall 檢驗,均值突變檢驗采用Pettitt 檢驗,方差突變采用單變異點檢測方法(AMOC)[23,24]。Man-Kendall 檢驗通過計算統(tǒng)計量Zs來進行趨勢統(tǒng)計的顯著性檢驗[25,26]。當Zs>0 時,表示時間序列呈增加趨勢,Zs<0 時,表示時間序列呈減少趨勢,如果Zs絕對值大于1.96,表明時間序列趨勢顯著(0.05 顯著水平)。Man-Kendall 檢驗中還計算了傾斜度Sen,表示趨勢的方向和大小。計算公式分別為:

式中:1<i<j<n,n為數(shù)據(jù)資料時間長度,如果Sen的值為正,表示呈現(xiàn)增加趨勢,如果Sen的值為負則表示呈現(xiàn)減小趨勢。

Pettitt 非參數(shù)變點檢驗假設一個水文時間序列y1,y2,…,yn中的突變點出現(xiàn)在t時刻,將t作為分界點,將該水文時間序列分為兩部分:y1,y2,…,yt和yt+1,yt+2,…,yn分別遵循分布Fy1(y)和Fy2(y),然后檢驗在顯著性水平α下突變點的顯著度,具體計算過程參照文獻[27,28]。

另外文章采用ADF 檢驗進行單位根檢驗,對于時間序列AR(p)模型,有

式中:c為漂移項;β、φi為參數(shù);εt為隨機誤差項,是服從獨立分布的白噪聲過程。

ADF 檢驗的原假設H0∶β=0,存在單位根,為非平穩(wěn)序列,從而可以根據(jù)ADF單位根檢驗結果判斷序列是否平穩(wěn)。

1.3.2 平穩(wěn)和非平穩(wěn)GEV模型

GEV模型適合用來處理AM序列(如年最大降水量、年極端最高和最低氣溫),其累積分布函數(shù)為:

式中:μ,σ,ξ分別為位置、尺度和形狀參數(shù)。

GEV模型T年重現(xiàn)水平為:

式中:P為不同重現(xiàn)期對應的概率。

由于在方差突變檢驗中,所有的極值時間序列均不顯著,且形狀參數(shù)一般不易改變,因此非平穩(wěn)GEV 分布模型僅設置為位置參數(shù)隨時間變化[μ(t)=ρ0+ρ1t],形式如下:

1.3.3 平穩(wěn)和非平穩(wěn)GPD模型

POT序列采用GPD分布進行模擬,其累計分布函數(shù)是:

式中:u,σ,ξ分別為分布的閾值、尺度和形狀參數(shù)。

假定某一極端事件的重現(xiàn)期為1/(1-p),則其對應的重現(xiàn)水平為:

由于閾值為95%分位數(shù),形狀參數(shù)一般不易改變,因此非平穩(wěn)GPD 分布模型定義如下[σ(t)=sin(2 π×(t-1960)/365.25]+cos[2π×(t-1960)/365.25)]:

對平穩(wěn)和非平穩(wěn)模型進行似然比檢驗(LR 檢驗),統(tǒng)計量為:

式中:L1為復雜模型最大似然值;L2為簡單模型最大似然值;LR近似符合卡方分布。

對于GEV 模型,分別構造了GEV0(無協(xié)變量)和GEV1(位置參數(shù)隨時間線性變化);對于GPD 構造GPD0和GPD1(尺度參數(shù)隨時間線性變化)。通過構造的復雜模型與簡單模型比較來檢驗加入非平穩(wěn)性后是否能夠顯著地提高模擬效率。

最后利用選擇的平穩(wěn)或非平穩(wěn)GPD和GEV模型,通過設置不同的重現(xiàn)期,得到每個站點2年、10年、20年、30年、50年、100年一遇的極端氣候事件重現(xiàn)水平,利用克里金插值得到重現(xiàn)水平空間分布圖。

2 結果分析

2.1 極端氣候事件時空變化趨勢及平穩(wěn)性分析

關于最高氣溫,淮河流域Tmaxgpd和Txx變化趨勢為0.07 d/a和0.01 ℃/a,且變化走勢相對一致,都為1960-1985呈減少趨勢,而在1985年以后呈增加趨勢[圖2(a),圖2(d)]。從空間變化趨勢來看,上述兩個指標大部分站點表現(xiàn)為增加趨勢,但表現(xiàn)為顯著變化的站點相對較少,分別為14%和21%[表1,圖3(a),圖3(d)]。對于最低氣溫,Tmingpd和Tnn均為顯著增加趨勢,變化幅度分別為0.16 d/a和0.04℃/a[圖2(b),圖2(e)],其變化幅度均高于最高氣溫,這也說明與最低氣溫相關的極值要高于與最高氣溫相關的極值,這與其他地區(qū)的研究結果一致[29]。空間變化趨勢具有與區(qū)域時間一致的變化趨勢,顯著增加的站點分別達到了59%和76%。Pregpd和Rx1day的區(qū)域變化趨勢分別為0.01 d/a 好0.11 mm/a,但不管是從區(qū)域上還是空間上均為不顯著變化。

圖2 淮河流域極端氣候事件區(qū)域時間變化趨勢Fig.2 Regional times series of climate extremes in the Huaihe River Basin

表1 極端氣候事件空間變化趨勢百分比 %Tab.1 Percentage showed different trends of extreme climate events

為進一步分析極端氣候指數(shù)的平穩(wěn)性,采用了ADF 檢驗和Pettitt 檢驗進行分析(圖4)。79%站點的Tmaxgpd和87%站點的Tmingpd的ADF 檢驗P值大于0.05,說明具有單位根,表現(xiàn)為非平穩(wěn);Txx和Tnn具有單位根的百分比為46%和36%,說明在此檢驗下表現(xiàn)為非平穩(wěn)的站點百分比不到一半;Pregpd和Rx1day僅有少部分站點具有單位根,說明在淮河流域降水極值表現(xiàn)為平穩(wěn)性變化。在Pettitt 檢驗下[圖4(b)],Tmingpd和Tnn變現(xiàn)為顯著變化的站點分別為71%和74%,說明大部分站點在此統(tǒng)計量下表現(xiàn)為非平穩(wěn)。其他指數(shù)表現(xiàn)為顯著變化的站點比例很少,均小于30%,且突變年份分布很廣[圖4(c)],說明大部分站點表現(xiàn)為非平穩(wěn)變化。從趨勢分析和相關平穩(wěn)性結果可以看出,Tmingpd和Tnn在淮河流域表現(xiàn)為非平穩(wěn)性,Pregpd和Rx1day表現(xiàn)為平穩(wěn)性,Tmaxgpd和Txx表現(xiàn)相對較為復雜,平穩(wěn)性和非平穩(wěn)性同時存在。

2.2 基于GEV模型的極端氣候事件變化特征

對AM 序列進行GEV 分析。首先利用平穩(wěn)性GEV 模型對所有站點的GEV 序列進行KS 擬合優(yōu)度檢驗[圖5(a)],發(fā)現(xiàn)幾乎所有的站點的P值均大于0.05,這說明在KS 統(tǒng)計量下所有AM序列均滿足平穩(wěn)性GEV模型。該研究結果表明在進行非平穩(wěn)性分析時,不能通過擬合優(yōu)度檢驗來探查數(shù)據(jù)的平穩(wěn)性。通過LR 檢驗來分析非平穩(wěn)模型是否對平穩(wěn)模型有改進[圖5(b)],對所有站點的AM 序列進行平穩(wěn)和非平穩(wěn)GEV 函數(shù)擬合(圖6)。結果表明Tnn80%的站點P值小于0.05,這說明大部分站點的非平穩(wěn)模型對平穩(wěn)模型有改進。而Txx和Rx1day,非平穩(wěn)模型改進的比率僅有30%和1%,這說明大部分站點滿足平穩(wěn)性。所以在后續(xù)GEV 分析中,Tnn采用非平穩(wěn)模型,而Txx和Rx1day采用平穩(wěn)性模型。

圖3 淮河流域極端氣候事件空間變化趨勢Fig.3 Spatial trend of extreme climate events in the Huaihe River Basin

圖4 極端氣候事件平穩(wěn)性檢驗及突變年份Fig.4 Stationary test and the year of abrupt change of extreme climate events

對AM 序列的2年(RTL2)、10年(RTL10)、20年(RTL20)、30年(RTL30)、50年(RTL50)和100年(RTL100)重現(xiàn)期重現(xiàn)水平進行求取,發(fā)現(xiàn)不用重現(xiàn)期下重現(xiàn)水平具有較高的相關系數(shù)(圖7),這說明不同重現(xiàn)期下極端氣候事件重現(xiàn)水平具有一致的空間分布特征,因此,后續(xù)研究僅分析30年重現(xiàn)期重現(xiàn)水平的空間分布(圖8)。淮河流域Rx1day30年重現(xiàn)水平西北向東南遞增[圖8(a)],這可能與東南季風的逐漸增強有關。Txx從西往東減少[圖8(b)],由于海洋對流域東部地區(qū)氣溫的調節(jié)作用,離海洋越近受海洋的影響越大,導致Txx自東往西逐漸升高。由于將Tnn視為非平穩(wěn)模型,重現(xiàn)水平會隨著時間而變化。為方便表示,將重現(xiàn)水平分為兩個時間段,1960-1989[圖8(c)]和1990-2018[圖8(d)],分別對兩個時間段重現(xiàn)水平求均值。結果顯示,Tnn-pre和Tnn-post具有一致的空間分布趨勢,重現(xiàn)水平從南向北遞增,這與緯度地帶性規(guī)律一致。但值得注意的是兩個時期的Tnn相差了近乎1~2 ℃,這說明在淮河流域Tnn的增溫趨勢非常明顯。

圖5 淮河流域極端氣候事件AM序列擬合優(yōu)度檢驗Fig.5 Goodness-of-fit test of extreme climate events of AM series in the Huaihe River Basin

圖6 基于平穩(wěn)和非平穩(wěn)模型的泗洪站Tnn(取絕對值)擬合優(yōu)度圖Fig.6 Goodness-of-fit test of absolute Tnn of Sihong station based on stationary and non-stationary models

2.3 基于GPD模型的極端氣候事件變化特征

采用平穩(wěn)和非平穩(wěn)性GPD 模型對所有站點POT 序列進行擬合,并利用LR 檢驗驗證模型是否有改進。結果表明,7%、78%和93%的非平穩(wěn)GPD 模型對降水、最高氣溫和最低氣溫有改進(圖9、圖10),這說明Pregpd表現(xiàn)為平穩(wěn)變化,Tmaxgpd和Tmingpd表現(xiàn)為非平穩(wěn)變化。因此,采用平穩(wěn)GPD 模型對Pregpd模型進行擬合,非平穩(wěn)模型對Tmaxgpd和Tmingpd進行擬合。同Rx1day空間分布特征一致,Pregpd空間部分從西北向東南遞降[圖11(a)],同樣反映了東亞季風從東南向西北逐步減弱的態(tài)勢。將Tmaxgpd和Tmingpd重現(xiàn)水平分為兩個時間段,1960-1989(pre)和1990-2008(post),分別對兩個時段重現(xiàn)水平求均值,分析其空間分布狀態(tài)。Tmaxgpd-post和Tmaxgpdpre反映了同樣的分布趨勢,從東部向西部遞增[圖11(b)、圖11(c)],這反映了海洋的熱力調節(jié)作用。前后兩個時期的重現(xiàn)水平相對變化較少,且有正有負[圖10(f)],這說明Tmaxgpd30年重現(xiàn)水平并不一定為增加趨勢,這可能與部分站點并不為非平穩(wěn)性有關。Tmingpd-pre和Tmingpd-post空間分布趨勢具有一致性,從南向北遞降[圖11(d)、圖11(e)]。相較于Tmaxgpd,Tmingpd-post比Tmingpd-pre數(shù)值都大,且大部分地區(qū)超過0.5 ℃以上[圖11(g)],這說明最低氣溫重現(xiàn)水平呈顯著增加趨勢。

圖7 30年重現(xiàn)水平與2、10、20、50、100年重現(xiàn)水平相關系數(shù)圖Fig.7 Correlation coefficients between the 30-year return level and the 2,10,20,50,and 100-year return levels

圖8 AM序列30年重現(xiàn)水平空間分布圖(其中Tnn取絕對值)Fig.8 Spatial distribution of the 30-year return level of AM extremes

3 結 論

基于淮河流域最新氣象站點日值氣溫和降水數(shù)據(jù),計算目前國際上常用的極端氣候事件指數(shù),分析極端氣候事件的時空分布趨勢。基于平穩(wěn)和非平穩(wěn)GEV和GPD極值模型,分析極端氣候事件在不同極值模型下的重現(xiàn)水平空間變化特征,主要結果如下。

(1)對于最低氣溫,超門限最低氣溫(Tmingpd)和1日最低氣溫(Tnn)呈顯著增加趨勢,且空間上亦表現(xiàn)為顯著增加趨勢。超門限最高氣溫(Tmaxgpd)、1日最高氣溫(Txx)在區(qū)域上表現(xiàn)為不顯著增加,且發(fā)生顯著變化的站點百分比相對較低。超門限降水(Pregpd)和1日最大降水量(Rx1day)不管是從區(qū)域上還是空間上均表現(xiàn)為不顯著變化。

圖9 淮河流域極端氣候事件POT序列LR檢驗Fig.9 LR test of extreme climate events of POT series in the Huaihe River Basin

(2)不同的平穩(wěn)性檢驗方法會得出一定差異的平穩(wěn)性檢驗結果,單位根檢驗下Tmaxgpd和Tmingpd表現(xiàn)為非平穩(wěn),Pregpd和Rx1day表現(xiàn)為平穩(wěn)性變化。在Pettitt檢驗下,Tmingpd和Tnn表現(xiàn)為非平穩(wěn)。

(3)AM 序列中,Tnn大部分站點非平穩(wěn)模型對平穩(wěn)模型有改進,Txx和Rx1day非平穩(wěn)模型改進的比率相對較低,說明在AM 序列中,Tnn宜采用非平穩(wěn)模型,而Txx和Rx1day采用平穩(wěn)模型。淮河流域Rx1day30年重現(xiàn)水平西北向東南遞增,Txx從西往東減少,Tnn絕對值從南向北遞增,且不同時期重現(xiàn)水平變化顯著。

圖10 基于平穩(wěn)和非平穩(wěn)GPD模型的泗洪Tmingpd擬合優(yōu)度圖Fig.10 Goodness-of-fit test of Tmingpd of Sihong station based on stationary and non-stationary models

圖11 POT 序列30年重現(xiàn)水平空間分布Fig.11 Spatial distribution of the 30-year return level of POT extremes

(4)POT 序列中,Tmaxgpd和Tmingpd非平穩(wěn)模型改進的比率超過78%,這兩個指數(shù)宜采用非平穩(wěn)模型,Pregpd改進比率較低,宜采用平穩(wěn)GPD 模型。Pregpd重現(xiàn)水平空間分布從西北向東南遞減,Tmaxgpd從東部向西部遞增,Tmingpd從南向北遞降。 □

猜你喜歡
趨勢模型
一半模型
趨勢
重要模型『一線三等角』
重尾非線性自回歸模型自加權M-估計的漸近分布
初秋唇妝趨勢
Coco薇(2017年9期)2017-09-07 21:23:49
3D打印中的模型分割與打包
SPINEXPO?2017春夏流行趨勢
“去編”大趨勢
FLUKA幾何模型到CAD幾何模型轉換方法初步研究
趨勢
汽車科技(2015年1期)2015-02-28 12:14:44
主站蜘蛛池模板: 久久频这里精品99香蕉久网址| 亚洲成a人片| 国产九九精品视频| 亚洲欧美在线精品一区二区| 熟女成人国产精品视频| 国产不卡网| 免费jizz在线播放| 亚洲精品天堂在线观看| 狠狠综合久久| 99热这里只有成人精品国产| 久草视频精品| 日韩毛片免费视频| 亚洲码在线中文在线观看| 国产尤物在线播放| Jizz国产色系免费| 日韩在线欧美在线| 老司国产精品视频91| 91久久国产综合精品女同我| 一级黄色网站在线免费看| 国产成人亚洲毛片| 免费看美女自慰的网站| 中文字幕波多野不卡一区| 波多野吉衣一区二区三区av| 中文字幕人成人乱码亚洲电影| 重口调教一区二区视频| 免费黄色国产视频| 国产亚洲精品自在线| 精品少妇人妻一区二区| 国产99免费视频| 无码人妻免费| 亚洲欧美在线综合一区二区三区| 一边摸一边做爽的视频17国产| 女同久久精品国产99国| 欧美一级一级做性视频| 国产www网站| 日韩在线网址| 亚洲高清无码久久久| 99视频在线看| 日韩一区二区在线电影| 91香蕉视频下载网站| 亚洲精品中文字幕无乱码| 午夜精品福利影院| 99久久亚洲综合精品TS| 成人毛片免费观看| 亚洲日本一本dvd高清| 国产精品林美惠子在线观看| 国产美女在线观看| 黄色网站在线观看无码| 国产成人啪视频一区二区三区| 在线日韩日本国产亚洲| 国产熟女一级毛片| 免费国产不卡午夜福在线观看| 国产第八页| 亚洲第一天堂无码专区| 欧美精品H在线播放| 亚洲国产亚综合在线区| 国产精品开放后亚洲| 国产欧美日韩在线在线不卡视频| 国产成本人片免费a∨短片| 国产成熟女人性满足视频| 成人福利免费在线观看| 午夜丁香婷婷| 狠狠色丁香婷婷综合| 国产高清在线观看91精品| 国产a网站| 色视频久久| 日韩免费中文字幕| 日韩视频精品在线| 狠狠综合久久| 欧美精品亚洲二区| 久久精品国产免费观看频道| 亚瑟天堂久久一区二区影院| 国产精品视频系列专区| 99re视频在线| 日韩福利视频导航| 久久国产乱子| 亚洲第一成年网| 亚洲综合二区| 亚洲精品片911| 伊人大杳蕉中文无码| 亚洲成a人片77777在线播放 | 久99久热只有精品国产15|