郭山川, 湯 傲, 李效順,2①, 卞正富, 侯湖平, 倪 衡

(1.中國礦業(yè)大學(xué)環(huán)境與測繪學(xué)院, 江蘇 徐州 221116; 2.中國礦業(yè)大學(xué)江蘇省資源環(huán)境信息工程重點實驗室, 江蘇 徐州 221116)

煤炭資源開發(fā)利用對礦區(qū)土地造成持續(xù)性擾動而引發(fā)的植被退化、地表形態(tài)變化、土壤肥力衰退與水土流失等土地?fù)p傷在土地資源承載力脆弱的干旱半干旱礦區(qū)尤為嚴(yán)重。與土地?fù)p毀相比,土地?fù)p傷著重于擾動效應(yīng)引起的土地功能性與結(jié)構(gòu)性衰變,其土地綜合整治方法包含土地復(fù)墾、生態(tài)修復(fù)和自然恢復(fù)等[1-3]。土地?fù)p傷機理主要有:開采沉陷過程導(dǎo)致的土壤豎向遷移運動,驅(qū)使沉陷盆地內(nèi)土壤的空間分布特征重建;開采、運輸、沉陷等過程中重金屬、銨態(tài)氮、速效磷等物質(zhì)發(fā)生轉(zhuǎn)移從而改變土壤化學(xué)特征;土壤微生物多樣性、數(shù)量和種群結(jié)構(gòu)在采礦脅迫下也呈現(xiàn)出顯著響應(yīng)特征。由此可見,礦區(qū)土地在采礦持續(xù)擾動下遭受長期功能性與結(jié)構(gòu)性損傷,導(dǎo)致其植被覆蓋與形態(tài)結(jié)構(gòu)發(fā)生趨勢性變化[3-6]。

傳統(tǒng)土地?fù)p傷監(jiān)測手段不僅費時、費力,且監(jiān)測點呈離散分布,難以進(jìn)行礦區(qū)全覆蓋監(jiān)測,亦不能揭示礦區(qū)土地?fù)p傷空間分布規(guī)律。近年來,國內(nèi)外學(xué)者利用3S技術(shù)和新興的合成孔徑雷達(dá)干涉測量技術(shù)(InSAR)對礦區(qū)土地利用類型、景觀格局、植被覆蓋和地表沉陷等指標(biāo)進(jìn)行了全方位監(jiān)測[7-11]。部分學(xué)者基于上述指標(biāo)分析了采礦擾動下礦區(qū)土地生態(tài)系統(tǒng)的時空演變規(guī)律,進(jìn)而討論了礦區(qū)生態(tài)演替驅(qū)動力、恢復(fù)力以及協(xié)同機制等[12-17]。然而對礦區(qū)土地?fù)p傷的測度研究還相對滯后,尤其是融合主被動遙感提取礦區(qū)土地?fù)p傷因子,建立損傷理論模型以定量分析西北礦區(qū)土地?fù)p傷的研究較為少見。同時礦區(qū)土地?fù)p傷的空間分布特征、時序演替規(guī)律以及土地修復(fù)范式的科學(xué)匹配等問題均需在土地?fù)p傷有效測度的基礎(chǔ)上開展進(jìn)一步研究。

鑒于此,以典型西北干旱半干旱地區(qū)烏海礦區(qū)為例,從采礦擾動過程分析識別關(guān)鍵土地?fù)p傷因子,利用主、被遙感提取土地覆被類型、地表形變以及植被覆蓋度等土地生態(tài)參量,構(gòu)建以植被損傷指數(shù)和地表形變指數(shù)為主控參量的土地?fù)p傷測度模型,實現(xiàn)礦區(qū)土地?fù)p傷的空間全覆蓋監(jiān)測,定量揭示其土地?fù)p傷空間分布規(guī)律,有助于西北礦區(qū)土地修復(fù)范式的精準(zhǔn)化研究。

1 研究區(qū)概況與數(shù)據(jù)介紹

1.1 研究區(qū)概況

烏海市是新興煤炭資源型城市,轄區(qū)內(nèi)礦產(chǎn)資源分布廣泛,因此研究區(qū)烏海礦區(qū)范圍以烏海市行政邊界為參考。烏海市位于內(nèi)蒙古自治區(qū)西南部(39.15°～39.52° N,106.36°～107.05° E),下轄海勃灣區(qū)、烏達(dá)區(qū)和海南區(qū)3個市轄區(qū)(圖1)。地勢呈東西高、中部低,3山成南北走向平行排列,中間形成2條平坦的谷地。地貌為構(gòu)造侵蝕中低山地、剝蝕丘陵區(qū)、山前堆積沖洪積扇區(qū)和黃河沖積堆積階地。烏海市屬于典型干旱半干旱地區(qū),年均溫9.6 ℃,年均降水量159.8 mm,年均蒸發(fā)量3 289 mm。區(qū)域內(nèi)累計查明煤炭資源儲量為28.8億t,是國家重要焦煤基地,持續(xù)大規(guī)模的煤炭資源開發(fā)給國家和地區(qū)帶來巨大的經(jīng)濟(jì)效益,同時也對區(qū)域生態(tài)環(huán)境帶來了負(fù)面影響。

圖1 烏海市示意圖Fig.1 The location of Wuhai City

1.2 數(shù)據(jù)準(zhǔn)備

采用數(shù)據(jù)主要包括遙感數(shù)據(jù)、合成孔徑雷達(dá)(SAR)數(shù)據(jù)、數(shù)字高程模型和行政區(qū)劃數(shù)據(jù)等(表1)。

表1主要數(shù)據(jù)來源

Table1Datasources

數(shù)據(jù) 采集時間數(shù)據(jù)來源備注 Landsat TM2000年8月22日美國地質(zhì)勘探局(USGS)云量0.34%,空間分辨率30 m Landsat TM2005年6月1日USGS云量0.12%,空間分辨率30 m Landsat TM2009年6月28日USGS云量0,空間分辨率30 m Landsat OLI2015年9月1日USGS云量0.08%,空間分辨率30 m 行政區(qū)劃圖2012年地方機構(gòu)比例尺1∶10萬 地形圖2012年地方機構(gòu)比例尺1∶10萬 礦區(qū)分布圖2012年地方機構(gòu)比例尺1∶10萬 礦山環(huán)境治理規(guī)劃評估2012年地方機構(gòu)比例尺1∶10萬 礦山地質(zhì)環(huán)境影響評價圖2012年地方機構(gòu)比例尺1∶10萬 Sentinel-1A雷達(dá)數(shù)據(jù)2015-01-13—2015-12-27(共9景)歐洲航天局(ESA)空間分辨率3.7×13.9 m SRTM數(shù)字高程模型2003年美國航空航天局(NASA)空間分辨率30 m

USGS、ESA、NASA數(shù)據(jù)來源網(wǎng)站分別為https:∥glovis.usgs.gov/、https:∥scihub.copernicus.eu/和http:∥srtm.csi.cgiar.org/。

由于季節(jié)及時間累積效應(yīng)差異等原因易使遙感影像信息產(chǎn)生偽變化,因此應(yīng)盡量采用同時相、同時間間隔的遙感影像作為數(shù)據(jù)源進(jìn)行處理分析。通過輻射校正、投影變化、幾何校正、數(shù)據(jù)融合和影像剪裁等方法對研究區(qū)2000年8月、2005年6月、2009年6月(2010年夏季遙感影像缺失)的Landsat TM和2015年9月的Landsat OLI遙感影像進(jìn)行預(yù)處理。SAR影像為2015年1月13日至2015年12月27日共9景3.7×13.9 m分辨率、干涉寬幅升軌模式(IW)的Sentinel-1A數(shù)據(jù)。

2 測度模型及方法

采礦活動對礦區(qū)土地的擾動形式主要有挖損、壓占、沉陷、地裂縫和滑坡等，這些擾動過程的持續(xù)擠壓式作用導(dǎo)致土地功能與結(jié)構(gòu)損傷,反映在地表植被覆蓋退化、土地覆被改變和地表形變等方面。因此,通過3S技術(shù)提取植被覆蓋與土地覆被信息,InSAR技術(shù)反演地表形變,并構(gòu)建損傷測度模型對礦區(qū)土地?fù)p傷進(jìn)行定量測度和評估[18-19]。

2.1 土地?fù)p傷測度模型

以地表形變指數(shù)和植被損傷指數(shù)為主控變量,同時引入土地覆被類型變量,從空間像元級尺度構(gòu)建西北礦區(qū)土地?fù)p傷測度模型,空間像元p土地?fù)p傷DL,p的測度模型為

DL,p=g(IND,p,TLC,p)+f(INVD,p,TLC,p)。

(1)

式(1)中,g(IND，p，TLC，p)為由像元p地表形變指數(shù)(normalized deformation index,NDI，IND)和土地覆被類型(land cover type,LCT，TLC)決定的土地結(jié)構(gòu)模塊損傷函數(shù);f(INVD，p，TLC，p)為由像元p植被損傷指數(shù)(normalized vegetation damage index,NVDI，INVD)和LCT決定的土地功能模塊損傷函數(shù),IND，p，INVD，p∈[0，1]。

礦區(qū)土地?fù)p傷在持續(xù)擾動過程中存在累積效應(yīng)和恢復(fù)力作用,在兩者的耦合作用下,損傷測度邊際函數(shù)呈先增大后減小特征,利用該特性建立形變變量和植被損傷變量與土地?fù)p傷的微分方程組,對微分方程組積分并進(jìn)行歐拉變換得:

(2)

(3)

式(2)～(3)中,α、β為與TLC，p相關(guān)的權(quán)重系數(shù),表示不同土地覆被類型對礦區(qū)土地?fù)p傷的敏感程度。結(jié)合烏海礦區(qū)實地調(diào)研資料和TD/T 1010—2015《土地利用動態(tài)遙感監(jiān)測規(guī)程》[20],將研究區(qū)分為林草地、耕地、水域、裸地、建設(shè)用地和采煤用地6種土地覆被類型,利用德爾菲法對不同土地覆被類型的土地?fù)p傷敏感程度進(jìn)行排序分析(表2)。

表2不同土地覆被類型的土地?fù)p傷敏感性分析

Table2Thesensitivityoflandusetypestolanddamageindex

土地覆被類型 土地?fù)p傷結(jié)構(gòu)損傷函數(shù)系數(shù)(α)功能損傷函數(shù)系數(shù)(β)林草地DB耕地CA水域FF裸地EC建設(shè)用地AD采煤用地BE

A～F的分值賦為5～0。

如圖2所示,測度模型輸入?yún)⒘繛楸磉_(dá)土地結(jié)構(gòu)損傷的形變參量,表達(dá)土地功能損傷的植被損傷參量,以及由土地覆被類型決定的系數(shù)α、β。礦區(qū)損傷測度模型由結(jié)構(gòu)損傷測度和功能損傷測度函數(shù)構(gòu)成,輸出結(jié)果為礦區(qū)內(nèi)像元p的土地?fù)p傷值。

圖2 烏海礦區(qū)土地?fù)p傷測度框架

2.2 損傷參量提取方法

利用遙感方法計算礦區(qū)歸一化差分植被指數(shù)(normalized difference vegetation index, NDVI,INDV)以獲取土地?fù)p傷測度模型中的植被參量,NDVI是遙感影像中近紅外波段的反射值和紅光波段的反射值之差與上述兩者之和的比值[21],其有效值域為[0,1],對地面覆蓋為云、水、雪等負(fù)NDVI值區(qū)域進(jìn)行零值處理。礦區(qū)下沉盆地土地植被覆蓋變化能一定程度上反映其土地生產(chǎn)功能損傷情況[22-23],因此植被損傷指數(shù)NVDI與NDVI函數(shù)關(guān)系可表達(dá)為INVD=1-INDV。

礦區(qū)地表形變參量信息提取采用小基線集合成孔徑雷達(dá)干涉測量(small baseline subset interferometric synthetic aperture radar, SBAS-InSAR)技術(shù),該技術(shù)以全天候、全天時、高分辨率和連續(xù)空間覆蓋的優(yōu)勢在地表形變測量中得到廣泛關(guān)注[24-25]。SBAS-InSAR作為重要時序InSAR技術(shù)之一,能夠有效減小一般InSAR易受失相干、大氣延遲效應(yīng)等因素引起的誤差,在礦區(qū)連續(xù)、劇烈形變場反演中具有獨一無二的優(yōu)勢[26-27]。對SBAS-InSAR提取的地表形變速率(ground deformation velocity, GDV,VGD)進(jìn)行無量綱歸一化處理,得到地表形變指數(shù)(IND):

(4)

另外,采用支持向量機算法(support vector machine, SVM)對烏海礦區(qū)遙感影像的光譜和紋理融合特征進(jìn)行監(jiān)督分類學(xué)習(xí),提取烏海礦區(qū)土地覆被類型。

3 結(jié)果與分析

3.1 植被提取結(jié)果

烏海礦區(qū)2000、2005、2009和2015年4個時相的植被覆蓋分布如圖3所示。從空間分布上看,整體上低植被覆蓋水平反映了烏海礦區(qū)突出的生態(tài)脆弱性,而且東南區(qū)域植被覆蓋水平顯著低于其他區(qū)域。

圖3 2000—2015年烏海礦區(qū)NDVI的空間分布

為進(jìn)一步分析烏海礦區(qū)植被覆蓋信息提取結(jié)果的時序變化特征,計算烏海礦區(qū)范圍內(nèi)每個地理單元的平均NDVI值,結(jié)果見圖4。可以發(fā)現(xiàn),各期NDVI值均小于0.05,說明烏海礦區(qū)植被覆蓋整體水平低;2000—2009年烏海礦區(qū)植被覆蓋呈現(xiàn)逐期衰退趨勢,期間植被覆蓋水平下降22.9%,標(biāo)準(zhǔn)差呈遞減趨勢;而2009—2015年烏海礦區(qū)植被覆蓋出現(xiàn)恢復(fù)反彈,反彈力度達(dá)10.8%,主要原因是地方政府和煤炭企業(yè)落實了生態(tài)文明建設(shè)政策方針[28-30],開展了系列土地整治、植樹造林和環(huán)境治理等工程,同期標(biāo)準(zhǔn)差值激增印證了研究區(qū)植被覆蓋恢復(fù)是土地綜合整治工程造成的空間離散式恢復(fù)。

統(tǒng)計烏海礦區(qū)不同植被覆蓋水平面積(表3),可以發(fā)現(xiàn),2000—2005年,NDVI水平為>0.6～1的區(qū)域面積減小45.8 km2,>0.3～0.6的區(qū)域面積減小12.5 km2,>0.1～0.3的區(qū)域面積增加86.2 km2;2005—2009年,NDVI水平為>0.3～1的區(qū)域面積略有下降,共減小1.5 km2;2009—2015年,NDVI水平為>0.6～1的區(qū)域面積增加52.6 km2,>0.3～0.6的區(qū)域增加21.1 km2,生態(tài)治理效果明顯。

3.2 沉陷反演結(jié)果

以時相為2015年9月22日和2015年11月9日的2景SAR圖像差分干涉圖為例(圖5)展示SBAS-InSAR干涉結(jié)果,可以發(fā)現(xiàn),干涉圖紋清晰說明大氣延遲、斑點噪聲等誤差相位通過空間低通和時間高通濾波得到有效削減,差分干涉結(jié)果質(zhì)量較好;干涉圖中呈現(xiàn)出系列沉陷漏斗,沉陷漏斗的位置、幾何形狀與礦區(qū)下沉盆地特征相似,說明區(qū)域在觀測期間的地表沉陷主要由采礦活動引起。

圖4 2000—2015年烏海礦區(qū)NDVI平均值變化

表32000—2015年烏海礦區(qū)不同NDVI級別的面積與比例

Table3TheareaandproportionofdifferentNDVIlevelsinthestudyareafrom2000to2015

NDVI值2000年2005年2009年2015年比例/%面積/km2比例/%面積/km2比例/%面積/km2比例/%面積/km2 0～0.166.31 163.364.71 135.447.0824.361.11 071.9 >0.1～0.324.3426.929.3513.147.1825.828.8504.4 >0.3～0.66.3110.55.698.05.799.46.9120.5 >0.6～13.053.30.47.40.34.63.357.1 總計1001 7541001 7541001 7541001 754

烏海礦區(qū)SBAS-InSAR地表形變監(jiān)測結(jié)果如圖6所示。從采煤沉陷地地表形變場空間分布來看,地表沉降場與礦區(qū)地理位置有強空間相關(guān)性,沉降場主要分布在煤炭資源開采區(qū),如烏海東部和中部山區(qū),以及黃河以西的烏達(dá)礦區(qū),而其余地方一般沒有沉陷;采煤工作面地理位置多分布于地形起伏多變、地質(zhì)條件復(fù)雜的山區(qū),導(dǎo)致開采下沉盆地等沉線表現(xiàn)形式為系列同心不規(guī)則橢圓群。

利用GIS統(tǒng)計烏海礦區(qū)采煤沉陷監(jiān)測結(jié)果,發(fā)現(xiàn)烏海礦區(qū)開采沉陷地面積為230.03 km2,占烏海礦區(qū)總面積的13.78%;年沉降速率為>200～800 mm·a-1的區(qū)域面積為44.98 km2,占沉陷地面積的19.55%;年沉降速率為>50～200 mm·a-1的區(qū)域面積為101.33 km2,占沉陷地面積的44.05%;年沉降速率為10～50 mm·a-1的區(qū)域面積為83.72 km2,占沉陷地面積的36.40%。

3.3 土地?fù)p傷測度

由圖6可見,烏海礦區(qū)土地?fù)p傷空間分布呈現(xiàn)出以高強度開采工作面為中心向四周輻射減弱的特征:土地?fù)p傷值為>4.5～8.0的高損傷區(qū)地理空間位置主要集中在東部山區(qū)、中部山區(qū)以及黃河以西等煤炭資源開采強度大的區(qū)域;土地?fù)p傷值為>1.0～4.5的中度損傷區(qū)圍繞在高損傷區(qū)周邊;土地?fù)p傷值為0～1.0的非損傷區(qū)由于擴散衰減顯著,其地理分布與礦區(qū)的空間相關(guān)性較弱。研究區(qū)共有586.57 km2土地在采礦擾動下遭到破壞,其中土地?fù)p傷值為>4.5～8.0的高損傷區(qū)面積為42.36 km2,占總面積的2.68%;中損傷區(qū)面積為544.21 km2,占比34.43%;研究區(qū)土地非損傷區(qū)面積為994.02 km2,占比62.89%。

4 討論與結(jié)論

筆者對采礦擾動形式、礦區(qū)土地響應(yīng)特征進(jìn)行分析,同時針對土地?fù)p傷過程中其結(jié)構(gòu)與功能的變化機制,構(gòu)建西北礦區(qū)土地?fù)p傷測度模型,融合主被動遙感提取模型關(guān)鍵參量,以烏海礦區(qū)為例獲取了該區(qū)域的土地?fù)p傷值,分析了其空間分布規(guī)律。以烏海市礦山環(huán)境治理規(guī)劃評估圖、礦山地質(zhì)環(huán)境影響評價圖、地面實測數(shù)據(jù)以及各礦區(qū)提供的礦山地質(zhì)環(huán)境治理數(shù)據(jù)為驗證數(shù)據(jù)對比分析筆者提出的土地?fù)p傷測度結(jié)果,發(fā)現(xiàn)測度結(jié)果與實際資料在空間上具有較強相關(guān)性,而且在時序上的關(guān)聯(lián)性也得到了很好的檢驗。

(1)采礦擾動下,烏海礦區(qū)在2000—2009年間植被覆蓋逐漸萎縮,2009—2015年期間植被覆蓋水平提高是由于地方進(jìn)行了局部區(qū)域的土地整治與修復(fù)工程,其標(biāo)準(zhǔn)差激增說明植被恢復(fù)過程是空間不連續(xù)的離散式事件。

圖5 烏海礦區(qū)局部差分干涉圖

(2)采礦擾動下,烏海礦區(qū)土地結(jié)構(gòu)破壞表現(xiàn)在區(qū)域內(nèi)出現(xiàn)大范圍地表沉陷場,其沉陷速率最大超800 mm·a-1,監(jiān)測到的沉陷漏斗與煤炭資源開采區(qū)有較強的空間相關(guān)性,其下沉等值線形式特征為系列不規(guī)則橢圓群。

(3)通過西北礦區(qū)土地?fù)p傷模型測算,發(fā)現(xiàn)土地?fù)p傷在空間分布上出現(xiàn)以開采工作面為中心向周邊擴散衰減的現(xiàn)象,說明采礦擾動是導(dǎo)致烏海礦區(qū)土地功能性和結(jié)構(gòu)性損傷的主因之一,另一方面也反映出烏海礦區(qū)土地修復(fù)工作刻不容緩。

(4)烏海礦區(qū)土地?fù)p傷測度結(jié)果與實際調(diào)研資料相比,具有較一致的空間相關(guān)性和時序關(guān)聯(lián)性,印證了該模型的有效性。土地?fù)p傷測度模型有助于西北礦區(qū)長時序、全覆蓋土地?fù)p傷監(jiān)測體系的搭建、土地?fù)p傷精準(zhǔn)修復(fù)策略的研究以及損傷-修復(fù)動態(tài)反饋機制的建立等,有利于“青山”與“金山”協(xié)同發(fā)展。

圖6 烏海礦區(qū)地表沉降速率、土地覆被和土地?fù)p傷值示意

筆者利用提出的土地?fù)p傷測度模型能夠有效應(yīng)用于烏海礦區(qū)土地?fù)p傷調(diào)查,對采煤塌陷地治理工作有一定的現(xiàn)實意義,但也存在不足之處:土地?fù)p傷提取結(jié)果缺乏大量實測數(shù)據(jù)進(jìn)行定量驗證;提供SAR影像的Sentinel-1A衛(wèi)星于2014年發(fā)射運行,形變監(jiān)測時相與植被監(jiān)測時相無法向前吻合,導(dǎo)致烏海礦區(qū)土地?fù)p傷監(jiān)測時相較為單一。后續(xù)的研究重點是針對上述2點問題展開進(jìn)一步分析,以更有力地指導(dǎo)礦區(qū)土地整治與生態(tài)修復(fù)工作。