具交叉对流和局部放牧Klausmeier模型的图灵斑图

Turing Patterns of the Klausmeier Model with Cross⁃Advection and Local Grazing

摘要

为了研究交叉对流对斑图形成的影响， 将缓坡上的水流携带植被种子的交叉对流项引入到Klausmeier模型中， 构建了含交叉对流和局部放牧项的植被-水模型。首先， 得到了常微分系统的正平衡点的存在性和稳定性； 其次， 基于常微分系统稳定的正平衡点， 利用线性化方法和特征值理论得到了对应反应扩散系统正平衡点的不稳定性， 从而得到图灵不稳定性； 最后， 通过数值模拟验证理论结果， 分析交叉对流项和局部放牧项的存在对分岔和斑图的影响。结果显示， 交叉对流项和局部放牧项的存在改变了系统发生分岔的难度， 交叉对流项的存在加快了稳定带状斑图的形成速度， 同时减缓了植被带的移动速度。

Abstract

In order to investigate the effects of cross-advection on pattern formation， the cross-advection term of water flow carrying vegetation seeds on slopes was introduced into the Klausmeier model. A new vegetation-water model with local grazing term and cross-advection term was established. Firstly， the existence and stability of the positive equilibrium for the ordinary differential system were obtained. Secondly， based on the stable positive equilibrium obtained， the instability of the positive equilibrium of the reaction diffusion system was obtained by linearization method and eigenvalue theory， and the Turing instability was obtained. Finally， numerical simulation was presented using Matlab to support the theoretical result， and the influence of cross-advection term and local grazing term on bifurcation and pattern was investigated. The results show that the difficulty of bifurcation is changed by the cross-advection term and local grazing term， and the cross-advection term makes the speed of the formation of stable stripe patterns more quickly， but the migration speed of vegetation bands more slowly.

Graphical abstract

关键词

反应扩散方程 / 图灵斑图 / 图灵不稳定性 / 交叉对流 / 局部放牧

Key words

reaction-diffusion equations / Turing patterns / Turing instability / cross-advection / local grazing

0 引　言

半干旱地区由于水资源有限， 植被破坏， 生态环境比较恶劣， 所以半干旱地区生态环境的平衡状态很容易受到破坏， 植被退化造成土壤荒漠化^［1］。例如， 宁夏盐池县常年自由放牧的天然草场^［2］， 属于中温带半干旱地区， 大地貌为缓坡丘陵， 有大面积的风沙土、 黄绵土等非地带性土壤。干旱和半干旱地区在全球陆地生态系统中占有重要地位， 其气候特征和生态环境受到多种自然和人为因素的影响， 了解这些区域的气候变化特征及其影响因素对于制定有效的生态保护和恢复措施具有重要意义^［3］。Turing提出的图灵不稳定性是指均匀稳定的平衡点在引入扩散作用后， 反而变得不稳定， 进而促使系统自发形成空间非均匀的有序结构（如生物中的花纹、 斑图等）的现象， 图灵不稳定性的讨论必须建立在“无扩散时平衡点稳定”的前提下， Turing 提出的图灵斑图理论为理解干旱-半干旱地区植被空间格局的形成提供了重要的理论基础， 同时， 图灵斑图的研究和应用有助于预测干旱地区植被的演替趋势， 更详细的内容见文献［4］。在生态系统研究领域， 植被-水模型作为关键工具， 用于揭示植被与水分之间复杂的相互作用关系。植被斑图^［5］是植被在空间上的分布模式， 如均匀分布、 斑块状分布、 条带状分布等， 其形成与变化受到多种因素的综合影响， 其中水分的传输与分配在植被斑图的形成过程中起着核心作用。

1999年， Klausmeier^［6］提出了植被-水模型

\left\{\begin{array}{l}\frac{\partial W}{\partial T}=K_{\mathrm{1}}-K_{\mathrm{2}}W{B}^{\mathrm{2}}-K_{\mathrm{3}}W+\nu \frac{\partial W}{\partial X_{\mathrm{1}}},\\ \frac{\partial B}{\partial T}=D_B\mathrm{\Delta }B+K_{\mathrm{4}}{K}_{\mathrm{2}}W{B}^{\mathrm{2}}-M_{\mathrm{0}}B,\end{array}\right.

式中： W为地表水； B为植物量； K_{\mathrm{1}}，K_{\mathrm{2}}，K_{\mathrm{3}}，K_{\mathrm{4}}，M_{\mathrm{0}}都是正常数； K_{\mathrm{1}}为降水量； K_{\mathrm{2}}为植物对水的消耗率； K_{\mathrm{3}}为水的蒸发速率； K_{\mathrm{4}}为植物的增长率； M_{\mathrm{0}}为植物的死亡率； D_B\mathrm{\Delta }B表示植物的扩散， 其中， D_B为扩散系数， \mathrm{\Delta }=\frac{{\partial }^{\mathrm{2}}}{\partial X_{\mathrm{1}}^{\mathrm{2}}}+\frac{{\partial }^{\mathrm{2}}}{\partial X_{\mathrm{2}}^{\mathrm{2}}}； \nu \frac{\partial W}{\partial X_{\mathrm{1}}}表示水的径流， 其中， \nu代表水沿着下坡流的速度。

近年来， 国内外学者以经典的Klausmeier模型为基础构造出适合不同生态系统规律的植被斑图模型^［7‐16］， 例如： Hardenberg等^［10］考虑了水的扩散和植被对水流动的影响， 将Klausmeier模型中的地表径流量\nu \frac{\partial W}{\partial X_{\mathrm{1}}}改为\nu \frac{\partial (W-\alpha B)}{\partial X_{\mathrm{1}}}来描述水的流动； Sherratt^［11］将水的扩散项引入Klausmeier模型； Consolo等^［12］将交叉对流引入Klausmeier模型中； Consolo^［14］在广义双曲Klausmeier模型框架下研究了平坦干旱环境中Turing植被斑图的形成； Guo等^{［15‐16］}在具有植被饱和吸水效应的Klausmeier模型中加入了交叉扩散项等。但是， Klausmeier模型并没有考虑土壤中水和微地形的变化， 如缓坡对于植被生物量和地表径流的相互作用， 这里将植被生物量和地表径流的相互作用称为交叉对流。实际上， 在缓坡上不但植被对水的流动有影响而且植被种子也会随径流流动^［17］。一般情况下， 植物的损失分为两种， 一种是自然因素导致的植物死亡（如干旱、 气候变化、 衰老等）， 另一种是人为因素（如树木乱砍滥伐、 城市化、 过度放牧、 草原退化等）。张鹏^［18］研究了具有交叉对流的改进的Klausmeier模型

\left\{\begin{array}{l}\frac{\partial v}{\partial t}=w{v}^{\mathrm{2}}-mv-a_v\frac{\partial w}{\partial x}+{\nabla }^{\mathrm{2}}v,\\ \frac{\partial w}{\partial t}=a-w-w{v}^{\mathrm{2}}-a_w\frac{\partial w}{\partial x}+d_w{\nabla }^{\mathrm{2}}w\mathrm{。}\end{array}\right.

在该模型中， 作者考虑了坡面水流对植被种子传播的影响， 将交叉对流引入了Klausmeier模型， 分析了交叉对流对斑图形成的影响， 但是， 模型中只考虑了植物死亡的自然因素， 没有考虑人为因素。本文研究将考虑自然因素的植物死亡和人为因素的过度放牧导致的植物死亡。Siero等^［19-21］考虑了放牧导致的植被损失， 并构造出了具有非局部项和局部项的过度放牧导致的植被损失函数， 考虑到地域限制， 放牧区域具有局限性， 因此， 本文研究一类具有交叉对流和局部放牧的植被-水模型

\left\{\begin{array}{l}\frac{\partial w}{\partial t}=d_{\mathrm{1}}\mathrm{\Delta }w+R-w{v}^{\mathrm{2}}-w-a_w\frac{\partial w}{\partial x},\\ x\in \Omega ,t>\mathrm{0},\\ \frac{\partial v}{\partial t}=d_{\mathrm{2}}\mathrm{\Delta }v+w{v}^{\mathrm{2}}-m_{\mathrm{0}}v-\frac{m_{\mathrm{1}}{v}^{\mathrm{2}}}{\mathrm{1}+v^{\mathrm{2}}}-a_v\frac{\partial w}{\partial x},\\ x\in \Omega ,t>\mathrm{0},\\ \frac{\partial w}{\partial \eta }=\frac{\partial v}{\partial \eta }=\mathrm{0},x\in \partial \Omega ,t>\mathrm{0},\\ w(x,\mathrm{0})=w_{\mathrm{0}}(x)\ge \mathrm{0},b(x,\mathrm{0})=b_{\mathrm{0}}(x)\ge \mathrm{0},\\ x\in \overline{\Omega }，\end{array}\right.（1）

式中： \Omega =(\mathrm{0},l_{\mathrm{1}})\times (\mathrm{0},l_{\mathrm{2}})\subset {\mathfrak{R}}^{\mathrm{2}}为带有光滑边界\partial \mathrm{\Omega }的有界区域； \mathrm{\Delta }w,\mathrm{\Delta }v分别为水的扩散项和植被的扩散项； \mathrm{\Delta }=\frac{{\partial }^{\mathrm{2}}}{\partial x^{\mathrm{2}}}+\frac{{\partial }^{\mathrm{2}}}{\partial y^{\mathrm{2}}}； d_{\mathrm{1}},d_{\mathrm{2}}为水的扩散系数和植被的扩散系数； R,m_{\mathrm{0}}分别为降水量和植被的自然死亡率； \frac{m_{\mathrm{1}}{v}^{\mathrm{2}}}{\mathrm{1}+v^{\mathrm{2}}}为局部放牧导致的植被损失函数； m_{\mathrm{1}}为放牧导致的植被损失率； R，m_{\mathrm{0}}，m_{\mathrm{1}}都为正常数； a_w\frac{\partial w}{\partial x}和a_v\frac{\partial w}{\partial x}为交叉对流项， a_w为水沿着缓坡的水流速度， a_v为植被随水流方向的对流系数。

1 常微分模型的解的存在性和稳定性

模型（1）对应的常微分系统为

\left\{\begin{array}{l}\frac{\mathrm{d}w}{\mathrm{d}t}=R-w{v}^{\mathrm{2}}-w,\\ \frac{\mathrm{d}v}{\mathrm{d}t}=w{v}^{\mathrm{2}}-m_{\mathrm{0}}v-\frac{m_{\mathrm{1}}{v}^{\mathrm{2}}}{\mathrm{1}+v^{\mathrm{2}}}\mathrm{。}\end{array}\right.（2）

显然， 模型（2）有一个半平凡解(w_{\mathrm{0}},v_{\mathrm{0}})=(R,\mathrm{0})， 下面给出两个正平衡解的存在性。

定理 1　假设R>\mathrm{2}{m}_{\mathrm{0}}+m_{\mathrm{1}}， 则系统（2）有两个正常数平衡点

\begin{array}{c}{E}_{\mathrm{1}}=(w_{\mathrm{1}},v_{\mathrm{1}})=\\ \left(\frac{R}{\mathrm{2}A}\left(A+\sqrt[]{A^{\mathrm{2}}-\mathrm{4}{m}_{\mathrm{0}}^{\mathrm{2}}}\right),\frac{\mathrm{2}{m}_{\mathrm{0}}}{A+\sqrt[]{A^{\mathrm{2}}-\mathrm{4}{m}_{\mathrm{0}}^{\mathrm{2}}}}\right),\end{array}

\begin{array}{c}{E}_{\mathrm{2}}=(w_{\mathrm{2}},v_{\mathrm{2}})=\\ \left(\frac{R}{\mathrm{2}A}\left(A-\sqrt[]{A^{\mathrm{2}}-\mathrm{4}{m}_{\mathrm{0}}^{\mathrm{2}}}\right),\frac{\mathrm{2}{m}_{\mathrm{0}}}{A-\sqrt[]{A^{\mathrm{2}}-\mathrm{4}{m}_{\mathrm{0}}^{\mathrm{2}}}}\right),\end{array}

式中： A=R-m_{\mathrm{1}}。

证明　令模型（2）的右端都为0， 可得

\left\{\begin{array}{l}R-w{v}^{\mathrm{2}}-w=\mathrm{0},\\ wv-m_{\mathrm{0}}-\frac{m_{\mathrm{1}}v}{\mathrm{1}+v^{\mathrm{2}}}=\mathrm{0}\mathrm{。}\end{array}\right.（3）

由式（3）的第一个式子可得w=\frac{R}{\mathrm{1}+v^{\mathrm{2}}}， 将此式代入式（3）的第二个式子可得v_{\mathrm{1,2}}=\frac{\mathrm{2}{m}_{\mathrm{0}}}{A\pm \sqrt[]{A^{\mathrm{2}}-\mathrm{4}{m}_{\mathrm{0}}^{\mathrm{2}}}}， (A=R-m_{\mathrm{1}})， 将v_{\mathrm{1,2}}代入到w=\frac{R}{\mathrm{1}+v^{\mathrm{2}}}可得w_{\mathrm{1,2}}=\frac{R}{\mathrm{2}A}(A\pm \sqrt[]{A^{\mathrm{2}}-\mathrm{4}{m}_{\mathrm{0}}^{\mathrm{2}}})。由R>\mathrm{2}{m}_{\mathrm{0}}+m_{\mathrm{1}}可得R-m_{\mathrm{1}}=A>\mathrm{2}{m}_{\mathrm{0}}>\mathrm{0}， 即A^{\mathrm{2}}-\mathrm{4}{m}_{\mathrm{0}}^{\mathrm{2}}>\mathrm{0}， 所以， w_{\mathrm{1}}>\mathrm{0}， v_{\mathrm{1}}>\mathrm{0}； 显然， A-\sqrt[]{A^{\mathrm{2}}-\mathrm{4}{m}_{\mathrm{0}}^{\mathrm{2}}}>\mathrm{0}， 所以， w_{\mathrm{2}}>\mathrm{0}， v_{\mathrm{2}}>\mathrm{0}。综上， 当R>\mathrm{2}{m}_{\mathrm{0}}+m_{\mathrm{1}}时， E_{\mathrm{1}}=(w_{\mathrm{1}},v_{\mathrm{1}})和E_{\mathrm{2}}=(w_{\mathrm{2}},v_{\mathrm{2}})为正平衡点。证毕。

定理 2　假设系统（2）的参数R,m_{\mathrm{0}},m_{\mathrm{1}}都是正常数， 那么有以下结论：

1） 半平凡解(R,\mathrm{0})是稳定的；

2） 若R>\mathrm{2}{m}_{\mathrm{0}}+m_{\mathrm{1}}， 则E_{\mathrm{1}}不稳定；

3） 若R>\mathrm{2}{m}_{\mathrm{0}}+m_{\mathrm{1}}， 且满足条件时， 平衡点E_{\mathrm{2}}是稳定的。

证明　系统（2）对应的雅可比矩阵得

\left[\begin{array}{cc}-v^{\mathrm{2}}-\mathrm{1}& -\mathrm{2}vw\\ v^{\mathrm{2}}& \mathrm{2}wv-m_{\mathrm{0}}-\frac{\mathrm{2}{m}_{\mathrm{1}}v}{(\mathrm{1}+v^{\mathrm{2}}{)}^{\mathrm{2}}}\end{array}\right]\mathrm{。}（4）

将半平凡解E_{\mathrm{0}}=(R,\mathrm{0})代入到雅可比矩阵得

\mathit{J}{|}_{E_{\mathrm{0}}}=\left[\begin{array}{cc}-\mathrm{1}& \mathrm{0}\\ \mathrm{0}& -m_{\mathrm{0}}\end{array}\right],

则有特征方程

{\left|\lambda \mathit{E}-\mathit{J}\right|}_{E_{\mathrm{0}}}=\left|\begin{array}{cc}\lambda +\mathrm{1}& \mathrm{0}\\ \mathrm{0}& \lambda +m_{\mathrm{0}}\end{array}\right|=

{\lambda }^{\mathrm{2}}+(\mathrm{1}+m_{\mathrm{0}})\lambda +m_{\mathrm{0}}=\mathrm{0},

故。因此， 半平凡解E_{\mathrm{0}}是稳定的。

下面证明结论2）和结论3）。

当v\ne \mathrm{0}时， 根据式（3）可将式（4）简化为

\left[\begin{array}{cc}-v^{\mathrm{2}}-\mathrm{1}& -\mathrm{2}vw\\ v^{\mathrm{2}}& m_{\mathrm{0}}+\frac{\mathrm{2}{m}_{\mathrm{1}}{v}^{\mathrm{3}}}{(\mathrm{1}+v^{\mathrm{2}}{)}^{\mathrm{2}}}\end{array}\right]\mathrm{。}

对应的特征方程为

\left|\lambda \mathit{E}-\mathit{J}\right|=\left|\begin{array}{cc}\lambda +v^{\mathrm{2}}+\mathrm{1}& \mathrm{2}vw\\ -v^{\mathrm{2}}& \lambda -m_{\mathrm{0}}-\frac{\mathrm{2}{m}_{\mathrm{1}}{v}^{\mathrm{3}}}{(\mathrm{1}+v^{\mathrm{2}}{)}^{\mathrm{2}}}\end{array}\right|=

\begin{array}{c}{\lambda }^{\mathrm{2}}+(v^{\mathrm{2}}+\mathrm{1}+m_{\mathrm{0}}-\frac{\mathrm{2}{m}_{\mathrm{1}}{v}^{\mathrm{3}}}{(\mathrm{1}+v^{\mathrm{2}}{)}^{\mathrm{2}}})\lambda +\\ (v^{\mathrm{2}}+\mathrm{1})(-m_{\mathrm{0}}-\frac{\mathrm{2}{m}_{\mathrm{1}}{v}^{\mathrm{3}}}{(\mathrm{1}+v^{\mathrm{2}}{)}^{\mathrm{2}}})+\mathrm{2}{v}^{\mathrm{3}}w=\\ {\lambda }^{\mathrm{2}}+a_{\mathrm{1}}\lambda +a_{\mathrm{2}}=\mathrm{0}，\end{array}

其中，

a_{\mathrm{1}}=v^{\mathrm{2}}+\mathrm{1}-m_{\mathrm{0}}-\frac{\mathrm{2}{m}_{\mathrm{1}}{v}^{\mathrm{3}}}{(\mathrm{1}+v^{\mathrm{2}}{)}^{\mathrm{2}}}，

a_{\mathrm{2}}=m_{\mathrm{0}}(v^{\mathrm{2}}-\mathrm{1}),

所以， \mathrm{d}\mathrm{e}\mathrm{t}\mathit{J}=a_{\mathrm{2}}，\mathrm{t}\mathrm{r}\mathit{J}=-a_{\mathrm{1}}。

下面判断E_{\mathrm{1}}的稳定性。

因为A=R-m_{\mathrm{1}}>\mathrm{2}{m}_{\mathrm{0}}， 所以可得A+\sqrt[]{A^{\mathrm{2}}-\mathrm{4}{m_{\mathrm{0}}}^{\mathrm{2}}}>\mathrm{2}{m}_{\mathrm{0}}， 即\frac{\mathrm{2}{m}_{\mathrm{0}}}{A+\sqrt[]{A^{\mathrm{2}}-\mathrm{4}{m_{\mathrm{0}}}^{\mathrm{2}}}}， 因此， ， 则\mathrm{d}\mathrm{e}\mathrm{t}\mathit{J}{|}_{E_{\mathrm{1}}}=m_{\mathrm{0}}(v_{\mathrm{1}}^{\mathrm{2}}， 所以E_{\mathrm{1}}不稳定。

接下来判断E_{\mathrm{2}}的稳定性。

由于， 可得v_{\mathrm{2}}=\frac{\mathrm{2}{m}_{\mathrm{0}}}{A-\sqrt[]{A^{\mathrm{2}}-\mathrm{4}{m}_{\mathrm{0}}^{\mathrm{2}}}}=\frac{A+\sqrt[]{A^{\mathrm{2}}-\mathrm{4}{m}_{\mathrm{0}}^{\mathrm{2}}}}{\mathrm{2}{m}_{\mathrm{0}}}>\mathrm{1}， 所以， \mathrm{d}\mathrm{e}\mathrm{t}\mathit{J}{|}_{E_{\mathrm{2}}}=m_{\mathrm{0}}(v_{\mathrm{2}}^{\mathrm{2}}-\mathrm{1})>\mathrm{0}。

令\mathrm{\Delta }=A^{\mathrm{2}}-\mathrm{4}{m}_{\mathrm{0}}^{\mathrm{2}}， \mathrm{t}\mathrm{r}\mathit{J}{|}_{E_{\mathrm{2}}}=-a_{\mathrm{1}}{|}_{E_{\mathrm{2}}}， 因为A=R-m_{\mathrm{1}}>\mathrm{2}{m}_{\mathrm{0}}且， 则

\begin{array}{c}{a}_{\mathrm{1}}{|}_{E_{\mathrm{2}}}=\frac{A(A+\sqrt[]{\mathrm{\Delta }})-\mathrm{2}{m}_{\mathrm{0}}^{\mathrm{3}}}{\mathrm{2}{m}_{\mathrm{0}}^{\mathrm{2}}}-\frac{\frac{\mathrm{2}{m}_{\mathrm{1}}(A+\sqrt[]{\mathrm{\Delta }})}{\mathrm{8}{m}_{\mathrm{0}}^{\mathrm{3}}}}{\frac{A^{\mathrm{2}}{(A+\sqrt[]{\mathrm{\Delta }})}^{\mathrm{2}}}{\mathrm{4}{m}_{\mathrm{0}}^{\mathrm{4}}}}=\\ \frac{A(A+\sqrt[]{\mathrm{\Delta }})-\mathrm{2}{m}_{\mathrm{0}}^{\mathrm{3}}}{\mathrm{2}{m}_{\mathrm{0}}^{\mathrm{2}}}-\frac{m_{\mathrm{1}}{m}_{\mathrm{0}}(A+\sqrt[]{\mathrm{\Delta }})}{A^{\mathrm{2}}}>\end{array}\frac{\mathrm{2}{m}_{\mathrm{0}}(\mathrm{2}{m}_{\mathrm{0}}+\sqrt[]{\mathrm{\Delta }})-\mathrm{2}{m}_{\mathrm{0}}^{\mathrm{3}}}{\mathrm{2}{m}_{\mathrm{0}}^{\mathrm{2}}}-\frac{m_{\mathrm{1}}{m}_{\mathrm{0}}(\mathrm{2}{m}_{\mathrm{0}}+\sqrt[]{\mathrm{\Delta }})}{\mathrm{4}{m}_{\mathrm{0}}^{\mathrm{2}}}=\frac{m_{\mathrm{0}}(\mathrm{8}-\mathrm{4}{m}_{\mathrm{0}}-\mathrm{2}{m}_{\mathrm{1}})+(\mathrm{4}-m_{\mathrm{1}})\sqrt[]{\mathrm{\Delta }}}{\mathrm{4}{m}_{\mathrm{0}}}>\mathrm{0},

所以， ， 故E_{\mathrm{2}}是稳定的。证毕。

2 图灵分岔分析

下面给出系统（1）在平衡点E_{\mathrm{2}}处发生图灵不稳定性的充分条件。

定理 3　假设R>\mathrm{2}{m}_{\mathrm{0}}+m_{\mathrm{1}}， 且

M+{\left(\frac{P+\sqrt[]{P^{\mathrm{2}}+Q^{\mathrm{2}}}}{\mathrm{2}}\right)}^{\frac{\mathrm{1}}{\mathrm{2}}}>\mathrm{0},

则平衡点E_{\mathrm{2}}关于ODE系统（2）是稳定的， 但关于PDE系统（1）是不稳定的， 系统（1）的平衡点E_{\mathrm{2}}是图灵不稳定的， 即系统（1）发生图灵分岔。其中，

\begin{array}{c}M=g_v+f_w-(d_{\mathrm{1}}+d_{\mathrm{2}})(l_{\mathrm{1}}^{\mathrm{2}}+l_{\mathrm{2}}^{\mathrm{2}})，\\ P=(g_v-f_w+(d_{\mathrm{1}}-d_{\mathrm{2}})(l_{\mathrm{1}}^{\mathrm{2}}+l_{\mathrm{2}}^{\mathrm{2}}{))}^{\mathrm{2}}-a_w^{\mathrm{2}}{l}_{\mathrm{1}}^{\mathrm{2}}+\\ \hspace{1em}\hspace{1em}\mathrm{4}{f}_v{g}_w，\\ Q=\mathrm{2}{a}_w{l}_{\mathrm{1}}(g_v-f_w+(d_{\mathrm{1}}-d_{\mathrm{2}})(l_{\mathrm{1}}^{\mathrm{2}}+l_{\mathrm{2}}^{\mathrm{2}}))-\\ \hspace{1em}\hspace{1em}\mathrm{4}{f}_v{a}_v{l}_{\mathrm{1}}，\end{array}\begin{array}{c}{f}_w=-v_{\mathrm{2}}^{\mathrm{2}}-\mathrm{1}，f_v=-\mathrm{2}{w}_{\mathrm{2}}{v}_{\mathrm{2}}，\\ g_w=v_{\mathrm{2}}^{\mathrm{2}}\hspace{1em}，g_v=m_{\mathrm{0}}+\frac{\mathrm{2}{m}_{\mathrm{1}}{v}^{\mathrm{3}}}{(\mathrm{1}+v^{\mathrm{2}}{)}^{\mathrm{2}}}\mathrm{。}\end{array}

证明　首先对(w_{\mathrm{2}}(x,y,t),v_{\mathrm{2}}(x,y,t))进行空间非均匀的扰动：

w=w_{\mathrm{2}}+\theta (x,y,t)，v=v_{\mathrm{2}}+\epsilon (x,y,t)，（5）

式中： \theta (x,y,t)和\epsilon (x,y,t)表示t时刻空间斑块(x,y)上的扰动变量。

将式（5）代入到系统（1）中并在E_{\mathrm{2}}处进行线性化， 可得

\left\{\begin{array}{l}\frac{\partial \theta }{\partial t}=d_{\mathrm{1}}\mathrm{\Delta }\theta +f_w\theta +f_v\epsilon -a_w\frac{\partial \theta }{\partial x},\\ x\in \Omega ,t>\mathrm{0},\\ \frac{\partial \epsilon }{\partial t}=d_{\mathrm{2}}\mathrm{\Delta }\epsilon +g_w\theta +g_v\epsilon -a_v\frac{\partial \theta }{\partial x},\\ x\in \Omega ,t>\mathrm{0},\\ \frac{\partial \theta }{\partial \eta }=\frac{\partial \epsilon }{\partial \eta }=\mathrm{0},x\in \partial \Omega ,t>\mathrm{0},\\ \theta (x,\mathrm{0})={\theta }_{\mathrm{0}}(x)\ge \mathrm{0},\epsilon (x,\mathrm{0})={\epsilon }_{\mathrm{0}}(x)\ge \mathrm{0},x\in \overline{\Omega }\mathrm{。}\end{array}\right.（6）

令\theta =\tilde{\theta }(t)e^{i(l_{\mathrm{1}}x+l_{\mathrm{2}}y)}，\epsilon =\tilde{\epsilon }(t)e^{i(l_{\mathrm{1}}x+l_{\mathrm{2}}y)}， 其中， (l_{\mathrm{1}},l_{\mathrm{2}})为波向量， (x,y)为方向向量。代入式（6）可得

\left\{\begin{array}{l}\frac{\partial \tilde{\theta }}{\partial t}=-d_{\mathrm{1}}(l_{\mathrm{1}}^{\mathrm{2}}+l_{\mathrm{2}}^{\mathrm{2}})\tilde{\theta }+f_w\tilde{\theta }+f_v\tilde{\epsilon }-a_w{l}_{\mathrm{1}}i\tilde{\theta },\\ x\in \Omega ,t>\mathrm{0},\\ \frac{\partial \tilde{\epsilon }}{\partial t}=-d_{\mathrm{2}}(l_{\mathrm{1}}^{\mathrm{2}}+l_{\mathrm{2}}^{\mathrm{2}})\tilde{\epsilon }+g_w\tilde{\theta }+g_v\tilde{\epsilon }-a_v{l}_{\mathrm{1}}i\tilde{\theta },\\ x\in \Omega ,t>\mathrm{0},\\ \frac{\partial \tilde{\theta }}{\partial \eta }=\frac{\partial \tilde{\epsilon }}{\partial \eta }=\mathrm{0},x\in \partial \Omega ,t>\mathrm{0},\\ \tilde{\theta }(x,\mathrm{0})={\tilde{\theta }}_{\mathrm{0}}(x)\ge \mathrm{0},\tilde{\epsilon }(x,\mathrm{0})={\tilde{\epsilon }}_{\mathrm{0}}(x)\ge \mathrm{0},\\ x\in \overline{\Omega }\mathrm{。}\end{array}\right.（7）

系统（7）对应的雅可比矩阵为

\left[\begin{array}{cc}{f}_w-a_w{l}_{\mathrm{1}}i-d_{\mathrm{1}}(l_{\mathrm{1}}^{\mathrm{2}}+l_{\mathrm{2}}^{\mathrm{2}})& f_v\\ g_w-a_v{l}_{\mathrm{1}}i& g_v-d_{\mathrm{2}}(l_{\mathrm{1}}^{\mathrm{2}}+l_{\mathrm{2}}^{\mathrm{2}})\end{array}\right]\mathrm{。}

对应的特征值为

{\lambda }_{\mathrm{1,2}}=\frac{\mathrm{1}}{\mathrm{2}}(g_v+f_w-a_w{l}_{\mathrm{1}}i-(d_{\mathrm{1}}+d_{\mathrm{2}})(l_{\mathrm{1}}^{\mathrm{2}}+l_{\mathrm{2}}^{\mathrm{2}}))\pm

\begin{array}{c}\frac{\mathrm{1}}{\mathrm{2}}[(g_v-f_w+(d_{\mathrm{1}}-d_{\mathrm{2}})(l_{\mathrm{1}}^{\mathrm{2}}+l_{\mathrm{2}}^{\mathrm{2}})+a_w{l}_{\mathrm{1}}{i)}^{\mathrm{2}}-\\ \mathrm{4}{f}_v(a_v{l}_{\mathrm{1}}i-g_w{)]}^{\mathrm{1}/\mathrm{2}}\mathrm{。}\end{array}

分离实部和虚部可得

\mathrm{R}\mathrm{e}({\lambda }_{\mathrm{1,2}})=\frac{\mathrm{1}}{\mathrm{2}}\left(M\pm {\left(\frac{P+\sqrt[]{P^{\mathrm{2}}+Q^{\mathrm{2}}}}{\mathrm{2}}\right)}^{\frac{\mathrm{1}}{\mathrm{2}}}\right)，

\mathrm{I}\mathrm{m}({\lambda }_{\mathrm{1,2}})=\frac{\mathrm{1}}{\mathrm{2}}\left(-a_w{l}_{\mathrm{1}}\pm {\left(\frac{-P+\sqrt[]{P^{\mathrm{2}}+Q^{\mathrm{2}}}}{\mathrm{2}}\right)}^{\frac{\mathrm{1}}{\mathrm{2}}}\right)，

M=g_v+f_w-(d_{\mathrm{1}}+d_{\mathrm{2}})(l_{\mathrm{1}}^{\mathrm{2}}+l_{\mathrm{2}}^{\mathrm{2}})，

P=(g_v-f_w+(d_{\mathrm{1}}-d_{\mathrm{2}})(l_{\mathrm{1}}^{\mathrm{2}}+l_{\mathrm{2}}^{\mathrm{2}}{))}^{\mathrm{2}}-a_w^{\mathrm{2}}{l}_{\mathrm{1}}^{\mathrm{2}}＋\mathrm{4}{f}_v{g}_w，

Q=\mathrm{2}{a}_w{l}_{\mathrm{1}}(g_v-f_w+(d_{\mathrm{1}}-d_{\mathrm{2}})(l_{\mathrm{1}}^{\mathrm{2}}+l_{\mathrm{2}}^{\mathrm{2}}))－\mathrm{4}{f}_v{a}_v{l}_{\mathrm{1}}\mathrm{。}

特征值实部的最大值为

\mathrm{m}\mathrm{a}\mathrm{x}\mathrm{R}\mathrm{e}(\lambda )=\frac{\mathrm{1}}{\mathrm{2}}\left(M+{\left(\frac{P+\sqrt[]{P^{\mathrm{2}}+Q^{\mathrm{2}}}}{\mathrm{2}}\right)}^{\frac{\mathrm{1}}{\mathrm{2}}}\right)\mathrm{。}

由定理条件可得： \mathrm{m}\mathrm{a}\mathrm{x}\mathrm{R}\mathrm{e}(\lambda )>\mathrm{0}， 对于PDE系统（1）是不稳定的， 对于ODE系统（2）是稳定的， 因此系统（1）的正平衡点E_{\mathrm{2}}是图灵不稳定的， 即系统（1）发生图灵分岔。

3 数值模拟

为了探究交叉对流项和局部放牧项对干旱半干旱放牧地区植被斑图分岔点及斑图形成的影响， 设定固定参数值， 即R=\mathrm{3}， m_{\mathrm{0}}=\mathrm{0.4}， m_{\mathrm{1}}=\mathrm{1.5}， 由于只考虑交叉对流和局部放牧项的影响， 因此将水和植被的扩散系数d_{\mathrm{1}}和d_{\mathrm{2}}都固定为1， 扰动模数k=l_{\mathrm{1}}^{\mathrm{2}}+l_{\mathrm{2}}^{\mathrm{2}}， 通过调整相关变量开展分析。

3.1 交叉对流和局部放牧对分岔点的作用

在假定x和y方向扰动相等， 即l_{\mathrm{1}}=l_{\mathrm{2}}且有局部放牧项(m_{\mathrm{1}}=\mathrm{1.5})的情况下， 对交叉对流项a_v\frac{\partial w}{\partial x}中a_v=\mathrm{0}和a_v\ne \mathrm{0}两种情况进行图灵分岔分析。由图 1（a） 可知， 当不存在交叉对流项(a_v=\mathrm{0})且a_w=\mathrm{3}时， 特征值实部的最大值\mathrm{m}\mathrm{a}\mathrm{x}(\mathrm{R}\mathrm{e}(\lambda ))始终小于0， 这表明系统未发生图灵分岔； 当a_w=\mathrm{7.65}时， \mathrm{m}\mathrm{a}\mathrm{x}(\mathrm{R}\mathrm{e}(\lambda ))=\mathrm{0}， 系统出现图灵分岔点； 当a_w>\mathrm{7.65}时， 系统会发生图灵分岔。

当a_w=\mathrm{3}且a_v\ne \mathrm{0}时， 由图 1（b） 可以看出， 随着a_v的增大， \mathrm{m}\mathrm{a}\mathrm{x}(\mathrm{R}\mathrm{e}(\lambda ))随之增大， 系统开始具备发生图灵分岔的条件， 当a_v=\mathrm{28.7}时， 系统出现图灵分岔点， 即当a_v>\mathrm{28.7}时， 系统发生图灵分岔。对比图 1（a） 和图 1（b） 中的红线可以发现， 交叉对流的存在使系统发生图灵分岔所需的a_w从7.65降到3， 交叉对流项的存在降低了水沿缓坡的流动速度。

接下来在交叉对流项a_v\frac{\partial w}{\partial x}存在的情况下， 对局部放牧项\frac{m_{\mathrm{1}}{v}^{\mathrm{2}}}{\mathrm{1}+v^{\mathrm{2}}}是否存在， 即m_{\mathrm{1}}=\mathrm{0}和m_{\mathrm{1}}\ne \mathrm{0}两种情况进行分析。图 1（b） 中局部放牧项存在(m_{\mathrm{1}}=\mathrm{1.5})， 由图 1（b） 可知， 当a_w=\mathrm{3}，a_v>\mathrm{28.7}时， 系统发生图灵分岔； 图 1（c） 中局部放牧项不存在(m_{\mathrm{1}}=\mathrm{0})， 由图1（c）可知， 当a_w=14.5， a_v>\mathrm{13.1}时， 系统发生图灵分岔。对比图 1（b） 和图 1（c）， 局部放牧项的存在使得图灵分岔所需的a_w从14.5降到3， 降低了水沿缓坡的流动速度。

进一步考虑x和y方向扰动存在差异的情况， 即l_{\mathrm{1}}\ne l_{\mathrm{2}}， 特征值实部的最大值在图 2 中给出， 其中图 2（a）， 图 2（c） 和图 2（e） 是特征值实部最大值在不同参数取值下随扰动模数k变化的立体图， 图2（b）， 图 2（d） 和图 2（f） 是特征值实部最大值在不同参数取值下随扰动模数k变化的平面图。

对比图 2（a）、 （b） 和图 2（c）、 （d）， 当a_w=\mathrm{15}， m_{\mathrm{1}}=1.5， 即局部放牧项存在时， 交叉对流项a_v\frac{\partial w}{\partial x}中a_v的增加稍微增大了\mathrm{m}\mathrm{a}\mathrm{x}(\mathrm{R}\mathrm{e}(\lambda ))的值， 促使斑图更易形成， 这个变化并不是特别明显， 因此需通过斑图形成的数值模拟做进一步验证； 对比图 2（c）、 （d）和图 2（e）、 （f）， 当交叉对流项存在时， 即a_w=15， a_v=20， 随着局部放牧导致的植被损失率m_{\mathrm{1}}的增加， \mathrm{m}\mathrm{a}\mathrm{x}(\mathrm{R}\mathrm{e}(\lambda ))的值变小， 在x和y方向扰动存在差异时， 斑图的空间差异减弱， 振幅变小。

3.2 交叉对流对斑图形成的影响

数值模拟在\mathrm{100}\times \mathrm{100}的网格上进行， 初始值设定为对空间均匀定态解施加非均匀随机扰动， 并采用周期性边界条件， 借助有限差分法进行迭代计算， 且固定m_{\mathrm{1}}=\mathrm{1.5}，a_w=\mathrm{15}。图 3 展示了植被从均匀定态解经扰动后逐渐演变为带状斑图的过程， 涵盖t=\mathrm{0}，t=\mathrm{5}，t=\mathrm{15}和t=\mathrm{100}四个关键时刻的植被生物量空间分布， 将竖直方向设为x方向， 即坡向， 水平方向设为y方向。

当a_v=\mathrm{0}时， 在t=\mathrm{0}时刻， 植被生物量均匀分布， 无明显斑块结构； 在t=\mathrm{5}时开始出现植被斑块， 但振幅较小， 空间分布较为零散； 到t=\mathrm{15}时， 小斑块逐渐聚合形成较大的斑块， 部分呈带状趋势， 但带状结构尚未清晰； 至t=\mathrm{100}时， 清晰的带状植被斑图形成， 斑图形态平缓， 并且垂直于x方向， 但振幅相对较小， 生物量分布差异不显著。

当a_v=\mathrm{20}时， 在t=\mathrm{0}时刻与a_v=\mathrm{0}初始状态类似， 生物量均匀分布； t=\mathrm{5}时， 斑块形成速度更快， 局部聚集现象比a_v=\mathrm{0}时更明显， 振幅开始增大； t=\mathrm{15}时， 斑块聚合速度加快， 振幅远大于a_v=\mathrm{0}同期的振幅， 斑图呈现“破碎-聚合”的动态过渡； t=\mathrm{100}时， 形成清晰且稳定的带状斑图， 振幅显著大于a_v=\mathrm{0}的情况且植被带的宽度更窄。

图 4 进一步展示了在(x,\mathrm{50})坐标方格上（y=\mathrm{50}的直线）的v和w的生物量的变化曲线（绿线代表v， 蓝线代表w）。图 4（a） 中的虚线（t=\mathrm{480}）和实线（t=\mathrm{500}）、 图 4（b） 中的虚线（t=\mathrm{180}）和实线（t=\mathrm{200}）均存在时间差\mathrm{\Delta }t=\mathrm{20}， 表明带状植被带沿x方向移动。对比可知， a_v=\mathrm{20}时植被带移动速度比a_v=\mathrm{15}时有所降低。

4 结论

本文将植物的自然死亡因素和植物的人为死亡因素“局部放牧项”同时引入到Klausmeier模型中， 研究了一类具有交叉对流和局部放牧的植被-水模型， 基于此模型， 得到了图灵斑图存在的充分条件， 并对分岔图和斑图进行了数值模拟。理论结果不仅包含了交叉对流相关的参数， 还加入了放牧的相关参数（如放牧导致的植被损失率）， 需要同时满足交叉对流和放牧的参数条件， 才会出现图灵不稳定性。

交叉对流模拟结果表明， 交叉对流项的存在总会降低水沿着缓坡的流动速度， 即植被具有固水性。当x和y方向扰动相等和交叉对流项存在时， 交叉对流项与局部放牧的共同作用， 降低了水沿缓坡的流动速度； 当x和y方向扰动不相等且交叉对流项存在时， 局部放牧项的存在使得斑图的空间差异减弱， 振幅变小。

斑图形成的模拟结果显示， 植被沿水流方向的对流系数的增大加快了斑图形成的速度并且增大了振幅， 改变了植被空间分布的均匀性， 支持了理论结果， 最后当斑图稳定后， 交叉对流项使得植被带的移动速度减慢。本研究为理解放牧干扰下和交叉对流存在时干旱区植被的动态及制定生态保护策略提供了理论支撑。

参考文献

[1] 蔡丽红.土地沙漠化原因及林业防沙治沙措施［J］. 农业灾害研究， 2023， 13（5）： 188-190.

[2] CAI Lihong. Causes of land desertification and measures of forest desertification prevention and control［J］. Journal of Agricultural Catastrophology， 2023， 13（5）： 188-190. （in Chinese）

[3] 肖绪培， 宋乃平， 王兴， 等. 放牧干扰对荒漠草原土壤和植被的影响［J］. 中国水土保持， 2013（12）： 19-23.

[4] XIAO Xupei， SONG Naiping， WANG Xing， et al. Effects of grazing disturbance to the soil and vegetation of desert grassland［J］. Soil and Water Conservation in China， 2013（12）： 19-23. （in Chinese）

[5] KÉFI S， GÉNIN A， GARCIA-MAYOR A， et al. Self-organization as a mechanism of resilience in dryland ecosystems［J］. Proceedings of the National Academy of Sciences， 2024， 121（6）： e2305153121.

[6] 欧阳颀. 反应扩散系统中的斑图动力学［M］. 上海： 上海科技教育出版社， 2000.

[7] 张红桃， 孙桂全. 植被系统的时空动力学研究进展［J］. 数学理论与应用， 2023， 43（2）： 1-15.

[8] ZHANG Hongtao， SUN Guiquan. Progress in spatiotemporal dynamics of vegetation systems［J］. Mathematical Theory and Applications， 2023， 43（2）： 1-15. （in Chinese）

[9] KLAUSMEIER C A. Regular and irregular patterns in semiarid vegetation［J］. Science， 1999， 284（5421）： 1826-1828.

[10] GILAD E， VON HARDENBERG J， PROVENZALE A， et al. A mathematical model of plants as ecosystem engineers［J］. Journal of Theoretical Biology， 2007， 244（4）： 680-691.

[11] HILLERISLAMBERS R， RIETKERK M， VAN DEN BOSCH F， et al. Vegetation pattern formation in semi-arid grazing systems［J］. Ecology， 2001， 82（1）： 50-61.

[12] SHNERB N M， SARAH P， LAVEE H， et al. Reactive glass and vegetation patterns［J］. Physical Review Letters， 2003， 90（3）： 038101.

[13] VON HARDENBERG J， MERON E， SHACHAK M， et al. Diversity of vegetation patterns and desertification［J］. Physical Review Letters， 2001， 87（19）： 198101.

[14] SHERRATT J A. Using wavelength and slope to infer the historical origin of semiarid vegetation bands［J］. Proceedings of the National Academy of Sciences， 2015， 112（14）： 4202-4207.

[15] CONSOLO G， VALENTI G. Secondary seed dispersal in the Klausmeier model of vegetation for sloped semi-arid environments［J］. Ecological Modelling， 2019， 402： 66-75.

[16] DE LA CRUZ M P C， SANTIESTEBAN D A， ÁLVAREZ L M M， et al. On a generalized Klausmeier model［J］. Mathematical Biosciences and Engineering， 2023， 20（9）： 16447-16470.

[17] CONSOLO G， CURRO C， VALENTI G. Turing vegetation patterns in a generalized hyperbolic Klausmeier model［J］. Mathematical Methods in the Applied Sciences， 2020， 43（18）： 10474-10489.

[18] GUO G H， QIN Q， CAO H， et al. Pattern formation of a spatial vegetation system with cross-diffusion and nonlocal delay［J］. Chaos， Solitons & Fractals， 2024， 181： 114622.

[19] GUO G H， ZHAO S， PANG D， et al. Stability and cross-diffusion-driven instability for a water-vegetation model with the infiltration feedback effect［J］. Zeitschrift für Angewandte Mathematik und Physik， 2024， 75（2）： 33.

[20] SACO P M， WILLGOOSE G R， HANCOCK G R. Eco-geomorphology and vegetation patterns in arid and semi-arid regions［J］. Hydrology and Earth System Sciences Discussions， 2006， 3（4）： 2559-2593.

[21] 张鹏. 基于改进Klausmeier模型的交叉对流对斑图形成的影响研究［J］. 应用数学进展， 2022， 11（6）： 3590-3598.

[22] ZHANG Peng. Research on the effects of cross-advection on pattern formation based on modified Klausmeier model［J］. Advances in Applied Mathematics， 2022， 11（6）： 3590-3598. （in Chinese）

[23] SIERO E. Nonlocal grazing in patterned ecosystems［J］. Journal of Theoretical Biology， 2018， 436： 64-71.

[24] SIERO E， SITEUR K， DOELMAN A， et al. Grazing away the resilience of patterned ecosystems［J］. The American Naturalist， 2019， 193（3）： 472-480.

[25] SIERO E. Resolving soil and surface water flux as drivers of pattern formation in Turing models of dryland vegetation： A unified approach［J］. Physica D： Nonlinear Phenomena， 2020， 414： 132695.

基金资助

陕西省自然科学基金资助项目(2021JQ-662)

国家自然科学基金资助项目(11626182)