太阳影子定位,2015数学建模国赛A题资料

对太阳影子定位算法探究

摘要

本文是对2013年全国大学生数学建模竞赛A题的解答.随着人们对数据挖掘的深入，如何确定视频的拍摄地点和拍摄日期已经成为视频数据分析的重要方面，太阳影子定位技术就是通过分析视频中物体的太阳影子变化，确定视频拍摄的地点和日期的一种方法，进而可以促进视频分析定位技术发展。

对于问题一，我们根据地球自转公转的自然规律，建立影子长度变化的数学模型，并且分析影子长度关于各个参数的变化规律。基于对问题的分析以及理论的学习研究，画出模拟概念图，然后计算相关量（如太阳高度角、赤纬角等）的表达式，并按照相关地理知识建立起模型，得到杆子影长与时间函数表达式，再将题目所给的数据代入求解方程，并用MATLAB作出曲线图，最后检验模型的准确性。

对于问题二，我们以问题一所求出的表达式和资料作为基础，继而利用球面天文学求算太阳视坐标的简化算法建立一模型直接求解出经度，纬度的估算值。再代进数据并用利用多项式拟合出更长的时间序列曲线，用函数的特征值（最低点）加上时角，时区计算相关知识，再推算出经度值。最后利用第一问模型，经度，加上曲线获得的几组影长数据联立求解出大致纬度，最后估算杆子所在的地区。

对于问题三，结合问题一问题二所建立的模型，将附件2，附件3的数据先画出散点图并以多项式拟合出两条相对完整的曲线，通过其曲线函数求得影长的最小值以及最小值所对应的时间求得经度，纬度，将经度和纬度代入赤纬角公式以及影长公式可求得相应的具体日期。

对于问题四，首先将视频数据利用MATLAB，并且编程处理视频得到每分钟一帧的图片，再把相关图片转化为灰度图矩阵，最后用语句转化为二值图（0为黑，1为白）。下一步把二值图集分析并且分析出杆子影长的变化规律，求出视频拍摄点经度，利用模型一求出纬度，即是位置。

关键字：影长位置 MATLAB编程多项式拟合最小二乘法二值图

一．问题重述

如何确定视频的拍摄地点和拍摄日期是视频数据分析的重要方面，太阳影子定位技术就是通过分析视频中物体的太阳影子变化，确定视频拍摄的地点和日期的一种方法。

1.建立影子长度变化的数学模型，分析影子长度关于各个参数的变化规律，并应用你们建立的模型画出2015年10月22日,北京时间9:00-15:00之间天安门广场,3米高的直杆的太阳影子长度的变化曲线。

2.根据某固定直杆在水平地面上的太阳影子顶点坐标数据，建立数学模型确定直杆所处的地点。将你们的模型应用于附件1的影子顶点坐标数据，给出若干个可能的地点。

3. 根据某固定直杆在水平地面上的太阳影子顶点坐标数据，建立数学模型确定直杆所处的地点和日期。将你们的模型分别应用于附件2和附件3的影子顶点坐标数据，给出若干个可能的地点与日期。

4．附件4为一根直杆在太阳下的影子变化的视频，并且已通过某种方式估计出直杆的高度为2米。请建立确定视频拍摄地点的数学模型，并应用你们的模型给出若干个可能的拍摄地点。

如果拍摄日期未知，你能否根据视频确定出拍摄地点与日期？

注：题目及数据附件都可到全国大学生数学建模竞赛官方网站（https://www.360docs.net/doc/0819212225.html,）下载.

二．模型的假设

1)假设地球在自转时角速度均匀；

2)假设地球在公转时的椭圆为规则椭圆；

3)假设题目涉及的太阳时计算为真太阳时；

4)假设在一天内的太阳高度角不变；

5)忽略大气折射对阳光传播的影响。

注：1真太阳日分为24真太阳时。这个时间系统称为真太阳时。真太阳时也称为视太阳时，简称视时。真太阳时即真太阳视圆面中心的时角加12小时。即：真太阳时=平太阳时+真平太阳时差

三．符号说明

表1

四．模型的建立与求解

4.1问题一

4.1.1问题1分析

题目要求描建立影子长度变化的数学模型，分析影子长度关于各个参数的变化规律。首先通过查阅相关地理、物理文献，找寻建模中可能遇到的概念、方法，再挖掘建立模型相关的方法，为之后的建模奠定坚实的基础。模型建立第一步,理清自变量（四个），因变量，选取参考系（参考面还有参考点），建立空间直角坐标系。再做出太阳随时间移动杆影的直观视图（日照杆影图），最后，利用太阳高度角、杆高等相关量与杆影的几何关系，构造对应函数，使之符合题意，同时根据相应的标准进行统计、分析和构建相关的数学模型。

4.1.2模型1的建立

通过对题目的探讨，一根垂直地面的杆子，研究它的影子长度随着日照光线的改变而改变，首先我们确定相应的参考系，以杆子所在地地平面α为参考面，以杆子的最低端为坐标原点，再得到以下的自变量：

表2

因变量：杆子在参考面上投影长度（可以从杆子顶点在参考面的投影求得）

4.1.3步骤一：求出建立模型必须的物理，地理相关量值量

在此基础上我们分析地球运动的特性规律与研究杆影与时间、日期、经度、

纬度、杆长的相关关系，首先做出日照杆影分析图1-1，如下:

图1-1

设杆子顶点在地平面上投影点相对于原点坐标为Q(x0,y0),另外杆子高度为H ，影子长度为L 。太阳高度角可求：

)arctan(1

L H

=θ; 从图1-1中分析，太阳与杆子相夹的角近似为太阳高度角，

太阳赤纬角：定义是又称太阳赤纬，是地球赤道平面与太阳和地球中心的连线之间的夹角。

由图1-1几何关系，查阅资料以及推理出计算公式如下： 因为题中涉及到的是真太阳日所以太阳赤纬角的计算公式为

)365k)

(284×36023.5sin(ε ；+= ①

k 为日期序号，如若3月22号,k 就为81。 时角计算公式：

12;-15n t , t β== ②

太阳高度角θ计算公式：

；βεεcos cosWcos sinWsin sin θ+= ③

根据三角函数知识以及图中几何关系得，太阳方位角计算公式：

;c o s

s i n c o s s i n θβ

εφ= ④

4.1.4步骤二：利用上面所得的数据求出Q 点的轨迹方程（即为影子轨迹）

图1-2

对图1-2进行分析，发现在某个地点某个时刻太阳照射杆子杆影的变化曲线图的大致规律：一条类似于抛物线的曲线对杆影的因变量和自变量进行量化分析得影子长度,推并推算相关公式得：

；φtan 0y0?=x ⑤ ；φ

cos x 0

L =

⑥ 然后将上述式子联立求解得Q 点的x 坐标的表达式：

；

φ

θ2

tan 1cot 0x +?±

=H ⑦ 代入数据得

⑧

]]]

))))

12(15cos(cos W cos sin W (sin cos(arcsin ))

12(15sin(cos n[[tan[arcsi 1))))

12(15cos(cos W cos sin nW (arcsin(cot H x0；-???+?-??+-???+??±

=n n n si εεεεε 由;00L 22y x +=

再以时间t 为自变量，可求得影子长度的变化曲线：

⑨]]]

))))

12(15cos(cos W cos sin W (sin cos(arcsin ))

12(15sin(cos n[[tan[arcsi 1))))

12(15cos(cos W cos sin nW (arcsin(cot H tan 1L 222

2

；）（-???+?-??+-???+??+=n n n si εεεεεφ

如果数据充足，可以利用MATLAB 做出Q 点的轨迹，以及影子长度随着自变量的变化而变化曲线图，其中纬度被包含在公式中，经度被时角所替代，地方时间n ，日期也嵌套在赤纬角里面，这样四个自变量就可以表示杆影的长度变化曲线。

4.1.5步骤三：计算取值，作曲线

计算画出2015年10月22日北京时间9:00-15:00之间天安门

广

场（北纬

39度54分26秒,东经116度23分29秒）3米高的直杆的太阳影子长度的变化曲线。

由步骤二可知道 φ

β2tan 1cot 0x +?±

=H , φtan 0y0?=x , 2200L y x +=.

是主要公式，联合步骤二的有关公式可以得出直杆的太阳影子长度随时间的变化曲线为

θφ

θ

φHcos tan 1cot tan 100t L 2222

2

2

=+?+=+=H y x ）（）（. ⑩

其中cos θ是时间的函数，因为H 、W 、J 、n 、日期已知，代入数据并用MATLAB

求得杆子顶点的影子长度的变化曲线如曲线图所示：

图1-3 检验准确性

分析图1-3，可见曲线直观上是一条类似于抛物线的曲线，表明3米高的杆子在

北京天安门广场9:00到15:00的杆影长度曲线如图1-3。

图1-4来自于·影端轨迹周年变化的实践与分析_以北温带地区为例_吴济廉

再根据相关文献中的依据，其中包含实验验证此现象的记载，表明结果的准确性

程度较高。

4.2问题二

4.2.1问题二的分析

题目要求根据某固定直杆在水平地面上的太阳影子顶点坐标数据，建立模型

确定该杆所在的地理方位。附件一中有的条件是影长，时间，日期，不明确竿子

高度H ，求所测地点的经纬度。

4.2.2模型的建立与检验：

设竿子长度为H （未知），以问题一所建立模型为基础，利用球面天文学求算太阳视坐标的简化算法对问题进行探究，建立此问题的数学模型，并且带进附件一中的数据，再利用数据拟合出一天之中的竿影长度随时间的变化图，得出最短阴影长度值，进而得到该点经度值。估算出测验地的地理未知，并检验。 首先我们要计算某一时间的天顶距S （是指在天体方位圈上，天体与天顶之间的交角称为天顶距---来自百度百科），可以从天顶的定位三角形,利用三角几何关系求

图2-1百度百科上的天顶角示意图 天顶距S 的代数关系式子：

;cos cos cos sin sinW cosS βεεW += ? 时角的精确计算方式：

;150a -J ?+++=r r f （忽略）β ?

a0为视赤经，f 为真太阳时（忽略），r 为地方时，ε为视赤纬；

)2623cos arctan(tan 0a ?=d , );2623sin sin(sin a ?=d rc ε

太阳的视天顶距为

;sin 9tan 0306-S S S S ''+'''''='

利用最小二乘法进行拟合求观测地点的经纬度，杆子的高度与影长之间的关系:

.tan

)

cos cos cos sin sinW (arccos tan H βεεW S L +='=

?

当影长一定时，假设杆子高度H 与经度纬度的函数，即是当影长一定时, 可把由影长而推算的物高H 看作为纬度和经度的函数

();W,J f H =

，??????

?'

-=??-+=??S H S W t H 3'3s i n /s i n c o s J c o s L J

s i n /)s i n c o s c o s W c o s (s i n L W

βεεε ? 结合第一问所得的L-t 图，解微分方程求解经度大约106°，再用多项式拟合求经度与之比较。

代入已知数据之后，用多项式拟合法拟合出出杆影随时间的变化而变化图像，如下:

最低点是： ），（）（a

b a

c a b ac a 44598724.1244,2b -2

2-=-

1275

.247519.31489.0y 2+-=x x

图4-1阴影部分为实际区，蓝色部分为预测区

拟合度的检验，利用MATLAB 进行regress 拟合度命令的多次分析，拟合度b 的数值的平均值为0.9998723，证明在区域内，曲线的拟合效果优良。

接下来就是分析计算经纬度值，由拟合出来的曲线可以计算杆子影子最短时的时间，拟合方程所代表式子为基础条件，求经度：

测量最短杆影，地球自转一周需要24小时,也就是1440分钟,地球自转360度,因此,可以算出每4分钟,地球自转一度。用上面的时间和正午12点之差,就可以求出你的经度与当地时区经度的差。最短竿影下的时间约为12.60时即是12:36

左右，查阅资料获得东八区的经度为120°,可知杆子所在地的经度为

4

15

1260.12120J 1?-±

=）（. 解得经度为：120°±9°E 即是111°E 或者129°E.

取最短影长为参考值，此时时间为12：36，日期为2015年4月8日，由第一问提取到的公式联立有：

θφ

θ

φHcos tan 1cot tan 100t L 2

222

2

2

=+?+=+=H y x ）（）（ . ? ；βεεc o s c o s W c o s s i n W s i n s i n θ+= ?

综合解得：纬度的可能取值是40.55°N 和18.25°N 。

在问题一建立的模型的基础上，通过对附件1所给的数据进行多项式拟合出时间序列的二次函数曲线，并通过验证拟合度的方式保证该函数曲线的可靠性。接着通过分析二次函数最小值点以求得最小影长，以及最小影长所对应的时间。将对应的时间代入求得经度。结合最小影长和经度可求得纬度。

4.2.3 问题二的求解

由所求的函数曲线可得出其最小影长所对应的时间为12.598，结合地球自转360度，即每4分钟地球自转一度，用所求时间与正午12点之差即可求得该地的经度与东八区的经度差，再由东八区120度可求得两个位置的经度分别为111.03°E 和128.97°E ，再通过模型一所建立的模型以及公式可求得纬度为40.55°N 和18.25°N ，由GPS 定位可得该地理位置在呼和浩特和海南三亚附近地区。

4.3 问题三

4.3.1问题3的分析

题目3与题目2较为相似，尽管题目3多出了另一个附件的数据，但是这两个附件的数据都没有给出确切测定的日期。然而纬度，经度的值与最小影长有关，即因通过求得最小影长去求得相应的经度和纬度。又因为影响赤纬角公式以及影长公式中有日期，即将求得的经度，纬度代入公式中进而求得日期。

4.3.2模型3建立

在前两问的基础上，我们仍然采用多项式拟合的方式对附件2，附件3的函数进行处理，同时通过求出拟合度的方式验证其求得函数的可靠性。为求得该经度，纬度，仍然是在第二问的基础上求得影长的最小值以及所对应的时间。由赤纬角以及影长公式可知其与具体日期有关，故将经度，纬度代入即可求得具体日期。

4.3.3模型3的求解

首先对附件2和附件3的数据利用MATLAB进行多项式拟合，

Y=0.0981x^2-2.9883x+23.3691

图4-1阴影部分为实际区，蓝色部分为预测区

Y=0.2964x^2-7.2839+48.2661

图4-2阴影部分为实际区，蓝色部分为预测区

然后对两个曲线拟合模型进行regress拟合度分析，都得到的拟合度匹配程度较高。

由附件2，附件3的函数曲线可求得最小的影长为15.23和12.29，并代入求得纬度为40.40°N，经度为80.24°E，地点在新疆阿瓦提县附近，日期为6月18日附近。

4.4 问题四

4.4.1问题4的分析

对问题四的分析中，题目中已知杆高为2，另外题目提供一个40分钟视频附件，我们能从视频得到影子长度随时间的变化规律，这样已知数据包含时间，日期，杆子长度，以及杆影随时间的变化曲线，根据这些已知，利用模型一公式，可以求出视频拍摄地的位置。重要的是如何分析提取视频里面的数据内容，因为视频中的是立体图，所得到的矩阵数列要经过修正才能代表杆子的长度，才能进行位置的估计。

4.4.2模型4建立

首先我们把视频导入MATLAB中,用附件四的视频导入MATLAB中，然后取每N帧大概每隔一分取一张图片，存入数据库矩阵A，再利用MATLAB编程将所得的所有彩色图片简化成二维矩阵的灰度图片以png格式保存（如图4-1左），下一

步将灰度图片化成二值图片（如4-2右），得到简化的矩阵组合，最后分析矩阵得到杆影与时间的大致变化曲线。

1111

n n nn a a A a a ??

??=?

?????.

分析过程如图

4-1

图

4-1

图4-2

二值矩阵的建立

为防止灰度值流失，Matlab 计算的灰度值会被限制在0-255之间。在模型计算中，为方便确定模型中杆长及杆影灰度矩阵，需将灰度值矩阵A 中杆长及杆影部分转化成二值矩阵，具体操作：

对于杆长的灰度矩阵ij a ：

；255=a 时，60<=a ；当0=a 时，60>a 当ij ij

对于随时间变化的杆影灰度矩阵ij b ：

；255=B 时，195>b 且160

灰度矩阵的处理

所以杆影长度与时间的函数曲线如图4-3所示 利用MATLAB 对灰度矩阵进行运算，求解出杆长、杆影二值化矩阵边缘，从而确定合适的。 边界定位模型

为了得到附件4中随时间变化的杆影长度，采用边界定位法。它是基于二值化矩阵的两极化特性：二值化后的杆影灰度图像的右上边界的255灰度值与左下边界的255灰度值即为杆影的边界，即杆影的长度为两灰度值之间的距离。 边界定位模型的具体步骤如下： 步骤一：右上边界模型的构建

将第i 张(1,2,

,40)i =图片的右上边界进行右上边界定位，即求第i 张

(1,2,,40)i =图片矩阵从第j 行(1,2,

,51)j =由下向上，第k 列(1,2,

,998)

k =由左向右开始，检测灰度值首先为255的元素，并提取其位置数据——(),i i x y 。 步骤二：左下边界模型的构建

将第i 张(1,2,

,40)i =图片的左下边界进行左下边界定位，即求第i 张

(1,2,,40)i =图片矩阵从第j 行(1,2,

,51)j =由上向下，第k 列(1,2,

,998)

k =由右向左开始，检测灰度值首先为255的元素，并提取其位置数据——(),i i x y 。

步骤三：杆影长度的确定

将提取出的第i 张(1,2,,40)i =图片的右上边界位置数据(),i i x y 与左下边界

位置数据(),i i x y 进行计算：

22()()(1,2,

,40)

i i i i i d x x y y i =-+-= (1)

4.4.3 模型4的求解

以方向北为基准方向，结合模型一，已知条件为时间，日期，；可求得最短杆影对应时间为，进而求出视频拍摄经度设L （t ）= 再利用相应时间与表7.1中数据按二次曲线进行拟合，使用最小二乘法确定拟合方程为

.

545.162886.20.0766t =）t （L 2+-t 其中，L 为视频中杆影长度，t 为当地时间。

因为杆影最短的时刻即为当地正午时刻，故t=14.939，所以初始经度为

.92.1051512tmin -E 120J E ?=??-?=）（（东经）

其中，J 为视频所在地初始经度，tmin 为视频所在地杆影最短时刻。

之后代进模型一求得纬度为北纬38.8640°N 或者41.325°N 。从GPS 搜索得拍摄地大约在宁夏回族自治区银川市。

4.5 额外问题

如果拍摄日期未知，你能否根据视频确定出拍摄地点与日期？ 如果缺少日期条件，采取建立基于遗传算法的太阳影子定位模型，可以对经纬度矩阵增加日期纬度，引入3个0-1变量，对应着该矩阵的经度、纬度列向量、日期列向量，考虑到可能有若干个地点和日期，以各0-1变量求和不超过3作为约束条件，以竿影长和实际竿影长的的相似度作为目标函数，拟建立基于太阳影子定位的目标优化模型。

五、结论

通过以上四个问题的分析解答探讨，本文研究了杆子的影子在太阳照射下的影长与时间，日期、杆长以及杆子所在地的位置五者之间的相互关系，进而研究分析视频中物体的太阳影子变化，确定视频拍摄的地点和日期或者反推，这几个的关系如图。

问题一中，我们根据所知道的杆长，日期，时间，以及位置，求得影长。再问题二、三则分别在未知杆长，未知杆长和日期的情况下求解出杆子的地理方位，而问题四则利用视频中的影长变化来推断拍摄点位置。将这个模型推广，就有利于视频定位系统的建立。只要知道其他的三个或更多量就可以计算另外的量，这对本文得到的模型方法应用型较强。

六、模型的改进与推广

分析几个模型后可以发现，模型有一定缺陷性，下面逐一分析：

其一，我们没有考虑大气对太阳光的折光率（大气层从下到上可分为对流层、平流层、电离层和磁层四层。无线电波在大气层中传播时，由于在各层中的传播速度变化产生影响）问题，这对相关量，所测得的目标角度、距离、高度都存在大气折射误差，如太阳高度角的计算是有相当大的影响的。改善的方法是利用折射系数对太阳光进行修正，然后再求相关角度，最后代进模型。

第二是在问题二、三中，我们利用多项式拟合的方法有相当大的不准确性，在数据提供范围外，预测曲线与实际曲线的拟合度不是太高，导致一定的误差性，对此我们利用多次拟合度分析，并用实际曲线贴近拟合实际曲线，尽量减少误差量。

关于模型的推广，此模型适用于太阳影子定位的相关计算中，将这个模型推广，就有利于视频定位系统的建立。只要知道其他的三个或更多量就可以计算另外的量，因此，此模型还可以应用于对已知视频定位，求物体的高度，求相关日期时间等等问题。

参考文献：

[1] 司守奎，孙玺菁:数学建模算法与应用[S]: 国防工业出版社,2013.

[2]何银涛，张梅黄华.垂直单轴跟踪光伏支架方阵间距计算界面设计.太阳能报刊.2015.

[3]郑鹏飞,林大钧,刘小羊,吴志庭.基于影子轨迹线反求采光效果的技术研究[].广东:华南理工大学20023 7.

[4] 林根石,利用太阳视坐标的计算进行物高测量与定位[A].南京大学学报,1991.

[5] 吴济廉,影端轨迹周年变化的实践与分析.地理教学报，2013,50(5):9-20.

[6]张丽梅,赵建立,乔立山.用二值矩阵表示研究格矩阵的{1}-广义逆和{1,2}-广义逆.中国

知网.

附录：

杆子影子随时间曲线图：

%% 太阳赤角求解

clc,clear;%清空环境

Ts=(31-21+1)+30+31+30+31+31+30+21;

%计算lmd经度

lmd=180+90*((Ts-186)/90);

lmd=pi*lmd/180;%将角度转化为弧度

%计算直射点纬度角度K正弦值

sink=0.39775*sin(lmd);

k=asin(sink);%输出太阳直射纬度k

%纬度化为弧度

fai=39+(54/60)+(26/3600);

fai=pi*fai/180;

%% 时角问题求解

x=9:0.1:15;

fai1=116+(23/60)+(29/3600);

st=x+(fai1-120)/15

w=15*(st-12);

%计算太阳高度

w=w/180*pi;

y=sin(fai)*sin(k)+cos(fai)*cos(k)*cos(w);

z=asin(y);%输出太阳高度角

l=3./tan(z);

x=9:0.1:15;

figure(1);

plot(x,y,'-b')%太阳角变化

xlabel('时间点');ylabel('太阳角度变化');

figure(2);

plot(x,l,'r-')%绘制杆长变化图像

xlabel('时间点');ylabel('直杆的太阳影子长度');

clc;

clear;

as=0.39795*cos(0.98563*(295-173));

a=asin(as)

b=0.7302;

t=9:0.01:15;

w=pi/12*(t-12);

c=sin(b)*sin(a)+cos(b)*cos(a)*cos(w);

l=asin(c);

plot(t,l);%太阳高度角；

s=3./tan(l);

plot(t,s)

xlabel('t');ylabel('L');

title('直杆的太阳影子长度随时间的变化曲线');

利用附件1曲线拟合并用regress进行检验:

clear

clc

t=xlsread('附件1-3','D4:D24');

x=xlsread('附件1-3','B4:B24');

y=xlsread('附件1-3','C4:C24');

l=sqrt(x.*x+y.^2);

p=polyfit(t,l,2)

x1=linspace(9,15);

f1=polyval(p,x1);

plot(x1,f1)

hold on

k=polyder(p)

scatter(t,l,5)

A=polyfit(t,l,2)

z=polyval(A,t);

[b,bint,r,rint,stats]=regress(l,[ones(21,1),t,t.^2]); stats

plot(t,l,'k+',t,z,'r')

xlabel('t');ylabel('L');

title('直杆的太阳影子长度随时间的变化曲线');

利用附件2曲线拟合并用regress进行检验:

clear

clc

t=((12+42./60):0.05:(13+42./60))';

x=xlsread('附件1-3',2,'B4:B24');

y=xlsread('附件1-3',2,'C4:C24');

l=sqrt(x.*x+y.^2);

p=polyfit(t,l,2)

x1=linspace(6,18);

f1=polyval(p,x1);

plot(x1,f1)

hold on

k=polyder(p)

scatter(t,l,5)

A=polyfit(t,l,2)

z=polyval(A,t);