全国大学生数学建模竞赛全国一等奖论文设计

2012高教社杯全国大学生数学建模竞赛

承诺书

我们仔细阅读了中国大学生数学建模竞赛的竞赛规则.

我们完全明白，在竞赛开始后参赛队员不能以任何方式（包括电话、电子邮件、网

上咨询等）与队外的任何人（包括指导教师）研究、讨论与赛题有关的问题。

我们知道，抄袭别人的成果是违反竞赛规则的, 如果引用别人的成果或其他公开的

资料（包括网上查到的资料），必须按照规定的参考文献的表述方式在正文引用处和参

考文献中明确列出。

我们郑重承诺，严格遵守竞赛规则，以保证竞赛的公正、公平性。如有违反竞赛规

则的行为，我们将受到严肃处理。

我们参赛选择的题号是（从A/B/C/D中选择一项填写）：

我们的参赛报名号为（如果赛区设置报名号的话）：

所属学校（请填写完整的全名）：

参赛队员 (打印并签名) ：1.

2.

3.

指导教师或指导教师组负责人 (打印并签名)：

日期：年月日赛区评阅编号（由赛区组委会评阅前进行编号）：

2012高教社杯全国大学生数学建模竞赛

编号专用页

赛区评阅编号（由赛区组委会评阅前进行编号）：

全国统一编号（由赛区组委会送交全国前编号）：全国评阅编号（由全国组委会评阅前进行编号）：

甲型H1N1流感的预测、控制和影响模型

摘要

甲型H1N1流感是全国乃至全球人们最受关注的传染病，它的传播速度快，对人们的身体健康危害极大。本文根据香港甲流疫情数据进行分析，对其传播的预测与控制进行研究并建出模型，并提出模型建立的关键和困难以及对卫生部门所采取的预防措施作出评定估计。

针对问题一，为了了解甲流的传播情况，先作出已确诊的病例散点图。根据散点图的情况，分别建立了马尔萨斯模型：

()t e t x 0175.08.1107=，阻滞增长模型：

()t e i t i λ-???

? ??-+=

11110，SIS 模型：

?????

?

---=)11(σλi i dt di ，SIR 模型： ()()?????????==-=-=0000s s i i s d ds N s d d i t i i t

i

λμλ， 以及SIR 模型的改进模型：???

??????

???

???-+=-+==-+=-=βεωωωβωωωεββεωβωωβ)()()(g s p dt d qi g dt di qi dt dr g g g s p gt dg s p dt ds

. 从SIR 模型的改进模型中，可以得出控制传染源、切断传播途径、保护易感人群、隔离

等措施进行预防和控制H1N1甲流的传播。

针对问题二，考虑H1N1对旅游经济的影响，对近几年香港接待海外游客的数据进行拟合，得出2009年后三个月的游客数目18.1984y3 , 26.7907y2, 25.5199y1===，进

而建立灰色预测模型:()()

()???????-=+-=+=--∧

2046820468))1((17669.2500.0124 0124.0)0(11e a b e a b x k x x dt

dx a ，并对其模型

进行了残差检验和关联度检验，从而较为准确的预测出2010的旅客人数为274.9568万人。

【关键词】 H1N1流感 马尔萨斯模型 Logistic 模型 SIR 模型 灰色预测法

一、问题重述

2009年3月底至4月中旬，由墨西哥、美国等地相继发生甲型H1N1流感(A/H1N1 influenza)疫情逐步迅速地蔓延到世界各地。甲型H1N1流感（简称甲流）是一种新型甲型流感病毒引起的急性呼吸道传染病。去年爆发期间全球数千万人染病，死亡人数超过16000人。截至去年12月21日，我国内地确诊110590例，死亡442人。由于甲流的传播速度快，对人们的身体健康危害大，因此得到世界卫生组织的重视和人们广泛的关注。

附件1是香港流感疫情的模拟数据；附件2是香港接待海外旅游人数的模拟数据。收集和阅读有关甲流的相关数据及文章，建立数学模型，解决如下问题：

问题一：对甲流的传播数学模型进行分析，特别地说明怎样才能建立一个真正能够预测以及能为预防和控制提供可靠、足够的信息的模型，这样做的困难在哪里？同时，对于卫生部门所采取的措施做出评论，如：提前或延后5天采取严格的隔离措施，对疫情传播所造成的影响做出估计（附件1提供的数据可供参考）。

问题二：收集甲流对经济某个方面影响的数据，建立相应的数学模型并进行预测（附件2提供的数据可供参考）。

二、问题分析

根据附件1香港疫情数据分析，我们初步观察到在对65天甲流传播情况包含了对已确诊病例、疑似病例、死亡人数累计量以及治愈出院人数累计量。依据这些数据，首先我们对香港疫情中的已确诊病例情况做出定量分析，运用Mtlab7.1编程得出了甲流传播速度情况的散点图。针对传染病的传播过程，首先，我们用()t x 表示时刻t 的病人人数，用λ表示每天每个病人有效接触的人数，考虑t 到t t ?+时刻病人人数的增加，建

立微分方程 x dt

dx

λ=，()00x x =，通过马尔萨斯模型求解得：()t e x t x λ0=。

接着在病人的有效接触人群中只有病人方可被传染为病人，因此要区分健康人和病人。那么我们再次对这些数据进行分析， 用常数λ表示日接触率；()t s 表示健康者；()t i 表示病人；用()t Ni 表示病人数。那么由此可知每天共有()()t i t Ns λ个健康者被感染。建立模型 Nsi dt di

N

λ=，()()1=+t i t s ，通过阻滞增长模型求解得：()??

???????? ??-+=-t e i t i λ111/10。 接着我们考虑当治愈后的健康者还可被感染变成病人的情况，我们用μ表示日治愈

率，

μ1表示平均传染期，建立模型 Ni Nsi dt

di

N -=λ。 对于问题二，首先我们利用2003年至2008年后7至9月份各个月份的平均值与2009年做差值，利用其差值进行拟合，利用Mtlab7.1求得2003年至2008年与2009年后三个月的差值为2.4468，-2.2407，0.6516，从而得到2009年后三个月香港海外旅游人数。接着同样运用Mtlab7.1编程对2003年到2009年香港海外旅游总人数进行了处理并假设 ()()}{1,196.697,297,326.3,250,292.229.2,217.0=k X ，再对其作一次性累加生成运算得到新的生成数列()()}{,1809 ,1612.3 ,1286.2 989.2 , 696.5 , ,446.5 229.21=k X ，紧接着对()

()k X

1作紧邻均值生成得出数据阵B 和数据向量n Y ，再对参数列T b a ],[=∧

α进行最

小二乘估计最后建立出了灰色模型（GM(1,1)模型）。我们又经过对GM(1,1)模型的残差检验和关联度检验，最终得出了预测结果。

三、符号说明

四、模型假设

1、假设已确诊人数作为主要的预测模型的指标，对于甲流感病情的预测没有影响。

2、假设所有的统计数据真实，没有遗漏现象。

3、假设与患者有效接触的易感染者（即未患过该病的健康者）均会被传染。

4、假设所考查人群的总数恒定，没有其他病源的输入和输出，不考虑总人口的出生率

和自然死亡率。

五、模型的建立与求解

5.1 对问题一建立模型与求解

5.1.1 已确诊病例散点图

根据问题一，由附件1（香港疫情数据）中的已确诊病例数据,用Mtlab7.1作出如下散点图（程序参见附件3）：

图1 散点图

从图1可看出，前25天（即5月20日至6月15日），甲流的传播速度增长幅度较大，而后四十天，甲流的传播速度持续增长，但增长速度趋于平缓。

5.1.2马尔萨斯模型（Malthusian 模型）

甲流传播预测模型类似于人口增长的预测模型，故首先采用马尔萨斯模型（Malthusian 模型）进行建模。设时刻t的病人人数()t x是连续、可微函数，并且每天每个病人有效接触（足以使人致病的接触）的人数为常数λ，考察t到t+t?病人人数的增加，则有

()()()t

=

-

?

+λ

t

x?

t

x

t

t

x

x个病人，即得微分方程

再设0

=

t时有

()?????==0

0x x x

dt

dx

λ 解之可得：

()t e x t x λ0= 其中， λ,0x 为常数。

根据香港疫情数据中的已确诊的病例数据散点图(图1)，考虑利用马尔萨斯模型

()t e x t x λ0=来预测甲流的传播情况。用matlab7.1求得8.11070=x 0175.0=λ。即得

马尔萨斯模型如下（程序参见附件4）： ()t e t x

0175.08.1107= 模型Ⅰ

图2 马尔萨斯拟合及预测图形

结果表明，随着t 的增加，病人人数()t x 无限增长。即马尔萨斯拟合及预测图线与香港疫情中的已确诊病例数据图线拟合程度较差，且对未来预测情况跟实际显然是不太

相符合的，因此暂不考虑用该模型进行数据预测。

对模型Ⅰ的结果分析：马尔萨斯模型是关于人口或种群增长的模型，它发现人口或种群成指数增长。即在模型I 中可引意为，患病人数随着时间得增长呈指数增长变化。但现实生活中，由于病人在有效接触的人群中，包含健康人和病人，而其中只有健康人才可以被传染为病人，因此在改进的模型中必须避免将健康人和病人混为一体这种情况，即要区别病人和健康人进行建模。 5.1.3 阻滞增长模型（Logistic 模型）

在疾病传播期内所考察地区的总人数N 不变，即不考虑生死，也不考虑迁移。人群分为易感染者和已感染者两类，以下简称健康者和病人。时刻t 这两类人在总人数中占得比例分别记作()t s 和()t i 。

假设病人每天的有效接触的平均人数是常数λ，λ成为日接触率。当病人与健康者有效接触时，使健康者受感染变为病人。

根据假设，每个病人每天可使()t s λ个健康者变为病人，因为病人数为()t Ni ，所以每天共有()()t i t Ns λ个健康者被感染，于是Nsi λ就是病人数Ni 的增加率，即有

Nsi dt

di

N

λ= 又因为

()()1=+t i t s

再记初始时刻()0=t 病人的比例为0i ，则

()()?????=-=001i i i i dt di

λ

解之得：

()t

e i t i λ-???

? ??-+=

11110 模型Ⅱ

用Mtlab7.1作出()t t i ~和

i dt

di

~的图形如下（程序参见附件5）

：

图3 Logistic 模型()t t i ~曲线 图4 Logistic 模型i dt

di

~曲线 模型Ⅱ结果分析：由图4可知，当5.0=i 时

dt di 达到最大值m

dt di ???

??，这个时刻为

???

? ??-=-11ln 01i t m λ 此时病人数增加得最快，预示着传染病的高潮的到来。m t 与λ成反比，由于日接触

率λ反应了该地区的卫生水平，λ越小卫生水平越高。所以改善保健设施、提高卫生水平可以延缓传染病高潮的到来。而当∞→t 时1→i ，即所有人终究将被传染，全变为病人，这显然与实际情况不符相。其中的原因是模型中没有考虑到病人是可以治愈的，人群中的健康者只能变成病人，而病人不会再变成健康者。下面模型中将讨论病人可以治愈的情况。

5.1.4 SIS 模型

由于病人被治愈后变成健康者，健康者还可以被感染再变成病人，那么由此得到需增加的条件为：每天被治愈的病人数占病人总数的比例为常数μ，称为日治愈率。病人治愈后成为仍可被感染的健康者。显然

μ

1

是这种传染病的平均传染期。 ?????=--=1)()(t i t s Ni

Nsi dt di

N μλ

记初始时刻()0=t 病人的比例为0i ，则

?????=--=0)0()1(i i i

i i dt di

μλ

设μ

λ

σ=

，则σ可表示整个传染病期内每个病人有效接触的平均人数，称为接触数。利 用σ，可得如下模型：

?????

?

---=)11(σλi i dt di 模型Ⅲ 根据模型Ⅲ，利用Mtlab7.1作出

i dt

di

~的图形，如下（程序参见附件6）

：

图5 SIS 模型的

i dt

di

~曲 图6 SIS 模型的t i ~曲

模型Ⅲ结果分析：不难看出，接触数1=σ是一个阈值。由图5可知道，随着病人所占的人数越多，那么在时间t 内病人的增长率就越大。当1>σ时()t i 的增减性取决于0i 的大小（见图6），单其极限值()σ

1

1-

=∞i 随着σ的增加而增加；当1≤σ时病人比例()

t i 越来越小，最终趋于0，这是由于传染期内经有接触从而使健康者变成的病人数不超过原来病人数的缘故。 5.1.5 SIR 模型

由于病人在治愈后有一定的免疫力，所以病愈的人既非健康者（易感染者），也非病人（已感染者），他们己经退出传染系统。人群分为健康者、病人、病愈与免疫的移

出者三类，即SIR 模型。这三类人在总人数N 中占得比例分别记作i(t)s(t),和r(t)。

???

??==++i t r N d

d N t r t i t s μ1

)()()(

记初始时刻的健康者和病人的比例分别是()000>s s 和()000>i i （设移出者的初始值

00=r ），则可得SIR 模型：

()()?

????????==-=-=0

000s s i i s d ds

N s d d i t

i i t i

λμλ 模型Ⅳ

由于模型Ⅳ无法直接求出)(t s 和)(t i 的值，故先作数值运算。设,02.0)0(,5.0,2===i μλ 96.0)0(=s ，用Mtlab7.1求解可得如下)(),0(t s i 图形和i s ~图形（程序参见附件7）

：

图7 )(),0(t s i 图形 图8 i s ~图形（相轨图）

模型Ⅳ结果分析：s ~ i 平面称为相平面，相轨线在相平面上的定义域为

(){}1,0,0,≤+≥≥=i s i s i s D

消去t d 可得：

00

,11i i s

d d s s s i =-==σ

利用积分特性可解得：()0

00ln

1

s s

s i s i σ

+

-+= ，在定义域D 内,该式表示的曲线即为相轨线,如图9所示.其中箭头表示了随着时间t 的增加)(t s 和)(t i 的变化趋向.

图9 SIR 模型的相轨线

根据图9，可分析)(),(t i t s 和)(t r 的变化情况如下：

1.不论初始条件00,i s 如何,病人消失将消失,即：0=∞i

2.最终未被感染的健康者的比例是s ∞.

3.若σ/10>s ,则)(t i 先增加,后减小.

4.若σ/10≤s ,则)(t i 单调减小至0.

5.1.6 模型Ⅳ（SIR 模型）的改进模型

由于在H1N1流行的过程中，各个地方（包括香港）都采取了一定的措施，一般是采取了隔离的制度，所以在模型模型Ⅳ（SIR 模型）的基础上进行改进。考虑到隔离人数比例g 和未隔离人人数比例w ，以及接触后没有及时隔离治疗的人数p ，从而建立如下改进模型：

???

????

?????

???-+=-+==-+=-=βεωωωβωωωεββεωβωωβ)()()(g s p dt d qi g dt di

qi

dt dr

g g g s p gt dg s p dt ds

模型Ⅴ 由于该模型的分析过程过于复杂，所以该模型在这里将不多做讨论。但从该模型中，

可以看出为预防和控制提供可靠的信息，比如：控制传染源、切断传播途径、保护易感人群、隔离等。

5.1.7 建立模型的关键和困难

建立模型的关键在于对模型进行动态的分析，当传染病发展到一定阶段时，在医疗水平提高、人员流动、出生率和死亡率以及以及采取防御传播措施等方面的影响促使传染率下降。此时仍用之前模型的误差会很大。在建立模型过程中有以下几个方面的困难：①对不同地区H1NI 的卫生知识的宣传程度，K 值取值不同；②对某一地区的不同地方的强化管理也不一样，K 值也就不一样；③防护措施不同、卫生条件不一等，都会影响到K 的取值。另外，本文模型大多假设种群总数为常数，且考虑的影响因素较少，但在实际问题中，由于疾病的复杂性往往涉及变动人口、年龄结构、隔离等多种因素的影响，致使模型的建立错综复杂。

5.1.8 对卫生部门采取的措施评价

经上网查询得知医学研究表明，从正式发病到治愈一般需一至两周，假定平均治愈时间为10天。假设新患者出现的数量与现有患者的数量成正比，也与现有易感者的数量成正比，即发病率是患者人数和易感者人数的双线性函数。则有：

?

?

???-=++-=+∑=t t t Z t M t Z t Z t M t M 0)()()()1()1()()1(ε?ελ 对其进行整理可得：

)())(1()

()

1(ε?ελ∑--=+t

t Z t M t M 模型Ⅵ 其中，)(t M 为t 时刻易感人群总数，)(t Z 为t 时刻新增病人数，)(ε?为病人从患病起经过ε时间仍为病人的概率（图中用p 表示）。

假设病人开始患病记为第1天，最迟到第10天病愈。那么病人从患病起经过ε时间仍为病人呈逐步递减的概率参数，如下图：

图10 概率—时间图

由图10可看出，如果病人发病后5天才开始隔离的话，病人仍患病的概率相当大（图10阴影区域D ），即病人在社会上与易感人群的接触率也相对较大。由模型Ⅴ可得：

1)())(1()

()

1(<--=+∑ε?ελt

t Z t M t M 即)()1(t M t M <+，说明易感人群总数将会以较大的数值递减，给疫情的控制带来更大的困难。所以，如果在病人发病前提前5天隔离的话，新增病人数将变得很小。

5.2 问题二的模型的建立与求解

问题二要我们收集甲流对经济某个方面影响的数据并建立相应的数学模型并进行预测，针对该问题二，我们充分利用附件二，建立甲流对旅游带来的经济影响，而旅游经济与游客数目成正比例关系，故建立预测游客数目模型来预测旅游经济。 5.2.1 香港接待海外旅游人数折线图

根据附件2，利用Excel2003作出2003年至2009年各个月份香港接待海外旅游人数的折线图，如下：

图11 香港接待海外旅游人数折线图

从图11可看出，2003至2008年各整年的海外游客人数的增长率相对稳定，2009年前三月份海外游客人数稳定，从四月份至六月份是因受到H1N1影响而急促下降，从七月份至九月份海外游客又逐步的上升，十月份至十二月份就是要预测的。 5.2.2 2009年后三个月预测

为了预测2009年后3个月的海外旅游人数，根据图11折线变化，利用2003年至2008年后7至9月份各个月份的平均值与2009做差值，利用其差值进行拟合，利用

即4468.21-27.9667=y ，2407.22-24.5500-=y ，6516.03-18.8500=y .所以2009年后三个月香港接待海外旅游人数分别为：18.1984y3 , 26.7907y2, 25.5199y1===（单位：万人）.

5.2.3 灰色预测模型

为了预测2010年香港接待海外旅游总人数，先分别计算出2003至2009年每年的

假设设()()}{1,196.697,297,326.3,250,292.229.2,217.0=k X .

5.2.3.1 GM(1,1)模型的建立

为了使其成为有规律的时间序列数据，对其作一次累加生成运算，即令

)72,1,()

()(1

)0()

1( ==∑=k t k X t X

t

n

从而得到新的生成数列()()}{,1809 ,1612.3 ,1286.2 989.2 , 696.5 , ,446.5 229.2

1=k X .对()()k X 1做紧邻均值生成. 则数据阵B 和数据向量n Y 为

()()()()[]

()()()()[]

()()()()[]

()()()()[]

()()()()[]

()()()()[]

??

??

??

????

??

??

???

???=????????????

??

??

?

?????????????????+-+-+-+-+-+-=1 1710.7-1 1449.3-1 1137.7-1 842.85-1 571.50- 1337.85-1762

116521154211432113221121211

1

1

1

1

1

1

1

1

1

1

1

x x x x x x x x x x x x B ()()()()()()()()()()()()?????????????

??

????

?=????????????????????=196.69326.1297292.7250217.3765432000000x x x x x x Y n 对参数列T b a ],[=∧

α进行最小二乘估计，可得

()

??

????=-0.3494- 0.1574- 0.0716 0.2882 0.4876 0.65930.0005- 0.0003- 0.0001- 0.0001 0.0003

0.00051

T

T B B B ()

?

?

?

???==??

????=-∧

250.7669 0.0124- 1

n T T Y B B

B b a α （其中，a 为发展系数，反映x 的发展趋势；b 为灰色作用量，反映数据间的变化关系. ） 从而可得出GM(1,1)模型：

()

()()???

????-=+-=+=--∧

2046820468))1((17669.2500.0124 0124.0)0(11e a b e a b x k x x dt

dx a 模型Ⅶ

其中，()2046820468))1((10124.0)0(-=+-=+-∧

e a

b

e a b x k x a 为时间响应函数形式。

5.2.3.2 GM(1,1)模型的残差检验

残差大小检验，即对模型值和实际值的残差进行逐点检验.

（1）根据预测公式，计算()

()k X 1∧，得

()

()}

(){6,1,0 91612.3,180.2,1286.2,,696.5,989

,446.5 229.21 ==∧k k X

（2）累减生成()()k X

0∧序列，72,1 =k ()

()}{5603.2059,271.4.8928,26861.6207,26,258.389,2

,255.1972 229.20=∧k X

原始序列： ()()}{1,196.697,297,326.3,250,292.229.2,217.0=k X （3）计算绝对残差和相对残差序列

绝对残差序列：()}{.38071,0.1775,01062,0.108,0.0336,0. 0,0.1744 0=? 相对残差：

}{%19.0%,054431.0%,036397.0%,036283.0%,013440.0%,080258.0,0=φ

GM(1,1)模型的残差检验结果：相对残差不超过0.19%，精确度高。

5.2.3.3GM(1,1)模型关联度检验

关联度是用来定量描述各变化过程之间的差别. 关联系数越大，说明预测值和实际值越接近.

（1）计算序列()

0x 与()

0∧X

的绝对残差序列()0?()k

()

}{.38071,0.1775,01062,0.108,0.0336,0. 0,0.1744 0=?

()(){}{}0.38071,0.1775,01062,0.108,0.0336,0. 0,0.1744m in m in 0==?k

()(){}{}3807.0.38071,0.1775,01062,0.108,0.0336,0. 0,0.1744m ax m ax 0

==?

k

（2）计算关联度

()()()()

)5.0,7,...,1()}

(max{)()}

(max{)}(min{)(0000==?+??+?=P k k P k k P k k η 0.9547 )(11

==∑=n

k i k i n r η

精度检验等级如下表：

精度等级 关联度

好（1级） 0.90≥ 合格（2级） 0.80≥ 勉强（3级） 0.70≥ 不合格（4级）

0.70<

GM(1,1)模型关联度检验结果：关联度为90.09547.0≥，精确度高。 5.2.3.4GM(1,1)模型求解

六、模型评价与推广

模型评价：

在建模前期，全面分析影响甲流疫情的各种因素，找出各因素之间的关系以及作用的时间段和范围，收集比较完整而准确的前期数据。在模型建立中我们采用了各种软件(如Mtlab，Excel等)进行求解，制图精确，计算结果较为准确。但在预测模型中，时间序列数据的时间间隔不是稳定的，这对模型的求解结果的准确性有一定的影响。本文所建立的控制模型忽略了人口流动、变化给该地区甲型H1N1流感带来的影响，从而模型预测结果会与实际情况有一定差距。

模型推广：

通过模型的分析可知，如果全社会的努力和投入的程度继续增加，即隔离措施的提早进行、隔离率增大、防疫药品的早日研发、公众的防御意识提高，可使得疫情周期缩短、患者人数逐步减少。实时监控甲流疫情走势，采集更多的数据以验证模型和改进模型，若有预料之外的干扰因素出现，应及时修正模型，重新预测其后期走势。

七、参考文献

[1] 姜启源谢金星叶俊，数学模型[M]，北京：高等教育出版社，2003（第三版）.

[2] 刘国卫 MATLAB程序设计与应用[M]，北京：高等教育出版社,2011.

[3] SARS传播的数学模型及应用，https://www.360docs.net/doc/599023775.html,/math/SARS123.doc 151K 2008-4-14.

[4] 马知恩周义仓王稳地等. 传染病动力学的数学建模与研究[M].北京，科学出版社，2004.2.

八、附录附件1：香港疫情的数据

附件2：香港接待海外旅游人数（单位：万人）

部分参考文献：

https://www.360docs.net/doc/599023775.html,/s2009/2813/s263659873/

https://www.360docs.net/doc/599023775.html,/html/qikan/dxxb/zsdxxbyxkxb/200910305/wz/20100315 100432420_515523.html

https://www.360docs.net/doc/599023775.html,/Scripts/viewArticle.jsp?articleID=1809092

附件3：

x1=[1:65];

y1=[339 482 588 693 774 877 988 1114 1199 1347 1440 1553 1636 1741 1803 1897 1960 2049 2136 2177 2227 2265 2304 2347 2370 2388 2405 2420 2434 2437 2444 2444 2456 2465 2490 2499 2504 2512 2514 2517 2520 2521 2522 2522 2522 2522 2522 2522 2523 2522 2522 2522 2523 2523 2522 2522 2522 2521 2521 2521 2521 2521 2521 2521 2521];