基于自然恢復的采煤沉陷區(qū)植物群落演替模擬系統(tǒng)開發(fā)

米家鑫，張紹良，楊永均，侯湖平，丁忠義

(1.中國礦業(yè)大學 公共管理學院，江蘇 徐州 221116；2.中國礦業(yè)大學 環(huán)境與測繪學院，江蘇 徐州 221116；3.中國礦業(yè)大學 教育部礦區(qū)生態(tài)修復重點實驗室,江蘇 徐州 221116)

煤炭開發(fā)戰(zhàn)略西移后，干旱半干旱區(qū)高強度地下開采引起了大面積的地表沉陷、水體破壞，嚴重脅迫著脆弱的生態(tài)系統(tǒng)，而大范圍的人工修復往往成本高、效率低、效果不明顯，因此一些學者提出了自然恢復的策略。盡管一些研究發(fā)現(xiàn)采煤沉陷區(qū)的植被的確具有自我恢復的能力，但是植被自然恢復過程時間長、效果和格局難以預測，這種不確定性極大地增加了采煤沉陷區(qū)植被自然恢復的風險，因此需要開發(fā)自然恢復效果評價技術(shù)，模擬礦區(qū)受損植物群落的自然演替過程成為了重要途徑。

植物群落的演替過程十分復雜，涉及到了多個生態(tài)水文過程，根據(jù)模擬模型的原理主要分為經(jīng)驗模型與機理模型。經(jīng)驗模型主要依賴于植物群落類型、與其生長演替相關(guān)要素間的相關(guān)關(guān)系，在將其擬合為定量的經(jīng)驗函數(shù)后進行模擬。例如物種分布模型、Holdridge生命地帶模型、BP神經(jīng)網(wǎng)絡(luò)模型等。但由于經(jīng)驗模型具有較強的區(qū)域性特征導致適用范圍較小，在植物群落演替模擬中研究較少。另一方面，機理模型則通過模擬植物群落演替過程中發(fā)生的生態(tài)水文變化過程，最終輸出植物群落的模擬結(jié)果。其中應(yīng)用較為廣泛的包括林窗模型和全球植被動態(tài)模型。林窗模型(Forest Gap Model)最早由BOTKIN與1972年提出，主要通過描述樹木的生長、死亡過程來模擬森林的演替過程，最終模型將會輸出模擬過程中不同樹木的數(shù)量、樹木的高度、森林的密度等參數(shù)的動態(tài)變化。全球植被動態(tài)模型(DGVM)，則主要從物質(zhì)交換和能量流動的角度模擬植物群落的生長演替過程，通過將氣象數(shù)據(jù)、CO濃度數(shù)據(jù)和土壤數(shù)據(jù)等輸入植被動力模塊以及生物地球物理-生物地球化學模塊后，輸出植物群落的類型和狀態(tài)?？傮w來看，已有的植物群落演替模擬系統(tǒng)中主要依靠降雨、光照等氣候數(shù)據(jù)模擬大尺度區(qū)域內(nèi)的物質(zhì)、能量流動來模擬植物群落的自然演替過程，而地形這類在小范圍區(qū)域內(nèi)存在空間差異的因素往往被忽視，沒有考慮地形差異對各類生態(tài)水文過程的影響。因此這些模擬系統(tǒng)大多難以反映模擬單元內(nèi)部植物群落在自然恢復過程中的空間分布及實現(xiàn)其可視化，更無法反映因采煤沉陷區(qū)地形變化而改變的植物群落演替過程及植被格局。

采煤沉陷區(qū)復雜地形條件下的水文生態(tài)過程是植被自然恢復模擬的難點。沉陷裂縫出現(xiàn)后改變了原有地形中坡度、坡向、高程等因素，從而影響著地表徑流、水分蒸發(fā)、水土流失等多個水文過程，進一步影響土壤養(yǎng)分積累、種子萌發(fā)、群落發(fā)展等過程。因此小尺度下生態(tài)水文過程的變化是影響采煤沉陷區(qū)植被自然恢復的主要因素。而元胞自動機因其能夠通過簡單的局部轉(zhuǎn)換規(guī)則來模擬十分復雜的空間結(jié)構(gòu)，常被用于生態(tài)水文過程模擬中，例如坡面水蝕模型、河道匯流模擬模型、城市內(nèi)澇模型等。一些學者也將其用于植被格局的模擬中，例如有學者基于元胞自動機模型模擬了高寒草甸退化過程中“禿斑”現(xiàn)象中出現(xiàn)的連通效應(yīng)；有的結(jié)合邏輯斯第增長模型模擬了不同資源條件形成的植被空間格局；還有的構(gòu)建了林火蔓延快速模型，模擬了火災發(fā)生時樹木的燃燒傳遞現(xiàn)象。然而，已有的基于元胞自動機的植被格局模擬研究大多使用植物群落與單一因素或單個生態(tài)水文過程之間的相關(guān)關(guān)系，作為演替的限制或驅(qū)動條件，卻沒有考慮復雜地形條件下地表徑流、水分蒸發(fā)、水土流失等多個水文過程相互作用關(guān)系以及這一過程中植被的變化，從而難以完整的模擬采煤沉陷區(qū)的植被自然恢復過程及最終形成的植被格局。

筆者基于植被自然恢復規(guī)律和元胞自動機模型開發(fā)了植物群落演替模擬系統(tǒng)，并以山西省大同市云岡區(qū)的煤峪口礦采煤沉陷區(qū)為例，模擬了植被自然恢復下的群落演替過程，以及演替形成的植被格局，可以為礦區(qū)植被自然恢復策略的實施提供案例參考和決策工具。

1 研究區(qū)概況

大同市云岡區(qū)位于山西省北部，地理位置坐標為東經(jīng)112.88°～113.19°，北緯40.17°～40.42°，如圖1所示。全區(qū)地處黃土高原，屬大陸性季風氣候，全年平均氣溫6.4 ℃，年平均降水量394.6 mm，主要集中在6—9月份。筆者選擇的沉陷裂縫區(qū)域位于大同煤業(yè)煤峪口礦沉陷裂縫區(qū)內(nèi)，總面積約為0.13 km，區(qū)域內(nèi)主要植被包括沙棘、狗尾草等灌木和草本植物。煤峪口礦年生產(chǎn)能力250萬t，自2012年至今，煤峪口煤礦主要對14-2號煤層進行開采，煤層平均厚度為2.11 m，采深在310～376 m，采煤方法以長壁大冒頂機械化綜合采煤為主。根據(jù)概率積分法的沉陷預計模型，煤峪口煤礦在2006—2018年間出現(xiàn)的下沉影響面積達到8.65 km(下沉大于10 mm)，沉陷達到2 000 mm以上的區(qū)域占到8.21%。

圖1 研究區(qū)位置Fig.1 Location of study area

2 植物群落演替模擬系統(tǒng)設(shè)計

2.1 模擬系統(tǒng)總體框架

基于干旱半干旱區(qū)氣候、水文、植被特征，本研究中的植物群落演替模擬系統(tǒng)由地表形變模塊、降雨徑流模塊和植被生長模塊組成，如圖2所示。在干旱半干旱礦區(qū)，當?shù)乇硇巫儼l(fā)生后，首先以地表沉陷和地表裂縫為代表的地表形變通過根系拉伸等直接引起植被的死亡。此外，地表形變也導致了區(qū)域內(nèi)高程、坡度和坡向的變化，并將對后續(xù)的水文過程產(chǎn)生影響。在所有生態(tài)過程和要素中，降雨既是地表植被水分的主要來源同時也是植物群落演替的主要限制因素。當降雨到達地表后，土壤入滲隨即發(fā)生將地表水轉(zhuǎn)化為土壤水；當單位時間的降雨量大于土壤入滲速率時一部分降雨轉(zhuǎn)化為地表徑流，并在重力作用下從高向低流出；高程決定了徑流的流向，坡度則影響著入滲速率和徑流速率，坡向則影響了土壤的蒸發(fā)速率。降雨結(jié)束后土壤蒸發(fā)隨即發(fā)生，土壤水分逐漸降低至田間持水率。當植被立地條件滿足不同類型的植被生長所必須的條件后，植被種子便有概率萌發(fā)。而在生長季結(jié)束之后，草本植被和一部分灌木植被死亡，一部分灌木植被和喬木植被則凋落，植被殘枝和凋落物進入土壤轉(zhuǎn)換為土壤養(yǎng)分；凋落的種子也將進入土壤種子庫，從而為次年植被提供養(yǎng)分和種子。這些過程將會不斷循環(huán)，最終促成了區(qū)域植物群落格局的形成。

圖2 植物群落演替模擬系統(tǒng)框架Fig.2 Framework of simulation model for spontaneoussuccession of plant communities

本次研究構(gòu)建的模擬系統(tǒng)主要基于元胞自動機(Cellular Automata，CA)模型，它是一種時間和空間上離散的動力學系統(tǒng)，主要由元胞空間、元胞鄰居、元胞狀態(tài)及轉(zhuǎn)化規(guī)則組成。

2.2 模擬系統(tǒng)元胞單元設(shè)計

在本研究中，為了兼容地形數(shù)據(jù)特選擇四邊形網(wǎng)格的劃分方式作為元胞空間(圖3(a))?？紤]到植物群落的劃分尺度，將每個元胞大小設(shè)定為現(xiàn)實世界中1 m的方塊，下文中植被覆蓋度、等狀態(tài)變量均以1 m大小元胞為最小空間單位。

圖3 四邊形網(wǎng)格元胞與摩爾型鄰居模型Fig.3 Quadrilateral grid cell and Mooreneighbor model

2.3 模擬系統(tǒng)元胞鄰居選擇

元胞鄰居定了元胞之間的傳遞關(guān)系及影響范圍。考慮到現(xiàn)實情況下地表徑流和植被生長過程中相鄰斑塊間的作用方式，本次研究選擇摩爾型作為此模型的元胞鄰居，如圖3(b)所示。每個元胞受到與其相鄰的8個鄰居元胞的作用，在不同的模塊中遵循不同的規(guī)則變化。

2.4 模擬系統(tǒng)元胞狀態(tài)參數(shù)選取

元胞狀態(tài)表示元胞在模型中的非空屬性集合。每個元胞都有特定的狀態(tài)值，并隨著模型的演化時刻發(fā)生改變。本次研究中每個元胞的狀態(tài)主要包括地形參數(shù)、入滲參數(shù)、徑流參數(shù)、蒸發(fā)參數(shù)、土壤參數(shù)、植被參數(shù)，以及屬于模擬系統(tǒng)的全局參數(shù)和降雨參數(shù)，見表1。

表1 模擬系統(tǒng)的元胞狀態(tài)參數(shù)

2.5 模擬系統(tǒng)元胞規(guī)則設(shè)計

轉(zhuǎn)換規(guī)則是元胞自動機的核心，其實質(zhì)是根據(jù)元胞當前時刻的狀態(tài)和鄰居狀態(tài)來決定下一個時刻該元胞的狀態(tài)的動力學函數(shù)。在本次研究中，主要包含降雨規(guī)則、入滲規(guī)則、徑流規(guī)則、蒸發(fā)規(guī)則、養(yǎng)分積累及流失規(guī)則、種子庫積累規(guī)則、植被生長規(guī)則。

2.5.1 根系損傷規(guī)則

地表形變發(fā)生時不僅引起了地表高程、坡度、坡向的變化，同時也導致了部分植被的死亡。考慮到不同類型植被的生長方式及對外界擾動的恢復力差異，設(shè)置草本、灌木和喬木植被對沉陷的忍耐閾值分別為1.5，1.0和0.5 m，即當沉陷深度大于1.5 m時，區(qū)域內(nèi)的草本植被即死亡，而灌木和喬木當沉陷超過1.0 m和0.5 m時死亡。在模擬模型中，將地表形變設(shè)置為模擬的第1天發(fā)生，死亡的植被將作為養(yǎng)分被土壤吸收，養(yǎng)分吸收規(guī)律與植被自然死亡時一致。

2.5.2 降雨規(guī)則

在真實的降雨場景中，單位時間落入地表的雨水呈現(xiàn)隨機分布規(guī)律。將雨滴橫截面作為單位面積并投影到其落入的地表，將秒作為單位時間并記錄每一秒內(nèi)落入地表的雨水數(shù)量，每一秒內(nèi)落入雨水的地表是隨機的，而單位時間內(nèi)含有降雨落入地表的面積與地表總面積的比值即降雨落入概率。考慮到模型的整體運行速度，將運行的步長即單位時間設(shè)定為分鐘，空間單位為1 m的元胞。在本次研究中，以云岡區(qū)過去50 a的降雨情況作為參考，分別確定不同季節(jié)的降雨量。由于云岡區(qū)降水70%集中在春夏季，且模擬系統(tǒng)僅將春夏季作為植被的生長期，因此僅將在春夏季時考慮降雨情況。其中設(shè)置春季為每年度的第59～151天，夏季為第151～243天，其余時間為秋冬季。具體參數(shù)見表2。

表2 生長季降雨參數(shù)設(shè)置

為了反映區(qū)域真實的降雨情況，通過區(qū)域的季度降雨總量、降雨天數(shù)、降雨時間、降雨強度、雨水落入概率以計算降雨強度，其具體公式為

=

(1)

(2)

式中，為每一步長內(nèi)元胞獲得的雨量;為雨強;為雨滴落入概率;為累計降雨量；為降雨天數(shù)；為降雨時間。

2.5.3 土壤入滲規(guī)則

土壤入滲是指地表水在重力勢、基質(zhì)勢等作用下運移、存儲變?yōu)橥寥浪倪^程。對于植被覆蓋的土壤來說，由于植被覆蓋能有效減緩地表水的流動，從而大大增加其穩(wěn)定入滲速率；坡度的增加則會減少來自重力的分量，從而引起入滲速率的減小，如圖4所示。而對于無植被覆蓋的土壤來說，由于其地表更容易受到降雨徑流的侵蝕從而引起土壤質(zhì)地的惡化，導致其土壤入滲速率的降低。因此針對裸地土壤增加了一個養(yǎng)分的調(diào)節(jié)參數(shù)，當土壤養(yǎng)分低于土壤初始條件時，其土壤入滲率將逐漸下降，以模擬現(xiàn)實情況下土壤侵蝕增加并逐漸沙化的現(xiàn)象。

圖4 降水在坡面的重力加速度分量Fig.4 Gravity acceleration component of precipitation on slope

=(191+0326e-0215)cos

(3)

(4)

式中，，分別為植被覆蓋和無植被覆蓋區(qū)域的土壤入滲速率；為入滲時間；為土壤養(yǎng)分含量；為土壤初始養(yǎng)分含量。

為了保證模型運行的連貫性，在此模型中僅考慮表土土壤水分的變化。入滲結(jié)束后，表土土壤的水分增加量即等于累計入滲量，而表土為地表20 cm以內(nèi)的土壤。當土壤水分達到土壤的飽和含水量時，土壤入滲過程將會暫停，多余的地表水會以徑流的形式流向周邊區(qū)域。通過對研究區(qū)域土壤的采樣化驗，測得植被覆蓋和無植被覆蓋土壤的飽和含水率分別0.35與0.30。為了便于降雨量、入滲量、含水量等不同量綱變量之間的轉(zhuǎn)化，筆者定義的含水量均統(tǒng)一為單位體積內(nèi)土壤水分的比例。由于降雨形成的地表水以及地表土壤的投影面積均為1 m，所以土壤含水量等于累計入滲量的深度與表土深度(20 cm)的比值。

2.5.4 徑流規(guī)則

當?shù)乇硭笥谕寥廊霛B速率時開始出現(xiàn)徑流。對于單個元胞來說，地表水輸入來源包括降雨和來自其他元胞的徑流，輸出方式則包括土壤入滲及流出至其他元胞的徑流。本文中降雨量的單位為mm，其含義為單位時間內(nèi)落入單位面積內(nèi)水量的深度。但對于地表徑流，由于坡度的存在使得地表徑流實際通過的面積將大于元胞的投影面積，如圖5所示。因此，從降雨中獲取的徑流量需要進行轉(zhuǎn)化，根據(jù)勾股定理可得，坡度為的坡面實際單位面積為1/cos，則徑流公式為

′=

(5)

=′cos

(6)

式中，′為單位元胞通過降雨獲得的理論徑流量；為考慮地表坡度后單位元胞上的實際徑流量。

圖5 不同坡度元胞的徑流狀態(tài)Fig.5 Runoff status of cells with different slopes

當徑流在一個元胞(,)內(nèi)形成后，如果元胞(,)并非與其相鄰的8個鄰居元胞中的高程最小元胞，那么徑流將在重力的作用下從該元胞流向相對于此元胞的最低點元胞(+1,+1)中。由于重力方向是固定的，所有每個元胞形成徑流的流出方向只有1個，徑流將始終流入高程相對較低的元胞，如圖6所示。

圖6 非相對最低點元胞徑流的流出方向Fig.6 Outflow direction of the runoff from non relativelowest cellular

當元胞(,)的高程小于周圍所有鄰居元胞時，此元胞則成為徑流的匯聚點，鄰居元胞形成的徑流都將流入元胞(,)中，且元胞(,)形成的徑流則不會流出。但是，當元胞(,)內(nèi)匯集的徑流超出其所能匯集的最大值時，元胞(,)的鄰居元胞則將停止流出。元胞(,)的鄰居元胞內(nèi)的徑流由于無法再流出，這些鄰居元胞也將逐漸匯集徑流，從而出現(xiàn)徑流由低向高增加的“溢出”現(xiàn)象，如圖7所示(圖中為匯聚點元胞的高程；為匯集點元胞周邊的最高點元胞的高程；Drop為匯聚點元胞與最高點元胞之間的落差)。

圖7 相對最低點元胞徑流的流出方向Fig.7 Outflow direction of the runoff from relative lowest cellular

在徑流過程中，元胞(,)的流出量與其徑流量以及徑流速率有關(guān)，而流入量則等于其鄰居元胞的流向元胞(,)的流出量總和。當元胞內(nèi)的徑流量小于元胞的徑流速率時，流出量等于此時的徑流量當元胞內(nèi)的徑流量大于元胞的徑流速率時，流出量則等于元胞的徑流速率。其中，元胞(,)的徑流速率等于單位時間內(nèi)通過元胞的最大徑流量，這一速率通過曼寧公式計算得到。曼寧公式是一個估測液體在開放管道(即明渠流)或非滿管流中平均速度的經(jīng)驗公式，具體公式為

(7)

式中,為轉(zhuǎn)換系數(shù)，國際標準制為1；為曼寧粗糙系數(shù)，通過查詢曼寧系數(shù)表確定；為水力半徑，計算方式是過水截面積與濕周長的比值；為水力坡線或是線性揚程損失的斜率。

在本次研究中，將每個元胞作為徑流流動的基本空間單元，則為通過元胞的的水深，為了便于降雨量、徑流量之間的水量交互與計算，對徑流量與降雨量均統(tǒng)一使用毫米作為單位，因此在計算最大徑流速率的公式中將徑流量直接用于水深計算；經(jīng)過查詢曼妮系數(shù)表后，規(guī)定裸土的粗糙系數(shù)為0.02，植被覆蓋土壤的粗糙系數(shù)為0.095；斜率則等于坡度的正切值tan。

2.5.5 蒸發(fā)規(guī)則

當土壤含水率高于土壤的田間持水率時則土壤水分開始蒸發(fā)，這一過程主要集中在無降雨期間。土壤的田間持水率主要受到覆蓋類型的影響，存在植被覆蓋的土壤其保水持水能力較強，而裸土則相對較弱。本研究中，測得植被覆蓋和無植被覆蓋土壤的田間持水率分別為0.22與0.15。蒸發(fā)速率在不同季節(jié)內(nèi)有地形位置的不同而有所差異。由于不同坡向的差異，陽坡區(qū)域的土壤由于光照時間較長而溫度較高，土壤的蒸發(fā)速率也較高。地表的植被覆蓋則會一定程度的降低地表土壤的溫度，從而降低土壤的蒸發(fā)速率。參考云岡區(qū)過去50 a的氣象數(shù)據(jù)，并將蒸發(fā)速率轉(zhuǎn)化為日蒸發(fā)量后，確定研究區(qū)不同季節(jié)、不同坡位及不同植被覆蓋土壤的蒸發(fā)速率見表3。

表3 不同坡向的角度范圍及蒸發(fā)速率

2.5.6 土壤養(yǎng)分積累及流失規(guī)則

土壤養(yǎng)分的積累主要來源于植被凋落物和殘枝分解。當植被在秋冬季凋落枝葉時，其凋落物將被土壤中的微生物分解并轉(zhuǎn)換為土壤養(yǎng)分；當植被死亡時，其殘留的根莖也將被土壤吸收轉(zhuǎn)化土壤養(yǎng)分其中，灌木和喬木植被包含植被死亡和凋落過程，而草本植被只通過死亡過程轉(zhuǎn)化養(yǎng)分。此外，由于灌木植被和喬木植被的冠幅較大，其凋落和死亡后的枝葉通常也將落入周邊區(qū)域，因此喬木和灌木植被的鄰居元胞也將在其凋落和死亡時獲得養(yǎng)分。本文使用土壤有機質(zhì)含量反映土壤養(yǎng)分情況，基于研究區(qū)土壤養(yǎng)分水平情況發(fā)育規(guī)律，規(guī)定不同植被凋落或死亡后土壤養(yǎng)分的增加情況見表4。為了避免局部區(qū)域養(yǎng)分過高不符合現(xiàn)實情況，同時設(shè)置土壤養(yǎng)分的最大值為40 g/kg。

表4 植被死亡或凋落時土壤養(yǎng)分增加量

另一方面，徑流對土壤的沖刷裹挾作用導致了土壤養(yǎng)分的改變。徑流通過土壤時，一部分土壤以及土壤中的有機質(zhì)、氮、磷、鉀等養(yǎng)分會在徑流的裹挾作用下融入到徑流中，從而造成土壤養(yǎng)分的流失。徑流的強度越大，其造成的養(yǎng)分流失則越明顯；降雨的強度越大，土壤養(yǎng)分的速率越快。在這一過程中，地表植被則會起到減弱徑流、穩(wěn)固水土的作用，因此植被覆蓋土壤的養(yǎng)分流失速率將明顯低于裸土的流失速率。以往的研究針對徑流引起的土壤養(yǎng)分流失積累大量預測模型，通過參考之前的研究，確定了本研究中不同類型土壤的養(yǎng)分流失速率(表5)。此外，分別設(shè)置5，30 g/kg作為土壤養(yǎng)分的最低值和最高值，為避免一些區(qū)域土壤養(yǎng)分持續(xù)流出或流入而造成養(yǎng)分過高或過低的情況。

表5 土壤養(yǎng)分流失速率

2.5.7 種子庫積累規(guī)則

植被在生長過程中都會產(chǎn)生一定數(shù)量的種子，這些種子大多會凋落在植被周圍并進入土壤種子庫中；相似的，周邊區(qū)域植被的種子也可能掉入其相鄰區(qū)域，從而滿足植被生長所需的土壤種子庫條件。此外，風力運輸?shù)扰既皇录瑯訒槟骋粎^(qū)域的土壤提供種子。在本次研究中，規(guī)定植被死亡后其種子數(shù)量增加1；而一天內(nèi)單個元胞因為偶然事件獲得種子的概率為0.000 4，偶然事件發(fā)生后其種子數(shù)量增加1。對于來自相鄰植被獲得種子，則統(tǒng)一使用邏輯判斷，即當周邊存在相同類型的植被時，則滿足植被生長的種子庫條件。

2.5.8 植被生長規(guī)則

不同植被生長所需的土壤水分、養(yǎng)分、種子庫以及種子萌發(fā)概率存在一定差異。植物群落的發(fā)展基本遵循草本—灌木—喬木的演替順序，且灌木植被所需的土壤養(yǎng)分條件一般高于草本植被，但在水分條件上無明顯差異，而喬木所需的土壤養(yǎng)分和水分條件均較高。此外為了模擬植物群落的“聚集”現(xiàn)象，同時考慮鄰居元胞植被覆蓋的影響?；谝陨弦?guī)則，在參考以往黃土丘陵區(qū)植被自然恢復的研究案例基礎(chǔ)上，確定了不同類型植被所需的土壤養(yǎng)分、水分和種子庫條件見表6。

表6 不同類型植被種子萌發(fā)條件

而對于不同類型植被的萌發(fā)概率來說，在一定范圍內(nèi)水分條件和養(yǎng)分條件越充足，其種子萌發(fā)的概率越高；鄰居元胞中存在灌木和喬木的數(shù)量越多，種子萌發(fā)的概率同樣越高。因此通過土壤實際的水分和養(yǎng)分與種子萌發(fā)所需的水分和養(yǎng)分的比值作為調(diào)節(jié)參數(shù)，以反映不同水分和養(yǎng)分的土壤在種子萌發(fā)條件上的差異。

2.6 模擬區(qū)域環(huán)境條件設(shè)置

此次研究中使用研究區(qū)已獲取的無人機影像作為模擬區(qū)域的起始地形條件，并將其植被格局作為模擬結(jié)果的驗證數(shù)據(jù)。研究區(qū)位于云岡礦區(qū)內(nèi)煤峪口礦開采范圍內(nèi)，開采時間在2006—2008年。無人機影像采集時間為2019-08-22 T 09:00—12:00，距沉陷裂縫形成時間約為11 a，影像獲取時區(qū)域內(nèi)植被生長茂盛。無人機飛行高度為100 m，云臺角度-90°，旁向重疊率65%，航向重疊率70%。獲取無人機影像后，使用Agisoft PhotoScan軟件進行影像預處理，并通過數(shù)字表面模型去除植被高度后獲得數(shù)字地形模型。然后將數(shù)字地形模型影像作為模擬區(qū)域的高程數(shù)據(jù)，并在ArcMap中通過柵格分析功能轉(zhuǎn)化為區(qū)域的坡度及坡向數(shù)據(jù)，如圖8所示。模擬區(qū)域的初始植被條件設(shè)置為無植被，而驗證數(shù)據(jù)則根據(jù)通過數(shù)字正射影像確定模擬區(qū)域?qū)嶋H的植被類型，如圖8所示。為了與植物群落演替模擬系統(tǒng)中的元胞空間匹配，在ArcMap中將高程、坡度、坡向及植被類型數(shù)據(jù)重采樣為1 m×1 m的柵格。考慮到模擬區(qū)域內(nèi)存在少量人為建造的道路，在ArcMap中通過矢量化構(gòu)建了掩膜圖層并在模擬結(jié)束后添加至模擬結(jié)果中。

圖8 沉陷裂縫區(qū)地形條件Fig.8 Terrains of subsidence area

2.7 模擬結(jié)果精度評價方法

由于植物群落自然演替的隨機性，比較不同景觀尺度下的植被格局進行精度評價較為合理。首先將模擬結(jié)果與驗證數(shù)據(jù)劃分為若干個大小的斑塊(圖9)，依次比較各個斑塊內(nèi)不同類型的植被覆蓋度。當模擬結(jié)果與驗證數(shù)據(jù)中對應(yīng)斑塊內(nèi)植被覆蓋度、草本覆蓋度及灌木覆蓋度相似，則說明此模擬系統(tǒng)能較好的模擬現(xiàn)實情況中的群落演替?？紤]到不同景觀尺度下植物群落格局的差異，本次研究中選擇了10×10，20×20,50×50,80×80,100×100 m的5種斑塊大小(P10，P20，P50，P80，P100)。

圖9 植被覆蓋類型及群落斑塊劃分情況Fig.9 Types and patches of plant communities

在本研究中，采用歸一化均方誤差反映模型的模擬精度。歸一化均方誤差是一種反映2張圖像之間相似度的指數(shù)，其值域為0到正無窮，接近0則說明實際值與預測值之間差異越小，模型的模擬精度越高，具體公式為

(8)

當模擬結(jié)果中的植被覆蓋度接近驗證數(shù)據(jù)的植被覆蓋度時，即停止模擬并將此時的植物群落格局記錄并與驗證數(shù)據(jù)進行比較。在實際情況中，植被覆蓋度為50.79%，草本覆蓋度為43.27%，灌木覆蓋度為7.51%，而喬木覆蓋度僅為0.01%。因此在進行模擬時僅將草本和灌木植被作為模擬結(jié)束和精度評價的標準，當模擬結(jié)果的草本植被覆蓋度為35%～55%。且灌木覆蓋度為7%～8%時即停止模擬，輸出此時的植被格局并分別計算不同尺度下的。為了減少單次模擬帶來的隨機性，設(shè)置模擬次數(shù)為50次，取每次模擬結(jié)果與驗證數(shù)據(jù)的NMSE平均值為最終值。

3 結(jié)果與討論

3.1 模擬結(jié)果與系統(tǒng)精度

植物群落的模擬結(jié)果與實際分布對比如圖10所示，其中圖10(a)～(c)為隨機選擇的植被自然恢復模擬結(jié)果，圖10(d)為實際的植被分布情況。從模擬結(jié)果與實際結(jié)果的植物群落整體來看，灌木植物群落均主要集中在沉陷裂縫區(qū)的低矮區(qū)域，且都表現(xiàn)出較為明顯的群落聚集現(xiàn)象。但在模擬結(jié)果中，草本植被在整個區(qū)域內(nèi)較為均勻分布，無明顯的無植被覆蓋區(qū)域。而實際情況中存在較為明顯的無植被覆蓋區(qū)域。實際調(diào)查發(fā)現(xiàn)，這些區(qū)域均為陡峭的邊坡，經(jīng)過長年的雨水沖刷以及成為裸露的巖石。這一差異主要由于此次模擬的初始土壤被設(shè)置為適宜草本植被的條件，但在實際情況下這些區(qū)域可能在自然恢復開始時便不適宜植被的生長，從而導致了模擬結(jié)果與實際情況在這些局部區(qū)域的差異。

圖10 驗證區(qū)域植物群落格局模擬結(jié)果Fig.10 Simulation results of vegetation pattern in stuyd area

從模擬結(jié)果與實際格局的來看，隨著斑塊尺度的增加模擬結(jié)果與實際情況的差異度隨之降低，表明模擬精度也隨之增加。總體模擬精度如圖11所示，其中，灌木植被的空間分布格局差異度最高，5種斑塊尺度下的在0.06～0.59；草本植被的總體較小，在0.05～0.14；而從總體的植被覆蓋度來看，植被格局的模擬結(jié)果與實際情況差異最小，其在0.04～0.12。這主要來自于不同的植被覆蓋度基數(shù)的差異，其中灌木覆蓋度基數(shù)最小，因此在劃分的板塊中較小的差異都將被放大，尤其考慮到植被的自然演替本身便是充滿隨機性的過程，使得灌木植物群落格局的模擬在較小尺度精度較低。但從整體效果來看，目前設(shè)計的模擬系統(tǒng)能夠還原現(xiàn)實世界中的植物群落格局。

圖11 不同尺度下模擬結(jié)果和實際情況的差異度Fig.11 NMSE between simulated and actual vegetationpatterns at different scales

3.2 植物群落演替的時間變化模擬結(jié)果

沉陷裂縫區(qū)的植物群落在自然恢復初期以草本植被為主，在自然恢復約3 a后開始出現(xiàn)灌木植被。圖12反映了沉陷裂縫區(qū)內(nèi)不同類型植被覆蓋度在自然恢復10 a間的動態(tài)變化，圖13展示了植被自然恢復模擬過程中年度生長期最后一天的植被覆蓋情況。在自然恢復初期，草本植物群落的發(fā)展速度較快，灌木植被由于對土壤養(yǎng)分和種子數(shù)量的要求較高，在自然恢復的前3 a間無灌木植被出現(xiàn)。但隨著草本植被凋落死亡并轉(zhuǎn)變?yōu)橥寥鲤B(yǎng)分，越來越多的區(qū)域內(nèi)滿足灌木的生長條件，因此在自然恢復5 a后灌木植被表現(xiàn)出明顯的增長趨勢。在本次研究中，群落演替模擬的停止標志為植被覆蓋度介于35%～55%，灌木覆蓋度介于7%～8%，50次模擬結(jié)果表明到達這一植被覆蓋條件評價平均需要10.46 a，這與以往研究中的植被自然恢復時間相似。在以往的自然恢復實例中發(fā)現(xiàn)，植物群落的自然恢復基本遵循草本群落→草本-灌木群落→草本-灌木-喬木群落的演替順序，其中灌木植被一般在自然恢復3～5 a后出現(xiàn)，而喬木植被的恢復需要十年乃至幾十年。本次研究的50次模擬結(jié)果中，僅2次模擬結(jié)果在自然恢復10 a后開始出現(xiàn)喬木植被，這與實際的喬木植被自然恢復規(guī)律基本一致。

圖12 模擬自然恢復過程中植被覆蓋度的時序變化Fig.12 Dynamic of vegetation coverage during succession

圖13 植被自然恢復過程中每一年度生長期末的植被覆蓋情況Fig.13 Vegetation coverage at the last day of growth peroid in each year duing spontaneous succession

3.3 植物群落演替的空間分布模擬結(jié)果

沉陷裂縫區(qū)內(nèi)低矮的溝壑區(qū)域由于水分與養(yǎng)分的匯集，成為了灌木植被的聚集區(qū)域。通過對比植被自然恢復模擬結(jié)果和沉陷裂縫區(qū)的地形條件可以發(fā)現(xiàn)，灌木植被主要集中分布在高程較低、坡度較小且坡向為北向的區(qū)域。由于本次研究中設(shè)置的初始土壤養(yǎng)分條件僅滿足草本植被的生長，因此在植被自然恢復初期不存在灌木植被。但隨著土壤獲取草本植被凋落物而積累養(yǎng)分的過程中，地表徑流通過裹挾區(qū)域內(nèi)的土壤養(yǎng)分，在重力作用下不斷向低矮區(qū)域匯集使其成為了養(yǎng)分和水分的富集區(qū)，因此高程較低的區(qū)域往往率先滿足了灌木植被的生長條件。而坡度較小的區(qū)域則由于土壤入滲速率、徑流速率較小，相較于坡度較大的土壤也夠儲存更多的水分和養(yǎng)分。坡向則主要影響著土壤的蒸發(fā)速率，這使得面向北方的陰坡更適宜植被的生長。因此不難發(fā)現(xiàn)，干旱、半干旱礦區(qū)沉陷裂縫的出現(xiàn)，在一定程度上反而促進了灌木植被的生長。沉陷裂縫一方面導致了部分區(qū)域成為了徑流和養(yǎng)分匯集的富集區(qū)，同時也增加了周邊區(qū)域的坡度，提高了徑流速率而降低了土壤入滲速率，加劇了干旱半干旱礦區(qū)的水土流失現(xiàn)象，逐漸形成了植物群落集中于低矮溝壑區(qū)域的空間分布格局。

4 結(jié) 論

(1)基于開發(fā)的群落演替模擬系統(tǒng)模擬了大同市云岡礦區(qū)一處沉陷裂縫區(qū)的自然恢復后形成的植被格局，模擬結(jié)果與實際植被格局間的誤差在0.04～0.59，其中灌木植被格局誤差在0.06～0.59，草本植被格局誤差在0.05～0.14，整體植被格局誤差在0.04～0.12，模擬結(jié)果與實際的植被空間分布規(guī)律一致。

(2)以裸地為初始條件的50次模擬結(jié)果表明，研究區(qū)通過自然恢復使草本覆蓋度達到35%～55%且灌木覆蓋度達到7%～8%水平時，平均需要10.46 a，模擬結(jié)果與實際的植被動態(tài)變化規(guī)律一致。

(3)模擬結(jié)果中草本植被在自然恢復前期快速發(fā)展，灌木植被在恢復3 a后開始出現(xiàn)，且在5 a后增長加速，但喬木植被在50次模擬中僅有2次在自然恢復10 a后出現(xiàn)，模擬結(jié)果與實際的植物群落演替規(guī)律一致。