太阳影子定位

太阳影子定位模型建立

摘要

本文讨论求解了在直杆影子随时间变化过程中，在知道日期、杆位置、影子坐标、时

间等参数条件中的某几个前提下，设计了确定型模型进行求解。

分析太阳方位与直杆影子关系，首先，将地球自转公转视为地球不动太阳动，利用立

体几何知识得出太阳高度角与影子长度关系。问题一的关键在于太阳高度角与日期、竿

位置、时间参数的关系。问题二中我们将立体平面化，把太阳与地球的运动关系转化为

平面上的角度关系，使模型简明直接。在模型求解时，我们把各解看为离散型随机变量，

对解进行权重处理，最后求得较精准的解。问题三，先结合前两题的模型预处理，再利

用matlab据最小二乘法原理，来对目标函数进行曲线拟合求解。对问题四中视频进行

分段截取照片处理，用photoshop软件测量影子长度与时间关系，再结合前几题模型与

求解方法，可求得结果。问题被函数化，模型简明直接，提高了确定性。

关键词：太阳高度角，立体平面化，权重处理，matlab曲线拟合

问题重述

确定视频的拍摄地点和拍摄日期是视频数据分析的重要方面，太阳影子定位技术就

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

通过影子长度变化建立数学模型，分析影子长度关于各个参数的变化规律，并应用

建立的模型画出某时间段某地某固定直杆的太阳影子长度的变化曲线。

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

处的地点。并利用模型对附件1的影子顶点坐标数据进行求解，求出若干个可能的地点。

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

处的地点和日期。将模型分别应用于附件的影子顶点坐标数据，求出若干个可能的地点

与日期。

根据一根直杆在太阳下的影子变化的视频，直杆的高度为2米。建立确定视频拍摄

地点的数学模型，并应用你们的模型给出若干个可能的拍摄地点。日期未知下再尝试求

解。

问题分析

根据影子变化来确定时间地点和时间，我们把地球自转看成太阳绕地球转，可以转

化为太阳方位与地球各地点和时间的关系问题。

对于问题一：可以把影子长度变化规律转化为光线与水平面夹角的变化规律。我们

根据地球自转公转规律和立体几何知识建立模型，且该模型得能体现光线和地面的夹角

与日期时间、地理位置的关系，最后通过matlab画出影子长度随时间的变化曲线。

对于问题二：问题二相比问题一缺少一个已知量，无法通过问题一中的模型来求解，

我们把太阳与直杆影子的关系转化为了平面角度关系，进而简便有效地求出杆的位置。

对于问题三：已知量较前一问更少，故我们先结合问题一和二建立的模型，再应用

matlab进行曲线拟合求得参数解。

对于问题四：通过提取视频特定帧，测量出杆的影子随时间变化的实际长度数据，

与问题三类似，结合模型用matlab曲线拟合求解即可，或者取多组数据用lingo软件求

解方程组。

模型假设

1.假设地球公转轨迹近似为圆。

2.忽略太阳光线进入大气层时的折射误差。

3.假设地面是水平的且直杆垂直地面。

4.忽略太阳直射点纬度一天内的变化。

5.假设所给数据准确可靠。

定义和符号说明

H：杆长

L：杆的影长

n：为从1月1号开始的天数

N：为当地经度



：为当地纬度

t：时间（小时）



：时角

ki

：区间均值

p

i

：解出现在某区间的频率

模型建立与求解

问题一：在水平地面上，垂直竖杆与在地面上的影长是一个直角三角形的两条直角

边，而太阳高度角α是杆长H的对角，运用正切公式可得：

L

H

tan

这样就把影长问题转化成了高度角问题。

一、构造高度角与各参数的关系图:

二、高度角与各参数关系模型构建

影响高度角的参数有太阳时角、赤纬角、杆的位置（经度、纬度）、时间。

时角，地球自转每小时15度，北京时间是以东经120度的为准的，故考虑时

差问题，可得:

)N120(12-t15

，

赤纬角



是地球赤道平面与太阳和地球中心的连线之间的夹角，也就是太阳直射点

的纬度。是因为地球公转形成的，因为地轴方向不变，所以赤纬角随地球在运行轨道上

的不同点具有不同的数值。赤纬角以年为周期，在南北纬23.45度范围内移动，考虑其

一年内运动规律，得赤纬角：



















365

n2842

sin45.23

，

结合当地纬度



，得太阳高度角与各参数间的关系模型（推算过程请见附件）：

coscoscossinsinsin

，

影子长度与各参数关系模型：

联合公式｛





coscoscossinsinsin

tan





L

H

可得：







coscoscossinsin

2^coscoscossinsin1







（

HL

三、根据模型求解：

将问题一中的地点位置，日期，杆长代入式中，用matlab软件画出了直杆影子

从9：00到15：00的长度变化曲线。（代码请见附件）。

影长曲线与时间的关系图如下图所示：

问题二：根据影子顶点数据来求地点，由于直杆高度无法确定，所以很难根据影子长度

来求太阳高度角问题，于是我们利用顶点坐标的轨迹来解决。

利用matlab将题中给出杆的影子顶点运动轨迹拟合了出来，如图所示，其中坐标原点为

杆的位置。可知影子发生了角度变化，以此来确定影子顶点与太阳的运动关系。

一、立体问题转化为平面问题

我们通过将太阳与地面直杆的立体方位关系转化为平面关系来建立模型，把地球的

自转看成太阳在直射点轨迹上方绕地球转，并且短时间内可近似看作太阳是作由东向西

的直线运动。并且存在两个关于直射点线（图中太阳运动轨迹）对称的位置，它们的影

子变化规律也是对称的。以下模型建立是视y轴指向北，x轴指向东（若y轴以南为正

方向，x轴以东为正方向，模型建立同理）。

如上图所示，此图可看作从太阳方向往地平面方向看去的俯视图，当太阳从A点运

动到B点时，杆影顶点从a移动到b，O点为杆的位置。由于太阳、杆顶、影子顶点永

远在一条直线上，据对顶角相等原理，可知太阳在时间段内移动过的经度差等于影子前

后转动的角度差。从而可以根据角度关系来确定杆位置。

二、模型建立

设北京时间t1~t2时刻太阳从A到B，可知太阳转过的经度为15（t2-t1），

由于北京时间12点时太阳所在的经度为东经120度，所以t1时刻太阳所在A点经

度为

120-15（t1-12）。

太阳运动轨迹为直射点轨迹，根据问题一中的公式:



















365

n2842

sin45.23

可知太阳运动路线的纬度。

影子从a点(x1,y1)到b点(x2,y2)，可知

1

1

tan

x

y



，

2

2

tan

x

y



。

假设经度差为N，纬度差为

，

由三角关系解方程组：｛

Nttx

y

Nx

y













)12(152

2

1

1





，

又因杆位置纬度



















365

n2842

sin45.23

+



,

经度N=15（t2-t1）+

N

,

解得杆位置与参数间的关系模型：

300)121(15

2121

21

)12(15





t

yxxy

yx

ttN



2112

21

)12(15

365

n2842

sin45.23

yxyx

yy

tt





















故若一直影子顶点坐标随时间变化的数据和日期，则可通过模型求出杆的位置。

三、模型求解

我们把题中给出的影子顶点坐标数据和日期代入模型进行求解，

每相邻两组坐标数据求一次解，利用matlab软件求得各组解的坐标分布图，如图所示：

由于各解的分布为离散型，所以我们对解进行权重处理。我们把各解的坐标看为离

散型随机变量，利用刻画离散型随机变量概率分布的方法来处理结果。把解出现在划分

的区间的个数作为该区间的频数，频率作为所占的权重p

i

（在某区间出现的频率越大，

则区间所占的权重越大），ki

为该区间的均值，最后结果为p

k

i

n

i

i





1

对经度处理：

区间频数频率

(110,111]10.05

(111,112]10.05

(112,113]30.15

(113,114]30.15

(114,115]30.15

(115,116]90.45

最后计算得经度N=p

k

i

i

i





6

1

=114.15

同理可求得纬度



=27.27，此处不赘述。

经查询，经纬度为114.15和27.27的地点位于江西吉安市安福县。

若影子顶点坐标y轴是以南为正方向，x轴以东为正方向。则该杆位置与上述求得

的位置关于直射点轨迹对称，计算得杆位置为经度N=114.15纬度



=—6.25。

位于印度尼西亚。

四、模型检验

将问题中的

问题三：问题三中已知的参数只有影子定点坐标和时间，所以无法利用问题二中的模型

来求杆的具体位置。我们联合问题一和问题二建立的模型（以下分别称为模型一、模型

二），然后用matlab对坐标数据进行曲线拟合，得出影长与时间的关系式，从而确定各

参数值。

一、模型结合

模型一：







coscoscossinsin

2^coscoscossinsin1







（

HL

模型二：

300)121(15

2121

21

)12(15





t

yxxy

yx

ttN

因为已经知道影子顶点坐标数据，故利用模型二先进行预处理，可以得到杆位置经

度N，由于日期未知，故模型二无法确定纬度。

把求得的经度代入模型一，可整理成函数关系式L),,(ntF，结合以下公式

)N120(12-t15



















365

n2842

sin45.23

输入影长与时间关系的数据利用matlab，据最小二乘法原理来进行曲线拟合，并且

关系式为给定的模型一，根据拟合结果确定参数



、n。则杆子位置的经纬度和日期就

可以求得。

二、模型求解

(一)附件二求解

（1)、利用模型二求经度

求题方法于问题二中已说明，此处不详解。将题给的附件二中坐标数据代入模

型二求解，整理得经度概率分布表：

区间频数

频率

(80,85]30.15

(85,90]60.3

(90,95]40.2

(95,100]50.25

(100,105]20.1

计算得经度N=p

k

i

i

i





5

1

=91.75度。

（2）利用模型一拟合曲线

首先，把坐标数据处理成影长随时间的关系，如下所示，为题中给出的附件2影长

与时间关系数据表格：

时间/h12.6813.6814.6815.6816.6817.6818.6819.6820.6821.6822.68

影长/m13.6814.6815.6816.6817.6818.6819.6820.6821.6822.6823.68

时间/h14.6815.6816.6817.6818.6819.6820.6821.6822.6823.68

影长/m15.6816.6817.6818.6819.6820.6821.6822.6823.6824.68

把求得的经度N=91.75影长与时间的关系式L=),,(ntF，用matlab软件拟合曲线，

进行参数估计，拟合结果如下图，

拟合结果得出了参数纬度



=—5.326度、天数n=357，

从而就确定了杆位置为经度N=91.75,纬度



=—5.326。日期为12月22日。

与问题二同理，赤纬角



=—21.8，存在关于太阳直射点轨迹的对称点



=—38.3，N=91.75。

（二）附件三求解

附件三数据求解过程与附件二同理，此处不赘述。（详情见附录）

求解结果为杆位置经度N=115.8，纬度



=5.6日期为12月20日。

问题四

首先，根据视频中可得在40分钟里直杆在太阳下影长变化规律。每隔两分钟截一张

图可得20张图片。利用photoshop软件中的度量工具可得图片中直杆的高度与影子的长

度可得一组数据。由于直杆的高度为2米，可得实际影长的长度数据如下：

（1）

时间/h8:558:578:599:019:039:059:079:099:119:13

影长/m2.402.382.352.322.292.272.242.202.182.15

时间/h

9:159:179:199:219:239:259:279:299:319:33

影长/m

2.122.102.072.052.011.991.971.941.911.89

日期已知，为7月13日，将天数n代入，用matlab软件对影长与时间的关系式进行

曲线拟合，确定关系式







coscoscossinsin

2^coscoscossinsin1







（

HL

得曲线图：

得参数经度N=-18.23度纬度



=21.86度

（2）

日期不确定，则有三个不确定的参数，三个未知量可用三个方程组求得解。用lingo

软件输入三组数据，即可确定三个参数。

误差分析：

1、不同时刻大气对光线的折射率不同；

2、地球表面不是完全水平面；

3、量取视频中影子长度时由于影子是斜的和视角问题，所量值非其实际值；

4、数据计算取小数保留；

5、曲线拟合存在误差；

模型评价

模型优点：

1、简洁性：模型原理简易，不需要太复杂的数学知识，易为大众所了解掌握。

2、确定性：模型基本都是由公式组成，有较强的严谨性，计算出的都是确定的解。

3、可靠性，模型建立在坚实可靠的数学基础上。

模型局限性：

1.计算公式过于繁琐，必要时得依赖计算机求解。

2.模型只是建立的假设条件只在短时间内成立。

3.应用曲线拟合只是大体上估计，存在较大偏差，只能拟合出大概解。

模型改进与推广：

模型建立条件是没有考虑大气层的折射问题，这可能在一定程度上影响模型准确度。

所以应该在这方面改进，使其更接近于现实。

参考文献：

【1】方荣生，《太阳能应用技术》，北京：中国农业机械出版社，1985。

【2】王松桂，《概率论与数理统计》，北京：科学出版社，2011。

【3】杨长青.正午太阳高度角模型”的设计、制作与运用.内蒙古师范大学地理科学

学院2011-08-01。

【4】王国安.太阳高度角和日出日落时刻太阳方位角一年变化范围的计算.河南省气

候中心2007-09-15。

【5】刘卫国，《MATLAB程序设计教程》，北京：中国水利水电出版社，2005。

附录

问题一相关代码：

b=[9:0.05:15];

c=15*(b-12)-(120-116.39139);

a=23.45*sin(2*3.1415926*(284+295)/365);

i=39.907222;

t=sind(i)*sind(a)+cosd(i)*cosd(a)*cosd(c);

x=asind(t);

d=3./tand(x);

plot(b,d)

问题二相关代码：

x=[1.03651.06991.10381.13831.17321.20871.24481.28151.31891.35681.3955

1.43491.47511.5161.55771.60031.64381.68821.73371.78011.8277

];

y=[0.49730.50290.50850.51420.51980.52550.53110.53680.54260.54830.5541

0.55980.56570.57150.57740.58330.58920.59520.60130.60740.6135

];

t=[14.714.7514.814.8514.914.951515.0515.115.1515.215.2515.315.35

15.415.4515.515.5515.615.6515.7

];

fori=1:20

e(i)=120-(t(i)-12)*15+15*(t(2)-t(1))*x(i)*y(i+1)/(y(i)*x(i+1)-x(i)*y(i+1));

n(i)=10.5109+15*(t(2)-t(1))*y(i)*y(i+1)/(y(i)*x(i+1)-x(i)*y(i+1));

end

scatter(e,n,'r*')

figure

scatter(n,e,'b*')

杆位置离散图：

问题三中附件二杆位置离散图：

相关代码：

x=[1.23521.20811.18131.15461.12811.10181.07561.04961.02370.998

0.97240.9470.92170.89650.87140.84640.82150.79670.77190.7473

0.7227];

y=[0.1730.1890.20480.22030.23560.25050.26530.27980.2940.3080.3218

0.33540.34880.36190.37480.38760.40010.41240.42460.43660.4484];

t=[12.6833333312.7333333312.7833333312.8333333312.88333333

12.9333333312.9833333313...13333333

13.1833333313.2333333313.2833333313.3333333313.38333333

13.4333333313.4833333313.5333333313.5833333313.63333333

13.68333333];

fori=1:20

区间频数频率

(110,111]10.05

(111,112]10.05

(112,113]30.15

(113,114]30.15

(114,115]30.15

(115,116]90.45

区间频数频率

(80,85]30.15

(85,90]60.3

(90,95]40.2

(95,100]50.25

(100,105]20.1

e(i)=120-(t(i)-12)*15-15*(t(2)-t(1))-15*(t(2)-t(1))*x(i+1)*y(i)/(y(i+1)*x(i)-x(i+1)*y(i));

n(i)=15*(t(2)-t(1))*y(i)*y(i+1)/(y(i+1)*x(i)-x(i+1)*y(i));

end

stem(e,n,'r*')

figure

stem(n,e,'b*')

y=[1.2472562051.222794591.1989214861.1754289641.152439573

1.129917471.107835481...044446265

1...9854909080.9667904940.948584735

0.9309278810.913751750.8971090510.8809737620.865492259

0.850504468

];

x=[12.6833333312.7333333312.7833333312.8333333312.88333333

12.9333333312.9833333313...13333333

13.1833333313.2333333313.2833333313.3333333313.38333333

13.4333333313.4833333313.5333333313.5833333313.63333333

13.68333333];

y=h*(1-(sind(a)*sind(23.45*sin(2*3.1415926*(284+n)/365))+cosd(a)*cosd(23.45*sin(2*3.14

15926*(284+n)/365))*cosd(15*(x-12)-120+91.75))^2)^(1/2)/(sind(a)*sind(23.45*sin(2*3.141

5926*(284+n)/365))+cosd(a)*cosd(23.45*sin(2*3.1415926*(284+n)/365))*cosd(15*(x-12)-1

20+91.75))

附件三中位置离散图：：

附件三相关代码：

x=[1.16371.22121.27911.33731.3961.45521.51481.5751.63571.697

1.75891.82151.88481.94882.01362.07922.14572.21312.28152.3508

2.4213];

y=[3.3363.32993.32423.31883.31373.30913.30483.30073.29713.2937

3.29073.28813.28593.2843.28243.28133.28053.28013.28013.2804

3.2812];

t=[13.1513.213.2513.313.3513.413.4513.513.5513.613.65

13.713.7513.813.8513.913.951414.0514.114.15];

fori=1:20

e(i)=120-(t(i)-12)*15+15*(t(2)-t(1))*x(i)*y(i+1)/(y(i)*x(i+1)-x(i)*y(i+1));

n(i)=15*(t(2)-t(1))*y(i)*y(i+1)/(y(i)*x(i+1)-x(i)*y(i+1));

end

stem(e,n,'k*')

figure

stem(n,e,'k*')

y=[3.5331421843.5467680293.5617976433.5781007153.595750783

3.614934283.6354259833.6572182723.6805411153.705167836

3.7312780253.7589179113.7880878883.8187010153.850809619

3.884585223.9199118283.9568759923.995534794..077863059];

x=[13.1513.213.2513.313.3513.413.4513.513.5513.6

13.6513.713.7513.813.8513.913.951414.0514.1

14.15];

y=h*(1-(sind(a)*sind(23.45*sin(2*3.1415926*(284+n)/365))+cosd(a)*cosd(23.45*sin(2*3.14

15926*(284+n)/365))*cosd(15*(x-12)-120+115.8))^2)^(1/2)/(sind(a)*sind(23.45*sin(2*3.141

5926*(284+n)/365))+cosd(a)*cosd(23.45*sin(2*3.1415926*(284+n)/365))*cosd(15*(x-12)-1

20+115.8))

问题四相关代码：

x=[8.9166666678.958.9833333339..059..1166666679.15

9.1833333339.2166666679.259.2833333339.3166666679.359.3833333339.416666667

9.459.4833333339.5166666679.55];

y=[2.4028776982.3836930462.3533173462.3245403682.2925659472.266986411

2.2414068752.203037572.1790567552.156674662.1262989612.10871303

2....9920063951.968025581.93764988

1.9136690651.889688249];

cftool

y=2*(1-(0.372*sind(a)+cosd(a)*0.928*cosd(15*(x-12)-(120-n)))^2)^(1/2)/(0.372*sind(a)+cosd(a)*0.928*c

osd(15*(x-12)-(120-n)))

附件三求解过程

（1)、利用模型二求经度

将题给的附件三中坐标数据代入模型二求解，整理得经度概率分布表：

计算得经度N=p

k

i

i

i





5

1

=115.8度。

（2）利用模型一拟合曲线

首先，把坐标数据处理成影长随时间的关系，

如下所示，为题中给出的附件三影长与时间关系数

据表格：

时间/h12.6813.6814.6815.6816.6817.6818.6819.6820.6821.6822.68

影长/m13.6814.6815.6816.6817.6818.6819.6820.6821.6822.6823.68

时间/h14.6815.6816.6817.6818.6819.6820.6821.6822.6823.68

影长/m15.6816.6817.6818.6819.6820.6821.6822.6823.6824.68

把求得的经度N=115.8影长与时间的关系式L=),,(ntF，用matlab软件拟合曲线，

进行参数估计，拟合结果如下图，

区间频数频率

(113.5,114.5]30.15

(114.5,115.5]50.25

(115.5,116.5]50.25

(116.5,117.5]70.35

拟合结果得出了参数纬度



=—5.326度、天数n=357，

从而就确定了杆位置为经度N=115.8,纬度



=—5.326。日期为12月22日。

与问题二同理，赤纬角



=—21.8，存在关于太阳直射点轨迹的对称点



=—38.3，N=115.8。