毛里塔尼亞海域短線竹?魚時空分布與海洋環境的關系

宋利明，徐雙泉，許 回，陳明銳，Ebango Ngando Narcisse

（1.上海海洋大學海洋科學學院，上海 201306；2.國家遠洋漁業工程技術研究中心，上海 201306；3.大洋漁業資源可持續開發省部共建教育部重點實驗室，上海 201306；4.遠洋漁業協同創新中心，上海 201306）

短線竹?魚（Trachurus trecae）為中底層魚類，一般活動在20～100m的近海底水域，有時也在中上層甚至近表層出現，適宜水溫為15～22℃，廣泛分布在東大西洋的摩洛哥至安哥拉海域［1］。毛里塔尼亞西鄰大西洋，處在幾內亞灣暖流和加納利寒流的交匯地帶，海域初級生產力比較高，擁有豐富的漁業資源［2］。FAO的統計結果顯示［3-4］，2007年毛里塔尼亞海域主要的小型中上層魚類中，短線竹?魚占總漁獲量的比重最多，為31%；2013年中東大西洋的短線竹?魚產量主要集中在毛里塔尼亞沿海，約占82%。因此對毛里塔尼亞海域短線竹?魚的研究十分重要。目前國內學者如牛明香等［5-6］研究了智利竹?魚（Trachurusmurphyi）的漁場與時空變動的關系，李顯森等［7］、化成君等［8］和汪金濤等［9］分別對智利竹?魚的資源分布特征、時空分布與環境因子的關系、漁情預報技術作了研究。RUMOLO等［10］通過同位素和胃含物分析，研究了西西里海峽大西洋竹?魚（Trachurus trachurus）的空間分布與攝食行為的關系，GUSHCHIN和CORTEN［11］研究了毛里塔尼亞海域大西洋竹?魚和短線竹?魚的攝食特性，ARKHIPOV［12］研究了毛里塔尼亞海域魚類浮游生物的季節和年際變化，班銳［13］利用廣義可加模型探討了時空和環境因子對毛里塔尼亞海域竹?魚漁場的影響，但毛里塔尼亞海域的竹?魚包括了短線竹?魚在內的多種竹?魚類，其具體的生活方式和時空分布特征存在差異［1］，且其研究中使用的環境數據來源于衛星遙感，分辨率較低。分位數回歸模型最早由KOENKER和BASSETT［14］提出，已在研究魚類的棲息地綜合指數方面得到了應用［15-16］。本文根據2017年在毛里塔尼亞海域測定的海洋環境數據以及短線竹?魚的漁獲率統計資料，用分位數回歸的方法探究短線竹?魚時空分布與海洋環境的關系，旨在進一步了解短線竹?魚在毛里塔尼亞海域的時空分布特征，為提高我國漁船在此海域的作業效率及毛里塔尼亞對短線竹?魚等小型中上層漁業資源養護與管理提供參考。

1 材料與方法

1.1 調查方法

1.1.1 調查船及基本信息

“福遠漁097”燈光圍網船總長49.95 m，型寬8.60 m，型深4.00 m。設計吃水3.30 m，設計排水量896.90 t，主機功率為828 kW，航速為6～11 kn，魚艙容量190m3，油艙容積120m3，淡水艙容積47m3。航區為國際I區。

1.1.2 調查時間、海域及站點分布

調查時間為2017年9月20日至12月31日。調查海域為：16°28′～17°59′W、17°50′～20°58′N。共119個調查站點（部分站點位置非常接近，圖中重疊在一起），隨機選擇75%的（90個）站點的數據用作建模，25%的（29個）站點數據用以驗證模型，站點分布見圖1。部分建模站點和驗證站點出現重合，原因是這些站點的下網位置幾乎相同。

1.1.3 調查的漁具漁法

“福遠漁097”所用圍網網具主體由網衣、包鉛繩、浮子、上綱和下綱組成，網口網衣拉直周長約1 643.4 m，網口至囊端拉直長度約116 m，上綱長度約809.84 m，下綱長度約761.44 m。

1.1.4 調查儀器、方法及內容

利用漁船上的全球定位儀記錄每天放網開始時的經緯度，同時記錄漁獲重量。調查站點預先確定，但實際調查位置與計劃站點存在一定偏差。每日放網后或起網前，使用加拿大型號為XR-620的多功能水質儀，測量調查站點海表面的葉綠素a濃度、海表面溫度以及鹽度數據。

1.2 數據處理方法

1.2.1 CPUE的計算方法

單位捕撈努力量漁獲量（catch per unit effort，簡稱CPUE）是用來表示魚類資源豐度的相對指數，是對某個漁業資源進行評估的基本內容之一［17］。船上所得到的直接數據為漁獲量（t），建立關系模型前，需要計算短線竹?魚的CPUE值，計算方法如下：

式（1）中，ci為在i站點的漁獲量（t），ni為i站點的網次數量。

1.2.2 海洋環境因子值的計算方法

海洋環境因子數據為各個站點10 m以淺所測得的各海洋環境因子數據的算術平均值［15］。

1.3 建立短線竹?魚CPUE預測模型

分位數回歸建模的具體方法可以參照文獻［15］。分位數回歸分析采用由美國地理調查局中陸生態科學研究中心開發的統計學軟件Blossom。

根據90個建模站點的CPUE、葉綠素a濃度（Ch）、海表面溫度（T）、海表面鹽度（S）及其交互作用項建立短線竹?魚CPUE與海洋環境因子的關系模型，其一般形式為：

圖1 調查站點分布圖Fig.1 Distribution map ofmodeling sites and validation sites

式（2）中，C為常數項，a-f為自變量的系數，ChT表示葉綠素a濃度和海表面溫度的交互作用項，ChS表示葉綠素a濃度和海表面鹽度的交互作用項，TS表示海表面溫度和海表面鹽度的交互作用項，ε為誤差項。i表示站點號數，如Chi表示在i站點的葉綠素a濃度。

1.4 棲息地綜合指數（IHI）計算方法

根據關系模型，得出建模站點和驗證站點的預測CPUE，按照下列公式計算得出短線竹?魚在不同站點的IHI值：

式（3）中，CPUEij是指j組（j=1為建模站點；j=2為驗證站點）i站點的CPUE預測值，CPUEmaxj是指j組中CPUE預測值中的最大值。

1.5 IHI模型預測能力的驗證和評價

將建模站點和驗證站點的環境數據值分別代入模型中，使用SPSS 23.0軟件，通過Wilcoxon（符號秩）對建模站點和驗證站點的CPUE實際值與預測值進行檢驗［18］，根據P值的大小判斷兩者間是否存在顯著性差異，以此評估模型的預測能力。通過計算短線竹?魚驗證站點的CPUE實際值與預測值之間的Spearman相關系數（秩相關系數）大小，結合雙尾檢驗對IHI模型的預測能力進行評價。

用Arcmap10.6作圖得出建模站點和驗證站點的IHI等值線與實測CPUE的分布圖。將29個驗證站點的實測CPUE值除以其中的最大值得出IHI′［式（4）］，用Excel作出IHI′和IHI的變化趨勢圖，定性評價模型的預測能力。

式（4）中，CPUEi2是指驗證站點i的實測CPUE，CPUEmax2是指驗證站點實測CPUE中的最大值。

2 結果與分析

2.1 分位數回歸得出的預測模型

利用Blossom建立CPUE與海洋環境的關系模型，將每次擬合得出的模型系數和對應的P值記錄下來，取分位數（用Q表示）最高的模型為最佳模型。結果如表1所示。

擬合得出的最佳預測模型如下：

由式（5）可知，影響短線竹?魚CPUE的海洋環境因素包括海表面溫度、葉綠素a濃度、海表面溫度與葉綠素a濃度的交互項，P值依次為0.028、0.004和0.032。

由自變量系數的大小可得，對CPUE影響最大的是海表面溫度，其次是葉綠素a濃度，再次是海表面溫度和葉綠素a濃度的交互作用項。

2.2 IHI模型驗證

對建模站點和驗證站點的CPUE實測值與預測值進行Wilcoxon（符號秩）檢驗，結果見表2，兩者的P值均高于顯著水平0.05，說明建模站點和驗證站點的CPUE實測值與預測值之間均無顯著性差異。

對29個驗證站點的實測CPUE和預測CPUE的相關系數和顯著性進行雙尾檢驗，得出實測CPUE與預測CPUE的Spearman相關系數達到0.832，相伴概率小于0.01，這說明顯著性水平為0.01時，實測CPUE和預測CPUE之間高度相關，模型可靠。

表1 不同分位數下各變量回歸系數秩得分檢驗結果Tab.1 Rank score test for regression parameters under different quantiles

表2 建模站點和驗證站點的實測CPUE與預測CPUE間W ilcoxon檢驗結果Tab.2 Result of W ilcoxon test between the observed and predicted CPUE ofmodeling sites and validation sites

2.3 IHI模型的定性評價

用Arcmap10.6作圖得出的建模站點的IHI等值線與實測CPUE分布見圖2。驗證站點的IHI等值線與實測CPUE分布見圖3。

本次調查期間，短線竹?魚的CPUE在0～3 t·次-1，因此在作圖和分析時，認為CPUE小于1 t·次-1時為CPUE的較低值，CPUE大于2 t·次-1時為CPUE的較高值。由圖2可觀察到，實測CPUE與IHI具有一定的相關性。IHI值在小于0.40時，實測CPUE值小于1 t·次-1；IHI值在0.40～0.56之間時，CPUE值隨IHI值的增大而相應增大，大部分在1～2 t·次-1之間，大于2 t·次-1的情況也有出現，但分布零散，具有偶然性；IHI值大于0.60時，CPUE較高值出現的頻率明顯增加，大于2 t·次-1的情況集中出現。

觀察圖2可得，CPUE較低值集中分布在16°15′～16°50′W、16°50′～18°50′N，對應的IHI值低于0.56；CPUE較高值集中分布在17°20′～17°45′W、20°15′～20°50′N，對應的IHI值大于0.60。

圖2 建模站點IHI等值線和實測CPUE分布Fig.2 IH I contour and CPUE distribution for modeling sites

圖3中，驗證站點的CPUE與IHI呈現出的相關性較為顯著。當IHI值低于0.40時，對應的海域為16°15′～17°40′W、18°05′～20°05′N，CPUE較低值在這一區域零散分布，大于1 t·次-1的只有1個，其余都在1 t·次-1以下。當IHI值大于0.50時，對應的海域為17°20′～17°45′W、20°15′～20°50′N，CPUE較高值集中出現，CPUE全部大于1 t·次-1，大于2 t·次-1的驗證站點也都在上述區域出現。這說明IHI對CPUE的預測具有良好的效果。

29個驗證站點的IHI′和IHI的變化趨勢見圖4。可以觀察到，IHI的變化趨勢與IHI′高度一致，基本呈現正相關，側面反映了預測CPUE的變化趨勢與實測CPUE十分吻合。這可說明，建立的模型能較好地預測CPUE。

圖3 驗證站點IHI等值線和實測CPUE分布Fig.3 IHI contour and CPUE distribution for validation sites

圖4 驗證站點IH I′和IH I的變化趨勢Fig.4 Trends of IHI′and IHI for validation sites

綜上，IHI值較高的海域對應的實測CPUE較高，因此基于分位數回歸擬合得出的關系模型能較好地反應CPUE與海洋環境因子的關系，預測能力良好。

3 討論

3.1 用分位數回歸模型及相關統計檢驗方法的依據

分位數回歸可提供不同分位數下的回歸方程，反映了自變量對因變量的某個特定分位數的邊際分布，因此能清楚地闡釋因變量的整個分配，甚至數據異質性的問題也能解決［19］。本文所使用的樣本數據時間跨度只有3個多月，樣本數據總量也只有119個，加上部分站點數據用作驗證模型，實際上用于建模的數據并不多，因此樣本數據并非呈正態分布，而是呈正偏態分布。當研究的數據集合呈現非正態分布時，選擇分位數回歸能夠具體地反映因變量的平均水平［20］。利用SPSS對樣本數據進行非參數檢驗，當樣本數據呈現正態分布時，可用t檢驗，非正態分布或無法確定分布狀態時，用Wilcoxon秩和檢驗［21］；進行相關性分析時，樣本數據呈正態分布用Pearson檢驗，非正態分布使用Spearman檢驗［22］；未知呈正相關和負相關，使用雙尾檢驗。

3.2 海洋環境因子的影響與時空分布

由基于分位數回歸得出的最佳關系模型可以看出，海洋環境因子對短線竹?魚CPUE的影響程度為：海表面溫度＞葉綠素a濃度＞海表面溫度與葉綠素a濃度的交互作用，海表面鹽度對CPUE的影響沒有在關系式中被表征出來。

漁場的形成與分布是由各種海洋環境因子的綜合作用共同決定的。在所有非生物環境因子對魚類資源分布的影響中，牛明香等［6］指出，水溫是對漁場形成最為重要的因素，海表面溫度的變化能夠影響竹?魚的分布和洄游等。本文的結論同樣為海表面溫度對短線竹?魚的CPUE影響最大，這也符合實際情況：毛里塔尼亞海域處在幾內亞灣暖流和加納利寒流的交匯地帶，有利于漁場形成，寒暖水團對應著溫度的變化，因此海表面溫度對短線竹?魚的影響尤其顯著。浮游生物是短線竹?魚等中上層魚類的食物來源，也是漁場形成的決定性因素，毛里塔尼亞海域為幾內亞灣暖流和加納利寒流冷暖水團的交匯處，攜帶大量的營養鹽，海域初級生產力較高，為浮游生物的繁殖創造了有利條件，因此葉綠素a濃度是形成漁場的另一個重要的海洋環境要素。需要指出的是，雖然關系模型顯示葉綠素a濃度是影響短線竹?魚CPUE的一個因素，但由于關系模型中其系數為-0.027 3，因此葉綠素a濃度對短線竹?魚CPUE的影響相對海表面溫度影響較小，這是因為葉綠素a雖然對于漁場的形成不可或缺，但形成漁場的關鍵還要取決于浮游生物的量和規模［23］，而除了葉綠素a濃度外，營養鹽、海流等其他海洋要素同樣影響著浮游生物的量和規模。牛明香等［5］利用廣義可加模型研究時空和環境因子對東南太平洋智利竹?魚漁場的影響時得出相似的結論，認為葉綠素a濃度雖然與資源密度有關聯，但對整個模型的貢獻程度相對較低，可作為輔助指標。RUMOLO等［10］通過同位素和胃含物分析西西里海峽大西洋竹?魚的空間分布與攝食行為的關系時指出，大西洋竹?魚的空間分布與溫度、光合作用的強度有關，此外與鹽度也有關系。本文的關系模型中，海表面鹽度對短線竹?魚CPUE的影響沒有被表征出來，可能是因為調查期間，作業地點的經緯度變化并不大，海表面鹽度值上下浮動很小，鹽度梯度相較溫度梯度低。

本文通過分位數回歸建模，得出關系式后建立IHI模型，并結合Arcmap作圖（圖2，圖3）分析，得出毛里塔尼亞短線竹?魚CPUE較高的區域為：17°20′～17°45′W、20°15′～20°50′N；CPUE較低的區域為16°15′～17°40′W、16°50′～18°50′N。這在資源豐度的時空分布上反映為9—12月短線竹?魚逐漸由南向北轉移，ARKHIPOV［12］在研究中指出，毛里塔尼亞海域的短線竹?魚等魚類在寒冷季節會向北聚集，在布朗角附近產卵，距離南部有一定距離，這從側面證明了本文研究結果的準確性。張鵬等［24］利用FAO報告和相關資料研究中東大西洋中上層小型魚類資源及其漁業現狀時總結得出短線竹?魚最密集分布區在21°00′N附近，這與本文的研究結果在緯度上大體一致，但也存在差異，原因在于研究所使用的樣本數據在經緯度上的分布和分辨率不同。

3.3 建議

本文認為IHI模型能較準確地預測短線竹?魚的CPUE。短線竹?魚IHI分布較高和資源豐度較高的區域在17°20′～17°45′W、20°15′～20°50′N，因此建議我國燈光圍網船在毛里塔尼亞海域9—12月作業時，盡量使下網區域集中在上述區域，以增加產量，提高經濟效益。短線竹?魚不僅是毛里塔尼亞海域小型中上層漁業的重要組成部分，也是金槍魚類等的主要食物，因此對其養護與管理十分重要。IHI分布較低的區域是16°15′～17°40′W、16°50′～18°50′N，由于短線竹?魚主要攝食浮游生物，因此在9—12月，可以投放少量流木為浮游生物的繁殖創造條件。

3.4 不足與展望

毛里塔尼亞海域的竹?魚主要包括大西洋竹?魚和短線竹?魚，它們有相似的生物學習性，但具體的棲息水層等存在差異，因此在研究竹?魚的時空分布時，深度也是應當考慮的變量之一。本文缺乏短線竹?魚棲息深度的數據，因此空間上未能說明短線竹?魚的垂直分布。此外，西北非海域上升流顯著［25］，魚類的分布還受溶解氧、營養鹽的分布等因素影響［26］，因此要更全面地探索短線竹?魚時空分布與海洋環境的關系，還需要考慮上述影響因子。本次調查時間跨度不大，有效數據偏少，因此之后的研究要盡可能增加觀測時長和站點數量。另外，如有可能，可以采集各水層的浮游生物量和短線竹?魚的胃含物等，為更深入地研究短線竹?魚的運動模式提供支撐。

致謝：本研究得到了毛里塔尼亞政府的許可，并得到宏東國際（毛里塔尼亞）漁業發展有限公司蘭平勇、林蘭英、曾志堅、趙松輝、朱國平先生和“福遠漁097”號船長及全體船員等大力支持，謹致謝意。