基于長度貝葉斯生物量估算法評估北部灣大頭白姑魚資源狀況

王淼娣 ，王雪輝，孫典榮，王躍中，陳新軍，杜飛雁，邱永松

（1. 上海海洋大學海洋科學學院，上海 201306； 2. 中國水產(chǎn)科學研究院南海水產(chǎn)研究所/廣東省漁業(yè)生態(tài)環(huán)境重點實驗室/農(nóng)業(yè)農(nóng)村部外海漁業(yè)開發(fā)重點實驗室，廣東 廣州 510300）

對漁業(yè)資源狀況進行科學評估，是有效管理漁業(yè)、緩解資源壓力的重要手段之一。傳統(tǒng)的漁業(yè)資源評估方法需要較多的參數(shù)輸入[1]，如剩余產(chǎn)量模型[2]、年齡結構模型[3]、延遲差分模型[4]等。輸入?yún)?shù)主要包括捕撈努力量、資源豐度指數(shù)、漁獲量、生長和年齡等。而上述參數(shù)的獲取，需投入較大的物力、經(jīng)費和專業(yè)技術人員，導致這些模型的運用和推廣受限。

長度頻率法是漁業(yè)資源數(shù)據(jù)缺乏下有效的評估手段之一，其優(yōu)點是長度組成數(shù)據(jù)易獲取[5]。近十幾年，基于長度的產(chǎn)卵潛力比評估[6]、基于單峰長度頻率數(shù)據(jù)[7]與多元回歸樹數(shù)據(jù)挖掘技術[8]等被廣泛應用到漁業(yè)種群評估中。隨著計算機技術的發(fā)展，F(xiàn)roese等[9]提出的一種新的基于長度的貝葉斯生物量估算法 (Length-based Bayesian biomass estimation method, LBB) 被應用。該方法可分析漁業(yè)種群參數(shù)和資源狀況，只需要輸入魚類的長度頻率數(shù)據(jù)，通過模型模擬即可獲得魚類的最適開捕規(guī)格、評估漁業(yè)資源的利用狀況 (如最大可持續(xù)產(chǎn)量的相對生物量與原始資源量的相對生物量的比值)。與其他漁業(yè)資源評估方法相比，LBB的主要優(yōu)點在于降低了初步資源評估的輸入數(shù)據(jù)要求[10]。

北部灣位于中國南海西北部，地處熱帶和亞熱帶，物種多樣、資源豐富，盛產(chǎn)魚、蝦、蟹及貝類，曾是我國的優(yōu)良漁場之一，也是中國和越南漁民的傳統(tǒng)作業(yè)場所[11]。為評估北部灣的漁業(yè)資源狀況，眾多學者對北部灣的漁業(yè)資源和主要經(jīng)濟種類的生物學做了較多研究。研究表明，北部灣漁業(yè)資源密度持續(xù)下降，魚類優(yōu)勢種更替明顯，壽命長、個體大的魚類種類減少[12]。進一步分析發(fā)現(xiàn)，捕撈強度的逐年增加是導致北部灣漁業(yè)資源衰退的主要原因[13]。此外，漁民非法捕撈行為和東南亞周邊國家的漁業(yè)政策，也是導致該海域漁業(yè)資源衰退的重要原因[11]。大頭白姑魚 (Pennahia macrocephalus)隸屬于鱸形目、石首魚科、白姑魚屬，是北部灣主要漁獲優(yōu)勢種之一，在南海北部魚類資源中具有重要地位，其產(chǎn)量約占白姑魚屬的30%[14]。

國內(nèi)外對大頭白姑魚的研究相對較少。在形態(tài)學方面，蘇新紅[15]、Yamada和Yamada[16]、Sasaki[17]對大頭白姑魚的形態(tài)進行了詳細描述分析，為其種屬的劃分提供了參考依據(jù)；在食物網(wǎng)方面，張其永等[18]通過對經(jīng)濟魚類的食物網(wǎng)研究發(fā)現(xiàn)大頭白姑魚的食物以底棲生物為主；在資源評估方面，多集中在白姑魚屬的其他種，如陳作志等[19]估算南海北部白姑魚 (P. argentata) 的生長和死亡參數(shù)，顏云榕等[20]和何雄波等[21]分別對斑鰭白姑魚 (P.pawak) 的生長、群體結構與分布進行分析。歷經(jīng)數(shù)十年的開發(fā)利用，作為北部灣主要優(yōu)勢漁獲的大頭白姑魚資源狀況目前尚不清楚。為此，本研究根據(jù)2006—2014年在北部灣開展的18個航次的大頭白姑魚生物學測量數(shù)據(jù)，利用LBB估算北部灣大頭白姑魚的漸進體長 (L∞)、相對自然死亡率(M/k)、相對捕撈死亡率 (F/k)、最適開捕體長(Lc_opt) 和相對生物量開發(fā)指標 (B/B0)，對大頭白姑魚的資源進行評估，以期為該魚種資源的持續(xù)利用與科學管理提供參考依據(jù)。

1 材料與方法

1.1 數(shù)據(jù)

本研究所用數(shù)據(jù)來自北部灣海域大頭白姑魚的生物學數(shù)據(jù)。采樣時間為2006—2014年，每年2個航次，共18個航次 (表1)。為便于表述，把1—2月作為冬季，7—8月作為夏季。每個航次布設調查站位52個，調查船為單拖漁船“北漁60011”，總噸位242 t，長度36.8 m，寬度6.8 m，主機功率為441 kW。采樣網(wǎng)具為404型底拖網(wǎng)，上綱長度為37.7 m，網(wǎng)口網(wǎng)目尺寸20 cm，網(wǎng)衣全長60.5 m，網(wǎng)囊目尺寸為3.9 cm[12]。

表1 采樣航次、時間與樣本量信息Table 1 Information of sampling cruise, sampling month and sampling size

1.2 種群參數(shù)估算

魚類體長數(shù)據(jù)以10 mm間距進行分組。采用LBB分析漁獲物的長度頻率數(shù)據(jù)，所有相關參數(shù)均采用馬爾科夫鏈-蒙特卡羅法，詳細推導過程見Froese等[9]，主要公式如下：

假設魚類生長遵循Von Bertalanffy生長方程[22]，

其中 Lt為t齡時的長度；k為達到 L∞的 速率； t0為魚的體長為0時的年齡。

其中 NL是長度為L的存活數(shù)量； NLstart是長度為Lstart的數(shù)量，具有完整的選擇范圍，且所有進入漁具的個體都被漁具保留；Z/k是總死亡率與生長率的比值。因為長度頻率數(shù)據(jù)不包含絕對豐度的信息，當?shù)仁?(2) 的兩邊除以它們各自的和時，等式不變。從右側分母的和中提取常量NLstart，然后在分子中消掉NLstart，剩下的2個要確定的參數(shù)是Z/k和L∞。

通過方程式 (4) 和 (5) 分別得出未開發(fā)群體生物量最大長度Lopt以及在給定的捕撈壓力下，使捕獲量與生物量最大化的第一次捕撈長度，即最適開捕體長Lc_opt：

單位補充量指數(shù)可以是Lc/L∞、F/k、M/k與捕撈強度 (F/M)的函數(shù)[23]，Lc指50%選擇性體長：

單位努力量指數(shù) (C PUE'/R) 由式 (6) 除以F/M得到，假設捕撈死亡率F與捕撈努力量成正比。由于CPUE與資源開發(fā)階段的生物量成正比，即相對 C PUE'/R與單位補充B'/R的開發(fā)生物量指標可由 (8) 計算得出[23]：

若不進行捕撈，則種群在開發(fā)階段的相對生物量為：

其中，B0'＞Lc為未捕撈生物量B0的可開發(fā)部分 (Lc)，種群B/B0的被開發(fā)部分的相對生物量指數(shù)為：

LBB估算在Bayesian Gibbs sampler軟件JAGS中實施[24]，并使用統(tǒng)計語言R執(zhí)行，代碼(R-code: LBB_20.R) 可從網(wǎng)站http://oceanrep.geomar.de/44832/下載。

2 結果

2.1 體長組成

2006—2014年冬季總體長范圍介于25～290 mm，優(yōu)勢體長范圍介于140～149 mm，平均體長為137 mm；夏季總體長范圍介于6.9～240 mm，優(yōu)勢體長范圍介于110～119 mm，平均體長為122 mm，優(yōu)勢體長范圍比冬季小(表2)?？傮w而言，冬季漁獲的個體比夏季漁獲的個體大。

表2 2006—2014年大頭白姑魚群體結構Table 2 Population structure of P. macrocephalus in summer and winter from 2006 to 2014

2.2 種群參數(shù)

估算的大頭白姑魚L∞介于213～307 mm，平均值為264.50 mm；Lc_opt介于114～226 mm，平均值為156.50；M/k介于0.15～1.74，平均值為1.39；F/k介于0.12～20，平均值為3.58；Z/k介于1.04～21.6，平均值為4.96；開發(fā)率 (E)介于0.12～0.96，平均值為 0.58 (表 3)。其中，L∞、Lc_opt、M/k、E變異系數(shù)均小于0.5，F(xiàn)/k、Z/k變異系數(shù)均大于0.5。

表3 2006—2014年根據(jù)LBB估算的大頭白姑魚的種群參數(shù)Table 3 Population parameters of P. macrocephalus estimated by LBB from 2006 to 2014

圖1是估算的2006—2014年大頭白姑魚種群參數(shù)的年際變化線性趨勢圖，其數(shù)值為各年冬季和夏季的平均值。從中可以看出，2006—2014年大頭白姑魚的L∞呈下降趨勢，年際變化無顯著性差異 (P=0.50＞0.05)；Lc_opt波動較小，略有上升趨勢(P=0.30＞0.05)。隨著大頭白姑魚F/k的明顯增加(P=0.07＞0.05)，E呈上升趨勢，但年際間變化不顯著 (P=0.11＞0.05)；M/k值在1.5上下波動，呈下降趨勢 (P=0.48＞0.05)。

圖1 2006—2014年北部灣大頭白姑魚種群參數(shù)的變化趨勢Figure 1 Trends of population parameters of P. macrocephalus in Beibu Gulf from 2006 to 2014

2.3 資源狀況

圖2-a中顯示2006—2014年北部灣大頭白姑魚的F/M呈現(xiàn)上升趨勢，且捕撈死亡率F大于自然死亡率M。捕撈強度的逐年增加使北部灣大頭白姑魚資源呈下降趨勢 (圖2-b)，且B/B0基本小于0.5，資源衰退明顯。

圖2 大頭白姑魚捕撈強度 (F/M) 與資源利用狀態(tài)Figure 2 Fishing intensity (F/M) and stock utilization status of P. macrocephalus in Beibu Gulf from 2006 to 2014

3 討論

北部灣地處熱帶-亞熱帶，棲息于該海域的魚類具有生長快、生命周期較短和個體較小的生物學特性。用傳統(tǒng)方法準確鑒定其年齡 (日齡) 較為困難，且容易產(chǎn)生較大的誤差。長度頻率法的運用解決了鑒定年齡的困難，克服了推算年齡和生長參數(shù)存在的主觀誤差，節(jié)省了大量的人力和物力。

3.1 群體體長組成

在漁業(yè)資源評估中，當技術有限、數(shù)據(jù)不充分或目標物種習性復雜等成為準確評估漁業(yè)資源的限制性因素時[25]，衍生出數(shù)據(jù)受限、數(shù)據(jù)有限或數(shù)據(jù)缺乏 (Data poor or data limited) 的漁業(yè)資源評估。本研究中，由于不需要任何除體長數(shù)據(jù)之外的信息輸入，這給一些數(shù)據(jù)缺乏的漁業(yè)資源評估提供了參考。但是，只用體長數(shù)據(jù)并不是最佳選擇，利用體長數(shù)據(jù)計算的體長頻率在漁獲物取樣的過程中，可能會出現(xiàn)體長頻率分布的偏差，樣本中小個體偏少或者大個體缺失都會對估算到的L∞造成影響[26]，其估計值與真實值存在差異。其次，個體的變異也對長度頻率數(shù)據(jù)估計種群參數(shù)存在影響。在個體增長差異較小且補充量統(tǒng)一情況下，廣義方法的估算值偏差較大[27]。本研究中年平均體長呈下降趨勢，表明大頭白姑魚呈體長變小的特點?？紤]到捕撈對體長分布的影響，有學者提出為了減少捕撈對一個種群體型分布的影響，捕撈死亡率不得超過任何體型且成年種群的平均自然死亡率[23]。此外，捕撈漁具的選擇性也會產(chǎn)生體長組成的偏差，從而影響評估結果。

3.2 不同方法估算的種群參數(shù)的比較

在南海北部，與大頭白姑魚同屬且較為常見的還有白姑魚、斑鰭白姑魚和截尾白姑魚 (P. anea)。其中，白姑魚與大頭白姑魚形態(tài)相似、個體大小相近，在近海?；鞐?，容易將二者混淆。表4為已估算的白姑魚屬不同種類的種群參數(shù)，L∞介于220.32～305.00 mm，不同方法的估算結果存在一定差異。本研究根據(jù)LBB估算的L∞為264 mm，介于顏云榕等[14]和Wang等[28]分別用耳石生長輪和體長頻率方法估算的L∞(分別為229.95和289 mm)之間。運用LBB估算的北部灣大頭白姑魚的E與Wang等[28]通過ELEFAN I技術估算的結果大體一致，分別為0.58與0.63。

表4 白姑魚屬種群參數(shù)的估算結果Table 4 Estimation results of population parameters of Pennahia

體長頻數(shù)分析法廣泛應用于熱帶-亞熱帶魚種的資源評估中。ELEFAN I技術、最大似然估計法和LBB都是在生長方程的基礎上發(fā)展而來的，不同的是為了模擬生長，ELEFAN I通過移動平均數(shù)重組大小頻率以減少不規(guī)則性，然后用模態(tài)級數(shù)分析擬合生長曲線[29]；LBB基于體長數(shù)據(jù)通過蒙特卡羅馬爾科夫鏈估算所有參數(shù)。ELEFAN I和LBB都是以體長為基礎的方法，但在實際漁業(yè)中，獲得毫無偏差的樣本十分困難[30]，造成偏差的原因可能與采樣、魚類的生物學特性或魚類的洄游相關。

3.3 大頭白姑魚資源開發(fā)狀態(tài)分析

陳丕茂[33]運用ELEFAN I技術根據(jù)1997—1999年南海北部 (包括北部灣) 白姑魚的生物學數(shù)據(jù)，建議南海北部灣白姑魚的Lc_opt為160 mm，此時調查的最小性成熟體長為158 mm。Wang等[28]估算2006—2007年北部灣大頭白姑魚的Lc_opt為174 mm。本研究根據(jù)2006—2014年的北部灣大頭白姑魚的長度頻率數(shù)據(jù)，運用LBB估算其Lc_opt為157 mm。此估算結果與陳丕茂[33]估算的160 mm接近，而略低于Wang等[28]估算的174 mm?？紤]到開捕規(guī)格越小，在實際工作中更容易推行，且此時的調查最小性成熟體長為130 mm[14]，即建議的開捕體長大于調查最小性成熟體長。另一方面，研究表明南海北部 (北部灣) 經(jīng)濟魚類的性成熟有進一步提前的趨勢。綜上，建議北部灣大頭白姑魚的Lc_opt為 157 mm。

以往對北部灣白姑魚屬自然死亡系數(shù)M的研究均估算其絕對的數(shù)值，LBB用相對自然死亡率M/k替代其絕對數(shù)據(jù)。因為M/k的變化小于參數(shù)本身，且該比值可以從生活史中近似得到[23]。因此，使用M/k代替單個參數(shù)M和k更便捷有利。Froese等[23]認為對于生長不確定的物種，M/k應介于1.0～2.0，如果M/k小于1.5，則認為自然死亡率相對較低。本研究估算的北部灣2006—2014年大頭白姑魚M/k的平均值為1.39，表明北部灣大頭白姑魚的自然死亡率較低，即捕撈死亡較高。從2007年開始，中國南部海域捕撈機動漁船數(shù)量比較穩(wěn)定，但是越南環(huán)北部灣的海洋機動船數(shù)量和功率急速增加，盡管中國出臺一系列相應政策控制近海捕撈，但是越南環(huán)北部灣作業(yè)漁船增加且沒有休漁政策，北部灣漁業(yè)仍處于較高捕撈強度狀態(tài)，導致漁業(yè)資源密度不斷下降[34]。F/M=2.58＞1也從側面反映了該魚種的過度捕撈[34]。其次，學者們一般認為B/B0＜0.5時資源過度開發(fā)[35]，從圖2-b可知除2008年，剩余年份B/B0均小于0.5，說明北部灣大頭白姑魚處于生物量極低的狀態(tài)。此外，Gulland[36]認為一般魚類的最適開發(fā)率為0.5，E介于0～0.5表示輕度開發(fā)，0.5～1.0表示過度開發(fā)，并據(jù)此來判斷某個海域漁業(yè)資源的開發(fā)程度。本研究的E=0.58＞0.5，以此為標準，大頭白姑魚已處于過度開發(fā)狀態(tài)。

3.4 LBB方法的優(yōu)缺點

基于長度的種群評估模型，一般有3個最直觀的優(yōu)點：1) 體長數(shù)據(jù)易獲??；2) 模型真實性更高；3) 充分利用了體長組成數(shù)據(jù)。LBB被認為是一種合理的資源評估方法[37]，可直接利用長度頻率評估數(shù)據(jù)信息較差的魚類資源。但是此種方法也存在一定局限性，由于只輸入長度頻率數(shù)據(jù)，LBB在資源評估中存在估算結果不理想的情況[35]。這可能是樣本量的大小對長度組成數(shù)據(jù)質量存在影響造成的[38]。由表3可知2014年種群參數(shù)指標波動較大，這可能與樣本少有關。其次該模型較為復雜，參數(shù)不易估算而且無法估算年齡結構。如果LBB和其他方法結合使用，評估效果可能會更準確。

本研究通過多個指標評價得出北部灣大頭白姑魚資源處于過度利用狀態(tài)，需制定人為干預方案來保持北部灣大頭白姑魚資源的可持續(xù)利用。從前人研究結果看，北部灣漁業(yè)資源的變動與環(huán)境因素顯著相關[39]，魚類生長常常受溫度、鹽度和其他環(huán)境因素影響。而本研究尚未考慮除體長之外影響漁業(yè)資源變動的其他因素。因此，本研究對大頭白姑魚的種群資源評估存在一定的局限性。其次，本研究采取較通用的10 mm間距進行體長頻率時間序列數(shù)據(jù)整理，不同分組組距構成的體長頻率數(shù)據(jù)來估算種群參數(shù)有明顯差異[40]。后續(xù)還應深入研究對比不同組距的估算結果，進而提高種群參數(shù)估算結果的可信度。