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

集合經(jīng)驗(yàn)?zāi)B(tài)分解和交叉小波在贛江流量與太陽(yáng)黑子活動(dòng)關(guān)系中的應(yīng)用

2018-03-16 01:52:48萬(wàn)智巍賈玉連洪祎君蔣梅鑫
關(guān)鍵詞:分析

萬(wàn)智巍, 賈玉連, 洪祎君, 蔣梅鑫

(1.江西師范大學(xué) 地理與環(huán)境學(xué)院 鄱陽(yáng)湖濕地與流域研究教育部重點(diǎn)實(shí)驗(yàn)室,江西 南昌 330022; 2.中國(guó)科學(xué)院地理科學(xué)與資源研究所 陸地表層格局與模擬院重點(diǎn)實(shí)驗(yàn)室, 北京 100101)

1 研究背景

河流的水文要素變化過(guò)程是氣候變化和下墊面要素交互作用的結(jié)果,具有復(fù)雜的時(shí)空變化演變規(guī)律[1]。近年來(lái),隨著全球極端水文、氣候事件頻發(fā),了解河川徑流時(shí)空變化規(guī)律成為社會(huì)的迫切需求[2-3]。地球作為一個(gè)復(fù)雜的開(kāi)放系統(tǒng),太陽(yáng)輻射是其重要的外部能量輸入和影響因素。研究表明,全球和區(qū)域降水在長(zhǎng)時(shí)間尺度上受到太陽(yáng)活動(dòng)的影響,但其具體的影響過(guò)程和相互關(guān)系在不同區(qū)域有著不同的表現(xiàn)形式和規(guī)律特征[4]。Neff等[5]的研究論證了全新世以來(lái),太陽(yáng)活動(dòng)與印度季風(fēng)之間復(fù)雜的耦合關(guān)系。Zhang等[6]的研究也表明,太陽(yáng)活動(dòng)的變化對(duì)東亞區(qū)域降水和河流流量有著復(fù)雜的影響機(jī)制。竇睿音等[7]的研究表明,太陽(yáng)黑子數(shù)與關(guān)中平原的整體干濕狀況的關(guān)系存在轉(zhuǎn)換現(xiàn)象,20世紀(jì)70年代之前二者為反向相關(guān),之后二者之間的關(guān)系轉(zhuǎn)換為正相關(guān)。

河流流量的變化是區(qū)域氣象水文活動(dòng)的綜合體現(xiàn),反映了流域范圍內(nèi)整體的干濕狀況,因此流量變化成為氣候變化區(qū)域響應(yīng)的指示器[8]。相關(guān)研究表明,盡管太陽(yáng)黑子活動(dòng)與河流流量之間存在相關(guān)性,但是由于氣象水文過(guò)程屬于非線性、非穩(wěn)定的復(fù)雜系統(tǒng),如利用傳統(tǒng)統(tǒng)計(jì)方法分析其相關(guān)性可能會(huì)受到原始序列中噪音的干擾[9],難以探究二者之間在不同時(shí)域和頻域上的復(fù)雜相關(guān)關(guān)系。因此,本文結(jié)合集合經(jīng)驗(yàn)?zāi)B(tài)分解(Ensemble Empirical Mode Decomposition, EEMD)[10]和交叉小波變換(Cross Wavelet Transform, CWT)[11],提出一種改進(jìn)的相關(guān)分析方法——EEMD-CWT綜合分析,并將其應(yīng)用于贛江長(zhǎng)時(shí)間流量序列和同期太陽(yáng)黑子數(shù)序列之間的相關(guān)性分析。文中結(jié)果一方面可以為贛江流域水資源和干濕狀況長(zhǎng)期變化規(guī)律和影響因素的探討提供基礎(chǔ)資料,另一方面可以為提高氣象水文系統(tǒng)與外部強(qiáng)迫之間的非線性關(guān)系研究提供參考和思路。

2 數(shù)據(jù)與方法

2.1 流量數(shù)據(jù)與太陽(yáng)黑子數(shù)據(jù)

贛江的流量數(shù)據(jù)來(lái)自1950-2015年南昌外洲水文站的實(shí)測(cè)資料。外洲水文站位于江西省南昌市范圍內(nèi)的贛江下游右岸(115.84°E、28.63°N),控制面積為80 948 km2,占整個(gè)贛江流域范圍的99%以上[12],因此外洲水文站的流量數(shù)據(jù)可以反映整個(gè)贛江流域的產(chǎn)流情況。外洲水文站設(shè)立于1949年10月,從1950年開(kāi)始有連續(xù)的流量觀測(cè)數(shù)據(jù),是目前贛江流域最長(zhǎng)的流量數(shù)據(jù)序列(圖1(a))。

太陽(yáng)黑子數(shù)據(jù)來(lái)自比利時(shí)皇家天文臺(tái)數(shù)據(jù)中心,該研究機(jī)構(gòu)是世界權(quán)威的太陽(yáng)活動(dòng)觀測(cè)中心,因此本文選取其發(fā)布的1950-2015年太陽(yáng)黑子相對(duì)數(shù)(圖1(b))作為反映太陽(yáng)活動(dòng)狀況的序列(數(shù)據(jù)下載網(wǎng)址為http://www.sidc.be/silso/datafiles)。

圖1 1950-2015年贛江年均流量和太陽(yáng)黑子數(shù)

由圖1可以觀察到,流量序列與太陽(yáng)黑子序列之間的相關(guān)關(guān)系并不明顯。因此,利用SPSS 20軟件,對(duì)流量數(shù)據(jù)與太陽(yáng)黑子數(shù)據(jù)進(jìn)行初步的相關(guān)分析。結(jié)果表明,二者之間的Pearson相關(guān)系數(shù)為-0.118(P=0.347)。

2.2 研究方法

2.2.1EEMD方法 經(jīng)驗(yàn)?zāi)B(tài)分解(EmpiricalModeDecomposition,EMD)是一種利用各信號(hào)周期不同的振幅頻率特征,將原始序列分解為一系列具有不同時(shí)間尺度本征模函數(shù)(IntrinsicModeFunction,IMF)的方法,被廣泛應(yīng)用于水文、氣象、地球物理等領(lǐng)域的多尺度時(shí)頻分析[13-15]。集合經(jīng)驗(yàn)?zāi)B(tài)分解(EEMD)則是EMD方法的改進(jìn)算法,通過(guò)在原始信號(hào)中加入噪聲減弱了模態(tài)之間的混疊,從而提高了分解的有效性[16-17]。

2.2.2CWT方法 交叉小波變換(CrossWaveletTransform,CWT)是一種結(jié)合交叉譜分析和小波分析的新方法,可以同時(shí)在時(shí)域和頻域上分析兩個(gè)信號(hào)之間的相關(guān)性[11]。交叉小波能量譜可以反映出兩個(gè)信號(hào)經(jīng)過(guò)小波分析之后具有相同能量譜的區(qū)域,因此可以揭示出兩個(gè)信號(hào)在不同時(shí)頻域上相互作用的程度。其具體計(jì)算方法為:

(1)

2.2.3EEMD-CWT綜合分析 本文提出的EEMD-CWT綜合分析的思路是:首先利用EEMD方法將原始信號(hào)分解為不同尺度的互相正交的本征模函數(shù)(IMF)分量,以實(shí)現(xiàn)信號(hào)的平穩(wěn)化,降低原始信號(hào)中的噪音影響。然后利用多窗譜分析(Multi-TaperMethod,MTM)[18]確定原始信號(hào)中所存在的各種周期。最后,選擇具有相同物理意義的周期尺度下的本征模函數(shù)(IMF)分量進(jìn)行交叉小波變換,得到二者的交叉小波能量譜和小波凝聚譜,并以此來(lái)分析兩個(gè)序列在不同時(shí)頻域上的相關(guān)性。

3 結(jié)果與分析

3.1 流量的多尺度特征和周期

對(duì)1950-2015年贛江流量序列進(jìn)行EEMD分解(圖2),共得到6個(gè)IMF分量,其中IMF1~5為不同時(shí)間尺度下流量序列的周期項(xiàng),IMF6為整個(gè)序列的長(zhǎng)期變化趨勢(shì)項(xiàng)。由圖2可以看出,IMF1分量的波動(dòng)振幅最大,頻率也最高,但是周期性較差。這一信號(hào)代表了原始序列中的主要變化因素,是水文序列非線性和非平穩(wěn)特征的體現(xiàn)。EMD方法的提出者——Huang等[18]的研究指出,如果趨勢(shì)項(xiàng)為單調(diào)上升或單調(diào)下降函數(shù)則說(shuō)明原始信號(hào)具有明顯的非平穩(wěn)性。由圖2可以清楚的觀察到IMF6為單調(diào)上升函數(shù),這也論證了贛江近66年來(lái)的流量序列具有明顯的非平穩(wěn)性。IMF2~5這4個(gè)分量的周期性和平穩(wěn)性已經(jīng)得到了很大程度上的改善,表現(xiàn)出了比較明顯的周期變化規(guī)律。并且隨著分解尺度的遞進(jìn),其振動(dòng)幅度和頻率逐漸下降,信號(hào)的平穩(wěn)性逐漸增強(qiáng)。如IMF4和IMF5已經(jīng)比較接近于正弦波。

(圖中縱坐標(biāo)為相對(duì)值)

通過(guò)對(duì)1950-2015年贛江流量序列的多窗譜(MTM)分析表明(圖3),近66年來(lái)贛江流量序列具有明顯的22、8、4.5和3.7a周期(通過(guò)α=0.05顯著性檢驗(yàn)),其中11a周期也比較明顯,但是未能通過(guò)相關(guān)檢驗(yàn)。譜值最大的22a周期即為太陽(yáng)活動(dòng)的Hale周期[19],也就是常見(jiàn)的太陽(yáng)黑子11a周期的2倍。另外,準(zhǔn)4~8a的周期很可能與ENSO的活動(dòng)周期有關(guān)[20]。

(虛線表示95%置信區(qū)間)

3.2 太陽(yáng)黑子的多尺度特征和周期

對(duì)1950-2015年太陽(yáng)黑子數(shù)序列進(jìn)行EEMD分解(圖4),共得到6個(gè)IMF分量,其中IMF1~5為不同時(shí)間尺度下太陽(yáng)黑子數(shù)序列的周期項(xiàng),IMF6為整個(gè)序列的長(zhǎng)期變化趨勢(shì)項(xiàng)。由圖4可以看出,太陽(yáng)黑子數(shù)的變化過(guò)程較為平穩(wěn),具有明顯的周期性特征。在IMF1分量中,除1950-1960年這段時(shí)間外,其余時(shí)間其振幅都較為平穩(wěn)。IMF2~5這4個(gè)分量具有非常明顯的周期變化特征,各個(gè)時(shí)間尺度下的變化過(guò)程都具有類似正弦波的結(jié)構(gòu)。IMF6趨勢(shì)項(xiàng)的變化過(guò)程顯示,1950-1970年為太陽(yáng)黑子數(shù)上升階段,1970-2015年為下降階段。

(圖中縱坐標(biāo)為相對(duì)值)

進(jìn)一步的多窗譜分析(MTM)表明(圖5),太陽(yáng)黑子數(shù)序列具有顯著的36、11、5.4和3.2a周期(通過(guò)α=0.05顯著性檢驗(yàn))。其中太陽(yáng)黑子數(shù)的11a周期表現(xiàn)最為明顯,其譜值明顯大于其他幾個(gè)周期。

(虛線表示95%置信區(qū)間)

3.3 流量與太陽(yáng)黑子交叉小波分析

現(xiàn)代研究表明[21-22],太陽(yáng)活動(dòng)異常對(duì)地球氣候系統(tǒng)、降水、流域流量的影響主要是通過(guò)影響環(huán)流和季風(fēng)強(qiáng)度實(shí)現(xiàn)的。如過(guò)去120年以來(lái)南亞和東亞地區(qū)太陽(yáng)輻射的增強(qiáng)通過(guò)影響沃克環(huán)流和哈德來(lái)環(huán)流,最終導(dǎo)致區(qū)域降水和流域流量的增加。因此可以通過(guò)太陽(yáng)活動(dòng)和流量之間的關(guān)系,進(jìn)一步揭示二者在不同尺度下的相互影響過(guò)程。

通過(guò)前文的分析可知,贛江流量序列與太陽(yáng)黑子數(shù)序列的相關(guān)系數(shù)僅為-0.118,且未能通過(guò)α=0.1的顯著性檢驗(yàn),不能直接利用未經(jīng)處理的原始序列進(jìn)行二者之間的相關(guān)分析。但是,經(jīng)過(guò)EEMD分解后,各IMF分量的平穩(wěn)性有了很大程度的改善,周期性特征也更加明顯。與此同時(shí),流量和太陽(yáng)黑子的多窗譜分析(MTM)表明,二者都具有明顯的太陽(yáng)活動(dòng)周期,因此可以選取具有與太陽(yáng)活動(dòng)周期相關(guān)的IMF分量進(jìn)行對(duì)比分析。以此為基礎(chǔ),經(jīng)過(guò)交叉小波分析得到贛江流量和太陽(yáng)黑子數(shù)之間的交叉小波能量譜和小波凝聚譜。

圖6為1950-2015年贛江流量序列IMF2分量與太陽(yáng)黑子數(shù)序列的IMF2分量對(duì)比。由圖6可知,贛江流量變化與太陽(yáng)黑子數(shù)序列之間的相關(guān)關(guān)系在不同時(shí)間段內(nèi)是不一樣的。其中,1950-1960年呈負(fù)相關(guān)、1961-1973年呈正相關(guān)、1974-1984年呈負(fù)相關(guān)、1985-1995年呈正相關(guān)、1996-2007年為負(fù)相關(guān)(部分時(shí)段相差1/4個(gè)周期)、2008-2015年為正相關(guān)(部分時(shí)段相差1/4個(gè)周期)。近66年來(lái),大體上二者以11a為周期,在負(fù)相關(guān)和正相關(guān)之間交替變換。

(圖中縱坐標(biāo)為相對(duì)值)

圖7為1950-2015年贛江流量與太陽(yáng)黑子數(shù)序列的交叉小波譜和小波凝聚譜,交叉小波能量譜中顏色越偏紅表示能量譜密度值越大。由圖7(a)可以觀察到在11a周期附近贛江流量與太陽(yáng)黑子數(shù)的相關(guān)關(guān)系最為顯著,其中黑色實(shí)線所圍合區(qū)域表示通過(guò)了95%置信度檢驗(yàn),灰色陰影區(qū)則表示可能會(huì)受到邊界效應(yīng)的影響。圖中→表示二者之間變化的位相一致(正相關(guān)),←表示位相相反(負(fù)相關(guān)),↑和↓則表示二者之間提前或落后1/4個(gè)周期。在圖7(a)中的11a周期尺度上,1950-2015年二者之間由負(fù)相關(guān)逐漸順時(shí)針轉(zhuǎn)為相差1/4個(gè)周期,再轉(zhuǎn)為正相關(guān)。小波凝聚譜(圖7(b))中相關(guān)系數(shù)較高的區(qū)域與交叉小波能量譜(圖7(a))中基本一致,也是在11a尺度上達(dá)到最大值。其中1950-2000年,二者的相關(guān)系數(shù)在0.6以上,2001-2015年二者的相關(guān)系數(shù)也在0.4以上。小波凝聚譜可以發(fā)現(xiàn)在32a周期尺度上,二者的相關(guān)系數(shù)值在整個(gè)研究時(shí)段都在0.8以上,但是考慮到本研究的長(zhǎng)度只有66a,而且由于邊界效應(yīng)的存在,這一相關(guān)關(guān)系可能并不具有實(shí)際意義。

圖7 1950-2015年贛江流量與太陽(yáng)黑子數(shù)序列的交叉小波譜和小波凝聚譜

4 結(jié) 論

本文基于贛江外洲水文站1950-2015年的長(zhǎng)時(shí)間流量序列和同期太陽(yáng)黑子數(shù)序列,結(jié)合EEMD及CWT綜合分析方法對(duì)流量與太陽(yáng)活動(dòng)之間的多尺度相關(guān)關(guān)系進(jìn)行了分析。研究結(jié)果表明:

(1)EEMD-CWT綜合分析可以有效利用EEMD方法將原始信號(hào)進(jìn)行分解,得到具有平穩(wěn)性的IMF分量,并借助多窗譜分析(MTM)識(shí)別出信號(hào)中的不同周期,選擇具有相同物理意義的贛江流量序列IMF分量和太陽(yáng)黑子IMF分量,再利用交叉小波變換(CWT)得到二者在時(shí)域和頻域上的相關(guān)關(guān)系變化過(guò)程。

(2)1950-2015年贛江流量序列的EEMD分解結(jié)果表明,該序列具有明顯的非平穩(wěn)性,其中IMF1分量的波動(dòng)振幅最大,頻率最高,代表了原始序列中的主要變化因素,是水文序列非線性和非平穩(wěn)特征的體現(xiàn)。

(3)1950-2015年太陽(yáng)黑子數(shù)序列的EEMD分解結(jié)果表明,太陽(yáng)黑子數(shù)的變化過(guò)程較為平穩(wěn),具有明顯的周期性特征。趨勢(shì)項(xiàng)的變化過(guò)程顯示,1950-1970年為太陽(yáng)黑子數(shù)上升階段,1970-2015年為下降階段。

(4)多窗譜分析(MTM)表明,近66年來(lái)贛江流量序列具有22、8、4.5和3.7a周期,太陽(yáng)黑子數(shù)序列具有36、11、5.4和3.2a周期。

(5)交叉小波分析(CWT)表明,在11a周期尺度上,1950-2015年贛江流量與太陽(yáng)黑子數(shù)之間由負(fù)相關(guān)逐漸順時(shí)針轉(zhuǎn)為相差1/4個(gè)周期,再轉(zhuǎn)為正相關(guān)。其中1950-2000年,二者的相關(guān)系數(shù)在0.6以上,2001-2015年相關(guān)系數(shù)在0.4以上。

[1]GATESWL.Anoverviewoftheresultsoftheatmosphericmodelintercomparisonproject(AMIPI)[J].BulletinoftheAmericanMeteorologicalSociety,1999,80(2):29-55.

[2] 夏 軍,談 戈. 全球變化與水文科學(xué)新的進(jìn)展與挑戰(zhàn)[J]. 資源科學(xué),2002,24(3):1-7.

[3] 范垂仁,夏 軍,張利平,等. 中國(guó)水旱災(zāi)害長(zhǎng)期預(yù)報(bào)[M]. 北京:中國(guó)水利水電出版社,2008.

[4] 張?zhí)m生,方修琦,任國(guó)玉. 全球變化[M]. 北京:高等教育出版社,2000.

[5]NEFFU,BURNSSJ,MANGINIA,etal.Strongcoherencebetweensolarvariabilityandthemonsooninomanbetween9and6KyrAgo[J].Nature,2001,411(6835):290-293.

[6]ZHANGP,CHENGH,EDWARDSRL,etal.Atestofclimate,sun,andculturerelationshipsfroman1810-Yearchinesecaverecord[J].Science,2008,322(5903):940-942.

[7] 竇睿音,延軍平. 關(guān)中平原太陽(yáng)黑子活動(dòng)周期與旱澇災(zāi)害的相關(guān)性分析[J]. 干旱區(qū)資源與環(huán)境,2013,27(8):76-82.

[8]CHENF,YUANYJ,WEIWS,etal.Tree-ringrecordedhydroclimaticchangeinTienshanMountainsduringthepast500years[J].QuaternaryInternational,2015,358:35-41.

[9]YEHJ,SHIEHJ,HUANGNE.Complementaryensembleempiricalmodedecomposition:anovelnoiseenhanceddataanalysismethod[J].AdvancesinAdaptiveDataAnalysis,2010,2(2):135-156.

[10]WUZhaohua,HUANGNE.EnsembleEmpiricalModeDecomposition:Anoise-assisteddataanalysismethod[J].AdvancesinAdaptiveDataAnalysis,2009,1(1):1-41.

[11]GRINSTEDA,MOOREJC,JEVREJEVAS.Applicationofthecrosswavelettransformandwaveletcoherencetogeophysicaltimeseries[J].NonlinearProcessesinGeophysics,2004,11(5/6):561-566.

[12] 顧朝軍,穆興民,高 鵬,等. 贛江流域徑流量和輸沙量的變化過(guò)程及其對(duì)人類活動(dòng)的響應(yīng)[J]. 泥沙研究,2016,41(3):38-44.

[13]JIF,WUZ,HUANGJ,etal.Evolutionoflandsurfaceairtemperaturetrend[J].NatureClimateChange,2014,4(6):462-466.

[14]CHENYangkang,MAJitao.RandomnoiseattenuationbyF-Xempiricalmodedecompositionpredictivefiltering[J].Geophysics,2013,79(3):81-91.

[15]BREAKERLC,LOORHR,CARROLLD.Trendsinseasurfacetemperatureoffthecoastofecuadorandthemajorprocessesthatcontributetothem[J].JournalofMarineSystems,2016,164:151-164.

[16]WUZ,HUANGNE,CHENX.Themulti-dimensionalensembleempiricalmodedecompositionmethod[J].AdvancesinAdaptiveDataAnalysi,2009,1(3):339-372.

[17]MANNME,LEESJM.Robustestimationofbackgroundnoiseandsignaldetectioninclimatictimeseries[J].ClimaticChange,1996,33(3):409-445.

[18]HUANGNE,SHENZ,LONGSR,etal.Theempiricalmodedecompositionandthehilbertspectrumfornonlinearandnon-stationarytimeseriesanalysis[J].ProceedingsoftheRoyalSocietyofLondonA(Mathematical,PhysicalandEngineeringSciences),1998,454(1971):903-995.

[19]RINDD.Thesun'sroleinclimatevariations[J].Science,2002,296(5568):673-677.

[20] 王鐘睿,錢永甫. 江淮梅雨的多尺度特征及其與厄爾尼諾和大氣環(huán)流的聯(lián)系[J]. 大氣科學(xué)學(xué)報(bào). 2004,27(3):317-325.

[21] 葛全勝,劉路路,鄭景云,等. 過(guò)去千年太陽(yáng)活動(dòng)異常期的中國(guó)東部旱澇格局[J]. 地理學(xué)報(bào),2016,71(5):707-717.

[22]BHATTACHARYYAS,NARASIMHAR.RegionaldifferentiationinmultidecadalconnectionsbetweenIndianmonsoonrainfallandsolaractivity[J].JournalofGeophysicalResearchAtmospheres,2007,112(24):177-180.

猜你喜歡
分析
禽大腸桿菌病的分析、診斷和防治
隱蔽失效適航要求符合性驗(yàn)證分析
電力系統(tǒng)不平衡分析
電子制作(2018年18期)2018-11-14 01:48:24
電力系統(tǒng)及其自動(dòng)化發(fā)展趨勢(shì)分析
經(jīng)濟(jì)危機(jī)下的均衡與非均衡分析
對(duì)計(jì)劃生育必要性以及其貫徹實(shí)施的分析
GB/T 7714-2015 與GB/T 7714-2005對(duì)比分析
出版與印刷(2016年3期)2016-02-02 01:20:11
中西醫(yī)結(jié)合治療抑郁癥100例分析
偽造有價(jià)證券罪立法比較分析
在線教育與MOOC的比較分析
主站蜘蛛池模板: 欧美日韩精品综合在线一区| 午夜成人在线视频| 无码中字出轨中文人妻中文中| 日本道综合一本久久久88| jizz亚洲高清在线观看| 成人免费视频一区| 毛片在线看网站| 日韩人妻无码制服丝袜视频| 黄片一区二区三区| 亚洲精品无码日韩国产不卡| 欧美自慰一级看片免费| 十八禁美女裸体网站| 国产香蕉在线视频| 毛片视频网址| 国产日韩精品欧美一区灰| 广东一级毛片| 亚洲一区二区黄色| 国产精品观看视频免费完整版| 99久久国产综合精品2023| 日韩成人在线网站| 国产在线观看人成激情视频| 国产日韩丝袜一二三区| 一边摸一边做爽的视频17国产| 欧美亚洲国产日韩电影在线| 日本精品视频| m男亚洲一区中文字幕| 亚洲美女一区| 九色国产在线| 国产激爽大片在线播放| 人妻91无码色偷偷色噜噜噜| 国产精品欧美激情| 久久精品国产国语对白| 国产亚洲精品97在线观看| 一本大道香蕉高清久久| 色综合网址| 婷婷激情五月网| 日韩无码视频播放| 黄色网页在线观看| 国产成在线观看免费视频| 波多野结衣中文字幕一区| 亚洲日韩精品综合在线一区二区| 美女免费黄网站| 成人福利在线视频| 老司国产精品视频91| 国产一区二区三区免费观看| 婷婷综合缴情亚洲五月伊| 亚洲精品国产自在现线最新| 欧美成人在线免费| 制服丝袜一区| 亚洲不卡网| 91精品啪在线观看国产| 毛片视频网址| 色老二精品视频在线观看| 久草国产在线观看| 午夜免费小视频| 国产又爽又黄无遮挡免费观看 | 精品国产免费观看| 97久久超碰极品视觉盛宴| 亚洲成人精品在线| 久久精品国产精品青草app| 呦女精品网站| av尤物免费在线观看| 亚洲首页国产精品丝袜| 污视频日本| 在线色综合| 欧美第一页在线| 东京热一区二区三区无码视频| 国产va在线观看| 国产精品无码AV中文| 午夜毛片免费观看视频 | 亚洲精品视频免费看| 欧美一级黄片一区2区| 亚洲一区二区约美女探花| 国产精品护士| 亚洲综合香蕉| 无码人妻热线精品视频| 亚洲自拍另类| 欧美三级日韩三级| 亚洲自拍另类| 毛片视频网址| 久久频这里精品99香蕉久网址| 在线精品视频成人网|