郭千慧子，湯伏全，2，馬 婷，柴成富

（1.西安科技大學 測繪科學與技術學院，陜西 西安 710054；2.國土資源部煤炭資源勘查與綜合利用重點實驗室，陜西 西安 710021）

0 引 言

西部黃土高原地區(qū)煤炭儲量豐富但生態(tài)環(huán)境十分脆弱。 受氣候環(huán)境因素和地形地質條件影響，區(qū)域地表植被覆蓋率低、水土流失現(xiàn)象嚴重［1］。 黃土礦區(qū)平地少坡地多，其特殊的地貌條件對農業(yè)種植活動具有明顯的限制作用，地下井工采煤方式造成的土地破壞，導致黃土糧礦復合區(qū)內優(yōu)質土地數(shù)量大幅減少，人－地矛盾日益突出［2－3］。 由采煤活動引起的地表土地沉陷、裂縫等損害嚴重影響著當?shù)剞r業(yè)生產活動［4－5］。 明確黃土礦區(qū)耕地面積及空間分布對區(qū)域農業(yè)資源管理、保護及可持續(xù)發(fā)展具有重要意義。

長期以來針對黃土礦區(qū)開采影響下農耕地的變化缺乏監(jiān)測數(shù)據支撐，礦區(qū)土地利用研究難以深入。遙感技術憑借快速、客觀、不受地面條件限制等優(yōu)點已成為監(jiān)測大范圍耕地種植信息的有效手段。 目前，基于單一時相遙感數(shù)據信息提取在地形破碎地區(qū)具有較大不確定性［6］，農作物生育期具有相對穩(wěn)定的時間變化規(guī)律，其特有的反射光譜特征和季節(jié)時相特征能夠通過時序遙感數(shù)據得以反映，相較于單時相數(shù)據，時序遙感數(shù)據在進行耕地提取時具有明顯的優(yōu)勢。 但由于天氣條件及傳感器重訪周期過長等因素，覆蓋作物完整生長發(fā)育期的時序遙感數(shù)據集構建困難，為獲取同時具備高時空分辨率的影像國內外研究學者展開了大量研究，其中增強型自適應反射率融合模型（ Enhanced Spatial and Temporal Adaptive Reflectance Fusion Model spatiotemporal fusion model，ESTARFM）因其融合精度高、細節(jié)效果好而被廣泛使用［7］。

為此，以彬長大佛寺煤礦為研究區(qū)，利用ES?TARFM 構建高時空分辨率NDVI 數(shù)據集，以該礦區(qū)耕地主要種植作物（冬小麥及夏玉米）為研究對象，結合相關物候信息，采用面向對象決策樹分類法實現(xiàn)礦區(qū)耕地種植面積進行的提取，并通過分析采煤沉陷與耕地利用演化情況的時空相關性，揭示采煤活動對礦區(qū)耕地資源利用的量化影響。

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

1.1 研究區(qū)概況

大佛寺煤礦地處隴東黃土高原，地貌類型以黃土梁塬、溝谷為主，屬典型半干旱黃土溝壑丘陵區(qū)，土地利用類型以林草地、耕地為主，其中耕地約占地全部土地類型的1／3，且主要分布于臺塬和溝谷平坦處，受地形因素影響明顯，綜合考慮Landsat 遙感影像分辨率以及研究區(qū)作物種植情況，所定義耕地僅包括每年耕種、收獲作物種植范圍。

1.2 地下采空區(qū)及地表沉陷影響區(qū)

該煤礦自2006 年投產，截止2019 年底，地下采煤工作面的總面積約12 km2，將該煤炭開采區(qū)域定義為采空區(qū)，根據開采沉陷理論，地下開采引起的地表沉陷區(qū)遠大于采空區(qū)面積，根據概率積分法預計模型及礦區(qū)地表移動參數(shù)估算，將采空區(qū)邊界向外500 m 確定為開采沉陷區(qū)，而開采沉陷區(qū)周圍的地下水位、大氣質量、地表徑流、土壤質量等，都會受到地下開采的間接影響［8－10］，故將沉陷區(qū)邊界以外500 m 范圍劃為間接影響區(qū)。 為了進行對比分析，將間接影響區(qū)以外500 m 范圍作為非開采影響區(qū)，如圖1 所示。

圖1 研究區(qū)地理位置及開采影響區(qū)劃分示意Fig.1 Schematic diagram of geographical location of study area and zoning affected by coal mining

1.3 遙感數(shù)據及預處理

選取作物生育期內連續(xù)的Landsat 影像，對缺失時相影像采用時空融合模型進行預測，參與模型的遙感數(shù)據選用NASA 提供的MOD13Q1 16 d 合成陸地專題數(shù)據集，獲取時間應與缺失的Landsat 數(shù)據相同或鄰近。 預處理首先在IDL 中實現(xiàn)Landsat 數(shù)據的批量輻射定標與大氣校正，然后利用MRT（MODIS Reproject Tools）對MODIS13Q1 數(shù)據集進行投影轉換和鑲嵌，并保證經重投影和重采樣之后的MODIS 數(shù)據應與對應時間的Landsat 影像在像元大小、數(shù)據分辨率、投影方式保持一致。

2 研究方法

2.1 時空數(shù)據融合模型

ESTARFM 模型首先假設高－低分辨率遙感影像傳感器間像元反射率的時間變化趨勢是線性的，通過擬合回歸移動窗口中與目標像素相似的所有像素，進而確定目標像素的大?。?1］，從而獲取具有高時空分辨率的結果影像。 圖2 為隨機選取的時空數(shù)據融合結果影像與真實影像對應樣本像素點的NDVI 值散點圖，從離散點的分布情況可知，融合影像與真實影像數(shù)據間呈線性相關，R2＝0.893 5，通過波段計算得出兩幅影像之間的方差低于0.01，相關系數(shù)達0.946 1，表明兩組數(shù)據具有較強的相關性。即在landsat 數(shù)據缺失的條件下，使用ESTARFM 模型獲取的融合影像能夠代替參與后續(xù)NDVI 時序數(shù)據集的構建。

圖2 真實影像與融合影像NDVI 散點圖Fig.2 Scatterplots of synthetic Landsat images with the original landsat images

2.2 面向對象多尺度分割

黃土高原特殊的地貌條件造成了礦區(qū)耕地分布在空間上的分異性。 考慮地物紋理及空間分布特征對提高分類精度具有重要意義，而確定分割對象大小是保證分類精度的前提［12］。 結合研究區(qū)的實際情況對分割效果圖進行目視判讀，最終確定3 個尺度下的分割結果。 首先在尺度85、形狀因子0.2、緊密度因子0.5 的參數(shù)條件下，區(qū)分水體及非水體；然后在尺度45、形狀因子0.3、緊密度因子0.2 的參數(shù)條件下，區(qū)分非耕地信息；最后在尺度30 下對分割結果進行掩膜獲取耕地邊界。

2.3 時間序列NDVI 數(shù)據集構建

基于相關資料和實地調查，確認研究區(qū)耕地主要種植冬小麥和夏玉米，2 種作物采取輪作模式耕種，在同一時間段內不存在交集。 通過對比農作物在不同生長發(fā)育期內的植被指數(shù)變化特性實現(xiàn)對具備良好耕作條件的在耕耕地的提取。 研究區(qū)主要作 物物候期見表1。

表1 研究區(qū)主要作物物候期Table 1 Phenological period of major crops in study area

為避免云雨、傳感器與地物間幾何特征等因素影響，本文采用TIMESAT3.2 程序包中的S－G 濾波方法對時間序列NDVI 曲線數(shù)據集進行平滑重構處理，以降低噪聲影響。 經預處理后的NDVI 時間序列數(shù)據能夠基本反映地類間的差異關系［14－15］。 圖3為滿足最優(yōu)特征的不同地類樣本數(shù)據集所提取的平均時間序列NDVI 變化曲線。 通過設置最優(yōu)閾值并建立決策樹規(guī)則模型來判別地物類型。

圖3 典型地物時序NDVI 曲線Fig.3 NDVI time series curve of typical features

分類結果精度評價采用目視解譯結合現(xiàn)場采樣，選取符合物候特征的耕地樣本點，其中40%用于確定NDVI 閾值，60%用于檢驗閾值精度，并將分類結果與監(jiān)督分類、決策樹分類進行對比。

3 結果與分析

3.1 分類結果對比分析

將面向對象時序決策樹分類結果與監(jiān)督分類、決策樹分類結果進行對比，如圖4 所示。 監(jiān)督分類法與決策樹分類結果中存在大量的“椒鹽現(xiàn)象”。 由于農作物長勢較好時，林草地也處于繁密階段，因此導致大量耕地被錯分為林地，工礦用地與建筑用地的時序變化曲線相似，也存在一定程度的誤分。

圖4 不同分類方法的分類結果對比Fig.4 Comparison of classification results of different classification methods

面向對象時序決策樹分類法不僅能考慮地形差異，按光譜、紋理信息將各地類進行聚合，還能識別耕地與其他地類的時序特征差異性，減少光譜相近像元的混分現(xiàn)象。 通過計算混淆矩陣進行精度評價，各分類方法的精度統(tǒng)計見表2。

表2 分類結果精度評定Table 2 Accuracy assessment of classification results

表2 中，面向對象決策樹法分類的總體精度及其Kappa 系數(shù)和耕地分類精度均高于監(jiān)督分類及決策樹分類方法，該方法適用于黃土礦區(qū)耕地變化信息的提取。

3.2 研究區(qū)耕地變化時空特征分析

以各影響區(qū)為統(tǒng)計單元，獲取年際耕地面積分布，見表3。 采用相對變化率A定量分析采動影響下耕地變化的空間差異［16］：

式中，sb為研究末期耕地面積；sa為研究初期耕地面積；T為研究時長。

由表3 可知，以2005 年礦井投產之前為基準數(shù)，截止2019 年，非開采影響區(qū)耕地面積基本保持穩(wěn)定，相對變化僅為0.45%。 開采影響區(qū)（包括間接影響區(qū)、沉陷區(qū)、采空區(qū)）內耕地面積變化大致可以分成，“減少—增加—減少”3 個階段。 第1 階段為2005—2013 年，是耕地面積快速減少階段，這是因為煤礦投產初期，采煤沉陷破壞了礦區(qū)原有的生態(tài)環(huán)境平衡狀態(tài)，此階段土地結構受損嚴重，耕地種植環(huán)境受到破壞。 第2 階段為2013—2015年，開采影響區(qū)耕地面積相對增加，經調查發(fā)現(xiàn)此階段政府出臺了各種耕地保護與生態(tài)補償政策，加之企業(yè)實施沉陷區(qū)治理工程，使得沉陷影響區(qū)耕地面積有所增加，但由于采煤沉陷效應的累積影響，開采影響區(qū)耕地面積減少量仍遠大于非影響區(qū)。 第3 階段為2015—2019 年，開采影響區(qū)耕地面積繼續(xù)減少。

表3 礦區(qū)耕地面積統(tǒng)計Table 3 Statistical data of cultivated land in coal mine

由于開采沉陷效應的滯后影響及相鄰工作面開采沉陷在空間上的疊加效應，此階段地表沉陷區(qū)范圍并非線性增加，從而開采影響區(qū)的耕地面積下降幅度較第1 階段有所減緩。

由表3 可知，采空區(qū)、沉陷區(qū)、間接影響區(qū)及非影響區(qū)耕地面積的相對變化率分別為－7.73%，－4.96%，－4.13%，－0.45%。 其中采空區(qū)耕地面積減少幅度最大，而非影響區(qū)耕地面積基本保持不變。分年度統(tǒng)計各影響區(qū)耕地面積變化率，如圖5 所示。采空區(qū)耕地面積變化率最大，距離采空區(qū)越遠時，其耕地面積變化率呈現(xiàn)遞減趨勢，這表明耕地面積減少幅度與采動影響程度呈正相關。

圖5 各影響區(qū)耕地面積年際間相對變化率Fig.5 Relative change rate of cultivated land

3.3 采煤沉陷對耕地變化的影響分析

計算研究區(qū)相關地形因子，獲取高程、坡度分級圖（圖6）并與土地利用分類結果進行疊加，可以發(fā)現(xiàn)研究區(qū)地處黃土丘陵溝壑區(qū)，地形復雜，耕地主要分布于海拔較高、地形平坦的黃土臺塬地區(qū)，少量位于臺塬與丘陵交界處坡度較緩的地帶，海拔低、地勢陡峭處，耕種難度較大，土地利用類型以林地、草地為主。

圖6 耕地利用地形空間分布Fig.6 Spatial distribution map of cultivated land use topography

基于分類結果，利用二元掩膜獲取年度耕地空間分布情況，如圖7 所示（圖7a—7g 為各年份耕地面積空間分布；圖7h 為耕地減少面積空間分布圖）。 由圖7 可以看出，不同年份的耕地在空間上的分布基本一致，但總覆蓋面積在逐年減少，且地形平坦地區(qū)耕地減少幅度小于坡度發(fā)生變化處。

圖7 年度耕地面積及減少分布Fig.7 Annual arable land area and spatial distribution of reduced area

為進一步分析地下開采與耕地面積變化的量化關系，采用GIS 空間技術，對年度耕地分布進行疊加分析，結合采煤擾動影響分區(qū)，獲取耕地減少面積的時空分布。

受采礦擾動直接影響，耕地面積減少的主要區(qū)域與相應的采空區(qū)及沉陷區(qū)域在時間和空間上存在明顯的對應關系，通過統(tǒng)計不同開采時段兩區(qū)域的耕地面積減少情況，發(fā)現(xiàn)其減少幅度基本相當，減少區(qū)域主要分布在采煤邊界兩側的沉陷區(qū)內，部分采煤工作面所對應地表覆蓋類型以林草地為主，所以即使位于擾動程度最大的影響區(qū)耕地面積減少幅度反而較低。

由于地表下沉存在超前影響和滯后效應，未開采工作面和已停采工作面地表仍會產生下沉變形，導致沉陷區(qū)及間接影響區(qū)耕地面積繼續(xù)減少。 而非影響區(qū)內部分耕地轉化為建設用地導致區(qū)域耕地面積發(fā)生減少。

繪出煤炭開采量與耕地減少面積隨年度變化的 關系曲線，如圖8 所示。 礦井投產初期，開采沉陷引起環(huán)境發(fā)生突變，導致耕地面積急劇減少。 隨著采空區(qū)范圍和煤炭開采量的線性增長，耕地面積減少的幅度則有所降低。

圖8 耕地累計減少面積與煤炭開采量相關Fig.8 Cumulative cultivated land change and coal mining volume correlation curve

4 討 論

4.1 采煤沉陷對土地破壞的影響

相較露天采煤，地下井工采煤方式對土地的影響看似良性，實則破壞巨大。 地下采煤引起的地表不均勻下沉和變形會破壞土地的平整性，影響正常的糧食種植活動，采煤沉陷使平坦的耕地因坡度變大而轉變?yōu)槠赂兀料葸吔缟戏揭蛩阶冃渭扒首冃渭卸纬傻乇砹芽p則導致耕地破碎程度加深，增加耕作難度。

4.2 采煤沉陷對耕地質量的影響

煤炭開采對礦區(qū)耕地質量上的影響主要有兩個方面，一是耕地土壤理化性質的改變，二是農作物品質及產量的降低。 沉陷邊界上方因水平變形及曲率變形集中而形成地表裂縫的形成，改變了土壤滲透性和孔隙度［17，19］，加重耕地水土和有機質流失。 由開采沉陷引起的土壤結構改變使得土壤肥力重新分布、土壤養(yǎng)分降低，造成耕地退化和作物生產力下降［19－21］，使得棄耕或季節(jié)性棄耕現(xiàn)象日益嚴重，采煤過程中產生的廢料占用、土地污染等，進一步導致耕地數(shù)量及農作物產量下降。

通過冬小麥距平值評估開采沉陷對耕地作物長勢的影響，如圖9 所示。 2005 年冬小麥長勢最好，2007 年小麥較2005 年長勢有所下降，但仍高于平均長勢。 小麥長勢低于平均值的情況首次在2009年采空區(qū)內出現(xiàn)，自此開采影響區(qū)的小麥長勢均低于平均值。 根據距平值變化幅度看采空區(qū)內農作物長勢受采煤擾動影響最大，沉陷區(qū)次之，間接影響區(qū)在采煤初期影響較大，而后減弱。 而非影響區(qū)農作物長勢保持平穩(wěn)波動，各階段受采煤擾動影響很小。

圖9 冬小麥NDVI 距平值Fig.9 The curve of NDVI anomaly

4.3 采煤過程中耕地恢復情況

隨著開采工作面的持續(xù)推進，在采空區(qū)面積保持線性增加的情況下，由于開采沉陷的疊加影響及自然力的恢復作用，沉陷區(qū)及間接影響區(qū)的采煤擾動影響趨于平緩，使得耕地面積變化趨于平緩。 據調查，2013—2015 年政府出臺了各種耕地保護以及煤礦開采區(qū)農業(yè)與生態(tài)補償?shù)南嚓P政策，加強了對礦區(qū)耕地的保護，加之企業(yè)實施沉陷區(qū)治理工程，使得沉陷影響區(qū)耕地面積反而有所增加，耕地質量甚至有所好轉。 從遙感時序數(shù)據反演的作物長勢指數(shù)變化的時間點上也證實了上述政策效應，這說明人為政策干預在短期內具有改善礦區(qū)生態(tài)環(huán)境質量的效果。 2015—2019 年，隨著地下開采區(qū)域的不斷增加，沉陷區(qū)及間接影響區(qū)在空間上出現(xiàn)重疊，伴隨著耕地保護政策和企業(yè)土地復墾工作的持續(xù)推進，此階段耕地面積雖然整體下降，但變化幅度小于煤礦投產的前期階段。

5 結 論

1）利用遙感時序數(shù)據可定量反演黃土高原礦區(qū)農耕地變化的時空特征。 利用ESTARFM 模型克服了遙感數(shù)據時空分辨率不足的局限性，融合后的時序影像與真實影像間的相關系數(shù)達0.946 1，具有很強的一致性；采用結合物候、光譜特征的面向對象決策樹法進行影像分類，能充分利用紋理及空間特征，避免產生“椒鹽現(xiàn)象”，提高了耕地面積提取的精度。 與傳統(tǒng)的監(jiān)督分類法和決策樹分類法相比，本文方法進行土地分類的總體精度分別提升了20.34%，7.78%，所提取耕地的精度達到93.3%。

2）研究區(qū)2005—2019 年耕地面積呈整體減少趨勢，不同時期耕地面積變化幅度有所不同，而政策與人工干預對耕地變化有明顯影響。 研究期內耕地總面積減少722 hm2，其中煤礦投產初期（2006—2009 年）減少幅度最大，達到329 hm2，隨著采空區(qū)面積的線性增大，各影響區(qū)內耕地面積減小幅度則趨于平緩。 2013—2015 年，受政府推出耕地保護和生態(tài)補償政策及企業(yè)開展土地復墾影響，研究區(qū)耕地面積有所增加。

3）不同影響區(qū)內耕地面積變化及作物長勢與采煤沉陷擾動程度相關，且呈現(xiàn)出明顯的空間分異特征。 各影響區(qū)耕地面積的相對變化率從采空區(qū)（－7.73%）、沉陷區(qū)（－4.96%）、間接影響區(qū)（－4.13%）至非影響區(qū)（－0.45%）依次遞減采空區(qū)內耕地面積在耕地總面積中的占比由2005 年煤礦投產前的32.86%減少至2019 年的19.87%，在開采工作面邊界附近地表變形劇烈區(qū)區(qū)域的農耕地變化最為顯著，揭示出黃土高原采煤沉陷區(qū)農耕地面積和作物長勢具有復雜的時空演化特征。