基于時間序列的微藻生長環境參數的預測模型

聶磊 蔡文濤 黃一凡 董正瓊

摘要：? 掌握微藻生長環境中的溫度、硝酸鹽濃度、氧氣濃度等參數的變化規律有利于提高生物產出量。本研究利用鈍頂螺旋藻全生長周期中的大量環境參數，引入時間序列預測方法，通過極大似然估計，得到線性回歸模型，同時引入狀態空間預測方法，將低維的線性模型映射到高維空間，以消除低維空間下模型的不準確性。為評估模型的有效性，采用了杜賓-瓦特森檢驗法、圖像法、均方根誤差、平均絕對誤差和最大誤差等評估方法和指標進行評估和驗證。預測結果顯示，自回歸滑動平均模型和自回歸滑動平均-卡爾曼濾波模型均可用于預測微藻生長的環境參數;與前者相比，后者預測誤差更小，模型擬合更精確，能更好地揭示鈍頂螺旋藻生長過程中環境參數的內在規律。

關鍵詞：? 鈍頂螺旋藻; 自回歸滑動平均模型; 自回歸滑動平均-卡爾曼濾波模型; 環境參數

中圖分類號：? 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建模。差分處理，可以去除序列的趨勢性，對序列進行一次差分處理后，得到的新序列通?？梢詽M足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] ?？柭鼮V波的基本公式如下：

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.

（責任編輯：陳海霞）