基于BIOMOD的黃河源區(qū)高原鼠兔潛在分布及其影響因子

杜嘉星，孫 義，向 波，陳建軍，秦 彧，侯秀敏，于紅妍，宜樹華

（1.冰凍圈科學國家重點實驗室 / 中國科學院西北生態(tài)環(huán)境資源研究院，甘肅 蘭州 730000；2.中國科學院大學，北京 100049；3.南通大學地理科學學院，江蘇 南通 226007；4.南通大學脆弱生態(tài)環(huán)境研究所，江蘇 南通 226007；5.重慶市氣候中心，重慶 401147；6.桂林理工大學測繪地理信息學院，廣西 桂林 541004；7.青海省草原總站，青海 西寧 810008）

高原鼠兔(Ochotona curzoniae)是分布在青藏高原及其周邊地區(qū)特有的小型草食哺乳動物，是高寒草甸生態(tài)系統(tǒng)的關鍵種[1]。棲息地海拔為3 000-5 100 m，喜歡居住在坡位較低、土壤疏松干燥、地形開闊、距離水源地較近、植被低矮的高山草甸和亞高山草甸[2-4]。高原鼠兔對其所處生態(tài)系統(tǒng)的影響主要取決于其種群密度。當種群密度過高時，其啃食牧草、挖掘鼠洞，常形成生物災害，能直接導致草地退化；而當種群密度在合理范圍內，能作為青藏高原中肉食性動物的重要食物來源，其所挖掘的洞道為鳥類、兩棲類動物提供了棲息場所，同時洞道的挖掘行為能增加入滲、減少徑流，增加土壤水分和碳含量，對生態(tài)系統(tǒng)食物網結構及其能量流通和物質循環(huán)具有重要意義[5-11]。

認識高原鼠兔的空間分布及其影響因子對了解草地退化的原因、“黑土灘”的形成和高原鼠兔對棲息地的選擇以及在生態(tài)系統(tǒng)中的作用有重要意義。傳統(tǒng)的高原鼠兔調查方法有直接計數(shù)法、標志重捕法、指數(shù)估計法[12]。這些方法雖然準確度較高，但費時費力，難以展開大范圍調查。無人機是利用無線電遙控設備和自備的程序控制裝置操縱的，或者由車載計算機完全或間歇地自主操作的不載人飛機[13]。無人機應用前景廣闊，如無人偵察機和靶機，以及航拍、農業(yè)植保、測繪等無人機[14-15]。在生物群落調查中，無人機具有成本低、效率高等特點，是進行生物群落快速調查的優(yōu)良工具。例如，Watts等[16]用無人機影像來估計美國短吻鱷(Alligator mississippiensis)的種群數(shù)量，Israel[17]用無人機裝載熱紅外相機來監(jiān)測西方狍(Capreolus capreolus)的活動，郭新磊等[18]在青藏高原地區(qū)用無人機獲取高原鼠兔的洞口數(shù)量。

物種分布模型 (species distribution models, SDMs)是利用物種的觀測數(shù)據(jù)與環(huán)境數(shù)據(jù)基于統(tǒng)計或理論推導得出的經驗模型[19]。物種數(shù)據(jù)可以是基于隨機或分層野外采樣的簡單存在、存在/缺失或物種豐富度觀測數(shù)據(jù)[20]。環(huán)境預測因子可以直接或間接影響物種，按照影響程度的大小排序，可以分為限制因子、擾動和資源[21-22]。目前，物種分布模型已經成為基礎生態(tài)學和生物地理學研究的重要工具，全球變化背景下物種的分布和氣候之間的關系、區(qū)域氣候變化對植物群落和功能的影響、生態(tài)系統(tǒng)功能群和關鍵種的監(jiān)測和預測、生態(tài)系統(tǒng)不同尺度多樣性的管理和保護、外來物種入侵區(qū)域的預測、生態(tài)系統(tǒng)的關鍵物種的潛在分布預測和保護區(qū)規(guī)劃等方面的研究上均有廣泛應用[23]。BIOMOD(BIOdiversity MODelling)是一個基于R語言的用于物種分布集成預測的計算機平臺，能夠處理一系列的模型中的方法不確定性和物種-環(huán)境關系的檢驗，由Wilfried Thuiller博士于2003年提出[24-25]。它可以使用多種不同方式建模，并對模型進行評估，從而使模型精度最大化并預測物種分布。目前該模型的最新版本為2016年發(fā)布的BIOMOD2(Version:3.3-7)。本研究應用BIOMOD物種分布集成模型集成平臺，結合無人機遙感技術(unmanned aerial vehicle remote sensing, UAVRS)，對高原鼠兔在黃河源區(qū)的潛在分布進行研究，揭示其空間分布的規(guī)律并探究影響其空間分布的限制因子。

1 材料與方法

1.1 研究區(qū)概況

黃河源區(qū)位于青藏高原東北部邊緣(95°50′-103°30′ E, 32°30′-36°10′ N)，橫跨青海、甘肅、四川三省，平均海拔 4 100 m，年平均氣溫-2.7 ℃，年平均降水517 mm。流域范圍達12.2萬km2，是黃河流域的重要產流區(qū)。研究區(qū)植被類型以高寒草甸和高山草原化草甸為主，優(yōu)勢種有矮生嵩草(Kobresia humilis)、垂穗披堿草(Elymus nutans)、青藏苔草(Carex moorcroftii)、異針茅(Stipa aliena)、高原早熟禾(Poa alpigena)等[26-27]。本研究區(qū)界定在河源至唐乃亥水文站的黃河流域范圍。

1.2 數(shù)據(jù)獲取

1.2.1 高原鼠兔分布數(shù)據(jù)

本研究所采用的高原鼠兔分布數(shù)據(jù)來自于2018年7-8月進行的野外調查，在研究區(qū)范圍內共布設了208個工作點，使用大疆創(chuàng)新科技有限公司 (DJ-Innovations, DJI)生產的 Phantom 3 Professional和Mavic Pro無人機進行航拍。由無人機生態(tài)環(huán)境研究網絡團隊自主研發(fā)的FragMAP軟件設置工作點及航線，控制無人機按照設定航線飛行并定點拍照[28]，航線的中心點設置與MODIS像元中心點相同，飛行高度為20 m，覆蓋地面范圍約為910 m2(35 m × 26 m)，單張照片分辨率為 4 000 像素 × 3 000像素。后期使用自主研發(fā)的Proposal Classifier軟件對照片進行處理，選取合適的閾值，通過自動識別和人工手動矯正相結合的方法對高原鼠兔洞口進行提取，并將處于同一航線上的洞口數(shù)量進行平均，結合當?shù)氐挠行Ф纯谙禂?shù)，獲取高原鼠兔的存在/不存在點。

1.2.2 環(huán)境數(shù)據(jù)

本研究共獲取了48種環(huán)境數(shù)據(jù)，包括以下幾種類型。

1)氣候數(shù)據(jù)：氣候數(shù)據(jù)來自中國科學院青藏高原研究所陽坤團隊[29-30]，包括年平均氣溫、氣溫周較差、等溫性、溫度季節(jié)性、最高氣溫、最低氣溫、氣溫年較差、最濕季平均氣溫、最干季平均氣溫、最暖季平均氣溫、最冷季平均氣溫、年降水量、最濕期降水量、最干期降水量、降水季節(jié)性、最濕季降水量、最干季降水量、最暖季降水量、最冷季降水量、年平均輻射、最大輻射、最小輻射、輻射季節(jié)性、最濕季輻射、最干季輻射、最暖季輻射、最冷季輻射、年平均濕度、最大濕度、濕度季節(jié)性、最濕季濕度、最干季濕度和最暖季濕度的0.1°格網數(shù)據(jù)，采用克里金插值的方法擴展到整個研究區(qū)，空間分辨率1 km，共33個柵格圖層。

2)植被數(shù)據(jù)：植被數(shù)據(jù)來自LAADS DAAC(ladsweb.modaps.eosdis.nasa.gov)的 MODIS/Terra MOD1-3A3植被指數(shù)產品，空間分辨率1 km，時間分辨率30 d。時間范圍為2008年11月-2018年10月(共11 年)，通過使用 MRT(Modis projection tool)軟件進行讀取并進行投影，在QGIS Desktop中進行鑲嵌和裁剪，最后通過像元統(tǒng)計得到NDVI年平均值、NDVI最大值、NDVI最小值、NDVI年較差共4個柵格圖層。

3)地形數(shù)據(jù)：數(shù)字高程數(shù)據(jù)(digital elevation model, DEM)為 SRTM 90 m DEM 數(shù)據(jù)，來自美國地質勘探局(www.usgs.gov)。使用QGIS Desktop對DEM數(shù)據(jù)提取坡度和坡向，得到海拔、坡度、坡向共3個柵格圖層。

4)土壤數(shù)據(jù)：土壤數(shù)據(jù)來自SoilGrids(www.soilgrids.org)[31]1 km分辨率的土層厚度、0.3-0.6 m土層有機碳儲量、0.3 m處土壤容重、0.3 m處黏土含量、0.3 m處粗屑體容積、0.3 m處淤泥含量、0.3 m處含沙量、0.3 m處土壤pH數(shù)據(jù)，共8個柵格圖層。

在使用上述環(huán)境數(shù)據(jù)構建模型過程中，為了減少因子自相關對模型精度的影響，對其進行Pearson 相關性分析，保留|r| ＞ 0.7 的環(huán)境變量對中的其中一個。經過篩選后剩下年平均氣溫、氣溫周較差、等溫性、溫度季節(jié)性、年降水量、最濕期降水量、最干期降水量、年平均輻射、輻射季節(jié)性、海拔、坡度、坡向、NDVI年平均、NDVI年較差、有機碳儲量、土壤容重、土壤pH共17種環(huán)境數(shù)據(jù)用于建模。

1.3 模型介紹

在BIOMOD中共提供了10種物種分布模型可供選擇，并且可以通過不同指標計算模型精度來評估不同模型在研究區(qū)的適用性，從而篩選出最優(yōu)的模型，包括廣義線性模型(generalized linear model,GLM)[32]、廣義增強回歸模型 (generalized boosted regression models, GBM)[33-34]、廣義加性模型 (generalized additive model, GAM)[35]、 分 類 樹 分 析 (classification tree analysis, CTA)[36]、人工神經網絡 (artificial neural network, ANN)[37]、表面分布區(qū)分室模型 (surface range envelope, SRE)[38]、彈性判別分析 (flexible discriminant analysis, FDA)[39]、多元自適應回歸樣條 (multiple adaptive regression splines, MARS)[40-41]、隨機森林 (random forest, RF)[42]和最大熵模型 (maximum entropy model,MaxEnt)[43]。

1.4 模型檢驗

為了評價模型的預測精度，將樣本劃分為兩個子集，一個子集包含全部樣本的70%，另一個包含30%。樣本的70%用于建模，30%用于檢驗，并引入真實技巧統(tǒng)計值 (true skill statistics, TSS)和 AUC(area under curve)兩種指標來對模型進行評估。

1.4.1 TSS

TSS = Sensitivity + Specificity - 1 = TPR - FPR。其中，Sensitivity為敏感性，Specificity為特異性，TPR(true positive rate)為真陽性率，F(xiàn)PR(false positive rate)為假陽性率[44]。TSS 的取值范圍在 (0, 1)區(qū)間，TSS的值越接近于1，表示真陽性率與假陽性率的差值越大，模型的效果就越好。

1.4.2 AUC

AUC來源于受試者工作特征曲線(receiver operating characteristic curve, ROC)，被定義為 ROC 曲線下與坐標軸圍成的面積[45]。ROC曲線是反映敏感性和特異性連續(xù)變量的綜合指標，ROC曲線上每個點反映著對同一信號刺激的感受性。AUC的取值范圍一般在(0.5, 1)區(qū)間內。越接近于1的模型表示預測效果越好，而越接近于0.5的模型表示越接近于隨機猜測，模型沒有預測價值。因此，AUC值越大的模型效果越好。

2 結果與分析

2.1 高原鼠兔分布

在實地調查中發(fā)現(xiàn)，如選用是否存在高原鼠兔洞口作為判別高原鼠兔是否存在的依據(jù)，則完全沒有洞口的樣本量非常少，故在本研究中選用是否存在1個或更多有效洞口作為判別依據(jù)。根據(jù)Yi等[46]在青藏高原地區(qū)不同草地類型的調查，該地區(qū)的有效洞口系數(shù)在0.22～0.42。根據(jù)泊松分布，10個洞口中存在1個或更多有效洞口的概率P＞ 0.96，可以認為10個洞口中至少存在1個活動的洞口。故對于洞口個數(shù)大于10的樣地認為有高原鼠兔存在，小于10則認為不存在。通過對黃河源區(qū)野外獲取的航拍照片進行洞口提取，共獲得了82個存在點數(shù)據(jù)和126個不存在點數(shù)據(jù)，并根據(jù)航拍工作點的定位進行繪圖(圖1)。

圖1 黃河源區(qū)采樣點Figure 1 Sampling points in the source region of the Yellow River Basin

2.2 模型精度評估

將10種模型每個重復運行多次，每次隨機選取70%的樣本進行建模，30%的樣本對建模的結果進行評估。10種模型中除GAM由于獲取估計參數(shù)困難運行失敗外，其余9種模型均成功運行，將9種模型運行多次得到的TSS與AUC兩個參數(shù)進行描述性統(tǒng)計(表1)，并繪制箱形圖(圖2)。在9種模型中，RF的整體表現(xiàn)最佳，TSS平均值為0.57，AUC平均值為0.83，變異系數(shù)僅為11%和5%，說明在多次運行的過程中模型能較準確地預測并能保持穩(wěn)定性。次優(yōu)的模型為GBM、MARS和FDA，這3種模型均能一定程度上預測高原鼠兔的分布。其余模型均表現(xiàn)一般，其中SRE最差，TSS為0.09，AUC 為 0.54，接近于隨機分類 (AUC = 0.5)，可以認為該模型基本無法在該地區(qū)準確預測。

2.3 環(huán)境因子重要性

使用模型內置的因子重要性計算函數(shù)，將包含全部環(huán)境因子的集合定義為ref數(shù)據(jù)集，而隨機剔除單個環(huán)境因子后的集合定義為shuffled數(shù)據(jù)集，使用兩個集合進行預測，并計算預測結果的簡單相關性(默認為Pearson相關性)。計算公式：Importance = 1 - cor(pred_ref, pred_shuffled)。剔除掉表現(xiàn)較差的SRE以及運行失敗的GAM，統(tǒng)計其余8種模型的各個環(huán)境因子的重要性求平均并按大小排序(表2)。

表1 9種模型TSS與AUC描述性統(tǒng)計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

2.4 預測結果

根據(jù)以上對模型精度的評估結果，采用準確度和穩(wěn)定性均較好的RF。在多次重復運行結果中選取預測效果最佳的一次結果 (TSS = 0.74, AUC = 0.92)，對結果進行繪圖(圖3)。

預測結果顯示，高原鼠兔分布概率較高的區(qū)域主要有瑪多縣扎陵湖鄉(xiāng)、瑪沁縣優(yōu)云鄉(xiāng)、同德縣尕巴松多鎮(zhèn)、河南縣優(yōu)干寧鎮(zhèn)以及達日縣與甘德縣交界處的上貢麻鄉(xiāng)，而在瑪曲縣東南部、若爾蓋縣西部與紅原縣北部三縣之間的山區(qū)、瑪多縣瑪查里鎮(zhèn)東南部以及同德縣和澤庫縣北部則分布概率較小。

3 討論與結論

本研究使用無人機于2018年7-8月在黃河源區(qū)進行野外調查，獲取了3 000余張航拍照片，從中提取了208個高原鼠兔的存在/不存在點數(shù)據(jù)，從氣候、植被、地形和土壤4類共48種環(huán)境因子篩選出17種作為解釋變量，使用BIOMOD集成模型平臺對該區(qū)域高原鼠兔的潛在分布進行預測。本研究方法相比較傳統(tǒng)的地面調查效率更高且更適合大范圍調查，預測的結果能為鼠害防治工作提供科學依據(jù)，同時能為今后的高原鼠兔研究提供新思路。

BIOMOD中集成的10種物種分布模型中表現(xiàn)最好的是RF，多次運行的TSS平均值為0.57，AUC平均值為0.83，變異系數(shù)分別為11%和5%，說明該模型在黃河源區(qū)能較為準確地預測高原鼠兔的分布，并且在多次運行中仍能保持較好的穩(wěn)定性。BIOMOD在此區(qū)域的應用能降低預測的不確定性和誤差，提高預測的精度。

高原鼠兔的分布受多種環(huán)境因子影響，重要性較高的幾個因子有海拔、NDVI年平均、坡度、土壤pH、年降水量、坡向以及年平均輻射。高原鼠兔的分布與海拔高度相關，其適合生長在高海拔低氧地區(qū)，王淯等[4]的研究表明，高原鼠兔多分布在海拔 3 000-5 100 m、地勢較為平坦且坡位較低的區(qū)域。高原鼠兔的分布同時受坡度和坡向的影響，衛(wèi)萬榮[47]的研究表明，84.2%高原鼠兔喜好將棲息地選擇在坡度在0°～5°的區(qū)域，而在6°～10°和11°～15°的區(qū)域分布分別只占13.2%和2.6%；對坡向的選擇更傾向于南或偏南向；在山地地貌上占總數(shù)的100%，在所有地貌上占41.3%。這有利于增加冬季對太陽輻射的獲取，提升洞內溫度。NDVI能夠在一定程度上反映植被的蓋度和生物量，馬波等[48]的研究認為高原鼠兔洞口數(shù)與NDVI存在顯著正相關的線性關系，高原鼠兔有選擇地利用NDVI較高、植被較好的生境。高原鼠兔的挖掘活動能夠改變當?shù)氐耐寥罓顩r，Guo等[49]的研究認為高原鼠兔的洞穴對土壤pH、含水量、含沙量等土壤理化性質均有影響，能增加土壤滲透率，并提高土壤養(yǎng)分。同時高原鼠兔的分布受降水量和輻射的影響，郭新磊[50]的研究表明，高原鼠兔洞口密度隨降水量和輻射的增高表現(xiàn)為先增大后減少，降水量在600～700 mm區(qū)間和輻射在240～250 J·m-2區(qū)間的地區(qū)高原鼠兔洞口密度達到最大。以上影響高原鼠兔空間分布的因子與前人的研究相一致，說明BIOMOD不僅能很好地的模擬高原鼠兔的空間分布，在影響因子的確定上也有著良好的效果。

圖3 基于RF的黃河源區(qū)高原鼠兔分布預測Figure 3 Prediction of the distribution of plateau pika, based on RF in the source region of the Yellow River Basin

高原鼠兔分布除了與環(huán)境因子有關外，還與人類活動密切相關，如放牧、害鼠防治、招鷹架的布置等。張興祿[51]在瑪曲高寒草甸植被區(qū)所做試驗證明輪牧制度可有效抑制高原鼠兔數(shù)量增長，而在連續(xù)放牧制度下，低放牧率下高原鼠兔的數(shù)量則明顯上升，而高放牧率下由于牛羊踐踏、采食過重，不適合高原鼠兔生存，數(shù)量則有所下降。梁杰榮等[52]在海北高寒草甸生態(tài)系統(tǒng)定位研究站調查了高原鼠兔的種群恢復過程，發(fā)現(xiàn)害鼠防治后殘鼠數(shù)量按Logistics曲線增長；當防治效果高達90%時，其數(shù)量約在兩年內也恢復至原水平。侯秀敏[53]在澤庫巴灘所做試驗表明，招鷹架的布置能有效吸引高原鼠兔的天敵，降低高原鼠兔的密度。當鷹架布置的間距為 500 m × 500 m 時，控制面積為每25 hm2一根，鷹架的利用率達到了98%，樣地內兩年間高原鼠兔的有效洞口數(shù)量降低了26.47%。

致謝：本研究使用的氣候數(shù)據(jù)由青藏高原多圈層模擬與數(shù)據(jù)同化中心 (Data Assimilation and Modeling Center for Tibetan Multi-spheres)開發(fā)，中國科學院青藏高原研究所陽坤研究員提供，特在此表示誠摯的謝意。