一种用于提高山区地貌倾斜摄影测量空三质量的方法

技术领域

本发明涉及一种针对山区地貌倾斜影像数据提高其空三质量的方法，具体涉及一种用于提高山区地貌倾斜摄影测量空三质量的方法。

背景技术

针对山区地貌的倾斜摄影测量数据进行空三，存在一些问题，包括：

1.若山区起伏大，在保证相同地面分辨率的情况下，遥感飞行器需要在不同绝对航高的航线上进行倾斜摄影，容易产生漏洞，使影像重叠率发生变化，不利于空三；在保证相同绝对航高的情况下，所获取的倾斜影像具有不同的地面分辨率，相邻影像上的同一地物大小不一致，不利于影像匹配；

2.山区容易产生方向和速度均不确定的风，使山上植物晃动，导致影像发虚，不利于特征点提取和影像匹配；另外，山中方向和速度均不确定的风会使飞行器产生晃动、倾角变大，甚至可能使飞行器撞山或倾覆，从而导致影像倾角变大、精度变差，不利于空三；

3.山区地貌容易产生阴影，不同时段、不同日期，影像的光影效果不同、颜色不同，不利于影像匹配；

4.山区一般有大面积的植物覆盖，植物属于低度纹理、几何简单的物体，同类植物在影像上近似，不易区分，不利于特征点提取和影像匹配。

因此，对于山区地貌的倾斜摄影测量，较之平地、丘陵，更容易产生空三问题，若空三质量差，后期的实景三维建模可能会出现断层、扭曲、拉伸、模糊、色彩不均匀及漏洞等问题，从而导致作业任务失败。

发明内容

本发明的目的是提供一种用于提高山区地貌倾斜摄影测量空三质量的方法，针对山区的倾斜影像数据进行空中三角测量，定量分析航测数据产品的质量，及时判断、选择最优化结果，解决因山区地貌产生的空三问题，因此，能提高工作效率，减少损失。

为了达到上述目的，本发明有如下技术方案：

本发明的一种用于提高山区地貌倾斜摄影测量空三质量的方法，有以下步骤：

1)根据地面分辨率要求及山区起伏情况确定每条航线的绝对航高；

2)通过航空遥感飞行获取山区的倾斜影像数据、GNSS数据及 IMU数据；

3)对倾斜影像数据、GNSS数据及IMU数据进行整理；

4)整体空三和先分块空三、再合并空三，得出相应的空三结果及空三报告；

5)比较不同条件下的空三结果及空三报告，选择最优结果，用于山区实景三维建模，以形成清晰精确的三维地理信息数据。

其中，步骤1)中每条航线的绝对航高包括若干个相同的高度，或包括若干个不同的高度，本发明既可以处理相同绝对航高的数据，也可以处理不同绝对航高的数据；根据步骤1)中的地面分辨率要求及山区起伏情况，能得到每条航线的相对航高和绝对航高，其中，相对航高与地面分辨率的关系式是：

式中，h是航线的相对航高，f是相机镜头焦距，GSD是地面分辨率，a是相机CCD阵列像元大小，f、a是常数值。

其中，绝对航高与相对航高的关系式是：

H＝h+h′

式中，H是航线的绝对航高，h是航线的相对航高，h′是山区的海拔高度。

若山区起伏小，即山区内部高差小，在相同地面分辨率的要求下，相对航高基本相同，这样所有航线就处于相同的绝对航高；

若山区起伏大，即山区内部高差大，在相同地面分辨率的要求下，航线就处于不同的绝对航高；

反之，若要求所有航线都处于相同的绝对航高，则起伏小的山区对应的数据具有相同的地面分辨率；起伏大的山区对应的数据具有不同的地面分辨率。

其中，所述步骤2)中的倾斜影像数据包括前视、后视、左视、右视及下视五个角度的彩色数码影像数据；所述步骤2)中的GNSS 数据包括经度、纬度及高程的数据，所述步骤2)中的IMU数据包括角元素Roll(Φ)、Pitch(Θ)、Heading(Ψ)；

其中，所述步骤3)中的整理是指：根据遥感飞行的天数，分别整理出每天获取的影像所对应的定位定姿数据以及所有影像所对应的定位定姿数据；

其中，所述步骤4)中的整体空三是指：针对获取的所有数据进行空三；先分块空三、再合并空三是指：分别对每天获取的数据进行空三，然后将各个空三结果合并，再对合并后的数据进行空三；整体空三是在不同时相、不同倾角的条件下进行；分块空三是在相同时相、相同倾角的条件下进行，因此分块空三比整体空三更加稳定，将分块空三合并后再进行空三是基于稳定空三结果的再次空三，其结果也较整体空三稳定；

整体空三和先分块空三、再合并空三包括SURF特征点检测、 SURF特征点描述及RANSAC精确匹配，具体步骤是：

1)SURF特征点检测：利用Hessian矩阵行列式的极大值来检测特征点，假设I是图像，X(x,y)是图像中的一点，尺度是σ，则点X处的Hessian矩阵为：

其中，L_xx(X,σ)是点X处的高斯二阶微分与图像I的卷积，其余各项含义相似；

将近似Hessian矩阵行列式作为X(x,y,σ)处的斑点响应，以加快运算效率，具体公式是：

Det(H_approx)＝D_xxD_yy-(0.9D_xx)²

计算图像中每一个点的斑点响应值，组成响应图像，通过比较某一点在本尺度空间和上下尺度空间的响应值的大小，来判断某一点是否为候选特征点，若响应值比26个邻域值都大或都小，则将该点作为最终的候选特征点，并计算其位置和尺度参数；

2)SURF特征点主方向分配：首先进行Haar小波运算，具体参数是：6s为半径，特征点为中心，边长为4s，获取点在x,y方向上的 Haar小波响应值，其中s为空间尺度；然后进行高斯加权运算，具体参数是：张角为π/3的扇形滑动窗口，步长为0.2弧度滑动窗口，对窗口内图像的Haar小波响应值dx,dy进行累加，获取矢量(m_ω,θ_ω)：

求取Haar小波响应值累加值在多个方向中值最大的一个方向，则将该方向作为特征点的主方向；

3)SURF特征点特征矢量生成：构建以20s为边长，特征点为中心，方向与特征点主方向保持一致，大小为4×4的规则子窗口。采用边长为2σ的Haar小波对图像进行处理，以获得x,y方向上的响应值 dx,dy，利用高斯加权计算每一个子窗口的响应值，以获得每一个子窗口的特征矢量：

υ_子窗口＝[∑dx∑dy∑|dx|∑|dy|]

一组描述子特征矢量共包含4×4×4＝64维特征矢量，可获取一个特征点的完整信息：空间尺度、坐标、64维矢量特征；

4)RANSAC算法精确匹配：根据左视图像像特征点匹配右视图像，右视图像特征点匹配左视图像，再进行筛选，如果两次匹配都能够成功，就把进行匹配点对存入新的数组，进行RANSAC模型估计，判断出正确的匹配点对，进行n次迭代计算，以得到最终的匹配点和转换矩阵。

其中，还包括以下步骤：

针对倾斜影像数据、GNSS数据及以角元素Roll(Φ)、Pitch(Θ)、 Heading(Ψ)作为IMU数据初始值的组合导航数据进行空三，分别迭代计算倾斜影像数据的位置信息及姿态信息，剔除残差、粗差，当达到最佳迭代次数时，若空三结果存在航线集与X-Y平面不平行，则不再继续迭代计算，且判定该结果失败；若所有航线集与X-Y平面平行，则不再继续迭代计算，且判定该结果可行；当没有达到最佳迭代次数时，若空三结果存在航线集与X-Y平面不平行，则继续进行迭代计算倾斜影像数据的位置信息及姿态信息，剔除残差、粗差，直到所有航线集与X-Y平面平行为止，且判定该结果可行；将可行的空三结果作为新的初始值继续空三，迭代计算倾斜影像数据的位置信息及姿态信息，剔除残差、粗差，观察每次计算后的连接点数量以及每次计算后航线集与X-Y平面是否平行，当连接点数量达到最多且航线集与 X-Y平面平行时，则选择该次空三结果作为新的初始值继续空三，迭代优化倾斜影像数据的位置信息及姿态信息，剔除残差、粗差，观察每次优化后的连接点数量以及每次计算后航线集与X-Y平面是否平行，当连接点数量达到最多且航线集与X-Y平面平行时，若GNSS数据没有采用RTK测量，则结束优化，得出空三结果及空三报告；若 GNSS数据采用了RTK测量，则选择该次空三结果作为新的初始值继续空三，迭代优化倾斜影像数据的位置信息及姿态信息，剔除残差、粗差，观察每次优化后的连接点数量以及每次计算后航线集与X-Y平面是否平行，当连接点数量达到最多且航线集与X-Y平面平行时，选择该次空三结果作为最终结果，比较两次空三最终结果，选择连接点多且精度高的用于实景三维建模。

其中，多次迭代计算、优化是为了剔除残差、粗差，使像片的位置、姿态正确，同时使连接点数足够多，确保建模所用到的三角形网格数足够多，从而提升三维模型的质量。

由于采取了以上技术方案，本发明的优点在于：

1本发明操作简单，可以使用不同的策略进行计算，然后从不同的结果中找到最优解，使成果质量最佳；

2方便快捷，可以定量分析，精度高，能提高工作效率，减少重复性工作及工作损失；

3能够解决山区地貌引起的问题，避免任务失败。

附图说明

图1为本发明的流程图；

图2为本发明中空三流程图；

图3为本发明中空三结果存在航线集与X-Y平面不平行示意图；

图4为本发明中空三结果正常示意图；

图中，1、航线集一；2、航线集二；

具体实施方式

以下实施例用于说明本发明，但不用来限制本发明的范围。

参见附图1-图4，本发明的一种用于提高山区地貌倾斜摄影测量空三质量的方法，有以下步骤：

1)根据地面分辨率要求及山区起伏情况确定每条航线的绝对航高；

2)通过航空遥感飞行获取山区的倾斜影像数据、GNSS数据及 IMU数据；

3)对倾斜影像数据、GNSS数据及IMU数据进行整理；

4)整体空三和先分块空三、再合并空三，得出相应的空三结果及空三报告；

5)比较不同条件下的空三结果及空三报告，选择最优结果，用于山区实景三维建模，以形成清晰精确的三维地理信息数据。

其中，步骤1)中每条航线的绝对航高包括若干个相同的高度，或包括若干个不同的高度，本发明既可以处理相同绝对航高的数据，也可以处理不同绝对航高的数据；根据步骤1)中的地面分辨率要求及山区起伏情况，能得到每条航线的相对航高和绝对航高，其中，相对航高与地面分辨率的关系式是：

式中，h是航线的相对航高，f是相机镜头焦距，GSD是地面分辨率，a是相机CCD阵列像元大小，f、a是常数值。

其中，绝对航高与相对航高的关系式是：

H＝h+h′

式中，H是航线的绝对航高，h是航线的相对航高，h′是山区的海拔高度。

若山区起伏小，即山区内部高差小，在相同地面分辨率的要求下，相对航高基本相同，这样所有航线就处于相同的绝对航高；

若山区起伏大，即山区内部高差大，在相同地面分辨率的要求下，航线就处于不同的绝对航高；

反之，若要求所有航线都处于相同的绝对航高，则起伏小的山区对应的数据具有相同的地面分辨率；起伏大的山区对应的数据具有不同的地面分辨率。

其中，所述步骤2)中的倾斜影像数据包括前视、后视、左视、右视及下视五个角度的彩色数码影像数据；所述步骤2)中的GNSS 数据包括经度、纬度及高程，所述步骤2)中的IMU数据包括角元素 Roll(Φ)、Pitch(Θ)、Heading(Ψ)；

其中，所述步骤3)中的整理是指：根据遥感飞行的天数，分别整理出每天获取的影像所对应的定位定姿数据以及所有影像所对应的定位定姿数据；

其中，所述步骤4)中的整体空三是指：针对获取的所有数据进行空三；先分块空三、再合并空三是指：分别对每天获取的数据进行空三，然后将各个空三结果合并，再对合并后的数据进行空三；整体空三是在不同时相、不同倾角的条件下进行；分块空三是在相同时相、相同倾角的条件下进行，因此分块空三比整体空三更加稳定，将分块空三合并后再进行空三是基于稳定空三结果的再次空三，其结果也较整体空三稳定；

整体空三和先分块空三、再合并空三包括SURF特征点检测、 SURF特征点描述及RANSAC精确匹配，具体步骤是：

1)SURF特征点检测：利用Hessian矩阵行列式的极大值来检测特征点，假设I是图像，X(x,y)是图像中的一点，尺度是σ，则点X处的Hessian矩阵为：

其中，L_xx(X,σ)是点X处的高斯二阶微分与图像I的卷积，其余各项含义相似；

将近似Hessian矩阵行列式作为X(x,y,σ)处的斑点响应，以加快运算效率，具体公式是：

Det(H_approx)＝D_xxD_yy-(0.9D_xx)²

计算图像中每一个点的斑点响应值，组成响应图像，通过比较某一点在本尺度空间和上下尺度空间的响应值的大小，来判断某一点是否为候选特征点，若响应值比26个邻域值都大或都小，则将该点作为最终的候选特征点，并计算其位置和尺度参数；

2)SURF特征点主方向分配：首先进行Haar小波运算，具体参数是：6s为半径，特征点为中心，边长为4s，获取点在x,y方向上的 Haar小波响应值，其中s为空间尺度；然后进行高斯加权运算，具体参数是：张角为π/3的扇形滑动窗口，步长为0.2弧度滑动窗口，对窗口内图像的Haar小波响应值dx,dy进行累加，获取矢量(m_ω,θ_ω)：

求取Haar小波响应值累加值在多个方向中值最大的一个方向，则将该方向作为特征点的主方向；

3)SURF特征点特征矢量生成：构建以20s为边长，特征点为中心，方向与特征点主方向保持一致，大小为4×4的规则子窗口。采用边长为2σ的Haar小波对图像进行处理，以获得x,y方向上的响应值 dx,dy，利用高斯加权计算每一个子窗口的响应值，以获得每一个子窗口的特征矢量：

υ_子窗口＝[∑dx∑dy∑|dx|∑|dy|]

一组描述子特征矢量共包含4×4×4＝64维特征矢量，可获取一个特征点的完整信息：空间尺度、坐标、64维矢量特征；

4)RANSAC算法精确匹配：根据左视图像像特征点匹配右视图像，右视图像特征点匹配左视图像，再进行筛选，如果两次匹配都能够成功，就把进行匹配点对存入新的数组，进行RANSAC模型估计，判断出正确的匹配点对，进行n次迭代计算，以得到最终的匹配点和转换矩阵。

其中，还包括以下步骤：

针对倾斜影像数据、GNSS数据及以角元素Roll(Φ)、Pitch(Θ)、 Heading(Ψ)作为IMU数据初始值的组合导航数据进行空三，分别迭代计算倾斜影像数据的位置信息及姿态信息，剔除残差、粗差，当达到最佳迭代次数时，若空三结果存在航线集与X-Y平面不平行，则不再继续迭代计算，且判定该结果失败；若所有航线集与X-Y平面平行，则不再继续迭代计算，且判定该结果可行；当没有达到最佳迭代次数时，若空三结果存在航线集与X-Y平面不平行，则继续进行迭代计算倾斜影像数据的位置信息及姿态信息，剔除残差、粗差，直到所有航线集与X-Y平面平行为止，且判定该结果可行；将可行的空三结果作为新的初始值继续空三，迭代计算倾斜影像数据的位置信息及姿态信息，剔除残差、粗差，观察每次计算后的连接点数量以及每次计算后航线集与X-Y平面是否平行，当连接点数量达到最多且航线集与 X-Y平面平行时，则选择该次空三结果作为新的初始值继续空三，迭代优化倾斜影像数据的位置信息及姿态信息，剔除残差、粗差，观察每次优化后的连接点数量以及每次计算后航线集与X-Y平面是否平行，当连接点数量达到最多且航线集与X-Y平面平行时，若GNSS数据没有采用RTK测量，则结束优化，得出空三结果及空三报告；若 GNSS数据采用了RTK测量，则选择该次空三结果作为新的初始值继续空三，迭代优化倾斜影像数据的位置信息及姿态信息，剔除残差、粗差，观察每次优化后的连接点数量以及每次计算后航线集与X-Y平面是否平行，当连接点数量达到最多且航线集与X-Y平面平行时，选择该次空三结果作为最终结果，比较两次空三最终结果，选择连接点多且精度高的用于实景三维建模。

其中，最佳迭代次数是6－10次。

其中，多次迭代计算、优化是为了剔除残差、粗差，使像片的位置、姿态正确，同时使连接点数足够多，确保建模所用到的三角形网格数足够多，从而提升三维模型的质量。

表1为本发明实施例的空三精度报告表：

由上表可知，空三精度基本不变，已经趋于稳定，可以看到“连接点数”已达到最大，这时需要观察航线集与X-Y平面是否平行，若平行，此次空三正常结束。

绝对高度：是指在地面或海面上空的飞机距地表或海面的垂直距离，亦称“海拔高度”。即距标准大气海平面的垂直距离，或者说以平均海平面做标准的高度。航行地图上标出的地形、地物的高度，都是按绝对高度计算的。

相对高度：飞行器高出某指定地点的垂直距离。

SURF：全称是Speed-up robust features，即加速鲁棒性特征算法，是一种高鲁棒性的局部特征点检测器。该算法可以用于计算机视觉领域的物体识别或者三维重建。

RANSAC：是Random Sample Consensus的缩写，它是根据一组包含异常数据的样本数据集，计算出数据的数学模型参数，得到有效样本数据的算法。RANSAC算法经常用于计算机视觉中。例如，在立体视觉领域中同时解决一对相机的匹配点问题及基本矩阵的计算。

X—Y平面：为平行于平均海平面的平面。

航线集：为若干张倾斜影像组成的集合。

倾斜摄影测量：倾斜摄影测量技术是国际测绘领域近些年发展起来的一项高新技术，它颠覆了以往正射影像只能从垂直角度拍摄的局限，通过在同一飞行平台上搭载多台传感器，同时从一个垂直、四个倾斜等五个不同的角度采集影像，将用户引入了符合人眼视觉的真实直观世界。航空倾斜影像不仅能够真实地反应地物情况，而且还通过采用先进的定位技术，嵌入精确的地理信息、更丰富的影像信息、更高级的用户体验，极大地扩展了遥感影像的应用领域，并使遥感影像的行业应用更加深入。由于倾斜影像为用户提供了更丰富的地理信息，更友好的用户体验，该技术在欧美等发达国家已经广泛应用于应急指挥、国土安全、城市管理、房产税收等行业。

实景三维建模：是指根据一系列二维相片，或者一组倾斜影像，自动生成高分辨的、带有逼真纹理贴图的三维模型。如果倾斜像片带有坐标信息，那么模型的地理位置信息也是准确的。这种模型效果逼真，要素全面，而且具有测量精度，不仅带给人身临其境之感还可用于测量学应用，是现实世界的真实还原。

IMU：即惯性测量单元，是测量物体三轴姿态角(或角速率)以及加速度的装置；

GNSS：即Global Navigation Satellite System的缩写，即全球导航卫星系统。早在20世纪90年代中期开始，欧盟为了打破美国在卫星定位、导航、授时市场中的垄断地位，获取巨大的市场利益，增加欧洲人的就业机会，一直在致力于民用全球导航卫星系统计划，称之为 Global Navigation Satellite System。该计划分两步实施：第一步是建立一个综合利用美国的GPS系统和俄罗斯的GLONASS系统的第一代全球导航卫星系统(当时称为GNSS-1，即后来建成的EGNOS)；第二步是建立一个完全独立于美国的GPS系统和俄罗斯的GLONASS系统之外的第二代全球导航卫星系统，即正在建设中的Galileo卫星导航定位系统。由此可见，GNSS从一问世起，就不是一个单一星座系统，而是一个包括GPS、GLONASS等在内的综合星座系统；差分GNSS指通过利用位置已知的参考GNSS接收机的额外数据来减小GPS系统或 GLONASS系统定位误差的一项技术。

组合导航数据：指卫星导航数据(GNSS数据)和惯性导航数据 (IMU数据)组合在一起的综合导航数据，包含了目标物的位置信息和姿态信息。

连接点(Tiepoints)：立体像对重叠范围内,同一物点在不同像片上的构像点称为同名像点，大量自动或手动生成的同名像点统称连接点。

影像匹配：通过一定的匹配算法在两幅或多幅影像之间识别同名点的过程。

空三：即空中三角测量，空中三角测量是立体摄影测量中，根据少量的野外控制点，在室内进行控制点加密，求得加密点的高程和平面位置的测量方法。其主要目的是为缺少野外控制点的地区测图提供绝对定向的控制点。空中三角测量一般分为两种:模拟空中三角测量即光学机械法空中三角测量；解析空中三角测量即俗称的电算加密。模拟空中三角测量是在全能型立体测量仪器(如多倍仪)上进行的空中三角测量。它是在仪器上恢复与摄影时相似或相应的航线立体模型，根据测图需要选定加密点，并测定其高程和平面位置。航空摄影测量中利用像片内在的几何特性，在室内加密控制点的方法。即利用连续摄取的具有一定重叠的航摄像片,依据少量野外控制点，以摄影测量方法建立同实地相应的航线模型或区域网模型(光学的或数字的)，从而获取加密点的平面坐标和高程。主要用于测地形图。

RTK：即实时动态差分法。这是一种新的常用的GPS测量方法，以前的静态、快速静态、动态测量都需要事后进行解算才能获得厘米级的精度，而RTK是能够在野外实时得到厘米级定位精度的测量方法，它采用了载波相位动态实时差分方法，是GPS应用的重大里程碑，它的出现为工程放样、地形测图，各种控制测量带来了新曙光，极大地提高了外业作业效率。RTK定位技术就是基于载波相位观测值的实时动态定位技术，它RTK能够实时地提供测站点在指定坐标系中的三维定位结果，并达到厘米级精度。在RTK作业模式下，基准站通过数据链将其观测值和测站坐标信息一起传送给流动站。流动站不仅通过数据链接收来自基准站的数据，还要采集GPS观测数据，并在系统内组成差分观测值进行实时处理，同时给出厘米级定位结果，历时不足一秒钟。流动站可处于静止状态，也可处于运动状态；可在固定点上先进行初始化后再进入动态作业，也可在动态条件下直接开机，并在动态环境下完成整周模糊度的搜索求解。在整周未知数解固定后，即可进行每个历元的实时处理，只要能保持四颗以上卫星相位观测值的跟踪和必要的几何图形，则流动站可随时给出厘米级定位结果。

Roll、Pitch、Heading：即侧滚角、俯仰角和旋偏角，是惯性测量单元(IMU)通常采用的、用以描述传感器姿态的角元素系统；

显然，本发明的上述实施例仅仅是为清楚地说明本发明所作的举例，而并非是对本发明的实施方式的限定。对于所属领域的普通技术人员来说，在上述说明的基础上还可以做出其它不同形式的变化或变动。这里无法对所有的实施方式予以穷举。凡是属于本发明的技术方案所引伸出的显而易见的变化或变动仍处于本发明的保护范围之列。