李 靜， 孫桂全， 靳 禎

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

引 言

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

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

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

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

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

1 動力學模型

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

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

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

2 動力學理論

2.1 平衡點及局部穩(wěn)定性

不考慮時滯和空間時，模型(2)轉變?yōu)槿缦碌某Ｎ⒎窒到y(tǒng)：

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

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

圖2 方程(4)存在唯一正實數(shù)根的示意圖：(a) 情形(Ⅰ)；(b) 情形(Ⅱ)(i)；(c) 情形(Ⅱ)(ii)(其中紅色實心點代表正實數(shù)根，圖(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隨著參數(shù) p的變化；(b) 在w-n 平面上零斜率線f1(n,w)=0(紫色/洋紅色)和f2(n,w)=0(橘色)形成的交點(黑色實心點)，其中 E0 和E u表示植被滅絕平衡點和植被生存平衡點，參數(shù)γ=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

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

對應的特征方程為

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

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

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

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

圖4 臨界時滯 τ00和植被生物量隨參數(shù)p 的變化情況： (a) τ00隨參數(shù) p的變化情況；(b) (p,n)平面上的單參數(shù)分支圖，參數(shù)σ=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 系統(tǒng)(2)的動力學行為

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

可以計算出

根據(jù)以上的分析內(nèi)容，可獲得下述的結果：

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

定理1 和2 的結果說明非空間和空間植被系統(tǒng)均可發(fā)生植被隨時間演化呈現(xiàn)周期振蕩模式的現(xiàn)象.

3 數(shù) 值 模 擬

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

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

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

圖5 非空間系統(tǒng)(2)解的數(shù)值模擬： (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隨著參數(shù)γ，b和 σ的變化情況.從圖6(a)可觀察到 τ00隨著參數(shù)γ的增加而快速減少，表明了植被對土壤水的吸收速率強度變大有利于植被與土壤水以周期振蕩的形式共存.從圖6(b)中觀察到τ00隨著參數(shù)b的增加而緩慢減少，說明蒸發(fā)率的增加會促進這種共存作用形式，但這種促進作用相比于土壤水的吸收速率是緩慢的.圖6(c)展示了 τ00隨著參數(shù) σ增加而變大，反映了功能反應強度的增加反而抑制這種周期振蕩行為.同時，3 個圖一致呈現(xiàn)出 τ00隨著降雨量的增加而減小，這一結果說明降雨量與植被土壤水的吸收速率和蒸發(fā)率的效應一樣.

圖6 對于p=1.0,1.3,1.6 ，τ 00隨著γ , b, σ 的變化情況，其他參數(shù)為σ=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 空間系統(tǒng)(2)的動力學行為

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

圖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和β，隨著參數(shù)p的變化Fig. 8 The variations of with p for different k and β values

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

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

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

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

圖10 周期振蕩模式的敏感性分析，參數(shù)σ=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 結論和討論

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

本文研究發(fā)現(xiàn)，降雨量、植被土壤水的吸收速率強度和蒸發(fā)率的增強均有利于非空間植被周期振蕩模式的發(fā)生，但植被與土壤水的功能反應強度的加大反而會產(chǎn)生抑制效應.另外，對具有空間擴散的系統(tǒng)，波數(shù)和植被根系引起的土壤水的擴散反而會抑制這種行為，表明了空間因素的引入使得植被系統(tǒng)不易發(fā)生周期振蕩模式現(xiàn)象，同時空間擴散的引入會影響其他參數(shù)對這種模式的影響.相比于無空間因素，影響效果會變得復雜.通過參數(shù)敏感性分析發(fā)現(xiàn)非空間參數(shù)均對和產(chǎn)生影響，但敏感性程度不一致，同時降雨量和植被的增長率影響尤為顯著；系統(tǒng)的非空間參數(shù)均可對非空間和空間振蕩模式的周期產(chǎn)生相同的影響效果，特別是植被增長率 γ和降雨量參數(shù)p的影響顯著；空間因素對依賴于空間的臨界時滯參數(shù) τ00有影響，敏感性程度為β>d2>d1，但對空間齊次振蕩模式的周期和振幅沒有影響.文中利用Hopf 分支理論研究干旱、半干旱地區(qū)植被生物量隨時間演化做周期振蕩模式的內(nèi)在機制，但未結合實際的干旱、半干旱地區(qū)的植被生物量數(shù)據(jù)進行討論，在后續(xù)工作中我們將對具體地區(qū)的植被展開相關的周期振蕩模式研究[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]. 應用數(shù)學和力學, 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］谷 雨萌, 黃明迪. 一類時間周期的時滯競爭系統(tǒng)行波解的存在性[J]. 應用數(shù)學和力學, 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.