2015国赛二等奖论文

太阳影子定位

摘要

太阳影子的变化以及地球太阳的对应变化规律广泛应用于定点、测时、测距等一系列操作中，本文通过分析物体太阳影子的变化来确定物体所处的地点及观测日期，从而实现太阳影子定位的功能。

针对问题一：通过分析直杆长度、太阳时角、赤纬角、真太阳时和日期序号等参数对影子长度的影响，建立影长变化模型，得到直杆的影子长度的变化与直杆长度、直杆所处地点的地理纬度、观测时间和日期相关，并得到相应的变化曲线（图4~8）。对于在北京时间9:00至15:00时间范围内3米高直杆进行计算并绘图，得到太阳影子长度的变化曲线图，是一个开口向上、最低点对应的时间在12时左侧的二次曲线（图9）。

针对问题二：通过最小二乘法对直杆影长和观测时间进行二次拟合，得到观测地点的真太阳时，计算出经度，在模型一的基础上建立地点确定模型，求出直杆的纬度，从而确定直杆所处的位置。根据附件1所提供的数据求得直杆位于海南（北纬19度30分50秒，东经110度53分18秒）和印度尼西亚（南纬2度52分26秒，东经110度53分18秒）。

针对问题三：通过分析经纬度、杆长和赤纬角对影长的影响，建立非线性回归模型，利用迭代的方法得到最优解，将该模型应用于附件2和附件3提供的影子顶点数据，得到附件2是长度为2米的直杆在日期为5月26或者7月21日，位于新疆（北纬39度53分31秒，东经79度46分26秒）；得到附件3是长度为3米的直杆，在日期为1月9日或者11月5日，位于湖北（北纬32度50分49秒，东经110度14分42秒）。

针对问题四：对附件4所给的视频每隔2分钟提取出一次灰度图片，采用阈值分割的方法处理图片，计算获得图中直杆影子顶点的坐标，运用问题三的模型，将7月13日作为约束条件，求出视频拍摄地点为内蒙古（北纬41度19分8秒，东经112度52分16秒）。当观测日期未知，去除约束条件，再次求得额视频拍摄地点为内蒙古（北纬41度14分56秒，东经111度42分29秒），拍摄日期为7月16日。

关键词：影子定位，最小二乘法，影长变化模型，地点确定模型，非线性回归模型

1.问题重述

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

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

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

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

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

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

2.问题分析

2.1问题一

题目要求建立影子长度变化的数学模型，且分析影子长度关于各个参数的变化规律。需要考虑影响直杆影子长度变化的多个因素，如地点和时间等。通过查找文献和整理，列出太阳高度角、赤纬角、太阳时角、真太阳时等方程式，以此建立模型，画图并分析模型中各参数对影子长度变化产生的影响。再通过所建立的模型，根据问题一提供的具体时间、地理位置和日期数据，求解出3米高的直杆的太阳影子长度的变化曲线。

2.2问题二

题目要求根据直杆的顶点坐标数据，建立数学模型确定直杆所处的地点。问题二在问题一的基础上已知直杆在某时间范围内影子顶点的坐标，转而求解直杆所处地点的经纬度，且直杆长度未知。首先通过拟合的方法确定直杆影长达到最短的时刻，进而确定经度，然后在问题一的基础上建立模型，通过已知的各个参数求出直杆所处地点的经纬度。

2.3问题三

题目要求根据直杆的顶点坐标数据，建立数学模型确定直杆所处的地点和日期。在问题二的基础上又多了一个未知参数日期，杆长仍然是一个未知的常数。该问题可以通过非线性回归的方法，建立模型，根据误差的大小，计算出符合题目要求的最优解，即为直杆所处地点的经纬度和观测日期。

2.4问题四

先对视频进行处理，计算得到直杆影长的坐标的多组数据，运用问题三的模型，通过约束日期和不约束日期分别对数据进行求解，得到两组答案，通过对比，判断可得到直杆所处地点的经纬度并评判所建模型的准确性。

3.模型假设

（1）假设地球的自转是匀速的。（2）假设太阳光是平行光。

（3）假设忽略所测各点的海拔高度对直杆影长产生的影响。（4）假设忽略太阳光通过大气层的折射对影子长度产生的影响。（5）假设测试当天的天气情况良好，天气不对直杆影子的测量产生影响。（6）假设太阳系是近似球体。

（7）假设忽略直杆底端的位置视频中直杆所处的位置即为摄影机拍摄的位置。

4.符号说明

tEωδ

北京时间

地球绕日公转时运动和转速变化而产生的时差真太阳时太阳时角太阳赤纬角直杆所处地点的纬度直杆所处地点的经度制定标准时间地区的经度

直杆所处地点经度和制定标准时间地区经度的差值太阳高度角直杆的长度

某固定时刻直杆的影长日期序号

任意x（或y）时刻的直杆影长任意x（或y）时刻的太阳高度角任意x（或y）时刻的太阳时角

LLsLcH0HnHx、Hyhx、hy

x、y

5.模型的建立与求解

5.1问题一

5.1.1模型的建立

地球绕太阳运行，采用相对运动的观点，以地球为参考点，可以看做是太阳的移动。某立于地面的直杆在太阳的照射下产生影子，影子的长度会随太阳的移动而产生变化，直杆影子轨

迹形成图如下：

图1直杆影子轨迹线形成图

建立影长变化模型，根据太阳高度角和勾股定理计算得到直杆的影子长度。首先计算太阳高度角h，根据参考文献[1]，得到太阳高度角的正弦值计算公式如下：

sinhsinsincoscoscos

（1）

式中，是直杆所处地点的纬度，是太阳赤纬角，其变化范围是23o27，是太阳时角，赤纬角和时角的位置如图2：

赤纬圈赤道面

地

球

太阳





时圈

图2赤纬角和太阳时角图

根据参考文献[2]，得到太阳赤纬角的计算公式如下：

23.45sin（360

284n

365

（2）

式中，n是日期序号，例如，当日期为1月1日时，n=1，以此类推，当日期为10月22日时，n=295。

根据参考文献[1]，得到太阳时角的计算公式如下：

（T12）15

式中，T是真太阳时（单位：min）。

（3）

太阳时与各地使用的标准时间并不一致，根据参考文献[1]，得到真太阳时T的计算公式如下：

Tt4（LLs）Et4LcE（4）

式中，t是北京时间（即东八区区时），L是直杆所处地点的经度，Ls是制定标准时间地区的经度，Lc是直杆所处地点经度和制定标准时间地区经度的差值（东半球取正，西半球取负），E是地球绕日公转时运动和转速变化而产生的时差（单位：min）。根据参考文献[3]，得到时差E的计算公式如下：

E9.87sin2-7.53cos-1.5sin

360o（n81）

364

根据上述公式计算，可以求出太阳高度角的正弦值sinh，再由勾股定理，如图3

：

（5）（6）

图3勾股定理求影长图

计算公式如下：

sinh

即可求得直杆的影长H。5.1.2模型的求解

H0H0H

（7）

分别以任意随机地点（东经115度，北纬40度）、任意日期（2015年5月25日）、任意杆长3米为例，通过MATLAB绘图，得到以下图片，分析影子长度随各参数的的变化规律：

（1）由图4可知，直杆影长与直杆本身的长度成正比，即直杆越长，影子越长。

影长/m

1234

5 杆长/m

678910

图4影长随日期序号的变化

（2）以3米的直杆为例，由图5可知，季节（年，月，日）的变化会引起日期序号n的改变。当处于夏至日时，太阳直射北回归线，直杆的影长达到最大值，当处于冬至日时，太阳直射南回归线，直杆的影长达到最小值，且杆长的变化以一年为周期，在最大值与最小值之间来回波动。

影长/m

日期序号/天

图5影长随日期序号的变化

（3）由图6可知，直杆的影长随直杆所处位置的纬度的变化不是正比或者反比，在某一固定的纬度，达到影长的最小值，并以该纬度为中心，向两边扩散时，影长变大。

影长/m

510

20253035 直杆所处地点的纬度/度

404550

图6影长随纬度的变化

（4）由图7可知，直杆的的影长在某一固定经度达到最小值，并以该经度为中心，向两边扩散时影长增加。

影长/m

直杆所处位置的经度/度

图7影长随经度的变化

（5）由图8可知，观测时间越靠近直杆所在地的真太阳时，直杆影长越短，越偏离真太阳

时，直杆影长越长。

影长/m

观测时间/n

图8影长随观测时间

综上可知，直杆的影子长度与直杆长度、直杆所处地点的地理经纬度、观测时间和日期均有关系。

针对2015年10月22日北京时间9:00-15:00之间天安门广场（北纬39度54分26秒,东经116度23分29秒）3米高的直杆的太阳影子长度的变化，通过MATLAB计算并画图，求出在北京时间9:00至15:00时间范围内，得到直杆影长随时间变化图及整点时刻的影子长度的数值，见表1和图9：

表1

北京时间（单位：h）直杆影子长度（单位：m）

影长（单位：m）

91011

1213

时间（单位：h）

1415

图9直杆的影长随时间变化图

根据参考文献[4]可知，具有时间一致对应关系的阴影轨迹，如果投射到水平平面上将得到一个二次曲线。再结合图4分析，可得出结论：直杆影长随时间的变化图的是一个开口向上、最低点对应的时间在12时左侧的二次曲线。

5.2问题二

5.2.1模型的建立

首先，采用MABLAB根据直杆影子的顶点坐标算出各时刻所对应的影子长度，由问题一的结论可知影子变化轨迹是一条二次曲线，将影长和对应的北京时间通过最小二乘法进行曲线的二次拟合，得到与原图像拟合度高的图像后，计算图像的顶点坐标，再由该坐标确定直杆影长达到最短的北京时间，通过公式（4）计算，可以得到直杆所处地点的经度。

其次，在问题一建立的影长变化模型基础上，建立地点确定模型，由于直杆处于同一位置且不同时刻的直杆长度不变，根据勾股定理和图3，可知算式如下：

H0HxtanhxHytanhy

（8）

将任意两个时刻得到的公式（8）联立，将杆长H0约去，再与公式（1）联立，可以得到下式：

sinhysinsincoscoscosy

（9）

sin2hy（sinsincoscoscosy）Hxtanhy



xHytanhxx

sin2hx（sinsincoscoscosx）

式中，Hx、Hy是在任意x或y时刻直杆的影长，hx、hy是在任意x或y时刻的太阳高度角，x、y是在任意x或y时刻的太阳时角，且xy。

由于经度通过计算已得到结果，根据观测日期、时间和任意时刻太阳影子顶点坐标数据，

用MATLAB计算，可求出直杆所处的地点的纬度，即得到直杆所处的地点。5.2.2模型的求解

采用MATLAB算法，根据附件1提供的21组x与y坐标数据，算出各个时刻所对应的直杆影长，见表2:

表2附录1中个时刻所对应的直杆影长

北京时间直杆影长/m北京时间直杆影长/m

14:421.149615:151.5402

14:451.182215:181.5799

14:481.215315:211.6201

14:511.249115:241.6613

14:541.283215:271.7033

14:57

1.318015:301.7462

15:001.353415:331.7901

15:031.389415:361.8350

15:061.426215:391.8809

15:091.463415:421.9279

15:121.5015

再用最小二乘法对所有时刻及其对应影长进行曲线的二次拟合，拟合后的误差为0.016，拟合度较高，与原图对比，得到拟合图像如下：

图10二次拟合图像

同时由MATLAB计算得到该曲线的三个参数0.1489，-3.7519和24.1275，可得到影长Hx

与北京时间t的方程：

Hx0.1489t23.7519t24.1275

求得方程顶点坐标（Hx，t）=（0.4929，12.5987），即当天直杆的最短影长为0.4929m，达到最短影长的北京时间是12:35:56，即该地的真太阳时，将此时刻带入问题一所建立的模型里，通过MATLAB计算，得到直杆所处地方的经度L为110.8883度，即东经110度53分18秒。

利用问题二建立的地点确定模型，经过MATLAB计算得到两个纬度1=19.514度，即北纬19度30分50秒，2=-2.874度，即南纬2度52分26秒，从而确定直杆所处地点为（北纬19度30分50秒，东经110度53分18秒）或（南纬2度52分26秒，东经110度53分18秒），如图11、图12所示：

图11直杆所处地点

图12直杆所处地点2

5.3问题三

5.3.1模型的建立

在问题一的基础上，分析直杆所处经纬度、杆长和观测日期对影长产生的影响，建立非线性回归模型：

H0sinhx

tanhxHxsin2hx

（9）

sinhxsinsincoscoscos

2(284n)2(284n)LLs

sinsin(23.45)coscos(23.45cos[15(t12)]

36536515

（10）

通过SPSS编程，采用线性回归的方法处理直杆影子顶点的坐标，求出拍摄位置。

对于所建模型，采用Levenberg-Marquard（简称LM）算法[5]进行对所求结果进行优化，能够求得较为准确的拍摄位置。5.3.2模型的求解

首先将附录二的多组数据导入SPSS，进行非线性回归分析，设置所求变量参数，其中x为纬度，y为经度，n为赤纬角，l为杆长，建立一个非线性规划模型关于x，y，n，l的方程式，通过编程对导入的多组数据进行迭代，得到最优解，如图13

所示：

图13附件二的经纬度求解

计算误差如图14

所示：

图14附件二的误差分析

由图14可知，计算误差极小，计算结果较为准确。

因此，由图13结果可知，在附件二的条件下，长度为2.030米的直杆位于东经79.744度（即79度44分38秒），北纬39.892度（即39度53分31秒），见图15。通过公式（2）和计算所得的赤纬角20.76，计算出日期序号为146.1719或202.3757天，则观测日期为5月26日或7月21日。

图15附件二的位置图

与附件二的计算方法相同，求得在附件三的条件下，是长度为3.036米的直杆，位于东经110.245度（即110度14分52秒），北纬32.847度（即32度50分49秒），日期序号为39.5384天或309.0172天，则观测日期为1月9日或11月5日。经纬度求解、误差分析及位置图如图16、17、18

：

图16

附件三的经纬度求解

图17附件三的误差分析

图18附件三的位置图

5.4问题四

5.4.1已知日期确定位置

在问题三的基础上，用MATLAB编程处理视频，附件4的视频可以分解为60000帧图片，每3000帧提取一张图片，共获得20张图片，例如图

19：

图198：55：00的灰度图像

通过阈值分割方法，得更为清晰的分割图像，如图20

：

图208：55：00的分割图像

运用MATLAB计算，得到影长关于时间变化的表格如下：

表3影长随时间的变化表

时间

8:55

8:57

8:59

9:01

9:03

9:05

9:07

9:09

9:11

9:13

9:15

9:17

9:19

9:21

9:23

9:25

9:27

9:29

9:31

9:33时间（10进制）影长8.91672.35228.952.31718.98332.28489.01672.25559.052.21739.08332.20559.11672.17629.152.14389.18332.11459.21672.0919.252.06469.28332.0449.31672.02889.351.99239.38331.96299.41671.93269.451.90039.48331.87989.51671.85049.551.8329

利用问题三的非线性规划模型，求解出附件四的日期序号为194，则观测日期为7月13日，由计算得知，赤纬角为21.826，用该条件对非线性规划方程进行约束，运用SPSS编程，得到如

下图结果：

图21附件四的经纬度求解—已知日期

可知视频拍摄地点为北纬41.319度（即41度19分8秒），东经112.871度（即112度52分16秒），拍摄位置如图22

：

图22附件四的位置图—未知日期

5.4.2未知日期确定位置

当拍摄视频的日期未知时，SPSS

求解结果如下：

图23附件四的经纬度求解—未知日期

由图23可知，视频拍摄地点为北纬41.249度（即41度14分56秒），东经111.708度（即111度42分29秒），拍摄位置位置如图24。赤纬角为21.579度，计算得到日期序列为197.64天，即拍摄日期为7月16日。

图24附件四的位置图—未知日期

6.模型评价

优点：

（1）运用非线性规划模型拟合分析多元未知参数时，更加简单方便；

（2）得到的数据较为准确，误差较小。

（3）通过几组简单的数据就可以确定模型中的较多参数，在现实中的应用较为广泛，实用性很强。

（4）非线性规划的方法在解决很多统计方面的问题都有很大的作用，有助于最优解的求解。

（5）模型较为易懂，在建立模型的过程中尽量保证了公式的简洁与明了，使模型更加清晰。缺点：

（1）未知量过多时，模型的精确度就会受到影响，模型就不能很好地进行判断，如第四问未知日期时计算所得的日期与视频提供的日期不吻合。

（2）拟合出的误差最小的为最优解，全局最优解只有一个，需要其他误差较小的地点来最终判断最终的结论。

（3）模型忽略了时差以及大气折射率的影响，可能会对结论产生误差。

7.参考文献

[1]中国气象局,地面气象观测规范[M],北京：气象出版社，2003：P133

[2]http://baike.baidu.com/link?url=SrlnklOPjmxXHYj4vGhC833-RRdv307s-WGRvT6N-tU91TR9_24rSAycbPLEiqjqn4JSr-tqeKjlXVZeOP4FYq，2015.09.11

[3]太阳能热利用原理与技术，

http://wenku.baidu.com/link?url=1vxj5fjD25jvTEViMZ2Fal1EaKEJyPSNtwEEjL7o16ljPti9qcvUFnkXRV9UKIadqcIFqOOMBHAhywWFuEryoiTMyxZbCDOd_f_dsTpnPGW，2015.09.11）

[4]武琳，基于太阳阴影轨迹的经纬度估计技术研究[D]，2010：P27

[5]张长胜，一种基于遗传算法和LM算法的混合学习算法[N]，吉林大学学报（理学版），46

（3）：P676，2008

附录：

附录一

使用软件：MATLABSPSS

附录二

程序源代码：

问题一：

%参数-杆长clc;clear;

%t=9:0.01:15%时间g=0:0.01:10;t=12;

n=145;%日期转换成天数

a1=23.45*sin(2*pi*(284+n)/365);%赤纬角B=2*pi*(n-81)/364;

e=9.87*sin(2*B)-7.53*cos(B)-1.5*sin(B);%时差e

H=t-((120-116.39)/15)+e/60;%真太阳时公式a2=(H-12)*15;

a=sind(a1)*sind(40)+cosd(40)*cosd(a1)*cosd(a2);%sinh的值

b=a./((1-a.^2).^0.5);%tanh的值

h=g./b;%影长plot(g,h);

xlabel('杆长/m');ylabel('影长/m');

%参数-天数clc;clear;

%t=9:0.01:15%时间t=12;

n=1:730;%日期转换成天数

a1=23.45*sin(2*pi*(284+n)/365);%赤纬角B=2*pi*(n-81)/364;

e=9.87*sin(2*B)-7.53*cos(B)-1.5*sin(B);%时差e

H=t-((120-110)/15)+e/60;%真太阳时公式a2=(H-12)*15;

a=sind(a1).*sind(40)+cosd(40)*cosd(a1).*cosd(a2);%sinh的值

b=a./((1-a.^2).^0.5);%tanh的值

h=3./b;%影长plot(n,h);

xlabel('日期序号/n');ylabel('影长/m');

%参数-时间clc;clear;

t=8:0.01:16;%时间%t=12;

n=145;%日期转换成天数

a1=23.45*sin(2*pi*(284+n)/365);%赤纬角B=2*pi*(n-81)/364;

e=9.87*sin(2*B)-7.53*cos(B)-1.5*sin(B);%时差e

H=t-((120-115)/15)+e/60;%真太阳时公式a2=(H-12)*15;

a=sind(a1)*sind(40)+cosd(40)*cosd(a1)*cosd(a2);%sinh的值

b=a./((1-a.^2).^0.5);%tanh的值

h=3./b;%影长plot(t,h);

xlabel('观测时间/n');ylabel('影长/m');

%参数-纬度clc;clear;

%t=8:0.01:16;%时间t=12;ww=5:50;

n=145;%日期转换成天数

a1=23.45*sin(2*pi*(284+n)/365);%赤纬角B=2*pi*(n-81)/364;

e=9.87*sin(2*B)-7.53*cos(B)-1.5*sin(B);%时差e

H=t-((120-115)/15)+e/60;%真太阳时公式a2=(H-12)*15;

a=sind(a1).*sind(ww)+cosd(ww).*cosd(a1)*cosd(a2);%sinh的值

b=a./((1-a.^2).^0.5);%tanh的值

h=3./b;%影长plot(ww,h);

xlabel('直杆所处地点的纬度/度');ylabel('影长/m');

%参数-经度clc;clear;

%t=8:0.01:16;%时间t=12;jj=70:130

n=145;%日期转换成天数

a1=23.45*sin(2*pi*(284+n)/365);%赤纬角B=2*pi*(n-81)/364;

e=9.87*sin(2*B)-7.53*cos(B)-1.5*sin(B);%时差e

H=t-((120-jj)./15)+e/60;%真太阳时公式a2=(H-12)*15;

a=sind(a1)*sind(40)+cosd(40)*cosd(a1).*cosd(a2);%sinh的值

b=a./((1-a.^2).^0.5);%tanh的值

h=3./b;%影长plot(jj,h);

xlabel('直杆所处位置的经度/度');ylabel('影长/m');

%3米杆的影长随时间变化曲线clc;clear;

t=9:0.01:15%时间

n=295;%日期转换成天数

a1=23.45*sin(2*pi*(284+n)/365);%赤纬角B=2*pi*(n-81)/364;

e=9.87*sin(2*B)-7.53*cos(B)-1.5*sin(B);%时差e

T=t-((120-116.3913889)/15)+e/60;%真太阳时公式a2=(T-12)*15;%时角

a=sind(a1)*sind(39.907222)+cosd(39.907)*cosd(a1)*cosd(a2);%太阳高度角h的sinh的值

b=a./((1-a.^2).^0.5);%tanh的值

H=3./b;%影长Hplot(t,H);

%title('直杆的影长随时间的变化');

xlabel('时间（单位：h）');ylabel('影长（单位：m）');

问题2：

（1）最小二乘法拟合

clc;clearall;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];t=[14.714.7514.814.8514.914.951515.0515.115.1515.215.2515.315.3515.415.4515.515.5515.615.6515.7];t1=9:0.01:16;%z=ones(1,21);z=(x.^2+y.^2).^0.5;%plot(t,z,'r*');%func=@(a,t)a(1)*(1-(sin(a(2))*sin(10.51)+cos(t+(a(3)-120)/15+0.523-12).*cos(a(2)).*cos(10.51)).^2).^0.5/cos(t+(a(3)-120)/15+0.523-12).*sin(a(2).*sin(10.51)+cos(a(2)*cos(10.51)));%A=lsqcurvefit(func,[1,1,1],t,z);p=polyfit(t,z,2);y1=polyval(p,t);y2=polyval(p,t1);e=y1-z;e1=abs(e);sum(e1);holdonplot(t1,y2);plot(t,z,'*');

xlabel('时间/t');ylabel('影长/m');holdoff

（2）求经度方程

clc;clear;symstx;

z=0.1489*t^2-3.7519*t+24.1275;%拟合方程%ezplot(z,[8,18]);a=3.7519/0.1489/2;

b=0.1489*a^2-3.7519*a+24.1275;%影长最小的时间

n=108;%日期转换成天数

a1=23.45*sin(2*pi*(284+n)/365);%赤纬角B=2*pi*(n-81)/364;

e=9.87*sin(2*B)-7.53*cos(B)-1.5*sin(B);%时差e

x=15*(12-a-e/60)+120;%该点的精度

d=((120-x)/15)+e/60;%求差值

（3）求解纬度方程

clcclearall;%symsx;%t=solve(1.9279*1.9279*(1/(sin(x)*sind(10.511)+cos(x)*cosd(10.511)*cosd((14.7-12.6162)*15))^2-1)-1.1496*1.1496*(1/(sin(x)*sind(10.511)+cos(x)*cosd(10.511)*cosd((15.7-12.6162)*15))^2-1))%t=solve(1.9279*(1/(sin(x)*0.1824+cos(x)*0.9832*cosd((14.7-12.6162)*15))^2-1)-1.1496*(1/(sin(x)*0.1824+cos(x)*0.9832*cosd((15.7-12.6162)*15))^2-1))%t=double(t)a=1.92792;b=1.66127;c=15.4;d=15.7;f=@(x)(a*a*(1/(sin(x)*sind(10.511)+cos(x)*cosd(10.511)*cosd((c-12.6162)*15))^2-1)-b*b*(1/(sin(x)*sind(10.511)+cos(x)*cosd(10.511)*cosd((d-12.6162)*15))^2-1));x=fsolve(f,0);

x1=x*180/pi

（4）Spss方程实例

l*COS(ARSIN(SIN(a*0.0174532)*SIN(n*0.0174532)+COS(a*0.0174532)*COS(n*0.0174532)*COS((15*t-300+b+0.1309)*0.0174532)))/SIN(ARSIN(SIN(a*0.0174532)*SIN(n*0.0174532)+COS(a*0.0174532)*COS(n*0.0174532)*COS((15*t-300+b+0.1309)*0.0174532)))

问题3：

Spss方程实例

l*COS(ARSIN(SIN(x*0.0174532)*SIN(n*0.0174532)+COS(x*0.0174532)*COS(n*0.0174532)*COS((15*t-300+y)*0.0174532)))/SIN(ARSIN(SIN(x*0.0174532)*SIN(n*0.0174532)+COS(x*0.0174532)*COS(n*0.0174532)*COS((15*t-300+y)*0.0174532)))

问题4：

（1）读取灰度图片程序

clc;clear;

video=VideoReader('Appendix4.avi');%读取图片

nFrames=video.NumberOfFrames;%帧的总数

vidHeight=video.Height;%得倒图片的高

vidWidth=video.Width;%得倒图片的宽fork=1350:3000:nFrames

Im=read(video,k);%读出第k帧I=rgb2gray(Im);

imshow(I);%显示每一帧

imwrite(I,strcat(num2str(k),'.jpg'),'jpg');%保存帧end

（2）灰度图片处理程序

clcclearBW=imread('1350.jpg');%BW=rgb2gray(BW);BW=im2bw(BW,0.8);%BW=edge(BW,'sobel');BW=~BW;BW=bwareaopen(BW,34600,8);

L=bwlabel(BW,4);imshow(BW);%imshow(L);STATS=regionprops(L,'MajorAxisLength');

（3）计算影长程序

clccleara=1497;b=886;c=((a-872)^2+(884-b)^2)^0.5;d=c*2/682

（4）spss方程实例

2*COS(ARSIN(SIN(x*0.0174532)*SIN(n*0.0174532)+COS(x*0.0174532)*COS(n*0.0174532)*COS((15*t-300+y)*0.0174532)))/SIN(ARSIN(SIN(x*0.0174532)*SIN(n*0.0174532)+COS(x*0.0174532)*COS(n*0.0174532)*COS((15*t-300+y)*0.0174532)))

2*COS(ARSIN(SIN(a*0.0174532)*SIN(n*0.0174532)+COS(a*0.0174532)*COS(n*0.0174532)*COS((15*t-300+b-1.3495)*0.0174532)))/SIN(ARSIN(SIN(a*0.0174532)*SIN(n*0.0174532)+COS(a*0.0174532)*COS(n*0.0174532)*COS((15*t-300+b-1.3495)*0.0174532)))

太阳影子定位

摘要

关键词：影子定位，最小二乘法，影长变化模型，地点确定模型，非线性回归模型

1.问题重述

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

2.问题分析

2.1问题一

2.2问题二

2.3问题三

2.4问题四

3.模型假设

（1）假设地球的自转是匀速的。（2）假设太阳光是平行光。

（7）假设忽略直杆底端的位置视频中直杆所处的位置即为摄影机拍摄的位置。

4.符号说明

tEωδ

北京时间

地球绕日公转时运动和转速变化而产生的时差真太阳时太阳时角太阳赤纬角直杆所处地点的纬度直杆所处地点的经度制定标准时间地区的经度

直杆所处地点经度和制定标准时间地区经度的差值太阳高度角直杆的长度

某固定时刻直杆的影长日期序号

任意x（或y）时刻的直杆影长任意x（或y）时刻的太阳高度角任意x（或y）时刻的太阳时角

LLsLcH0HnHx、Hyhx、hy

x、y

5.模型的建立与求解

5.1问题一

5.1.1模型的建立

迹形成图如下：

图1直杆影子轨迹线形成图

sinhsinsincoscoscos

（1）

式中，是直杆所处地点的纬度，是太阳赤纬角，其变化范围是23o27，是太阳时角，赤纬角和时角的位置如图2：

赤纬圈赤道面

地

球

太阳





时圈

图2赤纬角和太阳时角图

根据参考文献[2]，得到太阳赤纬角的计算公式如下：

23.45sin（360

284n

365

（2）

式中，n是日期序号，例如，当日期为1月1日时，n=1，以此类推，当日期为10月22日时，n=295。

根据参考文献[1]，得到太阳时角的计算公式如下：

（T12）15

式中，T是真太阳时（单位：min）。

（3）

太阳时与各地使用的标准时间并不一致，根据参考文献[1]，得到真太阳时T的计算公式如下：

Tt4（LLs）Et4LcE（4）

E9.87sin2-7.53cos-1.5sin

360o（n81）

364

根据上述公式计算，可以求出太阳高度角的正弦值sinh，再由勾股定理，如图3

：

（5）（6）

图3勾股定理求影长图

计算公式如下：

sinh

即可求得直杆的影长H。5.1.2模型的求解

H0H0H

（7）

（1）由图4可知，直杆影长与直杆本身的长度成正比，即直杆越长，影子越长。

影长/m

1234

5 杆长/m

678910

图4影长随日期序号的变化

影长/m

日期序号/天

图5影长随日期序号的变化

影长/m

510

20253035 直杆所处地点的纬度/度

404550

图6影长随纬度的变化

（4）由图7可知，直杆的的影长在某一固定经度达到最小值，并以该经度为中心，向两边扩散时影长增加。

影长/m

直杆所处位置的经度/度

图7影长随经度的变化

（5）由图8可知，观测时间越靠近直杆所在地的真太阳时，直杆影长越短，越偏离真太阳

时，直杆影长越长。

影长/m

观测时间/n

图8影长随观测时间

综上可知，直杆的影子长度与直杆长度、直杆所处地点的地理经纬度、观测时间和日期均有关系。

表1

北京时间（单位：h）直杆影子长度（单位：m）

影长（单位：m）

91011

1213

时间（单位：h）

1415

图9直杆的影长随时间变化图

5.2问题二

5.2.1模型的建立

H0HxtanhxHytanhy

（8）

将任意两个时刻得到的公式（8）联立，将杆长H0约去，再与公式（1）联立，可以得到下式：

sinhysinsincoscoscosy

（9）

sin2hy（sinsincoscoscosy）Hxtanhy



xHytanhxx

sin2hx（sinsincoscoscosx）

式中，Hx、Hy是在任意x或y时刻直杆的影长，hx、hy是在任意x或y时刻的太阳高度角，x、y是在任意x或y时刻的太阳时角，且xy。

由于经度通过计算已得到结果，根据观测日期、时间和任意时刻太阳影子顶点坐标数据，

用MATLAB计算，可求出直杆所处的地点的纬度，即得到直杆所处的地点。5.2.2模型的求解

采用MATLAB算法，根据附件1提供的21组x与y坐标数据，算出各个时刻所对应的直杆影长，见表2:

表2附录1中个时刻所对应的直杆影长

北京时间直杆影长/m北京时间直杆影长/m

14:421.149615:151.5402

14:451.182215:181.5799

14:481.215315:211.6201

14:511.249115:241.6613

14:541.283215:271.7033

14:57

1.318015:301.7462

15:001.353415:331.7901

15:031.389415:361.8350

15:061.426215:391.8809

15:091.463415:421.9279

15:121.5015

再用最小二乘法对所有时刻及其对应影长进行曲线的二次拟合，拟合后的误差为0.016，拟合度较高，与原图对比，得到拟合图像如下：

图10二次拟合图像

同时由MATLAB计算得到该曲线的三个参数0.1489，-3.7519和24.1275，可得到影长Hx

与北京时间t的方程：

Hx0.1489t23.7519t24.1275

图11直杆所处地点

图12直杆所处地点2

5.3问题三

5.3.1模型的建立

在问题一的基础上，分析直杆所处经纬度、杆长和观测日期对影长产生的影响，建立非线性回归模型：

H0sinhx

tanhxHxsin2hx

（9）

sinhxsinsincoscoscos

2(284n)2(284n)LLs

sinsin(23.45)coscos(23.45cos[15(t12)]

36536515

（10）

通过SPSS编程，采用线性回归的方法处理直杆影子顶点的坐标，求出拍摄位置。

对于所建模型，采用Levenberg-Marquard（简称LM）算法[5]进行对所求结果进行优化，能够求得较为准确的拍摄位置。5.3.2模型的求解

所示：

图13附件二的经纬度求解

计算误差如图14

所示：

图14附件二的误差分析

由图14可知，计算误差极小，计算结果较为准确。

图15附件二的位置图

：

图16

附件三的经纬度求解

图17附件三的误差分析

图18附件三的位置图

5.4问题四

5.4.1已知日期确定位置

在问题三的基础上，用MATLAB编程处理视频，附件4的视频可以分解为60000帧图片，每3000帧提取一张图片，共获得20张图片，例如图

19：

图198：55：00的灰度图像

通过阈值分割方法，得更为清晰的分割图像，如图20

：

图208：55：00的分割图像

运用MATLAB计算，得到影长关于时间变化的表格如下：

表3影长随时间的变化表

时间

8:55

8:57

8:59

9:01

9:03

9:05

9:07

9:09

9:11

9:13

9:15

9:17

9:19

9:21

9:23

9:25

9:27

9:29

9:31

下图结果：

图21附件四的经纬度求解—已知日期

可知视频拍摄地点为北纬41.319度（即41度19分8秒），东经112.871度（即112度52分16秒），拍摄位置如图22

：

图22附件四的位置图—未知日期

5.4.2未知日期确定位置

当拍摄视频的日期未知时，SPSS

求解结果如下：

图23附件四的经纬度求解—未知日期

图24附件四的位置图—未知日期

6.模型评价

优点：

（1）运用非线性规划模型拟合分析多元未知参数时，更加简单方便；

（2）得到的数据较为准确，误差较小。

（3）通过几组简单的数据就可以确定模型中的较多参数，在现实中的应用较为广泛，实用性很强。

（4）非线性规划的方法在解决很多统计方面的问题都有很大的作用，有助于最优解的求解。

（5）模型较为易懂，在建立模型的过程中尽量保证了公式的简洁与明了，使模型更加清晰。缺点：

（1）未知量过多时，模型的精确度就会受到影响，模型就不能很好地进行判断，如第四问未知日期时计算所得的日期与视频提供的日期不吻合。

（2）拟合出的误差最小的为最优解，全局最优解只有一个，需要其他误差较小的地点来最终判断最终的结论。

（3）模型忽略了时差以及大气折射率的影响，可能会对结论产生误差。

7.参考文献

[1]中国气象局,地面气象观测规范[M],北京：气象出版社，2003：P133

[2]http://baike.baidu.com/link?url=SrlnklOPjmxXHYj4vGhC833-RRdv307s-WGRvT6N-tU91TR9_24rSAycbPLEiqjqn4JSr-tqeKjlXVZeOP4FYq，2015.09.11

[3]太阳能热利用原理与技术，

http://wenku.baidu.com/link?url=1vxj5fjD25jvTEViMZ2Fal1EaKEJyPSNtwEEjL7o16ljPti9qcvUFnkXRV9UKIadqcIFqOOMBHAhywWFuEryoiTMyxZbCDOd_f_dsTpnPGW，2015.09.11）

[4]武琳，基于太阳阴影轨迹的经纬度估计技术研究[D]，2010：P27

[5]张长胜，一种基于遗传算法和LM算法的混合学习算法[N]，吉林大学学报（理学版），46

（3）：P676，2008

附录：

附录一

使用软件：MATLABSPSS

附录二

程序源代码：

问题一：

%参数-杆长clc;clear;

%t=9:0.01:15%时间g=0:0.01:10;t=12;

n=145;%日期转换成天数

a1=23.45*sin(2*pi*(284+n)/365);%赤纬角B=2*pi*(n-81)/364;

e=9.87*sin(2*B)-7.53*cos(B)-1.5*sin(B);%时差e

H=t-((120-116.39)/15)+e/60;%真太阳时公式a2=(H-12)*15;

a=sind(a1)*sind(40)+cosd(40)*cosd(a1)*cosd(a2);%sinh的值

b=a./((1-a.^2).^0.5);%tanh的值

h=g./b;%影长plot(g,h);

xlabel('杆长/m');ylabel('影长/m');

%参数-天数clc;clear;

%t=9:0.01:15%时间t=12;

n=1:730;%日期转换成天数

a1=23.45*sin(2*pi*(284+n)/365);%赤纬角B=2*pi*(n-81)/364;

e=9.87*sin(2*B)-7.53*cos(B)-1.5*sin(B);%时差e

H=t-((120-110)/15)+e/60;%真太阳时公式a2=(H-12)*15;

a=sind(a1).*sind(40)+cosd(40)*cosd(a1).*cosd(a2);%sinh的值

b=a./((1-a.^2).^0.5);%tanh的值

h=3./b;%影长plot(n,h);

xlabel('日期序号/n');ylabel('影长/m');

%参数-时间clc;clear;

t=8:0.01:16;%时间%t=12;

n=145;%日期转换成天数

a1=23.45*sin(2*pi*(284+n)/365);%赤纬角B=2*pi*(n-81)/364;

e=9.87*sin(2*B)-7.53*cos(B)-1.5*sin(B);%时差e

H=t-((120-115)/15)+e/60;%真太阳时公式a2=(H-12)*15;

a=sind(a1)*sind(40)+cosd(40)*cosd(a1)*cosd(a2);%sinh的值

b=a./((1-a.^2).^0.5);%tanh的值

h=3./b;%影长plot(t,h);

xlabel('观测时间/n');ylabel('影长/m');

%参数-纬度clc;clear;

%t=8:0.01:16;%时间t=12;ww=5:50;

n=145;%日期转换成天数

a1=23.45*sin(2*pi*(284+n)/365);%赤纬角B=2*pi*(n-81)/364;

e=9.87*sin(2*B)-7.53*cos(B)-1.5*sin(B);%时差e

H=t-((120-115)/15)+e/60;%真太阳时公式a2=(H-12)*15;

a=sind(a1).*sind(ww)+cosd(ww).*cosd(a1)*cosd(a2);%sinh的值

b=a./((1-a.^2).^0.5);%tanh的值

h=3./b;%影长plot(ww,h);

xlabel('直杆所处地点的纬度/度');ylabel('影长/m');

%参数-经度clc;clear;

%t=8:0.01:16;%时间t=12;jj=70:130

n=145;%日期转换成天数

a1=23.45*sin(2*pi*(284+n)/365);%赤纬角B=2*pi*(n-81)/364;

e=9.87*sin(2*B)-7.53*cos(B)-1.5*sin(B);%时差e

H=t-((120-jj)./15)+e/60;%真太阳时公式a2=(H-12)*15;

a=sind(a1)*sind(40)+cosd(40)*cosd(a1).*cosd(a2);%sinh的值

b=a./((1-a.^2).^0.5);%tanh的值

h=3./b;%影长plot(jj,h);

xlabel('直杆所处位置的经度/度');ylabel('影长/m');

%3米杆的影长随时间变化曲线clc;clear;

t=9:0.01:15%时间

n=295;%日期转换成天数

a1=23.45*sin(2*pi*(284+n)/365);%赤纬角B=2*pi*(n-81)/364;

e=9.87*sin(2*B)-7.53*cos(B)-1.5*sin(B);%时差e

T=t-((120-116.3913889)/15)+e/60;%真太阳时公式a2=(T-12)*15;%时角

a=sind(a1)*sind(39.907222)+cosd(39.907)*cosd(a1)*cosd(a2);%太阳高度角h的sinh的值

b=a./((1-a.^2).^0.5);%tanh的值

H=3./b;%影长Hplot(t,H);

%title('直杆的影长随时间的变化');

xlabel('时间（单位：h）');ylabel('影长（单位：m）');

问题2：

（1）最小二乘法拟合

xlabel('时间/t');ylabel('影长/m');holdoff

（2）求经度方程

clc;clear;symstx;

z=0.1489*t^2-3.7519*t+24.1275;%拟合方程%ezplot(z,[8,18]);a=3.7519/0.1489/2;

b=0.1489*a^2-3.7519*a+24.1275;%影长最小的时间

n=108;%日期转换成天数

a1=23.45*sin(2*pi*(284+n)/365);%赤纬角B=2*pi*(n-81)/364;

e=9.87*sin(2*B)-7.53*cos(B)-1.5*sin(B);%时差e

x=15*(12-a-e/60)+120;%该点的精度

d=((120-x)/15)+e/60;%求差值

（3）求解纬度方程

x1=x*180/pi

（4）Spss方程实例

问题3：

Spss方程实例

问题4：

（1）读取灰度图片程序

clc;clear;

video=VideoReader('Appendix4.avi');%读取图片

nFrames=video.NumberOfFrames;%帧的总数

vidHeight=video.Height;%得倒图片的高

vidWidth=video.Width;%得倒图片的宽fork=1350:3000:nFrames

Im=read(video,k);%读出第k帧I=rgb2gray(Im);

imshow(I);%显示每一帧

imwrite(I,strcat(num2str(k),'.jpg'),'jpg');%保存帧end

（2）灰度图片处理程序

clcclearBW=imread('1350.jpg');%BW=rgb2gray(BW);BW=im2bw(BW,0.8);%BW=edge(BW,'sobel');BW=~BW;BW=bwareaopen(BW,34600,8);

L=bwlabel(BW,4);imshow(BW);%imshow(L);STATS=regionprops(L,'MajorAxisLength');

（3）计算影长程序

clccleara=1497;b=886;c=((a-872)^2+(884-b)^2)^0.5;d=c*2/682

（4）spss方程实例

2015国赛二等奖论文

相关文章