[发明专利]一种基于影像协同分割与生态地理分区规则库的地表覆盖产品增量更新方法

权利要求书

1.一种基于影像协同分割与生态地理分区规则库的地表覆盖产品增量更新方法，包括协同分割提取两期影像的变化图斑；利用生态地理分区规则库识别剔除伪变化；地表覆盖产品的更新；具体步骤如下：

一、多时相遥感影像协同分割提取两期影像的变化图斑；

1)求得变化强度图

将两期原始遥感影像用常规方法进行预处理，包括几何预处理、辐射预处理；

变化强度图以图像的光谱特征为基础，变化强度值(CV_pixel)是在n维空间中分别求取不同时相同一像元的光谱值之差的绝对值后再求平均，即：

式中，BV_pqk(t₁)和BV_pqk(t₂)为像元(p，q)于时间t₁、t₂，在波段k的光谱值，k＝1，2，…，n，n为图像的波段数量；

2)求得变化特征项

以变化强度图和设定的阈值条件为基础，构建变化特征项E₁；

E₁＝∑_p∈PD_p(l_p) (2)

式中P为全部像元的集合；D_p(l_p)表示将像元p分配为目标和背景的代价；其中l_p∈{“obj”，“bkg”}，“obj”表示目标变化区域，“bkg”表示背景不变区域；

式中，t为设定的阈值，在本发明中t为距离变化强度图的均值(μ)若干倍的标准差(σ)处，即t＝μ+nσ，其中n为常数，同时保证t为正数；

3)求得图像特征项

首先需要构建综合图像特征项的函数E2，由两个时相的光谱和纹理特征组成的图像特征以及权重两部分组成：

E₂＝λ_t1E_t1+(1-λ_t1)E_t2 (5)

上式中E_t1和E_t2分别表示前一时相(t₁)的图像特征和后一时相(t₂)的图像特征；λ_t1表示t₁时相的图像特征的权重，(1-λ_t1)表示t₂时相的图像特征项的权重；

其中

E_ti＝∑_{p，q∈N}V_{p，q} (6)

上式中E_ti表示t_i(i＝1，2)时相的图像特征，V_{p，q}是像元p、q不连续的代价；为像元邻域；I_p、I_q为像元p、q的特征；(p，q)是p和q之间的欧式距离；σ²＝||I_p-I_q||²为归一化因子，表示在整幅图像上求均值；

如果|I_p-I_q|＜σ时，表示相邻像元之间的差异较小，则V_{p，q}的值会相对较大；如果|I_p-I_q|＞σ，相邻像元差异较大，V_{p，q}的值会相对较小；

图像特征选用光谱特征和纹理特征，即

E_ti＝λ_sE_tis+(1-λ_s)E_tiT (8)

式中，E_tiS表示t_i时相图像的光谱特征，E_tiT表示t_i时相图像的纹理特征；λ_s和(1-λ_S)分别表示光谱特征和纹理特征的权重，两个权重的取值范围均为[0，1]；

①计算光谱特征

对经过预处理的图像，首先根据公式(9)计算得到单时相图像的光谱特征，在本发明中采用所有光谱波段的均值来表示光谱特征，

式中，I_k(p)为像元p对应的第k波段的光谱值；n为所采用的光谱波段的总数量；N为所采用的光谱波段的最后一个波段数；

②计算纹理特征

对两时相图像分别进行主成分分析，得到第一主成分，然后采用等概率量化算法进行灰度级化简，计算得到四个纹理特征，分别为基于灰度共生矩阵的常规算法计算的角二阶矩(ASM)、对比度(CON)、相关性(COR)和熵(ENT)，即

E_T＝λ_ASME_ASM+λ_CONE_CON+λ_CORE_COR+λ_ENTE_ENT (10)

式中的λ为各个纹理特征的权重，四个纹理特征的权重之和应为1；

4)构建能量函数

变化检测能量函数形式为：

Energy_CD(l)＝λE₁+E₂ (11)

E₁表示变化特征项，E₂表示综合图像特征项，λ为变化特征的权重；完成纹理和光谱特征计算后即用公式(11)求得最终综合图像特征；

5)最小割/最大流方法求最小割

运用图论的原理，将图像的每个像素作为普通节点，设置两个特殊节点源点s和汇点t，以此来构成网络流图，再将计算完成的能量函数Energy_CD，作为构建的网络流图中的容量；寻找图中的最大流，图中的最大流使得图中的一组边达到饱和状态，从而将图中的普通节点分成相互独立的两个子集ν₁和ν₂，即为目标和背景，达到图像的最优分割结果，这组达到饱和的边的集合为图的最小割；

求取最大流的方法采用广度优先搜索，利用路径中边的数量表示距离，即两点间路径的长度，点的距离表示该节点到源点的最短路径的长度；根据图中普通节点到源点的最短路径的长度，用广度优先搜索方法将图中的普通节点分为不同的子集，当汇点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位代码是否与规则库中匹配，若是则为伪变化，否则为真变化，标记、剔除伪变化图斑；由于变化图斑可能跨越不同的生态分区，则需将图斑对应的全部分区逐个判断，实现结合生态地理分区数据库的去伪技术；去除伪变化图斑后，需将最终的变化图斑转换为栅格像素，转换时输出像素与输入矢量图斑边缘保持一致；像素存储地物类型代码；

三、地表覆盖产品更新

地表覆盖数据采用栅格形式存储，更新操作为利用经过上述步骤得到的变化图斑像素地表覆盖类型替换旧时相产品相应位置像素类型以生成新一期的地表覆盖产品。