CN102663815A - 一种基于水平集的lod2建筑物模型构建方法 - Google Patents

Abstract

一种基于水平集的LOD2建筑物模型构建方法，属于应用水平集算法对DSM数据的分割处理领域。为了解决现有的基于二维图像的建筑物模型构建方法中存在的顶面结构不够细致而导致构建精度不高的问题。具体过程为：提取建筑物轮廓掩膜Ω_m，选取DSM数据，配入到统一坐标系下；获取建筑物顶面数据T；获取建筑物顶面数据T的特征空间；进行多相水平集分割，获取子区域；提取每个子区域的点集，检测每个片元的边界点，获取建筑物每个基元的角点的图像坐标；建立建筑物顶面数据T的拓扑结构；根据航拍可见光图像提取建筑物表面的纹理数据，与建筑物不同基元相对应，完成LOD2建筑物模型的构建。用于LOD2级别的大规模建筑物的三维构建任务。

Description

一种基于水平集的LOD2建筑物模型构建方法

技术领域

本发明属于应用水平集算法对DSM数据的分割处理领域，具体涉及一种基于水平集的LOD2建筑物模型构建方法。

背景技术

随着对数字城市需求的日益提高，城市和郊区人造地物的高精度模型化能够使人们参与到城市的规划设计与建设，对经济的发展和避难措施的策划提出参考性的建议。游戏娱乐、汽车导航、电子商务等领域的飞速发展促进了人们对数字城市的迫切需求。三维建筑物模型通常采用多尺度的表达方式，以满足不同的应用需求。根据CityGML的定义，LOD1(Level Of Detail 1)是指仅有平面屋顶的建筑物模型，LOD1是第一级尺度，LOD2(LevelOf Detail 2)建筑物模型具有不同类型的建筑物屋顶和真实的表面纹理。目前，许多研究集中在LOD1建筑物模型的构建方法上，主要研究方向是对具有精细几何结构以及表面纹理的建筑物模型的构建方法。

在过去的几十年里，三维建筑物模型的构建技术可以使用多种数据源，尤其是DSM(Digital Surface Model，数字表面模型)的应用愈加广泛，通过多角度航拍图像或LiDAR(Light Detection And Ranging)数据处理所得的高精度DSM数据，其分辨率高达0.3m，可提供丰富的三维信息，从而构建出精细的屋顶几何结构。但是现有技术还没有从DSM数据中提取复杂建筑物模型的顶面基元以及拓扑结构、实现较高精度的LOD2建筑物模型构建的方法。

发明内容

本发明是为了解决现有的基于二维图像的建筑物模型构建方法中存在的顶面结构不够细致而导致构建精度不高的问题，提出的一种基于水平集的LOD2建筑物模型构建方法。

本发明的一种基于水平集的LOD2建筑物模型构建方法，其具体过程为：

步骤1、提取建筑物轮廓作为建筑物轮廓掩膜Ω_m，同时根据相对应的建筑物区域的地理坐标选取DSM数据，通过自动配准或人工配准的方式，将建筑物轮廓掩膜Ω_m和DSM数据配入到统一的坐标系下；

步骤2、根据步骤1选取的建筑物轮廓掩膜Ω_m和DSM数据获取建筑物顶面数据T；

步骤3、根据步骤2获得的建筑物顶面数据T获取建筑物顶面数据T的特征空间；

步骤4、根据步骤3获取的特征空间将建筑物区域进行多相水平集分割，获取子区域；

步骤5、根据多相水平集的定义分别提取步骤4获得的每个子区域的点集，然后检测建筑物每个片元的边界点，获取建筑物每个基元的角点的图像坐标；

步骤6、根据步骤5获得的的建筑物角点的图像坐标，与步骤1中选取的DSM数据中的地理坐标相对应，建立建筑物顶面数据T的拓扑结构；

步骤7、根据步骤6获取的建筑物顶面数据T的拓扑结构，同时根据相应的建筑物区域的航拍可见光图像提取建筑物表面的纹理数据，与建筑物的不同基元相对应，完成LOD2建筑物模型的构建。

本发明基于多相多通道水平集框架，将局部区域信息引入该框架，以DSM数据点的法向量作为特征向量进行水平集分割，分离出建筑物顶面的不同基元，从而实现顶面几何结构的构建。

本发明的优势在于能够从DSM数据中提取复杂建筑物模型的顶面基元以及拓扑结构，实现较高精度的LOD2建筑物模型构建。

本发明可用于LOD2级别的大规模建筑物的三维构建任务。

附图说明

图1是本发明的建筑物模型构建方法流程示意图；图2是多想水平集分割过程的流程示意图；图3是本发明所采用的多相水平集方法的示意图；图4和图5是本发明在水平集框架中引入局部化区域信息方法的示意图。

具体实施方式

具体实施方式一、结合图1说明本实施方式，一种基于水平集的LOD2建筑物模型构建方法，其具体过程为：

步骤1、提取建筑物轮廓作为建筑物轮廓掩膜Ω_m，同时根据相对应的建筑物区域的地理坐标选取DSM数据，通过自动配准或人工配准的方式，将建筑物轮廓掩膜Ω_m和DSM数据配入到统一的坐标系下；

步骤2、根据步骤1选取的建筑物轮廓掩膜Ω_m和DSM数据获取建筑物顶面数据T；

步骤3、根据步骤2获得的建筑物顶面数据T获取建筑物顶面数据T的特征空间；

步骤4、根据步骤3获取的特征空间将建筑物区域进行多相水平集分割，获取子区域；

步骤5、根据多相水平集的定义分别提取步骤4获得的每个子区域的点集，然后检测建筑物每个片元的边界点，获取建筑物每个基元的角点的图像坐标；

步骤6、根据步骤5获得的的建筑物角点的图像坐标，与步骤1中选取的DSM数据中的地理坐标相对应，建立建筑物顶面数据T的拓扑结构；

步骤7、根据步骤6获取的建筑物顶面数据T的拓扑结构，同时根据相应的建筑物区域的航拍可见光图像提取建筑物表面的纹理数据，与建筑物的不同基元相对应，完成LOD2建筑物模型的构建。

具体实施方式二、本实施方式与具体实施方式一的区别在于，本实施方式是对具体实施方式一中步骤2的进一步限定，所述步骤2中获取建筑物顶面数据T的方法为：用建筑物轮廓掩膜Ω_m与配准后的DSM数据进行对应点相乘，然后用获得的结果去除建筑物轮廓外点集，获得建筑物的顶面数据T。

具体实施方式三、本实施方式与具体实施方式一的区别在于，本实施方式是对具体实施方式一中步骤3的进一步限定，所述步骤3中获取建筑物顶面数据T的特征空间的具体过程为：对建筑物顶面数据T应用Delaunay三角形剖分算法计算每一个三角形的法向量，然后对法向量进行归一化处理，将建筑物顶面数据T中的每一个点p_i的法向量N_i用共用点p_i的相邻三角形的法向量的平均值来表示，结果作为建筑物顶面数据T的特征空间。

具体实施方式四、本实施方式与具体实施方式一的区别在于，本实施方式是对具体实施方式一中步骤4的进一步限定，所述步骤4中将建筑物区域进行多相水平集分割的具体过程为：

步骤4-1、在待分割数据I的定义域Ω内，以两条等间距、等半径的圆圈为初始曲线C₁和C₂，把定义域Ω划分为四个区域，分别将曲线C₁和C₂转化成水平集方程φ₁和φ₂的形式，令X表示任意点(x，y)，用符号距离函数表示一个点X的水平集方程φ₁(X)和φ₂(X)为：

其中，D(X，C₁)表示点X与曲线C₁之间的Euclidean距离，

其中，D(X，C₂)表示点X与曲线C₂之间的Euclidean距离；

步骤4-2、分别计算两个水平集方程φ₁(X)和φ₂(X)曲线的长度分量和面积分量，长度分量分别为：

length(C₁)＝∫_Ω|▽H(φ₁(X))|dX＝∫_Ωδ(φ₁(X))|▽φ₁(X)|dX

length(C₂)＝∫_Ω|▽H(φ₂(X))|dX＝∫_Ωδ(φ₂(X))|▽φ₂(X)|dX

面积分量分别为：

area(inside(C₁))＝∫_ΩH(φ₁(X))dX

area(inside(C₂))＝∫_ΩH(φ₂(X))dX

其中， $\mathrm{\&delta;}\left(\mathrm{\&phi;}\right)=\left\{\begin{array}{cc}0,& \left|\mathrm{\&phi;}\right|>\mathrm{\&epsiv;}\\ \frac{1}{2\mathrm{\&epsiv;}}(1+\cos \left(\frac{\mathrm{\&pi;\&phi;}}{\mathrm{\&epsiv;}}\right)),& \left|\mathrm{\&phi;}\right|<\mathrm{\&epsiv;}\end{array}\right.,$ $H\left(\mathrm{\&phi;}\right)=\left\{\begin{array}{cc}1,& \mathrm{\&phi;}>\mathrm{\&epsiv;}\\ 0,& \mathrm{\&phi;}<-\mathrm{\&epsiv;}\\ \frac{1}{2}(1+\frac{\mathrm{\&phi;}}{\mathrm{\&epsiv;}}+\frac{1}{\mathrm{\&pi;}}\sin \left(\frac{\mathrm{\&pi;\&phi;}}{\mathrm{\&epsiv;}}\right)),& \left|\mathrm{\&phi;}\right|<\mathrm{\&epsiv;}\end{array}\right.,$

φ为φ₁或φ₂；

步骤4-3、在多相水平集框架中引入局部化区域函数B(X，Y)，设置局部化区域函数B(X，Y)的半径参数r，以局部化区域参数B(X，Y)为局部区域掩膜，计算以X为中心、r为半径的圆圈内的特征空间的平均值：

$c_{++}=\frac{{\&Integral;}_{{\mathrm{\&Omega;}}_Y}B(X,Y)\&CenterDot;H\left({\mathrm{\&phi;}}_1\right)\&CenterDot;H\left({\mathrm{\&phi;}}_2\right)\&CenterDot;I\left(Y\right)\mathrm{dY}}{{\&Integral;}_{{\mathrm{\&Omega;}}_Y}B(X,Y)\&CenterDot;H\left({\mathrm{\&phi;}}_1\right)\&CenterDot;H\left({\mathrm{\&phi;}}_2\right)\mathrm{dY}}$

$c_{+-}=\frac{{\&Integral;}_{{\mathrm{\&Omega;}}_Y}B(X,Y)\&CenterDot;H\left({\mathrm{\&phi;}}_1\right)\&CenterDot;(1-H\left({\mathrm{\&phi;}}_2\right))\&CenterDot;I\left(Y\right)\mathrm{dY}}{{\&Integral;}_{{\mathrm{\&Omega;}}_Y}B(X,Y)\&CenterDot;H\left({\mathrm{\&phi;}}_1\right)\&CenterDot;(1-H\left({\mathrm{\&phi;}}_2\right))\mathrm{dY}}$

$c_{-+}=\frac{{\&Integral;}_{{\mathrm{\&Omega;}}_Y}B(X,Y)\&CenterDot;(1-H\left({\mathrm{\&phi;}}_1\right))\&CenterDot;H\left({\mathrm{\&phi;}}_2\right)\&CenterDot;I\left(Y\right)\mathrm{dY}}{{\&Integral;}_{{\mathrm{\&Omega;}}_Y}B(X,Y)\&CenterDot;(1-H\left({\mathrm{\&phi;}}_1\right))\&CenterDot;H\left({\mathrm{\&phi;}}_2\right)\mathrm{dY}}$

$c_{--}=\frac{{\&Integral;}_{{\mathrm{\&Omega;}}_Y}B(X,Y)\&CenterDot;(1-H\left({\mathrm{\&phi;}}_1\right))\&CenterDot;(1-H\left({\mathrm{\&phi;}}_2\right))\&CenterDot;I\left(Y\right)\mathrm{dY}}{{\&Integral;}_{{\mathrm{\&Omega;}}_Y}B(X,Y)\&CenterDot;(1-H\left({\mathrm{\&phi;}}_1\right))\&CenterDot;(1-H\left({\mathrm{\&phi;}}_2\right))\mathrm{dY}}$

I(Y)表示待分割数据；

步骤4-4、设置权重因子μ₁，μ₂，ν₁，ν₂，λ₁，λ₂，λ₃，λ₄，定义多相水平集能量泛函：

E_total＝μ₁·E₁+μ₂·E₂+ν₁·E₃+ν₂·E₄+λ₁·E₅+λ₂·E₆+λ₃·E₇+λ₄·E₈；

步骤4-5、以局部化区域函数B(X，Y)为局部区域掩膜，计算以X为中心、r为半径的圆圈内的能量分量，

水平集方程φ₁的能量泛函表示为：

$E_{\mathrm{total}}\left({\mathrm{\&phi;}}_1\right)={\&Integral;}_{{\mathrm{\&Omega;}}_X}\mathrm{\&delta;}\left({\mathrm{\&phi;}}_1\right){\&Integral;}_{{\mathrm{\&Omega;}}_Y}B(X,Y)[{\mathrm{\&lambda;}}_1{\left|\right|u-c_{++}\left|\right|}^{2}H\left({\mathrm{\&phi;}}_1\right)H\left({\mathrm{\&phi;}}_2\right)+{\mathrm{\&lambda;}}_2{\left|\right|u-c_{+-}\left|\right|}^{2}H\left({\mathrm{\&phi;}}_1\right)(1-H\left({\mathrm{\&phi;}}_2\right))$

$+{\mathrm{\&lambda;}}_3{\left|\right|u-c_{-+}\left|\right|}^2(1-H\left({\mathrm{\&phi;}}_1\right))H\left({\mathrm{\&phi;}}_2\right)+{\mathrm{\&lambda;}}_4{\left|\right|u-c_{--}\left|\right|}^2(1-H\left({\mathrm{\&phi;}}_1\right))(1-H\left({\mathrm{\&phi;}}_2\right))]\mathrm{dYdX}$

$+{\mathrm{\&mu;}}_1\&CenterDot;\mathrm{length}\left(C_1\right)+{\mathrm{\&mu;}}_2\&CenterDot;\mathrm{length}\left(C_2\right)+{\mathrm{\&nu;}}_1\&CenterDot;\mathrm{area}\left(\mathrm{inside}\left(C_1\right)\right)+{\mathrm{\&nu;}}_2\&CenterDot;\mathrm{area}\left(\mathrm{inside}\left(C_2\right)\right)$

水平集方程φ₂的能量泛函表示为：

$E_{\mathrm{total}}\left({\mathrm{\&phi;}}_2\right)={\&Integral;}_{{\mathrm{\&Omega;}}_X}\mathrm{\&delta;}\left({\mathrm{\&phi;}}_2\right){\&Integral;}_{{\mathrm{\&Omega;}}_Y}B(X,Y)[{\mathrm{\&lambda;}}_1{\left|\right|u-c_{++}\left|\right|}^{2}H\left({\mathrm{\&phi;}}_1\right)H\left({\mathrm{\&phi;}}_2\right)+{\mathrm{\&lambda;}}_2{\left|\right|u-c_{+-}\left|\right|}^{2}H\left({\mathrm{\&phi;}}_1\right)(1-H\left({\mathrm{\&phi;}}_2\right))$

$+{\mathrm{\&lambda;}}_3{\left|\right|u-c_{-+}\left|\right|}^2(1-H\left({\mathrm{\&phi;}}_1\right))H\left({\mathrm{\&phi;}}_2\right)+{\mathrm{\&lambda;}}_4{\left|\right|u-c_{--}\left|\right|}^2(1-H\left({\mathrm{\&phi;}}_1\right))(1-H\left({\mathrm{\&phi;}}_2\right))]\mathrm{dYdX}$

$+{\mathrm{\&mu;}}_1\&CenterDot;\mathrm{length}\left(C_1\right)+{\mathrm{\&mu;}}_2\&CenterDot;\mathrm{length}\left(C_2\right)+{\mathrm{\&nu;}}_1\&CenterDot;\mathrm{area}\left(\mathrm{inside}\left(C_1\right)\right)+{\mathrm{\&nu;}}_2\&CenterDot;\mathrm{area}\left(\mathrm{inside}\left(C_2\right)\right);$

步骤4-6、根据梯度下降法，分别计算步骤4-5获得的能量泛函E_total(φ₁)和E_total(φ₂)的演化方程：

$\frac{{\&PartialD;\mathrm{\&phi;}}_1}{\&PartialD;t}=\mathrm{\&delta;}\left({\mathrm{\&phi;}}_1\right)\{{\mathrm{\&mu;}}_1\&CenterDot;\mathrm{div}\left(\frac{{\&dtri;\mathrm{\&phi;}}_1}{\left|\&dtri;\mathrm{\&phi;}\right|}\right)-{\mathrm{\&nu;}}_1+{\&Integral;}_{{\mathrm{\&Omega;}}_Y}B(X,Y)\mathrm{\&delta;}\left({\mathrm{\&phi;}}_1\right)[({\mathrm{\&lambda;}}_3{\left|\right|u-c_{-+}\left|\right|}^2-{\mathrm{\&lambda;}}_1{\left|\right|u-c_{++}\left|\right|}^2)H\left({\mathrm{\&phi;}}_2\right)$

$-({\mathrm{\&lambda;}}_2{\left|\right|u-c_{+-}\left|\right|}^2-{\mathrm{\&lambda;}}_4{\left|\right|u-c_{--}\left|\right|}^2)(1-H\left({\mathrm{\&phi;}}_2\right))\left]\mathrm{dY}\right\}$

$\frac{{\&PartialD;\mathrm{\&phi;}}_2}{\&PartialD;t}=\mathrm{\&delta;}\left({\mathrm{\&phi;}}_2\right)\{{\mathrm{\&mu;}}_2\&CenterDot;\mathrm{div}\left(\frac{{\&dtri;\mathrm{\&phi;}}_2}{\left|{\&dtri;\mathrm{\&phi;}}_2\right|}\right)-{\mathrm{\&nu;}}_2+{\&Integral;}_{{\mathrm{\&Omega;}}_Y}B(X,Y)\mathrm{\&delta;}\left({\mathrm{\&phi;}}_2\right)[({\mathrm{\&lambda;}}_2{\left|\right|u-c_{+-}\left|\right|}^2-{\mathrm{\&lambda;}}_1{\left|\right|u-c_{++}\left|\right|}^2)H\left({\mathrm{\&phi;}}_1\right)$

$-({\mathrm{\&lambda;}}_3{\left|\right|u-c_{-+}\left|\right|}^2-{\mathrm{\&lambda;}}_4{\left|\right|u-c_{--}\left|\right|}^2)(1-H\left({\mathrm{\&phi;}}_1\right))\left]\mathrm{dY}\right\}$

步骤4-7、判断能量泛函是否达到了最小值，如果是则完成将待分割区域进行多相水平集分割的过程，将待分割区域分割成四个子区域，如果否，则返回步骤4-2，迭代更新水平集方程φ₁和φ₂。

本实施方式中，步骤4-4中设置权重因子的经验值为μ₁＝μ₂＝0.2，ν₁＝ν₂＝λ₁＝λ₂＝λ₃＝λ₄＝1。步骤4-5的作用是减少冗余信息的影响，并加速计算的速度。设置局部化区域函数B(X，Y)的半径参数r的经验值是r＝5。

具体实施方式五、本实施方式与具体实施方式四的区别在于，本实施方式是对具体实施方式四中步骤4-7的进一步限定，所述步骤4-7中判断能量泛函是否达到了最小值的依据为：通过设定阈值或设定迭代次数。

具体实施方式六、本实施方式与具体实施方式一的区别在于，本实施方式是对具体实施方式一中步骤5的进一步限定，所述步骤5中获取建筑物每个基元的角点的图像坐标具体方法为：

步骤5-1、令水平集方程φ₁＝0，φ₂＝0，获得建筑物顶面几何结构边缘的二值图；

步骤5-2、对步骤5-1获得的建筑物顶面几何结构边缘二值图应用Harris算子检测建筑物结构每个基元的角点的图像坐标，Harris算子为R＝det(C)-ktr²(C)，其中

表示图像点z的灰度在u的偏导，表示图像点z的灰度在v的偏导，I_uv(z)表示图像点z的灰度在u和v的二阶混合偏导。本实施方式中，

和I_uv(z)分别是图像点z的灰度在u和v的偏导以及二阶混合偏导，一般取0.04～0.06。

Claims (6)

1.一种基于水平集的LOD2建筑物模型构建方法，其特征在于：该方法的具体过程为：

步骤1、提取待构建模型的建筑物轮廓作为该建筑物轮廓掩膜Ω_m，同时根据相对应的建筑物区域的地理坐标选取DSM数据，通过自动配准或人工配准的方式，将建筑物轮廓掩膜Ω_m和DSM数据配入到统一的坐标系下；

步骤2、根据步骤1选取的建筑物轮廓掩膜Ω_m和DSM数据获取建筑物顶面数据T；

步骤3、根据步骤2获得的建筑物顶面数据T获取建筑物顶面数据T的特征空间；步骤4、根据步骤3获取的特征空间将建筑物区域进行多相水平集分割，获取子区域；

步骤5、根据多相水平集的定义分别提取步骤4获得的每个子区域的点集，然后检测建筑物每个片元的边界点，获取建筑物每个基元的角点的图像坐标；

步骤6、根据步骤5获得的建筑物角点的图像坐标，与步骤1中选取的DSM数据中的地理坐标相对应，建立建筑物顶面数据T的拓扑结构；

步骤7、根据步骤6获取的建筑物顶面数据T的拓扑结构，同时根据相应的建筑物区域的航拍可见光图像提取建筑物表面的纹理数据，与建筑物的不同基元相对应，完成LOD2建筑物模型的构建。

2.根据权利要求1所述的一种基于水平集的LOD2建筑物模型构建方法，其特征在于：所述步骤2中获取建筑物顶面数据T的方法为：用建筑物轮廓掩膜Ω_m与配准后的DSM数据进行对应点相乘，然后用获得的结果去除建筑物轮廓外点集，获得建筑物的顶面数据T。

3.根据权利要求1所述的一种基于水平集的LOD2建筑物模型构建方法，其特征在于：所述步骤3中获取建筑物顶面数据T的特征空间的具体过程为：对建筑物顶面数据T应用Delaunay三角形剖分算法计算每一个三角形的法向量，然后对所有法向量进行归一化处理，将建筑物顶面数据T中的每一个点p_i的法向量N_i用共用点p_i的相邻三角形的法向量的平均值来表示，结果作为建筑物顶面数据T的特征空间。

4.根据权利要求1所述的一种基于水平集的LOD2建筑物模型构建方法，其特征在于：所述步骤4中将建筑物区域进行多相水平集分割的具体过程为：

步骤4-1、在待分割数据I的定义域Ω内，以两条等间距、等半径的圆圈为初始曲线C₁和C₂，把定义域Ω划分为四个区域，分别将曲线C₁和C₂转化成水平集方程φ₁和φ₂的形式，令X表示任意点(x，y)，用符号距离函数表示一个点X的水平集方程φ₁(X)和φ₂(X)为：

其中，D(X，C₁)表示点X与曲线C₁之间的Euclidean距离，

其中，D(X，C₂)表示点X与曲线C₂之间的Euclidean距离；

步骤4-2、分别计算两个水平集方程φ₁(X)和φ₂(X)曲线的长度分量和面积分量，长度分量分别为：

length(C₁)＝∫_Ω|▽H(φ₁(X))|dX＝∫_Ωδ(φ₁(X))|▽φ₁(X)|dX

length(C₂)＝∫_Ω|▽H(φ₂(X))|dX＝∫_Ωδ(φ₂(X))|▽φ₂(X)|dX

面积分量分别为：

area(inside(C₁))＝∫_ΩH(φ₁(X))dX

area(inside(C₂))＝∫_ΩH(φ₂(X))dX

其中， $\mathrm{\&delta;}\left(\mathrm{\&phi;}\right)=\left\{\begin{array}{cc}0,& \left|\mathrm{\&phi;}\right|>\mathrm{\&epsiv;}\\ \frac{1}{2\mathrm{\&epsiv;}}(1+\cos \left(\frac{\mathrm{\&pi;\&phi;}}{\mathrm{\&epsiv;}}\right)),& \left|\mathrm{\&phi;}\right|<\mathrm{\&epsiv;}\end{array}\right.,$ $H\left(\mathrm{\&phi;}\right)=\left\{\begin{array}{cc}1,& \mathrm{\&phi;}>\mathrm{\&epsiv;}\\ 0,& \mathrm{\&phi;}<-\mathrm{\&epsiv;}\\ \frac{1}{2}(1+\frac{\mathrm{\&phi;}}{\mathrm{\&epsiv;}}+\frac{1}{\mathrm{\&pi;}}\sin \left(\frac{\mathrm{\&pi;\&phi;}}{\mathrm{\&epsiv;}}\right)),& \left|\mathrm{\&phi;}\right|<\mathrm{\&epsiv;}\end{array}\right.,$

φ为φ₁或φ₂；

步骤4-3、在多相水平集框架中引入局部化区域函数B(X，Y)，设置局部化区域函数B(X，Y)的半径参数r，以局部化区域参数B(X，Y)为局部区域掩膜，计算以X为中心、r为半径的圆圈内的特征空间的平均值：

$c_{++}=\frac{{\&Integral;}_{{\mathrm{\&Omega;}}_Y}B(X,Y)\&CenterDot;H\left({\mathrm{\&phi;}}_1\right)\&CenterDot;H\left({\mathrm{\&phi;}}_2\right)\&CenterDot;I\left(Y\right)\mathrm{dY}}{{\&Integral;}_{{\mathrm{\&Omega;}}_Y}B(X,Y)\&CenterDot;H\left({\mathrm{\&phi;}}_1\right)\&CenterDot;H\left({\mathrm{\&phi;}}_2\right)\mathrm{dY}}$

$c_{+-}=\frac{{\&Integral;}_{{\mathrm{\&Omega;}}_Y}B(X,Y)\&CenterDot;H\left({\mathrm{\&phi;}}_1\right)\&CenterDot;(1-H\left({\mathrm{\&phi;}}_2\right))\&CenterDot;I\left(Y\right)\mathrm{dY}}{{\&Integral;}_{{\mathrm{\&Omega;}}_Y}B(X,Y)\&CenterDot;H\left({\mathrm{\&phi;}}_1\right)\&CenterDot;(1-H\left({\mathrm{\&phi;}}_2\right))\mathrm{dY}}$

$c_{-+}=\frac{{\&Integral;}_{{\mathrm{\&Omega;}}_Y}B(X,Y)\&CenterDot;(1-H\left({\mathrm{\&phi;}}_1\right))\&CenterDot;H\left({\mathrm{\&phi;}}_2\right)\&CenterDot;I\left(Y\right)\mathrm{dY}}{{\&Integral;}_{{\mathrm{\&Omega;}}_Y}B(X,Y)\&CenterDot;(1-H\left({\mathrm{\&phi;}}_1\right))\&CenterDot;H\left({\mathrm{\&phi;}}_2\right)\mathrm{dY}}$

$c_{--}=\frac{{\&Integral;}_{{\mathrm{\&Omega;}}_Y}B(X,Y)\&CenterDot;(1-H\left({\mathrm{\&phi;}}_1\right))\&CenterDot;(1-H\left({\mathrm{\&phi;}}_2\right))\&CenterDot;I\left(Y\right)\mathrm{dY}}{{\&Integral;}_{{\mathrm{\&Omega;}}_Y}B(X,Y)\&CenterDot;(1-H\left({\mathrm{\&phi;}}_1\right))\&CenterDot;(1-H\left({\mathrm{\&phi;}}_2\right))\mathrm{dY}}$

I(Y)表示待分割数据；

步骤4-4、设置权重因子μ₁，μ₂，ν₁，ν₂，λ₁，λ₂，λ₃，λ₄，定义多相水平集能量泛函：

E_total＝μ₁·E₁+μ₂·E₂+ν₁·E₃+ν₂·E₄+λ₁·E₅+λ₂·E₆+λ₃·E₇+λ₄·E₈；

步骤4-5、以局部化区域函数B(X，Y)为局部区域掩膜，计算以X为中心、r为半径的圆圈内的能量分量，

水平集方程φ₁的能量泛函表示为：

$E_{\mathrm{total}}\left({\mathrm{\&phi;}}_1\right)={\&Integral;}_{{\mathrm{\&Omega;}}_X}\mathrm{\&delta;}\left({\mathrm{\&phi;}}_1\right){\&Integral;}_{{\mathrm{\&Omega;}}_Y}B(X,Y)[{\mathrm{\&lambda;}}_1{\left|\right|u-c_{++}\left|\right|}^{2}H\left({\mathrm{\&phi;}}_1\right)H\left({\mathrm{\&phi;}}_2\right)+{\mathrm{\&lambda;}}_2{\left|\right|u-c_{+-}\left|\right|}^{2}H\left({\mathrm{\&phi;}}_1\right)(1-H\left({\mathrm{\&phi;}}_2\right))$

$+{\mathrm{\&lambda;}}_3{\left|\right|u-c_{-+}\left|\right|}^2(1-H\left({\mathrm{\&phi;}}_1\right))H\left({\mathrm{\&phi;}}_2\right)+{\mathrm{\&lambda;}}_4{\left|\right|u-c_{--}\left|\right|}^2(1-H\left({\mathrm{\&phi;}}_1\right))(1-H\left({\mathrm{\&phi;}}_2\right))]\mathrm{dYdX}$

$+{\mathrm{\&mu;}}_1\&CenterDot;\mathrm{length}\left(C_1\right)+{\mathrm{\&mu;}}_2\&CenterDot;\mathrm{length}\left(C_2\right)+{\mathrm{\&nu;}}_1\&CenterDot;\mathrm{area}\left(\mathrm{inside}\left(C_1\right)\right)+{\mathrm{\&nu;}}_2\&CenterDot;\mathrm{area}\left(\mathrm{inside}\left(C_2\right)\right)$

水平集方程φ₂的能量泛函表示为：

$E_{\mathrm{total}}\left({\mathrm{\&phi;}}_2\right)={\&Integral;}_{{\mathrm{\&Omega;}}_X}\mathrm{\&delta;}\left({\mathrm{\&phi;}}_2\right){\&Integral;}_{{\mathrm{\&Omega;}}_Y}B(X,Y)[{\mathrm{\&lambda;}}_1{\left|\right|u-c_{++}\left|\right|}^{2}H\left({\mathrm{\&phi;}}_1\right)H\left({\mathrm{\&phi;}}_2\right)+{\mathrm{\&lambda;}}_2{\left|\right|u-c_{+-}\left|\right|}^{2}H\left({\mathrm{\&phi;}}_1\right)(1-H\left({\mathrm{\&phi;}}_2\right))$

$+{\mathrm{\&lambda;}}_3{\left|\right|u-c_{-+}\left|\right|}^2(1-H\left({\mathrm{\&phi;}}_1\right))H\left({\mathrm{\&phi;}}_2\right)+{\mathrm{\&lambda;}}_4{\left|\right|u-c_{--}\left|\right|}^2(1-H\left({\mathrm{\&phi;}}_1\right))(1-H\left({\mathrm{\&phi;}}_2\right))]\mathrm{dYdX}$

$+{\mathrm{\&mu;}}_1\&CenterDot;\mathrm{length}\left(C_1\right)+{\mathrm{\&mu;}}_2\&CenterDot;\mathrm{length}\left(C_2\right)+{\mathrm{\&nu;}}_1\&CenterDot;\mathrm{area}\left(\mathrm{inside}\left(C_1\right)\right)+{\mathrm{\&nu;}}_2\&CenterDot;\mathrm{area}\left(\mathrm{inside}\left(C_2\right)\right);$

步骤4-6、根据梯度下降法，分别计算步骤4-5获得的能量泛函E_total(φ₁)和E_total(φ₂)的演化方程：

$\frac{{\&PartialD;\mathrm{\&phi;}}_1}{\&PartialD;t}=\mathrm{\&delta;}\left({\mathrm{\&phi;}}_1\right)\{{\mathrm{\&mu;}}_1\&CenterDot;\mathrm{div}\left(\frac{{\&dtri;\mathrm{\&phi;}}_1}{|\&dtri;\mathrm{\&phi;}|}\right)-{\mathrm{\&nu;}}_1+{\&Integral;}_{{\mathrm{\&Omega;}}_Y}B(X,Y)\mathrm{\&delta;}\left({\mathrm{\&phi;}}_1\right)[({\mathrm{\&lambda;}}_3{\left|\right|u-c_{-+}\left|\right|}^2-{\mathrm{\&lambda;}}_1{\left|\right|u-c_{++}\left|\right|}^2)H\left({\mathrm{\&phi;}}_2\right)$

$-({\mathrm{\&lambda;}}_2{\left|\right|u-c_{+-}\left|\right|}^2-{\mathrm{\&lambda;}}_4{\left|\right|u-c_{--}\left|\right|}^2)(1-H\left({\mathrm{\&phi;}}_2\right))\left]\mathrm{dY}\right\}$

$\frac{{\&PartialD;\mathrm{\&phi;}}_2}{\&PartialD;t}=\mathrm{\&delta;}\left({\mathrm{\&phi;}}_2\right)\{{\mathrm{\&mu;}}_2\&CenterDot;\mathrm{div}\left(\frac{{\&dtri;\mathrm{\&phi;}}_2}{\left|{\&dtri;\mathrm{\&phi;}}_2\right|}\right)-{\mathrm{\&nu;}}_2+{\&Integral;}_{{\mathrm{\&Omega;}}_Y}B(X,Y)\mathrm{\&delta;}\left({\mathrm{\&phi;}}_2\right)[({\mathrm{\&lambda;}}_2{\left|\right|u-c_{+-}\left|\right|}^2-{\mathrm{\&lambda;}}_1{\left|\right|u-c_{++}\left|\right|}^2)H\left({\mathrm{\&phi;}}_1\right)$

$-({\mathrm{\&lambda;}}_3{\left|\right|u-c_{-+}\left|\right|}^2-{\mathrm{\&lambda;}}_4{\left|\right|u-c_{--}\left|\right|}^2)(1-H\left({\mathrm{\&phi;}}_1\right))\left]\mathrm{dY}\right\}$

步骤4-7、判断能量泛函是否达到了最小值，如果是则完成将待分割区域进行多相水平集分割的过程，将待分割区域分割成四个子区域，如果否，则返回步骤4-2，迭代更新水平集方程φ₁和φ₂。

5.根据权利要求4所述的一种基于水平集的LOD2建筑物模型构建方法，其特征在于：所述步骤4-7中判断能量泛函是否达到了最小值的依据为：通过设定阈值或设定迭代次数。

6.根据权利要求1所述的一种基于水平集的LOD2建筑物模型构建方法，其特征在于：所述步骤5中获取建筑物每个基元的角点的图像坐标具体方法为：

步骤5-1、令水平集方程φ₁＝0，φ₂＝0，获得建筑物顶面几何结构边缘的二值图；

步骤5-2、对步骤5-1获得的建筑物顶面几何结构边缘二值图应用Harris算子检测建筑物结构每个基元的角点的图像坐标，Harris算子为R＝det(C)-ktr²(C)，其中

表示图像点z的灰度在u的偏导，

表示图像点z的灰度在v的偏导，I_uv(z)表示图像点z的灰度在u和v的二阶混合偏导。本实施方式中，

和I_uv(z)分别是图像点z的灰度在u和v的偏导以及二阶混合偏导，一般取0.04～0.06。

Images