基于時(shí)序TM影像的淮南潘謝礦區(qū)采煤沉陷水域動(dòng)態(tài)監(jiān)測(cè)

王啟元 肖 武,2 李素萃 張文凱 張 偉

(1.中國(guó)礦業(yè)大學(xué)(北京)土地復(fù)墾與生態(tài)重建研究所，北京市海淀區(qū)，100083；2.浙江大學(xué)公共管理學(xué)院，浙江省杭州市，310058)

摘 要利用時(shí)序TM遙感影像數(shù)據(jù)，在對(duì)比單波段閾值法、譜間關(guān)系法、歸一化水體指數(shù)法和改進(jìn)歸一化水體指數(shù)法4種水體信息提取精度的基礎(chǔ)上，選擇改進(jìn)歸一化水體指數(shù)(MNDWI)法并結(jié)合GIS空間疊加技術(shù)，對(duì)潘謝礦區(qū)2005-2017年采沉陷水域進(jìn)行提取并分析其時(shí)空演變特征。結(jié)果顯示，2005-2017年，淮南潘謝礦區(qū)采煤沉陷水域面積增至4948.65 hm2，總增量達(dá)到4003.52 hm2，年均增長(zhǎng)333.63 hm2，年均增長(zhǎng)率為20.1%；分析研究區(qū)內(nèi)煤炭開采作為主要貢獻(xiàn)因素對(duì)沉陷水域面積變化影響及預(yù)測(cè)。通過線性擬合分析，預(yù)計(jì)到2020年潘謝礦區(qū)沉陷水域面積5968.12 hm2；2021-2030年沉陷水域面積年均增約為350 hm2，至2030其總量將達(dá)到9468.12 hm2。利用時(shí)序TM遙感影像對(duì)高潛水位礦區(qū)沉陷水體的時(shí)序動(dòng)態(tài)變化監(jiān)測(cè)，分析了其變化原因，為后期采煤沉陷區(qū)的土地復(fù)墾與治理提供了科學(xué)依據(jù)。

東部高潛水位礦區(qū)煤炭開采引起的地表沉陷積水，嚴(yán)重影響了當(dāng)?shù)氐霓r(nóng)業(yè)發(fā)展和居民生活。國(guó)內(nèi)礦區(qū)生態(tài)環(huán)境方面的研究大多集中在生態(tài)環(huán)境修復(fù)、沉降監(jiān)測(cè)以及環(huán)境影響評(píng)估等方面，對(duì)于礦區(qū)采煤沉陷積水的研究主要集中在沉陷積水動(dòng)態(tài)監(jiān)測(cè)以及范圍確定等方面?？焖偾揖珳?zhǔn)獲取沉陷水域時(shí)空變化信息是高潛水位煤糧復(fù)合區(qū)進(jìn)行環(huán)境整治與沉陷地復(fù)墾工作的前提。與傳統(tǒng)手段相比較，遙感監(jiān)測(cè)方法具有快速、精確、宏觀性強(qiáng)和低費(fèi)用等特點(diǎn)，在沉陷水域的長(zhǎng)時(shí)序動(dòng)態(tài)監(jiān)測(cè)上具有無可比擬的優(yōu)勢(shì)。

國(guó)內(nèi)外常見的水體信息提取方法有單波段閾值法、水體指數(shù)法、多波段譜間關(guān)系法等。Frazier等針對(duì)水體在近紅外波段表現(xiàn)為強(qiáng)吸收性而植被和干土壤在該波段具有強(qiáng)反射性特點(diǎn)，提出單波段閾值法對(duì)水體進(jìn)行提取。Kloiber等用TM影像的2波段加上3波段大于4波段加上5波段，提取水體信息。McFeeters提出歸一化水體指數(shù)(NDWI)，該指數(shù)提取水體時(shí)植被和土壤信息能夠被抑制。徐涵秋提出了改進(jìn)歸一化水體指數(shù)(MNDWI)來提取水體信息。彭蘇萍等基于多時(shí)相TM影像，提取了淮南礦區(qū)積水沉陷面積的擴(kuò)展變化信息；張瑞婭等將采區(qū)的原始地形歸入地表沉陷分析，確定了采煤沉陷的積水區(qū)域范圍；楊光華等利用高分遙感影像，對(duì)濟(jì)寧市采煤導(dǎo)致的塌陷積水耕地信息用人機(jī)交互解譯的方式進(jìn)行了提取測(cè)算；孟磊等用TM數(shù)據(jù)分析淮南潘謝礦區(qū)1976-2006年的水體變化情況并分析了其演變的景觀生態(tài)效應(yīng)。利用遙感技術(shù)手段對(duì)水體信息進(jìn)行提取，必須考慮到在光譜特征方面水體與其他地物的差異性，針對(duì)不同的研究區(qū)域的不同自然地理?xiàng)l件，在選擇水體提取方法時(shí)亦有所差異。本文以淮南潘謝礦區(qū)為研究區(qū)，比較4種遙感水體信息提取方法，結(jié)合礦區(qū)塌陷范圍實(shí)測(cè)資料及GIS空間疊加分析技術(shù)，實(shí)現(xiàn)研究區(qū)2005-2017年沉陷水域提取，揭示了潘謝礦區(qū)近十多年沉陷水域的時(shí)空演變特征極，并分析煤炭開采作為主要原因?qū)ζ浔O(jiān)測(cè)時(shí)段內(nèi)的影響及趨勢(shì)預(yù)測(cè)，對(duì)未來采煤沉陷水域變化趨勢(shì)預(yù)測(cè)以及生態(tài)環(huán)境治理具有現(xiàn)實(shí)意義。

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

1.1 研究區(qū)概況

潘謝礦區(qū)位于安徽省淮南市西北部，地形平坦開闊，地處淮河沖積平原，整體地勢(shì)西北高、東南低，地表高程一般為21.5～25.4 m，平均標(biāo)高在23.1 m左右。所在區(qū)域?qū)偌撅L(fēng)溫暖帶半濕潤(rùn)氣候，季節(jié)性表現(xiàn)為夏季炎熱，冬季寒冷。6～8月降雨量約占全年的40%，年平均降雨量865.5 mm，年平均蒸發(fā)量1610.44 mm，蒸發(fā)量大于降雨量，潮濕系數(shù)近似0.5。區(qū)內(nèi)溝渠縱橫，水利設(shè)施完善，當(dāng)?shù)貪撍宦裆罴s為1.5 m。因此地下煤炭開采沉陷后，地表極易形成積水，如圖1所示。截至2014年，淮南礦區(qū)累計(jì)原煤產(chǎn)量約6.7億t，根據(jù)生產(chǎn)規(guī)劃，到2020年，礦區(qū)累計(jì)原煤產(chǎn)量約11億t；2021-2030年，原煤產(chǎn)量穩(wěn)定在7000萬(wàn)t/a。

圖1 淮南潘榭礦區(qū)遙感影像

1.2 數(shù)據(jù)選擇及預(yù)處理

為了研究不同水體提取方法在本區(qū)域內(nèi)的精度與適用性，選取研究區(qū)內(nèi)丁集礦、顧橋礦和顧北礦作為實(shí)驗(yàn)區(qū)域?qū)Ρ确治?種水體提取方法中的最優(yōu)法，選用2016年9月2日Landsat8 OLI影像作為數(shù)據(jù)源，然后整體提取潘謝礦區(qū)6期沉陷水域信息，見表1。所有影像來源于地理空間數(shù)據(jù)云，空間分辨率為30 m，寬幅185 km×185 km，影像包含紅、綠、藍(lán)、近紅外和短波紅外等9個(gè)波段。由于研究所采用的遙感影像獲取時(shí)均已經(jīng)過系統(tǒng)的輻射校正與幾何校正等相關(guān)預(yù)處理工作，本文對(duì)數(shù)據(jù)的預(yù)處理工作主要包括輻射定標(biāo)、大氣校正以及影像裁剪等。2005-2017年6期遙感影像日期均在5月份左右，未到研究區(qū)多降水季節(jié)，降水量對(duì)于水體提取結(jié)果的影響不予考慮。

表1 研究區(qū)遙感數(shù)據(jù)基本情況

2 研究方法

基于TM影像的水體信息自動(dòng)提取方法，考慮影像光譜特性的提取方法目前應(yīng)用最多，其中單波段閾值法、多波段譜間關(guān)系法、水體指數(shù)法這3種方法最為廣泛。本文使用4種水體信息提取方法，選擇研究區(qū)內(nèi)丁集、顧橋、顧北3個(gè)礦作為實(shí)驗(yàn)區(qū)，對(duì)實(shí)驗(yàn)區(qū)內(nèi)所有水體進(jìn)行提取，經(jīng)過精度評(píng)價(jià)找到適合研究區(qū)的最優(yōu)水體信息提取方法，再利用優(yōu)選后的方法對(duì)研究區(qū)內(nèi)2005年、2008年、2009年、2013年、2015年、2017年共6期水體進(jìn)行分階段提取，分析煤炭開采對(duì)沉陷水域變化的影響以及趨勢(shì)預(yù)測(cè)。

2.1 TM影像水體信息提取方法

目前利用遙感影像進(jìn)行水體信息提取的方法主要有單波段(近紅外)閾值法、多波段譜間關(guān)系法和水體指數(shù)法等。單一的水體信息提取方法無法通用，需要根據(jù)研究數(shù)據(jù)情況以及區(qū)域地理?xiàng)l件擇優(yōu)選擇。

2.1.1 單波段閾值法

單波段閾值法利用近紅外波段水體的強(qiáng)吸收性以及植被和干土壤的強(qiáng)反射性特點(diǎn)，選取適當(dāng)?shù)拈撝祵⑺w信息提取出來。此方法原理簡(jiǎn)單，但水體與非水體之間存在的過渡區(qū)域容易被忽略，同時(shí)在區(qū)分細(xì)小水體和混淆在水體中的陰影上很困難，從而產(chǎn)生誤提、漏提的現(xiàn)象。利用水體信息和其他地物信息在遙感圖像不同波段上表現(xiàn)出來的不同亮度值來設(shè)定閾值，提取水體信息。依據(jù)實(shí)驗(yàn)，水體在TM5(近紅外)波段具備很強(qiáng)的吸收作用，亮度值表現(xiàn)較低，而在該波段其他地物亮度值和水體有明顯區(qū)別。使用ENVI5.2軟件，對(duì)2016年9月2日Landsat8 OLI影像近紅外波段灰度圖像使用密度分割法確定閾值，設(shè)定水體提取閾值，利用波段運(yùn)算工具，提取水體信息。

2.1.2 多波段譜間關(guān)系法

波譜間關(guān)系法可根據(jù)水體與地物波譜曲線的特征，找出其中的變化規(guī)律并設(shè)定邏輯判別式提取水體，這樣能夠?qū)⑺w與陰影較好地區(qū)分出來，而且提取精度較高。水體具備TM2(藍(lán))加TM3(綠)大于TM4(紅)加TM5(近紅外)的特征，經(jīng)過試驗(yàn)發(fā)現(xiàn)，設(shè)定一定的閾值來得到提取水體的效果更好，設(shè)定模型如下：

(1)

通過多次反復(fù)實(shí)驗(yàn)，選取K1及K2值，使用波段運(yùn)算工具，得到水體信息提取圖。

2.1.3 歸一化水體指數(shù)法

水體的反射率從可見光波段到短紅外波段逐級(jí)減弱，尤其在近紅外和短波熱紅外范圍內(nèi)顯示強(qiáng)烈的吸收性，因此可采用水體在可見光波段和近紅外波段(NIR)的光譜差異構(gòu)建水體指數(shù)，而植被在近紅外波段的反射率最強(qiáng)，采用綠波段(GREEN)和近紅外波段構(gòu)建指數(shù)可以有效地抑制植被和土壤信息，NDWI指數(shù)法如下：

(2)

同樣對(duì)于NDWI指數(shù)設(shè)定閾值上下限1，K2〗，可以更好地獲得水體信息提取效果。

2.1.4 改進(jìn)歸一化水體指數(shù)法

土壤、建筑物因素的干擾在NDWI指數(shù)中被忽略，而水體綠光和近紅外的波段波譜特征與土壤、建筑物幾乎一致。徐秋涵等利用中紅外波段(MIR)，提出MNDWI指數(shù)用于水體信息的提取，計(jì)算公式如下：

(3)

2.2 遙感影像與采煤沉陷邊界空間疊加

采用水體信息提取方法中精度最高的改進(jìn)歸一化水體指數(shù)法分別對(duì)研究區(qū)內(nèi)2005-2017年(2005年5月8日、2008年5月15日、2009年6月3日、2013年5月21日、2015年5月19日、2017年4月30日)共6期遙感影像分別進(jìn)行體提取，在Arcgis軟件中，使用空間疊加分析模塊，將水體提取結(jié)果與淮南潘榭礦區(qū)2015年地表沉降范圍疊加分析，實(shí)現(xiàn)對(duì)潘榭礦區(qū)采煤沉陷水域的提取。

3 結(jié)果與分析

3.1 水體信息提取方法比較及精度分析

經(jīng)過前期預(yù)處理，對(duì)實(shí)驗(yàn)區(qū)OLI遙感影像應(yīng)用各水體信息自動(dòng)提取方法，得到提取結(jié)果，見表2。

表2 不同方法水體信息提取結(jié)果對(duì)比(黑色為水體)

表2從左到右依次為單波段閾值法、多波段譜間關(guān)系法、NDWI指數(shù)法和MNDWI水體指數(shù)的水體信息提取結(jié)果，黑色部分為水體信息。結(jié)合原始影像目視解譯分析上述提取結(jié)果，MNDWI指數(shù)法對(duì)水體信息的提取效果最好，對(duì)水體斑塊完整性及連續(xù)性提取上表現(xiàn)出色，并且誤提現(xiàn)象最少；單波段閾值法對(duì)水體的提取效果次之，雖然也存在誤提現(xiàn)象，但在水體信息連續(xù)性提取上表現(xiàn)出色，相比與其他兩種自動(dòng)提取方法效果好；水體指數(shù)法效果和譜間關(guān)系法相對(duì)較差，水體信息邊界提取不連貫，存在較多誤提現(xiàn)象。

使用ENVI5.3軟件對(duì)實(shí)驗(yàn)區(qū)遙感原始影像目視解譯，均勻選擇100組樣本，包含像元5562個(gè)。經(jīng)統(tǒng)計(jì)得到水體提取混淆矩陣和Kappa統(tǒng)計(jì)，見表3。由表3可知，4種水體提取方法精度均達(dá)到90%以上，但MNDWI法提取研究區(qū)水體信息的用戶精度、總體精度以及Kappa相關(guān)系數(shù)為各方法中最高。

表3 不同水體方法提取結(jié)果精度評(píng)價(jià)

3.2 歷年水體信息提取

基于MNDWI指數(shù)法，用2005年、2008年、2009年、2013年、2015年和2017年TM/OLI影像潘謝礦區(qū)整體水域進(jìn)行提取，結(jié)果如圖2所示。

由圖2可以看到，2005-2017年，潘謝礦區(qū)地表全域水體擴(kuò)張明顯，尤其在2010年之后增速加快，區(qū)內(nèi)水體斑塊總量不斷增加，個(gè)體積水斑塊多數(shù)由小變大，且斑塊形狀變化巨大。

圖2 監(jiān)測(cè)年份潘謝礦區(qū)地表水域分布

3.3 采煤沉陷水域監(jiān)測(cè)

由于本文主要側(cè)重于分析開采沉陷導(dǎo)致的水體變化，而礦區(qū)內(nèi)仍存在大面積的自然水體，為了區(qū)分未受擾動(dòng)的天然水體與開采沉陷水體，實(shí)現(xiàn)對(duì)研究區(qū)內(nèi)采煤沉陷水域的提取監(jiān)測(cè)，將淮南潘謝礦區(qū)2015年采煤沉陷邊界與水域提取結(jié)果結(jié)合，利用Arcgis10.3軟件進(jìn)行空間疊加運(yùn)算，獲得2005-2017年共6期潘謝礦區(qū)沉陷區(qū)范圍內(nèi)的水域結(jié)果，如圖3所示，并在SPSS軟件中統(tǒng)計(jì)分析潘謝礦區(qū)監(jiān)測(cè)時(shí)段的沉陷水域面積變化情況，見表4。

分析結(jié)果表明，淮南潘謝礦區(qū)沉陷水域面積由2005年的945.13 hm2增加至2017年的4948.65 hm2，累積增長(zhǎng)量為4003.52 hm2，年均增長(zhǎng)333.63 hm2，年均增長(zhǎng)率為20.1%。沉陷水體總量明顯呈上升趨勢(shì)。其中，2005-2009年沉陷水域面積增加1459.3 hm2，增幅達(dá)154%，變化幅度相對(duì)最大；2009-2013年沉陷水域面積增加603.32 hm2，增幅為25.1%；2013-2017年沉陷水域面積增量為1940.9 hm2，增幅為64.5%，雖然增幅未至最大，但其沉陷水域面積的絕對(duì)增量在監(jiān)測(cè)時(shí)段內(nèi)為最大值。

表4 2005-2017年沉陷水域面積及變化值

綜合圖3和表4可知，2005年淮南潘謝礦區(qū)由于已采礦井?dāng)?shù)量較少，還未出現(xiàn)大范圍的采煤沉陷積水，只有潘一礦、潘二礦、潘三礦、顧北礦、張集礦出現(xiàn)小范圍的沉陷積水；2005-2009年，潘一礦、潘二礦、潘三礦、顧北礦、張集礦沉陷水域面積增加，其中謝集礦、張集礦擴(kuò)張最為明顯；2009-2013年期間，除朱集礦外，研究區(qū)內(nèi)各礦沉陷水域面積均有所變化，顧橋礦水域增加明顯，2013年其內(nèi)部產(chǎn)生一大面積沉陷積水區(qū)；2013-2017年，研究區(qū)境內(nèi)所有沉陷水域面積均擴(kuò)有增加，且分布均勻。

圖3 2005-2017年潘謝礦區(qū)沉陷水域分布

圖4 累積原煤產(chǎn)量與沉陷水域面積相關(guān)性分析

淮南潘謝礦區(qū)2005-2017年累積原煤產(chǎn)量與沉陷水域面積關(guān)系如圖4所示，二者相關(guān)系數(shù)r=0.94。隨著累積煤炭產(chǎn)量的增加以及礦區(qū)沉陷水域面積不斷擴(kuò)大，根據(jù)淮南礦業(yè)集團(tuán)生產(chǎn)規(guī)劃，到2020年潘謝礦區(qū)累積原煤產(chǎn)量預(yù)計(jì)達(dá)10.85億t；2020-2030年，原煤產(chǎn)量穩(wěn)定在7000萬(wàn)t/a，累計(jì)原煤產(chǎn)量達(dá)17.85億t。根據(jù)上述擬合結(jié)果預(yù)計(jì)到2020年潘謝礦區(qū)沉陷水域面積總量預(yù)計(jì)增至5968.12 hm2；2021-2030年期間，區(qū)內(nèi)沉陷水域面積年均增量預(yù)計(jì)為350 hm2，到2030年沉陷水域面積將達(dá)到9468.12 hm2。

4 結(jié)論

本研究在對(duì)比分析各種水體信息方法的基礎(chǔ)上，以淮南潘榭礦區(qū)為例，對(duì)研究區(qū)內(nèi)所有水體進(jìn)行了提取，結(jié)合潘謝礦區(qū)2015年采煤沉陷邊界，利用GIS空間疊加，實(shí)現(xiàn)對(duì)研究區(qū)2005-2017年地表沉陷水體的動(dòng)態(tài)監(jiān)測(cè)，分析了區(qū)內(nèi)煤炭開采對(duì)沉陷水域面積變化的影響并做出預(yù)測(cè)，得出以下結(jié)論：

(1)依據(jù)研究所用數(shù)據(jù)特點(diǎn)以及研究區(qū)地理?xiàng)l件，對(duì)比了單波段閾值法、譜間關(guān)系法、歸一化水體指數(shù)法及改進(jìn)歸一化水體指數(shù)法4種水體信息提取方法的優(yōu)劣，選擇提取效果最好、精度最高的改進(jìn)歸一化水體指數(shù)法作為潘謝礦區(qū)水體信息提取方法，其總體精度達(dá)96.03%，Kappa系數(shù)為0.88。

(2)2005-2017年潘謝礦區(qū)沉陷水域的監(jiān)測(cè)結(jié)果顯示，2005-2017年，區(qū)內(nèi)沉陷水域面積增至4948.65 hm2，總增量為4003.52 hm2，年均增長(zhǎng)333.63 hm2，年均增長(zhǎng)率為20.1%。其中，2005-2009年沉陷水域面積增加1459.3 hm2，增幅達(dá)154%，變化幅度相對(duì)最大，2013-2017年沉陷水域面積增量為1940.9 hm2，沉陷水域面積的絕對(duì)增量在監(jiān)測(cè)時(shí)段內(nèi)為最大值。

(3)通過對(duì)監(jiān)測(cè)時(shí)段內(nèi)的累積原煤產(chǎn)量與沉陷水域面積的相關(guān)性定量分析可知，隨著煤炭開采量的累積增加，潘謝礦區(qū)沉陷水域面積隨之不斷上升，二者具有很強(qiáng)的相關(guān)性，相關(guān)系數(shù)為0.94。依據(jù)生產(chǎn)規(guī)劃和擬合結(jié)果，預(yù)計(jì)到2020年，潘謝礦區(qū)沉陷水域面積總量將增至5968.12 hm2，2021-2030年，其沉陷水域面積年均增量預(yù)計(jì)為350 hm2，到2030年沉陷水域面積將達(dá)到9468.12 hm2。