CN103400151B - 一体化的光学遥感影像与gis自动配准与水体提取方法 - Google Patents

Abstract

一体化的光学遥感影像与GIS自动配准与水体提取方法：步骤1，对输入的光学遥感影像进行分割，获得水系初始分割结果；步骤2，进行局部水平集演化分割得到集合R；步骤3，匹配基础地理信息水系层矢量对象与集合R中对象，配准影像，如果配准成功则进入步骤4，否则返回步骤2；步骤4，通过缓冲区检测获取未变化水体、疑似新增水体、疑似变化水体，进一步滤除疑似水体对象，确认真正变化水体与未变化水体；步骤5，利用未变化水体对象中对应像元的多光谱值作为正样本，随机选择光学遥感影像上集合R以外区域中的像元多光谱值作为负样本，训练SVM分类器并验证过滤后的疑似水体对象集中各对象是否为水体，获得分割配准提取的水体结果。

Description

一体化的光学遥感影像与GIS自动配准与水体提取方法

技术领域

本发明是关于光学影像与矢量的异源数据配准及水体提取技术，具体涉及一体化的光学遥感影像与GIS自动配准与水体提取方法。

背景技术

水体分布信息的准确获取对于洪水监测与灾害评估、地理信息更新、水资源调查、流域综合治理、水利规划等领域具有重要的意义。由于遥感具有时效性、经济性等特点，使之成为提取水体最广泛、最有效的手段之一。

近年来高分辨率影像的迅猛发展给遥感信息提取提供了更多更好的数据支持，但同时也带来了巨大挑战。由于在遥感领域“同谱异物，同物异谱”的现象广泛存在，影像上存在有与水体表现特征类似的各种地物，如阴影、浓密植被、沥青路面等，加之受光照条件、地形因素、复杂环境、影像质量等影响，使得从遥感影像上自动提取水体信息是一个非常困难的问题。另一方面，在灾害发生或者其他应急响应情况下，水体信息提取在时效性上有很高的要求，现有的水体信息提取方法难以满足实际应用的需求。

现有的遥感影像水体提取方法主要分为以下四类：特征分割方法、模式分类方法、轮廓线检测方法以及多源数据综合分析方法。特征分割方法又包括自上而下的知识驱动和自下而上的数据驱动两类，知识驱动型依据地物特征构建先验模型和知识规则来指导分割，如水平集演化方法，数据驱动型根据数据自身的灰度分布特征进行分割，如单波段阈值法、特征指数法、色彩空间转换法等，其中特征指数法是通过多光谱波段运算获取反映水体与背景反差的指数，进而利用直方图阈值分割得到提取结果，典型的水体特征指数如归一化水体指数(NDWI)、改进的归一化差异水体指数(MNDWI)、新型水体指数(NWI)等，其模型简单，并且有不错的提取结果。模式分类方法首先对原图像进行特征提取，构造特征向量，然后利用神经网络、支撑向量机(SVM)等机器学习方法将原始图像中的每个像素点分类为水体与背景两类对象，从而实现对水体的提取。轮廓线检测方法通过边缘检测获得水体的岸线边缘，然后对边缘进行编组获取水体的主体区域，最后采用纹理跟踪、区域生长等方法对检测到的水体区域进行合并，该类方法由于充分考虑的地物的边缘特征，其提取结果的轮廓定位精度较高。多源数据综合分析方法充分利用水体在不同数据源中的特征，如SAR(合成孔径雷达)、GIS(地理信息系统)数据等，通过多源数据之间的相互验证确定最终水体提取结果，比依靠单一数据源的方法有更高的正确率。尽管以上方法在试验中取得了较好的提取结果，但是面对实际应用的需求时依然存在以下不足：第一，很多方法均是针对特定的试验数据取得了较好的效果，缺少对各类数据处理的普适性；第二，提取方法模型中的参数需要手动调节，或者模型训练样本需要人工选择，使得提取算法的自动化程度偏低；第三，对于多源数据综合分析的方法，其前提假设是多源数据已经过严格的配准，而实际处理过程中所面对的多源数据均是存在位置和坐标系的差异，需要手工的选择控制点并执行几何纠正处理。

综上可以看出，面对水体在遥感影像上的多样化表现和复杂背景干扰时，现有的水体提取方法都存在普适性差、人工干预多等问题。

发明内容

针对上述问题，本发明提出了一体化的光学遥感影像与GIS自动配准与水体提取方法。

本发明提供一体化的光学遥感影像与GIS自动配准与水体提取方法，包括以下步骤：

步骤1，对输入的光学遥感影像进行分割，获得水系初始分割结果；

步骤2，以每个水系初始分割结果中的每个对象为当前待处理对象构建局部缓冲区，进行局部水平集演化分割，得到新的分割结果，构成的分割结果对象集合记为R；

步骤3，基于形状曲线相似度匹配GIS数据中水系层已有的基础地理信息矢量数据集合V与步骤2所得集合R中对象，以GIS数据为基准配准光学遥感影像，如果配准成功则进入步骤4，如果配准不成功则将当前的分割结果中每个对象作为当前待处理对象返回步骤2，重新构建局部缓冲区并进行局部水平集演化分割；

步骤4，对经步骤3配准后的V与R做基于缓冲区的变化检测，得到未变化水体R_w0、疑似变化水体R_w1和疑似新增水体R_nw；然后分别对疑似变化水体R_w1和疑似新增水体R_nw做进一步的验证如下，

对疑似变化水体R_w1做基于骨架线相似性的检测，确定真正变化的水体和未变化的水体，未变化的水体并入未变化水体R_w0，得到集合R‘_w0，R_w1中保留的变化水体记为集合R‘_w1；

对疑似新增水体R_nw进行两重验证，包括首先将R_nw与GIS数据中道路和建筑层对象做叠加分析，拓扑相交的对象则非水体对象并从R_nw中滤除；然后对R_nw中剩余对象进行上下文纹理特征过滤，包括提取R_nw中对象周边背景纹理特征，并与疑似水体纹理特征做比较，如果纹理相似度超过预设纹理特征相似度阈值，则认为该疑似水体对象为阴影的半影区域并从R_nw中滤除，过滤后剩余的疑似水体对象构成集合R‘_nw；

步骤5，利用集合R‘_w0中对应像元的多光谱值作为正样本，随机选择光学遥感影像上集合R以外区域中的像元多光谱值作为负样本，训练SVM分类器，并采用训练所得SVM分类器验证步骤4过滤后剩余的集合R“_nw中各对象是否为水体，在R‘_nw中保留验证结果为水体的对象，从R‘_nw中取出验证结果为非水体的对象，最终保留在R‘_nw中的对象记为集合R“_nw，以集合R“_nw为分割配准提取的水体结果。

而且，步骤1中，判断输入的光学遥感影像是否具有NIR波段，有则采用NDWI分割，否则采用OTSU分割。

而且，步骤2实现方式为，设步骤1所得水系初始分割结果中的所有对象构成一个集合，以集合中的任一对象为当前待处理对象，进行以下步骤，

(a)对当前待处理对象，按预设缓冲距离d构建缓冲区，得到待处理影像区域；

(b)以(a)得到的待处理影像区域为数据域，以当前待处理对象为初始零水平集，表示为，其中(x,y)为像素点坐标，t为时间参数，设在曲线内部曲线外部引入Heaviside函数和一维Dirac函数建立模型的能量函数如下，

其中，Ω为待处理影像区域，u₁，u₂分别为零水平集内、外的平均灰度值，u为待处理影像区域某一点灰度值，为当前点梯度的模，固定参数α₁,α₂＞0,β≥0,γ≥0；令u₁，u₂保持不变，使能量函数最小，得到其对应的偏微分方程：

根据式(5)，采用C-V模型水平集演化方法分割提取水系的对象。

而且，步骤3中包括如下子步骤，

(a)设基础地理信息水系层矢量对象的集合为V＝{v₁，v₂，v₃，...v_m}，其中v_i为集合中的一个对象，i＝1,2...m，m为集合中对象总数；设步骤2所得集合R＝{r₁，r₂，r₃，...r_n}，其中r_j为集合中的一个对象，j＝1,2...n，n为集合中对象总数；

取出V中每一对象分别计算与R中各对象的形状相似度，生成一个二维形状相似度矩阵S_vr[m，n]，设V中某对象对应多边形A，R中某对象对应多边形B，按以下形状曲线相似度计算公式计算对象之间的形状相似度，

$\mathrm{pCor}\_\mathrm{shape}(A,B)=\max \left(\frac{\underset{k=1}{\overset{n_p}{\mathrm{\Σ}}}({d_k}^A-\stackrel{\‾}{d^A})(d_{\mathrm{lk}}^B-\stackrel{\‾}{d_l^B})}{\sqrt{{\underset{k=1}{\overset{n_p}{\mathrm{\Σ}}}(d_k^A-\stackrel{\‾}{d^A})}^2\underset{k=1}{\overset{n_p}{\mathrm{\Σ}}}{(d_{\mathrm{lk}}^B-\stackrel{\‾}{d_l^B})}^2}}\right),l=\mathrm{1,2}...n_p---\left(7\right)$

其中，为多边形A外边界的第k采样点到质心的距离，为所有的均值，表示多边形B由第l_k点到质心的距离，其中l_k＝l+k,l+k≤n_p; l_k＝l+k-n_p,l+k＞n_p，为所有的均值，n_p为多边形边界离散化后的点数，多边形A与B的边界点数相同，设其中某一个点为起始点，采用l和k标识相对于起始点的点对应的序号；

(b)将集合V中各对象分别对应的多边形按面积排序并从中依次选择面积最大的三个多边形，记为{v_i1，v_i2，v_i3}；

(c)从集合R中各对象分别对应的多边形中，依次选择形状曲线相似度最大的多边形与{v_i1，v_i2，v_i3}匹配，记为{r_j1，r_j2，r_j3}；

(d)计算多边形v_i1，r_j1，v_i2，r_j2，v_i3，r_j3对应的质心点在图像上的横纵坐标值

{(x_vi1，y_vi1)，(x_rj1，y_rj1);(x_vi2，y_vi2)，(x_rj2，y_rj2);(x_vi3，y_vi3)，(x_rj3，y_rj3)}

采用下式求解几何变换参数(a₀，a₁，a₂，b₀，b₁，b₂)，得到几何变换模型；

$\left\{\begin{array}{c}{x}_v=a_0+a_1\×x_r+a_2\×y_r\\ y_v=b_0+b_1\×x_r+b_2\×y_r\end{array}\right.---\left(8\right)$

其中，x_v,y_v,x_r,y_r为自变量，取值分别为

(x_vi1，y_vi1)，(x_rj1，y_rj1);(x_vi2，y_vi2)，(x_rj2，y_rj2);(x_vi3，y_vi3)，(x_rj3，y_rj3)；

(e)从集合V中剩余的m-3个对象中依次取一个对象，计算所取对象质心(x_vi，y_vi)与集合R中剩余n-3个对象质心基于式(8)变换后坐标的距离，设距离最小的质心为(x_tri，y_tri)，如果质心(x_vi，y_vi)与(x_tri，y_tri)的距离小于预设阈值D_thr，则将这两个质心作为一对匹配点，相应的两个对象为匹配对象；遍历完集合V中所有对象，计算成功匹配对象数量占集合V中对象总数的比例，超过预设阈值P则认为当前几何变换模型是稳定的；否则返回(c)，重新根据形状曲线相似度依次选择集合R中三个对象与集合V中已选出的对象匹配；如果遍历了集合R中所有多边形均无法与集合V中选出的多边形对应，则返回(b)，从集合V中按面积依次选择其他三个多边形，再次迭代与R中的多边形匹配，直到获取到稳定的几何变换模型；获取到稳定的几何变换模型后以GIS数据为基准配准光学遥感影像，配准成功；如果遍历结合V和R中所有对象仍无法获得稳定的几何变换模型，则配准不成功。

而且，步骤4中对经步骤3配准后的V与R做基于缓冲区的变化检测，得到未变化水体R_w0、疑似变化水体R_w1、疑似新增水体R_nw，包括如下子步骤，

(a)对V中对象以预设缓冲距离R_buffer构造缓冲区V_buffer；

(b)将R与V_buffer做叠加，基于式(9)分别以R中各对象为当前待检查的对象，统计在缓冲区外的长度占总长度的比例P_Line，

$P_{\mathrm{Line}}=\frac{L_{\mathrm{out}\_\mathrm{buffer}}}{L_{\mathrm{total}}}---\left(9\right)$

其中，L_{out_buffer}为当前待检查的对象边缘位于缓冲区外部分的长度累加和，L_total为当前待检查的对象边缘总长度，P_Line为当前待检查的对象边缘位于GIS对象边缘缓冲区外部分长度占总长度的比例。

(c)当P_Line大于预设阈值时，认为当前待检查的对象为疑似变化水体R_w1，当P_Line小于等于预设阈值时的时候，认为当前待检查的对象为未变化水体R_w0；R中与V在位置上不重叠的对象为疑似新增水体R_nw。

而且，步骤4中对R_nw中剩余对象进行上下文纹理特征过滤的实现方式如下，

(a)从R_nw的剩余对象中选择一个对象r_nwi，在其质心处构建一个矩形模板RECT_i，RECT_i的长度小于r_nwi最小外接矩形的宽度；

(b)根据预设步长绕质心旋转模板RECT_i，同时统计每次旋转后模板内像素的方差，得到关于对象r_nwi的角度-方差图，记角度-方差图中方差最小值对应的方向为r_nwi的纹理主方向Direct_i，对应的矩形模板为RECT_imin；

(c)以r_nwi质心为原点，Direct_i和-Direct_i为方向构造射线Line_i1和Line_i2，射线与r_nwi的边缘有两个交点Pt_i1,Pt_i2，构造以Pt_i1,Pt_i2为端点，Direct_i,-Direct_i为方向，RECT_i的宽度为宽，长宽比为5的矩形对象Rect_i1,Rect_i2；

(d)分别统计Rect_i1,Rect_i2,RECT_imin对应区域像素的灰度共生矩阵生成的自相关统计量，记为COR_i1,COR_i2,COR_imin；

(e)如果自相关统计量差值|COR_i1-COR_imin|,|COR_i2-COR_imin|中有一个值小于预设纹理特征相似度阈值，则认为对象r_nwi与周边纹理相近，r_nwi为阴影区域，从R_nw中剔除该对象；如果|COR_i1-COR_imin|,|COR_i2-COR_imin|均大于等于预设纹理特征相似度阈值，则保留该对象；

(f)返回(a)依次对R_nw中的所有剩余对象执行(a)-(e)，直到对R_nw中的所有剩余对象处理完毕。

本发明提出的方法提供一体化的光学遥感影像与GIS自动配准与水体提取方法，该方法基于NDWI或OTSU对水系进行初始分割，以初始分割结果作为初始轮廓进行局部水平集演化分割，对矢量对象与影像分割对象分别计算形状曲线相似度并排序匹配，如果匹配失败，则将当前演化结果作为初始轮廓再次进行局部水平集演化；通过迭代获取到精确匹配结果后，同时也获得水系提取结果，首先利用纹理特征对疑似水体做初始过滤，然后利用光谱特征训练SVM分类器，然后对水系提取结果做进一步的验证并获取最终水系提取结果。与传统的方法相比，本发明具有以下特点：

1)针对光学影像和水域矢量数据，提出了一种配准与目标地物提取的一体化方法，配准的精度与目标地物的提取精度相互影响，通过迭代实现不断精确化的过程。

2)提出基于改进的形状曲线相似度排序匹配方法，针对待匹配数据对象数量较小，但同名矢量多边形对象存在整体的相似性，同时顾及多边形对象之间相对位置约束，获取正确的匹配结果。

附图说明

图1为本发明实施例的流程图。

具体实施方式

通过分析水体在影像上的光谱、纹理、形状、上下文等特征以及水体变化规律可知，水体在遥感影像上具有显著的视觉特征，容易辨识；其次，水体具有恒稳性，即水域一般不会消失，即使受季节性影响范围和轮廓会发生一定的变化。基于这些特性，我们提出利用GIS矢量辅助全自动提取水体的方法。其核心思路是分割、配准、变化检测和提取的一体化，通过水平集迭代演化分割方法和改进形状曲线相似性方法实现水体遥感影像与GIS的自动配准，同时获得配准和水体分割结果，配准成功的水体可作为未变化的水体，未配准的水体可基于GIS矢量利用水平集演化得到精确水体边界，对于新增的水体及非水体对象，则利用已知水体来学习训练验证未知水体，同时考虑对象的上下文纹理特征过滤非水体对象。

本发明技术方案可采用计算机软件技术实现自动运行流程。以下结合附图和实施例详细说明本发明技术方案。如图1，实施例的技术方案的流程包括以下步骤：

步骤1，对输入的光学遥感影像分割，获得水系初始分割结果。可根据输入的光学遥感影像是否具有近红外波段分别采用归一化水指数(NormalDifferentialWaterIndex，简记NDWI)阈值分割和最大类间方差(OTSU)分割，获得水系初始分割结果。OTSU分割针对基于视觉显著图定位后的水体影像区域。

实施例中，先判断输入的光学遥感影像是否具有NIR波段，有则采用NDWI分割，否则采用OTSU分割。

现有技术中NDWI的计算表达式为：

NDWI＝(p(Green)-p(NIR))/(p(Green)+p(NIR))(1)

其中，p(Green)为绿波段灰度值，p(NIR)为近红外灰度值。

计算每个像素对应的NDWI值，生成一幅指数影像，然后取阈值为0对指数影像进行二值分割，影像中像元值大于0的像元构成的即为水体初步提取结果。

可采用Itti和Koch提出的视觉注意模型实现影像中水体显著区域定位，该模型中影像分解采用9种尺度，在每个尺度下分别获取1个亮度分量、4个颜色分量、4个方向分量的特征。

现有技术中OTSU是一种自适应的阈值确定方法，是按图像的灰度特性，将通过视觉注意模型确定的各个可能存在水体的影像块分别分成背景和目标两部分，使背景和目标之间的类间方差最大。其计算过程是一个遍历过程，计算灰度范围内的所有灰度值作为阈值时对应的类间方差，选择类间方差最大时对应的灰度值作为最终的阈值，以此对影像进行阈值分割，小于此阈值的像素构成的分割结果为水体初步提取结果；类间方差计算表达式为：

g＝w₀w₁(μ₀-μ₁)²(2)

其中，w₀为前景像素点数占整幅图像的比例，μ₀为其平均灰度值；w₁为背景像素点数占整幅图像的比例，μ₁为其平均灰度值。

步骤2，对水系初始分割结果中的每个对象构建局部缓冲区，进行局部水平集演化分割，具体实现方式如下，

设步骤1所得水系初始分割结果中的所有对象构成一个集合，即初始的水系分割结果对象集合，然后对初始水系分割结果对象集合中的每个对象分别执行以下操作，得到新的分割对象集合。

设初始水系分割结果对象集合中的任一对象为当前待处理对象，进行以下处理：

(a)对当前待处理对象，以预设缓冲距离d构建局部缓冲区，使得缓冲区外围轮廓所包围的影像区域(即待处理影像区域)中包括单一水系和背景；建议设缓冲距离d为当前待处理对象外接圆半径的十分之一所对应的像素数。

(b)以(a)得到的待处理影像区域为数据域，以当前待处理对象为初始零水平集，表示为其中(x,y)为像素点坐标，t为时间参数。零水平集是一个闭合曲线对应的函数，是在x，y，t都确定情况下的函数取值，设在曲线内部曲线外部引入Heaviside函数和一维Dirac函数这两个函数为现有技术C-V变分水平集分割方法中所定义，为便于实施参考起见提供如下：

式(4)中为模型能量，Ω为待处理影像区域，u₁，u₂分别为零水平集内、外的平均灰度值，u为待处理影像区域某一点灰度值，为当前点梯度的模，α₁,α₂＞0,β≥0,γ≥0为固定参数，一般建议取α₁＝α₂＝1,β＝0,γ＝1；令u₁，u₂保持不变，使能量函数最小，得到其对应的偏微分方程：

根据式(5)演化提取水系对象，该演化过程严格按照现有C-V模型水平集演化技术执行即可，本发明不予赘述。

步骤3，基于形状曲线相似度匹配计算基础地理信息水系层矢量对象与步骤2分割对象，并根据多边形对象的匹配结果配准影像，如果配准成功则进入步骤4，如果配准不成功，则将当前分割结果作为初始水系分割结果对象集合代入步骤2，重新构建缓冲区并进行演化分割。

实施例进一步提供改进的形状曲线相似度匹配方式以提高匹配效果，实现包括以下步骤：

(a)设GIS数据中水系层已有的基础地理信息矢量数据集合为V＝{v₁，v₂，v₃，...v_m}，其中v_i(i＝1,2...m)为集合中的一个矢量多边形对象，m为集合中对象总数；设步骤2所得影像分割结果对象集合(即疑似水体分割对象集)为R＝{r₁，r₂，r₃，...r_n}，其中r_j(j＝1,2...n)，n为集合中对象总数。

取出V中每一对象分别计算与R中各对象的形状相似度，生成一个二维形状相似度矩阵S_vr[m，n]，其中m，n分别为V与R中对象总数。通过计算矢量对象的形状曲线之间的相似度来得到形状相似度。设V中某对象对应多边形A，R中某对象对应多边形B，形状曲线相似度计算公式如下：

$\mathrm{pCor}\_\mathrm{shape}(A,B)=\max \left(\frac{\underset{k=1}{\overset{n_p}{\mathrm{\Σ}}}({d_k}^A-\stackrel{\‾}{d^A})(d_{\mathrm{lk}}^B-\stackrel{\‾}{d_l^B})}{\sqrt{{\underset{k=1}{\overset{n_p}{\mathrm{\Σ}}}(d_k^A-\stackrel{\‾}{d^A})}^2\underset{k=1}{\overset{n_p}{\mathrm{\Σ}}}{(d_{\mathrm{lk}}^B-\stackrel{\‾}{d_l^B})}^2}}\right),l=\mathrm{1,2}...n_p---\left(7\right)$

式中为多边形A外边界的第k采样点到质心的距离，为所有的均值，表示多边形B由第l_k点到质心的距离，其中l_k＝l+k,l+k≤n_p;l_k＝l+k-n_p,l+k＞n_p，为所有的均值，n_p为多边形边界离散化后的点数，多边形A与B的边界点数相同，设其中某一个点为起始点，采用l和k标识相对于起始点的点对应的序号。

(b)将基础地理信息矢量数据集合V中对象分别构成的多边形按面积排序并从中依次选择面积最大的三个多边形{v_i1，v_i2，v_i3}；。

(c)从影像分割结果对象集合R中对象分别构成的多边形中，依次选择形状曲线相似度最大的多边形{r_j1，r_j2，r_j3}与{v_i1，v_i2，v_i3}匹配。

(d)计算多边形的质心坐标

{(x_vi1，y_vi1)，(x_rj1，y_rj1);(x_vi2，y_vi2)，(x_rj2，y_rj2);(x_vi3，y_vi3)，(x_rj3，y_rj3)}，式中(x_vi1，y_vi1)、(x_rj1，y_rj1)、(x_vi2，y_vi2)、(x_rj2，y_rj2)、(x_vi3，y_vi3)、(x_rj3，y_rj3)分别为多边形v_i1，r_j2，v_i2，r_j2，v_i3，r_j3对应的质心点在图像上的横纵坐标值，

代入式(8)，求解几何变换参数(a₀，a₁，a₂，b₀，b₁，b₂);，得到几何变换模型；

$\left\{\begin{array}{c}{x}_v=a_0+a_1\×x_r+a_2\×y_r\\ y_v=b_0+b_1\×x_r+b_2\×y_r\end{array}\right.---\left(8\right)$

其中，x_v,y_v,x_r,y_r为自变量，取值分别为

(x_vi1，y_vi1)，(x_rj1，y_rj1);(x_vi2，y_vi2)，(x_rj2，y_rj2);(x_vi3，y_vi3)，(x_rj3，y_rj3)。

(e)从V中剩余的m-3个对象中依次取一个对象，计算所取对象质心(x_vi，_yvi)与R中剩余n-3个对象质心基于式(8)变换后坐标的距离，设距离最小的点为(x_tri，_ytri)，如果两者的距离小于预设阈值D_tr(一般建议取为10个像素)，则将这两个质心作为一对匹配点，相应两个对象为匹配对象；遍历完V中所有对象，计算成功匹配对象数量占V中数据对象总数的比例，超过预设阈值P(一般默认为50%)，则认为当前几何变换模型是稳定的，根据当前几何变换模型以GIS数据为基准，对输入的光学遥感影像做几何变换(R也随之变换)，即配准成功；否则返回(c)，重新依次选择R中三个对象与V中已选出的对象匹配，选择策略为与V中选出对象次相似的R中对象；如果遍历了R中所有多边形均无法与V中选出的多边形对应，则返回(b)，从V中按面积依次选择其他三个多边形组合，再次迭代与R中的多边形匹配，直到获取到稳定的几何变换模型。如果遍历结合V和R中所有对象仍无法获得稳定的几何变换模型，配准不成功，则将R作为初始轮廓代替步骤1所得水系初始分割结果，返回步骤2重新在影像上演化分割水系对象。直到获得稳定的几何变换模型，配准成功后进入步骤4。

步骤4，对经步骤3配准后的V与R做基于缓冲区的变化检测，得到未变化水体R_w0、疑似变化水体R_w1、疑似新增水体R_nw。

实施例的具体操作步骤包括首先获取疑似变化水体R_w1、未变化水体R_w0、疑似新增水体R_nw，采用缓冲区变化检测方法如下：

(a)对V中对象以缓冲距离R_buffer构造缓冲区得到V_buffer(可称为GIS对象边缘缓冲区)，缓冲距离R_buffer可由本领域技术人员预设，例如设为10个像素；

(b)将R与V_buffer做叠加，基于式(9)分别以R中各对象为当前待检查的对象，统计在缓冲区外的长度占总长度的比例，；

$P_{\mathrm{Line}}=\frac{L_{\mathrm{out}\_\mathrm{buffer}}}{L_{\mathrm{total}}}---\left(9\right)$

其中，L_{out_buffer}为当前待检查的对象边缘位于缓冲区外部分的长度累加和，L_total为当前待检查的对象边缘总长度，P_Line为当前待检查的对象边缘位于GIS对象边缘缓冲区外部分长度占总长度的比例。

(c)根据预设阈值进行分类，例如实施例设为0.2，当P_Line＞0.2的时候，认为当前待检查的对象为疑似变化水体R_w1，等待进一步的检验；当P_Line≤0.2的时候，认为当前待检查的对象为未变化水体R_w0；R中与V在位置上不重叠的对象为疑似新增水体R_nw。

疑似变化水体R_w1中的对象可能由于季节水位变化导致影像上水体对象变化，并不是变化检测期望的检测结果，因此需要对R_w1中的对象做进一步的基于GIS骨架线相似性的比对，来获取实际变化水体对象和未变化水体对象。实施例的具体操作步骤如下：

(a)利用现有技术中的形态学细化算法提取R_w1中各对象的骨架线；

(b)栅格化GIS对象(指V中与R_w1在位置上重叠的对象)为二值图像，同样利用现有形态学细化算法提取对象骨架线；

(c)以GIS对象骨架线为基准构建缓冲区(缓冲距离可由本领域技术人员预设)，检测R_w1对象骨架线与GIS对象骨架线缓冲区的拓扑关系，与上述缓冲区变化检测方法类似，统计缓冲区外骨架线长度占总长度的比例，同样，当P_Line＞0.2的时候，认为当前对象为变化水体保留，反之为未变化水体，并入集合R_w0，得到R‘_w0。R_w1中保留的变化水体集合记为R‘_w1。基于骨架线的变化检测主要避免将河流在不同时期水位变化导致的变化做为变化，而只将河道走向发生变化作为真正的变化。

之前分割配准用的是基础地理信息数据中的水系层的数据，影像经过配准后与基础地理信息其他图层的数据也是套合的，可以利用道路和建筑层来过滤部分阴影的干扰，进一步提高最终的提取检测精度。因此，实施例还对疑似新增水体R_nw进行两重验证，即首先将R_nw与GIS数据中道路和建筑层对象做叠加分析，拓扑相交的对象为非水体对象则从R_nw中滤除；然后对R_nw中剩余对象进行上下文纹理特征过滤，提取R_nw中对象周边背景纹理特征，并与疑似水体纹理特征做比较，如果纹理相似度超过预设纹理特征相似度阈值，则认为该疑似水体对象为阴影的半影区域并从疑似水体对象中滤除，否则保留。过滤后剩余的疑似水体对象构成疑似水体对象集R‘_nw。实施例的上下文纹理特征过滤具体操作步骤如下：

(a)从R_nw的剩余对象中选择一个对象r_nwi，在其质心处构建一个矩形模板RECT_i，RECT_i的长度小于r_nwi最小外接矩形的宽度，实施例中RECT_i的长宽比为5，RECT_i的质心与r_nwi的质心重合；

(b)预设步长，实施例以15度为步长，绕质心旋转模板RECT_i，同时统计每次旋转后模板内像素的方差，可得到关于对象r_nwi的角度-方差图，记角度-方差图中方差最小值对应的方向为r_nwi的纹理主方向Direct_i，对应的矩形模板为RECT_imin；

(c)以r_nwi质心为原点，Direct_i和-Direct_i为方向构造射线Line_i1和Line_i2，射线与r_nwi的边缘有两个交点Pt_i1,Pt_i2，构造以Pt_i1,Pt_i2为端点，Direct_i,-Direct_i为方向，RECT_i的宽度为宽，长宽比为5的矩形对象Rect_i1,Rect_i2；

(d)分别统计Rect_i1,Rect_i2,RECT_imin对应区域像素的灰度共生矩阵生成的自相关统计量，记为COR_i1,COR_i2,COR_imin；

(e)如果区域灰度共生矩阵的自相关统计量差值|COR_i1-COR_imin|,|COR_i2-COR_imin|中有一个值小于预设纹理特征相似度阈值，则认为对象r_nwi与周边纹理相近，r_nwi为阴影区域，从R_nw中剔除该对象；如果上述两个差值计算结果均大于等于预设纹理特征相似度阈值，则保留该对象；本领域技术人员可根据具体精度要求自行设定阈值，实施例中预设纹理特征相似度阈值为0.3；

(f)返回(a)依次对R_nw中的所有剩余对象执行(a)-(e)，直到对R_nw中的所有剩余对象处理完毕。

在上述过程中涉及到的灰度共生矩阵和其自相关统计量为现有技术，为便于实施参考起见，提供具体实现方式：假设待定分析的纹理图像为一矩形图像，该图像某一区域像元的灰度可以量化为N层，则对应该区域的灰度共生矩阵是一个N×N阶的矩阵，N一般取值为256。在矩阵中位置(i_c,j_c)处的元素代表灰度为i_c的像素离开某个固定偏移量处灰度为j_c，i_c,j_c的取值为0,1...N-1。这种现象出现的概率。偏移量可以用δ＝(Δx,Δy)来表示，其中Δx为水平方向的偏移量，Δy为垂直方法偏移量。矩阵中的(i_c,j_c)位置的元素记为p(i_c,j_c)。自相关(Correlation)用来描述共生矩阵中行或列元素之间的相似程度，计算公式如下：

$\mathrm{COR}=[\underset{i_c=0}{\overset{N-1}{\mathrm{\Σ}}}\underset{j_c=0}{\overset{N-1}{\mathrm{\Σ}}}{i}_c{j}_{c}p(i_c,j_c)-{\mathrm{\μ}}_1{\mathrm{\μ}}_2]/{\mathrm{\σ}}_1{\mathrm{\σ}}_2---\left(10\right)$

其中，μ₁、μ₂分别为灰度的期望；σ₁、σ₂分别为灰度的方差，计算公式如下：

$\left\{\begin{array}{c}{\mathrm{\μ}}_1=\underset{i_c=0}{\overset{N-1}{\mathrm{\Σ}}}{i}_c\underset{j_c=0}{\overset{N-1}{\mathrm{\Σ}}}p(i_c,j_c),{\mathrm{\σ}}_1^2=\underset{i_c=0}{\overset{N-1}{\mathrm{\Σ}}}{(i_c-{\mathrm{\μ}}_1)}^2\underset{j_c=0}{\overset{N-1}{\mathrm{\Σ}}}p(i_c,j_c)\\ {\mathrm{\μ}}_2=\underset{j_c=0}{\overset{N-1}{\mathrm{\Σ}}}{j}_c\underset{i_c=0}{\overset{N-1}{\mathrm{\Σ}}}p(i_c,j_c),{\mathrm{\σ}}_2^2=\underset{j_c=0}{\overset{N-1}{\mathrm{\Σ}}}{(j_c-{\mathrm{\μ}}_2)}^2\underset{i_c=0}{\overset{N-1}{\mathrm{\Σ}}}p(i_c,j_c)\end{array}\right.---\left(11\right)$

步骤5，利用集合R’_w0对应像元的多光谱值作为正样本，随机选择配准后的光学遥感影像上集合R以外的区域中像元多光谱值作为负样本，训练SVM分类器，并采用训练所得SVM分类器验证步骤4过滤后剩余的疑似水体对象集R‘_nw中各对象是否为水体，在R‘_nw中保留验证结果为水体的对象，从R‘_nw中取出验证结果为非水体的对象。最终保留在R‘_nw中的对象即为分割配准提取的水体结果，即新增水体对象R“_nw。SVM分类器为现有技术，本发明不予赘述。最终，可输出配准后的影像以及新增水体对象R“_nw，变换水体对象R‘_w0，未变换水体对象R‘_w1。

以上所述的具体实施例，对本发明的技术方案进行了进一步详细说明，仅是用以说明本发明的具体实施案例而已，并非用以限定本发明的可实施范围。熟悉本领域的技术人员在不违背本发明所指示的精神与原理下所完成的一切等效变形、替换或修饰，仍包含在本发明权利要求所限定的范围内。

Claims (5)

1.一体化的光学遥感影像与GIS自动配准与水体提取方法，包括以下步骤：

步骤1，对输入的光学遥感影像进行分割，获得水系初始分割结果；

步骤2，以每个水系初始分割结果中的每个对象为当前待处理对象构建局部缓冲区，进行局部水平集演化分割，得到新的分割结果，构成的分割结果对象集合记为R；

步骤3，基于形状曲线相似度匹配GIS数据中水系层已有的基础地理信息水系层矢量对象的集合V与步骤2所得集合R中对象，以GIS数据为基准配准光学遥感影像，如果配准成功则进入步骤4，如果配准不成功则将当前的分割结果中每个对象作为当前待处理对象返回步骤2，重新构建局部缓冲区并进行局部水平集演化分割；步骤3中包括如下子步骤，

(3.1)设基础地理信息水系层矢量对象的集合为V＝{v₁，v₂，v₃，…v_m}，其中v_i为集合中的一个对象，i＝1,2…m，m为集合中对象总数；设步骤2所得集合R＝{r₁，r₂，r₃，…r_n}，其中r_j为集合中的一个对象，j＝1,2…n，n为集合中对象总数；

取出V中每一对象分别计算与R中各对象的形状相似度，生成一个二维形状相似度矩阵S_vr[m，n]，设V中某对象对应多边形A，R中某对象对应多边形B，按以下形状曲线相似度计算公式计算对象之间的形状相似度，

$pCor\_shape(A,B)=max\left(\frac{\underset{k=1}{\overset{n_p}{\Σ}}({d_k}^A-\stackrel{\‾}{d^A})(d_{lk}^B-\stackrel{\‾}{d_l^B})}{\sqrt{\underset{k=1}{\overset{n_p}{\Σ}}{(d_k^A-\stackrel{\‾}{d^A})}^2\underset{k=1}{\overset{n_p}{\Σ}}{(d_{lk}^B-\stackrel{\‾}{d_l^B})}^2}}\right),l=1,2...n_p---\left(7\right)$

其中，为多边形A外边界的第k采样点到质心的距离，为所有的均值，表示多边形B由第l_k点到质心的距离，其中l_k＝l+k,l+k≤n_p；l_k＝l+k-n_p,l+k＞n_p，为所有的均值，n_p为多边形边界离散化后的点数，多边形A与B的边界点数相同，设其中某一个点为起始点，采用l和k标识相对于起始点的点对应的序号；

(3.2)将集合V中各对象分别对应的多边形按面积排序并从中依次选择面积最大的三个多边形，记为{v_i1，v_i2，v_i3}；

(3.3)从集合R中各对象分别对应的多边形中，依次选择形状曲线相似度最大的多边形与{v_i1，v_i2，v_i3}匹配，记为{r_j1，r_j2，r_j3}；

(3.4)计算多边形v_i1，r_j1，v_i2，r_j2，v_i3，r_j3对应的质心点在图像上的横纵坐标值

{(x_vi1，y_vi1)，(x_rj1，y_rj1)；(x_vi2，y_vi2)，(x_rj2，y_rj2)；(x_vi3，y_vi3)，(x_rj3，y_rj3)}

采用下式求解几何变换参数(a₀，a₁，a₂，b₀，b₁，b₂)，得到几何变换模型；

$\left\{\begin{array}{c}{x}_v=a_0+a_1\×x_r+a_2\×y_r\\ y_v=b_0+b_1\×x_r+b_2\×y_r\end{array}\right.---\left(8\right)$

其中，x_v,y_v,x_r,y_r为自变量，取值分别为

(x_vi1，y_vi1)，(x_rj1，y_rj1)；(x_vi2，y_vi2)，(x_rj2，y_rj2)；(x_vi3，y_vi3)，(x_rj3，y_rj3)；

(3.5)从集合V中剩余的m-3个对象中依次取一个对象，计算所取对象质心(x_vi，y_vi)与集合R中剩余n-3个对象质心基于式(8)变换后坐标的距离，设距离最小的质心为(x_tri，y_tri)，如果质心(x_vi，y_vi)与(x_tri，y_tri)的距离小于预设阈值D_thr，则将这两个质心作为一对匹配点，相应的两个对象为匹配对象；遍历完集合V中所有对象，计算成功匹配对象数量占集合V中对象总数的比例，超过预设阈值P则认为当前几何变换模型是稳定的；否则返回(3.3)，重新根据形状曲线相似度依次选择集合R中三个对象与集合V中已选出的对象匹配；如果遍历了集合R中所有多边形均无法与集合V中选出的多边形对应，则返回(3.4)，从集合V中按面积依次选择其他三个多边形，再次迭代与R中的多边形匹配，直到获取到稳定的几何变换模型；获取到稳定的几何变换模型后以GIS数据为基准配准光学遥感影像，配准成功；如果遍历结合V和R中所有对象仍无法获得稳定的几何变换模型，则配准不成功；

步骤4，对经步骤3配准后的V与R做基于缓冲区的变化检测，得到未变化水体R_w0、疑似变化水体R_w1和疑似新增水体R_nw；然后分别对疑似变化水体R_w1和疑似新增水体R_nw做进一步的验证如下，

对疑似变化水体R_w1做基于骨架线相似性的检测，确定真正变化的水体和未变化的水体，未变化的水体并入未变化水体R_w0，得到集合R‘_w0，R_w1中保留的变化水体记为集合R‘_w1；

对疑似新增水体R_nw进行两重验证，包括首先将R_nw与GIS数据中道路和建筑层对象做叠加分析，拓扑相交的对象则为非水体对象并从R_nw中滤除；然后对R_nw中剩余对象进行上下文纹理特征过滤，包括提取R_nw中对象周边背景纹理特征，并与疑似水体纹理特征做比较，如果纹理相似度超过预设纹理特征相似度阈值，则认为该疑似水体对象为阴影的半影区域并从R_nw中滤除，过滤后剩余的疑似水体对象构成集合R‘_nw；

步骤5，利用集合R‘_w0中对应像元的多光谱值作为正样本，随机选择光学遥感影像上集合R以外区域中的像元多光谱值作为负样本，训练SVM分类器，并采用训练所得SVM分类器验证步骤4过滤后剩余的集合R‘_nw中各对象是否为水体，在R‘_nw中保留验证结果为水体的对象，从R‘_nw中取出验证结果为非水体的对象，最终保留在R‘_nw中的对象记为集合R“_nw，以集合R“_nw为分割配准提取的水体结果。

2.根据权利要求1所述一体化的光学遥感影像与GIS自动配准与水体提取方法，其特征在于：步骤1中，判断输入的光学遥感影像是否具有NIR波段，有则采用NDWI分割，否则采用OTSU分割。

3.根据权利要求1或2所述一体化的光学遥感影像与GIS自动配准与水体提取方法，其特征在于：步骤2实现方式为，设步骤1所得水系初始分割结果中的所有对象构成一个集合，以集合中的任一对象为当前待处理对象，进行以下步骤，

(2.1)对当前待处理对象，按预设缓冲距离d构建缓冲区，得到待处理影像区域；

(2.2)以(2.1)得到的待处理影像区域为数据域，以当前待处理对象为初始零水平集，表示为其中(x,y)为像素点坐标，t为时间参数，设在曲线内部曲线外部引入Heaviside函数和一维Dirac函数建立模型的能量函数如下，

其中，Ω为待处理影像区域，u₁，u₂分别为零水平集内、外的平均灰度值，u为待处理影像区域某一点灰度值，为当前点梯度的模，固定参数α₁,α₂＞0,β≥0,γ≥0；令u₁，u₂保持不变，使能量函数最小，得到其对应的偏微分方程：

根据式(5)，采用C-V模型水平集演化方法分割提取水系的对象。

4.根据权利要求1或2所述一体化的光学遥感影像与GIS自动配准与水体提取方法，其特征在于：步骤4中对经步骤3配准后的V与R做基于缓冲区的变化检测，得到未变化水体R_w0、疑似变化水体R_w1、疑似新增水体R_nw，包括如下子步骤，

(4.11)对V中对象以预设缓冲距离R_buffer构造缓冲区V_buffer；

(4.12)将R与V_buffer做叠加，基于式(9)分别以R中各对象为当前待检查的对象，统计在缓冲区外的长度占总长度的比例P_Line，

$P_{Line}=\frac{L_{out\_buffer}}{L_{total}}---\left(9\right)$

其中，L_{out_buffer}为当前待检查的对象边缘位于缓冲区外部分的长度累加和，L_total为当前待检查的对象边缘总长度，P_Line为当前待检查的对象边缘位于GIS对象边缘缓冲区外部分长度占总长度的比例。

(4.13)当P_Line大于预设阈值时，认为当前待检查的对象为疑似变化水体R_w1，当P_Line小于等于预设阈值时的时候，认为当前待检查的对象为未变化水体R_w0；R中与V在位置上不重叠的对象为疑似新增水体R_nw。

5.根据权利要求4所述一体化的光学遥感影像与GIS自动配准与水体提取方法，其特征在于：步骤4中对R_nw中剩余对象进行上下文纹理特征过滤的实现方式如下，

(4.21)从R_nw的剩余对象中选择一个对象r_nwi，在其质心处构建一个矩形模板RECT_i，RECT_i的长度小于r_nwi最小外接矩形的宽度；

(4.22)根据预设步长绕质心旋转模板RECT_i，同时统计每次旋转后模板内像素的方差，得到关于对象r_nwi的角度-方差图，记角度-方差图中方差最小值对应的方向为r_nwi的纹理主方向Direct_i，对应的矩形模板为RECT_imin；

(4.23)以r_nwi质心为原点，Direct_i和-Direct_i为方向构造射线Line_i1和Line_i2，射线与r_nwi的边缘有两个交点Pt_i1,Pt_i2，构造以Pt_i1,Pt_i2为端点，Direct_i,-Direct_i为方向，RECT_i的宽度为宽，长宽比为5的矩形对象Rect_i1,Rect_i2；

(4.24)分别统计Rect_i1,Rect_i2,RECT_imin对应区域像素的灰度共生矩阵生成的自相关统计量，记为COR_i1,COR_i2,COR_imin；

(4.25)如果自相关统计量差值|COR_i1-COR_imin|,|COR_i2-COR_imin|中有一个值小于预设纹理特征相似度阈值，则认为对象r_nwi与周边纹理相近，r_nwi为阴影区域，从R_nw中剔除该对象；如果|COR_i1-COR_imin|,|COR_i2-COR_imin|上述两个差值计算结果均大于等于预设纹理特征相似度阈值，则保留该对象；

(4.26)返回(4.21)依次对R_nw中的所有剩余对象执行(4.22)-(4.25)，直到对R_nw中的所有剩余对象处理完毕。