全国大学生数学建模2015年国二a题

2015高教社杯全国大学生数学建模竞赛
承诺书
我们仔细阅读了中国大学生数学建模竞赛的竞赛规则.
我们完全明白，在竞赛开始后参赛队员不能以任何方式（包括电话、电子邮件、网上咨询等）与队外的任何人（包括指导教师）研究、讨论与赛题有关的问题。

我们知道，抄袭别人的成果是违反竞赛规则的, 如果引用别人的成果或其他公开的资料（包括网上查到的资料），必须按照规定的参考文献的表述方式在正文引用处和参考文献中明确列出。

我们郑重承诺，严格遵守竞赛规则，以保证竞赛的公正、公平性。

如有违反竞赛规则的行为，我们将受到严肃处理。

我们授权全国大学生数学建模竞赛组委会，可将我们的论文以任何形式进行公开展示（包括进行网上公示，在书籍、期刊和其他媒体进行正式或非正式发表等）。

我们参赛选择的题号是（从A/B/C/D中选择一项填写）： A
我们的参赛报名号为（如果赛区设置报名号的话）：
所属学校（请填写完整的全名）：
参赛队员(打印并签名) ：1.
2.
3.
指导教师或指导教师组负责人(打印并签名)：
日期：年月日赛区评阅编号（由赛区组委会评阅前进行编号）：
编号专用页
赛区评阅编号（由赛区组委会评阅前进行编号）：
全国统一编号（由赛区组委会送交全国前编号）：全国评阅编号（由全国组委会评阅前进行编号）：
太阳影子定位
摘要
本文研究了太阳影子定位问题，基于天球坐标系相关知识、球面几何理论以及相似度理论，对不同情况下的数据，建立了相应的数学模型并得到了最优的匹配地点与日期。

问题1中，利用球面三角形余弦定理给出了太阳高度角公式，并建立了影子长度变化的数学模型，定性的分析了影子长度关于时角、当地纬度以及赤纬角的变化规律：(1). 时角的绝对值越大，影子长度越大；(2). 在同一经度上（即时角一定），当地纬度与此时的太阳赤纬之差越大，影子长度越大；(3). 在同一纬度不同经度上，当地经度和此时太阳直射点所在的经度之差越大，影子长度越大。

用所建的模型，得到了2015年10月22日北京时间9:00-15:00之间天安门广场3米高的直杆的太阳影子长度的变化曲线。

问题2中，由太阳高度角、方位角、时角、赤纬角之间的关系，建立了杆的高度、太阳高度角、太阳方位角等变量间的方程组，推导出影子的长度与太阳高度角满足的方程，该方程不依赖于杆的高度、坐标系的选取。

利用附件1的数据并结合遍历相关知识，给出了最优匹配地点：最优地点为北纬24度，东经100度，（云南省沧市临翔区）；次优地点为北纬18度，东经110度（海南省三亚市）。

问题3中，在问题2模型的基础上，对未知日期引入了日期参数，导出了相应的优化模型。

由附件2的数据得出了最优匹配日期和地点：7月10日，北纬41度，东经80度，新疆维吾尔自治区阿克苏地区乌什县。

由附件3的数据得出了最优匹配日期和地点：3月02日，北纬37度，东经110度，陕西省延安市延川县。

问题4中，给出了杆的高度的近似值，我们修正了问题3中的模型，在优化模型中，引入了变量h。

对视频进行逐针读取，再灰度化处理，利用Photoshop 软件得到了影子坐标的和杆的坐标的像素，由它们之间的比例关系，得出了相应的坐标。

利用坐标数据和杆的近似高度，通过优化模型，经过Matlab遍历，得到了最优的拍摄地点：最优的为北纬40度，东经113度，山西省大同市南郊区；其它次优的为北纬44度，东经110度，蒙古；北纬39度，东经110度，陕西省榆林市神木县；北纬39度，东经111度，陕西省榆林市府谷县；北纬39度，东经112度，山西省忻州市神池县；北纬40度，东经114度，河北省张家口市阳原县；北纬41度，东经115度，河北省张家口市崇礼县。

如果拍摄日期未知，最优的拍摄地点为北纬40，东经115，河北省张家口市蔚县，日期为6月22日。

同时，我们也提出了分光理论，并建立了相应的模型和数学推导。

由于时间的限制，不能完善该模型。

针对问题二、三、四，提出了影子弧线的相似度模型。

由于提取数据的误差，可能使得我们的计算结果与实际地点有稍许偏差。

关键字:太阳影子定位，天球坐标系，遍历，时角，赤纬角，太阳高度角，太阳方位角
1. 问题重述
太阳影子定位技术就是通过分析视频中物体的太阳影子变化，确定视频拍摄的地点和日期的一种方法。

1.1问题一有两部分组成：
（1）：确定影子长度变化数学模型，分析影子长度关于各个参数的变化规律。

（2）：用已建立的模型画出2015年10月22日北京时间9:00-15:00之间天安门广场（北纬39度54分26秒,东经116度23分29秒）3米高的直杆的太阳影子长度的变化曲线。

1.2问题二有两部分组成：
（1）：根据某固定直杆在水平地面上的太阳影子顶点坐标数据，建立数学模型确定直杆所处的地点。

（2）：将模型应用于附件1的影子顶点坐标数据，给出若干个可能的地点。

1.3问题三有两部分组成：
（1）：根据某固定直杆在水平地面上的太阳影子顶点坐标数据，建立数学模型确定直杆所处的地点和日期。

（2）：将所建的模型分别应用于附件2和附件3的影子顶点坐标数据，给出若干个可能的地点与日期。

1.4问题有三部分组成：
（1）：根据附件4一根直杆在太阳下的影子变化视频，建立确定视频拍摄地点的数学模型。

（2）：应用已建立的模型给出若干个可能的拍摄地点。

（3）：如果拍摄日期未知，根据视频确定出拍摄地点与日期。

2. 模型假设
（1）：地球是球体；
（2）：直杆影子长度为直线段长度而非弧长；
（3）：到达地球的太阳光线为平行光线；
（4）：太阳高度角为太阳光的入射方向与入射点切平面（地平面）的夹角；（5）：假设地面无障碍物，直杆影子都投影在一个平面上；
（6）：假设附件4中拍摄角度对所提取的影长无影响；
（7）：假设一天当中太阳赤纬角不变；
3. 符号说明
注：规定北纬为正，南纬为负
4. 模型建立及求解
4.1问题一的分析
影子的长度主要依赖于两个因素，即直杆的长度和当地的太阳高度角。

由推导可得，太阳高度角又依赖于当时的太阳赤纬，时角和当地的纬度。

易知当地的纬度和时角；又可知赤纬角由日期确定，注意我们假设一天当中的赤纬角是不变的。

4.1.1问题一模型的建立及求解
直杆高度、太阳高度角及影子长度三者有如下关系：
= (4.1)
tan
r hα
其中，r为影子长度，h为直杆高度，α为太阳高度角，其关系如图4.1所示。

图4.1 直杆高度、太阳高度角及影子长度关系示意图
在天文三角形中，已知两边：测者余纬，天体极距PNB=2
π
δ-，它们的夹
角，即天体半圆时角ω，应用边的余弦公式可以求得第三边：天顶距z=2πα-,
如图4.2所示。

图4.2 天球坐标系中的相关参数
应用球面三角形边的余弦定理
[]
1,2,有
cosa cosbcosc sinbsinccos ω=+, (4.2) 又2a απ+=,2b δπ+=,2c φπ+=，所以
()()()cos 90cos 90cos 90οοοαδφ-=--
()()()sin
90
sin 90sin 90ο
οοδφω+---. (4.3)
即有
sin sin sin cos cos cos αφδφδω=+. (4.4)
带入(4.1)式，可得
r =
. (4.4) 其中，时角ω的计算公式为
[]
3：
()
()1512s t ω=-度 , (4.5)
太阳赤纬角δ
（即太阳直射点的维度）：
()()228423.45sin 365n πδ+⎛⎫
= ⎪
⎝⎭
度 (4.6)
其中n 为日期序号，例如，1月1日为1=n ，2月25日为56n =。

下面我们将定性分析影子长度r 关于各个参数（时角ω、纬度φ及赤纬角δ）的变化规律。

首先，分析在一天之中时角ω对影子长度的影响。

由假设可知赤纬角δ在一天之中是恒定的，又tan r
h α=可知，太阳高度角α越大，影子长度越小；又
由式(4.4)可得，时角ω的绝对值ω越大，α越小。

所以可知，ω越大，影子长度r 越大。

其次，分析纬度φ对影子长度的影响。

我们在同一经度上（即时角ω定）考虑此问题，当地纬度φ与此时的太阳赤纬之差越大，影子长度r 越大。

最后，分析赤纬角δ对影子长度的影响。

我们在同一纬度不同经度上考虑此问题，当地经度和此时太阳直射点所在的经度之差越大，影子长度r 越大。

应用我们所建立的模型(4.4)，基于Matlab plot 函数给出9:00—15:00的影子长度，见图(4.3)。

其中整点时刻的影子长度见表4.1.
表4.1 9:00—15:00间整点时刻直杆影子的长度
北京时间 9 10 11 12 13 14 15 影长(米)
7.8105
5.3652
4.2675
3.8585
3.9843
4.6874
6.2862
图4.3 9:00—15:00天安门广场直杆影子长度的变化曲线
4.2问题二的分析
首先建立一个坐标系（标准坐标系），附件数据所在的坐标系为旧坐标系，标准坐标系与旧坐标系之间有一个旋转角。

通过推导发现，影子长度与相应时刻太阳高度角的余切之比为定值，故可以直接利用旧坐标系下的影长数据，而避开了未知杆长、旋转角及标准坐标系下的坐标分量，从而简化了问题的复杂度和求解难度。

4.2.1问题二模型的建立及求解
设i α为i 时刻的太阳高度角，i β为i 时刻太阳方位角(太阳方位角为正北方向按顺时针旋转到太阳投影点所旋转的角度)
[]
4,i r 为i 时刻影子长度，(),x y 为旧坐标
系下的杆顶点投影的坐标（旧坐标系为附件中数据所在的坐标系），建立标准坐标系，以杆低端为原点，x 轴朝向正南，y 轴朝向正东。

(),x y 为标准坐标系下杆顶点投影的坐标。

在标准坐标系下，有
2
2
tan h
r
x y
α==+ (4.7)
所以
cot h α⋅= (4.8)
又i β为i 时刻太阳方位角，则i πβ-为旧坐标系和标准坐标系之间的旋转角度，故
cot sin()cot cos()i i i
i i i h y h x απβαπβ-=⎧⎨
-=⎩ (4.9)
即
sin()sin tan cos()cos i i i i i i i
y x πβββπββ-===--- (4.10) 也就是直杆影子顶点横纵坐标的比值与杆子的高度无关，仅仅与太阳的方位角有关。

因为
cot i h α== (4.11)
所以
i
h =
=恒定 (4.12)
由上面分析可知，新旧坐标系下的坐标，杆的高度，太阳高度角，太阳方位角，满足如下方程组：
cot cos()cot sin(cos sin sin cos cot )i i i i i i i x x y y h h x h y
απβαπβθθθθα⎡⎤⎡⎤⎡⎤
=⎢⎥⎢⎥⎢⎥--=⎧⎪-⎣⎦⎣⎦⎣⎦
==⎪
⎪=⎨⎪⎪
⎪⎩
(4.13) 上述方程组为含多个未知量的非线性的方程，求解非常困难。

但是由方程组可以得到如下关系式：
cot i α=
=
(4.14)
由此得到
cot cot i
j
j
i
r r αα==
(4.15) 若所求点恰好为附件所给数据的拍摄地点，必满足
()2
,cot cot 0i
j
j i i j
r
r αα⋅-⋅=∑ (4.16)
从而建立目标函数如下：
()2
,min cot cot i j j i i j
f r r αα=⋅-⋅∑ (4.17)
此目标函数和杆的高度、两个坐标系均无关，大大简化了模型的求解难度，加快了求解速度。

为了找出最佳匹配的点，我们需要在球面上，计算出一些点的某段时刻的高度角，带入目标函数，求出最佳匹配点。

在Matlab 环境下，基于遍历思想，求出最佳匹配点如下： 最优：北纬24度 东经100度 （云南省沧市临翔区） 次优：北纬18度 东经110度 （海南省三亚市）
基于目标函数(4.17)，给出影子长度随时间变化的规律，见图(4.4)，其中不同时刻影子的长度见表4.2.
图4.4附件1中影子长度的变化曲线
北京时间 12:41 12:44 12:47 12:50 12:53 12:56 12:59 13:02 13:05 13:08 13:11 影子长度 1.247256
1.222795
1.198921
1.175429
1.15244
1.129917
1.107835
1.086254
1.065081
1.044446
1.024264
北京时间 13:14 13:17 13:20 13:23 13:26 13:29 13:32 13:35 13:38 13:41 影子长度
1.00464
0.985491
0.96679
0.948585
0.930928
0.913752
0.897109
0.880974
0.865492
0.850504
4.3问题三的分析
在问题二中的目标函数里面，由于日期是已知的，即赤纬角是已知的；而问题三中，日期是未知的，即赤纬角是未知的。

故太阳高度角和当地经纬度、日期均有关，这就导致所建的目标函数还是日期的函数，故需要在问题二的基础上，还须对日期进行遍历。

4.3.1问题三模型的建立及求解
建立目标函数如下：
()2
,min cot cot sin sin sin cos cos cos i j j i i j f r r αααφδφδω
⎧=⋅-⋅⎪
⎨
⎪=+⎩∑ (4.18)
在Matlab 平台下，基于遍历思想，求出最佳匹配点如下： 附件二的最优解为：
日期191 7月10日 北纬41度 东经80度 新疆维吾尔自治区阿克苏地区乌什县；
基于目标函数(4.18)，给出影子长度随时间变化的规律，见图(4.5)，其中不同时刻影子的长度见表4.3
.
图4.5附件2中影子长度的变化曲线
北京
时间
12:41 12:44 12:47 12:50 12:53 12:56 12:59 13:02 13:05 13:08 13:11
影子
长度
1.247
256
1.222
795
1.198
921
1.175
429
1.152
44
1.129
917
1.107
835
1.086
254
1.065
081
1.044
446
1.024
264 北京
时间
13:14 13:17 13:20 13:23 13:26 13:29 13:32 13:35 13:38 13:41
影子
长度
1.004
64
0.985
491
0.966
79
0.948
585
0.930
928
0.913
752
0.897
109
0.880
974
0.865
492
0.850
504
日期61 3月02日北纬37度东经110度陕西省延安市延川县基于目标函数(4.18)，给出影子长度随时间变化的规律，见图(4.6)，其中不同时刻影子的长度见表4.4
图4.6附件2中影子长度的变化曲线
北
京
时
间
13:09 13:12 13:15 13:18 13:21 13:24 13:27 13:30 13:33 13:36 13:39
影
子
长
度
3.533
142
3.546
768
3.561
798
3.578
101
3.595
751
3.614
934
3.635
426
3.657
218
3.680
541
3.705
168
3.731
278
北
京
时
间
13:42 13:45 13:48 13:51 13:54 13:57 14:00 14:03 14:06 14:09
4.4问题四的分析
由于影子长度的数据未知，所以首先需要对视频进行逐针读取，再灰度化处理，利用Photoshop 软件提取影子长度数据。

在该问题中，由于给出了直杆高度的估计值，我们便可以进一步优化目标函数，找出最佳匹配点。

4.4.1问题四模型的建立及求解
问题四给出了直杆高度的估计值，基于问题二和问题三的分析，我们可以进一步优化目标函数为：
()2min cot sin sin sin cos cos cos i i i f r h ααφδφδω
⎧=⋅-⎪
⎨
⎪=+⎩∑ (4.19) 在Matlab 环境下，基于遍历思想，求出最佳匹配点如下：
最优：北纬40度 东经113度 山西省大同市南郊区
次优：北纬44度 东经110度 蒙古
北纬 39度 东经110度 陕西省榆林市神木县 北纬 39度 东经111度 陕西省榆林市府谷县 北纬 39度 东经112度 山西省忻州市神池县 北纬 40度 东经114度 河北省张家口市阳原县 北纬 41度 东经115度 河北省张家口市崇礼县
若拍摄日期未知，可以得到最优解，即可以确定拍摄地点和日期，日期
和地点为：
日期207 6月22日 北纬40度 东经115度 河北省张家口市蔚县
5. 模型的改进与评价
5.1问题一模型的改进与评价
原问题一中所建的模型的求解比较复杂，把地球看成了一个圆球体，而地球是一个椭球体，所以该模型舍去了一些高阶微小量，只是给出了一个近似模型。

我们需要注意，可以求得任何一地点在任何时刻的太阳直射点，一般情况下，该直射点和直杆处的经度和纬度是不同的，而我们假设太阳光线为平行光，基于此，可以把直杆处的光线沿着经线的切线方向和纬线的切线方向进行分解，分别得到两个方向的单位化向量，然后将经线切线的单位向量除以直杆处纬度和此时赤纬角的差绝对值的余弦，得到一向量；再将纬线切线的单位向量除以直杆处经度和直射点处的经度的差绝对值的余弦，得到另外一向量，分别将这两个向量加权合成光线向量，进一步可得影子的长度。

假设A 处为太阳直射点，B 处为直杆所在点，把直杆处的光线沿着经线的切线方向和纬线的切线方向进行分解,即分别在同一经度不同纬度下，同一纬度不同经度下考虑，如图5.1所示，假设北纬为正，南纬为负。

图5.1 直杆处光线分解示意图
首先考虑，在同一经度不同纬度下的情况，如图5.2所示。

直杆处纬度和赤纬之差的绝对值为=-γφδ，假设在B 点建立空间直角坐标系如下：如图5.3所示。

再考虑同一纬度不同经度的情况，如图5.4所示。

假设在D 点建立空间直角坐标系如下：如图5.5所示。

设1V 为单位向量，因为太阳光朝向Z 轴负半轴，所以()1V =0sin ,cos αα-，.设2V 为单位向量，因为太阳光朝向Z 轴负半轴,所以()2V =sin ,0,-cos ωω.
图5.2 直杆处光线经线切线分解示意图
图
5.3 B 处空间直角坐标系
图5.4 直杆处光线纬线切线分解示意图
则()11222112sin ,
sin ,cos cos V V V λλλωλαλαλω=+=--
2
12222
2
1
cos cos tan sin sin A λαλω
λωλα
--=
+其中，A 为V 与Z 轴夹角，tan r h A =. 合成光线向量如图5.6所示。

图5.5 D 处空间直角坐标系
()()2cos0sin ,sin0sin ,cos sin ,0,cos V ωωωωω=⋅⋅=
图5.6 C 处空间直角坐标系
5.2问题二、三、四模型的改进与评价
缺点：原模型由于所给数据为一小时内间隔的数据，所以模型不能很好的解决问题，只能得出最优解，而最优解并不一定是最佳答案
改进模型如下：同一位置的直杆，在地面上的投影顶点随着时间的变化形成
一条轨迹曲线，而对于不同长度的两根直杆L1和L2,他们对应影子顶点的轨迹曲线C1和C2是相似的，如图
5.7所示：
图5.7 影子的相似关系
我们做以下处理，把曲线相似性量度化，如下图所示两天曲线[]5：
图5.8 曲线相似关系
用234
1
1231
Z=,,,n
n
ob ob ob ob
ob ob ob ob
-
⎡⎤
⎢⎥
⎣⎦
代表
15
b b这条曲线特征,用
234
2
1231
,,n
n
oa oa oa oa
Z
oa oa oa oa
-
⎡⎤
=⎢⎥
⎣⎦
代表
15
a a曲线特征。

把两曲线相似性转化为矩阵12Z Z 相似性，用相关性系数来量化两矩阵的相关程度。

1Z n n h h -⎡⎤
=⎢⎥⎣⎦
,1n n h Z h φϕ
-⎡⎤=⎢⎥⎣⎦
,Z φϕ为经度为φ纬度为ϕ的地方的14：42至15：42这段时间内的某时刻影子长
度与上一时刻影子比值的得到的一维矩阵。

两个矩阵的相关系数为：
(Z,Z )(ZZ )EZ EZ Cov E d φϕφϕ
-⋅=
=
其中d φϕ表示经度φ维度ϕ的地方所得到的一维矩阵与附表所得到的基准矩阵的相关系数，由matlab 程序按遍历的索引方法对全球进行检索。

经度每隔1度，纬度每隔1度去求解，并建立一个二维矩阵M 去存放d φϕ这个二维矩阵行用φ，列用ϕ表示
1212
121
21
2
n n n n n n d d d d d M d d d φϕφϕφϕφϕφϕφϕφϕφϕ⎡⎤⎢⎥
⎢
⎥=⎢⎥⎢⎥⎣⎦
我们取出M 中前最大的元素，由大元素中定位这个点，为求全球各地在2015年4月18日14：42至15：42固定直杆在水平面上的太阳影长变化规律模型如下：
tan h
r α
=
其中α为太阳高度角，φ为当地纬度，δ为当时的太阳赤纬，ω为时角。

参考文献
[1] 杨世国，球面几何简介(Ⅱ)，中学数学教学，4：1-4，2005年4月。

[2] 周雪娟，球面三角形公式及其应用，浙江国际海运职业技术学院学报，4(1): 59-63, 2008年3月。

[3] 方荣生，太阳能应用技术，pp: 36-48，北京：中国农业出版社，1985年。

[4] 贺晓雷，于贺军，太阳方位角的公式求解及其应用，太阳能学报，29(1): 69-73,
2008年1月。

[5] 李秀敏,江卫华，相关系数与相关性度量，数学的实践与认识，36(12):188-192,2006年12月。