水文保留曲面下草原露天礦內排土場地貌近自然重塑

夏嘉南，李根生，卞正富，雷少剛，金慧霞,3

(1.中國礦業大學 公共管理學院，江蘇 徐州 221116；2.中國礦業大學 礦山生態修復教育部工程研究中心，江蘇 徐州 221116；3.溫州職業技術學院，浙江 溫州 325035)

露天開采作為我國草原礦區的重要開采方式，其內排土場占地面積隨開采工作面的不斷推進而增大。在此過程中，傳統地貌重塑方法往往將原始地貌轉變為“臺階式”或“錐狀堆積式”的人工堆疊地貌，使原有的地貌及水文特征消失，導致礦區重塑地貌與周邊自然地貌割裂。雖然在后續復墾設計中，內排土場通過輸水管網布設以滿足復墾區重建植被的生長需求，但未能實現礦區內外地貌水系的“無痕”拼接，從而在高強度降雨條件下，易引發邊界處嚴重表土水蝕等問題。相對人工規則地貌，自然穩定地貌作為長期演變的結果，具有較好的適宜性與系統性。因此，可嘗試在采排復過程中保留自然地貌及其水文特征以破解上述困境。

當前，GeoFluv模型作為主要的近自然地貌重塑方法，被廣泛應用于采損區地貌重塑，且塑形效果在美國墨西哥州和懷俄明州、西班牙、澳大利亞等礦區均得到了良好的檢驗，并在國內得到了初步的運用。然而，受模型框架影響，重塑區水文溝道位置需人為主觀設定，雖后續通過自然穩定溝道特征值加以修正，但整體上其地貌求解結果主觀性較高。NbS(Nature-based Solutions)理論謀求人與生態系統互惠，與“保護優先，自然恢復為主”的中國生態文明建設理念相契合。在上述研究及理論的指導下，團隊內諸多學者均提出師法自然穩定地貌的內排土場地貌重塑方法。其中，李恒等以內蒙古勝利礦為例，通過將周邊自然穩定地貌剖面形態應用至排土場邊坡，通過GeoWep模擬論證仿自然形態邊坡較傳統“臺階式”邊坡更具抗水蝕能力。陳航等以內蒙古勝利礦為例，通過提取周邊自然穩定地貌特征參數，在不考慮采復周期的情景下采用GeoFlue軟件構建研究區內排土場近自然重塑地貌，并基于CLiDE模擬驗證其表土抗水蝕能力優于傳統內排土場地貌重塑結果。然而，在上述研究中，內排土場采排復一體下復填區的可用土方量會隨采復周期變化。同時，相較于周邊自然穩定地貌統計結果，區域原始自然穩定地貌的特殊性(受區域氣候及相鄰區的長期作用影響)也未得到充分考慮。在前置研究中，筆者以新疆黑山露天礦為例，提出基于SAS(Surface Adjustment Surface)的內排土場地貌近自然重塑模型，強調通過局部調整的方式以盡可能保留原始地貌起伏特征，同時避免模型構建中水文溝道的人為主觀設計，且CLiDE水蝕模擬表明，NNDL_SAS(Near-natural Design Landform with Surface Adjustment Surface)相較于TDL(Traditional Design Landform)在10 a當地氣候條件下可減少約39.07%的表土水蝕量。然而，后續論證表明，上述模型求解形態局部調整可能會造成水系斷流、淤積等不可控問題。因此，需要一種既貼合周邊自然地貌，又保留原始地貌及水文特征的內排土場地貌重塑方法，以削弱前置研究中地貌設計結果水系斷流、淤積等問題。

為此，筆者在前置研究的基礎上，以新疆黑山露天礦為例，提出一種基于HPCS(Hydrological Preserved Curved Surface)的內排土場近自然地貌重塑模型，求解并以NNDL_HPCS(Near-natural Design Landform with Hydrological Preserved Curved Surface)為實驗組，NNDL_SAS,TDL和PNL(Pristine Natural Landform)為對照組，結合GIS空間分析軟件與CLiDE演化模型，模擬并討論NNDL_HPCS在地貌融合、土方調配和土壤抗水蝕3方面的改良效果。為相同及相似地區NbS理論下礦山近自然生態修復提供方法支持。

1 研究區概況與數據來源

1.1 研究區概況

研究區位于新疆黑山露天礦礦權范圍內，地理坐標：87.43°E—87.49°E，43.23°N—43.23°N，周邊地形以山間谷地為主，屬大陸半干旱及高寒氣候；區內表土以第四系為主，野生植被類型為典型草原植被。采區東西長約4 874 m、南北長約1 000 m；采坑預設邊界為812 m×1 000 m的矩形，內排作業由東至西(圖1)。

圖1 研究區采前原始DEM和真彩色衛星影像Fig.1 DEM of the natural surface of the study area before acquisition and RGB satellite image of the study area

依據研究區煤炭資源儲量核實報告，區域可采煤層為6,7，8,9，11,12-1,12-2,13-1和13-2號。其中13-2號煤層為全區可露采、其余為局部可露采(圖2)。內排土場臺階高度均設置為10 m。其中西幫為開采幫，幫坡角為20°；南北兩側邊幫幫坡角均為13°；東幫為內排幫，幫坡角為8°。

1.2 數據來源

研究區采前原始DEM(Digital Elevation Model)、煤層頂板DEM和煤層底板DEM均來源于前置研究，空間分辨率均為30 m×30 m；內排土場開采和復填子區空間分布數據與前置研究相同，均基于內排土場采復子區位置識別方法獲取；復填子區內的可用土方量數據通過可用土方量計算式求解獲??；復填區HPCS數據通過HPCS預構建模型獲??；復填子區NNDL_HPCS通過Matlab軟件基于土方平衡下區域坡度最緩原則定向篩選不同控制點系數下預構建的HPCS獲得；復填子區NNDL_SAS來源于前置研究成果中“近自然設計DEM”；研究區PNL數據來源于前置研究成果中“自然原始地貌”；研究區內排土TDL與前置研究成果中的“傳統設計DEM”一致(30°斜坡加平臺式“階梯狀”布局)；地貌融合效果中，各DEM相對于PNL的高程方差值通過ArcGIS空間疊加分析方法，逐柵格累計統計DEM與PNL間柵格高程差的平方值獲取；地貌融合效果中，自然水文保留率指標通過ArcGIS坡向模塊和空間疊加分析方法，提取并對比研究區NNDL_HPCS,NNDL_SAS和TDL三者圖像柵格坡向與PNL空間分布一致性獲取；土方調配中，NNDL_HPCS,NNDL_SAS和TDL的土方平均運距均由Matlab基于每立方土體空間直線位移距離統計獲取；最后運用CLiDE演化模型和Matlab軟件，分別模擬并統計研究區NNDL_HPCS,NNDL_SAS,TDL和PNL四者10 a的土壤水蝕量，以評估NNDL_HPCS的抗水蝕能力。

圖2 研究區各煤層厚度及最底層煤層底板標高空間分布Fig.2 Spatial distribution of the coal seam thickness and bottom plate elevation of the lowest coal seam in the study area

2 研究方法

2.1 基于HPCS的內排土場近自然地貌重塑模型

2.1.1 采復子區空間位置識別

依據露天礦采復周期，獲取開采子區m和復填子區f的水平投影空間位置(圖3)。

圖3 礦區采復周期n中采前和復填后豎直剖面示意((a)，(b))以及采復子區網格單元示意(c)Fig.3 Schematic diagram of the vertical section before and after mining in the mining area under the open-pit mining-refill cycle n， schematic diagram of the grid unit of the mining-refill sub-area

如圖3所示，依據露天礦實際采復周期及其時間先后順序，沿開采方向將露天礦區劃分為個子區。其中，開采子區和復填子區分別對應區域～m和區域～f。在露天礦內排構建的采復過程中，采坑每前進一個開采子區m對應修復后側一個復填子區f(圖3(a)，(b))，直至結束內排過程，礦坑由區域不斷迭代前進至采礦終止邊界。

2.1.2 可用土方量計算

基于輸入數據柵格空間分辨率，將開采及復填子區均細分為row行col列的網格單元組(圖3(c))。其中，網格單元大小與柵格分辨率一致。在采復周期中，開采子區m可用于復填子區f的土方量計算式為

(1)

式中，(m)為開采子區m可用于復填子區f的土方量，其值是開采子區m總巖土體積(最底層底板以上所有物質)與其可采礦層總物質體積之差；為土體修正系數，是煤層受剝離損失、夾矸、巖土膨脹等多因素綜合影響下實際土方量與剝離土方量之比；，分別為采復子區內網格單元的行、列序號，定義域分別為1～row與1～col；(m)為開采子區網格單元m(,)所在區域的采前地面平均標高；(m)為開采子區網格單元m(,)所在區域最底層礦層的底板平均標高；(m)為開采子區網格單元m(,)所在區域的可采礦層平均厚度；(m)為開采子區網格單元m(,)的水平投影面積，是行列分辨率的乘積，其余參數解釋同上。

2.1.3 復填區HPCS預構建

自然穩定地貌作為區域長時序演化的結果，其不僅形態穩定，且地表水文特征與上下流域間高度匹配。為構建與周邊自然地貌相銜接的重塑地貌，依據NbS理論，應對其水文特征進行保留，以實現重塑地貌與周邊自然地貌相融合。HPCS模型是在上述理論基礎上，通過對水文溝道和坡頂進行識別，令其中水文溝道和邊界處高程不變，以保證地貌重塑前后內外水系相通，匯水面關系不變；通過調整區域內各坡頂高程，以實現土方量平衡；并通過目標篩選，以獲取最優地貌設計形態的地表高程調整模型。在數值上，其是重塑地貌DEM較采前自然地貌DEM高程變化的空間分布數據。

在HPCS求解過程中，首先需構建一批形態各異的HPCS預構建模型，以為后續土方控制和坡度緩和優化提供可供操作的基底數據。故基于其模型定義與NbS理論，為盡可能保留復填子區f原始水文特征，通過ArcGIS水文分析模塊與Matlab軟件，創建圖4所示HPCS預構建模型。

圖4 HPCS預構建示意Fig.4 Schematic diagram of pre-construction of HPCS

如圖4所示，對原有自然穩定地貌進行概化，以減少NNDL_HPCS求解結果與PNL間的匯水面分布差異。其中，依據概化要素于原有地貌水文所起作用，劃分為分割匯水面的地形突出點(坡頂點)和匯水面內水流集中的溝道線(水系或者溝道)。其中，水文溝道網格提取結果通過ArcGIS水文分析模塊獲取，其最小匯水面積參數通過均值變點法(OTSU)確定。原始地貌的地形突出點則基于式(2)提?。?/p>

(2)

其中，(,)為復填子區網格單元(，)的取值，凸點取值為1，其余為0；area為以網格單元(，)為中點的*高程搜索窗口，area(2，2)為窗口中心網格單元(，)所對應的高程值，其余參數解釋同上。當搜索網格單元中心點單元area(2，2)高程為其窗口內最大值，定義該中心網格為凸點。通過Matlab遍歷填挖后復填子區PNL，識別凸點網格單元所在位置。在此基礎上，將復填子區凸點網格識別結果與水文溝道網格提取結果進行空間疊加，構建式(3)判別式對HPCS網格單元進行預賦值：

(3)

其中，為HPCS網格單元預賦值結果;(,)為其網格單元(,)的取值，凸點網格賦值為w，其余為0；w為凸點網格(，)所對應的無量綱系數，定義域為[0,1]，其余參數解釋同上。在此基礎上，基于式(4)所示模糊算法，預構建HPCS：

(4)

2.1.4 曲面土方量控制及坡度緩和優化

同SAS中的土方動態平衡方法，通過豎直伸縮變化使預構建的HPCS符合土方平衡要求，其計算式為

(5)

′(,)=(,)

(6)

′(f)=′(,)+(f)

(7)

式中，為伸縮變化系數，其值是復填子區f剝離的總巖土(包含礦層)體積減去由開采子區m獲取的膨脹后的可用土方量之后，與預構建的HPCS總體積比值的相反數；′(,)為預構建HPCS伸縮變換后網格單元(,)的值，是復填子區f最終設計結果(NNDL_HPCS)與PNL間高程值的變化；′(f)為NNDL_HPCS中復填子區網格單元f(，)所對應的地表標高，其值等于該處的自然原始標高與伸縮變換后的HPCS值之和。

在此基礎上，為篩選出區域表土抗水蝕能力較強(坡度平緩)且與周邊地貌融合的NNDL_HPCS，構建相鄰網格單元間坡角為對象的區域坡度評分標準。在坡度評分標準設置中，為避免區域不同坡度組合間的分值混淆，同時減少軟件運算中非必要的內存占用，以1°為步長，令相鄰步長間高坡度分值是低坡度分值乘上相鄰網格數(即8個)的值再加1。其計算式為

(8)

其中，()為坡度夾角(°)所對應的坡度評分值，其值域∈；為正整數，定義域在2°～90°。在區域坡度評分標準構建的基礎上，構建區域坡度評分模型，其計算式為

(9)

其中，(f)為復填子區f的區域坡度總評分值；p()為子區f內網格單元(,)沿某一朝向的坡度評分值;∑p為子單元(,)與周邊8個網格單元間坡角所對應評分值的總和。為確保最終篩選對象與周邊已知地貌無縫拼接，除復填子區f與礦坑鄰接處網格單元外，其余網格單元都需統計與周邊8個網格單元的坡度評分值。

將設為w的最小變化步長，通過上述HPCS預構建模型，由0至1遍歷獲取一批不同w組合下的預構建HPCS，進而獲取一組形態各異的復填子區′。在此基礎上，以區域坡度評分最小為目標，定向篩選復填子區′求解結果，最終獲取原有自然形態表面下，土方量動態平衡且坡度平緩的復填子區f內NNDL_HPCS。

2.2 水文保留率

坡向是地貌的重要表征參數，對地表徑流起著直接導向作用，同時受區域光照、溫度、雨量、風速等因子的綜合影響，間接決定區域植被類型及其生長狀況。基于灰色系統理論，以自然穩定地貌坡向空間特征為塑形導向。將水文保留率定義為重塑地貌與采前自然地貌坡向空間分布上一致性，即2者相同坡向區域的面積占總區域面積的比例。通過ArcGIS坡向組件，提取PNL,NNDL_HPCS,NNDL_SAS和TDL的坡向空間柵格數據。以PNL坡向空間柵格數據為基準，通過Matlab軟件逐柵格進行空間疊加分析，計算NNDL_HPCS,NNDL_SAS,TDL三者與PNL的區域坡向一致性，其計算式為

(10)

其中，為地貌相對于PNL的水文保留率，其定義域在[0,1]，若2者坡向空間一致性越高(即相同位置2者間的坡向夾角越小)，其取值越大；row，col分別為地貌坡向柵格數據的總行、總列數；為地貌的坡向空間柵格數據，其中(,)為其第行第列的柵格的坡向值，其取值大小符合羅馬盤規則；為PNL的坡向空間柵格數據，其中(,)是其第行第列柵格的坡向值，其取值大小也符合羅馬盤規則。在ArcGIS坡向模塊中，平臺坡向賦值為-1，考慮到平臺坡向的萬向性(可視作任意方向)，默認將平臺處的坡向認定為一致。

2.3 土方直線運距

基于露天采排復一體化理念，令開采工作面土方由上層臺階逐層向下剝離，后側復填工作面由煤層剝離后底板處逐層至上復填。其中，復填區當且僅當相鄰低處臺階的可用土方量(考慮膨脹系數)完全滿足復填所需時，才進行相鄰高處臺階的復填作業。由于實際采排復一體化工作中，土方運移受運輸通道等諸因素影響，其難以精確計算，故依據上述思路，按照研究區邊幫設計參數，以10 m為分段間隔，通過Matlab建立采復前后土方的空間對應索引，并依據其每m土體中心點的位移，計算其土方平均運距，以此簡要評價NNDL_HPCS,NNDL_SAS和TDL的土方調配特征。其計算式為

(11)

其中，為土方平均運距；為運移土方按時間的排序，取值為以1為步長從0至(m)的遞進值；(m,)，′(m,)分別為開采子區m中第運移的立方土體運移前后中心點所對應的三維坐標，其余參數解釋同上。通過統計各采復周期中各立方土體運移前后位移所對應的模長均值，簡要評價重塑地貌土方運移特征。

2.4 抗水蝕能力評價

CLiDE模型是基于CAESAR模型開發的景觀地貌演化模型。相比于傳統的侵蝕預測模型，在獲取侵蝕和沉積的時序數據的同時，還能獲取最終發育的地貌形態。并可模擬地表徑流-地貌-地下水3者間的交互作用和補給關系，極其適用于地貌修復效果評價。為評價內排土場設計結果的抗水蝕能力，將研究區NNDL_HPCS設為實驗組，NNDL_SAS,TDL和PNL設為對照組，套用前置研究設定參數(氣象數據來源于中國氣象網)構建CLiDE模型，模擬4者10 a的土壤侵蝕過程。其中，依據研究區典型草原植被類型，將土地利用類型統一設置為第5類(適用于半野生草地、粗糙的草地、中性草地等)，并基于區域植被以1 a生為主，將其植被恢復時間統一設置為1 a。并基于Matlab軟件統計各組地貌演變前后土方變化量，以評價其抗水蝕效果。其中為排除礦坑最終位置對土壤水蝕模擬結果的干擾，在PNL對應位置處設置相同大小的礦坑。其計算式為

第一，專業人才匱乏。當下養老服務存在的主要困境是養老專業人才隊伍的匱乏，導致服務專業水平不足，影響服務質量和服務成效。在養老服務提供方的人才隊伍中，主要包括志愿者、居委會工作人員，以及非專業出身的其他領域轉行而從事養老服務的人員等，在養老服務的能力和水平等方面都有待考量，這對于前期智慧養老服務模式的建設是一個亟待突破的限制性因素。更值得注意的是，智慧養老的運營模式不僅需要的是專業的養老服務人才，更是需要懂大數據技術、大數據工具的數據專業領域的養老服務復合型人才。

(12)

式中，為土壤水蝕總量；(f)，(f)分別為復填子區網格單元f(,)中地貌演化前后的地表標高均值；(f)為復填子區網格單元f(,)的水平投影面積，為行列分辨率的乘積。

3 結果與分析

3.1 數據準備結果及分析

研究區可采煤層共計9層，依據儲量報告，采用ArcGIS克里金插值模塊獲取研究區最底層煤層底板標高空間分布圖(圖2(j)，對應研究區13-2號煤層))及其各煤層厚度平面分布圖(圖2(a)～(i)，分別對應研究區6，7，8，9，11，12-1，12-2，13-1，13-2號煤層)。為避免數據空間分辨率差異造成的干擾，將影像空間分辨率統一設置為30 m×30 m，數據形式均為162列32行的柵格影像。

考慮到研究區可采煤層非惟一(9層)，且其厚度與底板標高存在差異，故依據三維地質資料，以13-2號煤層的底板標高為研究區煤層底板標高，并采用空間疊加方法疊加各煤層厚度數據，以獲取研究區煤層頂板標高等效數據，實現本研究地貌重塑模型于多礦層復雜情景下的應用。

3.2 近自然設計地貌構建結果及分析

依據2.1.1節采復子區空間位置識別方法，將研究區劃分為45個一一對應的采復子區，與已有研究保持一致。其中，采復子區網格行(row)列(col)數均為32行3列，其水平投影為東西長90 m，南北長960 m的矩形。通過2.1節內排土場近自然地貌重塑模型，將w的最小變化步長取0.01，土體修正系數取1.1，求解研究區NNDL_HPCS(圖5(b))。

圖5 研究區內排土場DEMFig.5 DEM of the internal dump in the study area

3.2.1 地貌融合結果與分析

在視覺效果上，相較于TDL(圖5(d))“階梯狀”地貌空間形態，NNDL_HPCS(圖5(b))和NNDL_SAS(圖5(c))兩者與PNL(圖5(a))較為相似，存在復雜的自然地貌紋理特征；在表土高程變化上，受土方平衡約束，NNDL_HPCS,NNDL_SAS,TDL三者較PNL高程變化幅度分別為-112～171,-205.90～50.03,-160.21～73.60 m，變化區間大小分別為NNDL_HPCS(283 m) > NNDL_SAS(250.93 m)>TDL(233.81 m)；在空間形變特征上，受“開采子區地表原始標高”、“復填子區煤層底板標高”、“復填子區煤層頂板等效標高”、“膨脹系數”4個要素的直接影響，本研究NNDL_HPCS與已有成果NNDL_SAS均較PNL發生空間上的不規則形變。例如，在東西方向上，復填子區～NNDL均較PNL有所抬升，而在其他區域均較原始地貌有所下降。與前兩者相似，TDL較PNL也具有明顯的空間形變特征(圖5(d)，6(d))，然而與前兩者不同，受斜坡設計影響，TDL與周邊原始地貌間“分割”明顯，具“斜坡-平臺-斜坡”3段式的地貌形態特征(圖5(d)，6(d))。

如圖6所示，以研究區西邊界為起點，以PNL為基準，由西向東逐列分別統計NNDL_HPCS,NNDL_SAS,TDL相對于PNL的高程方差。整體上，NNDL_HPCS,NNDL_SAS和TDL三者相對PNL的高程總方差分別為288 942.8,582 248.8和465 452.8 m。其中，NNDL_HPCS相對于PNL的高程方差最小，相較NNDL_SAS與TDL分別低50.37%與37.92%，在絕對高程上更接近于PNL，形態保留效果最優。隨距研究區西邊界距離變化，在0～810 m段，NNDL_HPCS,NNDL_SAS,TDL,PNL四者均位于礦坑位置，列相對高程方差值均為0；在810～1 800 m段和4 590～4 890 m段，NNDL_HPCS列相對高程方差遠高于NNDL_SAS和TDL，其均值分別為后2者均值的11.83倍與6.99倍，和17.16倍與6.43倍；在1 800～4 590 m段，NNDL_HPCS列相對高程方差值則遠低于NNDL_SAS和TDL，其均值分別為后兩者均值的0.13倍與0.18倍。由地質資料可知，在810～1 800 m段和4 590～4 890 m段，隨著前側開采子區煤層底板走低，可采煤層厚度變化較少，導致后側復填可用土方量不斷增加，NNDL_HPCS,NNDL_SAS和TDL對應復填子區均形成“凸包地形”；在1 800～4 590 m段，隨著前側開采子區煤層底板走低，可采煤層厚度增大，導致后側復填可用土方量穩中有減，NNDL_HPCS,NNDL_SAS和TDL對應復填子區均形成“凹陷地形”(圖6)。即在“凸包地形”段，NNDL_HPCS與PNL柵格絕對高程重合效果上遠弱于NNDL_SAS和TDL，在“凹陷地形”段則遠優于2者。

圖6 研究區PNL,NNDL_HPCS,NNDL_SAS,TDL三維示意(a)～(d)，及其復填子區物質等效厚度隨距研究區西邊界距離的 折線統計(e)和3者相對于PNL的列高程方差隨距研究區西邊界距離的折線統計(f)Fig.6 Schematic diagram of the PNL，NNDL_HPCS，NNDL_SAS,TDL(a)～(d)，the broken line statistics chart of the equivalent thickness of the material in the refill subregion with the distance from the western boundary(e) and they relative to the PNL column elevation variance with the distance from the western boundary in the study area(f)

圖7 研究區NNDL_HPCS,NNDL_SAS,TDL水文保留 率隨距研究區西邊界距離的折線統計Fig.7 Hydrological retention rate of NNDL_HPCS，NNDL_SAS， and TDL in the study area is broken-line statistics with the distance from the western boundary of the study area

如圖7所示，以研究區西邊界為起點，以PNL為基準，由西向東逐列分別統計NNDL_HPCS,NNDL_SAS,TDL相對于PNL的水文保留率。整體上，NNDL_HPCS,NNDL_SAS和TDL三者水文保留率分別為73.53%,70.68%和65.00%(不包括礦坑部分)。其中，NNDL_HPCS水文保留率最高，相較NNDL_SAS與TDL分別高2.85%與8.53%，在邊坡朝向上更接近于PNL，水文保留效果相較最優。隨距研究區西邊界距離變化上，與高程方差類似，在0～810 m段，NNDL_HPCS,NNDL_SAS和TDL三者均位于礦坑位置，水文保留率均為1；在810～1 800 m段和4 590～4 890 m段，NNDL_HPCS列水文保留率(59.19%和60.60%)低于NNDL_SAS和TDL，其均值分別較后2者均值低23.27%與9.31%，和11.84%與5.23%；在1 800～4 590 m段，NNDL_HPCS列水文保留率(79.44%)高于NNDL_SAS和TDL，其均值分別較后者均值高13.75%和16.52%。與上述高程方差類似，在“凸包地形”段，NNDL_HPCS與PNL柵格水文保留率弱于NNDL_SAS和TDL，在“凹陷地形”段則優于兩者(圖6)。

在土方運移上，依據采復子區土方對應關系，運用Matlab計算獲取NNDL_HPCS,NNDL_SAS及TDL三者的土方平均運距(PNL不存在采挖，故不進行計算)。計算結果表明，3者土方平均運距NNDL_HPCS(863.23 m/m)>NNDL_SAS(822.20 m/m)>TDL(821.71 m/m)。其中，NNDL_SAS土方平均運距較TDL僅增加約0.06%(0.49 m/m)，NNDL_HPCS則較前者約增加5.05%(41.03 m/m)，其施工成本勢必較前2者有所增加。

3.2.3 土壤抗水蝕能力求解結果與分析

地貌穩定性主要受氣候、地質、植被等諸多要素的影響。其組成物料強度直接影響邊坡的穩定性，是其抗侵蝕能力的主要指標。根據地質調查報告，確定研究區內排土場結構及巖土物理力學參數，并結合中國天氣網(http://www.weather.com.cn/)獲取研究區2007—2017年逐日氣象參數。通過CLiDE演化模型，分別模擬NNDL_HPCS,NNDL_SAS,TDL,PNL四者10 a的土壤水蝕過程。其中，為保證模擬結果的可比性，使NNDL_HPCS,NNDL_SAS和TDL三者的采復子區空間位置保持一致，最終輸入范圍如圖6所示。同時，為避免邊緣效應對模型模擬結果的影響，將研究區周邊1 000 m范圍共同納入運算(圖6)。最終通過Matlab軟件進行疊加對比分析，統計研究區內排土場NNDL_HPCS,NNDL_SAS,TDL和PNL四者10 a的土壤水蝕總量(表1)。

表1 研究區NNDL_HPCS,NNDL_SAS,TDL,PNL 土壤10 a總水侵蝕量Table 1 10-year total soil water erosion scale of the NNDL_HPCS，NNDL_SAS，TDLand PNL of the dump site

由表1可知，整體上，本研究內排土場10 a土壤水蝕總量在4.098 3×10～23.681 5×10m。其中，PNL土壤抗水蝕能力最強，模擬情景下其土壤10 a水蝕總量僅占4者平均水平(12.175 6×10m)的33.66%。在人工地貌設計下，本研究NNDL_HPCS和已有研究NNDL_SAS相較TDL分別可減少約39.07%和72.58%的土壤水侵總量，內排土場土壤抗水蝕能力提升明顯。并且，相較于已有研究中NNDL_SAS土壤抗水蝕效果，本研究NNDL_HPCS可進一步減少約54.99%的土壤水蝕總量。

如圖8，9所示，以研究區內排土場西邊界為起點，以PNL10 a土壤水蝕為基準(將其列統計值定義為對應列的0值)，由西向東逐列分別統計NNDL_HPCS,NNDL_SAS,TDL相對于PNL的10 a列土壤水蝕量。與水文保留率和絕對高程方差的統計結果不同，在810～1 800 m段，10 a土壤水蝕量中NNDL_SAS(-6.78×10m)

圖8 NNDL_HPCS,NNDL_SAS,TDL相對于PNL的土壤水蝕量隨距研究區西邊界距離變化Fig.8 Soil water intrusion of NNDL_HPCS，NNDL_SAS，and TDL relative to PNL varies with the distance from the west boundary of the study area

圖9 NNDL_HPCS,NNDL_SAS,TDL和PNL的 土壤水蝕量隨距研究區西邊界距離變化直方圖Fig.9 Histogram of soil water erosion of NNDL_HPCS， NNDL_SAS，TDL and PNL with the distance from the western boundary of the study area

即整體上，研究區PNL在全段均發生較低水平的土壤水蝕；NNDL_SAS土壤水蝕主要發生在初始段(4 590～4 890 m段)與中段(1 800～4 590 m段)，土壤堆積主要發生在終段(810～1 800 m段)；NNDL_HPCS土壤水蝕主要發生在起始段與終段，土壤堆積主要發生在中段；TDL則在全段均發生較高水平的土壤水蝕。

4 討 論

煤炭作為我國重要的能源礦產，其露采比例隨社會經濟的發展不斷增加。內排土場作為煤炭露采區的重要組成部分，其地貌形態不僅直接影響礦區人員設施及生產安全，也直接和間接影響區域生態環境系統穩定。受傳統復填方式影響，內排土場往往表現為“梯田式”人工堆墊地貌，從而使原地表要素及水文特征喪失，并在集中降雨時容易引發嚴重的表土水蝕等問題。諸多實踐與研究均表明，NbS理念下基于自然穩定特征的排土場地貌重塑具有重要的理論與實踐意義，其在有效提高區域表土穩定性的同時所需維護費用較低。與上述結論相同，團隊前置研究均表明，排土場近自然地貌的表土穩定性優于傳統“梯田式”地貌。且前置研究中，基于SAS的采排復一體下內排土場NNDL較TDL 10 a水蝕量可減少約39.07%，初步驗證了NNDL_SAS在理論上具備可行性。在此基礎上，本研究進一步提出了基于HPCS的內排土場地貌近自然重塑方法。與前者類似，本文所提方法綜合考慮了采排復一體化下土方平衡約束、區域坡度平滑控制與邊界地貌融合3方面約束，在地貌融合效果上不僅優于NNDL_SAS，還進一步提高了研究區的抗水蝕能力(10 a土壤水蝕量分別較NNDL_SAS和TDL減少約55%和72.58%)。然而與前者不同的是，本研究基于HPCS的內排土場地貌近自然重塑方法不再僅著眼于地貌重塑結果與周邊自然地貌間的銜接，而是致力于將復填可用土方再分配于原有的各個水文細小單元之中，在灰色系統理論支撐下通過保留或微調整此類要素以盡可能維護區域原有生態系統。

在露天礦采排復一體化過程中，已有坡面重塑方法能夠較好地應用于傳統內排土場復填區地貌重塑，并在縮短土方運距、減小占地面積、提升技術適用性等方向得到了充分發展。然而在地貌形態上，上述內排土場多呈現出“棱角分明”的“梯田式”空間布局，存在邊坡表土易發侵蝕，景觀明顯破碎化等問題。雖然后續部分學者提出在原有施工工藝的基礎上通過“弱化”排土場斜坡與平臺間“分界”的方式提高區域表土穩定性及其融合效果，但實際上仍是“去除棱角”的“臺階狀”地貌，未有效解決上述問題。并且，受采復周期、可采煤層厚度等諸要素的影響，復填區可用土方隨采復周期動態變化，各周期間內排形態存在差異。本研究和前置研究中，筆者均采用采復子區一一對應的方式簡化采復過程中復填可用土方的動態變化問題，并相較于TDL在土方運距小幅增加(0.06%～5.05%)的情況下大幅減少(39.07%～72.58%)表土水蝕量和獲得更高的地貌融合效果。然而值得注意的是，相較于傳統規則狀內排土場地貌，NNDL形態較為復雜。在實際施工作業中，為避免近自然重塑地貌不規則起伏所帶來的安全風險，應在排土場次高層平盤上運用卡車運輸加單斗挖掘的形式進行近自然地貌重塑，其余臺階均應嚴格遵照傳統內排土場地貌設計規范；同時應充分運用RTK等定位技術，確保排土場近自然重塑地貌特征要素(例如2.1.3節內溝道、坡定點)位置符合設計，通過“點線構面”的思路指導地貌近自然重塑。

受地貌發育影響，自然原始地貌并不總是穩定的(例如黃土高原沖刷地貌)。在特定條件下，傳統復墾后礦區或將改良原有地貌。由于本研究HPCS模型以PNL為學習對象，因此在使用該模型前首先需保證PNL地貌穩定(例如本研究PNL 10 a的土壤侵蝕量最小，僅為TDL的28.4%)。并且，本研究各樣本10 a土壤水蝕模擬量隨距西邊界距離變化的統計結果表明，在局部地區基于PNL的NNDL構建結果表土抗蝕能力甚至會劣于傳統方式(圖9中，1 800～4 590 m段，NNDL_SAS劣于TDL；4 590～4 890 m段，NNDL_HPCS劣于TDL)。因此，需充分論證NNDL,TDL和PNL三者的穩定性，從而確定有利于因地制宜的內排土場地貌重塑方式。

5 結 論

(1) 地貌設計上，基于HPCS的內排土場近自然地貌重塑模型可依據預設采復周期下采復子區采前地表、礦層頂板和礦層底板3者的DEM，自動求解復填可用土方量，獲取研究區NNDL_HPCS。

(2) 地貌融合上，NNDL_HPCS較NNDL_SAS和TDL視覺上更接近于PNL，其相對PNL的高程總方差較后2者分別低50.37%和37.92%，且其自然水文保留率較后2者分別高出2.85%和8.53%。

(3) 土方調配上，NNDL_HPCS土方運移成本小幅增加，其土方平均運距分別較NNDL_SAS和TDL分別高出約4.99%和5.05%。

(4)土壤抗水蝕上，NNDL_HPCS較NNDL_SAS 10 a可減少約55%的表土水蝕總量，較TDL減少約72.6%，其可顯著提高內排土場表土抗水蝕能力，與PNL相逼近。