基于改進的Lotka-Volterra種間競爭模型預測退化高寒草地人工恢復演替結果

吳曉慧,單熙凱,董世魁,*,高曉霞,許馭丹,張 曦,武勝男,胡 櫻,溫 璐,劉世梁,董全民,王文穎

1 北京師范大學環(huán)境學院, 北京 100875 2 北京師范大學天文系, 北京 100875 3 內蒙古大學生態(tài)與環(huán)境學院, 呼和浩特 010021 4 青海師范大學生命科學學院, 西寧 810008 5 青海大學畜牧獸醫(yī)科學院, 西寧 810016

青藏高原的草地面積約占全國草地總面積的1/3[1],提供了牧業(yè)生產基地和生態(tài)屏障等生態(tài)服務功能。但是,由于受到氣候條件的變化和人類活動的影響,該區(qū)50%—70%的天然草地已經退化[2]。為了有效遏制天然草地退化、提高草地生產能力,行之有效的辦法之一,就是在失去牧用價值,難以自我恢復的極度退化草地上建植和培育高產、優(yōu)質、穩(wěn)定的人工草地[2],降低天然草地放牧壓力,并通過人工草地植被逐漸恢復退化的天然草地。自20世紀80年代以來大面積的“黑土型”次生裸地的出現[3],人為干預下的高寒草甸恢復演替過程成為生態(tài)學家研究的熱點問題,目前的研究主要集中在人工恢復與其他恢復措施間的差異性[4]、播種牧草的優(yōu)化組合[5]及混播機理[6]、水肥調控[7]、群落生產力及穩(wěn)定性變化[8]等,但對人工草地植物群落種間競爭關系及其生態(tài)效應的研究較為少見。在人工草地的演替過程中,競爭是塑造植被群落結構和動態(tài)特征的動力[9]。尤其是在青藏高原的人工草地植物群落中,由于各組分生態(tài)位高度重疊,種間競爭尤為激烈,植物群落穩(wěn)定性難以維持,植物群落演替何時趨于穩(wěn)定難以預測,各組分之間的競爭成為控制群落演替的主導因子[2]。因此,對人工草地各植物組分的競爭效應進行研究,將有助于我們厘清植被恢復過程、預測植物群落演替的方向,可以為退化高寒草地的人工恢復管理提供科學依據。

目前,利用生態(tài)數學方法構建的諸多數學模型在草地植被研究中的應用越來越廣泛。其中,反應競爭-密度效應的數學表達模型以Logistic模型較為常見,但是該模型僅適用于實驗室控制條件的競爭試驗結果預測,對自然環(huán)境中的種間競爭關系預測較為困難。因此出現了許多改進模型,如Hallam模型[10]、Smith模型[11]等。從多物種競爭的角度出發(fā),Lotka-Volterra種間競爭模型(簡稱LVM)是Logistic模型的一種延伸,其應用較為廣泛。例如,王慧忠等[12]利用LVM求解牧草種群與放牧家畜種群的變化規(guī)律,林慧龍等[2]建立了適合高寒地區(qū)3組分禾草混播種間競爭機制的LVM模型。但是,由于青藏高原高寒地區(qū)植物群落演替的長期追蹤監(jiān)測較為困難,長時間序列的植物群落演替過程中的種間競爭關系鮮有報道。對于植物群落中種間關系的預測,也需要對傳統(tǒng)的LVM模型進行修正,才能應用于自然條件下植物群落演替過程中的種間競爭關系預測。

本研究旨在通過野外調查和模型模擬,建立適合黑土灘型高寒人工草地植物競爭效應的LVM修正模型,從牧業(yè)生產和生態(tài)功能兩個角度出發(fā),將高寒草地植物群落的組分植物進行歸類,結合不同恢復年限人工草地的實測組分植物的數量特征,通過計算機模擬研究改進后Lotka-Volterra系統(tǒng)的動力學行為,刻畫不同恢復年限高寒草地植物群落各組分的時間演替特征,分析高寒草地不同組分植物的競爭結局,預測退化高寒草地人工恢復所需的時間限值,從而為退化高寒草地的人工恢復管理提供科學依據。

1 研究方法

1.1 研究區(qū)域概況

研究區(qū)域位于青海省果洛藏族自治州瑪沁縣大武鎮(zhèn),屬于典型的高原山地類型。地理位置為34°21′—34°29′N,100°12′—100°29′E,平均海拔3800 m。氣候類型為高原大陸性氣候,春季干旱多風,冬季寒冷漫長,年均溫約-3.9℃,年均降水量423—565 mm,降水多集中在6月—9月份,占年降水量60.0%。牧草生長季約156 d,無絕對霜期[13]。全年日照時間2313—2607 h,太陽輻射強,晝夜溫差大。

土壤類型為高山草甸土,85%的草地為高寒草甸,植被多樣性程度高,為發(fā)展畜牧業(yè)提供了重要資源。但受到人類活動諸如超載放牧、過度開墾等影響,加之氣候變化帶來的影響,大部分天然草地出現不同程度的退化,黑土灘問題突出[3]。黑土灘退化草地給區(qū)域生態(tài)、生產、社會經濟等帶來負面效應[14]。已有研究結果表明,黑土灘退化草地使用圍欄封育自然恢復方法難以達到恢復標準,且最終會形成越來越穩(wěn)定的毒雜草群落[15],而人工草地措施在一定程度上促進了黑土灘退化草地植物群落的恢復演替[16]。為了解決草地退化問題,該地區(qū)建植了大片人工草地,目前已有不同建植年限的人工草地。

1.2 數據來源

于2017年7—8月進行野外采樣工作。在選取的研究區(qū)域內,選擇13塊不同建植年限的研究樣地,均為冬春草場。這13塊樣地均為禾本科——垂穗披堿草(Elymusnutans)或草地早熟禾(Poapratensis)單播的人工草地,建植期分別為17a、15a、13a、12a、11a、10a、9a、8a、7a、5a、3a、2a、0a(春季新翻撒播地)。我們采用空間分布代替時間演替[17]的方法研究人工草地恢復演替過程中各組分的動態(tài)競爭。樣地信息見表1。

表1 樣地基本信息

在上述13個樣地中,分別隨機選取5個1 m×1 m的樣方,統(tǒng)計各樣方中每一個物種的物種數量(多度)。

1.3 數據處理

本文選取物種的種群數量來表征物種競爭力。將樣地上的所有物種按不同的目的進行分組,基于統(tǒng)計的各樣方中每一個物種的種群數量,利用Python計算出每一組分的種群數量之和。采用3σ準則對每個樣地上的五個樣方的數據進行綜合,將數據取平均值作為該樣地上各組分的種群數量之和。

2 模型的建立

LVM是研究種群對某些公共資源的競爭動態(tài)的簡單模型,是經典的單種群增長的Logistic方程在種間競爭上的直接推廣,因此又被稱為種間競爭的Logistic方程[18]。已知LVM如下：

(1)

式中,N1、N2分別為兩個物種的種群數量；K1、K2分別為兩個物種的環(huán)境容納量；r1、r2分別為兩個物種的種群增長率；α為物種2對物種1的競爭系數,即每個N2個體所占用的空間相當于α個N1個體所占用空間；β為物種1對物種2的競爭系數,即每個N1個體所占用的空間相當于β個N2個體所占用空間。并且,N2種群中每個個體對N1種群的影響為α/K1,N1種群中每個個體對N2種群的影響為β/K2。

然而,LVM存在一定的局限性。模型中種群數量的改變與環(huán)境資源之間呈線性關系,不一定符合自然條件下的邏輯關系[19]。且實際情況中種群增長速度分布形式眾多,一般不對稱,與模型的預期存在較大的差別。只有當環(huán)境資源與種群所需要的資源供應數量相差很大時,才會出現方程中所描述的種群增長率隨種群密度的增加呈線性遞減的趨勢[19]。為了使模型對于種群競爭規(guī)律具有更大的包容性,借鑒張大勇[20]、李新運[21]等對Logistic方程提出的改進,我們對LVM進行非線性化,提出以下改進：

(2)

式中,i、j均為密度制約參數,取值范圍為0

2.1 模型解的存在唯一性

模型改進后,為求得微分方程組的數值解,首先需要證明該方程解的存在唯一性。為此,我們給出以下證明過程：

將方程(2)寫成矩陣形式即：

利用Lipschitz條件,得到以下證明：

因為

令

所以

因為

an+1-bn+1=(a-b)(an+an-1b+…+abn-1+bn)

令

g(a,b)=(an+an-1b+…+abn-1+bn)

所以

所以

因為

所以

同理可證

所以上式滿足Lipschitz條件,因此微分方程(2)存在唯一解,故改進后的模型在數學理論上可行。

2.2 模型的求解

在實際應用中,密度制約參數i、j的確定方法如下：

第二步：根據ODE算法求得微分方程(2)在不同i、j下的數值解；

第四步：循環(huán)Step 1—Step 3直至取遍所給定的所有i、j的組合；

第五步：畫出不同i和j的組合對應的相對誤差熱圖；

第六步：找出相對誤差最小的情況下所對應的i、j組合,并畫出此組合下的模型數值解N1(t),N2(t)趨勢圖。

3 模型的驗證

為了全面反映人工恢復草地植物的競爭格局,將人工草地植物物種組分分為以下幾組：(1)可食牧草與毒雜草；(2)栽培植物(垂穗披堿草、草地早熟禾)與原生植物；(3)先鋒植物與頂極植物[22]?；?3塊不同恢復年限人工草地的樣地實測數據,用Python進行樣條插值擴充數據量以保證擬合微分方程式系數時的穩(wěn)定性,并得到建植17年內的實際恢復演替圖。通過循環(huán)Step 1—Step 3得到確定的密度制約參數后,對以上改進的模型和LVM分別進行數值求解,得到改進的模型和LVM模擬的建植17年內的競爭演替圖。分別與實際恢復演替圖對比,得出改進模型對減小模型誤差的貢獻及其合理性。

需要注意的是,為了使相對誤差熱圖中誤差最小的部分(顏色最深處)得以突出體現,本文將該部分的誤差數值進行了適當調整。因此,相對誤差熱圖僅作對比參考以確定密度制約參數,具體的相對誤差數值在文中各個具體驗證部分展示。

3.1 可食牧草與毒雜草的競爭關系

對比圖1可食牧草與毒雜草模型的相對誤差熱圖中不同密度制約參數i和j下的誤差大小可以看出,在該分組中,i=0.5,j=1.5時,模型模擬值的平均相對誤差最小。利用Python求得計算機實驗模擬模型：

(3)

圖1 可食牧草與毒雜草模型相對誤差熱圖 Fig.1 Relative error thermography of edible forage and poisonous weed model

利用ODE算法求解微分方程(3),模擬建植17年內可食牧草與毒雜草的競爭情況并與實際恢復演替圖進行對比(圖2第一行)。同時對LVM進行求解,與實際恢復演替圖進行對比(圖2第二行),可以看出,改進后的模型模擬的結果比LVM模擬的結果更接近于實際恢復演替的情況。LVM模擬數值的平均相對誤差為30.98%,改進后模型模擬值的平均相對誤差為29.40%。其中,改進后模型對毒雜草的模擬值誤差僅為25.31%,優(yōu)于對可食牧草的模擬。模型改進后將平均相對誤差降低了5.10%。由此可見,對LVM進行非線性化的改進,可以實現可食牧草與毒雜草的競爭及演替結局的有效預測。

圖2 建植17年內可食牧草與毒雜草的模型預測與實際恢復演替對比圖Fig.2 A comparison between model predictions and actual restoration successes of edible forage and poisonous weed within 17 years上下兩行分別為改進后的模型和LVM的對比結果; 圖中帶有數據點的趨勢線為實際恢復演替圖,光滑的趨勢線為模型模擬的結果;X軸為恢復時間(年),Y軸為種群數量(株)

3.2 栽培植物與原生植物的競爭關系

對比圖3栽培植物與原生植物模型的相對誤差熱圖中不同密度制約參數i和j下的誤差大小,在該分組中,i=2.0,j=2.5時,模型模擬值相對誤差最小。利用Python求得計算機實驗模擬模型：

(4)

圖3 栽培植物與原生植物模型相對誤差熱圖 Fig.3 Relative error thermography of cultivated plants and native plants

同樣的,利用ODE算法求解微分方程(4),模擬建植17年內栽培植物與原生植物的競爭情況并與實際恢復演替圖進行對比(圖4第一行)。同時對LVM進行求解,與實際恢復演替圖進行對比(圖4第二行),結果表明,改進后的模型模擬的結果比LVM模擬的結果更接近于實際恢復演替,尤其是在第8年后的效果十分明顯。LVM模擬值的平均相對誤差為65.27%,改進后模型模擬值的平均相對誤差為55.24%。其中,改進后的模型對于原生植物的模擬誤差為30.33%,效果優(yōu)于栽培植物。模型改進后將相對誤差降低了15.37%。由此可見,對LVM進行非線性化的改進以預測栽培植物與原生植物未來的競爭結局并分析演替規(guī)律優(yōu)于LVM,對于較長時間尺度的預測是可行的,但仍需結合實際恢復演替圖綜合分析。

圖4 建植17年內栽培植物與原生植物的模型預測與實際恢復演替對比圖Fig.4 A comparison between model predictions and actual restoration successes of cultivated plants and native plants within 17 years上下兩行分別為改進后的模型和LVM的對比結果; 圖中帶有數據點的趨勢線為實際恢復演替圖,光滑的趨勢線為模型模擬的結果;X軸為恢復時間(年),Y軸為種群數量(株)

3.3 先鋒植物與頂極植物競爭關系

由圖5先鋒植物與頂極植物模型的相對誤差熱圖得出,在該分組中,密度制約參數i=0.5,j=1.0時,模型模擬值相對誤差最小。同樣利用Python求得計算機實驗模擬模型：

(5)

圖5 先鋒植物與頂極植物模型相對誤差熱圖 Fig.5 Relative error thermography of pioneer plants and climax plants

求解微分方程(5),模擬建植17年內先鋒植物與頂極植物的競爭情況并與實際恢復演替圖進行對比(圖6第一行)。同時對LVM進行求解,與實際恢復演替圖進行對比(圖6第二行)。結果表明,改進后的模型模擬的結果比LVM模擬的結果更接近于實際恢復演替。LVM模擬值的平均相對誤差為60.98%,改進后模型模擬值的平均相對誤差為52.88%。其中,對于頂極植物的模擬誤差為23.74%,效果優(yōu)于先鋒植物。模型改進后將相對誤差降低了13.28%。由此可見,對LVM進行非線性化的改進以預測先鋒植物與頂極植物未來的競爭結局并分析演替規(guī)律優(yōu)于LVM,對于較長時間尺度的趨勢預測是可行的。

圖6 建植17年內先鋒植物與頂極植物的模型預測與實際恢復演替對比圖Fig.6 A comparison between model predictions and actual restoration successes of pioneer plants and climax plants within 17 years上下兩行分別為改進后的模型和LVM的對比結果;圖中帶有數據點的趨勢線為實際恢復演替圖,光滑的趨勢線為模型模擬的結果;X軸為恢復時間(年),Y軸為種群數量(株)

4 競爭結果預測

4.1 可食牧草與毒雜草競爭結果

將模型模擬的預測年限延長至50年,得到自建植栽培植物后50年內可食牧草和毒雜草的競爭演替預測圖7。圖8為實際恢復演替圖。結合圖7、8可以看出,可食牧草在建植前7年內競爭力較強,呈直線上升,而7—10年急劇下降,在第10年達到最低,10年后呈現波動式上升；毒雜草在建植5年內種群數量不斷上升,5年后下降,后期呈小幅度波動式上升,總體競爭力弱于可食牧草。從畜牧生產的角度出發(fā),我們建議在建植7年后進行一定的人為干預，如施肥、滅鼠害以及適當去除雜草[23]等，以保證可食牧草持久的競爭力。從模型模擬預測圖7來看，建植第20年可食牧草與毒雜草的競爭開始逐漸趨于平衡，可恢復放牧生產，進行適度放牧。此外，通過Mathematica求解模型平衡點有：(0，0),(540.15，163.58),(0，816.33)。其中，滿足實際條件的穩(wěn)態(tài)平衡點為(540.15，163.58)。通過該平衡點對應預測演替數據得知在人工草地建植后第28年，可食牧草與毒雜草的競爭達到平衡。

圖7 可食牧草與毒雜草競爭演替預測圖 Fig.7 Competition succession prediction of edible forage and poisonous weed

圖8 可食牧草與毒雜草實際恢復演替圖Fig.8 Actual restoration succession of edible forage and poisonous weed

4.2 栽培植物與原生植物競爭結果

人工恢復草地中的栽培植物與原生植物的競爭結果預測如圖9所示。從實際恢復演替圖10可以看出，起初栽培植物的競爭力較弱，呈現一定的下降趨勢，但在第5年開始上升并與自然植物發(fā)生激烈的競爭，且二者之間的競爭比可食牧草與毒雜草之間的競爭更為激烈，尤其是在建植后5—8年之間栽培植物競爭力明顯強于原生植物。人工草地建植8年后，栽培植物的競爭力下降，原生植物的競爭力上升；人工草地建植10年后二者協同生長。人工草地建植20年后，二者之間的競爭逐漸趨于動態(tài)平衡。因此，人工草地建植最初2—4a，要進行施肥、除雜草等人工輔助措施。

模型預測平衡點有(0，0)、(0，±710.70)、(-262.60±128.02i，-390.39±247.50i)、(-211.65±224.46i，-371.22±316.86i)、(297.55±47.91i，-352.43±272.22i)、(±688.43，0)、(245.83，456.46)、(-189.48±147.01i，477.03±23.03i)。其中滿足實際條件的穩(wěn)態(tài)平衡點為(245.83，456.46)。通過該平衡點對應預測演替數據得知，在人工草地建植第41年時，栽培植物和原生植物的競爭可以達到平衡。

圖9 栽培植物與原生植物競爭演替預測圖 Fig.9 Competition succession prediction of cultivated plants and native plants

圖10 栽培植物與原生植物實際恢復演替圖 Fig.10 Actual restoration succession of cultivated plants and native plants

4.3 先鋒植物與頂極植物的競爭結果

先鋒植物與頂極植物的競爭結果預測如圖11所示。競爭預測圖11和實際恢復演替圖12可以看出,在種植人工草種后5年內,先鋒植物和頂極植物協同增長,但在第5年后,頂極植物呈現出強烈的競爭優(yōu)勢,先鋒植物競爭力直線下降甚至于第8年接近于0,先鋒植物幾乎被完全替代。因此,在高寒草甸種植栽培植物對先鋒植物的生長抑制起到了一定的作用。然而第8年后,先鋒植物呈現小幅度的波動式上升。相應的,頂極植物競爭力呈現下降趨勢,第10年后繼而呈現波動式上升。從模型預測演替圖來看,自第20年起,二者之間的競爭逐漸趨于動態(tài)平衡,這一點與前文提及的可食牧草與毒雜草、栽培植物與原生植物的演替趨勢基本保持一致。

圖11 先鋒植物與頂極植物競爭演替預測圖 Fig.11 Competition succession prediction of pioneer plants and climax plants

圖12 先鋒植物與頂極植物實際恢復演替圖 Fig.12 Actual restoration succession of pioneer plants and climax plants

模型預測平衡點有(0,0)、(597.95,59.44)、(0,2085.03)。其中,符合實際條件的穩(wěn)態(tài)平衡點為(597.95,59.44),通過該平衡點對應預測演替數據得知在建植45年時,先鋒植物與頂極植物的競爭達到平衡。在人工草地建植、輔助恢復的狀態(tài)下,頂級植物的競爭力明顯高于先鋒植物。但預測結果顯示,先鋒植物在演替過程中保持一定微弱的競爭力,并未完全被頂極植物替代。

5 結論與討論

退化草地人工措施恢復的實踐和研究起步相對較晚,人工措施下草地恢復過程中的牧草群落演替規(guī)律尚未被完全探明[4]?？紤]到人工草地的恢復過程較長,因此,構建模型模擬各組分之間的競爭效應成為預測人工恢復結局的有效方法。目前,非線性動力學方法已經成為分析生態(tài)系統(tǒng)演變的主要工具[24]?；谇叭藢ogistic模型的改進,我們對LVM模型進行了非線性化,并利用Lipschitz條件證明得到其數學數值解的存在唯一性。

本文利用改進的LVM模型對青藏高原退化高寒草地人工恢復過程中對植物競爭和演替的結果預測,可以得出以下結論：(1)隨著恢復年限的增加,栽培植物(垂穗披堿草和草地早熟禾)與原生植物在經歷一個激烈的競爭階段后逐漸趨于動態(tài)平衡后趨于穩(wěn)定,說明人工種植垂穗披堿草、草地早熟禾等本土植物,可以促進青藏高原地區(qū)“黑土灘”型退化高寒草甸的有效恢復。這與馮瑞章[25]、楊慧茹[26]等的研究結果相一致,人工草地的群落組成、多樣性、結構特征變化以及植物量動態(tài)指標變化基本上趨于穩(wěn)定,土壤理化性質和種子性狀的變化與群落演替相適應,是趨于穩(wěn)定的動態(tài)變化,總體上趨于良性的方向發(fā)展[26]；(2)在“黑土灘”型退化草地上, 通過多年生人工草地建設, 短期內可以恢復草地植被[27]。但是,王長庭等研究發(fā)現,雜類草(尤其是甘肅馬先蒿)在建植4年的人工草地群落中開始局部入侵,在5—6年時入侵速度、入侵面積達到高峰期[24]。周華坤等的研究也表明人工草地在建成后4年內植物群落由“生產穩(wěn)定性”急劇向“生態(tài)穩(wěn)定性”轉化,呈現出明顯的退化態(tài)勢[28]。本文研究則發(fā)現可食牧草、栽培植物和頂極植物均從第7年開始競爭力出現明顯下降的趨勢。并且,王啟蘭等表明人工草地建植后第7年,原群落優(yōu)勢種的豐富度明顯下降, 主要以雜類草為主, 群落處于退化階段[29]。因此,人工草地的恢復管理中,建議建植后第4年及7—10年間進行適度的人為干預,如施肥、滅鼠害以及適當去除有毒有害雜草等[30]；(3)高寒地區(qū)“黑土灘”型退化草地的恢復需要較長的時間,王啟蘭等人在研究中發(fā)現,至少在建植后第14年,群落呈相對穩(wěn)定的狀態(tài)[29]。本文各植物組分的競爭預測結果則表明,“黑土灘”型退化草地進行人工恢復至少需要20年以上,才會獲得較為穩(wěn)定的植物群落。同時,有大量研究表明,放牧干擾對混播人工草地的穩(wěn)定性[31]、土壤理化性質[32]、草地群落特征[33]和牧草營養(yǎng)[34]等方面具有顯著影響。因此,從畜牧生產和草地恢復的角度出發(fā),應在恢復20年時進行適度放牧。

通過對青藏高原退化高寒草地人工恢復過程中植物競爭及演替預測來看,改進后的LVM模型存在以下優(yōu)點：(1)將模型非線性化使之成為了具有自適應性的非線性制約種群競爭模型,對于種群的增長、競爭模式具有更強的包容性；(2)從誤差熱圖可以看出,改進后的競爭模型相對誤差降低了10%以上,更接近于實際恢復演替的情況；(3)對模型的動力學行為的研究有助于預測人工草地各組分間的競爭結果,從而判斷其恢復演替狀況。但是,改進的模型也具有一定的缺點：(1)與實際演替圖相比,模型對短期競爭效應的模擬能力較弱,因此在實際應用中,需將模型預測與實際監(jiān)測結果相結合綜合分析；(2)對于密度制約參數i和j的確定,由于受到計算機運行速度和算法優(yōu)化等因素的制約,對i和j的取值間隔不夠密集,應采用更多、更密集的i和j進行循環(huán)模擬以尋求更小的相對誤差。在未來的研究中需要進一步發(fā)展和完善LVM模型,使其更加準確、合理地預測高寒草地人工恢復演替的結果。