CN106780712B - 联合激光扫描和影像匹配的三维点云生成方法 - Google Patents

Abstract

本发明涉及一种联合激光扫描和影像匹配的三维点云生成方法，它包括如下步骤：1、激光LiDAR点云完整度评估；2、弱/无反射区域三维点云生成；3、线特征地物的三维点云生成；4、激光LiDAR点云和密集匹配点云融合。本发明充分利用激光扫描技术和影像密集匹配技术各自的优势，根据影像密集匹配技术，解决激光扫描在弱/无反射区域和线状地物区域表现较弱的问题；根据原始的LiDAR点云，极大程度上减少密集匹配的计算量，能够快速获得完整、高精度、密集的三维点云产品。

Description

联合激光扫描和影像匹配的三维点云生成方法

技术领域

本发明涉及三维点云获取技术领域，具体涉及一种联合激光扫描和影像匹配的三维点云生成方法。

背景技术

21世纪是一个信息化的时代，突出表现在将人类世界以信息化的方式进行显示、统计、分析和应用。美国前副总统戈尔在1998年的演讲中，提出了“数字地球”的概念，掀起了一股数字地球研究的浪潮。2005年，谷歌公司推出了GoogleEarth产品，将航片、卫片数据映射到虚拟的数字地球上，便于用户浏览和使用世界各地的地理信息产品。在全世界正从数字地球向智慧地球发展的过程中，我国也在推动“数字中国”、“智慧城市”等的快速发展，在为政府部门、企事业单位和社会公众等提供公共管理、突发事件应急、科学决策等服务方面具有重要意义。

数字表面模型(Digital Surface Model,DSM)是以数字化的形式，表达地球自然表面和人工地物的三维模型，是“智慧城市”重要的显示和分析工具。近十年来，随着“数字中国”、“智慧城市”的快速发展，对DSM生成的精度、速度和分辨率提出了越来越高的要求。激光扫描和立体影像密集匹配技术是目前两种主流的地形地物三维信息获取手段，通过将相机、LiDAR等传感器搭载于车辆、无人机、航空飞机乃至卫星等不同高度的遥感平台，间接或直接地快速获取测区地表的三维测绘信息。

LiDAR(Light Detection And Ranging，激光探测与测量点云)系统由全球定位系统(GPS)、惯性导航系统(IMU)以及激光器(Laser Scanner)共同组成，以主动方式快速获取地面模型及空间点云信息，具有速度快、时效性强、作业范围大等优点。但是LiDAR系统成本昂贵，且点云密度较低(与对应的影像分辨率比较而言)，无法很好地描述线状地物；并且由于自然地表和建筑物的材质问题，在某些区域无法形成反射点，无法很好地为智慧城市等三维重建应用服务。

影像密集匹配通过覆盖一个测区的多张立体影像，根据同名光线对对相交的原理，从影像的二维信息恢复整个测区的三维空间信息。与LiDAR激光点云相比，立体影像具有密集匹配点云密度大、平面精度高、粗差剔除技术完善、影像上的地物几何特征更明显、影像数据获取成本低等突出优势。但是影像匹配的结果取决于地物纹理情况以及影像辐射质量。影像密集匹配技术在纹理贫乏区域、重复纹理区域表现较弱。较差的影像辐射质量以及阴影区域也会严重影响匹配结果。影像密集匹配往往速度较慢，很难达到实时的要求，严重制约数字城市的重建速度。

由于上述两种技术各自的不足，单纯地采用激光扫描，或者单纯地采用影像密集匹配，均无法很好地描述测区地表的三维信息。为了满足高精度、高分辨率、高实时性的三维重建的需求，需要将两种技术相结合，充分发挥各自的优势，得到完整可靠的点云产品，能够满足快速、大范围、高精度、高分辨率三维重建的需求，为数字地球、智慧城市的发展提供技术支撑，为社会的可持续发展服务。

发明内容

本发明的目的在于提供一种联合激光扫描和影像匹配的三维点云生成方法，该方法充分发挥激光扫描技术和影像密集匹配技术各自的优势，利用激光扫描技术直接快速获取测区的三维点云，对于线状地物、无反射点区域(如水面)等扫描困难区域，采用影像密集匹配技术获取扫描困难区域的三维点云，最终快速生成稠密的高精度的三维点云，为智慧城市等应用服务。

为解决上述技术问题，本发明公开的一种联合激光扫描和影像匹配的三维点云生成方法，其特征在于，它包括如下步骤：

步骤1：根据原始激光探测与测量点云数据，生成数字表面模型，提取数字表面模型中的无效区域，将该无效区域中最小外接矩形的四个角点，分别投影到三维重建测区的影像上，生成一个三维重建测区影像上的四边形区域，该四边形区域即可认为是激光探测与测量点云数据中的无效区域在三维重建测区影像上对应的范围；

步骤2：将激光探测与测量点云数据中的无效区域在三维重建测区影像上对应的范围作为三维重建测区影像的密集匹配区域，根据该三维重建测区影像密集匹配区域构建影像金字塔，从金字塔顶层开始，采用半全局密集匹配方法(Semi-global Matching,SGM)进行由粗到精的分级匹配，获得密集匹配区域的二维视差图，并采用前方交会的方式，生成三维重建测区影像的密集匹配区域所对应的三维点云数据；

步骤3：采用Canny边缘检测算子，在三维重建测区影像上提取线特征，将每个线特征作为一个密集匹配区域，分别进行影像密集匹配，由于线特征只有一维，因此匹配方式采用沿着线特征方向的一维动态规划的优化方式，获取密集匹配区域的视差图，并采用前方交会的方式，生成线特征对应的三维点云数据；

步骤4：在原始激光探测与测量点云数据的基础上，采用点云融合技术将步骤2生成的三维重建测区影像密集匹配区域所对应的三维点云数据，以及步骤3生成的线特征对应的三维点云数据进行三维点云融合，生成融合激光探测与测量点云和密集匹配点云的三维重建点云。

本发明的有益效果：

本发明公开的一种联合激光扫描和影像匹配的三维点云生成方法，能够充分利用影像分辨率高、灰度特征明显的优势，解决LiDAR点云中无反射点区域、线状地物描述等重建问题，能够极大地减少影像密集匹配的计算复杂度，快速获取更加完整精确的点云产品，能够为智慧城市、数字城市等应用服务。

附图说明

图1为本发明的流程图；

图2为激光探测与测量点云中的无效区域；

图3为半全局密集匹配的扫描线方向；

图4为无效区域的密集匹配点云；

图5为激光探测与测量点云对线状地物的定位精度；

图6为线特征匹配的代价积聚方向。

具体实施方式

以下结合附图和具体实施例对本发明作进一步的详细说明：

本发明针对现有的主流三维点云获取技术：激光扫描和影像匹配各自的特点，充分发挥两者各自的优势，提出一种联合激光扫描和影像匹配的三维点云生成方法，能够解决激光扫描在弱/无反射区域、线状地物区域的三维点云获取困难问题，解决影像密集匹配时间复杂度高的问题，生成的三维点云具有完整度高、精度高、点云密度大、获取速度快等优势，如图1所示，它包括如下步骤：

步骤1：激光探测与测量点云完整度评估；

根据原始激光探测与测量点云数据，生成数字表面模型(Digital SurfaceModel，DSM)，提取数字表面模型中的无效区域，将该无效区域中最小外接矩形的四个角点，分别投影到三维重建测区的影像上，生成一个三维重建测区影像上的四边形区域，该四边形区域即可认为是激光探测与测量点云数据中的无效区域在三维重建测区影像上对应的范围；

步骤2：弱/无反射区域三维点云生成；

将激光探测与测量点云数据中的无效区域在三维重建测区影像上对应的范围作为三维重建测区影像的密集匹配区域，根据该三维重建测区影像的密集匹配区域构建影像金字塔，从金字塔顶层开始，采用半全局密集匹配方法进行由粗到精的分级匹配，获得密集匹配区域的二维视差图，并采用前方交会的方式，生成三维重建测区影像密集匹配区域对应的三维点云数据；

步骤3：线特征地物的三维点云生成；

采用Canny边缘检测算子，在三维重建测区影像上提取线特征，将每个线特征作为一个密集匹配区域，分别进行影像密集匹配，由于线特征只有一维，因此匹配方式采用沿着线特征方向的一维动态规划的优化方式，获取密集匹配区域的视差图，并采用前方交会的方式，生成线特征对应的三维点云数据；

步骤4：激光探测与测量点云和密集匹配点云融合；

在原始激光探测与测量点云数据的基础上，采用点云融合技术将步骤2生成的三维重建测区影像密集匹配区域对应的三维点云数据，以及步骤3生成的线特征对应的三维点云数据进行三维点云融合，生成融合激光探测与测量点云和密集匹配点云的三维重建点云。

上述技术方案的步骤1中，本发明的输入是激光探测与测量点云和测区的光学影像集合。为了保证影像密集匹配顺利实施，需要具有重叠度的影像至少两张，即能够组成至少一个立体像对。激光探测与测量是主流的获取地形地貌三维信息的手段之一，得到了广泛的应用。但是在水面等无法反射激光的区域，无法获得激光探测与测量点云，导致激光探测与测量点云中存在无效区域，如图2(a)中的河流区域所示。本文定义仅仅有极少量的，甚至不存在任何三维点的区域为无效区域。因此，需要对激光探测与测量点云进行完整度评估，找出无效区域，在后续处理中，采用影像密集匹配技术，对这些无效区域进行处理，生成三维点云。

上述技术方案的步骤1中，为了方便统计激光探测与测量点云中的无效区域，首先根据激光探测与测量点云，生成数字表面模型(DSM)，定义激光探测与测量点云的点间距为s，一般点间距s由工程文件提供。计算无效区域中最小外接矩形的四个角点的X、Y坐标的方式如下：

X_lb＝min{Xl_i|i＝1…t} Y_lb＝min{Yl_i|i＝1…t}

X_rt＝max{Xl_i|i＝1…t} Y_rt＝max{Yl_i|i＝1…t}

其中，t表示原始激光探测与测量点云数据中三维点的数目；(Xl_i,Yl_i)表示第i个激光探测与测量点的X、Y坐标；(X_lb,Y_lb)表示无效区域中最小外接矩形左下角的角点坐标；(X_rt,Y_rt)表示无效区域中最小外接矩形右上角的角点坐标，(X_lb,Y_rt)表示无效区域中最小外接矩形左上角的角点坐标，(X_rt,Y_lb)表示无效区域中最小外接矩形右下角的角点坐标；

将激光探测与测量点云数据的最小外接矩形范围定义为数字表面模型范围，数字表面模型的起点为无效区域中最小外接矩形左下角的角点，定义数字表面模型的正方形网格的大小为s_D x s_D，s_D＝σxs，其中，s表示激光探测与测量点云数据中点和点之间的平均间距，σ表示数字表面模型的网格大小与上述激光探测与测量点云数据中点和点之间平均间距之间的比值，σ要取大于1的值，一般取3，目的就是为了保证每个网格内都存在激光探测与测量点，如果不存在，则该网格可能处于激光探测与测量点云无效区域；

则数字表面模型的宽W_D和高H_D分别为：

W_D＝(int)(X_rt-X_lb)/s_D；H_D＝(int)(Y_rt-Y_lb)/s_D

其中，int表示取整操作，从而将数字表面模型定义为一个H_D x W_D大小的规则格网，该规则格网的起点为(X_lb,Y_lb)；如果数字表面模型的网格内不存在任何激光探测与测量点，则该数字表面模型的网格定义为无效区域；否则，该数字表面模型的网格定义为有效区域，在有效区域中，每个数字表面模型网格都存在至少一个激光探测与测量点，每个数字表面模型网格的高程为对应激光探测与测量点的最大高程，如下式所示：

Z(m,n)＝max(Z_i|(int)(X_i-X_lb)/s_D＝n；(int)(Y_i-Y_lb)/s_D＝m)

其中，(int)(X_i-X_lb)/s_D＝n；(int)(Y_i-Y_lb)/s_D＝m表示坐标为(X_i,Y_i)的激光探测与测量点落入数字表面模型第m行第n列的网格中；Z_i表示第i个激光探测与测量点的Z坐标；Z(m,n)表示数字表面模型网格(m,n)的高程，取落入网格中所有激光探测与测量点的最大高程；

在有效区域，每个DSM网格都存在有效的高程值；但是由于无效区域不存在LiDAR点，无效区域的网格同样也不存在有效高程。因此，可以根据DSM网格内是否存在有效高程，定义该区域是否为无效区域，如图2(b)中河流区域所示。根据所述无效区域内每个数字表面模型网格的X、Y坐标，定义无效区域的最小外接矩形，如下式所示：

N_lb＝min{n_i|i∈Ω} M_lb＝min{m_i|i∈Ω}

N_rt＝max{n_i|i∈Ω} M_rt＝max{m_i|i∈Ω}

其中，N_lb表示数字表面模型网格所有n列中的最小值，M_lb表示数字表面模型网格所有n列中的最大值，N_rt表示数字表面模型网格所有m行中的最小值，M_rt表示数字表面模型网格所有m行中的最大值，N_lb、M_lb、N_rt和M_rt共同定义了无效区域最小外接矩形的范围；Ω表示无效区域中数字表面模型网格的集合；n表示数字表面模型的网格的列数；m表示数字表面模型的网格的行数；具体外接矩形如图2(b)虚线框所示。

根据无效区域外接矩形四个角点的数字表面模型行列坐标，获得四个角点的物方空间三维坐标，如下式所示：

X_c＝X_lb+n x s_D Y_c＝Y_lb+m x s_D Z_c＝Z(m,n)

其中，(m,n)表示无效区域外接矩形角点的数字表面模型行列坐标；X_c、Y_c、Z_c分别表示无效区域外接矩形角点的物方三维坐标；

最后，通过共线方程，将无效区域外接矩形角点反投影到三维重建测区影像上得到无效区域在三维重建测区影像上的范围，生成无效区域的影像，即三维重建测区影像密集匹配区域，如下式所示：

其中，x_c、y_c表示无效区域外接矩形角点在三维重建测区影像上的坐标；x₀、y₀表示相机内像主点(相机摄影中心与像平面的垂线与像平面的交点)坐标参数，f表示相机的焦距参数；Xs、Ys、Zs表示相机的外方位线元素；a₁、a₂、a₃、b₁、b₂、b₃、c₁、c₂、c₃表示摄影坐标系与大地坐标系之间的旋转矩阵的九个元素，最终无效区域在影像上的显示范围如图2(c)所示。上述旋转矩阵由外方位元素计算得来的旋转矩阵，一般通过测量平差手段获得精确值。

上述技术方案的步骤2中，影像划分为有效区域和无效区域。有效区域已经存在对应的激光探测与测量点云，不需要影像匹配的方法重新生成三维点云；而无效区域不存在对应的激光探测与测量点云，需要借助影像密集匹配的方法来生成无效区域的三维点云。无效区域往往是水面等纹理贫乏区域，为了减少匹配歧义，提高匹配鲁棒性，需要对无效区域影像构建金字塔，本发明采用均值滤波的方式构建影像金字塔，将三维重建测区影像的密集匹配区域构建影像金字塔的方法为：定义S为影像金字塔重采样的尺度，上一级影像金字塔的一个像素，会对应下一级影像金字塔S×S的区域，影像金字塔构建方式如下式所示；

式中，gⁿ表示上一级第n层影像金字塔的像素灰度；

表示下一级第n-1层影像金字塔第i个像素的灰度，S²＝S×S；

在构建影像金字塔后，采用由粗到精的分级匹配策略，获得密集匹配区域的二维视差图，在影像金字塔顶层开始，采用半全局的密集匹配策略，每个像素受到八个方向扫描线的代价累积，如图3所示，在每条扫描线上，代价积聚方式为：

式中，L_r(p,d)表示像素p在当前扫描线对应视差d的累积代价；L_r(p-1,d)表示像素p-1在当前扫描线对应视差d的累积代价；L_r(p-1,d-1)表示像素p-1在当前扫描线对应视差d-1的累积代价；L_r(p-1,d+1)表示像素p-1在当前扫描线对应视差d+1的累积代价；m_kinL_r(p-1,k)表示像素p-1在当前扫描线对应视差k的累积代价中，最小的累积代价，其中k的绝对值大于1；

表示像素p-1在当前扫描线对应视差i的累积代价中，最小的累积代价；r表示扫描线的方向；C(p,d)表示当前像素p对应视差d的代价；p–1表示在当前路径上，像素p的前一个像素；P₁表示视差平滑项惩罚因子；P₂表示视差阶跃项惩罚因子；

单条扫描线的累积结果是不稳定的，容易产生“条纹”问题。为了使匹配结果更加稳定可靠，对无效区域影像进行8个方向的代价累积，最后对各个方向的累积结果相加，得到最终的代价积聚结果，如下式所示：

式中，S(p,d)表示对各个方向代价积聚结果相加后，得到的总体代价积聚结果，采用Winner Takes All(WTA)策略获取顶层影像的初始视差图，将初始视差图传递到下一级影像金字塔，约束下一级影像金字塔的视差搜索范围，然后再次采用半全局的密集匹配方法，获取视差图，并向下一级影像金字塔传递，直至计算到影像金字塔底层为止，根据视差图，可以快速获得立体影像的同名点坐标，如下式所示：

x_r＝x_l-d y_r＝y_l

式中，(x_l,y_l)、(x_r,y_r)分别表示左右影像的同名点坐标；d表示视差值，根据左右影像的同名点坐标，确定对应的三维点坐标：

B_u＝X_S2-X_S1B_v＝Y_S2-Y_S1B_w＝Z_S2-Z_S1

X＝X_S1+U₁＝X_S2+U₂

Y＝Y_S1+V₁＝Y_S2+V₂

Z＝Z_S1+W₁＝Z_S2+W₂

式中，X、Y、Z表示物方点的三维坐标；f表示相机的焦距参数；R_i(i＝1,2)表示旋转矩阵；X_s1、Y_s1、Z_s1、X_s2、Y_s2、Z_s2表示外方位线元素，U₁、V₁、W₁表示地面点在左影像像空间辅助坐标系中的坐标；U₂、V₂、W₂表示地面点在右影像像空间辅助坐标系中的坐标；u₁、v₁、w₁表示像点在左影像像空间辅助坐标系中的坐标；u₂、v₂、w₂表示像点在右影像像空间辅助坐标系中的坐标；B_u、B_v、B_w表示相机之间的基线分量；N₁、N₂表示点投影系数；只要无效区域在影像上是可见的，即可根据影像密集匹配技术恢复无效区域的三维点云，如图4所示。从图4可以看出，河流区域已经充满了密集匹配点云。

上述技术方案的步骤3，由于激光探测与测量点云本身密度较低，直接从LiDAR点云提取的建筑物边缘并不精确，产生边缘“偏移”问题，如图5所示。图5中，上面一条线条表示激光探测与测量点云提取建筑物边缘的投影结果；而下面一条直线表示影像上建筑物的边缘。可以看出，激光探测与测量点云中的建筑物边缘并不精确，存在一定的偏移。为了进一步提高三维点云的精度，有必要借助影像匹配技术，生成建筑物精确的“真实边缘”；

采用canny算子，在三维重建测区影像上提取线特征，将每个线特征作为一个密集匹配区域，分别进行影像密集匹配，具体方法为，沿着线特征的方向，采用一维动态规划方法，进行代价积聚，如下式所示：

式中，L_r(p,d)表示像素p在当前线特征方向上对应视差d的累积代价；r表示线特征的方向；C(p,d)表示当前像素p对应视差d的代价；p–1表示在当前路径上，像素p的前一个像素；

由于线特征一般存在正反两个方向，因此上式所示的代价积聚也存在两个方向，如图6所示，将两个方向的代价积聚结果累加，得到最终的代价积聚结果，采用WinnerTakes All策略获取每条线特征的视差值，计算对应的同名点对，并生成线特征对应的三维点云数据。

上述技术方案的步骤4中，通过影像密集匹配技术，可以生成无效区域和线状地物区域的三维点云。将密集匹配点云和LiDAR点云融合，可以生成更加完整、精度更高的三维点云产品。在无效区域匹配过程中，不可避免地会将一部分有效区域进行重复匹配。由于有效区域已经存在LiDAR点云，因此重复匹配的三维点是不需要融合的，在步骤2中的三维重建测区影像密集匹配区域对应的三维点云数据，以及步骤3中的线特征对应的三维点云数据进行三维点云融合前，需要判断哪些匹配点位于激光探测与测量点云有效区域，哪些匹配点位于激光探测与测量点云无效区域，如下式所示：

m＝(X_i-X_lb)/s_D n＝(Y_i-Y_lb)/s_D

其中，P(X_i,Y_i,Z_i)表示坐标为(X_i,Y_i,Z_i)的三维点；None表示该点不参与融合；Exist表示该点参与融合；Z_DSM表示该点所在数字表面模型网格的高程；(m,n)表示该点所在网格的行列号，valid表示有效值，即数字表面模型网格内存在激光点；invalid表示无效值，即数字表面模型网格内不存在任何激光点；(X_i,Y_i)表示三维点的平面坐标；(X_lb,Y_lb)表示数字表面模型的起点；S_D表示数字表面模型网格的大小。上式表明，当该匹配点落入的数字表面模型网格存在有效高程，则说明该匹配点是有效区域的匹配点，不参与融合过程；否则，说明该匹配点是无效区域的匹配点，参与融合。

本发明充分利用激光扫描技术和影像密集匹配技术各自的优势，生成完整、高精度、稠密的三维点云，在弱\无反射区域，线状地物区域有着可靠精细的三维描述，能够为数字城市、智慧城市等应用提供技术支撑。

本说明书未作详细描述的内容属于本领域专业技术人员公知的现有技术。

Claims (2)

1.一种联合激光扫描和影像匹配的三维点云生成方法，其特征在于，它包括如下步骤：

步骤1：根据原始激光探测与测量点云数据，生成数字表面模型，提取数字表面模型中的无效区域，将该无效区域中最小外接矩形的四个角点，分别投影到三维重建测区的影像上，生成一个三维重建测区影像上的四边形区域，该四边形区域即可认为是激光探测与测量点云数据中的无效区域在三维重建测区影像上对应的范围；

步骤2：将激光探测与测量点云数据中的无效区域在三维重建测区影像上对应的范围作为三维重建测区影像的密集匹配区域，根据该三维重建测区影像的密集匹配区域构建影像金字塔，从金字塔顶层开始，采用半全局密集匹配方法进行由粗到精的分级匹配，获得密集匹配区域的二维视差图，并采用前方交会的方式，生成三维重建测区影像密集匹配区域对应的三维点云数据；

步骤3：采用Canny边缘检测算子，在三维重建测区影像上提取线特征，将每个线特征作为一个密集匹配区域，分别进行影像密集匹配，由于线特征只有一维，因此匹配方式采用沿着线特征方向的一维动态规划的优化方式，获取密集匹配区域的视差图，并采用前方交会的方式，生成线特征对应的三维点云数据；

步骤4：在原始激光探测与测量点云数据的基础上，采用点云融合技术将步骤2生成的三维重建测区影像密集匹配区域对应的三维点云数据，以及步骤3生成的线特征对应的三维点云数据进行三维点云融合，生成融合激光探测与测量点云和密集匹配点云的三维重建点云；

所述步骤4中，在步骤2中的三维重建测区影像密集匹配区域对应的三维点云数据，以及步骤3中的线特征对应的三维点云数据进行三维点云融合前，需要判断哪些匹配点位于激光探测与测量点云有效区域，哪些匹配点位于激光探测与测量点云无效区域；

所述步骤2中，将三维重建测区影像的密集匹配区域构建影像金字塔的方法为：定义S为影像金字塔重采样的尺度，上一级影像金字塔的一个像素，会对应下一级影像金字塔S×S的区域，影像金字塔构建方式如下式所示；

式中，gⁿ表示上一级第n层影像金字塔的像素灰度；

表示下一级第n-1层影像金字塔第i个像素的灰度，S²＝S×S；

在构建影像金字塔后，采用由粗到精的分级匹配策略，获得密集匹配区域的二维视差图，在影像金字塔顶层开始，采用半全局的密集匹配策略，每个像素受到八个方向扫描线的代价累积，在每条扫描线上，代价积聚方式为：

式中，L_r(p,d)表示像素p在当前扫描线对应视差d的累积代价；L_r(p-1,d)表示像素p-1在当前扫描线对应视差d的累积代价；L_r(p-1,d-1)表示像素p-1在当前扫描线对应视差d-1的累积代价；L_r(p-1,d+1)表示像素p-1在当前扫描线对应视差d+1的累积代价；

表示像素p-1在当前扫描线对应视差k的累积代价中，最小的累积代价，其中k的绝对值大于1；

表示像素p-1在当前扫描线对应视差i的累积代价中，最小的累积代价；r表示扫描线的方向；C(p,d)表示当前像素p对应视差d的代价；p–1表示在当前路径上，像素p的前一个像素；P₁表示视差平滑项惩罚因子；P₂表示视差阶跃项惩罚因子；

对无效区域影像进行8个方向的代价累积，最后对各个方向的累积结果相加，得到最终的代价积聚结果，如下式所示：

式中，S(p,d)表示对各个方向代价积聚结果相加后，得到的总体代价积聚结果，采用Winner Takes All策略获取顶层影像的初始视差图，将初始视差图传递到下一级影像金字塔，约束下一级影像金字塔的视差搜索范围，然后再次采用半全局的密集匹配方法，获取视差图，并向下一级影像金字塔传递，直至计算到影像金字塔底层为止，根据视差图，可以快速获得立体影像的同名点坐标，如下式所示：

x_r＝x_l-d y_r＝y_l

式中，(x_l,y_l)、(x_r,y_r)分别表示左右影像的同名点坐标；d表示视差值，根据左右影像的同名点坐标，确定对应的三维点坐标：

B_u＝X_S2-X_S1 B_v＝Y_S2-Y_S1 B_w＝Z_S2-Z_S1

X＝X_S1+U₁＝X_S2+U₂

Y＝Y_S1+V₁＝Y_S2+V₂

Z＝Z_S1+W₁＝Z_S2+W₂

式中，X、Y、Z表示物方点的三维坐标；f表示相机的焦距参数；R_i(i＝1,2)表示旋转矩阵；X_s1、Y_s1、Z_s1、X_s2、Y_s2、Z_s2表示外方位线元素，U₁、V₁、W₁表示地面点在左影像像空间辅助坐标系中的坐标；U₂、V₂、W₂表示地面点在右影像像空间辅助坐标系中的坐标；u₁、v₁、w₁表示像点在左影像像空间辅助坐标系中的坐标；u₂、v₂、w₂表示像点在右影像像空间辅助坐标系中的坐标；B_u、B_v、B_w表示相机之间的基线分量；N₁、N₂表示点投影系数；只要无效区域在影像上是可见的，即可根据影像密集匹配技术恢复无效区域的三维点云；

所述步骤3，采用canny算子，在三维重建测区影像上提取线特征，将每个线特征作为一个密集匹配区域，分别进行影像密集匹配，具体方法为，沿着线特征的方向，采用一维动态规划方法，进行代价积聚，如下式所示：

式中，L_r(p,d)表示像素p在当前线特征方向上对应视差d的累积代价；r表示线特征的方向；C(p,d)表示当前像素p对应视差d的代价；p–1表示在当前路径上，像素p的前一个像素；

由于线特征一般存在正反两个方向，因此上式所示的代价积聚也存在两个方向，将两个方向的代价积聚结果累加，得到最终的代价积聚结果，采用Winner Takes All策略获取每条线特征的视差值，计算对应的同名点对，并生成线特征对应的三维点云数据；

所述步骤4中，在步骤2中的三维重建测区影像密集匹配区域对应的三维点云数据，以及步骤3中的线特征对应的三维点云数据进行三维点云融合前，需要判断哪些匹配点位于激光探测与测量点云有效区域，哪些匹配点位于激光探测与测量点云无效区域，如下式所示：

m＝(X_i-X_lb)/s_D n＝(Y_i-Y_lb)/s_D

其中，P(X_i,Y_i,Z_i)表示坐标为(X_i,Y_i,Z_i)的三维点；None表示该点不参与融合；Exist表示该点参与融合；Z_DSM表示该点所在数字表面模型网格的高程；(m,n)表示该点所在网格的行列号，valid表示有效值，即数字表面模型网格内存在激光点；invalid表示无效值，即数字表面模型网格内不存在任何激光点；(X_i,Y_i)表示三维点的平面坐标；(X_lb,Y_lb)表示数字表面模型的起点；S_D表示数字表面模型网格的大小，上式表明，当该匹配点落入的数字表面模型网格存在有效高程，则说明该匹配点是有效区域的匹配点，不参与融合过程；否则，说明该匹配点是无效区域的匹配点，参与融合。

2.根据权利要求1所述的联合激光扫描和影像匹配的三维点云生成方法，其特征在于：所述步骤1中，计算无效区域中最小外接矩形的四个角点的X、Y坐标的方式如下：

X_lb＝min{Xl_i|i＝1…t} Y_lb＝min{Yl_i|i＝1…t}

X_rt＝max{Xl_i|i＝1…t} Y_rt＝max{Yl_i|i＝1…t}

其中，t表示原始激光探测与测量点云数据中三维点的数目；(Xl_i,Yl_i)表示第i个激光探测与测量点的X、Y坐标；(X_lb,Y_lb)表示无效区域中最小外接矩形左下角的角点坐标；(X_rt,Y_rt)表示无效区域中最小外接矩形右上角的角点坐标，(X_lb,Y_rt)表示无效区域中最小外接矩形左上角的角点坐标，(X_rt,Y_lb)表示无效区域中最小外接矩形右下角的角点坐标；

将激光探测与测量点云数据的最小外接矩形范围定义为数字表面模型范围，数字表面模型的起点为无效区域中最小外接矩形左下角的角点，定义数字表面模型的正方形网格的大小为s_D x s_D，s_D＝σxs，其中，s表示激光探测与测量点云数据中点和点之间的平均间距，σ表示数字表面模型的网格大小与上述激光探测与测量点云数据中点和点之间平均间距之间的比值；

则数字表面模型的宽W_D和高H_D分别为：

W_D＝(int)(X_rt-X_lb)/s_D；H_D＝(int)(Y_rt-Y_lb)/s_D

其中，int表示取整操作，从而将数字表面模型定义为一个H_D x W_D大小的规则格网，该规则格网的起点为(X_lb,Y_lb)；如果数字表面模型的网格内不存在任何激光探测与测量点，则该数字表面模型的网格定义为无效区域；否则，该数字表面模型的网格定义为有效区域，在有效区域中，每个数字表面模型网格都存在至少一个激光探测与测量点，每个数字表面模型网格的高程为对应激光探测与测量点的最大高程，如下式所示：

Z(m,n)＝max(Z_i|(int)(X_i-X_lb)/s_D＝n；(int)(Y_i-Y_lb)/s_D＝m)

其中，(int)(X_i-X_lb)/s_D＝n；(int)(Y_i-Y_lb)/s_D＝m表示坐标为(X_i,Y_i)的激光探测与测量点落入数字表面模型第m行第n列的网格中；Z_i表示第i个激光探测与测量点的Z坐标；Z(m,n)表示数字表面模型网格(m,n)的高程，取落入网格中所有激光探测与测量点的最大高程；

根据所述无效区域内每个数字表面模型网格的X、Y坐标，定义无效区域的最小外接矩形，如下式所示：

N_lb＝min{n_i|i∈Ω} M_lb＝min{m_i|i∈Ω}

N_rt＝max{n_i|i∈Ω} M_rt＝max{m_i|i∈Ω}

其中，N_lb表示数字表面模型网格所有n列中的最小值，M_lb表示数字表面模型网格所有n列中的最大值，N_rt表示数字表面模型网格所有m行中的最小值，M_rt表示数字表面模型网格所有m行中的最大值，N_lb、M_lb、N_rt和M_rt共同定义了无效区域最小外接矩形的范围；Ω表示无效区域中数字表面模型网格的集合；n表示数字表面模型的网格的列数；m表示数字表面模型的网格的行数；

根据无效区域外接矩形四个角点的数字表面模型行列坐标，获得四个角点的物方空间三维坐标，如下式所示：

X_c＝X_lb+n x s_D Y_c＝Y_lb+m x s_D Z_c＝Z(m,n)

其中，(m,n)表示无效区域外接矩形角点的数字表面模型行列坐标；X_c、Y_c、Z_c分别表示无效区域外接矩形角点的物方三维坐标；

最后，通过共线方程，将无效区域外接矩形角点反投影到三维重建测区影像上得到无效区域在三维重建测区影像上的范围，生成无效区域的影像，即三维重建测区影像密集匹配区域，如下式所示：

其中，x_c、y_c表示无效区域外接矩形角点在三维重建测区影像上的坐标；x₀、y₀表示相机内像主点坐标参数，f表示相机的焦距参数；Xs、Ys、Zs表示相机的外方位线元素；a₁、a₂、a₃、b₁、b₂、b₃、c₁、c₂、c₃表示摄影坐标系与大地坐标系之间的旋转矩阵的九个元素。

Images