A201525003029_楊俊龍_蘇婷_李存課件

1、太陽影子定位摘要本文對太陽影子定位問題進行了深入探討，建立了太陽影子和時間以及變位參數(shù)之間關系的數(shù)學模型，應用了matlab進行求解。針對問題一：已知時間與地點，求解影長與時間的關系。通過影長與太陽高度角、太陽方位角的關系，建立參數(shù)模型，運用matlab畫出影子長度變化曲線，影子長度與太陽高度角、太陽方位角、時間之間的關系如表1。表1 2015年10月22日9:00-15:00天安門廣場直桿（3米）的太陽影子長度當?shù)貢r間9:0010:0011:0012:0013:0014:0015:00太陽高度角23.319531.00136.155837.99136.155831.00123.32太陽方位角

2、A40.021139.61939.375939.29439.375939.619240.021影子高度L13.8417573.8413643.8411283.8410493.8411283.8413643.841757針對問題二，已知日期n與時間t以及影子的觀測坐標P，經(jīng)度、緯度、桿長h未知。分析太陽高度角、太陽方位角A與參數(shù)h,的關系，求解理論坐標Q。根據(jù)運籌學理論，以各PQ2之和最小為目標函數(shù)，以h,的取值范圍為約束條件，建立非線性規(guī)劃模型。應用最小二乘參數(shù)估計法思想，運用窮舉法和放縮法，利用matlab求解，得出,的值。針對問題三，已知t與P，在問題二的基礎上增加一個未知參數(shù)日期序號n，

3、分析、A與參數(shù)h,n的關系，求解理論坐標Q。以各PQ2之和最小為目標函數(shù)，以h,n的取值范圍為約束條件，建立非線性規(guī)劃模型。應用最小二乘參數(shù)估計法，利用matlab求解，得出,的值。針對問題四，利用視頻幀實現(xiàn)經(jīng)緯度估計的方法，根據(jù)未經(jīng)過校準的影子的位置，得出所拍攝圖像的經(jīng)緯度信息，利用matlab求解，得到相應的經(jīng)緯度。關鍵詞：影子定位，參數(shù)方程，非線性規(guī)劃，最小二乘參數(shù)估計法一、 問題重述太陽影子定位技術就是通過分析視頻中物體的太陽影子變化，確定視頻拍攝的地點和日期的一種方法。問題一：（1）建立影子長度變化的數(shù)學模型，分析影子長度關于各個參數(shù)的變化規(guī)律。（2）應用建立的模型畫出2015年10

4、月22日北京時間9:00-15:00之間天安門廣場（北緯39度54分26秒,東經(jīng)116度23分29秒）3米高的直桿的太陽影子長度的變化曲線。問題二：根據(jù)某固定直桿在水平地面上的太陽影子頂點坐標數(shù)據(jù)，建立數(shù)學模型確定直桿所處的地點。將你們的模型應用于附件1的影子頂點坐標數(shù)據(jù)，給出若干個可能的地點。問題三：根據(jù)某固定直桿在水平地面上的太陽影子頂點坐標數(shù)據(jù)，建立數(shù)學模型確定直桿所處的地點和日期。將你們的模型分別應用于附件2和附件3的影子頂點坐標數(shù)據(jù)，給出若干個可能的地點與日期。問題四：（1）附件4為一根直桿在太陽下的影子變化的視頻，并且已通過某種方式估計出直桿的高度為2米。請建立確定視頻拍攝地點的數(shù)

5、學模型，并應用你們的模型給出若干個可能的拍攝地點。（2）如果拍攝日期未知，如何根據(jù)視頻確定出拍攝地點與日期？二、 問題分析2.1問題一的分析（1）建立影子長度變化的數(shù)學模型，分析影子長度關于各個參數(shù)的變化規(guī)律。從太陽影子長度隨著時間的變化而改變的角度出發(fā)，利用不同時刻的太陽高度角及太陽方位角來計算太陽影子長度。建立太陽影子頂點坐標的參數(shù)模型。（2）應用建立的模型，代入數(shù)據(jù)，運用matlab進行擬合，得到2015年10月22日北京時間9:00-15:00之間天安門廣場3米高的直桿的太陽影子長度的變化曲線。2.2問題二的分析已知直桿頂點坐標投影和時間（即已知日期序號和時間t），直桿高度固定，建立數(shù)

6、學模型確定直桿所處的地點。要根據(jù)某固定直桿在水平地面上的太陽影子頂點坐標數(shù)據(jù)，建立數(shù)學模型確定直桿所處的地點。要將太陽高度角和太陽方位角分別進行討論，考慮到實際觀測值P所在的坐標系和正北方向存在夾角的話會使表達式過于復雜，于是假設實際觀測值P所在的坐標系和正北方向存在夾角為零。根據(jù)太陽影子長度與時間、直桿長度、經(jīng)緯度之間的關系，建立S關于參數(shù)h,的非線性規(guī)劃模型，運用最小二乘參數(shù)估計法求解，利用matlab求解找到最佳的h,值，從而確定附件1直桿所在地點。2.3問題三的分析已知直桿頂點坐標投影和時間（年月日時分秒，即已知日期序號和時間t），直桿高度固定，建立數(shù)學模型確定直桿所處的地點。在問題二

7、的基礎上根據(jù)太陽影子長度與時間、直桿長度、經(jīng)緯度、日期之間的關系，建立 S 關于參數(shù)h,n的最小二乘擬合模型, 運用最小二乘參數(shù)估計法求解，利用matlab求解找到最佳的h,n值，從而確定附件2、附件3中直桿所在地點以及日期。2.4問題四的分析（1）根據(jù)高度為2米直桿在太陽下的影子變化的視頻，確定視頻拍攝地點的數(shù)學模型。利用視頻幀實現(xiàn)經(jīng)緯度估計的方法，根據(jù)未經(jīng)過校準的影子的位置，得出所拍攝圖像的經(jīng)緯度信息5，建立相應的灰色模型，利用matlab求解，得到相應的經(jīng)緯度，得到拍攝地點。（2）若拍攝日期未知，根據(jù)視頻確定出拍攝地點與日期。三、問題假設（1）假設不考慮大氣折射率；（2）假設不考慮海拔、

8、大氣壓強、氣候等因素；（3）假設直桿是細桿，不考慮桿的粗細對影長的影響；（4）假設所給數(shù)據(jù)真實可靠。四、符號定義和說明h桿高L1影長t北京時間ts太陽時太陽高度角A太陽方位角當?shù)鼐暥犬數(shù)亟?jīng)度當時的太陽赤緯w時角n日期序號1n365五、模型的建立5.1問題一的模型建立參數(shù)模型的建立思路：太陽方位角A太陽高度角時差ET經(jīng)度赤緯角緯度時角w太陽時ts時間t日期n 5.1.1太陽高度角（1）計算時角w1 w=15（ts-12） （1）其中ts為太陽時2 ts=t-4×120-+ET （2）其中ET為時角 ET=229.2（0.000075+0.001868cos-0.032077sin-0.

9、014615cos2-0.04089sin2 （3） =360n-13651n365 （4）（2）計算赤緯角1 =23.45sin2(284+n)365 度 （5）（3）確定太陽高度角3太陽高度角是太陽相對于地平線的高度角，這是以太陽視盤面的幾何中心和理想地平線所夾的角度。太陽高度角可以使用下面的算式，經(jīng)由計算得到很好的近似值： sin=sinsin+coscoscosw （6）5.1.2太陽方位角 4太陽方位角是太陽在方位上的角度，它通常被定義為從北方沿著地平線順時針量度的角。它可以利用下面的公式，經(jīng)由計算得到良好的近似值，使用正弦函數(shù)： sinA=-sinwcoscos （7）使用余弦函數(shù)

10、：當時角為負值時 (上午)，方位角的角度小于180度，時角為正值時 (下午)，方位角應該大于180度，即要取補角的值。 cosA=sincos-coswcossincos （8） cosA=sin-sinsincoscos （9）5.1.3建立參數(shù)方程定義圖1中y軸方向為正北方向，影子頂點坐標為x，yx，yL1hz y x A圖1 “立竿見影”坐標系由影子長度與太陽高度角的關系 L1=htan （10）可得出參數(shù)方程： x=sinAL1y=cosAL1 t為參數(shù) （11）根據(jù)建立的參數(shù)方程，運用matlab擬合出影子長度與時間的變化曲線。5.2問題二的模型建立假設太陽影子頂點坐標實際觀測值為P

11、Pxi,Pyi,太陽影子頂點坐標理論值為Q(QxiQyi),假設P所在的坐標系和正北方向的夾角為0°，則有 =hsh, A=hah, Qxi=L1sinAQyi=L1cosA (12) L1=htan (13)由上述三個式子得到： Qxi=htansinA1Qyi=htancosA1 (14) 化解得到： Qxi=-sinwcoshsin (15) Qyi=sincos-coswcossinsinhQyi=sin-sinsinsincosh （16）注：當時角為負值時 (上午)，方位角的角度小于180度，時角為正值時 (下午)，方位角應該大于180度，即要取補角的值。通過建立理論坐標

12、QQxi，Qyi與觀測坐標PPxi，Pyi的距離的平方和S最小的非線性規(guī)模型進行求解。即： minS=i=121（Pxi-Qxi）2+（Qxi-Qyi）2 (17)s.t.0360°-90°90° 2<h3 5.3問題三的模型建立在第二問的基礎上增加一個未知參數(shù)n ,建立非線性規(guī)劃模型： min S=i=121Qxi-Pxi2+Qyi-Pyi2 （18）s.t.0°360°-90°90°2h31n3655.4問題四的模型建立 要根據(jù)一根高度為2米的直桿在太陽下的影子變化的視頻，建立確定視頻拍攝地點。就是要利用視頻幀實現(xiàn)

13、經(jīng)緯度估計的方法，能夠根據(jù)未經(jīng)過校準的影子的位置，得出所拍攝圖像的經(jīng)緯度信息5。具體步驟5如下：（1） 獲取視頻中含有至少兩個影子軌跡的視頻幀；（2） 對每一幀圖像，檢測影子軌跡；（3） 找出滅點，計算滅點；（4） 由滅點擬合出地平線；（5） 在圖像沿著有垂直關系的物體畫出兩條直線，根據(jù)這兩條直線分別與地平線的交點坐標，這兩個坐標就是兩個相互垂直的滅點vx,vy的坐標。（6） 根據(jù)公式f=v0vx3vy2+u0vx1vy3+v0vx2vy3+u0vx3vy1-vx1vy1-vx2vy2vx3vy3-u02-v02計算出焦距f,式中，u0,v0為圖像中心點的坐標，根據(jù)公式g=10-u001-v0

14、-u0-v0u02+v02+f2計算出g（7） 在地平線上找到另一對滿足vxTgvy=0限制的點vp,vq；（8） 實現(xiàn)從圖像到經(jīng)過度量糾正的世界坐標的轉換；（9） 根據(jù)X=Hx還原出經(jīng)過度量糾正的世界坐標中的點，其中X是經(jīng)過度量糾正的世界坐標中的點，x是圖像坐標中的點；（10） 擬合影子軌跡；（11） 利用日晷原理，計算緯度；（12） 運用極值法求出經(jīng)度。（13） 根據(jù)得到的經(jīng)緯度給出可能的視頻拍攝地點。六、模型的求解6.1問題一的模型的求解曲線擬合根據(jù)參數(shù)方程模型，代入數(shù)據(jù)日期序號n=295，當?shù)貣|經(jīng)=116.8667，北緯=39.9072，桿高h=3米。以y軸為正北方，求解直桿的頂點坐標

15、（x，y），得到2015年10月22日北京時間9:00-15:00之間天安門廣場3米高的直桿的太陽影子的運動軌跡如圖2，影子長度隨時間和方位角變化的曲線如圖3，影子長度隨方位角的變化曲線如圖4，影子長度隨時間的曲線如圖5，不同時刻對應的影長如表1。圖2 直桿影子的運動軌跡（附錄1）圖3影子長度隨時間和方位角變化的曲線（附錄3）圖4 影子長度隨方位角的變化曲線（附錄3）圖5 直桿影子隨時間的變化曲線圖（附錄2）表1 2015年10月22日9:00-15:00天安門廣場直桿（3米）的太陽影子長度（附錄4）當?shù)貢r間9:0010:0011:0012:0013:0014:0015:00太陽高度角23.3

16、19531.00136.155837.99136.155831.00123.32太陽方位角A40.021139.61939.375939.29439.375939.619240.021影子高度L13.8417573.8413643.8411283.8410493.8411283.8413643.841757觀察圖形，分析數(shù)據(jù)：太陽高度越小，影子越長。中午12:00太陽高度最大，影子最短。早晚太陽高度最小，中午最大，早上9：00到中午12：00太陽高度慢慢變大，影子變短，中午12：00到下午15：00太陽高度慢慢變小，影子變長。此數(shù)據(jù)與日常生活中實際現(xiàn)象相符合，可見此模型是合理的。6.2問題二的

17、模型求解最小二乘法用最小二乘參數(shù)估計法確定參數(shù)。最小二乘參數(shù)估計法基本思想：根據(jù)Sh,t的關系表達式求出太陽影子的理論坐標，通過計算理論坐標與觀測坐標的距離，建立距離平方和最小模型。通過對,h進行等間距的窮舉最終求得理論坐標QQxi，Qyi與觀測坐標PPxi，Pyi的距離的平方和S，當S取得最小值時，此時，h的值即為所求的最佳值。即求解如下最小二乘擬合模型。算法描述：輸入：時間t、日期序號n輸出：經(jīng)度、緯度、桿高h的值Step1: 建立以y軸為正北方向的空間直角坐標系，以地面為水平面，直桿垂直于地面，建立理論坐標QQxi，Qyi的參數(shù)方程，根據(jù)太陽高度角、太陽方位角與影長的關系，求出理論值QQ

18、xi，Qyi關于，h的關系式。Step2:時間t1，t2，t3···t21，日期序號n=108，等間距窮舉（附錄5）經(jīng)度0,360°，緯度-90°,90°，桿高h2，3，根據(jù)理論值QQxi，Qyi關于，h的關系式，求得理論坐標QQxi，Qyi。Step3:根據(jù)求出的理論坐標QQxi，Qyi與附錄1中的觀測坐標 PiPxi，Pyi （i0,21）,求出距離平方和S=i=121（Pxi-Qxi）2+（Qxi-Qyi）2令S0初始值為，S<S0,令S0=S，循環(huán)找出S0的最小值，最小值S=S0。當S取得最小值時，求得，h的值。算法結束

19、運用窮舉法與放縮法，用matlab求解得出最優(yōu)解，如表2。取值范圍0：30：360-90：30：902：0.1：3Step_1所得值120032.9145取值范圍90：5：120-30:5:302：0.1：:3Step_2所得值120-1032.6843取值范圍115:0.1:125-15:0.1：-52：0.1:3Step_3所得值118.8-9.130.1086取值范圍118.7:0.01:118.9-9.2:0.01：-92：0.1：:3Step_4所得值118.83-92.80.0007表2窮舉法的縮放過程（附錄6）得到，h可能值為=118.83°，=±9°

20、;，h=2.8m。6.3問題三的模型求解在問題二的基礎上，增加一個未知參數(shù)n，求解如下非線性規(guī)劃模型：minSh,n=i=121（Pxi-Qxi）2+（Qxi-Qyi）2在matlab中根據(jù)n的取值范圍1，365作為限制條件，引入一個for循環(huán)，同理運用窮舉法和放縮法進行求解。七、模型的評價7.1模型的優(yōu)點（1）通過建立非線性規(guī)劃模型,在利用未知參數(shù)取值范圍的條件下，運用matlab編程得出了影子頂點坐標理論值與實際值的距離平方之和，最小二乘法求解最小值，該方法結果合理可信。（2）運用matlab繪制出圖形，用數(shù)形結合的方法來進行分析，使模型思路更加清晰，更有說服力。 （3）建立參數(shù)

21、方程，可直接求解與逆向求解，模型具有普遍性。（4）通過對參數(shù)范圍的縮小，減小了大量的運算，同時提高了精度，使目標函數(shù)的最小值達到一個很小的水平，可以認為實際值和理論值幾乎相等S00.0007.7.2模型的缺點 （1）由于未考慮實際觀測值P所在的坐標系和正北方向存在的夾角，導致問題二的求解出的值存在一定的誤差。（2）由于運用窮舉法，未知參數(shù)過多，等間距窮舉參數(shù)有限，運算量過大，存在誤差。（3）模型假設過于理想化，和實際情況存在一定的差異。八、模型推廣本文所建的模型可以作為太陽影子定位方面的參考，具有借鑒意義；可以通過該地經(jīng)緯度、時間、日期、桿高擬合影長隨時間的變化曲線，可通過該地的影長

22、頂點坐標與時間求解所處位置等。可以應用在野外求生時地理位置的預測、某地的經(jīng)緯度預測等九、參考文獻1 方榮生，太陽能應用技術.北京：中國農(nóng)業(yè)機械出版社，1985.2 邱國全，太陽時的計算.科普資訊，1998. 3 Wikipedia，/wiki/Solar_zenith_angle，2015,9,144 Wikipedia，/wiki/Solar_azimuth_angle，2015,9,14.5操曉春，曲彥齡，孫濟洲，武琳，郭曉杰，張煒.基于視頻中太陽影子軌跡的經(jīng)緯度估計方法.天津：中華人民共和國國

23、家知識產(chǎn)權局，2009.6地球在線，7韓中庚等.數(shù)學建模使用教程.北京：高等教育出版社，2012.十、附錄附錄1%problem_1程序如下：clcclear allclose allbjtime = 9:1:15;jingdu = 116+23/60+29/3600n = 295;tao = 360*(n-1)/365;ET = 229.2*(0.000075+0.001868*cosd(tao)-0.032077*sind(tao)-0.014615*cosd(2*tao)-0.04089*sind(2*tao);ts = bjtime - 4*(120-jingdu)+ ET; fai

24、= 39+54/60+26/3600sigma = 23.45 .* sin(2.*pi*(284+n)/365)ts1 = 9:1/6:12;w = 15.*(ts1-12)alpha = asind(sind(fai).*sind(sigma)+cosd(fai).*cosd(sigma).*cosd(w)A1 =acosd(sind(sigma)-sind(alpha).*sind(fai)./(cosd(alpha).*cosd(fai) -180ts2 = 12:1/6:15;A2 =180- acosd(sind(sigma)-sind(alpha).*sind(fai)./(co

25、sd(alpha).*cosd(fai)hold onL1 = 3./tand(alpha)x1 = -sind(A2).*3./tand(alpha);x2 = sind(A2).*3./tand(alpha);y1 = cosd(A1).*3./tand(alpha);y2 = cosd(A2).*3./tand(alpha);axis(-6,6,-1,6)plot(x1,y1,'*')plot(x2,y2,'*')plot(0,0,'ro')xlabel('x')ylabel('y') 附錄2程序如下：clc

26、clear allclose allbjtime = 9:1/6:15;jingdu = 116+23/60+29/3600n = 295;tao = 360*(n-1)/365;ET = 229.2*(0.000075+0.001868*cosd(tao)-0.032077*sind(tao)-0.014615*cosd(2*tao)-0.04089*sind(2*tao);ts = bjtime - 4*(120-jingdu)+ ET; fai = 39+54/60+26/3600sigma = 23.45 .* sin(2.*pi*(284+n)/365)ts1 = 9:1/6:15;

27、 w = 15.*(ts1-12)alpha = asind(sind(fai).*sind(sigma)+cosd(fai).*cosd(sigma).*cosd(w)L=3./tand(alpha)plot(bjtime,L)xlabel('北京時間')ylabel('影子長度')grid on附錄3時間、太陽方位角與影子長度的關系圖：將三維圖旋轉一定角度后可得到時間、太陽方位角與影子長度的關系圖：太陽方位角與影子長度的關系圖：clcclear allclose allbjtime = 9:1:15;jingdu = 116+23/60+29/3600n =

28、 295;tao = 360*(n-1)/365;ET = 229.2*(0.000075+0.001868*cosd(tao)-0.032077*sind(tao)-0.014615*cosd(2*tao)-0.04089*sind(2*tao);ts = bjtime - 4*(120-jingdu)+ ET; fai = 39+54/60+26/3600sigma = 23.45 .* sin(2.*pi*(284+n)/365)ts = 9:1/6:15;w = 15.*(ts-12).*pi./180alpha = asind(sind(fai).*sind(sigma)+cosd(

29、fai).*cosd(sigma).*cosd(w) A1=acosd(sind(sigma)-sin(alpha).*sin(fai)./(cos(alpha).*cos(fai) hold ongrid onts = 9:1/6:15;w = 15.*(ts-12)alpha = asind(sind(fai).*sind(sigma)+cosd(fai).*cosd(sigma).*cosd(w)L1 = 3./tand(alpha)plot3(ts,A1,L1,'r*-')xlabel('時刻')ylabel('方位角/度')zlabel

30、('影子高度/米')附錄4時間、太陽高度角、太陽方位角以及相應影子長度的數(shù)據(jù)如下表：時間9:009:109:209:309:409:5010:0010:1010:2010:30高度角23.319524.739126.10827.422428.678629.872931.00132.05933.04333.948A140.021139.943339.869839.800639.735839.675339.61939.56839.5239.477L16.959396.510756.121615.782065.484465.222884.99264.794.6124.4563時間10

31、:4010:5011:0011:1011:2011:3011:4011:5012:0012:10高度角34.771635.508636.155836.709737.167337.526137.78437.93937.99137.939A139.439239.405339.375939.351139.330739.314939.30439.29739.29439.297L14.3214.204514.105634.023393.957043.905993.86983.84823.84113.8482時間12:2012:3012:4012:5013:0013:1013:2013:3013:4013

32、:50高度角37.783937.526137.167336.709736.155835.508634.77233.94833.04332.059A139.303639.314939.330739.351139.375939.405339.43939.47739.5239.568L13.869823.905993.957044.023394.105634.204514.3214.45634.6124.79時間14:0014:1014:2014:3014:4014:5015:00高度角31.00129.872928.678627.422426.10824.739123.32A139.619239.

33、675339.735839.800639.869839.943340.021L14.992645.222885.484465.782066.121616.510756.9594附錄5%Problem_2clcclear allclose allbjtime = 14.70:0.05:15.70; n = 108;sigma = 23.45.*sind(2.*pi.*(284+n)./365);tao = 360*(n-1)/365;a=0;b=0;c=0;d=0;ET = 229.2*(0.000075+0.001868*cosd(tao)-0.032077*sind(tao)-0.01461

34、5*cosd(2*tao)-0.04089*sind(2*tao);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

35、0.5892 0.5952 0.6013 0.6074 0.6135;a=0;min=inf; %step_1for jingdu = 0:30:360 ts = bjtime - 4.*(120-jingdu)+ET; w = 15./(ts-12); for fai= -90:30:90 alpha = asind(sind(fai).*sind(sigma)+cosd(fai).*cosd(sigma).*cosd(w); A = acosd(sind(sigma).*cosd(fai)-cosd(w).*cosd(sigma).*sind(fai)./cosd(alpha); for

36、h = 2:0.1:3 x = sind(A).*h./tand(alpha); y = cosd(A).*h./tand(alpha); c=c+1; for i=1:21 D = (x(i)-X(i)2+(y(i)-Y(i)2; end s = sum(sum(D(1:1,:); if s<min min = s ; d = jingdu fai h; end end end endmin,d%step_2for jingdu = 90:5:150 ts = bjtime - 4.*(120-jingdu)+ET; w = 15./(ts-12); for fai= -30:5:30

37、 alpha = asind(sind(fai).*sind(sigma)+cosd(fai).*cosd(sigma).*cosd(w); A = acosd(sind(sigma).*cosd(fai)-cosd(w).*cosd(sigma).*sind(fai)./cosd(alpha); for h = 2:0.1:3 x = sind(A).*h./tand(alpha); y = cosd(A).*h./tand(alpha); c=c+1; for i=1:21 D = (x(i)-X(i)2+(y(i)-Y(i)2; end s = sum(sum(D(1:1,:); if s<