基于隨機(jī)森林的寒區(qū)奶牛舍環(huán)境因素與產(chǎn)奶量關(guān)系研究

丁 濤 蔣雷生,2 施正香,3 趙 洋,2 馬 慧

(1.中國農(nóng)業(yè)大學(xué)水利與土木工程學(xué)院， 北京 100083； 2.北京市供水管網(wǎng)與安全節(jié)能中心， 北京 100083；3.農(nóng)業(yè)農(nóng)村部設(shè)施農(nóng)業(yè)工程重點(diǎn)實(shí)驗(yàn)室， 北京 100083； 4.北京首農(nóng)畜牧發(fā)展有限公司， 北京 100083)

0 引言

奶牛舍環(huán)境對奶牛的生產(chǎn)性能具有決定性影響，當(dāng)牛舍小氣候變化幅度超出奶牛的適應(yīng)范圍時，奶牛的生產(chǎn)力和健康便會受到負(fù)面影響[1]。目前，評價(jià)牛舍環(huán)境的指標(biāo)大部分是基于熱環(huán)境的熱應(yīng)激評價(jià)方法[2]，如THI[3]、BGHI[4]、THVI[5]、HLI[6]以及不同的體感溫度表達(dá)式[7-9]，這些指數(shù)關(guān)注了溫度、相對濕度、風(fēng)速和熱輻射，卻忽略了光照強(qiáng)度及空氣質(zhì)量的影響，不適用于非熱應(yīng)激季節(jié)的奶牛舍環(huán)境評估。由于季節(jié)特殊性，冬季舍內(nèi)風(fēng)速低、有害氣體濃度高，特別是二氧化碳濃度和氨氣濃度對奶牛的健康和生產(chǎn)有重要影響。研究表明[10-12]，光照使產(chǎn)奶量最大化的明暗交替最佳節(jié)點(diǎn)是16 h光照，8 h黑暗，合理的光照能提高產(chǎn)奶量8%～10%、促進(jìn)產(chǎn)后發(fā)情和受孕。然而針對寒區(qū)的氣候特點(diǎn)，光照對于奶牛生產(chǎn)的影響沒有更加深入地研究。

隨機(jī)森林算法是數(shù)據(jù)挖掘分類算法的重要組成部分，在生物信息學(xué)、生態(tài)學(xué)、醫(yī)學(xué)等方面都被廣泛應(yīng)用[13-14]。針對畜禽舍環(huán)境的調(diào)控復(fù)雜性，國內(nèi)外已經(jīng)有采用模糊算法、神經(jīng)網(wǎng)絡(luò)等手段實(shí)現(xiàn)舍內(nèi)環(huán)境的智能調(diào)節(jié)[15-17]。有學(xué)者利用奶牛數(shù)據(jù)和機(jī)器學(xué)習(xí)算法預(yù)測奶牛泌乳早期的代謝狀態(tài)，其中隨機(jī)森林算法和支持向量機(jī)的表現(xiàn)最優(yōu)[18]。本文通過對奶牛舍環(huán)境參數(shù)以及產(chǎn)奶量的連續(xù)監(jiān)測，從不同角度分析其中的變化規(guī)律，建立基于環(huán)境因素的隨機(jī)森林回歸預(yù)測模型，并驗(yàn)證其在寒區(qū)奶牛舍的適用性。

1 材料與方法

1.1 數(shù)據(jù)集

試驗(yàn)于2018年9月至2019年3月在吉林省白城市首都農(nóng)業(yè)集團(tuán)白城牛場進(jìn)行。牛舍東西長度約250 m，南北跨度約27 m，屋頂為彩鋼夾芯板雙坡屋頂，奶牛采用全舍飼散欄飼養(yǎng)模式，舍內(nèi)奶牛約580頭，每日擠奶3班。牛舍四周窗戶除冬季外均開啟，屋頂小窗常年開啟，南面屋頂上每隔12 m均勻布置有采光板增加舍內(nèi)采光，牛舍內(nèi)無噴淋降溫設(shè)施。舍內(nèi)沿東西方向共有4排工字鋼，東西方向相鄰工字鋼間隔6 m。在中間兩排工字鋼上間隔12 m安裝36寸風(fēng)機(jī)，在躺臥區(qū)域上方間隔15 m安裝54寸風(fēng)機(jī)，舍內(nèi)共安裝有40臺36寸風(fēng)機(jī)、32臺54寸風(fēng)機(jī)。風(fēng)機(jī)在夏季全部24 h運(yùn)行，秋季白天部分風(fēng)機(jī)在高溫時段對稱開啟，晚上關(guān)閉。風(fēng)機(jī)運(yùn)行個數(shù)隨溫度降低而減少直至全部停止運(yùn)行，試驗(yàn)期間舍內(nèi)大部分時間都處于風(fēng)速較低狀態(tài)。

試驗(yàn)牛舍由于跨度不大、采光板的設(shè)置及非熱應(yīng)激時期太陽高度角較小，前期通過手持testo540型光照度測試儀測試奶牛活動區(qū)域，結(jié)果表明南北光照強(qiáng)度差異不大。牛舍光照為自然光照，舍內(nèi)在高度4 m處間隔6 m安裝15 W節(jié)能燈，共安裝80個節(jié)能燈，如圖1所示，夜間開啟節(jié)能燈便于人員工作。為了避免距離風(fēng)機(jī)過近數(shù)據(jù)不準(zhǔn)確，因此將傳感器安裝在相鄰兩個風(fēng)機(jī)中間的工字鋼上，距離地面1.8 m。舍內(nèi)各項(xiàng)設(shè)施布置均勻，舍外附近無其余建筑物，因此布置4個環(huán)境數(shù)據(jù)采集點(diǎn)可以較好地反映舍內(nèi)環(huán)境水平。測點(diǎn)分別距離奶牛舍東門、西門45 m及90 m，舍外布置1個采集點(diǎn)E，位于奶牛舍舍外東南角，如圖2所示。每個采集點(diǎn)可以采集溫度、相對濕度、風(fēng)速、二氧化碳濃度、氨氣濃度、光照強(qiáng)度。傳感器為武漢中科能慧公司生產(chǎn)，采樣間隔為10 min，具體參數(shù)見表1。選擇泌乳性能相近的30頭奶牛每日產(chǎn)奶量進(jìn)行分析。

1.2 數(shù)據(jù)處理

采集的數(shù)據(jù)往往存在部分缺失、不正確、含有噪聲或其他類型的不一致現(xiàn)象[19]。采用刪除法與插補(bǔ)法處理數(shù)據(jù)缺失的情況，同時利用Matlab進(jìn)行小波去噪以減少噪聲對數(shù)據(jù)的影響。考慮到環(huán)境因素的變化達(dá)到一定的閾值才會對畜禽的生產(chǎn)能力產(chǎn)生影響[20]，為了確定當(dāng)日環(huán)境數(shù)據(jù)的代表值，本文將數(shù)據(jù)從小到大進(jìn)行排列取不同百分位值及均值，研究其與產(chǎn)奶量的相關(guān)關(guān)系。

表1 試驗(yàn)傳感器參數(shù)Tab.1 Test sensor parameter list

1.3 機(jī)器學(xué)習(xí)

隨機(jī)森林算法通過自助法在原數(shù)據(jù)集中有放回地隨機(jī)抽取n個樣本集，每次抽取2/3樣本容量的數(shù)據(jù)作為袋內(nèi)數(shù)據(jù)，建立m棵決策樹構(gòu)建隨機(jī)森林，計(jì)算m棵決策樹回歸結(jié)果的平均值進(jìn)行預(yù)測[21]。該方法由于引入“雙隨機(jī)”的思想與組合投票取均值原則，因此對噪聲和異常值不敏感、不易陷入過擬合[22]。

1.3.1數(shù)據(jù)預(yù)處理

構(gòu)建機(jī)器學(xué)習(xí)模型時，使用樣本的平均誤差最小化作為模型的學(xué)習(xí)準(zhǔn)則，則少數(shù)類樣本的誤差在總體誤差中的比重會偏小，使得模型對少數(shù)類樣本的預(yù)測性能下降。本文采用合成少數(shù)類過采樣技術(shù)方法，利用隨機(jī)線性插值法，在兩個距離較近的少數(shù)樣本間創(chuàng)造新的少數(shù)類樣本來平衡樣本。由于各環(huán)境因子的取值方式對產(chǎn)奶量的影響程度尚不明確，因此將各環(huán)境因子的每日平均值、每日12個最大值平均、每日12個最小值平均以及第10、20、35、50、65、80、90百分位數(shù)組成基礎(chǔ)自變量集X，將產(chǎn)奶量作為基礎(chǔ)因變量集Y。對自變量和因變量分別采取平方根變換、平方變換、立方變換、自然對數(shù)變換、指數(shù)變換、倒數(shù)變換。經(jīng)過變換后形成7個訓(xùn)練集，每個訓(xùn)練集包含420個自變量。

1.3.2模型評價(jià)指標(biāo)

模型訓(xùn)練的目標(biāo)是實(shí)現(xiàn)對訓(xùn)練樣本的最佳擬合，是基于誤差最小化原則，確定模型最優(yōu)參數(shù)與結(jié)構(gòu)的過程。選用決定系數(shù)(R2)和均方根誤差(E)作為模型評價(jià)指標(biāo)，同時將兩者合并為指標(biāo)P以便對比模型性能，其計(jì)算公式為

(1)

E0——剔除某自變量前模型訓(xùn)練的均方根誤差

Ei——剔除某自變量后模型訓(xùn)練的均方根誤差

1.3.3模型參數(shù)優(yōu)選

隨機(jī)森林回歸模型的構(gòu)建需要確定決策樹數(shù)量，對420個自變量篩選降維，確定最佳訓(xùn)練集。根據(jù)各模型的袋外數(shù)據(jù)均方誤差估計(jì)隨決策樹數(shù)量增長的變化曲線確定決策樹數(shù)量為120棵[23]。基于袋外數(shù)據(jù)估計(jì)的特征重要性，逐步剔除重要性較小的特征，并迭代訓(xùn)練隨機(jī)森林模型，直到自變量的數(shù)量降低到目標(biāo)個數(shù)。最后利用R2和E對優(yōu)選后的模型進(jìn)行評價(jià)。

2 結(jié)果與討論

2.1 環(huán)境因素日均值與產(chǎn)奶量之間的關(guān)系

試驗(yàn)期間每月的數(shù)據(jù)變化情況見表2,目前研究表明奶牛舒適溫度區(qū)間為5～25℃，可以看出在9—11月舍內(nèi)溫度處于奶牛舒適的溫度范圍內(nèi)時，相對濕度、風(fēng)速、光照強(qiáng)度、二氧化碳濃度、氨氣濃度、產(chǎn)奶量均存在顯著性差異。9月和次年1月產(chǎn)奶量近似，期間除二氧化碳濃度外，其他因素均顯著變化。2月和3月產(chǎn)奶量近似，除二氧化碳濃度與風(fēng)速外，其他因素均存在顯著性差異。表明在部分時段里二氧化碳濃度不存在顯著性差異時，其他環(huán)境因素的變化對產(chǎn)奶量沒有產(chǎn)生明顯地促進(jìn)或抑制作用。

為了減少自變量從而準(zhǔn)確分析其他因素與產(chǎn)奶量的規(guī)律，引入溫濕指數(shù)(Temperature humidity index,THI)，溫濕指數(shù)計(jì)算式為

THI=0.8T+(H/100)(T-14.4)+46.4

(2)

式中T——干球溫度H——相對濕度

溫濕指數(shù)與產(chǎn)奶量的關(guān)系如圖3所示。由圖3可知，THI與產(chǎn)奶量間并沒有明顯的變化規(guī)律。分析其原因，在9、10月舍內(nèi)風(fēng)機(jī)隨溫度變化決定開啟數(shù)量，風(fēng)機(jī)運(yùn)行影響了舍內(nèi)氣體分布。當(dāng)溫度較低時，牛舍內(nèi)風(fēng)速處于較低水平，有害氣體濃度增加，光照強(qiáng)度也處于較高的水平。不同環(huán)境因素在不同時間段對奶牛的影響不同，因此在同一THI水平及不同THI水平時，產(chǎn)奶量動態(tài)變化。為了確定不同區(qū)間內(nèi)影響產(chǎn)奶量的重要因素，根據(jù)產(chǎn)奶量及THI劃分為不同階段進(jìn)行分析，結(jié)果見表3。

表2 試驗(yàn)期間各月環(huán)境數(shù)據(jù)Tab.2 Monthly environmental data during test period

由表3可知，在THI處于較低和中等水平時，高產(chǎn)奶量和中產(chǎn)奶量間光照強(qiáng)度、二氧化碳濃度存在顯著差異。中產(chǎn)奶量和低產(chǎn)奶量間光照強(qiáng)度不存在明顯差異，二氧化碳濃度較高時產(chǎn)奶量低。THI處于較高水平時，高產(chǎn)奶量和低產(chǎn)奶量樣本間光照強(qiáng)度存在顯著性差異，二氧化碳濃度不存在顯著性差異。表明當(dāng)二氧化碳濃度較低時產(chǎn)奶量較高，光照對產(chǎn)奶量有進(jìn)一步的促進(jìn)作用。

表3 不同THI與產(chǎn)奶量階段時光照強(qiáng)度與二氧化碳濃度關(guān)系Tab.3 Relationship between illumination intensity and carbon dioxide concentration at different THI and milk yield stage

光照能抑制褪黑激素的分泌，其分泌的交替變化能引發(fā)一系列的激素反應(yīng)，從而提升產(chǎn)奶量[24]。PETERS等[25]研究表明接受16 h光照比9～12 h光照奶牛產(chǎn)奶量高6.7%。ESPINOZA等[26]研究表明保持16 h光照相比8 h光照奶牛產(chǎn)奶量高2.2 kg/d。REKSEN等[27]研究結(jié)果表明在高緯度地區(qū)，夜晚接受光照比未接受光照的奶牛產(chǎn)奶量高0.5 kg/頭。MUTHURAMALINGAM等[28]研究了不同光照強(qiáng)度對褪黑激素的影響水平，當(dāng)夜晚光照強(qiáng)度為50 lx時，褪黑激素分泌會受到抑制。目前普遍認(rèn)為光照強(qiáng)度200 lx有助于生產(chǎn)，各類研究明確了光照對產(chǎn)奶量的促進(jìn)作用。由于各環(huán)境因素對奶牛的影響是動態(tài)的，本次試驗(yàn)主要分析正常生產(chǎn)狀況下多個環(huán)境因素與產(chǎn)奶量之間關(guān)系，確定不同階段各環(huán)境因素對產(chǎn)奶量的影響程度。夜間光照只用于照明，光照強(qiáng)度包含在每日環(huán)境因素日均值中，不影響本文分析結(jié)果。

東北地區(qū)晝夜溫差大，冬季溫度較低，牛舍的封閉措施導(dǎo)致二氧化碳濃度過高，繼而引起氧氣濃度降低，造成奶牛缺氧[29]。梅瑋等[30]研究表明牛舍中二氧化碳和氨氣濃度過高將不利于奶牛健康和生產(chǎn)性能。目前關(guān)于牛舍內(nèi)二氧化碳濃度的研究多是關(guān)于氣體分布規(guī)律方面的，包括不同牛舍類型、不同季節(jié)等方面，并以二氧化碳濃度1 500 mg/m3為標(biāo)準(zhǔn)進(jìn)行舍內(nèi)環(huán)境評價(jià)[31-33]。

考慮到隨著溫度的改變，人工管理會影響舍內(nèi)環(huán)境因素的分布規(guī)律。將觀測樣本劃分為高風(fēng)速與低風(fēng)速時各環(huán)境因素的相關(guān)系數(shù)見表4。可以看出，當(dāng)溫度較高時，風(fēng)機(jī)運(yùn)行的個數(shù)與舍內(nèi)溫度有關(guān)，因此風(fēng)速與溫度、相對濕度、二氧化碳濃度均存在明顯的相關(guān)性，各因素與產(chǎn)奶量間存在相關(guān)性，但不顯著。溫度與產(chǎn)奶量相關(guān)系數(shù)較高，此時評價(jià)舍內(nèi)環(huán)境時應(yīng)以溫度為主要因素。當(dāng)溫度較低時，舍內(nèi)風(fēng)速一直處于較為穩(wěn)定的低水平，溫濕度呈明顯的線性規(guī)律。二氧化碳濃度與產(chǎn)奶量相關(guān)系數(shù)為-0.532，光照強(qiáng)度與產(chǎn)奶量相關(guān)系數(shù)為0.720，二者均達(dá)到顯著水平。溫度、光照強(qiáng)度、二氧化碳濃度與產(chǎn)奶量的顏色映射圖如圖4所示，可以看出光照強(qiáng)度日均值250 lx與二氧化碳濃度日均值8×10-4可以明顯地劃分高低產(chǎn)奶量。因此建議在低溫時評價(jià)奶牛舍環(huán)境時要重點(diǎn)關(guān)注光照強(qiáng)度與二氧化碳濃度。

表4 高低風(fēng)速時各環(huán)境因素間相關(guān)系數(shù)Tab.4 Correlation coefficient between high and low wind speed and various environmental factors

2.2 環(huán)境因素非日均值與產(chǎn)奶量之間的關(guān)系

不同百分位值二氧化碳濃度、光照強(qiáng)度與產(chǎn)奶量間的相關(guān)系數(shù)見表5。可以看出二氧化碳濃度與產(chǎn)奶量的相關(guān)系數(shù)在所有百分位上均呈負(fù)相關(guān)，并且百分位數(shù)越小，兩者的負(fù)相關(guān)程度越高。相反，光照強(qiáng)度與產(chǎn)奶量的相關(guān)系數(shù)在所有百分位上均呈正相關(guān)，且隨著百分?jǐn)?shù)增大，兩者的正相關(guān)程度有增大的趨勢。光照強(qiáng)度第90百分位數(shù)、二氧化碳濃度第10百分位數(shù)與產(chǎn)奶量關(guān)系如圖5所示。可以看出，高產(chǎn)奶量樣本集中在光照強(qiáng)度第90百分位數(shù)大于800 lx且二氧化碳濃度第10百分位數(shù)低于6×10-4的區(qū)域。由于采樣間隔為10 min，大于等于第90百分位數(shù)的數(shù)據(jù)有15個，小于等于第10百分位數(shù)的數(shù)據(jù)有14個。因此為了提升產(chǎn)奶量，建議奶牛舍每天至少保證2.5 h不低于800 lx的光照時長。同時，建議控制二氧化碳濃度高于6×10-4的時長不超過2.33 h。

表5 光照強(qiáng)度與二氧化碳濃度非日均值與產(chǎn)奶量相關(guān)系數(shù)Tab.5 Correlation coefficient between non-daily mean values of illumination intensity and carbon dioxide concentration and milk yield

2.3 隨機(jī)森林回歸模型構(gòu)建

2.3.1自變量篩選降維

自變量篩選降維時要至少保留各個環(huán)境因素的一個自變量，剔除重要性較小的自變量后所訓(xùn)練的模型見表6。模型編號1～7分別表示因變量未作任何變換、平方根變換、平方變換、立方變換、自然對數(shù)變換、指數(shù)變換、倒數(shù)變換所訓(xùn)練的模型；自變量欄中的大寫字母F、W、S、C、A、G分別表示風(fēng)速、溫度、相對濕度、二氧化碳濃度、氨氣濃度、光照強(qiáng)度；ave、max、min分別代表日均值、每日12個最大值平均、每日12個最小值平均；數(shù)字代表其所在百分位數(shù)。可以看出，除了X和1/X兩種形式外，其余形式的自變量均在自變量篩選過程中被剔除。這可能是由于自變量在采用其他形式變換后線性關(guān)系不再明顯，不利于決策樹在分裂時找到合適的分裂點(diǎn)。同時可以看出在多個模型中光照強(qiáng)度、二氧化碳濃度、風(fēng)速、溫度出現(xiàn)次數(shù)較多，表明這4種環(huán)境因素是反映產(chǎn)奶量的重要自變量，其中光照尤為重要。

表6 篩選自變量后模型性能Tab.6 Model performance after filtering independent variables

除了第6個模型性能明顯劣于其他模型外，其余6個模型性能相近。說明除了對Y進(jìn)行自然指數(shù)變換，Y的其他變換形式對隨機(jī)森林回歸模型的構(gòu)建并無明顯影響。這可能是由于自然指數(shù)變換會導(dǎo)致Y的變換值急劇增大從而導(dǎo)致模型精度降低。因此，本文選用未對因變量進(jìn)行變換的模型進(jìn)行后續(xù)分析。

2.3.2確定最佳訓(xùn)練集

經(jīng)過自變量篩選降維后，剩余預(yù)測性能較好的6個模型中共有25個不重復(fù)的自變量，將其作為新模型的自變量，以未變換的Y作為因變量進(jìn)行性能評估。根據(jù)逐步回歸的方法對剩余25個自變量進(jìn)行逐步回歸，回歸結(jié)果見表7。可以看出，在剩余17個自變量時，模型性能達(dá)到最高值，隨后性能逐步降低，在剩余6個自變量時，模型性能低于25個自變量建立的預(yù)測模型，篩選到最后剩余G80變量。考慮到在保證一定模型精度的前提下，盡可能刪減自變量個數(shù)，因此選取7個自變量(G80、Fave、1/Cmin、Gmin、Smin、Gmax、1/A50)建立最佳模型。

表7 逐步回歸確定最佳模型參數(shù)Tab.7 Step by step regression to determine the best model process

2.3.3優(yōu)選模型預(yù)測效果

將完整樣本隨機(jī)分為4個小樣本集，選取1個作為驗(yàn)證集，其余作為訓(xùn)練集，進(jìn)行交叉驗(yàn)證評定模型的泛化能力。圖6為交叉驗(yàn)證過程中較好模型與較差模型的預(yù)測性能。全部模型對驗(yàn)證集預(yù)測值的平均決定系數(shù)R2為0.731 6，平均均方根誤差E為1.037 0 kg，說明模型預(yù)測值與實(shí)際值之間趨勢基本一致。兩種模型的性能差異主要來自個別樣本的較大誤差，尤其是對某些產(chǎn)奶量區(qū)間的預(yù)測精度較低，這可能是由于隨機(jī)抽取訓(xùn)練集導(dǎo)致某些產(chǎn)奶量區(qū)間的樣本量過少，模型對該區(qū)間預(yù)測能力降低所致。兩個模型雖然性能上有所差異，不過整體趨勢與實(shí)際值是一致的。

3 結(jié)論

(1)牛舍運(yùn)行過程中氣體分布會受到風(fēng)機(jī)干擾，在高風(fēng)速時舍內(nèi)部分環(huán)境因素彼此顯著相關(guān)，但與產(chǎn)奶量的相關(guān)性沒有達(dá)到顯著性水平，此時評價(jià)舍內(nèi)環(huán)境時主要以溫度、相對濕度、風(fēng)速為主。在低風(fēng)速時，光照強(qiáng)度日均值250 lx與二氧化碳濃度8×10-4可以明顯地劃分高低產(chǎn)奶量。

(2)二氧化碳濃度第10百分位數(shù)與光照強(qiáng)度第90百分位數(shù)與產(chǎn)奶量相關(guān)性最大，為了提升產(chǎn)奶量，建議奶牛舍每天至少保證不低于800 lx的光照時長2.5 h。同時，建議控制二氧化碳濃度高于6×10-4的時長不超過2.33 h。

(3)將產(chǎn)奶量作為預(yù)測指標(biāo)，利用隨機(jī)森林方法建立了基于環(huán)境因素的產(chǎn)奶量回歸模型。模型泛化能力檢驗(yàn)結(jié)果表明其具有良好的預(yù)測能力。