一类自扩散–交叉扩散离散捕食系统的斑图形成

The Pattern Formation of a Certain Kind of the Self- and Cross-Diffusion Discrete Predator-Prey System

摘要: 通过构建三链耦合映像格子模型，研究了Holling-II型时空离散反应自扩散–交叉扩散捕食系统的时空动力学。在分析并确定该时空离散捕食系统的均匀稳定状态后，对其进行图灵失稳分析，计算出了图灵失稳条件。在图灵失稳基础上进行的数值模拟，展现了空间密度分布斑图的自组织过程，发现捕食者与被捕食者空间密度斑图总会由无规律的随机分布演化出块状、圆盘状以及螺旋波状斑图。

Abstract: The spatiotemporal dynamics of Holling-II spatiotemporal discrete-time reactive self-diffusion and cross-diffusion predator-prey system was studied by constructing a three-chain coupled map lattice model. After analyzing and determining the homogeneous stable state of the spatio-temporal discrete predator-prey system, the Turing instability was analyzed and the Turing instability conditions were calculated. The numerical simulation on the basis of Turing instability showed the self-organizing process of spatial density distribution patterns, and it was found that the spatial density patterns of predator and prey always evolved into block, disk and spiral wave patterns from random distributions.

1. 引言

在捕食系统中，两个不同种群之间的互动方式为捕食，而捕食行为最主要的特征就是捕食者对被捕食者种群有直接的影响。在捕食系统中，捕食者捕捉被捕食者，同时被捕食者奋力地从捕食者的捕捉中逃脱。自Lotka和Volterra第一次对捕食系统进行研究以来，研究者们为了探讨捕食系统中的非线性行为，建立了许多动力学模型 [1] [2]，这极大地推进了种群动力学的发展。在捕食系统的理论模型研究中，必须考虑生态关系的空间组成和种群间的相互作用。基于耦合映像格子(CML)对捕食系统斑图动力学进行研究，有助于了解捕食系统中多模式复杂斑图的产生机理 [3]。

在反应自扩散–交叉扩散捕食系统中，存在着种群的自扩散过程，交叉扩散过程以及种间反应过程。前人的研究表明，交叉扩散能够产生稳定的捕食者和被捕食者斑图 [4]。但是，在多过程的情形下，尤其是多过程的发生时间和空间不匹配的情形下，双链CML可能不太适用。本研究在双链CML的基础上，根据反应自扩散–交叉扩散捕食系统中特有的生态学过程，拓展出三链CML以进行斑图研究。利用三链CML模型对种群自扩散–交叉扩散的空间运动形式进行研究具有创新意义，有助于进一步理解捕食系统中的自组织现象。

2. 三链CML的捕食模型构建

具有Holling-II型功能反应函数的捕食者–被捕食者模型可以表示如下

$\frac{\text{d}U}{\text{d}T}=\beta U\left(1-\frac{U}{K}\right)-\frac{EMUV}{A+U}$ (1a)

$\frac{\text{d}V}{\text{d}T}=V\left(\frac{MU}{A+U}-D\right)$ (1b)

其中U，V分别表示被捕食者与捕食者密度；T 表示时间； $\beta$，K ，E ，D ，M ，A 是正常数； $\beta$ 为被捕食者的内禀增长率；K为被捕食者的环境容纳量；E为被捕食者的相对损失；D为捕食者的死亡率。

为方便起见，令 $t=\beta T$， $u=U$， $v=EMV/\beta$，则方程(1)可表示为

$\frac{\text{d}u}{\text{d}t}=u\left(1-\frac{u}{k}\right)-\frac{uv}{a+u}$ (2a)

$\frac{\text{d}v}{\text{d}t}=v\left(\frac{bu}{a+u}-c\right)$ (2b)

其中 $a=A$， $k=K$， $b=M/\beta$， $c=D/\beta$。

反应扩散系统的研究表明，生物种群的空间扩散除了包含自扩散外，同时还可能受到其他生物种群空间运动的影响，即还存在交叉扩散 [5]。Bie等人提出了一个带强耦合交叉扩散的反应自扩散–交叉扩散捕食模型 [6]

$\frac{\partial u}{\partial t}-\Delta \left(d_{1}u+{\rho }_{12}uv\right)=u\left(1-\frac{u}{k}\right)-\frac{uv}{a+u}$ (3a)

$\frac{\partial v}{\partial t}-\Delta \left(d_{2}v+{\rho }_{21}uv\right)=v\left(\frac{bu}{a+u}-c\right)$ (3b)

$\frac{\partial u}{\partial t}=\frac{\partial v}{\partial t}=0$ (3c)

$u\left(x,0\right)=u_0\left(x\right)$， $v\left(x,0\right)=v_0\left(x\right)$ (3d)

其中d₁，d₂分别表示被捕食者与捕食者的自扩散系数；式(3c)表示捕食者与被捕食者在边界处的通量为0；非负常数ρ₁₂ 和ρ₂₁ 是交叉扩散常数，ρ₁₂表示被捕食者远离捕食者的趋势，而ρ₂₁表示捕食者远离庞大的被捕食者群体的趋势。在连续反应自扩散–交叉扩散模型中，交叉扩散会破坏均衡，产生不稳定性，当随机的自扩散不能产生稳定斑图时，交叉扩散可以产生稳定斑图。

上述反应自扩散-交叉扩散捕食系统是时空连续的，根据前人研究 [7] [8]，可以通过时间、空间离散化，将此时空连续模型发展为耦合映像格子模型。将二维空间划分为n × n 的网格，并对这些网格进行编号，使除边界外的每一个格子都与其周围4个格子相连。同时，将时间离散成一系列连续的时间段。用 $u_{\left(i,j,t\right)}$ 和 $v_{\left(i,j,t\right)}$ ( $i,j\in \left\{1,2,3,\cdots ,n\right\}$， $t\in Z^{+}$ ) 表示在时间t 内第(i, j) 网格中被捕食者与捕食者的密度。因为耦合映像格子的动力学包括扩散阶段和反应阶段，而在反应自扩散–交叉扩散捕食系统中，捕食者和被捕食者种群的空间扩散又包括两个部分，即自扩散和强耦合的交叉扩散 [9]。根据种群反应扩散中的自然生态学过程，耦合映像格子模型将反应自扩散–交叉扩散捕食系统的动力学分离成三段独立过程：种群自扩散，种群交叉扩散以及“反应”过程。首先对(3)中的自扩散空间项进行时空离散，则有

${u^{\prime }}_{\left(i,j,t\right)}=u_{\left(i,j,t\right)}+\frac{\tau }{h^2}{d}_1{\nabla }_d^2{u}_{\left(i,j,t\right)}$ (4a)

${v^{\prime }}_{\left(i,j,t\right)}=v_{\left(i,j,t\right)}+\frac{\tau }{h^2}{d}_2{\nabla }_d^2{v}_{\left(i,j,t\right)}$ (4b)

τ表示时间步长，h 表示空间步长， ${\nabla }_d^2$ 为离散拉普拉斯算子，且

${\nabla }_d^2{x}_{\left(i,j,t\right)}=x_{\left(i+1,j,t\right)}+x_{\left(i-1,j,t\right)}+x_{\left(i,j+1,t\right)}+x_{\left(i,j-1,t\right)}-4x_{\left(i,j,t\right)}$ (5)

其中x表示任意具备空间扩散的变量。

其次考虑捕食系统中两个种群的交叉扩散过程，对(3)的交叉扩散空间项进行时空离散，并令w = uv，得到

${u^{″}}_{\left(i,j,t\right)}={u^{\prime }}_{\left(i,j,t\right)}+\frac{\tau }{h^2}{\rho }_{12}{\nabla }_d^2{w^{\prime }}_{\left(i,j,t\right)}$ (6a)

${v^{″}}_{\left(i,j,t\right)}={v^{\prime }}_{\left(i,j,t\right)}+\frac{\tau }{h^2}{\rho }_{21}{\nabla }_d^2{w^{\prime }}_{\left(i,j,t\right)}$ (6b)

最后考虑捕食系统的反应阶段，对(3)中的非空间项进行时空离散，则反应阶段可表示为

$u_{\left(i,j,t+1\right)}={u^{″}}_{\left(i,j,t\right)}+\tau {u^{″}}_{\left(i,j,t\right)}\left(1-\frac{{u^{″}}_{\left(i,j,t\right)}}{k}\right)-\frac{\tau {u^{″}}_{\left(i,j,t\right)}{v^{″}}_{\left(i,j,t\right)}}{a+{u^{″}}_{\left(i,j,t\right)}}$ (7a)

$v_{\left(i,j,t+1\right)}={v^{″}}_{\left(i,j,t\right)}+\tau {v^{″}}_{\left(i,j,t\right)}\left(\frac{b{u^{″}}_{\left(i,j,t\right)}}{a+{u^{″}}_{\left(i,j,t\right)}}-c\right)$ (7b)

通过对此三链耦合映像格子模型进行分析，推导出在强交叉扩散影响下的捕食系统斑图形成条件，并进行相应的数值模拟。

3. 捕食系统的斑图形成条件

3.1. 不动点和稳定性

除了一些众所周知的英文缩写，如IP、CPU、FDA，所有的英文缩写在文中第一次出现时都应该给出其全称。文章标题中尽量避免使用生僻的英文缩写。

捕食系统的动力学分为空间同质性动态和空间异质性动态，其中空间同质性动态要求

${\nabla }_d^2{u}_{\left(i,j,t\right)}\equiv 0$， ${\nabla }_d^2{v}_{\left(i,j,t\right)}\equiv 0$， ${\nabla }_d^2{w}_{\left(i,j,t\right)}\equiv 0$ (8)

将式(8)代入式(4)~(7)，得到捕食系统的差分方程，写成映射形式，则可表示为

$\left(\begin{array}{c}u\\ v\end{array}\right)=\left(\begin{array}{c}u+\tau u\left(1-\frac{u}{k}\right)-\frac{\tau uv}{a+u}\\ v+\tau v\left(\frac{bu}{a+u}-c\right)\end{array}\right)$ (9)

根据不动点的定义求解映射(9)，可得到反应自扩散–交叉扩散捕食系统的三个不动点 $\left(u_0,v_0\right)$ ： $\left(0,0\right)$； $\left(u_1,v_1\right)$ ： $\left(k,0\right)$； $\left(u_2,v_2\right)$ ： $\left(\frac{ac}{b-c},\frac{1}{k}\left(a+\frac{ac}{b-c}\right)\left(k-\frac{ac}{b-c}\right)\right)$。 $\left(u_0,v_0\right)$ 表示被捕食者和捕食者都灭绝的情况； $\left(u_1,v_1\right)$ 表示捕食者灭绝，而被捕食者达到最大环境容纳量的情况； $\left(u_2,v_2\right)$ 代表捕食者和被捕食者的共存态。

不动点的稳定性可由映射(9)的Jacobian矩阵确定，其在任意点 $\left(u,v\right)$ 的Jacobian矩阵为

$J\left(u,v\right)=\left(\begin{array}{cc}1+\tau \left(1-\frac{2u}{k}-\frac{av}{{\left(a+u\right)}^2}\right)& -\frac{\tau u}{a+u}\\ \frac{\tau abv}{{\left(a+u\right)}^2}& 1+\tau \left(\frac{bu}{a+u}-c\right)\end{array}\right)$ (10)

由于捕食者和被捕食者都灭绝的情况对于捕食系统来说没有现实的生态学意义，所以本文不考虑不动点 $\left(u_0,v_0\right)$ 的稳定性。现在将不动点 $\left(u_1,v_1\right)$ 与 $\left(u_2,v_2\right)$ 分别代入Jacobian矩阵(10)，并计算其特征值。Nayfeh等的研究表明 [10]，如果某一不动点所对应的Jacobian矩阵的两个特征值的模全都小于1，则该不动点是稳定的；否则，若任意一个特征值的模大于1，那么该不动点就是不稳定的。不动点 $\left(u_1,v_1\right)$ 的Jacobian矩阵为

$J\left(u_1,v_1\right)=\left(\begin{array}{cc}1-\tau & -\frac{\tau k}{a+k}\\ 0& 1+\tau \left(\frac{bk}{a+k}-c\right)\end{array}\right)$ (11)

通过计算可得，矩阵 $J\left(u_1,v_1\right)$ 的两个特征值为 ${\lambda }_1=1-\tau$ 和 ${\lambda }_2=1+\tau \left(\frac{bk}{a+k}-c\right)$。可知当 $b>\frac{c\left(k+a\right)}{k}$ 时， ${\lambda }_2>1$，即不动点 $\left(u_1,v_1\right)$ 是不稳定的。

对于不动点 $\left(u_2,v_2\right)$，为方便起见，令 $u^{*}=\frac{ac}{b-c}$，则不动点 $\left(u_2,v_2\right)$ 的Jacobian矩阵为

$J\left(u_2,v_2\right)=\left(\begin{array}{cc}1+\frac{\tau }{k}\left(k-2u^{*}-\frac{a\left(k-u^{*}\right)}{a+u^{*}}\right)& -\frac{\tau u^{*}}{a+u^{*}}\\ \frac{\tau ab\left(k-u^{*}\right)}{k\left(a+u^{*}\right)}& 1+\tau \left(\frac{b{u}^{*}}{a+u^{*}}-c\right)\end{array}\right)$ (12)

通过计算可得，矩阵 $J\left(u_2,v_2\right)$ 的两个特征值为 ${\lambda }_{1,2}=\frac{p\pm \sqrt{p^2-4q}}{2}$，其中 $p=\text{tr}\left(J\left(u_2,v_2\right)\right)$，

$q=\det \left(J\left(u_2,v_2\right)\right)$。通过计算 $\left|{\lambda }_{1,2}\right|<1$ 可得到不动点 $\left(u_2,v_2\right)$ 稳定的条件为 $-1-q<p<1+q$ 和 $q<1$。

3.2. 图灵失稳条件

在没有空间种群运动时，捕食系统处于空间均匀稳定状态。然而在三链CML中，自扩散和交叉扩散的组合效应可能会破坏这种稳定状态，触发图灵不稳定性，导致捕食系统发生空间斑图自组织。为了进行图灵失稳分析，引入空间非均匀扰动在均匀稳定状态 $\left(u_2,v_2\right)$ 对捕食系统进行线性化。设

$u_{\left(i,j,t\right)}=u_2+{\tilde{u}}_{\left(i,j,t\right)}$， $v_{\left(i,j,t\right)}=v_2+{\tilde{v}}_{\left(i,j,t\right)}$ (13)

当 ${\tilde{u}}_{\left(i,j,t\right)}\ll u_2$， ${\tilde{v}}_{\left(i,j,t\right)}\ll v_2$ 时，将式(13)直接代入CML方程，得

$\left(\begin{array}{c}{\tilde{u}}_{\left(i,j,t+1\right)}\\ {\tilde{v}}_{\left(i,j,t+1\right)}\end{array}\right)=J\left(u_2,v_2\right)\left(\begin{array}{c}{\tilde{u}}_{\left(i,j,t\right)}+\frac{\tau }{h^2}{d}_1{\nabla }_d^2{\tilde{u}}_{\left(i,j,t\right)}+\frac{\tau }{h^2}{\rho }_{12}\left(c_1{\nabla }_d^2{\tilde{u}}_{\left(i,j,t\right)}+c_2{\nabla }_d^2{\tilde{v}}_{\left(i,j,t\right)}\right)+\frac{{\tau }^2}{h^4}{\rho }_{12}\left(c_1{d}_1{\nabla }_d^2{\nabla }_d^2{\tilde{u}}_{\left(i,j,t\right)}+c_2{d}_2{\nabla }_d^2{\nabla }_d^2{\tilde{v}}_{\left(i,j,t\right)}\right)\\ {\tilde{v}}_{\left(i,j,t\right)}+\frac{\tau }{h^2}{d}_2{\nabla }_d^2{\tilde{v}}_{\left(i,j,t\right)}+\frac{\tau }{h^2}{\rho }_{21}\left(c_1{\nabla }_d^2{\tilde{u}}_{\left(i,j,t\right)}+c_2{\nabla }_d^2{\tilde{v}}_{\left(i,j,t\right)}\right)+\frac{{\tau }^2}{h^4}{\rho }_{21}\left(c_1{d}_1{\nabla }_d^2{\nabla }_d^2{\tilde{u}}_{\left(i,j,t\right)}+c_2{d}_2{\nabla }_d^2{\nabla }_d^2{\tilde{v}}_{\left(i,j,t\right)}\right)\end{array}\right)$ (14)

其中， $c_1={\frac{\partial w^{\prime }}{\partial u}|}_{\left(u_2,v_2\right)}$， $c_2={\frac{\partial w^{\prime }}{\partial v}|}_{\left(u_2,v_2\right)}$。根据文献中描述的方法 [11]，可计算出算子 ${\nabla }_d^2$ 和 ${\nabla }_d^2{\nabla }_d^2$ 的特征值分别为

${\lambda }_{1i\bar{j}}=4\left({\sin }^2\frac{\left(\bar{i}-1\right)\pi }{n}+{\sin }^2\frac{\left(\bar{j}-1\right)\pi }{n}\right)$ (15a)

${\lambda }_{2i\bar{j}}=16\left({\sin }^4\frac{\left(\bar{i}-1\right)\pi }{n}+{\sin }^2\frac{\left(\bar{i}-1\right)\pi }{n}{\sin }^2\frac{\left(\bar{j}-1\right)\pi }{n}\right)$ (15b)

其中 $\bar{i}$ 和 $\bar{j}$ 代表波数。

假设 $X_{\bar{i}\bar{j}}^{ij}$ 是特征值 ${\lambda }_{i\bar{j}}$ 所对应的特征函数，则用 $X_{\bar{i}\bar{j}}^{ij}$ 分别乘式(14)的两端，再对所有的i和j求和得到

$\begin{array}{c}{\displaystyle \sum X_{\bar{i}\bar{j}}^{ij}}{\tilde{u}}_{\left(i,j,t+1\right)}=a_{11}{\displaystyle \sum X_{\bar{i}\bar{j}}^{ij}}{\tilde{u}}_{\left(i,j,t\right)}+a_{12}{\displaystyle \sum X_{\bar{i}\bar{j}}^{ij}}{\tilde{v}}_{\left(i,j,t\right)}\\ +a_{11}\frac{\tau }{h^2}\left(d_1{\displaystyle \sum X_{\bar{i}\bar{j}}^{ij}}{\nabla }_d^2{\tilde{u}}_{\left(i,j,t\right)}+{\rho }_{12}\left(c_1{\displaystyle \sum X_{\bar{i}\bar{j}}^{ij}}{\nabla }_d^2{\tilde{u}}_{\left(i,j,t\right)}+c_2{\displaystyle \sum X_{\bar{i}\bar{j}}^{ij}}{\nabla }_d^2{\tilde{v}}_{\left(i,j,t\right)}\right)\right)\\ +a_{12}\frac{\tau }{h^2}\left(d_2{\displaystyle \sum X_{\bar{i}\bar{j}}^{ij}}{\nabla }_d^2{\tilde{v}}_{\left(i,j,t\right)}+{\rho }_{21}\left(c_1{\displaystyle \sum X_{\bar{i}\bar{j}}^{ij}}{\nabla }_d^2{\tilde{u}}_{\left(i,j,t\right)}+c_2{\displaystyle \sum X_{\bar{i}\bar{j}}^{ij}}{\nabla }_d^2{\tilde{v}}_{\left(i,j,t\right)}\right)\right)\\ +a_{11}\frac{{\tau }^2}{h^4}{\rho }_{12}\left(c_1{d}_1{\displaystyle \sum X_{\bar{i}\bar{j}}^{ij}{\nabla }_d^2}{\nabla }_d^2{\tilde{u}}_{\left(i,j,t\right)}+c_2{d}_2{\displaystyle \sum X_{\bar{i}\bar{j}}^{ij}{\nabla }_d^2}{\nabla }_d^2{\tilde{v}}_{\left(i,j,t\right)}\right)\\ +a_{12}\frac{{\tau }^2}{h^4}{\rho }_{21}\left(c_1{d}_1{\displaystyle \sum X_{\bar{i}\bar{j}}^{ij}{\nabla }_d^2}{\nabla }_d^2{\tilde{u}}_{\left(i,j,t\right)}+c_2{d}_2{\displaystyle \sum X_{\bar{i}\bar{j}}^{ij}{\nabla }_d^2}{\nabla }_d^2{\tilde{v}}_{\left(i,j,t\right)}\right)\end{array}$ (16a)

$\begin{array}{c}{\displaystyle \sum X_{\bar{i}\bar{j}}^{ij}}{\tilde{v}}_{\left(i,j,t+1\right)}=a_{21}{\displaystyle \sum X_{\bar{i}\bar{j}}^{ij}}{\tilde{u}}_{\left(i,j,t\right)}+a_{22}{\displaystyle \sum X_{\bar{i}\bar{j}}^{ij}}{\tilde{v}}_{\left(i,j,t\right)}\\ +a_{21}\frac{\tau }{h^2}\left(d_1{\displaystyle \sum X_{\bar{i}\bar{j}}^{ij}}{\nabla }_d^2{\tilde{u}}_{\left(i,j,t\right)}+{\rho }_{12}\left(c_1{\displaystyle \sum X_{\bar{i}\bar{j}}^{ij}}{\nabla }_d^2{\tilde{u}}_{\left(i,j,t\right)}+c_2{\displaystyle \sum X_{\bar{i}\bar{j}}^{ij}}{\nabla }_d^2{\tilde{v}}_{\left(i,j,t\right)}\right)\right)\\ +a_{22}\frac{\tau }{h^2}\left(d_2{\displaystyle \sum X_{\bar{i}\bar{j}}^{ij}}{\nabla }_d^2{\tilde{v}}_{\left(i,j,t\right)}+{\rho }_{21}\left(c_1{\displaystyle \sum X_{\bar{i}\bar{j}}^{ij}}{\nabla }_d^2{\tilde{u}}_{\left(i,j,t\right)}+c_2{\displaystyle \sum X_{\bar{i}\bar{j}}^{ij}}{\nabla }_d^2{\tilde{v}}_{\left(i,j,t\right)}\right)\right)\\ +a_{21}\frac{{\tau }^2}{h^4}{\rho }_{12}\left(c_1{d}_1{\displaystyle \sum X_{\bar{i}\bar{j}}^{ij}{\nabla }_d^2}{\nabla }_d^2{\tilde{u}}_{\left(i,j,t\right)}+c_2{d}_2{\displaystyle \sum X_{\bar{i}\bar{j}}^{ij}{\nabla }_d^2}{\nabla }_d^2{\tilde{v}}_{\left(i,j,t\right)}\right)\\ +a_{22}\frac{{\tau }^2}{h^4}{\rho }_{21}\left(c_1{d}_1{\displaystyle \sum X_{\bar{i}\bar{j}}^{ij}{\nabla }_d^2}{\nabla }_d^2{\tilde{u}}_{\left(i,j,t\right)}+c_2{d}_2{\displaystyle \sum X_{\bar{i}\bar{j}}^{ij}{\nabla }_d^2}{\nabla }_d^2{\tilde{v}}_{\left(i,j,t\right)}\right)\end{array}$ (16b)

其中a₁₁，a₁₂，a₂₁，a₂₂是矩阵 $J\left(u_2,v_2\right)$ 的元素。令

${\bar{u}}_t={\displaystyle \sum X_{\bar{i}\bar{j}}^{ij}}{\tilde{u}}_{\left(i,j,t\right)}$， ${\bar{v}}_t={\displaystyle \sum X_{\bar{i}\bar{j}}^{ij}}{\tilde{v}}_{\left(i,j,t\right)}$ (17)

则有

$\left(\begin{array}{c}{\bar{u}}_{t+1}\\ {\bar{v}}_{t+1}\end{array}\right)=\left(\begin{array}{cc}{A}_{11}& A_{12}\\ A_{21}& A_{22}\end{array}\right)\left(\begin{array}{c}{\bar{u}}_t\\ {\bar{v}}_t\end{array}\right)$ (18)

其中

$A_{11}=a_{11}-\frac{\tau }{h^2}{\lambda }_{1i\bar{j}}\left(a_{11}{d}_1+a_{11}{\rho }_{12}{c}_1+a_{12}{\rho }_{21}{c}_1\right)-\frac{{\tau }^2}{h^4}{\lambda }_{2i\bar{j}}{c}_1{d}_1\left(a_{11}{\rho }_{12}+a_{12}{\rho }_{21}\right)$ (19a)

$A_{12}=a_{12}-\frac{\tau }{h^2}{\lambda }_{1i\bar{j}}\left(a_{12}{d}_2+a_{11}{\rho }_{12}{c}_2+a_{12}{\rho }_{21}{c}_2\right)-\frac{{\tau }^2}{h^4}{\lambda }_{2i\bar{j}}{c}_2{d}_2\left(a_{11}{\rho }_{12}+a_{12}{\rho }_{21}\right)$ (19b)

$A_{21}=a_{21}-\frac{\tau }{h^2}{\lambda }_{1i\bar{j}}\left(a_{21}{d}_1+a_{21}{\rho }_{12}{c}_1+a_{22}{\rho }_{21}{c}_1\right)-\frac{{\tau }^2}{h^4}{\lambda }_{2i\bar{j}}{c}_1{d}_1\left(a_{21}{\rho }_{12}+a_{22}{\rho }_{21}\right)$ (19c)

$A_{22}=a_{22}-\frac{\tau }{h^2}{\lambda }_{1i\bar{j}}\left(a_{22}{d}_2+a_{21}{\rho }_{12}{c}_2+a_{22}{\rho }_{21}{c}_2\right)-\frac{{\tau }^2}{h^4}{\lambda }_{2i\bar{j}}{c}_2{d}_2\left(a_{21}{\rho }_{12}+a_{22}{\rho }_{21}\right)$ (19d)

计算(18)的两个特征值，可得到

${\bar{\lambda }}_{\pm }=\frac{1}{2}\left(\left(A_{11}+A_{22}\right)\pm \sqrt{{\left(A_{11}-A_{22}\right)}^2+4A_{12}{A}_{21}}\right)$ (20)

由于(18)描述了所有位置空间异质扰动下捕食系统的动力学方程，所以，若(18)收敛，那么捕食系统将回到空间均匀的稳定共存状态；而当(18)发散时，则捕食系统的稳定状态将会被打破，发生图灵失稳，导致图灵斑图的自组织，形成图灵斑图。(18)的发散条件为 $\left|{\bar{\lambda }}_{+}\right|>1$ 或者 $\left|{\bar{\lambda }}_{-}\right|>1$，也即至少存在一组满足 $\left|{\lambda }_{i\bar{j}}\right|\ne 0$ 的 $\bar{i}$ 和 $\bar{j}$ 使得

$L=\underset{\bar{i},\bar{j}}{\max }\left({\bar{\lambda }}_{+},{\bar{\lambda }}_{-}\right)>1$ (21)

当条件(21)满足时，出现图灵失稳，捕食系统形成空间异质性斑图。

4. 数值模拟

根据以上理论结果，本节通过数值模拟呈现在强交叉扩散影响下的捕食系统空间密度分布斑图的形成。由于捕食者和被捕食者斑图空间结构的互补性，故只展现被捕食者的斑图。选取以下一组参数值进行数值模拟，a = 0.1，b = 0.5，c = 0.3，K = 2，d₁ = 0.1，d₂ = 0.2，ρ₁₂ = 0.02，ρ₂₁ = 0.01，τ = 0.01，h = 1，得到图1中的捕食者块状斑图形成。可以发现，当t = 100 时，空间斑图仍呈随机分布；从t = 800开始，被捕食者空间斑图开始出现密度较大的区域，捕食者在一定区域内密度也会增加；当t = 5000时，这一趋势变得更加明显，被捕食者密度区域面积分布更广，斑点面积增大；当 t = 8000 时，被捕食者种群密度进一步增大，占据大部分空间；当 t = 10000 时，由于被捕食者种群内部竞争以及捕食等原因，被捕食者高密度区域减少至块状斑图状态。

(a) t = 100 (b) t = 800 (c) t = 1000 (d) t = 5000 (e) t = 8000 (f) t = 10000

Figure 1. Pattern formation of block pattern image

图1. 块状斑图随时间的演变图像

(a) t = 100 (b) t = 500 (c) t = 800 (d) t = 1000 (e) t = 5000 (f) t = 10000

Figure 2. Pattern formation of disk pattern image

图2. 圆盘状斑图随时间的演变图像

选取参数a = 7，b = 10，c = 8，K = 41，d₁ = 2，d₂ = 2，ρ₁₂ = 0.7，ρ₂₁ = 0.9，τ = 1，h = 21，模拟得到图2的圆盘状斑图。当t = 100时，被捕食者的空间密度分布斑图呈斑点状；当t = 500~1000时，被捕食者种群开始聚集，高密度区域面积增大；当t = 5000~10000时，斑图呈圆盘分布，高密度区域集中在中央位置。

(a) t = 10 (b) t = 100 (c) t = 1000 (d) t = 2000 (e) t = 8000 (f) t = 10000

Figure 3. Pattern formation of spiral wave pattern image

图3. 螺旋波状斑图随时间的演变图像

当参数条件选择为a = 4，b = 8，c = 6，K = 19，d₁ = 1，d₂ = 1，ρ₁₂ = 0.8，ρ₂₁ = 0.7，τ = 1，h = 10，捕食系统呈现图3所示的螺旋波状斑图。从随机斑图出发，被捕食者斑图逐渐发展成不规则的块状和条状。当t > 2000时，被捕食者的空间密度分布斑图开始具有规律性，出现螺旋状和圆弧状。当t > 8000时，捕食系统稳定地呈现螺旋波状斑图。

NOTES

^*通讯作者。