荷斯坦奶牛泌乳曲线拟合分析

栗敏杰，任小丽，闫 磊，闫跃飞，白雪利，李静茹，李学钦，邓好好，张 震*，高腾云*

(1.河南农业大学 牧医工程学院，河南 郑州 450002；2.河南省奶牛生产性能测定中心，河南 郑州 450000；3.偃师市高龙镇隆鑫奶牛养殖场，河南 洛阳 471931；4洛阳生生乳业有限公司，河南 洛阳 471199)

泌乳曲线是奶牛在一个完整泌乳期内产奶量随时间变化的曲线。泌乳曲线能综合反映奶牛的泌乳能力、牧场营养供给水平、奶牛的利用年限及牧场的总体经济效益，因此针对泌乳曲线的研究在奶牛泌乳性能评定、饲养管理指导及品种选育方面都有重要意义[1-4]。从1967年Wood模型提出至今，泌乳曲线拟合模型已达20余种[5]，目前应用较多的是Wood模型[6]、逆多项式模型[7]、AS模型[8]及Wilmink模型等。已经有较多关于奶牛泌乳曲线的相关研究，如：Pal-liser等[9]采用动态的MOLLY模型，构建了新西兰奶牛的泌乳曲线；亓建刚等[10]利用Wood模型对不同产犊季节下1～3胎牛进行泌乳曲线拟合，分析不同产犊季节对各胎次泌乳曲线各参数的影响。Olori等[11]选择AS、IQP、ML、WIL、Wood等5种模型对英国488头初产荷斯坦奶牛进行了泌乳曲线拟合研究，其中AS模型拟合效果最好，Wood模型效果最差；付雪峰等[12]利用Wood、IQP、多项式回归、AS和WIL等5种模型对2843头新疆褐牛泌乳曲线进行拟合，AS的模型拟合效果最好；另王瑞军等[13]、马光辉等[14]对国内其他牛只品种也做过wood模型的拟合研究，拟合效果良好。不同生态类型区域奶牛的泌乳曲线存在差别，适宜的拟合模型及参数也不一样。中原地区是我国奶牛优势产区之一，并且是很有发展潜力的产区，目前还没有针对我国中原地区规模化牧场奶牛泌乳曲线的研究。本研究对中原地区两个规模化牧场测定间隔为5～30 d的日产奶量记录严格质控后拟合泌乳曲线，比较不同胎次及牧场下Wood模型和AS模型的拟合效果，以选择最优拟合模型，为305天产奶量计算和奶牛不同泌乳阶段营养管理提供依据。

1 材料与方法

1.1 数据来源及质控

本研究选择了2个牧场2006-2018年度1～6胎次的40003条日产奶量记录，按下列质控标准进行筛选：①胎次共4个水平，将4～6胎归类为≥4胎；②删除DIM<5 d或>305 d的记录；③一个完整泌乳期内最小泌乳天数<30 d、记录10次(包含)以上，记录间隔<30 d；④保留日产奶量5～80 kg·d-1的记录。得到374头荷斯坦奶牛的12829条日产奶量记录作为试验数据资料。

1.2 选择的奶牛泌乳曲线模型

本研究选择两种模型进行泌乳曲线拟合：

(1)Wood 不完全伽马函数模型[6]：Yt=atbe-ct。式中：t为泌乳天数(d)；Yt为该泌乳天数时的日产奶量；a、b、c为模型参数；a是泌乳潜力；b是产奶量在泌乳高峰后下降的速度；c是产奶量在泌乳高峰前上升的速度；e为自然常数。

在进行R-NLS非线性拟合时，各参数初始值采用亓建刚等[10]计算所得结果。拟合所得a、b、c等3个参数后，利用如下公式换算MPY、PY、Per:MPY=b/c；PY=a(b/c)be-b；Per=-(b+1)1n(c)。式中:MPY是达到泌乳高峰所需要的天数，即泌乳高峰日(d);PY为高峰产奶量(kg);e为自然常数;Per为泌乳持续力。

(2)AS模型(Ali and Schaeffer,1987)[8,11]：Yt=a+bt+ct2+dlogt+e(logt)2

式中：t为泌乳天数(d)；Yt为该泌乳天数时的日产奶量；a、b、c、d、e为模型参数；logt是以e为底；在进行R-NLS非线性拟合时，各参数初始值采用毛永江等[15]计算所得结果。

1.3 统计模型拟合优良性评价标准-AIC

1.4 统计分析

使用Excel 2013对数据进行初步筛选，利用SAS 9.4建立数据质控标准对数据进行质控，使用R v3.5.1软件中NLS模块进行非线性拟合，该模块运用Gauss-Newton迭代算法计算各模型参数。AIC由R v3.5.1软件中的logLik过程得到的极大似然值计算得到。

2 结果与分析

2.1 不同胎次、牧场日产奶量描述性统计量

由表1可见，头胎个体日产奶量均低于经产牛，产奶量在1号、2号牧场的不同胎次下变化趋势不同；1号牧场2胎次日产奶量最高，2号牧场3胎次产奶量最高；1号牧场1～4胎次的日产奶量平均值比2号牧场高5.15 kg。

表1 不同胎次、牧场日产奶量描述性统计量Table 1 改为Descriptive statistics of daily milk yield for different parities and farms

2.2 不同胎次、牧场泌乳曲线拟合

如表2所示，不同胎次、牧场的泌乳曲线二级参数表现不同。头胎个体泌乳高峰日和高峰产奶量均晚于或低于经产牛，头胎牛泌乳持续力高于经产牛。胎次相同时，1号牧场泌乳高峰日和高峰产奶量均高于2号牧场。胎次为1～3胎时，1号牧场泌乳持续力高于2号牧场。Wood模型拟合不同胎次、牧场的泌乳曲线如图1所示。

2.3 不同模型泌乳曲线拟合

由表3可见，随着胎次的增加，Wood模型和AS模型的AIC逐渐降低，即胎次越高，各模型的拟合效果越好。AS模型对于1胎次和3胎次拟合效果优于Wood模型，Wood模型对于2胎次拟合效果优于AS模型。4胎次时，AS模型适合于1号牧场，Wood模型适合2号牧场(图2)。

3 讨 论

相关研究结果表明，理想泌乳高峰日约在产后40～70 d[12,28]。除第1胎外，本研究中，2～4胎次的泌乳高峰日均在产后40～70 d内。1胎泌乳高峰日最晚，且高峰产奶量最低，与类似研究结果一致[15,29-31]，可能原因是个体在初次产奶时仍处于生长发育阶段，没有完全发挥泌乳潜力。

表2 不同胎次和牧场Wood模型群体泌乳曲线拟合参数Table 2 Lactation curve fitting parameters with Wood model for different parities and farms

高峰产奶量与胎次总产奶量呈强正相关。高峰产奶量每上升1 kg，胎次产奶量会上升200～300 kg[32]。如图1所示，1号牧场的2胎高峰产奶量最高，2胎以后逐渐降低，毛永江等[15]研究结果与之一致。实际生产中可通过分阶段分群饲养、添加饲料中精料、增加牛群饲喂次数及提高干物质采食量等措施来补充奶牛摄入能量进而提高高峰产奶量[28]，2号牧场各胎次高峰产奶量变化规律与1号牧场不同，牧场间生产水平不同、奶牛日粮中精粗料比例不同可能是两个牧场高峰产奶量变化规律不同的主要原因。

本研究中，泌乳高峰日越早、泌乳持续力越低，这与程郁昕等[33]的研究结果一致。泌乳持续力随着胎次的增加而逐渐减弱，初产牛的泌乳持续力最好，与贾先波[34]的研究相似。较强泌乳持续力得益于牧场科学的饲养管理，在达到泌乳高峰时，供给奶牛充足的营养，预防产后能量负平衡，从而增强牛群泌乳持续力、提高泌乳期产奶量。

在刘丽元等[35]对新疆昌吉地区荷斯坦奶牛泌乳曲线拟合的研究中指出，Wood模型相较于AS模型更适合于2胎泌乳曲线拟合，与本研究结果一致。随着胎次的增加，Wood模型和AS模型的AIC均逐渐降低，与亓建刚等[10]得出的北京地区奶牛1～3胎Wood模型拟合泌乳曲线的拟合度逐渐上升的结果相一致；但毛永江等[15]研究得出，南方地区奶牛1～3胎的wood模型和AS模型的拟合度逐渐降低，拟合效果逐渐减弱，与本研究结果相反，造成此差异的原因极有可能是不同生态类型区域的牧场所处的自然环境不同造成奶牛泌乳曲线变化趋势差别较大。

4 结 论

本研究通过拟合不同胎次的泌乳曲线得出，荷斯坦初产牛的泌乳高峰日晚于经产牛，高峰产奶量低于经产牛，但泌乳持续力高于经产牛；Wood模型和AS模型对泌乳曲线的拟合度随着胎次的增加均逐渐上升，且AS模型为1、3胎最佳拟合模型，Wood模型为2胎最佳拟合模型。在实际生产中，根据牧场历史产奶量数据拟合泌乳曲线，可以更准确、更精细化的指导牧场生产管理。