http://www.ma-xy.com
目录
第一章 太阳影子定位 1
1.1 题目要求 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
1.1.1 附件 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
1.2 太阳影子定位研究 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2
1.2.1 模型的假设 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2
1.2.2 符号说明 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2
1.2.3
问题一的分析与求解
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2
1.2.4 问题二的分析与求解 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
1.2.5 问题三的分析与求解 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
1
http://www.ma-xy.com
第一章 太阳影子定位
1.1 题目要求
A 题：太阳影子定位
如何确定视频的拍摄地点和拍摄日期是视频数据分析的重要方面，太阳影子定位技术就是通
过分析视频中物体的太阳影子变化，确定视频拍摄的地点和日期的一种方法。
1. 建立影子长度变化的数学模型，分析影子长度关于各个参数的变化规律，并应用你们建立
的模型画出 2015 年 10 月 22 日北京时间 9:00-15:00 之间天安门广场 (北纬 39 度 54 分 26 秒,
东经 116 度 23 分 29 秒)3 米高的直杆的太阳影子长度的变化曲线。
2. 根据某固定直杆在水平地面上的太阳影子顶点坐标数据，建立数学模型确定直杆所处的地
点。将你们的模型应用于附件 1 的影子顶点坐标数据，给出若干个可能的地点。
3. 根据某固定直杆在水平地面上的太阳影子顶点坐标数据，建立数学模型确定直杆所处的
地点和日期。将你们的模型分别应用于附件 2 和附件 3 的影子顶点坐标数据，给出若干个可能
的地点与日期。
4．附件 4 为一根直杆在太阳下的影子变化的视频，并且已通过某种方式估计出直杆的高度
为 2 米。请建立确定视频拍摄地点的数学模型，并应用你们的模型给出若干个可能的拍摄地点。
如果拍摄日期未知，你能否根据视频确定出拍摄地点与日期？
1.1.1 附件 1
说明：坐标系以直杆底端为原点，水平地面为 xy 平面。直杆垂直于地面。测量日期：2015
年 4 月 18 日。附件 1 数据如表 (1.1) 所示
1
http://www.ma-xy.com
1.2 太阳影子定位研究 第一章 太阳影子定位
表 1.1: 太阳影子定位：附件 1
t
i
(北京时间) l
i
影长 (米) t
i
(北京时间) l
i
影长 (米)
14:42 1.1496 15:15 1.5402
14:45 1.1822 15:18 1.5799
14:48 1.2153 15:21 1.6201
14:51 1.2491 15:24 1.6613
14:54 1.2832 15:27 1.7033
14:57 1.3180 15:30 1.7462
15:00 1.3534 15:33 1.7901
15:03 1.3894 15:36 1.8350
15:06 1.4262 15:39 1.8809
15:09 1.4634 15:42 1.9279
15:12 1.5015
1.2 太阳影子定位研究
1.2.1 模型的假设
1. 假设研究的物体是垂直于该地点的切平面；
2. 影子长度仅考虑切平面上的投影长，不考虑弧度；
3. 假设地球为球体；
4. 假设一天时间内，太阳直射点的纬度不变 (赤纬角在一天内保持不变)；
5. 不同年份的同一时间、同一地点的太阳高度角是相同的；
6. 假设地球绕太阳公转的轨道为圆形轨道。
1.2.2 符号说明
1.2.3 问题一的分析与求解
问题的分析
建立影子长度变化的数学模型，分析影子长度关于各个参数的变化规律，并应用你们建立的
模型画出 2015 年 10 月 22 日北京时间 9:00-15:00 之间天安门广场 (北纬 39 度 54 分 26 秒, 东
经 116 度 23 分 29 秒)3 米高的直杆的太阳影子长度的变化曲线。
第一问的主要目标是让我们确定要分析的量，并建立各个量之间的函数关系。
http://www.ma-xy.com 2 http://www.ma-xy.com
http://www.ma-xy.com
第一章 太阳影子定位 1.2 太阳影子定位研究
模型的建立与求解
我们要通过物体影子的数据 (各时刻的长度变化和方向变化) 来判定物体的位置 (经纬度) 及
拍照时间。某固定杆高的影子长度在不同地点和时点下应该是不同的。可以推测影长是关于杆
高、日期、时间、经度和纬度的函数，事实也确实如此。下面，我们找出各变量间具体的函数关
系。在建立函数关系式之前，我们对太阳直射点、赤纬角和时角等相关概念做了了解，并给出了
某日某时杆高为 h 的观测物在地球某点 O 处的太阳影长变化图，如图 (1.1) 所示
图 1.1: 太阳影子长度模拟图
假设研究的物体是垂直于该地点的切平面；并且影子长度仅考虑切平面上的影长，不考虑弧
度。由于地球半径很大，杆子在地球表面的影长 l
0
可以用杆子所在切平面上的影长 l 做近似代
替，由此可建立杆高 h(观测物的高度) 关于影长 l(真实影长 l
0
的估计值) 和太阳高度角 ρ 的关
系式
l =
h
tan
ρ
=
h
1 − sin
2
ρ
sin
ρ
其中：ρ 是太阳光线与拍摄地点切平面的交角，称为太阳高度角。我们设经度为 θ，纬度为 φ，年
份为 y，日期为 n，时刻为 t，杆高为 h，影长为 l。应该可以想到的是 ρ 是 φ, θ, n, t, y 的函数，
那么，这种函数关系到底是怎样的呢？为简单处理，假设不考虑 y 对 ρ 的影响，即不同年份的同
一时间、同一地点的太阳高度角是相同的。下面，我们来求解 ρ 或者 sin ρ。
假设地球是一球体，以其球心为坐标原点，以赤道面作为 xoy 面，以垂直 xoy 向北的方向
为 z 轴方向建立三维直角坐标系 oxyz，如图 (1.2) 所示
图 1.2: 地球坐标系
图中 w 为太阳时角，δ 为太阳直射点维度 (太阳赤纬)，O 为待求地点，φ 为 O 点的纬度。
http://www.ma-xy.com 3 http://www.ma-xy.com
http://www.ma-xy.com
1.2 太阳影子定位研究 第一章 太阳影子定位
我们假设一天内太阳直射点的维度不变，在给出 ρ 之前，有一些常识要说明
1. 太阳直射点的太阳高度角为 90
◦
；
2. 太阳赤纬是太阳直射点的纬度 δ，其一年内的变化情况如图 (1.3) 示意图
图 1.3: 太阳赤纬一年的变化情况
3. 太阳光线到达地球时时平行的 (这点需要假设)；
4. 以太阳直射点为圆心的圆上各地点的太阳高度角近似相等；
5. 自转改变太阳直射点的经度，公转改变维度。
我们知道，太阳高度角 ρ 是线 (太阳光线) 与面 (待求地点 O 的切平面) 的夹角，因此，我
们需要知道平面的法向量以及线的向量。O 点所在的切平面的单位法向量为
n = (cos φ cos w , cos φ sin w, sin φ)
而第 n 天 t 时刻的入射向量为
m = (−cos δ, 0, −sin δ)
则 O 点的太阳高度角 ρ(如图 (1.4) 所示) 求解为
sin ρ = −n · m = cos φ cos δ cos w + sin φ sin δ
其中：δ 为太阳赤纬，w 为太阳时角。
图 1.4: 太阳照射点的投影向量
http://www.ma-xy.com 4 http://www.ma-xy.com
http://www.ma-xy.com
第一章 太阳影子定位 1.2 太阳影子定位研究
关于赤纬角 δ 的计算方法有：Cooper 法、Spencer 法、Stome 法、Bpirges 法等；而关于太阳
时角 w 的计算方法亦有：Wloof 法、Spencer 法、Whillier 法、Lamm 法等。杜春旭等人在《一
种高精度太阳位置算法》一文中对上述算法进行了分类比较，显示在每天世界 0 时时，Bpirges
法和 Lamm 法求解的 δ 和 w 的误差较小，结果较为精确。现在我们引用 Spencer 在 1971 年给
出的赤纬角 δ(弧度) 的计算公式
δ = 0.006918 − 0.399912 cos(Γ) + 0.070257 sin(Γ) − 0.006758 cos(2Γ)
+ 0.00907 sin(2Γ) − 0.002697 cos(3Γ) + 0.00148 sin(3Γ)
其中：Γ = 2π(n − 1)/365，单位为弧度，n 表示一年中的第 n 天即日期，如 1 月 1 日计 n = 1，
12 月 31 日的计 n = 365。
太阳赤纬角：下面，我们尝试给出太阳赤纬 δ 的求解。设地球 P 绕太阳 S 公转的轨道为圆
形轨道。太阳位于圆形轨道的圆心，以太阳为原点，以圆形轨道面 (黄道平面) 为 xoy 面家里公
转坐标系，如图 (1.5) 所示
图 1.5: 公转坐标系
设一年有 T
∗
天，一天有 t
∗
小时，以春分 (3 月 21 日) 为第 0 天 (n = 0)，则第 n 天 t 时
刻，地球的坐标为
P (r cos
2πT
T
∗
, r sin
2πT
T
∗
, 0)
其中：r 为地球中心到太阳中心的距离，T = n +
t
t
∗
。
设第 n 天 t 时刻，太阳直射地球的维度 δ(太阳赤纬)，已知在春分 t
0
时刻太阳直射赤道，即
当 n = 0, t = t
0
时，δ = 0，此时的 T 为 T
0
t
0
t
∗
。而地球中心坐标为
P
0
(r cos
2πT
0
T
∗
, r sin
2πT
0
T
∗
, 0)
设北回归线的纬度为 δ
0
，则北极方向是 (0, 0, 1) 绕 SP
0
轴旋转 δ
0
后的方向。在地球绕太阳公转
过程中，北极方向始终保持不变，因而
−−→
P N =
sin
2πT
0
T
∗
sin δ
0
, −cos
2πT
0
T
∗
sin δ
0
, cos δ
0
−→
P S 与
−−→
P N 的夹角为
π
2
− δ，因此
cos(
π
2
− δ) =
−→
P S ·
−−→
P N
|
−→
P S||
−−→
P N|
http://www.ma-xy.com 5 http://www.ma-xy.com
http://www.ma-xy.com
1.2 太阳影子定位研究 第一章 太阳影子定位
即
sin δ = −cos
2πT
T
∗
sin
2πT
0
T
∗
sin δ
0
+ sin
2πT
T
∗
cos
2πT
0
T
∗
sin δ
0
= sin
2π
n +
t−t
0
t
∗
T
∗
sin δ
0
为简单，假设赤纬角在一天内保持不变，则上式可写为
sin δ = sin
2πn
T
∗
sin δ
0
太阳时角：定义太阳时角 w 在正午时为 0
◦
，每隔一小时增加 15
◦
，上午为正，下午为负。但
考虑到实际两地之间存在时差并对时角产生影响，且中国各地区均以北京时间为标准 (北京位于
120
◦
E，所在时区为东 8 区) 因此某地以北京时间为依据的时角 w(弧度化) 的计算公式为
w =
π(12 − t)15
◦
180
◦
+
π
180
◦
120
◦
− θ θ ∈ [0, 180
◦
E]
120
◦
+ θ θ ∈ [0, 60
◦
W ]
240
◦
− θ θ ∈ [60, 180
◦
W ]
其中：w 为时角，t 为时间 (以小时为单位，可为小数，如 13.326)， 为经度。
通过以太阳高度角赤纬角 δ 和太阳时角 w 为桥梁，搭建了影长 l 关于杆高 h、日期 n、时间
t、经度 θ 和纬度 φ 的具体函数关系式。依据此关系式，可以解决下列关于太阳影子定位的问题：
1. 给出具体的 φ、θ、n、h，观察 l 随 t 的变化；
2. 给出具体的 n、h、t，观察 l 随 φ、θ 的变化；
3. 给出一组 (t
i
, l
i
)i = 1, 2, . . . , k，(k 为观测个数，t
i
, l
i
皆可在具体的视频中求得)，在 n 已知
h 未知但固定的条件下，求出视频的拍摄位置 φ, θ；
4. 给出一组 (t
i
, l
i
)i = 1, 2, . . . , k，h 未知但固定的条件下，求出视频的拍摄位置 φ, θ 和拍摄
日期 n；
针对问题 1 和问题 2，可通过各参数间的函数关系式绘制图像进行分析，相对简单。问题 4
在问题 3 的基础上增加了一个参数 n 的求解，因此，我们将详细讨论问题 3 的求解，简要分析
问题 4。对模型求解之前，我们做如下规定：由于问题中的参数很多，均以如下表达式规范参数
y = f(x; α, β|C)
其中：y 表示输出变量 (可以是标量也可以是向量)，x 表示出入变量 (可以是标量也可以是向量)，
α, β 表示参数，C 表示某变量的固定值。
程序
Matlab 给出了一个求解影长变化的示例
¬
，可以作为参考。下面，我们给出 3m 长杆 9 点到
15 点太阳影子长度变化求解程序，其中 Rad2deg 函数是将角度转化为弧度，numdate 函数是用
于计算天数 n，求解程序如下：
¬
http://cn.mathworks.com/help/matlab/examples/live-editor-intteractive-narrative.html
http://www.ma-xy.com 6 http://www.ma-xy.com
http://www.ma-xy.com
第一章 太阳影子定位 1.2 太阳影子定位研究
1 phi = Rad2deg (39 ,5 4 ,26 ) ;%纬 度 （ 角 度 转 化 为 弧 度 ）
2 the ta = Rad2deg (11 6 ,23 , 29) ;%经 度
3 h = 3 ;%杆 高
4 n = numdate( 2015 , 10 ,22 ) ;%天 数 n
5 t = [ 9 : 0 . 2 : 1 5 ] ;
6 f o r i = 1 : le ngth ( t )
7 Gamma = 2* p i *(n−1)/3 65;%弧 度
8 de lta = 0.006918 − 0.399912* cos (Gamma) +0.070257* s i n (Gamma) . . .
9 −0.006758* cos (2*Gamma)+0.000907* sin (2*Gamma) − 0.002697* cos (3*Gamma)+0.00148*
sin (3*Gamma) ;%弧 度
10 w = (12− t ( i ) ) *15* pi /180+Rad2deg (3 , 36 , 31 ) ; %弧 度
11 sin_rho = s in ( phi ) * si n ( d elt a )+cos ( phi ) *c os (w) * cos ( de lta ) ;
12 l ( i ) = h* s qr t(1− sin_rho ^2)/ sin_rho ;
13 end
14 f i g u r e
15 p lot ( t , l , ’ r− ’ )
16 se t ( gca , ’ XTick ’ , [ 9 : 1 5 ] )
17 se t ( gca , ’ XTickLabel ’ ,{ ’ 9: 00 ’ , ’ 10:0 0 ’ , ’ 1 1:00 ’ , ’ 12: 00 ’ , ’ 13:00 ’ , ’ 1 4:00 ’ , ’ 15: 00 ’ })
18 xl a be l ( ’ 北 京 时 间 ’ )
19 yl a be l ( ’ 杆 影 子 长 /m’ )
20 t i t l e ( ’ 3m长 杆 9 点 到 15 点 太 阳 影 子 长 度 变 化 曲 线 图 ’ )
21
结果
3m 长杆 9 点到 15 点太阳影子长度变化曲线图如图 (1.6) 所示
图 1.6: 3m 长杆 9 点到 15 点太阳影子长度变化曲线图 (无插值)
1.2.4 问题二的分析与求解
问题的分析
根据某固定直杆在水平地面上的太阳影子顶点坐标数据，建立数学模型确定直杆所处的地
点。将你们的模型应用于附件 1 的影子顶点坐标数据，给出若干个可能 的地点。
http://www.ma-xy.com 7 http://www.ma-xy.com
http://www.ma-xy.com
1.2 太阳影子定位研究 第一章 太阳影子定位
通过抓帧工具获得某一视频中的 (t
i
, l
i
) 数据 k 组 (i = 1, 2, . . . , k)，且视频的拍摄日期 n 已
知，杆高 h 固定 (可分为已知和未知)，利用 k 组 (t
i
, l
i
) 数据求解视频的拍摄位置 (φ, θ)。针对此
问题，基本思想是利用最小二乘法求解 l = f(t; φ, θ|n, h) 中的参数 (φ, θ)，并给定输入变量 t 和
输出变量 l 的 k 组观测值 (t
i
, l
i
)(i = 1, 2, . . . , k) 进行求解，使总离差平方和最小的 (φ, θ) 即为所
求。然而这种方法并不适用于杆高 h 未知的情况，即使转化为 h = f(t, l; φ, θ|n) 仍无法求解，因
此需重新考虑。既然各变量关系已知，给出了 k 组 (t
i
, l
i
) 的观测值，那么每代入一组 (t
i
, l
i
) 都
会有一个杆高 h 的估计值 h
i
，相应的会得到 k 组 h
i
，由于真实杆高 h 未知但固定，可详细分析。
模型的建立与求解
方法 1：最小二乘模型
根据附件 1 中的点坐标 (x
i
, y
i
) 我们可以求出杆的影长
ˆ
l
i
(i = 1, 2, . . . , k)。当然，我们还可
以根据影长公式
l
i
=
h
√
tan ρ
i
(1.1)
我们可以求解 (φ, θ) 使真实值 l
i
与测量值
ˆ
l
i
的离差平方和最小 (最小二乘)，即
min
φ,θ
J(φ, θ) =
k
i=1
(l
i
−
ˆ
l
i
)
2
但是，上面计算影长的公式 (1.1) 的杆高 h 未知，这使得问题变得困难。其实，这里我们可
以任意设置一个杆高 h(只要有就可以了)，在计算的过程中杆高的影响会消失，我们仍能求解出
最佳的地点 (φ, θ)
∗
。
方法 2：比例模型
上面的分析中杆高 h 未知，我们可以想办法消去杆高 h。将各影长相除，我们研究影长变化
的改变。考虑相邻两个影长的变化，有
min
φ,θ
J(φ, θ) =
k−1
i=1
l
i+1
l
i
−
ˆ
l
i+1
ˆ
l
i
2
上面的模型中不再要求 h 已知。其实，我们还可以建立更一般的目标
min
φ,θ
J(φ, θ) =
k
j=i+1
k−1
i=1
l
j
l
i
−
ˆ
l
j
ˆ
l
i
2
方法 3：最小方差模型
最小方差模型是基于下面的思想：真正的视频拍摄位置 (φ, θ) 会使得 ∀i, j ∈ [1, k] 有 h
i
= h
j
成立，即杆高 h 的估计值 h
i
每次求解都相等 h
1
= h
2
= ··· = h
k
。越接近真实的 (φ, θ) 估计值
h
i
, h
j
之间的偏差越小，即
min
φ,θ
J(φ, θ) =
i,j
|h
i
− h
j
|
http://www.ma-xy.com 8 http://www.ma-xy.com
http://www.ma-xy.com
第一章 太阳影子定位 1.2 太阳影子定位研究
但上面目标处理起来并不是很方便，我们将其转化为下面的优化问题
min
φ,θ
1
k
k
i=1
(h
i
− E(h
i
))
2
= Var(h
i
) unknown h
1
k
k
i=1
(h
i
− h)
2
known h
其中：h 表示已知的杆高，h
i
为杆高的估计值，E(h
i
) 为杆高的估计值的期望，V ar(h
i
) 为杆高
的估计值的方差，k 为观测组数。
要对上述优化模型进行求解，提出的简单直观的想法是：遍历地球上所有的地点 (φ, θ)，每
个 (φ, θ) 都会有 h
i
，找到 min V ar(h
i
) 所在的 (φ, θ)
∗
即可。这种做法的求解结果虽然精确，但
计算量却异常大；且因为需要遍历更多的点所以精度越高计算量越大。基于这种思想，且考虑到
计算的精度和时间，设计了粗劣搜索算法 + 精确搜索算法对问题进行求解。
上面建立了 3 个不同的优化模型，下面我们来考虑相应的约束条件。首先，纬度应该在 ±90
◦
之间，即 −90
◦
< φ < 90
◦
；此外，直杆的影子只有在存在太阳直射的地方才能出现，所以要约束
ρ > 0。由此有最终的优化模型
min
φ,θ
J(φ, θ)
s.t.
− 90
◦
< φ < 90
◦
ρ
i
> 0
i = 1, 2, . . . , k
程序
根据影子长度的变化，还可以判断出当地时间是上午还是下午，影子长度单调递增是下午。
对于上述的 3 个优化模型，可以使用牛顿法、拟牛顿法或者遗传算法等进行求解。还可以采用搜
索算法进行求解，搜索法本质上是枚举，将优化变量按一定精度离散化，得到一系列离散值，对
所有这些离散值进行枚举，寻找其中的最优解及最优值。搜索算法编程较为容易，但工作量大，
为了减少搜索时间，可以采用变步长搜索法，即先用较大的精度对优化变量离散，搜索得到其中
的最优解，然后缩小精度，在最优解附近按更小精度离散化后再进行搜索。
问题二的求解程序如下，这里我们仅给出模型一：最小方差模型的求解程序
1 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%模 型 一 ： 最 小 方 差 模 型%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2 c lc , c l e a r
3 %部 分 量 输 入
4 n = numdate( 201 5 , 4 , 18) ;%日 期
5 f o r i =1:6
6 t ( i )=time2h (14 ,42+3*( i −1)) ;
7 end
8 f o r i =7:21
9 t ( i )=time2h (15 ,0+3*( i −7)) ;
10 end
11 l = [1. 1496 25826 1.182198976 1.215296955 1.249051052 1.28319534 1.317993149 . . .
http://www.ma-xy.com 9 http://www.ma-xy.com
http://www.ma-xy.com
1.2 太阳影子定位研究 第一章 太阳影子定位
12 1.353364049 1.389387091 1.426152856 1.463399853 1.501481622 1.540231817
1 .5 7 9 8 53 3 1 6 .. .
13 1.620144515 1.661270613 1.703290633 1.74620591 1.790050915 1.835014272 1.880875001
1.927918447 ] ;
14 a = [ t ’ , l ’ ] ;
15 k = s i z e ( a , 1) ;
16 %% 遍 历 式 粗 搜 索C
17 var_h = z ero s (90 , 180) ;
18 f o r phi = 1: 1 : 9 0 %纬 度
19 f o r the ta = 1: 1: 1 80 %经 度
20 f o r i = 1 : k
21 t = a( i , 1) ;
22 l = a( i , 2) ;
23 Gamma = 2* p i *(n − 1) /3 65; %弧 度
24 d el t a = 0.006918 − 0.399912* c os (Gamma) + 0.070257* s in (Gamma) . . .
25 −0.006758* cos (2*Gamma) + 0.000907* s in (2*Gamma) − 0.002697* cos (3*Gamma) +
0.00148* s i n (3*Gamma) ;%弧 度
26 w = ((12− t ) *15 + (120− theta ) ) * pi /1 80; %弧 度
27 sin_rho = si n ( Rad2deg( phi , 0 , 0) ) * s i n ( de lt a ) + cos (Rad2deg ( phi , 0 , 0) ) *
cos (w)*co s ( de l ta ) ;
28 h( i ) = l *sin_rho/ s q rt (1 − sin_rho ^2) ;
29 end
30 h ;
31 var_h( phi , theta ) = var (h) ;
32 end
33 end
34 best_h = min(min(var_h) ) ;
35 di sp ( ’ 粗 糙 区 域 的 经 纬 度 为 ： ’ )
36 [ cbest_phi , cbest_theta ] = f i nd ( var_h == best_h )
37 di sp ( ’ 粗 糙 区 域 的 弧 度 为 ： ’ )
38 cbe stx = Rad2deg ( cbest_phi , 0 , 0) %纬 度
39 cbe sty = Rad2deg ( cbest_theta , 0 , 0) %经 度
40 %% 遍 历 式 精 搜 索 S
41 ep s i lon = 0 . 00 1 ; %设 置 遍 历 搜 索 精 度
42 x = cbe stx − 0. 05 : ep s i lo n : cbestx + 0 . 0 5 ; %搜 索 范 围
43 y = cbe sty − 0. 05 : ep s i lo n : cbesty + 0 . 0 5 ;
44 svar_h = ze r os ( len gth (x ) , length (x ) ) ;
45 f o r i i i = 1 : length ( x)
46 f o r i i = 1 : length (y )
47 f o r i = 1 : k
48 t = a( i , 1) ;
49 l = a( i , 2) ;
50 phi = x ( i i i ) ; %纬 度
51 theta = y( i i ) ; %经 度
52 Gamma = 2* p i *(n − 1) /3 65; %弧 度
53 d el t a = 0.006918 − 0.399912* c os (Gamma) + 0.070257* s in (Gamma) . . .
54 −0.006758* cos (2*Gamma) + 0.000907* s in (2*Gamma) − 0.002697* cos (3*Gamma) +
0.00148* s i n (3*Gamma) ;%弧 度
55 w = ((1 2 − t ) *15) * pi/180+ (Rad2deg (120 , 0 , 0) − theta ) ; %弧 度
56 sin_rho = si n ( phi ) * s i n ( del ta ) + co s ( phi ) * cos (w) * cos ( d elt a ) ;
57 ch ( i ) = l *sin_rho / s qr t (1 − sin_rho ^2) ;
58 end
59 ch ;
http://www.ma-xy.com 10 http://www.ma-xy.com
http://www.ma-xy.com
第一章 太阳影子定位 1.2 太阳影子定位研究
60 svar_h ( i i i , i i ) = var ( ch ) ;
61 end
62 end
63 sbest_h = min(min( svar_h ) ) ;
64 di sp ( ’ 精 确 区 域 的 弧 度 为 ： ’ )
65 [ ix , i y ] = f i nd ( svar_h == sbest_h ) ;
66 s best x = x( i x )
67 s best y = y( i y )
68 di sp ( ’ 精 确 区 域 的 经 纬 度 为 ： ’ )
69 [ sbest_phi_jiao , sbest_phi_fen , sbest_phi_miao ] = Deg2rad ( sbe stx )
70 [ sbest_theta_jiao , sbest_theta_fen , sbest_theta_miao ] = Deg2rad ( s besty )
71 %% 设 置 初 始 种 群 的GA精 确 查 找
72 %产 生 初 始 种 群 ， 十 进 制 ， 由 上 面 程 序 产 生 的 一 些 粗 劣 的 区 域 和 一 些 随 机 生 成 的 种 群 组 成
73 %phi 的 取 值 为 0−90， th eta 的 取 值 为 0−180， 随 机 生 成 点 。
74 %
75 %依 据 cbestx ; cbes ty 构 建 十 进 制 初 始 种 群
76 N1 = 10 ; %种 群 大 小 为 100 = 10+90
77 n1 = 2; %变 量 个 数
78 Chrom01 = ze r os (N1, 2) ;
79 f o r i =1:N1
80 Chrom01( i , 1 ) = Rad2deg ( cbest_phi + round ( ( rand*2−1)*5) , rand *60 , rand *60) ;% 角 度
弧 度 化+round (( rand*2−1)*1)
81 Chrom01( i , 2 ) = Rad2deg ( cbest_theta + round (( rand*2−1)*5) , rand *60 , rand *60) ;% 角
度 弧 度 化+round ( ( rand*2−1)*1)
82 end
83 N2 = 90 ; %添 加 的 种 群 的 大 小
84 Chrom02 = ze r os (N2, 2) ;
85 f o r i =1:N2
86 Chrom02( i , 1) = cbestx + 0.1 * rand ;
87 Chrom02( i , 2) = cbesty + 0.1 * rand ;
88 end
89 Chrom0 = [ Chrom01 ; Chrom02 ] ; %100 行 2 列
90 % GA求 解
91 f i t n e s s f c n = @(Chrom) GA_objfun_wen2(Chrom, a ,n) ;
92 A= [ ] ; b = [ ] ;
93 LB = [ 0 , 0 ] ;
94 UB = [ pi /2 , pi ] ;
95 nonlcon = [ ] ;
96 IntCon = [ ] ;
97 % optimoptions 的 设 置 h ttps :/ / cn . mathworks . com/help/gads/ ge net ic−algorithm−op tio ns . html
#f14223
98 o ptions = optimoptions ( ’ ga ’ , ’ PopulationType ’ , ’ doubleVector ’ , . . .
99 ’ Popula tionS ize ’ , 100 , ’ I ni ti al Po pu la ti on Ma tr ix ’ , Chrom0 , . . .
100 ’ PlotFcn ’ , { @ga plotbe stf }) ;
101 GA_xy = ga ( f i t n e s s f c n , n1 , A, b , [ ] , [ ] , LB, UB, nonlcon , IntCon , options ) ;
102 di sp ( ’GA精 确 区 域 的 弧 度 为 ： ’ )
103 GA_xy
104 di sp ( ’GA精 确 区 域 的 经 纬 度 为 ： ’ )
105 [ GAbest_phi_jiao , GAbest_phi_fen , GAbest_phi_miao ] = Deg2rad (GA_xy( 1) )
106 [ GAbest_theta_jiao , GAbest_theta_fen , GAbest_theta_miao ] = Deg2rad (GA_xy( 2) )
107 %% 不 设 置 初 始 种 群 的GA查 找
108 o ptions = optimoptions ( ’ ga ’ , ’ PopulationType ’ , ’ doubleVector ’ , . . .
109 ’ Popula tionS ize ’ , 500 , ’ Con straintTol erance ’ ,1 e −3, ’ FunctionTolerance ’ ,1 e − 10 , .. .
http://www.ma-xy.com 11 http://www.ma-xy.com
http://www.ma-xy.com
1.2 太阳影子定位研究 第一章 太阳影子定位
110 ’ PlotFcn ’ , { @ga plotbe stf }) ;
111 GA_xy_c = ga ( f i t n e s s f c n , n1 , A, b , [ ] , [ ] , LB, UB, nonlcon , IntCon , options ) ;
112 di sp ( ’GA精 确 区 域 的 弧 度 为 ： ’ )
113 GA_xy_c
114 di sp ( ’GA精 确 区 域 的 经 纬 度 为 ： ’ )
115 [ GAbest_phi_jiao , GAbest_phi_fen , GAbest_phi_miao ] = Deg2rad (GA_xy_c(1 ) )
116 [ GAbest_theta_jiao , GAbest_theta_fen , GAbest_theta_miao ] = Deg2rad (GA_xy_c( 2) )
117 %% GA查 找 其 他 可 能 的 位 置 ( 经 纬 度 扩 展 )
118 LB = [− pi /2 , 0 ] ;
119 UB = [ pi /2 , pi ] ;
120 GA_maybe_xy = ga ( f i t n e s s f c n , n1 , A, b , [ ] , [ ] , LB, UB, nonlcon , IntCon , opt ion s ) ;
121 di sp ( ’GA精 确 区 域 的 弧 度 为 ： ’ )
122 GA_maybe_xy
123 di sp ( ’GA精 确 区 域 的 经 纬 度 为 ： ’ )
124 [ GAbest_phi_jiao , GAbest_phi_fen , GAbest_phi_miao ] = Deg2rad (GA_maybe_xy( 1) )
125 [ GAbest_theta_jiao , GAbest_theta_fen , GAbest_theta_miao ] = Deg2rad (GA_maybe_xy( 2) )
126
结果
我们选取 2015 年全国大学生数学建模竞赛 A 题中附件 1 的数据 (见表1.1) 进行求解。已知
视频拍摄日期为 2015 年 4 月 18 日，即 n = 108，杆高 h 未知但固定，时间 t 和影长 l 的一系
列观测值 (t
i
, l
i
) 如表 (1.1) 所示，共计 21 个观测值，即 k = 21，并且视频真实的拍摄地点为真
实拍摄地 18.3
◦
N, 109.5
◦
E，这可以将求解后的地点与之比较。
利用上面的最小方差法对问题进行求解，结果如表 (1.2) 所示。
表 1.2: 太阳影子问题二结果表
变量 粗搜索 精搜索 GA(有初始种群) GA(无初始种群) 可能地点 真实拍摄地
纬度 20
◦
N 19
◦
22
′
11.0871
′′
N 19
◦
20
′′
57.2972
′′
N 19
◦
20
′
41.9352
′′
N −80
◦
N 18.3
◦
N
经度 108
◦
E 108
◦
51
′
33.9721
′′
E 108
◦
52
′
54.4503
′′
E 108
◦
53
′
13.0378
′′
E 49
◦
E 109.5
◦
E
注：杆高 h 未知、日期 n = 108 天 (2015 年 4 月 18 日)，影长 l
i
和时间 t
i
已有观测值.
遗传算法进化图如图 (1.7a) 所示，不设置初始种群的遗传算法进化图如图 (1.7b) 所示
(a) (b)
图 1.7: 遗传算法迭代图
http://www.ma-xy.com 12 http://www.ma-xy.com
http://www.ma-xy.com
第一章 太阳影子定位 1.2 太阳影子定位研究
1.2.5 问题三的分析与求解
问题的分析
根据某固定直杆在水平地面上的太阳影子顶点坐标数据，建立数学模型确定直杆所处的地点
和日期。将你们的模型分别应用于附件 2 和附件 3 的影子顶点坐标数据，给出若干个可能 的地
点与日期。
和第二问不同的是，这里不仅要求杆的地点 φ, θ，还要求拍摄日期 n。这里我们仍然可以使
用第二问中的三个优化模型，不过这里会增加一个新的优化变量 (日期 n)，并且这个参数是整数
型变量，与经纬度的实值型变量不同。对于这一点，我们要重新讨论。
如果采用搜索法求解，只需要对日期也进行搜索即可：得到 365 天的最优解及最优值，然后
在得到的 365 个最优值中选取最优解 (这种方法好像有些问题)。如果采用牛顿法，不应将日期 n
作为实值型变量，可将优化算法与搜索法结合，对日期进行枚举，对每一天采用优化算法求解。
模型的建立与求解
单纯采用太阳影子长度进行求解的效率不是很高，因此，我们需要更多的影子信息。注意到，
根据附件中的坐标信息 x
i
, y
i
，我们还可以求出影子的角度。但是附件数据并未说明坐标轴是如
何选取的，这使得问题变得困难。像未知高度 h 一样，我们可以将方位角进行相除，得到方位角
变化情况。仅考虑影子角度的相对变化 (即影子角度的改变量)，这会方便我们建立坐标系 xOy，
因为不需要过多关注 x 和 y 轴的真实方向，建立影子角度变化模拟图如图 (1.8)
图 1.8: 太阳影子角度变化模拟图
图 (1.8) 中的 h 代表杆高，O 为杆高地点，l
k
为第 k 个影长观测值，(x
i
, y
i
) 表示第 i(i =
1, 2, . . . , k) 个影长在 xOy 坐标系的顶点坐标，ˆγ
i
为第 i 个影子角度的大小 (影角) 的观测/估计，
其计算公式为
ˆγ
i
= arctan
y
i
x
i
上面给出了影角的观测值求解方法，下面，我们给出影角的理论公式。考虑图 (1.4) 的情况，
http://www.ma-xy.com 13 http://www.ma-xy.com
http://www.ma-xy.com
1.2 太阳影子定位研究 第一章 太阳影子定位
杆 h 的影长公式为
l =
h
tan ρ
太阳光在地面上的投影向量为
l = m − (m ·n)n
其长度平方为
|
l|
2
= |m − (m ·n)n|
2
= |m|
2
− |m ·n|
2
= 1 − sin
2
ρ = cos
2
ρ
O 点的正北方向为 O 指向北极的方向在 O 点出切平面的单位投影向量，也就是 z 轴正想在切
平面上的单位投影向量，即
N =
z − (z ·n)n
|z − (z · n)n|
由于
z −(z ·n)n = (−sin φ cos φ cos w, −sin φ cos φ sin w, cos
2
φ)
其长度平方为
|z|
2
− |z ·n|
2
= 1 − sin
2
φ = cos
2
φ
因此
N = (−cos w sin φ, −sin w sin φ, cos φ)
于是，正北方向 N 旋转至太阳影子方向的角度 (影角)γ 满足
N ·
l = cos ρ cos γ = cos δ cos w sin φ − sin δ cos φ
即
cos γ =
cos δ cos w sin φ − sin δ cos φ
cos ρ
注意，这里定义的影角与地学中定义的太阳方位角互为补角：余弦相差一个负号，正弦相同。结
合 sin δ 的计算公式，上式可以写为
cos
γ
=
sin ρ sin φ − sin δ
cos ρ cos φ
另一方面，地球的正东方向为
E = N × n = (−sin w, cos w, 0)
http://www.ma-xy.com 14 http://www.ma-xy.com
http://www.ma-xy.com
第一章 太阳影子定位 1.2 太阳影子定位研究
因而
N ·
l = |E||
l|cos(
π
2
− γ) = cos ρ sin γ = E · m = sin w cos δ
即
sin γ =
sin w cos δ
cos ρ
同时给出 sin γ 和 cos γ 是为了避免出现影角 γ 多解的情况。
和问题二中的处理方法相似，考虑影角相对变化以抵消 x 轴方向，建立影角的最小二乘目标
J(n, φ, θ) =
k
i=1
[(γ
i+1
− γ
i
) − (ˆγ
i−1
− ˆγ
i
)]
2
将上述目标和影长的最小二乘目标相结合，设置目标权重 c，有
min
n,φ,θ
J
(
n, φ, θ
) =
k
i=1
l
i+1
l
i
−
ˆ
l
i+1
ˆ
l
i
2
+ c[(γ
i
+1
− γ
i
) − (ˆγ
i
−
1
− ˆγ
i
)]
2
s.t.
− 90
◦
< φ < 90
◦
ρ
i
> 0
i = 1, 2, . . . , k
程序
结果
http://www.ma-xy.com 15 http://www.ma-xy.com
