不同模態下氣候變化對智利竹筴魚補充量的影響

張暢，李綱，陳新軍,2,3,4,5*

( 1.上海海洋大學 海洋科學學院，上海 201306；2.青島海洋科學與技術試點國家實驗室 海洋漁業科學與食物產出過程功能實驗室，山東 青島 266237；3.上海海洋大學 大洋漁業資源可持續開發教育部重點實驗室，上海 201306；4.上海海洋大學 國家遠洋漁業工程技術研究中心，上海 201306；5.上海海洋大學 農業農村部大洋漁業開發重點實驗室，上海 201306)

1 引言

智利竹筴魚（Trachurus murphyi）是東南太平洋重要的經濟種類，于20世紀70年代開始大規模開發和利用，1995年達到歷年最高產量，為496萬t，其后多年漁業資源出現衰退，2013年產量僅35萬t，為近30 a以來的最低產量，2019年捕撈產量又恢復到63.8萬t[1]。作為沿海國家和遠洋漁業國家的重要捕撈對象之一，其資源評估與管理工作受到世界各國和南太平洋區域漁業管理組織的極大關注，補充量作為研究種群資源量變動的一個重要因素，也是資源評估管理中不可忽略的重要因素。

通常，漁業資源的補充量除了取決于資源本身的再生產機制，還受到生態環境的制約[2]。在目前已有研究中，主要通過兩個方面考慮環境對補充量的影響：一是使用一個或多個已知環境變量，構建短期的補充量預測模型[3]；二是考慮長期環境氣候變化對補充量的影響，即考慮模態變動[4]。自然的模態變動可能是劇烈且突然的，與海洋氣候聯系密切，它有別與捕撈因素的影響，在漁業管理中有著十分重要的意義[5]。近年來，基于貝葉斯模型平均法構建的親體-補充量模型在漁業中得到較好應用[6-7]，該方法不僅能分析各環境因子對資源補充量是否具有長期、穩健的解釋能力，還能用于補充量預測。為此，本研究擬通過貝葉斯模型平均法，分析在不同模態中的環境因子對智利竹筴魚的親體補充量影響及其內在變化規律，為其資源評估和科學管理提供參考。

2 材料與方法

2.1 漁業與環境數據

（1）漁業數據。1971-2017年智利竹筴魚親體量與補充量源自南太平洋區域漁業管理組織（South Pacific Reginal Fisheries Management Organization, SPR FMO）評估報告[1,8]，以1齡魚的資源量作為補充量，以參與繁殖的群體資源量作為親體量。

（2）海洋環境數據?？紤]智利竹筴魚種群結構復雜，結合已有文獻研究[9-10]，選取目前已知的產卵期（11 月至翌年2 月）最大產卵場（35°～45°S，90°～75°W）（圖1）的海表面溫度(Sea Surface Temperature, SST)、海表面鹽度（Sea Surface Salinity, SSS）、海表面高度（Sea Surface Height, SSH）的平均值作為影響補充量的潛在影響因子，SST數據來源于https://oceanwatch.pifsc.noaa.gov/erddap/griddap/index.html?page=1&itemsPer-Page=1000，SSH與SSS數據來源于 http://apdrc.soest.hawaii.edu/las_ofes/v6/dataset?catitem=71。海洋尼諾指數 （Oceanic Ni?o Index, ONI） 以 Ni?o3.4 區 中 每 年11月至翌年2月的SSTA平均值來表示，時間范圍為1971-2017年，數據來源于https://ggweather.com/ONI/oni.htm。太平洋年代際振蕩(Pacific Decadal Oscillation,PDO）可直接造成太平洋及周邊地區氣候的年代際變化，對中上層魚類資源造成影響，為此，本研究選取每年11月至翌年2月的平均太平洋年代際振蕩，作為表征其資源量變化的環境因子，時間范圍為1971-2017年，數據來源于http://research.jisao.washington.edu/pdo/PDO.latest。

圖1 智利竹筴魚產卵場Fig.1 The spawning ground of Chilean jack mackerel

2.2 模態檢驗

基于連續t檢驗的模態變動分析(Sequential T-test Analysis of Regime Shifts, STARS），進行補充量模態分析[11]。當模態改變時，可通過模態躍變指數（Regime Shift Index, RSI）反映，其計算公式為

式中，diff為兩個連續且存在顯著差異的格局模態平均值的差值；j為新模態的起始年份；L為所選取的截斷長度；t為置信水平p下，2L-2自由度的t分布值；為方差；RSI為模態躍變指數；當模態向上轉變，則=xi-， 當模態向下轉變，則=-xi， 則=±diff；xi第i個變量；為L年內變量x的平均值；為新模態中變量x的平均值。根據相關文獻和參數計算設定標準[11-12]，設置顯著性檢驗置信水平p等于0.1，截斷長度L等于20，huber權重參數H等于2。

2.3 補充量模型構建

選用Ricker模型構建基于環境因子的親體-補充量模型[13]。其公式為

模型兩邊取對數后可轉化為線性形式：

式中，R/S為繁殖率，又稱為再生長率；α、β、β0為模型中的待定參數；R為補充量；S為親體量；X表示影響因子；n為影響因子的數量。本次研究中變量因子包括親體量、SST、SSS、SSH、ONI和 PDO，共 6個變量，可構建k=26個模型，所有模型構成模型集合M={M1，M2，···，Mk}。

2.4 貝葉斯模型平均法

2.4.1 模型參數估計

貝葉斯模型平均法（Bayesian Model Averaging，BMA）的基本原理為，根據已有的主觀信息指定先驗分布，解釋變量的后驗包含概率（Posterior Inclusion Probabilities，PIP）大小可作為選擇解釋變量的客觀標準，以模型后驗概率（Posterior Model Probabilities，PMP）為權重對2.3節中所有單項模型進行加權平均，形成包含所有信息的集合M，用于被解釋變量的預測[14]。

為了不對模型的優劣加入主觀色彩，模型先驗設為均勻先驗分布（即相同的模型先驗概率），參數先驗選擇單位信息先驗。利用馬爾科夫鏈蒙特卡洛（Markov Chain Monte Carb，MCMC）方法在模型空間中抽樣，計算迭代10 000次，拋棄前2 000次的計算結果。

模型后驗概率PMP表達式為

各因子的后驗包含概率PIP的表達式為

參數 β的后驗分布為

式中，y為ln(R/S)；X表示影響因子；為模型Mγ的先驗分布概率；為觀測數據的概率總和。

2.4.2 補充量預測

貝葉斯模型平均法不僅可用于解釋變量的選擇，還可用于預測被解釋變量。不同模態中，不同環境因子以不同的組合形式擬合多個親體補充量模型，以PMP為權重對所有可能的模型進行加權平均，通過BMA綜合多個模型進行集合預報，求得最后兩年的預報結果[17]。為便于考慮模態變動對親體-補充量的影響，將根據模態分析結果，設計多個基于BMA的親體-補充量預測模型用于對比分析。

式中，yt為t年ln(R/S)預測值；yt,i為t年第i個模型的ln(R/S)預測值；為第i個模型的后驗概率。

3 結果

3.1 模態檢驗

1971-2017 年，智利竹筴魚補充量發生了顯著的年際變化，大體上呈現先上升再降低，再上升降低的連續波動過程，具體表現為補充量在1971-1985年處于波動上升階段，且在1985年達到最高值，在隨后的幾年里其補充量持續下降，1995年后再次上升，于2000年達到第2個高峰值，隨后補充量下降，在2002-2015年維持在較低的范圍內（圖2）。STARS分析結果表明，1971-2017年間智利竹筴魚資源補充量可分為5個模態：第1模態為1971-1980年，該模態中補充量處于上升階段，補充量平均水平相對較低，多個國家和地區對該漁業進行了開發和利用；第2模態為1981-1990年，該模態中的補充量平均水平最高；第3模態為1991-2001年，該階段補充量明顯的起伏變化，但其平均值仍保持在較高的水平；第4模態為2002-2015年，該模態中的補充量平均水平最低；第5模態為2016-2017年，為新的模態階段，與前一個模態相比，該階段補充量水平有明顯的上升，但該模態時間序列短，故在后續分析中不予討論（圖2，圖3）。

圖2 1971-2017年智利竹筴魚補充量模態檢測Fig.2 Regime shift detection in recruitment for different values of Chilean jack mackerel during 1971 to 2017

圖3 1971-2017年智利竹筴魚補充量躍變指數Fig.3 Regime shift index values of recruitment for Chilean jack mackerel during 1971 to 2017

3.2 氣候變化特征

由圖4可知，在第1模態中，厄爾尼諾累計值持續下降。在第2模態中，正負序列（不同氣候事件）交替的頻次高，變動趨勢不穩定。在第3模態中，連續4 a為正數序列，厄爾尼諾累計值持續上升，經歷連續兩年減少，1997年迎來強厄爾尼諾，厄爾尼諾指數高達1.38，其后4 a為負數序列，累計值持續下降。在該模態中，厄爾尼諾指數累計值有明顯的上升下降過程，變動趨勢明顯。在第4模態中，2002-2006年間整體趨勢上升，2007-2014年整體趨勢下降，但連續增長或下降的過程沒有第3模態明顯。分析發現，發生模態變化的年份都為厄爾尼諾年（圖4）。

圖4 1971-2015年間厄爾尼諾指數的年際變化Fig.4 Annual changes of El Ni?o index during 1971 to 2015

由圖5可知，在第1模態中，PDO指數累計值持續下降。在第2模態中，1981-1988年PDO指數累計值持續上升，其后兩年持續下降。在第3模態中，1991年PDO指數驟降至-1.78，整體趨勢變動不明顯，但增減變化的頻率更高。在第4模態中，2003年PDO指數上升至1.86，其后小幅增長，2008-2014年整體趨勢下降，2015年又上升至2.25。PDO指數的累計值在3個模態中都經歷了增長再下降的過程，但在第2模態、第4模態中該過程更明顯、更規律，在第3模態中增減變化的頻率更高。分析發現，發生模態變動的年份前后，PDO指數都有明顯的變化（圖5）。

圖5 1971-2015年間PDO指數的年際變化Fig.5 Annual changes of PDO index during 1971 to 2015

3.3 不同模態中影響因子的篩選分析

在不考慮模態變動的親體-補充量關系中，親體量的PIP為99%，在環境因子中，SSH的PIP最高，為76%。第1模態中，親體量的PIP占34%，PIP最高的環境因子為SSS，其PIP值為46%；第2模態中，親體量的 PIP為 86%，El Ni?o的 PIP值為 49%，為該模態中PIP值最高的環境因子；第3模態中，親體量的PIP為100%，PIP最高的環境因子為PDO，其值為41%；第4模態中，親體量的PIP為98%，PIP最高的環境因子為 El Ni?o，其值為 33%（表1）。

表1 不同影響因子在不同模態階段的PIP（%）Table 1 PIP (%) of different impact factors on different regimes

3.4 預測模型

在不考慮模態變動的親體-補充模型中，6個因子各自組合形成多個模型，表2列舉PMP最高的前5個模型。模型1為基于SSH的親體-補充量模型，后驗概率最高，為36%。5個模型集合預報的累計后驗概率達到了70%。

表2 基于BMA的PMP最高的前5個模型的后驗概率（%）Table 2 Posterior distribution (%) of top five models with highest PMP of all models from BMA for no shift model

考慮模態變動對親體-補充量的影響，依次去除第1模態、第2模態、第3模態，設計3個基于BMA的親體補充量模型 (表3)：（1）方案 a，由第 2模態、第3模態、第4模態數據（1981-2015年）構建親體補充量模型；（2）方案 b，基于第 3模態、第 4模態（1991-2015年）構建親體補充量模型；（3）方案 c，由第 4模態（2002-2015年）構建親體補充量模型。3個方案中，方案a中最佳模型的后驗概率最高為24%，各個方案的累計后驗概率基本相同(表3)。

表3 不同方案中基于BMA的PMP最高的前5個模型的后驗概率（%）Table 3 Posterior distribution (%) of top five models with highest PMP of all models from BMA for different program

由圖6可知，從預測結果看，BMA可通過概率分布定量評價預報的不確定性（圖6）。不區分模態的預測結果最接近真實值，優于方案a、b、c，方案c最為保守，預測值明顯小于其他方案。但其他方案的預測概率約為60%～80%，方案c的預測概率大于80%，優于其他方案。

4 分析與討論

4.1 氣候變化與補充量模態變動之間的聯系

南太平洋智利竹筴魚分布極其廣泛，東起智利、秘魯和厄瓜多爾沿岸，西至新西蘭、澳大利亞沿海，在 30°～60°S、78°～177°W 廣闊海域均有密集的智利竹筴魚群體，該區域也稱為“竹筴魚帶”[18]。該海域大尺度環流、中尺度海流以及各種海洋環境要素的變化活躍，智利竹筴魚的補充量受到環境要素、氣候事件的影響明顯[19-20]。變動的生態環境也可能會引起竹筴魚補充量的模態變動，從而給資源評估和管理帶來不可忽略的影響[21]。

本研究基于STARS對智利竹筴魚補充量進行模態分析，結果表明其存在多個模態，且在不同模態中，占主導地位的環境氣候因子并不相同。在第1模態中，親體量對補充量變動的解釋概率明顯偏低，各影響因子的解釋概率也明顯存在異常，可能的原因是該階段漁業資源變化更多的受到捕撈因素的影響[22]，并不適用于分析補充量與環境的關系。其余模態與傳統環境因子分析（不區分模態變動）相比，在環境因子的解釋概率排序上有明顯區別。在傳統環境因子分析中，SSH占有極高的解釋概率，而考慮模態變動的分析結果認為，各環境因子在不同的時期存在不同的排序。從統計結果分析，SSH的PIP明顯偏大且優于其他環境因子，極少有環境因子的解釋概率超過60%[23]，從目前已有的智利竹筴魚與環境關系的研究結果分析，也更偏向于區分模態的分析結果[24-25]。傳統環境因子分析往往將長時間序列的資源變動歸因于一個關鍵因子上，而忽略了資源變動受到多個環境因子的共同影響，在不同時期占主導地位的環境因子也會發生改變。要深入了解漁業資源對環境氣候的響應以及魚類動態變動及生態系統調節機制等內容，模態分析是不可忽略的。

在不同環境時期，智利竹筴魚的補充量水平會發生改變，基于這些不同時期的平均補充量水平來擬合親體-補充量關系和制定管理方案，可能不利于漁業的可持續發展[26]。在智利竹筴魚的評估報告[1]中，選取兩個基于不同時間段親體-補充量關系（1970-2015年、2000-2015年），并分別放入漁業資源評估模型，其結果存在較大的區別。研究認為，基于2000-2015年數據的分析結果更為保守，可用于最終的資源評估與管理。此外，SPRFMO評估報告[1]中選取的時間段（2000-2015年）與本研究中的第4模態相近，基于第4模態的預測結果也明顯比其余方案更保守（圖6）。但從預測結果看，4種預測方案中，不考慮模態變動的方案最接近真實值，基于第4模態的結果偏差最大，其原因可能是用于預測的2016年和2017年補充量水平高于第4模態補充量水平，與2015年以前有較大的差異（圖2），因此只基于第4模態（低水平）進行預測其結果自然會偏小。目前，在SPRFMO智利竹筴魚評估報告中，數據段的選取更多的依賴于人為經驗，但本文基于模態檢驗的科學方法將智利竹筴魚劃分為多個模態進行分析，能為數據段的選取及資源評估提供了科學依據。

圖6 基于BMA的ln(R/S)預測概率Fig.6 Predictive probabilities for ln(R/S) based on BMA

4.2 厄爾尼諾和PDO事件對智利竹筴魚補充量的影響

在第2模態、第4模態中，厄爾尼諾皆為解釋概率最高的環境因子，其指數的正負序列（不同氣候事件）交替頻次高，持續處于某一氣候事件的時間短，與第3模態中的表現有明顯的區別。這可能是由于不同氣候事件高頻的交替發生，這種不穩定的氣候環境對智利竹筴魚資源補充量產生了更大的影響。研究認為，智利竹筴魚在面對不同氣候事件時，存在不同的應對機制。以厄爾尼諾和拉尼娜為例，在厄爾尼諾期間，海表面溫度、海表面高度、風速、鹽度等眾多環境要素都會發生改變，智利竹筴魚仔稚魚分布、洄游路徑等多方面產生變化[27]。同時，由于捕食者減少、不同性質海水擴展等原因，智利竹筴魚棲息范圍會有所擴展[9]。Arcos等[28]研究認為，厄爾尼諾對智利竹筴魚仔稚魚肥育場的影響往往會持續3～4 a的時間。在拉尼娜期間，氣候條件與厄爾尼諾存在顯著不同，海水的降溫可能導致智利附近海域的幼魚進入厄瓜多爾和智利的專屬經濟區[29]。

PDO是一種反映從年際到年代際的氣候變率強信號，是太平洋重要的氣候系統之一，對太平洋乃至全球的大氣環流有著重要影響，也是氣候預測的關鍵因子[30-31]。除了對各地區氣候的影響，PDO對海洋魚類種群的興衰也起著調節作用[32]。PDO與許多中上層魚類的補充量有著直接的聯系，PDO正負相位時期、冷暖時期對應的補充量都有明顯的變化[33]。Dioses等[34]研究認為，年際間PDO變化與智利竹筴魚的緯度隔離有一定關系。分析認為，1977-1999年為PDO暖期，1999-2015年為PDO冷期[35]。冷暖期交替發生在本研究中的第3模態時期，這也解釋了該模態中PDO年際間變動更為不穩定，而不穩定的年際間PDO變化對智利竹筴魚補充量產生了更大的影響。

本研究認為，在發生模態變動年份前后，均發生了厄爾尼諾事件且PDO指數有明顯的躍變，這表明智利竹筴魚補充量與厄爾尼諾和PDO氣候事件存在緊密聯系，PDO年代際冷暖期交替與厄爾尼諾現象等可能是誘發智利竹筴魚補充量發生模態轉變的重要因素。Espino[36]研究也表明，PDO、厄爾尼諾、智利竹筴魚三者間的變動存在緊密聯系。一方面，不同尺度間的氣候變化存在交互作用，年代際PDO變化會影響厄爾尼諾現象，而物理環境中這些低頻的變化可能會導致海洋生態系統的改變[37-38]，從而導致智利竹筴魚補充量的模態改變。另一方面，PDO和厄爾尼諾會以不同的方式和不同的量級對漁業資源造成影響[14]。在今后的研究中，需要結合物理海洋、氣候學等多學科，加強了解智利竹筴魚補充量對氣候環境變化的響應機制和機理。

4.3 其他環境因素對智利竹筴魚補充量的影響

海表面高度是用于檢測和表征中尺度渦流的重要物理量[39-40]。渦流系統對智利海域葉綠素的聚集起著重要作用，冬季高濃度葉綠素海域多數是因此形成[41]。同時，渦流和上升流之間的耦合使得高初級生產力區域擴張，持續時間延長，浮游生物隨之聚集，為智利竹筴魚仔魚提供連續的食物來源，有利于其生長[42-43]。此外，渦流的垂直影響深度完全覆蓋仔魚生活水層，部分仔魚被夾帶入渦流內會被運到外海，在適當條件下可能會促進產生新的仔、稚魚棲息地[44-45]。在早期生活史階段，渦流等中尺度海流變化過程對智利竹筴魚空間分布、生長以及資源補充量等有著不可忽略的影響。

在本研究中，溫鹽的PIP基本維持在20%～30%，這表明兩者對智利竹筴魚補充量也具有一定的影響。東南太平洋海域擁有眾多不同溫鹽特質的水團，占主導地位的水團改變及強弱改變都會導致竹筴魚分布情況的改變[46]。當高溫（溫度大于18.5℃）、高鹽（鹽度大于34.9）亞熱帶表層水團南移并覆蓋近岸海域時，魚群分布更為集中，漁獲量增加；當低溫（約12℃）、低鹽（約34.25）的亞南極中層水團強于赤道表層水時，智利竹筴魚分布在更深的水層中，漁獲量明顯減少。秘魯海域的智利竹筴魚常聚集在表溫為11～28℃的水域，智利和公海海域智利竹筴魚聚集在表溫為9～20℃的水域, 當海表溫上升時，智利竹筴魚適宜棲息地有明顯的南移趨勢。與其他魚類不同，鹽度對智利竹筴魚的影響更為直接[47]，張敏等[29]研究認為，較穩定的鹽度將有利于智利竹筴魚集群的穩定。

模態檢驗方法的參數設置不僅需要遵循模型使用的前提條件，同時還需考慮智利竹筴魚漁業資源本身變動規律及其特性，任意降低截斷時間長度可能會出現更多的模態[48]，時間序列過短的模態不利于漁業資源評估和管理，選擇一個恰當合理的時間序列長度在不同模態的漁業資源評估是極為重要的。氣候環境變化對智利竹筴魚早期生活史階段的影響是不可忽略的，合理的劃分模態，有助于探究不同階段的環境與補充量的影響模式。其次，利用模態變動的方法將長期氣候變化對補充量的影響納入考慮，有助于提高對魚類動態及其生態系統調節機制的理解，制定更合理的管理策略[49]。同時，貝葉斯模型平均法將補充量影響機制歸因于多個環境因子共同作用的結果，分析每個因子不同的解釋能力，還能以概率分布的形式描述補充量預報的準確性。模態分析與貝葉斯模型平均法的結合適用于多信息集合，可在今后的研究中對更多的環境氣候因素進行分析，促進物理海洋、氣候、生物等多學科交叉，建立起科學的智利竹筴魚補充量預報模型，為后續的漁業資源評估和管理提供基礎。