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

基于GIS和MaxEnt模型分析氣候變化背景下紫果云杉的潛在分布區

2022-04-18 08:23:42曹家豪高本強
西北植物學報 2022年3期
關鍵詞:物種模型

劉 婷,曹家豪,3*,齊 瑞,李 波,高本強

(1 甘肅省白龍江林業科學研究所,蘭州730070;2 甘肅白龍江森林生態系統國家定位觀測研究站,甘肅舟曲 746300;3 蘭州大學 生命科學學院,蘭州730070)

紫果云杉是中國特有樹種,材質優良,生長快,是長江、黃河上游水源涵養區重要的建群種和森林群落演替的頂級種,在水源涵養、防止水土流失、保護生物多樣性等方面具有重要的作用。紫果云杉的海拔分布范圍廣,耐陰性很強,是生態環境脆弱、氣候條件惡劣的高海拔地區造林綠化的首選樹種,在維護流域生態安全和保障流域生態高質量發展等方面具有十分重要的作用。因此,開展在各項環境因子的制約條件下紫果云杉分布范圍的預測模擬工作并進行適宜性評價,了解制約其生長的主導環境因子,定量描述影響其生長、生存的主要環境因子,能夠對開展紫果云杉資源調查、掌握其分布狀況、實施資源保護提供有價值的參考,為在適生區開展人工栽培及播種育苗工作、有效推進造林綠化提供一定的理論依據,對長江、黃河上游水源涵養區生態環境保護工作具有十分重要的意義。

物種分布模型(species distribution models, SDMs)通常利用存在的發生數據和假定影響其分布的環境變量預測物種潛在分布的地理范圍[1-2],是生態學領域研究物種與環境變量之間的關系及物種適生性的重要工具[3],目前已有多種物種分布模型在模擬物種分布區得到了廣泛應用,其中最大熵(the maximum entropy, MaxEnt)通過已知樣本點和對應環境變量具有最大熵值時擬合物種的分布區,在環境變量的影響下預測物種存在的概率[4],不受樣本量的影響,能夠生成物種響應曲線,定量分析適宜生境的環境因子,是物種分布模型中應用最廣泛[5]、預測精度較高[6]的模型。目前MaxEnt在梭梭[7]、獨葉草[8]、冬青樹[9]、青海云杉[10]、白櫟[11]等植物的潛在分布區預測及適宜性評價等眾多相關研究所采用。

目前,針對紫果云杉的研究主要集中在種群結構特征[12]及物種組成[13]、種子萌發[14]、育苗技術[15]、徑向生長對氣候的響應[16]等方面,另外,王婧如等利用物種分布模型的研究側重于紫果云杉與其親本種生態位的環境差異[17],李寧寧等側重于在氣候變化下紫果云杉與麗江云杉分布區的變化情況[18],二者均使用第五次國際耦合模式比較計劃(CMIP5)不同大氣環流模型中代表性濃度路徑情景(representative concentration pathways, RCPs)模擬未來紫果云杉的潛在分布,而針對主導紫果云杉分布的環境因子及潛在分布區變化的定量分析卻鮮見報道。另外,第六次國際耦合模式比較計劃(CMIP6)對未來的氣候模擬數據已更新,逐漸應用于模型預測[19],且與CMIP5提供的代表濃度路徑(RCPs)只考慮CO2濃度與氣候關系的氣候情景模式不同,其使用共享社會經濟路徑(SSPs)預測未來不同氣候政策下溫室氣體的排放情景,綜合考慮了社會經濟、土地利用對區域氣候變化發展的影響,預測結果與真實值更接近[20]。本研究擬利用數字化標本信息及野外調查采集分布點數據,結合最大熵模型與GIS技術,利用22個涵蓋溫度、降水、地形3方面的環境因子,使用CMIP6提供的未來2個時期3種氣候情景模式數據,分析不同時期紫果云杉的分布區變化,繪制影響紫果云杉生長的主要環境因子的響應曲線,定量描述適宜紫果云杉生長、生存的環境條件,以期更好地開展紫果云杉林資源調查、野生撫育及種植區劃工作,為紫果云杉的保護與管理提供理論支撐。

1 材料和方法

1.1 紫果云杉分布數據及處理方法

本文使用本課題組在甘肅省甘南藏族自治州卓尼縣、夏河縣、臨潭縣、碌曲縣、迭部縣及四川省阿壩藏族羌族自治州若爾蓋縣的30個紫果云杉固定監測樣地的分布數據,以及中國數字植物標本館(Chinese Virtual Herbarium,CVH)、全球生物多樣性信息庫(Global Biodiversity Information Facility,GBIF)的記載數據。將中國數字植物標本館檢索得到的紫果云杉分布數據去除重復樣本、地理位置不明確及人工引種栽培的記錄,通過BIGEMAP地圖下載器獲取分布點的經緯度坐標信息,共獲得83處分布數據;從GBIF檢索到348條紫果云杉分布數據,剔除無明確坐標信息的數據后共獲得164處紫果云杉分布數據,加上野外調查取樣的30處分布數據,共獲得277條紫果云杉有效分布數據。

將獲得的紫果云杉分布數據利用ArcGIS 10.2軟件進行緩沖區分析,對277處分布數據進行篩選,刪除重復、空間關聯性較大的分布數據,以免造成過擬合模擬的影響。由于選取的環境變量柵格數據空間分辨率為2.5 arc-minutes(~5 km),故設置2.5 km的緩沖半徑,在距離<5 km的分布點中只取唯一分布點,最終得到有效紫果云杉分布樣點108個,保存為.csv格式待用。

1.2 環境數據及處理方法

氣候因子對植物的地理分布有非常重要的影響,本研究使用的當前(1970-2000年)及未來2個時期(2041-2060年、2081-2100年)的氣候數據均下載自WorldClim全球氣候數據庫(https://www. worldclim.org),坐標系為WGS84,空間分辨率為2.5 arc-minutes(~5 km),包含19個氣候變量(bio1~bio19,表1所示)。未來2個時期的氣候數據選擇第6次國際耦合模式比較計劃(CMIP6)BCC-CSM2-MR共享社會經濟路徑(SSPs)中代表了可持續發展(ssp126)、局部發展(ssp370)、常規發展(ssp585)3種情景模式的數據。將3個時間段的氣候數據進行處理,提取出19個氣候變量柵格數據裁剪為研究區范圍并轉換為.asc格式待用。

海拔與氣溫﹑氧氣含量、氣壓等關系密切,坡度、坡向體現了太陽光照條件[21],對植物生長影響顯著,本文使用的海拔(altitude)、坡度(slope)、坡向(aspect)3種地形變量是通過地理空間數據云STRM30m高程數據結合GIS技術提取而來。考慮到選取的環境數據的精度為2.5 arc-minutes,故將3種地形變量進行重采樣,轉換成與環境變量一致的空間分辨率,裁剪為研究區范圍并轉換為.asc格式待用。

為減少22個環境因子間的自相關性干擾,本研究通過Pearson相關系數(r)檢驗各變量之間的多重共線性,一組高度交叉相關的變量(|r|>0.8)將保留貢獻率較大的變量納入模型分析中,最后選取了8個環境變量用于最大熵模型分析,如表1所示,這些變量分別是年平均降水量(bio12)、最干月降水量(bio14)、降水量變異系數(bio15)、最冷季度平均溫度(bio11)、年均溫變化范圍(bio7)、平均日溫差(bio2)、等溫性(bio3)、海拔(alt)。

表1 最大熵模型中各個環境變量的貢獻率分布數據

1.3 MaxEnt模型預測分析

本研究使用最大熵模型MaxEnt 3.4.4模擬預測紫果云杉的分布區。將篩選出的108處紫果云杉分布數據與8個環境氣候變量導入MaxEnt軟件中,設置25%的分布點數據作為測試數據,75%為訓練數據,重復運行10次,使用受試者工作特征曲線(receiver operating characteristic curve, ROC曲線)評價模型預測結果的精度,使用刀切法Jackknife確定環境變量的貢獻。

2 結果與分析

2.1 預測精度

受試者工作特征曲線下面積(area under receiver operating characteristic curve, AUC)表示模型模擬精度[11,22],通常AUC值介于0~1.0之間[23],AUC值為0.50表示模型的預測結果并不優于隨機結果,當AUC值≥0.9時預測結果精度較高,且越接近1,表示模型預測的結果越準確[24],可以較為準確地反映物種分布區的范圍。

本次研究MaxEnt模型重復運行10次后訓練數據AUC值的平均值為0.972(標準差±0.009),測試數據平均值為0.967(標準差±0.009),如圖1所示。另外,用于模型預測模擬的108處紫果云杉分布數據全部位于分布區內,且70%的分布點數據分布于最適生區,21%位于較適生區,9%位于低適生區,模擬結果與紫果云杉實際分布數據具有較好的一致性。本次模型對紫果云杉分布區的預測效果非常好,具有很高的可信度。

圖1 基于Maxent模型預測紫果云杉分布的受試者工作特征曲線下面積

2.2 當前氣候下紫果云杉在中國的分布

如圖2所示,紫果云杉分布區主要集中在以祁連山脈-橫斷山脈以西,青藏高原東緣地區的四川西北部、甘肅南部、青海東南部、西藏東部等地,其他地區偶有零星分布。

圖2 紫果云杉的地理分布情況及當前氣候環境下(1970-2000)的潛在分布區

利用地理信息系統空間分析模塊提取各省份紫果云杉最適生區、較適生區、低適生區面積(表2所示)。當前氣候環境下,中國適宜紫果云杉生長的面積為77.553 8 萬km2,僅占中國陸地面積的約8%。最適生區面積為4.093 7 萬km2,主要分布在四川西北部、甘肅南部;較適生區面積16.734 4 萬km2,四川、甘肅、青海、西藏地區分布較廣;低適生區面積56.725 7 萬km2,分布范圍最廣,占中國適生區面積的73%以上。

表2 當前氣候條件下(1970-2000)紫果云杉適生區面積統計

2.3 影響紫果云杉分布的主要環境變量

通過MaxEnt模型Jackknife模塊分析當前氣候環境下不同環境變量對紫果云杉分布預測的影響權重(圖3),Without Variable表示缺失該變量時,模型的正規化訓練增益值,With Only Variable表示僅有該變量時模型的正規化訓練增益值,With All Variable 表示全部變量參與運算時,模型的正規化訓練增益值,當僅有單變量參與模型時,海拔(alt)和年平均降水量(bio12)得分值超過1,當缺失單個變量時,海拔(alt)和年平均降水量(bio12)收益降低最大,表明這二者是最重要的變量。結合表1中各環境變量的貢獻率,貢獻率較高的環境變量為海拔alt(44.0%)、年平均降水量bio12(32.1%)、最干月降水量bio14(8.7%)、最冷季度平均溫度bio11(5%)、降水量變異系數bio15(4.8%),而年均溫變化范圍bio7、等溫性bio3和平均日溫差bio2貢獻率較低,分別為2.8%、1.3%和1.2%,其中降水累計貢獻率為45.6%,溫度累計貢獻率為10.3%,由此可見,降水、海拔對紫果云杉的分布影響起主導作用,其次是溫度。

圖3 Maxent模型中對環境變量重要性的刀切法檢驗

選取降水變量中的年平均降水量bio12、溫度變量中的最冷季度平均溫度bio11及海拔alt進行環境變量響應曲線分析,探討各環境變量對紫果云杉的分布之間的生態學聯系。根據本次適生區預測結果,本次研究主要關注存在概率大于0.4時環境變量的變化情況。如圖4所示,海拔(alt)、年平均降水量(bio12)、最冷季度平均溫度(bio11)均具有顯著的峰值。

如圖4,A所示,海拔在小于1 000 m和大于5 000 m時紫果云杉存在概率極低(小于0.05),當海拔在2 600~4 200 m時紫果云杉的存在概率大于0.4,當海拔在3 600 m時,紫果云杉存在概率達到最大,為0.68;如圖4,B所示,年均降水量響應曲線反映紫果云杉對年均降水量的耐受范圍明顯小于海拔及最冷季平均溫度,從另一方面說明了年均降水量是限制紫果云杉分布的主要環境變量,當年均降水量在590~820 mm時紫果云杉存在概率大于0.5,最適宜紫果云杉生長,當年均降水量達到712 mm左右時紫果云杉存在概率最大為0.71;如圖4,C所示,最冷季平均溫度為-9.8~-1.4 ℃時紫果云杉最適宜生長,當最冷季平均溫度為-7 ℃時紫果云杉存在概率最大為0.61。

圖4 主要氣候因子的響應曲線

利用GIS技術提取108處紫果云杉分布樣點年均溫(bio1)、1~12月降水量與溫度數據(圖5,數據來自WorldClim全球氣候數據庫),發現紫果云杉分布樣點年均溫在2.8 ℃左右,11月至次年3月溫度均在0 ℃以下。紫果云杉通常10月結果,成熟后的種子經過11月至次年3月的低溫貯藏,長時間的積雪覆蓋和低溫環境更容易促進種子萌發,而低海拔的紫果云杉種群經歷了相對短的積雪覆蓋和低溫環境,種子萌發能力較低,可能導致幼苗成功建植和存活的概率較低。

圖5 108處紫果云杉分布樣點氣溫、降水月均變化

2.4 未來氣候情景下紫果云杉的潛在適生區

本次研究基于當前紫果云杉的分布預測結果,對未來不同氣候情景下的適生區分布范圍進行了預測,如圖6所示。在3種氣候情景下,紫果云杉在未來2041-2060年、2081-2100年的分布范圍與當前相比呈現向周邊擴張的趨勢,擴張范圍主要分布在當前適宜生境的西、北方向。在3種未來氣候情景下,2081-2100年的分布范圍比2041-2060年更大,未來2個時期新疆博樂市、西藏日喀則新增了一部分適宜生境,陜西、山西、河南、湖北等地新增的部分潛在分布區與當前環境下的零星分布連接,形成了該地區未來氣候條件下的大面積潛在適生區。與當前和2041-2060年的潛在分布區相比,2081-2100年四川南部、云南地區潛在適生區向南擴張,其中尤以ssp585氣候情景下的潛在分布格局最為明顯。另外,未來氣候條件下紫果云杉適宜生境丟失較少,ssp126氣候情景下僅有2041-2060年西藏噶爾、日喀則、林芝地區零星減少,ssp370氣候情景下2041-2060年、2081-2100年在西藏南部地區減丟失生境最多,ssp585氣候情景下2081-2100年日喀則、林芝地區零星減少。

圖6 未來不同氣候情景下紫果云杉潛在適生區預測

統計不同氣候情景下潛在適生區面積及面積變化,如表3所示。在相同氣候情景下,2081-2100年各適生區面積最大,面積增長百分比遠高于2041-2060年。相對當前氣候情景下的適生區面積,未來2041-2060年、2081-2100年潛在分布區面積將不斷增加,尤其是在ssp370氣候情景下2081-2100年潛在分布區面積最大,為109.140 4 萬km2,其中最適生區面積增長164.60%,增長幅度最大,低適生區面積增長百分比為40.73%,但增長面積最大,是該時期潛在分布區面積增長的主要原因。

表3 氣候變化下紫果云杉潛在適生區面積與當前相比面積變化百分比

在3種氣候情景下,持續發展(ssp126)模式下溫度和降水增長幅度最小(圖7),紫果云杉分布范圍與面積增長百分比最小,局部發展(ssp370)模式下降水量增幅最大,紫果云杉潛在分布區面積增長最大,這一現象也說明了降水是主導紫果云杉分布的重要因素。

圖7 未來氣候情景下108處紫果云杉分布樣點bio11和bio12變化情況

3 討 論

3.1 MaxEnt模型模擬預測結果的可靠性

本次研究模擬的紫果云杉集中分布范圍與李寧寧等人的研究[18]結果一致,垂直分布范圍與中國植物志記載的2 600~3 800 m及陳育峰提出紫果云杉在大渡河中、上游海拔3 000 m~3 800 m或更高的區域分布[25]相吻合,年降水量與最冷季平均溫度的閾值研究與李賀等的研究結果[26]具有較高的一致性,并且本次預測結果經過ROC曲線精度檢驗,AUC值高于0.9,表明本次模型對紫果云杉分布區的預測效果較好,具有很高的可信度,但預測生成的部分分布區與紫果云杉實際分布區域仍舊存在一定差異。本次預測結果顯示,紫果云杉低適生區分布面積較廣,尤其是新疆博樂市西北部與伊寧市北部、山西、河南、湖北等地有零星低適生區分布,此預測范圍大于目前已知的紫果云杉的分布范圍。

MaxEnt模型通過物種分布的不完整信息找到物種分布規律的最大熵[6],從而對物種的存在概率進行預測,確定物種的分布區。一般而言,樣本數量及樣本在空間上的分布狀況直接關系到模型模擬預測結果的精確度與可靠性[23],樣本量越多模型預測結果越準確,但是研究表明MaxEnt模型規則化的程序在小樣本量模擬過程中能夠抵消過度擬合的趨勢,表現出高質量的預測結果[27];從樣本分布范圍來看,樣本點在空間上的分布差異會影響預測模擬過程的擬合狀況,包括水平尺度上的分布與垂直分布,樣本的不均等分布會導致預測結果過度擬合,造成模擬預測的物種分布遠大于物種實際分布范圍,樣本點分布均勻且與研究區域尺度一致,可以避免由樣本問題導致的預測結果的偏差。另外,選取的參加模型模擬預測過程的環境因子是預測結果可靠與否的關鍵。植物在地理上產生分布差異受多種影響因子制約,如溫度、降水、地形、土壤、微生物、水文狀況及其他外在因素,在模型模擬中導入的已去除自相關關系的環境因子越多,建立的約束條件越多,預測結果就越準確[24],若在環境因子導入模型前未去除具有自相關關系的環境變量則會導致模擬結果過度擬合的情況,預測結果將遠大于實際分布。

本次研究在選取的樣本點在導入模型預測前已經過數量及空間分布的篩選,避免了樣本點的選取對預測結果帶來的偏差,但仍舊產生預測范圍大于實際分布的原因主要可能有:(1)本次研究主要選取經過相關性篩選、包含溫度、降水、地形的3種環境因子導入模型對分布區進行模擬,未考慮土壤、微生物、水文條件等其他環境因子,導致分布范圍大于實際分布范圍;(2)林木在實際生長過程中還受到種間關系及一些外在的不確定因素影響,如種間競爭、林木病蟲害、森林火災、人為擾動[17],及物種遷徙、擴散、自然更新能力等因素的影響[10],模型預測模擬時未考慮到這部分影響導致模擬結果過大;(3)林木的實際分布還受制于土地利用類型,模型模擬預測結果是基于氣候條件及地形條件形成的滿足林木生長的基本條件,在實際分布中,林木主要集中分布于林地,而有些高原山地雖然滿足林木生長的氣候及地形條件,但實際卻是耕地或建設用地等其他土地利用類型,并不適于林木生長。總體而言,本次預測結果具有很高的可信度,模擬預測生成的最適生區及較適生區分布范圍對以紫果云杉為主的森林生態系統的生物多樣性保護具有重要意義。

3.2 影響紫果云杉分布和適生性的關鍵因子

根據MaxEnt模型Jackknife、環境因子響應曲線的分析可知降水(主要是年均降水量)及海拔是限制紫果云杉地理空間分布差異的主要環境因子,與李寧寧等的研究結果認為影響紫果云杉分布的主要氣候因子是最暖季均溫和最暖季降水量不同。根據本次研究Pearson相關分析,海拔與最暖季均溫顯著相關(r= 0.934,P<0.01),年均降水量與最暖季降水量顯著相關(r= 0.898,P<0.01),造成主導紫果云杉分布差異的環境因子不同可能主要是由于在模型預測初期篩選環境因子時保留的環境因子不同導致的。

物種響應曲線描述了環境變量與物種存在概率之間的關系,展現了物種分布區對環境的耐受能力[5],適宜紫果云杉生長的海拔閾值為2 600~4 200 m,峰值為3 600 m。海拔作為影響植物生長的綜合環境因子,溫度、水分、光照等環境因子均會隨著海拔的升高而受到影響,進而對植物分布范圍產生影響。已有研究表明海拔是影響紫果云杉種子萌發能力的主要原因[14],而種子萌發是影響植物幼苗存活及未來生長的關鍵階段,是導致植物在地理空間分布上存在差異的重要原因。王桔紅等的研究表明紫果云杉種子萌發率與萌發速率隨海拔升高而增大,且低溫貯藏是打破紫果云杉種子休眠的有效方法之一,紫果云杉種子最大萌發率為海拔3 500 m采集、低溫貯藏160 d的種子[14]。紫果云杉集中分布于高海拔地區,種子需要長時間的積雪覆蓋和低溫環境來萌發,在長期環境選擇的壓力下,逐漸演變為在特殊環境條件下的一種生態適應[28],使得低海拔地區種子萌發育苗成功率較低,造成紫果云杉對高海拔具有獨特適應性的重要原因。

本研究生成的響應曲線中,降水量生態幅較其他兩種環境變量窄,表明紫果云杉對降水量的耐受范圍較小,地理分布范圍受降水量限制較大。紫果云杉作為耐陰性極強、淺根性的樹種,4月底氣溫回暖,降水量逐漸增多,紫果云杉進入生長季,氣溫處于7~12 ℃之間,降水量在80~110 mm之間(圖5),適宜的溫度使植物光合作用、蒸騰作用加強,充足的降水能夠保證土壤中供給樹木的有效水分充足,利于紫果云杉的徑向生長,若在該時期氣溫過低則蒸騰作用、根系生長受到影響,抑制植物養分吸收,降水量過高則導致土壤通氣不良,抑制根對水分的吸收、運轉和疏導,阻礙根系呼吸和養分吸收,光合作用受到抑制,不利于樹木生長。另外,紫果云杉分布范圍內,當年11月至次年3月氣溫在-1~-7 ℃之間,降水量在20 mm以下,較高的最低氣溫能夠有效降低紫果云杉根系、枝葉等遭受凍害的風險,從而較少的消耗樹內能量,最大限度保持生長活性[29]。若在該時期若氣溫過低、降水量過高,則會使凍土情況嚴重,地下部分的根系很難從土壤中吸收水分,致使植物水分失衡,枝條干枯;或會凍裂主干,樹皮脫落或卷折,招致病蟲害,削弱樹木的生存能力;或導致凍拔現象使根系外露,嚴重時導致樹木倒伏死亡。

3.3 未來氣候變化情景下紫果云杉的潛在適生區變化

在未來氣候變化情景下,紫果云杉的分布范圍主要向西、北方向延伸,適宜生境減少的區域主要在西藏南部,適生區面積總體增加,且2081-2100年面積增長遠高于2041-2060年。分析108處紫果云杉分布點未來3種氣候情景下2041-2060年、2081-2100年最冷季度平均溫度(bio11)、年平均降水量(bio12)數據,在未來氣候情景下,氣候變暖,相較于當前氣候環境,最冷季度平均溫度、年平均降水量不斷升高,結合前文研究,認為這可能是導致紫果云杉適生區擴張的主要原因。已有研究表表明,氣候變暖可能導致高緯度濕潤中緯度更加干旱[22,30],這可能是紫果云杉適宜生境向西、北方向延伸,西藏南部適宜生境有部分減少的重要原因。

4 結 論

紫果云杉在當前氣候環境條件下分布區主要集中在青藏高原東緣地區的四川、甘肅、青海、西藏等地,四川西北部、甘肅南部、青海東部和西藏東部林芝、昌都地區是最適宜紫果云杉生長的地區。影響紫果云杉分布的主要環境變量為降水量和海拔,適宜紫果云杉生長的海拔、年降水量和最冷季平均溫度的閾值分別為2 600~4 200 m、590~820 mm和-9.8~-1.4 ℃,峰值分別為3 600 m、712 mm和-7 ℃。未來氣候變化情景下,紫果云杉的分布范圍主要向西、北方向延伸,西藏南部小部分生境丟失,總體呈現適宜生境面積增加的趨勢,2018-2100年局部發展模式下適生區面積達到最大,為109.1404 萬km2。

猜你喜歡
物種模型
物種大偵探
物種大偵探
一半模型
吃光入侵物種真的是解決之道嗎?
英語世界(2023年10期)2023-11-17 09:18:18
重要模型『一線三等角』
重尾非線性自回歸模型自加權M-估計的漸近分布
回首2018,這些新物種值得關注
電咖再造新物種
汽車觀察(2018年10期)2018-11-06 07:05:26
3D打印中的模型分割與打包
FLUKA幾何模型到CAD幾何模型轉換方法初步研究
主站蜘蛛池模板: 国产亚洲欧美在线人成aaaa| 国产一国产一有一级毛片视频| a级毛片免费播放| 亚洲精品人成网线在线| 亚洲福利视频网址| 麻豆精品视频在线原创| 久久精品欧美一区二区| 91久久青青草原精品国产| 美女扒开下面流白浆在线试听 | 9啪在线视频| 国产白浆在线观看| 日本AⅤ精品一区二区三区日| 91精品小视频| 国产黄网永久免费| 麻豆精品在线视频| 国产精品一区二区无码免费看片| 丁香五月婷婷激情基地| 国产一级裸网站| 99re免费视频| 一本色道久久88综合日韩精品| 国产成人成人一区二区| 国产精品丝袜视频| 欧美日韩国产精品va| 成人中文字幕在线| 国产一区二区三区在线精品专区| 又黄又湿又爽的视频| 欧美精品成人一区二区在线观看| 青青草原国产免费av观看| 无遮挡国产高潮视频免费观看| a级毛片一区二区免费视频| 亚洲视频无码| 在线观看国产一区二区三区99| 亚洲成a人片| 四虎永久免费地址在线网站 | 国产毛片基地| 国产福利一区在线| 日韩a级片视频| 久久大香香蕉国产免费网站| 色婷婷成人网| 国产又大又粗又猛又爽的视频| 色综合手机在线| 国产成人精品一区二区| 亚洲不卡av中文在线| 欧美国产在线一区| 亚欧成人无码AV在线播放| 国产精品尤物铁牛tv| 婷婷色在线视频| 亚洲成人播放| 网友自拍视频精品区| 无码高潮喷水在线观看| 91在线一9|永久视频在线| a国产精品| 午夜日本永久乱码免费播放片| 91久久国产成人免费观看| 无码AV日韩一二三区| 日韩欧美国产精品| 狠狠v日韩v欧美v| 久久成人18免费| 91丨九色丨首页在线播放 | 青青国产在线| 国产精品私拍99pans大尺度| 91激情视频| 免费国产一级 片内射老| 亚洲欧美精品日韩欧美| 国禁国产you女视频网站| 亚洲欧美综合另类图片小说区| 一本大道东京热无码av| 青青草原国产| 国产三级国产精品国产普男人| 亚洲高清中文字幕| 不卡无码网| 免费看a毛片| 日韩精品专区免费无码aⅴ| 九九热精品视频在线| 午夜福利网址| 特级毛片免费视频| 欧美精品H在线播放| 在线精品自拍| 露脸一二三区国语对白| 在线一级毛片| 亚洲国产成熟视频在线多多| 国产午夜看片|