基于優化模型的萊州灣三疣梭子蟹(Portunus trituberculatus)棲息地適宜性研究*

張孝民 王秀霞 李少文 楊艷艷 徐炳慶 李 凡

(山東省海洋資源與環境研究院 山東省海洋生態修復重點實驗室 山東煙臺 264006)

棲息地適宜性是對棲息地特征與物種在該條件下生存質量的定量描述(易雨君等, 2013), 是當前海洋生態保護領域重要的研究內容之一, 對于生態養護、資源評估、種群動態研究都具有重要的理論和現實意義(鄒易陽等, 2016)。在人類活動不斷加劇和氣候變化的大背景下, 全球生物多樣性正處于持續的快速喪失中(徐煒等, 2016), 其中棲息地的減少已經成為物種和局部種群絕滅的重要原因(鄧文洪, 2009; Keinathet al, 2017)。三疣梭子蟹(Portunus trituberculatus)生命周期短、繁殖力強、經濟價值高(林群等, 2015), 捕撈數量占全世界經濟蟹類的四分之一左右(Liuet al, 2013)。萊州灣位于渤海南部, 是赫赫有名的“萊州梭子蟹”的產地, 更有黃河和小清河等眾多河流入海, 灣內基礎餌料豐富, 水溫鹽度適宜, 是包括蟹類在內的眾多經濟性海洋生物的產卵和棲息場所。近幾十年來, 受環境污染、過度捕撈等人類活動的影響, 萊州灣海洋生境遭到破壞, 眾多生物的棲息地出現破碎化現象, 部分海區呈現荒漠化趨勢。與其他海洋生物一樣, 三疣梭子蟹的棲息地也在不斷減少, 資源量持續衰退(盧曉等, 2018)。棲息地適宜性評價在內陸河流及湖泊生物資源中運用較為廣泛, 在海洋中的運用起步較晚(龔彩霞等, 2011), 近年來研究集中在大洋魚類生物資源, 對于近海生物資源有較大的應用空間, 國內有關三疣梭子蟹棲息地適宜性的研究尚未見報道, 亟待開展相關研究。

HSI 模型是進行棲息地適宜性定量評價的經典方法, 該方法最早由美國地理調查局國家濕地研究中心魚類與野生生物署于20 世紀80 年代初提出(U.S.Fish and Wildlife Service, 1981)。GAM 模型常用來分析生物豐度和環境變量的非線性關系(Woodet al,2002), 可以在HSI 模型中選擇環境變量。BRT 模型基于決策樹理論的學習方法, 可有效計算出變量因子對模型建立的貢獻大小(De’ath, 2007), 生成多重回歸樹, 獲取變量的權重分布。

本文構建了四種HSI 模型, 包括普通HSI 模型;GAM 優化HSI 模型; BRT 優化HSI 模型; GAM 和BRT優化HSI 模型, 對萊州灣三疣梭子蟹棲息地適宜性進行研究, 以期摸清資源分布狀況, 明確適宜的環境特征, 全面了解其棲息地生存質量, 為三疣梭子蟹棲息地保護和萊州灣生態修復提供幫助。

1 材料與方法

1.1 數據來源和預處理

三疣梭子蟹數據來源于2010～2020 年夏季底拖網調查, 調查區域為萊州灣海域(119°00′～120°15′E,37°12′～37°40′N), 調查站位每年20 個(圖1)。調查船只功率260 kW, 調查網具為單船底拖網, 網口周長30.6 m, 囊網網目20 mm, 拖曳時網口寬度約8 m。

圖1 萊州灣三疣梭子蟹調查站位Fig.1 Sampling stations of P. trituberculatus in the Laizhou Bay

環境數據利用美國 ORI-ON520M-01A 便攜式YSI 水質分析儀現場測定, 內容為水深(depth,D)、底層水溫(bottom sea temperature, BST)、底層葉綠素(bottom sea chlorophyll, BSC)、底層鹽度(bottom sea salinity, BSS)。大型底棲生物(Macrozoobenthos, Mzb)取樣使用開口面積為0.05 m2箱式采泥器采集, 每個站位取2 個平行樣, 使用0.05 mm 的網篩分選大型底棲動物, 樣品的保存、處理、計數和稱量等均按《海洋調查規范》GB/T 12763.6—2007 (中華人民共和國國家質量監督檢驗檢疫總局等, 2008)進行。

以調查網具拖速3 kn、拖網時間1 h 為基準, 根據各站位實際拖網時間對調查數據進行標準化處理,將其換算為單位時間的生物量(biomass, Bio) (g/h)和個數(individual, Ind) (ind./h), 以表征三疣梭子蟹的資源豐度。

1.2 研究方法

1.2.1 變量因子的選擇 GAM 是一種非參數化的多元線性回歸模型, 可以用來表示響應變量和解釋變量的非線性關系, 本文利用GAM 模型選擇HSI 模型中的變量。GAM 模型通常表示為:

式中,g為鏈接函數,μi=E(Yi),Xi為第i個響應變量的解釋變量,β為模型估計參數,Yi為第i個響應變量,fi為平滑函數。將變量因子逐漸加入GAM 模型中, 選取AIC 值最小的模型作為最佳模型(Damalaset al,2007)。AIC 值的計算(Reid, 2003)如下:

式中,m為模型中參數的個數,n為觀察值(數據樣本)的個數, RSS 為殘差平方和。

1.2.2 變量因子權重分析 BRT 模型基于變量因子的相對貢獻確定其在HSI 模型中的權重(Changet al, 2010; Yiet al, 2016), BRT 方法是在運算過程中隨機抽取一定量的數據, 分析自變量對因變量的影響程度, 而剩余數據用來對擬合結果進行檢驗, 對生成的多重回歸樹取平均值輸出, 從而得出自變量對因變量的相對重要性(Keinathet al, 2017)。本文使用R 4.0.3 中“GBM”包的“GBM.STEP”功能實現BRT 模型,其中抽樣率(Train. Fraction)設置為 0.8, 重復循環計算1 000 次, 獲取每個環境變量的相對貢獻, 即環境變量的權重分布。

1.2.3 HSI 模型 首先, 利用樣條平滑函數建立變量因子與適宜性指數(suitability index, SI)的關系模型,分別以生物量Bio 和個數Ind 建立SI, 調查數據中Bio和Ind 最小值為0, 三疣梭子蟹在調查海域資源密度不均, 最大值和最小值量級差距較大, 直接以Bio 和Ind 數據建模會降低SI 指數精度, 本文對環境變量進行分組, 以各組距范圍內Bio 和Ind 的總和建立每個組距下的SI。SI 的取值范圍為0～1, 一般認為0.7～1之間是適宜的棲息環境范圍(Brooks, 1997; Xueet al,2017), 計算公式如下:

式中, SIi為i環境變量的適宜性指數, SIi,Bio為以生物量Bio 為基礎建立的適宜性指數, SIi,Ind為以個數Ind為基礎建立的適宜性指數, Bioij為i環境變量的j組距的Bio, Bioi,max為i環境變量的最大Bio, Indij為i環境變量的j組距的Ind, Indi,max為i環境變量的最大Ind。

HSI 模型是由SI 模型構建的綜合模型, 取值范圍為0～1, 代表棲息地適宜性“差”到“好”, 本文采用算數平均法(arithmetic mean model, AMM)和幾何平均法(geometric mean model, GMM)來構建HSI 模型, 計算公式分別為:

1.2.4 HSI 模型的比較和驗證 采用2010～2018 年數據構建HSI 模型, 包括四種類型, 分別為: 普通HSI 模型; GAM 優化HSI 模型(GAM 對變量進行篩選后建立的HSI 模型); BRT 優化HSI 模型(BRT 對變量權重進行優化建立的HSI模型); GAM和BRT優化HSI模型。通過Pearson 檢驗2020 年實測值HSI 和模型預測值HSI 的相關性, 采用ArcGIS 10.2 普通Kriging 插值法繪制2020 年萊州灣三疣梭子蟹HSI 地圖。

2 結果

2.1 變量因子的選擇和權重分析

以生物量表征資源豐度的GAM 模型(表1), 依次選取了D、BST、Mzb 密度和BSS 組成AIC 值最小的模型, 累計解釋方差為49.3%,D和BST 與生物量極顯著相關(P＜0.01), Mzb 密度和BSS 與生物量顯著相關(P＜0.05)。以尾數表征資源豐度的GAM 模型, 依次選取了D、BST、BSS 和Mzb 密度組成AIC 值最小的模型, 累計解釋方差為43.5%,D、BST 和BSS與生物量極顯著相關(P＜0.01), Mzb 密度與生物量顯著相關(P＜0.05)。

表1 以生物量和尾數表征資源密度構建的GAM 模型Tab.1 P. trituberculatus abundance shown in the GAM model constructed by using biomass and mantissa

由圖2 可知, 以生物量和尾數表征資源密度構建的BRT 模型變量因子的相對貢獻順序相同, BSS 的相對貢獻最大, 分別為26.9% (生物量)和23.5% (尾數),BSC 的相對貢獻最小, 分別為 12.2% (生物量)和11.4% (尾數), 其他因子的相對貢獻大小依次為BST、D和Mzb 密度, 區間范圍為19.0%～23.3%。

圖2 以生物量和尾數表征資源密度構建的提升回歸樹模型Fig.2 P. trituberculatus abundance shown in the BRT model constructed in biomass and mantissa

2.2 變量因子與生物量和尾數的適宜性指數關系模型

從變量因子的適宜性指數曲線(圖3)可以看出,以生物量和尾數表征資源豐度擬合的三疣梭子蟹適宜棲息環境范圍(SI＞0.7)基本相同。BST 的適宜范圍為27.2～28.5 °C, BSS 的適宜范圍為25.5～30.5,D的 適宜范圍為 9.0～10.5 m, Mzb 密度的適宜范圍為2 900～3 100 個/m2, BSC 濃度的適宜范圍為 13.5～14.5 mg/m3。

圖3 以生物量和尾數表征資源豐度與環境因子的適宜性指數曲線Fig.3 Fitted suitability index (SI) curves based on the relationship between P. trituberculatus abundance in terms of biomass and mantissa and environmental factors

2.3 HSI 模型的比較和驗證

以生物量表征資源豐度的GAM 和BRT 優化算術平均模型, 與 2019～2020 實際調查的生物量數據Pearson 相關系數最高, 為0.653, 表明預測值與實測值有較強的相關性, 最優模型具有較強的預測能力。總體看來, 優化后的HSI 模型均好于普通HSI 模型,BRT 優化的模型略好于GAM 優化的模型, GAM 和BRT 優化的模型好于其他3 種模型(表2)。此外, 生物量表征資源豐度的模型好于尾數表征資源豐度的模型,算術平均模型的相關系數高于幾何平均模型。

2.4 萊州灣三疣梭子蟹棲息地空間分布

根據GAM 和BRT 優化的算術平均HSI 模型預測結果可以看出, 伏季休漁結束前三疣梭子蟹主要分布在萊州灣南部、東南部和東北部海域, 中部和西南部海域分布較少, HSI 值較低(圖4)。模型預測結果與2020 年實際調查的三疣梭子蟹分布情況大體相符,以生物量表征資源豐度的HSI 模型預測的三疣梭子蟹最適海域(HSI＞0.7)范圍明顯大于以尾數表征資源豐度的HSI 模型, 以生物量表征資源豐度的HSI 模型對資源豐度較高區域的預測更加準確。

圖4 2020 年萊州灣三疣梭子蟹棲息地優化模型預測結果與實際調查情況分布圖Fig.4 Spatial distribution of observed abundance of P. trituberculatus in 2020 overlaid with the predicted HSI values based on optimized HSI models

3 討論

3.1 HSI 模型的優化

生物的棲息地受多種因子的影響, 不同因子的影響大小不同, 合理選擇影響因子和設置各因子的權重是優化HSI 模型的兩種主要方法。本研究分別采用GAM 模型對三疣梭子蟹影響因子進行選擇和BRT模型設置因子權重, 結果(表2)顯示BRT 優化模型好于GAM 優化模型, 說明在對HSI 模型優化的過程中,權重設置比因子選擇更加重要, 與Zhang 等(2020)研究結果一致。通過查閱文獻資料基本可以囊括研究對象棲息地的大部分影響因子, 加上專家知識的判斷篩選重要因子, 可以比較準確的完成因子選擇的過程(龔彩霞等, 2011), 但權重設置受不同研究海域和不同研究對象的影響較大, 專家無法獲得足夠的信息對不同因子權重準確設置(Tianet al, 2009), 本研究利用BRT 模型設置因子權重, 與非優化模型相比更具科學性和針對性, 將因子影響程度大小定量表示, 這也是導致權重設置比因子選擇優化效果更好的原因。但仍不能忽視GAM 模型對因子選擇的重要性, 綜合利用GAM 模型和BRT 模型才能構建最佳的HSI 模型。

表2 不同HSI 模型預測值與2019～2020 年實際調查數據的Pearson 相關檢驗結果Tab.2 Results of the Pearson correlation in different HSI model-predicted values and the actual survey data in 2019～2020

研究對象資源豐度指標的選取也是影響HSI 模型建模效果的重要因素, 主要的資源豐度指標包括漁獲量(范秀梅等, 2020)、漁獲尾數、作業次數、單位捕撈努力量漁獲量(catch per unit effort, CPUE) (劉瑜等, 2021)、生物量和漁獲密度(朱承之等, 2020)等,相關研究大都采用單一指標描述資源豐度, 本文比較了生物量和尾數表征三疣梭子蟹資源豐度對構建HSI 模型的影響, 發現以生物量表征資源豐度的模型好于尾數, 在GAM 模型的累計解釋率也相對較高,可能的原因是三疣梭子蟹個體大小差距較大, 部分站位大型個體較多, 用尾數表征降低了該區域的資源豐度, 但在適宜性指數曲線擬合得到的適宜棲息環境范圍基本相同。范秀梅等(2020)分別基于漁獲量和作業次數表征日本鯖魚的資源豐度建立HSI 模型,發現兩種方法在預報漁場上精度相差不大, 原因是日本鯖漁船作業次數和漁獲量通常成正比關系, 作業次數反映了漁船的集中程度和漁民對產量的滿意度(Bordalo-Machado, 2006), 作業次數多表明該海域有很高的漁獲量。Tian 等(2009)比較了基于CPUE 和作業次數建立的西北太平洋柔魚HSI 模型, 結果表明,基于CPUE 的HSI 模型高估了柔魚的最佳棲息地范圍,基于作業次數的HSI 模型能更好地定義柔魚的棲息地范圍。因此, 在選取資源豐度指標時, 應充分考慮研究對象的生物學以及作業方式等特性, 并采用多種指標表征資源豐度, 以獲得更好的HSI 模型。此外,HSI 模型的計算方法中, 算術平均法AMM 普遍具有較高的效果(陳新軍等, 2009), 本文研究結果也證明了AMM 的優勢, 但相關研究表明GMM 對SI 的極大值和極小值較敏感(Vinagreet al, 2006)。

3.2 三疣梭子蟹棲息地的影響因素

萊州灣三疣梭子蟹冬季在渤海10～30 m 水深較深海域越冬, 4 月開始生殖洄游, 到3～5 m 深的港灣或河口等淺海附近海域產卵, 8 月大多數三疣梭子蟹度過產卵期, 進入索餌育肥期, 是三疣梭子蟹資源量最高的時期, 9 月伏季休漁結束后, 高強度的捕撈會使資源量迅速下降, 人為因素成為主導三疣梭子蟹棲息地分布的決定因素, 環境因素影響力下降(Tianet al, 2009), 因此8 月伏季休漁結束前是研究三疣梭子蟹棲息地適宜性的最佳時期。本文研究了影響萊州灣三疣梭子蟹的部分環境因子和生物因子, GAM 模型和BRT 模型結果均表明BST、BSS 和D對其棲息地影響較大, 與吳強等(2016a)研究結果相似。三疣梭子蟹是廣鹽性蟹類, 但其 BSS 的最適范圍為25.5～30.5 (圖3), 夏季雨水豐富, 萊州灣淡水入海量增加, 近海鹽度變化成為影響三疣梭子蟹棲息地的最主要因素(圖2)。三疣梭子蟹最適棲息地的BST 范圍為27.2～28.5 °C, 水溫不僅是決定三疣梭子蟹的季節性洄游的主要因素, 夏季是索餌育肥的重要時期,水溫過高會降低其體內消化酶的活性, 影響攝食行為(廖永巖等, 2008), 進而對棲息地分布產生影響。萊州灣大部分區域水深小于10 m, 最大水深18 m, 水深較淺的近岸水域是三疣梭子蟹產卵棲息的重要場所, 水深的變化會引起水溫、鹽度和營養鹽等的改變,對棲息地分布產生影響, 本文研究發現伏季休漁結束前三疣梭子蟹在萊州灣的最適水深為9.0～10.5 m。相關研究表明三疣梭子蟹分布與底棲生物的關系密切(吳強等, 2016b), 本文研究認為 Mzb 密度為 2 900～3 100 個/m2時是三疣梭子蟹最適范圍, 底棲生物是大型甲殼類的重要食物來源, 三疣梭子蟹與日本蟳、口蝦蛄和中國明對蝦之間存在較強的食物競爭關系(李凡等, 2021), 導致底棲生物密度成為限制其棲息地分布的影響因素, 夏季萊州灣餌料食物較為豐富, 底棲生物與水溫和鹽度的影響相比較小, 在食物缺乏季節或海域, 底棲生物的影響力將增強。

3.3 夏季萊州灣三疣梭子蟹棲息地空間分布

結合實際調查數據以及HSI 模型的預測結果, 夏季萊州灣整體表現為近岸海域比離岸海域更適合三疣梭子蟹棲息, 這與徐炳慶等(2021)的研究結果一致,8～9 月份伏季休漁結束前水溫較高, 三疣梭子蟹尚未進入越冬洄游階段, 集中分布于近岸海域。本研究調查船只為底拖網漁船, 無法在水深低于5 m 海區設置站位, 如采用小型船只或其他作業方式對淺海三疣梭子蟹資源進行調查, 將使近岸資源好于離岸資源的趨勢更加明顯。2020 年萊州灣西南部近岸海域在HSI 模型預測中棲息地適宜性較差, 主要原因為7211站位和7214 站位底層平均水溫30 °C, 底層平均鹽度23.8, 均處于適宜性指數曲線SI 較低的位置, 導致HSI 值較小, 成為不適宜三疣梭子蟹棲息的海域。但在2020 年的實際調查中, 7214 站位生物量較高, 為2 196 g/h, 可能的原因是站位設置中環境數據是與三疣梭子蟹資源數據同步, 近岸海域溫鹽環境變化較大, 站位設置稀疏使得環境數據的獲取出現偶然性,在以后的調查研究中, 應當對近岸海域環境監測站位設置更加緊密, 以獲取更加準確地數據, 提高棲息地預測的精確性。

4 結論

本研究運用GAM 模型和BRT 模型在建模過程中對HSI 模型進行了優化, 優化模型對適宜棲息地預測效果均好于未優化模型, 優化模型對萊州灣三疣梭子蟹棲息地適宜性有較強的預測能力。優化模型顯示BST、BSS 和D對三疣梭子蟹棲息地影響較大, 萊州灣南部、東南部和東北部是夏季三疣梭子蟹的適宜棲息海域。