2000—2020年環渤海地區植被NDVI時空演變及其與人口密度的相關性

徐 勇, 趙 純, 竇世卿, 郝文強, 鄭志威, 靖娟利

(桂林理工大學 測繪地理信息學院, 廣西 桂林 541006)

1 研究區概況

環渤海地區以“C”字形環抱渤海，包括三省兩市，即遼寧省、山東省、河北省、天津市以及北京市，位于34°22′—43°28′N，113°26′—125°52′E。氣候類型為溫帶季風氣候，主要植被帶為暖溫帶落葉闊葉林帶。地形主要以平原、丘陵和山地為主，整體上地勢呈東南低西北高的分布格局。環渤海地區是中國最為發達的經濟區之一，基礎設施完善，自然資源豐沛，工業基礎扎實，科技實力雄厚。2019年環渤海地區年末總人口約2.6億(2020年統計年鑒)，約占全國總人口的18.38%，為中國人口最為稠密的地區之一，尤其是北京市和天津市，不僅是中國政治、經濟和文化中心，2020年人口密度分別高達1 334和1 158人/km2，分別居全國第2和第3位，而山東省和河北省為中國人口大省，人口總量分居全國第2和第6(第7次人口普查數據)。

2 數據來源與研究方法

2.1 數據來源

本文所使用 MODIS NDVI數據來源于NASA MODIS陸地產品組發布的植被指數產品MOD13A3，該產品時間分辨率為每月，空間分辨率為1 km，時間跨度從2000年2月至2020年12月，該數據已經經過大氣校正，去除了水、云、重氣溶膠和云影的影響。本文通過對得到的每期MOD13A3數據集進行影像鑲嵌、重采樣、投影轉換、數據格式轉換、最大值合成等預處理，生成2000—2020年覆蓋中國全境的年最大值植被NDVI時間序列，然后通過裁剪得到環渤海地區植被年最大合成NDVI時間序列。2000年1月植被NDVI數據缺失。由于1月為冬季，植被NDVI較低，故本文使用最大值合成法得到年最大合成植被NDVI，缺失月份不會對本文研究結果產生影響。人口密度數據來源于WorldPop發布的Population Density數據，時間分辨率為每年，空間分辨率為1 km，時間跨度為2000—2020年。人口密度數據需要經過鑲嵌、投影轉換、圖像裁剪和重采樣等預處理，裁剪得到2000—2020年環渤海地區人口密度時間序列。

2.2 研究方法

2.2.1 Theil-Sen Median趨勢分析與Mann-Kendall顯著性檢驗 本文采用Theil-Sen Median趨勢分析和Mann-Kendall顯著性檢驗[15]，分析環渤海地區2000—2010年、2010—2020年和2000—2020年植被NDVI時空演變特征及其變化的顯著性。Theil-Sen Median趨勢分析是一種非參數統計的趨勢分析方法，可以有效去除離散值和異常值對趨勢結果的影響，常被用于探究長時間序列數據的變化趨勢。

Theil-Sen median趨勢分析計算公式為：

(1)

式中：NDVIi，NDVIj分別表示第i和第j年的植被NDVI值；βNDVI代表植被NDVI的變化斜率,當βNDVI>0時,表示在研究時間段內植被NDVI呈上升趨勢,植被覆蓋有所增加; 當βNDVI=0時,表示在研究時間段內植被NDVI基本保持不變,植被覆蓋較為穩定;當βNDVI<0時,表示在研究時間段內植被NDVI呈下降趨勢,植被覆蓋有所減少。

Mann-Kendall顯著性檢驗是穩健的非參數統計方法，不需要方差正態性假設，在一定程度上避免異常值對分析結果的影響。本文利用Mann-Kendall顯著性檢驗來判斷植被NDVI時間序列變化趨勢的顯著性。具體公式為：

給定NDVI時間序列：

{NDVIi},i=2000，2001，2002，…,2020

(2)

計算檢驗統計量S：

(3)

符號函數sgn：

(4)

計算方差：

(5)

定義標準化統計量Z：

(6)

按照Mann-Kendall顯著性檢驗的檢驗統計量Z值，將植被NDVI變化趨勢分成4個等級：顯著下降、不顯著下降、不顯著上升和顯著上升。當|z|>1.96即置信水平為95%[16]下，為顯著變化，當0<|z|<1.96，為不顯著變化。

(7)

(8)

2.2.3 格蘭杰因果分析 為檢驗X組時間序列是否為Y組時間序列的原因，采用格蘭杰因果關系檢驗法驗證其因果關系[18]。本文設Y為年際植被NDVI即Y={yt}，X為年際人口密度即X={xt}，假設兩者之間是否具有因果解釋能力，并設定滯后階，進行格蘭杰檢驗分析人口密度能否影響植被NDVI。其回歸模型為：

(9)

式中：βt為時間序列X的趨勢項；γt為時間序列Y的趨勢項;n為滯后階數;α為常數項;εt為殘差項。

假定原假設H0：人口密度不是引起植被NDVI的格蘭杰原因(H0:?t∈(1,2,…,n),γt=0)；其余假設H1：人口密度是引起植被NDVI的格蘭杰原因〔H1:?t∈(1,2,…,n),γt≠0〕。當不考慮X的情況下(即假設βt=0)，設置回歸模型的殘差平方和RSR;當考慮X即人口密度的情況下(即假設βt≠0)，設置殘差平方和RSU。構建F統計值：

(10)

選取置信水平為95%，若F≥F0.05(n,N-2n-1)，即拒絕原假設H0，認為人口密度是引起植被NDVI的格蘭杰原因，反之則保持原假設。

式中：F為統計值； RSR為沒有X約束下的殘差平方和； RSU為在X約束下的殘差平方和；n為待估參數的個數；N為樣本容量。

2.2.4 重心轉移模型和標準差橢圓 為分析2000—2020年環渤海地區植被NDVI和人口密度遷移特征，本次研究運用屬性重心遷移模型[19]進行植被NDVI和人口密度演化模式的定量化表達，具體公式為:

(11)

(12)

標準差橢圓從數據范圍、數據集聚性、變化方向和形狀等多個角度，呈現地理要素的空間格局分布特征。方向角是正北方向與順時針旋轉的橢圓長軸形成的夾角，代表數據空間分布的主趨勢方向，長半軸表示數據在主趨勢方向上集聚性，短半軸反映數據在此趨勢方向上集聚性。主要計算公式為：

tanθ=

(13)

(14)

(15)

3 結果與分析

3.1 植被NDVI的時空演變特征

3.1.1 NDVI的時間變化趨勢 由圖1可知，整體上，2000—2020年環渤海地區植被NDVI呈波動上升趨勢，上升速率為0.022/10 a。2000—2020年環渤海地區植被NDVI均值為0.728，最小值出現在2000年，為0.666，最大值出現在2011年，為0.748，且2009和2014年植被NDVI值明顯低于相鄰年份。2000—2010年環渤海地區植被NDVI呈上升趨勢，上升速率為0.057/10 a， 2000—2010年環渤海地區植被NDVI均值為0.716，略低于2000—2020年環渤海地區植被NDVI平均值。相較于2000—2010年和2010—2020年環渤海地區植被NDVI上升趨勢較為平緩，上升速率僅為0.000 1/10 a。2010—2020年環渤海地區平均植被NDVI值最高，為0.738。綜上所述，整體上，2000—2020年環渤海地區植被覆蓋呈改善態勢，改善趨勢明顯，且2000—2010年植被NDVI改善程度強于2010—2020年，這與徐勇等[15]和張鵬騫等[20]的研究結果一致。自2000年以來，針對環渤海地區采取了林草地植被保護、小流域生態治理和高原地區水土保持等一系列林業生態工程，且取得了良好的效果，加之為籌備2008年奧運會，北京市開展了一系列綠化工程，短期內明顯提升了區域植被NDVI值。2009年植被NDVI值明顯低于2008和2010年，這與李卓等[21]的研究結果一致，是由于2009—2010年環渤海地區遭受了歷史極嚴重的旱災，導致整體植被NDVI值較低。

圖1 環渤海地區2000—2020年植被NDVI變化趨勢

3.1.2 植被NDVI的空間變化特征 本文利用Theil-Sen Median趨勢分析法得到2000—2010年、2010—2020年以及2000—2020年環渤海地區植被NDVI變化趨勢，并根據Mann-Kendall顯著性檢驗，檢驗其變化趨勢的顯著性。由表1和圖2可知，2000—2010年和2010—2020年環渤海地區植被NDVI變化呈明顯空間異質性，變化斜率分別在-0.076～0.072/a和-0.087～0.070/a之間。2000—2010年植被NDVI呈上升趨勢的區域面積占83.83%，其中，呈顯著上升的面積占30.26%，均勻分布在研究區內，無明顯集現象；呈顯著下降的面積僅占1.35%，零星分布在各大中城市中心及其周邊地區。而2010—2020年植被NDVI呈上升趨勢的面積為50%，其中，呈顯著上升的面積僅占9.85%，主要分布在北京和河北省西北部以及遼寧省東部；呈顯著下降的面積占5.71%，零星分布在遼寧省中部、天津以及河北和山東部分地區。由上可知，2000—2010年環渤海地區植被NDVI以上升為主，上升區域面積明顯高于下降區域面積，而2010—2020年環渤海地區植被NDVI上升趨勢明顯減緩，植被NDVI呈上升趨勢面積等于呈下降趨勢面積。整體上，2000—2020年植被NDVI變化斜率在-0.034～0.037/a之間，植被NDVI整體呈上升趨勢。植被NDVI上升區域面積遠大于下降區域面積。植被NDVI呈上升趨勢的面積占76.37%，下降趨勢面積占23.68%。植被NDVI顯著上升區域面積占總面積的45.18%，主要分布在遼寧省東部和西部、北京和河北西北部、山東省中部以及河北省和山東省交界處。這與Zheng等[5]和Jiang等[22]和研究結果一致，以上地區植被NDVI上升，主要得益于環渤海地區多年來實施退耕還林，京津風沙源治理工程和三北防護林工程等項目，使得地表植被覆蓋程度得到提升，植被生態穩定性以及抵御風險的能力得到加強[6,8-9]。植被NDVI呈顯著下降趨勢區域占總面積的6.42%，主要分布在遼寧省中部、北京市城市中心和南部、天津市南部、河北省西南部至東北部沿海地區、山東省少部分地區，呈斑塊狀零星分布，這與Zheng等[5]、Jiang等[22]、孫濤等[23]和趙安周等[24]的研究結果一致，由于北京、天津、石家莊等城市周邊的工業用地急劇擴張，耕地、有林地等生態空間大幅度縮減，人類活動對植被NDVI的變化產生負向影響[11,25-26]。

表1 環渤海地區植被NDVI變化趨勢面積比例統計

3.2 植被NDVI與人口密度變化的相關性

圖2 環渤海地區2000—2020年植被NDVI變化趨勢及顯著性檢驗

表2 環渤海地區人口密度冷熱點區域面積比例統計 %

3.2.2 植被NDVI與人口密度變化的相關性 為檢驗環渤海地區人口密度和植被NDVI的因果關系，本文利用格蘭杰因果關系檢驗法對兩者的相關關系進行驗證。該方法可判斷兩個時間序列之間是否存在因果性[30]。采用E-G兩步法[31]對人口密度和植被NDVI兩者的時間序列進行平穩性驗證，由協整檢驗一般規律可得，兩者的時間序列較平穩，存在協整關系進行格蘭杰因果檢驗。滯后期取1～6，選擇最優解作為結果進行分析。當p值<0.05時，拒絕原假設，即認定兩者之間存在單向格蘭杰原因[32]。根據表3可知，滯后期為2 a時，環渤海地區人口密度是植被NDVI的格蘭杰原因，但植被NDVI不是人口密度的格蘭杰原因，山東省人口密度是植被NDVI的格蘭杰原因，但植被NDVI不是人口密度的格蘭杰原因。滯后期為1 a時，北京市、河北省和遼寧省的人口密度是植被NDVI的格蘭杰原因，但植被NDVI不是人口密度的格蘭杰原因。滯后期為4 a時，天津市人口密度是植被NDVI的格蘭杰原因，但植被NDVI不是人口密度的格蘭杰原因。綜上所述，環渤海地區人口密度對植被NDVI存在單向格蘭杰原因，人口密度的改變是引起植被NDVI變化的格蘭杰原因，而植被NDVI的變化不是造成人口密度改變的格蘭杰原因。為探究植被NDVI變化與人口密度變化的空間耦合關系，本文將植被NDVI變化趨勢結果與人口密度冷熱點區進行疊置分析，得到植被NDVI與人口密度的空間耦合特征和數量關系，結果如圖4所示。2000—2010年各冷熱點區植被NDVI呈上升趨勢的面積遠大于呈下降趨勢的面積。植被NDVI呈顯著上升的面積占比隨人口密度變化熱點區到冷點區依次上升，其在冷點區面積占比最大，為37.11%，在熱點區面積占比最小，為13.83%。而植被NDVI呈顯著下降的面積占比在各冷點區面積占比與顯著上升面積占比成相反，顯著下降的面積占比在熱點區最大，為31.16%，在冷點區最小，為9.17%。由些可見，2000—2010年環渤海地區的植被覆蓋整體呈改善趨勢，其中冷點區改善程度最為明顯，次冷點區其次。

圖3 環渤海地區人口密度冷熱點區域分布

由圖4可見，2010—2020年植被NDVI呈顯著上升的面積占比隨人口集聚程度的下降呈現先減少后升高的趨勢，其中次熱點區面積占比最小，為3.97%，冷點區面積占比最大，為14.02%。同樣，植被NDVI呈顯著下降的面積占比呈現先減少后升高的趨勢，其中次冷點區面積占比最小，為4.81%，熱點區面積占比最大，為7.02%。相較于2000—2010年，2010—2020年植被NDVI變化趨勢的面積占比變化規律不顯著，且植被NDVI顯著上升的面積占比有所減少，植被NDVI顯著下降的面積占比增加。2000—2020年冷點區、次冷點區和次熱點區的植被NDVI呈上升趨勢變化的面積大于呈下降趨勢面積。植被NDVI呈顯著上升區域面積占比排序為：冷點區>次冷點區>次熱點區>熱點區，其中在冷點區的面積占比為65.14%，熱點區的面積占比為19.94%。由此可知，植被NDVI呈顯著上升的面積占比隨人口集聚程度的下降呈上升趨勢。植被NDVI顯著下降變化趨勢的面積占比排序為：熱點區>次熱點區>次冷點區>冷點區，其中熱點區的面積占比為18.50%，冷點區的面積占比為1.88%，植被NDVI呈顯著下降的面積占比隨人口集聚程度的下降呈現降低的趨勢。2000—2020年環渤海地區植被NDVI在次熱點區、次冷點區和冷點區呈改善趨勢，其中冷點區改善程度最為明顯，但熱點區的植被NDVI呈退化趨勢。綜上所述，環渤海地區人口密度變化是引起植被NDVI改變的格蘭杰原因，且對區域植被生長具有雙重作用。環渤海地區植被NDVI呈下降趨勢的面積主要分布在熱點區和次熱點區，原因是熱點區和次熱點區主要位于各大中型城市及其周邊地區，以上地區人口基數大，人口增長快，其中北京和天津的人口增長最為顯著，其他省會城市次之。另外，以上地區城市化程度高，如2019年北京、天津、遼寧和山東的城市化率在分別為86.6%，83.5%，68.1%和61.5%，高于全國平均水平，相較于2000年，城市化率分別增長了9.1%，12.3%，13.9%和23.5%(中國統計年鑒2001—2020)；其次，熱點區和次熱點區城市的第二、三產業對勞動力需求較大，對流動人口具有較強的吸引力，導致熱點區和次熱點區城市人口集聚程度較高，城市擴張，人類活動較為頻繁，對植被覆蓋造成破壞，從而抑制植被生長。環渤海地區植被NDVI呈上升趨勢的部分主要集中在冷點區，冷點區的地理位置缺乏優越區位條件，基礎設施較為落后，導致人口向北京市、天津市兩個直轄市以及石家莊、邯鄲、濟南、青島和沈陽等經濟核心城市流動。冷點區的經濟開發活動較低，對自然環境的負向干擾程度較弱，為了保障環渤海地區經濟的可持續發展，冷點區主要推行“三北”防護林工程、京津風沙源治理工程以及退耕還林還草工程等[6,8-9]，對植被生長起促進作用。

表3 人口密度和植被NDVI兩個時間序列的Granger檢驗結果

圖4 環渤海地區人口密度冷熱點區域植被NDVI變化顯著性面積比例統計

3.3 植被NDVI與人口密度重心遷移

3.3.1 植被NDVI重心遷移 由圖5可知，2000—2020年環渤海地區植被NDVI標準差橢圓的重心經度在118°39′54″—118°46′13″之間，重心緯度在39°03′40″—39°10′52″之間，表明環渤海地區植被NDVI呈“東北—西南”走向，且在“東北—西南”空間格局進一步擴大。長半軸標準差由4.890 km下降至4.846 km，而短半軸標準差由2.628 km上升至2.668 km，標準差橢圓變化較小，表明植被NDVI年際變化較為平穩。總體上，2000—2020年環渤海地區植被NDVI重心往西北方向遷移。研究區內植被NDVI的改變以及重心遷移受人類活動因子的影響[33]。主要有兩個方面原因：一方面，與環渤海地區，尤其是沿海地區，經濟發展和城市擴張緊密相關。隨著改革開放，以京津冀為經濟中心帶動輻射周圍城市經濟發展，城市用地的急劇擴張，人口密度顯著增加，導致耕地、草地和林地等轉為城市建設用地，人類活動頻繁，環渤海地區沿海地區的植被覆蓋率下降顯著，植被NDVI重心向西北方向遷移。另一方面，三北防護林工程和環渤海地區農業生態補償建設，實施封山育林、退耕還林等山區保護措施，退耕還林等舉措初見成效，增加了環渤海地區北部的森林覆蓋率，提高水土涵養能力，尤其是環渤海地區東北和西北部地區，植被覆蓋得到明顯的改善[34-35]，因此，環渤海地區植被NDVI重心向西北方向遷移。

3.3.2 人口密度重心遷移 2000—2020年環渤海地區人口密度標準差橢圓形成東北—西南的空間分布格局，標準差橢圓參數變化較為穩定，重心經度由118°05′38″逐步下降至117°58′30″，而重心緯度由38°12′00″逐步上升至38°17′21″，方位角由59.48°逐步擴大至61.32°。長半軸標準差從4.254 km縮小至4.061 km，表明環渤海地區人口密度在“東北—西南”方向上呈收斂趨勢，空間集聚性輕微上升，而短半軸標準差由2.405 km擴大至2.419 km，表明2000—2020年的環渤海地區人口密度在“西北—東南”方向上呈擴張趨勢，空間集聚性輕微下降。由圖6可知，2000—2020年環渤海地區人口密度變化較為平穩，重心逐步向西北方向遷移。影響環渤海地區人口密度重心遷移的主要因素是本地人口基數大以及外部人口的遷入[36]，使得人口密度逐年上升。遼中南經濟集聚區、京津唐工業區的城市經濟發展良好，城市化和工業化水平高，帶動了產業經濟由東南沿海逐步向西北內陸擴張，而產業發展和城市建設需更多人口彌補勞動力不足，吸引外來人口向該地區遷移，導致環渤海地區的人口重心向西北內陸遷移。2000—2020年的環渤海地區植被NDVI重心與人口密度重心均往西北方向遷移，表明人口重心遷移方向和NDVI重心遷移方向呈現一定相關性，人口密度重心遷移軌跡比較均一，單方向往西北遷移，植被NDVI重心整體往西北遷移，但遷移軌跡較為復雜，無明顯遷移模型，振蕩幅度較大，不與人口密度重心遷移軌跡呈明顯的線性關系。原因是除人口密度外，植被生長還與區域地形、氣候條件以及林業生態工程等息息相關，已有研究[11]表明，環渤海地區植被與降水、氣溫以及相對濕度呈明顯的正向效應，近年來，環渤海地區年際氣溫和降水呈波動上升趨勢，溫暖濕潤的環境對區域植被生長具有促進作用。此外，環渤海地區積極開展退耕還林還草工程、京津風沙源治理工程和三北防護林工程，以上工程不僅抵消人類活動對植被生長帶來的負向影響，還提高了環渤海地區植被覆蓋程度，增加了區域植被NDVI[37]。綜上所述，環渤海地區人口密度與植被NDVI重心變化具有一致性，但人口密度與植被NDVI變化沒有明顯的線性相關性，人口密度是影響環渤海地區植被NDVI的因素，但該地區的植被NDVI同時還受氣候變化和生態林業工程等因素的影響[6,8-9,11]。

圖5 環渤海地區2000—2020年植被NDVI重心轉移空間變化

圖6 環渤海地區2000—2020年人口密度重心轉移空間變化

4 結 論

(1) 2000—2020年環渤海地區植被NDVI整體呈上升趨勢，上升斜率為0.022/10 a。2000—2010年植被NDVI上升速率高于2010—2020年植被NDVI上升速率。2000—2020年植被NDVI呈上升趨勢的面積大于呈下降趨勢的面積，其中顯著上升區域面積占總面積的45.18%，主要分布于研究區的西北部和東北部。

(2) 2000—2020年環渤海地區人口密度變化以冷點區為主，主要從遼寧省東北部延伸至河北省西北部，以及山東省東南地區。相較于2000—2010年，2010—2020年人口密度變化的冷點區、次冷點區以及熱點區均有所減少，主要表現為次熱點區增加。

(3) 環渤海地區人口密度變化是引起植被NDVI變化的格蘭杰原因，而植被NDVI的變化不是造成人口密度改變的格蘭杰原因。且人口密度變化對植被NDVI變化的負向效應要強于正向效應，植被NDVI呈顯著下降趨勢面積主要集中于熱點區，呈顯著上升面積主要位于冷點區。

(4) 2000—2020年的環渤海地區植被NDVI重心與人口密度重心均往西北方向遷移，表明人口重心遷移方向和NDVI重心遷移方向呈現一定相關性，人口密度重心遷移軌跡比較均一，但植被NDVI重心遷移軌跡較為復雜，不與人口密度重心遷移軌跡呈明顯的線性關系。

本文的研究結果基本符合一般規律，表明利用人口密度分析研究對植被NDVI的研究方法切實可行，但仍存在不足之處。首先考慮的人為因素較少，人口密度作為單一因子進行分析導致的結論可能過于片面，后期的研究可以加入地區GDP、人口勞動力等人為因子[38]，從多角度分析人類活動對植被NDVI的影響。其次未采用最優滯后期模型計算人口密度與植被NDVI的最優滯后階，導致格蘭杰因果關系方法的檢驗結果比較粗略，在一定程度上影響評價結果的精度。因此，后續的研究應從適當增加研究因子、采用先進數學模型計算相關性等方面對人類活動與植被NDVI的相關性進行深入研究，使得研究結果與結論更加準確與科學。