基于GAM 的黃海南部越冬小黃魚資源豐度與環境因子關系

李國東，李冬佳，熊 瑛，仲霞銘，宋大德，康中杰，湯建華，施金金，吳 磊

（1.江蘇省海洋水產研究所，江蘇南通 226007；2.上海海洋大學海洋科學學院，上海 201306）

小黃魚（Larimichthyspolyactis）隸屬鱸形目，石首魚科，黃魚屬，系暖溫性底層魚類［1］，廣泛分布于中國渤海、黃海、東海及朝鮮半島西岸海域［2］。小黃魚在我國漁業中占有重要地位，其資源曾經歷興盛期、衰退期和恢復期幾個階段［3-5］，海洋捕撈產量從20世紀80年代末期的2萬t已上升到2020年的近30萬t［6-7］。但與之相伴隨的是小黃魚低齡化、小型化和性成熟提前等特征日益凸顯［3，8］。小黃魚是一種洄游性魚類，一年四季在產卵場、索餌場和越冬場之間有規律地做季節性洄游［9］。劉效舜等［10］研究表明，黃海存在3個主要小黃魚種群：黃海北部-渤海種群、黃海中部種群以及黃海南部種群。春季黃海南部小黃魚種群從沙外漁場遷徙到大沙漁場中西部，秋季則開始向東遷移到外海越冬場，在冬季抵達沙外漁場西部水域進行越冬［11］。魚類越冬洄游是魚類種群動態的重要組成部分，與海洋環境因子間有著密切聯系。近年來，受氣候變化和人類捕撈活動的影響，小黃魚生長洄游習性發生了變化，越冬場位置有一定程度偏移［12-13］。因此，分析越冬小黃魚資源豐度對環境變化的響應，對于研究小黃魚種群補充和漁業資源管理具有重要意義。

核密度估計（kernel density estimation，KDE）是一種用于計算要素在其周圍鄰域中密度的非參數方法［14］。相較于其他空間分析方法，核密度估計法需要參數較少，受人為主觀因素影響較小，因此是探究漁業資源空間分布模式的有效工具［15］。廣義可加模型（generalized additivemodel，GAM）是一個能很好擬合高維數據中因變量與自變量之間非線性關系的生態位模型［16］，已被廣泛應用于漁業資源分布與環境關系相關研究［17-18］。然而，關于越冬期間黃海南部小黃魚種群資源分布與環境的關系尚少見研究報道。另外，目前近海漁業資源分布與環境因子關系的研究主要是利用調查站點的捕撈漁獲數據及在該站位使用CTD等儀器測定的環境數據，如海水溫度、海水鹽度、葉綠素a濃度等［19-21］。隨著遙感衛星技術的發展，海洋遙感環境數據精度不斷提高。與使用CTD等儀器測量站點數據相比，遙感環境數據具有覆蓋面廣、人工成本低和自動化獲取等諸多優勢，促進了漁場分析預報等應用研究的發展［22］。

本文利用2010—2019年冬季黃海南部小黃魚漁獲數據，基于核密度估計法分析越冬期間小黃魚資源分布情況，并結合相應時空的海洋衛星遙感監測數據，運用GAM探究黃海南部越冬小黃魚資源豐度與環境因子的關系，以期為預測小黃魚越冬漁場、制定小黃魚資源合理捕撈和科學養護政策提供科學依據。

1 材料與方法

1.1 數據來源

本文小黃魚捕撈數據來源于江蘇省海洋水產研究所2010—2019年冬季（12月—2月）在黃海南部進行的帆張網商業漁船捕撈日志，主要包含捕撈小黃魚所在漁區、小黃魚產量以及所布網次數據。依據小黃魚捕撈樣本數據所在漁區（定義為10′×10′的空間分辨率）進行了地圖可視化處理（圖1）。

圖1 2010—2019年捕撈越冬小黃魚分布位置Fig.1 Distribution locations of overw intering Larim ichthys polyactis fishing data during 2010—2019

依據小黃魚資源分布與環境因子關系相關研究［20-21］，初步選擇水深（water depth，Depth）、底層溫度（sea bottom temperature，SBT）、底層鹽度（sea bottom salinity，SBS）和葉綠素a濃度（chlorophyll-a concentration，Chla）作為環境因子。其中，水深數據來源于英國海洋學數據中心（https：／／www.gebco.net／data＿and＿products／gridded＿bathymetry＿data），空間分辨率為36″，數據更新時間為2017年12月；葉綠素a濃度數據來源于NOAAOceanWatch（https：／／oceanwatch.pifsc.noaa.gov／），空間分辨率為3′，時間分辨率為月；底層溫度和底層鹽度數據來源于中科院大氣物理所網站（http：／／www.ocean.iap.ac.cn），空間分辨率為30′，海水分層測定的垂直分辨率約為5～10 m，時間分辨率為月。

1.2 數據處理與分析

1.2.1 標準化處理

由于2010—2019年商業捕撈漁船數據及其攜帶的網具規格不一，需要對所有漁獲量數據進行標準化處理，剔除網具規格差異導致的誤差，因此采用單位捕撈努力量漁獲量（catch per unit effort，CPUE）表示小黃魚越冬資源豐度，計算方法如下：

式（1）中，C為小黃魚漁獲量（kg），E為該漁獲量所投入的捕撈努力量，S為所使用的帆張網網具迎流面積（m2），N為投網次數。

由于捕撈數據和環境因子數據的空間分辨率不完全一致，且它們之間呈獨立狀態，需利用ArcGIS10.2軟件進行處理。處理過程如下：通過捕撈漁獲位置生成相應捕撈點集，基于此點集對各環境柵格數據采用最鄰近分配法進行采樣，并通過空間位置信息連接捕撈點集的資源豐度和各環境要素數據，最后進行點集合并及統計分析。

1.2.2 核密度估計法分析

本文基于計算得到的小黃魚產量數據，運用核函數［14］進行分析，能直觀反映出黃海南部越冬小黃魚相對資源豐度的時空分布情況。

式（2）中，K［］是核函數；h為帶寬；n是在帶寬范圍內的已知點數目；d是數據的維度。

1.2.3 影響因子篩選

本文選取小黃魚越冬資源豐度作為因變量。除了初選的環境因子外，亦加入時間（年份）作為初始解釋變量，探究黃海南部越冬小黃魚資源豐度的年際變動情況。本研究采用條件數k［23］和方差膨脹因子（variance inflation factor，VIF）［24］對初始因子進行多重共線性檢驗，從而剔除存在共線性的因子。

設x1，x2，x3，…，xp是選取的初始因子，X1，X2，X3，…，Xp是經過中心化和標準化的向量，記X=［x1，x2，x3，…，xp］，設λ為XTX的一個特征值，φ為對應的特征向量，其長度為1，即φTφ=1。度量多重共線性嚴重程度的一個重要指標是矩陣XTX的條件數k：

式中，λmax（XTX）、λmin（XTX）分別表示矩陣的最大、最小特征值。條件數表征了XTX特征值差異的大小。若條件數k＜100，則表示初始因子之間存在較小的多重共線性程度；若100≤k≤1 000，則表示多重共線性為中等或略強程度；若k＞1 000，則表示多重共線性十分嚴重［23］。條件數k的計算可通過R語言kappa函數實現。

此外，VIF的平方根表示變量回歸參數的置信區間能膨脹與模型無關的預測變量的程度，可通過使用R語言中car包中的VIF函數確定VIF值。一般認為VIF＞4，即存在多重共線問題［24］。

1.2.4 GAM模型的構建

本研究以小黃魚資源豐度為響應變量，以初步篩選的相關因子作為解釋變量，利用GAM模型對小黃魚越冬資源豐度與相關環境因子的關系進行擬合。GAM模型的一般形式如下［25］：

式（4）中，Y為小黃魚資源豐度，即每1 m2網具迎流面積的漁獲量，單位為g·m-2；α是適合函數的截距；ε是服從正態分布的隨機誤差項；xj表示解釋變量，即各站位的影響因子；fi（xj）是各環境變量xj的任意單變量函數，通過樣條平滑函數來估計。模型的誤差分布估計為Gaussian分布，連接函數為自然對數。

依照赤池信息準則（Akaike information criterion，AIC），在AIC值最小的因子預測函數上依次加入其他因子，進而得到AIC值最小的多因子預測模型［25］。其中，AIC值越低，模型的擬合效果越好；方差解釋率越高，模型擬合效果越好。AIC的計算過程如下：

式中，k是參數個數，L是似然函數。

2 結果與分析

2.1 越冬期間小黃魚資源分布

黃海南部越冬期間小黃魚資源分布呈現一定的聚集性，主要分布在32°24′～34°00′N、123°00′～125°00′E區域。而各月份小黃魚資源分布情況有所不同（圖2～圖4），其中12月小黃魚資源核密度高值區域主要集中在32°48′～33°48′N、123°00′～124°10′E，即黃海南部水深45～60 m海域。1月小黃魚資源逐漸遷移至水深55～75 m海域，集中在32°24′～34°00′N、123°24′～125°00′E區域，呈現東南-西北分布。2月小黃魚資源分布區域與1月份基本相同，但資源核密度高值區域主要集中在大沙漁場中北部，且該區域的越冬小黃魚資源核密度與1月相比明顯增加。

圖2 12月黃海南部越冬小黃魚資源核密度分布Fig.2 Kernel density distribution of Larim ichthys polyactis resources in December

圖3 1月黃海南部越冬小黃魚資源核密度分布Fig.3 Kernel density distribution of Larim ichthys polyactis resources in January

圖4 2月黃海南部越冬小黃魚資源核密度分布Fig.4 Kernel density distribution of Larim ichthys polyactis resources in February

2.2 小黃魚越冬洄游影響因子

通過條件數k和VIF對初始因子進行多重共線性檢驗（表1），結果顯示初始因子的條件數小于100，VIF值均小于4，即初始因子的多重共線性程度均較小，因此構建GAM時將年份、底層溫度、底層鹽度、葉綠素a濃度以及水深作為本研究的解釋變量。

表1 初始因子的多重共線性檢驗Tab.1 Multicollinearity test of influencing factors

通過GAM模型擬合小黃魚越冬資源豐度與各影響因子間的關系，利用AIC準則進行因子篩選（表2）。根據AIC值最小原則得到最優GAM模型，模型的最終表達式為：

表2 最優GAM 模型中各影響因子的相關參數Tab.2 Relevant parameters of each factor in the optimal GAM model

式（6）中，Y為小黃魚越冬資源豐度，α為截距，s（Year）為年份因子效應，s（SBT）為底層溫度效應，s（SBS）為底層鹽度效應，s（Chla）為葉綠素a濃度效應，s（Depth）為水深效應。

小黃魚最優GAM模型及各影響因子相關參數表明，影響黃海南部小黃魚資源越冬洄游的變量為水深、底層溫度、底層鹽度、葉綠素a濃度和年份。模型累積偏差解釋率為54.4%。其中在各海洋環境因子中，對小黃魚越冬洄游影響最大的因子為水深（34.3%），最小的為葉綠素a濃度（4.6%）。

2.3 小黃魚越冬對環境因子的響應

本文建立GAM模型時加入年份因子作為初始解釋變量，結果顯示2010—2017年黃海南部越冬小黃魚資源豐度的年際變動較小，而2017年后小黃魚資源豐度有所上升（圖5）。在選取的影響因子中，水深對小黃魚的資源豐度影響明顯最大（表2）。隨著水深的增加，小黃魚資源豐度呈現先減少后增加的趨勢。其中，水深在30～40 m范圍內，小黃魚資源豐度隨水深增加而減少；在40～70 m水深范圍內，小黃魚相對資源豐度總體呈現隨著水深增加而明顯增加趨勢。底層溫度對小黃魚越冬資源豐度影響波動較大，但整體上資源豐度隨著溫度升高而振蕩上升。底層鹽度范圍在32.8～33.5時，小黃魚資源豐度隨著鹽度的增加而明顯增加；當鹽度大于33.5時，小黃魚資源豐度一直較為穩定。葉綠素a濃度范圍在1.0～2.3 mg·m-3時，小黃魚資源豐度隨著葉綠素a濃度增加而明顯增加；當葉綠素a濃度大于2.3 mg·m-3時，小黃魚資源豐度隨著葉綠素a濃度增加而緩慢減少。

圖5 基于GAM 模型的小黃魚資源豐度對環境因子響應情況Fig.5 Response of Larim ichthys polyactis abundance to environmental factors based on GAM model

3 討論

小黃魚越冬洄游期間資源豐度高值區域一定程度上可看作為小黃魚的越冬場。越冬場對于補充小黃魚種群十分關鍵，因為越冬場匯集了一年中的大部分小黃魚成魚［4］。本研究表明，2010—2019年黃海南部越冬小黃魚群體資源豐度高值區域在空間分布上均呈現一定的聚集性（圖2～圖4），主要分布在32°20′～34°00′N、123°00′～125°00′E區域，與之前研究結果基本相同［10-13］。此外，越冬期間各月份小黃魚資源分布情況有所不同，其中12月小黃魚資源核密度高值主要集中在32°48′～33°48′N、123°00′～124°10′E海域。隨著海溫的逐漸降低，1月和2月黃海南部小黃魚資源繼續向外海遷移，即集中在32°20′～34°00′N、123°24′～125°00′E區域，呈現東南-西北向分布，與冬季黃海暖流在黃海南部的主軸位置和流向一致［26］。冬季黃海南部海域受黃海暖流的影響明顯，黃海暖流是對馬暖流和陸架水混合而成，它主要是在濟州島西側海域鋒區中衍生出來的具有高溫高鹽特征的海流［26］。因此該海域為黃海南部小黃魚群體越冬洄游提供了合適的場所。本研究為探尋黃海南部小黃魚越冬群體提供了重要的空間信息，據此可以通過保護一部分小黃魚成魚以維持小黃魚的合理種群補充水平。

海洋水體環境是海洋魚類賴以生存的空間，因此可通過對海洋水體環境要素的信息獲取、越冬小黃魚資源的時空分布分析研究小黃魚越冬的環境偏好，進而可為小黃魚越冬漁場海況預報提供重要支撐。GAM建模結果顯示，黃海南部小黃魚越冬洄游主要受水深和底層溫度影響。隨著水深的增加，小黃魚資源豐度呈現先減少后增加的趨勢（圖5）。其中水深在30～40 m范圍內，小黃魚資源豐度隨水深增加而減少；在40～70 m水深范圍內，小黃魚相對資源豐度總體呈現隨著水深增加而明顯增加趨勢。這與陳新軍［27］研究結果基本一致，即越冬期的小黃魚在黃海區域主要分布在55～75 m。海水鹽度的顯著變化是支配魚類行為的一個重要因素，對魚類滲透壓甚至整個魚類生活史產生重要影響［28］。鹽度與魚類行動的關系主要表現在間接方面，通過水團、海流等起到了影響作用，如暖水性魚類隨著暖流的高溫高鹽進行洄游［27］。本研究擬合結果顯示，底層鹽度范圍在32.8～33.5時，小黃魚資源豐度隨著鹽度的增加而明顯增加；當鹽度大于33.5時，小黃魚資源豐度一直較為穩定（圖5）。海水鹽度變化不是孤立的物理變化現象，它會隨著水系、水流等運動而變化，這對于探測小黃魚越冬分布具有一定意義，后期可通過了解鹽度分布及其與海流運動的關系來推測小黃魚漁場的可能位置。海洋生物受海水溫度變化刺激所產生的行為是主動選擇最適的溫度環境，以使其體溫維持在一定的溫度范圍內。因此海水溫度直接影響和控制著海洋生物種群分布及其洄游和繁殖過程，是影響魚類生態習性的主要環境因子之一［29］。小黃魚屬于暖溫性洄游魚類，其生活習性受溫度等水文環境的影響顯著［19］，越冬洄游主要是由于環境溫度降低所引起的，因此魚群洄游到漁場時間、資源豐度以及漁期長短均與漁場水溫有密切關系。近10年數據擬合結果表明，底層溫度對小黃魚越冬資源豐度影響波動情況較大，但整體上資源豐度隨著溫度增加而振蕩增加（圖5），這是因為屬暖溫性的小黃魚越冬洄游傾向于水溫較高的海域。海洋遙感反演葉綠素a濃度在一定程度上代表該海域的初級生產力水平的高低，常用于估算海洋浮游生物等魚類餌料的時空分布，因此葉綠素a濃度對小黃魚資源分布也會產生重要影響。研究結果顯示，葉綠素a濃度范圍在1.0～2.3 mg·m-3時，小黃魚資源豐度隨著葉綠素a濃度增加而明顯增加；當葉綠素a濃度大于2.3 mg·m-3時，小黃魚資源豐度隨著葉綠素a濃度增加而緩慢減少，表明黃海南部小黃魚越冬洄游資源豐度與葉綠素a濃度并非是完全正相關關系。這可能與在越冬期間魚類通常減少或停止攝食、而主要依靠索餌期間體內所積累的營養來供應機體能量消耗有關［27］。所以在越冬洄游時期，葉綠素a及餌料的分布和變動在一定程度上難以支配小黃魚越冬洄游行為。本文主要利用葉綠素a濃度數據進行GAM建模，葉綠素a濃度與其他環境因子相比，其解釋率僅為4.6%，即對黃海南部小黃魚越冬洄游產生影響較小。但后期可通過對比分析索餌期與越冬期海水葉綠素濃度時空動態特征以及此背景下小黃魚資源豐度情況，深入探究海水葉綠素濃度與小黃魚索餌和越冬洄游資源豐度之間的時空關系，為預測小黃魚秋季索餌漁場、越冬漁場位置和制定相關漁業管理政策提供重要參考。

本文建立GAM模型時加入時間（年份）因子作為初始解釋變量，可為探究黃海南部越冬小黃魚資源豐度的年際變動情況、分析小黃魚越冬適宜生境動態變化提供科學依據。氣候環境變化和人類捕魚活動差異可造成魚類資源的時空分布變化［30］，并改變其生活史特征和遷移模式［31-32］。探究海洋環境變化與魚類行為的關系對深入了解和掌握海洋魚類生活史過程、開發和利用魚類資源極為重要。所以后期可基于小黃魚洄游群體資源豐度與環境因子關系，結合海洋衛星環境監測數據對小黃魚越冬洄游海域和時間、漁期長短等進行預測。