杜嘉星,孫 義,向 波,陳建軍,秦 彧,侯秀敏,于紅妍,宜樹(shù)華
(1.冰凍圈科學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室 / 中國(guó)科學(xué)院西北生態(tài)環(huán)境資源研究院,甘肅 蘭州 730000;2.中國(guó)科學(xué)院大學(xué),北京 100049;3.南通大學(xué)地理科學(xué)學(xué)院,江蘇 南通 226007;4.南通大學(xué)脆弱生態(tài)環(huán)境研究所,江蘇 南通 226007;5.重慶市氣候中心,重慶 401147;6.桂林理工大學(xué)測(cè)繪地理信息學(xué)院,廣西 桂林 541004;7.青海省草原總站,青海 西寧 810008)
高原鼠兔(Ochotona curzoniae)是分布在青藏高原及其周邊地區(qū)特有的小型草食哺乳動(dòng)物,是高寒草甸生態(tài)系統(tǒng)的關(guān)鍵種[1]。棲息地海拔為3 000-5 100 m,喜歡居住在坡位較低、土壤疏松干燥、地形開(kāi)闊、距離水源地較近、植被低矮的高山草甸和亞高山草甸[2-4]。高原鼠兔對(duì)其所處生態(tài)系統(tǒng)的影響主要取決于其種群密度。當(dāng)種群密度過(guò)高時(shí),其啃食牧草、挖掘鼠洞,常形成生物災(zāi)害,能直接導(dǎo)致草地退化;而當(dāng)種群密度在合理范圍內(nèi),能作為青藏高原中肉食性動(dòng)物的重要食物來(lái)源,其所挖掘的洞道為鳥(niǎo)類(lèi)、兩棲類(lèi)動(dòng)物提供了棲息場(chǎng)所,同時(shí)洞道的挖掘行為能增加入滲、減少?gòu)搅鳎黾油寥浪趾吞己?,?duì)生態(tài)系統(tǒng)食物網(wǎng)結(jié)構(gòu)及其能量流通和物質(zhì)循環(huán)具有重要意義[5-11]。
認(rèn)識(shí)高原鼠兔的空間分布及其影響因子對(duì)了解草地退化的原因、“黑土灘”的形成和高原鼠兔對(duì)棲息地的選擇以及在生態(tài)系統(tǒng)中的作用有重要意義。傳統(tǒng)的高原鼠兔調(diào)查方法有直接計(jì)數(shù)法、標(biāo)志重捕法、指數(shù)估計(jì)法[12]。這些方法雖然準(zhǔn)確度較高,但費(fèi)時(shí)費(fèi)力,難以展開(kāi)大范圍調(diào)查。無(wú)人機(jī)是利用無(wú)線電遙控設(shè)備和自備的程序控制裝置操縱的,或者由車(chē)載計(jì)算機(jī)完全或間歇地自主操作的不載人飛機(jī)[13]。無(wú)人機(jī)應(yīng)用前景廣闊,如無(wú)人偵察機(jī)和靶機(jī),以及航拍、農(nóng)業(yè)植保、測(cè)繪等無(wú)人機(jī)[14-15]。在生物群落調(diào)查中,無(wú)人機(jī)具有成本低、效率高等特點(diǎn),是進(jìn)行生物群落快速調(diào)查的優(yōu)良工具。例如,Watts等[16]用無(wú)人機(jī)影像來(lái)估計(jì)美國(guó)短吻鱷(Alligator mississippiensis)的種群數(shù)量,Israel[17]用無(wú)人機(jī)裝載熱紅外相機(jī)來(lái)監(jiān)測(cè)西方狍(Capreolus capreolus)的活動(dòng),郭新磊等[18]在青藏高原地區(qū)用無(wú)人機(jī)獲取高原鼠兔的洞口數(shù)量。
物種分布模型 (species distribution models, SDMs)是利用物種的觀測(cè)數(shù)據(jù)與環(huán)境數(shù)據(jù)基于統(tǒng)計(jì)或理論推導(dǎo)得出的經(jīng)驗(yàn)?zāi)P蚚19]。物種數(shù)據(jù)可以是基于隨機(jī)或分層野外采樣的簡(jiǎn)單存在、存在/缺失或物種豐富度觀測(cè)數(shù)據(jù)[20]。環(huán)境預(yù)測(cè)因子可以直接或間接影響物種,按照影響程度的大小排序,可以分為限制因子、擾動(dòng)和資源[21-22]。目前,物種分布模型已經(jīng)成為基礎(chǔ)生態(tài)學(xué)和生物地理學(xué)研究的重要工具,全球變化背景下物種的分布和氣候之間的關(guān)系、區(qū)域氣候變化對(duì)植物群落和功能的影響、生態(tài)系統(tǒng)功能群和關(guān)鍵種的監(jiān)測(cè)和預(yù)測(cè)、生態(tài)系統(tǒng)不同尺度多樣性的管理和保護(hù)、外來(lái)物種入侵區(qū)域的預(yù)測(cè)、生態(tài)系統(tǒng)的關(guān)鍵物種的潛在分布預(yù)測(cè)和保護(hù)區(qū)規(guī)劃等方面的研究上均有廣泛應(yīng)用[23]。BIOMOD(BIOdiversity MODelling)是一個(gè)基于R語(yǔ)言的用于物種分布集成預(yù)測(cè)的計(jì)算機(jī)平臺(tái),能夠處理一系列的模型中的方法不確定性和物種-環(huán)境關(guān)系的檢驗(yàn),由Wilfried Thuiller博士于2003年提出[24-25]。它可以使用多種不同方式建模,并對(duì)模型進(jìn)行評(píng)估,從而使模型精度最大化并預(yù)測(cè)物種分布。目前該模型的最新版本為2016年發(fā)布的BIOMOD2(Version:3.3-7)。本研究應(yīng)用BIOMOD物種分布集成模型集成平臺(tái),結(jié)合無(wú)人機(jī)遙感技術(shù)(unmanned aerial vehicle remote sensing, UAVRS),對(duì)高原鼠兔在黃河源區(qū)的潛在分布進(jìn)行研究,揭示其空間分布的規(guī)律并探究影響其空間分布的限制因子。
黃河源區(qū)位于青藏高原東北部邊緣(95°50′-103°30′ E, 32°30′-36°10′ N),橫跨青海、甘肅、四川三省,平均海拔 4 100 m,年平均氣溫-2.7 ℃,年平均降水517 mm。流域范圍達(dá)12.2萬(wàn)km2,是黃河流域的重要產(chǎn)流區(qū)。研究區(qū)植被類(lèi)型以高寒草甸和高山草原化草甸為主,優(yōu)勢(shì)種有矮生嵩草(Kobresia humilis)、垂穗披堿草(Elymus nutans)、青藏苔草(Carex moorcroftii)、異針茅(Stipa aliena)、高原早熟禾(Poa alpigena)等[26-27]。本研究區(qū)界定在河源至唐乃亥水文站的黃河流域范圍。
1.2.1 高原鼠兔分布數(shù)據(jù)
本研究所采用的高原鼠兔分布數(shù)據(jù)來(lái)自于2018年7-8月進(jìn)行的野外調(diào)查,在研究區(qū)范圍內(nèi)共布設(shè)了208個(gè)工作點(diǎn),使用大疆創(chuàng)新科技有限公司 (DJ-Innovations, DJI)生產(chǎn)的 Phantom 3 Professional和Mavic Pro無(wú)人機(jī)進(jìn)行航拍。由無(wú)人機(jī)生態(tài)環(huán)境研究網(wǎng)絡(luò)團(tuán)隊(duì)自主研發(fā)的FragMAP軟件設(shè)置工作點(diǎn)及航線,控制無(wú)人機(jī)按照設(shè)定航線飛行并定點(diǎn)拍照[28],航線的中心點(diǎn)設(shè)置與MODIS像元中心點(diǎn)相同,飛行高度為20 m,覆蓋地面范圍約為910 m2(35 m × 26 m),單張照片分辨率為 4 000 像素 × 3 000像素。后期使用自主研發(fā)的Proposal Classifier軟件對(duì)照片進(jìn)行處理,選取合適的閾值,通過(guò)自動(dòng)識(shí)別和人工手動(dòng)矯正相結(jié)合的方法對(duì)高原鼠兔洞口進(jìn)行提取,并將處于同一航線上的洞口數(shù)量進(jìn)行平均,結(jié)合當(dāng)?shù)氐挠行Ф纯谙禂?shù),獲取高原鼠兔的存在/不存在點(diǎn)。
1.2.2 環(huán)境數(shù)據(jù)
本研究共獲取了48種環(huán)境數(shù)據(jù),包括以下幾種類(lèi)型。
1)氣候數(shù)據(jù):氣候數(shù)據(jù)來(lái)自中國(guó)科學(xué)院青藏高原研究所陽(yáng)坤團(tuán)隊(duì)[29-30],包括年平均氣溫、氣溫周較差、等溫性、溫度季節(jié)性、最高氣溫、最低氣溫、氣溫年較差、最濕季平均氣溫、最干季平均氣溫、最暖季平均氣溫、最冷季平均氣溫、年降水量、最濕期降水量、最干期降水量、降水季節(jié)性、最濕季降水量、最干季降水量、最暖季降水量、最冷季降水量、年平均輻射、最大輻射、最小輻射、輻射季節(jié)性、最濕季輻射、最干季輻射、最暖季輻射、最冷季輻射、年平均濕度、最大濕度、濕度季節(jié)性、最濕季濕度、最干季濕度和最暖季濕度的0.1°格網(wǎng)數(shù)據(jù),采用克里金插值的方法擴(kuò)展到整個(gè)研究區(qū),空間分辨率1 km,共33個(gè)柵格圖層。
2)植被數(shù)據(jù):植被數(shù)據(jù)來(lái)自LAADS DAAC(ladsweb.modaps.eosdis.nasa.gov)的 MODIS/Terra MOD1-3A3植被指數(shù)產(chǎn)品,空間分辨率1 km,時(shí)間分辨率30 d。時(shí)間范圍為2008年11月-2018年10月(共11 年),通過(guò)使用 MRT(Modis projection tool)軟件進(jìn)行讀取并進(jìn)行投影,在QGIS Desktop中進(jìn)行鑲嵌和裁剪,最后通過(guò)像元統(tǒng)計(jì)得到NDVI年平均值、NDVI最大值、NDVI最小值、NDVI年較差共4個(gè)柵格圖層。
3)地形數(shù)據(jù):數(shù)字高程數(shù)據(jù)(digital elevation model, DEM)為 SRTM 90 m DEM 數(shù)據(jù),來(lái)自美國(guó)地質(zhì)勘探局(www.usgs.gov)。使用QGIS Desktop對(duì)DEM數(shù)據(jù)提取坡度和坡向,得到海拔、坡度、坡向共3個(gè)柵格圖層。
4)土壤數(shù)據(jù):土壤數(shù)據(jù)來(lái)自SoilGrids(www.soilgrids.org)[31]1 km分辨率的土層厚度、0.3-0.6 m土層有機(jī)碳儲(chǔ)量、0.3 m處土壤容重、0.3 m處黏土含量、0.3 m處粗屑體容積、0.3 m處淤泥含量、0.3 m處含沙量、0.3 m處土壤pH數(shù)據(jù),共8個(gè)柵格圖層。
在使用上述環(huán)境數(shù)據(jù)構(gòu)建模型過(guò)程中,為了減少因子自相關(guān)對(duì)模型精度的影響,對(duì)其進(jìn)行Pearson 相關(guān)性分析,保留|r| > 0.7 的環(huán)境變量對(duì)中的其中一個(gè)。經(jīng)過(guò)篩選后剩下年平均氣溫、氣溫周較差、等溫性、溫度季節(jié)性、年降水量、最濕期降水量、最干期降水量、年平均輻射、輻射季節(jié)性、海拔、坡度、坡向、NDVI年平均、NDVI年較差、有機(jī)碳儲(chǔ)量、土壤容重、土壤pH共17種環(huán)境數(shù)據(jù)用于建模。
在BIOMOD中共提供了10種物種分布模型可供選擇,并且可以通過(guò)不同指標(biāo)計(jì)算模型精度來(lái)評(píng)估不同模型在研究區(qū)的適用性,從而篩選出最優(yōu)的模型,包括廣義線性模型(generalized linear model,GLM)[32]、廣義增強(qiáng)回歸模型 (generalized boosted regression models, GBM)[33-34]、廣義加性模型 (generalized additive model, GAM)[35]、 分 類(lèi) 樹(shù) 分 析 (classification tree analysis, CTA)[36]、人工神經(jīng)網(wǎng)絡(luò) (artificial neural network, ANN)[37]、表面分布區(qū)分室模型 (surface range envelope, SRE)[38]、彈性判別分析 (flexible discriminant analysis, FDA)[39]、多元自適應(yīng)回歸樣條 (multiple adaptive regression splines, MARS)[40-41]、隨機(jī)森林 (random forest, RF)[42]和最大熵模型 (maximum entropy model,MaxEnt)[43]。
為了評(píng)價(jià)模型的預(yù)測(cè)精度,將樣本劃分為兩個(gè)子集,一個(gè)子集包含全部樣本的70%,另一個(gè)包含30%。樣本的70%用于建模,30%用于檢驗(yàn),并引入真實(shí)技巧統(tǒng)計(jì)值 (true skill statistics, TSS)和 AUC(area under curve)兩種指標(biāo)來(lái)對(duì)模型進(jìn)行評(píng)估。
1.4.1 TSS
TSS = Sensitivity + Specificity - 1 = TPR - FPR。其中,Sensitivity為敏感性,Specificity為特異性,TPR(true positive rate)為真陽(yáng)性率,F(xiàn)PR(false positive rate)為假陽(yáng)性率[44]。TSS 的取值范圍在 (0, 1)區(qū)間,TSS的值越接近于1,表示真陽(yáng)性率與假陽(yáng)性率的差值越大,模型的效果就越好。
1.4.2 AUC
AUC來(lái)源于受試者工作特征曲線(receiver operating characteristic curve, ROC),被定義為 ROC 曲線下與坐標(biāo)軸圍成的面積[45]。ROC曲線是反映敏感性和特異性連續(xù)變量的綜合指標(biāo),ROC曲線上每個(gè)點(diǎn)反映著對(duì)同一信號(hào)刺激的感受性。AUC的取值范圍一般在(0.5, 1)區(qū)間內(nèi)。越接近于1的模型表示預(yù)測(cè)效果越好,而越接近于0.5的模型表示越接近于隨機(jī)猜測(cè),模型沒(méi)有預(yù)測(cè)價(jià)值。因此,AUC值越大的模型效果越好。
在實(shí)地調(diào)查中發(fā)現(xiàn),如選用是否存在高原鼠兔洞口作為判別高原鼠兔是否存在的依據(jù),則完全沒(méi)有洞口的樣本量非常少,故在本研究中選用是否存在1個(gè)或更多有效洞口作為判別依據(jù)。根據(jù)Yi等[46]在青藏高原地區(qū)不同草地類(lèi)型的調(diào)查,該地區(qū)的有效洞口系數(shù)在0.22~0.42。根據(jù)泊松分布,10個(gè)洞口中存在1個(gè)或更多有效洞口的概率P> 0.96,可以認(rèn)為10個(gè)洞口中至少存在1個(gè)活動(dòng)的洞口。故對(duì)于洞口個(gè)數(shù)大于10的樣地認(rèn)為有高原鼠兔存在,小于10則認(rèn)為不存在。通過(guò)對(duì)黃河源區(qū)野外獲取的航拍照片進(jìn)行洞口提取,共獲得了82個(gè)存在點(diǎn)數(shù)據(jù)和126個(gè)不存在點(diǎn)數(shù)據(jù),并根據(jù)航拍工作點(diǎn)的定位進(jìn)行繪圖(圖1)。
圖1 黃河源區(qū)采樣點(diǎn)Figure 1 Sampling points in the source region of the Yellow River Basin
將10種模型每個(gè)重復(fù)運(yùn)行多次,每次隨機(jī)選取70%的樣本進(jìn)行建模,30%的樣本對(duì)建模的結(jié)果進(jìn)行評(píng)估。10種模型中除GAM由于獲取估計(jì)參數(shù)困難運(yùn)行失敗外,其余9種模型均成功運(yùn)行,將9種模型運(yùn)行多次得到的TSS與AUC兩個(gè)參數(shù)進(jìn)行描述性統(tǒng)計(jì)(表1),并繪制箱形圖(圖2)。在9種模型中,RF的整體表現(xiàn)最佳,TSS平均值為0.57,AUC平均值為0.83,變異系數(shù)僅為11%和5%,說(shuō)明在多次運(yùn)行的過(guò)程中模型能較準(zhǔn)確地預(yù)測(cè)并能保持穩(wěn)定性。次優(yōu)的模型為GBM、MARS和FDA,這3種模型均能一定程度上預(yù)測(cè)高原鼠兔的分布。其余模型均表現(xiàn)一般,其中SRE最差,TSS為0.09,AUC 為 0.54,接近于隨機(jī)分類(lèi) (AUC = 0.5),可以認(rèn)為該模型基本無(wú)法在該地區(qū)準(zhǔn)確預(yù)測(cè)。
使用模型內(nèi)置的因子重要性計(jì)算函數(shù),將包含全部環(huán)境因子的集合定義為ref數(shù)據(jù)集,而隨機(jī)剔除單個(gè)環(huán)境因子后的集合定義為shuffled數(shù)據(jù)集,使用兩個(gè)集合進(jìn)行預(yù)測(cè),并計(jì)算預(yù)測(cè)結(jié)果的簡(jiǎn)單相關(guān)性(默認(rèn)為Pearson相關(guān)性)。計(jì)算公式:Importance = 1 - cor(pred_ref, pred_shuffled)。剔除掉表現(xiàn)較差的SRE以及運(yùn)行失敗的GAM,統(tǒng)計(jì)其余8種模型的各個(gè)環(huán)境因子的重要性求平均并按大小排序(表2)。
表1 9種模型TSS與AUC描述性統(tǒng)計(jì)Table 1 Descriptive statistics of the TSS and AUC of nine models
圖2 9種模型TSS與AUC比較Figure 2 Comparison of the TSS and AUC of nine models
表2 8種模型環(huán)境因子重要性Table 2 Importance of the environment factors in eight models
根據(jù)以上對(duì)模型精度的評(píng)估結(jié)果,采用準(zhǔn)確度和穩(wěn)定性均較好的RF。在多次重復(fù)運(yùn)行結(jié)果中選取預(yù)測(cè)效果最佳的一次結(jié)果 (TSS = 0.74, AUC = 0.92),對(duì)結(jié)果進(jìn)行繪圖(圖3)。
預(yù)測(cè)結(jié)果顯示,高原鼠兔分布概率較高的區(qū)域主要有瑪多縣扎陵湖鄉(xiāng)、瑪沁縣優(yōu)云鄉(xiāng)、同德縣尕巴松多鎮(zhèn)、河南縣優(yōu)干寧鎮(zhèn)以及達(dá)日縣與甘德縣交界處的上貢麻鄉(xiāng),而在瑪曲縣東南部、若爾蓋縣西部與紅原縣北部三縣之間的山區(qū)、瑪多縣瑪查里鎮(zhèn)東南部以及同德縣和澤庫(kù)縣北部則分布概率較小。
本研究使用無(wú)人機(jī)于2018年7-8月在黃河源區(qū)進(jìn)行野外調(diào)查,獲取了3 000余張航拍照片,從中提取了208個(gè)高原鼠兔的存在/不存在點(diǎn)數(shù)據(jù),從氣候、植被、地形和土壤4類(lèi)共48種環(huán)境因子篩選出17種作為解釋變量,使用BIOMOD集成模型平臺(tái)對(duì)該區(qū)域高原鼠兔的潛在分布進(jìn)行預(yù)測(cè)。本研究方法相比較傳統(tǒng)的地面調(diào)查效率更高且更適合大范圍調(diào)查,預(yù)測(cè)的結(jié)果能為鼠害防治工作提供科學(xué)依據(jù),同時(shí)能為今后的高原鼠兔研究提供新思路。
BIOMOD中集成的10種物種分布模型中表現(xiàn)最好的是RF,多次運(yùn)行的TSS平均值為0.57,AUC平均值為0.83,變異系數(shù)分別為11%和5%,說(shuō)明該模型在黃河源區(qū)能較為準(zhǔn)確地預(yù)測(cè)高原鼠兔的分布,并且在多次運(yùn)行中仍能保持較好的穩(wěn)定性。BIOMOD在此區(qū)域的應(yīng)用能降低預(yù)測(cè)的不確定性和誤差,提高預(yù)測(cè)的精度。
高原鼠兔的分布受多種環(huán)境因子影響,重要性較高的幾個(gè)因子有海拔、NDVI年平均、坡度、土壤pH、年降水量、坡向以及年平均輻射。高原鼠兔的分布與海拔高度相關(guān),其適合生長(zhǎng)在高海拔低氧地區(qū),王淯等[4]的研究表明,高原鼠兔多分布在海拔 3 000-5 100 m、地勢(shì)較為平坦且坡位較低的區(qū)域。高原鼠兔的分布同時(shí)受坡度和坡向的影響,衛(wèi)萬(wàn)榮[47]的研究表明,84.2%高原鼠兔喜好將棲息地選擇在坡度在0°~5°的區(qū)域,而在6°~10°和11°~15°的區(qū)域分布分別只占13.2%和2.6%;對(duì)坡向的選擇更傾向于南或偏南向;在山地地貌上占總數(shù)的100%,在所有地貌上占41.3%。這有利于增加冬季對(duì)太陽(yáng)輻射的獲取,提升洞內(nèi)溫度。NDVI能夠在一定程度上反映植被的蓋度和生物量,馬波等[48]的研究認(rèn)為高原鼠兔洞口數(shù)與NDVI存在顯著正相關(guān)的線性關(guān)系,高原鼠兔有選擇地利用NDVI較高、植被較好的生境。高原鼠兔的挖掘活動(dòng)能夠改變當(dāng)?shù)氐耐寥罓顩r,Guo等[49]的研究認(rèn)為高原鼠兔的洞穴對(duì)土壤pH、含水量、含沙量等土壤理化性質(zhì)均有影響,能增加土壤滲透率,并提高土壤養(yǎng)分。同時(shí)高原鼠兔的分布受降水量和輻射的影響,郭新磊[50]的研究表明,高原鼠兔洞口密度隨降水量和輻射的增高表現(xiàn)為先增大后減少,降水量在600~700 mm區(qū)間和輻射在240~250 J·m-2區(qū)間的地區(qū)高原鼠兔洞口密度達(dá)到最大。以上影響高原鼠兔空間分布的因子與前人的研究相一致,說(shuō)明BIOMOD不僅能很好地的模擬高原鼠兔的空間分布,在影響因子的確定上也有著良好的效果。
圖3 基于RF的黃河源區(qū)高原鼠兔分布預(yù)測(cè)Figure 3 Prediction of the distribution of plateau pika, based on RF in the source region of the Yellow River Basin
高原鼠兔分布除了與環(huán)境因子有關(guān)外,還與人類(lèi)活動(dòng)密切相關(guān),如放牧、害鼠防治、招鷹架的布置等。張興祿[51]在瑪曲高寒草甸植被區(qū)所做試驗(yàn)證明輪牧制度可有效抑制高原鼠兔數(shù)量增長(zhǎng),而在連續(xù)放牧制度下,低放牧率下高原鼠兔的數(shù)量則明顯上升,而高放牧率下由于牛羊踐踏、采食過(guò)重,不適合高原鼠兔生存,數(shù)量則有所下降。梁杰榮等[52]在海北高寒草甸生態(tài)系統(tǒng)定位研究站調(diào)查了高原鼠兔的種群恢復(fù)過(guò)程,發(fā)現(xiàn)害鼠防治后殘鼠數(shù)量按Logistics曲線增長(zhǎng);當(dāng)防治效果高達(dá)90%時(shí),其數(shù)量約在兩年內(nèi)也恢復(fù)至原水平。侯秀敏[53]在澤庫(kù)巴灘所做試驗(yàn)表明,招鷹架的布置能有效吸引高原鼠兔的天敵,降低高原鼠兔的密度。當(dāng)鷹架布置的間距為 500 m × 500 m 時(shí),控制面積為每25 hm2一根,鷹架的利用率達(dá)到了98%,樣地內(nèi)兩年間高原鼠兔的有效洞口數(shù)量降低了26.47%。
致謝:本研究使用的氣候數(shù)據(jù)由青藏高原多圈層模擬與數(shù)據(jù)同化中心 (Data Assimilation and Modeling Center for Tibetan Multi-spheres)開(kāi)發(fā),中國(guó)科學(xué)院青藏高原研究所陽(yáng)坤研究員提供,特在此表示誠(chéng)摯的謝意。