白思琦，鄒曉榮, ，張 鵬，丁 鵬

（1.上海海洋大學(xué)海洋科學(xué)學(xué)院，上海 201306；2.國家遠(yuǎn)洋漁業(yè)工程技術(shù)研究中心/大洋漁業(yè)資源可持續(xù)開發(fā)省部共建教育部重點(diǎn)實(shí)驗(yàn)室/農(nóng)業(yè)農(nóng)村部大洋漁業(yè)資源環(huán)境科學(xué)觀測實(shí)驗(yàn)站，上海 201306）

智利竹?魚 (Trachurus murphyi) 是一種具有較高商業(yè)價(jià)值的高度洄游性經(jīng)濟(jì)魚種[1-4]，主要分布在東南太平洋的中上層海域，包括厄瓜多爾和秘魯?shù)妊睾蛞约爸抢麑俳?jīng)濟(jì)區(qū)向西延伸的部分外海水域[5-8]。其產(chǎn)量在相關(guān)沿岸國家及一些遠(yuǎn)洋漁業(yè)國家的捕撈產(chǎn)業(yè)中具有舉足輕重的地位[9-12]。智利竹?魚漁業(yè)自20 紀(jì)70 年代被大規(guī)模開發(fā)以來，產(chǎn)量多年占據(jù)世界前列[13]，1995 年達(dá)到高峰 (4.96×107t)。近20 年來，智利竹?魚資源量由于全球氣候長期變化和人類捕撈強(qiáng)度的增大而不斷降低，漁場的年際變動(dòng)愈加顯著[14]。

漁場資源變動(dòng)與海洋環(huán)境息息相關(guān)[15-21]，對(duì)此國內(nèi)外學(xué)者利用不同模型研究海洋環(huán)境因子對(duì)智利竹?魚漁場時(shí)空分布的影響，Laura 等[16]通過神經(jīng)網(wǎng)絡(luò)模型、汪金濤等[22]利用主成分和BP 神經(jīng)網(wǎng)絡(luò)方法、晉偉紅[23]利用偏最小二乘法計(jì)算模型、陳春光[24]和Li 等[25]通過建立不同權(quán)重的棲息地適應(yīng)性指數(shù)模型、楊香帥等[14]采用灰色關(guān)聯(lián)度分析方法研究了智利竹?魚漁場時(shí)空變動(dòng)與海洋環(huán)境因子間的關(guān)系，這些模型的全局性研究會(huì)忽略資源分布與環(huán)境因子的特殊區(qū)域。

海洋生態(tài)系統(tǒng)具有高度復(fù)雜性，魚類資源由于受到多種海洋環(huán)境因子的交互作用產(chǎn)生動(dòng)態(tài)響應(yīng)機(jī)制，出現(xiàn)空間上的變動(dòng)，且該變動(dòng)具有空間異質(zhì)性[26]。目前對(duì)智利竹?魚資源的研究均偏向于靜態(tài)，而地理權(quán)重回歸模型 (Geographically weighted regression,GWR) 將采樣點(diǎn)間的空間距離考慮在內(nèi)，對(duì)每個(gè)空間位置擬合一個(gè)局部回歸方程，得出各環(huán)境因子與單位捕撈努力量漁獲量 (CPUE) 關(guān)系的局部回歸系數(shù)，便于發(fā)現(xiàn)資源分布的“熱點(diǎn)區(qū)域”[27-29]。Windle 等[28]首次利用GWR 模型分析漁業(yè)數(shù)據(jù)，發(fā)現(xiàn)環(huán)境因子對(duì)漁業(yè)資源的空間分布有異質(zhì)性影響，目前該模型已廣泛應(yīng)用于陸上資源空間分布的研究，在漁業(yè)資源空間分布研究中也有較為廣闊的應(yīng)用前景。本研究運(yùn)用GWR 模型對(duì)影響智利竹?魚種群生物資源豐度的環(huán)境因子進(jìn)行對(duì)比分析，以期為我國智利竹?魚漁船的生產(chǎn)提供參考。

1 材料與方法

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

智利竹?魚捕撈活動(dòng)數(shù)據(jù)由中國遠(yuǎn)洋漁業(yè)協(xié)會(huì)大型加工拖網(wǎng)項(xiàng)目技術(shù)組提供，數(shù)據(jù)來源于大型拖網(wǎng)漁船“開富號(hào)”和“開裕號(hào)”，兩船為相同類型的中層拖網(wǎng)船，使用的漁具規(guī)格為1 632 ×444.00 m(277.40 m)，網(wǎng)板面積為12 m2。船長、型深、船寬、總噸位、主機(jī)功率分別為91.10 m、9.20 m、20 m、7 671 t、5 920 kW。捕撈日志記錄數(shù)據(jù)包括作業(yè)日期、作業(yè)水深、起放網(wǎng)時(shí)間、作業(yè)位置以及漁獲量等。作業(yè)區(qū)域?yàn)?3°W—92°W、19°S—46°S。由于每年9 月到翌年3 月為生產(chǎn)淡季，出于對(duì)數(shù)據(jù)充足性的考慮，時(shí)間范圍選擇為2012—2018 年的4—8 月。去除漁撈日志中捕撈數(shù)據(jù)記錄有誤的數(shù)據(jù)，有效數(shù)據(jù)為1 880 個(gè)。環(huán)境數(shù)據(jù)選取海面溫度 (Sea surface temperature,SST)、葉綠素a濃度 (Chlorophyllaconcentration,Chl-a)、海面溫度異常值 (Sea surface temperature anomaly,SSTA)和海面高度異常值 (Sea surface height anomaly,SSHA)及作業(yè)水深，海面溫度異常值是指實(shí)際海面溫度按照時(shí)間序列減去該時(shí)間序列的平均海面溫度所得到的海面溫度差值，海面高度異常值是指實(shí)際海面高度按照時(shí)間序列減去該時(shí)間序列的平均海面高度所得到的海面高度差值。SST、Chl-a、SSTA 和SSHA 來源于美國國家大氣局 (NOAA)Oceanwatch 網(wǎng)站 (www.oceanwatch.pifsc.noaa.gov/)，作業(yè)水深為實(shí)測數(shù)據(jù)，時(shí)間范圍為2012—2018 年4—8 月，空間分辨率為0.125°×0.125°，時(shí)間分辨率為周。

1.2 數(shù)據(jù)處理及分析

1.2.1 CPUE 計(jì)算 CPUE (t·h?1) 是用于表示智利竹?魚資源密度的主要指標(biāo)[14]，計(jì)算公式為：

1.2.2 漁場重心計(jì)算 分別統(tǒng)計(jì)各個(gè)年度每一經(jīng)度和緯度的漁獲量，對(duì)各周作業(yè)的漁場重心進(jìn)行計(jì)算，并分析比較漁場空間分布的變化，以0.5°×0.5°的空間分辨率按周進(jìn)行統(tǒng)計(jì)，作業(yè)漁場的經(jīng)緯度重心計(jì)算公式[30]為：

式中X表示漁場經(jīng)度重心，Xi表示第i漁區(qū)的經(jīng)度中心，Ci表示第i漁區(qū)的漁獲量，Ei表示第i漁區(qū)的捕撈努力量；Y表示漁場緯度重心，Yj表示第j漁區(qū)的緯度中心，Cj表示第j漁區(qū)的漁獲量，Ej表示第j漁區(qū)的捕撈努力量。

1.2.3 空間自相關(guān)分析 用于表示空間自相關(guān)性方法眾多，Moran's I 指數(shù)為國際中最廣泛使用的方法，該指數(shù)用于衡量空間要素的相互關(guān)系，數(shù)值范圍為 [?1,1]，大于0 表示呈正相關(guān)，反之則為負(fù)相關(guān)，經(jīng)檢驗(yàn)等于0 則表明不存在空間相關(guān)性[31-33]。該方法表達(dá)式為：

1.2.4 地理權(quán)重回歸模型 GWR 模型是將地理學(xué)第一定律深度融入局部空間統(tǒng)計(jì)方法的模型，通過對(duì)獨(dú)立抽樣的捕撈站點(diǎn)分別進(jìn)行模型的回歸分析計(jì)算，得到與空間位置相對(duì)應(yīng)的空間回歸系數(shù)，隨著空間位置不同而變化的參數(shù)估計(jì)量化表征空間關(guān)系異質(zhì)性特征[29]。若智利竹?魚CPUE 與環(huán)境因子之間的關(guān)系存在空間異質(zhì)性，則GWR 模型中變量的關(guān)系可能存在空間位置函數(shù)，因此采用反映局部空間分異的GWR 模型能在不同局部空間進(jìn)行參數(shù)估計(jì)[34]。GWR 模型是全局線性回歸模型的擴(kuò)展，根據(jù)每個(gè)空間位置周圍一定范圍內(nèi)的作業(yè)點(diǎn)，都可得到一個(gè)局部回歸模型的結(jié)果數(shù)據(jù)[27-28,35-36]。與傳統(tǒng)線性回歸模型的區(qū)別在于GWR 模型中用于擬合局部回歸模型各采樣點(diǎn)的權(quán)重值不是假定相等的，而是與距離回歸點(diǎn)的距離成反比，采樣點(diǎn)權(quán)重值的確定常常采用高斯和雙平方權(quán)重函數(shù)[35]。GWR模型公式[34]表示為：

式中yi表示智利竹?魚CPUE，βj(j=0、1、2、3、4、5) 表示回歸系數(shù)；(ui,vi) 表示地理空間中第i個(gè)回歸點(diǎn)的空間坐標(biāo)；ESSTi表示地理空間中第i個(gè)回歸點(diǎn)的海面溫度；EChlai表示地理空間中第i個(gè)回歸點(diǎn)的Chl-a濃度；ESSTAi表示地理空間中第i個(gè)回歸點(diǎn)的海面溫度異常值；ESSHAi表示地理空間中第i個(gè)回歸點(diǎn)的海面高度異常值；Edepthi表示地理空間中第i個(gè)回歸點(diǎn)的作業(yè)水深。

1.2.5 顯著性檢驗(yàn) 對(duì)地理權(quán)重回歸模型進(jìn)行顯著性檢驗(yàn)和模型誤差項(xiàng)空間自相關(guān)的顯著性檢驗(yàn)，常用的檢驗(yàn)方法為赤池信息準(zhǔn)則 (Akaike information criterion,AIC) 檢驗(yàn)[37]、蒙特卡羅檢驗(yàn)[38]和F檢驗(yàn)[39-40]。本文選取AIC 檢驗(yàn)，AIC檢驗(yàn)多用于比較不同模型的顯著性，當(dāng)兩個(gè)模型AIC 值相差超過3 時(shí)，較低AIC 值的模型具有更好的擬合效果。表達(dá)式為：

式中n是樣本數(shù)，σ是誤差項(xiàng)估計(jì)的標(biāo)準(zhǔn)離差，tr(S)是地理加權(quán)回歸模型的投影矩陣S 軌跡，是帶寬的函數(shù)。

本研究使用ArcGIS 10.6 軟件來進(jìn)行數(shù)據(jù)空間自相關(guān)性分析、GWR 模型建立及模型結(jié)果的顯著性檢驗(yàn)。

2 結(jié)果

2.1 漁場重心時(shí)空變化

2.1.1 CPUE 及漁場重心月際變化 智利竹?魚漁場2012—2018 年CPUE 月際變動(dòng)見圖1，智利竹?魚漁場2012—2018 年4—8 月CPUE 總體基本呈現(xiàn)為逐月下降的趨勢，僅2012、2013、2014年基本呈上升趨勢；8 月CPUE 波動(dòng)最小，2015 年數(shù)值最大 (10.64 t·h?1)，2017 年數(shù)值最小 (6.11 t·h?1)；4 月CPUE 波動(dòng)最大，2017 年數(shù)值最大 (14.08 t·h?1)，2013 年數(shù)值最小 (3.76 t·h?1)。智利竹?魚漁場2012—2018 年月際變動(dòng)見圖2，漁場逐月向北移動(dòng)，4 月漁場主要分布在78°W—87°W、40°S—45°S；5 月主要分布在78°W—90°W、40°S—45°S；6 月主要分布在78°W—88°W、38°S—44°S；7 月主要分布在76°W—90°W、33°S—41°S；8 月主要分布在73°W—88°W、19°S—39°S，且開始分為兩個(gè)區(qū)域，一部分逐漸向西北偏移，一部分向東北偏移；8 月漁場位置達(dá)到最北端，5 月達(dá)到最南端。

圖1 2012—2018年月平均單位捕撈努力量漁獲量變化圖Figure 1 Monthly CPUE change from 2012 to 2018

圖2 漁場重心月際變化Figure 2 Monthly gravity change in fishing ground

圖3 漁場重心年際變化Figure 3 Anually gravity change in fishing grounds

2.1.2 CPUE 及漁場重心年際變化 智利竹?魚漁場2012—2018 年CPUE 年際變動(dòng)見圖1，2012—2018 年年際CPUE 基本呈逐年上升的趨勢，2012 年CPUE 值最小(4.82 t·h?1)，2015 年CPUE 值最高(10.94 t·h?1)。智利竹?魚漁場2012—2018 年年際變動(dòng)見圖3，2012 年漁場主要分布在75°W—88°W、27°S—45°S；2013 年主要分布在77°W—81°W、36°S—45°S；2014 年主要分布在75°W—82°W、27°S—44°S；2015 年主要分布在77°W—90°W、34°S—45°S；2016 年主要分布在74—86°W、28°S—44°S；2017 年主要分布在73°W—84°W、20°S—45°S；2018 年主要分布在73°W—90°W、22°S—44°S；2017 年漁場位置達(dá)到最北端，2015 年達(dá)到最南端，2015 年漁場位置最為集中。

2.2 空間自相關(guān)性分析結(jié)果

Moran's I 指數(shù)結(jié)果顯示CPUE 與SST、Chl-a、SSHA、SSTA、作業(yè)水深 (Depth) 均具有顯著的正空間自相關(guān)性 (Moran's I 指數(shù)>0，表1)，且該變量空間分布格局較為清晰，可以進(jìn)行GWR 模型的建立。

表1 不同環(huán)境因子地理權(quán)重回歸模型回歸效果的比較Table 1 Comparison of regression effects of GWR models with different environmental factors

2.3 CPUE 與環(huán)境因子關(guān)系的空間異質(zhì)性

GWR 模型結(jié)果顯示智利竹?魚CPUE 與SST、Chl-a、SSTA、SSHA 和作業(yè)水深在空間分布上具有異質(zhì)性 (圖4)?；貧w系數(shù)值的大小表明GWR 模型中該環(huán)境因子對(duì)漁場時(shí)空分布影響的強(qiáng)弱，正值表明該環(huán)境因子對(duì)CPUE 存在正向影響，負(fù)值表明存在反向控制作用。SST、Chl-a、SSTA、SSHA、作業(yè)水深的局部回歸系數(shù)值都存在正負(fù)符號(hào)的改變，即不同區(qū)域環(huán)境因子對(duì)CPUE 的影響不同，數(shù)值的大小表明了這5 個(gè)環(huán)境因子在空間上均對(duì)CPUE 的影響顯著；整體來看，基于GWR 模型的CPUE 與SSTA 回歸系數(shù)比Chl-a、SSTA、SSHA、作業(yè)水深大，表明SSTA 對(duì)于東南太平洋智利竹?魚漁場分布的影響大于其他4 個(gè)因子；CPUE與SST、Chl-a、SSTA、SSHA、作業(yè)水深的回歸系數(shù)大部分為正值，即SST、Chl-a、SSTA、SSHA、作業(yè)水深對(duì)CPUE 均具有正向影響。

圖4 地理權(quán)重回歸模型評(píng)估的單位捕撈努力量漁獲量與海面溫度、葉綠素 a 濃度、海面溫度異常值、海面高度異常值和作業(yè)水深關(guān)系的局部回歸系數(shù)值Figure 4 Local regression coefficients of SST,Chl-a,SSTA,SSHA and gear depth by GWR model

自然環(huán)境條件具有顯著的空間差異，因此CPUE 與SST、Chl-a、SSTA、SSHA、作業(yè)水深也是呈現(xiàn)顯著的空間異質(zhì)性。擬合優(yōu)度 (R2) 和局部擬合優(yōu)度 (校正的R2) 越大則表示對(duì)變量的解釋率越大，校正的R2計(jì)算是按照分子和分母的自由度對(duì)它們進(jìn)行正規(guī)化，具有對(duì)模型中變量數(shù)進(jìn)行補(bǔ)償?shù)男Ч?。SSTA 基于GWR 模型回歸的R2值最大，為0.74，SSHA 的校正R2值最大，為0.59，SSTA 的校正的R2值次之，為0.58 (表1)。結(jié)果表明SSTA 對(duì)東南太平洋智利竹?魚漁場時(shí)空分布影響最大。

2.4 AIC 檢驗(yàn)結(jié)果比較

AIC 檢驗(yàn)見表1，結(jié)果顯示，5 種環(huán)境因子中SST 基于GWR 模型回歸的AIC 值最小 (1 022.08)，雖然其擬合優(yōu)度較低，但擬合精度最好。

3 討論

3.1 東南太平洋智利竹?魚漁場時(shí)空變動(dòng)

本研究發(fā)現(xiàn)2012—2018 年4—8 月智利竹?魚漁場分布在73°W—90°W、20°S—45°S，2017年漁場位置達(dá)到最北端，2015 年達(dá)到最南端且漁場位置最為集中，為77°W—90°W、34°S—45°S。魏廣恩等[41]研究認(rèn)為漁場位置可以用漁場重心來判斷。西風(fēng)漂流帶為45°S—60°S，鄒曉榮[42]研究認(rèn)為，西風(fēng)漂流帶為智利竹?魚南部的漁場屏障，本文研究結(jié)果與其一致。在漁場月際變化趨勢上，漁場逐月向北移動(dòng)，與秘魯寒流從南向北的流向有關(guān)，8 月漁場分為兩部分，一部分逐漸向西北偏移，一部分向東北偏移，該變動(dòng)可能與秘魯寒流在該區(qū)域分成兩部分，一部分流入南赤道暖流，一部分流向東澳大利亞暖流有關(guān)。2012—2018 年4—8 月CPUE 基本呈逐月下降的趨勢，8 月CPUE 最小，4 月CPUE 最高，表明在生產(chǎn)旺季中4 月為生產(chǎn)最高，可為今后東南太平洋智利竹?魚捕撈生產(chǎn)提供科學(xué)的管理依據(jù)。

21 世紀(jì)初，由于環(huán)境變化和捕撈導(dǎo)致資源驟減，2006 年成立了南太平洋區(qū)域漁業(yè)管理組織，限制各國的捕撈水平，旨在恢復(fù)智利竹?魚資源量，在時(shí)間年際變化趨勢上，智利竹?魚漁場2012—2018 年CPUE 基本呈逐年上升的趨勢，說明了南太平洋區(qū)域漁業(yè)管理組織總捕撈許可量制度管理初見成效。有研究認(rèn)為厄爾尼諾和拉尼娜現(xiàn)象會(huì)引發(fā)環(huán)境發(fā)生異常改變，進(jìn)而導(dǎo)致海洋環(huán)境生產(chǎn)力的變化，資源分布產(chǎn)生差異[14]，研究表明2015 年為厄爾尼諾年[12]，受其影響，CPUE 值最高，且探捕報(bào)告調(diào)查表明生產(chǎn)時(shí)間和產(chǎn)量均增加，可能與厄爾尼諾現(xiàn)象有關(guān)。

3.2 東南太平洋智利竹?魚CPUE 與環(huán)境因子關(guān)系

本文研究了2012—2018 年4—8 月環(huán)境因子對(duì)東南太平洋智利竹?魚漁場時(shí)空分布的異質(zhì)性影響，結(jié)果表明，環(huán)境因子中SST、Chl-a、SSHA、SSTA、作業(yè)水深存在顯著的時(shí)空差異性，表明了智利竹?魚漁場時(shí)空分布與該5 種環(huán)境因子存在異質(zhì)性關(guān)系，由于復(fù)雜的海洋環(huán)境以及較大的研究范圍與海洋環(huán)境因子的相互作用導(dǎo)致CPUE 與環(huán)境因子存在空間異質(zhì)性關(guān)系。結(jié)果表明SSTA 對(duì)于東南太平洋智利竹?魚漁場時(shí)空分布影響最大，說明厄爾尼諾和拉尼娜現(xiàn)象更能影響漁場的時(shí)空分布；CPUE 與SST、Chl-a、SSTA、SSHA、作業(yè)水深的回歸系數(shù)多為正值，即該5 種環(huán)境因子對(duì)CPUE具有正向影響，在78°W—85°W、35°S—45°S 正向關(guān)系尤其顯著，從整體來看，基于GWR 模型回歸系數(shù)的對(duì)比，海面溫度和SSTA 對(duì)漁場時(shí)空變動(dòng)的影響大于其他環(huán)境因子，相關(guān)空間格局的結(jié)果表明，漁場從南到北，自東向西移動(dòng)。牛明香等[9]基于GAM 模型將各個(gè)環(huán)境因子綜合起來定量評(píng)估其重要程度及影響，但GAM 模型得出的全局統(tǒng)一規(guī)律可能會(huì)忽略資源分布與環(huán)境因子的特殊區(qū)域，GWR 模型的使用，考慮了數(shù)據(jù)空間關(guān)系的局部特征，可研究漁場的動(dòng)態(tài)分布，便于發(fā)現(xiàn)智利竹?魚資源空間分布的“熱點(diǎn)”海域，目前GWR 模型多用于陸地資源的研究，未來其在海洋研究上具有廣闊的應(yīng)用前景，也能對(duì)未來漁業(yè)資源的開發(fā)和利用提供參考依據(jù)。

本研究僅分析了部分環(huán)境因子，還有許多環(huán)境因子例如海水鹽度、溶解氧、海流和食物來源等，均對(duì)CPUE 產(chǎn)生一定影響，之后可將其融入到模型中，增加模型精確度，得到更為客觀的分析結(jié)果；本文數(shù)據(jù)為我國大型拖網(wǎng)漁船的商業(yè)捕撈數(shù)據(jù)，缺少他國漁船捕撈數(shù)據(jù)，在今后的研究中需要對(duì)結(jié)果進(jìn)行進(jìn)一步驗(yàn)證；GWR 模型不適合預(yù)測調(diào)查海域以外的資源分布；在經(jīng)濟(jì)、自然資源管理、醫(yī)學(xué)等多個(gè)領(lǐng)域的研究趨向于平面型，在地理上僅部分研究可呈立體型，例如不同海拔具有不同溫度、降水、土壤環(huán)境等，與海洋環(huán)境相似，但部分海洋生物具有垂直移動(dòng)特性，使用該模型時(shí)需要考慮具體情況或改變相關(guān)影響因子。