CN109443359B - 一种地面全景影像的地理定位方法 - Google Patents

Abstract

本发明一种地面全景影像的地理定位方法，先采用一种对空间组织方位关系整体性描述的全方位特征量，通过在卫星影像地图网格点处描述向量的检索完成粗定位，然后由粗定位确定若干匹配线作为高度信息失真时的控制线，最后在满足相对高度顺序不变的约束条件下，通过粗定位引入迭代初值，实现外方位元素的精确解算。不同于多帧地面影像的光束平差解算，该方法仅需要中高分辨率的卫星影像数据和场景表面相对高程信息，无需任何其它定位装置(如GPS)，也无需添加人工控制点，即可实现单帧地面全景影像的精确外参解算，非常适合于天、地这样视场角、尺度差异巨大的影像配准。

Description

一种地面全景影像的地理定位方法

所属技术领域

本发明属于遥感信息系统技术领域，涉及一种地面全景影像的地理定位方法。

背景技术

基于视觉的地面影像定位(即外方位元素解算)在城市三维增强现实、移动机器人等领域有重要价值，在无GPS信号或无线电等手机定位信号情况下，可通过地面影像与卫星/机载遥感影像数据的融合配准，获得当前影像在地理空间中的标定，完成精确的定位定向。这种仅通过视觉影像信息进行定位定向的方法，不仅避免使用复杂的传统地面移动测量装置(配备严格标定的LiDAR和全景相机、精准的GPS、IMU测位测姿单元)，而且这种空天地数据的融合，还将有助于发挥不同视角观测优势，为场景三维重建提供全面信息，而不再受传统测绘领域DEM/DSM产品无侧面纹理信息的限制。由于全景影像可提供全方位信息，较普通影像更容易提取稳定的几何相关测量信息，因此能提供单帧影像更可靠的外方位解算。

传统基于视觉的定位方法是利用标志点，搜索其在数据库或地形图中的位置，实现外方位解算。标志点通常需要手工挑选，或是采用特征点匹配方法得到。但由于特征点在超大视角更迭时已经变得面目全非，即使是人眼在面对地面观察和卫星俯视观察时也很难给出匹配点，所以常规特征点配准方法通常无法处理如此大视角变化。利用高程模型进行配准的定位方法，不受植被四季变化、光照变化的影响，定位稳定，文献“DynamicProgramming and Skyline Extraction in Catadioptric Infrared Images”用动态规划分割提取天际线，和高程模型匹配，估计相机的姿态和位置。文献“Large Scale VisualGeo-localization of Images in Mountainous Terrain”提取山区的天际线，在DEM数据库中查找与之对应的变换，可在4万平方公里定位在1公里以内。这类方法受高程模型精度限制，多适用于大场景粗定位，而且适合山区或开阔地区，对城市区域并不适合。文献“Registering Ground and Satellite Imagery for Visual Localization”提出了一种基于building footprint orientation(FPO)描述向量的全景影像定位方法，适用于城市环境，但它仅能给出粗略定位，精度不高。

发明内容

本发明解决的技术问题是：克服现有技术的不足，提供一种基于卫星影像配准的单帧全景影像地理定位方法，该方法无需人工添加任何控制点，仅需要中高分辨率的卫星影像数据和不精确的场景表面相对高程信息，即可精确定位单帧地面全景影像的空间方位信息。

本发明的技术解决方案是：一种地面全景影像的地理定位方法，包括如下步骤：

1)提取卫星影像地图网格点处的全方位特征量；

2)提取地面全景影像全方位特征量；

3)对地面全景影像进行粗定位，得到地面全景影像的大致位置和大致姿态；

4)对地面全景影像进行基于线匹配的地面全景影像精确地理定位。

所述步骤1)中进行卫星影像地图网格点全方位特征量提取时，选择不同分辨率的光学卫星遥感图像，或选用正射校正后的航空影像。

所述步骤1)的具体过程为：

11)针对卫星影像地图，提取建筑物顶面轮廓，然后将建筑物顶面轮廓经直线段拟合，提取出建筑顶面的外侧轮廓直线段；

12)在卫星影像地图中均匀布置网格点，计算每个网格点处的全方位特征量。

所述步骤12)中，所述全方位特征量为V×D矩阵，其中V是视线方向的量化数目，D是深度位数，该特征量记录了不同视线方向上的视线与步骤11)获得的建筑顶面的外侧轮廓直线段的夹角，并分成D个不同深度。

所述步骤2)的具体过程为：

21)对地面全景影像提取直线段，进一步提取消失点，得到建筑物垂直立面的直线段簇和建筑物平行地面的直线段簇；

22)利用建筑物垂直立面直线段组成的消失点确定水平地面法向量，进而确定从地面全景影像成像坐标系到地面水平面坐标系的校正角度α,β；

23)通过计算消失点，获得建筑物平行水平地面的直线段簇的空间方向，将建筑物平行水平地面的直线段簇映射到水平地面，获得地面全景影像的全方位特征量。

所述全方位特征量为V×D矩阵，其中V是视线方向的量化数目，D是深度位数，该特征量记录了不同视线方向上的视线与步骤21)获得的建筑物平行水平地面的直线段簇的夹角，并分成D个不同深度。

所述步骤22)中计算得到校正角度α,β的具体公式为：

G_{N_c}＝[cosβ·cosα,sinβ,cosβ·sinα]

其中，G_{N_c}为已知的水平地面法向量。

所述步骤3)粗定位的具体过程为：将全方位特征量循环上移，每次移动一行，将每次移动后的全方位特征量与卫星影像地图中每个网格点处的全方位特征量进行做差比较，选取差值最小所对应的网格点和循环上移量，并根据网格点和循环上移量推算出地面全景影像的大致位置和大致姿态γ₀。

所述步骤4)中进行基于线匹配的地面全景影像精确地理定位时，目标函数的迭代求解采用梯度下降法，或高斯牛顿法，或Levenberg-Marquardt法。

所述步骤4)匹配线的确定过程如下：

a1)对卫星影像地图上的每条建筑外轮廓直线段L_m向地面全景影像大致位置的全景成像面引出射线，找到由射线围成的地面全景影像建筑物平行地面的直线段簇集合{l_k}；

a2)经坐标转换将每条直线l_k的空间方向变换到卫星影像地图坐标系下，剔除与L_m方向夹角大于5°的线段；

a3)由步骤a2)得到的剔除后的剩余直线段簇的相对高度，进一步剔除非建筑物顶面的线段，得到卫星影像地图和地面全景影像中的匹配线。

所述步骤4)精确定位的具体过程为：引入权重矩阵W，构造log-barrier目标函数：

其中，Δγ是待求角度调整量，T_x,T_y是待求解的全景相机在卫星影像地图上的位置，Z_{1_m},…,Z_{s_m}是未知的S个线对的空间高度，τ_i是比例项，

是基于地面全景影像模式的匹配线Plucker投影关系，μ是值大于0的约束项参数，以满足空间直线相对高度顺序不变的约束；权重矩阵W_3S×3S＝diag(w₁I_3×3,…,w_sI_3×3)，若第i对匹配线是正确的，则W_i为1，若被判断为误匹配，则w_i为0；

粗定位提供迭代求解的一个有效初始值，权重矩阵中w_i初始值全设置为1，在每次迭代过程中，进行粗差筛查，若单项匹配线的

值是所有匹配线

平均值的3倍且大于设定阈值，则该项对应权值w_i设为0；利用log-barrier目标函数求解方法得到真值

从而得到精确的外方位参数

其中T_z是已知的全景相机真实高度。

本发明的优点在于：

(1)、本发明不同于基于特征点的常规配准方法，先利用全景影像的全方位特征量获得其在卫星影像地图全方位特征量库中最接近的网格点，完成粗定位和线匹配，然后利用匹配线完成相对高度约束下的方向角和位置的精确求解，解决了天、地大视角差异下的影像配准问题；

(2)、本发明不同于常规的地面多点成像光束平差法，在无任何人工控制点，不依赖精确的高程模型时，仅给定较高分辨率的卫星影像和相对高度顺序，即可实现单帧全景影像的精确外方位元素解算，同时估算了匹配线的空间高度。

附图说明

图1为本发明一种地面全景影像的地理定位方法工作流程图；

图2是卫星影像地图上某网格点处全方位特征量的灰度编码显示；

图3是地面全景影像坐标系到卫星影像地图坐标系的转换关系示意图；

图4是已知空间方向的直线段在水平面上的投影示意图。

具体实施方式

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

本发明提出了一种地面全景影像的地理定位方法，该方法分粗定位和精定位两个阶段，前者通过在卫星影像地图网格点处全方位特征量的检索确定大致位姿；后者利用匹配线作为高度信息失真时的控制线，构建约束最小化问题，求解精确外方位元素。相较于现有配准方法，该方法适合于天、地这样视场角、尺度差异巨大的情况，而且极大降低了对高程精度的依赖，可实现单帧全景影像的无控制点自动精确定位。如图1所示为本发明一种地面全景影像的地理定位方法工作流程图，本发明共分成四个步骤完成，具体如下：

步骤(一)：提取卫星影像地图网格点处的全方位特征量

针对卫星影像地图，提取建筑物顶面轮廓，然后经直线段拟合，提取出建筑顶面轮廓直线段。卫星影像地图由N个建筑立面f₁，f₂，...，f_N组成，其中每个

是由立面投影到水平地面(在卫星图像中的像素坐标)

以及立面相对高度h_{r_m}组成。

在卫星影像地图中均匀布置网格点C_m，计算每个网格点处的全方位特征量。全方位特征量F是一个V×D矩阵，其中V是视线方向的量化数目，D是深度位数。对于从网格点C_m引出的视线方向θ，以此视线为中心考察其两侧的5条视线，对由这些视线与投影立面的夹角进行聚类，然后将具有最多成员数量的D个聚类的均值分配给此视线。如果少于D个聚类，则将此角度的深度赋值设置为-1。这里，我们取V＝72，D＝7。图2上标出了每个网格点，并展示了其中一个网格点处的全方位特征量。在从0到360度的每个观察方向上，最多可识别D个立面方向(0到180度之间)，每个环代表不同深度值处的立面角度，由近及远立面高度依次增高，若条目是白色，表示该视线方向在此立面高度没有标识。

步骤(二)：提取地面全景影像全方位特征量

对地面全景影像提取直线段，进一步提取消失点。由于全景成像面涵盖了消失点，所以可以利用建筑物垂直立面直线段组成的消失点来确定水平地面法向量G_{N_c}，从而确定从地面全景影像成像坐标系到地面水平面坐标系的转换矩阵R_c→hor，满足：

其中，

是地面全景影像成像坐标系下的空间点坐标，

是在水平地面坐标系下的坐标，坐标系符合右手定则，坐标系的z轴垂直水平地面向下，x轴指向正南方。

为了准确确定水平地面的法向量，由于消失点成对出现，对所有建筑物垂直立面直线段簇的交点进行k-means聚类，设置成两大类，选择数量最多一簇的中心即为水平地面的法向量G_{N_c}。由法向量可以进一步确定转换矩阵R_c→hor，见式(3)。

G_{N_c}＝[cosβ·cosα,sinβ,cosβ·sinα] (2)

其中，β是G_{N_c}向量与其在x-z面上投影的夹角，α是x-z面上投影与x轴的夹角。R_axis作用是将旋转后的坐标系与水平地面坐标系的x″-y″-z″方向对齐。坐标系转换关系如图3所示。

通过计算消失点，可以获得直线段的空间方向，然后将这些已知空间方向的直线段映射到水平地面，具体如下：图4中，L_c是在地面全景影像成像坐标系下的空间直线，l_c是其在全景成像面上的投影线段(端点p_s和端点p_e)，由于空间直线的三维方向可以在计算消失点时同时得到，因此可确定空间点p_E，过p_E和p_s的直线在三角形Δp_sCp_e所在的平面上，并且平行于L_c在此平面上的投影。采用此方法，可以获得每条空间直线在地面全景影像成像坐标系下的对应直线段，从而可以计算出它们在水平地面，以C为坐标原点，坐标系x″-y″平面上的投影。按照步骤(一)类似的方法，计算出地面全景影像的全方位特征量。由于不知道相对高度，对一个方向我们选取D个与卫星影像地图上建筑顶面轮廓直线段的夹角，不足D个的补-1。

步骤(三)：地面全景影像粗定位

由于全方位特征量不具有旋转不变性，所以该向量要循环上移逐一和卫星影像地图网格点库中的向量进行比对，选取误差最小的所在的地图网格点和循环上移量推算其粗略的位置和姿态。将全方位特征量每向上循环移动一行，相当于绕z轴旋转dγ＝360°/72＝5°，会生成72种不同起始方向时的描述向量，将它们与卫星影像地图中每个网格点处的描述向量比较，距离量按式(6)计算，选取令距离差最小的网格点位置和旋转角度。

其中，

表示地面全景影像全方位特征量上移

行的结果，g_(u,v)表示网格点(X_{u_m},Y_{v_m})的全方位特征量，(m,n)表示描述矩阵中第m行，第n列的元素，q是0或1的异或操作值，当

和g_(u,v)的第m行全为-1时，或

和g_(u,v)的第m行全不为-1时，取1。e是一个常数，由当前所有有效行的平均差值决定，K是一个比例系数，对于无效行，指定的误差应比平均误差大，所以一般选择大于1的值，经验选择1.5～3之间。上移

行等同于绕z轴旋转

此时

通过上述方法，找到距离拍摄点最近的网格点作为粗定位的初始值

和最合适的夹角作为初始方位元素γ₀。

步骤(四)：基于线匹配的地面全景影像精确地理定位。

粗定位可以匹配得到卫星影像地图上建筑顶面轮廓直线段和它们在地面全景影像中的投影线段。从而进一步由线匹配解算精确外方位参数。匹配线的确定过程如下：1)对卫星影像地图上的每条建筑外轮廓直线段L_m向地面全景影像大致位置的全景成像面引出射线，找到由射线围成的地面全景影像建筑物平行地面的直线段簇集合{l_k}；2)由每条直线l_k的空间方向经下式(9)可以得到卫星影像地图坐标系下的方向，剔除与L_m方向夹角大于5度的线段；3)由剔除后的剩余直线段簇的相对高度再进一步剔除非建筑物顶面的线段，得到卫星影像地图和地面全景影像中的匹配线。

Rot＝R_z(γ₀)·R_axis·R_y(α)·R_x(β) (9)

卫星影像地图上的空间直线由空间点齐次坐标P_{1_m}＝(X_{1_m},Y_{1_m},Z_{1_m},1)^T和P_{2_m}＝(X_{2_m},Y_{2_m},Z_{2_m},1)^T确定。由于是建筑物顶面线段，一般平行于地面，两个点的高度是一致的，但由于高程信息误差较大，高度作为待求解量。空间直线的Plucker坐标表示为

分别代表了卫星影像地图坐标系原点和空间直线确定平面的法向量和空间线段方向向量。地面全景影像的直线段和卫星影像地图的空间直线满足下式

其中，

是3×4的投影矩阵，τ是非零比例项，

是由地面全景影像中心和全景成像面直线段投影确定的平面法向量，由直线段投影两个端点的地面全景影像坐标系坐标P_{1_c}＝(X_{1_c},Y_{1_c},Z_{1_c})^T和P_{2_c}＝(X_{2_c},Y_{2_c},Z_{2_c})^T可以得到

而此时的投影矩阵P描述为

P＝[Rot^-1|-Rot^-1·T_m] (13)

其中，T_m是地面全景影像在卫星影像地图坐标系下的坐标(T_x,T_y,T_z)^T，T_z由相机成像高度决定，可以精确得知。对于S对线对，通过

可以构建3S个方程，有3+2S个未知数，包括1个姿态角Δγ，2个待定位的位置元素T_x,T_y，以及所有卫星影像地图上空间直线段在z方向的高度信息Z_m，和比例项τ。当n≥3时，即可解算。令

由于已经过了粗定位，R_z(Δγ)可以用小角近似为一个线性关系，则每对匹配线有下式：

(Y_{1_m}-Y_{2_m})(Z_m-T_z)＝τ(l_x-Δγl_y)

(X_{2_m}-X_{1_m})(Z_m-T_z)＝τ(Δγl_x+l_y)

只要S≥3，就可以解算，忽略Δγ²的小角项，通过化简合并，可以得到如下方程组：

其中，b_i，c_i，d_i，f_i，g_i，h_i，k_i，u_i，v_i是根据X_{1_mi}，X_{2_mi}，Y_{1_mi}，Y_{2_mi}，T_z解算的中间量，此外，所有匹配线的相对高度有个大小排序，所以约束条件为：

Z_{1_m}≤Z_{2_m}≤…≤Z_{S_m}

这里构造了log-barrier目标函数，见式(16)，

其中，

是方程组(15)的形式，μ是值大于0的约束项参数，权重矩阵W_3S×3S＝diag(w₁I_3×3,…,w_sI_3×3)，若第i对匹配线是正确的，则w_i为1，若被判断为误匹配，则w_i为0。

粗定位可以提供迭代求解的一个有效初值，权重矩阵中w_i初值全设置为1，在每次迭代过程中，进行粗差筛查，若某匹配线的

误差是所有匹配线

平均误差的3倍且大于某一阈值，则该匹配线对应的w_i设为0。利用常规的log-barrier目标函数求解方法可以得到真值

从而得到精确的外方位参数

其中T_z是已知的全景相机真实高度。

在上述一种地面全景影像的地理定位方法中，步骤(一)中进行卫星影像地图采样点全方位特征量提取时，可以选择不同分辨率的光学卫星遥感图像，也可以选用正射校正后的航空影像，只要地面轮廓直线段清晰便于提取即可。步骤(四)中进行基于线匹配的地面全景影像精确地理定位时，目标函数的迭代求解可以采用梯度下降法、高斯牛顿法、Levenberg-Marquardt(LM)法，并且不限于这些方法。

以上所述，仅为本发明的较佳实施例，并非对本发明做任何限制。凡是根据发明技术实质对以上实施例所做的任何简单修改、变更及等效变化，仍属于本发明技术方案的保护范围。

Claims (8)

1.一种地面全景影像的地理定位方法，其特征在于包括如下步骤：

1)提取卫星影像地图网格点处的全方位特征量；

2)提取地面全景影像全方位特征量；

3)对地面全景影像进行粗定位，得到地面全景影像的大致位置和大致姿态；

4)对地面全景影像进行基于匹配线的地面全景影像精确地理定位；

所述步骤4)中进行基于匹配线的地面全景影像精确地理定位时，目标函数的迭代求解采用梯度下降法，或高斯牛顿法，或Levenberg-Marquardt法；

所述步骤4)匹配线的确定过程如下：

a1)对卫星影像地图上的每条建筑外轮廓直线段L_m向地面全景影像大致位置的全景成像面引出射线，找到由射线围成的地面全景影像建筑物平行地面的直线段簇集合{l_k}；

a2)经坐标转换将每条直线l_k的空间方向变换到卫星影像地图坐标系下，剔除与L_m方向夹角大于5°的线段；

a3)由步骤a2)得到的剔除后的剩余直线段簇的相对高度，进一步剔除非建筑物顶面的线段，得到卫星影像地图和地面全景影像中的匹配线；

所述步骤4)精确定位的具体过程为：引入权重矩阵W，构造log-barrier目标函数：

其中，Δγ是待求角度调整量，T_x,T_y是待求解的全景相机在卫星影像地图上的位置，Z_{1_m},…,Z_{S_m}是未知的S个线对的空间高度，τ_i是比例项，

是基于地面全景影像模式的匹配线Plucker投影关系，μ是值大于0的约束项参数，以满足空间直线相对高度顺序不变的约束；权重矩阵W_3S×3S＝diag(w₁I_3×3,…,w_sI_3×3)，若第i对匹配线是正确的，则w_i为1，若被判断为误匹配，则w_i为0；

粗定位提供迭代求解的一个有效初始值，权重矩阵中w_i初始值全设置为1，在每次迭代过程中，进行粗差筛查，若单项匹配线的

值是所有匹配线

平均值的3倍且大于设定阈值，则该项对应权值w_i设为0；利用log-barrier目标函数求解方法得到真值

从而得到精确的外方位参数

其中T_z是已知的全景相机真实高度。

2.根据权利要求1所述的一种地面全景影像的地理定位方法，其特征在于：所述步骤1)中进行卫星影像地图网格点全方位特征量提取时，选择不同分辨率的光学卫星遥感图像，或选用正射校正后的航空影像。

3.根据权利要求1或2所述的一种地面全景影像的地理定位方法，其特征在于：所述步骤1)的具体过程为：

11)针对卫星影像地图，提取建筑物顶面轮廓，然后将建筑物顶面轮廓经直线段拟合，提取出建筑顶面的外侧轮廓直线段；

12)在卫星影像地图中均匀布置网格点，计算每个网格点处的全方位特征量。

4.根据权利要求3所述的一种地面全景影像的地理定位方法，其特征在于：所述步骤12)中，所述全方位特征量为V×D矩阵，其中V是视线方向的量化数目，D是深度位数，该特征量记录了不同视线方向上的视线与步骤11)获得的建筑顶面的外侧轮廓直线段的夹角，并分成D个不同深度。

5.根据权利要求1所述的一种地面全景影像的地理定位方法，其特征在于：所述步骤2)的具体过程为：

21)对地面全景影像提取直线段，进一步提取消失点，得到建筑物垂直立面的直线段簇和建筑物平行地面的直线段簇；

22)利用建筑物垂直立面直线段组成的消失点确定水平地面法向量，进而确定从地面全景影像成像坐标系到地面水平面坐标系的校正角度α,β；

23)通过计算消失点，获得建筑物平行水平地面的直线段簇的空间方向，将建筑物平行水平地面的直线段簇映射到水平地面，获得地面全景影像的全方位特征量。

6.根据权利要求5所述的一种地面全景影像的地理定位方法，其特征在于：所述全方位特征量为V×D矩阵，其中V是视线方向的量化数目，D是深度位数，该特征量记录了不同视线方向上的视线与步骤21)获得的建筑物平行水平地面的直线段簇的夹角，并分成D个不同深度。

7.根据权利要求5所述的一种地面全景影像的地理定位方法，其特征在于：所述步骤22)中计算得到校正角度α,β的具体公式为：

G_{N_c}＝[cosβ·cosα,sinβ,cosβ·sinα]

其中，G_{N_c}为已知的水平地面法向量。

8.根据权利要求1所述的一种地面全景影像的地理定位方法，其特征在于：所述步骤3)粗定位的具体过程为：将全方位特征量循环上移，每次移动一行，将每次移动后的全方位特征量与卫星影像地图中每个网格点处的全方位特征量进行做差比较，选取差值最小所对应的网格点和循环上移量，并根据网格点和循环上移量推算出地面全景影像的大致位置和大致姿态γ₀。

Images