夏嘉南，李根生，卞正富，雷少剛，3，宮傳剛，李 恒

（1.中國礦業(yè)大學(xué) 礦山生態(tài)修復(fù)教育部工程研究中心，江蘇 徐州 221116；2.中國礦業(yè)大學(xué) 公共管理學(xué)院，江蘇 徐州 221116；3. 山東省煤田地質(zhì)局第二勘探隊(duì) 山東省采煤塌陷地和采空區(qū)治理工程研究中心，山東 濟(jì)寧 272000）

0 引 言

內(nèi)排土場作為露天礦區(qū)重要的組成部分，其占地面積隨采復(fù)工作面推進(jìn)不斷增大[1]。在內(nèi)排過程中，傳統(tǒng)地貌重塑方法往往將原始微起伏地貌景觀轉(zhuǎn)變?yōu)槠脚_(tái)和邊坡相間的“梯田式”人工堆疊景觀，使原始地表要素及水系特征消失[2]。雖然在后續(xù)復(fù)墾設(shè)計(jì)中，內(nèi)排土場內(nèi)布設(shè)輸水管網(wǎng)以滿足復(fù)墾植被生長需求，但整體上未能實(shí)現(xiàn)其與周邊地貌融合，易引發(fā)土地退化等問題[3-4]。因此，需嘗試在采復(fù)過程中保留原始自然地貌特征以破解上述困境。

近年來，不少專家學(xué)者就內(nèi)排土場規(guī)劃及內(nèi)排方案開展了研究。例如，王東等[5]以內(nèi)排土場空間利用最大化為目標(biāo)，基于采場工作線及邊坡角、內(nèi)排土場邊坡角、基底傾角對排土線的影響規(guī)律，構(gòu)建了采場工作幫、內(nèi)排土場、端幫和基底的空間幾何模型；栗嘉彬[6]以技術(shù)最佳適用性為目標(biāo)，根據(jù)物料用途，構(gòu)建了露天礦表土最優(yōu)采集周期、最佳覆土厚度等求解模型；OGGERI等[7]以最小復(fù)填成本為目標(biāo)，根據(jù)場地、運(yùn)輸?shù)戎T因素，提出了多情境下露天礦采復(fù)處置方法。然而上述研究中，未充分考慮重塑地貌結(jié)果與周邊自然地貌起伏貼合情況，導(dǎo)致內(nèi)排土場往往被重塑為平臺(tái)加斜坡的“簡單梯田式”形態(tài)，人工地貌形態(tài)明顯，易引發(fā)嚴(yán)重的景觀破碎化問題[8-10]。相較于人工規(guī)則地貌，自然地貌作為長期演變的結(jié)果，具有高度的區(qū)域適宜性與穩(wěn)定性，可作為修復(fù)區(qū)地貌重塑的參考對象[3,10-11]。且國外實(shí)踐證明[12-13]，基于自然的地貌重塑能更好地適應(yīng)當(dāng)?shù)厮畾鈼l件，在有效提高區(qū)域抗水蝕能力的同時(shí)僅需要較少的維護(hù)費(fèi)用[14]。因此，迫切需要一種貼合周邊自然地貌的內(nèi)排土場重塑方法，用于解決內(nèi)排土場設(shè)計(jì)結(jié)果與周邊自然地貌間貼合的問題[15]。

基于此，筆者以新疆黑山露天礦為例，提出一種采排復(fù)一體下基于調(diào)整曲面的內(nèi)排土場近自然地貌重塑模型，并結(jié)合地理空間分析技術(shù)與MATLAB軟件，構(gòu)建內(nèi)排土場近自然地貌重塑結(jié)果。并通過CLiDE景觀演化模型，模擬了研究區(qū)10 a間原始地貌、傳統(tǒng)設(shè)計(jì)地貌及近自然設(shè)計(jì)地貌三者的土壤水蝕情況，對比評價(jià)了近自然設(shè)計(jì)地貌抗水蝕效果，為相同及相似地區(qū)內(nèi)排土場地貌形態(tài)設(shè)計(jì)提供借鑒。

1 研究區(qū)概況與數(shù)據(jù)來源

1.1 研究區(qū)概況

研究區(qū)位于新疆托克遜縣西北約90 km處，地理坐標(biāo)：87°26′04″～87°29′40″E，43°13′48″～43°14′01″N，周邊地貌以山間谷地為主；年均降水量約160.3 mm，年均氣溫約1.6 ℃，冬寒夏涼，屬大陸干旱及高寒氣候。區(qū)內(nèi)表土以第四系為主，野生植被類型為典型草原植被。

依據(jù)地質(zhì)勘探資料，選取吉布拉克煤礦探礦權(quán)范圍以南為研究區(qū)（圖1）。區(qū)域東西長約4 874 m、南北長約1 000 m。其中，預(yù)設(shè)采坑邊界東西長約812 m、南北長約1 000 m，采礦作業(yè)由東至西；區(qū)域可采煤層分別為6、7、8、9、11、12-1、12-2、13-1和13-2號煤層，其中13-2號煤層為全區(qū)可露采、其余為局部可露采。

圖1 研究區(qū)示意Fig.1 Schematic of the study area

1.2 數(shù)據(jù)來源

研究區(qū)采前自然原始數(shù)字高程數(shù)據(jù)（DEM）來源于地理空間數(shù)據(jù)云，拍攝時(shí)間2011年12月，條帶號87，行編號43，空間分辨率30 m×30 m；研究區(qū)煤層底板、頂板標(biāo)高及水平坐標(biāo)數(shù)據(jù)來源于《新疆托克遜縣黑山礦區(qū)黑山露天煤礦資源儲(chǔ)量核實(shí)報(bào)告》。此外，基于上述報(bào)告通過ArcGIS克里金插值模塊獲取煤層底板及頂板空間分布數(shù)據(jù)；內(nèi)排土場采、復(fù)子區(qū)空間分布數(shù)據(jù)基于內(nèi)排土場采復(fù)子區(qū)位置識別方法獲??；復(fù)填子區(qū)內(nèi)的可用土方量數(shù)據(jù)通過可用土方量計(jì)算式求解獲??；復(fù)填區(qū)地表調(diào)整曲面數(shù)據(jù)通過地表調(diào)整曲面預(yù)構(gòu)建模型獲??；復(fù)填子區(qū)地表近自然DEM通過MATLAB軟件基于動(dòng)態(tài)土方量一致、且區(qū)域斜率較緩的原則定向篩選不同B樣條控制點(diǎn)位置下的地表調(diào)整曲面獲得；研究區(qū)自然原始DEM數(shù)據(jù)通過ArcGIS空間疊加提取獲得；研究區(qū)內(nèi)排土場傳統(tǒng)設(shè)計(jì)DEM數(shù)據(jù)基于土方量動(dòng)態(tài)平衡下“30°斜坡加平臺(tái)”的設(shè)計(jì)思路獲??；最后運(yùn)用CLiDE演化模型與MATLAB軟件，分別模擬并統(tǒng)計(jì)研究區(qū)原始DEM、傳統(tǒng)設(shè)計(jì)DEM及近自然設(shè)計(jì)DEM三者10 a的土壤水蝕量，以評估近自然整形結(jié)果的抗侵蝕能力。

2 研究方法

2.1 內(nèi)排土場近自然地貌重塑模型構(gòu)建

2.1.1 采復(fù)子區(qū)空間位置識別

依據(jù)露天礦采復(fù)周期，獲取圖2所示各周期下開采子區(qū)amn和復(fù)填子區(qū)afn水平投影空間位置。

圖2 露天礦采復(fù)周期下開采及復(fù)填子區(qū)空間位置提取示意Fig.2 Schematic of extraction of spatial location of the mining sub-area and the refilling sub-area under the open-pit mining-refill cycle

如圖2所示，依據(jù)露天礦實(shí)際采復(fù)周期，沿開采方向?qū)⒙短斓V區(qū)按采復(fù)周期時(shí)間先后順序劃分為n個(gè)子區(qū)。其中，開采子區(qū)和復(fù)填子區(qū)分別對應(yīng)區(qū)域amn和區(qū)域afn。在露天礦內(nèi)排構(gòu)建的采復(fù)過程中，采坑每前進(jìn)一個(gè)開采子區(qū)amn對應(yīng)修復(fù)后側(cè)一個(gè)復(fù)填子區(qū)afn，直至內(nèi)排過程結(jié)束，礦坑由區(qū)域Ab移至區(qū)域Aa。

2.1.2 可用土方量計(jì)算

為便于后續(xù)土方量計(jì)算，將開采及復(fù)填子區(qū)均細(xì)分為rowc行colc列的網(wǎng)格單元組（圖2）。在采復(fù)周期n中，開采子區(qū)amn可用于對應(yīng)復(fù)填子區(qū)afn的土方量，計(jì)算式如下：

（1）

式中:V（amn）為開采子區(qū)amn可用于對應(yīng)復(fù)填子區(qū)afn的土方量，其值是開采子區(qū)afn最下層可采礦層底板以上的總物質(zhì)體積減去區(qū)域內(nèi)所有可采礦層總物質(zhì)體積；i,j分別為開采及復(fù)填子區(qū)內(nèi)網(wǎng)格單元的行、列序號，取取值分別在1～rowc和1～colc之間；H（amnij）為開采子區(qū)amn網(wǎng)格單元（i,j）所在區(qū)域的采前地面標(biāo)高；hd（amnij）為開采子區(qū)amn網(wǎng)格單元（i,j）所在區(qū)域的最底層礦層底板標(biāo)高；ht（amnij）為開采子區(qū)amn網(wǎng)格單元（i,j）所在區(qū)域的可采礦層總厚度；S（amnij）為開采子區(qū)amn子單元（i,j）的水平投影面積，其余參數(shù)解釋同上。

2.1.3 復(fù)填區(qū)地表調(diào)整曲面預(yù)構(gòu)建

依據(jù)區(qū)域原有自然穩(wěn)定地貌，構(gòu)建圖3所示調(diào)整曲面模型，以盡可能保留復(fù)填子區(qū)afn原有地貌起伏特征。

圖3 地表調(diào)整曲面預(yù)構(gòu)建示意Fig.3 Schematic of pre-built adjustment surface

如圖3所示，將反映復(fù)填子區(qū)最終地貌設(shè)計(jì)形態(tài)較其原始自然地貌空間高程變化的曲面定義為地表調(diào)整曲面。為保證地表調(diào)整曲面能使復(fù)填子區(qū)表面與周邊自然地貌“無縫”融合，選取邊界處高程變化小，中間變化大，且與區(qū)域自然地貌相似的連續(xù)放緩曲面。參照“船體”構(gòu)建思路，采用連續(xù)性較優(yōu)的B樣條曲線[16]為“龍骨”，三角函數(shù)為側(cè)方“支架”，構(gòu)建復(fù)填子區(qū)afn所對應(yīng)的地表調(diào)整曲面。圖中虛線箭頭為開采方向；x,y,z三軸正方向分別對應(yīng)開采方向、水平垂直開采方向及豎直方向；fb和ft分別對應(yīng)B樣條曲線函數(shù)和三角函數(shù)。其中，B樣條曲線fb位于復(fù)填子區(qū)afn水平投影中心線位置，采用控制點(diǎn)構(gòu)建，其計(jì)算式如下：

fb（j）=Bezier（v,tb,ts,te）

（2）

（3）

其中:fb（j）為復(fù)填子區(qū)afn第j列的B樣條曲線值，定義域在0～1；Bezier為豎直平面上的2維B樣條曲線函數(shù)，形態(tài)由起始點(diǎn)0,y1、終止點(diǎn)1,yNC和中間NC-2個(gè)控制點(diǎn)共同控制；NC為Bezier函數(shù)中起點(diǎn)和終點(diǎn)在內(nèi)所有控制點(diǎn)的個(gè)數(shù)，取值為>2的整數(shù)；v為NC行2列的矩陣，由上至下行數(shù)據(jù)分別對應(yīng)起點(diǎn)、控制點(diǎn)1至控制點(diǎn)NC-2、終點(diǎn)的坐標(biāo)；tb為起點(diǎn)所對應(yīng)網(wǎng)格單元的列序號；ts為Bezier函數(shù)取值間隔，是終點(diǎn)與起點(diǎn)所對應(yīng)網(wǎng)格單元的列序號之差，取值為大于0的正數(shù)；te是終點(diǎn)所對應(yīng)網(wǎng)格單元的列序號；x1～xNC-2分別為控制點(diǎn)1至控制點(diǎn)NC-2在中心線上以起點(diǎn)為原點(diǎn)的水平投影距離，其中x1≤x2≤…≤xNC-2；y1～yNC分別是起點(diǎn)、控制點(diǎn)1至控制點(diǎn)NC-2、終點(diǎn)的預(yù)設(shè)高程值，定義域皆在0～1之間，其余參數(shù)解釋同上。在確定復(fù)填子區(qū)B樣條曲線函數(shù)形態(tài)的基礎(chǔ)上，通過三角函數(shù)確定地表調(diào)整曲面，其計(jì)算式如下：

（4）

式中:ft（i,j）為復(fù)填子區(qū)afn地表調(diào)整曲面網(wǎng)格單元（i,j）所對應(yīng)的預(yù)設(shè)高程值，取值在0～1；π為圓周率，其余參數(shù)解釋同上。

2.1.4 曲面土方量控制及坡度緩和優(yōu)化

為保證采復(fù)周期中可用土方量動(dòng)態(tài)平衡，通過豎直伸縮變化使地表調(diào)整曲面符合土方平衡要求，其計(jì)算式如下：

（5）

（6）

（7）

在此基礎(chǔ)上，為篩選出區(qū)域整體坡度較緩的復(fù)填子區(qū)地表高程模型，以提高其表土抗水蝕能力，構(gòu)建以網(wǎng)格單元與周邊網(wǎng)格單元間坡角為對象的區(qū)域坡度評分標(biāo)準(zhǔn)，其計(jì)算式如下：

（8）

其中：fp（a）為坡度夾角a所對應(yīng)的坡度評分值，其值為大于0的整數(shù)；au為正整數(shù)，其取值在2°～90°，其余參數(shù)解釋同上。在區(qū)域坡度評分標(biāo)準(zhǔn)構(gòu)建的基礎(chǔ)上，構(gòu)建區(qū)域坡度評分模型，其計(jì)算式如下：

（9）

式中：p（afn）為復(fù)填子區(qū)afn的區(qū)域坡度總評分；fpij（a）為子區(qū)afn內(nèi)網(wǎng)格單元（i,j）沿某一朝向的坡度評分值；∑fpij（a）為子單元（i,j）與周邊8個(gè)網(wǎng)格單元間坡角所對應(yīng)評分值的總和。

依據(jù)上述評分模型，通過MATLAB軟件，調(diào)整2.1.3節(jié)中B樣條曲線起點(diǎn)、控制點(diǎn)1至控制點(diǎn)NC-2、終點(diǎn)的坐標(biāo)位置，遍歷構(gòu)建一組形態(tài)各異的復(fù)填子區(qū)地表高程模型。并以區(qū)域坡度總評分最小為目標(biāo)，定向篩選復(fù)填子區(qū)Ht′求解結(jié)果，最終獲取原有自然形態(tài)表面下，土方量動(dòng)態(tài)平衡且坡度平緩的復(fù)填子區(qū)afn地表高程近自然設(shè)計(jì)形態(tài)。

2.2 設(shè)計(jì)結(jié)果抗水蝕能力評價(jià)

CLiDE模型是基于CAESAR模型開發(fā)的景觀地貌演化模型[17]。與傳統(tǒng)的侵蝕預(yù)測模型相比，其不僅能獲取侵蝕和沉積的時(shí)序數(shù)據(jù)，還能獲取最終發(fā)育的地貌形態(tài)。同時(shí)嵌入了Lisflood-FP流體力學(xué)相關(guān)模塊與地下水分布式水文模塊，可以模擬地表徑流-地貌-地下水三者間的交互作用和補(bǔ)給關(guān)系[18]，故極其適用于地貌修復(fù)效果評價(jià)[10]。為評價(jià)內(nèi)排土場設(shè)計(jì)結(jié)果的抗水蝕能力，將研究區(qū)內(nèi)排土場近自然設(shè)計(jì)結(jié)果為試驗(yàn)組，原始自然地貌及傳統(tǒng)復(fù)填方法下內(nèi)排土場設(shè)計(jì)結(jié)果為對照組，結(jié)合實(shí)地調(diào)查情況構(gòu)建CLiDE模型，模擬其10 a的土壤侵蝕過程。并基于MATLAB軟件統(tǒng)計(jì)各組地貌演變前后土方變化量，以評價(jià)其抗水蝕效果。

3 結(jié)果與分析

3.1 數(shù)據(jù)準(zhǔn)備結(jié)果及分析

研究區(qū)可采煤層共計(jì)9層，依據(jù)儲(chǔ)量報(bào)告，采用ArcGIS克里金插值模塊獲取研究區(qū)13-2號煤層底板標(biāo)高平面分布圖及其各煤層厚度平面分布圖（圖4）。為避免數(shù)據(jù)空間分辨率差異造成的干擾，將影像空間分辨率統(tǒng)一設(shè)置為30 m×30 m，數(shù)據(jù)形式均為162列32行的柵格影像。

圖4 13-2號煤層底板標(biāo)高平面分布與研究區(qū)各煤層厚度平面分布Fig.4 Plane distribution map of floor elevation of No. 13-2 coal seam and plane distribution map of the thickness of each coal seam

3.2 近自然設(shè)計(jì)地貌構(gòu)建結(jié)果及分析

依據(jù)2.1.1節(jié)采復(fù)子區(qū)空間位置識別方法，在30 m×30 m柵格影像數(shù)據(jù)下，按采復(fù)周期將研究區(qū)劃分為45個(gè)一一對應(yīng)的采復(fù)子區(qū)。其中，采復(fù)子區(qū)網(wǎng)格行數(shù)為32，列數(shù)為3，即復(fù)填子區(qū)和開采子區(qū)均設(shè)計(jì)為水平投影東西長90 m，南北長960 m的水平矩形區(qū)域，以使采復(fù)工作前后采坑的水平投影大小不變（圖5）。

圖5 研究區(qū)開采子區(qū)及復(fù)填子區(qū)空間位置示意Fig.5 Schematic of the spatial location of the mining sub-area and the refill sub-area

在此基礎(chǔ)上，依據(jù)2.1節(jié)所示內(nèi)排土場地表近自然設(shè)計(jì)方法，通過MATLAB軟件以0.01為控制點(diǎn)數(shù)值最小變化間距，膨脹系數(shù)k取均值1.1，求解研究區(qū)內(nèi)排土場地貌近自然設(shè)計(jì)DEM結(jié)果（圖6）。其中，考慮到礦區(qū)可采煤層非惟一（9層），且厚度及底板標(biāo)高存在差異，依據(jù)研究區(qū)三維地質(zhì)資料，以13-2號煤層的底板標(biāo)高為研究區(qū)煤層底板標(biāo)高，并采用空間疊加方法疊加各煤層厚度數(shù)據(jù)，獲取研究區(qū)煤層頂板標(biāo)高等效數(shù)據(jù)，以實(shí)現(xiàn)本研究內(nèi)排土場近自然地貌重塑模型于多礦層復(fù)雜情景下的應(yīng)用。

圖6 研究區(qū)內(nèi)排土場DEM圖Fig.6 DEM of the internal dump in the study area

如圖6所示，視覺效果上，研究區(qū)內(nèi)排土場近自然設(shè)計(jì)DEM與原始自然地貌及紋理特征較為相似，且高程值均在+2 000.26～+2 959.00 m；表土高程上，受土方平衡約束，近自然設(shè)計(jì)DEM較自然原始DEM存在差異，變化幅度在-205.90～50.03 m。一方面，受研究區(qū)“開采子區(qū)地表原始標(biāo)高”、“復(fù)填子區(qū)煤層底板標(biāo)高”和“復(fù)填子區(qū)煤層頂板等效標(biāo)高”的直接影響，近自然設(shè)計(jì)結(jié)果較原始地貌發(fā)生不規(guī)則形變。例如，在東西方向上，復(fù)填子區(qū)af1和af14近自然設(shè)計(jì)DEM較原始地貌有所抬升，而在其他區(qū)域較原始地貌有所下降。另一方面，受近自然地貌重塑模型中地表調(diào)整曲面模型影響（2.1.3節(jié)，地表調(diào)整曲面突出部分集中在區(qū)域中部），使得近自然設(shè)計(jì)DEM較原始地貌間形變也集中在區(qū)域中部。

與前者相似，相較于區(qū)域自然原始DEM，傳統(tǒng)設(shè)計(jì)結(jié)果也具有中部形變明顯的空間立體特征（圖6c、圖7c），然而與前者不同，受斜坡設(shè)計(jì)影響，內(nèi)排土場傳統(tǒng)設(shè)計(jì)DEM與周邊原始地貌間地貌“分割”明顯。例如，在同一位置由北至南分別做原始DEM、近自然設(shè)計(jì)DEM及傳統(tǒng)設(shè)計(jì)DEM三者的豎直剖面，在相同可用土方量控制下，相較于前兩者地貌“波動(dòng)式”形態(tài)（圖7a、圖7b），傳統(tǒng)設(shè)計(jì)DEM具有明顯的“斜坡-平臺(tái)-斜坡”三段式的地貌形態(tài)特征（圖7c）。

在土方運(yùn)移上，依據(jù)采復(fù)子區(qū)土方對應(yīng)關(guān)系，運(yùn)用MATLAB計(jì)算獲取近自然設(shè)計(jì)及傳統(tǒng)設(shè)計(jì)下土方平均運(yùn)距。其中，近自然設(shè)計(jì)DEM與傳統(tǒng)設(shè)計(jì)DEM兩者土方平均運(yùn)距分別為822.20 m/m3與821.71 m/m3，前者較后者僅高出約0.06%，兩者土方平均運(yùn)距極為接近（相差0.26 m/m3）。

3.3 設(shè)計(jì)結(jié)果抗水蝕能力評價(jià)

地貌穩(wěn)定性主要受氣候、地質(zhì)、植被等諸要素影響。其組成物料強(qiáng)度直接影響邊坡的穩(wěn)定性，是其抗侵蝕能力的主要指標(biāo)。依據(jù)地質(zhì)調(diào)查報(bào)告，確定研究區(qū)內(nèi)排土場地層結(jié)構(gòu)及巖土物理力學(xué)參數(shù)，并結(jié)合中國天氣網(wǎng)獲取研究區(qū)2007—2017年氣候參數(shù)。通過CLiDE演化模型，分別模擬3.2節(jié)所示10 a間內(nèi)排土場自然原始DEM、傳統(tǒng)設(shè)計(jì)DEM及近自然設(shè)計(jì)DEM的演變過程。其中，為保證模擬結(jié)果具有可比性，將研究區(qū)內(nèi)排土場傳統(tǒng)設(shè)計(jì)DEM的采復(fù)子區(qū)空間位置與3.2節(jié)近自然設(shè)計(jì)DEM設(shè)定為一致，其最終設(shè)計(jì)結(jié)果如圖7c所示。同時(shí)，為避免模型邊緣效應(yīng)影響，將研究區(qū)周邊1 000 m范圍共同納入運(yùn)算（圖7）。最后通過MATLAB疊加及對比分析，統(tǒng)計(jì)研究區(qū)內(nèi)排土場自然原始DEM、傳統(tǒng)設(shè)計(jì)DEM及近自然設(shè)計(jì)DEM下10 a的土壤水蝕總量（表1）。

如表1所示，整體上，本研究內(nèi)排土場10 a土壤水蝕量在4.098 3×105～23.681 5×105m3。其中，自然原始DEM土壤抗水蝕能力最強(qiáng)，模擬情境下其土壤10 a水蝕總量僅占樣本平均水平（14.069 5×105m3）的29.13%。在人工地貌設(shè)計(jì)下，本研究近自然設(shè)計(jì)DEM相較于傳統(tǒng)設(shè)計(jì)DEM可減少約39.07%的土壤水蝕量，其土壤水蝕總量逼近樣本平均水平。

為進(jìn)一步分析排土場各種地貌設(shè)計(jì)方式下對區(qū)域整體土壤水蝕的影響，將各DEM演化模擬結(jié)果與其初始DEM求差，獲取研究區(qū)內(nèi)排土場表土高度的空間變化影像（圖8），其中正數(shù)為該區(qū)域表土高度增高高度，負(fù)數(shù)為降低高度。

圖7 研究區(qū)內(nèi)排土場自然原始、傳統(tǒng)設(shè)計(jì)及近自然設(shè)計(jì)下DEM影像及橫向剖面示意Fig.7 DEM image and cross-sectional schematic of the natural primitive, traditional design and near-natural design of the internal dump

表1 研究區(qū)內(nèi)排土場自然原始、傳統(tǒng)設(shè)計(jì)及近自然設(shè)計(jì)下區(qū)域土壤10 a總水侵蝕量

圖8 研究區(qū)內(nèi)排土場演化前后自然原始、傳統(tǒng)設(shè)計(jì)及近自然設(shè)計(jì)下DEM表層高度變化Fig.8 The surface height changes of the natural original DEM, the traditional design DEM, and the near-natural design DEM before and after the evolution of the dump

如圖8所示，整體上，各樣本土壤在10 a的演變模擬中均發(fā)生了不同程度的運(yùn)移，使區(qū)域原有DEM發(fā)生不同程度的改變。在表土高度變化上，各DEM樣本高程變化區(qū)間均在-0.62～2.79 m之間，其變化幅度由大到小依次為：傳統(tǒng)設(shè)計(jì)DEM（3.39 m）、近自然設(shè)計(jì)DEM（3.35 m）和自然原始DEM（3.06 m）。排序結(jié)果與各樣本土壤10 a水蝕總量排序一致。在空間分布上，各樣本的土壤堆積及流失區(qū)域均有所不同。其中，傳統(tǒng)設(shè)計(jì)DEM表土流失及堆積區(qū)與研究區(qū)邊界高度相關(guān)，且主要沿邊界（圖8研究區(qū)邊界）呈向內(nèi)規(guī)則分布，視覺上與周邊自然地貌間具有明顯割裂，景觀及水文融合性較低。與前者不同，近自然設(shè)計(jì)DEM表土堆積及流失區(qū)空間分布與自然原始DEM較為相似，兩者均呈現(xiàn)出高度復(fù)雜的“離散化”空間分布。在周邊自然地貌融合上，兩者與周邊自然地貌間無明顯“痕跡”間隔，表土流失及堆積區(qū)域空間分布與研究區(qū)邊界無明顯相關(guān)，景觀及水文銜接性較好。并且，以自然原始DEM表土流失區(qū)（演化后高程下降區(qū)域）和堆積區(qū)（演化后高程上升區(qū)域）空間位置為參照的空間疊加分析表明，相較傳統(tǒng)設(shè)計(jì)DEM（空間重疊率：54.31%），近自然設(shè)計(jì)DEM土壤流失及堆積區(qū)具有更高的空間重疊率（57.27%）。

4 討 論

內(nèi)排土場作為露天礦山的重要組成部分，其占地面積隨工作面的推進(jìn)不斷增大。然而，受傳統(tǒng)復(fù)填方式影響，內(nèi)排土場區(qū)域原始微起伏地貌往往被替換為斜坡加平臺(tái)的“梯田式”人工堆墊地貌，使原有地表要素及水文特征喪失，并在集中降雨條件下，易引發(fā)表土水蝕等問題[2-4]。為此，專家學(xué)者提出依據(jù)臨近自然穩(wěn)定地貌特征以實(shí)現(xiàn)排土場地貌重塑[3,8,10]。例如，有學(xué)者以勝利一號礦排土場邊坡為例，提出以仿自然“反S形”邊坡替代傳統(tǒng)斜坡，減少區(qū)域表土侵蝕[10]；同樣，有學(xué)者以勝利一號礦為例，提出依據(jù)內(nèi)排土場周邊自然地貌參數(shù)，基于區(qū)域整體土方平衡，構(gòu)建內(nèi)排土場近自然地貌結(jié)果，實(shí)現(xiàn)景觀內(nèi)外融合，且相較傳統(tǒng)設(shè)計(jì)地貌具備更好的表土穩(wěn)定性[3]。與上述研究相似，本研究同樣以自然地貌為學(xué)習(xí)對象，且內(nèi)排土場近自然設(shè)計(jì)地貌較傳統(tǒng)設(shè)計(jì)地貌具有更強(qiáng)的抗水蝕能力（10年土壤水蝕減少約39.07%）。然而不同的是，上述研究強(qiáng)調(diào)以區(qū)域地貌統(tǒng)計(jì)特征為依據(jù)，本研究則強(qiáng)調(diào)以原始自然穩(wěn)定地貌為依據(jù)，即在排土場塑形過程中盡可能保留原始自然穩(wěn)定地貌特征，以實(shí)現(xiàn)礦區(qū)內(nèi)外地貌相連。

在露天礦采復(fù)過程中，受土體膨脹系數(shù)、可采煤層標(biāo)高、采復(fù)周期及采前地表標(biāo)高等諸要素影響，使得各采復(fù)周期間后方復(fù)填可用土方量存在差異；且由于復(fù)填區(qū)形態(tài)各異，造成內(nèi)排土場復(fù)填地貌存在走高、走低等諸多情況。在近自然設(shè)計(jì)中，已有坡面近自然重塑方法雖然能夠較好地應(yīng)用于傳統(tǒng)內(nèi)排土場復(fù)填區(qū)地貌設(shè)計(jì)，然而其實(shí)際形態(tài)仍以去除邊角的“梯田式”地貌為主[8,10]。后者提出的內(nèi)排土場近自然塑形方法雖具有良好的景觀融合效果，并考慮到開采前后土方平衡，然而設(shè)計(jì)中未考慮采復(fù)周期影響[3]，易造成采復(fù)周期中復(fù)填區(qū)實(shí)際可用土方量無法或超出原有設(shè)計(jì)表面等問題。受采復(fù)周期影響，本研究近自然設(shè)計(jì)地貌也呈現(xiàn)出部分區(qū)域較自然原始表面抬高、沉降等高程形變特征。然而與已有方法不同的是，本研究內(nèi)排土場近自然地貌重塑方法通過采復(fù)子區(qū)一一對應(yīng)的方式解決采復(fù)過程中土方動(dòng)態(tài)變化的問題，并相較于傳統(tǒng)設(shè)計(jì)地貌，近自然設(shè)計(jì)地貌與周邊自然地貌無明顯分界，模擬構(gòu)建上可實(shí)現(xiàn)采后地形與周邊自然地形的融合。

值得注意的是，在自然地貌中，區(qū)域原始地貌并不總是穩(wěn)定的[19]。在一定條件下，傳統(tǒng)復(fù)墾后礦區(qū)或?qū)⒏欣趨^(qū)域生態(tài)恢復(fù)[20]。由于本研究近自然地貌重塑方法以采前原始地貌為學(xué)習(xí)對象，故在構(gòu)建內(nèi)排土場地貌形態(tài)前，需充分論證原有自然地貌是否穩(wěn)定（例如，本研究自然原始地貌10年水蝕量最小，僅為傳統(tǒng)設(shè)計(jì)地貌的28.40%），以選擇有利于區(qū)域生態(tài)恢復(fù)的內(nèi)排土場設(shè)計(jì)方法。同時(shí)，受地下煤層剝離影響，本研究設(shè)計(jì)結(jié)果無法在形態(tài)上與原始自然地貌保持一致，并受采后松散土體空間不規(guī)則壓實(shí)沉降[21]，設(shè)計(jì)結(jié)果可能遠(yuǎn)離預(yù)期。故在后續(xù)研究中，可進(jìn)一步從力學(xué)視角，優(yōu)化本研究基于調(diào)整曲面的內(nèi)排土場近自然地貌重塑模型。

5 結(jié) 論

1）研究區(qū)內(nèi)排土場近自然設(shè)計(jì)結(jié)果表明，提出的近自然地貌重塑方法可在土方動(dòng)態(tài)平衡的前提下，依據(jù)采復(fù)子區(qū)采前地表、礦層頂板及底板數(shù)字高程數(shù)據(jù)，獲取可用地表標(biāo)高設(shè)計(jì)結(jié)果；在土方調(diào)配上，近自然設(shè)計(jì)地貌土方平均運(yùn)距為822.20 m/m3，較傳統(tǒng)設(shè)計(jì)地貌高0.06%，兩者土方平均運(yùn)距極為接近；且在視覺效果上，與周邊自然景觀“無痕”拼接，對區(qū)域開采前后地貌改變較少。

2）研究區(qū)10 a土壤水蝕模擬結(jié)果表明，本研究內(nèi)排土場近自然整形結(jié)果較傳統(tǒng)“梯田式”減少約39.07%的土壤流失量，表土抗水蝕能力較優(yōu)；在土壤演化運(yùn)移中，本研究近自然設(shè)計(jì)地貌與區(qū)域原始地貌在土壤堆積及流失上空間重疊率達(dá)57.27%，較傳統(tǒng)設(shè)計(jì)結(jié)果高2.96%，且視覺上呈“離散化”分布，與周邊自然地貌銜接較好。