技术领域
[0001] 本
发明涉及一种蓝藻水华预测方法，属于水环境预测技术领域，具体地说，是指一种基于非线性动力学时序模型的蓝藻水华暴发预测预警方法。
背景技术
[0002] 随着经济社会的发展，
水体富营养化已经成为一个全球性的重大水环境问题。水体富营养化是指在人类活动的影响下，
生物所需的氮、磷等营养物质大量进入湖泊、河湖、海湾等缓流水体，引起藻类及其他浮游生物迅速繁殖、水体溶解
氧量下降、鱼类及其他生物大量死亡的现象。藻类水华是水体富营养化的典型表现，是我国河湖水环境污染的主要问题之一。藻类暴发性繁殖，致使水质恶化、缺氧、产生腥臭等异味物质，甚至产生藻毒素并通过食物链对人、畜和水生生物造成毒害，继而破坏河流生态的
稳定性，严重影响城市供水和
饮用水安全，危害人类健康，制约区域经济的可持续发展，给人类生活带来极其恶劣的影响，水华问题已被列入全球性的水污染问题。因此，深入研究藻类水华生成过程，对藻类水华暴发这一非常规突发事件进行有效模拟和预测，对促进水环境保护和技术进步具有重要意义。
[0003] 目前蓝藻水华预测方法研究主要包括机理驱动模型和数据驱动模型。机理驱动模型包含单一营养盐模型、浮游
植物生态模型、生态动力学模型等，这些模型中含有许多生态变量和待定参数，通过描述水华发生的机理过程进行建模，在认识水华和治理水华
发面发挥重要作用；而数据驱动模型包括数理统计模型和
人工智能模型，这些模型从大量的实测数据中通过智能
算法等技术挖掘隐藏于系统的相关信息，能很好的模拟非线性过程，尤其适用于机理不明确的高维非线性系统，在水华预测中得到了运用。上述两类建模方法具有各自的优势，但是也存在一些问题，机理驱动模型多是针对影响水华形成的生物、化学和物理因素等方面展开研究，但是水华的发生仍然是一个机理不太清楚的非线性问题，还需要进行深入的研究，同时，模型中大量的待定参数也限制了模型的适用性；而数据驱动模型缺乏机理
支撑，无法合理解释各种影响因素与水华生成的因果关系，难以保证模型的准确性。
[0004] 在对水华预测的研究中，现有的机理驱动模型大多将蓝藻生长动力学系统视为时不变系统，其模型参数被设为不随时间变化的恒定值，因而忽略了实际水体中水温、光照等影响因素随时间的变化对蓝藻生长率和死亡率等模型参数的影响，无法有效描述实际蓝藻生长时变系统的动力学特征和动态过程，从而难以实现蓝藻水华的准确预测。现有的数据驱动模型采用了多种智能方法，比如
人工神经网络、
支持向量机、模糊理论、
决策树等，对水华影响因子的主观取舍会给准确预测水华带来不稳定的因素，并且会存在局部最小问题、泛化能力较差、预测
精度有限等问题。
[0005] 机理驱动建模和数据驱动建模具有各自的优势和局限性，因此，如何将两者结合，建立一个预测精度高的水华
预测模型，是蓝藻水华研究领域中亟待解决的问题，对促进水环境保护和技术进步具有重要意义。
发明内容
[0006] 本发明目的是解决现有的蓝藻水华预测精度不高、忽略部分影响因素对蓝藻生长时变系统的影响等问题，针对实际蓝藻生长时变系统中采集到的叶绿素a浓度、氮磷营养盐浓度、水温和光照等实测数据，考虑光照、水温等影响因素随时间变化对蓝藻生长的影响，将蓝藻生长率作为时变参数，建立带有双营养盐循环的蓝藻生长非线性动力学时序模型，采用数值算法和智能
进化算法相结合对蓝藻生长非线性动力学时序模型中的定常参数进行优化率定，通过对光照、水温等影响因素进行多元时间序列建模，实现对时变参数及蓝藻生物量的预测，并采用分岔理论和中心流形理论对蓝藻生长时变系统进行非线性动力学分析，得到蓝藻水华暴发的条件，进而实现对水华暴发行为的预警。本发明提出的基于蓝藻生长非线性动力学时序模型的蓝藻水华预测方法，不仅确定了蓝藻水华暴发的条件，也提高了水华预测精度，为环保部
门提供有效的参考依据，为水环境治理提供了治理决策。
[0007] 本发明提供的基于蓝藻生长非线性动力学时序模型的蓝藻水华预测方法，包括以下五个步骤：
[0008] 步骤一、蓝藻生长双营养盐循环非线性动力学建模；
[0009] 蓝藻生长与氮、磷营养盐负荷相互作用和关系是揭示水华形成机理的主要途径，可通过建立氮磷营养盐与蓝藻生物量之间的关系来表示水华形成过程。叶绿素a浓度是表征水体中蓝藻生物量的最直接指标，可采用叶绿素a浓度作为反映蓝藻水华生成的表征因素，叶绿素a浓度的变化率与蓝藻生长率、死亡率密切相关，其中蓝藻生长率主要由水温、光照等水华影响因素决定，根据Monod方程(莫诺方程)得到的氮磷双营养盐与蓝藻生长率的关系，建立叶绿素a浓度与氮磷双营养盐之间的非线性动力学模型，模拟蓝藻生长过程中叶绿素a浓度与氮磷双营养盐之间的变化。
[0010] 步骤二、蓝藻生长非线性动力学时序建模；
[0011] 考虑到蓝藻生长率在实际环境中并不是不变的定常参数，而是随时间变化的时变参数，将步骤一中建立的非线性动力学模型中蓝藻生长率作为受光照、水温等因素的影响随时间变化的时变参数，建立蓝藻生长非线性动力学时序模型。
[0012] 步骤三、蓝藻生长非线性动力学时序模型参数优化率定；
[0013] 蓝藻生长非线性动力学时序模型包含定常参数和时变参数。对于定常参数，结合实际水体中采集到的数据进行优化率定。由于建立的非线性动力学时序模型含有三个非线性微分方程，属于三变量非线性微分方程组，且定常参数的数量较多，采用数值算法和智能进化算法中的
遗传算法相结合，进行三变量非线性微分方程组的定常参数优化率定。对于时变参数，由于蓝藻生长率与光照、水温等影响因素有关，且水温和光照等影响因素的数值按时间顺序变化，同时考虑多种影响因素间的相关性，将其监测数据视为多元时间序列。从而实现对蓝藻生长非线性动力学时序模型预测。
[0014] 步骤四、蓝藻水华暴发非线性动力学条件确定；
[0015] 针对已经确定定常参数和时变参数的蓝藻生长非线性动力学时序模型，采用分岔理论研究蓝藻生长时变系统的平衡点及其稳定性。由于蓝藻生长非线性动力学时序模型为三变量非线性微分方程组，其非线性动力学分析复杂，需要在分岔理论研究
基础上，采用中心流形理论进行
降维处理，以进一步研究蓝藻生长时变系统
临界点的分岔行为，确定各个范围内蓝藻生长时变系统的稳定性，得到系统稳定时蓝藻生长率的范围，确定蓝藻水华暴发的非线性动力学条件。
[0016] 步骤五、蓝藻水华暴发预警；
[0017] 根据步骤四中的水华暴发非线性动力学条件，得到蓝藻生长时变系统稳定时蓝藻生长率的变化范围，再根据预测的水温、光照等影响参数随时间变化对蓝藻生长率的影响，判断蓝藻生长率是否超过蓝藻的稳定生长率范围，即是否存在蓝藻水华暴发的
风险，从而对蓝藻水华暴发行为进行风险预警。并且结合实际水环境中水华发生的灾害点，与叶绿素a浓度的预测值对比，即可预测蓝藻水华暴发的时刻,从而对水华暴发灾害预警。
[0018] 本发明的优点在于：
[0019] 1、本发明建立了蓝藻生长双营养盐循环非线性动力学模型，在此基础上，考虑了实际水体中水温、光照等影响因素随时间的变化对蓝藻生长率的影响，构建蓝藻生长非线性动力学时序模型，对蓝藻生长时变系统的描述更符合实际。
[0020] 2、本发明在对蓝藻生长非线性动力学时序模型定常参数优化率定时，采用数值方法与智能进化算法相结合的方法求取最优参数，使模型具有适用性和实践应用性。克服现有研究不适用于对三变量非线性微分方程组参数率定的问题。
[0021] 3、本发明在蓝藻水华预测时，结合机理驱动和数据驱动预测的优势，将两种预测方法相结合，先对影响因子时序建模预测，预测影响因素随时间的变化，并且反映到蓝藻生长非线性动力学时序模型的时变参数中，进而对蓝藻的叶绿素a浓度预测。本发明既能反映复杂水环境下各种影响因素对藻类水华形成的影响，又能发挥数据驱动搜索建模的优势，可适用于不同湖库水体环境要求，具有可扩展性。
[0022] 4、本发明提出动力学临界点和实际水环境灾害点相结合的方法对蓝藻水华的暴发进行预警的方法。根据分岔理论和中心流形理论研究蓝藻生长时变系统稳定性转变的临界点，在此临界点由于系统的振荡作用，当叶绿素a浓度的最大值大于或等于实际水环境的灾害点，表示会发生水华暴发。本发明的方法可及时预测水华暴发这一突发性水质灾害，为决策部门制定水华的防治对策提供依据。
附图说明
[0023] 图1是本发明基于非线性动力学时序模型的蓝藻水华预测方法的
流程图；
[0024] 图2是遗传算法结合数值算法的参数率定流程；
[0025] 图3是原始建模数据与参数率定后的模型拟合结果对比图；
[0026] 图4是蓝藻生长率在G0＝29.92附近的分岔图；
[0027] 图5是第二个平衡点随生长率(G>G0)变化图；
[0028] 图6～图9是G＝38时的系统相轨迹、叶绿素a浓度、氮浓度、磷浓度的时间历程图；
[0029] 图10～图13是G＝25时的系统相轨迹、叶绿素a浓度、氮浓度、磷浓度的时间历程图；
[0030] 图14～图17是G＝50时的系统相轨迹、叶绿素a浓度、氮浓度、磷浓度的时间历程图；
[0031] 图18是水温时间序列预测图；
[0032] 图19是光照时间序列预测图；
[0033] 图20是叶绿素、氮、磷浓度时变相轨迹图；
[0034] 图21是时变叶绿素预测值；
[0035] 图22是时变氮营养盐预测值；
[0036] 图23是时变磷营养盐预测值；
[0037] 图24是时变动力学预测值、源数据、时间序列预测值对比图；
[0038] 图25是蓝藻生长率随时间变化图。
具体实施方式
[0039] 下面将结合附图对本发明作进一步的详细说明。
[0040] 本发明提供一种基于非线性动力学时序模型的蓝藻水华预测方法，如图1所示流程，具体步骤如下：
[0041] 步骤一、蓝藻生长双营养盐循环非线性动力学建模；
[0042] 蓝藻生物量可通过叶绿素a的含量指标来衡量，通过建立氮磷营养盐与叶绿素a浓度之间的关系来表示蓝藻生长过程。叶绿素a的浓度变化率与蓝藻的生长率、净损失率密切相关，经过对蓝藻生长率的影响因素分析，发现蓝藻生长率主要由光照、水温等影响因素决定，因此将蓝藻生长率由光照、水温的作用函数表示。根据营养盐限制条件下藻类生长率的关系，建立营养盐增长函数时使其符合Monod方程，建立的蓝藻生长双营养盐循环非线性动力学模型如下所示：
[0043]
[0044] 其中，ca为t时刻叶绿素a的浓度，N、P分别为t时刻营养盐总氮、总磷的浓度；N0、P0分别为营养盐总氮、总磷在初始时刻的浓度；G为蓝藻生长率，T1为藻类的环境损失率；T2为氮营养盐的耗损率，gN为蓝藻对氮营养盐的吸收率，dN为藻类代谢分解对氮营养盐的贡献率；T3为磷营养盐的耗损率，gP为蓝藻对磷营养盐的吸收率，dP为藻类代谢分解对磷营养盐的贡献率。g(N)、g(P)分别为营养盐总氮、总磷增长函数，表达式为：
[0045]
[0046] 其中，KN、KP分别代表蓝藻生长对氮和磷的半饱和浓度常数。在实际水体中，蓝藻生长率G由水温T、光照I等影响因素决定，以水温T和光照I为例，其关系如下：
[0047]
[0048] 其中，GT(T),GI(I)分别为水温和光照对蓝藻生长率的影响函数，gmax为蓝藻最大生长率，kI为光照半饱和浓度。
[0049] 步骤二、蓝藻生长非线性动力学时序建模；
[0050] 在实际水体中，由于水温和光照等影响因素随时间变化而变化，由这些影响因素决定的蓝藻生长率也同样随时间而变化，并非如公式(1)中所定义的恒定不变。因此，蓝藻生长双营养盐循环动力学系统是时变系统，本发明在公式(1)基础上，引入时变的蓝藻生长率G(t)，构建蓝藻生长双营养盐循环时变参数非线性动力学模型：
[0051]
[0052] 其中，时变的蓝藻生长率G(t)随时间变化的公式为：
[0053]
[0054] 其中，T(t),I(t)分别为水温和光照随时间变化的t时刻的数值。
[0055] 步骤三、蓝藻生长非线性动力学时序模型参数优化率定；
[0056] 对于定常参数，根据实际水体监测数据，结合遗传算法与数值算法，对蓝藻生长双营养盐循环非线性动力学模型即公式(1)中的定常参数(T1、T2、gN、dN、T3、gP、dP、KN、KP、gmax、kI)进行优化率定。定常参数优化率定流程如图2所示，具体步骤为：
[0057] (1)设定初始化条件及初始化种群个体。确定个体数目、最大遗传代数、待优化参数个数、代沟、适应度
阈值。采用多参数级联浮点数编码。随机产生若干种不同的参数组合作为初始种群，构成待优化参数空间。所述的待优化参数包括T1、T2、gN、dN、T3、gP、dP、KN、KP、gmax、kI。
[0058] (2)适应度值评价。种群在进化过程中，通过适应度值来评价个体的优劣，建立适应度函数：
[0059]
[0060] 其中，F为适应度值，n为采集水环境数据的总天数；cat为叶绿素a在t时刻的叶绿素a浓度真实值，ca(t)为叶绿素a在t时刻的叶绿素a浓度函数值；Nt为营养盐在t时刻的总氮浓度真实值，N(t)为营养盐在t时刻的总氮浓度函数值；Pt为营养盐在t时刻的总磷浓度真实值，P(t)为营养盐在t时刻的总磷浓度函数值。
[0061] 由于多元微分方程组公式(1)的结构复杂，无法得到ca(t)和N(t)的解析表达式，需采用数值算法计算，本发明采用数值算法中的4阶龙格库塔算法，则其数值积分表达分别为：
[0062]
[0063] 其中，参数qi(i＝1,2,3,4)具体表达式如下：
[0064] q1＝f1(ca(t),N(t),P(t))
[0065] q2＝f1(ca(t)+q1×h/2,N(t)+q1×h/2,P(t)+q1×h/2)
[0066] q3＝f1(ca(t)+q2×h/2,N(t)+q2×h/2,P(t)+q2×h/2)
[0067] q4＝f1(ca(t)+q3×h,N(t)+q3×h,P(t)+q3×h)
[0068] f1(x,y,z)＝T1x+Gg(y)g(z)x
[0069] 其中，h为步长，T1是藻类的环境损失率，g(y)、g(z)分别与g(N)、g(P)的函数形式相同。x,y,z为函数f1(x,y,z)的自变量。
[0070] 同理，参数ri(i＝1,2,3,4)具体表达式如下：
[0071] r1＝f2(ca(t),N(t),P(t))
[0072] r2＝f2(ca(t)+r1×h/2,N(t)+r1×h/2,P(t)+r1×h/2)
[0073] r3＝f2(ca(t)+r2×h/2,N(t)+r2×h/2,P(t)+r2×h/2)
[0074] r4＝f2(ca(t)+r3×h,N(t)+r3×h,P(t)+r3×h)
[0075] f2(x,y,z)＝T2(N0-y)-gNg(y)g(z)ca+dNx
[0076] 参数si(i＝1,2,3,4)具体表达式如下：
[0077] s1＝f3(ca(t),N(t),P(t))
[0078] s2＝f3(ca(t)+s1×h/2,N(t)+s1×h/2,P(t)+s1×h/2)
[0079] s3＝f3(ca(t)+s2×h/2,N(t)+s2×h/2,P(t)+s2×h/2)
[0080] s4＝f3(ca(t)+s3×h,N(t)+s3×h,P(t)+s3×h)
[0081] f3(x,y,z)＝T3(P0-z)-gPg(y)g(z)ca+dPx
[0082] 根据数值方法计算获得各个时刻叶绿素a、氮、磷浓度值，将得到的叶绿素a、氮、磷浓度值与实际监测的叶绿素a、氮、磷浓度真实值代入适应度函数公式(6)中，获得种群各个个体的适应度值F。根据公式(6)可知，当适应度值越小，表明数值方法的计算值越接近实际监测值，说明该个体优良，会参与下一代种群的繁殖。
[0083] (3)选择、交叉、变异、重插入运算。
[0084] (4)遗传算法终止。根据个体的适应度值不断重复第(2)、(3)步形成新个体，终止条件是重复次数达到最大的遗传代数或参数组合的最小适应度值小于或等于适应度阈值。取达到终止条件时参数组合适应度值的最小值作为蓝藻生长非线性动力学时序模型的适应度值，该参数组合为蓝藻生长非线性动力学时序模型的最优参数组合。
[0085] 对于时变参数，由于水温T(t)和光照I(t)等影响因素的函数值按时间顺序变化，同时考虑多种影响因素间的相关性，将其监测数据视为多元时间序列。根据多元时间序列模型预测水温T(t)和光照I(t)等影响因素随时间变化的数值，将其预测结果代入公式(5)得到时变参数蓝藻生长率G(t)的预测值。将蓝藻生长率G(t)预测值代入到蓝藻生长非线性动力学时序模型中，根据龙格库塔算法对叶绿素a浓度、氮、磷营养盐浓度预测的计算公式(7)，即可实现对蓝藻生长非线性动力学时序模型的预测。
[0086] 步骤四、蓝藻水华暴发非线性动力学条件确定；
[0087] 蓝藻生长时变系统是高维动力学系统，分岔行为较复杂，本发明对蓝藻水华暴发非线性动力学条件确定的方法是将蓝藻生长率视为时变参数，首先采用分岔理论确定蓝藻生长非线性动力学时序模型的分岔集和分岔点，再采用中心流形方法对蓝藻生长非线性动力学时序模型降维得到约化方程，具体研究临界局部稳定性和分岔的类型，确定蓝藻水华暴发的临界点，并且提出了临界点满足的公式条件。
[0088] (1)首先采用分岔理论确定蓝藻生长非线性动力学时序模型的分岔集和分岔点。具体方法为：令方程组(4)中三个方程的右边都等于零，解方程组得到蓝藻生长非线性动力* * * * * *
学时序模型的平衡点Q(ca ,N ,P)，ca、N、P分别为平衡点处的叶绿素a浓度、总氮浓度和总磷浓度，为了便于分析进行坐标平移，令，
[0089]
[0090] 坐标平移后的方程组变换为：
[0091]
[0092] 分离线性项和高次项得到：
[0093]
[0094] 其中 k为线性部分的雅克比矩阵， 为高次项。
[0095] 通过求解方程(8)，得到蓝藻生长非线性动力学时序模型方程线性部分的雅克比矩阵：
[0096]
[0097] 高次项：
[0098]
[0099]
[0100]
[0101] 其中，h.o.t表示展开式中三次及其以上的高次项。
[0102] 由方程(1)可知，蓝藻生长时变系统的平衡点Q(ca*,N*,P*)有两种情况：
[0103] (I)(0,N0,P0),N0,P0分别为初始时刻两种营养盐的浓度；
[0104] (II)(ca*,N*,P*)，其中ca*,N*,P*均不为0。
[0105] 下面根据在不同情况下求解雅克比矩阵的特征值考虑蓝藻生长时变系统的稳定性状况。
[0106] (I)对于第一个平衡点(0,N0,P0)，其物理意义对应于初始状态，蓝藻生物量的浓度为0，此时的雅克比矩阵为：
[0107]
[0108] 由分岔理论，蓝藻生长时变系统的稳定性取决于对应的雅克比矩阵(10)的特征值。第一个平衡点的三个特征值分别为 -T2，-T3。由于T2、T3分别是氮、磷营养盐的耗损率，均为正数，当第一个特征值的 特征值为
正数，状态是稳定的；当特征值 为负数，状态是不稳定的，所以
当 时，特征值为0，蓝藻生长时变系统处于
临界状态，令
则G0为分岔点。
[0109] (II)对于第二个平衡点(ca*,N*,P*)，是三个变量都存在的情况。其中ca*,N*,P*,G(t)，满足的关系式为：
[0110]
[0111] 将关系式(11)中的后两个式子代入第一个式子可得到只含有叶绿素a浓度和蓝藻生长率的一元二次方程，求解所述的一元二次方程可以得到关于叶绿素a浓度的两个根，对应的平衡态为(ca1*,N1*,P1*)和(ca2*,N2*,P2*)，对应的雅克比矩阵为：
[0112]
[0113] 由于在参数表示下，无法显式表达雅克比矩阵(12)的三个特征值。但是可根据雅克比矩阵的特征值方程求出蓝藻生长时变系统的分岔点，即当特征值出现零根时对应的蓝藻生长率，设为G1。经过计算得到特征值方程为：λ3+aλ2+bλ+c＝0
[0114] 其中
[0115]
[0116] 当特征根出现0根时，满足c＝0，与方程式(11)联立，即可求出ca*,N*,P*,G(t)，令此时的G(t)值表示为G1，则G1也是分岔点。
[0117] 通过上述分析可知蓝藻生长时变系统存在两个分岔点，所以可以将分岔集分成三个范围(0，G0)、(G0，G1)和大于G1的范围，并且当蓝藻生长率取不同量值时，蓝藻生长时变系统的第一个平衡点无变化，第二个平衡点会随蓝藻生长率变化。
[0118] (2)确定蓝藻生长时变系统的分岔点后，为确定系统稳定时蓝藻生长率的范围以及蓝藻暴发的临界点，需要研究各个范围内蓝藻生长时变系统的稳定性和临界点系统的分岔情况，具体分析如下：
[0119] 第一个范围(0，G0)：根据蓝藻生长双营养盐非线性动力学模型对第一个平衡点的分析可知，当蓝藻生长率处于(0，G0)范围时，蓝藻生长时变系统处于稳定状态，会稳定在第一个平衡点。当生长率 时，特征值为0，蓝藻生长时变系统处于临界状态，需要研究蓝藻生长时变系统在分岔点G0附近的局部稳定性。蓝藻生长时变系统处于临界状态时，非线性系统的局部稳定性不能采用传统的线性化稳定性原理来分析，而是决定于系统的非线性项，且非线性系统的临界局部稳定性经常是与系统动态行为的分岔行为紧密相关。而对于高维动力学系统，分岔行为较复杂，需要采用降维措施化为低维方程组再研究，本发明采用中心流形方法对蓝藻生长时变系统降维得到约化系统，确定蓝藻生长率为G0附近系统的分岔行为，具体步骤为：
[0120] 为了得到约化系统，令 通过变换即可知Gc在0附近的分岔行为。此时蓝藻生长时变系统的雅克比矩阵为：
[0121]
[0122] 三个特征值是0，-T2，-T3，对应的
特征向量组成的线性变换矩阵为：
[0123]
[0124] 通过线性变换(u,v,w)T＝M(x,y,z)T，则原来的微分方程(8)可以转化为：
[0125] φ'＝Dφ+M-1f(Mφ)。
[0126] 其中,D＝diag(0-T2-T3),φ＝(x,y,z)T,并且将蓝藻生长率视为分岔变量时，展开得到：
[0127]
[0128] 根据中心流形理论，蓝藻生长时变系统的局部动力学特性主要由中心流形上的流决定，该系统有二维中心流形，设，
[0129]
[0130] 其中，O((x+Gc)3)是函数展开式的高次项。
[0131] 为了确定方程组(15)中的系数a1,a2,a3,b1,b2,b3，将(15)代入微分方程组(14)中关于求解y、z的后两个方程，比较同次幂系数后得到系数值，再代入(14)中关于求解x的第一个方程，得到中心流形上的方程即约化方程，并且将生长率参数转化为G(t)，可得到约化方程如下：
[0132]
[0133] 其中，
[0134]
[0135] 为了分析的方便，对上面所求的约化方程用x'＝f(x,G(t))＝lx+mx2+nx3表示，分别用字母l,m,n表示系数，由于a1＜0,b1＜0，判断约化方程(16)的系数的正负性可得m＞0、n＜0。对约化方程分析可知：方程f(x,G(t))＝0有两条解曲线，对于解曲线x＝0，导数值Dxf(x,G(t))＝a，故当a＜0即生长率 时，平衡点渐进稳定，当G(t)＞G0时，不稳定；对于另一条解曲线是在对称轴为 的抛物线，当解存在时，需要满足：
[0136] 并且在时，两个解都大于0，但是Dxf(x,G(t))＞0表明平衡状态不稳定，当蓝藻生长率G(t)＞G0时，两个解一正一负，对应的平衡态设为(ca1*,N1*,P1*)和(ca2*,N2*,P2*)，并令ca1*＞0,ca2*＜0，对于ca2*无实际研究意义，因此只考虑ca1*，并且导数值Dxf(x,G(t))＜0，说明G(t)＞G0时平衡状态(ca1*,N1*,P1*)稳定。
[0137] 通过中心流形对蓝藻生长率为G0附近系统的分岔行为分析，说明G0点是超临界叉形分岔点，不是蓝藻水华的暴发点。分析得到当蓝藻生长率在(0，G0)范围内时，蓝藻生长时变系统会稳定在第一个平衡点，当蓝藻生长率在G0和大于G0的一段范围时，蓝藻生长时变系统会稳定在第二个平衡点。
[0138] 继续对第二个平衡点分析，对于有意义的平衡点(ca1*,N1*,P1*)，在分岔集求解时可知第二个平衡点对应的雅克比矩阵的特征值不能显式表达，但是可以求出特征值出现0根时对应的蓝藻生长率G1，即G1满足：
[0139]
[0140] 当生长率处于(G0，G1)范围，特征值均小于0，平衡点会稳定在第二个平衡点，当大于G1后，会出现特征值大于0的情况，蓝藻生长时变系统不稳定。说明G1点是蓝藻水华暴发的临界点。临界点G1的取值满足方程组(17)。
[0141] 通过上述分析，当蓝藻生长率在(0，G1)范围内，蓝藻生长时变系统是稳定的，即在(0，G0)时，蓝藻生长时变系统会稳定在第一个平衡点，当蓝藻生长率处于(G0，G1)时，蓝藻生长时变系统会稳定在第二个平衡点；当生长率大于G1时，蓝藻生长时变系统不稳定，发生蓝藻水华暴发行为。
[0142] 步骤五、蓝藻水华暴发预警；
[0143] 当蓝藻生长率作为时变参数随时间变化时，蓝藻生长时变系统的第一个平衡点无变化，第二个平衡点会随时间变化。由于平衡点不再是一个固定值，而是随时间变化，系统的稳定性也随时间变化。根据步骤三中对系统分岔集和分岔点的分析可知：当蓝藻的生长率大于G1时，蓝藻生长时变系统不稳定，有发生水华暴发行为的风险，其中蓝藻的生长率是由光照、水温等影响因素决定，故蓝藻水华暴发风险的条件公式是：
[0144] G(t)＝GT(T(t))·GI(I(t))＞G1 (18)
[0145] 其中，G1由方程组(17)求解得出。通过时间序列模型预测到水温和光照值，代入蓝藻生长率函数，当预测出蓝藻生长率大于蓝藻生长率的暴发点G1时，有发生水华暴发的风险。
[0146] 再根据预测到的蓝藻生长率，结合蓝藻生长非线性动力学时序模型，可预测蓝藻生物量以及营养盐浓度值，当预测的蓝藻生物量大于实际水环境的水华灾害点阈值(设为caT)时，即满足条件：
[0147] ca(t+1)＝ca(t)+h×(q1+2q2+2q3+q4)/6＞caT (19)
[0148] 表示发生水华暴发行为，根据预测到的蓝藻水华发生时刻，提前为蓝藻水华暴发行为预警。
[0150] 以江苏省太湖流域数据为例采用本发明所提方法进行蓝藻水华预测和蓝藻水华预警。数据主要包括叶绿素a浓度、总氮浓度、总磷浓度、水温、光照等数据。为方便分析，对叶绿素a浓度、总氮、总磷等的原始监测数据进行了标准化及异常点剔除等预处理，如图3中无标记的曲线所示。
[0151] 根据2010年至2011年太湖共400天的1个表征因子(叶绿素a浓度)和四个影响因子(总氮浓度、总磷浓度、水温、光照)的实测数据，实测数据类型及单位具体见表1。
[0152] 表1实测数据类型及单位
[0153]
[0154] 步骤一、蓝藻生长双营养盐循环非线性动力学建模；
[0155] 见具体实施方式中的公式(1)-(3)。
[0156] 步骤二、蓝藻生长非线性动力学时序建模；
[0157] 见具体实施方式中的公式(4)-(5)。
[0158] 步骤三、蓝藻生长非线性动力学时序模型参数优化率定；
[0159] 对于定常参数，采用本发明提出的数值算法与智能
优化算法相结合的参数优化率定方法对公式(1)中的定常参数进行率定，率定结果如表2。
[0160] 表2蓝藻生长非线性动力学时序模型定常参数率定结果
[0161]参数 约束范围 率定结果
蓝藻最大生长率gmax 0.01～100 24.7
光照半饱和浓度kI 0.01～10 0.08
藻类的环境损失率T1 0.01～30 21.6
氮营养盐的耗损率T2 0.01～20 10.3
磷营养盐的耗损率T3 0.01～1 0.09
蓝藻对氮的吸收率gN 0.01～110 2.68
藻类代谢分解对氮的贡献率dN 0.01～100 57
蓝藻对磷的吸收率gP 0.01～1 0.014
藻类代谢分解对磷的贡献率dP 0.01～1 0.23
蓝藻生长对氮的半饱和常数KN 0.01～10 3
蓝藻生长对磷的半饱和常数KP 0.01～5 0.51
氮浓度初始值N0 0.01～50 23
磷浓度初始值P0 0.01～10 3
[0162] 将定常参数率定后的蓝藻生长非线性动力学时序模型对叶绿素a浓度的拟合结果如图3中带“×”形标注曲线所示，其中营养盐浓度和水温、光照均采用的实际监测数据，从图3中实际数据与拟合数据对比可见定常参数率定的结果较为理想。
[0163] 对于时变参数，为了对蓝藻浓度和氮磷营养盐浓度预测，需要首先对光照、水温等影响因素预测，所以先对光照、水温建立二元时间序列模型，并且将二元时间序列模型按趋势项、周期项、随机项进行特征项分解，对光照和水温进行预测，水温和光照的时序数据预测结果如图4和图5的带“×”形标注曲线所示。
[0164] 将水温和光照的时序数据预测值代入时变参数蓝藻生长率公式中，则蓝藻生长时变系统相轨迹如图6，其叶绿素a浓度和氮磷营养盐浓度的时间历程如图7、图8、图9所示。
[0165] 将蓝藻生长时变系统的叶绿素a浓度时间历程(预测值)与其监测数据、直接采用时间序列预测方法进行对比，曲线如图10中所示，表明蓝藻生长时变系统比直接采用时间序列方法的预测结果精度更高。
[0166] 步骤四、蓝藻水华暴发非线性动力学条件确定；根据理论分析，首先计算第一个平衡点下的叉形分岔点G0，带入率定参数值可计算G0＝29.92。再根据中心流形理论和分岔理论研究G0点的分岔情形。经线性变换令Gc＝G(t)-G0＝G(t)-29.92，即转换至研究Gc在零附近时蓝藻生长时变系统的分岔。此时，三个特征值为0、-10.3、-0.09，三个特征值对应的特征向量构成的特征矩阵为：
[0167]
[0168] 将时变参数蓝藻生长率视为变量时，原来的微分方程则转化为：
[0169]
[0170] 其中f1,f2,f3均为关于f(-0.619x,0.965x+y,x+z)的函数。
[0171] 根据中心流形定理，设，
[0172]
[0173] 为了确定两个函数的系数值，分别将两个函数带入代入微分方程组的后两个式子，
[0174] 得到：
[0175] 可知蓝藻生长时变系统有二维中心流形，计算后得到中心流形上的方程即约化方程为：
[0176]
[0177] 将生长率参数转化为G，可得到：
[0178]
[0179] 约化方程的状态变量随参数变化的分岔图如图11所示。根据蓝藻生长率G0在29.92附近的分岔图可知：当蓝藻的生长率小29.92时，对于第一个平衡点是稳定的，大于
29.92时，第一个平衡点不稳定，当蓝藻生长时变系统存在微小的
波动时，就会稳定在第二个平衡点处，故蓝藻生长率在29.92处属于叉形分岔。
[0180] 当蓝藻的生长率大于29.92时，继续对第二个平衡点分析，求解蓝藻生长时变系统的分岔点G1，通过对方程组(9)的求解，得到G1＝48.85。随着在生长率的变化，第二个平衡点的变化图如图12所示，在蓝藻生长率为48.05时，特征根出现0根，即是在第二个平衡点下蓝藻生长时变系统由稳定状态转化为不稳定状态的临界点。当蓝藻的生长率大于48.85时，蓝藻生长时变系统不稳定，会发生水华暴发现象。因此对蓝藻生长时变系统的整体分析可知，当蓝藻的生长率处于(0，29.92)时，系统会稳定在第一个平衡点，当蓝藻的生长率处于(29.92，48.05)区间中，系统会稳定在第二个平衡点，大于48.05时，系统不稳定，会发生水华暴发行为。
[0181] 下面通过选取不同的生长率，说明不同生长率下蓝藻生长时变系统的稳定性。当生长率取为38时，先分析系统的第一个的平衡点(0,23,3)，雅克比矩阵对应的三个特征值分别为7.13、-10.3、-0.09，特征值含有正数，所以第一种平衡点是不稳定的；对第二个的平衡点分析，当生长率为38时，分析叶绿素a的浓度值为正值时的情况，通过计算得到平衡点(2.86,22.85,0.92)，对应的雅克比矩阵的三个特征值为-0.2093+1.1882i、-0.2093-1.1882i、-11.1336，表明蓝藻生长时变系统是振荡稳定的。G＝38时的蓝藻生长时变系统的相轨迹如图13，相应的叶绿素a浓度和氮、磷营养盐浓度的时间历程如图14、图15、图16所示。
[0182] 由图13～图16可见，当营养盐积累到一定程度，同时水温、光照等影响因素满足一定条件后，蓝藻(叶绿素a浓度)生长速度较快，即蓝藻生长时变系统产生振荡。随着蓝藻生长对营养盐的大量消耗，当营养盐消耗到不足以维持蓝藻的生长时叶绿素a浓度也开始降低，由于死亡藻类的降解和藻类的代谢，营养盐的浓度又会开始上升，如此反复，最终会稳定在某一水平上。所以公式(1)能较好的解释蓝藻的生长情况，在一定程度上符合实际情况。
[0183] 当蓝藻的生长率为25时，第一个平衡点(0,23,3)对应的矩阵的三个特征值分别为-2.69、-10.3、-0.09，三个特征值均为负值，表示蓝藻生长时变系统会稳定在第一个平衡点，不会发生蓝藻水华暴发行为。蓝藻生长时变系统的相轨迹如图17，相应的叶绿素a浓度和氮、磷营养盐浓度的时间历程如图18、图19、图20所示，叶绿素a浓度曲线会维持在零附近，说明系统是稳定的。
[0184] 当蓝藻的生长率为50时，第一个平衡点(0,23,3)对应的矩阵的三个特征值分别为16.2，第二个平衡点中有意义的平衡点值是(45,34,0.3)，对应的矩阵的三个特征值分别为-7.5530、0.6598、20.1382，特征根含有正数，表示蓝藻生长时变系统是不稳定的，会发生蓝藻水华暴发行为。蓝藻生长时变系统的相轨迹如图21，相应的叶绿素a浓度和氮、磷营养盐浓度的时间历程如图22、图23、图24所示，叶绿素a浓度曲线在突变至较大数值，表明系统不稳定，直接发生水华暴发现象。
[0185] 通过上面选取的不同生长率下的蓝藻生物量和营养盐的时间历程图，展现了蓝藻生长时变系统的不同的动力学行为。
[0186] 步骤五、蓝藻水华暴发预警；
[0187] 由于蓝藻的生长率是随光照、水温等影响因素随时间变化的，导致蓝藻生长时变系统的第二个平衡点也随时间变化，所以蓝藻生长时变系统的稳定性也随时间会发生改变。根据步骤三中的对分岔集和分岔点的讨论，可知当生长率大于48.05时，会存在发生蓝藻水华暴发的风险，用影响因素表示是：
[0188] G(t)＝GT(T(t))·GI(I(t))＞48.05
[0189] 通过时间序列预测的光照、水温值代入生长率函数，得到的蓝藻生长率随时间变化图如图25所示，图中圆圈标记处对应的蓝藻生长率为48.05，说明在第313天会发生蓝藻水华暴发现象。
[0190] 根据预测的蓝藻生长率，代入蓝藻生长非线性动力学时序模型，根据步骤二中的数值方法可预测蓝藻生物量的变化，当预测的蓝藻生物量大于或等于实际水环境水华灾害点(阈值为40)时，表明会发生水华暴发现象。根据图7中蓝藻叶绿素浓度随时间变化图，在第313天蓝藻叶绿素浓度预测值为43，大于40，确实发生蓝藻水华暴发行为。
[0191] 本发明建立的蓝藻生长非线性动力学时序模型不仅可以预测水环境中蓝藻生物量和营养盐的浓度，还可判别水华暴发的风险和发生，实现了对水华暴发的预测和预警，为提前做好水华防范工作提供了新方法。
