太陽(yáng)影子定位

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

2、何原理計(jì)算太陽(yáng)高度角、赤緯角、地理緯度、時(shí)角和太陽(yáng)方位角之間的關(guān)系，計(jì)算出目標(biāo)地的地理緯度。同時(shí)選取北京正午十二時(shí)作為一個(gè)參考值，聯(lián)系擬合的曲線得出目標(biāo)地太陽(yáng)影子最短的時(shí)刻，與參考值進(jìn)行對(duì)比得出經(jīng)度。經(jīng)緯度的雙鎖定確定直桿所在地的具體位置大致是廣東省茂州市。針對(duì)問(wèn)題三，建立與第二問(wèn)類似的目標(biāo)規(guī)劃模型，但是加上日期因素，為了提高遍歷速度，使用粒子群算法減小時(shí)間復(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、廣西河池市。針對(duì)問(wèn)題四，對(duì)視頻進(jìn)行處理以得到必須的數(shù)據(jù)，由于拍攝時(shí)角度的存在，得到的影長(zhǎng)并非實(shí)際長(zhǎng)度，而是投影長(zhǎng)度，所以采用基于Hough變換和透視變換的圖像矯正法，對(duì)斜視圖像矯正，得出實(shí)際長(zhǎng)度，然后將得到的數(shù)據(jù)代入第二問(wèn)中的模型，得出的結(jié)果是視頻拍攝地點(diǎn)在內(nèi)蒙古包頭市境內(nèi)；在拍攝日期未知的情況下，將圖像變換得出的實(shí)際影長(zhǎng)代入第三問(wèn)中的基于粒子群算法的規(guī)劃模型里，可得出日期為6月11日和7月13日. 對(duì)于模型的推廣，根據(jù)某地的陽(yáng)光信息與地理信息，可以應(yīng)用到建筑行業(yè)日照設(shè)計(jì)，設(shè)計(jì)住宅布置、建筑間距、房屋朝向等方面，甚至還可以應(yīng)用到農(nóng)業(yè)大棚養(yǎng)殖安排蔬菜的種植位置等。 關(guān)鍵詞：二次函數(shù)擬合 最小二乘

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

21、我們采用圖標(biāo)形式進(jìn)行展現(xiàn)：（圖6）圖6：目標(biāo)地直桿影長(zhǎng)隨時(shí)間變化圖標(biāo)示意圖利用最小二乘法的思想，通過(guò)曲線擬合，得出目標(biāo)直桿影子長(zhǎng)度隨時(shí)間的變化關(guān)系滿足下式：（圖7） (17)圖7：直桿影長(zhǎng)擬合曲線并因此式得出擬合曲線圖像，通過(guò)擬合函數(shù)關(guān)系找到使得目標(biāo)直桿影長(zhǎng)達(dá)到最短的時(shí)間，即要求找到使得Y值達(dá)到最小的t值，將得到的t值代入式(3)，從而得到此地的經(jīng)度為東經(jīng)111度9秒。 在求解目標(biāo)地緯度的時(shí)候，同樣應(yīng)用最小二乘法思維，分別擬合得出直桿影子頂點(diǎn)的X值隨時(shí)間的函數(shù)關(guān)系以及直桿影長(zhǎng)頂點(diǎn)的Y值隨時(shí)間的函數(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)兩式，得出此時(shí)的影子頂點(diǎn)坐標(biāo)，根據(jù)圖5算出此時(shí)的太陽(yáng)方位角，并因此得出sin，聯(lián)立(14)式求出此地的地理緯度。為北緯22度11分24秒。綜上分析，通過(guò)該模型所最終確定的目標(biāo)地的經(jīng)緯度為東經(jīng)111度9秒，北緯22度11分24秒。經(jīng)查找，該地大概位于廣東省茂名市高州市G65包茂高速旁。圖10 分析結(jié)果所示地點(diǎn)5.3 問(wèn)題三模型建立與求解 5.3.1 模型建立第三問(wèn)將日期設(shè)為未知量，在前兩問(wèn)模型的基礎(chǔ)上，要求求解桿子所在位置的經(jīng)緯度與日期。由地理知識(shí)可以得出，日期影響的唯一因素為赤緯角，即太陽(yáng)直

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

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

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

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

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

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

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

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

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

32、則OP連線L與水平線相互垂直。設(shè)此時(shí)的直桿影子頂點(diǎn)為R，連接PR得輔助線A并延伸。由地平線的性質(zhì)可知，此時(shí)從P點(diǎn)引出的任意兩條直線都是平行的，即直線A與直線L相互平行，則A與B必有一交點(diǎn)，記為Q。則RQO為直角三角形。利用霍夫變換提取平行直線，當(dāng)圖像中存在很多條直線的時(shí)候，霍夫變換直接使用效果并不好，對(duì)算法進(jìn)行改進(jìn)，舍棄直線寬度在1以上的直線。同時(shí)，由于照片平面中物體因?yàn)槌上裢敢曌儞Q而產(chǎn)生形變，需要對(duì)其進(jìn)行透視逆變換使其恢復(fù)成原來(lái)樣子。已知圖片空間中二維圖像經(jīng)過(guò)透視變化變成另一個(gè)圖像，過(guò)程符合下列公式： 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ù)。為了提高算法的精確度，對(duì)圖像進(jìn)行預(yù)處理，以達(dá)到減少計(jì)算難度的目的。利用MATLAB對(duì)其進(jìn)行圖像提取，并且對(duì)視頻中提取產(chǎn)生圖片進(jìn)行灰度處理和高斯濾波處理。可以發(fā)現(xiàn)，進(jìn)行高斯濾波處理后的圖像噪點(diǎn)明顯減少，減少了計(jì)算難度。繼續(xù)通過(guò)MATLAB進(jìn)行矯正處理，提取出實(shí)際情況下的影子長(zhǎng)度變化。經(jīng)過(guò)一系列數(shù)據(jù)處理，我們得到了不同時(shí)刻下的直桿影子長(zhǎng)度變換情況，代入到問(wèn)題二的經(jīng)緯網(wǎng)模型中去，即可得到視頻拍攝地點(diǎn)。5.4.3 模型求解對(duì)經(jīng)緯度進(jìn)行離散化處理后，對(duì)其進(jìn)行編程求解：已知日期

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

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

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

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

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

39、)y坐標(biāo)(米) 影長(zhǎng)（米）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附件二：北京時(shí)間x坐標(biāo)(米)y坐標(biāo)(米)影長(zhǎng)（米）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附件三：北京時(shí)間x坐標(biāo)(米)y坐標(biāo)(米)影長(zhǎng)（米）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 第一問(wèn)影子長(zhǎng)度變化的程序：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');第四問(wèn)：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*&