建立模型计算太阳影子长度和梯度下降法拟合太阳影子数据定位全国大学生数学建模竞赛a题.docx
《建立模型计算太阳影子长度和梯度下降法拟合太阳影子数据定位全国大学生数学建模竞赛a题.docx》由会员分享,可在线阅读,更多相关《建立模型计算太阳影子长度和梯度下降法拟合太阳影子数据定位全国大学生数学建模竞赛a题.docx(25页珍藏版)》请在冰豆网上搜索。
![建立模型计算太阳影子长度和梯度下降法拟合太阳影子数据定位全国大学生数学建模竞赛a题.docx](https://file1.bdocx.com/fileroot1/2022-12/30/4fa07263-e624-400e-8c8b-fd2f25d7197e/4fa07263-e624-400e-8c8b-fd2f25d7197e1.gif)
建立模型计算太阳影子长度和梯度下降法拟合太阳影子数据定位全国大学生数学建模竞赛a题
建立模型计算太阳影子长度和梯度下降法拟合太阳影子数据定位
1.摘要
本文建立了一个较为理想的几何立体模型来计算某个日期,某个时间,某个纬度的太阳高度角,从而计算太阳照射直杆的影长。
并且通过对一系列按时间变化记录的数据以及部分已知参数,利用梯度下降法进行曲线拟合样本数据点,以此求出相应的未知的参数,算出数据记录时所在的可能的地点经纬度和日期。
关键词:
几何立体模型拟合曲线梯度下降法
2.正文
1.建立太阳光照射地球,地球上垂直地面杆的影长模型
假设太阳光线完全平行,地球是完全标准的球体,地球公转轨道是标准圆形,且随日期变化匀速运转。
1.1建立太阳光线和地球上某一纬度β的上某点X切面所形成的太阳高度角模型
O为地球球心。
设太阳光线和地球的赤道面形成的角度为∠EOC=α,纬度即∠BOA=β,根据纬度的定义显然OB是β纬度上X点的切面(即X点地面)的法向量方向。
X点的经度和太阳垂直照射点所在的经度的夹角即∠AOC为Θ(可以由该点的真实时间计算得到,地球每小时自转15°,Θ=(该点真实时间-12)*15)。
欲求太阳在X点的高度角,相当于求X点地面法向量方向OB和太阳光线的夹角∠BOE的余角。
假设太阳高度角为γ,则sinγ=cos∠BOE。
如图作AC⊥OC,作EC和BA垂直于赤道面,BD⊥CE于D。
设OA长为z,则OB=z/cosβ,
OE=zcosΘ/cosα,
BD=AC=zsinΘ,
CD=AB=ztanβ,
CE=zcosΘtanα。
DE=CE-CD=zcosΘtanα-ztanβ。
BE²=BD²+DE²=(zsinΘ)²+(zcosΘtanα-ztanβ)²。
根据余弦定理可得,
经化简可得:
cos∠BOE=sinαsinβ+cosαcosβcosΘ。
即sinγ=sinαsinβ+cosαcosβcosΘ
用反三角函数即可知太阳高度角γ。
1.2建立太阳光线和地球赤道面所成角度α随日期变化的模型
如图,如果不考虑地球自转,则随着公转,太阳直射点将在地球上画出一个圆,这个圆所在的平面即为回归面。
设A点为夏至日(6月22日)的直射点,此时∠AOB=N=23°26′,随着日期推移,直射点移到点C,则此时∠COD=α正是我们所要求的太阳光线和地球赤道面所成的角,而∠AOC则是日期相对6月22日(公转导致)移动所形成的夹角,所以∠AOC=相对6月22日天数差距/365*π(闰年为366),其中比6月22日小则天数差距为负值。
易证,赤道面和回归面的二面角大小为N=23°26′,平面AOB⊥赤道面,平面COD⊥赤道面。
所以该问题可以转化为如下图所示的情形。
即一个二面角∠A’OB’大小为N,垂直于二面角∠A’OB’其中一个面B’OD’的面C’OD’,求C’OD’被这个二面角所夹出来的∠C’OD’(即∠α)大小。
图中△A’OB’,△A’OC’,△C’OD’为直角三角形。
C’D’=A’B’=OA’sinN,
OC’=OA’/cos∠A’OC’,
∴sinα=C’D’/OC’=sinN*cos∠A’OC’
用反三角函数即可求得α。
1.3建立某太阳高度角γ下,杆影长的模型
由于杆影长大部分时候远小于地球的半径,所以不考虑地球地面的弯曲,当做地面为平面。
则影长shadowlength=L*cotγ。
其中L为杆长。
1.4北京时间和某一经度longitude真实时间的关系。
由于北京时间是以东经120°为基准计算的,所以该经度longitude上的真实时间应该为realtime=(longitude-120)/15+beijingtime(小时),其中单位为小时和度,除以15因为地球每小时自转15°。
2.画出2015年10月22日北京时间9:
00-15:
00之间天安门广场(北纬39度54分26秒,东经116度23分29秒)3米高的直杆的太阳影子长度的变化曲线。
根据上述模型,在MATLAB上编码得到该变化曲线。
编码参照附录
3.根据某固定直杆在水平地面上的太阳影子顶点坐标数据,建立数学模型确定直杆所处的地点。
将你们的模型应用于附件1的影子顶点坐标数据,给出若干个可能的地点。
附件1得到不同时刻的影长数据,我们采用梯度下降的方法,去拟合出一条曲线尽可能符合附件1的数据。
根据模型,影长函数
shadowlength(L,α,β,Θ)=L*cotγ
=L*cot(arcsin(sinαsinβ+cosαcosβcosΘ))
α可由日期2015年4月18日计算得到。
每个数据样本都可以计算出对应的Θ。
设数据样本有n个,第i个样本计算得到的Θ记为Θi,第i个样本的影子长度记为Lexami。
首先定义一个代价函数:
拟合即找到使代价函数达到尽可能小的值的参数L和β。
让cost对L和β分别求偏导数记为∂cost/∂L和∂cost/∂β。
(具体求导由MATLAB完成,在代码文件gradientdescent.m中实现)
初始化L,β为随机值,进行多次迭代每次迭代都修正参数L和β的值,以减小cost的值。
每次迭代修正方法:
L:
=L-stepl*∂cost/∂L;
β:
=β-stepb*∂cost/∂β;
则cost将随着L和β的修正逐渐变小,即拟合程度变好。
关于梯度下降法原理:
详见[1]。
其中stepl和stepb是控制修正速度的值,在机器学习中也叫学习速率,太大容易造成略过最优点,太小造成修正缓慢,效率低下。
此次计算中选用速率都为0.03,效果较好。
迭代次数为40000次。
并且为了防止出现局部最优解,试用了多次不同初始值。
最终得到的拟合结果,每个样本点影长差距远小于1毫米。
并且考虑了经度对时间的影响,试用了不同经度来拟合曲线。
得到了拟合结果类似如下:
由于差距很小,所以两条线几乎合并在一起,放大后如下图:
最终将拟合后得到的一些可能的经纬度制成下图:
如图,在东经100°~125°,北纬23.55°都是可能的地点。
杆长为2.6073米。
3.根据某固定直杆在水平地面上的太阳影子顶点坐标数据,建立数学模型确定直杆所处的地点和日期。
将你们的模型分别应用于附件2和附件3的影子顶点坐标数据,给出若干个可能的地点与日期。
由于数据的进一步缺失,没有交代日期,所以依旧采用梯度下降的方法拟合数据。
但是多了一个变量需要拟合就是太阳光和赤道面的夹角α。
对比上一题多了一个需拟合的变量,同样地求出cost关于α的偏导数∂cost/∂α(具体看MATLAB代码文件gradientdescent.m)。
迭代过程中增加每次对α的修正:
α:
=α-stepa*∂cost/∂α。
由于多出了变量,导致拟合的可能性变多,所以对结果进行筛选,需要满足cost<0.008²且拟合最终得到的参数|α|<23°26′;|β|<66°34′,因为太阳光和赤道面所成角的范围所致。
考虑到附件2中影子在北京时间12点41分后直到13点41分仍慢慢减短,说明其只可能处于经度小于95度的范围内(这些地点的真实时间还没到12点,正在向12点靠近,所以影子在缩短)。
所以在东经50~95°范围每隔5°进行一次拟合。
拟合结果:
附件2可能的地点与日期以及杆长:
(7月1日,50°E,42.71°N,杆长0.846917m),
(7月13日,55°E,40.19°N,杆长0.90609m),
(7月24日,60°E,37.15°N,杆长0.980429m)
。
。
。
。
。
。
还有一些,具体看MATLAB代码的运行结果。
经度越靠近90度拟合效果越差。
同样地,附件3影子在13点09后影长数据不断增大,所以考虑拟合经度范围在东经110~120°范围每隔1°进行一次拟合(原本在105~150°每隔5°进行一次拟合,只有在110~120°之间有两组拟合成功数据,所以改为110~120°)。
拟合结果:
附件3可能的地点与日期以及杆长:
(9月18日,111°E,37.95°N,杆长3.059654m)
(9月7日,112°E,42.47°N,杆长3.027828m)
(8月27日,113°E,46.27°N,杆长3.005987m)
。
。
。
。
。
。
还有一些,具体看MATLAB代码的运行结果。
拟合较好的是在111°E到113°E之间。
4.附件4为一根直杆在太阳下的影子变化的视频,并且已通过某种方式估计出直杆的高度为2米。
请建立确定视频拍摄地点的数学模型,并应用你们的模型给出若干个可能的拍摄地点。
如果拍摄日期未知,你能否根据视频确定出拍摄地点与日期?
视频中给出了时间即相当于给出Θ,给出了杆长L,给出了日期相当于给出了α,所以唯一需要的就是纬度β。
用尺子粗略地测量了视频中影子的长度,得到以下数据:
杆长11厘米。
影长的变化情况:
8:
55、影长13厘米;8:
57、12.9厘米;8:
59、12.8厘米9:
01、12.6厘米;9:
03、12.4厘米;9:
05、12.2厘米;9:
07、12厘米;9:
09、11.8厘米;9:
11、11.7厘米;9:
13、11.6厘米;9:
15、11.5厘米;9:
17、11.3厘米;9:
19、11.2厘米;9:
21、11厘米;9:
23、10.8厘米;9:
25、10.7厘米;9:
27、10.55厘米;9:
29、10.4厘米;9:
31、10.3厘米;9:
33、10.2厘米
将数据输入程序,在经度为90°~125°对曲线拟合,得到β的情况,即可能的地点如下:
Possibleplace:
longitude:
90,latitude:
52.383329
Possibleplace:
longitude:
91,latitude:
50.940848
Possibleplace:
longitude:
92,latitude:
49.560709
Possibleplace:
longitude:
93,latitude:
48.241217
Possibleplace:
longitude:
94,latitude:
46.980504
Possibleplace:
longitude:
95,latitude:
45.776579
Possibleplace:
longitude:
96,latitude:
44.627373
Possibleplace:
longitude:
97,latitude:
43.530781
Possibleplace:
longitude:
98,latitude:
42.484696
Possibleplace:
longitude:
99,latitude:
41.487046
Possibleplace:
longitude:
100,latitude:
40.535827
Possibleplace:
longitude:
101,latitude:
39.629140
Possibleplace:
longitude:
102,latitude:
38.765238
Possibleplace:
longitude:
103,latitude:
37.942608
Possibleplace:
longitude:
104,latitude:
37.160098
Possibleplace:
longitude:
105,latitude:
36.417195
Possibleplace:
longitude:
106,latitude:
35.714660
Possibleplace:
longitude:
107,latitude:
35.056272
Possibleplace:
longitude:
108,latitude:
34.454905
Possibleplace:
longitude:
109,latitude:
33.964138
Possibleplace:
longitude:
110,latitude:
34.130393
Possibleplace:
longitude:
111,latitude:
40.322309
Possibleplace:
longitude:
112,latitude:
44.391538
Possibleplace:
longitude:
113,latitude:
47.150567
Possibleplace:
longitude:
114,latitude:
49.289467
Possibleplace:
longitude:
115,latitude:
51.049417
Possibleplace:
longitude:
116,latitude:
52.547724
Possibleplace:
longitude:
117,latitude:
53.851881
Possibleplace:
longitude:
118,latitude:
55.004924
Possibleplace:
longitude:
119,latitude:
56.036277
Possibleplace:
longitude:
120,latitude:
56.967101
Possibleplace:
longitude:
121,latitude:
57.813207
Possibleplace:
longitude:
122,latitude:
58.586770
Possibleplace:
longitude:
123,latitude:
59.297392
Possibleplace:
longitude:
124,latitude:
59.952799
Possibleplace:
longitude:
125,latitude:
60.559310
其中数据都为东经和北纬。
如果没有给出日期,则相当于对两个参数α和β进行梯度下降逼近。
依然可以计算出可能的地点(纬度β和经度决定)和日期(太阳和赤道面所成角α决定)。
参考文献:
[1]woxincd,梯度下降法,
网址:
访问时间:
2015年9月13日
附录:
使用软件:
MATLAB2012a
以下A1、A2、A3、A4为各个小题的解决程序脚本,其余为脚本调用的各个函数。
A1.m
clearall;clc;
%betaislatitude.
beta=((39*60+54)*60+26)/(180*60*60)*pi;
%alphaissunverticleangleonequatorialplane.
alpha=alphaondate(10,22,0);
%bjtime=[8:
0.05:
16];
bjtime=[9:
0.05:
15];
%thetaisthetimeangle.
theta=thetafuncontime(BeijingChange2realTime((116+(23+29/60)/60),bjtime));
%gammaissunverticleangleonearth.
gamma=asin(sinr(alpha,beta,theta));
shadowlength=cot(gamma)*3;
plot(bjtime,shadowlength);
xlabel('time');
ylabel('lengthofstick''sshadow');
A2.m
clearall;clc;
%alphaissunverticleangleonequatorialplane.
alpha=alphaondate(4,18,0);
%inputexample
x=[1.0365
1.0699
1.1038
1.1383
1.1732
1.2087
1.2448
1.2815
1.3189
1.3568
1.3955
1.4349
1.4751
1.516
1.5577
1.6003
1.6438
1.6882
1.7337
1.7801
1.8277
];
y=[0.4973
0.5029
0.5085
0.5142
0.5198
0.5255
0.5311
0.5368
0.5426
0.5483
0.5541
0.5598
0.5657
0.5715
0.5774
0.5833
0.5892
0.5952
0.6013
0.6074
0.6135
];
bjtime=[14+42/60:
3/60:
15+42/60]';
%size(bjtime)
shadowlexample=(x.^2+y.^2).^(1/2);
%size(lexample)
beta=60/180*pi;%40
l=2;%31.52
longtitude=100;%120116
i=1;
stepb=0.03;
stepl=0.03;%(0.03,0.1)(0.03,0.07)alternative
forlontitude=100:
5:
125
%initializationdone.
realtime=BeijingChange2realTime(longtitude,bjtime);
theta=thetafuncontime(realtime);
[betal]=Descent2(l,alpha,beta,theta,shadowlexample,stepb,stepl,40000);
latitude(i)=(beta/pi)*180;%recordthedifferentlatitudecombinedwithdifferentlontitude
i=i+1;
l%lshowthelengthofstick.
shadowlength=l*cot(asin(sinr(alpha,beta,theta)));
figure,plot(bjtime,shadowlength);%showthemodelpredictedlengthofshadow.
xlabel('beijingtime');
ylabel('shadowlength');
holdon;
plot(bjtime,shadowlexample,'r');
end
%shadowlexample%showthelengthmeasurementofshadow.
figure,plot([100:
5:
125],latitude);%showthepossibleplace.
xlabel('longtitude');
ylabel('latitude');
A3.m
clearall;clc;
x1=[-1.2352
-1.2081
-1.1813
-1.1546
-1.1281
-1.1018
-1.0756
-1.0496
-1.0237
-0.998
-0.9724
-0.947
-0.9217
-0.8965
-0.8714
-0.8464
-0.8215
-0.7967
-0.7719
-0.7473
-0.7227
];
y1=[0.173
0.189
0.2048
0.2203
0.2356
0.2505
0.2653
0.2798
0.294
0.308
0.3218
0.3354
0.3488
0.3619
0.3748
0.3876
0.4001
0.4124
0.4246
0.4366
0.4484
];
shadowlexample1=(x1.^2+y1.^2).^(1/2);
bjtime1=[12+41/60:
3/60:
13+41/60]';
x2=[1.1637
1.2212
1.2791
1.3373
1.396
1.4552
1.5148
1.575
1.6357
1.697
1.7589
1.8215
1.8848
1.9488
2.0136
2.0792
2.1457
2.2131
2.2815
2.3508
2.4213
];
y2=[3.336
3.3299
3.3242
3.3188
3.3137
3.3091
3.3048
3.3007
3.2971
3.2937
3.2907
3.2881
3.2859
3.284
3.2824
3.2813
3.2805
3.2801
3.2801
3.2804
3.2812
];
shadowlexample2=(x2.^2+y2.^2).^(1/2);
bjtime2=[13+09/60:
3/60:
14+09/60]';
alpha1=0/180*pi;
alpha2=0/180*pi;
beta1=60/180*pi;
beta2=45/180*pi;
l1=1;
l2=3;
longtitude=120;%longtitude2shouldbelargerthan(120-41/60*15)degree.
%Becausethedataisturnlarger,therealtimeontheplace..
%shouldbeafter12o'clock.
%whilelongitude1shouldbesmallerthan(120-15-41/60*15)
%forlongtitude=:
:
%initializationdone.
forlongtitude1=50:
5:
95
realtime1=BeijingChange2realTime(longtitude1,bjtime1);
theta1=thetafuncontime(realtime1);
stepa=0.003;
stepb=0.003;
stepl=0.003;
%%{
%[alpha1alpha2betal]=Descent3(l,alpha1,alpha2,beta,theta1,theta2,shadowlexample1,shadowlexample2,stepa,stepb,stepl,100000);
[alpha1beta1l1goodflag]=Descent3(l1,alpha1,beta1,theta1,shadowlexample1,stepa,stepb,s