CN107784666B - 基于立体影像的地形地物三维变化检测和更新方法 - Google Patents

基于立体影像的地形地物三维变化检测和更新方法 Download PDF

Info

Publication number
CN107784666B
CN107784666B CN201710946704.1A CN201710946704A CN107784666B CN 107784666 B CN107784666 B CN 107784666B CN 201710946704 A CN201710946704 A CN 201710946704A CN 107784666 B CN107784666 B CN 107784666B
Authority
CN
China
Prior art keywords
image
representing
dimensional
point cloud
points
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
Application number
CN201710946704.1A
Other languages
English (en)
Other versions
CN107784666A (zh
Inventor
周刚
黄旭
胡堃
蔡刚山
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
WUHAN ENGINEERING SCIENCE & TECHNOLOGY INSTITUTE
Original Assignee
WUHAN ENGINEERING SCIENCE & TECHNOLOGY INSTITUTE
Priority date (The priority date 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 date listed.)
Filing date
Publication date
Application filed by WUHAN ENGINEERING SCIENCE & TECHNOLOGY INSTITUTE filed Critical WUHAN ENGINEERING SCIENCE & TECHNOLOGY INSTITUTE
Priority to CN201710946704.1A priority Critical patent/CN107784666B/zh
Publication of CN107784666A publication Critical patent/CN107784666A/zh
Application granted granted Critical
Publication of CN107784666B publication Critical patent/CN107784666B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00—Image analysis
    • G06T7/30—Determination of transform parameters for the alignment of images, i.e. image registration
    • G06T7/33—Determination of transform parameters for the alignment of images, i.e. image registration using feature-based methods
    • G06T7/344—Determination of transform parameters for the alignment of images, i.e. image registration using feature-based methods involving models
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/10—Image acquisition modality
    • G06T2207/10004—Still image; Photographic image
    • G06T2207/10012—Stereo images
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/10—Image acquisition modality
    • G06T2207/10028—Range image; Depth image; 3D point clouds

Landscapes

  • Engineering & Computer Science (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Image Processing (AREA)
  • Length Measuring Devices By Optical Means (AREA)
  • Other Investigation Or Analysis Of Materials By Electrical Means (AREA)

Abstract

本发明涉及一种基于立体影像的地形地物三维变化检测和更新方法，它包括如下步骤：1、根据输入的立体影像和三维点云，进行影像和点云的配准；2、根据立体影像的像方一致性约束，检测点云中的变化区域，并剔除变化区域中的点云；3、采用立体影像密集匹配算法，重新生成变化区域的三维点云，达到更新点云的目的。本发明能够以较低的成本，解决三维地形更新的问题，能够满足大范围地形测绘、智慧城市、智能交通等应用。

Description

基于立体影像的地形地物三维变化检测和更新方法
技术领域
本发明涉及三维变化检测和更新技术领域，具体涉及一种基于立体影像的地形地物三维变化检测和更新方法。
背景技术
目前，大范围自然地表和人工地物的三维建模技术已经得到了极为广泛的应用，如智慧城市、国防建设、智能交通、虚拟电商、文物保护、增强现实等。现有主流的三维建模技术包括激光扫描技术LiDAR(Light Detection And Ranging)和影像密集匹配技术。LiDAR系统由全球定位系统(GPS)、惯性导航系统(IMU)以及激光器(Laser Scanner)共同组成，以主动方式快速获取地面模型及空间点云信息，具有速度快、时效性强、作业范围大等优点。但是LiDAR系统成本昂贵，不利于三维模型的更新。影像密集匹配通过覆盖一个测区的多张立体影像，根据同名光线对对相交的原理，从影像的二维信息恢复整个测区的三维空间信息。与LiDAR激光点云相比，影像密集匹配技术的重建速度较慢，但是建模成本低廉。
随着国内城市化进程的大规模开展，建筑物、道路等人工地物信息发生着日新月异的变化，因此，对智慧城市、智能交通、地形测绘等三维应用提出了实时更新的需求。但是由于LiDAR技术成本昂贵，因此不适合大范围地形的更新。
发明内容
本发明的目的在于提供一种基于立体影像的地形地物三维变化检测和更新方法，该方法首先采用特征匹配的方式，实现影像和三维点云的配准，其次根据立体影像的像方一致性约束，实现三维变化区域的检测和剔除；最后采用立体影像密集匹配方法，对变化区域进行更新。该方法能够以低廉的成本，实现三维变化检测和更新，能够为大范围地形测绘、智慧城市、智能交通等应用服务。
为解决上述技术问题，本发明公开的一种基于立体影像的地形地物三维变化检测和更新方法，其特征在于，它包括如下步骤：
步骤1：根据输入的立体影像数据和三维点云数据，利用以下方式进行影像和点云的配准；
对三维点云数据采用降维处理，将三维点云数据重采样成2.5维的数字表面模型；计算三维点云平面范围内最小外接矩形的四个角点的水平方向坐标X，垂直方向坐标Y，即：
Xlb＝min{Xi|i＝1...t}Ylb＝min{Yi|i＝1...t}
Xrt＝max{Xi|i＝1...t}Yrt＝max{Yi|i＝1...t}
其中，t表示三维点云中三维点的数目；(Xi,Yi)表示第i个三维点的X、Y坐标；(Xlb,Ylb)表示三维点云平面范围内最小外接矩形左下角的角点坐标；(Xrt,Yrt)表示三维点云平面范围内最小外接矩形右上角的角点坐标；(Xlb,Yrt)表示三维点云平面范围内最小外接矩形左上角的角点坐标；(Xrt,Ylb)表示三维点云平面范围内最小外接矩形右下角的角点坐标；min表示取最小值；max表示取最大值；
将三维点云的最小外接矩形范围定义为数字表面模型的范围，数字表面模型的起点坐标为(Xlb,Ylb)，则数字表面模型的宽WD和高HD分别为：
WD＝(int)(Xrt-Xlb)/sD；HD＝(int)(Yrt-Ylb)/sD
其中，sD表示数字表面模型网格的大小，int表示取整操作；WD表示数字表面模型的宽；HD表示数字表面模型的高；Xrt,Xlb,Yrt,Ylb表示最小外接矩形的角点坐标；
从而将数字表面模型定义为一个HDxWD大小的规则格网；该格网的起点坐标为(Xlb,Ylb)；在每个数字表面模型的网格中都存在多个三维点，每个网格的高程为对应三维点的最大高程，如下式所示：
Z(m,n)＝max(Zi|(int)(Xi-Xlb)/sD＝n；(int)(Yi-Ylb)/sD＝m)
其中，Zi表示第i个三维点的高程方向坐标Z；Z(m,n)表示网格(m,n)的高程，取落入网格中所有三维点的最大高程；Xrt,Xlb,Yrt,Ylb表示三维点云平面范围内最小外接矩形的角点坐标；sD表示数字表面模型网格的大小；int表示取整操作；m,n表示数字表面模型中网格的行、列号；max表示取最大值；(int)(Xi-Xlb)/sD＝n；(int)(Yi-Ylb)/sD＝m表示坐标为(Xi,Yi)的三维点落入数字表面模型中第m行第n列的网格中；
然后，采用加速鲁棒特征匹配算子来提取数字表面模型与影像之间的特征匹配点，特征匹配点包括数字表面模型上的三维点和对应影像上的像点，其中数字表面模型上的三维点云和影像上的像点一一对应的关系；
然后，利用三维点云和对应影像上的像点得到确定三维点云和影像之间投影关系的后方交会模型，根据确定三维点云和影像之间投影关系的后方交会模型，实现三维点云和影像的配准；将拍摄上述影像的相机的位置和姿态表示为三个线元素Xs,Ys,Zs，以及三个角元素phi,omega,kappa；三个线元素和三个角元素统称为影像的外方位元素；三维点云和影像配准的实质是，在三维点云坐标系下，计算影像的外方位元素；
将特征匹配的三维点作为控制点，将特征匹配的影像点作为对应的像点，代入确定三维点云和影像之间投影关系的后方交会模型：
Figure BDA0001431853840000031
Figure BDA0001431853840000032
式中，(vx,vy)表示模型的残差；ΔXs、ΔYs、ΔZs、Δphi、Δomega、Δkappa表示相机外方位元素的改正数；x、y表示真实的特征匹配像点坐标；x0、y0表示将外方位元素初值代入成像模型后，得到的特征匹配像点坐标近似值，其中，成像模型表示为：
Figure BDA0001431853840000041
Figure BDA0001431853840000042
式中，(x,y)表示真实的特征匹配像点坐标；f表示相机焦距；(X,Y,Z)表示特征匹配控制点坐标；a1、a2、a3、b1、b2、b3、c1、c2、c3表示由角元素初值计算的旋转矩阵元素；
Figure BDA0001431853840000043
表示相机线元素的初值；
采用最小二乘平差迭代的方式计算相机外方位元素的改正数，在每一次迭代中，外方位元素的计算公式为：
Figure BDA0001431853840000044
phi＝phi0+Δphi omega＝omega0+Δomega kappa＝kappa0+Δkappa
式中，Xs、Ys、Zs表示相机的外方位线元素，phi、omega、kappa表示相机的外方位角元素，即相机在物方空间的朝向；ΔXs、ΔYs、ΔZs、Δphi、Δomega、Δkappa表示相机外方位元素的改正数；
Figure BDA0001431853840000045
phi0、omega0、kappa0表示相机外方位元素的初值；
步骤2：根据立体影像的像方一致性约束，检测点云中的变化区域，并剔除变化区域中的点云，具体方法如下:
首先，根据所述相机内、外方位元素，将立体影像重采样成核线立体影像；采用成像模型，将所有的三维点投影到核线立体影像上，生成一系列同名像点，成像模型表示为：
Figure BDA0001431853840000046
Figure BDA0001431853840000047
式中，(xl,yl)、(xr,yr)表示核线立体影像上的同名像点的两个坐标；f表示相机焦距；(X,Y,Z)表示点云的三维坐标；
Figure BDA0001431853840000051
表示左核线影像的相机外方位线元素；
Figure BDA0001431853840000052
表示右核线影像的相机外方位线元素；
Figure BDA0001431853840000053
表示由左核线影像角元素计算的旋转矩阵元素；
Figure BDA0001431853840000054
表示由右核线影像角元素计算的旋转矩阵元素；
然后，计算核线立体影像上的同名像点的两个坐标(xl,yl)、(xr,yr)之间影像特征相似性的相关系数：
Figure BDA0001431853840000055
式中，Cov(xl,yl；xr,yr)表示核线立体影像上的同名像点的两个坐标(xl,yl)、(xr,yr)之间影像特征相似性的相关系数；m表示所述相关系数窗口的半径；i,j表示相关系数窗口内的各个像素距离相关系数窗口中心像素的偏移量；g表示左核线影像的灰度值；g'表示右核线影像的灰度值；
Figure BDA0001431853840000056
表示相关系数窗口在左核线影像上的平均灰度值；
Figure BDA0001431853840000057
表示相关系数窗口在右核线影像上的平均灰度值；
将核线立体影像上的同名像点的两个坐标(xl,yl)、(xr,yr)之间影像特征相似性的相关系数Cov(xl,yl；xr,yr)与预设的相关系数阈值进行比较，如果核线立体影像上的同名像点的两个坐标(xl,yl)、(xr,yr)之间影像特征相似性的相关系数Cov(xl,yl；xr,yr)小于上述相关系数阈值，则认为该同名像点以及对应的三维点位于变化区域，删该同名的两个像点，以及同名的两个像点对应的三维点；
步骤3：采用立体影像密集匹配算法，重新生成变化区域的三维点云，达到更新点云的目的，具体方法为：
在上述变化区域中采用半全局密集匹配算法，生成变化区域的核线影像同名点；然后根据前方交会法，利用变化区域中的核线影像同名点生成变化区域的三维点云，实现更新点云。
本发明的有益效果：
本发明能够设计一种基于立体影像的地形地物三维变化检测和更新，能够解决传统激光扫描LiDAR的三维检测和更新成本昂贵的问题；在保证三维模型精度的情况下，能够以低廉的成本，实现三维变化检测和更新，最终为地理测绘、智慧城市等应用服务。
附图说明
图1为本发明的流程图；
图2为三维点云和影像之间的配准结果；
图3为三维变化检测结果；
图4为三维变化区域的更新结果；
具体实施方式
以下结合附图和具体实施例对本发明作进一步的详细说明：
本发明的一种基于立体影像的地形地物三维变化检测和更新方法，如图1所示，它包括如下步骤：
步骤1：根据输入的立体影像数据和三维点云数据，利用以下方式进行影像和点云的配准；
三维点云是已有的三维产品，可以根据激光扫描LiDAR获得，也可以根据影像密集匹配获得。立体影像可以根据多类遥感平台获得，包括卫星影像、航拍影像、倾斜影像、无人机影像、近景影像等等。
由于相机安置误差、POS系统误差以及其他系统误差的原因，影像与三维点云之间并不能直接套合，如图2中(a)所示；因此需要实现影像和点云之间的配准；
对三维点云数据采用降维处理，将三维点云数据重采样成2.5维的数字表面模型；计算三维点云平面范围内最小外接矩形的四个角点的水平方向坐标X，垂直方向坐标Y，即：
Xlb＝min{Xi|i＝1...t}Ylb＝min{Yi|i＝1...t}
Xrt＝max{Xi|i＝1...t}Yrt＝max{Yi|i＝1...t}
其中，t表示三维点云中三维点的数目；(Xi,Yi)表示第i个三维点的X、Y坐标；(Xlb,Ylb)表示三维点云平面范围内最小外接矩形左下角的角点坐标；(Xrt,Yrt)表示三维点云平面范围内最小外接矩形右上角的角点坐标；(Xlb,Yrt)表示三维点云平面范围内最小外接矩形左上角的角点坐标；(Xrt,Ylb)表示三维点云平面范围内最小外接矩形右下角的角点坐标；min表示取最小值；max表示取最大值；
将三维点云的最小外接矩形范围定义为数字表面模型的范围，数字表面模型的起点坐标为(Xlb,Ylb)，则数字表面模型的宽WD和高HD分别为：
WD＝(int)(Xrt-Xlb)/sD；HD＝(int)(Yrt-Ylb)/sD
其中，sD表示数字表面模型网格的大小，int表示取整操作；WD表示数字表面模型的宽；HD表示数字表面模型的高；Xrt,Xlb,Yrt,Ylb表示最小外接矩形的角点坐标；
从而将数字表面模型定义为一个HDxWD大小的规则格网；该格网的起点坐标为(Xlb,Ylb)；在每个数字表面模型的网格中都存在多个三维点，每个网格的高程为对应三维点的最大高程，如下式所示：
Z(m,n)＝max(Zi|(int)(Xi-Xlb)/sD＝n；(int)(Yi-Ylb)/sD＝m)
其中，Zi表示第i个三维点的高程方向坐标Z；Z(m,n)表示网格(m,n)的高程，取落入网格中所有三维点的最大高程；Xrt,Xlb,Yrt,Ylb表示三维点云平面范围内最小外接矩形的角点坐标；sD表示数字表面模型网格的大小；int表示取整操作；m,n表示数字表面模型中网格的行、列号；max表示取最大值；(int)(Xi-Xlb)/sD＝n；(int)(Yi-Ylb)/sD＝m表示坐标为(Xi,Yi)的三维点落入数字表面模型中第m行第n列的网格中；
然后，采用加速鲁棒特征匹配算子(Speeded-Up Robust Features，SURF，该算子能够抵抗尺度、旋转等变形的影响，该算子会为数字表面模型上的点，自动寻找影像上对应的像点)来提取数字表面模型与影像之间的特征匹配点，特征匹配点包括数字表面模型上的三维点和对应影像上的像点，其中数字表面模型上的三维点云和影像上的像点一一对应的关系；
在实际工程作业中，本发明采用开源函数库OpenCV中的SurfFeature Detector类来实现Surf特征点的提取和匹配；
然后，利用三维点云和对应影像上的像点得到确定三维点云和影像之间投影关系的后方交会模型，根据确定三维点云和影像之间投影关系的后方交会模型，实现三维点云和影像的配准，如图2中(b)所示；将拍摄上述影像的相机的位置和姿态表示为三个线元素Xs,Ys,Zs，以及三个角元素phi,omega,kappa；三个线元素和三个角元素统称为影像的外方位元素；三维点云和影像配准的实质是，在三维点云坐标系下，计算影像的外方位元素；
将特征匹配的三维点作为控制点，将特征匹配的影像点作为对应的像点，代入确定三维点云和影像之间投影关系的后方交会模型：
Figure BDA0001431853840000081
Figure BDA0001431853840000082
式中，(vx,vy)表示模型的残差；ΔXs、ΔYs、ΔZs、Δphi、Δomega、Δkappa表示相机外方位元素的改正数；x、y表示真实的特征匹配像点坐标；x0、y0表示将外方位元素初值代入成像模型后，得到的特征匹配像点坐标近似值，其中，成像模型表示为：
Figure BDA0001431853840000083
Figure BDA0001431853840000084
式中，(x,y)表示真实的特征匹配像点坐标；f表示相机焦距；(X,Y,Z)表示特征匹配控制点坐标；a1、a2、a3、b1、b2、b3、c1、c2、c3表示由角元素初值计算的旋转矩阵元素；
Figure BDA0001431853840000091
表示相机线元素的初值；
相机线元素和角元素的初值，可以由遥感平台搭载的POS系统提供，也可以根据自由网平差手段来提供。其中，自由网平差可以通过开源软件VisualSFM来实现；
根据外方位元素的初值和后方交会模型，采用最小二乘平差迭代的方式计算相机外方位元素的改正数，在每一次迭代中，外方位元素的计算公式为：
Figure BDA0001431853840000092
phi＝phi0+Δphi omega＝omega0+Δomega kappa＝kappa0+Δkappa
式中，Xs、Ys、Zs表示相机的外方位线元素，phi、omega、kappa表示相机的外方位角元素，即相机在物方空间的朝向；ΔXs、ΔYs、ΔZs、Δphi、Δomega、Δkappa表示相机外方位元素的改正数；
Figure BDA0001431853840000093
phi0、omega0、kappa0表示相机外方位元素的初值。这个初值一般是不准确的，所以我们采用后方交会模型，将初值进行精化，得到准确的外方位元素；
将当前迭代结果，作为下一次迭代的初值，继续代入后方交会模型进行平差处理；一般迭代三次即可得到高精度的相机外方位元素；
步骤2：根据立体影像的像方一致性约束，检测点云中的变化区域，并剔除变化区域中的点云，图3中(a)表示遥感影像，(b)表示对应的三维点云，(c)表示将三维变化区域剔除的结果；比较图3中(a)和(b)可以看出，该区域已经发生三维变化；需要检测变化的区域，并剔除变化区域的三维点云；
具体方法如下:
首先，根据所述相机内、外方位元素，将立体影像重采样成核线立体影像；采用成像模型，将所有的三维点投影到核线立体影像上，生成一系列同名像点，成像模型表示为：
Figure BDA0001431853840000101
Figure BDA0001431853840000102
式中，(xl,yl)、(xr,yr)表示核线立体影像上的同名像点的两个坐标；f表示相机焦距；(X,Y,Z)表示点云的三维坐标；
Figure BDA0001431853840000103
表示左核线影像的相机外方位线元素；
Figure BDA0001431853840000104
表示右核线影像的相机外方位线元素；
Figure BDA0001431853840000105
表示由左核线影像角元素计算的旋转矩阵元素；
Figure BDA0001431853840000106
表示由右核线影像角元素计算的旋转矩阵元素；
然后，计算核线立体影像上的同名像点的两个坐标(xl,yl)、(xr,yr)之间影像特征相似性的相关系数：
Figure BDA0001431853840000107
式中，Cov(xl,yl；xr,yr)表示核线立体影像上的同名像点的两个坐标(xl,yl)、(xr,yr)之间影像特征相似性的相关系数；m表示所述相关系数窗口的半径；i,j表示相关系数窗口内的各个像素距离相关系数窗口中心像素的偏移量；g表示左核线影像的灰度值；g'表示右核线影像的灰度值；
Figure BDA0001431853840000108
表示相关系数窗口在左核线影像上的平均灰度值；
Figure BDA0001431853840000109
表示相关系数窗口在右核线影像上的平均灰度值(相关系数表示同名像点在核线影像上的像方一致性。当相关系数大，说明同名像点在左右核线影像上的灰度特征是一致的，因此对应的三维点没有发生变化；当相关系数小，说明同名像点在左右核线影像上的灰度特征不一样，对应的三维点发生变化)；
将核线立体影像上的同名像点的两个坐标(xl,yl)、(xr,yr)之间影像特征相似性的相关系数Cov(xl,yl；xr,yr)与预设的相关系数阈值进行比较，如果核线立体影像上的同名像点的两个坐标(xl,yl)、(xr,yr)之间影像特征相似性的相关系数Cov(xl,yl；xr,yr)小于上述相关系数阈值，则认为该同名像点以及对应的三维点位于变化区域，删该同名的两个像点，以及同名的两个像点对应的三维点，如图3中(c)所示；
步骤3：采用立体影像密集匹配算法，重新生成变化区域的三维点云，达到更新点云的目的，具体方法为：
在上述变化区域中采用半全局密集匹配算法(Semi-global Matching,SGM)，半全局密集匹配算法SGM可以采用开源函数库OpenCV的类StereoSGBM来实现，生成变化区域的核线影像同名点；然后根据前方交会法，利用变化区域中的核线影像同名点生成变化区域的三维点云，实现更新点云。
上述技术方案中，前方交会公式为：
Figure BDA0001431853840000111
Figure BDA0001431853840000112
Figure BDA0001431853840000113
Figure BDA0001431853840000114
Figure BDA0001431853840000115
Figure BDA0001431853840000116
Figure BDA0001431853840000117
式中，(xl,yl)、(xr,yr)表示核线立体影像上的同名像点的两个坐标；X1、Y1、Z1表示变化区域的三维点坐标；f表示相机的焦距参数；
Figure BDA0001431853840000121
表示由左核线影像角元素计算的旋转矩阵元素；
Figure BDA0001431853840000122
表示由右核线影像角元素计算的旋转矩阵元素；
Figure BDA0001431853840000123
表示左核线影像的相机外方位线元素；
Figure BDA0001431853840000124
表示右核线影像的相机外方位线元素；U1、V1、W1表示地面点在左核线影像像空间辅助坐标系中的坐标；U2、V2、W2表示地面点在右核线影像像空间辅助坐标系中的坐标；u1、v1、w1表示像点在左核线影像像空间辅助坐标系中的坐标；u2、v2、w2表示像点在右影像像空间辅助坐标系中的坐标；Bu、Bv、Bw表示相机之间的基线分量；N1、N2表示点投影系数。
上述技术方案中，所述步骤1中，数字表面模型网格的大小sD与三维点云的平面分辨率一致。
上述技术方案中，所述步骤1中，需要将当前平差的结果作为下一次迭代平差的初值，继续代入后方交会模型进行平差处理；迭代三次即可得到高精度的相机外方位元素。
上述技术方案中，所述步骤2中，相关系数的阈值设为0.8。即若三维点的相关系数小于0.8，则剔除该三维点；若三维点的相关系数大于0.8，则不剔除该三维点；
上述技术方案中，所述相关系数窗口的半径m为6个像素。
本发明设计一种基于立体影像的地形地物三维变化检测和更新方法，能够以低廉的成本，快速实现高精度的三维地形变化检测和更新，能够为智慧城市、智能交通、国土测绘等三维应用提供技术支撑。
本说明书未作详细描述的内容属于本领域专业技术人员公知的现有技术。

Claims (6)

1.一种基于立体影像的地形地物三维变化检测和更新方法，其特征在于，它包括如下步骤：
步骤1：根据输入的立体影像数据和三维点云数据，利用以下方式进行影像和点云的配准；
对三维点云数据采用降维处理，将三维点云数据重采样成2.5维的数字表面模型；计算三维点云平面范围内最小外接矩形的四个角点的水平方向坐标X，垂直方向坐标Y，即：
Xlb＝min{Xi|i＝1...t}Ylb＝min{Yi|i＝1...t}
Xrt＝max{Xi|i＝1...t}Yrt＝max{Yi|i＝1...t}
其中，t表示三维点云中三维点的数目；(Xi,Yi)表示第i个三维点的X、Y坐标；(Xlb,Ylb)表示三维点云平面范围内最小外接矩形左下角的角点坐标；(Xrt,Yrt)表示三维点云平面范围内最小外接矩形右上角的角点坐标；(Xlb,Yrt)表示三维点云平面范围内最小外接矩形左上角的角点坐标；(Xrt,Ylb)表示三维点云平面范围内最小外接矩形右下角的角点坐标；min表示取最小值；max表示取最大值；
将三维点云的最小外接矩形范围定义为数字表面模型的范围，数字表面模型的起点坐标为(Xlb,Ylb)，则数字表面模型的宽WD和高HD分别为：
WD＝(int)(Xrt-Xlb)/sD；HD＝(int)(Yrt-Ylb)/sD
其中，sD表示数字表面模型网格的大小，int表示取整操作；WD表示数字表面模型的宽；HD表示数字表面模型的高；Xrt,Xlb,Yrt,Ylb表示最小外接矩形的角点坐标；
从而将数字表面模型定义为一个HD x WD大小的规则格网；该格网的起点坐标为(Xlb,Ylb)；在每个数字表面模型的网格中都存在多个三维点，每个网格的高程为对应三维点的最大高程，如下式所示：
Z(m,n)＝max(Zi|(int)(Xi-Xlb)/sD＝n；(int)(Yi-Ylb)/sD＝m)
其中，Zi表示第i个三维点的高程方向坐标Z；Z(m,n)表示网格(m,n)的高程，取落入网格中所有三维点的最大高程；Xrt,Xlb,Yrt,Ylb表示三维点云平面范围内最小外接矩形的角点坐标；sD表示数字表面模型网格的大小；int表示取整操作；m,n表示数字表面模型中网格的行、列号；max表示取最大值；(int)(Xi-Xlb)/sD＝n；(int)(Yi-Ylb)/sD＝m表示坐标为(Xi,Yi)的三维点落入数字表面模型中第m行第n列的网格中；
然后，采用加速鲁棒特征匹配算子来提取数字表面模型与影像之间的特征匹配点，特征匹配点包括数字表面模型上的三维点和对应影像上的像点，其中数字表面模型上的三维点云和影像上的像点一一对应的关系；
然后，利用三维点云和对应影像上的像点得到确定三维点云和影像之间投影关系的后方交会模型，根据确定三维点云和影像之间投影关系的后方交会模型，实现三维点云和影像的配准；将拍摄上述影像的相机的位置和姿态表示为三个线元素Xs,Ys,Zs，以及三个角元素phi,omega,kappa；三个线元素和三个角元素统称为影像的外方位元素；三维点云和影像配准的实质是，在三维点云坐标系下，计算影像的外方位元素；
将特征匹配的三维点作为控制点，将特征匹配的影像点作为对应的像点，代入确定三维点云和影像之间投影关系的后方交会模型：
Figure FDA0002788500680000021
Figure FDA0002788500680000022
式中，(vx,vy)表示模型的残差；ΔXs、ΔYs、ΔZs、Δphi、Δomega、Δkappa表示相机外方位元素的改正数；x、y表示真实的特征匹配像点坐标；x0、y0表示将外方位元素初值代入成像模型后，得到的特征匹配像点坐标近似值，其中，成像模型表示为：
Figure FDA0002788500680000031
Figure FDA0002788500680000032
式中，(x,y)表示真实的特征匹配像点坐标；f表示相机焦距；(X,Y,Z)表示特征匹配控制点坐标；a1、a2、a3、b1、b2、b3、c1、c2、c3表示由角元素初值计算的旋转矩阵元素；
Figure FDA0002788500680000033
表示相机线元素的初值；
采用最小二乘平差迭代的方式计算相机外方位元素的改正数，在每一次迭代中，外方位元素的计算公式为：
Figure FDA0002788500680000034
phi＝phi0+Δphi omega＝omega0+Δomega kappa＝kappa0+Δkappa
式中，Xs、Ys、Zs表示相机的外方位线元素，phi、omega、kappa表示相机的外方位角元素，即相机在物方空间的朝向；ΔXs、ΔYs、ΔZs、Δphi、Δomega、Δkappa表示相机外方位元素的改正数；
Figure FDA0002788500680000035
phi0、omega0、kappa0表示相机外方位元素的初值；
步骤2：根据立体影像的像方一致性约束，检测点云中的变化区域，并剔除变化区域中的点云，具体方法如下:
首先，根据相机内、外方位元素，将立体影像重采样成核线立体影像；采用成像模型，将所有的三维点投影到核线立体影像上，生成一系列同名像点，成像模型表示为：
Figure FDA0002788500680000036
Figure FDA0002788500680000037
式中，(xl,yl)、(xr,yr)表示核线立体影像上的同名像点的两个坐标；f表示相机焦距；(X,Y,Z)表示点云的三维坐标；
Figure FDA0002788500680000038
表示左核线影像的相机外方位线元素；
Figure FDA0002788500680000041
表示右核线影像的相机外方位线元素；
Figure FDA0002788500680000042
表示由左核线影像角元素计算的旋转矩阵元素；
Figure FDA0002788500680000043
表示由右核线影像角元素计算的旋转矩阵元素；
然后，计算核线立体影像上的同名像点的两个坐标(xl,yl)、(xr,yr)之间影像特征相似性的相关系数：
Figure FDA0002788500680000044
式中，Cov(xl,yl；xr,yr)表示核线立体影像上的同名像点的两个坐标(xl,yl)、(xr,yr)之间影像特征相似性的相关系数；m表示所述相关系数窗口的半径；i,j表示相关系数窗口内的各个像素距离相关系数窗口中心像素的偏移量；g表示左核线影像的灰度值；g'表示右核线影像的灰度值；
Figure FDA0002788500680000045
表示相关系数窗口在左核线影像上的平均灰度值；
Figure FDA0002788500680000046
表示相关系数窗口在右核线影像上的平均灰度值；
将核线立体影像上的同名像点的两个坐标(xl,yl)、(xr,yr)之间影像特征相似性的相关系数Cov(xl,yl；xr,yr)与预设的相关系数阈值进行比较，如果核线立体影像上的同名像点的两个坐标(xl,yl)、(xr,yr)之间影像特征相似性的相关系数Cov(xl,yl；xr,yr)小于上述相关系数阈值，则认为该同名像点以及对应的三维点位于变化区域，删该同名的两个像点，以及同名的两个像点对应的三维点；
步骤3：采用立体影像密集匹配算法，重新生成变化区域的三维点云，达到更新点云的目的，具体方法为：
在上述变化区域中采用半全局密集匹配算法，生成变化区域的核线影像同名点；然后根据前方交会法，利用变化区域中的核线影像同名点生成变化区域的三维点云，实现更新点云。
2.根据权利要求1所述的基于立体影像的地形地物三维变化检测和更新方法，其特征在于：前方交会公式为：
Figure FDA0002788500680000051
Figure FDA0002788500680000052
Figure FDA0002788500680000053
Figure FDA0002788500680000054
Figure FDA0002788500680000055
Figure FDA0002788500680000056
Figure FDA0002788500680000057
式中，(xl,yl)、(xr,yr)表示核线立体影像上的同名像点的两个坐标；X1、Y1、Z1表示变化区域的三维点坐标；f表示相机的焦距参数；
Figure FDA0002788500680000058
表示由左核线影像角元素计算的旋转矩阵元素；
Figure FDA0002788500680000059
表示由右核线影像角元素计算的旋转矩阵元素；表示旋转矩阵；
Figure FDA00027885006800000510
表示左核线影像的相机外方位线元素；
Figure FDA00027885006800000511
表示右核线影像的相机外方位线元素；U1、V1、W1表示地面点在左核线影像像空间辅助坐标系中的坐标；U2、V2、W2表示地面点在右核线影像像空间辅助坐标系中的坐标；u1、v1、w1表示像点在左核线影像像空间辅助坐标系中的坐标；u2、v2、w2表示像点在右影像像空间辅助坐标系中的坐标；Bu、Bv、Bw表示相机之间的基线分量；N1、N2表示点投影系数。
3.根据权利要求1所述的基于立体影像的地形地物三维变化检测和更新方法，其特征在于：所述步骤1中，数字表面模型网格的大小sD与三维点云的平面分辨率一致。
4.根据权利要求1所述的基于立体影像的地形地物三维变化检测和更新方法，其特征在于：所述步骤1中，需要将当前平差的结果作为下一次迭代平差的初值，继续代入后方交会模型进行平差处理；迭代三次即可得到高精度的相机外方位元素。
5.根据权利要求1所述的基于立体影像的地形地物三维变化检测和更新方法，其特征在于：所述步骤2中，相关系数的阈值设为0.8。
6.根据权利要求1所述的基于立体影像的地形地物三维变化检测和更新方法，其特征在于：所述相关系数窗口的半径m为6个像素。
CN201710946704.1A 2017-10-12 2017-10-12 基于立体影像的地形地物三维变化检测和更新方法 Active CN107784666B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201710946704.1A CN107784666B (zh) 2017-10-12 2017-10-12 基于立体影像的地形地物三维变化检测和更新方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201710946704.1A CN107784666B (zh) 2017-10-12 2017-10-12 基于立体影像的地形地物三维变化检测和更新方法

Publications (2)

Publication Number Publication Date
CN107784666A CN107784666A (zh) 2018-03-09
CN107784666B true CN107784666B (zh) 2021-01-08

Family

ID=61434650

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201710946704.1A Active CN107784666B (zh) 2017-10-12 2017-10-12 基于立体影像的地形地物三维变化检测和更新方法

Country Status (1)

Country Link
CN (1) CN107784666B (zh)

Families Citing this family (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN111091594B (zh) * 2019-10-17 2023-04-11 如你所视(北京)科技有限公司 多点云平面融合方法及装置
CN113840127B (zh) * 2021-08-12 2024-02-27 长光卫星技术股份有限公司 一种卫星视频影像获取水域自动掩膜处理dsm的方法

Citations (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101840517A (zh) * 2010-04-27 2010-09-22 武汉大学 一种基于影像配准的控制点影像库匹配方法及其装置
CN102411778A (zh) * 2011-07-28 2012-04-11 武汉大学 一种机载激光点云与航空影像的自动配准方法
CN102831646A (zh) * 2012-08-13 2012-12-19 东南大学 一种基于扫描激光的大尺度三维地形建模方法
CN104049245A (zh) * 2014-06-13 2014-09-17 中原智慧城市设计研究院有限公司 基于LiDAR点云空间差异分析的城市建筑物变化检测方法
CN104063898A (zh) * 2014-06-30 2014-09-24 厦门大学 一种三维点云自动补全方法
CN105783878A (zh) * 2016-03-11 2016-07-20 三峡大学 一种基于小型无人机遥感的边坡变形检测及量算方法
CN106683173A (zh) * 2016-12-22 2017-05-17 西安电子科技大学 一种基于邻域块匹配提高三维重建点云稠密程度的方法
CN106780712A (zh) * 2016-10-28 2017-05-31 武汉市工程科学技术研究院 联合激光扫描和影像匹配的三维点云生成方法

Family Cites Families (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102928846B (zh) * 2012-10-24 2014-07-02 华南理工大学 小型无人直升机超低空激光雷达数字地形测绘系统及方法
US9846975B2 (en) * 2016-02-18 2017-12-19 Skycatch, Inc. Generating filtered, three-dimensional digital ground models utilizing multi-stage filters
CN106951860A (zh) * 2017-03-20 2017-07-14 河南腾龙信息工程有限公司 一种基于点云的三维数据智能识别方法

Patent Citations (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101840517A (zh) * 2010-04-27 2010-09-22 武汉大学 一种基于影像配准的控制点影像库匹配方法及其装置
CN102411778A (zh) * 2011-07-28 2012-04-11 武汉大学 一种机载激光点云与航空影像的自动配准方法
CN102831646A (zh) * 2012-08-13 2012-12-19 东南大学 一种基于扫描激光的大尺度三维地形建模方法
CN104049245A (zh) * 2014-06-13 2014-09-17 中原智慧城市设计研究院有限公司 基于LiDAR点云空间差异分析的城市建筑物变化检测方法
CN104063898A (zh) * 2014-06-30 2014-09-24 厦门大学 一种三维点云自动补全方法
CN105783878A (zh) * 2016-03-11 2016-07-20 三峡大学 一种基于小型无人机遥感的边坡变形检测及量算方法
CN106780712A (zh) * 2016-10-28 2017-05-31 武汉市工程科学技术研究院 联合激光扫描和影像匹配的三维点云生成方法
CN106683173A (zh) * 2016-12-22 2017-05-17 西安电子科技大学 一种基于邻域块匹配提高三维重建点云稠密程度的方法

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
大规模城市三维重建中点云配准及平面提取研究;王瑞岩;《信号处理》;20150630;第31卷(第6期);第660-668页 *

Also Published As

Publication number Publication date
CN107784666A (zh) 2018-03-09

Similar Documents

Publication Publication Date Title
CN106780712B (zh) 联合激光扫描和影像匹配的三维点云生成方法
Wu et al. Integration of aerial oblique imagery and terrestrial imagery for optimized 3D modeling in urban areas
Bosch et al. A multiple view stereo benchmark for satellite imagery
KR100912715B1 (ko) 이종 센서 통합 모델링에 의한 수치 사진 측량 방법 및장치
Yang et al. Image-based 3D scene reconstruction and exploration in augmented reality
KR101165523B1 (ko) 다중 소스의 지리 정보를 이용한 지리공간 모델링 시스템 및 관련 방법
CN111060924B (zh) 一种slam与目标跟踪方法
CN111486855A (zh) 一种具有物体导航点的室内二维语义栅格地图构建方法
CN109255808B (zh) 基于倾斜影像的建筑物纹理提取方法和装置
GB2506411A (en) Determination of position from images and associated camera positions
Maurer et al. Tapping into the Hexagon spy imagery database: A new automated pipeline for geomorphic change detection
CN112270698B (zh) 基于最邻近曲面的非刚性几何配准方法
CN112184786B (zh) 一种基于合成视觉的目标定位方法
Cosido et al. Hybridization of convergent photogrammetry, computer vision, and artificial intelligence for digital documentation of cultural heritage-a case study: the magdalena palace
CN112862966B (zh) 地表三维模型构建方法、装置、设备及存储介质
CN112288637A (zh) 无人机航拍图像快速拼接装置及快速拼接方法
CN112767461A (zh) 激光点云与序列全景影像自动配准方法
CN113566793A (zh) 基于无人机倾斜影像的真正射影像生成方法和装置
CN114677435A (zh) 一种点云全景融合要素提取方法和系统
CN113947638A (zh) 鱼眼相机影像正射纠正方法
CN116883604A (zh) 一种基于天、空、地影像的三维建模技术方法
Bybee et al. Method for 3-D scene reconstruction using fused LiDAR and imagery from a texel camera
CN107784666B (zh) 基于立体影像的地形地物三维变化检测和更新方法
Zhao et al. Alignment of continuous video onto 3D point clouds
CN112767459A (zh) 基于2d-3d转换的无人机激光点云与序列影像配准方法

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant