種內競爭時滯對植被周期振蕩模式的影響研究*

李 靜， 孫桂全， 靳 禎

（1. 山西財經大學 應用數學學院， 太原 030006；2. 中北大學 理學院，太原 030051；3. 山西大學 復雜系統研究所， 太原 030006；4. 疾病防控的數學技術與大數據分析山西重點實驗室， 太原 030006）

引 言

干旱地區約占地球陸地總面積的41%，與干旱地區人口占世界總人口比例(38%)相近.由于人為因素和氣候變化，干旱和半干旱地區土地退化的面積預計將在未來幾十年內急劇增加[1].植被作為這種脆弱生態系統的主體，其生長規律、空間分布模式以及生物量覆蓋比例，可作為生態系統惡化和向好的早期指標[2-4].植被被稱為生態系統工程師，不僅對物種和生境的豐富度有重要影響，還對生態環境的改善和維持具有重要作用[5].從2002 年初國家在西部地區開始啟動退耕還林工程，到習總書記關于生態文明建設方面提出“綠水青山就是金山銀山”的科學論斷，進一步表明了對植被系統的研究變得更加必要和迫切.

植被的生存離不開水資源，土壤水分過低會抑制植被正常的生長發育.同時，植被的生長在一定程度上也會促進土壤水分分布結構發生變化.在干旱、半干旱地區，土壤水分已成為影響植被系統恢復和重建的主要限制因素.因而，利用數學動力學模型來研究植被與土壤水的作用機制引起了生態學家和數學家的廣泛關注[6-12].1999 年，Klausmeier[6]首次利用反應擴散方程構建了一個包含植被和土壤水的動力學模型去解釋非洲尼亞美附近地區植被形成的大尺度空間條紋狀結構，發現這種空間分布模式主要是由土壤水的流動和植被對水資源的競爭導致的.Von Hardenberg 等[7]為了研究水資源有限地區的植物根系吸水而引起植被斑塊對水的競爭作用，建立了一個在植被演化方程中考慮Holling-Ⅱ型的植被-土壤水反應擴散對流模型，解釋了水資源有限地區觀察到的點狀、條狀和帶狀分布模式產生的機制.Sun 等[11]解釋了土壤水擴散強度的降低會導致植被斑圖結構發生相變且經歷間隙狀、條狀和點狀分布結構，并且植被的增長依賴于反饋強度的臨界值.

目前植被-土壤水模型的數學分析主要集中在空間分布模式形成理論的研究上，從而去解釋自然界中各種植被空間分布結構的形成機制.事實上，自然界中還可觀察到植被的生物量會隨時間做周期振蕩現象.例如，根據植被覆蓋度1 公里的月數據發現甘肅省2010 年1 月至2018 年12 月植被覆蓋度演變隨年度推移呈現周期變化趨勢(圖1)，其中每一年的12 個柱體代表從1 月到12 月的植被覆蓋度.然而，這種現象無法利用空間分布模式形成理論進行研究，需要尋求新的理論方法.因此，本文將從另一種理論角度出發去解釋植被系統中存在的這類周期振蕩現象.

圖1 從2010 年1 月到2018 年12 月甘肅省植被覆蓋度隨年變化的柱狀圖Fig. 1 Time series of vegetation coverages in Gansu Province from Jan., 2010 to Dec., 2018

本文考慮到植被的種內競爭存在于幼年期與成年期的植被之間，因此建立具有種內競爭時滯的植被-土壤水模型，從新的理論角度研究了植被隨時間演化呈現周期振蕩現象的內在演化機制.本文的內容結構安排如下：第1 節構建具有種內競爭時滯的植被-土壤水動力學模型；第2 節展示對應的常微分系統有植被滅絕平衡點和植被存在平衡點，并分析出系統在唯一的植被存在平衡點附近發生Hopf 分支的條件，給出相應的臨界分支參數的表達式；第3 節基于數值模擬展示出系統在唯一的植被存在平衡點處產生了一族空間齊次的周期解，并發現種內競爭時滯影響空間周期振蕩模式的發生，同時找到對臨界分支參數、周期和振幅有顯著影響的關鍵參數；第4 節對全文內容進行了總結和討論.

1 動力學模型

不同于以往的經典植被模型[6]，文獻[7]考慮到植被的飽和吸水效應提出了一個在植被方程中引入Holling-Ⅱ型作用項的植被-土壤水反應擴散對流模型：

其中，n(x,t)，w(x,t)分別為t時刻x位置的植被和土壤水的密度，參數γ表示植被對土壤水的吸收率， σ表示功能反應參數， μ表示植被的死亡速率，p表示降雨量，b表示蒸發率，d1，d2分別表示植被和土壤水的擴散率，?2=?2/?x2表示一維空間上的Laplace 算子， β表示滲透層中植被根系引起的土壤水的擴散率.項γw/(1+σw)描述了植被的飽和吸水效應，?(1?bn)w表示土壤水的蒸發.利用數值模擬方法研究發現，模型產生了廣泛分布在干旱和半干旱地區的點狀、條狀和帶狀空間分布模式.進一步預測了植被系統發生了從低降雨量的裸地態轉變為高降雨量的均勻植被態的過程，其間產生不同的穩定狀態共存現象.

接下來，我們將主要從動力學行為分析的角度對模型(2)進行研究.

2 動力學理論

2.1 平衡點及局部穩定性

不考慮時滯和空間時，模型(2)轉變為如下的常微分系統：

表1 系統(3)正平衡點的存在情況Table 1 Existence properties of positive equilibrium points for system (3)

(Ⅰ) 當?≤0時，由f(0)=α0>0和limw→+∞f(w)=?∞，易知三次函數f(w)和w軸之間僅有一個位于w右半軸的交點(見圖2(a))，則f(w)=0有唯一一個正實數根w?>0.因此，系統(3)有一個唯一的正平衡點E?=(n?,w?).

圖2 方程(4)存在唯一正實數根的示意圖：(a) 情形(Ⅰ)；(b) 情形(Ⅱ)(i)；(c) 情形(Ⅱ)(ii)(其中紅色實心點代表正實數根，圖(a)中亮藍色和藍色曲線分別表示?=0和?<0的情形)Fig. 2 Schematic diagrams of the existence of the unique positive real root for eq. (4): (a) case (Ⅰ); (b) case (Ⅱ)(i); (c) case (Ⅱ)(ii), where the red dot represents the positive real root, the light blue and the blue curves of fig. (a) represent the cases of ?=0 and ?<0, respectively

圖3 平衡點的存在情況和零斜率線圖： (a) 植被生物量n隨著參數 p的變化；(b) 在w-n 平面上零斜率線f1(n,w)=0(紫色/洋紅色)和f2(n,w)=0(橘色)形成的交點(黑色實心點)，其中 E0 和E u表示植被滅絕平衡點和植被生存平衡點，參數γ=1.6，σ=0.8 ，μ=0.2，b=0.2Fig. 3 The existence of equilibria and zero-slope diagrams: (a) vegetation biomass n vs. p; (b) the intersection (black solid points) formed by zero-slope line f1(n,w)=0 (purple/magenta) and f2(n,w)=0 (orange) on w-n plane, where E0 and Eu are the vegetation extinction equilibrium points and the vegetation existence equilibrium points, γ=1.6，σ=0.8 ，μ=0.2，b=0.2

通過計算，可推出系統(3)在植被滅絕平衡點E0=(0,p)處的Jacobi 矩陣為

對應的特征方程為

容易得到特征值為λ1=?1和λ2=γp/(σp+1)?μ.若條件

2.2 非空間系統(2)的動力學行為

通過以上的分析可得到如下結論.

定理1 如果條件(C2)成立，

圖4 臨界時滯 τ00和植被生物量隨參數p 的變化情況： (a) τ00隨參數 p的變化情況；(b) (p,n)平面上的單參數分支圖，參數σ=0.8，μ=0.2，b=0.2，γ=1.6 ，τ=3Fig. 4 The changes of the critical delay τ00 and the vegetation biomass with parameter p: (a) the variation of τ00 with p; (b) bifurcation diagram in (p,n) plane for σ=0.8 , μ=0.2, b=0.2 , γ=1.6 and τ=3

2.3 系統(2)的動力學行為

當系統(2)不考慮時滯時，即k≠0和τ=0，對應的特征方程為

可以計算出

根據以上的分析內容，可獲得下述的結果：

定理2 如果條件(C3)和(C4)同時成立，則

定理1 和2 的結果說明非空間和空間植被系統均可發生植被隨時間演化呈現周期振蕩模式的現象.

3 數 值 模 擬

這一節通過對系統(2)在一維空間上執行一系列的數值模擬去驗證第2 節的理論結果.此外，進一步細致分析系統參數對臨界時滯值、周期振蕩模式的周期和振幅的影響程度.

3.1 非空間系統(2)的動力學行為

本小節通過數值模擬驗證定理1 的理論結果，并研究降雨量p對臨界時滯值 τ00、周期震蕩模式的周期和振幅的影響.主要關注非空間系統(2)發生Hopf 分支的行為，因此假設d1=d2=β=0.選擇參數σ=0.8， μ=0.2，b=0.2，γ=1.6，p=1.0，計算出Eu=(0.557,0.762)， τ00=2.41. 當τ=1<τ00時，圖5(a)展示了非空間系統(2)的解最終趨于唯一的植被存在平衡點Eu，表明Eu是局部漸近穩定的；當τ=3>τ00時，圖5(b)展示了植被生物量和土壤水的密度隨時間演化呈周期振蕩現象，表明非空間系統(2)在Eu處產生了一族周期解，并通過相圖(c)呈現極限環(藍色線和紅色點分別代表軌線和Eu).

圖5 非空間系統(2)解的數值模擬： (a) τ=1；(b) τ=3；(c) 相圖Fig. 5 Numerical simulation of the solution of non-spatial system (2): (a) τ=1; (b) τ=3; (c) the phase diagram

圖6 展示了臨界時滯 τ00隨著參數γ，b和 σ的變化情況.從圖6(a)可觀察到 τ00隨著參數γ的增加而快速減少，表明了植被對土壤水的吸收速率強度變大有利于植被與土壤水以周期振蕩的形式共存.從圖6(b)中觀察到τ00隨著參數b的增加而緩慢減少，說明蒸發率的增加會促進這種共存作用形式，但這種促進作用相比于土壤水的吸收速率是緩慢的.圖6(c)展示了 τ00隨著參數 σ增加而變大，反映了功能反應強度的增加反而抑制這種周期振蕩行為.同時，3 個圖一致呈現出 τ00隨著降雨量的增加而減小，這一結果說明降雨量與植被土壤水的吸收速率和蒸發率的效應一樣.

圖6 對于p=1.0,1.3,1.6 ，τ 00隨著γ , b, σ 的變化情況，其他參數為σ=0.8 ，μ=0.2，b=0.2 ，γ=1.6Fig. 6 The variations of τ 00 with γ , b, σ for p=1.0,1.3,1.6, with other parameters of σ=0.8 , μ=0.2, b=0.2, γ=1.6

3.2 空間系統(2)的動力學行為

本小節將通過數值模擬驗證定理2 的理論結果，主要關注降雨量參數p對臨界時滯、齊次空間周期振蕩模式的周期和振幅的影響.選取參數σ=0.8， μ=0.2，b=0.2， γ=1.6，p=1.0，d1=0.02，d2=0.2， β=2.計算出常數穩態解Eu=(0.557,0.762)，=4.64.當τ=1<時，圖7(a)展示了空間系統(2)的解最終趨于唯一的穩態解Eu，表明Eu是局部漸近穩定的.當τ=4.8>τ02時，圖7(b)、(c)展示了空間系統(2)的解最終趨于空間齊次的周期解，表明穩態解Eu失去穩定性，即從Eu處分支出一族漸進穩定的空間齊次周期解.為了更清晰地展示出周期振蕩行為，從圖7(b)中截取時間1200～1400 部分作為圖7(c)，通過將圖7(c)投影到時間和空間坐標面形成圖7(d)，呈現一種時空周期行為.

圖7 n(x,t)的時空演化圖：(a) τ=1；(b) τ=4.8；(c) 截取圖(b)形成的部分圖；(d) 圖(c)的二維時空圖Fig. 7 The time series of n(x,t): (a) τ=1; (b) τ=4.8; (c) the partial graph cut out of fig. (b); (d) the 2D spatiotemporal graph of fig. (c)

圖8 對于不同的k和β，隨著參數p的變化Fig. 8 The variations of with p for different k and β values

圖9 對于不同的參數 p ，隨著d 2，β ，γ ，b ，σ 的變化情況Fig. 9 The variations of with d 2, β , γ , b, σ for different p values

3.1 和3.2 小節的數值模擬結果已表明在特定的參數取值下非空間系統和空間系統都能發生周期振蕩模式.為了研究參數取值變化能否保持這種模式，以及辨識出臨界時滯值、周期振蕩模式的振幅和周期對參數變化的敏感性程度，我們將對兩種系統作參數局部敏感性分析[16]，其中局部敏感性系數(Slocal)的計算公式為

這里，Vpar表示兩個系統中參數的值，?par是對參數值做兩個方向的微小擾動(+/?10%)，VNR和?NR分別對應于Vpar 和?par的數值模擬結果，如τ00/τ02、振幅和周期的模擬結果.為了評估參數的影響程度，作以下規定：敏感性系數大于2 時，此參數為高敏感參數；敏感性系數介于1 和2 之間為中度敏感參數，否則為低敏感參數.

在圖10 中，第一行和第二行分別描述了非空間系統和空間系統對應的τ00/τ02、振幅和周期對每個參數的敏感性程度，可以看出參數γ始終是高敏感參數而b是低敏感性參數，表明植被對土壤水的吸收率對植被系統的震蕩模式影響最顯著，但系統對參數b的 魯棒性最強.從圖10(a)、(c)、(d)和(f)中可觀察到 σ 和μ對 τc和周期來說是低敏感參數，但對于振幅來說是中度敏感性參數.此外，由圖10(d)～(f)可知， τ02受土壤水的擴散速率d2和植被根系引起的土壤水的擴散速率β的影響較大，而受植被的擴散速率d1的影響小，但這三個參數對振幅和周期沒有任何影響，表明空間擴散因素會影響周期振蕩模式的發生.此外，降雨量p的影響比較復雜，對 τ02、非空間周期、空間振幅是中度敏感性參數，對非空間振幅是高敏感參數，對 τ00、空間周期是低敏感參數.

圖10 周期振蕩模式的敏感性分析，參數σ=0.8 ，μ=0.2，b=0.2 ，γ=1.6 ，d1=0.02 ，d2=0.2，β=2，k=2Fig. 10 Sensitivity analysis of periodic oscillation patterns with parameters of σ=0.8 , μ=0.2, b=0.2 , γ=1.6, d1=0.02, d2=0.2, β=2, k=2

4 結論和討論

考慮到干旱、半干旱地區的幼年植被與成年植被之間存在競爭土壤水資源的現象，本文構建了一個具有種內競爭時滯的植被-土壤水動力學模型.利用動力系統理論分析出對應的常微分系統存在一個植被滅絕平衡點E0，并給出唯一植被存在平衡點出現的判定條件.隨后給出非空間系統在唯一的植被存在平衡點Eu處發生Hopf 分支的條件并計算出獨立于空間擴散的臨界時滯值(j=1,2,···).同時可計算出依賴于空間擴散的臨界時滯值，發現當時滯值大于依賴于空間擴散的臨界時滯時，原來考慮空間擴散的系統在Eu處可以產生空間齊次的Hopf 分支周期解.數值模擬進一步驗證了植被隨時間的推移呈現周期振蕩模式，為解釋干旱、半干旱地區植被生物量的周期振蕩演化提供了新的理論視角，也可以為干旱、半干旱地區植被系統的保護和治理提供理論依據.

本文研究發現，降雨量、植被土壤水的吸收速率強度和蒸發率的增強均有利于非空間植被周期振蕩模式的發生，但植被與土壤水的功能反應強度的加大反而會產生抑制效應.另外，對具有空間擴散的系統，波數和植被根系引起的土壤水的擴散反而會抑制這種行為，表明了空間因素的引入使得植被系統不易發生周期振蕩模式現象，同時空間擴散的引入會影響其他參數對這種模式的影響.相比于無空間因素，影響效果會變得復雜.通過參數敏感性分析發現非空間參數均對和產生影響，但敏感性程度不一致，同時降雨量和植被的增長率影響尤為顯著；系統的非空間參數均可對非空間和空間振蕩模式的周期產生相同的影響效果，特別是植被增長率 γ和降雨量參數p的影響顯著；空間因素對依賴于空間的臨界時滯參數 τ00有影響，敏感性程度為β>d2>d1，但對空間齊次振蕩模式的周期和振幅沒有影響.文中利用Hopf 分支理論研究干旱、半干旱地區植被生物量隨時間演化做周期振蕩模式的內在機制，但未結合實際的干旱、半干旱地區的植被生物量數據進行討論，在后續工作中我們將對具體地區的植被展開相關的周期振蕩模式研究[17-18]，進而因地制宜地給出影響植被生物量的關鍵因素.

參考文獻（ References ） ：

［1］E IGENTLER L, SHERRATT J A. Metastability as a coexistence mechanism in a model for dryland vegetation patterns[J].Bulletin of Mathematical Biology, 2019, 81（7）: 2290-2322.

［2］P RINGLE R M, DOAK D F, BRODY A K, et al. Spatial pattern enhances ecosystem functioning in an African savanna[J].PLoS Biology, 2010, 8（5）: e1000377.

［3］G OWDA K, CHEN Y, IAMS S, et al. Assessing the robustness of spatial pattern sequences in a dryland vegetation model[J].Proceedings of the Royal Society A, 2016, 472: 20150893.

［4］G ETZIN S, YIZHAQ H, BELL B, et al. Discovery of fairy circles in Australia supports self-organization theory[J].Proceedings of the National Academy of Sciences of the United States of America, 2016, 113（13）:3551-3556.

［5］G ILAD E, VON HARDENBERG J, PROVEBZALE A, et al. Ecosystem engineers: from pattern formation to habit creation[J].Physical Review Letters, 2004, 93（9）: 098105.

［6］K LAUSMEIER C A. Regular and irregular patterns in semiarid vegetation[J].Science, 1999, 284（5421）: 1826-1828.

［7］V ON HARDENBERG J, MERON E, SHACHAK M, et al. Diversity of vegetation patterns and desertification[J].Physical Review Letters, 2001, 87（19）: 198101.

［8］L IU Q X, JIN Z, LI B L. Numerical investigation of spatial pattern in a vegetation model with feedback function[J].Journal of Theoretical Biology, 2008, 254（2）: 350-360.

［9］B ONACHELA J A, PRINGLE R M, SHEFFER E, et al. Termite mounds can increase the robustness of dryland ecosystems to climatic change[J].Science, 2015, 347（6222）: 651-655.

［10］T ARNITA C E, BONACHEL J A, SHEFFER E, et al. A theoretical foundation for multi-scale regular vegetation patterns[J].Nature, 2017, 541（7637）: 398-401.

［11］S UN G Q, WANG C H, CHANG L L, et al. Effects of feedback regulation on vegetation patterns in semi-arid environments[J].Applied Mathematical Modelling, 2018, 61: 200-215.

［12］L I J, SUN G Q, JIN Z. Interaction of time delay and spatial diffusion induce the periodic oscillation of the vegetation system[J].Discrete and Continuous Dynamical Systems(Series B), 2022, 27（4）: 2147-2172.

［13］曹 建智, 譚軍, 王培光. 一類具有時滯的云杉蚜蟲種群模型的Hopf分岔分析[J]. 應用數學和力學, 2019, 40（3）: 332-342. （CAO Jianzhi, TAN Jun, WANG Peiguang. Hopf bifurcation analysis of a model for spruce budworm populations with delays[J].Applied Mathematics and Mechanics, 2019, 40（3）: 332-342.（in Chinese））

［14］谷 雨萌, 黃明迪. 一類時間周期的時滯競爭系統行波解的存在性[J]. 應用數學和力學, 2020, 41（6）: 658-668. （GU Yumeng, HUANG Mingdi. Existence of periodic traveling waves for time-periodic Lotka-Volterra competition systems with delay[J].Applied Mathematics and Mechanics, 2020, 41（6）: 658-668.（in Chinese））

［15］歐 陽頎. 非線性科學與斑圖動力學導論[M]. 北京: 北京大學出版社, 2010. （OUYANG Qi.Introduction to Nonlinear Science and Pattern Dynamics[M]. Beijing: Peking University Press, 2010. （in Chinese））

［16］I NGALLS B. Sensitivity analysis: from model parameters to system behaviour[J].Essays in Biochemistry, 2008,45: 177-194.

［17］K éFI S, RIETKERK M, ALADOS C L, et al. Spatial vegetation patterns and imminent desertification in Mediterranean arid ecosystems[J].Nature, 2007, 449（7159）: 213-217.

［18］B ASTIAANSEN R, JA?BI O, DEBLAUWE V, et al. Multistability of model and real dryland ecosystems through spatial self-organization[J].Proceedings of the National Academy of Sciences, 2018, 115（44）: 11256-11261.