太陽影子定位數(shù)學建模論文

1、太陽影子定位摘 要 太陽影子定位對視頻拍攝分析至關(guān)重要，本文通過建立幾何模型、太陽高度角模型和遍歷模型，繪制出了影長變化曲線，確定了視頻拍攝的地點和日期，解決了太陽影子定位問題。針對問題一，基于地球是球形的假設，建立幾何模型和坐標系，用與影長相關(guān)的五個參數(shù)表示出了桿頂點坐標和光線向量。由點、線、面間的關(guān)系求得直線（光線向量所在直線）和平面（過原點且與桿垂直的平面）的方程，聯(lián)立方程組求得了影子頂點坐標，并用向量的模表示出了影子的長度，由Matlab軟件繪制出了影長變化曲線。最后，基于地球的真實形狀，用太陽高度角模型對模型進行檢驗，驗證了它的合理性。針對問題二，由于它恰好有兩個條件的確定性與問題一

2、相反，所以我們采用問題一精度更高的太陽高度角模型倒推求解。由勾股定理求出21組影長，以影長為縱坐標，北京時間為橫坐標作圖，得到的最低點同時對應北京時間和當?shù)貢r間12:00，根據(jù)問題一中的時差關(guān)系式，反推出當?shù)氐慕?jīng)度，再用遍歷法，求出了合適的緯度，由經(jīng)緯度確定了地點為蒙古自治區(qū)鄂爾多斯市。針對問題三，由于它與問題二的區(qū)別僅是日期未知，所以只需求日期，地點用問題二模型來求。我們采用逆向思維，反推出求日期需要知道五個角，然后正向求解。先通過方位角與影軸角、時角和太陽高度角的關(guān)系建立兩個等式求出方位角，再利用赤緯角與方位角的關(guān)系求出赤緯角，最后利用赤緯角與日期的關(guān)系，建立遍歷模型求出了日期。最后采用問

3、題二的模型求解地點。針對問題四，我們通過處理圖像求得了影長。在日期已知時，我們通過繪制影長-時間圖求解出經(jīng)度為，通過求解高度角與其三個參數(shù)的值，建立太陽高度角模型，通過總關(guān)系式求出了緯度為。 在日期未知時，我們給出了遍歷模型求解緯度，由于經(jīng)度與日期無關(guān)，所以依舊采用日期已知時的方法來求解。本模型考慮了地球不是規(guī)則球體的因素，引入了修正值，使結(jié)果更加可靠，且后續(xù)問題參照前面的模型來求解，使問題大大簡化。關(guān)鍵字：太陽影子定位；幾何模型；太陽高度角模型；遍歷模型1. 問題重述視頻數(shù)據(jù)分析需要確定視頻的拍攝地點和拍攝日期兩大方面，太陽影子定位技術(shù)就是其中一種確定方法，它通過分析視頻中物體太陽影子的變化

4、情況來反推拍攝視頻的地點和日期。為了解決太陽影子的定位問題，我們要解決以下四個問題。1.通過設定中間參數(shù)，建立影子長度變化的數(shù)學模型，并分析出影子長度關(guān)于各個參數(shù)的變化規(guī)律。應用建立的模型畫出2015年10月22日北京時間9:00-15:00之間天安門廣場（北緯39度54分26秒,東經(jīng)116度23分29秒）3米高的直桿的太陽影子長度的變化曲線。2.由直桿在水平面上的影子頂點坐標數(shù)據(jù)，建立數(shù)學模型，確定出直桿所處的地點。將模型應用于附件1的影子頂點坐標數(shù)據(jù)，給出若干個符合條件的地點。3.由直桿在水平面上的影子頂點坐標數(shù)據(jù)，建立數(shù)學模型，確定直桿所處的地點和日期。將模型分別應用于附件2和附件3的影

5、子頂點坐標數(shù)據(jù)，給出若干個符合條件的地點與日期。4.建立確定視頻拍攝地點的數(shù)學模型，并將其應用于附件4，求出若干個符合條件的拍攝地點。附件4為一根直桿在太陽下的影子變化的視頻截圖，并且已通過某種方式估計出直桿的高度為2米。在拍攝日期未知的情況下，判斷是否可以根據(jù)視頻確定出拍攝地點與日期，若可以則給出模型和方法。2模型假設1.假設地球是一個規(guī)則的球形；2.假設附件所給數(shù)據(jù)真實可信；3.假設本文引用的定理公式準確度足夠高。3通用符號說明序號符號符號說明1直桿的長度2影長3經(jīng)度4緯度5太陽赤緯：地球赤道平面與太陽和地球的中心的連線之間的夾角6日期：從1月1日開始計時，每天累計+17太陽高度角8時角4

6、問題一模型的建立與求解4.1問題分析問題一要我們建立影子長度變化的數(shù)學模型，并分析影長關(guān)于各個參數(shù)的變化規(guī)律。由于在基本假設中已經(jīng)假設了地球是一個規(guī)則的球形，我們打算建立幾何模型，采用建立坐標系的方法來求解直桿影子的向量，擬將此坐標系的原點設在坐標原點，赤道設為xoz平面，地球自轉(zhuǎn)軸設為y軸。由于影子的長度與直桿長度、經(jīng)度、緯度和太陽光線直射角度有關(guān)，而太陽直射角度又與日期有關(guān)，我們打算設定五個參數(shù)分別表示它們。我們擬用這五個參數(shù)表示出直桿頂點坐標和太陽光線向量，可能由點、線、面間的關(guān)系表示出太陽光線向量所在直線的方程、過原點且與直桿垂直的平面的方程，聯(lián)立兩個方程可能得到影子頂點的坐標。此坐標

7、到原點的距離就是直桿影子的長度，我們打算用向量的模來表示影長1。根據(jù)影長的公式可能分析出各參數(shù)對影長的影響。我們擬以影長為縱坐標，時間為橫坐標作圖來繪制曲線。4.2模型的建立我們以地心為坐標原點O，設直桿所在經(jīng)線與赤道的交點為B點，以OA所在直線在X軸，地球自轉(zhuǎn)軸為y軸，z軸過圓心且與xoy平面垂直。直桿的影子變化與直桿長度、經(jīng)度、緯度和太陽光線直射角度有關(guān)，而太陽直射角度又與日期有關(guān)，所以我們設定了桿長、直桿所在經(jīng)度、直桿所在緯度、太陽赤緯、時角五個參數(shù)。為了便于求解，我們將直桿向量平移到圓心，此時直桿底端位于圓心，設直桿的頂端A坐標為。在坐標圖中由簡單的幾何關(guān)系可以得到A的坐標與四個參數(shù)間

8、的關(guān)系： （1）其中 （2）T是觀測日期，它從2015年1月1日開始計時，每天+1。 （3）由四個參數(shù)的定義可以得到太陽光線向量為 （4）由式（4）得向量所在直線的斜率為，過點A且斜率為的直線即為從太陽發(fā)射經(jīng)A照射到平面S的光線所在的直線，它的點斜式方程為 （5）設過原點且與OA垂直的平面為平面S，平面S的點法式方程表示為： （6）其中為平面S內(nèi)已知點的坐標，即原點坐標；為平面的法向量即。帶入上述坐標值，可以將式（6）寫為 （7）直線與平面S的交點即為直桿影子的頂點,聯(lián)立方程（5）和（7）即可得到頂點的坐標。求投影點坐標的Matlab程序和運行結(jié)果見附錄。由運行結(jié)果我們得出投影點的坐標為：直桿

9、影子的長度即為向量的模， （8）4.3模型的求解（1）影長與參數(shù)間關(guān)系分析（ki為常數(shù)）我們以其它參數(shù)為常數(shù)，桿長為變量，對影長與桿長間的關(guān)系進行分析，得到可以得出：在其他條件一樣的情況下，影長與桿長呈線性關(guān)系，影長隨桿長的增加而增加。我們以其它參數(shù)為常數(shù)，緯度為變量，對影長和緯度間關(guān)系進行分析，近似得到可以看出影長與緯度近似呈二次函數(shù)關(guān)系，且存在一個最低點，隨緯度的增加影長先減小后增加，最低點對應緯度為0.我們以其它參數(shù)為常數(shù)，日期為變量，對影長和日期間關(guān)系進行分析，近似得到可以看出影長與日期呈二次函數(shù)關(guān)系，有一個最低點，對應冬（夏）至日，隨日期的增大，影長先減小后增大。 同理分析，影長隨時

10、刻的增加呈先減小后增加的趨勢。（2）坐標圖的繪制本題要我們求2015年10月22日北京時間9:00-15:00之間天安門廣場（北緯39度54分26秒,東經(jīng)116度23分29秒）3米高的直桿的太陽影子長度的變化曲線，我們由題目可以讀出直桿長度、經(jīng)度、緯度和日期，由模型一可以求出太陽赤緯。表1 問題一已知條件表符號數(shù)值直桿的長度2米直桿所在經(jīng)度東經(jīng)116度23分29秒直桿所在緯度北緯39度54分26秒日期2015年10月22日，即295天將上述已知條件帶入公式（2）和（3），用Matlab軟件求出和 并將坐標和參數(shù)的具體數(shù)值帶入向量的模的公式（8），然后繪圖，程序見附錄1。我們得到下圖： 圖1 問

11、題一影長-北京時間圖由圖得，在北京時間9:00-15:00間，影子長度呈先減小后增加的趨勢，二者符合二次函數(shù)關(guān)系，圖像存在一個最低點，此點對應北京時間12:00，由于直桿位于東八區(qū)，所以12:00也是最低點對應的當?shù)貢r間。4.4模型的檢驗由于模型一是基于地球是一個規(guī)則的球形的假設來計算的，存在誤差，為了消除誤差，我們考慮到地球的實際形狀，建立了太陽高度角模型2，來檢驗模型一。直桿、直桿影子和太陽光線可以構(gòu)成一個直角三角形，令太陽高度角為，直桿長為，直桿影長為，則它們之間的關(guān)系為：太陽高度角與四個因素有關(guān)，分別是太陽赤緯、緯度、和時角有關(guān)，它們之間的關(guān)系為：太陽赤緯的計算公式為：其中，T為觀測日

12、期，從每年1月1日開始計算，每天數(shù)值累計+1。本題要我們求2015年10月22日對應T=295，將T帶入上式得：緯度的計算公式為時角的計算公式為 （9）其中為真太陽時；為北京時間；為北京時間與當?shù)貢r間的時差；為天數(shù)，從1月2日開始計時；為修正參數(shù)，由于地球不是一個規(guī)則的球形，所以需要修正。是當?shù)鼐暥龋捎诒本r間指東京的地方時間，有緯度差，而每經(jīng)度變化，時間相差1小時。用上述模型求出9:00時的影長為5.67米，而模型一求出此時的影長為5.13米，偏差為-9.53%，偏差不大，所以模型一合理可信。5問題二模型的建立與求解5.1問題分析 考慮到問題二恰好有兩個條件的確定性與問題一相反，我們打算采

13、用問題一的模型倒推經(jīng)緯度，來確定地點。我們擬通過勾股定理對附件1中的影子的坐標進行處理，求21組直桿影子在對應北京時間的長度。以影長為縱坐標，北京時間為橫坐標作圖，可能得到一個最低點，這個點對應一個北京時間和當?shù)貢r間（12:00）。我們準備根據(jù)問題一模型中的北京時間與當?shù)貢r間的關(guān)系式，反推出當?shù)氐慕?jīng)度。對于緯度的確定，我們打算采用遍歷的方法，將連續(xù)變化的緯度離散化，可能找出直桿的方差符合要求的緯度。經(jīng)度和緯度求出后，排除出不合適的地點，正確地點就可以確定了。5.2模型的建立（1）影長的求法：用勾股定理求出21組直桿影子在對應北京時間的長度。（2）經(jīng)度的求法：由問題一知影長與時間呈拋物線，以影長

14、為縱坐標，北京時間為橫坐標作圖，可以得到一個最低點，這個點對應一個北京時間，此點同時對應當?shù)貢r間12:00。由于經(jīng)度每差，時間相差一小時，我們得到北京時間與當?shù)貢r間的關(guān)系為北京時間，為當?shù)貢r間。由于已知，所以經(jīng)度可求。（3）太陽高度角的求法：太陽高度角公式為根據(jù)時角公式上述各字母意義與問題一相同。將求出的數(shù)值帶入上式，即可求出太陽高度角與緯度的關(guān)系式。（4）緯度的求法：我們用遍歷法來確定緯度。以北緯為正，南緯為負，緯度區(qū)間為。將連續(xù)變化的緯度離散成每一緯度一個值進行遍歷，確定出一個可能的緯度區(qū)間。每一個可能的緯度與經(jīng)度結(jié)合可以確定出一個當?shù)貢r間。由問題一可以得到直桿與其影子的關(guān)系為由于對一根桿

15、來說，它的高度一定，所以我們求這21組直桿長度的方差，對應方差最小的緯度是最優(yōu)解。為組直桿長度的平均值。進行方差檢驗求出可能區(qū)間。然后在可能區(qū)間內(nèi)，細化步長，按上述方法再次求解，直到得出符合要求的結(jié)果為止。我們確定目標函數(shù)如下：綜上所述，問題二的模型為 5.3模型的求解（1）經(jīng)度的求解求解經(jīng)度的Matlab程序見附錄2，運行結(jié)果圖如下：圖2 問題二影長-時間圖由上圖得影長與時間呈二次函數(shù)關(guān)系，由系數(shù)可求最低點對應的時間為12:29，經(jīng)度為110.23度。（2）緯度的求解編寫Matlab程序及運行結(jié)果見附錄2，由結(jié)果我們得到方差較小的，滿足條件的可能區(qū)域為：表2 問題二緯度求解表緯度343536

16、桿長1.296371848121611.294834288103711.29322806743802緯度373839桿長1.289811295328511.28800160666087緯度404142桿長緯度434445桿長1.280098230862491.277958666344341.27575463247034由上表我們可以得到可能的緯度區(qū)間為北緯。我們再次運用遍歷法，將此區(qū)間進一步細化，選擇直桿方差最小的一組數(shù)據(jù)對應的緯度，確定此緯度為北緯40度。綜上所述，該圖像拍攝的地點為東經(jīng)110.23度，北緯40度，此地點對應內(nèi)蒙古自治區(qū)鄂爾多斯市。6問題三模型的建立與求解6.1問題分析問題三

17、在問題二的基礎上增加了日期變量，我們考慮先建立模型求出日期，地點的求解與問題二一致。此問題就轉(zhuǎn)化為了建立模型求日期問題。由于日期與赤緯角存在函數(shù)關(guān)系，我們打算通過求赤緯角來確定日期。我們擬引入新的變量方位角3，通過幾何關(guān)系，可能得出方位角與赤緯角的關(guān)系。我們打算用太陽高度角和時角來表示方位角。我們準備引入影長與x軸的夾角來表示方位角，此夾角對日期的變化率等于方位角隨日期的變化率，由此可建立一個等式。由于方位角又可以由太陽高度角和時角來表示，由方位角的兩種表示方法可以建立等式方程。通過這種方法可能求解出日期。6.2模型的建立日期T的求法：我們通過赤緯角求日期，赤緯角通過新變量方位角來確定。太陽方

18、位角是指太陽光線在地平面上的投影與當?shù)亟?jīng)線的夾角，可近似看作是豎立在地面上的直線在陽光下的陰影與正南方的夾角s， 對于方位角的確定，我們通過定義影長與x軸的夾角來表示。太陽方位角表示為，i表示附件2第i組數(shù)據(jù)。由于和關(guān)于時間的變化率相等。我們得到方位角還可以由太陽高度角和時角來表示由問題一的公式（9）可以求出時角t。由于日期是連續(xù)的，不利于我們逐點分析，所以我們采用遍歷法對日期以1天為步長進行遍歷。方法類似問題二，這里不再贅述。綜上所述，我們建立模型如下： 6.3模型的求解（1）日期的求解求解日期的Matlab程序見附錄3，由運行結(jié)果我們得到：附件2 數(shù)據(jù)對應的日期為172和355，分別對應6

19、月21日和12月21日，即冬至和夏至日。附件3 對應的日期為172和355，分別對應6月21日和12月21日。（2）經(jīng)度的求解圖3 附件2影長-時間圖由上圖，我們得到：對于附件2,我們得到曲線的三個系數(shù)分別為0.000245，-0.0253,1.2725，由系數(shù)我們求出最低點對應的北京時間為15:16，經(jīng)度差為49度，所以所求地點的位置位于東經(jīng)度71度。圖4 附件3影子-時間圖由上圖，我們得到：對于附件3,我們得到曲線的三個系數(shù)分別為0.0007411，0.0108, 3.5227，由系數(shù)我們求出最低點對應的北京時間為12:47，經(jīng)度差為11.75度，所以所求地點的位置位于東經(jīng)度108.25度

20、。（3）緯度的求解緯度的求解同問題二，附件2可能的緯度為表3 附件2緯度求解表緯度日期6月21日12月21日附件3可能的緯度為表4 附件3緯度求解表緯度日期6月21日12月21日（4）地點和日期的求解結(jié)果我們將緯度與經(jīng)度結(jié)合，在地圖中尋找地點，由于存在某些位于海上的點，不符合我們的要求，所以我們將其舍棄。最后得到符合要求的地點。對于附件2，在6月21日，符合的地點為Jalal-Abad，對應東經(jīng)71度，北緯42度；在12月21日，符合的地點為Karagandy Region，對應東經(jīng)71度，北緯48度。對于附件3，在6月21日，符合的地點為Burytia,Sibirsky federating

21、 okruh，對應東經(jīng)108度，北緯53度；在12月21日，符合的地點為Erdene,Tov，對應東經(jīng)108度，南緯48度。7問題四模型的建立與求解7.1問題分析 此題給出的是直桿在太陽下的影子變化的視頻截圖，并且已通過某種方式估計出直桿的高度為2米，我們打算用Matlab來處理圖像，應該可以根據(jù)桿長求出影長。當日期已知時，此題與問題二相似，我們擬采用問題二的作圖法，結(jié)合時差與經(jīng)度的關(guān)系式來求解經(jīng)度。對于緯度的求解，由于在問題一的檢驗模型中我們知道了太陽高度角與包含緯度在內(nèi)的四個變量之間的總關(guān)系式，所以我們打算建立太陽高度角模型，準備求解出高度角與其他三個變量，應該可以通過總關(guān)系式得出緯度來。

22、 當日期未知時，經(jīng)度的求法不涉及日期，所以上一小問的求法依舊可以用。而緯度的求解時，由于上一小問用到了太陽赤緯，它隨日期而變，所以我們需要采用其它辦法求解赤緯，若求出赤緯，則此題與上一題變得相同，可以采用上一小問的方法求解。結(jié)合問題三，我們考慮采用遍歷模型，對緯度進行遍歷，尋找合適的緯度，求出赤緯角。7.2模型的建立日期已知時（1）經(jīng)度的求法：由問題一知影長與北京時間呈拋物線，以影長為縱坐標，北京時間為橫坐標作圖，可以得到一個最低點，這個點對應一個北京時間，此點同時對應當?shù)貢r間12:00。由北京時間與當?shù)貢r間的關(guān)系式 我們可以求出經(jīng)度。（2）緯度的求法：由問題一、二，我們得到如下關(guān)系：太陽高度

23、角與太陽赤緯、緯度、時角的關(guān)系為： 時角的計算公式為：太陽赤緯的計算公式為：太陽高度角與直桿長、影長的關(guān)系：我們通過圖像處理，根據(jù)上述公式可以得出、和的大小來，由太陽高度角與太陽赤緯、緯度、時角的關(guān)系式，我們可以求出緯度的大小。結(jié)合經(jīng)緯度，可求出地點。7.2.2日期未知時經(jīng)度的求解方法同7.2.1。緯度的求解：參考上一小問的求解過程，我們使用如下關(guān)系式： 將緯度以為步長進行遍歷，遍歷范圍為，由每一個緯度，我們可以對應求出一個太陽赤緯角，計算方差值我們保留方差值小的赤緯，由此確定出可能的區(qū)域。再對此區(qū)域進行遍歷，將步長細化十倍，直到方差足夠小為止，我們就可以得到一個合適的緯度區(qū)間。綜合分析緯度與

24、經(jīng)度，我們可以得到合適的地點。7.3模型的求解由附件4求得15張圖片的影長為表5 附件4影長表序號123影長2.389743589743592.32991452991453序號456影長2.264957264957272.22222222222222序號789影長序號101112影長2.037606837606842.003418803418801.95726495726496序號131415影長1.935042735042741.890598290598291.85982905982906結(jié)合上表數(shù)據(jù)，我們以影長為縱坐標，北京時間為橫坐標作圖圖5 問題四影長-時間圖我們由此圖的最低點求得經(jīng)度

25、為東經(jīng)102.165度。然后由附錄4求解緯度的Matlab程序，求得緯度為北緯58.1度。結(jié)合地圖我們得到此處為Irkutsk Oblast,Sibirsky federalny okruh.8模型的評價8.1模型的優(yōu)點1.針對幾何模型，通過建立坐標系，合理的將影子長度問題轉(zhuǎn)化為了向量模的問題，思路清晰，便于理解，運用了點、線、面之間的關(guān)系，很好的解決了影長的求算問題。2.針對太陽高度角模型，它基于地球的真實形狀來考慮，引入了修正值，可信度十分高。3.針對遍歷模型，它將連續(xù)變化的緯度離散化，逐點求解，極大的縮小了緯度區(qū)間，再通過細化步長，進一步縮小緯度區(qū)間，給出了合理的緯度。8.2模型的缺點由

26、于附件所給的數(shù)據(jù)太少，使我們的作圖過程可用信息太少，圖的準確度達不到極高的水平。8.3模型的推廣 本文的太陽高度角模型可以推廣到建筑領域，對于建筑的采光問題十分有效。同時，我們可以將直桿換成其他任何難以求高度的物體，利用本文模型來求算其高度。參考文獻1 高等數(shù)學（下）M，北京：高等教育出版社，38-39，2007.62 胡毅華，楊旭龍，劉媛萍，太陽影子定位模型的構(gòu)建J, 洛陽師范學院學報，Vol.34, No.11:13-18，2015.113 鄭鵬飛，基于影子軌跡線反求采光效果的技術(shù)研究J，華東理工大學學報: 自然學科版，Vol.36，No .3 : 458-463，2010附錄附錄11.求

27、投影點坐標的Matlab程序（由于Matlab程序不能識別和，我們用w和c分別代替）syms x y z w r c l t;eq1=(x-l*cos(w)*cos(t)*tan(c)-(y-l*sin(w);eq2=z-l*cos(w)*sin(t);eq3=l*x*cos(w)*cos(t)+l*y*sin(w)+l*z*cos(w)*sin(t);x y z=solve(eq1,eq2,eq3,x,y,z)2.運行結(jié)果x=-(l*sin(w)2+l*cos(w)2*sin(t)2-l*tan(c)*cos(t)*cos(w)*sin(w)/(cos(t)*cos(w)+tan(c)*si

28、n(w)y=-(cos(w)*(l*tan(c)*cos(t)2*cos(w)-l*cos(t)*sin(w)+l*tan(c)*cos(w)*sin(t)2)/(cos(t)*cos(w) + tan(c)*sin(w)z =l*cos(w)*sin(t)3.畫圖像的Matlab程序w=40/180;l=3;t=-0.8:0.01:0.8;c=0;for i=1:161x(i)=-(l*sin(w)2+l*cos(w)2*sin(t(i)2-l*tan(c)*cos(t(i)*cos(w)*sin(w)/(cos(t(i)*cos(w) + tan(c)*sin(w);y(i)=-(cos(

29、w)*(l*tan(c)*cos(t(i)2*cos(w)-l*cos(t(i)*sin(w)+l*tan(c)*cos(w)*sin(t(i)2)/(cos(t(i)*cos(w) + tan(c)*sin(w);z(i)=l*cos(w)*sin(t(i);H(i)=(x(i)2+y(i)2+z(i)2)0.5;endT=9:0.0375:15;plot(T,H,'r-')ylabel('高度/米');xlabel('時間')附錄21.求經(jīng)度的Matlab程序for i=1:21 y(i)=(S1(i,1)2+S1(i,2)2)0.5;endx=1:21;A=polyfit(x,y,2);y1=polyval(A,x);2.求緯度的Matlab程序B=0.048;for c=-90:90 for i=1:21 t(i)=0.262*(2.7+0.05*i); end for i=1:21 A1=sin(B)*sin(c/180)+cos(B)*cos(c/180)*cos(t(i); af=asin(A1); H(i)=