太阳影子定位
摘 要
本文考虑了影响直杆影长的五个因素——直杆所在地经度、纬度、日期（太阳积日）、时间、直杆高度五个参数，引入太阳高度角、时角、太阳方位角、赤纬四个辅助变量，利用三角函数关系分析了五个参数与影长变化的关系，建立出影长变化的数学模型，并在两个或者三个变量未知的情况下，基于最小二乘准则，采取网格搜索的方法，逐步逼近求得最优解。 问题一，为建立影长l 变化的数学模型，需要求出模型中的直杆所在地经度λ、纬度ϕ、日期（太阳积日）N 、时间t 和直杆高度H 五个参数与影长l 的关系l =f (λ, ϕ, N , t , H ) 。我们引入了太阳高度角这一辅助量，利用太阳高度角θ和经度λ、纬度ϕ和时间t 的关系和太阳高度角与影长l 的关系，建立出影长变化的数学模型l =f (λ, ϕ, N , t , H ) （见第5页公式5）。模型建立后，选择其中一个参数作为变量，固定其他四个参数，代入模型中，即可求得影长l 关于各参数的变化规律，所求规律具体可见正文部分第5-7页。问题中还需要2015年10月22日北京时间9:00-15:00之间天安门广场（北纬39度54分26秒, 东经116度23分29秒）3米高的直杆的太阳影子长度的变化曲线，将题目中的参数代入即可求出在该参数条件下的影长变化曲线，结果见第7页图5。
问题二，附件1给出21组影子顶点坐标x 、y 和日期N 、时间t 数据。直杆高度H 和影子顶点坐标所在的坐标系未知，我们要求出直杆所处地点地理经度λ、纬度ϕ。首先将顶点坐标x 、y 数据转化成影长l 数据，然后基于最小二乘准则，采用网格搜索的方式求出最优解。为缩小搜索范围，在网格搜索前采用二次拟合加时角公式的方法估算出经度大致值。求取的最终结果是H =1. 8m ，λ=111︒01'25"E ，ϕ=19︒57'30"N ，地理位置大约在海南省海口市附近。
问题三，与题目二不同的是，日期N 这一变量也是未知的，因此可以直接套用问题二的解决方法，只是在网格搜索的时候增加一个日期变量。附件2所求结果是H =1. 9m ，λ=78 01'23''E ，ϕ=40 01'15''N ，N =6月4日，7月28日，在新疆乌鲁木齐附近；附
件3所求结果是H =3. 0m ，λ=10957'02''E ，ϕ=2710'42''N ，N =1月14日，12月
16日，在湖南省长沙市附近。
问题四，已知视频数据和日期N 、时间t 、直杆高度H ，我们要求的是直杆所在地的经度λ、纬度ϕ。首先将视频等时间间隔截取成21张图片，使用CorelDRAW 软件对这21张图片进行预处理，手动读取出21组相关点坐标数据，将坐标数据转化成影长数据，并对这些数据进行校正。在已知直杆高度H 的条件下，采取与问题二相同的解决方法求解。最
终结果是λ=10900'39''E ，ϕ=4155'39''N ，在内蒙古包头市附近。如果拍摄日期未知，
我们认为可以求解，具体的分析过程见17-18页。
关键字 影子定位 顶点坐标 网格搜索 太阳高度角 最小二乘准则
一、问题提出
在大数据时代，随着海量数据的出现，人们需要在短暂的时间内从大量的数据中提取出有效的信息，因此视频数据分析技术也越来越重要。如何确定视频的拍摄地点和拍摄日期是视频数据分析的重要方面，这在生产生活的各个领域都有着极为广泛的应用。太阳影子定位技术就是通过分析视频中物体的影子变化，确定视频拍摄的地点和日期的一种方法。太阳每天东升西落，世间万物影子也随之变化，其中影子的长度、方位与记录时间、物体地理坐标、物体本身长度之间存在密不可分的关系。因此，研究它们之间的变化规律成为太阳影子定位技术的核心问题，更确切的说，通过影子长度或者影子顶点坐标的变化确定拍摄日期和地点是此次建模需要解决的核心问题。需要解决的具体问题有以下几个：
1. 建立影子长度变化的数学模型，分析影子长度关于各个参数的变化规律，并应用建立的模型画出2015年10月22日北京时间9:00-15:00之间天安门广场（北纬39度54分26秒, 东经116度23分29秒）3米高的直杆的太阳影子长度的变化曲线。
2. 根据某固定直杆在水平地面上的太阳影子顶点坐标数据，建立数学模型确定直杆所处的地点。将模型应用于附件1的影子顶点坐标数据，给出若干个可能的地点。
3. 根据某固定直杆在水平地面上的太阳影子顶点坐标数据，建立数学模型确定直杆所处的地点和日期。将模型分别应用于附件2和附件3的影子顶点坐标数据，给出若干个可能的地点与日期。
4．附件4提供了一段某直杆在太阳下影子变化的视频，并且已通过某种方式估计出直杆的高度为2米。请建立确定视频拍摄地点的数学模型，并应用模型给出若干个可能的拍摄地点。
如果拍摄日期未知，能否根据视频确定出拍摄地点与日期？
几个问题环环相扣，密不可分，在求解过程中，通过建立影子变化的数学模型，解决根据视频中影子顶点坐标变化或者影子长度变化来确定拍摄地点和日期的问题。
二、问题分析
由题意可知，该问题主要是根据视频中直杆影子变化的数据来确定拍摄日期和地点，从而达到太阳影子定位的目的。事实上，影子变化与时间、物体地理经度、纬度和物体本身高度息息相关，但是又无法找到它们之间的直接关系，因此将它们联系起来，并建立模型是解决问题的核心。此外，如何比较准确地从视频中读取出影子长度也是本次建模的难点所在。 关于问题一，主要问题是建立影子长度l 变化的模型，这也是本次建模的核心。事实上，模型中涉及的参数有时间t 、日期N 、地理经度λ、纬度ϕ和直杆高度H 。为了建立它们与影长l 之间的关系，需要找到一个中间量，而由于影长与直杆高度、太阳高度角之间存在
一个三角函数关系，而太阳高度角又是时间t 、日期N 、地理经度λ、纬度ϕ的函数，因此可以将两者联系起来，进而求得影子长度变化模型l =f (λ, ϕ, N , t , H ) 。
关于问题二，已知的附件1数据是影子顶点坐标数据、日期和时间，要求直杆高度和地理经度、纬度。模型中有三个变量是未知的，因此如何准确的找到一组解，使得其代入模型后能够很好的与真实影长拟合，成为本问题的关键。因为变量过多，可以利用最小二乘准则，采取网格搜索的方法筛选出最优解。为提高搜索效率，可以通过其他约束条件事先缩小变量的搜索范围。
关于问题三，附件2、3数据是影子顶点坐标数据、日期和时间，要求直杆高度和地理经度、纬度和日期。问题三与问题二的区别仅在于缺少一个日期变量，其余则完全相同，于是问题解决同样可以利用最小二乘准则，采取网格搜索的方法进行求解，事先依然可以通过约束条件缩小变量搜索范围。
关于问题四，附件4所给的数据由顶点坐标变成了视频，因此本问的难点在于如何将视频数据读成顶点坐标数据。影长数据可以通过从视频中截取图片，从图片读取的方式获得，但是由于图片中实物存在一定宽度，需要对图像进行预处理，而且由于拍摄角度问题，还需要对影长数据进行校正。得到影长数据后，这个问题就变成已知影长、日期时间、直杆高度，求地理经度、纬度的问题，与问题二类似，且已知直杆高度，求解更为简单。而最后一个问题只是增加一个变量日期，与问题三类似。
三、模型假设
1. 假设太阳光垂直入射地表。
2. 假设地球为均匀球体，入射表面为水平面，无地形起伏。
3. 由于我们所求影长时间在9:00-15:00之间，此时太阳高度角较大，因此可以忽略假设大气层、电离层、磁层对太阳光的影响。
四、定义与符号说明
五、模型的建立与求解
（一）问题一模型的建立与求解
1. 影子长度变化的数学模型
要建立影子长度变化的模型，必须先分析出该模型涉及到的参数变量，但其实问题中已经通过提问2015年10月22日北京时间9:00-15:00之间天安门广场（北纬39度54分26秒, 东经116度23分29秒）3米高的直杆的太阳影子长度的变化曲线给出了这些参数变量，即日期N 、时间t 、地理经纬度λ、ϕ、直杆高度H ，再加上影长l 就凑足了模型涉及到的所有变量。建立影子长度l 与地理经纬度λ、ϕ、直杆高度H 、地方时t 之间的关系模型需要借助太阳高度角θ这个变量。
（1）太阳高度角θ的求解
通过查阅相关资料得出太阳高度角θ的表达式为
sin θ=sin ϕsin δ+cos ϕcos δcosh (1)
θ是太阳高度角，δ是太阳赤纬，h 是时角，通过查阅相关资料得出太阳赤纬δ表示式分别为【1】：
δ=0. 006918-0. 399912⋅cos β+0. 01025⋅sin β-0. 006758⋅cos 2β+0. 000907⋅sin 2β式中β=2π⋅N /365 (2) 时角h 的表达式是【1】
h =⋅t +λ- (3)
为减小误差，时角的计算使用经度作为变量，而不是在某个时区内固定使用一个时角。
（2）影子长度l 与直杆高度H 和太阳高度角θ之间的关系
如图1，影子长度l 与直杆高度H 和太阳高度角θ
之间的数学关系为
l =H /tan θ (4)
图1 太阳高度角示意图
（3）模型的求取
由公式（1）-（4）可得出影子长度l 与地理经纬度λ、ϕ、地方时t 和直杆高度H 之间的模型函数：
l =H /{tan{a sin[sinϕsin δ+cos δcos ϕt +λ-]}},
δ=0. 006918-0. 399912⋅cos β+0. 01025⋅sin β-0. 006758⋅cos 2β+0. 000907⋅sin 2β
, β=2π⋅N /365 (5)
2. 影长l 关于各参数的变化规律
我们固定其中四个参数为常量，只留下一个参数作为变量，得到影子长度与该变量之间的变化曲线，分析其中的规律。
①影子长度l 与日期N 之间的关系
我们固定λ=116°E ，ϕ=40°N ，H =3m ，t =12h ，日期选择在5月份、8月份和11月份，分别绘制三个月的影子长度变化曲线图（如图2）。从图中可以看出，在固定地点λ=116 E ，ϕ=40 N 处5月份正午12时影子长度随日期的增加而减少，3m 长的直杆其影子变化范围在1~1.5m之间，而在同样地点，8月份和11月份正午12h 影子长度随日期的增加而增加，8月份影子变化范围在1.1-1.7m 之间，11月份影子变化范围在2.9-4.3m 之间。由于我们所选取的地点位于北半球、东半球，分析可得，5月份时，太阳相对地球的位置从赤道向北回归线移动，距离选点越来越近，高度角越来越大，因此正午时分的影长也越来越小，而8、11月份时，太阳相对地球的位置从赤道向南回归线移动，距离选点越来越远，高度角越来越小，正午时分的影长也越来越长。11月份的影长明显要大于5月份的影长，这是因为夏至日在6月21日左右，也就是太阳在6月21日左右直射赤道，而5月份距离夏至日较11月份距离夏至日明显要近，因此前者较后者距离太阳也就更近，高度角更大，影长更小。同理8月份与5月份的影长则要接近一些。除此之外，我们还可以看出，影长与天数大致呈线性相关关系。
图2 5月份（左）、8月份（中）、11月份（右）影子长度变化曲线图
②影子长度l 与经度λ之间的关系
我们固定地方时t 为2015年10月22日正午12时，杆长H =3m ，给定经度范围为100 ~130 E ，绘制出纬度在20 N ，40 N ，60 N 处，该经度范围内影子长度的变化曲
线，如图3。从图中可以看出，在我们选取的经度范围内，三个纬度都对应在λ=120E 左
右影长达到极小值，该经线以西，影长变大，以东，影长也变大，并且关于λ=120 E 左右对称，整体呈一条圆滑的下凹曲线。根据时角公式计算可得，正午12时整λ=120 E 时时角达到最小值0，此时高度角最大，因此影长最小，而且当时间t 确定时时角关于λ的函数是一个简单的带有绝对值的一次线性函数，因此必然是关于最小值这一点对称的。由此我们也可知，某一时刻影长关于经度曲线是一条对称曲线，且对称经度为此时刻时角为0的经线。
图3 经度100 ~130 E 范围不同纬度影子长度变化曲线图
③影子长度l 与纬度ϕ之间的关系
我们固定日期、时间为2015年10月22日正午12时，杆长H =3m ，给定纬度范围为40 S ~40 N ，绘制出当λ=60 E ，λ=120 E , λ=160 E ，该纬度范围内影子长度的变化曲线（如图4）。从图中可以看出在我们选取的纬度范围内，λ=120E 时，ϕ=12 S 左右影长达到极小值0，说明在这一时刻太阳直射λ=120E ，ϕ=8 S 地点处，往南往北影子长度均增加。而在其他两个经度线上影长也是存在一个极小值，但并不为0，从极小值向两侧影长数据是增加的。
图4 纬度40°S~40°N 范围影子长度变化曲线图
④影子长度l 与时间t 之间的关系
我们固定日期为2015年10月22日，杆长H=3m，选择三个不同地点分别是λ=120E ， ϕ=40 S 、λ=90 E ，ϕ=10 S 、λ=140 E ，ϕ=60 S ，绘制当地一段时间内影长变化曲线图。如图5。我们可以清晰的看出三个地点的曲线都是关于某一时刻对称的，而且存在一个极小值，并且三个地点的对称点也是极小点对应的北京分别为12:00,14:00,10:42。从时角公式的角度考虑，在每一条经度线上都存在一个时刻，这个时刻对应的时角值为0，即太阳高度角最大，影子最小。因此曲线必然存在一个极小值，这个极小值对应的时刻也就是时角为0时对应的时刻，并且时角是带有绝对值的简单一次线性函数，因此影长曲线也必然
是关于该极小值是对称的。
图5 2015年10月22日三个不同地点（从左到右依次为λ=120°E ，ϕ=40°S 、λ=90°E ，ϕ=10°
S 、λ=140°E ，ϕ=60°S ）影长随时间变化曲线图
3. 北京时间2015年10月22日9:00-15:00之间天安门广场（北纬39度54分26秒，东经116度23分29秒）3米高的直杆的太阳影子长度的变化曲线
题目中要求应用建立的模型画出2015年10月22日北京时间9:00-15:00之间天安门广场（3954'26''N ，11623'29''E ）3米高的直杆的太阳影子长度的变化曲线。已知的参数
变量为λ、ϕ、H 、N ，因此此问题即是固定这几个参数变量，求解影子长度l 与地方时t 之间的关系曲线，将已知参数代入模型，可以得到影长l 关于时间t 的数据结果，绘制成图6。
从图6中我们可以看出，该时间段内影长在3.3-6.6m 之间变化，在12时15分左右达到极小值，向两侧处影长增加，曲线整体呈一条下凹的曲线，并且仔细看不难发现，在9:30-15:00范围内该曲线关于t=12:15对称，类似于一条抛物线，分析可知，在之前的模型
分析中计算时角时我们代入了λ这一变量，并且通过计算得出λ=116E 时在t=12:16处时
角为0，也就是在此时刻高度角最大，因此影长也最小，而且当经度λ确定时时角关于t 的函数是一个简单的带有绝对值的一次线性函数，因此必然是关于最小值这一点对称的。这就不难理解该曲线在t=12:15时刻处达到最小，且关于该点对称了。
图6 北京时间9:00-15:00之间天安门广场3米高的直杆的太阳影子长度的变化曲线图
（二）问题二模型的建立与求解
题目中有三个未知量：直杆高度H 、地理经度λ、纬度ϕ，已知量有影子顶点坐标x 、y 、时间t 、日期N 。此题目可以将顶点坐标数据转化成影长数据，那么问题就等同于已知l 、t 、N 求H 、λ、ϕ。我们利用最小二乘准则，采取网格搜索的方法将这三个变量在规定的误差范围内筛选出一组最优解。由于变量较多，可以通过时角公式加二次拟合的方法缩小经度λ的搜索范围。
1. 根据时角公式估计经度λ值
由时角公式
h =⋅t +λ-
可以看出时角h 是时间t 和经度λ的函数，因此我们只要估计出时角为0时所对应的时间t 0，即可估计经度λ的大概值了。在之前的分析中我们得到，影长随时间变化曲线极小值对应的是高度角最大的值，而高度角最大的值对应的是时角为0的值，因此我们可以拟合影长随时间变化曲线，找到极小值对应的时间t 0，这个时间对应的就是该地时角为0的时间，将其代入时角公式，令时角为0，就可以找到当地经度λ了，λ的计算公式如下：
λ=300-15t 0 (6)
将影长随时间变化曲线绘制成图，再利用matlab 中的多项式拟合功能，通过反复尝试发现二次拟合曲线与真实影长曲线符合最好，因此我们选用二次拟合得到的结果，得到影长曲线关于时间的表达式：
l =0. 1489t 2-3. 7519t +24. 1275（t 单位为h ）
得到9:42-15:42之间时长6个小时的影长曲线，计算二次曲线极小值在t 0=12.6h处，如图7从图中以及实际数据中也可以得到曲线在12:36分左右达到最小值，也就是该地点的零时角时刻为12:36，再根据时角计算公式计算可以得到该地经度在东经111度左右。
图7 二次多项式拟合曲线图与真实影长曲线图
2. 网格搜索确定直杆高度H 、地理经度λ、地理纬度ϕ
利用问题一中建立的影子长度模型（5）, 基于最小二乘准则，其表达式为：
σ=min ∑(l 模型-l 实际) （7）
i =1212
采取网格搜索的方法，初步圈定直杆高度H 的范围为0-5m ，步长为1m ，纬度ϕ的范围为
90 S ~90 N ，步长为1°，λ的搜索范围为101E ~121E ，步长为1°。这样我们就得
到由若干个元素组成的一维数组，每个元素由（H ，ϕ，λ）组成。网格搜索流程图如图
8。利用数据中的时间值，以及数组中的任一元素，代入到长度模型中，计算得到21个影长值，再利用顶点坐标数据x 、y 计算得到真实影长值，将其与真实影长值相减计算得到一个均方误差值，将数组中每一个元素都照此操作，选择均方误差较小的元素范围作为下一次搜索的初始范围，以此类推，当数组元素达到我们要求的精度范围时，即可得到准确的直杆高度H 和地理纬度ϕ值，三次搜索的范围分别为H=0:1:5，ϕ=90 S :1:90 N ；H =1:0.1:3，ϕ=22 S :1:45 N ；H =1.5:0.05:2.5，ϕ=19 N :1/3600:21 N ，最终得到H =1.80m，λ=111°01′25″，ϕ=19°57′30″N ，将此元素代入模型得到的均方误差为6.00109×10-6，如图9，为H =1.80m，λ=111°01′25″，ϕ=19°57′30″N 时所得到的影长与真实影长数据之间的拟合结果，拟合情况良好，符合精度要求。通过经纬度查询发现该地点大致在海口市附近。
.
图8 算法流程图
图9 模型计算影长与真实影长拟合效果图
（三）问题三的求解
题目中有四个未知量：直杆高度H 、地理经度λ、纬度ϕ、日期N 。已知量有影子顶点坐标x 、y 、时间t 。与问题二不同的是，本问多一个变量日期N ，而模型则完全相同。问题二的解决方式是利用最小二乘准则加网格搜索的方法，由于网格搜索本身的特点使得其能够解决多变量问题，因此在这我们同样可以采取网格搜索的方法，N 也作为变量代到模型中去求最优解。而经度λ同样可以采取时角公式加二次拟合先确定下大致范围来。
1. 附件2的求解
（1）根据根据时角公式估计经度λ值
将影长随时间变化曲线绘制成图，再利用matlab 中的多项式拟合功能，得到影长曲线关于时间的表达式，根据表达式计算和读图求得极小值对应的时刻，也就是零时角时刻，通过时角计算公式
h =⋅t +λ-
即可先大致确定出当地经度。
二次拟合表达式是
l =0. 0981t 2-2. 9850t +23. 3192（t 单位为h ）
得到9:00-15:00之间时长6个小时的影长曲线如图10，零时角时刻在t =15.2075h 处，计算可以得到该地经度大致在72°E 左右。
图10 二次多项式拟合曲线图与拟合影长曲线图
（2）网格搜索确定地理纬度ϕ和直杆高度H 和日期N
目前已知的参数有时刻t ，影子顶点坐标xy 以及λ的大致值。利用问题一建立的影子长度模型（5），给定地理经度λ、纬度ϕ、直杆高度H 和日期N 的初始搜索范围，其中λ的搜索范围为之前求出的大概值左右各10°，比如，附件二中求出的经度大概在72度左右，那么经度的初始搜索范围则为62E ~82E ，步长为1°，ϕ的搜索范围为90S ~90N ，搜索步长为1°，H 的搜索范围为0-5m ，步长为1m ，日期N 的搜索范围为1-365，步长为1，将这几个变量和已知参数代入模型中，规定一个均方误差，在变量取值范围内遍历搜索，找到满足该均方误差的几组解，每个解都由地理经度λ、纬度ϕ、直杆高度H 和日期N 组成，从实际程序运行结果来看，第一步的搜索中可以将H 确定到米量级，再将均方误差缩小，并缩小H 的搜索范围，步长变为0.1m ，其他不变，就能将H 确定到分米量级，其他参数依此类推，最终我们可以得到一组最优解。值得一提的是，最终得到的解对应两个日期N ，这是因为太阳从赤道到北回归线再回到赤道向南回归线再回到赤道这一年的周期性变化中，必然有两天的太阳位置完全相同，因此在每个地点必然存在两天，这两天的影子变化完全相同，并且这两个日期是关于夏至日或者冬至日对称的。
所求解的均方误差σ=2. 5462⨯10，最终得到解为：H =1.9m，λ=7801'23''E ， -7
ϕ=40 01'15''N ，N =155、210。将日期N 换算成以年月日为单位，附件二得到的最终解为：H =1.9m，λ=78°01′23″S ，ϕ=40°01′15″，日期为6月4日，7月28日。地点在乌鲁木齐附近。图11为该解代入模型后计算得到的影长值与真实影长值拟合效果图。
图11 模型计算影长与真实影长拟合效果图
2. 附件3的求解
（1）根据时角公式估计经度λ值
二次拟合表达式是
l =0. 2964t 2-7. 5507t +51. 5639（t 单位为h ）
得到9:00-15:00之间时长6个小时的影长曲线如图12，零时角时刻在t=12.7373h处，计算可以得到该地经度大概为东经109度左右。
图12 二次多项式拟合曲线图与拟合影长曲线图
（2）网格搜索确定地理纬度ϕ和直杆高度H 和日期N
均方误差σ=2. 5115⨯10，最终得到的解为：ϕ=27°λ=109°57′02″E ，H =3.0m，
10′42″N ，N =350、13。将日期N 换算成以年月日为单位，附件3得到的最终解为：
-7
H =3.0m，λ=109°57′02″E ，ϕ=27°10′42″N ，N =1月14日，12月16日。地点在长沙市附近。图13为该解代入模型后计算得到的影长值与真实影长值拟合效果图。
图13 模型计算影长与真实影长拟合效果图
（四）问题四的求解
题目的已知量有直杆高度H 、一段视频数据、日期N 、时间t ，未知量有影子长度l 、地理经度λ、纬度ϕ。本问的难点在于将影子长度从视频数据中获取出来，一旦影子长度已知，那么问题就和问题二完全相同了。视频是由一系列照片组成的，我们需要首先将照片截取下来，然后从照片中读出影子变化的数据。该数据可以是坐标，可以直接是影长，但不管怎样我们都需将其转化成影长l ，这样才能应用到我们建立的模型中去，而求取λ、ϕ的方式则与问题二完全相同，采取网格搜索的方式。
1. 图像预处理
根据之前处理的数据，我们在8:54:06-9:34:06时间段内，选择时间间隔为2分钟，使用视频播放器等时间间隔截取出21张照片。由于照片中的实物无论直杆还是影子都存在一定宽度，如果在原始照片上直接读取数据，必然会带来了一定的误差。所以我们将照片在CorelDRAW 软件中打开，利用CorelDRAW 软件中位图选项下的模式转换功能将图片转换成黑白线条图，并选择一定的阈值使得影子和直杆清晰的显现出来，如图所示，上面的图是将阈值调低得到的图，放大后就可以清晰的找到杆顶的位置，下面的图是将阈值调高得到的图，可以看到影子变成一条直线，这样我们就可以找到影子顶点的位置，并且影子还是由两条不同的直线叠加在一起，这说明这根2m 的直杆是由两段直杆绑在一起组成的，这也与实际模型相符。
图14 阈值较低时图像预处理后得到的图
图15 阈值较高时图像预处理后得到的图
2. 照片中顶点坐标数据的读取
将预处理后的图像在CorelDRAW 软件中打开，由于直接读取长度较为麻烦，因此我们选择读取坐标数据。需要读取的坐标数据有直杆底部的x 、y 坐标、直杆顶点的x 、y 坐标、影子顶点的x 、y 坐标。读取过程中将照片完全叠放在一起，这样就能保证所有照片的直杆底部都位于一个点上，而我们所需要的就是这一个点的坐标数据作为所有照片中直杆底部的坐标。将鼠标放在需要读取的点上即可直接在软件下方状态栏读取数据。得到的数据如表1.
3. 影长的求取与校正
因为直杆高度是已知的，所以可以通过照片直杆与影长的比值将实际影长求取出来。又因为题目中直杆是由两段直杆绑在一起的组成，所以读取的直杆底部横坐标与直杆顶部横坐标是存在一定差距的，又因为假设直杆与地面垂直，而且给出的数据也是高度而不是长度，故计算直杆高度时采用直杆顶部纵坐标减去直杆底部纵坐标的方法，那么影子长度l 与坐标数据之间的关系是：
l =(H ⋅y 1-y 0) /((x 1-x 0) 2+(y 1-y 0) 2) (m)
式中x 1、y 1是直杆顶部横纵坐标，x 0、y 0是直杆底部横纵坐标，H =2m，x 、y 分别为影子顶点横纵坐标。
影长的校正：根据计算公式得到的影长在表1中给出，进行校正。如图16所示，实际所见的实际影长与图片影长有如下关系：假设相片影长与杆在一个平面内，杆高的两米所在的比例尺与图片影长比例尺可以认为相同。基于此可以推导出实际影长与图片影长的关系为实际影长=图片影长/cos(θ) ，θ为实际影长与图片影长的夹角。对应不同的θ角度实际影长与图片影长的关系就不同。据观测，在视频的最后θ角接近0度。整个视频期间θ角变化不大，cos(0︒) =1 因此实际影长校正后仍等于图片影长。
求取的数据以及计算的结果如下：
表1 影子顶点坐标求取数据表
图16 图片影长与真实影长关系示意图
4. 根据时角公式估计经度λ值
利用曲线二次拟合和时角公式计算得到λ的估值。
二次拟合表达式是
l =0. 1264t 2-3. 1632t +20. 5621（t 单位为h ）
得到8:00-15:00之间时长7个小时的影长曲线如图17，零时角时刻在t =12.5127h处，计算可以得到该地经度大致在东经111度左右。
图17 二次拟合影长曲线与真实曲线拟合效果图
5. 网格搜索确定地理纬度ϕ
将影长数据、直杆高度、日期时间都代入模型中，采取网格搜索，求得最优解是：λ=109°00′39″E ，ϕ=41°55′39″N ，在包头市附近，均方误差是σ=4. 4941⨯10-5。图18是将所求解代入模型求取的影长值与真实影长值拟合的效果图。
图18 模型影长与真实影长拟合效果图
关于“如果拍摄日期未知，能否根据视频确定出拍摄地点与日期”这一问题，我们考虑过后，认为此问题是已知影长数据、时间、直杆高度，求取地理经度、纬度和日期的问题，是可以求解的，具体求解步骤如下：
①图像预处理；
②照片中坐标数据读取；
③影长数据的求取；
④根据时角公式估计经度λ值；
⑤网格搜索确定地理纬度ϕ和日期N 。
前四个步骤与问题四的求解是完全相同的，只是在第五步网格搜索的过程中增加一个变量日期。
六、模型的评价
我们利用直杆高度、影长和高度角之间的三角函数关系，并通过高度角这一辅助量将影长与地理经纬度、日期、时间联系起来，最终得到影子长度变化的模型。模型的优点在于简单易懂，比较容易用程序实现，运行速度也较快，但是模型也遇到一些问题，我们做了相关考虑。
（1）影子坐标数据转换成影长数据
题目中给出的都是影子坐标数据，但我们都是将其转化成影长数据再进行处理，可能会丢失掉一些信息，因此我们在求解过程中又建立了一个方位角模型，其大致原理如图19所示，其中影子顶点坐标x 、y 与方位角之间的关系式是
x
/y =tan(A -α)
图19 角度示意图
而太阳方位角A 又是地理经度、纬度、时间、日期的函数，这样我们就可以建立出影子横纵坐标x 、y 与地理经度λ、纬度ϕ、时间t 、日期N 之间的模型。求解过程中利用最小二乘准则，采取网格搜索的方法求出最优解，如图20是我们求出的解代入模型后求得的x /y 比值与真实x /y 拟合效果图，效果并不是很好，并且均方误差也较大，实际求得的纬度也与影长模型求得的纬度存在一定差距。我们猜想这是因为方位角公式本身的不准确性造成的，因为目前得到的方位角公式基本都是经验所得，存在一定的误差，因此我们最终选取了影长模型求得的解。
图20 模型x /y 比值与真实x /y 比值拟合效果图
（2）图片影长的校正
在从图片中读取影子数据时，由于拍摄的平面与影子和杆所在的平面存在一个角度，而如果该角度过大，将会影响真实影长值的获取，因此我们需进行校正，但是分析可得在视频记录时间段内，该夹角并不大，十分接近0，因此这一校正实际上是可以忽略的。
七、参考文献
[1] 王国安， 米鸿涛， 邓天宏等， 太阳高度角和日出日落时刻太阳方位角一年变化范围的计算，气象与环境科学，30(增刊) ：161-164，2007
[2]于贺军，气象用太阳赤纬和时差计算方法研究，气象水文海洋仪器，3：50-53，2006
[3] 贺晓雷，于贺军，李建英等，太阳方位角的公式求解及其应用，太阳能学报，29(1)：69-73，2008
八、附录
附录一-------------------------------------------------------------------------影子长度求取的matlab 程序 附录二----------------------------------------------------------------网格筛选（粗筛）的matlab 程序 附录三----------------------------------------------------------------网格筛选（粗筛）的matlab 程序 附录四----------------------------------------------------------------------经度拟合缩小范围matlab 程序 附录五-------------------------------------------------------------------------太阳方位角计算matlab 程序 附录六-------------------------------------------------------------------------------网格筛选结果数据表格
附录一-------------------------------------------------------------------------影子长度求取的matlab 程序 MATLAB （R2013a ）编程代码
clear all ;clc
l=3; %竿子高
Ts=datenum( '2015-01-01' );
Te=datenum( '2015-12-31' );
Ts1=datenum( '2015-01-02' );
Tt=datenum( '2015-04-18' );
Tn=Ts1-Ts;
T1=Ts:Tn:Te;
for i=1:length(T1)
if T1(i)==Tt;
N=i-1; %求取积日
end
end
tstart=9; %开始时间
tend=15; %截止时间 （小时制）
tn=1/20; %每小时20个点
t=tstart:tn:tend;
st=2*pi*N/365;
c=0.006918-0.399912*cos(st)+0.010257*sin(st)-0.006758*cos(2*st)+0.000907*sin(2*st);%太阳赤经计算
wei=39+54/60+26/60/60; %弧度转换角度
jing=116+23/60+29/60/60;
s=abs(t*15+jing-300); %时角
gao=asin(sin(wei/180*pi)*sin(c)+cos(wei/180*pi)*cos(c)*cos(s/180*pi)); %太阳高度角
x=l./tan(gao);
plot(x);%画出模拟影长
hold on ;
plot (sj4,'r' );
hold on ; %画出实际影长
附录二----------------------------------------------------------------网格筛选（粗筛）的matlab 程序
sx1=zeros(50,180,20,365);
jie=zeros(1000,4);
n=1;
for i=1:1:20
for j=1:1:90
for k=1:365
for p=1:20
l=i/10; %杆长
wei=j; %纬度
jing=60+p; %经度
tstart=12+41/60; %开始时间
tend=13+41/60; %截止时间 （小时制）
tn=1/20; %每小时20个点
N=k-1;
t=tstart:tn:tend;
st=2*pi*N/365;
c=0.006918-0.399912*cos(st)+0.010257*sin(st)-0.006758*cos(2*st)+0.000907*sin(2*st);
s=abs(t*15+jing-300);
gao=asin(sin(wei/180*pi)*sin(c)+cos(wei/180*pi)*cos(c)*cos(s/180*pi));
x=l./abs(tan(gao)); %影长
for m=1:21
sx1(i,j,p,k)=sx1(i,j,p,k)+(x(m)-yc(m))^2; %误差 end
if sx1(i,j,p,k)
jie(n,1)=jie(n,1)+i;
jie(n,2)=jie(n,2)+j;
jie(n,3)=jie(n,3)+p;
jie(n,4)=jie(n,4)+k;
n=n+1;
end
end
end
end
end
sx1=sx1/21;
jie1=jie(1:100,:); %选取误差小的对应数据值
附录三----------------------------------------------------------------网格筛选（粗筛）的matlab 程序
sx2=zeros(7200,7200);
for i=1:1:7200
for j=1:1:7200
wei=j/3600+39; %纬度
jing=i/3600+77; %经度
t=tt;
st=2*pi*N/365;
c=0.006918-0.399912*cos(st)+0.010257*sin(st)-0.006758*cos(2*st)+0.000907*sin(2*st);
s=abs(t*15+jing-300);
gao=asin(sin(wei/180*pi)*sin(c)+cos(wei/180*pi)*cos(c)*cos(s/180*pi));
x=l./abs(tan(gao)); %影长
for m=1:21
sx2(i,j)=sx2(i,j)+(x(m)-yc(m))^2;
end
end
end
sx2=sx2/21;
xiao=min(min(sx2)); %最小误差的获得
[jjing,jwei]=find(sx2==xiao); %最精确位置的获得
jwei=39+jwei/3600;
jjing=77+jjing/3600 ;
附录四----------------------------------------------------------------------经度拟合缩小范围matlab 程序
tt=12.6833:0.05:13.6833;
a=polyfit(tt,yc,2);
ty=9:0.05:17;
ni=a(1)*ty.^2+a(2)*ty+a(3); %应用二次多项式拟合效果最好
%ni=a(1)*ty.^3+a(2)*ty.^2+a(3)*ty+a(4);
%ni=a(1)*ty.^4+a(2)*ty.^3+a(3)*ty.^2+a(4)*ty+a(5);
plot(ty,ni,'r' )
hold on ;
plot(tt,yc)
hold on ;
jing=a(2)/(a(1)*2)*15+300;
附录五-------------------------------------------------------------------------太阳方位角计算matlab 程序
ww=zeros(5,180);
xg=zeros(450,7200);
for jia=1:1:450
for we=1:1:7200
tstart=14+42/60; %开始时间
tend=15+42/60; %截止时间 （小时制）
tn=1/20; %每小时20个
l=1.8;
jiao=jia/10; %夹角
wei=16+we/3600; %纬度
N=108-1;
jing=111; %经度
t=tstart:tn:tend;
st=2*pi*N/365;
c=0.006918-0.399912*cos(st)+0.010257*sin(st)-0.006758*cos(2*st)+0.000907*sin(2*st);
s=abs(t*15+jing-300);
gao=asin(sin(wei/180*pi)*sin(c)+cos(wei/180*pi)*cos(c)*cos(s/180*pi)); %高度角
x=l./tan(gao); %影长
A=acos((sin(gao)*sin(wei/180*pi)-sin(c))./(cos(gao)*cos(wei/180*pi))); %太阳方位角的计算
fang=tan(A-jiao/180*pi); %x/y的比例
for k=1:1:21
xg(jiao*10,we)=xg(jiao*10,we)+(tan(A(k)-jiao/180*pi)-sj3(k))^2; end
% xg(l,wei+90)=(x(21)-x(1))/1-(sj4(21)-sj4(1))/1;
% for k=1:21
% ww(l,wei+90)=ww(l,wei+90)+(x(k)-sj4(k))^2;
%end
end
end
ww=ww.^0.5;
ww=ww/21;
bb=min(min(ww)); %得到最小值
%ba=min(min(xg));
[hedd,sdd]=find(xg==ba); %找到最小值
附录六-------------------------------------------------------------------------------网格筛选结果数据表格
第二问网格筛选结果数据表
第三问附件2网格筛选结果数据表
第四问网格筛选结果数据表