基于自然邊坡模型的草原煤礦排土場坡形重塑

李 恒，雷少剛，黃云鑫，程 偉，鄧 彪，張周愛

(1.中國礦業大學 礦山生態修復教育部工程研究中心，江蘇 徐州 221116; 2.神華北電勝利能源有限公司，內蒙古 錫林浩特 026015; 3.神華寶日希勒能源有限公司，內蒙古 呼倫貝爾 021500)

采礦活動會對生態環境造成極大的破壞，主要體現在改變了礦區原有地表形態和自然景觀，從而引發一系列的土地與環境問題，因此將受采礦擾動的區域一定程度地恢復到擾動前的生態服務功能，使之能與周圍的自然環境相協調是十分重要的[1-2]。在土地復墾與生態重建過程中，地形重塑、土壤重構和植被重建被認為是3個十分重要的環節，其中地形重塑為基礎[3]。但以往的土地復墾方案更重視土壤重構[4-5]和植被重建[6]方面，而地形重塑受制于理論和技術層面的問題，往往會被研究者所輕視。在這種情況下，排土場往往被重塑為簡單的階梯式地形，未能與當地自然環境相融合，不僅養護成本較高，且缺乏長期穩定性，從而產生較為嚴重的水土流失和景觀破碎現象。因此尋找合理的地形重塑理論依據和方法顯得十分重要。

近自然地形重塑作為一種仿造當地自然形態進行地形重塑模型構建的理論，在礦區復墾活動中被國內外諸多學者應用。國外學者和研究單位在結合水文地質學[7-8]、自然地理學[9]、系統動力學等相關理論和原則的基礎上，以塑造不會輕易變動，且與周圍未損毀土地相適應的近自然地形為目的[10-11]，從流域角度開發出了GeoFluv[12]，GeoWEPP[13-15]等模型進行廢棄礦區的地形重塑。國內學者也逐漸意識到近自然地形重塑的優勢，并應用河流動力學[16-17]和流域地貌學[18-20]的相關理論和分析重新構建項目區的流域格局。現有近自然地形重塑的相關研究主要著眼于以水文為基礎的損毀地形整體構建[21]，在排土場等工地形的坡面形態方面研究較少。

近自然地形重塑的主要理論依據是“師法自然”，即學習受當地氣候、土壤、植被、水文條件影響下的自然地表形態，并應用其重塑礦區損毀地形。筆者以內蒙古勝利礦區為研究區，為設計出適宜勝利礦區自然條件的排土場邊坡模型，將研究區內未受采礦擾動的自然地表形態的關鍵特征參數化并進行統計分析，分析其特征參數之間的相關性，在此基礎上構建出與當地環境相適應，土壤損失量較小的近自然邊坡模型。

1 研究區概況

研究區勝利礦區(圖1)位于內蒙古自治區錫林浩特市西北部，錫林河以西的緩坡丘陵干草原與河谷沖積、湖積平原的過渡帶,西部屬緩坡丘陵,東部為河谷、湖積平原。氣候類型屬中溫帶半干旱、干旱大陸性季風氣候[22]。勝利礦區排土場共5層臺階，每層臺階高15 m，邊坡角為33°，排土場整體坡角為17°。由于排土場簡單的采用階梯式方法構建，長期受到雨水的沖刷侵蝕，土壤受侵蝕強度遠大于周邊自然環境，由此導致排土場邊坡已經出現溝蝕、垮塌等現象。

圖1 研究區區位Fig.1 Location of the study area

2 研究方法

2.1 DEM數據獲取

數據獲取時間為4月，進行邊坡特征參數提取的參照區位于勝利礦區內的自然區域，其中陡坡提取區域位于勝利礦區西北部，緩坡提取區域位于勝利礦區西部。在研究區未受采礦擾動的自然區域內，使用大疆精靈3無人機獲取高精度的地面數字高程模型(DEM)數據。數據范圍為3 000 m×3 000 m，像元大小為0.068 m×0.068 m。

2.2 邊坡特征參數提取及邊坡模型構建

在Arcmap10.2中對DEM數據進行填洼后進行坡向分析和獲取等高線兩部分處理，進行這兩部分處理的原因是坡向分析后的圖像可較為清晰的顯示出分水嶺，而等高線分布可以顯示出參照區域的坡度。以此為基準獲取自然地形部分的邊坡剖面線，通過邊坡剖面線提取可知勝利礦區自然地形部分的主要邊坡類型為反S型邊坡，共提取出52組反S型邊坡剖面線，其中包括28組陡坡部分的邊坡剖面線及24組緩坡部分的邊坡剖面線。

反S型邊坡以拐點為分界處分為上側的凸起部分和下側的凹陷部分。反S型邊坡與其他類型邊坡有部分共同的邊坡特征參數，包括坡高、坡長、凸面曲率和凹面曲率，但反S型邊坡與其他類型邊坡的不同之處在于其凸起部分和凹陷部分存在拐點，因此這兩部分分別在水平面和鉛垂面上有不同的占比。為此筆者提出了凸面水平占比和凸面豎直占比兩個新的邊坡特征參數，通過這兩個參數可以確定拐點在邊坡上的位置，進而可以確定反S型邊坡的具體形狀。

根據剖面線數據計算反S型邊坡的特征參數，包括坡高、坡長、凸面曲率、凹面曲率、凸面水平占比、凸面垂直占比。其中，坡高指邊坡坡面在鉛垂面上的投影距離，坡長指邊坡坡面在水平面上的投影距離，凸/凹面曲率指邊坡凸起/凹陷部分形成的圓弧的曲率，凸面水平占比指邊坡上側凸起部分在水平面上的投影長度占坡長的百分比，凸面豎直占比指邊坡凸起部分在鉛垂面上的投影長度占坡高的比例，根據這些邊坡特征參數可確定反S型邊坡的坡形。分析邊坡特征參數之間的相關性，利用統計學方法構建勝利礦區自然邊坡的邊坡模型。

2.3 邊坡土壤損失量對比

利用土壤水蝕預報模型WEPP對構建的自然邊坡模型和勝利礦區原排土場邊坡模型進行土壤水蝕分析。WEPP模型在進行模擬時需要氣象、土壤及植被覆蓋情況等基礎數據。其中氣象數據以研究區2016年的氣象數據為基礎，使用氣候發生器CLIGEN生成模型所需要的氣象數據文件。土壤數據包括沙粒含量、黏粒含量、有機質含量、石礫含量、陽離子交換量及土壤容重，將這些數據輸入WEPP中即可生成土壤數據文件。由于本次研究主要考慮反S型邊坡模型與原排土場邊坡模型之間土壤損失量的差異，為控制變量先假設邊坡坡面的植被覆蓋程度為0，即全部為裸地的情況。在確定了氣象、土壤、植被覆蓋情況等基礎數據后，通過對比不同坡高、不同模擬年份條件下反S型邊坡模型和原排土場模型的土壤損失量和土壤流失曲線圖，以驗證構建的近自然邊坡模型的適用性。

3 研究結果

3.1 邊坡特征參數相關性分析

為得到勝利礦區自然邊坡特征參數之間的擬合公式，首先需要分析這些邊坡特征參數之間的相關性。在利用SPSS進行相關性分析之前，由于選取的邊坡特征參數有可能存在偏離一般規律的數值，因此首先要對偏離值進行去除。在去除了偏離值之后，根據坡高和坡長的關系制作散點圖(圖2)，根據圖像可看出，勝利礦區緩坡部分與陡坡部分有著不同的規律，因此需要對兩部分分別進行分析。

圖2 坡高對坡長的散點Fig.2 Scatter plot and fitting curves of slope height to slope length

陡坡部分的邊坡剖面線坡度在22%～46%，對陡坡部分的邊坡特征參數進行相關性分析的結果發現，坡高與坡長、凸面水平占比與凸面豎直占比在0.01水平上顯著相關;坡長與凸面曲率、坡長與凹面曲率在0.05水平上顯著相關。

緩坡部分的邊坡剖面線坡度在6%～10%，對緩坡部分的邊坡特征參數進行相關性分析的結果發現，坡高與坡長、凸面曲率與凹面曲率、凸面水平占比與凹面豎直占比分別在0.01水平上顯著相關。

3.2 邊坡模型構建

3.2.1陡坡部分邊坡模型構建

根據陡坡部分邊坡特征參數之間的相關性分析，選取坡高對坡長、坡長對凸面曲率和凹面曲率、凸面水平占比對凸面豎直占比這4組相關特征參數進行曲線擬合，通過這種擬合得到的公式組，可以得到一個以坡高為自變量的近自然邊坡模型。在進行曲線擬合之前，為檢驗模型的準確性，將提取兩組數據作為驗證數據。在曲線擬合的過程中(圖3)，通過分析R2值，參數的T檢驗和方差的F檢驗后，分別得到了4組相關特征參數的擬合式:

s=e3.829+0.024d

(1)

qt=e2.530-0.007s

(2)

qa=e2.039-0.011s

(3)

zt=e-2.219+2.673pt

(4)

式中，d為坡高;s為坡長;qt為凸面曲率;qa為凹面曲率;pt為凸面水平占比;zt為凸面豎直占比。

圖3 陡坡部分邊坡特征參數擬合Fig.3 Slope characteristic parameter fitting diagram of steep part

其中由于凸面水平占比和凸面豎直占比與其他邊坡特征參數沒有相關性，而最終的邊坡模型只能有坡高一個自變量，為獲取特征規律應盡量采取與自然邊坡特征參數相近的數值，因此求取實測數據凸面水平占比的均值為46.82%，根據式(4)計算凸面豎直占比的參考值為38.01%。由此構建出陡坡部分的自然邊坡模型。將驗證數據的坡高代入模型中，所得到的坡長、凸面曲率和凹面曲率數據與實測數據的相差幅度均低于10%，可說明模型的準確性。

3.2.2緩坡部分邊坡模型構建

根據緩坡部分邊坡特征參數的相關性分析，分別選取坡高對坡長、凸面曲率對凹面曲率、凸面水平占比對凹面豎直占比這3組相關特征參數進行曲線擬合(圖4)。同樣提取出兩組數據作為驗證，用上述方法擬合出了3組特征參數的擬合式:

s=e3.782+0.103d

(5)

qa=e0.24-0.443qt

(6)

zt=e-1.187+0.736pt

(7)

同樣對實測數據的凸面曲率和凸面水平占比求取均值，并按照式(6)，(7)計算凹面曲率和凸面豎直占比的數值。經計算的凸面曲率參考值為1.41×10-3，凹面曲率為0.69×10-3，凸面水平占比為52.39%，凸面豎直占比為45.26%。

3.3 反S型邊坡模型與原排土場邊坡模型對比

考慮到實際情況，由于在進行排土場邊坡重塑時不可能按照緩坡的坡度進行放坡處理，而陡坡部分邊坡模型的坡度與原排土場邊坡相近，為此在對比分析時使用陡坡部分邊坡模型與原排土場邊坡模型進行分析。在進行水蝕分析的時候，為確保分析結果的準確性，應當對除坡型因素外的其他變量進行控制，主要包括邊坡高度，土壤類型，氣候數據，植被覆蓋以及模擬年份。為了使分析結果更加符合原排土場的實際情況，在進行水蝕分析模擬時，應當以一個或者數個臺階的高度作為基準考慮。在綜合考慮邊坡設計需要的基礎上，筆者最終選取邊坡高度分別為15，30，45和60 m來構建不同的邊坡模型。而在模擬年份這一變量中，由于邊坡形狀或許會因為常年受到侵蝕而發生改變，導致土壤損失量發生變化，因此需要劃分不同的模擬年份。本次分析將以1，10和50 a分別進行水蝕分析。分析結果見表1。

圖4 緩坡部分邊坡特征參數擬合Fig.4 Slope characteristic parameter fitting diagram of gentle part

表1 不同邊坡的特征參數及年均土壤損失量
Table 1 Characteristic parameters annual soil loss of each slope

參數原排土場邊坡反S型邊坡坡高/m1530456015304560坡長/m56.86113.72170.58227.4465.9694.54135.50194.22邊坡特征參數凸面曲率/10-3————7.916.484.863.22凹面曲率/10-3————3.722.721.730.91凸面水平占比/%————46.8246.8246.8246.82凸面豎直占比/%————38.0138.0138.0138.011 a23.9348.5072.1495.3917.0031.8642.2750.64土壤損失量/(kg·m-2)10 a32.0765.4696.40127.9021.6741.1854.0559.9150 a35.0071.78106.30141.2023.3544.3258.7065.65

由表1可知，在模擬坡高為15 m時，反S型邊坡1，10，50 a的土壤損失量分別為原排土場邊坡的71.04%,67.57%,66.71%；在模擬坡高為30 m時，反S型邊坡1,10,50 a的土壤損失量分別為原排土場邊坡的65.69%,62.91%,61.74%；在模擬坡高為45 m時，反S型邊坡1,10,50 a的土壤損失量分別為原排土場邊坡的58.59%,56.07%,55.22%；在模擬坡高為60 m時，反S型邊坡1,10,50 a的土壤損失量分別為原排土場邊坡的53.09%,46.84%,46.49%。

在對各個坡型的年均土壤損失量進行了對比分析后，為了對不同坡型土壤流失的分布情況進行對比，還需要根據兩種坡型的土壤流失曲線圖進行分析，圖5分別為模擬年份為1 a，邊坡高度為30 m的條件下兩種坡型的土壤流失曲線圖。其中紅色的曲線表示坡面剖面線，綠線表示坡面上相應點的侵蝕情況，灰色的面積圖形在Y軸上的數值為實際上的土壤流失量或土壤沉積量。

圖5 原排土場邊坡和反S型邊坡的土壤流失曲線Fig.5 Soil loss curves of the slope of the original dump and anti-s slope

根據原排土場邊坡模型的土壤流失曲線圖(圖5(a))可知，原排土場邊坡的土壤流失從第1坡面開始，而在第2平臺上有土壤顆粒的沉積，由此可看出平臺設置對邊坡保土能力有一定的影響。但在第2坡面位置，土壤損失量比第1坡面位置提升了一倍左右，而最終在第2坡面底部土壤損失量達到最大，最大分離量為116 kg/m2。

根據反S型邊坡模型的土壤流失曲線圖(圖5(b))中可知，邊坡的土壤損失量從坡頂開始增長，在到達邊坡長度50 m左右的時候在47 kg/m2處上下波動并保持到坡底。最大分離量為48.4 kg/m2。而在邊坡長度為51.5 m處的位置土壤損失量有一個驟減，之后繼續增加。由于此位置接近反S型邊坡凸起部分和凹陷部分的拐點位置，因此推測出現此峰值的原因是此位置剖面曲率趨于平緩，地表徑流流速減慢，導致地表徑流中攜帶的土壤顆粒得以沉降，在此處部分堆積使得土壤損失量減少。

4 討 論

4.1 反S型邊坡與其他類型邊坡的對比

由于自然地形的復雜性，在勝利礦區未受采礦擾動的自然區域內對邊坡剖面線進行提取時，獲得的坡型并不只有反S型邊坡一種。陡坡部分在提取過程中總共獲得60條邊坡剖面線(圖1)，其中反S型邊坡28條，凸型邊坡21條，S型邊坡6條，凹型邊坡4條，直線型邊坡1條，反S型邊坡和凸型邊坡分別占邊坡提取總數的46.67%和35%。緩坡部分在提取過程中總共獲得55條邊坡剖面線，其中反S型邊坡24條，凸型邊坡15條，S型邊坡3條，凹型邊坡6條，直線型邊坡7條，反S型邊坡和凸型邊坡分別占邊坡提取總數的43.63%和27.28%。由于后3種邊坡在本次提取中獲得的數量過少，可以認為其并非為勝利礦區自然邊坡的主要坡型。而在對凸型邊坡進行特征參數提取和邊坡模型構建并進行水蝕分析后，發現該邊坡模型的土壤損失量高于反S型邊坡，因此采用反S型邊坡模型對排土場邊坡進行設計。

在對陡坡部分反S型邊坡進行模型構建時，統計實測數據的凸面水平占比并求取均值為46.82%，根據凸面水平占比和凸面豎直占比的擬合式(4)，計算得出凸面豎直占比為38.01%，即凹面部分的豎直占比為61.99%，十分接近61.80%。這表明陡坡部分的反S型邊坡模型，其凸面向凹面過渡的拐點在坡高方向上的投影位置在黃金分割點左右。據此可以推測，自然界的反S型邊坡，其凹面部分和凸面部分在豎直方向上的拐點位于黃金分割點附近，蘊含著豐富的美學價值。

4.2 反S型邊坡的保土能力和長期穩定性

保土能力是評價邊坡模型優劣的重要依據。邊坡的保土能力越強，土壤損失量越少，其穩定性就越高。隨著坡高的增加，無論是原排土場邊坡模型還是反S型邊坡模型的年均土壤損失量都有一定程度的提升。原排土場邊坡模型在坡高分別為30,45,60 m時的1 a年均土壤損失量分別比15 m時的1 a年均土壤損失量增加了102.67%,201.46%,298.62%，而反S型邊坡模型同等條件下的增加百分比分別為87.41%,148.65%,197.88%，相較于原排土場邊坡模型分別減少了15.26%,52.82%,100.74%，并且原排土場邊坡模型的年均土壤損失量增幅差異不大，而反S型邊坡模型的年均土壤損失量增幅呈現遞減趨勢。這表明隨著坡高的增加，反S型邊坡模型的保土能力相較于原排土場邊坡模型越來越強。

另一方面，隨著模擬年份的增加，兩種邊坡模型的年均土壤損失量都有一定程度的提升，這可能是因為邊坡在經歷數年的侵蝕導致邊坡形狀發生改變后，其保土能力有所下降，因此邊坡模型的長期穩定性也是一個重要的評價依據。隨著模擬年份的增加，邊坡模型在不同坡高條件下的年均土壤損失量增加百分比相差不多，如原排土場模型在坡高為15,30,45,60 m的10 a年均土壤損失量相較于1 a年均土壤損失量增加百分比分別為34.02%,34.97%,33.63%,34.08%，故可以用均值代替。原排土場邊坡模型各坡高條件下的10 a年均土壤損失量相較于1 a年均土壤損失量提升了34.17%左右，50 a年均土壤損失量相較于1 a年均土壤損失量提升了47.41%左右;而反S型邊坡模型分別提升了28.2%和38.44%左右。這表明反S型邊坡模型在長時間尺度下的保土能力優于原排土場邊坡模型。

4.3 近自然邊坡模型的應用和局限

地面坡度越小，地表越平緩，其受侵蝕的影響效果也就越低[23]，從這種角度而言使用緩坡部分的邊坡模型貌似是最佳選擇。但考慮到實際情況，現有排土場的坡角為17°，若把坡度放緩到緩坡部分(大多<10%)，那將占用坡底部分極多的土地，與實際情況不符。而陡坡部分的反S型邊坡模型，其坡角與原排土場的坡角相近，但其保土能力和長期穩定性的提升效果十分顯著，故使用陡坡部分的反S型邊坡模型對勝利礦區排土場進行設計是可行的。

本文構建邊坡模型的方法不僅適用于排土場邊坡，還可用于自然地形和設計地形的銜接區域。長久以來，各近自然重塑研究都著眼于以流域為主的整體地形構建[24]，而構建形成的地形往往與自然地形不能很好的進行銜接。本文使用的邊坡模型構建方法為解決設計地形與自然地形的銜接問題提供了一種思路，即根據研究區周邊自然地形的邊坡特征參數，采用統計學方法構建出自然地形和設計地形的過渡區，使得設計地形能較好的與自然地形進行銜接，與當地的流域地貌相協調。

另一方面，本文構建邊坡模型的方法也可用于礦區復墾規劃中的排土場設計。根據開采前獲取的埋深、采深、采厚等地質信息估算外排土場的土方量，計算松散系數和下沉系數后設計外排土場高度，利用本文所構建的近自然邊坡模型對外排土場邊坡進行優化設計。在對內排土場進行設計時，由于采掘區域的煤層埋深和采深一般具有空間異質性，可以利用本文構建的近自然邊坡模型將內排土場重塑為近自然地形而非傳統復墾模式中的平地，并對內、外排土場與周圍地貌的銜接區域進行優化，將具有一定高度差的銜接區域設計成近自然的反S型邊坡，使其能與周圍的自然環境相協調。

本文對邊坡模型保土能力的評價是基于WEPP土壤水蝕預報模型進行的。該模型的優點是能在短時間內模擬不同時間尺度下的土壤侵蝕過程，相較于野外不同時段實測來進行對比的方法而言，簡化了研究過程并降低了研究難度。但應用WEPP模型進行評價的方法存在一定的局限性。一方面，WEPP模型是基于美國本土環境開發的水蝕模型，因此在很多文獻中提到過WEPP模型的模擬過程中有部分經驗參數，雖然通過氣候、土壤等數據能減弱這些經驗參數的影響，但還是有可能造成模擬結果上的誤差。另一方面，WEPP模型是基于水力侵蝕進行的預報模型，對以水力侵蝕為主的勝利礦區而言適用性較好，但對于風力侵蝕、凍融侵蝕影響效果遠大于水力侵蝕的區域，WEPP模型的預測結果可能與實際結果有偏差。

5 結 論

(1)勝利礦區自然區域內邊坡類型主要為反S型邊坡。通過對邊坡形態進行分析，提出了反S型邊坡不同于其他類型邊坡的兩個邊坡特征參數，即凸面水平占比和凸面豎直占比。分別對陡坡部分和緩坡部分的反S型邊坡的邊坡特征參數進行相關性分析發現，陡坡部分邊坡特征參數的坡高與坡長、凸面水平占比與凸面豎直占比、坡長與凸面曲率、坡長與凹面曲率顯著相關;緩坡部分邊坡特征參數的坡高與坡長、凸面曲率與凹面曲率、凸面水平占比與凹面水平占比顯著相關。

(2)在對邊坡特征參數進行擬合后分別得出陡坡部分和緩坡部分的反S型邊坡模型，并發現陡坡部分的反S型邊坡模型，其凸面和凹面在豎直方向上的占比符合黃金比例。考慮到實際情況，確定了以陡坡部分反S型邊坡模型與原排土場邊坡模型進行對比分析。

(3)在WEPP模型中輸入勝利礦區的氣候、土壤、植被覆蓋等數據，對不同坡高、不同模擬年份條件下反S型邊坡模型和原排土場邊坡模型進行了水蝕分析，得到了兩種邊坡模型的土壤損失量及土壤流失曲線圖。結果表明隨著坡高和模擬年份的增加，反S型邊坡模型的保土能力和長期穩定性明顯優于原排土場邊坡模型，驗證了反S型邊坡模型在排土場邊坡設計中的適用性。

(4)本文基于自然邊坡特征參數構建近自然邊坡模型的方法不僅適用于排土場坡形設計，也為解決礦區土地復墾活動中設計地形與自然地形的銜接問題提供了一種思路，可作為其他地區進行地形重塑的理論依據。