CN114254707A - 一种基于GlobeLand30的历史地表覆盖快速重建方法 - Google Patents
一种基于GlobeLand30的历史地表覆盖快速重建方法 Download PDFInfo
- Publication number
- CN114254707A CN114254707A CN202111580258.XA CN202111580258A CN114254707A CN 114254707 A CN114254707 A CN 114254707A CN 202111580258 A CN202111580258 A CN 202111580258A CN 114254707 A CN114254707 A CN 114254707A
- Authority
- CN
- China
- Prior art keywords
- surface coverage
- earth surface
- clustering
- image
- obtaining
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Pending
Links
- 238000000034 method Methods 0.000 title claims abstract description 56
- 238000004422 calculation algorithm Methods 0.000 claims abstract description 42
- 239000000203 mixture Substances 0.000 claims abstract description 32
- 238000012549 training Methods 0.000 claims abstract description 30
- 238000009826 distribution Methods 0.000 claims abstract description 24
- 230000003595 spectral effect Effects 0.000 claims abstract description 17
- 238000004364 calculation method Methods 0.000 claims abstract description 15
- 230000008859 change Effects 0.000 claims abstract description 14
- 238000005457 optimization Methods 0.000 claims abstract description 14
- 238000001514 detection method Methods 0.000 claims abstract description 9
- 238000007637 random forest analysis Methods 0.000 claims abstract description 8
- 230000000877 morphologic effect Effects 0.000 claims abstract description 7
- 230000008569 process Effects 0.000 claims abstract description 6
- 239000013598 vector Substances 0.000 claims description 46
- 238000001228 spectrum Methods 0.000 claims description 19
- 230000006870 function Effects 0.000 claims description 14
- XLYOFNOQVPJJNP-UHFFFAOYSA-N water Substances O XLYOFNOQVPJJNP-UHFFFAOYSA-N 0.000 claims description 11
- 238000002310 reflectometry Methods 0.000 claims description 10
- 238000005070 sampling Methods 0.000 claims description 8
- 238000000354 decomposition reaction Methods 0.000 claims description 5
- 238000002156 mixing Methods 0.000 claims description 3
- 238000010606 normalization Methods 0.000 claims description 2
- 230000011218 segmentation Effects 0.000 claims description 2
- 238000003066 decision tree Methods 0.000 claims 1
- 238000012804 iterative process Methods 0.000 claims 1
- 238000000605 extraction Methods 0.000 abstract description 4
- 238000012360 testing method Methods 0.000 description 5
- 238000007635 classification algorithm Methods 0.000 description 3
- 230000000694 effects Effects 0.000 description 3
- 238000010801 machine learning Methods 0.000 description 3
- 230000008901 benefit Effects 0.000 description 2
- 230000015572 biosynthetic process Effects 0.000 description 2
- 239000002131 composite material Substances 0.000 description 2
- 238000010276 construction Methods 0.000 description 2
- 238000003064 k means clustering Methods 0.000 description 2
- 238000003786 synthesis reaction Methods 0.000 description 2
- 101000701876 Homo sapiens Serpin A9 Proteins 0.000 description 1
- 238000007476 Maximum Likelihood Methods 0.000 description 1
- 238000004458 analytical method Methods 0.000 description 1
- 230000001419 dependent effect Effects 0.000 description 1
- 102000050111 human SERPINA9 Human genes 0.000 description 1
- 238000011835 investigation Methods 0.000 description 1
- 238000013507 mapping Methods 0.000 description 1
- 239000000463 material Substances 0.000 description 1
- 238000005259 measurement Methods 0.000 description 1
- 238000013508 migration Methods 0.000 description 1
- 230000005012 migration Effects 0.000 description 1
- 238000012545 processing Methods 0.000 description 1
- 238000011160 research Methods 0.000 description 1
- 238000010187 selection method Methods 0.000 description 1
- 238000006467 substitution reaction Methods 0.000 description 1
- 230000009466 transformation Effects 0.000 description 1
- 238000000844 transformation Methods 0.000 description 1
- 230000000007 visual effect Effects 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F18/00—Pattern recognition
- G06F18/20—Analysing
- G06F18/24—Classification techniques
- G06F18/243—Classification techniques relating to the number of classes
- G06F18/24323—Tree-organised classifiers
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F18/00—Pattern recognition
- G06F18/20—Analysing
- G06F18/21—Design or setup of recognition systems or techniques; Extraction of features in feature space; Blind source separation
- G06F18/214—Generating training patterns; Bootstrap methods, e.g. bagging or boosting
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F18/00—Pattern recognition
- G06F18/20—Analysing
- G06F18/23—Clustering techniques
- G06F18/232—Non-hierarchical techniques
- G06F18/2321—Non-hierarchical techniques using statistics or function optimisation, e.g. modelling of probability density functions
- G06F18/23213—Non-hierarchical techniques using statistics or function optimisation, e.g. modelling of probability density functions with fixed number of clusters, e.g. K-means clustering
Landscapes
- Engineering & Computer Science (AREA)
- Data Mining & Analysis (AREA)
- Theoretical Computer Science (AREA)
- Physics & Mathematics (AREA)
- Bioinformatics & Cheminformatics (AREA)
- Artificial Intelligence (AREA)
- Life Sciences & Earth Sciences (AREA)
- Bioinformatics & Computational Biology (AREA)
- Computer Vision & Pattern Recognition (AREA)
- Evolutionary Biology (AREA)
- Evolutionary Computation (AREA)
- General Engineering & Computer Science (AREA)
- General Physics & Mathematics (AREA)
- Probability & Statistics with Applications (AREA)
- Image Analysis (AREA)
Abstract
本发明涉及一种基于GlobeLand30的历史地表覆盖快速重建方法，包括以下步骤：1）GlobeLand30的有效斑块提取：利用形态学运算去除错误斑块，获得可以正确表达地表覆盖信息空间分布的斑块。2）伪样本推选：计算影像光谱特征，以优化后的斑块作为聚类单元，采用无监督算法推选出伪样本集。3）全局样本优化：对不同地物类型分别构建高斯混合模型，通过求解模型的过程优化全局样本。4）地表覆盖重建：根据第三步结果获得训练样本集，判断是否需要进行变化检测，将样本集输入随机森林，重建目标历史时相的地表覆盖分类图。本算法考虑了地表覆盖产品中信息的有效性和不确定性，通过有效降低不确定性快速重建历史地表覆盖。
Description
技术领域
本发明涉及一种基于GlobeLand30的历史地表覆盖快速重建方法，属于遥感智能信息提取技术领域。
背景技术
地表覆盖是描述地球表面构成的关键指标，对全球陆表水分、能量与物质循环有重要影响。遥感影像监督分类解译地表覆盖已成为主流的技术方法，通过对遥感影像进行目视解译，手动选取训练样本，基于监督分类算法生成多期地表覆盖分类结果。这类方法得益于监督分类方法有效性，在遥感地表覆盖分类中得到了广泛的应用。监督分类结果对训练样本非常依赖，训练样本获取环节对上述方法尤其重要。样本选取是一个依赖专家知识且需要花费一定时间与人力的过程，这个过程需要高分影像作为参考或者实地调查获得地表真实数据。
因此，对于历史影像地表覆盖重建问题，传统基于机器学习的遥感影像分类方法存在以下局限性：1)理想的监督分类结果依赖充分的高质量训练样本，需要大量的时间与人力的投入；2)对于历史影像，用于标记参考的高分影像非常缺乏，难以覆盖整个研究区域，使得训练样本集不能很好地代表整体研究区的地物分布。针对以上问题，提出了一种快速、无监督的历史地表覆盖快速重建方法，充分利用已有地表覆盖产品，实现历史影像训练样本的自动获取。构建一种无监督样本迁移方法嵌入地表覆盖分类框架中，快速准确地重建历史地表覆盖。
发明内容
本发明要解决的技术问题是：克服历史地表覆盖分类任务中训练样本获取难度大、成本高的不足，提出一种无监督样本迁移的方法，该方法利用存档地表覆盖产品数据，通过产品斑块的几何约束、影像光谱特征的属性约束以及地物的分布约束充分降低地表覆盖产品的时空不确定性，自动快速生产高质量历史影像训练样本，结合变化检测与分类算法，快速准确地完成历史地表覆盖的重建。
为了解决以上技术问题，本发明提供的一种基于GlobeLand30的历史地表覆盖快速重建方法，包括以下步骤：
第一步、准备地表覆盖产品GlobeLand30、遥感影像Landsat数据以及DEM数据，选择感兴趣区域并生成感兴趣区域的矢量文件，对GlobeLand30、DEM与遥感影像Landsat数据进行裁剪；
第二步、基于形态学开运算构建对地表覆盖产品土地斑块的优化从而获得斑块单元：将地表覆盖产品GlobeLand30按照类别划分为若干二值图层，同时对每一层二值图像进行窗口大小递增的形态学开运算，针对不同窗口大小下的计算结果，当前背景像元占比大于用户设定阈值时，对二值图层进行合并，获得优化后的斑块单元；
第三步、计算影像归一化差值光谱向量NDSV：通过遍历所有波段反射率的两两组合，对每一对波段组合计算归一化差值指数，获得对应的归一化差值光谱向量NDSV，所有的归一化差值光谱向量NDSV构成了影像光谱特征集；
第四步、将每一个斑块单元作为聚类计算单元，采用K-means算法进行无监督聚类，采用K-means++算法逐个初始化每个斑块单元的聚类中心：对于任意斑块单元，从中随机选择一个像元的归一化差值光谱向量NDSV作为第一个簇的聚类中心，当目前有k个聚类中心时，计算当前斑块单元内其他像元的归一化差值光谱向量NDSV到第k个聚类中心μk的欧式距离，并根据该欧式距离设置对应像元的抽样权重，根据抽样权重随机抽取第k+1个聚类中心，重复该过程直到簇的数量满足用户预设要求；
第五步、采用方差比准则VRC进行每个斑块单元聚类簇数的寻优：定义一个聚类簇数集合，对于任意斑块单元，将所有像元的归一化差值光谱向量NDSV输入K-means算法进行聚类，根据聚类结果计算方差比准则VRC结果，通过遍历簇数集合内的所有元素，获得最大方差比准则VRC对应的簇数Kj作为当前斑块单元内聚类的最优簇数；
第六步、根据第四、第五步方法，对所有斑块单元执行K-means方法，根据最优簇数聚类结果获得全局影像伪样本集：将第三步获得的光谱影像特征集作为K-means输入特征，遍历第二步获得的所有斑块单元，对各个斑块单元内的所有像元使用第四步方法获得聚类中心，采用第五步方法获得最优簇数，在此基础上获得聚类结果，令斑块单元内像元数量占比最大的簇继承斑块单元的类别信息，获得伪样本集；
第七步、根据地表覆盖产品确定地表覆盖类别体系，对每一个地类类别选取对应的关键特征，构建多个二类高斯混合模型：计算全局影像的改进归一化水体指数MNDWI作为水体类型的关键特征，计算非水体区域的归一化植被指数NDVI分别作为林地、草地、人工地表类型的关键特征，计算植被区域的坡度作为耕地类型的关键特征，将每一个关键特征的特征值分布视为两个高斯分布的混合，构建如下高斯混合模型：
其中
为概率密度函数，y表示关键特征的特征值，α1、α2分别为第一和第二个高斯分布的混合系数、μ1、μ2分别为两个分布的均值向量，∑1、∑2分布为两个分布的协方差向量，p(y|μ1,∑1)与p(y|μ2,∑2)分别为第一与第二个高斯分布密度函数。以上{(αj,μj,∑j)|j＝1,2}统称为混合成分参数；
第八步、采用迭代期望最大化EM算法估算每个高斯混合模型中的混合成分参数，完成对高斯混合模型的分解，从而完成对伪样本集的优化，获得训练样本集；
第九步、利用变化检测与分类方法完成目标历史时相的地表覆盖重建：判断目标时相与地表覆盖产品是否属于同一年份，如果属于同一年份，则将第八步中优化后的训练样本输入随机森林RF分类器完成训练，获得目标时相地表覆盖分类结果；如果属于不同年份，则采用双时相影像变化检测方法获得不变区域，传递不变区域样本标签，结合目标时相影像形成训练样本，训练样本输入随机森林RF分类器完成训练，获得目标历史时相地表覆盖分类结果。
本发明采用的历史地表覆盖重建方法提取属于算法创新，算法核心思想是利用存档地表覆盖产品中的信息，同时考虑产品信息的有效性和不确定，从地表覆盖产品中土地斑块信息的角度来初步降低产品的地表覆盖信息不确定性。本发明引入改进的K-means算法，在斑块单元的几何约束和光谱特征的属性约束下，快速划分获得伪样本集。本发明引入高斯混合模型，利用不同地类在关键特征空间中的分布约束，通过混合模型分解优化伪样本集，得到可靠的训练样本集。本发明基于两种基础聚类算法，从地表覆盖产品斑块的几何约束、影像光谱特征的属性约束以及地类特征空间的分布约束三个角度有效降低产品信息的不确定性，自动获得可靠的历史时相训练样本集。本发明提出的算法是一种快速、自动的历史地表覆盖重建方法，为遥感影像智能解译和地表覆盖动态制图问题提供了解决思路。所有步骤均由Matlab编程实现，除了设置算法必要参数外无需人工参与。
算法的主要创新在以下三方面：
1、所述第二步中，构建了形态学开运算为基础构建的斑块单元优化方法，通过逐个类别的二值图像开运算达到去除原始产品中细小斑块与错误斑块联合的目的，得到更加符合算法需求的斑块单元。
2、所述第六步中，提出了以改进K-means算法为基础的伪样本推选方法，算法为第二步优化后的土地斑块为聚类单元，以第三步计算的NDSV特征为聚类特征，实现从局部单元显著降低产品不确定性，同时伪样本推选算法计算单元相互独立，采用并行计算保证算法运行效率。
3、所述第七步中，构建了基于高斯混合模型的全局样本优化方法，通过选择与地表覆盖类型相关的关键特征，从特征空间数据分布的角度对伪样本进行全局优化，获得可靠的训练样本集。
本发明方法在多个历史年份的大尺度地表覆盖重建中样本提取效果好，分类精度高，相较于原始地表覆盖产品更加准确，对中高分辨率历史地表覆盖重建与新影像的地表覆盖更新起到很好的借鉴效果。
附图说明
结合附图对本发明作进一步的说明。
图1为基于GlobeLand30的历史地表覆盖快速重建方法流程图。
图2为伪样本推选方法流程图。
图3为基于高斯混合模型的全局伪样本优化流程图。
图4(a)为2000年中值合成影像。
图4(b)为2000年地表覆盖分类结果。
图5(a)为1995年中值合成影像。
图5(b)为1995年地表覆盖分类结果。
具体实施方式
根据附图详细阐述本发明的技术路线和具体操作步骤。
本发明实例研究区域为太湖流域地区，地表覆盖产品选择GlobeLand30系列，产品数据获取时间为2000年，为了方便描述，后续将产品简称为GlobeLand30V20000，影像获取时间分别为1995年与2000年。
本实施案例基于GlobeLand30的历史地表覆盖快速重建方法(图1)，包括以下步骤：
第一步：下载GlobeLand30 V2000数据产品(空间分辨率为30m)，根据试验区矢量文件裁剪数据产品。通过谷歌地球引擎(Google Earth Engine)获取Landsat中值合成的反射率数据(空间分辨率为30m)，该数据产品在GEE的目录为"LANDSAT/LT05/C02/T1_L2"，该目录下的Landsat数据附带CFMask算法生成的云、阴影、雪等地物的掩膜文件，基于掩膜文件在GEE中获取1995、2000年的全部晴空观测数据，通过中值合成的方式输出影像，根据试验区矢量文件裁剪影像，下载数字高程模型(Digital Elevation Model，DEM)数据，根据试验区矢量文件裁剪。
第二步：利用Matlab编程完成对土地斑块的优化方法。读取栅格格式的GlobeLand30 V2000，按照类别栅格值将GlobeLand30 V2000分解成多个二值图层，设置结构元形状为正方形，领域大小范围为[3,21]，步长为2。利用Matlab图像处理工具箱中的imopen函数完成形态学开运算，在完成每个窗口值形态学开运算后，将二值图像合并，计算背景像元比例，记为ξ。为了保证在优化斑块的同时尽可能保留产品信息，将ξ设置阈值为5％，当背景像元比例高于5％(用户设定阈值)时，停止优化算法，将当前窗口值下的合并结果作为优化后的斑块单元。
第三步：利用Matlab编程完成归一化差值光谱向量NDSV的提取，读取2000年年度中值合成的Landsat 5影像，遍历所有波段的两两组合，对每一对组合的两个波段计算归一化差值指数，获得对应的归一化差值光谱向量NDSV，所有的归一化差值光谱向量NDSV形成对各个地类更具判别性的影像光谱特征集，归一化差值光谱向量NDSV具体的计算公式入下所示：
其中yd代表影像中第d个像元的光谱特征向量，
为第d个像元第bi个波段的反射率，
为第d个像元第bj个波段的反射率，i、j的取值范围为[1，B]。
第四步：将每一个斑块单元作为聚类计算单元，采用K-means算法进行无监督聚类，采用K-means++算法逐个初始化每个斑块单元的聚类中心。对于任意斑块单元pj，根据第三步NDSV计算方法，获得其对应的光谱特征集
从
中随机选择一个归一化差值光谱向量作为第一个簇的聚类中心，记为μ1；然后计算
中每一个归一化差值光谱向量到聚类中心的欧式距离，在
中根据抽样权重随机选择的下一个聚类中心μ2，通过归一化差值光谱向量到μ1的欧式距离来地设置对应的抽样权重，对于第n个特征向量yn而言，抽样权重
的计算方式如下：
当获得第k个聚类中心时，从
中继续按照权重随机选取下一个聚类中心，对于属于簇
的特征向量yn而言，其权重的计算方式为：
重复k→k+1的过程直到簇的数量满足预设要求。
第五步：利用Matlab编程实现方差比准则VRC算法，自动寻优每一个斑块单元的最佳聚类簇数。以任意斑块单元pj为例，利用K-means将光谱特征集
聚类为Kj个簇，在聚类结果的基础上，获得当前聚类结果下的总体簇间方差SSB(overall between-clustervariance)，当前聚类结果下的总体簇间方差SSB的计算方式如下：
其中μk是簇
的均值向量，Kj为簇的数量，μ为数据集的整体均值向量。
是簇
内的像元数量。总体簇间方差SSB越大表示簇与簇之间的相似度越低。
计算当前聚类结果下的总体簇内方差SSW(overall within-cluster variance)，通过下面的公式计算：
其中yn为特征向量，μk是簇
的均值向量，Kj为簇的数量，而总体簇内方差越小则表示簇内数据之间相似度越高。
方差比准则通过计算总体簇间方差和总体簇内反差的比值来衡量聚类效果的好坏，即保证簇内数据相似性和簇间数据的差异性，具体衡量指标计算如下：
其中N为特征向量数量，Kj为簇的数量，设置最优簇数的取值范围为1到10。
第六步：利用Matlab编程实现伪样本推选算法(图2)，提出的伪样本推选算法需要三类输入数据，首先将第二步中优化后的斑块单元进行顺序编码，得到斑块单元的编码栅格数据，利用Matlab读取；其次读取归一化差值光谱向量NDSV的特征影像，最后读取地表覆盖产品栅格数据。按照编码信息顺序逐个获取单个斑块单元，作为伪样本推选算法的几何约束，即改进K-means方法的聚类单元，同时获取该斑块单元对应的NDSV特征影像子集和地表覆盖产品类别信息，将NDSV特征作为属性约束，即改进K-means方法的聚类特征。对每个独立的斑块单元，执行改进的K-means算法，在完成初始化跟簇数寻优后，利用Matlab统计与机器学习工具箱中的kmeans函数，选择输入数据为斑块单元内的NDSV特征，聚类簇数为寻优后的结果，最大迭代次数设置为1000，完成K-means聚类。算法通过定义的初始的Kj个聚类中心，通过计算
中的每个归一化差值光谱向量yn到每个聚类中心欧式距离，将yn划分到最相似聚类中心的簇中。通过不断调整聚类中心，最小化平方误差：
其中，μk是簇
的均值向量。
当K-means算法收敛时，获取当前簇划分结果，将占比最高的簇保留，同时将地表覆盖产品类别信息赋予该簇。通过并行计算对试验区所有斑块单元完成伪样本推选算法，获得整个试验区的伪样本推选结果。
第七步：为各个地类选取关键特征，选择改进的归一化水体指数(ModifiedNormalized Difference Water Index,MNDWI)作为水体与其他地类之间的判别信息。MNDWI通过绿光波段反射率与中红外波段的反射率的归一化差值获得：
其中ρMIR与ρGreen分别代表中红外波段与绿光波段的反射率。
选择归一化植被指数(Normalized Difference Vegetation Index,NDVI)作为林地、草地以及人工地表的关键指标，NDVI通过近红外波段与红光波段的归一化差值指数获得：
其中ρNIR与ρRed分别代表近红外波段与红光波段的反射率。
选择坡度作为耕地的关键指标，坡度基于DEM数据计算获得，通过ArcGIS软件执行获得，首先选择空间选取工具，打开“表面”子栏，点击“坡度”工具，输入DEM计算获得。
利用Matlab编程实现多层二类高斯混合模型的构建(图3)，通过将整体影像的MNDWI视为两类高斯分布的混合，在求解二类高斯混合模型的基础上，通过比较簇均值将簇标记为水体和非水体两类，全局影像的非水体部分继续分别通过NDVI和坡度特征构建高斯混合模型。其中NDVI构建的高斯混合模型同时用于林草地和人工地表伪样本的优化，耕地的伪样本则通过坡度构建的高斯混合模型进行分解优化。二类高斯混合模型的构建采用Matlab统计与机器学习工具箱中的fitgmdist函数编程完成。
高斯混合模型公式如下：
其中
为概率密度函数，y表示关键特征的特征值，α1、α2分别为第一和第二个高斯分布的混合系数、μ1、μ2分别为两个分布的均值向量，∑1、∑2分布为两个分布的协方差向量，p(y|μ1,∑1)与p(y|μ2,∑2)分别为第一与第二个高斯分布密度函数。以上{(αj,μj,Σj)|j＝1,2}统称为混合成分参数。
第八步：利用Matlab编程完成基于高斯混合模型分解的全局伪样本优化，获得可靠的训练样本集。采用EM算法完成对第七步中构建的多个二类高斯混合模型的分解，首先采用马氏距离作为权重，使用权重采样的方式完成EM算法的初始化，然后实行EM算法的E步，利用当前划分结果的混合成分参数，估算数据集中的元素
属于第一和第二个高斯混合模型的后验概率，具体计算方式为：
其中
为隶属j的后验概率，{(αj,μj,∑j)|j＝1,2}为混合成分参数，在此基础上执行M步，通过最大似然估计求解调整混合分布参数，需要最大化的对数似然函数为：
其中
为当前数据集，Ni为数据集中元素的数量。
设置EM算法的停止条件为似然函数的变化不再显著，本实施例中具体参数为ΔLL≤10-6。当完成对多个高斯混合模型的分解后，将伪样本集中的错误标签的样本去除，获得最终的训练样本结果。
第九步：利用Matlab编程实现变化检测与分类算法，获得地表覆盖重建结果。对于2000年份的地表覆盖重建，直接将第八步获得的训练样本与2000年的影像数据输入随机森林(Random Forest)分类器，RF分类器采用开源代码实现(Classification andregression based on a forest of trees using random inputs,based on Breiman(2001)<DOI:10.1023/A:1010933404324>.)，其中RF的基分类器数目设置为500，RF的特征子集的维数设置为输入特征维数的平方根。利用RF获得2000年地表覆盖重建结果(图4)。
对于1995年份的地表覆盖重建任务，考虑到地表覆盖变化的影响，采用变化矢量分析法结合大津法阈值分割，检测1995-2000年间的不变区域。变化检测的算法通过Matlab编程实现，首先计算变化强度信息：
其中ρd为第个像元的变化强度信息，
与
分别为该像元的第一个时相和第二时相在n波段的反射率值，N为影像波段总数。然后采用大津法实现对二值化强度图像，获得变化与不变区域。大津法通过遍历图像中的所有灰度级，计算不同灰度级阈值下的类间方差
的集合，将集合内的极大值对应的阈值T作为最佳阈值。对于灰度级为{0,1,2,3,…,l,…L-1}的图像，最佳阈值T*为:
其中T为阈值，L为图像灰度级数量，
为类间方差，根据最佳阈值将强度图像二值化，获得不变区域，将不变区域的训练样本标签从2000年影像传递到1995年的影像中，形成1995年的训练样本，输入RF分类器完成1995的地表覆盖重建(图5)。
除上述实施例外，本发明还可以有其他实施方式。凡采用等同替换或等效变换形成的技术方案，均落在本发明要求的保护范围。
Claims (8)
1.一种基于GlobeLand30的历史地表覆盖快速重建方法，包括以下步骤：
第一步、准备地表覆盖产品GlobeLand30、遥感影像Landsat数据以及DEM数据，选择感兴趣区域并生成感兴趣区域的矢量文件，对GlobeLand30、DEM与遥感影像Landsat数据进行裁剪；
第二步、基于形态学开运算构建对地表覆盖产品土地斑块的优化从而获得斑块单元：将地表覆盖产品GlobeLand30按照类别划分为若干二值图层，同时对每一层二值图像进行窗口大小递增的形态学开运算，针对不同窗口大小下的计算结果，当前背景像元占比大于用户设定阈值时，对二值图层进行合并，获得优化后的斑块单元；
第三步、计算影像归一化差值光谱向量NDSV：通过遍历所有波段反射率的两两组合，对每一对波段组合计算归一化差值指数，获得对应的归一化差值光谱向量NDSV，所有的归一化差值光谱向量NDSV构成了影像光谱特征集；
第四步、将每一个斑块单元作为聚类计算单元，采用K-means算法进行无监督聚类，采用K-means++算法逐个初始化每个斑块单元的聚类中心：对于任意斑块单元，从中随机选择一个像元的归一化差值光谱向量NDSV作为第一个簇的聚类中心，当目前有k个聚类中心时，计算当前斑块单元内其他像元的归一化差值光谱向量NDSV到第k个聚类中心μk的欧式距离，并根据该欧式距离设置对应像元的抽样权重，根据抽样权重随机抽取第k+1个聚类中心，重复该过程直到簇的数量满足用户预设要求；
第五步、采用方差比准则VRC进行每个斑块单元聚类簇数的寻优：定义一个聚类簇数集合，对于任意斑块单元，将所有像元的归一化差值光谱向量NDSV输入K-means算法进行聚类，根据聚类结果计算方差比准则VRC结果，通过遍历簇数集合内的所有元素，获得最大方差比准则VRC对应的簇数Kj作为当前斑块单元内聚类的最优簇数；
第六步、根据第四、第五步方法，对所有斑块单元执行K-means方法，根据最优簇数聚类结果获得全局影像伪样本集：将第三步获得的光谱影像特征集作为K-means输入特征，遍历第二步获得的所有斑块单元，对各个斑块单元内的所有像元使用第四步方法获得聚类中心，采用第五步方法获得最优簇数，在此基础上获得聚类结果，令斑块单元内像元数量占比最大的簇继承斑块单元的类别信息，获得伪样本集；
第七步、根据地表覆盖产品确定地表覆盖类别体系，对每一个地类类别选取对应的关键特征，构建多个二类高斯混合模型：计算全局影像的改进归一化水体指数MNDWI作为水体类型的关键特征，计算非水体区域的归一化植被指数NDVI分别作为林地、草地、人工地表类型的关键特征，计算植被区域的坡度作为耕地类型的关键特征，将每一个关键特征的特征值分布视为两个高斯分布的混合，构建如下高斯混合模型：
其中
为概率密度函数，y表示关键特征的特征值，α1、α2分别为第一和第二个高斯分布的混合系数、μ1、μ2分别为两个分布的均值向量，∑1、∑2分布为两个分布的协方差向量，p(y|μ1,∑1)与p(y|μ2,∑2)分别为第一与第二个高斯分布密度函数。以上{(αj,μj,∑j)|j＝1,2}统称为混合成分参数；
第八步、采用迭代期望最大化EM算法估算每个高斯混合模型中的混合成分参数，完成对高斯混合模型的分解，从而完成对伪样本集的优化，获得训练样本集；
第九步、利用变化检测与分类方法完成目标历史时相的地表覆盖重建：判断目标时相与地表覆盖产品是否属于同一年份，如果属于同一年份，则将第八步中优化后的训练样本输入随机森林RF分类器完成训练，获得目标时相地表覆盖分类结果；如果属于不同年份，则采用双时相影像变化检测方法获得不变区域，传递不变区域样本标签，结合目标时相影像形成训练样本，训练样本输入随机森林RF分类器完成训练，获得目标历史时相地表覆盖分类结果。
2.根据权利要求1所述的基于GlobeLand30的历史地表覆盖快速重建方法，其特征在于：第二步所述的土地斑块优化方法，其中形态学结构元形状为正方形，领域大小范围为[3,21],步长为2，背景像素比例ξ设置阈值为5％。
3.根据权利要求1所述的基于GlobeLand30的历史地表覆盖快速重建方法，其特征在于：第三步所述的归一化差值光谱向量NDSV特征计算公式如下：
其中yd代表影像中第d个像元的光谱特征向量，
为第d个像元第bi个波段的反射率，
为第d个像元第bj个波段的反射率，i、j的取值范围为[1，B]。
4.根据权利要求1所述的基于GlobeLand30的历史地表覆盖快速重建方法，其特征在于：第三步所述的方差比准则设置簇数的取值范围为1到10的自然数。
5.根据权利要求1所述的基于GlobeLand30的历史地表覆盖快速重建方法，其特征在于：第六步所述改进的K-means算法聚类单元设置为第二步的斑块单元，聚类特征设置为第三步计算的归一化差值光谱向量NDSV特征。
6.根据权利要求1所述的基于GlobeLand30的历史地表覆盖快速重建方法，其特征在于：第八步中，通过对迭代期望最大化EM算法中两次似然函数的变化量设置阈值，当似然函数不再显著变化时，停止迭代过程，完成对高斯混合模型的分解，将伪样本集中从错误标签去除，获得训练样本集；其中，似然函数变化量阈值设置为ΔLL≤10-6。
7.根据权利要求1所述的基于GlobeLand30的历史地表覆盖快速重建方法，其特征在于：第九步中随机森林RF分类器参数设置为决策树的数量Ntree设置为500,基分类器的特征子集维度Mtry设置为输入分类器的特征集维数的平方根。
8.根据权利要求1所述的基于GlobeLand30的历史地表覆盖快速重建方法，其特征在于：第九步中变化检测算法设置为在变化矢量分析强度图像结果基础上采用大津法进行阈值分割。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202111580258.XA CN114254707A (zh) | 2021-12-22 | 2021-12-22 | 一种基于GlobeLand30的历史地表覆盖快速重建方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202111580258.XA CN114254707A (zh) | 2021-12-22 | 2021-12-22 | 一种基于GlobeLand30的历史地表覆盖快速重建方法 |
Publications (1)
| Publication Number | Publication Date |
|---|---|
| CN114254707A true CN114254707A (zh) | 2022-03-29 |
Family
ID=80796691
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202111580258.XA Pending CN114254707A (zh) | 2021-12-22 | 2021-12-22 | 一种基于GlobeLand30的历史地表覆盖快速重建方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN114254707A (zh) |
Cited By (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN115205675A (zh) * | 2022-06-29 | 2022-10-18 | 中国科学院地理科学与资源研究所 | 一种基于多源遥感数据的森林变化驱动力分类方法 |
| CN115346120A (zh) * | 2022-08-16 | 2022-11-15 | 中国科学院空天信息创新研究院 | 一种草地地上生物量及其固碳量遥感估算方法 |
| CN116912202A (zh) * | 2023-07-13 | 2023-10-20 | 中国中医科学院眼科医院 | 一种医用高值耗材管理方法和系统 |
| CN117577563A (zh) * | 2024-01-16 | 2024-02-20 | 东屹半导体科技(江苏)有限公司 | 一种半导体划片机的优化控制方法及系统 |
-
2021
- 2021-12-22 CN CN202111580258.XA patent/CN114254707A/zh active Pending
Cited By (7)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN115205675A (zh) * | 2022-06-29 | 2022-10-18 | 中国科学院地理科学与资源研究所 | 一种基于多源遥感数据的森林变化驱动力分类方法 |
| CN115205675B (zh) * | 2022-06-29 | 2023-05-16 | 中国科学院地理科学与资源研究所 | 一种基于多源遥感数据的森林变化驱动力分类方法 |
| CN115346120A (zh) * | 2022-08-16 | 2022-11-15 | 中国科学院空天信息创新研究院 | 一种草地地上生物量及其固碳量遥感估算方法 |
| CN116912202A (zh) * | 2023-07-13 | 2023-10-20 | 中国中医科学院眼科医院 | 一种医用高值耗材管理方法和系统 |
| CN116912202B (zh) * | 2023-07-13 | 2024-01-30 | 中国中医科学院眼科医院 | 一种医用高值耗材管理方法和系统 |
| CN117577563A (zh) * | 2024-01-16 | 2024-02-20 | 东屹半导体科技(江苏)有限公司 | 一种半导体划片机的优化控制方法及系统 |
| CN117577563B (zh) * | 2024-01-16 | 2024-04-12 | 东屹半导体科技(江苏)有限公司 | 一种半导体划片机的优化控制方法及系统 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN114254707A (zh) | 一种基于GlobeLand30的历史地表覆盖快速重建方法 | |
| Petropoulos et al. | Support vector machines and object-based classification for obtaining land-use/cover cartography from Hyperion hyperspectral imagery | |
| Du et al. | Mapping wetland plant communities using unmanned aerial vehicle hyperspectral imagery by comparing object/pixel-based classifications combining multiple machine-learning algorithms | |
| CN103679675B (zh) | 一种面向水质定量遥感应用的遥感影像融合方法 | |
| CN102982338B (zh) | 基于谱聚类的极化sar图像分类方法 | |
| CN112101271A (zh) | 一种高光谱遥感影像分类方法及装置 | |
| CN105608473B (zh) | 一种基于高分辨率卫星影像的高精度土地覆盖分类方法 | |
| CN110309780A (zh) | 基于bfd-iga-svm模型的高分辨率影像房屋信息快速监督识别 | |
| Luo et al. | Comparison of machine learning algorithms for mapping mango plantations based on Gaofen-1 imagery | |
| CN113033714B (zh) | 多模态多粒度遥感影像面向对象全自动机器学习方法及系统 | |
| Zhang et al. | Mapping freshwater marsh species in the wetlands of Lake Okeechobee using very high-resolution aerial photography and lidar data | |
| CN109359525B (zh) | 基于稀疏低秩的判别谱聚类的极化sar图像分类方法 | |
| Zafari et al. | A multiscale random forest kernel for land cover classification | |
| Bagan et al. | Improved subspace classification method for multispectral remote sensing image classification | |
| Sharma et al. | Land cover classification: a comparative analysis of clustering techniques using Sentinel-2 data | |
| Kar et al. | Classification of multispectral satellite images | |
| Vignesh et al. | A novel multiple unsupervised algorithm for land use/land cover classification | |
| Jenicka | Land Cover Classification of Remotely Sensed Images | |
| Li et al. | Measuring detailed urban vegetation with multisource high-resolution remote sensing imagery for environmental design and planning | |
| Acar et al. | Performance Assessment of Landsat 8 and Sentinel-2 Satellite Images for the Production of Time Series Land Use/Land Cover (Lulc) Maps | |
| Bai et al. | Kernel low-rank entropic component analysis for hyperspectral image classification | |
| Selvaraj et al. | Assessment of object-based classification for mapping land use and land cover using google earth | |
| Kavzoglu et al. | Agricultural crop type mapping using object-based image analysis with advanced ensemble learning algorithms | |
| Sekertekin | Mapping of surface water resources: a comparative analysis of eight image classification methods | |
| Dumeur et al. | Paving the way toward foundation models for irregular and unaligned Satellite Image Time Series |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination |