劉 雅, 王 晶, 栗小東, 王迎賓
(浙江海洋大學水產學院, 浙江 舟山 316022)
三疣梭子蟹(Portunustrituberculatus)隸屬甲殼綱十足目梭子蟹科梭子蟹屬,為廣溫廣鹽性種類,是中國最重要的海洋經濟蟹類之一[1]。主要分布于日本、朝鮮、菲律賓、馬來群島、紅海及中國的渤海、黃海、東海和南海。在東海三疣梭子蟹主要分布在大沙漁場、長江口漁場和舟山漁場[2],其漁獲量約占全國產量的50%[3]。東海北部海域是東海最大的三疣梭子蟹產區,其漁獲量占東??倽O獲量的40%以上[4]。國外學者對三疣梭子蟹的相關研究較少,Ariyama H[5]基于1984—2008年日本大阪灣區域三疣梭子蟹的漁獲量數據,探究了過于頻繁的漁業開發和偶發的人類活動,以及氣候事件對三疣梭子蟹補充量造成的影響。我國有關三疣梭子蟹的研究主要集中于養殖領域,包括三疣梭子蟹的人工育苗[6]、池塘養殖[7]、常見病害防控[8]、配合飼料開發[9-10]以及優良品種選育[11]等。對野生三疣梭子蟹資源的相關研究相對有限,主要集中在資源分布、增殖容量和資源評估方面。如徐炳慶等[12]通過輪廓系數法對萊州灣三疣梭子蟹資源的時空分布特征進行了研究;栗小東等[13]基于機器學習方法,對東海北部海域三疣梭子蟹的資源分布進行了預測;林群等[14]構建了黃河口鄰近海域的Ecopath模型,初步評估了三疣梭子蟹的增殖容量;俞存根[4]對東海中北部重要的經濟蟹類及其現存資源量和最大可持續產量(Maximum sustainable yield, MSY)進行了分析;Wang等[15]基于傳統Schaefer模型對三疣梭子蟹的最大可持續產量及資源開發情況開展了評估。以上在對三疣梭子蟹進行資源評估的研究中并未考慮模型中相關參數的不確定性。在資源評估過程中,如果不考慮不確定性,將對漁業資源評估結論、漁業生產和生產管理造成負面影響。因此,實現漁業科學管理和合理開發的前提是對資源量的正確評估[16]。
本研究依據2001—2020年東海北部海域三疣梭子蟹的漁獲量數據和單位捕撈努力量漁獲量(Catch per unit effort, CPUE)數據,運用Schaefer剩余產量模型,將參數基準方案(參數的先驗分布均為均勻分布)與敏感性分析方案(參數與基準方案相同,但采用不同先驗分布)進行對比,對東海北部海域三疣梭子蟹的資源進行評估,并設定收獲率作為三疣梭子蟹資源管理策略且分析管理策略可能存在的風險。研究結果有助于掌握當前三疣梭子蟹的漁業資源變動規律,并基于目前三疣梭子蟹資源狀況對未來資源的動態變化趨勢進行預估。
東海北部海域三疣梭子蟹漁獲量數據來源于《中國漁業統計年鑒》[17],時間段為2001—2020年(見圖1)。除了資源量以外,其他一些因素也會導致CPUE的變動,因此需要對CPUE進行標準化,以剔除這些因素的影響。本研究使用了近20年的漁業生產數據,其中包含了3種三疣梭子蟹生產船型,同時考慮主要環境因子變化的影響,因此在使用廣義線性模型(Generalized linear model,GLM)對CPUE數據進行標準化時[18],考慮了年份、漁船類型和平均水溫3個因素的影響。具體模型如下:
lm(ln CPUE~factor(Year)+factor(Effort)+
factor(Vessel)+factor(Temp), data=RData)。
(1)
式中:lm是R語言中廣義線性模型的函數命令;ln CPUE代表對CPUE取對數;Year代表年份;Effort代表捕撈努力量;Vessel代表3種不同的船型,分別是刺網、蟹籠、單拖網;Temp代表底層溫度。
Schaefer非平衡剩余產量模型的離散時間形式如下:
(2)
(3)
式中:Bt為t年的生物量;r為內稟自然增長率;K為環境容納量;Ct-1為t-1年的漁獲量;q為可捕系數。其中r、K、q為模型的參數。It為三疣梭子蟹資源豐度的相對指數,εt為誤差項。
在進行資源評估的過程中,一般會假設剩余產量模型中初始年的生物量B0為K,通過減少模型估算的參數數量來達到降低模型評估難度的目的[19]。但這種假設成立的前提條件是能夠確定漁業開始時前一年的資源量[20]。由于中國對三疣梭子蟹的捕撈活動開始較早,因此將初始一年的資源量視為K并不合理。本研究的B0值參考Wang等[15]的研究結果。
假設CPUE正比于生物量,且觀測誤差服從對數正態分布,則似然函數的表達式為:
(4)
式中:L和ln均表示自然對數;I表示單位捕撈努力量漁獲量;θ表示隨機變量;Π表示連乘符號;q表示可捕系數;Bt表示t年的生物量;σ表示標準差。
由于本研究采用的漁獲量和CPUE數據的時間序列為20年,這對標準差(σ)的估算來說存在一定的不穩定性,因此在本研究將其視為恒定值(σ=0.2)[19]。
在貝葉斯理論中,參數的后驗分布會受到先驗分布設定的影響[21]。因此,要對先驗分布進行敏感性檢驗來驗證樣本數據提供的信息是否充足。具體方法是比較各參數在不同先驗分布下得到的后驗分布的差異大小,來判斷樣本數據提供的信息是否充足。若結果差異較大,則表明數據沒有提供足夠的信息[18]。
因為較難搜集到三疣梭子蟹相關參數的先驗信息,為了結果的合理性,故采用均勻分布作為模型參數先驗概率分布的基準方案。此方案K的基準先驗分布假設服從均勻分布(50×104t,550×104t),下邊界設置為50×104t,此值大于2001—2020年數據中的最大漁獲量。為了盡量減少上邊界值對后驗概率的影響,上邊界值被設定為550×104t[22]。對于可捕系數q,獲取其有信息的先驗分布具有較大難度,在研究中一般將這種無信息的先驗分布設置為均勻分布或對數均勻分布。為足以覆蓋內稟增長率r的合理范圍,r的上、下邊界區間一般包括與評估種類相似的物種r的估計值[23]。本研究為了驗證樣本數據是否提供了足夠的信息,用敏感性分析方案來替代基準方案中的參數先驗分布,比較兩種情況下得到的參數后驗分布差異大小。
根據三疣梭子蟹的相關研究[15],參數r、K、q基準方案和敏感性分析方案的先驗分布設定見表1。

表1 剩余產量模型參數r、K、q的先驗概率分布
貝葉斯方法不僅考慮了評估模型的不確定性,而且考慮了參數的不確定性,通過將經驗信息與似然函數結合在一起,得到模型參數的后驗分布,所以可以得到比傳統分析方法更可靠的結果。因此,貝葉斯方法被認為是目前國際漁業資源研究領域較有代表性的方法,且適用于缺乏某些特定數據的情況下研究漁業資源動態[24]。
本研究中模型參數r、K、q的后驗概率分布是通過MCMC(Markov chain monte carlo algorithm)算法測算的。此算法共迭代50 000次,將前5 000次的結果舍棄,后每隔40次儲存剩下的45 000次的運算結果。采用此算法計算時,各參數的初始值設定分別為:r=0.15,K=37 706.4,q=0.005[15]。
生物學參考點是從生物學角度設定的指標,包括目標參考點(Target reference points,TRP)和限制參考點(Limit reference point,LRP)。目標參考點是漁業管理為了持續獲得某一漁獲量所設定的指標,限制參考點指的是在漁業管理中為了避免出現種群的生物量過低導致種群不能被可持續利用的情況而設置的指標[25]。本研究涉及到的生物參考點有FMSY、BMSY、F0.1和MSY。FMSY和BMSY分別指漁業達到MSY水平時對應的捕撈死亡系數和生物量,F0.1表示平衡漁獲量和捕撈死亡系數關系曲線最大斜率的10%對應的捕撈死亡系數。通過計算可得:
(5)
F0.1=0.45r,
(6)
(7)
(8)
在本研究中,捕撈死亡系數的目標參考點Ftar設定為F0.1,限制參考點Flim設定為FMSY;BMSY設定為資源量的目標參考點Btar,BMSY/4設定為限制參考點Blim。可以通過以上參考點來判斷當前的漁業資源狀況。當捕撈死亡系數大于Flim時,表明該資源正在遭受過度捕撈(Overfishing),相反則表明該資源未遭受過度捕撈;當資源量小于Blim時,表明該資源已經遭到破壞,處于過度捕撈狀態(Overfished),相反則表明該資源未處于過度捕撈狀態。
本文采用不同的收獲率(Harvest rate,即捕撈死亡率,規定了漁業年捕撈漁獲量在資源量中的固定占比)來作為三疣梭子蟹資源管理策略。收獲率分別設定為0.01、0.02、0.03、0.04、0.05、0.06、0.07和0.08。因此,未來t年的漁獲量的計算公式如下:
Ct=hiBteε,ε?N(0,12)。
(9)
式中:Ct為t年的漁獲量;hi為收獲率;Bt為t年的生物量;ε為誤差項。
資源評估過程的復雜性會導致評估結果存在較多不確定性,因此通過建立指標來對不同的管理策略進行風險評估是十分必要的。本研究假設實施管理策略的時間段為2021—2030年。建立的指標有:
(1)管理結束時的生物量B2030,即2030年資源量的期望值。
(2)管理策略實施后2030年漁獲量的期望值C2030。
(3)2030年三疣梭子蟹生物量與BMSY之比的期望值(B2030/BMSY)。
(4)2030年三疣梭子蟹生物量的衰減率(Depletion),即實施某一管理策略后,2030年三疣梭子蟹生物量與K的比例(B2030/K)。
(5)2030年三疣梭子蟹生物量大于2020年三疣梭子蟹生物量的概率P(B2030>B2020),表示管理措施實施后資源較管理措施實施前資源恢復的概率。
(6)2030年三疣梭子蟹生物量大于BMSY的概率P(B2030>BMSY),表示資源恢復到健康水平的概率。
(7)2030年三疣梭子蟹生物量小于BMSY/4的概率P(B2030 在1.5中,我們得出剩余產量模型參數r、K、q的后驗概率分布后,可以通過反復的模擬計算得到指標概率分布,具體步驟見圖2。 圖2 指標概率分布計算流程圖 圖3為不同方案模型的后驗概率分布?;鶞史桨讣僭Or、K、q的先驗概率均為均勻分布,但由圖3的曲線可知r、K、q的后驗概率明顯不同于假設的分布類型。同時,由兩圖對比可知這兩種方案的參數后驗分布非常相似,這表明模型參數后驗分布的結果受本研究采用的漁業數據的影響遠大于先驗分布的影響,即數據提供了較豐富的信息。 圖3 基準方案和敏感性分析方案的后驗概率分布 模型估計的參數估計值見表2,結果表明兩方案下內稟增長率r的估計值分別為基準方案0.192,敏感性分析方案0.194,即后驗分布均值無明顯差異;基準方案和敏感性分析方案的環境最大容納量K的估計值分別為503.38×104t和504.23×104t;兩方案下可捕系數q的估計值相等。 表2 模型相關參數的估計量 基準方案和敏感性分析方案估算的三疣梭子蟹生物學和漁業學參考點見表3?;鶞史桨赶鹿浪愕腗SY為24.07×104t,MSY對應的生物量BMSY為251.70×104t,此時的捕撈死亡系數為0.096。此方案評估的F0.1為0.086。在敏感性分析方案下,估算的MSY為24.33×104t,BMSY為252.12×104t,此時的捕撈死亡系數為0.097。此方案評估的F0.1為0.087。由這兩種方案估算得到的結果可知,F2020的均值小于FMSY的均值且B2020的均值大于BMSY的均值,這表明當前東海北部海域三疣梭子蟹資源未遭受過度捕撈且未處在過度捕撈的狀態(見表3)。 表3 兩種方案下估算的生物學參考點 圖4的4個相位,紅色區域代表已過度捕撈且正遭受過度捕撈,而綠色區域代表未出現過度捕撈,黃色區域代表已過度捕撈或正在過度捕撈[26]。從圖4分析可知,2001—2020年,捕撈死亡系數基本低于FMSY且生物量大都在BMSY之上,這表明了近年來三疣梭子蟹資源狀況良好,基本不存在過度捕撈的狀況。 (藍線展示了B/BMSY和F/FMSY在2001—2020年的變化。The blue line shows the interannual variation of B/BMSY and F/FMSY between 2001 and 2020.) 基準方案和敏感性分析方案在不同的收獲率下的管理決策和風險分析的指標結果見表4。兩種方案的各項指標值均比較接近,在基準方案下,隨著收獲率的增大,2030年漁獲量期望值總體呈上升趨勢,但當收獲率從0.05增大到0.07時,2030年最大漁獲量期望值的變化不顯著。在敏感性分析方案下,收獲率從0.01增大到0.06時,2030年最大漁獲量期望值呈上升趨勢,但隨著收獲率的繼續增大,2030年最大漁獲量期望值卻逐漸減小。綜合兩種方案的結果可知,如果要獲取2021—2030年最大持續產量,應將收獲率控制在0.06左右,此時對應的資源量為23×104t左右。 表4 兩種方案下管理決策和風險分析的指標 當收獲率為0.06時,風險分析表明資源狀況并不十分樂觀。在基準方案下,B2030/BMSY的概率期望值為0.96,敏感性分析方案為0.99。對于B2030>BMSY,基準方案和敏感性分析方案的概率期望值分別為0.53和0.56。對于B2030 從漁業管理的角度出發,通過分析表4可知,當收獲率為0.04時,基準方案和敏感性分析方案中B2030>B2020的概率期望值均大于0.6,B2030>BMSY的概率期望值均大于0.8,且B2030 貝葉斯方法是利用統計數據(如漁獲量和CPUE等)與先驗概率分布相結合,從而推斷出后驗概率分布的統計學方法[27]。先驗概率分布可以依據以往的研究結果,也可以參考研究者的主觀判斷。由于先驗概率分布會對資源評估結果產生影響,選擇了不合理的先驗概率分布可能會導致錯誤的結果,所以對先驗概率分布的推斷目前仍有爭議[19]。因此,需要選擇正確的先驗概率分布,以避免后驗概率分布結果完全被先驗概率所主導[28]。 在本研究中,三疣梭子蟹資源評估方面的參考資料較少,信息不足,只能依據已有的研究結果設定模型參數K、q、r的先驗概率分布[15]。因此,為了減少先驗分布對后驗分布的影響,在基準方案中,我們選擇了范圍較大的均勻分布,同時設置了敏感性分析方案進行比較分析。從后驗分布來看,參數K、q、r與先驗的均勻分布相比差異很大,說明本研究采用的數據提供了充足的信息;敏感性分析方案下的后驗分布與基準方案下的后驗分布十分相似,說明資源評估的結果受到模型參數先驗概率分布的影響十分有限。 當前東海北部海域三疣梭子蟹資源狀況相對樂觀,在近年來并未遭受到過度捕撈且三疣梭子蟹的資源狀況良好。本研究得出的捕撈死亡系數F較小,可能是由于三疣梭子蟹適應力強、存活率較高且增殖放流數量大,使得三疣梭子蟹資源量維持在較高水平,因此估算出的捕撈死亡系數F較小。本研究得出的三疣梭子蟹的MSY約為24×104t,這一結果不同于Wang等[15]研究結果。Wang等[15]研究顯示,假設三疣梭子蟹的捕撈強度和放流數量保持在2015年的水平,三疣梭子蟹的MSY約為14.2×104t。造成結果不同的原因可能是采用了不同的評估方法或所用數據的時間跨度不同。Wang等[15]使用的是傳統的非平衡剩余產量模型,評估所使用的數據是2001—2015年間的漁業統計數據。當采用相同時間跨度的數據(2001—2015年)時,使用貝葉斯方法得到的MSY為19×104t,與Wang等[15]14×104t的研究結果差異明顯縮小。具體分析可得,2001—2015年的年產量均值約為10.4×104t,而采用2001—2020年的漁業統計數據算出的年產量均值約為12.1×104t。2015年的年產量約為18.5×104t,與本文評估得到的MSY值更接近。對于Wang等[15]未使用的時間段即2016—2020年,東海北部海域三疣梭子蟹的年產量均值約為17.0×104t,明顯大于2001—2015年的均值,且2001—2015年CPUE均值也偏小。對數據趨勢分析可得,2013—2020年的CPUE均值明顯大于2001—2012年,且2013年之后三疣梭子蟹年產量相較于之前明顯增加,2015年之后的產量雖略有下降,但總體仍保持較高水平。因為本文采用的數據時間跨度為20年,故后期研究海域三疣梭子蟹統計數據的變化可能會對模型產生較大影響,導致評估結果(如MSY)數據變大。由此可見,采用不同時間段的數據做出的評估結果可能會產生差異,影響主要來自漁獲量和捕撈努力量的年際變化趨勢。因此,全面且準確的長時間序列數據對于評估MSY及相關的生物學參考點至關重要。 漁業資源管理的主要目標是通過合理的調控與科學的管理,從漁業資源中獲得最佳的可持續利益。這就要求在對漁業管理策略進行選擇時,應從多角度考慮,要在保證漁業資源可持續開發利用狀態與獲得高漁獲量之間進行全面權衡。綜合兩種方案下的管理策略和風險指標分析結果可知,若要獲取最大產量,應將收獲率控制在0.06。但風險分析表明,此時的資源狀況并不樂觀。若將收獲率設定為0.06開展捕撈活動,三疣梭子蟹資源會出現過度捕撈現象,甚至有可能會出現資源崩潰的情況。所以,將收獲率設置為0.06對資源的可持續發展來說存在很大風險。 從資源保護的角度來看,將三疣梭子蟹資源的收獲率控制在0.03時,資源量會維持在較高的水平且不會出現資源崩潰的狀況。從經濟收益的角度來看,將收獲率控制在0.05時,基準方案和敏感性分析方案中2030年最大漁獲量期望值與結果中最大值非常接近,且出現資源崩潰的可能性較低。但在大多數漁業已經出現衰退趨勢的現狀下,我們應采取使資源狀況更安全的策略作為最佳管理策略。對利益及其存在的風險進行權衡后得出,0.04的收獲率在未來10年的管理期內是較合理且保守的管理策略。 本研究分基準和敏感性分析兩種方案,使用了與貝葉斯方法相結合的Schaefer剩余產量模型,對東海北部海域三疣梭子蟹的資源狀況進行了評估,并估算了三疣梭子蟹資源在2030年恢復到健康水平的概率和崩潰的概率。研究結果可為三疣梭子蟹資源養護和可持續利用提供有利支持,為制定科學合理的漁業管理策略提供參考。 但是,本研究仍然存在一些不足和需改善之處。第一,模型參數的不確定性,本研究假設了CPUE數據的標準差σ為0.2,該值的假設對評估結果造成的影響需要進一步研究與分析。第二,本研究采用的生產數據時間段為2001—2020年,因此模型參數的結果只能較客觀地反映這一時期的漁業狀態。由于三疣梭子蟹屬于易受環境影響的短生命周期物種,且漁業統計年鑒中對該種的統計也存在誤差,因此,基于2001—2020年間估算出的模型參數對其未來的漁業資源變化進行的預測會存在一定的不確定性。第三,三疣梭子蟹是東海北部海域重要的增殖放流種類之一,每年幼蟹放流數量達到幾千萬尾,而本研究并未將該影響納入模型評估當中,這無疑會對資源評估結果產生一定的影響。今后的研究中,我們將繼續積累更多的統計和調查數據,進一步完善評估模型,在評估過程中充分考慮環境因子和增殖放流等因素的影響,并嘗試使用生態系統模型對三疣梭子蟹的MSY開展評估,為合理開發、有效保護三疣梭子蟹資源提供科學參考。
2 結果
2.1 后驗概率分布和期望值


2.2 三疣梭子蟹資源狀況分析


2.3 管理策略分析

3 討論與分析
3.1 先驗和后驗概率分布
3.2 當前資源狀況評價及未來管理策略分析
3.3 研究的不足和展望