利用后向預(yù)報方法分析印度洋黃鰭金槍魚資源評估模型?

官文江, 吳佳文, 曹友華

(1. 上海海洋大學(xué)海洋科學(xué)學(xué)院，上海 201306； 2. 上海海洋大學(xué)大洋漁業(yè)資源可持續(xù)開發(fā)省部共建教育部重點(diǎn)實驗室， 上海 201306)

印度洋黃鰭金槍魚(Thunnusalbacares)是商業(yè)金槍魚漁業(yè)的重要捕撈目標(biāo)魚種之一，主要被法國、西班牙的圍網(wǎng)漁業(yè)，印度尼西亞、中國臺灣的冰鮮延繩釣漁業(yè)，中國臺灣、日本、韓國的超低溫延繩釣漁業(yè)及馬爾代夫、印度等國的手釣、竿釣、刺網(wǎng)等漁業(yè)所利用[1]。近年，隨著印度洋海盜活動減少、安全形式好轉(zhuǎn)，印度洋黃鰭金槍魚漁獲量又得以大幅提升，年平均產(chǎn)量超過了4.0×105t。盡管我國捕撈的印度洋黃鰭金槍魚較少，但在我國的印度洋金槍魚漁業(yè)中，其產(chǎn)量僅次于大眼金槍魚(Thunnusobesus)與長鰭金槍魚(Thunnusalalunga)。印度洋黃鰭金槍魚的資源狀況也日益受到我國學(xué)者的關(guān)注[1-5]。

當(dāng)前，剩余產(chǎn)量模型(Bayesian Biomass Dynamic Model，BBDM)[1]、年齡結(jié)構(gòu)模型(Statistical Catch at Age，SCAA)[6]及合成模型(Stock Synthesis，SS3)[7-8]被用于印度洋黃鰭金槍魚的資源評估。盡管SS3模型的評估結(jié)果被用于確定印度洋黃鰭金槍魚的資源狀態(tài)[9]，但上述各模型的評估結(jié)果均存在較大不確定性[1,9]。特別是這些模型的診斷與選擇，主要采用的依據(jù)是模型對觀察數(shù)據(jù)(主要是資源豐度指數(shù)及年齡或體長組成數(shù)據(jù))的擬合度。但對漁業(yè)管理而言，如總可捕量(Total Allowable Catch，TAC)的確定等，模型的預(yù)測能力也非常重要。在當(dāng)前的漁業(yè)資源評估中，相關(guān)學(xué)者并未足夠重視模型的預(yù)測能力，很少評價模型的預(yù)測能力、并將其作為評價漁業(yè)資源評估與管理質(zhì)量的依據(jù)[10]。為此，本文以印度洋黃鰭金槍魚資源評估為例，利用后向預(yù)報方法評價模型的預(yù)測能力，以分析各模型漁業(yè)資源評估與管理質(zhì)量，為科學(xué)選擇漁業(yè)資源評估模型、避免漁業(yè)管理風(fēng)險提供科學(xué)依據(jù)。

1 材料與方法

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

用于印度洋黃鰭金槍魚資源評估的數(shù)據(jù)包括漁獲量數(shù)據(jù)與標(biāo)準(zhǔn)化CPUE(Catch Per Unit Effort)數(shù)據(jù)。其中，標(biāo)準(zhǔn)化CPUE數(shù)據(jù)分別來自中國臺灣與日本延繩釣漁業(yè)。漁獲量數(shù)據(jù)的時間跨度為1950—2014年，中國臺灣標(biāo)準(zhǔn)化CPUE數(shù)據(jù)的時間跨度為1980—2012年，而日本標(biāo)準(zhǔn)化CPUE數(shù)據(jù)的時間跨度為1963—2014年[1]。由于1972年以前的日本標(biāo)準(zhǔn)化CPUE數(shù)據(jù)存在問題[1,7]，因此，本文僅使用了1972年及以后的數(shù)據(jù)。上述數(shù)據(jù)均來自印度洋金槍魚委員會(Indian Ocean Tuna Commission，IOTC)網(wǎng)站(http://www.iotc.org/meetings/17th-working-party-tropical-tunas-wptt17)。

1.2 評估模型及配置

1.2.1 評估模型 本文使用貝葉斯剩余產(chǎn)量模型對印度洋黃鰭金槍魚進(jìn)行資源評估，該模型由軟件JABBA(Just Another Bayesian Biomass Assessment)實現(xiàn)參數(shù)估計[11]。

1.2.2 模型參數(shù)的先驗設(shè)置

(1)K，φ，qJ及qT的先驗設(shè)置 根據(jù)文獻(xiàn)[1]，假設(shè)K服從1.06×106～1.69×107t的均勻分布，記作U(1.06×106, 1.69×107)，φ服從中值為0.90，變異系數(shù)為0.001的對數(shù)正態(tài)分布，記作L(0.9, 0.001)；假設(shè)qJ與qT均服從0至100的均勻分布[11]，記作U(0, 100)。

(3)r的先驗設(shè)置 根據(jù)文獻(xiàn)[1]，假設(shè)r服從中值與變異系數(shù)分別為 0.46與0.22及0.75與0.15的對數(shù)正態(tài)分布，分別記作L(0.46,0.22)、L(0.75,0.15)。

根據(jù)m及r先驗設(shè)置的不同，共有10個模型(見表1)。

1.2.4 模型收斂診斷 模型收斂診斷采用Gelman-Rubin 統(tǒng)計量，并以1.1為閾值，即Gelman-Rubin 統(tǒng)計量小于1.1，則認(rèn)為模型收斂[13]。本文僅對收斂模型進(jìn)行分析。

表1 不同模型的先驗設(shè)置

注：*表示在實際計算中取1.001，因為1.0會使模型分母為0[11]；其他參數(shù)的先驗分布在各模型中均相同，因此表1未列出。

Note:*denote the 1.001was used in practical calculations, because the 1.0 make the denominator in the model become zero[11]; the priors for the other parameters are same for the 10 models and are not shown in Table 1.

模型擬合效果結(jié)合偏差信息準(zhǔn)則(Deviance Information Criterion, DIC)[14]與泰勒圖[15]進(jìn)行分析評價，其中DIC由式(1)計算。

(1)

D(s,θ)=-2ln(L(s|θ))。

(2)

式中：ln為自然對數(shù)；L為似然函數(shù)。

1.3 后向預(yù)報及預(yù)測能力評價

1.3.1 后向預(yù)報 為執(zhí)行后向預(yù)報，本文將1950—2005年期間的數(shù)據(jù)用于估計參數(shù)，2006—2014年的漁獲量數(shù)據(jù)用于后向預(yù)報、并預(yù)測日本及中國臺灣延繩釣漁業(yè)CPUE的值，依JABBA軟件[11]，其計算方程如下：

(3)

(4)

(5)

對式(5)多次迭代的結(jié)果取中值，則得到日本及中國臺灣延繩釣漁業(yè)CPUE的預(yù)測值。

(6)

(7)

(8)

當(dāng)Ei越接近0時，模型預(yù)測能力越好[10]。

(9)

為更好比較不同數(shù)據(jù)對觀測的預(yù)測效果，需去除不同觀測方差影響，因此，式(9)等式兩邊同時除以σi,O，得:

(10)

(11)

(12)

由式(10)可知，當(dāng)模型預(yù)測CPUE方差等于標(biāo)準(zhǔn)化CPUE方差(此時為1)，且相關(guān)系數(shù)為1時，E′i為0，模型預(yù)測能力最好。

1.4 管理風(fēng)險分析

在漁業(yè)管理中，通常根據(jù)捕撈控制規(guī)則及評估模型投影的結(jié)果確定最終TAC。為評估管理風(fēng)險，本文利用不同時段(即1950—2005年與1950—2014年)數(shù)據(jù)估計的參數(shù)及式(3)與式(4)，分別以2006年為起始年，進(jìn)行投影計算TAC，以比較不同時段數(shù)據(jù)、不同模型估計的TAC的差異。確定TAC的捕撈控制規(guī)則是：在2014年，當(dāng)年資源量(B2014)大于最大可持續(xù)產(chǎn)量下的資源量(BMSY)、并且當(dāng)年捕撈死亡系數(shù)(F2014)小于最大可持續(xù)產(chǎn)量下的捕撈死亡系數(shù)(FMSY)的概率大于G。G通常設(shè)為50%和60%。

2 結(jié)果

2.1 模型擬合效果

由圖1可知，模型預(yù)測的CPUE與日本標(biāo)準(zhǔn)化CPUE具有較高的相關(guān)系數(shù)(大于0.8)，與日本標(biāo)準(zhǔn)化CPUE所具有的方差(圖1為1.0)相比，模型預(yù)測CPUE的方差偏小(小于0.5)；模型預(yù)測的CPUE與中國臺灣標(biāo)準(zhǔn)化CPUE的相關(guān)系數(shù)較小(小于0.45)，但模型預(yù)測CPUE的方差偏大(大于0.9)。模型預(yù)測CPUE與日本標(biāo)準(zhǔn)化CPUE的去中心均方根誤差(即圖1實心點(diǎn)至空心點(diǎn)[1.0, 0] 的距離)小于0.8，而模型預(yù)測CPUE與中國臺灣標(biāo)準(zhǔn)化CPUE的去中心均方根誤差大于1.0，因此所有模型均對日本標(biāo)準(zhǔn)化CPUE擬合較好。

圖1 不同模型對CPUE擬合效果的泰勒圖

隨形狀參數(shù)值的增加，模型預(yù)測CPUE與日本標(biāo)準(zhǔn)化CPUE的相關(guān)系數(shù)變化較小，但模型預(yù)測CPUE的方差不斷減少，使其去中心均方根誤差逐漸變大。中國臺灣標(biāo)準(zhǔn)化CPUE與模型預(yù)測CPUE具有類似關(guān)系，但其去中心均方根誤差隨模型預(yù)測CPUE方差變小而逐漸變小(見圖1)。因此，隨m的增加，模型對日本標(biāo)準(zhǔn)化CPUE的擬合效果變差，而對中國臺灣標(biāo)準(zhǔn)化CPUE的擬合效果則逐漸變好。

若采用DIC選擇模型，則模型S6或S7應(yīng)為最佳模型(見表2)。

表2 不同模型的DIC估計的去中心均方根誤差、預(yù)測標(biāo)準(zhǔn)差與相關(guān)系數(shù)

注：*表示P<0.05。*denoteP<0.05.

2.2 模型預(yù)測能力分析

當(dāng)采用模型S1～S10進(jìn)行預(yù)測時，所有模型預(yù)測的CPUE均大于日本標(biāo)準(zhǔn)化CPUE，并且當(dāng)模型為S1、S2、S6和S7時，大部分日本標(biāo)準(zhǔn)化CPUE數(shù)據(jù)點(diǎn)不在模型預(yù)測CPUE的95%預(yù)測區(qū)間內(nèi)(見圖2)。模型預(yù)測CPUE與日本標(biāo)準(zhǔn)化CPUE均呈負(fù)相關(guān)關(guān)系，甚至是顯著的負(fù)相關(guān)關(guān)系(見表2)。模型預(yù)測CPUE與中國臺灣標(biāo)準(zhǔn)化CPUE關(guān)系略好(見表2、圖2)，均呈正相關(guān)關(guān)系，但不顯著(見表2)。同時，由表2可知，隨m的增加，模型預(yù)測能力有逐漸變好的趨勢，即去中心均方根誤差逐漸變小(見表2)。并且，當(dāng)r先驗設(shè)為L(0.46, 0.22)時(即模型S1～S5)，其模型預(yù)測效果好于r先驗設(shè)為L(0.75, 0.15)時(即模型S6～S10)的模型預(yù)測效果，模型S5的預(yù)測能力最好，但其DIC最大(見表2)。

(點(diǎn)為標(biāo)準(zhǔn)化CPUE，實線為模型預(yù)測CPUE，始于2006年的虛線表示95%預(yù)測區(qū)間。The points was observed CPUE； The solid lines was predicted CPUE； The dash lines began at 2006 denoted 95% prediction interval.)

圖2 觀測CPUE與模型預(yù)測CPUE

Fig.2 Observed CPUE versus predicted CPUE

2.3 模型估計的TAC比較

使用相同模型、不同時段數(shù)據(jù)估計的TAC存在較大差異，如當(dāng)使用模型S6及1950—2005年數(shù)據(jù)時，在F2014BMSY的概率大于50%和60%的情況下、其TAC分別為4.77×105和4.58×105t，而使用相同模型及1950—2014年數(shù)據(jù)時，該值分別為3.69×105和3.66×105t，TAC分別減少了23%與21%。使用1950—2005年數(shù)據(jù)時，采用DIC選擇的最佳模型為S6或S7，而在使用1950—2014年數(shù)據(jù)時，采用DIC選擇的最佳模型為S9或S10，四模型計算的TAC差異更大(見表3)，這使?jié)O業(yè)管理存在巨大風(fēng)險。

表3 不同模型估計的TAC

注：*模型分別使用了1950至2005年的漁獲量數(shù)據(jù)及1980—2005年、1972—2005年的中國臺灣與日本標(biāo)準(zhǔn)化CPUE數(shù)據(jù)；**模型分別使用了1950至2014年的漁獲量數(shù)據(jù)及1980至2014年、1972至2014年的中國臺灣與日本標(biāo)準(zhǔn)化CPUE數(shù)據(jù);***表示模型具有最小DIC。

Note：*denote the catch data from 1950 to 2005 and the standardized CPUE from 1980 to 2005 from China Taiwan or from 1972 to 2005 from Japan were used in the model;**denote the catch data from 1950 to 2014 and the standardized CPUE from 1980 to 2014 from China Taiwan or from 1972 to 2014 from Japan were used in the model;***denote the model with minimum DIC.

3 討論

3.1 標(biāo)準(zhǔn)化CPUE的選擇及先驗、形狀參數(shù)的設(shè)置

標(biāo)準(zhǔn)化CPUE的選擇及參數(shù)先驗的設(shè)置對資源評估結(jié)果有重要影響，但不是本文討論的重點(diǎn)，為此，本文僅參考了官文江等[1]的結(jié)果。本文沒有考慮r均勻無信息分布的情形，因為使用該先驗分布時，模型無法收斂。

在使用剩余產(chǎn)量模型時，由于形狀參數(shù)較難被估計，因此在漁業(yè)資源評估實踐中主要使用Fox模型(m為1)或Schaefer模型(m為2)。同時，由于Schaefer模型的產(chǎn)量函數(shù)具有對稱性，這使得其評估結(jié)果較難被接受，因此造成Fox模型更合理的假象，但事實可能并非如此。為此，本文增加了多種形狀參數(shù)的設(shè)置，以展示形狀參數(shù)的重要影響。本文結(jié)果表明，m對參數(shù)、生物參考點(diǎn)的估計有重要影響(見表3)，也影響模型的預(yù)測能力(見表2、圖2)。同時，使用1950—2005年數(shù)據(jù)與使用1950—2014年數(shù)據(jù)估計的最佳模型的m不同，可能暗示m不應(yīng)該作為常數(shù)處理。此外，Winker等[11]認(rèn)為m的值可由最大可持續(xù)產(chǎn)量下的產(chǎn)卵生物量(SSBMSY)與初始條件下的產(chǎn)卵生物量(SSB0)的比值確定。由于SSBMSY受選擇系數(shù)等具有時變特點(diǎn)的參數(shù)影響，從這個角度看，m也應(yīng)該具有時變特點(diǎn)。

3.2 后向預(yù)報方法可用于評價漁業(yè)資源評估模型的質(zhì)量

在當(dāng)前漁業(yè)資源評估中，模型診斷與選擇主要依賴于模型對觀測數(shù)據(jù)的擬合度[10]，如采用DIC或均方根誤差選擇模型，但模型對觀測數(shù)據(jù)擬合的好壞，并不一定能決定模型正確與否。如本文，當(dāng)使用1950—2005年數(shù)據(jù)時，所有模型均能較好擬合日本標(biāo)準(zhǔn)化CPUE數(shù)據(jù)，但模型預(yù)測的2006—2014年的CPUE均與相應(yīng)時段日本標(biāo)準(zhǔn)化CPUE呈負(fù)相關(guān)關(guān)系，模型預(yù)測能力較差。同樣，采用DIC選擇的最佳模型S6或S7，其預(yù)測能力也最差；并且當(dāng)使用1950—2014年數(shù)據(jù)時，采用DIC選擇的最佳模型為S9或S10，這表明隨著使用數(shù)據(jù)時段的不同，最佳模型缺少一致性(見表2、表3)。而模型預(yù)測能力較差、最佳模型缺少一致性等均表明模型可能沒有正確模擬種群的演化動態(tài)、漁業(yè)資源評估結(jié)果可能存在問題。

由于缺少對漁業(yè)資源量進(jìn)行直接觀測的有效手段，漁業(yè)資源評估與預(yù)測結(jié)果很難得到驗證，而后向預(yù)報方法為評價漁業(yè)資源評估模型的預(yù)測能力與評估質(zhì)量提供了有效手段，從而有利于漁業(yè)資源評估模型的正確診斷與選擇。如本文，利用后向預(yù)報方法，可評價不同模型的預(yù)測能力，評價DIC選擇模型的穩(wěn)定性，從而可在一定程度上判斷模型對種群演化動態(tài)模擬是否正確，模型評估結(jié)果是否存在問題。

3.3 后向預(yù)報方法可揭示評估結(jié)果的不確定性及其引起的漁業(yè)管理風(fēng)險

在當(dāng)前金槍魚漁業(yè)管理中，主要采用的捕撈控制規(guī)則是：經(jīng)過若干年 (如10年) 的模型投影后，采用的TAC必須使捕撈死亡系數(shù)小于FMSY，資源量(或產(chǎn)卵生物量)大于BMSY的概率大于一個確定數(shù)，這個數(shù)通常是50%或60%。就本文而言，當(dāng)使用1950—2005年數(shù)據(jù)時，DIC選擇的最佳模型S6的TAC分別為4.77×105t (50%)和4.58×105t (60%)，而使用1950—2014年數(shù)據(jù)時，最佳模型S9的TAC分別為3.22×105t (50%)和3.18×105t (60%)，這兩者估計的TAC相差1.4×105t以上，這將使?jié)O業(yè)管理存在極大風(fēng)險。因此，通過后向預(yù)報方法不僅可發(fā)現(xiàn)DIC選擇的最佳模型不能提供最佳預(yù)測、DIC選擇的最佳模型不能保持時間上的一致性，并可揭示這些模型提供的TAC存在巨大差異，這將使?jié)O業(yè)管理者更加了解漁業(yè)資源評估的不確定性及其可能引起的漁業(yè)管理風(fēng)險。當(dāng)需要利用這些剩余產(chǎn)量模型為印度洋黃鰭金槍魚的漁業(yè)管理提供TAC時，后向預(yù)報方法提供的結(jié)果將提醒漁業(yè)管理者必須考慮這些不確定性。如本文，漁業(yè)管理者應(yīng)將F2014BMSY的概率提高，至少應(yīng)大于60%以上，以增加漁業(yè)管理的穩(wěn)健性、確保TAC能使?jié)O業(yè)管理目標(biāo)得以實現(xiàn)。