CN103093466A - 基于LiDAR点云和影像的建筑物三维变化检测方法 - Google Patents

Abstract

本发明提出一种基于LiDAR点云和影像的建筑物三维变化检测方法，通过LiDAR数据滤波、插值生成高精度DSM，然后对DSM进行差值得到候选变化区域，然后通过相关公式将候选变化区域对应到二维影像上，再通过边缘提取、直线特征匹配、边缘差值等提取影像上的变化区域。本发明的优点是能够充分利用LiDAR点云数据提供的高程信息，同时结合二维影像信息进行人工建筑物的三维变化检测，显著提高了变化检测的精度、广度和深度。

Description

基于LiDAR点云和影像的建筑物三维变化检测方法

技术领域

本发明属于测绘科学与技术领域，涉及一种基于LiDAR点云和影像数据的建筑物三维变化检测方法，主要应用于人工建筑物变化信息提取与处理等领域。

背景技术

变化检测的数据源有遥感影像，正射影像，矢量地图和LiDAR（激光雷达）点云数据等。其中，LiDAR不受日照和天气条件的限制，能全天候对地观测，快速获取低成本、高精度的地表的三维坐标，为准确提取建筑物特征提供了良好的数据源。特别地，相比于其他方法，LiDAR点云数据可提供人工建筑物高精度的高程信息。但是由于机载激光雷达点云数据缺乏纹理信息，同时数据密度分布不均匀并且在三维空间中呈不规则性、不连续分布。因此直接利用LiDAR点云数据进行人工建筑物三维变化检测比较困难，尤其是在建筑物形状复杂、建筑物边界细节较多时，很难对建筑物进行准确、快速的变化检测。通过航空摄影获取的航空影像数据空间分辨率高，同时数据分布连续，具有十分丰富纹理信息，但是影像数据不能提供高精度的高程信息。故将LiDAR点云数据和影像数据结合进行建筑物的三维变化检测可以达到不同数据源的优势互补。两种数据源的融合研究对于人工建筑物高效率、高精度的三维变化检测，对于LiDAR数据后处理技术发展以满足实际应用都有十分重要的研究价值。现有的变化检测方法有：（1）基于代数运算的变化检测方法，包括图像差分、图像比值、图像回归方法、图像植被指数差分、变化矢量分析和背景相减等方法；（2）基于图像变换的变化检测方法，包括主成分分析方法(PCA)、穗帽变换方法(KT)、典型相关分析等方法；（3）基于图像分类的变化检测方法，包括分类后比较方法，多时相图像直接分类方法（也称光谱/时相分类）；（4）基于对象的变化检测方法；（5）基于统计模型的变化检测方法；（6）基于小波变换的变化检测方法。尽管上述各种变化检测方法采用的变化检测基本单元不同，变换检测的策略差别很大，但它们都只使用二维影像的平面信息，没有充分顾及高程变化，而人工建筑物一个十分重要的变化就表现在高度上的变化上。因此迫切需要寻求考虑高程变化的变化检测方法。

发明内容

本发明的目的是提供一种顾及高程信息变化的人工建筑物三维变化检测方法，它能够克服上述现有变化检测方法技术的不足，满足人工建筑物三维变化检测应用的需求。

本发明的技术方案为一种基于LiDAR点云和影像的建筑物三维变化检测方法，包括以下步骤，

步骤1，对两个不同时期的LiDAR点云数据分别进行滤波和插值处理，生成两个不同时期的DSM；

步骤2，求取步骤1所得两个不同时期的DSM的差值，得到DSM差值影像；

步骤3，对步骤2所得DSM差值影像处理，获取DSM差值影像上的侯选变化区域；

步骤4，根据DSM和影像坐标对应关系，解算出步骤3所得DSM差值影像上的侯选变化区域在两个不同时期相应原始航空影像上分别相对应的候选变化区域；

步骤5，对步骤4所得两个不同时期相应航空影像上候选变化区域分别进行边缘特征提取，得到两幅边缘图像；

步骤6，对步骤5所得两幅边缘图像分别进行直线特征提取和匹配，得到两幅直线特征提取和匹配后的边缘图像；

步骤7，对步骤6所得两幅直线特征提取和匹配后的边缘图像进行差值运算，同时剔除不匹配的直线及直线区域外的边缘，得到候选变化区域中的变化图像。

而且，步骤4所述DSM和影像坐标对应关系如下式，

$x=-f\frac{a_1(X-\mathrm{Xs})+b_1(Y-\mathrm{Ys})+c_1(Z-\mathrm{Zs})}{a_3(X-\mathrm{Xs})+b_3(Y-\mathrm{Ys})+c_3(Z-\mathrm{Zs})}$

$y=-f\frac{a_2(X-\mathrm{Xs})+b_2(Y-\mathrm{Ys})+c_2(Z-\mathrm{Zs})}{a_3(X-\mathrm{Xs})+b_3(Y-\mathrm{Ys})+c_3(Z-\mathrm{Zs})}$

其中，a₁、a₂、a₃、b₁、b₂、b₃、c₁、c₂、c₃为影像外方位元素生成的9个旋转矩阵元素；(X_s,Y_s,Z_s)为影像的三个外方位线元素；(X,Y,Z)为DSM上的地面点物方坐标；(x,y)为影像上的像方坐标。

而且，步骤6所述对步骤5所得两幅边缘图像分别进行直线特征提取和匹配，匹配实现方式如下，

设两个不同时期分别记为T1时期和T2时期，设在T2时期影像的某一个候选变化区域内的直线表示为集合S_n，在T1时期影像上的同一候选变化区域内的直线表示为集合S_o，

对于集合S_n中的任意一条直线(ρ，θ)，采用以下公式计算到集合S_o中所有直线的距离，

$D=\sqrt{{(\mathrm{\ρ}-{\mathrm{\ρ}}^{\′})}^2+{(\mathrm{\θ}-{\mathrm{\θ}}^{\′})}^2}$

其中，D表示集合S_n上某直线(ρ，θ)与集合S_o上某直线(ρ′，θ')之间的距离；

如果计算所得距离中的最小者小于某个给定的阈值T，则认为该距离所对应的两条直线是匹配的，否则认为集合S_o中没有对应的直线与直线(ρ，θ)匹配。

本发明的优点是不仅充分利用了图像匹配检测的特点(直线匹配)，而且充分利用了人工建筑物高度变化这一重要特点。通过对LiDAR数据的处理分析，利用生成的DSM图像来更加准确地检测出人工建筑物的变化。该方法在整个变化检测过程中不需要人工进行干预，在时间上能满足更新速度快的数据库和地理信息***的要求，而且可以进行定量的分析，同时也能对检测出的变化区域进行变化性质的判断。

附图说明

图1是本发明实施例的基于LiDAR数据和影像的建筑物三维变化检测流程图。

具体实施方式

以下结合附图和实施例详细说明本发明技术方案。

参见图1，本发明提供人工建筑物三维变化检测方法，通过LiDAR数据滤波、插值生成高精度DSM，然后对DSM进行差值得到候选变化区域，然后通过相关公式将候选变化区域对应到二维影像上，再通过边缘提取、直线特征匹配、边缘差值等提取影像上的变化区域。输入数据为两期机载激光点云和两期航空影像。具体实施时可采用计算机软件技术实现自动运行。

实施例的具体实现方法包括以下步骤：

步骤1，高精度DSM生成：对两个不同时期的LiDAR点云数据分别进行滤波和插值处理，生成两幅DSM（数字表面高程模型）。

实施例对新旧两个不同时期的原始LiDAR点云数据进行处理，分别记为T1时期LiDAR数据和T2时期LiDAR数据。首先对两期原始LiDAR点云进行滤波以剔除粗差点（极高点、极低点和噪声点），然后对剔除粗差点的点云数据进行三角网插值（TIN插值），同时采用一定的配色方案，不同颜色代表不同的高程值，进而生成两个不同时期的高精度DSM。T1时期LiDAR数据和T2时期LiDAR数据的对应DSM分别记为T1时期DSM和T2时期DSM。

步骤2，生成DSM差值影像：根据高程求取两个不同时期的DSM的差值，得到DSM差值影像。

实施例利用高程的变化，对生成的两期高精度DSM进行差值运算，得到DSM差值影像。

步骤3，DSM候选变化区域生成：对DSM差值影像进行二值化处理、边缘跟踪等处理手段，获取DSM差值影像上的侯选变化区域。

实施例设置合适的阈值，对DSM差值影像进行二值化处理得到二值图像，然后对该二值图像进行形态学边界跟踪得到DSM差值影像上的候选变化区域。

步骤4，航空影像候选变化域生成：根据DSM和影像坐标对应关系，解算出步骤3所得DSM差值影像上的侯选变化区域在两个不同时期相应原始航空影像上分别相对应的候选变化区域。

为便于实施参考，提供涉及的DSM和影像坐标转化数学公式如式（1）所示：

$x=-f\frac{a_1(X-\mathrm{Xs})+b_1(Y-\mathrm{Ys})+c_1(Z-\mathrm{Zs})}{a_3(X-\mathrm{Xs})+b_3(Y-\mathrm{Ys})+c_3(Z-\mathrm{Zs})}$

$y=-f\frac{a_2(X-\mathrm{Xs})+b_2(Y-\mathrm{Ys})+c_2(Z-\mathrm{Zs})}{a_3(X-\mathrm{Xs})+b_3(Y-\mathrm{Ys})+c_3(Z-\mathrm{Zs})}---\left(1\right)$

其中，a₁、a₂、a₃、b₁、b₂、b₃、c₁、c₂、c₃为影像外方位元素生成的9个旋转矩阵元素；(X_s,Y_s,Z_s)为影像的三个外方位线元素；(X,Y,Z)为DSM上的地面点物方坐标；(x,y)为影像上的像方坐标。

实施例通过公式（1）中DSM和影像坐标投影换算关系，解算出DSM中侯选变化区域在两期相应原始航空影像上相对应的区域，即可得到T1时期影像数据和T2时期影像数据的候选变化区域。

步骤5，航空影像候选变化区域边缘提取：对步骤4所得两个不同时期相应航空影像上候选变化区域分别进行边缘特征提取，得到两幅边缘图像。

由于Canny算子可以在检测精度和抑制噪声方面取得很好的平衡，故实施例采用Canny算子对两期航空影像候选变化区域进行边缘特征提取，得到T1时期边缘图像和T2时期边缘图像。Canny算子具体实现为现有技术。

步骤6，边缘图像线特征提取与匹配：对步骤5所得两幅边缘图像分别进行直线特征提取和匹配，得到两幅直线特征提取和匹配后的边缘图像（可简称为匹配后边缘图像）。

由于噪声和光照不均等因素，边缘图像中有的边缘点不连续，必须通过边缘连接将它们转化为有意义的边缘。实施例使用Hough变换提取边缘图像的直线特征，然后按照公式（2）计算两期边缘图像待匹配直线间距离D，如果D小于给定的阈值T，则认为直线匹配成功，否则认为直线匹配失败。Hough变换具体实现为现有技术。

为便于实施参考起见，提供实施例使用Hough变换提取新旧影像候选变化区域的直线特征，进而进行直线特征匹配，具体匹配策略如下：

设在T2时期影像的某一个候选变化区域内的直线表示为集合S_n：

s_n={(ρ₀,θ₀),(ρ₁,θ₁)......(ρ_n,θ_n)}

其中，ρ_i(i＝0,1,2......n)表示T2时期影像候选变化区域中第i条直线到原点的垂直距离，θ_i(i＝0,1,2……n)表示x轴到第i条直线的垂线的角度，θ_i∈[-90°,90°]。

在T1时期影像上的同一候选变化区域内的直线表示为集合S_o：

$S_o=\{({\mathrm{\ρ}}_0^{\′},{\mathrm{\θ}}_0^{\′}),({\mathrm{\ρ}}_1^{\′},{\mathrm{\θ}}_1^{\′})......({\mathrm{\ρ}}_n^{\′},{\mathrm{\θ}}_n^{\′})\}$

其中，ρ'_i(i＝0,1,2......n)表示T1时期影像候选变化区域中第i条直线到原点的垂直距离，θ′_i(i＝0,1,2......n)表示x轴到第i条直线的垂线的角度，θ'_i∈[-90°,90°]。

计算集合S_n上某直线(ρ，θ)与集合S_o上某直线(ρ′，θ')之间的距离D：

$D=\sqrt{{(\mathrm{\ρ}-{\mathrm{\ρ}}^{\′})}^2+{(\mathrm{\θ}-{\mathrm{\θ}}^{\′})}^2}---\left(2\right)$

对于集合S_n中的任意一条直线，计算到集合S_o中所有直线的距离D，如果所有距离中的最小距离小于某个给定的阈值T，则认为该距离所对应的两条直线是匹配的，否则认为集合S_o中没有对应的直线与直线(ρ，θ)匹配。

步骤7，变化区域图像生成：对步骤6所得两幅直线特征提取和匹配后的边缘图像进行差值运算，同时删除匹配失败的直线及直线区域外的边缘，得到候选变化区域中的变化图像，即变化区域图像。

本文中所描述的具体实施例仅仅是对本发明精神作举例说明。本发明所属技术领域的技术人员可以对所描述的具体实施例做各种各样的修改或补充或采用类似的方式替代，但并不会偏离本发明的精神或者超越所附权利要求书所定义的范围。

Claims (3)

1.一种基于LiDAR点云和影像的建筑物三维变化检测方法，其特征在于：包括以下步骤，

步骤1，对两个不同时期的LiDAR点云数据分别进行滤波和插值处理，生成两个不同时期的DSM；

步骤2，求取步骤1所得两个不同时期的DSM的差值，得到DSM差值影像；

步骤3，对步骤2所得DSM差值影像处理，获取DSM差值影像上的侯选变化区域；

步骤4，根据DSM和影像坐标对应关系，解算出步骤3所得DSM差值影像上的侯选变化区域在两个不同时期相应原始航空影像上分别相对应的候选变化区域；

步骤5，对步骤4所得两个不同时期相应航空影像上候选变化区域分别进行边缘特征提取，得到两幅边缘图像；

步骤6，对步骤5所得两幅边缘图像分别进行直线特征提取和匹配，得到两幅直线特征提取和匹配后的边缘图像；

步骤7，对步骤6所得两幅直线特征提取和匹配后的边缘图像进行差值运算，同时剔除不匹配的直线及直线区域外的边缘，得到候选变化区域中的变化图像。

2.根据权利要求1所述基于LiDAR点云和影像的建筑物三维变化检测方法，其特征在于：步骤4所述DSM和影像坐标对应关系如下式，

$x=-f\frac{a_1(X-\mathrm{Xs})+b_1(Y-\mathrm{Ys})+c_1(Z-\mathrm{Zs})}{a_3(X-\mathrm{Xs})+b_3(Y-\mathrm{Ys})+c_3(Z-\mathrm{Zs})}$

$y=-f\frac{a_2(X-\mathrm{Xs})+b_2(Y-\mathrm{Ys})+c_2(Z-\mathrm{Zs})}{a_3(X-\mathrm{Xs})+b_3(Y-\mathrm{Ys})+c_3(Z-\mathrm{Zs})}$

其中，a₁、a₂、a₃、b₁、b₂、b₃、c₁、c₂、c₃为影像外方位元素生成的9个旋转矩阵元素；(X_s,Y_s,Z_s)为影像的三个外方位线元素；(X,Y,Z)为DSM上的地面点物方坐标；(x,y)为影像上的像方坐标。

3.根据权利要求1或2所述基于LiDAR点云和影像的建筑物三维变化检测方法，其特征在于：步骤6所述对步骤5所得两幅边缘图像分别进行直线特征提取和匹配，匹配实现方式如下，

设两个不同时期分别记为T1时期和T2时期，设在T2时期影像的某一个候选变化区域内的直线表示为集合S_n，在T1时期影像上的同一候选变化区域内的直线表示为集合S_o，

对于集合S_n中的任意一条直线(ρ，θ)，采用以下公式计算到集合S_o中所有直线的距离，

$D=\sqrt{{(\mathrm{\ρ}-{\mathrm{\ρ}}^{\′})}^2+{(\mathrm{\θ}-{\mathrm{\θ}}^{\′})}^2}$

其中，D表示集合S_n上某直线(ρ，θ)与集合S_o上某直线(ρ′，θ')之间的距离；

如果计算所得距离中的最小者小于某个给定的阈值T，则认为该距离所对应的两条直线是匹配的，否则认为集合S_o中没有对应的直线与直线(ρ，θ)匹配。

Images