基于貝葉斯狀態(tài)空間產(chǎn)量模型的大西洋黃鰭金槍魚資源評估

田志盼，田思泉,2,3，戴黎斌，麻秋云,2,3*

(1.上海海洋大學 海洋科學學院，上海 201306；2.國家遠洋漁業(yè)工程技術(shù)研究中心，上海 201306；3.大洋漁業(yè)資源可持續(xù)開發(fā)教育部重點實驗室，上海 201306)

1 引言

黃鰭金槍魚（Thunnus albacares）是高度洄游類的魚種，廣泛分布于三大洋的熱帶和亞熱帶海域，資源量相對豐富且價值較高，是金槍魚漁業(yè)中重要的經(jīng)濟魚種。隨著人類需求的逐漸上升，其受到的捕撈威脅也不斷增加，因此進行科學的資源評估進而制定合理的養(yǎng)護管理措施，是實現(xiàn)漁業(yè)可持續(xù)開發(fā)的基礎(chǔ)。

目前，國際上對公海金槍魚實施管轄的是各類區(qū)域性漁業(yè)管理組織（Regional Fisheries Management Organizations,RFMOs)，大西洋黃鰭金槍魚由國際大西洋金槍魚養(yǎng)護委員會（International Commission for Conservation of Atlantic Tunas,ICCAT）進行管理。ICCAT 在對黃鰭金槍魚的資源評估中，主要使用的評估模型有非平衡剩余產(chǎn)量模型（A Stock Production Model Incorporating Covariates,ASPIC）[1]、年齡結(jié)構(gòu)產(chǎn)量模型（Age-Structured Production Model,ASPM）[2]、實際種群分析（Virtual Population Analysis,VPA）[3]和資源綜合模型（Stock Synthesis III,SS3）[4]等，每個模型評估的資源狀態(tài)不盡相同，ICCAT 綜合各模型認為，當前資源處于資源型過度捕撈而無捕撈型過度捕撈狀態(tài)。其中，剩余產(chǎn)量模型，相較其他模型來說，對漁業(yè)數(shù)據(jù)需求較低（僅需漁獲量和種群豐度指數(shù)），且能得到最大可持續(xù)產(chǎn)量（Maximum Sustainable Yield,MSY）等參考點信息，故是各類漁業(yè)資源評估中使用最廣泛的模型之一[5]。

除了觀測誤差，漁業(yè)資源種群動態(tài)中還存在環(huán)境變化等因素產(chǎn)生的過程誤差，而ICCAT 當前使用的ASPIC 等產(chǎn)量模型無法估計過程誤差，由該類誤差產(chǎn)生的不確定性難以被考慮在內(nèi)。JABBA（Just Another Bayesian Biomass Assessment）是一種基于貝葉斯方法的狀態(tài)空間產(chǎn)量模型[6]，其中貝葉斯框架可以通過合理的先驗信息來降低模型中的不確定性[7-8]，狀態(tài)空間建模則可以同時估計過程誤差和觀測誤差[9-11]。為此，本文嘗試采用JABBA 模型來評估大西洋黃鰭金槍魚的資源狀況，研究狀態(tài)空間建模對該資源進行評估的適用性，以期為該重要魚種的科學研究和漁業(yè)管理提供更多基礎(chǔ)資料和參考信息。

2 材料與方法

2.1 數(shù)據(jù)來源

本文利用的1950?2017 年漁獲量數(shù)據(jù)來自ICCAT 數(shù)據(jù)庫，漁獲量在1990 年達到最高的19.36×104t，2017 年為13.53×104t（圖1）。單位捕撈努力量漁獲量（Catch Per Unit Effort，CPUE）數(shù)據(jù)來自ICCAT 黃鰭金槍魚資源評估會議報告和CPUE 標準化研究報告[12-16]，共計8 個延繩釣船隊的標準化CPUE 數(shù)據(jù)（表1）。雖然各船隊根據(jù)各自漁業(yè)和數(shù)據(jù)情況采用了不同的CPUE 標準化方法（日本、委內(nèi)瑞拉和美國為廣義線性模型，而烏拉圭1、烏拉圭2、巴西、中國臺北1 和中國臺北2 則為廣義線性混合模型），但應(yīng)ICCAT 相關(guān)工作組的建議，本文資源評估模型在建模時納入了所有8 個CPUE 數(shù)據(jù)。

2.2 JABBA 模型

JABBA（版本為v1.1[6]）中運行的Pella-Tomlinson 剩余產(chǎn)量函數(shù)形式如下：

式中，SP為剩余產(chǎn)量；r為種群的內(nèi)稟增長率；K為平衡狀態(tài)時的未開發(fā)資源生物量（即環(huán)境容納量）；B為資源量；m為形狀參數(shù)。

獲得MSY 時的生物量BMSY和捕撈死亡率FMSY分別為

圖1 大西洋黃鰭金槍魚1950?2017 年的年漁獲量Fig.1 The annual catch of Atlantic yellowfin tuna from 1950 to 2017

表1 大西洋黃鰭金槍魚各延繩釣船隊標準化CPUE 數(shù)據(jù)Table 1 Standardized CPUE for each longline fleet of Atlantic yellowfin tuna

捕撈死亡率F定義為

MSY 定義為

B/BMSY＜1 表示當前種群已發(fā)生資源型過度捕撈，F(xiàn)/FMSY＞1 表示種群正遭受捕撈型過度捕撈。如果m=2，則SP函數(shù)為Scheafer 形式；如果m趨近于1，則為Fox 形式。在JABBA 中，采用默認設(shè)置值BMSY/K=0.4，由此得出m=1.2。

過程方程定義如下：

式中，y為年份，Py為y年B與K的比值；ηy為過程誤差，且為過程方差，服從逆伽馬分布（inverse-gamma (4,0.01)）；Cf,y為y年船隊f的漁獲量。

JABBA 中觀測方程定義如下：

式中，qi為豐度指數(shù)i的可捕性系數(shù)；εy,i為觀測誤差，且為觀測方差，包含固定項和預(yù)測項，預(yù)測項服從無信息的逆伽馬分布（inverse-gamma(0.001,0.001)）。

本研究中各參數(shù)的先驗分布設(shè)置如下：BMSY/K=0.4；σfix=0.2；初始資源消耗率B1950/K服從對數(shù)正態(tài)分布，其中值和變異系數(shù)分別為1.0 和0.1；r和K的先驗信息參考Matsumoto 等[1]的研究結(jié)果：假設(shè)r服從0.14～0.34 的均勻分布，K服從139.2×104～265.8×104的均勻分布；可捕性系數(shù)q為無信息均勻分布。

因Schaefer 的對稱形式不符合黃鰭金槍魚種群動態(tài)變化情況[17]，本文只考慮選擇Fox 和Pella-Tomlinson 函數(shù)。根據(jù)不同的CPUE 數(shù)據(jù)和剩余產(chǎn)量函數(shù)，預(yù)實驗共設(shè)置了S1?S8 共8 種方案進行分析（表2）。當均方根誤差（Root Mean Squared Error，RMSE）或偏差信息準則（Deviation Information Criteria，DIC）較小時，說明模型擬合效果較好[18]。

當選擇所有CPUE 數(shù)據(jù)并使用Pella-Tomlinson 函數(shù)時，得到S1；在S1 預(yù)實驗基礎(chǔ)上，去掉擬合效果差的CPUE 數(shù)據(jù)得到了S2；在S2 基礎(chǔ)上，考慮到ICCAT在CPUE 數(shù)據(jù)方面的建議[12]—因CPUE 標準化當中未考慮目標魚種的變化，認為日本延繩釣標準化CPUE數(shù)據(jù)應(yīng)該從1976 年開始，舍棄之前的數(shù)據(jù)得到S3；S3 的CPUE 數(shù)據(jù)不變，剩余產(chǎn)量函數(shù)選擇Fox 得到S4；在S4 基礎(chǔ)上，考慮對CPUE 數(shù)據(jù)的敏感性，依次去掉1 條CPUE 數(shù)據(jù)后得到S5?S8 共4 種方案。

表2 大西洋黃鰭金槍魚JABBA 模型S1?S8 方案設(shè)置Table 2 Different scenarios (S1?S8) of Atlantic yellowfin tuna in JABBA

2.3 回溯性分析

隨著漁業(yè)數(shù)據(jù)逐年增加到資源評估中，模型估算結(jié)果可能因為出現(xiàn)系統(tǒng)性偏差而導(dǎo)致持續(xù)高估或低估的問題稱為回溯性問題（Retrospective Problem，RP）。RP 誤差的強度主要由Mohn[19]定義的ρ來衡量：

式中，y1、y2分別為整個數(shù)據(jù)的起始年和結(jié)束年，y1：y表示利用y1到y(tǒng)年的數(shù)據(jù)進行模型估計；X為某一估計的模型參數(shù)（如資源生物量或捕撈努力量等）。

如果ρ趨于0，則表明不存在RP；ρ大于0，則存在正RP，即同一年某參數(shù)短時間序列的估計值大于整個時間序列的估計值，反之則為負RP[20]。

2.4 敏感性分析

本文通過敏感性分析，研究了種群關(guān)鍵參數(shù)K和r的先驗分布以及漁獲量數(shù)據(jù)的誤報比例對評估結(jié)果的影響，進而探討模型的穩(wěn)健性。本文分別研究K和r無信息的先驗分布和有信息的先驗分布（表3）。20 世紀90 年代中后期ICCAT 數(shù)據(jù)收集上報過程才更加規(guī)范[21-22]，這段時間前后數(shù)據(jù)的可信度存在差異。鑒于此，本文假設(shè)1950?1994 年間漁獲量數(shù)據(jù)存在不同程度的誤報問題，即上報漁獲量占實際漁獲量的比例分別設(shè)為70%、80%、90%、110%、120%和130%共6 種情況。

2.5 預(yù)測分析

ICCAT 在2016 年大西洋黃鰭金槍魚資源評估會議中，預(yù)測分析顯示，漁獲量低于12×104t 時能使種群到2024 年一直保持健康狀態(tài)，所以將其總允許可捕量（Total Allowable Catch，TAC）設(shè)定為11×104t[22]。因此本研究以11×104t 為基礎(chǔ)，設(shè)置8.80×104t (80%)、9.35×104t (85%)、9.90×104t (90%)、10.45×104t (95%)、11.00×104t (100%)、12.10×104t (110%)和13.20×104t(120%)共7 個TAC 指標，假設(shè)2018 年漁獲量為2015?2017 年的平均值，并以2019 年為起始年，預(yù)測2019?2027 年的種群動態(tài)變化，并以生物量B＞BMSY及種群處于健康狀態(tài)等的概率評價TAC 指標的管理效果。

3 結(jié)果

3.1 模型結(jié)果

各方案下得到了JABBA 模型的CPUE 指數(shù)趨勢（圖2）和擬合優(yōu)度（表4）。在S1 方案下，URU1、URU2、BR 和TAI2 的擬合效果非常差（圖2a），存在許多異常值且RMSE 極高（表4），而去除異常值后，擬合效果有較大改善，RMSE 大幅降低（圖2b）。S3使用JAP_RE 后，擬合效果略有改善，S4 改用Fox 函數(shù)后RMSE 基本不變但DIC 降低。S5、S6、S8 下的擬合效果稍有提升，S7 則變差（表4）。綜上所述，且鑒于S4 方案涵蓋了更多船隊的CPUE 信息，本文將S4 方案作為基礎(chǔ)模型來提供資源評估結(jié)果。

表3 大西洋黃鰭金槍魚JABBA 模型中K 和r 有信息和無信息的先驗分布設(shè)定以及后驗分布Table 3 The informative and non-informative prior and posterior distributions for K and r in the JABBA for Atlantic yellowfin tuna

基礎(chǔ)模型所有參數(shù)后驗分布均左右對稱且在合理的范圍內(nèi)，說明模型收斂并得到了可靠的結(jié)果（圖3）。本文求得大西洋黃鰭金槍魚的MSY 為13.7×104t，BMSY為65.2×104t，B2017略高于BMSY，F(xiàn)2017略低于FMSY（表5）。

3.2 資源狀態(tài)

隨著漁業(yè)開發(fā)程度的增加，種群由初始（1950 年）的健康狀態(tài)逐漸進入資源捕撈過度的狀態(tài)（1997 年前后），隨后逐漸恢復(fù)，2017 年資源有65%的概率既沒有處于資源型過度捕撈狀態(tài)，也沒有處于捕撈型過度捕撈狀態(tài)，資源狀態(tài)健康（圖4）。20 世紀90 年代和21 世紀初期，種群處于過度捕撈狀態(tài)（相對捕撈死亡率F/FMSY＞1，而相對生物量B/BMSY＜1（圖5））。

以2019 年為起始年，在7 個不同的TAC 目標下，預(yù)測分析顯示，2019?2027 年資源量均保持增長的趨勢（圖6）。當TAC 為8.8×104t 時，生物量增長最快，隨著TAC 變大，資源量增長速度放緩。風險分析結(jié)果顯示（表6至表8），TAC為11×104t時，2024年B＞BMSY和資源健康的概率均為84.6%，F(xiàn)＞FMSY的概率為6.2%；當TAC 為13.2×104t 時，2024 年B＞BMSY的概率降低到了69.2%（表6），但F＞FMSY的概率明顯增大（28.9%，表7）。

敏感性分析結(jié)果表明，當K為無信息先驗時，K的估計值略有減小，r估計值略有增大；當r無信息先驗時，K的估計值略有增大，r的估計值略有減小（表3）。F2017隨漁獲量少報程度增大而減小，B2017則與之相反；F2017隨多報程度增大而減小，B2017則與之相反（表9）。但不同誤報率情況下，B2017/BMSY都大于1，F(xiàn)2017/FMSY都小于1，且資源處于健康狀態(tài)的概率變化較小。回溯性分析結(jié)果表明，當數(shù)據(jù)逐年減少至2012 年，B、B/BMSY估計值略有減小，F(xiàn)、F/FMSY估計值略有增大，但差別極小（圖7）。計算得到B、B/BMSY、F、F/FMSY的ρ值分別為?0.360、0.448、?0.183、0.296。

4 討論

本文通過貝葉斯狀態(tài)空間的剩余產(chǎn)量模型，在JABBA 中評估了1950?2017 年大西洋黃鰭金槍魚的資源狀況。當前種群處于沒有過度捕撈的健康狀態(tài)，在ICCAT 當前的TAC 養(yǎng)護管理措施下，2024 年能夠達成其保持種群健康狀態(tài)的養(yǎng)護管理目標。研究結(jié)果表明，當使用美國、委內(nèi)瑞拉、日本（去除1976 年以前）、中國臺北1993?2014 年4 個CPUE 數(shù)據(jù)及Fox函數(shù)時，JABBA 模型的擬合效果最佳，評估結(jié)果對參數(shù)K、r的先驗分布和1950?1994 年間的漁獲量誤報不太敏感，且模型不存在明顯的回溯性誤差。

1994 年ICCAT 成立工作組對大西洋黃鰭金槍魚進行評估[22]，之后分別在2000 年、2003 年、2008 年、2011 年和2016 年都進行了資源評估。在2016 年的資源評估中，ASPIC 模型評估認為2014 年大西洋黃鰭金槍魚處于資源型過度捕撈狀態(tài)，但沒有遭受捕撈型過度捕撈[1]；SS3 和VPA 模型認為其處于上述兩種過度捕撈狀態(tài)[4,12]；而ASPM 模型則表明其均不處于過度捕撈狀態(tài)[2]。綜合上述模型，ICCAT 認為種群處于資源型過度捕撈而未遭受捕撈型過度捕撈的狀態(tài)[12]。本研究與之產(chǎn)生差異的原因可能是由狀態(tài)空間建模與上述幾種模型的結(jié)構(gòu)差異及使用的先驗信息不同所導(dǎo)致的。本研究得到大西洋黃鰭金槍魚的環(huán)境容納量K為178×104t，內(nèi)稟增長率r為0.210，與同樣基于剩余產(chǎn)量理論構(gòu)建的ASPIC 模型的結(jié)果相似[1]，說明評估結(jié)果較為可信；與ASPIC 模型相比，評估的種群狀態(tài)脫離了資源型過度捕撈，這可能是因為近兩年捕撈死亡率的下降，使黃鰭金槍魚資源有機會得到部分恢復(fù)。

圖2 大西洋黃鰭金槍魚JABBA 模型S1?S8 方案的CPUE 指數(shù)趨勢Fig.2 Time-series of input CPUE of Atlantic yellowfin tuna and predicted CPUE of S1?S8 scenarios in JABBA

表4 大西洋黃鰭金槍魚JABBA 模型S1?S8 方案的擬合效果Table 4 Goodness of fitting of S1?S8 scenarios in JABBA for Atlantic yellowfin tuna

圖3 大西洋黃鰭金槍魚JABBA 基礎(chǔ)模型參數(shù)先驗分布（深色）和后驗分布（淺色）Fig.3 Priors (dark) and posteriors (light) of parameters of base case in JABBA for Atlantic yellowfin tuna

相對ICCAT 當前的資源評估模型而言，JABBA作為剩余產(chǎn)量模型的一種，其結(jié)果可靠性較高但無法充分利用魚類的生物學數(shù)據(jù)，與ICCAT 使用的其他剩余產(chǎn)量模型如ASPIC 等相比，JABBA 可以估計過程誤差，對形狀參數(shù)的估計更自由，但JABBA 的基本假設(shè)為漁獲量不存在誤差，這一點有待改進。此外，在v1.1 版本的JABBA 中，在選擇Pella-Tomlinson 產(chǎn)量函數(shù)時，模型無法將m作為未知參數(shù)直接估算，必須通過假定BMSY與K的關(guān)系得到，因此下一步我們將探尋JABBA 的更高版本，以研究m的估算問題。JABBA 模型參數(shù)設(shè)定自由，擬合快速，當前已有用其評估大西洋劍魚（Xiphias gladius）[23]、大西洋大眼金槍魚（Thunnus obesus）[24]、大西洋藍槍魚（Makaira nigric-ans）[25]等的研究，相信其在RFMOs 的資源評估中將發(fā)揮越來越重要的作用。

表5 大西洋黃鰭金槍魚JABBA 基礎(chǔ)模型參數(shù)后驗估計值及其95%置信區(qū)間Table 5 Posterior estimates and 95% confidence intervals of parameter of base case in JABBA for Atlantic yellowfin tuna

近年來大西洋黃鰭金槍魚的漁獲量為13×104t 左右，預(yù)測分析結(jié)果表明，當前TAC（11×104t）管理措施對其種群的養(yǎng)護是有效的，可以實現(xiàn)ICCAT 的管理目標，而在當前的漁獲量水平下，種群生物量仍能保持一定的速度增長。近年來，我國的漁獲量僅為0.05×104t 左右[26]，占總漁獲量比例較小，且都來自于延繩釣漁業(yè)的兼捕漁獲，因此我國的大西洋黃鰭金槍魚漁業(yè)仍有一定的開發(fā)空間。

大西洋黃鰭金槍魚漁業(yè)主要有圍網(wǎng)、延繩釣和餌釣3 種，其中圍網(wǎng)漁業(yè)占總漁獲量的70%左右，且圍網(wǎng)漁業(yè)主要在東部大西洋作業(yè)[26]。而黃鰭金槍魚的生長分為兩個階段[27-28]，幼魚階段生長較為緩慢，成魚階段生長快速，且黃鰭金槍魚幼魚主要在大西洋東部完成早期生活史[29]。當前圍網(wǎng)漁業(yè)漁獲量上升[26]，造成的黃鰭金槍魚幼魚死亡率偏高[12]，可能導(dǎo)致補充量不足，種群內(nèi)稟增長率降低，致使剩余產(chǎn)量減少，可以考慮適當限制圍網(wǎng)漁業(yè)的捕撈投入，以更好地養(yǎng)護大西洋黃鰭金槍魚資源。

圖4 1950?2017 年大西洋黃鰭金槍魚JABBA 模型基礎(chǔ)模型資源開發(fā)狀態(tài)變化圖Fig.4 Kobe phase plot showing estimated trajectories(1950?2017) of B/BMSYand F/FMSYof Atlantic yellowfin tuna of base case in JABBA

圖5 大西洋黃鰭金槍魚JABBA 基礎(chǔ)模型1950?2017 年F/FMSY和B/BMSY趨勢Fig.5 F/FMSYand B/BMSYof Atlantic yellowfin tuna from 1950 to 2017 of base case in JABBA

圖6 不同TAC 目標下的大西洋黃鰭金槍魚JABBA 模型基礎(chǔ)模型B/BMSY預(yù)測（2019?2027 年）Fig.6 Future projection (2019?2027) of B/BMSYof Atlantic yellowfin tuna of base case in JABBA under different TACs

剩余產(chǎn)量模型將種群所有個體生命史的動態(tài)變化過程進行了高度綜合，模型具有參數(shù)少、所需數(shù)據(jù)相對簡單的特點，而形狀參數(shù)較難準確估計且容易導(dǎo)致資源評估的失敗[30]，因此在模型擬合結(jié)果相差不大時，本研究最終選擇了較簡單的Fox 而放棄形狀參數(shù)不易估計的Pella-Tomlinson 產(chǎn)量函數(shù)形式。貝葉斯方法把經(jīng)驗判斷、前人的研究結(jié)果與現(xiàn)有數(shù)據(jù)相結(jié)合[8,31]，后驗概率分布由先驗概率分布和模型數(shù)據(jù)共同決定，但如果所用的數(shù)據(jù)不包含足夠的信息，那么后驗概率分布可能完全由先驗概率主導(dǎo)和控制[32]。因此在使用貝葉斯資源評估方法時，對后驗概率分布與先驗概率分布進行比較分析顯得尤為重要[33-34]。本研究中的敏感性分析顯示，K、r的后驗分布對先驗分布是否有信息并不敏感，說明數(shù)據(jù)為模型的貝葉斯方法提供了足夠的信息。

回溯性誤差在漁業(yè)資源評估中比較普遍，誤差過大可能導(dǎo)致漁業(yè)管理的失敗[35]。對基礎(chǔ)模型的回溯性分析中，對生物量的估計過低而對捕撈死亡率估計過高，這可能是由近年來黃鰭金槍魚漁獲量下降導(dǎo)致。4 個參數(shù)的ρ值均趨近于0，結(jié)合圖形繪制結(jié)果可以表明不存在明顯的回溯性誤差，這可能是由于狀態(tài)空間建模不僅給出了傳統(tǒng)模型的點估計值，同時能量化觀測誤差和過程誤差的不確定性，從而避免了一定的回溯性問題[20,36]。

本研究表明，早期漁獲量數(shù)據(jù)誤報率會對資源量和捕撈死亡率的結(jié)果產(chǎn)生一定影響，而種群狀態(tài)并沒有明顯改變，不會影響對種群健康狀態(tài)的判斷。一般來說漁獲量數(shù)據(jù)以少報居多，本研究表明此時資源評估的結(jié)果將更加樂觀。但本研究未考慮其他時間段內(nèi)漁獲量數(shù)據(jù)失真問題，而近期漁獲量數(shù)據(jù)對當前資源狀態(tài)的判斷有更大影響，此外，下一步的研究還應(yīng)考慮近期數(shù)據(jù)的隨機誤差等情況[37]。

表6 不同TAC 目標下大西洋黃鰭金槍魚2019?2027 年B>BMSY的概率Table 6 The probability that B>BMSYof Atlantic yellowfin tuna under different TAC targets in 2019?2027

表7 不同TAC 目標下大西洋黃鰭金槍魚2019?2027 年F>FMSY的概率Table 7 The probability that F>FMSYof Atlantic yellowfin tuna under different TAC targets in 2019?2027

表8 不同TAC 目標下大西洋黃鰭金槍魚2019?2027 年處于健康狀態(tài)的概率Table 8 The probability that the Atlantic yellowfin tuna is in healthy status under different TAC targets in 2019?2027

表9 不同漁獲量誤報比例下大西洋黃鰭金槍魚JABBA 基礎(chǔ)模型評估資源狀態(tài)Table 9 Stock status of Atlantic yellowfin tuna in different mis-reported rates of catches of base case in JABBA

圖7 大西洋黃鰭金槍魚JABBA 基礎(chǔ)模型B、B/BMSY、F、F/FMSY的回溯性分析Fig.7 Retrospective analysis of B,B/BMSY,F,F/FMSYof base case in JABBA of Atlantic yellowfin tuna

致謝：感謝漁業(yè)資源和生態(tài)系統(tǒng)量化評估與管理研究室趙蓬蓬等師兄在論文修改方面的幫助。