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

錫林郭勒盟植被物候枯黃期對干濕變化的時間多尺度響應

2022-01-21 02:11:50佟斯琴
中國環境科學 2022年1期
關鍵詞:研究

呂 達,包 剛,2,3*,佟斯琴,2,雷 軍

錫林郭勒盟植被物候枯黃期對干濕變化的時間多尺度響應

呂 達1,包 剛1,2,3*,佟斯琴1,2,雷 軍1

(1.內蒙古師范大學地理科學學院,內蒙古 呼和浩特 010022;2.內蒙古師范大學,內蒙古自治區遙感與地理信息系統重點實驗室,內蒙古 呼和浩特 010022;3.寧夏大學,西北土地退化與生態恢復省部共建國家重點實驗室培育基地,寧夏 銀川 750021)

基于全球監測與建模組(GIMMS)1982~2015年第三代歸一化植被指數(NDVI) GIMMS NDVI 3g數據集和氣象觀測數據,采用累計NDVI的Logistic曲線曲率極值法提取錫林郭勒草原植被枯黃期,并結合不同時間尺度標準化降水蒸散指數(SPEI)分析植被枯黃期對干濕變化的響應特征.結果表明:1982~2015年錫林郭勒草原植被平均枯黃期主要集中在265~280d,由東北向西南推遲,整體呈現提前趨勢(-0.14d/a).年SPEI傾向率總體呈弱下降趨勢(-0.10d/10a),其中上升趨勢主要分布在東北部草甸草原,下降趨勢主要分布在西南荒漠草原.近34a來,錫林郭勒草原干旱發生頻率高,輕、中旱面積占比較大,干旱突變始于20世紀90年代.時間多尺度SPEI影響EOS空間格局變化差異顯著,錫林郭勒盟植被EOS受月尺度SPEI-1正面影響較大,短時間尺度干濕環境向好將延遲植被枯黃期;SPEI-3、SPEI-6、SPEI-9及SPEI-12與EOS負相關像元比例較高,植被受干濕環境長期變化的影響提前進入枯黃期.時滯分析表明,錫林郭勒盟SPEI-1的滯后效應主要發生在季前2個月,長時間尺度SPEI(SPEI-3、SPEI-6和SPEI-9)累積至6月的滯后效應尤為顯著,特別在SPEI-3顯著相關面積達到最大.

錫林郭勒盟;枯黃期;標準化降水蒸散指數;氣候響應

秋季植被枯黃期(EOS)是植被生長季的結束日期,它是陸地生態系統對全球氣候變化響應的重要指標之一,其變化與陸地水熱條件密不可分,并通過調節陸地和大氣能量交換對氣候產生明顯的反饋作用[1].研究表明,由于全球氣候變暖,北半球中高緯度大多地區的植被EOS呈推遲趨勢,從而延長了植被年生長天數,生長期增加[2].干旱是擾動生態系統的主要因素,自20世紀下半葉以來,干旱事件的頻度、時間和影響范圍大幅增加,引起土壤水分脅迫的負面影響,阻礙植被生長和生理功能[3],長時間影響著生態系統碳匯增量.目前廣泛使用干燥度(drought index)來表征氣候干濕狀況,標準化降水蒸散指數(SPEI)[4]融合帕爾默干旱指數[5]對蒸散的敏感性和標準化降水指數[6]多時間尺度特征等優點,在區域尺度上能夠較好反映干濕狀況[7-8],被廣泛用于分析干旱對植被動態的影響.

研究表明,植被的枯黃與干濕變化關系密切.如Guo等[9]研究發現,氣候變化與植被枯黃期的變化密切相關,干旱環境下草本植物更易受到影響.Kang等[10]利用SPEI論證了中國北方半干旱生態系統植被物候和生產力對干旱擾動的響應,發現秋季干旱在降低植被生產力的同時,也顯著提前草地和貧瘠/稀疏植被土地的枯黃期.然而,黃文琳等[11]在內蒙古荒漠草原區的研究表明,秋季氣候環境干燥將推遲草原植被EOS,同年氣候環境濕潤則會導致EOS天數提前.EOS對環境干濕水平的響應機制具有復雜且不對稱的時間效應,特別受時滯和累積作用的綜合影響(即植被在連續的前幾個月中受到氣候條件的顯著影響)[12].近期的研究表明,植被早期的枯黃可能不是由當前水熱條件所驅動,而是先前的干濕變化驅動所致.干旱事件發生后,植被通常需要長時間來修復受損根系并恢復之前的生長能力[13].如Cui等[14]對加拿大草地的研究發現,枯黃期前干旱主導區域變化面積占比為26%~44%,在植被生長的關鍵階段累積降水量減少,最終導致草地植被生長季度縮短.以往研究大多是在單一時間尺度(如某月、季節或年尺度)分析干濕變化對植被物候的影響,而對物候在時間多尺度干濕變化的時滯和累積效應有所忽視,導致草地植被EOS及其驅動機制的解析很不確定.

錫林郭勒草原地處中國北方干旱半干旱氣候區,生態環境脆弱,極易受到干旱等氣象災害的影響,屬蒙古高原南段環境變化敏感地帶[15].近年來,在氣候變化、土地利用變化和社會經濟轉型的綜合影響下,錫林郭勒草原植被生長結構、組成、功能等均發生變化,草地沙化及鹽漬化面積呈進一步擴大的態勢[16].此外,水分的虧缺抑制植被養分的積累并加速葉片衰老,影響EOS的時間變化[17].鑒于此,為了更好地了解EOS對干濕變化的時間響應特征,本文以錫林郭勒盟為研究區,基于第三代全球檢測與建模組(GIMMS)的歸一化植被指數(NDVI)GIMMS NDVI 3g數據,利用時間序列諧波分析法(HANTS)對數據進行預處理,采用累計NDVI的 Logistic曲線曲率極值法提取EOS.并結合同期氣溫降水數據,采用多時間尺度SPEI指數進行干旱程度和干旱時間的量化,分析EOS對干濕變化的時間多尺度響應特征.研究結果可為錫林郭勒盟草原植被保護、沙化防治提供參考.

1 材料與方法

1.1 研究區概況

圖1 研究區植被類型及氣象臺站分布

錫林郭勒盟位于內蒙古自治區中部(43°02′~ 44°52′N, 115°13′~117°06′E),是中國最大的天然牧場和畜牧業基地,可利用草地面積達1.96×105km2,占全盟總面積的96.55%[15](圖1).該區屬中溫帶半干旱、干旱大陸性季風氣候,四季分明,具有干旱少雨、風多風大的氣候特點,區域年平均溫度在0~3℃之間波動,自東向西逐漸升高,降水的空間分布特征與溫度相反,年降水量300~380mm,雨季多集中在7~9月,蒸發量在1500~2700mm之間[18].錫林郭勒盟地貌類型主要為高平原,地勢南高北低,平均海拔800~1800m,南接陰山北麓,是京津冀北部的天然屏障[19].近幾十年間,在過度打草放牧、采礦填挖、旅游開發等不合理資源利用和多變的氣候要素聯合驅動下,區域植被退化、生態承載力下降、沙塵暴肆虐、次生災害頻發,嚴重危及京津乃至華北等地區的生態系統安全[20].自2000年啟動實施了退牧還草、風沙源治理等生態工程以來,區域生態有所改善.

1.2 數據來源及預處理

NDVI時間序列數據為美國國家航空航天局(NASA)提供的1982~2015年GIMMS NDVI 3g數據集產品(http://glcf.umd.edu/data/gimms), GIMMS NDVI 3g數據來自NOAA衛星的AVHRR傳感器,時間和空間分辨率分別為15d和0.0833°,每年共24幅數據影像,覆蓋全球范圍.該數據集經過相關校正,最大程度地消除了衛星傳感器、大氣對數據產生的誤差,具有較好的數據一致性和可信度,已廣泛應用于識別植被長勢和物候動態研究[21-22].本文為得到錫林郭勒盟植被枯黃期數據集,在此基礎上對原始數據做格式轉換、掩膜等預處理,以對應各月尺度的干旱數值.

氣象數據下載自中國氣象數據共享網(http: //data.cma.cn),包括錫林郭勒盟15個氣象站點(圖1)1982~2015年月平均溫度和總降水量數據,主要采用均值替代法對個別臺站的缺測值及異常數據進行插補處理[23].為提高插值精度,本研究利用基于部分薄板樣條(PTPSS)函數模型專業氣候插值軟件ANUSPLIN,引入各氣象站點數字高程模型數據(DEM)作為協變量,并調用splina和lapgrd程序對34a月平均氣溫和月總降水量數據進行插值處理,獲得與遙感數據匹配的逐月氣象要素柵格數據[24].

1.3 研究方法

1.3.1 NDVI數據的平滑與物候參數的提取 為進一步減少噪聲和保持數據的周期性,利用HANTS時間序列平滑法[25]對數據進行平滑重構(圖2a),提高時序數據質量,使其能更好地反映植被動態變化.HANTS方法的核心是最小二乘法的迭代擬合和傅里葉變換的分解重構.其表達式為:

式中:A0為諧波的余項;N為時間序列長度;Aj、wj、qj為各諧波的振幅、頻率、初相位;m為諧波個數,等于-1.采用Hou等[26]提出的累計NDVI的Logistic曲線曲率極值法提取EOS.該方法是HANTS平滑后的年內累積NDVI時間序列進行擬合(式2),計算Logistic擬合曲線的曲率(式3,4),將曲率變化率極小值點對應時間定義為EOS的開始日期(圖2b).

1.3.2 標準化降水蒸散指數(SPEI) 本文利用標準化降水指數(SPEI)量化了錫林郭勒盟1982~2015年干濕變化特征及趨勢.該指數根據降水量與潛在蒸散量(PET)的差值偏離平均狀態的程度計算水分盈虧[27],根據SPEI值大小,能夠判斷出區域的干濕水平,因其多時間尺度的特點,特別適于監測干旱及其對植被長勢的影響.本研究根據研究區地理特性[28],并參考GB/T 20481—2006《氣象干旱等級標準》[29]將SPEI劃分為5個等級,詳見表1.具體計算步驟如下:

第1步,計算潛在蒸散.Vicente-Serrano采用的是Thornthwaite方法.

式中:T為月平均溫度,℃;為年熱量指數;為常數.

第2步,計算逐月降水量與潛在蒸散發量的差值:

式中:表示1~12月;D為降水與蒸散的差值;P為月降水量,mm;PET為月蒸散量,mm/月.

在不同時間尺度下構建水平衡累積序列:

式中:為時間尺度,一般為月;為計算次數.

第3步,對D數據序列進行正態化.由于原始數據序列D中可能存在負值,所以SPEI指數采用3參數的log-Logistic概率分布,得到概率分布的累積函數:

式中:參數、分別為尺度參數、形狀參數和位置參數,如式9,10,11.

對累積概率密度()進行標準正態化,當累積概率£0.5時:

當>0.5時:

本文為探索干濕變化對植被枯黃期時滯及累積效應,時間多尺度干濕特征表達如下,主要分析1,3,6,9和12個月時間尺度的SPEI,分別記作SPEI- 1、SPEI-3、SPEI-6、SPEI-9和SPEI-12.SPEI-時間尺度表示累積月份向前延續-1個月,其中SPEI-1表示月干濕程度,可以反映枯黃期當月旱澇細微變化,SPEI-3表示季節干濕程度,SPEI-6和 SPEI-9表示累積6個月和9個月的干濕程度, SPEI-12表示年際干濕程度[30].例如9月份的SPEI-12是根據上年10月到9月共12個月的降水量與蒸散量插值聚合得到,用以反映區域植被枯黃期前的整年干濕狀況.此外,文中干旱站次比表示發生不同等級干旱的站點數占研究區所有站點數的比值,用于量化干旱影響范圍的大小.

1.3.3 趨勢分析和相關分析方法 采用一元線性回歸趨勢分析方法計算EOS和SPEI 指數傾向率,以綜合表征研究區多年空間格局演變特征.其傾向率計算公式為:

表1 SPEI干濕等級劃分

式中:R表示EOS和SPEI之間的相關系數,cor代表皮爾遜相關函數,為不同時間尺度(1,3,6,9和12月),max表示EOS對應的各尺度干旱中的最大相關系數.

2 結果分析

2.1 植被EOS和SPEI空間分布及變化趨勢

圖3 1982~2015年錫林郭勒盟平均枯黃期、枯黃期傾向率及SPEI傾向率

表2 在氣象站點計算的多年平均枯黃期

2.1.1 植被EOS與SPEI空間變化特征 由圖3a和表2可以看出,研究區植被平均枯黃期主要集中在265~280d(95.52%),總體上由東北向西南推遲,即9月下旬至10月上旬.枯黃期較早的地區主要分布在蘇尼特左旗中部、東烏珠穆沁旗北部草甸草原區、西烏珠穆沁旗、正藍旗、阿巴嘎旗大部地區,約發生于267d左右(68.29%).而蘇尼特右旗等荒漠草原枯黃期較晚,約發生于270d后(31.71%).植被平均枯黃期由小到大依次跨越草甸草原、典型草原再到荒漠草原逐漸變化,從東北向西南平均枯黃期呈推遲趨勢,總體與該地的植被類型分布相對應(圖1).同樣,薩日蓋等[31]、劉紅海等[32]、烏日汗等[33]利用不同遙感數據源在內蒙草地研究表明,草原生態區枯黃期早于荒漠生態區,與本研究結果基本一致.

從34a錫林郭勒盟枯黃期總體變化趨勢看(圖3b),大部分像元平均變化趨勢為-0.14d/a,整體上呈現提前趨勢.約占研究區89.8%的像元呈現提前趨勢,主要分布在錫林浩特市、西烏珠穆沁旗和東烏珠穆沁旗的接壤地帶等大部分研究區域.10.2%的像元則呈現推遲趨勢,零星分布在阿巴嘎旗以北、太仆寺旗以南、蘇尼特右旗中東部等區域.其中在中東部地區枯黃期的提前趨勢變化較大,整體提前-0.21d/a以上.結合圖3a和圖3b可以發現,錫林郭勒盟枯黃期變化趨勢的空間格局差異性明顯,除太仆寺旗南部及蘇尼特右旗小部分區域外,整體上枯黃期提前趨勢呈現小面積多點狀零星分布且數值波動性較大,表明研究區枯黃期變化可能受多種因素共同的作用,除溫度和降水外可能還存在著更復雜的響應機理.

由圖3c可知,錫林郭勒盟SPEI變化趨勢為-0.47~0.19d/10a,由東北向西南逐漸遞減,有較明顯的水平地帶性規律,區域上總體呈弱下降趨勢(-0.10d/10a),下降和上升趨勢分別占總面積的59.25%和43.75%.上升趨勢平均變化速率為0.17d/ 10a,主要分布在東烏珠穆沁旗和西珠穆沁旗的東北地區,此區域水熱條件較好,即草甸草原總體上越來越濕潤.下降趨勢平均變化率為-0.38d/10a,主要分布在蘇尼特左旗與二連浩特市、蘇尼特右旗交界處,表明荒漠草原總體上越來越干旱,且干旱化速率超過了草甸草原.值得注意的是研究區中部典型草原SPEI增長速率為-0.18d/10a,說明近34a此區域干旱有逐漸增加趨勢,且空間分布上干旱速率由北向南逐漸上升.

圖4 1982~2015年錫林郭勒盟不同時間尺度SPEI動態特征

2.1.2 時間多尺度的SPEI變化特征 從圖4可以看出,錫林郭勒盟地區1982~2015年時間多尺度的SPEI值波動敏感性明顯不同.干旱指數的波動隨時間尺度的增大,干濕變化越平緩,只有在連續降雨、干旱或高溫等環境下才會使之發生改變.相反,時間尺度越小,受月氣溫降水的影響,干濕波動頻率變化越顯著.由于受短期干旱的影響,SPEI-1沿0值劇烈波動,反映出錫林郭勒盟各月旱澇的頻繁交替. SPEI-3和SPEI-6的波動周期相對較長,波動減弱,顯示出錫林郭勒盟干濕季節性變化規律.SPEI-12波動相對穩定,可以較為清晰地反映出34a錫林郭勒盟水分多寡的年際變化特征.1982~1999年期間結果表明,錫林郭勒盟的干旱發生頻率較低,持續時間短.2000年后錫林郭勒盟進入干旱頻發期,特別在2002, 2006, 2008, 2010和2012年開始出現了重度干旱,并且干旱事件持續了較長的時間.

為進一步探究干旱影響范圍的大小,分別對研究區34a干旱站次比進行分析統計(圖5a).結果表明錫林郭勒盟年干旱站次增加趨勢明顯,干旱面積總體上表現為:輕度干旱>中度干旱>重度干旱>極端干旱.1982~1998年干旱站次比為34%,區域受干旱影響范圍較小,其中1983年和1995年區域干旱占比均超過了66%,僅1990年發生一次全區域干旱.2000年起干旱頻率呈現遞增趨勢,尤其在2001~2003年干旱面積達到峰值,中度以上干旱約占總面積的60%.2006年干旱范圍占總面積的85.6%,其中重度干旱、極端干旱分別占總面積的29.8%和13.3%; 2008年干旱面積以重度干旱為主,相較于2006年增加了23.5%,其中中度干旱面積減少了13.3%;2010年受干旱影響范圍占總面積的86.6%,干旱面積比例呈減少趨勢,極端干旱面積減弱,說明錫林郭勒盟生態逐步回緩,可能與當地治沙、休牧等生態措施有關[15,19].

根據先前對SPEI值的變化分析,在不同的時間尺度上,本研究表明1999年左右的SPEI值發生了很大的變化,采用Mann-Kendall突變檢驗法對SPEI值進行突變檢驗(圖5b).錫林郭勒盟SPEI值呈下降、上升、再緩慢下降的變化趨勢,波動在-1.58~1.38之間,整體趨勢呈現不顯著減小(-0.036/a).UF值在1982~1989年間為正,之后緩慢上升,其間UF曲線整體未超出置信區間,表明研究區濕潤化趨勢不顯著.在顯著性水平1.96的閾值內,UF和UB的交點位于1999年,說明1999年SPEI發生了突變,錫林郭勒盟自1999年旱情加劇,由先前的濕潤轉為大面積干旱.2000年后UF曲線整體小于0,并在2013年超出了置信區間,干旱化趨勢顯著.

2.2 錫林郭勒盟時間多尺度干濕變化對植被EOS的影響

短時間尺度的SPEI與EOS相關性越大表明植被EOS受干旱的響應越迅速;相應地,長時間尺度的SPEI與EOS相關性越大則表明植被枯黃期受干旱累積效應越明顯,由于不同時間尺度的干旱累積,最大響應空間分布具有明顯的地帶分異性.由圖6a可以看出,近34a單月尺度(SPEI-1)的干濕水平下,SPEI-1大面積推遲了植被EOS.具體而言,SPEI-1與EOS的最大相關系數像元分布主要呈正相關(均>75%),其中多數地區介于0.3~0.6之間,像元值最大超過0.7.二者相關性在整個研究區大部區域反映較強,僅在東烏珠穆沁旗北部、蘇尼特左旗西部和蘇尼特右旗接壤地帶反映較弱,可能是因為錫林郭勒盟以北小部分區域為濕地生態系統,而中西部地區則多為荒漠生態系統,區域土壤含水量變異性小且較為穩定[35],SPEI-1響應不敏感.

圖5 1982~2015年錫林郭勒盟干旱影響范圍和SPEI突變年份M-K檢驗

UF表示標準正態分布,即正向時間序統計量曲線,UB為逆序統計量曲線

從圖6a與圖6b對比可知,隨著干旱累積時間尺度的增加,最大相關性面積由先前的正相關轉為負相關,且SPEI-3與EOS負相關面積(74.89%)遠高于正相關面積(25.11%),表明EOS前干濕作用累積時間尺度的增加對植被枯黃具有提前作用.尤其在阿巴嘎旗東南部、蘇尼特左旗、正鑲白旗及正藍旗大部分地區像元一直維持著負相關高值態勢(79.01%) (圖6b),該地區植被類型多為典型草原,其水熱條件較好,植被生長受干濕環境響應較強.這印證了先前的相關研究[13],植被的提前枯黃可能受先前的水熱條件影響[22],長期的干旱通過碳水交換改變了秋季物候的形成,并隨著時間尺度的增加,植被提前枯黃.

半年累積尺度上(圖6c),EOS與SPEI-6最大相關系數負相關面積較3個月尺度增加15.5%,但相關性減弱,僅在阿巴嘎旗西部東烏珠穆沁旗中部及西烏珠穆沁旗中部17.5%的區域呈正相關,表明錫林郭勒盟植被生長季節(4~9月)的水熱變化大面積提前了植被的生長.EOS與SPEI-9的最大相關系數像元空間分布與SPEI-6的空間格局逐漸趨于一致(圖6b),其中72.27%的區域呈負相關,且負相關均在-0.3~-0.4之間,總體上,區域內長期的干旱累積致使植被較早地進入枯黃期.已有研究顯示,12個月尺度的SPEI通常與地下水分含量及徑流補給密切相關[36],從圖6e可以看出,EOS與SPEI-12最大相關系數空間分布與枯黃期的平均分布相似,較大的相關系數分布在海拔相對較高的地區.

圖6 1982~2015年錫林郭勒盟EOS與各尺度SPEI最大相關系數的空間分布

a、b、c、d、e 分別表示EOS與SPEI-1、SPEI-3、SPEI-6、SPEI-9和SPEI-12的相關性

為更加全面揭示錫林郭勒盟植被EOS對干旱響應程度,將各相關系數頻率做三維統計處理(圖7),隨著SPEI累積時間尺度的增加,錫林郭勒盟EOS與SPEI的相關頻率幅度發生改變,正相關頻數峰值由先前的7.2附近(SPEI-1)波動逐漸降到了2.2附近(SPEI-6),隨著累積尺度的增加,峰值響應逐漸趨于平穩.正相關系數呈負增長的趨勢,由先前的0.45逐漸降到了0.2.負相關峰值則由先前的2.8附近(SPEI-1)波動逐漸上升到6.8附近(SPEI-6).由圖7可知,錫林郭勒盟SPEI對EOS正相關遠大于負相關的時間尺度為1個月,其次是12個月;負相關則主要集中在3, 6, 9個月時間尺度.表明錫林郭勒盟植被枯黃期受1個月尺度的干濕狀況影響敏感,隨著時間尺度的增加,植被枯黃期對環境干濕響應敏感性降低.

2.3 錫林郭勒盟植被EOS受時間多尺度干濕變化的時滯響應

圖7 錫林郭勒盟EOS與不同時間SPEI最大相關系數頻率分布

從上文可以看出不同時間尺度干濕變化對錫林郭勒盟植被EOS具有一定的累積效應,由于植被生理過程受環境水熱條件的制約,在不同時間尺度干濕變化影響下植被EOS存在滯后性.本文選取植被EOS前1~9月各尺度SPEI的不同滯時與EOS進行相關性分析,得到了最大響應滯后期的空間格局(圖8)及其顯著相關面積(圖9),以探討各尺度干濕變化對植被EOS的影響.

由圖8a所示,1個月尺度的干濕程度對植被枯黃期影響最大,其中滯后2個月的SPEI-1與EOS呈正相關(69.15%),且通過顯著檢驗的區域面積占比為57.56%(圖9).呈正相關的區域主要分布在阿巴嘎旗、錫林浩特市、東烏珠穆沁旗及西烏珠穆沁旗等草甸草原地區,而研究區北部的草甸草原地區約滯后4個月.一般1個月尺度的SPEI可以反映當月氣象干濕程度,從不同季前長度來看(圖9),隨著滯后時間的推移,區域干濕程度與植被枯黃期顯著性逐漸減弱,表明錫林郭勒盟植被EOS的干濕程度對物候期推遲具有重要作用.另一方面EOS與滯后2個月的SPEI-1正相關,而與滯后4個月的SPEI-1負相關,表明植被EOS前2個月水熱條件向好對植被生長具有延遲作用,而滯后4個月的氣象干濕程度加劇了植被受水分脅迫的負面作用,導致其提前結束生長[37].

圖8 1982~2015年錫林郭勒盟EOS與各尺度SPEI最大響應滯后期的空間分布

a、b、c、d、e 分別表示EOS與SPEI-1、SPEI-3、SPEI-6、SPEI-9和SPEI-12的相關性,負值代表負最大相關,正值代表正最大相關

3和6個月尺度的干濕程度對植被EOS影響較大(圖8b,圖8c),其中滯后4個月的SPEI-3和SPEI-6與EOS均呈負相關,達到顯著水平的區域面積占比分別為42.82%和36.32%(圖9).以阿巴嘎旗、蘇尼特左旗、正鑲白旗和正藍旗地區達顯著負相關比例最高,約占總面積的69.41%,表明4~6月或1~6月的干濕水平對植被枯黃期的提前具有滯后響應.這可能與植被生長對耗水量的需求有關[38],3、6個月尺度的SPEI常用于檢測生態農業干旱,植被前期的生長消耗土壤較多水分,加之錫林郭勒盟逐年低降水量和高蒸發量帶來的干旱化程度增加[18],導致植被秋季缺水枯黃期提前.與SPEI-6不同的是,滯后1個月的SPEI-3與EOS正相關面積占比為30.88%,主要分布在東烏珠穆沁旗和西烏珠穆沁旗南部的典型草原.緣于7~9月份錫林郭勒盟進入到雨季,季前溫度和降水的增加緩解了植被水分虧缺,進一步減輕了土壤缺水帶來的負面效應,從而延緩了植被的枯黃期[39-40].

9和12個月尺度的干濕程度對植被EOS影響較小(圖8d,圖8e),SPEI-9與EOS最大滯后期是季前4個月,且均呈負相關(55.41%),達到顯著水平的區域面積比例縮小為25.11%;SPEI-12與EOS最大滯后期則是季前1個月,且均呈正相關(31.83%),達到顯著水平面積減小到3.72%(圖9).9、12個月長時間尺度的SPEI通常代表水文干旱指數,水文干旱對植被生長的滯后性不同于氣象干旱敏感[41],從圖直觀反映出植被對干旱的響應較遲緩,且由于研究區大部植被類型均為多年生草本植被,長時間尺度的SPEI對植被EOS的滯后效應不明顯.

圖9 1982~2015年錫林郭勒盟季前不同尺度SPEI與EOS相關系數顯著相關面積比

3 討論

通過對1982~2015年錫林郭勒盟植被GIMMSNDVI 3g數據的分析結果表明,錫林郭勒盟植被枯黃期呈微弱提前趨勢(-0.14d/a),平均枯黃期空間格局表現為由北向南逐漸增大的分布特征,這與前人在內蒙古的秋季物候研究結果較為一致[31].由于截取時間段、研究區選取和物候提取方式的不同,本文結果與個別學者的研究結論存有一定偏差,如Kang等[10]分析了1982~2012年中國北方干旱區秋季枯黃期變化為-0.2d/a,提前趨勢略高于本研究;在數據處理上,前人采用的是閾值法和擬合法,而本文利用年內累計NDVI的Logistic曲線曲率極值法提取枯黃期,該方法經多學者驗證在內蒙古地區有較好的物候擬合效果[31-33].近34a來,不同時間尺度的SPEI顯示錫林郭勒地區自20世紀80年代,干旱事件頻發且持續時間較長(圖4),空間上總體表現為從東北部草甸到西南部荒漠增加趨勢(圖3c),其時空變化與現有研究基本一致[34].錫林郭勒盟地處中國北方干旱半干旱地帶,區域氣候的暖干化趨勢加劇了植被的退化[42],從而導致地表水和熱量條件的改變和地表徑流的增加,這在一定程度上加大了地區的干旱化程度[43].

錫林郭勒盟多時間尺度的SPEI影響植被EOS在空間格局變化上存在明顯差異.較短時間尺度的SPEI-1與EOS在空間上正相關面積占比較大(圖6a),表明短期內區域干濕水平的向好將延緩植被的枯黃,干濕氣候的改善一定程度上緩解水分脅迫的負面效應,減輕了秋季霜凍對植被生理活動的損害[39],從而進一步推遲植被秋季的生長進程[22].此外,受不同植被分布影響,北部及西南部常年濕潤或干旱地區,零星分布著多年生草本和耐旱植被,植被EOS的敏感性隨著環境生態穩定性的增加而降低[44].然而,長時間尺度的SPEI(SPEI-3、SPEI-6、SPEI-9和SPEI-12)與EOS的響應則與之相反(圖6b、c、d、e),隨著累積時間尺度的增加,植被較早地進入枯黃期.這印證了先前的相關研究[13],植被的提前枯黃可能受先前的水熱條件影響[22],盡管氣溫和降水對研究區植被生長起到了促進作用,但長期干旱導致土壤水分改變影響秋季物候的形成.從時間多尺度的干濕變化與枯黃期最大相關系數頻率分布特征來看(圖7),二者在SPEI-1尺度正相關性最強,SPEI-12尺度次之,負相關性最大出現在SPEI-6.表明干濕變化在植被生長的各個階段共同形成了對秋季枯黃期的累積效應.與本文結果類似,在中國北方干旱半干旱區發現,植被物候對氣候變化高度敏感,秋季前3個月的干濕變化將提前植被枯黃期[10],這可能與植被內在的調節機制有關[45],水熱活動強烈促使葉片氣孔較早地關閉[46],進一步加劇植被蒸騰速率和光化學作用帶來的水分流失.另一方面,植被呼吸維持較高的速率而加速碳退化[47],縮短生長周期以應對水分蒸發流失的抵抗力[48],因此,相較于正常年份植被更早停止生長.

通過結合不同時間尺度SPEI對植被EOS進行滯后效應分析得知(圖8, 9),錫林郭勒盟SPEI-1的滯后響應主要發生在EOS前2個月(圖8a),與大區域尺度研究結果相似,如Peng等[49]在北半球監測到46.2%的植被EOS表現出一定的滯后效應,除當月的干濕變化對枯黃期影響較小外,尤其在季前2個月草原和灌叢上檢測到干旱對植被枯黃影響顯著,同降水的時滯效應一致,說明前一階段的水分利用對下一階段的植被生長發育具有間接影響[50].而在3個月時間尺度(SPEI-3)下,本研究結果顯示,除7~9月的干濕水平對植被枯黃期的顯著影響較大外,4~6月的干濕水平對植被枯黃期的顯著影響甚至超過了7~9月.這可能與植被水分吸收和下滲消耗產生的滯后效應有關[50],生長季前期大量的水分利用降低了土壤水含量,且研究區植被多以淺根植物為主[29],土壤水分虧缺無法用于植物碳的合成,導致EOS的提前.隨時間尺度的增加,累積至6月的干濕狀況對植被枯黃期影響減弱,此時SPEI與EOS呈現出較弱的相關性.

此外,研究中利用遙感數據獲取植被枯黃期除其固有的分辨率限制條件外,仍缺乏對物候長期有效的地面驗證工作.同時,本研究中側重于利用降水和潛在蒸散量來計算SPEI,但氣候模型預測表明,多個極端事件(例如,嚴重干旱造成的大范圍變暖)可能比單個干旱事件發生的強度更高[51],因此,有必要利用田間試驗或模型來收集地面物候觀測資料來探索極端氣候環境對植被生長后期的影響,以揭示植被枯黃期對干濕變化的響應機制.另外,考慮到人類活動(如生態恢復計劃,城市擴張和農業灌溉),其他非氣候因素(如二氧化碳肥,氮和磷的沉積)以及自然干擾(例如野火,植物病蟲害)等綜合因素對植被枯黃期空間格局變化產生一定的影響[52],對于后續更準確預測和研究植被動態變化尤為重要.

4 結論

4.1 1982~2015年錫林郭勒盟植被枯黃期主要集中在265~280d,時間波動性較大,89.8%像元呈現提前趨勢(-0.14d/a).空間分布上,錫林郭勒盟東北部草甸草原植被枯黃時間早于西南部荒漠草原.年SPEI值呈弱下降趨勢,下降速率為-0.10d/10a干旱化趨勢整體表現為南部高,北部低,具有較明顯干濕地帶性分布規律,1999年SPEI值發生突變.

4.2 錫林郭勒盟植被枯黃期對短時間尺度干濕變化最為敏感,并且枯黃期前兩個月的SPEI-1對植被枯黃期具有顯著推遲作用.但隨著時間尺度增加的影響,植被枯黃期對環境干濕敏感性降低,植被枯黃期提前,且空間格局逐漸趨于一致.

4.3 受時間多尺度尺度干濕變化的影響,研究區69.15%的植被枯黃期表現出一定的滯后效應,主要分布于阿巴嘎旗、錫林浩特市、東烏珠穆沁旗等地,響應程度集中在枯黃期前2個月和4個月.植被枯黃期與時間多尺度SPEI在大部分區域均呈不同程度的負相關,總體表現為,滯后4個月的SPEI-3、SPEI-6和SPEI-9對植被枯黃期響應程度最高,且顯著相關面積占比均最大.

[1] R D Koster, G K Walker, G J Collatz, et al. Hydroclimatic controls on the means and variability of vegetation phenology and carbon uptake [J]. Journal of Climate, 2014,27(14):5632-5652.

[2] Irene Garonna, Rogier de Jong, Allard J.W. de Wit,et al . Strong contribution of autumn phenology to changes in satellite-derived growing season length estimates across Europe (1982~2011) [M]. England: Global Change Biology, 2014:3457-3470.

[3] Mcdowell N, Pockman W T, Allen C D, et al. Mechanisms of plant survival and mortality during drought: Why do some plants survive while others succumb to drought? [J]. New Phytologist, 2010,178(4): 719-739.

[4] Begueria-Portugues S, Vicente-Serrano S M, Marta Angulo-Martinez, et al. The standardized precipitation-evapotranspiration index (SPEI): A multiscalar drought index [C]//Zürich: The EMS annual meeting, 2010.

[5] Palmer W C. Meteorological Drought [N]. Research paper, 1965- 2(58).

[6] Mckee T B, Doesken N J, Kleist J. The relationship of drought frequency and duration to time scales[C]. California: Eighth Conference on Applied Climatology, 1993.

[7] Chen F, Yuan Y, Yu S, et al. A 225-year long drought reconstruction for east Xinjiang based on Siberia larch (Larix sibirica) tree-ring widths: Reveals the recent dry trend of the eastern end of Tien Shan [J]. Quaternary International, 2015,358(2):42-47.

[8] Mayank S, Ashish K P, Amalava B, et al. Tree-ring based reconstruction of winter drought since 1767 CE from Uttarkashi, Western Himalaya [J]. Quaternary International, 2017, 479(1):58-69.

[9] Guo M, Wu C, Peng J, et al. Identifying contributions of climatic and atmospheric changes to autumn phenology over mid-high latitudes of Northern Hemisphere [J]. Global and Planetary Change, 2021,197 (6237):103396.

[10] Kang W P, Wang T, Liu S L. The response of vegetation phenology and productivity to drought in semi-arid regions of Northern China [J]. Remote Sensing, 2018,10(5):727.

[11] 黃文琳,張 強,孔冬冬,等.1982~2013年內蒙古地區植被物候對干旱變化的響應 [J]. 生態學報, 2019,39(13):4953-4965.

Huang W L, Zhang Q, Kong D D, et al. Response of vegetation phenology to drought in Inner Mongolia from 1982 to 2013 [J]. Acta Ecologica Sinica, 2019,39(13):4953-4965.

[12] Ding Y, Li Z, Peng S. Global analysis of time-lag and accumulation effects of climate on vegetation growth [J]. International Journal of Applied Earth Observation and Geoinformation, 2020,92:102179.

[13] Hu Z M, Yu G R, Fan J W, et al. Effects of drought on ecosystem carbon and water processes: A review at different scales [J]. Progress in Geography, 2006,25(6):12-20.

[14] Cui T, Lawrence M, Guo X. Grassland phenology response to drought in the Canadian Prairies [J]. Remote Sensing, 2017,9(12):1258.

[15] 胡云鋒.內蒙古錫林郭勒生態系統綜合監測與評估 [M]. 北京:中國環境出版社, 2013:1-33.

Hu Y F. Comprehensive monitoring and evaluation of xilin Guole ecosystem in Inner Mongolia Inner Mongolia [M]. Beijing: China Environmental Science Press, 2013:1-33.

[16] Batu N C, Claas N, Hu Y, et al. Land-use change and land degradation on the Mongolian Plateau from 1975 to 2015—A case study from Xilingol, China [J]. Land Degradation and Development, 2018,29(6): 1595-1606.

[17] Wu D, Zhao X, Zhou T, et al. Diverse responses of global vegetation to climate changes: Spatial patterns and time-lag effects [C]. Washington DC: AGU Fall Meeting Abstracts, 2014:3250-3531.

[18] 杭玉玲,包 剛,包玉海,等.2000~2010年錫林郭勒草原植被覆蓋時空變化格局及其氣候響應[J]. 草地學報, 2014,22(6):1194-1204.

Hang Y L, Bao G, Bao Y H, et al. Spatiotemporal changes of vegetation coverage in Xilin Gol Grassland and its responses to climate change during 2000~2010 [J]. Acta Agrestia Sinica, 2014,22(6):1194-1204.

[19] Tang M F, Wu D, Fu X, et al. An assessment of ecological carrying capacity of Xilingol, Inner Mongolia [J]. The International Journal of Sustainable Development and World Ecology, 2017,24(5):408-414.

[20] Batu N C, Ralf Wieland; Tobia Lakes; et al. Identifying drivers of land degradation in Xilingol, China, between 1975 and 2015 [J]. Land Use Policy, 2019,83(4):543-559.

[21] Piao S L, Liu Z, Wang T, et al. Weakening temperature control on the interannual variations of spring carbon uptake across northern lands [J]. Nature Climate Change, 2017,7(5):359-363.

[22] Bao G, Jin Hgjlt, Tong S Q, et al. Autumn phenology and its covariation with climate, spring phenology and annual peak growth on the Mongolian Plateau [J]. Agricultural and Forest Meteorology, 2021,298-299(8):108-312.

[23] 王曉峰,胡春艷,衛 偉,等.基于SPI的渭北黃土高原干旱時空特征 [J]. 生態環境學報, 2016,25(3):415-421.

Wang Xiaofeng, Hu Chunyan, Wei Wei, et al. Temporal and spatial characteristics of drought based on standardized precipitation index in Weibei Loess Plateau [J]. Ecology and Environmental Sciences, 2016, 25(3):415-421.

[24] 趙 芬.基于CASA模型的錫林郭勒盟草地凈初級生產力遙感估算與驗證 [D]. 北京:中國農業科學院, 2015.

Zhao F. Remote sensing estimation and validation of net primary production in the Xilingol grassland based on CASA model [D]. Beijing: Chinese Academy of Agricultural Sciences Dissertation, 2015.

[25] 包 剛,覃志豪,包玉海,等.1982~2006年蒙古高原植被覆蓋時空變化分析 [J]. 中國沙漠, 2013,33(3):918-927.

Bao G, Qin Z H, Bao Y H, et al. Spatial-temporal changes of vegetation cover in Mongolian Plateau during 1982~2006 [J]. Journal of Desert Research, 2013,33(3):918-927.

[26] Hou X, Shuai G, Zheng N, et al. Extracting grassland vegetation phenology in North China based on cumulative spot-vegetation NDVI data [J]. International Journal of Remote Sensing, 2014,35(9/10): 3316-3330.

[27] Vicente-Serrano S M, Beguería S, López-Moreno J I. A multi-scalar drought index sensitive to global warming: The standardized precipitation evapotranspiration index [J]. Journal of Climate, 2010, 23:1696-1718.

[28] Yu X Y, He X Y, Zheng H F, et al. Spatial and temporal analysis of drought risk during the crop-growing season over northeast China [J]. Natural Hazards, 2014,71(1):275-289.

[29] GB/T 20481-2006 氣象干旱等級[S].

GB/T 20481-2006 Grades of meteorological drought [S].

[30] Vicenteserrano S M, Beguería S, Lópezmoreno J I, et al. A new global 0.5 gridded dataset (1901~2006) of a multiscalar drought index: Comparison with current drought index datasets based on the Palmer drought severity index [J]. Journal of Hydrometeorology, 2010,11(4): 1033 1043.

[31] 薩日蓋,包 剛,包玉海,等.內蒙古植被枯黃期變化及其與氣候和植被生產力的關系 [J]. 應用生態學報, 2020,31(6):1898-1908.

Sa R G, Bao G, Bao Y H, et al. Variation in the vegetation fade stage and its relationships with climate and vegetation productivity in Inner Mongolia, China [J]. Chinese Journal of Applied Ecology, 2020,31(6): 1898-1908.

[32] 劉紅海,王 宏,李曉兵,等.內蒙古溫帶典型草原物候期的變化研究[J]. 北京師范大學學報(自然科學版), 2016,52(4):512-517,406,537.

Liu H H, Wang H, Li X B, et al. Monitoring phenological phase of typical temperate steppe in Inner Mongolia [J]. Journal of Beijing Normal University (Natural Science), 2016,52(4):512-517,406,537.

[33] 烏日汗,紅 雨,包 剛.2001~2016年內蒙古植被物候變化及其對生產力的影響[J]. 草地學報, 2019,27(6):1685-1693.

Wu R H, Hong Y, Bao G. The change of vegetation phenology and its impacts on vegetation [J]. Acta Agrestia Sinica, 2019,27(6):1685- 1693.

[34] 趙汝冰,肖如林,萬華偉,等.錫林郭勒盟草地變化監測及驅動力分析[J]. 中國環境科學, 2017,37(12):4734-4743.

Zhao R B, Xiao R L, Wan W H, et al. Grassland change monitoring and driving force analysis in Xilingol League [J]. China Environmental, 2017,37(12):4734-4743.

[35] 張巧鳳.錫林郭勒草原干旱災害監測與風險評估研究[D]. 北京:中國農業科學院, 2016.

Zhang QF. Study on drought disaster monitoring and risk assessment in Xilingol grassland [D]. Beijing: Chinese Academy of Agricultural Sciences Dissertation, 2016.

[36] Geneva. World meteorological organization [EB/OL]. https://World Meteorological Organisation – Geneva (mfa.ee), 1951,16(1):241-243.

[37] Che M L, Chen B Z, John L, et al. Spatial and temporal variations in the end date of the vegetation growing season throughout the Qinghai- Tibetan Plateau from 1982 to 2011 [J]. Agr Forest Meteorol, 2014, 189-190(5):81-90.

[38] Jackson Robert B, Jobbáagy Esteban G, Avissar Roni, et al. Trading water for carbon with biological carbon sequestration [J]. Science, 2005,310(5756):1944–1947.

[39] Liu Q, Fu Y H, Zeng Z Z, et al. Temperature, precipitation, and insolation effects on autumn vegetation phenology in temperate China. [J]. Global Change Biology, 2016,22(2):644–655.

[40] Yang Y T, Guan H D, Shen M G, et al. Changes in autumn vegetation dormancy onset date and the climate controls across temperate ecosystems in China from 1982 to 2010 [J]. Global Change Biology, 2015,21(2):652-665.

[41] Beguería S, Vicente-Serrano S M, Reig F, et al. Standardized precipitation evapotranspiration index (SPEI) revisited: Parameter fitting, evapotranspiration models, tools, datasets and drought monitoring [J]. International Journal of Climatology, 2014, 34(10):3001-3023.

[42] Ma M, Zhang S W, Wei B C. Temporal and spatial pattern of grassland degradation and its determinants for recent 30Years in xilingol [J]. Chinese Journal of Grassland, 2017,39(4):86-93.

[43] 苗百嶺,梁存柱,王 煒,等.植被退化對典型草原地表徑流的影響[J]. 水土保持學報, 2008,22(2):10-14.

Miao B L, Liang C Z, Wang W, et al. Effects of vegetation on degradation surface runoff of typical steppe [J]. Journal of Soil & Water Conservation, 2008,22(2):10-14.

[44] Hou J, Du LT, Liu K, et al. Characteristics of vegetation activity and its responses to climate change in desert/grassland biome transition zones in the last 30years based on GIMMS3g [J]. Theoretical and Applied Climatology, 2019,136(3/4):915-928.

[45] Elisa D, Wim V E, Maurice D P, et al. Influence of environmental factors light, CO2, temperature, and relative humidity on stomatal opening and development: A review [J]. Agronomy, 2020, 10(12): 1975-1975.

[46] Pillitteri L J, Torii K U. Mechanisms of stomatal development [J]. Annual Review of Plant Biology, 2012,63(1):591– 614.

[47] Lei T J, Feng J, Zheng C Y, et al. Review of drought impacts on carbon cycling in grassland ecosystems [J]. Frontiers of Earth Science, 2020,14(24):462-478.

[48] Xu Z Z, Zhou G S, Shimizu H. Plant responses to drought and rewatering [J]. Plant Signaling & Behavior, 2010,5(6):649–654.

[49] Peng J, Wu C Y, Zhang X Y, et al. Satellite detection of cumulative and lagged effects of drought on autumn leaf senescence over the Northern Hemisphere [J]. Global Change Biology, 2019,25(6):2174-2188.

[50] Wu X C, Liu H Y, Li X Y, et al. Differentiating drought legacy effects on vegetation growth over the temperate Northern Hemisphere [J]. Global Change Biology, 2018,24(1):504-516.

[51] Knapp A, Carroll C J W, Denton E, et al. Differential sensitivity to regional-scale drought in six central US grasslands [J]. Oecologia, 2015,177(4):949-957.

[52] Sillmann J, Roeckner E. Indices for extreme events in projections of anthropogenic climate change [J]. Climatic Change, 2008,86(1):83- 104.

Response of phenological vegetation wilting period to multi-scale drying-wetting changes in Xilingol.

LYU Da1,BAO Gang1,2,3*, TONG Si-qin1,2, LEI Jun1

(1.College of Geographical Science, Inner Mongolia Normal University, Hohhot 010022, China;2.Inner Mongolia Key Laboratory of Remote Sensing & Geography Information System, Inner Mongolia Normal University, Hohhot 010022, China;3.Breeding Base for State Key Laboratory of Land Degradation and Ecological Restoration in Northwest China, Ningxia University, Yinchuan 750021, China)., 2022,42(1):323~335

Based on the third generation《Global Inventory Monitoring and Modelling Studies》(GIMMS) and the《Normalized Difference Vegetation Index》(NDVI) (1982~2015), this study extracted vegetation EOS (End of Growing Season) data in Xilingol grassland from 1982 to 2015, and analyzed its responses to multi-scale drying-wetting changes derived from the Standardized Precipitation Evapotranspiration Index (SPEI). The results showed that during 1982~2015, the EOS of Xilingol vegetation was 265~280 d, and delayed from northeast to southwest, with a decreasing trend of -0.14 d/a. Generally, the interannual variation in spatially averaged SPEI tended to decrease at -0.10d/10a. In addition, an increasing trend was mainly observed in the northeastern meadows, and a decreasing trend was observed in the southwestern desert grasslands. Over the past 34 years, the frequency of drought in Xilingol has increased, with the large areas experiencing the light to mid-droughts. The sudden change to frequent droughts occurred in the 1990s. There were significant spatial differences in the impacts of multiple SPEI time scales on EOS. In detail, the positive effect of monthly SPEI-1(Index for one month data) on vegetation EOS in the Xilingol was the greatest, suggesting that grassland vegetation EOS could be delayed in a favorable dry-wet environment over a short time scale. The SPEI-3, SPEI-6, SPEI-9, and SPEI-12(Index for many months data), however, had relatively large negative correlations with EOS, indicating that the long-term changes in the drying conditions could advance grassland EOS in Xilingol. The time-lag effect of SPEI-1 in the Xilingol mainly occurs two months before EOS. The long-term SPEI-3, SPEI-6, and SPEI-9 (especially SPEI-3), which are accumulated to June, have particularly considerable lag effects.

Xilingol;end of season;standardized precipitation evapotranspiration index;climate response

X87

A

1000-6923(2022)01-0323-13

呂 達(1995-),男,內蒙古呼和浩特人,內蒙古師范大學碩士研究生,主要研究方向為資源與環境遙感.發表論文1篇.

2021-05-12

國家自然科學基金資助項目(42061070);內蒙古自治區自然科學基金資助項目(2021MS04014);內蒙古高??萍加⒉彭椖?NJYT- 18-A11)

* 責任作者, 副研究員, baogang@imnu.edu.cn

猜你喜歡
研究
FMS與YBT相關性的實證研究
2020年國內翻譯研究述評
遼代千人邑研究述論
視錯覺在平面設計中的應用與研究
科技傳播(2019年22期)2020-01-14 03:06:54
關于遼朝“一國兩制”研究的回顧與思考
EMA伺服控制系統研究
基于聲、光、磁、觸摸多功能控制的研究
電子制作(2018年11期)2018-08-04 03:26:04
新版C-NCAP側面碰撞假人損傷研究
關于反傾銷會計研究的思考
焊接膜層脫落的攻關研究
電子制作(2017年23期)2017-02-02 07:17:19
主站蜘蛛池模板: 亚洲一区二区视频在线观看| 日本精品中文字幕在线不卡| 国产a在视频线精品视频下载| 激情五月婷婷综合网| 天堂亚洲网| 欧美一级特黄aaaaaa在线看片| 特级毛片免费视频| 婷婷久久综合九色综合88| 黄色一级视频欧美| 毛片视频网址| 在线日韩日本国产亚洲| 亚洲三级成人| 精品国产电影久久九九| 国产精品一区二区久久精品无码| 日本免费福利视频| 欧美日韩国产在线人成app| 亚洲 欧美 中文 AⅤ在线视频| 一区二区日韩国产精久久| 免费在线看黄网址| 国产18页| 欧美第二区| 中国黄色一级视频| 精品少妇三级亚洲| 国模在线视频一区二区三区| 国产日产欧美精品| 日本午夜影院| 久久精品亚洲专区| 亚洲aⅴ天堂| 亚洲码在线中文在线观看| 亚洲VA中文字幕| 亚洲欧洲日产国产无码AV| 黄片在线永久| 亚洲高清国产拍精品26u| 老司国产精品视频91| 粉嫩国产白浆在线观看| 美女被躁出白浆视频播放| 国产肉感大码AV无码| 国产综合欧美| 91在线丝袜| 日本久久网站| 亚洲日本中文字幕乱码中文 | 久久大香伊蕉在人线观看热2| 国产人碰人摸人爱免费视频| 91www在线观看| 国产91视频免费观看| 欧美一级黄色影院| 国语少妇高潮| 狠狠色丁香婷婷综合| 久久男人资源站| 亚洲欧美自拍视频| 亚洲综合18p| 亚洲女同欧美在线| 好吊色妇女免费视频免费| 波多野结衣亚洲一区| 免费一极毛片| 四虎精品国产AV二区| 亚洲h视频在线| 波多野结衣在线se| 波多野结衣一级毛片| 另类欧美日韩| 午夜国产在线观看| 色首页AV在线| 中文字幕不卡免费高清视频| 日本午夜三级| 亚洲第一视频网| 26uuu国产精品视频| 波多野结衣一二三| 8090成人午夜精品| 欧美a在线看| 亚洲日韩在线满18点击进入| 国产精品亚洲欧美日韩久久| 欧美国产在线一区| 久久综合丝袜长腿丝袜| 97超碰精品成人国产| igao国产精品| 五月激情婷婷综合| 国产精品视频久| 久久无码av三级| 亚洲第一天堂无码专区| 亚洲欧美一区二区三区麻豆| 国产在线视频福利资源站| 国产在线拍偷自揄观看视频网站|