野生动物的分布受到其栖息区域景观空间格局的显著影响^[1]。景观空间格局既是景观异质性本身的具体体现, 又是野生动物栖息地生境特征的有效反映^[1]。因此, 了解景观空间格局特征对野生动物空间利用的影响, 对物种的科学保护和有效管理有着重要的意义^[2], 已有的景观格局与野生动物分布关系的研究多借助于景观指数对景观格局进行量化, 然后采用传统的线性回归模型等方法定量分析景观格局对野生动物分布的影响^[3-9]。然而, 生态学关系的空间异质性或空间非均匀性(即变量间的关系或结构会随着地理位置的变化而改变), 是空间分析中一个普遍存在的现象^[10]。传统线性回归模型, 如普通最小二乘法(Ordinary Least Square, OLS)模型, 将自变量与因变量的关系视为全局不变, 得到的回归参数是整个研究区域内的平均值, 从而忽略了两者的空间变异性所导致的关系变化, 并因此造成回归模型有效性降低和预测结果失效等问题^{[10, 11]}。

地理加权回归(Geographically weighted regression, GWR)模型是传统回归分析方法的扩展, 是探索空间关系异质性方面有效工具之一^[10]。GWR模型是一种典型的局部空间回归模型, 通过将全局参数分解成局部参数进行估计, 对每一个空间位置点的参数进行估计时考虑了非平稳性的关系, 因此能深刻解释地理空间数据的某类指标和影响因子之间的空间非平稳性关系, 这是传统OLS模型无法比拟的^[12-13]。近年来, 许多学者对GWR模型进行了研究, 主要集中在应用GWR模型分析社会环境因子对区域经济、区域房价和区域污染等方面的研究^[14-18], 如关伟和郝金连^[14]运用OLS和GWR模型分析东北地区旅游经济和旅游产业因子、消费因子、投资因子之间的关系, 发现GWR模型拟合优度比OLS模型有显著提高。陈辉等^[17]采用GWR模型构建了我国区域范围内近地表的PM_2.5遥感反演模型, 结果表明GWR模型既能体现PM_2.5时空分布的全局变化特征, 又能体现局部空间异质性。陈强等^[18]基于GWR模型评估了土地利用对地表水质的影响, 其结果也证明GWR模型在预测精度和处理空间自相关过程中都要优于OLS模型。然而, 在野生动物栖息地选择与利用研究领域, GWR方法尚未得到应用。

大熊猫(Ailuropoda melanoleuca)是中国乃至世界野生动物保护的旗舰物种, 是中国特有珍稀濒危物种^{[19, 20]}。近年来应用景观生态学原理与方法对大熊猫的研究很多, 内容涉及生境适宜性评价、潜在栖息地分布模拟、景观连通性分析及栖息地景观破碎化等^[21-26]。这些研究主要是利用景观格局指数描述大熊猫栖息地特征, 或者使用传统线性回归方法分析景观格局对大熊猫的影响, 很少有研究涉及从整体到局部探讨环境因子对大熊猫空间分布的影响。本文结合前人的研究结果, 采用GWR方法探讨环境异质性对秦岭大熊猫空间分布的影响, 并通过与传统的OLS分析方法进行比较, 检验GWR模型在大熊猫栖息地空间分析中是否可用, 以期提高大熊猫-环境空间关系的科学理解, 为进一步开展大熊猫等野生动物的栖息地选择与利用精准量化分析和评价提供参考。

1 材料与方法 1.1 研究区域概况

本研究区域介于北纬33°7′—34°6′、东经106°21′—108°54′和海拔376—3770 m之间, 属于我国传统南北地理分界线的秦岭中段地区, 辖陕西省17个县部分区域, 且包含19个自然保护区, 面积约为2万km²(图 1)。该研究区属于季风气候, 是北亚热带和暖温带的过渡区, 是长江和黄河的分水岭, 同时包含秦岭山地、低山丘陵和平原等一系列地貌类型, 是我国生物多样性和生态系统类型最为丰富的地区之一。秦岭大熊猫主要分布在秦岭中段南坡, 北坡和西段有少量分布^[27]。总体上为两块相隔较远的区域, 其西端的宁强县和甘肃康县、四川青川连成一片；佛坪、太白、洋县、周至、宁陕、留坝、城固、户县和镇安9县的分布区相隔较近或连成一片, 分布区域横跨长江、黄河水系^[27]。

1.2 数据来源及处理

本研究所采用的秦岭大熊猫分布数据来源于1999—2002年全国第三次大熊猫调查数据。根据秦岭大熊猫平均家域面积为10.62 km^2[28], 为避免空间自相关影响模型预测, 剔除距离小于2 km的冗余痕迹点, 最终保留267个大熊猫分布痕迹点进入分析(图 1红点所示)。鉴于OLS与GWR建模要求使用连续型数值变量, 我们利用ArcGIS 10.2软件的核密度(Kernel Density)分析功能估计了每一个大熊猫分布点所在位置的大熊猫分布相对密度。

研究区的数字高程模型(DEM)数据源于美国SRTM(Shuttle Radar Topography Mission, SRTM)30 m精度数据；公路和居民点GIS数据源于中国科学院资源环境科学数据中心的国家基础矢量数据集数据；土地覆盖/利用数据源于中国土地利用现状遥感监测数据库2000年土地利用30 m精度数据(http://www.resdc.cn), 该数据集的土地利用类型包括耕地、林地、草地、水域、居民地和未利用土地等6个一级类型和包括有林地、灌木林、疏林地等25个二级类型, 数据分类精度均可达到85%以上^[29]。参考以往大熊猫生态学方面的研究^[21-26], 本研究选取了有林地、农田、灌草地等三类与大熊猫分布相关的土地覆盖类型进入下一步景观空间格局量化分析。

1.3 环境变量筛选

本研究选取了以下3个方面的环境因素来分析大熊猫空间分布与环境异质性之间的关系：1)地理特征(包括海拔和坡度)；2)人类影响(包括距离道路的距离和距离居民点的距离)；3)景观格局特征(包括最大斑块指数、边缘密度、斑块面积、连接度指数和聚合度指数)^[30]。每一个大熊猫分布点处的人类影响因素通过ArcGIS 10.2软件中的距离量算功能计算距离道路和居民点的最近距离。大熊猫分布点处的景观指数值是在Fragstats 4.0软件中选择User provided points模块, 设置搜索半径为3 km计算完成的^[30]。为避免景观指数值的空间自相关, 本文利用R 3.3.2软件对所提取景观格局指数值进行因子分析, 通过变量间的相关关系剔除信息严重冗余变量, 将特征值>0.8的变量保留, 特征值< 0.3的变量移除^[30], 最终选择农田斑块面积(Cropland-CA)、有林地斑块面积(Forest-CA)和灌草地斑块面积(Grassland-CA)等独立的景观变量作为模型变量参与下一步的模型运算。最终进入建模分析的环境变量见表 1。

1.4 模型建立与评价

本研究首先选用经典的OLS进行建模。通过OLS建模, 既可以初步了解变量间的相互关系, 也可以评估OLS对空间数据的建模效果, 进而确定是否选择GWR进一步建模。本文所有的模型分析在GWR 4.0软件中完成。

1) OLS建模。OLS是一种对因变量和解释变量之间相互关系进行研究的统计方法, 是所有空间回归分析的起点, 假设线性回归关系满足全局空间平稳条件, 其公式如下：

(1)

式中, Y_i是第i点因变量的值, β₀为截距, X_ik为第k个解释变量在第i点的值, β_k为第k个解释变量的斜率或回归系数, ε_i为残差。

2) GWR建模。GWR是Fotheringham等在传统OLS模型基础上将数据的地理位置加入回归参数中, 同时考虑了相邻点的空间权重, 允许局部参数估计的地学统计方法^{[31, 32]}。其公式如下：

(2)

式中, (μ_i, ν_i)是第i点的空间位置, β_k(μ_i, ν_i)是连续函数β_k(μ, ν)在(μ_i, ν_i)处的值。

由于大熊猫的分布在空间上呈集聚分布, 故GWR模型采用调整型的空间核和Gaussian权重函数, 以AICc作为衡量标准, 通过黄金分割搜索(Golden section search)方法来确定最佳带宽, 也就是说, 使得AICc最小的带宽即为最佳带宽。

3) 模型评价。为了比较OLS模型和GWR模型的拟合效果, 选择修正的AIC信息准则(Akaike, 1973)作为模型评价指标, 其公式如下：

(3)

式中, n表示观测值数量；σ表示误差项估计的标准差；tr(S)表示GWR模型S矩阵的迹, 是带宽的函数。一般情况下, 具有较小AIC_C值的模型拟合效果相对较好。若两个模型间的AIC_C差>3, 则表明两个模型之间具有明显差异, 具有较小AIC_C值的模型拟合度更优。此外, 还可以通过计算模型的残差标准差(Sigma)来比较OLS模型和GWR模型, Sigma值越小, 表示模型的拟合效果越好。

2 结果与分析 2.1 最小二乘法(OLS)模型分析

以大熊猫相对分布密度为因变量建立OLS回归方程结果表明(表 2、图 2), 大熊猫分布相对密度与森林面积占比和灌草地面积占比呈显著相关关系, 与海拔和距离道路距离也呈现出一定的相关性。各解释变量的方差膨胀因子(VIF)均小于7.5, 表明各变量间共线性对模型的影响不显著。OLS模型的R²为0.152, 表明模型的拟合度不理想。OLS模型残差的全局Moran′s I系数为0.16, Z-score为10.44, 说明残差项还存在空间自相关性。以上说明, OLS对数据空间信息的提取还不完全, 不能很好地解释大熊猫空间分布密度与环境影响因素的空间非均匀性。

2.2 地理加权回归(GWR)模型分析

使用与OLS模型相同的变量构建GWR模型, 依据最小AIC值确定最佳带宽为107获得各变量系数分布(表 3、图 3), 可以看到, 森林面积占比和灌草地面积占比等变量的参数估计值在OLS方法下是正值, 但在GWR方法下则有正有负, 说明各环境变量对大熊猫空间分布的影响呈现出空间不均匀性。

GWR模型与OLS模型相比(表 4), GWR模型的R²达到了0.468, 且与OLS模型的AICc之差远大于3(ΔAICc=24.401), 说明GWR模型的拟合度明显高于OLS模型；GWR模型残差的方差分析结果显示, GWR模型的残差平方和(Residuals SS)下降了84.088, 亦表明GWR模型的拟合效果比OLS模型有很大改善(Residuals F = 2.368)；GWR模型残差的Moran′s I系数为0.09, Z-score为6.08, 说明残差不存在空间自相关。综上, 结果表明GWR模型比OLS模型模拟的效果要好。

2.3 影响因子空间分异特征举例分析

GWR模型是局部模型, 在每个样本数据点都有一组局部的参数估计, 而不像OLS模型那样是一个全局或平均意义上的估计值。因此, 我们对GWR模型的局部参数估计值使用IDW方法进行空间插值, 得到各环境变量的局部参数估计的空间格局图(图 4), 以此展示各环境变量对大熊猫分布影响的空间分异特征。

如图 4所示, 距居民点距离对大熊猫相对分布密度的影响作用自西向东呈递增趋势, 形成以凤县、勉县、留坝县、城固县、太白县和洋县西北部为主的低值区和以周至县和佛坪县为主的高值区；距道路距离对大熊猫相对分布密度的影响作用从中心向四周呈圈层递增趋势, 形成以太白县东南部、周至县西南部、佛坪县西北部和洋县东北部为主的低值区和以凤县、城固县、留坝县、勉县和略阳县为主的高值区格局；而有林地面积占比对大熊猫分布相对密度的影响作用从中心向四周呈圈层递增趋势, 形成以周至县南部和佛坪县北部为主的低值区和以城固县为主的高值区格局, 表明有林地斑块面积占比对周至县南部和佛坪县北部的大熊猫核心分布区域的相对分布密度的影响作用不明显。

3 讨论

本研究分别应用OLS模型和GWR模型分析了环境异质性与大熊猫空间分布之间的关系, 并将两者的模拟结果进行比较, 证明GWR模型对于刻画大熊猫空间分布与景观格局之间的关系具有独特的空间非平稳性和尺度依存特性优势。OLS模型中多采用全局空间回归模型量化大熊猫分布与环境变量之间的关系, 得到的回归参数估计是在整个研究区域内的平均值, 不能反映物种-环境关系的空间异质特征；而GWR模型允许将全局参数的估计分解成局部参数进行估计, 在模拟精度和空间特征提取上都显著优于OLS模型, 模型的结果可反映环境变量的影响具有空间非均匀性, 这是传统OLS模型无法实现的。

参考已有大熊猫与环境异质性之间关系的研究^[21-26], 本文的GWR分析结果进一步说明环境异质性对大熊猫空间利用的影响呈现出明显的空间非平衡性, 不同区域大熊猫空间分布与主导因子发生作用的特征尺度也不同。例如, 有林地斑块面积对大熊猫相对分布密度的影响作用从中心向四周呈圈层递增趋势, 形成以周至县南部和佛坪县北部为主的低值区和以城固县为主的高值区格局, 表明有林地斑块面积对周至县南部和佛坪县北部的大熊猫核心分布区域的相对分布密度的影响作用不明显, 而对城固县大熊猫相对分布密度的影响作用要大于其他区域。距离道路距离对大熊猫相对分布密度的影响作用从中心向四周呈圈层递增趋势, 形成以太白县东南部、周至县西南部、佛坪县西北部和洋县东北部为主的低值区和以凤县、城固县、留坝县、勉县和略阳县为主的高值区格局, 说明道路对凤县、城固县、留坝县、勉县和略阳县大熊猫相对分布密度的影响作用要大于其他区域, 这就要求我们保护区管理规划者要格外注意对不同区域采用整体和局部相结合的办法, 做到因地制宜, 提升保护措施有效性。

虽然GWR模型在理解野生动物空间分布与环境异质性之间的空间关联具有很大优势, 但模型本身也存在一些不足之处。首先, 由于回归模型对环境变量数依赖性较大, 当环境变量数较少时, 难以建立合理的回归模型, 对因变量的预测不准确；而当环境变量数较多时, 又会产生多重共线性, 同样使得结果出现不确定性。其次, GWR模型不支持类型变量的加入, 如土地利用类型和保护地类型等, 而野生动物空间分布却受这些类型要素的影响较大, 忽视重要的类型变量会降低模型预测精度。如何将类别变量合理的应用到GWR模型中还需要进一步研究。此外, 本研究仅分析了大熊猫空间分布与部分环境因子之间的回归关系, 实际上, 除了这些环境因子之外诸多环境因子也影响着大熊猫的空间分布, 如气候、水域和食物等。下一步工作将把这些环境因子融合到GWR模型中, 以期得到更为全面客观的分析结果。

4 结论

本文基于地理加权回归模型(GWR)和普通最小二乘法线性回归模型(OLS)探究了环境异质性对秦岭大熊猫空间利用的影响, 并对模型模拟结果进行了比较。主要结论如下：

(1) GWR模型可解释环境异质性对大熊猫空间利用影响的46.8%, 优于OLS模型的15.2%。

(2) GWR模型适用于野生动物分布与环境变量之间关系的建模分析, 能有效的刻画野生动物分布与环境变量之间的空间异质性, 可为探究物种-环境关系的空间异质特征提供一种新的思路和方法。

(3) 大熊猫空间分布具有明显的空间异质性特征, 不同环境因子对大熊猫空间分布影响的作用不同, 现有大熊猫空间分布格局受多重因素综合影响。

致谢: 感谢陕西省林业局、陕西周至老县城国家自然保护区以及陕西佛坪国家自然保护区给予的技术支持。