唐 政 單秀娟 關(guān)麗莎 金顯仕2, 胡芷君 韓青鵬,4

(1. 上海海洋大學(xué)海洋科學(xué)學(xué)院 上海 201306；2. 青島海洋科學(xué)與技術(shù)試點(diǎn)國(guó)家實(shí)驗(yàn)室海洋漁業(yè)科學(xué)與食物 產(chǎn)出過(guò)程功能實(shí)驗(yàn)室 青島 266071；3. 中國(guó)水產(chǎn)科學(xué)研究院黃海水產(chǎn)研究所 農(nóng)業(yè)農(nóng)村部海洋漁業(yè)可持續(xù)發(fā)展 重點(diǎn)實(shí)驗(yàn)室 山東省漁業(yè)資源與生態(tài)環(huán)境重點(diǎn)實(shí)驗(yàn)室 青島 266071；4. 中國(guó)海洋大學(xué)水產(chǎn)學(xué)院 青島 266003)

為了探究變量在不同空間尺度上的分布格局和變化規(guī)律，通常要求采樣信息具有空間上的連續(xù)性，但是受成本和技術(shù)等限制，對(duì)于高山、地層和海洋等難以到達(dá)或范圍較大的研究區(qū)域，只能獲取具有一定分辨率的分散樣本點(diǎn)數(shù)據(jù)(Li et al, 2014)，空間插值法可以利用有限的數(shù)據(jù)，對(duì)整個(gè)研究區(qū)域的變量進(jìn)行 預(yù)測(cè)，減少采樣調(diào)查的費(fèi)用和難度，因此，被廣泛應(yīng)用于水文、氣象、土壤、礦產(chǎn)、農(nóng)業(yè)和生態(tài)等方面 (孫銘帥等, 2017; 岳文澤等, 2005)?？臻g插值法可以分為非地統(tǒng)計(jì)法、地統(tǒng)計(jì)法和二者結(jié)合的方法(Li et al, 2014)。受自身算法、數(shù)據(jù)本身屬性以及樣本點(diǎn)數(shù)量等多種因素的影響，不同情況下空間插值法的精度存在一定差異(朱會(huì)義等, 2004)，例如，孫慧等(2017)研究發(fā)現(xiàn)，Cubist 方法預(yù)測(cè)的土壤Cd 含量結(jié)果精度較高；Chen 等(2016)發(fā)現(xiàn)，不同采樣時(shí)間對(duì)漁業(yè)資源密度空間插值的精度存在差異；Luo 等(2011)對(duì)多種方法估算的月降水量的空間分布進(jìn)行比較發(fā)現(xiàn)，協(xié)同克里格法(Ordinary cokriging, COK)的精度要高于泰森多邊形法、反距離加權(quán)法(Inverse distance weighted, IDW)和普通克里格法(Ordinary kriging, OK)。因此，在實(shí)際應(yīng)用時(shí)，有必要先對(duì)不同插值方法進(jìn)行比較分析，選擇最優(yōu)的方法。

物種多樣性在不同尺度上的空間分布特征是生態(tài)學(xué)研究中的重要問(wèn)題(李圣法, 2005; 馬斌等, 2008)。游泳動(dòng)物作為水生生物群落的重要組成部分，其多樣性空間變化可以在一定程度上反映物種多樣性的變動(dòng)趨勢(shì)。目前，只在陸地生物多樣性空間分布中發(fā)現(xiàn)不同空間插值法的比較研究(Hernandez- Stefanoni et al, 2006)，雖然空間插值法在海洋生物多樣性空間分布中已有應(yīng)用(Stelzenmüller et al, 2009；Tittensor et al, 2010; 張敏等, 2011; 呂振波等, 2012)，但是尚未見(jiàn)不同方法優(yōu)劣的判別研究。因此，本研究比較分析了5 種常用的空間插值法[包括IDW、徑向基函數(shù)法(Radial basis function, RBF)、局部多項(xiàng)式法(Local polynomial interpolation, LPI)、OK 和COK]對(duì)4 種多樣性指數(shù)[即物種種類數(shù)(S)、Shannon-Weaver多樣性指數(shù)(H’)、Pielou 均勻度指數(shù)(J)和Margalef豐富度指數(shù)(d)]的插值結(jié)果，選出最適合海洋游泳動(dòng)物多樣性插值的方法，以期為海洋生態(tài)學(xué)研究提供參考。

1 數(shù)據(jù)與方法

1.1 數(shù)據(jù)來(lái)源

數(shù)據(jù)來(lái)自2017年5月黃海和東海北部漁業(yè)資源底拖網(wǎng)調(diào)查，調(diào)查站位見(jiàn)圖1。調(diào)查船為中國(guó)水產(chǎn)科學(xué)研究院黃海水產(chǎn)研究所“北斗”號(hào)，平均拖速為3.0 kn，拖網(wǎng)時(shí)間為0.5～1 h，調(diào)查網(wǎng)具網(wǎng)目尺寸為836目× 20 cm，囊網(wǎng)網(wǎng)目為2.4 cm，網(wǎng)具總長(zhǎng)度為83.2 m。漁獲物鑒定到種，并記錄每一個(gè)種類的尾數(shù)和重量。數(shù)據(jù)分析前，將原始數(shù)據(jù)標(biāo)準(zhǔn)化為1 h的漁獲尾數(shù)和漁獲重量。

1.2 物種多樣性指數(shù)

采用物種種類數(shù)(S)、Shannon-Weaver 多樣性指數(shù)(H′)、Pielou 均勻度指數(shù)(J)和Margalef 豐富度指數(shù)(d)對(duì)游泳動(dòng)物多樣性進(jìn)行分析，公式如下：

圖1 黃海和東海北部調(diào)查站位圖 Fig.1 Survey stations in the Yellow Sea and north East China Sea

式中，S 為每個(gè)站位游泳動(dòng)物種類數(shù)，Pi為第i種游泳動(dòng)物重量占總重量的比例，N 為所有游泳動(dòng)物的總重。

1.3 空間插值法

1.3.1 反距離加權(quán)插值 IDW 基于相近相似原理(Chaplot et al, 2006)，假設(shè)每一個(gè)樣本點(diǎn)的值都會(huì)對(duì)鄰近區(qū)域產(chǎn)生影響，距離未知點(diǎn)越近的樣本點(diǎn)對(duì)其影響程度越大。公式如下(Shepard, 1968)：

以統(tǒng)計(jì)學(xué)上的無(wú)偏性和最優(yōu)性為條件，結(jié)合拉格朗日原理，可推導(dǎo)出權(quán)重值和半方差函數(shù)之間關(guān)系的公式(王政權(quán), 1999)。

1.3.5 協(xié)同克里格插值 COK 是利用主變量之間的相關(guān)性及其與輔助變量的交互相關(guān)性，優(yōu)化插值的精度(姜勇等, 2006)，常用于主變量數(shù)據(jù)難以獲取或直接用主變量插值的結(jié)果不準(zhǔn)確的情況(李潤(rùn)林等, 2013)。公式如下(Myers, 1982)：

1.4 數(shù)據(jù)的正態(tài)分布檢驗(yàn)和轉(zhuǎn)換

進(jìn)行OK 和COK 之前，數(shù)據(jù)必須符合正態(tài)分布，因此，對(duì)上文計(jì)算的4 種物種多樣性指數(shù)數(shù)據(jù)進(jìn)行Shapiro-Wilk 檢驗(yàn)(Chen et al, 2016)。若不滿足正態(tài)分布，則進(jìn)行數(shù)據(jù)轉(zhuǎn)換。因?yàn)閿?shù)據(jù)結(jié)構(gòu)存在差異，本研究采用BOX-COX 和平方2 種轉(zhuǎn)換方式(于洋等, 2015; 章清等, 2016)。

1.5 效果評(píng)價(jià)

采用jackknife 檢驗(yàn)比較不同空間插值法的預(yù)測(cè)效果。該檢驗(yàn)的主要步驟為：首先，對(duì)于已有的n 個(gè)樣本，依次去除其中1 個(gè)樣本點(diǎn)，利用剩余的n?1 個(gè)樣本作為訓(xùn)練集對(duì)去除點(diǎn)進(jìn)行估算，重復(fù)n 次，得到所有樣本點(diǎn)的預(yù)測(cè)值；其次，計(jì)算預(yù)測(cè)值與樣本實(shí)際觀測(cè)值的決定系數(shù)(Coefficient of determination, R2)、平均絕對(duì)誤差(Mean absolute error, MAE)和均方根誤差(Root mean square error, RMSE)。R2反映樣本值和預(yù)測(cè)值之間的相關(guān)程度，R2越接近1，表示該空間插值法的預(yù)測(cè)值與樣本值的變化趨勢(shì)越相似，反之則偏差越大；MAE 可以反映預(yù)測(cè)值可能的誤差范圍；MAE和RMSE 都可以反映插值結(jié)果的精度，且都是越小越好(林忠輝等, 2002; 董志南等, 2015; 于洋等, 2015)。公式如下：

物種多樣性指數(shù)計(jì)算及其描述性統(tǒng)計(jì)、MAE及R2的計(jì)算、數(shù)據(jù)的正態(tài)分布檢驗(yàn)和轉(zhuǎn)換采用R3.5.2完成；半方差函數(shù)的擬合采用GS+7.0軟件完成，因GS+7.0軟件無(wú)法計(jì)算具有多個(gè)輔助變量條件下的半方差函數(shù)，因此，COK的半變異函數(shù)中的C0、C和C0/C 3個(gè)參數(shù)是未知的，其變程的計(jì)算是由ArcGIS10.2完成；4種多樣性指數(shù)之間的相關(guān)性判斷采用SPSS17.0完成；坐標(biāo)轉(zhuǎn)換、空間插值和交叉驗(yàn)證采用ArcGIS10.2完成。

2 結(jié)果

2.1 數(shù)據(jù)的統(tǒng)計(jì)描述

4 種多樣性指數(shù)的描述性統(tǒng)計(jì)見(jiàn)表1，Shpiro-Wilk正態(tài)檢驗(yàn)發(fā)現(xiàn)，S、J 和d 屬于偏態(tài)分布(P<0.05)。根據(jù)數(shù)據(jù)特點(diǎn)，對(duì)S 和d 進(jìn)行BOX-COX 轉(zhuǎn)換，對(duì) J 進(jìn)行平方轉(zhuǎn)換，轉(zhuǎn)換前對(duì)數(shù)據(jù)進(jìn)行加1 處理。轉(zhuǎn) 換后，數(shù)據(jù)均通過(guò)Shpiro-Wilk 正態(tài)檢驗(yàn)(P>0.1)，符合正態(tài)分布。此外，H′符合正態(tài)分布(P>0.05)，不做轉(zhuǎn)換。

表1 物種多樣性指數(shù)描述性統(tǒng)計(jì)和Shapiro-Wilk 檢驗(yàn)的P 值 Tab.1 Descriptive statistics of five species diversity indices and P-values of the Shapiro-Wilk test

2.2 協(xié)同克里格法輔助變量的選擇

表2 列舉了物種多樣性指數(shù)之間的相關(guān)性，其中，S、H′和d 之間具有極顯著的正相關(guān)關(guān)系。J 與H′呈極顯著正相關(guān)關(guān)系，與d 呈顯著正相關(guān)，而與S 相關(guān)性不顯著。將具有極顯著相關(guān)的物種多樣性指數(shù)作為輔助變量，當(dāng)存在多個(gè)具有顯著相關(guān)的輔助變量

時(shí)，根據(jù)ArcGIS10.2 交叉驗(yàn)證的結(jié)果初步進(jìn)行篩選。S 和J 最適合的輔助變量均為H′和d；H′最適合的輔助變量為S 和J；d 最適合的輔助變量為S、H′和J。

2.3 克里格插值半方差函數(shù)分析

根據(jù)GS+7.0 軟件的擬合結(jié)果(表3)，不同物種多樣性指數(shù)的半方差函數(shù)適合不同的模型。S 的空間分布為中等空間相關(guān)性(C0/C=0.43)，其他3 個(gè)指數(shù)均為強(qiáng)空間相關(guān)性(C0/C<0.25)，符合空間插值的前提。不同物種多樣性指數(shù)對(duì)應(yīng)的變程差別較大：S 和d 的變程較大，說(shuō)明其在大范圍內(nèi)的空間自相關(guān)性較強(qiáng)，H′和J 的變程都較小，說(shuō)明在小尺度范圍內(nèi)二者的空間自相關(guān)性較強(qiáng)。

2.4 不同空間插值法的比較

對(duì)于S，5 種空間插值法按預(yù)測(cè)值與樣本值的R2大小排序?yàn)長(zhǎng)PI

表3 物種多樣性指數(shù)的半方差函數(shù)及相關(guān)參數(shù) Tab.3 Semivariance models and the corresponding parameters of species diversity indices

4 種多樣性指數(shù)在樣本點(diǎn)的大小分布與預(yù)測(cè)的空間分布基本一致，在非采樣區(qū)域的曲線平滑度都低于采樣區(qū)域(圖2)。S 和d 整體的空間分布格局相似，都呈現(xiàn)出明顯的東南部高，北部次之，中部最低的格局。H′和J 也是中部最低，但是在南部和北部海域差異不明顯。4 種多樣性指數(shù)在32°N 以南海域高值區(qū)的范圍大致相同。

3 討論與分析

3.1 黃海和東海北部海域游泳動(dòng)物多樣性空間分布

分析COK的預(yù)測(cè)圖，在大尺度上，2017年春季黃海和東海北部游泳動(dòng)物的S和d大致呈現(xiàn)明顯的南高北低的格局，而H′和J的南北差異不大。與臺(tái)灣海峽游泳動(dòng)物多樣性(宋普慶等, 2012)以及相鄰海域魚(yú)類群落多樣性(程濟(jì)生, 2000)相比，除了d具有相同的空間格局，后2種多樣性指數(shù)空間分布格局均有較大差別。在小尺度上，34°～35.5°N，121°～124°E附近海域的S和d明顯比周圍海域低，而H′和J在此范圍內(nèi)沒(méi)有明顯的低值區(qū)。在30°～32°N范圍內(nèi)的海域中，4種多樣性指數(shù)均有高值區(qū)，而且西部海域均低于東部海域。在劉勇等(2017)的研究中，30°～32°N海域底層漁業(yè)生物的S和d也存在高值區(qū)，而且S、H′和d三者也呈現(xiàn)東部海域高于西部海域的格局，但是J在124°～ 126°E之間海域存在大范圍的低值區(qū)，與本研究結(jié)果有較大差別。

游泳動(dòng)物多樣性指數(shù)空間分布變化特征可能與環(huán)境、優(yōu)勢(shì)種、餌料豐度及分布、棲息地狀況、人為活動(dòng)和魚(yú)類洄游行為等多種因素有關(guān)(Jin et al, 1996; 戴芳群等, 2020; 金顯仕等, 2009; 程濟(jì)生等, 2004; 朱鑫華等, 1994)。從水文環(huán)境方面看，因?yàn)槎鄻有灾笖?shù)與環(huán)境因子的相關(guān)關(guān)系可以量化，而且數(shù)據(jù)容易獲取，常受到人們關(guān)注。已有研究表明，春季黃海中南部的S和d與海底水溫顯著正相關(guān)(P<0.05)，與深度、海表溫度和鹽度關(guān)系不顯著(P>0.05)(單秀娟等, 2014)，說(shuō)明物種多樣性的空間分布可能受海底水溫影響程度較大；而H’和J與深度、海底和海表的水溫和鹽度的關(guān)系不密切(單秀娟等, 2014)。受海底地形、海流等因素影響，在春、夏和秋季，黃海中北部海域存在一個(gè)穩(wěn)定且大范圍的冷水團(tuán)(李洪波等, 2006; 于非等, 2006; 王家棟等, 2011)。而在2017年春季，35°N，122°～124°E附近海域的海底溫度最高為10℃，低于附近海域的溫度(孔凡洲等, 2018)。這可能是導(dǎo)致4種多樣性指數(shù)在大尺度上以及34°～35.5°N，122°～124°E附近海域呈現(xiàn)獨(dú)特的分布特點(diǎn)和差異的原因。對(duì)于30°～32°N范圍內(nèi)的東海北部海域，多樣性指數(shù)較低的原因可能是靠近長(zhǎng)江口受淡水沖擊導(dǎo)致附近海域環(huán)境條件不穩(wěn)定，海水濁度較高，不利于生物生活(Whitfield, 1994; 劉淑德等, 2009, 徐勇等, 2014)。從繁殖特性方面看，春季是很多游泳動(dòng)物繁殖的季節(jié)，繁殖群體會(huì)洄游到環(huán)境適宜的海域產(chǎn)卵。如位于32°～34°N、122.5°E以西海域的呂四漁場(chǎng)是銀鯧(Pampus argenteus)、小黃魚(yú)(Larimichthys ployactis)和黃鯽(Setipinna tenuifilis)等種類主要產(chǎn)卵場(chǎng)(陳新軍, 2004; 劉勇等, 2017)。在30°～32.5°N，123°～126°E之間的海域存在多個(gè)重要漁場(chǎng)，例如，舟山漁場(chǎng)、江外漁場(chǎng)和舟外漁場(chǎng)等，是多種重要漁業(yè)種類的產(chǎn)卵場(chǎng)(劉勇等, 2007、2014; 林龍山等, 2008)。影響游泳動(dòng)物空間分布的原因復(fù)雜多變，對(duì)于其他小尺度上的多樣性指數(shù)空間特征形成原因，還需要進(jìn)一步研究和分析。

表4 不同插值方法的交叉驗(yàn)證結(jié)果 Tab.4 Cross-validation diagnostic of different interpolation methods

3.2 空間插值法分析

插值結(jié)果的精度取決于空間插值模型對(duì)變量空間變異性和相關(guān)性的解釋度(朱會(huì)義等, 2004)。本研究中，4 種多樣性指數(shù)的COK jackknife 檢驗(yàn)結(jié)果都比較好，解釋度高于其他方法，這與很多研究結(jié)果一致(盧月明等, 2017; 李潤(rùn)林等, 2013; 石朋等, 2005; 岳文澤等, 2005)。但是，在日降水量(宋麗瓊等, 2008)和熱帶雨林植物多樣性(Hernandez-Stefanoni et al, 2006)的空間插值法比較研究中，COK 的精度要低于OK。根據(jù)王平等(2014)對(duì)不同輔助變量的COK 插值結(jié)果比較，發(fā)現(xiàn)當(dāng)主變量和輔助變量之間的相關(guān)性大于0.45 時(shí)，COK 精度才會(huì)明顯高于OK?？梢?jiàn)，輔助變量的選擇是COK 插值的關(guān)鍵步驟，合適的輔助變量可以優(yōu)化COK。目前，最簡(jiǎn)單的方法是直接將與主變量相關(guān)性較高的因子作為輔助變量。如果主變量有多個(gè)顯著相關(guān)的因子，可以根據(jù)ArcGIS 軟件交叉驗(yàn)證的結(jié)果進(jìn)行初步篩選。但是，當(dāng)與主變量顯著相關(guān)的因子數(shù)量過(guò)多時(shí)，這種篩選方式會(huì)導(dǎo)致工作量增加，而且單一的輔助變量可能對(duì)模型優(yōu)化程度有限(章清等, 2016)。因此，基于主成分分析的COK 受到廣泛關(guān)注。這種方法是利用主成分分析，將多個(gè)與主變量顯著相關(guān)的因子賦予不同的權(quán)重，再合并成綜合指標(biāo)作為輔助變量(盧月明等, 2017)，在魚(yú)類資源豐度和土壤要素的空間分布研究中都取得了較好的結(jié)果(方學(xué)燕等, 2017; 章清等, 2016)。

圖2 協(xié)同克里格插值(COK)預(yù)測(cè)出的物種種類數(shù)(S)、Shannon-Weaver 多樣性指數(shù)(H′)、 Pielou 均勻度指數(shù)(J)和Margalef 豐富度指數(shù)(d)的空間分布圖 Fig.2 Spatial distribution of number of species (S)、Shannon-Weaver diversity index (H′)、 Pielou′s evenness index (J) and Margalef′s richness index (d) by ordinary cokriging (COK)

3.3 展望

本研究采用了多個(gè)輔助變量，主變量與輔助變量之間的相關(guān)性很高，因此，相對(duì)于OK，COK 的精度有了極大提高。研究表明，S、H′、J 和d 之間存在較高的相關(guān)性，這種現(xiàn)象在不同季節(jié)和年份同樣存在(單秀娟等, 2014)，因此，COK 不僅可以用于分析和預(yù)測(cè)春季黃海和東海北部游泳動(dòng)物多樣性空間分布，還可能適用于不同時(shí)間維度。在今后的研究中，可以結(jié)合與游泳動(dòng)物多樣性有顯著關(guān)系的環(huán)境因子，并基于主成分分析法減少輔助變量的數(shù)量，嘗試進(jìn)一步提高COK 的精度。