胡 曉，明安遠(yuǎn)，李士玉，趙傳華，段建國

(1.山東農(nóng)業(yè)大學(xué) 資源與環(huán)境學(xué)院，山東 泰安 271018；2.山東農(nóng)業(yè)大學(xué) 信息科學(xué)與工程學(xué)院，山東 泰安 271018)

兗州地處魯西南平原的東北部，境內(nèi)煤炭資源豐富，為全國八大煤炭生產(chǎn)基地之一，煤炭資源對兗州的經(jīng)濟(jì)發(fā)展起到了重要的作用[1-2]。然而在煤礦生產(chǎn)開采過程中也帶來了許多負(fù)面影響，煤層采空引起地面沉降和塌陷是本地區(qū)最突出的地質(zhì)災(zāi)害[3]。由于地下潛水位高，地面塌陷導(dǎo)致地下水位容易上升到地表標(biāo)高以上，形成大面積的常年或季節(jié)性積水區(qū)，造成了水資源空間分布的紊亂[4-5]。此外兗州礦區(qū)的煤糧復(fù)合面積大，地面塌陷對耕地資源的利用和保護(hù)也產(chǎn)生了較大威脅。因此，客觀分析兗州耕地和水域資源的時(shí)空變化情況，可為我國其它高潛水位煤礦區(qū)城市的耕地和水域資源的利用保護(hù)提供借鑒。

Landsat TM/OLI遙感影像因具有較好的時(shí)序連續(xù)性、數(shù)據(jù)獲取免費(fèi)等特點(diǎn)，被廣泛地應(yīng)用于耕地和水域信息的提取。在耕地提取方面，牛魯燕等[6]以Landsat8 OLI為數(shù)據(jù)源對山東省的耕地進(jìn)行了提取，與詳查數(shù)據(jù)相比，精度達(dá)到了91.8%；滿衛(wèi)東等[7]為揭示1990～2013年東北地區(qū)的耕地變化規(guī)律，利用Landsat TM/ETM+/OLI遙感影像，提取東北地區(qū)1990年、2000年和2013年耕地信息，輔以氣候、地形、社會經(jīng)濟(jì)等數(shù)據(jù)，分析耕地時(shí)空變化特征及其驅(qū)動因素；趙曉麗等[8]基于Landsat TM/ETM+、CBERS1/2和HJ-1等遙感影像，依靠全數(shù)字人機(jī)交互判讀系統(tǒng)對中國1987～2010年的耕地變化進(jìn)行了全面分析，揭示了中國近30年來耕地變化的時(shí)空特征，分析了其變化的主要原因；Zuo等[9]利用多時(shí)相Landsat TM/ETM+遙感數(shù)據(jù)分析了20世紀(jì)80年代末至2010年6個(gè)時(shí)期西北干旱和半干旱區(qū)耕地的時(shí)空演變特征；歐陽玲等[10]利用 2000年、2013年Landsat TM/ETM+/OLI遙感數(shù)據(jù)和2351個(gè)樣點(diǎn)土壤實(shí)測數(shù)據(jù)，對吉林省中西部耕地分布的時(shí)空變化和特征進(jìn)行了分析；Wang等[11]利用1999年和2011年Landsat TM/ETM+影像研究了延川縣耕地轉(zhuǎn)化時(shí)空格局對“退耕還林工程”的響應(yīng)。在水域提取方面，李晶等[12]在比較改進(jìn)歸一化差異水體指數(shù)法、單波段閾值法、譜間關(guān)系法、K-T變換等4種水體提取方法的精度及優(yōu)缺點(diǎn)基礎(chǔ)上，采用基于閾值分割改進(jìn)的歸一化差異水體指數(shù)法，在Landsat TM/OLI影像上提取了兗州煤田1990～2014年的水體信息并分析了其時(shí)空變化特征；Lobo等[13]利用Landsat MSS/TM/OLI遙感影像分析了1973～2013年亞馬遜水域受金礦開采活動影響的時(shí)空變化情況；Jain等[14]對比單波段密度分割法、歸一化差異水體指數(shù)法、譜間關(guān)系法和非監(jiān)督分類法對TM數(shù)據(jù)水體信息提取的精度。

通過前人研究可以發(fā)現(xiàn)，利用Landsat TM/OLI遙感影像開展耕地和水域時(shí)空變化的研究雖然較多，但是鮮有針對礦區(qū)城市耕地和水域時(shí)空變化的研究，多數(shù)研究針對耕地或者水域進(jìn)行單一研究，將二者結(jié)合進(jìn)行共同研究的還較少。此外，礦區(qū)城市耕地和水域的時(shí)空變化特點(diǎn)、二者變化關(guān)系、煤炭開采是否是二者變化的主要驅(qū)動力等問題還有待進(jìn)一步研究。

1 材料與方法

1.1 研究區(qū)概況

兗州隸屬于山東省濟(jì)寧市，地理坐標(biāo)介于東經(jīng)116°35′21″～116°53′36″，北緯35°23′31″～35°43′17″，總面積535 km2。境內(nèi)氣候?qū)倥瘻貣|亞季風(fēng)區(qū)大陸性氣候，四季分明，降水較為充沛，呈春季干燥多風(fēng)、夏季炎熱多雨、秋季溫和涼爽、冬季干冷等氣候特點(diǎn)。土壤以褐土、潮土、砂姜黑土為主，農(nóng)業(yè)發(fā)展歷史悠久，主要的糧食作物包括小麥、玉米、花生、大豆等。境內(nèi)有泗河、洸府河、白馬河、南泉河水系，含一級支流14條，二級支流4條。煤炭是兗州最具優(yōu)勢的礦產(chǎn)資源，面積約占全區(qū)國土面積的37.06%，兗州是濟(jì)寧市也是山東省重要的煤炭生產(chǎn)基地之一。

1.2 數(shù)據(jù)來源及預(yù)處理

影像數(shù)據(jù)：本文所使用的Landsat TM/OLI遙感影像是從美國地質(zhì)調(diào)查局(USGS)下載的經(jīng)LEDAPS(The Landsat Ecosystem Disturbance Adaptive Processing System)程序處理的地表反射率數(shù)據(jù)，影像列號122，行號035，空間分辨率為30 m。數(shù)據(jù)選擇時(shí)盡量選擇降水量、氣溫、月份相近，云量及陰影較少的影像，本研究均選取了5月份的影像,如表1所示。因上述數(shù)據(jù)已經(jīng)過系統(tǒng)的輻射校正、幾何校正等相關(guān)預(yù)處理，本文對遙感影像所做的預(yù)處理工作包括影像的波段合成、影像裁剪等。

表1 遙感影像基本信息

其它數(shù)據(jù)：土地利用現(xiàn)狀圖來源于兗州國土資源局，比例尺為1: 50 000，時(shí)間為2000年、2006年、2013年和2018年；濟(jì)寧市采煤塌陷地現(xiàn)狀圖來源于濟(jì)寧市采煤塌陷地治理中心；降水量數(shù)據(jù)來源于“中國氣象數(shù)據(jù)網(wǎng)”(http://data.cma.cn/)的中國地面降水月值0.5°×0.5°格點(diǎn)數(shù)據(jù)集(V2.0)；煤炭開采量數(shù)據(jù)來源于兗州煤業(yè)股份有限公司(http://www.yanzhoucoal.com.cn)的年度報(bào)告以及部分調(diào)查數(shù)據(jù)。

1.3 耕地提取

在耕地提取方面，最常用的方法主要是監(jiān)督分類[15-16]，本文使用ENVI 5.3軟件對4期影像中的耕地進(jìn)行提取。首先根據(jù)土地利用現(xiàn)狀圖，將兗州的土地類型分成六類：耕地、林地、草地、水域、建筑用地、未利用地[17]；其次定義各土地類型的訓(xùn)練樣本，樣本的可分離性均大于1.9，屬于合格樣本；最后采用最大似然進(jìn)行分類提取。

為評價(jià)耕地提取結(jié)果，在研究區(qū)內(nèi)隨機(jī)選取了200個(gè)驗(yàn)證樣本點(diǎn)，其中，耕地樣本點(diǎn)100個(gè)，非耕地樣本點(diǎn)100個(gè)，利用隨機(jī)樣本對兩種方法分別提取的4期影像耕地精度進(jìn)行評價(jià)。

1.4 水域提取

改進(jìn)的歸一化差異水體指數(shù)(MNDWI)是在歸一化差異水體指數(shù)(NDWI)模型的基礎(chǔ)上，利用短波紅外波段代替NDWI中的近紅外波段創(chuàng)建的新型水體指數(shù)[18]，計(jì)算公式見式(1)和(2)。它增強(qiáng)了水體和建筑物的反差，提高了水域信息的提取精度，不僅適用于城市范圍，也同樣適用于礦區(qū)[19]。

(1)

(2)

1.5 變化檢測

利用ArcGIS 10.6中的統(tǒng)計(jì)分析工具分別計(jì)算各期影像的耕地和水域面積，通過將4期影像的面積做差值運(yùn)算等，得到年際間耕地和水域面積時(shí)間變化及變化速率；將各期提取影像進(jìn)行疊置分析，得到年際間耕地和水域的空間變化，并應(yīng)用SPSS等軟件分析其變化特征與煤炭開采、降水量因素的關(guān)系，以揭示兗州耕地和水域變化的驅(qū)動機(jī)制。

2 結(jié)果與分析

2.1 耕地提取與時(shí)空變化分析

2.1.1 提取結(jié)果

采用監(jiān)督分類提取的兗州耕地信息如圖1所示，提取精度如表2所示。

表2 耕地提取精度

圖1 監(jiān)督分類法提取耕地結(jié)果

2.1.2 耕地時(shí)空變化分析

統(tǒng)計(jì)4期影像提取的耕地面積及變化率如表3所示。

將提取耕地面積做差值運(yùn)算并疊加，得到了兗州年際間耕地空間變化，如圖2所示。時(shí)間變化方面，兗州的耕地面積從2000年的480.19 km2減少至2018年的412.43 km2，減少了67.76 km2，年均減少3.76 km2，年減少率為0.84%。其中，2006～2013年耕地減少較大，年均減少率1.05%；2003～2018年耕地年均減速相對最大為1.15%。

表3 各觀測年耕地面積及變化率

空間變化方面，從2000～2006年，耕地減少的區(qū)域主要集中在東部主城區(qū)周圍、東南部的興隆莊、大統(tǒng)、楊村、鮑店、田莊等礦區(qū)以及西南部的許廠礦區(qū)，兗州北部和西部區(qū)域耕地面積呈現(xiàn)增加趨勢，如圖2(a)所示；而從2006～2013年，耕地減少的范圍進(jìn)一步擴(kuò)大，并延伸到了西部新驛礦區(qū)，增加的耕地主要集中在東南部礦區(qū)，如圖2(b)所示；2013～2018年，伴隨著土地復(fù)墾、塌陷地治理等工作的全面開展，兗州境內(nèi)幾大煤礦區(qū)的耕地面積呈現(xiàn)小幅增加，如圖2(c)所示。

圖2 年際間耕地空間變化

2.2 水域提取與時(shí)空變化分析

2.2.1 水域提取結(jié)果

采用MNDWI水體指數(shù)提取的兗州水域如圖3所示，并利用隨機(jī)建立的10個(gè)水域感興趣區(qū)域(ROI)對提取精度進(jìn)行評定，由于水體在影像上容易識別，而且數(shù)量較少，經(jīng)過統(tǒng)計(jì)總體分類精度都高于90.00%，Kappa系數(shù)均高于0.9，這表明水體的提取較為準(zhǔn)確。

圖3 各觀測年兗州地表水域分布

2.2.2 水域時(shí)空變化分析

統(tǒng)計(jì)4期影像提取的水域面積及變化率如表4所示，并將提取水域面積做差值運(yùn)算并疊加，得到了兗州年際間水域空間變化，如圖4所示。

表4 各觀測年水域面積及變化率

時(shí)間變化方面，兗州的水域面積從2000年的2.77 km2增加至2018年的7.24 km2，增加了6.37 km2，年均增加0.35 km2，年增加率為0.07%。其中，2000～2006年增加了4.47 km2，水域面積增加較大，年均增加率0.75%，年均增速相對最大；2006～2013年增加了0.78 km2，2003～2018年水域面積增加了1.12 km2，年均增加率降低至0.03%。

空間變化方面，從圖4可以看出兗州天然水體較少，2000～2018年間增加的水體主要集中在兗州東南部，為興隆莊、楊村、鮑店、田莊等大中型煤礦區(qū)所在地。2000～2006年兗州東南部水域出現(xiàn)大面積增加，如圖4(a)所示，這說明煤炭長時(shí)間、大面積的開采引發(fā)地面塌陷，致使?jié)撍簧?，形成了積水塌陷地，伴隨著開采活動的劇烈，積水區(qū)的面積和深度進(jìn)一步增大[20-21]。2006～2013年水域面積有所減少，如圖4(b)所示，而水域面積的增加速率從0.17%降低至0.03%。2013～2018年兗州新增積水面積較少，這說明塌陷后主要增加的是原水體的深度。

圖4 年際間水域空間變化

2.3 耕地、水域變化驅(qū)動力分析

本文從煤炭開采和降水量兩個(gè)方面分析其對耕地和水域變化的貢獻(xiàn)度。本文中降水量，指影像獲取時(shí)間相應(yīng)年的年降水量，分別為529.9 mm，766.2 mm，580.3 mm，686.7 mm；煤炭累計(jì)開采量為同一觀測期幾年內(nèi)主要煤礦的煤炭累積開采量，分別為2 746萬t、5 605萬t、7 380萬t、10 590萬t。

兗州年降水量和水域面積之間的相關(guān)系數(shù)為r=0.602，顯著性P>0.4，相關(guān)程度不明顯，這說明降水量不是兗州水域面積增加的主要原因。相比較而言，觀測期兗州煤炭開采量與水域面積的相關(guān)系數(shù)r=0.922，顯著性P<0.07，呈顯著正相關(guān)，如圖5(a)所示，即煤炭開采量的增加，導(dǎo)致塌陷積水區(qū)面積增加，是影響兗州水域面積變化的主要原因。

圖5 累計(jì)開采量與面積的相關(guān)性

觀測期兗州煤炭開采量與耕地面積的相關(guān)系數(shù)為r=-0.967，顯著性P<0.03，呈顯著負(fù)相關(guān)，如圖5(b)所示，即煤炭開采量的增加，是導(dǎo)致兗州耕地減少的主要因素。究其原因，受采動影響的地表從原有標(biāo)高向下沉降，地面坡度發(fā)生變化，改變了地表原有的形態(tài)，引起了高低、坡度及水平位置的變化，逐漸形成一個(gè)比采空區(qū)大的多的塌陷盆地。特別是在高潛水位地區(qū)，地表下沉后抬高潛水位，使地下水埋深下降，在華東地區(qū)地表一般下沉1 m左右，盆地內(nèi)部就可形成季節(jié)性積水或常年積水，被水淹的土壤逐漸喪失了其耕種功能，而且受積水影響，周邊地區(qū)土壤的沼澤化和鹽漬化現(xiàn)象加劇，致使農(nóng)作物枯死和減產(chǎn)、絕產(chǎn)，耕地面積大量減少。

3 結(jié) 論

本文以東部高潛水位煤礦區(qū)城市—兗州為例，利用Landsat TM/OLI遙感影像提取了兗州2000、2006、2013、2018年的耕地和水域信息，分析了2000～2018年間的耕地和水域時(shí)空變化信息。從2000～2018年，兗州耕地總體呈減少趨勢，減少了67.76 km2，主要集中在煤礦區(qū)；而兗州的水域面積則從2000年的2.77 km2增加至2018年的7.24 km2，增加了6.37 km2，主要集中在兗州東南部興隆莊、鮑店、楊村等礦區(qū)。兗州4個(gè)觀測期的水域面積和年降水量之間相關(guān)程度不明顯；而煤炭開采量與耕地面積的相關(guān)系數(shù)為-0.967，與水域面積相關(guān)系數(shù)為0.922，結(jié)果表明煤炭開采是研究區(qū)耕地面積減少、水域面積增加的主要驅(qū)動力。

本研究客觀分析了2000年以來兗州礦區(qū)的耕地和水域變化，并證實(shí)了煤炭開采是二者變化的主要驅(qū)動力。但是由于2012年的遙感影像無法獲取，導(dǎo)致定性分析的年份間隔不一致，加之受限于遙感影像空間分辨率、缺乏持續(xù)觀測數(shù)據(jù)等因素的制約，本研究仍存在一些不足。以后研究應(yīng)盡量使用高分辨率的遙感影像，采用面向?qū)ο蟮姆椒ㄩ_展相關(guān)研究。