基于ArcGIS ModelBuilder的LUCC转移矩阵自动计算模型研究

王育红

(江苏师范大学 地理测绘与城乡规划学院，江苏 徐州 221116)

自1995年，“国际地圈生物圈计划(IGBP)”和“全球环境变化人文计划(HDP)”两大国际科学计划组织共同拟定发表《土地利用/覆被变化科学研究计划》以来，土地利用覆被/变化(Land Use/Cover Change，LUCC)一直是世界各国学术界持续关注和广泛研究的热点。人们从不同目的出发构造多种LUCC分析与评价模型[1]，如空间统计模型、系统动力学模型、元胞自动机模型、智能体模型以及综合模型等[2]。

转移矩阵作为其中最为核心的基础性模型，更是被众多专家学者所使用。例如，史培军等利用转移矩阵分析了深圳市1980—1994年间的土地利用变化空间过程[3]；段增强等基于转移矩阵分析挖掘了北京市海淀区1991—2001年间的土地利用变化信息[4]；谢金开等通过转移矩阵研究了乌鲁木齐河下游地区1989— 2007年间的土地利用/覆盖时空动态变化规律[5]；廖谌婳等凭借转移矩阵对中老缅泰交界地区1990—2010年间的土地利用变化特征进行了分析[6]；任东风等运用转移矩阵对于彰武县2008—2018年间的土地沙漠化问题进行了研究[7]。

从文献资料上看，现有的大多数研究都是直接给出研究区的转移矩阵结果，对其计算方法与过程往往缺乏深入详细的阐述。针对这一情况，文中在界定总结转移矩阵概念及其计算原理的基础上，描述ArcGIS平台下人工手动计算转移矩阵的过程及工具，并采用ModelBuilder技术通过集成现有工具设计实现一个转移矩阵自动批量计算工具模型。

1 转移矩阵及其计算原理

LUCC转移矩阵来源于系统分析中对系统状态与状态之间转移的定量描述，可以从数量和结构上，综合描述某一地区在期初和期末限定的研究时期内，各种土地利用/覆被类型之间的总体变化转移情况。在实际应用中，转移矩阵有两种形式：即以变化转移面积绝对数量为统计指标的转移面积矩阵(通常称为转移矩阵)和以变化转移面积相对数量为统计指标的转移概率矩阵。

1.1 定义与性质

假设某一研究区的土地利用/覆盖类型(简称地类)在研究期初(T0)为m种，其编码值分别为1，2，3，…，m；在研究期末(T1)为n种，其编码值分别为1，2，…，n。如果n=m，表示研究时期内地类种类没有变化；如果n>m，表示研究期末地类种类增加；如果n

表1 地类转移面积矩阵

表2 地类转移概率矩阵

1.1.1 地类转移面积矩阵及其基本性质

1)aii(i≤min(m,n))表示研究期间地类i没有发生变化转移的土地面积。

2)aij(i≠j,1≤i≤m,1≤j≤n)表示研究期初地类i在研究期末向地类j变化转移的土地面积。

3)每一行的面积值之和(记作Ai)表示研究期初地类i对应的土地总面积，即：

(1)

4)每一列的面积值之和(记作Ai)表示研究期末地类j对应的土地总面积，即：

(2)

5)由于在期初和期末研究区的土地总面积保持不变，即式(3)恒成立。

(3)

1.1.2 地类转移概率矩阵及其基本性质

2)每一行的转移面积比例值之和为1，即式(4)恒成立。

(4)

为全面反映地类变化转移数量和程度，建议将两种转移矩阵结合起来使用，避免误判或低估某些地类的变化转移情况。例如，地类变化转移绝对面积值相对较小而相对变化转移比例(概率)值却很大的情况。但在土地利用/覆盖空间数据集采用经纬度地理坐标系，不便于计算真实面积值的情况下，应使用转移概率矩阵。通过转移矩阵可以进一步挖掘获取地类转移变化的净变化量、交换变化量、总变化量及其持续性、随机性、系统性等特征[8-9]，转移概率矩阵还可用于地类未来构成及数量的估算与预测[10]。有关这方面的具体模型与方法，这里不再详述，请参考相应文献资料。

1.2 基本计算原理

由于根据转移面积矩阵可计算派生出转移概率矩阵，因此本节及以后内容将围绕转移面积矩阵(简称转移矩阵)讨论相关问题。目前，转移矩阵的计算主要利用遥感影像解译分类后所获得的土地利用/覆被栅格数据集，其核心思想是地类代码融合。为便于叙述与理解，文中采用具体的示意性数据对其进行解释。

假设研究区由8×8个像元组成，反映研究期初与期末土地利用/覆被状况的栅格数据集分别为R0和R1，两期都有5种地类，分别用代码1～5表示，前后两期的具体像元地类代码值分别如图1所示，那么根据R0和R1计算转移矩阵将由如下两大步骤组成：

1)采用地图代数按照式(5)对R0地类代码进行放大，并将放大处理结果与R1相加，形成反映地类变化转移情况的栅格数据，即:

R01=R0×s+R1.

(5)

其中，s为初地类代码放大倍数，其值取决于期末地类代码最大值的位数(b)，两者之间的关系为s=10b。图1中R1数据集相应的b为1，故放大倍数s为10，因此可得如图1所示的R01。

如果土地分类栅格数据集的像元地类代码信息以文本名称、标准编码等形式表示，计算前应按相应的映射规则将其转换成连续的整数型代码。

图1 基于栅格数据的转移面积矩阵计算核心思想示意图

2)对R01中不同编码值的像元个数进行统计，根据空间分辨率(栅格像元所代表的实地面积)计算不同编码对应的土地总面积，并对统计表格进行相应转换处理以生成转移矩阵。假设图1中像元的空间分辨率为10 m×10 m，则对R01进行统计,所得表格和最终生成的转移矩阵可分别为表3和表4的形式。

表3 栅格数据R01统计表

表4 转移矩阵计算结果

2 基于ArcGIS的计算方法与模型

ArcGIS是全球著名GIS技术与软件提供商—美国ESRI公司发布的全系列GIS软件平台的统称，主要包括桌面GIS、嵌入GIS、移动GIS以及服务端GIS等多种子系统。据ARC咨询公司2019年调查研究报告的数据显示，ArcGIS产品占全球GIS市场的45%以上[11]。初步估计，目前ArcGIS产品也至少占国内GIS市场的25%以上。基于ArcGIS在国内的应用广泛性，本节进一步讨论基于ArcGIS桌面系统的转移矩阵具体计算方法、步骤及工具。

2.1 手动计算方法

根据上述原理，使用ArcGIS桌面系统中的ArcMap模块即可通过人工手动方式实现转移矩阵的计算，具体计算过程主要包括如下步骤。

2.1.1 生成栅格数据集R01

打开位于“空间分析工具箱地图代数工具集”下的“栅格计算器”工具，在该工具中的计算表达式窗口中输入形如“R0”*s+“R1”的表达式，并指定输出栅格数据集(即R01)存储格式和位置，确认无误后单击“确定”生成R01。其中，s为代码放大系数。

2.1.2 添加和计算字段值

打开R01的属性表，在“ObjectID”、“Value”与“Count”现有字段的基础上，添加“期初代码”(整型)、“期末代码”(整型)和“面积”(浮点型)，并在相应字段依次选择右键菜单中的“字段计算器”，分别使用形如:Left([Value], Len([Value]) -b)、Right([Value],b)、[Count]*r的VB脚本语言计算表达式计算表中记录3个字段的相应值。其中，b为期末地类代码最大值位数；r为R01的空间分辨率(面积值)。

2.1.3 复制和透视属性表

打开位于“数据管理工具箱表工具集”下的“复制行”工具，通过该工具将处理后的R01属性表，复制备份到指定存储位置。打开位于“数据管理工具箱表工具集”下的“数据集透视表”工具，在该工具中选择刚复制的表作为输入表，并分别将输入字段、透视表字段和值字段三参数选择设置为“期初代码”“期末代码”和“面积”字段，最后指定输出表位置单击“确定”即完成转移矩阵的计算。

由于ArcGIS表格不支持数值型字段名及单元格的合并，所生成的转移矩阵还难以做到与上述的定义形式完全一致。如果需要可通过修改字段别名，做到除“期初(T0)地类”“期初(T1)地类”字样之外的非常一致。

通过上述描述可以看出，利用ArcMap手动计算转移矩阵仍是一项比较繁琐的工作，人机交互频繁，所需工具及参数设置多，稍有不慎将出现错误，进而影响计算效率及结果的准确性。为解决这一问题，作者进一步采用Modelbuilder技术对上述计算步骤与工具进行组合与封装，构建形成一个通用型自动批量计算工具模型。

2.2 自动计算模型

ModelBuilder是ArcGIS桌面系统中一个用来创建、编辑和管理模型的子程序，也可将其看作一种可视化编程语言。ModelBuilder可直接利用ArcGIS系统已有的各类工具，按一定的规则将其有机组合串联起来，最终形成一个流程化的模型程序，从而实现高效的分析计算[12-14]。

根据上述手动计算流程，可通过如下主要步骤在ModelBuilder模型编辑窗口中创建实现转换矩阵自动计算模型工具。

2.2.1 添加输入变量

为提高所建模型工具的通用性，添加两个栅格图层类型和一个整数类型的变量，其名称分别为“期初栅格数据集”“期末栅格数据集”和“期末地类代码最大值位数”。

2.2.2 添加、设置、连接现有处理工具

通过拖拽方式依次从ArcMap系统工具箱里将“计算值”(位于数据管理工具箱常规工具集下，用于计算放大倍数)、“栅格计算器”“添加字段”“计算字段”(与前一个工具一样都位于数据管理工具箱字段工具集下，都需要添加3个)，“获取栅格属性”(位于数据管理工具箱栅格工具集下，需要添加两个分别用于获取单个像元的宽度值和高度值)、“复制行”“数据透视表”等工具添加到模型编辑窗口中。

在添加工具之后，依次打开相应处理工具通过选择已有变量及相关工具输出结果或手动输入的方式设置处理工具所需参数，并将工具连接起来形成工作流。为提高模型易于阅读和理解的程度，还需进一步修改工具及其输出结果的名称。

2.2.3 设置输出变量及模型属性

为方便处理结果的管理、查看与对比，除“数据透视表”工具输出的最终处理结果(即转移矩阵)外，还将“栅格计算器”工具输出的中间处理结果(即研究期内地类变化转移栅格数据集)设置为输出变量。

在验证模型能无误运行后，在模型属性对话框的“参数”选项卡中调整变量参数显示顺序，将“期末地类代码最大值位数”的过滤器参数设置为“值列表”，并添加1、2、3等3个可选值。另外，在“常规”选项卡中，为名称、标签、描述等参数设置准确的信息，以方便工具的查找和使用。

最后所得到的模型工具内部工作流程及外部运行界面效果分别如图2与图3所示。

图2 模型内部工作流程

3 初步实验与效率对比

为了对比分析上述计算方法与工具的效率性能，笔者以前期相关研究积累的南京、郑州、徐州3地在不同时期的6幅土地利用分类栅格数据集为例进行初步实验。其中，南京市包括城市用地、非城市用地以及未定义用地3种地类，郑州、徐州两地包括水体、林地、农田、菜地、旱地、建设用地6种地类。笔者邀请3位本科生使用这些数据在各自的笔记本电脑上，分别采用手动和自动方式对相应时期内的转移矩阵进行10次计算。

图3 模型外部运行界面效果图

表5给出实验数据的基本特征及每人的平均计算用时。从表中可以看出，通过模型自动计算用时平均可比人工手动计算缩短10倍以上。另外，该工具只需要选择输入5个参数，避免常规手动计算方法用户全程参与、工具及参数设置多、人机交互频繁、易出错的不足。

表5 实验数据集基本特征及计算用时统计表

4 结束语

文中在界定总结转移矩阵概念、性质及其计算原理的基础上，详细阐述ArcGIS桌面系统下的人工手动计算方法与自动计算模型创建过程，为土地资源管理、规划、科研以及相关开发人员提供有益的理论参考与技术支持。随着全国土地利用调查及地理国情监测等国家重大战略工程的常态化实施，我国土地利用/覆被矢量正以前所未有的速度快速增长，亟需开展深层次的挖掘分析。针对这一需要，研究创建基于矢量数据的转移矩阵自动计算方法和模型工具将是下一步研究工作的重点。