)其中，表示北京時間矩陣，表示北京時間下測量的影子長度矩陣。根據(jù)附表2中數(shù)據(jù)通過計算得到表3：表3：計算后附表2數(shù)據(jù)北京時間坐標坐標影子長度12:41-1.23520.1731.24725620512:44-1.20810.1891.2227945912:47-1.18130.20481.19892148612:50-1.15460.2203153-1.12810.2356156-1.10180.25051.1299174712:59-1.07560.26531.1078354813:02-1.04960.27981.08625420613:05-1.02370.2941.06508107213:08-0.9980.3081.04444626513:11-0.97240.32181.02426412613:14-0.9470.33541.00464031413:17-0.92170.34880.98549090813:20-0.89650.36190.96679049413:23-0.87140.37480.94858473513:26-0.84640.38760.93092788113:29-0.82150.40010.9137517513:32-0.79670.41240.89710905113:35-0.77190.42460.88097376213:38-0.74730.43660.86549225913:41-0.72270.44840.850504468根據(jù)附表3中數(shù)據(jù)通過計算得到表3：表4：計算后附表3數(shù)據(jù)北京時間坐標坐標影子長度13:091.16373.3363.53314218413:121.22123.32993.54676802913:151.27913.32423.56179764313:181.33733.31883.57810071513:211.3963.31373.59575078313:241.45523.30913.6149342813:271.51483.30483.63542598313:301.5753.30073.65721827213:331.63573.29713.68054111513:361.6973.29373.70516783613:391.75893.29073.73127802513:421.82153.28813.75891791113:451.88483.28593.78808788813:481.94883.2843.81870101513:512.01363.28243.85080961913:542.07923.28133.8845852213:572.14573.28053.91991182814:002.21313.28013.95687599214:032.28153.28013.9955347914:062.35083.28044.03575083514:092.42133.28124.077863059由于多了一個未知量間隔時間，采用模型三的解法無法高效進行計算，遂采用牛頓迭代法進行求解：Step1：開始，編寫自定義函數(shù)mulNewton.m。Step2：準備，選定初始近似值，計算，。其中，是影子長度的關(guān)系式。Step3：迭代，按照公式迭代一次，得到新的近似值，計算，。Step4：控制，如果滿足。則終止迭代，以作為所求的根；否則轉(zhuǎn)Step5。此時是允許誤差，而 Step5：修改，如果迭代次數(shù)達到預(yù)定指定的次數(shù)10000，或者則方法失??；否則以代替轉(zhuǎn)Step3繼續(xù)迭代。按照上面計算步驟，通過MATLAB計算得到最優(yōu)解：表5：附件2和附件3最優(yōu)解時間經(jīng)緯度直桿高度附件25月20日新疆2.01米附件39月21日內(nèi)蒙古4.16米5.4問題四求解附錄4是一段視頻，我們能從視頻中看出直桿影子長度在隨時間變化，視頻中已經(jīng)給出拍攝日期。我們首先用MATLAB使用函數(shù)“VideoReader”讀取視頻，由于附件1中時間間隔是3分鐘，所以我們按照3分鐘提取一幀圖像，使用“imwrite“保存圖片，一共得到15張圖片。圖11MATLAB處理視頻得到圖片由于特征不明顯，程序無法直接對影子進行識別，所以，我們通過Photoshop對所有圖片的對比度進行改善，圖12Photoshop處理后的圖片然后對圖片進行灰度處理和圖像二值化：圖13MATLAB圖像處理后的圖片再使用函數(shù)測量圖片上直桿高度和影子長度，接著通過直桿實際高度和圖片上直桿高度和影子長度的比例算出實際影子長度，見表6表6：視頻處理后數(shù)據(jù)北京時間影子長度（像素）影子長度（米）8:54726.20372.12038:57717.18762.0949:00702.12012.059:03686.10512.00329:06675.08931.97119:09658.07451.92149:12643.04181.87759:156311.84239:186221.81619:216041.76359:245921.72859:275781.68769:305641.64679:335501.60589:345441.5883根據(jù)問題二的模型三，代入表中數(shù)據(jù)，利用MATLAB計算得到最優(yōu)解：呼和浩特。5.5問題五求解拍攝時間未知，時間間隔作為未知量，根據(jù)問題三的模型四，代入表中數(shù)據(jù)，利用MATLAB計算得到最優(yōu)解：呼和浩特,日期為6月18日。六、模型評價與推廣1.優(yōu)點：1、在影子長度計算模型中，利用黃經(jīng)推導和日期推導太陽赤緯，按季節(jié)分段時進行計算，確保推算的相對精確。2、在一元方程求解模型中，擬合圖像得出經(jīng)度，計算簡單。3、在多元非線性超定方程組模型中，能夠計算多組解。4、模型具有較強的適用性。2.缺點：1、在影子長度計算模型中，按季節(jié)分段時進行計算太陽赤緯，公式可能因為時間間隔不同而不同，增加計算復雜度。2、在一元方程求解模型中，擬合圖像得出經(jīng)度，由于擬合是二次函數(shù)，所以存在一定誤差，精確度不高。3、在多元非線性超定方程組模型中，求解算法復雜，計算時間較長。3.推廣：1、實際應(yīng)用中，在進行影子長度計算的數(shù)學模型中，考慮海拔、地心引力、折射率、壓強、溫度等因素能得到更加精確的地理位置。2、可以按照算法做成一個手機軟件，通過手機測量執(zhí)機人的影子長度和方位等信息反饋給終端設(shè)備，來確定他所在的位置，可以增加手機的防盜功能，也有助于尋找走失的老人和小孩。3、可以根據(jù)模型做成手機導航軟件，直接通過影子導航，代替手機地圖導航。七、參考文獻[1]費云翔，對太陽高度角的了解及其計算方法，中小企業(yè)管理與科技旬刊，116卷117期，2008年[2]林根石，利用太陽視坐標的計算進行物高測量與定位，南京林業(yè)大學學報，第15卷第3期，1991年9月[3]蔣洪力，太陽直射點緯度的數(shù)學推導和分析，數(shù)學通報，第46卷第9期，2007年[4]賀曉雷，太陽方位角的公式求解及其應(yīng)用，太陽能學報，第29卷第1期，2008年1月[5]劉浩

韓晶，MATLAB一點通，北京，飛電子工業(yè)出版社，2013年1月附錄源程序：第一問：1.作出天安門廣場影子隨時間變化：draw2.m%h太陽高度角%b當?shù)鼐暥?a太陽直射緯度%c太陽時角n=215;%北京時間t=9:0.05:15;%求太陽直射緯度a=-asin(0.39775*sin((n-186)*pi/180));%當?shù)氐木暥萣=(((26/60)+54)/60+39)*pi/180;%太陽時角c=(15*(t+0.05-12))*pi/180;%太陽高度角h=asin(sin(b)*sin(a)+cos(b)*cos(a)*cos(c));%影子長度s=3./tan(h);plot(t,s)xlabel('時間');ylabel('影子的長度');title('天安門廣場影子隨時間變化');2．畫出影子長度隨著緯度，經(jīng)度，間隔天數(shù)，時間的變化圖像。draw1.m%畫出l=cot(arcsin(cos(W)));影子長度隨緯度的變化W=3:0.1:53;l=cot(asin(cos(W*pi/180)));plot(W,l)title('影子長度隨緯度的變化');xlabel('緯度');ylabel('影子的長度');%y=x;x3=[1234567];y3=[1234567];plot(x3,y3,'b')title('影子長度隨直竿高度的變化');xlabel('直桿高度');ylabel('影子的長度');%l=cot(pi/3*(120-N));影子長度隨經(jīng)度的變化N=73*pi/180:0.1:135*pi/180;l=cot(pi/3.*(120-N));plot(N,l)title('影子長度隨經(jīng)度的變化');xlabel('經(jīng)度');ylabel('影子的長度');%%%%l=cot(asin(cos(asin(0.39775*sin(pi/186*n1)))));影子長度隨間隔天數(shù)的變化n1=0:1:180;l=cot(asin(cos(asin(0.39775*sin(pi/186*n1)))));plot(n1,l)title('影子長度隨間隔天數(shù)的變化');xlabel('間隔天數(shù)');ylabel('影子的長度');%%l=cot(asin(cos(asin(0.39775*sin(pi/186*n1)))));影子長度隨時間的變化t1=7:0.1:17;l=cot(asin(cos(pi/12*(t1-12))));plot(t1,l)title('影子長度隨時間的變化');xlabel('北京時間');ylabel('影子的長度');%%%第二問：求解緯度，經(jīng)度Jianmo2_1.m%x1(1)高度;x1(2)緯度;x1(3)經(jīng)度x=[1.03651.06991.10381.13831.17321.20871.24481.28151.31891.35681.39551.43491.47511.51601.55771.60031.64381.68821.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];fun=inline('sqrt(x1(1).^2/(tan(asin(sin(x1(2)*pi/180)*0.3336+0.7168*cos(x1(2)*pi/180)*cos(15*(15.2-(120-x1(3))*0.067)-12))))^2-x.^2)','x1','x');r1=1.91;r2=20.48;r3=111.08;x1=nlinfit(x,y,fun,[r1r2r3])第三問：求解經(jīng)度，緯度Jianmo3.m%H高度x2(1)%W緯度x2（2）%N經(jīng)度x2(3);%n天數(shù)x2(4);fun=inline('(sqrt((x2(1).*x2(1)/(tan(asin(sin(x2(2)*pi/180)*sin(sin((2*pi*(284+x2(4)))/365))))+cos(x2(2)*pi/180)*cos(sin((2*pi*(284+x2(4)))/365))*cos(15*((pi/180*15.2-(120-x(3))*4)-12)))^2)-x.^2))','x2','x');x2=nlinfit(x,y,fun,[15.801360.227856.4981-13])%牛頓迭代法function[x,m]=mulBroynewton(F,x0,A,eps)ifnargin==2A=eye(length(x0));eps=1.0e-4;elseifnargin==3eps=1.0e-1;endendx0=transpose(x0);Fx=subs(F,findsym(F),x0);x=x0-A\Fx;m=1;tol=1;whiletol>epsx0=x;Fx=subs(F,findsym(F),x0);x=x0-A\Fx;y=x-x0;Fr=subs(F,findsym(F),x0);x=x0-A\Fx;y=x-x0;Fr=subs(F,findsym(F),x);z=Fr-Fx;A1=A+(z-A*y)*transpose(y)/norm(y);A=A1;m=m+1;if(m>100000)disp('迭代次數(shù)太多，可能不收斂');return;endtol=norm(x-x0);end第四問：求解視頻上的測量位置：jianmo4.m%W緯度x(1)%N經(jīng)度x(2)%圖片上的影子長度單位（像素）length=[726.2037,717.1876,702.1201,686.1051,675.0893,...658.0745,643.0418,631.0000,622.0000,604.0000...592.0000,578.0000,564.0000,550.0000,544.0000];%桿長單位（像素）l=685;%影子的實際長度單位（米）length_yz=(2/l).*length;%時間time4=[8.90008.95009.00009.05009.10009.15009.20009.25009.30009.35009.40009.45009.50009.55009.5667];%fun=inline('0.2492*sin(W*pi/180)+0.9685*cos(W*pi/180)*cos(15*pi/180*(t+(120-N)*4-12))-H/sqrt(length_yingzi.^2+H.^2)','x1','x');fun=inline('sqrt(4/(0.2492*sin(x1(1)*pi/180)+0.9685*cos(x1(1)*pi/180)*cos(15*pi/180*(time4(1:,)+(120-x1(2))*4-12))).^2-4)','x1','time4');x1=nlinfit(time4,length_yz,fun,[12])jianmo5.m%W緯度x(1)%N經(jīng)度x(2)%圖片上的影子長度單位（像素）lengt