基于全局搜索收敛性太阳影子的定位

基于全局搜索收敛性太阳影子的定位

摘 要

本文通过对太阳影子定位技术的研究，针对题目中的各问题分别建立模型，解决了 由给出的太阳下影子的变化规律确定其所在地理位置和时间.

问题一要求我们分析影子的长度变化和哪些因素有关，据此建立模型.并用所建模 型确定题目中给定时间和长度的天安门广场上的直杆影子长度的变化曲线.我们分析影 子的变化主要与太阳的高度角有关，于是我们建立了以太阳高度角为核心的影长模型并 用Matlab 模拟出它影长随时间的变化轨迹.

问题二给出了已知日期的杆影运动轨迹，要求我们求出可能的地点.我们经过讨论发 现在预设了已知坐标方向与正北方向的偏差角度的情况下，可以知道太阳方位角函数与 时角的函数的差值是一个常量.我们建立了一个以方位角和时角之和方差最小为目标的 单目标规划模型.我们在搜寻合适经纬度坐标时，首先限制搜索范围为该时刻的地球表面 可接受光的区域中.由于全局搜索时任意一个经纬坐标对应的函数是连续的，所以我们可 以进行离散点择优，所以我们设计了全局搜索算法.通过该算法我们得出了最优点位于 19°N，110°E（海南黎族苗族自治县附近）.

问题三要求我们根据给出的影子顶点坐标数据，推测出可能的测量地点与日期.为此 我们将问题二中原本为常量的赤纬角变异为一个关于日期的函数，同样使用全局搜索搜 索一年中的最优解，从而求出附件二的最优解位于41°N，79°E(新疆维吾尔自治区阿 苏克地区)，时间为 5 月 26 日或第 7 月 18 日；附件三的最优解 3 月 2 日或 10 月 8 日取 得最优值，其位置坐标为38°N，110°E(陕西延川县)，另一个最优解在 28°N,110°

E(湖南湘西).

问题四要求我们根据视频中图像的变化判断该视频的可能拍摄点. 我们首先先获取 视频帧，并且对每一帧图像检测出影子的轨迹点；然后确定多个灭点，并拟合出地平线；

拟合互相垂直的灭点，计算出仿射纠正和投影纠正矩阵. 进而还原出经过度量纠正的世 界坐标.然后我们又利用了问题二中的模型，反推视频拍摄点的经纬度.得到在 40°N，

117°E（在呼和浩特）点处的拟合效果最小. 所以我们认为该视频的最佳可能拍摄地为 呼和浩特市.

关键词：太阳高度角；赤纬角；目标规划；Matlab；全局搜索；纠正矩阵

1问题的重述 1. 1背景知识

现在的互联网时代，时代的进步催生技术的产生.对处理图片和视频越来越重要了，

如何确定视频的拍摄地点和拍摄日期是视频数据分析的重要方面，它不但可以应用到图 像的识别上还可以用在图像的追踪.太阳影子定位技术就是通过分析视频中物体的太阳 影子变化，确定视频拍摄的地点和日期的一种方法.

1 . 2要解决的问题

在四个问题中我们首先应当建立影长与哪些因素相关的模型，并需要把该模型应用 到题目中给出的天安门广场上北纬 39 度 54 分 26 秒,东经 116 度 23 分 29 秒）一个 3 米 高的直杆上，求出其影子长度的变化曲线.

问题二中要求我们根据某固定直杆在水平地面上的太阳影子顶点坐标数据，建立数 学模型确定直杆所处的地点.并且把该模型应用到附件 1 中给出的数据，并给出几个合 理的直杆所在地理位置.

问题三中要我们根据某固定直杆在水平地面上的太阳影子顶点坐标数据，建立数学 模型确定直杆所处的地点和日期.这个问题要求我们通过所建模型给出直杆的地点和日 期.并且要应用到附件 2.3 给出的和附件 1 相类似的数据上.推测该直杆所在位置和时间.

问题四给我们提供的是一个并无明显特征的时长 40 多分钟的一直杆的影子在晴天 随时间变化的视频.问题要求我们根据视频影像内容和直杆的长度建立能确定视频拍摄 地点的数学模，并通过所建模型给出若干个可能的拍摄地点.

2模型的假设

1.地球的自转速度忽略不计；

2.略地球的公转为椭圆轨道；

3.不考虑地形因素，海拔因素；

4.假设太阳光源为面光源.

3名词解释与符号说明 3.1 名词解释

1. 日出时间：是指太阳每天从东方地平线升起的时间.

2. 太阳直射点：太阳直射点是地球表面太阳光射入角度为90 度的地点，它是地心 与日心连线和地球球面的交点.太阳直射点所在的经线的地方时为正午 12 时.

3. 时角：日面中心的时角，即从观测点天球子午圈沿天赤道量至太阳所在时圈的 角距离.

4. 太阳高度角：太阳高度角，是指某地太阳光线与该地作垂直于地心的地表切线

的夹角.

5. 太阳方位角：太阳方位角即太阳所在的方位，指太阳光线在地平面上的投影与 当地经线的夹角，可近似地看作是竖立在地面上的直线在阳光下的阴影与正南方的夹角.

方位角以目标物正北方向为零，顺时针方向逐渐变大，其取值范围是0——360°.因此太 阳方位角一般是以目标物的 北方向为起始方向，以太阳光的入射方向 为终止方向，按 顺时针方向所测量的角度 .

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

3.2 符号说明

符号 意义

h ^{太阳高度角}

 赤纬角，即太阳直射纬度



A 太阳方位角

t 时角

Bzt 北京正午时间

Bt 北京时间

4 模型的建立与求解

建立影子长度变化模型

4.1.1 太阳直射点的纬度在空间直角坐标系中参数方程

假设地球绕太阳公转的速度是匀速的，以地球为参照物，将地球绕太阳公转运动转 换为太阳绕地球做匀速旋转运动，如图1 所示。

图 1 太阳“绕”地球旋转

以地心为原点，黄道面为xOy 平面建立平面直角坐标系，假设太阳“绕”地球旋转轨 迹为圆，设夏至日太阳直射点坐标为(1,0)，秋分日为(0, ，冬至日为1) (1,0)，春分日 为(0,1)，如图2 所示。

图 2 太阳直射点的平面直角坐标系

把365 天平均分成 360 份，每一天占 365

360 ，设s 表示观测当天距离 1 月 1 日的天数，

夏至日距离1 月 1 日的天数为 175 天，根据建立的平面直角坐标系可得目标当天太阳直

射点在平面直角坐标系的坐标方程为











 



  



 



365 175 2

365 sin 2

365 ) 175 2

365 cos(2







 y s

x s

(1)

以地心为原点，黄道面为xOy 平面，过原点垂直于黄道面的直线为 z 轴建立空间直 角坐标系，如图3 所示。

图 3 太阳直射点的空间直角坐标系

则太阳直射点在空间直角坐标系的参数坐标方程为



















 



  



 



0

365 175 2

365 sin 2

365 ) 175 2

365 cos(2

z y s x s









(2)

4.1.2 观测地点地理纬度在空间直角坐标系中参数方程

以地心为原点，黄道面为xOy 平面，过原点垂直于黄道面的直线为 z 轴建立空间直 角坐标系，如图4 所示。

图 4 天安门广场所在纬度的空间直角坐标系

假设地球半径为1，则天安门广场所在纬度的轨迹参数方程为



















5 . 66 sin sin 5 . 66 cos cos cos

cos sin

5 . 66 cos sin 5 . 66 sin cos cos

















z y x

(,^为参数) (3)

4.1.3 观测地点的太阳高度角

设太阳直射点的纬度为为₁，观测地点地理纬度为₂，当纬度在南半球时记为负，

当纬度在北半球时记为正，则太阳直射点和地心的连线与赤道面的夹角为₁ ，观测地 点和地心的连线与赤道面的夹角为₂ 。观测地点的太阳高度角为，当观测地点与太 阳直射点不在同一半球时，可以得到 90 ₁₂ ，如图5 所示；当观测地点与太 阳直射点在同一半球时，可以得到 90 ₂ ₁ ，以北半球为例，如图6 所示。

图5 不在同一半球 图6 在同一半球

因此目标地点的太阳高度角可以由OA与 OB 的夹角₁ ₂ 表出，O 表示地

心，B 表示目标地点，A 表示太阳直射点，如图 7 所示。

图 7

4.1.4 影子长度变化曲线

根据已经求得的目标地点的太阳高度角计算得到直杆的影子，如图8 所示。

图 8

由图可知，影子长度

 arctan l h 。

4.2 模型求解

4.2.1 北京时间转换为目标地点的地方时

北京时间是指东经120 度的当地时，经度每隔 15 度，地方时相差一小时。当北京 时间12 点时，天安门广场(北纬 39 度 54 分 26 秒,东经 116 度 23 分 29 秒)的地方时不是 12 点整，而是 11 点 45 分 31 秒，太阳高度角有所不同，因此考虑将北京时间转换为目 标地点的地方时来计算。设天安门广场当地时与北京时间的差值为t，则





4  ^{' ''} 



t

得到北京时间9:00—15:00 天安门广场的地方时为 8:45—14:45。

4.2.2 确定参数 的取值范围

当天安门广场的地方时为0 点时， 0 ；地方时为12 点时， 180 ，因此地方

时从0 点开始，每过 1 小时，^增加15 ，由天安门广场的地方时为 8:45—14:45 得到参 数的取值范围为131.25 —171.25 。

4.2.3 直杆影子长度变化曲线

图 5 直杆影子长度变化曲线图

我们可以得到直杆影子最长为7.1518m,最短为 3.7788m。

4.4 模型检验

4.4.1 检验太阳直射点的纬度

由太阳直射点的回归运动规律可知，太阳直射点在地表移动的轨迹图接近于正弦函 数，如图1 所示，

图 6 太阳直射点在地表移动的轨迹图 1

对正弦函数解析式进行变形得到太阳直射点纬度的变化公式1

 

365 80 sin2

5 .

23  

 s

L 

(4) (其中 s 表示目标当天距离 2015 年 1 月 1 日的天数，80 表示春分日距离 2015 年 1 月 1 日的天数。)求得 2015 年 10 月 22 日太阳直射点的纬度为南纬 12.47°。

由于太阳直射点在地表移动的轨迹图并不完全符合正弦函数，为更精确地得到太阳

直射点的纬度，以秋分日为起点，对图1 进行调整得到图 2，

图 7 太阳直射点在地表移动的轨迹图 2

利用公式

 

365 166 sin2

5 .

23  



 s

L 

(5) 其中s 表示目标当天距离 2015 年 1 月 1 日的天数，166 表示秋分日距离 2015 年 1 月 1 日的天数。

计算可得2015 年 10 月 22 日太阳直射点的纬度为南纬 11.25°。

根据影子长度变化模型可得2015 年 10 月 22 日太阳直射点的纬度为南纬 11.61°，

误差最大为0.36°，因此影子长度变化模型中对太阳直射点的纬度计算方法可行。

4.4.2 检验影子长度变化曲线(程序见附件二)

根据文献[1]，通过解天文球面三角方程可以得到太阳高度角计算公式如下：

t cos cos sin sin

sin

sin       (6) 其中表示太阳高度角， 表示太阳赤纬，表示观测地地理纬度，t 表示太阳时角(以 当地正午为0 ，上午为负，每小时 15 ，下午为正，每小时15 )

根据公式(6)可得太阳赤纬为南纬11.61 ，即 11.61 ，计算可得观测地地理纬度 9

.

39

 ，太阳时角t 的范围为 8:45—14:45 的太阳时角，则直杆影子长度变化曲线为

可以得到直杆影子最长为7.6187m,最短为 3.7735m。

比较两者的结果，如表1

图中的曲线一表示由影子长度变化模型得到的影子长度变化曲线，曲线二表示由太阳高 度角计算公式得到的影子长度变化曲线，从图中我们可以发现两条曲线基本吻合，因而 本文所建立的影子长度变化模型具有可行性。

上述方法是利用相关变量的参数方程解决了问题一. 将参数方程整合成一个表达式 如下：

图9：太阳高度角与杆长和影长之间的关系 太阳高度角计算公式：

sinhsin sin  cos cos cos  t ① 从示意图中我们可得到如下关系：

3 1 sin2

tanh sinh

L h

x   ②

日出时间计算公式为：

360( 9)

180 cos[ tan( 23.31 ) tan ] 365.24

15

arc ^ n^ 

  

③ 算法流程图：

将当地坐标北纬39 度 54 分 26 秒,东经 116 度 23 分 29 秒代入③式，可求得当地北 京时间9：00-15：00 的时角 t 的变化范围t ( 44.5 , 45.5 ) , 将变化范围代入①式可以求 得这段时间内的太阳高度变化量，再根据②式即可求出影长.

我们用Matlab 软件编写程序对影子长度随时间的变化模拟，得出下图：

-50 -40 -30 -20 -10 0 10 20 30 40 50 3.5

4 4.5 5 5.5 6 6.5 7

t时 角 (对 应 时 间 从 9点 到 15点 )

杆的影长

图10：杆的影长和时间的关系

图5 中显示杆长最小的时间并非是 12:00，经计算该点的正午时间为 11:58 .近乎相 同.而且该模拟图形并非完全对称的，这也正印证了我们的猜想.在误差允许范围内结果 是可以接受的.

由①②③可以看出，影响影子长度变化的主要变量有纬度和日期决定.

由于时间来不及，我们没有对两种作比较.

4.2 问题二的分析与求解

对此问题我们建立了两个模型分别求解并对比，并对二者结果进行择优选择:

4. 2. 1模型一：利用方位角和时角差为常数而建立的最优化模型.

该问题需要我们解决如何由影子的顶点坐标数据判断直杆所在的地理位置.由于我 们并不知道在测量该数据时的地理方向相吻合，所以我们假设其X 坐标与正南北线的夹 角_.

图11：坐标修正后的各角度之间的关系

注：

Y

从图11 可得：

( , ) ( )

A   t   

因为偏移角是不变的，所以在该地的任意时刻，其方位角的函数与时角的函数之差 应是一个常数，也就是说我们求得目标地点的方位角与时角之差为一个常数.

由此关系式我们建立了一个单目标规划模型，其约束条件为：

cos sin sin cosh

sin sin sin cos cos cos 15 [ 120 ]

. . 15

( 180,180) ( 90,90) (14 : 700,15 : 700)

A t

h t

t Bt Bzt

s t

t



   







  



 

 

    



  

  

 



要使得到的目标地点数据与真实地点的数据非常吻合，我们需要将两者的方差进行 比较，越符合的数据其方差越接近0.所以我们建立了一个以方位角和时角之和方差最小 为目标的单目标规划模型.其目标函数为：

min(var[ ( , ) ( )]) Y  A   t

我们用MATLAB 软件该模型进行编程求.算法流程图如下：

运行程序后得到的图形如下所示：

0 20 40 60 80 100

0 50 100 150 200

0 1 2 3 4

x 10⁷

经 度 纬 度

方差的倒数

图12：模型一经纬度与方差倒数的关系

由于优化的结果用方差很小，所以我们选用方差的倒数，即高度越高越符合要

求，示中我们得到在 17°N，111°E（海南黎族苗族自治县附近）及其附近的点处的方 位角和太阳高度角之和的方差最小，即是该直杆最合理的位置.

解决问题的过程中我们发现此问题的解决方法并不是唯一的，近乎相同的约束条件 而目标函数不同也很合理，为了比较两模型之间的优劣我们又对模型二进行问题的求解.

模型二：利用同一地点连续时间内太阳高度角正切值变化的唯一性建立的最优化模 型.

由测量得到的( yx, )我们可以求出影长X ，假设杆长为 L ，则可以得到：

x

 l tanh 所以对于任意一个已测点都有：

1 1 1

2 2

2

tanh tanh

L

x x

L x

x

 

由此我们就得到了这段时间内该处太阳高度角正切tanh的变化率，日常生活经验告 诉我们：在杆长相同的情况下，由于经纬度的不同，不同地点在同时刻的太阳角高度的 变化曲线是不相同的，所以tanh的变化率对应唯一的地点.即需要求得的目标地点，由 此我们得到约束条件：

sin sin sin cos cos cos . . 120

15 [ ]

15

h t

s t t Bt Bzt

   



 





    



要使目标地点和实际地点一致，必须要求其_tanh的变化率在这一时间段内是一致的，

对此我们将目标地点的变化率与实际的变化率做方差比较，越小越符合要求，所以目标 函数：

20 2 2

1

(tan _i tan _i)

i

Y h x







使用MATLAB 软件编写算法进行求解，算法的框图和模型一的类似，程序见附 录，运行程序后的如下图：

0 20 40 60 80 100 0

50 100 150

2000 5 10 15

x 10⁶

经 度 纬 度

方差和的倒数

图13：模型二经纬度与方差倒数的关系

由图示我们得到在19°N，110°E（海南黎族苗族自治县附近）处的方位角和太阳 高度角之和的方差最小，即该处和该处的临近点是该直杆最合理的位置.

三、问题三的分析与求解

由于这个问题我们不知道时间和地点，如果我们知道时间的话，问题三就可以 简化为问题二，从而得到最优解，为了简化题目，我们引进了赤纬角δ关于日期变化的 函数y(n)，其推导过程如下：

夏半年：从春分日(3 月 21 日前后)到秋分日(9 月 3 日前后)太阳在黄道上运转 180°.

设从春分日开始，视太阳运行了n 天，则 n 天运行了λ经度：

(180* ) 186 n

  将其带入到太阳直射点纬度的计算公式中：

arcsin[0.39775 sin(180 ) ] 186 n

   (n[0,186])

冬半年：自科分日（9 月 23 日前后）至冬至日（12 月 22 日前后）太阳在黄道上运 转90°.设从春分日开始，视太阳运行了 n 天，则 n 天运行了

180 (n 186)

     再带入到太阳直射点纬度的计算公式中得：

sin 3.9775 sin[180 (n186) ] arcsin[0.39775 sin(n 186) ]

     (n[186,276])

自冬至日至次年春分日太阳在黄道上运转90°，设从春分日开始，视太阳行了 n 天，则n 天运行了_经度：

sin 0.39775 sin[270 90 ( 276) ] 89 n

     

带入到太阳直射点纬度的计算公式中得：

arcsin{0.39775 cos[90 ( 276)] } 89 n

     (n[276,365])

将上题中的赤纬角常数用这个函数代替，设立一年中每天的循环，以同样的方 式运行程序，可得到365 天中每天的运行结果，选取每个运行结果中的最优解，则一共 有365 个最优解，为了更直观的表明两个模型的相似度，我们同时比较两模型的图形:

0 50 100 150 200 250 300 350 400

0 0.5 1 1.5 2 2.5 3 3.5 4 4.5x 10⁷

天 数

该天最佳方差和的倒数

0 50 100 150 200 250 300 350 400

0 2 4 6 8 10 12 14x 10⁷

天 数

方差的倒数

图14：附件 2 求得当天最优值与天数的二维表

可以看出在第146 天(5 月 26 日左右)和第 200 天(7 月 18 日左右)处方差最小，故可 以认定此时间为要求得的时间，查看这两个日期的方差表如下图所示：

0

20 40

60

60 80 100 120 1400

1 2 3 4 5

x 10⁷

经 度 纬 度

方差和的倒数

0

20 40

60

60 80 100 120 140

0 0.5 1 1.5 2

x 10⁶

经 度 纬 度

方差的倒数

图15:附件二中符合条件的两个日期的经度纬度方差三维图

图像在41°N，79°E(新疆维吾尔自治区阿苏克地区)取得最优值，故可判定附件二 所描述的地点位于5 月 26 日或者 7 月 18 日左右的新疆阿苏克地区.

用 同 样 的 方 法 处 理 附 件 三 ， 得 到 图 像 如 下 图 所 示 ：

0 50 100 150 200 250 300 350 400 0

1 2 3 4 5 6 7 8x 10⁷

天 数

方差和的倒数

0 50 100 150 200 250 300 350 400

0 2 4 6 8 10 12 14 16x 10⁷

天 数

方差的倒数

图16：附件3 求得当天最优值与天数的二维表

0

20 40

60

60 80 100 120 1400

1 2 3 4 5

x 10⁷

经 度 纬 度

方差和的倒数

0

20

40

60

60 80 100 120 1400

1 2 3 4

x 10⁷

经 度 纬 度

方差倒数

图17:附件三中符合条件的两个日期的经度纬度方差三维图

由上图可知道图像在第61 天(3 月 2 日)和第 281 天(10 月 8 日)取得最优值，其位置 坐标为38°N，110°E，即陕西延川县.

四、问题四的分析与求解

看完问题给出的的这段视频我们发现在视频中唯一的明显的变化的就是直杆的影 子，由于影子的变化较为缓慢，我们首先对视频做了压缩处理.然后我们用 Matlab 软件 对视频进行了桢抓取工作.由于影子的变化较为缓慢，我们就做了简单的处理我们仅仅抓 取了视频中20 个间隔相等的桢作为处理样本.并且我们假设摄像机在摄像时它的

镜头是水平的.

在确定好图片后，我们用 Photoshop 对图片进行了灰度处理，其中处理后的图片如 下：

我们对给出的数据进行简单的拟合，如图为拟合轨迹.直杆的影子轨迹已经得 出来，那么我们就可以直接代入到问题二的模型中，由影子的轨迹直接推出直杆可能地 理位置.

图 18：影子轨迹的拟合曲线

表1：图片对应质感的顶点的修正后的二维坐标

消隐点与消隐线：

透视投影的一个主要特征是可以延伸到无穷远的物体在照片平面中可以获得它的有限延伸范围.

如，一条无穷远的直线可以投影到一条某消隐点为端点的直线.通常情况下，消隐点的确定可以通过 照片中一组平行的线段计算获得.通过消隐线的确定我们就可以得到相对于照片平面平行的平面了.

针对于本题视频后的两个方形的箱子恰好为我们得出消隐点和消隐线提供了方便.我们绘出的 点如下图所示

图19：消隐点与消隐线和水平线

将平面上的二维的点转换为空间上也可以确定的点.

[ , ]*

n

n n

u

v k r t i

x y z

 

   

    

   

    

n

n n

u v p i

x y z

 

   

    

   

    

, , n

n n

u

v p

i

x F T y

z

 

   

    

   

    

顶点 1 2 3 4 5 6 7 8 9 10

x 967 946 946 936 925 916 909 896 887 878 y 47 47 48 50 51 52 52 53 53 57 顶点 11 12 13 14 15 16 17 17 19 20 x 869 863 851 841 832 824 816 805 797 791 y 58 59 59 60 60 61 63 63 64 64

由上述关系式我们得到了一组影长数据：

表2： 所求的影长

于是我们可以分别推出我们抽取出来的图像的直杆上顶点随时间的变化规律，随即 我们就可得出影子的轨迹，从而我们就可以将数据直接带入到问题二的模型及算法当中，

得出最优的可能位置.

将拟合后的轨迹和模型二中的方法相对比，进行全局搜索，我们对任意一个经纬度 点和要求进行匹配.得到的坐标为得到在 40°N，117°E（在呼和浩特）点处的拟合效 果最小.所以我们认为该视频的最佳可能拍摄地为呼和浩特市.

组数 1 2 3 4 5 6 7 8 9 10

影长 2.5848 2.5728 2.5598 2.5470 2.5344 2.5221 2.5101 2.4988 2.4868 2.4755

组数 11 12 13 14 15 16 17 18 19 20

影长 2.4645 2.4537 2.4431 2.4327 2.4226 2.4127 2.4029 2.3935 2.3842 2.3751

5 模型的误差分析与检验

5.1 误差分析

在我们处理问题的时候我们处理问题时是基于很多假设前提下，如第一问中我们认为在地球围 绕太阳转动的过程中赤纬角是恒定不变的.而实际上它不但会变化还会直接影响方位角。

5.2 典型相关分析的检验

1， 误差来源：由于是离散点的全局搜索会导致搜索结果可能导致收敛于局部最优，从而导致偏差.

0 20 40 60 80 100

0 50 100 150 200

0 1 2 3 4

x 10⁷

经 度 纬 度

方差的倒数

图20：线性拟合前（左）线性你和后（右）

从多次拟合结果来看，拟合后的函数和离散点的收敛性一样，基本上认为函数的峰值和全局搜 索的一致.

2， 第四问中模型数据测量的误差：

误差来源：由于在测量图片影长时，采用的是比较传统的测量方式，即matla 图片上手动测量工具.

使测量范围在±5 内随机变化，得到多次的所得 21 组实际影长的变化：

0 5 10 15 20 25

2.35 2.4 2.45 2.5 2.55 2.6

组 数

实际影长

1 图21：测量数据有微小变化时对应的影长变化

6 模型的评价

在问题一中在解决太阳高度角的问题上，我们选择使用太阳高度角和时角之间的联 系作为突破口.巧妙的将经纬度和太阳时角联系在一起.

在问题二的处理当中我们并非简单的使用一种方法或模型去解决问题，我们使用了 完全不同的思路构造出不同的单目标规划.亮点在于我们将两者之间的处理结果进行了 对比，发现两者结果近似相等，这就更加印证了我们结论的正确性.

在论文当中我们经常的用到利用时间来推算太阳直射点纬度，而地球公转的周期是 一个回归年，现行功力的历年是历日的整数被，它和回归年并不相等.这直接导致了日 期和黄经并不完全相对应.而我们的模型在处理时间时是按季节分段时行计算，确保了 推算的相对准确.

在问题四中，我们应用了图像学经常使用的技术，对视频进行了简单且少失真处理.

并且将由摄像机拍摄的二维图像经过处理转换为真实的三维坐标.减少了误差.

论文中最明显的缺点莫过于经常使用的近似处理问题,导致精度不太高。

在第四问数据的处理过程中有部分关键数据严重依赖人工完成，具有较大的偶然性。

所以得到的结果可能会有较大的偏差。

在第四问的处理中我们在处理直杆并不是竖直的放在水平面上的问题上我们完全 的依赖于透视变换。这种结果的准确性我们并不知道。

参考文献

[1] 蒋洪力，太阳直射点纬度的数学推倒与分析.数学通报[J],第 46 卷 39-40,2007 年。

[2] 武琳，基于太阳阴影轨迹的经纬度估计技术研究，天津大学[J],201012。

[3] 袁信，张珣，基于太阳方位算法的非 GPS 定位系统的研究与设计，中国水运[J] 第 11 卷 68-70,2011 年 11 月。

［４］卓金武 《Matlab 在数学建模中的应用》，北京市海淀区，北京航空航天大学出 版社，2011 年 4 月一版。

附 录

模型一所用代码

function t=sj(bt,jd,bzt)%bt 北京时间 bzt 北京正午时间 jd 经度 l=length(bt);

jd1=jd*pi/180;

t=15.*(bt+ones(1,l)*(bzt-(120-jd)/15));

end

function x=jsxl(year,m,d)

a=[31 28 31 30 31 30 31 31 30 31 30 31];

x=0;ii=1;

if((rem(year,4)==0|rem(year,100)~=0)&(rem(year,400)==0)) yearl=366;a(2)=29;

else

yearl=365;

end

while ii<=m-1;

31+

x=x+a(ii);

ii=ii+1;

end x=x+d;

function y=time2num(h,m,s) m=m/60;

s=s/3600;

y=h+m+s;

function y=cwj(n)%计算赤纬角

x=n-80;

% if x<0

% x=365+x;

% end if x<=186

y=asin(0.397775*sin(x*(180/186)*pi/180))*180/pi;

else if x<=276

y=-asin(0.397775*sin((x-182)*pi/180))*180/pi;

else

y=-asin(0.397775*cos(90/89*(n-276)*pi/180))*180/pi;

end end end clc clear

%第一问的脚本文件

tianshu=jsxl(2015,10,22);

cwd=cwj(tianshu);

st=sunuptime(39.9,tianshu)-0.15;

t=-44.5:0.1:55.5;

sinh=Sinh(cwd,39.9,t);

x=3.*sqrt(1-sinh.^2)./sinh;

plot(t,x);

第二问的脚本文件

% 注：首先将附件 1 的 x,y 的坐标放入 x1,y1;

%%%%%首先将附件的 x,y 的坐标放入 x1,y1;

clc

tanshu=jsxl(2015,4,18);%计算 2015 4 18 号的序列天数 bt1=time2num(14,42,0);%将 14:42 换算为实数

bt=bt1:0.05:bt1+1;%因为每隔 3min 中即为 0.05 小时 for jd=0:180;

for wd=0:90;

h=length(jd);

l=length(wd);

bzt=time2num(11,59,26);%北京正午时间 11:59:26 根据第一题程序带入 4 月 18 号可得 t=sj(bt,jd,bzt);%求时角

zwd=10.9;%据上题可得

sinh=Sinh(zwd,wd,t);%太阳高度正弦 a=cosfwj(sinh,t,zwd);%求方位角的正弦 dd1=asin(a)';%dd1 为方位角

dd2=atan(x1./y1);%dd2 为目标角度 var(dd1-dd2);

jz(jd+1,wd+1)=var(dd1-dd2);

end end

%画图

[x,y]=meshgrid(0:90,0:180);

mesh(x,y,1./jz);

第三问 附件 2 的脚本文件

clc

% 注：首先将附件 2 的 x,y 的坐标放入 x2,y2;

bt1=time2num(12,41,0);%将附件二 12:41:00 时间换算为实数 bt=bt1:0.05:bt1+1;%因为每隔 3min 中即为 0.05 小时

for n=1:365;

zwd(n)=cwj(n);

end

for jj=1:365;

for jd=70:135;

for wd=0:55;

bzt=time2num(11,59,26);%北京正午时间 11:59:26 根据第一题程序带入 4 月 18 号可得 t=sj(bt,jd,bzt);%求时角

sinh=Sinh(zwd(jj),wd,t);%太阳高度正弦 a=cosfwj(sinh,t,zwd(jj));%求方位角的正弦 dd1=asin(a)';%dd1 为方位角

dd2=atan(x2./y2);%dd2 为目标角度 var(dd1+dd2);

jz(jd-70+1,wd+1)=var(dd1-dd2);

jz3{jj}=jz;

end end end

%

for ii=1:365

da(ii)=min(min(jz3{ii}));

jl=ii;

end

[x,y]=meshgrid(0:55,70:135);

mesh(x,y,1./jz{146});

% switch floor(n/90)+1;

% case 1;

clc

%第三问 2 的脚本文件

% 注：首先将附件 2 的 x,y 的坐标放入 x3,y3;

bt1=time2num(13,09,0);%将附件三 13:09:00 时间换算为实数 bt=bt1:0.05:bt1+1;%因为每隔 3min 中即为 0.05 小时

for n=1:365;

zwd(n)=cwj(n);

end

for jj=1:365;

for jd=70:135;

for wd=0:55;

bzt=time2num(11,59,26);%北京正午时间 11:59:26 根据第一题程序带入 4 月 18 号可得 t=sj(bt,jd,bzt);%求时角

sinh=Sinh(zwd(jj),wd,t);%太阳高度正弦 a=cosfwj(sinh,t,zwd(jj));%求方位角的正弦 dd1=asin(a)';%dd1 为方位角

dd2=atan(x3./y3);%dd2 为目标角度 var(dd1+dd2);

jz(jd-70+1,wd+1)=var(dd1-dd2);

jz3{jj}=jz;

end end end

%

for ii=1:365

da(ii)=min(min(jz3{ii}));

jl=ii;

end

[x,y]=meshgrid(0:55,70:135);

mesh(x,y,1./jz{61});

% switch floor(n/90)+1;

% case 1;

模型二所用代码：

bt1=time2num(12,41,0);%将 14:42 换算为实数 bt=bt1:0.05:bt1+1;%因为每隔 3min 中即为 0.05 小时 dd2=sqrt(x1.^2+y1.^2);

for ii=1:20

bz(ii)=dd2(ii+1)./dd2(ii);

end

for n=1:365;

zwd(n)=cwj(n);

end

for jj=1:365;

for jd=75:135;

for wd=0:55;

bzt=time2num(11,59,26);%北京正午时间 11:59:26 根据第一题程序带入 4 月 18 号可得 t=sj(bt,jd,bzt);%求时角

% zwd=11.9;%据上题可得

sinh=Sinh(zwd(jj),wd,t);%太阳高度正弦 dd1=asin(sinh);%dd1 为高度角

for ii=1:20

bz1(ii)=tan(dd1(ii))./tan(dd1(ii+1));

end

sum(sqrt((bz-bz1).^2));

jz(jd-75+1,wd+1)=sum((bz-bz1).^2);

jz3{jj}=jz;

end end end