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

基于地理加權(quán)回歸模型探究環(huán)境異質(zhì)性對(duì)秦嶺大熊貓空間利用的影響

2020-06-11 13:40:52薛瑞暉于曉平李東群葉新平
生態(tài)學(xué)報(bào) 2020年8期
關(guān)鍵詞:景觀影響模型

薛瑞暉,于曉平,李東群,葉新平,*

1 陜西師范大學(xué)生命科學(xué)學(xué)院, 西安 710119 2 陜西師范大學(xué)易科泰無人機(jī)遙感生態(tài)研究中心, 西安 710119 3 陜西周至老縣城國家級(jí)自然保護(hù)區(qū)管理局, 西安 710400

野生動(dòng)物的分布受到其棲息區(qū)域景觀空間格局的顯著影響[1]。景觀空間格局既是景觀異質(zhì)性本身的具體體現(xiàn),又是野生動(dòng)物棲息地生境特征的有效反映[1]。因此,了解景觀空間格局特征對(duì)野生動(dòng)物空間利用的影響,對(duì)物種的科學(xué)保護(hù)和有效管理有著重要的意義[2],已有的景觀格局與野生動(dòng)物分布關(guān)系的研究多借助于景觀指數(shù)對(duì)景觀格局進(jìn)行量化,然后采用傳統(tǒng)的線性回歸模型等方法定量分析景觀格局對(duì)野生動(dòng)物分布的影響[3- 9]。然而,生態(tài)學(xué)關(guān)系的空間異質(zhì)性或空間非均勻性(即變量間的關(guān)系或結(jié)構(gòu)會(huì)隨著地理位置的變化而改變),是空間分析中一個(gè)普遍存在的現(xiàn)象[10]。傳統(tǒng)線性回歸模型,如普通最小二乘法(Ordinary Least Square,OLS)模型,將自變量與因變量的關(guān)系視為全局不變,得到的回歸參數(shù)是整個(gè)研究區(qū)域內(nèi)的平均值,從而忽略了兩者的空間變異性所導(dǎo)致的關(guān)系變化,并因此造成回歸模型有效性降低和預(yù)測結(jié)果失效等問題[10,11]。

地理加權(quán)回歸(Geographically weighted regression,GWR)模型是傳統(tǒng)回歸分析方法的擴(kuò)展,是探索空間關(guān)系異質(zhì)性方面有效工具之一[10]。GWR模型是一種典型的局部空間回歸模型,通過將全局參數(shù)分解成局部參數(shù)進(jìn)行估計(jì),對(duì)每一個(gè)空間位置點(diǎn)的參數(shù)進(jìn)行估計(jì)時(shí)考慮了非平穩(wěn)性的關(guān)系,因此能深刻解釋地理空間數(shù)據(jù)的某類指標(biāo)和影響因子之間的空間非平穩(wěn)性關(guān)系,這是傳統(tǒng)OLS模型無法比擬的[12-13]。近年來,許多學(xué)者對(duì)GWR模型進(jìn)行了研究,主要集中在應(yīng)用GWR模型分析社會(huì)環(huán)境因子對(duì)區(qū)域經(jīng)濟(jì)、區(qū)域房價(jià)和區(qū)域污染等方面的研究[14-18],如關(guān)偉和郝金連[14]運(yùn)用OLS和GWR模型分析東北地區(qū)旅游經(jīng)濟(jì)和旅游產(chǎn)業(yè)因子、消費(fèi)因子、投資因子之間的關(guān)系,發(fā)現(xiàn)GWR模型擬合優(yōu)度比OLS模型有顯著提高。陳輝等[17]采用GWR模型構(gòu)建了我國區(qū)域范圍內(nèi)近地表的PM2.5遙感反演模型,結(jié)果表明GWR模型既能體現(xiàn)PM2.5時(shí)空分布的全局變化特征,又能體現(xiàn)局部空間異質(zhì)性。陳強(qiáng)等[18]基于GWR模型評(píng)估了土地利用對(duì)地表水質(zhì)的影響,其結(jié)果也證明GWR模型在預(yù)測精度和處理空間自相關(guān)過程中都要優(yōu)于OLS模型。然而,在野生動(dòng)物棲息地選擇與利用研究領(lǐng)域,GWR方法尚未得到應(yīng)用。

大熊貓(Ailuropodamelanoleuca)是中國乃至世界野生動(dòng)物保護(hù)的旗艦物種,是中國特有珍稀瀕危物種[19,20]。近年來應(yīng)用景觀生態(tài)學(xué)原理與方法對(duì)大熊貓的研究很多,內(nèi)容涉及生境適宜性評(píng)價(jià)、潛在棲息地分布模擬、景觀連通性分析及棲息地景觀破碎化等[21-26]。這些研究主要是利用景觀格局指數(shù)描述大熊貓棲息地特征,或者使用傳統(tǒng)線性回歸方法分析景觀格局對(duì)大熊貓的影響,很少有研究涉及從整體到局部探討環(huán)境因子對(duì)大熊貓空間分布的影響。本文結(jié)合前人的研究結(jié)果,采用GWR方法探討環(huán)境異質(zhì)性對(duì)秦嶺大熊貓空間分布的影響,并通過與傳統(tǒng)的OLS分析方法進(jìn)行比較,檢驗(yàn)GWR模型在大熊貓棲息地空間分析中是否可用,以期提高大熊貓-環(huán)境空間關(guān)系的科學(xué)理解,為進(jìn)一步開展大熊貓等野生動(dòng)物的棲息地選擇與利用精準(zhǔn)量化分析和評(píng)價(jià)提供參考。

1 材料與方法

1.1 研究區(qū)域概況

本研究區(qū)域介于北緯33°7′—34°6′、東經(jīng)106°21′—108°54′和海拔376—3770 m之間,屬于我國傳統(tǒng)南北地理分界線的秦嶺中段地區(qū),轄陜西省17個(gè)縣部分區(qū)域,且包含19個(gè)自然保護(hù)區(qū),面積約為2萬km2(圖1)。該研究區(qū)屬于季風(fēng)氣候,是北亞熱帶和暖溫帶的過渡區(qū),是長江和黃河的分水嶺,同時(shí)包含秦嶺山地、低山丘陵和平原等一系列地貌類型,是我國生物多樣性和生態(tài)系統(tǒng)類型最為豐富的地區(qū)之一。秦嶺大熊貓主要分布在秦嶺中段南坡,北坡和西段有少量分布[27]。總體上為兩塊相隔較遠(yuǎn)的區(qū)域,其西端的寧強(qiáng)縣和甘肅康縣、四川青川連成一片;佛坪、太白、洋縣、周至、寧陜、留壩、城固、戶縣和鎮(zhèn)安9縣的分布區(qū)相隔較近或連成一片,分布區(qū)域橫跨長江、黃河水系[27]。

1.2 數(shù)據(jù)來源及處理

本研究所采用的秦嶺大熊貓分布數(shù)據(jù)來源于1999—2002年全國第三次大熊貓調(diào)查數(shù)據(jù)。根據(jù)秦嶺大熊貓平均家域面積為10.62 km2[28],為避免空間自相關(guān)影響模型預(yù)測,剔除距離小于2 km的冗余痕跡點(diǎn),最終保留267個(gè)大熊貓分布痕跡點(diǎn)進(jìn)入分析(圖1紅點(diǎn)所示)。鑒于OLS與GWR建模要求使用連續(xù)型數(shù)值變量,我們利用ArcGIS 10.2軟件的核密度(Kernel Density)分析功能估計(jì)了每一個(gè)大熊貓分布點(diǎn)所在位置的大熊貓分布相對(duì)密度。

研究區(qū)的數(shù)字高程模型(DEM)數(shù)據(jù)源于美國SRTM(Shuttle Radar Topography Mission,SRTM)30 m精度數(shù)據(jù);公路和居民點(diǎn)GIS數(shù)據(jù)源于中國科學(xué)院資源環(huán)境科學(xué)數(shù)據(jù)中心的國家基礎(chǔ)矢量數(shù)據(jù)集數(shù)據(jù);土地覆蓋/利用數(shù)據(jù)源于中國土地利用現(xiàn)狀遙感監(jiān)測數(shù)據(jù)庫2000年土地利用30 m精度數(shù)據(jù)(http://www.resdc.cn),該數(shù)據(jù)集的土地利用類型包括耕地、林地、草地、水域、居民地和未利用土地等6個(gè)一級(jí)類型和包括有林地、灌木林、疏林地等25個(gè)二級(jí)類型,數(shù)據(jù)分類精度均可達(dá)到85%以上[29]。參考以往大熊貓生態(tài)學(xué)方面的研究[21- 26],本研究選取了有林地、農(nóng)田、灌草地等三類與大熊貓分布相關(guān)的土地覆蓋類型進(jìn)入下一步景觀空間格局量化分析。

1.3 環(huán)境變量篩選

本研究選取了以下3個(gè)方面的環(huán)境因素來分析大熊貓空間分布與環(huán)境異質(zhì)性之間的關(guān)系:1)地理特征(包括海拔和坡度);2)人類影響(包括距離道路的距離和距離居民點(diǎn)的距離);3)景觀格局特征(包括最大斑塊指數(shù)、邊緣密度、斑塊面積、連接度指數(shù)和聚合度指數(shù))[30]。每一個(gè)大熊貓分布點(diǎn)處的人類影響因素通過ArcGIS 10.2軟件中的距離量算功能計(jì)算距離道路和居民點(diǎn)的最近距離。大熊貓分布點(diǎn)處的景觀指數(shù)值是在Fragstats 4.0軟件中選擇User provided points模塊,設(shè)置搜索半徑為3 km計(jì)算完成的[30]。為避免景觀指數(shù)值的空間自相關(guān),本文利用R 3.3.2軟件對(duì)所提取景觀格局指數(shù)值進(jìn)行因子分析,通過變量間的相關(guān)關(guān)系剔除信息嚴(yán)重冗余變量,將特征值>0.8的變量保留,特征值<0.3的變量移除[30],最終選擇農(nóng)田斑塊面積(Cropland-CA)、有林地斑塊面積(Forest-CA)和灌草地斑塊面積(Grassland-CA)等獨(dú)立的景觀變量作為模型變量參與下一步的模型運(yùn)算。最終進(jìn)入建模分析的環(huán)境變量見表1。

1.4 模型建立與評(píng)價(jià)

本研究首先選用經(jīng)典的OLS 進(jìn)行建模。通過OLS 建模,既可以初步了解變量間的相互關(guān)系,也可以評(píng)估OLS 對(duì)空間數(shù)據(jù)的建模效果,進(jìn)而確定是否選擇GWR進(jìn)一步建模。本文所有的模型分析在GWR 4.0軟件中完成。

表1 本研究模型分析中所用的環(huán)境因素及變量

1)OLS建模。OLS是一種對(duì)因變量和解釋變量之間相互關(guān)系進(jìn)行研究的統(tǒng)計(jì)方法,是所有空間回歸分析的起點(diǎn),假設(shè)線性回歸關(guān)系滿足全局空間平穩(wěn)條件,其公式如下:

Yi=β0+∑kβkXik+εi

(1)

式中,Yi是第i點(diǎn)因變量的值,β0為截距,Xik為第k個(gè)解釋變量在第i點(diǎn)的值,βk為第k個(gè)解釋變量的斜率或回歸系數(shù),εi為殘差。

2)GWR建模。GWR是Fotheringham等在傳統(tǒng)OLS模型基礎(chǔ)上將數(shù)據(jù)的地理位置加入回歸參數(shù)中,同時(shí)考慮了相鄰點(diǎn)的空間權(quán)重,允許局部參數(shù)估計(jì)的地學(xué)統(tǒng)計(jì)方法[31,32]。其公式如下:

Yi=β0(μi,νi)+∑kβk(μi,νi)Xik+εi

(2)

式中,(μi,νi)是第i點(diǎn)的空間位置,βk(μi,νi)是連續(xù)函數(shù)βk(μ,ν)在(μi,νi)處的值。

由于大熊貓的分布在空間上呈集聚分布,故GWR模型采用調(diào)整型的空間核和Gaussian權(quán)重函數(shù),以AICc 作為衡量標(biāo)準(zhǔn),通過黃金分割搜索(Golden section search)方法來確定最佳帶寬,也就是說,使得AICc最小的帶寬即為最佳帶寬。

3)模型評(píng)價(jià)。為了比較OLS模型和GWR模型的擬合效果,選擇修正的AIC信息準(zhǔn)則(Akaike,1973)作為模型評(píng)價(jià)指標(biāo),其公式如下:

(3)

式中,n表示觀測值數(shù)量;σ表示誤差項(xiàng)估計(jì)的標(biāo)準(zhǔn)差;tr(S)表示GWR模型S矩陣的跡,是帶寬的函數(shù)。一般情況下,具有較小AICC值的模型擬合效果相對(duì)較好。若兩個(gè)模型間的AICC差>3,則表明兩個(gè)模型之間具有明顯差異,具有較小AICC值的模型擬合度更優(yōu)。此外,還可以通過計(jì)算模型的殘差標(biāo)準(zhǔn)差(Sigma)來比較OLS模型和GWR模型,Sigma值越小,表示模型的擬合效果越好。

2 結(jié)果與分析

2.1 最小二乘法(OLS)模型分析

以大熊貓相對(duì)分布密度為因變量建立OLS 回歸方程結(jié)果表明(表2、圖2),大熊貓分布相對(duì)密度與森林面積占比和灌草地面積占比呈顯著相關(guān)關(guān)系,與海拔和距離道路距離也呈現(xiàn)出一定的相關(guān)性。各解釋變量的方差膨脹因子(VIF)均小于7.5,表明各變量間共線性對(duì)模型的影響不顯著。OLS模型的R2為0.152,表明模型的擬合度不理想。OLS模型殘差的全局Moran′sI系數(shù)為0.16,Z-score 為10.44,說明殘差項(xiàng)還存在空間自相關(guān)性。以上說明,OLS對(duì)數(shù)據(jù)空間信息的提取還不完全,不能很好地解釋大熊貓空間分布密度與環(huán)境影響因素的空間非均勻性。

表2 OLS模型參數(shù)估計(jì)及檢驗(yàn)結(jié)果

***、**、*分別表示通過1%、5%、10%水平下的顯著性檢驗(yàn)

2.2 地理加權(quán)回歸(GWR)模型分析

使用與OLS模型相同的變量構(gòu)建GWR模型,依據(jù)最小AIC值確定最佳帶寬為107獲得各變量系數(shù)分布(表3、圖3),可以看到,森林面積占比和灌草地面積占比等變量的參數(shù)估計(jì)值在OLS方法下是正值,但在GWR方法下則有正有負(fù),說明各環(huán)境變量對(duì)大熊貓空間分布的影響呈現(xiàn)出空間不均勻性。

圖2 OLS模型影響因子回歸系數(shù)Fig.2 OLS model influence factor regression coefficient

圖3 GWR模型影響因子回歸系數(shù)Fig.3 GWR model influence factor regression coefficient

表3 GWR模型參數(shù)估計(jì)及檢驗(yàn)結(jié)果

GWR模型與OLS模型相比(表4),GWR模型的R2達(dá)到了0.468,且與OLS模型的AICc之差遠(yuǎn)大于3(ΔAICc=24.401),說明GWR模型的擬合度明顯高于OLS模型;GWR模型殘差的方差分析結(jié)果顯示,GWR模型的殘差平方和(Residuals SS)下降了84.088,亦表明GWR 模型的擬合效果比OLS 模型有很大改善(ResidualsF= 2.368);GWR模型殘差的Moran′sI系數(shù)為0.09,Z-score 為6.08,說明殘差不存在空間自相關(guān)。綜上,結(jié)果表明GWR模型比OLS模型模擬的效果要好。

表4 GWR模型和OLS模型的比較

2.3 影響因子空間分異特征舉例分析

GWR模型是局部模型,在每個(gè)樣本數(shù)據(jù)點(diǎn)都有一組局部的參數(shù)估計(jì),而不像OLS 模型那樣是一個(gè)全局或平均意義上的估計(jì)值。因此,我們對(duì)GWR模型的局部參數(shù)估計(jì)值使用IDW方法進(jìn)行空間插值,得到各環(huán)境變量的局部參數(shù)估計(jì)的空間格局圖(圖4),以此展示各環(huán)境變量對(duì)大熊貓分布影響的空間分異特征。

如圖4所示,距居民點(diǎn)距離對(duì)大熊貓相對(duì)分布密度的影響作用自西向東呈遞增趨勢,形成以鳳縣、勉縣、留壩縣、城固縣、太白縣和洋縣西北部為主的低值區(qū)和以周至縣和佛坪縣為主的高值區(qū);距道路距離對(duì)大熊貓相對(duì)分布密度的影響作用從中心向四周呈圈層遞增趨勢,形成以太白縣東南部、周至縣西南部、佛坪縣西北部和洋縣東北部為主的低值區(qū)和以鳳縣、城固縣、留壩縣、勉縣和略陽縣為主的高值區(qū)格局;而有林地面積占比對(duì)大熊貓分布相對(duì)密度的影響作用從中心向四周呈圈層遞增趨勢,形成以周至縣南部和佛坪縣北部為主的低值區(qū)和以城固縣為主的高值區(qū)格局,表明有林地斑塊面積占比對(duì)周至縣南部和佛坪縣北部的大熊貓核心分布區(qū)域的相對(duì)分布密度的影響作用不明顯。

圖4 第三次大熊貓調(diào)查GWR模型影響因子回歸系數(shù)空間分布Fig.4 The third giant panda survey GWR model influence factor regression coefficient spatial distribution A:距居民點(diǎn)距離對(duì)大熊貓相對(duì)分布密度的影響;B:距道路距離對(duì)大熊貓相對(duì)分布密度的影響;C:有林地斑塊面積對(duì)大熊貓相對(duì)分布密度的影響

3 討論

本研究分別應(yīng)用OLS模型和GWR模型分析了環(huán)境異質(zhì)性與大熊貓空間分布之間的關(guān)系,并將兩者的模擬結(jié)果進(jìn)行比較,證明GWR模型對(duì)于刻畫大熊貓空間分布與景觀格局之間的關(guān)系具有獨(dú)特的空間非平穩(wěn)性和尺度依存特性優(yōu)勢。OLS模型中多采用全局空間回歸模型量化大熊貓分布與環(huán)境變量之間的關(guān)系,得到的回歸參數(shù)估計(jì)是在整個(gè)研究區(qū)域內(nèi)的平均值,不能反映物種-環(huán)境關(guān)系的空間異質(zhì)特征;而GWR模型允許將全局參數(shù)的估計(jì)分解成局部參數(shù)進(jìn)行估計(jì),在模擬精度和空間特征提取上都顯著優(yōu)于OLS模型,模型的結(jié)果可反映環(huán)境變量的影響具有空間非均勻性,這是傳統(tǒng)OLS模型無法實(shí)現(xiàn)的。

參考已有大熊貓與環(huán)境異質(zhì)性之間關(guān)系的研究[21- 26],本文的GWR分析結(jié)果進(jìn)一步說明環(huán)境異質(zhì)性對(duì)大熊貓空間利用的影響呈現(xiàn)出明顯的空間非平衡性,不同區(qū)域大熊貓空間分布與主導(dǎo)因子發(fā)生作用的特征尺度也不同。例如,有林地斑塊面積對(duì)大熊貓相對(duì)分布密度的影響作用從中心向四周呈圈層遞增趨勢,形成以周至縣南部和佛坪縣北部為主的低值區(qū)和以城固縣為主的高值區(qū)格局,表明有林地斑塊面積對(duì)周至縣南部和佛坪縣北部的大熊貓核心分布區(qū)域的相對(duì)分布密度的影響作用不明顯,而對(duì)城固縣大熊貓相對(duì)分布密度的影響作用要大于其他區(qū)域。距離道路距離對(duì)大熊貓相對(duì)分布密度的影響作用從中心向四周呈圈層遞增趨勢,形成以太白縣東南部、周至縣西南部、佛坪縣西北部和洋縣東北部為主的低值區(qū)和以鳳縣、城固縣、留壩縣、勉縣和略陽縣為主的高值區(qū)格局,說明道路對(duì)鳳縣、城固縣、留壩縣、勉縣和略陽縣大熊貓相對(duì)分布密度的影響作用要大于其他區(qū)域,這就要求我們保護(hù)區(qū)管理規(guī)劃者要格外注意對(duì)不同區(qū)域采用整體和局部相結(jié)合的辦法,做到因地制宜,提升保護(hù)措施有效性。

雖然GWR模型在理解野生動(dòng)物空間分布與環(huán)境異質(zhì)性之間的空間關(guān)聯(lián)具有很大優(yōu)勢,但模型本身也存在一些不足之處。首先,由于回歸模型對(duì)環(huán)境變量數(shù)依賴性較大,當(dāng)環(huán)境變量數(shù)較少時(shí),難以建立合理的回歸模型,對(duì)因變量的預(yù)測不準(zhǔn)確;而當(dāng)環(huán)境變量數(shù)較多時(shí),又會(huì)產(chǎn)生多重共線性,同樣使得結(jié)果出現(xiàn)不確定性。其次,GWR模型不支持類型變量的加入,如土地利用類型和保護(hù)地類型等,而野生動(dòng)物空間分布卻受這些類型要素的影響較大,忽視重要的類型變量會(huì)降低模型預(yù)測精度。如何將類別變量合理的應(yīng)用到GWR模型中還需要進(jìn)一步研究。此外,本研究僅分析了大熊貓空間分布與部分環(huán)境因子之間的回歸關(guān)系,實(shí)際上,除了這些環(huán)境因子之外諸多環(huán)境因子也影響著大熊貓的空間分布,如氣候、水域和食物等。下一步工作將把這些環(huán)境因子融合到GWR模型中,以期得到更為全面客觀的分析結(jié)果。

4 結(jié)論

本文基于地理加權(quán)回歸模型(GWR)和普通最小二乘法線性回歸模型(OLS)探究了環(huán)境異質(zhì)性對(duì)秦嶺大熊貓空間利用的影響,并對(duì)模型模擬結(jié)果進(jìn)行了比較。主要結(jié)論如下:

(1)GWR模型可解釋環(huán)境異質(zhì)性對(duì)大熊貓空間利用影響的46.8%,優(yōu)于OLS模型的15.2%。

(2)GWR模型適用于野生動(dòng)物分布與環(huán)境變量之間關(guān)系的建模分析,能有效的刻畫野生動(dòng)物分布與環(huán)境變量之間的空間異質(zhì)性,可為探究物種-環(huán)境關(guān)系的空間異質(zhì)特征提供一種新的思路和方法。

(3)大熊貓空間分布具有明顯的空間異質(zhì)性特征,不同環(huán)境因子對(duì)大熊貓空間分布影響的作用不同,現(xiàn)有大熊貓空間分布格局受多重因素綜合影響。

致謝:感謝陜西省林業(yè)局、陜西周至老縣城國家自然保護(hù)區(qū)以及陜西佛坪國家自然保護(hù)區(qū)給予的技術(shù)支持。

猜你喜歡
景觀影響模型
一半模型
是什么影響了滑動(dòng)摩擦力的大小
景觀別墅
哪些顧慮影響擔(dān)當(dāng)?
火山塑造景觀
重要模型『一線三等角』
包羅萬象的室內(nèi)景觀
重尾非線性自回歸模型自加權(quán)M-估計(jì)的漸近分布
3D打印中的模型分割與打包
擴(kuò)鏈劑聯(lián)用對(duì)PETG擴(kuò)鏈反應(yīng)與流變性能的影響
中國塑料(2016年3期)2016-06-15 20:30:00
主站蜘蛛池模板: 人妻丰满熟妇αv无码| 激情無極限的亚洲一区免费| 久久99久久无码毛片一区二区| 99久视频| 国产丝袜无码精品| 青青草原偷拍视频| 亚洲中文字幕国产av| 香蕉视频国产精品人| 精久久久久无码区中文字幕| 国产精品一老牛影视频| 天天色天天操综合网| 色婷婷综合在线| 人妻中文久热无码丝袜| 精品无码视频在线观看| 在线国产资源| 国产95在线 | 国产精品亚洲天堂| 一区二区三区在线不卡免费| 国产在线拍偷自揄观看视频网站| 不卡网亚洲无码| 99re这里只有国产中文精品国产精品 | 欧美精品三级在线| 青草免费在线观看| 欧美激情视频在线观看一区| 香蕉久人久人青草青草| 91精品啪在线观看国产60岁| 91毛片网| 国产主播福利在线观看| 成人在线不卡| 久久这里只精品国产99热8| 国产欧美成人不卡视频| 人妻精品全国免费视频| 国产精品三级专区| 国产原创第一页在线观看| 国产精品三级专区| 国产精品国产三级国产专业不| 91无码人妻精品一区二区蜜桃| 国产9191精品免费观看| 无码精品一区二区久久久| 亚洲欧洲天堂色AV| 亚洲动漫h| 狠狠干欧美| 国产白浆在线观看| 国产精品9| 国产精品人人做人人爽人人添| 国产女人综合久久精品视| 成人午夜网址| 日韩麻豆小视频| 67194成是人免费无码| 成人在线亚洲| 波多野结衣视频网站| 国产高清免费午夜在线视频| 国产成人精品免费视频大全五级| 人妻夜夜爽天天爽| 亚洲三级电影在线播放| 国产精品大白天新婚身材| 99久久精品国产综合婷婷| 日本a∨在线观看| 国产女人喷水视频| 中字无码av在线电影| 国产精品成人观看视频国产| 亚洲黄色成人| 2021亚洲精品不卡a| 免费在线观看av| 国产精品三级专区| 日韩欧美中文字幕一本| 欧美日韩精品在线播放| 国产一二视频| 香蕉综合在线视频91| 成人av手机在线观看| 免费在线不卡视频| 狠狠躁天天躁夜夜躁婷婷| 亚洲欧美日韩天堂| 国产精品lululu在线观看| 色婷婷亚洲综合五月| 欧美一级高清免费a| 国产黑丝视频在线观看| 性色一区| 欧美日本中文| 在线va视频| 欧美人在线一区二区三区| 98超碰在线观看|