1995年，“国际地圈生物圈计划”和“全球环境变化人文计划”两大国际科学计划组织共同拟定发表了《土地利用/覆被变化科学研究计划》，标志着土地利用/覆被变化(land use/cover change，LUCC)问题全球研究序幕的正式拉开^[1-3]。20多年来，LUCC问题始终是全球变化与可持续发展领域的关注热点和研究核心。不同学科的众多专家学者依据不同研究范式^[4]，围绕LUCC何地发生、何时发生、为何发生、如何发生等基本问题积极开展研究工作，在时空模拟、驱动机制、影响效应等方面提出了多种模型与方法，而转移矩阵则是其中最为常用的基础性模型之一^[5-9]。

目前，LUCC转移矩阵的计算生成大多采用遥感影像栅格数据，栅格化计算方法一般由数据预处理、地类解译分类、像元叠置运算和统计生成转移矩阵4个步骤组成^[10-11]。该方法原理简单、容易理解，具有支持叠加及统计分析的成熟模块工具较多、便于实施等优点，被广泛使用。但也存在如下三点不足：

1) 遥感影像一般仅包含土地利用/覆被要素的“地类”属性信息，所以只能用来计算生成地类转移矩阵，不能描述反映土地利用/覆被要素其他属性的变化转移情况。

2) 由于目前缺乏真正实用的全要素自动分类方法，遥感影像土地利用/覆被类型的识别与解译仍普遍需要大量的人工参与，周期长、成本高。另外，由于缺乏系统的标准规范指导和严格的质量控制措施，解译结果主观性强、准确度低。

3) 该方法以虚设的栅格像元作为统计单元，要通过合并处理方能将其转化为土地利用/覆被要素，受像元大小的影响转化结果与真实要素之间往往存在较大误差，进一步降低了计算结果的准确性。

针对上述方法的不足以及中国土地利用/覆被矢量数据日益增多的实际情况，本文首先提出了广义转移矩阵的概念，在描述挖掘及其内涵的基础上，设计了基于矢量数据的自动计算生成方法，并采用集成二次开发方式对相应工具软件进行了程序编码实现。

1 LUCC广义转移矩阵

LUCC转移矩阵来源于系统分析中对系统状态与状态之间转移的定量描述。目前，转移矩阵主要通过土地利用/覆被要素的位置、形态、面积等信息来概括描述其“地类”属性的变化情况，不仅能够直接反映研究区域内“地类”属性变化的组成结构、转移方向与转移面积，而且蕴含着变化率、变化速度、变化趋势等多种潜在信息。而事实上，除了位置、形态、面积、地类这些基本属性外，土地利用/覆盖要素还具有坡度、坡度等级、所有权人、权属性质、使用权人、使用类型等其他属性，理应利用转移矩阵对这些属性的变化情况进行分析，以更加全面地反映挖掘土地利用/覆盖要素变化性质和规律。为此，本文提出了广义转移矩阵(general transfer matrix，GTM)的概念。

1.1 形式化定义

所谓LUCC广义转移矩阵就是能够反映描述土地利用/覆被要素各种属性变化的转移矩阵。设F为选定的一个土地利用/覆盖要素的属性，在研究期初和研究期末其取值值域分别为{f₁，f₂…f_m}和{f₁，f₂…f_n}，那么该属性F在研究期内的广义转移矩阵可用表 1表示。

如果研究期初和研究期末F的值域相同(即m = n)，则GTM为一个方阵。表 1中，a_ii(i=1，2…m)表示研究期内属性值f_i保持不变的对应土地面积；a_ij(i≠j，j=1，2…n)表示研究期初属性值f_i向研究期末属性值f_j变化转移的对应土地面积；$A_{i+}=\sum\limits_{j=1}^{n} a_{i j} $表示研究期初属性值f_i对应的土地总面积；$ A_{+j}=\sum\limits_{i=1}^{m} a_{i j}$表示研究期末属性f_j对应的土地总面积；$ A_{{总}}=\sum\limits_{i=1}^{m} A_{i+}=\sum\limits_{j=1}^{n} A_{+j}$表示研究区域内土地的总面积。由于土地利用/覆盖要素的固位连续及不可消失性，研究期初和研究期末研究区域的土地利用总面积保持不变；$R_{i} =\sum\limits_{j=1, j \neq i}^{n} a_{i j}=A_{i+}-a_{i i}$表示研究期末属性值f_i变为其他值的对应土地面积；$ N_{j}=\sum\limits_{i=1, i \neq j}^{m} A_{i j}=A_{+j}-A_{j j}$表示研究期末其他属性值变为f_j的对应土地面积；$ R_{总}=\sum\limits_{i=1}^{m} R_{i}$表示研究期内属性F所有取值发生转出变化的对应土地总面积；$N_{总}= \sum\limits_{j=1}^{n} N_{j}$表示研究期内属性F所有取值发生转入变化的对应土地总面积。R_总和N_总表示的都是转移矩阵中前m行n列元素之和减去主对角线各元素之和，因此R_总 = N_总，表示整个研究区域内属性F的转出量和转入量相等。

1.2 其他隐含信息挖掘

除上述基本信息外，对广义转移矩阵做进一步的分析挖掘，还可以定义如下量化度量指标来更加全面地描述和揭示土地利用/覆被变化。

1) 净变化量(net change，NC)：是研究期初和研究期末属性值f_i对应土地面积之差的绝对值，也是属性值f_i转出变化与转入变化对应土地面积之差的绝对值，计算公式为：

$ \begin{array}{l} {\rm{N}}{{\rm{C}}_i} = \left| {{A_{i + }} - {A_{ + i}}} \right| = \\ \left| {\left( {{A_{i + }} - {a_{ii}}} \right) - \left( {{A_{ + i}} - {a_{ii}}} \right)} \right| = \\ |{R_i} - {N_i}| = \\ \max \left( {{R_i}, {N_i}} \right) - \min \left( {{R_i}, {N_i}} \right) \end{array} $

(1)

由于土地利用/覆盖空间区位的固定性与独特性，仅靠净变化量并不能真实反映土地利用/覆被要素的动态演变过程，往往低估其变化程度^[12]。例如，在一个研究期内，一块远郊区工业用地被整理为耕地，而另一块近郊区同等面积的耕地则被新开发为工业区。此时土地利用/覆被属性(耕地)数量没有发生变化，即净变化量为0，可在空间位置上却发生了实际变化。

2) 交换变化量(swap change，SW)：同时考虑了属性值f_i在一定区域转为其他值，而在另外区域又有其他值转为f_i的情况，它通过定量与定位相结合的方式弥补了净变化量的不足，计算公式为：

$ \mathrm{SC}_{i}=2 \min \left(R_{i}, N_{i}\right) $

(2)

3) 总变化量(total change，TC)：是净变化量与交换变化量之和，也是属性值f_i发生转入与转出变化的对应土地面积之和，即：

$ \begin{array}{l}{\mathrm{TC}_{i}=\mathrm{NC}_{i}+\mathrm{SC}_{i}=} \\ {\quad \max \left(R_{i}, N_{i}\right)+\min \left(R_{i}, N_{i}\right)=} \\ {\quad R_{i}+N_{i}}\end{array} $

(3)

4) 数量动态度(quantitative dynamic degree，QDD)：通过研究期内属性值f_i对应土地面积数据量的变化率来表示，其计算公式为：

$ \begin{array}{c}{\mathrm{QDD}_{i}=\frac{\left(A_{+i}-A_{i+}\right)}{A_{i+}} \times \frac{1}{T} \times 100 \%=} \\ {\frac{\left(N_{i}-R_{i}\right)}{A_{i+}} \times \frac{1}{T} \times 100 \%}\end{array} $

(4)

式中，T = t₂-t₁为研究时段长。当T设定为年时，QDD_i为研究期内属性值f_i的年变化率。该指标简明扼要，但不能反映土地利用/覆被动态变化的空间过程和相关属性，不能识别土地变化的“热点”或“敏感”区域^[12]。

5) 转出动态度(rolling-out dynamic degree，RoDD)：通过研究期内属性值f_i发生转出变化的对应土地面积变化率来表示，其计算公式如式(5)所示。该指标仅考虑了属性转出这一单向变化过程，容易低估忽略转出慢、但转入(增长)快的土地利用/覆被属性(如城市建设用地)。

$ \mathrm{RoDD}_{i}=\frac{R_{i}}{A_{i+}} \times \frac{1}{T} \times 100 \% $

(5)

6) 转入动态度(rolling-in dynamic degree，RiDD)：通过研究期内属性值f_i发生转入变化的对应土地面积变化率来表示，其计算公式如式(6)所示。显然，RiDD = QDD + RoDD。

$\operatorname{RiDD}_{i}=\frac{N_{i}}{A_{i+}} \times \frac{1}{T} \times 100 \% $

(6)

7) 双向动态度(bidirectional dynamic degree，BDD)：通过研究期内属性值f_i发生转出和转入变化的对应土地面积变化率来表示，其计算公式如式(7)所示。显然，BDD_i = RoDD_i + RiDD_i。

$ \mathrm{BDD}_{i}=\frac{R_{i}+N_{i}}{A_{i+}} \times \frac{1}{T} \times 100 \% $

(7)

8) 综合动态度(total dynamic degree，TDD)：通过研究区域内属性F发生转出或转入变化对应土地总面积的变化率来表示，其计算公式为：

$ \begin{array}{l} {\rm{TDD}} = \frac{{{R_总}}}{{{A_总}}} \times \frac{1}{T} \times 100\% = \\ \frac{{{N_总}}}{{{A_总}}} \times \frac{1}{T} \times 100\% \end{array} $

(8)

2 LUCC广义转移矩阵自动生成

随着全国土地利用与变更调查、全国地理国情普查与监测等项目的组织实施与常态化推进，目前我国各级土地资源管理与测绘生产服务部门已经或即将积累拥有大量不同时期的土地利用/覆被矢量数据。这些数据通常以快照序列的形式加以存储管理，每一快照都准确详细地记录了某一时刻土地利用/覆被要素的空间位置、几何图形以及各种属性信息，为广义转移矩阵的计算生成提供了充分、精准的数据基础。

然而，目前却缺乏直接利用这些数据快速挖掘生成广义转移矩阵的有效方法，所涉及的匹配、对比、识别、计算、汇总等操作大多需要在常用GIS软件工具的支持下由人工手动完成，人机交互频繁、效率低、易出错，难以用于大数据量的计算。为此，本文提出了一种基于矢量数据的广义转移矩阵自动生成方法。

2.1 数据预处理

为保证计算过程及结果的可审、准确、高效，首先应该对所获得的原始矢量数据进行备份，然后再对数据副本进行预处理，以供计算生成满足需要的转移矩阵。

对于多时态土地利用/覆被矢量数据而言，所涉及的预处理主要包括坐标系统转换、数据格式转换、模式结构同化、类别编码转换、升级、数据再分类等内容^[13]。在本研究中，还应包括数据集空间及属性索引建立与更新、融合组成同一要素的多图元等内容，这些对提高后续处理效率至关重要。

2.2 匹配关系建立

根据§1.1的形式化定义可知，广义转移矩阵的计算生成是以单个土地利用/覆被要素的属性变化为前提的。只有明确了单个要素的属性变化情况，才能对所有要素属性变化情况进行汇总统计, 以生成具体的转移矩阵，而单个要素变化情况的识别则需要通过要素匹配操作，来建立研究期初与研究期末土地利用/覆被矢量数据集中同一要素之间的对应匹配关系。

为了克服常用人工判别匹配的不足，人们已研究提出了多种自动要素匹配方法。这些方法的基本策略是计算比较待匹配要素在位置、图形、拓扑、属性等上的相近性，如果相近性达到某种程度，则认为它们是匹配对应要素^[14-15]。根据这一基本策略并考虑到属性匹配执行效率相对较高的特点，通过如下5个步骤来匹配发现不同时期土地利用/覆被矢量数据集之间的对应匹配要素。

1) 从研究期初土地利用/覆被矢量数据集中提取一个要素记录r_i，其系统唯一标识属性值为sid，如果存在用户定义唯一标识属性且取值为uid，则执行第2)步；否则，执行第3)步。

2) 从研究期末土地利用/覆被矢量数据集中查询检索用户定义唯一标识属性与r_i.uid相同的要素记录，如果查到要素记录r_j，则认为r_i与r_j对应，并计算二者间的重叠面积，连同其系统唯一标识值r_i.sid、r_j.sid记录在“要素匹配变化表”中，返回1)步重新提取下一个要素记录；否则，执行第3)步。

3) 从研究期末土地利用/覆被矢量数据集中查询检索位置形状与r_i完全重合的要素记录，如果查到要素记录r_j，则认为r_i与r_j对应，并计算二者间的重叠面积(实质上为r_i或r_j的面积)，连同其系统唯一标识值r_i.sid、r_j.sid记录在“要素匹配变化表”中，返回1)步重新提取下一个要素记录；否则，执行第4)步。

4) 从研究期末土地利用/覆被矢量数据集中查询检索与r_i重叠面积大于0的要素记录，如果查到要素记录集R={r_j}，|E|≥1，则认为r_i与R中的所有要素对应，并将r_i.sid与R中所有记录的系统唯一标识值以及彼此间的重叠面积值记录在“要素匹配变化表”中，返回1)步重新提取下一个要素记录。

5) 如果研究期初土地利用/覆被矢量数据集中的要素已被全部提取一次，则终止匹配过程。

2.3 要素变化识别

对矢量数据集而言，通过双向匹配操作一般可以获得基数为1:0、1:1、1:m、m:1、m:n、0:1的6种对应关系，依据这些对应关系可以识别发现未变化以及消失、位移、缩小、扩张、分解、合并、聚合、新增等多种类型的要素变化^[16-17]。在本文研究中，由于面状土地利用/覆被矢量数据的空间固定连续无缝性以及匹配操作的单向性(期初要素向期末要素)，因此，要素匹配对应关系实际上仅有1:1和1:m两种类型。

根据上述两种匹配关系，通过进一步的分析比较可得到相应要素的几何图形变化以及属性变化。对于几何图形的变化，本文制定了如下具体识别规则：

1) 对于1:1匹配，如果两个要素的面积相同，则认为研究期初要素没发生几何图形变化; 如果不相同，则认为研究期初要素发生了扩张变化。

2) 对于1:m匹配，如果研究期初要素面积与对应研究期末要素总面积相同，则认为研究期初要素发生了分解变化; 如果不相同，则认为研究期初要素可能发生侵蚀、分解、扩张等综合变化。

对于属性的变化，其识别规则如下：

1) 对于1:1匹配，如果两个要素的字符型属性取值不相同，则认为研究初期要素发生了属性变化。

2) 对于1:1匹配，如果两个要素的数值型属性取值之差的绝对值大于预设阈值，则认为研究期初要素发生属性变化。该规则可以降低数值型属性取值精度、误差对变化识别结果的影响。

3) 对于1:m匹配，将其分解为m个1:1形式的匹配，然后根据上述两个规则来判断识别其属性的变化。

通过上述规则，对所有匹配要素进行对比分析，并将变化结果记录在要素匹配变化表中。

2.4 矩阵设置与生成

对所得到的单个要素变化信息做进一步的处理计算汇总，可灵活准确地生成土地利用/覆被要

素在某一属性上的总体变化转移矩阵，具体步骤如下：

1) 选择设置拟研究分析的土地利用/覆被要素的属性，统计该属性在研究期初和研究期末的唯一取值，根据取值个数设置转移矩阵的行数与列数，并将属性取值准确记录在矩阵相应位置。

2) 以“期初属性值”+“→”+“期末属性值”的形式依次组合所有属性值，在要素匹配变化表的相应属性列中查询这些属性值组合。如果查询结果非空，则对所有返回记录包含的重叠面积信息进行累加求和，所得结果即为该属性组合所示变化情况对应的土地总面积，将其记录在转移矩阵对应位置；如果查询结果为空，则表示实际中没有这种变化情况，对应的土地总面积为0，将其记录在转移矩阵对应位置。

3) 根据第2)步计算结果及§2.1中所列相关公式，进一步计算完善转移矩阵中“总计”、“转入”、“转出”行与列中的相应信息，最终生成具体完整的转移矩阵。

3 试验验证与分析

根据上述方法与策略，采用Visual C# 2010与ArcObjects 10集成二次开发方式设计实现了基于矢量数据的土地利用/覆被变化广义转移矩阵计算生成软件工具(简称General TM)，其主界面如图 1所示。在该工具支持下，分别采用武汉市汉南区、琼海市大路镇、文昌市潭牛镇、长沙县干杉镇、长沙县黄兴镇共5个县区、乡镇的第二次土地利用现状及变更调查地类图斑矢量数据，在配置为AMD 2.90 GHz CPU、8.00 GB内存、64位Windows操作系统的个人PC上，对所提方法进行了试验验证。

3.1 总体性能分析

表 2列出了所选5对数据集的基本特征以及在不同处理计算阶段所用的时间。

在试验过程中，没有对后4对数据集进行同要素多图元融合处理，而是直接将面积相对较小的图元视为要素，以模拟和验证较多要素量的运行环境要素和运行效率。要素匹配均以几何重叠为依据；均选择一级地类(由二级地类派生而来)、二级地类、权属性质、权属单位4个属性作为变化检测内容；矩阵计算生成4个属性对应的转移矩阵。

通过与人工手动计算进行对比，本文方法执行效率平均提高10倍以上。用户可以通过表格、图形等形式灵活便捷地对中间及最终处理结果进行检查核对，加上计算结果间的自动交叉验证，能够有效避免错算、漏算情况的发生，准确度可达98%以上。

3.2 部分结果分析

限于篇幅，这里仅以武汉市汉南区试验为例，对转移矩阵挖掘提取过程及相关结果做进一步的分析与说明。

1) 要素匹配结果。以几何重叠作为匹配条件，共发现1:1匹配824对，1:n匹配16对(其中，1:2匹配12对，1:3匹配2对，1:4匹配2对)，将1:n匹配转为1:1匹配后，共有1:1匹配862对。

2) 变化检测结果。比较匹配发现的862对要素在4个所选属性上的取值，共发现一级地类变化23对；二级地类变化23对；权属性质变化3对；权属单位变化0对。要素匹配及变化检测结果记录形式如图 1所示(面积单位：m²)。

3) 转移矩阵计算结果。依次对要素匹配变化表中记录的匹配要素属性值和重叠面积进行汇总统计，可计算生成4个转移矩阵，这里仅给出一级地类转移矩阵和权属性质转移矩阵，分别如表 3、表 4所示。

4) 其他度量指标计算结果。根据转移矩阵及§2.2讨论的相关公式，可挖掘计算出不同属性取值的其他变化特征。

4 结语

本文定义了LUCC广义转移矩阵的概念，设计实现了利用两期快照式矢量数据自动挖掘生成广义转移矩阵的方法与软件工具，在县乡尺度上验证了其可行性与有效性，为高效充分利用我国日益增多的多时态土地利用/覆被矢量数据，全方位多属性分析土地利用/覆被变化特征和规律提供了技术支撑和数据依据。

在后续研究中，计划从3个方面对现有工作进行扩充完善：①结合其他数据资料，对已获得试验结果做进一步的分析，更为深入和准确地认识研究区土地利用/覆被变化的动因、机制和过程；② 扩展所提方法对不同结构组织形式数据源的支持，着重研究面向增量式时态矢量数据的转移矩阵计算方法；③采用并行计算技术进一步优化现有方法的运行效率，以适应市、省尺度上大数据量的挖掘计算处理。