CN103017653B - 一种球面全景影像与三维激光扫描点云的配准及测量方法 - Google Patents

Abstract

本发明提供了一种球面全景影像与三维激光扫描点云的配准及测量方法，首先用单位球体的三维结构模拟球面全景影像，将球面全景影像的像素坐标转换为三维直角坐标，然后计算出球面全景影像在点云数据中的姿态参数以及球面全景影像的光心在点云数据中的位置坐标将球面全景影像在点云数据中进行配准；再通过计算得到球面全景影像的深度图；最后当需要在球面全景影像上测量时，将鼠标移动到需要测量的目标点位置，即可计算得出目标点的三维坐标。本发明提供的配准及测量方法不仅能够将球面全景影像与三维激光扫描点云进行高度配准，同时还能够对目标点实现快速测量。

Description

一种球面全景影像与三维激光扫描点云的配准及测量方法

技术领域

本发明涉及一种将球面全景影像与三位激光扫描点云进行配准以及配准后根据图像进行快速测量的方法，属于光电测绘领域。

背景技术

三维激光点云是物体表面的高精度三维坐标，通过点云直接可以获取地物坐标及特征信息，但没有纹理和颜色信息，表现不直观。球面全景影像数据具有形象的纹理和颜色信息，可以作为三维激光点云数据的纹理补充。三维激光点云数据和球面全景影像数据各有优势和特长，倘若将三维激光点云和全景影像数据结合使用，那么既可以通过三维激光点云可以找到对应影像像素，进行点云着色，反之以球面全景影像作为展现载体，通过影像素坐标可以找到对应三维激光点云，进行地物量测。但是要实现上述目的必须解决两者之间的配准问题，该问题目前是研究的热点。同时，如何在在配准后的球面全景视图中根据任意像素快速获取的点云坐标也是亟待解决的问题，由于海量的点云数据不适合于互联网的应用，想要解决球面全景影像快速测量问题必须解决点云数据的存储难题。

发明内容

本发明提供了一种球面全景影像与三维激光扫描点云的配准及测量方法，解决了上述背景技术中的难题，该方法不仅能够将球面全景影像与三维激光扫描点云进行高度配准，同时还能够对目标点实现快速测量。

实现本发明上述目的所采用的技术方案为：

一种球面全景影像与三维激光扫描点云的配准及测量方法，包括以下步骤：

（1）、用单位球体的三维结构模拟球面全景影像，将球面全景影像的像素坐标（u，v）首先转换为球面极坐标 ，然后将球面全景影像的球面极坐标转换为三维直角坐标(x₀，y₀，z₀)；

（2）、在球面全景影像中选取三个以上的分别位于不同平面上的点作为控制点，同时与之对应的在点云数据中选取上述控制点，根据控制点在球面全景影像中的坐标以及在点云数据中的坐标，即可计算得出球面全景影像在点云数据中的姿态参数(φ，ω，κ)，以及球面全景影像的光心在点云数据中的位置坐标(Xc，Yc，Zc)，从而将球面全景影像在点云数据中进行配准；

（3）、以球面全景影像的光心位置坐标(Xc，Yc，Zc)为原点，搜索空间范围内的 点云，将搜索到的每一个点至光心的距离值d均转化为RGB值，将转化后的RGB值以图片格式存储，得到球面全景影像的深度图；

（4）、当需要在球面全景影像上测量时，将鼠标移动到需要测量的目标点位置，此时读取目标点的像素坐标（u，v）以及球面极坐标 ，根据姿态参数(φ，ω，κ)以及像素坐标（u，v）计算目标点在球面全景影像的深度图中的像素坐标（u′，v′），读取此像素中的RGB值，并将其转化为距离值d，最后根据距离值d和角度坐标 ，即可计算得出目标点的三维坐标(x,，y，z)。

步骤（1）中像素坐标（u，v）转换为球面极坐标 的公式为：

，

式中，dx为球面全景影像的纵向像素大小，dy为球面全景影像的横向像素大小；球面极坐标 转换为三维直角坐标(x₀，y₀，z₀)的计算公式为：

。

步骤（2）中姿态参数(φ，ω，κ)以及球面全景影像的光心位置坐标(Xc，Yc，Zc)的计算方法为：根据点云数据中的点、球面全景影像中的点以及球面全景影像球心三点共线，得到以下关系式：

其中R为旋转矩阵， 

转换后得：

其中：

，

，

其中a₁、a₂、a₃、b₁、b₂、b₃、c₁、c₂、c₃为旋转矩阵的元素；

将转换后的公式展开，得：

设： ,

则： ,

,

又设： ，

其中α、β为观测值，其相应的改正数为 v_α、 v_β；将上式用泰勒公式展开成线性形式，即可得到线性化的误差方程如下：

，

其中，（α）、（β）为使用待定参数的近似值带入公式中求出的α、β的计算值，ΔXc、ΔYc、ΔZc、Δφ、Δω、Δκ为未知数Xc、Yc、Zc、φ，ω，κ的改正数；

然后采用以下公式计算α、β对各个未知数的偏导数：

再将误差方程转化为矩阵形式：

其中，

根据最小二乘平差原理，列出法方程式：

则法方程的解为：

从而求出待测未知数的近似值的改正数矩阵H的6元素ΔXc、ΔYc、ΔZc、Δφ、Δω、Δκ；

最后对计算结果进行重复迭代计算，每次迭代时用未知数近似值与上次迭代计算的改正数之和作为新的近似值，重复计算过程，求出新的改正数，依此反复趋近，直至改正数小于设定的阈值，最后得到未知数的解，为：

其中Xc⁰为Xc的近似值，ΔXc¹为第一次计算的改正数，ΔXc²为第二次计算的改正数，其他数据均依此类推。

步骤（3）中的具体步骤为：首先预设一张球面全景影像的深度图，根据角度坐标 以及像素坐标将该预设深度图划分为网格状，并开辟保存每个网格中保存的点云距离的数据组；然后以球面全景影像的光心位置坐标(Xc，Yc，Zc)为原点，搜索空间范围内的点云，根据搜索到的点以及光心坐标，计算点的极坐标 ，以及到光心的距离值d；根据点的极坐标 ，计算确定点所对应深度图中的具体网格，并更新该网格中的点的距离值d；最后将所有距离值d均转化为RGB值，将转化后的RGB值一一存储于对应网格中，即可得到球面全景影像的深度图。

当有多个点同时对应到一个网格中时，选取多个点中的距离值最小者。

本发明提供的球面全景影像与三维激光扫描点云的配准及测量方法有以下优点：将三维激光点云和全景影像数据配准后相结合，既可以通过三维激光点云可以找到对应影像像素，进行点云着色，反之以球面全景影像作为展现载体，通过影像素坐标可以找到对应三维激光点云，进行地物量测。同时本发明中将大量的点云数据转化为深度图，以图片的形式进行存储，从而大大的减少了数据量，从而实现快速测量。

具体实施方式

下面结合具体实施例对本发明做详细具体的说明。

本实施例中采用以下方法，首先用单位球体的三维结构模拟球面全景影像，将球面全景影像的像素坐标（u，v）首先转换为球面极坐标 ，转换公式如下：

，

式中，dx为球面全景影像的纵向像素大小，dy为球面全景影像的横向像素大小；

然后将球面全景影像的球面极坐标转换为三维直角坐标(x₀，y₀，z₀)；

然后在球面全景影像中选取三个以上的分别位于不同平面上的点作为控制点，同时与之对应的在点云数据中选取上述控制点，根据控制点在球面全景影像中的坐标以及在点云数据中的坐标，即可计算得出球面全景影像在点云数据中的姿态参数( ，ω，κ)，以及球面全景影像的光心在点云数据中的位置坐标(Xc，Yc，Zc)，从而将球面全景影像在点云数据中进行配准；具体的计算步骤如下：首先根据点云数据中的点、球面全景影像中的点以及球面全景影像球心三点共线，得到以下关系式：

其中R为旋转矩阵， 

转换后得：

其中：

，

，

其中a₁、a₂、a₃、b₁、b₂、b₃、c₁、c₂、c₃为旋转矩阵的元素；

将转换后的公式展开，得：

设： ,

则： ,

,

又设： ，

其中α、β为观测值，其相应的改正数为v_α、 v_β； Xc、Yc、Zc、 ，ω，κ为待定的未知数，由于上式中观测值与未知数之间是非线性函数关系，为了便于计算，需把非线性函数表达式用泰勒公式展开成线性形式。将上式用泰勒公式展开后，即可得到线性化的误差方程如下：

，

其中，（α）、（β）为使用待定参数的近似值带入公式中求出的α、β的计算值，ΔXc、ΔYc、ΔZc、 、Δω、Δκ为未知数Xc、Yc、Zc、 ，ω，κ的改正数；

然后采用以下公式计算α、β对各个未知数的偏导数：

，求得的α、β对各个未知数的偏导数如下举例所示， ， 。

再将误差方程转化为矩阵形式：

其中，

根据最小二乘平差原理，列出法方程式：

则法方程的解为：

从而求出待测未知数的近似值的改正数矩阵H的6元素ΔXc、ΔYc、ΔZc、Δφ、Δω、 Δκ；

由于线性化的过程中各系数取自泰勒级数展开式的一次项，且未知数的初值一般比较粗略，因此计算需要迭代进行。每次迭代时用未知数近似值与上次迭代计算的改正数之和作为新的近似值，重复计算过程，求出新的改正数，依此反复趋近，直至改正数小于设定的阈值，最后得到未知数的解，为：

其中Xc⁰为Xc的近似值，ΔXc¹为第一次计算的改正数，ΔXc²为第二次计算的改正数，其他数据均依此类推。

进行配准后需要生成球面全景影像的深度图，具体方法如下：首先预设一张球面全景影像的深度图，根据角度坐标以及像素坐标将该预设深度图划分为网格状，并开辟保存每个网格中保存的点云距离的数据组；然后以球面全景影像的光心位置坐标(Xc，Yc，Zc)为原点，搜索空间范围内的点云，根据搜索到的点以及光心坐标，计算点的极坐标 ，以及到光心的距离值d；根据点的极坐标 ，计算确定点所对应深度图中的具体网格，并更新该网格中的点的距离值d；最后将所有距离值d均转化为RGB值，将转化后的RGB值一一存储于对应网格中，即可得到球面全景影像的深度图。

最后当需要在球面全景影像上测量时，将鼠标移动到需要测量的目标点位置，此时读取目标点的像素坐标（u，v）以及球面极坐标 ，根据姿态参数( ，ω，κ)以及像素坐标（u，v）计算目标点在球面全景影像的深度图中的像素坐标（u′，v′），读取此像素中的RGB值，并将其转化为距离值d，最后根据距离值d和角度坐标 ，即可计算得出目标点的三维坐标(x,，y，z)。

Claims (4)

1.一种球面全景影像与三维激光扫描点云的配准及测量方法，其特征在于包括以下步骤：

(1)、用单位球体的三维结构模拟球面全景影像，将球面全景影像的像素坐标(u，v)首先转换为球面极坐标然后将球面全景影像的球面极坐标转换为三维直角坐标(x₀，y₀，z₀)；

(2)、在球面全景影像中选取大于三个的分别位于不同平面上的点作为控制点，同时与之对应的在点云数据中选取上述控制点，根据控制点在球面全景影像中的坐标以及在点云数据中的坐标，即可计算得出球面全景影像在点云数据中的姿态参数以及球面全景影像的光心在点云数据中的位置坐标(Xc，Yc，Zc)，从而将球面全景影像在点云数据中进行配准；

(3)、以球面全景影像的光心位置坐标(Xc，Yc，Zc)为原点，搜索空间范围内的点云，将搜索到的每一个点至光心的距离值d均转化为RGB值，将转化后的RGB值以图片格式存储，得到球面全景影像的深度图；

步骤(3)中的具体步骤为：首先预设一张球面全景影像的深度图，根据角度坐标以及像素坐标将该预设深度图划分为网格状，并开辟保存每个网格中保存的点云距离的数据组；然后以球面全景影像的光心位置坐标(Xc，Yc，Zc)为原点，搜索空间范围内的点云，根据搜索到的点以及光心坐标，计算点的极坐标以及到光心的距离值d；根据点的极坐标计算确定点所对应深度图中的具体网格，并更新该网格中的点的距离值d；最后将所有距离值d均转化为RGB值，将转化后的RGB值一一存储于对应网格中，即可得到球面全景影像的深度图；

(4)、当需要在球面全景影像上测量时，将鼠标移动到需要测量的目标点位置，此时读取目标点的像素坐标(u，v)以及球面极坐标根据姿态参数以及像素坐标(u，v)计算目标点在球面全景影像的深度图中的像素坐标(u′，v′)，读取此像素中的RGB值，并将其转化为距离值d，最后根据距离值d和角度坐标即可计算得出目标点的三维坐标(x，y，z)。

2.根据权利要求1所述的球面全景影像与三维激光扫描点云的配准及测量方法，其特征在于：步骤(1)中像素坐标(u，v)转换为球面极坐标的公式为：

式中，dx为球面全景影像的纵向像素大小，dy为球面全景影像的横向像素大小；球面极坐标转换为三维直角坐标(x₀，y₀，z₀)的计算公式为：

3.根据权利要求1所述的球面全景影像与三维激光扫描点云的配准及测量方法，其特征在于：步骤(2)中姿态参数以及球面全景影像的光心位置坐标(Xc，Yc，Zc)的计算方法为：根据点云数据中的点、球面全景影像中的点以及球面全景影像球心三点共线，得到以下关系式：

$\left[\begin{array}{c}{X-X}_c\\ {Y-Y}_c\\ {Z-Z}_c\end{array}\right]=R\left[\begin{array}{c}{x}_0\\ y_0\\ z_0\end{array}\right]$

其中R为旋转矩阵，

转换后得： $\left[\begin{array}{c}{x}_0\\ y_0\\ z_0\end{array}\right]={\mathrm{kR}}^T\left[\begin{array}{c}{X-X}_c\\ {Y-Y}_c\\ {Z-Z}_c\end{array}\right]$

其中：

$R=\left[\begin{array}{ccc}{a}_1& a_2& a_3\\ b_1& b_2& b_3\\ c_1& c_2& c_3\end{array}\right]=\left[\begin{array}{ccc}\cos \mathrm{\φ}& 0& -\sin \mathrm{\φ}\\ 0& 1& 0\\ \sin \mathrm{\φ}& 0& \cos \mathrm{\φ}\end{array}\right]\left[\begin{array}{ccc}1& 0& 0\\ 0& \cos \mathrm{\ω}& -\sin \mathrm{\ω}\\ 0& \sin \mathrm{\ω}& \cos \mathrm{\ω}\end{array}\right]\left[\begin{array}{ccc}\cos \mathrm{\κ}& -\sin \mathrm{\κ}& 0\\ \sin \mathrm{\κ}& \cos \mathrm{\κ}& 0\\ 0& 0& 1\end{array}\right],$

$R=\left[\begin{array}{ccc}\cos \mathrm{\φ}\cos \mathrm{\κ}-\sin \mathrm{\φ}\sin \mathrm{\ω}\sin \mathrm{\κ}& -\cos \mathrm{\φ}\sin \mathrm{\κ}-\sin \mathrm{\φ}\sin \mathrm{\ω}\cos \mathrm{\κ}& -\sin \mathrm{\φ}\cos \mathrm{\ω}\\ \cos \mathrm{\ω\κ}\sin & \cos \mathrm{\ω\κ}\cos & -\sin \mathrm{\ω}\\ \sin \mathrm{\φ}\cos \mathrm{\κ}+\cos \mathrm{\φ}\sin \mathrm{\ω}\sin \mathrm{\κ}& -\sin \mathrm{\φ}\sin \mathrm{\κ}+\cos \mathrm{\φ}\sin \mathrm{\ω}\cos \mathrm{\κ}& \cos \mathrm{\φ}\cos \mathrm{\ω}\end{array}\right],$

$R^T=\left[\begin{array}{ccc}{a}_1& b_1& c_1\\ a_2& b_2& c_2\\ a_3& b_3& c_3\end{array}\right]=\left[\begin{array}{ccc}\cos \mathrm{\φ}\cos \mathrm{\κ}-\sin \mathrm{\φ}\sin \mathrm{\ω}\sin \mathrm{\κ}& \cos \mathrm{\ω}\sin \mathrm{\κ}& \sin \mathrm{\φ}\cos \mathrm{\κ}+\cos \mathrm{\φ}\sin \mathrm{\ω}\sin \mathrm{\κ}\\ -\cos \mathrm{\φ}\sin \mathrm{\κ}-\sin \mathrm{\φ}\sin \mathrm{\ω}\cos \mathrm{\κ}& \cos \mathrm{\φ}\cos \mathrm{\ω}& -\sin \mathrm{\φ}\sin \mathrm{\κ}+\cos \mathrm{\φ}\sin \mathrm{\ω}\cos \mathrm{\κ}\\ -\sin \mathrm{\φ}\cos \mathrm{\ω}& -\sin \mathrm{\ω}& \cos \mathrm{\ω}\cos \mathrm{\κ}\end{array}\right]$

其中a₁、a₂、a₃、b₁、b₂、b₃、c₁、c₂、c₃为旋转矩阵的元素；

将转换后的公式展开，得： $\begin{array}{c}{x}_0=k[a_1\left({X-X}_c\right)+b_1\left({Y-Y}_c\right)+c_1\left({Z-Z}_c\right)]\\ y_0=k[a_2\left({X-X}_c\right)+b_2\left({Y-Y}_c\right)+c_2\left({Z-Z}_c\right)]\\ z_0=k[a_3\left({X-X}_c\right)+b_3\left({Y-Y}_c\right)+c_3\left({Z-Z}_c\right)]\end{array}$

设： $\begin{array}{c}{X}_0=a_1\left({X-X}_c\right)+b_1\left({Y-Y}_c\right)+c_1\left({Z-Z}_c\right)\\ Y_0=a_2\left({X-X}_c\right)+b_2\left({Y-Y}_c\right)+c_2\left({Z-Z}_c\right)\\ Z_0=a_3\left({X-X}_c\right)+b_3\left({Y-Y}_c\right)+c_3\left({Z-Z}_c\right)\end{array},$

则： $\begin{array}{c}{x}_0={\mathrm{kX}}_0\\ y_0={\mathrm{kY}}_0\\ z_0={\mathrm{kZ}}_0\end{array},$

$\begin{array}{c}{x}_0/z_0=X_0/Z_0\\ y_0/z_0=Y_0/Z_0\end{array},$

又设： $\begin{array}{c}\mathrm{\α}=x_0/z_0=X_0/Z_0\\ \mathrm{\β}=y_0/z_0=Y_0/Z_0\end{array},$

其中α、β为观测值，其相应的改正数为ν_α、ν_β；将上式用泰勒公式展开成线性形式，即可得到线性化的误差方程如下：

$\left\{\begin{array}{c}{v}_{\mathrm{\α}}=\left(\mathrm{\α}\right)-\mathrm{\α}+\frac{\∂\mathrm{\α}}{{\∂X}_c}\mathrm{\ΔXc}+\frac{\∂\mathrm{\α}}{{\∂Y}_c}\mathrm{\ΔYc}+\frac{\∂\mathrm{\α}}{{\∂Z}_c}\mathrm{\ΔZc}+\frac{\∂\mathrm{\α}}{\∂\mathrm{\φ}}\mathrm{\Δ\φ}+\frac{\∂\mathrm{\α}}{\∂\mathrm{\ω}}\mathrm{\Δ\ω}+\frac{\∂\mathrm{\α}}{\∂\mathrm{\κ}}\mathrm{\Δ\κ}\\ v_{\mathrm{\β}}=\left(\mathrm{\β}\right)-\mathrm{\β}+\frac{\∂\mathrm{\β}}{{\∂X}_c}\mathrm{\ΔXc}+\frac{\∂\mathrm{\β}}{{\∂Y}_c}\mathrm{\ΔYc}+\frac{\∂\mathrm{\β}}{{\∂Z}_c}\mathrm{\ΔZc}+\frac{\∂\mathrm{\β}}{\∂\mathrm{\φ}}\mathrm{\Δ\φ}+\frac{\∂\mathrm{\β}}{\∂\mathrm{\ω}}+\mathrm{\Δ\ω}+\frac{\∂\mathrm{\β}}{\∂\mathrm{\κ}}\mathrm{\Δ\κ}\end{array}\right.,$

其中，(α)、(β)为使用待定参数的近似值带入公式中求出的α、β的计算值，ΔXc、ΔYc、ΔZc、Δω、Δκ为未知数Xc、Yc、Zc、φ，ω，κ的改正数；

然后采用以下公式计算α、β对各个未知数的偏导数：

${\mathrm{\α}}^{\′}={\left(\frac{X_0}{Z_0}\right)}^{\′}=\frac{X_0^{\′}{Z}_0-X_0{Z}_0^{\′}}{Z_0^2}=\frac{1}{Z_0}(X_0^{\′}-\mathrm{\α}{Z}_0^{\′})$

${\mathrm{\β}}^{\′}={\left(\frac{Y_0}{Z_0}\right)}^{\′}=\frac{Y_0^{\′}{Z}_0-Y_0{Z}_0^{\′}}{Z_0^2}=\frac{1}{Z_0}(Y_0^{\′}-\mathrm{\β}{Z}_0^{\′})$

再将误差方程转化为矩阵形式：V＝AH-L

其中，

V＝[v_α,v_β]^T

$A=\left[\begin{array}{c}\frac{\∂\mathrm{\α}}{{\∂X}_c},\frac{\∂\mathrm{\α}}{{\∂Y}_c},\frac{\∂\mathrm{\α}}{{\∂Z}_c},\frac{\∂\mathrm{\α}}{\∂\mathrm{\φ}},\frac{\∂\mathrm{\α}}{\∂\mathrm{\ω}},\frac{\∂\mathrm{\α}}{\∂\mathrm{\κ}}\\ \frac{\∂\mathrm{\β}}{{\∂X}_c},\frac{\∂\mathrm{\β}}{{\∂Y}_c},\frac{\∂\mathrm{\β}}{{\∂Z}_c},\frac{\∂\mathrm{\β}}{\∂\mathrm{\φ}},\frac{\∂\mathrm{\β}}{\∂\mathrm{\ω}},\frac{\∂\mathrm{\β}}{\∂\mathrm{\κ}}\end{array}\right]$

H＝[Δx,Δy,Δz,Δφ,Δω,Δκ]^T

L＝[α-(α),β-(β)]^T

根据最小二乘平差原理，列出法方程式：

A^TAH＝A^TL

则法方程的解为：

H＝(A^TA)^-1A^TL

从而求出待测未知数的近似值的改正数矩阵H的6元素ΔXc、ΔYc、ΔZc、Δω、Δκ；

最后对计算结果进行重复迭代计算，每次迭代时用未知数近似值与上次迭代计算的改正数之和作为新的近似值，重复计算过程，求出新的改正数，依此反复趋近，直至改正数小于设定的阈值，最后得到未知数的解，为：

X_c＝X_c ⁰+ΔX_c ¹+ΔX_c ²+Λ

Y_c＝Y_c ⁰+ΔY_c ¹+ΔY_c ²+Λ

Z_c＝Z_c ⁰+ΔZ_c ¹+ΔZ_c ²+Λ

φ＝φ⁰+Δφ¹+Δφ²+Λ

ω＝ω⁰+Δω¹+Δω²+Λ

κ＝κ⁰+Δκ¹+Δκ²+Λ

其中Xc⁰为Xc的近似值，ΔXc¹为第一次计算的改正数，ΔXc²为第二次计算的改正数，其他数据均依此类推。

4.根据权利要求1所述的球面全景影像与三维激光扫描点云的配准及测量方法，其特征在于：当有多个点同时对应到一个网格中时，选取多个点中的距离值最小者。