吳鳳敏，鄭稚棚，余 靜，梁 星，張靈犀，楊光譜，余 洋，趙珍妮

（1. 重慶市地理信息和遙感應(yīng)用中心，重慶 401147）

目前，對(duì)于植被遙感動(dòng)態(tài)監(jiān)測(cè)有大量研究[1-9]，其中Landsat 數(shù)據(jù)應(yīng)用最為廣泛，該類數(shù)據(jù)可以免費(fèi)下載，且適用于大尺度范圍的長(zhǎng)時(shí)間序列監(jiān)測(cè)[10-22]。基于此，研究以重慶市中梁山片區(qū)為實(shí)驗(yàn)區(qū)，利用Landsat8 數(shù)據(jù)對(duì)礦山植被長(zhǎng)時(shí)間序列情況進(jìn)行分析，同時(shí)結(jié)合高分?jǐn)?shù)據(jù)對(duì)礦山植被短時(shí)間序列進(jìn)行動(dòng)態(tài)監(jiān)測(cè)，并對(duì)于不同分辨率數(shù)據(jù)的植被監(jiān)測(cè)結(jié)果及關(guān)聯(lián)性進(jìn)行分析，希望能夠?yàn)楹罄m(xù)礦山修復(fù)治理過(guò)程中植被狀況動(dòng)態(tài)監(jiān)管提供技術(shù)支撐。

1 研究區(qū)概況及研究方法

1.1 研究區(qū)概況

中梁山片區(qū)地理坐標(biāo)為106°23′20″～106°24′25″E，29°26′10″～29°31′50″N，總面積為28.63 km2，最高點(diǎn)位于中部山頂，最高高程為695.8 m，最低點(diǎn)位于山腳，最低高程為270 m。研究區(qū)內(nèi)礦山有70個(gè)，總面積為191.36 hm2，位于中梁山背斜之南傾伏端，以低山地貌為主。研究選擇中梁山片區(qū)為示范區(qū)域，主要由于該區(qū)域露天礦山數(shù)量較多，存在已完成修復(fù)治理、正在修復(fù)治理和關(guān)閉未治理等多種類型，在礦山植被狀況動(dòng)態(tài)監(jiān)測(cè)研究中具有一定代表性（圖1）。

圖1 中梁山片區(qū)研究區(qū)域

研究利用landsat8 數(shù)據(jù)對(duì)2012—2021 年中梁山片區(qū)露天礦山植被指數(shù)進(jìn)行變化分析，利用高分?jǐn)?shù)據(jù)進(jìn)行2017—2021 年短時(shí)間序列變化分析。Landsat8 數(shù)據(jù)來(lái)源于地理空間數(shù)據(jù)云，高分?jǐn)?shù)據(jù)包括高分1號(hào)和高分6 號(hào)數(shù)據(jù)，來(lái)源于重慶市地理信息和遙感應(yīng)用中心。為盡量保證植被監(jiān)測(cè)年度生長(zhǎng)一致，Landsat8數(shù)據(jù)主要采用夏天時(shí)期影像， 數(shù)據(jù)獲取時(shí)間分別為2013-08-19、2014-07-21、2015-07-08、2016-07-11、2017-08-14、2018-05-13、2019-03-29、2020-09-07、2021-06-06，根據(jù)1∶5 000DEM數(shù)據(jù)計(jì)算該區(qū)域平均高程為0.37 km。高分?jǐn)?shù)據(jù)數(shù)據(jù)獲取從2017—2021年，部分?jǐn)?shù)據(jù)因?yàn)椴煌耆采w實(shí)驗(yàn)區(qū)需要進(jìn)行拼接，每景數(shù)據(jù)平均高程同樣根據(jù)1∶5 000DEM進(jìn)行計(jì)算（表1）。

表1 高分衛(wèi)星數(shù)據(jù)使用總體情況表

1.2 主要技術(shù)流程

研究采用ENVI5.3 軟件首先對(duì)Landsat8 數(shù)據(jù)和高分?jǐn)?shù)據(jù)進(jìn)行預(yù)處理，分別進(jìn)行植被指數(shù)計(jì)算，并對(duì)礦山植被狀況進(jìn)行現(xiàn)狀及變化分析，最后對(duì)兩者的關(guān)聯(lián)性進(jìn)行詳細(xì)對(duì)比分析（圖2）。研究選擇NDVI 作為礦山植被監(jiān)測(cè)指數(shù)，由于礦山生態(tài)修復(fù)過(guò)程中植被生長(zhǎng)變化較大，NDVI 能夠較好地反映植被指數(shù)從低值到高值區(qū)間的動(dòng)態(tài)變化，與植被分布密度呈線性相關(guān)。數(shù)據(jù)預(yù)處理主要包括輻射定標(biāo)、大氣校正、正射校正等流程。輻射定標(biāo)處理主要指將影像的灰度（D）值轉(zhuǎn)化成輻射亮度值或表觀反射率；大氣校正主要是將輻射定標(biāo)后的表觀反射率或輻射亮度值轉(zhuǎn)化為地表反射率的過(guò)程；正射校正需基于地形高程模型（DEM）來(lái)對(duì)影像中所有的像元執(zhí)行地形變形的校正，讓圖像滿足正射投影的要求，可根據(jù)影像自帶的RPC參數(shù)和ENVI中的DEM信息，對(duì)影像數(shù)據(jù)進(jìn)行正射校正。

圖2 Landsat8數(shù)據(jù)和高分?jǐn)?shù)據(jù)處理及關(guān)聯(lián)性分析技術(shù)流程

采用ENVI軟件進(jìn)行輻射定標(biāo)時(shí)，Landsat8數(shù)據(jù)可以直接采用自帶模型Radiometric Calibration進(jìn)行輻射定標(biāo)，高分?jǐn)?shù)據(jù)需要網(wǎng)上下載并安裝China Satellites Support插件，同時(shí)根據(jù)中國(guó)資源衛(wèi)星應(yīng)用中心發(fā)布的輻射定標(biāo)參數(shù)（http://www.cresda.com/CN/index.shtml）。由于NDVI計(jì)算取值范圍為[-1.0，1.0]，研究對(duì)該范圍外的作為異常值進(jìn)行剔除。根據(jù)已有學(xué)者研究成果，裸土NDVI 理論上接近于0，但由于受眾多因素影響，其取值范圍在-0.1～0.2之間。

2 Landsat植被動(dòng)態(tài)監(jiān)測(cè)

2.1 現(xiàn)狀情況分析

參考已有研究，按照無(wú)植被覆蓋區(qū)[-1.0，0.2），低植被覆蓋區(qū)[0.2，0.3），中等植被覆蓋區(qū)[0.3，0.4），高植被覆蓋區(qū)[0.4，1.0）進(jìn)行分級(jí)。結(jié)果顯示，2021年中梁山片區(qū)礦山高植被覆蓋區(qū)面積為81.09 hm2，占中梁山片區(qū)礦山總面積的42.38%；中等植被覆蓋區(qū)面積為29.35 hm2，占比15.34%；低植被覆蓋區(qū)面積為30.72%，占比16.05%，表明中梁山片區(qū)礦山植被整體情況較好，但存在部分礦山還在修復(fù)治理過(guò)程中，植被指數(shù)較小或者為推土情況無(wú)植被覆蓋（表2）。高植被覆蓋區(qū)中，NDVI 以0.4～0.5 區(qū)間為主，該區(qū)間礦山面積占高植被覆蓋區(qū)面積的50%以上，＞0.6的面積占比為46.88%。以單個(gè)礦山來(lái)看，中梁山片區(qū)70 個(gè)礦山中，平均NDVI 超過(guò)0.4 的礦山數(shù)量有36 個(gè)，平均NDVI 為0.3～0.4 的礦山為19 個(gè)，平均NDVI 為中等的礦山共有55 個(gè)，說(shuō)明隨著礦山修復(fù)治理工作的不斷開(kāi)展，礦山的植被指數(shù)狀況總體較好，大部分礦山植被恢復(fù)較好，部分礦山由于正處于修復(fù)治理過(guò)程中，可能存在部分區(qū)域推填土或者剛種植植被，植被指數(shù)較低。

表2 2013—2021年中梁山片區(qū)礦山植被指數(shù)NDVI分級(jí)表

2.2 年際變化分析

研究以每隔1 a NDVI值對(duì)礦山范圍內(nèi)數(shù)據(jù)進(jìn)行對(duì)比分析，2013—2021年中梁山片區(qū)植被指數(shù)NDVI呈現(xiàn)先減小后增加趨勢(shì)（圖3、4），主要波谷在2017年左右，由于礦山開(kāi)采時(shí)期植被破壞較大，植被指數(shù)較低，在修復(fù)治理過(guò)程中，也可能存在部分區(qū)域推填土破壞植被的情況；隨著礦山不斷修復(fù)，前期的植被指數(shù)從較低的覆蓋指數(shù)逐漸增加，最終修復(fù)完成后植被不斷恢復(fù)，植被指數(shù)不斷增加，表明礦山生態(tài)修復(fù)治理取得了較好的成績(jī)。其中，礦山高植被覆蓋區(qū)（NVDI＞0.4）從2013 年的104.81 hm2下降至2017 年的42.17 hm2，減少面積為62.64 hm2；從2017—2021 年面積又逐漸增加，2021年的高植被覆蓋區(qū)面積達(dá)到81.09 hm2。低植被覆蓋區(qū)（0.2＜NDVI≤0.3） 具有同樣的規(guī)律，2013—2017 年面積減少12.64 hm2，2017 年到2021 年面積增加17.54 hm2。2017—2021年間高植被覆蓋區(qū)增加的趨勢(shì)較低植被覆蓋區(qū)更大，說(shuō)明這個(gè)期間為礦山植被恢復(fù)的主要時(shí)期，最開(kāi)始以低植被覆蓋的治理初期隨著時(shí)間的推移，植被逐漸變好，并且變?yōu)楦咧脖桓采w區(qū)較多。

圖3 2013年和2021年中梁山片區(qū)植被指數(shù)NDVI分布圖

研究發(fā)現(xiàn)，存在個(gè)別礦山植被指數(shù)到2021年總體情況較低，比如中梁山片區(qū)萌特礦山，經(jīng)實(shí)地調(diào)查，該礦山于2019年6月開(kāi)始工程治理，目前正處于修復(fù)治理過(guò)程中，項(xiàng)目周期為66 個(gè)月，到2025 年為止，存在部分區(qū)域因修復(fù)治理過(guò)程中對(duì)周邊土地的踩踏以及修復(fù)治理最開(kāi)始種植小樹(shù)苗在影像上不明顯，因此總體上植被指數(shù)較低。

圖4 中梁山片區(qū)礦山2013—2021年植被指數(shù)變化情況

2.3 變化斜率分析

研究采用一元線性回歸方法模擬每個(gè)柵格的NDVI變化趨勢(shì)，重點(diǎn)分析礦山范圍內(nèi)不同修復(fù)治理區(qū)域植被變化情況。對(duì)于NDVI 時(shí)間序列數(shù)據(jù)，每個(gè)像元對(duì)應(yīng)有若干年的時(shí)間序列數(shù)值，對(duì)這些數(shù)值進(jìn)行線性擬合，獲取直線的斜率，像元斜率表明了在該時(shí)間序列中該像素所代表的植被指數(shù)的演化趨勢(shì)。主要計(jì)算公式如下：

式中，k為斜率；xi為具體年份；yi為每年的NDVI值。

2017—2021 年中梁山片區(qū)礦山范圍內(nèi)斜率＞0 的像元總面積為85.15 hm2，占礦山總面積的44.49%，表明這些區(qū)域內(nèi)植被逐漸變好；斜率＜0 的像元面積為106.21 hm2，表明這些區(qū)域植被狀況有變差的趨勢(shì)。參考已有研究[13]，將NDVI 變化程度分為7 級(jí)（表3），圖5 為以30 m×30 m 像元大小，統(tǒng)計(jì)中梁山片區(qū)礦山范圍內(nèi)所有像元的斜率值分布情況。中梁山片區(qū)礦山范圍內(nèi)像元明顯改善面積有55.22 hm2，占像元總面積的28.85%；輕微改善及中度改善的面積合計(jì)為26.14 hm2，占比13.66%，存在一般以上斜率為負(fù)值，植被情況較差的情況，主要原因是中梁山片區(qū)面積最大的礦山重慶萌特礦山處于修復(fù)治理過(guò)程中。雖然有一部分區(qū)域在種植樹(shù)木，但由于數(shù)據(jù)還在恢復(fù)中，植被指數(shù)變化較小，且治理過(guò)程中存在部分推填土情況，因此植被指數(shù)總體變小。

表3 2013—2021年中梁山片區(qū)礦山像元NDVI斜率分級(jí)（Landsat8數(shù)據(jù)）

圖5 中梁山片區(qū)礦山范圍內(nèi)單個(gè)像元斜率分布情況

研究通過(guò)ArcGIS 工具將計(jì)算斜率值與中梁山片區(qū)礦山像元進(jìn)行空間關(guān)聯(lián)，可以看出中梁山片區(qū)中部及東部礦山修復(fù)情況較好，植被斜率大部分為正數(shù)，右下角萌特礦山和中部部分礦山目前還處于修復(fù)治理過(guò)程中或者關(guān)閉未治理情況，植被覆蓋較差，如圖6所示。

圖6 2013—2021年中梁山片區(qū)礦山范圍內(nèi)像元植被指數(shù)斜率分布圖

3 高分?jǐn)?shù)據(jù)植被動(dòng)態(tài)監(jiān)測(cè)

3.1 現(xiàn)狀情況分析

研究按照與Landsat8 數(shù)據(jù)同樣分析方法對(duì)高分?jǐn)?shù)據(jù)計(jì)算的植被指數(shù)進(jìn)行分析，由于分辨率較高（8 m），礦山植被指數(shù)監(jiān)測(cè)結(jié)果較landsat8 數(shù)據(jù)更精細(xì)。中梁山片區(qū)整體植被情況較好，大部分區(qū)域均為高植被覆蓋區(qū)域（NDVI ＞0.4），70個(gè)礦山總體上位于中高植被覆蓋區(qū)域范圍內(nèi)，2021 年有植被覆蓋區(qū)域面積為150.49 hm2，其中高植被覆蓋區(qū)域94.84 hm2，占礦山總面積的49.56%。與landsat8 計(jì)算結(jié)果相比，高植被覆蓋區(qū)面積較大，說(shuō)明高分?jǐn)?shù)據(jù)對(duì)于植被指數(shù)的響應(yīng)更敏感。中等植被覆蓋區(qū)面積為29.48 hm2，占比為15.41%（圖7）。

圖7 2021年中梁山片區(qū)高分?jǐn)?shù)據(jù)NDVI監(jiān)測(cè)分布圖

3.2 變化斜率分析

由于高分?jǐn)?shù)據(jù)只有2017—2021 年5 a 的數(shù)據(jù)，因此進(jìn)行短時(shí)間序列斜率變化分析（表4，圖8），圖8為以8 m×8 m像元大小，統(tǒng)計(jì)中梁山片區(qū)礦山范圍內(nèi)所有像元的斜率值分布情況。結(jié)果表明中梁山片區(qū)礦山范圍內(nèi)NDVI斜率＞0的像元面積為78.94 hm2，占礦山總面積的41.25%；斜率＜0的像元面積為112.42 hm2，總體上超過(guò)礦山總面積的一半，NDVI 斜率的分布情況與Landsat8 數(shù)據(jù)分析結(jié)論基本一致，僅部分分級(jí)數(shù)值略有差異。

圖8 中梁山片區(qū)礦山范圍內(nèi)單個(gè)像元斜率分布情況（高分?jǐn)?shù)據(jù)）

表4 2017—2021年中梁山片區(qū)礦山像元NDVI斜率分級(jí)（GF數(shù)據(jù)）

從分級(jí)情況看，礦山范圍內(nèi)嚴(yán)重退化區(qū)域面積最大，為108.51 hm2，主要是由于中梁山片區(qū)萌特礦山正在修復(fù)治理過(guò)程中造成。明顯改善區(qū)域面積居第二，占礦山總面積的39.45%，表明其他大部分礦山區(qū)域的植被都是變好的，說(shuō)明礦山修復(fù)治理取得了良好的效果。

4 植被指數(shù)關(guān)聯(lián)性分析

研究選擇2019 年數(shù)據(jù)進(jìn)行關(guān)聯(lián)性分析，利用Landsat8 （2019-03-29） 和高分?jǐn)?shù)據(jù)（2019-07-28）分別計(jì)算NDVI，得到兩者空間分布圖。結(jié)果顯示中梁山片區(qū)范圍內(nèi)高分?jǐn)?shù)據(jù)在高植被覆蓋區(qū)域分級(jí)粒度更細(xì)，NDVI 大于0.8 在圖上分布較為明顯，小于0.2區(qū)域面積總體較少，主要由于高分?jǐn)?shù)據(jù)分辨率為8 m，而Landsat8分辨率為30 m，在像元尺度上部分高值和低值區(qū)域進(jìn)行了一定綜合，因此在中間部分值的像元數(shù)量較多（圖9，表5）。

表5 中梁山片區(qū)礦山2019年Landsat8和高分?jǐn)?shù)據(jù)計(jì)算NDVI分級(jí)情況對(duì)比

圖9 2019年Landsat8和高分?jǐn)?shù)據(jù)計(jì)算NDVI空間分布圖對(duì)比

通過(guò)對(duì)中梁山片區(qū)礦山范圍內(nèi)的兩種影像計(jì)算的DNVI 進(jìn)行計(jì)算得出，高分?jǐn)?shù)據(jù)有植被覆蓋區(qū)域（NDVI ＞0.2）為82.68 hm2，明顯小于Landsat8計(jì)算出的有植被覆蓋區(qū)域面積（132.93 hm2），分析造成趨勢(shì)情況不一致的原因有以下幾種：一是分辨率原因?qū)е虏糠諰andsat8 像元綜合了周?chē)腘DVI 值較高的植被區(qū)域?qū)е孪裨獌?nèi)植被指數(shù)整體較高；二是由于兩類NDVI數(shù)據(jù)在進(jìn)行正射糾正時(shí)均采用ENVI自帶的DEM數(shù)據(jù)，可能存在部分像元位置有一定偏移；三是兩類NDVI 數(shù)據(jù)的時(shí)相不一致，對(duì)于植被生長(zhǎng)期情況的計(jì)算結(jié)果本身也有一定差異。

研究利用中梁山片區(qū)范圍對(duì)2 種數(shù)據(jù)進(jìn)行植被指數(shù)關(guān)聯(lián)性分析，通過(guò)ArcGIS 中Spatial Analyst Tools 工具進(jìn)行計(jì)算。計(jì)算結(jié)果顯示，通過(guò)Landsat8 和高分?jǐn)?shù)據(jù)計(jì)算NDVI 相關(guān)系數(shù)為0.739 3，表明2 種遙感影像在計(jì)算結(jié)果上呈現(xiàn)較高的正相關(guān)性。同時(shí)，研究以Landsat8 像元為網(wǎng)格（30 m × 30 m），計(jì)算高分?jǐn)?shù)據(jù)NDVI 平均值，分析空間上同一像元NDVI 值差異情況，結(jié)果表明，高分?jǐn)?shù)據(jù)與Landsat8 數(shù)據(jù)計(jì)算NDVI值呈現(xiàn)花瓣?duì)钕嚓P(guān)特征，兩者在高值區(qū)和低值區(qū)呈聚攏狀態(tài)，在中間值區(qū)呈離散狀態(tài)，總體上表現(xiàn)為中間寬兩頭窄的花瓣形。以線性關(guān)系式表達(dá)兩者特征，決定系數(shù)R2=0.692，表明兩者雖然對(duì)于同一像元的值存在差異，但總體上還是具有較好的線性相關(guān)性（圖10）。

圖10 2019年Landsat8和高分?jǐn)?shù)據(jù)計(jì)算NDVI相關(guān)性分析

5 結(jié)論與討論

通過(guò)對(duì)Landsat8 數(shù)據(jù)和高分?jǐn)?shù)據(jù)進(jìn)行預(yù)處理，對(duì)礦山修復(fù)治理過(guò)程中植被狀況進(jìn)行時(shí)間序列的動(dòng)態(tài)分析，并對(duì)2 種數(shù)據(jù)計(jì)算處理的植被指數(shù)進(jìn)行關(guān)聯(lián)性分析，希望通過(guò)礦山植被遙感監(jiān)測(cè)能夠?qū)ξ磥?lái)礦山修復(fù)治理監(jiān)管提供基礎(chǔ)支撐。本研究仍然存在一定不足之處，需要進(jìn)一步探討。