实现图像上像素点与实际位置的GPS对应

作者有话说

这篇随笔是基于我自己完成的一个项目，这个项目虽然看起来较为简单，但是由于我本身不是学这个方向的，因此在做的过程中还是遇到了一些大大小小的问题。经过仔细研究并多次调试代码，终于把这个问题的原理弄懂了。下面我将详细介绍该问题的解决过程，并在随笔末尾附上所有相关代码，希望有兴趣的可以一起交流学习。

目录

问题介绍  

无人机的出现给拍照录像带来了极大的便利。如下图所示，无人机停留在高为2m的高空中，其所装配的相机的拍摄方向与垂直方向的夹角为60°，利用相机参数求解所拍摄图像上所有像素点的大地坐标，实现图像上所有像素点与实际位置的GPS对应。另外，相机拍摄时的经纬度为（118.8675992，32.032575），相机的焦距为16mm，成像像素长和宽分别为4.65µm和3.9µm。

相关名词介绍

解决思路及过程

考虑到地球是椭球形状的，而经纬度是大地坐标，不方便计算，所以先利用高斯投影正算公式将大地平面坐标转换为高斯平面坐标，高斯平面坐标又和高度构成一个三维空间坐标。根据相机成像原理，把像素正中心坐标作为一个空间坐标，把相机中心坐标作为一个空间坐标，过这两点的直线与高斯平面的交点就是这个像素点所对应的实际位置。

下图反映了由图片像素点得到对应实际位置的过程，考虑到成像屏的尺寸相对于拍摄的实物来说非常小，所以我们把这两部分分开画图。图1和图2结合起来表示由像素点对应到实际位置的原理，即过像素点（彩色方格）与相机中心（空心圆点）的射线（彩线）与高斯平面的交点即为像素点对应的实际位置（射线终点）。

图1  光线从相机中心到像素点示意图

图2  光线从实际位置到相机中心示意图（斜视）

图3  光线从实际位置到相机中心示意图（俯视）

1.已知大地坐标$\left ( B,L \right )$及中央子午线经度$L_{0}$，计算高斯平面坐标$\left ( x,y \right )$，公式如下：

$ x=X+\frac{N}{2}sin\left ( B \right )cos\left ( B \right )l^{2}+\frac{N}{24}sin\left ( B \right )cos^{3}\left ( B \right )\left ( 5-t^{2}+9\eta ^{2}+4\eta ^{4} \right )l^{4}+\frac{N}{720}sin\left ( B \right )cos^{5}\left ( B \right )( 61-58t^{4}+270\eta ^{2}-330\eta ^{2}t^{2})l^{6}$

$y=Ncos\left ( B \right )l+\frac{N}{6}cos^{3}\left ( B \right )\left ( 1-t^{2}+\eta ^{2} \right )l^{3}+\frac{N}{120}cos^{5}\left ( B \right )\left ( 5-18t^{2}+t^{4}+14\eta ^{2}-58\eta ^{2}t^{2} \right )l^{5}$

其中，$B$为纬度，$l=L-L_{0}$，单位为弧度，$N=\frac{a}{\sqrt{1-e^{2}sin^{2}\left ( B \right )}}$，为卯酉圈曲率半径，$t=tan\left ( B \right )$，$\eta ^{2}=e^{2}cos^{2}\left ( B \right )$，$e=\frac{\sqrt{a^{2}-b^{2}}}{b}$为第二偏心率，$a$为旋转椭球长半轴，$b$为短半轴，$X$为子午线弧长。

2.已知高斯平面坐标$\left (x,y  \right )$及指定中央子午线经度$L_{0}$，计算大地坐标$\left ( B,L \right )$:

$B=B_{f}-\frac{t_{f}}{2M_{f}N_{f}}y^{2}+\frac{t_{f}}{24M_{f}N_{f}^{3}}\left (5+3t_{f}^{2}+\eta _{f}^{2}-9t_{f}^{2}\eta _{f}^{2}  \right )y^{4}+\frac{t_{f}}{720M_{f}N_{f}^{5}}\left ( 61+90t_{f}^{2}+45t_{f}^{4} \right )y^{6}$

$L=L_{0}+\frac{1}{N_{f}cos\left ( B_{f} \right )}y-\frac{1}{6N_{f}^{3}cos\left ( B_{f} \right )}\left ( 1+2t_{f}^{2}+\eta _{f}^{2} \right )y^{3}+\frac{1}{120N_{f}^{5}cos\left ( B_{f} \right )}\left ( 5+28t_{f}^{2}+24t_{f}^{4}+6\eta _{f}^{2}+8t_{f}^{2}\eta _{f}^{2} \right )y^{5}$

其中，$N_{f}=\frac{a}{\sqrt{1-e^{2}sin\left ( B_{f} \right )}}$，$M_{f}=\frac{a\left ( 1-e^{2} \right )}{\sqrt{\left (1-e^{2}sin\left ( B_{f} \right )  \right )^{3}}}$，$\eta _{f}^{2}=e^{2}cos^{2}\left ( B_{f} \right )$，$t_{f}=tan\left ( B_{f} \right )$，$B_{f}$为根据子午线弧长$X$反算的底点纬度。

一些值得注意的细节

保证一定的精度：因为像素尺寸非常小，他们所对应的实际位置的经纬度差别也是很小的，我刚处理的时候没有考虑到这个问题，所以导致得到的结果很多都是一样的（结果被截断了）。所以我后来将输出数据的位数设定为10为小数，才看出其差别。

结果

 经过编程计算，我们把得到经纬度数据分别根据像素大小画出一张曲面来观察结果的合理性，如下图（图4、图5）所示：

图4：经度                                         图5：纬度

根据图4和图5，可以看出经纬度变化的连续性和对称性，而且数值又徘徊在相机经纬度坐标周围，因此是合理的。

参考文献链接

matlab大地测量高斯投影正反算程序设计实验 - 百度文库https://wenku.baidu.com/view/121aedfdee06eff9aff807a3.html#

MATLAB代码

main.m

dms2degree.m

dms2rad.m

GaussianMapDirect.m

GaussianMapInverse.m

ImagePointCoordinates.m

 ImagePointMapping.m

latitude2meridian.m

meridian2latitude.m