-
利用倾斜摄影平台快速获取多角度影像,已成为当前获取空间数据的有效方式之一,尤其是在三维城市建设等领域有着独特的优势。倾斜摄影测量系统所获取的影像不同于传统的航空影像,它存在同视影像基线较短、几何变形大(摄影倾角大)、色差大(逆光、顺光摄影)、海量数据(镜头数量增加、影像重叠度增加)等问题,这些问题给影像匹配、定向与空三等内业处理带来了一系列的挑战。目前,国内对倾斜航空摄影测量技术的研究尚处于早期阶段,在倾斜影像的信息挖掘与行业应用等方面还需进一步地探索与试验。对倾斜影像进行后期处理时,需解决的主要技术有影像预处理、多视倾斜影像匹配、多视倾斜影像空三解算、基于倾斜影像的地形提取、真正射影像构建、基于单张影像的空间信息提取及基于影像的建筑物模型三维重建[1, 2]等。
国内外学者在基于影像进行三维模型重建这一领域进行了大量相关研究[1-3]。
在倾斜影像匹配方面,文献[4]提出了一种具有仿射不变性的倾斜影像匹配方法。文献[5]利用影像初始外方位元素和数字高程模型数据对影像进行几何纠正,利用特征检测(features from accelerated segment test,FAST)和独立二元鲁棒初级特征(binary robust independent elementary features,BRIEF)描述子对纠正影像进行最邻近匹配并利用随机一致性约束算法(random sample consensus,RANSAC)[6]剔除误匹配,再利用三种空间关系约束进一步剔除存在的少量误匹配,最后将匹配点对反算到原始影像上,取得了较好的匹配效果。文献[7]利用影像内方位元素和初始外方位元素对影像进行几何纠正,在纠正影像上提取Harris角点并利用尺度不变特征变换(scale invariant feature transform,SIFT)描述子进行特征描述,将特征点反算到原始影像上,对原始影像上的特征点进行匹配时,采用基本矩阵F、单应矩阵H引导,在局部范围内进行显著性匹配,并利用归一化互相关测度约束剔除少量误匹配,最终得到的匹配点对分布均匀、正确率高。
在影像空三处理方面,文献[8, 9]提出了连续像对相对定向、光束法区域网平差模型和定位定向系统(position orientation system,POS)辅助光束法区域网平差模型。
在影像密集匹配方面,半全局匹配算法(semi-global matching,SGM)[10]及其改进算法tSGM[11]和基于面元的多视立体匹配算法(patch-based multi-view stereo matching,PMVS)[12]由于其良好的匹配性能,应用十分广泛。
在基于密集三维点云构建三维三角网方面,文献[13, 14]提出了一种图割法构网算法,文献[15, 16]提出了一种基于隐函数求解的泊松表面重构算法。在基于影像进行纹理重建方面,文献[2]提出了一种三维纹理重建方法,实现了最优纹理选择和纹理最佳缝合;文献[17, 18]实现了城市建筑物纹理的自动提取及稠密建筑物纹理的遮挡处理。
本文在倾斜影像匹配方面,针对文献[5]存在的两个问题:(1) 由于利用RANSAC算法剔除误匹配时,要保证粗匹配点对中包含至少近50%的正确同名点对[4, 19]。即在“利用FAST检测子和BRIEF描述子对纠正影像进行最邻近匹配,并利用RANSAC方法剔除误匹配”的过程中,对于特殊情况,最邻近匹配中包含的正确同名点对低于50%而导致RANSAC剔除后保留点集为误匹配点集;(2) 对于城市地区大倾角倾斜影像,RANSAC可能剔除了大量正确匹配点对,导致匹配点对分布不均匀,影响空三处理精度[7]。本文采用一种基于透视变换和改进匹配策略的倾斜影像匹配方法得到大量分布均匀的同名点对;针对多视倾斜影像光束法平差方法未知参数多,可能导致平差解算不稳定的问题,推导并实现了带有相对姿态参数的倾斜影像光束法平差模型。
通过一套典型的国产SWDC-5平台[20]获取的倾斜影像数据,对本文提出的倾斜影像自动空中三角测量处理方法进行空三实验,实验结果表明,本文方法空三处理精度较高,验证了方法的可行性和正确性。将本文方法解算得到的空三成果应用到物方面元多视立体匹配、泊松构网和三维纹理快速重建中,计算机可以自动得到城市三维表面模型且效果逼真。
HTML
-
基于倾斜摄影测量的三维城市模型重建的线性结构关系如图 1所示。
1) 在相机标定场,对倾斜相机平台进行检校,得到相机标定结果和相机之间的相对姿态参数。文献[7, 19, 21]对SWDC-5倾斜相机平台和平台参数、同一摄站的影像关系进行了详细介绍。
2) 利用倾斜相机平台获取倾斜影像数据和下视影像的粗略外方位元素。利用相机标定结果对倾斜影像进行畸变校正[22]。
3) 利用下视相机的相机参数、下视影像及下视影像的粗略外方位元素进行空中三角测量,得到下视影像的精确外方位元素。
4) 利用下视影像的精确外方位元素、下视相机与各侧视相机之间的相对姿态参数,计算同一摄站各张侧视影像的粗略外方位元素。
5) 利用倾斜影像数据及各张影像的内、外方位元素进行倾斜影像匹配,得到多视倾斜影像匹配结果。
6) 利用倾斜影像匹配结果及每张影像的内、外方位元素进行POS辅助倾斜影像光束法平差,得到倾斜影像的精确外方位元素。
7) 利用倾斜影像数据及每张影像精确的内、外方位元素进行密集匹配,得到密集三维点云(并计算点云的法向量)。
8) 利用已定向的密集三维点云构建三维三角网,得到三维三角网表面模型(白模)。
9) 利用倾斜影像及其精确的内、外方位元素来确定物方与像方的几何关系,并对三维三角网表面模型进行自动纹理映射,最后得到城市三维表面模型。
-
本文采用一种具有视点不变性的倾斜影像快速匹配算法[7],并对其作了两点改进: (1) 由透视矩阵P计算单应矩阵HP时,文献[7]将地面平均高程值直接取为0,对于特殊情况,比如两摄影点的高程相差较大时,影像纠正得不够准确。因此,本文在多张下视影像上标出一个地面上的同名点,通过前方交会计算出地面的高程值,作为地面平均高程面;(2) 为了进一步提高空三精度,在特征提取阶段,融入了互补不变性特征,对纠正影像提取Harris角点[23, 24]和高斯差分(difference of Gaussian,DoG)特征点[20]。
倾斜影像匹配的基本原理如下。
1) 利用下视影像的外方位元素、下视相机与各侧视相机之间的相对姿态计算同一摄站各张侧视影像的粗略外方位元素。
2) 利用相机参数和倾斜影像的外方位元素计算影像的透视投影矩阵P [6, 7](物方到像方的关系)。
3) 取测区地面平均高程值,令P=[p1 p2 p3 p4],计算物方平面到像方的关系矩阵HP[25]:
4) 对物方平面进行缩放,并统一纠正影像的分辨率,得到近似的正射影像,即纠正影像。所有纠正影像之间不存在尺度和旋转问题。
5) 对纠正影像提取Harris角点[23, 24]和DoG特征点[8, 20],并用SIFT描述子[8, 20]描述这些关键点。将特征点反算到原影像上,并剔除在原影像边缘附近50像素内的特征点。
6) 匹配时,为了使匹配点对分布均匀,提高正确匹配点对数量;并减小搜索区域,提高匹配效率,本文利用基本矩阵F、单应矩阵H引导,仅在局部范围内进行显著性匹配(比值提纯法[8, 20]),即在同名核线附近的10个像素内(由F矩阵的容差
,其中xi为左片中的待匹配点;I′i=F×xi表示待匹配点xi在右片中对应的核线; a′和b′为I′i的前两个元素;计算出右片中对应同名点x′i的取值范围[6]),且在理论同名点附近的200个像素内(按照文献[6]由H矩阵的容差εh=‖x′i-H×xi‖<200,计算出点x′i的取值范围)的区域内进行显著性匹配,再利用归一化互相关测度约束[26]来剔除误匹配,得到匹配点对:式中,ai、bi分别为匹配点对中两个描述向量的元素;a、b分别为两个描述向量中各元素的平均值;γ(a,b)>0.75。
-
平差原理为,首先根据下视影像的外方位元素以及相机之间的相对姿态推算出侧视影像的粗略外方位元素,倾斜影像匹配后,将所有影像参数作为平差参数整体进行光束法平差。平差模型[8, 9]为:
其中,(X,Y,Z)表示物点坐标;(x,y)表示像点坐标;f为相机焦距;Rimu表示惯性测量单元(inertial measurement unit,IMU)坐标到物方空间坐标系之间的旋转矩阵;Rmis表示像空间坐标系到IMU坐标系之间的旋转矩阵;(Xstation,Ystation,Zstation)表示GPS获取的摄站坐标;(X0,Y0,Z0)、(X1,Y1,Z1)表示GPS摄站坐标漂移系统误差改正参数;t表示GPS获取摄站坐标的时间;λ为比例系数。
对于下视影像,Rimu中三个角元素的初值等于下视影像外方位元素中的角元素。对于侧视影像,Rimu中三个角元素的初值等于侧视影像外方位元素中的角元素。对于下视影像,(Xstation,Ystation,Zstation)的初值等于下视影像外方位元素中的线元素;对于侧视影像,(Xstation,Ystation,Zstation)的初值等于侧视影像外方位元素中的线元素。通过前方交会计算(X,Y,Z)的初值,令(X0,Y0,Z0)、(X1,Y1,Z1)及Rmis中三个角元素的初值为0。通过该平差模型的法方程求出加密点物方坐标的改正数以及其他参数的改正数(用于计算影像的外方位元素)。由RimuRmis求影像外方位元素中的角元素,由
求影像外方位元素中的线元素。 -
1) 平差原理。多视倾斜影像光束法平差方法是一种很严密的平差方法。但是在实际的摄影测量任务中,一次性需要处理的影像数量很多,因此该方法的未知参数数量过多,可能导致平差解算不稳定。
如果把“相机的相对姿态在空中飞行任务中为刚体”这一条件考虑在内,将其他相机相对于下视相机的相对姿态作为平差参数参与运算,这样在列误差方程时,同时刻曝光的影像只有下视影像的外方位元素参与平差,大大减少了影像外方位元素的参数。然后根据平差后所得到的下视影像外方位元素和相对姿态来转换同时刻曝光下的侧视影像外方位元素,本文称之为带有相对姿态的倾斜影像光束法平差。
2) 平差模型。垂直航空摄影GPS/IMU联合平差的数学模型[8, 9]同式(2)。针对倾斜相机影像,需要将侧视影像的像空间坐标系转换到下视影像的像空间坐标系,才可利用式(2)。将下视相机的像空间坐标系作为一个基准像空间坐标系,那么,相对姿态参数即为各视相机的像空间坐标系到这个基准坐标系的转换参数。通过相对姿态将各视相机像空间坐标转换到这个基准下的关系式为:
其中,R′rel为相对姿态角元素得到的旋转矩阵;(X′rel,Y′rel,Z′rel)为相对姿态的线元素;(xstation,ystation,zstation)为基准坐标系下的像空间坐标。
将式(3)代入式(2)可得带有相对姿态参数的数学模型:
事实上,本文已经得到带有相对姿态参数的倾斜影像光束法联合平差数学模型,但是在表达侧视相机像空间坐标系与IMU坐标系之间相对关系的时候,出现了两个变换矩阵:Rmis为下视影像像空间坐标系与IMU坐标系之间的变换矩阵;R′rel为侧视相机与下视相机的相对姿态的变换矩阵。这里完全可以用一个相对姿态表达侧视相机像空间坐标系与IMU坐标系之间的关系。那么,侧视影像像空间坐标系到IMU坐标系之间的变换矩阵可以表示为RmisR′rel,相对线元素为Rmis[X′rel Y′rel Z′rel]T。因此,将侧视相机像空间坐标系到IMU坐标系之间的角参数,即mis参数,转换到相对姿态中,用Rrel、Xrel、Yrel、Zrel表示。带有相对姿态参数的数学模型可以重新表示为:
利用Rimu、Rrel、(Xrel,Yrel,Zrel)、(Xstation,Ystation,Zstation)、(X,Y,Z)、(X0,Y0,Z0)、(X1,Y1,Z1)的初值,通过该平差模型的法方程求出加密点(物点)坐标以及其他参数的改正数。这样就可以计算出影像的外方位元素和加密点的物方坐标。
2.1. 倾斜影像匹配
2.2. 倾斜影像光束法平差
2.2.1. 多视倾斜影像光束法平差模型
2.2.2. 带有相对姿态参数的倾斜影像光束法平差模型
-
实验摄影平台为SWDC-5倾斜相机,摄取区域为贵阳市金阳新区(该地区以浅丘地形为主)。对应实验数据的飞行航高为600~1 000 m,航带内间距为150 m,航带间间距为500 m。飞行航线、五相机影像及航带关系示意图如图 2所示。
由于该套倾斜影像数据的像片数量较多,本文选取该套数据中6条航带(20,21,22,23,24,25)的部分影像进行实验,其中下视、前视各13张,后视15张,左视、右视各16张,共73张影像。这73张影像对应的区域为下视影像(13张)的两条航带(航带23和22),对应同时刻曝光的下视影像为65张,如表 1所示。
相机名 影像名 单视影像张数/张 对应同时刻曝光下视影像名 对应航带 A相机 07025AR0031~36;05020AR0031~40 16 07025ER0031~36;05020ER0031~40 25;20 B相机 07023BR0036~42;10022BR0035~42 15 07023ER0036~42;10022ER0035~42 23;22 C相机 07024CR0043~50;05021CR0031~38 16 07024ER0043~50;05021ER0031~38 24;21 D相机 07023DR0023~32;10022DR0027~32 13 07023ER0023~32;10022ER0027~32 23;22 E相机 07023ER0031~37;10022ER0032~37 13 07023ER0031~37;10022ER0032~37 23;22 总张数/张 73 65 Table 1. Relationship Between Experimental Image and Its Bottom-View Image at the Same Camera Station
-
由SWDC-5各相机参数、下视影像外方位元素、下视相机与侧视相机的相对姿态参数,采用一种基于透视变换的改进倾斜影像匹配算法对有重叠区的倾斜影像数据进行匹配,得到连接点文件,图 3(a)是对73张倾斜影像的匹配结果。本文选取由5个不同相机对同一区域拍摄的5张影像,其匹配效果如图 3(b)所示,图 3中显示的匹配点为2度及以上的连接点。由图 3可知,本文改进倾斜影像匹配算法得到的匹配点对分布均匀且较为密集。以某一张下视影像为例,理论最高连接点度数为38度,实际最高连接点度数为22度,10度以上的连接点有696个,5度以上的连接点有7 123个,实际效果如图 4所示。
对于73张实验影像,参与空三的连接点有1 016 575个(像点个数为3 082 798个),点云效果如图 5所示。由图 5可知,参与空三的点云数量密集,分布均匀,这为光束法平差提供了良好的连接点数据。
-
本文使用多视倾斜影像光束法平差模型(设为平差方法1)和带有相对姿态参数的倾斜影像光束法平差模型(设为平差方法2)分别对上述73张影像进行自由网平差实验。
其中,平差方法1将73张倾斜影像的外方位元素作为平差参数进行光束法平差;平差方法2将73张倾斜影像对应同时刻曝光65张下视影像的外方位元素,以及4个倾斜相机相对于下视相机的相对姿态参数作为平差参数进行光束法平差。
由于本文实验数据缺少控制点,因此在比较两种倾斜影像光束法平差方法的精度时,采取了以下办法:分别在下视影像与侧视影像人工刺取相同点(共设置了25个检查点),对下视影像中同名点进行前方交会计算其物方坐标(设为W1),对侧视影像中同名点进行前方交会计算其物方坐标(设为W2),通过坐标差值(W1-W2在三个方向上的分量ΔX、ΔY、ΔZ)和距离差值(W1-W2)来比较自由网精度。两种平差方法的实验结果如表 2所示,表 2中的ΔX、ΔY、ΔZ表示侧视与下视相同点的坐标差值。对于两种倾斜影像光束法平差方法,25个侧视与下视相同点的距离差值比较结果见图 6。
平差方法1 平差方法2 迭代次数/次 3 3 单位权中误差/μm 2.740 505 2.743 400 像点最大残差/μm 11.2 11.8 像点平均残差/μm 1.610 909 1.612 831 ΔX最大值/m 0.518 332 0.523 301 ΔX平均值/m 0.117 331 0.118 542 ΔY最大值/m 0.233 792 0.236 428 ΔY平均值/m 0.094 365 0.101 304 ΔZ最大值/m 0.500 458 0.506 545 ΔZ平均值/m 0.160 64 0.164 913 Table 2. Comparison of Aerotriangulation Accuracy for Two Bundle Adjustment Models (Pixel Size: 6 μm)
Figure 6. Coordinate Differences Between Side-view Image Forward Intersection and Bottom-view Image Forward Intersection for 25 Checkpoints (Blue: the Result of Bundle Adjustment Model 1,Red: the Result of Bundle Adjustment Model 2)
由表 2可知,利用带相对姿态参数的倾斜影像光束法平差方法对该套数据进行空三处理后,空三单位权中误差为2.743 400/6=0.46像素,像点平均残差为1.612 831/6=0.27像素。这一结果表明,针对SWDC-5倾斜相机平台,使用带相对姿态参数的倾斜影像光束法平差方法能够获得较高的空三精度。关联表 2和图 6可知,两种倾斜影像光束法平差方法的自由网平差精度相当且25个检查点的坐标差值和距离差值基本吻合。由于多视倾斜影像光束法平差模型是一种非常严密的倾斜影像光束法平差,也进一步表明了本文提出的带相对姿态参数的倾斜影像光束法平差方法的可行性和稳健性。
实验结果分析如下。
1) 两种倾斜影像平差方法平差参数数量的比较。对于一个大测区,假定摄站数量为N,多视倾斜影像光束法平差方法有6N × 5 = 30N 个外方位未知参数,而带相对姿态参数的倾斜影像光束法平差方法仅有6N + 24个外方位未知参数,平差未知参数的数量大大减少。
2) 精度和适用范围分析。空三实验结果显示:两种倾斜影像平差方法的空三精度都比较高,且25个检查点的坐标差值和距离差值基本吻合,因此当相机之间相对姿态较为稳定且各视相机曝光延迟小的情况下,可以选用带相对姿态参数的光束法平差方法。这样可以减少未知参数,使得解算更加稳定,且可以保证比较高的空三精度,能够满足一般的应用需求。不过,由于存在一定曝光延迟和飞行中气流等的影响,导致SWDC-5侧视相机与下视相机之间的相对姿态参数在飞行中有一定的变化,特别是后视相机B与下视相机E、前视相机D与下视相机E的相对姿态参数变化较大(具体比较数据参见表 3。事实上,相对姿态结算值还包含了带相对姿态平差模型的系统误差),使得带相对姿态参数的倾斜影像光束法平差模型精度略低于多视倾斜影像光束法平差模型精度。需要强调的是,如果各视相机曝光延迟大,或由于实际空中作业时的相对姿态受气压、温度等的影响较大(即五相机的相对姿态在空中飞行时不能看作严格的刚体),那么就需要选用多视倾斜影像光束法平差模型。两种平差模型的优点和适用范围见表 4。
相对姿态参数 A相机-E相机 C相机-E相机 D相机-E相机 B相机-E相机 初始值 差值 初始值 差值 初始值 差值 初始值 差值 X/m 0.146 -0.051 4 -0.111 0.058 5 -0.027 0.242 8 0.045 -0.236 9 Y/m -0.004 -0.003 9 -0.014 0.027 2 -0.121 0.268 2 0.149 0.186 3 Z/m 0.003 -0.022 9 0.037 -0.006 9 0.018 -0.293 8 0.026 -0.333 4 φ/(°) -44.800 2 0.000 0 44.913 33 0.000 0 -1.437 79 -0.004 3 0.230 507 -0.010 7 w/(°) -0.581 43 0.000 0 0.499 547 0.000 0 45.147 42 0.038 7 -44.607 8 -0.033 4 κ/(°) 90.153 45 0.000 0 -89.844 4 0.000 0 1.183 133 0.015 7 -179.615 0.004 9 注:差值=结算值-初始值。 Table 3. Comparison of the Relative Attitude Parameters Between the Two Side- and Bottom- View Cameras Before and After Bundle Adjustment
优点 适用范围 平差方法1 解算精度高,适用范围广 相对姿态变化,各视相机曝光延迟大 平差方法2 未知数少,解算更加稳定 相对姿态稳定,且各视相机曝光延迟小 Table 4. Advantages and Adoption Scopes of the Two Bundle Adjustment Models
-
从73张影像数据中选取15张影像进行三维重建实验。利用平差方法2解算出的空三成果计算每张影像的精确P矩阵[6, 7],再使用PMVS算法[12]进行多视立体匹配,得到密集三维点云如图 7所示。
利用泊松重建方法[15, 16]对已定向的三维密集点云进行构网,得到三维三角网表面模型(白模),其侧视图如图 8所示。
利用精确P矩阵计算物方(即城市三维三角网表面模型)与像方的几何投影关系,获得与地物面对应的影像纹理信息,然后通过遮挡检测技术剔除有遮挡现象的影像纹理信息,最后通过一定的最优算法选择出地物面的理想纹理,实现城市建筑物纹理的自动提取。使用文献[18]中设计的纹理映射方法对三维三角网表面模型自动贴上纹理(对白模的中间部分进行了自动纹理映射),最终的三维表面模型效果如图 9。
该方法无需人工干预,能够快速地得到三维城市模型,其纹理要比手工建模更自然、真实,能够满足一定程度的视觉需求。
3.1. 数据说明
3.2. 倾斜影像自动空三效果
3.2.1. 倾斜影像匹配效果
3.2.2. 倾斜影像空三结果
3.3. 在城市真三维模型重建中的应用
-
本文提出了一种适用于SWDC-5倾斜影像的自动空中三角测量方法,首先设计了一种改进的倾斜影像匹配方法,其次针对多视倾斜影像光束法平差方法未知参数多的问题,提出了带有相对姿态参数的倾斜影像光束法平差模型,并分析了其优点和适用范围。针对SWDC-5倾斜相机平台获取的实验数据,本文提出的倾斜影像自动空中三角测量方法的单位权中误差在0.5像素以内,能够满足一定的生产需求。倾斜影像自动空三成果在城市真三维模型重建中的应用表明,运用本文方法进行城市三维模型重建自动化程度比较高,成本低廉,纹理自然真实,效果逼真,为大规模的城市地区三维表面模型重建提供了良好的参考意义。
下一步的研究重点如下。
1) 改善密集匹配算法的精度和效率。比如,改进PMVS算法的种子点传播扩散过程[27],优化匹配约束条件。对密集点云进行后期处理,包括点云去噪和RANSAC平面拟合等。
2) 提高三维构网的精度。对于泊松表面重建:(1)改善由采样点法向量计算叶子结点法向量的空间插值方法;(2)提高移动立方体算法(marching cubes,MC)[28]由表面重建函数到其三维三角网表面模型的精度,比如,在计算立方体某条边上的等值点坐标与法向量时,用非线性插值方法代替线性插值方法,例如,可以使用二点三次Hermite插值,既保证函数的连续性,也可以保证插值的光滑度;(3)优化和简化三维三角网模型。