土地利用/覆被变化(land use and land cover change, LUCC)是人类社会经济活动行为与自然生态过程交互和链接的纽带, 直接反映人类活动与自然环境的相互作用强度与影响特征^[1-2]。国际全球环境变化人文因素计划(International Human Dimensions Programme, IHDP)和国际地圈生物圈计划(International Geosphere-Biosphere Programme, IGBP)于1995年联合提出了LUCC计划, 成为土地利用/覆被变化研究领域的历史性转折点; 于2006年提出的全球土地计划(Global Land Plan, GLP)将LUCC研究推进入快速发展阶段。在国际综合评估中, LUCC研究可概括为变化过程、变化驱动力机制和模型开发与应用3个核心问题^[3-4]。其中, 土地利用变化过程监测是LUCC研究的基础, 多是基于连续的多期土地利用数据, 分析区域土地利用变化的特征与机制; 土地利用变化情景分析和模拟则是LUCC研究的热点, 国内外学者根据研究的需要选用不同的模型对流域、矿区、城市等类型区域的土地利用变化进行模拟, 常用的方法有马尔可夫(Markov)模型、系统动力学(Systems Dynamics, SD)模型、元胞自动机(cellular automata, CA)模型、土地利用转化及其效应(conversion of land use and its effects, CLUE)模型、多智能体(Multi-agent System, MAS)模型以及混合模型等^[5-8]。

随着人类活动与自然环境之间相互作用的加剧, 越来越多的研究聚焦于人地关系矛盾尖锐、生态环境脆弱、人口与城市化快速发展的地区^[9]。海岸带是陆地和海洋之间的过渡地带, 是典型的生态交错带和脆弱区, 也是开发利用强度最高的区域之一^[10]。在我国, 主要是20世纪90年代以来, 随着社会经济持续快速发展, 沿海区域率先进入快速城市化和工业化阶段^[11], 与此相应, 海岸带土地利用呈现出新的变化特征和趋势, 因而成为新时期区域发展研究的重要部分。黄河三角洲—莱州湾海岸带区域包含河口三角洲和海湾海岸带2大地貌类型, 在长期的海陆交互作用基础上, 加上经济快速发展的刺激, 致使土地利用类型复杂多样、变化迅速, 是海岸带土地利用变化研究的代表性和热点区域。已有很多学者针对黄河三角洲和莱州湾区域的土地利用/覆盖变化开展了研究, 例如:叶庆华等对黄河三角洲土地利用时空进行复合变化图谱分析^[12]; 张成扬和赵智杰对黄河三角洲土地利用/覆盖时空变化特征与驱动因素进行了定量分析^[13]; 张宝雷等基于复杂网络对黄河三角洲进行土地利用变化过程研究^[14]; 齐善忠和蔡文华对莱州湾海岸带土地利用变化及环境灾害进行了研究^[15]等。这些研究较好地揭示了黄河三角洲、莱州湾海岸带土地利用变化的历史特征和基本规律, 揭示并指出建设用地面积和比例快速增加、自然湿地有增有减但总面积和比例不断减少是土地利用/覆盖变化的最基本特征, 河口三角洲发育、工业化、城市化、围填海和农业发展等是土地利用变化的主要原因, 但将黄河三角洲和莱州湾海岸带作为一个整体并对该区域土地利用演化过程及未来趋势进行分析和模拟的研究则未见报道, 因此, 在前人研究进展的基础上, 本文将黄河三角洲和莱州湾海岸带作为一个整体, 强调陆海兼顾基于遥感技术监测土地利用变化的历史特征, 并利用Markov方法和基于CA技术的未来土地利用模拟(Future Land Use Simulation, FLUS)模型对土地利用的未来变化趋势进行多情景分析和模拟, 以期为黄河三角洲—莱州湾海岸带区域土地资源的可持续利用及海岸带综合管理决策提供科学依据。

1 研究区域概况

黄河三角洲—莱州湾海岸带地处渤海南部、山东半岛北部(图 1), 包括黄河三角洲和莱州湾海岸带2个地貌单元, 介于117°52′~120°74′E、36°28′~38°28′N。西北部为黄河三角洲, 自1855年黄河改道以来快速淤积, 河口湿地面积不断增加, 成为世界上增长最快的三角洲; 东南部为莱州湾海岸带, 有支脉河、小清河、弥河、白浪河、虞河、潍河、胶莱河等注入, 海湾东部濒临胶东半岛西北部山地, 沿其山前堆积平原发育砂质海岸。黄河三角洲—莱州湾海岸带处于海陆交互的生态过渡带, 生态环境较为脆弱, 抗干扰能力弱, 存在咸水倒灌、海岸侵蚀、盐碱化、风暴潮等灾害。本文以陆域市县行政边界和海域–6 m等深线确定研究区的边界, 包括无棣县、沾化县、河口区、利津县、垦利县、东营市辖区、广饶县、寿光市、潍坊市辖区、昌邑市、平度市、莱州市、招远市和龙口市以及–6 m等深线覆盖的浅海水域, 总面积约2.86万km²。该区域涉及黄河三角洲高效生态经济区和山东半岛蓝色经济区2大战略规划, 经济快速发展, 人口不断增加, 工农业较为发达, 人类活动越来越强烈。

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

基于2000年1︰10万中国土地利用数据库^[2]和Landsat TM卫星影像, 通过分类系统调整和图斑修改, 获得2000年中国海岸带土地利用数据, 进而, 通过动态更新依次获得2005年、2010年和2015年中国海岸带土地利用数据^[16-17]。该数据集将海岸带土地利用分为8个一级类型、24个二级类型^[18], 从中提取黄河三角洲—莱州湾海岸带区域4个年份的数据用于本研究; 为便于研究以及突出区域特征, 将数据重分类为耕地、林地、草地、建设用地、内陆水体、滨海湿地、浅海水域、人工湿地和未利用地共9个类型。土地利用分布与变化受到多重因素的驱动和影响, 根据已有研究的经验总结^[19-20]和研究区实际及数据可获取性, 情景分析与模拟选取6个驱动力因子, 数据源及其处理如表 1所示。

2.2 研究方法 2.2.1 土地利用变化幅度

土地利用变化幅度主要反映不同土地利用类型面积总量的变化^[22], 表达式为:

$ R = \frac{{{S_2} - {S_1}}}{{{S_1}}} \times 100\% , $

(1)

式中, R表示研究时段内某一地类变化幅度; S₁、S₂分别表示时段初期及末期某一土地利用类型面积。

2.2.2 土地利用变化速率

以单一土地利用动态度定量表示一定时间范围内某种土地利用类型变化的速度^[23]。正值表示该土地利用类型面积增加, 反之则减小; 具体衡量时使用绝对值, 绝对值越高, 表明该土地利用类型转换数量多, 变化相对剧烈, 反之则相对稳定。表达式为:

$ V = R \times \frac{1}{T}, $

(2)

式中, V表示研究时段内某一地类动态度; T表示研究时段。

2.2.3 土地利用程度综合指数

土地利用程度综合指数反映土地利用的广度和深度, 用于定量表征土地利用的集约化水平。将土地利用程度按照土地自然综合体在社会因素影响下的状态分为1级(未利用地、浅海水域)、2级(林地、草地、内陆水体、滨海湿地)、3级(耕地、人工湿地)和4级(建设用地)4种土地利用级别, 分别赋予1、2、3、4不同的分级指数, 并以土地利用类型面积比例作为权重, 加和得到100~400之间连续分布的土地利用程度综合指数, 体现人与自然交互的综合效应^[24], 表达式为:

$ I = 100 \times \sum\limits_{i = 1}^n {{A_i}} \times {C_i}, I \in [100, 400], $

(3)

式中, I为土地利用程度综合指数; A_i为第i级土地利用程度分级的分级指数; C_i为第i级土地利用程度分级的面积百分比; n为土地利用程度的分级数。

2.2.4 土地利用转移矩阵

土地利用转移矩阵来源于系统分析中对系统状态与状态转移的定量描述^[25], 在土地利用变化分析中被广泛使用, 可以直接反映出土地的来源与去向, 全面显示区域土地利用的变化特征, 表达式为:

$ {\mathit{\boldsymbol{P}}_{ij}} = \left[ {\begin{array}{*{20}{c}} {{P_{11}}}& \cdots &{{P_{1n}}}\\ \vdots & \ddots & \vdots \\ {{P_{n1}}}& \cdots &{{P_{nn}}} \end{array}} \right], $

(4)

式中, P_ij为状态转移矩阵, 矩阵元素是土地利用转移面积或转移比例。

2.2.5 Markov模型

Markov过程是一种具有“无后效性”的特殊随机运动过程^[26]。Markov模型是一种基于转移概率的数学统计模型, 土地利用的动态演化过程具有较明显的马尔科夫性质, 因而可以利用土地利用类型的转移概率来对土地利用的面积进行预测^[27]。设S_t是一随机过程, 将某一时刻的土地利用类型对应于Markov过程的可能状态, 土地利用类型之间相互转换的面积比例作为转移概率。其公式为:

$ {\mathit{\boldsymbol{S}}_{t + 1}} = {P_{ij}} \cdot {\mathit{\boldsymbol{S}}_t}, $

(5)

式中, S_t、S_t+1分别为t、t+1时刻的系统状态; P_ij为转移概率矩阵。

2.2.6 FLUS模型

FLUS模型是用于模拟人类活动与自然影响下的土地利用变化以及未来土地利用情景的模型^[28]。该模型源自CA, 但做了较大的改进:采用能够有效处理非线性关系的神经网络算法(Artificial Neural Network, ANN)从一期土地利用数据与包含人为活动与自然效应的多种驱动力因子获取各类用地类型的适宜性概率。ANN的优点是特别适用于模拟复杂的非线性系统, 显著优于多准则判断等传统分析方法^[29]; 从一期土地利用分布数据中采样, 能较好的避免误差传递的发生^[30]; 基于CA模型进行用地空间配置, 模拟过程采用基于轮盘赌选择的自适应惯性竞争机制, 能有效处理多种土地利用类型在自然作用与人类活动共同影响下发生相互转化时的不确定性与复杂性, 具有较高的模拟精度并且能获得与现实土地利用分布非常相似的结果^[31]。模型具体使用BP-ANN方法, 是一种多层前馈神经网络, 分为输入层、隐藏层和输出层3层, 用于训练和评估每个栅格土地利用类型发生的概率^[32]。

基于自适应惯性机制的元胞自动机, 对邻域影响的计算公式为:

$ \mathit{\Omega }_{p, k}^t = \frac{{\sum\nolimits_{N \times N} {{\mathop{\rm con}\nolimits} } \left( {c_p^{t - 1} = k} \right)}}{{N \times N - 1}} \times {w_k}, $

(6)

式中, Ω_{p, k}^t表示邻域影响程度, $\sum\nolimits_{N \times N} {{\rm{con}}\left( {c_p^{t - 1} = k} \right)} $表示在N×N的Moore邻域窗口, 上一次迭代结束后第k种地类的栅格总数, w_k为各类用地的邻域作用的权重。

自适应惯性系数用于判断特定土地利用类型的发展趋势与实际需求是否存在较大差距, 如存在, 则在下一次迭代中调整该土地利用的发展趋势, 从而动态控制该土地利用类型的数量^[31], 其表达式为:

$ I_k^t = \left\{ {\begin{array}{*{20}{l}} {I_k^{t - 1}}&{{\rm{ if }}\left| {D_k^{t - 2}} \right| \le \left| {D_k^{t - 1}} \right|}\\ {I_k^{t - 1} \times \frac{{D_k^{t - 2}}}{{D_k^{t - 1}}}}&{{\rm{ if }}0 > D_k^{t - 2} > D_k^{t - 1}}\\ {I_k^{t - 1} \times \frac{{D_k^{t - 1}}}{{D_k^{t - 2}}}}&{{\rm{ if }}D_k^{t - 1} > D_k^{t - 2} > 0} \end{array}} \right. $

(7)

式中, I_k^t表示k种地类在t时刻的自适应惯性系数, D^t–1_k、D^t–2_k分别为t–1、t–2时刻第k种用地类型的栅格数与需求数量之差。

2.2.7 情景分析与模拟

情景分析是对未来态势进行探测的常用方法, 考虑多方面的因素从而预测多种可能情景, 提供更加科学的参考依据^[33]。根据研究区土地利用和经济社会发展的现状特征, 设置趋势延续(Business as Usual, BAU)、经济社会发展(Social and Economic Development, SED)和生态保护优先(Ecological Protection Priority, EPP)三种情景, 将目标年份设定为2025年, 利用2005—2015年间的转移概率矩阵, 基于Markov数学方法实现3种情景的总量计算:直接预测得到BAU情景, 修改转移概率矩阵和地类转换限制矩阵, 进而预测得到SED和EPP情景。多情景空间模拟使用FLUS模型, 选取均匀采样模式, 基于2015年土地利用现状和土地利用变化驱动因子, 将各因素进行归一化, 提取0.1%的栅格样本进行ANN样本训练及适应性概率计算。

3 结果与分析 3.1 土地利用数量变化特征

统计4个年份的土地利用面积, 如表 2所示; 以2015年为例展示土地利用的宏观格局特征, 如图 2所示; 土地利用面积净变化及变化速率特征分别如表 3、图 3所示。

由表 2和图 3可知, 黄河三角洲—莱州湾海岸带区域土地利用类型丰富:主要类型是耕地, 约占50%;其次为建设用地, 占12.53%, 与耕地镶嵌, 分布较为离散; 浅海水域面积位列第3, 占11.60%;人工湿地的面积比例为8.55%, 主要连片分布在海岸带区域; 草地面积比例是5.19%, 主要分布于黄河三角洲及东部山地区域; 滨海湿地、内陆水体、林地和未利用地的面积比例较小, 林地主要分布在东南部山地, 未利用地主要分布在沿海区域, 包括盐碱地、沼泽地等。黄河三角洲—莱州湾海岸带土地利用由海向陆总体呈现出由“浅海水域→滨海湿地→人工湿地→陆域多类型混合区域”过渡的宏观格局特征。

由表 3、图 3可以看出, 建设用地和人工湿地面积净增加趋势显著, 相反, 草地和滨海湿地面积净减少趋势显著, 这表明黄河三角洲-莱州湾海岸带区域的城镇化、工业化以及养殖业等处于快速发展的过程中。各地类变化幅度的绝对值介于0.78%~ 49.94%之间, 类型间的差异性以及时间阶段的差异性均较为显著: 2000—2005年, 人工湿地增长和草地减少最显著; 人工湿地的变化幅度高达49.94%, 变化速率为9.99%, 而草地则分别为–24.32%和–4.86%; 2005—2010年, 建设用地和人工湿地的增长幅度在25%左右, 而滨海湿地的减少幅度则高达21.87%; 2010—2015年, 建设用地和内陆水体分别增加了9.29%和15.04%, 而滨海湿地和未利用地则分别减少了6.89%和11.21%。综合三个时段, 人工湿地和建设用地增长幅度大, 相反, 草地、滨海湿地和未利用地减少幅度大。

各个地类不同时段的变化趋势差异显著:耕地和林地在第1时段内小幅增加, 后2个时段内则小幅度减少; 草地持续减少, 变化幅度由24.32%降低至4.03%;建设用地持续扩张, 变化幅度在9%~25%;内陆水体在前2个时段以9%的幅度缩减, 但在第3时段扩增了15.04%;滨海湿地面积减少的幅度在6.89%~21.87%;浅海水域15 a间持续小幅度缩减; 人工湿地在第1时段变化幅度为49.94%, 第2时段降为25.07%, 第3时段骤降至1.98%, 表明前2个阶段是海岸带盐业和养殖业大发展的时期。综合各地类近15 a的变化幅度, 耕地、林地和浅海水域变化幅度相对较小, 而其他地类较为剧烈; 3个时段相比, 第1时段整体的变化幅度较为剧烈, 后2个时段则逐渐降低, 土地利用趋于稳定。

计算各个年份的土地利用程度综合指数, 由于1级地类(未利用地、浅海水域)和2级地类(林地、草地、内陆水体、滨海湿地)面积总体减少, 3级地类(耕地、人工湿地)和4级地类(建设用地)面积总体增加, 2000年、2005年、2010年和2015年土地利用程度综合指数分别为255.79、261.93、266.10、268.82, 呈较为显著的递增趋势。

3.2 土地利用面积转移特征

2000—2005年、2005—2010年、2010—2015年三个时段的土地利用面积转移矩阵如表 4所示。各个时段耕地主要转变为建设用地, 而其获得补充的主要来源有草地、内陆水体和未利用地。林地与草地互换以及向耕地、建设用地的转变较为显著, 但不同阶段的转换方向及规模差异较大。草地主要转为耕地、人工湿地, 其获得补充的主要来源是耕地、内陆水体、滨海湿地和未利用地。建设用地和人工湿地的扩张显著, 前者主要占用耕地、滨海湿地和人工湿地, 后者主要占用滨海湿地。内陆水体在第1时段大量变为人工湿地, 后2个时段主要转变为耕地、草地、人工湿地和建设用地。滨海湿地大量转为人工湿地和建设用地。浅海水域在前2个时段主要转变为滨海湿地(河口沉积、水域变浅), 在第3时段主要转变为建设用地、滨海湿地和未利用地。上述不同地类之间相互转移的可能性及规模特征体现了自然-经济社会耦合系统对土地利用变化的影响和制约, 总体上, 受经济发展驱动, 土地利用转移呈现出由自然类型转向人工类型以及人工化程度和集约化水平梯级递升的特征。

3.3 土地利用多情景模拟

至2025年土地利用多情景分析与模拟的结果如表 5和图 4所示。

3.3.1 BAU情景

在BAU情景下直接基于2005—2015年土地利用概率转移矩阵建立Markov模型, 所以延续2005— 2015年的驱动力和驱动机制, 模拟结果表明:与2015年相比, 至2025年耕地、林地、草地和湿地的减少态势以及建设用地和人工湿地的扩张态势较为显著, 土地利用程度综合指数将达到275.08, 土地利用的破碎化程度也将愈加显著。

3.3.2 SED情景

该情景以持续满足快速的经济社会增长对土地资源的需求为基本特征, 其转移概率矩阵是对2005—2015年转移概率矩阵进行一定的修改而获得:提高耕地转为建设用地和人工湿地的概率以及草地转为建设用地的概率; 降低林地转向草地的概率; 限制建设用地转变为其他类型; 提高滨海湿地转为人工湿地、建设用地的概率; 提高未利用地转为人工湿地的概率。模拟结果表明:至2025年, 建设用地和人工湿地将大幅扩张, 滨海湿地则大幅缩减, 与2015年现状相比, 3种类型的变化率分别达38.52%、27.61%和–25.27%;格局方面, 沿海区域土地利用的人工化和破碎化特征极为显著; 土地利用程度综合指数将达到278.15, 在3种情景中是最高的, 而且比2015年高出3.47%。

3.3.3 EPP情景

该情景以减缓自然型地类转向人工型地类、保证生态用地的面积、优化生态结构为导向, 其转移概率矩阵是对2005—2015年转移概率矩阵进行如下修改而获得:耕地转为林地的概率加倍, 转向建设用地的概率减半; 林地转向耕地和草地的概率减半, 且限制转为建设用地和人工湿地; 草地转向耕地和建设用地的概率减半, 转为林地的概率加倍, 限制转为人工湿地; 限制建设用地转为其他地类; 内陆水体转为草地的概率减半, 限制转为建设用地; 滨海湿地不得转为建设用地和人工湿地; 浅海水域转为建设用地的概率减半; 人工湿地转为草地、浅海水域和湿地的概率加倍; 未利用地转出概率加倍。模拟结果表明:该情景下, 至2025年土地利用相对于前两种情景变化幅度较小, 草地、林地、湿地等生态用地面积萎缩的速率明显减缓, 未利用地得到较高程度的开发利用, 土地利用程度综合指数为274.92;总体上, 该情景较好地达到了生态保护的目标。

4 结论与讨论

利用遥感技术、土地利用变化指数模型等方法, 对黄河三角洲—莱州湾海岸带2000—2015年间的土地利用变化进行定量分析, 进而, 利用Markov方法和FLUS模型对至2025年的土地利用结构和格局进行多情景分析与模拟。主要结论如下:

(1) 黄河三角洲—莱州湾海岸带区域耕地面积最大, 其次是浅海水域、湿地(内陆水体、滨海湿地、人工湿地)以及建设用地; 由海向陆土地利用总体呈现出由“浅海水域→滨海湿地→人工湿地→陆域多类型混合区域”过渡的宏观格局特征。

(2) 2000—2015年, 耕地、林地、草地、湿地等生态型地类的面积和占比显著减少, 而建设用地和人工湿地的面积与占比则显著增加, 土地利用的集约化水平和土地利用程度综合指数不断提高; 三个5 a时段土地利用变化的数量特征、主导类型、变化速率等存在较大差异, 区域整体土地利用的变化幅度逐渐降低, 表明土地利用变化的驱动因素和驱动机制存在一定的阶段性特征, 各种规划和管控性政策措施的作用显现, 尤其是在研究期末, 大规模围填海等发展趋势已经明显“降温”。

(3) 土地利用多情景分析与模拟结果表明, 至2025年, 土地利用仍以耕地、建设用地和浅海水域为主, 建设用地和人工湿地的扩张仍然较为显著, 但不同情景的扩张强度及途径有较大的差异: BAU情景和SED情景之间的共性较为显著, 但SED情景的变化强度和速率等更为显著, 因而也对应最高的土地利用程度综合指数; EPP情景充分发掘了区域丰富的未利用地资源的潜力和优势, 同时强调了建设用地的优化利用, 因而耕地、草地、林地、湿地等生态型用地的萎缩程度在3种情景中是最低的, 总体上较好的促进了生态保护和经济发展之间的平衡。

黄河三角洲—莱州湾海岸带是中国海岸带中非常具有典型性和代表性的区域, 主要是20世纪90年代以来, 该区域城镇化、工业化以及农业现代化的进程加快, 加之气候变化、海平面上升、海咸水入侵等自然因素的影响, 海岸带土地利用进入快速和急剧变化的时期, 突出表现为耕地大量流失、建设用地急剧扩张、大规模围填海以及湿地资源锐减等特征。而且, 进入21世纪以来, 经济发展等人类活动因素的影响在广度和深度上均愈来愈显著, 已逐渐成为黄河三角洲—莱州湾海岸带土地利用变化的主控因素:海岸带优越的地理区位吸引着大量人口、产业在此聚集, 集聚效应导致土地资源配置和需求量急剧上升, 需要大量新增土地作为人口容纳和产业发展的承载空间, 主要通过大量占用耕地和生态用地而得到满足, 但由于严格的耕地占补平衡政策的约束作用, 又刺激了大量滨海湿地资源的开垦和占用; 而且, 沿海居民生活水平提高、海洋经济发展迅速、产业结构调整等也大大推动了海岸带区域的围填海热潮、房地产业和旅游业的发展, 从而进一步加剧湿地资源的流失和破坏。

黄河三角洲—莱州湾海岸带在土地利用变化过程中, 资源环境约束问题逐渐显现, 尤其是自然岸线和湿地破坏所带来的一系列环境和生态效应问题, 生态保护和经济社会发展之间的矛盾日趋尖锐。这种问题及趋势在我国其他海岸带区域也较显著, 因此, 中央及地方相关政府部门也通过出台有针对性的政策法规等途径加强对沿海区域的管控, 例如, 2016年6月国家海洋局印发《关于全面建立实施海洋生态红线制度的意见》, 2018年7月国务院印发《关于加强滨海湿地保护严格管控围填海的通知》等, 这些政策的严厉性大大超出了以往, 有望在海岸带生态保护方面取得显著的成效。未来时期, 针对海岸带区域, 建议重视陆海统筹、城乡统筹原则, 进一步优化土地利用结构和格局, 充分挖掘未利用土地的资源优势, 尤其是前期大规模围填海阶段已经成陆但尚处于闲置状态的土地资源, 合理划分生产、生活和生态空间, 分类开发和管控, 有效促进经济发展和生态保护之间的协调与平衡。