基于物元可拓法的高潛水位煤礦區土地損毀程度評價

李 超，崔瑞豪，張 曦，張曦元，馬晉明，王培俊

（1.永煤集團復墾公司，河南 永城 476600；2.中國礦業大學環境與測繪學院，江蘇 徐州 221116；3.西北農林科技大學資源環境學院，陜西 楊凌 712100；4.中國礦業大學公共管理學院（應急管理學院），江蘇 徐州 221116）

中國煤炭資源95%以上采用井工開采方式，煤炭開采造成大面積的土地下沉，并導致房屋坍塌、基礎設施受損、水土流失加重、生態環境惡化，農田遭到了嚴重的破壞，農民的生產生活受到嚴重的影響[1-2]。高潛水位采煤沉陷區具有地勢平坦、潛水位埋深小、可采煤層數量多、煤層厚度大和地表下沉系數大等特點，致使地表容易塌陷積水，農田減產甚至絕產，土地損毀嚴重[3-4]。因此，為了更好地保護生態環境，提高復墾效率，相關學者提出了“邊采邊復”技術[5-6]。由于“邊采邊復”技術是在地面沉陷前或沉陷過程中采取的復墾措施，從而促進了“末端治理”向“源頭和過程控制與治理”的轉變[7]，極大地提高了土地復墾率，減少了對土地的破壞。

土地損毀程度評價是復墾方案編制的重要組成部分，為土地復墾提供基礎數據，其結果直接影響復墾的效果及受損農民補償標準的制定[8-9]。國外的煤炭資源大多采用露天礦開采方式，不會產生采煤沉陷地。因此，國外對于采煤沉陷區土地損毀程度評價的研究很少[10]。隨著中國對土地復墾工作重視程度不斷提升，對采煤沉陷區土地損毀程度評價的研究也越來越多[11-12]。部分學者對評價指標體系及其分級等進行了探討，常用的評價方法有模糊綜合評價法、極限條件法、指數法等[13-16]，并利用GIS 手段優化評價過程[17]，實現評價結果的可視化。但如何客觀有效地構建采煤沉陷區土地損毀程度評價指標體系，并完善其評價方法是重中之重。程琳琳等[18]使用可拓評價方法，建立高潛水位采煤沉陷區土地損毀評價指標的形式化模型，通過分析事物變換的規律，降低了多因素評價中的主觀片面性，大大提高了評價結果真實性。該方法已廣泛應用于其他領域的評價[19]，陳秋計等[20]初步驗證了該方法在礦區土地復墾適宜性評價、采煤沉陷區土地損毀程度評價中的可行性，并取得了良好的效果。

在多數采煤沉陷區土地損毀程度評價中，積水情況屬于必須要考慮的內容，多數土地損毀評價中對當地積水情況不清楚，對積水范圍和積水時間不明確，這嚴重影響了土地損毀程度評價的準確性。因此，借鑒HE 等[21]、趙艷玲等[22]對華東地區沉陷區積水的識別方法，使用Landsat 8 遙感數據生成年內水頻率指數（AWFI），以代表沉陷區的年度積水量，并定量描述積水情況，作為評價基礎數據，使得評價結果更為準確。

目前，礦區土地損毀程度評價研究蓬勃發展，但對于東部高潛水位平原礦區的土地損毀評價研究較少。以永城市的城郊礦區和新橋礦區為例，在梳理已有研究的基礎上，構建評價指標體系，劃分評價單元，采用熵權法確定指標權重，針對高潛水位礦區土地損毀的特點，將極限條件法和可拓評價方法相結合，對過去數年內開采造成的沉陷區土地損毀程度進行評價，并使用無人機影像數據繪制植被長勢進行評價驗證[23-24]，為該區域土地復墾方案編制、受損農民補償標準的制定等提供理論支持，為其他高潛水位采煤沉陷區土地損毀程度的評價提供借鑒。

1 研究區概況及數據來源

1.1 研究區概況

永城市屬暖溫帶半濕潤半干旱大陸性季風氣候，夏季炎熱多雨，冬季干燥寒冷，四季氣候變化分明。永城礦區整體土地約70%屬于耕地，易積水形成湖泊，耕地破壞面積大，沉陷區鹽漬化趨勢加重。

城郊礦區位于黃淮沖積平原的東部，地勢較為平坦開闊，總體地勢西北高，東南低，地面標高32～34 m，相對高差2 m 左右。據永城市氣象站資料，降雨量多集中于6 月—8 月，占全年降雨量的50%；夏季多東南-南風，冬季多西北-北風，年平均風速3.4 m/s，最大風速20 m/s。井田內地表水系多為縱橫交錯的人工排澇溝渠，最大河流為沱河，由西北流向東南，橫貫全井田。該河流為季節性河流，雨季流量增大，旱季流量劇減。新橋礦區處于華北平原與黃淮沖積平原的交接部位，地勢平坦，自西向東略有傾斜。地面標高31～33 m，相對高差2 m 左右。降雨、水系等情況與城郊礦區類似。城郊礦區和新橋礦區地理區位和土地利用分類情況如圖1 所示。

圖1 城郊礦區與新橋礦區地理區位和土地利用分類情況Fig.1 Geographical location and land use classification of Chengjiao Mining Area and Xinqiao Mining Area

1.2 數據來源

數據資料主要包括基礎礦山數據、實地調研數據以及遙感影像數據：①商丘市永城市的清華2017土地利用數據，主要獲取農田土地利用位置，以確定試驗田和采樣點在大致位置；②商丘市永城市2020年地形圖數據；③Landsat 8 TM 影像用以繪制年積水頻率圖；④無人機影像數據；⑤煤炭資源開發利用方案、土地復墾方案、環境監測報告等。

使用的遙感傳感器是基于大疆PSDK 定制開發的MS600 PRO 無人機多光譜相機，相機鏡頭可以獲取藍光波段（波長450 nm、帶寬27 nm）、綠光波段（波長555 nm、帶寬27 nm）、紅光波段（波長660 nm、帶寬22 nm）、紅邊波段（波長710 nm、帶寬10 nm）和近紅外波段（波長840 nm、帶寬15 nm）等波段影像。搭載MS600 PRO 的遙感平臺為四旋翼無人機大疆M210-RTK，航向重疊度及旁向重疊度均設置為75%，航高90 m，航速1.9 m/s，三次飛行累計航點數81 個，航線長度約13.6 km，飛行范圍約133 749 m2，預計飛行時間67 min，共計獲取像片8 388 張，覆蓋蓮藕、芡實、水稻的主要種植范圍。

空間數據包括土地利用數據、DEM（數字高程模型）、礦區分布圖及Landsat 衛星影像數據等。表1展示了數據源、分辨率和用途。通過裁剪、光柵化和重采樣對這些數據進行預處理，并對收集統計到的數據進行梳理檢查，確保單位的規范統一；在QGIS 中實現各專題地圖的空間投影坐標系一致。

表1 數據源、分辨率及用途Table 1 Data source, resolution and purpose

2 土地損毀評價

2.1 評價過程

2.1.1 評價指標的選取

根據《農用地質量分等規程》（GB/T 28407—2012）、《土地復墾方案編制規程》（TD/T 1031.3—2011）、相關研究文獻和土地復墾方案案例，初步篩選了地表形變的7 個備選指標，即積水情況、下沉深度、裂縫密度、傾斜變形、曲率、水平移動和水平變形等。結合實地考察，剔除多項爭議指標。研究區待評價單元主要種植水稻、芡實和蓮藕等水生植物，需要從時間角度考慮積水情況，因此加入了時間因素層；東部高潛水位礦區的裂縫與西部露天煤礦的地裂縫情況不同，常出現在水下，且較為細小難以觀測，因此去掉裂縫密度指標；但又因為裂縫是實際存在的，所以需要水平位移指標進行表示。最終確定的高潛水位平原型采煤沉陷區土地損毀程度評價指標體系見表2。

表2 采煤沉陷土地損毀程度評價指標體系Table 2 Evaluation index system of damage degree of coal mining subsidence land

2.1.2 指標因子獲取

積水時間指標主要通過Landsat 數據生成計算獲取，并在實地考察時驗證準確性。使用年內水頻率指數表示年積水情況，在GEE 平臺上通過歸一化差分植被指數（NDVI）和增強植被指數（EVI）來降低植被對水體分類的干擾，水體分類主要是使用歸一化修正差分水指數（MNDWI）[25]，具體計算見式（1）～式（3）。一個像素點滿足MNDWI>NDVI或MNDWI>EVI和EVI<0.1 條件時將被歸類為水體，通過計算像素點有水的影像數量與全年的影像數量的值來表示積水情況。相關研究表明水體分為季節性水（0.3≤AWFI<0.75）和永久水（AWFI≥0.75），（0，0.3）范圍為未積水或無效數據區。積水結果如圖2 所示。

圖2 研究區積水分布圖Fig.2 Water distribution map of study area

式 中： ρblue、 ρgreen、 ρred、 ρnir和 ρswir分別為Landsat 8 影像的藍色波段（450～520 nm）、綠色波段（520～600 nm）、紅色波段（630～690 nm）、近紅外波段（760～900 nm）和短波紅外波段（1 550～1 750 nm）。

地表形變的空間指標獲取主要通過MSPS 軟件獲取，如下沉深度、傾斜變形、水平變形和曲率。根據城郊礦區和新橋礦區實測地表移動參數，結合開采實踐及上覆巖層物理力學特性，確定礦區內工作面采用概率積分法進行地表移動變形預計的參數為：下沉系數q=0.70、水平移動系數b=0.33、主要影響角正切tgβ=2.0、最大下沉角θ=90°?0.6α、拐點偏移距s=0.05H。地表移動變形預計的參數結果見表3。

表3 工作面開采地表移動變形各參數最大值Table 3 Maximum parameters of surface movement and deformation in working face mining

2.1.3 指標因子分級

1）劃分評價等級。依據積水、地形變化和土壤肥力等方面差異，將高潛水位采煤沉陷區土地損毀程度評價的結果統一分為三個等級（表4），即輕度損毀、中度損毀和重度損毀，以便于該區域復墾標準及受損農民補償標準的制定及評價結果的判斷。

表4 損毀程度等級分級Table 4 Classification of damage degree

2）因子的分級標準。依據《農用地質量分等規程》（GB/T 28407—2012）、《土地復墾方案編制規程》（TD/T 1031.3—2011），結合已有研究的土地損毀程度的分級標準，根據研究區的實際情況，得到各因子的損毀程度等級劃分標準，見表5。

表5 研究區損毀程度等級劃分Table 5 Classification of damage degree in study area

2.1.4 確定指標權重

評價指標權重的確定是建模評價的關鍵環節，權重賦予合理與否直接關系到評價結果的可信度。為有效避免指標權重賦值的主觀隨意性，使用熵權法對各評價指標賦予權重，具體步驟見式（4）～式（6）。

熵值是一種物理計量單位，熵值越大說明數據越混亂，攜帶的信息越少，效用值越小，因而權重也越小。由表6 可知，積水時間熵值最小、權重最高符合一般土地損毀評價中的賦值規律。下沉值的信息熵值最大、信息效用值最小，主要是因為評價單元的下沉值較為接近。

表6 熵權值結果Table 6 The results of entropy weight

2.1.5 采用可拓理論評價損毀程度

1）數據預處理。評價指標對損毀程度的影響有正效應或負效應，即部分指標因子的數值越小損毀等級越高，部分指標數值越大損毀等級越高（如下沉深度），無法直接進行比較，因此將各評價指標進行歸一化無量綱處理，定量化的指標可以根據式（7）進行歸一化處理。

式中：X為指標的邊界值；xˉ為對應指標無量綱化后的值；Xmax和Xmin為該指標的最大值和最小值。定性指標需結合經驗以及當地的實際資料來賦值，最后按照定量指標的方式進行歸一化處理。

2）確定評價指標的節域及經典域。系統節域物元模型見式（8）。

式中：P為采煤沉陷區土地損毀等級評價的全部集合，即評價指標，是對應于Ci（i=1,2,···,n）的全部取值范圍，也就是P的節域系統的經典域物元模型。

各指標的分級（j=1,2,3）分別表示輕度損毀、中度損毀和重度損毀，計算見式（9）。

式中，VPi為對應的指標Ci（i=1,2,···,n）的取值范圍，經典域為 ?aPn,bPn?。

3）建立待評價對象的物元模型和關聯函數，計算關聯度。以物元表示分析得到的評價數據。式（10）為評價單元對應因子的具體數值，是待評價對象。

關聯度是事物間或者因素間關聯性的量度，這里描述的是各評價指標對于不同損毀程度的影響水平，計算公式見式（11）。

式中：Kji為第i個評價因子隸屬于評價等級j的關聯度； ρ（vi,vo ji） 為 點vi與區間vo ji的距；的模。其中， ρ（vi,vo ji）可以通過式（12）確定。

式中：y為實域（–∞，+∞）上的一點；Y=（a，b）為實域上的區間。

4）確定損毀等級。損毀等級計算見式（13）。

式中：Kj（Po）為評價單元對不同土地損毀程度等級的符合程度。依據式（14）確定土地損毀程度等級。

2.2 評價單元的劃分

評價單元是具有特定土地特性和土地質量的土地評價和制圖的基本單元。依據《農用地分等定級規程》，參考已有研究成果，得出常用的評價單元劃分方法有圖斑法、疊置法、地塊法及動態網格法。首先選取疊置法將永城市土地利用現狀圖、DEM、開采沉陷預計的沉陷范圍和礦權邊界圖進行疊加，其次在此基礎上提取所有建筑用地、水體和農用地作為評價對象，用以表示整個采煤沉陷區地物類別，再次將其與年內水頻率指數圖疊加并使用緩沖區分析得出34 個評價單元，如圖3（a）和圖3（b）所示。水稻玉米田評價單元編號為1～10、芡實評價單元編號為11～21、蓮藕池評價單元編號為22～34。最后將獲取的評價單元根據位置狀況相互連接繪制成圖。

圖3 土地利用圖和年內水頻率圖疊加獲取的評價單元圖Fig.3 Evaluation unit diagram obtained by superposition of land use map and annual water frequency map

3 結果與討論

3.1 損毀程度結果

芡實地、水稻玉米田和蓮藕池各評價單元相應指標的實測值見表7，土地損毀情況如圖4～圖6 所示。通過上述計算得出芡實地評價單元中11 點位、15 點位、16 點位、19 點位、20 點位、21 點位為中度損毀，其他評價單元為輕度損毀。芡實地總面積為444.7 hm2，中度損毀面積為135.7 hm2，占比30.51%，即有接近三分之一的土地會有季節性積水，地形產生形變，土壤肥力下降。如圖4（b）所示芡實地評價單元中中度損毀分為三塊，分別是11 點位，位于西北處邊界，面積最小只有9.0 hm2；15 點位和16 點位位置接近所以合并單元，其面積為48.6 hm2，位于西南邊界處；19 點位～21 點位位于東部，面積最大為78.0 hm2。如圖5（b）所示水稻玉米田整體屬于輕度損毀狀態，總面積為747.2 hm2，屬于三個研究區中土地損毀情況最輕的等級，即無積水或短期積水，地形未發生明顯改變，土壤肥力無影響。蓮藕池總損毀面積為562.6 hm2，中度損毀面積為151.8 hm2，占比26.98%。如圖6（b）所示23 點位、24 點位、27 點位、28 點位、31 點位、34 點位為蓮藕池的中度損毀地區。其中，23 點位位于蓮藕池的東北角，損毀面積為51.0 hm2；24 點位、27 點位和34 點位位于其西南角，總面積為45.7 hm2。

表7 研究區評價指標實測值（部分）Table 7 Measured value of evaluation index in study area（part）

圖4 芡實地土地損毀總體情況和放大圖Fig.4 Overall situation and enlarged map of land damage in gorgon nut land

圖5 水稻玉米田土地損毀總體情況和放大圖Fig.5 Overall situation and enlarged map of land damage in rice corn field

圖6 蓮藕池土地損毀總體情況和放大圖Fig.6 Overall situation and enlarged map of land damage in lotus root pond

3.2 結果驗證

研究中的土地損毀評價主要通過積水和地形的改變進行判斷，因此使用土地肥力條件進行驗證。土壤肥力的驗證主要體現在農作物長勢上，圖4～圖6 中無人機影像位置與損毀程度評價不同程度類型同時重疊的僅有新橋礦區的蓮藕池，因此以其為例進行結果驗證。將無人機多光譜影像進行預處理和指數運算，繪制代表蓮藕長勢的NDVI圖，如圖7所示。輕度損毀地區即圖7（a）和圖7（b）中圓點位置，放大后可明顯觀察出該區域大部分蓮藕處于長勢良好的狀態。中度損毀地區即圖7（a）和圖7（c）中菱形位置，長勢呈現明顯的淺色即長勢較差。該結果較符合中度損毀定義中土壤肥力下降的直觀表現，表明土地損毀程度評價結果較為準確合理。

圖7 蓮藕長勢NDVI 反演結果Fig.7 NDVI inversion results of lotus root growth

3.3 討論

積水情況是土地損毀評價的重要影響因素，但不是唯一影響因素。對高潛水位礦區土地損毀程度評價，通常的做法是通過估計積水的范圍與程度來確定損毀程度。但該經驗方法既具有主觀片面性與不準確性，又不能定量表示損毀范圍，因此需要引入新的評價方法。因此將基于遙感影像的積水時間數據引入評價體系中，從而構建全新的高潛水位采煤沉陷區土地損毀程度評價指標體系，使用基于極限條件法與可拓物元法的評價方法，以及基于熵權法的權重確定方法，提出了采煤沉陷土地損毀指標體系，使評價結果更為客觀、科學和合理。將可拓方法應用于礦區土地損毀評價，把各評價因子定量化，盡量減少人為因素的影響，且使用物元可拓模型評價出的土地損毀程度與依據積水情況的經驗做法已經有了直觀上的差異，克服了多因素識別評價中的主觀片面性，大大提高了評價結果真實性、科學性和合理程度。

4 結 論

通過構建的高潛水位平原型采煤沉陷區土地損毀程度評價指標體系方法，對城郊礦區和新橋礦區進行土地損毀程度評價，得到以下結論。

1）研究區整體土地損毀程度較輕，其中，蓮藕池中度損毀面積最大，為151.8 hm2，占蓮藕池損毀總面積的26.98%；水稻玉米田損毀程度最低，整體屬于輕度損毀程度。

2）評價指標體系方法的評價結果與單憑積水情況的評價方法有明顯不同，引入遙感等方法后，使得研究區損毀情況更加直觀、客觀、科學、合理。

3）通過無人機影像繪制植被長勢可以直觀地區分出長勢的不同，進而初步得出土壤肥力不同的區域，即中度損毀區域植被長勢差于輕度損毀區域。