摘要: 本文利用影子随着时间的变化，建立了太阳影子定位的模型。首先对视频中提取的图像进行裁剪、灰度化以及二值化的预处理，可以大大降低图像处理的计算量；其次利用透视变换，求出视频中实际的杆影长，对于受测量误差影响较大的异常数据，我们给予剔除；然后利用天文学中太阳运行轨迹相关理论，我们给出了太阳高度角，实际杆影长等的公式；最后，在已知杆长和未知杆长的情况下，把模型定位问题转化成求解经纬度和日期的非线性优化问题。通过比较，求解结果较好地吻合了实际的拍摄地点和日期。
Abstract:
With the length variation of pole’s shadow, we construct a localization model. Firstly, in order to reduce the work of images processing, we preprocess the extracted images, including cutting, gray processing, binarization processing, etc. Secondly, applying perspective transform, we get the actual length of pole’s shadow in video. For some images, measurement error results in some abnormal data of pole’s length, consequently, we remove them. Thirdly, using solar paths’ related theory in astronomy, we obtain formulas of the solar elevation angle and the shadow’s length. Lastly, whether the real length of pole in video is known or not, we transform the localization model into solving the longitude, latitude, and date of the nonlinear optimization problem. By comparison, the results of the solution are in good agreement with the actual shooting location and date.
1. 引言
随着数据业务和多媒体业务的扩展，人们对定位的需求日益增大，从最开始15世纪的路基无线电导航到现代的GPS定位、超声波定位和蓝牙定位等，定位技术获得了快速发展。为弥补传统定位技术在信号接收、网络环境和设备性能等方面的局限性 [1] ，本文结合计算机视觉 [2] 这一研究热点，建立了利用太阳影子变化的定位模型。通过任意一段太阳下拍摄的视频，太阳影子定位模型能利用影子的变化推算出视频拍摄地点与日期。目前国内对太阳影子定位技术的研究与记载甚少，2015年全国大学生建模竞赛A题涉及相关思路，但研究者的定位模型在图像识别和模型建立方面都存在不足，例如胡毅华 [3] 、黄小卉 [4] 等测量影长时未对视频图像做透视变换图像校正，朱海洋 [5] 建立模型时忽略地球绕日的椭圆轨道影响，降低了模型精度。
本文首先利用Matlab软件提取视频中的图像，经过裁剪、灰度化以及二值化降低运算量后，利用透视变换校正图像，从而消除了照相机拍摄带来的视觉偏差，随后基于Hough变换的直线检测技术测算出视频中不同时刻下杆子的影长；然后利用影长与赤纬角、纬度和时角等参数的理化关系建立非线性优化模型的过程中，本文先从已知视频里杆长的情况出发，求解出视频拍摄的地点和日期，其次为提高模型的适用性，将原模型拓展到杆长未知的情况下，建立新的模型同样能求解出视频拍摄的地点与日期；最后讨论了拍摄地的大气折射率对结果的影响。
2. 透视投影校正
2.1. 图像预处理
本文以2015年全国大学生建模竞赛A题提供视频为例 [6] ，将时长为40分钟的视频导入Matlab软件，共61,022帧。若提取图像时间间隔过短时，提取的相邻图像中影子长度变化将不明显，同时计算复杂度变大，导致难以定位，而提取图像时间间隔过长时，样本数据变少，同样难以定位，因此提取图像时间间隔不宜过短或过长。本文每隔1525帧即1分钟提取视频图像，得到41张图像。由于提取的图像里面包含了大量的无效信息，我们首先对图像进行剪裁，只保留有效信息，见图1。
其次为了降低图像处理的计算量，我们对剪裁的图像先后进行灰度化和二值化处理。依据人眼对不同颜色的敏感度，将彩色图像的R、G、B三个颜色分量进行加权平均求和，可得到灰度化后的图像。本文采用Matlab中的rgb2函数进行图像灰度化 [7] ，其默认的亮度方程为
(1)
其中
，
，
为图像各像素点的红、绿、蓝像素值。
然后对灰度化后的图像进行二值化处理，我们选择普遍采用的最大类间法 [8] 。通过最大类间法，我们确定了一个二值化的阈值
，当图像上的灰度值大于等于该阈值
时，像素点灰度设置为255，否则为0。该方法能在不影响图像整体和局部分布的基础上呈现明显的黑白效果。最后利用形态学的开、闭操作剔除一些杂点和去除小面积形状，最终得到图2。
2.2. 利用透视变换计算影长
以降噪后图像的左上角为坐标原点建立平面直角坐标系，利用Hough变换 [9] ，检测出降噪图像中的杆影所在线段及端点坐标，由于照相机拍摄得到的物像存在视觉透视，因此不能直接根据降噪后的图像计算杆影长。降噪后的图像中，Matlab能识别出杆子底座的4个顶点及杆末端的坐标分别为A(14,137)，B(20,113)，C(123,113)，D(122,137)和E(70,125)，见图3。已知底座底面是正方形，AD平行于横坐标轴，
故在透视变化下，线段AD长度不变，从而得到A’(14,137)，B’(14,29)，C’(122,29)以及D’(122,137)。利用四对透视变换前后对坐标即可求得透视变换矩阵 [10] ，本文直接调用Matlab的cp2tform函数 [11] 计算出透视变换的转换矩阵Tform，进而求出不同时刻透视变换后的E'和F'的坐标，然后计算E’，F’的欧式距离便得到所需要的实际杆影长。
2.3. 影长数据处理
经过透视变换后，我们可得到41张图片中的实际杆影长。在视频时间段内，杆影长是逐渐减少的，
如图4所示。因为短期内影长变化微小扩大了误差的影响，导致图4中存在部分异常点，为了保证模型的可靠性，我们需要剔除这些异常点。最终保留41个影长中的前10个与后12个，共22个数据作为有效杆影长数据，用于后续的非线性优化求解。
3. 模型建立与求解
3.1. 杆长已知日期未知求解地点
在现实生活中要在视频拍摄地寻找一个已知高度的物体并非难事，并且通常能轻易获知当时日期，部分相机拍摄时也会在视频上体现拍摄日期，但是也会经常遇到搜集的视频没有拍摄时间记载的情况，因此建立一个杆长已知，拍摄日期未知，求解拍摄地经纬度的模型具有一定现实意义。本例视频中杆长为2米，拍摄日期为7月13日，我们在假设日期未知的情况下建立模型，模型的目标函数和约束条件需涉及太阳高度角、赤纬角、时角、时差、本地太阳时等相关天文物理参数。
首先给出太阳高度角的公式 [12] ：
(2)
其中
为高度角，
为赤纬角，
为纬度，
为太阳时角，而赤纬角
近似公式 [13] 为：
(3)
其中
为自1月1日算起的累计日。
时角
的计算公式 [14] 为：
Figure 4. Change of shadow length with time
图4. 影长随时间的变化
(4)
其中
为本地太阳时，
为正表示偏东，
为负表示偏西。通常
可用本地时
代替，但本地时与本地太阳时因为地球自转和人为调整等因素存在差异，故本文利用调整公式计算精度更高的
：
(5)
上式中本地时
可由视频中提供的拍摄时对应的北京时间
和拍摄地经度
转换得到
(6)
而时差
的近似公式 [15] 为：
(7)
(8)
根据太阳影长和折射后的太阳高度角的几何关系可以得到太阳影子的长度
为
(9)
在已知杆长
为2米，得到实际的太阳影长
是关于本地时
，经度
，纬度
和累积日
的函数，即
(10)
因此模型定位问题可转化为求解非线性优化问题，其中目标函数为：
(11)
其中
为论文2.3部分得到的不同时刻
对应的实际杆影长，共22个数据，
为由公式(10)得到的理论值，即：
(12)
而关于
，
，
的约束条件为：
使用变步长搜索法、拟牛顿法以及lingo软件都可以解出以上非线性优化问题，最终得最优解为：
，
即拍摄地点为内蒙古鄂尔多斯(E109.9727˚, N39.0051˚)，比较吻合实际拍摄地点，而对应日期为7月12日，与真正的日期相差一天，关于日期结果的相对误差为
，也比较接近拍摄日期。结果的误差主要来自于剔除部分异常的杆影长数据，测量有一定的误差等，这些最终都对结果有稍微的影响。
3.2. 杆长和日期未知求地点
当拍摄的视频没有杆长的信息，或者视频从他处获得，持有者对视频中对象了解较少时，需要重新建立模型。上述情况对应于杆长和时间信息未知，确定拍摄地点和时间的模型。
由太阳高度角公式(2)可计算出正切值
，也是关于本地时
，经度
，纬度
和累积日
的函数，即：
(13)
函数实际上应该等于杆长与影长之比，由本文2部分，利用透视变换我们可以求得实际杆长，以及筛选得到的22个杆影长数据，于是便得到相应的22个的
。由此杆长和日期未知的模型定位问题便转化为下列非线性优化模型，其中目标函数为：
(14)
约束条件为：
求解得到：
这与实际拍摄地点内蒙古自治区鄂尔多斯的经纬度非常吻合，而求解的对应日期仍然为7月12日，与实际的日期也比价接近。与3.1部分类似，剔除部分异常的杆影长数据，测量本身存在一定的误差等导致对结果有稍微的影响。
4. 模型的进一步讨论
为了进一步提高模型精度，考虑大气折射率
对太阳高度角的影响，原太阳高度角为
，经过折射后，太阳高度角为角
，如图5。
若设大气层的折射率为
，则可得：
(15)
而大气层折射率近似的计算公式 [16] 为：
(16)
其中
为该点的大气压，
为温度，
为相对温度。根据中国天气网，鄂尔多斯在2015年7月13日气温为18℃~29℃，视频拍摄时为8点54分到9点34分，
取20℃，气压为895.6 hpa，
为37%，得
，说明当天鄂尔多斯地区大气折射作用非常弱，可以忽略不计。
5. 结论和评价
本文利用影子随着时间的变化，建立了太阳影子定位的模型。在已知杆长和未知杆长的情况下，把模型定位问题转化成求解拍摄地点的经纬度和拍摄日期的非线性优化问题。通过比较，求解的结果较好地吻合了实际的拍摄地点和日期。
在模型中，视频图像是否经过透视校正，本地太阳时与时差的模拟计算公式选择以及是否考虑不同地区大气折射率的差异都将对模型结果的精确程度产生影响，本文综合考虑了上述影响因素，得到了较为理想的结果。另外测量误差，剔除异常数据，求解优化模型等环节也不可避免地会影响最终的结果。
基金项目
北京高等学校“青年英才计划”(YETP0769)，国家自然科学基金(61571002, 61370193)。