基于時序TM影像的淮南潘謝礦區采煤沉陷水域動態監測

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

(1.中國礦業大學(北京)土地復墾與生態重建研究所，北京市海淀區，100083；2.浙江大學公共管理學院，浙江省杭州市，310058)

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

東部高潛水位礦區煤炭開采引起的地表沉陷積水，嚴重影響了當地的農業發展和居民生活。國內礦區生態環境方面的研究大多集中在生態環境修復、沉降監測以及環境影響評估等方面，對于礦區采煤沉陷積水的研究主要集中在沉陷積水動態監測以及范圍確定等方面?？焖偾揖珳诗@取沉陷水域時空變化信息是高潛水位煤糧復合區進行環境整治與沉陷地復墾工作的前提。與傳統手段相比較，遙感監測方法具有快速、精確、宏觀性強和低費用等特點，在沉陷水域的長時序動態監測上具有無可比擬的優勢。

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

1 研究區及數據源

1.1 研究區概況

潘謝礦區位于安徽省淮南市西北部，地形平坦開闊，地處淮河沖積平原，整體地勢西北高、東南低，地表高程一般為21.5～25.4 m，平均標高在23.1 m左右。所在區域屬季風溫暖帶半濕潤氣候，季節性表現為夏季炎熱，冬季寒冷。6～8月降雨量約占全年的40%，年平均降雨量865.5 mm，年平均蒸發量1610.44 mm，蒸發量大于降雨量，潮濕系數近似0.5。區內溝渠縱橫，水利設施完善，當地潛水位埋深約為1.5 m。因此地下煤炭開采沉陷后，地表極易形成積水，如圖1所示。截至2014年，淮南礦區累計原煤產量約6.7億t，根據生產規劃，到2020年，礦區累計原煤產量約11億t；2021-2030年，原煤產量穩定在7000萬t/a。

圖1 淮南潘榭礦區遙感影像

1.2 數據選擇及預處理

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

表1 研究區遙感數據基本情況

2 研究方法

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

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

目前利用遙感影像進行水體信息提取的方法主要有單波段(近紅外)閾值法、多波段譜間關系法和水體指數法等。單一的水體信息提取方法無法通用，需要根據研究數據情況以及區域地理條件擇優選擇。

2.1.1 單波段閾值法

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

2.1.2 多波段譜間關系法

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

(1)

通過多次反復實驗，選取K1及K2值，使用波段運算工具，得到水體信息提取圖。

2.1.3 歸一化水體指數法

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

(2)

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

2.1.4 改進歸一化水體指數法

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

(3)

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

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

3 結果與分析

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

經過前期預處理，對實驗區OLI遙感影像應用各水體信息自動提取方法，得到提取結果，見表2。

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

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

使用ENVI5.3軟件對實驗區遙感原始影像目視解譯，均勻選擇100組樣本，包含像元5562個。經統計得到水體提取混淆矩陣和Kappa統計，見表3。由表3可知，4種水體提取方法精度均達到90%以上，但MNDWI法提取研究區水體信息的用戶精度、總體精度以及Kappa相關系數為各方法中最高。

表3 不同水體方法提取結果精度評價

3.2 歷年水體信息提取

基于MNDWI指數法，用2005年、2008年、2009年、2013年、2015年和2017年TM/OLI影像潘謝礦區整體水域進行提取，結果如圖2所示。

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

圖2 監測年份潘謝礦區地表水域分布

3.3 采煤沉陷水域監測

由于本文主要側重于分析開采沉陷導致的水體變化，而礦區內仍存在大面積的自然水體，為了區分未受擾動的天然水體與開采沉陷水體，實現對研究區內采煤沉陷水域的提取監測，將淮南潘謝礦區2015年采煤沉陷邊界與水域提取結果結合，利用Arcgis10.3軟件進行空間疊加運算，獲得2005-2017年共6期潘謝礦區沉陷區范圍內的水域結果，如圖3所示，并在SPSS軟件中統計分析潘謝礦區監測時段的沉陷水域面積變化情況，見表4。

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

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

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

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

圖4 累積原煤產量與沉陷水域面積相關性分析

淮南潘謝礦區2005-2017年累積原煤產量與沉陷水域面積關系如圖4所示，二者相關系數r=0.94。隨著累積煤炭產量的增加以及礦區沉陷水域面積不斷擴大，根據淮南礦業集團生產規劃，到2020年潘謝礦區累積原煤產量預計達10.85億t；2020-2030年，原煤產量穩定在7000萬t/a，累計原煤產量達17.85億t。根據上述擬合結果預計到2020年潘謝礦區沉陷水域面積總量預計增至5968.12 hm2；2021-2030年期間，區內沉陷水域面積年均增量預計為350 hm2，到2030年沉陷水域面積將達到9468.12 hm2。

4 結論

本研究在對比分析各種水體信息方法的基礎上，以淮南潘榭礦區為例，對研究區內所有水體進行了提取，結合潘謝礦區2015年采煤沉陷邊界，利用GIS空間疊加，實現對研究區2005-2017年地表沉陷水體的動態監測，分析了區內煤炭開采對沉陷水域面積變化的影響并做出預測，得出以下結論：

(1)依據研究所用數據特點以及研究區地理條件，對比了單波段閾值法、譜間關系法、歸一化水體指數法及改進歸一化水體指數法4種水體信息提取方法的優劣，選擇提取效果最好、精度最高的改進歸一化水體指數法作為潘謝礦區水體信息提取方法，其總體精度達96.03%，Kappa系數為0.88。

(2)2005-2017年潘謝礦區沉陷水域的監測結果顯示，2005-2017年，區內沉陷水域面積增至4948.65 hm2，總增量為4003.52 hm2，年均增長333.63 hm2，年均增長率為20.1%。其中，2005-2009年沉陷水域面積增加1459.3 hm2，增幅達154%，變化幅度相對最大，2013-2017年沉陷水域面積增量為1940.9 hm2，沉陷水域面積的絕對增量在監測時段內為最大值。

(3)通過對監測時段內的累積原煤產量與沉陷水域面積的相關性定量分析可知，隨著煤炭開采量的累積增加，潘謝礦區沉陷水域面積隨之不斷上升，二者具有很強的相關性，相關系數為0.94。依據生產規劃和擬合結果，預計到2020年，潘謝礦區沉陷水域面積總量將增至5968.12 hm2，2021-2030年，其沉陷水域面積年均增量預計為350 hm2，到2030年沉陷水域面積將達到9468.12 hm2。