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

基于GIS的淮南市土壤Cu含量克里格插值方法比較研究

2015-11-12 13:34:12郭旻欣喻根王曉輝等
湖北農業科學 2015年20期
關鍵詞:模型

郭旻欣+喻根+王曉輝+等

x摘要:利用ArcGIS地統計學中克里格空間插值方法,對分布類型為對數正態分布、具有二階多項式趨勢的淮南市土壤重金屬Cu含量的空間結構及分布特征進行了研究。基于經過對數變換的Cu含量分析數據,考慮是否進行分析數據的趨勢移除,綜合比較了三類克里格法(普通克里格、簡單克里格和泛克里格)共6個具體插值方法的分析結果。結果表明,土壤Cu元素含量數據經過二階趨勢移除,采用簡單克里格方法插值效果最優。在11種可選的半變異函數模型中,孔穴效應模型插值效果優于高斯、指數等其他模型。研究區土壤Cu含量屬中等程度的空間相關,說明Cu含量空間變異受隨機性因素和結構性因素共同影響。淮南市土壤Cu含量除受土壤內在因素影響外,可能來源于人類的采礦活動和對化肥、有機肥和農藥的過度使用。

關鍵詞:土壤;銅;地統計學;空間結構;克里格插值;淮南市

中圖分類號:X53 文獻標識碼:A 文章編號:0439-8114(2015)20-4993-06

DOI:10.14088/j.cnki.issn0439-8114.2015.20.018

Study on Comparison of Kriging Interpolation Methods of Soil Cu Content in Huainan City Based on GIS

GUO Min-xin1,YU Gen1,WANG Xiao-hui2,LIU Yu-juan1,GU Xiao-yang1

(1.School of Resources and Environmental Engineering, Hefei University of Technology, Hefei 230009, China;

2.Research Academy of Environmental Sciences of Anhui, Hefei 230009, China)

Abstract: To identify the spatial structure and distribution of soil heavy metal Cu within Huainan city, the Cu data with lognormal distribution, second-order polynomial trend were analyzed and processed using Kriging Geostatistical Analyst tool in ArcGIS. Based on the Cu content data of logarithmic transformation and considering data trend removal issues, a total of 6 interpolation methods of three categories of Kriging (ordinary, simple and universal Krigings) were compared. The results showed that using data with second trend removal, simple Kriging interpolation method provided optimal results. Among 11 kinds of available semivariogram models, hole effect model performed best. Cu content in soil was auto-correlated to a moderate degree, which indicated that the concentration of Cu in the soils was influenced by the random and structural factors. Apart from soil intrinsic factor, soil heavy metal Cu content of Huainan city may be originated from human mining activities and the excessive use of fertilizers and pesticides.

Key words: soil; copper; geostatistics; spatial structure; Kriging; Huainan city

土壤是人類食物生長的地方,同時也是人類活動的基本場所,土壤環境質量與人類健康密切相關。淮南市是中國億噸煤基地、華東火電基地和煤化工基地的“三大基地”,經濟以重工業為主,煤礦業為支柱產業。近年來,由于煤炭生產和工礦企業的迅速發展,使得該區的土壤環境遭受了不同程度的污染[1]。因此,準確模擬土壤重金屬的空間分布,合理評價土壤重金屬的污染程度,及時有效地采取防護、修復措施,對減少重金屬污染,保障人們的身體健康具有重要意義。

對于淮南市重金屬污染狀況已有學者進行一定研究。崔龍鵬等[2]對有百年開采歷史的淮南礦區土壤污染的調查發現,采礦歷史越長礦井周圍土壤中重金屬的富集越高。李海霞等[3]對淮南某廢棄礦區污染場的土壤重金屬污染潛在生態風險進行了評價。王興明等[4]對淮南新莊孜煤矸石山附近土壤中4種重金屬元素(Zn、Pb、Cd和Cu)分布特征進行了研究。王曉輝等[5]對淮南礦區6種重金屬空間分布及污染來源進行了討論。這些研究范圍主要集中在淮南市某些局部礦區,其空間尺度很小。研究重點在淮南礦區土壤重金屬含量的空間分布及污染來源,對于影響土壤重金屬含量分布的半變異函數理論模型、最優的空間插值方法沒有進行具體研究。因此,很有必要對淮南全市范圍內土壤重金屬含量的空間結構及分布特征進行系統的研究。endprint

國外學者[6-8]對地統計學中的插值方法進行了比較研究,發現距離反比權重法(IDW)比三次樣條插值效果要好,克里格法比距離反比權重法好。在國內應用于土壤重金屬空間分布研究的方法多為普通克里格法[9-14],對簡單克里格和泛克里格方法研究未見較多的報道。針對土壤重金屬的空間結構研究,以GS+軟件進行半變異函數理論模型擬合居多[10,13,14],這種方法僅僅比較了高斯、指數等4種常見半變異函數模型,不能保證重金屬含量空間插值的準確性。因此,本研究以淮南市土壤重金屬Cu為例,應用3種克里格插值方法,對11種半變異函數模型插值結果進行分析比較,選取最優半變異函數模型和插值方法,對淮南市土壤重金屬Cu的空間結構和空間分布特征進行研究,對認識淮南市土壤污染特征和環境管理具有一定的意義。

1 材料與方法

1.1 研究區概況

淮南市位于安徽省中北部,淮河之濱,全市總面積約2 596 km2。位于東經116°21′21″~117°11′59″,北緯32°32′45″~33°0′24″,東與滁州市毗鄰,南與合肥市接壤,西南與六安市相連,西北與阜陽市、亳州市交界,東北與蚌埠市相交。淮河由西向東橫穿淮南市,淮河北岸為地勢平坦的淮北平原,淮河南岸為丘陵地區。土壤類型復雜多樣,淮河以北的平原地區,基質是古河流沉積物,主要為砂礓黑土和黃土。淮河沿岸的灣地為潮土類土壤。灣地與丘陵之間的崗地,基質為下蜀系黃土,土壤主要為黃棕壤[1]。

1.2 樣品采集

按照《全國土壤污染狀況調查技術規定》的布點要求和編碼規則,對淮南市0~20 cm表層土壤進行8 km×8 km網格采樣。在50 m×50 m采樣區內以梅花形分點采樣后,將各分點樣品等質量混勻,用“四分法”棄取,保留相當于3 kg干土壤的土樣量。另一部分的樣品采集以礦井為中心由密漸疏向周圍放射狀布設采樣點,在矸石山堆放處根據具體情況適當加密布點,采集0~20 cm表層土壤樣,兩部分共布設采樣點200個(圖1)。Cu元素含量測定方法采用X射線熒光光譜法(XRF)。

1.3 基本原理

1.3.1 克里格法 克里格法最早是由法國地理數學家Matheron和南非礦山工程師Krige提出的。克里格方法是將被插值的要素當作一個區域化變量來對待,區域化變量隨所在區域位置的改變而連續地變化。因此,彼此離得近的點之間相關性比較強,距離遠的點之間相關性較弱[15]。克里格內插法是根據無偏估計和方差最小兩項原則來確定加權系數,其中關鍵的函數是半方差函數。在地統計分析中通過下面公式來對點進行預測:

Z(s0)=∑■■?姿iZ(si) (1)

式中,Z(s0)為預測點的數值,?姿i為第i個采樣點的權重,Z(si)為第i個采樣點的數值,n為估計預測點的采樣點個數。

為了使預測結果是無偏估計,N個采樣點的權重和需要為1。通過下面計算公式求出采樣點的權重。

?祝*?姿=g (2)

?酌11 … ?酌1N 1■ ?塤 ■ ■?酌N1 … ?酌NN 1 1 … 1 0*?姿1■?姿Nm =?酌10■?酌N01 (3)

式中,?祝為所有采樣點之間對應的半方差函數值組成的矩陣,λ為權重矩陣,g為采樣點和預測點對應的半方差函數值組成的向量。?酌ij為第i個和第j個采樣點之間對應的半方差函數值,λi為第i個采樣點的權重,m為約束預測無偏性變量,?酌i0為第i個采樣點和預測點對應的半方差函數值。將向量g乘以?祝的逆矩陣可以求出向量λ,根據公式(1)算出預測點的值[16]。

1.3.2 交叉驗證法 通常在比較各種地統計學插值方法時,主要是對其半方差函數模型進行檢驗,但到目前為止還沒有一個行之有效的直接檢驗方法。地統計學者多采用一種間接的方法,即將克里格插值方法與半方差函數模型相結合進行檢驗,這種檢驗方法被稱為交叉驗證法。設Z(xi)和Z*(xi)分別為實測值和預測值,?滓(xi)為點xi處的預測標準方差。則它們的平均預測誤差ME(Mean Error)、均方根誤差RMSE(Root-Mean-Square Error)、平均標準誤差ASE(Average Standard Error)、平均標準化預測誤差MSE(Mean Standardized Error)、均方根標準預測誤差RMSSE(Root-Mean-Square Standardized Error)可分別表示為:

ME=■■[Z(xi)-Z*(xi)] (1)

RMSE=■ (2)

ASE=■ (3)

MSE=■■■ (4)

RMSSE=■ (5)

根據以下4個原則選擇最優的克里格插值法及半變異函數模型:①預測是無偏的,則平均預測誤差ME應接近于零;②預測值盡可能地接近測量值,均方根誤差RMSE和平均標準化預測誤差MSE越小越好;③預測的不確定性是有效的,如果平均標準誤差ASE接近于均方根誤差RMSE,則表示正確地評價了預測的變異性。或者,如果預測標準誤差是有效的,那么均方根標準預測誤差RMSSE應接近于1;④預測值與測量值分布的擬合直線的斜率接近于1[16-19]。

1.4 數據處理

數據統計分析和正態分布檢驗采用SPSS 19.0軟件完成,半變異函數理論模型擬合和克里格空間插值采用ArcGIS10.0完成。

2 結果與分析

2.1 土壤重金屬Cu含量特征

2.1.1 Cu元素含量統計特征值 淮南市土壤重金屬Cu含量統計分析結果見表1。從偏度和峰度系數可知,土壤重金屬Cu基本服從對數正態分布(偏度接近于0,峰度接近于3)。淮南市土壤重金屬Cu的平均含量為46.55 mg/kg,顯著高于淮南市“十一五”背景值26.1 mg/kg,說明淮南市土壤重金屬Cu具有明顯的累積效應。變異系數反映了各采樣點之間的平均變異程度。表1中Cu的變異系數為101%,屬強變異性,說明土壤中Cu受人為活動干擾顯著[10,11,20]。2.1.2 Cu元素含量趨勢效應 Cu元素含量趨勢如圖2所示,X軸表示東西方向,Y軸表示南北方向,Z軸為采樣點Cu元素含量,按東西(綠線)和南北(藍線)兩個方向投影到與地圖平面正交的平面上。探索兩個方向存在的趨勢,如果兩條曲線平直,則說明全局趨勢不存在,如果擬合曲線可以用二次函數、指數函數等數學公式表示,則表明存在全局趨勢。由圖2所示,淮南市土壤重金屬Cu的含量分布具有明顯的趨勢效應,XOZ和YOZ平面上的曲線形狀均為拋物線型,表明Cu元素含量數據存在二階多項式趨勢[19]。endprint

2.2 土壤重金屬Cu的各種估值方法比較

2.2.1 克里格插值方法比較 克里格法對服從正態分布的數據插值精度更高,如果數據不服從正態分布會使半變異函數產生比例效應,抬高基臺值和塊金值。因為本試驗研究區內數據為偏態分布,為避免半變異函數分析產生比例效應,在進行變異函數分析前首先對原始數據進行了對數變換,使數據接近正態分布[9,21]。

基于經過對數變換后的Cu含量分析數據,假定半變異函數均擬合為高斯模型情況下,插值的趨勢效應分別選0階(無趨勢)和2階(二階多項式),綜合比較了3類克里格方法(普通克里格、簡單克里格和泛克里格)共6個具體插值方法的分析結果(表2)。由表2可知,Cu含量數據經過二階趨勢移除后的插值結果要優于0(無趨勢);對于3種克里格法,平均預測誤差ME最小、RMSSE最接近于1的是簡單克里格法,MSE最小以及ASE與RMSE最接近的是普通克里格法和簡單克里格法。綜合比較各種誤差大小和插值效果,二階趨勢移除的簡單克里格法比其余5種方法要好。

2.2.2 半變異函數模型的比較 為了比較不同模型的半方差函數模型擬合的精度情況,基于經過對數變換的Cu含量分析數據,假定在經過二階趨勢移除條件下,均采用簡單克里格插值方法,擬合的半方差函數理論模型分別為球形模型(Circular model)、球狀模型(Spherical model)、四球模型(Tetraspherical model)、五球模型(Pentaspherical model)、指數模型(Exponential model)、高斯模型(Gaussian model)、有理二次模型(Rational Quadratic model)、孔穴效應模型(Hole Effect model)、K貝塞爾模型(K Bessel model)、J貝塞爾模型(J Bessel model)和穩定模型(Stable model)等11種半變異函數模型。

由表3可知,在11種半變異函數模型中,MSE最接近于0、RMSSE最接近于1以及ASE與RMSE最接近的均是孔穴效應模型(Hole Effect model),孔穴效應模型的平均預測誤差ME值也接近于0,斜率接近于1。綜合考慮各種參數半變異函數模型選擇孔穴效應模型,平均誤差最小,插值效果最好。

2.3 土壤重金屬Cu含量空間結構特征

利用內插效果最優的孔穴效應模型對試驗半變異函數進行了擬合,重金屬Cu的理論模型對試驗半變異函數擬合程度比較好。塊金值與基臺值的比值是反映變量空間異質性的指標,表示區域化變量的空間相關程度,可以反映影響因素中結構性因素(自然因素)和隨機性因素(人為因素)的作用。對照1994年Cambardella等[22]的劃分方案,塊金值和基臺值的比值小于25%表明變量的空間變異以結構性變異為主,空間相關性強;大于75%表明空間變異以隨機性變異為主,空間相關性較弱;處于兩者之間說明空間相關性中等。土壤重金屬Cu元素含量擬合半變異函數模型的參數以及空間相關性分級見表4。

土壤重金屬Cu元素半變異函數的塊金值與基臺值比值為31.2%,表明Cu元素具有中等空間相關性,說明Cu元素含量的富集是由結構性因素和隨機性因素共同作用的結果。結構性因素如氣候、母質、地形、土壤類型等自然因素可以導致土壤元素有強的空間相關性,而隨機性因素如施肥、農藥使用、耕作措施、礦業開采、企業生產等各種人為活動使得土壤元素的空間相關性減弱,變異性增強[10]。

變程表示隨機變量在空間上的自相關尺度,它與觀測以及取樣尺度有關。在變程范圍內,變量才有空間自相關性。由Cu元素半變異函數圖面(圖3)可見,淮南市土壤Cu元素含量空間變異特征無方向性差異,其主變程和次變程均為4.054 km。Cu的變程較小,說明其受采礦、施肥和噴灑農藥等人為活動的隨機性因素的影響較大,從而導致其在相對較小范圍內存在相關關系。

2.4 土壤重金屬Cu含量空間分布特征

根據擬合效果最好的孔穴效應理論模型,采用簡單克里格法,對經過對數變換和二階趨勢移除后的Cu含量數據進行插值,進而得到淮南市土壤重金屬Cu含量的空間分布圖(圖4)。由圖4可見,重金屬Cu含量范圍為6.39~170.91 mg/kg,含量高值區為43.10~170.91 mg/kg,主要分布在向陽橋、蒼溝北部、架河、牛家橋及劉魏橋;含量中值區為33.19~43.10 mg/kg,主要分布在界河、田家庵港、幸福堤及窯河紀念碑;含量低值區為6.39~33.19 mg/kg,分布較為分散。

由圖4可知,Cu元素含量高值區較為集中,主要分布在淮南北部的礦區。煤礦開采產生的煤矸石在堆放過程中,受到風化、淋溶等因素的影響,Cu元素從煤矸石中析出產生污染,這與王興明等[4]和孫賢斌等[14]研究的淮南礦區Cu元素污染來源一致。Cu含量中值區主要分布在淮南東部農業集中區附近,農業基地農藥和化肥的投入量很大,直接使用未經處理的畜禽糞便等有機肥料較為普遍。在這種高強度農業活動的影響下,土壤中Cu含量普遍偏高。國內外有研究表明農業活動是土壤Cu元素的主要來源,隨著農田耕種歷史的延長,表層土壤的Cu含量呈增加趨勢。農藥、化肥和有機肥中Cu含量較高,在高強度農業活動的影響下,農田土壤Cu富集,并最終導致農田的重金屬污染[23-25]。由于Cu、Zn、As等微量元素飼料添加劑的普遍使用[26],造成集約化養殖場的畜禽糞便中重金屬含量超標[27,28]。畜禽糞便在農田中的使用,可能成為淮南東部農田土壤重金屬Cu污染的重要來源。

3 小結

淮南市土壤重金屬Cu含量的統計結果表明,土壤Cu的平均含量顯著高于淮南市“十一五”的土壤背景值,具有明顯的累積效應。

Cu含量原始數據經過二階多項式趨勢移除比不經過趨勢移除的插值效果要好。相對普通克里格和泛克里格法,簡單克里格法對極大值的評估更加合理,能夠準確地評價淮南市土壤重金屬Cu元素污染的區域。Cu含量半變異函數模型符合孔穴效應模型,屬中等程度的空間相關,空間分異為隨機性因素和結構性因素共同作用的結果,但主要受人類活動等隨機性因素影響。endprint

淮南市土壤重金屬Cu元素污染高值區分布在淮南北部的礦區,污染中值區分布在淮南東部農業化集中區,主要原因是人類的采礦活動和農田中化肥、畜禽糞便等有機肥料和農藥的過度使用所致。因此,需要加強礦區和農田土壤的預防、監控和治理,對重點污染源和高污染風險區域進行治理和修復。

參考文獻:

[1] 姚尚和.安徽省淮南市礦區土壤污染現狀及風險評估[D].合肥:合肥工業大學,2011.

[2] 崔龍鵬,白建峰,史永紅,等.采礦活動對煤礦區土壤中重金屬污染研究[J].土壤學報,2004,41(6):896-904.

[3] 李海霞,胡振琪,李 寧,等.淮南某廢棄礦區污染場的土壤重金屬污染風險評價[J].煤炭學報,2008,33(4):423-426.

[4] 王興明,董眾兵,劉桂建,等.Zn,Pb,Cd,Cu在淮南新莊孜煤礦矸石山附近土壤和作物中分布特征[J].中國科學技術大學學報,2012,42(1):17-25.

[5] 王曉輝,楊 晨.基于GIS和地統計學的淮南礦區土壤重金屬含量與空間分布研究[J].長江流域資源與環境,2014,23(S):60-65.

[6] GOTWAY C A,FERGUSON R B, HERGERT G W, et al. Comparison of kriging and inverse-distance methods for mapping soil parameters[J]. Soil Sci Soc Am J, 1996, 60(4):1237-1247.

[7] KRAVCHENKO A , BULLOCK D G. A comparative study of interpolation methods for mapping soil properties[J]. Agronomy Journal,1999,91(3):393-400.

[8] BRODSKY L, VANEK V, BAZALOVA M, et al. The differences in the interpolation methods for mapping spatial variability of soil properties[J]. Rostlinna Vyroba,2001,47(12):529-535.

[9] 張朝生,章 申,何建邦.長江水系沉積物重金屬含量空間分布特征研究-地統計學方法[J].地理學報,1997,52(2):184-192.

[10] 曹會聰,王金達,張學林.吉林省農田黑土中Cd、Pb、As含量的空間分布特征[J].環境科學,2006,27(10):2117-2122.

[11] 鐘曉蘭,周生路,李江濤,等.長江三角洲地區土壤重金屬污染的空間變異特征-以江蘇省太倉市為例[J].土壤學報,2007, 44(1):33-39.

[12] 李保杰,顧和和,紀亞洲.基于地統計的礦業城市土壤重金屬污染研究-以徐州市為例[J].江蘇農業科學,2011,39(3):524-526.

[13] 呂建樹,張祖陸,劉 洋,等.日照市土壤重金屬來源解析及環境風險評價[J].地理學報,2012,67(7):971-984.

[14] 孫賢斌,李玉成.淮南大通煤礦廢棄地土壤重金屬空間分布及變異特征[J].地理科學,2013,33(10):1238-1244.

[15] 黃杏元,馬勁松.地理信息系統概論[M].北京:高等教育出版社,2008.101-103.

[16] KALPANA H K, PRAMILA A. Geostatistical analyst for deciding optimal interpolation strategies for delineating compact zones[J]. International Journal of Geosciences, 2011, 2(4):585-596.

[17] 胡克林,李保國,呂貽忠,等.非平穩型區域土壤汞含量的各種估值方法比較[J].環境科學,2004,25(3):132-137.

[18] 史 舟,李 艷.地統計學在土壤學中的應用[M].北京:中國農業出版社,2006.89-91.

[19] 鄭學成.冶煉渣場周邊土壤重金屬污染空間分布研究[D]. 重慶:重慶大學,2011.

[20] 郭廣慧,張航程,彭 穎.基于GIS的宜賓城市土壤Pb含量空間分布特征及污染評價[J].環境科學學報,2011,31(1):164-171.

[21] 夏學齊,陳 駿,廖啟林,等.南京地區表土鎘汞鉛含量的空間統計分析[J].地球化學,2006,35(1):95-102.

[22] CAMBARDELLA C A, MOORMAN T B, NOVAK J M, et al. Field-scale variability of soil properties in central iowa soils[J]. Soil Science Society of America Journal,1994, 58(5):1501-1511.

[23] 張 民,龔子同.我國菜園土壤中某些重金屬元素的含量與分布[J].土壤學報,1996,33(1):85-93.

[24] 陳同斌,黃銘洪,黃煥忠,等.香港土壤中的重金屬含量及其污染現狀[J].地理學報,1997,52(3):228-236.

[25] STANNERS D, BOURDEAU P. Europes environment: the Dobris assessment[M]. Copenhagen: European Environment Agency,1995.146-171.

[26] 張 勇,朱宇旌.飼料與飼料添加劑[M].北京:化學工業出版社,2008.100-105.

[27] NICHOLSON F A, SMITH S R, ALLOWAY B J, et al.An inventory of heavy metals inputs to agricultural soils in England and Wales[J]. Science of the Total Environment, 2003, 311(1-3):205-219.

[28] XIONG X, LI Y X, LIN W, et al. Copper content in animal manures and potential risk of soil copper pollution with animal manure use in agriculture[J]. Resources, Conservation and Recycling, 2010, 54(11):985-990.endprint

猜你喜歡
模型
一半模型
一種去中心化的域名服務本地化模型
適用于BDS-3 PPP的隨機模型
提煉模型 突破難點
函數模型及應用
p150Glued在帕金森病模型中的表達及分布
函數模型及應用
重要模型『一線三等角』
重尾非線性自回歸模型自加權M-估計的漸近分布
3D打印中的模型分割與打包
主站蜘蛛池模板: 亚洲午夜综合网| 伊人色婷婷| 色婷婷亚洲综合五月| 亚洲一区第一页| 狠狠色噜噜狠狠狠狠色综合久| 午夜影院a级片| 精品视频在线观看你懂的一区| 国产理论精品| 久久久久国产一级毛片高清板| h视频在线播放| 国内精品九九久久久精品| 午夜视频免费一区二区在线看| 午夜国产精品视频黄| 国产欧美日韩精品综合在线| 欧美全免费aaaaaa特黄在线| 欧美成人国产| 欧美日韩中文字幕在线| 内射人妻无码色AV天堂| 国产一区二区在线视频观看| 五月丁香伊人啪啪手机免费观看| 亚洲三级电影在线播放| 五月婷婷丁香综合| 免费av一区二区三区在线| 亚洲无码高清一区| 亚洲国产精品无码久久一线| 精品少妇人妻av无码久久| 欧美国产日韩在线| 成人va亚洲va欧美天堂| AV色爱天堂网| 午夜无码一区二区三区在线app| 免费国产无遮挡又黄又爽| 国产区福利小视频在线观看尤物| 欧美国产视频| 国产大片喷水在线在线视频 | 中国成人在线视频| 久久精品无码专区免费| 亚洲综合婷婷激情| 国产91精品最新在线播放| 亚洲 欧美 中文 AⅤ在线视频| 看国产毛片| 天天躁夜夜躁狠狠躁躁88| 精品视频福利| 国产JIZzJIzz视频全部免费| 一级毛片在线免费视频| 欧美精品在线免费| 欧美精品影院| 国产91无码福利在线| 精品一区二区三区无码视频无码| 久久久噜噜噜久久中文字幕色伊伊 | 欧美日韩中文字幕二区三区| 91久久夜色精品国产网站| 一本一本大道香蕉久在线播放| 美女毛片在线| 亚欧成人无码AV在线播放| 91 九色视频丝袜| 亚洲大尺度在线| 亚洲欧美自拍中文| 色婷婷在线播放| 国内精品手机在线观看视频| 欧美亚洲国产精品久久蜜芽| 亚洲熟女偷拍| 少妇精品久久久一区二区三区| 91免费精品国偷自产在线在线| 国内精品自在自线视频香蕉| 精品少妇人妻一区二区| 91精品伊人久久大香线蕉| 日韩精品久久无码中文字幕色欲| 国产欧美网站| 国产网站免费看| 亚洲 欧美 中文 AⅤ在线视频| 国产精品女熟高潮视频| 国产无人区一区二区三区| 国产激情影院| 人妻丝袜无码视频| 国产91无毒不卡在线观看| 国产午夜小视频| 性色一区| 精品国产aⅴ一区二区三区| 好吊色国产欧美日韩免费观看| 亚洲男人的天堂在线| 亚洲二区视频| 免费一级成人毛片|