蔡懿軒， 楊 剛?， 張成梁， 顧清敏

(1.北京林業(yè)大學(xué)信息學(xué)院，100083，北京； 2.國(guó)家林業(yè)草原林業(yè)智能信息處理工程技術(shù)研究中心，100083，北京；3.輕工業(yè)環(huán)境保護(hù)研究所，100089，北京； 4.國(guó)家能源集團(tuán)寧夏煤業(yè)有限責(zé)任公司羊場(chǎng)灣礦區(qū)，751410，寧夏靈武)

礦區(qū)廢棄地的生態(tài)修復(fù)中，地形重塑是第一步、也是最核心的部分。地形構(gòu)建不合理將導(dǎo)致礦區(qū)嚴(yán)重的水土流失，使后期的植被重建受阻，從而導(dǎo)致極大的維持和養(yǎng)護(hù)成本。近年來(lái)，近自然地形重塑的理念和方法逐漸被提出并應(yīng)用。大自然具有自我更新和再生能力，其自身原始的地形特征為生態(tài)重建提供最重要的參考依據(jù)。近自然地形重塑就是以礦區(qū)周邊未被擾動(dòng)的自然地形為參照來(lái)進(jìn)行地形設(shè)計(jì)。因此，充分了解當(dāng)?shù)匚幢粩_動(dòng)的自然區(qū)域的地形、水文等特征已成為近自然地形重塑的先決條件。而在此過(guò)程中，存在一系列有待探討的問(wèn)題：1)參照區(qū)域的選擇問(wèn)題，即選擇哪塊兒區(qū)域作為近自然地形重塑的參照區(qū)域；2)參數(shù)的選擇問(wèn)題，即選擇哪些地形特征參數(shù)進(jìn)行計(jì)算和分析；3)參數(shù)的計(jì)算。

近年來(lái)，國(guó)內(nèi)外許多學(xué)者對(duì)近自然地形重塑方法進(jìn)行探討[1-3]：Terrence等[4]針對(duì)礦區(qū)廢棄地的坡面和溝道治理，以集水區(qū)為單元討論地形重塑方法。楊翠霞等[2]探討露天開采礦區(qū)廢棄地近自然地形重塑的理論體系和技術(shù)方法，并結(jié)合Geofluv模型重建周口店采石場(chǎng)廢棄地。Bugosh等[5]基于流域地貌學(xué)理論提出地形重塑模型GeoFluv。但這些工作重點(diǎn)放在地形的流域水力侵蝕作用上，且大多僅考慮流域尺度的參數(shù)計(jì)算。其對(duì)于參照區(qū)域選擇、參數(shù)選擇等問(wèn)題尚未有系統(tǒng)工作，還未形成成熟的解決方案。本研究以寧夏寧東羊場(chǎng)灣礦區(qū)一區(qū)排矸場(chǎng)為實(shí)例，討論近自然地形重塑中的參照區(qū)域確定、特征參數(shù)選取及參數(shù)計(jì)算問(wèn)題，以給出可行解決方案，其分析結(jié)果可為該礦區(qū)排矸場(chǎng)近自然地形重塑提供參照目標(biāo)，并為該地區(qū)地理分析、環(huán)境評(píng)價(jià)及生態(tài)環(huán)境治理等提供基礎(chǔ)數(shù)據(jù)。

1 研究區(qū)概況

羊場(chǎng)灣礦區(qū)位于寧夏回族自治區(qū)靈武市寧東鎮(zhèn)境內(nèi)(圖1a)，為我國(guó)寧東能源化工基地的主力礦井之一。該地區(qū)位于毛烏素沙漠的邊緣地帶，屬溫帶干旱氣候區(qū)，典型的大陸性氣候，多年平均降雨量212 mm，年平均蒸發(fā)量2 862.2 mm，研究區(qū)全年盛行西北風(fēng)，年平均風(fēng)速3.1 m/s，全年>17 m/s的大風(fēng)超過(guò)50 d。研究區(qū)地勢(shì)較為平緩，屬低緩剝蝕殘丘地貌。研究區(qū)內(nèi)土壤主要以淡灰鈣土為主，部分地區(qū)有明顯土壤沙質(zhì)化現(xiàn)象，地表抗侵蝕能力弱；常年大風(fēng)，風(fēng)力對(duì)地表土壤的搬運(yùn)作用和塑形作用明顯。

羊場(chǎng)灣礦區(qū)平均年排矸量可達(dá)38萬(wàn)t。排矸場(chǎng)位于礦區(qū)東南處，主要依據(jù)地形自上而下傾倒的方式進(jìn)行排矸，目前已形成約10 m高差的陡坡，并繼續(xù)向東側(cè)推進(jìn)。矸石山坡度達(dá)40°，邊坡易發(fā)生侵蝕現(xiàn)象；在極端氣象條件下，易發(fā)生山體滑坡、塌方。隨著矸石堆放的增多，對(duì)當(dāng)?shù)赝恋刭Y源、大氣環(huán)境、水資源等有強(qiáng)烈擾動(dòng)。

2 研究方法

2.1 數(shù)據(jù)準(zhǔn)備

本研究采用資源三號(hào)衛(wèi)星拍攝5 m分辨率DEM數(shù)據(jù)為基礎(chǔ)數(shù)據(jù)(圖1b)。從資源三號(hào)衛(wèi)星的一景數(shù)據(jù)中選取礦區(qū)周邊區(qū)域87.33 km2的數(shù)據(jù)作為試驗(yàn)區(qū)，數(shù)據(jù)獲取時(shí)間為2019年8月，地圖數(shù)據(jù)范圍為： E 106°32′～106°41′、N 37°57′～38°03′。圖中待重塑區(qū)劃定范圍為羊場(chǎng)灣礦區(qū)一區(qū)排矸場(chǎng)區(qū)域。

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

2.2 參照區(qū)域的選取

在近自然地形重塑時(shí)，需根據(jù)礦區(qū)待重塑區(qū)域的情況選擇與之有類似特征、對(duì)其重塑有參照價(jià)值的未干擾區(qū)域作為參照目標(biāo)。

一般而言，流域侵蝕作用對(duì)地貌形成影響顯著，與礦區(qū)的地質(zhì)安全息息相關(guān)。研究區(qū)屬西北干旱荒漠區(qū)，年平均降雨量很低，但其降雨非常集中，偶發(fā)性暴雨對(duì)地面沖蝕作用非常強(qiáng)烈，流域侵蝕作用依然是影響其地表形態(tài)的最重要因素。因此，在進(jìn)行參照區(qū)域的選取和分析時(shí)，依據(jù)流域的鄰近性、相似性等因素，并以流域?yàn)檫x取單位進(jìn)行參照區(qū)域的確定；同時(shí)兼顧氣候、土壤、巖性等其他自然因素。

具體而言，首先通過(guò)水文分析得到相關(guān)區(qū)域的流域信息，在此基礎(chǔ)上依據(jù)如下原則進(jìn)行參照區(qū)域選取：

1) 鄰近性原則。選取待重塑區(qū)的鄰近流域?yàn)閰⒄諈^(qū)域；

2) 近似性原則。選擇與待重塑區(qū)在水文、氣候、土壤、植被和巖性等相似的流域區(qū)域；

3) 穩(wěn)定性原則。選擇處于內(nèi)應(yīng)力和外應(yīng)力基本平衡的穩(wěn)定狀態(tài)的參照區(qū)域。

在此基礎(chǔ)上，還需對(duì)所選取流域的相關(guān)特征進(jìn)行篩查，剔除明顯偏離整體樣本數(shù)據(jù)的“異?！绷饔颉Ｆ渲胁捎孟渚€圖方法[6]進(jìn)行異常流域的剔除。

2.3 參數(shù)的選擇

與地形特征相關(guān)的參數(shù)有許多[7-9]，但需重點(diǎn)選取對(duì)近自然地形重塑有參考意義的參數(shù)進(jìn)行計(jì)算和分析。本研究主要參考楊翠霞等[10]給出的近自然地形重塑指標(biāo)及美國(guó)GeofluvTM模型[11]中的地形計(jì)算參數(shù)。楊翠霞等[10]綜合前人工作，選取流域面積、流域圓度、溝道密度等7項(xiàng)指標(biāo)參數(shù)來(lái)刻畫流域地貌的多方面形態(tài)特征，這些參數(shù)與礦區(qū)地形重塑設(shè)計(jì)密切相關(guān)，因此選取這些參數(shù)作為自然流域地貌特征的指標(biāo)。

上述流域地貌特征重點(diǎn)在于流域宏觀形態(tài)的把握，而在進(jìn)行近自然溝道設(shè)計(jì)時(shí)，本研究重點(diǎn)參考了美國(guó)GeoFluvTM模型中有關(guān)溝道形態(tài)和局部地表形態(tài)的若干參數(shù)。GeoFluvTM模型是Bugosh等[5]開發(fā)的用于廢棄礦區(qū)地形重建的數(shù)學(xué)模型。該模型可根據(jù)給定的溝道參數(shù)和地形參數(shù)來(lái)模擬生成近自然地形。

流域地貌特征是以流域?yàn)閱卧M(jìn)行地形特征的描述，其測(cè)量尺度一般為數(shù)km2到數(shù)百km2。但地表更小尺度內(nèi)的局部形狀特征同樣影響自然地形的穩(wěn)定性，對(duì)近自然地形重塑具有重要參照作用。為此，選取坡度、坡向、局部地形起伏度、局部地表粗糙度等4項(xiàng)參數(shù)來(lái)刻畫局部地表形態(tài)特征。其中坡度、坡向?yàn)槠旅嫣卣?；局部地形起伏度、局部地表粗糙度分別刻畫地形在豎直方向、水平方向上的變化特征。

綜上所述，擬將測(cè)量參數(shù)歸為3類，共14項(xiàng)參數(shù)：1)流域地貌形態(tài)：溝道級(jí)別、流域面積、主溝道長(zhǎng)度、溝道密度、流域圓度、流域高程差、平均坡度；2)溝道形態(tài)：蜿蜒度、溝道比降、分水嶺到溝道源頭最小距離；3)局部地表形態(tài)：坡度、坡向、局部地形起伏度、局部地表粗糙度。

2.4 參數(shù)的計(jì)算

擬計(jì)算的地形參數(shù)涉及地形多種類型的形態(tài)特征，需采用不同的計(jì)算方法，且有計(jì)算上的先后次序。為此我們制定系統(tǒng)的計(jì)算路線(圖2)：首先進(jìn)行水文分析并確定參照區(qū)域；其次分2個(gè)分支進(jìn)行參數(shù)計(jì)算，左側(cè)分支進(jìn)行流域地貌特征和溝道形態(tài)特征的計(jì)算，右側(cè)分支對(duì)局部地表形態(tài)特征進(jìn)行計(jì)算。

圖2 參數(shù)計(jì)算的技術(shù)路線Fig.2 Technical route of parameter calculation

在“匯流量閾值的確定”步驟中首先擬合出溝道密度與匯流量閾值間的關(guān)系曲線，其次采用均值變點(diǎn)分析法[12]找到拐點(diǎn)來(lái)確定河網(wǎng)匯流量閾值。均值變點(diǎn)法是目前最佳分析窗口計(jì)算中的主流方法，相比人工作圖法等方法，其分析效果更符合實(shí)際。在3.3.3節(jié)中局部地形起伏度計(jì)算中，采用均值變點(diǎn)法來(lái)確定最佳分析窗口尺寸。

“分水嶺到溝道源頭最小距離”這項(xiàng)參數(shù)來(lái)自于GeoFluv模型[11]。當(dāng)降水發(fā)生后，坡面上部(靠近分水嶺的部分)不會(huì)立刻形成侵蝕的溝道，而是需要經(jīng)過(guò)一定距離的流量匯集，對(duì)坡面形成足夠的侵蝕力形成溝道。參數(shù)“分水嶺到溝道源頭最小距離”就是為反映其現(xiàn)象。該參數(shù)與降水量、降水強(qiáng)度、地表粗糙度、土壤的抗沖性和滲透性都有關(guān)系，是地形重塑中近自然溝道設(shè)計(jì)的重要參考值[11]。首先基于地形位置指數(shù)[13]來(lái)識(shí)別地形類別，由此確定分水嶺區(qū)域；其次針對(duì)每個(gè)子流域確定入水口；之后利用ArcGIS測(cè)量分水嶺到溝道源頭最小距離值。

3 結(jié)果與分析

3.1 流域分析及參照區(qū)域的確定

3.1.1 匯流量閾值分析及河網(wǎng)提取 以2000柵格為間隔，提取不同匯流量閾值下的河網(wǎng)，并計(jì)算其溝道密度，得到溝道密度與匯流量閾值間的關(guān)系曲線(圖3)。

圖3 溝道密度與匯流量閾值擬合關(guān)系Fig.3 Fitting relationship between gully density and catchment flow threshold

從圖3中看出，隨著匯流量閾值的增大，溝道密度逐漸遞減趨于平緩。其擬合曲線的相關(guān)性系數(shù)R趨于1，符合冪函數(shù)關(guān)系。根據(jù)均值變點(diǎn)分析法[12]，計(jì)算原始樣本的統(tǒng)計(jì)量S與樣本分段后的統(tǒng)計(jì)量Si，S與Si之間最大差值對(duì)應(yīng)的點(diǎn)即為變點(diǎn)。觀察S與Si的差值變化曲線(圖4)，當(dāng)閾值為4 000時(shí)差值達(dá)到最大，因此確定河網(wǎng)提取的最佳匯流量閾值為4 000。

利用ArcGIS得到閾值4 000河網(wǎng)提取結(jié)果(圖5)。其中河網(wǎng)中河道總長(zhǎng)度為228.04 km，河道數(shù)目為438條。

圖5 河網(wǎng)提取結(jié)果(匯流量閾值為4 000)，粉紅色框區(qū)域則是3.2節(jié)分析后所選定的參照區(qū)域Fig.5 River network extraction results (catchment flow threshold is 4 000), the pink framed area is the reference area selected after the analysis in Section 3.2

3.1.2 流域分析及參照流域的確定 根據(jù)圖5給出的待重塑區(qū)及其周邊的溝道狀況，對(duì)溝道進(jìn)行分級(jí)[14]。圖中排矸場(chǎng)西側(cè)存在部分人工區(qū)域，因此主要從排矸場(chǎng)東側(cè)的未擾動(dòng)的自然流域中選擇參照流域進(jìn)行分析。

根據(jù)鄰近性原則，重點(diǎn)考慮與待重塑區(qū)鄰近并相交的2個(gè)流域區(qū)域：待重塑區(qū)北部南向北的流域(圖6中紫色流域1)與待重塑區(qū)南部從東向西的流域(圖6中綠色流域2)。待重塑區(qū)北部所覆蓋的Ⅰ級(jí)溝道由南向北匯合成1個(gè)Ⅱ級(jí)溝道，并與右側(cè)的另一個(gè)由東南向西北流的支流匯入流域1的主溝道(圖中黃色虛線溝道)。而待重塑區(qū)南部的Ⅰ級(jí)溝道由北向南匯合成1個(gè)Ⅱ級(jí)溝道，并匯入到流域2主溝道(圖中的藍(lán)綠色虛線溝道)。待重塑區(qū)所覆蓋的都是Ⅰ級(jí)溝道，并與流域1和流域2相聯(lián)通。由于位置鄰近，且這2個(gè)自然流域與待重塑區(qū)在氣候、土壤、植被和巖性等方面都具有很大的相似性，依據(jù)相似性原則，選擇2個(gè)流域內(nèi)含Ⅰ級(jí)溝道的流域作為參照流域(圖6)。

圖6 分級(jí)水文網(wǎng)與Ⅰ級(jí)子流域劃分示意圖Fig.6 Schematic diagram of classification of hierarchical hydrological network and level Ⅰ sub-watershed

圖6為初步確定的29個(gè)Ⅰ級(jí)子流域區(qū)域。其中子流域1、3和4與待重塑區(qū)域相交，因此剔除不作為參照流域。將除上述3個(gè)子流域外的26個(gè)Ⅰ級(jí)流域作為初步選定的參照流域，26個(gè)子流域的6項(xiàng)參數(shù)計(jì)算結(jié)果見表1。

這26個(gè)子流域中，有些子流域與其他流域形態(tài)特征差別很大，對(duì)數(shù)據(jù)統(tǒng)計(jì)造成干擾；為此通過(guò)箱型圖方法[6]對(duì)子流域6項(xiàng)特征參數(shù)進(jìn)行統(tǒng)計(jì)，剔除“異?！绷饔?，最終篩選出21個(gè)子流域作為參照流域。在此基礎(chǔ)上，將這21個(gè)子流域區(qū)域所屬的2條主流域(流域1和流域2)所轄區(qū)域作為參照區(qū)域(圖5中的粉紅色框線區(qū)域)進(jìn)行后續(xù)局部地表形態(tài)特征的計(jì)算和分析。

3.2 流域地貌形態(tài)特征及溝道形態(tài)特征

采用SPSS軟件對(duì)表1中21個(gè)參照流域的6項(xiàng)流域地貌特征參數(shù)進(jìn)行S-W檢驗(yàn)，其顯著性>0.05。表明參數(shù)分布都較為集中，符合正態(tài)分布，可作為地形重塑參考目標(biāo)。21個(gè)流域的參數(shù)均值為：平均子流域面積0.25 km2；平均主溝道長(zhǎng)度0.54 km；平均溝道密度2.29；平均流域圓度0.27；流域高程差32.33 m；平均流域坡度4.61 °。

表1 初步選定的26個(gè)Ⅰ級(jí)子流域地貌形態(tài)特征Tab.1 Morphological characteristics of the 26 level Ⅰ sub-watersheds initially selected

進(jìn)一步地，對(duì)參照流域內(nèi)的23條Ⅰ級(jí)溝道(6、24號(hào)子流域內(nèi)分別有2條 Ⅰ 級(jí)溝道)溝道形態(tài)特征參數(shù)進(jìn)行測(cè)量(表2)。

根據(jù)河流分類標(biāo)準(zhǔn)[15]，當(dāng)蜿蜒度S=1.0時(shí)，為順直河段；S< 1.2時(shí)為低度蜿蜒，S=[1.2，1.4]為中度蜿蜒，S>1.4為高度蜿蜒。從表2中得出目標(biāo)流域內(nèi)溝道的蜿蜒度范圍為[1.030，1.316]，其中57%的溝道蜿蜒度<1.2；35%的溝道蜿蜒度在[1.2，1.3]之間，說(shuō)明大部分溝道屬于中低度蜿蜒溝道。參照流域的溝道比降在[0.001 0，0.045 1]之間，其中52%的溝道比降<0.02；30%的溝道比降在[0.020 0，0.030 0]之間；除11號(hào)溝道，溝道比降都<0.040 0。

表2 Ⅰ級(jí)溝道形態(tài)參數(shù)Tab.2 Level Ⅰ gully shape parameters

對(duì)這3項(xiàng)參數(shù)數(shù)據(jù)進(jìn)行S-W檢驗(yàn)，3項(xiàng)參數(shù)顯著性> 0.05，表明數(shù)據(jù)分布集中，服從正態(tài)分布。此外，對(duì)3項(xiàng)參數(shù)進(jìn)行箱線圖分析，發(fā)現(xiàn)無(wú)異常值。因此，可將此3項(xiàng)參數(shù)的均值作為近自然地形重塑設(shè)計(jì)的參考目標(biāo)。

3.3 局部地表形態(tài)特征

3.3.1 坡度 利用ArcGIS對(duì)參照區(qū)域進(jìn)行坡度分析，根據(jù)我國(guó)水力侵蝕強(qiáng)度分級(jí)指標(biāo)中面蝕分級(jí)標(biāo)準(zhǔn)，將坡度分為6個(gè)等級(jí)生成面蝕分級(jí)標(biāo)準(zhǔn)下的坡度分級(jí)圖(圖7)。從圖中可知，自然地形中大部分區(qū)域?yàn)椤? °的不易發(fā)生面蝕的平地和微坡面，面積占72.74%；而>15 °的陡坡區(qū)域僅占1.42%。

3.3.2 坡向 坡向利用ArcGIS對(duì)參照區(qū)域進(jìn)行坡向分析，生成坡向分級(jí)占比圖(圖7b)。根據(jù)圖中的占比數(shù)據(jù)統(tǒng)計(jì)可知，自然地形的主要坡向?yàn)楸?、東、西、東北、西北，共占區(qū)域面積的70.41%。若按照陰坡、陽(yáng)坡、半陰坡、半陽(yáng)坡對(duì)坡向重分類，比例分別為30.18%、19.34%、24.45%和26.04%；陰坡面積明顯大于陽(yáng)坡。

圖7 局部地表形態(tài)特征參數(shù)分布圖(括號(hào)中的百分?jǐn)?shù)為該等級(jí)區(qū)域占整體面積的比例)Fig.7 Distribution map of local surface morphological characteristic parameters (the percentage in parentheses after the level is the percentage of the area of the level in the overall area)

研究區(qū)的主風(fēng)向?yàn)槲鞅憋L(fēng)，而參照區(qū)域中偏北、偏西坡向所占比例較大。從圖中可看出許多子溝道的分水嶺呈東北- 西南走向，因而形成許多垂直于西北風(fēng)向的坡地。這種垂直于主力風(fēng)向的山脊地形對(duì)防風(fēng)固沙有一定作用，這對(duì)于地形重塑有一定參考意義。

3.3.3 局部地形起伏度 局部地形起伏度指局部窗口內(nèi)地形高程最大值與最小值之差。為確定地形起伏度最佳分析窗口，依次計(jì)算網(wǎng)格2×2、3×3……10×10下的平均地形起伏度，通過(guò)分析地形起伏度與網(wǎng)格面積間的關(guān)系曲線，采用均值變點(diǎn)分析法確定7×7網(wǎng)格為最佳分析窗口，在該網(wǎng)格窗口下的地形起伏度情況見圖7c。從圖中可知，比例最大的為(1,2] m起伏度的區(qū)域(35.01%)；≤3 m的地形區(qū)域占72.1%；而>5 m的區(qū)域僅占7%。說(shuō)明整體地形起伏度較小，地勢(shì)平坦。

3.3.4 局部地表粗糙度 局部地表粗糙度可反映地形在水平方向上的復(fù)雜度，是衡量地表侵蝕程度的重要量化指標(biāo)之一。在ArcGIS中計(jì)算得到分級(jí)地表粗糙度(圖7d)。從圖中可見，自然地形96.63%區(qū)域的地表粗糙度在[1.000,1.019]之間，粗糙度較低，表明侵蝕程度較小。

4 結(jié)論

1)流域侵蝕作用依然是影響西北干旱荒漠區(qū)地表形態(tài)的最重要因素，因此在參照區(qū)域選擇上，提出以流域?yàn)檫x取單位，依據(jù)流域的鄰近性、相似性，并兼顧氣候、土壤、巖性等其他自然因素的策略。

2)從流域地貌形態(tài)、溝道形態(tài)和局部地表形態(tài)3方面出發(fā)，選取14項(xiàng)參數(shù)作為刻畫近自然地形特征的因子。這3方面的參數(shù)綜合近自然地形重塑對(duì)不同尺度(流域尺度、局部地形小尺度)以及不同地形構(gòu)建單元(流域、溝道、坡面)的參數(shù)需求，以滿足近自然地形重塑的設(shè)計(jì)之需。

3)基于研究區(qū)流域分析，選取21個(gè)I級(jí)子流域作為參照流域。統(tǒng)計(jì)得到這21個(gè)參照子流域的流域地貌特征參數(shù)值及溝道特征參數(shù)值的分布比較集中，顯著性> 0.05，表明這些參數(shù)值有代表性，可作為地形重塑的參照目標(biāo)。