CN101126812A - 高分辨率遥感影像的分割和分类及变化检测一体化方法 - Google Patents

Abstract

本发明公开了一种基于多水平集演化的高分辨率遥感影像分割和分类及变化检测一体化方法。其特征是：(1)影像预处理(辐射、配准及特征滤波)；(2)建立多水平集演化的分割及分类模型，用配准后GIS数据自动确定各水平集函数的初始轮廓，对第1时相影像分割和分类；(3)仍采用(2)中所述模型，优化各水平集函数的初始轮廓，对第2至第T时相影像均采用增量式分割及分类；(4)以分割后对象为单元，比较第i与第i+1两个相邻时相影像分类结果，确定变化区域；(5)返回至(3)，直到完成所有T个时相影像的分割和分类及变化检测。优点：与传统面向像元的K均值法相比，分类和检测精度均有所提高，适于序列遥感影像的变化检测，在灾害监测和国土资源调查等方面有广阔应用。

Description

高分辨率遥感影像的分割和分类及变化检测一体化方法

技术领域

本发明属于计算机遥感影像信息提取技术领域，涉及一种基于多水平集曲面演化的高空间分辨率遥感影像分割和分类及变化检测一体化方法。

背景技术

遥感影像信息提取是计算机与遥感等领域的一个基本而又重要的研究内容。但随着传感器技术，卫星通讯技术、计算机技术的迅猛发展，遥感影像愈来愈呈现出“三多”的特点，即多传感器、多种空间或光谱分辨率、多时相，这使得遥感影像自动信息提取更加困难。

近年来，高空间分辨率(空间分辨率小于或等于5米，以下简称高分辨率)卫星，如IKONOS，QuickBird，以及国内的军事卫星、小卫星等均具有很高的空间分辨率，它为使用者提供了丰富的地物信息，但也给自动信息提取任务带来了更大的困难。随着高分辨率影像在土地资源调查、城市规划、灾害监测等方面的应用日益广泛，传统的面向像元的分割、分类及变化检测方法已不再适用，其表现的不足主要在于：(1)影像中噪声及过多的细节常常会导致分割、分类及变化检测结果中出现“椒盐”现象，其输出结果无法代表一个有意义的目标，与实际地物不相吻合；(2)它未充分考虑高分辨率影像除光谱之外的纹理，结构、几何及地物间的空间关系等特征，分类精度低；(3)其分割和分类结果不便于直接矢量化，不便于与地理信息数据建立联系。前人工作已证明了面向对象的信息提取技术是缓解这些不足的有效手段，相关文献有：G..J.Hay，G.Castilla，“Object-based image analysis：strengths，Weakness，Opportunities and Threats(SWOT)，”OBIA2006：International Archive of Photogrammetry，Remote Sensing and Spatial Information Sciences，July，2006.http://www.commission4.isprs.org/obia06/papers.htm。商业软件包eCognition中的多尺度分割及在分割基础上的分类正是这类技术的代表，它改善了传统分割及分类结果中的“椒盐”现象，提高了分类精度。此外，分水岭(watershed)、均值漂移(meanshift)等多尺度分割基础上的分类也是处理高分辨率遥感影像的有效途径。但这类方法中，分割方法的性能对分类精度有很大影响。

变化检测是遥感影像信息提取技术领域的一个重要内容，它涉及影像辐射校正、配准、分割、分类识别等多个环节，误差传播严重，现在的变化检测方法大致分两类：

●影像代数运算方法。它对辐射校正及配准精度要求很高，不适用于多时相异质传感器影像的变化检测问题。

●分类后比较法。该方法中，分类是关键，根据分类技术不同又将其分为面向像元和面向对象的变化检测方法，前者分类精度低，分类结果不便于GIS数据库更新，且对配准误差要求非常高，一般小于0.5个像元；后者分类精度高，对配准误差要求相对较低，误差值视影像空间分辨率及变化检测目标而定。此外，这类方法大多针对有限多个时相数(2-10个时相)的遥感影像。

对于高分辨遥感影像，采用面向对象的思想进行影像分割或分类或变化检测是一种行之有效的途径，参见文献：G.Willhauck，Comparison of object oriented classificationtechniques and standard image analysis for the use of change detection between SPOTmultispectral satellite images and aerial photos.ISPRS，Vol.XXXIII，Amsterdam，2000，以及Niemeyer，I.&M.J.Canty，Pixel-Based and Object-Oriented Change DetectionAnalysis Using High-Resolution Imagery.Proc.25th Symposium on Safeguards andNuclear Material Managment，Stockholm，Sweden，13-15 May 2003。基于以上思想，本专利公开了一种新的面向对象的高分辨率遥感影像信息提取方法，它是一种基于多水平集曲面演化的分割和分类及变化检测一体化方法。

水平集理论(level set theory)是处理封闭运动界面随时间演化过程中几何拓扑变化的有效计算工具。其基本原理是：通过水平集函数Φ(x，y，t)表示的三维曲面在图像内、外力的作用下，改变其拓扑形状，当曲面演化速度趋于零时，该水平集函数的零水平集隐含表示的平面曲线即可表示目标或区域边界，相关文献有：S.Osher，J.A Sethian，“Fronts Propagatingwith Curvature Dependent Speed：Algorithms Based on Hamilton-Jacobi Formulations，”Journalof Computational Physics，vol.79，pp.12-49，1988和S.Osher，R.Fedkiw.Level Set Methods andDynamic Implicit Surfaces.New York：Springer-Verlag New York Inc.，ch.3，pp.25-37，2003。它将二维的图像分割问题转化成三维的曲面演化问题，通过求解与时间t有关的偏微分方程(组)的来实现图像分割。它是一种具有拓扑自由变化的主动轮廓方法，可实现全局能量最优准则下的图像划分。其优点在于：不仅考虑图像的区域统计特征或目标梯度特征，而且利用了曲线的几何特性(曲率、法向力)，并可融合目标的先验知识等，最后将这些信息作为曲线演化动力集成在水平集函数的演化方程中。

自2000年以来，水平集理论逐渐成为国、内外流体力学、材料学、数学、计算机视觉、图像处理等领域学者研究的热点。其中，在医学图像中的研究及应用颇多，相关文献有：L.Vese，T.F.Chan，“A Multiphase Level Set Framework for Image Segmentation Using theMumford and Shah Model，”International Journal of Computer Vision，vol.50，no.3，pp.271-293，July，2002，及文献H.K Zhao，T.F.Chan，B.Merriman，S.Osher，“Variational Level SetApproach to Multiphase Motion，”Journal of Computational Physics，no.0167，pp.179-195，February，1996。这证明了水平集理论在图像分割和分类中的可行性，但以上方法中存在其不足，如分类数受水平集个数制约或存在模糊分割现象等。

不同于医学图像，遥感影像具有成像平台的多样性、成像环境的复杂性以及多波段、多时相、高空间分辨率等特性，因此，提出行之有效的信息提取方法一直是人们关注的热点。多水平集理论在遥感影像分割或分类等方面研究不多。相关文献有：A.R.Mansouri，A.Mitiche，C.Va′zquez，“Multiregion competition：A level set extension of regioncompetition to multiple region image partitioning，”Computer Vision and ImageUnderstanding，no.101，pp.137-150，2006中，作者主要是针对医学、视频图像提出的多水平集演化的分割方法，但未涉及高空间分辨率或多光谱影像，且未充分考虑遥感影像的成像环境复杂，地物类型繁多等特点。文献：O.Besbesy，Z.Belhadj，N.Boujema，“AdaptiveSatellite Images Segmentation by Level Set Multiregion Competition，”INRIA，Technical Report，no.5855，March，2006基于A.R.Mansouri的工作提出了一个自适应的分类模型；文献：Ball，J.E.Bruce，L.M.Level set segmentation of remotely sensed hyperspectral images.IGARSS 2005.vol.8，page(s)：5638-5642中，作者针对高光谱影像，建立了多水平集的分类模型，但它未充分利用遥感影像的纹理，几何等更多的地物特征，准确性尚有待改善。以上相关文献中的方法均是监督性的，在分割或分类前需要样本学习和参数估计，增加了工作量；

在变化检测的应用方面，佛罗里达州立大学的Vivien Mallet等人基于单水平集演化的思想开发了Multivac软件包，用于对火势蔓延过程进行建模以及模拟硅纳米材料的生长，相关文献见http://vivienmallet.net/fronts/index.php。此外，文献Dell′Acqua，F.Gamba，P.Prevedini，P.，Level-set based extraction and tracking of meteorological objects in satellite images，IGARSS 2000.vol.2，page(s)：627-629中，利用水平集理论的动态演化特性设计了中低分辨率气象卫星数据的云目标跟踪系统，以上工作证实了水平集理论在众多时相遥感影像中单目标变化检测方面的可应用性，但针对土地利用状况调查等任务中涉及的多类目标的变化监测，以上方法无效，需采用多水平集(水平集函数个数大于1)分割或分类方法。

在图像处理中，分割和分类的概念相似。在计算机视觉等领域，一般称图像分割。由于现有的水平集分割研究者主要来自于应用数学、计算机视觉及医学图像处理领域，因此在该领域的相关文献中常常使用图像分割这个术语，不涉及分类的概念；而在模式识别及遥感领域，分割和分类是两个不同的过程，概念上也有差别，前者是指将影像划分为多个同质区域，不关心同质区的个数及每个同质区的地物属性，而后者是将影像划分成不同类别的地物，一般用户很关心所分的类别数及每一类的地物属性。水平集演化理论在遥感影像的分割及分类应用研究中，部分作者将分类过程用分割术语来表示，如作者A.R.Mansouri et.al.以及作者O.Besbesy et.al.发表的文献，实质上，他们仅实现了影像分类的过程。本发明是面向遥感领域，它区分了这两个概念，所提出的方法不仅实现了分类，而且给出了图像分割结果，但分割和分类结果的表达方法不同，前者是用各水平集函数的符号将整个影像划分成不同类别后，再逐个标定各类别所包含的所有闭合区域，形成同质区；而后者仅仅是用各水平集函数的符号将整个影像划分成不同类别。

发明内容

本发明所要解决的问题是：提供一种基于多水平集曲面演化的高分辨率遥感影像分割和分类及变化检测的一体化方法，该方法以区域为对象通过平面闭合曲线簇的演化来实现影像的分割和分类，是一种新的面向对象的分割和分类及变化检测思想，有效缓解了传统的面向像元的信息提取方法中的“椒盐”现象。与传统K均值分类后的变化检测方法相比，其分类及变化检测精度较高。

本发明提供的技术方案是：一种基于多水平集曲面演化的高分辨率遥感影像分割和分类及变化检测的一体化方法，包括以下步骤：

一、高分辨率遥感影像的预处理：

a)对各时相高分辨率遥感影像分别进行辐射校正；并对辐射校正后的各时相高分辨率遥感影像进行配准，然后再与同一地区的GIS数据配准；

b)对辐射校正和配准之后的高分辨率遥感影像提取纹理特征及光谱特征，组成m维特征矢量影像

，其中，m是提取的特征矢量影像的维数，且m≥1；u_k表示第k个分量的特征影像，1≤k≤m；然后利用矢量总变分的多尺度异向扩散技术进行滤波；

二、对预处理之后的特征矢量影像建立分割和分类模型：

首先，建立多水平集演化方程组，即偏微分方程组：

$\left\{\begin{array}{c}\frac{\∂{\mathrm{\Φ}}_1}{\∂t}(x,y)=[-{\mathrm{\ξ}}_1(x,y)+{\mathrm{\φ}}_1(x,y)+\mathrm{\μ}{\mathrm{Cur}}_{{\mathrm{\Φ}}_1}]\left|\right|\▿{\mathrm{\Φ}}_1(x,y)\left|\right|\\ \·\\ \·\\ \·\\ \frac{\∂{\mathrm{\Φ}}_i}{\∂t}(x,y)=[-{\mathrm{\ξ}}_i(x,y)+{\mathrm{\φ}}_i(x,y)+\mathrm{\μCu}{r}_{{\mathrm{\Φ}}_i}]\left|\right|\▿{\mathrm{\Φ}}_i(x,y)\left|\right|\\ \·\\ \·\\ \·\\ \frac{{\∂\mathrm{\Φ}}_{N-1}}{\∂t}(x,y)=[-{\mathrm{\ξ}}_{N-1}(x,y)+{\mathrm{\φ}}_{N-1}(x,y)+\mathrm{\μ}{\mathrm{Cur}}_{{\mathrm{\Φ}}_{N-1}1}]\left|\right|\▿{\mathrm{\Φ}}_{N-1}(x,y)\left|\right|\end{array}\right.---\left(1\right)$

式(1)中，

是隐含表示N-1维水平集函数矢量，N是高分辨率遥感影像包含的地物类别个数；函数ξ_i(x，y)是二维影像域Ω中第i个类别C_i的统计特征，计算公式为：

${\mathrm{\ξ}}_i(x,y)=\frac{1}{m}\underset{k=1}{\overset{m}{\mathrm{\Σ}}}{(u_k(x,y)-c_{\mathrm{ik}})}^4,$ 1≤i≤N    (2)

式(2)中，c_ik是二维影像域Ω中第k维特征矢量影像上第i个类别C_i的均值，计算公式为：

$c_{\mathrm{ik}}=\frac{{\∫}_{C_i}{u}_k(x,y)\mathrm{dxdy}}{{\∫}_{C_i}\mathrm{dxdy}}---\left(3\right)$

式(3)中，积分符号∫_Ci表示积分域是第i个类别C_i包含的像元集合；式(1)中，函数φ_i(x，y)的计算公式为：

${\mathrm{\φ}}_i(x,y)={\mathrm{\ξ}}_{i+1}(x,y){\mathrm{\χ}}_{\{{\mathrm{\Φ}}_{i+1}(x,y,t)>0\}}(x,y)+{\mathrm{\ξ}}_{i+2}(x,y){\mathrm{\χ}}_{\{{\mathrm{\Φ}}_{i+1}(x,y,t)\≤0\}}(x,y){\mathrm{\χ}}_{\{{\mathrm{\Φ}}_{i+2}(x,y,t)>0\}}(x,y)+\·\·\·$

$+{\mathrm{\ξ}}_{N-1}(x,y){\mathrm{\χ}}_{\{{\mathrm{\Φ}}_{i+1}(x,y,t)\≤0\}}(x,y)\·\·\·{\mathrm{\χ}}_{\{{\mathrm{\Φ}}_{N-2}(x,y,t)\≤0\}}(x,y){\mathrm{\χ}}_{\{{\mathrm{\Φ}}_{N-1}(x,y,t)>0\}}(x,y)$

$+{\mathrm{\ξ}}_N(x,y){\mathrm{\χ}}_{\{{\mathrm{\Φ}}_{i+1}(x,y,t)\≤0\}}(x,y)\·\·\·{\mathrm{\χ}}_{\{{\mathrm{\Φ}}_{N-2}(x,y,t)\≤0\}}(x,y){\mathrm{\χ}}_{\{{\mathrm{\Φ}}_{N-1}(x,y,t)\≤0\}}(x,y)---\left(4\right)$

式(1)中，Cur_Φi是第i个水平集函数的平均曲率，计算公式为： ${\mathrm{Cur}}_{{\mathrm{\Φ}}_i}=\▿\·\frac{{\▿\mathrm{\Φ}}_i}{\left|\right|{\▿\mathrm{\Φ}}_i\left|\right|},$

为梯度算子；参数μ是平均曲率的权值；X_{Φi(x，y，t)}(x，y)是第i个水平集函数的示性函数，定义为：

给定约束条件：

$C_i(x,y,t)=\underset{1\≤i\≤N-1}{\min }\left\{{\mathrm{\ξ}}_i(x,y,t)\right\}---\left(5\right)$

式(5)中，C_i(x，y，t)表示第t次迭代中，ξ_i(x，y，t)值最小的第i个水平集函数对应的像元(x，y)的类标号；当多水平集演化过程满足收敛条件时，类标号C_i的变化处于静止状态；

带有约束条件(5)的偏微分方程组(1)就是高分辨率遥感影像的分割和分类模型，该模型的解(即影像分割和分类结果)采用以下方式来表示：第1至第N-1类表示为：C_i(x，y，t)＝{(x，y)∈Ω|Φ_i(x，y，t)＞0}，其中，1≤i≤N-1，第N类表示为： $C_N(x,y,t)=\{(x,y)\∈\mathrm{\Ω}|\underset{i=1}{\overset{N-1}{\∩}}{\mathrm{\Φ}}_i(x,y,t)\≤0\},$ 而分割结果就是所有水平集函数的零水平

集围成的区域；

三、第1个时相高分辨率遥感影像的分割和分类：

a)利用有限差分法对分割和分类模型进行离散化，边界条件是Neumann边值条件；输入式(1)-(3)中参数：影像包含的地物类别个数N取正整数，且N≥2，这样，分割和分类模型的所需的水平集函数个数为N-1；平均曲率的权值μ＝(0.1)^p×255²，p为整数，取值范围为0～5；行、列方向上的空间采样间隔Δx＝1，Δy＝1，即一个像元；迭代时间步长Δt依赖于空间采样间隔及水平集函数的演化速度，它的取值必须满足CFL条件；

b)在对第1时相的高分辨率遥感影像水平集分割和分类中，分割和分类模型的初始条件设置方法如下：在配准后的该地区GIS数据辅助下，设置每一个水平集函数的初始轮廓，具体方法：利用GIS数据中每一个地物类别的空间位置及属性来画一个初始圆，并将每个初始圆进行符号距离变换作为各个水平集函数的初始条件，采用步骤三中a)所述的边界条件和离散化方式及参数值来实现第1时相影像的水平集分割和分类，分割和分类结束时得到分割后的平面闭合曲线簇；

四、第2至第T时相的高分辨率遥感影像采用增量式分割和分类：对第2至第T时相的高分辨率遥感影像仍采用步骤二中所述的分割和分类模型，但采用新的初始条件及新的影像作用域，新的初始条件设置方法如下：将第j时相高分辨率遥感影像分割后的平面闭合曲线簇作为与其相邻的第j+1时相高分辨率遥感影像分割和分类前的初始轮廓；仍采用步骤三中a)所述的边界条件和离散化方式及参数值来实现第j+1时相的高分辨率遥感影像的分割和分类，分割和分类结束时得到分割后的平面闭合曲线簇；其中：T为变化检测影像的时相总数；j为整数变量，表示第j时相的高分辨率遥感影像，取值范围为是1≤j≤T-1；

五、第j与第j+1两个相邻时相高分辨率遥感影像的变化检测：从j＝1开始，当第j与第j+1时相高分辨率遥感影像的分割和分类过程结束后，开始对这两个相邻时相高分辨率遥感影像进行变化检测，具体方法是：对比这两个相邻时相高分辨率遥感影像的分类结果，以其中任意一个时相高分辨率遥感影像的分类结果为参考，在另一时相高分辨率遥感影像的分类结果中，以分割后的对象为比较单元，若任意一个对象的全部像元或部分像元的类标号发生了改变，则表明该对象的地物类型发生了变化，并输出所有变化像元的空间位置及类别属性；

六、返回第四步，循环执行步骤四至六的过程，直到完成了所有的T个时相影像的分割和分类及变化检测任务。

对大小为512×512以上的T个不同时相的高分辨率遥感影像分割和分类及变化检测可采用以下步骤：按步骤一进行预处理；再采用小波多分辨率分解；按照上述步骤三至六进行分割和分类及变化检测，直到完成了所有的T个不同时相影像的分割和分类及变化检测任务。

在没有同一地区的GIS矢量图的情况下，上述步骤二中所述的分割和分类模型的初始条件设置方法为：在影像作用域范围内自动均匀地画任意半径的初始圆，再进行符号距离变换来初始化各个水平集函数。

本发明建立了一种基于区域的多水平集演化模型来实现面向对象的影像分割和分类及变化检测方法，它的优点表现在：

a)它利用水平集框架的多源信息可融合性，建立了融合光谱、纹理、几何等多种特征的多波段的高分辨率遥感影像分割和分类模型；它是对A.R.Mansouri提出方法的简化，参数少且无需样本学习及参数估计过程，优化了分割和分类速度，对大小为512×512的四波段的QuickBird影像，其分割和分类过程耗费时间平均为77秒；

b)实现了面向对象的分割及分类，避免了面向像元方法的“椒盐”现象，缓解了传统的分水岭，区域增长方法中表现的过分割或欠分割现象，分类精度提高了2％-10％，相比传统的K均值分类后的变化检测方法，其检测精度提高1％-8％，且信息提取结果易用于GIS数据库的更新；

c)它利用水平集曲面在演化过程中表现出的分裂、合并等自由拓扑改变特性，同时实现了高分辨率遥感影像的分割与分类，改善了现有方法中分割与分类独立进行的模式，减少了误差传播途径，不需要分类后处理；

d)它凭借水平集曲面演化的连续性及动态性，在变化检测的分类环节中采用了增量式分割和分类的思想；它通过优化模型的初始化条件来加速单个时相的高分辨率遥感影像的分割和分类过程，与未优化初始条件的变化检测过程相比，这种增量式分割和分类技术使整个变化检测过程所耗费的时间缩短了5％-15％；

e)本发明是面向对象的分割及分类后变化检测方法，对配准误差具有一定的鲁棒性。

本发明在军事目标检测中，能有效降低虚检及误检率；在灾情监测、土地利用变化检测方面，该方法灵活，模型参数较少；相比传统的K均值分类后的变化检测方法，对序列遥感影像的检测精度提高1％-8％，且信息提取结果易用于地理信息数据库更新。

与现有的K均值以及eCognition软件包中的多分辨率分割方法相比，该方法对512×512以上的大幅遥感影像多类地物的信息提取中仍存在速度方面的不足，但在方程的数值求解方面还有优化的潜力，因此在遥感信息提取技术领域有很大的应用前景。

附图说明

图1为本发明的主程序运行流程图。

图2为本发明的分割和分类结束时区域及类别标定示意图。

图3为本发明的增量式分割和分类示意图。

图4为本发明的分类后面向对象的变化检测示意图。

具体实施方式

参见图1、图2、图3和图4，本发明提供一种基于多水平集演化的分割和分类及变化检测一体化方法。设影像的时相数为T(T为大于1的正整数)，则对第1至第T时相的分割和分类及变化检测实施流程参见附图1，包括如下步骤：

1、原始影像的预处理

a)通常，绝对辐射校正所需参数难以获取，因此，对同质传感器的不同时相遥感影像采用常见的相对辐射校正，如直方图匹配法、暗集一亮集法等，这两种方法及更多方法参见文献：丁丽霞，周斌，王人潮；遥感监测中5种相对辐射校正方法研究．浙江大学学报(农业与生命科学版).2005，31(3)：269～276；对异质传感器的不同时相遥感影像，在不具备绝对辐射校正参数的情况下，略去辐射校正操作；

b)对辐射校正后的高分辨率遥感影像运用多项式法或有理函数法等进行影像之间的互配准，具体操作方法参见文献：孙家柄，舒宁，关泽群，遥感原理、方法和应用，北京：测绘出版社，1997；以及文献：张永生、巩丹超等著，高分辨率遥感卫星应用--成像模型、处理算法及应用技术，北京：测绘出版社，2005；

c)对辐射校正和配准之后的不同时相影像，采用Gabor滤波器或非线性结构张量法提取纹理特征[具体操作方法参见文献：赵银娣，马尔可夫模型的扩展及其在高分辨率影像分类中的应用，博士学位论文，武汉大学，2006，以及Thomas Brox，Joachim Weickerta.A TV flow based local scale estimate and its application to texture discrimination，Journalof Visual Communication and Image Representation，2006，17(5)：1053-1073]，再综合原始影像的所有或特征选择后的光谱波段(特征选择具体操作方法参见文献：孙家柄，舒宁，关泽群，遥感原理、方法和应用，北京：测绘出版社，1997)，一同构成m维特征矢量空间，其中，m是提取的特征矢量的维数，且m≥1；再利用矢量总变分流进行特征空间滤波，改善区域内部的异质性，增强区域间的特征差异，并保持重要的区域边界，具体操作方法参见文献：P.Blomgren,T.F.Chan，“Color TV：Total VariationMethods for Restoration of Vector-Valued Images,”IEEE Transactions on ImageProcessing，1998，7(3)：304-309。

2、建立多水平集演化的影像分割和分类模型

a)建立关于平面闭合曲线簇的能量方程

假设位于二维影像作用域Ω上的一个以弧长s、时间t为参数的闭合曲线簇

1＜i＜N-1，建立曲线簇 ${\stackrel{\→}{\mathrm{\γ}}}_i(x(s,t),y(s,t)):\left[\mathrm{0,1}\right]\→\mathrm{\Ω},$ 设区域

分别表示曲线

的内、外区域，曲线簇的能量方程可表示为式(1)，该方程的具体涵义及建立方法参见文献：A.R.Mansouri，A.Mitiche，C.Va′zquez，“Multiregioncompetition：A level set extension of region competition to multiple region imagepartitioning,”Computer Vision and Image Understanding，no.101，pp.137-150，October，2006。

$E\left[{\left\{{\stackrel{\→}{\mathrm{\γ}}}_i\right\}}_{i=1}^{N-1}\right]={\∫}_{R_{{\stackrel{\→}{\stackrel{\→}{\mathrm{\γ}}}}_1}}{\mathrm{\ξ}}_1(x,y)\mathrm{dxdy}+{\∫}_{R_{\stackrel{\→}{\mathrm{\γ}}2}}{\mathrm{\ξ}}_2(x,y)\mathrm{dxdy}+...$

$+{\∫}_{R_{\stackrel{\→}{\mathrm{\γ}}1}^c\∩R_{\stackrel{\→}{\mathrm{\γ}}2}^c\∩...\∩R_{\stackrel{\→}{\mathrm{\γ}}k-1}^c\∩R_{\stackrel{\→}{\mathrm{\γ}}k}}{{\mathrm{\ξ}}_k}_{}(x,y)\mathrm{dxdy}+...$

$+{\∫}_{R_{\stackrel{\→}{\mathrm{\γ}}1}^c\∩R_{\stackrel{\→}{\mathrm{\γ}}2}^c\∩...\∩R_{\stackrel{\→}{\mathrm{\γ}}N-2}^c\∩R_{\stackrel{\→}{\mathrm{\γ}}N-1}^c}{\mathrm{\ξ}}_N(x,y)\mathrm{dxdy}$

(1)

式(1)中，N是预先定义的影像包括的地物类别个数；最后一项是曲线簇总长度的惩罚项，μ是惩罚因子，μ是常数且μ＞0，这个惩罚项保证了分割边界的光滑性，这个惩罚因子也可视为分割尺度参数，此值越大，分割对象面积越大。其值μ＝(0.1)^p×255²，p为整数，取值范围为0～5；

能量方程(1)利用了N-1维参数化的曲线簇

将影像域Ω划分为N个类别，即， $\mathrm{\Ω}=\underset{i=1}{\overset{N}{\∪}}{C}_i,$ C_i表示第i个类别，每个类别都是一个像元集合。由于高分辨率影像分割和分类过程是一个不断优化求解的过程，因此，C_i包含的像元个数及空间位置均随时间不断变化，当能量达到全局最小或局部最小时，类标号C_i的变化处于静止状态，此时分割和分类过程结束。尽管该能量方程(1)形式与作者A.R.Mansouri所述方法中的能量方程相同，但对其具体的计算表达式已进行了改进。改进之处主要在于：

i)ξ_i(x，y)函数的涵义及计算方面：在本方法中，ξ_i(x，y)函数反映了区域

的统计特征，它是整个曲线簇演化的主要驱动力，它的计算公式如式(2)示所示。改进之后的能量公式(1)主要针对多波段高分辨率遥感影像的分割和分类，如IKONOS全色与多光谱融合后的影像；它不仅考虑了高分辨率遥感影像的光谱特征，还充分利用了其重要的空间纹理特征，使分割和分类结果更准确；此外，ξ_i(x，y)的计算无需类概率分布假设，从而省去了繁琐的样本学习及参数估计过程，使能量方程(1)参数更少，计算更简便；重要的是，ξ_i(x，y)函数的这一计算方法避免了协方差矩阵及其求逆过程，采用了幂为4的数据置信度项，加快了曲线簇的演化。而且。注意：以下阐述中，能量方程(1)均表示本发明中改进后的能量方程(1)。

${\mathrm{\ξ}}_i(x,y)=\frac{1}{m}\underset{k=1}{\overset{m}{\mathrm{\Σ}}}{(u_k(x,y)-c_{\mathrm{ik}})}^4,$ 1≤i≤N    (2)

式(2)中，

是滤波后的特征矢量影像，m是提取的特征矢量影像的维数，u_k表示第k个分量的特征影像，1≤k≤m；c_ik是二维影像域中第k维特征矢量影像上第i个类别C_i所包含的所有像元的均值，随着迭代次数的增加，c_ik值趋于实际区域的均值，其计算公式(3)为：

$c_{\mathrm{ik}}=\frac{{\∫}_{C_i}{u}_k(x,y)\mathrm{dxdy}}{{\∫}_{C_i}\mathrm{dxdy}}---\left(3\right)$

式(3)中，积分符号∫_Ci表示积分域是第i个类别C_i包含的像元集合；

ii)本发明给出了防模糊分割的约束条件，其表达式如式(4)：

$C_i(x,y,t)=\underset{1\≤i\≤N-1}{\min }\left\{{\mathrm{\ξ}}_i(x,y,t)\right\}---\left(4\right)$

式(4)的含义是：在每个迭代步中，选择ξ_i(x，y)，1≤i≤N-1值最小的第i个水平集函数对应的类标号作为t时刻像元(x，y)最终的类标号C_i。当多水平集演化过程满足收敛条件时，类标号C_i的变化处于静止状态。这个约束条件消除了本方法对高分辨率遥感影像分割和分类结果中的模糊分割现象。

b)对能量方程(1)运用一阶变分法推导出关于曲线簇的演化方程。

由于具体实施方式中的第2步中给出的能量方程(1)形式上与作者A.R.Mansouri提出的能量方程相同，因此曲线簇演化方程的推导都利用了一阶变分原理最小化每一曲线的能量泛函，再利用梯度下降法得出曲线簇的演化方程。具体操作过程参见文献：A.R.Mansouri，A.Mitiche，C.Va′zquez，“Multiregion competition：A level set extension ofregion competition to multiple region image partitioning，”Computer Vision and ImageUnderstanding，no.101，pp.137-150，2006和文献：A.R.Mansouri，A.Mitiche，C.Va′zquez，“Multiregion Competition：A Level Set Extension of Region Competition to MultipleRegion Partitioning of Images and Image Sequences，”Computer Vision and ImageUnderstanding，2004。这样，得出曲线簇的演化方程，即式(6)所示的偏微分方程组：

$\left\{\begin{array}{c}\frac{\∂{\stackrel{\→}{\mathrm{\γ}}}_1}{\∂t}(x,y)=[-{\mathrm{\ξ}}_1\left({\stackrel{\→}{\mathrm{\γ}}}_1(x,y)\right)+{\mathrm{\φ}}_1\left({\stackrel{\→}{\mathrm{\γ}}}_1(x,y)\right)+{\mathrm{\μCur}}_1(x,y)]{\stackrel{\→}{n}}_1(x,y)\\ \·\\ \·\\ \·\\ \frac{\∂{\stackrel{\→}{\mathrm{\γ}}}_i}{\∂t}(x,y)=[-{\mathrm{\ξ}}_i\left({\stackrel{\→}{\mathrm{\γ}}}_i(x,y)\right)+{\mathrm{\φ}}_i\left({\stackrel{\→}{\mathrm{\γ}}}_i(x,y)\right)+{\mathrm{\μCur}}_i(x,y)]{\stackrel{\→}{n}}_i(x,y)\\ \·\\ \·\\ \·\\ \frac{\∂{\stackrel{\→}{\mathrm{\γ}}}_{N-1}}{\∂t}(x,y)=[-{\mathrm{\ξ}}_{N-1}\left({\stackrel{\→}{\mathrm{\γ}}}_{N-1}(x,y)\right)+{\mathrm{\φ}}_{N-1}\left({\stackrel{\→}{\mathrm{\γ}}}_{N-1}(x,y)\right)]+{\mathrm{\μCur}}_{N-1}(x,y)]{\stackrel{\→}{n}}_{N-1}(x,y)\end{array}\right.---\left(5\right)$

式(5)中，

Cur_i分别表示曲线

的单位外法矢量和平均曲率。而

${\mathrm{\φ}}_i(x,y)={\mathrm{\ξ}}_{i+1}(x,y){\mathrm{\χ}}_{{\stackrel{\→}{R}}_{{\stackrel{\→}{\mathrm{\γ}}}_{i=1}}}(x,y)+{\mathrm{\ξ}}_{i+2}(x,y){\mathrm{\χ}}_{{\stackrel{\→}{R}}_{{\stackrel{\→}{\mathrm{\γ}}}_{i+1}^c}}(x,y){\mathrm{\χ}}_{{\stackrel{\→}{R}}_{{\stackrel{\→}{\mathrm{\γ}}}_{i+2}}}(x,y)+\·\·\·$

$+{\mathrm{\ξ}}_{N-1}(x,y){\mathrm{\χ}}_{{\stackrel{\→}{R}}_{{\stackrel{\→}{\mathrm{\γ}}}_{i+1}^c}}(x,y)\·\·\·{\mathrm{\χ}}_{{\stackrel{\→}{R}}_{{\stackrel{\→}{\mathrm{\γ}}}_{N-2}^c}}(x,y){\mathrm{\χ}}_{{\stackrel{\→}{R}}_{{\stackrel{\→}{\mathrm{\γ}}}_{N-1}}}(x,y)$

$+{\mathrm{\ξ}}_N(x,y){\mathrm{\χ}}_{{\stackrel{\→}{R}}_{{\stackrel{\→}{\mathrm{\γ}}}_{i+1}^c}}(x,y)\·\·\·{\mathrm{\χ}}_{{\stackrel{\→}{R}}_{{\stackrel{\→}{\mathrm{\γ}}}_{N-2}^c}}(x,y){\mathrm{\χ}}_{{\stackrel{\→}{R}}_{{\stackrel{\→}{\mathrm{\γ}}}_{N-1}^c}}(x,y)$

c)水平集演化方程的导出

水平集理论对如方程(5)所示的曲线求解问题具有以下优点：曲线拓扑自由改变时无需参数化，而且平面闭合曲线可用水平集曲面的零水平集来表示。设平面闭合曲线簇 

对应的水平集函数矢量为

，平面闭合曲线簇与水平集函数矢量的每一分量一一对应。因此，平面闭合曲线

可用水平集函数

的零水平集 ${\stackrel{\→}{\mathrm{\Φ}}}_i=0$ 表示，而

内、外部区域

和

分别对应于区域 ${\stackrel{\&RightArrow;}{\mathrm{\&Phi;}}}_i>0$ 和 ${\stackrel{\&RightArrow;}{\mathrm{\&Phi;}}}_i<0$ ，则将方程(5)转化成为水平集函数表示的偏微分方程组(6)

$\left\{\begin{array}{c}\frac{{\&PartialD;\mathrm{\&Phi;}}_1}{\&PartialD;t}(x,y)=[-{\mathrm{\&xi;}}_1(x,y)+{\mathrm{\&phi;}}_1(x,y)+{{\mathrm{\&mu;Cur}}_{\mathrm{\&Phi;}}}_1]\left|\right|{\&dtri;\mathrm{\&Phi;}}_1(x,y)\left|\right|\\ \&CenterDot;\\ \&CenterDot;\\ \&CenterDot;\\ \frac{{\&PartialD;\mathrm{\&Phi;}}_i}{\&PartialD;t}(x,y)=[-{\mathrm{\&xi;}}_i(x,y)+{\mathrm{\&phi;}}_i(x,y)+{{\mathrm{\&mu;Cur}}_{\mathrm{\&Phi;}}}_i]\left|\right|{\&dtri;\mathrm{\&Phi;}}_i(x,y)\left|\right|\\ \&CenterDot;\\ \&CenterDot;\\ \&CenterDot;\\ \frac{{\&PartialD;\mathrm{\&Phi;}}_{N-1}}{\&PartialD;t}(x,y)=[-{\mathrm{\&xi;}}_{N-1}(x,y)+{\mathrm{\&phi;}}_{N-1}(x,y)+\mathrm{\&mu;}{\mathrm{Cur}}_{{\text{\&Phi;}}_{N-1}1}]\left|\right|{\&dtri;\mathrm{\&Phi;}}_{N-1}(x,y)\left|\right|\end{array}\right.---\left(6\right)$

式(6)中：函数φ_i(x，y)表达式为：

${\mathrm{\&phi;}}_i(x,y)={\mathrm{\&xi;}}_{i+1}(x,y){\mathrm{\&chi;}}_{\{{\mathrm{\&Phi;}}_{i+1}(x,y,t)>0\}}(x,y)+{\mathrm{\&xi;}}_{i+2}(x,y){\mathrm{\&chi;}}_{\{{\mathrm{\&Phi;}}_{i+1}(x,y,t)\&le;0\}}(x,y){\mathrm{\&chi;}}_{\{{\mathrm{\&Phi;}}_{i+2}(x,y,t)>0\}}(x,y)+\&CenterDot;\&CenterDot;\&CenterDot;$

$+{\mathrm{\&xi;}}_{N-1}(x,y){\mathrm{\&chi;}}_{\{{\mathrm{\&Phi;}}_{i+1}(x,y,t)\&le;0\}}(x,y)\&CenterDot;\&CenterDot;\&CenterDot;{\mathrm{\&chi;}}_{\{{\mathrm{\&Phi;}}_{N-2}(x,y,t)\&le;0\}}(x,y){\mathrm{\&chi;}}_{\{{\mathrm{\&Phi;}}_{N-1}(x,y,t)>0\}}(x,y)$

$+{\mathrm{\&xi;}}_N(x,y){\mathrm{\&chi;}}_{\{{\mathrm{\&Phi;}}_{i+1}(x,y,t)\&le;0\}}(x,y)\&CenterDot;\&CenterDot;\&CenterDot;{\mathrm{\&chi;}}_{\{{\mathrm{\&Phi;}}_{N-2}(x,y,t)\&le;0\}}(x,y){\mathrm{\&chi;}}_{\{{\mathrm{\&Phi;}}_{N-1}(x,y,t)\&le;0\}}(x,y)$

上式中，Cur_Φi表示第i个水平集函数的平均曲率，其表达式为： ${\mathrm{Cur}}_{{\mathrm{\&Phi;}}_i}=\&dtri;\&CenterDot;\frac{{\&dtri;\mathrm{\&Phi;}}_i}{\left|\right|{\&dtri;\mathrm{\&Phi;}}_i\left|\right|}$ (为梯度算子)；χ_{Φi(x，y，t)}(x，y)是第i个水平集函数的示性函数，定义为：

带有约束条件(4)的偏微分方程组(6)就是高分辨率遥感影像的分割和分类模型。该模型的解(即影像分割和分类结果)采用以下方式来表示：第1至第N-1类表示为：C_i(x，y，t)＝{(x，y)∈Ω|Φ_i(x，y，t)＞0}，其中，1≤i≤N-1，第N类表示为： $C_N(x,y,t)=\{(x,y)\&Element;\mathrm{\&Omega;}|\underset{i=1}{\overset{N-1}{\&cap;}}{\mathrm{\&Phi;}}_i(x,y,t)\&le;0\},$ 而分割结果就是所有水平集函数的零水平

集围成的区域；

这种表示方法不同于作者A.R.Mansouri所述方法。此外，A.R.Mansouri所述方法主要应用于计算机视觉领域，作者用分割这一术语来表达遥感领域中影像分类的概念，所以实质上作者的方法是对遥感影像实施了分类，但未涉及分割。

3、第1时相影像的分割和分类

影像的分割和分类就是求解带有约束条件(4)的偏微分方程组(6)的数值解的过程。式(6)是一个带有初、边值条件的连续意义下的偏微分方程组。离散的二维图像中，采用数值解法来逼近真实解，其数值化过程如下：

a)采用有限差分法这种简单而快速的离散化方式，其中，空间偏导是用中心差分法逼近，时间偏导采用前向欧拉差分法逼近；行、列方向上的空间采样间隔Δx＝1，Δy＝1，即一个像元；迭代时间步长Δt依赖于空间采样间隔及水平集函数的演化速度，它的取值必须满足CFL条件，以保持演化过程的稳定性；采用二阶精度的水平集函数重初始化方法，确保较高的准确性及适中的计算代价。水平集函数被定义为一个符号距离函数，作为偏微分方程组(6)的初始条件，目的是保持水平集曲面的光滑性及零水平集求解的方便，其初始化方法如式(7)所示：

式(7)中，d(x，y)是欧氏距离函数；边界条件是Neumann边值条件，目的是使水平集曲面在影像边界Ω处停止演化，其表达式为式(8)：

$\frac{{\&PartialD;\mathrm{\&Phi;}}_i(x,y,t)}{\&PartialD;{\stackrel{\&RightArrow;}{n}}_i}=0$ 当(x，y，t)∈Ω×[0，+∞]    (8)

式(8)表示：在影像边界处水平集函数的法向导数为零。关于欧氏距离函数，CFL条件约束下Δt的取值及水平集函数重初始化方法的具体操作过程参见文献：S.Osher，R.Fedkiw.Level Set Methods and Dynamic Implicit Surfaces.New York：Springer-VerlagNew York Inc.，ch.3，pp.25-37，2003；

b)偏微分方程组(6)的初始条件的设置是非常重要的，它在影像中的空间位置及拓扑形状对曲线簇的演化速度有很大的影响。具体地说，初始轮廓的位置愈接近于最终的分割及分类结果，则所需的时间代价越少；而初始轮廓原则上可为任意形状的平面闭合曲线，但其拓扑形状愈接近于待分割的目标或区域形状，则曲线簇收敛所需时间越短，即分割和分类过程耗费的时间越少。传统的水平集分割和分类方法中均采用人工方式设定水平集函数的初始轮廓，自动化程度低。GIS数据，如地形图、专题图，它包含了地物目标的物理信息、形状大小、空间关系等，在较大程度上反映了地物的分布现状，它可辅助遥感影像进行分割和分类，更详细内容参见文献：王莹，刘敏莺，黄文骞；GIS对遥感影像分类判读的辅助作用；海洋测绘，2002，22(3)：12-15。因此，对于第1时相的高分辨率遥感影像初始条件的设置利用已配准的同一地区的GIS数据，将其视为类别的先验知识来指导初始轮廓的分布，以此来提高分割和分类的自动化程度及准确性。具体方法是：利用GIS数据中已有地物类别的空间位置及属性来画一个初始圆，并将每个初始圆进行符号距离变换作为各个水平集函数的初始条件，符号距离变换的具体方法参见文献S.Osher，R.Fedkiw.Level Set Methods and DynamicImplicit Surfaces.New York：Springer-Verlag New York Inc.，ch.3，pp.25-37，2003。在没有同一地区GIS数据的情况下，对第一时相的高分辨率遥感影像分割和分类的初始条件设置方法为：在影像作用域范围内自动均匀地画任意半径的初始圆，再进行符号距离变换来初始化各个水平集函数，这样，自动化程度提高了，但本方法对初始轮廓的空间位置较敏感，分割和分类精度可能会低于未采用GIS数据辅助初始化时的精度。

c)分割和分类收敛条件：分割和分类过程就是在初、边值条件(7)及约束条件(4)下，迭代地求解与人工时间有关的偏微分方程组(6)，当平面闭合曲线簇收敛，分割和分类过程结束，收敛条件及涉及到的参数参见文献：H.K Zhao，T.F.Chan，B.Merriman，S.Osher，“Variational Level Set Approach to Multiphase Motion，”Journal of ComputationalPhysics，no.0167，pp.179-195，February，1996。

d)分割和分类结果的表达：带有约束条件[公式(4)]的偏微分方程组(6)的解是一组位于影像区域及类边界的曲线簇；采用具体实施内容第2步c)中的阐述的模型解的表示方式对第1时相的高分辨率遥感影像的分割和分类结果进行表示。但采用这种表示方式的前提是偏微分方程组(6)必须具备公式(4)所述的约束条件，这个约束条件的作用是避免模糊分割和分类现象，即避免在像元(x，y)处，Φ_i(x，y，t)＞0且Φ_s(x，y，t)＞0(1≤i，s≤N-1，i≠s)。这一表示方式的含义是：当曲线簇收敛时，若像元(x，y)处第i个水平集函数的符号为正，则该像元的类标号为C_i，但第N类是由影像域Ω中不属于第1至第N-1类中任何一类的像元构成的集合。而最终的分割结果就是所有水平集函数的零水平集所在的曲线围成的区域。见附图2所示，三个水平集函数将影像域Ω分为四个类别，即N＝4，分割和分类结果中，类标号为C₁所在的左上区域中包含着一个类标号为C₂的同质区A₁，类标号为C₂所在的右上区域中包含着一个类标号为C₃的同质区A₂，类标号为C₃所在的左下区域中包含着一个类标号为C₄的同质区A₃，C₄所在的右下区域未包含任何同质区。类似于A₁、A₂、A₃这样的同质区就是分割后的对象，当N取其它值时，分割结果识别方法与此相同。以上阐述表明了多水平集演化模型可同时获得分割和分类结果。最后，根据步骤2建立的模型，采用步骤3中所述的初始和边值条件及离散化方式进行数值求解，来实现第1时相影像的分割和分类，并获得位于分割后区域边界的平面闭合曲线簇(最终曲线簇)的空间位置。

4、第2至第T时相高分辨率遥感影像采用增量式分割和分类

a)对第2至第T时相的高分辨率遥感影像仍采用步骤2中所述的分割和分类模型，但采用新的初始条件及新的影像作用域，新的初始条件设置方法如附图3所示，它表示：将第j时相高分辨率遥感影像分割后的平面闭合曲线簇作为与其相邻的第j+1时相高分辨率遥感影像分割和分类前的初始轮廓；例如，完成第1时相高分辨率遥感影像的分割和分类(其初始化是利用GIS数据辅助完成的)后，将分割后的平面闭合曲线簇的空间位置直接作为第2个时相的高分辨率遥感影像分割和分类前的初始轮廓，当完成第2时相高分辨率遥感影像的分割和分类过程后，将分割后的平面闭合曲线簇空间位置又作为第3个时相的高分辨率遥感影像分割和分类前的初始轮廓，以此类推，对第3时相高分辨率遥感影像之后的各时相影像，其分割和分类前的初始化方式均采用这样的方式；与步骤3中所述的初始化方法相比，这种方法使各水平集函数的初始轮廓空间位置更接近于最终的分割及分类结果，减少了第j+1时相影像分割和分类所需的时间。仍采用步骤3中所述的边界条件、收敛条件、离散化方式、及参数值求解步骤2中所述的分割和分类模型，直到达到新的收敛状态，以实现第j+1时相的高分辨率遥感影像的分割和分类。对第2至第T时相的高分辨率遥感影像的分割和分类结果仍采用具体实施内容第2步c)中的阐述的模型解的表示方式。在以上阐述中，T为变化检测影像的时相总数，其值取决于实际应用中变化检测任务使用的影像时相个数，但必须满足T≥2，即，至少需要两个不同时相的影像才能使用本方法完成变化检测任务，但时相数越大，累积误差越大，导致变化检测精度降低。此外，由于高分辨率遥感影像数据量大，以及获取困难的问题，因此这个T值不能太大。一般取2～10。以上阐述中，j为整数变量，表示第j时相的高分辨率遥感影像，取值范围为是1≤j≤T-1。

b)当第j与第j+1两个相邻时相影像的类别数增加或减少时，需要在原有水平集函数类别数基础上增加新的水平集函数或减少已有的某个水平集函数，增加的水平集函数初始位置及减少的水平集函数类别数可通过以下两种方式确定：1)根据先验知识定；2)通过这两个相邻时相影像的简单差分，变化区域作为新增类别的先验知识来指导新增的水平集函数的初始轮廓的设置。

5、相邻两个时相影像的变化检测

当完成了第j与第j+1两个相邻时相影像的分割和分类后，开始这两时相影像的变化检测。例如，T＝3时，获得了第1与第2两个相邻时相影像的分割和分类结果后，开始对这两个影像进行变化检测，当完成了第3时相的影像的分割和分类，开始对第2与第3两个相邻时相的影像进行变化检测；参见附图4所示，具体的检测方法是：以其中任意一个影像的分类结果为参考，对比这两相邻时相影像的分类结果，若另一个影像的分类结果中某个分割对象全部或部分的类标号发生了改变，则表明该对象的地物类型发生了变化，并输出所有变化像元的空间位置及类别属性。

6、返回第4步，循环执行步骤4至6的过程，直到完成了所有的T个时相影像的分割和分类及变化检测任务。

Claims (3)

1.一种高分辨率遥感影像的分割和分类及变化检测一体化方法，其特征包括以下步骤：

一、高分辨率遥感影像的预处理：

a)对各时相高分辨率遥感影像分别进行辐射校正；并对辐射校正后的各时相高分辨率遥感影像进行配准，然后再与同一地区的GIS数据配准；

b)对辐射校正和配准之后的高分辨率遥感影像提取纹理特征及光谱特征，组成m维特征矢量影像

其中，m是提取的特征矢量影像的维数，且m≥1；u_k表示第k个分量的特征影像，1≤k≤m；然后利用矢量总变分的多尺度异向扩散技术进行滤波；

二、对预处理之后的特征矢量影像建立分割和分类模型：

首先，建立多水平集演化方程组，即偏微分方程组：

$\left\{\begin{array}{c}\frac{\&PartialD;{\mathrm{\&Phi;}}_1}{\&PartialD;t}(x,y)=[-{\mathrm{\&xi;}}_1(x,y)+{\mathrm{\&phi;}}_1(x,y)+\mathrm{\&mu;Cu}{r}_{{\mathrm{\&Phi;}}_1}]\left|\right|\&dtri;{\mathrm{\&Phi;}}_1(x,y)\left|\right|\\ .\\ .\\ .\\ \frac{\&PartialD;{\mathrm{\&Phi;}}_i}{\&PartialD;t}(x,y)=[-{\mathrm{\&xi;}}_i(x,y)+{\mathrm{\&phi;}}_i(x,y)+\mathrm{\&mu;Cu}{r}_{{\mathrm{\&Phi;}}_i}]\left|\right|\&dtri;{\mathrm{\&Phi;}}_i(x,y)\left|\right|\\ .\\ .\\ .\\ \frac{\&PartialD;{\mathrm{\&Phi;}}_{N-1}}{\&PartialD;t}(x,y)=[-{\mathrm{\&xi;}}_{N-1}(x,y)+{\mathrm{\&phi;}}_{N-1}(x,y)+\mathrm{\&mu;Cu}{r}_{{\mathrm{\&Phi;}}_{N1}1}]\left|\right|\&dtri;{\mathrm{\&Phi;}}_{N-1}(x,y)\left|\right|\end{array}\right.---\left(1\right)$

式(1)中，

是隐含表示N-1维水平集函数矢量，N是高分辨率遥感影像包含的地物类别个数；函数ξ_i(x，y)是二维影像域Ω中第i个类别C_i的统计特征，计算公式为：

${\mathrm{\&xi;}}_i(x,y)=\frac{1}{m}\underset{k=1}{\overset{m}{\mathrm{\&Sigma;}}}{(u_k(x,y)-c_{\mathrm{ik}})}^4,1\&le;i\&le;N---\left(2\right)$

式(2)中，c_ik是二维影像域Ω中第k维特征矢量影像上第i个类别C_i的均值，计算公式为：

$c_{\mathrm{ik}}=\frac{{\&Integral;}_{C_i}{u}_k(x,y)\mathrm{dxdy}}{{\&Integral;}_{C_i}\mathrm{dxdy}}---\left(3\right)$

式(3)中，积分符号∫_Ci表示积分域是第i个类别C_i包含的像元集合；式(1)中，函数φ_i(x，y)的计算公式为：

${\mathrm{\&phi;}}_i(x,y)={\mathrm{\&xi;}}_{i+1}(x,y){\mathrm{\&chi;}}_{\{{\mathrm{\&Phi;}}_{i+1}(x,y,t)>0\}}(x,y)+{\mathrm{\&xi;}}_{i+2}(x,y){\mathrm{\&chi;}}_{\{{\mathrm{\&Phi;}}_{i+1}(x,y,t)\&le;0\}}(x,y){\mathrm{\&chi;}}_{\{{\mathrm{\&Phi;}}_{i+2}(x,y,t)>0\}}(x,y)+...$

$+{\mathrm{\&xi;}}_{N-1}(x,y){\mathrm{\&chi;}}_{\{{\mathrm{\&Phi;}}_{i+1}(x,y,t)\&le;0\}}(x,y)...{\mathrm{\&chi;}}_{\{{\mathrm{\&Phi;}}_{N-2}(x,y,t)\&le;0\}}(x,y){\mathrm{\&chi;}}_{\{{\mathrm{\&Phi;}}_{N-1}(x,y,t)>0\}}(x,y)---\left(4\right)$

$+{\mathrm{\&xi;}}_N(x,y){\mathrm{\&chi;}}_{\{{\mathrm{\&Phi;}}_{i+1}(x,y,t)\&le;0\}}(x,y)...{\mathrm{\&chi;}}_{\{{\mathrm{\&Phi;}}_{N-2}(x,y,t)\&le;0\}}(x,y){\mathrm{\&chi;}}_{\{{\mathrm{\&Phi;}}_{N-1}(x,y,t)\&le;0\}}(x,y)$

式(1)中，Cur_Φi是第i个水平集函数的平均曲率，计算公式为： ${\mathrm{Cur}}_{{\mathrm{\&Phi;}}_i}=\&dtri;\&CenterDot;\frac{\&dtri;{\mathrm{\&Phi;}}_i}{\left|\right|{\&dtri;\mathrm{\&Phi;}}_i\left|\right|},$ 为梯度算子；参数μ是平均曲率的权值；χ_{Φi(x，y，t)}(x，y)是第i个水平集函数的示性函数，定义为：

给定约束条件：

$C_i(x,y,t)=\underset{1\&le;i\&le;N-1}{\min }\left\{{\mathrm{\&xi;}}_i(x,y,t)\right\}---\left(5\right)$

式(5)中，C_i(x，y，t)表示第t次迭代中，ξ_i(x，y，t)值最小的第i个水平集函数对应的像元(x，y)的类标号；当多水平集演化过程满足收敛条件时，类标号C_i的变化处于静止状态；

带有约束条件(5)的偏微分方程组(1)就是高分辨率遥感影像的分割和分类模型，该模型的解采用以下方式来表示：第1至第N-1类表示为：C_i(x，y，t)＝{(x，y)∈Ω|Φ_i(x，y，t)＞0}，其中，1≤i≤N-1，第N类表示为： $C_N(x,y,t)=\{(x,y)\&Element;\mathrm{\&Omega;}|\underset{i=1}{\overset{N-1}{\mathrm{\&cap;}}}{\mathrm{\&Phi;}}_i(x,y,t)\&le;0\},$ 而分割结果就是所有水平集函数的零水平集围成的区域；

三、第1个时相高分辨率遥感影像的分割和分类：

a)利用有限差分法对分割和分类模型进行离散化，边界条件是Neumann边值条件；输入式(1)-(3)中参数：影像包含的地物类别个数N取正整数，且N≥2，这样，分割和分类模型的所需的水平集函数个数为N-1；平均曲率的权值μ＝(0.1)^p×255²，p为整数，取值范围为0～5；行、列方向上的空间采样间隔Δx＝1，Δy＝1，即一个像元；迭代时间步长Δt依赖于空间采样间隔及水平集函数的演化速度，它的取值必须满足CFL条件；

b)在对第1时相的高分辨率遥感影像水平集分割和分类中，分割和分类模型的初始条件设置方法如下：在配准后的该地区GIS数据辅助下，设置每一个水平集函数的初始轮廓，具体方法：利用GIS数据中每一个地物类别的空间位置及属性来画一个初始圆，并将每个初始圆进行符号距离变换作为各个水平集函数的初始条件，采用步骤三中a)所述的边界条件和离散化方式及参数值来实现第1时相影像的水平集分割和分类，分割和分类结束时得到分割后的平面闭合曲线簇；

四、第2至第T时相的高分辨率遥感影像采用增量式分割和分类：对第2至第T时相的高分辨率遥感影像仍采用步骤二中所述的分割和分类模型，但采用新的初始条件及新的影像作用域，新的初始条件设置方法如下：将第j时相高分辨率遥感影像分割后的平面闭合曲线簇作为与其相邻的第j+1时相高分辨率遥感影像分割和分类前的初始轮廓；仍采用步骤三中a)所述的边界条件和离散化方式及参数值来实现第j+1时相的高分辨率遥感影像的分割和分类，分割和分类结束时得到分割后的平面闭合曲线簇；其中：T为变化检测影像的时相总数；j为整数变量，表示第j时相的高分辨率遥感影像，取值范围为是1≤j≤T-1；

五、第j与第j+1两个相邻时相高分辨率遥感影像的变化检测：从j＝1开始，当第j与第j+1时相高分辨率遥感影像的分割和分类过程结束后，开始对这两个相邻时相高分辨率遥感影像进行变化检测，具体方法是：对比这两个相邻时相高分辨率遥感影像的分类结果，以其中任意一个时相高分辨率遥感影像的分类结果为参考，在另一时相高分辨率遥感影像的分类结果中，以分割后的对象为比较单元，若任意一个对象的全部像元或部分像元的类标号发生了改变，则表明该对象的地物类型发生了变化，并输出所有变化像元的空间位置及类别属性；

六、返回第四步，循环执行步骤四至六的过程，直到完成了所有的T个时相影像的分割和分类及变化检测任务。

2.根据权利要求1所述的方法，其特征包括：对大小为512×512以上的T个不同时相的高分辨率遥感影像分割和分类及变化检测采用以下步骤：按步骤一进行预处理；再采用小波多分辨率分解；按照步骤三至六进行分割和分类及变化检测，直到完成了所有的T个不同时相影像的分割和分类及变化检测任务。

3.根据权利要求1或2所述的方法，其特征在于：在没有同一地区的GIS数据的情况下，步骤二中所述的分割和分类模型的初始条件设置方法为：在影像作用域范围内自动均匀地画任意半径的圆，再对这些圆进行符号距离变换来初始化各个水平集函数。

Images