第37 卷第14 期
2007 年7 月
数学的实践与认识
MA THEMA T ICS IN PRACT ICE AND THEORY
Vo l137　No114　
July, 2007　
确定Lotka-Volterra 生态系统模型高精度参数的研究
唐静波, 　胡智渊, 　张彦琼
指导教师: 任力锋
(中南大学生物医学工程研究院, 湖南长沙　410078)
摘要: 　研究确定Lo tka2Vo lterra 生态系统模型的高精度参数估计问题. 利用周期性, 先对测量数据进行预
处理; 然后用三种不同的方法构造了误差函数, 进行非线性最小二乘法参数估计; 再用计算机仿真对其进行
验证. 结果表明该方法能够有效地解决高精度参数估计中消除测量数据误差的问题.
关键词: 　非线性最小二乘法; 参数估计; 周期性; 平均值法
1　问题提出
收稿日期: 2007203201
　　在“神舟六号”载人航天宇宙飞船、人造地球卫星等航天器围绕地球在轨运行的过程中,
必须建立高精度的数学模型和高精度地估计模型中的大批参数. 因为只有这样的数学模型
才能解决实际问题, 而不会出现差之毫厘, 结果却失之千里的情况. 为此人们要设法利用长
期积累的丰富的观测资料, 高精度确定这些重要的参数.
由于航天器的问题太复杂, 下面仅考虑较简单的确定高精度参数问题.
假设有一个生态系统, 其中含有两种生物, 即: A 生物和B 生物, 其中A 生物是捕食者,
B 生物是被捕食者. 假设t 时刻捕食者A 的数目为x ( t) , 被捕食者B 数目为y ( t) , 它们之间
满足以下变化规律:
x ′( t) = x ( t) [a1 + a2y ( t) ]
y ′( t) = y ( t) [a3 + a4x ( t) ]
　　初始条件为:
x ( t0) = a5
y ( t0) = a6
(1)
其中ak (1 F k F 6) 为模型的待定参数.
该模型为著名的Lo tka2Vo lterra 方程[ 1—2 ]. 通过对此生态系统的观测, 可以得到相关的
观测数据. 为了说明问题, 本文仅以DA TA 1. TXT , DA TA 4. TXT 为例[ 3 ].
问题1. 1　在观测数据无误差的情况下, 至少需要多少组观测数据, 才能确定参数ak (1
F k F 6) ? 有关数据见数据文件:DA TA 1. TXT
问题1. 2　在观测数据x、y、t 均有误差的情况下, 试建立数学模型, 确定参数ak (1 F k
F 6) 在某种意义下的最优解, 并与仿真结果比较. 有关数据见数据文件: DA TA 4. TXT
2　模型求解
Lo tka2Vo lterra 模型的假设:
1) 假设是一个包括捕食与被捕食两个种群的生态系统, 捕食者靠被捕食者而生存, 系
统与外界没有种群交换关系.
2) 假设种群自身的增殖变化速率与t 时刻种群的数量成比率(马尔萨斯定律).
2. 1　方程通解法
对于问题1. 1, 可由Lo tka2Vo lterra 方程求得方程的通解, 然后利用其通解, 构建误差函
数E
1, 由于参数ak 越精确, E
1 的值越接近0. 由此可以转化为对函数M in E
1 的非线性最小值
无约束优化模型. 对于此非线性最小值无约束优化模型的求解, 我们采用非线性最小二乘法
对参数a1～ a4 和参数C 进行估计[ 4—5 ]. 另外根据测量资料DA TA 1. TXT 得知x ( t0) = a5 =
10, y ( t0) = a6 = 60.
Lo tka2Vo lterra 方程通解:
a1 ln y - a3 ln x + a2y - a4x = C (2)
　　误差函数:
E
1 = (a1 ln y - a3 ln x + a2y - a4x - C) 2 (3)
M inE
1 = Σ
N
i= 0
[a1 ln y i ( t) - a3 ln x i ( t) + a2y i ( t) - a4x i ( t) - C ]2 (4)
　　MA TLAB 实现过程:
1) 用Lo tka2Vo lterra 方程通解, 构建误差函数E
1 (3) ,M inE
1 (4).
2) 用非线性最小二乘法对参数a1～ a4 和参数C 进行估计(M at lab fm in search 函数).
3) 用估计得到的参数对Lo tka2Vo lterra 方程进行求解(M at lab ode45 函数).
4) 将仿真结果与原始数据进行比较.
仿真结果与原始数据非常吻合, 见图1. 图中横坐标表示捕食者x ( t) , 纵坐标表示被捕
食者y ( t) ; 原始数据用红色‘3 ’表示, 仿真数据用蓝色实线表示.
我们做了一个统计, 分别使用5 组、4 组、3 组数据进行参数估计的结果见表1. 结果表
示: 至少需要4 组观测数据可以确定参数, 但计算结果不够精确, 建议至少使用5 组数据.
图1　DA TA 1. TXT 数据与仿真数据图2　DA TA 4. TXT 数据与仿真数据
对比图(通解法) 对比图(通解法)
对于问题1. 2, 由于测量资料DA TA 4. TXT 的x ( t) , y ( t) 和时间变量t 均有误差, 故需要对
这些数据进行预处理. 我们首先对原始数据进行插值: ti = 0: 0. 01: 15, 求出ti 时刻对应的x ( t) ,
y ( t) , 以消除时间变量t 的误差. 然后利用周期解定理(见2. 2 的定理1) , 对插值后的数据按周期
T 进行分段, 共9 段(每段长162 个数据点, 共1458 个数据点). 将1 至9 段中对应时间点的x ( t) ,
y ( t) 加起来, 求其平均值x
q
i ( t)、y
q
i ( t) , i = 1, 2, ⋯, 162, 以减少数据x ( t) , y ( t) 的误差.
78 数　学　的　实　践　与　认　识 37 卷
表1　用5 组、4 组、3 组数据对Lo tka2Vo lterra 方程进行参数估计统计表
数据组5 组4 组3 组
模
型
参
数
a1 - 0. 23239389370315 - 0. 25718506931796 - 0. 60526159001990
a2 0. 02323002030168 0. 02581146824179 0. 27384481371427
a3 1. 39566452999218 1. 54218980966463 1. 45825956312549
a4 - 0. 11628571769750 - 0. 12853240689509 - 0. 15877718474521
SS
残差平方和
1. 85807e- 007 2. 05723e- 007
M S= SSön
平均残差平方和
3. 7161400e- 008 5. 143075e- 008
仿真曲线和原始数据
偏差很大, 故不推荐
只使用3 组数据.
利用误差函数E
1 和函数M inE
1, 采用非线性最小二乘法对参数a1～ a4 和参数C 进行估
计, 参数a5, a6 可由测量资料DA TA 4. TXT 中直接得到. 仿真结果见图2.
x i ( t) =
Σ
9
j= 1
x ( j - 1) ×162+ i ( t)
9
, 　y i ( t) =
Σ
9
j= 1
y ( j - 1) ×162+ i ( t)
9 　( i = 1, 2, ⋯, 162) (5)
2. 2　平均值法
定理1　设x ( t) , y ( t) 是Lo tka2Vo lterra 方程(1) 的一个周期解, 具有周期T > 0. 定义x
和y 的平均值为
x
q
=
1
T∫T
0
x ( t) dt, 　y
q
=
1
T∫T
0
y ( t) dt, (6)
　　那么, x
q
= - a3
öa
4 及y
q
= - a1
öa
2, 换句话说, x ( t) 和y ( t) 的平均值就是平衡值(图3).
图3　Lo tka2Vo lterra 模型的轨迹
证明　用x ( t) 去除(1) 的第一个等式的两边, 得到x ′
x
= a1 + a2y , 于是
1
T∫T
0
x ′
x
dt =
1
T∫T
0
[a1 + a2y ]dt (7)
　　令∫T
0
x ′
x
dt = ln x (T ) - ln x (0) , 因为x (T ) = x (0) , 于是这个等式等于零. 因此
14 期 唐静波, 等: 确定Lo tka2Vo lterra 生态系统模型高精度参数的研究 79
1
T∫T
0
a2y ( t) dt =
1
T∫T
0
- a1dt = - a1, (8)
　　于是, y
q
= - a1
öa
2. 类似地, 用y ( t) 去除(1) 的第二个等式的两边并从0 到T 进行积分,
可得x
q
= - a3
öa
4.
对于问题1. 2, 由于测量资料DA TA 4. TXT 的x ( t) , y ( t) 和时间变量t 均有误差, 故需
要对这些数据进行预处理. 我们首先对原始数据进行插值: ti = 0: 0. 01: 15, 求出ti 时刻对
应的x ( t) , y ( t) , 以消除时间变量t 的误差. 然后利用周期解定理(定理1) , 对插值后的数据
按周期T 进行分段(每段长162 个数据点, 共9 段, 共1458 个数据点) , 并分别将1 至9 段内
的x ( t) , y ( t) 的162 个数据点加起来求平均值x
q
i ( t) , y
q
i ( t) , i = 1, 2, ⋯9, 以减少系统变量
x ( t) , y ( t) 的误差.
利用周期解定理(定理1) 中的种群的平均值模型, 构造误差函数E
2 (11) , 并采用非线性
最小二乘法对参数a1～ a4 进行估计, 仿真结果见图4.
x i ( t) =
Σ
162
j= 1
x ( i- 1) ×162+ j ( t)
162
, 　y i ( t) =
Σ
162
j= 1
y ( i- 1) ×162+ j ( t)
162 　( i = 1, 2, ⋯, 9) (9)
　　种群的平均值模型:
x
q
= - a3
öa
4
y
q
= - a1
öa
2
　] a3 + a4x
q
= 0
a1 + a2y
q
= 0
(10)
　　误差函数:
E
2 = [a1 + a2y
q
i ( t) ]2 + [a3 + a4x
q
i ( t) ]2 (11)
M inE
2 = Σ
N
i= 1
{[a1 + a2y
q
i ( t) ]2 + [a3 + a4x
q
i ( t) ]2} (12)
图4　DA TA 4. TXT 数据与仿真
数据对比图(平均值法)
图5　DA TA 4. TXT 数据与仿真数据
对比图(四阶龙格2库塔法)
2. 3　四阶龙格2库塔法
由于已知在DA TA 4. TXT 中测量数据x ( t) , y ( t) , t 均有误差, 故首先对原始数据进行
插值, 求出ti 时刻( ti = 0: 0. 001: 15) 对应的x ( t) , y ( t) , 在这里时间t 的增量, 即步长h 取
01001. 然后根据四阶龙格2库塔法公式, 求出Lo tka2Vo lterra 方程中x ( t) , y ( t) 的增量d x ,
80 数　学　的　实　践　与　认　识 37 卷
d y , 并利用d x , d y 构建误差函数E
3 (14) , 用非线性最小二乘法对参数a1～ a4 进行参数估
计, 以减少系统变量x ( t) , y ( t) 和时间t 的误差. 仿真结果见图5.
d x n= x n+ 1- x n= h
6
(k1+ 2k2+ 2k 3+ k4)
k 1= x n (a1+ a2y n)
k 2= x n+ h
2 a1+ a2 y n+ h
2 k1
k 3= x n+ h
2 a1+ a2 y n+ h
2 k2
k 4= (x n+ h) [a1+ a2 (y n+ hk3) ]
　
dy n= y n+ 1- y n= h
6
(k1+ 2k2+ 2k 3+ k4)
k1= y n (a3+ a4x n)
k2= y n+ h
2 a3+ a4 x n+ h
2 k 1
k3= y n+ h
2 a3+ a4 x n+ h
2 k 2
k4= (y n+ h) [a3+ a4 (x n+ hk 3) ]
(13)
误差函数:
E
3 = [d x n - d x
3
n ]2 + [d y n - d y
3
n ]2 (14)
M inE
3 = Σ
N
i= 0
{[d x n - d x
3
n ]2 + [d y n - d y
3
n ]2} (15)
　　d x n , d y n: 根据(13) 计算得到增量
d x
3
n , d y
3
n : 对DA TA 4. TXT 数据进行插值, 然后计算得到的增量.
3　模型验证
我们对所有模型验证的方式都采用把原始数据与计算机仿真数据进行比较、计算残差
平方和及平均残差平方和的方法. 计算结果显示, 残差平方和都在可接受误差范围内.
4　模型评价和推广
本文根据Lo tka2Vo lterra 模型的周期性, 先对测量数据进行预处理以消除测量数据的
误差, 然后用三种不同的方法构造误差函数进行非线性最小二乘法参数估计. 计算机仿真的
结果表明该方法能够有效地解决高精度参数估计中消除测量数据误差的问题.
参考文献:
[ 1 ]　W ILL IAM E LUCA S. 微分方程模型[M ]. 长沙: 国防科技大学出版社, 1998.
[ 2 ]　徐克学. 生物数学[M ]. 北京: 科学出版社, 2001.
[ 3 ]　h ttp: ööwww. shumo. com 第三届全国研究生数学建模竞赛B 题.
[ 4 ]　刘庆上, 赵捍东, 王芳. 用优化方法求解最小二乘法问题[J ]. 弹箭与制导学报, 2004, 24 (3) : 83—86.
[ 5 ]　张志涌. 精通MA TLA TB6. 5 版[M ], 北京: 北京航空航天大学出版社, 2003.
Study on Determ ination of High Accuracy Parameter in
Lotka-Volterra Artif ic ial Ecosystem Model
TAN G J ing2bo, 　HU Zh i2yuan, 　ZHAN G Yan2qiong
A dviso r: REN L i2feng
( Inst itute of B iomedical Engineering of Cent ral South U niversity, Changsha Hunan 410078, Ch ina)
14 期 唐静波, 等: 确定Lo tka2Vo lterra 生态系统模型高精度参数的研究 81
Abstract: 　To study the est imat ion p roblem on determ inat ion of h igh accuracy parameter in
Lo tka2Vo lterra art ificial eco system model. A t first, the measured datum w ere p ret reated by
periodicity, and then const ruct th ree different erro r funct ion w h ile p rocessing nonlinear least
square method parameter est imat ion; finally, to verify the results by computer simulat ion. In
conclusion, it is effect ive to elim inate erro rs of the measured datum in h igh accuracy parameter
est imat ion by th ismethod.
Keywords: 　nonlinear least squaremethod; parameter est imat ion; periodicity; average number
method
期刊简介
本刊主要刊登数学的最新的理论成果, 及其在工业、农业、环境
保护、军事、教育、科研、经济、金融、管理、决策等工程技术、自然科学
和社会科学中的应用成果、方法和经验. 主要任务是沟通数学工作
者与其他科技工作者之间的联系, 推动应用数学在我国的发展, 为四
化建设作贡献.
主要栏目: 数学建模、管理科学、工程、应用、问题研究、知识与进
展、学科介绍、方法介绍、高等数学园地、数学史、研究简报、书刊评
介、简讯.
注: ① 投稿一式二份, 原稿自己保存, 编辑部不退还投稿的稿件.
② 编辑部不接收网上投稿.
82 数　学　的　实　践　与　认　识 37 卷