郭虹揚，史明昌，楊建英，陳春陽

（1. 北京林業(yè)大學(xué) 水土保持學(xué)院，北京 100083；2. 北京林業(yè)大學(xué) 教育部水土保持與荒漠化防治重點實驗室，北京 100083）

白洋淀及大清河水系是雄安新區(qū)水資源安全保障體系的重要部分，但氣候干旱、上游土壤流失嚴(yán)重、森林水源涵養(yǎng)功能低下，導(dǎo)致白洋淀干淀、上游斷流，難以為新區(qū)發(fā)展提供保障[1]。土地或巖石的裸露、造林樹種的不適宜以及森林空間結(jié)構(gòu)不良等是森林水源涵養(yǎng)功能低下的主要原因[2-3]。早期植樹造林時，主要根據(jù)區(qū)劃及經(jīng)驗進行判斷，導(dǎo)致油松Pinus tabulaeformis等部分防護林違背適宜性原則，致使其逐漸衰退為低效林，未達預(yù)期效果[4-5]。張春霞等[6]以海拔、坡度、坡向為主要立地因子對陜西黃陵油松進行立地質(zhì)量評價，賴文豪等[7]通過地形、土壤、氣象等立地因子分析了油松的適宜性情況?？梢姡延醒芯恐饕性谛×饔虻刃〕叨确秶鷥?nèi)，難以為當(dāng)?shù)貐^(qū)域防護林的營造提供支撐，且各立地因子精度較低，如僅用個別站點數(shù)據(jù)直接代表整個區(qū)域，或克里金插值(Kriging)的氣象因子站點數(shù)據(jù)較少，不夠精準(zhǔn)化，進而導(dǎo)致適宜性分析結(jié)果無法貼合生產(chǎn)實際。本研究基于植物生態(tài)學(xué)特性，考慮了地形、土壤地質(zhì)、氣象等因子，且將各因子落實到30 m分辨率，采用了空間限制性因子疊加和最大熵模型(MaxEnt)，對較大尺度下油松的精準(zhǔn)適宜性空間分布進行了研究，以期為白洋淀大清河流域造林規(guī)劃提供科學(xué)依據(jù)。

1 研究區(qū)概況與研究方法

1.1 研究區(qū)概況

大清河地處海河流域中部 (38°10′～40°10′N，113°39′～117°34′E)，西起太行山，東臨渤海灣，北臨永定河，南至子牙河。河系跨山西、河北、北京、天津4省市，總面積43 060 km2，其中山區(qū)18 659 km2，丘陵平原24 401 km2，分別占河系總面積的43.33%和56.67%。大清河流域地處溫帶半干旱大陸性季風(fēng)氣候區(qū)，多年平均氣溫為6.58～13.26 ℃，多年平均降水量為365～626 mm，年內(nèi)降水量極為不均，多年平均積溫為858～4 664 ℃，多年平均無霜期為73～217 d。研究區(qū)內(nèi)優(yōu)勢植物有華北落葉松Larix principisrupprechtii、油松、側(cè)柏Platycladus orientalis、楊樹Populusspp.、刺槐Robinia pseudoacacia、臭椿Ailanthus altissima、繡線菊Spiraea salicifolia、荊條Vitex negundovar.heterophylla、胡枝子Lespedeza bicolor等。

1.2 數(shù)據(jù)來源

本研究數(shù)據(jù)分為地形、土壤地質(zhì)、氣象等3類，其中：地形數(shù)據(jù)來源于地理空間數(shù)據(jù)云(http://www.gscloud.cn/)；土壤地質(zhì)數(shù)據(jù)來源于中國寒區(qū)旱區(qū)科學(xué)數(shù)據(jù)中心(http://westdc.westgis.ac.cn)和全國地質(zhì)資料館(http://ngac.org.cn/Map/List)；氣象數(shù)據(jù)為站點數(shù)據(jù)，來源于中國氣象數(shù)據(jù)網(wǎng)(http://data.cma.cn/)各數(shù)據(jù)集。站點數(shù)量為79個，其中流域內(nèi)36個，流域外43個，各地理環(huán)境變量如表1。

表1 地理環(huán)境變量數(shù)據(jù)Table 1 Geographical environment data

在白洋淀大清河流域上游山區(qū)，以油松為優(yōu)勢種，隨機選取不同高程、坡度、坡向、坡位、土壤類型、齡組的樣方進行調(diào)查，高程按照中山、低山區(qū)間結(jié)合太行山植物垂直地帶性進行考慮[8]，記錄樹高、胸徑、郁閉度、植被覆蓋度、林下灌草種類及蓋度和地理位置等信息。在樣方調(diào)查時，山脈、高大植物會遮擋全球定位系統(tǒng)(GPS)衛(wèi)星信息接收[9]，故地理位置信息通過手持GPS儀結(jié)合衛(wèi)星地圖比對位置并標(biāo)注后自動獲取，其精確度要高于手持GPS儀直接獲取的位置數(shù)據(jù)，樣方分布如圖1所示。根據(jù)植物學(xué)和植物地理學(xué)等原理，在篩選油松可能適宜性空間單元時，考慮人為干預(yù)較少的樣方和海拔100 m以上地區(qū)，因此，本研究調(diào)查樣方主要分布于流域西部山區(qū)，共計190個30 m×30 m油松樣方數(shù)據(jù)(圖1)。

圖1 白洋淀大清河流域油松樣方分布示意圖Figure 1 Distribution of P. tabulaeformis quadrats in Baiyangdian Daqing River Basin

1.3 地理環(huán)境數(shù)據(jù)處理方法

采用多元線性回歸克里金法(multiple linear regression Kriging, MLRK)對氣象站點數(shù)據(jù)進行插值，得到連續(xù)氣象空間分布數(shù)據(jù)，其主要公式如下[10-11]：

在進行氣象數(shù)據(jù)插值時，將高程、坡度、坡向數(shù)據(jù)標(biāo)準(zhǔn)化后作為輔助變量，利用逐步回歸法篩選后進行回歸擬合，計算多元線性回歸在氣象站點位置的殘差，再利用普通克里金法對回歸殘差進行插值，將每個預(yù)測點的確定部分和殘差插值結(jié)果部分相加得到各點的預(yù)測值[12]。本研究在插值時考慮到研究區(qū)面積較大且跨越不同氣候區(qū)，故將研究區(qū)按照氣候區(qū)劃分為2個區(qū)，即半濕潤大區(qū)與半干旱大區(qū)，分別擴寬30 km進行插值。為提升流域邊緣結(jié)果的準(zhǔn)確性，站點除選取流域內(nèi)的36個站點外，另選取流域周圍43個站點納入插值過程，插值完成后再進行拼接得到連續(xù)的氣象要素空間分布數(shù)據(jù)，分辨率為30 m，根據(jù)精度檢驗，各氣象數(shù)據(jù)MLRK插值結(jié)果交叉驗證值R2均大于0.9，插值結(jié)果達到準(zhǔn)確水平。

1.4 精準(zhǔn)適宜性研究方法

1.4.1 油松可能適宜空間單元篩選 根據(jù)油松的生態(tài)習(xí)性，利用限制性因子疊加[13]，得到油松生長的可能適宜單元。油松是中國特有樹種，主產(chǎn)華北、西北；生于海拔100～2 600 m，為研究區(qū)中、低山區(qū)重要的造林樹種。喜光，適干冷氣候，產(chǎn)區(qū)極端低溫-25 ℃，年降水量400～750 mm，耐干旱，也能生長在年降水量250 mm的部分山區(qū)；在酸性、中性、鈣質(zhì)土壤中均能生長，極耐干瘠，石質(zhì)荒山亦能生長，不耐水澇和鹽堿[14-15]?？紤]到以上生態(tài)習(xí)性，空間疊加模型為：

式(2)中：PSG表示油松可能適宜單元，當(dāng)PSG=1時，為油松可能適宜空間單元，當(dāng)PSG=0，為不適宜；當(dāng)高程≥100 m時，D取值為1，反之則取值為0；當(dāng)極端最低氣溫≥-25 ℃時，T取值為1，反之則取值為0；S表示土壤類型，S為草甸土、沼澤土及城區(qū)、島嶼、鹽土?xí)r，取值為0，其余土壤類型為1；R表示巖石類型，R為水域時，取值為0，其余類型取值為1。由于大清河流域降水量符合油松降水量限制條件，故在疊加模型中暫不予考慮。

1.4.2 油松適宜性計算方法 油松適宜性計算主要采用最大熵模型(MaxEnt)。PHILLIPS等[16]基于最大熵原理編寫了MaxEnt生態(tài)位模型軟件，被廣泛用來預(yù)測各物種的潛在適宜分布區(qū)[17-19]。本研究使用MaxEnt 3.4.1，將分辨率為30 m的24個地理環(huán)境變量作為模型的輸入?yún)?shù)，地形變量如圖2，土壤地質(zhì)變量如圖3，氣象變量如圖4。

圖2 地形變量示意圖Figure 2 Terrain variables

圖3 土壤地質(zhì)變量示意圖Figure 3 Soil geological variables

圖4 氣象變量示意圖Figure 4 Meteorological variables

2 結(jié)果與分析

2.1 最大熵模型精度分析

經(jīng)MaxEnt軟件計算，得到受試者工作曲線(ROC)圖，該圖反映模型靈敏度，可評價預(yù)測模型的準(zhǔn)確程度，ROC曲線下面積(AUC)越大，評價模型越準(zhǔn)確[20]。由圖5可見：基于MaxEnt模型和地理環(huán)境變量構(gòu)建的油松適宜性預(yù)測模型的訓(xùn)練數(shù)據(jù)集AUC為0.869，大于0.8，達到很準(zhǔn)確水平；檢驗數(shù)據(jù)集AUC為0.817，大于0.8，達到很準(zhǔn)確水平。綜合模型準(zhǔn)確性，說明油松適宜性預(yù)測模型有效性較高，可用于白洋淀大清河流域油松精準(zhǔn)適宜性空間分布研究。

圖5 ROC曲線圖Figure 5 ROC curve

刀切法(Jackknife)[21-22]可預(yù)測各地理環(huán)境變量對油松適宜性空間分布影響的得分情況，可確定各個環(huán)境變量的重要程度。從圖6可見：多年年平均積溫、多年年平均最大日降水量、高程、多年年平均氣溫標(biāo)準(zhǔn)差依次是對油松預(yù)測貢獻率最高的地理環(huán)境變量，為流域內(nèi)油松精準(zhǔn)適宜性預(yù)測提供了豐富信息。此外，當(dāng)忽略巖石類型變量時，會最大程度上降低地理環(huán)境變量總得分，因此，它具有最多其他變量中沒有的信息，也需要綜合考慮。

圖6 刀切法地理環(huán)境變量貢獻圖Figure 6 Jackknife environmental variable contribution graph

通過MaxEnt模型得到各地理環(huán)境變量與發(fā)生概率之間的反饋曲線，該曲線可反映單個地理環(huán)境變量值與發(fā)生概率之間的響應(yīng)關(guān)系。從圖7A可以看出：隨著多年年平均積溫的升高，分布概率呈先增大后減小的趨勢，以分布概率0.5為界，積溫適宜區(qū)間為2 800～3 800 ℃。由圖7B看出：隨著高程升高，分布概率先增大后減小，以分布概率0.5為界，高程適宜區(qū)間為500～1 200 m。由圖7C看出：隨著多年年平均最大日降水量升高，分布概率先增大后減小，以分布概率0.5為界，年平均最大日降水量適宜區(qū)間為70～110 mm。由圖7D看出：隨著多年年平均氣溫標(biāo)準(zhǔn)差升高，分布概率先增大后減小，以分布概率0.5為界，年平均氣溫標(biāo)準(zhǔn)差適宜區(qū)間為0.31～0.42 ℃。

圖7 主導(dǎo)地理環(huán)境變量反饋曲線Figure 7 Leading environmental variable feedback curve

2.2 適宜性結(jié)果分析

參考李培琳[23]、麻亞鴻[24]的適宜性分級標(biāo)準(zhǔn)，確定概率在0～0.05為不適宜區(qū)，0.05～0.30為低適宜區(qū)，0.30～0.50為中適宜區(qū)，0.50～1.00為高適宜區(qū)。白洋淀大清河流域油松高適宜區(qū)面積為1 924 km2，中適宜區(qū)為3 477 km2，低適宜區(qū)為7 621 km2，不適宜區(qū)為30 038 km2，中、高適宜區(qū)面積達到流域總面積的12.54%，中、高適宜區(qū)主要分布在大清河流域上游山區(qū)的西北部，在靈丘縣、淶源縣、淶水縣、易縣、房山區(qū)、阜平縣等區(qū)縣。油松精準(zhǔn)適宜性空間分布及其細節(jié)放大圖如圖8，分辨率為30 m。表2為油松高適宜區(qū)面積大于20 km2的縣，可見，高適宜區(qū)主要集中在淶源縣、靈丘縣、易縣、淶水縣、清苑縣、涿鹿縣，其中淶源縣面積最大，為462.60 km2，占流域總面積的1%，這是由于這些縣符合油松適宜條件的區(qū)域較多，在利用油松林進行綠化時可著重考慮。

圖8 油松精準(zhǔn)適宜性分布及細節(jié)放大圖Figure 8 Accurate suitability distribution and detailed enlarged drawing of P. tabulaeformis

表2 各縣高適宜區(qū)面積及比例Table 2 Area and proportion of high suitable area in each county

3 討論

油松是華北地區(qū)主要造林樹種，由于早期造林時，沒有考慮油松適宜性或僅根據(jù)氣候區(qū)劃種植了大面積的人工油松林，使得油松生長緩慢，或衰退現(xiàn)象明顯[25-26]。近年來，針對油松適宜性的研究僅停留在小尺度，且數(shù)據(jù)和研究成果精度不高，難以為油松防護林營造提供科學(xué)支撐。

本研究從油松生態(tài)習(xí)性的角度，得出影響油松適宜分布的主導(dǎo)地理環(huán)境變量有積溫、最大日降水量、高程、年平均氣溫標(biāo)準(zhǔn)差，貢獻率均大于20%，高適宜區(qū)分布在積溫2 800～3 800 ℃、高程500～1 200 m、年平均最大日降水量70～110 mm、年平均氣溫標(biāo)準(zhǔn)差0.31～0.42 ℃的區(qū)域內(nèi)。張春霞等[6]認為：油松中低海拔(900～1 100 m)立地質(zhì)量為中等以上水平，與本研究高程最大值接近，最小值有差異，是由于其調(diào)查樣方集中于900 m以上。魏安琪[27]認為：坡向是影響油松適宜性最大的因子，但在本研究中坡向不是主導(dǎo)因子，這是由于該研究僅考慮了地形因子未考慮氣象因子，且本研究分析后認為坡向通過影響區(qū)域局部水熱條件而對氣象因子產(chǎn)生影響，這在氣候因子中有所體現(xiàn)。

在油松林管護過程中，針對不適宜導(dǎo)致的油松低效林，可間伐及間種其他適宜植物，通過逐漸演替由新植物種代替油松成為優(yōu)勢種；其次，在油松防護林營造時，需根據(jù)上述條件判斷該地塊是否適宜油松生長。本研究由于研究區(qū)范圍較大，僅考慮了土壤類型對油松生長的影響，沒有考慮土壤理化性質(zhì)[28]、土壤厚度[21]等，在防護林建設(shè)中需綜合考慮。此外，在以地理環(huán)境為角度進行油松適宜性研究時，未考慮林業(yè)經(jīng)營管理、經(jīng)濟因素、其他人為因素等。在防護林建設(shè)中，整地方式、苗木品質(zhì)、森林采伐等經(jīng)營管理技術(shù)均對油松有較大影響；從經(jīng)濟角度看，造林成本、管護投入等因素也會造成一定影響[29-30]，也需綜合考慮。