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

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

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

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

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

本文應用四川省第1高山——貢嘎山來研究竹類資源分布與生態模型尺度的關系。采用貢嘎山竹類調查分布位點和WorldClim生物氣候數據，利用集合模型構建了貢嘎山區3種竹類的分布預測，以期為保護區分區規劃和生物多樣性保護工作的依據。

1 研究區域

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

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

2 材料與方法

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

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

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

3 結果

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

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

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

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

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

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

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

4 討論

4.1 氣候異質性

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

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

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

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

4.2 物種特征

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

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

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

4.3 模型的適用性

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

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