本发明是一种基于stella系统动力学构建的模拟在人类投食影响下大天鹅数量变化的方法，特别是一种适用于模拟预测开放投食的保护区的大天鹅数量。
背景技术：
大天鹅是一种冬候鸟，属于雁形目鸭科天鹅属，国家ⅱ级保护野生动物，是湿地的生态健康指示器。在我国，大天鹅的主要越冬地有山东荣成沿海、东营、埕口，黄河中下游湿地、青海湖、新疆等。随着人类活动能力范围的扩大和不断的改变自然环境，人为干扰或已经成为鸟类的主要威胁之一。由于湿地受人为于扰程度的日趋加剧，湿地水鸟正在更大范围内受到严重影响，甚至在以往人迹罕至的南极洲，海鸟的种群数量由于发展快速的旅游和研究者的活动而急剧下降。而在一些保护区内，人类的投食和观赏行为是主要的人为干扰，这两种行为对鸟类数量有影响，鸟类经常把人类作为潜在的捕食者，大量的研究探讨了鸟类对人为干扰的反应，如人为干扰会使鸟类心脏频率加快，增加报警或防御行为，人为干扰还能够增加鸟类被捕食的风险。在人类活动较多的区域，频繁的人为干扰导致鸟类的存活或繁殖成功率下降。因此，为更加深刻、直观的研究鸟类数量与人为干扰之间的作用机理，除了许多实际工作以外还需要构建一些模型来模拟人为干扰对鸟类数量的影响。
关于鸟类数量模型可分为三个模块：鸟类模块、游客模块、食物模块。人为干扰对自然保护区鸟类数量影响模拟分析软件是通过模拟游客投喂食物数量及环境提供食物数量对自然保护区栖息鸟类数量影响，进而研究人为干扰对鸟类数量影响的内在机制，模拟不同程度干扰下对自然保护区鸟类的影响。通过鸟类食物来源和数量变化，研究人为干扰对自然保护区鸟类数量的影响，对自然保护区管理及鸟类保护有着重要意义。
技术实现要素：
本发明的目的是，提供一种基于stella系统动力学软件构建的适用于开放投食的保护区的大天鹅数量的方法，用于模拟在不同人为干扰下保护区越冬大天鹅数量变化。以更加有效地预测大天鹅数量变化趋势，为大天鹅保护提供更加科学合理的参考数据和建议。
为了实现上述技术目的，本发明采用的技术方案是：
一种人类投食影响下大天鹅数量模拟方法，适用于开放投食的天鹅湖保护区，包括如下几个步骤：
s1、通过统计大天鹅所吃的食物种类、以及大天鹅所吃食物种类中各种食物所占比例，构建大天鹅数量模型的食物模块，所述食物模块用于模拟大天鹅的食物来源，对大天鹅整个越冬期不同种类食物数量进行动态模拟；
s2、通过统计大天鹅越冬时期天鹅湖保护区内游客数量变化，及每位游客平均投食的量，构建大天鹅数量模型的游客模块，所述游客模块用于模拟大天鹅越冬期游客数量变化，游客数量受到鸟类数量的影响，对大天鹅整个越冬期游客数量进行动态模拟；
s3、通过统计天鹅湖保护区内大天鹅数量、行为、以及大天鹅所吃食物的总量及比例，构建大天鹅数量模型的鸟类模块，所述鸟类模块用于模拟大天鹅越冬期大天鹅数量变化；
s4、根据天鹅湖保护区内食物、游客数量及投食对大天鹅数量的影响，将食物模块、游客模块和鸟类模块进行耦合，组合完成人类投食影响下大天鹅数量模型；
s5、输入模型初始值，运行模型，得到运行结果，完成人类投食影响下大天鹅数量的模拟。
所述天鹅湖保护区内，大天鹅所吃的食物种类包括自然环境中的大叶藻、周边农田的冬小麦苗和游客投喂的玉米。
所述步骤s2中统计游客数量变化具体是：统计天鹅湖保护区内游客投食的食物种类、投食量，得到天鹅湖保护区内大天鹅越冬时期游客数量变化。
所述步骤s3中统计天鹅湖保护区内大天鹅行为包括大天鹅飞入和飞离的数量。
通过stella软件对所述人类投食下大天鹅数量模型进行构件，其中，
人类投食下大天鹅数量模型中鸟类模块公式为：
x(t)＝x(t–dt) (xin–xout)·dt，其中，
x(t)为天鹅湖保护区内大天鹅数量(只)，是鸟类模块随时间变化的状态变量；x为天鹅湖大天鹅数量(只)；
xin为每周飞入天鹅湖保护区的大天鹅数量(只/周)，通过以下公式获得：
xin＝rxin xinli–x·0.0001，其中，
rxin为大天鹅每周飞来的增长速率(只/周)；xinli为食物限制大天鹅飞入因素，通过以下公式获得：xinli＝0.0000001·(ywheat·90 yzostera·8 ycorn·2)，该公式是鸟类模块与食物模块耦合公式，以下公式为食物模块食物种类及数量变化公式，其中，ywheat为大天鹅食物之冬小麦量(kg)；
yzostera为大天鹅食物之大叶藻量(kg/周)；
ycorn为大天鹅食物之玉米量(kg/周)；
xout为每周飞离天鹅湖保护区的大天鹅数量(只/周)，通过以下公式获得：
xout＝rxout x·0.0001 xoutli
其中，rxout为大天鹅每周飞离的增长速率(只/周)；
xoutli为食物限制大天鹅飞离因素，通过以下公式获得：
xoutli＝(ywheat·0.5 yzostera·0.3 ycorn·0.2)/1000000
ywheat(t)为大天鹅食物之冬小麦量(kg)，是食物模块随时间变化的状态变量，通过以下公式获得，
ywheat(t)＝ywheat(t–dt) (yin_wheat–yout_wheat)·dt，其中，
yin_wheat为大天鹅食物之冬小麦输入情况(kg/周)；
yout_wheat为大天鹅食物之冬小麦消耗情况(kg/周)；
yzostera(t)为大天鹅食物之大叶藻量(kg/周)，是食物模块随时间变化的状态变量，通过以下公式获得，
yzostera(t)＝yzostera(t–dt) (yin_zostera–yout_zostera)·dt，其中，
yin_zostera为大天鹅食物之大叶藻量输入情况(kg/周)；
yout_zostera为大天鹅食物之大叶藻量消耗情况(kg/周)；
ycorn(t)为大天鹅食物之玉米量(kg/周)，是食物模块随时间变化的状态变量，通过以下公式获得，
ycorn(t)＝ycorn(t–dt) (yin_corn–yout_corn)·dt，其中，
yin_corn为大天鹅食物之玉米输入情况(kg/周)；
yout_corn为大天鹅食物之玉米消耗情况(kg/周)；
ycorn＝b
yin_wheat＝ryin
ryin为大天鹅食物增长率(kg/周)，通过以下公式获得，
ryin＝b/1000，其中，
b为游客喂食的玉米量(g)，通过以下公式获得，该公式为食物模块与游客模块耦合公式，
b＝a·z，其中，
a为每位游客投食的玉米量(g/人)；z为天鹅湖景区游客数量(人)；
yout_wheat＝ryout·0.9
ryout为大天鹅食物消耗率(kg/周)；
yout_zostera＝ryout·0.05
yout_corn＝ryout·0.05
ryout＝x·2·7
z(t)为天鹅湖景区游客数量(人)，是食物模块随时间变化的状态变量，通过以下游客模块游客数量模拟公式获得，
z(t)＝z(t–dt) (zin–zout)·dt，其中，
zin为天鹅湖景区游客输入量(人/周)，通过以下公式获得，
zin＝rzin x·0.001，其中，
rzin为天鹅湖景区游客数量增长率(人/周)；
zout为天鹅湖景区游客输出量(人/周)，通过以下公式获得，
zout＝rzout z/400，其中，
rzout为天鹅湖景区游客数量减少率(人/周)。
步骤s5具体是：输入模型初始值，包括食物、游客数量和大天鹅数量数据，在stella软件界面设置运行时长及步长，本模型中运行时长为大天鹅整个越冬期10月至次年4月，共24周，步长为0.25周，再运行模型，得到运行结果，完成人类投食影响下大天鹅数量的模拟。
本发明的有益效果：
1.通过大量数据收集工作和不同模块及模块间生态作用的耦合，可以更好的综合生态系统中多个要素的作用，更加有效地预测大天鹅数量变化趋势，为大天鹅保护提供更加科学合理的参考数据和建议。
2.本发明针对开放投食的大天鹅保护区大天鹅和游客数量进行分析、模拟及预测，对开放投食的保护区大天鹅数量预测更有针对性，对鸟类生态系统研究有着长期且深远的影响。
附图说明
图1——大天鹅数量模型流程图；
图2——大天鹅数量模型框架示意图；
图3——大天鹅数量模型高层结构图；
图4——大天鹅数量模型在stella软件中结构图；
图5——大天鹅数量模型模拟结果图。
具体实施方式
下面结合具体实施例以及说明书附图对本发明的技术方案作进一步说明。
如图1和如图3所示，本发明有三部分组成：食物模块、游客模块和鸟类模块。食物模块主要模拟大天鹅在保护区越冬时期三种食物来源、数量及变化。游客模块模拟在大天鹅越冬期保护区内游客数量及投食数量变化。鸟类模块模拟在游客和食物影响下，整个越冬期的大天鹅数量变化过程，模块输出鸟类数量变化数据。
一种通过构建人类投食影响下大天鹅数量模型来模拟大天鹅数量的方法，如图2示意图所示，其特征步骤如下：
(1)通过研究荣成大天鹅食性得到大天鹅主要食物来源为：冬小麦苗、大叶藻和玉米，在stella软件中设置3种食物来源及大天鹅食物占比，构建大天鹅数量模型的食物模块；
(2)通过观察游客投食行为和游客数量估算，得到大天鹅越冬时期(10月至次年4月)荣成天鹅湖游客数量变化，及每位游客平均投食玉米的量，在stella软件中构建大天鹅数量模型的游客模块；
(3)通过研究荣成天鹅湖数量、行为等，大天鹅数量受食物总量影响，其中包括自然环境中的大叶藻、周边农田的冬小麦苗和游客投喂的玉米，在stella软件中构建大天鹅数量模型的游客模块；
(4)根据食物、游客数量及投食对大天鹅数量的影响，将食物模块、游客模块和鸟类模块耦合，组合完成人类投食影响下大天鹅数量模型，如图4所示；
所述食物模块模拟大天鹅主要食物来源，包括：自然环境中的大叶藻、周边农田的冬小麦苗和游客投喂的玉米，对大天鹅整个越冬期(10月至次年4月)不同种类食物数量进行动态模拟，影响鸟类模块食物供给；
所述游客模块模拟大天鹅越冬期(10月至次年4月)游客数量变化，游客数量受到鸟类数量的影响，鸟类数量越多，吸引的游客越多，游客数量越多，投喂的食物(玉米)越多，根据这个作用关系，对大天鹅整个越冬期(10月至次年4月)游客数量进行动态模拟，影响鸟类模块和食物模块；
所述鸟类模块模拟大天鹅越冬期(10月至次年4月)大天鹅数量变化，鸟类数量变化受到食物供给以及游客数量的影响，游客数量增加后，投食的玉米量大，可以满足更多大天鹅的实物需求，但同时游客干扰会破坏大天鹅对安全生境的需求；大天鹅数量增加可以吸引更多游客，同时对自然环境大叶藻和周边农田的冬小麦造成压力；当食物数量下降后，大天鹅会去寻找其他食物来源甚至改变越冬栖息地；越冬期大天鹅数量受到多种因素影响，需要耦合三个模块进行综合的模拟分析。
人类投食下大天鹅数量模型中相关参数赋值为：
x＝0
ywheat＝4786238.75
yzostera＝224000
yin_wheat＝0
yin_zostera＝0
z＝1000
a＝100
rxout＝graph(time)(0.00,0.00),(4.80,0.00),(9.60,50.0),(14.4,100),(19.2,200),(24.0,400)
rzin＝graph(time)(0.00,500),(4.80,600),(9.60,600),(14.4,400),(19.2,250),(24.0,200)
rzout＝graph(time)(0.00,0.00),(4.80,50.0),(9.60,150),(14.4,700),(19.2,700),(24.0,800)
rxin＝graph(time)(0.00,200),(4.80,150),(9.60,70.0),(14.4,0.00),(19.2,0.00),(24.0,0.00)
x为天鹅湖大天鹅数量(只)；
ywheat为大天鹅食物之冬小麦量(kg)，食物的种类和数量是通过stella模型中矩阵进行赋值，wheat代表冬小麦苗量；
yzostera为大天鹅食物之大叶藻量(kg)，食物的种类和数量是通过stella模型中矩阵进行赋值，zoster代表大叶藻量；
yin_wheat为大天鹅食物之冬小麦量输入情况(kg/周)，根据食物数量和种类的矩阵赋值，食物输入变化也按照不同种类进行输入，wheat代表冬小麦苗量；
yin_zostera为大天鹅食物之大叶藻量输入情况(kg/周)，根据食物数量和种类的矩阵赋值，食物输入变化也按照不同种类进行输入，zostera代表大叶藻量；
z为天鹅湖景区游客数量(人)；
a为每位游客投食的玉米量(g/人)；
rxout为大天鹅每周飞离的增长速率(只/周)，该参数是通过stella软件中图形输入法进行赋值，以(4.80,0.00)为例，4.80代表时间(week)，0.00代表游客数量增长率；
rzin为天鹅湖景区游客数量增长率(人/周)，该参数是通过stella软件中图形输入法进行赋值，以(4.80,600)为例，4.80代表时间(week)，600代表飞离增长率；
rzout为天鹅湖景区游客数量减少率(人/周)，该参数是通过stella软件中图形输入法进行赋值，以(4.80,50.0)为例，4.80代表时间(week)，50.0代表游客数量减少率；
rxin为大天鹅每周飞来的增长速率(只/周)，该参数是通过stella软件中图形输入法进行赋值，以(4.80,150)为例，4.80代表时间(week)，150代表飞来的增长速率；
人类投食下大天鹅数量模型中相关公式为：
x(t)＝x(t–dt) (xin–xout)·dt
xin＝rxin xinli–x·0.0001
xout＝rxout x·0.0001 xoutli
ywheat(t)＝ywheat(t–dt) (yin_wheat–yout_wheat)·dt
yzostera(t)＝yzostera(t–dt) (yin_zostera–yout_zostera)·dt
ycorn(t)＝ycorn(t–dt) (yin_corn–yout_corn)·dt
ycorn＝b
yin_wheat＝yin
yout_wheat＝yout·0.9
yout_zostera＝yout·0.05
yout_corn＝yout·0.05
z(t)＝z(t–dt) (zin–zout)·dt
zin＝rzin x·0.001
zout＝rzout z/400
b＝a·z
xinli＝0.0000001·(ywheat·90 yzostera·8 ycorn·2)
xoutli＝(ywheat·0.5 yzostera·0.3 ycorn·0.2)/1000000
ryin＝b/1000
ryout＝x·2·7
x(t)为天鹅湖大天鹅数量(只)，是鸟类模块随时间变化的状态变量；
xin为每周飞来天鹅湖的大天鹅数量(只/周)；
xout为每周飞离天鹅湖的大天鹅数量(只/周)；
ywheat(t)为大天鹅食物之冬小麦量(kg)，是食物模块随时间变化的状态变量；
yzostera(t)为大天鹅食物之大叶藻量(kg/周)，是食物模块随时间变化的状态变量；
ycorn(t)为大天鹅食物之玉米量(kg/周)，是食物模块随时间变化的状态变量；
b为游客喂食的玉米量(g)；
yin_wheat为大天鹅食物之大玉米输入情况(kg/周)，根据食物数量和种类的矩阵赋值，食物输入变化也按照不同种类进行输入，corn代表游客投食玉米量；
ryin为大天鹅食物增长率(kg/周)；
yout_wheat为大天鹅食物之冬小麦量输出情况(kg/周)，根据食物数量和种类的矩阵赋值，食物输入变化也按照不同种类进行输入，wheat代表冬小麦苗量；
ryout为大天鹅食物消耗率(kg/周)；
z(t)为天鹅湖景区游客数量(人)，是食物模块随时间变化的状态变量；
zin为天鹅湖景区游客输入量(人/周)；
zout为天鹅湖景区游客输出量(人/周)；
xinli为食物限制大天鹅飞入因素；
xoutli为食物限制大天鹅飞离因素；
(5)输入模型初始值，包括食物、游客数量和大天鹅数量等数据，在stella软件界面设置运行时长及步长，本模型中运行时长为大天鹅整个越冬期10月至次年4月，共24周(week)，步长为0.25周，再运行模型，得到运行结果，完成人类投食影响下大天鹅数量的模拟。
模拟实例：输入2016年山东荣成大天鹅保护区的食物数据，冬小麦苗和大叶藻初始质量分别为4786238.75kg、224000kg和，大天鹅初始数量0只；模拟结果如下表；
参数值为相关文献和实地数据测定所得，假设大天鹅可觅食小麦范围为整个成山镇(天鹅湖所在区域)，成山镇小麦种植面积为1765公顷，换算单位后为271.175g·m–2，所以成山镇可供天鹅觅食的冬小麦的量为4786238.75kg。其中，种植面积是根据2007年威海统计年鉴中荣成市夏粮种植面积换算而来，换算方法为按照成山镇占荣成市的面积比换算种植面积。经过查询，天鹅湖面积为5000km2，假设湖中大叶藻生长状况一致，大叶藻冬季生物量为44.8g·m–2，可供大天鹅觅食的大叶藻质量为224000kg。大天鹅食物中的玉米主要来源于参观游客投食，根据观察，每个游客投食量在100g左右，随着人数增加，投食量增加。模拟结果如图4所示，经过实测数据验证，准确度高。
根据图5可以看到大天鹅数量呈现先增加后减少的趋势，分别是越冬前期数量增加，越冬中期数量达到最高，越冬后期数量逐渐减少。
以上所述，仅是本发明的较佳实例而已，并非对本发明作其他形式的限制，任何熟悉本专业的技术人员可能利用以上揭示的技术内容加以变更或改型为等同变化的等效实施例。但是凡是未脱离本发明技术方案内容，依据本发明的技术实质对以上实施例所作的任何简单修改、等同变化与改型，仍属于本发明技术方案的保护范围。
1.一种人类投食影响下大天鹅数量模拟方法，适用于开放投食的天鹅湖保护区，其特征在于，包括如下几个步骤：
s1、通过统计大天鹅所吃的食物种类、以及大天鹅所吃食物种类中各种食物所占比例，构建大天鹅数量模型的食物模块，所述食物模块用于模拟大天鹅的食物来源，对大天鹅整个越冬期不同种类食物数量进行动态模拟；
s2、通过统计大天鹅越冬时期天鹅湖保护区内游客数量变化，及每位游客平均投食的量，构建大天鹅数量模型的游客模块，所述游客模块用于模拟大天鹅越冬期游客数量变化，游客数量受到鸟类数量的影响，对大天鹅整个越冬期游客数量进行动态模拟；
s3、通过统计天鹅湖保护区内大天鹅数量、行为、以及大天鹅所吃食物的总量及比例，构建大天鹅数量模型的鸟类模块，所述鸟类模块用于模拟大天鹅越冬期大天鹅数量变化；
s4、根据天鹅湖保护区内食物、游客数量及投食对大天鹅数量的影响，将食物模块、游客模块和鸟类模块进行耦合，组合完成人类投食影响下大天鹅数量模型；
s5、输入模型初始值，运行模型，得到运行结果，完成人类投食影响下大天鹅数量的模拟。
2.根据权利要求1所述的人类投食影响下大天鹅数量模拟方法，其特征在于，所述天鹅湖保护区内，大天鹅所吃的食物种类包括自然环境中的大叶藻、周边农田的冬小麦苗和游客投喂的玉米。
3.根据权利要求1所述的人类投食影响下大天鹅数量模拟方法，其特征在于，所述步骤s2中统计游客数量变化具体是：统计天鹅湖保护区内游客投食的食物种类、投食量，得到天鹅湖保护区内大天鹅越冬时期游客数量变化。
4.根据权利要求1所述的人类投食影响下大天鹅数量模拟方法，其特征在于，所述步骤s3中统计天鹅湖保护区内大天鹅行为包括大天鹅飞入和飞离的数量。
5.根据权利要求1所述的人类投食影响下大天鹅数量模拟方法，其特征在于，通过stella软件对所述人类投食下大天鹅数量模型进行构件，其中，
人类投食下大天鹅数量模型中鸟类模块公式为：
x(t)＝x(t–dt) (xin–xout)·dt，其中，
x(t)为天鹅湖保护区内大天鹅数量(只)，是鸟类模块随时间变化的状态变量；x为天鹅湖大天鹅数量(只)；
xin为每周飞入天鹅湖保护区的大天鹅数量(只/周)，通过以下公式获得：
xin＝rxin xinli–x·0.0001，其中，
rxin为大天鹅每周飞来的增长速率(只/周)；xinli为食物限制大天鹅飞入因素，通过以下公式获得：xinli＝0.0000001·(ywheat·90 yzostera·8 ycorn·2)，该公式是鸟类模块与食物模块耦合公式，以下公式为食物模块食物种类及数量变化公式，其中，ywheat为大天鹅食物之冬小麦量(kg)；
yzostera为大天鹅食物之大叶藻量(kg/周)；
ycorn为大天鹅食物之玉米量(kg/周)；
xout为每周飞离天鹅湖保护区的大天鹅数量(只/周)，通过以下公式获得：
xout＝rxout x·0.0001 xoutli
其中，rxout为大天鹅每周飞离的增长速率(只/周)；
xoutli为食物限制大天鹅飞离因素，通过以下公式获得：
xoutli＝(ywheat·0.5 yzostera·0.3 ycorn·0.2)/1000000
ywheat(t)为大天鹅食物之冬小麦量(kg)，是食物模块随时间变化的状态变量，
通过以下公式获得，
ywheat(t)＝ywheat(t–dt) (yin_wheat–yout_wheat)·dt，其中，
yin_wheat为大天鹅食物之冬小麦输入情况(kg/周)；
yout_wheat为大天鹅食物之冬小麦消耗情况(kg/周)；
yzostera(t)为大天鹅食物之大叶藻量(kg/周)，是食物模块随时间变化的状态变量，通过以下公式获得，
yzostera(t)＝yzostera(t–dt) (yin_zostera–yout_zostera)·dt，其中，
yin_zostera为大天鹅食物之大叶藻量输入情况(kg/周)；
yout_zostera为大天鹅食物之大叶藻量消耗情况(kg/周)；
ycorn(t)为大天鹅食物之玉米量(kg/周)，是食物模块随时间变化的状态变量，通过以下公式获得，
ycorn(t)＝ycorn(t–dt) (yin_corn–yout_corn)·dt，其中，
yin_corn为大天鹅食物之玉米输入情况(kg/周)；
yout_corn为大天鹅食物之玉米消耗情况(kg/周)；
ycorn＝b
yin_wheat＝ryin
ryin为大天鹅食物增长率(kg/周)，通过以下公式获得，
ryin＝b/1000，其中，
b为游客喂食的玉米量(g)，通过以下公式获得，该公式为食物模块与游客模块耦合公式，
b＝a·z，其中，
a为每位游客投食的玉米量(g/人)；z为天鹅湖景区游客数量(人)；
yout_wheat＝ryout·0.9
ryout为大天鹅食物消耗率(kg/周)；
yout_zostera＝ryout·0.05
yout_corn＝ryout·0.05
ryout＝x·2·7
z(t)为天鹅湖景区游客数量(人)，是食物模块随时间变化的状态变量，通过以下游客模块游客数量模拟公式获得，
z(t)＝z(t–dt) (zin–zout)·dt，其中，
zin为天鹅湖景区游客输入量(人/周)，通过以下公式获得，
zin＝rzin x·0.001，其中，
rzin为天鹅湖景区游客数量增长率(人/周)；
zout为天鹅湖景区游客输出量(人/周)，通过以下公式获得，
zout＝rzout z/400，其中，
rzout为天鹅湖景区游客数量减少率(人/周)。
6.根据权利要求5所述的人类投食影响下大天鹅数量模拟方法，其特征在于，步骤s5具体是：输入模型初始值，包括食物、游客数量和大天鹅数量数据，在stella软件界面设置运行时长及步长，本模型中运行时长为大天鹅整个越冬期10月至次年4月，共24周，步长为0.25周，再运行模型，得到运行结果，完成人类投食影响下大天鹅数量的模拟。
技术总结