近年来，因土地利用开发不合理和城市无序蔓延扩张，城郊区陆地生态系统结构和功能遭到严重破坏，致使区域水文过程、生物地球化学循环、生物多样性发生了巨大变化[1]，产生了工业“三废”污染、森林面积减少、土地景观破碎化等一系列生态环境问题。研究城郊区土地景观变化驱动力，可以为制定生态环境管理和生态风险防控措施提供科学依据，对缓解城郊区生态环境问题和促进区域土地持续利用、城镇合理扩张、经济良性发展具有重要意义。
目前，常用景观变化驱动力识别方法为传统统计分析法，如典型相关分析法[2]、Logistic回归模型分析法[3]、多元线性回归模型分析法[4]等。这类方法应用前提是数据必须符合统计上独立且均匀分布的假设条件[5]，而土地景观格局及驱动因子往往不独立，存在空间依赖关系，应用传统统计分析法进行驱动因子识别可能会产生偏差，所以有必要应用能体现数据空间相关性的空间回归模型进行土地景观格局变化驱动力分析。空间回归模型假设空间统计数据具有空间多维特征和时空相关性，能利用研究对象空间分布信息，较好地揭示土地景观格局变化影响因素及其时空分布。该模型不要求数据独立，可以充分分析数据空间属性，广泛应用于社会经济、农业环境、公共卫生等领域[6-7]以及空间分布预测模拟[8]、影响因素识别分析[9-10]等方面，但在土地景观格局变化驱动力分析中应用还比较少见。
龙泉驿属四川省会成都市辖区，是典型大城市近郊区，21世纪以来，该区在加快推进工业化和城镇化过程中，也面临诸如森林植被遭受破坏、耕地资源减少等生态环境问题，对其进行土地景观格局变化驱动力研究具有较强代表性和研究价值。有鉴于此，本研究以龙泉驿为研究区，从自然驱动力和人文驱动力2个方面构建土地景观格局变化驱动力指标体系，应用空间回归模型对—年、—年2个时期土地景观格局变化驱动力进行分析，以期能为区域生态安全格局构建和制定生态风险防控措施提供依据。
龙泉驿区位于成都平原东缘、龙泉山脉西侧，地处30°27′52″～30°43′23″N、°08′19″～°27′09″E，幅员面积.4km2。地质构造为成都断陷带与龙泉山隆褶带间的构造断块，地势由东南逐渐向西北微倾，最高海拔.3m，最低海拔.8m。有平坝、丘陵、山地等地貌类型，其中平坝.54km2、山地.38km2、丘陵21.48km2。属四川盆地中亚热带湿润气候，年平均气温16.5℃、降水量.4mm、日照時数h、蒸发量.7mm。成土母质为侏罗纪、白垩纪、第四纪成土母质，主要土壤类型为水稻土、黄泥土、紫色新冲积土和紫色土。龙泉驿属成都第二经济圈，是国务院批准的成都城市向东发展主体区，也是成都国家级经济技术开发区所在地，重点发展汽车产业，年全区GDP实现.6亿元，位居全省第一。
以研究区2期（年5月2日、年5月6日）TM影像和1期（年8月13日）OLI影像（来自USGS网站）为基础数据，年1∶[KG-*3]15万比例的土地利用现状图、年1∶[KG-*3]1万比例的土地利用现状图、年10—11月景观野外调查成果为辅助数据，应用QUEST遥感影像分类方法[11]解译得研究区土地景观类型图（图1）。自然驱动力指标中，—年、—年间年均降雨量和气温空间分布图根据龙泉驿气象局提供的降雨、气温数据，通过空间插值获取；高程分布图采用ASTERGDEMV2数据（来自NASA网站）；坡度图、坡向图均基于DEM计算得到；土壤有机质含量空间分布图根据龙泉驿区测土配方施肥项目测得的土壤有机质含量数据，通过空间插值获得。人文驱动力指标中，x7、x8、x9、x21、x23、x24、x25、x27、x28、x29、x35、x36数据来自《龙泉驿统计年鉴（—年）》和成都统计信息网，x10、x13、x14、x15、x17、x31、x32、x33、x34数据由龙泉驿统计局、农林局和各街道（镇、乡）提供，x11、x12、x18、x19、x20、x22、x26、x30数据根据已有数据计算得到。在此基础上，按照孙才志等采用的指标栅格化方法[3]，将各人文驱动力指标转化为栅格数据。各栅格数据分辨率统一为30m。
（1）建立空间回归分析基本单元网格。以0.6km为分析尺度，运用ArcGIS10.0构建景观格局变化驱动力分析单元网格矢量数据（图2-a）。研究区共涉及网格个。endprint（2）计算单元网格景观面积变化率。以—年农田景观面积变化率为例阐述景观面积变化率计算过程。首先，经重分类把农田景观变化图中变化区赋值为“1”、未变区赋值为“0”（图2-b）；然后，将单元网格矢量图和农田景观变化图叠加结果属性表导入MicrosoftExcel中统计各网格农田变化区、未变化区面积；最后，按式（1）计算各单元网格农田景观面积变化率。批量计算各指标在各单元网格中的均值。
（4）土地景观格局变化与驱动力指标的相关分析。应用SPSS19.0对各单元网格对应土地景观面积变化率、驱动指标均值进行Z-score标准化处理，计算土地景观面积变化率与驱动指标之间的皮尔森相关系数，删除相关系数未通过显著性检验和相关系数小于0.2的指标，筛选出与土地景观格局变化相关的指标。（5）应用普通最小二乘法（OLS）线性回归模型进行土地景观格局变化驱动力分析。应用OLS线性回归模型[9]分析驱动指标对土地景观格局变化的影响，以回归系数显著性水平、方差扩大因子VIF10为依据构建OLS线性回归模型，筛选出无多重共线性问题的指标。
（6）OLS线性回归模型残差、因变量、自变量空间相关分析与判断。应用OpenGeoDa计算1次Rook邻接、1次Queen邻接权重矩阵[16]下OLS线性回归模型残差、因变量、自变量MoransI值，根据P值和Z值判断MoransI显著性[17-18]。若模型残差、因变量、自变量存在空间自相关，则必须选择恰当空间回归模型进行驱动力分析。（7）空间回归模型选择。常见空间回归模型包括空间滞后模型（SLM）和空间误差模型（SEM），通用表达式为：式中：y为因变量；X为自变量；μ为模型残差；β为自变量空间回归系数；ε为白噪声；W1为反映因变量自身空间趋势的权重矩阵；
W2为反映模型残差空间趋势的权重矩阵；ρ为空间滞后项W1y系数，称为空间自回归系数；λ为空间误差系数。当ρ≠0，λ=0，时，为空间滞后模型，方程为y=ρWyβX+μ；当ρ=0，λ≠0时，为空间误差模型，方程为y=βX+λWμ+ε。根据拉格朗日乘数（LM-lag和LM-error）、稳健拉格朗日算子（RobustLM-lag和RobustLM-error）显著性确定空间回归模型类型[7-8]。（8）空间回归模型评估。采用确定系数（R2）、赤池信息准则（AIC）和施瓦茨准则（SC）和对数似然值（LG）模型评估指标[10]，对空间回归模型擬合结果进行评估，选取效果最好的进行驱动力分析。
各阶段OLS线性回归模型在1次Rook邻接、1次Queen邻接权重矩阵下的残差MoransI值相差不大（图3）。—年OLS线性回归模型在1次Rook邻接、1次Queen邻接权重矩阵下的残差MoransI值分别介于0.～0.、0.～0.之间，—年分别介于0.～、0.～0.之间，且各阶段OLS线性回归模型残差MoransI值显著性检验指标P0.01、Z1.96，其残差存在显著空间自相关性，表明OLS线性回归模型对土地景观格局变化驱动力的分析结果可靠性差。此外，各阶段OLS线性回归模型在1次Rook邻接、1次Queen邻接权重矩阵下的因变量、自变量MoransI值都较大，且彼此差异较小（图3）。基于1次Rook邻接（1次Queen邻接）权重矩阵计算得到的—年景观面积变化率及对应驱动指标MoransI指数最小值为0.（0.）、最大值为0.（0.），—年最小值为0.（0.），最大值为0.（0.4），且各MoransI指数显著性检验指标P0.01、Z1.96。表明—年各景观面积变化率及对应驱动指标存在显著空间自相关性，必须选择恰当空间回归模型进行土地景观格局变化驱动力分析。
空间回归模型评估指标R2、LG值均大于OLSM，AIC、SC值均小于OLSM（表2），表明空间回归模型拟合效果总体优于OLSM，因为空间回归模型残差MoransI值接近为零且均小于OLSM，基本排除残差空间自相关性，且考虑了变量空间自相关性[19]。农田景观在2个阶段1次Queen邻接权重矩阵对应空间滞后模型SLMQueen的R2和LG值最大、AIC和SC值最小，所以把SLMQueen拟合结果作为农田景观格局变化驱动力最终结果（表3）；果园景观在—年1次Rook邻接权重矩阵对应空间滞后模型SLMRook的R2和LG值最大、AIC和SC值最小，在—年1次Queen邻接权重矩阵对应空间误差模型SEMQueen的R2和LG值最大、AIC和SC值最小，所以分别把SLMRook和SEMQueen拟合结果作为—年、—年果园景观格局变化驱动力最终估计结果（表3）；
森林景观在—年1次Rook邻接权重矩阵对应空间误差模型SEMRook的R2和LG值最大、AIC和SC值最小，在—年SLMQueen的R2和LG值最大、AIC和SC值最小，所以分别把SEMRook和SLMQueen拟合结果作为—年、—年森林景观格局变化驱动力最终估计结果（表3）；城乡人居及工矿、交通运输、水体景观在2个阶段SLMRook的R2和LG值皆最大、AIC和SC值皆最小，所以把SLMRook拟合结果作为2个阶段3种景观格局变化驱动力最终估计结果（表3）。endprint表中，P、Z分别为MoransI值显著性检验统计量P值、Z值；SLM、SEM和OLSM分别为空间滞后模型、空间误差模型和普通最小二乘法线性回归模型；Queen、Rook分别表示1次Queen邻接、1次Rook邻接权重矩阵。
（1）农田景观格局变化驱动力分析。第一阶段（—），影响指标为粮食总产量、农业人口密度、农村用电、糧食单产、年末大牲畜存栏数、坡度和人口自然增长率；第二阶段（—），影响指标为农业人口、农业人口密度、坡度和年末大牲畜存栏数（表3）。2个阶段中，除农业人口密度、坡度、年末大牲畜存栏数3个指标相同外，其余指标各不相同，影响程度排前3位的是农业人口、农业人口密度和粮食总产量，表明人口状况是—年农田景观格局变化主要驱动力，这是由于人口数量、密度的变化可能导致粮食需求增大而扩大耕种面积、调整农业种植结构或导致城镇扩展而侵占大量农地，从而加剧了农田景观格局变化。
（2）果园景观格局变化驱动力分析。第一阶段，影响指标为人口密度、农业机械总动力、地方财政收入、耕地有效灌溉面积（表3）；第二阶段，第一产业产值占GDP比重、第三产业产值占GDP比重回归系数均未通过5%水平显著性检验（表3），表明产业结构变化在—年间对果园景观变化影响不大，这是因为这期间研究区已从农业区转变为以发展汽车产业为主的工业区，农业产业结构调整基本完成，加之果园景观主要分布在山区，受坝区汽车产业发展用地扩张影响较小，所以产业结构变化对果园景观影响不大；其余自变量回归系数均通过5%水平显著性检验（表3），是该阶段果园景观格局变化影响指标。2个阶段仅人口密度1个共有驱动指[CM（25]标，影响程度排前几位的是人口密度、农业人口密度、农业机械总动力、农村用电、地方财政收入和城镇居民人均可支配收入，表明人口状况、科技水平和经济发展是—年期间果园景观格局变化主要驱动力。
（3）森林景观格局变化驱动力分析。第一阶段，影响指标为坡度、土壤有机质含量、高程、非农业人口和化肥使用量；第二阶段，影响指标为粮食总产量、坡度、坡向、粮食播种面积、高程、年均降水量和农业人口密度（表3）。2个阶段只有高程、坡度2个共有驱动指标，影响程度排前几位的是高程、坡度、坡向、土壤有机质含量和粮食总产量，表明地形、土壤等自然驱动因子是—年间森林景观格局变化主要驱动力，由于在海拔较低、坡度较缓地带，土壤肥力较高、水热条件较好，土地开发利用较容易，处在此地带的森林景观最先受到人为影响而发生改变，但随着时间推移，易开发区域开发殆尽，一些高程较低、坡度相对较大区域的森林景观也受到人类干扰而发生变化。
（4）城乡人居及工矿景观格局变化驱动力分析。第一阶段，影响指标为坡度、第一产业产值占GDP比重、粮食播种面积，而且均为正效应影响（表3），表明坡度、第一产业产值占GDP比重、粮食播种面积较大区域，城乡人居及工矿景观变化越大，因为—年间，此类指标较大区域主要是农村，而这期间正是该区农村农民住房大量建设期，农村建设用地空间变化幅度较大，所以呈现出该部分区域城乡人居及工矿景观格局变化明显的特点。第二阶段，影响指标为坡度和年均降水量，其中坡度与之呈正相关（表3），表明坡度越大区域城乡人居及工矿景观变化越大，因为该区农村地形坡度普遍大于城市地区，—年，伴随着城镇化、工业化进程的加快，农村人口不断向城镇聚集，城镇、工业园区不断向农村蔓延扩张，在这双重因素驱动下使得研究区坡度较大区域城乡人居及工矿景观呈现出变化越大的特点。
区生产总值、高程2个自变量回归系数通过5%水平显著性检验，是交通运输景观格局变化主要影响指标（表3）。其中，地区生产总值回归系数为0.，远大于高程回归系数0.，表明经济发展水平是—年期间交通运输景观格局变化主要驱动力，这是因为经济越发达区域，城市化、工业化水平相对越高，地方财政势力相对越强，交通路网建设投入也较大，因而交通运输景观变化相对较大。
（6）水体景观格局变化驱动力分析。第一阶段，影响指标为地方财政收入、人口自然增长率和综合城镇化率；第二阶段，影响指标为农业人口密度和土壤有机质含量（表3）。2个阶段无一共同驱动指标，排在首位的指标分别是地方财政收入和农业人口密度，表明经济发展、人口状况等人文驱动因子是—年期间水体景观格局变化主要驱动力，这是由于这期间区域人口快速增长（—年总人口增幅达30.36%）、经济加速发展（—年GDP年均增长1.32倍）的压力对水体景观变化影响越发深刻，聚集经济效益使得经济水平越高区域，产业和人口聚集越大，土地需求量亦越大，势必导致这些区域部分水体转变为耕地或建设用地，呈现出变化较其他区域大的特点。
（1）OLS线性回归模型残差、自变量、因变量在1次Rook邻接、1次Queen邻接权重矩阵下均存在显著空间自相关性，说明研究区土地景观格局变化不仅与相关驱动因子有关而且还与邻近区域土地景观格局变化相关，忽略空间相关性的OLS线性回归分析结果存在偏差。空间回归模型在2种权重矩阵下的残差MoransI值均接近零，基本排除残差空间自相关性影响，拟合效果总体上优于OLS线性回归模型，这是由于空间回归模型引入空间权重矩阵，能充分挖掘数据空间特性，使分析效果更好[20]。
（2）各阶段各土地景观格局变化驱动指标差异较大，有部分土地景观在各阶段有少量共同驱动指标，但影响程度均不尽相同，这表明土地景观格局变化驱动因子受时间尺度影响较大，同一土地景观格局变化驱动因子会随时间推移而发生不同程度的变化，同一驱动因子对土地景观格局变化的影响力也会随时间变化而发生改变。因此，在进行景观格局变化驱动力分析时不能简单地将短期影响因子作为长期变化驱动力。endprint
（3）—年，人文驱动因子是研究区农田、果园、交通运输、水体景观格局变化主要驱动力，其中农田景观格局主要受人口状况影响，果园景观格局主要受人口状况、科技水平、经济发展影响，交通运输景观格局主要受经济发展影响，水体景观格局主要受经济发展、人口状况影响；自然驱动因子则是研究区森林、城乡人居及工矿景观格局变化主要驱动力，其中森林景观格局主要受地形、土壤等驱动因子影响，城乡人居及工矿景观格局主要受地形驱动因子影响。人文驱动因子对研究区景观格局变化的影响程度总体上大于自然驱动因子。
（4）鉴于部分数据缺失和部分指标较难量化，本研究在驱动指标选择时未能充分考虑水文、自然干扰、政策和文化因子，土地景观格局变化驱动指标体系还有待完善。此外，有研究证实局部空间回归模型（地理加权回归模型）在疾病空间数据影响因素筛选中比全局空间回归模型更可靠[19]，本研究限于篇幅仅对空间滞后模型、空间误差模型2种常用全局空间回归模型在土地景观格局变化驱动力分析中的应用进行了探讨，下一步可尝试将地理加权回归模型应用于景观格局变化驱动力分析中。