建立模型计算太阳影子长度和梯度下降法拟合太阳影子数据定位-全国大学生数学建模竞赛A题论文

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

2、。设太阳光线和地球的赤道面形成的角度为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 = z / cos，OE=z cos / cos，BD=AC=z sin，CD=AB=z tan，CE=z cos t

3、an。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点为夏至日6月22日的直射点，此时AOB=N=2326，随着日期推移，直射点移到点C，那么此时COD=正是我们所要求的太阳光线和地

4、球赤道面所成的角，而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 = AB = OA sinN，OC=OA / cosAOC，sin=CD / OC = sinN * cosAOC用反三角函数即可求

5、得。1.3 建立某太阳高度角下，杆影长的模型由于杆影长大局部时候远小于地球的半径，所以不考虑地球地面的弯曲，当做地面为平面。那么影长shadowlength = L * cot 。其中L为杆长。1.4 北京时间和某一经度longitude真实时间的关系。由于北京时间是以东经120为基准计算的，所以该经度longitude上的真实时间应该为realtime =longitude - 120/ 15 + beijingtime (小时)，其中单位为小时和度，除以15因为地球每小时自转15。画出2021年10月22日北京时间9:00-15:00之间天安门广场北纬39度54分26秒,东经116度23分

6、29秒3米高的直杆的太阳影子长度的变化曲线。根据上述模型，在MATLAB上编码得到该变化曲线。编码参照附录3.根据某固定直杆在水平地面上的太阳影子顶点坐标数据，建立数学模型确定直杆所处的地点。将你们的模型应用于附件1的影子顶点坐标数据，给出假设干个可能的地点。附件1得到不同时刻的影长数据，我们采用梯度下降的方法，去拟合出一条曲线尽可能符合附件1的数据。根据模型，影长函数shadowlength L，=L * cot=L * cotarcsinsinsin + coscoscos可由日期2021年4月18日计算得到。每个数据样本都可以计算出对应的。设数据样本有n个，第i个样本计算得到的记为i，第

7、i个样本的影子长度记为Lexam i。首先定义一个代价函数：拟合即找到使代价函数到达尽可能小的值的参数L和。让cost对L和分别求偏导数记为cost / L 和cost / 。具体求导由MATLAB完成，在代码文件gradientdescent.m中实现初始化 L ，为 随机值，进行屡次迭代每次迭代都修正参数L和的值，以减小cost的值。每次迭代修正方法：L := L - stepl * cost / L；:= - stepb * cost / ；那么cost将随着L和的修正逐渐变小，即拟合程度变好。关于梯度下降法原理：详见1。其中stepl 和 stepb是控制修正速度的值，在机器学习中也叫

8、学习速率，太大容易造成略过最优点，太小造成修正缓慢，效率低下。此次计算中选用速率都为0.03，效果较好。迭代次数为40000次。并且为了防止出现局部最优解，试用了屡次不同初始值。最终得到的拟合结果，每个样本点影长差距远小于1毫米。并且考虑了经度对时间的影响，试用了不同经度来拟合曲线。得到了拟合结果类似如下：由于差距很小，所以两条线几乎合并在一起，放大后如下列图：最终将拟合后得到的一些可能的经纬度制成下列图：如图，在东经100125，北纬23.55都是可能的地点。杆长为2.6073米。根据某固定直杆在水平地面上的太阳影子顶点坐标数据，建立数学模型确定直杆所处的地点和日期。将你们的模型分别应用于附

9、件2和附件3的影子顶点坐标数据，给出假设干个可能的地点与日期。由于数据的进一步缺失，没有交代日期，所以依旧采用梯度下降的方法拟合数据。但是多了一个变量需要拟合就是太阳光和赤道面的夹角。比照上一题多了一个需拟合的变量，同样地求出cost关于的偏导数cost / 具体看MATLAB代码文件gradientdescent.m。迭代过程中增加每次对的修正：:= - stepa * cost / 。由于多出了变量，导致拟合的可能性变多，所以对结果进行筛选，需要满足cost0.008且拟合最终得到的参数|2326；| 23+26/60 | alpha1d 66 + 34/60 | beta1d 23+26

10、/60 | alpha2d 66 + 34/60 | beta2d M(month) day = day-M(month); month = month + 1;end month = month - 1; end Descent2.mfunction b l = Descent2( l,a,b,e,shadowlexample,stepb,stepl,itertimes )%UNTITLED5 Summary of this function goes here% Detailed explanation goes herefor i = 1:itertimes cost = costFun

11、c(l,a,b,e,shadowlexample); % to show cost in the descent process when debugging . %and deciding the parameters of gradient descent. dl da db = gradientdescent(l,a,b,e,shadowlexample); b = b - db*stepb; l = l - dl*stepl;end end Descent3.mfunction a b l goodflag = Descent3( l,a,b,e,shadowlexample,step

12、a,stepb,stepl,itertimes )%UNTITLED5 Summary of this function goes here% Detailed explanation goes herefor i = 1:itertimes cost = costFunc(l,a,b,e,shadowlexample); % to show cost in the descent process when debugging . %and deciding the parameters of gradient descent. dl da db = gradientdescent(l,a,b

13、,e,shadowlexample); b = b - db*stepb; l = l - dl*stepl; a = a - da*stepa; end if(cost 0.0082) goodflag = 1; else goodflag = 0; end end gradientdescent.mfunction gradl grada gradb = gradientdescent( l,a,b,e,lexample ) %length,alpha,beta,theta.theta is calculated by timeexample.%UNTITLED2 Summary of t

14、his function goes here% Detailed explanation goes hereexamplenumber = size(lexample,1);shadowlength = l*cot(asin(sinr(a,b,e); dl = (1 - (sin(a)*sin(b) + cos(a)*cos(b)*cos(e).2).(1/2). ./(sin(a)*sin(b) + cos(a)*cos(b)*cos(e); % dshadowlength/dl gradl = (1/examplenumber) * sum( 2*(shadowlength - lexam

15、ple) .* dl , 1); % dcost/dl da = - (l*(cos(a)*sin(b) - cos(b)*cos(e)*sin(a). ./(1 - (sin(a)*sin(b) + cos(a)*cos(b)*cos(e).2).(1/2). - (l*(cos(a)*sin(b) - cos(b)*cos(e)*sin(a).*(1 - (sin(a)*sin(b) + cos(a)*cos(b)*cos(e).2).(1/2). ./(sin(a)*sin(b) + cos(a)*cos(b)*cos(e).2; % dshadowlength/da grada = (

16、1/examplenumber) * sum( 2*(shadowlength - lexample) .* da , 1); % dcost/da db = - (l.*(cos(b)*sin(a) - cos(a)*cos(e)*sin(b). ./(1 - (sin(a)*sin(b) + cos(a)*cos(b)*cos(e).2).(1/2) . - (l*(cos(b)*sin(a) - cos(a)*cos(e)*sin(b).*(1 - (sin(a)*sin(b) + cos(a)*cos(b)*cos(e).2).(1/2). ./(sin(a)*sin(b) + cos

17、(a)*cos(b)*cos(e).2; % dshadowlength/db gradb = (1/examplenumber) * sum( 2*(shadowlength - lexample) .* db , 1); % dcost/db end sinr.mfunction sinr = sinr( a,b,e )%alpha,beta,theta%alpha sun vertical angle on equatorial planebeta latitude angletheta time angle% %UNTITLED Summary of this function goes here% Detailed explanation goes here%cosr = (cos(e)2/cos(a)2 + 1/cos(b)2 - (tan(b) - cos(e)*tan(a)2 - sin(e)2 ). /(2 * cos(e)/cos(a) * 1/cos(b)%cosr = (cos(e)2*cos(b)2 + cos(a)2 -(sin(b)*cos(a)-cos(e)*sin(a)*cos(b)2 - si