2015数学建模A题太阳影子定位论文

太阳影子定位

摘 要

本文主要讨论分析并建立太阳影子定位的数学模型，利用所建模型及所给数据确定事件发生的位置，运用多种数学方法研究物体影子长度与当地经纬度、时间及日期之间的联系，具有一定的实际意义。主要利用MATLAB 对数据进行分析和处理。

针对问题一：利用日期求出当天太阳直射点的纬度，然后利用太阳直射点的纬度，当地的纬度及时角求出太阳高度角，由三角函数关系得直杆影子长度，利用MATLAB 对所求数据曲线拟合构建影子长度变化模型，分析影子长度关于各个参数的变化规律。利用所求模型求出问题一中直杆的影子长度及其变化曲线见图5.1.3.3。

针对问题二：建立时间与影子长度之间的模型，利用MATLAB 曲线拟合并对求出的曲线函数分析，由时间差得出当地的经度。方法一根据相似三角形求出直杆高度，通过三角函数关系以及求太阳直射点纬度模型得直杆所在地的经纬度

()''''''1824107,4670??东经北纬或()'

'''''1824107,484420??东经北纬。方法二利用太阳方位角和太阳高度角模型求出直杆所在地的经纬度()

''''''1824107,1628??东经南纬。

针对问题三：利用相似三角形求出直杆的长度，通过曲线拟合求出正午最短的影子长度，根据三角函数关系求出正午太阳高度角和太阳高度角的方程计算当地纬度和太阳直射点的纬度，再利用问题二构建的模型求经度，最后求出附件二直杆所在地的经纬度为（北纬''''''221668,485239??东经），日期为5月24日，附件三直杆所在地的经纬度为（南纬'''151051?，东经'''2519105?），日期为3月24日。

针对问题四：利用MATLAB 对视频进行处理，求出任意时刻影子长度，通过时间差求出视频拍摄位置的经度。利用问题一求出太阳直射点的纬度和太阳高度角随时间的变化，继而求出直杆所在地可能的经纬度），东经（北纬’''2923218626

.38??或()'

''29231213626.38??，东经北纬。通过曲线拟合出当地正午时最短的影子长度，利用三角函数求出正午太阳高度角和求太阳高度角的方程，计算当地纬度、太阳直射点的经纬度（南纬'''151051?，东经'''2519105?），日期为3月24日。

关键词：太阳高度角 太阳方位角 太阳直射点 曲线拟合 相似三角形 MATLAB

一问题重述

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

请分析题目，试建立数学模型讨论下列问题：

1.建立影子长度变化的数学模型，分析影子长度关于各个参数的变化规律，并应用你们建立的模型画出2015年10月22日北京时间9:00-15:00之间天安门广场（北纬39度54分26秒,东经116度23分29秒）3米高的直杆的太阳影子长度的变化曲线。

2.根据某固定直杆在水平地面上的太阳影子顶点坐标数据，建立数学模型确定直杆所处的地点。将你们的模型应用于附件1的影子顶点坐标数据，给出若干个可能的地点。

3.根据某固定直杆在水平地面上的太阳影子顶点坐标数据，建立数学模型确定直杆所处的地点和日期。将你们的模型分别应用于附件2和附件3的影子顶点坐标数据，给出若干个可能的地点与日期。

4.附件4为一根直杆在太阳下的影子变化的视频，并且已通过某种方式估计出直杆的高度为2米。请建立确定视频拍摄地点的数学模型，并应用你们的模型给出若干个可能的拍摄地点。

如果拍摄日期未知，你能否根据视频确定出拍摄地点与日期？

二问题分析

问题一中，由地理知识可知直杆高度、影子长度与太阳高度角存在三角函数关系，因此要求出太阳高度角随时刻的变化就能求出影子长度的变化。而太阳高度角的变化与太阳直射纬度、当地纬度和时角有关，因此要先求出该日期的太阳直射纬度。依次求解建立模型，并利用MATLAB进行数据拟合求出直杆影子变化曲线。最后利用建立的模型求出直杆的影子长度及变化曲线。

问题二中，法一利用附件1所给的影子顶点坐标位置以及北京时间和日期，建立影子长度L与时间M之间的数学模型，运用MATLAB对影子长度和北京时间进行曲线拟合并对拟合出来的曲线进行求导，得出当地正午太阳高度时直杆的最短影子长度及此时的北京时间。利用地理知识根据北京时间与当地正午十二点之间的时间差值得出当地的经度。接下来求纬度，利用相似三角形求出直杆的高度为2.76米，利用影长与直杆之间的三角函数关系求出正午太阳高度角H，根据问题一中的模型求太阳直射纬度，利用正午太阳高度角计算公式计算当地纬度,从而得出可能的经纬度位置。法二中利用太阳高度

角和太阳方位角避开直杆高度，并利用问题一求太阳高度角的模型求出所在地的纬度。

问题三中利用问题二中求经度的模型来对附件2、3中所给的影子长度和时间对所在位置的经度进行求解，利用相似三角形求出直杆的高度。将附件中的数据使用MATLAB进行曲线拟合得出当地正午时最短的影子长度，根据三角函数关系求出某一时刻的太阳高度角，利用问题一中建立的太阳高度角模型和正午太阳高度角计算公式求出当地的纬度和太阳直射点的纬度。利用问题一求太阳直射点的纬度模型求当天的日期。

问题四中，第一小问，首先利用视频中所给的日期，利用问题一建立的模型求出太阳直射点的纬度，再利用MATLAB对所给视频进行灰度处理，求出影长，利用三角函数关系，将影长数据代入求出太阳高度角随时间的变化，最后利用问题一中建立的太阳高度角模型的公式求得当地的纬度。第二小问，用MATLAB进行曲线拟合得出当地正午时最短的影子长度，根据三角函数关系求出某一时刻的太阳高度角，利用问题一中建立的太阳高度角模型和正午太阳高度角计算公式求出当地的纬度和太阳直射点的纬度，并通过问题一求太阳直射点的纬度模型求当天的日期。另外本文以北京时间为标准，利用当地与北京的时间差确定视频拍摄的经度。

三、模型假设

1.假设模型中一天中太阳直射纬度的变化忽略不计；

2.假设不考虑固定直杆的海拔高度；

3.假设地球自转时全球各地线速度一致；

四、符号说明

h太阳高度角；

H正午太阳高度角；

l直杆杆长；

L直杆影子的长度；

?地理纬度；

t时角；

δ太阳直射点的纬度；

A太阳方位角。

五、模型的建立与求解

5.1 问题一：模型的建立与求解

5.1.1太阳直射点纬度δ[1]

如图5.1.1.1,O 为地球,S 为太阳,δ为太阳直射点纬度,SOC ∠=δ;λ是黄经度(太阳在黄道上由春分点γ自西向东运行到S 点所转过的角度,即SOA ∠);γO 是黄道平面和天赤道平面的交线(棱),OS AS OD 、、在黄道面内,OC AC OE O AS O OD 、、;,γγ⊥⊥分别是OS AS OD 、、在天赤道面的射影,它们的垂足分别是;C C E 、、AOE ∠是黄赤交角(二面角的平面角),'''212623?=∠=∠SAC DOE 。在SAO Rt ?中,λsin OS AS =;在OCS Rt ?中，δsin ?=OS SC ，在ACS Rt ?中，''''''212623sin sin 212623sin ???=??=λOS AS SC '''232623sin sin sin ???=?λδOS OS ，所以，'''232623sin sin sin ??=λδ，得到太阳直射点纬度δ的计算公式：

()λδsin 397775.0arcsin =

图5.1.1.1

地球公转的周期是一个回归年（365.2422），现行公历的历年是历日的整数倍，它和回归年并不精准相等；另外由于复杂的历史演变过程，以及一些人为的原因，上半年和下半年、冬半年和夏半年以及各季、各月之间的天数也并不完全相同，因此用日期来推算太阳直射点纬度(天文上称赤纬,一年中变动在'2623?±范围内)时，需按季节分段时行计算，以确保推算的相对精确。

因此，把全年的日期分成夏半年、冬半年以及自冬至日到春分日三个阶段来进行推

导太阳直射点纬度δ。

)1.夏半年：

从春分日（3月21日前后）到秋分日（9月23日前后）约186天，在此时段内视太阳在黄道上运转?180。设从春分日开始，视太阳运行了d 天，则d 天运行了β经度：

d ??=186

180β 把上式带入()βδsin 397775.0arcsin =中得：

?????

???? ?????=d 186180s i n 39775.0a r c s i n δ []()1862,1,0 ∈d ； )2.冬半年：

自秋分日到冬至日（12月22日）前后总计90天，在此阶段上运行了?90，设从春分开始，视太阳运行了d 天，则d 天运行了β经度：

()18690

90180-??+?=d β 把上式带入()βδsin 39775.0arcsin =中得：

()[]?-?-=186sin 39775.0arcsin d δ []()276

187,186 ∈d ； )3.自冬至日到次年春分日总计89天，在此阶段内也假设运转了?90。设从春分日开始，假设太阳运转了d 天，则d 天运行了β经度：

()27690

90270-??+?=n β 把上式带入()βδsin 39775.0arcsin =中得：

()?

???

????????-???-=2768990cos 39775.0arcsin n δ []()365277,276 ∈n ； 所以：

=δ []()[][]()[]

?????

??????∈????????????-???-∈?-?-∈????????? ?????365277,2762768990cos 39775.0arcsin 276188,187,186186sin 39775.0arcsin 1863,2,1,0186180sin 39775.0arcsin d n d d d d

5.1.2太阳高度角h [2] 我们知道中午太阳高度角最高;在北半球夏季比冬季的太阳高度角要高;低纬度地区要比中、高纬度地区的太阳高度角高。说明太阳高度角的变化是随时间、地理纬度和太阳直射点而变化的。太阳高度角计算公式是一个多元函数方程。首先研究真太阳时(真太阳时=地方平均太阳时+时差)正午的太阳高度角计算,见图5.1.2.1,AGFED 表示以M 为测站中心的南北天子午圈,ME 为赤道截面,F 为天顶。

图5.1.2.1

在真太阳时正午时,太阳正位于当地子午面上(即O ), δ?+-=90h ，由此式可见,在固定的测站上,?为常数,真太阳时正午的太阳高度角只决定于太阳倾角的变化,在一年内真太阳时正午的太阳高度角变化在'2623?±之内。

计算任意时刻的太阳高度角,见图5.1.2.2,因为是任意时刻,太阳不恰好处在当地子午面上,O 为太阳移动的位置,AGFED 为天子午圈,BE 为赤道截面,BF 为天顶距,LG 为赤道天顶距。赤道截面与天顶的交角为纬度角,故有FE ?=,则GF ?-=90。太阳与赤道截面的交角为太阳倾角则OL=δ,故OG=δ-90。太阳与地平面的交角为太阳高度角,则OB=h, FO=90-h 。在GOF 球面三角形中,根据球面三角公式：

A c b c b a cos sin sin cos cos cos +=

则有：

()()()()()t h cos 90sin 90sin 90cos 90cos 90cos δ?δ?--+--=-

将此式简化得出：

t cos cos cos sin sin sinh δ?δ?+=

t 是太阳位置与当地子午面的偏角,即时角(所谓时角,是指太阳所在的时圈与通过南点的时圈构成的夹角,单位为度。自天球北极看,顺时针方向为正,逆时针方向为负。时角表示太阳的方位,因为天球在一天24h 内旋转360°,所以每小时旋转15°，当地时间12点时的时角为零）。

图5.1.2.2

5.1.3求影子长度变化

如图5.1.2.2，利用三角函数，根据L

l tanh =即可求出指定时刻直杆的影子长度。

图5.1.2.2影子长度与太阳高度角的关系

由上述求解过程分析得出影子长度的变化规律：

1) 与经纬度有关：

就某一天来看，太阳直射纬度所在的位置，日影最短，为一圆点。在直射纬度以北、以南的地区，正午日影随着正午太阳高度缩小而逐渐变长。

2) 与时间有关：

就某个地点来看，一年中正午太阳高度增大时，日影逐渐缩短；正午太阳高度达最

大时，日影最短；正午太阳高度减小时，日影逐渐增长；正午太阳高度达最小时，日影最长。

例如：①6月22日，太阳直射北回归线，北回归线及其以北各地的正午太阳高度达到全年最大，其日影也达到全年最短；

②6月22日~12月22日，在太阳直射点向南移动过程中，北回归线及其以北各地的正午太阳高度逐渐减小，那么其日影逐渐增长；

③12月22日，太阳直射南回归线，北回归线及其以北各地的正午太阳高度达到全年最小，其日影也达到全年最长。

④12月22日~6月22日，在太阳直射点向北移动过程中，北回归线及其以北各地的正午太阳高度逐渐增大，那么其日影逐渐缩短。

5.1.3利用模型求解

通过上述模型来对2015年10月22日北京时间9:00-15:00之间天安门广场（北纬39度54分26秒,东经116度23分29秒）3米高的直杆的太阳影子长度的变化进行求解。每过15分钟计算一下本时刻的太阳高度角，数据见表5.1.3.1。

表5.1.3.1太阳高度角随时刻的变化

利用MATLAB对表5.1.3.1中的数据进行拟合求出直杆的影子长度的曲线变化见图5.1.3.3，详细代码见附录1。

图5.1.3.3时刻与直杆影长之间的关系

5.2问题二：模型的建立与求解

5.2.1所在地经度的求解

本文首先利用附件一中所给的影子顶点坐标求出影子长度，建立横坐标为时间M与纵坐标为影子长度L之间的函数关系，利用MATLAB整合出曲线方程以及曲线拟合图像，如图5.2.1.1：

图5.2.1.1

曲线方程为：

=M

M

L

-

.02+

752

.3

1489

24013

正午太阳高度最大，直杆影子长度最短，由于图像横坐标为北京时间，因此曲线最低点对应的北京时间为当地正午十二时。对所求曲线方程进行求导，求出最低点的北京时间和最短的影子长度：

752.32978.0'-=M L

令0'=L ,得出M=12.5991，0.4952l =。

利用地理知识经度每相差?15在时间上相差1小时得出当地的经度与北京的经度相差?9865.8，由求导得出的结果可知北京时间比当地时间快，因此该地应在北京的西面，所以得出当地经度为：

?=-4049.1079865.83914.116

5.2.2所在地纬度的求解

法一：

由于直杆高度和地球最外层距离地面的距离都是固定的，利用直杆影长的变化和经过相同时间地球自转距离的差值，进而用相似三角形求出直杆的高度，原理如图

5.2.2.1。

图5.2.2.1

相似三角形的比值公式如下：

2

121P P X L L l -=- 21L L -为确定时间内影子变化长度；

2

1P P -为确定时间内地球自转距离的差值[4]； X 为地球最外层距离地球表面的距离[5] 。 根据比值求出固定直杆的高度m l 7226.260

346510306.06

=????=。

利用直角三角形的角度关系L l H =

tan 求出当地正午太阳高度角H ，代入正午最短的直杆影子长度，由L

H 3arctan =最终求得正午太阳高度角?=6914.79H 。 通过问题一中求太阳高度角的计算公式t cos cos cos sin sin sinh δ?δ?+=，将正午的时角0代入，利用三角函数和角公式B A B A B A os sin sin cos cos )(c ?+?=-可推导出正午太阳高度角计算公式纬差-?=90H ，纬差即为当地纬度与太阳直射点纬度之差，从而求出纬度位置，利用题目一所建立的日期与太阳直射点纬度之间的关系模型来求出4月18日太阳直射点纬度，最终求得当地纬度:

?=-?-?6914.794381.1090x

????=?=7467.201295.02

1x x 最后得出直杆两个可能的位置:

()''''''1824107,4670??东经北纬，()'

'''''1824107,484420??东经北纬 法二：

利用太阳高度角求出所在地的纬度[3]

太阳方位角即太阳所在的方位，指太阳光线在地平面上的投影与当地子午线的夹角，可近似地看作是竖立在地面上的直线在阳光下的阴影与正南方的夹角。方位角以目标物正北方向为零，顺时针方向逐渐变大，其取值范围是0--360°。因此太阳方位角一般是以目标物的北方向为起始方向，以太阳光的入射方向为终止方向，按顺时针方向所测量的角度，如图5.2.2.1：

图5.2.2.1

根据“天球理论”，可以导出太阳方位角与太阳直射点纬度、太阳高度角以及观测

时刻的关系如下:

cosh

sin cos sin t A δ= 由图5.2.2.1可以得出：

()()????

?????+=-=-?-=+=）（）（

2cos 1cos 90sin sin sin 22y x x A θθθπα 由上述方程组可知：

x

y x t -+=2

2sin cos cosh δ 又有t cos cos cos sin sin sinh δ?δ?+=，将附件一中的北京时间转化为当地的时角代入方程，并利用题目一所建立的日期与太阳直射点纬度之间的关系模型来求出4月18日太

阳直射点的纬度，最后由MATLAB 解得：?-=0378

.8? ，所以求出直杆所在位置为()'

'''''1824107,1628??东经南纬。

5.3问题三：模型的建立与求解

5.3.1所在地经度的求解

本文首先利用附件2和附件3中所给的影子顶点坐标求出影子长度，建立横坐标为时间M 与纵坐标为影子长度L 之间的函数关系，利用MATLAB 分别拟合出附件2、附件3数据的曲线方程以及图像，如图5.3.1.1和5.3.1.2，代码见附录2：

图5.3.1.1

图5.3.1.2

附件2拟合出来的曲线方程为：

32.23985.209814.012

11+-=M M L

附件3拟合出来的曲线方程为： 56.51551.72964.022

22+-=M M L

正午太阳高度最大，直杆影子长度最短，由于图像横坐标为北京时间，因此曲线最低点对应的北京时间为当地正午十二时。对上述所求曲线方程分别进行求导，求出最低点的北京时间和最短的影子长度： 985.219628.01'

1-=M L

551.75928.02'2-=M L

令0'1=L ,得出2079.151=M ，0.6223L 1=。

令0'2=L ,得出7379.122=M ， 3.4809L 2=。

利用地理知识经度每相差?15在时间上相差1小时得出当地的经度与北京的经度分别相差'''7748?和'''4411?，由求导得出的结果可知北京时间比当地时间快，因此两地都应在北京的西面，所以得出两地经度分别为：

'''''''''22166877482923116?=?-?

'''''''''251910544112923116?=?-? 5.3.2所在地纬度与日期的求解

对附件2，利用问题二的相似三角形模型求出固定直杆的长度，所以

m l 5.260

346510306.06

=????=。 利用直杆的高度与当地最短的影子长度的比值得出正午太阳高度角L H 5.2arctan =，求得?=02.76H 。取附件二中某一时刻的影子顶点的坐标并求出影子长度，运用三角函数公式求出该时刻的太阳高度角，联立方程组：

??

???±-?=+=?δδ?δ?90cos cos cos sin sin sinh H t 分两种情况讨论：若当地纬度与太阳直射点的纬度在同一半球，纬差为?δ-；若当地纬度与太阳直射点的纬度在不同半球，纬差为?δ+。用MATLAB 求解得当地纬度?=87.39?，根据问题一建立的求太阳直射点的纬度的模型将δ代入三个阶段，筛选出符合条件的日期为5月24日。当地的经纬度为（北纬''''''221668,485239??东经）。

同理，对附件三进行求解得，太阳直射点的纬度δ为?3.1430，当地经纬度为（南纬'''151051?，东经'''2519105?），日期为3月24日。

5.4问题四：模型的建立与求解

5.4.1模型假设

1) 假设直杆底座厚度忽略不计；

2) 假设忽略风速对直杆高度、影子长度的影响；

3) 假设拍摄过程中摄像头位置及角度的变化忽略不计。

5.4.2第一小问模型的建立与求解

1) 所在地纬度的求解

本文首先利用视频中所给的日期，利用问题一中建立的求太阳直射点纬度的模型进

行求解，根据公式?????

???? ?????=d 186180sin 39775.0arcsin δ求出太阳直射点的纬度'''364943?=δ。

对用MATLAB 进行灰度处理的21张图片进行影子长度的测量，测量数据见表

5.4.2.1，利用影子的长度变化数据根据三角函数求出对应时刻的太阳高度角，根据问题

一建立的求太阳高度角的数学公式t cos cos cos sin sin sinh δ?δ?+=，利用MATLAB [6] 对数据进行求解得出当地的纬度?=8626.38?。

表5.4.2.1 2) 所在地经度的求解

我们以北京时间为标准，求出北京与当地的时间差，从而根据两地的经度差求出当地的经度。假设将此2m 的直杆放在北京，利用达到相同影子长度的时刻求时间差，选取某一时刻的影子长度，求出该时刻的太阳高度角，再利用视频中已知的日期求出当天太阳直射点的纬度，利用问题一建立的模型t cos cos cos sin sin sinh δ?δ?+=求出时角?±=2.50t ，由两地的时角差计算出两地的时间差，继而根据经度每隔?15时间相差一小时可求出分别位于北京东面和西面的两地经度为''''''2923121292321??和。

因此，得出视频大体拍摄位置为），东经（北纬’''2923218626

.38??和()'

''29231213626.38??，东经北纬。 5.4.3问题4.2模型的建立与求解

利用MATLAB 对上述影长与时间进行曲线拟合，图像如图5.4.2.2得出的表达式如下：

16.16201.207359.02+-=M M L ；

图5.4.2.2

当地正午十二时时，直杆影子的长度最短，根据曲线求得最短影长L=0.3450。利用三角函数关系，求出影子最短时的正午太阳高度角?=2127.80H 。再联立方程组：

??

???±-?=+=?δδ?δ?90cos cos cos sin sin sinh H t 筛选出合适的?=?=9779.47651.14δ?，。将太阳直射点的纬度δ代入问题一建立的

模型中，得出日期为3月26日。

经度的求解同问题四第一小问中的求法，可得两地经度为''''''2923121292321??和。因此，视频中直杆的拍摄位置为（南纬'''544514?，东经'''292321?）和（南纬'''544514?，东经'''2923121?）。

综上所述，如果拍摄日期未知，可以根据视频确定出若干地点与日期，不同的一个日期能够确定一个不同的经纬度位置满足视频中时间与影子长度的对应关系。

六、模型的评价

6.1模型评价

6.1.1问题一模型评价

问题一建立的几何模型，运用三角函数的数学知识和地理知识，用简单的方法去解决复杂的问题，简明易懂。本模型主要解决了物体影子长度与所在地的经纬度、时间以及日期之间的关系，该模型基于严密的数学推导，求解过程严谨，结果可信度高，说服力强。

6.1.2问题二模型评价

问题二模型是基于问题一的进一步理论推导，模型理论严谨。

6.1.3问题三模型平价

问题三是在第二问的基础上进一步研究，频繁利用一二题中建立的模型，利用MATLAB等软件进行求解。

6.1.4问题四模型评价

问题四模型综合了问题二、三模型中的情况，利用前三问的结论进行求解。

6.2模型改进

由于利用相似三角形进行求解采用地球的自转与距离有所误差，对视频进行处理时因为摄像机角度原因和直杆底座厚度引起测量数据的不准，以及忽略了风速对影子长度的影响，需要改进。

七参考文献

[1] 蒋洪力.太阳直射点纬度的数学推导和分析.数学通报.2007

[2] 李玉海.太阳高度角及其计算.气象.1977.052.01

[3] 贺晓类.太阳方位角的公式求解及其应用.太阳能学报.2008年01期

[4] 百度百科.地球自转的速度. https://www.wendangku.net/doc/9d10887471.html,/q/1363977514068699

[5] 百度百科.地球最外表距离地球表面的距离. https://www.wendangku.net/doc/9d10887471.html,/doc/9525544-9869677.html

[6] 张志涌.MATLAB教程.北京航空航天大学出版社.2010.8.