数学建模论文太阳影子定位模型

太阳影子定位的数学模型

摘要

太阳影子的变化与人类生活紧密相连，本文分别针对影子长度变化问题和太阳定位问题建立模型并求解，最后分析了模型的优缺点并给出了改进方案。

针对问题一，首先通过查阅资料，确定影响影长变化的直接参数是杆长和太阳高度角，并得到杆长、影长、太阳高度角的几何关系。然后基于太阳高度角、太阳赤纬、当地的北京时间、观测地纬度、观测地经度之间的关系式，建立关于杆长、太阳赤纬等五个参数的影长变化模型。最后利用matlab编程，代入题中所给参数数据，实现模型的计算。

针对问题二，仅考虑已知测量日期的太阳影子定位。首先假定影子顶点坐标系x轴的正方向指向正南方，根据太阳方位角的定义计算出实际方位角。然后，基于经度、纬度、太阳方位角间的关系式，根据实际方位角和太阳方位角之间的匹配程度建立模型，并将其划归为线性规划模型。最后，穷举所有匹配的经纬度求出最优解，并结合人工选择给出可能的地点。

针对问题三，考虑未知测量日期的太阳影子定位，首先对问题二所建立的模型进行匹配修改，在问题二中的模型添入了太阳赤纬作为决策变量，并根据太阳赤纬取值范围添加约束条件，建立了新的规划模型。然后，为了降低枚举法的时间复杂度，通过最小二乘法拟合影长求得时差，基于时差和经度的关系式确定出经度，减少了枚举变量。最后枚举了所有纬度、太阳赤纬，确定了测量地点和日期。

针对问题四，首先类比问题二中模型建立的方法，建立根据杆长、影长数据确定拍摄地点的模型。基于经度、纬度、太阳高度角间的关系式，改变问题二中模型的目标函数、约束条件建立新的规划模型。求解模型时，首先利用FFDshow将视频分解为若干图像（每两分钟一张），再将图像进行数字化处理，转换为二值图像后得到各图像中影长数据，最后代入已提取的影长数据，利用枚举法在经纬度的取值范围内搜索出可能的地点。

关键词：太阳高度角、太阳方位角、规划模型、搜索算法、枚举法优化、图像数字化处理

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

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

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

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

4．附件4为一根直杆在太阳下的影子变化的视频，并且已通过某种方式估计出直杆的高度为2米。需建立确定视频拍摄地点的数学模型，并应用该模型给出若干个可能的拍摄地点。考虑在拍摄日期未知的情况下，能否根据视频确定出拍摄地点与日期。

2.1问题一的分析

由于影子长度变化受地球自转、地球公转、黄赤交角等多个参数的影响，建立影子长度变化模型时需要抓住关键参数，而影响影长变化最直接的参数是杆长和太阳高度角（与太阳赤纬、观测地经纬度、时角相关），因此选择杆长、太阳高度角作为关键参数，结合各参数间的关系式，进行合理的假设和建模。采取控制参数法，每次改变模型中的某一个参数，其余参数固定不变，可分析出模型中各个参数对影长变化的影响。将题中的杆长、经纬度等已知参数带入模型，利用matlab即可绘制出在9:00-15:00影长变化的图像。

2.2问题二的分析

建立数学模型来根据太阳影子顶点坐标数据确定直杆地点，即确定所在地点的经纬度。因此，可选择经度、纬度作为决策变量，利用太阳方位角与实际方位角（可由太阳影子顶点坐标确定）的关系式确定目标函数，利用经度、纬度、太阳方位角间的关系式以及经纬度的取值范围确定约束条件，建立合理的规划模型。在求解模型时，只需合理地规定经度、纬度的变动步长，利用枚举法在经纬度的取值范围内搜索出最优解。

2.3问题三的分析

针对问题三，基于问题二中所建立的太阳影子定位模型，分析可发现日期仅与太阳赤纬相关，因此，在问题二的模型中添加太阳赤纬作为决策变量，并根据太阳赤纬的取值范围添加约束条件，建立出新模型。求解模型时，为降低枚举的时间复杂度，考虑减少枚举变量。因此，可通过最小二乘法对影长（可由影子顶点坐标算得）数据进行拟合求得时差，基于时差、经度的关系式将经度确定，最后只须枚举纬度和太阳赤纬，即可搜索出最优解。

2.4问题四的分析

基于问题二所建立的模型，通过太阳高度角和实际太阳高度角（由影长、杆长确定）的关系式，改变模型中的目标函数、约束条件，即可建立新的规划模型来确定视频拍摄地点。求解模型的方法与问题二中基本相同，关键在于如何从视频中提取影长数据，首先利用FFDshow将视频分解为若干图像，再将图像进行数字化处理，转换为二值图像后得到各图像中影长数据。

三、基本假设与符号说明

3.1基本假设

1.将太阳光近似看作平行光；

2.固定直杆附近无高大建筑物遮挡阳光；

3.固定直杆以及影子的直径可忽略；

4.不考虑大气折射对太阳光的影响；

5.假定太阳照射稳定，不受云层等自然因素的影响；

6.假定问题二中的坐标系x轴的正方向指向正南方。

3.2符号说明符号

h符号说明太阳高度角

太阳赤纬

观测地纬度（北纬为正，南纬为负）

观测地经度（东经为正，西经为负）

时角

杆长

影长

真太阳时（地方时）

当地的北京时间（简称北京时间）

时差（真太阳时与北京时间的时差）

经度差（观测地与北京的经度差）单位tLlSTmmhhhT

i第i组纬度

第i组经度

实际的太阳方位角（由影子顶点坐标x,y确定）

第i组的太阳方位角

第i组太阳赤纬

实际的太阳高度角（由影长l、杆长L确定）

iAsAiihs

四、模型建立与求解

4.1问题一模型的建立与求解

4.1.1模型的建立

对于地球上的某个地点，太阳高度角是指太阳光的入射方向和地平面之间的夹角，一根直立的杆其影子的长度在一天中随太阳高度角的变化而不断变化，且杆长L、影长l与太阳高度角h满足图1的几何关系。

L光线

图1杆长、影长与太阳高度角间的几何关系示意图h(1)由上图的几何关系，根据正切函数可得到杆长L、影长l与太阳高度角h之间的关系式：

lLcoth

(2)又由于太阳高度角随着时角t、太阳的赤纬（与太阳直射点纬度相等）和观测地纬度（太阳赤纬与观测地纬度都是北纬为正，南纬为负）的变化而变化，参考百度百科中有关太阳高度角的文献，可得太阳高度角的计算公式：[1]

sinhsinsincoscoscost

根据时角的相关资料可知，规定真太阳时ST12时的时角为零，前后每

隔一小时增加3602415，由此可求得时角t的表达式：[2]

t15ST12

杆长L、影长l与太阳高度角h的关系式中需求得太阳高度角的余切值coth，因此首先可根据反三角函数的求出太阳高度角h，进而根据三角函数间的关系求出其余切值coth：

harcsin[sinsincoscoscos(15ST12)]

cothcot{arcsin[sinsincoscoscos(15ST12)]}

(3)由于需要建立影长l关于太阳赤纬、观测地纬度、观测地经度、北京时间T、杆长L的模型，因此需要将上式中的真太阳时ST以北京时间T替代，

参阅资料了解到真太阳时ST与北京时间T的关系式如下：

120STTTT1515

表示观测地与北京的经度差。

（4）通过（1）——（3）中得到的关系式，可建立影长l随着太阳赤纬、观测地纬度、观测地经度、北京时间T、杆长L变化的模型：[3]表示观测地与北京的时差，

lLcot{arcsin[sinsincoscoscos(15T

该模型中影长l为因变量，其余参数均为自变量。12012)]}15

4.1.2模型的分析

根据控制参数的思想，分析所建模型中影长l关于各参数的变化规律。利用matlab编写程序（见附录一），依次改变太阳赤纬、观测地纬度、北京时间T、杆长L等参数，分析出影长l与各个参数的变化规律如下：

（1）固定观测地纬度39.9072、北京时间T10、经度120、杆长L3，仅有太阳赤纬变化时，影长l随着太阳赤纬的增大而减小，见图2。

影长与太阳赤纬的关系

l（-25-20-15-10-50510

q（太阳赤纬单位：°）152025

图2影长与太阳赤纬的关系图

（2）固定太阳赤纬10.4、北京时间T10、经度120、杆长L3，仅有观测地纬度变化时，影长l随着观测地纬度的变化如图3。

影长与观测地纬度的关系

5.5

4.5

3.5

3l（2.5

1.5-50-40-30-20-1001020

p（观测地纬度单位：°）304050

图3影长与观测地的关系图

（3）固定太阳赤纬10.4、观测地纬度39.9072、经度120、杆长L3，仅有北京时间T变化时，影长l随着时刻T的变化如图4。

影长与时间的关系

l（910111213

T（北京时间单位：h）1415

图4影长与北京时间的关系图

（4）固定太阳赤纬10.4、观测地纬度39.9072、经度120、北京时间T10，仅有杆长L变化时，影长l随着杆长L的增长而增长,见图5。

影长与杆长的关系

l（6

002468

L（杆长单位：m）1012

图5影长与杆长的关系图

4.1.3模型的求解

将题中已知的参数代入所建模型，画出太阳影子长度在北京时间9:00-15:00之间的变化曲线。

（1）直杆长L3；

（2）地理纬度39.9072；

（3）查阅赤纬DE表（见附录二）可知2015年10月22日的太阳赤纬：

10.4

（4）已知地理经度116.3914，可计算经度差：

120120116.39143.6086

将（1）——（4）中数据带入模型，得到以北京时间T为自变量，影长l为因变量的关系式：

l3cot{arcsin[sin39.9072sin(10.4)

cos39.9072cos(10.4)cos(15T12.2406)]}

利用matlab编写程序（见附录一），可绘制影子长度的变化曲线如图6。

影长与时间的关系7.5

6.5

5.5

l（4.5

3.5910111213

T（时间单位：h）1415

图69:00-15:00影长变化图

4.2问题二模型的建立与求解

4.2.1模型的准备

（1）太阳方位角：太阳方位角即太阳所在的方位，指太阳光线在地平面上的投影与当地经线的夹角，可近似地看作是竖立在地面上的直线在阳光下的阴影与正南方的夹角[4]。假设坐标系x轴的正方向指向正南方，可根据某固定直杆在水平地面上的太阳影子顶点坐标数据求出实际方位角A：

Asarctan

[5]sxy根据天体理论，可以导出太阳方位角A与太阳高度角h、太阳赤纬、时

角t的关系如下：

Aiarcsin(cossintii1n)coshi

（2）范数：范数是在广义长度意义下，对函数、向量或矩阵的一种度

量定义。使用范数可以测量两个函数、向量或矩阵之间的距离。对于向量[6]

Xx1,

,xn经常使用的一种范数是：

X

4.2.2模型的建立

（1）目标分析：太阳方位角Ai与实际方位角A间的距离越小，则纬度i、经度i与直杆所处的地点的经、纬度越接近。因此，确定目标为太阳方位角Ai与实际方位角As间距离最小，得到目标函数：s

minAiAs

（2）约束条件：

根据经、纬度取值范围得到不等式约束：

90i90 (i1n)

180i180 (i1n)

基于模型的准备，可知太阳方位角

关系式：

Aiarcsin(cossinti) coshiAit与太阳高度角hi、太阳赤纬、时角i的基于问题一，可知时角ti、经度i的关系式：

ti15T

可知太阳高度角hi、纬度120i1215i、经度i的关系式：

120i12)]15hiarcsin[sinisincosicoscos(15T

上述三个关系式，即为模型中的等式约束。（3）综上可得规划模型：

minAiAs

cossintis.. tAarcsin() (i1n)icoshi120i ti15T12 (i1n)15 harcsin[sinsincoscoscos(15T120i12)] (i1n)iii15 90i90 (i1n)

 180i180 (i1n)

4.2.3模型的求解

（1）计算太阳赤纬：

附件1中给出了测量日期，查阅赤纬DE表（见附录二）可知测量日期2015年4月18日的太阳赤纬：10.2

（2）计算实际方位角As：

基于模型的准备可知顶点坐标xj,yj与实际方位角的Asj的关系式，将附件1中21组顶点坐标xj,yj代入关系式，可求每组坐标对应的实际方位角Asj：

Aarctan

则实际方位角A可表示为：

xjyj

(j121)

AsA1sA21

（3）计算太阳方位角Ai：

T j121将附件1中21组北京时间j代入太阳方位角与太阳赤纬、北京

T A

时间、经度、纬度间的关系式，可求出每组北京时间j所对应的太阳方位角i,j：

120i

cos10.2sin15Tj12

15Ai,jarcsin (j121)

120icosarcsinsinisin10.2cosicos10.2cos15Tj12

15

则太阳方位角可表示为：

AiAi,1Ai,21

（4）计算目标函数：

基于模型的准备可知，使用范数可以测量示为：

AiA

As,Ai

间的距离，则目标函数可表



A

（5）利用Matlab编写程序（见附录一），合理设置经、纬度的步长，得到最优解。

为得出多个可能地点，设置经、纬度的不同步长，运行程序可得到以下结果（表1）：

步长(纬度，经度)

（0.5,1）

经纬度(纬度，经度)

（-34,139）

地点澳大利亚

（0.2,0.2）（4,0.1）（3,0.1）（4，0.3）

表1

（-89.2140.2）（6，-67.3）（-3,37.6）（46,19.8）

不同搜索步长下确定的地点

南极洲委内瑞拉肯尼亚塞尔维亚

4.3问题三模型的建立与求解4.3.1模型的建立

添加太阳赤纬i作为决策变量，查阅资料

23.5的范围内移动。

[7]

可得，太阳赤纬i在23.5到

基于问题二中所建立的模型，得到问题三的模型：

minAiAs

cosisinti

s.. tAarcsin(i1n)icoshi120i

12 (i1n) ti15T

15

120i

harcsin[sinsincoscoscos(15T12)] (i1n)iiiii15

 90i90 (i1n) 180180 (i1n)

i

 23.5i23.5 (i1n)

i为第i组地理纬度，i为第i组地理经度，i为第i组太阳赤纬，A为第i组

太阳方位角，T为北京时间，A为实际方位角。

4.3.3模型的求解

（1）计算当地经度：为降低枚举的时间复杂度，考虑确定经度，以此减少枚举变量。确定经度的步骤（以附件2为例）如下。

a.计算影长l:

根据x,y坐标和影长l的几何关系，利用勾股定理得到坐标与影长的关系式：

ljxj2yj2

代入附件2中的x,y坐标，利用excel编写函数求解出影长ll1l21b.计算当地正午时间Tm：

基于问题一，可知影长l关于北京时间T的函数为二次函数。利用matlab软件编写二次拟合程序（见附件一），代入附件2中的21组Tj,lj，可得函数表达式与图像（图7）：

l0.09814T22.985T23.32

拟合附件2曲线

l （91011

12131415

T （时间单位：h）

161718

图7影长与北京时间的关系图

根据拟合出的函数可求得当地的正午时间

Tm15.2079c.计算经度：

基于问题一，可知经度与时差的关系：

120

12Tm

代入（2）中求得的正午时间Tm可求得经度：71.8860

同理可得，附件3中108.9315（附件3中影长l关于北京时间T的函数表达式：l0.2964T27.551T51.56）

（2）计算实际方位角As：将附件2、附件3中21组顶点坐标

实际方位角j：

xj,yj

代入关系式，可求每组坐标对应的

Aarctan

则实际方位角A可表示为：

xjyj

(j121)

AsA1sA21

（3）计算太阳方位角Ai：

将附件2、附件3中21组北京时间

Tj j121

和在（1）中求得的经度

依次代入太阳方位角与太阳赤纬、北京时间、经度、纬度间的关系式，可求出每

T A

组北京时间j所对应的太阳方位角i,j：

120

cosisin15Tj12

15Ai,jarcsin (j121)

120

cosarcsinsinisinicosicosicos15Tj12

15

则太阳方位角Ai可表示为：

AiAi,1Ai,21

（4）计算目标函数：

AiA



A

（5）求解附件2的测量地点和日期：

利用Matlab编写程序（见附录一），代入已求出附件2中的经度、影子顶点坐标数据，设置纬度i的步长为0.5、太阳赤纬i的步长为0.1，运行程序可得到结果：

24.96.2

查阅地图、赤纬表可得，测量地点为广西壮族自治区河池市罗城仫佬族自治县204省道测量时间为平年的4月7日或闰年的4月6日。

（5）求解附件3的测量地点和日期：

利用Matlab编写程序（见附录一），代入附件3中的的经度、影子顶点坐标数据，设置纬度i的步长为0.2、太阳赤纬i的步长为0.1，运行程序可得到结果：

11

闰年的8月22日。

4.4问题四模型的建立与求解4.4.1模型的建立

12.1

查阅地图、赤纬表可得，测量地点为印度，测量时间为平年的8月23日或

（1）目标分析：太阳高度角hi与实际高度角hs间的距离越小，则纬度i、经度i与直杆所处的地点的经、纬度越接近。因此，确定目标为太阳高度角hi与实际高度角hs间距离最小，得到目标函数：

minhihs

（2）约束条件分析：

根据经、纬度取值范围得到不等式约束：

90i90 (i1n)180i180 (i1n)

基于模型的准备，可知太阳高度角hi与太阳赤纬、地理纬度i、时角ti的关系式：

hiarcsin[sinisincosicoscosti]

t基于问题一，可知时角i、经度i的关系式：

ti15T

120i

1215

综上，可知太阳高度角hi、纬度i、经度i的关系式：

hiarcsin[sinisincosicoscos(15T

（3）综上可得规划模型：

120i

12)]15

minhihs120is.. thiarcsin[sinisincosicoscos(15T12)] (i1n)15

 90i90 (i1n) 180i180 (i1n)

i为第i组地理纬度，i为第i组地理经度，为太阳赤纬，h为第i组太阳

方位角，T为北京时间，hs为实际高度角。

4.4.2模型的求解

（1）计算影长li：

a.利用FFDshow软件，将视频每2分钟分解共得到21张图像。

b.利用matlab软件编写程序（见附录一），对得到的图像进行数字化处理，转换为二值图像后，分别得到图像中影子起点的坐标x1,y1、终点坐标x2,y2。利用坐标间的距离公式：

l

值，可求得21组影长数据。

（2）计算太阳赤纬：

视频中给出了测量日期，查阅赤纬DE表（见附录二）可知测量日期2015年7月13日的太阳赤纬：22.1

（3）计算实际高度角hs：

基于问题一，可知影长lj、杆长L与实际方位角的hsj间的关系式：

ljLcothsj

将杆长L2以及从视频中提取的21组影长lj代入上述关系式，可求每组影长lj对应的实际高度角hsj：

hsjarccot

则实际高度角hsj可表示为：

ljL

(j121)

hsh1sh21

（4）计算太阳高度角hi：

将视频中的21组北京时间Tj等式约束，可求出每组北京时间Tj所对应的太阳高度角hi,j：

hi,jarcsin[sinisincosicoscos(15Tj则太阳高度角hi可表示为：

120i

12)] (j121)15

titi,1ti,21

（5）计算目标函数：

基于模型的准备可知，使用范数可以测量ts,ti间的距离，则目标函数可表示

为：

tit



t



（6）利用Matlab编写程序（见附录一），设置经度的步长为0.5度、纬度的步长为0.5度，得到最优解为：纬度50.9，经度96.7，地点为俄罗斯。

五、模型的评价与改进

5.1模型的评价

计算简便，结果合理

优点

具有适用性，可根据影子定位模型推导出确定日期的模型能与实际紧密联系,结合实际情况对所提出的问题进行求解利用枚举法求解，比较直观，易于理解未考虑坐标轴实际朝向，模型不具有推广性

缺点

利用枚举法求解，运算量比较大，解题效率不高进行图像的数字化处理时，存在一定误差

5.2模型的改进

现有模型的算法对于相关参数比较敏感，为了提高算法的适用性，可以尝试引入模拟退火算法、遗传算法、蚁群算法等智能算法。

参考文献

[1]百度百科，太阳高度角_百度百科，https://www.baidu.com/link?url=Ow1IVmyA1ceNTRid=afbe1c08000714c[1**********]f3a3bd，2015年9月11日

[2]百度百科，时角_百度百科，，2015年9月11日

[3]百度知道，地方时计算公式，http://zhidao.baidu.com/link?url=k-LZkgh_VV11dMIgM6，2015年9月11日

[4]百度百科，太阳方位角_百度百科，mw7UsOsr6s78j-_eB5pPHUzbbyNzj-XM5cuBbJ7DHMitA5XtTylNpEgbjbIfEDCcuHl01A0qOha，2015年9月12日

[5]贺晓雷等，太阳方位角的公式求解及其应用，太阳能学报，第29卷第1期：p60，2018年

[6]张韵华等，数值计算方法与算法，北京：科学出版社，2006

[7]百度百科，太阳赤纬_百度百科，http://baike.baidu.com/link?url=NqXYAQsPP852Jj，2015年9月11日

附录

附录一：模拟所用的程序

（1）影子长度关于地方时（时间）变化的变化规律

%t为地方时(时角)，h为太阳赤纬，T时间，l影长，q太阳赤纬，p观测地纬度,L杆长%固定杆长、太阳赤纬、观测地纬度，考虑影子长度关于地方时（时间）变化的变化规律q=-10.4;L=3;

%q太阳赤纬%p观测地纬度%m杆长

p=39.9072;

T=[9:0.25:15];%步长15分钟=0.25小时t=abs(T-12).*15;

h=asind(sind(p)*sind(q)+cosd(p)*cosd(q).*cosd(t));l=L./tand(h);plot(T,l)

%求影长

%求太阳高度角

%画出影长y轴与时间x轴的关系

title('影长与时间的关系')xlabel('T（北京时间单位：h）')ylabel('l（影长单位：m）')（2）影子长度关于杆长变化的变化规律

%固定地方时(上午10点)、太阳赤纬、观测地纬度，考虑影子长度关于杆长变化的变化规律t=30;

%地方时

%q太阳赤纬%p观测地纬度

%m杆长

%求太阳赤纬

q=-10.4;

p=39.9072;

L=[0:0.5:12];

h=asind(sind(p)*sind(q)+cosd(p)*cosd(q).*cosd(t));l=L./tand(h);plot(L,l)

%求影长

%画出影长y轴与时间x轴的关系

title('影长与杆长的关系')xlabel('L（杆长单位：m）')ylabel('l（影长单位：m）')

（3）影子长度关于太阳赤纬变化的变化规律

%固定地方时(上午10点)、杆长、观测地纬度，考虑影子长度关于太阳赤纬变化的变化规律t=30;L=3;

%地方时

%p观测地纬度

%q太阳赤纬

%求太阳赤纬

%m杆长

p=39.9072;

q=[-23.5:2:25.5];

h=asind(sind(p)*sind(q)+cosd(p)*cosd(q).*cosd(t));l=L./tand(h);

%求影长

plot(q,l)%画出影长y轴与时间x轴的关系

title('影长与太阳赤纬的关系')xlabel('q（太阳赤纬单位：°）')ylabel('l（影长单位：m）')

（4）影子长度关于观测地纬度变化的变化规律

%固定地方时(上午10点)、杆长、太阳赤纬，考虑影子长度关于观测地纬度变化的变化规律t=30;L=3;q=-10.4;

%地方时%m杆长%q太阳赤纬

%p观测地纬度

%求太阳高度角

p=[-45:1:45];

h=asind(sind(p)*sind(q)+cosd(p)*cosd(q).*cosd(t));l=L./tand(h);plot(p,l)

%求影长

%画出影长y轴与时间x轴的关系

title('影长与观测地纬度的关系')xlabel('p（观测地纬度单位：°）')ylabel('l（影长单位：m）')（5）太阳影子定位模型的求解代码

%方位角1A1，方位角2A2，误差值r，影子端点坐标（X,Y）,i=1~21,j=1~21AW=[];BE=[];FW=[];FE=[];A1_t=[];A2_t=[];STOR_A12=[];STOR_W=[];STOR_E=[];suma1=0;suma2=0;X_A1=0;X_A2=0;X_A12=0;

min_x_a12=1000;CX=[];

%T存放当地北京时刻，单位小时

T=[14.7014.7514.8014.8514.9014.9515.0015.0515.1015.1515.2015.2515.3015.3515.4015.4515.5015.5515.6015.6515.70];

%存放影子端点坐标

X=[1.03651.06991.10381.13831.17321.20871.24481.28151.31891.35681.39551.43491.47511.5161.55771.60031.64381.68821.73371.78011.8277];

Y=[0.49730.50290.50850.51420.51980.52550.53110.53680.54260.54830.55410.55980.56570.57150.57740.58330.58920.59520.60130.60740.6135];q=10.2;

forW=-90:0.2:90

%纬度，北纬为正，南纬为负

%经度，东经为正，西经为负

%求当地时角（西经）

%求太阳赤

forE=-180:0.2:180

fori=1:1:21

t=abs((T(i)-(120-E)/15-12))*15;

纬h

h=asind(sind(W)*sind(q)+cosd(W)*cosd(q)*cosd(t));

A2=asind(cosd(q)*sind(t)/cosd(h));%求从太阳高度角和赤纬和时角

得到的方位角A2

A1=atand(X(i)/Y(i));a1=A1;

A1_t=[AW,a1];AW=A1_t;a2=A2;

A2_t=[BE,a2];BE=A2_t;end

forj=1:1:21

suma1=(A1_t(j))^2+suma1;suma2=(A2_t(j))^2+suma2;end

X_A1=sqrt(suma1);X_A2=sqrt(suma2);X_A12=abs(X_A1-X_A2);if(X_A12

min_x_a12=X_A12;FAL_W=W;FAL_E=E;endAW=[];BE=[];

%求从影子端点坐标得到的方位角

A1_t=[];A2_t=[];suma1=0;suma2=0;X_A1=0;X_A2=0;X_A12=0;endendFAL_WFAL_E

（5）问题三模型关于附件2的求解代码

%方位角1A1，方位角2A2，误差值r，影子端点坐标（X,Y）,i=1~21,j=1~21AW=[];BE=[];FW=[];FE=[];QT=[];A1_t=[];A2_t=[];STOR_A12=[];STOR_W=[];STOR_E=[];suma1=0;suma2=0;X_A1=0;X_A2=0;X_A12=0;

min_a12_step=1000;CX=[];

%T存放当地北京时刻，单位小时

T=[14.7014.7514.8014.8514.9014.9515.0015.0515.1015.1515.2015.2515.3015.3515.4015.4515.5015.5515.6015.6515.70];

%存放影子端点坐标

X=[1.03651.06991.10381.13831.17321.20871.24481.28151.31891.35681.39551.43491.47511.5161.55771.60031.64381.68821.73371.78011.8277];

Y=[0.49730.50290.50850.51420.51980.52550.53110.53680.54260.5483

0.55410.55980.56570.57150.57740.58330.58920.59520.60130.60740.6135];

%拟合结果：附件二中最小的影长对应时间为tmin=15.2079；附件三中最小的影长对应的时间为tmin=12.7379；E=120-(tmin-12)*15%E2=71.8860E3=108.9315E=71.8860;forW=-90:0.5:90

%纬度，北纬为正，南纬为负

%赤纬q

%求当地时角（西经）

%求太

forq=-23.4:0.1:23.4

fori=1:1:21

t=abs((T(i)-(120-E)/15-12))*15;

阳赤纬h

h=asind(sind(W)*sind(q)+cosd(W)*cosd(q)*cosd(t));

A2=asind(cosd(q)*sind(t)/cosd(h));%求从太阳高度角和赤纬和

时角得到的方位角A2

A1=atand(X(i)/Y(i));a1=A1;

A1_t=[AW,a1];AW=A1_t;a2=A2;

A2_t=[BE,a2];BE=A2_t;end

forj=1:1:21

suma1=(A1_t(j))^2+suma1;suma2=(A2_t(j))^2+suma2;end

X_A1=sqrt(suma1);X_A2=sqrt(suma2);X_A12=abs(X_A1-X_A2);if(X_A12

min_a12_step=X_A12;FAL_W=W;FAL_E=E;endAW=[];BE=[];A1_t=[];A2_t=[];suma1=0;

%求从影子端点坐标得到的方位角

suma2=0;X_A1=0;X_A2=0;X_A12=0;end

endFAL_WFAL_E

（6）问题三模型关于附件3的求解代码

min_a12_step=1000;CX=[];

%T存放当地北京时刻，单位小时

T=[14.7014.7514.8014.8514.9014.9515.0015.0515.1015.1515.2015.2515.3015.3515.4015.4515.5015.5515.6015.6515.70];

%存放影子端点坐标

X=[1.03651.06991.10381.13831.17321.20871.24481.28151.31891.35681.39551.43491.47511.5161.55771.60031.64381.68821.73371.78011.8277];

Y=[0.49730.50290.50850.51420.51980.52550.53110.53680.54260.54830.55410.55980.56570.57150.57740.58330.58920.59520.60130.60740.6135];

%拟合结果：附件二中最小的影长对应时间为tmin=15.2079；附件三中最小的影长对应的时间为tmin=12.7379；E=120-(tmin-12)*15%E2=71.8860E3=108.9315

E=108.9315;forW=-90:0.2:90

%纬度，北纬为正，南纬为负

%赤纬q

%求当地时角（西经）

%求太

forq=-23.4:0.1:23.4

fori=1:1:21

t=abs((T(i)-(120-E)/15-12))*15;

阳赤纬h

h=asind(sind(W)*sind(q)+cosd(W)*cosd(q)*cosd(t));

A2=asind(cosd(q)*sind(t)/cosd(h));%求从太阳高度角和赤纬和

时角得到的方位角A2

A1=atand(X(i)/Y(i));a1=A1;

A1_t=[AW,a1];AW=A1_t;a2=A2;

A2_t=[BE,a2];BE=A2_t;end

forj=1:1:21

suma1=(A1_t(j))^2+suma1;suma2=(A2_t(j))^2+suma2;end

X_A1=sqrt(suma1);X_A2=sqrt(suma2);X_A12=abs(X_A1-X_A2);if(X_A12

min_a12_step=X_A12;FAL_W=W;FAL_E=E;endAW=[];BE=[];A1_t=[];A2_t=[];suma1=0;suma2=0;X_A1=0;X_A2=0;X_A12=0;end

end

%求从影子端点坐标得到的方位角

FAL_WFAL_E

（7）附件2的拟合代码a=[12.6833333312.73333333

12.7833333312.8333333312.88333333

12.9333333312.9833333313.0333333313.0833333313.1333333313.1833333313.2333333313.2833333313.3333333313.3833333313.4333333313.4833333313.5333333313.5833333313.6333333313.68333333];

b=[1.2472562051.222794591.1989214861.1754289641.1524395731.129917471.107835481.0862542061.0650810721.0444462651.0242641261.0046403140.9854909080.9667904940.9485847350.9309278810.913751750.8971090510.8809737620.8654922590.850504468];Cftoolx=9:0.5:18;

y=0.0981.*x.^2-2.985.*x+23.32;plot(x,y,'*')title('拟合曲线')

xlabel('T（时间单位：h）')ylabel('l

（影子长度单位：m）')

（8）附件3的拟合代码

a=[13.1513.213.2513.313.3513.413.4513.513.5513.613.6513.713.7513.813.8513.913.951414.0514.114.15];b=[3.533142184

3.546768029

3.561797643

3.578100715

3.595750783

3.614934283.6354259833.6572182723.6805411153.7051678363.7312780253.7589179113.7880878883.8187010153.8508096193.884585223.9199118283.9568759923.995534794.0357508354.077863059];Cftoolx=9:0.5:18;

y=0.2964.*x.^2-7.551.*x+51.56;plot(x,y,'*')

title('拟合附件3曲线')xlabel('T（时间单位：h）')ylabel('l

（影子长度单位：m）')

（9）将图像转化为二值图像的代码%求灰度信息

imshow(read(obj,1))f1=read(obj,60001);imshow(f1);

gray=rgb2gray(f1);imshow(gray);f2=read(obj,1);

g2=rgb2gray(f2);imshow(g2);t=g2-grayimshow(t);

img=edge(f2,'sobel');（10）模型四的求解代码：

%T存放当地北京时刻，单位小时

T=[14.7014.7514.8014.8514.9014.9515.0015.0515.1015.1515.2015.2515.3015.3515.4015.4515.5015.5515.6015.6515.70];

%l存放不同时刻的影长

l=[1.03651.06991.10381.13831.17321.20871.24481.28151.31891.35681.39551.43491.47511.5161.55771.60031.64381.68821.73371.78011.8277];

L=2;%杆长，单位米

q=22.1;%赤纬由日期(2015\7\13)得到min_step=1000;AW=[];BE=[];h1_t=[];h2_t=[];sumh1=0;sumh2=0;X_h1=0;X_h2=0;X_h12=0;

forW=-90:0.5:90%纬度，北纬为正，南纬为负

forE=-180:0.5:180

fori=1:1:21

h1=acotd(l(i)/L);

当地时角

h2=asind(sind(W)*sind(q)+cosd(W)*cosd(q)*cosd(t));

当地时角、经度、纬度、赤纬求太阳高度角h2

a=h1;

h1_t=[AW,a];AW=h1_t;b=h2;

h2_t=[BE,b];BE=h2_t;

%由

%由影长、杆长求太阳高度角h1

%由当地北京时间、当地经度求

t=abs((T(i)-(120-E)/15-12))*15;

end

forj=1:1:21

sumh1=(h1_t(j))^2+sumh1;sumh2=(h2_t(j))^2+sumh2;end

X_h1=sqrt(sumh1);X_h2=sqrt(sumh2);X_h12=abs(X_h1-X_h2);if(X_h12

min_step=X_h12;Fin_W=W;Fin_E=E;endAW=[];BE=[];h1_t=[];h2_t=[];sumh1=0;sumh2=0;X_h1=0;X_h2=0;X_h12=0;end

endFin_WFin_E

附录二：赤纬DE表

赤纬DE表（单位：度）（经度=120度，1992年）（12时0分）

日期平年闰年[***********][***********][**************]31

[***********][***********][**************]31

1月-23.1-23.1-23.0-22.9-22.8-22.7-22.6-22.5-22.4-22.2-22.1-21.9-21.8-21.6-21.5-21.3-21.1-20.9-20.7-20.5-20.3-20.1-19.9-19.7-19.4-19.2-18.9-18.7-18.4-18.2-17.9

2月-17.6-17.3-17.1-16.8-16.5-16.2-15.9-15.6-15.3-15.0-14.6-14.3-14.0-13.6-13.3-13.0-12.6-12.3-11.9-11.6-11.2-10.9-10.5-10.1-9.8-9.4-9.0-8.7-8.3

3月-8.3-7.9-7.5-7.1-6.8-6.4-6.0-5.6-5.2-4.8-4.4-4.0-3.6-3.2-2.8-2.5-2.1-1.7-1.3-0.9-0.5-0.10.30.71.11.51.92.32.73.13.53.8

4月3.84.24.65.05.45.86.16.5-6.97.97.68.08.48.79.19.59.810.2-10.510.911.211.611.912.312.612.913.213.613.914.214.5

5月14.514.815.115.415.716.016.316.516.917.117.417.717.918.218.418.718.919.119.419.619.820.020.220.420.620.821.021.221.321.521.6

6月21.821.922.122.222.322.522.622.722.822.922.923.023.123.223.223.323.223.423.423.423.423.423.423.423.423.423.423.423.323.323.3

7月23.223.223.123.022.922.922.822.722.622.522.322.222.121.921.821.621.521.321.221.020.820.620.420.220.019.819.619.419.218.918.718.5

8月18.518.218.017.717.417.216.916.616.416.115.815.515.214.914.614.314.013.713.313.012.712.412.011.711.411.010.710.310.09.69.38.9

9月8.98.68.27.87.57.16.76.46.05.65.24.94.54.13.73.32.92.62.21.81.41.00.60.2-0.2-0.5-0.9-1.3-1.7-2.1-2.5

10月-2.5-2.9-3.3-3.6-4.0-4.4-4.8-5.2-5.6-6.0-6.3-6.7-7.1-7.5-7.8-8.2-8.6-9.0-9.3-9.7-10.1-10.4-10.8-11.1-11.5-11.8-12.2-12.5-12.9-13.2-13.513.9

11月-13.9-14.2-14.5-14.8-15.1-15.5-15.8-16.1-16.4-16.6-16.9-17.217.5-17.8-18.0-18.3-18.6-18.8-19.1-19.3-19.5-19.8-20.020.2-20.4-20.6-20.8-21.0-21.2-21.4-21.5

12月-21.5-21.7-21.8-22.0-22.1-22.3-22.4-22.5-22.6-22.7-22.8-22.9-23.0-23.1-23.2-23.2-23.3-23.3-23.4-23.4-23.4-23.4-23.4-23.4-23.4-23.4-23.4-23.4-23.3-23.3-23.2

注：⑴用月份、日期查表，闰年1、2月份与平年同，从3月1日开始查闰年一行。

⑵一般情况（即不符合1992年、12时、120°E）查此表时，最大误差不大于0.03°。