基于分布滯后非線性模型的蟹塘水質受投喂量影響的關系分析

金 晶， 毛 星， 戴紅君， 劉 楊， 任 妮

(江蘇省農業(yè)科學院信息中心，江蘇 南京 210014)

中華絨螯蟹(Eriocheirsinensis)，俗名大閘蟹，是中國淡水養(yǎng)殖的重要水產品種，因其風味獨特、營養(yǎng)豐富而受到消費者的喜愛[1-2]。2020年全國中華絨螯蟹養(yǎng)殖面積超過6.7×105hm2，總產量達到7.8×105t[3-4]。中華絨螯蟹養(yǎng)殖產業(yè)已成為漁業(yè)養(yǎng)殖的重要收入來源，對鄉(xiāng)村振興戰(zhàn)略的實施具有積極意義[5]。中華絨螯蟹養(yǎng)殖對環(huán)境要求較高，在高密度養(yǎng)殖的背景下，投喂方式直接影響?zhàn)B殖對象的產量和品質，并可能對養(yǎng)殖區(qū)域的生態(tài)環(huán)境造成威脅[6-7]。因此，明確中華絨螯蟹的飼料投喂對養(yǎng)殖水質的影響，能夠降低養(yǎng)殖風險，對水質的及時調節(jié)和中華絨螯蟹的標準化、綠色化生產具有重要意義。

飼料投喂對中華絨螯蟹個體生長和養(yǎng)殖水體水質有重要影響。目前，中華絨螯蟹養(yǎng)殖主要用雜魚、玉米、小麥、大豆和油餅等傳統(tǒng)飼料投喂[8]，但由于供應不穩(wěn)定、營養(yǎng)不均衡、攜帶病原體等原因，這種投喂方式會引起池塘水質惡化、中華絨螯蟹疾病暴發(fā)等嚴重問題[9-11]。當投喂量不足時，中華絨螯蟹的生長受限，經濟效益降低[12]；當投喂量過高時，中華絨螯蟹無法全部攝食，殘留的餌料會腐敗變質，造成水體污染和中華絨螯蟹生長受限等負面影響[13]。水溫、溶解氧質量濃度、pH值等作為評價水體質量的重要指標[4]，對中華絨螯蟹養(yǎng)殖有重要影響。其中，針對溶解氧的研究最為廣泛，其與中華絨螯蟹生存和飼料利用息息相關[14-15]。也有研究關注葉綠素、pH值和水溫等其他參數(shù)[16-18]。另外，許多水質參數(shù)之間存在相互作用[19]，例如氮、磷和有機質的養(yǎng)分負荷量過多時，會增大需氧量和有毒代謝物濃度，導致溶解氧含量降低和水體富營養(yǎng)化[20]。

目前，已有研究者提出了基于溶解氧、pH值等參數(shù)的投喂策略[21-22]，分析了不同飼料配比對中華絨螯蟹生長的影響[10-11]，但關于投喂對水質影響的研究相對較少。王海候等[6]分析了螺螄投喂量對中華絨螯蟹產量的影響，發(fā)現(xiàn)投喂量過多或過少均對中華絨螯蟹生長不利。針對對蝦養(yǎng)殖的研究則發(fā)現(xiàn)，過高的投喂頻率不但無法顯著提高對蝦的飼料利用效率，反而污染了水體環(huán)境[7]；饒毅等[23]認為生物絮團技術顯著改善了草魚養(yǎng)殖水質。上述研究多是通過控制試驗的方法分析不同投喂策略、投喂頻率對養(yǎng)殖水質的影響，但當前中國的中華絨螯蟹養(yǎng)殖仍處在經驗養(yǎng)殖階段，普通農戶的投喂策略以經驗判斷為主，實際生產環(huán)境下的飼料投喂對養(yǎng)殖水質的影響尚不明晰。此外，中華絨螯蟹投喂行為對水質的影響存在一定的滯后效應，而現(xiàn)有研究大多未考慮滯后效應。Gasparrini等[24]以廣義線性模型和廣義加法模型等傳統(tǒng)模型的思想為基礎，闡述了分布滯后非線性模型(Distributed lag non-linear models, DLNM)理論，適用于投喂與水質變化之間的滯后效應研究。

本研究基于2020年江蘇省宜興市高塍鎮(zhèn)中華絨螯蟹水產養(yǎng)殖試驗示范基地中2個試驗池塘的水質、氣象和養(yǎng)殖投喂數(shù)據(jù)，在考慮滯后效應的基礎上，采用分布滯后非線性模型，以溶解氧、pH值、氨氮等6個指標表征養(yǎng)殖水質情況，分析飼料投喂后養(yǎng)殖水質變化，以探討實際生產環(huán)境下中華絨螯蟹飼料投喂量對養(yǎng)殖水質的影響，為中華絨螯蟹養(yǎng)殖過程中水體環(huán)境的及時調節(jié)和中華絨螯蟹養(yǎng)殖產業(yè)的可持續(xù)發(fā)展提供參考。

1 材料與方法

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

中華絨螯蟹是江蘇省宜興市的特色水產養(yǎng)殖品種，具有養(yǎng)殖面積大、產量高、效益好等特點。本研究所用數(shù)據(jù)采集自2020年江蘇省宜興市高塍鎮(zhèn)中華絨螯蟹水產養(yǎng)殖試驗示范基地的2個試驗池塘，分別由2個養(yǎng)殖戶經營管理(圖1)，其中1號池塘面積為2.17 hm2，2號池塘面積為0.94 hm2。

蟹塘水質數(shù)據(jù)來自池塘中布設的傳感器。在各池塘選取4個位置，取距離水底30 cm、60 cm和90 cm處作為布設點，每個池塘分別布設12個傳感器(圖1，表1)，用于監(jiān)測溶解氧質量濃度、葉綠素質量濃度、濁度、氨氮質量濃度、水溫、pH和電導率。數(shù)據(jù)采集間隔時間為10 min，2020年中華絨螯蟹生長季數(shù)據(jù)量總計1 602 400條。氣象數(shù)據(jù)來源于試驗池塘附近建立的小型氣象站，用于監(jiān)測大氣溫度、大氣壓強、空氣濕度、雨量、太陽輻射、風速和風向，數(shù)據(jù)采集間隔時間為1 min，數(shù)據(jù)量總計481 456條。

中華絨螯蟹投喂數(shù)據(jù)為2個試驗池塘的養(yǎng)殖戶通過“與蟹同行”小程序上傳獲得，每條數(shù)據(jù)均有專業(yè)人員進行人工審核，確保數(shù)據(jù)的可靠性。1號池塘共有201條投喂記錄，2號池塘共有207條投喂記錄，涉及投喂的飼料類型、投喂時間、投喂量、飼料的采購價格、采購渠道，以及人工審核情況。

1.2 數(shù)據(jù)融合與處理

由于投喂數(shù)據(jù)在采集時已進行人工審核，因此數(shù)據(jù)清洗主要針對氣象和水質數(shù)據(jù)。根據(jù)各類設備采集數(shù)據(jù)的參數(shù)范圍(表2)，判斷氣象觀測和水質監(jiān)測時序數(shù)據(jù)中的異常值，并識別可能由于網絡波動、設備異常等原因導致的長時間數(shù)值唯一且異常的情況，將數(shù)據(jù)集中的異常值作為缺失值對待。由于水產傳感器的穩(wěn)定性較差，而水質指標往往受到氣象、投喂、增氧泵設備運轉等多種因素影響，直接對缺失值進行填補可能會對模型結果帶來較大不確定性，因此研究中只選擇連續(xù)、沒有缺失的數(shù)據(jù)集進行模型構建。在后續(xù)研究中，氣象和水質觀測數(shù)據(jù)均通過取平均值的方法統(tǒng)一至以天作為時間分辨率。

表2 傳感器和小型氣象站設備型號及參數(shù)

1.3 DLNM模型的構建

選取溶解氧質量濃度、葉綠素質量濃度、濁度、氨氮質量濃度、pH和電導率6個指標表征水體環(huán)境情況，對每個池塘中12個傳感器監(jiān)測的每組水質數(shù)據(jù)分別進行分析?？紤]到觀測的7個氣象因素之間可能存在多重共線性，且各水質指標并非均與所有氣象因素相關，為了保證DLNM模型構建時參數(shù)選擇的合理性，先建立各水質參數(shù)與氣象因素的最小二乘線性回歸模型，即：

yindex=β1xtemperature+β2xpressure+β3xradiation+β4xwindspeed+β5xwinddirection+β6xprecipitation+β7xhumidity+β0+e

式中，yindex表示水質參數(shù)，xtemperature表示大氣溫度，xpressure表示氣壓，xradiation表示太陽輻射，xwindspeed表示風速，xwinddirection表示風向，xprecipitation表示雨量，xhumidity表示大氣濕度，β0～β7表示回歸參數(shù)，e表示殘差。

通過逐步回歸方法，基于AIC準則選取氣象因素組合，并基于方差膨脹系數(shù)進行共線性檢驗，確保方差膨脹系數(shù)小于10，在排除多重共線性的同時保證模型具有較好的解釋性，進而確定引入各個DLNM模型的氣象指標。

在DLNM模型構建中，根據(jù)已有研究結果，投喂量數(shù)據(jù)的基函數(shù)選擇自然樣條函數(shù)，滯后維度的基函數(shù)選擇多項式函數(shù)，氣象參數(shù)選用自由度為3的自然樣條函數(shù)進行擬合[25-26]。通過分別擬合7～40期不同滯后天數(shù)的DLNM模型并計算改進的赤池信息準則(Quasi-Akaike information criterion，QAIC)和改進的貝葉斯信息準則(Quasi-Bayesian information criterion，QBIC)2個指標，選擇其中指標值最小時所對應的滯后期作為DLNM模型的最大滯后天數(shù)。以溶解氧參數(shù)為例，基本模型如下：

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是當日的溶解氧監(jiān)測數(shù)值，α是模型截距項，β1為日投喂量的交叉矩陣系數(shù)，F(xiàn)eedt,l為利用每日投喂量建立的交叉基，ns為自然樣條，3為自由度。氣象指標為篩選出的與溶解氧相關且不存在多重共線性的變量。建模過程中，參考谷少華等[27]的方法，各DLNM模型均選取投喂量的中位數(shù)作為參考點，分別以投喂量的99%和1%分位數(shù)表征高和低投喂量情景。

數(shù)據(jù)的處理分析均在R 4.1.0環(huán)境下進行。線性回歸分析基于lm程序包，廣義線性模型構建基于glm程序包，DLNM模型構建基于dlnm程序包[28]，時間序列數(shù)據(jù)的處理基于tseries和lubridate程序包。

2 結果與分析

2.1 氣象參數(shù)的選擇

根據(jù)最小二乘回歸模型擬合結果，選擇引入各水質參數(shù)DLNM建模過程的氣象指標(表3)。

表3 各水質參數(shù)相關的氣象指標

2.2 飼料投喂量對蟹塘水質的影響

2.2.1 飼料投喂與DLNM建模數(shù)據(jù)概況 2020年生長季中，2個中華絨螯蟹養(yǎng)殖試驗池塘的主要投喂飼料類型包括冰魚、蛋白質飼料和玉米，其中1號池塘全生長季3種飼料的投喂量占比分別為87.8%、9.0%和3.2%；2號池塘3種飼料投喂量占比分別為72.0%、23.2%和4.8%，另有少量的復合型益生菌。

在DLNM模型建模過程中，不同觀測點引入建模的數(shù)據(jù)量各不相同。每個池塘所觀測的6項水質指標共可生成21組數(shù)據(jù)。1號池塘分析用最小連續(xù)觀測量為55條，最大值為113條；2號池塘最小連續(xù)觀測量為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號觀測點)監(jiān)測了所有的水質指標，因此以該點對比分析各水質指標受投喂量的影響。對于溶解氧質量濃度，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%置信區(qū)間：0.001～1.228 mg/L)，該正效應持續(xù)至第25階，效應值為0.498 mg/L(95%置信區(qū)間：0.024～0.972 mg/L)，其余滯后期不顯著。當2號池塘投喂量由參考值降為9 kg/hm2時，在10～19階滯后具有使溶解氧質量濃度降低的顯著負效應，其中最大負效應為-1.063 mg/L(95%置信區(qū)間：-1.917～-0.208 mg/L)，但在高投喂量情景下各滯后期效應均不顯著。

陰影和虛線區(qū)間部分為95%置信區(qū)間。圖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%置信區(qū)間：-0.264～-0.076)；高投喂量情景下，滯后9～36 d均顯示出負效應，最大負效應為-0.096(95%置信區(qū)間：-0.159～-0.032)。1號池塘無論是在特定滯后期或是特定投喂量情景下效應均不顯著。

陰影和虛線區(qū)間部分為95%置信區(qū)間。圖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

在高和低投喂量情景下，不同水質參數(shù)對投喂的響應模式不同(表4)。對于溶解氧質量濃度，在高投喂量情景下，各池塘在全滯后期無顯著效應；在低投喂量情景下，1號池塘在36期滯后中的4期(11.1%)具有顯著正效應，2號池塘則在38期滯后中的10期(26.3%)具有顯著負效應。對于pH值，2號池塘在2種情景下均呈現(xiàn)顯著負效應。對于氨氮質量濃度，1號池塘在高投喂量情景下出現(xiàn)5期顯著負效應，在低投喂量情景下2個池塘均在部分滯后期表現(xiàn)為顯著正效應，表明當投喂量較低時，水體的氨氮質量濃度相對于參考值有顯著上升。2號池塘水體濁度和電導率受投喂行為的影響相對明顯，高、低投喂量情景下分別有11期(32.4%)和15期(44.1%)使水體濁度顯著上升，分別有26期(66.7%)和14期(35.9%)使電導率顯著下降。2種情景下飼料投喂量對葉綠素質量濃度變化均無顯著效應。

表4 各池塘0-2觀測點水質在高、低投喂量情景下的滯后效應分析

由于在各池塘布設的12個觀測點均有對溶解氧質量濃度的監(jiān)測，選用模型結果更顯著的2號池塘分析不同觀測點溶解氧質量濃度在各投喂量下的滯后效應(圖4)。1-3和3-3 2個觀測點模型結果均不顯著，無統(tǒng)計學意義。其余10個觀測點中，由于飼料投喂對溶解氧質量濃度帶來的顯著正效應較少，負效應較多。當投喂量降低至9 kg/hm2時，7個觀測點部分滯后期的溶解氧質量濃度指標呈顯著負效應，且多出現(xiàn)在滯后期前中期(如0-1點)和后期(如3-1點)。其中，1-2號點位的10階滯后負效應最大，使得溶解氧質量濃度下降1.731 mg/L(95%置信區(qū)間：-2.647～-0.816 mg/L)。投喂量升高至181 kg/hm2時，9個觀測點部分滯后期的溶解氧質量濃度指標表現(xiàn)出顯著負效應，該負效應多出現(xiàn)在前期和后期(如0-1點)。其中，滯后期數(shù)最多的點位為0-1，在5～10期和30～38期共15期滯后，使溶解氧質量濃度顯著降低。在所有觀測點中，2-2點位的4階滯后期使溶解氧質量濃度下降的效應最大，下降了1.781 mg/L(95%置信區(qū)間：-2.916～-0.647 mg/L)。總體而言，同一深度不同位置及同一位置不同深度的溶解氧質量濃度滯后效應模式也有所不同，變化趨勢沒有明顯規(guī)律。

陰影和虛線區(qū)間部分為95%置信區(qū)間。圖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

3 結論與討論

本研究采用分布滯后非線性模型，在考慮滯后效應的基礎上，基于2020年宜興市中華絨螯蟹水產養(yǎng)殖試驗示范基地的2個試驗池塘數(shù)據(jù)，分析了生產環(huán)境下飼料投喂量對蟹塘水質的影響。研究結果可為生產環(huán)境下蟹塘水質的影響因素分析及水質調節(jié)策略的及時制定提供支撐。

研究發(fā)現(xiàn)，在高、低兩種投喂量下，飼料投喂與蟹塘溶解氧質量濃度、pH、氨氮質量濃度、電導率、濁度的變化之間存在長期滯后效應，但與葉綠素質量濃度變化之間的關系不顯著。在高投喂量情景下，相對于參考值，總體表現(xiàn)出氨氮質量濃度、電導率、pH值在部分滯后期顯著減少、濁度顯著上升的趨勢，其中溶解氧質量濃度在75%的觀測點顯著減少，滯后期內最大負效應達-1.781 mg/L(95%置信區(qū)間：-2.916～-0.647 mg/L)。已有研究者指出，中華絨螯蟹飼料投喂量較高時，飼料無法充分利用，殘餌、糞便等使得水體濁度顯著上升，隨著有機物質分解和中華絨螯蟹攝食行為的增加，耗氧量和二氧化碳的排放量有所增加，導致溶解氧質量濃度和pH值顯著降低[13,20]，這與本研究結論相一致。1號池塘的氨氮質量濃度和2號池塘的電導率在后期顯著減小，濁度在部分滯后期也表現(xiàn)出顯著的負效應，這可能是由于養(yǎng)殖戶在投喂量過大時會定期進行池塘改底等措施。試驗池塘的用藥記錄顯示，在中華絨螯蟹養(yǎng)殖過程中會使用特效芽孢桿菌、底居安、碧水安等藥物降解水體中的有機物、氨氮等有害物質，從而使池塘氨氮質量濃度在投喂后15～20 d、電導率在投喂后30～40 d、濁度在投喂后26～31 d左右顯著下降。在低投喂量情景下，相對于參考值，氨氮質量濃度和濁度在部分滯后期顯著增加，pH值顯著減小。這與餌料不足時，螃蟹為尋覓食物在池塘底部淤泥處大量活動有關[29]。另外，由于投喂量減少，殘餌量相對減少，因此2號池塘的電導率在多個滯后期呈顯著負效應。1號池塘的溶解氧質量濃度在部分滯后期顯著上升，但2號池塘在7個觀測點均出現(xiàn)對溶解氧質量濃度的顯著負效應, 滯后期內最大負效應達-1.731 mg/L(95%置信區(qū)間：-2.647～-0.816 mg/L)；pH值的負效應多出現(xiàn)在滯后的后期，在滯后期前期雖然估計效應值并不顯著，但相比于高投喂量情景效應值較小。這可能是由于投喂量減少后，前期溶解氧質量濃度上升、pH值增高，養(yǎng)殖戶基于較好的水質情況調整了投喂策略，增加了投喂量，使得后期又出現(xiàn)水質下降的趨勢。

研究發(fā)現(xiàn)，同一深度不同位置及同一位置不同深度的溶解氧質量濃度受投喂量影響的滯后效應模式有所不同，多數(shù)觀測點滯后期在30 d以上。這與溶解氧在水體不同位置質量濃度分布有所不同有關?，F(xiàn)有研究結果表明，養(yǎng)殖池塘溶解氧質量濃度受到多種生態(tài)因素影響，并在水體中呈現(xiàn)出三維分布特征[30-31]。具體表現(xiàn)為在池塘同一水平面上不同位置的溶解氧質量濃度有所不同，在同一垂直面上不同水層的溶解氧質量濃度也不同[4, 32]，這與本研究在不同觀測點觀測到的滯后效應模式有所不同的結果相符合。

本研究結果可輔助生產環(huán)境下蟹塘水質調節(jié)策略的制定。對于養(yǎng)殖戶，在投喂大量飼料后，宜著重關注投喂后5～10 d和20 d后的溶解氧質量濃度下降情況，以及投喂9 d后的pH值下降情況；在投喂少量飼料后，宜關注投喂后8～12 d和30 d后的溶解氧質量濃度下降情況、投喂后10～18 d pH值下降情況、投喂后10 d內氨氮質量濃度和投喂后5～20 d池塘濁度變化，及時通過開啟增氧泵設備、改變投喂策略、池塘改底等方式調節(jié)養(yǎng)殖水質。本研究也存在一定的局限性，通過DLNM模型的構建，可以給出投喂和水質參數(shù)之間暴露-反應關系的總結，但由于未報告各滯后期、投喂量對應影響效應的統(tǒng)計顯著性[28]，因此本研究采用了特定投喂量和滯后期進行分析。另外，由于中華絨螯蟹養(yǎng)殖水質受到氣象、生態(tài)、養(yǎng)殖經驗、投喂、用藥等多方面影響，且同一池塘不同位置的水質變化情況有所不同，未來還需要融合更多類型的數(shù)據(jù)并將其集成在分析過程中，進一步明確蟹塘養(yǎng)殖水質的復雜變化規(guī)律。