包妮沙 李秋玥 楊天鴻 劉善軍 齊 跡

(1.東北大學(xué)資源與土木工程學(xué)院,遼寧 沈陽 110819;2.遼寧工程技術(shù)大學(xué)測繪與地理科學(xué)學(xué)院,遼寧 阜新 123000)

新疆地區(qū)礦產(chǎn)資源豐富,有色金屬礦和煤炭資源儲量均位于我國前列。人類經(jīng)濟(jì)活動與資源環(huán)境保護(hù)矛盾已成為制約干旱半干旱區(qū)域可持續(xù)發(fā)展的關(guān)鍵問題[1],新疆作為典型的干旱半干旱地區(qū)生態(tài)環(huán)境脆弱,礦產(chǎn)資源開發(fā)對區(qū)域生態(tài)環(huán)境有很大沖擊,露天開采尤為嚴(yán)重。據(jù)2020年新疆地質(zhì)環(huán)境監(jiān)測院調(diào)查統(tǒng)計:新疆地區(qū)因采礦破壞土地面積達(dá)9.07×104km2;水土流失面積達(dá)103×104km2,占全國28.1%;近2/3的土地遭受荒漠化危害。露天開采全生命周期下生態(tài)環(huán)境的動態(tài)監(jiān)測是指導(dǎo)綠色礦山建設(shè)和生態(tài)恢復(fù)的重要基礎(chǔ)。

傳統(tǒng)的礦區(qū)生態(tài)環(huán)境調(diào)查多以地面調(diào)查、采樣、化學(xué)測試方法獲取生態(tài)指標(biāo),實(shí)現(xiàn)礦區(qū)生態(tài)環(huán)境監(jiān)測,難以滿足多尺度、快速動態(tài)監(jiān)測需求。遙感技術(shù)憑借周期短、實(shí)時性、動態(tài)性和信息豐富等特點(diǎn),通過地面目標(biāo)反射、輻射特性表征地表生態(tài)環(huán)境,不僅能全面真實(shí)地記錄地表覆蓋信息,揭示礦區(qū)邊開采、邊復(fù)墾過程中生態(tài)環(huán)境的陸面演變規(guī)律,而且對露天開采過程中引起的突發(fā)環(huán)境問題能夠快速響應(yīng)和及時預(yù)警,以揭示受采礦等人類活動影響的生態(tài)環(huán)境變化。多波段遙感數(shù)據(jù)衍生的地表生態(tài)狀況指數(shù),如地表溫度(LST)、歸一化植被指數(shù)(NDVI)等已被廣泛應(yīng)用于礦區(qū)生態(tài)環(huán)境監(jiān)測中[2-4],反映露天開采引起的地表溫度、植被覆蓋的變化,為分析礦區(qū)地表生態(tài)擾動提供了更多可用信息。

單一生態(tài)因子遙感反演指標(biāo)只能表征生態(tài)系統(tǒng)一個方面的變化,難以應(yīng)用于復(fù)雜的生態(tài)系統(tǒng)和揭示礦山環(huán)境的改變[5]。因此,可以識別更多特征的多生態(tài)因子耦合得到的遙感模型或指數(shù)成為研究的熱點(diǎn)[6-7]。徐涵秋[8]于2013年提出了遙感生態(tài)指數(shù)(RSEI),通過主成分分析(PCA)耦合了綠度、濕度、干度和熱度4個生態(tài)因子,是近年來應(yīng)用最廣泛的遙感模型之一,能夠應(yīng)用于礦區(qū)的生態(tài)累積效應(yīng)和生態(tài)環(huán)境質(zhì)量監(jiān)測。隨著應(yīng)用場景和目的不同,許多學(xué)者嘗試對其方法進(jìn)行了改進(jìn),劉英等[9]認(rèn)為RSEI中PCA耦合方式忽略4個生態(tài)因子間的弱線性或非線性的影響,因此采用非線性的核主成分分析(KPCA)進(jìn)行改進(jìn),構(gòu)建了非線性的遙感生態(tài)指數(shù)(NRSEI)。宋美杰等[10]認(rèn)為RSEI中僅用到PCA第一主成分,存在信息利用不充分的問題,因此采用前3個主成分分量來集成MRSEI。李晶等[11]認(rèn)為RSEI對黃河流域水土侵蝕易發(fā)區(qū)的適用性尚不明確,除了以上4個生態(tài)因子外,加入土壤侵蝕度因子,耦合5個生態(tài)因子對山西省礦區(qū)進(jìn)行了生態(tài)環(huán)境監(jiān)測。ZHENG等[12]認(rèn)為RSEI在監(jiān)測長時序生態(tài)環(huán)境質(zhì)量變化方面不穩(wěn)定,因此在PCA之前,采用標(biāo)準(zhǔn)化處理方法代替歸一化將4個生態(tài)因子的值轉(zhuǎn)換至一個公共域,構(gòu)建了標(biāo)準(zhǔn)化遙感生態(tài)指數(shù)(RSEIs)。以上改進(jìn)均是對4個生態(tài)因子去量綱或者耦合方式進(jìn)行調(diào)整,但是2022年RSEI指數(shù)的提出者徐涵秋教授在《武漢大學(xué)學(xué)報》《遙感技術(shù)與應(yīng)用》等刊物上撰文[13-14]針對以上改進(jìn)提出了質(zhì)疑,指出這些改進(jìn)指數(shù)普遍存在算法論證不充分、精度驗(yàn)證不合理等不足,如MRSEI指數(shù)僅增加信息量而忽略了PCA各分量的生態(tài)意義,會造成結(jié)果失真;RSEIs指數(shù)中標(biāo)準(zhǔn)化方法難以避免異常值的影響。

由于RSEI是基于綠度、濕度、干度和熱度4個生態(tài)因子耦合的模型,KARIMI等[15]指出干度指標(biāo)的計算是通過裸土指數(shù)和建筑指數(shù)實(shí)現(xiàn),采礦活動的增加會改變自然表面反射率、電導(dǎo)率、土壤pH值等土壤物理化學(xué)特性,導(dǎo)致自然表面向不透水面轉(zhuǎn)化,因此裸土和受采礦活動影響土地最重要的區(qū)別是表面不透水率百分比,干度指標(biāo)難以有效區(qū)分采礦活動引起的地表變化。對此,KARIM等[15]選擇了亮度指標(biāo)代替干度指標(biāo),通過代數(shù)運(yùn)算方式耦合4個生態(tài)因子,提出了一種新的指數(shù),即LSESCI,并選擇了14個不同類型的研究區(qū),對比LSESCI與RSEI指數(shù)的應(yīng)用效果,發(fā)現(xiàn)相對于RSEI指數(shù),LSESCI指數(shù)在礦區(qū)生態(tài)環(huán)境質(zhì)量監(jiān)測方面的性能較優(yōu)。

由于生態(tài)環(huán)境背景及礦種類型存在差異,導(dǎo)致采礦活動誘發(fā)的生態(tài)環(huán)境問題具有明顯的區(qū)域異質(zhì)性?；谝陨仙鷳B(tài)因子構(gòu)建的指數(shù)模型均在單一礦區(qū)進(jìn)行對比,在不同景觀條件和露天開采模式下表征礦區(qū)生態(tài)環(huán)境的差異性及穩(wěn)定性尚不明確。遙感模型的穩(wěn)定性與普適性,是利用多源遙感數(shù)據(jù)實(shí)現(xiàn)礦區(qū)多要素參數(shù)同步反演與構(gòu)建監(jiān)測技術(shù)體系的基礎(chǔ),同時也是科學(xué)認(rèn)識和評價礦區(qū)生態(tài)環(huán)境狀況的關(guān)鍵[16]。新疆全區(qū)面積166 萬km2[17],涉及數(shù)十個生態(tài)功能區(qū)。針對礦區(qū)所在區(qū)域的生態(tài)環(huán)境特點(diǎn),綜合考慮背景環(huán)境中土壤、植被以及巖礦類型,選擇區(qū)域尺度科學(xué)、有效、合理的地表生態(tài)狀況指數(shù),揭示露天開采過程中對區(qū)域生態(tài)環(huán)境的擾動性,使其評價結(jié)果具有穩(wěn)定性和可比性,是新疆地區(qū)綠色礦山建設(shè)及評價的關(guān)鍵基礎(chǔ)。本研究選擇了位于高寒高海拔山前沖洪積扇區(qū)的金寶鐵礦區(qū)、荒漠戈壁區(qū)的黑山煤礦區(qū)以及荒漠綠洲區(qū)的烏拉根鋅礦區(qū),利用2005—2020年長時序Landsat衛(wèi)星影像為數(shù)據(jù)源,提取濕度、熱度、干度、綠度、亮度5個生態(tài)因子,構(gòu)建LSESCI、RSEI、MRSEI及RSEIs 4個指數(shù)模型,從指數(shù)模型與各生態(tài)因子相關(guān)性、土地利用及地形因素影響3個方面,闡明各指數(shù)在表征干旱半干旱區(qū)露天礦地表生態(tài)狀況性能差異性和穩(wěn)定性,評價其在不同礦區(qū)的應(yīng)用效果和適用性,并選擇適宜的指數(shù)揭示露天開采過程中地表生態(tài)環(huán)境質(zhì)量的時空變化特征。

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

1.1 研究區(qū)概況

本研究選擇了新疆位于山前沖洪積扇區(qū)、荒漠戈壁區(qū)、荒漠綠洲區(qū)的3座典型露天礦進(jìn)行分析(圖1)。

圖1 研究區(qū)分布示意Fig.1 Distribution of the study area

金寶鐵礦區(qū)礦山產(chǎn)能約400 萬t/a,屬于大型露天金屬礦山。礦區(qū)位于阿勒泰低中山區(qū)山坡地貌,以山地為主且其中鑲嵌著山間盆地、谷地,地勢復(fù)雜,屬于寒旱牧區(qū)。年平均氣溫1.9 ℃,年平均降水量271.31 mm,年平均蒸發(fā)量1 477.2 mm,海拔946～1 570 m,屬大陸性中溫帶干旱氣候。

黑山煤礦區(qū)屬于準(zhǔn)東煤田,探明儲量469 億t,設(shè)計生產(chǎn)能力為1 000 萬t/a,屬特大型露天煤礦。礦區(qū)位于沖積戈壁平原荒漠地帶,以戈壁為主,兼有沙地、裸巖石質(zhì)地等荒漠用地類型,也有少量荒漠草地、灌叢等荒漠植被,地形整體平坦開闊。區(qū)內(nèi)年平均氣溫為1.2～2.3 ℃,年均降水量約215.3 mm,年均蒸發(fā)量1 099.2 mm,海拔2 322～2 845 m,屬于典型的荒漠礦區(qū)。

目前,烏拉根鋅礦區(qū)已探明鉛鋅金屬儲量超過500 萬t,預(yù)測的鉛鋅金屬資源量在1 000萬t以上,開采方式為露天—地下聯(lián)合開采。礦區(qū)位于荒漠綠洲區(qū),土壤類型為棕漠土,鹽漬化敏感,植被以天然牧草地為主。研究區(qū)平均氣溫7.4 ℃,年平均降水量181.6 mm,年平均蒸發(fā)量為2 600.5 mm,為降水量的14.3倍,海拔1 932～2 493 m。

1.2 數(shù)據(jù)源

本研究遙感數(shù)據(jù)源為GEE平臺2005年至2020年6—10月Landsat TM和OLI地表反射率數(shù)據(jù)集,共計145景影像,其中金寶鐵礦區(qū)55景,黑山煤礦區(qū)55景,烏拉根鋅礦區(qū)35景(圖2)?；贕EE平臺,利用CFMask算法對影像進(jìn)行去云處理。對同一年中涉及的所有月份數(shù)據(jù)提取生態(tài)因子,最終通過取均值的方式,表征這一年的生態(tài)環(huán)境質(zhì)量狀況,從而減少因?yàn)榧竟?jié)性差異而降低生態(tài)環(huán)境質(zhì)量時空變化特征監(jiān)測的準(zhǔn)確性。

圖2 研究區(qū)影像數(shù)量統(tǒng)計Fig.2 Statistics of the number of images in the study area

土地利用數(shù)據(jù)參考2020年ESRI全球10 m分辨率的土地覆蓋數(shù)據(jù),由ArcGIS Living Atlas 中的Web服務(wù)訪問下載,總體分類精度可達(dá)85%。按照該數(shù)據(jù)的分類體系和定義,研究區(qū)涉及的土地利用類型有建筑、裸地、牧場及林地。裸地指礦場、沙漠和沙丘等,牧場指低植被覆蓋度、可暴露出土壤和巖石的開闊區(qū)域,林地指較高大的植被,建筑指人造建筑物[18]。結(jié)合礦區(qū)典型地物和土地利用特點(diǎn)[16],將裸地重新劃分為開采區(qū),牧場劃分為裸土,林地劃分為植被。數(shù)字高程數(shù)據(jù)來源于地理空間數(shù)據(jù)云 (http:∥www.gscloud.cn/) ASTER DEM數(shù)據(jù),空間分辨率為30 m。

2 研究方法

本研究首先提取濕度、熱度、干度、綠度、亮度生態(tài)因子,進(jìn)而耦合構(gòu)建LSESCI、RSEI、MRSEI以及RSEIs指數(shù);其次結(jié)合土地利用數(shù)據(jù)、DEM數(shù)據(jù),針對礦區(qū)所在區(qū)域生態(tài)環(huán)境特點(diǎn),通過相關(guān)性、緩沖區(qū)及滑動濾波等綜合分析指數(shù)間的性能差異;最后選擇最優(yōu)指數(shù)利用Sen+Mann-Kendall時間序列分析法揭示露天礦生態(tài)環(huán)境變化趨勢和演變規(guī)律,如圖3所示。

圖3 研究流程Fig.3 Study process

2.1 地表生態(tài)狀況指數(shù)構(gòu)建方法

2.1.1 生態(tài)因子提取

2.1.1.1 濕 度

纓帽變換[19]是KAUTH與THOMAS于1976年根據(jù)多光譜遙感中植物、水體等地物的分布特點(diǎn)提出的一種經(jīng)驗(yàn)線性變換。衛(wèi)星遙感影像經(jīng)過纓帽變換之后生成的前3個分量與地面景物密切相關(guān),第1分量為反映地物總體反射率強(qiáng)弱的亮度指數(shù);第2分量為反映地表植被覆蓋率等信息的綠度指數(shù);第3分量為反映土壤等水分狀態(tài)的濕度指數(shù)。本研究中濕度采用纓帽變換中的濕度分量來代表,其表達(dá)式為[20-21]

式中,Blue、Green、Red、NIR、SWIR1、SWIR2為藍(lán)、綠、紅、近紅外、短波紅外1和短波紅外2波段,分別對應(yīng)TM影像中的第1、2、3、4、5、7波段和OLI影像中的第2、3、4、5、6、7波 段。對 于TM影 像[22],系 數(shù)ωi(i=1,2,…,6)取值分別為0.031 5、0.202 1、0.310 2、0.159 4、-0.680 6、-0.610 9;對 于OLI影 像[23],ωi(i=1,2,…,6)取值分別為0.151 1、0.197 3、0.328 3、0.340 7、-0.711 7、-0.455 9。

2.1.1.2 綠 度

以應(yīng)用最廣泛的歸一化植被指數(shù)NDVI表征RSEI中的綠度指標(biāo),公式為

式中,Red、NIR為紅和近紅外波段,分別對應(yīng)TM影像中的第3、4波段和OLI影像中的第4、5波段。

采用纓帽變換中的綠度分量來表征LSESCI中的綠度指標(biāo),公式為

式中,Blue、Green、Red、NIR、SWIR1、SWIR2為藍(lán)、綠、紅、近紅外、短波紅外1和短波紅外2波段,分別對應(yīng)TM影像中的第1、2、3、4、5、7波段和OLI影像中的第2、3、4、5、6、7波段。對于TM影像,系數(shù)ωi(i=1,2,…,6)取值分別為-0.160 3、-0.281 9、-0.493 4、0.794 0、0.000 2、-0.144 6;對于OLI影像,系數(shù)ωi(i=1,2,…,6)取值分別為-0.294 1、-0.243 0、-0.542 4、0.727 6、0.071 3、-0.160 8。

2.1.1.3 熱 度

熱度指標(biāo)由經(jīng)比輻射率校正的地表溫度來代表,將遙感影像熱紅外波段反演為亮溫再對其進(jìn)行比輻射率校正,獲得地表溫度LST[24-25]:

式中,L為TM/TIRS熱紅外波段的輻射值;DN為灰度值;gain和bias為波段的增益與偏置值,可從影像的頭文件中獲得;T為傳感器處溫度值;K1和K2為定標(biāo)參數(shù),TM6波段K1=607.76 W/(m2·sr·μm),K2= 1 260.56 K, TIRS10波 段K1= 774.89 W/(m2·sr·μm),K2=1 321.08 K;λ為TM6波段中心波長(11.5μm)和TIRS10波段的中心波長(10.9 μm);ρ=1.438×10-2mK;ε為地表比輻射率,由NDVI閾值法估算[26]。

2.1.1.4 干 度

干度指標(biāo)NDSI通過裸土指數(shù)SI和建筑指數(shù)IBI合成得到[27-28],公式為

式中,Blue、Green、Red、NIR、SWIR1為藍(lán)、綠、紅、近紅外和短波紅外1波段,分別對應(yīng)TM影像中的第1、2、3、4、5波段和OLI影像中的第2、3、4、5、6波段。

2.1.1.5 亮 度

采用纓帽變換中的亮度分量來表征LSESCI中的亮度指標(biāo),公式為

式中,Blue、Green、Red、NIR、SWIR1、SWIR2為藍(lán)、綠、紅、近紅外、短波紅外1和短波紅外2波段,分別對應(yīng)TM影像中的第1、2、3、4、5、7波段和OLI影像中的第2、3、4、5、6、7波段。對于TM影像,系數(shù)ωi(i=1,2,…,6)取 值 分 別 為 0.204 3、0.415 8、0.552 4、0.574 1、0.312 4、0.230 3;對于OLI影像,系數(shù)ωi(i=1,2,…,6)取值分別為0.302 9、0.278 6、0.473 3、0.559 9、0.508 0、0.187 2。

2.1.2 地表生態(tài)狀況指數(shù)構(gòu)建

如表1所示,RSEI指數(shù)涉及綠度、濕度、熱度和干度4個生態(tài)因子,采用歸一化方式對4個因子去量綱,再利用PCA第1主成分耦合;MRSEI在RSEI基礎(chǔ)上對耦合方式進(jìn)行了改進(jìn),利用PCA前3個主成分進(jìn)行生態(tài)因子耦合;RSEIs的改進(jìn)方式主要表現(xiàn)為用標(biāo)準(zhǔn)化處理方式代替歸一化方式對各生態(tài)因子去量綱;LSESCI指數(shù)采用亮度指標(biāo)代替RSEI指數(shù)中干度指標(biāo),并且耦合方式與前3個指數(shù)不同,通過代數(shù)運(yùn)算來進(jìn)行生態(tài)因子耦合。

表1 不同地表生態(tài)狀況指數(shù)構(gòu)建Table 1 Construction of different land surface ecological status indices

2.2 地表生態(tài)狀況指數(shù)性能評估方法

皮爾森相關(guān)系數(shù)是計算變量之間線性相關(guān)程度的統(tǒng)計量,也是研究多個生態(tài)時間序列之間關(guān)系的主要方法之一,并且被廣泛應(yīng)用于生態(tài)指數(shù)評價[8,15,29]。本研究選用皮爾森相關(guān)系數(shù)計算不同指數(shù)模型與各生態(tài)因子之間的相關(guān)系數(shù),以表征其綜合代表程度和各生態(tài)因子的相關(guān)性,公式為

式中,Xi和Yi分別表示地表生態(tài)狀況指數(shù)及生態(tài)因子各像素值;和分別表示地表生態(tài)狀況指數(shù)及生態(tài)因子各像素值均值。

在各指數(shù)與地形相關(guān)性分析中,以開采區(qū)域?yàn)橹行?創(chuàng)建50、100、150 m緩沖區(qū),分析隨著緩沖區(qū)范圍的變化,其地形和各地表生態(tài)狀況指數(shù)的變化規(guī)律;進(jìn)一步沿開采方向繪制剖面線,分析高程和地表生態(tài)狀況指數(shù)在開采方向上變化規(guī)律。為了更好地反映剖面線上的主要變化趨勢,對根據(jù)像素值生成的曲線進(jìn)行滑動濾波處理,即在可變窗口下對所選取的數(shù)據(jù)求其平均值;再利用求得的平均值作為窗口中心,按照點(diǎn)距向前移動,重復(fù)對數(shù)組中的數(shù)據(jù)求其平均值。濾波器窗口寬度為50,公式為

式中,xi為第i個像素值;yn為由N個像素值得到的第n個像素值;N為滑動濾波器窗口寬度。

以金寶鐵礦區(qū)LSESCI指數(shù)為例,滑動濾波前后對比如圖4所示。

圖4 金寶鐵礦區(qū)LSESCI指數(shù)滑動濾波前后對比Fig.4 Comparison of LSESCI index before and after sliding filtering in Jinbao Iron Mining area

2.3 時空變化特征分析

本研究選擇Sen+Mann-Kendall方法對反演得到的生態(tài)環(huán)境時間序列的變化規(guī)律進(jìn)行分析。Sen[30]趨勢度是一種穩(wěn)健的非參數(shù)計算變化率方法,通過計算時間序列的中值,能夠有效降低噪聲的干擾,但是不能判斷序列趨勢的顯著性。因此需引入Mann-Kendall[31-32]方法來實(shí)現(xiàn)對長時間序列變化趨勢的顯著性檢驗(yàn)。Sen+Mann-Kendall法可以很好地減少噪聲的干擾,對數(shù)據(jù)分布無要求且對異常值不敏感,是時空分析中變化趨勢判斷的重要方法。Sen趨勢度計算公式為

式中,QLSES為時間序列數(shù)據(jù)趨勢的變化量;mean·()為取中值函數(shù);n為研究時間節(jié)點(diǎn)總個數(shù),1＜i＜j＜n;LSESi、LSESj分別為時間序列的地表生態(tài)狀況指數(shù)值。根據(jù)研究區(qū)分布特征,通常當(dāng)QLSES＞0.000 5時地表生態(tài)狀況指數(shù)值呈上升趨勢;當(dāng)≤0.000 5時趨勢穩(wěn)定不變;當(dāng)QLSES＜-0.000 5時地表生態(tài)狀況指數(shù)值呈下降趨勢。

Mann-Kendall趨勢檢驗(yàn)統(tǒng)計量S的計算公式為

當(dāng)n≥10時,S近似服從標(biāo)準(zhǔn)正態(tài)分布,使用檢驗(yàn)統(tǒng)計量Z進(jìn)行趨勢檢驗(yàn):

式中,VAR·()表示方差運(yùn)算。采用顯著性水平α=0.05進(jìn)行顯著性檢驗(yàn),即時,礦區(qū)生態(tài)環(huán)境質(zhì)量呈顯著變化趨勢;時,礦區(qū)生態(tài)環(huán)境質(zhì)量發(fā)生不顯著變化[33]。

3 結(jié)果與討論

3.1 地表生態(tài)狀況指數(shù)適用性分析

3.1.1 地表生態(tài)狀況指數(shù)與各生態(tài)因子之間的相關(guān)性分析

研究區(qū)內(nèi),各地表生態(tài)狀況指數(shù)與生態(tài)因子之間的相關(guān)性分析結(jié)果(表2)均通過顯著性檢驗(yàn)(p＜0.05),平均相關(guān)系數(shù)均達(dá)到0.6以上,能綜合代表各生態(tài)因子信息。其中,綠度和濕度與RSEI、RSEIs、MRSEI呈正相關(guān)關(guān)系,與LSESCI指數(shù)呈負(fù)相關(guān);溫度、亮度及干度與各地表生態(tài)狀況指數(shù)的相關(guān)性關(guān)系則相反。LSESCI指數(shù)與綠度平均相關(guān)性最高,最大負(fù)相關(guān)值為0.87,高于其他指數(shù)。RSEI及MRSEI指數(shù)與溫度平均相關(guān)性最高,最大負(fù)相關(guān)值為0.82,高于LSESCI指數(shù)。RSEIs指數(shù)與干度平均相關(guān)性最高達(dá)0.94,高于RSEI與MRSEI指數(shù)。

表2 地表生態(tài)狀況指數(shù)與生態(tài)因子相關(guān)性分析Table 2 Correlation analysis of land surface ecological status indices and ecological factors

由于LSESCI與RSEI指數(shù)對各生態(tài)因子的耦合方式不同,兩者與各生態(tài)因子相關(guān)性數(shù)據(jù)形態(tài)分布差異較大,2020年金寶鐵礦區(qū)LSESCI指數(shù)與溫度、亮度相關(guān)性分布及RSEI指數(shù)與溫度、干度相關(guān)性分布如圖5所示。由圖5可知:LSESCI指數(shù)與溫度和亮度呈正相關(guān)且相關(guān)性較高,數(shù)值集中在0.65～0.90;RSEI指數(shù)與溫度、干度呈負(fù)相關(guān),數(shù)據(jù)分布均勻且線性規(guī)律明顯,但是相關(guān)性數(shù)值高的樣本點(diǎn)明顯少于LSESCI生態(tài)因子。

圖5 金寶鐵礦區(qū)LSESCI、RSEI指數(shù)與生態(tài)因子相關(guān)性Fig.5 Correlation between LSESCI,RSEI indices and ecological factors in Jinbao Iron Mining Area

2020年黑山煤礦區(qū)不同地表生態(tài)狀況指數(shù)與溫度、干度相關(guān)性分布如圖6所示。由圖6可知:RSEI和MRSEI指數(shù)與各生態(tài)因子相關(guān)性分布差異較大。RSEI指數(shù)與各生態(tài)因子的相關(guān)性較高且數(shù)值主要集中在0.60～0.82,表現(xiàn)在圖中三點(diǎn)分布均勻且密切,但是MRSEI指數(shù)與各生態(tài)因子相關(guān)性之間差異較大,線性規(guī)律不明顯,說明對PCA前3個主成分加和的過程降低了與各生態(tài)因子之間的相關(guān)性。

圖6 黑山煤礦區(qū)RSEI、MRSEI指數(shù)與生態(tài)因子相關(guān)性Fig.6 Correlation between RSEI,MRSEI indices and ecological factors in Heishan Coal Mining Area

2020年烏拉根鋅礦區(qū)不同地表生態(tài)狀況指數(shù)與溫度、干度相關(guān)性分布如圖7所示。由圖7可知:RSEI和RSEIs指數(shù)與各生態(tài)因子相關(guān)性分布差異較大。RSEI指數(shù)與各生態(tài)因子的相關(guān)性數(shù)值主要集中在0.6以上且數(shù)據(jù)分布密集,但是RSEIs指數(shù)與各生態(tài)因子相關(guān)性較低,數(shù)值主要集中在0.5以下,說明對各生態(tài)因子采用標(biāo)準(zhǔn)化代替歸一化的方式去量綱,存在降低地表生態(tài)狀況指數(shù)與各生態(tài)因子之間相關(guān)性的情況。

圖7 烏拉根鋅礦區(qū)RSEI、RSEIs指數(shù)與生態(tài)因子相關(guān)性Fig.7 Correlation between RSEI,RSEIs and Ecological Factors in Wulagen Zinc Mining Area

3.1.2 不同土地利用類型下地表生態(tài)狀況指數(shù)分布

土地覆蓋類型的改變是驅(qū)動地表生態(tài)環(huán)境質(zhì)量變化的重要原因之一[34],研究區(qū)土地利用類型中的開采區(qū)覆蓋了采坑、剝離區(qū)及未復(fù)墾排土場,生態(tài)環(huán)境質(zhì)量處于“差”及“較差”等級。裸土指研究區(qū)周邊原地貌裸露土石、戈壁、荒漠,是干旱半干旱生態(tài)系統(tǒng)中的典型生態(tài)要素,生態(tài)環(huán)境質(zhì)量處于“中等”等級。由于研究區(qū)處于生態(tài)脆弱區(qū),植被稀少、有裸土、基巖、戈壁等,因此,需探究地表生態(tài)狀況指數(shù)區(qū)分不同露天礦開采區(qū)和原地貌戈壁、荒漠的能力,這是客觀評價該區(qū)域露天開采區(qū)對生態(tài)環(huán)境影響的基礎(chǔ)。

如圖8所示,金寶鐵礦區(qū)原地貌主要的土地覆蓋類型為開采區(qū)、裸土、建筑及低覆蓋度植被。該礦區(qū)LSESCI指數(shù)獲得生態(tài)環(huán)境質(zhì)量分布與土地覆蓋類型分布形態(tài)最吻合,開采區(qū)和裸土的平均LSESCI值之差最大為0.295,其次為MRSEI指數(shù),RSEI和RSEIs得到的開采區(qū)和裸土差值小于0.15。因此露天開采周邊低覆蓋植被及裸土等均被劃分為生態(tài)環(huán)境質(zhì)量“差”和“較差”等級,面積分別為7.50 km2、7.16 km2,遠(yuǎn)大于開采區(qū)面積3.97 km2,低估了研究區(qū)的生態(tài)環(huán)境質(zhì)量。

圖8 2020年金寶鐵礦不同地表生態(tài)狀況指數(shù)空間分布Fig.8 Spatial distribution of different land surface ecological status indices in Jinbao Iron Mine in 2020

黑山煤礦區(qū)原地貌主要為戈壁和沙漠,地形相對平坦,土壤類型為荒漠化鹽堿土壤,粗粒含量高。植被分布稀疏,覆蓋率通常為5%～10%。露天開采對該區(qū)域的影響主要為開采區(qū)和排土場占地。RSEI和RSEIs得到不同生態(tài)環(huán)境質(zhì)量分級的空間分布與土地利用數(shù)據(jù)基本吻合(圖9)。由圖9可知:RSEI能夠很好地區(qū)分原地貌為荒漠化鹽堿土壤、植被覆蓋度低土壤與露天開采區(qū),開采區(qū)和原地貌土地利用的RSEI值相差接近0.2,而LSESCI相差值不到0.1,存在將周邊戈壁、荒漠與開采區(qū)、排土場一同評估為生態(tài)環(huán)境質(zhì)量“差”等級的可能性,無法客觀反映露天開采的影響。RSEI指數(shù)不僅獲得的研究區(qū)“較差”和“差”等級面積與開采區(qū)、排土場和剝離區(qū)面積接近,而且能夠反映出排土場區(qū)域生態(tài)環(huán)境質(zhì)量為“良”等級的土地,即排土場生態(tài)恢復(fù)及復(fù)墾情況。

圖9 2020年黑山煤礦不同地表生態(tài)狀況指數(shù)空間分布Fig.9 Spatial distribution of different land surface ecological status indices in Heishan Coal Mine in 2020

烏拉根鋅礦位于荒漠綠洲區(qū),礦區(qū)生態(tài)問題以荒漠化和鹽漬化脅迫為主,土壤為棕漠土,開采對環(huán)境的影響除了土地挖損,還包括廢石堆、尾礦,以及土壤鹽漬化,如圖10所示。RSEI指數(shù)能夠通過生態(tài)環(huán)境質(zhì)量“差”等級反映露天開采對土地的破壞、占用,通過生態(tài)環(huán)境質(zhì)量“較差”等級反映開采區(qū)土壤鹽漬化問題,LSESCI及RSEIs指數(shù)存在低估整體生態(tài)環(huán)境質(zhì)量的情況(圖10(d)、圖10(f))。如圖10(c)所示,RSEI指數(shù)區(qū)分周邊原地貌荒漠化土壤與露天開采區(qū)的能力最強(qiáng),開采區(qū)和原地貌土地利用RSEI差值為0.167,遠(yuǎn)高于其他指數(shù)。

圖10 2020年烏拉根鋅礦不同地表生態(tài)狀況指數(shù)空間分布Fig.10 Spatial distribution of different land surface ecological status indices in Wulagen Zinc Mine in 2020

綜上,LSESCI指數(shù)區(qū)分金寶鐵礦開采區(qū)和裸土的能力最強(qiáng),RSEI指數(shù)區(qū)分黑山煤礦區(qū)及烏拉根鋅礦開采區(qū)和裸土的能力最強(qiáng),能夠客觀反映露天開采對區(qū)域環(huán)境的影響程度。構(gòu)建LSESCI與RSEI指數(shù)的生態(tài)因子中最顯著差異是LSESCI中含有基于纓帽變換獲得的亮度、綠度和濕度組合信息,而RSEI中不含有亮度指標(biāo)。因此,開采區(qū)和裸土的反射率差異是造成兩者生態(tài)環(huán)境質(zhì)量反演結(jié)果不同的重要原因。比較不同露天礦開采區(qū)與裸土亮度均值差異大小發(fā)現(xiàn),金寶鐵礦開采區(qū)與裸土亮度均值之差為0.174,大于黑山煤礦區(qū)和烏拉根鋅礦區(qū)(圖11)?？梢?對于開采區(qū)和周圍自然地物亮度差異小的礦區(qū),RSEI指數(shù)比LSESCI指數(shù)能更好地反演露天開采對環(huán)境的影響。

3.1.3 不同地形條件下地表生態(tài)狀況指數(shù)分布時空特征

通過將研究區(qū)內(nèi)不同地表生態(tài)狀況指數(shù)與高程進(jìn)行疊加,并以開采區(qū)域?yàn)橹行?創(chuàng)建50、100、150 m緩沖區(qū),分析緩沖區(qū)內(nèi)不同地形條件下地表生態(tài)狀況指數(shù)分布特點(diǎn),通過對沿開采方向地表生態(tài)狀況指數(shù)的變化程度進(jìn)行滑動濾波,來揭示地表生態(tài)狀況指數(shù)對露天開采影響程度的表征能力。

金寶鐵礦區(qū)隨著緩沖區(qū)范圍變大,海拔升高,不同地表生態(tài)狀況指數(shù)對生態(tài)環(huán)境質(zhì)量的反演結(jié)果均表現(xiàn)為,距離開采區(qū)遠(yuǎn)區(qū)域的生態(tài)環(huán)境質(zhì)量優(yōu)于距離開采區(qū)近區(qū)域(圖12)。其中,距離開采區(qū)50 m區(qū)域LSESCI值為0.549,高于距離開采區(qū)150 m區(qū)域的LSESCI值(0.506),較其他地表生態(tài)狀況指數(shù)差異最為明顯。沿礦區(qū)開采方向,LSESCI指數(shù)像素波動最為劇烈,且與RSEI指數(shù)像素波動規(guī)律一致。

圖12 金寶鐵礦地表生態(tài)狀況指數(shù)與地形因子疊加分析Fig.12 Superposition analysis of land surface ecological status indices and topographic factors in Jinbao Iron Mine

黑山煤礦區(qū)地形主要呈現(xiàn)西高東低,緩沖區(qū)范圍增加,海拔升高。距離開采區(qū)50 m區(qū)域的RSEI值為0.558,隨著緩沖區(qū)范圍增加,距離開采區(qū)150 m區(qū)域的RSEI值為0.604,可以通過RSEI來明確露天開采的影響范圍。RSEIs和MRSEI指數(shù)隨著緩沖區(qū)范圍的擴(kuò)大差異很小,如圖13所示。通過沿礦區(qū)開采方向地表生態(tài)狀況指數(shù)像素值的分析,MRSEI指數(shù)明顯低于其他指數(shù),低估了生態(tài)環(huán)境質(zhì)量。

圖13 黑山煤礦地表生態(tài)狀況指數(shù)與地形因子疊加分析Fig.13 Superposition analysis of land surface ecological status indices and topographic factors in Heishan Coal Mine

烏拉根鋅礦區(qū)隨著緩沖區(qū)范圍增加,海拔降低。距離開采區(qū)50 m區(qū)域RSEI值為0.608,隨著緩沖區(qū)范圍增加,距離開采區(qū)150 m區(qū)域RSEI值為0.639,較其他地表生態(tài)狀況指數(shù)差異最明顯(圖14)。RSEIs過低估計烏拉根鋅礦區(qū)生態(tài)環(huán)境質(zhì)量,RSEI、MRSEI、LSESCI指數(shù)對生態(tài)環(huán)境質(zhì)量反演結(jié)果均表現(xiàn)為,距離開采區(qū)遠(yuǎn)區(qū)域生態(tài)環(huán)境質(zhì)量優(yōu)于距離開采區(qū)近區(qū)域。其中,RSEI指數(shù)沿礦區(qū)開采方向地表生態(tài)狀況指數(shù)像素值波動最劇烈,能更好地監(jiān)測開采區(qū)周圍土地受采礦活動擾動導(dǎo)致的生態(tài)環(huán)境質(zhì)量變化。

圖14 烏拉根鋅礦地表生態(tài)狀況指數(shù)與地形因子疊加分析Fig.14 Superposition analysis of land surface ecological status indices and topographic factors in Wulagen Zinc Mine

3.2 生態(tài)環(huán)境質(zhì)量時空變化特征動態(tài)監(jiān)測

根據(jù)上述分析,本研究選擇LSESCI指數(shù)分析金寶鐵礦區(qū)生態(tài)環(huán)境質(zhì)量,選擇RSEI指數(shù)分析黑山煤礦區(qū)及烏拉根鋅礦區(qū)生態(tài)環(huán)境質(zhì)量,相關(guān)結(jié)果如圖15所示。

由圖15分析可知:

(1)2005—2020年,金寶鐵礦區(qū)生態(tài)環(huán)境質(zhì)量處于“差”和“較差”等級面積從2.50 km2增加至5.15 km2,處于“中等”等級面積占比由66.57%降至49.12%,呈明顯下降趨勢。如圖15(b1)所示,大部分區(qū)域生態(tài)環(huán)境質(zhì)量為穩(wěn)定不變,生態(tài)質(zhì)量顯著下降區(qū)域主要為采坑和排土場,生態(tài)環(huán)境顯著上升區(qū)域包括礦區(qū)內(nèi)部的排土場復(fù)墾、原地貌草場休牧以及圍欄,其中2020年該礦通過植被恢復(fù)等方式對排土場邊坡治理的面積共計0.25 km2。

圖15 各礦山生態(tài)環(huán)境質(zhì)量時空特征動態(tài)監(jiān)測Fig.15 Dynamic monitoring of spatio-temporal characteristics of ecological environment quality in different mines

(2)2005—2020年,黑山煤礦區(qū)生態(tài)環(huán)境質(zhì)量處于“差”和“較差”等級面積從32.86 km2增加至50.03 km2,主要分布于開采區(qū)及排土場,生態(tài)環(huán)境質(zhì)量等級處于“良”和“中等”面積變化較小,區(qū)內(nèi)復(fù)墾工程以防止水土流失、邊坡穩(wěn)定為主,無植被恢復(fù)工程,外排土場復(fù)墾面積約為1.20 km2。

(3)2005—2020年,烏拉根鋅礦區(qū)生態(tài)環(huán)境質(zhì)量處于“差”和“較差”等級面積變化不大,主要由于開采方式為露天—地下聯(lián)合開采,此外有部分廢石場和尾礦庫。如圖15(b3)所示,復(fù)墾工程主要集中在內(nèi)排土場和外排土場,表現(xiàn)為生態(tài)環(huán)境顯著上升,面積約0.83 km2。

4 結(jié)論與討論

本研究基于濕度、熱度、干度、綠度、亮度生態(tài)因子耦合構(gòu)建了LSESCI、RSEI、MRSEI以及RSEIs指數(shù),結(jié)合土地利用數(shù)據(jù)、DEM數(shù)據(jù),通過相關(guān)性分析、緩沖區(qū)分析、滑動濾波等方法分析了指數(shù)間性能差異,并利用Sen+Mann-Kendall時間序列分析法揭示了露天礦生態(tài)環(huán)境變化趨勢。所得結(jié)論如下:

(1)各地表生態(tài)狀況指數(shù)在干旱半干旱區(qū)露天礦與各生態(tài)因子之間的平均相關(guān)系數(shù)達(dá)0.6以上,RSEI指數(shù)對各生態(tài)因子擬合程度最好。金寶鐵礦區(qū)LSESCI指數(shù)與各生態(tài)因子相關(guān)系數(shù)集中在0.65～0.90,但數(shù)據(jù)分布形態(tài)不均勻。黑山煤礦區(qū)MRSEI指數(shù)與各生態(tài)因子之間相關(guān)性差異較大,數(shù)據(jù)分散,線性規(guī)律不明顯。烏拉根鋅礦區(qū)RSEIs指數(shù)與各生態(tài)因子相關(guān)性較低,數(shù)值主要集中在0.5以下。

(2)通過不同指數(shù)反演獲得的生態(tài)環(huán)境質(zhì)量結(jié)果與土地利用數(shù)據(jù)進(jìn)行疊加分析發(fā)現(xiàn),金寶鐵礦開采區(qū)和裸土的平均LSESCI值之差最大為0.295,LSESCI指數(shù)對于亮度均值差異大的土地覆蓋類型的生態(tài)環(huán)境質(zhì)量反演結(jié)果更準(zhǔn)確。與地形數(shù)據(jù)疊加分析表明,各地表生態(tài)狀況指數(shù)表征不同露天礦開采影響范圍的能力不同,LSESCI指數(shù)在金寶鐵礦區(qū)表現(xiàn)最好,RSEI指數(shù)在黑山煤礦區(qū)和烏拉根鋅礦區(qū)表現(xiàn)最好。

(3)2005—2020年,金寶鐵礦區(qū)由于露天開采導(dǎo)致的生態(tài)環(huán)境質(zhì)量下降區(qū)域面積為2.65 km2,通過排土場復(fù)墾,生態(tài)環(huán)境質(zhì)量上升區(qū)域面積達(dá)到0.55 km2,占下降面積的20.75%。黑山煤礦區(qū)由于露天開采導(dǎo)致的生態(tài)環(huán)境質(zhì)量下降區(qū)域面積為17.17 km2,主要為開采區(qū)和排土場,外排土場復(fù)墾面積約為1.20 km2,占下降面積的6.99%。烏拉根鋅礦區(qū)生態(tài)環(huán)境質(zhì)量下降區(qū)域面積為2.40 km2,復(fù)墾面積約0.83 km2,復(fù)墾區(qū)域生態(tài)環(huán)境質(zhì)量顯著提高,占下降面積的34.58%。

本研究的主要工作是通過綜合考慮礦區(qū)的類型、背景環(huán)境、地形地貌等條件,對典型地表生態(tài)狀況指數(shù)在干旱半干旱露天礦區(qū)生態(tài)環(huán)境監(jiān)測適用性方面取得了新的認(rèn)識。但是,定量遙感監(jiān)測的準(zhǔn)確性受到數(shù)據(jù)分辨率、輻射傳輸過程、模型精度等多個因素的影響,未來還需要開展更多的實(shí)地驗(yàn)證工作,包括光譜測量、地面調(diào)查等,使其生態(tài)因子更具有代表性,并在耦合方法選擇上充分考慮生態(tài)因子的統(tǒng)計學(xué)意義及地學(xué)意義,提高礦區(qū)生態(tài)環(huán)境反演和監(jiān)測的精度及準(zhǔn)確度。