聂磊 蔡文涛 黄一凡 董正琼

摘要：  掌握微藻生长环境中的温度、硝酸盐浓度、氧气浓度等参数的变化规律有利于提高生物产出量。本研究利用钝顶螺旋藻全生长周期中的大量环境参数，引入时间序列预测方法，通过极大似然估计，得到线性回归模型，同时引入状态空间预测方法，将低维的线性模型映射到高维空间，以消除低维空间下模型的不准确性。为评估模型的有效性，采用了杜宾-瓦特森检验法、图像法、均方根误差、平均绝对误差和最大误差等评估方法和指标进行评估和验证。预测结果显示，自回归滑动平均模型和自回归滑动平均-卡尔曼滤波模型均可用于预测微藻生长的环境参数;与前者相比，后者预测误差更小，模型拟合更精确，能更好地揭示钝顶螺旋藻生长过程中环境参数的内在规律。

关键词：  钝顶螺旋藻; 自回归滑动平均模型; 自回归滑动平均-卡尔曼滤波模型; 环境参数

中图分类号：  Q949.2    文献标识码： A    文章编号：  1000-4440（2021）05-1183-07

Prediction model of microalgae growth environment parameters based on time series

NIE Lei， CAI Wen-tao， HUANG Yi-fan， DONG Zheng-qiong

（College of Mechanical Engineering， Hubei University of Technology/Hubei Key Laboratory of Manufacture Quality Engineering， WuHan 430068， China）

Abstract：  Mastering the variation law of parameters such as temperature， nitrate concentration and oxygen concentration in the growth environment of microalgae is beneficial to improve the biological output. In this study， a large number of environmental parameters in the full growth cycle of  Spirulina platensis  were used， and the time series prediction method was introduced to obtain a linear regression model by maximum likelihood estimation， while a state space prediction method was introduced to map the low-dimensional linear model to the high-dimensional space to eliminate the inaccuracy of the model in the low-dimensional space. Evaluation methods and indicators such as the Durbin-Watson test， image method， root mean square error， mean absolute error and maximum error were used to assess the validity of the model. The prediction results showed that both the autoregressive moving average model and the autoregressive moving average-Kalman filter model could be used to predict the environmental parameters of microalgae growth. Compared with the autoregressive moving average， the autoregressive moving average-Kalman filter model has smaller prediction error and more accurate model fitting， which can better reveal the intrinsic patterns of environmental parameters during the growth of  Spirulina platensis .

Key words：   Spirulina platensis ; autoregressive moving average model; autoregressive moving average-Kalman filter model; environmental parameters

在過去的二十年里，微藻细胞作为蛋白质、维生素、必需氨基酸和脂肪酸的重要来源而受到关注  [1-5] 。螺旋藻是一种丝状蓝藻，它由蓝绿色的圆柱形细胞的丝状螺旋毛状体组成  [6] ，可以在许多其他生物不适宜的环境中定居  [7] ，是微藻培养的典型藻种。螺旋藻被开发为“未来的食物”，因为它具有比其他藻类更高效地合成高质量浓缩食物的惊人能力。最值得注意的是，螺旋藻蛋白质含量超过50%，所有必需氨基酸都处于完美平衡状态。它的光合转化率为8%～10%，相比之下，大豆等陆生植物只有3%。

户外跑道池系统是当下商业藻类生产最常用的方法  [8-10] ，它试图最大限度地将太阳能和无机化学物质转化成高价值产品。目前，利用户外跑道池系统培养藻类的运营成本较高，其主要原因在于最佳种群密度、最适温度难以维持以及过饱和光照度易引起高溶解氧水平降低等问题。针对上述问题，已有许多研究结果表明微藻的生长速度受到周围环境参数的影响。例如， Baky等  [11] 比较了螺旋藻在不同氮浓度下的生长，结果表明，在氮限制条件下，螺旋藻的总类胡萝卜素和总脂肪含量最高。Li等  [12] 采用CO  2 吸收和微藻（CAMC）转化混合系统，研究不同含氮量对微藻培养工艺性能的影响，发现当 NH  4 HCO  3  与NaNO  3 的比例为 1∶ 4时，复合工艺的碳利用率可达40.45%，高于传统的微藻CO  2 固定工艺（约 10%～ 30%）。李姿等  [13] 在恒温条件下，比较了4个不同生长期内的光照度和CO  2 体积分数的最适条件，得到的钝顶螺旋藻最大生物量和固碳速率分别为4.126  g/L 和58.231  mg/（L·h） 。赵嶝科等  [14] 研究了Ni  2+ 对螺旋藻生长的影响，试验结果表明，质量浓度0.5  μg/ml 的Ni  2+ 可促進螺旋藻生长。掌握跑道池系统参数变化规律，预测未来参数数值，是实现微藻培养过程优化的重要前提。

本研究以钝顶螺旋藻为研究对象，采用时间序列模型对其培养过程中的温度、硝酸盐浓度和氧气浓度等变量的变化进行预测，随后对时间序列模型进行卡尔曼滤波，以消除单一模型预测的局限性，揭示微藻整个生长周期中环境变量的变化规律，为其延迟期、对数期前期、对数期后期和稳定期中环境变量的调控提供依据，为发现低成本、高产出的微藻培养条件参数组合打下基础。

1 材料与方法

1.1 试验材料及数据来源

1.1.1 试验材料  试验材料为从美国模式培养物集存库（American type culture collection，ATCC）中获得的蓝藻S.platensis 29408;户外跑道池中接种1.5 L钝顶螺旋藻和200.0 L培养基。

1.1.2 试验仪器  四叶片浆轮，以6.875  r/min 速度旋转;Hobo数据记录器;15 ml Falcon管;Dionex离子色谱仪DX-120;克拉克式氧微电极与Unisense皮安计。

1.1.3 试验环境  试验地点在美国加利福尼亚州马斯堡，室外跑道水池系统是在屋顶温室中构建的，在整个试验过程中，空气温度是唯一的控制因子。

1.1.4 试验数据的测定  温度数据的测定：温度数据（以摄氏度为单位）使用漂浮在每个滚道中的Hobo数据记录器收集，每5 min测量1次，数据每周从数据记录器下载。硝酸盐浓度的测定：硝酸盐浓度数据收集每周3次，从每个滚道每天收集的25 ml样本中取10 ml以9  g 的速度在15 ml Falcon离心管中离心10 min，倒出上清液并储存在30 ml闪烁小瓶中，储存在 -20 ℃ 的冷冻室中，备用。使用型号为DX-120的Dionex离子色谱仪分析500 μl来自闪烁瓶的硝酸盐样品。氧气浓度的测定：从每个滚道收集2份50 ml样品来测量氧气浓度，每周3次。首先，将克拉克式微电极耦合到Unisense皮安计上，放入样品试管中，记录下电极的电流，并采用两点法对电极进行校准  [15] ，然后分别向2份样品中通入氮气和空气约3 min，同时记录电流，最后对照水饱和空气溶氧表获得此时的氧气含量。

1.2 参数预测方法

1.2.1 时间序列模型的建立  时间序列分析是一种处理随时间变化的相关数据的数学方法  [16-18] 。该方法的基本思想是：在不同的时间，同一变量的观测值之间存在着一定的联系，通过分析观测值的序列，找出它所反映的事物发展过程和趋势，然后进行类推或延伸，以达到预测的目的  [19-22] 。时间序列分析中最常用的自回归滑动平均模型（Autoregressive moving average，ARMA）将平稳的有色噪声序列视为与时刻有关的序列，其一般形式为：

x k=φ 1x  k-1 +φ 2x  k-2 +…+φ px  k-p +ε k-θ 1ε  k-1 - θ 1ε  k-1 -… -θ qε  k-q   （1）

其中， φ i <1（ i =1，2，…， p ），为回归系数; θ i <1 （ i = 1，2，…， q ），为滑动平均系数; p 、 q 为ARMA（ p ，  q ）模型的阶次; ε k 为误差。当 q 为0时，该模型是自回归模型AR（ p ），当 p 为0时，该模型是滑动平均模型MA（ q ）。建立时序模型的步骤如图1所示：

（1） 平稳性检验

在对时间序列进行ARMA建模之前，要对序列进行平稳性检验。本研究采用ADF（Augmented dickey-fuller）检验，也称单根检验，是一种时序模型建立常用的检验方法。

（2） 差分

在实际测量中，由于测量目标不断变化，测量得到的数据是不平稳的，不能对其直接进行ARMA建模。差分处理，可以去除序列的趋势性，对序列进行一次差分处理后，得到的新序列通常可以满足ARMA建模平稳性的需求。

（3） 零均值化

对于录入的试验数据，选取若干个数据组成新的数据序列，并将其零均值化，处理数据速度快，能较好地符合建立ARMA模型所需的条件。

（4） 序列的自相关函数（ACF）和偏自相关函数（PACF）

由ACF和PACF的图像可以观察数据是否平稳，同时根据图像是否拖尾，来分别判断模型中的自回归部分 p 和滑动平均部分 q 的最大值。二者计算公式如下：

ρ ^   k=n=k t=1 （x t-x - ）（x  t+k -x - ）n t=1 （x t-x - ） 2   （2）

φ ^    kk =  D ^   k  D ^     （3）

ρ ^   k 為第 k 个自相关函数，  x  t  为序列中第 t 个数据， x  - 为序列均值， x 为序列中第 t+k 个数据， D  ^ 为自相关函数与偏自相关函数互相转换的系数矩阵，  D  ^   k是根据克拉默法则得到的矩阵。公式（3）中，

D ^ =

1   ρ ^   1  …  ρ ^    k-1

ρ ^   1  1  …  ρ ^    k-2

ρ ^    k-1   ρ ^    k-2  …  1   ;

D ^  k=

1   ρ ^   1  …  ρ ^    1

ρ ^   1  1  …  ρ ^    2

ρ ^    k-1   ρ ^    k-2  …   ρ ^    k    。

（5） 模型的定阶

赤池信息准则 （Akaike information criterion，AIC）和贝叶斯信息准则 （Bayesian information criterion，BIC）是ARMA建模常用的定阶方法，本研究选取BIC准则作为拟合模型的定阶准则。二者的计算方法如下：

BIC（p，q） =lg （σ 2  （p，q） ）+ p+q n  lg n  （4）

AIC（p，q） =lg （σ 2  （p，q） ）+ 2 n （p+q）  （5）

公式（4）和（5）中， p 、 q 为ARMA模型的阶数， σ 2  （p，q）  为模型残差的方差， n 为样本个数。

（6） 模型检验

模型定阶后指定模型结构，采用最小二乘法计算模型参数，对得到的模型需检验其可用性，即检验模型的残差是否为白噪声，是否具有相关性，若模型残差为白噪声且不具有相关性，则模型可用;否则，模型不可用。

1.2.2 卡尔曼滤波模型的建立  考虑到单一ARMA模型预测结果存在着滞后、不精确等问题，本研究采用卡尔曼滤波方法对模型预测的结果进行修正，以时序模型建模得到的线性方程为卡尔曼滤波模型的状态方程，以状态方程为依据建立卡尔曼滤波模型的预测部分和修正部分，从而实现对时序模型的优化。

卡尔曼滤波器（Kalman filter，KF）是一种以线性随机系统的状态方程和量测方程为基本方程，用状态方程进行递推，通过计算线性最小方差估计，滤除观测数据及系统中噪声和干扰的影响，从而得到最优自回归数据的处理算法  [23-24] 。卡尔曼滤波的基本公式如下：

X（k）=AX（k-1）+BU（k）+W（k）  （6）

Z（k）=HX（k）+V（k）  （7）

公式（6）和（7）中， X（k ）是 k 时刻系统的状态，  U（k ） 是 k 时刻系统的控制量。 A 、 B 是系统的系数，对多变量系统模型，它们是矩阵。 Z（k ）是 k 时刻系统的测量值， H 是测量系统系数，对于多变量系统，它是矩阵。 W（k ）和 V（k ）分别是模型和测量过程的噪声。

进一步根据公式（6）和（7）推导出如下的卡尔曼滤波方程：

X（k|k-1）=AX（k-1|k-1）+BU（k）  （8）

P（k|k-1）=AP（k-1|k-1）A T+Q  （9）

X（k|k）=X（k|k-1）+Kg（k）[Z（k）-HX（k|k-1）]  （10）

Kg（k）=P（k|k-1）H T/[HP（k|k-1）H T+R]  （11）

P（k|k）=[1-Kg（k）H]P（k|k-1）  （12）

公式（8）中， X（k|k -1）是利用上一状态预测的结果， X（k -1| k -1）是上一状态的最优结果， U（k ）为控制量， P（k|k -1）是 X（k|k -1）对应的协方差，  P（k -1| k -1） 是 X（k -1| k -1）对应的协方差， Q 是系统过程的噪声方差， R 是测量过程的噪声方差，  Kg（k） 为 k 时刻的卡尔曼增益， I 为全1矩阵。获取ARMA模型后，系统的状态方程和测量方程可表示为：

X 1（k+1） X 2（k+1）    X p（k+1） =  φ 1 φ 2 … φ p  1  0  0  0  0  1  0  0  0  0  …  1  X 1（k） X 2（k）   X p（k） + 1 θ 1 θ 2 θ q ω（k+1）  （13）

Z（k+1）= 1 0    0  X 1（k+1） X 2（k+1）    X p（k+1） +v（k+1）  （14）

1.2.3 模型的评估与验证  首先，采用残差图分析和QQ图（分位数图示法）分析对ARMA模型进行正态性检验，然后利用相关系数、偏自相关系数图像法和杜宾-瓦特森（Durbin-watson，DW）法检验其相关性，最后采用均方根误差（ RMSE ）、平均绝对误差（ MAE ）和误差极差对整体的建模精度开展评估。

DW=   n k （e k-e  k-1 ） 2   n k e 2 k   （15）

误差极差=max （x k-x ^  k）  （16）

MAE= 1 nn k=1 |（x k-x ^  k）|  （17）

RMSE=  1 nn k=1 （x k-x ^  k） 2   （18）

其中： e  k  是 k 时刻的残差， e  k    -1 是 k -1时刻的残差， x  k  是 k 时刻序列实际值， x ^  k 是 k 时刻序列预测值， n 是样本数量， DW 为杜宾-瓦特森统计量的值。

模型的相关性可通过计算 DW 值进行评价，数值越接近2，说明模型不具有相关性，若其接近0或4，则说明模型分别有负相关和正相关的性质，模型不可用。误差极差评价模型误差的最大波动情况， MAE 评价模型整体的波动， RMSE 衡量预测准确度可以说明模型预测值的离散程度。

2 结果与分析

2.1 自回归滑动平均-卡尔曼滤波模型的预测结果

選取培养微藻时的温度、硝酸盐浓度和氧气浓度3种环境参数的数据作为ARMA模型的输入，得到不同参数条件下ARMA模型表达式，以该表达式作为卡尔曼滤波模型的基本公式，最终得到自回归滑动平均-卡尔曼滤波模型（ARMA-KF）。拟合得到3种环境参数的ARMA模型阶数及参数（表1），2种模型的预测结果如图2、图3和图4所示。其中，ARMA-KF模型所对应的测量数据与拟合数据的匹配程度最高，故可利用该模型预测微藻生长过程中环境参数的变化规律。

2.2 模型的评估与验证

2.2.1 模型的验证  考虑到ARMA模型对混合模型的精度影响很大，对ARMA模型的相关性和正态性分别进行了检验，结果如图5、图6和图7所示。由图5、图6、图7中可以看出，标准化残差十分平稳，且QQ图中大多数残差落在图中虚线上，可以判断模型基本接近正态分布，由残差的自相关函数图与偏自相关函数可初步得出模型无相关性。通过DW方法检验模型自相关性得到检验值分别为1.97、2.06、2.02，说明模型无相关性。此时，原始数据中有用的信息提取完毕，可对序列进行预测。

2.2.2 模型的评估  利用训练数据建模的过程中，ARMA模型得到的拟合曲线存在着一定的误差，且明显滞后于实际曲线。将ARMA模型与卡尔曼滤波结合后，曲线的3个误差指标下降了约60%（表2）。说明ARMA-KF模型很适合运用于钝顶螺旋藻生长环境指标的预测。

3 结 论

环境因素对微藻的生长有着重要的影响，也是限制微藻大量生产的主要因素。本研究建立的ARMA模型能够对微藻生长的环境因素如营养物质浓度、氧气浓度等进行预测，但其存在局限性，即预测结果存在一定的滞后，而ARMA-KF模型能很好地消除ARMA模型的局限性，且其预测结果优于单一的ARMA模型，大大地提高了模型的可用性。本研究模型可运用于微藻生产的各项重要环境指标参数的预测，为其生长过程的调控提供理论依据，同时也为其他微生物的生长过程预测提供了新思路。本模型对于变化较为缓慢的线性环境参数较为适用，对于环境参数变化较为剧烈的非线性参数适用性不强，下一步可采用自适应的预测方法，将线性和非线性预测模型组合起来，进一步掌握环境参数的变化规律。

参考文献：

[1]  ZHAO X L， KATRIN S， GERO C， et al. Hydrothermal carbonization of Spirulina platensis and Chlorella vulgaris combined with protein isolation and struvite production[J]. Bioresource Technology Reports， 2019， 1（6）：159-167.

[2] IMAR D S， REGIN M Q， ALHADJ D M， et al. Vitamin a status in healthy women eating traditionally prepared spirulina （ Dihé ） in the Chad Lake area[J]. PLoS One， 2018， 13（1）： e0191887.

[3] 闫春宇，胡冰涛，王素英. 常压室温等离子体诱变对螺旋藻中氨基酸成分的影响[J]. 食品与发酵工业， 2017， 43（1）： 60-65.

[4] CHEN G S， CAI Y， SU Y Y， et al. Effects of Spirulin a algae as a feed supplement on nutritional value and flavour components of silkie hens eggs[J]. Journal of Animal Physiology and Animal Nutrition， 2019， 103（5）：1408-1417.

[5] 罗光宏，雷 婕，张喜峰，等. 碳酸二甲酯三相分配体系萃取分离螺旋藻多糖[J]. 食品与发酵工业， 2020， 46（17）： 196-203.

[6] 朱孝晨. 钝顶螺旋藻藻蓝蛋白与多糖的制备及生物活性研究[D]. 烟台：烟台大学， 2018.

[7] 王雪飞. 光生物反应器中市政污水培养钝顶螺旋藻的条件优化[D]. 无锡：江南大学， 2017.

[8] CARVALHO A P， MEIRELES L A， MALCATA F X. Microalgal reactors： a review of enclosed system designs and performances[J]. Biotechnology Progress， 2006， 22（6）：1490-1506.

[9] NORSKER N H， BARBOSA M J， VERMUE M H， et al. Microalgal production--a close look at the economics[J]. Biotechnology Advances， 2011， 29（1）：24-27.

[10] 余 灿. 城市污水户外跑道池培养微藻的生长促进方法研究[D]. 哈尔滨：哈尔滨工业大学，2015.

[11] BAKY H E， BAROTY G E， MOSTAFA E M .  Optimization growth of spirulina （ Arthrospira ） platensis in photobioreactor under varied nitrogen concentration for maximized biomass， carotenoids and lipid contents.[J]. Recent Patents on Food， 2020， 11（1）：40-48.

[12] LI S H， SONG C F， LI M D， et al. Effect of different nitrogen ratio on the performance of CO  2  absorption and microalgae conversion （CAMC） hybrid system[J]. Bioresource Technology， 2020， 306： 123126.

[13] 李 姿，任洪艳，周骞骞，等. 光强和CO  2 体积分数对不同生长时期的钝顶螺旋藻生长和固碳的交互影响[J]. 食品与生物技术学报， 2014， 33（8）： 827-836.

[14] 赵嶝科，刘 慧，张少斌，等. Ni  2+ 对钝顶螺旋藻生长及藻胆蛋白含量的影响[J]. 水产科学， 2020， 39（5）：1-5.

[15] 李日清.在线pH仪在污水处理系统中的使用及维护[J].山西电子技术，2011，4（5）：35-36.

[16] 李 桐，董维红，张琦琛，等.基于时间序列模型的黑龙江省粮食水足迹分析与预测[J].排灌机械工程学报，2020，38（11）：1152-1159.

[17] 甄晓菊，张雪红，吴国明，等.基于Sentinel-2A NDVI时间序列数据的冬小麦识别[J].江苏农业科学，2019，47（16）：239-245.

[18] 江显群，陈武奋，邵金龙，等. 大坝变形预报模型应用[J]. 排灌机械工程学报，2019， 37（10）：870-874，920.

[19] 吳 婕. 数据驱动的雷达健康预测技术研究[D].北京：中国电子科技集团公司电子科学研究院， 2019.

[20] 李 波，林 聪，刘清蝉，等. 基于时序建模的光纤电流互感器随机噪声卡尔曼滤波方法[J]. 电机与控制学报， 2017， 21（4）： 83-88，.

[21] 张瑞国，李春雨，丁志宏，等. 基于时间序列模型的雷达数据随机误差建模与补偿[J]. 现代雷达， 2016， 38（11）： 49-52.

[22] 唐 俊，赵成萍，周新志. 基于EVI-RBF的玉米长势监测及产量预测[J].江苏农业学报，2020，36（3）：577-583.

[23] 孙振华. 基于动态数据驱动的生物氧化槽进气量预测研究[D]. 乌鲁木齐：新疆大学， 2018.

[24] 修春波，任 晓，李艳晴，等. 基于卡尔曼滤波的风速序列短期预测方法[J]. 电工技术学报， 2014， 29（2）： 253-259.

（责任编辑：陈海霞）