1982—2020年黃河流域植被變化特征及驅(qū)動(dòng)因素

王紫荊，徐夢(mèng)珍，胡宏昌，張向萍

(1. 清華大學(xué)水利水電工程系，北京 100084；2. 黃河水利科學(xué)研究院，河南 鄭州 450003)

黃河流域生態(tài)環(huán)境脆弱，資源環(huán)境高負(fù)載，是中國(guó)人地矛盾最為突出的區(qū)域之一。對(duì)黃河流域人地關(guān)系空間分異及演化規(guī)律的認(rèn)識(shí)，是制定流域可持續(xù)發(fā)展策略的基礎(chǔ)。然而，黃河流域是長(zhǎng)期受強(qiáng)人類活動(dòng)干擾的復(fù)雜巨系統(tǒng)，流域各地差異極大，深入理解其物質(zhì)能量傳遞與轉(zhuǎn)換過(guò)程頗具挑戰(zhàn)[1-2]。植被在陸地生態(tài)系統(tǒng)的碳、水以及能量循環(huán)中起關(guān)鍵性作用[3-4]，其變化被認(rèn)為是環(huán)境變化的指紋[5]，植被變化及其對(duì)流域的綜合影響是當(dāng)前黃河流域研究關(guān)注的關(guān)鍵問(wèn)題[6-7]。黃河流經(jīng)的黃土高原是中國(guó)植被綠度增加趨勢(shì)最突出的地區(qū)[8-9]，生態(tài)工程帶來(lái)的植被改善有效控制侵蝕，增加碳匯，改善了當(dāng)?shù)丨h(huán)境，是黃河進(jìn)入21世紀(jì)以來(lái)水沙驟減的重要原因[10-11]。然而，植被改善也導(dǎo)致蒸散發(fā)的增加，水資源消耗顯著增加，甚至有研究認(rèn)為黃土高原上恢復(fù)的植被水平已接近水資源可持續(xù)利用可承受的極限[12]。針對(duì)1982—2016年黃淮海流域植被指數(shù)年序列與徑流關(guān)系的研究表明，黃河流域植被改善對(duì)徑流的影響有顯著的空間差異[6]。除生態(tài)工程外，農(nóng)業(yè)集約化發(fā)展同樣會(huì)帶來(lái)顯著的植被指數(shù)增長(zhǎng)[9]。黃河流域是中國(guó)重要的糧食產(chǎn)區(qū)，其農(nóng)業(yè)活動(dòng)對(duì)流域植被變化的影響尤為突出[13]。在嚴(yán)峻的水資源短缺情勢(shì)下，流域生態(tài)工程與糧食生產(chǎn)間存在用水沖突[7，14]。然而，已有研究多將植被綠度的增加歸于生態(tài)工程，忽略了黃河流域農(nóng)業(yè)集約化的影響。綜上，在黃河流域植被變化分析中，區(qū)分生態(tài)恢復(fù)與農(nóng)業(yè)生產(chǎn)活動(dòng)的影響，可提升對(duì)黃河流域復(fù)雜巨系統(tǒng)的認(rèn)識(shí)，為流域生態(tài)保護(hù)與高質(zhì)量發(fā)展以及山水林田湖草沙一體化治理目標(biāo)提供支持。

已有研究多使用基于遙感的植被指數(shù)年序列趨勢(shì)分析，來(lái)指示黃河流域生態(tài)工程下的植被變化[7]。植被物候是植物長(zhǎng)期適應(yīng)環(huán)境的季節(jié)性變化而形成的生長(zhǎng)發(fā)育節(jié)律[15]。遙感物候?qū)W挖掘氣候變化的生物指示，在平滑植被指標(biāo)年內(nèi)序列后，用基于導(dǎo)數(shù)、閾值或是模型擬合方法提取植被序列年內(nèi)特征[16-17]。提取的年內(nèi)指標(biāo)不僅受氣候變化影響，也會(huì)直接反饋人類活動(dòng)下植被覆蓋、物種組成的變化[17]。因此，遙感物候分析方法可提供挖掘黃河流域植被信息、細(xì)化植被變化特征的新思路。在大尺度、強(qiáng)人類干擾的黃河流域植被分析中，使用逐日植被遙感數(shù)據(jù)，應(yīng)用遙感物候提取方法，有望獲得更為精細(xì)的流域植被時(shí)空變化特征。

本文在植被年際變化分析的基礎(chǔ)上，選擇Savitzky-Golay濾波與雙logistic模型結(jié)合的方法提取植被年內(nèi)變化特征，以標(biāo)準(zhǔn)化的方法挖掘黃河流域各區(qū)域植被信息；通過(guò)機(jī)器學(xué)習(xí)聚類方法，歸納1982—2020年黃河流域植被變化特征，探究其驅(qū)動(dòng)因素；基于植被變化與黃河水資源利用情況，探討黃河流域生態(tài)保護(hù)與高質(zhì)量發(fā)展對(duì)策。

1 研究區(qū)概況與數(shù)據(jù)來(lái)源

1.1 研究區(qū)概況

黃河發(fā)源于青藏高原，流經(jīng)9省(自治區(qū))，在山東省注入渤海，流域地勢(shì)西高東低，流域面積約79.5萬(wàn)km2。流域氣候?qū)儆诖箨懶詺夂?，東南部半濕潤(rùn)，中部半干旱，西北部干旱，各地氣候差異大、地貌類型多樣、生境復(fù)雜，孕育了豐富的生態(tài)系統(tǒng)。本研究將黃河流域外的下游引黃灌區(qū)納入研究區(qū)，并根據(jù)黃河干流7個(gè)關(guān)鍵水文站點(diǎn)將黃河流域分為7個(gè)區(qū)域(圖1)。黃河源區(qū)(唐乃亥以上，S1)平均海拔在4 000 m以上，主要植被類型為高寒草甸；中部(唐乃亥—花園口，S2—S6)地處內(nèi)蒙古高原和黃土高原，海拔在1 000～2 000 m，水土流失嚴(yán)重，植被自西北向東南大致可分為4個(gè)植被帶，依次為中溫帶荒漠草原、干草原、森林草原及暖溫帶夏綠闊葉林；東部(花園口—利津，S7)海拔不超過(guò)100 m，由黃河下游沖積平原和魯中丘陵組成，主要植被類型為栽培植被。

圖1 研究區(qū)分區(qū)、植被類型及土地利用類型Fig.1 Map of division，vegetation types，and land use types of the study area

1.2 數(shù)據(jù)來(lái)源

本文選擇1982—2020年NOAA氣候數(shù)據(jù)記錄(Climate Data Record，CDR)歸一化植被指數(shù)(Normalized Difference Vegetation Index，NDVI)數(shù)據(jù)產(chǎn)品分析黃河流域植被變化特征[18]。該產(chǎn)品為逐日數(shù)據(jù)，空間分辨率為0.05°×0.05°，時(shí)間跨度長(zhǎng)，且相比于黃河流域相關(guān)研究多使用的GIMMS NDVI 3g時(shí)空分辨率均大幅提升。此外，本文補(bǔ)充了植被類型、耕地類型、土地利用信息以分析植被變化驅(qū)動(dòng)因素。植被類型數(shù)據(jù)選擇國(guó)家自然資源和自然條件的基本圖件《1∶1 000 000中國(guó)植被圖集》，土地利用數(shù)據(jù)選擇2020年1 km分辨率中國(guó)多時(shí)期土地利用遙感監(jiān)測(cè)數(shù)據(jù)集(CNLUCC)[19]，2項(xiàng)數(shù)據(jù)可從中國(guó)科學(xué)院資源環(huán)境與數(shù)據(jù)中心獲得。耕地類別來(lái)源于美國(guó)國(guó)家航空航天局(NASA)提供的全球1 km糧食安全支持分析數(shù)據(jù)GFSAD1KCM，該產(chǎn)品基于2017—2012年信息生成[20]。

2 研究方法

2.1 年內(nèi)特征提取方法

本文利用Savitzky-Golay濾波重構(gòu)遙感數(shù)據(jù)[21-22]，利用雙Logistic函數(shù)擬合獲得年內(nèi)特征值。Savitzky-Golay濾波一般方程如下：

(1)

將重構(gòu)的NDVI序列擬合至雙Logistic函數(shù)來(lái)獲取年內(nèi)序列特征，擬合方程如下[23-24]：

(2)

式中：f(t)為特定日期t的INDV；v1為擬合的全年INDV的最小值；v2為擬合的全年INDV最大值與最小值之差；參數(shù)m和n用于確定NDVI增大階段和減小階段的總體斜率和基本相位，m1、n1和m2、n2分別是一對(duì)參數(shù)。該模型6個(gè)參數(shù)通過(guò)Levenberg-Marquardt算法求解。擬合模型曲線斜率最大值對(duì)應(yīng)的時(shí)間為生長(zhǎng)期開(kāi)始時(shí)間(ISOS)，斜率最小值對(duì)應(yīng)的時(shí)間為生長(zhǎng)期結(jié)束時(shí)間(IEOS)，兩者之差為生長(zhǎng)期持續(xù)時(shí)間(IGSL)。該方法不依賴主觀經(jīng)驗(yàn)，具有自適應(yīng)性，適用于多種植被類型，可應(yīng)用于黃河流域(圖2)。

圖2 植被年內(nèi)特征提取方法示意Fig.2 Examples of vegetation intra-annual indices extraction

2.2 植被變化分析及特征歸納方法

3 結(jié)果與討論

3.1 黃河流域植被指標(biāo)提取結(jié)果

本文對(duì)研究區(qū)1982—2020年間逐年逐像元的植被指標(biāo)進(jìn)行了提取，共計(jì)42 403個(gè)像元，其中有效識(shí)別41 570個(gè)，識(shí)別率達(dá)98%，表明本文選擇的年內(nèi)指標(biāo)提取方法可靠有效。圖3呈現(xiàn)了2個(gè)時(shí)期(1980s、2010s)的植被指標(biāo)提取結(jié)果，選擇各像元序列1982—1990年的中位數(shù)代表1980s的水平，2011—2020年的中位數(shù)代表2010s的水平。圖3中，4個(gè)指標(biāo)呈現(xiàn)了直觀的時(shí)空分異，體現(xiàn)了所選指標(biāo)對(duì)植被時(shí)空變化特征的表征能力?？臻g上大致趨勢(shì)為由東南至西北INDV，max減小，ISOS推遲，IEOS提前，IGSL縮短。時(shí)間變化上，從1980s到2010s，INDV，max普遍增長(zhǎng)，以中游S4和S5(頭道拐—三門峽)尤為突出。2個(gè)階段年內(nèi)特征ISOS、IEOS及IGSL同樣有較大的變化，但流域內(nèi)變化趨勢(shì)不統(tǒng)一，如中上游S3、S4(蘭州—龍門)大范圍ISOS推遲，IGSL縮短；而在下游S7(花園口以下)趨勢(shì)則完全相反，ISOS提前，IGSL延長(zhǎng)。

表1 植被聚類特征Table 1 Characteristics of vegetation clusters

圖3 黃河流域植被指標(biāo)提取結(jié)果Fig.3 Results of vegetation indices extraction

3.2 黃河流域植被變化趨勢(shì)分析

3.2.1分區(qū)趨勢(shì)分析

圖4 分區(qū)INDV，max統(tǒng)計(jì)及年變化序列Fig.4 Annual regional statistics and their variation of maximum NDVI(INDV，max)

圖5 分區(qū)植被年內(nèi)指標(biāo)統(tǒng)計(jì)及年變化序列Fig.5 Annual regional statistics and their variation of phenological indices

3.2.2逐像元趨勢(shì)分析

4項(xiàng)植被指標(biāo)逐像元計(jì)算的Kendall秩相關(guān)系數(shù)分布如圖6所示。在4項(xiàng)指標(biāo)中，INDV，max的變化總體最強(qiáng)，整個(gè)黃河流域表現(xiàn)為普遍的增加趨勢(shì)，這與已有的黃河流域植被研究結(jié)果一致[8]。細(xì)節(jié)上，81%的像元表現(xiàn)了顯著增加趨勢(shì)，僅1%的像元顯著減小(圖6(a))。7個(gè)分區(qū)中，生態(tài)工程的主要實(shí)施區(qū)S4和S5(頭道拐—三門峽)西北部植被改善顯著，S4中INDV，max檢出顯著增加的像元達(dá)到了97%。除了該地區(qū)外，可以觀察到黃河上游河套灌區(qū)、中游汾渭灌區(qū)表現(xiàn)了不弱于黃土高原中部的INDV，max增長(zhǎng)趨勢(shì)，這些歷史悠久的黃河灌區(qū)與生態(tài)工程關(guān)聯(lián)較小，植被變化主要由農(nóng)業(yè)活動(dòng)驅(qū)動(dòng)。然而，現(xiàn)有研究多忽略了上述農(nóng)業(yè)因素，僅考慮單一的生態(tài)工程對(duì)植被恢復(fù)的影響[6-7]。

圖6 1982—2020年黃河流域植被指標(biāo)時(shí)空變化趨勢(shì)Fig.6 Temporal and spatial trend of vegetation indices from 1982 to 2020 in the Yellow River basin

3.3 黃河流域植被變化特征及驅(qū)動(dòng)因素

為了進(jìn)一步探究黃河流域植被1982—2020年間變化的特征和時(shí)空異質(zhì)性，以及其驅(qū)動(dòng)因素，在基于上述植被趨勢(shì)分析結(jié)果的基礎(chǔ)上，補(bǔ)充植被類別和耕地信息，用k均值聚類方法歸納植被變化特征，分析不同植被變化的驅(qū)動(dòng)因素，揭示人類活動(dòng)對(duì)黃河流域的改變。

聚類結(jié)果如圖7所示。通過(guò)列舉，k=7為聚類結(jié)果最優(yōu)。相比于k=6，增加的聚類像元數(shù)目與其他類別相當(dāng)，且主要分布于黃土高原中部過(guò)渡帶，是生態(tài)工程的重點(diǎn)區(qū)域；而將k增至8時(shí)，增加的聚類樣本數(shù)量小，且較為分散。為了便于區(qū)分，綜合主要分布區(qū)特征與土地利用為聚類結(jié)果命名。圖7中7個(gè)餅圖為分區(qū)內(nèi)各聚類的占比，表1列出了各聚類的植被變化特征、主要植被類型與土地利用類型。中國(guó)植被的顯著改變主要由人類活動(dòng)主導(dǎo)[9]，在黃河流域具體為生態(tài)工程與農(nóng)業(yè)集約化，因此，本文得到的7個(gè)聚類可再簡(jiǎn)單歸為兩大組，分別由生態(tài)工程與農(nóng)業(yè)活動(dòng)主導(dǎo)。

圖7 黃河流域植被變化特征聚類結(jié)果Fig.7 Results of vegetation change clustering in the Yellow River basin

本文將分布于黃河流域重點(diǎn)農(nóng)業(yè)生產(chǎn)區(qū)的2個(gè)聚類命名為中上游耕地類和下游耕地類。中上游耕地類INDV，max增長(zhǎng)顯著，下游耕地類ISOS顯著提前。中上游耕地類INDV，max增加顯著主要是由于灌溉、肥料等農(nóng)業(yè)技術(shù)發(fā)展下干旱半干旱區(qū)種植能力提升。圖3(a)中，1980s下游耕地INDV，max明顯高于中上游耕地，而圖3(b)中，2010s重點(diǎn)耕地區(qū)域INDV，max水平無(wú)明顯差異。年內(nèi)指標(biāo)變化差異主要是農(nóng)業(yè)集約化下黃河主要灌區(qū)確定的種植結(jié)構(gòu)與熟制差異。耕地類在汾渭灌區(qū)(位于S5)的分界線與春玉米和夏玉米種植帶的分界線接近[27]。下游耕地(S7)冬小麥和夏玉米的規(guī)?；沟肐NDV曲線轉(zhuǎn)為顯著的雙峰型(圖2(b))，ISOS提前，IGSL延長(zhǎng)。中上游耕地類氣候條件不能支撐一年兩熟，識(shí)別的IGSL較下游短。特別地，1980年以來(lái)河套灌區(qū)IGSL縮短顯著，這主要是由于農(nóng)業(yè)節(jié)水下種植結(jié)構(gòu)的改變。作為節(jié)水的重點(diǎn)區(qū)域，河套灌區(qū)耗水量大的小麥種植面積極速萎縮，轉(zhuǎn)為葵花和玉米，在2018年河套灌區(qū)小麥播種面積僅占9.29%[28-29]。替代小麥的葵花和玉米生長(zhǎng)期比小麥短，且生長(zhǎng)高峰期晚于小麥[30]，因而河套灌區(qū)表現(xiàn)為ISOS推遲，IGSL縮短。

3.4 啟示與建議

水資源短缺是黃河流域的最大矛盾[31]。黃河流域植被變化及其對(duì)水資源的影響是黃河流域生態(tài)保護(hù)與高質(zhì)量發(fā)展的關(guān)鍵科學(xué)問(wèn)題。本文可為解釋黃河流域植被變化對(duì)徑流的作用有顯著的空間差異性提供支持。研究結(jié)果表明，1982—2020年農(nóng)業(yè)集約化對(duì)黃河流域植被的影響不弱于生態(tài)工程。農(nóng)業(yè)是黃河水資源使用的主要行業(yè)，根據(jù)2020年《黃河水資源公報(bào)》，2020年黃河地表水農(nóng)業(yè)耗水量為231.01億m3，占總耗水量的65.3%。在山水林田湖草沙一體化治理的要求下，黃河流域的“變綠”不能僅強(qiáng)調(diào)生態(tài)工程，在關(guān)注“林”與“草”的同時(shí)，也要重視“田”這一關(guān)鍵環(huán)境要素。經(jīng)典的INDV年序列分析不能區(qū)分農(nóng)業(yè)集約化與生態(tài)工程。本文使用遙感物候?qū)W方法和長(zhǎng)序列植被逐日數(shù)據(jù)提取的植被年內(nèi)指標(biāo)對(duì)植被組成變化敏感，是解析黃河流域復(fù)雜植被變化的有效途徑。

水資源短缺背景下，黃河流域汾渭平原、河套平原和黃淮海平原是《國(guó)家主體功能區(qū)規(guī)劃》劃定的全國(guó)重要糧食主產(chǎn)區(qū)，是保障中國(guó)糧食安全的關(guān)鍵地區(qū)[13]。黃河上游河套灌區(qū)通常被認(rèn)為是黃河流域節(jié)水潛力所在，但根據(jù)本文結(jié)果，河套灌區(qū)提取的生長(zhǎng)期顯著縮短，表明種植結(jié)構(gòu)發(fā)生了巨大轉(zhuǎn)變，高耗水的小麥種植比例大幅下降，這有效地減少了農(nóng)業(yè)用水。然而與此同時(shí)，農(nóng)業(yè)節(jié)水下河套灌區(qū)平均地下水埋深下降了0.5 m以上，部分地區(qū)出現(xiàn)地下水漏斗，對(duì)自然植被的水分支撐作用減弱[1]。引黃灌溉是維持干旱區(qū)綠洲生態(tài)健康的重要水源，且隨著未來(lái)荒漠化治理的繼續(xù)推進(jìn)，干旱半干旱區(qū)生態(tài)工程的耗水量將持續(xù)增加，上游灌區(qū)進(jìn)一步節(jié)水可能威脅當(dāng)?shù)厣鷳B(tài)。因此，對(duì)于上游的干旱灌溉綠洲區(qū)，需要加強(qiáng)對(duì)農(nóng)業(yè)節(jié)水生態(tài)影響的動(dòng)態(tài)監(jiān)控，以防過(guò)度節(jié)水破壞生態(tài)環(huán)境。

與此同時(shí)，黃河中下游耕地提取的生長(zhǎng)期顯著延長(zhǎng)，冬小麥夏玉米的輪種模式規(guī)?；１M管農(nóng)業(yè)節(jié)水技術(shù)在不斷進(jìn)步，但在高種植強(qiáng)度下很難大幅縮減該地區(qū)農(nóng)業(yè)用水。另一方面，黃河流域生態(tài)工程深刻地改變了流域植被條件，生態(tài)工程對(duì)流域水循環(huán)的作用有待進(jìn)一步探究。在氣候變化背景下，未來(lái)氣候不確定性增加。2022年中國(guó)就出現(xiàn)了嚴(yán)重的旱情，糧食安全受到威脅[32]。綜上，黃河水-糧食-生態(tài)沖突下，須加快推進(jìn)流域外調(diào)水以緩解用水壓力[31-33]，重視農(nóng)業(yè)用水保障，優(yōu)化用水分配，預(yù)測(cè)未來(lái)生態(tài)用水的變化，核算用水效益，在用水沖突中做好權(quán)衡，確保水、糧食與生態(tài)安全[34]。

4 結(jié) 論

本文選擇1982—2020年NOAA CDR NDVI逐日數(shù)據(jù)集，在傳統(tǒng)植被NDVI年序列分析基礎(chǔ)上補(bǔ)充使用遙感物候?qū)W方法提取的植被日序列年內(nèi)特征，用k均值聚類方法歸納植被變化特征，分析植被變化驅(qū)動(dòng)因素，最后給出基于植被的黃河流域管理建議。

1982—2020年黃河流域植被變化廣泛且顯著。黃河流域生態(tài)工程驅(qū)動(dòng)下各地植被大幅改善，NDVI年最大值均顯著提升，但植被年內(nèi)特征變化不同，濕潤(rùn)半濕潤(rùn)區(qū)提取的生長(zhǎng)期持續(xù)時(shí)間延長(zhǎng)，干旱半干旱區(qū)提取的生長(zhǎng)期持續(xù)時(shí)間縮短，不同氣候條件下的生態(tài)工程影響不能一概而論。黃河流域農(nóng)業(yè)集約化下，干旱半干旱的中上游耕地種植水平大幅提升，NDVI年最大值增加，其中上游河套灌區(qū)種植結(jié)構(gòu)改變，高耗水的小麥轉(zhuǎn)為葵花和玉米，提取的生長(zhǎng)期持續(xù)時(shí)間顯著縮短；水熱條件較好的下游耕地冬小麥夏玉米的輪種模式規(guī)?；C(jī)械化，提取的生長(zhǎng)期持續(xù)時(shí)間顯著延長(zhǎng)。黃河流域耕地呈現(xiàn)的植被變化顯著性與黃土高原中部生態(tài)工程重點(diǎn)實(shí)施區(qū)相當(dāng)，農(nóng)業(yè)集約化是黃河流域植被研究中不可忽略的重要影響因素。黃河流域生態(tài)工程與農(nóng)業(yè)生產(chǎn)存在突出的用水矛盾，需重視農(nóng)業(yè)用水優(yōu)化，核算生態(tài)用水效益，并加快推進(jìn)流域外調(diào)水以緩解用水壓力，確保水-糧食-生態(tài)安全。