程　強(qiáng)　周海洋　單桂林　李夢華　BUESCHER W　孫宇瑞

(1.中國農(nóng)業(yè)大學(xué)信息與電氣工程學(xué)院， 北京 100083； 2.波恩大學(xué)農(nóng)業(yè)工程研究所， 波恩 53115)

0　引言

目前，國內(nèi)外使用較廣的青飼料貯存方式包括：圓柱形裹包式、方形裹包式、青貯槽式、塔倉式、窖倉式等。在這幾種方式中，無論從減少收獲過程中青飼料損失考慮，還是從保證發(fā)酵質(zhì)量與安全存儲等方面綜合衡量，裹包式青貯飼料都是最優(yōu)的選擇。裹包式青貯飼料在西方發(fā)達(dá)國家已經(jīng)非常普遍，我國也在日漸普及過程中。然而，裹包式貯存也有一些缺陷，比如，因其自身表面積/體積比很大，在運(yùn)輸過程中塑封膜極易破損。氧氣的侵入破壞了無氧貯存環(huán)境，各類好氧微生物大量繁殖，消耗飼料的營養(yǎng)成分，產(chǎn)生熱量，引起溫度異常上升，伴隨著霉菌的滋生，導(dǎo)致青貯飼料腐敗變質(zhì)，產(chǎn)生有害的揮發(fā)物對環(huán)境造成影響[1]。在變質(zhì)過程中，氧氣與溫度的變化是衡量各類好氧微生物活動的關(guān)鍵參數(shù)，也是反映與評價青貯飼料貯藏質(zhì)量的重要指標(biāo)[2-6]。

國際學(xué)者構(gòu)建了氧氣與溫度動態(tài)模型用于預(yù)測青貯飼料有氧脅迫下的變質(zhì)過程，并借助實(shí)驗(yàn)觀測方法檢驗(yàn)?zāi)Ｐ皖A(yù)測的準(zhǔn)確性與適用性。美國學(xué)者PITT等[7]和MUCK等[8]建立了槽式青貯玉米飼料有氧變質(zhì)過程的數(shù)學(xué)模型。隨后，使用破壞性取樣法抽取氣體，再用氣相色譜儀測定青貯玉米飼料內(nèi)部氧氣濃度變化，并用這些觀測數(shù)據(jù)檢驗(yàn)所構(gòu)建的模型的適用性。RUXTON等[9-10]建立了裹包式青貯牧草飼料有氧變質(zhì)過程的數(shù)學(xué)模型，WILLIAMS等[11]對該模型進(jìn)行了驗(yàn)證。在后續(xù)研究中，其他學(xué)者運(yùn)用類似方法測量有氧脅迫下青貯飼料內(nèi)部氧氣與溫度變化[12]，但這些觀測方法有一個共同的缺陷：抽取氣體樣本將會擾動飼料內(nèi)部的氣體分布，造成實(shí)驗(yàn)結(jié)果出現(xiàn)偏差。此外，抽氣分析法每天只能取少數(shù)幾個氣體樣本，無法實(shí)時記錄氧氣濃度變化，難以突破連續(xù)觀測技術(shù)瓶頸，致使青貯飼料有氧變質(zhì)風(fēng)險預(yù)測的研究進(jìn)展緩慢。隨著科學(xué)技術(shù)的不斷進(jìn)步，SUN等[13]運(yùn)用一種基于Clark氧電極原理的溶解氧傳感器成功實(shí)現(xiàn)了青貯玉米飼料內(nèi)部氧氣與溫度的原位實(shí)時觀測，并對PITT等[7]提出的預(yù)測模型進(jìn)行了評價[14]。而對于青貯牧草，其成分、顆粒大小、干質(zhì)量、菌落類型與數(shù)量等理化生性質(zhì)與青貯玉米差異明顯，基于青貯玉米對象提出的原位觀測方法與預(yù)測模型是否適用于青貯牧草還有待驗(yàn)證。程洪等[15]雖然測量了氧脅迫下裹包式青貯牧草飼料內(nèi)部溫升分布，但缺少氧氣濃度的測量數(shù)據(jù)，也未涉及模型檢驗(yàn)。為此，本文嘗試運(yùn)用Clark氧電極傳感器原位觀測裹包式青貯牧草內(nèi)部氧氣與溫度變化，建立有氧變質(zhì)過程中青貯牧草內(nèi)部氧動態(tài)、溫升以及微生物繁殖耦合模型。由于裹包式青貯飼料的對稱性，將三維模型簡化成二維模型，最后根據(jù)實(shí)驗(yàn)觀測數(shù)據(jù)與模型預(yù)測數(shù)據(jù)的對比結(jié)果評價模型的適用性。

1　材料與方法

1.1　模型建立

青貯飼料有氧變質(zhì)過程可包括氧氣擴(kuò)散、好氧微生物生長、營養(yǎng)物質(zhì)分解、生化反應(yīng)產(chǎn)熱與傳導(dǎo)。當(dāng)裹包青貯的塑封膜破損或青貯飼料開窖后，氧氣侵入，根據(jù)菲克擴(kuò)散定律，氧氣在青貯飼料中的擴(kuò)散過程可運(yùn)用二維拋物線型偏微分方程來描述[7,10]

(1)

式中ψ——氧氣濃度，%

D——氧氣在空氣中的擴(kuò)散率，0℃時D=0.064 1 m2/h

Φ——飼料孔隙度

ρ——飼料容重，kg/m3

τ——多孔介質(zhì)中孔隙的彎曲度，與孔隙度(容重)有關(guān)，本文取常數(shù)0.667

Yg——實(shí)際增長率，取0.5

Rs——用于呼吸作用的氧氣體積分?jǐn)?shù)，取0.747 m3/kg

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

μg——微生物增長率，h-1，受溫度和氧氣濃度影響

式(1)右邊第2項(xiàng)為源/匯項(xiàng)，反映了微生物分解營養(yǎng)物質(zhì)消耗的氧氣，這一過程會產(chǎn)生熱量，熱量的二維傳導(dǎo)過程可描述為

(2)

其中

cp=cp,w(1.451dm+1)

(3)

式中T——青貯飼料溫度，K

cp——青貯飼料比熱容，J/(kg·K)

k——青貯飼料熱導(dǎo)率，W/(m·K)

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

α——熱量轉(zhuǎn)化系數(shù)

dm——干物質(zhì)含量，%

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

式(2)右邊第2項(xiàng)為源/匯項(xiàng)，反映了微生物在有氧環(huán)境下分解營養(yǎng)物質(zhì)產(chǎn)生的熱量，用化學(xué)方程式描述

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

(4)

模型參數(shù)的測定對于模型預(yù)測精度非常重要。式(1)中微生物生長模型可描述為

(5)

其中

μg=μgmaxfTfψ

(6)

(7)

式中μg-μm——微生物凈生長率

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

μm——酵母菌實(shí)際維持率

fT——微生物生長率隨溫度變化的函數(shù)[16-17]

Km——常數(shù)，取0.001 45[7]

ψ0——空氣中的氧氣濃度，取21%

fψ——微生物生長率隨氧氣濃度變化的函數(shù)

由于酵母菌的繁殖是青貯飼料有氧變質(zhì)初期的主導(dǎo)因素，因此本文主要模擬變質(zhì)初期酵母菌繁殖過程。營養(yǎng)成分的消耗(L，%)計(jì)算公式為

(8)

飼料熱導(dǎo)率k、孔隙度Φ和最大容重ρmax為

k=kw(1-Φ)+k0Φ

(9)

(10)

(11)

式中kw——水熱導(dǎo)率k0——空氣熱導(dǎo)率模型求解過程詳見文獻(xiàn)[7,10,14]。

1.2　Clark氧電極基本原理

薄膜氧電極是為測量水中所溶解的氧氣量而研制的極譜電極，其本質(zhì)上是一個電化學(xué)電池。該電極由Clark于1953年最先研制的，也稱為Clark氧電極，其結(jié)構(gòu)如圖1所示。Clark氧電極由嵌在不導(dǎo)電原料上的銀電極和鉑電極組成。其中，銀極作為陽極或參比電極制成圓環(huán)形狀，面積盡可能大一些，以降低電極表面電流密度，減少陽極的極化現(xiàn)象，使其電極電位不受外加電壓的影響。鉑極作為陰極或指示電極，通常是圓點(diǎn)狀，置于銀電極的中間，電化學(xué)反應(yīng)就產(chǎn)生在鉑電極上。用15～20 μm的聚四氟乙烯薄膜包裹住氧電極的表層，把電解質(zhì)(飽和氯化鉀溶液)充入氧電極與薄膜之間。因?yàn)樗腥芙獾难鯕饽芡ㄟ^乙烯膜而電解質(zhì)溶液不能進(jìn)入，從而消除被測溶液中存在的各類離子發(fā)生反應(yīng)的影響。

根據(jù)亨利定律，電解質(zhì)溶液中的溶解氧含量與氣液界面處的氧分壓呈比例關(guān)系[13]，即

p=kHCDO

(12)

式中p——待測氣體中的氧分壓，Pa

CDO——電解質(zhì)溶液中的溶解氧含量，mol/L

kH——亨利常數(shù)，對于氧氣而言，在298 K下取770 Pa·L/mol

由于亨利常數(shù)受溫度影響，因此需要借助van’t Hoff方程進(jìn)行校正

(13)

式中，C1值為1 700 K(對于氧氣而言)。將式(13)代入式(12)中即可得到待測氣體/液體中的氧分壓。本文使用商業(yè)化的Clark氧電極(型號：DO-10，測量范圍：0～20 mg/L，電極輸出范圍：0～2.5 V，精度：±3%滿量程；氧電極內(nèi)部集成熱敏電阻：MF5，范圍0～50℃，分辨率：0.1℃，北京聯(lián)創(chuàng)思源測控技術(shù)有限公司)觀測裹包式青貯飼料內(nèi)部高溫高濕環(huán)境下氧氣與溫度的變化。關(guān)于氧電極的標(biāo)定方法以及二氧化碳對氧電極的輸出影響詳見文獻(xiàn)[18]。

圖1　Clark氧電極內(nèi)部示意圖Fig.1　Schematic diagram of Clark oxygen electrodes1.電流計(jì)　2.氯化鉀電解液　3.橡膠圈　4.聚乙烯或聚四氟乙烯薄膜　5.鉑極　6.銀極　7.電源

1.3　實(shí)驗(yàn)方案設(shè)計(jì)

圖2　實(shí)驗(yàn)設(shè)計(jì)與傳感器安裝位置Fig.2　Experiment design and sensor location

如圖2所示，所測試的裹包式青貯飼料由德國CLAAS公司的Rollant455型打包機(jī)打包成圓柱形，底面直徑約120 cm，高約120 cm，飼料種類為黑麥草。裹包式青貯飼料樣本在打包時打包機(jī)參數(shù)設(shè)置相同，其粒徑大小、體積與鮮容重非常接近，最大限度消除了其他物理因素對實(shí)驗(yàn)結(jié)果的影響。首先將氧電極和溫度傳感器(DS18B20型)原位植入未開封的青貯裹包飼料內(nèi)部[14]，探頭離飼料表面深度均為20 cm。溫度傳感器分別記錄距離破損表面約20、60、100 cm處的數(shù)據(jù)，氧傳感器記錄距離破損表面約20 cm處的數(shù)據(jù)。傳感器植入后立即封住植入?yún)^(qū)域，以確保青貯飼料內(nèi)部的無氧環(huán)境。待青貯飼料內(nèi)部因安裝傳感器而引入的微量氧氣消耗完之后，將青貯飼料圓柱面一側(cè)開封，使其暴露在空氣中，模擬破損或飼喂過程空氣侵入方式，數(shù)據(jù)采集器(自制多通道采集器，微處理器：MSP430F149，14位A/D，Texas Instruments, Dallas, Texas, USA)記錄開封之后氧氣、飼料不同位置溫度以及環(huán)境溫度的變化。

圖3　開封后5 d內(nèi)的觀測數(shù)據(jù)Fig.3　Observation data of five days after opening plastic cover of bale silage

同時對2種不同干物質(zhì)含量(DM-1和DM-2)的青貯牧草進(jìn)行監(jiān)測，每個水平測試3個裹包，觀測數(shù)據(jù)取平均值，共6個裹包飼料樣本。實(shí)驗(yàn)在室內(nèi)環(huán)境下進(jìn)行，共經(jīng)歷了5 d，室溫保持在(27.3±1.1)℃，風(fēng)速可以忽略。在實(shí)驗(yàn)開始階段，每個裹包用取樣器收集N個樣本(N≥3)，每個樣本質(zhì)量大于30 g，用于分析微生物種類及含量初值，同時確定不同位置飼料鮮容重，具體分析方法詳見文獻(xiàn)[14]，分析結(jié)果如表1所示。由于裹包式青貯飼料長期水平放置，其自重改變了容重的縱向分布，上部飼料(DM-1：381 kg/m3；DM-2：402 kg/m3)較下部飼料容重變小(DM-1：597 kg/m3；DM-2：643 kg/m3)。模型仿真算法采用的是非穩(wěn)態(tài)顯式差分法[19]，最小間距設(shè)定為0.1 m，最小時間間隔為0.067 h，在Matlab(Mathworks, Inc., Version 7.0)開發(fā)平臺進(jìn)行仿真計(jì)算。

表1　青貯飼料參數(shù)分析結(jié)果Tab.1　Analysis of related parameters in silage

2　結(jié)果與討論

2.1　氧氣與溫度原位觀測動態(tài)變化與分析

裹包式青貯飼料開封后氧氣侵入與溫度變化過程如圖3所示(圖中顯示的數(shù)據(jù)曲線是3個樣本的平均值)。從整體上看，氧動態(tài)變化過程中峰值的出現(xiàn)表明氧氣擴(kuò)散與微生物消耗達(dá)到了瞬態(tài)平衡，在峰值之前氧氣擴(kuò)散占據(jù)主導(dǎo)，之后微生物耗氧占據(jù)優(yōu)勢，達(dá)到動態(tài)平衡后穩(wěn)定在某個值[7,14]。圖3a反映了高干物質(zhì)含量飼料樣本的觀測數(shù)據(jù)，圖3b是低干物質(zhì)含量飼料樣本的觀測數(shù)據(jù)。相比之下不難看出：首先，氧氣達(dá)到的峰值有些差異，高干物質(zhì)含量飼料的峰值(15.1%)略大于低干物質(zhì)含量飼料的峰值(14.2%)。其次在耗氧速率方面，高干物質(zhì)含量的飼料的耗氧速率(約2.5 d)僅僅略低于低干物質(zhì)含量(約2.3 d)，整體差異不顯著。最后，在耗氧初始階段，高干物質(zhì)含量的飼料樣本耗氧速率較慢，大約經(jīng)過1.2 d后斜率變大，耗氧速率開始加快；而低干物質(zhì)含量的飼料樣本在開始階段耗氧速率較快，經(jīng)過0.5 d后斜率變小，速率降低。其可能原因在于低干物質(zhì)含量的青貯飼料濕度較高，更易于酵母菌的存活與繁殖，并且水分活性及可溶性碳水化合物(WSC)含量也高，因而耗氧速率較高[20]。

高干物質(zhì)含量飼料溫度增幅整體上比低干物質(zhì)含量飼料的溫度增幅高。其原因在于低干物質(zhì)含量的飼料水分含量較高，本身的比熱容也較大，在釋放同等數(shù)量級熱量的條件下，低干物質(zhì)含量溫度上升幅度較小。從局部分析可看出，不同位置(上部距開封處20、60、100 cm，下部距開封處20、100 cm)的青貯飼料溫度上升幅度也有差異。首先，對比上、下部分20、100 cm的溫升幅度，上部較下部更高，這是由于裹包式青貯飼料長期水平放置導(dǎo)致其自重改變了容重的縱向分布，上部較下部飼料容重變小，孔隙度變大，氧氣侵入更多導(dǎo)致產(chǎn)熱更多，溫度上升幅值更大[14]。其次，飼料上部不同位置(20、60、100 cm)的溫度變化呈現(xiàn)出距離開封處越近溫度上升越顯著的規(guī)律，即溫升幅度由淺至深依次遞減，表明在開封處附近的氧氣濃度以及微生物的分解產(chǎn)生的生化反應(yīng)最為劇烈，隨著氧氣侵入到飼料深處，導(dǎo)致裹包內(nèi)部的飼料有氧變質(zhì)。最后，在20 cm處的氧氣消耗完全后，該位置的溫度不再繼續(xù)上升，但在100 cm處的溫度出現(xiàn)第2次上升，只是增速放緩，這是由于開封處附近的飼料變質(zhì)過程產(chǎn)生的熱量傳導(dǎo)進(jìn)入深層飼料所致，這也不利于未受氧脅迫的深層飼料的存儲。

2.2　模型仿真值與測量值的比較

圖4給出了原位觀測與模型模擬的溫度和氧氣的動態(tài)變化過程，其中，圖4a是干物質(zhì)含量(DM-1)為39.0%的青貯牧草飼料溫度和氧氣變化，圖4b是干物質(zhì)含量(DM-2)為31.4%的青貯牧草飼料溫度和氧氣變化。從圖中可以看出，模擬值的變化趨勢整體上與觀測值相符合，溫升幅度均呈現(xiàn)由淺入深依次遞減的規(guī)律，氧氣變化也是先升后降最終穩(wěn)定在某個值，但局部區(qū)域的溫升與氧動態(tài)過程仍有差異。首先，如圖4a所示，60 cm的溫度預(yù)測值在2.5 d前出現(xiàn)低估，而在2.5 d后則出現(xiàn)高估。同時，圖4b所示的20 cm處的預(yù)測值也出現(xiàn)類似情況。其可能原因在于酵母菌遇到氧氣后被激活并開始快速繁殖，其繁殖速率受到溫度、氧氣等各種因素影響，各因素之間相互作用復(fù)雜，無法用簡單的指數(shù)模型來模擬微生物的生長，因此出現(xiàn)了預(yù)測的早期生長速率偏低，與實(shí)際有差異。其次，圖4a所示100 cm與圖4b中60 cm和100 cm的上部溫度預(yù)測值均出現(xiàn)低估，而在圖4c所示100 cm下部和圖4b所示20 cm下部的溫度則出現(xiàn)明顯的高估。其原因在于飼料本身內(nèi)部的容重分布(即孔隙分布)不均也會對氧氣擴(kuò)散系數(shù)造成影響，從而影響到局部微生物的繁殖速率。在本研究中，實(shí)際的容重值比模型輸入?yún)?shù)值低，則出現(xiàn)低估，反之則出現(xiàn)高估。最后，圖4給出了兩種干物質(zhì)含量的青貯飼料氧氣濃度在20 cm上部的預(yù)測值與觀測值的對比。對于高干物質(zhì)含量的青貯飼料(圖4a)，氧氣濃度在升降過程中兩者吻合較好，到了平穩(wěn)階段，預(yù)測值出現(xiàn)了高估(<3%)；而對于低干物質(zhì)含量的青貯飼料，預(yù)測值均出現(xiàn)高估，尤其在升降過程中預(yù)測誤差較大(<9%)，同時與相同位置的溫度預(yù)測值的低估(<4℃)相對應(yīng)，其原因也是由裹包式青貯飼料內(nèi)部容重分布不均造成。此外，在青貯飼料暴露在空氣3 d后，觀察到因有氧變質(zhì)消耗氧氣，產(chǎn)生了大量二氧化碳和水。產(chǎn)生的水填充了飼料中的部分孔隙，致使氧氣擴(kuò)散速率下降，多余水分因重力作用在飼料垂直方向流動，并滲出飼料底部，也很好地解釋了圖4中氧氣濃度在平穩(wěn)階段均出現(xiàn)高估的原因。

圖4　開封后5 d的原位觀測數(shù)據(jù)與模型模擬數(shù)據(jù)對比Fig.4　Comparison between measurements and simulations for five days after opening plastic cover of bale silage

3　結(jié)束語

運(yùn)用基于Clark電極原理的溶解氧傳感器成功實(shí)現(xiàn)了高溫高濕環(huán)境下裹包式青貯牧草飼料有氧脅迫過程內(nèi)部氧動態(tài)與溫度變化的原位觀測。在此基礎(chǔ)上，構(gòu)建了基于偏微分-常微分方程系統(tǒng)的有氧脅迫下青貯飼料內(nèi)部氧動態(tài)、溫升以及微生物繁殖耦合模型，并檢驗(yàn)了模型的適用性。通過分析對比實(shí)驗(yàn)原位觀測值與模型預(yù)測值的差異發(fā)現(xiàn)：模型預(yù)測結(jié)果與實(shí)驗(yàn)測定結(jié)果總體上有較好的一致性，兩者差異的原因可能是青貯飼料內(nèi)部容重(孔隙)分布不均造成的。此外，有氧變質(zhì)后期產(chǎn)生較多的水從而改變了氣體擴(kuò)散率、比熱容等初始參數(shù)，導(dǎo)致后期預(yù)測精度降低。

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魏炳泉,馬奎元. 青貯飼料的二次發(fā)酵與防止[J]. 畜禽業(yè), 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, 等. 破損青黃貯裹包飼料內(nèi)部溫度動態(tài)研究[J/OL]. 農(nóng)業(yè)機(jī)械學(xué)報, 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高應(yīng)才. 數(shù)學(xué)物理方程及其數(shù)值解法[M]. 北京: 高等教育出版社, 1983.

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