两期土地利用图摆在面前一张是几年前的一张是最新的任务只有一个算出这片区域的土地利用动态图。这个需求在 ArcGIS 里几乎是最常见的一类活儿——说它简单流程确实不复杂说它坑多也确实能让人来回折腾一整天。我做土地利用变化分析这些年见过太多人卡在第一步两期数据的分类编码对不上叠加出来一片乱码算出来的动态度要么大得离谱要么小数点后全是零。先把这件事说清楚土地利用动态图这个词在实际项目里其实是个笼统叫法它往往同时包含三块内容——各类地类的单一土地利用动态度每年变化百分之几、全区的综合土地利用动态度整体活跃程度以及支撑这两者的土地利用转移矩阵和变化图斑分布图。前两个是数字后两个是图和表四样东西凑齐才算把动态这件事讲明白。这篇文章面向的是手上有两期或多期土地利用矢量数据、会用 ArcGIS 基本工具、但没系统做过变化分析的人。如果你正在做土地资源调查、生态评估、国土空间规划前期分析或者毕业论文里要用到土地利用动态度这个指标下面这套流程可以直接照着走。我会把每一步为什么这么做、参数怎么定、哪里最容易翻车都讲透公式和算例一并给出。1. 先把土地利用动态图这五个字拆开看1.1 动态度、变化图、转移矩阵三个东西三条路线很多人一上来就问动态度在哪个工具里算这个问题本身就问偏了。ArcGIS 里没有哪个按钮叫计算动态度动态度是一个统计指标它的输入是各地类的面积输出是一个百分比数字ArcGIS 负责的是把面积统计出来剩下的四则运算可能要在属性表或者 Excel 里完成。具体拆开看单一土地利用动态度针对某一个地类比如耕地算它在一段时间内平均每年变化了多少个百分点。输入是这一类地类的期初面积和期末面积。综合土地利用动态度针对整个研究区把所有地类之间的转移面积加总除以期初总面积衡量区域整体的土地利用活跃程度。土地利用转移矩阵一张 n×n 的表格行是期初地类列是期末地类交叉格子里是对应的面积对角线是没变的部分。变化图斑图把所有发生了地类改变的多边形挑出来做成专题图直观展示变化发生在哪里。这四样东西的关系是转移矩阵是底层数据动态度是从矩阵里提炼出的指标变化图斑图是空间表达。路线上有矢量叠加和栅格叠加两条矢量精度高、边界贴合实际栅格速度快、适合大范围快速评估。我的习惯是研究区不大比如一个县以内或者要求图斑边界精确走矢量跨市跨省、动辄几千万个图斑走栅格。提示动态度这个词在不同文献里的公式写法不完全一样尤其是综合动态度有的分母带系数 2有的不带。写报告前先确定你参考的那篇文献用的是哪种别中途换公式导致前后数据对不上。1.2 数据口径不统一后面全白干真正常见的翻车点不在计算而在数据准备。土地利用数据的来源五花八门可能是某次国土调查的成果可能是遥感解译的分类结果也可能是从别人那里拷来的历史数据。这些数据的分类体系、编码位数、坐标系、边界范围经常不一致。我遇到过最典型的一次两期数据的耕地编码一期是01另一期是0101看着都是耕地但字段一比对全是错位叠加之后原本没变化的耕地被识别成了变化图斑动态度算出来高得离谱。还有一次是两期数据的行政边界范围差了半条街相交之后边缘出现一条几百米宽的空白带所有地类的面积都少算了一块。所以在动手之前必须先做三件事确认两期数据用的是同一套分类体系。如果一期是国标二级类、一期是一级类先归并到统一的层级再做。归并的原则是向上归并二级类合成一级类不会丢信息反过来则不行。确认编码字段类型一致都是文本型或者都是数值型位数补齐。文本01和数值 1 在连接和分组统计时是两个东西。确认研究区边界统一用同一份行政边界去裁剪两期数据或者先用相交工具取两期的公共范围之后所有统计都基于这个公共范围。这三件事做完可能只花你二十分钟但能省下一整天的返工。2. 数据准备把两期数据磨成能算的状态2.1 地类编码对齐与融合去碎两期数据的分类对齐之后下一步是融合Dissolve。为什么要融合因为很多土地利用数据的图斑碎得离谱同一块连片的耕地可能被切成七八个多边形中间还夹杂着几十平方米的细碎图斑。直接拿去做叠加分析输出要素数量会呈指数级增长跑起来慢结果也难看。融合的操作路径是按地类编码字段做 Dissolve在统计字段里保留面积字段求和其实面积后面重新算也行。融合前建议先删掉面积小于某个阈值的碎图斑比如小于 0.1 公顷的图斑根据研究精度决定。这一步要慎重因为删多了会影响总面积我的做法是把碎图斑合并到相邻的最大图斑里而不是直接删除这样总面积不变。还有一种情况是同一个地类编码下存在重复的、相互重叠的图斑。这种数据在叠加时会产出大量零面积的碎屑必须先处理掉。处理办法是跑一遍拓扑检查里的不能重叠规则把重叠部分找出来修正。ArcGIS Pro 里还有个提速技巧如果两期数据都很大用Pairwise Intersect代替普通的 Intersect多核并行速度能快好几倍结果基本一致。这个工具在 Analysis Tools Pairwise Overlay 下面很多人不知道它的存在。2.2 拓扑检查与几何修复的实操套路土地利用数据要做叠加分析几何质量必须过关否则报错或者出怪结果。我固定会做这几项检查检查项工具/规则常见问题处理方式自相交拓扑规则不能自相交多边形边界自己打结用修复几何工具自动修重叠压盖拓扑规则不能重叠两块地类互相压盖手动编辑或用联合后重建空隙拓扑规则不能有空隙图斑之间有条缝补充缝隙多边形或融合尖锐角几何检查或专用插件夹角过小导致后续叠加报错删除节点或局部重建伪结点、悬挂点拓扑规则线要素未闭合编辑捕捉修正尖锐角这个东西平时不起眼但在做大批量叠加的时候特别要命。角度过小的节点会让叠加算法在某些位置产生异常细长的多边形面积算出来是个极小的正数或者干脆是零。市面上有一些专门检查尖锐角的插件如果不装插件也可以自己用字段计算器算相邻边的夹角或者干脆把夹角小于某个阈值比如 1 度的节点在编辑时手动拉开。注意修复几何Repair Geometry能自动处理自相交、空几何、环方向错误这类问题但它不会处理重叠和空隙那两类必须靠拓扑和人工编辑。还有一点容易被忽略地理坐标系下的数据不要直接做叠加分析。虽然工具能跑但容差单位是度精度完全不可控叠加结果可能出现肉眼可见的偏移。2.3 坐标系与面积口径为什么一定要投影坐标系这是我最想强调的一条。只要涉及面积计算必须用投影坐标系且最好是等面积投影。用地理坐标系经纬度算面积ArcGIS 虽然会给结果但那是在球面上按测地线方式算的单位是平方米没错可一旦数据范围跨了几个纬度带误差会明显放大。更麻烦的是很多人在叠加分析之后直接读 Shape_Area 字段而这个字段是在数据本身的坐标系下算的如果原数据是地理坐标系Shape_Area 的单位是平方度这个数没有任何物理意义。我的标准做法是根据研究区位置选一个合适的投影坐标系。中国范围内常用的是 CGCS2000 的高斯-克吕格投影按 3 度带或 6 度带分带如果研究区跨带就换成 Albers 等面积投影中央经线和双标准纬线按研究区范围设定。两期数据投影到同一个坐标系投影参数完全一致。面积字段统一用公顷因为土地利用动态度、转移矩阵的常规表达单位就是公顷。平方米转公顷除以 10000直接用 Calculate Geometry 时选 HECTARES 就行。选 3 度带还是 6 度带也有讲究。3 度带变形更小适合精度要求高的县级尺度研究6 度带带号少适合省级以上范围。如果你不确定用 Albers 等面积投影最保险因为动态度和面积直接挂钩等面积投影能保证面积不变形。投影做完建议先算一遍两期数据的各地类总面积各自核对一下看看是不是和研究区公布的统计面积对得上。这一步是体检如果面积对不上后面所有计算都不可信。对不上的原因通常是边界裁剪不一致、碎图斑没处理干净、或者分类编码有重复。3. 单一动态度与综合动态度公式、参数与算例3.1 单一土地利用动态度怎么算单一土地利用动态度的公式不复杂但每一个参数的含义要弄清楚$$K_i \frac{U_b - U_a}{U_a} \times \frac{1}{T} \times 100%$$其中$K_i$ 是第 i 类地类的单一土地利用动态度单位是百分比每年$U_a$ 是该地类在研究期初的面积单位公顷$U_b$ 是该地类在研究期末的面积单位公顷$T$ 是研究时段长度单位年。这个公式的物理含义很直白先算整个研究期内这一类地类变化了多少比例再除以年数得到平均每年的变化率。正数说明面积增加负数说明面积减少绝对值越大说明变化越剧烈。举个算例。某研究区 2010 年耕地面积 10240.5 公顷2020 年耕地面积 8632.7 公顷研究时段 $T 10$ 年。代入公式$$K_{耕地} \frac{8632.7 - 10240.5}{10240.5} \times \frac{1}{10} \times 100% \frac{-1607.8}{10240.5} \times 0.1 \times 100% \approx -1.57%$$也就是说这十年里耕地年均减少约 1.57%。这个数字放在报告里配上一句耕地年均减少 1.57%主要流向建设用地和林地说服力就出来了。关于 $T$ 的取值有个细节。如果期初数据是 2010 年末的期末是 2020 年末的那 $T$ 就是 10 年整。但如果期初是 2010 年 6 月、期末是 2020 年 3 月严格算的话 $T$ 应该是 9.75 年左右。实际项目里大多数人是按年份差直接取整数这个做法可以接受但如果时段很短比如两年T 取整带来的误差就不能忽略了这时候建议按实际月份折算。3.2 综合动态度分母到底取什么综合土地利用动态度的公式版本特别多这里给一个在中文文献里最常见、也最容易被审稿人接受的写法$$LC \frac{\sum_{i1}^{n} \Delta LU_{i-j}}{\sum_{i1}^{n} LU_i} \times \frac{1}{T} \times 100%$$参数含义$LC$ 是综合土地利用动态度$\Delta LU_{i-j}$ 是第 i 类地类转为其他地类的面积之和的绝对值注意这里取的是转出面积不包括自身不变的对角线部分$\sum LU_i$ 是研究期初所有地类的面积总和也就是研究区总面积$T$ 同样是时段长度。关于分母还有一种写法是分母乘 2理由是变化面积被转出方和转入方各算了一次要除以 2 去重。这两种写法都有人用关键是在报告里写清楚用的是哪一个并且前后一致。我个人的习惯是用不带系数 2 的版本因为它更容易解释分子是全区发生转出的总面积分母是全区总面积比值就是整体上有多少比例的用地发生了改变再除以年数。继续用算例。假设研究区期初总面积 20000 公顷十年间各地类转出面积之和为 2430 公顷则$$LC \frac{2430}{20000} \times \frac{1}{10} \times 100% 1.215%$$意思是这片区域平均每年有约 1.2% 的用地发生了类型转换。这个数字可以和同区域的其他时期对比也可以和相邻区域对比用来判断哪个区域的土地利用更活跃。注意算综合动态度时分子千万别把转移矩阵对角线上的未变化面积加进去。对角线是没变的部分加进去分母就成了总面积乘以年数结果会小得莫名其妙。这个错我见过不止一次。3.3 字段计算器实操与算例复核有了公式落地到 ArcGIS 里就是在属性表里做字段计算。假设你已经把各地类的期初面积和期末面积汇总到一张表里表里有字段 A_start期初面积公顷、A_end期末面积公顷、Years时段年数那么在字段计算器里新建一个 DOUBLE 型字段 K用 Python 解析器写# 单一土地利用动态度%每年 (!A_end! - !A_start!) / !A_start! / !Years! * 100注意 ArcGIS Desktop 用的是 Python 2 语法ArcGIS Pro 用的是 Python 3两者在字段计算器里的表达式写法基本一致但涉及字符串处理时有差异。上面这个纯数值计算两边通用。如果想在脚本里批量跑用 arcpy 更省事。下面是矢量路线里从叠加到汇总统计的完整思路代码import arcpy arcpy.env.workspace rD:\LUCC\LUCC.gdb arcpy.env.overwriteOutput True # 1. 两期数据相交保留两期地类编码字段 arcpy.analysis.Intersect( [LU2010, LU2020], intersect_1020, join_attributesALL, cluster_tolerance0.001 Meters ) # 2. 添加面积字段按公顷计算 arcpy.management.AddField(intersect_1020, AREA_HA, DOUBLE) arcpy.management.CalculateGeometryAttributes( intersect_1020, [[AREA_HA, AREA]], area_unitHECTARES ) # 3. 按两期地类编码分组统计面积之和 arcpy.analysis.Statistics( intersect_1020, stat_1020, [[AREA_HA, SUM]], [DLBM_2010, DLBM_2020] )这段代码跑完stat_1020 表里就是转移矩阵的原始数据每一行是一个期初地类—期末地类组合加上对应的面积之和。把它导到 Excel 做透视表行放 DLBM_2010列放 DLBM_2020值放 SUM_AREA_HA一张标准的转移矩阵就出来了。这里有个关键细节相交之后两期的同名字段会自动加后缀比如两期都叫 DLBM输出可能是 DLBM 和 DLBM_1名字不直观还容易搞混。稳妥的做法是在相交之前就把字段重命名好比如统一改成 DLBM10 和 DLBM20或者用相交工具的字段映射Field Map手动指定输出字段名。矩阵算出来之后我建议做一次交叉复核把矩阵每一行的所有格子加起来应该等于该地类期初的总面积每一列加起来应该等于该地类期末的总面积。如果哪一行对不上说明叠加过程中有几何丢失回去检查修复几何和容差设置。4. 转移矩阵与变化图斑动态度背后的支撑数据4.1 矢量路线相交叠加加汇总统计矢量路线的核心就是上面那段 arcpy 脚本的逻辑手工操作用户可以这样走Intersect相交Analysis Tools Overlay Intersect输入两期数据输出类型默认 ALL容差如果数据是投影坐标系且经过了拓扑修复可以先用默认或者设成 0.001 米。Add Geometry Attributes添加几何属性勾选 AREA单位选 HECTARES或者用 Calculate Geometry 逐字段算。Summary Statistics汇总统计统计字段选面积求和分组字段选两期的地类编码。导出 Excel 做透视表。这四步里最容易出问题的是第二步和第三步之间。如果相交输出的要素数量很大几十万以上计算几何属性会很慢建议先用 ArcGIS Pro 的并行处理环境变量提高效率或者改用栅格路线。还有一个提升精度的技巧在相交之前对两期数据分别做一次融合。融合能把同一地类的相邻图斑合并大幅减少相交输出的要素数量同时保留面积信息。融合会丢失原有图斑的独立边界信息但对于算转移矩阵来说只需要地类之间的面积关系不影响结果。这个操作能让相交环节的速度提升数倍。4.2 栅格路线Combine 与 Tabulate Area如果数据是栅格格式或者矢量太大想转成栅格快速评估有两条现成的路第一条是 Tabulate Area制表面积这个是 Spatial Analyst 里的工具能直接输出两期栅格之间的交叉面积表输出本身就是一张近似转移矩阵的表非常省事import arcpy from arcpy.sa import * arcpy.CheckOutExtension(Spatial) # 输入栅格、地类编码字段、输出表、像元大小或参考投影 out TabulateArea( LU2010_img, Value, LU2020_img, Value, tab_area, 30 ) out.save(tab_area_table)Tabulate Area 的输出单位取决于输入的坐标系统务必先投影到等面积坐标系否则输出的面积数字是错的。栅格的分辨率决定精度30 米分辨率的栅格每个像元代表 0.09 公顷小图斑会被吞掉这一点要在报告里说明。第二条是 Combine组合把两期栅格组合成一个新栅格每个唯一的组合值对应一种转移类型然后统计各组合值的像元数乘以像元面积也能得到转移矩阵。Combine 的好处是逻辑清晰、易于后续做变化图坏处是组合值的数量可能很多属性表巨大。# 假设两期栅格的地类编码都在 1-20 之间 combined Raster(LU2010_img) * 100 Raster(LU2020_img) combined.save(combine_1020) # 建立属性表后即可按 Value 统计像元数 arcpy.management.BuildRasterAttributeTable(combine_1020)这里有个小技巧编码乘以 100 再相加相当于把期初编码放在了百位以上期末编码放在了个位和十位后期解析组合值的时候直接拆位就能还原出转移类型不用去查对照表。提示栅格路线算出的转移矩阵边界位置会受像元分辨率影响比如一条 30 米宽的地类边界在栅格里可能被压成一条线或者干脆消失。矢量路线的面积统计更贴合实际图斑两者结果通常有差异这是正常的报告里注明方法即可不要为了对齐两个数字去硬改数据。4.3 变化图斑提取与专题制图出图数字算完了图还得做。转移矩阵和动态度在空间上怎么体现核心是变化图斑提取。矢量数据的做法是在相交结果里筛选出 DLBM10 不等于 DLBM20 的要素这些就是变化图斑。筛选出来的图层按从什么变成什么这个组合字段做唯一值符号化比如耕地转建设用地用红色林地转耕地用绿色一张变化专题图就出来了。栅格数据的做法类似用组合值筛选出不等于自身编码的那部分或者用 Con 函数做条件运算# 提取发生变化的区域变化区域为 1未变化为 0 changed Con(Raster(LU2010_img) ! Raster(LU2020_img), 1, 0) changed.save(changed_area)除了变化图斑本身还有一类图很受欢迎变化强度分布图。做法很简单用 Create Fishnet 生成一个规则格网比如 1 公里 × 1 公里然后用格网去裁剪变化图斑统计每个格网内的变化面积占格网总面积的比例最后按比例做分级设色。这张图能直观看出哪些区域变化密集比单纯的变化图斑图更有信息量在论文和汇报里都很吃香。制图这一步有几个经验地类符号色彩要有区分度别全用相近的绿读者一眼看不出区别如果地类符号太密导致图面糊成一团考虑做图例放大的局部小图而不是硬往主图里塞变化箭头或者流向图用 Excel 的条件格式热力图或者第三方图表工具画ArcGIS 里画桑基图比较费劲不必强求出图分辨率 300 dpi布局里带上比例尺、指北针和图例这是基本要求。5. 常见问题与排查技巧实录5.1 面积对不上、比例异常的几类原因做这行做久了你会发现结果不对的原因基本就那么几类整理成速查表方便对照现象最可能的原因排查方法解决方式动态度大得离谱超过 10%两期地类编码错位比对编码字段的取值集合统一编码体系后重算所有地类动态度都是 0分类字段类型不一致检查字段是文本还是数值统一类型或补齐位数总面积比实际少一块两期边界范围不一致叠加后看是否有空白带统一用同一份边界裁剪面积数字带有小数点很多位但明显偏小地理坐标系下算的面积查看数据坐标系投影到等面积坐标系重算转移矩阵行列和不等叠加时几何丢失逐行逐列求和核对修复几何、调小容差重跑小地类完全消失栅格分辨率太粗查看最小图斑面积提高分辨率或改走矢量这里单独说说编码错位这个坑。有一种很隐蔽的情况是两期数据编码本身没错但字段名不同比如一期叫 DLBM、一期叫 DLMC你在分组统计时选错了字段结果统计出来的是一个空表或者全是对角线。还有一种是编码里混了全角字符和半角字符肉眼看着一模一样实际不相等这种只能靠导出唯一值列表逐一比对来发现。5.2 叠加出碎屑、运行卡死的处理办法相交运算跑完之后输出要素数量暴增十倍二十倍属性表卡得打不开这是非常常见的场景。原因通常是两期数据的边界几何在局部不一致相交时在边界处产生了大量细长的碎屑多边形数据本身有重叠、自相交等几何问题容差设置过小导致本应合并的节点被拆开。处理办法按优先级排先融合后相交。这一步最能减少要素数量效果立竿见影。适当放大容差。投影坐标系下容差设成 0.001 米到 0.01 米之间通常够用设太小反而容易出碎屑。相交后按面积筛选。把面积小于阈值的碎屑多边形删掉或者按面积排序后合并到相邻大图斑。这一步会影响总面积所以删除前先算一遍删除的面积总量如果占比很小比如千分之一以内可以接受。改用 Pairwise Intersect。多核并行速度明显提升。大范围研究改用栅格路线。几千万个图斑的矢量叠加基本跑不动转栅格是唯一现实的选择。关于运行时进度框长时间不动这件事先别急着强制结束。相交分析在后台生成几何索引遇到大数据停顿几分钟是正常的。如果超过半小时还没动静检查一下是不是磁盘空间不足或者输出路径用了中文和特殊字符。这些看似无关的细节实际非常容易导致工具静默失败。5.3 结果表达与交付的几个细节数据算完了最后一步是让结果能被别人看懂。表格部分转移矩阵建议做成带条件格式的热力图对角线深色、非对角线按面积大小渐变一眼就能看出主要转移方向。表格里面积单位统一到公顷或平方公里小数点保留两位足够保留太多位反而显得不专业。图表部分各地类动态度做成柱状图正负用不同颜色区分综合动态度按时期做成折线图多期数据一对比趋势就出来了。如果只有两期数据折线图意义不大用柱状图对比各地类动态度更合适。图件部分至少出三张图——期初土地利用图、期末土地利用图、变化图斑分布图三张图的配色、图例位置、比例尺保持统一放在报告里才成体系。数据部分交付时把中间过程数据也打包尤其是相交结果和转移矩阵原始表。别人复核你的结果时最想知道的就是你这个数怎么来的有中间数据沟通成本能降一半。最后说一个我自己踩过的坑。有一次做多期动态度为了省事把三期数据两两相交算了两组动态度结果报告里前后两组数用的地类编码版本不一样被审的人一眼看出来。从那以后我养成了一个习惯所有输入数据、字段名、坐标系、分类版本先写进一个数据说明文档里算之前对照一遍。这个文档花十分钟写能避免后面几天的返工。土地利用动态分析这件事工具操作只占三成七成功夫在数据准备和口径统一上。公式背得再熟源数据是脏的结果照样不能用。我现在拿到数据的第一反应不是打开工具而是先看字段、看坐标系、看范围、看分类四样都对上了再动手心里才踏实。
ArcGIS两期土地利用动态图制作：转移矩阵与动态度计算实战
两期土地利用图摆在面前一张是几年前的一张是最新的任务只有一个算出这片区域的土地利用动态图。这个需求在 ArcGIS 里几乎是最常见的一类活儿——说它简单流程确实不复杂说它坑多也确实能让人来回折腾一整天。我做土地利用变化分析这些年见过太多人卡在第一步两期数据的分类编码对不上叠加出来一片乱码算出来的动态度要么大得离谱要么小数点后全是零。先把这件事说清楚土地利用动态图这个词在实际项目里其实是个笼统叫法它往往同时包含三块内容——各类地类的单一土地利用动态度每年变化百分之几、全区的综合土地利用动态度整体活跃程度以及支撑这两者的土地利用转移矩阵和变化图斑分布图。前两个是数字后两个是图和表四样东西凑齐才算把动态这件事讲明白。这篇文章面向的是手上有两期或多期土地利用矢量数据、会用 ArcGIS 基本工具、但没系统做过变化分析的人。如果你正在做土地资源调查、生态评估、国土空间规划前期分析或者毕业论文里要用到土地利用动态度这个指标下面这套流程可以直接照着走。我会把每一步为什么这么做、参数怎么定、哪里最容易翻车都讲透公式和算例一并给出。1. 先把土地利用动态图这五个字拆开看1.1 动态度、变化图、转移矩阵三个东西三条路线很多人一上来就问动态度在哪个工具里算这个问题本身就问偏了。ArcGIS 里没有哪个按钮叫计算动态度动态度是一个统计指标它的输入是各地类的面积输出是一个百分比数字ArcGIS 负责的是把面积统计出来剩下的四则运算可能要在属性表或者 Excel 里完成。具体拆开看单一土地利用动态度针对某一个地类比如耕地算它在一段时间内平均每年变化了多少个百分点。输入是这一类地类的期初面积和期末面积。综合土地利用动态度针对整个研究区把所有地类之间的转移面积加总除以期初总面积衡量区域整体的土地利用活跃程度。土地利用转移矩阵一张 n×n 的表格行是期初地类列是期末地类交叉格子里是对应的面积对角线是没变的部分。变化图斑图把所有发生了地类改变的多边形挑出来做成专题图直观展示变化发生在哪里。这四样东西的关系是转移矩阵是底层数据动态度是从矩阵里提炼出的指标变化图斑图是空间表达。路线上有矢量叠加和栅格叠加两条矢量精度高、边界贴合实际栅格速度快、适合大范围快速评估。我的习惯是研究区不大比如一个县以内或者要求图斑边界精确走矢量跨市跨省、动辄几千万个图斑走栅格。提示动态度这个词在不同文献里的公式写法不完全一样尤其是综合动态度有的分母带系数 2有的不带。写报告前先确定你参考的那篇文献用的是哪种别中途换公式导致前后数据对不上。1.2 数据口径不统一后面全白干真正常见的翻车点不在计算而在数据准备。土地利用数据的来源五花八门可能是某次国土调查的成果可能是遥感解译的分类结果也可能是从别人那里拷来的历史数据。这些数据的分类体系、编码位数、坐标系、边界范围经常不一致。我遇到过最典型的一次两期数据的耕地编码一期是01另一期是0101看着都是耕地但字段一比对全是错位叠加之后原本没变化的耕地被识别成了变化图斑动态度算出来高得离谱。还有一次是两期数据的行政边界范围差了半条街相交之后边缘出现一条几百米宽的空白带所有地类的面积都少算了一块。所以在动手之前必须先做三件事确认两期数据用的是同一套分类体系。如果一期是国标二级类、一期是一级类先归并到统一的层级再做。归并的原则是向上归并二级类合成一级类不会丢信息反过来则不行。确认编码字段类型一致都是文本型或者都是数值型位数补齐。文本01和数值 1 在连接和分组统计时是两个东西。确认研究区边界统一用同一份行政边界去裁剪两期数据或者先用相交工具取两期的公共范围之后所有统计都基于这个公共范围。这三件事做完可能只花你二十分钟但能省下一整天的返工。2. 数据准备把两期数据磨成能算的状态2.1 地类编码对齐与融合去碎两期数据的分类对齐之后下一步是融合Dissolve。为什么要融合因为很多土地利用数据的图斑碎得离谱同一块连片的耕地可能被切成七八个多边形中间还夹杂着几十平方米的细碎图斑。直接拿去做叠加分析输出要素数量会呈指数级增长跑起来慢结果也难看。融合的操作路径是按地类编码字段做 Dissolve在统计字段里保留面积字段求和其实面积后面重新算也行。融合前建议先删掉面积小于某个阈值的碎图斑比如小于 0.1 公顷的图斑根据研究精度决定。这一步要慎重因为删多了会影响总面积我的做法是把碎图斑合并到相邻的最大图斑里而不是直接删除这样总面积不变。还有一种情况是同一个地类编码下存在重复的、相互重叠的图斑。这种数据在叠加时会产出大量零面积的碎屑必须先处理掉。处理办法是跑一遍拓扑检查里的不能重叠规则把重叠部分找出来修正。ArcGIS Pro 里还有个提速技巧如果两期数据都很大用Pairwise Intersect代替普通的 Intersect多核并行速度能快好几倍结果基本一致。这个工具在 Analysis Tools Pairwise Overlay 下面很多人不知道它的存在。2.2 拓扑检查与几何修复的实操套路土地利用数据要做叠加分析几何质量必须过关否则报错或者出怪结果。我固定会做这几项检查检查项工具/规则常见问题处理方式自相交拓扑规则不能自相交多边形边界自己打结用修复几何工具自动修重叠压盖拓扑规则不能重叠两块地类互相压盖手动编辑或用联合后重建空隙拓扑规则不能有空隙图斑之间有条缝补充缝隙多边形或融合尖锐角几何检查或专用插件夹角过小导致后续叠加报错删除节点或局部重建伪结点、悬挂点拓扑规则线要素未闭合编辑捕捉修正尖锐角这个东西平时不起眼但在做大批量叠加的时候特别要命。角度过小的节点会让叠加算法在某些位置产生异常细长的多边形面积算出来是个极小的正数或者干脆是零。市面上有一些专门检查尖锐角的插件如果不装插件也可以自己用字段计算器算相邻边的夹角或者干脆把夹角小于某个阈值比如 1 度的节点在编辑时手动拉开。注意修复几何Repair Geometry能自动处理自相交、空几何、环方向错误这类问题但它不会处理重叠和空隙那两类必须靠拓扑和人工编辑。还有一点容易被忽略地理坐标系下的数据不要直接做叠加分析。虽然工具能跑但容差单位是度精度完全不可控叠加结果可能出现肉眼可见的偏移。2.3 坐标系与面积口径为什么一定要投影坐标系这是我最想强调的一条。只要涉及面积计算必须用投影坐标系且最好是等面积投影。用地理坐标系经纬度算面积ArcGIS 虽然会给结果但那是在球面上按测地线方式算的单位是平方米没错可一旦数据范围跨了几个纬度带误差会明显放大。更麻烦的是很多人在叠加分析之后直接读 Shape_Area 字段而这个字段是在数据本身的坐标系下算的如果原数据是地理坐标系Shape_Area 的单位是平方度这个数没有任何物理意义。我的标准做法是根据研究区位置选一个合适的投影坐标系。中国范围内常用的是 CGCS2000 的高斯-克吕格投影按 3 度带或 6 度带分带如果研究区跨带就换成 Albers 等面积投影中央经线和双标准纬线按研究区范围设定。两期数据投影到同一个坐标系投影参数完全一致。面积字段统一用公顷因为土地利用动态度、转移矩阵的常规表达单位就是公顷。平方米转公顷除以 10000直接用 Calculate Geometry 时选 HECTARES 就行。选 3 度带还是 6 度带也有讲究。3 度带变形更小适合精度要求高的县级尺度研究6 度带带号少适合省级以上范围。如果你不确定用 Albers 等面积投影最保险因为动态度和面积直接挂钩等面积投影能保证面积不变形。投影做完建议先算一遍两期数据的各地类总面积各自核对一下看看是不是和研究区公布的统计面积对得上。这一步是体检如果面积对不上后面所有计算都不可信。对不上的原因通常是边界裁剪不一致、碎图斑没处理干净、或者分类编码有重复。3. 单一动态度与综合动态度公式、参数与算例3.1 单一土地利用动态度怎么算单一土地利用动态度的公式不复杂但每一个参数的含义要弄清楚$$K_i \frac{U_b - U_a}{U_a} \times \frac{1}{T} \times 100%$$其中$K_i$ 是第 i 类地类的单一土地利用动态度单位是百分比每年$U_a$ 是该地类在研究期初的面积单位公顷$U_b$ 是该地类在研究期末的面积单位公顷$T$ 是研究时段长度单位年。这个公式的物理含义很直白先算整个研究期内这一类地类变化了多少比例再除以年数得到平均每年的变化率。正数说明面积增加负数说明面积减少绝对值越大说明变化越剧烈。举个算例。某研究区 2010 年耕地面积 10240.5 公顷2020 年耕地面积 8632.7 公顷研究时段 $T 10$ 年。代入公式$$K_{耕地} \frac{8632.7 - 10240.5}{10240.5} \times \frac{1}{10} \times 100% \frac{-1607.8}{10240.5} \times 0.1 \times 100% \approx -1.57%$$也就是说这十年里耕地年均减少约 1.57%。这个数字放在报告里配上一句耕地年均减少 1.57%主要流向建设用地和林地说服力就出来了。关于 $T$ 的取值有个细节。如果期初数据是 2010 年末的期末是 2020 年末的那 $T$ 就是 10 年整。但如果期初是 2010 年 6 月、期末是 2020 年 3 月严格算的话 $T$ 应该是 9.75 年左右。实际项目里大多数人是按年份差直接取整数这个做法可以接受但如果时段很短比如两年T 取整带来的误差就不能忽略了这时候建议按实际月份折算。3.2 综合动态度分母到底取什么综合土地利用动态度的公式版本特别多这里给一个在中文文献里最常见、也最容易被审稿人接受的写法$$LC \frac{\sum_{i1}^{n} \Delta LU_{i-j}}{\sum_{i1}^{n} LU_i} \times \frac{1}{T} \times 100%$$参数含义$LC$ 是综合土地利用动态度$\Delta LU_{i-j}$ 是第 i 类地类转为其他地类的面积之和的绝对值注意这里取的是转出面积不包括自身不变的对角线部分$\sum LU_i$ 是研究期初所有地类的面积总和也就是研究区总面积$T$ 同样是时段长度。关于分母还有一种写法是分母乘 2理由是变化面积被转出方和转入方各算了一次要除以 2 去重。这两种写法都有人用关键是在报告里写清楚用的是哪一个并且前后一致。我个人的习惯是用不带系数 2 的版本因为它更容易解释分子是全区发生转出的总面积分母是全区总面积比值就是整体上有多少比例的用地发生了改变再除以年数。继续用算例。假设研究区期初总面积 20000 公顷十年间各地类转出面积之和为 2430 公顷则$$LC \frac{2430}{20000} \times \frac{1}{10} \times 100% 1.215%$$意思是这片区域平均每年有约 1.2% 的用地发生了类型转换。这个数字可以和同区域的其他时期对比也可以和相邻区域对比用来判断哪个区域的土地利用更活跃。注意算综合动态度时分子千万别把转移矩阵对角线上的未变化面积加进去。对角线是没变的部分加进去分母就成了总面积乘以年数结果会小得莫名其妙。这个错我见过不止一次。3.3 字段计算器实操与算例复核有了公式落地到 ArcGIS 里就是在属性表里做字段计算。假设你已经把各地类的期初面积和期末面积汇总到一张表里表里有字段 A_start期初面积公顷、A_end期末面积公顷、Years时段年数那么在字段计算器里新建一个 DOUBLE 型字段 K用 Python 解析器写# 单一土地利用动态度%每年 (!A_end! - !A_start!) / !A_start! / !Years! * 100注意 ArcGIS Desktop 用的是 Python 2 语法ArcGIS Pro 用的是 Python 3两者在字段计算器里的表达式写法基本一致但涉及字符串处理时有差异。上面这个纯数值计算两边通用。如果想在脚本里批量跑用 arcpy 更省事。下面是矢量路线里从叠加到汇总统计的完整思路代码import arcpy arcpy.env.workspace rD:\LUCC\LUCC.gdb arcpy.env.overwriteOutput True # 1. 两期数据相交保留两期地类编码字段 arcpy.analysis.Intersect( [LU2010, LU2020], intersect_1020, join_attributesALL, cluster_tolerance0.001 Meters ) # 2. 添加面积字段按公顷计算 arcpy.management.AddField(intersect_1020, AREA_HA, DOUBLE) arcpy.management.CalculateGeometryAttributes( intersect_1020, [[AREA_HA, AREA]], area_unitHECTARES ) # 3. 按两期地类编码分组统计面积之和 arcpy.analysis.Statistics( intersect_1020, stat_1020, [[AREA_HA, SUM]], [DLBM_2010, DLBM_2020] )这段代码跑完stat_1020 表里就是转移矩阵的原始数据每一行是一个期初地类—期末地类组合加上对应的面积之和。把它导到 Excel 做透视表行放 DLBM_2010列放 DLBM_2020值放 SUM_AREA_HA一张标准的转移矩阵就出来了。这里有个关键细节相交之后两期的同名字段会自动加后缀比如两期都叫 DLBM输出可能是 DLBM 和 DLBM_1名字不直观还容易搞混。稳妥的做法是在相交之前就把字段重命名好比如统一改成 DLBM10 和 DLBM20或者用相交工具的字段映射Field Map手动指定输出字段名。矩阵算出来之后我建议做一次交叉复核把矩阵每一行的所有格子加起来应该等于该地类期初的总面积每一列加起来应该等于该地类期末的总面积。如果哪一行对不上说明叠加过程中有几何丢失回去检查修复几何和容差设置。4. 转移矩阵与变化图斑动态度背后的支撑数据4.1 矢量路线相交叠加加汇总统计矢量路线的核心就是上面那段 arcpy 脚本的逻辑手工操作用户可以这样走Intersect相交Analysis Tools Overlay Intersect输入两期数据输出类型默认 ALL容差如果数据是投影坐标系且经过了拓扑修复可以先用默认或者设成 0.001 米。Add Geometry Attributes添加几何属性勾选 AREA单位选 HECTARES或者用 Calculate Geometry 逐字段算。Summary Statistics汇总统计统计字段选面积求和分组字段选两期的地类编码。导出 Excel 做透视表。这四步里最容易出问题的是第二步和第三步之间。如果相交输出的要素数量很大几十万以上计算几何属性会很慢建议先用 ArcGIS Pro 的并行处理环境变量提高效率或者改用栅格路线。还有一个提升精度的技巧在相交之前对两期数据分别做一次融合。融合能把同一地类的相邻图斑合并大幅减少相交输出的要素数量同时保留面积信息。融合会丢失原有图斑的独立边界信息但对于算转移矩阵来说只需要地类之间的面积关系不影响结果。这个操作能让相交环节的速度提升数倍。4.2 栅格路线Combine 与 Tabulate Area如果数据是栅格格式或者矢量太大想转成栅格快速评估有两条现成的路第一条是 Tabulate Area制表面积这个是 Spatial Analyst 里的工具能直接输出两期栅格之间的交叉面积表输出本身就是一张近似转移矩阵的表非常省事import arcpy from arcpy.sa import * arcpy.CheckOutExtension(Spatial) # 输入栅格、地类编码字段、输出表、像元大小或参考投影 out TabulateArea( LU2010_img, Value, LU2020_img, Value, tab_area, 30 ) out.save(tab_area_table)Tabulate Area 的输出单位取决于输入的坐标系统务必先投影到等面积坐标系否则输出的面积数字是错的。栅格的分辨率决定精度30 米分辨率的栅格每个像元代表 0.09 公顷小图斑会被吞掉这一点要在报告里说明。第二条是 Combine组合把两期栅格组合成一个新栅格每个唯一的组合值对应一种转移类型然后统计各组合值的像元数乘以像元面积也能得到转移矩阵。Combine 的好处是逻辑清晰、易于后续做变化图坏处是组合值的数量可能很多属性表巨大。# 假设两期栅格的地类编码都在 1-20 之间 combined Raster(LU2010_img) * 100 Raster(LU2020_img) combined.save(combine_1020) # 建立属性表后即可按 Value 统计像元数 arcpy.management.BuildRasterAttributeTable(combine_1020)这里有个小技巧编码乘以 100 再相加相当于把期初编码放在了百位以上期末编码放在了个位和十位后期解析组合值的时候直接拆位就能还原出转移类型不用去查对照表。提示栅格路线算出的转移矩阵边界位置会受像元分辨率影响比如一条 30 米宽的地类边界在栅格里可能被压成一条线或者干脆消失。矢量路线的面积统计更贴合实际图斑两者结果通常有差异这是正常的报告里注明方法即可不要为了对齐两个数字去硬改数据。4.3 变化图斑提取与专题制图出图数字算完了图还得做。转移矩阵和动态度在空间上怎么体现核心是变化图斑提取。矢量数据的做法是在相交结果里筛选出 DLBM10 不等于 DLBM20 的要素这些就是变化图斑。筛选出来的图层按从什么变成什么这个组合字段做唯一值符号化比如耕地转建设用地用红色林地转耕地用绿色一张变化专题图就出来了。栅格数据的做法类似用组合值筛选出不等于自身编码的那部分或者用 Con 函数做条件运算# 提取发生变化的区域变化区域为 1未变化为 0 changed Con(Raster(LU2010_img) ! Raster(LU2020_img), 1, 0) changed.save(changed_area)除了变化图斑本身还有一类图很受欢迎变化强度分布图。做法很简单用 Create Fishnet 生成一个规则格网比如 1 公里 × 1 公里然后用格网去裁剪变化图斑统计每个格网内的变化面积占格网总面积的比例最后按比例做分级设色。这张图能直观看出哪些区域变化密集比单纯的变化图斑图更有信息量在论文和汇报里都很吃香。制图这一步有几个经验地类符号色彩要有区分度别全用相近的绿读者一眼看不出区别如果地类符号太密导致图面糊成一团考虑做图例放大的局部小图而不是硬往主图里塞变化箭头或者流向图用 Excel 的条件格式热力图或者第三方图表工具画ArcGIS 里画桑基图比较费劲不必强求出图分辨率 300 dpi布局里带上比例尺、指北针和图例这是基本要求。5. 常见问题与排查技巧实录5.1 面积对不上、比例异常的几类原因做这行做久了你会发现结果不对的原因基本就那么几类整理成速查表方便对照现象最可能的原因排查方法解决方式动态度大得离谱超过 10%两期地类编码错位比对编码字段的取值集合统一编码体系后重算所有地类动态度都是 0分类字段类型不一致检查字段是文本还是数值统一类型或补齐位数总面积比实际少一块两期边界范围不一致叠加后看是否有空白带统一用同一份边界裁剪面积数字带有小数点很多位但明显偏小地理坐标系下算的面积查看数据坐标系投影到等面积坐标系重算转移矩阵行列和不等叠加时几何丢失逐行逐列求和核对修复几何、调小容差重跑小地类完全消失栅格分辨率太粗查看最小图斑面积提高分辨率或改走矢量这里单独说说编码错位这个坑。有一种很隐蔽的情况是两期数据编码本身没错但字段名不同比如一期叫 DLBM、一期叫 DLMC你在分组统计时选错了字段结果统计出来的是一个空表或者全是对角线。还有一种是编码里混了全角字符和半角字符肉眼看着一模一样实际不相等这种只能靠导出唯一值列表逐一比对来发现。5.2 叠加出碎屑、运行卡死的处理办法相交运算跑完之后输出要素数量暴增十倍二十倍属性表卡得打不开这是非常常见的场景。原因通常是两期数据的边界几何在局部不一致相交时在边界处产生了大量细长的碎屑多边形数据本身有重叠、自相交等几何问题容差设置过小导致本应合并的节点被拆开。处理办法按优先级排先融合后相交。这一步最能减少要素数量效果立竿见影。适当放大容差。投影坐标系下容差设成 0.001 米到 0.01 米之间通常够用设太小反而容易出碎屑。相交后按面积筛选。把面积小于阈值的碎屑多边形删掉或者按面积排序后合并到相邻大图斑。这一步会影响总面积所以删除前先算一遍删除的面积总量如果占比很小比如千分之一以内可以接受。改用 Pairwise Intersect。多核并行速度明显提升。大范围研究改用栅格路线。几千万个图斑的矢量叠加基本跑不动转栅格是唯一现实的选择。关于运行时进度框长时间不动这件事先别急着强制结束。相交分析在后台生成几何索引遇到大数据停顿几分钟是正常的。如果超过半小时还没动静检查一下是不是磁盘空间不足或者输出路径用了中文和特殊字符。这些看似无关的细节实际非常容易导致工具静默失败。5.3 结果表达与交付的几个细节数据算完了最后一步是让结果能被别人看懂。表格部分转移矩阵建议做成带条件格式的热力图对角线深色、非对角线按面积大小渐变一眼就能看出主要转移方向。表格里面积单位统一到公顷或平方公里小数点保留两位足够保留太多位反而显得不专业。图表部分各地类动态度做成柱状图正负用不同颜色区分综合动态度按时期做成折线图多期数据一对比趋势就出来了。如果只有两期数据折线图意义不大用柱状图对比各地类动态度更合适。图件部分至少出三张图——期初土地利用图、期末土地利用图、变化图斑分布图三张图的配色、图例位置、比例尺保持统一放在报告里才成体系。数据部分交付时把中间过程数据也打包尤其是相交结果和转移矩阵原始表。别人复核你的结果时最想知道的就是你这个数怎么来的有中间数据沟通成本能降一半。最后说一个我自己踩过的坑。有一次做多期动态度为了省事把三期数据两两相交算了两组动态度结果报告里前后两组数用的地类编码版本不一样被审的人一眼看出来。从那以后我养成了一个习惯所有输入数据、字段名、坐标系、分类版本先写进一个数据说明文档里算之前对照一遍。这个文档花十分钟写能避免后面几天的返工。土地利用动态分析这件事工具操作只占三成七成功夫在数据准备和口径统一上。公式背得再熟源数据是脏的结果照样不能用。我现在拿到数据的第一反应不是打开工具而是先看字段、看坐标系、看范围、看分类四样都对上了再动手心里才踏实。
