李志明，王 乾,余國(guó)寶，陳緒文*

(1.四川貢嘎山國(guó)家級(jí)自然保護(hù)區(qū)石棉管理處，四川 石棉 625400;2.中國(guó)科學(xué)院成都生物研究所，四川 成都 610041)

生物多樣性是區(qū)域性生態(tài)—進(jìn)化過程的外在表現(xiàn)，是生態(tài)系統(tǒng)功能的基礎(chǔ)，是維系人類社會(huì)生存和發(fā)展的根本[1]。目前全球生物多樣性在以前所未有的速度消失，即使作為旗艦性保護(hù)物種的大熊貓，也受到棲息地生境破碎化的不良影響。四川盆地西緣山地是大熊貓主要分布區(qū)域，也是箭竹類竹種的分化中心和現(xiàn)代分布中心[4]，大熊貓保護(hù)很大程度上得益于這里適宜的生境條件和竹類資源[2]。箭竹類作為一種林下和灌叢中的木本植物資源，不僅是當(dāng)?shù)厣锒鄻有缘闹匾M成部分，還是重要的經(jīng)濟(jì)植物，其種類和生存狀態(tài)對(duì)生態(tài)系統(tǒng)平衡有著重要的作用[3]。竹類主要靠無性克隆繁殖，營(yíng)養(yǎng)期長(zhǎng)達(dá)幾十年，擴(kuò)散能力較差[5]。而區(qū)域內(nèi)竹子一旦開花多會(huì)同步進(jìn)行，并隨之造成大片竹林死亡，導(dǎo)致大熊貓出現(xiàn)食物短缺甚至死亡。因此竹類分布對(duì)環(huán)境反應(yīng)具有滯后性和周期性，環(huán)境波動(dòng)過大容易導(dǎo)致局部滅絕。

氣候是決定物種分布的頂層驅(qū)動(dòng)力，其次為生境類型和生態(tài)過程，尺度從大到小通過層層篩選，物種的實(shí)際分布范圍要小于氣候適宜區(qū)[9]。熊貓保護(hù)區(qū)天然竹類主要分布在林下，從上層林冠類型到地面基質(zhì)狀況，從種間競(jìng)爭(zhēng)到演替階段以及干擾歷史無不對(duì)其分布起著約束作用。氣候變暖速率超過物種擴(kuò)散速率也會(huì)造成竹類生境范圍和資源的大幅減少[6～8]。因此詳細(xì)研究竹類分布與氣候和環(huán)境的關(guān)系特征對(duì)于了解竹類分布狀況，掌握承載力，及早發(fā)現(xiàn)敏感區(qū)域，開展科學(xué)合理的保護(hù)措施有重要意義。目前關(guān)于竹類的空間分布研究?jī)H有少數(shù)描述性的資料發(fā)表，沒有具體分布位點(diǎn)，不能確定相應(yīng)的環(huán)境因素，或者涉及保密只給出整個(gè)熊貓保護(hù)區(qū)的分析結(jié)果[7]，這種結(jié)果涉及范圍大，氣候、植被等環(huán)境因素的空間異質(zhì)性較強(qiáng)，難以反映物種的環(huán)境需求。

本文應(yīng)用四川省第1高山——貢嘎山來研究竹類資源分布與生態(tài)模型尺度的關(guān)系。采用貢嘎山竹類調(diào)查分布位點(diǎn)和WorldClim生物氣候數(shù)據(jù)，利用集合模型構(gòu)建了貢嘎山區(qū)3種竹類的分布預(yù)測(cè)，以期為保護(hù)區(qū)分區(qū)規(guī)劃和生物多樣性保護(hù)工作的依據(jù)。

1 研究區(qū)域

貢嘎山是四川省第1高山，也是青藏高原東部的最高山，整個(gè)山區(qū)東部為高山峽谷地貌，向西則與川西高原相接，東西氣候差異顯著[10]，動(dòng)植物資源極為豐富，植被類型多樣，空間分布格局顯著[11]。河谷地帶為亞熱帶干旱河谷氣候區(qū)，山頂則為冰雪永凍帶，東坡山勢(shì)陡峭，降水豐富，植被茂密，植被垂直帶譜豐富，類型多樣，物種豐富。山體整體處于中亞熱帶向北亞熱帶過渡區(qū)域，基帶性植被基本穩(wěn)定，氣候梯度由山體地形地貌格局造成，南北環(huán)境差異不顯著，由于整體范圍不大，可以認(rèn)為環(huán)境因素的響應(yīng)在尺度上是不變的，因此適合分析物種的環(huán)境需求。

貢嘎山西部為山地向高原過渡的山原地貌區(qū)，氣候干冷、輻射強(qiáng)烈，植被以亞高山針葉林和高山灌草叢為主。貢嘎山自然保護(hù)區(qū)環(huán)繞貢嘎山主峰及其周邊山峰，包括石棉、瀘定、康定和九龍4個(gè)縣的部分地區(qū)，總面積約 4 091 km2。貢嘎山位于橫斷山脈東部，在地理環(huán)境上為水熱轉(zhuǎn)換地帶，是東西向和南北向上的生物區(qū)系過渡地帶，具有重要的科研價(jià)值。貢嘎山同時(shí)也是大熊貓3個(gè)種群(小相嶺、大相嶺和邛崍山)的連接區(qū)域，對(duì)保持熊貓遺傳多樣性的完整性起著關(guān)鍵作用。而竹類作為熊貓的食物在廊道的質(zhì)量上的作用尤為重要。我們2017年的調(diào)查發(fā)現(xiàn)整個(gè)分布區(qū)有竹子種類8種，其中天然竹類6種，包括箭竹屬3種，巴山木竹屬兩種，玉山竹屬1種。我們選用其中分布最廣的3個(gè)竹種做空間分布研究，即豐實(shí)箭竹(Fargesiaferax)、冷箭竹(Bashaniafangiana)和峨熱竹(Bashaniaspanostachya)。

2 材料與方法

竹類分布位點(diǎn)數(shù)據(jù)來自于帶GPS定位功能的相機(jī)拍攝的竹類微距照片，照片拍攝于2017年4月和9月的貢嘎山竹類調(diào)查工作中。共拍攝兩萬余張照片，從中選出1500個(gè)竹類位點(diǎn)。去除重疊的位點(diǎn)后，共篩選出320個(gè)竹種分布點(diǎn)，各竹種的分布位點(diǎn)在預(yù)測(cè)圖中標(biāo)出，以資參考。調(diào)查區(qū)域基本涵蓋貢嘎山區(qū)可通行的區(qū)域。

氣候數(shù)據(jù)來自于WorldClim數(shù)據(jù)庫[12]的分辨率約1 km的空間插值生成的生物氣象數(shù)據(jù)，共包括19個(gè)變量，由于數(shù)據(jù)之間可能存在共線性問題，造成模型預(yù)測(cè)效果較差，因此需要事先將共線性較強(qiáng)的數(shù)據(jù)剔除。依據(jù)成對(duì)Pearson相關(guān)系數(shù)矩陣對(duì)這19個(gè)生物氣候變量進(jìn)行選擇，去除平均相關(guān)系數(shù)絕對(duì)值大于0.7的變量[7]，最后保留6個(gè)變量，即bio2、bio3、bio4、bio12、bio15、bio17變量，其中前3個(gè)變量與溫度有關(guān)，后3個(gè)與降水有關(guān)。所有分析使用biomod2包[13]在CRAN R軟件中進(jìn)行。Biomod2包中包括了多種模型分析，其中集合模型(Ensemble Modeling)將模型中的10個(gè)模型(GLM,GBM,GAM,CTA,ANN,SRE,FDA,MARS,RF,MAXENT)作為子模型進(jìn)行物種分布建模，然后利用投票方式，以票數(shù)過半來決定物種的在每個(gè)空間子區(qū)上的分布情況。這10個(gè)模型包括廣義線性模型、廣義增強(qiáng)模型、廣義加性模型、分類樹分析、人工神經(jīng)網(wǎng)絡(luò)、生物氣候模型、柔性判別分析、多重適應(yīng)性回歸樣條、隨機(jī)森林以及最大熵模型等。這些模型各自在某一方面表現(xiàn)較強(qiáng)，綜合的結(jié)果通常較單一模型表現(xiàn)更優(yōu)。

調(diào)查的竹類分布位點(diǎn)數(shù)據(jù)被分成兩部分，其中80%用于模型訓(xùn)練，其余部分用于模型質(zhì)量評(píng)估，采用5折交叉驗(yàn)證，計(jì)算模型預(yù)測(cè)精度、敏感性和專一性，TSS(true skill statitic)[14]進(jìn)行模型比較。氣候數(shù)據(jù)為柵格數(shù)據(jù)，空間分辨率30″，在貢嘎山區(qū)約890m，坐標(biāo)系為WGS84經(jīng)緯度格式。通過竹類分布位點(diǎn)計(jì)算對(duì)應(yīng)柵格點(diǎn)上的各氣象要素的值分布情況，然后應(yīng)用各種概率模型計(jì)算最優(yōu)的反應(yīng)曲線。最后反演到整個(gè)區(qū)域，找到相應(yīng)的分布范圍。

3 結(jié)果

單個(gè)模型中以廣義線性模型表現(xiàn)最差，隨機(jī)森林模型表現(xiàn)最好，TSS值大于0.85，說明模型的敏感性和專一性都較高。集合模型通過對(duì)10種模型的預(yù)測(cè)投票，以少數(shù)服從多數(shù)的方法，博采眾家之長(zhǎng)，得到預(yù)測(cè)結(jié)果，TSS值均大于0.90，表現(xiàn)效果超出單個(gè)模型的最好效果。

預(yù)測(cè)的竹類基本沿著溝谷分布，包含著實(shí)際分布點(diǎn)位，各竹種之間呈現(xiàn)一定的空間分布格局(圖1)。豐實(shí)箭竹分布海拔較低，集中分布在瀘定縣磨西鎮(zhèn)和石棉縣草科鄉(xiāng)地段，南部九龍縣洪壩鄉(xiāng)也有一定的分布?？刀h城附近有零星的分布。沿貢嘎山主峰經(jīng)線以西區(qū)域均無豐實(shí)箭竹分布。冷箭竹分布區(qū)域較為狹窄，主要分布在磨西鎮(zhèn)和康定縣城附近海拔較高的有林地段，從康定折多塘至貢嘎山主峰一線以西地段均無冷箭竹分布。峨熱竹分布較為廣泛，除西南部和東部海拔較低處外，基本都有分布，其中草科鄉(xiāng)以西的溝谷中、康定縣城附近、洪壩鄉(xiāng)附近分布較多，在保護(hù)區(qū)之外康定縣朋布西和六巴的東部地區(qū)分布較廣。

圖1 使用集體模型構(gòu)建的貢嘎山區(qū)3種竹類的空間分布預(yù)測(cè)圖(A為豐實(shí)箭竹，B為冷箭竹，C為峨熱竹)，白色方塊構(gòu)成的區(qū)域?yàn)轭A(yù)測(cè)分布區(qū)，紅色叉為實(shí)際調(diào)查分布位點(diǎn)，灰白色線包圍區(qū)域?yàn)樨暩律阶匀槐Ｗo(hù)區(qū)范圍，圖片背景為海拔高程圖。Fig.1 The predicted distribution map of three bamboo species using the ensemble model(A:Fargesia ferax; B:Bashania fangiana; C:Bashania spanostachya).The area composed of white squares are predicted distribution area of the specific species.The crosses are actual present records of the species.The natural reserve of the Gongga Mountains are surrounded by the pale white lines.The background is a digital elevation map.

實(shí)際調(diào)查竹類的海拔分布以豐實(shí)箭竹最低，中位數(shù)海拔為 2 200 m；其次為峨熱竹，中位數(shù)海拔為 3 100 m；最高的是冷箭竹，中位數(shù)海拔為 3 400 m(圖2)。3個(gè)竹種中均存在海拔分布特別低的異常值，但竹種分布海拔之間t統(tǒng)計(jì)檢驗(yàn)存在顯著的區(qū)別 (p<0.001)。

圖2 調(diào)查竹種的實(shí)際海拔分布箱形圖。Fig.2 The boxplot of elevational distribution of the surveyed bamboo species.

但是預(yù)測(cè)區(qū)域的海拔分布與實(shí)際調(diào)查分布有一定的區(qū)別，例如豐實(shí)箭竹的海拔中位數(shù)為 2 340 m，比實(shí)際提高140 m。預(yù)測(cè)的海拔最低值 1 050 m，最高值達(dá)到 3 770 m，遠(yuǎn)遠(yuǎn)超過實(shí)際調(diào)查的海拔高度范圍。峨熱竹預(yù)測(cè)區(qū)的海拔中位數(shù)為 3 910 m，比調(diào)查的中位數(shù)提高810 m，海拔幅度在 1 740 m～4 630 m之間，也遠(yuǎn)超調(diào)查的范圍。冷箭竹預(yù)測(cè)區(qū)海拔中位數(shù)為 3 320 m，與實(shí)際調(diào)查的中位數(shù)最為接近，但是預(yù)測(cè)的海拔范圍(2 170 m～4 630 m)與實(shí)際調(diào)查情況差異也很大。

在保護(hù)區(qū)內(nèi)，各竹種的預(yù)測(cè)分布面積以峨熱竹最大，達(dá)到231.7 km2，其次為豐實(shí)箭竹，達(dá)到170.6 km2，冷箭竹面積最小，為134.9 km2。保護(hù)區(qū)中這3種竹類的分布面積只占13%。

4 討論

4.1 氣候異質(zhì)性

貢嘎山東西坡降水呈現(xiàn)不同的模式[10,11]，東坡受東南季風(fēng)影響強(qiáng)烈，隨著海拔升高，降水量逐步增多。由于貢嘎山地形抬升，形成地形雨，在迎風(fēng)坡的一定海拔范圍內(nèi)降水隨海拔增加而增加。其增加幅度從瀘定縣城到磨西鎮(zhèn)為每百米131 mm，從磨西鎮(zhèn)到海螺溝索道(海拔 3 000 m)為67 mm，總降水量達(dá)到約 2 000 mm[15]。貢嘎山西部地區(qū)受西南季風(fēng)和東南季風(fēng)的雙重影響，降水模式為雙峰形式，最大降水量約出現(xiàn)在6月和9月。貢嘎山西南部的九龍縣城降水量雖大于貢嘎山北部的康定縣城，但也只有同海拔的東部區(qū)域的一半。除了降水之外，同海拔地區(qū)的溫度也有顯著的區(qū)別，東部顯著低于西部地區(qū)。海螺溝年均溫4.2℃，九龍縣為10.7℃。貢嘎山東坡林線位置多在海拔 3 600 m左右，而九龍縣湯古瓦灰山林線海拔在 4 300 m左右，兩者相差700 m。而林線的位置被認(rèn)為具有相同的最熱月溫度或者生長(zhǎng)季溫度[16]，因此東部氣溫的垂直遞減速率顯著大于西部。

我們用相關(guān)系數(shù)法去除了大部分氣候變量，只保留了共線性較弱的6個(gè)變量，這樣可以增強(qiáng)模型的穩(wěn)健性和移植性，對(duì)附近區(qū)域的竹類分布預(yù)測(cè)有指示意義。例如我們發(fā)現(xiàn)康定朋布西和六巴以東區(qū)域預(yù)測(cè)有大片的峨熱竹分布，將來的實(shí)地探查可以驗(yàn)證模型的可靠性。

這種水熱分布上的顯著區(qū)別在山勢(shì)陡峭、地形復(fù)雜的東部區(qū)域往往造成復(fù)雜的地形氣候模式，用分辨率粗糙的氣候模型模擬會(huì)掩蓋許多特殊氣候模式。目前用于預(yù)測(cè)未來氣候變化的數(shù)據(jù)空間分辨率較為粗糙，基本采用WorldClim的約1km精度數(shù)據(jù)，其精度受限于站點(diǎn)密度和地形特征[10]，這就造成了顯著的取樣偏差。在平原和丘陵地區(qū)使用不會(huì)造成太大的誤差，但是在川西高山峽谷區(qū)由于氣象站點(diǎn)稀少、地形起伏劇烈以及地形氣候效應(yīng)較為突出的地方，使用粗分辨率數(shù)據(jù)往往造成較大的誤差。例如海拔每升高1 km，溫度大約下降6.5℃[16]，按照45°的山坡坡度計(jì)算，1km水平分辨率的溫度數(shù)據(jù)在峽谷區(qū)的誤差就已經(jīng)超過了IPCC估計(jì)的未來幾十年的溫度升高幅度。在這種尺度下的模擬研究可能會(huì)夸大變化范圍，誤導(dǎo)生物多樣性的研究和保護(hù)。而且這種預(yù)測(cè)區(qū)域所覆蓋地形海拔不可避免會(huì)超過實(shí)際的分布海拔，其預(yù)測(cè)的可信度降低。因此粗分辨率數(shù)據(jù)只能給出平均的環(huán)境特征，不能反映鑲嵌分布的面積較小的特殊生境特征。

在物種分布的海拔預(yù)測(cè)上，使用粗分辨率的氣候數(shù)據(jù)導(dǎo)致海拔出現(xiàn)較大偏差，特別是在地形陡峭的地段，預(yù)測(cè)范圍往往超出竹類可能分布的海拔上限。因此粗分辨率數(shù)據(jù)只能估計(jì)平原等與平均狀態(tài)較為接近的區(qū)域，在所研究的竹類中冷箭竹和豐實(shí)箭竹的實(shí)際海拔分布的中位數(shù)與預(yù)測(cè)海拔的中位數(shù)分布最接近。但是和峨熱竹的海拔中位數(shù)差異達(dá)到810 m，這種情況可能和西北部大片出現(xiàn)的新預(yù)測(cè)分布區(qū)的海拔較高有很大關(guān)系。

4.2 物種特征

竹子的分布和植被也有一定的關(guān)系。竹子通常分布在稀疏的林下或者林窗內(nèi)，要求一定的光照，因此其分布狀況和森林的演替階段有一定的關(guān)聯(lián)。密林下或者人工林下由于林內(nèi)光線較弱，竹類通常較難生存。冷箭竹和峨熱竹生長(zhǎng)在高海拔的針闊混交林內(nèi)，特別是樺木云冷杉林下，竹子常是灌木層的優(yōu)勢(shì)種，呈散生狀態(tài)，高度可達(dá)2 m左右。豐實(shí)箭竹常生長(zhǎng)在針闊混交林和落葉闊葉林內(nèi)，竹叢密集，竹體高大，通常高度在3 m以上。在高山灌叢和高山草地中未見竹子分布。在高海拔段的稀疏冷云杉林中，香柏、紅杉、柳樹、高山櫟、杜鵑等物種的蓋度較大，冷箭竹和峨熱竹通常呈低矮的禾草狀，并且較為稀疏，當(dāng)?shù)厮追Q“麥秧子”，不注意很難識(shí)別。

竹子對(duì)溫度和光照的需求是有范圍限制的，在高海拔地段，通常一個(gè)因素就可以決定其生存。我們?cè)诰琵垳耪{(diào)查發(fā)現(xiàn)，高海拔段陽坡雖然光照充分，溫度也相對(duì)較高，在和有竹類分布的陰坡同一海拔高度的陽坡或者半陽坡上卻沒有竹類分布。模型中如果加入植被和坡向等數(shù)據(jù)可以進(jìn)一步縮小預(yù)測(cè)偏差。

物種的生存除了要滿足氣候條件外(基礎(chǔ)生態(tài)位)，還需要滿足土壤、植被、生物互作等因素，因此實(shí)際分布區(qū)小于基礎(chǔ)生態(tài)位[9]。而且生物氣候關(guān)系在外推時(shí)還有關(guān)系發(fā)生變化的可能，距離越遠(yuǎn)，這種變化的危險(xiǎn)性越大。收集高質(zhì)量的物種分布數(shù)據(jù)、采用更高精度的氣候數(shù)據(jù)、加入更多的環(huán)境變量對(duì)于準(zhǔn)確預(yù)測(cè)物種實(shí)際生態(tài)位會(huì)起到積極的推動(dòng)作用。

4.3 模型的適用性

本文采用集合模型來預(yù)測(cè)物種分布，包括10種模型，每種模型都將數(shù)據(jù)拆成兩份，其中一份用于訓(xùn)練，另一份用于驗(yàn)證，其結(jié)果可以通過敏感性和專一性的測(cè)量來進(jìn)行有效性判別。數(shù)據(jù)的劃分是采用隨機(jī)抽樣的方法生成的，對(duì)多次抽樣生成的數(shù)據(jù)集進(jìn)行模擬，然后求平均值，得到該種模型相應(yīng)的分布預(yù)測(cè)。相比于單用某一模型，集合模型通過模型投票決定最終的預(yù)測(cè)結(jié)果，這樣可以防止模型的過擬合以及彌補(bǔ)模型結(jié)構(gòu)差異帶來的缺陷，從而獲得一個(gè)較為穩(wěn)定的預(yù)測(cè)結(jié)果。竹類分布海拔存在一些異常值，與未檢測(cè)到相反，這些是確定存在的物種分布點(diǎn)，而未檢測(cè)到卻是實(shí)際存在，但是沒有發(fā)現(xiàn)的情況。因此可能存在海拔分布的連續(xù)性，但是由于沒有調(diào)查或者調(diào)查了卻沒有發(fā)現(xiàn)而造成海拔分布的間斷性。集合模型將這些異常值多做拋棄處理，可以看到檢測(cè)點(diǎn)處沒有預(yù)測(cè)分布。這種問題可能是取樣偏差造成的，均勻布設(shè)樣地開展取樣工作可以減少這種問題。

集合模型的基礎(chǔ)模型都是屬于相關(guān)模型，通過關(guān)聯(lián)找到植物與氣候之間的反應(yīng)曲線，進(jìn)而預(yù)測(cè)到全區(qū)域適宜生態(tài)位。我們的研究只是考慮了氣象因素，而這種數(shù)據(jù)集的相關(guān)性隨著區(qū)域發(fā)生改變，例如Li等[7]在整個(gè)熊貓分布區(qū)的分布研究中，發(fā)現(xiàn)有8個(gè)氣象因素滿足要求，而我們則只找到6個(gè)。說明前者環(huán)境異質(zhì)性更強(qiáng)，大區(qū)域中變量方差變大，導(dǎo)致相關(guān)性降低。但是局部區(qū)域變量之間的相關(guān)性依然很強(qiáng)，這會(huì)導(dǎo)致模型預(yù)測(cè)效果變差。因此將研究區(qū)域細(xì)分成相對(duì)均勻的小塊做模型分析，能夠更加準(zhǔn)確地預(yù)測(cè)物種的分布。