基于种群动态膜系统的大熊猫种群模型

田 昊，张葛祥*，荣海娜，Mario J．P REZ-JIM NEZ，Luis VALENCIA-CABRERA，陈 鹏，侯 蓉，齐敦武

(1．西南交通大学电气工程学院，成都610031; 2．塞维利亚大学计算机科学与人工智能系，塞维利亚41012，西班牙;3．成都大熊猫繁育研究基地，成都610081; 4．四川省濒危野生动物保护生物学重点实验室，成都610081)

(*通信作者电子邮箱zhgxdylan@126．com)

0 引言

大熊猫是世界自然基金会的形象大使和当之无愧的世界生物多样性保护的旗舰物种［1－6］。在研究人员从事与大熊猫活动相关的工作中，需要对大熊猫的生活习性、种群数量、种群年龄结构组成以及引起种群数量变化的各种因素进行分析和预测［7－9］。这些研究对大熊猫的保护工作具有重要意义［10－11］。

在生态系统种群动态性建模中最广泛使用的方法是常微分方程法和偏微分方程法［12－13］，相关科研人员曾广泛使用莱斯利矩阵对不同大熊猫种群进行动态预测［14－15］。常微分数学方程通过简化的数学分析很容易用简单的公式描述复杂系统，但该方法无法捕捉到空间动态性和随机效应。偏微分方程建模方法虽克服了常微分方程法建模方法中不能捕获生态系统空间动态性的缺点，但建模过程中存在非常复杂的数学分析。随着生态系统复杂性不断增加，采用常微分方程与偏微分方程方法建立的模型会变得更加难以解决问题，以至于这种建模框架在许多情况下不适用。

目前，关于大熊猫种群动态的研究文献还比较少。文献［16－17］利用计算机仿真分析了大熊猫种群生存力;文献［18－19］分析讨论了大熊猫种群动态与它们主要食物来源竹子之间的相互影响作用;文献［20－21］通过计算机仿真研究了栖息地对大熊猫种群动态的影响;文献［22－23］研究了汶川地震对大熊猫种群的影响;在文献［24－25］中，相关研究人员分析了气候变化对大熊猫种群及其保护政策的威胁和挑战。这些研究都无法很好地反映大熊猫种群的空间动态性和随机效应，国外研究人员已经将膜计算数据建模应用于生态系统种群动态性的研究中［26－31］。在国内，本文的实验室团队已经使用概率膜系统对成都大熊猫繁育研究基地的圈养大熊猫生态系统进行建模研究［32－33］，该模型粗略地考虑了大熊猫生命周期中的繁殖、死亡、进食、救护行为，繁殖模块设为每年固定数量的新生大熊猫，并且假定大熊猫为单胞胎繁殖，救护模块中一年实际救护数量最多为1只个体。根据对谱系数据更深入的研究发现大熊猫双胞胎现象较为普遍，并且2005—2008年大熊猫的出生率较之后的年份明显高出很多，查证大熊猫繁殖期后设置合理的大熊猫繁殖年龄，根据处于繁殖期的雌性大熊猫数量作为主要因素之一设计进化规则。救护模块中加入同一年内出现两只个体的可能性，最为重要的是随着野放现象的出现，在计算模型中需要考虑加入野放模块。而且根据概率膜系统的定义，其无法对两个大熊猫种群之间的相互作用与交流进行数据建模，如若后期需要对多个种群进行建模研究，此时概率膜系统则无法实现。

在此提出使用种群动态膜系统对中国大熊猫保护研究中心的圈养大熊猫进行数据建模，新模型从个体行为上来模拟中国大熊猫保护研究中心圈养大熊猫的生态系统，更详细地分析讨论了2005—2016年的大熊猫谱系数据，充分考虑了可能会影响该种群动态性的相关因素。该模型更完整地体现出中国大熊猫保护研究中心圈养大熊猫生态系统的特点，更加符合该生态系统的实际情况，同时适用于接下来对中国大熊猫保护研究中心圈养大熊猫生态系统的研究。种群动态膜系统(Population Dynamics P system，PDP system)比微分方程法更复杂，但是使用的语言更接近专家语言(生态学家)。通俗地说，用PDP系统建模方法通过对象集和一些规则来控制它们的进化和相互影响来模拟大熊猫种群现实生活中的进化过程。在复杂系统建模中对象是最重要的组成部分(大熊猫、空间资源、食物等)，它们同时可以根据一系列的重写规则来描述系统的动态性。每一个相关规则采用概率函数试图捕捉自然中固有的随机性。

1 谱系数据分析处理

本文以中国动物园协会已发布的全球圈养大熊猫谱系数据为依据，从中搜集统计中国大熊猫保护研究中心大熊猫种群数据，并以其作为研究对象。根据对大熊猫谱系数据中年龄结构、繁殖率和死亡率的分析考量，以及与大熊猫专家研讨得到的相关意见，在大熊猫从出生到死亡的整个生命过程中，依照其不同时期的行为特点将其分为幼年期、亚成年期、青年期、中年期、中老年期和老年期共6个年龄阶段，并确定繁殖期所处的年龄范围。以6个年龄阶段为基础，追踪每只大熊猫的行为，统计确定大熊猫在各阶段的食物种类、进食数量、死亡率以及在繁殖期内的大熊猫繁殖率。

同时，由于人为因素的干扰，在此需要考虑野外大熊猫救护行为对大熊猫种群动态性产生影响的因素，根据往年救护大熊猫的数量、年龄和性别的统计和分析，确定每年可能救护的大熊猫数量和年龄分布。野外的雌性大熊猫由于生存环境等原因无法同时抚养双胞胎幼体，之前都是考虑大熊猫是单胞胎的哺乳动物，但在圈养大熊猫生态系统中随着生存环境的改善和人工繁殖技术的发展，双胞胎个体都能得到良好的成长环境，故在此需要把双胞胎的可能性纳入模型。通过与中国大熊猫保护研究中心专家讨论研究发现，不同时间段的大熊猫出生率有较大差别，2005—2008年出生率较高，之后随着大熊猫种群数量的增长、年龄结构的稳定，出生率也随之下降并稳定下来，故本文首次在种群动态膜系统的计算模型中加入分段式的进化规则，对2008年前后两个时期的大熊猫种群参数进行不同的参数设置，使其更符合中国大熊猫保护研究中心圈养大熊猫生态系统的种群特点。圈养大熊猫除了保护种群这个重要原因，更是希望通过野放让大熊猫回归自然，复壮野外小种群。故随着野化技术的成熟，从2012年开始选择合适的大熊猫进行野化放归。

模型中所需的参数如下所示:

i表示圈养大熊猫的性别，其中i=1表示雄性大熊猫，i=2表示雌性大熊猫。

ki，1表示圈养大熊猫年龄达到亚成年阶段;ki，2表示圈养大熊猫年龄达到成年阶段;ki，3表示圈养大熊猫年龄达到中年阶段;ki，4，1表示圈养大熊猫年龄达到中老年阶段;ki，4，2表示圈养大熊猫年龄达到老年阶段;ki，5表示圈养大熊猫的最大年龄;ki，6表示圈养大熊猫处于幼年阶段的死亡率;ki，7表示圈养大熊猫处于亚成年阶段的死亡率;ki，8表示圈养大熊猫处于成年阶段的死亡率;ki，9表示圈养大熊猫处于中年阶段的死亡率;ki，10表示圈养大熊猫处于中老年阶段的死亡率;ki，11表示圈养大熊猫处于老年阶段的死亡率;ki，12表示圈养大熊猫进入繁殖期的年龄;ki，13表示圈养大熊猫结束繁殖期的年龄。

g1表示繁殖期大熊猫生育单胞胎的概率;g2表示繁殖期大熊猫生育双胞胎的概率;g3表示一年中供给圈养大熊猫的竹笋量(kg);g4表示一年中供给圈养大熊猫的竹子量(kg);g5表示一年中供给圈养大熊猫的其他食物(例如水果、牛奶等)量(kg)。

fi，1表示每只幼年大熊猫一年中需要消耗的竹笋总量;fi，2表示每只幼年大熊猫一年中需要消耗的竹子总量;fi，3表示每只幼年大熊猫一年中需要消耗的其他食物总量;fi，4表示每只亚成年大熊猫一年中需要消耗的竹笋总量;fi，5表示每只亚成年大熊猫一年中需要消耗的竹子总量;fi，6表示每只亚成年大熊猫一年中需要消耗的其他食物总量;fi，7表示每只成年大熊猫和老年大熊猫一年中需要消耗的竹笋总量;fi，8表示每只成年大熊猫和老年大熊猫一年中需要消耗的竹子总量;fi，9表示每只成年大熊猫和老年大熊猫一年中需要消耗的其他食物总量。

cmin表示每年救护大熊猫的最少数量;cmax表示每年救护大熊猫的最大数量;cmaxage表示救护大熊猫的最大年龄;pcc表示一年中救护c只大熊猫的概率;pgi表示救护大熊猫性别为i的概率;paj表示救护大熊猫年龄为j的概率;pdd表示一年中野放d只大熊猫的概率;psi表示野放大熊猫性别为i的概率;pyj表示野放大熊猫年龄为j的概率。

2 种群动态膜系统

2．1 种群动态膜系统定义

膜计算是自然计算的新分支，由欧洲科学院院士、罗马尼亚科学院院士Paun提出［34］。膜计算受生物细胞的启发，旨在从细胞组织的结构和功能中抽象出计算模型［35－36］。膜系统(计算模型)由膜结构、初始对象集和进化规则集三要素构成，膜系统中膜结构分隔了对象和规则所处的区域，对象按照所在区域的规则不断进化，进化规则以极大并行的方式执行。膜系统运行开始前的状态称为初始格局，膜系统自初始格局开始根据不同的进化规则从一个格局进化到下一个格局，直到没有规则可以执行时系统终止，此时的膜系统状态称为终止格局，同时终止格局时可以指定膜系统里不同区域中的对象数目作为膜系统的输出结果。多环境概率膜系统是在膜系统的基础上引入概率函数而得到的一种新的膜系统概念，种群动态膜系统是多环境膜系统的一个变种［37］。

一个度为(q，m)，同时满足q，m≥1且T≥1种群动态膜系统可形式化定义为:

其中:

G=(V，S) 是一个 V={e1，e2，…，em} 的有向图。

Γ和Σ是字母表。Γ是包含有限对象的非空字母表，由种群动态膜系统所有区域中的所有对象组成;Σ代表字母表Γ中可以出现在所有环境中的对象组成的字母表。

Τ表示生态系统进化过程被仿真的年数。

RE表示有限规则集，形如

其中:x，y1，y2，…，yh∈ Σ，(ej，ejl) ∈ S，1 ≤l≤ h，1 ≤ k≤ n，且p，p'是定义域为{1，2，…，T}的可计算的函数。

u是集合{1，2，…，q}×{0，+，－}中元素单射标号的有q个节点的根树，如果膜标号为(i，α)那么这个膜将表示为同时该膜的标号为i控制电荷为α，树的根节点设为1作为关联标签。

对于每一个r∈R和1≤j≤m，fr，j是一个定义域为{1，2，…，T}的可计算函数。

对于每一个 i，j(1 ≤ i≤ q，1 ≤ j≤ m)，Mi，j是字母表 Γ中的有限多重集。

对于每一个j，1≤j≤m，Ej是字母表Σ中的有限多重集。

2．2 计算模型设计

整个模型由初始模块、繁殖模块、救护模块、进食模块、死亡模块、野放模块和更新模块6个部分组成。

计算模型考虑使用度为(2，1)时间单元为T的种群动态膜系统模型

定义如下:

G是一个空表。

工作字母表Γ为:

其中:Xi，j，y表示繁殖模块之前性别为i年龄为j的大熊猫个体，Yi，j，y表示进入死亡模块性别为i年龄为j的大熊猫个体，Zi，j，y表示存活性别为 i年龄为 j的大熊猫个体，Wi，j，y表示进食模块之后性别为i年龄为j的大熊猫个体，Ci，j，y表示救护性别为i年龄为j的大熊猫个体，y表该个体目前所处的仿真周期。对象S，B，O表示不同类型的食物:分别为竹笋、竹子和其他食物。对象F，A1，G1是辅助符号分别表示如下:F是一个符号用来作为在每个周期初始产生新的食物，A1是一个符号用来作为触发救护个体的行为;G1是一个符号用来作为触发野放个体的行为，下标1表示目前为第一个仿真周期。

T表示在进化生态系统中仿真的年数。

u=［［］1］2表示采用两层膜结构。

M1，1和M2，1是Γ中的有限多重集，描述区域u中的初始对象:

qi，j表示性别为i年龄为j的初始大熊猫数量

M2，1={F A1G1}

1)初始规则。

为大熊猫提供食物

2)繁殖规则。

未达到繁殖期年龄的个体

3 实验仿真结果

为了验证上述设计的中国大熊猫保护研究中心圈养大熊猫种群动态膜系统计算模型的有效性和正确性，以2005—2016年全球大熊猫谱系数据为依据，分别统计了各年份中国大熊猫保护研究中心圈养大熊猫的数量，用实际数据和实验仿真结果进行对比。把2005年各年龄段的大熊猫数量作为初始参数输入，2006—2016年各年龄段的大熊猫数量作为结果输出。仿真实验所用计算机为戴尔 Inspiron N5110，2．10 GHz，3．16 GB，使用 MeCoSim 仿真软件作为仿真平台。

实验中的参数设置为:仿真周期cycles设置为11，实验中一个仿真周期代表大熊猫生态系统中的一个自然年，周期运行步数 steps为5，为保证实验结果的可靠性仿真次数simulations设置为1000次。

大熊猫生态系统通常涉及参数由于自然变化以及许多因素影响大熊猫个体的健康、死亡率和生物学现象，故所研究的过程本质上是随机的。在种群定投膜系统中，进化规则以概率方式执行，可以有效处理含随机性、不确定性的数据。部分进化规则的关键参数设置如表1～3所示。

表2 大熊猫出生率及其繁殖期年龄Tab．2 Birth ratio and reproductive age of giant panda

表3 野放大熊猫数量及年龄Tab．3 Quantity and age of released giant panda

不同年份雄性大熊猫和雌性大熊猫数量变化的实验仿真结果与实际数据对比分别如图1所示。

图1 雌雄大熊猫数量实验结果与实际数据对比Fig．1 Comparison of experimental and real results with respect to the numbers of female and male giant pandas

不同年份大熊猫总体数量变化实验仿真结果与实际数据进行对比，结果如图2所示。

图2 大熊猫总体数量变化Fig．2 Changes in the overall number of giant pandas

从图1～2的仿真结果可以看出，本文设计的种群动态膜系统模型仿真结果与真实数据相比较，雄性大熊猫仿真结果与实际数据的相对误差最大不超过±8．95%，大部分控制在±2．5%以内;雌性大熊猫仿真结果与实际数据的相对误差最大不超过±7．81%，大部分控制在±2．4%以内;全体大熊猫仿真结果与实际数据的相对误差最大不超过±4．13%，大部分控制在±2．7%以内。同时由于种群基数较小导致相对误差较为明显，但依然远小于±10%的相对误差标准。同时与之前采用概率膜系统设计的简单模型相比，该模型实验结果误差更小，拟合度更高，故本文设计的大熊猫种群动态膜系统较好地模拟了中国大熊猫保护研究中心大熊猫生态系统种群动态变化趋势。

4 结语

本文提出了一种研究大熊猫种群动态的新方法，种群动态膜系统模型能够有效预测种群未来的变化趋势，可以及时评估种群发展的合理性，设计出一个基于大熊猫谱系数据的具有双层膜结构、对象集和一系列进化规则的种群动态膜系统模拟实际生态系统中大熊猫固有的随机进化过程。这是首次搜集和分析如此完整的谱系数据，并以此为基础研究大熊猫的种群动态。同时有专门开发的仿真软件MeCoSim来协助设计该计算模型并仿真验证所提出的计算模型。它还可以进行基于仿真实验结果的决策，诊断种群异常的原因，分析外界因素对种群动态性的影响，为大熊猫基地的工作人员提供决策参考，避免盲目决策。