李曉婷，楊麗帆，鄒友峰,3，馬 超,3

(1.河南理工大學(xué) 測繪與國土信息工程學(xué)院，河南 焦作 454003; 2.中國地質(zhì)大學(xué) 地理與信息工程學(xué)院，湖北 武漢 430074; 3.河南理工大學(xué) 自然資源部礦山時空信息與生態(tài)修復(fù)重點實驗室，河南 焦作 454003)

煤礦開采對塌陷區(qū)地表植被影響及恢復(fù)的程度、范圍和時間，一直是學(xué)術(shù)界爭論的焦點問題之一[1-2]。由于受迫植被生理生態(tài)指標(biāo)的甄別需要大范圍、長時間的野外調(diào)查和實驗室工作，相關(guān)研究進展緩慢。自20世紀(jì)70年代初以來，遙感已經(jīng)幫助我們極大地提高了對我們所居住的星球的認(rèn)知[3]，基于遙感手段的很多研究已經(jīng)提高了我們對植被響應(yīng)外在條件變化的理解[4]。近10 a來，基于星載、機載和地面平臺的多光譜、高光譜遙感技術(shù)的植被脅迫研究，一直是植物生態(tài)毒理學(xué)、病理學(xué)研究的熱點，國內(nèi)外關(guān)于重金屬脅迫如銅和鉛[5-7]、鎘[8-9]、鋅、砷、錳[10-11]等，水熱條件脅迫如降水[12]、干旱[13-15]、高溫[16]、凍害[17]等，以及農(nóng)藥脅迫[18]，交通脅迫[19]，鹽堿脅迫[20]，病蟲害脅迫[21-22]和酸性礦井水脅迫[23]導(dǎo)致的植被葉綠素、胡蘿卜素的熒光特征、反射光譜特征、紅邊響應(yīng)特征研究屢見報導(dǎo)。

開采沉陷對地表植被的直接損傷主要是根系損傷和土壤物理性質(zhì)的變化，沉陷造成土壤表面裂縫生成，土壤持水能力減弱，進而形成斷根脅迫和干旱脅迫[24];長期而言，更會造成土壤氮、磷、鉀等養(yǎng)分元素?fù)p失，廢矸廢水廢氣排放引起植被健康狀況惡化，導(dǎo)致光合作用能力變差，生產(chǎn)力下降[25-26]。開采沉陷造成植物生長層受到損傷，對植被的影響是持久的、長期的，植被作為陸地生態(tài)系統(tǒng)的主體，是聯(lián)結(jié)土壤、大氣和水分等要素的自然紐帶，是礦區(qū)生態(tài)環(huán)境變化最直觀的“指示器”。而歸一化植被指數(shù)NDVI是公認(rèn)的表征植被變化最有效的生態(tài)指標(biāo)之一，基于連續(xù)一致的NDVI時間序列數(shù)據(jù)，沉陷區(qū)植被生長狀況及其動態(tài)變化研究取得了眾多成果[27-31]。但由于植被監(jiān)測的多光譜遙感大多缺少紅邊波段，鮮有其他遙感生態(tài)指標(biāo)應(yīng)用于礦區(qū)植被生態(tài)脅迫的范疇，單一指標(biāo)表征礦區(qū)生態(tài)走勢尚顯薄弱。進入21世紀(jì)，以哨兵2號(Sentinel-2)衛(wèi)星為代表的一系列紅邊波段傳感器研究取得突破性進展，作為植被生長變化的最敏感的波段，一系列創(chuàng)新性紅邊指數(shù)被提出[32-36]。有研究表明，利用Sentinel-2A紅邊波段估算植被葉綠素時，其決定系數(shù)R2最高達到了0.99[37];利用Sentinel-2A數(shù)據(jù)估算的大豆紅邊葉綠素指數(shù)(CIre)與其葉綠素濃度具有非常高的相關(guān)性(R2= 0.94)[38];使用Sentinel-2A多光譜衛(wèi)星數(shù)據(jù)提取的紅邊葉綠素指數(shù)(CIre)、紅邊比值指數(shù)(SRre)和葉綠素敏感指數(shù)(MTCI)在冬小麥全生長季內(nèi)對生物量的變化保持高靈敏性[39]。

目前，基于紅邊波段的高敏感指數(shù)的采礦脅迫遙感研究尚屬空白，開采沉陷影響地表植被的機理有待深入研究。筆者以地處黃河流域中上游的鄂爾多斯市泊江海子礦為工程示范，針對113101 首采工作面地表覆被，開展一個生長周期的采礦脅迫實驗研究。利用歸一化植被指數(shù)(NDVI)和改進紅邊歸一化植被指數(shù)(Red Edge Normalized Difference Vegetation Index,NDVIre)、歸一化差值紅邊(Normalized Difference Red Edge,NDre)、修改的葉綠素反射率指數(shù)(Modified Chlorophyll Absorption in Reflectance Index,MCARI)、修改的紅邊簡單比值指數(shù)(Modified Red Edge Simple Ratio,MSRre)和三角葉綠素指數(shù)(Triangular Chlorophyll Index,TCI)共5種高敏感紅邊指數(shù)，診斷煤炭地下開采導(dǎo)致生長層生態(tài)損傷的因果關(guān)系，定量分析黃河流域脆弱生態(tài)區(qū)受迫植被退化程度與恢復(fù)過程，旨在闡明采礦與復(fù)墾活動對自然生態(tài)的擾動與植被生態(tài)系統(tǒng)修復(fù)的時空關(guān)系，期望豐富煤炭開采生態(tài)損傷與生態(tài)環(huán)境協(xié)同保護評價指標(biāo)體系，為黃河流域脆弱生態(tài)環(huán)境保護與生態(tài)重建提供一定的科學(xué)依據(jù)。

1 研究區(qū)與數(shù)據(jù)

1.1 研究區(qū)

泊江海子礦地處黃河流域中上游的黃土-風(fēng)積沙區(qū)，庫布齊沙漠以南，毛烏素沙地以北，資源同屬于我國十四大煤炭基地的神東基地(圖1)。行政隸屬內(nèi)蒙古鄂爾多斯市東勝區(qū)泊江海子鎮(zhèn)。地理經(jīng)度109.246°～109.417°，地理緯度39.838°～39.879°。井田東西長約20.16 km，南北寬約5.94 km，面積115.62 km2。113101工作面呈正南正北傾向長壁布置，于2016-12-01聯(lián)合試采。工作面地表為荒漠化地區(qū)，有少量農(nóng)田和陳家社、瓦窯社村。研究區(qū)氣候特點為高原大陸性干旱-半干旱氣候，日照時長較長，干燥，風(fēng)沙大，降水少[40]。

圖1 研究區(qū)地理地貌Fig.1 Geography and geomorphology in the study area

1.2 數(shù) 據(jù)

哨兵二號(Sentinel-2)A，B星同處一條軌道上，相差180°。其中Sentinel-2A于2015年6月發(fā)射。Sentinel-2A重訪周期為10 d。

載荷多光譜成像儀(Multi-Spectral Imager,MSI)的13個波段覆蓋了可見光-近紅外到短波紅外，掃描帶寬度290 km，具有3種分辨率(10，20和60 m)，見表1。

表1 Sentinel-2A衛(wèi)星多光譜波段信息Table 1 Multispectral band information of Sentinel-2A satellite

Sentinel-2A特有的3個紅邊波段(即B5，B6，B7)具有高空間分辨率、高光譜分辨率、寬幅寬、高重訪頻率特征，可為監(jiān)測地表植被的生長狀況提供更加有利的數(shù)據(jù)支持?？捎糜陉懙?、海洋、極地以及全球變化監(jiān)測，還可應(yīng)對洪澇、森林災(zāi)害和地震等突發(fā)情況[41]。

Sentinel-2A影像每月1次覆蓋研究區(qū)域。主要選擇天氣晴朗且云量少的高質(zhì)量影像。由于7月和9月云量覆蓋較大，影像質(zhì)量較差，因此未選用這2期數(shù)據(jù)。11期數(shù)據(jù)信息見表2。

表2 研究采用Sentinel-2A數(shù)據(jù)Table 2 Data of Sentinel-2A used in study

2 方 法

2.1 數(shù)據(jù)預(yù)處理

Sentinel-2A的Level-1C產(chǎn)品數(shù)據(jù)是經(jīng)過輻射校正和幾何精校正的大氣表觀反射率產(chǎn)品，沒有進行大氣校正。采用歐空局提供的Sentinel-2A/B大氣校正工具Sen2cor(Sentinel-2 Level-2A Atmospheric Correction Processor)實施高精度大氣校正[42]。Sen2cor可以將L1C產(chǎn)品處理成12個波段的L2A產(chǎn)品，同時將大氣上層表觀反射率轉(zhuǎn)化為大氣下層地表反射率(圖2)。

圖2 典型地物大氣校正前后的光譜曲線Fig.2 Spectral curves of typical surface features before and after Atmospheric Correction

圖2為典型地物(植被、水體、建筑物、裸地)大氣校正后效果，藍波段B2經(jīng)大氣校正后相比于大氣上層表觀反射率TOA顯著下降，紅邊、近紅外波段大幅增加，其光譜曲線與典型地物的光譜曲線的基本吻合。

大氣校正完成后，因不同波段的空間分辨率不同，利用Sentinel Application Platform(SNAP)對L2A數(shù)據(jù)采用最鄰近像元法進行重采樣，將所有的波段重采樣為20 m分辨率，并輸出為ENVI格式。最終去除B10波段(1 375 nm)，得到相同空間分辨率的12個波段。Sentinel-2A的1C和2A 級產(chǎn)品的反射率會乘上一個固定的系數(shù)。該系數(shù)默認(rèn)為10 000(從頭文件的QUANTIFICATION Value語句中查得)[43]。將重采樣后的結(jié)果導(dǎo)入ENVI進行波段融合(Layer Stacking)，得到一個多波段文件。將其各波段的灰度值DN除以10 000，還原為地表反射率。計算公式為

ρλ=Qcal/10 000

(1)

其中，ρλ為λ波段的地表反射率;Qcal為影像以16 bit量化的DN值。正常的DN值應(yīng)該在0～1，上式計算后如果存在DN<0,DN>1，可將小于0賦值為0，大于1的賦值為1。

2.2 植被紅邊指數(shù)

植被指數(shù)(VI)是一種利用遙感影像來監(jiān)測地表植被狀態(tài)的方法。紅光和紅外區(qū)間的反射率會隨綠色植被的不同覆蓋而變化。組合這2個波段信息可以用來區(qū)分土壤和植被。一定條件下，不同的植被指數(shù)在可以用來定量說明植被在某一階段的生長狀況。通過建立VI，綜合有關(guān)的光譜信號，使與植被有關(guān)信息增強，而使與植被無關(guān)信息減弱。

目前，在科學(xué)文獻中發(fā)布了超過150種植被指數(shù)模型，主要篩選出歸一化植被指數(shù)(NDVI)及5種紅邊植被指數(shù)—改進紅邊歸一化植被指數(shù)(Red Edge Normalized Difference Vegetation Index,NDVIre)、歸一化差值紅邊(Normalized Difference Red Edge,NDre)、修改的葉綠素反射率指數(shù)(Modified Chlorophyll Absorption in Reflectance Index,MCARI)、修改的紅邊簡單比值指數(shù)(Modified Red Edge Simple Ratio,MSRre)和三角葉綠素指數(shù)(Triangular Chlorophyll Index,TCI)。通過IDL批量計算11個期次的NDVI,NDVIre,NDre,MCARI,MSRre和TCI指數(shù)。NDVI和5種紅邊植被指數(shù)表達式，見表3。

表3 植被指數(shù)計算公式及對應(yīng)波段Table 3 Calculation formula and corresponding bands of Vegetation Index

由于Sentinel-2A的波段的中心波長與所選用植被指數(shù)對應(yīng)的中心波長不完全重合，因此利用與Sentinel-2A相鄰近波段進行波段運算見表3。

2.3 剖線與信息提取

113101工作面是泊江海子礦首采工作面，開采時間2016-12-01—2017-08-31，其沉陷區(qū)域隨開采由南向北移動。參考InSAR結(jié)果，沿113101工作面中心位置設(shè)計1條從北走向南(N—S走向)的剖面線，并在工作面東西兩側(cè)各設(shè)計1條平行的剖面線，剖面線長約6.6 km，每1條剖線包含334個像元。從西向東分別命名為對比區(qū)1(CK1)、沉陷區(qū)(S1)、對比區(qū)2(CK2)(圖3)，其中，圖3(a)中藍色箭頭為開采方向，圖3(b)中黃色箭頭為柵格排序方向，334個像元從北向南編號為1，2，…，334。

圖3 InSAR剖面線與光學(xué)影像疊加分析Fig.3 Superposition analysis of InSAR section line and optical image

2.4 趨勢分析與檢驗

根據(jù)其研究區(qū)剖線對應(yīng)的植被指數(shù)，利用長時序的遙感影像數(shù)據(jù)，建立沉陷區(qū)內(nèi)外植被指數(shù)隨時間的變化曲線。通過比較其植被指數(shù)在一時間段內(nèi)的各種變化檢測指標(biāo)的生長期曲線的差異，獲取研究區(qū)植被的物候信息。主要采用剖面分析法和非對稱高斯函數(shù)擬合法。

2.4.1剖面分析

剖面分析法主要為了更精確的分析沉陷區(qū)內(nèi)外植被分布情況，是多時相趨勢分析的數(shù)據(jù)基礎(chǔ)。根據(jù)植被指數(shù)NDVI的特點，當(dāng)旺季植被指數(shù)NDVI<0，則為水體;00.6為灌叢及林地[44]。

2.4.2回歸分析

(1)線性回歸。在剖面數(shù)據(jù)的基礎(chǔ)上，建立自變量與因變量之間的相關(guān)關(guān)系，從而建立一個具有良好相關(guān)性的回歸方程，即函數(shù)表達式。主要采用最小二乘法對數(shù)據(jù)進行一元線性回歸，回歸系數(shù)k主要是用來反映植被指數(shù)長期來的變化趨勢。

y=a+kxi+εi

(2)

式中，a為常數(shù)項;k為回歸系數(shù);εi為殘差。

(2)二元Logistic回歸。二元Logistic回歸是因變量為二分類的廣義線性回歸模型。二分響應(yīng)變量取值為0或者1，0表示事件未發(fā)生，1表示事件發(fā)生，正好用來檢驗自變量(如NDVI)是否受開采沉陷影響。我們設(shè)剖面上受采動影響的像元值為1，則未受采動影響的像元值為0，構(gòu)建如下logistic回歸模型[45]:

(3)

式中，Y為二分響應(yīng)變量;P為二分類變量取特定值時的概率;α0為常數(shù)項；Xi(i=1，2，…，m)為剖線上第i個像元的NDVI值;βi為自變量Xi的回歸系數(shù)。

2.4.3高斯擬合

高斯函數(shù)是表示連續(xù)隨機變量的概率密度函數(shù)。它通常用于表達自然物候的周期性演變，可以表示植被枯—榮—枯(或榮—枯—榮)過程[46]。函數(shù)(單峰型)的一維表達式為

(4)

式中，y0為基線的偏移量;A為鐘形曲線下側(cè)的積分面積(NDVI等效生物量);w大約等于峰高一半處寬度的0.849倍(當(dāng)峰高為1/2時，鐘形曲線的寬度),該值可以表示植被生長繁茂時的生長期的長度;x0為峰值位置。

2.4.4變化幅度

變化幅度指的是沉陷區(qū)變化值與非沉陷區(qū)變化值相比，增長或者下降的比率。函數(shù)表達式為

(5)

式中，V為變化幅度，當(dāng)V為正值時為增長幅度，為負(fù)值時為下降幅度；δ1為沉陷區(qū)差值；δ2為非沉陷區(qū)差值。

2.4.5決定系數(shù)

決定系數(shù)R2，又稱擬合優(yōu)度。R2越高，模型的擬合效果越好，自變量對因變量的解釋程度就越好，自變量本身所造成的變化百分比也就越高。R2的取值范圍(0，1)，沒有單位。若觀測點直接落在擬合曲線上，則R2=1，此時擬合曲線可以解釋因變量的所有變化。

(6)

(7)

其中,n為樣本數(shù)量;p為特征數(shù)量。即樣本為n個(x1，x2，x3，…,xp,y)。

3 結(jié)果與分析

3.1 時序NDVI剖面

獲取研究區(qū)2016年10月至2017年10月的植被指數(shù)時序變化剖面圖(S1，CK1，CK2)有3個作用:① 反映剖面植被類型及植被覆被情況;② 認(rèn)識剖面植被的物候特征及生物量(和值、均值)波動情況;③ 分析各剖面相關(guān)性，形成定量化比較的數(shù)據(jù)基礎(chǔ)。

NDVI是使用最廣泛的植被指數(shù)，已被證明是監(jiān)測植被狀況的最佳指標(biāo)之一。以此為例繪制S1，CK1，CK2時間序列剖面，如圖4所示，其中，黃色線段部分為開采沉陷影響區(qū)間，其余像元為非沉陷部分。

沉陷區(qū)剖面線(S1)的NDVI值時序圖(圖4(b))，該剖線上從12月開始形成沉陷，沉陷區(qū)域由南向北移動。其中，S為該期該剖線NDVI總和，M為均值。黃色線段部分為開采沉陷影響區(qū)間，其余像元為非沉陷部分。根據(jù)生長旺季(2017年8月)NDVI閾值，可知研究區(qū)地表覆被主要中等覆被類型的草地和濕草地;每年的4—10月為植被的生長期，S1的NDVI和值較高，但在S1開采前后沒有明顯變化，應(yīng)是季節(jié)變化掩蓋了采動影響;3條剖面NDVI同一季節(jié)相關(guān)性較高，說明地表覆被類型相近。

圖4 NDVI剖線時序Fig.4 Sequence diagram of NDVI section line

3.2 時序植被指數(shù)均值

將3條剖線(S1，CK1，CK2)上的植被指數(shù)的每1期均值(ΣVI/334)做一個散點圖，按照每月1期的時間間隔獲取該研究區(qū)2016年10月至2017年10月的植被指數(shù)時序變化(圖5)?？梢娫谝粋€生命周期內(nèi)的植被指數(shù)呈單峰高斯分布，因此采用高斯函數(shù)擬合曲線。

圖5 6種植被指數(shù)均值的高斯分布曲線Fig.5 Gaussian distribution curves of mean value of 6 vegetation indexes

根據(jù)變化趨勢可知，研究區(qū)所處的荒漠化草地在4月之后進入返青期，7，8月份植被的生長最為旺盛，大約在8月初達到峰值，10月之后進入衰退期。沉陷區(qū)與對比區(qū)的植被生長狀況有所差異，但差異不大。

表4 4—8月間6種植被指數(shù)的線性回歸Table 4 Linear regression of six Vegetation Indices from April to August

均值回歸結(jié)果表明，開采沉陷造成沉陷區(qū)植被的增長速度低于對比區(qū)，但絕對影響量有多少，有待開展增量研究。

3.3 時序植被指數(shù)差值

為了定量化表達沉陷前后植被指數(shù)的變化，設(shè)計了差值回歸方法。主要是針對采動影響，分析采動期和非采動期兩個時期植被指數(shù)的變化。即用采動時期的所有像元減去非采動時期相對應(yīng)像元的植被指數(shù)，求出采動引起的VI增量。

考慮到剖線上的植被分布狀況有所差異，僅利用S1剖面線的數(shù)據(jù)，由InSAR獲得的時間序列沉陷區(qū)可知，第1次衛(wèi)星過境時(2016-10-19)未發(fā)生沉陷，以這一期數(shù)據(jù)作為基準(zhǔn)數(shù)據(jù)，結(jié)合InSAR技術(shù)獲得的沉陷區(qū)具體位置，剖面線上落入沉陷區(qū)的像元位置見表5，其余的為非沉陷像元。

表5 沉陷區(qū)剖面線上涉及沉陷的像元位置Table 5 Pixel position of subsidence on the section line of subsidence area

用每一期沉陷區(qū)和非沉陷區(qū)像元對應(yīng)植被指數(shù)VI的均值，減去2016-10-19那一期對應(yīng)位置像元植被指數(shù)VI的均值，得到植被指數(shù)增量δVI，分析沉陷區(qū)和非沉陷區(qū)增量的時序變化趨勢。

以衛(wèi)星數(shù)據(jù)接收時間為橫軸，繪制植被指數(shù)增量的散點圖并采用高斯回歸(圖6)?？梢?種植被植數(shù)增量δNDVI，δNDVI，re，δND，re,δMCARI,δTCI和δMSR,re在沉陷區(qū)增量均小于非沉陷區(qū)增量。這表明受沉陷脅迫的影響，沉陷區(qū)植被的長勢要比非沉陷區(qū)差。這與3.2節(jié)結(jié)論一致。通過研究沉陷區(qū)相對于非沉陷區(qū)的下降幅度，發(fā)現(xiàn)6種植被指數(shù)增量δNDVI，δNDVI，re，δND，re,δMCARI,δTCI和δMSR,re平均下降幅度分別為24.1%，18.8%，25.9%，32.8%，28.4%，23.1%，沉陷區(qū)的植被的增長速率均低于非沉陷區(qū)的增長速率。

圖6 6種植被指數(shù)差值時序變化分析Fig.6 Time series analysis of difference of 6 Vegetation Indexes

對6種植被指數(shù)增量(δVI)時間序列進行了高斯回歸，δNDVI，δNDVI，re，δND，re,δMCARI,δTCI和δMSR,re回歸方程的決定系數(shù)分別為0.954 91，0.961 01，0.945 52，0.974 5，0.982 11和0.963 32。易見，當(dāng)NDVI中的近紅外波段替換為紅邊波段變?yōu)镹DVIre后，決定系數(shù)R2有明顯改善。

3.4 Logistic 回歸分析

為了確認(rèn)各植被指數(shù)是否受開采沉陷影響，筆者構(gòu)建了6種植被指數(shù)每一期的Logistic回歸模型，然后進行邏輯斯蒂回歸分析。檢驗情況見表6。其中顯著性水平Sig又稱為P值，當(dāng)P<0.05時，表明該模型是顯著的。

表6顯示，采動對各植被指數(shù)響應(yīng)在非生長季沒有顯著影響，即在非生長季其統(tǒng)計學(xué)意義不大;而在生長季會有顯著影響。且除MCARI，TCI以外的4種植被指數(shù)模型在植被生長較為旺盛的4，5，6，8月份均是顯著的，在植被生長最旺盛的8月份，6個植被指數(shù)均顯著相關(guān)。表明在生長季，采礦活動對植被指數(shù)有顯著影響。

表6 6種植被指數(shù)響應(yīng)采動影響的Logistic回歸的顯著水平Table 6 Significant level of Logistic regression of six vegetation indexes in response to mining

4 結(jié) 論

(1)在植被一個生命周期內(nèi)，6種植被指數(shù)值曲線呈明顯單峰變化。研究區(qū)地表的荒漠化草地在4月啟動返青期，7，8月份進入生長旺季，10月之后進入衰退期。

(2)對沉陷區(qū)(S1)、對比區(qū)1(CK1)、對比區(qū)2(CK2)的4—8月的6種植被指數(shù)均值進行線性回歸分析，沉陷區(qū)的回歸方程的斜率83%小于對比區(qū)(非沉陷區(qū))，沉陷區(qū)植被的增長速度慢于非沉陷區(qū)。

(3)沉陷區(qū)和非沉陷區(qū)1 a內(nèi)植被指數(shù)增量(δVI)的時序分析表明，沉陷區(qū)增量均小于非沉陷區(qū)增量。表明植被受沉陷脅迫，沉陷區(qū)植被的長勢要比非沉陷區(qū)差。采礦活動使6種植被指數(shù)增幅減少18.8%～32.8%。

(4)在植被生長期，二元Logistic回歸證實6種植被指數(shù)與采礦活動相關(guān)或顯著相關(guān)，采礦活動對植被指數(shù)有顯著影響。

黃河流域中上游的黃土-風(fēng)積沙礦區(qū)，氣候干旱，生態(tài)脆弱，植被紅邊指數(shù)對采礦脅迫極為敏感。這為研究黃河流域中上游煤礦生態(tài)修復(fù)的理論與方法，實現(xiàn)黃河流域高質(zhì)量發(fā)展，提供了敏感性極佳的生態(tài)指標(biāo)。本文的研究結(jié)果，回答了開采沉陷對地表植被有沒有影響？且部分回答了有多少影響？的問題。但對于影響多久？和長期影響規(guī)律如何？等問題的解答，仍有待長期實踐和深入研究。

致謝感謝中國-意大利合作培養(yǎng)博士、邊山大學(xué)訪問學(xué)者劉培副教授對英文摘要的潤色。