CN114283335A - 一种历史时期遥感识别精度验证制备方法 - Google Patents

Abstract

本发明提出了一种历史时期遥感识别精度验证制备方法，其步骤为：根据地表目标遥感识别结果分布图布，采用随机抽样方法布设样本框；对于每个样本框，基于空间分辨率优于1米的高清卫星影像对地物边界进行矢量化，并结合历史资料、影像特征、识别结果分布图综合解译确认每个地物矢量多边形的类型属性；将解译得到的地物矢量多边形数据转成栅格数据，其空间分辨率与识别结果分布图的空间分辨率一致；最后采用混淆矩阵精度验证方法开展精度验证，并使用精度数据对目标地物面积数据进行修正。本发明提出的遥感识别精度验证制备流程，具有操作简便、高效，精度数据客观、科学的优势，解决了历史时期精度验证样本无法实地调研获取的难题。

Description

一种历史时期遥感识别精度验证制备方法

技术领域

本发明涉及遥感目标识别技术领域，特别是指一种历史时期遥感识别精度验证制备方法。

背景技术

历史时期地表目标遥感识别是遥感目标识别技术领域的常见应用。受历史实地调查无法逆转再现的影响，准确获得历史时期真实地表目标的空间分布成为遥感目标识别结果精度验证的瓶颈之一。此外，现行技术领域缺少精度验证的技术规范，大多采用随机抽取像元点的方式开展精度验证，易造成精度验证客观性、代表性、科学性存疑。

发明内容

针对上述技术背景中存在的不足，本发明提出了一种历史时期遥感识别精度验证制备方法，解决了历史时期真实地表目标的空间分布数据获取困难的技术问题。

本发明的技术方案是这样实现的：

一种历史时期遥感识别精度验证制备方法，其步骤如下：

S1、根据研究区内历史时期X年的遥感识别结果设计a×a的格网覆盖研究区，然后采用随机抽样方法确定M个验证样本框；

S2、根据X年遥感识别结果中的像元分布对步骤S1中M个验证样本框的边界进行修正，以确保验证样本框的边界与X年遥感识别结果中像元边界重叠；

S3、采集第i个验证样本框内的空间分辨率优于1米的X年的高清卫星影像，利用目视解译方法对步骤S2中修正后的验证样本框内的各种地物进行矢量化，得到地物矢量多边形，并且将矢量误差控制在1米以内；

S4、结合历史资料、影像特征、遥感识别结果信息综合解译确认第i个验证样本框内的地物矢量多边形的类别属性；

S5、将步骤S4中得到的带有类别属性的第i个验证样本框的矢量数据在ArcGIS软件中转换成栅格数据，并设置栅格数据的空间分辨率与遥感识别结果的空间分辨率一致，得到地表真实验证样本数据；

S6、循环执行步骤S3至步骤S5，得到M个地表真实验证样本数据；

S7、利用步骤S2中修正后的M个验证样本框对X年遥感识别结果进行裁剪，得到M个参与精度验证的识别结果样本数据；

S8、将步骤S6中得到的地表真实验证样本数据和S7中得到的识别结果样本数据进行一一比对，构造混淆矩阵，计算识别精度；

S9、统计遥感识别结果中研究目标的面积数据，并基于步骤S8中得到的识别精度对面积数据进行修正。

所述步骤S1中M个验证样本框的确定方法为：验证样本框的总面积占研究区总面积的0.1％～1％，验证样本框的边长大于遥感识别结果中像元的边长的30倍，采用随机抽样方法确定M个验证样本框。

所述根据X年遥感识别结果中的像元分布对步骤S1中M个验证样本框的边界进行修正的方法为：对于第i个验证样本框的四条边界，如果其第j条边界与X年遥感识别结果中像元边界不重叠，即第j条边界穿过像元，则以与第j条边界平行且距离第j条边界最近的像元边界为基准，平移第j条边界，使之与基准像元边界重叠，完成第j条边界的修正；如果第j条边界与X年遥感识别结果中像元边界重叠，则不再修正第j条边界；依此完成M个验证样本框边界的修正。

所述步骤S8中的混淆矩阵的构造方法为：对于第i组地表真实验证样本数据和识别结果样本数据，以地表真实验证样本数据为基准，对比识别结果样本数据，构造第i组地表真实验证样本数据的像元混淆矩阵，如表1所示；表1中x表示识别结果样本数据中的类别，x＝1,2,┄,m，y表示地表真实验证样本数据中的类别，y＝1,2,┄,m，m是地表真实验证样本数据和识别结果样本数据中的不同类别的总数；N_xy,i表示第i组的识别结果样本数据中的类别x识别为地表真实验证样本数据中的类别y的像元个数，N_i表示第i组的地表真实验证样本数据中的像元总个数；

表1第i组地表真实验证样本数据的像元混淆矩阵

将M组地表真实验证样本数据和识别结果样本数据的像元混淆矩阵进行累加，即：N_xy＝N_xy,1+N_xy,2+┄+N_xy,i，i＝1,2,┄,M，x＝1,2,┄,m，y＝1,2,┄,m，如表2所示；N_xy表示所有识别结果样本数据中的类别x识别为地表真实验证样本数据中的类别y的像元个数，N表示所有地表真实验证样本数据中的像元总个数；

表2 M组地表真实验证样本数据的总体像元混淆矩阵

根据公式1和2以及表2计算得到类别比例估算混淆矩阵，如表3所示；

P_xy＝W_x×(N_xy/N_x·) (1)；

W_x＝N_·x/N (2)；

式中，P_xy表示识别结果样本数据中的类别x识别为地表真实验证样本数据中的类别y的像元个数占识别结果样本数据中的类别x像元总个数的比例，其中W_x表示地表真实验证样本数据中的类别x的像元个数占总像元个数的比例，N_x·表示识别结果样本数据中的类别x的像元总个数，N_.x表示地表真实验证样本数据中的类别x的像元总个数；

表3类别比例估算混淆矩阵

据公式3计算第k个类别的使用者精度，根据公式4计算第k个类别的生产者精度，根据公式5计算总体识别精度，根据公式6计算kappa系数；

UA_k＝P_kk/P_k· (3)；

PA_k＝P_kk/P_.k (4)；

其中，UA_k为第k个类别的使用者精度，PA_k为第k个类别的生产者精度，OA为总体识别精度，kappa为kappa系数，k＝1,2,┄,m。

所述统计遥感识别结果中研究目标的面积数据的方法为：首先将X年遥感识别结果的栅格数据转换为矢量数据，然后将矢量数据的地理投影转换为等积投影，最后基于ArcGIS软件的统计功能计算研究目标的面积。

所述基于步骤S8中得到的识别精度对面积数据进行修正的方法为：

A′_k＝A_k×UA_k/PA_k (7)；

其中，A_k为第k类目标类别的识别面积，A′_k为修正后的面积。

与现有技术相比，本发明产生的有益效果为：

1)本发明提出了一种历史时期真实地表目标样本的空间分布数据的制备方法，为遥感识别精度验证技术领域提供了理论基础；

2)本发明提出了遥感识别精度验证的技术规范体系，保障了遥感识别精度的客观性、科学性。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。

图1为本发明的流程图；

图2为本发明的地表真实验证样本数据分布图。

具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有付出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

如图1所示，本发明实施例提供了一种历史时期遥感识别精度验证制备方法，其步骤如下：

S1、根据研究区内历史时期X年的遥感识别结果设计a km×a km的格网覆盖研究区，然后采用随机抽样方法确定M个验证样本框；验证样本框的总面积占研究区总面积的0.1％～1％，验证样本框的边长大于遥感识别结果中像元的边长的30倍，采用随机抽样方法确定M个验证样本框。当验证样本框的数量M大于100时，可以适当提高样本框的边长以减少样本框的数量。研究区内2019年的遥感识别结果的空间分辨率为30米，根据验证样本框的边长应大于遥感识别结果中像元的边长的30倍的原则，设计1km×1km的格网覆盖研究区；根据验证样本框总面积占研究区总面积的0.1％～1％的原则，采用随机抽样方法确定50个验证样本框的空间分布。

S2、根据X年遥感识别结果中的像元分布对步骤S1中M个验证样本框的边界进行修正，以确保验证样本框的边界与X年遥感识别结果中像元边界重叠；根据2019年遥感识别结果中像元分布状况修正步骤S1中得到的50个验证样本框的具体边界，以确保验证样本框的边界与2019年遥感识别结果中像元边界重叠，修正后的样本框的大小为30像元×33像元；具体修正方案为：对于第i个验证样本框的四条边界，如果其第j条边界与2019年遥感识别结果中像元边界不重叠，即第j条边界穿过像元，则以与第j条边界平行且距离第j条边界最近的像元边界为基准，平移第j条边界，使之与基准像元边界重叠，完成第j条边界的修正；如果第j条边界与2019年遥感识别结果中像元边界重叠，则不再修正第j条边界；依此完成50个验证样本框边界的修正。

S3、采集第i个验证样本框内的空间分辨率优于1米的2019年的高清卫星影像，利用目视解译方法对步骤S2中修正后的验证样本框内的各种地物进行矢量化，得到地物矢量多边形，并且将矢量误差控制在1米以内；

S4、结合历史资料、影像特征、遥感识别结果信息综合解译确认第i个验证样本框内的地物矢量多边形的类别属性；

S5、将步骤S4中得到的带有类别属性的第i个验证样本框的矢量数据在ArcGIS软件中转换成栅格数据，并设置栅格数据的空间分辨率与遥感识别结果的空间分辨率一致，即均为30米，得到地表真实验证样本数据；样本内的像元个数为30×33个像元。

S6、循环执行步骤S3至步骤S5，得到50个地表真实验证样本数据；

S7、利用步骤S2中修正后的50个验证样本框对2019年遥感识别结果进行裁剪，得到50个参与精度验证的识别结果样本数据；

S8、将步骤S6中得到的地表真实验证样本数据和S7中得到的识别结果样本数据进行一一比对，构造混淆矩阵，计算识别精度；

混淆矩阵的构造方法为：

对于第i组地表真实验证样本数据和识别结果样本数据，以地表真实验证样本数据为基准，对比识别结果样本数据，构造第i组地表真实验证样本数据的像元混淆矩阵，如表1所示；表1中x表示识别结果样本数据中的类别，x＝1,2,┄,m，y表示地表真实验证样本数据中的类别，y＝1,2,┄,m，m是地表真实验证样本数据和识别结果样本数据中的不同类别的总数；N_xy,i表示第i组的识别结果样本数据中的类别x识别为地表真实验证样本数据中的类别y的像元个数，N_i表示第i组的地表真实验证样本数据中的像元总个数。

表1第i组地表真实验证样本数据的像元混淆矩阵

将M组地表真实验证样本数据和识别结果样本数据的像元混淆矩阵进行累加，即：N_xy＝N_xy,1+N_xy,2+┄+N_xy,i，i＝1,2,┄,M，x＝1,2,┄,m，y＝1,2,┄,m，如表2所示；N_xy表示所有识别结果样本数据中的类别x识别为地表真实验证样本数据中的类别y的像元个数，N表示所有地表真实验证样本数据中的像元总个数。

表2 M组地表真实验证样本数据的总体像元混淆矩阵

根据公式1和2以及表2计算得到类别比例估算混淆矩阵，如表3所示。

P_xy＝W_x×(N_xy/N_x.) (1)；

W_x＝N_.x/N (2)；

式中，P_xy表示识别结果样本数据中的类别x识别为地表真实验证样本数据中的类别y的像元个数占识别结果样本数据中的类别x像元总个数的比例，其中W_x表示地表真实验证样本数据中的类别x的像元个数占总像元个数的比例，N_x·表示识别结果样本数据中的类别x的像元总个数，N_·x表示地表真实验证样本数据中的类别x的像元总个数。

表3类别比例估算混淆矩阵

据公式3计算第k个类别的使用者精度(User Accuracy，UA_k)，根据公式4计算第k个类别的生产者精度(Production Accuracy，PA_k)，根据公式5计算总体识别精度(OverallAccuracy，OA)，根据公式6计算kappa系数。

UA_k＝P_kk/P_k· (3)；

PA_k＝P_kk/P_·k (4)；

其中，UA_k为第k个类别的使用者精度，PA_k为第k个类别的生产者精度，OA为总体识别精度，kappa为kappa系数，k＝1,2,┄,m。

具体实例中，对于第1组验证样本和结果样本，以验证样本为基准，对比结果样本，构造第1组验证样本的像元混淆矩阵，如表4所示。表4中x表示结果样本中的类别，x＝1,2,3，y表示验证样本中的类别，y＝1,2,3，其中，本次实施例中的识别结果中共包括三类地物类型，即第1类地物是冬小麦，第2类地物是冬油菜，第3类地物是其他地物。N_xy,1表示第1组的结果样本中的类别x识别为验证样本中的类别y的像元个数。N₁表示第1组的验证样本中的像元总个数。其中，N_11,1＝149，N_12,1＝10，N_13,1＝2，N_21,1＝11，N_22,1＝762，N_23,1＝4，N_31,1＝1，N_32,1＝2，N_33,1＝49，N₁·,₁＝161，N₂·,₁＝777，N₃·,₁＝52，N·_1,1＝161，N·_2,1＝774，N·_3,1＝55，N₁＝990。

表4第1组验证样本的像元混淆矩阵

依次得到50组验证样本的像元混淆矩阵，然后将这50组混淆矩阵进行累加，得到所有验证样本的总体像元混淆矩阵，即：N_xy＝N_xy,1+N_xy,2+┄+N_xy,i，i＝1,2,┄,50，x＝1,2,3，y＝1,2,3，如表5所示。N_xy表示所有结果样本中的类别x识别为验证样本中的类别y的像元个数。N表示所有验证样本中的像元总个数。其中，N₁₁＝16976，N₁₂＝302，N₁₃＝280，N₂₁＝279，N₂₂＝13476，N₂₃＝55，N₃₁＝411，N₃₂＝99，N₃₃＝17622，N_1·＝17558，N_2·＝13810，N_3·＝18132，N_·1＝17666，N_·2＝13877，N_·3＝17957，N＝49500。

表5 50组验证样本的总体像元混淆矩阵

根据公式1和2以及表5计算得到类别比例估算混淆矩阵，如表6所示。本次实施例中W₁＝0.3569，W₂＝0.2803，W₃＝0.3628，P₁₁＝0.3451，P₁₂＝0.0061，P₁₃＝0.0057，P₂₁＝0.0057，P₂₂＝0.2736，P₂₃＝0.0011，P₃₁＝0.0082，P₃₂＝0.0020，P₃₃＝0.3526，P_1·＝0.3569，P_2.＝0.2803，P_3·＝0.3628，P_.1＝0.3589，P_·2＝0.2817，P_·3＝0.3594。

表6类别比例估算混淆矩阵

本次实施例中UA₁＝0.9669，UA₂＝0.9758，UA₃＝0.9719，PA₁＝0.9613，PA₂＝0.9712，PA₃＝0.9811，OA＝0.9712，kappa＝0.9565。

S9、统计遥感识别结果中研究目标的面积数据，并基于步骤S8中得到的识别精度对面积数据进行修正。首先将X年遥感识别结果的栅格数据转换为矢量数据，然后将矢量数据的地理投影转换为等积投影，最后基于ArcGIS软件的统计功能计算研究目标的面积。

所述基于步骤S8中得到的识别精度对面积数据进行修正的方法为：

A′_k＝A_k×UA_k/PA_k (7)；

其中，A_k为第k类目标类别的识别面积，A′_k为修正后的面积。

基于ArcGIS软件的统计功能计算冬小麦、冬油菜的面积分别为7062.44km²、3755.71km²；对冬小麦、冬油菜的面积数据进行修正，修正后的冬小麦、冬油菜的面积分别为7103.14km²、3773.65km²。

为了验证本发明的效果，本发明以2019年湖南省安乡县部分区域内的冬小麦和冬油菜遥感识别结果为实验对象，开展了其识别结果的精度验证实验。

以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (6)

1.一种历史时期遥感识别精度验证制备方法，其特征在于，其步骤如下：

S1、根据研究区内历史时期X年的遥感识别结果设计a×a的格网覆盖研究区，然后采用随机抽样方法确定M个验证样本框；

S2、根据X年遥感识别结果中的像元分布对步骤S1中M个验证样本框的边界进行修正，以确保验证样本框的边界与X年遥感识别结果中像元边界重叠；

S3、采集第i个验证样本框内的空间分辨率优于1米的X年的高清卫星影像，利用目视解译方法对步骤S2中修正后的验证样本框内的各种地物进行矢量化，得到地物矢量多边形，并且将矢量误差控制在1米以内；

S4、结合历史资料、影像特征、遥感识别结果信息综合解译确认第i个验证样本框内的地物矢量多边形的类别属性；

S5、将步骤S4中得到的带有类别属性的第i个验证样本框的矢量数据在ArcGIS软件中转换成栅格数据，并设置栅格数据的空间分辨率与遥感识别结果的空间分辨率一致，得到地表真实验证样本数据；

S6、循环执行步骤S3至步骤S5，得到M个地表真实验证样本数据；

S7、利用步骤S2中修正后的M个验证样本框对X年遥感识别结果进行裁剪，得到M个参与精度验证的识别结果样本数据；

S8、将步骤S6中得到的地表真实验证样本数据和S7中得到的识别结果样本数据进行一一比对，构造混淆矩阵，计算识别精度；

S9、统计遥感识别结果中研究目标的面积数据，并基于步骤S8中得到的识别精度对面积数据进行修正。

2.根据权利要求1所述的历史时期遥感识别精度验证制备方法，其特征在于，所述步骤S1中M个验证样本框的确定方法为：验证样本框的总面积占研究区总面积的0.1％～1％，验证样本框的边长大于遥感识别结果中像元的边长的30倍，采用随机抽样方法确定M个验证样本框。

3.根据权利要求1或2所述的历史时期遥感识别精度验证制备方法，其特征在于，所述根据X年遥感识别结果中的像元分布对步骤S1中M个验证样本框的边界进行修正的方法为：对于第i个验证样本框的四条边界，如果其第j条边界与X年遥感识别结果中像元边界不重叠，即第j条边界穿过像元，则以与第j条边界平行且距离第j条边界最近的像元边界为基准，平移第j条边界，使之与基准像元边界重叠，完成第j条边界的修正；如果第j条边界与X年遥感识别结果中像元边界重叠，则不再修正第j条边界；依此完成M个验证样本框边界的修正。

4.根据权利要求1所述的历史时期遥感识别精度验证制备方法，其特征在于，所述步骤S8中的混淆矩阵的构造方法为：对于第i组地表真实验证样本数据和识别结果样本数据，以地表真实验证样本数据为基准，对比识别结果样本数据，构造第i组地表真实验证样本数据的像元混淆矩阵，如表1所示；表1中x表示识别结果样本数据中的类别，x＝1,2,┄,m，y表示地表真实验证样本数据中的类别，y＝1,2,┄,m，m是地表真实验证样本数据和识别结果样本数据中的不同类别的总数；N_xy,i表示第i组的识别结果样本数据中的类别x识别为地表真实验证样本数据中的类别y的像元个数，N_i表示第i组的地表真实验证样本数据中的像元总个数；

表1 第i组地表真实验证样本数据的像元混淆矩阵

将M组地表真实验证样本数据和识别结果样本数据的像元混淆矩阵进行累加，即：N_xy＝N_xy,1+N_xy,2+┄+N_xy,i，i＝1,2,┄,M，x＝1,2,┄,m，y＝1,2,┄,m，如表2所示；N_xy表示所有识别结果样本数据中的类别x识别为地表真实验证样本数据中的类别y的像元个数，N表示所有地表真实验证样本数据中的像元总个数；

表2 M组地表真实验证样本数据的总体像元混淆矩阵

根据公式1和2以及表2计算得到类别比例估算混淆矩阵，如表3所示；

P_xy＝W_x×(N_xy/N_x) (1)；

W_x＝N_x/N (2)；

式中，P_xy表示识别结果样本数据中的类别x识别为地表真实验证样本数据中的类别y的像元个数占识别结果样本数据中的类别x像元总个数的比例，其中W_x表示地表真实验证样本数据中的类别x的像元个数占总像元个数的比例，N_x·表示识别结果样本数据中的类别x的像元总个数，N_·x表示地表真实验证样本数据中的类别x的像元总个数；

表3 类别比例估算混淆矩阵

据公式3计算第k个类别的使用者精度，根据公式4计算第k个类别的生产者精度，根据公式5计算总体识别精度，根据公式6计算kappa系数；

UA_k＝P_kk/P_k· (3)；

PA_k＝P_kk/P_·k (4)；

其中，UA_k为第k个类别的使用者精度，PA_k为第k个类别的生产者精度，OA为总体识别精度，kappa为kappa系数，k＝1,2,┄,m。

5.根据权利要求4所述的历史时期遥感识别精度验证制备方法，其特征在于，所述统计遥感识别结果中研究目标的面积数据的方法为：首先将X年遥感识别结果的栅格数据转换为矢量数据，然后将矢量数据的地理投影转换为等积投影，最后基于ArcGIS软件的统计功能计算研究目标的面积。

6.根据权利要求4或5所述的历史时期遥感识别精度验证制备方法，其特征在于，所述基于步骤S8中得到的识别精度对面积数据进行修正的方法为：

A′_k＝A_k×UA_k/PA_k (7)；

其中，A_k为第k类目标类别的识别面积，A′_k为修正后的面积。

Images