本科毕业设计建立模型计算太阳影子长度和梯度下降法拟合太阳影子数据定位全国大学生数学建模竞赛a题论文.docx

1、本科毕业设计建立模型计算太阳影子长度和梯度下降法拟合太阳影子数据定位全国大学生数学建模竞赛a题论文建立模型计算太阳影子长度和梯度下降法拟合太阳影子数据定位1摘要本文建立了一个较为理想的几何立体模型来计算某个日期，某个时间，某个纬度的太阳高度角，从而计算太阳照射直杆的影长。并且通过对一系列按时间变化记录的数据以及部分已知参数，利用梯度下降法进行曲线拟合样本数据点，以此求出相应的未知的参数，算出数据记录时所在的可能的地点经纬度和日期。关键词：几何立体模型 拟合曲线 梯度下降法 2正文1.建立太阳光照射地球，地球上垂直地面杆的影长模型 假设太阳光线完全平行，地球是完全标准的球体，地球公转轨道是标准圆

2、形，且随日期变化匀速运转。1.1 建立太阳光线和地球上某一纬度的上某点X切面所形成的太阳高度角模型O为地球球心。设太阳光线和地球的赤道面形成的角度为EOC = ，纬度即BOA=，根据纬度的定义显然OB是纬度上X点的切面（即X点地面）的法向量方向。X点的经度和太阳垂直照射点所在的经度的夹角即AOC为（可以由该点的真实时间计算得到，地球每小时自转15，=（该点真实时间-12） * 15 ）。欲求太阳在X点的高度角，相当于求X点地面法向量方向OB和太阳光线的夹角BOE的余角。假设太阳高度角为，则sin = cos BOE。如图作ACOC，作EC和BA垂直于赤道面，BDCE于D。设OA长为z，则OB

3、= z / cos，OE=z cos / cos，BD=AC=z sin，CD=AB=z tan，CE=z cos tan。DE = CE - CD =z cos tan - z tan。BE = BD + DE=（z sin）+（z cos tan - z tan）。根据余弦定理可得，经化简可得：cosBOE = sinsin + coscoscos。即sin = sinsin + coscoscos用反三角函数即可知太阳高度角。1.2 建立太阳光线和地球赤道面所成角度随日期变化的模型如图，如果不考虑地球自转，则随着公转，太阳直射点将在地球上画出一个圆，这个圆所在的平面即为回归面。设A点为夏

4、至日（6月22日）的直射点，此时AOB=N=2326，随着日期推移，直射点移到点C，则此时COD=正是我们所要求的太阳光线和地球赤道面所成的角，而AOC则是日期相对6月22日（公转导致）移动所形成的夹角，所以AOC=相对6月22日天数差距/365 * (闰年为366)，其中比6月22日小则天数差距为负值。易证，赤道面和回归面的二面角大小为N=2326，平面AOB赤道面，平面COD赤道面。所以该问题可以转化为如下图所示的情形。即一个二面角AOB大小为N，垂直于二面角AOB其中一个面BOD的面COD，求COD被这个二面角所夹出来的COD（即）大小。图中AOB，AOC，COD为直角三角形。CD =

5、AB = OA sinN，OC=OA / cosAOC，sin=CD / OC = sinN * cosAOC用反三角函数即可求得。1.3 建立某太阳高度角下，杆影长的模型由于杆影长大部分时候远小于地球的半径，所以不考虑地球地面的弯曲，当做地面为平面。则影长shadowlength = L * cot 。其中L为杆长。1.4 北京时间和某一经度longitude真实时间的关系。由于北京时间是以东经120为基准计算的，所以该经度longitude上的真实时间应该为realtime =（longitude - 120）/ 15 + beijingtime (小时)，其中单位为小时和度，除以15因为

6、地球每小时自转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（sinsi

7、n + coscoscos）可由日期2015年4月18日计算得到。每个数据样本都可以计算出对应的。设数据样本有n个，第i个样本计算得到的记为i，第i个样本的影子长度记为Lexam i。首先定义一个代价函数：拟合即找到使代价函数达到尽可能小的值的参数L和。让cost对L和分别求偏导数记为cost / L 和cost / 。（具体求导由MATLAB完成，在代码文件gradientdescent.m中实现）初始化 L ，为 随机值，进行多次迭代每次迭代都修正参数L和的值，以减小cost的值。每次迭代修正方法：L := L - stepl * cost / L；:= - stepb * cost /

8、；则cost将随着L和的修正逐渐变小，即拟合程度变好。关于梯度下降法原理：详见1。其中stepl 和 stepb是控制修正速度的值，在机器学习中也叫学习速率，太大容易造成略过最优点，太小造成修正缓慢，效率低下。此次计算中选用速率都为0.03，效果较好。迭代次数为40000次。并且为了防止出现局部最优解，试用了多次不同初始值。最终得到的拟合结果，每个样本点影长差距远小于1毫米。并且考虑了经度对时间的影响，试用了不同经度来拟合曲线。得到了拟合结果类似如下：由于差距很小，所以两条线几乎合并在一起，放大后如下图：最终将拟合后得到的一些可能的经纬度制成下图：如图，在东经100125，北纬23.55都是可

9、能的地点。杆长为2.6073米。3.根据某固定直杆在水平地面上的太阳影子顶点坐标数据，建立数学模型确定直杆所处的地点和日期。将你们的模型分别应用于附件2和附件3的影子顶点坐标数据，给出若干个可能的地点与日期。由于数据的进一步缺失，没有交代日期，所以依旧采用梯度下降的方法拟合数据。但是多了一个变量需要拟合就是太阳光和赤道面的夹角。对比上一题多了一个需拟合的变量，同样地求出cost关于的偏导数cost / （具体看MATLAB代码文件gradientdescent.m）。迭代过程中增加每次对的修正：:= - stepa * cost / 。由于多出了变量，导致拟合的可能性变多，所以对结果进行筛选，

10、需要满足cost0.008且拟合最终得到的参数|2326；|6634，因为太阳光和赤道面所成角的范围所致。考虑到附件2中影子在北京时间12点41分后直到13点41分仍慢慢减短，说明其只可能处于经度小于95度的范围内（这些地点的真实时间还没到12点，正在向12点靠近，所以影子在缩短）。所以在东经5095范围每隔5进行一次拟合。拟合结果：附件2 可能的地点与日期以及杆长：（7月1日，50E，42.71N，杆长0.846917m），（7月13日 ，55E，40.19N，杆长0.90609m），（7月24日，60E，37.15N，杆长0.980429m）。还有一些，具体看MATLAB代码的运行结果。经

11、度越靠近90度拟合效果越差。同样地，附件3影子在13点09后影长数据不断增大，所以考虑拟合经度范围在东经110120范围每隔1进行一次拟合（原本在105150每隔5进行一次拟合，只有在110120之间有两组拟合成功数据，所以改为110120）。拟合结果:附件3可能的地点与日期以及杆长：（9月18日，111E，37.95N，杆长3.059654m）（9月7日，112E，42.47N，杆长3.027828m）（8月27日，113E，46.27N，杆长3.005987m）。还有一些，具体看MATLAB代码的运行结果。拟合较好的是在111E到113E之间。4.附件4为一根直杆在太阳下的影子变化的视频，

12、并且已通过某种方式估计出直杆的高度为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.

13、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厘米将数据输入程序，在经度为90125对曲线拟合，得到的情况，即可能的地点如下：Possible place:longitude:90,latitude:52.383329Possible place:longitude:91,latitude:50.940848Possible place:longitude:92,latitude:49.560709Possible place

14、:longitude:93,latitude:48.241217Possible place:longitude:94,latitude:46.980504Possible place:longitude:95,latitude:45.776579Possible place:longitude:96,latitude:44.627373Possible place:longitude:97,latitude:43.530781Possible place:longitude:98,latitude:42.484696Possible place:longitude:99,latitude:4

15、1.487046Possible place:longitude:100,latitude:40.535827Possible place:longitude:101,latitude:39.629140Possible place:longitude:102,latitude:38.765238Possible place:longitude:103,latitude:37.942608Possible place:longitude:104,latitude:37.160098Possible place:longitude:105,latitude:36.417195Possible p

16、lace:longitude:106,latitude:35.714660Possible place:longitude:107,latitude:35.056272Possible place:longitude:108,latitude:34.454905Possible place:longitude:109,latitude:33.964138Possible place:longitude:110,latitude:34.130393Possible place:longitude:111,latitude:40.322309Possible place:longitude:112

17、,latitude:44.391538Possible place:longitude:113,latitude:47.150567Possible place:longitude:114,latitude:49.289467Possible place:longitude:115,latitude:51.049417Possible place:longitude:116,latitude:52.547724Possible place:longitude:117,latitude:53.851881Possible place:longitude:118,latitude:55.00492

18、4Possible place:longitude:119,latitude:56.036277Possible place:longitude:120,latitude:56.967101Possible place:longitude:121,latitude:57.813207Possible place:longitude:122,latitude:58.586770Possible place:longitude:123,latitude:59.297392Possible place:longitude:124,latitude:59.952799Possible place:lo

19、ngitude:125,latitude:60.559310其中数据都为东经和北纬。如果没有给出日期，则相当于对两个参数和进行梯度下降逼近。依然可以计算出可能的地点（纬度和经度决定）和日期（太阳和赤道面所成角决定）。参考文献：1woxincd，梯度下降法，网址：访问时间：2015年9月13日附录：使用软件：MATLAB 2012a以下A1、A2、A3、A4为各个小题的解决程序脚本，其余为脚本调用的各个函数。A1.mclear all;clc;%beta is latitude.beta = (39*60+54)*60+26)/(180*60*60) * pi; %alpha is sun ve

20、rticle angle on equatorial plane.alpha = alphaondate(10,22,0); %bjtime = 8:0.05:16;bjtime = 9:0.05:15;%theta is the time angle.theta = thetafuncontime(BeijingChange2realTime(116+(23+29/60)/60),bjtime); %gamma is sun verticle angle on earth.gamma = asin(sinr(alpha,beta,theta); shadowlength = cot(gamm

21、a) * 3; plot(bjtime,shadowlength);xlabel(time);ylabel(length of sticks shadow);A2.mclear all; clc; %alpha is sun verticle angle on equatorial plane.alpha = alphaondate(4,18,0); %input examplex = 1.03651.06991.10381.13831.17321.20871.24481.28151.31891.35681.39551.43491.47511.5161.55771.60031.64381.68

22、821.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; bjtime = 14+42/60:3/60:15+42/60;%size(bjtime)shadowlexample = (x.2 + y.2).(1/2);%size(lexample) beta = 60/180 * pi;%40l = 2;%3 1.5 2longtitude = 100

23、;%120 116i=1; stepb = 0.03;stepl = 0.03;%(0.03,0.1) (0.03,0.07)alternative for lontitude = 100:5:125%initialization done. realtime = BeijingChange2realTime(longtitude,bjtime); theta = thetafuncontime(realtime); beta l = Descent2(l,alpha,beta,theta,shadowlexample,stepb,stepl,40000); latitude(i) = (be

24、ta/pi) * 180; %record the different latitude combined with different lontitude i=i+1; l % l show the length of stick. shadowlength = l*cot(asin(sinr(alpha,beta,theta); figure,plot(bjtime,shadowlength);% show the model predicted length of shadow. xlabel(beijingtime); ylabel(shadowlength); hold on; pl

25、ot(bjtime,shadowlexample,r); end %shadowlexample % show the length measurement of shadow . figure,plot(100:5:125,latitude);% show the possible place. xlabel(longtitude); ylabel(latitude); A3.mclear all; 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.92

26、17-0.8965-0.8714-0.8464-0.8215-0.7967-0.7719-0.7473-0.7227; y1 = 0.1730.1890.20480.22030.23560.25050.26530.27980.2940.3080.32180.33540.34880.36190.37480.38760.40010.41240.42460.43660.4484; shadowlexample1 = (x1.2 + y1.2).(1/2); bjtime1 = 12+41/60:3/60:13+41/60; x2 = 1.16371.22121.27911.33731.3961.45

27、521.51481.5751.63571.6971.75891.82151.88481.94882.01362.07922.14572.21312.28152.35082.4213; y2 = 3.3363.32993.32423.31883.31373.30913.30483.30073.29713.29373.29073.28813.28593.2843.28243.28133.28053.28013.28013.28043.2812; shadowlexample2 = (x2.2 + y2.2).(1/2); bjtime2 = 13+09/60:3/60:14+09/60; alph

28、a1 = 0/180 * pi;alpha2 = 0/180 * pi; beta1 = 60/180 * pi;beta2 = 45/180 * pi; l1 = 1;l2 = 3; longtitude = 120; %longtitude2 should be larger than (120 - 41/60 * 15)degree. %Because the data is turn larger ,the realtime on the place. %should be after 12 oclock. %while longitude1 should be smaller tha

29、n (120 - 15 - 41/60*15) %for longtitude = : %initialization done. for longtitude1 = 50 : 5 : 95 realtime1 = BeijingChange2realTime(longtitude1,bjtime1); theta1 = thetafuncontime(realtime1); stepa = 0.003; stepb = 0.003; stepl = 0.003; % %alpha1 alpha2 beta l = Descent3(l,alpha1,alpha2,beta,theta1,theta2,shadowlexample1,shadowlexample2,stepa,stepb,stepl,100000 ); alpha1 beta1 l1 goodflag = Descent3(l1,alpha1,beta1,theta1,shadowlexample1,