A201528009009_黃悅悅_楊萬里_董對林課件

1、太陽影子定位摘要：通過太陽影子來確定時間是人類歷史上常用的方法，本文主要討論了不同日期，不同時間，不同地理位置上固定直桿影子的變化規(guī)律，涉及太陽高度角，地球時角，經(jīng)緯度等天文知識。針對問題一，固定日期時間和經(jīng)緯度地點(diǎn)3米高直桿影子變化規(guī)律，首先根據(jù)天文知識、幾何知識和題意建立了以東方向?yàn)閤軸，北方向?yàn)閥軸，直桿為z軸的坐標(biāo)系，推導(dǎo)出根據(jù)日期時間經(jīng)緯度來計(jì)算影子長度公式，影子端點(diǎn)坐標(biāo)隨時間變換模型公式（）.然后根據(jù)模型，運(yùn)用MATLAB軟件計(jì)算出該問題的數(shù)據(jù)，并繪出曲線圖。針對問題二，給出影子測量坐標(biāo)要求根據(jù)模型來計(jì)算直桿所在經(jīng)緯度，雖然是問題一的逆問題，但涉及公式不完全可逆，首先利用測量影子坐

2、標(biāo)，用二次曲線擬合影子曲線，運(yùn)用擬合函數(shù)來計(jì)算地方時正午12點(diǎn)的影子坐標(biāo)，利用該數(shù)據(jù)，簡化了太陽高度計(jì)算公式，容易得到緯度計(jì)算式。然后擬合影子長度與時間的函數(shù)關(guān)系，根據(jù)影子最短時是地方時12點(diǎn)，算出該地方時與北京時間的時差，從而計(jì)算出與北京的經(jīng)度差，最終導(dǎo)出經(jīng)度。針對問題三，拍攝的日期，經(jīng)緯度都是未知，只能根據(jù)影子坐標(biāo)數(shù)據(jù)來計(jì)算，但計(jì)算原理與問題二的相似，先來計(jì)算出該問題的經(jīng)度，再由于日期為1-365，提出了窮舉計(jì)算模型（12），來計(jì)算緯度。針對問題四，要根據(jù)視頻中直桿影子變化來確定拍攝的地點(diǎn)，時間，日期，還有攝像機(jī)的位置，這問題首先要從視頻抓取圖片，在圖片中提取影子坐標(biāo)數(shù)據(jù)，要從平面數(shù)據(jù)轉(zhuǎn)換

3、到三維數(shù)據(jù)，隸屬于射影變化，因此建立了射影變換模型，并設(shè)計(jì)了計(jì)算的流程圖。關(guān)鍵詞：太陽高度； 太陽赤緯角；直桿影長；經(jīng)緯度；射影變換；1、 問題重述如何確定視頻的拍攝地點(diǎn)和拍攝日期是視頻數(shù)據(jù)分析的重要方面，太陽影子定位技術(shù)就是通過分析視頻中物體的太陽影子變化，確定視頻拍攝的地點(diǎn)和日期的一種方法。問題一： 建立模型，分析影子長度關(guān)于各個參數(shù)的變化規(guī)律，并用自己建立的模型來求解2015年10月22日北京時間9:00-15:00之間天安門廣場（北緯39度54分26秒,東經(jīng)116度23分29秒）3米高的直桿的太陽影子長度的變化曲線。問題二： 根據(jù)某固定直桿在水平地面上旗桿影子頂端坐標(biāo)數(shù)據(jù)，建立模型確定

4、直桿所處的地點(diǎn)。并用你們的模型和附件1的數(shù)據(jù)來求出該直桿的可能地點(diǎn)。問題三： 根據(jù)某固定直桿在水平地面上的太陽影子頂點(diǎn)坐標(biāo)數(shù)據(jù)，建立數(shù)學(xué)模型確定直桿所處的地點(diǎn)和日期。將你們的模型分別應(yīng)用于附件2和附件3的影子頂點(diǎn)坐標(biāo)數(shù)據(jù)，給出若干個可能的地點(diǎn)與日期。問題四：附件4為一根直桿在太陽下的影子變化的視頻，并且已估計(jì)出直桿的高度為2米。要建立確定視頻拍攝地點(diǎn)的數(shù)學(xué)模型，并應(yīng)用所建模型給出若干個可能的拍攝地點(diǎn)。進(jìn)一步如果拍攝日期未知，你能否根據(jù)視頻確定出拍攝地點(diǎn)與日期？二、問題分析 2.1對問題一的分析 針對問題一，運(yùn)用地球與太陽之間的公轉(zhuǎn)和自轉(zhuǎn)知識，進(jìn)行逆向運(yùn)用找出太陽光投射到地表上直桿的影子變化軌跡

5、。然后運(yùn)用MATLAB軟件進(jìn)行分析，再結(jié)合天文數(shù)學(xué)方法，進(jìn)而確定影子長度的變化規(guī)律。最后運(yùn)用影子長度變化的數(shù)學(xué)模型畫出2015年10月22日北京時間9:00-15:00之間天安門廣場（北緯39度54分26秒，東經(jīng)116度23分29秒）3米高的直桿的太陽影子長度的變化曲線。地球上日照的變化, 是由地球自轉(zhuǎn)和繞太陽公轉(zhuǎn)引起的。一天中, 地球自轉(zhuǎn)一周360°,(地球每小時自轉(zhuǎn)15°, 稱之為時角(R), 太陽位置也隨時間變化, 因此可獲得不同時間竿頂落影點(diǎn)P, 諸P點(diǎn)形成了一天的竿影軌跡線。1桿影頂端的軌跡與當(dāng)?shù)亟?jīng)緯度、太陽高度角以及時角（C）存在著函數(shù)關(guān)系，所以只要利用幾何知識

6、運(yùn)算出這幾者之間的關(guān)系，就解決了該問題。2.2問題二的分析：從問題一中得知，我們可以通過時角、經(jīng)緯度、太陽高度就可以描繪出旗桿影子長度的變化曲線，當(dāng)影子在平面xy平面上運(yùn)動時，若以東西和南北分別以y軸和x軸、桿為z軸建立空間直角坐標(biāo)系，就可以描繪影子頂端xy平面上的運(yùn)動關(guān)系，相反，當(dāng)然可以通過變化的坐標(biāo)點(diǎn)以及一些參變量就可以求出直桿所在的地點(diǎn)。2.3問題三的分析：本問題和問題二屬于同一類問題，其實(shí)還是問題一的反求過程，所以先準(zhǔn)備用附件2和附件3所給的數(shù)據(jù)進(jìn)行最小二乘法進(jìn)行擬合來計(jì)算出直桿的頂端在地上形成影子的端點(diǎn)劃過的軌跡。再通過擬合得到的方程來計(jì)算直桿所在的地點(diǎn)與日期。2.4問題四的分析：本

7、問題和前三個不同之處在于第一個問題是利用幾何關(guān)系來計(jì)算直桿太陽影子長度變化的曲線，二、三問是已知變化的曲線，來求日期和高桿所在地的具體經(jīng)緯度，而本題通過一個40分鐘的視頻，通過閱讀視頻的內(nèi)容。已知在該視頻中，直桿隨著時間的推移，高桿的影子也在發(fā)生著變化。由于攝像機(jī)在錄制視頻的過程中，圖像會發(fā)生畸變現(xiàn)象，所以在處理高桿與影子的函數(shù)關(guān)系之前，還必須對畸變的圖像進(jìn)行矯正，將圖像校正后再進(jìn)行處理，來確定該視頻的錄制地點(diǎn)與日期。三、符號與變量說明 符號 說明A 太陽高度角B 觀測地地理緯度C 太陽赤緯角 t 當(dāng)?shù)貙?shí)際時間t0 當(dāng)?shù)厝粘鰰r間t1 當(dāng)?shù)厝章鋾r間J 經(jīng)度N 日期序列數(shù)即積日Y 年份H 旗桿的高

8、度L 旗桿的影長R 太陽的是時角E 日角Tz 時區(qū) W 緯度四、模型假設(shè)1. 假設(shè)地球的自轉(zhuǎn)是勻速的2. 地球近似為一個形狀規(guī)則的球體3. 忽略水蒸氣，透明污染物對光的折射4. 附件所給出的所有點(diǎn)都真實(shí)有效五、模型的建立與求解5.1.問題一的求解5.11相關(guān)概念知識1太陽高度角（A） 對于地球上的某個地點(diǎn)，太陽高度角是指太陽光的入射方向和地平面之間的夾角，專業(yè)上講太陽高度角就是指某地太陽光線與地心相連的地表切線的夾角，記為A。太陽高度也叫高度角。2太陽赤緯角（C）太陽赤緯角，是地球赤道平面與太陽和地球中心的連線之間的夾角，記為W。赤緯角一年周期，在+2326與-2326的范圍內(nèi)移動，成為季節(jié)的

9、標(biāo)志。3太陽時角(R)日面中心的時角，即從觀測點(diǎn)天球子午圈沿天赤道量至太陽所在時圈的角距離，記為R.太陽時角的計(jì)算公式： R=15×（t-12）（1）其中t為真太陽時，以24小時計(jì)。4.太陽高度角的計(jì)算公式： （2） A為太陽高度角，其中B為所觀測地的地理緯度，C為太陽赤緯角，R為太陽時角5.太陽赤緯角的計(jì)算公式： (3) (4) 其中N為積日，E為日角 (5) （6）6.日出時間（t0）的計(jì)算方法： （7)其中J為經(jīng)度，Tz為時區(qū)，W為所觀測地緯度，N為積日。45.1.2模型建立及求解如圖1以直桿的位置為坐標(biāo)原點(diǎn)，東西方向?yàn)閥軸,在南北方向?yàn)閤軸，以xoy平面為地平面，直立方向?yàn)閦

10、軸建立空間直角坐標(biāo)系圖1 坐標(biāo)系從圖中可以看出，直桿OP=h，,影長為，所以根據(jù)三角關(guān)系有，點(diǎn)的坐標(biāo)如下： (9)對（9）式代入（6）式，得到影子端點(diǎn)的坐標(biāo)（x，y）與時間t的關(guān)系 （）根據(jù)（1-9）式和問題一給出的已知條件（2015年10月22日北京時間9:00-15:00之間天安門廣場（北緯39度54分26秒,東經(jīng)116度23分29秒）3米高的直桿）進(jìn)行求解，用Matlab軟件求解得到的坐標(biāo)（計(jì)算程序見附件），其坐標(biāo)滿足二次曲線，繪出圖如下：圖2 影子端點(diǎn)規(guī)律圖與參考文獻(xiàn)1提供的圖作對比，發(fā)現(xiàn)結(jié)論與之吻合由于影子端點(diǎn)的運(yùn)動規(guī)律的曲線已經(jīng)描繪出來了，它的坐標(biāo)軸與L和時間t之間存在著二次曲線的

11、關(guān)系，所以我們采用算出來的影長L和t最小二乘法二次擬合（程序見附件）得到影長L和時間t之間的函數(shù)關(guān)系為： （10）圖3 影長L與時間t的曲線關(guān)系5.2 問題二建模及求解針對第二問已知影子端點(diǎn)的坐標(biāo)（x，y），在此基礎(chǔ)上分析數(shù)據(jù)可知，直桿影子子?xùn)|西方向的投影為y,在南北方向的投影為x，根據(jù)數(shù)據(jù)曲線擬合，預(yù)測地方時中午12點(diǎn)影子數(shù)據(jù)，來計(jì)算緯度，在此基礎(chǔ)上固定緯度和時間，來計(jì)算經(jīng)度。5.2.1緯度的計(jì)算根據(jù)的坐標(biāo)（x，y）的數(shù)據(jù)，用Matlab做最小二乘法曲線擬合（程序見附件）,公式如下 （11）在此基礎(chǔ)上預(yù)測中午12點(diǎn)，影子端點(diǎn)的坐標(biāo)（x，y），由于此時R為太陽時角為R=0，=1，根據(jù)（2）式，

12、合并之后有， ，所以緯度角 （12） 5.2.2經(jīng)度的計(jì)算方法一：在第一問我們得出，如果利用本地的太陽時，則在中午12時直桿的影長達(dá)到最短，我們可以通過擬合出來的圖像來觀察來求經(jīng)度，因?yàn)楦郊薪o出的時間都為北京時間，只要我們算出在那個時間點(diǎn)直桿的影子最短，根據(jù)它跟北京中午的時間差就可以推算出直桿所在的經(jīng)度。根據(jù)根據(jù)的坐標(biāo)（x，y）的數(shù)據(jù)，擬合得到影子長度L與時間t之間的函數(shù) （13）當(dāng) 時，L有最值，也就是影長最短與北京時差與北京經(jīng)度之差：該地經(jīng)度 ： （14）方法二： 根據(jù)（9）式可得 ，由反三角函數(shù)得出 由公式（6）得日出時間 (15)令Q=，由于日期已知，積日知道，上述緯度B計(jì)算出來，Q

13、為常數(shù)，因?yàn)?，設(shè)，根據(jù)（10）式有 （16） 5.2.3計(jì)算結(jié)果 x=0時，y=0.139, C=10.6305 因?yàn)镃在0到90°之間，所以真正的=C-3=1.2057; A=1.5245;聯(lián)立公式（12）解得°°根據(jù)附件1數(shù)據(jù)，擬合出影子長度與時間之間的函數(shù)關(guān)系Linear model Poly2: f(x) = p1*x2 + p2*x + p3Coefficients (with 95% confidence bounds): p1 = 0.1489 (0.1435, 0.1543) p2 = -3.752 (-3.915, -3.589) p3 = 24

14、.13 (22.89, 25.37) 5.3.針對第三問的建模與求解對于問題三，根據(jù)實(shí)際測的影子端點(diǎn)的坐標(biāo)（x，y）的數(shù)據(jù)，不知道具體日期和測驗(yàn)點(diǎn)的經(jīng)緯度。所以問題要比問題二更難，根據(jù)對問題二計(jì)算經(jīng)度的方法，利用附件2的數(shù)據(jù)，擬合出多項(xiàng)式 Linear model Poly2: f(x) = p1*x2 + p2*x + p3Coefficients (with 95% confidence bounds): p1 = 0.09814 (0.09726, 0.09902) p2 = -2.988 (-3.011, -2.965) p3 = 23.37 (23.22, 23.52) 計(jì)算出附件2

15、的實(shí)驗(yàn)點(diǎn)所在位置的經(jīng)度 附件3數(shù)據(jù)擬合多項(xiàng)式Linear model Poly2: f(x) = p1*x2 + p2*x + p3Coefficients (with 95% confidence bounds): p1 = 0.2964 (0.2927, 0.3002) p2 = -7.551 (-7.654, -7.448) p3 = 51.56 (50.86, 52.27)計(jì)算采用窮舉法來探求誤差最小的經(jīng)緯度和日期。日期可以遍歷，所以設(shè)積日N從1到365變動賦值，然后在計(jì)算太陽赤緯角C，得到C(N)序列，在問題二的基礎(chǔ)上，來計(jì)算每一天對應(yīng)的緯度角B(N）序列，同樣根據(jù)問題二計(jì)算經(jīng)度的方

16、法，計(jì)算經(jīng)度J(N).最后根據(jù)N，B(N),C(N)來計(jì)算影子端點(diǎn)的坐標(biāo)（X(N),Y(N)）的數(shù)據(jù)，然后與測量數(shù)據(jù)（x，y）差距最小的那個日期N。 （17） 5.4針對第四問的建模一根直桿在太陽下的影子變化的視頻，并且估計(jì)出直桿的高度為2米。要建立確定視頻拍攝地點(diǎn)的數(shù)學(xué)模型，并應(yīng)用模型給出若干個可能的拍攝地點(diǎn)。對于上述問題，視頻中的物體長度位置與現(xiàn)實(shí)世界中的對應(yīng)物品長度位置是一種射影變換關(guān)系，可以通過抓取視頻的數(shù)據(jù)，建立射影模型，來通過數(shù)據(jù)還原原來真實(shí)物體位置及長度，從而使問題轉(zhuǎn)化為問題二或者問題三的模式。設(shè)視頻中抓取的圖片中物品坐標(biāo)為（x，y），現(xiàn)實(shí)世界中的物品在三維坐標(biāo)系的坐標(biāo)為（X,Y

17、,Z），下面探討轉(zhuǎn)換原理，是一種矩陣變換：1.透過射影幾何法將二維影像坐標(biāo)轉(zhuǎn)換到三維地面坐標(biāo) （18）2. 矩陣為二維影像坐標(biāo)與三維地面坐標(biāo)間的轉(zhuǎn)換矩陣，其內(nèi)容可以透過11個參數(shù)表示， （19）3. 將2-2式中的齊次坐標(biāo)轉(zhuǎn)換到歐氏坐標(biāo)，可得一公式2-3 （20）其中是待定系數(shù)，可以通過多張抓拍照片數(shù)據(jù)來計(jì)算該系數(shù)，然后得到照片中對應(yīng)的物體在三維空間中的坐標(biāo)。5. 計(jì)算流程圖6六模型的評價與推廣模型的優(yōu)點(diǎn)：1. 根據(jù)經(jīng)緯度可以提供不同時節(jié)日出及日落的時間。2. 根據(jù)時間的變化，確定確定太陽光照射的區(qū)域。3. 根據(jù)拍攝的照片確定出該地的經(jīng)度、緯度等地理知識 。4. 將浩瀚的大海區(qū)域化。模型的改進(jìn)

18、：1. 在進(jìn)行角度測量時角度引起誤差的可能性較大。2. 在該模型中，角度的轉(zhuǎn)化很繁瑣，因此出錯的可能很高。3. 用最小二乘法的二次擬合的時候，有些試驗(yàn)點(diǎn)不能夠完全通過，只是一個變化趨勢。參考文獻(xiàn)1鄭鵬飛等，基于影子軌跡反求采光效果的技術(shù)研究J，上海2002.3.72鄭鵬飛等,基于物體影子反求特定時空J(rèn)，華東理工大學(xué)學(xué)報(bào)，上海2002373汪和平，影子與季節(jié)緯度的關(guān)系M，數(shù)學(xué)園地246309.4 太陽高度角_百度百科 5賀曉雷，于賀軍等 太陽方位角的公式求解及其應(yīng)用J，太陽能學(xué)報(bào)，第19卷 第1期。北京100816武琳，基于太陽陰影軌跡的經(jīng)緯度估計(jì)技術(shù)研究D，天津大學(xué)計(jì)算機(jī)科學(xué)與技術(shù)學(xué)院，201

19、0.12附 件1、題一求解程序%數(shù)學(xué)建模A題第一問程序%A為太陽高度角%B為觀測地地理緯度%C為太陽赤緯角%t為當(dāng)?shù)貙?shí)際時間%J為經(jīng)度%N為日期序列數(shù)即積日%Y為年份%h為旗桿的高度%L為旗桿的影長%r為太陽時角%E為日角%Tz為時區(qū)t=9:0.05:15J=116.3914;N=108;Y=2015;h=3;Tz=8;W=39.9072;t0=24*(180+Tz*15-J-acos(-tan(-23.4*cos(2*pi*(N+9)/365)*pi/180)*tan(W*pi/180)*180/pi)/360t1=24*(1+(Tz*15-J)/180)-t0;r=(t-t0)*pi/12

20、;R=(t-12)*pi/12;N0=79.6764+0.2422*(Y-1985)-floor(Y-1985)/4);T0=N-N0;E=2*pi*T0/365.2422;C=0.3723+23.2567*sin(E)+0.1149*2*sin(E)*cos(E)-0.1712*sin(3*E)-0.758*cos(E)+0.3656*cos(2*E)+0.0201*cos(3*E);A=asin(sin(W)*sin(C)+cos(W)*cos(C)*cos(R);L=h./tan(A) x=h./tan(A).*cos(r)y=h./tan(A).*sin(r) plot(x,y)2、用

21、來擬合L和t的函數(shù)關(guān)系的函數(shù)function fitresult, gof = createFit18(t, L) % Fit: 'untitled fit 1'.xData, yData = prepareCurveData( t, L ); % Set up fittype and options.ft = fittype( 'poly2' ); % Fit model to data.fitresult, gof = fit( xData, yData, ft ); % Plot fit with data.figure( 'Name', 'untitled fit