转载 2015A国赛优秀论文

1
太阳影子定位
摘要
本题意在通过分析实际测量的竖直杆影子随时间变化的坐标或视频数据,确
认竖直杆所在的经纬度以及数据对应的记录日期。
第一问要求分析画出已知杆长、地理位置和日期的竖直杆的太阳影子变化曲
线,并分析影子长度关于各个参数的变化规律。本文利用天文学公式,建立了该
竖直杆太阳影子关于时间变化的计算求解模型,得出影长随时间的变化曲线如下。
此后,本文使用控制变量法,分别对经纬度、日期、杆长及当日时间这五个参数
对影子长度的影响展开讨论。经过误差分析并改进算法后,得到在11 点59 分24
秒,取到最短的影长为3.6638 米
第二问要求分析某竖直杆的太阳影子顶点坐标数据,确定其地点。本文沿用
第一问的天文学公式,以每一对时刻理论与实测的影长之比的差的平方和最小为
首要目标,以每一对时刻理论与实测方位角差值之差的平方和最小为次要目标,
建立多目标规划模型,采用分层求解法,遍历求解得到最能满足条件的位置为北
纬19 度,东经109 度,位于我国海南省。
第三问要求分析两组某竖直杆的太阳影子坐标数据,确定其地点及数据记录
日期。与第二问类似,在添加日期这一变量后,建立与第二问相同的多目标规划
模型,采用分层求解法遍历求解。得到附件二和附件三数据对应的最能满足条件
的位置分别为:
纬度 经度 日期 大致地点
附件二 40 79 5 月25 日 新疆图木舒克市
附件三 33 106 10 月31 日 陕西,汉中广元之间
第四问要求分析竖直杆的影子变化视频,分别在日期已知和未知的情况下求
出其所在地点与拍摄日期。本文首先利用动态追踪技术,找出影子顶点在每一时
刻的像素坐标,计算出每时刻的太阳高度角。进而利用天文学公式中太阳高度角
与日期、纬度、经度的确定关系,分别在日期已知和未知的情况下建立计算求解
模型,考虑到模型存在系统误差,难以得到完美解,故将问题转换为非线性规划
问题遍历求解。得到最能满足条件的位置与日期为:
纬度 经度 日期 大致地点
日期确定 43 115 7 月13 日 内蒙古锡林郭勒盟西南
日期未定 44 113 6 月25 日 内蒙古锡林郭勒盟西
本文亦对文中使用的天文学公式的误差及所求出的结果进行了讨论,在给出
最优解之外还讨论了其他解存在的可能性。
关键词:计算求解模型、控制变量法、多目标规划、分层求解法、误差分析
本资料由:数模之家(春哥带你学建模)收集并分享,需要更多数学建模资料添加微信公众号:shumohome
数模之家
2
1. 问题重述
1.1 问题背景
现代科技的发展使得人们能够更为方便地记录高质量的视频文件。在分析视
频材料时,有时需要确定视频的拍摄地点及日期,而利用天文学知识,对视频物
体中的太阳影子变化进行分析是确定视频拍摄地点及日期的一种方法。
1.2 相关信息
本题给出了三组某处某固定直杆在水平地面上太阳影子的顶点坐标数据,每
组数据皆包含 21 个等间距时刻(北京时间)直杆太阳影子的x, y顶点坐标。其
中,第一组数据还额外包括了数据对应的日期。坐标系以直杆底端为原点,水平
地面为x, y平面,直杆垂直于地面。
除三组数据以外,本题亦给出了一根直杆在太阳下的影子变化视频,视频中
清晰的记录了2015 年7 月13 日8 点54 分06 秒到同日9 点34 分36 秒某地高
为2 米的直杆的太阳影子变化过程。以上信息中,各直杆的地理位置皆未知。
1.3 需要解决的问题
1) 建立太阳影子变化的数学模型,分析影子长度关于各个参数的变化规律并应
用这一模型画出2015 年10 月22 日北京时间9:00-15:00 之间天安门广场(北
纬39 度54 分26 秒,东经116 度23 分29 秒)3 米高的直杆的太阳影子长度
变化曲线。
2) 建立数学模型,根据题目提供的附件一中的影子顶点坐标数据,给出该直杆
所在的可能地点。
3) 建立数学模型,根据题目提供的附件二、三种影子顶点坐标数据,给出对应
直杆所在的可能的地点。
4) 分析题目提供的视频,确定视频拍摄地点的数学模型,分别讨论在拍摄日期
已知和未知两种情况下,能否确定视频的拍摄地点与日期,并给出该视频可
能的拍摄地点。
2. 模型假设
假设一:本文研究的所有对象所在的地面皆为平地。
假设二:本文研究的所有对象所处地海拔为0。
本资料由:数模之家(春哥带你学建模)收集并分享,需要更多数学建模资料添加微信公众号:shumohome
数模之家
春哥给大家定期发资料啦,所有符合要求的,会被标记“定期发资料”,然后春
给会定期发布内部课讲义、各种数学建模资料、以及春哥独家整理的内部资料。
资料主要包括:优秀论文、优秀论文框架分析、排版格式问题汇总、各种论文
模版、各种数学建模教材PDF 版、各种算法源代码以及源代码使用说明以及历
年国赛赛题分析课件等资料。
(1)VIP 学员、传课学员(VIP 学员有获奖保证)
(2)关注微信公众号:shumohome,并回复“分享者”按照要求操作即可。
注:资料包括春哥分析整理的资料等,因为资料性价比非常高,所以请大家不要
分享给别人。
有其他疑问可以添加QQ 群:598397784 咨询
郑老师
郑老师被学员亲切的称为春哥,精通数学建模思路分析,数学建模思维培养,建
模论文写作等方面。所指导过的学员70%以上获得了省级二等奖及以上的奖项,

春哥带你学建模的所有课程按照数学建模学习四步走进行授课:
(1)排版与格式——在国赛评阅中,评阅用时非常短,阅卷者会首先看一篇论
文的排版与格式是否正确。虽然说评阅过程摘要最重要,但是本课程会提供摘要
模板,可以直接在比赛中掌握。
(2)论文与摘要——论文的重要程度,不言而喻,如何写论文,是非常重要的,
本课程从两个方面讲解如何写数学建模论文,一个是完全没有思路如何编凑论文,
一个是如何写出一篇完美的论文。
(3)解题思路——依照讲师整理的算法模型简介,以及优秀论文,对数学建模
赛题进行分析,让学员学会对于各种类型题目都能迅速分析出最佳算法。
(4)方法积累——讲师整理各类算法及其在数学建模中的需要掌握程度,让你
在最短的时间掌握一种算法,例如两分钟掌握层次分析在数学建模中的应用
3
3. 符号说明
符号 符号说明
i 
时刻i 的太阳高度角
 太阳赤纬
d 代表天数,当日期为 1 月 1 日时,d=1
 目标对象所在地的纬度。为正数时代表北纬,负数为南纬
l 目标对象所在位置经度。为正数时代表东经,负数为西经
i t 时刻i 下研究对象对应时角
L 直杆长度
yi L 直杆在时刻i 时的影长

yi L 直杆在时刻i 的理论影长
i 
直杆影子在时刻i 被测得的方位角

i 
直杆影子在时刻i 的理论方位角值
1 f 每一对理论的影长之比的差的平方和
2 f 方位角差值之差的平方和
ls t 本地太阳时间
lt t 本地时
i0 t i 时刻i 时子午线当地标准时间
E 时差
4. 问题分析
问题(一)
第一问要求建立影子长度变化的数学模型,分析影子长度关于各个参数的变
化规律,并画出题目指定的直杆的太阳影子长度变化曲线。由于本题要求研究太
阳影子以定位,结合实际生活体验可知,对于一个自身长度固定的竖直杆,其影
本资料由:数模之家(春哥带你学建模)收集并分享,需要更多数学建模资料添加微信公众号:shumohome
数模之家
4
子长度不仅仅与其所在地理位置有关,还与地球自转情况、地球绕太阳公转的情
况有关。因地球自转与公转分别导致了昼夜交替及四季变化,故影子长度应该还
与对应的当天时刻和日期有关。
考虑到这些因素,可进一步搜索相关的文献,建立数学模型表现影长与地理
位置、时间的关系,通过研究关系式以得出题目要求求解的答案。
问题(二)
第二问要求根据某固定直杆在水平地面上的太阳影子顶点坐标数据,建立数
学模型确定直杆所处的地点。本问相对于第一问,直杆长度未知,但因知道太阳
影子顶点在一个连续时间段内的数据,故相当于知道了此时间段内影子的影长以
及其角度变化情况。经过第一问的分析与解答可知,影长与直杆长度、当天日期、
时间、经纬度此5 个因素有关;本问中,日期、时间皆为已知量,而特殊地,当
其他因素恒定时,直杆长度与影长呈正比关系,故可以研究两个时刻的影长比值,
消除直杆长度这一未知量。这一比值仅仅与直杆所在的经纬度有关,因而能够达
到减少未知数的目的。
除了分析影长以外,角度也是一个不可忽略的因素。虽然影长信息已经足以
解出未知数,但再增加角度信息,能够提高方程组冗余度,使得模型的稳健性和
抗噪性得以提升。
用来表示角度的量应当能够体现影子的方位,参考资料可知,因方位角能够
表现影子在地平面上的方向,使用方位角[1]这一指标较为合适。在本题中,虽然
每个时刻的方位角的值难以得到,但是两个时刻的方位角差却与坐标系自身无关,
可以根据两个时刻的影子顶点在同一坐标系中的相对位置计算。
因而,综合考虑影长和角度两个指标,以直杆所在经纬度为x, y时,计算所
得两时刻的影子长度比与对应时刻测得的影子长度比的差值最小为第一个目标
函数,以计算所得两时刻方位角差值与测得的方位角差值的差值最小为第二个目
标,建立多目标规划模型,找出可行的直杆对应的经纬度值。
问题(三)
第三问与第二问类似,不同的是日期未知,要求根据某固定直杆在水平地面
上的太阳影子顶点坐标数据,建立数学模型确定直杆所处的地点与日期。与第二
问相同,可综合考虑影长和角度两个指标,以直杆所在经纬度与日期分别为x, y, z
时,计算所得两时刻的影子长度比与对应时刻测得的影子长度比的差值最小为第
一个目标函数,以计算所得两时刻方位角差值与测得的方位角差值的差值最小为
第二个目标,建立多目标规划模型,找出可行的直杆对应的经纬度值与日期。
本资料由:数模之家(春哥带你学建模)收集并分享,需要更多数学建模资料添加微信公众号:shumohome
数模之家
5
问题(四)
第四问要求分析题目提供的视频,确定视频拍摄地点的数学模型,分别讨论
在拍摄日期已知和未知两种情况下,能否确定视频的拍摄地点与日期。因此,需
要先对附件四的视频进行处理,找出正确的太阳高度角在拍摄时间段内各个时刻
的值。
因第一、二、三问建模过程中,确定经纬度并不一定需要直杆长度以及影长
这两个量同时已知,需要的仅仅是由之计算得到的太阳高度角在不同时刻的值,
因而得到此信息后可将问题转换为与前三问类似的问题。在前三问的建模过程中,
使用了关联太阳高度角、日期、纬度、经度、及时间的等式,故可利用该等式建
立计算求解模型求解。
5. 数据分析
5.1 附件四视频处理
本题附件四是一根直杆在太阳下的影子变化视频,视频中清晰的记录了2015
年7 月13 日8 点54 分06 秒到同日9 点34 分36 秒某地高为2 米的直杆的太阳
影子变化过程。直杆的地理位置皆未知。为了使用该视频中的信息解答第四问,
需要对视频做一定的处理,找出各个时刻直杆影子尖端部分的变化。
图 5.1.1 附件 4 视频截图
图 5.1.1 为附件4 的视频截图。因已知竖直杆垂直于地面,观察该视频发现:
竖直杆的中心位于画面中央,可认为摄像机的轴线是与直杆垂直的。随着时间的
推移,可以用肉眼观察到影子的长度变短,因此可以推断该视频拍摄在当地正午
时间前,则可以作为左上角的时间是该地所在时区的时间的佐证。
此外,由于在视频所记录的时间里,直杆影子转动的角度十分小,而在图示
画面中,可认为影子所在的直线是和摄像机的视线垂直的,因此可以认为影子长
本资料由:数模之家(春哥带你学建模)收集并分享,需要更多数学建模资料添加微信公众号:shumohome
数模之家
6
度在透视矫正前后变化不大,因而不需要在处理视频的时候对视频做透视矫正,
即认为影子和直杆构成的平面垂直于摄像机的视线。
图 5.1.2 附件4 处理示意图
利用 Adobe After Effect 软件读取该视频,采用软件自带的动态追踪技术读取
视频每一帧影子顶部的像素坐标C、直杆顶端的坐标A 以及影子上的任一其他
点的像素坐标 B由此,即可得到同一坐标系内的两条向量BC, AC,进而解得两
向量的夹角大小,用于之后的各项运算。
图5.1.3 动态追踪操作示意图
动态追踪 C 点的具体操作如图所示,图中有一大一小两个白边方框,小方框
本资料由:数模之家(春哥带你学建模)收集并分享,需要更多数学建模资料添加微信公众号:shumohome
数模之家
7
中是追踪的目标即影子顶点及周围的地面。由于影子顶点具有固定的形状与较深
的颜色,小白框内的像素的颜色信息有一种特征明显的分布,故可对每一帧图像
寻找最符合的像素区域来进行追踪。而另一方面,由于移动是连续的,为提高追
踪精度,大方框规定了下一帧目标点出现的范围。对追踪参数进行适当调试,可
得出图示的轨迹。
实际操作中,应将选择较为靠近直杆底座的点作为B 点,如此便不需要每次
都移动B 点,而可将之大致固定在一个位置,便于计算。在6.1.1 的分析中将体
现,ACB即为太阳高度角。
6. 模型的建立与求解
问题(一)的求解
6.1.1 模型的分析
为建立合适的数学模型,首先须对垂直直杆在太阳照耀下成影情况进行分析。
可构建示意图如下:
图 6.1.1.1 竖直直杆成影示意图
如图 6.1.1.1 所示,平面XOY 代表大地,此处可将大地看作是平坦的。OA为
竖直直杆,其中O点为直杆和大地的接触点,OB 为竖直杆OA的影子。黄色有向
线代表了由太阳射来的光线。由此构成了竖直直杆成影的示意图。由于直杆长度
OA已知为 3 米,故只需知道ABO的值,即可求出影子OB的长度。其中,ABO
为太阳高度角。
参考有关太阳高度角的资料[2],可以得到太阳高度角的求解公式如下:
  arcsin sin sin  cos cos cos t 
其 中 ,  为太阳高度角,  为太阳赤纬,  的计算方式为[3] :
 
360
=23.45 sin 81
365
 d
 
   
 
,d 代表天数,当日期为 1 月 1 日时,d=1。 为竖
直杆所在地的纬度,t为时角。同时,由于AO  BO,故
2 2
sin
y
L
L L
 

, , y L L
分别为直杆长度和影长。
本资料由:数模之家(春哥带你学建模)收集并分享,需要更多数学建模资料添加微信公众号:shumohome
数模之家
8
通常,时角可由本地时换算得到。然而因本地时常常因地球自转及人为调整
影响,而与该地真实的本地太阳时不相同,故本文采用精度更高的调整后计算时
角的公式。参考相关文献可知[3],=15  12 ls t  t  ,ls t 为本地太阳时间,
60 ls lt t t

  。
lt t 为本地时。本地时的计算方法为0 15 lt
l
t  t  ,其中,l为目标所在位置经度,当
l 为负时代表所在位置为西经,为正时即东经。0 t 为子午线当地标准时间。E 为
时差(EoT)。值得注意的是,此处的时差并非生活中使用的区时时差,而是定义
为:E  9.87sin 2B  7.53cos B 1.5sin B,  
360
81
365
B  d  。由此,即可得出
精度较高的时角数值,其计算方式为:
     
  0
9.87sin 2 7.53cos 1.5sin 360
15 12 , 81
15 60 365
l B B B
t t B d
   
       
 
经过以上分析,可发现,当目标经纬度、日期、时间都已知时,可以求出相
应的太阳高度角。又因杆长已知,故可求出相应时刻的影长。因题目要求绘制出
2015 年10 月22 日北京时间9:00-15:00(即子午线标准时间1:00-7:00)之间天安
门广场(北纬39 度54 分26 秒,东经116 度23 分29 秒)3 米高的直杆的太阳影
子长度的变化曲线,故可写出太阳此地此日太阳影子长度与时间的关系式,从而
在给定的定义域中,画出此直杆影子长度的变化曲线。
6.1.2 模型的建立
根据 6.1.1 中的分析,可建立有关影长y L 与子午线当地标准时间0 t 的计算求
解模型如下:
其中 y L 为直杆的影长, L 是直杆长度,为3 米。 为太阳赤纬,计算方式
为:  
360
=23.45 sin 81
365
 d
 
   
 
,d 代表天数,当日期为 1 月 1 日时,d=1。
为竖直杆所在地的纬度,竖直杆在北半球时, 为正,否则为负。t 为时角,根
据6.1.1 中的分析,可以得出:
     
  0
9.87sin 2 7.53cos 1.5sin 360
15 12 , 81
15 60 365
l B B B
t t B d
   
       
 
 
 
2
0
1
1
sin sin cos cos cos
1,7
y L L
t
t
   
 
   
    

本资料由:数模之家(春哥带你学建模)收集并分享,需要更多数学建模资料添加微信公众号:shumohome
数模之家
9
0 t 为子午线当地标准时间,l 为目标所在位置经度,当竖直杆经度为东经
时,l 为正,否则为负。
可以看出,由于日期、杆长、经纬度皆为已知,所建立的模型实际上为y L
关于t0,定义域为1,7,解析式为
 2
1
1
sin sin cos cos cos
y L L
    t
 
   
    

函数关系式。
综上所述,第一问的模型为:
 
 
 
     
 
2 0
0
1
1 , 1,7
sin sin cos cos cos
360
=23.45 sin 81
365
9.87sin 2 7.53cos 1.5sin
15 12
15 60
360
81
365
y L L t
t
d
l B B B
t t
B d
   

 
    
    
  
    
  
    
      
  

  

6.1.3 模型的求解
采用 MATLAB 编程求解,可以得出直杆影子长度的变化曲线如下:
图 6.1.3.1 直杆的太阳影子长度变化曲线
图 6.1.3.1 表示了2015 年10 月22 日北京时间9:00-15:00(即子午线标准时
间1:00-7:00)之间天安门广场(北纬39 度54 分26 秒,东经116 度23 分29 秒)
3 米高的直杆的太阳影子长度的变化曲线,横轴北京时间(小时),纵轴为直杆影
本资料由:数模之家(春哥带你学建模)收集并分享,需要更多数学建模资料添加微信公众号:shumohome
数模之家
10
长(米)。
需要注意的是,由于本题中计算太阳赤纬、时角误差修正的公式并非完全精
确,因此本模型仍然存在一定的系统误差。在本问解答结束后,将会对本模型的
误差进行分析。
6.1.4 结果分析
根据本文的计算结果,直杆的太阳影子长度变化曲线如图6.1.3.1 所示,其在
11 点58 分30 秒时,直杆影长达到最小值3.8410 米;影长在9 点到11 点58 分
30 秒时随时间变短,此后随时间增长。
决定影子长度的除了当天的时间之外,还有杆长、日期、经度、纬度这四个
因素。由于杆长与影子长度呈正比例关系,故本文不多加讨论,而着重在日期、
经度、纬度三个方面给出讨论。为方便表现影子长度关于各个因素的关系,须控
制变量,作图如下:
图 6.1.4.1 直杆的太阳影子长度在一年中的变化曲线
图 6.1.4.1 是1 月1 日至12 月31 日北京时间12:00(即子午线标准时间4:00)
北纬39 度54 分26 秒,东经116 度23 分29 秒(上图)与南纬39 度54 分26 秒,
东经116 度23 分29 秒(下图)3 米高的直杆的太阳影子长度的变化曲线,横轴
为日期(1 代表1 月1 日,2 代表1 月2 日,由此类推),纵轴为直杆影长(米)。
可以看出,对处在北半球的直杆而言从1 月1 日到夏至日左右(6 月22 日),
影子的长度随着时间推移而变短,夏至日后,随着时间推移而变长,在冬至日左
右(12 月22 日)达到最长后又开始减小。而对于处在南半球的直杆,情况恰好
本资料由:数模之家(春哥带你学建模)收集并分享,需要更多数学建模资料添加微信公众号:shumohome
数模之家
11
完全相反。
图 6.1.4.2 修正前与修正后直杆影长与经度关系曲线
图 6.1.4.2 体现的是6 月22 日北京时间12:00(即子午线标准时间4:00)北纬
39 度54 分26 秒,处于180 度到本初子午线的3 米高的直杆的太阳影子长度的变
化曲线,横轴为经度(负数为西经,正数为东经),纵轴为直杆影长(米)。
第一幅图未经过任何处理,为修正前变化曲线。可以看出,在西经120 度左
右到东经10 度左右的区域,子午线标准时间为4:00 时太阳尚未升起或已经落下,
因此对其讨论直杆影长毫无意义。而在此二经度附近的区域,太阳刚刚升起或落
下,若当地大地较为平坦,则影子长度极大,会影响对其他经度区域影子长度的
观察,因此本文将影子长度大于30 米的区域统一设定为30 米后,绘制出修正后
的曲线。
本资料由:数模之家(春哥带你学建模)收集并分享,需要更多数学建模资料添加微信公众号:shumohome
数模之家
12
观察两幅图可以发现,在东经120 度附近,直杆太阳影子最短,而以120 度
为中心向东西延伸,影子长度不断变长,直到达到晨昏线分界线后,影子消失。
图 6.1.4.3 修正前与修正后直杆影长与纬度关系曲线
图 6.1.4.3 体现的是6 月22 日北京时间12:00(即子午线标准时间4:00)东经
116 度23 分29 秒,处于南纬90 度到北纬90 度的3 米高的直杆的太阳影子长度
的变化曲线,横轴为纬度(负数为南纬,正数为北纬),纵轴为直杆影长(米)。
第一幅图未经过任何处理,为修正前变化曲线。可以看出,在南纬84 度左右
到南纬90 度,北京时间12:00 时为黑夜,因此对其讨论直杆影长毫无意义。而
在此南纬84 度附近的区域,太阳刚刚升起,若当地大地较为平坦,则影子长度
极大,会影响对其他经度区域影子长度的观察,因此本文将影子长度大于30 米
的区域统一设定为30 米后,绘制出修正后的曲线。
观察两幅图可以发现,在北回归线处,直杆太阳影子最短,而以北回归线为
中心向东西延伸,影子长度不断变长,直到达到南纬84 度左右,影子消失。
综上所述,题目要求绘出的天安门广场3 米高竖直杆在10 月22 日9 至12
点的曲线如图6.1.3.1 所示。观察竖直杆与当日时间、日期、经度、纬度的关系可
以总结规律如下:
表 6.1.4.4 影子长度关于各个参数的变化规律


当日
时间
日期 经度 纬度









正午
前逐
渐变
短,
正午
北半球的直杆在去年
12 月22 日至6 月22
日逐渐变短,6 月22
日至当年12 月22 日
在有日照的区域
内,某一经线上
的太阳影子最
短,以此经线为
在有日照的区域
内,某一纬度上
的太阳影子最
短,以此纬线为
本资料由:数模之家(春哥带你学建模)收集并分享,需要更多数学建模资料添加微信公众号:shumohome
数模之家
13





后逐
渐变

逐渐变长,南半球直
杆反之。
中心,向东西两
侧影子逐渐变长
中心,向南北两
侧影子逐渐变长
误差分析与改进
由于第一问求解时所用的公式会在二、三、四问的建模过程中沿用,故于此
处,有必要对使用的模型进行误差分析。
于第一问的计算求解模型建立过程中,其核心公式为:太阳高度角的计算公
式,即:  arcsin sin sin  cos cos cos t 。各项中,太阳赤纬 以及时角t
的修正E 都是近似式,因而本模型的主要系统误差来自于此两项。由于 与E 皆
为仅仅与日期d 有关的变量,故可以选择一个精度较高的天文计算器[4],将之计
算的d 为1 至365 时的 和E 值作为真实值,而将本文模型中计算的 与E 值作
为计算值,做误差分析。
本文作为参考的计算器为MIDC SPA CALCULATOR[5]。对比计算得出的误
差如下图所示:
图 6.1 误差对比
如图 6.1 所示,上方图为E 值在取1 至365 时采用本文以及真实结果的比较,
下方图为太阳赤纬 的计算值和理论值的比较。其中,蓝色线为本文方案的计算
本资料由:数模之家(春哥带你学建模)收集并分享,需要更多数学建模资料添加微信公众号:shumohome
数模之家
14
值,棕色线为改进后结果,X 点为SPA 计算器算出的值,即真实值。可以看出,
本文的计算值与真实值之间存在一定偏差,因此,可以采用更为精确的模型加以
修正。
参考其它文献[6][7][8],得到更为精确的太阳赤纬 与时差E 的计算方式如下:
  
180sin 0.98565  2
arcsin 0.39779cos 0.98565 10 1.914
D
 D

   
       
 
 
E=720C C,其中:
   
 
tan( )
arctan
cos(23.44 )
180
1.914sin 2
10
360
365.24
B
A
C
B A W D
A W D
W
     
     
      


    

  

 


D为天数。当日期为1 月1 日时,D为0。
修正前后,本文模型计算得到的误差与SPA 计算器的比较如下表:
表 6.2 修正前后时差与太阳赤纬误差对比
最大残差绝对值 残差平方和 平均残差绝对值
修正前 修正后 修正前 修正后 修正前 修正后
时差 1.3206 0.2221 147.462 6.2765 0.4977 0.1153
太阳赤纬 1.1492 0.2221 115.0509 0.6689 0.4399 0.0377
可见,经过精度更高的公式修正之后,时差和太阳赤纬的误差得到了较为明
显的改善。然而,由于在之后的模型中,时差和太阳赤纬是用来计算太阳高度角
以及方位角的中间步骤,因此还须检测修正前后太阳高度角和方位角的误差变化。
表 6.3 修正前后太阳高度角与方位角误差对比
平均残差绝对值(单位为角度)
修正前 修正后
太阳高度角 0.3411 0.0222
方位角 0.3203 0.1760
从表 6.3 可以发现,修正后对于太阳高度角的误差减少效果较为明显。然而
可以看出,即使不加以修正,现在使用的公式产生的误差亦并非很大。
本资料由:数模之家(春哥带你学建模)收集并分享,需要更多数学建模资料添加微信公众号:shumohome
数模之家
15
图 6.4 修正前后太阳高度角与方位角残差分布图
图 6.4 为修正前后太阳高度角与方位角残差分布图。横轴以小时为单位,从
左到右为1 月1 日至12 月31 日,纵轴为残差。蓝色为修正前曲线,红色为修正
后曲线。可见,修正前的残差最大亦不足2 度,虽然修正后降低误差效果明显,
但是原本的误差也在可接受的范围之内。
考虑到修正的计算方法复杂,会影响模型效率,因而本文将在明确模型存在
系统误差的前提下,在之后几问沿用第一问中使用的天文学公式。
将修正后的公式代入第一问,可以得到影长在北京时间9 至15 点的最小的
值为3.6638 米,对应时间为11 点59 分24 秒。
问题(二)的求解
6.2.1 模型的分析
本问相对于第一问,不知道直杆的长度,但因知道太阳影子顶点在一个连续
时间段内的坐标数据,故相当于知道了此时间段内等间隔时刻影子的影长以及其
角度变化(即方向)。因此,合理的经纬度(直杆位置)应不仅仅能够拟合直杆
的影长及其变化,还应当拟合影子方位角的变化。
本资料由:数模之家(春哥带你学建模)收集并分享,需要更多数学建模资料添加微信公众号:shumohome
数模之家
16
由于特定时刻,特定经纬度下的影子太阳高度角和方位角皆可以计算,且同
一地点的同一直杆在一天中两不同时刻的影长只与太阳高度角有关,故可以计算
不同时刻的理论影长比及方位角差。因而,可以依照问题分析(二)中提出的思
路,构造以直杆所在经纬度为x, y时,计算所得两时刻的影子长度比与对应时刻
测得的影子长度比的差值最小为第一个目标函数,以计算所得两时刻方位角差值
与测得的方位角差值的差值最小为第二个目标的多目标规划模型,利用遍历求出
尽可能精确的解。
图 6.2.1.1 竖直直杆在不同时刻成影示意图
如图 6.2.1.1 所示,该图表现了竖直直杆在同一天中两个时刻的成影示意图。
该示意图表示的情况已过当地时间正午,故OB 为较早时的成影,OC 为较晚时
的成影。与第一问相同,B,C 分别为影子 OB、OC的太阳高度角。而第一
问中,参考有关太阳高度角的资料[1],可以得到太阳高度角的求解公式为
  arcsin sin sin  cos cos cos t ,且 0,
2


 
 
 
。因此,时刻1 2 t , t 的理论影
长之比即为

1 1 1

2 2 2

cot cot

cot cot
y
y
L L
L L
 
 
 ,实际影长之比可以通过所给数据计算,即为
1
2
y
y
L
L

由于方位角是在地平坐标系上的角[1],竖直直杆是垂直于大地的,因此可以
将面XOY 看作是地平坐标系,因而BOC即为时刻i, j方位角 1 2  , 的差值。查阅
资料可知,方位角的计算公式为:
sin cos cos sin cos 
arccos
cos
    t


  
  
 
,因而
可以计算出两时刻方位角理论计算值的差值的绝对值’ ’
1 2   。对于实际的方位
角差值,因已经知道了向量OB,OA的坐标,故 1 2 arccos
OB OA
OB OA
 

  。
通过计算所得到的理想经纬度应该能够使得实际和理论的影长之比以及方
本资料由:数模之家(春哥带你学建模)收集并分享,需要更多数学建模资料添加微信公众号:shumohome
数模之家
17
位角差值之差尽可能为0,因而,可建立多目标多目标规划模型,找出满足条件
的经纬度。
需要注意的是,因为附件中一共提供了 21 组数据,为了使一对数据的差别
尽可能大,从而减少误差,选择每一个时刻与其之后的第十个时刻为一对(例如,
时刻1,即14:42 的数据应与时刻11,即15:12 时的数据成为一对),得到共十组
数据,找出符合条件的经纬度使得每一对理论的影长之比的差的平方和接近于0,
且方位角差值之差的平方和接近于0,由此建立多目标规划模型。
6.2.2 模型的建立
根据模型分析及问题分析,建立已实际和理论的影长之比以及方位角差值之
差尽可能为0 二者为目标的多目标规划模型。目标函数为:
 
     
2

1 ’
,
2
’ ’
2
,
min ,
min ,
yi yi
i j yj yj
i j i j
i j
L L
f l
L L
f l

    
 
      
 
      
 


 和l 分别为直杆所在的纬度和经度。目标函数是两个关于经纬度的二元函
数。其中, 1 f 为每一对理论与实测的影长之比的差的平方和, , yi yj L L 分别为时刻
i, j测得的影长, ’ ’ , yi yj L L 为时刻i, j的理论影长。 2 f 为每一对理论与实测方位角差
值之差的平方和, ,i j   分别为时刻i, j 测得的方位角, ’ ’ ,i j   分别为时刻i, j 的理
论方位角值。
目标函数的约束条件如下:
本资料由:数模之家(春哥带你学建模)收集并分享,需要更多数学建模资料添加微信公众号:shumohome
数模之家
18
 
 
 
     
 



0

cot

cot
0,
2
arcsin sin sin cos cos cos
sin cos cos sin cos
arccos
cos
360
=23.45 sin 81
365
9.87sin 2 7.53cos 1.5sin
15 12
15 60
360
81
365
{1
yi i
yj j
k
k k
k
k
k
k k
L
L
t
t
d
l B B B
t t
B d
i




    
   



 
 
 
 
  
  
 
 
   
 
   
     
 
 

 
 
, 2,3, 4,5,6, 7,8,9,10}, 10
90 ,90
180 ,180
j i
l





















   
    

    
条件一给出了理论影长之比的计算方式, k  为k 时刻的太阳高度角, k 可取
i, j。条件二给出了 k  的约束条件,即k  须为锐角,仅仅在这种情况下,当地时
间才为白天,讨论影子才有意义。条件三给出了k  的计算方式: 为太阳赤纬,
 为竖直杆所在地的纬度,竖直杆在北半球时, 为正,否则为负。k t 为时刻k
的时角。条件四给出了理论方位角的计算公式。
条件五、六、七给出了 、k t 以及一个相关量B 的求解方式。其中d 代表天
数,当日期为 1 月 1 日时,d=1,本问中d 108;l为直杆所在经度, k 0 t 为时刻
k 的格林威治标准时间。
条件八给出了i, j对应所指的时刻。最后两个条件给出了经纬度的定义域。
综上所述,本文建立的多目标规划模型为:
 
     
2

1 ’
,
2
’ ’
2
,
min ,
min ,
yi yi
i j yj yj
i j i j
i j
L L
f l
L L
f l

    
 
      
 
      
 


本资料由:数模之家(春哥带你学建模)收集并分享,需要更多数学建模资料添加微信公众号:shumohome
数模之家
19
 
 
 
     
 



0

cot

cot
0,
2
arcsin sin sin cos cos cos
sin cos cos sin cos
arccos
cos
360
=23.45 sin 81
365
9.87sin 2 7.53cos 1.5sin
15 12
15 60
360
81
365
{1
yi i
yj j
k
k k
k
k
k
k k
L
L
t
t
d
l B B B
t t
B d
i




    
   



 
 
 
 
  
  
 
 
   
 
   
     
 
 

 
 
, 2,3, 4,5,6, 7,8,9,10}, 10
90 ,90
180 ,180
j i
l





















   
    

    
6.2.3 模型的求解
由于在所给的数据中,影子长度的变化较为明显,而角度的变化则相比之下
更为微小,因此即使是在所用公式精度相同的情况下,以角度为首要标准来解答
也会产生更大的误差。因此,为较为精确的定位,本文采用分层求解法解该多目
标规划问题。
以   1 min f , l(即影子长度比之差)为首要优化目标,在完成此步优化之后,
得到多组 , l 使得   1 f , l   ,进而以第一步的优化结果作为约束条件,即
  1 f , l  ,以   2 min f , l 作为优化目标再次优化,得出答案。
通过 MATLAB 编程遍历求解,得出最为符合条件的地点为北纬19 度,东经
109 度,位于我国海南省境内:
本资料由:数模之家(春哥带你学建模)收集并分享,需要更多数学建模资料添加微信公众号:shumohome
数模之家
20
图 6.2.3.1 北纬19 度,东经109 度定位(谷歌地图)
由于本文采用的计算方法自身存在一定的系统误差,因而此多目标规划得到
的答案并非精确的答案。于结果分析中,本文会进一步给出各经纬度下的误差分
析,指出其他可能点所在的范围。
6.2.4 结果分析
得到最优结果之后,可以分别以 f1 , l 值为因变量,以纬度、经度为自变量
作图分析。如本文之前所解释,角度的计算引起的偏差可能较大,因此结果分析
部分将着重于影子长度的理论值和测量值吻合情况,这一做法亦与本文求解多目
标规划模型的思路一致。所得图像如下:
图 6.2.4.1 不同经纬度的理论影长与实际测得影长的吻合情况
图 6.2.4.1 反映的是不同经纬度的理论影长与实际测得的影长的吻合情况,图
中,横轴为经度,经度为负即代表西经,为正即东经;纵轴为纬度,纬度为正时
即北纬,为负时即南纬。颜色越蓝,越深则代表此经纬度的直杆吻合情况越好;
颜色越黄,越浅即代表此经纬度的直杆吻合情况越差。
为了更好的描述   1 f , l 的相对大小,可对   1 f , l 的数据值做离差标准化,将
本资料由:数模之家(春哥带你学建模)收集并分享,需要更多数学建模资料添加微信公众号:shumohome
数模之家
21
离差标准化后的值命名为误差程度,该值越低代表误差越小。
根据该图可以看出,除了北纬19 度,东经109 度以外,还有其他的地点的
误差程度也相对较好。因此,可以再遍历所有的地点,搜寻各个局部的最优中心:
即以此点为中心,其邻近点的误差程度皆劣于该点。找到这些局部最优中心后,
再筛选出误差程度低于0.002 的点,可以得到七个中心,这些中心也是直杆可能
的所在地:
表 6.2.4.2 直杆可能所在地经纬度
纬度 经度 误差程度 大致地点
24 100 0.0012 云南
-3 104 0.0017 印度尼西亚,巨港
21 106 0.0008 越南,河内
20 108 0.0019 海南西北
19 109 0 海南
18 110 0.0006 海南东南
17 111 0.0018 西沙群岛
表 6.2.4.2 中,纬度为正代表北纬,反之代表南纬;经度为正代表东经,反
之为西经。可以看出,北纬19 度,东经109 度也包含在这7 个点中,其误差程
度亦是最小的。误差程度越小代表:依照本题模型,直杆处在该地的可能性越
大。可以看出:直杆处在海南(北纬19 度,东经109 度)的可能性是最大的,
其次是在海南东南的北纬18,东经109 度以及越南河内(北纬21 度,东经106
度)。但因模型自身或存在系统误差,测量影子坐标时亦可能出现错误,因此,
仍然不能排除直杆位于其他各地点的可能性。
以上各点在地图上的相对位置如下:
图 6.2.4.3 直杆可能所在地在地图上的位置(谷歌地图截图)
如图 6.2.4.3 所示,上述7 个可能的位置都已经用星星标出,其中,加上红色
本资料由:数模之家(春哥带你学建模)收集并分享,需要更多数学建模资料添加微信公众号:shumohome
数模之家
22
图钉的点为北纬 19 度,东经109 度。
问题(三)的求解
6.3.1 模型的分析
第三问的其他条件与第二问一致,只不过相比于第二问,日期是一个未知值。
因此,可以沿用第二问的多目标规划模型,不过须将原来的二元目标函数更改为
含有日期(几月几日)、纬度和经度三个未知数的三元函数,再建立与第二问类
似的多目标规划模型,采用分层求解法遍历求解。
6.3.2 模型的建立
根据模型分析,参照第二问的模型,建立已实际和理论的影长之比以及方位
角差值之差尽可能为0 二者为目标的多目标规划模型。目标函数为:
 
     
2

1 ’
,
2
’ ’
2
,
min , ,
min , ,
yi yi
i j yj yj
i j i j
i j
L L
f l d
L L
f l d

    
 
      
 
      
 


 和l 分别为直杆所在的纬度和经度,d 为日期,当日期为1 月1 日时,d 的
值为1。目标函数是两个关于经纬度与日期的二元函数。其中, 1 f 为每一对理论
的影长之比的差的平方和, , yi yj L L 分别为时刻i, j测得的影长, ’ ’ , yi yj L L 为时刻i, j
的理论影长。2 f 为方位角差值之差的平方和, ,i j   分别为时刻i, j测得的方位角,
’ ’ ,i j   分别为时刻i, j的理论方位角值。
目标函数的约束条件如下:
本资料由:数模之家(春哥带你学建模)收集并分享,需要更多数学建模资料添加微信公众号:shumohome
数模之家
23
 
 
 
     
 

1 1

2 2

0

cot

cot
0,
2
arcsin sin sin cos cos cos
sin cos cos sin cos
arccos
cos
360
=23.45 sin 81
365
9.87 sin 2 7.53cos 1.5sin
15 12
15 60
360
81
365
{1
y
y
k
k k
k
k
k
k k
L
L
t
t
d
l B B B
t t
B d
i




    
   



 
 
 
 
  
  
 
 
   
 
   
     
 
 

   
,2,3,4,5,6, 7,8,9,10}, 10
90 ,90 , 180 ,180
[1,365]
j i
l
d





















   
         

 
条件一给出了理论影长之比的计算方式, k  为k 时刻的太阳高度角, k 可取
i, j。条件二给出了 k  的约束条件,即k  须为锐角。条件三给出了k  的计算方式:
 为太阳赤纬,  为竖直杆所在地的纬度,竖直杆在北半球时, 为正,否则为
负。k t 为时刻k 的时角。条件四给出了理论方位角的计算公式。
条件五、六、七给出了 、k t 以及一个相关量B 的求解方式。其中d 代表天
数, l 为直杆所在经度, k 0 t 为时刻k 的格林威治标准时间。
条件八给出了i, j 对应所指的时刻。最后三个条件分别给出了经纬度和日期
的定义域。
综上所述,本文建立的多目标规划模型为:
 
     
2

1 ’
,
2
’ ’
2
,
min , ,
min , ,
yi yi
i j yj yj
i j i j
i j
L L
f l d
L L
f l d

    
 
      
 
      
 


本资料由:数模之家(春哥带你学建模)收集并分享,需要更多数学建模资料添加微信公众号:shumohome
数模之家
24
 
 
 
     

1 1

2 2

0

cot

cot
0,
2
arcsin sin sin cos cos cos
sin cos cos sin cos
arccos
cos
. . 360 =23.45 sin 81
365
9.87sin 2 7.53cos 1.5sin
15 12
15 60
360
8
365
y
y
k
k k
k
k
k
k k
L
L
t
t
s t d
l B B B
t t
B d




    
   



 
 
 
 
  
  
 
      
 
   
     
 
   
   
1
{1,2,3,4,5,6, 7,8,9,10}, 10
90 ,90 , 180 ,180
[1,365]
i j i
l
d





















    
         

 
6.3.3 模型的求解
与第二问类似,采用分层求解法,利用MATLAB 编程遍历求解,可得到根
据附件二、三计算出相应位置相应的日期为:
表、图6.3.3.1 第三问最可能答案及在地图上的大概位置(谷歌地图)
纬度 经度 日期 大致地点
附件二 40 79 5/25 新疆图木舒克市
附件三 33 106 10/31 陕西,汉中广元之间
纬度为正即为北纬,经度为正代表为东经。图中两星标点即为上述两点。其
本资料由:数模之家(春哥带你学建模)收集并分享,需要更多数学建模资料添加微信公众号:shumohome
数模之家
25
中,带图钉的点为附件二的位置。由于年份不可确定,只能确定日期,因而对应
的日期具体年份不可知。同时,因为本题的计算方式存在一定的误差,因而满足
目标函数最好的解并不一定是真实的解,进一步分析会在结果分析中给出。
6.3.4 结果分析
得到最优结果之后,与第二问类似,可以分别对特定的纬度与经度 , l 以其
对应的最小的  1 f , l,d 值为因变量,以纬度、经度为自变量作图分析。因附件二
数据和附件三的数据在分析过程中没有区别,故本部分仅仅以附件二为样本,进
行结果分析。所得图像如下:
图 6.3.4.1 不同经纬度的理论影长与实际测得影长的误差值
图 6.3.4.1 反映的是不同经纬度的理论影长与实际测得的影长的吻合情况,图
中,横轴为经度,经度为负即代表西经,为正即东经;纵轴为纬度,纬度为正时
即北纬,为负时即南纬。竖轴为误差值,该值越小,代表理论影长与实际影长的
误差越小,即目标吻合情况越好。
为了更好的描述  1 f , l,d 的相对大小,可对   1 f , l,d 的数据值做离差标准
化,更易于绘图。此外,为了便于观察,本文将误差值大于 0.001 的点的误差值
统一设为0.001,由此在保留观察误差值趋势变化的同时,便于读者观察。
根据该图可以看出,当遍历地点逐渐靠近到东经50 度,南北回归线范围前
后时,误差明显减小,由此从该图可以看出,较为准确的定位应该在此范围内。
为了更进一步的从图像上看出结果,可以将误差值大于0.000001 的点的误差值
全部规定为0.000001,修正后得到的图如下:
本资料由:数模之家(春哥带你学建模)收集并分享,需要更多数学建模资料添加微信公众号:shumohome
数模之家
26
图 6.3.4.2 修正后不同经纬度的理论影长与实际测得影长的误差值
从图 6.3.4.2 中可以明显的看出,存在两簇地理位置的点的误差很小,因本题
所采用的模型自身存在一定的误差,所以尽管各点对第一个目标的满足情况仍然
有差异,但是不能排除这些点为正确地点的可能。因而,图中的结果进一步体现
了选择若干的可能点作为可能所在的位置的理由,并展示了这些结果的地理位置
之间直观的相对关系。对附件二中的点取误差值低于10-7,附件三中的点取误差
值小于10-9,得出的可能解如下:
表 6.3.4.3 第三问可能答案
附件二 附件三




日期 误差值




日期 误差值
78 -41 12/3 3.22E-08 110 -36 4/14 7.23E-09
82 -41 1/7 6.56E-08 111 -34 8/17 1.96E-10
77 -40 11/24 3.94E-08 112 -33 8/2 7.32E-09
82 -40 1/17 8.91E-08 112 -32 7/29 1.19E-09
76 -39 11/17 3.59E-08 110 -31 5/31 9.22E-09
82 -38 1/30 8.85E-08 112 -31 7/21 2.81E-09
83 -38 1/29 8.93E-08 110 -30 6/3 8.39E-09
75 -37 11/6 6.08E-08 107 30 11/22 8.75E-09
82 -37 2/4 8.77E-08 108 30 12/1 8.58E-09
80 36 8/10 7.88E-08 109 30 12/10 9.94E-09
78 37 5/7 5.73E-08 113 30 1/20 6.35E-09
78 38 5/12 2.83E-08 107 31 11/20 6.88E-09
81 38 7/31 5.19E-08 113 31 1/16 5.93E-09
81 39 7/25 9.54E-08 106 32 11/4 4.67E-09
79 40 5/25 0 106 33 10/31 0
81 40 7/19 2.11E-08 106 34 10/25 1.19E-10
79 41 6/2 7.30E-08 114 34 2/10 5.94E-09
本资料由:数模之家(春哥带你学建模)收集并分享,需要更多数学建模资料添加微信公众号:shumohome
数模之家
27
80 41 6/4 6.93E-08 106 35 10/18 9.48E-09
81 41 7/10 2.97E-08
上述的点为第三问中附件二和附件三的直杆可能的所在地及对应数据记录
的日期,注意年份无法确定。误差值越小的可认为直杆在该位置的可能越大。
问题(四)的求解
6.4.1 模型的分析
在数据分析部分已经得出了视频所给的时间间隔中各个时刻的太阳高度角
信息。由于已知关系式  arcsin sin sin  cos cos cos t ,且各个时刻的太阳
高度角及时间皆为已知量,故可以组成若干个方程形成的方程组,建立计算求解
模型求解。
6.4.2 模型的建立
根据模型分析,参照第二、三问的模型,当日期已知时,建立的模型为:
arcsin sin sin cos cos cos  i i i i        t
日期未知时,建立的模型为:
arcsin sin sin cos cos cos  i i       t
其中,  
360
=23.45 sin 81
365
 d
 
   
 
,当日期未知时,d 是一个变量。d 为日期,
当日期为1 月1 日时,d 的值为1。 为直杆所在的纬度,为未知数。t 为时角,
其计算方式为:
     
0
9.87sin 2 7.53cos 1.5sin
15 12
15 60 k k
l B B B
t t
   
     
 
,其
中,  
360
81
365
B  d  ,l为直杆所在的经度,是一个未知数。此处的i代表该方程
为利用第i 时刻数据所列的方程。
因而可以看出,当日期已知时,待求解方程组中有 与l 两个未知数;日期未
知时,有 与l 和d 三个未知数。
综上所述,所建立的计算求解模型为:
当日期已知时:
arcsin sin sin cos cos cos  i i i i        t
本资料由:数模之家(春哥带你学建模)收集并分享,需要更多数学建模资料添加微信公众号:shumohome
数模之家
28
日期未知时:
arcsin sin sin cos cos cos  i i       t
其中:
 
     
 
0
360
=23.45 sin 81
365
9.87sin 2 7.53cos 1.5sin
15 12
15 60
360
81
365
k k
d
l B B B
t t
B d

  
    
  
    
      
  

  

6.4.3 模型的求解
利用 MATLAB 编程,考虑到天文学公式自身存在系统误差,故转换为非线
性规划模型求解。
在日期已知的情况下,以   2
min sin sin cos cos cos sin i i i i
i
      t   为
目标函数, 在日期未知的情况下, 以
  2
min sin sin cos cos cos sin i i
i
      t   为目标函数,遍历求解,可得到日
期确定与未确定的情况下,计算出的相应经纬度、日期和大致地点如下:
表、图6.4.3.1 第四问最可能答案及在地图上的大概位置(谷歌地图)
纬度 经度 日期 大致地点
日期确定 43 115 7/13 内蒙古锡林郭勒盟西南
日期未定 44 113 6/25 内蒙古锡林郭勒盟西
纬度为正即为北纬,经度为正代表为东经。图中两星标点即为上述两点。其
中,带图钉的点为日期未定时的位置。年份不可确定。虽然本文在日期已知和未
本资料由:数模之家(春哥带你学建模)收集并分享,需要更多数学建模资料添加微信公众号:shumohome
数模之家
29
知的情况下都给出了答案,但是因为本题的计算方式存在一定的误差,因而满足
目标函数最好的解并不一定是真实的解,进一步分析会在结果分析中给出。
6.4.4 结果分析
得到结果之后,与第三问类似,可以分别对特定的纬度与经度 , l,以其对应
的最小的离差标准化后的本问目标函数值为因变量,以纬度、经度为自变量作图
分析。因当日期确定的情况下该部分的分析与日期未确定时类似,故本部分主要
以日期未确定时的情况为样本,进行结果分析。所得图像如下:
图 6.4.4.1 不同经纬度的理论影长与实际测得影长的误差值
图 6.4.4.1 反映的是不同经纬度的理论影长与实际测得的影长的吻合情况,图
中,横轴为经度,经度为负即代表西经,为正即东经;纵轴为纬度,纬度为正时
即北纬,为负时即南纬。竖轴为误差值,该值越小,代表理论影长与实际影长的
误差越小,即目标吻合情况越好。
为了更好的描述目标函数值的相对大小,可对目标函数值做离差标准化,更易
于绘图。
本资料由:数模之家(春哥带你学建模)收集并分享,需要更多数学建模资料添加微信公众号:shumohome
数模之家
30
根据该图可以看出,当地点在东经110 度左右,西经150 度左右,误差明显
减小,由此从该图可以看出,较为准确的定位应该在此范围内。为了更进一步的
从图像上看出结果,可以将误差值大于0.001 的点的误差值全部规定为0.001,
修正后得到的图如下:
图 6.4.4.2 修正后不同经纬度的理论影长与实际测得影长的误差值
从图 6.4.4.2 中可以明显的看出,存在两簇地理位置的点的误差很小,其位置
分别在东经110 度左右,北纬40 度、南纬40 度左右。因本题所采用的模型自身
存在一定的误差,所以尽管各点对第一个目标的满足情况仍然有差异,但是不能
排除这些点正确的可能。因而,图中的结果进一步体现了选择若干的可能点作为
可能所在的位置的理由,并展示了这些结果的地理位置之间直观的相对关系。
对于日期确定的情况,可作类似的图如下:
图 6.4.4.3 日期确定时修正后不同经纬度的理论影长与实际测得影长的误差值
图 6.4.4.3 给出了日期确定的情况下不同经纬度的理论影长与实际测得的影
长的吻合情况。横轴为经度,纵轴为纬度,等高线状图代表着离差标准化后误差
本资料由:数模之家(春哥带你学建模)收集并分享,需要更多数学建模资料添加微信公众号:shumohome
数模之家
31
值相等的点。越蓝代表误差越小。同样的,可以看出在东经110 度左右,北纬40
度、南纬10 度左右存在两簇地理位置的点的误差很小。
对日期不定的情况下取误差值低于10-3 的点,附件三中的点取误差值小于
0.05 的点,得出的可能解如下:
表 6.3.4.3 第四问可能答案
日期确定 日期未定
经度 纬度 误差值 经度 纬度 日期 误差值
115 43 0 111 -44 12/5 0.000904
116 46 0.009332 112 -44 12/21 0.000711
114 39 0.016437 113 -44 12/28 0.000268
127 -6 0.017337 111 -43 11/26 3.73E-05
128 -7 0.022651 114 -43 1/4 0.000304
126 -5 0.027358 115 -43 1/9 0.000495
117 48 0.028334 120 -42 1/30 0.000835
133 -12 0.029911 121 -41 2/4 0.000927
129 -8 0.030939 118 40 5/1 0.000947
132 -11 0.032562 115 42 5/15 0.000955
130 -9 0.035856 116 42 5/11 0.000939
131 -10 0.036016 117 42 7/26 0.00046
134 -13 0.036652 113 43 5/30 0.000466
118 50 0.047725 114 43 7/7 0.000311
124 -2 0.04793 113 44 6/25 0
114 44 7/5 0.00074
上述的点为第四问中日期确认时与日期不确认时直杆可能的所在地及对应
数据记录的日期,注意年份无法确定。误差值越小的可认为直杆在该位置的可能
越大。
7. 模型的评价、改进及推广
7.1 模型评价
7.1.1 模型优点
本文中的模型具有以下优点:
1. 模型结构简单,在精度要求不高的情况下可以采纳。
2. 求解时可采用遍历求解,在现有情况下提供了较为精确的答案。
3. 采用多目标规划模型,并合理的根据模型采用的计算公式和实际情况决定
目标的优先顺序,考虑更为周全。
本资料由:数模之家(春哥带你学建模)收集并分享,需要更多数学建模资料添加微信公众号:shumohome
数模之家
32
4. 在提供现有模型给出的最可能解的前提下,模型可给出其他可能解,方便
读者参考选择。
7.1.2 模型缺点
本文中的模型主要缺点如下:
1. 没有考虑诸如海拔等复杂的现实因素。
2. 模型本身存在一定的误差,无法给出十分准确的答案。
3. 第四问模型缺少通用性,因没有对视频进行透视矫正,所用方法并不适用
于一般的视频。
7.1.3 模型改进
1. 在条件允许的情况下,多考虑如海拔等实际因素。
2. 利用误差分析的结果对现有的模型做出修正,提高精准度。
3. 在第四问模型中加以透视矫正,增强模型的通用性。
7.1.4 模型推广
1) 用于分析同一组照片拍摄的地点及时间。
2) 检验视频或照片是否被为伪造品。
8. 参考文献
[1] PV Education, Azimuth Angle, http://www.pveducation.org/pvcdrom/propertiesof-
sunlight/azimuth-angle, 2015/9/12
[2] PV Education, Elevation Angle, http://www.pveducation.org/pvcdrom/propertiesof-
sunlight/elevation-angle, 2015/9/12
[3] PV Education, The Sun’s Position,
http://www.pveducation.org/pvcdrom/properties-of-sunlight/suns-position , 2015/9/12
[4] Reda I, Andreas A. Solar Position Algorithm for Solar Radiation Applications. 2003.
[5] National Renewable Energy Laboratory, MIDC SPA Calculator,
http://www.nrel.gov/midc/solpos/spa.html, 2015/9/13
[6] Wikipedia, Position of the Sun. (2015, August 29),
https://en.wikipedia.org/w/index.php?title=Position_of_the_Sun&oldid=678477646,
2015/9/13
[7] C Johnson. Sundial Time Correction - Equation of Time, http://mbsoft.
com/public3/equatime.html, September 13, 2015
[8] Wikipedia, Equation of time. (2015, May 27),
https://en.wikipedia.org/w/index.php?title=Equation_of_time&oldid=664315122,
September 13, 2015
本资料由:数模之家(春哥带你学建模)收集并分享,需要更多数学建模资料添加微信公众号:shumohome
数模之家
33
附录
附录一:所用程序
1. MATLAB R2015a
2. Adobe After Effect CC
3. Adobe Photoshop CS
附录二:MATLAB 源程序代码
第一问解答代码(Q1):
clear all; clc;
L = 3;
LT = 9:0.001:15;
GMT=LT-8;
DateString = ‘22-Oct-2015’;
DateString2 = ‘1-Jan-2015’;
formatIn = ‘dd-mmm-yyyy’;
d= datenum(DateString,formatIn) -
datenum(DateString2,formatIn) + 1 ;
LONGITUDE=116.3914;
LATITUDE=39.0972;
dGMT = LT - GMT;
LSTM = 15.*dGMT;
B = (360./365).*(d - 81);
EoT = 9.87.*sin(2.*B./180.*pi) -
7.53.*cos(B./180.*pi) - 1.5.*sin(B./180.*pi);
TC = 4.*(LONGITUDE - LSTM) + EoT;
LST = LT + TC./60;
HRA = 15.*(LST - 12);
DEC = 23.45.*sin(B./180.*pi);
elevation =
asin(sin(DEC./180.*pi).*sin(LATITUDE./180.*pi) +
cos(DEC./180.*pi).*cos(LATITUDE./180.*pi).*cos(HR
A./180.*pi));
azimuth =
acos((sin(DEC./180.*pi).*cos(LATITUDE./180.*pi) -
本资料由:数模之家(春哥带你学建模)收集并分享,需要更多数学建模资料添加微信公众号:shumohome
数模之家
34
cos(DEC./180.*pi).*sin(LATITUDE./180.*pi).*cos(HR
A./180.*pi))./cos(elevation));
%azimuth = azimuth(elevation>0);
%elevation = elevation(elevation>0);
for i=1:length(azimuth)
if LST(i)>12 || HRA(i)>0
azimuth(i) = 2*pi - azimuth(i);
end
end
%disp(elevation./pi.*180);
%disp(azimuth./pi.*180);
%polar(azimuth,L./tan(elevation),’-‘);
plot(LT,L./tan(elevation),’-‘);
xlabel(‘±±¾©Ê±¼ä’);
ylabel(‘Ö±¸ËÓ°³¤’);
title(‘Ö±¸ËμÄÌ«ÑôÓ°×Ó³¤¶Èμı仯ÇúÏß’);
第二问解答程序(Q2_tol)
data = xlsread(‘dataQ2’ , ‘sheet1’ );
len = sqrt(data(: , 1) .^2 + data(: , 2) .^ 2);
left = 1 : ceil(length(len)/2);
right = ceil(length(len)/2) : (length(len) );
%left = 1 : length(len)-1;
%right = 2 : length(len);
bleft = left;
bright = right;
rateT = len(left) ./ len(right) ;
Xl = data(bleft , 1);
Xr = data(bright , 1);
Yl = data(bleft , 2);
Yr = data(bright , 2);
angleT =acos( (Xl .* Xr + Yl .* Yr) ./
( len(bleft) .* len(bright) ) );
longRange = - 180: 180;
latiRange = -90 : 90 ;
ErrorT = zeros(length(longRange) ,
length(latiRange));
本资料由:数模之家(春哥带你学建模)收集并分享,需要更多数学建模资料添加微信公众号:shumohome
数模之家
35
ErrorTA = zeros(length(longRange) ,
length(latiRange));
DateString = ‘11-Sep-2015’;
DateString2 = ‘1-Jan-2015’;
formatIn = ‘dd-mmm-yyyy’;
d= datenum(DateString,formatIn) -
datenum(DateString2,formatIn) + 1 ;
L = 3;
LT = 14.7 :0.05:15.7;
%LT = T;
GMT=LT- 8;
dGMT = LT - GMT;
LSTM = 15.*dGMT;
B = (360./365).*(d - 81);
EoT = 9.87.*sin(2.*B./180.*pi) -
7.53.*cos(B./180.*pi) - 1.5.*sin(B./180.*pi);
for i =1 : length(longRange)
%disp(i);
for j =1 : length( latiRange)
LONGITUDE =longRange(i);
LATITUDE = latiRange(j);
TC = 4.*(LONGITUDE - LSTM) + EoT;
LST = LT + TC./60;
HRA = 15.*(LST - 12);
DEC = 23.45.*sin(B./180.*pi);
elevation =
asin(sin(DEC./180.*pi).*sin(LATITUDE./180.*pi) +
cos(DEC./180.*pi).*cos(LATITUDE./180.*pi).*cos(HR
A./180.*pi));
azimuth =acos
(sin(DEC./180.*pi).*cos(LATITUDE./180.*pi) -
cos(DEC./180.*pi).*sin(LATITUDE./180.*pi).*cos(HR
A./180.*pi))./cos(elevation);
azimuth(elevation < 0) = 0;
angleP = abs(azimuth(bleft) -
azimuth(bright));
%azimuth = azimuth(elevation>0);
%elevation = elevation(elevation>0);
本资料由:数模之家(春哥带你学建模)收集并分享,需要更多数学建模资料添加微信公众号:shumohome
数模之家
36
%azimuth(LST > 12) = 2*pi - azimuth(LST >
12);
lenP = L ./tan(elevation);
lenP(elevation < 0) = 0;
rateP =lenP(left) ./ lenP(right ) ;
ErrorT(i , j) = sqrt(sum(( rateTrateP’
) .^ 2));
ErrorTA(i , j) = sqrt(sum((anglePangleT’
) .^ 2));
end
end
Tmin = min(min(ErrorT));
tol = 1.01*Tmin;
%ErrorT((ErrorT > tol)) = tol;
ErrorT((ErrorT > 1)) = 1;
ErrorTA((ErrorTA > 1)) = 1 ;
ErrorTA((ErrorT == tol)) = 1 ;
%ErrorT = ( ErrorT - min(min(ErrorT))) /
(max(max(ErrorT)) - min(min(ErrorT)));
%ErrorTA = ( ErrorTA - min(min(ErrorTA))) /
(max(max(ErrorTA)) - min(min(ErrorTA)));
%ErrorTT = ErrorT + ErrorTA;
[X , Y] = meshgrid(longRange , latiRange);
%surf(X , Y , ErrorTT’);
figure;
contour(X , Y , ErrorT(selectX , selectY)’);
%figure;
%surf(X , Y , ErrorT’);
%figure;
%contour(X , Y , ErrorTA’);
%[a , b] = find(ErrorTT ==min(min(ErrorTT)));
第三问解答程序(Q3)
clc;
clear all;
data = xlsread(‘dataQ2’ , ‘sheet2’ );
xTemp = data(: , 1);
yTemp = data(: , 2);
T = (12 + 41/60 ): 1/20 : (13 + 41/60 );
len = sqrt(xTemp .^ 2 + yTemp .^ 2);
本资料由:数模之家(春哥带你学建模)收集并分享,需要更多数学建模资料添加微信公众号:shumohome
数模之家
37
%rateT = len(1 : ceil(length(len)/2)) ./
len(ceil(length(len)/2) : length(len) ) ;
rateT = len(1 : ( length(len) - 1)) ./ len(2 :
length(len) ) ;
%angleT =acos( (data(1 : ceil(length(len)/2) ,
1) .* data(ceil(length(len)/2) :
(length(len) ) , 1) + …
% data(1 : ceil(length(len)/2) , 2) .*
data(ceil(length(len)/2) : (length(len) ) ,
2) )./ …
%( len(1 : ceil(length(len)/2)) .*
len(ceil(length(len)/2) :
%length(len) )) )
%angleT = acos((data(1 : length(len) - 1 , 1) .*
data(2 : length(len) , 1) + …
%(data(1 : length(len) - 1 , 2) .* data(2:
(length(len) ) , 2))) ./…
%(len(1 : length(len) - 1) .* len(2 :
length(len) )) );
angleT = acos((xTemp(1 : length(len) - 1 ) .*
xTemp(2 : length(len) ) + …
(yTemp(1 : length(len) - 1 ) .* yTemp(2:
(length(len) ) ))) ./…
(len(1 : length(len) - 1) .* len(2 :
length(len) )) );
%angleT =atan(yTemp (2 : length(len) ) ./ xTemp
(2 : length(len) ) ) - …
% atan(yTemp (1 : length(len) - 1) ./ xTemp
(1 : length(len) - 1) );
longRange = -180 : 180;
latiRange = -90: 90 ;
dateRange = 1 : 365 ;
ErrorT = zeros(length(longRange) ,
length(latiRange) , length(dateRange));
%ErrorTA = zeros(length(longRange) ,
length(latiRange));
for i =1 : length(longRange)
disp(i);
for j =1 : length( latiRange)
for k = 1 : length( dateRange)
本资料由:数模之家(春哥带你学建模)收集并分享,需要更多数学建模资料添加微信公众号:shumohome
数模之家
38
d = k ;
LONGITUDE =longRange(i);
LATITUDE = latiRange(j);
L = 3;
LT = T;
GMT=LT- 8;
dGMT = LT - GMT;
LSTM = 15.*dGMT;
B = (360./365).*(d - 81);
EoT = 9.87.*sin(2.*B./180.*pi) -
7.53.*cos(B./180.*pi) - 1.5.*sin(B./180.*pi);
TC = 4.*(LONGITUDE - LSTM) + EoT;
LST = LT + TC./60;
HRA = 15.*(LST - 12);
DEC = 23.45.*sin(B./180.*pi);
elevation =
asin(sin(DEC./180.*pi).*sin(LATITUDE./180.*pi) +
cos(DEC./180.*pi).*cos(LATITUDE./180.*pi).*cos(HR
A./180.*pi));
%azimuth =acos
(sin(DEC./180.*pi).*cos(LATITUDE./180.*pi) -
cos(DEC./180.*pi).*sin(LATITUDE./180.*pi).*cos(HR
A./180.*pi))./cos(elevation);
%azimuth(elevation < 0) = 0;
% angleP = abs(azimuth(1 :
ceil(length(len)/2)) -
azimuth(ceil(length(len)/2) : length(len)));
% angleP = abs((azimuth(1 : length(len) -
1 ) - azimuth(2 : length(len))));
%azimuth = azimuth(elevation>0);
%elevation = elevation(elevation>0);
% azimuth(LST > 12) = 2*pi - azimuth(LST >
12);
lenP = L ./tan(elevation);
lenP(elevation < 0) = 0;
% rateP =lenP(1 : ceil(length(len)/2)) ./
lenP(ceil(length(len)/2) : length(len) ) ;
rateP =lenP(1 : length(len) - 1) ./
lenP(2: length(len) );
ErrorT(i , j , k) = (sum(( rateTrateP’
) .^ 2));
% ErrorTA(i , j) = sqrt(sum((angleP-
本资料由:数模之家(春哥带你学建模)收集并分享,需要更多数学建模资料添加微信公众号:shumohome
数模之家
39
angleT’ ) .^ 2));
end
end
end
ErrorT((ErrorT > 1)) = 1 ;
%ErrorTA((ErrorTA >0.5)) = 0.5 ;
%ErrorT = ( ErrorT - min(min(ErrorT))) /
(max(max(ErrorT)) - min(min(ErrorT)));
%ErrorTA = ( ErrorTA - min(min(ErrorTA))) /
(max(max(ErrorTA)) - min(min(ErrorTA)));
%ErrorTT = ErrorT + ErrorTA;
RangeX = 1 : 361 ;
RangeY = 1 : 181;
[X , Y] = meshgrid(longRange(RangeX) ,
latiRange(RangeY));
min(ErrorT() , [] , 3);
[temp , test]= min(ErrorT() , [] , 3);
temp= ( temp - min(min(temp))) / (max(max(temp))
- min(min(temp)));
%surf(X , Y , temp’);
A = temp(RangeX , RangeY);
A(A > 0.000001) = 0.000001;
%contour(X , Y , A’)
%surf(X , Y , ErrorTT’);
%surf(X , Y , A’);
%surf(X , Y , ErrorTA’);
acc = 1 * 10 ^( - 6);
[a , b] = find(A < acc);
a = a - 181 ;
b = b - 91;
c = [a b test( find(A < acc)) A(find(A <
acc)) ];
第四问解答程序(Q4,日期已知)
clc ;
clear all;
dataraw = xlsread(‘keyframeData’ , ‘sheet1’ );
top = [891.2 205] ;
bottom = [891.2 885];
L = bottom(2) - top(2) ;
本资料由:数模之家(春哥带你学建模)收集并分享,需要更多数学建模资料添加微信公众号:shumohome
数模之家
40
len =sqrt((dataraw(: , 1) - top(1)) .^2 +
(dataraw(: , 2) - top(2)) .^2);
startTime= 8 + 54/60 + 7/3600;
endTime = 9 + 34/60 + 46/3600;
step = (endTime - startTime) / (length(dataraw) -
1);
elevationTest = L ./ len ;
T = startTime : step : endTime;
longRange = -180: 180;
latiRange = -90 : 90 ;
ErrorT = zeros(length(longRange) ,
length(latiRange));
%ErrorTA = zeros(length(longRange) ,
length(latiRange));
DateString = ‘13-July-2015’;
DateString2 = ‘1-Jan-2015’;
formatIn = ‘dd-mmm-yyyy’;
d= datenum(DateString,formatIn) -
datenum(DateString2,formatIn) + 1 ;
for i =1 : length(longRange)
disp(i);
for j =1 : length( latiRange)
LONGITUDE =longRange(i);
LATITUDE = latiRange(j);
LT = T;
GMT=LT- 8;
dGMT = LT - GMT;
LSTM = 15.*dGMT;
B = (360./365).*(d - 81);
EoT = 9.87.*sin(2.*B./180.*pi) -
7.53.*cos(B./180.*pi) - 1.5.*sin(B./180.*pi);
TC = 4.*(LONGITUDE - LSTM) + EoT;
LST = LT + TC./60;
HRA = 15.*(LST - 12);
DEC = 23.45.*sin(B./180.*pi);
elevation =
(sin(DEC./180.*pi).*sin(LATITUDE./180.*pi) +
cos(DEC./180.*pi).*cos(LATITUDE./180.*pi).*cos(HR
A./180.*pi));
%azimuth =acos
(sin(DEC./180.*pi).*cos(LATITUDE./180.*pi) -
cos(DEC./180.*pi).*sin(LATITUDE./180.*pi).*cos(HR
本资料由:数模之家(春哥带你学建模)收集并分享,需要更多数学建模资料添加微信公众号:shumohome
数模之家
41
A./180.*pi))./cos(elevation);
%azimuth(elevation < 0) = 0;
%angleP = abs(azimuth(1 :
ceil(length(len)/2)) -
azimuth(ceil(length(len)/2) : length(len)));
%angleP = abs(azimuth(1 : length(len) - 1 )
- azimuth(2 : length(len)));
%azimuth = azimuth(elevation>0);
%elevation = elevation(elevation>0);
% azimuth(LST > 12) = 2*pi - azimuth(LST >
12);
% lenP = L ./tan(elevation);
% lenP(elevation < 0) = 0;
% rateP =lenP(1 : ceil(length(len)/2)) ./
lenP(ceil(length(len)/2) : length(len) ) ;
%rateP =lenP(1 : length(len) - 1) ./
lenP(2: length(len) );
ErrorT(i , j) = sqrt(sum(( elevationelevationTest’
) .^ 2));
%ErrorTA(i , j) = sqrt(sum((anglePangleT’
) .^ 2));
end
end
ErrorT((ErrorT > 1)) = 1 ;
%ErrorTA((ErrorTA >0.5)) = 0.5 ;
ErrorT = ( ErrorT - min(min(ErrorT))) /
(max(max(ErrorT)) - min(min(ErrorT)));
%ErrorTA = ( ErrorTA - min(min(ErrorTA))) /
(max(max(ErrorTA)) - min(min(ErrorTA)));
%ErrorTT = ErrorT + ErrorTA;
[X , Y] = meshgrid(longRange , latiRange);
%surf(X , Y , ErrorTT’);
%surf(X , Y , ErrorT’);
%surf(X , Y , ErrorTA’);
count = 0 ;
for i =2 : length(longRange) - 1
for j =2 : length( latiRange) - 1
if( ErrorT(i , j) - ErrorT(i , j + 1) < 0
&&ErrorT(i , j) - ErrorT(i , j - 1) < 0 &&
ErrorT(i , j) - ErrorT(i + 1 , j ) < 0 …
&& ErrorT(i , j) - ErrorT(i-1 ,
本资料由:数模之家(春哥带你学建模)收集并分享,需要更多数学建模资料添加微信公众号:shumohome
数模之家
42
j ) < 0 )
count = count + 1 ;
c(: , count ) = [i-181 j-91 ErrorT(i ,
j)];
end
end
end
contour(X ,Y , ErrorT’);
%[a , b] = find(ErrorTT ==min(min(ErrorTT)));
第四问解答程序(Q4_2,日期未知)
clc ;
clear all;
dataraw = xlsread(‘keyframeData’ , ‘sheet1’ );
top = [891.2 205] ;
bottom = [891.2 885];
L = bottom(2) - top(2) ;
len =sqrt((dataraw(: , 1) - top(1)) .^2 +
(dataraw(: , 2) - top(2)) .^2);
startTime= 8 + 54/60 + 7/3600;
endTime = 9 + 34/60 + 46/3600;
step = (endTime - startTime) / (length(dataraw) -
1);
elevationTest = L ./ len ;
T = startTime : step : endTime;
longRange = -180: 180;
latiRange = -90 : 90 ;
dateRange = 1 : 365;
ErrorT = zeros(length(longRange) ,
length(latiRange) , length(dateRange));
%ErrorTA = zeros(length(longRange) ,
length(latiRange));
for i =1 : length(longRange)
disp(i);
for j =1 : length( latiRange)
for k = 1 : length( dateRange)
d = k ;
LONGITUDE =longRange(i);
LATITUDE = latiRange(j);
LT = T;
GMT=LT- 8;
本资料由:数模之家(春哥带你学建模)收集并分享,需要更多数学建模资料添加微信公众号:shumohome
数模之家
43
dGMT = LT - GMT;
LSTM = 15.*dGMT;
B = (360./365).*(d - 81);
EoT = 9.87.*sin(2.*B./180.*pi) -
7.53.*cos(B./180.*pi) - 1.5.*sin(B./180.*pi);
TC = 4.*(LONGITUDE - LSTM) + EoT;
LST = LT + TC./60;
HRA = 15.*(LST - 12);
DEC = 23.45.*sin(B./180.*pi);
elevation =
(sin(DEC./180.*pi).*sin(LATITUDE./180.*pi) +
cos(DEC./180.*pi).*cos(LATITUDE./180.*pi).*cos(HR
A./180.*pi));
%azimuth =acos
(sin(DEC./180.*pi).*cos(LATITUDE./180.*pi) -
cos(DEC./180.*pi).*sin(LATITUDE./180.*pi).*cos(HR
A./180.*pi))./cos(elevation);
%azimuth(elevation < 0) = 0;
%angleP = abs(azimuth(1 :
ceil(length(len)/2)) -
azimuth(ceil(length(len)/2) : length(len)));
%angleP = abs(azimuth(1 : length(len) - 1 )
- azimuth(2 : length(len)));
%azimuth = azimuth(elevation>0);
%elevation = elevation(elevation>0);
% azimuth(LST > 12) = 2*pi - azimuth(LST >
12);
% lenP = L ./tan(elevation);
% lenP(elevation < 0) = 0;
% rateP =lenP(1 : ceil(length(len)/2)) ./
lenP(ceil(length(len)/2) : length(len) ) ;
%rateP =lenP(1 : length(len) - 1) ./
lenP(2: length(len) );
ErrorT(i , j , k) = sqrt(sum(( elevationelevationTest’
) .^ 2));
%ErrorTA(i , j) = sqrt(sum((anglePangleT’
) .^ 2));
end
end
end
ErrorT((ErrorT > 30)) = 30 ;
本资料由:数模之家(春哥带你学建模)收集并分享,需要更多数学建模资料添加微信公众号:shumohome
数模之家
44
RangeX = 1 : 361 ;
RangeY = 1 : 181;
[X , Y] = meshgrid(longRange(RangeX) ,
latiRange(RangeY));
[temp , test]= min(ErrorT() , [] , 3);
temp= ( temp - min(min(temp))) / (max(max(temp))
- min(min(temp)));
A = temp(RangeX , RangeY);
A(A > 0.1) =0.1;
surf(X , Y , A’);
acc = 1 * 10 ^( - 6);
[a , b] = find(A < acc);
a = a - 181 ;
b = b - 91;
c = [a b test( find(A < acc)) A(find(A <
acc)) ];
太阳赤纬误差分析(Declination_acc)
figure;
d = 0:364;
DEC2 = -asin( 0.39779 * cosd(0.98565*(d+10) +
1.914 * sind(0.98565*(d-2)) ))/pi*180;
B = (360./365).*(d - 82);
DEC1 = 23.45.*sind(B);
%B = (360./365).*(d + 10);
%DEC1 = -23.44.*cosd(B);
spa = dlmread(‘Data/SPA_DEC_EoT_2015/data.txt’);
spa = spa(:,1);
plot(0:5:364,spa(1:5:end),’kx’,d,DEC1,d,DEC2);
legend(‘SPA Data’,’Method 1’,’Method 2’);
disp(max(abs(DEC1’-spa)));
disp(max(abs(DEC2’-spa)));
disp(sum((DEC1’-spa).^2));
disp(sum((DEC2’-spa).^2));
本资料由:数模之家(春哥带你学建模)收集并分享,需要更多数学建模资料添加微信公众号:shumohome
数模之家
45
disp(mean(abs(DEC1’-spa)));
disp(mean(abs(DEC2’-spa)));
时差误差分析(EoT_acc)
D = 0:364;
W = 360/365.24;
A = W.*(D+10);
B = A+1.914.sind(W.(D-2));
C = (A-atand(tand(B)./cosd(23.44)))/180;
EoT2 = 720.*(C-round(C));
B = (360./365).*(D - 80);
EoT1 = 9.87.*sin(2.*B./180.*pi) -
7.53.*cos(B./180.*pi) - 1.5.*sin(B./180.*pi);
spa = dlmread(‘Data/SPA_DEC_EoT_2015/data.txt’);
spa = spa(:,2);
plot(0:5:364,spa(1:5:end),’kx’,d,EoT1,d,EoT2);
legend(‘SPA Data’,’Method 1’,’Method 2’);
disp(max(abs(EoT1’-spa)));
disp(max(abs(EoT2’-spa)));
disp(sum((EoT1’-spa).^2));
disp(sum((EoT2’-spa).^2));
disp(mean(abs(EoT1’-spa)));
disp(mean(abs(EoT2’-spa)));
总误差分析(Q1_err)
clear all;
LT = 0:23;
GMT=LT-8;
d=1:365;
LONGITUDE=116.3914;
LATITUDE=39.0972;
%convert to matrix
d = repmat(d’,1,24);
本资料由:数模之家(春哥带你学建模)收集并分享,需要更多数学建模资料添加微信公众号:shumohome
数模之家
46
LT = repmat(LT,365,1);
dGMT = 8;
LSTM = 15.*dGMT;
B = (360./365).*(d - 81);
EoT = 9.87.*sin(2.*B./180.*pi) -
7.53.*cos(B./180.*pi) - 1.5.*sin(B./180.*pi);
TC = 4.*(LONGITUDE - LSTM) + EoT;
LST = LT + TC./60;
HRA = 15.*(LST - 12);
DEC = 23.45.*sin(B./180.*pi);
D = d-1;
W = 360/365.24;
A = W.*(D+10);
B = A+1.914.sind(W.(D-2));
C = (A-atand(tand(B)./cosd(23.44)))/180;
EoT2 = 720.*(C-round(C));
TC2 = 4.*(LONGITUDE - LSTM) + EoT;
LST2 = LT + TC./60;
HRA2 = 15.*(LST - 12);
%DEC2 = -23.44.*cosd(A);
DEC2 = -asin( 0.39779 * cosd(0.98565*(D+10) +
1.914 * sind(0.98565*(D-2)) ))/pi*180;
elevation =
asin(sin(DEC./180.*pi).*sin(LATITUDE./180.*pi) +
cos(DEC./180.*pi).*cos(LATITUDE./180.*pi).*cos(HR
A./180.*pi));
azimuth =
acos((sin(DEC./180.*pi).*cos(LATITUDE./180.*pi) -
cos(DEC./180.*pi).*sin(LATITUDE./180.*pi).*cos(HR
A./180.*pi))./cos(elevation));
elevation2 =
asin(sin(DEC2./180.*pi).*sin(LATITUDE./180.*pi) +
cos(DEC2./180.*pi).*cos(LATITUDE./180.*pi).*cos(H
RA2./180.*pi));
azimuth2 =
acos((sin(DEC2./180.*pi).*cos(LATITUDE./180.*pi)
本资料由:数模之家(春哥带你学建模)收集并分享,需要更多数学建模资料添加微信公众号:shumohome
数模之家

47

cos(DEC2./180.*pi).*sin(LATITUDE./180.*pi).*cos(H
RA2./180.*pi))./cos(elevation2));
azimuth(LST>12) = 2*pi - azimuth(LST>12);
azimuth(LST<0) = 2*pi - azimuth(LST<0);
azimuth2(LST>12) = 2*pi - azimuth2(LST>12);
azimuth2(LST<0) = 2*pi - azimuth2(LST<0);
elevation = 90-elevation./pi.*180;
elevation =
reshape(elevation’,numel(elevation),1);
azimuth = (azimuth)./pi.*180;
azimuth = reshape(azimuth’,numel(azimuth),1);
elevation2 = 90-elevation2./pi.*180;
elevation2 =
reshape(elevation2’,numel(elevation2),1);
azimuth2 = (azimuth2)./pi.*180;
azimuth2 = reshape(azimuth2’,numel(azimuth2),1);
spa = dlmread(‘Data/SPA_DEC_EoT_2015/angle.txt’);
spa_ele = spa(:,1);
spa_azi = spa(:,2);
%special cases
TMP = abs(azimuth-spa_azi)>3;
id = find(TMP);
azimuth(id) = 360-azimuth(id);
azimuth2(id) = 360-azimuth2(id);
disp(mean(abs(elevation-spa_ele)));
disp(mean((elevation2-spa_ele)));
disp(mean(abs(azimuth-spa_azi)));
disp(mean(abs(azimuth2-spa_azi)));
本资料由:数模之家(春哥带你学建模)收集并分享,需要更多数学建模资料添加微信公众号:shumohome
数模之家

  • 3
    点赞
  • 0
    评论
  • 1
    收藏
  • 一键三连
    一键三连
  • 扫一扫,分享海报

©️2021 CSDN 皮肤主题: 大白 设计师:CSDN官方博客 返回首页
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、C币套餐、付费专栏及课程。

余额充值