張 鑫,尹文萍,謝 菲,樊 輝,*,陳 飛

1 云南大學國際河流與生態(tài)安全研究院,昆明 650091 2 云南省國際河流與跨境生態(tài)安全重點實驗室,昆明 650091 3 國家林業(yè)和草原局亞洲象研究中心,昆明 650216 4 國家林業(yè)和草原局昆明勘察設(shè)計院,昆明 650216

生境退化與破碎化加劇正導致全球生物多樣性的快速喪失[1—3]。生境適宜性評價通過綜合多個與物種相關(guān)的氣候、地形、地表覆被與人為干擾因子,量化物種的現(xiàn)有生境并預測其潛在生境,是物種保護和生境管理與規(guī)劃的基礎(chǔ)[4]。生境適宜性評價應用廣泛,目前,常用的生境適宜性評價方法包括生境適宜性指數(shù)(Habitat Suitability Index,HSI)模型、生態(tài)位因子分析(Ecological Niche Factor Analysis,ENFA)模型、最大熵(Maximum Entropy,MaxEnt)模型、生境質(zhì)量(Integrated Valuation of Ecosystem Services and Trade-offs,InVEST)模型等[5—8]。以往的生境適宜性研究大多根據(jù)經(jīng)驗選擇生境因子,依賴于物種出現(xiàn)點數(shù)據(jù)賦予相應權(quán)重來衡量生境適宜性,遙感大數(shù)據(jù)在生境評價中的作用尚未得到充分挖掘。近年來,遙感大數(shù)據(jù)制圖的迅速發(fā)展為生境適宜性評價研究提供了新的數(shù)據(jù)源[9],如土地覆被、植被、氣候、水體、人類干擾等遙感數(shù)據(jù)產(chǎn)品對目標地物較為全面、一致的多元信息描述[10],為更準確、精細地辨識大區(qū)域范圍內(nèi)物種適宜生境提供了可能。因此,構(gòu)建遙感大數(shù)據(jù)驅(qū)動的生境適宜性評價體系,對物種生境做出快速、準確的評價,是當前亟待解決的問題,對野生動物保護、生物多樣性維持以及野生動物突發(fā)事件應急響應具有重要意義。

亞洲象(Elephasmaximus)是生態(tài)系統(tǒng)的旗艦物種,在生態(tài)系統(tǒng)中發(fā)揮著重要作用[11—12]。我國亞洲象主要分布在云南省南部西雙版納州、普洱市和臨滄市等地[13]。但自2020年3月以來,原棲息于西雙版納州國家級自然保護區(qū)勐養(yǎng)子保護區(qū)的“北移象群”從保護區(qū)出走,近數(shù)十年來首次進入元江流域,并持續(xù)向元江流域上游方向遷移。我國亞洲象家域范圍生境適宜性已有廣泛研究[8, 14—15],但因亞洲象數(shù)十年未曾在鄰近區(qū)域出現(xiàn),故鮮有研究涉及亞洲象家域范圍外的生境適宜性問題。此次“北移象群”突然遷入元江流域,該流域內(nèi)是否存在適宜生境,其空間格局如何？是否能支撐該象群的長期生存？遷移路線是否為最適生境等問題尚不清楚,難以支撐遷移過程中亞洲象的保護和風險防范應急等。因此,本文基于遙感大數(shù)據(jù)驅(qū)動的生境適宜性指數(shù)(HSI)模型,采用薈萃分析挖掘并篩選出對亞洲象分布具有重要影響的生境因子,定量評價并對比分析元江-李仙江流域亞洲象生境適宜性格局,以期為遷移中亞洲象群的保護與風險防范應急調(diào)控、亞洲象后備生境規(guī)劃等提供科學支撐。

1 研究區(qū)概況

元江-李仙江流域位于東經(jīng)100°06′—103°53′,北緯22°27′—25°32′,處于云南省中部、東南部,海拔范圍76—3344 m(圖1)。流域總面積57376.39 km2,其中元江、李仙江流域面積分別為23611.39 km2(占41.15%)和33764.00 km2(占58.85%)。流域涵蓋中部高原溫和區(qū)、元江炎熱地區(qū)和熱帶季風林區(qū)三個氣候區(qū)[21],年降水量731—2307 mm[22],年均溫度14.8—23.8℃[23]。流域內(nèi)土地利用/土地覆被類型主要為林地、耕地、草地和灌木地,而水體、濕地、建設(shè)用地、橡膠林和裸地分布較少,其中橡膠林主要分布在元江、李仙江和藤條江下游。

圖1 研究區(qū)位置及地形、亞洲象出現(xiàn)點和亞洲象活動路線Fig.1 Location and topography of the study area, Asian elephant occurrence points and Asian elephant migration routes國界線來源于中國科學院資源環(huán)境科學與數(shù)據(jù)中心https://www.resdc.cn/；亞洲象出現(xiàn)點來自新聞媒體報道[16—20],活動路線依據(jù)出現(xiàn)點推測

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

2.1 數(shù)據(jù)來源與處理

使用的數(shù)據(jù)包括：研究區(qū)2020年土地利用/土地覆被、地形、植被類型、歸一化差異植被指數(shù)、降水量、河流、道路等數(shù)據(jù)。2020年土地利用/土地覆被數(shù)據(jù)來源于http://www.globallandcover.com/,空間分辨率為30 m,由于該分類系統(tǒng)沒有對橡膠林單獨分類,本文基于2019年云量小于15%的地表反射率OLI影像,采用基于橡膠林物候的歸一化燃燒指數(shù)(Normalized Burn Ratio,NBR)差值法提取了流域范圍內(nèi)的橡膠林分布[24—25],然后依據(jù)亞洲象對不同土地利用類型的偏好(最小值為0,最大值為10)重新分類土地利用/土地覆被類型[26]；地形數(shù)據(jù)來源于http://gdex.cr.usgs.gov/gdex/航天雷達地形測繪任務(wù)的數(shù)字高程模型(Shuttle Radar Topography Mission Digital Elevation Model,SRTM DEM),空間分辨率為30 m；植被類型數(shù)據(jù)來自中國科學院資源環(huán)境科學與數(shù)據(jù)中心(https://www.resdc.cn/data.aspx?DATAID=122),依據(jù)亞洲象對不同植被類型的偏好作重分類[26]；歸一化差異植被指數(shù)數(shù)據(jù)來自https://developers.google.com/earth-engine/datasets/catalog/,取2013—2018年逐年最大值合成后時序NDVI的中值；降水量數(shù)據(jù)來自https://www.worldclim.org/,為1970—2000年降水量平均值；河流數(shù)據(jù)提取自云南省1∶25萬地形圖；道路數(shù)據(jù)來自全國地理信息資源目錄服務(wù)系統(tǒng)(https://www.webmap.cn/main.do?method=index)。對上述所有數(shù)據(jù)進行地理編碼后,利用最鄰近法(Nearest Neighbor)重采樣至30 m空間分辨率,并進行相應處理(表1)。

表1 生境評價因子及處理方式Table 1 Habitat evaluation factors and their processing methods

1.2 研究方法

1.2.1生境評價因子選擇

評價因子選擇對于模型運行結(jié)果準確性至關(guān)重要[27]。亞洲象對生境的選擇利用主要受氣候、地形、地表覆被、干擾等因素的影響[28—29]。本文利用Web of Science和CNKI數(shù)據(jù)庫檢索出亞洲象生境適宜性相關(guān)文獻。其中,Web of Science的文獻檢索式為：(TS="habitat suitability" OR TS="habitat use" OR TS="habitat index" OR TS="habitat select") AND (TS="Asian elephant" OR TS="wild elephant"),CNKI數(shù)據(jù)庫的文獻檢索式為：(生境評價+生境適宜性+生境選擇+棲息地適宜性+棲息地評價+棲息地選擇)*亞洲象。根據(jù)檢索到的文獻統(tǒng)計生境評價因子詞頻,共檢索到亞洲象生境適宜性評價文獻37篇,涉及生境適宜性評價因子142個。參考亞洲象習性,依據(jù)統(tǒng)計結(jié)果選取出現(xiàn)頻率較高的生境評價因子。

為了消除評價因子間的多重共線性對模型預測結(jié)果的不利影響[30],采用Pearson相關(guān)系數(shù)(Pearson Correlation Coefficient)和方差膨脹因子(Variance Inflation Factor, VIF)獨立性檢驗篩選生境評價因子[31]。當評價因子間的相關(guān)性大于0.7時,變量間存在共線性[32]；當評價因子的方差膨脹因子大于10時,變量間存在多重共線性[33]。在計算VIF時,逐步剔除VIF值最高的評價因子,直到剩余評價因子的VIF值均小于10。

2.2.2生境適宜性模型

通過計算HSI[34]衡量亞洲象生境適宜性。計算公式為：

(1)

式中,fi為第i個評價因子,wi為第i個因子的權(quán)重,n為評價因子的數(shù)量。評價因子權(quán)重計算有多種方法[34—36],本文綜合層次分析法(Analytic Hierarchy Process,AHP)和主成分分析法(Principal Components Analysis,PCA)兩種方法,對兩者的計算結(jié)果加權(quán)平均,得到綜合權(quán)重。AHP是一種多準則的決策分析方法[37],通過比較兩個影響因素之間的相對重要性來構(gòu)建矩陣[8, 15, 38—39],再由矩陣的最大特征根及其對應的特征向量,計算出各因素的權(quán)重；PCA首先將原變量重新組合成主成分,然后根據(jù)原變量的總方差和主成分的特征值計算原變量的權(quán)重[35]。加權(quán)平均上述兩種方法所得權(quán)重,既能克服AHP的主觀影響,也能有效削減指標因子間相關(guān)性的影響,從而提高模型結(jié)果的可靠性[40]。

采用專家經(jīng)驗法重分類生境適宜性結(jié)果,將亞洲象生境分為最適生境(0.75—1)、較適生境(0.6—0.75)、邊際生境(0.5—0.6)和非生境(0—0.5)4個等級[8, 15, 39],并計算不同緯度平均生境適宜性和2021年“北移象群”活動路線沿程生境適宜性及其趨勢。

2.2.3生境景觀格局分析

采用斑塊數(shù)(Number of Patches, NP)、斑塊密度(Patch Density, PD)、最大斑塊指數(shù)(Largest Patch Index, LPI)、平均斑塊面積(Mean Patch Area, AREA_MN)、斑塊分維數(shù)(Fractal Dimension Index, FRAC)和斑塊連接度指數(shù)(Connectance Index, CONNECT)6個指數(shù)來表征亞洲象生境景觀格局特征,這些指數(shù)由FRAGSTATS 4.2(https://www.umass.edu/landeco/research/fragstats/downloads/fragstats_downloads.html)計算得出。其中,NP表示景觀中某種類型斑塊的數(shù)量；PD表示景觀中每平方公里的斑塊數(shù)量,代表景觀整體斑塊分化程度；LPI表示某種類型的最大斑塊面積占整個景觀的總面積的百分比；AREA_MN表示所有某種類型斑塊的平均斑塊的大小；FRAC取值范圍為[1,2],值越接近1,表示斑塊形狀越接近歐幾里得幾何形狀[39]；CONNECT取值范圍[0,100],0表示閾值距離內(nèi)生境間不存在連接狀態(tài)、生境破碎化程度高,100代表閾值距離內(nèi)所有生境處于連接狀態(tài)。

3 結(jié)果與分析

3.1 元江-李仙江流域亞洲象生境格局

圖2表明,依據(jù)使用頻次高低排序,初步篩選出海拔、坡度、到水源距離、到居住點的距離、土地利用/土地覆被等17個亞洲象生境適宜性評價因子。

圖2 亞洲象生境適宜性評價因子使用頻次Fig.2 Frequency of Habitat suitability evaluation factors

由圖3和圖4可見,經(jīng)相關(guān)性分析和方差膨脹因子檢驗篩選,森林頻率和地表粗糙度兩個因子被剔除,最終得到15個亞洲象生境適宜性評價因子。圖5為基于開源遙感數(shù)據(jù)產(chǎn)品量化出的各生境適宜性評價因子空間格局,各因子的空間相似性較低；圖6顯示了結(jié)合AHP和PCA所得各生境適宜性評價因子的權(quán)重,其中土地利用/土地覆被的權(quán)重最高(0.1708),灌草叢頻率的權(quán)重最低(0.0159)。

圖3 亞洲象生境適宜性評價因子相關(guān)性Fig.3 Correlation of habitat suitability evaluation factors

圖4 亞洲象生境適宜性評價因子方差膨脹因子：剔除森林頻率和地表粗糙度前后Fig.4 Variance inflation factor of habitat suitability evaluation factors：Before and after removing the factors of forest frequency and surface roughness

圖5 亞洲象生境適宜性評價因子空間格局Fig.5 Spatial patterns of habitat suitability evaluation factors for Asian elephants

圖6 亞洲象生境適宜性評價因子權(quán)重Fig.6 Weights of habitat suitability evaluation factors for Asian elephants

由圖7可知,研究區(qū)內(nèi)亞洲象最適生境面積為3271.27 km2(5.70%),較適和邊際生境面積分別為19278.97 km2(33.62%)和16360.20 km2(28.52%)。分流域來看,元江流域內(nèi)亞洲象的最適生境面積為590.78 km2(1.75%),主要分布于元江流域下游中-越邊境附近；較適生境和邊際生境面積分別為8175.83 km2(24.21%)和12412.95 km2(36.76%)。李仙江流域內(nèi)亞洲象最適生境面積為2680.49 km2(11.35%),主要分布于李仙江流域下游段,且相對集中,流域中游段僅有少量最適生境分布；較適生境和邊際生境面積分別為11112.15 km2(47.06%)和3947.25 km2(16.72%)。圖7進一步表明,隨緯度升高,由南到北亞洲象生境適宜性總體上呈逐漸下降趨勢。2021年“北移象群”活動路線由西南向東北橫貫流域中游,沿程生境適宜性亦呈下降趨勢。

圖7 元江-李仙江流域亞洲象生境分布、不同緯度平均生境適宜性及生境適宜性隨活動路線變化Fig.7 The habitat distribution of Asian elephants in the Yuanjiang-Lixianjiang River Basin, the average habitat suitability along the latitude gradients and migration routes

3.2 生境景觀格局

由圖8可知,元江流域內(nèi)三種亞洲象生境類型的NP、LPI、AREA表現(xiàn)為邊際生境>較適生境>最適生境,PD和CONNECT表現(xiàn)為最適生境>較適生境>邊際生境；而李仙江流域內(nèi)三種亞洲象生境類型的NP和PD表現(xiàn)為邊際生境>最適生境>較適生境,LPI和AREA表現(xiàn)為較適生境>最適生境>邊際生境。元江流域亞洲象最適生境的NP、LPI、AREA和FRAC均小于李仙江流域,而PD和CONNECT高于李仙江流域。

圖8 元江-李仙江流域亞洲象生境景觀格局Fig.8 Landscape indices of Asian elephants′ habitats in the Yuanjiang-Lixian River basin

4 討論與結(jié)論

4.1 討論

(1)遙感大數(shù)據(jù)驅(qū)動的亞洲象生境適宜性評價：遙感大數(shù)據(jù)制圖的迅速發(fā)展為大尺度生境研究提供了豐富的數(shù)據(jù)源[9, 41—42]?；谝延衼喼尴笊尺m宜性研究,采用薈萃分析挖掘并篩選出使用頻率較高的生境因子,有效避免了生境因子選取的主觀性干擾,并利用具有較好一致性、較高空間分辨率與精度且開源的遙感數(shù)據(jù)產(chǎn)品量化生境適宜性評價因子,構(gòu)建了一套遙感大數(shù)據(jù)驅(qū)動的亞洲象生境適宜性評價體系,快速評價了元江-李仙江流域的亞洲象生境適宜性。與亞洲象生境適宜性評價常用的ENFA模型和MaxEnt模型相比[5, 39],該方法不依賴亞洲象出現(xiàn)點數(shù)據(jù),能有效規(guī)避人為景觀(如耕地、居民點和城鎮(zhèn)等)中的出現(xiàn)點導致的適宜生境誤判[15, 27]以及出現(xiàn)點觀測數(shù)據(jù)獲取周期長、獲取難度大等問題。土地利用/土地覆被是影響亞洲象生境適宜性的重要因子。橡膠林種植作為影響我國亞洲象生境適宜性的主要土地利用/土地覆被類型[39],在現(xiàn)有全球尺度土地利用/土地覆被數(shù)據(jù)產(chǎn)品(如GlobeLand30)多被劃歸為林地[9]。針對此,本研究基于2019年Landsat 8 OLI合成影像,采用基于橡膠林獨特物候特征的NBR差值法提取研究區(qū)橡膠林種植范圍,總體精度達86.05%,為流域亞洲象適宜生境評價提供更為準確的數(shù)據(jù)支撐。此外,現(xiàn)有全球尺度土地利用/土地覆被數(shù)據(jù)集存在分類精度偏低、產(chǎn)品間空間一致性較差等問題[43—44]。但隨著土地利用/土地覆被遙感制圖技術(shù)的發(fā)展,全球尺度土地利用/土地覆被產(chǎn)品將愈加豐富[9],未來可整合數(shù)據(jù)融合等技術(shù),充分發(fā)揮多源遙感數(shù)據(jù)生產(chǎn)的土地利用/土地覆被數(shù)據(jù)集的各自優(yōu)勢特點,提高數(shù)據(jù)集精度[43],更好地服務(wù)于亞洲象生境適宜性評價。遙感大數(shù)據(jù)驅(qū)動的適宜生境評價不僅極大地拓寬了傳統(tǒng)生境適宜性評價的數(shù)據(jù)來源,提升了大尺度適宜生境制圖的能力[45],且具有較好的時效性,在野生動物突發(fā)事件應急處置與調(diào)控中極具潛力。而在本次亞洲象遷移應急響應中,主要依賴于無人機遙感技術(shù)監(jiān)測象群動態(tài)[46],而高空間分辨率、近實時的“空-天”遙感未能得到有效應用。未來迫切需要健全和完善現(xiàn)有地面野外觀測網(wǎng)絡(luò)[47],整合高分辨率國產(chǎn)資源衛(wèi)星數(shù)據(jù)[48]、高分數(shù)據(jù)[49],以及無人機遙感[50—51]數(shù)據(jù)等,建立近實時“空-天-地”一體化的亞洲象監(jiān)測與突發(fā)事件應急響應平臺[52—53],為亞洲象保護與應急管理及防控提供更全面、實時、精準的決策支持。

(2)流域亞洲象適宜生境格局及其種群承載力：元江流域面積(23611.39 km2)遠小于李仙江流域(33764.00 km2),但其較適和最適生境的NP、PD高于李仙江流域,而AREA卻小于李仙江流域。因此,較于李仙江流域,元江流域亞洲象生境破碎化程度更高。兩個流域FRAC均接近于1,CONNECT均處于較低水平,表明兩個流域生境均曾受到一定程度的人類活動干擾[8, 39],而元江流域較適生境和最適生境LPI小于李仙江流域,表明前者受到人類活動的干擾更為強烈[54]。研究區(qū)內(nèi)最適生境斑塊主要位于李仙江流域下游段,其最適生境中最大斑塊面積為629.93 km2,而元江流域最適生境中最大斑塊面積為僅為56.94 km2,遠小于亞洲象種群最小家域面積(105 km2)[55],表明元江流域生境無法支撐該象群的長期生存與種群維持。隨著該象群向流域上游段遷移,受海拔、降水等因素的疊合影響,天然闊葉林分布減少,針葉林和灌草叢增多,象群生存所需的隱蔽場所和食物來源更為匱乏[8],這可能是該象群進入元江流域后遷移速度加快的主要原因。目前,由于該象群逐漸遠離亞洲象原棲息地,若長期滯留,極可能形成極小種群效應[56],對其遺傳多樣性產(chǎn)生威脅。此外,元江流域上游段人口密度較大,遠高于我國亞洲象生存可耐受的20人/km2最大人口壓力閾值[57]?！氨币葡笕骸钡耐蝗贿w入,已對當?shù)鼐用竦纳c財產(chǎn)安全構(gòu)成極大威脅[58—59],人象沖突潛在風險劇增。為規(guī)避此風險,亟需結(jié)合亞洲象生境適宜性、亞洲象行為生物學機制與規(guī)律[60],盡快引導象群回歸適宜生境區(qū)。本研究僅分析了元江-李仙江流域亞洲象生境現(xiàn)狀,未回溯該流域亞洲象適宜生境歷史格局及變化動態(tài)。因我國亞洲象種群數(shù)量持續(xù)增長對適宜生境需求勢必上升,未來可結(jié)合長時序遙感大數(shù)據(jù),開展對亞洲象適宜生境變化時空過程及驅(qū)動機制研究[39],為我國亞洲象生境恢復及后備生境創(chuàng)建提供科學支撐。

4.2 結(jié)論

采用薈萃分析挖掘并篩選出亞洲象生境評價因子,結(jié)合相關(guān)性分析和VIF等方法診斷生境適宜性評價因子獨立性,并基于開源遙感數(shù)據(jù)產(chǎn)品量化生境因子,構(gòu)建HSI模型,探究了元江-李仙江流域亞洲象生境適宜性及其景觀格局,得出以下結(jié)論：

(1)元江-李仙江流域亞洲象生境適宜性主要表現(xiàn)為從下游到上游逐漸降低,最適生境較少且主要分布于流域下游段。

(2)元江流域亞洲象生境適宜性低于李仙江流域,最適生境分布更少且破碎化更嚴重,無法滿足亞洲象種群最小家域面積需求,難以獨立維持“北移象群”的長期生存。

(3)沿“北移象群”遷移活動路線,土地利用/覆被類型差異較大,天然闊葉林沿程分布減少,生境適宜性呈下降趨勢。