太阳影子定位技术江西师范大学科学技术学院商江山_袁星_蒋梦婷
A题太阳影子定位
摘要：正由于太阳光线的照射,物体会形成一个影子,影子的朝向、长短等会随着地球的自转运动位置的不同,而时刻发生变化。面对以此作为载体，通过太阳影子的定位技术，分析视频中物体的太阳影子变化规律，从而确定视频的拍摄地点和日期。
问题一：影子长度模型的建立。本论文运用了太阳高度角、太阳赤纬和太阳时角公式，再结合麦克劳林公式，分析影子长度关于各个参数的变化规律，建立影子长度变化的数学模型。
当时间为2015年10月22日北京时间9:00—15:00时，天安门广场3米高的直杆的影子变化曲线则可以根据地点的经纬度值（北纬39度54分26秒,东经116度23分29秒）和地球的自转速度为465.18m/s绘出。
问题二：直杆位置的确定。首先根据影子顶点的坐标数据，参考解析几何三角形的勾股定理建立数学表达式，通过计算得到对应的影子长度值。得出物体处于北半球，利用matlab求方差公式计算出最小方差值，因此得出参考拍摄位置的东经与北纬的数值，即东经108.6000，北纬19.3000，所以直杆的拍摄地点为海南。
问题三：思路和问题二基本相似，先利用勾股定理求得影长，但由于不知太阳照射的日期，从而无法直接对所求数据进行处理，所以结合交比不变性定理算法原理，采用直线网络模板，使用最小二乘法对数据进行处理，从而估算出附件2与附件3所处的地理位置，它们分别为东经76，北纬33和东经109，北纬46，所以附件2的拍摄地点为新疆喀什，时间为2015年3月19日；附件3的拍摄地点为内蒙古，时间为2015年3月18日。
问题四：根据已知影子顶点的坐标采用Matlab软件对图像进行拟合分析和拟合公式得出影子长度最短所对应的当地时间得出经度。再次把本题的模型建立转化成问题二和问题三的模型。由问题三的纬度模型，把问题四得到影子长度和时间对应关系转化成附件二中的影子顶点坐标，最后把得出的坐标应用于纬度模型从而得到对应的纬度，由以上所求得的经纬度数值，得出视频的拍摄地点与日期。关键字：太阳运动定位技术影子变化规律
太阳时角麦克劳林公式勾股定理
最小二乘法太阳高度角太阳赤纬比不变性定理算法原理
Ⅰ问题重述与分析
1.1问题重述
1、由于地球围绕太阳的自转，每时每刻物体在太阳下的影子长度也是时刻在变化着，分析影子长度关于各个参数的变化规律，建立影子变化模型。那么2015年10月22日北京时间9:00—15:00之间天安门广场（北纬39度54分26秒,东经116度23分29秒）3米高的直杆的太阳影子长度变化曲线是什么样的呢？
2、已知固定2米的直杆在水平地面上太阳影子顶点坐标数据与拍摄时间为2015年4月18日，建立数学模型确定直杆所处位置。
3、根据某2米的固定直杆在水平地面上太阳影子顶点坐标数据，建立数学模型，进而确定附件2和附件3中直杆所处的地点和日期。
4、附件4为一根直杆在太阳下的影子变化的视频，并且已通过某种方式估计出直杆的高度为2米。请建立确定视频拍摄地点的数学模型，并应用你们的模型给出若干个可能的拍摄地点。
如果拍摄日期未知，你能否根据视频确定出拍摄地点与日期？
1.2、问题的分析：
针对问题一：
如图1，是直杆在太阳照
射下不同时间段影子的长短变
化图形，我们可以采用立体几
何的特点得直杆影子长度和太
阳高度角的数学关系式，也就
是通过立体几何的特点建立相
应的数学模型。首先已知题1
中的日期参数，可得赤纬角表达式，再依据时间段的
要求，结合地球自转速度，则可以知道太阳高度角的
变化曲线，最后通过分析地点的经纬度得到方位角变
化曲线，进而画出直杆的太阳影子长度的变化曲线。
针对问题二：
从附件1中可以知道某直杆的在水平地面上的太阳影子顶点坐标数据和拍摄的日期为2015年4月18日，利用三角形勾股定理知识求得对应的影子长度a，如表（一）所述。通过对地球北半球经纬度的扫描，结合Matlab中的方差公式，可求得不同地点的方差值，那么方差越小，位置越准确。又因为杆子是固定的，杆长为2米，所以在不考虑杆子与地面倾斜造成的误差的前提下，得出杆子所处的位置为在东经108.6000，北纬
19.3000，所以直杆的拍摄地点为海南。
针对问题三：
因为只知道直杆在水平地面上的太阳影子顶点坐标数据，却不知拍摄的日期，所以需要先计算出对应的影子长度如表（二）和表（三），然后根据影子长度的变化规律，随着一天中时间的推移，影子的长度越来越短，推算出是平年，再采用matlab中的循环语句得出拍摄的日期。最后重复题二的解题步骤得出拍摄的地点。
得出的结果为：1附件2中直杆的地理位置为76，北纬33，地处新疆喀什，时间为2015年3月19日。2附件3中直杆的地理位置为109，北纬46，地处内蒙古，时间为2015年3月18日。
针对问题四：视频帧拟合技术
从视频中，我们可以大概观察到影子的长度变化越来越慢，这与我们的生活常识相一致，我们的观察能力是有限度的，这需要我们借助计算机分析，从而的到更多的数据。我们利用MATLAB软件，在一定的时间隔内在视频里面抽取视频画面，获取自然图像序列或者视频帧,并且对每一帧图像检测出影子的轨迹点；然后确定多个灭点,并拟合出地平线；拟合互相垂直的灭点,计算出仿射纠正和投影纠正矩阵；进而还原出经过度量纠正的世界坐标；再拟合出经过度量纠正世界坐标中的影子点的轨迹,并且参考日晷设计,利用相似关系估计出纬度；利用二次曲线的极值点计算时差,从而有效地恢复出所拍摄图像的经纬度信息。本发明复杂度较低,精度较高,是一种利用自然图像序列或视频帧实现经纬度估计的方法,能够根据未经过校准的影子的位置,得出所拍摄图像的经纬度信息。
Ⅱ模型的假设与符号的说明
2.1模型的假设
假设1：假设地球的公转对直杆影子的变化不产生影响。
假设2：假设空气中二氧化碳，水蒸气等各种因素对太阳光进入大气层时发生折射的影响忽略不计。
假设3：不考虑地壳板块运动对拍摄影子顶点坐标的影响。
假设4：忽略不计人工测量时，数据采集近似认为是无穷小。
假设5：假设测量期间天气状况平稳，对测量的结构影响甚小。
2.2符号说明
i；经度
j；纬度
A；太阳方位角
a；杆子的影子长度
std；方差
N；从一月一号算起的天数
y；年
m；月
d；日
H；高度角
delta；弧度
k：天数
Ⅲ模型的建立与求解
3.1、模型一的建立：基于给定参数的几何模型
本题需要求解的是影子长度关于各个参数的变化规律。分析题目可知影响影子长度的因素有：杆子所处的位置和太阳的相对位置、杆子的长度等。太阳高度角h、太阳赤纬和太阳时角t是表示太阳和杆子相对位置的参数，计算公式为：太阳高度角h： sinh=sinsin+coscoscost
（：太阳赤纬，：实际地理纬度，t为当地太阳时角，都是弧度制）太阳赤纬：=0.3723+23.2567sin+0.1149sin2-0.1712sin3-
0.758cos+0.3656cos2+0.0201cos3
（N-N0）其中：=365.2422
(N为天数，自1月日起，）1
（ N0=796764+0.2422（年份-1985）-int(（年份-1985）4))
太阳时角t：t(12w)180(120)
其中：w为北京时间，为当地经度
通过公式的计算，建立了影子长度变化的数
学模型一。若要求得天安门前杆子影子变化曲线，
只需要将北京的日期时间，经纬度带入即可得到
方位变化曲线，高度变化曲线以及影子变化曲线。
如图2所示。
3.2、影子长度的求解
运用立体几何知识，以直杆底端为中心点作为原点建立平面直角坐标系，将附件1、附件2和附件3中影子顶点的坐标数据描于坐标系中，运用三角形的勾股定理得到不同坐标下影子的长度值aixiyi，分别制成表（一）、表（二）、2
2
表（三）。
3.3纬度模型的建立
因为不同的节气，太阳的高
度变化规律不同，且一天中由于
黄赤交角的存在，太阳直射点的
南北移动，引起太阳高度的变化。
从题目二中得知拍摄的时间为
2015年4月18日，故可确定太
阳是直射北半球的。再根据太阳
高度角概念经纬度的变化规律估
算出拍摄地点的经度与纬度值如
表（四），
从而确定所处的位置。
表（四）
附件
东经（度）
北纬（度）
地点
时间附件1108.600019.3000海南2015年4月18日附件27633新疆喀什2015年3月19日附件310946内蒙古2015年3月18日变化规律：1纬度变化规律：由于太阳直射点所在经纬度向南北两侧递减。可推知与太阳直射点的纬度相差一度，正午太阳高度角就减小一度。
2
季节变化规律：太阳直射点移来时渐增，移去时渐减（太阳直射点相对某地所在纬线而言）。
图四为附件4视频中不同时间点物体的太阳影子变化图。根据对问题四的分析思路，从而得到物体的东经与北纬的数值。
3.4最小二乘法的曲线拟合
数据拟合的具体做法是：对给定数据xi,yi（i=1,2,
类中，求ripxiyi，使误差（i=1,2,
21）,在取定的函数21）的平方和最小，即
rpxy2
iii
i1i1nn2min
21）的距离平方和为从集合意义上讲，就是寻求与给定点xi,yi（i=1,2,
最小的曲线ypx，函数px称为拟合函数或最小二乘解，求拟合函数px的方法称为曲线拟合的最小二乘法。
根据附件1影子的坐标，根据勾股定理易得影子的长度在一个小时的变化
1.149625826
1.353364049
1.579853316
1.8350142721.1821989761.2152969551.2490510521.283195341.3179931491.3893870911.4261528561.4633998531.5014816221.5402318171.6201445151.6612706131.7032906331.746205911.7900509151.8808750011.927918447
分析这21组数据，我们断定拍摄地点在北半球，我们采用MATLAB经纬度自动搜索经度从0开始步长为1到180为止纬度，纬度从0到90搜索，将每一组得到的影长数据与给定的数据进行拟合，数据相关程度越高，就与拍摄地点越接近。然后缩小范围，得出准确的经纬度坐标。
Ⅳ误差分析及模型评价与改进
4.1误差的分析
在测量中存在许多影响影子的长度因素，针对这些因素，主要从以下几个方面进行考虑。1一是根据经纬度确定杆子与太阳的相对位置和确定太阳的高度角时，由于测量数据的近似值，会造成经纬度值计算的偏差；2二是空气中二氧化碳，水蒸气等各种因素会使太阳光进入大气层时发生折射，那么太阳光线与地平面的交角会与实际的交角有误差；
4.2模型评价与改进
对于本论文中的四个问题分别建立了近似数学模型，在初步计算的数值中，选取经度较准确的几组数据，根据经度与纬度不同取值范围，进行统计分析。通过所得的数据可以明显的看出，在纬度相同的情况下，太阳的高度角则相同。第三问所建立的数学模型是最准确的，从而验证了我们问题求解的合理性结果都比较准确，但是始终存在一些不能排除误差，在此我们认为是系统误差，于此我们
针对问题3中确定拍摄日期问题模型利用蒙特卡洛随机模拟进行改进模型。
我们在三维空间产生随机点，通过已建立的空间坐标系拟合影长与地理位置坐标是否符合。当所求得方差值最小时，则可求得准确的位置经度与纬度。
蒙特卡洛随机模拟如下程序表：
Ⅵ参考文献
[1]《对建筑日照计算中太阳赤纬交公式的探讨》2011陈晓勇
[2]文献《数学及其应用》彭祖贈
[3]《MATLAB函数速查手册》李玉莉等编著
[4]《MATLAB仿真及电子信息应用》王亚芳
[5]姜启源.数学模型（第二版）.北京：高等教育出版社，1993.
[6]叶其孝.大学生数学建模竞赛辅导教材.长沙：湖南教育出版社，1993.
[7]阮沈勇，王永利.MATLAB程序设计.北京：电子工业出版社，2004.
[8]吕林根，许子道。解析几何.高等教育出版社，2006年5月.
附件
Ⅰ）直杆在太阳投影下影子的长度问题二：附件1中直杆的影长a=[1.1496258261.1821989761.2152969551.3179931491.3533640491.3893870911.5014816221.5402318171.5798533161.7032906331.746205911.7900509151.927918447];
1.2490510521.4261528561.6201445151.8350142721.283195341.4633998531.6612706131.880875001
问题三：附件2中直杆的影长a=[1.247256201.222794591.1989214861.152439571.129917471.107835481.086254201.0444462651.024264121.004640310.966790490.9485847350.9309278810.897109050.880973760.865492250.850504468];
1.175428964
1.0650810720.985490900.91375175
附件3中直杆的影长a=[3.5331421843.5467680293.5617976433.578100715
3.5957507833.614934283.6354259833.6572182723.6805411153.7051678363.7312780253.7589179113.7880878883.8187010153.8508096193.884585223.9199118283.9568759923.995534794.0357508354.077863059];
Ⅱ）程序
第一问：clcclearN=295;
%日期
b=2*pi*(N-1)/365;
delta=(0.006918-0.399912*cos(b)+0.070257*sin(b)-0.006758*cos(2*b)+
0.000907*sin(2*b)-0.002697*cos(3*b)+0.00148*sin(3*b));
T=-3:0.01:3;
%上午9点到下午3点
%赤纬角
fori=1:length(T)t=T(i)*15*pi/180;
%时间
%纬度
%高度角
weidu=39.943*pi/180;
h=abs(asin(sin(weidu)*sin(delta)+cos(delta)*cos(weidu)*cos(t)));
h=min(h,pi-h);H(i)=h*180/pi;
a=asin((cos(delta)*sin(t))/cos(h));A(i)=a*180/pi;end
plot(T+12,A,'r')%方位角变化曲线holdon
plot(T+12,H,'c')%高度角变化曲线L=3;%杆长
Lc=L*ones(1,length(T))./tan(H/180*pi);figure
plot(T+12,Lc);%杆长变化规律
%方位角
第二问：
程序一：
functionn=rqi(y,m,d)
if(mod(y,4)==0&&mod(y,100)~=0||mod(y,400)==0)z=1;else
z=0;end
m0=[0,31,z+59,z+90,z+120,z+151,z+181,z+212,z+243,z+273,z+304,z+334];%几月几号的天数n=m0(m)+d;程序二：function
[h]=hs(y,m,d,jd,wd,t0);
%判断闰年
n=rqi(y,m,d);fi=wd/180*pi;
gdj=23.45*sin((2*pi*(284+n))/365)/180*pi;t0=t0+(jd-120)/15;
t=15*(t0-12)/180*pi;
hh=sin(gdj).*sin(fi)+cos(gdj).*cos(fi).*cos(t);h=asin(hh);程序三：
Clearclc
a=[1.1496258261.1821989761.2152969551.2490510521.28319534
1.3179931491.3533640491.3893870911.4261528561.4633998531.5014816221.5402318171.5798533161.6201445151.6612706131.7032906331.746205911.7900509151.8350142721.8808750011.927918447];j0=0;
fori=0:1:180
forj=0:1:90
i0=0;j0=j0+1;
fort=14.7:0.05:15.7i0=i0+1;
b(i0)=cot(jiaodua);end
c(j0)=std(a./b);ifall(c(j0)
ii=i;jj=j;bb=a./b;
end
end
endjjiibb;
最度终的经纬度：
jj=19.3000
ii=108.6000
第三问：
附件2程序一：
functionn=rqi(y,m,d)
if(mod(y,4)==0&&mod(y,100)~=0||mod(y,400)==0)z=1;else
z=0;end
m0=[0,31,z+59,z+90,z+120,z+151,z+181,z+212,z+243,z+273,z+304,z+334];%几月几号的天数n=m0(m)+d;
%判断闰年
程序二：function
[h]=hs(y,m,d,jd,wd,t0);
n=rqi(y,m,d);fi=wd/180*pi;
gdj=23.45*sin((2*pi*(284+n))/365)/180*pi;t0=t0+(jd-120)/15;t=15*(t0-12)/180*pi;
hh=sin(gdj).*sin(fi)+cos(gdj).*cos(fi).*cos(t);h=asin(hh);
程序三：
clearclc
a=[1.247256201.222794591.1989214861.1754289641.152439571.129917471.107835481.086254201.0650810721.0444462651.024264121.004640310.985490900.966790490.9485847350.9309278810.913751750.897109050.880973760.865492250.850504468
];
j0=0;
for
k=70:1:85
fori=0:1:180
forend
end
endjjiikkbb;
%最度终的经纬度：
jj=
33
ii=
76
kk=
78
j=0:1:90
i0=0;j0=j0+1;
fort=12.7:0.05:13.7
jiaodua=drawyc(2015,2,k,i,j,t);i0=i0+1;
b(i0)=cot(jiaodua);end
c(j0)=std(a./b);ifall(c(j0)
ii=i;jj=j;kk=k;bb=a./b;
end
附件3
程序一：
functionn=rqi(y,m,d)
if(mod(y,4)==0&&mod(y,100)~=0||mod(y,400)==0)z=1;else
z=0;end
m0=[0,31,z+59,z+90,z+120,z+151,z+181,z+212,z+243,z+273,z+304,z+334];%几月几号的天数n=m0(m)+d;
%判断闰年
程序二：function
[h]=hs(y,m,d,jd,wd,t0);
n=rqi(y,m,d);fi=wd/180*pi;
gdj=23.45*sin((2*pi*(284+n))/365)/180*pi;t0=t0+(jd-120)/15;t=15*(t0-12)/180*pi;
hh=sin(gdj).*sin(fi)+cos(gdj).*cos(fi).*cos(t);h=asin(hh);
程序三：
clearclc
a=[3.5331421843.5467680293.5617976433.578100715
3.5957507833.614934283.6354259833.6572182723.6805411153.7051678363.7312780253.7589179113.7880878883.8187010153.8508096193.884585223.9199118283.9568759923.995534794.0357508354.077863059];j0=0;
fork=70:1:85
fori=65:1:180
for
end
end
endjjiikkbb;
%最终的经纬度：
jj=
46
ii=109
kk=
77
j=0:1:90i0=0;j0=j0+1;
fort=13.15:0.05:14.15
jiaodua=drawyc(2015,2,k,i,j,t);i0=i0+1;
b(i0)=cot(jiaodua);end
c(j0)=std(a./b);ifall(c(j0)
ii=i;jj=j;kk=k;bb=a./b;end
Ⅲ）模拟地球经纬度的程序
[jin,wei]=meshgrid(0:180,0:90);jin=jin*pi/180;wei=wei*pi/180;
mesh(cos(wei).*cos(jin),cos(wei).*sin(jin),sin(wei))