高潛水位煤礦區(qū)開采擾動的長時序過程遙感監(jiān)測與影響評價

肖 武，陳文琦，何廳廳,趙艷玲,胡振琪

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

中國東部高潛水位地區(qū)包含十四大煤炭基地中的5個，即兩淮基地、魯西(兗州)基地、河南基地、冀中基地、蒙東基地(東北部分)，可采煤炭儲量 189 億t左右，含煤面積超過6 955 km。該區(qū)域80%以上為高產(chǎn)糧田，采煤沉陷后地表易積水，導(dǎo)致厚沖積層形成的優(yōu)質(zhì)耕地變?yōu)樗颍r(nóng)田生態(tài)系統(tǒng)嚴(yán)重受損，直接威脅國家糧食安全和生態(tài)安全。煤炭開采導(dǎo)致地表沉陷積水，是中國東部高潛水位平原煤礦區(qū)土地生態(tài)變化的主要特征之一，掌握煤炭開采影響下的地表水變化，有助于定量評估煤炭開采對土地、生態(tài)與社會的綜合影響效應(yīng)。

實地調(diào)查是監(jiān)測地面變化與沉陷積水和影響評估最準(zhǔn)確的方法，但是受到沉陷積水空間分布大和過程時間長的限制，無法回溯歷史不同開采階段煤炭開采對土地利用的定量變化與影響。由于遙感影像更新頻率高、影像容易獲取以及大面積監(jiān)控等特征，遙感技術(shù)已廣泛應(yīng)用于森林退化、城市擴展、水體變化、耕地拋荒等大范圍資源變化監(jiān)測。對于礦區(qū)來說，礦區(qū)土地變化監(jiān)測可以分為露天開采和地下開采，露天礦區(qū)主要集中在利用植被指數(shù)監(jiān)測采礦活動對地表的擾動,如李晶等利用時序NDVI實現(xiàn)了對草原礦區(qū)植被擾動的監(jiān)測，賈鐸等與李恒凱等分別基于SSA-Mann Kendall與多源時序NDVI對草原露天礦區(qū)與稀土礦區(qū)的土地損毀與恢復(fù)過程進行了分析。對地下開采的礦區(qū)，多是對不同時期進行人工解譯，然后比較解譯結(jié)果來監(jiān)測采礦活動對地表的擾動與影響，比如郝成元和楊志茹、劉培等分別利用MODIS數(shù)據(jù)與CA_Markov模型對潞安礦區(qū)NPP與徐州礦區(qū)地表熱環(huán)境進行了分析，此外，雷達數(shù)據(jù)也被廣泛的應(yīng)用于地下開采礦區(qū)的地表形變觀測。然而，這種多期影像人工解譯分類的方法工作量大且分類誤差極易累積，選擇特定年份某一時間節(jié)點的數(shù)據(jù)，往往只能代表在這一瞬時狀態(tài)地表的植被與土地利用狀態(tài)，無法表征某一時間階段內(nèi)的特征，這給非線性且高時空異質(zhì)性的事件監(jiān)測帶來困難。在研究過程中，由于受到天氣條件比如云量與衛(wèi)星過境時間等多因素影響，各年之間的數(shù)據(jù)也難以選擇在每一年的同一時間進行，這也給監(jiān)測結(jié)果的可信度與可比較性帶來挑戰(zhàn)。對于高潛水位礦區(qū)地面擾動遙感監(jiān)測來說，如何準(zhǔn)確快速的獲取地表積水的時空規(guī)律是關(guān)鍵，李晶等在比較改進歸一化差異水體指數(shù)法、單波段閾值法、譜間關(guān)系法、K-T變換4種水體提取方法的精度及優(yōu)缺點基礎(chǔ)上，采用基于閾值分割的改進歸一化差異水體指數(shù)法提取了兗州煤田1990—2014年的水體信息并分析了其時空變化特征。但是地面水體變化受到多種因素的影響(人工挖掘魚塘、降水量等)，在缺乏地下采礦信息的情況下，如何根據(jù)單一的遙感影像數(shù)據(jù)，分辨采煤沉陷導(dǎo)致的地面積水與人為活動導(dǎo)致的地面挖損水體，以及消除由于水文年際變化(豐水年與枯水年)導(dǎo)致的噪音，成為遙感提取采煤沉陷水體的重大阻礙與難題。受氣候條件以及采煤引起的地表沉陷程度的影響，沉陷區(qū)域的積水面積在年內(nèi)和年間是變化的。因此，用單景影像代表年內(nèi)的積水面積是不合理的，通過比較多景影像檢測礦區(qū)年間積水面積變化，存在一定的局限性。

隨著云計算技術(shù)與時序遙感變化監(jiān)測算法等相關(guān)技術(shù)發(fā)展，利用可獲取的海量遙感衛(wèi)星過境數(shù)據(jù)，為實現(xiàn)精細化與高時空分辨率的地面變化檢測提供了可能。Google 公司推出的 Google Earth Engine (GEE) 平臺當(dāng)前已收集 MODIS,Landsat,Sentinel 等常用遙感數(shù)據(jù)集，可利用在線或離線的編程方式獲取和處理共享數(shù)據(jù)，并基于強大的谷歌云平臺，利用云計算進行遙感數(shù)據(jù)分析與處理，避免了傳統(tǒng)遙感分析模式帶來的數(shù)據(jù)下載、預(yù)處理等繁瑣過程，GEE平臺在數(shù)據(jù)分析與處理方面體現(xiàn)出強大的能力，逐步應(yīng)用于城鎮(zhèn)擴張、生態(tài)環(huán)境質(zhì)量動態(tài)監(jiān)測、植被覆蓋度動態(tài)監(jiān)測等方面。而長時間序列變化監(jiān)測的方法概括為4種主要類型：光譜變量、影像分類、基于光譜軌跡的分析及數(shù)據(jù)融合的方法。其中，ZHU等提出的基于光譜軌跡的CCDC(Continue Change Detection and Classification)算法能夠通過監(jiān)測時間序列中的斷點信息識別突變事件，最初用在監(jiān)測土地覆被變化，目前在遙感影像時間序列變化檢測領(lǐng)域得到廣泛應(yīng)用，包括城市不透水面變化信息提取、草原景觀動態(tài)變化等，但對于煤礦區(qū)沉陷水體的動態(tài)變化識別仍有待進一步研究。

此外，煤礦開采除了導(dǎo)致沉陷積水的產(chǎn)生，該類地區(qū)面臨的另外一個難題是對周圍社會生態(tài)環(huán)境也帶來顯著影響，尤其對于該區(qū)域內(nèi)優(yōu)質(zhì)耕地的農(nóng)業(yè)生態(tài)系統(tǒng)造成嚴(yán)重的損害。煤礦開采在實現(xiàn)經(jīng)濟效益的同時也要考慮社會-生態(tài)的可持續(xù)發(fā)展問題，開采沉陷造成周邊農(nóng)田作物遭受水淹直接影響了作物產(chǎn)量，導(dǎo)致當(dāng)?shù)剞r(nóng)民的利益受到損害，需要識別沉陷影響區(qū)域來確定對農(nóng)作物賠償和征地復(fù)墾范圍。傳統(tǒng)的方法較多使用實地測量來確定，基于下沉參數(shù)的確定進行沉陷預(yù)計，并用水準(zhǔn)測量或GPS測量進行開采沉陷觀測，從而監(jiān)測煤礦開采沉陷識別影響范圍，但是該辦法的缺點在于效率低，成本高，且缺乏客觀的判別標(biāo)準(zhǔn)。當(dāng)前已有學(xué)者通過測度植被信息、植被碳儲量等方式評估礦區(qū)開采對農(nóng)作物生長的影響。對高潛水位礦區(qū)而言，地形改變與潛水位的相對變化，直接導(dǎo)致開采沉陷區(qū)的表層土壤含水量上升。土壤水分在礦區(qū)高強度開采過程中也會受到極大影響，進而導(dǎo)致生態(tài)環(huán)境惡化。土壤水分作為作物生長的基礎(chǔ)條件，與土壤的其他理化性質(zhì)等共同形成適宜作物生長的土壤環(huán)境。土壤長時間處于水分過飽和狀態(tài)，會導(dǎo)致土壤肥力下降，作物受漬害脅迫，直接影響作物生長甚至危及作物存活。農(nóng)田大量減產(chǎn)甚至絕產(chǎn)，耕作系統(tǒng)遭到破壞，對當(dāng)?shù)鼐用竦纳a(chǎn)生活產(chǎn)生極大影響。由于水分在土壤中的滲透作用存在距離衰減規(guī)律，通過分析開采擾動帶來的外緣土壤水分空間變化量化農(nóng)業(yè)用地受損程度與范圍，能夠有效評估農(nóng)業(yè)生產(chǎn)的影響邊界。

因而，為了監(jiān)測高潛水位煤礦開采的動態(tài)擾動以及量化擾動帶來的農(nóng)業(yè)生產(chǎn)效用影響，筆者嘗試基于GEE云計算平臺與CCDC算法，運用Landsat系列遙感光譜數(shù)據(jù)識別地表水體，監(jiān)測每期影像中水體動態(tài)變化中的突變事件從中提取沉陷水體和回填進行土地復(fù)墾的事件，進而分析突變年份斑塊形態(tài)特征，利用形態(tài)學(xué)方法提取沉陷積水斑塊。在此基礎(chǔ)上，運用土壤濕度檢測指數(shù)(Soil Moisture Monitoring Index,SMMI)、 地表水指數(shù)(Land Surface Water Index,LSWI)、可見光-短波紅外干旱指數(shù)(Visible and Shortwave Infrared Drought Index,VSDI)3個遙感指數(shù)表征區(qū)域內(nèi)土壤水分的空間分布格局，分析沉陷水體周圍土壤水分隨距離的空間變化規(guī)律，識別煤礦開采沉陷對耕地生產(chǎn)影響的邊界，以此量化開采擾動帶來的農(nóng)業(yè)生產(chǎn)效用的影響，從而為高潛水位采煤塌陷區(qū)沉陷積水的時序遙感監(jiān)測與沉陷影響范圍評估提供新思路。

1 研究區(qū)概況

兗州煤田位于山東省濟寧市境內(nèi)，位于兗州、曲阜、鄒城和任城4縣區(qū)交界處，煤田南北跨度26.0 km，東西跨度17.6 km，總面積375.4 km(圖1)。礦區(qū)煤層在300 m以下，已開采煤層厚度大約在1～10 m，大多數(shù)為緩傾斜煤層，適用大型機械開采。礦區(qū)煤炭資源豐富，探明儲量30.14億t，可采儲量17.12億t，自20世紀(jì)90年代開始進入大規(guī)模開發(fā)時期，目前礦區(qū)內(nèi)有北宿、興隆莊、鮑店、楊村、東灘、南屯等大中型礦井。該地區(qū)均采取豎井開拓、地下開采的方式開展作業(yè)，主要使用長壁采煤法、條帶式采煤法、水平分層采煤法等采礦方法。由于地區(qū)的潛水位較高，在長期的煤炭開采作業(yè)下，地面形成大面積沉陷積水區(qū)，礦區(qū)生態(tài)系統(tǒng)受到嚴(yán)重影響，大量優(yōu)質(zhì)耕地遭到破壞，對地區(qū)居民的生產(chǎn)、生活產(chǎn)生了極大影響。

圖1 研究區(qū)概況

研究區(qū)內(nèi)地形平坦，水系遍布，主要流經(jīng)河流有泗河、白馬河、沙河、泥河等。氣候溫和，屬于溫帶大陸性季風(fēng)氣候。四季分明，降雨集中，雨熱同期，平均氣溫為13.6～14.9 ℃，平均降水量在580～747 mm。地區(qū)土壤肥沃，光照充足，耕性良好，主要種植作物有小麥、玉米等，是重要的糧食產(chǎn)區(qū)。

2 數(shù)據(jù)與方法

2.1 研究所采用的數(shù)據(jù)及預(yù)處理

為了識別采煤沉陷與復(fù)墾的動態(tài)過程，需要獲取長時序的遙感影像，因而本研究選取了Google Earth Engine(GEE)平臺上提供的Landsat 5,7,8的大氣地表反射率TOA數(shù)據(jù)，空間分辨率為30 m，如圖2所示，在1986—2017年內(nèi)研究區(qū)共有1 430景可用數(shù)據(jù)。為了能夠更加精細準(zhǔn)確地反映地表水體與土壤水分的空間分布，筆者運用了空間分辨率為10 m的Sentinel-2遙感影像數(shù)據(jù)，并獲取更豐富的光譜信息。所有圖像已經(jīng)過幾何校正處理，利用GEE平臺提供的云掩膜算法移除高于50%的高云量數(shù)據(jù)，將掩膜后所有數(shù)據(jù)作為可用數(shù)據(jù)，再進行影像裁剪等預(yù)處理工作。

圖2 研究所采用Landsat影像數(shù)據(jù)概況

2.2 研究方法

筆者以識別開采擾動以及影響范圍評估為目的，首先利用Landsat數(shù)據(jù)生成每期影像水體指數(shù)軌跡數(shù)據(jù)，運用CCDC算法動態(tài)擬合每一像素的軌跡曲線，識別具有突變特征的像素。在此基礎(chǔ)上，監(jiān)測沉陷積水和回填后土地復(fù)墾的區(qū)域并識別年份，運用形態(tài)學(xué)方法識別煤礦開采導(dǎo)致地表沉陷積水范圍，并進行精度驗證。最后對于提取的沉陷積水斑塊進一步構(gòu)建緩沖區(qū)，分析緩沖區(qū)土壤水分空間分布的距離規(guī)律，從而識別開采擾動的影響邊界，研究具體技術(shù)路線如圖3所示。

圖3 研究技術(shù)流程

..基于CCDC算法的沉陷水體動態(tài)監(jiān)測

煤礦開采會造成地表覆被的顯著變化，而高潛水位煤礦區(qū)的特征在于開采所導(dǎo)致的地表沉陷可能會形成大量的沉陷水體。眾多學(xué)者研究表明，運用遙感光譜信息能夠有效提取地表水體，其中調(diào)整歸一化水體指數(shù)(MNDWI)，改進了NDWI提取水體的能力，對水體識別更為敏銳，能夠?qū)⑺w與植被、土壤和建筑物等較好區(qū)分，被廣泛運用于水體的識別提取研究中，因而筆者運用MNDWI指數(shù)來表征區(qū)域內(nèi)水體的變化，有

(1)

其中，為影像中的綠光波段；為影像中的中紅外波段。MNDWI>0時，該像元為水體，MNDWI<0時為非水體。

而沉陷水體并非在某個時間點突然出現(xiàn)，通過時序遙感數(shù)據(jù)能夠觀察MNDWI指數(shù)的變化，從而監(jiān)測區(qū)域內(nèi)開采擾動帶來的沉陷水體變化。高潛水位礦區(qū)地表沉陷引起水體變化包括3個階段。第1階段，開采前MNDWI指數(shù)較低且穩(wěn)定(說明：本研究排除煤炭開采區(qū)域地表屬于水域的情況，這種情況地表沉陷會增加水域深度，但對環(huán)境的影響有限，故研究仍有意義)。第2階段，開采引起地表下沉，塌陷程度不足不會出現(xiàn)積水。隨著塌陷程度的加重，將逐步出現(xiàn)積水，MNDWI指數(shù)隨之上升。第3階段，對于出現(xiàn)的沉陷積水，如果沒有展開復(fù)墾修復(fù)工程，MNDWI指數(shù)將繼續(xù)維持高位，若通過修復(fù)工程對水體進行填平修復(fù)，MNDWI指數(shù)將降到低值，但是回填后土地復(fù)墾區(qū)域再次塌陷可能導(dǎo)致指數(shù)回升。高潛水礦區(qū)開采沉陷積水軌跡也會受降水、遙感影像質(zhì)量等影響，但是煤炭開采沉陷和土地復(fù)墾工作是礦區(qū)地表水體變化軌跡的主導(dǎo)因素。煤礦擾動積水變化所產(chǎn)生的MNDWI指數(shù)的變化過程如圖4所示。

圖4 2種擾動情境下的MNDWI指數(shù)變化示意

為了監(jiān)測開采擾動的動態(tài)過程，筆者運用CCDC算法監(jiān)測地表水體的時序變化。CCDC算法是一種基于時間序列的動態(tài)變化監(jiān)測和分類方法，該算法將每個像素所有可用的Landsat觀測數(shù)據(jù)均運用于估計時間序列模型，并使用這些模型來預(yù)測未來的觀測數(shù)據(jù)。如果新的連續(xù)觀測值超出預(yù)期范圍，則標(biāo)記一個斷點并將序列進行分割，然后從斷點開始重新模擬一個新的時間序列模型，因而變化前后會生成2個時間序列。將這一過程重復(fù)直至擬合所有觀測值，最終識別研究時段內(nèi)所有水體的突變事件。CCDC算法可以統(tǒng)計在時間序列中發(fā)現(xiàn)的斷點信息，包括突變事件出現(xiàn)總量、突變事件光譜指數(shù)的變化幅度與斷點出現(xiàn)的時間，通過檢測長時序的開采擾動事件為沉陷水體識別提供支撐。

..沉陷水體識別與制圖

根據(jù)CCDC算法對MNDWI軌跡的分割結(jié)果，可以獲得一系列的MNDWI片段。為了準(zhǔn)確地識別沉陷水體，選擇需要確定變化幅度閾值來準(zhǔn)確識別水體。通過實驗的方法，筆者選擇100個沉陷水體作為樣本點，設(shè)置參數(shù)在[0.1,0.5]區(qū)間內(nèi)，并以0.05為間隔，測算監(jiān)測沉陷水體準(zhǔn)確性最佳的閾值?；趯嶒灲Y(jié)果判斷：MNDWI指數(shù)增加幅度大于0.2，并且突變后的指數(shù)大于0，則被認(rèn)定為出現(xiàn)沉陷水體。與此同時，筆者用同樣的方法判定水體修復(fù)的突變情況。即認(rèn)為出現(xiàn)恢復(fù)水體的條件為：MNDWI指數(shù)減少幅度大于0.2，并且突變后的MNDWI指數(shù)小于0。在識別沉陷水體與恢復(fù)水體的基礎(chǔ)上，進一步判斷擾動與修復(fù)時間。如果所識別的相同變化方向的斷點僅出現(xiàn)一次，那么將這一突變時間點作為沉陷水體的出現(xiàn)年份，對于恢復(fù)水體也是如此。當(dāng)出現(xiàn)多個連續(xù)同向突變的片段，將擾動的最短時間點作為沉陷水體出現(xiàn)時間，將最長的片段作為恢復(fù)水體的時間。

對于沉陷積水和人工池塘的區(qū)分方式，考慮到地下煤炭開采空間形式分為縱向多煤層開采和橫向單煤層多工作面開采，對地表的擾動空間形態(tài)表現(xiàn)為“單中心”向外擴張和沿開采走向“帶狀”蔓延，“帶狀”形態(tài)強于“單中心”，擾動時間表現(xiàn)為多周期連續(xù)發(fā)生。沉陷積水年份圖斑呈現(xiàn)“時間多樣、空間上大面積連續(xù)出現(xiàn)且輪廓不規(guī)整等特征”。工程水體人工挖掘周期短，面積相對較小且不隨時間擴充，多以“孤島”形式存在。筆者采取形態(tài)學(xué)識別方式，利用圖斑的“時間一致性”和“分維度指數(shù)”剔除面積<4個單位像元的斑塊并提取年份<2的水體斑塊，計算提取斑塊的分維度，分維度≥0.5作為工程水體，從而獲取范圍內(nèi)沉陷積水區(qū)斑。

..數(shù)據(jù)驗證與分析

筆者利用集成概率積分法的沉陷預(yù)計軟件(Mining Subsidence Prediction System，MSPS)預(yù)測礦區(qū)的地表沉陷區(qū)域，計算結(jié)果與積水年份斑塊疊加，驗證本文所用方法在識別采煤擾動區(qū)域的位置精度。為了驗證沉陷積水年份和修復(fù)年份圖斑的精度，考慮到難以獲取高空間分辨率的公開遙感數(shù)據(jù)，筆者以年為單位進行突變時間的精度驗證。筆者在1990—2017年每年選擇40個樣點，經(jīng)歷沉陷積水和回填恢復(fù)的樣點比例為5∶3。利用Google Earth上時序影像數(shù)據(jù)交互目視標(biāo)定樣點的沉陷積水年份和恢復(fù)年份，其中700個樣點經(jīng)歷沉陷積水，375個樣點經(jīng)歷沉陷積水恢復(fù)(1993年之前未發(fā)現(xiàn)回填恢復(fù)區(qū)域，選取樣點的時間段定為1993—2017年)，共1 075個樣點。通過將樣本標(biāo)簽與算法識別結(jié)果進行比較，計算采礦擾動沉陷積水和復(fù)墾恢復(fù)檢測的用戶精度(User’s Accuracy，UA)、生產(chǎn)者精度(Producer’s Accuracy，PA)、總體精度(Overall Accuracy，OA)和Kappa系數(shù)，驗證變化檢測及精度。

..土壤水分的遙感提取

傳統(tǒng)的土壤水分獲取方法主要是針對局部區(qū)域的實驗測量，但是遙感技術(shù)的發(fā)展提供了大區(qū)域尺度監(jiān)測土壤水分的新方式。基于遙感數(shù)據(jù)的土壤水分反演方式已有較為成熟的研究，眾多學(xué)者通過高光譜的豐富信息構(gòu)建估算模型模擬土壤水分的空間分布。研究表明遙感影像的波段信息包括可見光、近紅外波段以及微波波段等，與地表土壤水分具有較強的相關(guān)性。土壤的反射率特征受到土壤中含水量的顯著影響，光學(xué)遙感影像由于較高的空間分辨率，能夠較為精細反映區(qū)域內(nèi)土壤水分的相對情況，因而筆者選取SMMI,LSWI,VSDI三個遙感指數(shù)反演土壤水分。劉英等基于SWIR和NIR波段構(gòu)造了SMMI指數(shù)，并在實驗中表明在礦區(qū)地區(qū)土壤水分變化測度也有著可靠表現(xiàn)。同時已有研究表明XIAO等提出的LSWI指數(shù)是通過選取水分的敏感光譜構(gòu)建植被指數(shù)反饋植物葉片的含水量進而表征土壤水分情況。VSDI由ZHANG等提出，并通過驗證表明與實測土壤水分高度相關(guān)。LSWI指數(shù)與VSDI指數(shù)對土壤水分呈現(xiàn)正相關(guān)的關(guān)系，也就是隨著土壤水分的增加，指數(shù)越高，而SMMI指數(shù)與土壤水分呈負相關(guān)關(guān)系，3個指數(shù)的具體計算式為

(2)

(3)

VSDI=1-(-)+(-)

(4)

式中，為影像中的近紅外波段；為影像中的近紅外波段；為影像中的紅光波段；為影像中的藍光波段。

..開采沉陷對農(nóng)作物影響邊界分析

高潛水位煤礦經(jīng)歷開采擾動后，地表沉陷以及生成的眾多積水區(qū)成為了對農(nóng)作物生長造成持續(xù)損害的影響源，對周圍的土壤水分產(chǎn)生極大影響。但是由于水分在土壤中的滲透作用有限，通過分析土壤水分的變化可以識別影響邊界。

由于小面積的沉陷水體易受季節(jié)性影響，影響評估結(jié)果的穩(wěn)定性，因而筆者根據(jù)2017年的開采沉陷水體空間分布，以10 ha為閾值，選取單個積水區(qū)域面積大于10 ha的28個沉陷水體為研究對象，分析伴隨沉陷水體向外距離的延伸，周圍土壤水分的變化趨勢。為了避免季節(jié)性影響，基于2017年Sentinel-2的合成影像，選取的礦區(qū)內(nèi)沉陷水體區(qū)域為開采擾動中心向外做射線，以30 m為間隔距離對每個水體做緩沖區(qū)分析。考慮到區(qū)域的總體范圍以及沉陷水體之間的間距較小，為了避免水體之間相互影響，共構(gòu)建了750 m范圍的25個緩沖區(qū)。根據(jù)3種指數(shù)反演土壤水分的結(jié)果，計算不同距離緩沖區(qū)內(nèi)土壤水分的均值，隨著距積水區(qū)域的距離增加，基于距離衰減規(guī)律分析開采沉陷對農(nóng)作物生長的影響程度變化，可以確定煤炭開采對周邊農(nóng)業(yè)生態(tài)系統(tǒng)的影響邊界。

然而，緩沖區(qū)內(nèi)土壤水分仍受到人類活動以及復(fù)雜因素影響，為了排除其他影響因素的干擾，本研究將緩沖區(qū)進一步做如下處理：① 由于建設(shè)用地斑塊具有不透水的特性，緩沖區(qū)內(nèi)建設(shè)用地圖斑的土壤水分值不具備研究意義，因而剔除緩沖區(qū)內(nèi)建設(shè)用地斑塊；② 通過遙感影像目視觀察，鑒于研究區(qū)西面有河流經(jīng)過并在河岸建設(shè)堤壩，一定程度組織了沉陷水體對土壤水分的滲透，因而將西面沿河一側(cè)沉陷水體緩沖區(qū)不納入研究范圍。③ 對于研究區(qū)其他較小面積的水體，由于相鄰的積水空間內(nèi)也存在一定的交互作用，為了排除無關(guān)變量的干擾，去除研究區(qū)內(nèi)的其余水體數(shù)據(jù)。

3 研究結(jié)果

3.1 結(jié)果精度驗證

積水年份圖斑與采煤沉陷預(yù)計軟件計算的沉陷區(qū)域疊加結(jié)果顯示(圖5)，積水年份圖斑與沉陷區(qū)域重疊度達到98.12%。積水年份和修復(fù)年份識別的總體精度分別為85%,77%，Kappa系數(shù)分別為84%,76%。對生產(chǎn)者精度與用戶精度而言，積水年份的總體精度均高于修復(fù)年份。2者均處于一定程度的波動，生產(chǎn)者精度主要在73%～88%，用戶精度主要在74%～92%。年份被識別錯誤的情況只出現(xiàn)在個別年份，如1997年和2006年，這可能是由于遙感數(shù)據(jù)的質(zhì)量不高?？傮w而言，積水年份驗證精度高于修復(fù)年份，是因為地表沉陷積水后，回填后土地復(fù)墾成本太高，導(dǎo)致出現(xiàn)沉陷積水區(qū)域未進行回填后土地復(fù)墾，或者只進行小范圍回填后土地復(fù)墾的情況，筆者運用Landsat數(shù)據(jù)集的空間分辨率為30 m，可能忽視了小范圍的土地復(fù)墾現(xiàn)象，增加了修復(fù)年份識別的難度。

圖5 精度驗證分析

3.2 沉陷積水與水體恢復(fù)的時序特征

自煤礦開采活動開展以來，研究區(qū)域逐漸受到開采帶來的擾動，從結(jié)果來看，研究區(qū)自1990年起開始出現(xiàn)沉陷水體(圖6)?？傮w而言，1990—2017年出現(xiàn)的開采沉陷積水總面積共計3 021.08 ha。從空間上看，沉陷積水主要出現(xiàn)在興隆莊礦區(qū)，沉陷水體面積達到936.14 ha。鮑店礦區(qū)與東灘礦區(qū)也出現(xiàn)了大量的沉陷積水現(xiàn)象，其余礦區(qū)沉陷水體現(xiàn)象較少。積水圖斑呈現(xiàn) “倒錐狀”、“橢圓形盆狀”等形狀，積水年份屬性呈現(xiàn)“時間連續(xù)、空間片狀”等特征。從時間上看，年度積水面積呈現(xiàn)梯度變化趨勢，且呈現(xiàn)先增長后減少的趨勢，可以大致分為1990—2000,2001—2011,2011—2017三個階段。其中，在1990—2000年，沉陷積水的面積變化較為平穩(wěn)，且沉陷水體出現(xiàn)面積較小，該階段總計出現(xiàn)沉陷水體755.87 ha，占比25.02%。而2001—2011年，沉陷水體面積沉陷顯著增長趨勢，該階段沉陷水體共計1 780.09 ha，占比58.92%。其中2006年沉陷積水最為嚴(yán)重，面積最大值達到417.64 ha。而2011年之后，沉陷水體面積迅速縮減并且趨于平緩，年平均積水面積為81.17 ha。

圖6 1990—2017年兗州礦區(qū)沉陷積水年份分布

沉陷水體面積的變化趨勢一定程度上可以反映中國煤炭開采周期的時序特征，即伴隨中國經(jīng)濟的發(fā)展，煤炭行業(yè)在20世紀(jì)初年呈現(xiàn)出迅猛發(fā)展的勢頭，但與此同時，開采帶來的不可逆的生態(tài)問題也讓人們意識到生態(tài)保護的重要性，2010年以后隨著逐步進入深部開采以及重視生態(tài)修復(fù)，開采擾動的范圍相對下降。

礦區(qū)沉陷水體的回填復(fù)墾活動從1993年開始出現(xiàn)，累計恢復(fù)面積為888.37 ha(圖7)。沉陷水體土地復(fù)墾工作在空間上主要呈現(xiàn)出2種特征：其一是由點到面，從小范圍的水體修復(fù)出發(fā)，將相鄰小面積的修復(fù)區(qū)域拼接從而實現(xiàn)區(qū)域的生態(tài)修復(fù)；其二呈現(xiàn)由外圍向內(nèi)修復(fù)的特征，出現(xiàn)圓環(huán)狀。早期沉陷水體的回填復(fù)墾力度較為薄弱，恢復(fù)面積較小，呈現(xiàn)零星的空間分布。1993—2005年共回填復(fù)墾125.49 ha，占比14.13%。從2006年起，沉陷水體的回填復(fù)墾面積呈現(xiàn)顯著的增加趨勢，可以看到在2013—2015持續(xù)地大范圍進行生態(tài)修復(fù)工作，3 a年間充填水體后復(fù)墾的面積達到318.88 ha，占比35.89%。

圖7 1993—2017年兗州礦區(qū)沉陷積水修復(fù)年份分布

3.3 開采沉陷造成的影響范圍確定

通過3種土壤水分的反演，并構(gòu)建緩沖區(qū)進行基于距離的土壤水分統(tǒng)計分析。由圖8可以看出，3種曲線均出現(xiàn)相似的變化趨勢，一定程度上能夠反映土壤水分的空間分異。也就是從開采沉陷中心往外，隨著距離的增加，土壤水分均呈現(xiàn)遞減的趨勢。其中，可以看出，從沉陷積水區(qū)向外120 m范圍內(nèi)，3個指數(shù)均呈現(xiàn)強烈的變化，土壤水分驟減，可以看出，在這一距離內(nèi)，開采擾動對土壤水分的造成了嚴(yán)重影響，沉陷積水隨著土壤滲透作用，顯著提高了土壤水分含量，可能導(dǎo)致植物吸水過多而死亡，進而可能對該區(qū)域內(nèi)農(nóng)業(yè)生產(chǎn)和生態(tài)環(huán)境也帶來較大損害，是開采擾動的嚴(yán)重影響區(qū)。在120～300 m，土壤水分仍呈現(xiàn)出減少趨勢，但變化幅度較之前平緩很多，說明仍處于開采擾動的影響范圍中，只是程度較輕。因而，該區(qū)域的農(nóng)業(yè)生產(chǎn)一定程度上也難以維持正常水平，未來仍要重視這一區(qū)域的生態(tài)修復(fù)工作。而在沉陷水體的300 m范圍之外，3種指數(shù)的均處于較為平緩的狀態(tài)，土壤水分的變化幅度較小，開采沉陷造成的影響接近于無。但是由于在與開采沉陷源距離較遠，可能受到外界因素的干擾，例如種植作物的不同、與河流的距離縮短，可能呈現(xiàn)出略有波動的情況。

圖8 與沉陷水體不同距離的土壤水分變化

4 結(jié)論與展望

(1)兗州礦區(qū)開采30 a來共出現(xiàn)沉陷水體3 021.08 ha，呈現(xiàn)先增加后減少的趨勢。1990年起礦區(qū)開始出現(xiàn)沉陷水體，但早期開采強度不高，沉陷水體面積較少。沉陷水體大量出現(xiàn)在2001—2011年，并在2006年達到峰值，在此階段出現(xiàn)了高強度的煤礦開采工作，而在這一階段之后又恢復(fù)較低的開采強度。除了控制開采強度，礦區(qū)也同步對沉陷水體進行恢復(fù)回填，1993—2017年已完成恢復(fù)回填水體888.37 ha，其中自2007年起，沉陷水體恢復(fù)工程開始大范圍開展。目前國家已逐漸意識到自然資源的稀缺性與生態(tài)修復(fù)的重要性，對礦區(qū)開采和復(fù)墾都提出了嚴(yán)格的要求，為資源的可持續(xù)利用和地區(qū)的可持續(xù)發(fā)展尋找新路徑。

(2)運用LSWI,SMMI,VSDI三種遙感指數(shù)表征土壤水分，以開采產(chǎn)生的沉陷水體為中心沿射線向外均呈現(xiàn)隨距離的增加，擾動效用遞減的變化趨勢。其中，距離沉陷水體0～120 m為嚴(yán)重影響區(qū)，120～300 m為中度影響區(qū)，300 m以外可視作幾乎無影響區(qū)。從沉陷水體擾動角度分析高潛水位礦區(qū)開采的影響范圍，可以進一步為礦區(qū)損毀復(fù)墾以及農(nóng)業(yè)生態(tài)修復(fù)提供支撐。但是本研究僅基于遙感指數(shù)反演分析土壤水分的變化趨勢，并沒有準(zhǔn)確量化出真正的土壤水分情況，未來可以依據(jù)土壤水分實測數(shù)據(jù)進行更準(zhǔn)確的分析。

(3)在GEE平臺上，識別的沉陷積水和回填后土地復(fù)墾水體的總體準(zhǔn)確度>77%，具有較高的魯棒性。本研究驗證了長時序Landsat數(shù)據(jù)在監(jiān)測礦區(qū)地表沉陷積水的實用性，繪制的地表沉陷積水變化年份圖可為土地修復(fù)工作提供科學(xué)依據(jù)。但是本研究仍存在一定局限性：① 由于本研究采用了30 m分辨率的遙感數(shù)據(jù)，可能對水體識別的精確程度產(chǎn)生一定影響，難以識別細碎的開采沉陷積水及土地復(fù)墾現(xiàn)象；② 研究中運用了1986—2017年可獲得的所有Landsat遙感影像，部分年份由于氣候原因(例如云、陰影等)影像質(zhì)量不佳，獲取不到足夠的影像信息，可能導(dǎo)致沉陷積水及土地復(fù)墾水體識別的準(zhǔn)確性降低；③ 對于長時序遙感影像的變化監(jiān)測研究，驗證精度存在年份間的波動情況，這是由于基于像素的年度影像分析可能存在識別的地表水體變化年份與相鄰年份混淆的情況，同時 “椒鹽”現(xiàn)象也造成了一定的噪聲干擾。

下一步研究將考慮利用多源遙感影像數(shù)據(jù)分析礦區(qū)開采沉陷水體與復(fù)墾水體，從而更加精細地識別沉陷水體的動態(tài)變化，以便將該方法應(yīng)用于其他高潛水位礦區(qū)的地表沉陷積水監(jiān)測，為大范圍地表沉陷積水制圖提供有力支撐。