金 晶， 毛 星， 戴红君， 刘 杨， 任 妮
(江苏省农业科学院信息中心，江苏 南京 210014)
中华绒螯蟹(Eriocheirsinensis)，俗名大闸蟹，是中国淡水养殖的重要水产品种，因其风味独特、营养丰富而受到消费者的喜爱[1-2]。2020年全国中华绒螯蟹养殖面积超过6.7×105hm2，总产量达到7.8×105t[3-4]。中华绒螯蟹养殖产业已成为渔业养殖的重要收入来源，对乡村振兴战略的实施具有积极意义[5]。中华绒螯蟹养殖对环境要求较高，在高密度养殖的背景下，投喂方式直接影响养殖对象的产量和品质，并可能对养殖区域的生态环境造成威胁[6-7]。因此，明确中华绒螯蟹的饲料投喂对养殖水质的影响，能够降低养殖风险，对水质的及时调节和中华绒螯蟹的标准化、绿色化生产具有重要意义。
饲料投喂对中华绒螯蟹个体生长和养殖水体水质有重要影响。目前，中华绒螯蟹养殖主要用杂鱼、玉米、小麦、大豆和油饼等传统饲料投喂[8]，但由于供应不稳定、营养不均衡、携带病原体等原因，这种投喂方式会引起池塘水质恶化、中华绒螯蟹疾病暴发等严重问题[9-11]。当投喂量不足时，中华绒螯蟹的生长受限，经济效益降低[12]；
当投喂量过高时，中华绒螯蟹无法全部摄食，残留的饵料会腐败变质，造成水体污染和中华绒螯蟹生长受限等负面影响[13]。水温、溶解氧质量浓度、pH值等作为评价水体质量的重要指标[4]，对中华绒螯蟹养殖有重要影响。其中，针对溶解氧的研究最为广泛，其与中华绒螯蟹生存和饲料利用息息相关[14-15]。也有研究关注叶绿素、pH值和水温等其他参数[16-18]。另外，许多水质参数之间存在相互作用[19]，例如氮、磷和有机质的养分负荷量过多时，会增大需氧量和有毒代谢物浓度，导致溶解氧含量降低和水体富营养化[20]。
目前，已有研究者提出了基于溶解氧、pH值等参数的投喂策略[21-22]，分析了不同饲料配比对中华绒螯蟹生长的影响[10-11]，但关于投喂对水质影响的研究相对较少。王海候等[6]分析了螺蛳投喂量对中华绒螯蟹产量的影响，发现投喂量过多或过少均对中华绒螯蟹生长不利。针对对虾养殖的研究则发现，过高的投喂频率不但无法显著提高对虾的饲料利用效率，反而污染了水体环境[7]；
饶毅等[23]认为生物絮团技术显著改善了草鱼养殖水质。上述研究多是通过控制试验的方法分析不同投喂策略、投喂频率对养殖水质的影响，但当前中国的中华绒螯蟹养殖仍处在经验养殖阶段，普通农户的投喂策略以经验判断为主，实际生产环境下的饲料投喂对养殖水质的影响尚不明晰。此外，中华绒螯蟹投喂行为对水质的影响存在一定的滞后效应，而现有研究大多未考虑滞后效应。Gasparrini等[24]以广义线性模型和广义加法模型等传统模型的思想为基础，阐述了分布滞后非线性模型(Distributed lag non-linear models, DLNM)理论，适用于投喂与水质变化之间的滞后效应研究。
本研究基于2020年江苏省宜兴市高塍镇中华绒螯蟹水产养殖试验示范基地中2个试验池塘的水质、气象和养殖投喂数据，在考虑滞后效应的基础上，采用分布滞后非线性模型，以溶解氧、pH值、氨氮等6个指标表征养殖水质情况，分析饲料投喂后养殖水质变化，以探讨实际生产环境下中华绒螯蟹饲料投喂量对养殖水质的影响，为中华绒螯蟹养殖过程中水体环境的及时调节和中华绒螯蟹养殖产业的可持续发展提供参考。
1.1 数据采集
中华绒螯蟹是江苏省宜兴市的特色水产养殖品种，具有养殖面积大、产量高、效益好等特点。本研究所用数据采集自2020年江苏省宜兴市高塍镇中华绒螯蟹水产养殖试验示范基地的2个试验池塘，分别由2个养殖户经营管理(图1)，其中1号池塘面积为2.17 hm2，2号池塘面积为0.94 hm2。
蟹塘水质数据来自池塘中布设的传感器。在各池塘选取4个位置，取距离水底30 cm、60 cm和90 cm处作为布设点，每个池塘分别布设12个传感器(图1，表1)，用于监测溶解氧质量浓度、叶绿素质量浓度、浊度、氨氮质量浓度、水温、pH和电导率。数据采集间隔时间为10 min，2020年中华绒螯蟹生长季数据量总计1 602 400条。气象数据来源于试验池塘附近建立的小型气象站，用于监测大气温度、大气压强、空气湿度、雨量、太阳辐射、风速和风向，数据采集间隔时间为1 min，数据量总计481 456条。
中华绒螯蟹投喂数据为2个试验池塘的养殖户通过“与蟹同行”小程序上传获得，每条数据均有专业人员进行人工审核，确保数据的可靠性。1号池塘共有201条投喂记录，2号池塘共有207条投喂记录，涉及投喂的饲料类型、投喂时间、投喂量、饲料的采购价格、采购渠道，以及人工审核情况。
1.2 数据融合与处理
由于投喂数据在采集时已进行人工审核，因此数据清洗主要针对气象和水质数据。根据各类设备采集数据的参数范围(表2)，判断气象观测和水质监测时序数据中的异常值，并识别可能由于网络波动、设备异常等原因导致的长时间数值唯一且异常的情况，将数据集中的异常值作为缺失值对待。由于水产传感器的稳定性较差，而水质指标往往受到气象、投喂、增氧泵设备运转等多种因素影响，直接对缺失值进行填补可能会对模型结果带来较大不确定性，因此研究中只选择连续、没有缺失的数据集进行模型构建。在后续研究中，气象和水质观测数据均通过取平均值的方法统一至以天作为时间分辨率。
表2 传感器和小型气象站设备型号及参数
1.3 DLNM模型的构建
选取溶解氧质量浓度、叶绿素质量浓度、浊度、氨氮质量浓度、pH和电导率6个指标表征水体环境情况，对每个池塘中12个传感器监测的每组水质数据分别进行分析。考虑到观测的7个气象因素之间可能存在多重共线性，且各水质指标并非均与所有气象因素相关，为了保证DLNM模型构建时参数选择的合理性，先建立各水质参数与气象因素的最小二乘线性回归模型，即：
yindex=β1xtemperature+β2xpressure+β3xradiation+β4xwindspeed+β5xwinddirection+β6xprecipitation+β7xhumidity+β0+e
式中，yindex表示水质参数，xtemperature表示大气温度，xpressure表示气压，xradiation表示太阳辐射，xwindspeed表示风速，xwinddirection表示风向，xprecipitation表示雨量，xhumidity表示大气湿度，β0～β7表示回归参数，e表示残差。
通过逐步回归方法，基于AIC准则选取气象因素组合，并基于方差膨胀系数进行共线性检验，确保方差膨胀系数小于10，在排除多重共线性的同时保证模型具有较好的解释性，进而确定引入各个DLNM模型的气象指标。
在DLNM模型构建中，根据已有研究结果，投喂量数据的基函数选择自然样条函数，滞后维度的基函数选择多项式函数，气象参数选用自由度为3的自然样条函数进行拟合[25-26]。通过分别拟合7～40期不同滞后天数的DLNM模型并计算改进的赤池信息准则(Quasi-Akaike information criterion，QAIC)和改进的贝叶斯信息准则(Quasi-Bayesian information criterion，QBIC)2个指标，选择其中指标值最小时所对应的滞后期作为DLNM模型的最大滞后天数。以溶解氧参数为例，基本模型如下：
yDissolve_oxgen_t=α+β1Feedt,l+ns(Humidity,3)+ns(Temperature,3)+ns(Pressure,3)+ns(WindSpeed,3)+ns(WindDirection,3)+ns(Radiation,3)
式中，t表示观察时间(d),yDissolve_oxgen_t是当日的溶解氧监测数值，α是模型截距项，β1为日投喂量的交叉矩阵系数，Feedt,l为利用每日投喂量建立的交叉基，ns为自然样条，3为自由度。气象指标为筛选出的与溶解氧相关且不存在多重共线性的变量。建模过程中，参考谷少华等[27]的方法，各DLNM模型均选取投喂量的中位数作为参考点，分别以投喂量的99%和1%分位数表征高和低投喂量情景。
数据的处理分析均在R 4.1.0环境下进行。线性回归分析基于lm程序包，广义线性模型构建基于glm程序包，DLNM模型构建基于dlnm程序包[28]，时间序列数据的处理基于tseries和lubridate程序包。
2.1 气象参数的选择
根据最小二乘回归模型拟合结果，选择引入各水质参数DLNM建模过程的气象指标(表3)。
表3 各水质参数相关的气象指标
2.2 饲料投喂量对蟹塘水质的影响
2.2.1 饲料投喂与DLNM建模数据概况 2020年生长季中，2个中华绒螯蟹养殖试验池塘的主要投喂饲料类型包括冰鱼、蛋白质饲料和玉米，其中1号池塘全生长季3种饲料的投喂量占比分别为87.8%、9.0%和3.2%；
2号池塘3种饲料投喂量占比分别为72.0%、23.2%和4.8%，另有少量的复合型益生菌。
在DLNM模型建模过程中，不同观测点引入建模的数据量各不相同。每个池塘所观测的6项水质指标共可生成21组数据。1号池塘分析用最小连续观测量为55条，最大值为113条；
2号池塘最小连续观测量为63条，最大值为168条，符合大样本要求。1号池塘日投喂量参考点为55 kg/hm2，99%分位点为120 kg/hm2，1%分位点为7 kg/hm2。2号池塘日投喂量参考点为64 kg/hm2，99%分位点为181 kg/hm2，1%分位点为9 kg/hm2。
2.2.2 饲料投喂量对蟹塘各水质指标变化的影响 由于各池塘0号点位中层(即0-2号观测点)监测了所有的水质指标，因此以该点对比分析各水质指标受投喂量的影响。对于溶解氧质量浓度，1号池塘均值为7.96 mg/L，标准差为2.60 mg/L，DLNM模型滞后期为36 d；
2号池塘溶解氧质量浓度均值为7.66 mg/L，标准差为2.08 mg/L，滞后期为38 d，表明投喂量对溶解氧质量浓度的影响具有长期滞后效应(图2)。在30阶滞后情况下，相对于参考值，1号和2号池塘不同投喂量对于溶解氧质量浓度的变化影响均不显著(P>0.05)。当1号池塘投喂量由参考值降为7 kg/hm2时，在22阶滞后期使溶解氧质量浓度升高了0.614 mg/L(95%置信区间：0.001～1.228 mg/L)，该正效应持续至第25阶，效应值为0.498 mg/L(95%置信区间：0.024～0.972 mg/L)，其余滞后期不显著。当2号池塘投喂量由参考值降为9 kg/hm2时，在10～19阶滞后具有使溶解氧质量浓度降低的显著负效应，其中最大负效应为-1.063 mg/L(95%置信区间：-1.917～-0.208 mg/L)，但在高投喂量情景下各滞后期效应均不显著。
阴影和虚线区间部分为95%置信区间。图2 基于DLNM模型的各池塘0-2观测点投喂量与溶解氧质量浓度变化的滞后效应Fig.2 Lag effects of feeding amount of feedstuff on the mass concentration variation of dissolved oxygen in each crab pond based on DLNM model at 0-2 observation point
对于pH值，1号池塘平均值为7.94，标准差为0.23；
2号池塘平均值为8.40，标准差为0.48。在30阶滞后情况下，相对于参考值，2号池塘仅在投喂量为65 kg/hm2及76～141 kg/hm2时对pH值变化的效应不显著，在其余投喂量下均有显著负效应(图3)。对于2号池塘，低投喂量情景下滞后40 d的负效应值最大，为-0.17(95%置信区间：-0.264～-0.076)；
高投喂量情景下，滞后9～36 d均显示出负效应，最大负效应为-0.096(95%置信区间：-0.159～-0.032)。1号池塘无论是在特定滞后期或是特定投喂量情景下效应均不显著。
阴影和虚线区间部分为95%置信区间。图3 基于DLNM模型的各池塘0-2观测点投喂量与pH值变化的滞后效应Fig.3 Lag effects of feeding amount of feedstuff on the pH value in each crab pond based on DLNM model at 0-2 observation point
在高和低投喂量情景下，不同水质参数对投喂的响应模式不同(表4)。对于溶解氧质量浓度，在高投喂量情景下，各池塘在全滞后期无显著效应；
在低投喂量情景下，1号池塘在36期滞后中的4期(11.1%)具有显著正效应，2号池塘则在38期滞后中的10期(26.3%)具有显著负效应。对于pH值，2号池塘在2种情景下均呈现显著负效应。对于氨氮质量浓度，1号池塘在高投喂量情景下出现5期显著负效应，在低投喂量情景下2个池塘均在部分滞后期表现为显著正效应，表明当投喂量较低时，水体的氨氮质量浓度相对于参考值有显著上升。2号池塘水体浊度和电导率受投喂行为的影响相对明显，高、低投喂量情景下分别有11期(32.4%)和15期(44.1%)使水体浊度显著上升，分别有26期(66.7%)和14期(35.9%)使电导率显著下降。2种情景下饲料投喂量对叶绿素质量浓度变化均无显著效应。
表4 各池塘0-2观测点水质在高、低投喂量情景下的滞后效应分析
由于在各池塘布设的12个观测点均有对溶解氧质量浓度的监测，选用模型结果更显著的2号池塘分析不同观测点溶解氧质量浓度在各投喂量下的滞后效应(图4)。1-3和3-3 2个观测点模型结果均不显著，无统计学意义。其余10个观测点中，由于饲料投喂对溶解氧质量浓度带来的显著正效应较少，负效应较多。当投喂量降低至9 kg/hm2时，7个观测点部分滞后期的溶解氧质量浓度指标呈显著负效应，且多出现在滞后期前中期(如0-1点)和后期(如3-1点)。其中，1-2号点位的10阶滞后负效应最大，使得溶解氧质量浓度下降1.731 mg/L(95%置信区间：-2.647～-0.816 mg/L)。投喂量升高至181 kg/hm2时，9个观测点部分滞后期的溶解氧质量浓度指标表现出显著负效应，该负效应多出现在前期和后期(如0-1点)。其中，滞后期数最多的点位为0-1，在5～10期和30～38期共15期滞后，使溶解氧质量浓度显著降低。在所有观测点中，2-2点位的4阶滞后期使溶解氧质量浓度下降的效应最大，下降了1.781 mg/L(95%置信区间：-2.916～-0.647 mg/L)。总体而言，同一深度不同位置及同一位置不同深度的溶解氧质量浓度滞后效应模式也有所不同，变化趋势没有明显规律。
阴影和虚线区间部分为95%置信区间。图4 基于DLNM模型的2号池塘各观测点投喂量与溶解氧质量浓度变化的滞后效应Fig.4 Lag effects of feeding amount of feedstuff on the mass concentration of dissolved oxygen in pond 2 based on DLNM model at each observation point
本研究采用分布滞后非线性模型，在考虑滞后效应的基础上，基于2020年宜兴市中华绒螯蟹水产养殖试验示范基地的2个试验池塘数据，分析了生产环境下饲料投喂量对蟹塘水质的影响。研究结果可为生产环境下蟹塘水质的影响因素分析及水质调节策略的及时制定提供支撑。
研究发现，在高、低两种投喂量下，饲料投喂与蟹塘溶解氧质量浓度、pH、氨氮质量浓度、电导率、浊度的变化之间存在长期滞后效应，但与叶绿素质量浓度变化之间的关系不显著。在高投喂量情景下，相对于参考值，总体表现出氨氮质量浓度、电导率、pH值在部分滞后期显著减少、浊度显著上升的趋势，其中溶解氧质量浓度在75%的观测点显著减少，滞后期内最大负效应达-1.781 mg/L(95%置信区间：-2.916～-0.647 mg/L)。已有研究者指出，中华绒螯蟹饲料投喂量较高时，饲料无法充分利用，残饵、粪便等使得水体浊度显著上升，随着有机物质分解和中华绒螯蟹摄食行为的增加，耗氧量和二氧化碳的排放量有所增加，导致溶解氧质量浓度和pH值显著降低[13,20]，这与本研究结论相一致。1号池塘的氨氮质量浓度和2号池塘的电导率在后期显著减小，浊度在部分滞后期也表现出显著的负效应，这可能是由于养殖户在投喂量过大时会定期进行池塘改底等措施。试验池塘的用药记录显示，在中华绒螯蟹养殖过程中会使用特效芽孢杆菌、底居安、碧水安等药物降解水体中的有机物、氨氮等有害物质，从而使池塘氨氮质量浓度在投喂后15～20 d、电导率在投喂后30～40 d、浊度在投喂后26～31 d左右显著下降。在低投喂量情景下，相对于参考值，氨氮质量浓度和浊度在部分滞后期显著增加，pH值显著减小。这与饵料不足时，螃蟹为寻觅食物在池塘底部淤泥处大量活动有关[29]。另外，由于投喂量减少，残饵量相对减少，因此2号池塘的电导率在多个滞后期呈显著负效应。1号池塘的溶解氧质量浓度在部分滞后期显著上升，但2号池塘在7个观测点均出现对溶解氧质量浓度的显著负效应, 滞后期内最大负效应达-1.731 mg/L(95%置信区间：-2.647～-0.816 mg/L)；
pH值的负效应多出现在滞后的后期，在滞后期前期虽然估计效应值并不显著，但相比于高投喂量情景效应值较小。这可能是由于投喂量减少后，前期溶解氧质量浓度上升、pH值增高，养殖户基于较好的水质情况调整了投喂策略，增加了投喂量，使得后期又出现水质下降的趋势。
研究发现，同一深度不同位置及同一位置不同深度的溶解氧质量浓度受投喂量影响的滞后效应模式有所不同，多数观测点滞后期在30 d以上。这与溶解氧在水体不同位置质量浓度分布有所不同有关。现有研究结果表明，养殖池塘溶解氧质量浓度受到多种生态因素影响，并在水体中呈现出三维分布特征[30-31]。具体表现为在池塘同一水平面上不同位置的溶解氧质量浓度有所不同，在同一垂直面上不同水层的溶解氧质量浓度也不同[4, 32]，这与本研究在不同观测点观测到的滞后效应模式有所不同的结果相符合。
本研究结果可辅助生产环境下蟹塘水质调节策略的制定。对于养殖户，在投喂大量饲料后，宜着重关注投喂后5～10 d和20 d后的溶解氧质量浓度下降情况，以及投喂9 d后的pH值下降情况；
在投喂少量饲料后，宜关注投喂后8～12 d和30 d后的溶解氧质量浓度下降情况、投喂后10～18 d pH值下降情况、投喂后10 d内氨氮质量浓度和投喂后5～20 d池塘浊度变化，及时通过开启增氧泵设备、改变投喂策略、池塘改底等方式调节养殖水质。本研究也存在一定的局限性，通过DLNM模型的构建，可以给出投喂和水质参数之间暴露-反应关系的总结，但由于未报告各滞后期、投喂量对应影响效应的统计显著性[28]，因此本研究采用了特定投喂量和滞后期进行分析。另外，由于中华绒螯蟹养殖水质受到气象、生态、养殖经验、投喂、用药等多方面影响，且同一池塘不同位置的水质变化情况有所不同，未来还需要融合更多类型的数据并将其集成在分析过程中，进一步明确蟹塘养殖水质的复杂变化规律。