专利名称：基于视频中太阳影子轨迹的经纬度估计方法
技术领域：
本发明属于图像处理和地理信息系统技术领域，具体涉及一种经纬度估计方法。
背景技术：
经纬度的估计在军事和民用中都有广泛的应用，在法医学，安全领域，以及航海等领域也有着巨大的应用潜力。本发明使用未进行相机校准的视频中的影子轨迹来估计照片采集位置的经纬度信息。其输入是用固定的针孔相机拍摄的含有真实世界垂直关系且在同一场景中含有两个影子轨迹的视频帧，如图1所示。本发明技术的输出是采集照片位置的经纬度信息。
本发明中使用到的背景技术有(1)日晷的原理。投影日晷1是观测日影记时的仪器。投影日晷不设置指时针，仅在地平面依地理纬度的不同绘制不同扁率的椭圆，在其上刻划时间线，并将长轴指向正东西方向，南北向的短轴上则需刻上日期，指示立竿测量时刻的正确方向。(2)Fei Lu和操晓春等2提出的从影子轨迹恢复出地平线的原理。
目前基于图像的经纬度估计方法主要有James Hays和Alexei A.Efros2提出的基于单个图像的经纬度估计方法。该技术采用了一种纯数据驱动的场景匹配技术。首先建立一个有地理标志的图像信息数据库，对于测试图像集手动的移除黑白图像，艺术图像，受噪音影响特别大的图像和涉及隐私等的图像，然后运用颜色直方图，纹理直方图，线性特征等技术进行场景匹配，最后基于数据驱动的地理信息估计。该方法需要超过600百万张有地理标记的GPS图像，收集比较困难，而且这600百万不能代表整个地球上的图像样本。该方法计算复杂度较高，在允许200km以内定位精度时只有约16％的查全率和查准率。Sunkavalli4提出一种基于颜色变化的经纬度估计方法。该技术通过对场景进行光度分析，间接的得到太阳位置，进而在不需要卫星数据和太阳在相机视场内不可见的情况下，通过使用一般几何限制得到经纬度信息。但是这种方法限制条件很多，例如要求线性传感器，要求朗伯平面，不考虑屏幕像素点上光的相互作用，仅仅考虑太阳光和环境光，而且建模复杂，误差较大。
主要参考文献 1http://en.wikipedia.org/wiki/Sundial. 2Fei Lu and Xiaochun Cao.Camera Calibration from Two Shadow Trajectories.Proceedings of the18th International Conference on Pattern Recognition，2006.(Fei Lu和操晓春，基于两条影子轨迹的相机校准，2006年第18届模式识别国际会议) 3James Hays and Alexei A.Efros.Estimating geographic information from a single image.ProcIEEE Int.Conf.on Computer Vision and Pattern Recognition(CVPR)，Florida，USA，June，2008.(James Hays和Alexei A.Efros，基于单张图片的地理信息估计，2008年美国福罗里达州，IEEE计算机视觉与模式识别国际会议) 4Kalyan Sunkavalli and Fabiano Romeiro.What do color changes reveal about an outdoor scene.Proc IEEE Int.Conf.on Computer Vision and Pattern Recognition(CVPR)，Florida，USA，June，2008.(Kalyan Sunkavalli和Fabiano Romeiro.户外场景的颜色变化可以揭示什么，2008年美国福罗里达州，IEEE计算机视觉与模式识别国际会议)
发明内容
本发明的目的在于克服现有技术的上述不足，提供一种计算复杂度较低，精度较高，是一种利用自然图像序列或视频帧实现经纬度估计的方法，利用该方法，能够根据未经过校准的影子的位置，得出所拍摄图像的经纬度信息。
为此，本发明采用如下的技术方案 一种基于视频中太阳影子轨迹的经纬度估计方法，包括下列步骤 (1)获取同一场景中含有至少两个影子轨迹的自然图像序列或视频帧； (2)对每一帧图像，检测影子轨迹； (3)同一个影子轨迹里不同的两个点的所确定的直线与对应时刻另一个影子轨迹中的两个点所确定的直线的交点即是一个灭点，计算灭点； (4)由各个灭点拟合出地平线； (5)在图像中沿着有垂直关系的物体画出两条直线，计算出这两条直线分别与地平线的交点坐标，这两个坐标就是两个互相垂直的灭点vx，vy坐标； (6)根据公式计算焦距f，式中，u0，v0为相机主点即图像中心点的坐标，根据公式计算ω； (7)在地平线上找到另一对满足限制的点vp，vvertical； (8)按照下列方法实现从图像到经过度量纠正的世界坐标的转换 1)任意取两个点M，N，直线vxM，vyM，vpN，vverticalN的斜率分别记做k1，k2，k3，k4。
2)以
为半径，以
为圆心得到一个圆其中
为直线vxM，vyM的夹角； 3)以
为半径，以
为圆心得到另一个圆其中
为直线vpN，vverticalN的夹角； 4)计算上述两个圆的交点α，β。
5)计算H＝AP，其中，(l1，l2，l3是地平线L∞的三个分量)， (9)根据公式X＝Hx还原出经过度量纠正的世界坐标，其中X是经过度量纠正的世界坐标中的点，x是图像坐标中的点； (10)拟合影子轨迹； (11)基于日晷原理，计算纬度φ； (12)按照下列方法计算经度以太阳时间12点所对应的视频帧为中心，在步骤(10)所拟合的影子轨迹上寻找到极值点，并且记录下该点所对应的时刻t，通过以下公式来计算时间差值TT＝|t-12|+teot-tdst，其中，teot是时差，tdst是夏令时，如果t大于12点，那么经度L为L＝Lmeridian-T/4，如果如果t小于12点，那么经度L为L＝Lmeridian+T/4，其中，Lmeridian是该时区中央经线。
作为优选实施方式，上述的步骤(2)中运用基于最大流算法的图分割方法，得到影子区域的二进制掩膜，对这个区域求主轴，主轴与二进制掩膜边缘的交点即是影子的轨迹点；步骤(4)中采用随机采样一致性算法拟合地平线；步骤(11)中按照下列方法计算纬度φ设日晷上X轴指向东西方向，Y轴指向南北方向，水平距离X和垂直距离Y分别定义为X＝sinh和
其中，h是时角，其大小为h＝(T24-12)×15，T24是以24小时表示的标准时间，φ是所要估算的纬度值，根据方程计算φ，其中，δ是太阳的倾角，A是太阳的方位角。
本发明提供一种自然图像序列或视频帧中基于影子轨迹的经纬度估计方法，从影子的轨迹来恢复出地平线，进而用恢复出来的地平线来消除透视畸变。本发明采用逆转日晷设计，输入是未经过校准的影子的位置，输出是所拍摄图像的经纬度信息。Sunkavalli的方法的输入之一要求提供三个互相正交的灭点，而本发明的方法从影子的轨迹本身恢复出几何限制，而且本发明仅仅需要采集影子的轨迹。
本发明将经纬度的估计应用到视频帧中。仅仅从影子的轨迹来恢复出地平线，进而用恢复出来的地平线来消除透视变形，并不要求投影的物体在图像中可见。本发明相比James Hays和Alexei A.Efros提出的基于单个图像的经纬度估计方法，不需要收集大量的有地理标记的GPS图像。相比Sunkavalli提出的基于颜色变化的经纬度估计方法，不需要提供三个互相正交的灭点，不需要非线性最小化，因此计算复杂度较低，但却有较高的精度，且比他们的方法简单。
图1含有影子轨迹的视频帧序列。
图2本发明的经纬度估计方法的总流程图。
图3标注了经过追踪得到的不同时间的影子点的图像，图中，白色实线为拟合的两个物体阴影的主轴，白色的星号标记为计算的影子点。
图4灭点的计算原理图。
图5(a)为经过随机采样一致性算法拟合出来的地平线，其上的黑色十字号表示的点为灭点。
图5(b)标注了影子点的视频图像，图中，黑色的点为图像中观测的影子点。
图6通过日晷计算经度原理图。
图7标注了随时间不断变化的影子轨迹点以拟合出来的二次曲线的图像，图中白色的十字号表示的是随时间不断变化的影子轨迹点，白色的曲线是由轨迹点拟合出来的二次曲线。
具体实施例方式 本发明的技术方案如下首先获取自然图像序列或者视频帧，并且对每一帧图像运用基于最大流(MAX FLOW)算法的图分割(graph cut)技术检测出阴影区域的二进制掩膜，然后找到影子区域二进制掩膜的主轴，主轴与二进制掩膜区域边缘的交点，就是所要求的影子的轨迹点。
然后，同一个影子轨迹里不同的两个点的所确定的直线与对应时刻另一个影子轨迹中的两个点所确定的直线的交点就是一个灭点，这样可以确定多个灭点。运用随机采样一致性算法，剔除掉异常值点后拟合出地平线。拟合互相垂直的灭点vx，vy，计算出仿射纠正和投影纠正矩阵从而得到只经过相似变换的单应矩阵。进而可以还原出经过度量纠正的世界坐标。
最后，利用最小二乘法实现二次曲线的拟合，拟合出经过度量纠正世界坐标中的影子点的轨迹，并且参考日晷设计，利用相似关系估计出纬度。利用二次曲线的极值点计算时差，从而有效地恢复出所拍摄图像的经纬度信息。
下面结合附图和实施例对本发明做进一步详述。
参见图1，本发明的基于影子轨迹的经纬度估计方法，包括下列步骤 步骤1获取视频帧 获取同一场景中含有至少两个影子轨迹的自然图像序列或者同一场景中含有至少两个影子轨迹的视频帧。图1为含有影子轨迹的视频帧序列。
步骤2检测影子轨迹 本实施例使用如下半自动的检测方法对每一帧图像，运用基于最大流(MAXFLOW)算法的graph cut(图分割)技术，得到影子区域的二进制掩膜，对这个区域求主轴，主轴与二进制掩膜边缘的交点就是所要的影子的轨迹点。还可以采用现有的图像处理技术中常见的基于阈值的图像分割方法检测影子轨迹。本实施例采用全局最优的方法检测影子轨迹，因而检测效果要更好些。
图3为经过追踪得到的不同时间的影子点。图中，白色实线为拟合的两个物体阴影的主轴，白色的星号标记为计算的影子点。
步骤3计算灭点 由于太阳光是平行的，同一个影子轨迹里不同的两个点的所确定的直线与对应时刻另一个影子轨迹中的两个点所确定的直线的交点就是一个灭点。图4为灭点的计算原理图。图中，3D点ti在不同时刻j在地面上的投影为Sij。
如图4，影子轨迹1上任何两个时刻所收集到的影子点所确定的直线与影子轨迹2上对应时刻所收集到的影子点所确定的直线是平行的，因此我们可以求得灭点。原理是这样的图4中ti是世界坐标中的点，bi是ti距离地面最近的点，即tibi与地面垂直。而Sij表示不同时刻j影子点的位置。于是在3D世界坐标中有由于太阳距离地球大约1.52*1012米，我们可以假设太阳光是平行的，故b1S11和b2S21平行，b1S12和b2S22是平行的。垂直地面的两个物体t1b1和t2b2也是平行的。因此可以得到Δb1s11s12和Δb2s21s22相似，进而有直线S11S12和S21S22平行。可见bi在此步的计算当中并不需要，因此，在图像中不要求物体一定可见，并且同一个影子轨迹里不同的两个点所确定的直线与对应时刻另一个影子轨迹中的两个点所确定的直线的交点就是一个灭点。
步骤4计算地平线L∞ 假设每条影子的轨迹上都有n个时刻对应的n个点，那么用步骤3我们可以求出Cn2个灭点，但是由于噪音的影响这些点中可能有异常值点，会影响拟合出来的地平线。因此，我们对所得到的Cn2个灭点运用随机采样一致性算法，剔除掉异常值点后再拟合地平线。其中随机采样一致性算法的思想如下 1.随机的选择2个点。
2.实例化一条直线。
3.计算出满足直线模型的内点的个数，记为ni。
4.如果ni的值大于阈值，则程序终止，返回所得到的模型。阈值通常是依据经验值取得。在本发明中，我们阈值取的是5。
5.重复N次测试，直到所得到的模型的ni的值最大，N可以灵活定义。本发明中N取的是100。
因此从随着时间变化的影子轨迹中我们可以还原出由各个灭点所拟合出的直线，即地平线L∞。计算结果如图5(a)所示。
图5(b)为标注了影子点的视频图像，图中，黑色的点为图像中观测的影子点。
步骤5计算灭点vx，vy 3D世界中有垂直关系的物体在图像中是随处可见的，例如图像中的两面互相垂直的墙，图像中与地面垂直的杆等。在图像中沿着在3D世界中有垂直关系的物体画出两条直线。计算出这两条直线分别与地平线的交点坐标，这两个坐标就是两个互相垂直的灭点vx，vy。
步骤6计算焦距f 步骤5计算的vx，vy为一对正交灭点，对绝对二次曲线的图像ω提供了一个线性的约束由于ω＝K-TK-1，其中为相机内部参数矩阵，f为相机焦距，λ为纵横比，u0，v0为相机主点即图像中心点的坐标，γ为偏斜系数。通常情况下，常用相机内部参数K的纵横比λ取值为1，偏斜系数γ为0。绝对二次曲线的图像ω只依赖于变量f。因此，通过vx，vy我们可以计算焦距f和绝对二次曲线的图像ω如下 步骤7在地平线上找到另一对满足限制的点 找到这样点的求解方式有很多，其中一种求解方式如下 1.求解第一个点。因为vx，vy都在地平线上，所以满足如下线性约束 由(1)(2)可以得出 令vp＝vx+vy，由(3)式可知vp是在地平线上的。vp为所求的第一个点。
2.求解满足约束的第二个点。首先求出一条与vp垂直的直线Lp，Lp满足约束Lp＝ωvp。Lp与地平线L∞的交点就是另一个满足条件的点vvertical。
步骤8从图像到经过度量纠正的世界坐标的转换 1.任意取两个点M，N。直线vxM，vyM，vpN，vverticalN的斜率分别记做k1，k2，k3，k4。
2.以
为半径，以
为圆心可以在虚平面内得到一个圆。其中
为直线vxM，vyM的夹角，本发明中为90°。第一个圆的方程为 3.以
为半径，以
为圆心可以在虚平面内得到另一个圆。其中
为直线vpN，vverticalN的夹角，本发明中为90°。第二个圆的方程为 4.解两个方程(1)(2)得到虚平面内两个圆的交点α，β。
5.H＝AP。
其中，(l1，l2，l3是地平线L∞的三个分量)， 步骤9还原出经过度量纠正的世界坐标 根据公式X＝Hx，其中X是经过度量纠正的世界坐标中的点，x是图像坐标中的点。
步骤10利用最小二乘法实现影子轨迹的拟合以及纬度的估计 参见图6，图中的F点是静态物体距离地面最近的点即脚点，S是影子轨迹点，P是影子轨迹点在过F点与X轴平行的直线上的投影点。显然，ΔFSP和ΔFMN相似。于是根据三角形相似的比例关系可以推导出经度。
根据天文知识，一天当中静态物体影子轨迹近似一个二次曲线。利用最小二乘法拟合出影子轨迹的二次曲线。而投影日晷通过观测投影方向来计时，其影子轨迹是一个椭圆。如图6所示通过日晷计算经度原理图日晷上X轴指向东西方向，Y轴指向南北方向。水平距离X和垂直距离Y分别定义为 X＝sinh(8)
其中，h是时角，其大小为h＝(T24-12)×15，T24是以24小时表示的标准时间。φ是所要估算的纬度值。投影物体脚点的位置跟投影日晷中心的关系为E＝tanδcosφ，其中δ是太阳的倾角。由于ΔFSP和ΔFMN相似，所以可以得 其中A是太阳的方位角。式子(10)中，只有一个未知数φ，它就是所要估算的纬度值。
步骤11计算经度 以太阳时间12点所对应的视频帧为中心，在步骤10所拟合的二次曲线上找到极值点，并且记录下该点所对应的时刻t。因此时间差值T可以通过以下公式来计算 T＝|t-12|+teot-tdst(11) 其中，teot是时差，tdst是夏令时。由于经度每隔15度，地方时相差1小时。因此，如果t大于12点，那么经度L为 L＝Lmeridian-T/4(12) 如果如果t小于12点，那么经度L为 L＝Lmeridian+T/4(13) 其中，Lmeridian是该时区中央经线。
图7标注了随时间不断变化的影子轨迹点以拟合出来的二次曲线的图像，图中白色的十字号表示的是随时间不断变化的影子轨迹点，白色的曲线是由轨迹点拟合出来的二次曲线。
权利要求
1.一种基于视频中太阳影子轨迹的经纬度估计方法，包括下列步骤
(1)获取同一场景中含有至少两个影子轨迹的自然图像序列或视频帧；
(2)对每一帧图像，检测影子轨迹；
(3)同一个影子轨迹里不同的两个点的所确定的直线与对应时刻另一个影子轨迹中的两个点所确定的直线的交点即是一个灭点，计算灭点；
(4)由各个灭点拟合出地平线；
(5)在图像中沿着有垂直关系的物体画出两条直线，计算出这两条直线分别与地平线的交点坐标，这两个坐标就是两个互相垂直的灭点vx，vy坐标；
(6)根据公式计算焦距f，式中，u0，v0为相机主点即图像中心点的坐标，根据公式计算ω；
(7)在地平线上找到另一对满足限制的点vp，vvertical；
(8)按照下列方法实现从图像到经过度量纠正的世界坐标的转换
1)任意取两个点M，N，直线vxM，vyM，vpN，vverticalN的斜率分别记做k1，k2，k3，k4。
2)以
为半径，以
为圆心得到一个圆其中
为直线vxM，vyM的夹角；
3)以
为半径，以
为圆心得到另一个圆其中
为直线vpN，vverticalN的夹角；
4)计算上述两个圆的交点α，β。
5)计算H＝AP，其中，(l1，l2，l3是地平线L∞的三个分量)，
(9)根据公式X＝Hx还原出经过度量纠正的世界坐标，其中X是经过度量纠正的世界坐标中的点，x是图像坐标中的点；
(10)拟合影子轨迹；
(11)基于日晷原理，计算纬度φ；
(12)按照下列方法计算经度以太阳时间12点所对应的视频帧为中心，在步骤(10)所拟合的影子轨迹上寻找到极值点，并且记录下该点所对应的时刻t，通过以下公式来计算时间差值TT＝|t-12|+teot-tdst，其中，teot是时差，tdst是夏令时，如果t大于12点，那么经度L为L＝Lmeridian-T/4，如果如果t小于12点，那么经度L为L＝Lmeridian+T/4，其中，Lmeridian是该时区中央经线。
2.根据权利要求1所述的基于视频中太阳影子轨迹的经纬度估计方法，其特征在于，步骤(2)中运用基于最大流算法的图分割方法，得到影子区域的二进制掩膜，对这个区域求主轴，主轴与二进制掩膜边缘的交点即是影子的轨迹点。
3.根据权利要求1所述的基于视频中太阳影子轨迹的经纬度估计方法，其特征在于，步骤(4)中采用随机采样一致性算法拟合地平线。
4.根据权利要求1所述的基于视频中太阳影子轨迹的经纬度估计方法，其特征在于，步骤(11)中按照下列方法计算纬度φ设日晷上X轴指向东西方向，Y轴指向南北方向，水平距离X和垂直距离Y分别定义为X＝sinh和
其中，h是时角，其大小为h＝(T24-12)×15，T24是以24小时表示的标准时间，φ是所要估算的纬度值，根据方程计算φ，其中，δ是太阳的倾角，A是太阳的方位角。
全文摘要
本发明属于图像处理和地理信息系统技术领域，涉及一种基于视频中太阳影子轨迹的经纬度估计方法先获取自然图像序列或者视频帧，并且对每一帧图像检测出影子的轨迹点；然后确定多个灭点，并拟合出地平线；拟合互相垂直的灭点，计算出仿射纠正和投影纠正矩阵；进而还原出经过度量纠正的世界坐标；再拟合出经过度量纠正世界坐标中的影子点的轨迹，并且参考日晷设计，利用相似关系估计出纬度；利用二次曲线的极值点计算时差，从而有效地恢复出所拍摄图像的经纬度信息。本发明复杂度较低，精度较高，是一种利用自然图像序列或视频帧实现经纬度估计的方法，能够根据未经过校准的影子的位置，得出所拍摄图像的经纬度信息。
文档编号G01C1/00GK101493322SQ20091006781
公开日2009年7月29日 申请日期2009年2月4日 优先权日2009年2月4日
发明者操晓春, 曲彦龄, 孙济洲, 琳 武, 郭晓杰, 炜 张 申请人:天津大学