基于改进稳定映射法的土地利用/覆被变化轨迹分析

2020-08-27王世东冯正英张合兵

农业机械学报 2020年8期

王世东冯正英余洋,2 张合兵

(1.河南理工大学测绘与国土信息工程学院，焦作 454000； 2.中国地质环境监测院，北京 100081)

0 引言

近年来，国内外学者较多利用转移矩阵[1-2]和土地利用动态变化模型[3-5]对土地利用/覆被动态变化进行研究，分析土地利用和覆被的数量变化、程度变化及空间分布。上述研究大多是基于单期或相邻两期数据的对比揭示在多时间序列上土地利用变化的时空规律。构建土地利用变化轨迹模型不仅可以得到其时空变化特征，还可以通过各种轨迹分析方法[6-10]分析其转移和流向[11]，从轨迹过程中寻求变化规律，揭示时间序列上土地利用/覆被变化的时空动态演变特征[12]。

土地利用/覆被变化研究是一个长期动态过程，以上研究较多只是刻画了一个相对“连续”时间尺度[13-14]的变化，若两期间隔时间较长，势必会割裂其连续变化，破坏变化规律中的过程完整性。因此，在遥感影像的支持下，研究基于连续时间尺度的土地利用/覆被变化轨迹，对于精细的土地利用动态监测及土地资源的可持续利用具有重要意义。当时间序列数据较多时，土地利用/覆被变化轨迹更加繁杂，轨迹提取和时空演变分析愈加困难，因此需要一种针对多时间序列数据的轨迹提取和分析方法。目前，稳定映射变化轨迹法已广泛应用于各领域的轨迹变化中，但该方法在土地利用变化方面的研究较少，尤其是针对多时间序列土地利用/覆被变化的研究。本文结合ArcGIS技术，基于长时序遥感数据和改进稳定映射法的土地利用/覆被变化(Land use/cover change，LUCC)轨迹分析方法，研究时相数与稳定映射变化轨迹STD指标的对应关系，为长时序连续时间尺度的土地利用变化轨迹提取和变化轨迹类型划分提供借鉴。

丹江流域水资源丰富，位于河南省淅川县的丹江口水库是南水北调中线工程的主要水源地[15]，自2003年南水北调工程实施以来，不仅改变了丹江口库区水资源的自然地理分布，同时受人类活动影响，丹江口库区的人口密度增大、土地开发利用强度增加、环境资源过度消耗等也使得丹江口库区的土地利用/覆被变化更加显著[16]。鉴于此，本文以丹江流域(河南段)为例，选用2002—2017年16期Landsat TM/OLI及HJ-1A CCD遥感影像数据，提取连续时间尺度的土地利用变化图谱，建立改进稳定映射法的LUCC轨迹分析方法，提取变化轨迹，并对轨迹类型进行划分，同时结合土地利用动态度等定量模型和景观指数，分析研究区16年间土地利用和覆被的转移流向总体特征和时空变化规律。

1 研究区概况

丹江流域(河南段)位于河南省西南部，地理坐标为东经110°35′～111°58′，北纬32°30′～34°1′，南接湖北省，西临陕西省，总面积约为8 419 km2。丹江干流自淅川县荆紫关入境，汇水总面积占整个研究区的3.11%，流经河南省洛阳市、南阳市和三门峡市，包含三门峡市的卢氏县，洛阳市的栾川县部分区域，南阳市的西峡县、淅川县、内乡县以及邓州市部分区域，其地理位置如图1所示。

图1 研究区地理位置图Fig.1 Geographical location map of study area

研究区具有典型的流域特色，整体地形起伏较大，地势由西北向东南倾斜，西北部环绕伏牛山脉，东南部依次为岗地和冲击平原区，属亚热带向暖温带过渡的季风性气候。研究区景观分布呈规律性，林地主要分布于研究区西部和北部的山地丘陵地带，大部分为块状山林，且受河流切割使得地貌起伏跌宕；耕地主要分布于东南部的平原地带；水域主要位于西南部淅川县境内的部分丹江口库区，是南水北调中线工程的主要水源地。2003年南水北调中线工程正式开工，2005年丹江口大坝开始实施加高加固，正常蓄水水位由157 m提升至170 m，库容由1.745×1010m3扩大至2.905×1010m3，至2014年底南水北调中线工程正式通水。南水北调中线工程实施及建成后，优化了我国水资源配置，有力支撑了水资源的社会经济发展，但对于丹江口库区的生态质量、土地资源可持续利用及社会经济发展都产生了重要影响。因此，以丹江流域(河南段)为例，基于长时序遥感数据，研究土地利用/覆被变化轨迹，揭示土地利用时空演变规律，为研究区土地可持续利用、生态环境建设以及水源地保护提供科学依据。

2 研究方法

2.1 数据来源与处理

采用单一遥感数据源往往难以支持长时间序列对土地变化轨迹分析的需求，为获得长时序多时相土地利用/覆被信息，选用2002—2017年共16期遥感影像数据，分别为Landsat TM(2002—2011年)、HJ-1A CCD(2012年)、Landsat OLI(2013—2017年)，其中选用Landsat TM影像的1～5波段和第7波段、HJ-1A CCD相机波段和Landsat OLI影像的1～7波段，从而保证所用影像的分辨率均为30 m，其他数据包括研究区2015年土地利用变更数据库和道路水系等。

为保证长时序遥感影像对空间叠加及土地利用分类的准确性，需要对影像进行辐射定标、校正、波段组合等预处理[17]，从而减少光照、大气等对不同时相影像造成的误差，同时也可以提高土地利用/覆被信息的解译精度。本研究在对长时序遥感影像进行预处理后，采用统计特征分析、主成分分析和相关分析等方法，利用相关系数矩阵、特征向量和特征值，按照最佳波段组合原则[18]，得到适合研究区土地利用分类的最佳波段组合：Landsat TM(5、4、3波段)、HJ-1A CCD(1、2、3波段)、Landsat OLI(7、5、4波段)。

根据土地利用一级类型划分体系和丹江流域(河南段)实际土地利用特征，将土地利用类型划分为水域、裸地、耕地、林地、草地、建设用地6类。结合土地利用变更数据库、Google历史影像及实地勘探，每期影像均选择两套样本点，分别用于分类和精度验证[19]。采用支持向量机分类和人工目视解译法，获取研究区16年间土地利用/覆被分类结果，如图2所示。经验证，分类总体精度不小于90.43%，Kappa系数达到了0.84，满足分类精度的要求。影像数据分类完成后，对地类属性做重分类处理，将水域、裸地、耕地、林地、草地、建设用地依次重编码为1、2、3、4、5、6，为后续研究区多时相土地利用变化图谱的提取提供数据准备。

图2 2002—2017年研究区土地利用分类图Fig.2 Land use status maps of study area

2.2 土地利用变化图谱的提取

土地利用变化轨迹的建立等同于多时相遥感数据土地空间单元变化图谱的提取，史培军等[1]提出了土地利用变化空间表达方法，变化图谱提取公式为

Ci×j=10nAki×j+10n-1A(k+1)i×j+10n-2A(k+2)i×j+…+ 10n-nA(k+n)i×j

(1)

式中Aki×j——k时刻研究区土地利用图谱

i——研究区行数

j——研究区列数

n——时相数

Ci×j——n+1个时相的土地利用变化图谱

变化图谱表现了由k时刻到k+n时刻的土地利用变化类型和空间分布，该方法适用于土地利用类型小于10的情况。

较多学者均采用土地利用变化空间表达方法，并结合ArcGIS技术的栅格运算功能，构建长时序相对“连续”的土地利用/覆被变化轨迹，时相数均为4～6个[20-22]。在本研究中，需要构建丹江流域(河南段)2002—2017年共16个时相的土地利用变化轨迹，变化图谱在影像栅格中表示为16位的属性值，例如某空间单元初始时相的土地利用类型为水域，并从始至终未发生变化，则其变化图谱的栅格属性值为1 111 111 111 111 111。但由于栅格运算功能只能使用系统字段进行计算，该字段值类型为长整型，取值范围为-2 147 483 648～2 147 483 648，最多只能满足10个时相的土地利用变化图谱的表达。因此，在本研究中将16个时相的土地利用分类结果通过转换工具转换为矢量，使得每一个土地空间单元对应一个矢量点，采用空间连接和字段计算器功能，代入式(1)，再将矢量点转换为栅格，此时栅格数据中参与运算的字段类型为字符串，字段值即为研究区16个时相的土地利用变化图谱信息。

2.3 改进稳定映射变化轨迹分析

AMY等[23]基于6时相23个随机地点的土地利用变化结果，通过计算研究区长时序多时相变化图谱中的STD指标，即研究时段内各土地空间单元所经历演变过程的变化次数(Turnover，T)、相似性(Similarity，S)和多样性(Diversity，D)来综合反映土地利用轨迹变化过程，从而提出了土地利用/覆被长时序多时相稳定映射变化轨迹分析方法。其中，S表示某土地空间单元在研究时段内经历的相同土地利用类型的数量；T表示研究时段内该土地空间单元在相邻时间段发生变化的次数；D表示某土地空间单元在研究时段内经历的不同土地利用类型的数量。

STD指标是稳定映射变化轨迹模型的基础，是划分研究区土地利用变化轨迹类型的关键。参考AMY等[23]提出的土地利用变化轨迹的STD判定原则，本研究结合长时序遥感数据变化图谱信息，推导计算出稳定映射STD指标与时相数的关系式，并根据研究区的实际情况，提出适合研究区的基于改进稳定映射法的LUCC变化轨迹分析方法。

2.3.1STD指标计算

土地利用变化轨迹的种类受研究时相数和土地利用类型数量的影响，令研究时相数为n，土地利用类型数为m，则可能出现的变化轨迹有mn种。本研究参考土地利用类型一级划分体系，将研究区划分为6种土地利用类型，选用连续的16个时相的土地利用结果，则可能出现的变化轨迹有616种，在应用稳定映射变化轨迹STD判定方法时，各土地空间单元演变过程中经历变化次数、多样性和相似性3个指标的计算更加复杂。由稳定映射STD指标的原理可知，长时序土地利用变化轨迹中，研究时段内某土地空间单元在相邻时间段发生变化的次数决定着该土地空间单元所经历的相同和不同土地利用类型的数量，即其所经历演变过程的变化次数T决定着该条变化轨迹上土地利用类型的相似性和多样性。而随着研究时段内时相数的增多，其变化轨迹中所经历演变过程的变化次数也越来越多，则土地利用类型就具有更高的相似性和多样性。本研究发现，时相数量的变化与变化次数、多样性、相似性3个指标的计算具有一定的对应关系，由此推导出在土地利用类型为6的情况下，时相数与STD指标的计算关系式(表1)，并通过随机抽样验证了其正确性。表1中，Ceil为向右取整函数。

表1 长时序土地利用变化轨迹的STD指标Tab.1 STD index for long time series land use change track

将本研究中n=16代入表1的范围中，可以计算出研究区土地利用变化过程中，经历不同变化次数的土地利用变化轨迹所包含的地类多样性和相似性，从而得到研究区16个时相的土地利用稳定映射STD指标(表2)，为下一步土地利用变化轨迹类型划分提供数据基础。

表2 研究区16个时相的土地利用变化STD指标Tab.2 Changes in soil utilization STD index at 16 in study area

2.3.2改进的土地利用变化轨迹分析方法

结合丹江流域(河南段)实际土地利用变化情况，本研究提出了一种基于长时序遥感的土地利用变化轨迹的STD判定方法(表3)。将一级变化轨迹类型划分为稳定型、渐变型、非连续渐变型、循环型和波动型5类。由于研究时序较长，研究区稳定土地空间单元会存在因土地利用规划占地或分类误差等影响而在某一年产生突变，但却长期保持该种土地利用类型的情况，因此本研究将稳定型土地利用变化轨迹划分为持续稳定型和长期保持型两类二级过程；同时，针对渐变型土地利用变化轨迹，又根据变化发生时期的早晚将其划分为早期变化后期稳定、中期变化和早期稳定后期变化3类二级过程；其次，根据研究区实际土地演变情况和STD指数，将循环型土地利用类型又划分为AB类循环和ABC类循环。最后，根据土地利用类型，将持续稳定型和长期保持型分别划分为持续稳定林地、持续稳定草地、持续稳定水域、持续稳定耕地、持续稳定裸地、持续稳定建设用地和长期保持林地、长期保持草地、长期保持水域、长期保持耕地、长期保持裸地、长期保持建设用地6类三级过程。

2.4 土地利用变化轨迹类型分析方法

传统的稳定映射变化轨迹分析方法是根据研究时段内区域土地利用变化的轨迹类型，分析其不同轨迹类型的面积变化和分布特征。近年来，国内外学者在分析其轨迹类型的基础上，结合定量分析模型和景观指数进行了更深层次的土地利用变化轨迹的演变特征分析。定量分析模型可以反映研究时段内区域的单一土地利用类型和综合土地利用类型的变化量、变化幅度和变化状态等特征，景观指数可以计算景观单元时空变化并反映景观格局结构等空间分布特征，其中景观水平上的景观格局分析能够反映出整个研究区的景观格局的演变规律[24-26]。本文将研究区土地利用变化轨迹划分为5种轨迹类型，结合其土地利用轨迹类型特征，选用单一土地利用动态度和综合动态度两个定量模型以及景观水平上的面积(CA)、面积占比(PLAND)、斑块密度(PD)、最大斑块指数(LPI)、散布与并列指数(IJI)和聚集度(AI)6种常用的景观指数，将变化轨迹相同的斑块看作同一种景观类型，利用ENVI软件和Fragstats软件，对研究区土地利用变化轨迹类型从动态变化总体特征和时空变化规律两方面进行系统分析。

表3 研究区土地利用变化轨迹的STD判定方法Tab.3 STD determining method of change track of soil use in study area

3 结果与分析

3.1 土地利用/覆被动态变化总体特征

在2002—2017年间，研究区各土地利用类型的组成结构和面积均发生了显著变化。从土地利用类型的组成结构来看，研究区地势从西北到东南依次为山地、丘陵、垄岗和平原，研究区的土地利用状况主要以林地和耕地为主；从土地利用类型的面积变化来看(表4)，随着研究区可持续发展进程加快以及南水北调中线工程的实施，2002—2017年间，研究区土地利用变化总体趋势主要以耕地向建设用地转化以及水域面积不断增加为主。具体表现为：①耕地面积在2004—2011年间持续减少，但2004年耕地面积较2003年增长约453.40 km2，主要原因为1999—2003年间，城镇、市区的迅速扩张和乡镇企业的突出发展占用了大量耕地，国家在该阶段相继出台了一系列保护耕地的法律法规，通过土地开发、复垦、整理等措施增加耕地面积[27]。2004年后，耕地面积逐年稳定减少。②林草地面积变化趋势较为复杂，2002—2004年呈现下降趋势，但在2006年后林草地面积有小幅度上升，主要因为该阶段“退耕还林还草”政策对其有较大的影响。 ③建设用地面积稳步增加，由2002年的106.35 km2增加至2017年的416.68 km2，共增加了310.33 km2。④水域面积总体呈上升趋势，原因在于2003年底南水北调中线工程开始实施，为使调水工程顺利进行，2005年开始对丹江口大坝进行加高加固处理，至2014年库区面积有了较为明显的增加。在2010—2012年期间，水域面积有所减少，与该时期丹江上游季节性降水和丹江口的人工调蓄有一定关系[28-29]。

表4 研究区土地利用类型面积统计Tab.4 Statistics of land use type area in study area km2

对比连续16年间研究区各土地利用类型的年变化率和综合土地利用动态度(图3)可以看出：①2002—2008年各土地利用类型的年变化率较大，这期间内研究区的土地利用类型之间的转变尤为频繁。②建设用地在研究时段内基本保持正向增长趋势，最大年增长率为51.00%，而耕地基本保持负向减少。③年变化率的计算受前一年度土地利用类型面积的影响，因此可以看到2005—2006年草地的年增长率为127.00%，主要原因为2003—2005年出台的一系列“复垦耕地”政策导致草地面积减少。④综合土地利用空间动态度曲线定量描述了研究区土地利用类型的变化幅度和活跃程度。由图3可以看出，研究初期丹江流域(河南段)各土地利用类型转换频繁，其中2008年土地利用变化最为剧烈；2009—2015年研究区综合土地利用动态度均在8%以下，其土地利用类型变化幅度较小；2016年后综合土地利用动态度有明显增加，表明这两年研究区土地利用类型变化较为活跃。

图3 研究区2002—2017年土地利用/覆被定量分析折线图Fig.3 Quantitative analysis of land use/covers in study area in 2002—2017

3.2 土地利用/覆被变化轨迹类型空间分布

基于改进稳定映射法的变化轨迹分析方法，按照土地利用变化轨迹动态过程将研究区变化轨迹划分为稳定型、渐变型、非连续渐变型、循环型和波动型5个一级类，并根据研究区实际土地利用情况划分二级和三级轨迹类型，然后结合6个常用景观指数，从空间构型上对研究区土地利用变化轨迹作深入分析(表5)。其构成特征如下：

(1)稳定型

此类轨迹类型表示自研究初期至研究末期，研究区土地利用类型从未发生过变化，或在研究期内长期保持某种地类的图谱单元。本研究中稳定型面积占比为61.80%。此种轨迹类型可以细分为持续稳定型和长期保持型。

持续稳定过程表示在连续16年间一直保持同一种土地利用类型的演变过程。结合图4a可以看出，在该种轨迹类型中，持续稳定林地面积最大，占研究区总面积的47.85%，主要分布于研究区北部和西北部的山地丘陵地带，其次为持续稳定耕地和持续稳定水域，面积占比分别为4.85%和2.03%，分别位于研究区东南部的平原地带和南部的丹江口库区。同时，分析最大斑块指数(LPI)和聚集度(AI)可知，持续稳定林地的最大斑块景观占比最高，为14.01%，聚集度为94.02%，说明持续稳定林地的景观破碎度较低，斑块聚集程度较高，林地的土地利用状况较为稳定；草地的最大斑块占比最小，趋近于0，聚集度为32.45%，说明该种地类景观破碎度高，主要以细碎斑块为主，土地利用变化活跃。裸地为过渡性土地利用类型，其面积占比很小，持续稳定裸地和长期保持裸地类型占比趋近于0，因此不对其进行景观指数分析。

长期保持过程是指某土地空间单元在研究期内土地利用类型原则上未发生变化，但由于研究时相较多，可能因为土地利用规划占地或分类结果存在误差等原因导致某一时刻地类发生一次变化，但仍长期保持该种地类的演变过程。结合图4b，此类轨迹斑块面积占研究区总面积的6.47%，主要集中在持续稳定型轨迹斑块附近，并在持续稳定型斑块周围离散分布。因此，长期保持型轨迹的AI值和整体水平AI值都偏低。

(2)渐变型

渐变型土地利用轨迹变化过程表示在研究区内土地利用类型只发生一次变化，即从一类转换为另一类的过程。该类型图谱单元变化次数为1，多样性为2，占研究区总面积的7.32%。根据变化时期的早晚，又将2002—2007年期间发生的变化归为早期变化后期稳定，2008—2011年期间发生的变化归为中期变化，将2012—2017年期间发生的变化归为早期稳定后期变化。结合图4c和表6可以看出，早期变化后期稳定的面积占比为3.86%，主要分布在研究区中部、北部和西南部，转化为林地、草地的面积较多；中期变化的面积占比为0.45%，主要转化为草地和建设用地；早期稳定后期变化的面积占比为3.01%，主要分布在研究区东南部平原地带和丹江流域支流的流经区域，对应的土地利用流向分别为耕地向草地、建设用地和耕地向水域的转化。

表5 研究区土地利用变化轨迹景观指数统计Tab.5 Land-use change track landscape index statistics

图4 研究区土地利用变化轨迹的空间分布Fig.4 Spatial distribution maps of land-use change trajectories

表6 渐变型土地利用类型转化面积统计Tab.6 Land use type conversion area statistics km2

(3)非连续渐变型

非连续渐变型是表示2种或3种土地利用类型的周期性转化，但趋势不明显的土地利用轨迹过程，占研究区总面积的8.86%。本研究提取出面积占比超过0.30%的变化轨迹如图4d所示，其中非连续耕地林地渐变和非连续耕地草地渐变共占此类型面积的87.33%。结合表5中非连续渐变型AI值为18.66%，说明该种轨迹类型由较多细碎小斑块组成，各个地类之间转换频繁。

(4)循环型

循环型轨迹变化过程是描述2种或3种土地利用类型间的相互转化，由于研究时序较长，时相较多，结合研究区实际土地利用状况，将该轨迹类型细分为两种地类间的循环转化和3种地类间的循环转化。经统计发现，耕地草地循环、耕地水域循环和林地草地循环的面积在该类循环中面积占比为83.50%，因此主要针对这3种循环转化进行分析(图4e)。其中耕地草地循环面积占比为3.79%，从空间分布来看，主要分布在稳定型耕地和稳定型草地的交错区域；耕地水域循环主要分布在淅川县丹江口库区支流的径流区域，面积占比为0.49%；林地草地循环面积占比为2.06%，受季节和区域影响，主要分布在研究区山地的山脊、山谷以及丘陵和平原相接的区域。

(5)波动型

波动型土地利用轨迹变化过程具有更高的多样性，表达不同地类间的频繁转化，占研究区总面积的14.43%。同时，因为波动型轨迹类型较多，本研究中仅将面积占比大于0.03%的轨迹类型选取出来。在此类轨迹变化过程中，图谱单元多样性为3的轨迹面积占该类总面积的64.20%，即在波动型的轨迹变化过程中，主要是3种地类间的频繁转化。其中，涉及到耕地和草地的变化轨迹居多，说明在研究区连续16年的土地利用变化中，耕地草地的变化是较为活跃的。

综上所述，在2002—2017年期间，研究区土地利用流向主要是耕地转化为建设用地；林草转化和耕草转化较为频繁；同时，水域面积明显增大，与南水北调中线工程的实施和丹江口库区大坝的增高紧密相关。从研究区土地利用变化总体特征和时空变化规律可以看出，16年来土地利用变化总体趋势主要以耕地向建设用地转化和水域面积不断增加为主，其中耕地面积由2002年的2 474.54 km2减少至2017年的1 562.60 km2，建设用地面积和水域面积分别由2002年的106.35 km2和218.60 km2增加至2017年的416.68 km2和400.31 km2，林草地在16年间面积变化较为复杂，呈现先减少后上升的趋势，主要受该阶段“退耕还林还草”政策的影响。综合分析16年间土地利用变化率可以得出，建设用地在研究时段内基本保持正向增长，这与16年间社会经济的快速发展和城镇、市区的迅速扩张有密切关系，同时2008年研究区综合土地利用动态度呈增加趋势，表明该年研究区土地利用类型的变化最为剧烈。研究区2002—2017年间稳定型轨迹面积占比最大，主要分布在研究区北部山区、丘陵地带林地、丹江口库区和东南部耕地等区域，达到61.80%，说明研究区总体土地利用状况稳定，符合研究区土地利用实际情况。而地类转化频繁的区域主要位于丘陵与平原的过渡地带，林地、耕地和草地出现了显著的相互转化现象。

4 讨论

提取土地利用/覆被变化轨迹，以研究土地利用类型的转移和流向，揭示土地利用变化规律，从而分析时间序列中土地利用/覆被的时空演变特征[30-31]。但随着研究时间序列数据的增多，土地利用/覆被变化轨迹更加繁杂，轨迹提取更加困难。本研究以丹江流域(河南段)为研究区，基于长时序遥感数据，提取土地利用变化图谱信息，并以改进的稳定映射变化轨迹分析方法将变化轨迹重分类，结合土地利用动态度等定量模型和景观格局指数，揭示研究区2002—2017年间土地利用类型的变化轨迹和时空演变规律，研究结果为土地资源可持续利用和水源地生态环境保护提供科学依据，同时也为南水北调工程后期进一步扩大引汉规模，并预计在2030年完成年均调水量达130亿m3的目标提供数据参考与决策支持。