CN104484668B - 一种无人机多重叠遥感影像的建筑物轮廓线提取方法 - Google Patents
一种无人机多重叠遥感影像的建筑物轮廓线提取方法 Download PDFInfo
- Publication number
- CN104484668B CN104484668B CN201510025503.9A CN201510025503A CN104484668B CN 104484668 B CN104484668 B CN 104484668B CN 201510025503 A CN201510025503 A CN 201510025503A CN 104484668 B CN104484668 B CN 104484668B
- Authority
- CN
- China
- Prior art keywords
- mrow
- msub
- mtd
- mtr
- point
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
- 238000000034 method Methods 0.000 title claims abstract description 53
- 238000004422 calculation algorithm Methods 0.000 claims abstract description 29
- 238000000605 extraction Methods 0.000 claims abstract description 21
- 238000012545 processing Methods 0.000 claims abstract description 5
- 238000001914 filtration Methods 0.000 claims description 32
- 239000011159 matrix material Substances 0.000 claims description 19
- 230000000877 morphologic effect Effects 0.000 claims description 11
- 239000013598 vector Substances 0.000 claims description 9
- 230000003287 optical effect Effects 0.000 claims description 6
- 230000009466 transformation Effects 0.000 claims description 6
- 239000000284 extract Substances 0.000 claims description 5
- 230000011218 segmentation Effects 0.000 claims description 5
- 238000004364 calculation method Methods 0.000 claims description 4
- 230000001133 acceleration Effects 0.000 claims description 3
- 230000001174 ascending effect Effects 0.000 claims description 3
- 238000013507 mapping Methods 0.000 claims description 3
- 238000007781 pre-processing Methods 0.000 claims description 3
- 238000012163 sequencing technique Methods 0.000 claims description 3
- 230000004927 fusion Effects 0.000 claims description 2
- 238000009795 derivation Methods 0.000 claims 1
- 238000001514 detection method Methods 0.000 abstract description 2
- 230000003139 buffering effect Effects 0.000 abstract 2
- 238000012217 deletion Methods 0.000 abstract 1
- 230000037430 deletion Effects 0.000 abstract 1
- 238000012544 monitoring process Methods 0.000 description 3
- 230000000694 effects Effects 0.000 description 2
- 238000004519 manufacturing process Methods 0.000 description 2
- 238000004458 analytical method Methods 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 238000011161 development Methods 0.000 description 1
- 230000018109 developmental process Effects 0.000 description 1
- 238000010586 diagram Methods 0.000 description 1
- 238000011156 evaluation Methods 0.000 description 1
- 238000007726 management method Methods 0.000 description 1
Landscapes
- Measurement Of Radiation (AREA)
- Image Processing (AREA)
- Image Analysis (AREA)
Abstract
一种无人机多重叠遥感影像的建筑物轮廓线提取方法，包括利用空三结合密集匹配的方法生成三维点云，并对点云进行滤波处理，从其中检测出建筑物。对检测的建筑删除墙面后，从建筑物顶面信息提取建筑物粗轮廓。建筑物粗轮廓作为缓冲区叠加拼接影像上，利用建筑物粗轮廓作为形状先验信息，在缓冲区内用水平集算法进行演化，最后得到建筑物精确轮廓。本发明充分利用了多重叠影像生成的点云三维信息，同时结合高分辨率遥感影像的高精度几何信息，不但显著提高了建筑物轮廓提取的精度，而且降低了方法的复杂度。
Description
技术领域
本发明涉及遥感影像应用技术领域，尤其是涉及一种无人机遥感影像建筑物轮廓线提取方法。
背景技术
建筑物是城市中一种重要的地理空间要素，它在城市规划与管理、城市发展与变化以及灾害检测与评估等应用领域具有重要的意义。建筑物轮廓线提取是城市基础地理信息***的建立和更新的一个重要步骤。无人机作为一种新型遥感监测平台，飞行操作智能化程度高，可按预定航线自主飞行、摄像，实时提供遥感监测数据和低空视频监控，具有机动性强、便捷、成本低等特点，其所获取的高分辨率重叠的遥感数据具有抗干扰能力强，成像范围大等特点，使之成为建筑物轮廓线提取有效的数据来源之一。
高分辨率遥感影像的中包含了大量丰富的信息，建筑物轮廓提取往往受到各种其它的地物干扰，比如建筑物和非建筑物区分，建筑物周围树木的遮挡，道路边线的影响等等。因此，仅仅靠单一的影像进行建筑物轮廓提取，技术难度很大。建筑物轮廓提取不仅需要依靠遥感二维信息的提取与分析，而且还需要结合建筑物三维信息，所以二维和三维信息互为融合和补充将更加有利于遥感影像中建筑物轮廓的提取。目前利用高分辨率遥感影像进行建筑物轮廓提取的典型方法包括以下几种：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×wk+1获得下次滤波所需的窗口大小，且该窗口的大小小于预设的滤波窗口最大值，再进行一次形态学滤波；最后当窗口大于预设窗口，结束滤波；其中，k为迭代次数，w前一次滤波窗口的大小。
所述步骤(3.2)包括以下过程：
由于由影像密集匹配生成的点云具有颜色信息，因此利用颜色不变量理论对绿色植被进行过滤；设点云中每个点的坐标为(x,y,z)颜色三个通道为(R,G,B)，颜色不变量的对于植被的阈值为Tg,利用绿色和蓝色通道定义的颜色不变量公式为：
其中，Ig(x,y,z)、Ib(x,y,z)表示点云在(x,y,z)点的绿色和蓝色分量值；
ψg(x,y,z)表示在(x,y,z)点的颜色不变量值；当ψg＜Tg时，表示该点为植被点；当ψg＞Tg时，表示该点为非植被点。
所述步骤四，利用基于平面拟合区域增长法分割步骤三中的点云，得到每个建筑物的点云区域VBi，具体包括以下步骤：
(4.1)把点云中的点分为两类：如果点的邻域中有一个点缺省，这类点属于边界点；否则这类点属于内部点；
(4.2)设某个内部点p0(x0,y0,z0),p0的八邻域点为：{p1,p2,…,p8}，对这9个点利用最小二乘法拟合出平面方程(2)，具体实现如下：
z＝ax+by+c (2)
对于点p0和邻域内的8个点：(xi,yi,zi),i＝0,1,2,…,7；拟合计算式(2)平
面方程，则使(3)式最小；
要使得S最小，应满足：
由此可得下列正规方程组：
(4.3)根据公式(4)求得任何内部点在8邻域内的方差和SSD：
其中M是p0和它的八邻域点的集合，hk和zk分别是观测高程值和拟合平面高程值；
(4.4)对点云中每个点按照SSD值进行升序排序，取最小的SSD值的点作为种子点；
(4.5)在种子点邻域内求每个点到种子点所在平面的距离h，如果h<hT,hT为高程阈值，则合并到同一区域，并对该区域进行标记；否则，这为非面片点；
(4.6)当邻域点全部遍历完，从新的SSD值中找未处理的点作为新的种子点，反复(4.5)的操作，直到所有的点遍历结束。
所述步骤五，对步骤四中的每个建筑物的点云区域VBi删除建筑物的墙面点，通过求取顶面边界点最后得到建筑物的粗轮廓信息，具体包括以下步骤：
(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)方程如下：
直接求取此方程比较困难，因此利用测绘学中的距离交汇算法得：
其中，S2＝(x1-x2)2+(y1-y2)2；
c)从点集合p2中求出所有点到p0的距离l，如果l>α，那么p1和p2i是边界点；如果l<α，转入下一步d)；
d)对p2中其它点重复a)b)c)三步，直到p2中所有点全部判断结束。
所述步骤六，将步骤五中的得到每个建筑物的顶面点云轮廓点叠加在影像上，得到影像的上的筑物缓冲区ΩBi；具体过程如下：
设投影矩阵P3×4表示单视图几何的已知的内外方位元素，计算P3×4矩阵：
其中，P即投影矩阵P3×4，f为影像的焦距，x0和y0为光轴距离光心在水平方向和垂直方向的偏心距；Xs、Ys、Zs为相机中心在世界坐标系中的坐标，RT表示3×3的旋转矩阵；
投影公式为：
利用式(8)和式(9)计算点云上的点到影像上的投影；x,y,z表示物点在相机坐标系下的坐标，X,Y,Z表示物点在世界坐标系下的坐标。
所述步骤七，把步骤六叠加到影像上的点进行边界连接，得到的建筑物粗轮廓作为缓冲区叠加在拼接影像上，同时利用建筑物粗轮廓的形状作为先验信息，在缓冲区内用水平集算法演化出建筑物精确轮廓；具体步骤如下：
(7.1)对叠加在影像上的点进行边界跟踪，得到封闭的轮廓；对该轮廓进行扩大，得到影像上建筑物轮廓提取的缓冲区ΩBi；
(7.2)根据(7.1)得到影像区域的轮廓作为初始水平集，根据(7.1)得到影像区域的轮廓形状作为先验信息在缓冲区内部的遥感影像中进行局部水平集演化，得到遥感影像的建筑物的精细轮廓；具体实施如下：
首先，预设影像上每个建筑物的缓冲区为ΩBi，(7.1)中获取轮廓的形状为ФPi，ФPi作为先验形状，ФSi作为目标分割轮廓，基于形状约束的能量泛函定义如下：
Eto＝E(c1,c2,φSi)+Esh(φSi,φP1) (10)
目标的初始轮廓和先验形状均利用运动目标区域进行表达：
现有技术基于水平集C-V模型的能量函数为：
式中，H(φ)是Heaviside函数，其形式为：
δ(φ)为Dirac函数，其形式为u0(x,y)是待处理影像区域某一点灰度值，为当前点梯度的模，系数λ1，λ2＞0，μ，v≥0为固定参数，一般取λ1＝λ2＝1，μ＝0，v＝1，(12)式对应的偏微分方程：
其中，参数c1、c2根据下式得到，
根据式(13)，采用C-V水平集演化方法分割提取建筑物区域，作为影像的
建筑物分割对像viSeg，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×wk+1(其中k为迭代次数，w前一次滤波窗口的大小)获得下次滤波所需的窗口大小，且该窗口的大小小于预设的滤波窗口最大值，在进行一次形态学滤波；最后当窗口大于预设窗口，结束滤波。
(2)对步骤(1)中的地面点中利用颜色不变量对植被进行滤波。由于由影像密集匹配生成的点云具有颜色信息，因此可以利用颜色不变量理论对绿色植被进行过滤。设点云中每个点的坐标为(x,y,z)颜色三个通道为(R,G,B)，颜色不变量的对于植被的阈值为Tg,利用绿色和蓝色通道定义的颜色不变量公式为：
其中，Ig(x,y,z)、Ib(x,y,z)表示点云在(x,y,z)点的绿色和蓝色分量值。ψg(x,y,z)表示在(x,y,z)点的颜色不变量值。当ψg＜Tg时，表示该点为植被点；当ψg＞Tg时，表示该点为非植被点。
(3)基于建筑物的特点，例如高度不会低于2米，面积不应该小于35平米。利用阈值过滤掉非建筑物目标点。
步骤四，利用基于平面拟合区域增长法分割步骤三中的点云，得到每个建筑物的点云区域VBi。具体流程如下：
(1)把点云中的点分为两类：如果点的邻域中有一个点缺省，这类点属于边界点；否则这类点属于内部点。
(2)设某个内部点p0(x0,y0,z0),p0的八邻域点为：{p1,p2,…,p8}。对这9个点利用最小二乘法拟合出平面方程(2)，具体实现如下：
z＝ax+by+c (2)
其中，a、b、c表示式(2)平面方程的参数；对于点p0和邻域内的8个点：
(xi,yi,zi),i＝0,1,2,…,7；拟合计算式(2)平面方程，则使(3)式最小。
要使得S最小，应满足：
由此可得下列正规方程组：
(3)根据公式(4)求得任何内部点在8邻域内的方差和SSD：
其中M是p0和它的八邻域点的集合，hk和zk分别是观测高程值和拟合平面高程值。
(4)对点云中每个点按照SSD值进行升序排序，取最小的SSD值的点作为种子点。
(5)在种子点邻域内求每个点到种子点所在平面的距离h，如果h<hT,合并到同一区域，并对该区域进行标记。否则，这为非面片点。本领域技术人员可自行预设相应阈值。
(6)当邻域点全部遍历完，从新的SSD值中找未处理的点作为新的种子点，反复(5)的操作，直到所有的点遍历结束。
步骤五，对步骤四中的每个建筑物的点云区域VBi删除建筑物的墙面点，通过求取顶面边界点最后得到建筑物的粗轮廓信息。具体实施过程如下：
(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)方程如下：
直接求取此方程比较困难，因此利用测绘学中的距离交汇算法得：
其中，S2＝(x1-x2)2+(y1-y2)2。
c)从点集合p2中求出所有点到p0的距离l，如果l>α，那么p1和p2i是边界点；如果l<α，转入下一步d)；d)对p2中其它点重复a)b)c)三步，直到p2中所有点全部判断结束。
步骤六，将步骤五中的得到每个建筑物的顶面点云轮廓点叠加在影像上，得到影像的上的筑物缓冲区ΩBi。由于步骤五中的建筑物顶面边界点是由影像多视立体匹配重建获得，因此把这些点重新投影回到影像上只需要利用二维到三维的投影矩阵作逆运算即可。具体流程如下：
设投影矩阵P3×4表示单视图几何的已知的内外方位元素，计算P3×4矩阵：
其中，P即投影矩阵P3×4，f为影像的焦距，x0和y0为光轴距离光心在水平方向和垂直方向的偏心距。Xs、Ys、Zs为相机中心在世界坐标系中的坐标，RT表示3×3的旋转矩阵。
投影公式为：
利用式(8)和式(9)计算点云上的点到影像上的投影。x,y,z表示物点在相机坐标系下的坐标，X,Y,Z表示物点在世界坐标系下的坐标。
步骤七，把步骤六叠加到影像上的点进行边界连接，得到的建筑物粗轮廓作为缓冲区叠加在拼接影像上，同时利用建筑物粗轮廓的形状作为先验信息，在缓冲区内用水平集算法演化出建筑物精确轮廓。具体过程如下：
(1)对叠加在影像上的点进行边界跟踪，得到封闭的轮廓。对该轮廓进行扩大，得到影像上建筑物轮廓提取的缓冲区ΩBi。
(2)根据(1)得到影像区域的轮廓作为初始水平集，根据(1)得到影像区域的轮廓形状作为先验信息在缓冲区内部的遥感影像中进行局部水平集演化，得到遥感影像的建筑物的精细轮廓。具体实施如下：
首先，预设影像上每个建筑物的缓冲区为ΩBi，(1)中获取轮廓的形状为ФPi，ФPi作为先验形状，ФSi作为目标分割轮廓，基于形状约束的能量泛函定义如下：
Eto＝E(c1,c2,φSi)+Esh(φSi,φP1) (10)
目标的初始轮廓和先验形状均可以利用运动目标区域进行表达：
现有技术基于水平集C-V模型的能量函数为：
式中，H（φ）是Heaviside函数，其形式为：
δ（φ）为Dirac函数，其形式为u0（x,y）是待处理影像区域某一点灰度值，为当前点梯度的模，系数λ1，λ2＞0，μ，v≥0为固定参数，一般建议取λ1＝λ2＝1，μ＝0，v＝1，(12)式对应的偏微分方程：
其中，参数c1、c2根据下式得到，
根据式(13)，采用C-V水平集演化方法分割提取建筑物区域，作为影像的
建筑物分割对像viSeg，C-V水平集为现有技术，本发明不予赘述。
为距离dBuffer1在灾后遥感影像上相应位置处为vi建立局部缓冲区viBuf，缓冲区viBuf轮廓以内的影像区域为当前待处理区域Ω。dBuffer1可由本领域技术人员根据情况自行预先设定。
Claims (3)
1.一种无人机多重叠遥感影像的建筑物轮廓线提取方法，其特征在于，包括以下步骤：
步骤一，利用空三对无人机遥感影像进行平差，同时利用GPU加速后的PMVS算法对影像密集匹配，最后得到精度高的密集彩色点云；
步骤二，对平差后的无人机遥感影像进行拼接；
步骤三，对彩色点云进行滤波；
先利用改进的形态学滤波算法进行地面和非地面分离，然后利用颜色不变量对地面点中的植被滤除，最后利用高程和面积作为阈值滤除非建筑物；
步骤四，利用区域增长算法检测点云中的建筑物；
步骤五，删除建筑物的墙面，通过对顶面边界拟合最后得到建筑物的粗轮廓信息；
步骤六，利用步骤五得到的建筑物粗轮廓叠加在拼接影像上，形成建筑物轮廓提取的缓冲区；
步骤七，同时利用建筑物粗轮廓的形状作为先验信息，在缓冲区内用水平集算法演化出建筑物精确轮廓；
所述步骤三包括以下步骤：
(3.1)利用改进的形态学滤波对点云的地面点和非地面点进行分离；
(3.2)对步骤(3.1)中的地面点利用颜色不变量对植被进行滤波；
(3.3)基于建筑物的特点，利用阈值过滤掉非建筑物目标点；
所述步骤(3.1)包括以下过程：
首先取任意一个点和它的邻域点组成一个固定大小的窗口，通过形态学的开运算检测出窗口内的最低点，如果窗口内的点的高程值与最低点高程之差在阈值范围内，表示该点为地面点，以此取出点云中所有点进行滤波；
其次根据y＝2×wk+1获得下次滤波所需的窗口大小，且该窗口的大小小于预设的滤波窗口最大值，再进行一次形态学滤波；最后当窗口大于预设窗口，结束滤波；
其中，k为迭代次数，w为前一次滤波窗口的大小；
所述步骤(3.2)包括以下过程：
由于由影像密集匹配生成的点云具有颜色信息，因此利用颜色不变量理论对绿色植被进行过滤；设点云中每个点的坐标为(x,y,z)，颜色三个通道为(R,G,B)，颜色不变量对于植被的阈值为Tg,利用绿色和蓝色通道定义的颜色不变量公式为：
<mrow>
<msub>
<mi>&psi;</mi>
<mi>g</mi>
</msub>
<mrow>
<mo>(</mo>
<mi>x</mi>
<mo>,</mo>
<mi>y</mi>
<mo>,</mo>
<mi>z</mi>
<mo>)</mo>
</mrow>
<mo>=</mo>
<mfrac>
<mn>4</mn>
<mi>&pi;</mi>
</mfrac>
<mo>&times;</mo>
<mi>a</mi>
<mi>r</mi>
<mi>c</mi>
<mi>t</mi>
<mi>a</mi>
<mi>n</mi>
<mrow>
<mo>(</mo>
<mfrac>
<mrow>
<msub>
<mi>I</mi>
<mi>g</mi>
</msub>
<mrow>
<mo>(</mo>
<mi>x</mi>
<mo>,</mo>
<mi>y</mi>
<mo>,</mo>
<mi>z</mi>
<mo>)</mo>
</mrow>
<mo>-</mo>
<msub>
<mi>I</mi>
<mi>b</mi>
</msub>
<mrow>
<mo>(</mo>
<mi>x</mi>
<mo>,</mo>
<mi>y</mi>
<mo>,</mo>
<mi>z</mi>
<mo>)</mo>
</mrow>
</mrow>
<mrow>
<msub>
<mi>I</mi>
<mi>g</mi>
</msub>
<mrow>
<mo>(</mo>
<mi>x</mi>
<mo>,</mo>
<mi>y</mi>
<mo>,</mo>
<mi>z</mi>
<mo>)</mo>
</mrow>
<mo>+</mo>
<msub>
<mi>I</mi>
<mi>b</mi>
</msub>
<mrow>
<mo>(</mo>
<mi>x</mi>
<mo>,</mo>
<mi>y</mi>
<mo>,</mo>
<mi>z</mi>
<mo>)</mo>
</mrow>
</mrow>
</mfrac>
<mo>)</mo>
</mrow>
<mo>-</mo>
<mo>-</mo>
<mo>-</mo>
<mrow>
<mo>(</mo>
<mn>1</mn>
<mo>)</mo>
</mrow>
</mrow>
其中，Ig(x,y,z)、Ib(x,y,z)表示点云在(x,y,z)点的绿色和蓝色分量值；ψg(x,y,z)表示在(x,y,z)点的颜色不变量值；当ψg＜Tg时，表示该点为植被点；当ψg＞Tg时，表示该点为非植被点；
所述步骤四，利用基于平面拟合区域增长法分割步骤三中的点云，得到每个建筑物的点云区域VBi，具体包括以下步骤：
(4.1)把点云中的点分为两类：如果点的邻域中有一个点缺省，这类点属于边界点；否则这类点属于内部点；
(4.2)设某个内部点p0(x0,y0,z0),p0的八邻域点为：{p1,p2,…,p8}，对这9个点利用最小二乘法拟合出平面方程(2)，具体实现如下：
z＝ax+by+c (2)
对于点p0和邻域内的8个点：(xi,yi,zi),i＝0,1,2,…,7；拟合计算式(2)平面方程，则使(3)式最小；
<mrow>
<mi>S</mi>
<mo>=</mo>
<munderover>
<mo>&Sigma;</mo>
<mrow>
<mi>i</mi>
<mo>=</mo>
<mn>0</mn>
</mrow>
<mn>7</mn>
</munderover>
<msup>
<mrow>
<mo>(</mo>
<mi>a</mi>
<mi>x</mi>
<mo>+</mo>
<mi>b</mi>
<mi>y</mi>
<mo>+</mo>
<mi>c</mi>
<mo>-</mo>
<mi>z</mi>
<mo>)</mo>
</mrow>
<mn>2</mn>
</msup>
<mo>-</mo>
<mo>-</mo>
<mo>-</mo>
<mrow>
<mo>(</mo>
<mn>3</mn>
<mo>)</mo>
</mrow>
</mrow>
要使得S最小，应满足：其中，a、b、c表示式(2)平面方程的参数，和分别表示S对a、b和c求导；
由此可得下列正规方程组：
<mrow>
<mfenced open = "|" close = "|">
<mtable>
<mtr>
<mtd>
<mrow>
<mo>&Sigma;</mo>
<msup>
<msub>
<mi>x</mi>
<mi>i</mi>
</msub>
<mn>2</mn>
</msup>
</mrow>
</mtd>
<mtd>
<mrow>
<mo>&Sigma;</mo>
<msub>
<mi>x</mi>
<mi>i</mi>
</msub>
<msub>
<mi>y</mi>
<mi>i</mi>
</msub>
</mrow>
</mtd>
<mtd>
<mrow>
<mo>&Sigma;</mo>
<msub>
<mi>x</mi>
<mi>i</mi>
</msub>
</mrow>
</mtd>
</mtr>
<mtr>
<mtd>
<mrow>
<mo>&Sigma;</mo>
<msub>
<mi>x</mi>
<mi>i</mi>
</msub>
<msub>
<mi>y</mi>
<mi>i</mi>
</msub>
</mrow>
</mtd>
<mtd>
<mrow>
<mo>&Sigma;</mo>
<msup>
<msub>
<mi>y</mi>
<mi>i</mi>
</msub>
<mn>2</mn>
</msup>
</mrow>
</mtd>
<mtd>
<mrow>
<mo>&Sigma;</mo>
<msub>
<mi>y</mi>
<mi>i</mi>
</msub>
</mrow>
</mtd>
</mtr>
<mtr>
<mtd>
<mrow>
<mo>&Sigma;</mo>
<msub>
<mi>x</mi>
<mi>i</mi>
</msub>
</mrow>
</mtd>
<mtd>
<mrow>
<mo>&Sigma;</mo>
<msub>
<mi>y</mi>
<mi>i</mi>
</msub>
</mrow>
</mtd>
<mtd>
<mi>n</mi>
</mtd>
</mtr>
</mtable>
</mfenced>
<mfenced open = "(" close = ")">
<mtable>
<mtr>
<mtd>
<mi>a</mi>
</mtd>
</mtr>
<mtr>
<mtd>
<mi>b</mi>
</mtd>
</mtr>
<mtr>
<mtd>
<mi>c</mi>
</mtd>
</mtr>
</mtable>
</mfenced>
<mo>=</mo>
<mfenced open = "(" close = ")">
<mtable>
<mtr>
<mtd>
<mrow>
<mo>&Sigma;</mo>
<msub>
<mi>x</mi>
<mi>i</mi>
</msub>
<msub>
<mi>z</mi>
<mi>i</mi>
</msub>
</mrow>
</mtd>
</mtr>
<mtr>
<mtd>
<mrow>
<mo>&Sigma;</mo>
<msub>
<mi>y</mi>
<mi>i</mi>
</msub>
<msub>
<mi>z</mi>
<mi>i</mi>
</msub>
</mrow>
</mtd>
</mtr>
<mtr>
<mtd>
<mrow>
<mo>&Sigma;</mo>
<msub>
<mi>z</mi>
<mi>i</mi>
</msub>
</mrow>
</mtd>
</mtr>
</mtable>
</mfenced>
<mo>,</mo>
<mi>i</mi>
<mo>=</mo>
<mn>0</mn>
<mo>,</mo>
<mn>1</mn>
<mo>,</mo>
<mn>...7</mn>
<mo>;</mo>
<mi>n</mi>
<mo>=</mo>
<mn>8</mn>
<mo>-</mo>
<mo>-</mo>
<mo>-</mo>
<mrow>
<mo>(</mo>
<mn>4</mn>
<mo>)</mo>
</mrow>
</mrow>
(4.3)根据公式(4)求得任何内部点在8邻域内的方差和SSD：
<mrow>
<mi>S</mi>
<mi>S</mi>
<mi>D</mi>
<mo>=</mo>
<munder>
<mo>&Sigma;</mo>
<mrow>
<mo>(</mo>
<msub>
<mi>p</mi>
<mi>k</mi>
</msub>
<mo>)</mo>
<mo>&Element;</mo>
<mi>M</mi>
</mrow>
</munder>
<msup>
<mrow>
<mo>(</mo>
<msub>
<mi>z</mi>
<mi>k</mi>
</msub>
<mo>-</mo>
<msub>
<mi>h</mi>
<mi>k</mi>
</msub>
<mo>)</mo>
</mrow>
<mn>2</mn>
</msup>
<mo>-</mo>
<mo>-</mo>
<mo>-</mo>
<mrow>
<mo>(</mo>
<mn>5</mn>
<mo>)</mo>
</mrow>
</mrow>
其中M是p0和它的八邻域点的集合，hk和zk分别是观测高程值和拟合平面高程值；
(4.4)对点云中每个点按照SSD值进行升序排序，取最小的SSD值的点作为种子点；
(4.5)在种子点邻域内求每个点到种子点所在平面的距离h，如果h<hT,hT为高程阈值，则合并到同一区域，并对该区域进行标记；否则，这为非面片点；
(4.6)当邻域点全部遍历完，从新的SSD值中找未处理的点作为新的种子点，反复(4.5)的操作，直到所有的点遍历结束；
所述步骤五，对步骤四中的每个建筑物的点云区域VBi删除建筑物的墙面点，通过求取顶面边界点最后得到建筑物的粗轮廓信息，具体包括以下步骤：
(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)方程如下：
<mrow>
<mfenced open = "{" close = "">
<mtable>
<mtr>
<mtd>
<mrow>
<msup>
<mrow>
<mo>(</mo>
<mi>x</mi>
<mn>1</mn>
<mo>-</mo>
<mi>x</mi>
<mn>0</mn>
<mo>)</mo>
</mrow>
<mn>2</mn>
</msup>
<mo>+</mo>
<msup>
<mrow>
<mo>(</mo>
<mi>y</mi>
<mn>1</mn>
<mo>-</mo>
<mi>y</mi>
<mn>0</mn>
<mo>)</mo>
</mrow>
<mn>2</mn>
</msup>
<mo>=</mo>
<msup>
<mi>&alpha;</mi>
<mn>2</mn>
</msup>
</mrow>
</mtd>
</mtr>
<mtr>
<mtd>
<mrow>
<msup>
<mrow>
<mo>(</mo>
<mi>x</mi>
<mn>2</mn>
<mo>-</mo>
<mi>x</mi>
<mn>0</mn>
<mo>)</mo>
</mrow>
<mn>2</mn>
</msup>
<mo>+</mo>
<msup>
<mrow>
<mo>(</mo>
<mi>y</mi>
<mn>2</mn>
<mo>-</mo>
<mi>y</mi>
<mn>0</mn>
<mo>)</mo>
</mrow>
<mn>2</mn>
</msup>
<mo>=</mo>
<msup>
<mi>&alpha;</mi>
<mn>2</mn>
</msup>
</mrow>
</mtd>
</mtr>
</mtable>
</mfenced>
<mo>-</mo>
<mo>-</mo>
<mo>-</mo>
<mrow>
<mo>(</mo>
<mn>6</mn>
<mo>)</mo>
</mrow>
</mrow>
直接求取此方程比较困难，因此利用测绘学中的距离交汇算法得：
<mrow>
<mfenced open = "{" close = "">
<mtable>
<mtr>
<mtd>
<mrow>
<mi>x</mi>
<mn>0</mn>
<mo>=</mo>
<mi>x</mi>
<mn>1</mn>
<mo>+</mo>
<mfrac>
<mrow>
<mi>x</mi>
<mn>2</mn>
<mo>-</mo>
<mi>x</mi>
<mn>1</mn>
</mrow>
<mn>2</mn>
</mfrac>
<mo>+</mo>
<mi>H</mi>
<mrow>
<mo>(</mo>
<mi>y</mi>
<mn>2</mn>
<mo>-</mo>
<mi>y</mi>
<mn>1</mn>
<mo>)</mo>
</mrow>
</mrow>
</mtd>
</mtr>
<mtr>
<mtd>
<mrow>
<mi>y</mi>
<mn>0</mn>
<mo>=</mo>
<mi>y</mi>
<mn>1</mn>
<mo>+</mo>
<mfrac>
<mrow>
<mi>y</mi>
<mn>2</mn>
<mo>-</mo>
<mi>y</mi>
<mn>1</mn>
</mrow>
<mn>2</mn>
</mfrac>
<mo>+</mo>
<mi>H</mi>
<mrow>
<mo>(</mo>
<mi>x</mi>
<mn>2</mn>
<mo>-</mo>
<mi>x</mi>
<mn>1</mn>
<mo>)</mo>
</mrow>
</mrow>
</mtd>
</mtr>
</mtable>
</mfenced>
<mo>-</mo>
<mo>-</mo>
<mo>-</mo>
<mrow>
<mo>(</mo>
<mn>7</mn>
<mo>)</mo>
</mrow>
</mrow>
2
其中，S2＝(x1-x2)2+(y1-y2)2；
c)从点集合p2中求出所有点到p0的距离l，如果l>α，那么p1和p2i是边界点；如果l<α，转入下一步d)；
d)对p2中其它点重复a)b)c)三步，直到p2中所有点全部判断结束；
所述步骤六，将步骤五中得到的每个建筑物的顶面点云轮廓点叠加在影像上，得到影像上的建筑物轮廓提取的缓冲区ΩBi；具体过程如下：
设投影矩阵P3×4表示单视图几何的已知的内外方位元素，计算P3×4矩阵：
<mrow>
<mi>P</mi>
<mo>=</mo>
<mfenced open = "[" close = "]">
<mtable>
<mtr>
<mtd>
<mrow>
<mo>-</mo>
<mi>f</mi>
</mrow>
</mtd>
<mtd>
<mn>0</mn>
</mtd>
<mtd>
<msub>
<mi>x</mi>
<mn>0</mn>
</msub>
</mtd>
</mtr>
<mtr>
<mtd>
<mn>0</mn>
</mtd>
<mtd>
<mi>f</mi>
</mtd>
<mtd>
<msub>
<mi>y</mi>
<mn>0</mn>
</msub>
</mtd>
</mtr>
<mtr>
<mtd>
<mn>0</mn>
</mtd>
<mtd>
<mn>0</mn>
</mtd>
<mtd>
<mn>1</mn>
</mtd>
</mtr>
</mtable>
</mfenced>
<msup>
<mi>R</mi>
<mi>T</mi>
</msup>
<mfenced open = "[" close = "]">
<mtable>
<mtr>
<mtd>
<mn>1</mn>
</mtd>
<mtd>
<mn>0</mn>
</mtd>
<mtd>
<mn>0</mn>
</mtd>
<mtd>
<mrow>
<mo>-</mo>
<msub>
<mi>X</mi>
<mi>s</mi>
</msub>
</mrow>
</mtd>
</mtr>
<mtr>
<mtd>
<mn>0</mn>
</mtd>
<mtd>
<mn>1</mn>
</mtd>
<mtd>
<mn>0</mn>
</mtd>
<mtd>
<mrow>
<mo>-</mo>
<msub>
<mi>Y</mi>
<mi>s</mi>
</msub>
</mrow>
</mtd>
</mtr>
<mtr>
<mtd>
<mn>0</mn>
</mtd>
<mtd>
<mn>0</mn>
</mtd>
<mtd>
<mn>1</mn>
</mtd>
<mtd>
<mrow>
<mo>-</mo>
<msub>
<mi>Z</mi>
<mi>s</mi>
</msub>
</mrow>
</mtd>
</mtr>
</mtable>
</mfenced>
<mo>-</mo>
<mo>-</mo>
<mo>-</mo>
<mrow>
<mo>(</mo>
<mn>8</mn>
<mo>)</mo>
</mrow>
</mrow>
其中，P即投影矩阵P3×4，f为影像的焦距，x0和y0为光轴距离光心在水平方向和垂直方向的偏心距；Xs、Ys、Zs为相机中心在世界坐标系中的坐标，RT表示3×3的旋转矩阵；
投影公式为：
<mrow>
<mfenced open = "[" close = "]">
<mtable>
<mtr>
<mtd>
<mrow>
<mi>x</mi>
<mi>z</mi>
</mrow>
</mtd>
</mtr>
<mtr>
<mtd>
<mrow>
<mi>y</mi>
<mi>z</mi>
</mrow>
</mtd>
</mtr>
<mtr>
<mtd>
<mi>z</mi>
</mtd>
</mtr>
</mtable>
</mfenced>
<mo>=</mo>
<mi>P</mi>
<mfenced open = "[" close = "]">
<mtable>
<mtr>
<mtd>
<mi>X</mi>
</mtd>
</mtr>
<mtr>
<mtd>
<mi>Y</mi>
</mtd>
</mtr>
<mtr>
<mtd>
<mi>Z</mi>
</mtd>
</mtr>
<mtr>
<mtd>
<mn>1</mn>
</mtd>
</mtr>
</mtable>
</mfenced>
<mo>-</mo>
<mo>-</mo>
<mo>-</mo>
<mrow>
<mo>(</mo>
<mn>9</mn>
<mo>)</mo>
</mrow>
</mrow>
利用式(8)和式(9)计算点云上的点到影像上的投影；x,y,z表示物点在相机坐标系下的坐标，X,Y,Z表示物点在世界坐标系下的坐标；
所述步骤七，把步骤六叠加到影像上的点进行边界连接，得到的建筑物粗轮廓叠加在拼接影像上作为缓冲区，同时利用建筑物粗轮廓的形状作为先验信息，在缓冲区内用水平集算法演化出建筑物精确轮廓；具体步骤如下：
(7.1)对叠加在影像上的点进行边界跟踪，得到封闭的轮廓；对该轮廓进行扩大，得到影像上建筑物轮廓提取的缓冲区ΩBi；
(7.2)根据(7.1)得到影像区域的轮廓作为初始水平集，根据(7.1)得到影像区域的轮廓形状作为先验信息在缓冲区内部的遥感影像中进行局部水平集演化，得到遥感影像的建筑物的精细轮廓；具体实施如下：
首先，预设影像上每个建筑物的缓冲区为ΩBi，(7.1)中获取轮廓的形状为ФPi，ФPi作为先验形状，ФSi作为目标分割轮廓，基于形状约束的能量泛函定义如下：
Eto＝E(c1,c2,φSi)+Esh(φSi,φP1)(10)
目标的初始轮廓和先验形状均利用运动目标区域进行表达：
<mrow>
<msub>
<mi>E</mi>
<mrow>
<mi>s</mi>
<mi>h</mi>
</mrow>
</msub>
<mrow>
<mo>(</mo>
<msub>
<mi>&phi;</mi>
<mrow>
<mi>s</mi>
<mi>i</mi>
</mrow>
</msub>
<mo>,</mo>
<msub>
<mi>&phi;</mi>
<mrow>
<mi>P</mi>
<mi>i</mi>
</mrow>
</msub>
<mo>)</mo>
</mrow>
<mo>=</mo>
<msub>
<mo>&Integral;</mo>
<msub>
<mi>&Omega;</mi>
<mrow>
<mi>B</mi>
<mi>i</mi>
</mrow>
</msub>
</msub>
<msup>
<mrow>
<mo>(</mo>
<msub>
<mi>&phi;</mi>
<mrow>
<mi>s</mi>
<mi>i</mi>
</mrow>
</msub>
<mo>(</mo>
<mrow>
<mi>x</mi>
<mo>,</mo>
<mi>y</mi>
</mrow>
<mo>)</mo>
<mo>-</mo>
<msub>
<mi>&phi;</mi>
<mrow>
<mi>P</mi>
<mi>i</mi>
</mrow>
</msub>
<mo>(</mo>
<mrow>
<mi>x</mi>
<mo>,</mo>
<mi>y</mi>
</mrow>
<mo>)</mo>
<mo>)</mo>
</mrow>
<mn>2</mn>
</msup>
<mi>d</mi>
<mi>x</mi>
<mi>d</mi>
<mi>y</mi>
<mo>-</mo>
<mo>-</mo>
<mo>-</mo>
<mrow>
<mo>(</mo>
<mn>11</mn>
<mo>)</mo>
</mrow>
</mrow>
基于水平集C-V模型的能量函数为：
<mrow>
<mtable>
<mtr>
<mtd>
<mrow>
<mi>E</mi>
<mrow>
<mo>(</mo>
<mi>&phi;</mi>
<mo>,</mo>
<msub>
<mi>c</mi>
<mn>1</mn>
</msub>
<mo>,</mo>
<msub>
<mi>c</mi>
<mn>2</mn>
</msub>
<mo>)</mo>
</mrow>
<mo>=</mo>
<msub>
<mi>&lambda;</mi>
<mn>1</mn>
</msub>
<msub>
<mo>&Integral;</mo>
<mi>&Omega;</mi>
</msub>
<msup>
<mrow>
<mo>(</mo>
<msub>
<mi>u</mi>
<mn>0</mn>
</msub>
<mo>(</mo>
<mrow>
<mi>x</mi>
<mo>,</mo>
<mi>y</mi>
</mrow>
<mo>)</mo>
<mo>-</mo>
<msub>
<mi>u</mi>
<mn>1</mn>
</msub>
<mo>)</mo>
</mrow>
<mn>2</mn>
</msup>
<mi>H</mi>
<mrow>
<mo>(</mo>
<mi>&phi;</mi>
<mo>)</mo>
</mrow>
<mi>d</mi>
<mi>x</mi>
<mi>d</mi>
<mi>y</mi>
<mo>+</mo>
</mrow>
</mtd>
</mtr>
<mtr>
<mtd>
<mrow>
<msub>
<mi>&lambda;</mi>
<mn>2</mn>
</msub>
<msub>
<mo>&Integral;</mo>
<mi>&Omega;</mi>
</msub>
<msup>
<mrow>
<mo>(</mo>
<msub>
<mi>u</mi>
<mn>0</mn>
</msub>
<mo>(</mo>
<mrow>
<mi>x</mi>
<mo>,</mo>
<mi>y</mi>
</mrow>
<mo>)</mo>
<mo>-</mo>
<msub>
<mi>u</mi>
<mn>2</mn>
</msub>
<mo>)</mo>
</mrow>
<mn>2</mn>
</msup>
<mrow>
<mo>(</mo>
<mn>1</mn>
<mo>-</mo>
<mi>H</mi>
<mo>(</mo>
<mi>&phi;</mi>
<mo>)</mo>
<mo>)</mo>
</mrow>
<mi>d</mi>
<mi>x</mi>
<mi>d</mi>
<mi>y</mi>
<mo>+</mo>
</mrow>
</mtd>
</mtr>
<mtr>
<mtd>
<mrow>
<mi>&mu;</mi>
<msub>
<mo>&Integral;</mo>
<mi>&Omega;</mi>
</msub>
<mi>H</mi>
<mrow>
<mo>(</mo>
<mi>&phi;</mi>
<mo>)</mo>
</mrow>
<mi>d</mi>
<mi>x</mi>
<mi>d</mi>
<mi>y</mi>
<mo>+</mo>
<mi>v</mi>
<msub>
<mo>&Integral;</mo>
<mi>&Omega;</mi>
</msub>
<mi>&delta;</mi>
<mrow>
<mo>(</mo>
<mi>&phi;</mi>
<mo>)</mo>
</mrow>
<mo>|</mo>
<mrow>
<mo>&dtri;</mo>
<mi>&phi;</mi>
</mrow>
<mo>|</mo>
<mi>d</mi>
<mi>x</mi>
<mi>d</mi>
<mi>y</mi>
</mrow>
</mtd>
</mtr>
</mtable>
<mo>-</mo>
<mo>-</mo>
<mo>-</mo>
<mrow>
<mo>(</mo>
<mn>12</mn>
<mo>)</mo>
</mrow>
</mrow>
3
式中，H(φ)是Heaviside函数，其形式为：δ(φ)为Dirac函数，其形式为u0(x,y)是待处理影像区域某一点灰度值，为当前点梯度的模，系数λ1，λ2＞0，μ，ν≥0为固定参数，一般取λ1＝λ2＝1，μ＝0，ν＝1，(12)式对应的偏微分方程：
<mrow>
<mfrac>
<mrow>
<mo>&part;</mo>
<mi>&phi;</mi>
</mrow>
<mrow>
<mo>&part;</mo>
<mi>t</mi>
</mrow>
</mfrac>
<mo>=</mo>
<mi>&delta;</mi>
<mrow>
<mo>(</mo>
<mi>&phi;</mi>
<mo>)</mo>
</mrow>
<mo>&CenterDot;</mo>
<mo>&lsqb;</mo>
<mi>&mu;</mi>
<mo>&dtri;</mo>
<mo>&CenterDot;</mo>
<mfrac>
<mrow>
<mo>&dtri;</mo>
<mi>&phi;</mi>
</mrow>
<mrow>
<mo>|</mo>
<mrow>
<mo>&dtri;</mo>
<mi>&phi;</mi>
</mrow>
<mo>|</mo>
</mrow>
</mfrac>
<mo>-</mo>
<mi>v</mi>
<mo>-</mo>
<msub>
<mi>&lambda;</mi>
<mn>1</mn>
</msub>
<msup>
<mrow>
<mo>&lsqb;</mo>
<msub>
<mi>u</mi>
<mn>0</mn>
</msub>
<mrow>
<mo>(</mo>
<mi>x</mi>
<mo>,</mo>
<mi>y</mi>
<mo>)</mo>
</mrow>
<mo>-</mo>
<msub>
<mi>c</mi>
<mn>1</mn>
</msub>
<mo>&rsqb;</mo>
</mrow>
<mn>2</mn>
</msup>
<mo>+</mo>
<msub>
<mi>&lambda;</mi>
<mn>2</mn>
</msub>
<msup>
<mrow>
<mo>&lsqb;</mo>
<msub>
<mi>u</mi>
<mn>0</mn>
</msub>
<mrow>
<mo>(</mo>
<mi>x</mi>
<mo>,</mo>
<mi>y</mi>
<mo>)</mo>
</mrow>
<mo>-</mo>
<msub>
<mi>c</mi>
<mn>2</mn>
</msub>
<mo>&rsqb;</mo>
</mrow>
<mn>2</mn>
</msup>
<mo>&rsqb;</mo>
<mo>-</mo>
<mo>-</mo>
<mo>-</mo>
<mrow>
<mo>(</mo>
<mn>13</mn>
<mo>)</mo>
</mrow>
</mrow>
其中，参数c1、c2根据下式得到，
<mfenced open = "" close = "">
<mtable>
<mtr>
<mtd>
<mrow>
<msub>
<mi>c</mi>
<mn>1</mn>
</msub>
<mo>=</mo>
<mfrac>
<mrow>
<msub>
<mo>&Integral;</mo>
<mi>&Omega;</mi>
</msub>
<msub>
<mi>u</mi>
<mn>0</mn>
</msub>
<mrow>
<mo>(</mo>
<mi>x</mi>
<mo>,</mo>
<mi>y</mi>
<mo>)</mo>
</mrow>
<mi>H</mi>
<mrow>
<mo>(</mo>
<mi>&phi;</mi>
<mo>)</mo>
</mrow>
<mi>d</mi>
<mi>x</mi>
<mi>d</mi>
<mi>y</mi>
</mrow>
<mrow>
<msub>
<mo>&Integral;</mo>
<mi>&Omega;</mi>
</msub>
<mi>H</mi>
<mrow>
<mo>(</mo>
<mi>&phi;</mi>
<mo>)</mo>
</mrow>
<mi>d</mi>
<mi>x</mi>
<mi>d</mi>
<mi>y</mi>
</mrow>
</mfrac>
</mrow>
</mtd>
<mtd>
<mrow>
<msub>
<mi>c</mi>
<mn>2</mn>
</msub>
<mo>=</mo>
<mfrac>
<mrow>
<msub>
<mo>&Integral;</mo>
<mi>&Omega;</mi>
</msub>
<msub>
<mi>u</mi>
<mn>0</mn>
</msub>
<mrow>
<mo>(</mo>
<mi>x</mi>
<mo>,</mo>
<mi>y</mi>
<mo>)</mo>
</mrow>
<mo>&lsqb;</mo>
<mn>1</mn>
<mo>-</mo>
<mi>H</mi>
<mrow>
<mo>(</mo>
<mi>&phi;</mi>
<mo>)</mo>
</mrow>
<mo>&rsqb;</mo>
<mi>d</mi>
<mi>x</mi>
<mi>d</mi>
<mi>y</mi>
</mrow>
<mrow>
<msub>
<mo>&Integral;</mo>
<mi>&Omega;</mi>
</msub>
<mo>&lsqb;</mo>
<mn>1</mn>
<mo>-</mo>
<mi>H</mi>
<mrow>
<mo>(</mo>
<mi>&phi;</mi>
<mo>)</mo>
</mrow>
<mo>&rsqb;</mo>
<mi>d</mi>
<mi>x</mi>
<mi>d</mi>
<mi>y</mi>
</mrow>
</mfrac>
</mrow>
</mtd>
</mtr>
</mtable>
</mfenced>
根据式(13)，采用C-V水平集演化方法分割提取建筑物区域，作为影像的建筑物分割对像viSeg。
2.根据权利要求1所述的一种无人机多重叠遥感影像的建筑物轮廓线提取方法，其特征在于，所述步骤一包括以下步骤：
(1.1)利用先验信息对无人机多重叠遥感影像进行预处理：
(1.2)在步骤(1.1)的基础上进行空三摄影测量，利用空三勾网，求出每张影像的外方位元素，并进行光束法的整体平差；
(1.3)根据影像分组，在步骤(1.2)的基础上利用GPU加速的PMVS算法进行快速的密集匹配，生成密集的三维点云，所生成的三维点云作为三维高程数据。
3.根据权利要求1所述的一种无人机多重叠遥感影像的建筑物轮廓线提取方法，其特征在于，所述步骤二包括以下步骤：
(2.1)特征提取：利用SIFT进行影像的特征提取；
(2.2)影像配准：先进行粗配准，利用k-d树搜索匹配的特征点；然后进行精配准，粗配准往往出现错误的匹配点，因此利用RANSAC算法剔除错误的匹配点；通过影像的配准，得到影像之间的变换矩阵；
(2.3)影像的拼接：通过(2.2)得到的变换矩阵进行影像的拼接；
(2.4)影像的融合：拼接后，利用双线性插值算法进行影像的融合。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201510025503.9A CN104484668B (zh) | 2015-01-19 | 2015-01-19 | 一种无人机多重叠遥感影像的建筑物轮廓线提取方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201510025503.9A CN104484668B (zh) | 2015-01-19 | 2015-01-19 | 一种无人机多重叠遥感影像的建筑物轮廓线提取方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN104484668A CN104484668A (zh) | 2015-04-01 |
| CN104484668B true CN104484668B (zh) | 2017-11-10 |
Family
ID=52759209
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201510025503.9A Active CN104484668B (zh) | 2015-01-19 | 2015-01-19 | 一种无人机多重叠遥感影像的建筑物轮廓线提取方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN104484668B (zh) |
Families Citing this family (48)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN104809689B (zh) * | 2015-05-15 | 2018-03-30 | 北京理工大学深圳研究院 | 一种基于轮廓的建筑物点云模型底图配准方法 |
| CN105046201B (zh) * | 2015-06-19 | 2018-05-08 | 东南大学 | 一种基于形态学对建筑物图像快速识别的方法 |
| CN105184854B (zh) * | 2015-08-24 | 2018-01-16 | 北京麦格天宝科技股份有限公司 | 针对地下空间扫描点云成果数据的快速建模方法 |
| CN105808926B (zh) * | 2016-03-02 | 2017-10-03 | 中国地质大学（武汉） | 一种基于gpu并行加速的预条件共轭梯度区域网平差方法 |
| CN105844707B (zh) * | 2016-03-15 | 2019-04-16 | 南京大学 | 基于城墙断面从LiDAR点云数据自动提取古城墙数据的方法 |
| CN106023312B (zh) * | 2016-05-13 | 2019-02-22 | 南京大学 | 基于航空LiDAR数据的三维建筑物模型自动重建方法 |
| CN106023224A (zh) * | 2016-05-30 | 2016-10-12 | 天水师范学院 | 一种中药材显微图像的pcnn自动分割方法 |
| CN106056614A (zh) * | 2016-06-03 | 2016-10-26 | 武汉大学 | 一种地面激光点云数据的建筑物分割与轮廓线提取方法 |
| CN106203335B (zh) * | 2016-07-11 | 2019-05-28 | 厦门大学 | 基于三维点云的标志牌可视度评价方法 |
| CN107633523B (zh) * | 2016-07-18 | 2021-04-16 | 巧夺天宫（深圳）科技有限公司 | 基于点云的提取建筑特征线方法和*** |
| CN106709897B (zh) * | 2016-12-28 | 2019-11-26 | 武汉大学 | 基于梯度领域的正射影像间最优拼接线寻找方法及*** |
| CN107240128B (zh) * | 2017-05-09 | 2020-09-11 | 北京理工大学 | 一种基于轮廓特征的x线片和彩色照片配准方法 |
| CN107103603A (zh) * | 2017-05-18 | 2017-08-29 | 北京道亨时代科技有限公司 | 一种倾斜测量场景的地物提取方法 |
| CN107393004A (zh) * | 2017-07-17 | 2017-11-24 | 北京数字绿土科技有限公司 | 一种获取输电线走廊中建筑物拆迁量的方法及装置 |
| EP3451203A1 (en) * | 2017-08-30 | 2019-03-06 | Dassault Systèmes | Computer-implemented method for computing an envelope for a building complying with shadow duration requirements |
| CN107514993B (zh) * | 2017-09-25 | 2019-11-05 | 同济大学 | 基于无人机的面向单体建筑建模的数据采集方法及*** |
| CN107741233A (zh) * | 2017-11-10 | 2018-02-27 | 邦鼓思电子科技（上海）有限公司 | 一种三维室外地图的构建方法 |
| CN108088422B (zh) * | 2018-01-24 | 2020-12-15 | 成都纵横自动化技术股份有限公司 | 一种确定真实重叠率的方法 |
| CN108571930A (zh) * | 2018-04-09 | 2018-09-25 | 刘默然 | 一种利用无人机的房屋测量方法 |
| CN108627112A (zh) * | 2018-05-09 | 2018-10-09 | 广州市杜格科技有限公司 | 车辆销轴距动态测量方法 |
| CN109242855B (zh) * | 2018-07-19 | 2020-08-11 | 中国科学院自动化研究所 | 基于多分辨率三维统计信息的屋顶分割方法、***及设备 |
| CN109146990B (zh) * | 2018-08-08 | 2023-02-24 | 广州市城市规划勘测设计研究院 | 一种建筑轮廓的计算方法 |
| CN109584262A (zh) * | 2018-12-15 | 2019-04-05 | 中国科学院深圳先进技术研究院 | 基于遥感影像的云检测方法、装置及电子设备 |
| CN109671109B (zh) * | 2018-12-25 | 2021-05-07 | 中国人民解放军61540部队 | 密集点云生成方法及*** |
| CN110097588B (zh) * | 2019-04-22 | 2021-01-15 | 西安交通大学 | 一种航发叶片陶瓷型芯点云模型的修型边缘提取方法 |
| CN110111349B (zh) * | 2019-04-22 | 2020-12-08 | 西安交通大学 | 一种基于点云的非刚体复杂构件高精度边缘提取方法 |
| CN112146564B (zh) * | 2019-06-28 | 2022-04-15 | 先临三维科技股份有限公司 | 三维扫描方法、装置、计算机设备和计算机可读存储介质 |
| CN110415337B (zh) * | 2019-07-12 | 2021-03-23 | 广州大学 | 一种基于影像直线特征的建筑物屋顶面重建方法和装置 |
| CN112241661A (zh) * | 2019-07-17 | 2021-01-19 | 临沂大学 | 一种结合机载LiDAR点云数据和航空影像的城市地物精细化分类方法 |
| CN110570428B (zh) * | 2019-08-09 | 2023-07-07 | 浙江合信地理信息技术有限公司 | 一种从大规模影像密集匹配点云分割建筑物屋顶面片的方法及*** |
| CN110543872B (zh) * | 2019-09-12 | 2023-04-18 | 云南省水利水电勘测设计研究院 | 一种基于全卷积神经网络的无人机影像建筑物屋顶提取方法 |
| CN110910387B (zh) * | 2019-10-09 | 2022-03-04 | 西安理工大学 | 一种基于显著性分析的点云建筑物立面窗户提取方法 |
| CN110879978A (zh) * | 2019-11-13 | 2020-03-13 | 广西中煤地质有限责任公司 | 一种无人机倾斜摄影三维模型的建筑轮廓线自动提取方法 |
| CN111209920B (zh) * | 2020-01-06 | 2022-09-23 | 桂林电子科技大学 | 一种复杂动态背景下飞机检测方法 |
| CN111652241B (zh) * | 2020-02-17 | 2023-04-28 | 中国测绘科学研究院 | 融合影像特征与密集匹配点云特征的建筑物轮廓提取方法 |
| CN111523391B (zh) * | 2020-03-26 | 2021-11-05 | 上海刻羽信息科技有限公司 | 建筑物的识别方法、***、电子设备及可读存储介质 |
| CN111383355B (zh) * | 2020-04-03 | 2023-06-27 | 如你所视(北京)科技有限公司 | 三维点云补全方法、装置和计算机可读存储介质 |
| CN111322950B (zh) * | 2020-04-17 | 2021-08-17 | 易思维（杭州）科技有限公司 | 利用线结构光传感器定位圆柱***置的方法及其用途 |
| CN111814715B (zh) * | 2020-07-16 | 2023-07-21 | 武汉大势智慧科技有限公司 | 一种地物分类方法及装置 |
| CN112132029B (zh) * | 2020-09-23 | 2023-07-11 | 中国地震局地震预测研究所 | 一种面向地震应急响应的无人机遥感影像快速定位方法 |
| CN112232248B (zh) * | 2020-10-22 | 2023-04-07 | 中国人民解放军战略支援部队信息工程大学 | 一种多线LiDAR点云数据平面特征的提取方法及装置 |
| CN112748700A (zh) * | 2020-12-18 | 2021-05-04 | 深圳市显控科技股份有限公司 | 数控代码生成方法、装置、电子设备及存储介质 |
| CN113296543B (zh) * | 2021-07-27 | 2021-09-28 | 成都睿铂科技有限责任公司 | 航拍航线规划方法及*** |
| CN114782826B (zh) * | 2022-06-20 | 2022-11-18 | 绵阳天仪空间科技有限公司 | 一种灾后建筑物的安全监测***及方法 |
| CN115546421B (zh) * | 2022-12-02 | 2023-03-24 | 中南大学 | 一种点线双域互增强建筑几何轮廓线重建方法 |
| CN116229339B (zh) * | 2023-02-06 | 2024-07-02 | 南京航空航天大学 | 一种基于语义分割和由粗到精策略的船闸水位检测方法 |
| CN117274605B (zh) * | 2023-11-20 | 2024-03-01 | 北京飞渡科技股份有限公司 | 一种从无人机拍摄的照片中提取水域轮廓的方法及装置 |
| CN118212366A (zh) * | 2024-05-21 | 2024-06-18 | 中国科学院空天信息创新研究院 | 基于多重遥感影像的运动目标三维重建方法及其装置 |
Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| WO2009131108A1 (ja) * | 2008-04-23 | 2009-10-29 | 株式会社パスコ | 建物屋根輪郭認識装置、建物屋根輪郭認識方法、及び建物屋根輪郭認識プログラム |
| CN102930540A (zh) * | 2012-10-26 | 2013-02-13 | 中国地质大学（武汉） | 城市建筑物轮廓检测的方法及*** |
| CN103020342A (zh) * | 2012-12-04 | 2013-04-03 | 南京大学 | 一种从地面LiDAR数据中提取建筑物轮廓和角点的方法 |
| CN104091369A (zh) * | 2014-07-23 | 2014-10-08 | 武汉大学 | 一种无人机遥感影像建筑物三维损毁检测方法 |
-
2015
- 2015-01-19 CN CN201510025503.9A patent/CN104484668B/zh active Active
Patent Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| WO2009131108A1 (ja) * | 2008-04-23 | 2009-10-29 | 株式会社パスコ | 建物屋根輪郭認識装置、建物屋根輪郭認識方法、及び建物屋根輪郭認識プログラム |
| CN102930540A (zh) * | 2012-10-26 | 2013-02-13 | 中国地质大学（武汉） | 城市建筑物轮廓检测的方法及*** |
| CN103020342A (zh) * | 2012-12-04 | 2013-04-03 | 南京大学 | 一种从地面LiDAR数据中提取建筑物轮廓和角点的方法 |
| CN104091369A (zh) * | 2014-07-23 | 2014-10-08 | 武汉大学 | 一种无人机遥感影像建筑物三维损毁检测方法 |
Also Published As
| Publication number | Publication date |
|---|---|
| CN104484668A (zh) | 2015-04-01 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN104484668B (zh) | 一种无人机多重叠遥感影像的建筑物轮廓线提取方法 | |
| US20220028163A1 (en) | Computer Vision Systems and Methods for Detecting and Modeling Features of Structures in Images | |
| CN110264567B (zh) | 一种基于标记点的实时三维建模方法 | |
| CN109345574B (zh) | 基于语义点云配准的激光雷达三维建图方法 | |
| Zhou et al. | Seamless fusion of LiDAR and aerial imagery for building extraction | |
| CN106595659A (zh) | 城市复杂环境下多无人机视觉slam的地图融合方法 | |
| CN107063228A (zh) | 基于双目视觉的目标姿态解算方法 | |
| CN111998862B (zh) | 一种基于bnn的稠密双目slam方法 | |
| Cosido et al. | Hybridization of convergent photogrammetry, computer vision, and artificial intelligence for digital documentation of cultural heritage-a case study: the magdalena palace | |
| CN112465849B (zh) | 一种无人机激光点云与序列影像的配准方法 | |
| KR101549155B1 (ko) | 라이다 자료를 활용한 구조물의 직선경계 추출방법 | |
| Alcantarilla et al. | Large-scale dense 3D reconstruction from stereo imagery | |
| Zakharov et al. | Automatic building detection from satellite images using spectral graph theory | |
| Wu et al. | A new stereo dense matching benchmark dataset for deep learning | |
| Chaloeivoot et al. | Building detection from terrestrial images | |
| CN112767459A (zh) | 基于2d-3d转换的无人机激光点云与序列影像配准方法 | |
| Brunet et al. | AI4GEO: A Path From 3D Model to Digital Twin | |
| Boerner et al. | Brute force matching between camera shots and synthetic images from point clouds | |
| Wang et al. | Automated mosaicking of UAV images based on SFM method | |
| CN107784666B (zh) | 基于立体影像的地形地物三维变化检测和更新方法 | |
| Verhoeven | Getting computer vision airborne: using structure from motion for accurate orthophoto production | |
| Chen et al. | True orthophoto generation using multi-view aerial images | |
| Previtali et al. | An automatic multi-image procedure for accurate 3D object reconstruction | |
| Zhou et al. | Occlusion detection for urban aerial true orthoimage generation | |
| Bai et al. | Application of unmanned aerial vehicle multi-vision image 3D modeling in geological disasters |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| C06 | Publication | ||
| PB01 | Publication | ||
| C10 | Entry into substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |