Precise spatial distribution of suitability of Pinus tabulaeformis in Daqing River Basin, Baiyangdian

大清河地处海河流域中部(38°10′~40°10′N，113°39′~117°34′E)，西起太行山，东临渤海湾，北临永定河，南至子牙河。河系跨山西、河北、北京、天津4省市，总面积43 060 km²，其中山区18 659 km²，丘陵平原24 401 km²，分别占河系总面积的43.33%和56.67%。大清河流域地处温带半干旱大陆性季风气候区，多年平均气温为6.58~13.26 ℃，多年平均降水量为365~626 mm，年内降水量极为不均，多年平均积温为858~4 664 ℃，多年平均无霜期为73~217 d。研究区内优势植物有华北落叶松Larix principis-rupprechtii、油松、侧柏Platycladus orientalis、杨树Populus spp.、刺槐Robinia pseudoacacia、臭椿Ailanthus altissima、绣线菊Spiraea salicifolia、荆条Vitex negundo var. heterophylla、胡枝子Lespedeza bicolor等。

本研究数据分为地形、土壤地质、气象等3类，其中：地形数据来源于地理空间数据云(http://www.gscloud.cn/)；土壤地质数据来源于中国寒区旱区科学数据中心(http://westdc.westgis.ac.cn)和全国地质资料馆(http://ngac.org.cn/Map/List)；气象数据为站点数据，来源于中国气象数据网(http://data.cma.cn/)各数据集。站点数量为79个，其中流域内36个，流域外43个，各地理环境变量如表1。

在白洋淀大清河流域上游山区，以油松为优势种，随机选取不同高程、坡度、坡向、坡位、土壤类型、龄组的样方进行调查，高程按照中山、低山区间结合太行山植物垂直地带性进行考虑^[8]，记录树高、胸径、郁闭度、植被覆盖度、林下灌草种类及盖度和地理位置等信息。在样方调查时，山脉、高大植物会遮挡全球定位系统(GPS)卫星信息接收^[9]，故地理位置信息通过手持GPS仪结合卫星地图比对位置并标注后自动获取，其精确度要高于手持GPS仪直接获取的位置数据，样方分布如图1所示。根据植物学和植物地理学等原理，在筛选油松可能适宜性空间单元时，考虑人为干预较少的样方和海拔100 m以上地区，因此，本研究调查样方主要分布于流域西部山区，共计190个30 m×30 m油松样方数据(图1)。

Figure 1.  Distribution of P. tabulaeformis quadrats in Baiyangdian Daqing River Basin

采用多元线性回归克里金法(multiple linear regression Kriging, MLRK)对气象站点数据进行插值，得到连续气象空间分布数据，其主要公式如下^[10-11]：

式(1)中：$\hat z\left( {{s_0}} \right)$表示预测位置$ {{{{s}}_{\bf{0}}}}$点的插值结果，$ \displaystyle\sum \limits_{{{k}} = {\bf{0}}}^{{\rho }} {\hat {{\beta }}_{{k}}}{{{q}}_{{k}}}\left( {{{{s}}_{\bf{0}}}} \right) $是通过回归拟合的确定性部分；$ \displaystyle\sum \limits_{{{i}} = 1}^{{n}} {{{\lambda }}_{{i}}}{{e}}\left( {{{{s}}_{{i}}}} \right)$是普通克里金(OK)对回归残差的插值结果部分；$ {{k}}$表示回归拟合时的位置序号，$ \; \rho$表示空间位置总数，$ \; {\hat {{\beta }}_{{k}}}$是回归模型的系数，当k=0时，$ \; {\hat{{\beta }}_{\bf{0}}}$是截距；$i $表示回归残差插值时的位置序号，$ n$表示空间位置总数，${{{q}}_{{k}}}\left( {{{{s}}_{\bf{0}}}} \right) $是预测位置$ {{{{s}}_{\bf{0}}}}$点辅助变量的值，$ {{{\lambda }}_{{i}}}$是由回归残差的空间相关性结构确定的普通克里金插值的权重；$ {{e}}\left( {{{{s}}_{{i}}}} \right)$是位置$ {{{{s}}_{{i}}}}$处的残差。

在进行气象数据插值时，将高程、坡度、坡向数据标准化后作为辅助变量，利用逐步回归法筛选后进行回归拟合，计算多元线性回归在气象站点位置的残差，再利用普通克里金法对回归残差进行插值，将每个预测点的确定部分和残差插值结果部分相加得到各点的预测值^[12]。本研究在插值时考虑到研究区面积较大且跨越不同气候区，故将研究区按照气候区划分为2个区，即半湿润大区与半干旱大区，分别扩宽30 km进行插值。为提升流域边缘结果的准确性，站点除选取流域内的36个站点外，另选取流域周围43个站点纳入插值过程，插值完成后再进行拼接得到连续的气象要素空间分布数据，分辨率为30 m，根据精度检验，各气象数据MLRK插值结果交叉验证值R²均大于0.9，插值结果达到准确水平。

根据油松的生态习性，利用限制性因子叠加^[13]，得到油松生长的可能适宜单元。油松是中国特有树种，主产华北、西北；生于海拔100~2 600 m，为研究区中、低山区重要的造林树种。喜光，适干冷气候，产区极端低温−25 ℃，年降水量400~750 mm，耐干旱，也能生长在年降水量250 mm的部分山区；在酸性、中性、钙质土壤中均能生长，极耐干瘠，石质荒山亦能生长，不耐水涝和盐碱^[14-15]。考虑到以上生态习性，空间叠加模型为：

式(2)中：$ {P_{{\rm{SG}}}}$表示油松可能适宜单元，当$ {P_{{\rm{SG}}}}=1$时，为油松可能适宜空间单元，当$ {P_{{\rm{SG}}}}=0$，为不适宜；当高程≥100 m时，D取值为1，反之则取值为0；当极端最低气温≥−25 ℃时，T取值为1，反之则取值为0；S表示土壤类型，S为草甸土、沼泽土及城区、岛屿、盐土时，取值为0，其余土壤类型为1；R表示岩石类型，R为水域时，取值为0，其余类型取值为1。由于大清河流域降水量符合油松降水量限制条件，故在叠加模型中暂不予考虑。

油松适宜性计算主要采用最大熵模型(MaxEnt)。PHILLIPS等^[16]基于最大熵原理编写了MaxEnt生态位模型软件，被广泛用来预测各物种的潜在适宜分布区^[17-19]。本研究使用MaxEnt 3.4.1，将分辨率为30 m的24个地理环境变量作为模型的输入参数，地形变量如图2，土壤地质变量如图3，气象变量如图4。

Figure 2.  Terrain variables

Figure 3.  Soil geological variables

Figure 4.  Meteorological variables

经MaxEnt软件计算，得到受试者工作曲线(ROC)图，该图反映模型灵敏度，可评价预测模型的准确程度，ROC曲线下面积(AUC)越大，评价模型越准确^[20]。由图5可见：基于MaxEnt模型和地理环境变量构建的油松适宜性预测模型的训练数据集AUC为0.869，大于0.8，达到很准确水平；检验数据集AUC为0.817，大于0.8，达到很准确水平。综合模型准确性，说明油松适宜性预测模型有效性较高，可用于白洋淀大清河流域油松精准适宜性空间分布研究。

Figure 5.  ROC curve

刀切法(Jackknife)^[21-22]可预测各地理环境变量对油松适宜性空间分布影响的得分情况，可确定各个环境变量的重要程度。从图6可见：多年年平均积温、多年年平均最大日降水量、高程、多年年平均气温标准差依次是对油松预测贡献率最高的地理环境变量，为流域内油松精准适宜性预测提供了丰富信息。此外，当忽略岩石类型变量时，会最大程度上降低地理环境变量总得分，因此，它具有最多其他变量中没有的信息，也需要综合考虑。

Figure 6.  Jackknife environmental variable contribution graph

通过MaxEnt模型得到各地理环境变量与发生概率之间的反馈曲线，该曲线可反映单个地理环境变量值与发生概率之间的响应关系。从图7A可以看出：随着多年年平均积温的升高，分布概率呈先增大后减小的趋势，以分布概率0.5为界，积温适宜区间为2 800~3 800 ℃。由图7B看出：随着高程升高，分布概率先增大后减小，以分布概率0.5为界，高程适宜区间为500~1 200 m。由图7C看出：随着多年年平均最大日降水量升高，分布概率先增大后减小，以分布概率0.5为界，年平均最大日降水量适宜区间为70~110 mm。由图7D看出：随着多年年平均气温标准差升高，分布概率先增大后减小，以分布概率0.5为界，年平均气温标准差适宜区间为0.31~0.42 ℃。

Figure 7.  Leading environmental variable feedback curve

参考李培琳^[23]、麻亚鸿^[24]的适宜性分级标准，确定概率在0~0.05为不适宜区，0.05~0.30为低适宜区，0.30~0.50为中适宜区，0.50~1.00为高适宜区。白洋淀大清河流域油松高适宜区面积为1 924 km²，中适宜区为3 477 km²，低适宜区为7 621 km²，不适宜区为30 038 km²，中、高适宜区面积达到流域总面积的12.54%，中、高适宜区主要分布在大清河流域上游山区的西北部，在灵丘县、涞源县、涞水县、易县、房山区、阜平县等区县。油松精准适宜性空间分布及其细节放大图如图8，分辨率为30 m。表2为油松高适宜区面积大于20 km²的县，可见，高适宜区主要集中在涞源县、灵丘县、易县、涞水县、清苑县、涿鹿县，其中涞源县面积最大，为462.60 km²，占流域总面积的1%，这是由于这些县符合油松适宜条件的区域较多，在利用油松林进行绿化时可着重考虑。

Figure 8.  Accurate suitability distribution and detailed enlarged drawing of P. tabulaeformis

油松是华北地区主要造林树种，由于早期造林时，没有考虑油松适宜性或仅根据气候区划种植了大面积的人工油松林，使得油松生长缓慢，或衰退现象明显^[25-26]。近年来，针对油松适宜性的研究仅停留在小尺度，且数据和研究成果精度不高，难以为油松防护林营造提供科学支撑。

本研究从油松生态习性的角度，得出影响油松适宜分布的主导地理环境变量有积温、最大日降水量、高程、年平均气温标准差，贡献率均大于20%，高适宜区分布在积温2 800~3 800 ℃、高程500~1 200 m、年平均最大日降水量70~110 mm、年平均气温标准差0.31~0.42 ℃的区域内。张春霞等^[6]认为：油松中低海拔(900~1 100 m)立地质量为中等以上水平，与本研究高程最大值接近，最小值有差异，是由于其调查样方集中于900 m以上。魏安琪^[27]认为：坡向是影响油松适宜性最大的因子，但在本研究中坡向不是主导因子，这是由于该研究仅考虑了地形因子未考虑气象因子，且本研究分析后认为坡向通过影响区域局部水热条件而对气象因子产生影响，这在气候因子中有所体现。

在油松林管护过程中，针对不适宜导致的油松低效林，可间伐及间种其他适宜植物，通过逐渐演替由新植物种代替油松成为优势种；其次，在油松防护林营造时，需根据上述条件判断该地块是否适宜油松生长。本研究由于研究区范围较大，仅考虑了土壤类型对油松生长的影响，没有考虑土壤理化性质^[28]、土壤厚度^[21]等，在防护林建设中需综合考虑。此外，在以地理环境为角度进行油松适宜性研究时，未考虑林业经营管理、经济因素、其他人为因素等。在防护林建设中，整地方式、苗木品质、森林采伐等经营管理技术均对油松有较大影响；从经济角度看，造林成本、管护投入等因素也会造成一定影响^[29-30]，也需综合考虑。