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

磴口縣地下水埋深時空變化特征

2017-06-09 22:01:35李寧岳德鵬于強張啟斌馬歡
南水北調與水利科技 2017年3期
關鍵詞:模態

李寧+岳德鵬+于強+張啟斌+馬歡

摘要:選取荒漠綠洲區磴口縣1988年-2013年17個觀測站逐月水位埋深數據,運用kernel K-means及經驗模態分解(EMD)方法,探索26年來研究區地下水埋深時空變化特征。結果表明:17個測站分為三個聚類中心,第一聚類中心包括6個測站,地下水平均埋深最大。第二聚類中心包括4個測站,地下水平均埋深次之。第三聚類中心包括7個測站,地下水平均埋深最小;26年來第一和第二聚類中心地下水埋深呈增大趨勢,增大幅度分別為0.014 m、0.26 m。第三聚類中心地下水埋深呈減小趨勢,減小幅度為0.08 m;三個聚類中心地下水埋深年內變化趨勢基本相同。

關鍵詞:經驗模態分解;kernel K-means;磴口縣;地下水埋深;時空變化

中圖分類號:P642 文獻標識碼:A 文章編號:1672-1683(2017)03-0049-06

Abstract:This paper selected the monthly observed groundwater depth data during 1988-2013 from 17 monitoring wells at Dengkou County in the desert oasis region,and explored the temporal and spatial variation characteristics of the groundwater depth in study area during the 26 years by using kernel K-means and Empirical Mode Decomposition method.The results indicated the following findings:Firstly,17 monitoring wells can be divided into 3 clusters.The first cluster contains 6 monitoring wells with the largest average depth of groundwater.The second cluster contains 4 monitoring wells with the second largest average depth of groundwater.The third cluster contains 7 monitoring wells with the smallest average depth of groundwater.Secondly,during the past 26 years,the groundwater depths of the first and second clusters tend to increase,respectively up by 0.014 m and 0.26 m.The groundwater depth of the third cluster tends to decrease,down by 0.08 m.Thirdly,the groundwater depths of the three clusters basically show the same annual variation tendency.

Key words:EMD;kernel K-means;Dengkou County;groundwater depth;temporal and spatial variation

干旱地區水資源是制約區域經濟發展的關鍵因素。干旱地區的水資源特性使得地下水成為干旱地區極為重要的水源之一,如何合理、高效地利用地下水資源成為解決干旱區生態、環境問題的核心所在[1]。研究區磴口縣位于烏蘭布和東北緣,屬于典型的荒漠綠洲交錯區,生態環境十分脆弱,黃河水側滲使得該區地下水資源較為豐富,地下水資源已經成為維持該區經濟、社會以及綠洲可持續發展的重要因素[2]。因此掌握研究區地下水埋深時空變化規律、變化趨勢和變化幅度,進而能對當地地下水資源進行合理利用、區域開發戰略制定及生態環境保護提供科學支撐。

地下水動態變化是地下水對各種自然和人為影響因素的綜合響應,是地下水水資源量變的最直接反映[3]。近年來諸多學者對地下水動態變化研究較多,如徐永亮[4]對額濟納三角洲年內地下水水位進行動態分類,分析不同波動類型的地下水位動態變化特征,趙潔[5]分析了過去20年黑河中游地下水空間變異規律,楊懷德[6]對民勤綠洲地下水年際尺度埋深變化進行趨勢性分析,并探究了各相關因子與地下水埋深的敏感程度,肖彩虹[7]對烏蘭布和沙漠東北部人工綠洲七年來不同灌溉區地下水水位時空動態變化進行了研究,但對研究區磴口縣地下水動態變化研究較少且時間尺度短。

以往在研究地下水埋深動態變化特征時,大多采用傳統的描述性統計分析,而地下水埋深變化是一個非線性、非平穩的過程,伴隨著各種尺度的震蕩,由于傳統方法本身的局限,并不能真實有效地提取出地下水埋深變化的自然變率[8]。經驗模態分解(EMD)是一種新的非線性、非平穩時間序列分析方法。它與其他的時頻分解方法(如小波分析)相比,完全擺脫了傅里葉變化的束縛,可根據時間序列局部時變特征進行自適應地時頻分解,去掉持久穩固信號的噪聲,得到極高的時頻分辨率,非常適合對非平穩、非線性時間序列的分析[9-10]。

本文基于荒漠綠洲區磴口縣1988年-2013年17個長觀孔逐月觀測水位埋深數據,采用kernel K-means聚類方法和經驗模態分解法(EMD),分析了26年來研究區地下水埋深動態變化規律,以期為該區水資源與社會、經濟和生態的可持續發展提供參考依據。

1 研究區域及資料來源

研究區磴口縣地處烏蘭布和沙漠與河套平原的結合部,地理坐標為106°9′-107°10′E,40°9′-40°57′N,面積4 166.6 km2,南與鄂爾多斯市杭景旗及鄂托克旗隔河相望,西同阿拉善盟阿拉善左旗毗鄰,北與巴彥淖爾市烏拉特后旗相連,東與巴彥淖爾市臨河區接壤[11-12]。該區域屬溫帶大陸性季風氣候,干旱少雨、風大沙多,盛行西南風,受高空西風環流控制,平均相對濕度47%,多年平均降水量為142.7 mm,沙區不足100 mm,平均蒸發量2 372.1 mm,為降水的16.6倍,水資源豐富,地下水資源儲量5.258億m3 ,大小湖泊46處,水域面積3.61萬畝。磴口縣地處河套黃溉上游,攔河閘控制著整個河套的灌溉,因此引黃灌溉較其他旗縣條件優越,縣內絕大多數耕地可引黃灌溉。黃河流經磴口縣52 km,年徑流量310億m3,年均流量在580~1 600 m/s之間,河套灌區水利大動脈總干渠及烏審干渠橫穿縣境而過,黃河水側滲豐富,同時由于古地理環境及黃河改道,使磴口縣地下水資源十分豐富,地下水埋深2~9 m,單井出水量80~120 m/h,地下水資源分布狀態較為穩定,埋深淺、極易開采[13-14]。

本研究地下水埋深數據和位置資料來自內蒙古磴口縣黃河管理局,包括1988年-2013年逐月地下水埋深數據和17個觀測站的空間坐標。

2 研究方法

2.1 kernel K-means方法

每一個測站點地下水埋深數據共記錄312個月,即有312個維度。通過探索17個測站的逐月觀測數據,發現一些站許多數據記錄交織在一起,并不總是高于或者低于其他站,這很難用一個線性的平面來區分。kernel K-means是一個很好的聚類方法,擅長分離高維非線性分離數據。對于一個給定的集群數,它能通過最大限度的減少集群內的平方和,并把它作為每個站和相應的集群中心之間的平方歐氏距離的總和,最終找到最佳的適當分區[15]。它類似于經典的k-均值,但所有的操作都在一個特征空間中進行。kernel K-means方法將原始數據映射到一個更高維的特征空間,其中的數據可以很容易地分離成線性。在一個核特征空間中,原始數據(x1,…,x n)成為(φ(xi ),…,φ(xN )),其中φ(·)是映射函數。然后,數據可以通過以下步驟進行聚類。

(1)初始化特征空間中的聚類中心mk(k=1,…,K)。

(2)計算每個站φ(xi )(i=1,…,N)與聚類中心的距離k,根據最近鄰原則將每個站分配到相應類中

(4)重復步驟(2)和(3),直到連續n次迭代E值穩定為止。

然而,映射函數φ(·)的顯式表達式很難在沒有足夠的先驗信息的情況下確定。核函數提供了有效的映射函數替代,并且可以直接用于在特征空間中提供內積,而不需要顯式映射函數[16]。在本研究中使用高斯徑向基函數,它是一個常用的核函數。對應的特征空間是具有無限維度的希爾伯特空間,其中原始維空間中的非線性分離數據可以使用線性分類器容易地分離:

式中:B表示總距離平方距離之間的總和;W是群內總和的平方距離;k是聚類中心數;n為測站的數量。 CH值大表示是良好的聚類,其中聚類中心間的不相似性應該很大,而聚類中心內的相異性應該很小。選擇由CH最大化的k值作為最終聚類中心數。我們使用統計軟件R2.15運行聚類方法和包kernlab。

2.2 EMD方法

經驗模態分解法(EMD)是一種自適應時間序列數據分析模型,它的提出是對傅里葉分析為基礎的線性和穩態譜分析的一個重大突破[17]。在EMD分析中,原始數據被分解成一系列的模式,而不需要事先知道。與傅立葉和小波分解相比,EMD有自己的優勢。傅立葉分析可以將數據轉換成具有不同頻率的正弦和余弦函數的組合,而小波分析需要小波函數。由于基礎函數諧波特性的限制,這些類型的分解包含了許多雜散分量。經驗模態分解法(EMD)將數據分解成若干本征模態函數(IMF)和一個殘余信號,殘留底層即是數據的趨勢[18]。

本征模態函數(IMF)需要滿足兩個條件:首先,極值的數目和在整個數據跨度過零點的數目之間的差應小于或等于1。第二,由局部極大值和局部最小值確定的包絡線所限定的包絡線的平均值是在任何點為零[19]。本研究使用篩選方法計算IMF和地下水埋深殘差。

步驟如下。

(1)初始化r(t)=x(t),i=0,k=1,其中x(t)是時間序列地下水埋深向量,計算局部極大值和時間序列r(t)的局部極小值。

(2)計算通過內插極大值上包絡線emax(t)和下包絡emin(t)的經內插的局部最小值。

式中:K是的本征模態函數的個數,r是原始數據的趨勢。第一本征模態函數(IMF1)是具有最高頻率的分量,所有本征模態函數分量進行統計學零假設測試[20]。在應用 EMD 方法時無法避免邊界問題,因為構成上下包絡線時,信號數據序列的兩端會出現發散現象,并且這種發散會不斷逐漸向內“污染”,使所得結果嚴重失真。本文在MATLAB軟件中運行Rilling[21]等人編寫的EMD代碼,采用鏡像對稱延伸方法進行包絡擬合。

3 結果分析

3.1 1988年-2013年地下水埋深變化統計特征

通過對研究區17點測站點地下水埋深數據進行描述性統計分析可得(表1),測得地下水埋深最小的是巴11測站點,為0.09 m;最大的是巴10測站點,為4.30 m。地下水埋深平均值巴4測站點最大,為2.65 m;巴6測站點最小,為0.94 m。地下水埋深變化幅度最大的巴10測站點,變化幅度為3.81 m;最小的為巴6測站點,變化幅度為1.45 m。根據變異系數(Cv)的大小可知,17個測站點都屬于中等變異性,其中巴9、巴11變異性最大,巴5、巴16變異性次之,巴1點變異性最小。

3.2 站點聚類的空間分布

利用kernel K-means方法對磴口縣17個地下水埋深觀測站進行聚類,從1到8不同的集群數來計算CH值,發現其達到最大值時集群數為3,表明此時各類內差異最小,類間的差異最大。因此17個觀測站被劃分為3個聚類中心。第一聚類中心包含6個站:巴1,巴7,巴9,巴10,巴13,巴14。第二聚類中心包含4個站:巴4,巴8,巴15,巴17。第三聚類中心包括7個站:巴2、巴3、巴5、巴6、巴11、巴12、巴16。

第一個聚類中心中地下水埋深最大,主要分布于磴口縣北部及磴口縣城附近。

第二個聚類中心中地下水埋深次之,分布在磴口縣中部,在第一聚類中心和第三聚類中心之間。

第三個聚類中心中地下水埋深最小,分布在磴口縣東北部和西南部。

3.3 年際及年內變化趨勢

利用EMD方法對三個聚類中心進行趨勢性變化分析(圖3)。IMF1具有最高的頻率和最大的振幅,而IMF6具有最低的頻率和最小的振幅。當所有的周期從原始數據中去掉時,即得到整個時間跨度的趨勢。三個聚類中心變化范圍分別為2.2~2.214 m、 1.80~2.06 m、1.19~1.27 m。從三條趨勢線可以得出兩個重要的不同。一個是變化幅度不同。第二聚類中心變化幅度最大,為0.26 m;第三聚類中心次之,變化幅度為0.08 m,第一聚類中心變化幅度最小,僅為0.014 m,基本保持穩定。二是各個聚類中心趨勢線的不同。第一聚類中心單調遞增。第二聚類中心先緩慢增加而后趨于穩定。而第三聚類中心先增加,到1994年時達到最大,而后緩慢減小,到2008年之后,基本穩定不變。

對三個聚類中心26年來每個月地下水埋深觀測值取平均(圖4),從圖中可以非常明顯看出:綠洲區內三個聚類中心地下水埋深雖然不同,但年內變化特征基本相似:從2月到6月,地下水埋深逐漸減??;從6月一直持續到9月,地下水埋深開始逐漸增大,從9月到11月地下水埋深又逐漸減小,達到全年地下水埋深最??;從11月開始地下水埋深又開始逐漸增大,到第二年2月份左右達到全年最大。河套灌區年內有三次灌溉期,分別是夏灌(4月-6月,灌水3次)、秋灌(7月-9月,灌水3次)、秋澆(10月-11月,灌水1次)。灌溉回滲導致3月-6月、9月-11月地下水埋深變小,地下水蒸發、蒸騰以及人類的生活用水的增加使6月-9月地下水埋深逐漸變大。

3.4 時間趨勢分析

各觀測站地下水埋深數據經過EMD方法處理得到年變化趨勢圖(圖5),所有分解的趨勢圖都是非線性的。分析17個站的年變化趨勢圖,巴9是變化最大的站(2.3 m),其次是巴10站(1.3 m)、巴4站(1.2 m)。根據趨勢線的變化情況,分為四個類型(圖5)。第一類(I)是單調上升的。地下水埋深逐漸變大,或短時間內保持不變。這類型包括7個站(巴2、巴3、巴4、巴7、巴10、巴15、巴17)。第二類(II)單調遞減的。這類型地下水埋深逐漸變小,包括3個站(巴5、巴9、巴11)。第三類(III)先下降后上升再下降。地下水埋深先變小,而后逐漸變大最后再變小,包括5個站(巴1、巴6、巴8、巴16、巴14)。第四類(IV)先上升后下降。地下水埋深早期逐漸變大而后又變小,包括兩個站(巴12、巴13)。 4 結論與討論

(1)從三個聚類中心趨勢線分析可知,研究區地下水平均埋深大于1.8 m的測站地區,地下水埋深逐漸增大,地下水平均埋深小于1.8 m的測站地區,

地下水埋深從1988年-1993年緩慢變大而后逐漸減小,到2008左右保持穩定。近些年隨著磴口縣人口的不斷增長,加之不合理的利用開采地下水,造成地下水埋深有增大趨勢。由于當地政府加強對濕地的保護,使靠近納林湖等湖泊的測站點地下水埋深小、基本保持穩定。

(2)綠洲區地下水埋深年內變化特征: 地下水埋深變小期發生在2月-6月、9月- 11月,埋深變大期發生在6月-9月、11月-2月,其中最大值出現在2月份,最小值出現在11月份,造成這種變化的主要原因是河套地區一年三次的灌溉以及人類和動植物的活動。

(3)各個測站點年際變化趨勢呈現四種不同類型,分別為單調上升、單調遞減、先下降后上升再下降、先上升后下降。從空間上分析17個測站,地下水埋深逐漸增加的站大都位于人口密集的地區,比如巴10點位于巴彥淖爾鎮附近,離縣城最近,其地下水埋深增加幅度最大,其次是巴4點,位于納林套海農場。由于對濕地的保護,巴5、巴9、巴11測站地下水埋深均有不同程度的減小。

(4)磴口縣處于荒漠綠洲交錯帶,生態環境脆弱,一旦破壞難以恢復。降雨稀少,蒸發強烈,水資源大都由黃河側滲的地下水進行補給,因此合理利用地下水資源,進一步加強對濕地的保護成為重中之重。

參考文獻(References):

[1] 陳亞寧,李稚,范煜婷,等.西北干旱區氣候變化對水文水資源影響研究進展[J].地理學報,2014,69(9):1295-1304.(CHEN Ya-ning,LI Zhi,FAN Yu-ting,et al.Research progress on the impact of climate change on water resources in the arid region of Northwest China[J].Acta Geographica Sinica,2014,69(9):1295-1304.(in Chinese))

[2] 李慧卿,李慧勇,李清河,等.烏蘭布和沙漠東北緣地下水位對白刺生長的影響[J].林業科學,2011,47(11):25-30.(LI Hui-qing,LI Hui-yong,LI Qing-he,et al.Impact of groundwater depth on growth of Nitraria tangutorum on the northeast of ulanbuh desert[J].Scientia Silvae Sinicae,2011,47(11):25-30.(in Chinese))

[3] 唐克旺,唐蘊,李原園,等.地下水功能區劃體系及其應用[J].水利學報,2012,43(11):1349-1356.(TANG Ke-wang,TANG Yun,LI Yuan-yuan,et al.Technical system of functional division for groundwater and its application in China[J].Journal of Hydraulic Engineering,2012,43(11):1349-1356.(in Chinese))

[4] 徐永亮,于靜潔,王平,等.額濟納三角洲地下水位年內動態變化特征分類分析[J].干旱區資源與環境,2013,27 (4):135-140.(XU Yong-liang,YU Jing-jie,WANG Ping,et al.Classification of annual dynamic of groundwater depth in Ejina delta[J].Journal of Arid Land Resources and Environment,2013,27 (4):135-140.(in Chinese))DOI:10.13448/j.cnki.jalre.2013.04.035

[5] 趙潔,徐宗學,周劍.黑河中游過去20年地下水位空間變異性分析[J].干旱區資源與環境,2011,25(8):172-178.(ZHAO Jie,XU Zong-xue,ZHOU Jian.Spatial variation of ground water in middle stream of The Heihe River Basin during past two decades[J].Journal of Arid Land Resources and Environment,2011,25(8):172-178.(in Chinese))DOI:10.13448/j.cnki.jalre.2011.08.027

[6] 楊懷德,馮起,郭小燕.1999-2013年民勤綠洲地下水埋深年際變化動態及影響因素[J].中國沙漠,2017,37(3):1-9.(YANG Huai-de,FENG Qi,GUO Xiao-yan.Variation of groundwater depth and its influence factors in the Minqin oasis in 1999-2013[J].Journal of Desert Research,2017,37(3):1-9.(in Chinese))DOI:10.7522/j.issn.1000-694X.2016.00077

[7] 肖彩虹,郝玉光,辛智鳴,等.烏蘭布和沙漠東北部人工綠洲地下水位動態變化研究[J].防護林科技,2013,12(12):1-3.(Xiao Cai-hong,Hao yu-guang,Xin Zhi-ming,et al.Dynamic Changes of Groundwater Level at Artificial Oasis in Northeastern Ulan Buh Desert[J].Protection Forest Science and Technology,2013,12(12):1-3(in Chinese))DOI:10.13601/j.issn.10055215.2013.12.013

[8] 張殷欽,胡偉,劉俊民.基于R/S分析法的地下水位動態變化趨勢分析[J].中南大學學報:自然科學版,2012,43(12):4912-4916.(ZHANG Yin-qin,HU Wei,LIU Jun-min.Groundwater level regime variation trend on basis of rescaled range analysis[J].Journal of Central South University:Science and Technology,2012,43(12):4912-4916.(in Chinese))

[9] 孫嫻,林振山.經驗模態分解下中國氣溫變化趨勢的區域特征[J].地理學報,2007,62(11):1132-1141.(SUN Xian,LIN Zhen-shan.The regional features of temperature variation trends over China by empirical mode decomposition method[J].Acta Geographica Sinica,2007,62(11):1132-1141.(in Chinese))

[10] Huang J,Xia Z,Zhou S,et al.Spatio-temporal analysis of malaria vectors in national malaria surveillance sites in China[J].Parasite Vector,2015,8(1):146.DOI:10.1186/s13071-015-0741-5

[11] 劉芳,郝玉光,張景波,等.磴口縣土地景觀變化遙感監測研究[J].內蒙古農業大學學報:自然科學版,2009(3):112-116.(LIU Fang,HAO Yu-guang,ZHANG Jing-bo,et al.Study on landscape changes for various land utilization types by using remote sensing in Dengkou county[J].Journal of Inner Mongolia Agricultural University:Science and Technology,2009(3):112-116.(in Chinese))

[12] 高儒,溫俊,周敏,等.淺談磴口縣水資源現狀及保護利用對策[J].內蒙古水利,2013(3):44-45.(GAO Ru,WEN Jun,ZHOU Min,et al.Protection and utilization of water resources in Dengkou county[J].Inner Mongolia Water Resources,2013(3):44-45.(in Chinese))

[13] 于強,岳德鵬,郝玉光,等.磴口縣荒漠綠洲景觀特征及時空變化研究[J].干旱區資源與環境,2016,30(4):178-183.(YU Qiang,YUE De-peng,HAO Yu-guang,et al.The spatial and temporal variations of oasis desert landscape in Dengkou County[J].Journal of Arid Land Resources and Environment,2016,30(4):178-183(in Chinese))DOI:10.13448 /j.cnki.jalre.2016.133

[14] 磴口編纂委員會.磴口縣志(1988-2010)[M].內蒙古文化出版社,2014.(Dengkou Compilation Committee.Dengkou county annals(1988-2010)[M].Hulun Buir:Inner Mongolia Culture Press,2014.(in Chinese))

[15] Hu M G,Jia L,Wang J F,et al.Spatial and temporal characteristics of particulate matter in Beijing,China using the Empirical Mode Decomposition method[J].Science of the Total Environment,2013,458-460:70-80.DOI:10.1016/j.scitotenv.2013.04.005

[16] 趙睿.擴散距離下的kernel K-means聚類算法的改進[D].清華大學,2011.(ZHAO Rui.Clustering with kernel K-means and diffusion distance[D].Tsinghua University,2011.(in Chinese))

[17] 孫小龍,劉耀煒,晏銳.經驗模態分解法在地下水資料處理中的應用[J].大地測量與地球動力學,2011(2):80-83.(SUN Xiao-long,LIU Yao-wei,YAN Rui.Application of empirical mode decompostion method to ground water data[J].Journal of Geodesy and Geodynamics,2011(2):80-83(in Chinese))DOI:10.14075/ j.jgg.2011.02.034

[18] 劉慧婷,倪志偉,李建洋.經驗模態分解方法及其實現[J].計算機工程與應用,2006,42(32):44-47.(LIU Hui-ting,NI Zhi-wei,LI Jian-yang.Empirical mode decomposition method and its implementation[J].Computer Engineering and Applications,2006,42(32):44-47.(in Chinese))

[19] 王婷.EMD算法研究及其在信號去噪中的應用[D].哈爾濱:哈爾濱工程大學,2010.(WANG Ting.Research on EMD algorithm and its application in signal denoising[D].Harbin:Harbin Engineering University,2010.(in Chinese))

[20] Wang D,Tang L,Liu H,et al.Application of structural equation models for elucidating the ecological drivers of anopheles sinensis in the three gorges reservoir[J].PLoS One,2013,8(7):e68766.DOI:10.1371/journal.pone.0068766

[21] Rilling G,Flandrin P,Goncalves P.On empirical mode decomposition and its algorithms[C].Proceedings of IEEE-EURASIP workshop on nonlinear signal and image processing,2003.

猜你喜歡
模態
基于BERT-VGG16的多模態情感分析模型
跨模態通信理論及關鍵技術初探
一種新的基于模態信息的梁結構損傷識別方法
工程與建設(2019年1期)2019-09-03 01:12:12
多跨彈性支撐Timoshenko梁的模態分析
車輛CAE分析中自由模態和約束模態的應用與對比
國內多模態教學研究回顧與展望
基于HHT和Prony算法的電力系統低頻振蕩模態識別
由單個模態構造對稱簡支梁的抗彎剛度
計算物理(2014年2期)2014-03-11 17:01:39
利用源強聲輻射模態識別噪聲源
日版《午夜兇鈴》多模態隱喻的認知研究
電影新作(2014年1期)2014-02-27 09:07:36
主站蜘蛛池模板: 狂欢视频在线观看不卡| 成年女人a毛片免费视频| 国产亚洲精久久久久久无码AV| 丁香婷婷综合激情| 91亚洲免费| 亚洲综合色婷婷| 国产一区二区三区在线观看视频 | 国产欧美日韩综合在线第一| 亚洲欧美成人网| 青青草原国产免费av观看| 精品无码日韩国产不卡av| 久久semm亚洲国产| 亚洲天堂视频在线观看免费| 久久91精品牛牛| 国产色婷婷| 91青青草视频在线观看的| 在线视频一区二区三区不卡| 欧美97色| 精品免费在线视频| 精品亚洲欧美中文字幕在线看| 超碰aⅴ人人做人人爽欧美 | 热这里只有精品国产热门精品| 国产高清在线观看91精品| 国产成人精品一区二区秒拍1o| 欧美日本激情| 91麻豆国产视频| 久久不卡精品| 国产麻豆精品手机在线观看| 亚洲视频一区| 亚洲天堂精品在线| 国产丝袜啪啪| 精品福利视频网| 婷婷综合色| 波多野结衣一二三| 青草精品视频| 成人无码区免费视频网站蜜臀| 超清人妻系列无码专区| 国产成人综合在线观看| 亚洲日本www| 欧美在线视频a| 香蕉蕉亚亚洲aav综合| 国产情侣一区| 日韩小视频在线观看| 国产免费一级精品视频 | 亚洲精品无码抽插日韩| 久久久波多野结衣av一区二区| 欧美成人A视频| 国产自无码视频在线观看| 亚洲欧美在线精品一区二区| 色综合久久久久8天国| 亚洲国产第一区二区香蕉| 99在线观看精品视频| 一级毛片不卡片免费观看| 香蕉视频在线精品| 在线观看网站国产| 在线观看视频一区二区| 日韩美一区二区| 久久人人97超碰人人澡爱香蕉| 精品国产自在现线看久久| 亚洲中文字幕在线精品一区| 久久久91人妻无码精品蜜桃HD| 亚洲人成成无码网WWW| 性喷潮久久久久久久久| 欧美成人午夜视频免看| 她的性爱视频| 亚洲日韩久久综合中文字幕| 91网址在线播放| 欧亚日韩Av| 青青草国产在线视频| 国产精品毛片一区视频播| 欧美人与性动交a欧美精品| 成人一级免费视频| 午夜天堂视频| 欧美午夜视频在线| 国产精品第| 99激情网| 中文字幕无码中文字幕有码在线 | 中文一级毛片| 欧美日韩在线第一页| 亚洲欧美另类中文字幕| 国产区成人精品视频| 小说 亚洲 无码 精品|