CN108089850B - 一种基于影像协同分割与生态地理分区规则库的地表覆盖产品增量更新方法 - Google Patents

一种基于影像协同分割与生态地理分区规则库的地表覆盖产品增量更新方法 Download PDF

Info

Publication number
CN108089850B
CN108089850B CN201810001776.3A CN201810001776A CN108089850B CN 108089850 B CN108089850 B CN 108089850B CN 201810001776 A CN201810001776 A CN 201810001776A CN 108089850 B CN108089850 B CN 108089850B
Authority
CN
China
Prior art keywords
change
image
ecological
pseudo
rule
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.)
Expired - Fee Related
Application number
CN201810001776.3A
Other languages
English (en)
Other versions
CN108089850A (zh
Inventor
朱凌
石若明
谢振雷
陈旭
孙杨
剌怡璇
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Beijing University of Civil Engineering and Architecture
Original Assignee
Beijing University of Civil Engineering and Architecture
Priority date (The priority date 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 date listed.)
Filing date
Publication date
Application filed by Beijing University of Civil Engineering and Architecture filed Critical Beijing University of Civil Engineering and Architecture
Priority to CN201810001776.3A priority Critical patent/CN108089850B/zh
Publication of CN108089850A publication Critical patent/CN108089850A/zh
Application granted granted Critical
Publication of CN108089850B publication Critical patent/CN108089850B/zh
Expired - Fee Related legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F8/00—Arrangements for software engineering
    • G06F8/20—Software design
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F18/00—Pattern recognition
    • G06F18/20—Analysing
    • G06F18/24—Classification techniques
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00—Image analysis
    • G06T7/0002—Inspection of images, e.g. flaw detection
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00—Image analysis
    • G06T7/10—Segmentation; Edge detection
    • G06T7/136—Segmentation; Edge detection involving thresholding
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06V—IMAGE OR VIDEO RECOGNITION OR UNDERSTANDING
    • G06V10/00—Arrangements for image or video recognition or understanding
    • G06V10/40—Extraction of image or video features
    • G06V10/46—Descriptors for shape, contour or point-related descriptors, e.g. scale invariant feature transform [SIFT] or bags of words [BoW]; Salient regional features
    • G06V10/462—Salient features, e.g. scale invariant feature transforms [SIFT]
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/10—Image acquisition modality
    • G06T2207/10032—Satellite or aerial image; Remote sensing
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/30—Subject of image; Context of image processing
    • G06T2207/30181—Earth observation

Landscapes

  • Engineering & Computer Science (AREA)
  • Theoretical Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • General Engineering & Computer Science (AREA)
  • Data Mining & Analysis (AREA)
  • Software Systems (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Evolutionary Computation (AREA)
  • Evolutionary Biology (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Artificial Intelligence (AREA)
  • Multimedia (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Quality & Reliability (AREA)
  • Image Processing (AREA)

Abstract

一种基于影像协同分割与生态地理分区规则库的地表覆盖产品增量更新方法，包括协同分割提取两期影像的变化图斑；利用生态地理分区规则库识别剔除伪变化；地表覆盖产品的更新；具体步骤如下：一、多时相遥感影像协同分割技术；二、基于生态地理分区规则库的变化图斑的去伪技术；三、地表覆盖产品更新。本发明可以解决利用多时相卫星影像，采用协同分割技术提取地表覆盖增量，即地表覆盖的变化，之后对变化图斑进行分类，再结合全球生态地理分区规则库对分类增量图斑进行去伪变化，提高增量检测结果的可靠度，最后用增量图斑更新旧时相地表覆盖产品，获得新的地表覆盖产品的技术问题。

Description

一种基于影像协同分割与生态地理分区规则库的地表覆盖产 品增量更新方法
技术领域
本发明涉及基于遥感影像的地表变化检测、地表覆盖产品更新领域。
背景技术
目前全球地表覆盖产品及更新情况如下：
1、美国地质调查局(USGS)的方法生产的，Global Land Cover CharacteristicsDatabase (GLCC)数据库中的地表覆盖产品，IGBP－DIS(International Geosphere-Biosphere Programme Data and Information System，简称IGBP-DIS)产品主要是通过非监督分类1公里分辨率的AVHRR的NDVI产品而获得。影像是1992年4月至1993年3月的AVHRR遥感数据。IGBP-DIS数据分类系统采取IGBP的分类体系【1】。
2、马里兰大学(UMD)地理系利用与IGBP-DIS相同的影像数据数据，包括AVHRR卫星影像5个波段影像和NDVI数据，采用决策树监督分类方法制作了另一个地表覆盖产品，产品包含12个地表覆盖类。产品的分辨率为8公里和1公里。【2】
3、欧盟土地覆被数据GLC2000。基于SPOT4卫星VEGETATION-1传感器1999年11月至2000年12月影像，分类方法主要采用非监督分类的方法。全球产品是通过混合各区域产品，以LCCS图例概括各区域图例，最终得到GLC2000产品【3】。
4、GLCNMO
全球制图国际指导委员会的全球制图计划生产了全球1公里分辨率的地表覆盖数据集GLCNMO(Global Land Cover by National Mapping Organizations)，产品采用LCCS分类体系定义了20类地表覆盖类型，其中14类是通过监督分类方法获得，剩余的6类包urban,tree open,mangrove,wetland,snow/ice,water是分别独立分类的。基础数据源是2003年的7波段1公里分辨率的MODIS影像。【4】
以上四种产品都只推出了一期产品，之后未见更新或新一期产品推出。
5、Globcover 2005和2009产品。300米分辨的数据产品Globcover是由欧洲空间局(ESA)发起。Globcover发布了两期产品：Globcover2005和2009。
数据生产是采用MERIS FRS level 1B影像，根据UN LCCS分类系统对影像进行分类.【5】
值得一提的是，GlobCover 2009产品分类时主要参考GlobCover 2005(V2.2)产品。确定类型既考虑时空类的类型，还考虑GlobCover2005的类型。通过行业专家制定了一些规则以帮助确定类型。
GlobCover 2005和2009采用了相同的分类算法、相同的分类体系。但是，Bicheron【5】等专门对这两期产品进行了对比，发现两期产品类型的空间分布存在着很大的不同。这些不同并不是全部因为地表覆盖的变化，而被解释为分类的不稳定引起的。大部分错误的变化是发生在相似的地物类型之间。特意指出的是，GlobCover 2009地表覆盖产品不能与2005产品进行直接对比，也不能用于任何变化检测。认为地表变化率比产品的分类错误率通常还要低，所以不适宜与老的地表覆盖产品对比变化。
6、MODIS地表覆盖产品。波士顿大学(BU)利用中低分辨率的MODIS卫星影像生成的地表覆盖产品MODIS Land Cover Type(MCD12Q1)，简称MLCT，产品的更新周期为每年更新，2001-2012年的每年产品有两种分辨率：5'x 5'分辨率，0.5°x 0.5°分辨率。采用监督分类的算法对每年影像进行分类得到产品。由于影像分辨率低，不同年份的产品由于受到混合像元的影响，在地物交界地带的分类结果不一致，此外还存在物候变化、火灾及干旱、虫灾等灾害而使得不同年度产品之间可比较性受到影响。MODIS collection 5产品开发了一种稳定分类结果的算法以减少年与年产品间比较时的伪变化。每个像素加了与基类有关的后验概率(posterior probability)约束值，如果像素分类结果与前一年不同，只有在新类型的后验概率值高于以前的后验概率值的时候类型才会改变。但这种方式可能造成在地表变化区域传递不正确的类型，而使得地表不能更新，因此产品生成时在连续三年的数据上进行操作，这样既能准确更新地表变化，还能使减少10％伪变化。然而，比较不同年度的地表覆盖产品MLCT得到的地表变化还是高于地表真实的变化。因此得出的结论是地表覆盖的变化通过MLCT分类后产品直接相减也就是分类后比较方法是不适宜的【6】。
7.ESA-CCI(European Space Agency-Climate Change Initiative)
第一阶段的CCI地表覆盖产品包含1998-2002,2003-2007，2008-2012三个时间段的产品。2014年开始的CCI第二阶段产品中的地表覆盖产品，目前包括全球1990至2015每年的300米分辨率地表覆盖图。分类系统采用FAO的LCCS系统。CCI-LC产品的关键在于其随着时间的推移一致性。每年的分类结果不是独立生产的而是都源于一个独特的由整个2003到2012MERIS FR和RR的影像生成的基准地表覆盖图。独立于基准地表覆盖图，利用1992至1999的AVHRR、1999至2013SPOT-VGT、2013-2015的PROBA-V时间序列影像检测地表覆盖的变化，最后一步是从基准地表覆盖图和各期变化图向后和向前推出从1992到2015绘制24年度的地表覆盖产品【7】。
由于全球变化研究的不断发展，1km以及300m的全球地表覆盖遥感数据产品已经不足以满足现在研究变化的需要，主要表现有分类体系不能满足地球系统模式需求、时相局限性较大以及全球变化研究与地球模式发展需要有更高的空间分辨率、更全面的地表覆盖信息和更符合全球变化需要的地表覆盖数据产品。
8、我国依托于863重点项目“全球地表覆盖遥感制图与关键技术研究”，研制了分辨率为30m的全球地表覆盖遥感数据产品GlobeLand30，包括2000和2010期全球十大类地物的地表覆盖产品，产品总体精度达到80％以上。【8】GlobeLand 30两期产品分别对影像进行分类获得，两期产品不建议直接比较提取地表的变化。目前，GlobeLand 30正在准备研制新一期产品，方法拟采用了利用两期影像提取地表变化，再更新产品的方式，即增量更新。增量更新(Incremental Updating)是指在进行更新操作时,只对己发生变化的目标区域迸行更新,己更新过或者未发生变化的区域则不会更新,具有更新实时快捷、处理的数据量少、便于传输的特点。
拟采用的增量更新方法CVAPS【9】首先利用对两期影像进行监督分类，进而利用分类同时得到的像元归属不同类型的后验概率分析像元发生变化的可能性。但本方法实质上还是一种先分类后比较的变化检测方法，分类本身的精度会对后续提取变化像元产生影响，由于分类的错误率往往高于地表的变化率，所以即使达到80％的分类精度也难以有效识别小于20％的地表变化。
除以上八种全球地表覆盖产品外，以下介绍一些区域地表覆盖产品及其更新的方法。
9、美国NLCD
由USGS为主导研制的NLCD(The National Land Cover Database)已经发布了4期30米分辨率的美国地表覆盖产品NLCD1992，NLCD 2001,NLCD 2006和NLCD 2011。NLCD 2006增加了从2001到2006的土地覆盖变化数据。NLCD 2006是第一个全美范围内的评价每个像素的地表覆盖变化的数据。
目前NLCD系列是五年的更新周期，以满足其成员和美国国内用户的需要。从2006版起，地
表覆盖产品的生成是利用增量更新的方法来实现。
实现影像光谱变化检测使用的算法选用了四个有互补效果的指数：NormalizedBurn Ratio归一化燃烧指数(NBR)；Normalized Difference vegetation Index归一化植被差异指数(NDVI)；Change Vector变化矢量(CV)；Relative Change Vector相对变化适量(RCV)。每个像素计算各个指数，之后两期相减计算差异。
最后一步处理是要保证新产品于前期产品趋势一致，这样才能支持与前期产品进行变化比较和叠加分析。随着新产品NLCD2011的发布，早期的产品NLCD2001、NLCD2006都进行了更新又重新发布(re-version)。更新版的早期产品修正了与新产品间的不一致性。这一步是手工完成的。【10】
10、欧洲CORINE
从上世纪80年代开始，欧洲就开始了CORINE(环境信息协调统筹)项目，建立了欧洲环境系统。其中地表覆盖产品CLC(CORINE Land Cover)利用卫星影像为主要信息源，采用自下而上的生产方式，即各成员国负责生产自己国家范围的数据产品，之后集成生成欧洲范围的地表覆盖产品。目前CLC产品包括CLC1990、CLC2000、CLC2006和CLC2012。从CLC2000版开始，出现了CLCC(Land Cover Change)产品CLCC(1990,2000)以反映10年间地表覆盖的变化。当时大部分国家的方法是首先将CLCC1990升级制作出CLC2000,再将两产品叠置产生CLCC产品。然而由于两期产品的MMU(最小制图单元)的不一致，CLCC产品中包含大量的噪声和伪变化。有些国家采用增量更新的方式，即首先检测出1990-2000的地表变化，再生成新一期的产品，这种方法后来被EEA推荐采用。
CLC从2006版CLCC产品作为基本产品之一。生产方法都是采用双时相卫星影像作为数据源(如CLCC2006采用SPOT4/5和IRS P6)，通过人工目视解译或机助半自动解译。变化图斑提取要参照顾及前期产品图斑生成以避免在生成新一期产品时产生矛盾。同时，前期产品的错误也在此过程中得到改正【11】。
11、德国DeCOVER
DeCOVER作为GMES(Global Monitoring for Environment and Security)计划的延伸，是德国服务于国内的、基于遥感影像的土地覆盖、利用系统。DeCOVER系统采用高分辨率卫星影像作为原始数据，产品比例尺大，变化检测的算法是采用基于对象的方法，通过比较两个时相的影像，并用前一时相的地表覆盖产品为约束，通过多尺度分割获得图斑，两时相影像间的diff_norm被选作为分割指标，用多尺度分割算法根据经验确定分割参数。值得一提的是，在分割过程中集成了DeCOVER前一时相的产品作为分割的约束条件，这样就保持了前一时相产品的边界不会被破坏，与分割结果及变化检测的结果边界一致。
更新过程中需要前一时相产品、两期影像参与。前期产品的边界被后期继承，错误无法被发现和改正。【12】
12、澳大利亚碳核算系统NCAS-LCCP(National Carbon Accounting-Land CoverChange Project)
利用Landsat的MSS,TM和ETM+影像进行澳大利亚25米分辨率自1972年以来的15个年份的多年生植被制图。方法是采用长时间序列的影像，监督分类各个不同年份影像估计类型参数，再通过空间/时间模型，这个模型是联合分析时间序列中的某像元的类型来减少分类错误。分类时首先较为精确地确定一个基准图像的分类，以后验概率为其他年份图像分类的依据。空间/时间模型联合分析某像元不同年份的分类结果及周围像素的类型以提高分类的精度。【13】
综上所述，目前地表覆盖产品更新方法主要包括各期影像分别分类的方法和增量更新的方法。实践证明增量更新的方法可获得更高的精度以及使各期产品可以相互对比提取变化。增量的提取是通过影像变化检测的方法。基于遥感影像的变化检测方法主要包括基于像素和基于对象两类方法【14】。基于像素的变化检测算法以像元作为基本单位，利用像元光谱特征判断变化的发生，并不考虑空间特征的变化。使用基于像元的变化检测方法会不可避免的产生“椒盐”噪声，而基于对象变化检测是以图像分割为基础，将同质性像元集合最为一个对象，确定最佳分割尺度来进行变化检测。
基于对象变化检测的图斑结果是以矢量形式表现的，这种形式相对于基于像素的栅格形式来说，有效避免了“椒盐”噪声，变化图斑的边界明确，一些零散图斑相对于栅格形式的零散像元来说更加容易处理。由于基于对象方法所产生的对象与各时相的图像特征有关，对象的几何特征会随着时间的推移而改编，所以两个时相对象边界会存在不一致，难以建立多时相对象之间的对应关系【15】。用栅格形式表达图斑，其准确度相较于矢量形式来说更高，更能表达出一些较为微小的变化。基于对象变化检测的图斑结果是以多边形矢量形式表现的，但是对于一些非人工、分布零散的地物如森林、灌木、草地来说，不同时期的多边形图斑叠加往往难以准确反映地表变化范围。而用栅格形式表达图斑更为灵活，其准确度相较于矢量形式来说更高，更能表达出一些较为微小的变化。
基于像元和基于对象的方法均需要依赖某一指数或某一特征值的阈值来进行变化区域与非变化区域的评判标准，阈值的选取会直接影响结果的精度。
国外学者Rother等【16】在2006年首次提出了协同分割的概念，在基于马尔科夫随机场的能量函数中同时包含了空间一致性和图像对共同部分直方图匹配的全局约束。通过置信域图割的方法实现能量函数的最小化，使得分割结果达到最优。优点是具有广泛的普遍性，但是该模型能量函数的优化只考虑了一个图像对的情况。Joulin等则将协同分割的应用拓展到分割一组图像，其考虑到非监督分类方法能够将一张图像分割为前景和背景区域，将归一化割和对象识别中的内核方法相结合，提出了基于判别聚类框架的图像协同分割方法，其目标是使用一个训练好的监督分类器在一组图像上标记前景/背景标签，以使图像达到最优分割【17】。为了实现协同分割能够处理大量数据的目标，Joulin等在判别聚类框架的基础上，将光谱聚类和判别聚类在能量函数中结合，利用能量最小化方法，将协同分割方法拓展到能够处理多个类别，同时分割图像的数量也有明显提高【18】。
我国学者孟凡满等【19】考虑提出的多组协同分割框架，不仅能够发现每一个图像组内部的图像协同信息，而且能够在不同的图像组中传递这些信息，来获得更准确的对象先验知识。袁敏等【15】在基于图割的能量函数中包含不同时相图像的变化信息，将变化检测与图像分割联系在一起，在能量函数的最小化的过程中，同时完成图像分割与变化检测。这种方法将运动图像的协同分割与地表覆盖的变化过程联系起来，拓展了协同分割的应用领域。
协同分割思想的提出将图像分割转化为能量函数的最小化问题，将不再依赖阈值的选取。同时协同分割的变化检测方法能直接得到边界准确、空间对应的多时相变化对象，解决了多时相对象边界不一致问题。
相比于基于像素以及基于对象这两种变化检测方法，协同分割方法的结果和基于像素方法的结果一样是栅格的形式，但由于算法中与周边像素的关系而能够很好的避免椒盐噪声的出现，集成了栅格与矢量数据两种类型的特点。
以往的地表覆盖分类和变化检测多是从遥感影像出发，分析其光谱、形状、纹理等特征因子，但遥感影像反映的只是地表的瞬时状态，而且均是基于像素级别的，所以存在很多错误及不确定因素。针对全球范围内复杂多样的地类，存在大量同物异谱和异物同谱的现象，同时由于季候时相的原因，在变化检测的过程中都会存在大量的伪变化。比如在灌溉期，水田由于灌水灌溉，会与水体有非常相似的光谱特征，如果单从影像的光谱特征来做变化检测的话很容易判断此处发生了变化，但实际这属于一种伪变化。如果想要得到高精度的地表分类及变化检测结果，仅仅从光谱上来判断是很难的。想要提高变化检测的精度，就要对检测后的变化图斑进行伪变化识别与剔除。
目前伪变化一般是利用人工来识别，不仅浪费时间、耗费大量资源，而且由于识别人员经验不足很容易造成错判。伪变化识别最好的办法就是请遥感行业专家来识别和判断，利用领域专家的知识，目前地表覆盖变化检测需要寻找一种机器替代专家来判断的可行方法，即建立专家系统。专家系统构建中最重要的就是规则库，要建立规则库来辅助变化检测，就需要找到可用的与地表覆盖相关的辅助数据及相关知识，其中生态地理分区由于其全球性、分区内部地类稳定性、地物变化规律性和信息量大等特点，可以用来构建规则库辅助变化检测。
目前全球范围的遥感变化检测规则库技术还没有人展开过研究，对生态地理分区的应用也仅限于利用某些特定区域的生态分区，而且以对植被生长和变化的判断为主，很少应用于全球地表覆盖变化检测中，也没有相应的系统被整体归纳、总结与开发。
生态地理区域系统反映了包括气候条件、地形条件、水分条件、土壤和植被等在内的自然要素的空间格局，并且体现了这些条件与资源环境的匹配关系。目前对于全球范围变化检测结果的自动判断工作还没有人展开过研究，对生态地理分区的应用也仅限于利用某些特定区域的分区，而且以对植被生长和变化的判断为主，很少有对全球生态地理分区进行整体归纳与总结。
发明内容
本发明的目的是提供一种基于影像协同分割与生态地理分区规则库的地表覆盖产品增量更新方法，以解决利用多时相卫星影像，采用协同分割技术提取地表覆盖增量，即地表覆盖的变化，之后对变化图斑进行分类，再结合全球生态地理分区规则库对分类增量图斑进行去伪变化，提高增量检测结果的可靠度，最后用增量图斑更新旧时相地表覆盖产品，获得新的地表覆盖产品的技术问题。
为了实现上述发明目的，本发明所述的技术方案如下：
一种基于影像协同分割与生态地理分区规则库的地表覆盖产品增量更新方法，包括协同分割提取两期影像的变化图斑；利用生态地理分区规则库识别剔除伪变化；地表覆盖产品的更新；具体步骤如下：
一、多时相遥感影像协同分割提取两期影像的变化图斑；
1)求得变化强度图
为了保证变化强度图的质量，需将两期原始遥感影像进行常规方法预处理(几何预处理、辐射预处理)；
变化强度图以图像的光谱特征为基础，变化强度值(CVpixel)是在n维空间中分别求取不同时相同一像元的光谱值之差的绝对值后再求平均，即：
Figure BDA0001537012810000071
式中，
Figure BDA0001537012810000072
和
Figure BDA0001537012810000073
为像元(p,q)于时间t1、t2，在波段k的光谱值，k＝1,2,…,n，n为图像的波段数量；
2)求得变化特征项
以变化强度图和设定的阈值条件为基础，构建变化特征项E1；
E1＝∑p∈PDp(lp) (2)
式中P为全部像元的集合；Dp(lp)表示将像元p分配为目标和背景的代价；其中lp∈{“obj”，“bkg”}，“obj”表示目标(变化区域)，“bkg”表示背景(不变区域)；
Figure BDA0001537012810000074
Figure BDA0001537012810000075
式中，t为设定的阈值，在本发明中t为距离变化强度图的均值(μ)若干倍的标准差(σ)处，即t＝μ+nσ，其中n为常数，同时保证t为正数；
3)求得图像特征项
首先需要构建综合图像特征项的函数(公式5)，由两个时相的光谱和纹理特征组成的图像特征以及权重两部分组成:
E2＝λt1Et1+(1-λt1)Et2 (5)
上式中Et1和Et2分别表示前一时相(t1)的图像特征和后一时相(t2)的图像特征；λt1表示t1时相的图像特征的权重，(1-λt1)表示t2时相的图像特征项的权重；
其中
Eti＝∑{p，q∈N}V{p，q} (6)
Figure BDA0001537012810000081
上式中Eti表示ti(i＝1,2)时相的图像特征，V{p，q}是像元p、q不连续的代价；
Figure BDA0001537012810000083
为像元邻域；Ip、Iq为像元p、q的特征；(p,q)是p和q之间的欧式距离；σ2＝<||Ip-Iq||2>为归一化因子，<>表示在整幅图像上求均值；
如果|Ip-Iq|<σ时，表示相邻像元之间的差异较小，则V{p，q}的值会相对较大；如果|Ip-Iq|>σ，相邻像元差异较大，V{p，q}的值会相对较小；
图像特征选用光谱特征和纹理特征，即
Eti＝λSEtiS+(1-λS)EtiT (8)
式中，EtiS表示ti时相图像的光谱特征，EtiT表示ti时相图像的纹理特征；λs和(1-λs)分别表示光谱特征和纹理特征的权重，两个权重的取值范围均为[0,1]；
①计算光谱特征
对经过预处理的图像，首先根据公式(9)计算得到单时相图像的光谱特征，在本发明中采用所有光谱波段的均值来表示光谱特征，
Figure BDA0001537012810000082
式中，Ik(p)为像元p对应的第k波段的光谱值；n为所采用的光谱波段的总数量；；N为所采用的光谱波段的最后一个波段数
②计算纹理特征
对两时相图像分别进行主成分分析，得到第一主成分，然后采用等概率量化算法进行灰度级化简，计算得到四个纹理特征，分别为基于灰度共生矩阵的采用常规算法计算的角二阶矩(ASM)、对比度(CON)、相关性(COR)和熵(ENT)，即
ET＝λASMEASM+λCONECON+λCORECOR+λENTEENT (10)
式中的λ为各个纹理特征的权重，四个纹理特征的权重之和应为1。
4)构建能量函数
变化检测能量函数形式为：
Energy_CD(l)＝λE1+E2 (11)
E1表示变化特征项，E2表示综合图像特征项，λ为变化特征的权重。完成纹理和光谱特征计算后即用公式(11)求得最终综合图像特征；
5)最小割/最大流方法求最小割
运用图论的原理，将图像的每个像素作为普通节点，设置两个特殊节点源点s和汇点t，以此来构成网络流图，再将计算完成的能量函数Energy_CD，作为构建的网络流图中的容量；寻找图中的最大流，图中的最大流使得图中的一组边达到饱和状态，从而将图中的普通节点分成相互独立的两个子集
Figure BDA0001537012810000091
和
Figure BDA0001537012810000092
即为目标和背景，达到图像的最优分割结果，这组达到饱和的边的集合为图的最小割；
求取最大流的方法采用广度优先搜索，利用路径中边的数量表示距离，即两点间路径的长度，点的距离表示该节点到源点的最短路径的长度；根据图中普通节点到源点的最短路径的长度，用广度优先搜索方法将图中的普通节点分为不同的子集，当汇点t进入到某一子集中后，分层结束；分层结束后，根据分层结果搜索最短路径，然后在最短路径中进行流量更新；汇点到源点的初识距离为2，每次循环汇点到源点的长度+1，以此类推进行计算直到找到图中的最大流/最小割；
求得最大流/最小割后，对分割结果进行二值化，完成图像的协同分割算法；
6)并行计算
协同分割算法对于大尺寸图像设计了并行计算的方法；利用MATLAB的parallelcomputing toolbox中的parfor并行方法，对for语句进行并行计算；
首先将大尺寸的图像裁切成数个小尺寸的图像，用matlabpool开启多个并行计算池，再用parfor语句一次读取多个图像同时进行计算和存储。完成所有小图像的计算后，再将这些小图像按照其原始图像的位置进行拼接，得到最后的大尺寸图像的结果；
7)变化图斑像素分类
利用常规的监督分类方法对协同分割提取的变化像素相应的两期影像像素进行地物分类；
二、利用生态地理分区规则库识别剔除伪变化；
在协同分割提取变化像素基础上，基于生态地理规则库来识别剔除伪变化；
具体包括：
1)生态地理分区规则库框架设计
采用国际生态学界普遍认同的世界基金组织为自然保护目的建立的全球生态分区作为全球生态地理分区规则库的基础框架，该生态分区把全球分为8个生态地理分区和14个生物群落；采用面向对象的框架形式，设计规则库的基本组成部分，此规则库自上而下建立不同层间的派生和继承关系，以此来构造规则库；本规则库采用左右两支并行的方式，按照左右分别向下一层展开，前三层左右互不交叉，到最后一层按照各分区对应的地理位置、生物群落和各种自然属性继承其父类的规则；第一层，左支设计为生态地理类，右支设计为自然地理类，本层是TOP层，所有其他对象都是它们两个的子类；第二层，左支为8个地理分区和14个生物群落作为生态地理类的子类，右支设计为海拔、坡度、NDVI、温度、水分；第三层，左支为大地理生态分区，即地理分区和生物群落的交叉继承类，不同分区存在的生物群落种类不同，由于每一个大地理生态分区范围过大，此处均取生态分区的平均值；本层规则库分别存储其相应的伪变化规则；第四层为最底层，是由小地理生态分区组成，分别存储属性信息及本小分区内特有的伪变化规则，及由第三层继承而来的伪变化规则。
2)规则的建立和表达
伪变化规则的表达方式采用产生式方式，即前提和结论的形式并映射为数据库中的表格，采用常见的对象-关系型数据库存储；
按照生态地理分区规则库框架结构，逐层建立伪变化规则；规则库始于框架模型的第三层，左支大地理生态分区各分区的伪变化规则，主要包括本分区伪变化规则和季节时相原因造成的伪变化规则；伪变化规则采用6位代码XXXXXX表示，前3位XXX为变化前类别代码，后3位XXX为变化后类别代码，代码采用表1中GlobeLand 30地物类型代码表示，各分区伪变化规则表如表2，伪变化规则存储在pseudochange字段；季节时相原因造成的伪变化规则如表3，各时期伪变化规则存储在pseudochange字段；
第三层右支：不同属性，即高度、坡度、NDVI值、温度和水分条件下存在的伪变化规则，表格的设计包括ID号、属性类型和伪变化规则字段；
分区规则库模型中的第四层为各小生态分区伪变化规则，包括从上层继承的规则和不能继承的特有规则；各分区内不存在继承关系的特有规则，主要包括本小分区伪变化规则和季节时相原因造成的伪变化规则；
3)规则库管理系统的功能设计
作为可视化显示、对规则库中的规则表进行查询、增加、删除和修改操作，以及实现伪变化识别标记和剔除的系统，伪变化判断方法采用正向推理，首先将协同分割提取出的变化像素栅格数据利用常用软件(如ArcGIS软件中的“栅格转面”工具)按照地物类型代码转换为矢量图斑，转换时输出矢量边缘与输入栅格的单元保持一致。矢量图斑属性需存储地物类型代码。将变化图斑根据图斑变化前后地表覆盖类型表达为6位代码，前3位XXX为变化前类别代码，后3位XXX为变化后类别代码，根据地图坐标范围判断图斑所在小生态地理分区，调用相应生态分区的伪变化规则库，即分区规则库模型第四层，包括从第三层继承的伪变化规则和不能继承的特有规则，即本小分区伪变化规则和季节时相原因造成的伪变化规则，判断6位代码是否与规则库中匹配，若是则为伪变化，否则为真变化，标记、剔除伪变化图斑；由于变化图斑可能跨越不同的生态分区，则需将图斑对应的全部分区逐个判断，实现结合生态地理分区数据库的去伪技术；去除伪变化图斑后，需将最终的变化图斑利用常用软件(如ArcGIS软件中的“面转栅格”工具)转换为栅格像素，转换时输出像素与输入矢量图斑边缘保持一致。像素存储地物类型代码。
三、地表覆盖产品更新
地表覆盖数据一般采用栅格形式存储，更新操作为利用经过上述步骤得到的变化图斑地表覆盖类型替换旧时相产品相应位置像素类型以生成新一期的地表覆盖产品。
本发明的优点与积极效果：
使用本发明提出的基于协同分割的遥感影像变化检测算法，可以提取出同一地区不同时相的遥感影像上的变化图斑，改变了基于像元和基于对象的变化检测方法中以阈值条件作为变化和非变化的唯一标准，有效的避免了椒盐噪声的出现。
本发明针对原始的利用遥感影像分类和变化检测方法存在的问题，提出了一种基于生态地理分区规则库来实现伪变化识别的方法。
生态地理分区规则库能实现的功能有：
(1)能够使用空间关系数据库加载数据，并可以实现查询结果的可视化。实现数据在数据库和系统之间的导入与导出，使数据能够在本系统中进行编辑、查询和删除等功能。
(2)能够以不同的比例尺、不同的范围显示生态分区矢量图、栅格和tif格式等影像图和辅助数据。可以实现图形在显示区的缩放、移动和漫游等。
(3)提供多种查询功能。包括属性查询、条件查询和点查询。进行空间查询之后，选择的生态分区在图形显示区被高亮显示，其相对应的属性表显示在下方。
(4)能够实现输入变化图斑后自动判断生态分区，并连接数据库调取属性表及相应伪变化规则，可以实现伪变化的识别和标注。数据库的连接直接写在了后台程序里，在系统界面不做展示，当执行的操作需要调用数据库时，再与数据库进行连接，如果只是对图层进行操作则默认不连接数据库。
附图说明
图1是本发明的增量更新技术路线图。
图2是本发明的生态分区规则库框架图。
图3是本发明的伪变化推理流程图。
具体实施方式
本发明所述的一种基于影像协同分割与生态地理分区规则库的地表覆盖产品增量更新方法，如图1所示，包括协同分割提取两期影像的变化图斑；利用生态地理分区规则库识别剔除伪变化；地表覆盖产品的更新；具体步骤如下：
一、多时相遥感影像协同分割技术
1)求得变化强度图
为了保证变化强度图的质量，需将两期原始遥感影像进行常规方法预处理(几何预处理、辐射预处理)。
变化强度图以图像的光谱特征为基础，变化强度值(CVpixel)是在n维空间中分别求取不同时相同一像元的光谱值之差的绝对值后再求平均，即：
Figure BDA0001537012810000121
式中，
Figure BDA0001537012810000122
和
Figure BDA0001537012810000123
为像元(p,q)于时间t1、t2，在波段k的光谱值，k＝1,2,…,n，n为图像的波段数量；
2)求得变化特征项
以变化强度图和设定的阈值条件为基础，构建变化特征项E1；
E1＝∑p∈PDp(lp) (2)
式中P为全部像元的集合；Dp(lp)表示将像元p分配为目标和背景的代价；其中lp∈{“obj”，“bkg”}，“obj”表示目标(变化区域)，“bkg”表示背景(不变区域)；
Figure BDA0001537012810000124
Figure BDA0001537012810000125
式中，t为设定的阈值，在本发明中t为距离变化强度图的均值(μ)若干倍的标准差(σ)处，即t＝μ+nσ，其中n为常数，同时保证t为正数，在本发明实验中采用了n＝+0.5；
3)求得图像特征项
首先需要构建综合图像特征项的函数(公式5)，由两个时相的光谱和纹理特征组成的图像特征以及权重两部分组成:
E2＝λt1Et1+(1-λt1)Et2 (5)
上式中Et1和Et2分别表示前一时相(t1)的图像特征和后一时相(t2)的图像特征；λt1表示t1时相的图像特征的权重，(1-λt1)表示t2时相的图像特征项的权重；
其中
Eti＝∑{p，q∈N}V{p，q} (6)
Figure BDA0001537012810000131
上式中Eti表示ti(i＝1,2)时相的图像特征，V{p，q}是像元p、q不连续的代价；
Figure BDA0001537012810000133
为像元邻域；Ip、Iq为像元p、q的特征；(p,q)是p和q之间的欧式距离；σ2＝<||Ip-Iq||2>为归一化因子，<>表示在整幅图像上求均值；
如果|Ip-Iq|<σ时，表示相邻像元之间的差异较小，则V{p，q}的值会相对较大；如果|Ip-Iq|>σ，相邻像元差异较大，V{p，q}的值会相对较小；
图像特征选用光谱特征和纹理特征，即
Eti＝λSEtiS+(1-λS)EtiT (8)
式中，EtiS表示ti时相图像的光谱特征，EtiT表示ti时相图像的纹理特征；λs和(1-λs)分别表示光谱特征和纹理特征的权重，两个权重的取值范围均为[0,1]；
①计算光谱特征
对经过预处理的图像，首先根据下式计算得到单时相图像的光谱特征，在本发明中采用所有光谱波段的均值来表示光谱特征，
Figure BDA0001537012810000132
式中，Ik(p)为像元p对应的第k波段的光谱值；n为所采用的光谱波段的总数量；；N为所采用的光谱波段的最后一个波段数
②计算纹理特征
对两时相图像分别进行主成分分析，得到第一主成分，然后采用等概率量化算法进行灰度级化简，计算得到四个纹理特征，分别为基于灰度共生矩阵的采用常规算法计算的角二阶矩(ASM)、对比度(CON)、相关性(COR)和熵(ENT)，即
ET＝λASMEASM+λCONECON+λCORECOR+λENTEENT (10)
式中的λ为各个纹理特征的权重，四个纹理特征的权重之和应为1。
4)构建能量函数
变化检测能量函数形式为：
Energy_CD(l)＝λE1+E2 (11)
E1表示变化特征项，E2表示综合图像特征项，λ为变化特征的权重。完成纹理和光谱特征计算后即用公式(11)求得最终综合图像特征；
5)最小割/最大流方法求最小割
运用图论的原理，将图像的每个像素作为普通节点，设置两个特殊节点源点s和汇点t，以此来构成网络流图，再将计算完成的能量函数Energy_CD，作为构建的网络流图中的容量；寻找图中的最大流，图中的最大流使得图中的一组边达到饱和状态，从而将图中的普通节点分成相互独立的两个子集
Figure BDA0001537012810000141
和
Figure BDA0001537012810000142
即为目标和背景，达到图像的最优分割结果，这组达到饱和的边的集合为图的最小割；
求取最大流的方法采用广度优先搜索，利用路径中边的数量表示距离，即两点间路径的长度，点的距离表示该节点到源点的最短路径的长度；根据图中普通节点到源点的最短路径的长度，用广度优先搜索方法将图中的普通节点分为不同的子集，当汇点t进入到某一子集中后，分层结束；分层结束后，根据分层结果搜索最短路径，然后在最短路径中进行流量更新；汇点到源点的初识距离为2，每次循环汇点到源点的长度+1，以此类推进行计算直到找到图中的最大流/最小割；
求得最大流/最小割后，对分割结果进行二值化，完成图像的协同分割算法；
6)并行计算
协同分割算法对于大尺寸图像设计了并行计算的方法；利用MATLAB的parallelcomputing toolbox中的parfor并行方法，对for语句进行并行计算。
首先将大尺寸的图像裁切成数个小尺寸的图像，用matlabpool开启多个并行计算池，再用parfor语句一次读取多个图像同时进行计算和存储。完成所有小图像的计算后，再将这些小图像按照其原始图像的位置进行拼接，得到最后的大尺寸图像的结果。
7)变化图斑像素分类
利用常规的监督分类方法对协同分割提取的变化像素相应的两期影像像素进行地物分类。本专利以我国生产的GlobeLand 30地表覆盖产品分类体系为例，如表1所示。用户可根据具体应用需要选择其他分类体系。
二、基于生态地理分区规则库的变化图斑的去伪技术；
在协同分割提取变化像素基础上，基于生态地理规则库来识别剔除伪变化的方法以进一步提高精度；
具体包括：
1)生态地理分区规则库框架设计
采用国际生态学界普遍认同的世界基金组织为自然保护目的建立的全球生态分区作为全球生态地理分区规则库的基础框架，该生态分区把全球分为8个生态地理分区和14个生物群落。基于这两个基础图层，全球共划分867个生态区。8大地理分区为欧亚大陆地区(PA)、南美洲地区(NT)、东亚东南亚地区(IM)、太平洋地区(OC)、北美洲地区(NA)、非洲地区(AT)、大洋洲地区(AA)、南极洲地区(AN)；14个生物群落为热带及亚热带湿润阔叶林(01)、热带和亚热带干旱阔叶林(02)、热带亚热带针叶森林(03)、温带阔叶混交林(04)、温带针叶林(05)、寒带森林\针叶林(06)、热带和亚热带草原，稀树草原和灌丛(07)、温带草原，稀树草原和灌丛(08)、淹水的草原和稀树草原(09)、山地草原和灌丛(10)、苔原(11)、地中海森林，林地和灌丛(12)、沙漠与旱生灌丛(13)、红树林(14)。
采用面向对象的框架形式，设计规则库的基本组成部分，见图2，在8个生态地理分区和14个生物群落的基础上，对每一分区考虑了温度、降水量、NDVI、高程、坡度等属性因素；此规则库采用面向对象的方式，自上而下建立不同层间的派生和继承关系，以此来构造规则库；本规则库采用左右两支并行的方式，按照左右分别向下一层展开，前三层左右互不交叉，到最后一层按照各分区对应的地理位置、生物群落和各种自然属性继承其父类的规则。第一层，左支设计为生态地理类，右支设计为自然地理类，本层是TOP层，所有其他对象都是它们两个的子类。第二层，左支为8个地理分区和14个生物群落作为生态地理类的子类，右支设计为海拔、坡度、NDVI、温度、水分。第三层，左支为大地理生态分区，即地理分区和生物群落的交叉继承类，不同分区存在的生物群落种类不同，本层左支一共有64个，命名为IM01…PA01等。右支按照不同生态分区的自然条件不同，把高程(E)、坡度(S)、NDVI值(N)、温度(T)和水分(W)分别按照其不同值来划分。高程分为：E1低海拔(<2000m)、E2中海拔(2000m-4000m)、E3亚高海拔(4000m-6000m)、E4高海拔(6000m-7000m)和E5级高海拔(>7000m)；坡度分为：S1低坡度(<5°)、S2中坡度(5°-30°)和S3大坡度(>30°)；温度分为：T1热带、T2亚热带、T3温带、T4寒带，按公知的纬度划分；水分分为：W1干旱区(降雨量<200mm)、W2半干旱区(降雨量200mm-400mm)、W3半湿润区(降雨量400mm-800mm)、W4湿润区(降雨量>800mm)；NDVI值分为：N1(-1-0)、N2(0)、N3(0-1)。由于每一个大地理生态分区范围过大，此处均取生态分区的平均值。本层规则库分别存储其相应的伪变化规则。第四层为最底层，是由867个小地理生态分区组成，分别存储属性信息及本小分区内特有的伪变化规则，及由第三层继承而来的伪变化规则。
2)规则的建立和表达
伪变化规则的表达方式采用产生式方式，即前提和结论的形式并映射为数据库中的表格，采用常见的对象-关系型数据库存储。
按照生态地理分区规则库框架结构，逐层建立伪变化规则。规则库始于框架模型的第三层，左64个大地理生态分区各分区的伪变化规则，主要包括本分区伪变化规则和季节时相原因造成的伪变化规则。伪变化规则采用6位代码XXXXXX表示，前3位XXX为变化前类别代码，后3位XXX为变化后类别代码，代码采用表1中GlobeLand 30地物类型代码表示，耕地为010、林地为020、草地为030、灌木为040、湿地为050、水体为060、苔原为070、人造代表为080、裸地为090、冰川和永久积雪为100。各分区伪变化规则表如表2，伪变化规则存储在pseudochange字段。参见表2所示。
季节时相原因造成的伪变化规则分为ST1代表丰水期；ST2表示枯水期；ST3表示灌溉期；ST4表示收割期；ST5表示植物生长期；ST6表示植物落叶期；ST7表示冰冻期。各时期伪变化规则存储在pseudochange字段，如表3。
第三层右支：不同属性(高度、坡度、NDVI值、温度和水分)条件下存在的伪变化规则。表的设计包括ID号、属性名称、属性的取值范围和伪变化规则字段。表4为不同高程条件下的伪变化规则表，高程分为ID、elevation、pseudochange三个字段，ID字段名为1、2、3、4、5，对应elevation分为E1、E2、E3、E4、E5五个等级，pseudochange字段中存储伪变化规则。表5为不同坡度条件下的伪变化规则表，坡度分为ID、slope、pseudochange三个字段，ID字段名为1、2、3，对应slope分为S1、S2、S3三个等级，pseudochange字段中存储伪变化规则。表6为不同NDVI值条件下的伪变化规则表，NDVI值分为ID、NDVI、pseudochange三个字段，ID字段名为1、2、3、4、5，对应NDVI值分为N1、N2、N3、N4、N5五个等级，pseudochange字段中存储伪变化规则。表7为不同温度条件下的伪变化规则表，共ID、temperature、pseudochange三个字段，ID字段名为1、2、3、4，对应温度值分为T1、T2、T3、T4四个等级，pseudochange字段中存储伪变化规则。表8为不同水分条件下的伪变化规则表，水分分为ID、wet、pseudochange三个字段，ID字段名为1、2、3、4，对应wet分为W1、W2、W3、W4四个等级，pseudochange字段中存储伪变化规则。
分区规则库模型中的第四层为各小生态分区伪变化规则，包括从上层继承的规则和不能继承的特有规则。各分区内不存在继承关系的特有规则，主要包括本小分区伪变化规则和季节时相原因造成的伪变化规则。
各小分区特有伪变化规则仍采用表格2的方式存储。季节时相原因造成的伪变化规则存储方式与表3相同。
各表格pseudochange字段的具体伪变化规则6位代码需要用户根据应用需求采集行业专家知识自行填入。本发明提供一些示例，参见表9、10、11、12、13所示。
3)规则库管理系统的功能设计
规则库管理系统主要作用是规则库可视化显示、对规则库中的规则表进行查询、增加、删除和修改操作，以及实现伪变化识别标记和剔除。伪变化判断方法采用正向推理，首先将协同分割提取出的变化像素栅格数据利用常用软件(如ArcGIS软件中的“栅格转面”工具)按照地物类型代码转换为矢量图斑，转换时输出矢量边缘与输入栅格的单元保持一致。矢量图斑属性需存储地物类型代码。如图3，将变化图斑根据图斑变化前后地表覆盖类型表达为6位代码，前3位XXX为变化前类别代码，后3位XXX为变化后类别代码，根据地图坐标范围判断图斑所在小生态地理分区，调用相应生态分区的伪变化规则库，即分区规则库模型第四层，包括从第三层继承的伪变化规则和不能继承的特有规则(本小分区伪变化规则和季节时相原因造成的伪变化规则)。判断6位代码是否与规则库中匹配，若是则为伪变化，否则为真变化，标记、剔除伪变化图斑。由于变化图斑可能跨越不同的生态分区，则需将图斑对应的全部分区逐个判断，实现结合生态地理分区数据库的去伪技术。去除伪变化图斑后，需将最终的变化图斑利用常用软件(如ArcGIS软件中的“面转栅格”工具)转换为栅格像素，转换时输出像素与输入矢量图斑边缘保持一致。像素存储地物类型代码。
三、地表覆盖产品更新
地表覆盖数据一般采用栅格形式存储，更新操作为利用经过上述步骤得到的变化图斑地表覆盖类型替换旧时相产品相应位置像素类型以生成新一期的地表覆盖产品。表1：GlobeLand 30地表覆盖系统一级分类体系
Figure BDA0001537012810000181
表2:各分区伪变化规则
id Eco-region name pseudochange
1 IM01
…
64 PA01
表3:季侯时相对应的伪变化规则
id season_time pseudochange
1 ST1
2 ST2
3 ST3
4 ST4
5 ST5
6 ST6
7 ST7
表4：高程对应的伪变化规则
id elevation pseudochange
1 E1
2 E2
3 E3
4 E4
5 E5
表5：坡度对应的伪变化规则
id slope pseudochange
1 S1
2 S2
3 S3
表6：NDVI值对应的伪变化规则
id ndvi pseudochange
1 N1
2 N2
3 N3
4 N4
5 N5
表7：温度对应的伪变化规则
id temperature pseudochange
1 T1
2 T2
3 T3
4 T4
表8：水分对应的伪变化规则
id wet pseudochange
1 W1
2 W2
3 W3
4 W4
表9：高程对应的伪变化规则示例
Figure BDA0001537012810000201
Figure BDA0001537012810000211
表10：坡度对应的伪变化规则示例
Figure BDA0001537012810000212
表11：NDVI值对应的伪变化规则示例
id ndvi pseudochange
1 N1 020050 030060 040080 050090 010020 040070
2 N2 050060 040050 050080 020060 010040 020050
3 N3 010020 040050 080090 060050 060020 060030
4 N4 070010 070050 070030 050080 020090 080050
5 N5 020050 040050 090060 080040 070050 090050
表12：温度对应的伪变化规则示例
Figure BDA0001537012810000213
Figure BDA0001537012810000221
表13：水分对应的伪变化规则示例
Figure BDA0001537012810000222
参考文献
【1】T.R.Loveland,B.C.Reed,J.F.Brown,D.O.Ohlen,Z.Zhu,L.Yang &J.，W.Merchant(2000)Development of a global land cover characteristics databaseand IGBP DISCover from 1km AVHRR data,International Journal of RemoteSensing,21:6-7,1303-1330,
DOI:10.1080/014311600210191
【2】Hansen MC,DeFries RS,Townshend JRG,Sohlberg R(2000)Global landcover classification at 1km spatial resolution using a classification treeapproach.Int J Remote Sens 21(6/7):1331–1364
【3】Bartholomé,E.,&Belward,A.S.(2005).GLC2000:a new approach to globalland cover mapping from Earth observation data.International Journal ofRemote Sensing,26,1959-1977
【4】Tateishi R,Uriyangqai B,Al-Bilbisi H,Ghar MA,Tsend-Ayush J,Kobayashi T,Kasimu A et al(2011)Production of global land cover data–GLCNMO.Int J Digit Earth 4(1):22–49
【5】Sophie Bontemps,Pierre Defourny,Eric Van Bogaert,Olivier Arino,Vasileios Kalogirou,Jose Ramos Perez，GLOBCOVER 2009Products Description andValidation Report，2011
【6】Friedl,M.A.,Sulla-Menashe,D.,Tan,B.,Schneider,A.,Ramankutty,N.,Sibley,A.,&Huang,X.(2010).MODIS Collection 5global land cover:Algorithmrefinements and characterization of new datasets.Remote Sensing ofEnvironment,114,168-182
【7】Land Cover CCI PRODUCT USER GUIDE VERSION 2.0，
http://maps.elie.ucl.ac.be/CCI/viewer/index.php
【8】Chen Jun,Yifang Ban，Songnian Li，China:Open access to Earth land-cover map，Nature 514,434(23October 2014)doi:10.1038/514434c published online22October 2014
【9】陈军，陈晋，廖安平等，全球地表覆盖遥感制图，科学出版社，2016.211-228.
【10】Jin,S.,Yang,L.,Danielson,P.,Homer,C.,Fry,J.,&Xian,G.(2013).Acomprehensive change detection method for updating the National Land CoverDatabase to circa 2011.Remote Sensing of Environment,132,159-175
【11】Büttner,G.(2014).CORINE Land Cover and Land Cover ChangeProducts,18,55-74
【12】Buck,O.(2010).DeCOVER 2–the German GMES extension to support landcover systems:status and outlook[13]Mora,B.,Tsendbazar,N.-E.,Herold,M.,&Arino,O.(2014).Global Land Cover Mapping:Current Status and Future Trends,18,11-30【13】Caccetta,P.A.,Furby,S.L.,O'Connell,J.,Wallace,J.F.,&Wu,X.(2007).Continental monitoring:34years of land cover change using Landsatimagery.32nd International Symposium on Remote Sensing of Environment,June25–29,2007,San José,Costa Rica.
【14】Masroor,H.,et al.,Change detection from remotely sensed images:From pixel-based to object-based approaches.ISPRS Journal of Photogrammetryand Remote Sensing,2013.80:p.91-106.
【15】袁敏,肖鹏峰,冯学智,等.基于协同分割的高分辨率遥感图像变化检测[J].南京大学学报自然科学,2015(5):1039-1048.
【16】Rother C,Minka T,Blake A,et al.Cosegmentation of Image Pairs byHistogram Matching-Incorporating a Global Constraint into MRFs[C]//ComputerVision and Pattern Recognition,2006IEEE Computer Society Conference on.IEEEXplore,2006:993-1000.
【17】Joulin A,Bach F,Ponce J.Discriminative clustering for image co-segmentation[C]//Computer Vision and Pattern Recognition.IEEE,2010:1943-1950.
【18】Joulin.A.,Multi-class cosegmentation[C]//IEEE Conference onComputer Vision and Pattern Recognition.IEEE Computer Society,2012:542-549.
【19】Meng F,Cai J,Li H.Cosegmentation of multiple image groups[J].Computer Vision&Image Understanding,2016,146:67-76.

Claims (1)

1.一种基于影像协同分割与生态地理分区规则库的地表覆盖产品增量更新方法，包括协同分割提取两期影像的变化图斑；利用生态地理分区规则库识别剔除伪变化；地表覆盖产品的更新；具体步骤如下：
一、多时相遥感影像协同分割提取两期影像的变化图斑；
1)求得变化强度图
将两期原始遥感影像用常规方法进行预处理，包括几何预处理、辐射预处理；
变化强度图以图像的光谱特征为基础，变化强度值(CVpixel)是在n维空间中分别求取不同时相同一像元的光谱值之差的绝对值后再求平均，即：
Figure FDA0002696079640000011
式中，BVpqk(t1)和BVpqk(t2)为像元(p，q)于时间t1、t2，在波段k的光谱值，k＝1，2，…，n，n为图像的波段数量；
2)求得变化特征项
以变化强度图和设定的阈值条件为基础，构建变化特征项E1；
E1＝∑p∈PDp(lp) (2)
式中P为全部像元的集合；Dp(lp)表示将像元p分配为目标和背景的代价；其中lp∈{“obj”，“bkg”}，“obj”表示目标变化区域，“bkg”表示背景不变区域；
Figure FDA0002696079640000012
Figure FDA0002696079640000013
式中，t为设定的阈值，在本发明中t为距离变化强度图的均值(μ)若干倍的标准差(σ)处，即t＝μ+nσ，其中n为常数，同时保证t为正数；
3)求得图像特征项
首先需要构建综合图像特征项的函数E2，由两个时相的光谱和纹理特征组成的图像特征以及权重两部分组成：
E2＝λt1Et1+(1-λt1)Et2 (5)
上式中Et1和Et2分别表示前一时相(t1)的图像特征和后一时相(t2)的图像特征；λt1表示t1时相的图像特征的权重，(1-λt1)表示t2时相的图像特征项的权重；
其中
Eti＝∑{p，q∈N}V{p，q} (6)
Figure FDA0002696079640000021
上式中Eti表示ti(i＝1，2)时相的图像特征，V{p，q}是像元p、q不连续的代价；
Figure FDA0002696079640000023
为像元邻域；Ip、Iq为像元p、q的特征；(p，q)是p和q之间的欧式距离；σ2＝<||Ip-Iq||2>为归一化因子，<>表示在整幅图像上求均值；
如果|Ip-Iq|＜σ时，表示相邻像元之间的差异较小，则V{p，q}的值会相对较大；如果|Ip-Iq|＞σ，相邻像元差异较大，V{p，q}的值会相对较小；
图像特征选用光谱特征和纹理特征，即
Eti＝λsEtis+(1-λs)EtiT (8)
式中，EtiS表示ti时相图像的光谱特征，EtiT表示ti时相图像的纹理特征；λs和(1-λS)分别表示光谱特征和纹理特征的权重，两个权重的取值范围均为[0，1]；
①计算光谱特征
对经过预处理的图像，首先根据公式(9)计算得到单时相图像的光谱特征，在本发明中采用所有光谱波段的均值来表示光谱特征，
Figure FDA0002696079640000022
式中，Ik(p)为像元p对应的第k波段的光谱值；n为所采用的光谱波段的总数量；N为所采用的光谱波段的最后一个波段数；
②计算纹理特征
对两时相图像分别进行主成分分析，得到第一主成分，然后采用等概率量化算法进行灰度级化简，计算得到四个纹理特征，分别为基于灰度共生矩阵的常规算法计算的角二阶矩(ASM)、对比度(CON)、相关性(COR)和熵(ENT)，即
ET＝λASMEASM+λCONECON+λCORECOR+λENTEENT (10)
式中的λ为各个纹理特征的权重，四个纹理特征的权重之和应为1；
4)构建能量函数
变化检测能量函数形式为：
Energy_CD(l)＝λE1+E2 (11)
E1表示变化特征项，E2表示综合图像特征项，λ为变化特征的权重；完成纹理和光谱特征计算后即用公式(11)求得最终综合图像特征；
5)最小割/最大流方法求最小割
运用图论的原理，将图像的每个像素作为普通节点，设置两个特殊节点源点s和汇点t，以此来构成网络流图，再将计算完成的能量函数Energy_CD，作为构建的网络流图中的容量；寻找图中的最大流，图中的最大流使得图中的一组边达到饱和状态，从而将图中的普通节点分成相互独立的两个子集ν1和ν2，即为目标和背景，达到图像的最优分割结果，这组达到饱和的边的集合为图的最小割；
求取最大流的方法采用广度优先搜索，利用路径中边的数量表示距离，即两点间路径的长度，点的距离表示该节点到源点的最短路径的长度；根据图中普通节点到源点的最短路径的长度，用广度优先搜索方法将图中的普通节点分为不同的子集，当汇点t进入到某一子集中后，分层结束；分层结束后，根据分层结果搜索最短路径，然后在最短路径中进行流量更新；汇点到源点的初始距离为2，每次循环汇点到源点的长度+1，以此类推进行计算直到找到图中的最大流/最小割；
求得最大流/最小割后，对分割结果进行二值化，完成图像的协同分割算法；
6)并行计算
协同分割算法对于大尺寸图像设计了并行计算的方法；利用MATLAB的parallelcomputing toolbox中的parfor并行方法，对for语句进行并行计算；
首先将大尺寸的图像裁切成数个小尺寸的图像，用matlabpool开启多个并行计算池，再用parfor语句一次读取多个图像同时进行计算和存储；完成所有小图像的计算后，再将这些小图像按照其原始图像的位置进行拼接，得到最后的大尺寸图像的结果；
7)变化图斑像素分类
利用常规的监督分类方法对协同分割提取的变化像素相应的两期影像像素进行地物分类；
二、利用生态地理分区规则库识别剔除伪变化；
在协同分割提取变化像素基础上，基于生态地理规则库来识别剔除伪变化；
具体包括：
1)生态地理分区规则库框架设计
采用国际生态学界普遍认同的世界基金组织为自然保护目的建立的全球生态分区作为全球生态地理分区规则库的基础框架，该生态分区把全球分为8个生态地理分区和14个生物群落；采用面向对象的框架形式，设计规则库的基本组成部分，此规则库自上而下建立不同层间的派生和继承关系，以此来构造规则库；本规则库采用左右两支并行的方式，按照左右分别向下一层展开，前三层左右互不交叉，到最后一层按照各分区对应的地理位置、生物群落和各种自然属性继承其父类的规则；第一层，左支设计为生态地理类，右支设计为自然地理类，本层是TOP层，所有其他对象都是它们两个的子类；第二层，左支为8个地理分区和14个生物群落作为生态地理类的子类，右支设计为海拔、坡度、NDVI、温度、水分；第三层，左支为大地理生态分区，即地理分区和生物群落的交叉继承类，不同分区存在的生物群落种类不同，由于每一个大地理生态分区范围过大，此处均取生态分区的平均值；本层规则库分别存储其相应的伪变化规则；第四层为最底层，是由小地理生态分区组成，分别存储属性信息及本小分区内特有的伪变化规则，及由第三层继承而来的伪变化规则；
2)规则的建立和表达
伪变化规则的表达方式采用产生式方式，即前提和结论的形式并映射为数据库中的表格，采用常见的对象-关系型数据库存储；
按照生态地理分区规则库框架结构，逐层建立伪变化规则；规则库始于框架模型的第三层，左支大地理生态分区各分区的伪变化规则，主要包括本分区伪变化规则和季节时相原因造成的伪变化规则；伪变化规则采用6位代码XXXXXX表示，前3位XXX为变化前类别代码，后3位XXX为变化后类别代码，代码采用GlobeLand 30地物类型代码表示，各分区伪变化规则存储在pseudochange字段；
各时期伪变化规则存储在pseudochange字段；
第三层右支：不同属性，即高度、坡度、NDVI值、温度和水分条件下存在的伪变化规则，表格的设计包括ID号、属性类型和伪变化规则字段；
分区规则库模型中的第四层为各小生态分区伪变化规则，包括从上层继承的规则和不能继承的特有规则；各分区内不存在继承关系的特有规则，主要包括本小分区伪变化规则和季节时相原因造成的伪变化规则；
3)规则库管理系统的功能设计
作为可视化显示、对规则库中的规则表进行查询、增加、删除和修改操作，以及实现伪变化识别标记和剔除的系统，伪变化判断方法采用正向推理，首先将协同分割提取出的变化像素栅格数据按照地物类型代码转换为矢量图斑，转换时输出矢量边缘与输入栅格的单元保持一致；矢量图斑属性需存储地物类型代码；将变化图斑根据图斑变化前后地表覆盖类型表达为6位代码，前3位XXX为变化前类别代码，后3位XXX为变化后类别代码，根据地图坐标范围判断图斑所在小生态地理分区，调用相应生态分区的伪变化规则库，即分区规则库模型第四层，包括从第三层继承的伪变化规则和不能继承的特有规则，即本小分区伪变化规则和季节时相原因造成的伪变化规则，判断6位代码是否与规则库中匹配，若是则为伪变化，否则为真变化，标记、剔除伪变化图斑；由于变化图斑可能跨越不同的生态分区，则需将图斑对应的全部分区逐个判断，实现结合生态地理分区数据库的去伪技术；去除伪变化图斑后，需将最终的变化图斑转换为栅格像素，转换时输出像素与输入矢量图斑边缘保持一致；像素存储地物类型代码；
三、地表覆盖产品更新
地表覆盖数据采用栅格形式存储，更新操作为利用经过上述步骤得到的变化图斑像素地表覆盖类型替换旧时相产品相应位置像素类型以生成新一期的地表覆盖产品。
CN201810001776.3A 2018-01-02 2018-01-02 一种基于影像协同分割与生态地理分区规则库的地表覆盖产品增量更新方法 Expired - Fee Related CN108089850B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201810001776.3A CN108089850B (zh) 2018-01-02 2018-01-02 一种基于影像协同分割与生态地理分区规则库的地表覆盖产品增量更新方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201810001776.3A CN108089850B (zh) 2018-01-02 2018-01-02 一种基于影像协同分割与生态地理分区规则库的地表覆盖产品增量更新方法

Publications (2)

Publication Number Publication Date
CN108089850A CN108089850A (zh) 2018-05-29
CN108089850B true CN108089850B (zh) 2020-11-24

Family

ID=62181541

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201810001776.3A Expired - Fee Related CN108089850B (zh) 2018-01-02 2018-01-02 一种基于影像协同分割与生态地理分区规则库的地表覆盖产品增量更新方法

Country Status (1)

Country Link
CN (1) CN108089850B (zh)

Families Citing this family (10)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109934154B (zh) * 2019-03-08 2021-06-01 北京科技大学 一种遥感影像变化检测方法及检测装置
CN110837875B (zh) * 2019-11-18 2022-07-05 国家基础地理信息中心 地表覆盖数据质量异常判断方法及装置
CN111611968B (zh) * 2020-05-29 2022-02-01 中国科学院西北生态环境资源研究院 一种遥感图像的处理方法以及遥感图像处理模型
CN111860204A (zh) * 2020-06-29 2020-10-30 成都数之联科技有限公司 基于语义分割技术的多时相遥感影像变化检测方法和介质
CN111897907B (zh) * 2020-08-12 2024-03-01 国家基础地理信息中心 地表覆盖数据伪变化信息过滤方法、装置及存储介质
CN113487483A (zh) * 2021-07-05 2021-10-08 上海商汤智能科技有限公司 影像分割网络的训练方法和装置
CN114003755B (zh) * 2021-10-25 2022-10-11 中国自然资源航空物探遥感中心 多源卫星分景影像数据组织存储与检索方法、系统及设备
CN114119575B (zh) * 2021-11-30 2022-07-19 二十一世纪空间技术应用股份有限公司 一种空间信息变化检测方法及系统
CN114972440B (zh) * 2022-06-21 2024-03-08 江西省国土空间调查规划研究院 用于国土调查的es数据库图斑对象链式追踪方法
CN115147726B (zh) * 2022-09-05 2023-03-24 清华大学 城市形态图的生成方法、装置、电子设备和可读存储介质

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US6571242B1 (en) * 2000-07-25 2003-05-27 Verizon Laboratories Inc. Methods and systems for updating a land use and land cover map using postal records
CN104851113A (zh) * 2015-04-17 2015-08-19 华中农业大学 多分辨率遥感影像的城市植被自动提取方法
CN107273813A (zh) * 2017-05-23 2017-10-20 国家地理空间信息中心 基于高分卫星遥感数据的地理空间要素提取系统

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US6571242B1 (en) * 2000-07-25 2003-05-27 Verizon Laboratories Inc. Methods and systems for updating a land use and land cover map using postal records
CN104851113A (zh) * 2015-04-17 2015-08-19 华中农业大学 多分辨率遥感影像的城市植被自动提取方法
CN107273813A (zh) * 2017-05-23 2017-10-20 国家地理空间信息中心 基于高分卫星遥感数据的地理空间要素提取系统

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
"地表覆盖遥感产品更新完善的研究动向";陈军 等;《遥感学报》;20160925;第991-1001页 *

Also Published As

Publication number Publication date
CN108089850A (zh) 2018-05-29

Similar Documents

Publication Publication Date Title
CN108089850B (zh) 一种基于影像协同分割与生态地理分区规则库的地表覆盖产品增量更新方法
Hermosilla et al. Land cover classification in an era of big and open data: Optimizing localized implementation and training data selection to improve mapping outcomes
Khandelwal et al. An approach for global monitoring of surface water extent variations in reservoirs using MODIS data
Coulter et al. Classification and assessment of land cover and land use change in southern Ghana using dense stacks of Landsat 7 ETM+ imagery
Zhang et al. Object-oriented method for urban vegetation mapping using IKONOS imagery
Guan et al. Integration of orthoimagery and lidar data for object-based urban thematic mapping using random forests
Lucas et al. The earth observation data for habitat monitoring (EODHaM) system
Ardila et al. Markov-random-field-based super-resolution mapping for identification of urban trees in VHR images
Pérez-Hoyos et al. A methodology to generate a synergetic land-cover map by fusion of different land-cover products
Schiewe Segmentation of high-resolution remotely sensed data-concepts, applications and problems
Phinn et al. Monitoring the composition of urban environments based on the vegetation-impervious surface-soil (VIS) model by subpixel analysis techniques
Pérez-Hoyos et al. Conventional and fuzzy comparisons of large scale land cover products: Application to CORINE, GLC2000, MODIS and GlobCover in Europe
Nabil et al. Assessing factors impacting the spatial discrepancy of remote sensing based cropland products: A case study in Africa
Kiptala et al. Land use and land cover classification using phenological variability from MODIS vegetation in the Upper Pangani River Basin, Eastern Africa
Chen et al. Detecting changes in high-resolution satellite coastal imagery using an image object detection approach
Salih Classification and mapping of land cover types and attributes in Al-Ahsaa Oasis, Eastern Region, Saudi Arabia using Landsat-7 data
Toosi et al. Citrus orchard mapping in Juybar, Iran: Analysis of NDVI time series and feature fusion of multi-source satellite imageries
Guo et al. Mozambique flood (2019) caused by tropical cyclone idai monitored from sentinel-1 and sentinel-2 images
Zhu et al. Exploiting cosegmentation and Geo-Eco zoning for land cover product updating
Qu et al. Mapping large area tea plantations using progressive random forest and Google Earth Engine
Danoedoro et al. Preliminary study on the use of digital surface models for estimating vegetation cover density in a mountainous area
Korting GEODMA: a toolbox integrating data mining with object-based and multi-temporal analysis of satellite remotely sensed imagery
Mudau et al. Towards development of a national human settlement layer using high resolution imagery: a contribution to SDG reporting
Sedona et al. An automatic approach for the production of a time series of consistent land-cover maps based on long-short term memory
OAP OBJECT· BASED RE-CLASSIFICATION OF HIGH RESOLUTION DIGITAL IMAGERY FOR URBAN LAND-USE MONITORING

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
GR01 Patent grant
GR01 Patent grant
CF01 Termination of patent right due to non-payment of annual fee

Granted publication date: 20201124

Termination date: 20220102

CF01 Termination of patent right due to non-payment of annual fee