[发明专利]一种无人机多重叠遥感影像的建筑物轮廓线提取方法

摘要

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

主权项

一种无人机多重叠遥感影像的建筑物轮廓线提取方法，其特征在于，包括以下步骤：步骤一，利用空三对无人机遥感影像进行平差，同时利用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,利用绿色和蓝色通道定义的颜色不变量公式为：ψg(x,y,z)=4π×arctan(Ig(x,y,z)-Ib(x,y,z)Ig(x,y,z)+Ib(x,y,z))---(1)]]>其中，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=Σi=07(ax+by+c-z)2---(3)]]>要使得S最小，应满足：其中，a、b、c表示式(2)平面方程的参数，和分别表示S对a、b和c求导；由此可得下列正规方程组：Σxi2ΣxiyiΣxiΣxiyiΣyi2ΣyiΣxiΣyinabc=ΣxiziΣyiziΣzi,i=0,1,...7;n=8---(4)]]>(4.3)根据公式(4)求得任何内部点在8邻域内的方差和SSD：SSD=Σ(pk)∈M(zk-hk)2---(5)]]>其中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)方程如下：(x1-x0)2+(y1-y0)2=α2(x2-x0)2+(y2-y0)2=α2---(6)]]>直接求取此方程比较困难，因此利用测绘学中的距离交汇算法得：x0=x1+x2-x12+H(y2-y1)y0=y1+y2-y12+H(x2-x1)---(7)]]>其中，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=-f0x00fy0001RT100-Xs010-Ys001-Zs---(8)]]>其中，P即投影矩阵P3×4，f为影像的焦距，x0和y0为光轴距离光心在水平方向和垂直方向的偏心距；Xs、Ys、Zs为相机中心在世界坐标系中的坐标，RT表示3×3的旋转矩阵；投影公式为：xzyzz=PXYZ1---(9)]]>利用式(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)目标的初始轮廓和先验形状均利用运动目标区域进行表达：Esh(φsi,φPi)=∫ΩBi(φsi(x,y)-φPi(x,y))2dxdy---(11)]]>基于水平集C‑V模型的能量函数为：E(φ,c1,c2)=λ1∫Ω(u0(x,y)-u1)2H(φ)dxdy+λ2∫Ω(u0(x,y)-u2)2(1-H(φ))dxdy+μ∫ΩH(φ)dxdy+v∫Ωδ(φ)|▿φ|dxdy---(12)]]>式中，H(φ)是Heaviside函数，其形式为：δ(φ)为Dirac函数，其形式为u0(x,y)是待处理影像区域某一点灰度值，为当前点梯度的模，系数λ1，λ2＞0，μ，ν≥0为固定参数，一般取λ1＝λ2＝1，μ＝0，ν＝1，(12)式对应的偏微分方程：∂φ∂t=δ(φ)·[μ▿·▿φ|▿φ|-v-λ1[u0(x,y)-c1]2+λ2[u0(x,y)-c2]2]---(13)]]>其中，参数c1、c2根据下式得到，c1=∫Ωu0(x,y)H(φ)dxdy∫ΩH(φ)dxdyc2=∫Ωu0(x,y)[1-H(φ)]dxdy∫Ω[1-H(φ)]dxdy]]>根据式(13)，采用C‑V水平集演化方法分割提取建筑物区域，作为影像的建筑物分割对像viSeg。