<ol id="ebnk9"></ol> 
夏嘉南，李根生，卞正富，雷少剛，3，宮傳剛，李 恒
（1.中國礦業大學 礦山生態修復教育部工程研究中心，江蘇 徐州 221116；
2.中國礦業大學 公共管理學院，江蘇 徐州 221116；
3. 山東省煤田地質局第二勘探隊 山東省采煤塌陷地和采空區治理工程研究中心，山東 濟寧 272000）
內排土場作為露天礦區重要的組成部分，其占地面積隨采復工作面推進不斷增大[1]。在內排過程中，傳統地貌重塑方法往往將原始微起伏地貌景觀轉變為平臺和邊坡相間的“梯田式”人工堆疊景觀，使原始地表要素及水系特征消失[2]。雖然在后續復墾設計中，內排土場內布設輸水管網以滿足復墾植被生長需求，但整體上未能實現其與周邊地貌融合，易引發土地退化等問題[3-4]。因此，需嘗試在采復過程中保留原始自然地貌特征以破解上述困境。
近年來，不少專家學者就內排土場規劃及內排方案開展了研究。例如，王東等[5]以內排土場空間利用最大化為目標，基于采場工作線及邊坡角、內排土場邊坡角、基底傾角對排土線的影響規律，構建了采場工作幫、內排土場、端幫和基底的空間幾何模型；
栗嘉彬[6]以技術最佳適用性為目標，根據物料用途，構建了露天礦表土最優采集周期、最佳覆土厚度等求解模型；
OGGERI等[7]以最小復填成本為目標，根據場地、運輸等諸因素，提出了多情境下露天礦采復處置方法。然而上述研究中，未充分考慮重塑地貌結果與周邊自然地貌起伏貼合情況，導致內排土場往往被重塑為平臺加斜坡的“簡單梯田式”形態，人工地貌形態明顯，易引發嚴重的景觀破碎化問題[8-10]。相較于人工規則地貌，自然地貌作為長期演變的結果，具有高度的區域適宜性與穩定性，可作為修復區地貌重塑的參考對象[3,10-11]。且國外實踐證明[12-13]，基于自然的地貌重塑能更好地適應當地水氣條件，在有效提高區域抗水蝕能力的同時僅需要較少的維護費用[14]。因此，迫切需要一種貼合周邊自然地貌的內排土場重塑方法，用于解決內排土場設計結果與周邊自然地貌間貼合的問題[15]。
基于此，筆者以新疆黑山露天礦為例，提出一種采排復一體下基于調整曲面的內排土場近自然地貌重塑模型，并結合地理空間分析技術與MATLAB軟件，構建內排土場近自然地貌重塑結果。并通過CLiDE景觀演化模型，模擬了研究區10 a間原始地貌、傳統設計地貌及近自然設計地貌三者的土壤水蝕情況，對比評價了近自然設計地貌抗水蝕效果，為相同及相似地區內排土場地貌形態設計提供借鑒。
研究區位于新疆托克遜縣西北約90 km處，地理坐標：87°26′04″～87°29′40″E，43°13′48″～43°14′01″N，周邊地貌以山間谷地為主；
年均降水量約160.3 mm，年均氣溫約1.6 ℃，冬寒夏涼，屬大陸干旱及高寒氣候。區內表土以第四系為主，野生植被類型為典型草原植被。
依據地質勘探資料，選取吉布拉克煤礦探礦權范圍以南為研究區（圖1）。區域東西長約4 874 m、南北長約1 000 m。其中，預設采坑邊界東西長約812 m、南北長約1 000 m，采礦作業由東至西；
區域可采煤層分別為6、7、8、9、11、12-1、12-2、13-1和13-2號煤層，其中13-2號煤層為全區可露采、其余為局部可露采。
圖1 研究區示意Fig.1 Schematic of the study area
研究區采前自然原始數字高程數據（DEM）來源于地理空間數據云，拍攝時間2011年12月，條帶號87，行編號43，空間分辨率30 m×30 m；
研究區煤層底板、頂板標高及水平坐標數據來源于《新疆托克遜縣黑山礦區黑山露天煤礦資源儲量核實報告》。此外，基于上述報告通過ArcGIS克里金插值模塊獲取煤層底板及頂板空間分布數據；
內排土場采、復子區空間分布數據基于內排土場采復子區位置識別方法獲??；
復填子區內的可用土方量數據通過可用土方量計算式求解獲??；
復填區地表調整曲面數據通過地表調整曲面預構建模型獲??；
復填子區地表近自然DEM通過MATLAB軟件基于動態土方量一致、且區域斜率較緩的原則定向篩選不同B樣條控制點位置下的地表調整曲面獲得；
研究區自然原始DEM數據通過ArcGIS空間疊加提取獲得；
研究區內排土場傳統設計DEM數據基于土方量動態平衡下“30°斜坡加平臺”的設計思路獲??；
最后運用CLiDE演化模型與MATLAB軟件，分別模擬并統計研究區原始DEM、傳統設計DEM及近自然設計DEM三者10 a的土壤水蝕量，以評估近自然整形結果的抗侵蝕能力。
2.1.1 采復子區空間位置識別
依據露天礦采復周期，獲取圖2所示各周期下開采子區amn和復填子區afn水平投影空間位置。
圖2 露天礦采復周期下開采及復填子區空間位置提取示意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所示，依據露天礦實際采復周期，沿開采方向將露天礦區按采復周期時間先后順序劃分為n個子區。其中，開采子區和復填子區分別對應區域amn和區域afn。在露天礦內排構建的采復過程中，采坑每前進一個開采子區amn對應修復后側一個復填子區afn，直至內排過程結束，礦坑由區域Ab移至區域Aa。
2.1.2 可用土方量計算
為便于后續土方量計算，將開采及復填子區均細分為rowc行colc列的網格單元組（圖2）。在采復周期n中，開采子區amn可用于對應復填子區afn的土方量，計算式如下：
（1）
式中:V（amn）為開采子區amn可用于對應復填子區afn的土方量，其值是開采子區afn最下層可采礦層底板以上的總物質體積減去區域內所有可采礦層總物質體積；
i,j分別為開采及復填子區內網格單元的行、列序號，取取值分別在1～rowc和1～colc之間；
H（amnij）為開采子區amn網格單元（i,j）所在區域的采前地面標高；
hd（amnij）為開采子區amn網格單元（i,j）所在區域的最底層礦層底板標高；
ht（amnij）為開采子區amn網格單元（i,j）所在區域的可采礦層總厚度；
S（amnij）為開采子區amn子單元（i,j）的水平投影面積，其余參數解釋同上。
2.1.3 復填區地表調整曲面預構建
依據區域原有自然穩定地貌，構建圖3所示調整曲面模型，以盡可能保留復填子區afn原有地貌起伏特征。
圖3 地表調整曲面預構建示意Fig.3 Schematic of pre-built adjustment surface
如圖3所示，將反映復填子區最終地貌設計形態較其原始自然地貌空間高程變化的曲面定義為地表調整曲面。為保證地表調整曲面能使復填子區表面與周邊自然地貌“無縫”融合，選取邊界處高程變化小，中間變化大，且與區域自然地貌相似的連續放緩曲面。參照“船體”構建思路，采用連續性較優的B樣條曲線[16]為“龍骨”，三角函數為側方“支架”，構建復填子區afn所對應的地表調整曲面。圖中虛線箭頭為開采方向；
x,y,z三軸正方向分別對應開采方向、水平垂直開采方向及豎直方向；
fb和ft分別對應B樣條曲線函數和三角函數。其中，B樣條曲線fb位于復填子區afn水平投影中心線位置，采用控制點構建，其計算式如下：
fb（j）=Bezier（v,tb,ts,te）
（2）
（3）
其中:fb（j）為復填子區afn第j列的B樣條曲線值，定義域在0～1；
Bezier為豎直平面上的2維B樣條曲線函數，形態由起始點0,y1、終止點1,yNC和中間NC-2個控制點共同控制；
NC為Bezier函數中起點和終點在內所有控制點的個數，取值為>2的整數；
v為NC行2列的矩陣，由上至下行數據分別對應起點、控制點1至控制點NC-2、終點的坐標；
tb為起點所對應網格單元的列序號；
ts為Bezier函數取值間隔，是終點與起點所對應網格單元的列序號之差，取值為大于0的正數；
te是終點所對應網格單元的列序號；
x1～xNC-2分別為控制點1至控制點NC-2在中心線上以起點為原點的水平投影距離，其中x1≤x2≤…≤xNC-2；
y1～yNC分別是起點、控制點1至控制點NC-2、終點的預設高程值，定義域皆在0～1之間，其余參數解釋同上。在確定復填子區B樣條曲線函數形態的基礎上，通過三角函數確定地表調整曲面，其計算式如下：
（4）
式中:ft（i,j）為復填子區afn地表調整曲面網格單元（i,j）所對應的預設高程值，取值在0～1；
π為圓周率，其余參數解釋同上。
2.1.4 曲面土方量控制及坡度緩和優化
為保證采復周期中可用土方量動態平衡，通過豎直伸縮變化使地表調整曲面符合土方平衡要求，其計算式如下：
（5）
（6）
（7）
在此基礎上，為篩選出區域整體坡度較緩的復填子區地表高程模型，以提高其表土抗水蝕能力，構建以網格單元與周邊網格單元間坡角為對象的區域坡度評分標準，其計算式如下：
（8）
其中：fp（a）為坡度夾角a所對應的坡度評分值，其值為大于0的整數；
au為正整數，其取值在2°～90°，其余參數解釋同上。在區域坡度評分標準構建的基礎上，構建區域坡度評分模型，其計算式如下：
（9）
式中：p（afn）為復填子區afn的區域坡度總評分；
fpij（a）為子區afn內網格單元（i,j）沿某一朝向的坡度評分值；
∑fpij（a）為子單元（i,j）與周邊8個網格單元間坡角所對應評分值的總和。
依據上述評分模型，通過MATLAB軟件，調整2.1.3節中B樣條曲線起點、控制點1至控制點NC-2、終點的坐標位置，遍歷構建一組形態各異的復填子區地表高程模型。并以區域坡度總評分最小為目標，定向篩選復填子區Ht′求解結果，最終獲取原有自然形態表面下，土方量動態平衡且坡度平緩的復填子區afn地表高程近自然設計形態。
CLiDE模型是基于CAESAR模型開發的景觀地貌演化模型[17]。與傳統的侵蝕預測模型相比，其不僅能獲取侵蝕和沉積的時序數據，還能獲取最終發育的地貌形態。同時嵌入了Lisflood-FP流體力學相關模塊與地下水分布式水文模塊，可以模擬地表徑流-地貌-地下水三者間的交互作用和補給關系[18]，故極其適用于地貌修復效果評價[10]。為評價內排土場設計結果的抗水蝕能力，將研究區內排土場近自然設計結果為試驗組，原始自然地貌及傳統復填方法下內排土場設計結果為對照組，結合實地調查情況構建CLiDE模型，模擬其10 a的土壤侵蝕過程。并基于MATLAB軟件統計各組地貌演變前后土方變化量，以評價其抗水蝕效果。
研究區可采煤層共計9層，依據儲量報告，采用ArcGIS克里金插值模塊獲取研究區13-2號煤層底板標高平面分布圖及其各煤層厚度平面分布圖（圖4）。為避免數據空間分辨率差異造成的干擾，將影像空間分辨率統一設置為30 m×30 m，數據形式均為162列32行的柵格影像。
圖4 13-2號煤層底板標高平面分布與研究區各煤層厚度平面分布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
依據2.1.1節采復子區空間位置識別方法，在30 m×30 m柵格影像數據下，按采復周期將研究區劃分為45個一一對應的采復子區。其中，采復子區網格行數為32，列數為3，即復填子區和開采子區均設計為水平投影東西長90 m，南北長960 m的水平矩形區域，以使采復工作前后采坑的水平投影大小不變（圖5）。
圖5 研究區開采子區及復填子區空間位置示意Fig.5 Schematic of the spatial location of the mining sub-area and the refill sub-area
在此基礎上，依據2.1節所示內排土場地表近自然設計方法，通過MATLAB軟件以0.01為控制點數值最小變化間距，膨脹系數k取均值1.1，求解研究區內排土場地貌近自然設計DEM結果（圖6）。其中，考慮到礦區可采煤層非惟一（9層），且厚度及底板標高存在差異，依據研究區三維地質資料，以13-2號煤層的底板標高為研究區煤層底板標高，并采用空間疊加方法疊加各煤層厚度數據，獲取研究區煤層頂板標高等效數據，以實現本研究內排土場近自然地貌重塑模型于多礦層復雜情景下的應用。
圖6 研究區內排土場DEM圖Fig.6 DEM of the internal dump in the study area
如圖6所示，視覺效果上，研究區內排土場近自然設計DEM與原始自然地貌及紋理特征較為相似，且高程值均在+2 000.26～+2 959.00 m；
表土高程上，受土方平衡約束，近自然設計DEM較自然原始DEM存在差異，變化幅度在-205.90～50.03 m。一方面，受研究區“開采子區地表原始標高”、“復填子區煤層底板標高”和“復填子區煤層頂板等效標高”的直接影響，近自然設計結果較原始地貌發生不規則形變。例如，在東西方向上，復填子區af1和af14近自然設計DEM較原始地貌有所抬升，而在其他區域較原始地貌有所下降。另一方面，受近自然地貌重塑模型中地表調整曲面模型影響（2.1.3節，地表調整曲面突出部分集中在區域中部），使得近自然設計DEM較原始地貌間形變也集中在區域中部。
與前者相似，相較于區域自然原始DEM，傳統設計結果也具有中部形變明顯的空間立體特征（圖6c、圖7c），然而與前者不同，受斜坡設計影響，內排土場傳統設計DEM與周邊原始地貌間地貌“分割”明顯。例如，在同一位置由北至南分別做原始DEM、近自然設計DEM及傳統設計DEM三者的豎直剖面，在相同可用土方量控制下，相較于前兩者地貌“波動式”形態（圖7a、圖7b），傳統設計DEM具有明顯的“斜坡-平臺-斜坡”三段式的地貌形態特征（圖7c）。
在土方運移上，依據采復子區土方對應關系，運用MATLAB計算獲取近自然設計及傳統設計下土方平均運距。其中，近自然設計DEM與傳統設計DEM兩者土方平均運距分別為822.20 m/m3與821.71 m/m3，前者較后者僅高出約0.06%，兩者土方平均運距極為接近（相差0.26 m/m3）。
地貌穩定性主要受氣候、地質、植被等諸要素影響。其組成物料強度直接影響邊坡的穩定性，是其抗侵蝕能力的主要指標。依據地質調查報告，確定研究區內排土場地層結構及巖土物理力學參數，并結合中國天氣網獲取研究區2007—2017年氣候參數。通過CLiDE演化模型，分別模擬3.2節所示10 a間內排土場自然原始DEM、傳統設計DEM及近自然設計DEM的演變過程。其中，為保證模擬結果具有可比性，將研究區內排土場傳統設計DEM的采復子區空間位置與3.2節近自然設計DEM設定為一致，其最終設計結果如圖7c所示。同時，為避免模型邊緣效應影響，將研究區周邊1 000 m范圍共同納入運算（圖7）。最后通過MATLAB疊加及對比分析，統計研究區內排土場自然原始DEM、傳統設計DEM及近自然設計DEM下10 a的土壤水蝕總量（表1）。
如表1所示，整體上，本研究內排土場10 a土壤水蝕量在4.098 3×105～23.681 5×105m3。其中，自然原始DEM土壤抗水蝕能力最強，模擬情境下其土壤10 a水蝕總量僅占樣本平均水平（14.069 5×105m3）的29.13%。在人工地貌設計下，本研究近自然設計DEM相較于傳統設計DEM可減少約39.07%的土壤水蝕量，其土壤水蝕總量逼近樣本平均水平。
為進一步分析排土場各種地貌設計方式下對區域整體土壤水蝕的影響，將各DEM演化模擬結果與其初始DEM求差，獲取研究區內排土場表土高度的空間變化影像（圖8），其中正數為該區域表土高度增高高度，負數為降低高度。
圖7 研究區內排土場自然原始、傳統設計及近自然設計下DEM影像及橫向剖面示意Fig.7 DEM image and cross-sectional schematic of the natural primitive, traditional design and near-natural design of the internal dump
表1 研究區內排土場自然原始、傳統設計及近自然設計下區域土壤10 a總水侵蝕量
圖8 研究區內排土場演化前后自然原始、傳統設計及近自然設計下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的演變模擬中均發生了不同程度的運移，使區域原有DEM發生不同程度的改變。在表土高度變化上，各DEM樣本高程變化區間均在-0.62～2.79 m之間，其變化幅度由大到小依次為：傳統設計DEM（3.39 m）、近自然設計DEM（3.35 m）和自然原始DEM（3.06 m）。排序結果與各樣本土壤10 a水蝕總量排序一致。在空間分布上，各樣本的土壤堆積及流失區域均有所不同。其中，傳統設計DEM表土流失及堆積區與研究區邊界高度相關，且主要沿邊界（圖8研究區邊界）呈向內規則分布，視覺上與周邊自然地貌間具有明顯割裂，景觀及水文融合性較低。與前者不同，近自然設計DEM表土堆積及流失區空間分布與自然原始DEM較為相似，兩者均呈現出高度復雜的“離散化”空間分布。在周邊自然地貌融合上，兩者與周邊自然地貌間無明顯“痕跡”間隔，表土流失及堆積區域空間分布與研究區邊界無明顯相關，景觀及水文銜接性較好。并且，以自然原始DEM表土流失區（演化后高程下降區域）和堆積區（演化后高程上升區域）空間位置為參照的空間疊加分析表明，相較傳統設計DEM（空間重疊率：54.31%），近自然設計DEM土壤流失及堆積區具有更高的空間重疊率（57.27%）。
內排土場作為露天礦山的重要組成部分，其占地面積隨工作面的推進不斷增大。然而，受傳統復填方式影響，內排土場區域原始微起伏地貌往往被替換為斜坡加平臺的“梯田式”人工堆墊地貌，使原有地表要素及水文特征喪失，并在集中降雨條件下，易引發表土水蝕等問題[2-4]。為此，專家學者提出依據臨近自然穩定地貌特征以實現排土場地貌重塑[3,8,10]。例如，有學者以勝利一號礦排土場邊坡為例，提出以仿自然“反S形”邊坡替代傳統斜坡，減少區域表土侵蝕[10]；
同樣，有學者以勝利一號礦為例，提出依據內排土場周邊自然地貌參數，基于區域整體土方平衡，構建內排土場近自然地貌結果，實現景觀內外融合，且相較傳統設計地貌具備更好的表土穩定性[3]。與上述研究相似，本研究同樣以自然地貌為學習對象，且內排土場近自然設計地貌較傳統設計地貌具有更強的抗水蝕能力（10年土壤水蝕減少約39.07%）。然而不同的是，上述研究強調以區域地貌統計特征為依據，本研究則強調以原始自然穩定地貌為依據，即在排土場塑形過程中盡可能保留原始自然穩定地貌特征，以實現礦區內外地貌相連。
在露天礦采復過程中，受土體膨脹系數、可采煤層標高、采復周期及采前地表標高等諸要素影響，使得各采復周期間后方復填可用土方量存在差異；
且由于復填區形態各異，造成內排土場復填地貌存在走高、走低等諸多情況。在近自然設計中，已有坡面近自然重塑方法雖然能夠較好地應用于傳統內排土場復填區地貌設計，然而其實際形態仍以去除邊角的“梯田式”地貌為主[8,10]。后者提出的內排土場近自然塑形方法雖具有良好的景觀融合效果，并考慮到開采前后土方平衡，然而設計中未考慮采復周期影響[3]，易造成采復周期中復填區實際可用土方量無法或超出原有設計表面等問題。受采復周期影響，本研究近自然設計地貌也呈現出部分區域較自然原始表面抬高、沉降等高程形變特征。然而與已有方法不同的是，本研究內排土場近自然地貌重塑方法通過采復子區一一對應的方式解決采復過程中土方動態變化的問題，并相較于傳統設計地貌，近自然設計地貌與周邊自然地貌無明顯分界，模擬構建上可實現采后地形與周邊自然地形的融合。
值得注意的是，在自然地貌中，區域原始地貌并不總是穩定的[19]。在一定條件下，傳統復墾后礦區或將更有利于區域生態恢復[20]。由于本研究近自然地貌重塑方法以采前原始地貌為學習對象，故在構建內排土場地貌形態前，需充分論證原有自然地貌是否穩定（例如，本研究自然原始地貌10年水蝕量最小，僅為傳統設計地貌的28.40%），以選擇有利于區域生態恢復的內排土場設計方法。同時，受地下煤層剝離影響，本研究設計結果無法在形態上與原始自然地貌保持一致，并受采后松散土體空間不規則壓實沉降[21]，設計結果可能遠離預期。故在后續研究中，可進一步從力學視角，優化本研究基于調整曲面的內排土場近自然地貌重塑模型。
1）研究區內排土場近自然設計結果表明，提出的近自然地貌重塑方法可在土方動態平衡的前提下，依據采復子區采前地表、礦層頂板及底板數字高程數據，獲取可用地表標高設計結果；
在土方調配上，近自然設計地貌土方平均運距為822.20 m/m3，較傳統設計地貌高0.06%，兩者土方平均運距極為接近；
且在視覺效果上，與周邊自然景觀“無痕”拼接，對區域開采前后地貌改變較少。
2）研究區10 a土壤水蝕模擬結果表明，本研究內排土場近自然整形結果較傳統“梯田式”減少約39.07%的土壤流失量，表土抗水蝕能力較優；
在土壤演化運移中，本研究近自然設計地貌與區域原始地貌在土壤堆積及流失上空間重疊率達57.27%，較傳統設計結果高2.96%，且視覺上呈“離散化”分布，與周邊自然地貌銜接較好。
<ol id="ebnk9"></ol>