太陽影子定位

1、太陽影子定位摘要對于太陽影子的定位問題，本文結(jié)合地理學(xué)與天文學(xué)的相關(guān)知識，建立了適用于不同因素影響的太陽影子定位模型，實(shí)現(xiàn)了對于拍攝地點(diǎn)和日期的快速確定。針對問題一，建立直桿影子長度變化的數(shù)學(xué)模型，通過分析影子長度關(guān)于經(jīng)度，緯度等參數(shù)的變化規(guī)律以及與太陽高度角、太陽赤緯和太陽時角等中間因素存在的聯(lián)系，建立起數(shù)學(xué)模型。根據(jù)所給的日期、時間、地點(diǎn)、桿高等數(shù)據(jù)，使用模型繪制位于天安門，2015年10月22日的直桿影子長度的變化曲線，得出最小長度出現(xiàn)在正午12時，最短長度約為3.68米。針對問題二，建立平面坐標(biāo)，根據(jù)給定的平面影子頂點(diǎn)坐標(biāo)，用最小二乘法思想進(jìn)行曲線擬合，得出影長隨時間的變化規(guī)律。通過幾

2、何原理計(jì)算太陽高度角、赤緯角、地理緯度、時角和太陽方位角之間的關(guān)系，計(jì)算出目標(biāo)地的地理緯度。同時選取北京正午十二時作為一個參考值，聯(lián)系擬合的曲線得出目標(biāo)地太陽影子最短的時刻，與參考值進(jìn)行對比得出經(jīng)度。經(jīng)緯度的雙鎖定確定直桿所在地的具體位置大致是廣東省茂州市。針對問題三，建立與第二問類似的目標(biāo)規(guī)劃模型，但是加上日期因素，為了提高遍歷速度，使用粒子群算法減小時間復(fù)雜度。計(jì)算得出的結(jié)果：附件二地點(diǎn)坐標(biāo)（80.34°E,33.32°N），8月13日或4月27日，約在西藏山南地區(qū)；附件三地點(diǎn)坐標(biāo)（113.21°E，23.54°N），12月26日或12月12日，約在

3、廣西河池市。針對問題四，對視頻進(jìn)行處理以得到必須的數(shù)據(jù)，由于拍攝時角度的存在，得到的影長并非實(shí)際長度，而是投影長度，所以采用基于Hough變換和透視變換的圖像矯正法，對斜視圖像矯正，得出實(shí)際長度，然后將得到的數(shù)據(jù)代入第二問中的模型，得出的結(jié)果是視頻拍攝地點(diǎn)在內(nèi)蒙古包頭市境內(nèi)；在拍攝日期未知的情況下，將圖像變換得出的實(shí)際影長代入第三問中的基于粒子群算法的規(guī)劃模型里，可得出日期為6月11日和7月13日. 對于模型的推廣，根據(jù)某地的陽光信息與地理信息，可以應(yīng)用到建筑行業(yè)日照設(shè)計(jì)，設(shè)計(jì)住宅布置、建筑間距、房屋朝向等方面，甚至還可以應(yīng)用到農(nóng)業(yè)大棚養(yǎng)殖安排蔬菜的種植位置等。 關(guān)鍵詞：二次函數(shù)擬合 最小二乘

4、法 粒子群算法 圖形矯正法 目錄一、問題重述1.1問題的背景1.2問題的提出二、模型假設(shè)三、符號說明四、問題分析4.1 問題一分析4.2問題二分析4.3問題三分析4.4 問題四分析五、模型建立與求解5.1 問題一模型建立與求解5.1.1 概念引入 5.1.1.1 太陽高度角5.1.1.2 太陽赤緯5.1.1.3 太陽時角5.1.2 模型建立 5.1.3 模型求解5.2 問題二模型建立與求解5.2.1 概念引入5.2.2 模型建立5.2.2.1 緯度求解模型5.2.2.2 經(jīng)度求解模型5.2.3 模型求解5.3問題三模型建立與求解5.3.1 模型的建立與求解5.3.2粒子群算法目標(biāo)規(guī)劃模型5.3

5、.3粒子群算法的簡介5.3.4粒子群算法模型的求解過程5.3.5 模型的分析比較5.4問題四模型建立與求解 5.4.1 概念引入5.4.2 模型建立 5.4.3 模型求解5.4.4 結(jié)果分析六、誤差分析與模型檢驗(yàn)6.1 誤差分析6.2 模型的檢驗(yàn)七、模型評價(jià)八、參考文獻(xiàn)附錄一、問題重述1.1問題的背景 對于直立在地面的直桿，其影長與方向會隨著日光方向的變化而改變，而且與一天中的時刻有關(guān)，因此對太陽光影的變化進(jìn)行研究以獲得時間和地理信息是可行的，特別是使用對圖像的處理方法，獲得拍攝地點(diǎn)的各種信息，對采光，建筑等工程提出改進(jìn)，具有現(xiàn)實(shí)應(yīng)用的意義。1.2問題的提出如何確定視頻的拍攝地點(diǎn)和拍攝日期是視

6、頻數(shù)據(jù)分析的重要方面，太陽影子定位技術(shù)就是通過分析視頻中物體的太陽影子變化，確定視頻拍攝的地點(diǎn)和日期的一種方法。1.建立影子長度變化的數(shù)學(xué)模型，分析影子長度關(guān)于各個參數(shù)的變化規(guī)律，并應(yīng)用你們建立的模型畫出2015年10月22日北京時間9:00-15:00之間天安門廣場（北緯39度54分26秒,東經(jīng)116度23分29秒）3米高的直桿的太陽影子長度的變化曲線。2.根據(jù)某固定直桿在水平地面上的太陽影子頂點(diǎn)坐標(biāo)數(shù)據(jù)，建立數(shù)學(xué)模型確定直桿所處的地點(diǎn)。將你們的模型應(yīng)用于附件1的影子頂點(diǎn)坐標(biāo)數(shù)據(jù)，給出若干個可能的地點(diǎn)。3. 根據(jù)某固定直桿在水平地面上的太陽影子頂點(diǎn)坐標(biāo)數(shù)據(jù)，建立數(shù)學(xué)模型確定直桿所處的地點(diǎn)和日

7、期。將你們的模型分別應(yīng)用于附件2和附件3的影子頂點(diǎn)坐標(biāo)數(shù)據(jù)，給出若干個可能的地點(diǎn)與日期。4附件4為一根直桿在太陽下的影子變化的視頻，并且已通過某種方式估計(jì)出直桿的高度為2米。請建立確定視頻拍攝地點(diǎn)的數(shù)學(xué)模型，并應(yīng)用你們的模型給出若干個可能的拍攝地點(diǎn)。如果拍攝日期未知，你能否根據(jù)視頻確定出拍攝地點(diǎn)與日期？二、模型假設(shè)假設(shè)一：假設(shè)地球?yàn)橐粋€球體，忽略實(shí)際橢圓體對解題的影響； 假設(shè)二：忽略海拔、天氣等因素對測量和計(jì)算造成的影響； 假設(shè)三：假設(shè)陽光照射到地球上時為平行光； 假設(shè)四：假設(shè)直桿本身理想筆直，無彎曲情況；假設(shè)五：忽略行星軌道差異； 假設(shè)六：太陽影子的長度變化只與太陽照射有關(guān)，不考慮大氣對太陽

8、的折射；三、符號說明符號含義解釋 太陽高度角太陽相對于地面觀察點(diǎn)的仰角 太陽時角子午圈沿赤道方向至太陽時圈的角度緯度角觀察點(diǎn)所在的緯度 赤緯角太陽直射點(diǎn)的緯度 太陽方位角從正北方起順時針到太陽時圈的角度（m,n） 正視坐標(biāo)第四問視頻經(jīng)變換后從俯視角度觀察呈現(xiàn)的坐標(biāo)值四、問題分析 太陽影子定位技術(shù)就是利用投影物體在太陽光的照射，根據(jù)影子長度的變化規(guī)律，確定該地的所處緯度位置、所處的時間以及所處的位置坐標(biāo)，進(jìn)而求得準(zhǔn)確地點(diǎn)。4.1 問題一分析問題一要求建立反映影子長度變化的數(shù)學(xué)模型，引入太陽高度角、太陽赤緯和太陽時角概念，根據(jù)三者之間的關(guān)系，利用三角函數(shù)關(guān)系和幾何關(guān)系，建立不同時刻太陽下影子長度的

9、變化的數(shù)學(xué)模型，并根據(jù)已給的日期、時間、地點(diǎn)、桿高數(shù)據(jù)，應(yīng)用所建立的模型畫出影子長度變化曲線。4.2 問題二分析問題二給出了一組數(shù)據(jù)，要求建立關(guān)于直桿所處的地點(diǎn)的數(shù)學(xué)模型，并應(yīng)用于附件1中影子頂點(diǎn)坐標(biāo)數(shù)據(jù)，得出直桿所在地點(diǎn)?；谥睏U影子與時間之間的關(guān)系，建立時間與直桿影子頂點(diǎn)坐標(biāo)的數(shù)學(xué)模型。同時利用地理上的時區(qū)計(jì)算方法，來進(jìn)一步求出該地的經(jīng)緯度，進(jìn)而通過經(jīng)緯度的雙重鎖定確定所求的地點(diǎn)。     對于緯度求解，我們在問題一中已經(jīng)有了相應(yīng)的說明，找出太陽方位角、太陽高度角、時角和太陽赤緯角之間的關(guān)系即可求出。根據(jù)經(jīng)線的劃分原則，則可知地球自轉(zhuǎn)每分鐘

10、為0.25度，然后結(jié)合問題一中的方法計(jì)算出經(jīng)度。根據(jù)所建立的模型，利用最小二乘法的思維，利用MATLAB軟件擬合出直桿影長與時間的關(guān)系圖像。 4.3 問題三分析對于問題三，可視為已知投影頂點(diǎn)的坐標(biāo)和投影時刻，要求得出直桿所在位置的經(jīng)緯度與拍攝的日期。首先，將附件2與3中的影長和時刻數(shù)據(jù)帶入問題一求得的影長變化函數(shù)中，得到方程組。其次該問題是一個雙目標(biāo)最優(yōu)化問題，可以建立以影長差值最小，太陽方位角差值最小為目標(biāo)的雙目標(biāo)最優(yōu)化模型對附件2,3 的數(shù)據(jù)進(jìn)行求解，得出位置的經(jīng)緯度與拍攝的日期，其中約束條件考慮經(jīng)緯度范圍。由此，可建立雙目標(biāo)最優(yōu)化模型，利用粒子群算法優(yōu)化遍歷過程，并對模型的正確性進(jìn)行檢驗(yàn)

11、。4.4 問題四分析對問題四所給的視頻附件進(jìn)行一定的預(yù)處理，由于給出視頻的清晰度不能滿足要求，還需要對圖片進(jìn)行灰度處理，以達(dá)到減少計(jì)算量和提高計(jì)算精確度的目的。利用MATLAB對處理后的數(shù)據(jù)進(jìn)行分析，與問題二的模型相互聯(lián)系，得出直桿影子長度變換規(guī)律，繼而得出可能的拍攝地點(diǎn)。 五、模型建立與求解5.1 問題一模型建立與求解 5.1.1 概念引入 5.1.1.1 太陽高度角對地球上所選取的某個地點(diǎn)來說，太陽高度角是指太陽光照射的方向和地平面之間的夾角，從專業(yè)術(shù)語上講太陽高度角是指某一個地方太陽光線與通過該地與地心相連的地表切線的夾角。太陽高度角隨著地方時間和太陽的赤緯的變化而變化。而且當(dāng)太陽高度角

12、為90度時，此時太陽輻射強(qiáng)度是最大的；當(dāng)太陽斜射在地面時，太陽輻射強(qiáng)度就變小。5.1.1.2 太陽赤緯太陽赤緯又稱赤緯角，是地球赤道平面與太陽和地球中心的連線之間的夾角。赤緯角以年為周期，在+23°26與-23°26的范圍內(nèi)移動，成為季節(jié)的標(biāo)志。（圖1） 圖1 赤緯角的定義 5.1.1.3 太陽時角日面中心的時角，即從觀測點(diǎn)天球子午圈沿天赤道量至太陽所在時圈的角距離。以地球?yàn)槔诘厍蛏?，同一時刻，對同一經(jīng)度，不同緯度的人來說，太陽對應(yīng)的時角是相同的。單位時間地球自轉(zhuǎn)的角度定義為時角w，規(guī)定正午時角為0，上午時角為負(fù)值，下午時角為正值。地球自轉(zhuǎn)一周360度，對應(yīng)的時間為24

13、小時，即每小時相應(yīng)的時角為15度。5.1.2 模型建立 引入了太陽高度角、太陽赤緯、太陽時角的相應(yīng)概念，根據(jù)三者之間的定義，我們可以得出太陽高度角隨著太陽赤緯和太陽時角的變化而變化的結(jié)論。根據(jù)三者的空間幾何聯(lián)系，建立合適的模型來反映影子長度變化。（下圖2） 圖2：高度角、赤緯角、時角關(guān)系圖已知太陽高度角隨太陽赤緯和太陽時角的變化而變化，三者之間存在以下的運(yùn)算關(guān)系： (1)其中為太陽高度角，為直射點(diǎn)地理緯度，為太陽赤緯角，為太陽時角。通過查閱資料，我們發(fā)現(xiàn)，因太陽赤緯值日變化很小，所以一年內(nèi)任意一天的赤緯角可以通過下式進(jìn)行計(jì)算： (2)其中，N為日數(shù)，自每年第一天計(jì)算。 (3)其中，T為時刻，自

14、每日凌晨計(jì)算。將查閱到的太陽赤緯角以及太陽時角的的計(jì)算方法代入太陽高度角的計(jì)算公式，可以得到以下的方程: (4) 根據(jù)題目要求畫出北京天安門廣場的3米高的直桿的太陽影子直射簡易圖： 圖3：太陽影子直射圖陽光與地面的夾角即為此地的太陽高度角，AC為目標(biāo)桿，桿長為3m，在三角形ABC中，由正切定理可知： (5)要想求出直桿AC的影長BC隨時間的變化規(guī)律，既要求出該地太陽直射角在指定時間內(nèi)的變化規(guī)律，即可得出影長： (6)由（4）式可知太陽直射角： (7) (8)5.1.3 模型求解已知目標(biāo)對象為2015年10月22日的北京天安門廣場（北緯39度54分26秒，東經(jīng)116度23分29秒）。根據(jù)以上信息

15、可得地理緯度為39度54分26秒，日數(shù)N為295天，并且根據(jù)目前已有信息可求得該地點(diǎn)赤緯度數(shù)。以時刻為單位，并且聯(lián)立(1)-(8)式，求得各個時刻的影長數(shù)據(jù)，通過編寫程序以及MATLAB軟件的應(yīng)用，畫出對應(yīng)的影長圖像，從而得出太陽直射的直桿影子長度變化曲線。（圖4） 圖4：直桿影長變化圖根據(jù)MATLAB所得圖像，可知2015年10月22日北京時間9:00-15:00之間天安門廣場（北緯39度54分26秒,東經(jīng)116度23分29秒）3米高的直桿的太陽影子長度的變化曲線為一開口向上的拋物線，直桿影子長度隨時間呈現(xiàn)出先減小后增大的規(guī)律，且最小長度出現(xiàn)在正午12時，最短長度約為3.68米。5.2 問題

16、二模型建立與求解5.2.1 概念引入太陽方位角：即太陽所在的方位，指太陽光線在地平面上的投影與當(dāng)?shù)刈游缇€的夾角，可近似地看作是豎立在地面上的直線在陽光下的陰影與正南方的夾角。方位角以目標(biāo)物正北方向?yàn)榱悖槙r針方向逐漸變大，其取值范圍是0360°。因此太陽方位角一般是以目標(biāo)物的 北方向?yàn)槠鹗挤较?，以太陽光的入射方?為終止方向，按順時針方向所測量的角度。5.2.2 模型建立已知目標(biāo)桿各個時刻于水平地面上在太陽下的太陽影子頂點(diǎn)坐標(biāo)，可以得出相對應(yīng)時刻的太陽影長，進(jìn)而得到在指定時間內(nèi)的目標(biāo)桿太陽影子長度變化規(guī)律，聯(lián)系幾何原理以及太陽光射角度等變量，得出直桿在該地太陽影子長度的最小值，聯(lián)系一

17、定的天文知識，計(jì)算出該地的經(jīng)度和緯度，繼而查找出該地所在。（圖5） 圖5;太陽直射下直桿影子投射示意圖 以直桿所在點(diǎn)為目標(biāo)原點(diǎn)，以東西方向?yàn)閄軸，南北方向?yàn)閅軸，建立XY平面坐標(biāo)系。模擬出太陽直射下直桿影子在平面上的投射，此時影子頂點(diǎn)為S點(diǎn)，其在XY坐標(biāo)軸上的投影分別為Xs，Ys。（圖中角為太陽高度角，角為太陽方位角）5.2.2.1 緯度求解模型日升日落，同一地點(diǎn)一天內(nèi)太陽高度角是不斷變化的。日出日落時角度都為0，正午時太陽高度角最大，時角為0，以上的公式可以簡化為： (9)由兩角和與差的三角函數(shù)公式，可得 (10)此時，發(fā)現(xiàn)存在兩種情況，因此對太陽位于天頂以北的地區(qū)和太 陽位于天頂以南的地區(qū)

18、分別進(jìn)行討論：對于太陽位于天頂以北的地區(qū)而言，=90º()；對于太陽位于天頂以南的地區(qū)而言，=90º()；二者合并，因?yàn)闊o論是()還是()，都是為了求當(dāng)?shù)鼐暥扰c太陽直射緯度之差，不會是負(fù)的，因此都等于它的絕對值，所以得到正午太陽高度角計(jì)算公式： (11)因此，只需要知道該地的太陽高度角和太陽赤緯即可得出該地的地理緯度，太陽赤緯角可以由式（2）計(jì)算得出。引入太陽方位角，要得出太陽高度角的值，經(jīng)查閱查閱推導(dǎo)得出太陽方位角、太陽高度角、太陽赤緯角以及太陽時角之間滿足下列關(guān)系： (12) 經(jīng)過（11）（12）式的聯(lián)立，可以得到以下的等式： (13) (14)在上列處理中，為太陽高度

19、角，為直射點(diǎn)地理緯度，為太陽赤緯度, 為太陽方位角，為太陽時角。聯(lián)立(2)(3)(14)即可得出該地的地理緯度。5.2.2.2 經(jīng)度求解模型地球自傳一周需要24小時，即將全球被劃分成24個時區(qū)，可知地球每小時自傳15度，即可推出每分鐘地球自轉(zhuǎn)了0.25度，根據(jù)這個數(shù)據(jù)既然可以計(jì)算出不同經(jīng)度地方的時間差，反之亦可從時間差推導(dǎo)出經(jīng)度差。通過問題一的求解，我們已經(jīng)知道，直桿的影子最短發(fā)生在一天中正午時，通過選取已知緯度和正午時刻的地方作為參考地，利用正午時刻時間差即可求出所求地的經(jīng)度值。建立目標(biāo)地點(diǎn)對比點(diǎn)，選取北京正午時間t。和經(jīng)度。作為參考點(diǎn)（其中t。為正午12:00, 。為東經(jīng)120 º

20、;），所求地正午時刻t1與t。對比計(jì)算出兩地的時間差，再通過數(shù)學(xué)計(jì)算把差值換算成相應(yīng)的度數(shù)，即可以得出所求地的經(jīng)度值1。兩者有如下關(guān)系： （ t1<t。） (15) （ t1>t。） (16)為了得出附件1所提供的數(shù)據(jù)反映的可能地點(diǎn)，需要求出該地的經(jīng)緯度。因此還要對附件1的數(shù)據(jù)進(jìn)行模擬得出目標(biāo)直桿達(dá)到影長最短的時刻，即可求出該地的經(jīng)度數(shù)。5.2.3 模型求解為了得出問題中直桿所處的地點(diǎn)，我們需要先對附件1直桿的影子頂點(diǎn)坐標(biāo)數(shù)據(jù)進(jìn)行預(yù)處理，對給出的數(shù)據(jù)在應(yīng)用到XY坐標(biāo)系的基礎(chǔ)上，通過對兩坐標(biāo)軸投影點(diǎn)的數(shù)據(jù)處理，得出各個時刻的影長，由此得出影長隨時間的變化規(guī)律，為了較為直觀的體現(xiàn)出來，

21、我們采用圖標(biāo)形式進(jìn)行展現(xiàn)：（圖6）圖6：目標(biāo)地直桿影長隨時間變化圖標(biāo)示意圖利用最小二乘法的思想，通過曲線擬合，得出目標(biāo)直桿影子長度隨時間的變化關(guān)系滿足下式：（圖7） (17)圖7：直桿影長擬合曲線并因此式得出擬合曲線圖像，通過擬合函數(shù)關(guān)系找到使得目標(biāo)直桿影長達(dá)到最短的時間，即要求找到使得Y值達(dá)到最小的t值，將得到的t值代入式(3)，從而得到此地的經(jīng)度為東經(jīng)111度9秒。 在求解目標(biāo)地緯度的時候，同樣應(yīng)用最小二乘法思維，分別擬合得出直桿影子頂點(diǎn)的X值隨時間的函數(shù)關(guān)系以及直桿影長頂點(diǎn)的Y值隨時間的函數(shù)關(guān)系 ：（圖8，9） (18) (19)圖8：X坐標(biāo)擬合曲線圖9：Y坐標(biāo)擬合曲線使用同樣的方法，在

22、MATLAB軟件的運(yùn)行下得出兩函數(shù)的擬合曲線圖，分別找到使函數(shù)值達(dá)到最小的t值，將其代入(18)(19)兩式，得出此時的影子頂點(diǎn)坐標(biāo)，根據(jù)圖5算出此時的太陽方位角，并因此得出sin，聯(lián)立(14)式求出此地的地理緯度。為北緯22度11分24秒。綜上分析，通過該模型所最終確定的目標(biāo)地的經(jīng)緯度為東經(jīng)111度9秒，北緯22度11分24秒。經(jīng)查找，該地大概位于廣東省茂名市高州市G65包茂高速旁。圖10 分析結(jié)果所示地點(diǎn)5.3 問題三模型建立與求解 5.3.1 模型建立第三問將日期設(shè)為未知量，在前兩問模型的基礎(chǔ)上，要求求解桿子所在位置的經(jīng)緯度與日期。由地理知識可以得出，日期影響的唯一因素為赤緯角，即太陽直

23、射點(diǎn)的緯度，有以下兩個計(jì)算某日赤緯角的公式可供選擇：(1) =23.5°sin()式中，n為當(dāng)日日期序號，1月1日時n=1,以此類推，日期每增加一天，即n加1。(2) =23.5°-（日期-6月22日）×（23.5°×4/365）此公式可以大致計(jì)算一年中某年太陽直射點(diǎn)所在的緯度值。 計(jì)算結(jié)果若是正值，則為北緯；若為負(fù)值，則為南緯；R為某日日期，R-6月22日為該日與6月22日相差的天數(shù)，（23°26*4/365）為太陽直點(diǎn)一日內(nèi)移動的緯度距離（假設(shè)其移動是勻速的）。直射點(diǎn)經(jīng)度即太陽高度最大（太陽上中天）的經(jīng)線，地方時12：00的經(jīng)線；

24、建立目標(biāo)函數(shù)：(1)太陽方位角的差值平方和最小，要求求出的每個時刻方位角與附件2，3中對應(yīng)時刻方位角差值平方和最小。 (2)影長差值平方和最小,要求求出的每個時刻影長與附2，3中對應(yīng)時刻影長差值平方和最小。 (20)該目標(biāo)函數(shù)的約束條件：(1)影長要求符合問題一所提出的五個基本因素。(2)影長非負(fù)。(3)太陽直射點(diǎn)（赤緯角）在南北緯23.5°之間移動。(4)經(jīng)度的范圍在東西方向各180°范圍之內(nèi)。 (21)對于經(jīng)緯度進(jìn)行取點(diǎn)處理，進(jìn)行編程求解： 表1：對于桿子所處地點(diǎn)的日期與誤差分析數(shù)據(jù)表附件二附件三經(jīng)度79.34°E112.21°E緯度34.45

25、76;N23.14°N桿高2.01m3.11m太陽方向殘差比2.3%2.1%影長殘差比1.2%0.42%推算日期8/13或4/2612/26或12/145.3.2 粒子群算法目標(biāo)規(guī)劃模型在第三問中，加入的新變量日期，增加了模型的可行域，對于遍歷算法的實(shí)現(xiàn)增加了時間復(fù)雜度，因此，我們引入粒子群算法對求解過程進(jìn)行優(yōu)化，降低模型求解的時間復(fù)雜度。5.3.3 粒子群算法的簡介假設(shè)有一個D維的目標(biāo)搜索空間，有n個微粒組成了一個粒子群，其中每個微粒都用一個D維的向量描述，將它的空間位置表示為: , i=1,2,n;這可看做目標(biāo)優(yōu)化問題中的一個解，帶入適應(yīng)度函數(shù)計(jì)算出適應(yīng)度值可以衡量微粒的優(yōu)劣；第

26、i個微粒的飛行速度也是一個D維度的向量，記為 ；i=1,2,n;第i個微粒所經(jīng)歷過的具有最好適應(yīng)值的位置稱為個體的歷史最好位置，記為 ;i=1,2,n;整個微粒群所經(jīng)歷過的最好位置稱為全局歷史最好位置，記為 ；i=1,2,n;粒子群的進(jìn)化方程可描述為： (22) 其中，i,j 分別表示第i個微粒的第j個維度，t代表第t代，是兩個加速常量，通常的取值范圍是(0,2)，是兩個相互獨(dú)立的隨機(jī)函數(shù)。從上式可以看出，可以調(diào)節(jié)微粒去往自身周邊的最好位置，可以調(diào)節(jié)微粒去往整個粒子群所能找到的最好位置。將這種算法引用到本問中。算法中的“位置”就是本文中不同日期的所有離散的坐標(biāo)值點(diǎn)：，即“粒子”就是存放并尋找解

27、的動態(tài)變量。為了能夠?qū)ふ逸^為合適的優(yōu)化解，并且避免計(jì)算量的增大，本模型中設(shè)置了20個“粒子”來進(jìn)行優(yōu)化解搜索，即令i=20。由于本問中直桿坐標(biāo)和日期都是未知的，因此將粒子設(shè)為儲存三個參數(shù)： (23)5.3.4 粒子群算法模型的求解過程(1) 首先要對所有的粒子進(jìn)行初始化，即讓粒子中的不同參數(shù)可行范圍內(nèi)的隨機(jī)值，如令。這個過程相當(dāng)于在整個“位置”區(qū)域內(nèi)的某些點(diǎn)散布了尋找最優(yōu)解的“種子”。在何處散布了”種子“與其后來能夠?qū)ふ业降摹白顑?yōu)解”密切相關(guān)。當(dāng)散布位置不理想時，“最優(yōu)解”很可能會是局部最優(yōu)解，而非全局最優(yōu)解。(2) 其次計(jì)算每個粒子的適應(yīng)值。本問的適應(yīng)度值可以使用目標(biāo)函數(shù)來進(jìn)行計(jì)算，但由于最

28、優(yōu)解是令目標(biāo)函數(shù)最小，為最優(yōu)解的適應(yīng)度值應(yīng)該盡量大。所以對適應(yīng)度值Q做定義為： (24)(3) 然后記錄這些值，對每個粒子，將其適應(yīng)度與個體最好歷史記錄和全局最好記錄作比較。如果能夠比最好記錄更好就將此值設(shè)置為最好記錄。(4) 接下來，參照式（22），對每個粒子的去向進(jìn)行調(diào)整。每次調(diào)整時，到要對，進(jìn)行重新隨機(jī)取值，保證粒子行為的隨機(jī)性：既能有向最優(yōu)解靠近的趨勢，也能存在尋找其他更優(yōu)解的可能。(5) 最后，查看迭代次數(shù)。為了保證找到的最優(yōu)解，需要機(jī)遇粒子足夠多的行動時間尋找解。這里令迭代次數(shù)為1000。 對經(jīng)緯度進(jìn)行離散化處理后，選取步長為0.01，對其進(jìn)行MATLAB編程求解：表2：對于桿子所

29、處地點(diǎn)的日期與誤差分析數(shù)據(jù)表附件二附件三經(jīng)度80.34°E113.21°E緯度33.32°N23.54°N桿高2.24m3.01m太陽方向殘差比1.7%3.1%影長殘差比0.41%1.15%推算日期8/13或4/2712/26或12/125.3.5 模型的分析比較將兩個模型的結(jié)果對比，發(fā)現(xiàn)結(jié)果較為接近，附件1所得出的結(jié)果顯示桿子處在西藏山南地區(qū)，附件2所得出的結(jié)果顯示桿子所在地為廣西河池市內(nèi)。且得到的日期關(guān)于夏至日6/22與冬至日12/22對稱，這是由于太陽的運(yùn)動是關(guān)于夏至日以及冬至日對稱的。同時在太陽方向角殘差比以及影長殘差比這兩個參量上，粒子群算法的

30、誤差并沒有顯著提升，最大值為3.1%，在合理范圍之內(nèi)；在時間復(fù)雜度上，粒子群算法的求解時間大大縮短，所以當(dāng)未知參量增加時，使用粒子群算法是較合適的算法。5.4問題四模型建立與求解 5.4.1 概念引入透視投影：透視投影即使用中心投影法，可以使得延伸到無窮遠(yuǎn)的物體在照片平面中獲得其有效的延伸范圍。例如，在現(xiàn)實(shí)世界坐標(biāo)中存在兩平行鐵軌，則可以投影為照片平面中的相交直線，且其交點(diǎn)存在于兩平行直線的消隱點(diǎn)方向。在水平地面上，消隱點(diǎn)無限延伸后，即為地平線。并據(jù)此引出地平線性質(zhì)：在平面透視圖中，從地平線上引出的任意兩條直線都是平行的。5.4.2 模型建立在題目附件的視頻中，我們發(fā)現(xiàn)，在假設(shè)攝像機(jī)自身不隨周

31、圍氣流而發(fā)生抖動的基礎(chǔ)上，由于攝像機(jī)在拍攝時的角度問題，使得從視頻中所獲取的影長是真實(shí)影長向攝像機(jī)方向上的投影值，其值必定小于等于真實(shí)影長。由于從視頻中無法判斷影子隨時間的真實(shí)轉(zhuǎn)動角度，我們需要通過進(jìn)行一定的透視變換來還原出影子的真實(shí)長度及其變換角度。對視頻進(jìn)行分析，建立相應(yīng)的示意圖。根據(jù)前面的相應(yīng)分析，認(rèn)定攝像機(jī)水平拍攝，且圖中左下角的兩條綠色直線在現(xiàn)實(shí)中處于水平面中且相互平行，將其延伸在遠(yuǎn)處找到消隱點(diǎn)，且消隱點(diǎn)在水平方線無限延伸后就是地平線。（圖11） 圖11 視頻分析模型 確立直桿立點(diǎn)為O，在視角方向，令地平線與直桿的交點(diǎn)為P，P點(diǎn)即為直桿向正前方的消隱點(diǎn)。同時從底部引出一水平輔助線B，

32、則OP連線L與水平線相互垂直。設(shè)此時的直桿影子頂點(diǎn)為R，連接PR得輔助線A并延伸。由地平線的性質(zhì)可知，此時從P點(diǎn)引出的任意兩條直線都是平行的，即直線A與直線L相互平行，則A與B必有一交點(diǎn)，記為Q。則RQO為直角三角形。利用霍夫變換提取平行直線，當(dāng)圖像中存在很多條直線的時候，霍夫變換直接使用效果并不好，對算法進(jìn)行改進(jìn)，舍棄直線寬度在1以上的直線。同時，由于照片平面中物體因?yàn)槌上裢敢曌儞Q而產(chǎn)生形變，需要對其進(jìn)行透視逆變換使其恢復(fù)成原來樣子。已知圖片空間中二維圖像經(jīng)過透視變化變成另一個圖像，過程符合下列公式： m=(ax+by+c)/(ox+py+1) n=(dx+ey+f)/(qx+py+1) (

33、25)式中，斜視圖像點(diǎn)的像素坐標(biāo)為(x,y）,正視圖像的像素坐標(biāo)為(m,n），a,b,c,d,e,f,o,p,q均為透視變換參數(shù)。為了提高算法的精確度，對圖像進(jìn)行預(yù)處理，以達(dá)到減少計(jì)算難度的目的。利用MATLAB對其進(jìn)行圖像提取，并且對視頻中提取產(chǎn)生圖片進(jìn)行灰度處理和高斯濾波處理。可以發(fā)現(xiàn)，進(jìn)行高斯濾波處理后的圖像噪點(diǎn)明顯減少，減少了計(jì)算難度。繼續(xù)通過MATLAB進(jìn)行矯正處理，提取出實(shí)際情況下的影子長度變化。經(jīng)過一系列數(shù)據(jù)處理，我們得到了不同時刻下的直桿影子長度變換情況，代入到問題二的經(jīng)緯網(wǎng)模型中去，即可得到視頻拍攝地點(diǎn)。5.4.3 模型求解對經(jīng)緯度進(jìn)行離散化處理后，對其進(jìn)行編程求解：已知日期

34、未知日期經(jīng)度東經(jīng)110.70經(jīng)度東經(jīng)109.76緯度南緯42.31緯度南緯42.66桿高2m桿高2m影長殘差比1.32%影長殘差比2.45%日期日期6/11或7/135.4.4 結(jié)果分析通過表格對比發(fā)現(xiàn)當(dāng)把日期作為未知量對模型進(jìn)行求解時，可得地點(diǎn)坐標(biāo)與已知日期所求地點(diǎn)坐標(biāo)經(jīng)度上基本無變化，在內(nèi)蒙古包頭市內(nèi)，緯度相差3.4度，變化幅度為1.8%，由于數(shù)據(jù)的采集和預(yù)處理不可避免的造成了誤差，可以接受此誤差。所求日期為6/11，視頻顯示日期為7/13，剛好關(guān)于夏至日6/22對稱，說明模型可靠。6、 誤差分析與模型檢驗(yàn)6.1 誤差分析本文中涉及了大量的角度，由于對角度進(jìn)行計(jì)算時得到結(jié)果的小數(shù)點(diǎn)位數(shù)是較

35、多的，因此還需要對數(shù)據(jù)進(jìn)行預(yù)先處理，舍棄結(jié)果的小數(shù)點(diǎn)后四位，這樣一來，對于數(shù)據(jù)的準(zhǔn)確性還是有一點(diǎn)影響。不失一般性而且由于同一天中太陽赤緯角的變化很小，因此本文解答過程中進(jìn)行處理時，將太陽赤緯角當(dāng)作一天內(nèi)并沒有變化，這樣也會對所要求的結(jié)果產(chǎn)生一些影響。6.2 模型的檢驗(yàn) 本題中所提供的附件本身就是對所建立的模型的一個檢驗(yàn)，給模型提供一些檢驗(yàn)數(shù)據(jù)，從理論上來講，如果所建立的數(shù)學(xué)模型合理且正確，那么最終求出的地點(diǎn)和時間也必定符合題目的要求。七、模型評價(jià)本文所建立的數(shù)學(xué)模型和實(shí)際情況相符合，對于問題有合理的猜想、假設(shè)、計(jì)算以及檢驗(yàn)，具有一定的指導(dǎo)作用，其中在計(jì)算直桿太陽影子長度的變化曲線時，結(jié)合較為大

36、眾所接受的幾何原理以及現(xiàn)實(shí)中“立竿見影”的原理，建立了太陽位置和影子長度的變化的數(shù)學(xué)模型。并且在此基礎(chǔ)上，可以進(jìn)行進(jìn)一步研究，例如可以推廣到求解在不同的時間條件下或者不同位置時，太陽對直桿影子的長度變化的影響。 整體來說，本論文中所建立的幾個數(shù)學(xué)模型具有一般推廣性，不僅僅可以求解更為復(fù)雜的太陽影子為題，同時，利用太陽影子長度可以推算出當(dāng)?shù)氐慕?jīng)緯度、時間、日期和相對位置的思想，可以推廣到建筑行業(yè)日照設(shè)計(jì)，設(shè)計(jì)住宅布置、建筑間距、房屋朝向等現(xiàn)實(shí)應(yīng)用方面，甚至還可以應(yīng)用到農(nóng)業(yè)大棚養(yǎng)殖方面，利用不同高度植物采光度不同，合理安排蔬菜的種植位置。 在實(shí)際緯度求解過程中，對于數(shù)據(jù)的處理，我們采用了最小二乘法

37、的思想來進(jìn)行曲線擬合，得出直桿影長的變化與時間的圖像，從而找到目標(biāo)地點(diǎn)正午的時間，后來經(jīng)過簡潔分析，發(fā)現(xiàn)用多項(xiàng)式擬合的方法處理曲線會更加精確，但由于時間問題并未實(shí)施。 同時，本文在求解過程中，設(shè)置了很多假設(shè)，例如不考慮大氣對陽光的折射，實(shí)際上應(yīng)該重點(diǎn)考慮當(dāng)?shù)氐拇髿獾恼凵鋵Ξ?dāng)?shù)靥柕恼凵渎?，因?yàn)椴煌牡胤降奶柟饩€折射率是不同的，根據(jù)所考慮的問題，還需要建立氣壓與太陽折射率的數(shù)學(xué)模型，但是因?yàn)闆]有時間，沒有完成模型的修正?？傊瑢τ诓煌那闆r，需要要對模型做適當(dāng)?shù)男拚脱a(bǔ)充，針對實(shí)際情況選擇合適的參數(shù)，然后利用這些模型中的核心思維，可以比較靈活的去解決所要求解的問題。 八、參考文獻(xiàn)1于賀軍，氣象

38、用太陽赤緯和時差計(jì)算方法研究，氣象水文海洋儀器，2006年3期：P23-452張波，王宇欣，李天穎，楊樹國，基于遺傳算法的太陽影子定位，科技展望，2016年12期：P78-903徐威，曹小鴿，MATLAB在太陽影子定位技術(shù)中的應(yīng)用，通訊世界，2016年2期：P67-704張夢，李媛媛，張敏，何智恒，胡傳仟，基于太陽影子定位枚舉法模型的研究，科技創(chuàng)新與應(yīng)用，2016年6期：P30-355霍俊爽 楊柳 吳兆薇 張春轉(zhuǎn) 王柳行，由物體的影子變化定位其位置的數(shù)學(xué)模型研究，決策與信息，2015年12期：P25-266Reaganlee，太陽影子定位模型的分析， 附錄影長數(shù)據(jù)處理附件一：北京時間x坐標(biāo)(米

39、)y坐標(biāo)(米) 影長（米）14:421.03650.49731.149614:451.06990.50291.182214:481.10380.50851.215314:511.13830.51421.249114:541.17320.51981.283214:571.20870.52551.318015:001.24480.53111.353415:031.28150.53681.389415:061.31890.54261.426215:091.35680.54831.463415:121.39550.55411.501515:151.43490.55981.540215:181.4751

40、0.56571.579915:211.5160.57151.620115:241.55770.57741.661315:271.60030.58331.703315:301.64380.58921.746215:331.68820.59521.790115:361.73370.60131.835015:391.78010.60741.880915:421.82770.61351.9279附件二：北京時間x坐標(biāo)(米)y坐標(biāo)(米)影長（米）12:41-1.23520.1731.247312:44-1.20810.1891.222812:47-1.18130.20481.198912:50-1.15

41、460.22031.175412:53-1.12810.23561.152412:56-1.10180.25051.129912:59-1.07560.26531.107813:02-1.04960.27981.086313:05-1.02370.2941.065113:08-0.9980.3081.044413:11-0.97240.32181.024313:14-0.9470.33541.004613:17-0.92170.34880.985513:20-0.89650.36190.966813:23-0.87140.37480.948613:26-0.84640.38760.930913

42、:29-0.82150.40010.913813:32-0.79670.41240.897113:35-0.77190.42460.881013:38-0.74730.43660.865513:41-0.72270.44840.8505附件三：北京時間x坐標(biāo)(米)y坐標(biāo)(米)影長（米）13:091.16373.336 3.5331 13:121.22123.3299 3.5468 13:151.27913.3242 3.5618 13:181.33733.3188 3.5781 13:211.3963.3137 3.5958 13:241.45523.3091 3.6149 13:271.51

43、483.3048 3.6354 13:301.5753.3007 3.6572 13:331.63573.2971 3.6805 13:361.6973.2937 3.7052 13:391.75893.2907 3.7313 13:421.82153.2881 3.7589 13:451.88483.2859 3.7881 13:481.94883.284 3.8187 13:512.01363.2824 3.8508 13:542.07923.2813 3.8846 13:572.14573.2805 3.9199 14:002.21313.2801 3.9569 14:032.28153

44、.2801 3.9955 14:062.35083.2804 4.0358 14:092.42133.2812 4.0779 第一問影子長度變化的程序：clcclearclose allsyms xy_longI=imread('sunlight.png');rotI=imrotate(I,90,'crop');fig1=imshow(rotI);BW1=rgb2gray(rotI);BW2=im2bw(BW1,0.95);BW3=double(BW2);BW=edge(BW3,'canny');figure, im

45、show(BW);H,theta,rho=hough(BW);figure,imshow(imadjust(mat2gray(H),'XData',theta,'YData',rho,.    'InitialMagnification','fit');xlabel('theta(degrees)'),ylabel('rho');axis on, axis normal, hold on;colormap(ho

46、t)P=houghpeaks(H,50,'threshold',ceil(0.3*max(H(:);x=theta(P(:,2);y=rho(P(:,1);plot(x,y,'s','color','black');lines=houghlines(BW,theta,rho,P,'FillGap',18,'MinLength',180);figure, imshow(rotI), hold onmax_len=0;for k=1:length(lines)&#

47、160; xy=lines(k).point1;lines(k).point2;  plot(xy(:,1),xy(:,2),'LineWidth',2,'color','green');    plot(xy(1,1),xy(1,2),'x','LineWidth',2,'Color','yellow');  plot(xy(2,1),xy(2,2),'x','Line

48、Width',2,'Color','red');    len=norm(lines(k).point1-lines(k).point2);  if (len>max_len)      max_len=len;      xy_long=xy;  endendplot(xy_long(:,1),xy_long(:,2),'LineW

49、idth',2,'Color','red');第四問：clcclearlength=2.424024024 2.391081081 2.357837838 2.324834835 2.315825826  2.276546547 2.246546547 2.207417417 2.176066066  2.162462462 2.121861862 2.105255255 2.09039039 2.0541

50、44144 2.027087087 2.003063063 1.966996997 1.94 1.912942943 1.900900901 1.876876877 1.858858859;Length=length;sOLen=sum(Length);result=;t=8.90166666666667:0.0333333333333332:9.60166666666667;min=2;bc=0.1;lon=105:bc:150;beginl=-30*ones(1,size(lon,2);endl=70*ones

51、(1,size(lon,2);n=194;delta=23.45*sin(2*pi*(284+n)/365);for i=1:size(lon,2)    ts=t-(120-lon(1,i)/15;    omega=15*(ts-12);    beginl(1,i);    endl(1,i);    for lat=beginl(1,i):bc:endl(1,i)

52、60;       phi=lat;         h=asin(sin(phi*pi/180)*sin(delta*pi/180)+cos(phi*pi/180)*cos(delta*pi/180)*cos(omega*pi/180); l=2./tan(h);  if sum(abs(l-length)<min      mi

53、n=abs(l-length);      ganH=mean(Length.*tan(h);      result=lon(l,i) lat ganH abs(l-length) n;       end  clear h;   end iendplot(result(:,1),result(:,2),'r*&