基于随机森林的寒区奶牛舍环境因素与产奶量关系研究

丁 涛 蒋雷生,2 施正香,3 赵 洋,2 马 慧

(1.中国农业大学水利与土木工程学院， 北京 100083； 2.北京市供水管网与安全节能中心， 北京 100083；3.农业农村部设施农业工程重点实验室， 北京 100083； 4.北京首农畜牧发展有限公司， 北京 100083)

0 引言

奶牛舍环境对奶牛的生产性能具有决定性影响，当牛舍小气候变化幅度超出奶牛的适应范围时，奶牛的生产力和健康便会受到负面影响[1]。目前，评价牛舍环境的指标大部分是基于热环境的热应激评价方法[2]，如THI[3]、BGHI[4]、THVI[5]、HLI[6]以及不同的体感温度表达式[7-9]，这些指数关注了温度、相对湿度、风速和热辐射，却忽略了光照强度及空气质量的影响，不适用于非热应激季节的奶牛舍环境评估。由于季节特殊性，冬季舍内风速低、有害气体浓度高，特别是二氧化碳浓度和氨气浓度对奶牛的健康和生产有重要影响。研究表明[10-12]，光照使产奶量最大化的明暗交替最佳节点是16 h光照，8 h黑暗，合理的光照能提高产奶量8%～10%、促进产后发情和受孕。然而针对寒区的气候特点，光照对于奶牛生产的影响没有更加深入地研究。

随机森林算法是数据挖掘分类算法的重要组成部分，在生物信息学、生态学、医学等方面都被广泛应用[13-14]。针对畜禽舍环境的调控复杂性，国内外已经有采用模糊算法、神经网络等手段实现舍内环境的智能调节[15-17]。有学者利用奶牛数据和机器学习算法预测奶牛泌乳早期的代谢状态，其中随机森林算法和支持向量机的表现最优[18]。本文通过对奶牛舍环境参数以及产奶量的连续监测，从不同角度分析其中的变化规律，建立基于环境因素的随机森林回归预测模型，并验证其在寒区奶牛舍的适用性。

1 材料与方法

1.1 数据集

试验于2018年9月至2019年3月在吉林省白城市首都农业集团白城牛场进行。牛舍东西长度约250 m，南北跨度约27 m，屋顶为彩钢夹芯板双坡屋顶，奶牛采用全舍饲散栏饲养模式，舍内奶牛约580头，每日挤奶3班。牛舍四周窗户除冬季外均开启，屋顶小窗常年开启，南面屋顶上每隔12 m均匀布置有采光板增加舍内采光，牛舍内无喷淋降温设施。舍内沿东西方向共有4排工字钢，东西方向相邻工字钢间隔6 m。在中间两排工字钢上间隔12 m安装36寸风机，在躺卧区域上方间隔15 m安装54寸风机，舍内共安装有40台36寸风机、32台54寸风机。风机在夏季全部24 h运行，秋季白天部分风机在高温时段对称开启，晚上关闭。风机运行个数随温度降低而减少直至全部停止运行，试验期间舍内大部分时间都处于风速较低状态。

试验牛舍由于跨度不大、采光板的设置及非热应激时期太阳高度角较小，前期通过手持testo540型光照度测试仪测试奶牛活动区域，结果表明南北光照强度差异不大。牛舍光照为自然光照，舍内在高度4 m处间隔6 m安装15 W节能灯，共安装80个节能灯，如图1所示，夜间开启节能灯便于人员工作。为了避免距离风机过近数据不准确，因此将传感器安装在相邻两个风机中间的工字钢上，距离地面1.8 m。舍内各项设施布置均匀，舍外附近无其余建筑物，因此布置4个环境数据采集点可以较好地反映舍内环境水平。测点分别距离奶牛舍东门、西门45 m及90 m，舍外布置1个采集点E，位于奶牛舍舍外东南角，如图2所示。每个采集点可以采集温度、相对湿度、风速、二氧化碳浓度、氨气浓度、光照强度。传感器为武汉中科能慧公司生产，采样间隔为10 min，具体参数见表1。选择泌乳性能相近的30头奶牛每日产奶量进行分析。

1.2 数据处理

采集的数据往往存在部分缺失、不正确、含有噪声或其他类型的不一致现象[19]。采用删除法与插补法处理数据缺失的情况，同时利用Matlab进行小波去噪以减少噪声对数据的影响。考虑到环境因素的变化达到一定的阈值才会对畜禽的生产能力产生影响[20]，为了确定当日环境数据的代表值，本文将数据从小到大进行排列取不同百分位值及均值，研究其与产奶量的相关关系。

表1 试验传感器参数Tab.1 Test sensor parameter list

1.3 机器学习

随机森林算法通过自助法在原数据集中有放回地随机抽取n个样本集，每次抽取2/3样本容量的数据作为袋内数据，建立m棵决策树构建随机森林，计算m棵决策树回归结果的平均值进行预测[21]。该方法由于引入“双随机”的思想与组合投票取均值原则，因此对噪声和异常值不敏感、不易陷入过拟合[22]。

1.3.1数据预处理

构建机器学习模型时，使用样本的平均误差最小化作为模型的学习准则，则少数类样本的误差在总体误差中的比重会偏小，使得模型对少数类样本的预测性能下降。本文采用合成少数类过采样技术方法，利用随机线性插值法，在两个距离较近的少数样本间创造新的少数类样本来平衡样本。由于各环境因子的取值方式对产奶量的影响程度尚不明确，因此将各环境因子的每日平均值、每日12个最大值平均、每日12个最小值平均以及第10、20、35、50、65、80、90百分位数组成基础自变量集X，将产奶量作为基础因变量集Y。对自变量和因变量分别采取平方根变换、平方变换、立方变换、自然对数变换、指数变换、倒数变换。经过变换后形成7个训练集，每个训练集包含420个自变量。

1.3.2模型评价指标

模型训练的目标是实现对训练样本的最佳拟合，是基于误差最小化原则，确定模型最优参数与结构的过程。选用决定系数(R2)和均方根误差(E)作为模型评价指标，同时将两者合并为指标P以便对比模型性能，其计算公式为

(1)

E0——剔除某自变量前模型训练的均方根误差

Ei——剔除某自变量后模型训练的均方根误差

1.3.3模型参数优选

随机森林回归模型的构建需要确定决策树数量，对420个自变量筛选降维，确定最佳训练集。根据各模型的袋外数据均方误差估计随决策树数量增长的变化曲线确定决策树数量为120棵[23]。基于袋外数据估计的特征重要性，逐步剔除重要性较小的特征，并迭代训练随机森林模型，直到自变量的数量降低到目标个数。最后利用R2和E对优选后的模型进行评价。

2 结果与讨论

2.1 环境因素日均值与产奶量之间的关系

试验期间每月的数据变化情况见表2,目前研究表明奶牛舒适温度区间为5～25℃，可以看出在9—11月舍内温度处于奶牛舒适的温度范围内时，相对湿度、风速、光照强度、二氧化碳浓度、氨气浓度、产奶量均存在显著性差异。9月和次年1月产奶量近似，期间除二氧化碳浓度外，其他因素均显著变化。2月和3月产奶量近似，除二氧化碳浓度与风速外，其他因素均存在显著性差异。表明在部分时段里二氧化碳浓度不存在显著性差异时，其他环境因素的变化对产奶量没有产生明显地促进或抑制作用。

为了减少自变量从而准确分析其他因素与产奶量的规律，引入温湿指数(Temperature humidity index,THI)，温湿指数计算式为

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

(2)

式中T——干球温度H——相对湿度

温湿指数与产奶量的关系如图3所示。由图3可知，THI与产奶量间并没有明显的变化规律。分析其原因，在9、10月舍内风机随温度变化决定开启数量，风机运行影响了舍内气体分布。当温度较低时，牛舍内风速处于较低水平，有害气体浓度增加，光照强度也处于较高的水平。不同环境因素在不同时间段对奶牛的影响不同，因此在同一THI水平及不同THI水平时，产奶量动态变化。为了确定不同区间内影响产奶量的重要因素，根据产奶量及THI划分为不同阶段进行分析，结果见表3。

表2 试验期间各月环境数据Tab.2 Monthly environmental data during test period

由表3可知，在THI处于较低和中等水平时，高产奶量和中产奶量间光照强度、二氧化碳浓度存在显著差异。中产奶量和低产奶量间光照强度不存在明显差异，二氧化碳浓度较高时产奶量低。THI处于较高水平时，高产奶量和低产奶量样本间光照强度存在显著性差异，二氧化碳浓度不存在显著性差异。表明当二氧化碳浓度较低时产奶量较高，光照对产奶量有进一步的促进作用。

表3 不同THI与产奶量阶段时光照强度与二氧化碳浓度关系Tab.3 Relationship between illumination intensity and carbon dioxide concentration at different THI and milk yield stage

光照能抑制褪黑激素的分泌，其分泌的交替变化能引发一系列的激素反应，从而提升产奶量[24]。PETERS等[25]研究表明接受16 h光照比9～12 h光照奶牛产奶量高6.7%。ESPINOZA等[26]研究表明保持16 h光照相比8 h光照奶牛产奶量高2.2 kg/d。REKSEN等[27]研究结果表明在高纬度地区，夜晚接受光照比未接受光照的奶牛产奶量高0.5 kg/头。MUTHURAMALINGAM等[28]研究了不同光照强度对褪黑激素的影响水平，当夜晚光照强度为50 lx时，褪黑激素分泌会受到抑制。目前普遍认为光照强度200 lx有助于生产，各类研究明确了光照对产奶量的促进作用。由于各环境因素对奶牛的影响是动态的，本次试验主要分析正常生产状况下多个环境因素与产奶量之间关系，确定不同阶段各环境因素对产奶量的影响程度。夜间光照只用于照明，光照强度包含在每日环境因素日均值中，不影响本文分析结果。

东北地区昼夜温差大，冬季温度较低，牛舍的封闭措施导致二氧化碳浓度过高，继而引起氧气浓度降低，造成奶牛缺氧[29]。梅玮等[30]研究表明牛舍中二氧化碳和氨气浓度过高将不利于奶牛健康和生产性能。目前关于牛舍内二氧化碳浓度的研究多是关于气体分布规律方面的，包括不同牛舍类型、不同季节等方面，并以二氧化碳浓度1 500 mg/m3为标准进行舍内环境评价[31-33]。

考虑到随着温度的改变，人工管理会影响舍内环境因素的分布规律。将观测样本划分为高风速与低风速时各环境因素的相关系数见表4。可以看出，当温度较高时，风机运行的个数与舍内温度有关，因此风速与温度、相对湿度、二氧化碳浓度均存在明显的相关性，各因素与产奶量间存在相关性，但不显著。温度与产奶量相关系数较高，此时评价舍内环境时应以温度为主要因素。当温度较低时，舍内风速一直处于较为稳定的低水平，温湿度呈明显的线性规律。二氧化碳浓度与产奶量相关系数为-0.532，光照强度与产奶量相关系数为0.720，二者均达到显著水平。温度、光照强度、二氧化碳浓度与产奶量的颜色映射图如图4所示，可以看出光照强度日均值250 lx与二氧化碳浓度日均值8×10-4可以明显地划分高低产奶量。因此建议在低温时评价奶牛舍环境时要重点关注光照强度与二氧化碳浓度。

表4 高低风速时各环境因素间相关系数Tab.4 Correlation coefficient between high and low wind speed and various environmental factors

2.2 环境因素非日均值与产奶量之间的关系

不同百分位值二氧化碳浓度、光照强度与产奶量间的相关系数见表5。可以看出二氧化碳浓度与产奶量的相关系数在所有百分位上均呈负相关，并且百分位数越小，两者的负相关程度越高。相反，光照强度与产奶量的相关系数在所有百分位上均呈正相关，且随着百分数增大，两者的正相关程度有增大的趋势。光照强度第90百分位数、二氧化碳浓度第10百分位数与产奶量关系如图5所示。可以看出，高产奶量样本集中在光照强度第90百分位数大于800 lx且二氧化碳浓度第10百分位数低于6×10-4的区域。由于采样间隔为10 min，大于等于第90百分位数的数据有15个，小于等于第10百分位数的数据有14个。因此为了提升产奶量，建议奶牛舍每天至少保证2.5 h不低于800 lx的光照时长。同时，建议控制二氧化碳浓度高于6×10-4的时长不超过2.33 h。

表5 光照强度与二氧化碳浓度非日均值与产奶量相关系数Tab.5 Correlation coefficient between non-daily mean values of illumination intensity and carbon dioxide concentration and milk yield

2.3 随机森林回归模型构建

2.3.1自变量筛选降维

自变量筛选降维时要至少保留各个环境因素的一个自变量，剔除重要性较小的自变量后所训练的模型见表6。模型编号1～7分别表示因变量未作任何变换、平方根变换、平方变换、立方变换、自然对数变换、指数变换、倒数变换所训练的模型；自变量栏中的大写字母F、W、S、C、A、G分别表示风速、温度、相对湿度、二氧化碳浓度、氨气浓度、光照强度；ave、max、min分别代表日均值、每日12个最大值平均、每日12个最小值平均；数字代表其所在百分位数。可以看出，除了X和1/X两种形式外，其余形式的自变量均在自变量筛选过程中被剔除。这可能是由于自变量在采用其他形式变换后线性关系不再明显，不利于决策树在分裂时找到合适的分裂点。同时可以看出在多个模型中光照强度、二氧化碳浓度、风速、温度出现次数较多，表明这4种环境因素是反映产奶量的重要自变量，其中光照尤为重要。

表6 筛选自变量后模型性能Tab.6 Model performance after filtering independent variables

除了第6个模型性能明显劣于其他模型外，其余6个模型性能相近。说明除了对Y进行自然指数变换，Y的其他变换形式对随机森林回归模型的构建并无明显影响。这可能是由于自然指数变换会导致Y的变换值急剧增大从而导致模型精度降低。因此，本文选用未对因变量进行变换的模型进行后续分析。

2.3.2确定最佳训练集

经过自变量筛选降维后，剩余预测性能较好的6个模型中共有25个不重复的自变量，将其作为新模型的自变量，以未变换的Y作为因变量进行性能评估。根据逐步回归的方法对剩余25个自变量进行逐步回归，回归结果见表7。可以看出，在剩余17个自变量时，模型性能达到最高值，随后性能逐步降低，在剩余6个自变量时，模型性能低于25个自变量建立的预测模型，筛选到最后剩余G80变量。考虑到在保证一定模型精度的前提下，尽可能删减自变量个数，因此选取7个自变量(G80、Fave、1/Cmin、Gmin、Smin、Gmax、1/A50)建立最佳模型。

表7 逐步回归确定最佳模型参数Tab.7 Step by step regression to determine the best model process

2.3.3优选模型预测效果

将完整样本随机分为4个小样本集，选取1个作为验证集，其余作为训练集，进行交叉验证评定模型的泛化能力。图6为交叉验证过程中较好模型与较差模型的预测性能。全部模型对验证集预测值的平均决定系数R2为0.731 6，平均均方根误差E为1.037 0 kg，说明模型预测值与实际值之间趋势基本一致。两种模型的性能差异主要来自个别样本的较大误差，尤其是对某些产奶量区间的预测精度较低，这可能是由于随机抽取训练集导致某些产奶量区间的样本量过少，模型对该区间预测能力降低所致。两个模型虽然性能上有所差异，不过整体趋势与实际值是一致的。

3 结论

(1)牛舍运行过程中气体分布会受到风机干扰，在高风速时舍内部分环境因素彼此显著相关，但与产奶量的相关性没有达到显著性水平，此时评价舍内环境时主要以温度、相对湿度、风速为主。在低风速时，光照强度日均值250 lx与二氧化碳浓度8×10-4可以明显地划分高低产奶量。

(2)二氧化碳浓度第10百分位数与光照强度第90百分位数与产奶量相关性最大，为了提升产奶量，建议奶牛舍每天至少保证不低于800 lx的光照时长2.5 h。同时，建议控制二氧化碳浓度高于6×10-4的时长不超过2.33 h。

(3)将产奶量作为预测指标，利用随机森林方法建立了基于环境因素的产奶量回归模型。模型泛化能力检验结果表明其具有良好的预测能力。