CN109726679A - 一种遥感分类误差空间分布制图方法 - Google Patents

Abstract

本发明属于遥感信息不确定性可视化表达领域，涉及一种遥感分类误差空间分布制图方法，在遥感分类误差空间分布制图模型中加入最大后验概率、后验概率信息熵作为模型指标因子。这两个指数有效地表达了遥感图像像元由于同物异谱、异物同谱现象引起的分类误差的概率。通过加入两个新的光谱指数结合空间景观特征指数来构建遥感分类误差空间分布制图模型并制作遥感分类误差空间分布图。本专利提出的遥感分类误差制图模型与技术具有表达精度高，稳定性强，适用范围更广等特点。本专利的方法适用于各种景观特征下的遥感图像分类误差空间分布制图，尤其适用于需要表达分类结果中由于同物异谱、异物同谱引起的分类误差的空间分布。

Description

一种遥感分类误差空间分布制图方法

技术领域

本发明属于遥感信息不确定性可视化表达领域，特别涉及遥感分类误差空间分布可视化表达领域。

背景技术

土地覆盖等专题地图已经广泛地服务于全球或区域不同尺度的生态、农业、环境等模型等多种后续研究和应用中，为全球气候变化、农业、生态、环境等政策的制定提供了重要的依据。遥感分类是从遥感影像中获取土地覆盖类型等信息的一种重要方法。然而，由于遥感影像成像过程中的尺度效应、异物同谱、混合像元等问题以及分类方法本身的缺陷等因素导致由遥感分类得到的专题产品与实际地面相比存在一定的误差。这种误差直接或间接地增加了以遥感专题分类成果数据为数据源的后续研究结论和应用的不确定性，甚至影响了决策制定的可靠性。因此，准确、定量地表达遥感分类误差空间分布至关重要。

当前，定量表达遥感分类误差空间分布主要有以下三种方法：

(1)误差矩阵法

误差矩阵又称混淆矩阵，是一种基于随机样本用于比较分类结果与地面参考数据之间符合程度的矩阵(见表1)。矩阵的行代表遥感影像分类得到的分类结果图，列表示地面参考数据中的类型。误差矩阵中某个元素q_mn表示遥感分类图像上属于第m类且在地面参考数据中数据属于第n类的像元个数。通过构建该误差矩阵能够计算出总体精度(OA)，Kappa系数，用户精度(UA)和生产者精度(PA)等反映遥感分类准确度的指标。该方法的优势在于用简单的指标有效地表达了遥感分类图的精度，易于理解；但是不能表达遥感分类误差在空间上的分布。

(2)空间插值法

空间插值法的基本假设是遥感分类误差在空间上不是相互独立的，存在一定的空间自相关性。该方法的基本做法是首先对遥感分类图随机抽取一定量的样本，然后将样本分类结果与参考数据进行对比，正确的标记为“1”，错误的标记为“0”；然后将样本的“0-1”值插值成全图的遥感分类误差空间分布图。空间插值采用的数学方法一般有克里金方插值法、反距离加权插值法等。空间插值方法的优势在于通过样本间的相关性和简单的计算方法就得到了误差的空间分布图；劣势在于空间插值法中每一个待插值的像元须在样本的有效控制范围之内，从而造成该方法对样本的数量和在空间分布上的要求极高。

(3)函数回归法

函数回归模型用于表达遥感分类误差空间分布的理论假设是像元分类误差与像元的空间景观结构等地物分布性质高度相关。例如，地块破碎的区域在遥感影像上更容易发生错误的分类现象。基本思路是寻找合适的空间景观特征指标，然后建立这些指标与分类误差之间的函数关系式，最后对每个像元的误差概率进行计算和制图。该方法的关键在于空间景观特征指数的选择。常用的空间景观指数主要有空间异质性和图斑斑块面积等。相比空间插值法，函数回归模型降低了对样本数量的依赖。但是，目前作为回归模型的变量均来自表达像元与其邻接像元之间空间关系的景观格局指数；缺乏像元本身光谱相似性等信息导致分类误差的指标，使得最终表达的遥感分类误差空间分布图主要表达了地块边缘部分的误差，表达的精度存在一定的偏差。

综上所述，以上传统的空间插值法对样本的数量和空间分布有严格的要求，增加了野外调查的成本；函数回归法只考虑了空间景观特征方面的指标，导致表达的分类误差空间分布的一般特点是地物边缘分类误差多，斑块内部像元的误差小，难以反映遥感普遍的“异物同谱”、“同物异谱”、“分类归属”等问题导致的像元误差，致使表达的遥感分类误差空间分布与实际偏差很大。

发明内容

本发明要解决的技术问题是：找到能够反映遥感分类误差在空间上分布的像元光谱特征因子指标，并结合传统的空间景观特征因子指标；然后以这些因子指标建立遥感分类误差空间分布表达模型进行误差空间分布制图，提高遥感分类误差空间分布图表达的准确度与精确性。

本发明解决其技术问题所采用的技术方案是：

一种遥感分类误差空间分布制图方法，包含以下步骤：

步骤1)：遥感图像分类；

步骤2)：随机抽选样本并标记样本分类正确与否；

步骤3)：模型因子指标选择与归一化计算；

步骤4)：构建误差概率制图模型；

步骤5)：遥感分类误差空间分布制图；

具体地，上述步骤1)至步骤5)具体为：

步骤1)：对一副遥感影像利用机器学习算法进行分类，得到土地覆盖类型空间分布图以及每个像元属于各类型的后验概率图；

步骤2)：随机抽样本并标记样本的分类正确性。对遥感分类图像随机抽取1套独立的、一定比例(数量)的样本；选择参考影像或通过地面调查，通过比较分类图像与参考影像(或地面调查)数据，标记每个样本是否分类正确，正确的标记为“1”，错误的标记为“0”；

步骤3)：依据分类图和后验概率图，计算表达遥感分类误差分布模型的5个因子指标，主要包括3个景观格局指数(空间异质性、平均斑块大小、像元所在斑块面积)和2个光谱特征指数(最大后验概率、后验概率信息熵)；

步骤4)：构建误差概率回归模型，模型自变量为步骤3)中的5个因子指标，因变量为步骤2)中的分类是否正确的值，值为“1”或“0”，制图模型采用Logistic回归模型；利用步骤2)中的样本数据解算模型中的未知参数；

步骤5)：遥感分类误差空间分布制图，将模型应用于遥感分类图像中的每个像元，依据每个像元的5个模型指标计算每个像元的分类误差概率，然后在Arc GIS软件中进行遥感分类误差空间分布制图。

其中，步骤1)中机器学习分类算法一般选择最大似然法、SVM支持向量机、光谱角、平行六面体等算法；分类的类型数量至少在2类(含)以上。

具体地，步骤2)中的一定数量的样本，一般为总体分类像元数量的0.5％-2％，遥感分类图像中每个类型的样本总数最低不能少于120个。

具体地，步骤2)中的参考数据来源如果是地面调查，则一般采用GPS在野外采集实际地物的类型；参考数据的来源如果是高分影像，那么高分影像的空间分辨率至少小于等于遥感分类图像分辨率的1/5(例如，分类影像的空间分辨率是30米，那么高分影像的空间分辨率至少是6米)。

具体地，步骤3)中的空间异质性(HET)的计算公式如下，

其中，y_i为分类图像上第i个像元的类型，y_k表示与y_i邻接的像元的类型。当y_i与y_k类型相同时，Ψ(y_i＝y_k)取值为0；当y_i与y_k类型不同时，Ψ(y_i＝y_k)取值为1。z表示邻域关系，四邻域时z＝4，八邻域时z＝8。HET值代表了像元类别与周围像元类别之间的差异性。H值越高表示像元的空间结构越复杂，分类的误差概率越高。

具体地，步骤3)中像元所在斑块面积的计算公式如下，

A_j＝(PatchArea_j|i∈Patch_j)

其中，i∈Patch_j表示图像上第i个像元属于斑块j，判定第i个像元属于斑块j的原则是第i个像元的分类类型与斑块j(Patch_j)中每一个像元的类型相同，并且与斑块j中至少有一个像元相邻。统计第j个斑块的面积为PatchArea_j。

具体地，步骤3)中平均斑块大小计算公式如下，

其中，PatchArea_i是第i个斑块的面积，所有斑块的面积的总和，n为所有斑块的数量。

具体地，步骤3)中最大后验概率(P_max(k))计算的方法，如下：

P_max(k)＝Max(P(ω_i|x_k))

其中，x_k是遥感图像上第k个像元的光谱向量；P(ω_i|x_k)是经过机器学习遥感分类后x_k属于第ω_i类的后验概率；Max(P(ω_i|x_k))是x_k属于不同类的后验概率中的最大值。最大后验概率越大，属于第ω_i类的可能性越高，在分类过程中发生错误的可能性越低。

具体地，步骤3)中后验概率信息熵1₂(k)的计算方法，如下：

其中，c代表遥感图像经过机器学习分类后的类型数量，ω_i是c个类型中的第i个类型；x_k是遥感图像上第k个像元的光谱向量；P(ω_i|x_k)是经过机器学习分类后x_k属于第ω_i类的后验概率；l6gP(ω_i|x_k)是P(ω_i|x_k)的对数值，H₂(k)是后验概率信息熵。如果H₂(k)的值较大,则说明各类的后验概率较接近，在类别归属过程发生错误判断的可能性较高。当各类后验概率相等时，实现正确分类识别是最困难的，此时描述可分性的信息熵值最大。

具体地，步骤3)中的归一化公式如下：

其中，X_i ^*为某个因子指标归一化后的值，取值范围为[0,1]，X_i为某个因子指标归一化前的值，X_min为某个因子指标的最小值，X_max为某个因子指标的最大值。

具体地，步骤4)中图像上任意一像元发生错误率的可能性与上述模型参数的定量模型如下：

公式中，p_error(c)为自变量，代表单像元的错误率，在建立Logistics回归模型时为步骤2)中的分类是否正确的值，值为“1”或“0”。x_k(c)为回归变量，当K＝5时，x_k(c)分别是空间异质性(HET)、像元所在斑块面积(A_j)、最大后验概率(P_max(k))、平均斑块大小(A_average)、最大后验概率(P_max(k))和后验概率信息熵H_c(k)。β_k是回归系数，其获取的方式是通过步骤2)中样本数据利用统计软件(例如，SPSS、SAS等)进行回归解算。这些样本中已知数据为自变量(9_error(c))和回归变量(空间异质性(HET)、像元所在斑块面积(A_j)、最大后验概率(P_max(k))、平均斑块大小(A_average)、最大后验概率(P_max(k))和后验概率信息熵H₂(k))。

本发明的有益效果主要有以下2个方面：

本发明提供了一种遥感分类误差空间分布制图方法，本发明的方法提高了分类误差空间分布表达的精确性，尤其是在遥感分类斑块内部像元误差表达精度方面有显著的提升。试验结果表明：在景观相对破碎的情况下，本发明技术生产的遥感分类误差分布图精细地描述了破碎地块边缘和地块内部的误差空间分布图。在景观相对破碎的情况下，本发明的评估值(s)为0.21，函数回归法的评估值为0.45，空间插值法的评估值为0.64。在景观相对规整的情况下，本发明的评估值(s)为0.27，函数回归法的评估值为0.39，空间插值法的评估值为0.61。评估值越小表明，遥感分类误差空间分布与实际的分布差异越小。这表明本专利生产的遥感分类误差空间分布图在景观相对破碎、相对规整的情况下，好于目前普遍采用的空间插值方法和函数回归法。同时相比传统的空间插值方法，该方法对野外样本的需求量较少，极大地减少了实际工作量，丰富了遥感精度验证关于误差表达方式的技术体系。

本发明方法主要针对遥感分类误差在空间分布上表达模型的指标因子选择方面进行了创新，以此能够增强制作遥感分类误差空间分布图的精确性。本发明的核心创新点在于提出了从景观特征和光谱特征两个角度选择分类误差的影响因子；建立单像元分类误差与景观格局指数(空间异质性、平均斑块大小、像元所在斑块面积)和光谱特征指数(最大后验概率、后验概率信息熵)之间的函数模型，并且对比分析了加入光谱特征前后误差空间分布表达的精度。在实际应用中，主要应用于土地覆盖/利用、农业遥感分类后野外样本抽选、农作物面积估算、全球气候模型中关于不确定性研究等方面的输入数据。

附图说明

下面结合附图和实施例对本发明进一步说明。

图1为本发明总体技术路线图。本发明提供的一种遥感分类误差空间分布制图方法的总体步骤是遥感分类、随机抽选样本与正确度标记、模型因子计算、分类误差概率模型构建、遥感分类误差空间分布制图、模型评估。

图2(a)为北京市大兴地区2011年7月26日的一景Landsat TM 5遥感影像(空间分辨率：30米)，空间范围为22.25km²，图像中心点坐标为116°21'55.64"E39°40'45.034"N。

图2(b)为在进行标准化预处理(投影转换、辐射定标、大气校正等)之后，采用最大似然法分类后的土地覆盖分类图。主要分为裸地、水体、居民地、耕地、其他植被5类。

图2(c)为遥感分类误差进行误差空间分布制图的野外样本分布图。

图2(d)为大兴地区土地覆盖遥感分类误差空间分布图。

图3(a)为北京市通州地区2011年7月26日的一景Landsat TM 5遥感影像(空间分辨率：30米)，空间范围为22.25km²，图像中心点坐标为116°42'19.91"E39°45'46.967"N。

图3(b)为在进行标准化预处理(投影转换、辐射定标、大气校正等)之后，采用最大似然法分类后的土地覆盖分类图。主要分为裸地、水体、居民地、耕地、其他植被5类。

图3(c)为遥感分类误差进行误差空间分布制图的野外样本分布图。

图3(d)为通州地区土地覆盖遥感分类误差空间分布图。

具体实施方式

现在结合附图对本发明作进一步详细的说明。这些附图均为简化的示意图，仅以示意方式说明本发明的基本结构，因此其仅显示与本发明有关的构成。

实施案例1：应用于景观相对破碎区域土地覆盖遥感分类误差空间分布制图场景，试验区域以北京市大兴地区以为例，采用简单随机抽样方法抽选的样本量为800个，平均每个类别至少有150个。模型指标中的空间异质性采用8邻域计算方法。建立的模型的表达式为：采用以下评估公式对遥感分类误差空间分布图进行评估：

公式中，Y_i为样本中第i样本标记的分类正确与否值(“1”或“0”)，Y_Pi为第i样本对应的步骤5)计算的遥感分类误差概率值。s值越小代表遥感分类误差制图模型精度越高。一般来说，s值小于0.4则表示模型精度高；s值在0.4-0.6表示模型精度一般；s值大于0.6则表示模型精度低。

试验结果表明：在景观相对破碎的情况下(图2(a))，本发明技术生产的遥感分类误差分布图精细地描述了破碎地块边缘和地块内部的误差空间分布图(图2(d))。表1表明，在景观相对破碎的情况下，本发明的评估值(s)为0.21，函数回归法的评估值为0.45，空间插值法的评估值为0.64。评估值越小表明，遥感分类误差空间分布与实际的分布越相近。这表明本专利生产的遥感分类误差空间分布图在景观相对破碎的情况下，好于目前普遍采用的空间插值方法和函数回归法。

实施例2：应用于景观像对规整区域土地覆盖遥感分类误差空间分布制图场景，试验区以北京市通州地区为例。采用简单随机抽样方法抽选的样本量为800个，平均每个类别至少有150个。模型指标中的空间异质性采用8邻域计算方法。建立的模型的表达式为：，采用以下评估公式对遥感分类误差空间分布图进行评估：

公式中，Y_i为样本中第i样本标记的分类正确与否值(“1”或“0”)，Y_Pi为第i样本对应的步骤4)计算的错误概率值，s值小于0.4则表示模型精度高；s值在0.4-0.6表示模型精度一般；s值大于0.6则表示模型精度低。

试验结果表明：在景观相对规整的情况下(图3(a))，本发明技术生产的遥感分类误差分布图精细地描述了规整地块边缘和地块内部的误差空间分布图(图3(d))。表1表明，在景观相对规整的情况下，本发明的评估值(s)为0.27，函数回归法的评估值为0.39，空间插值法的评估值为0.61。评估值越小表明，遥感分类误差空间分布与实际的分布越相近。这表明本专利生产的遥感分类误差空间分布图在景观相对规整的情况下，好于目前普遍采用的空间插值方法和函数回归法。

表1本专利技术和空间插值方法对比

Claims (10)

1.一种遥感分类误差空间分布制图方法，包含以下步骤：

步骤1)：遥感图像分类；

步骤2)：随机抽选样本并标记样本分类正确与否；

步骤3)：模型因子指标选择与归一化计算；

步骤4)：构建误差概率制图模型；

步骤5)：遥感分类误差空间分布制图。

2.如权利要求1所述的一种遥感分类误差空间分布制图方法，其特征在于：

所述步骤1)至步骤5)具体为：

步骤1)：对一副遥感影像利用机器学习算法进行分类，得到土地覆盖类型空间分布图以及每个像元属于各类型的后验概率图；

步骤2)：对遥感分类图像随机抽取1套独立的、一定比例的样本；选择参考影像或通过地面调查，通过比较分类图像与参考影像或地面调查数据，标记每个样本是否分类正确，正确的标记为“1”，错误的标记为“0”；

步骤3)：依据分类图和后验概率图，计算表达遥感分类误差分布模型的5个因子指标，包括空间异质性、平均斑块大小、像元所在斑块面积、最大后验概率和后验概率信息熵；然后按照归一化公式将所有指标的值归至[0,1]范围内；

步骤4)：构建误差概率回归模型，模型自变量为步骤3)中的5个因子指标，因变量为步骤2)中的分类是否正确的值，值为“1”或“0”，制图模型采用Logistic回归模型；利用步骤2)中的样本数据解算模型中的未知参数；

步骤5)：遥感分类误差空间分布制图，将模型应用于遥感分类图像中的每个像元，依据每个像元的5个模型指标计算每个像元的分类误差概率，然后在Arc GIS软件中进行遥感分类误差空间分布制图。

3.如权利要求2所述的一种遥感分类误差空间分布制图方法，其特征在于：所述步骤2)中一定比例的样本为总体分类像元数量的0.5％-2％，并且遥感分类图像中每个类型的样本数不少于120个。

4.如权利要求2所述的一种遥感分类误差空间分布制图方法，其特征在于：步骤3)中最大后验概率P_max(k)的计算方法如下：

P_max(k)＝Max(P(ω_i|x_k))

其中，x_k是遥感图像上第k个像元的光谱向量；P(ω_i|x_k)是经过机器学习遥感分类后x_k属于第ω_i类的后验概率；Max(P(ω_i|x_k))是x_k属于不同类的后验概率中的最大值。

5.如权利要求2所述的一种遥感分类误差空间分布制图方法，其特征在于：步骤3)中后验概率信息熵H_c(k)的计算方法如下：

其中，c代表遥感图像经过机器学习分类后的类型数量，ω_i是c个类型中的第i个类型；x_k是遥感图像上第k个像元的光谱向量；P(ω_i|x_k)是经过机器学习分类后x_k属于第ω_i类的后验概率；logP(ω_i|x_k)是P(ω_i|x_k)的对数值。

6.如权利要求2所述的一种遥感分类误差空间分布制图方法，其特征在于：步骤3)中的空间异质性HET的计算公式如下：

其中，y_i为分类图像上第i个像元的类型，y_k表示与y_i邻接的像元的类型，当y_i与y_k类型相同时，Ψ(y_i＝y_k)取值为0；当y_i与y_k类型不同时，Ψ(y_i＝y_k)取值为1；z表示邻域关系，四邻域时z＝4，八邻域时z＝8。

7.如权利要求2所述的一种遥感分类误差空间分布制图方法，其特征在于：步骤3)中像元所在斑块面积的计算公式如下，

A_j＝(PatchArea_j|i∈Patch_j)

其中，i∈Patch_j表示图像上第i个像元属于斑块j，判定第i个像元属于斑块j的原则是第i个像元的分类类型与斑块j(Patch_j)中每一个像元的类型相同，并且与斑块j中至少有一个像元相邻，统计第j个斑块的面积为PatchArea_j。

8.如权利要求2所述的一种遥感分类误差空间分布制图方法，其特征在于：步骤3)中平均斑块大小计算公式如下，

其中，PatchArea_i是第i个斑块的面积，所有斑块的面积的总和，n为所有斑块的数量。

9.如权利要求2所述的一种遥感分类误差空间分布制图方法，其特征在于：步骤4)中，误差概率回归模型如下：

其中，A_k为回归系数，x_k(c)为自变量，p_9rror(c)为表达的单像元的错误率。

10.如权利要求2所述的一种遥感分类误差空间分布制图方法，其特征在于：步骤2)中参考影像的空间分辨率高于分类图像的空间分辨率。