海岸带是复杂、动态的地球表层自然系统^[1], 同时也是人类经济社会活动密集分布区。相比于内陆地区, 海岸带地区的海洋性特征和高强度的人类活动使得土地资源的利用出现时空双维的高强度复杂变化, 海岸带地区城市空间不断扩张^[2-4]、围填海开发规模不断增大^[5]。而这种变化也必将反作用于人类经济社会活动和海岸带生态环境系统, 在促进沿海地区经济社会发展的同时, 也造成滨海湿地萎缩、海岸带环境质量降低^[6-8]和生态系统服务功能下降^[9-11]。因此, 开展海岸带土地利用变化研究对海岸带生态环境保护^{[7, 12-13]}、海岸带资源统筹规划与管理^[14-17]具有重要意义。

冀津海岸带地区位于渤海西岸, 是国内主要的平原型潮滩分布区^[18], 同时也是目前国内围填海连片开发地区、港口集中分布区及滨海盐场集中分布区。近40年来, 尤其是2000年以来, 冀津沿海地区经济社会取得显著发展, 土地资源及海洋资源开发强度不断提高, 土地利用特征也出现相应的变化。现有土地利用变化研究^[19-23]多基于中国科学院土地资源分类系统^[24], 以行政区为单元进行分析, 对于冀津沿海地区的特殊用地类型分析不足, 对于土地利用时空变化特征的刻画有待深入, 研究时间段相对较短。为解决上述问题, 本文在现有土地资源分类系统的基础上, 参考相关文献研究及冀津地区土地资源利用现状特征, 对海岸带土地利用分类系统进行相应调整。在此基础上, 充分利用GIS工具, 对改革开放以来冀津地区海岸带土地利用时空变化及强度进行分析, 以期为冀津地区海岸带国土资源陆海统筹管理及生态环境保护提供数据支撑与决策参考。

1 研究区域

本文研究区域为河北省、天津市两地海岸带地区(图 1)。为综合考虑统计单元完整性和划分科学性, 陆域部分确定为冀津两地沿海县市区, 海域部分根据国际湿地公约^[25]延伸至水深6 m的区域, 总面积达17 308.64 km², 共涉及12个县市区。研究区域2018年GDP总值达11 156.52亿元^[26-27], 人口共计760.24万人, 分别占冀津两地总量的24.33%、8.34%。研究区范围内主要入海水系为滦河和海河水系; 气候为温带海洋性季风气候, 多年平均降水量介于500 mm~650 mm; 陆地地貌以滨海平原为主, 海拔多在10 m以下, 仅在研究区域北侧有丘陵和山地分布; 近海滩涂广布, 曹妃甸及秦皇岛市区附近海域深水岸线资源优良。

2 数据与方法 2.1 数据来源与处理

本文所使用的主要数据为中国科学院土地利用现状遥感监测数据库(www.resdc.cn)所提供的1980年、1990年、2000年、2010年、2018年共5期土地利用矢量数据。该数据通过人机交互式目视判读的方式构建^[28], 土地利用分类系统为中国科学院土地资源分类系统^[24], 主要包含耕地、林地、草地、内陆水体、建设用地、未利用地等6类, 该系统考虑了遥感解译特征及内陆地区普遍存在的用地类型。

为更好反映研究地区滩涂宽阔、围填海及城镇建设强度大的状况, 本文在分类系统的基础上, 参考Feng等^[29]、邸向红等^[30]的研究成果并考虑实际开发利用特征, 引入养殖盐田用地、待利用地、滨海滩涂、浅海海域等4类。养殖与盐田用地主要指用于渔业养殖及制盐等沿海低强度涉水开发活动的用地类型, 原多划分为内陆水体或建设用地, 是研究区域浅海及水域重要的利用类型。待利用地是指处于由其他用地类型向建设用地的过渡阶段但尚闲置待建设的土地类型, 原多划分为未利用地, 对于识别土地开发效率和城镇扩张趋势有重要意义, 因此本文将其单独识别并进行分析。滨海滩涂指理论最低潮面到海岸线之间的土地类型, 原属于未利用地, 是研究区域分布广阔的特色土地类型。浅海海域指理论最低潮面以下及水深6 m以上的海域。养殖与盐田用地、待利用地的解译标志见表 1。基于改进的土地利用分类体系, 采用对研究时段Landsat影像的目视解译方法, 对中国科学院土地利用现状遥感监测数据进行核实, 划分出养殖与盐田用地及待利用地, 并进一步补充滨海滩涂及浅海海域图斑, 形成冀津海岸带地区土地利用数据(图 2)。结合WorldView等更高精度影像及野外经验判别, 遥感核实结果符合土地利用实际开发状态。

2.2 土地利用变化及其驱动力分析方法

本文采用土地利用数量动态变化分析方法、土地利用强度分析方法及土地利用类型转换分析方法来反映改革开放以来冀津海岸带土地利用时空变化特征, 并使用地理加权回归的方法对土地利用分布及其变化驱动因素进行分析。

土地利用数量动态变化分析方法: 采用土地利用年变化率及土地利用动态度指数。土地利用年变化率(K_i)^[31]反映了指定土地利用类型i的年变化速率, 土地利用动态度S^{[28, 32]}可反映研究区域内土地利用动态变化的总体情况。土地利用变化的计算方法为:

${K_i} = \left( {\sum\nolimits_j^n {\frac{{\Delta {S_{i, j}}}}{{{S_i}}}} } \right) \times \frac{1}{{{t_1} - {t_0}}} \times 100\% , $

(1)

$S = \left( {\frac{{\sum\nolimits_{ij}^n {\Delta {S_{i - j}}} }}{{{S_a}}}} \right) \times \frac{1}{{{t_1} - {t_0}}} \times 100\% , $

(2)

其中$ {{S}}_{{i}} $S_i是研究时段开始时间第i类土地利用类型总面积, $ {∆{S}}_{{i}, {j}} $ΔS_{i, j}是研究时段内第i类土地利用类型与其他土地利用类型j相互转换后的净面积变化, S_i-j是研究时段内由i类型转变为其他土地类型的面积的总和, $ {{S}}_{{a}} $S_a为研究区土地利用总面积, t₁—t₀为研究时段时长。

土地利用强度指数变化反映了土地由低级利用类型向高级利用类型的演变^[33], 具体计算方法可见式(3)。本文在此基础上, 为反映土地利用强度时间变化特征, 引入土地利用强度变化量, 见式(4)。

$ L=100 \times \sum\limits_{i=1}^{n} A_{i} \times C_{i}, $

(3)

$ L_{s}=L_{t_{1}}-L_{t_{0}}, $

(4)

其中, L是土地利用程度综合指数, A_i是第i级的土地利用程度分级指数, C_i是第i级的土地利用程度分级面积百分比。L_s是土地利用强度变化速率, L_t1和 L_t0分别是研究末期和初期的土地利用强度。对于土地利用分级指数, 庄大方等^[34]基于土地利用极限给出了各土地利用类型的程度的分级指数, 其后部分文献[29, 35]均基于研究区特征对该指数有适当调整。本文基于Feng等^[29]调整结果和所划定的土地利用类型, 对土地利用强度进行分级赋值, 结果见表 2。

土地利用转移矩阵分析方法是常见的土地利用变化分析方法, 最早可见于Singh^[36], 转移矩阵分析可以直接反映出土地的来源与去向, 全面显示区域土地利用的变化特征^[37], 计算方法见式(5)。

$\boldsymbol{P}_{i j}=\left[\begin{array}{ccc} P_{11} & \cdots & P_{1 n} \\ \vdots & \ddots & \vdots \\ P_{n 1} & \cdots & P_{n n} \end{array}\right], $

(5)

式中, P_ij为状态转移矩阵, 其中元素为不同土地利用类型转移的面积或比例。

地理加权回归分析方法(GWR)是由Fotheringham^[38]等基于局部回归分析和变参数的设想提出的处理空间异质性的主要手段之一, 计算方法见式(6)。

${y_i} = {\beta _0}\left( {{u_i}, {v_i}} \right) + \sum\limits_{k = 1}^p {{\beta _k}\left( {{u_i}, {v_i}} \right)} {x_{ik}} + \varepsilon, $

(6)

其中, y_i是第i点的预测值, (u_i, v_i)是第i点的位置, β₀(u_i, v_i)是第i点的回归常数, β_k(u_i, v_i)是第k个变量在第i点的回归参数, x_ik是第i点第k个变量的值, p是变量总数, ε是不确定性误差。β_k(u_i, v_i)由周围一定样本数量根据权重核函数确定。本文选择应用二次核函数(W_ij)确定每一个样本点的回归参数, 见式(7)。

${W_{ij}} = \left\{ {\begin{array}{*{20}{c}} {{{\left[ {1 - {{\left( {\frac{{{d_{ij}}}}{b}} \right)}^2}} \right]}^2}}&{\left( {{d_{ij}} < b} \right)} \\ 0&{\left( {{d_{ij}} > b} \right)} \end{array}} \right., $

(7)

其中, d_ij是参与局部回归的样本点与该点的距离。b为参与局部归回的样本点的最远距离, 又称为带宽, 可人为设定, 也可选择自适应带宽并由模型迭代计算确定。本文选择自适应带宽, 并以AIC_c准则作为模型计算带宽确定方法, 具体及可见文献[39]。

3 研究结果 3.1 土地利用空间分布及动态变化

改革开放初期, 土地面积较大且分布广泛的类型为耕地、浅海海域、养殖盐田用地、滩涂、建设用地, 其初始占比分别达46.76%、15.56%、8.76%、8.71%、7.59%; 未利用地集中分布在北大港、南大港及曹妃甸湿地等地区, 林地及草地主要分布在北部燕山山地, 内陆水体主要沿河流分布, 待利用地占比极小且仅分布在天津港附近。1980年以来, 耕地、海域、滩涂、未利用地一直处于减少趋势, 其中未利用地和滩涂降幅最大, 年变化率达–1.60%及–1.29%, 耕地面积缩减量最大, 达1 056.56 km²。养殖盐田用地面积在研究初期出现快速增加, 2000年之后开始不断减少, 累计变化率达1.07%。建设用地面积年变化率达2.46%, 且累计增量为各类型中最高, 达1 227.15 km²。待利用地2000年以来增幅明显, 新增区域集中分布于黄骅、滨海新区及曹妃甸; 研究期内年变化率最高, 达409.16%, 且其累计增量仅次于建设用地。其他利用类型中, 内陆水体处于持续减少状态, 林地及草地近年来由减少转为增加, 年变化率相对较低。

土地利用动态度分析显示, 近40 a以来冀津海岸带地区总动态度为0.48%。不同时段的动态度差异显著, 2000—2010年达0.83%, 为各时段最高; 1990—2000年最低, 为0.33%; 其他年份介于0.5%~0.6%。

为进一步了解土地利用动态度变化的空间差异, 分析研究区域1 km×1 km格网内土地利用动态度, 并依据其动态度值域范围及自然断裂法, 将各时段动态度结果分为无变化(0%)、极慢(1%~3%)、慢速(3%~6%)、快速(6%~8%)、极快(8%~10%)5个级别。结果见图 3。极快变化区集中分布在距离海岸带3~8 km范围内。其中1980—1990年极快变化区占2.90%, 主要分布在沧州沿海、滦河口沿海地区。1990—2000年极快变化区面积仅占1.28%, 零星分布于海河河口、南大港湿地及曹妃甸沿海。2000—2010年研究区域极快变化区面积高达4.18%, 广泛分布于滨海新区、曹妃甸等地区及乐亭、黄骅、秦皇岛市区沿海。2010年以来, 研究区域极快变化区下降至1.99%, 零星分布在海海地区及滨海新区、曹妃甸区、丰南区等地区。

3.2 土地利用类型转移

近40 a来, 冀津海岸带地区累计发生利用方式转移的面积占研究区域的18.22%, 各土地类型之间的转换情况见表 3。主要的土地利用主要转移特征可汇总为3大类: 耕地、养殖盐田用地、滩涂转为建设用地(6.20%), 主要位于滨海新区、曹妃甸、乐亭、秦皇岛市区及黄骅; 耕地、未利用地、滩涂转为养殖盐田用地(5.61%), 主要位于黄骅、海兴、曹妃甸; 滩涂、海域、养殖盐田用地转为待利用地(3.80%), 主要位于滨海新区、曹妃甸、黄骅。

土地利用类型转移特征具有时段差异性(图 4)。2000—2010年土地利用类型发生转移的面积最大, 达8.22%, 主要为耕地、养殖盐田转建设用地及滩涂、海域转待利用地。1980—1990年土地利用类型发生转移的面积占比达5.25%, 主要为耕地、滩涂及未利用地转为养殖盐田用地。其他时间段土地利用类型发生转移的面积均低于5%, 其中1990—2000年主要的转移类型为耕地转为养殖盐田用地, 2010—2018年主要的转移类型包括耕地、养殖盐田用地、待利用地转为建设用地及养殖盐田用地、海域转为待利用地。

临近海岸线的部分区域发生过两次或以上的土地利用类型转换(图 4)。其中发生过两次土地利用转换的区域总面积达453.11 km², 占比达2.63%, 主要转换类型包括耕地→养殖盐田用地→建设用地、未利用地→养殖盐田用地→建设用地、未利用地→耕地→养殖盐田用地、滩涂→养殖盐田用地→待利用地、未利用地→养殖盐田用地→待利用地。发生过3次土地利用转换的区域主要位于曹妃甸和乐亭, 总面积达32.23 km², 转换类型为滩涂、未利用地、耕地→养殖与盐田用地→待利用地→建设用地。

3.3 土地利用强度及其变化

通过土地利用强度指数及其变化分析方法, 分别分析了研究区域各时期总体及分格网的土地利用强度及其变化, 结果见图 5。就全部研究区域而言, 各时期土地利用强度分别为247.67、253.88、258.21、274.86、283.27, 总体呈现出持续上升的趋势, 其中2000年后利用强度出现显著增加。

借助3.2节生成的空间格网, 进一步分析了研究区域土地利用强度变化时空差异及变化特征(图 5)。结果显示, 研究区域土地利用强度变化显著的区域在陆域集中于滨海新区、秦皇岛市区, 在海域范围内集中于滨海新区、曹妃甸、乐亭、黄骅等地区, 其中土地利用强度出现3个层级及以上的变化区域占研究区域的3.93%, 主要位于黄骅、滨海新区、曹妃甸围填海区域及乐亭沿海陆域, 港口或临海工业区的建设使这些区域出现从难利用土地向最高强度的建设用地变化。

时间变化上, 1980—2000年, 土地利用强度变化高值区面积仅占0.01%, 且多属于零星分布, 主要的强度变化表现为从难利用土地向中等利用强度变化。2000年以来, 土地利用强度变化高值区明显增加, 强度变化级别上表现为从难利用、中等强度向高强度变化; 其中2000—2010年土地利用强度显著变化区集中于曹妃甸、滨海新区及黄骅港围填海区域, 2010年以来集中分布于滨海新区。

3.4 主要土地利用类型二级分类变化

建设用地、待利用地、养殖及盐田用地是最主要的土地利用转移目标。在中科院土地利用分类系统的基础上, 结合土地利用实际特征, 对主要土地利用类型进行二级分类调整及分析。将建设用地划分为城镇用地、农村居民点用地、工矿用地、陆域交通用地、港口堆场用地, 将养殖与盐田用地划分为养殖用地、盐田用地, 将待利用地划分为待利用土地和待利用围海。各二级分类累计增幅、年变化率及时空变化情况见表 4及图 6。

结果显示, 待利用土地、工矿用地、养殖用地及城镇用地在研究期内显著增长, 增量均在300 km²以上; 待利用地、港口堆场用地、工矿用地增长较快, 年变化率均在10%以上, 其中待利用地甚至接近400%。时间变化上, 近年来养殖用地及盐田用地面积近年来都出现下降, 农村居民点用地相对平稳, 面积增长有限; 城镇用地2000—2010年出现显著增加, 2010年以来增量出现回落。工矿用地、待利用土地2000年以来仍处于持续增长之中。

空间变化上, 沧州市城镇用地、农村居民点用地、交通用地、及盐田用地年变化率居首位, 其中盐田用地净增量为各地市最高。秦皇岛市城镇用地、养殖用地年变化率居首位, 其中城镇用地净增量为各地市最高。唐山市工矿用地、交通用地年变化率居研究区域首位, 港口用地、陆域交通用地、养殖用地增量为各地市最高。天津市工矿用地、待利用土地增量为各地市最高。

4 讨论 4.1 海岸带土地利用变化驱动因素定量分析

前述研究显示, 研究区域各类型土地显著变化时段为2000—2010年。受篇幅及自变量数据所限, 本文仅讨论2000—2010年陆域部分土地利用变化驱动因素。以数值连续且有意义的土地利用强度为因变量, 选择可以量化且弱共线性的参数为自变量, 使用地理加权回归模型(GWR)对驱动因素进行量化分析。自变量包括地形, 与河流、城市中心、县城中心、海岸线、港口、工业区、开放型交通线(国道与省道)、封闭型交通线(铁路与高速公路)的距离, 人口^[40], GDP^[41], 共计12个自变量。受人口、GDP等自变量数据格式所限, 本文将研究区域划分为1′为单位的网格, 计算各网格的因变量及自变量值并输入模型。结果显示, 2000年及2010年模拟结果的最优带宽均为72(即使用72个周围样本点参与局部回归), 调整后R²分别达0.837及0.836, 显示模型拟合效果较好, 各影响因素均值及中值见表 5。

表 5显示, 除回归常数之外, 回归系数绝对值较高的影响因子包括地形、与城市中心距离、GDP、人口及与港口距离, 回归系数绝对值的中值均在0.5以上, 显示了中心城市、社会经济、人口分布、港口建设对土地利用强度的强大作用力以及自然条件对于开发利用活动的限制。地形、回归常数、与城市中心距离、与工业区距离、与开放型及封闭性交通线距离、与港口距离以负相关为主, 与河流距离、人口、GDP以正相关为主。

2000年及2010年, 各驱动因素的回归系数出现了一定变化。地形、GDP、回归常数等因素的驱动效应相对减弱, 一方面显示了技术进步逐渐克服地形的限制, 另一方面也显示了2000年以来新城区规划、建设与经济发展水平之间不匹配性区域增多。人口、与港口距离等因子的驱动效应显著增强, 显示了人口聚集、港口建设是影响土地利用变化的刚性因素。

对于驱动效应显著增强的驱动因子, 将P≤0.05的网格区落于图内(图 7)。结果显示, 相比于2000年, 2010年人口、与港口距离、与城市中心距离、与工业区距离等因子所能解释的区域有所增多: 主城区以外地区的人口数量对制约土地利用强度的正相关效应增强, 天津港、京唐港和黄骅港对土地利用强度的驱动效应及影响面积不断扩大。

4.2 海岸带土地利用变化政策性驱动因素

自然资源管理及生态环境保护政策的实施是近年来土地利用变化的重要影响因素。农耕地保护方面, 2009年国家开始实施永久基本农田保护制度^[42], 有效保护了基本农田, 降低了耕地的年动态度变化幅度。滨海滩涂保护方面, 2011年国家发改委及国家海洋局出台了围填海计划管理办法^[43], 对各沿海省份每年的新增围填海面积进行指标限定, 2017年国家海洋局出台《关于进一步加强渤海生态环境保护工作的意见》^[44], 暂停渤海海域围填海; 2018年初国务院文件《关于加强滨海湿地保护严格管控围填海的通知》^[45]进一步作出在全国范围内严禁新增围填海的决定, 这些决定和政策使得2010年以来冀津沿海地区滩涂转换为其他用地类型的面积显著减少。

土地财政依赖性强、土地开发后继乏力是2000年来待利用地激增的主要原因。土地出让是政府部门的重要财政收入来源。近年来, 地方政府为增加财政收入, 往往编制不符合当地发展条件和规律的城镇规划, 并借此多征地、多占地, 将大量土地提前纳入储备计划^[46], 或在沿海地区超前推动大规模填海造陆。而由于部分地区产业发展、城市发展及人口增长并没有随时间变化出现城镇规划所设计的增长情景, 如曹妃甸围填海区、黄骅围填海区、滨海新区的汉沽及大港等地, 导致城市土地开发后继乏力, 储备土地及填海土地大量闲置。

5 结论

改革开放以来, 冀津海岸带地区土地利用动态度达0.48%, 土地利用强度也由247.67上升至283.27, 土地利用逐渐向养殖盐田用地、待利用地及建设用地转换。土地利用空间变化主要出现在滨海新区、曹妃甸等地区及乐亭、黄骅、秦皇岛市区沿海, 且集中在距离海岸带3 km~8 km范围内。时间变化方面土地利用动态度峰值、强度增长峰值区出现在2000—2010年。海岸带土地利用变化主要受到地形、中心城市、GDP、人口、港口、政策等因素影响, 其中人口与GDP为强正相关因子, 地形、与中心城市距离、与港口距离为强负相关因子。相比于2010年, 地形、GDP、回归常数等因素的驱动效应相对减弱, 而人口、与港口距离等因子的驱动效应显著增强。政策性因素方面, 土地财政依赖性强、过于超前的规划及后续开发乏力导致待利用地迅速增长, 而2010年以来基本农田保护制度、滨海湿地保护及严禁围填海的政策的实施有效减弱了耕地及滩涂面积的下降趋势。

为推动冀津海岸带地区土地资源高效及可持续利用, 保护沿海地区海洋生态环境, 当前应当进一步提高待利用地使用效率, 进一步加强耕地及滨海滩涂资源保护, 加强土地利用高强度变化区域的生态环境监测与修复。加强对囤积土地等行为的打击力度, 加强待利用地的政府收储管理, 鼓励新增建设项目优先布局在待利用土地, 对长期闲置的待利用土地试行梯级开发或就地建设为临时绿地。严格执行耕地及滨海湿地保护政策, 控制永久基本农田、生态红线区及自然保护区附近的土地利用变化率及强度增加趋势, 在重要河口、重要湿地分布区推行退养还海还滩政策。对于由难利用地或低强度利用土地直接转变为建设用地的区域, 如曹妃甸、滨海新区、黄骅及乐亭沿海地区, 应当加强陆域和海域生态环境的调查监测, 识别富营养化、重金属污染、湿地退化等问题, 及早对破损的生态系统和受污染的环境进行修复。