貢嘎山竹類資源分布模型模擬研究

李志明，王 乾,余國寶，陳緒文*

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

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

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

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

1 研究區(qū)域

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

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

2 材料與方法

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

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

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

3 結(jié)果

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

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

圖1 使用集體模型構(gòu)建的貢嘎山區(qū)3種竹類的空間分布預測圖(A為豐實箭竹，B為冷箭竹，C為峨熱竹)，白色方塊構(gòu)成的區(qū)域為預測分布區(qū)，紅色叉為實際調(diào)查分布位點，灰白色線包圍區(qū)域為貢嘎山自然保護區(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.

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

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

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

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

4 討論

4.1 氣候異質(zhì)性

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

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

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

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

4.2 物種特征

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

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

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

4.3 模型的適用性

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

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