基于近红外漫反射光谱的规模化奶牛场粪水氮磷定量分析及模型构建

时间：2024-05-24

赵润，牟美睿，王鹏，孙迪，刘海学，张克强*，杨仁杰*

（1.农业农村部环境保护科研监测所，天津 300191；2.天津农学院农业分析测试中心，天津 300384；3.天津农学院工程技术学院，天津 300384）

近年来奶牛场规模化程度的提升造成环境压力凸显，特别是夏季挤奶车间和牛舍降温用水量加大、雨污混流等造成粪水量陡增而降低了收运和处理效率，大量粪水和沼液的出路成为现阶段阻碍奶牛场正常运转乃至奶业绿色转型发展的首要难题。国内外实践经验表明，还田农用是奶牛场粪水和沼液的必然出路，然而缺少粪水中氮磷等成分的快速检测方法是造成还田难的技术瓶颈。有研究表明，时空差异、圈舍类型、养殖规模、清粪方式、粪污处理工艺等因素均会不同程度地影响粪污中氮磷含量变化[1-2]，但通过常规采样监测来获取粪水中氮磷含量动态变化的做法存在时效性低、成本高、适用范围有限及准确度低等问题，因此有必要建立满足规模化奶牛场现场影响因素条件下粪水中氮磷含量的快速检测方法。近期我国密集出台的政策性文件中也明确指出要健全畜禽粪污还田利用和检测方法标准体系，制定畜禽粪水、沼液还田利用技术规范，因此建立适合现场的快速检测方法对指导奶牛场粪水安全和科学还田及防治环境污染意义重大。

近红外光谱作为一种快速、便捷、可实现在线检测的方法已被应用于畜禽粪污中组分的检测[3-11]。王晓燕等[12]采用近红外光谱对鸡粪中的氮磷钾进行定量分析，对该方法的可行性进行了验证。Chen等[13]、严旭等[14]、黄光群等[15]和渠清博等[16]综述了近红外光谱技术在预测动物粪污中干物质、氮和磷等含量方面的应用，指出近红外光谱技术作为监测畜禽粪污化学成分的有效手段具有很大的潜力。Jancewicz等[17]基于近红外光谱对65种不同日粮的牛粪样品进行分析，结果表明近红外光谱能较好地预测奶牛粪便中大部分营养物质的化学组成和消化率。同时，近红外光谱技术还被应用于猪粪便中氮的全肠道表观消化率研究[18-19]，结果表明近红外光谱技术能准确预测粪便中的化学成分。但目前发表的论文大多围绕单一环节、单一尺度畜禽粪便样品的组分和含量进行近红外光谱定量化预测的可行性分析，未见面向规模化奶牛场粪水处理全过程环节来建立符合现场复杂、复合影响因素条件下近红外光谱快速检测方法及其全尺度综合模型的报道。这也是本研究中的主要技术难点。

本文选取天津市23家典型种养结合模式的规模化奶牛场，开展动态复合影响因素下全局、全过程近红外漫反射光谱快速检测方法与模型构建，为实现满足现场条件下的快速准确定量化预测，破解规模化奶牛场粪水还田难和推进奶业绿色发展提供技术支撑。

1 材料与方法

1.1 研究对象概况

1.1.1 典型奶牛场筛选

基于对天津市85家正常生产的规模化奶牛场情况摸底，分别从全市5个奶牛场集中分布地区（武清区、滨海新区、静海区、宁河区、北辰区）中选取23家采用种养结合模式的典型规模化奶牛场，其共性特点包括：饲养品种均为中国荷斯坦奶牛，全群自繁自养，按照不同牛群生长阶段划分区域，80%的奶牛场具有10年以上的经营历史，常年稳定运转。泌乳牛多为卧床式养殖，卧床垫料为固液筛分后的干牛粪，育成牛为散栏式养殖，泌乳牛舍中产生的粪水全部进入处理系统。不同之处在于：总存栏规模在400～5000头不等；清粪方式和粪水处理工艺组合各不相同，如清粪方式包括干清粪、干清粪-水冲粪、水冲粪等，清粪设备包括刮粪板、铲车、吸粪车等，粪水收储设施包括集粪沟、集粪渠、收集池等，粪水处理方式包括固液筛分、厌氧消化、好氧发酵等，粪水处理设施包括集污池、分离池、调节池、沉淀池、氧化塘等，粪水处理设备包括USR、UASB、CSTR、PFR等厌氧反应器；部分奶牛场育成牛的粪便进入粪水处理系统，部分不进；部分奶牛场挤奶车间（奶厅）的废水进入处理系统，部分不进。选取的23家典型规模化奶牛场综合了上述情况特点，涵盖了可能对粪水氮磷含量测定结果造成影响的现场要素。

1.1.2 样品采集

选取滨海新区某代表性规模化奶牛场作为纵向采样对象，以上述5个地区23家规模化奶牛场作为横向采样对象，参照《规模化奶牛场环境监测技术规程》（DB 12/T 655—2016）[20]，采样布点从牛舍内和奶厅粪水收储设施开始，到粪水还田前的末级贮存设施为止，遍及收集池、集粪沟、集污池、分离池、调节池、沉淀池、贮存池和氧化塘等所有粪水流经设施。纵向和横向采样对象的粪水处理工艺路线和采样位点分布分别如图1和图2所示。收集池为奶厅废水收储设施，集粪沟为牛舍粪水汇集点，集污池为奶牛场所有粪水交汇点，分离池为固液筛分后的粪水暂存设施，调节池为进入沼气工程前匀浆调质的设施，沉淀池为厌氧发酵后分离沼液和沼渣的设施，贮存池和氧化塘为粪水还田前的贮存设施。

参照《畜禽养殖污水采样技术规范》（GB/T 27522—2011）[21]和《规模化奶牛场环境监测技术规程》（DB 12/T 655—2016）[20]，于2018年12月下旬采集上述23家规模化奶牛场的粪水样品，用自制不锈钢材质提桶、水舀等工具在各设施采样点垂直液面下10～20 cm处随机采集3个位点的水样，并在大水桶中用水舀混合搅匀后取500 mL置于集水瓶中，送回实验室即时检测。连续6 d采集滨海新区某典型奶牛场的粪水样品共计57个，采集23家奶牛场的粪水样品共计138个。

1.2 仪器与分析条件

1.2.1 近红外漫反射光谱的采集

图1 滨海新区某奶牛场粪水处理工艺流程和采样位点分布图Figure 1 Flow chart of slurry treatment route within sampling distribution in one dairy farm of Binhai District

图2 23家规模化奶牛场粪水处理总工艺流程和采样位点分布图Figure 2 Flow chart of slurry treatment route within sampling distribution in 23 large-scale dairy farms

实验采用美国PerkinElmer公司的傅里叶变换近红外光谱仪，InGaAs检测器，扫描范围为4000～12 000 cm-1。将集水瓶中待测粪水样品充分摇匀后，用3 mL一次性胶头滴管取集水瓶中部2～3 mL样品装入样品杯中，并放置在积分球旋转样品台上（仪器自带积分球附件），以积分球内置参比为背景，分别采集每一个样品的近红外漫反射光谱。光谱扫描参数为：分辨率8 cm-1，扫描间隔2 cm-1，扫描次数64次。

1.2.2 总氮和总磷的测定

按照《水质凯氏氮的测定》（GB 11891—1989）[22]中规定的方法测定粪水中总氮（TN）的含量，仪器选用全自动凯氏定氮仪（Foss kjeltec 8400型，丹麦）。按照《水质总磷的测定钼酸铵分光光度法》（GB 11893—1989）[23]中规定的方法测定粪水中总磷（TP）的含量，仪器选用可见光分光光度计（722E型，中国）。对23家规模化奶牛场采集的138个粪水样品中TN和TP的含量进行测定，统计结果如表1所示，表明受养殖规模、清粪方式、粪水处理工艺等多种复合因素影响，不同奶牛场粪水处理全过程各环节粪水中TN和TP含量变化较大。

表1 23家规模化奶牛场粪水样品TN和TP含量统计结果Table 1 Statistical results of TN and TP in slurry samples of 23 large-scale dairy farms

1.3 建模方法选择

1.3.1 建模样品选择

原始光谱和化学值的准确性在很大程度上决定了近红外光谱分析结果的可靠性。由于待分析的样品包括不同规模化奶牛场牛舍集粪沟、挤奶车间粪水收集池、场区集污池到氧化塘全过程环节的样品，且每个环节的样品中都包含着养殖规模、牛群阶段、清粪方式、粪水处理工艺等众多影响因素信息，从采样到实验室光谱和化学值的测定结果可能会存在显著性差异，因此在建模前需要对异常样品进行剔除。本研究中采用蒙特卡洛交叉验证方法（Monte Carlo cross validation，MCCV）对异常样品进行剔除，图3为随机采样1000次，横向138个样本（TN模型）的预测均方根误差（RMSEP）统计分布图（均值MEAN和标准偏差STD）。设定阈值范围：均值为1500，标准偏差为300，剔除7、15、52、54、69、71和126共7个异常样品。采用同样的方法，从横向采集的138个样本（TP模型）中剔除了8个异常样品；纵向采集的51个样本（TN模型）中未发现异常样品；从纵向采集的51个样本（TP模型）中剔除了2个异常样品。

将剔除异常样品之后的样本集，分为校正集和预测集。校正集样品用来建立定标的数学模型，预测集样品用来验证模型的稳定性和准确性。为了使所建的定量分析模型能对规模化奶牛场粪水处理全过程环节粪水中TN和TP含量进行准确预测，校正集样品中必须包含不同奶牛场各个环节的代表性样品。本研究采用K-S（Kennard-Stone）方法来选择校正集和预测集样品，其结果如表2所示。

1.3.2 建模算法选择

在对现场规模化奶牛场各环节粪水中TN和TP含量进行检测时，模型算法越复杂，模型的适用性、稳定性和预测性越容易受到外界因素的影响。偏最小二乘法（PLS）是光谱多元校正最常用的一种方法，被广泛应用于近红外、拉曼、荧光等波谱定量模型的建立，是建立光谱定量校正模型的通用方法，并在现场检测中得到实际应用[24-25]。因此，本研究基于PLS算法来建立适合天津市规模化奶牛场粪水处理全过程环节中TN和TP含量的现场快速定量分析的全局数学模型。

1.3.3 建模预处理方法选择

图3 均值-方差分布图Figure 3 Mean-variance distribution

表2 K-S法划分的校正集和预测集样本信息Table 2 Classification information of samples in calibration set and prediction set

由于试验用粪水样品来源于规模化奶牛场粪水处理过程的各环节，其中包含固态粪渣等悬浮颗粒物质，是一种复杂的分析体系，对光具有较强的散射性。因此，原始的近红外漫反射光谱数据不仅包括了样本的化学信息，还包括了许多外界干扰信息，需要采用适合的预处理方法消除干扰因素的影响，以提高模型的可靠性、适用性和稳定性。对原始漫反射近红外光谱进行不同方法预处理，并采用PLS算法建立多变量回归模型，采用交叉验证均方根误差（RMSECV）来选择最佳建模的因子数（主成分数）。表3给出了不同预处理方法PLS回归结果。从决定系数R2、校正均方根误差（RMSEC）和预测均方根误差（RMSEP）来看，对于TN的近红外PLS回归模型最佳预处理方法为归一化方法；而对于TP的最佳预处理方法为SG平滑+基线校正。

1.3.4 建模波段选择

考虑到采用全波段4001个（12 000～4000 cm-1）变量建模，这些变量中包含与待测组分无关的冗余信息，同时又降低了建模效率。因此，本文采用间隔偏最小二乘法来选择建立粪水中TN和TP含量数学模型所需的有效波数变量。通过计算，选取4400～8800 cm-1（2201个波数变量）和4000～8000 cm-1（2001个波数变量）范围分别建立定量分析粪水中TN和TP含量的数学模型。

本文采用自行编写的PLS Matlab代码建立定量分析模型，采用The Unscrambler 9.7软件对光谱数据进行预处理，所有计算均采用MatlabR2017a软件工具（Mathwork Inc.）完成。

2 结果与讨论

2.1 规模化奶牛场粪水近红外漫反射光谱特性

图4为天津市23家规模化奶牛场粪水处理全过程环节采集的138个样品在4000～12 000 cm-1范围的近红外漫反射光谱。从中可观察到在5158、6938、8362 cm-1和10 232 cm-1处存在明显的吸收峰，这些峰都是由粪水中的水吸收所引起，其中5158 cm-1附近的波带主要是O-H反对称伸缩振动和弯曲振动的合频吸收引起[26]；6938 cm-1为O-H对称伸缩振动和反对称伸缩振动合频吸收引起；8362 cm-1附近波带为O-H反对称和对称伸缩振动，以及弯曲振动的合频吸收引起；10 232 cm-1处附近波带为倍频O-H对称伸缩振动和反对称伸缩振动的合频吸收引起。值得一提的是，虽然待分析样品来自养殖规模、清粪方式、粪污处理工艺组合和粪水流经环节等各不相同的奶牛场，但其近红外漫反射光谱在全波长范围内整体轮廓相似，仅在强度上存在差别，表明通过近红外漫反射光谱建立适合现场快速定量分析的数学模型具有可行性。

图4 奶牛场粪水样品的近红外漫反射光谱Figure 4 Near infrared diffuse reflection spectrum of slurry samples in dairy farms

表3 不同预处理方法PLS回归结果Table 3 PLS regression results of different pretreatment

2.2 粪水处理全过程环节模型的建立

2.2.1 主成分分析

为明确同一奶牛场纵向粪水处理路径中各环节粪水样品之间的关系和性质组成的变化，以及对定量分析模型的影响，对滨海新区某典型规模化奶牛场51个粪水样品的近红外漫反射光谱进行主成分分析。在主成分得分图上，样品之间的距离越近，表明这些样品的性质、组成越相似；距离越远，表明样品间性质、组分变化越大。因此，通过得分图上样品分布，可推断出粪水样品有机组分随处理环节的变化差异明显。

图5为前两个主成分的得分图，其中第一主成分PC1（Principal component 1）解释84%的光谱总变量，第二主成分PC2解释14%的光谱总变量。图中粉色椭圆所包含的是集污池样品，蓝色和红色椭圆所包含的是氧化塘样品，黑色椭圆所包含的是奶厅附近收集池样品。集污池样品主要分布在得分图的右上角区域，氧化塘样品主要分布在右下角区域，奶厅附近收集池样品主要分布在左下角区域，因此从得分图上可明确样品来源的环节，以及样品之间的组分、浓度等特性关联。同时可观察到，奶厅附近收集池样品与集污池、氧化塘样品距离较远，表明样品中有机物组成随处理环节发生了较大变化，这与奶厅收集池中的废水浓度远低于牛舍集粪沟中的粪水浓度；奶厅废水与牛舍粪水混合后水体中的物质组成和浓度值变化较大；混合水依次流经集污池、分离池、氧化塘的过程中水质组分和浓度呈梯度变化等现场情况相符。2.2.2定量分析模型的建立

图5 51个样品前两个主成分得分图Figure 5 Scores plot for 51 samples on first two principal components

对1.3所处理的光谱数据，建立定量分析纵向奶牛场粪水处理过程中各环节粪水样品TN和TP含量的PLS数学模型，交叉验证均方根误差（RMSECV）分别为181.97和4.72。将预测集未知样品光谱数据（18×2201）输入所建立的PLS模型中，得到各样品对应的TN含量。图6a为所建立的TN模型对预测集18个未知样品（8个收集池样品、6个氧化塘样品、4个奶厅附近收集池样品）TN含量预测结果与实际含量的线性拟合，其拟合关系为：C预测=1.02C实际+50.4，拟合相关系数R为0.96，预测均方根误差RMSEP为187.80。同理将预测集未知样品光谱数据（17×2001）输入所建立的PLS模型中，得到各样品对应的TP含量。图6b为所建立的TP模型对预测集17个未知样品（7个收集池样品，7个氧化塘样品，3个奶厅附近收集池样品）TP含量预测结果与实际含量的线性拟合，其拟合关系为：C预测=0.64C实际+7.3，拟合相关系数R为0.91，预测均方根误差RMSEP为3.59。上述结果表明，基于近红外漫反射光谱实现纵向同一奶牛场粪水处理过程中各环节样品TN和TP的检测是可行的。

图6 纵向PLS模型预测含量与实际含量的线性拟合Figure 6 Linear fitting between predicted and practical content using vertical PLS model

2.3 粪水处理全局模型的建立

2.3.1 主成分分析

为明确不同奶牛场粪水处理过程中粪水样品之间的关联和特性差异，以及对于定量分析模型的影响，对天津市23家规模化奶牛场在不同的饲料投入（饲料配方和饲喂量）、饲喂方式、清粪方式、处理工艺及环节等因素作用下共138个粪水样品的近红外漫反射光谱进行主成分分析。图7为前两个主成分的得分图，前两个主成分分别解释光谱总变量的64%和33%。从整体样品分布来看，同一奶牛场的多数样品距离接近，图中仅给出了5个椭圆来说明所包含的样品来源于同一奶牛场，其中包含样品数最多的为2.2节分析的滨海新区某一奶牛场（纵向采样对象），94号和95号样品来源于另外一家奶牛场，该场仅采集了这两个样品。上述结果说明：相对于纵向（同一奶牛场）不同环节，横向（不同奶牛场）采集的样品之间差异性更大。

2.3.2 定量分析模型的建立

图7 138个样品前两个主成分得分图Figure 7 Scores plot for 138 samples on first two principal components

对1.3所处理的光谱数据，建立定量分析横向不同奶牛场粪水处理过程中各环节样品TN和TP含量的PLS数学模型，其RMSECV分别为355.74和8.13。图8a为所建立的TN模型对预测集41个未知样品TN含量预测结果与实际含量的线性拟合，其拟合关系为：C预测=0.91C实际-0.84，拟合相关系数R为0.96，预测均方根误差RMSEP为238.59。图8b为所建立的TP模型对预测集38个未知样品TP含量预测结果与实际含量的线性拟合，其拟合关系为：C预测=0.8C实际+7.12，拟合相关系数R为0.91，预测均方根误差RMSEP为6.56。上述结果表明，基于近红外漫反射光谱实现横向不同奶牛场粪水处理过程中各环节样品TN和TP的检测也是完全可行的。

2.4 影响因素分析

从纵向和横向样品得分图可以看出，不同环节和不同奶牛场的粪水分布在不同区域。因此，定标模型中不同环节和不同奶牛场粪水样品的分布对定量分析模型的适用性、稳定性和预测能力都存在影响。对于纵向和横向定量分析模型，定标模型中应包括全程各环节代表性的样品，这样所建立的模型才能实现全程和全局未知样品的准确预测，上述建立的全程和全局模型取得的理想预测结果恰好验证这一点。同时注意到，图5第一主成分提取了84%的总变量信息，而图7第一主成分仅提取了64%的总变量信息，进一步表明了相对于纵向样品，横向样品变异性更大。

图8 横向PLS模型预测含量与实际含量的线性拟合Figure 8 Linear fitting between predicted and practical content using horizontal PLS model