全国大学生数学建模竞赛国家一等奖论文 A题 太阳影子定位模型的分析.docx

全国大学生数学建模竞赛国家一等奖论文 A题 太阳影子定位模型的分析.docx

《全国大学生数学建模竞赛国家一等奖论文 A题 太阳影子定位模型的分析.docx》由会员分享，可在线阅读，更多相关《全国大学生数学建模竞赛国家一等奖论文 A题 太阳影子定位模型的分析.docx（44页珍藏版）》请在冰点文库上搜索。

全国大学生数学建模竞赛国家一等奖论文A题太阳影子定位模型的分析

太阳影子定位模型

摘要

针对太阳影子定位问题，本文结合地理学和天文学的相关知识，建立了不同数据类型下的太阳影子定位模型，实现了视频拍摄地点和日期的快速精准确定。

对于问题一，首先从地理学角度，基于地理坐标，直杆长度，时间这三个影响影子长度的参数，计算出时角，赤纬角，太阳高度角，进而给出了影子长度与三个参数之间的关系式。

结果显示，影长对日期和时刻都呈现出先减小后增大的趋势；对杆长呈正比关系增长；对经度呈现先急剧增长到峰值再突变为0，而后

突变到峰值后再急剧下降；对纬度呈缓慢上升趋势。

然后，根据附件1中提供的数据，画出了天安门广场上直杆的太阳影子分布曲线图。

对于问题二，基于问题一中对影响影子长度因素的分析，根据地理学知识建

2021

∑ii

∑归i归i

立双目标规划模型，确立目标函数分别为：

min|∆A-∆A'|,min|S

-S'|。

i=1

i=1

然后在约束条件下对杆子的地点坐标应用网格逼近算法优化求解，得出最符合题目所提供数据的杆子地理位置为：

（19.1︒E,108.71︒N）——海南东方市境内，此时，杆长为2.03米，太阳方向角残差比为1.8%，影长残差比为0.9%，误差均很小。

对于问题三，首先建立了与问题二相似的目标规划模型，由于日期未知，模型求解的时间复杂度较高。

为提高计算速度，引入了粒子群算法。

分别对附件2和3中的数据进行分析，确定出的地点坐标分别为（80.51︒E，32.13︒N），（110.20︒E，24.83︒N）和（81.43︒E，32.24︒N），（111.56︒E，23.68︒N），附件2为西藏阿

里，日期为8/14或4/29，附件3为广西梧州市，日期为12/27或12/14。

可以发现，两种算法的结果极为接近，但粒子群算法计算时间要远小于网格逼近算法。

对于问题四，首先对视频数据进行采集和预处理，由于视频拍摄角度的存在，从视频中直接得到的影长并不是实际长度，而是其投影长度，这里采用基于Hough变换和透视变换的图像矫正法，对斜视图像进行矫正，得出实际影长。

然后将得到的数据带入问题二的模型中，给出视频拍摄地点为（110.70︒E，42.31︒N）

——内蒙包头市境内；在拍摄日期未知的情况下，将变化而来的实际影长代入问题三的基于粒子群算法的目标规划模型，求解出视频拍摄地点为（109.76︒E，42.66︒N）——内蒙包头市境内，拍摄日期为6/11或7/13。

对于模型的推广，根据物体采集到的太阳地理信息进行计算，可以应用到求建筑物群合理间距问题。

关键词：

双目标规划粒子群算法Hough变换透视变换

一、问题的重述

1.1问题的背景

“日长影移”是生活中人人熟知的自然现象，这个词说明地面上的影子变化与太阳活动有着密切的联系。

而古代智慧的先民就利用了这个现象制作了日晷，是最早且最精确的计时工具之一。

在图像信息充斥的当代，如何通过图像数据获得图像拍摄时的相关信息是图像分析学科的重要课题，而利用太阳光影变换获得时间和地理信息，是非常方便可靠的。

1.2问题的提出

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

针对上述背景和应用需求，提出以下问题：

1.建立影子长度变化的数学模型，分析影子长度关于各个参数的变化规律，并应用你们建立的模型画出2015年10月22日北京时间9:

00-15:

00之间天安

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

2.根据某固定直杆在水平地面上的太阳影子顶点坐标数据，建立数学模型确定直杆所处的地点。

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

3.根据某固定直杆在水平地面上的太阳影子顶点坐标数据，建立数学模型确定直杆所处的地点和日期。

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

4．如果已有一根直杆在太阳下的影子变化的视频，并且已通过某种方式估计出直杆的高度为2米。

请建立确定视频拍摄地点的数学模型，并应用你们的模型给出若干个可能的拍摄地点。

若拍摄日期未知，能否根据视频确定出拍摄地点与日期。

二、问题的分析

2.1问题一的分析

题目要求在固定地点，给定日期和杆长的条件下，求解出直杆投影长度的变化曲线。

对于水平地面上的垂直直杆，直杆长度与影子的比值即为太阳高度角的正切值，因此需要知道此时间段内的太阳高度角变化。

查阅资料[6]可得，太阳高度角与当地地方时、经纬度密切相关，因此知道上述两个量就可确定直杆的变化过程。

2.2问题二的分析

题目要求根据影子的变化情况和给出的日期求出直杆的位置，实际上是问题

一的逆求解过程。

这里杆长和地点都是未知量，逆求解是非常困难的，于是将问题二转化为双目标规划问题。

当太阳方位角与影长的实际值与理论值差值的绝对值之和达到最小时，所得经纬度即为杆子的地点坐标。

2.3问题三的分析

与问题二不同的是，该问中日期是个未知量。

首先考虑沿用上一问的模型。

由于日期未知，所以要考虑一年365天的所有情况，这将大大增加运算时间。

从减少运算量的角度考虑，有必要改进算法。

考虑引入现代优化算法之一的粒子群算法，将所有解视为粒子所要去的位置，由于适应度与目标函数相联系，选取合理的适应度函数，期望提高计算效率。

2.4问题四的分析

本问的数据由视频给出，那么首先要对视频进行数据预处理。

由于视频时间较长，所以不考虑将视频逐帧分析，而是每隔一段时间对视频进行获取分析，对图像所给出的影长和影子角度进行测量。

考虑到视频拍摄角度的存在，从视频中直接得到的影长并不是实际长度，而是其投影长度，因而要对斜视图像进行矫正，得出实际影长。

然后再分别利用问题二和问题三的模型进行求解。

三、问题的假设

1.每年的太阳活动情况是相同的，均为“恒星年”。

2.地球是一个完美的球形，不考虑海拔、地球扁率的影响。

3.无光线衍射造成的影子减淡现象。

5.在小尺度考虑直杆投影问题时，地表是绝对水平的。

6.不考虑地球公转的影响。

7.题目所给的数据是真实的，可靠的。

四、符号说明

符号

说明

θ

太阳高度角

L

水平地面上直杆长度

S

水平地面上直杆影子长度

h

当地地方时时角

δ

当地太阳赤纬

当地纬度值

γ

当地经度值

t

当地地方时

t0

北京时间

n

当日日期序号

P

当地经纬度坐标

A

当地太阳高度角

T

时间（包含时刻与日期）

β

离散化后的方位朝向

εi

每组数据影长与x坐标轴的夹角

xi，yi

题目附件提供的第i组直杆坐标

d

日期

注：

其他符号将在下文中给出具体说明。

五、模型的建立与求解

5.1问题一

5.1.1模型的建立

本模型结合相关地理学知识，对影子的变化情况进行分析描述。

下面将明确一些地理学定义，以及重新定义一些本模型需要用到的参数。

太阳高度角，也称太阳高度，是指某地的太阳光线与当地地平面的所交的最小线面角，这是以太阳视盘面的几何中心和理想地平线所夹的角度。

在水平地面上，直杆长度与影长的比值即为太阳高度角的的正切值：

tanθ=L

（1）

S

图1太阳方位角示意图其中L为杆长，S为影长。

通过查阅相关资料[6]，得知太阳高度角的计算公式为：

sinθ=coshcosδcosφ+sinδsinφ

（2）其中θ为太阳高度角，h为地方时时角，δ为当时的太阳赤纬，φ为当地纬

度。

以一个地方太阳升到最高的地方的时间为正午12时，将连续两个正午12时

之间等分为24个小时，所成的时间系统，称为地方时。

地球上每一个地点都有其相应的地方时。

由于题目只提供了当地时间的北京时间，因此在计算地方时时角时，要先将北京时间换算为当地地方时t：

γ-120︒

t=t0+

15︒

（3）

其中t0为北京时间，γ为当地经度。

根据某地地方时，可以换算出当地的地方时时角。

地方时时角h即为当地与

子午线之间相差的角度：

h=15︒⨯（t-12）（4）

太阳的赤纬等于太阳入射光与地球赤道之间的角度，由于地球自转轴与公转平面之间的角度基本不变，因此太阳的赤纬随季节不同而有周期性变化。

太阳赤纬的最高度数为23︒26'，夏至时太阳的赤纬为+23︒26'，冬至时太阳的赤纬为

-23︒26'。

春分和秋分时太阳的赤纬为0︒。

由于地球公转轨道的偏心率非常低，可以看作是一个圆圈，太阳赤纬δ可用下面这个公式来计算：

δ=23.45︒sin⎛2π（284+n）⎫（5）

ç365⎪

⎝⎭

其中n为当日日期序号，1月1日时n=1，以此类推得10月22日n=295。

联立式子

（2）-（5）得到方程组：

⎧L

tanθ=

S

⎪

⎪

⎪sinθ=coshcosδcosφ+sinδsinφ

⎪⎪

⎨t=t0+

⎪

γ-120︒

15︒

（6）

⎪h=15︒⨯（t-12）

⎪⎛+⎫

⎪δ=23.45︒sinç2π（284n）⎪

⎪⎩

求解上述方程组，得：

⎝365⎭

S=L/tan（arcsin（cos（15︒⨯（t-12））cos（δ）cosφ+sin（δ）sinφ））（7）可见，影子长度的变化与当地地理位置P（ϕ,φ），直杆长度L，时间T（t,d）

这三个参数有关。

5.1.2模型的求解

首先计算题目所给条件下的h，δ与φ，再将上述参数值代入

（2）式，得出从9:

00-15:

00的太阳高度角随时间的变化（具体值见附录）。

相应地，可求得直杆影长数据（具体见附录）。

从结果中挑选出几个比较重要的时间点，将相应结果制成下表以供参考：

表1影子长度分布

北京时间

太阳高度角/度

影长/m

9:

00

21.18

7.74

12:

00

37.88

3.85

12:

14

37.99

3.84

15:

00

25.35

6.33

将影长随时间变化的情况用MATLAB绘制成图像：

图2影子长度变化曲线

对于影子长度关于各个参数的变化规律，从时间，地点和杆长三个方面对其

进行分析，分析某一个参数时，将其他变量看做常量，只改变一个未知量。

我们将时间参数分为日期和时刻这两种情况进行运算求解，影长关于时刻的变化在图2中已经给出；对于地点，将其分为经度γ和纬度φ这两个参数进行运算求解，求解后的变化规律图如下：

图3影长与日期关系曲线图4影长与杆长关系曲线

图5影长与经度关系曲线图6影长与纬度关系曲线

5.1.3结果分析

从图3中可以发现，影长随着时间的增加，呈现先减小后增大的趋势，影长最小点出现在12:

14，这是由于北京时刻为120︒E的地方时，换算到116︒E附近时，会产生时差，显然是符合常理的。

从图4中可以发现，影长与杆长呈正比关系，这是由于tanθ=L，从而验证

S

了模型的准确性。

从图5中可以发现，影长随时间变化出现了两个明显的“脉冲”。

这是因为当太阳逼近地平线时，影长变化的速度非常快，影子也是最长的。

而太阳一旦在地平线以下，就不存在影子，影子长度也为0。

从图6中可以发现，图5中的曲线随纬度增加，总趋势是增加的，这是随着纬度的增加太阳高度角减小，从而导致了影长的增加。

5.2问题二——基于网格逼近算法的双目标规划模型

5.2.1数据预处理

分析题目所给的数据（附件1），可以发现这些数据不仅能计算出影长，还能够得知影子方位角的变化，但是附件没有给出坐标系的方位朝向，故假设坐标系x轴的方位朝向角为β（0︒≤β<360︒），取正北方向时，x轴的方位朝向角为0︒。

由此，可以得出21个时间点影子的太阳方位角。

我们以每五个时间点为间隔，选取部分呈现在下表中：

表3影子方位角数据表

北京时间

x坐标/米

y坐标/米

方位角/度

14:

42

1.0356

0.4973

β+25.63︒

14:

57

1.2087

0.5255

β+23.49︒

15:

12

1.3955

0.5541

β+21.65︒

15:

27

1.6033

0.5833

β+20.02︒

15:

42

1.8277

0.6135

β+18.55︒

5.2.2模型的建立

首先引入太阳方位角的定义。

太阳方位角是太阳在方位上的角度，它通常被定义为从北方沿着地平线顺时针量度的角：

cosA=sinδ⋅cosφ-cosh⋅cosδ⋅sinφ

cosα

（8）

上述公式可以用来计算近似的太阳方位角，不过因为公式是使用余弦函数，所以方位角永远是正值，因此，角度永远被解释为小于180度，而必须依据时角

来修正。

当时角为负值时（上午），方位角的角度小于180度，时角为正值时（下

午），方位角应该大于180度，即要取补角的值，故作如下修正：

⎧arccos⎛sinδ⋅cosφ-cosh⋅cosδ⋅sinφ⎫,

h<0

⎪çcosα⎪

⎨

A=⎪⎝⎭

⎪360︒-arccos⎛sinδ⋅cosφ-cosh⋅cosδ⋅sinφ⎫,h≥0

（9）

⎪çcosα⎪

⎩⎝⎭

由上述推导可得，方位角A与坐标系方位朝向β的关系为：

Ai=β+εi,

i=1,2,,21（10）

其中εi为每组数据影长与x坐标轴的夹角：

i

ε=arctanyi,

i=1,2,,21（11）

xi

将上一问题的某些参数关系与式（8）-（11）列成方程组：

⎧

⎪εi

⎪

=arctanyi,

xi

i=1,2,,21

⎪Ai=β+εi,

⎪

i=1,2,,21

⎛+⎫

⎪ç

⎪δ=23.45︒sin

2π（284n）

365⎪

（12）

⎨⎝⎭

⎪h=15︒⨯（t-12）

⎪

⎪θ=cos（15︒⨯（t-12））cos（δ）cosφ+sin（δ）sinφ

⎪

⎪t=t0+

⎪γ-120︒

⎩15︒

将上述求太阳高度角的方程化为如下关系式：

cos（A）=

（sin（δ）⋅cosφ-cos⎡︒⨯（t+

⎣

γ-120︒

15︒

-12）⎤⨯cos（δ）⋅sinφ）

⎥⎦

（13）

⎛⎛γ-120︒⎫⎫

cos（arcsin（cosç15︒⨯çt0+-12⎪⎪cos（δ）cosφ+sin（δ）sinφ））

⎝⎝15︒⎭⎭

下面对各个参量进行影响因素的分析：

γ-120︒

对于时角来说，由t=t0+

15︒

，h=15︒⨯（t-12）式子可知，时角的变化与

⎛2π（284+n）⎫

当地经度γ的变化以及时间t0有关；对于赤纬角，由δ=23.45︒sinç⎪

⎝365⎭

式子可知，赤纬角的变化与当地日期有关，而在此模型中，日期为定值，故可认为赤纬角δ为常量；对于太阳高度角，由sinθ=coshcosδcosφ+sinδsinφ式子可知，太阳高度角的变化，与时角的大小，纬度的大小有关，即与经纬度的大小有关，我们可以将其理解为地点的变化；对于太阳方位角，由式子（8）可知，太阳方位角的变化，与时角的大小，纬度的大小以及太阳高度角有关，即与经纬度的大小有关，同样可以将其理解为地点参数的变化。

综上，我们对参量的影响因素进行总结，可以得到下表：

表4影响因素分析

参量

影响因素

时角（h）

经度γ，时间t0

赤纬角（δ）

常量

太阳高度角（θ）

地点坐标（γ,φ），时间t0

太阳方位角（A）

地点坐标（γ,φ），时间t0

可以发现，在式（13）中，方程右端的未知量只有地点坐标（γ,φ）和北京时间

t0，故可以将其简化为如下形式：

F（γ,φ,t0）=cos（A）（14）

式中，映射F表示方程右端关于地点坐标（γ,φ）这两个自变量的函数。

因此，对于21个不同时刻对应的太阳方位角Ai（i=1,2,...,21），我们可以得出

21个不同的方程，将其列成方程组如下所示：

⎧F（γ,φ,t0）=cos（A1）

⎪

⎪F（γ,φ,t0）=cos（A2）

⎨

⎪

⎪⎩F（γ,φ,t0）=cos（A21）

（15）

由式（9），我们可以根据地点坐标（γ,φ）计算得出此处不同时刻的太阳方位角

Ai，将其去前一项作差，即可得到太阳方位角的差值∆Ai；由式（7），我们可以根

''

i

据地点坐标（γ,φ）算出此处不同时刻下的影子长度S'，对其进行归一化处理，可

以得到：

Si-

21

∑Si

'

归iS'

21i=1

-S'

（16）

imax

imin

i

对此我们可采用遍历算法对不同经纬度（γ,φ）下的太阳方位角的插值∆A'以

及归一化后的影子长度S归i

优解。

进行求解，从而与实际值作差进行比较，进而得出最

故此，我们建立多目标规划模型：

⎧

20

⎪∑ii

⎨

目标函数：

⎪

min|∆A-∆A'|

i=1

21

（17）

⎪

⎪min|S

-S'|

∑归i归i

⎩i=1

⎧-180︒≤γ≤180︒

⎨

约束条件：

⎪-90︒≤φ≤90︒

⎩

⎪0︒≤θ≤90︒

（18）

式中，∆Ai表示实际的太阳方位角，S归i表示归一化后的实际影长。

5.2.1模型的求解

为了方便计算机编程求解，我们对经纬度进行离散化处理，选取步长为0.01，对其进行MATLAB编程求解：

表5定点坐标数据

经度

纬度

杆高

108.71︒

19.12︒

2.03m

表6误差分析表

太阳方向角残差

影长残差

太阳方向角残差比

影长残差比

0.026

0.014

1.8%

0.9%

即杆子所处地点坐标为（19.1︒E,108.71︒N），在海南东方市境内。

5.2.1结果分析

通过误差分析表可以发现，在经纬度坐标为（108.71︒E,19.12︒N）时，计算所得的影子长度残差比为0.9%，太阳方向角残差比为1.8%，这个数值是较小的，说明得到的经纬度坐标（108.71︒E,19.12︒N）是较为精确的。

5.3问题三

5.3.1模型的建立

不难发现，此问在上一问的基础上将日期变为了未知量，在其他输入参数不变的情况下，要求模型输出杆子所处的地点和日期。

由表3可以看出，日期影响的参量为赤纬角，其他参量并不受日期的影响，太阳赤纬公式为：

δ=23.45︒sin⎛2π（284+n）⎫（19）

ç365⎪

⎝⎭

式中，n为当日日期序号，1月1日时n=1，以此类推，日期每增加一天，对n进行加1即可。

对此，我们采用遍历算法，对日期参数n及经纬度坐标（γ,φ）进行遍历求解，建立与上一问相似的多目标规划模型，从而求出最优解：

⎧

20

⎪∑ii

⎨

目标函数：

⎪

min|∆A-∆A'|

i=1

21

（20）

⎪

⎪min|S

-S'|

∑归i归i

⎩i=1

⎧-180︒≤γ≤180︒

⎪-90︒≤φ≤90︒

⎨0︒≤θ≤90︒

⎪⎩L>0

（21）

式中，∆Ai表示实际的太阳方位角，S归i表示归一化后的实际影长。

5.3.2模型的求解

同样，对经纬度进行离散化处理后，选取步长为0.01，对其进行MATLAB编程求解：

表7杆子所处地点和日期及误差分析数据表

附件二

附件三

经度

80.51︒E

经度

110.20︒E

纬度

32.13︒N

纬度

24.83︒N

杆高

2.04m

杆高

3.10m

太阳方向角残差比

2.0%

太阳方向角残差比

2.3%

影长残差比

1.1%

影长残差比

0.35%

日期

8/14或4/29

日期

12/27或12/14

5.3.3模型的改进——基于粒子群算法的目标规划模型

在实际求解中，发现由于加入新的未知量日期，使得模型的可行域大大增加，增加了遍历算法的时间复杂度，因此，我们引入粒子群算法[9]对目标函数进行求解，从而降低模型求解的时间复杂度。

假定有一个D维的目标搜索空间，有n个微粒组成了一个粒子群，其中每个微粒都用一个D维的向量描述，将它的空间位置表示为mi=（mi1,mi2,,miD）,i=1,2,,n；这可看做目标优化问题中的一个解，代入适应度函数计算出适应度值可以衡量微粒的优劣；第i个微粒的飞行速度也是一个D维的向量，记为vi=（vi1,vi2,,viD）；第i个微粒所经历过的具有最好适应值的位置称为个体历史最好位置，记为pi=（pi1,pi2,,piD）；整个