建立模型計算太陽影子長度和梯度下降法擬合太陽影子數(shù)據定位-全國大學生數(shù)學建模競賽A題論文

1、建立模型計算太陽影子長度和梯度下降法擬合太陽影子數(shù)據定位摘要本文建立了一個較為理想的幾何立體模型來計算某個日期，某個時間，某個緯度的太陽高度角，從而計算太陽照射直桿的影長。并且通過對一系列按時間變化記錄的數(shù)據以及局部參數(shù)，利用梯度下降法進行曲線擬合樣本數(shù)據點，以此求出相應的未知的參數(shù)，算出數(shù)據記錄時所在的可能的地點經緯度和日期。關鍵詞：幾何立體模型 擬合曲線 梯度下降法 正文建立太陽光照射地球，地球上垂直地面桿的影長模型 假設太陽光線完全平行，地球是完全標準的球體，地球公轉軌道是標準圓形，且隨日期變化勻速運轉。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用反三角函數(shù)即可知太陽高度角。1.2 建立太陽光線和地球赤道面所成角度隨日期變化的模型如圖，如果不考慮地球自轉，那么隨著公轉，太陽直射點將在地球上畫出一個圓，這個圓所在的平面即為回歸面。設A點為夏至日6月22日的直射點，此時AOB=N=2326，隨著日期推移，直射點移到點C，那么此時COD=正是我們所要求的太陽光線和地

4、球赤道面所成的角，而AOC那么是日期相對6月22日公轉導致移動所形成的夾角，所以AOC=相對6月22日天數(shù)差距/365 * (閏年為366)，其中比6月22日小那么天數(shù)差距為負值。易證，赤道面和回歸面的二面角大小為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用反三角函數(shù)即可求

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.根據某固定直桿在水平地面上的太陽影子頂點坐標數(shù)據，建立數(shù)學模型確定直桿所處的地點。將你們的模型應用于附件1的影子頂點坐標數(shù)據，給出假設干個可能的地點。附件1得到不同時刻的影長數(shù)據，我們采用梯度下降的方法，去擬合出一條曲線盡可能符合附件1的數(shù)據。根據模型，影長函數(shù)shadowlength L，=L * cot=L * cotarcsinsinsin + coscoscos可由日期2021年4月18日計算得到。每個數(shù)據樣本都可以計算出對應的。設數(shù)據樣本有n個，第i個樣本計算得到的記為i，第

7、i個樣本的影子長度記為Lexam i。首先定義一個代價函數(shù)：擬合即找到使代價函數(shù)到達盡可能小的值的參數(shù)L和。讓cost對L和分別求偏導數(shù)記為cost / L 和cost / 。具體求導由MATLAB完成，在代碼文件gradientdescent.m中實現(xiàn)初始化 L ，為 隨機值，進行屢次迭代每次迭代都修正參數(shù)L和的值，以減小cost的值。每次迭代修正方法：L := L - stepl * cost / L；:= - stepb * cost / ；那么cost將隨著L和的修正逐漸變小，即擬合程度變好。關于梯度下降法原理：詳見1。其中stepl 和 stepb是控制修正速度的值，在機器學習中也叫

8、學習速率，太大容易造成略過最優(yōu)點，太小造成修正緩慢，效率低下。此次計算中選用速率都為0.03，效果較好。迭代次數(shù)為40000次。并且為了防止出現(xiàn)局部最優(yōu)解，試用了屢次不同初始值。最終得到的擬合結果，每個樣本點影長差距遠小于1毫米。并且考慮了經度對時間的影響，試用了不同經度來擬合曲線。得到了擬合結果類似如下：由于差距很小，所以兩條線幾乎合并在一起，放大后如下列圖：最終將擬合后得到的一些可能的經緯度制成下列圖：如圖，在東經100125，北緯23.55都是可能的地點。桿長為2.6073米。根據某固定直桿在水平地面上的太陽影子頂點坐標數(shù)據，建立數(shù)學模型確定直桿所處的地點和日期。將你們的模型分別應用于附

9、件2和附件3的影子頂點坐標數(shù)據，給出假設干個可能的地點與日期。由于數(shù)據的進一步缺失，沒有交代日期，所以依舊采用梯度下降的方法擬合數(shù)據。但是多了一個變量需要擬合就是太陽光和赤道面的夾角。比照上一題多了一個需擬合的變量，同樣地求出cost關于的偏導數(shù)cost / 具體看MATLAB代碼文件gradientdescent.m。迭代過程中增加每次對的修正：:= - stepa * cost / 。由于多出了變量，導致擬合的可能性變多，所以對結果進行篩選，需要滿足cost0.008且擬合最終得到的參數(shù)|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