全 文 ：第 35 卷第 17 期
2015年 9月
生 态 学 报
ACTA ECOLOGICA SINICA
Vol.35,No.17
Sep.,2015
http: / / www.ecologica.cn
基金项目:中央高校基本科研业务费专项资金资助(2014FZA5008)
收稿日期:2013鄄11鄄05; 摇 摇 网络出版日期:2014鄄11鄄03
*通讯作者 Corresponding author.E鄄mail: zhaoxd@ cjlu.edu.cn
DOI: 10.5846 / stxb201311052669
赵晓东, 张宏建, 周洪亮.基于混沌理论的河流藻类生长特性分析———以德国易北河为例.生态学报,2015,35(17):5585鄄5596.
Zhao X D, Zhang H J, Zhou H L.Analysis of the growth characteristics of river algae based on chaos theory: a case study of Elbe River.Acta Ecologica
Sinica,2015,35(17):5585鄄5596.
基于混沌理论的河流藻类生长特性分析
———以德国易北河为例
赵晓东1,2,*, 张宏建1, 周洪亮1
1 工业控制技术国家重点实验室, 浙江大学, 杭州摇 310007
2 中国计量学院计量测试工程学院, 杭州摇 310018
摘要:针对藻类生长具有高度非线性特征和实际采样样本间隔稀疏的问题,采用了混沌理论对采样序列的混沌特征量进行估
计。 采用 C鄄C方法估计时间序列的嵌入维和延迟时间,采用 G鄄P 算法对关联维进行估计,并采用小数据量法估计最大 Lyapunov
指数,最终可实现对最长可预测时长的估计。 以易北河为例,对易北河水体叶绿素 a 1997年至 2001年间各年的观测序列进行
了混沌分析,分析结果表明,各年的叶绿素 a观测序列均具有低维混沌特性,关联维 D= 2.75—4.02,各年的叶绿素 a序列的最长
预测时间变化范围为 8.01—18.94 d,平均为 13.98d(约 2周)。 采用同样方法对 5a易北河连续日径流量时间序列分析表明,该
径流量也具有低维混沌特性(最大 Lyapunov指数 姿1 = 0.0125),径流量的最长预测时间估计约为 80 d。 气候因素的混沌特性对
藻类生长表现出的混沌特征的影响可能要大于径流量等水文因素的影响。
关键词:藻类生长; 混沌; 关联维数; Lyapunov指数
Analysis of the growth characteristics of river algae based on chaos theory: a case
study of Elbe River
ZHAO Xiaodong1,2,*, ZHANG Hongjian1, ZHOU Hongliang1
1 State Key Laboratory of Industrial Control Technology, Zhejiang University, Hangzhou 310007, China
2 College of Metrological and Measurement Engineering, China Jiliang University, Hangzhou 310018, China
Abstract: Algal growth at equilibrium is not sustainable in a river or lake experiencing point or non鄄point source pollution,
which is likely to change with season, location, and human activity. Chlorophyll a, as a common indictor, is an important
reference point for water resource management. It can be affected by the biodiversity of species, as a result of their
oscillations, and chaotic fluctuations. Algal growth has many highly nonlinear characteristics. The characteristics that differ
among species, however, are all susceptible to change under external disturbance. It is difficult to provide a comprehensive
and detailed description of nonlinear algal growth. Variations in chlorophyll a tend to maintain certain regularity; for
example, seasonal variation and the 24 hour cycle, which also display self鄄similarity. However, it is difficult to observe
similar variations at different times from the studied sampling series. These features are similar to aspects of chaotic motion,
such as boundedness, ergodicity, and inherent randomness. The sampling of the indictor chlorophyll a is typically performed
on an hourly, weekly, or even monthly basis. With higher sampling frequency, the sampling series of chlorophyll a in the
field becomes more unstable and appears to be more chaotic. Therefore, this paper aims to study the variations in a
chlorophyll a series sampled from the field, rather than constructing a theoretical model to recapitulate the field data. While
http: / / www.ecologica.cn
there is extensive research on algae in many lakes or rivers, few studies discuss the prediction of algal growth times in
aquatic environments. The aspects of algal growth times are partially addressed in this study.
The variation characteristics of the algal data series were analyzed using chaos theory. The characters of the one鄄
dimensional time series were recovered by reconstruction into the multi鄄dimensional phase space, using phase space
reconstruction. The reconstruction parameters, namely the embedding dimension m and time delay 子, were estimated using
the correlation integral method (C鄄C method). The correlation dimension, D, is the basic mathematical description of the
strange attractor, which is the main characteristic of a chaotic system. D was calculated using the Grassberger鄄Procaccia
algorithm (G鄄P algorithm). Only the largest Lyapunov exponent, 姿1, was estimated through the small data method to
evaluate the diffusion degree of the phase trajectory. The reciprocal of 姿1 is the upper bound of the deterministic prediction
time in the chaotic system, which is designated as the Lyapunov time t0 . This property indicates that the system is
unpredictable beyond t0 . In this paper, the chaotic characteristics of hourly chlorophyll a concentrations and the daily runoff
time series of the Elbe River over a five year period (1997—2001) were analyzed. It was found that both sequences had the
properties of low鄄dimension chaos with 姿1>0, D = 2.75—4.02 for the chlorophyll a series, and D = 1.84 for the runoff
series. The average value of t0(14 days) was estimated for the five鄄year chlorophyll a sequence data. These findings were
remarkably close to the current biggest day鄄to鄄day weather forecast time (two to three weeks) . A much longer value of t0 for
the runoff series, for the same period, was 80 days. This result indicated that compared to the weather factors, the runoff
factor is clearly weaker in affecting the chaotic characteristics of chlorophyll a.
Key Words: algal growth; chaos; correlation dimension; Lyapunov exponent
水体藻种浓度指标实际上是多藻种浓度指标的综合体现,受限于藻种种群间生长与竞争特性。 Huisman
等[1]给出了多藻种在不同营养条件下的数值模拟,藻种生长显示出具有混沌特性,藻种显示出平衡条件下的
多样性。 文中只针对恒定营养盐输入进行了数值模拟,并未给出试验验证结果。 自 1978 年,Connell 等[2]提
出非平衡营养条件下的 “中等干扰假说冶后,对水生环境非平衡营养输入条件下引发的藻种多样性进行了大
量的数值模拟和实验研究[2鄄8]。 国内研究者在对两种常见水华藻类进行 “脉冲冶营养盐输入条件下的实验室
培养的研究过程中,发现具有相同初始生长条件的藻种在非稳态营养输入条件下显示出截然不同的竞争特
性,在生物量上表现出交替占优的特征[9]。 对于实际河流藻类生长观测序列而言,河流下游某观测点获得的
藻类生长数据与上游藻类生长具有相关性。 影响其相关性的因素除沿途期间的气象因素之外,河流径流量也
是主要影响因素之一[10]。 可见,河流中藻类这种生长不确定特征除营养盐输入因素影响之外,还包括水文和
气象条件等因素的影响。 在影响藻类生长的各种因子中,早在 1987 年,Tsonis 和 Elsner 就证明了气候系统是
混沌系统[11],对于河流径流量以及降雨量等水文数据的研究也表明具有混沌特征[12鄄15]。
在对地表水体中藻类生长机理进行描述时,水体网格的划分细密程度,各“藕联冶数据的测量精度都会对
最终藻类预测产生影响,而当该系统具有混沌特性时,上述计算或测量误差都会对藻类长期预测结果产生极
大影响,最终导致不可预测。 另外,由于藻类指标实际采样样本间隔稀疏(理论上应无穷小)以及藻类生长本
身具有高度非线性特征,有必要对藻类实测时间序列是否具有混沌特性进行分析。 目前,尽管各种预测方法
都针对特定湖库和河流能给出了较为成功的应用,但对于究竟能最大预测多长时间国内外并未有文献给出这
方面的报道。 因此,对于河流藻类预测而言,有必要对其观测序列是否混沌进行研究,实现对其可预测时长的
估计。
本文对混沌理论和所采用的经典的混沌分析算法进行了说明,并以易北河为例,对叶绿素 a 观测序列和
日径流量的时间序列进行了混沌特征量估计,以此判断是否该时间序列具有混沌特征。 并在此基础上,对该
时间序列的最大可预测时间进行了估计,并指出了可能成为影响藻类生长混沌特性的主要外部因素。
6855 摇 生摇 态摇 学摇 报摇 摇 摇 35卷摇
http: / / www.ecologica.cn
1摇 分析方法
1.1摇 Takens定理和相空间重构理论
Takens[16]证明了可以找到一个合适的嵌入维,在这个嵌入维空间里可以把有规律的轨迹恢复出来。 该
定理表明,可以将单变量时间序列重构成一个相空间,只要嵌入维数足够高,就可构建出一个拓扑等价的相空
间。 由于混沌系统的策动因素是相互影响的,因而在时间上先后产生的数据点也是相关的,Packard 等[17]提
出了采用原系统中某变量的延迟坐标来重构相空间。 因此重构相空间时,只对一个分量在固定时延点上进行
维向量拓展至三维甚至更高维的空间,就可以重构一个拓扑等价的相空间。
设 x( t),t= 1,2,…,N为观测时间序列。 重构相空间选择嵌入维数 m,时间延迟设为 子,则由 x( t)构造出
一组新的向量序列 X( ti)= {x( ti),x( ti +子),…,x( ti +(m-1) 子)} T,i = 1,2,…,M,M =N-(m-1) 子,如式(1)
所示:
x( t1) x( t2) x( ti) x( tN - (m - 1)子)
x( t1 + 子) x( t2 + 子) x( ti + 子) x( tN - (m - 2)子)
x( t1 + 2子) x( t2 + 2子) x( ti + 2子) x( tN - (m - 3)子)
... ... ... ...
x( t1 + (m - 1)子) x( t2 + (m - 1)子) x( ti + (m - 1)子) x( tN)
X( t1) X( t2) X( ti) X( tM)
(1)
即 X( t)= {X( t1), X( t2),…, X( ti),…, X( tM)}
式(1)表明,时间延迟 子和嵌入维数 m的选择是相空间重构的主要内容,直接影响到相空间重构质量和
其奇异吸引子特征量的描述。
1.2摇 时间延迟的选取
对于无限长、无噪声的时间序列,时间延迟 子的选择原则无限制,但对于有限时间序列,相空间的特征量
则依赖于 子。 Taken[16]认为时间延迟和嵌入维数两者理论上相互独立无关,两者可以单独确定。 也有研究者
认为两者的选择应互相关联[18鄄19],延迟时间的选取应依赖于延迟时间窗口,即 子w =(m-1) 子。
本文主要采用 C鄄C方法进行时间延迟估计,并采用互信息算法计算结果作为对比。 互信息量方法是估
计重构相空间时间延迟的一种有效方法。 该方法利用互信息函数的第一次极小值来确定相空间重构的时间
延迟,Fraser等给出了互信息计算的递归算法[20]。 Kim等人[19]提出了 C鄄C方法,即通过嵌入时间序列的关联
积分构造统计量来代表非线性时间序列的相关性,对时间延迟 子 和时间窗口 子w同时进行估计。 该方法效果
和互信息量法一致,但能有效减少互信息法的计算量,并能保持其非线性特征,且对小数据组可靠,抗干扰能
力较强。
1.3摇 嵌入维数的确定
相空间嵌入维数需要足够大,以便能刻画出该系统的奇异吸引子。 奇异吸引子是描述混沌系统时间序列
的重要特征,其维数一般低于相空间的维数。 关联维数 D常用于描述奇异吸引子的维数,嵌入维数的选取通
常要求 m逸2D+1。 本文采用饱和关联维法(G鄄P 算法) 求取关联维数。
Grassberger和 Proacaccia[21]依据是嵌入定理和相空间重构理论提出的求取关联维的 G鄄P 算法,只需针对
实验观测数据就可以得到吸引子维数,目前采用较多:
对于嵌入维数为 m的一组向量序列维{X i} i =1,2,…,N,其中 X i ={xi 1, xi 2,…, xim}为参考点,分别计算其它
N-1个点与 X i的距离,统计落在以 X i点为中心且距离小于正数 r 的点的个数。 向量之间的距离一般采用两
种方法,2鄄范数和肄 鄄范数。 两种方法的计算距离的公式是拓扑等价的,后者计算量较小[22]。 式(2)所示采用
的是肄 鄄范数计算,即以两个向量的最大分量差作为距离:
椰X i - X j椰 = max1臆k臆m椰xik - x jk椰 (2)
7855摇 17期 摇 摇 摇 赵晓东摇 等:基于混沌理论的河流藻类生长特性分析———以德国易北河为例 摇
http: / / www.ecologica.cn
满足该条件的矢量对数在种配对中所占比例被称为关联积分:
Cm ( )r =
1
N2移i屹j H r - 椰X i - X j( )椰 (3)
其中表示 Heaviside阶跃函数:
( )H x =
0, x 臆 0
1, x >{ 0 (4)
令 dm为 m维空间中的奇异吸引子最大伸展距离,r逸dm时 Cm( r)= N(N-1) / N2 = (N-1) / N,当 N寅肄时,
Cm( r)抑1。 实际上,Cm( r)反映了奇异吸引子中各点之间距离的分布概率,则有:
Cm ( )r =
r
d
æ
è
ç
ö
ø
÷
max
D m,( )r
( r臆dmax) (5)
根据分形理论的标度不变性,当 r位于无标度区间时,D(m,r)是与 m和 r有关的常数,称为关联指数,是
关联维数 D的确定逼近。 取小距离 r1和 r2,则有:
Cm r( )2
Cm r( )1
=
r2
d
æ
è
ç
ö
ø
÷
max
D m,r( )2
/
r1
d
æ
è
ç
ö
ø
÷
max
D m,r( )1
(6)
对式 6两边取对数,则有:
lnCm r( )2 - lnCm r( )1 = ln
r2
d
æ
è
ç
ö
ø
÷
max
D m,r( )2
- ln
r1
d
æ
è
ç
ö
ø
÷
max
D m,r( )1
(7)
当值∣ r2-r1∣很小时,可认为 D(m,r2) 抑 D(m,r1)。 则有:
D m,( )r = lim
r2寅r1
lnCm r( )2 - lnCm r( )1
lnr2 - lnr1
= lim
驻lnr寅0
驻lnCm ( )r
驻lnr
=
dlnCm ( )r
dlnr
(8)
由式(8)可知,D(m,r)是 lnCm( r)—lnr曲线的斜率。 当 r寅0时,可获得关联维数 D的近似值:
D 抑
lnCm ( )r
lnr
(9)
确定关联维的关键是先确定无标度区间。 实际动力系统的分形不同于数学上的分形所具有的在无穷尺
度上的自相似或自仿射性,而是近似地或统计意义上存在自相似,这种自相似仅存在于一定的尺度变化范围,
一旦超出了这个尺度变化范围,其自相似性就不复存在,这个尺度变化范围就是分形的无标度区。 根据标度
区间不变性定理,无标度区间内,Cm( r)与 r 成指数关系,如式(9)。 标量的取值应包含在分形的无标度区之
内,其序列长度不应少于 20,以保证估计到的关联维可靠且满足相关性检验[23]。 近些年,针对无标度区间的
确定方法已展开多方面的研究[24鄄27]。 这些确定的方法有效性和针对性虽然不相同,但其基本步骤均是:1) 改
变嵌入维 m值,获得 lnCm( r)鄄lnr的曲线簇;2) 用最小二乘法拟合每条曲线最佳直线段(无标度区),获得该
直线段的斜率作为嵌入维数 m对应的关联指数;3) 随着 m的增加,关联指数会出现饱和现象,此时最佳直线
段对应的斜率即为该混沌序列对应的关联维数 D。
1.4摇 最大 Lyapunov指数的估计方法
混沌系统可通过其吸引子的宏观特征量表征。 宏观特征量主要包括:关联维,Kolmogorov 熵,Lyapunov 指
数等。 其中,Lyapunov指数的时间数列是根据序列本身所计算出来的客观规律,是量化初始闭轨道的指数发
散和估计方法的混沌量,反映了动力系统整体混沌量水平。 因此,Lyapunov 指数计算对预测混沌时间序列而
言显得尤为重要。
在判断实际动力系统是否为混沌系统时,通常只估计最大 Lyapunov 指数 姿1。 小数据量法是只计算最大
Lyapunov指数的一种方法,是通过对基本轨道上每个点的最近邻近点的平均发散速率估计出最大 Lyapunov
指数。 该算法需要给定嵌入维数,时间延迟以及序列的平均周期。 考虑到小数据量算法对嵌入维数,时间延
迟,数据序列长度以及信噪比的敏感程度较低,具有较好的鲁棒性[28],本文采用了该算法实现对最大
8855 摇 生摇 态摇 学摇 报摇 摇 摇 35卷摇
http: / / www.ecologica.cn
Lyapunov指数的求取。 但为避免计算过程中参数选取的主观性,采用了 C鄄C 方法对嵌入维数和时间延迟同
时估计。 其平均周期的计算是将序列 FFT 变换获得的每个功率值先与其对应的周期加权,再最终求加权平
均计算获得[29鄄30]。 本文只对时间序列的最大 Lyapunov 指数进行求取,该指数是否大于零可用来判别一个时
间序列是否为混沌系统。
1.5摇 Lyapunov时间计算
最大 Lyapunov指数 姿1表示的是当前两相邻轨道距离最大分量的指数发散程度:
啄 ( )x t = 啄x ( )0 e姿1t (10)
式中,啄x(0)表示两相邻轨道初始距离, 啄x( t)表示经过时间后的轨道距离。 混沌运动并不是真正的随机,他
服从确定性的规律,这表明在一定的临界距离变化范围内是可以预测的,而超过这个范围时,轨道运动就不可
预测了[27, 31]。 距离变化如式(11)所示,达到临界值所经历的时间如式(12)所示:
C = 啄
( )x t
啄x ( )0
= e姿1t (11)
t0 =
1
姿1
lnC (12)
通常认为轨迹分离达到原间距的数倍或十几倍时,轨道就不确定了。 此时可认为 lnC抑1,则:
t0 抑
1
姿1
(13)
式(13)所示的时间为 Lyapunov时间,也称为最大预测时间。 当 姿1 越大,预测时间越短,系统运动的可预测性
越差。
2摇 分析结果
本文主要对易北河 1997年—2001年的叶绿素 a观测序列进行了混沌特性分析。 易北河全长 1094km,为
欧洲第四大内陆河,流域面积达 148268 km2。 流域面积主要覆盖捷克和德国,叶绿素 a 的测点位置为
Geesthacht Weir(N53毅25.486忆 E10毅20.665忆),并采用在线叶绿素监测仪(10鄄AU鄄005)每小时采样 1次。 本文采
用 C鄄C算法[19]进行时间延迟 子的确定,采用 G鄄P 算法[21]对数据序列的关联维数进行估计 D,进而确定嵌入
维数 m,并采用 Rosenstein方法[27]计算了叶绿素 a观测序列的最大 Lyapunov指数 姿1,对最大预测时间 t0进行
了估计。
2.1摇 时间延迟估计
叶绿素 a样本选取每年(1997—2001)3月 14日到 9 月 30 日之间的数据序列,每年样本长度均为 4824。
时间延迟如图 1所示。 图 1中 驻軈S( t)表示统计量,P表示互信息量。
每年的叶绿素 a观测数据序列通过 C鄄C方法获得的时间延迟(图 1)。 时间延迟 子 = t子s,子s为采样时间间
隔,取值为 1,则 子= t。 曲线第一个极小值所对应的 t为时间延迟,图中所示时间延迟为 8或 12,该迟延与互信
息法所确定的时间基本一致(图 1)。 互信息法中时间延迟对应曲线第一个极小值。
时间延迟 子比较如表 1所示。 采用互信息方法和 C鄄C方法获得 子差异较小,表明这两种方法在 子的判断
上比较接近。 本文选用 C鄄C方法获得的 子值进行关联维数的计算。
2.2摇 关联维数估计
针对 1997—2001年的叶绿素 a观测序列采取了 C鄄C方法获得时间延迟 子 (表 1),并采用 G鄄P 算法获取
了对应各年份的关联维数 D,如图 2所示,其中标度 r起始值分别参考空间向量距离的最小值和最大值,其步
长设定为 30。 无标度区间的确定方法是采用目前常用的视觉识别方法。
图 2所示,每年的叶绿素 a的观测序列均具有低维分形维。 通过 C鄄C 算法获得的最佳嵌入维数与通过
G鄄P 算法中获得的饱和嵌入维数相接近,如表 2所示。 饱和嵌入维数是指最佳直线段斜率刚达到饱和时所对
应的嵌入维数。 当增加超过饱和嵌入维数时,其直线段斜率变化趋于稳定,由 G鄄P 算法可知,对应的最佳直
9855摇 17期 摇 摇 摇 赵晓东摇 等:基于混沌理论的河流藻类生长特性分析———以德国易北河为例 摇
http: / / www.ecologica.cn
线段斜率值即为关联维数。
图 1摇 叶绿素 a观测序列时间延迟 子 (子= t)
Fig.1摇 The time delay of the sequence of hourly chlorophyll a observation (子= t)
表 1摇 叶绿素 a观测序列时间延迟不同算法计算结果比较*
Table 1摇 Comparison of the time delay of chlorophyll a sequence
年份
Year
互信息方法
Mutual information method
时间延迟 子
Time delay
时间延迟 子
Time delay
C鄄C方法
C鄄C method
延迟时间窗口 子w
Time delay window
嵌入维数 m
Embedding dimension
1997 12 12 545 46
1998 12 8 233 30
1999 14 8 320 41
2000 10 12 537 45
2001 13 12 514 43
摇 摇 *样本长度 N= 4824
表 2摇 各年关联维数和嵌入维数列表
Table 2摇 List of the correlation dimension and embedding dimension in each year
年份 Year 关联维数Correlation dimension
最小嵌入维数
Minimum embedding
dimension
最佳嵌入维数
Optimum embedding
dimension
饱和嵌入维数
Saturated embedding
dimension
1997 4.00—4.02 9 46 50
1998 3.50—3.52 7 30 38
1999 3.10—3.11 7 41 44
2000 2.75—2.76 6 45 44
2001 3.35—3.36 7 43 44
2.3摇 最大 Lyapunov指数估计
选取 C鄄C方法计算所得的嵌入维数 m和时间延迟 子,如表 1,根据 Rosenstein 算法[28],结果如图 3 所示。
图中横坐标表示为离散时间步长,纵坐标 ln
3所示,在 i= 20—100的范围内,曲线表现出较好的直线特性。 选取该区域,采用最小二乘法做出回归直线,
该直线的斜率即为最大 Lyapunov指数 姿1。 Rosenstein算法中根据时间序列的平均周期限制相空间中最近邻
点的短暂分离,时间序列平均周期的算法参考 Rathje等[29鄄30]提出的算法。 每年叶绿素 a浓度时间序列的平均
周期如图 3中 tp值所示。
如表 3所示,各年的叶绿素 a 观测序列的最大 Lyapunov 指数均为非零数,表明其时间序列具有混沌特
0955 摇 生摇 态摇 学摇 报摇 摇 摇 35卷摇
http: / / www.ecologica.cn
性,Lyapunov指数接近于零值,表明该时间序列混沌特性较弱。 5a的叶绿素 a 观测序列最大预测时间最大约
为 19 d,最短约为 8 d,平均约为 14 d (两周)。
表 3摇 不同年份最大 Lyapunov指数和最大预测时间列表
Table 3摇 List of the largest Lyapunov exponent and the maximum prediction time in each year
年份 Year
最大 Lypunov指数 姿1
The largest Lyapunov exponent
最大预测时间 t0 / h
The Lyapunov time
最大预测时间 t0 / d
The Lyapunov time
1997 0.0022 454.55 18.94
1998 0.0052 192.31 8.01
1999 0.0038 263.16 10.96
2000 0.0028 357.14 14.88
2001 0.0025 400.00 16.67
1955摇 17期 摇 摇 摇 赵晓东摇 等:基于混沌理论的河流藻类生长特性分析———以德国易北河为例 摇
http: / / www.ecologica.cn
图 2摇 叶绿素 a观测序列关联维数
Fig.2摇 The correlation dimension D of hourly chlorophyll a observation sequence
图 3摇 不同年份叶绿素 a观测序列 Lyapunov指数(Rosenstein算法)
Fig.3摇 The largest Lyapunov index of the chlorophyll a sequence estimated by the small鄄data method
2.4摇 河流日流量时间序列混沌特性分析
采用与叶绿素 a观测序列同样的分析方法,对易北河 1997年 3月 3日至 2001 年 10 月 31 日之间的日径
2955 摇 生摇 态摇 学摇 报摇 摇 摇 35卷摇
http: / / www.ecologica.cn
流量时间序列进行混沌特性分析,流量采样点位置为 Neu Darcha((N53毅14.207忆 E10毅53.314忆),采样设备为多
普勒测流仪,每日采样一次。 首先采用了 C鄄C方法同时对延迟时间窗口和时间延迟同时估计。 C鄄C 方法中的
最小值对应时间延迟窗口值,的第一个零点或者的第一个极小值对应着时间延迟(子 = t)。 由图 4 可知,tw =
60,子= 17,其中的第一个零点对应值为 子= 60。 通常取较小值作为时间延迟,避免使得非线性模型的拟合关系
复杂,降低难度和减少计算量[24]。 C鄄C方法后的时间延迟与互信息法计算结果相接近,后者 子 = 15,对应曲线
的第一个极小值。
图 4摇 1997—2001年易北河日径流量时间序列时间延迟
Fig.4摇 The time delay of the daily runoff sequence (1997—2001)
1997—2001年易北河日径流量时间序列饱和关联维数为 9,如图 5所示。 而 C鄄C方法获得的 tw和 子,根据
可得 tw =(m-1) 子,m= 5。 G鄄P 算法获得的关联维数 D= 1.84,则日径流量时间序列的嵌入维数 m逸2D+1。 同
样,C鄄C算法获得的嵌入维数也满足此条件。
图 5摇 1997—2001年易北河日径流量时间序列关联维数
Fig.5摇 The correlation dimension of the daily runoff sequence (1997—2001)
根据 Rosenstein算法,取 m= 5,子= 17,tp = 26(平均周期),计算结果如图 6,由图中观察可知,在 0—100 的
范围内具有直线上升趋势,采用最小二乘法对该范围进行直线拟合,获得直线段斜率,即估计最大 Lyapunov
指数为 姿1 = 0.0125,估计最大预测时间为 80d。 非零数的指数估计和关联维数表明日径流量同样具有混沌
特性。
3摇 结果讨论
采用两种方法即 C鄄C方法和互信息法进行了时间序列时间延迟的判定,从叶绿素 a观测序列和日径流量
时间序列的分析结果可知,两种方法各自获得的时间延迟比较接近,如表 1,图 4所示。 在重构相空间时,子 不
3955摇 17期 摇 摇 摇 赵晓东摇 等:基于混沌理论的河流藻类生长特性分析———以德国易北河为例 摇
http: / / www.ecologica.cn
摇 图 6摇 1997—2001 年易北河日径流量时间序列 Lyapunov 指数
(Rosenstein算法)
Fig. 6 摇 The largest Lyapunov index of daily runoff sequence
(1997—2001)
取最佳值只会影响关联维数的计算,并不会影响重构吸
引子如实反映系统的动力学性质[27]。 子 取值太小,迫
使空间轨迹沿同一方向挤压,信息被压缩,获取困难。
反之,导致相邻时刻的动力学状态变化剧烈,使得简单
的系统运动变得看起来很复杂。 采用 C鄄C 方法获得的
子依赖于时间延迟窗口 子w,子w反映了数据依赖的最大时
间,而互信息方法仅仅反映了第一次局部最大时间,不
能估计 子w。 另外,评价重构吸引子的质量是其几何上
的重复性和不相关性,重构吸引子应具有较低的重复性
和较强的相关性。 因此,子w是一种用于估计维数的更好
的量。 尽管 C鄄C 方法是利用统计结果获得,但通过和
互信息对比表明,在确定序列时间延迟上仍具有较好的
效果,而且该方法能同时确定嵌入窗口和嵌入维数,为
较好的估计 Lyapunov指数提供了前提条件。
为更好的验证 C鄄C 方法计算嵌入维数的可靠性,
采用了 G鄄P 算法对嵌入维数进行估计。 计算结果表明,通过 G鄄P 算法获得的饱和嵌入维数与 C鄄C 方法获得
嵌入维数相近,如图 2和图 5所示。 这进一步验证了 C鄄C方法的可靠性。
根据 G鄄P 算法获得关联维数过程中,对无标度区的估计采用了视觉识别法。 无标度区的判断尽管有很
多种方法[25鄄26, 32],但在无法获知理论关联维数的情况下,视觉识别法仍可作为一种方便快捷的估计方法。 如
图 2所示,当 m取较小值时,其直线段(无标度区)较为明显,直线区间范围较大。 随着 m的增加,其直线段范
围发生变化,并有逐渐变小的趋势,如图 2中的曲线包围范围所示,这与文献[23]描述结果一致。 该方法获得
关联维数的精度对判定系统混沌特性以及确定关联维数并无显著影响[33],而且由 C鄄C 方法确定的嵌入维数
均比获得由关联维数确定的最小嵌入维数大得多(>2D+1),确保展开的相空间与原系统的动力学特性拓扑
等价。 同样,叶绿素 a序列的低关联维数估计值表明藻类生长在较低维空间产生自相似,但藻类生长的复杂
性可能使得需要采取高的嵌入维数值才能完全展现藻类的生长轨迹特征。
时间序列的数据长度应能保证其奇异吸引子可以被围绕多次,这样就能够对该序列的吸引子进行较好的
统计采样。 Smith等[34]认为 Nmin>42m,显然条件过于苛刻。 洪时中等[35]对此进行了重新理论推导,认为序列
最小长度应满足 Nmin> 2·( 27.5 ) D。 林振山等[36]则认为数据不必过多,满足 Nmin >260 即可。 实际上,本
文采用的单年叶绿素 a的时间序列长度 4824,嵌入维数范围为 30—45,关联维数为 2.75—4.02,符合洪时中等
人提出的最小数据长度要求。 而对于单年的日径流量而言,最少的数据序列只有 245,无法满足林振山等人
的最小数据量要求。 因此,将连续 5a 的数据作为连接起来作为一个时间序列,总长为 1704。 该序列获得的
嵌入维数为 9,关联维数为 1.84,符合洪时中等[35]提出的样本要求。 由上述分析可知,采用的时间序列符合混
沌序列的分析要求。
采用 Rosenstein算法对最大 Lyapunov指数进行了估计,先后对叶绿素 a和日径流量时间序列进行了指数
估计。 结果表明,两种时间序列的最大 Lyapunov指数估计值均较小,如表 2,图 6 所示,这表明该序列的邻近
吸引子轨道离散程度较低,可以进行较长时期的预测。 其中,5a 的日净流量的最大预测时间估计为 80 d,5a
的叶绿素 a最大预测时间平均为 14 d,前者要远大于后者。 可见,尽管河流径流量对叶绿素 a有影响,但可能
并不是对后者产生混沌特性影响的主要因素。
河流水体中的藻中生长的另一个重要因素是天气条件,而气候系统已被证实是具有混沌特性。 目前天气
预报可以被较好预测的时间约为 15 d,和叶绿素 a观测序列的最大预测时间相吻合。 相比于水文条件,气候
条件可能是对叶绿素 a的混沌特性产生影响的一个主要因素。
4955 摇 生摇 态摇 学摇 报摇 摇 摇 35卷摇
http: / / www.ecologica.cn
在 Elbe河流中的硅藻几乎总是优势藻种,其限制性营养盐为硅。 因此,叶绿素 a 所反映的主要是硅藻浓
度,始终作为优势藻,使得由于藻种间的竞争导致的混沌现象似乎不可能发生,另外,其限制性营养盐硅输入
主要是来自于 Elbe河源头的 Elbe砂岩山的硅砂岩土层,而沿途补充的营养源多为氮磷,因此,由于营养盐输
入变化导致的藻种浓度出现混沌特征可以忽略。 这也说明了气象因素要相对于其他因素更有可能造成藻种
叶绿素 a浓度的混沌变化。
4摇 总结
采用了混沌理论对易北河叶绿素 a 1997 年至 2001 年间各年的观测序列进行了混沌分析,分析结果表
明,各年的叶绿素 a观测序列均具有低维混沌特性,同样在对 5a 易北河连续日径流量时间序列分析表明,该
径流量也具有低维混沌特性。 通过对两种序列的最大可预测时间比对表明径流量等水文因子可能并不是影
响叶绿素 a混沌特性的主要因素,相反气象因子可能占据主导地位。
致谢:Mr. Scharfe博士提供德国易北河历史数据,特此致谢。
参考文献(References):
[ 1 ]摇 Huisman J, Weissing F J. Biodiversity of plankton by species oscillations and chaos. Nature, 1999, 402(6760): 407鄄410.
[ 2 ] 摇 Connell J H. Diversity in tropical rain forests and coral reefs. Science, 1978, 199(4335): 1302鄄1310.
[ 3 ] 摇 Sala S, Vighi M. GIS鄄based procedure for site鄄specific risk assessment of pesticides for aquatic ecosystems. Ecotoxicology and Environmental Safety,
2008, 69(1): 1鄄12.
[ 4 ] 摇 Malmaeus J M, Blenckner T, Markensten H, Persson I. Lake phosphorus dynamics and climate warming: a mechanistic model approach. Ecological
Modelling, 2006, 190(1 / 2): 1鄄14.
[ 5 ] 摇 Hu W P, J覬rgensen S E, Zhang F B. A vertical鄄compressed three鄄dimensional ecological model in Lake Taihu, China. Ecological Modelling, 2006,
190(3 / 4): 367鄄398.
[ 6 ] 摇 饶群, 芮孝芳. 完全混合系统总磷随机模型研究. 水科学进展, 2002, 13(1): 21鄄25.
[ 7 ] 摇 Xu F L, J覬rgensen S E, Tao S, Li B G. Modeling the effects of ecological engineering on ecosystem health of a shallow eutrophic Chinese lake
(Lake Chao) . Ecological Modelling, 1999, 117(2 / 3): 239鄄260.
[ 8 ] 摇 Hassan H, Hanaki K, Matsuo T. A modeling approach to simulate impact of climate change in lake water quality: phytoplankton growth rate
assessment. Water Science and Technology, 1998, 37(2): 177鄄185.
[ 9 ] 摇 赵晓东, 潘江, 李金页, 陶晓磊, 庞坤. 铜绿微囊藻和斜生栅藻非稳态营养盐限制条件下的生长竞争特性. 生态学报, 2011, 31(13):
3710鄄3719.
[10] 摇 Scharfe M, Callies U, Bl觟cker G, Petersen W, Schroeder F. A simple Lagrangian model to simulate temporal variability of algae in the Elbe River.
Ecological Modelling, 2009, 220(18): 2173鄄2186.
[11] 摇 Tsonis A A, Elsner J B. The weather attractor over very short timescales. Nature, 1988, 333(6173): 545鄄547.
[12] 摇 Radhakrishnan P, Dinesh S. An alternative approach to characterize time series data: case study on Malaysian rainfall data. Chaos, Solitons and
Fractals, 2006, 27(2): 511鄄518.
[13] 摇 丁晶, 王文圣, 赵永龙. 长江日流量混沌变化特性研究———域相空间嵌入维数的确定. 水科学进展, 2003, 14(4): 412鄄416.
[14] 摇 丁晶, 王文圣, 赵永龙. 长江日流量混沌变化特性研究———玉相空间嵌入滞时的确定. 水科学进展, 2003, 14(4): 407鄄411.
[15] 摇 Sivakumar B. Rainfall dynamics at different temporal scales: a chaotic perspective. Hydrology and Earth System Sciences, 2001, 5(4): 645鄄652.
[16] 摇 Takens F. Detecting strange attractors in turbulence. Lecture Notes in Mathematics, 1981, 898: 366鄄381.
[17] 摇 Packard N H, Crutchfield J P, Farmer J D, Shaw R S. Geometry from a time series. Physical Review Letters, 1980, 45(9): 712鄄716.
[18] 摇 Kugiumtzis D. State space reconstruction parameters in the analysis of chaotic time series鄄the role of the time window length. Physica D: Nonlinear
Phenomena, 1996, 95(1): 13鄄28.
[19] 摇 Kim H S, Eykholt R, Salas J D. Nonlinear dynamics, delay times, and embedding windows. Physica D: Nonlinear Phenomena, 1999, 127(1 / 2):
48鄄60.
[20] 摇 Fraser A M, Swinney H L. Independent coordinates for strange attractors from mutual information. Physical Review A, 1986, 33(2): 1134鄄1140.
[21] 摇 Grassberger P, Procaccia I. Characterization of strange attractors. Physical Review Letters, 1983, 50(5): 346鄄349.
5955摇 17期 摇 摇 摇 赵晓东摇 等:基于混沌理论的河流藻类生长特性分析———以德国易北河为例 摇
http: / / www.ecologica.cn
[22]摇 王安良, 杨春信. 评价奇怪吸引子分形特征的 Grassberger鄄Procaccia算法. 物理学报, 2002, 51(12): 2719鄄2729.
[23] 摇 党建武, 黄建国. 基于 G. P 算法的关联维计算中参数取值的研究. 计算机应用研究, 2004, 21(1): 48鄄51.
[24] 摇 韩敏. 混沌时间序列预测理论与方法. 北京: 中国水利水电出版社, 2007.
[25] 摇 党建武, 施怡, 黄建国. 分形研究中无标度区的计算机识别. 计算机工程与应用, 2003, 39(12): 25鄄27.
[26] 摇 巫兆聪. 分形分析中的无标度区确定问题. 测绘学报, 2002, 31(3): 240鄄244.
[27] 摇 吕金虎, 陆君安, 陈士华. 混沌时间序列分析及其应用. 武汉: 武汉大学出版社, 2002.
[28] 摇 Rosenstein M T, Collins J J, De Luca C J. A practical method for calculating largest Lyapunov exponents from small data sets. Physica D: Nonlinear
Phenomena, 1993, 65(1 / 2): 117鄄134.
[29] 摇 Rathje E M, Abrahamson N A, Bray J D. Simplified frequency content estimates of earthquake ground motions. Journal of Geotechnical and
Geoenvironmental Engineering, 1998, 124(2): 150鄄159.
[30] 摇 Rathje E M, Faraj F, Russell S, Bray J D. Empirical relationships for frequency content parameters of earthquake ground motions. Earthquake
Spectra, 2004, 20(1): 119鄄144.
[31] 摇 刘式达, 梁福明, 刘式适, 辛国君. 自然科学中的混沌和分形. 北京: 北京大学出版社, 2003.
[32] 摇 陈敏. 一种简单的 GP 算法无标度区识别方法. 计算机与信息技术, 2007, (11): 35鄄36, 39鄄39.
[33] 摇 姬翠翠, 朱华, 江炜. 混沌时间序列关联维数计算中无标度区间识别的新方法. 科学通报, 2010, 55(31): 3069鄄3076.
[34] 摇 Smith L A. Intrinsic limits on dimension calculations. Physics Letters A, 1988, 133(6): 283鄄288.
[35] 摇 洪时明, 洪时中. 用 Grassberger鄄Procaccia方法计算吸引子维数的基本限制. 物理学报, 1994, 43(8): 1228鄄1233.
[36] 摇 林振山. 长期预报的相空间理论和模式. 北京: 气象出版社, 1993.
6955 摇 生摇 态摇 学摇 报摇 摇 摇 35卷摇