CN111079635B - 一种基于变化向量的有林地变化检测方法 - Google Patents

Abstract

本发明公开了一种基于变化向量的有林地变化检测方法，在进行变化检测的过程中，通过在不同变化向量角度范围设置不同的阈值，能够有效去除伪变化信息，获取到更准确的有林地变化区域及其变化类型。

Description

一种基于变化向量的有林地变化检测方法

技术领域

本发明公开了一种基于变化向量的有林地变化检测方法，属于变化检测技术领域。

背景技术

近年来，随着经济的快速发展和人口的快速增加，我国的城市规模不断扩张。对建设用地迅猛增加的需求在一定程度上造成了我国森林资源的快速变化，因此及时掌握森林资源的消长变化情况能够为我国森林保护管理和生态建设提供可靠支撑。随着遥感技术的发展，基于遥感图像进行有林地变化检测已经成为快速准确获取有林地变化情况的有效手段。

遥感变化检测方法主要包括直接比较法和分类后比较法两类，分类后比较法的结果会受到遥感图像分类精度的影响且分类过程费时费力，而直接比较法运算速度快且能有效获取土地覆盖变化信息。在直接比较法中，变化向量是一种较为常用的方法，利用图像全部波段或部分波段，在得到变化区域的同时获取其变化类型信息。

目前基于变化向量的变化检测方法主要技术实现方案为：

(1)输入同一区域不同时期的两幅遥感图像，这两幅图像已经过辐射标定、大气校正、地形校正和图像配准等预处理。假定获取时间为t₁的遥感图像为X₁，获取时间为t₂的遥感图像为X₂；

(2)确定两幅图像参与变化检测的波段，此时应保证两幅图像各自参与的波段数相同且对应波段的波长范围相同。假定两幅图像参与的波段数均为B，令X_t，b表示多光谱图像X_t(t＝1,2)在第b波段的所有像元光谱值的集合，则X₁＝(X_1，1,X_1，2,…,X_1，B)^T，X₂＝(X_2，1,X_2，2,…,X_2，B)^T；

(3)通过图像X₁和X₂中的每个对应空间位置(即对应像元)的光谱向量相减获取多光谱差分图像X_D，如公式1所示，令X_D，b表示多光谱差分图像X_D在第b波段的分量，即所有像元的光谱变化向量的集合，则X_D＝(X_D，1,X_D，2,…,X_D，B)^T＝(X_2，1-X_1，1,X_2，2-X_1，2,…,X_2，B-X_1，B)^T；

X_D＝X₂-X₁ (1)

(4)令ρ表示像元光谱变化向量的强度，令θ表示像元光谱变化向量的角度。则在极坐标系中，ρ的计算如公式2所示，θ的计算如公式3所示：

(5)首先确定像元的光谱变化向量强度的阈值T，如果像元对应的ρ≥T，则该像元为变化像元，反之为未变化像元。然后在变化像元集合中分出属于不同角度范围的子像元集合，每一个子像元集合代表了一种变化类型。如图1所示，图中灰色区域为未变化像元所在区域，黑色区域中的像元代表了有着某种变化类型的变化区域。

在利用变化向量进行变化检测的过程中，所有类型的土地覆盖变化信息采用同一阈值获取，可能会造成伪变化信息的增加和真变化信息的漏检。同时，通过变化向量获取的土地覆盖变化结果中包含多种土地覆盖变化类型，要得到有林地的变化信息还需对结果进行进一步的筛选。

发明内容

本发明的目：针对上述现有技术存在的问题和不足，本发明的目的是提供一种从遥感图像中提取有林地变化信息的方法，主要解决获取具有有林地变化类型属性的变化区域的问题。

技术方案：本发明公开了一种基于变化向量的有林地变化检测方法，包括以下步骤：

步骤1：获取待检测区域不同时期的两幅遥感图像，并对该两幅遥感图像进行预处理；

步骤2：根据两幅遥感图像参与变化检测的波段，得到两幅遥感图像在该波段的所有像元的光谱向量；

步骤3：将两幅遥感图像中的每个对应像元的光谱向量相减得到各像元的变化向量；

步骤4：在极坐标系中，采用式(2)计算得到各像元的变化向量强度ρ，采用式(3)计算得到各像元的变化向量角度θ；

步骤5：根据步骤3计算得到的各像元的变化向量，设置其对应的变化强度阈值T，将变化向量强度值大于该变化强度阈值T的像元判定为候选变化像元；

步骤6：根据聚类数对步骤5得到的候选变化像元进行聚类，得到聚类结果，聚类结果所在的角度范围为变化类型对应的候选角度范围；该聚类数由待检测区域上存在的变化类型数决定；

步骤7：在每个候选角度范围下，设置其对应的变化强度阈值T’，将变化向量强度值大于该变化强度阈值T’的像元判定为最终变化像元，根据该像元所属的候选角度范围确定其变化类型。

进一步的，所述步骤5中，根据步骤3计算得到的各像元的变化向量，采用最大期望算法计算得到变化强度阈值T。

进一步的，所述步骤6中，利用k-means聚类方法对步骤5得到的候选变化像元进行聚类。

进一步的，所述步骤7中，在每个候选角度范围下，采用大津算法计算得出其对应的变化强度阈值T’。

进一步的，在步骤7之后还包括：

步骤8：逐一判断每个候选角度范围下的最终变化像元是否均匀随机分布在待检测区域上，将均匀随机分布在待检测区域上的最终变化像元去除。

有益效果：本发明在进行变化检测的过程中，通过在不同变化向量角度范围设置不同的阈值，能够有效去除伪变化信息，获取到更准确的有林地变化区域及其变化类型。

附图说明

图1为现有技术中变化检测极坐标示意图；

图2为本发明的有林地变化检测方法示意图；

图3为本实施例的实验区及各区域影像；

图4为本实施例的基于变化向量的有林地变化检测方法的各区域变化检测结果。

具体实施方式

现结合附图和实施例进一步阐述本发明的技术方案。

本发明为一种改进的变化向量分析法，旨在通过在不同变化向量角度范围设置不同的阈值，获取更准确的有林地变化区域及其变化类型，具体步骤如下：

步骤1：确定候选变化像元。首先计算出图像上各像元的变化向量，其中，变化向量的强度用ρ表示(由公式2计算得到)，变化向量的角度用θ表示(由公式3计算得到)。在极坐标系下，未变化的像元集中分布在坐标原点附近，它们的ρ相对偏小，因此通过强度阈值T可以区分变化像元和非变化像元。在这里我们假定变化像元和非变化像元均服从高斯分布，因此，图像上所有的像元服从由两个高斯分布组合而成的混合高斯分布。根据贝叶斯理论，可以采用最大期望算法计算出变化强度阈值T。当像元的ρ值大于T时，此像元被判定为候选变化像元。

步骤2：确定候选的角度范围。在极坐标系下，具有相同变化类型的变化像元有着相近的变化向量角度θ，因此根据θ的分布情况，可以得到各变化类型对应的角度范围。在本实施例中，利用k-means聚类方法对候选变化像元进行聚类，聚类数由图像上存在的变化类型数决定，得到的聚类结果所在的角度范围即为变化类型对应的候选的角度范围。

步骤3：确定最终的变化像元及其变化类型，基于步骤2中得到的候选角度范围，在每一个角度范围下，基于属于这个角度范围内的像元的变化向量强度ρ，采用大津算法计算得出变化强度阈值T’，当像元的ρ值大于T’时，此像元被判定为最终变化像元，同时根据该像元所属的角度范围确定其变化类型。如果某个角度范围下的变化像元均匀随机分布在影像上，则该角度范围下的变化像元均为伪变化像元，将其从最终变化像元中去除。

本发明在进行变化检测的过程中，通过在不同变化向量角度范围设置不同的阈值，能够有效去除伪变化信息，获取到更准确的有林地变化区域及其变化类型。

附图2中，T代表步骤1确定的强度阈值，α₁、α₂、α₃用来区分不同的角度范围，T₁、T₂、T₃和T₄是每个角度范围下的强度阈值，同时，图中四种不同的颜色代表图像上的不同变化类型。

实施例：

本实施例在安徽省黄山市和河北省承德市选取3个区域作为实验区，其中区域1位于黄山，区域2和区域3位于承德。具体情况如图3所示。

以区域1为例，简述利用本发明获得变化区域及类型的操作步骤：

(1)确定影像参与变化向量计算的波段：区域1位于安徽省黄山市，对应两幅影像的获取时间分别是2017年10月9日和2018年10月4日。根据这两幅图像各波段的波长覆盖范围，选择绿光波段(中心波长为560.0nm)、红光波段(中心波长为664.5nm)、红边1波段(中心波长为703.9nm)和近红外波段(中心波长为835.1nm)组成光谱向量，参与区域1的变化向量计算；

(2)计算区域1的变化向量：根据公式1，将区域1对应的2017年10月9日的影像和2018年10月4日的影像中的每个对应空间位置(即对应像元)的光谱向量相减获取多光谱差分图像X_D；

(3)确定区域1的候选变化像元：根据步骤(2)计算得到的变化向量，采用最大期望算法计算得到区域1的强度阈值为711，当像元的ρ值大于711时，此像元被判定为候选变化像元；

(4)确定区域1的候选角度范围：由于区域1中仅存在着有林地到裸地、建设用地这一种变化类型，因此利用k-means聚类方法对步骤(3)得到候选变化像元进行聚类，将其分为2类，由此得到两个角度范围[0°,82°)和[82°,180°]。

(5)确定区域1最终的变化像元及其变化类型：基于步骤(4)中得到的候选角度范围，分别在[0°,82°)和[82°,180°]这两个角度范围下，采用大津算法计算得出变化强度阈值T。其中，[0°,82°)下的变化强度阈值为1327，当像元的ρ值大于1327时，此像元被判定为最终变化像元，同时根据该像元所属的角度范围确定其变化类型为有林地到裸地、建设用地。[82°,180°]下的强度阈值为357，当像元的ρ值大于357时，此像元被判定为最终变化像元，但[82°,180°]下的变化像元均匀随机分布在影像上，这些变化像元均为伪变化像元，将其从最终变化像元中去除。由此得到区域最终的变化像元，它们的变化类型为有林地到裸地、建设用地。

区域2、3的操作步骤与区域1相同，最终的变化检测结果如图4所示。其中，区域1中检测出的深灰色变化区域对应有林地到裸地、建设用地，区域2中检测出的深灰色变化区域对应有林地到裸地、建设用地，浅灰色变化区域对应裸地、建设用地到水体，区域3中检测出的深灰色变化区域对应有林地到火烧、采伐迹地，浅灰色变化区域对应有林地到云，灰度介于深灰和浅灰的变化区域(已在图中标注出)对应有林地到云影。

Claims (5)

1.一种基于变化向量的有林地变化检测方法，包括以下步骤：

步骤1：获取待检测区域不同时期的两幅遥感图像，并对该两幅遥感图像进行预处理；

步骤2：根据两幅遥感图像参与变化检测的波段，得到两幅遥感图像在该波段的所有像元的光谱向量；

步骤3：将两幅遥感图像中的每个对应像元的光谱向量相减得到各像元的变化向量；

步骤4：在极坐标系中，采用式(2)计算得到各像元的变化向量强度ρ，采用式(3)计算得到各像元的变化向量角度θ；

式中，B为参与的波段数，b表示第b波段，X_D，b表示多光谱差分图像X_D在第b波段的分量，X_D＝(X_D，1，X_D，2，...，X_D，B)^T＝(X_2，1-X_1，1，X_2，2-X_1，2，...，X_2，B-X_1，B)^T，X_2，b表示图像X₂在第b波段的分量，X_1，b表示图像X₁在第b波段的分量；

其特征在于：还包括以下步骤：

步骤5：根据步骤3计算得到的各像元的变化向量，设置其对应的变化强度阈值T，将变化向量强度值大于该变化强度阈值T的像元判定为候选变化像元；

步骤6：根据聚类数对步骤5得到的候选变化像元进行聚类，得到聚类结果，聚类结果所在的角度范围为变化类型对应的候选角度范围；该聚类数由待检测区域上存在的变化类型数决定；

步骤7：在每个候选角度范围下，设置其对应的变化强度阈值T’，将变化向量强度值大于该变化强度阈值T’的像元判定为最终变化像元，根据该像元所属的候选角度范围确定其变化类型。

2.根据权利要求1所述的一种基于变化向量的有林地变化检测方法，其特征在于：所述步骤5中，根据步骤3计算得到的各像元的变化向量，采用最大期望算法计算得到变化强度阈值T。

3.根据权利要求1所述的一种基于变化向量的有林地变化检测方法，其特征在于：所述步骤6中，利用k-means聚类方法对步骤5得到的候选变化像元进行聚类。

4.根据权利要求1所述的一种基于变化向量的有林地变化检测方法，其特征在于：所述步骤7中，在每个候选角度范围下，采用大津算法计算得出其对应的变化强度阈值T’。

5.根据权利要求1所述的一种基于变化向量的有林地变化检测方法，其特征在于：在步骤7之后还包括：

步骤8：逐一判断每个候选角度范围下的最终变化像元是否均匀随机分布在待检测区域上，将均匀随机分布在待检测区域上的最终变化像元去除。

Images