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

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

Info

Publication number
CN107784666A
CN107784666A CN201710946704.1A CN201710946704A CN107784666A CN 107784666 A CN107784666 A CN 107784666A CN 201710946704 A CN201710946704 A CN 201710946704A CN 107784666 A CN107784666 A CN 107784666A
Authority
CN
China
Prior art keywords
mrow
msubsup
msub
mtd
mtr
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.)
Granted
Application number
CN201710946704.1A
Other languages
English (en)
Other versions
CN107784666B (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 and technology research institute
Original Assignee
Wuhan engineering science and technology research 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 and technology research institute filed Critical Wuhan engineering science and technology research 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

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

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；三个线元素和三个角元素统称为影像的外方位元素；三维点云和影像配准的实质是，在三维点云坐标系下，计算影像的外方位元素；
将特征匹配的三维点作为控制点，将特征匹配的影像点作为对应的像点，代入确定三维点云和影像之间投影关系的后方交会模型：
式中，(vx,vy)表示模型的残差；ΔXs、ΔYs、ΔZs、Δphi、Δomega、Δkappa表示相机外方位元素的改正数；x、y表示真实的特征匹配像点坐标；x0、y0表示将外方位元素初值代入成像模型后，得到的特征匹配像点坐标近似值，其中，成像模型表示为：
式中，(x,y)表示真实的特征匹配像点坐标；f表示相机焦距；(X,Y,Z)表示特征匹配控制点坐标；a1、a2、a3、b1、b2、b3、c1、c2、c3表示由角元素初值计算的旋转矩阵元素；表示相机线元素的初值；
采用最小二乘平差迭代的方式计算相机外方位元素的改正数，在每一次迭代中，外方位元素的计算公式为：
phi＝phi0+Δphi omega＝omega0+Δomega kappa＝kappa0+Δkappa
式中，Xs、Ys、Zs表示相机的外方位线元素，phi、omega、kappa表示相机的外方位角元素，即相机在物方空间的朝向；ΔXs、ΔYs、ΔZs、Δphi、Δomega、Δkappa表示相机外方位元素的改正数；phi0、omega0、kappa0表示相机外方位元素的初值；
步骤2：根据立体影像的像方一致性约束，检测点云中的变化区域，并剔除变化区域中的点云，具体方法如下:
首先，根据所述相机内、外方位元素，将立体影像重采样成核线立体影像；采用成像模型，将所有的三维点投影到核线立体影像上，生成一系列同名像点，成像模型表示为：
式中，(xl,yl)、(xr,yr)表示核线立体影像上的同名像点的两个坐标；f表示相机焦距；(X,Y,Z)表示点云的三维坐标；表示左核线影像的相机外方位线元素；表示右核线影像的相机外方位线元素；表示由左核线影像角元素计算的旋转矩阵元素；表示由右核线影像角元素计算的旋转矩阵元素；
然后，计算核线立体影像上的同名像点的两个坐标(xl,yl)、(xr,yr)之间影像特征相似性的相关系数：
式中，Cov(xl,yl；xr,yr)表示核线立体影像上的同名像点的两个坐标(xl,yl)、(xr,yr)之间影像特征相似性的相关系数；m表示所述相关系数窗口的半径；i,j表示相关系数窗口内的各个像素距离相关系数窗口中心像素的偏移量；g表示左核线影像的灰度值；g'表示右核线影像的灰度值；表示相关系数窗口在左核线影像上的平均灰度值；表示相关系数窗口在右核线影像上的平均灰度值；
将核线立体影像上的同名像点的两个坐标(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；三个线元素和三个角元素统称为影像的外方位元素；三维点云和影像配准的实质是，在三维点云坐标系下，计算影像的外方位元素；
将特征匹配的三维点作为控制点，将特征匹配的影像点作为对应的像点，代入确定三维点云和影像之间投影关系的后方交会模型：
式中，(vx,vy)表示模型的残差；ΔXs、ΔYs、ΔZs、Δphi、Δomega、Δkappa表示相机外方位元素的改正数；x、y表示真实的特征匹配像点坐标；x0、y0表示将外方位元素初值代入成像模型后，得到的特征匹配像点坐标近似值，其中，成像模型表示为：
式中，(x,y)表示真实的特征匹配像点坐标；f表示相机焦距；(X,Y,Z)表示特征匹配控制点坐标；a1、a2、a3、b1、b2、b3、c1、c2、c3表示由角元素初值计算的旋转矩阵元素；表示相机线元素的初值；
相机线元素和角元素的初值，可以由遥感平台搭载的POS系统提供，也可以根据自由网平差手段来提供。其中，自由网平差可以通过开源软件VisualSFM来实现；
根据外方位元素的初值和后方交会模型，采用最小二乘平差迭代的方式计算相机外方位元素的改正数，在每一次迭代中，外方位元素的计算公式为：
phi＝phi0+Δphi omega＝omega0+Δomega kappa＝kappa0+Δkappa
式中，Xs、Ys、Zs表示相机的外方位线元素，phi、omega、kappa表示相机的外方位角元素，即相机在物方空间的朝向；ΔXs、ΔYs、ΔZs、Δphi、Δomega、Δkappa表示相机外方位元素的改正数；phi0、omega0、kappa0表示相机外方位元素的初值。这个初值一般是不准确的，所以我们采用后方交会模型，将初值进行精化，得到准确的外方位元素；
将当前迭代结果，作为下一次迭代的初值，继续代入后方交会模型进行平差处理；一般迭代三次即可得到高精度的相机外方位元素；
步骤2：根据立体影像的像方一致性约束，检测点云中的变化区域，并剔除变化区域中的点云，图3中(a)表示遥感影像，(b)表示对应的三维点云，(c)表示将三维变化区域剔除的结果；比较图3中(a)和(b)可以看出，该区域已经发生三维变化；需要检测变化的区域，并剔除变化区域的三维点云；
具体方法如下:
首先，根据所述相机内、外方位元素，将立体影像重采样成核线立体影像；采用成像模型，将所有的三维点投影到核线立体影像上，生成一系列同名像点，成像模型表示为：
式中，(xl,yl)、(xr,yr)表示核线立体影像上的同名像点的两个坐标；f表示相机焦距；(X,Y,Z)表示点云的三维坐标；表示左核线影像的相机外方位线元素；表示右核线影像的相机外方位线元素；表示由左核线影像角元素计算的旋转矩阵元素；表示由右核线影像角元素计算的旋转矩阵元素；
然后，计算核线立体影像上的同名像点的两个坐标(xl,yl)、(xr,yr)之间影像特征相似性的相关系数：
式中，Cov(xl,yl；xr,yr)表示核线立体影像上的同名像点的两个坐标(xl,yl)、(xr,yr)之间影像特征相似性的相关系数；m表示所述相关系数窗口的半径；i,j表示相关系数窗口内的各个像素距离相关系数窗口中心像素的偏移量；g表示左核线影像的灰度值；g'表示右核线影像的灰度值；表示相关系数窗口在左核线影像上的平均灰度值；表示相关系数窗口在右核线影像上的平均灰度值(相关系数表示同名像点在核线影像上的像方一致性。当相关系数大，说明同名像点在左右核线影像上的灰度特征是一致的，因此对应的三维点没有发生变化；当相关系数小，说明同名像点在左右核线影像上的灰度特征不一样，对应的三维点发生变化)；
将核线立体影像上的同名像点的两个坐标(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来实现，生成变化区域的核线影像同名点；然后根据前方交会法，利用变化区域中的核线影像同名点生成变化区域的三维点云，实现更新点云。
上述技术方案中，前方交会公式为：
式中，(xl,yl)、(xr,yr)表示核线立体影像上的同名像点的两个坐标；X1、Y1、Z1表示变化区域的三维点坐标；f表示相机的焦距参数；表示由左核线影像角元素计算的旋转矩阵元素；表示由右核线影像角元素计算的旋转矩阵元素；表示左核线影像的相机外方位线元素；表示右核线影像的相机外方位线元素；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表示最小外接矩形的角点坐标；
从而将数字表面模型定义为一个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；三个线元素和三个角元素统称为影像的外方位元素；三维点云和影像配准的实质是，在三维点云坐标系下，计算影像的外方位元素；
将特征匹配的三维点作为控制点，将特征匹配的影像点作为对应的像点，代入确定三维点云和影像之间投影关系的后方交会模型：
<mfenced open = "" close = ""> <mtable> <mtr> <mtd> <mrow> <msub> <mi>v</mi> <mi>x</mi> </msub> <mo>=</mo> <mfrac> <mrow> <mo>&amp;part;</mo> <mi>x</mi> </mrow> <mrow> <mo>&amp;part;</mo> <mi>X</mi> <mi>s</mi> </mrow> </mfrac> <mi>&amp;Delta;</mi> <mi>X</mi> <mi>s</mi> <mo>+</mo> <mfrac> <mrow> <mo>&amp;part;</mo> <mi>x</mi> </mrow> <mrow> <mo>&amp;part;</mo> <mi>Y</mi> <mi>s</mi> </mrow> </mfrac> <mi>&amp;Delta;</mi> <mi>Y</mi> <mi>s</mi> <mo>+</mo> <mfrac> <mrow> <mo>&amp;part;</mo> <mi>x</mi> </mrow> <mrow> <mo>&amp;part;</mo> <mi>Z</mi> <mi>s</mi> </mrow> </mfrac> <mi>&amp;Delta;</mi> <mi>Z</mi> <mi>s</mi> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <mo>+</mo> <mfrac> <mrow> <mo>&amp;part;</mo> <mi>x</mi> </mrow> <mrow> <mo>&amp;part;</mo> <mi>p</mi> <mi>h</mi> <mi>i</mi> </mrow> </mfrac> <mi>&amp;Delta;</mi> <mi>p</mi> <mi>h</mi> <mi>i</mi> <mo>+</mo> <mfrac> <mrow> <mo>&amp;part;</mo> <mi>x</mi> </mrow> <mrow> <mo>&amp;part;</mo> <mi>o</mi> <mi>m</mi> <mi>e</mi> <mi>g</mi> <mi>a</mi> </mrow> </mfrac> <mi>&amp;Delta;</mi> <mi>o</mi> <mi>m</mi> <mi>e</mi> <mi>g</mi> <mi>a</mi> <mo>+</mo> <mfrac> <mrow> <mo>&amp;part;</mo> <mi>x</mi> </mrow> <mrow> <mo>&amp;part;</mo> <mi>k</mi> <mi>a</mi> <mi>p</mi> <mi>p</mi> <mi>a</mi> </mrow> </mfrac> <mi>&amp;Delta;</mi> <mi>k</mi> <mi>a</mi> <mi>p</mi> <mi>p</mi> <mi>a</mi> <mo>-</mo> <mrow> <mo>(</mo> <mi>x</mi> <mo>-</mo> <msup> <mi>x</mi> <mn>0</mn> </msup> <mo>)</mo> </mrow> </mrow> </mtd> </mtr> </mtable> </mfenced>
<mfenced open = "" close = ""> <mtable> <mtr> <mtd> <mrow> <msub> <mi>v</mi> <mi>y</mi> </msub> <mo>=</mo> <mfrac> <mrow> <mo>&amp;part;</mo> <mi>y</mi> </mrow> <mrow> <mo>&amp;part;</mo> <mi>X</mi> <mi>s</mi> </mrow> </mfrac> <mi>&amp;Delta;</mi> <mi>X</mi> <mi>s</mi> <mo>+</mo> <mfrac> <mrow> <mo>&amp;part;</mo> <mi>y</mi> </mrow> <mrow> <mo>&amp;part;</mo> <mi>Y</mi> <mi>s</mi> </mrow> </mfrac> <mi>&amp;Delta;</mi> <mi>Y</mi> <mi>s</mi> <mo>+</mo> <mfrac> <mrow> <mo>&amp;part;</mo> <mi>y</mi> </mrow> <mrow> <mo>&amp;part;</mo> <mi>Z</mi> <mi>s</mi> </mrow> </mfrac> <mi>&amp;Delta;</mi> <mi>Z</mi> <mi>s</mi> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <mo>+</mo> <mfrac> <mrow> <mo>&amp;part;</mo> <mi>y</mi> </mrow> <mrow> <mo>&amp;part;</mo> <mi>p</mi> <mi>h</mi> <mi>i</mi> </mrow> </mfrac> <mi>&amp;Delta;</mi> <mi>p</mi> <mi>h</mi> <mi>i</mi> <mo>+</mo> <mfrac> <mrow> <mo>&amp;part;</mo> <mi>y</mi> </mrow> <mrow> <mo>&amp;part;</mo> <mi>o</mi> <mi>m</mi> <mi>e</mi> <mi>g</mi> <mi>a</mi> </mrow> </mfrac> <mi>&amp;Delta;</mi> <mi>o</mi> <mi>m</mi> <mi>e</mi> <mi>g</mi> <mi>a</mi> <mo>+</mo> <mfrac> <mrow> <mo>&amp;part;</mo> <mi>y</mi> </mrow> <mrow> <mo>&amp;part;</mo> <mi>k</mi> <mi>a</mi> <mi>p</mi> <mi>p</mi> <mi>a</mi> </mrow> </mfrac> <mi>&amp;Delta;</mi> <mi>k</mi> <mi>a</mi> <mi>p</mi> <mi>p</mi> <mi>a</mi> <mo>-</mo> <mrow> <mo>(</mo> <mi>y</mi> <mo>-</mo> <msup> <mi>y</mi> <mn>0</mn> </msup> <mo>)</mo> </mrow> </mrow> </mtd> </mtr> </mtable> </mfenced>
式中，(vx,vy)表示模型的残差；ΔXs、ΔYs、ΔZs、Δphi、Δomega、Δkappa表示相机外方位元素的改正数；x、y表示真实的特征匹配像点坐标；x0、y0表示将外方位元素初值代入成像模型后，得到的特征匹配像点坐标近似值，其中，成像模型表示为：
<mrow> <mi>x</mi> <mo>=</mo> <mo>-</mo> <mi>f</mi> <mfrac> <mrow> <msub> <mi>a</mi> <mn>1</mn> </msub> <mrow> <mo>(</mo> <mi>X</mi> <mo>-</mo> <msubsup> <mi>X</mi> <mi>S</mi> <mn>0</mn> </msubsup> <mo>)</mo> </mrow> <mo>+</mo> <msub> <mi>b</mi> <mn>1</mn> </msub> <mrow> <mo>(</mo> <mi>Y</mi> <mo>-</mo> <msubsup> <mi>Y</mi> <mi>S</mi> <mn>0</mn> </msubsup> <mo>)</mo> </mrow> <mo>+</mo> <msub> <mi>c</mi> <mn>1</mn> </msub> <mrow> <mo>(</mo> <mi>Z</mi> <mo>-</mo> <msubsup> <mi>Z</mi> <mi>S</mi> <mn>0</mn> </msubsup> <mo>)</mo> </mrow> </mrow> <mrow> <msub> <mi>a</mi> <mn>3</mn> </msub> <mrow> <mo>(</mo> <mi>X</mi> <mo>-</mo> <msubsup> <mi>X</mi> <mi>S</mi> <mn>0</mn> </msubsup> <mo>)</mo> </mrow> <mo>+</mo> <msub> <mi>b</mi> <mn>3</mn> </msub> <mrow> <mo>(</mo> <mi>Y</mi> <mo>-</mo> <msubsup> <mi>Y</mi> <mi>S</mi> <mn>0</mn> </msubsup> <mo>)</mo> </mrow> <mo>+</mo> <msub> <mi>c</mi> <mn>3</mn> </msub> <mrow> <mo>(</mo> <mi>Z</mi> <mo>-</mo> <msubsup> <mi>Z</mi> <mi>S</mi> <mn>0</mn> </msubsup> <mo>)</mo> </mrow> </mrow> </mfrac> </mrow>
<mrow> <mi>y</mi> <mo>=</mo> <mo>-</mo> <mi>f</mi> <mfrac> <mrow> <msub> <mi>a</mi> <mn>2</mn> </msub> <mrow> <mo>(</mo> <mi>X</mi> <mo>-</mo> <msubsup> <mi>X</mi> <mi>S</mi> <mn>0</mn> </msubsup> <mo>)</mo> </mrow> <mo>+</mo> <msub> <mi>b</mi> <mn>2</mn> </msub> <mrow> <mo>(</mo> <mi>Y</mi> <mo>-</mo> <msubsup> <mi>Y</mi> <mi>S</mi> <mn>0</mn> </msubsup> <mo>)</mo> </mrow> <mo>+</mo> <msub> <mi>c</mi> <mn>2</mn> </msub> <mrow> <mo>(</mo> <mi>Z</mi> <mo>-</mo> <msubsup> <mi>Z</mi> <mi>S</mi> <mn>0</mn> </msubsup> <mo>)</mo> </mrow> </mrow> <mrow> <msub> <mi>a</mi> <mn>3</mn> </msub> <mrow> <mo>(</mo> <mi>X</mi> <mo>-</mo> <msubsup> <mi>X</mi> <mi>S</mi> <mn>0</mn> </msubsup> <mo>)</mo> </mrow> <mo>+</mo> <msub> <mi>b</mi> <mn>3</mn> </msub> <mrow> <mo>(</mo> <mi>Y</mi> <mo>-</mo> <msubsup> <mi>Y</mi> <mi>S</mi> <mn>0</mn> </msubsup> <mo>)</mo> </mrow> <mo>+</mo> <msub> <mi>c</mi> <mn>3</mn> </msub> <mrow> <mo>(</mo> <mi>Z</mi> <mo>-</mo> <msubsup> <mi>Z</mi> <mi>S</mi> <mn>0</mn> </msubsup> <mo>)</mo> </mrow> </mrow> </mfrac> </mrow>
式中，(x,y)表示真实的特征匹配像点坐标；f表示相机焦距；(X,Y,Z)表示特征匹配控制点坐标；a1、a2、a3、b1、b2、b3、c1、c2、c3表示由角元素初值计算的旋转矩阵元素；表示相机线元素的初值；
采用最小二乘平差迭代的方式计算相机外方位元素的改正数，在每一次迭代中，外方位元素的计算公式为：
<mfenced open = "" close = ""> <mtable> <mtr> <mtd> <mrow> <mi>X</mi> <mi>s</mi> <mo>=</mo> <msubsup> <mi>X</mi> <mi>S</mi> <mn>0</mn> </msubsup> <mo>+</mo> <mi>&amp;Delta;</mi> <mi>X</mi> <mi>s</mi> </mrow> </mtd> <mtd> <mrow> <mi>Y</mi> <mi>s</mi> <mo>=</mo> <msubsup> <mi>Y</mi> <mi>S</mi> <mn>0</mn> </msubsup> <mo>+</mo> <mi>&amp;Delta;</mi> <mi>Y</mi> <mi>s</mi> </mrow> </mtd> <mtd> <mrow> <mi>Z</mi> <mi>s</mi> <mo>=</mo> <msubsup> <mi>Z</mi> <mi>S</mi> <mn>0</mn> </msubsup> <mo>+</mo> <mi>&amp;Delta;</mi> <mi>Z</mi> <mi>s</mi> </mrow> </mtd> </mtr> </mtable> </mfenced>
phi＝phi0+Δphi omega＝omega0+Δomega kappa＝kappa0+Δkappa
式中，Xs、Ys、Zs表示相机的外方位线元素，phi、omega、kappa表示相机的外方位角元素，即相机在物方空间的朝向；ΔXs、ΔYs、ΔZs、Δphi、Δomega、Δkappa表示相机外方位元素的改正数；phi0、omega0、kappa0表示相机外方位元素的初值；
步骤2：根据立体影像的像方一致性约束，检测点云中的变化区域，并剔除变化区域中的点云，具体方法如下:
首先，根据所述相机内、外方位元素，将立体影像重采样成核线立体影像；采用成像模型，将所有的三维点投影到核线立体影像上，生成一系列同名像点，成像模型表示为：
<mfenced open = "{" close = ""> <mtable> <mtr> <mtd> <msub> <mi>x</mi> <mi>l</mi> </msub> <mo>=</mo> <mo>-</mo> <mi>f</mi> <mfrac> <mrow> <msubsup> <mi>a</mi> <mn>1</mn> <mi>l</mi> </msubsup> <mrow> <mo>(</mo> <mi>X</mi> <mo>-</mo> <msubsup> <mi>X</mi> <mi>S</mi> <mi>l</mi> </msubsup> <mo>)</mo> </mrow> <mo>+</mo> <msubsup> <mi>b</mi> <mn>1</mn> <mi>l</mi> </msubsup> <mrow> <mo>(</mo> <mi>Y</mi> <mo>-</mo> <msubsup> <mi>Y</mi> <mi>S</mi> <mi>l</mi> </msubsup> <mo>)</mo> </mrow> <mo>+</mo> <msubsup> <mi>c</mi> <mn>1</mn> <mi>l</mi> </msubsup> <mrow> <mo>(</mo> <mi>Z</mi> <mo>-</mo> <msubsup> <mi>Z</mi> <mi>S</mi> <mi>l</mi> </msubsup> <mo>)</mo> </mrow> </mrow> <mrow> <msubsup> <mi>a</mi> <mn>3</mn> <mi>l</mi> </msubsup> <mrow> <mo>(</mo> <mi>X</mi> <mo>-</mo> <msubsup> <mi>X</mi> <mi>S</mi> <mi>l</mi> </msubsup> <mo>)</mo> </mrow> <mo>+</mo> <msubsup> <mi>b</mi> <mn>3</mn> <mi>l</mi> </msubsup> <mrow> <mo>(</mo> <mi>Y</mi> <mo>-</mo> <msubsup> <mi>Y</mi> <mi>S</mi> <mi>l</mi> </msubsup> <mo>)</mo> </mrow> <mo>+</mo> <msubsup> <mi>c</mi> <mn>3</mn> <mi>l</mi> </msubsup> <mrow> <mo>(</mo> <mi>Z</mi> <mo>-</mo> <msubsup> <mi>Z</mi> <mi>S</mi> <mi>l</mi> </msubsup> <mo>)</mo> </mrow> </mrow> </mfrac> </mtd> </mtr> <mtr> <mtd> <msub> <mi>y</mi> <mi>l</mi> </msub> <mo>=</mo> <mo>-</mo> <mi>f</mi> <mfrac> <mrow> <msubsup> <mi>a</mi> <mn>2</mn> <mi>l</mi> </msubsup> <mrow> <mo>(</mo> <mi>X</mi> <mo>-</mo> <msubsup> <mi>X</mi> <mi>S</mi> <mi>l</mi> </msubsup> <mo>)</mo> </mrow> <mo>+</mo> <msubsup> <mi>b</mi> <mn>2</mn> <mi>l</mi> </msubsup> <mrow> <mo>(</mo> <mi>Y</mi> <mo>-</mo> <msubsup> <mi>Y</mi> <mi>S</mi> <mi>l</mi> </msubsup> <mo>)</mo> </mrow> <mo>+</mo> <msubsup> <mi>c</mi> <mn>2</mn> <mi>l</mi> </msubsup> <mrow> <mo>(</mo> <mi>Z</mi> <mo>-</mo> <msubsup> <mi>Z</mi> <mi>S</mi> <mi>l</mi> </msubsup> <mo>)</mo> </mrow> </mrow> <mrow> <msubsup> <mi>a</mi> <mn>3</mn> <mi>l</mi> </msubsup> <mrow> <mo>(</mo> <mi>X</mi> <mo>-</mo> <msubsup> <mi>X</mi> <mi>S</mi> <mi>l</mi> </msubsup> <mo>)</mo> </mrow> <mo>+</mo> <msubsup> <mi>b</mi> <mn>3</mn> <mi>l</mi> </msubsup> <mrow> <mo>(</mo> <mi>Y</mi> <mo>-</mo> <msubsup> <mi>Y</mi> <mi>S</mi> <mi>l</mi> </msubsup> <mo>)</mo> </mrow> <mo>+</mo> <msubsup> <mi>c</mi> <mn>3</mn> <mi>l</mi> </msubsup> <mrow> <mo>(</mo> <mi>Z</mi> <mo>-</mo> <msubsup> <mi>Z</mi> <mi>S</mi> <mi>l</mi> </msubsup> <mo>)</mo> </mrow> </mrow> </mfrac> </mtd> </mtr> </mtable> </mfenced>
<mfenced open = "{" close = ""> <mtable> <mtr> <mtd> <msub> <mi>x</mi> <mi>r</mi> </msub> <mo>=</mo> <mo>-</mo> <mi>f</mi> <mfrac> <mrow> <msubsup> <mi>a</mi> <mn>1</mn> <mi>r</mi> </msubsup> <mrow> <mo>(</mo> <mi>X</mi> <mo>-</mo> <msubsup> <mi>X</mi> <mi>S</mi> <mi>r</mi> </msubsup> <mo>)</mo> </mrow> <mo>+</mo> <msubsup> <mi>b</mi> <mn>1</mn> <mi>r</mi> </msubsup> <mrow> <mo>(</mo> <mi>Y</mi> <mo>-</mo> <msubsup> <mi>Y</mi> <mi>S</mi> <mi>r</mi> </msubsup> <mo>)</mo> </mrow> <mo>+</mo> <msubsup> <mi>c</mi> <mn>1</mn> <mi>r</mi> </msubsup> <mrow> <mo>(</mo> <mi>Z</mi> <mo>-</mo> <msubsup> <mi>Z</mi> <mi>S</mi> <mi>r</mi> </msubsup> <mo>)</mo> </mrow> </mrow> <mrow> <msubsup> <mi>a</mi> <mn>3</mn> <mi>r</mi> </msubsup> <mrow> <mo>(</mo> <mi>X</mi> <mo>-</mo> <msubsup> <mi>X</mi> <mi>S</mi> <mi>r</mi> </msubsup> <mo>)</mo> </mrow> <mo>+</mo> <msubsup> <mi>b</mi> <mn>3</mn> <mi>l</mi> </msubsup> <mrow> <mo>(</mo> <mi>Y</mi> <mo>-</mo> <msubsup> <mi>Y</mi> <mi>S</mi> <mi>r</mi> </msubsup> <mo>)</mo> </mrow> <mo>+</mo> <msubsup> <mi>c</mi> <mn>3</mn> <mi>l</mi> </msubsup> <mrow> <mo>(</mo> <mi>Z</mi> <mo>-</mo> <msubsup> <mi>Z</mi> <mi>S</mi> <mi>r</mi> </msubsup> <mo>)</mo> </mrow> </mrow> </mfrac> </mtd> </mtr> <mtr> <mtd> <msub> <mi>y</mi> <mi>r</mi> </msub> <mo>=</mo> <mo>-</mo> <mi>f</mi> <mfrac> <mrow> <msubsup> <mi>a</mi> <mn>2</mn> <mi>r</mi> </msubsup> <mrow> <mo>(</mo> <mi>X</mi> <mo>-</mo> <msubsup> <mi>X</mi> <mi>S</mi> <mi>r</mi> </msubsup> <mo>)</mo> </mrow> <mo>+</mo> <msubsup> <mi>b</mi> <mn>2</mn> <mi>r</mi> </msubsup> <mrow> <mo>(</mo> <mi>Y</mi> <mo>-</mo> <msubsup> <mi>Y</mi> <mi>S</mi> <mi>r</mi> </msubsup> <mo>)</mo> </mrow> <mo>+</mo> <msubsup> <mi>c</mi> <mn>2</mn> <mi>r</mi> </msubsup> <mrow> <mo>(</mo> <mi>Z</mi> <mo>-</mo> <msubsup> <mi>Z</mi> <mi>S</mi> <mi>r</mi> </msubsup> <mo>)</mo> </mrow> </mrow> <mrow> <msubsup> <mi>a</mi> <mn>3</mn> <mi>r</mi> </msubsup> <mrow> <mo>(</mo> <mi>X</mi> <mo>-</mo> <msubsup> <mi>X</mi> <mi>S</mi> <mi>r</mi> </msubsup> <mo>)</mo> </mrow> <mo>+</mo> <msubsup> <mi>b</mi> <mn>3</mn> <mi>r</mi> </msubsup> <mrow> <mo>(</mo> <mi>Y</mi> <mo>-</mo> <msubsup> <mi>Y</mi> <mi>S</mi> <mi>r</mi> </msubsup> <mo>)</mo> </mrow> <mo>+</mo> <msubsup> <mi>c</mi> <mn>3</mn> <mi>r</mi> </msubsup> <mrow> <mo>(</mo> <mi>Z</mi> <mo>-</mo> <msubsup> <mi>Z</mi> <mi>S</mi> <mi>r</mi> </msubsup> <mo>)</mo> </mrow> </mrow> </mfrac> </mtd> </mtr> </mtable> </mfenced>
式中，(xl,yl)、(xr,yr)表示核线立体影像上的同名像点的两个坐标；f表示相机焦距；(X,Y,Z)表示点云的三维坐标；表示左核线影像的相机外方位线元素；表示右核线影像的相机外方位线元素；表示由左核线影像角元素计算的旋转矩阵元素；表示由右核线影像角元素计算的旋转矩阵元素；
然后，计算核线立体影像上的同名像点的两个坐标(xl,yl)、(xr,yr)之间影像特征相似性的相关系数：
<mrow> <mi>C</mi> <mi>o</mi> <mi>v</mi> <mrow> <mo>(</mo> <msub> <mi>x</mi> <mi>l</mi> </msub> <mo>,</mo> <msub> <mi>y</mi> <mi>l</mi> </msub> <mo>;</mo> <msub> <mi>x</mi> <mi>r</mi> </msub> <mo>,</mo> <msub> <mi>y</mi> <mi>r</mi> </msub> <mo>)</mo> </mrow> <mo>=</mo> <mfrac> <mrow> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>i</mi> <mo>=</mo> <mo>-</mo> <mi>m</mi> </mrow> <mi>m</mi> </munderover> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>j</mi> <mo>=</mo> <mo>-</mo> <mi>m</mi> </mrow> <mi>m</mi> </munderover> <mrow> <mo>(</mo> <mi>g</mi> <mo>(</mo> <mrow> <msub> <mi>x</mi> <mi>l</mi> </msub> <mo>+</mo> <mi>j</mi> <mo>,</mo> <msub> <mi>y</mi> <mi>l</mi> </msub> <mo>+</mo> <mi>i</mi> </mrow> <mo>)</mo> <mo>-</mo> <mover> <mi>g</mi> <mo>&amp;OverBar;</mo> </mover> <mo>)</mo> </mrow> <mo>&amp;CenterDot;</mo> <mrow> <mo>(</mo> <msup> <mi>g</mi> <mo>&amp;prime;</mo> </msup> <mo>(</mo> <mrow> <msub> <mi>x</mi> <mi>r</mi> </msub> <mo>+</mo> <mi>j</mi> <mo>,</mo> <msub> <mi>y</mi> <mi>r</mi> </msub> <mo>+</mo> <mi>i</mi> </mrow> <mo>)</mo> <mo>-</mo> <mover> <msup> <mi>g</mi> <mo>&amp;prime;</mo> </msup> <mo>&amp;OverBar;</mo> </mover> <mo>)</mo> </mrow> </mrow> <msqrt> <mrow> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>i</mi> <mo>=</mo> <mo>-</mo> <mi>m</mi> </mrow> <mi>m</mi> </munderover> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>j</mi> <mo>=</mo> <mo>-</mo> <mi>m</mi> </mrow> <mi>m</mi> </munderover> <msup> <mrow> <mo>(</mo> <mi>g</mi> <mo>(</mo> <mrow> <msub> <mi>x</mi> <mi>l</mi> </msub> <mo>+</mo> <mi>j</mi> <mo>,</mo> <msub> <mi>y</mi> <mi>l</mi> </msub> <mo>+</mo> <mi>i</mi> </mrow> <mo>)</mo> <mo>-</mo> <mover> <mi>g</mi> <mo>&amp;OverBar;</mo> </mover> <mo>)</mo> </mrow> <mn>2</mn> </msup> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>i</mi> <mo>=</mo> <mo>-</mo> <mi>m</mi> </mrow> <mi>m</mi> </munderover> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>j</mi> <mo>=</mo> <mo>-</mo> <mi>m</mi> </mrow> <mi>m</mi> </munderover> <msup> <mrow> <mo>(</mo> <msup> <mi>g</mi> <mo>&amp;prime;</mo> </msup> <mo>(</mo> <mrow> <msub> <mi>x</mi> <mi>r</mi> </msub> <mo>+</mo> <mi>j</mi> <mo>,</mo> <msub> <mi>y</mi> <mi>r</mi> </msub> <mo>+</mo> <mi>i</mi> </mrow> <mo>)</mo> <mo>-</mo> <mover> <msup> <mi>g</mi> <mo>&amp;prime;</mo> </msup> <mo>&amp;OverBar;</mo> </mover> <mo>)</mo> </mrow> <mn>2</mn> </msup> </mrow> </msqrt> </mfrac> </mrow>
式中，Cov(xl,yl；xr,yr)表示核线立体影像上的同名像点的两个坐标(xl,yl)、(xr,yr)之间影像特征相似性的相关系数；m表示所述相关系数窗口的半径；i,j表示相关系数窗口内的各个像素距离相关系数窗口中心像素的偏移量；g表示左核线影像的灰度值；g'表示右核线影像的灰度值；表示相关系数窗口在左核线影像上的平均灰度值；表示相关系数窗口在右核线影像上的平均灰度值；
将核线立体影像上的同名像点的两个坐标(xl,yl)、(xr,yr)之间影像特征相似性的相关系数Cov(xl,yl；xr,yr)与预设的相关系数阈值进行比较，如果核线立体影像上的同名像点的两个坐标(xl,yl)、(xr,yr)之间影像特征相似性的相关系数Cov(xl,yl；xr,yr)小于上述相关系数阈值，则认为该同名像点以及对应的三维点位于变化区域，删该同名的两个像点，以及同名的两个像点对应的三维点；
步骤3：采用立体影像密集匹配算法，重新生成变化区域的三维点云，达到更新点云的目的，具体方法为：
在上述变化区域中采用半全局密集匹配算法，生成变化区域的核线影像同名点；然后根据前方交会法，利用变化区域中的核线影像同名点生成变化区域的三维点云，实现更新点云。
2.根据权利要求1所述的基于立体影像的地形地物三维变化检测和更新方法，其特征在于：前方交会公式为：
<mfenced open = "" close = ""> <mtable> <mtr> <mtd> <mrow> <mfenced open = "[" close = "]"> <mtable> <mtr> <mtd> <msub> <mi>u</mi> <mn>1</mn> </msub> </mtd> </mtr> <mtr> <mtd> <msub> <mi>v</mi> <mn>1</mn> </msub> </mtd> </mtr> <mtr> <mtd> <msub> <mi>w</mi> <mn>1</mn> </msub> </mtd> </mtr> </mtable> </mfenced> <mo>=</mo> <mfenced open = "(" close = ")"> <mtable> <mtr> <mtd> <msubsup> <mi>a</mi> <mn>1</mn> <mi>l</mi> </msubsup> </mtd> <mtd> <msubsup> <mi>a</mi> <mn>2</mn> <mi>l</mi> </msubsup> </mtd> <mtd> <msubsup> <mi>a</mi> <mn>3</mn> <mi>l</mi> </msubsup> </mtd> </mtr> <mtr> <mtd> <msubsup> <mi>b</mi> <mn>1</mn> <mi>l</mi> </msubsup> </mtd> <mtd> <msubsup> <mi>b</mi> <mn>2</mn> <mi>l</mi> </msubsup> </mtd> <mtd> <msubsup> <mi>b</mi> <mn>3</mn> <mi>l</mi> </msubsup> </mtd> </mtr> <mtr> <mtd> <msubsup> <mi>c</mi> <mn>1</mn> <mi>l</mi> </msubsup> </mtd> <mtd> <msubsup> <mi>c</mi> <mn>2</mn> <mi>l</mi> </msubsup> </mtd> <mtd> <msubsup> <mi>c</mi> <mn>3</mn> <mi>l</mi> </msubsup> </mtd> </mtr> </mtable> </mfenced> <mfenced open = "[" close = "]"> <mtable> <mtr> <mtd> <msub> <mi>x</mi> <mi>l</mi> </msub> </mtd> </mtr> <mtr> <mtd> <msub> <mi>y</mi> <mi>l</mi> </msub> </mtd> </mtr> <mtr> <mtd> <mrow> <mo>-</mo> <mi>f</mi> </mrow> </mtd> </mtr> </mtable> </mfenced> </mrow> </mtd> <mtd> <mrow> <mfenced open = "[" close = "]"> <mtable> <mtr> <mtd> <msub> <mi>u</mi> <mn>2</mn> </msub> </mtd> </mtr> <mtr> <mtd> <msub> <mi>v</mi> <mn>2</mn> </msub> </mtd> </mtr> <mtr> <mtd> <msub> <mi>w</mi> <mn>2</mn> </msub> </mtd> </mtr> </mtable> </mfenced> <mo>=</mo> <mfenced open = "(" close = ")"> <mtable> <mtr> <mtd> <msubsup> <mi>a</mi> <mn>1</mn> <mi>r</mi> </msubsup> </mtd> <mtd> <msubsup> <mi>a</mi> <mn>2</mn> <mi>r</mi> </msubsup> </mtd> <mtd> <msubsup> <mi>a</mi> <mn>3</mn> <mi>r</mi> </msubsup> </mtd> </mtr> <mtr> <mtd> <msubsup> <mi>b</mi> <mn>1</mn> <mi>r</mi> </msubsup> </mtd> <mtd> <msubsup> <mi>b</mi> <mn>2</mn> <mi>r</mi> </msubsup> </mtd> <mtd> <msubsup> <mi>b</mi> <mn>3</mn> <mi>r</mi> </msubsup> </mtd> </mtr> <mtr> <mtd> <msubsup> <mi>c</mi> <mn>1</mn> <mi>r</mi> </msubsup> </mtd> <mtd> <msubsup> <mi>c</mi> <mn>2</mn> <mi>r</mi> </msubsup> </mtd> <mtd> <msubsup> <mi>c</mi> <mn>3</mn> <mi>r</mi> </msubsup> </mtd> </mtr> </mtable> </mfenced> <mfenced open = "[" close = "]"> <mtable> <mtr> <mtd> <msub> <mi>x</mi> <mi>r</mi> </msub> </mtd> </mtr> <mtr> <mtd> <msub> <mi>y</mi> <mi>r</mi> </msub> </mtd> </mtr> <mtr> <mtd> <mrow> <mo>-</mo> <mi>f</mi> </mrow> </mtd> </mtr> </mtable> </mfenced> </mrow> </mtd> </mtr> </mtable> </mfenced>
<mfenced open = "" close = ""> <mtable> <mtr> <mtd> <mrow> <msub> <mi>B</mi> <mi>u</mi> </msub> <mo>=</mo> <msubsup> <mi>X</mi> <mi>S</mi> <mi>r</mi> </msubsup> <mo>-</mo> <msubsup> <mi>X</mi> <mi>S</mi> <mi>l</mi> </msubsup> </mrow> </mtd> <mtd> <mrow> <msub> <mi>B</mi> <mi>v</mi> </msub> <mo>=</mo> <msubsup> <mi>Y</mi> <mi>S</mi> <mi>r</mi> </msubsup> <mo>-</mo> <msubsup> <mi>Y</mi> <mi>S</mi> <mi>l</mi> </msubsup> </mrow> </mtd> <mtd> <mrow> <msub> <mi>B</mi> <mi>w</mi> </msub> <mo>=</mo> <msubsup> <mi>Z</mi> <mi>S</mi> <mi>r</mi> </msubsup> <mo>-</mo> <msubsup> <mi>Z</mi> <mi>S</mi> <mi>l</mi> </msubsup> </mrow> </mtd> </mtr> </mtable> </mfenced>
<mfenced open = "" close = ""> <mtable> <mtr> <mtd> <mrow> <msub> <mi>N</mi> <mn>1</mn> </msub> <mo>=</mo> <mfrac> <mrow> <msub> <mi>B</mi> <mi>u</mi> </msub> <msub> <mi>w</mi> <mn>2</mn> </msub> <mo>-</mo> <msub> <mi>B</mi> <mi>w</mi> </msub> <msub> <mi>u</mi> <mn>2</mn> </msub> </mrow> <mrow> <msub> <mi>u</mi> <mn>1</mn> </msub> <msub> <mi>w</mi> <mn>2</mn> </msub> <mo>-</mo> <msub> <mi>u</mi> <mn>2</mn> </msub> <msub> <mi>w</mi> <mn>1</mn> </msub> </mrow> </mfrac> </mrow> </mtd> <mtd> <mrow> <msub> <mi>N</mi> <mn>2</mn> </msub> <mo>=</mo> <mfrac> <mrow> <msub> <mi>B</mi> <mi>u</mi> </msub> <msub> <mi>w</mi> <mn>1</mn> </msub> <mo>-</mo> <msub> <mi>B</mi> <mi>w</mi> </msub> <msub> <mi>u</mi> <mn>1</mn> </msub> </mrow> <mrow> <msub> <mi>u</mi> <mn>1</mn> </msub> <msub> <mi>w</mi> <mn>2</mn> </msub> <mo>-</mo> <msub> <mi>u</mi> <mn>2</mn> </msub> <msub> <mi>w</mi> <mn>1</mn> </msub> </mrow> </mfrac> </mrow> </mtd> </mtr> </mtable> </mfenced>
<mfenced open = "" close = ""> <mtable> <mtr> <mtd> <mrow> <mfenced open = "[" close = "]"> <mtable> <mtr> <mtd> <msub> <mi>U</mi> <mn>1</mn> </msub> </mtd> </mtr> <mtr> <mtd> <msub> <mi>V</mi> <mn>1</mn> </msub> </mtd> </mtr> <mtr> <mtd> <msub> <mi>W</mi> <mn>1</mn> </msub> </mtd> </mtr> </mtable> </mfenced> <mo>=</mo> <msub> <mi>N</mi> <mn>1</mn> </msub> <mfenced open = "[" close = "]"> <mtable> <mtr> <mtd> <msub> <mi>u</mi> <mn>1</mn> </msub> </mtd> </mtr> <mtr> <mtd> <msub> <mi>v</mi> <mn>1</mn> </msub> </mtd> </mtr> <mtr> <mtd> <msub> <mi>w</mi> <mn>1</mn> </msub> </mtd> </mtr> </mtable> </mfenced> </mrow> </mtd> <mtd> <mrow> <mfenced open = "[" close = "]"> <mtable> <mtr> <mtd> <msub> <mi>U</mi> <mn>2</mn> </msub> </mtd> </mtr> <mtr> <mtd> <msub> <mi>V</mi> <mn>2</mn> </msub> </mtd> </mtr> <mtr> <mtd> <msub> <mi>W</mi> <mn>2</mn> </msub> </mtd> </mtr> </mtable> </mfenced> <mo>=</mo> <msub> <mi>N</mi> <mn>2</mn> </msub> <mfenced open = "[" close = "]"> <mtable> <mtr> <mtd> <msub> <mi>u</mi> <mn>2</mn> </msub> </mtd> </mtr> <mtr> <mtd> <msub> <mi>v</mi> <mn>2</mn> </msub> </mtd> </mtr> <mtr> <mtd> <msub> <mi>w</mi> <mn>2</mn> </msub> </mtd> </mtr> </mtable> </mfenced> </mrow> </mtd> </mtr> </mtable> </mfenced>
<mrow> <msub> <mi>X</mi> <mn>1</mn> </msub> <mo>=</mo> <msubsup> <mi>X</mi> <mi>S</mi> <mi>l</mi> </msubsup> <mo>+</mo> <msub> <mi>U</mi> <mn>1</mn> </msub> <mo>=</mo> <msubsup> <mi>X</mi> <mi>S</mi> <mi>r</mi> </msubsup> <mo>+</mo> <msub> <mi>U</mi> <mn>2</mn> </msub> </mrow>
<mrow> <msub> <mi>Y</mi> <mn>1</mn> </msub> <mo>=</mo> <msubsup> <mi>Y</mi> <mi>S</mi> <mi>l</mi> </msubsup> <mo>+</mo> <msub> <mi>V</mi> <mn>1</mn> </msub> <mo>=</mo> <msubsup> <mi>Y</mi> <mi>S</mi> <mi>r</mi> </msubsup> <mo>+</mo> <msub> <mi>V</mi> <mn>2</mn> </msub> </mrow>
<mrow> <msub> <mi>Z</mi> <mn>1</mn> </msub> <mo>=</mo> <msubsup> <mi>Z</mi> <mi>S</mi> <mi>l</mi> </msubsup> <mo>+</mo> <msub> <mi>W</mi> <mn>1</mn> </msub> <mo>=</mo> <msubsup> <mi>Z</mi> <mi>S</mi> <mi>r</mi> </msubsup> <mo>+</mo> <msub> <mi>W</mi> <mn>2</mn> </msub> </mrow>
式中，(xl,yl)、(xr,yr)表示核线立体影像上的同名像点的两个坐标；X1、Y1、Z1表示变化区域的三维点坐标；f表示相机的焦距参数；表示由左核线影像角元素计算的旋转矩阵元素；表示由右核线影像角元素计算的旋转矩阵元素；表示旋转矩阵；表示左核线影像的相机外方位线元素；表示右核线影像的相机外方位线元素；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 true CN107784666A (zh) 2018-03-09
CN107784666B 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)

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN111091594A (zh) * 2019-10-17 2020-05-01 贝壳技术有限公司 多点云平面融合方法及装置
CN113840127A (zh) * 2021-08-12 2021-12-24 长光卫星技术有限公司 一种卫星视频影像获取水域自动掩膜处理dsm的方法

Citations (11)

* 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 东南大学 一种基于扫描激光的大尺度三维地形建模方法
CN102928846A (zh) * 2012-10-24 2013-02-13 华南理工大学 小型无人直升机超低空激光雷达数字地形测绘系统及方法
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 武汉市工程科学技术研究院 联合激光扫描和影像匹配的三维点云生成方法
CN106951860A (zh) * 2017-03-20 2017-07-14 河南腾龙信息工程有限公司 一种基于点云的三维数据智能识别方法
US20170243404A1 (en) * 2016-02-18 2017-08-24 Skycatch, Inc. Generating filtered, three-dimensional digital ground models utilizing multi-stage filters

Patent Citations (11)

* 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 东南大学 一种基于扫描激光的大尺度三维地形建模方法
CN102928846A (zh) * 2012-10-24 2013-02-13 华南理工大学 小型无人直升机超低空激光雷达数字地形测绘系统及方法
CN104049245A (zh) * 2014-06-13 2014-09-17 中原智慧城市设计研究院有限公司 基于LiDAR点云空间差异分析的城市建筑物变化检测方法
CN104063898A (zh) * 2014-06-30 2014-09-24 厦门大学 一种三维点云自动补全方法
US20170243404A1 (en) * 2016-02-18 2017-08-24 Skycatch, Inc. Generating filtered, three-dimensional digital ground models utilizing multi-stage filters
CN105783878A (zh) * 2016-03-11 2016-07-20 三峡大学 一种基于小型无人机遥感的边坡变形检测及量算方法
CN106780712A (zh) * 2016-10-28 2017-05-31 武汉市工程科学技术研究院 联合激光扫描和影像匹配的三维点云生成方法
CN106683173A (zh) * 2016-12-22 2017-05-17 西安电子科技大学 一种基于邻域块匹配提高三维重建点云稠密程度的方法
CN106951860A (zh) * 2017-03-20 2017-07-14 河南腾龙信息工程有限公司 一种基于点云的三维数据智能识别方法

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
王瑞岩: "大规模城市三维重建中点云配准及平面提取研究", 《信号处理》 *

Cited By (4)

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

Also Published As

Publication number Publication date
CN107784666B (zh) 2021-01-08

Similar Documents

Publication Publication Date Title
CN110570466B (zh) 三维实景点云模型的生成方法和装置
KR100912715B1 (ko) 이종 센서 통합 모델링에 의한 수치 사진 측량 방법 및장치
CN102506824B (zh) 一种城市低空无人机系统生成数字正射影像图的方法
US7944547B2 (en) Method and system of generating 3D images with airborne oblique/vertical imagery, GPS/IMU data, and LIDAR elevation data
CN106780712B (zh) 联合激光扫描和影像匹配的三维点云生成方法
Xie et al. Study on construction of 3D building based on UAV images
Haala et al. Extracting 3D urban models from oblique aerial images
Maurer et al. Tapping into the Hexagon spy imagery database: A new automated pipeline for geomorphic change detection
CN112270698B (zh) 基于最邻近曲面的非刚性几何配准方法
CN113358091B (zh) 一种利用三线阵立体卫星影像生产数字高程模型dem的方法
Cosido et al. Hybridization of convergent photogrammetry, computer vision, and artificial intelligence for digital documentation of cultural heritage-a case study: the magdalena palace
CN110889899A (zh) 一种数字地表模型的生成方法及装置
CN112288637A (zh) 无人机航拍图像快速拼接装置及快速拼接方法
CN110986888A (zh) 一种航空摄影一体化方法
CN107784666A (zh) 基于立体影像的地形地物三维变化检测和更新方法
Gao et al. Automatic geo-referencing mobile laser scanning data to UAV images
CN112785686A (zh) 一种基于大数据的林区地图构建方法及可读存储介质
CN116883604A (zh) 一种基于天、空、地影像的三维建模技术方法
Majid et al. Three-dimensional mapping of an ancient cave paintings using close-range photogrammetry and terrestrial laser scanning technologies
Seong et al. UAV Utilization for Efficient Estimation of Earthwork Volume Based on DEM
Previtali et al. An automatic multi-image procedure for accurate 3D object reconstruction
Jun et al. Research on oblique aerial photography and its application in detection of hidden geological hazards
Zhang et al. Tests and performance evaluation of DMC images and new methods for their processing
Gabara et al. Kortowo test field for testing photogrammetric products accuracy–design and first evaluation
Wu et al. Building Facade Reconstruction Using Crowd-Sourced Photos and Two-Dimensional Maps

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