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

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

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

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

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

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

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

1 研究區(qū)概況

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

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

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

2.2 研究方法

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

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

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

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

(1)

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

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

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

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

..沉陷水體識(shí)別與制圖

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

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

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

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

..土壤水分的遙感提取

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

(2)

(3)

VSDI=1-(-)+(-)

(4)

式中，為影像中的近紅外波段；為影像中的近紅外波段；為影像中的紅光波段；為影像中的藍(lán)光波段。

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

高潛水位煤礦經(jīng)歷開(kāi)采擾動(dòng)后，地表沉陷以及生成的眾多積水區(qū)成為了對(duì)農(nóng)作物生長(zhǎng)造成持續(xù)損害的影響源，對(duì)周?chē)耐寥浪之a(chǎn)生極大影響。但是由于水分在土壤中的滲透作用有限，通過(guò)分析土壤水分的變化可以識(shí)別影響邊界。

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

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

3 研究結(jié)果

3.1 結(jié)果精度驗(yàn)證

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

圖5 精度驗(yàn)證分析

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

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

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

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

礦區(qū)沉陷水體的回填復(fù)墾活動(dòng)從1993年開(kāi)始出現(xiàn)，累計(jì)恢復(fù)面積為888.37 ha(圖7)。沉陷水體土地復(fù)墾工作在空間上主要呈現(xiàn)出2種特征：其一是由點(diǎn)到面，從小范圍的水體修復(fù)出發(fā)，將相鄰小面積的修復(fù)區(qū)域拼接從而實(shí)現(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)顯著的增加趨勢(shì)，可以看到在2013—2015持續(xù)地大范圍進(jìn)行生態(tài)修復(fù)工作，3 a年間充填水體后復(fù)墾的面積達(dá)到318.88 ha，占比35.89%。

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

3.3 開(kāi)采沉陷造成的影響范圍確定

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

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

4 結(jié)論與展望

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

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

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

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