一种无人机多重叠遥感影像的建筑物轮廓线提取方法

摘要

一种无人机多重叠遥感影像的建筑物轮廓线提取方法，包括利用空三结合密集匹配的方法生成三维点云，并对点云进行滤波处理，从其中检测出建筑物。对检测的建筑删除墙面后，从建筑物顶面信息提取建筑物粗轮廓。建筑物粗轮廓作为缓冲区叠加拼接影像上，利用建筑物粗轮廓作为形状先验信息，在缓冲区内用水平集算法进行演化，最后得到建筑物精确轮廓。本发明充分利用了多重叠影像生成的点云三维信息，同时结合高分辨率遥感影像的高精度几何信息，不但显著提高了建筑物轮廓提取的精度，而且降低了方法的复杂度。

说明书

技术领域

本发明涉及遥感影像应用技术领域，尤其是涉及一种无人机遥感影像建筑物 轮廓线提取方法。

背景技术

建筑物是城市中一种重要的地理空间要素，它在城市规划与管理、城市发展 与变化以及灾害检测与评估等应用领域具有重要的意义。建筑物轮廓线提取是城 市基础地理信息系统的建立和更新的一个重要步骤。无人机作为一种新型遥感监 测平台，飞行操作智能化程度高，可按预定航线自主飞行、摄像，实时提供遥感 监测数据和低空视频监控，具有机动性强、便捷、成本低等特点，其所获取的高 分辨率重叠的遥感数据具有抗干扰能力强，成像范围大等特点，使之成为建筑物 轮廓线提取有效的数据来源之一。

高分辨率遥感影像的中包含了大量丰富的信息，建筑物轮廓提取往往受到各 种其它的地物干扰，比如建筑物和非建筑物区分，建筑物周围树木的遮挡，道路 边线的影响等等。因此，仅仅靠单一的影像进行建筑物轮廓提取，技术难度很大。 建筑物轮廓提取不仅需要依靠遥感二维信息的提取与分析，而且还需要结合建筑 物三维信息，所以二维和三维信息互为融合和补充将更加有利于遥感影像中建筑 物轮廓的提取。目前利用高分辨率遥感影像进行建筑物轮廓提取的典型方法包括 以下几种：1)基于单一的高分辨率遥感影像建筑物轮廓线提取。虽然高分辨率的 遥感影像具有清晰的建筑物轮廓信息，但是人造的建筑物和非建筑物难以区分开 来，另外建筑物周围的树木遮挡也对建筑物的轮廓产生一定的干扰，因此这类方 法具有一定的局限性。2)基于阴影辅助下的建筑物轮廓线提取。虽然在阴影辅助 下进行建筑物轮廓提取间接利用了建筑物的高度信息，但是阴影的提取不具有一 定的普适性，而且利用阴影求得建筑物高度的需要相关的参数较多，因此此类方 法很难满足实际的需要。3)基于Lidar和遥感影像的建筑物轮廓线提取。虽然 这类方法既利用了Lidar的三维信息，又融合了影像的高精度几何轮廓信息，通 过两种数据优劣的互为补充来提取建筑物轮廓信息。但是这类方法存在是Lidar 和遥感影像的高精度配准困难，而且Lidar数据获取的成本也较为昂贵。4)基于 立体航空影像的建筑物轮廓线提取。虽然这类方法利用立体匹配获得了三维信 息，同时利用了影像的高精度的二维信息，通过两类信息的互补进行建筑物轮廓 信息提取。但是这类方法的问题是立体像对幅面较小，对于提取大范围的城区建 筑物轮廓有一定的影响。因此需要迫切寻找一种数据易获取、提取方法自动化程 度高、提取结果相对精确高且符合实际生产需要的方法。

发明内容

本发明充分利用了多重叠影像生成的点云三维信息，同时结合遥感影像自身 的高精度信息，显著提高了建筑物轮廓提取的精度。

本发明的技术方案为一种无人机多重叠遥感影像的建筑物轮廓线提取方法， 包括以下步骤：

步骤一，利用空三对无人机遥感影像进行平差，同时利用GPU加速后的PMVS 算法对影像密集匹配，最后得到精度高的密集彩色点云；

步骤二，对平差后的无人机遥感影像进行拼接；

步骤三，对彩色点云进行滤波；先利用改进的形态学滤波算法进行地面和非地面 分离，然后利用颜色不变量对地面点中的植被滤除，最后利用高程和面积作为阈 值滤除非建筑物；

步骤四，利用区域增长算法检测点云中的建筑物；

步骤五，删除建筑物的墙面，通过对顶面边界拟合最后得到建筑物的粗轮廓信息；

步骤六，利用步骤三得到的建筑物粗轮廓作为叠加在拼接影像上，形成建筑物轮 廓提取的缓冲区；

步骤七，同时利用建筑物粗轮廓的形状作为先验信息，在缓冲区内用水平集算法 演化出建筑物精确轮廓。

所述步骤一包括以下步骤：

(1.1)利用先验信息对多视重叠无人机遥感影像进行预处理：

(1.2)在步骤(1.1)的基础上进行空三摄影测量，利用空三勾网，求出每张影 像的外方位元素，并进行光束法的整体平差；

(1.3)根据影像分组，在步骤(1..2)的基础上利用现有技术中GPU加速的PMVS 算法进行快速的密集匹配，生成密集的三维点云，所重建的点云作为三维高程数 据。

所述步骤二包括以下步骤：

(2.1)特征提取：利用SIFT进行影像的特征提取；

(2.2)影像配准：先进行粗配准，利用k-d树搜索匹配的特征点；然后进行 精配准，粗配准往往出现错误的匹配点，因此利用RANSAC算法剔除错 误的匹配点；通过影像的配准，得到影像之间的变换矩阵；

(2.3)影像的拼接：通过(2.2)得到的变换矩阵进行影像的拼接；

(2.4)影像的融合：拼接后，利用双线性插值算法进行影像的融合。

所述步骤三包括以下步骤：

(3.1)利用改进的形态学滤波对点云的地面点和非地面点进行分离；

(3.2)对步骤(3.1)中的地面点中利用颜色不变量对植被进行滤波；

(3.3)基于建筑物的特点，利用阈值过滤掉非建筑物目标点。

所述步骤(3.1)包括以下过程：

首先取任意一个点和它的邻域点组成一个固定大小的窗口，通过形态学的开 运算检测出窗口内的最低点，如果窗口内的点的高程值与最低点高程之差在阈值 范围内，表示该点为地面点，以此取出点云中所有点进行滤波；

其次根据y＝2×w^k+1获得下次滤波所需的窗口大小，且该窗口的大小小于预 设的滤波窗口最大值，再进行一次形态学滤波；最后当窗口大于预设窗口，结束 滤波；其中，k为迭代次数，w前一次滤波窗口的大小。

所述步骤(3.2)包括以下过程：

由于由影像密集匹配生成的点云具有颜色信息，因此利用颜色不变量理论对 绿色植被进行过滤；设点云中每个点的坐标为(x,y,z)颜色三个通道为(R,G,B)， 颜色不变量的对于植被的阈值为T_g,利用绿色和蓝色通道定义的颜色不变量公式 为：

${\psi }_g(x,y,z)=\frac{4}{\pi }\times \arctan \left(\frac{I_g(x,y,z)-I_b(x,y,z)}{I_g(x,y,z)+I_b(x,y,z)}\right)---\left(1\right)$

其中，I_g(x,y,z)、I_b(x,y,z)表示点云在(x,y,z)点的绿色和蓝色分量值；

ψ_g(x,y,z)表示在(x,y,z)点的颜色不变量值；当ψ_g＜T_g时，表示该点为植被 点；当ψ_g＞T_g时，表示该点为非植被点。

所述步骤四，利用基于平面拟合区域增长法分割步骤三中的点云，得到每个 建筑物的点云区域V_Bi，具体包括以下步骤：

(4.1)把点云中的点分为两类：如果点的邻域中有一个点缺省，这类点属于边 界点；否则这类点属于内部点；

(4.2)设某个内部点p₀(x₀,y₀,z₀),p₀的八邻域点为：{p₁,p_2,…,p₈}，对这9个点 利用最小二乘法拟合出平面方程(2)，具体实现如下：

z＝ax+by+c   (2)

对于点p₀和邻域内的8个点：(x_i,y_i,z_i),i＝0,1,2,…,7；拟合计算式(2)平

面方程，则使(3)式最小；

$S=\underset{i=0}{\overset{7}{\Sigma }}{(\mathrm{ax}+\mathrm{by}+c-z)}^2---\left(3\right)$

要使得S最小，应满足：

由此可得下列正规方程组：

$\left(\begin{array}{ccc}\Sigma {x_i}^2& \Sigma x_i{y}_i& \Sigma x_i\\ \Sigma x_i{y}_i& \Sigma {y_i}^2& \Sigma y_i\\ \Sigma x_i& \Sigma y_i& n\end{array}\right)\left(\begin{array}{c}a\\ b\\ c\end{array}\right)=\left(\begin{array}{c}\Sigma x_i{z}_i\\ \Sigma y_i{z}_i\\ \Sigma z_i\end{array}\right)i=\mathrm{0,1},...7;n=8---\left(4\right)$

(4.3)根据公式(4)求得任何内部点在8邻域内的方差和SSD：

$\mathrm{SSD}=\underset{\left(p_k\right)\in M}{\Sigma }{(z_k-h_k)}^2---\left(5\right)$

其中M是p₀和它的八邻域点的集合，h_k和z_k分别是观测高程值和拟合平面 高程值；

(4.4)对点云中每个点按照SSD值进行升序排序，取最小的SSD值的点作为种 子点；

(4.5)在种子点邻域内求每个点到种子点所在平面的距离h，如果h<h_T,h_T为 高程阈值，则合并到同一区域，并对该区域进行标记；否则，这为非面片点；

(4.6)当邻域点全部遍历完，从新的SSD值中找未处理的点作为新的种子点， 反复(4.5)的操作，直到所有的点遍历结束。

所述步骤五，对步骤四中的每个建筑物的点云区域V_Bi删除建筑物的墙面点， 通过求取顶面边界点最后得到建筑物的粗轮廓信息，具体包括以下步骤：

(5.1)根据步骤四分割的面片，如果面片的法向量平行于地面，说明是墙面面 片；如果面片的法向量垂直地面，说明是屋顶面片；因此根据法向量删除墙面面 片；

(5.2)根据(5.1)中得到的顶面点云，利用alpha-shape算法得到点云的边界；具 体方法是：

a)从顶面点云中取出任意一点p1，从剩余点中搜索距离小于等于2α的点集 合p2，а为圆的半径，设p2＝{p21,p22,p23,…,p2n}；

b)从p2中任取出一点p2i，利用公式(6)(7)求出过p1和p2i点的圆心p0；已 知两点(x1,y1)、(x2,y2)和圆的半径α，求该圆的圆心(x0,y0)方程如下：

$\left(\begin{array}{c}{(x1-x0)}^2+{(y1-y0)}^2={\alpha }^2\\ {(x2-x0)}^2+{(y2-y0)}^2={\alpha }^2\end{array}\right)---\left(6\right)$

直接求取此方程比较困难，因此利用测绘学中的距离交汇算法得：

$\left(\begin{array}{c}x0=x1+\frac{x2-x1}{2}+H(y2-y1)\\ y0=y1+\frac{y2-y1}{2}+H(x2-x1)\end{array}\right)---\left(7\right)$

其中，S²＝(x1-x2)²+(y1-y2)²；

c)从点集合p2中求出所有点到p0的距离l，如果l>α，那么p1和p2i是 边界点；如果l<α，转入下一步d)；

d)对p2中其它点重复a)b)c)三步，直到p2中所有点全部判断结束。

所述步骤六，将步骤五中的得到每个建筑物的顶面点云轮廓点叠加在影像 上，得到影像的上的筑物缓冲区Ω_Bi；具体过程如下：

设投影矩阵P_3×4表示单视图几何的已知的内外方位元素，计算P_3×4矩阵：

$P=\left(\begin{array}{ccc}-f& 0& x_0\\ 0& f& y_0\\ 0& 0& 1\end{array}\right)R^T\left(\begin{array}{cccc}1& 0& 0& -X_s\\ 0& 1& 0& -Y_s\\ 0& 0& 1& -Z_s\end{array}\right)---\left(8\right)$

其中，P即投影矩阵P_3×4，f为影像的焦距，x₀和y₀为光轴距离光心在水平 方向和垂直方向的偏心距；X_s、Y_s、Z_s为相机中心在世界坐标系中的坐标，R^T表 示3×3的旋转矩阵；

投影公式为：

$\left(\begin{array}{c}\mathrm{xz}\\ \mathrm{yz}\\ z\end{array}\right)=P\left(\begin{array}{c}X\\ Y\\ Z\\ 1\end{array}\right)---\left(9\right)$

利用式(8)和式(9)计算点云上的点到影像上的投影；x,y,z表示物点在相 机坐标系下的坐标，X,Y,Z表示物点在世界坐标系下的坐标。

所述步骤七，把步骤六叠加到影像上的点进行边界连接，得到的建筑物粗轮 廓作为缓冲区叠加在拼接影像上，同时利用建筑物粗轮廓的形状作为先验信息， 在缓冲区内用水平集算法演化出建筑物精确轮廓；具体步骤如下：

(7.1)对叠加在影像上的点进行边界跟踪，得到封闭的轮廓；对该轮廓进行扩 大，得到影像上建筑物轮廓提取的缓冲区Ω_Bi；

(7.2)根据(7.1)得到影像区域的轮廓作为初始水平集，根据(7.1)得到影像区 域的轮廓形状作为先验信息在缓冲区内部的遥感影像中进行局部水平集演化，得 到遥感影像的建筑物的精细轮廓；具体实施如下：

首先，预设影像上每个建筑物的缓冲区为Ω_Bi，(7.1)中获取轮廓的形状为 Ф_Pi，Ф_Pi作为先验形状，Ф_Si作为目标分割轮廓，基于形状约束的能量泛函定 义如下：

E_to＝E(c₁,c₂,φ_Si)+E_sh(φ_Si,φ_P1) (10)

目标的初始轮廓和先验形状均利用运动目标区域进行表达：

$E_{\mathrm{sh}}({\phi }_{\mathrm{si}},{\phi }_{\mathrm{Pi}})={\int }_{{\Omega }_{\mathrm{Bi}}}{({\phi }_{\mathrm{si}}(x,y)-{\phi }_{\mathrm{Pi}}(x,y))}^2\mathrm{dxdy}---\left(11\right)$

现有技术基于水平集C-V模型的能量函数为：

$\left(\begin{array}{c}E(\phi ,c_1,c_2)={\lambda }_1{\int }_{\Omega }{(u_0(x,y)-u_1)}^{2}H\left(\phi \right)\mathrm{dxdy}+\\ {\lambda }_2{\int }_{\Omega }{(u_0(x,y)-u_2)}^2(1-H\left(\phi \right))\mathrm{dxdy}+\\ \mu {\int }_{\Omega }H\left(\phi \right)\mathrm{dxdy}+v{\int }_{\Omega }\delta \left(\phi \right)|▿\phi |\mathrm{dxdy}\end{array}\right)---\left(12\right)$

式中，H(φ)是Heaviside函数，其形式为：$H\left(\phi \right)=\left(\begin{array}{cc}1,& \phi \ge 0\\ 0,& \phi <0\end{array}\right),$

δ(φ)为Dirac函数，其形式为u₀(x,y)是待处理影像区域某一点 灰度值，为当前点梯度的模，系数λ₁，λ₂＞0，μ，v≥0为固定参数，一 般取λ₁＝λ₂＝1，μ＝0，v＝1，(12)式对应的偏微分方程：

$\frac{\partial \phi }{\partial t}=\delta \left(\phi \right)·[\mu ▿·\frac{▿\phi }{|▿\phi |}-v-{\lambda }_1{[u_0(x,y)-c_1]}^2+{\lambda }_2{[u_0(x,y)-c_2]}^2]---\left(13\right)$

其中，参数c₁、c₂根据下式得到，

$\left(\begin{array}{cc}{c}_1=\frac{{\int }_{\Omega }{u}_0(x,y)H\left(\phi \right)\mathrm{dxdy}}{{\int }_{\Omega }H\left(\phi \right)\mathrm{dxdy}}& c_2=\frac{{\int }_{\Omega }{u}_0(x,y)[1-H\left(\phi \right)]\mathrm{dxdy}}{{\int }_{\Omega }[1-H\left(\phi \right)]\mathrm{dxdy}}\end{array}\right)$

根据式(13)，采用C-V水平集演化方法分割提取建筑物区域，作为影像的

建筑物分割对像v_iSeg，C-V水平集为现有技术。

与现有技术相比，本发明具有以下优点和有益效果：

(1)虽然是单一的数据源，但是可以生成了两类数据信息的，达到了两类数据 源提取建筑物轮廓的效果。降低了方法的复杂度，也节约了生产成本。

(2)利用影像生成三维点云的三维信息获取粗建筑物的轮廓，同时也利用高分 辨影像的高精度轮廓信息，两者互为补充和融合，为建筑物轮廓的提取提高的 自动化程度和精度。

(3)利用三维点云信息提取的建筑物粗轮廓作为水平集演化的先验形状信息和 建筑物轮廓提取的缓冲区，保证了利用水平集演化建筑物轮廓的速度和精度。

附图说明

图1为本发明的流程图；

图2为步骤一中A点的影像在航带图中的分组图。

具体实施方式

本发明提出了一种无人机遥感影像的建筑物轮廓线提取方法，该方法先利用 摄影测量中的空三算法结合计算机视觉中的多视几何立体重建快速生成带有地 理坐标的三维点云，再通过三维点云分割提取出建筑物的粗略轮廓信息，将建筑 物的粗略轮廓信息建立成缓冲区叠加在高分辨率影像上，然后利用粗略轮廓作为 先验形状信息，在缓冲区内利用水平集进行演化迭代，得到建筑物精确的几何轮 廓。由于三维点云由影像密集匹配生成，因此点云和影像之间不存在的配准困难。 以下结合附图和实施例详细说明本发明技术方案，流程图如图1所示，实施例的 技术方案流程包括以下步骤：

步骤一，利用空三对无人机遥感影像进行平差，同时利用GPU加速后的PMVS 算法对影像密集匹配，最后得到精度高且带有地理坐标的密集彩色点云。具体实 施如下：

(1)利用先验信息对多视重叠无人机遥感影像进行预处理：

无人机航拍的相邻影像之间有一定的重叠度。由于航拍的数据量非常大，直 接进行三维重建，一方面无法得到较好的重建效果，另外一方面会使得重建的计 算量大，重建时间较长。因此，利用已有POS信息和航带先验信息对影像进行分 组。由于实施中无人机影像的航向重叠度是80％，旁向重叠度是35％，那么对于 某张影像应该和同一航带的连续4张影像以及航带间的连续两张影像分为一组。 如图2所示A点的影像在航带图中的分组，黑色矩形虚线框住部分为与A影像分 在同一组的影像。

(2)在步骤(1)的基础上进行空三摄影测量，利用空三勾网，求出每张影像的 外方位元素，并进行光束法的整体平差。本步骤实现可采用现有技术，本发明不 予赘述。

(3)根据影像分组，在步骤(2)的基础上利用现有技术中GPU加速的PMVS算 法进行快速的密集匹配，生成密集的三维点云，所重建的点云作为三维高程数据。 步骤二，对平差后的无人机遥感影像进行拼接。具体流程如下：

(1)特征提取：利用SIFT进行影像的特征提取。

(2)影像配准：先进行粗配准，利用k-d树搜索匹配的特征点；然后进行精配准， 粗配准往往出现错误的匹配点，因此利用RANSAC算法剔除错误的匹配点。通 过影像的配准，可以得到影像之间的变换矩阵。

(3)影像的拼接：通过(2)得到的变换矩阵进行影像的拼接。

(4)影像的融合：拼接后，利用双线性插值算法进行影像的融合。

步骤三，对彩色点云进行滤波。具体实施如下：

(1)利用改进的形态学滤波对点云的地面点和非地面点进行分离。其运算过程是： 首先取任意一个点和它的邻域点组成一个固定大小的窗口，通过形态学的开 运算检测出窗口内的最低点，如果窗口内的点的高程值与最低点高程之差在 阈值范围内，表示该点为地面点，以此取出点云中所有点进行滤波；其次根 据y＝2×w^k+1(其中k为迭代次数，w前一次滤波窗口的大小)获得下次滤 波所需的窗口大小，且该窗口的大小小于预设的滤波窗口最大值，在进行一 次形态学滤波；最后当窗口大于预设窗口，结束滤波。

(2)对步骤(1)中的地面点中利用颜色不变量对植被进行滤波。由于由影像密集 匹配生成的点云具有颜色信息，因此可以利用颜色不变量理论对绿色植被进 行过滤。设点云中每个点的坐标为(x,y,z)颜色三个通道为(R,G,B)，颜色 不变量的对于植被的阈值为T_g,利用绿色和蓝色通道定义的颜色不变量公式 为：

${\psi }_g(x,y,z)=\frac{4}{\pi }\times \arctan \left(\frac{I_g(x,y,z)-I_b(x,y,z)}{I_g(x,y,z)+I_b(x,y,z)}\right)---\left(1\right)$

其中，I_g(x,y,z)、I_b(x,y,z)表示点云在(x,y,z)点的绿色和蓝色分量值。 ψ_g(x,y,z)表示在(x,y,z)点的颜色不变量值。当ψ_g＜T_g时，表示该点为植被 点；当ψ_g＞T_g时，表示该点为非植被点。

(3)基于建筑物的特点，例如高度不会低于2米，面积不应该小于35平米。利用 阈值过滤掉非建筑物目标点。

步骤四，利用基于平面拟合区域增长法分割步骤三中的点云，得到每个建筑物的 点云区域V_Bi。具体流程如下：

(1)把点云中的点分为两类：如果点的邻域中有一个点缺省，这类点属于边界点； 否则这类点属于内部点。

(2)设某个内部点p₀(x₀,y₀,z₀),p₀的八邻域点为：{p₁,p_2,…,p₈}。对这9个点利用 最小二乘法拟合出平面方程(2)，具体实现如下：

z＝ax+by+c   (2)

其中，a、b、c表示式(2)平面方程的参数；对于点p₀和邻域内的8个点：

(x_i,y_i,z_i),i＝0,1,2,…,7；拟合计算式(2)平面方程，则使(3)式最小。

$S=\underset{i=0}{\overset{7}{\Sigma }}{(\mathrm{ax}+\mathrm{by}+c-z)}^2---\left(3\right)$

要使得S最小，应满足：

由此可得下列正规方程组：

$\left(\begin{array}{ccc}\Sigma {x_i}^2& \Sigma x_i{y}_i& \Sigma x_i\\ \Sigma x_i{y}_i& \Sigma {y_i}^2& \Sigma y_i\\ \Sigma x_i& \Sigma y_i& n\end{array}\right)\left(\begin{array}{c}a\\ b\\ c\end{array}\right)=\left(\begin{array}{c}\Sigma x_i{z}_i\\ \Sigma y_i{z}_i\\ \Sigma z_i\end{array}\right)i=\mathrm{0,1},...7;n=8---\left(4\right)$

(3)根据公式(4)求得任何内部点在8邻域内的方差和SSD：

$\mathrm{SSD}=\underset{\left(p_k\right)\in M}{\Sigma }{(z_k-h_k)}^2---\left(5\right)$

其中M是p₀和它的八邻域点的集合，h_k和z_k分别是观测高程值和拟合平面 高程值。

(4)对点云中每个点按照SSD值进行升序排序，取最小的SSD值的点作为种子点。

(5)在种子点邻域内求每个点到种子点所在平面的距离h，如果h<h_T,合并到同 一区域，并对该区域进行标记。否则，这为非面片点。本领域技术人员可自 行预设相应阈值。

(6)当邻域点全部遍历完，从新的SSD值中找未处理的点作为新的种子点，反复 (5)的操作，直到所有的点遍历结束。

步骤五，对步骤四中的每个建筑物的点云区域V_Bi删除建筑物的墙面点，通过求 取顶面边界点最后得到建筑物的粗轮廓信息。具体实施过程如下：

(1)根据步骤四分割的面片，如果面片的法向量平行于地面，说明是墙面面片； 如果面片的法向量垂直地面，说明是屋顶面片。因此根据法向量删除墙面面 片。

(2)根据(1)中得到的顶面点云，利用alpha-shape算法得到点云的边界。具体方法 是：a)从顶面点云中取出任意一点p1，从剩余点中搜索距离小于等于2α(α 为圆的半径)的点集合p2,设p2＝{p21,p22,p23,…,p2n}；b)从p2中任取出一点 p2i，利用公式(6)(7)求出过p1和p2i点的圆心p0；已知两点(x1,y1)、(x2,y2) 和圆的半径α，求该圆的圆心(x0,y0)方程如下：

$\left(\begin{array}{c}{(x1-x0)}^2+{(y1-y0)}^2={\alpha }^2\\ {(x2-x0)}^2+{(y2-y0)}^2={\alpha }^2\end{array}\right)---\left(6\right)$

直接求取此方程比较困难，因此利用测绘学中的距离交汇算法得：

$\left(\begin{array}{c}x0=x1+\frac{x2-x1}{2}+H(y2-y1)\\ y0=y1+\frac{y2-y1}{2}+H(x2-x1)\end{array}\right)---\left(7\right)$

其中，S²＝(x1-x2)²+(y1-y2)²。

c)从点集合p2中求出所有点到p0的距离l，如果l>α，那么p1和p2i是边 界点；如果l<α，转入下一步d)；d)对p2中其它点重复a)b)c)三步，直到p2 中所有点全部判断结束。

步骤六，将步骤五中的得到每个建筑物的顶面点云轮廓点叠加在影像上，得到影 像的上的筑物缓冲区Ω_Bi。由于步骤五中的建筑物顶面边界点是由影像多视立体 匹配重建获得，因此把这些点重新投影回到影像上只需要利用二维到三维的投影 矩阵作逆运算即可。具体流程如下：

设投影矩阵P_3×4表示单视图几何的已知的内外方位元素，计算P_3×4矩阵：

$P=\left(\begin{array}{ccc}-f& 0& x_0\\ 0& f& y_0\\ 0& 0& 1\end{array}\right)R^T\left(\begin{array}{cccc}1& 0& 0& -X_s\\ 0& 1& 0& -Y_s\\ 0& 0& 1& -Z_s\end{array}\right)---\left(8\right)$

其中，P即投影矩阵P_3×4，f为影像的焦距，x₀和y₀为光轴距离光心在水平 方向和垂直方向的偏心距。X_s、Y_s、Z_s为相机中心在世界坐标系中的坐标，R^T表 示3×3的旋转矩阵。

投影公式为：

$\left(\begin{array}{c}\mathrm{xz}\\ \mathrm{yz}\\ z\end{array}\right)=P\left(\begin{array}{c}X\\ Y\\ Z\\ 1\end{array}\right)---\left(9\right)$

利用式(8)和式(9)计算点云上的点到影像上的投影。x,y,z表示物点在相机坐 标系下的坐标，X,Y,Z表示物点在世界坐标系下的坐标。

步骤七，把步骤六叠加到影像上的点进行边界连接，得到的建筑物粗轮廓作为缓 冲区叠加在拼接影像上，同时利用建筑物粗轮廓的形状作为先验信息，在缓冲区 内用水平集算法演化出建筑物精确轮廓。具体过程如下：

(1)对叠加在影像上的点进行边界跟踪，得到封闭的轮廓。对该轮廓进行扩大， 得到影像上建筑物轮廓提取的缓冲区Ω_Bi。

(2)根据(1)得到影像区域的轮廓作为初始水平集，根据(1)得到影像区域的轮廓 形状作为先验信息在缓冲区内部的遥感影像中进行局部水平集演化，得到遥 感影像的建筑物的精细轮廓。具体实施如下：

首先，预设影像上每个建筑物的缓冲区为Ω_Bi，(1)中获取轮廓的形状为Ф_Pi， Ф_Pi作为先验形状，Ф_Si作为目标分割轮廓，基于形状约束的能量泛函定义如 下：

E_to＝E(c₁,c₂,φ_Si)+E_sh(φ_Si,φ_P1) (10)

目标的初始轮廓和先验形状均可以利用运动目标区域进行表达：

$E_{\mathrm{sh}}({\phi }_{\mathrm{si}},{\phi }_{\mathrm{Pi}})={\int }_{{\Omega }_{\mathrm{Bi}}}{({\phi }_{\mathrm{si}}(x,y)-{\phi }_{\mathrm{Pi}}(x,y))}^2\mathrm{dxdy}---\left(11\right)$

现有技术基于水平集C-V模型的能量函数为：

$\left(\begin{array}{c}E(\phi ,c_1,c_2)={\lambda }_1{\int }_{\Omega }{(u_0(x,y)-u_1)}^{2}H\left(\phi \right)\mathrm{dxdy}+\\ {\lambda }_2{\int }_{\Omega }{(u_0(x,y)-u_2)}^2(1-H\left(\phi \right))\mathrm{dxdy}+\\ \mu {\int }_{\Omega }H\left(\phi \right)\mathrm{dxdy}+v{\int }_{\Omega }\delta \left(\phi \right)|▿\phi |\mathrm{dxdy}\end{array}\right)---\left(12\right)$

式中，H（φ）是Heaviside函数，其形式为：$H\left(\phi \right)=\left(\begin{array}{cc}1,& \phi \ge 0\\ 0,& \phi <0\end{array}\right),$

δ（φ）为Dirac函数，其形式为u₀（x,y）是待处理影像区域某一点 灰度值，为当前点梯度的模，系数λ₁，λ₂＞0，μ，v≥0为固定参数，一 般建议取λ₁＝λ₂＝1，μ＝0，v＝1，(12)式对应的偏微分方程：

$\frac{\partial \phi }{\partial t}=\delta \left(\phi \right)·[\mu ▿·\frac{▿\phi }{|▿\phi |}-v-{\lambda }_1{[u_0(x,y)-c_1]}^2+{\lambda }_2{[u_0(x,y)-c_2]}^2]---\left(13\right)$

其中，参数c₁、c₂根据下式得到，

$\left(\begin{array}{cc}{c}_1=\frac{{\int }_{\Omega }{u}_0(x,y)H\left(\phi \right)\mathrm{dxdy}}{{\int }_{\Omega }H\left(\phi \right)\mathrm{dxdy}}& c_2=\frac{{\int }_{\Omega }{u}_0(x,y)[1-H\left(\phi \right)]\mathrm{dxdy}}{{\int }_{\Omega }[1-H\left(\phi \right)]\mathrm{dxdy}}\end{array}\right)$

根据式(13)，采用C-V水平集演化方法分割提取建筑物区域，作为影像的

建筑物分割对像v_iSeg，C-V水平集为现有技术，本发明不予赘述。

为距离d_Buffer1在灾后遥感影像上相应位置处为v_i建立局部缓冲区v_iBuf，缓冲区 v_iBuf轮廓以内的影像区域为当前待处理区域Ω。d_Buffer1可由本领域技术人员根据 情况自行预先设定。