面向土地覆盖变化信息局域精度预测的逻辑回归协变量选用的直接策略与组合策略

张景雄，万 月，梅莹莹

(1. 武汉大学 测绘学院，湖北 武汉 430079；2. 武汉大学 遥感信息工程学院，湖北 武汉 430079；3.地球空间信息技术协同创新中心，湖北 武汉 430079)

0 引 言

反映地表覆盖类型的空间分布及其动态变化特征的土地覆盖（或称地表覆盖）信息产品日益积累，其在资源环境规划和管理中具有越来越重要的作用。得益于时间维度信息的富集，地表覆盖变化信息产品对于地理国情监测、生态多样性等具有重要意义[1]。精确的地表覆盖变化信息有助于人们更好地理解地表过程，开展地球系统的科学研究，并为地理国情监测和诸多地学应用提供信息服务。地理国情监测通过对地理环境进行感知、统计和分析，服务于国计民生[2-3]。

有关土地覆盖或类似的专题分类信息精度的研究也称为真实性验证，在这一领域国内外学者做出了许多重要工作[4]。1991年，Congalton[5]提出了分类精度的度量指标PCC(percentage correctly classified)和误差矩阵。为了将精度指标应用到具体区域，即从全部问题域到局部，McGwire[6]提出报告用户精度时应考虑具体小区域的精度。为计算用户精度，De Wit[7]使用设置参考点的方式，以荷兰为研究区域进行真实性验证。为了避免设置参考点方式对精度计算的限制，一些学者探讨了基于模型的精度空间分布估计方法。Steele等[8]使用克里金方法从样本点空间内插得到了误分类概率分布图。Smith等[9]使用逻辑回归方法评价了地块大小和邻域异质性对单一像素分类精度的影响。2004年，Van oort等[10]考虑了地表类别、邻域异质性、邻域同质性、地块大小、熵和蔓延度等多种因素，分别计算了对应的逻辑回归模型参数，并评价了各因素对地表精度估计的影响程度。

土地覆盖变化信息的精度研究也越来越受到重视[11]。Burnicki[12]于2011年提出了一种方法，通过综合分析时间序列图和景观结构来估计像素级的分类错误概率空间分布。2016年，陈军等[13]基于景观指数，对精度验证使用的抽样方法加以改进，分别基于验证区域、区内地类(土地景观的空间异质性)和抽样格网三个层级，为大范围精度验证提供了定量化的自适应抽样方法。Zhang and Mei (2016)[14]应用逻辑回归方法估计了湖北省洪湖市某子区2012～2013地表覆盖变化的像素级的分类精度。

本文对地表覆盖变化信息的局域精度预测方法作了进一步的探究。前已述及，变化类别的像素级的正确分类概率可以通过逻辑回归来计算,选用基于像素点所在3×3邻域的类别的空间发生模式所计量的景观指数和基于区内地类邻域的景观指数（如面积）作为回归建模的协变量。这里，景观指数可以有两种方式作为协变量，一是直接基于待验证的地表覆盖变化地图，二是组合利用变化前后两时相对应的地表覆盖地图。本研究通过实验，比较协变量选用的直接和组合策略的精度预测效果，评判两种策略的预测效果和适用性。

论文第一部分说明研究区域和实验数据。第二部分描述逻辑回归估计局域土地覆盖精度的基本原理。第三部分分析实验结果。最后是结论和展望。

1 研究区域与研究数据

1.1 研究区域

本文选择武汉市中心城区约10 km×10 km的区域作为研究区域，经度范围为114°1′40″E～114°5′40″E，纬度范围为30°3′35″N～30°8′39″。选取的遥感影像为ETM+数据，像素尺寸为30 m×30 m，研究区域包括313×312个像素，时间为2011年12月和2013年8月两期影像。

1.2 类别体系

考虑到研究区域的景观特点，覆盖类型大致包括植被（农田、草地、树木）、裸土（未开发土壤、河岸/湖岸滩地）、不透水表面（道路、桥梁、公共设施等建设用地、建筑物、居民地）和水体（水稻田、湖泊、河流、鱼塘等）。因此，为选择代表性和分离性较强的地表覆盖类别，选用分类体系（不透水表面、水体、裸土、植被）在后文中分别以英文首字母I、W、S、V代替。

1.3 土地覆盖变化信息提取

利用分类后处理的方法（post-classification）获得土地覆盖变化地图，如图1所示。考虑到2011年11月～2013年8月研究区域可能发生的地理景观变化，去除了不可能发生的(变化)类别和所占像素点过少的类别，比如水体->不透水表面，保留共9种类别（包括4种未变化类别和5种变化类别）。4类未变化类别为5种变化类别包括：IS：不透水表面->裸地；IV：不透水表面->植被；SI：裸地->不透水表面；SV：裸地->植被；VI：植被->不透水表面。

图1 土地覆盖与土地覆盖变化图：(a)-(b)时相1与时相2的土地覆盖图，(c)土地覆盖变化Fig.1 Maps showing: (a)-(b) land cover at Time 1 and Time 2,respectively, and (c) land cover change

1.4 采样方案与样本数据采集

根据空间采样原理[15]，在地表覆盖变化图中抽取样本数据，包括训练样本和测试样本。前者用于建立逻辑回归模型，后者用于比较两种策略预测结果好坏。为了能准确表达用户精度，训练样本和测试样本的选取均须反映土地覆盖变化类别的分布。本文采用基于分层随机抽样的原则，抽取符合要求的训练样本：

式中，p表示期望达到的精度，Z为正态分布95%置信区间的标准差，Z=1.96。现取p=77.12% p=0.7712，设定允许误差E=0.05，由（1）式计算得到N为304。但在实际应用中，每类别的样本数不能过少，根据经验设置最小类别数为20，调整后的训练样本总数为322，根据同样的抽样原理，抽取同样数量的测试样本。表1给出训练样本和测试样本中各个变化类别的个数，图2是训练样本和测试样本的分布。

表1 样本点数量Tab.1 The number of sample points

图2 训练和测试样本分布图Fig.2 Distribution map of training and testing sample points

2 方法

2.1 景观指数计算

采用下述景观指数：

CLASS: 像素所属类别，使用二值变量来表示。变化类别共有9种类别，用8个二值变量表示，因为有一冗余；如类别1可表示为[10 000 000]；单时相影像上分别有4种类别，各自用3个二值变量表示，则组合变量由6个二值变量表示。如若time1的类别为1，time2的类别为4，它们的组合可表示为[100 000]。

HOM(homogeneity): 局域同质性，像素所在3×3邻域内与中心像素类别相同的像素个数。

HET(heterogeneity): 局域异质性，即像素所在3×3邻域所包含的不同类别数量。

LPS(log10(Patch size)): patch size为像素所在相连区域包含的元素个数（此处相连为4邻域相连），再对其取对数。此处的计算方法利用了计算机形态学相关知识，首先从分类图中分离出各类别的二值图像，再使用函数bwconncomp分别计算各二值图像相连区域的像素个数和各斑块包含的像素地址列表，最后确定各元素所属斑块面积。本部分操作使用matlab完成。

鉴于本文研究的方法是面向局域的精度预测，上述景观指数的量化的空间尺度包括单像素、像素所在3×3邻域以及像素所在相连通的区域，类似于Zhang and Mei (2016)[14]。

2.2 逻辑回归

建立分类精度与各景观指数之间的逻辑回归模型，从而得到逐像素的分类正确概率，可显示为分类精度表面。回归关系式可由式（2）表示，式中xi为各协变量，βi为对应于各协变量的系数，β0为截距。本文选择CLASS、HOM、HET、LPS等景观指数作为协变量，P为回归精度。

选取统计量Deviance来进行模型选择。它通过比较包含或不包含待检验协变量的模型的对数似然函数，计算待评估协变量对模型的有用性。当样本量较大时，Deviance近似满足自由度为上述待检验变量个数（相对于不包含这些变量的基准模型，这里的待检验变量个数是待检验模型所增加的变量个数）的χ2分布。

2.3 直接策略与组合策略模型介绍

为了探究土地覆盖变化信息局部精度与景观指数之间的定量关系，可考虑前述的两种协变量选用策略。直接策略是直接利用所测试的地表覆盖变化地图的CLASS、HOM、HET、LPS作为协变量，见表2。组合策略利用两个单时相地表覆盖地图的协变量的组合，即CLASS1、CLASS2、HOM1、HOM2、HET1、HET2、LPS1、LPS2，见表3。

表2 直接策略候选模型描述Tab.2 Descriptions of candidate models in the direct strategy

表3 组合策略候选模型描述Tab.3 Descriptions of candidate models in the composite strategy

在表2中模型0不包含任何协变量，只有截距。模型1a、1b、1c、1d中分别包含CLASS、HOM、HET、LPS。模型2a、2b、2c在模型1a的基础上，分别增加协变量HOM、HET、LPS。模型3a，3b则相比于模型2a，分别增加HET、LPS。表3与表2类似。

2.4 模型验证

本文选择以下评价指标来比较模型的预测精度：均值ME(mean error), 平均绝对误差MAE(mean absolute error)均方根误差RMSE(root mea〈n square error)。i(xj)代表测试样本点xj实际精度，(xj)为点xj估计值。

3 实验结果与分析

3.1 逻辑回归建模

根据表2、表3分别建立直接策略和组合策略的逻辑回归模型；利用1.4节抽取的训练样本估计回归模型的系数，直接策略和组合策略逻辑回归模型系数分别见表4、表5。表4、表5中的每一行分别对应表2、表3中对应项建立模型。以表4为例，第一行对应模型0，β0为截距。第二行对应模型1a, β1～β8为协变量CLASS中每一位二值变量的系数。

表4 直接策略候选模型系数Tab.4 Estimated regression coefficients for candidate models, the direct strategy

表5 组合策略候选模型系数Tab.5 Estimated regression coefficients for candidate models,the composite strategy

3.2 模型选择

为了比较不同模型拟合的优劣，根据（3）式选取统计量Deviance，采用似然比检验法对模型进行显著性检验，用于选择含有最多有效协变量的模型。直接策略和组合策略逻辑回归模型选择的分析过程分别见表6、表7。

表6 直接策略模型选择的x2检验Tab.6 Chi-square tests for model selection,the direct strategy

表7 组合策略模型选择的x2检验Tab.7 Chi-square tests for model selection, the composite strategy

由表6可知，模型2a、2b、2c、2d中分别相对于模型0增加了CLASS、HOM、HET、LPS4个变量，这4个模型相对于基准模型的改善作用可以分别用Deviance差值来衡量。由表6可以看出，CLASS的增加对改善模型的作用最大，且满足自由度为8，α=0.01的χ2检验。因此，增加CLASS作为优选模型的第一个协变量。同理，相对于模型1a的比较，增加HOM作为第二个协变量。在模型3a、3b相对于模型2a的比较中，可看出变量HET、LPS的增加均对改善模型不具显著作用。因此可得直接策略中协变量CLASS、HOM对改善模型有显著作用，故直接策略中优选模型为2a。

同样，由表7可知，模型2a，2b，2c,2d中分别相对于模型0增加了CLASS1,2、HOM1,2、HET1,2、LPS1,2这几个变量，这4个模型相对于基准模型的改善作用分别可以用Deviance差值来衡量。可看出，CLASS1,2的增加对改善模型的作用最大，且CLASS1,2的增加满足自由度为6，α=0.01的χ2检验。因此增加CLASS1,2作为优选模型的第一组协变量。同理，相对于模型1a的比较过程，增加HOM1,2作为第二组协变量。在模型3a、3b相对于模型2a的比较中，可看出变量HET1,2、LPS1,2的增加均对改善模型不具显著作用。因此，组合策略协变量CLASS1,2、HOM1,2对改善模型有显著作用，故得到组合策略优选模型为2a。

3.3 策略比较

本部分根据2.4节方法对两策略精度估计结果进行定量评价。首先，分别根据直接策略2a模型和组合策略2a模型计算表示变化类型的分类精度的概率值。然后，基于322个测试样本点，由式（4），（5），（6）计算精度预测值与样本真值之间的平均误差ME(mean error)、平均绝对误差MAE(mean absolute error)、均方根误差RMSE(root mean square error)，见表8。

表8 直接策略与组合策略模型2a预测精度比较Tab.8 Comparison accuracy estimation by model 2a's of the direct and composite strategies

RMSE、MAE、ME3个指标值愈小，模型预测精度愈高。由表8可看出，两策略下土地覆盖变化信息的预测精度差别不大。由于直接策略相对于组合策略计算量较小，故推荐选择直接策略。

4 结束语

本文比较了预测地表覆盖变化精度的逻辑回归方法中使用协变量的两种策略，探究两者在预测精度上是否有显著差异。直接策略直接使用待验证的地表覆盖变化地图中提取的CLASS、HOM、HET、LPS等景观指数作为协变量，组合策略以所涉及的两个时相的土地覆盖分类图的景观指数的组合（即CLASS1、CLASS2、HOM1、HOM2、HET1、HET2、LPS1、LPS2）作为协变量。以似然比检验法对两种策略可选模型进行了显著性检验，分别获得两种策略下含有最多有效协变量的模型：前者包含CLASS和HOM，后者包括CLASS1、CLASS2、HOM1和HOM2。通过独立的测试样本，计算比较了上述最优模型的ME、MAE和RMSE等评价指标，发现两策略下的模型预测能力差异较小。但是，直接基于土地覆盖变化地图进行局域精度预测建模和计算的效率更高。因此，基于运算效率的考虑，推荐在实际精度预测工作中使用直接策略。

本文的研究成果可以应用于地表覆盖变化图的局域精度估计。后续研究可探讨逻辑回归与地质统计方法的联合使用，用以深入探讨地理研究领域广泛存在的空间-时间依赖特性，精化土地覆盖变化信息精度预测方法。其次，应研究发展面向土地覆盖变化的精度验证的空间抽样方法[13]，结合地图精度-土地景观格局的耦合关系，优化抽样方案，满足精度验证的可靠性要求和特定的信息服务和应用目的。另外，由于不同传感器具有不同的空间分辨率，当像素尺寸增大或减小时，景观指数将随之变化，最优逻辑回归模型及其预测值也将不同；探究逻辑回归预测精度与所研究变化信息的尺度效应将具有重要意义。