破损青贮牧草氧气与温度原位观测及变质预测模型研究

程　强　周海洋　单桂林　李梦华　BUESCHER W　孙宇瑞

(1.中国农业大学信息与电气工程学院， 北京 100083； 2.波恩大学农业工程研究所， 波恩 53115)

0　引言

目前，国内外使用较广的青饲料贮存方式包括：圆柱形裹包式、方形裹包式、青贮槽式、塔仓式、窖仓式等。在这几种方式中，无论从减少收获过程中青饲料损失考虑，还是从保证发酵质量与安全存储等方面综合衡量，裹包式青贮饲料都是最优的选择。裹包式青贮饲料在西方发达国家已经非常普遍，我国也在日渐普及过程中。然而，裹包式贮存也有一些缺陷，比如，因其自身表面积/体积比很大，在运输过程中塑封膜极易破损。氧气的侵入破坏了无氧贮存环境，各类好氧微生物大量繁殖，消耗饲料的营养成分，产生热量，引起温度异常上升，伴随着霉菌的滋生，导致青贮饲料腐败变质，产生有害的挥发物对环境造成影响[1]。在变质过程中，氧气与温度的变化是衡量各类好氧微生物活动的关键参数，也是反映与评价青贮饲料贮藏质量的重要指标[2-6]。

国际学者构建了氧气与温度动态模型用于预测青贮饲料有氧胁迫下的变质过程，并借助实验观测方法检验模型预测的准确性与适用性。美国学者PITT等[7]和MUCK等[8]建立了槽式青贮玉米饲料有氧变质过程的数学模型。随后，使用破坏性取样法抽取气体，再用气相色谱仪测定青贮玉米饲料内部氧气浓度变化，并用这些观测数据检验所构建的模型的适用性。RUXTON等[9-10]建立了裹包式青贮牧草饲料有氧变质过程的数学模型，WILLIAMS等[11]对该模型进行了验证。在后续研究中，其他学者运用类似方法测量有氧胁迫下青贮饲料内部氧气与温度变化[12]，但这些观测方法有一个共同的缺陷：抽取气体样本将会扰动饲料内部的气体分布，造成实验结果出现偏差。此外，抽气分析法每天只能取少数几个气体样本，无法实时记录氧气浓度变化，难以突破连续观测技术瓶颈，致使青贮饲料有氧变质风险预测的研究进展缓慢。随着科学技术的不断进步，SUN等[13]运用一种基于Clark氧电极原理的溶解氧传感器成功实现了青贮玉米饲料内部氧气与温度的原位实时观测，并对PITT等[7]提出的预测模型进行了评价[14]。而对于青贮牧草，其成分、颗粒大小、干质量、菌落类型与数量等理化生性质与青贮玉米差异明显，基于青贮玉米对象提出的原位观测方法与预测模型是否适用于青贮牧草还有待验证。程洪等[15]虽然测量了氧胁迫下裹包式青贮牧草饲料内部温升分布，但缺少氧气浓度的测量数据，也未涉及模型检验。为此，本文尝试运用Clark氧电极传感器原位观测裹包式青贮牧草内部氧气与温度变化，建立有氧变质过程中青贮牧草内部氧动态、温升以及微生物繁殖耦合模型。由于裹包式青贮饲料的对称性，将三维模型简化成二维模型，最后根据实验观测数据与模型预测数据的对比结果评价模型的适用性。

1　材料与方法

1.1　模型建立

青贮饲料有氧变质过程可包括氧气扩散、好氧微生物生长、营养物质分解、生化反应产热与传导。当裹包青贮的塑封膜破损或青贮饲料开窖后，氧气侵入，根据菲克扩散定律，氧气在青贮饲料中的扩散过程可运用二维抛物线型偏微分方程来描述[7,10]

(1)

式中ψ——氧气浓度，%

D——氧气在空气中的扩散率，0℃时D=0.064 1 m2/h

Φ——饲料孔隙度

ρ——饲料容重，kg/m3

τ——多孔介质中孔隙的弯曲度，与孔隙度(容重)有关，本文取常数0.667

Yg——实际增长率，取0.5

Rs——用于呼吸作用的氧气体积分数，取0.747 m3/kg

C——微生物浓度，kg/kg

μg——微生物增长率，h-1，受温度和氧气浓度影响

式(1)右边第2项为源/汇项，反映了微生物分解营养物质消耗的氧气，这一过程会产生热量，热量的二维传导过程可描述为

(2)

其中

cp=cp,w(1.451dm+1)

(3)

式中T——青贮饲料温度，K

cp——青贮饲料比热容，J/(kg·K)

k——青贮饲料热导率，W/(m·K)

ke——呼吸作用释放能量，取1.594×107J/kg

α——热量转化系数

dm——干物质含量，%

cp,w——水比热容，J/(kg·K)

式(2)右边第2项为源/汇项，反映了微生物在有氧环境下分解营养物质产生的热量，用化学方程式描述

C6H12O6+6O2→6CO2+6H2O+2 870 kJ/mol

(4)

模型参数的测定对于模型预测精度非常重要。式(1)中微生物生长模型可描述为

(5)

其中

μg=μgmaxfTfψ

(6)

(7)

式中μg-μm——微生物净生长率

μd——微生物死亡率[7,14]

μm——酵母菌实际维持率

fT——微生物生长率随温度变化的函数[16-17]

Km——常数，取0.001 45[7]

ψ0——空气中的氧气浓度，取21%

fψ——微生物生长率随氧气浓度变化的函数

由于酵母菌的繁殖是青贮饲料有氧变质初期的主导因素，因此本文主要模拟变质初期酵母菌繁殖过程。营养成分的消耗(L，%)计算公式为

(8)

饲料热导率k、孔隙度Φ和最大容重ρmax为

k=kw(1-Φ)+k0Φ

(9)

(10)

(11)

式中kw——水热导率k0——空气热导率模型求解过程详见文献[7,10,14]。

1.2　Clark氧电极基本原理

薄膜氧电极是为测量水中所溶解的氧气量而研制的极谱电极，其本质上是一个电化学电池。该电极由Clark于1953年最先研制的，也称为Clark氧电极，其结构如图1所示。Clark氧电极由嵌在不导电原料上的银电极和铂电极组成。其中，银极作为阳极或参比电极制成圆环形状，面积尽可能大一些，以降低电极表面电流密度，减少阳极的极化现象，使其电极电位不受外加电压的影响。铂极作为阴极或指示电极，通常是圆点状，置于银电极的中间，电化学反应就产生在铂电极上。用15～20 μm的聚四氟乙烯薄膜包裹住氧电极的表层，把电解质(饱和氯化钾溶液)充入氧电极与薄膜之间。因为水中溶解的氧气能通过乙烯膜而电解质溶液不能进入，从而消除被测溶液中存在的各类离子发生反应的影响。

根据亨利定律，电解质溶液中的溶解氧含量与气液界面处的氧分压呈比例关系[13]，即

p=kHCDO

(12)

式中p——待测气体中的氧分压，Pa

CDO——电解质溶液中的溶解氧含量，mol/L

kH——亨利常数，对于氧气而言，在298 K下取770 Pa·L/mol

由于亨利常数受温度影响，因此需要借助van’t Hoff方程进行校正

(13)

式中，C1值为1 700 K(对于氧气而言)。将式(13)代入式(12)中即可得到待测气体/液体中的氧分压。本文使用商业化的Clark氧电极(型号：DO-10，测量范围：0～20 mg/L，电极输出范围：0～2.5 V，精度：±3%满量程；氧电极内部集成热敏电阻：MF5，范围0～50℃，分辨率：0.1℃，北京联创思源测控技术有限公司)观测裹包式青贮饲料内部高温高湿环境下氧气与温度的变化。关于氧电极的标定方法以及二氧化碳对氧电极的输出影响详见文献[18]。

图1　Clark氧电极内部示意图Fig.1　Schematic diagram of Clark oxygen electrodes1.电流计　2.氯化钾电解液　3.橡胶圈　4.聚乙烯或聚四氟乙烯薄膜　5.铂极　6.银极　7.电源

1.3　实验方案设计

图2　实验设计与传感器安装位置Fig.2　Experiment design and sensor location

如图2所示，所测试的裹包式青贮饲料由德国CLAAS公司的Rollant455型打包机打包成圆柱形，底面直径约120 cm，高约120 cm，饲料种类为黑麦草。裹包式青贮饲料样本在打包时打包机参数设置相同，其粒径大小、体积与鲜容重非常接近，最大限度消除了其他物理因素对实验结果的影响。首先将氧电极和温度传感器(DS18B20型)原位植入未开封的青贮裹包饲料内部[14]，探头离饲料表面深度均为20 cm。温度传感器分别记录距离破损表面约20、60、100 cm处的数据，氧传感器记录距离破损表面约20 cm处的数据。传感器植入后立即封住植入区域，以确保青贮饲料内部的无氧环境。待青贮饲料内部因安装传感器而引入的微量氧气消耗完之后，将青贮饲料圆柱面一侧开封，使其暴露在空气中，模拟破损或饲喂过程空气侵入方式，数据采集器(自制多通道采集器，微处理器：MSP430F149，14位A/D，Texas Instruments, Dallas, Texas, USA)记录开封之后氧气、饲料不同位置温度以及环境温度的变化。

图3　开封后5 d内的观测数据Fig.3　Observation data of five days after opening plastic cover of bale silage

同时对2种不同干物质含量(DM-1和DM-2)的青贮牧草进行监测，每个水平测试3个裹包，观测数据取平均值，共6个裹包饲料样本。实验在室内环境下进行，共经历了5 d，室温保持在(27.3±1.1)℃，风速可以忽略。在实验开始阶段，每个裹包用取样器收集N个样本(N≥3)，每个样本质量大于30 g，用于分析微生物种类及含量初值，同时确定不同位置饲料鲜容重，具体分析方法详见文献[14]，分析结果如表1所示。由于裹包式青贮饲料长期水平放置，其自重改变了容重的纵向分布，上部饲料(DM-1：381 kg/m3；DM-2：402 kg/m3)较下部饲料容重变小(DM-1：597 kg/m3；DM-2：643 kg/m3)。模型仿真算法采用的是非稳态显式差分法[19]，最小间距设定为0.1 m，最小时间间隔为0.067 h，在Matlab(Mathworks, Inc., Version 7.0)开发平台进行仿真计算。

表1　青贮饲料参数分析结果Tab.1　Analysis of related parameters in silage

2　结果与讨论

2.1　氧气与温度原位观测动态变化与分析

裹包式青贮饲料开封后氧气侵入与温度变化过程如图3所示(图中显示的数据曲线是3个样本的平均值)。从整体上看，氧动态变化过程中峰值的出现表明氧气扩散与微生物消耗达到了瞬态平衡，在峰值之前氧气扩散占据主导，之后微生物耗氧占据优势，达到动态平衡后稳定在某个值[7,14]。图3a反映了高干物质含量饲料样本的观测数据，图3b是低干物质含量饲料样本的观测数据。相比之下不难看出：首先，氧气达到的峰值有些差异，高干物质含量饲料的峰值(15.1%)略大于低干物质含量饲料的峰值(14.2%)。其次在耗氧速率方面，高干物质含量的饲料的耗氧速率(约2.5 d)仅仅略低于低干物质含量(约2.3 d)，整体差异不显著。最后，在耗氧初始阶段，高干物质含量的饲料样本耗氧速率较慢，大约经过1.2 d后斜率变大，耗氧速率开始加快；而低干物质含量的饲料样本在开始阶段耗氧速率较快，经过0.5 d后斜率变小，速率降低。其可能原因在于低干物质含量的青贮饲料湿度较高，更易于酵母菌的存活与繁殖，并且水分活性及可溶性碳水化合物(WSC)含量也高，因而耗氧速率较高[20]。

高干物质含量饲料温度增幅整体上比低干物质含量饲料的温度增幅高。其原因在于低干物质含量的饲料水分含量较高，本身的比热容也较大，在释放同等数量级热量的条件下，低干物质含量温度上升幅度较小。从局部分析可看出，不同位置(上部距开封处20、60、100 cm，下部距开封处20、100 cm)的青贮饲料温度上升幅度也有差异。首先，对比上、下部分20、100 cm的温升幅度，上部较下部更高，这是由于裹包式青贮饲料长期水平放置导致其自重改变了容重的纵向分布，上部较下部饲料容重变小，孔隙度变大，氧气侵入更多导致产热更多，温度上升幅值更大[14]。其次，饲料上部不同位置(20、60、100 cm)的温度变化呈现出距离开封处越近温度上升越显著的规律，即温升幅度由浅至深依次递减，表明在开封处附近的氧气浓度以及微生物的分解产生的生化反应最为剧烈，随着氧气侵入到饲料深处，导致裹包内部的饲料有氧变质。最后，在20 cm处的氧气消耗完全后，该位置的温度不再继续上升，但在100 cm处的温度出现第2次上升，只是增速放缓，这是由于开封处附近的饲料变质过程产生的热量传导进入深层饲料所致，这也不利于未受氧胁迫的深层饲料的存储。

2.2　模型仿真值与测量值的比较

图4给出了原位观测与模型模拟的温度和氧气的动态变化过程，其中，图4a是干物质含量(DM-1)为39.0%的青贮牧草饲料温度和氧气变化，图4b是干物质含量(DM-2)为31.4%的青贮牧草饲料温度和氧气变化。从图中可以看出，模拟值的变化趋势整体上与观测值相符合，温升幅度均呈现由浅入深依次递减的规律，氧气变化也是先升后降最终稳定在某个值，但局部区域的温升与氧动态过程仍有差异。首先，如图4a所示，60 cm的温度预测值在2.5 d前出现低估，而在2.5 d后则出现高估。同时，图4b所示的20 cm处的预测值也出现类似情况。其可能原因在于酵母菌遇到氧气后被激活并开始快速繁殖，其繁殖速率受到温度、氧气等各种因素影响，各因素之间相互作用复杂，无法用简单的指数模型来模拟微生物的生长，因此出现了预测的早期生长速率偏低，与实际有差异。其次，图4a所示100 cm与图4b中60 cm和100 cm的上部温度预测值均出现低估，而在图4c所示100 cm下部和图4b所示20 cm下部的温度则出现明显的高估。其原因在于饲料本身内部的容重分布(即孔隙分布)不均也会对氧气扩散系数造成影响，从而影响到局部微生物的繁殖速率。在本研究中，实际的容重值比模型输入参数值低，则出现低估，反之则出现高估。最后，图4给出了两种干物质含量的青贮饲料氧气浓度在20 cm上部的预测值与观测值的对比。对于高干物质含量的青贮饲料(图4a)，氧气浓度在升降过程中两者吻合较好，到了平稳阶段，预测值出现了高估(<3%)；而对于低干物质含量的青贮饲料，预测值均出现高估，尤其在升降过程中预测误差较大(<9%)，同时与相同位置的温度预测值的低估(<4℃)相对应，其原因也是由裹包式青贮饲料内部容重分布不均造成。此外，在青贮饲料暴露在空气3 d后，观察到因有氧变质消耗氧气，产生了大量二氧化碳和水。产生的水填充了饲料中的部分孔隙，致使氧气扩散速率下降，多余水分因重力作用在饲料垂直方向流动，并渗出饲料底部，也很好地解释了图4中氧气浓度在平稳阶段均出现高估的原因。

图4　开封后5 d的原位观测数据与模型模拟数据对比Fig.4　Comparison between measurements and simulations for five days after opening plastic cover of bale silage

3　结束语

运用基于Clark电极原理的溶解氧传感器成功实现了高温高湿环境下裹包式青贮牧草饲料有氧胁迫过程内部氧动态与温度变化的原位观测。在此基础上，构建了基于偏微分-常微分方程系统的有氧胁迫下青贮饲料内部氧动态、温升以及微生物繁殖耦合模型，并检验了模型的适用性。通过分析对比实验原位观测值与模型预测值的差异发现：模型预测结果与实验测定结果总体上有较好的一致性，两者差异的原因可能是青贮饲料内部容重(孔隙)分布不均造成的。此外，有氧变质后期产生较多的水从而改变了气体扩散率、比热容等初始参数，导致后期预测精度降低。

1HAFNER S D, MONTES F, ROTZ C A, et al. Ethanol emission from loose corn silage and exposed silage particles[J]. Atmospheric Environment, 2010, 44(34): 4172-4180.

2BECHTEL H E, ATKESON F W, HUGHES J S. Brown silage from atlas sorgo-chemical composition, and apparent digestibility as determined by feeding to dairy cows[J]. Journal of Animal Science, 1943, 2(4): 295-303.

3MÜLLER C E, PAULY T M. Storage of small bale silage and haylage-influence of storage period on fermentation variables and microbial composition [J]. Grass and Forage Science, 2007, 62(3): 274-283.

4魏炳泉,马奎元. 青贮饲料的二次发酵与防止[J]. 畜禽业, 2014(6): 29-30.

5DRIEHUIS F. Silage and the safety and quality of dairy foods: a review[J]. Agricultural and Food Science, 2013, 22(1): 16-34.

7PITT R E, MUCK R E. A diffusion model of aerobic deterioration at the exposed face of bunker silos[J]. Journal of Agricultural Engineering Research, 1993, 55(1): 11-26.

8MUCK R E, PITT R E. Aerobic deterioration in corn silage relative to the silo face[J]. Transactions of the ASAE, 1994, 37(3): 735-743.

9RUXTON G D, GIBSON G J. A critical introduction to mathematical modelling of the aerobic deterioration of silage[J]. Journal of the Science of Food and Agriculture, 1993, 62(3): 209-217.

10RUXTON G D, GIBSON G J. A mathematical model of the aerobic deterioration of big-hale silage and its implications for the growth ofListeriamonocytogenes[J]. Grass and Forage Science, 1995, 50(4): 331-344.

11WILLIAMS A G, HOXEY R P, LOWE J F. Changes in temperature and silo gas composition during ensiling, storage and feeding-out grass silage[J]. Grass and Forage Science, 1997, 52(2): 176-189.

12MCENIRY J, FORRISTAL P D, O’KIELY P. Gas composition of baled grass silage as influenced by the amount, stretch, colour and type of plastic stretch-film used to wrap the bales, and by the frequency of bale handling[J]. Grass and Forage Science, 2011, 66(2): 277-289.

13SUN Y, LI M, CHENG Q, et al. Tracking oxygen and temperature dynamics in maize silage-novel application of a Clark oxygen electrode[J]. Biosystems Engineering, 2015, 139: 60-65.

14SUN Y, LI M, ZHOU H, et al. In situ measurements and simulation of oxygen diffusion and heat transfer in maize silage relative to the silo surface[J]. Computers and Electronics in Agriculture, 2017, 137: 1-8.

15程洪, 孙宇瑞, BUESCHER W, 等. 破损青黄贮裹包饲料内部温度动态研究[J/OL]. 农业机械学报, 2015, 46(4): 194-199.http:∥www.j-csam.org/ch/reader/view_abstract.aspx?file_no=20150429&flag=1&journal_id=jcsam.DOI:10.6041/j.issn.1000-1298.2015.04.029.

CHENG Hong, SUN Yurui, BUESCHER W, et al. Temperature dynamics in bale silage with damaged cover [J/OL]. Transactions of the Chinese Society for Agricultural Machinery, 2015, 46(4): 194-199. (in Chinese)

16BORREANI G, TABACCO E. The relationship of silage temperature with the microbiological status of the face of corn silage bunkers[J]. Journal of Dairy Science, 2010, 93(6): 2620-2629.

17VAN DERLINDEN E, VAN IMPE J F. Modeling growth rates as a function of temperature: model performance evaluation with focus on the suboptimal temperature range[J]. International Journal of Food Microbiology, 2012, 158(1): 73-78.

18SHAN G, SUN Y, LI M, et al. An assessment of three different in situ oxygen sensors for monitoring silage production and storage[J]. Sensors, 2016, 16(1): 91.

19高应才. 数学物理方程及其数值解法[M]. 北京: 高等教育出版社, 1983.

20MUCK R E. Factors influencing silage quality and their implications for management[J]. Journal of Dairy Science, 1988, 71(11): 2992-3002.